aeons: approximating the end of nested sampling

Previous literature (Akrami et al., 2010; Feroz et al., 2011) has explored the potential for nested sampling to be used as a global maximiser, given its ability to handle multi-modalities. In particular, the latter authors emphasised that posterior samplers such as nested sampling find the bulk of the mass, not the maximum of the distribution, but that this can be remedied by tightening the termination criterion. We now use the machinery we have developed to put this statement on a more quantitative footing.

Let us take the current iteration to be the termination point with likelihood $\log\mathcal{L}_{\mathrm{f}}$ and prior volume $X_{\mathrm{f}}$, so that

Note that we have assumed that prior effects are negligible (so $1\approx\infty$), and that $\epsilon\ll 1$ so that the denominator is approximately the accumulated evidence. Computing this for (11), we find the answer in terms of lower incomplete gamma functions

Taking the $X_{\mathrm{f}}\ll(\sqrt{2}\sigma)^{d}$ limit (almost certainly valid at termination) we find

We thus have an expression relating $\mathcal{L}_{\mathrm{f}}$ at termination to the termination fraction $\epsilon$. This becomes yet more pleasing in the large $d$ limit, since $\epsilon^{2/d}\to 1$, we find via a Stirling approximation:

In the event that we retain $\epsilon$, we replace $\frac{d}{2e}\to\frac{d}{2e}\epsilon^{2/d}$, allowing one to battle the $\frac{d}{2e}$ term exponentially as dimensions increase.

Putting this together, taking $\mathcal{L}_{i}$ in (18) to be $\mathcal{L}_{\mathrm{f}}$ and combining this with (23), we find

showing that in general nested sampling will finish at a contour $d/2e$ away from the maximum log-likelihood. The final set of $n_{i}$ live points gets you $\log n_{i}/e$ closer, with a chance of getting $\sim\pi/\sqrt{6}e=0.472$ closer still by statistical fluctuation.

Making the traditional termination criterion stricter therefore has limited returns in high-dimensions, if it is ultimately still based on the remaining evidence. However, nested sampling still shrinks around the maximum exponentially, so provided a good alternative termination criterion is chosen, it will get there in reasonable time.

To quantify this statement, consider the number of iterations $\Delta i$ required to get from the posterior bulk to the true maximum, which we will now calculate. At the posterior, we have roughly speaking

The prior volume $\log X_{\delta}$ at which the likelihood is within $\delta$ of the maximum can then be found by inverting Eq. 11:

which corresponds to an additional

iterations after the bulk is reached, where $n$ is the harmonic mean of the number of live points. One would therefore expect a general distribution to take $\mathcal{O}(\tfrac{nd}{2}\log\tfrac{d}{2})$ iterations to get from the usual nested sampling stopping point to within an $e$-fold of the maximum. A rule-of-thumb termination criterion could therefore be to run for at least $\tfrac{nd}{2}\log\tfrac{d}{2}$ iterations after the posterior is reached.

A summary of the distances between the notable points at the end of a run is shown in Fig. 4.

As shown in the previous section, midway through the run nearly all of the posterior mass is concentrated at a single point. However, this does not capture the structure of the posterior that has been explored and all of the information it provides.

We have the potential to fix this because nested sampling is invariant to monotonic transformations, so we can transform the likelihood as $\mathcal{L}\to\mathcal{L}^{\beta}$ without loss of information by trivially re-weighting the samples. Increasing $\beta$ worsens the situation, while $\beta\to 0$ simply gives back the prior. There is, on the other hand, a significant intermediate range which makes the samples look like a posterior centred at the present contour, which will allow us to recover the structure of the samples. A schematic of the procedure is shown in Fig. 5.

At this point it is relevant to note the correspondence between Bayesian inference and statistical mechanics, from which the above transform is derived. If one equates the parameters to microstates $i$, the negative log-likelihood to the microstate energy $E_{i}$, and the prior to the density of states $g_{i}$, then the posterior as given by the generalised Bayes’ rule is the canonical ensemble

at the inverse temperature $\beta=1/T$. As noted by Habeck (2015), thermal algorithms such as thermodynamic integration (Gelman & Meng, 1998) get the evidence by evolving through a series of canonical ensembles via some temperature schedule, but nested sampling instead maintains microcanonical ensembles, which are the iso-likelihood contours. Instead of using temperature (Kirkpatrick et al., 1983; Swendsen & Wang, 1986) or energy (Wang & Landau, 2001) as a control parameter, nested sampling chooses a series of ensembles with constant relative volume entropy $\Delta\log X$, which allows the algorithm to handle phase transitions (Baldock et al., 2016).

Because the temperature of a microcanonical ensemble is a derived property rather than a parameter, there is some freedom in its definition. Returning to our original motivation, we make the connection that the temperature is the re-weighting $\mathcal{L}\to\mathcal{L}^{\beta}$ which centres the ensemble around the current energy. We now present several temperatures that achieve this aim, each of which one can plausibly consider to be the current temperature of a nested sampling run.

The obvious candidate is the microcanonical temperature $\partial S/\partial E$, where the volume entropy is $\log X$ and the energy is as usual $-\log\mathcal{L}$. This gives the density of states; as discussed in Skilling’s original paper,

is the $\beta$ at which $\mathcal{L}^{\beta}X$ peaks at $\log X_{i}$, if we assume differentiability, which is exactly the intuition we were aiming for to put the ensemble bulk at the current contour.

Its value can be easily obtained via finite difference of the $\log\mathcal{L}$ and $\log X$ intervals, albeit subject to an arbitrary window size for the differencing. Indeed, material science applications (Baldock et al., 2017) use this estimator to monitor the ‘cooling’ progress of nested sampling, with a window size of 1000 iterations.

Another temperature considered by Habeck (2015) is that at which the current energy (i.e. $-\log\mathcal{L}_{i}$) is the average energy of the entire ensemble. One can obtain it by inverting

to get the ‘canonical’ temperature $\beta_{\mathrm{C}}$. While $\beta_{\mathrm{M}}$ is derived from (the gradient of) a single contour, this temperature uses the entire ensemble. It has the desirable property that it rises monotonically with compression, in analogy to a monotonic annealing schedule.

We furthermore propose a temperature $\beta_{\mathrm{B}}$ that is obtained via Bayesian inference, which returns a distribution rather than a point estimate. Since each value of $\beta$ leads to a different likelihood $\mathcal{L}^{\beta}$, one can consider the posterior distribution as a function of $\log X$ to be conditioned on $\beta$. We can therefore write

What we would really like is the distribution of $\beta$ at the present iteration, so the natural step is to invert this via Bayes’ rule;

As with all Bayesian analyses, the distribution of $\beta$ is fixed up to a prior, which we choose to be uniform in $\beta$. The obtained temperatures are consistent with the previous two choices, which may seem oddly coincidental. In fact, closer inspection reveals that large values of $P(\beta\mid\log X_{i})$ are the temperatures with a large value of the posterior at the present contour, normalised by the corresponding evidence. Thus the Bayesian temperature uses the same idea as the microcanonical one, except it accounts for the spread in the result.

Fig. 6 shows the three temperatures as a function of compression for two cases, one containing a phase transition and one without. They are consistent in both cases when there is a single dominant phase, but differ during a phase transition. The canonical temperature is the only one that rises monotonically with compression.

One should keep in mind that despite the above theoretical reasoning, our introduction of the likelihood transformation was ultimately motivated by our wish to utilise the extra degree of freedom it provides. As we will see below, we recommend choosing the exact definition depending on what is useful for the problem at hand.

Plots of samples dimensionality against compression also draw attention to the directions in which the samples are constrained throughout the run. As a concrete example, consider an elongated Gaussian in a unit hypercube prior with $\mu=\bm{0}$ and $\Sigma=\mathrm{diag}\left(10^{-3},10^{-3},10^{-3},10^{-6},10^{-6},10^{-6}\right)$, for which the dimensionality estimates are plotted in Fig. 8. Alongside is a view of the distribution of live points across the prior for two directions with different scales, which shows the level to which those parameters have been constrained at different times.

A feature of nested sampling made apparent here is that parameters with high variance are initially ‘hidden’. Compression occurs in the direction which is most likely to have a sample of higher likelihood, and initially it is much easier to find a better point along the direction of a parameter that is poorly constrained. Lower variance parameters are constrained much later, and before that happens it appears as though those parameters have a uniform distribution.

It is important to appreciate that at lower compression the samples truly lie in a lower-dimensional space, rather than some artefact of the way we view them. Anticipating the full dimensionality of the space is therefore just as impossible as that associated with a slab-spike geometry, so in this sense such geometries contain a ‘compressive phase transition’.

An alternative approach used in Ultranest (Buchner, 2021) derives a progress bar based on the fraction of the accumulated integral compared to the remaining integral, approximated as

This has a several shortcomings. First, run-time is proportional to compression rather than accumulation of the integral, since it takes just as long to traverse the width of the bulk as it does any other width. Second, because of the point-like nature of the posterior mid-run, the remaining integral approximated as such holds nearly all of the evidence, so the relative fraction of the accumulated and remaining evidence is almost zero for most of the run. Finally, approximation (53) is always an underestimate, because as previously found the maximum live point is generally nowhere near the true maximum. This approach can however be useful in the low dimensions appropriate for Ultranest when the live points are always near the maximum, but in general is less reliable.

Seasoned watchers of nested sampling runs might be curious how the method compares to simply extrapolating the increments of evidence to roughly estimate when the evidence converges. We do this for a spherical Gaussian and compare it to our method. At an intermediate stage of the run, the most recent outputs might look something like that shown in the first two columns of the table in Fig. 9(a). Extrapolating those data to a linear and exponential profile yields endpoint estimates plotted in the graph to the right.

The linear extrapolation is clearly an underestimate, since it fails to account for the long tail of the nonlinear profile. The increments are also not exactly exponential, since the exponential fit leads to a large over-prediction. The predicted endpoint over the course of a run for $d=16$, $\sigma=0.01$, as shown in Fig. 10, shows the same result. One might expect an average to be more accurate, but this tends to be biased towards the exponential prediction, and there is no obvious choice of weighting that would fix this.

More importantly, we find that for real likelihoods which have an element of noise the extrapolation often diverges, for instance when the increments do not monotonically decrease. Directly extrapolating the evidence increments is therefore far less stable than the previous method, and generally not a reliable method for prediction.

Predictions for spherical Gaussians of various dimensions are shown in Fig. 11 as a benchmark for when fitting a Gaussian distribution is exact. In all endpoint prediction plots, the shading indicates the $1-2\sigma$ uncertainties. All were run with $n_{i}=500$ except for one $n_{i}=2000$ for comparison, with each Gaussian having a width of $\sigma=0.01$. The correct endpoint is recovered to within standard error at all points except the very beginning, when the parameters have hardly been constrained.

We note that as with nested sampling in general, increasing $n_{i}$ improves the resolution and reliability of the inferences, which can be seen from the middle two plots.

We also observe the effect of elongating the Gaussian, using the same example as Section 3.4. Fig. 12 shows a step-like trend similar to the inferred dimensionalities, reflecting the fact that the full dimensionality is undetectable at lower compression factors. The endpoint for a likelihood whose remaining three directions are completely unconstrained coincides with our predictions at early iterations, showing that the two cases are indistinguishable.

One case that might be expected to cause problems is the pathological Cauchy distribution, which is far from a Gaussian. Fig. 13 shows the predictions for a likelihood of the form

choosing $d=10$ and allowing $\gamma$ to vary. The correct estimate is obtained to within standard error by about halfway, but before that is inaccurate. The key limitation is that the estimate is wrong early on, not because the compression is anisotropic, or because there is a phase transition; but rather as a limitation of the reducing the likelihood to a Gaussian via the BMD, which is itself less stable for a Cauchy.

Nevertheless, the right order of magnitude is obtained at all times, so this remains sufficient for most use-cases. The Cauchy is also a pathological case, and the same problem does not in practice appear for more realistic cases, as we shall see next.