PolyChord: nested sampling for cosmology

TL;DR

PolyChord is a new nested sampling algorithm designed for high-dimensional cosmological parameter estimation, capable of handling multiple modes and exploiting parameter hierarchies with efficient parallelization.

Contribution

It introduces a novel nested sampling method that uses slice sampling, handles multi-modal posteriors, and exploits parameter hierarchies, improving efficiency in high-dimensional cosmological analyses.

Findings

Efficiently samples high-dimensional, multi-modal posteriors.

Exploits parameter hierarchies for faster convergence.

Parallelized implementation using openMPI.

Abstract

PolyChord is a novel nested sampling algorithm tailored for high dimensional parameter spaces. In addition, it can fully exploit a hierarchy of parameter speeds such as is found in CosmoMC and CAMB. It utilises slice sampling at each iteration to sample within the hard likelihood constraint of nested sampling. It can identify and evolve separate modes of a posterior semi-independently and is parallelised using openMPI. PolyChord is available for download at: http://ccpforge.cse.rl.ac.uk/gf/project/polychord/

Figures

Equations

$$\text{P}(\theta |\mathcal{D},\mathcal{M})\equiv {\mathcal{P}}_{\mathcal{M}}(\theta )=\frac{{\mathcal{L}}_{\mathcal{M}}(\theta ){\pi }_{\mathcal{M}}(\theta )}{{\mathcal{Z}}_{\mathcal{M}}},$$

$${\mathcal{Z}}_{\mathcal{M}}=\int {\mathcal{L}}_{\mathcal{M}}(\theta ){\pi }_{\mathcal{M}}(\theta )d\theta .$$

$$\text{P}({\mathcal{M}}_i|\mathcal{D})=\frac{\text{P}(\mathcal{D}|{\mathcal{M}}_i)\text{P}({\mathcal{M}}_i)}{\text{P}(\mathcal{D})}=\frac{{\mathcal{Z}}_i{\pi }_i}{{\sum }_j{\mathcal{Z}}_j{\pi }_j},$$

$$X_i=t_i{X}_{i-1},P(t_i)=n{t}_i^{n-1},X_0=1$$

$$\mathcal{Z}\approx \sum _{i=1}{\mathcal{L}}_i(X_{i-1}-X_i).$$

$$T^{-1}\mathbf{x}=\mathbf{y}$$

$$\text{Speedup}(n_{\mathrm{procs}})=n_{\mathrm{live}}\log \left[1+\frac{n_{\mathrm{procs}}}{n_{\mathrm{live}}}\right],$$

Full text

††pubyear: 2015††pagerange: PolyChord: nested sampling for cosmology–References††pubyear: 2014

1 Introduction

Over the past two decades, the quantity and quality of astrophysical and cosmological data has increased substantially. In response to this, Bayesian methods have been increasingly adopted as the standard inference procedure.

Bayesian inference consists of parameter estimation and model comparison. Parameter estimation is generally performed using Markov-Chain Monte-Carlo (MCMC) methods, such as the Metropolis Hastings (MH) algorithm and its variants (Mackay, 2003). In order to perform model comparison, one must calculate the evidence; a high-dimensional integration of the likelihood over the prior density (Sivia & Skilling, 2006). MH methods cannot compute this on a usable timescale, hindering the use of Bayesian model comparison in cosmology and astroparticle physics.

A contemporary methodology for computing evidences and posteriors simultaneously is provided by Nested Sampling (Skilling, 2006). This has been successfully implemented in the now widely adopted algorithm MultiNest  (Feroz & Hobson, 2008; Feroz et al., 2009, 2013). Modern cosmological likelihoods now involve a large number of parameters, with a hierarchy of speeds. MultiNest struggles with high-dimensional parameter spaces, and is unable to take advantage of this separation of speeds. PolyChord aim to address these issues, providing a means to sample high-dimensional spaces across a hierarchy of parameter speeds.

This letter is a brief overview of our algorithm, which is now in use in several cosmological applications (Planck Collaboration XX, 2015). It will be followed in the near future by a longer and more extensive paper.

PolyChord is available for download from the link at the end of the paper.

2 Bayesian Inference

Given some dataset $\mathcal{D}$, one may use a model $\mathcal{M}$ with parameters $\theta$ to calculate a likelihood ${\text{P}(\theta|\mathcal{D},\mathcal{M})\equiv\mathcal{L}_{\mathcal{M}}(\theta)}$. This may be inverted using Bayes’ theorem to find the posterior distribution on the parameters:

[TABLE]

where ${\text{P}(\theta|\mathcal{M})\equiv\pi_{\mathcal{M}}(\theta)}$ is our prior degree of belief on the values of the parameters and ${\text{P}(\mathcal{D}|\mathcal{M})\equiv\mathcal{Z}_{\mathcal{M}}}$ is the evidence or marginal likelihood, calculated as:

[TABLE]

One may use the evidence $\mathcal{Z}_{\mathcal{M}}$ to perform model comparison, comparing $\mathcal{Z}_{\mathcal{M}}$ to the evidences of other competing models $\{\mathcal{M}_{1},\mathcal{M}_{2},\ldots\}$. Applying Bayes theorem to model comparison yields:

[TABLE]

where $\text{P}(M_{i})\equiv\pi_{i}$ is the prior probability of a model $\mathcal{M}_{i}$ (usually taken as uniform). Since the evidences $\mathcal{Z}_{i}$ play a key role in Bayesian model comparison, so their calculation is vital for this aspect of inference. For a full discussion of Bayesian inference, we recommend Sivia & Skilling (2006) and part IV of Mackay (2003).

3 Nested Sampling

Skilling (2006) constructed an ingenious method of sampling which enables one to estimate the evidence integral (2) efficiently.

One begins by drawing $n$ “live points” uniformly from the prior. At iteration $i$, the point with the lowest likelihood $\mathcal{L}_{i}$ is replaced by a new live point drawn uniformly from the prior with the constraint that its likelihood $\mathcal{L}>\mathcal{L}_{i}$. One may estimate the fraction of the prior $X_{i}$ enclosed by the iso-likelihood contour $\mathcal{L}_{i}$, since it will be the lowest volume of $n$ volumes drawn uniformly from $[0,X_{i-1}]$. $X_{i}$ thus satisfies:

[TABLE]

We may then calculate the evidence by quadrature:

[TABLE]

and the mean and variance of the random variable $\mathcal{Z}$ may be inferred using equation (4). The posterior mass contained within the live points may be estimated as ${\mathcal{Z}_{\mathrm{live}}\approx\langle\mathcal{L}\rangle_{\mathrm{live}}X_{\mathrm{live}}}$, and convergence is reached when $\mathcal{Z}_{\mathrm{live}}$ is some small fraction of $\mathcal{Z}$ (Keeton, 2011).

In addition, the pool of dead points may be used to construct a set of posterior samples for use in parameter estimation. For further details on the nested sampling procedure, see Sivia & Skilling (2006).

4 Sampling within an iso-likelihood contour

The difficult step in nested sampling is sampling from the prior subject to the hard likelihood constraint $\mathcal{L}>\mathcal{L}_{i}$. The prior may be sampled using inverse transform sampling, so that sampling effectively occurs in the unit hypercube (Feroz et al., 2009). How one samples in this space under the hard likelihood constraint is where variations of the algorithm differ.

The MultiNest algorithm samples by using the live points to construct a set of intersecting ellipsoids which together aim to enclose the likelihood contour, and then samples by rejection sampling within it. Whilst being an excellent algorithm for small number of parameters, any rejection sampling algorithm has an exponential scaling with dimensionality.

Galilean Sampling (Feroz & Skilling, 2013; Betancourt, 2011) samples by an MCMC-chain based approach using reflection off iso-likelihood contours. This however suffers from the need to tune parameters, and the requirement of likelihood gradients.

Diffusive nested sampling (Brewer et al., 2009) is an alternative and promising variation on \NAT@swafalse\NAT@partrue\NAT@fullfalse\NAT@citetpskilling2006 algorithm, which utilises MCMC to explore a mixture of nested probability distributions. Since it is MCMC based, it scales well with dimensionality. In addition, it can deal with multimodal and degenerate posteriors, unlike traditional MCMC. It does however have multiple tuning parameters.

We find that \NAT@swafalse\NAT@partrue\NAT@fullfalse\NAT@citetpNealSlice slice sampling is well suited to the problem of drawing points from within an iso-likelihood contour. In the one-dimensional case, he suggests the sampling procedure detailed in Figure 1. In the $N$-dimensional case, he suggests multiple methods, the most simple being to perform one dimensional slice sampling of each parameter in turn (in a random order). This method has been implemented in nested sampling by Aitken & Akman (2013) to great effect on some simple examples. However, like Gibbs and MH sampling, this methodology struggles with degenerate distributions. PolyChord aims to address this issue.

5 The PolyChord Algorithm

PolyChord implements several novel features compared to \NAT@swafalse\NAT@partrue\NAT@fullfalse\NAT@citetpSystemsBio slice-based nested sampling. We utilise slice sampling in a manner that uses the information in the live points to deal with degeneracies in the contour. We also implement a general clustering algorithm which identifies and evolves separate modes in the posterior semi-independently. The algorithm is written in FORTRAN95 and parallelised using openMPI. It is optimised for the case where the dominant cost is the generation of a new live point. This is frequently the case in astrophysical applications, either due to high dimensionality, or to costly likelihood evaluation. In addition, it has the option of implementing fast-slow parameters, which is extremely effective in its combination with CosmoMC. This is termed CosmoChord, which may be downloaded from the link at the end of the paper.

Our generic $N$-dimensional slice-sampling procedure is detailed in Figure 2.

5.1 Dealing with degenerate contours

To solve degenerate cases, we make a linear transformation to “whiten” the space, so that in general the contour has size $\sim\mathcal{O}(1)$ in every direction. Uniform sampling is preserved under affine transformations, so this strategy is valid.

In order to find such a linear transformation, note that at every step of the algorithm one has some knowledge of the dimensions of the contour, provided by the set of live points and all of the inter-chain points. We take the covariance matrix of these points as a reasonable descriptor of the correlation. The Cholesky decomposition $T$ of the covariance matrix $\Sigma=TT^{T}$ is a good choice of transformation so that:

[TABLE]

transforms from $\bf x$ in the unit hypercube to $\bf y$ in the sampling space. Working in the sampling space, a reasonable approach that mixes well is to start from a random live point as the seed, choose a randomly oriented orthonormal basis, and sample along these directions in a random order. For difficult distributions, one may require multiple iterations of this procedure to ensure that the final point is fully decorrelated from the start point. This is illustrated further in Figure 2. Within the sampling space, the initial step size may be chosen as $w=1$.

This procedure has the advantage of being dynamically adaptive, and requires no tuning parameters. However, this “whitening” process is ineffective for pronounced curving degeneracies.

5.2 Clustering

Multi-modal posteriors are a challenging problem for any sampling algorithm. “Perfect” nested sampling (i.e. the entire prior volume enclosed by the iso-likelihood contour is sampled uniformly) in theory solves multi-modal problems as easily as uni-modal ones. In practice however, there are two issues.

First, one is limited by the resolution of the live points. If a given mode is not populated by enough live points, it runs the risk of “dying out”. Indeed, a mode may be entirely missed if the algorithm is not sufficiently resolved. In many cases, this problem can be alleviated by increasing the number of live points.

Second, and more importantly for PolyChord, the sampling procedure may not be appropriate for multi-modal problems. We “whiten” the unit hypercube using the covariance matrix of live points. For far-separated modes, the covariance matrix will not approximate the dimensions of the contours, but instead falsely indicate a high degree of correlation. It is therefore essential for our purposes to have PolyChord recognise and treat modes appropriately.

This methodology splits into two distinct parts; (i) recognising that clusters are there, and (ii) evolving the clusters semi-independently.

5.2.1 Recognition of clusters

Any cluster recognition algorithm can be substituted at this point. One must take care that this is not run too often, or one runs the risk of adding a large overhead to the calculation. In practice, checking for clustering every $\sim\mathcal{O}(n_{\mathrm{live}})$ iterations is sufficient. We encourage users of PolyChord to experiment with their own preferred cluster recognition, in addition to that provided and described below.

It should be noted that the live points of nested sampling are amenable to most cluster recognition algorithms for two reasons. First, all clusters should have the same density of live points in the unit hypercube, since they are uniformly sampled. Second, there is no noise (i.e. outside of the likelihood contour there will be no live points). Many clustering algorithms struggle when either of these two points are not satisfied.

We therefore choose a relatively simple variant of the $k$-nearest neighbours algorithm to perform cluster recognition. If two points are within one another’s $k$-nearest neighbours, then these two points belong to the same cluster. We iterate $k$ from $2$ upwards until the clustering becomes stable (the cluster decomposition does not change from one $k$ to the next). If sub-clusters are identified, then this process is repeated on the new sub-clusters.

5.2.2 Evolving the clusters semi-independently

An important novel feature comes from what one does once clusters are identified.

First, when spawning from an existing live point, the whitening procedure is now defined by the covariance matrix of the live points within that cluster. This solves the issue detailed above.

Second, PolyChord would naively spawn live points into a mode with a probability proportional to the number of live points in that mode. In fact, what it should be doing is to spawn in proportion to the volume fraction of that mode. These should be the same, but the difference between these two ratios will exhibit random-walk like behaviour, and can lead to biases in evidence calculations, or worse, cluster death. Instead, one can keep track of an estimate of the volume in each cluster, using an approached based on (4), and choose the mode to spawn into in proportion to that estimate. This methodology will be fully documented in the more extensive paper.

Thus, the point to be killed off is still the global lowest-likelihood point, but we control the spawning of the new live point into clusters by using our estimates of the volumes of each cluster. We call this ‘semi-independent’, because it retains global information, whilst still treating the clusters as separate entities.

When spawning within a cluster, we determine the cluster assignment of the new point by which cluster it is nearest to. It does not matter if clusters are identified too soon; the evidence calculation will remain consistent.

5.3 Parallelisation

Currently, PolyChord is parallelised by openMPI using a master-slave structure. One master process takes the job of organising all of the live points, whilst the remaining ${n_{\mathrm{procs}}-1}$ “slave” processes take the job of finding new live points.

When a new live point is required, the master process sends a random live point and the Cholesky decomposition to a waiting slave. The slave then, after some work, signals to the master that it is ready and returns a new live point and the inter-chain points to the master.

A point generated from an iso-likelihood contour $\mathcal{L}_{i}$ is usable as a new live point for an iso-likelihood contour $\mathcal{L}_{j}>\mathcal{L}_{i}$, providing it is within both contours. One may keep slaves continuously active, and discard any points returned which are not usable. The probability of discarding a point is proportional to the volume ratio of the two contours, so if too many slaves are used, then most will be discarded. The parallelisation goes as:

[TABLE]

and is illustrated in Figure 3.

5.4 Fast-slow parameters and CosmoChord

In cosmological applications, likelihoods exhibit a hierarchy of parameters in terms of calculation speed (Lewis, 2013). The effect of this is that the likelihoods may be quickly recalculated if one changes certain subsets of the parameters. In PolyChord it is very easy to exploit such a hierarchy. Our transformation to the sampling space is laid out so that if parameters are ordered from slow to fast, then this hierarchy is automatically exploited: A Cholesky decomposition, being a upper-triangular skew transformation, mixes each parameter only with faster parameters. Further to this, one may use the fast directions to extend the chain length by many orders of magnitude. This helps to ensure an even mixing of live points.

6 PolyChord in action

6.1 Gaussian Likelihood

As an example of the efficacy of PolyChord, we compare it to MultiNest on a Gaussian likelihood in $D$ dimensions. In both cases, convergence is defined as when the posterior mass contained in the live points is $10^{-2}$ of the total calculated evidence. We set $n_{\mathrm{live}}=25D$, so that the evidence error remains constant with $D$. MultiNest was run in its default mode with importance nested sampling and expansion factor $e=0.1$. Whilst constant efficiency mode has the potential to reduce the number of MultiNest evaluations, the low efficiencies required in order to generate accurate evidences negate this effect.

With these settings, PolyChord produces consistent evidence and error estimates with an error $\sim 0.4$ log units (Figure 5). Using importance nested sampling, MultiNest produces estimates that are within this accuracy.

Figure 4 shows the number of likelihood evaluations $N_{\mathcal{L}}$ required to achieve convergence as a function of dimensionality $D$.

Even on a simple likelihood such as this, PolyChord shows a significant improvement over MultiNest in scaling with dimensionality. PolyChord at worst scales as ${N_{\mathcal{L}}\sim\mathcal{O}(D^{3})}$, whereas MultiNest has an exponential scaling which emerges in higher dimensions.

However, we must point out that a good rejection algorithm like MultiNest will always win in low dimensions. MultiNest is also extremely effective at navigating pronounced curving degeneracies in low dimensions, whereas PolyChord must use very long chains in order to navigate such contours.

6.2 CosmoChord

PolyChord’s real strength lies in its ability to exploit a fast-slow hierarchy common in many cosmological applications. We have successfully implemented PolyChord within CosmoMC, termed CosmoChord. The traditional Metropolis-Hastings algorithm is replaced with nested sampling. This implementation is available to download from the link at the end of the paper.

This combination has been effectively implemented in multiple cosmological applications in the latest Planck paper describing constraints on inflation (Planck Collaboration XX, 2015), including application to a $37$-parameter reconstruction problem ($4$ slow, $19$ semi-slow, $14$ fast).

7 Conclusions

We have introduced PolyChord, a novel nested sampling algorithm tailored for high dimensional parameter spaces. It is able to fully exploit a hierarchy of parameter speeds such as is found in CosmoMC and CAMB (Lewis & Bridle, 2002; Lewis et al., 2000). It utilises slice sampling at each iteration to sample within the hard likelihood constraint of nested sampling. It can identify and evolve separate modes of a posterior semi-independently and is parallelised using openMPI. We demonstrate its efficacy on a toy problem. A further more detailed paper will follow imminently.

Acknowledgements

We would like to thank Farhan Feroz for numerous helpful discussions during the inception of the PolyChord algorithm.

Download Link

PolyChord is available for download from: http://ccpforge.cse.rl.ac.uk/gf/project/polychord/