Occasionally the finalization step produces a worse ELBO than the previous steps

[pipme]

In finilization step the elbo mean drops from 0.04 to -0.01.

To reproduce, add np.random.seed(43) in the beginning of tests/vbmc/test_vbmc_optimize.pyand run test_vbmc_multivariate_normal.

[lacerbi]

Thoughts:

• The finalization step does "variational boosting" which means it first decides what the "best" iteration is from a bracket of past iterations (after warmup), then takes the vp and gp from that chosen "best" iteration and runs variational inference on them with typically 50 components.
• In this specific case, the difference in Mean[ELBO] is small but still, I agree that in theory it shouldn't happen that the mean ELBO worsens after variational boosting - adding more components should always improve the mean ELBO. However, since we are optimizing the ELBO using stochastic optimization, it could happen in practice.
• The issue here might be in how the "best" iteration is chosen for starting the variational boosting. The algorithm for selecting the best iteration follows some heuristics, based on the Mean and Std of the ELBO, the iteration number (later iterations are favored with respect to earlier iterations), the stability index (here reported as "Convergence"), etc.
• I don't remember seeing similar behavior in MATLAB, but I wouldn't put my money on that (it's been shown multiple times that the MATLAB version has some minor issues that were found while porting pyvbmc).
• In short, while I think it's somewhat likely that the issue could be at this stage (variational boosting and/or selection of the best iteration), it could either be (1) an implementation bug in pyvbmc or (2) an intrinsic problem with the heuristic choice (which would also be in the MATLAB version). I have a slightly higher posterior over (1), but also the second option is not so unlikely.
• Alternative option (also somewhat likely, in fact more likely now that I think about it): it's not an issue of the choice of iteration in variational boosting, but a problem with stochastic optimization. Anecdotally, our homemade (and tweaked) implementation of Adam in Python seems to behave slightly differently from our MATLAB version, for reasons unknown (it's a problem mentioned by Mikko and never resolved, see here under "Things to investigate"). It's possible that this difference becomes important when optimizing a large number of variables (as when doing variational boosting with 50 components).

Having said this, checking whether the issue is with the choice of iteration (as opposed to somewhere else, like during optimization) should be quite easy.