zeus: A Python implementation of Ensemble Slice Sampling for efficient Bayesian parameter inference

We introduce zeus, a well-tested Python implementation of the Ensemble Slice Sampling (ESS) method for Bayesian parameter inference. ESS is a novel Markov chain Monte Carlo (MCMC) algorithm specifically designed to tackle the computational challenges posed by modern astronomical and cosmological analyses. In particular, the method requires only minimal hand--tuning of 1-2 hyper-parameters that are often trivial to set; its performance is insensitive to linear correlations and it can scale up to 1000s of CPUs without any extra effort. Furthermore, its locally adaptive nature allows to sample efficiently even when strong non-linear correlations are present. Lastly, the method achieves a high performance even in strongly multimodal distributions in high dimensions. Compared to emcee, a popular MCMC sampler, zeus performs 9 and 29 times better in a cosmological and an exoplanet application respectively.


INTRODUCTION
Over the past few decades the volume of astronomical and cosmological data has increased substantially. In response to that, a variety of astrophysical models have been developed to explain the plethora of observations. Markov chain Monte Carlo (MCMC) has been established as the standard procedure of inferring the model parameters subject to the available data in a Bayesian framework. Within the Bayesian context, the object that quantifies the probability distribution of the model parameters given the data and model M is the posterior distribution P ( ) ≡ ( | , M) which is defined using Bayes's theorem: where L ( ) ≡ ( | , M) is the likelihood function, ( ) ≡ ( |M) is the prior distribution of the model parameters , and Z ≡ ( |M) is the, so called, Bayesian model evidence or marginal likelihood and in this context can be treated as a simple normalisation constant. MCMC does not in general require knowing the value of the model evidence and it only depends on the ability to evaluate the unnormalised posterior distribution for arbitrary values of . MCMC methods can then be used to generate (Markov) chains of samples from the posterior distribution. Those samples can be used to calculate integrals (e.g. parameter uncertainties, marginal distributions etc.) that are paramount for modern astronomical and cosmological analyses.
The most commonly used MCMC methods are variants of the ★ E-mail: minas.karamanis@ed.ac.uk Metropolis-Hastings (MH) algorithm (Metropolis et al. 1953;Hastings 1970). MH consists of two steps. First, given the last sample in the chain, a new sample is proposed and then the Metropolis criterion determines whether or not that new sample should be accepted and thus added to the chain. The resulting chain is Markovian in the sense that each sample is proposed based only on the previous sample. The purpose of the Metropolis acceptance criterion is to bias the chain so that the time spent in a region of the parameter space would be proportional to the posterior probability in that region. In other words, the stationary distribution of the Markov chain is the target distribution i.e. the posterior distribution. For a detailed introduction to MCMC methods we direct the reader to MacKay (2003) and for an intuitive introduction to Bayesian inference to Jaynes (2003).
Arguably, the most difficult part of the MH algorithm is the proposal step. There are many ways of choosing a new sample and the efficiency of the method depends on this choice. By far the simplest one is the use of a normal (Gaussian) distribution, centred around the previous sample to generate the new proposed sample. The resulting method is often called Random Walk Metropolis algorithm and its performance is highly sensitive to the ( + 1)/2 elements that form its covariance matrix. Those elements generally need to be chosen a priori or be adaptively tuned. More efficient methods utilise the gradient of the target distribution (Betancourt 2017) or an ensemble of parallel and communicating chains (Gilks et al. 1994;Ter Braak 2006;ter Braak & Vrugt 2008;Goodman & Weare 2010).
Out of the methods mentioned in the previous paragraph we will focus our attention on the last one, the ensemble or population MCMC variety. The reason is simple: the Random Walk Metropolis algorithm requires a great amount of tuning (or a priori knowledge) for it to perform efficiently and even then there is no guarantee that the proposal covariance matrix is optimal for the whole parameter space. On the other hand, gradient based methods, although very powerful, are in general unsuitable for astronomical applications in which the models that are used are almost always not differentiable.
One benefit of ensemble MCMC over its alternatives is that the ensemble of parallel chains (also known as walkers) collectively sample the posterior, thus information about their distribution can be shared and used to make better educated proposals. Other advantages include the lack of hand-tuning of hyper-parameters and their capacity for parallel implementation. For the aforementioned reasons, ensemble MCMC methods have dominated astronomical analyses. The most common ones are affine-invariant ensemble sampling (AIES) (Goodman & Weare 2010) and differential evolution MCMC (DEMC) (Ter Braak 2006;ter Braak & Vrugt 2008), both implemented in the popular Python package emcee (Foreman-Mackey et al. 2013. In this paper we introduce zeus, a stable and well-tested Python implementation of Ensemble Slice Sampling (ESS) (Karamanis & Beutler 2020). ESS is a method based on the ensemble MCMC paradigm, with the crucial difference being that its proposals are performed via Slice Sampling updates (Neal 2003) instead of Metropolis-Hastings ones. As we will thoroughly demonstrate in Section 3, this subtle difference leads to substantial improvements in terms of sampling efficiency and robustness. zeus is a user-friendly tool that does not require any hand-tuning or preliminary runs and can scale up to 1000s of CPUs without any extra effort from the user.
zeus has been used in various astronomical and cosmological analyses, including cosmological tests of gravity (Tamosiunas 2020), relativistic effects and primordial non-Gaussianity (Wang et al. 2020), 21cm intensity mapping (Umeh et al. 2021), and has been implemented as part of the CosmoSIS package (Zuntz et al. 2015).
zeus is open source software that is publicly available at https: //github.com/minaskar/zeus under the GPL-3 Licence. Detailed documentation and examples on how to get started are available at https://zeus-mcmc.readthedocs.io.

ENSEMBLE SLICE SAMPLING
zeus is a Python implementation of the Ensemble Slice Sampling (ESS) method presented in Karamanis & Beutler (2020). Here we will provide a high-level description of the method and will refer to the accompanying paper for more details about the underlying algorithmic structure and mathematics.
ESS combines the ensemble MCMC paradigm with slice sampling. Since the use of slice sampling in astronomical parameter inference is rare we will start by explaining its function and how it differs from Metropolis updates. Then we will move on to discuss how it can be efficiently combined with ensemble MCMC.

Slice sampling
Slice sampling is based on the idea that sampling from a distribution with density ( ) is equivalent to uniform sampling from the area under the plot of ( ) ∝ ( ). To this end, we introduce an auxiliary variable , called height, such that the joint distribution ( , ) is uniform over the region = {( , ) : 0 < < ( )}. To sample from the marginal distribution ( ), we first sample from ( , ) and then we marginalise by dropping the value of each sample.
x y 0 0 Figure 1. Illustration of the univariate slice sampling update. Given the current sample 0 , a value 0 is uniformly sampled along the vertical slice (0, ( 0 )) (dashed line) thus defining the initial point (blue). An interval ( , ) is uniformly positioned horizontally around ( 0 , 0 ) and it is expanded in steps of size − until both its ends are outside the slice. The new sample is generated by repeatedly sampling (uniformly) from the interval ( , ) until a sample (green star) is found inside the slice. Samples outside of the slice (red star) are rejected and they are instead used to shrink ( , ).
(iii) Draw the new sample 1 uniformly from ∩ .
To construct the interval (step ii), Neal (2003) introduced the stepping-out procedure that works by randomly positioning an interval of length around the sample 0 (i.e. blue dot in Figure 1) and then expanding it in steps of size until both its ends (i.e. and ) are outside the slice. To obtain 1 (i.e. green star in Figure 1) we then use the shrinking procedure in which candidates are sampled uniformly from until a point inside the slice is found. Samples outside of the slice are used to shrink the interval . The two procedures are shown in Figure 1. The length scale is the only free hyperparameter of slice sampling and although its choice can reduce or increase the computational cost of the method it generally does not affect its mixing properties (e.g. convergence rate, autocorrelation time, etc.). zeus utilises a stochastic optimization algorithm similar to Tibbits et al. (2014) and based on the Robbins & Monro (1951) optimisation scheme in order to tune to its optimal value (see Section 3.1 of Karamanis & Beutler 2020 for more details).
It is important to note here that for multimodal target distributions there is no guarantee that the approximate slice would cross any of the other modes. In particular, if the initial estimate of the length scale is low then the probability of missing the other peaks, assuming that they are located far away, is also low. As we will show in Section 3, unlike simple slice sampling, ESS and thus zeus does not suffer from this effect.

Walkers, moves and parallelism
The slice sampling update described in the previous paragraphs is a univariate update scheme. For it to be used to sample from mul-X X Y k X l m Figure 2. The figure illustrates the differential move in the context of Ensemble Slice Sampling. The walker to be updated is shown in red. Two walkers, and , (blue) are uniformly selected from the complementary ensemble (grey). The approximate slice (dotted line) is constructed parallel to the two walkers and using the stepping-out procedure. The new position (green) of is sampled using the shrinking procedure along the approximate slice. tivariate target distributions it needs to be generalised accordingly. Perhaps the simplest such generalisation in a multivariate setting is the use of slice sampling to sample along each coordinate axis in turn (i.e. component-wise slice sampling) or to sample along randomly selected directions in parameter space (MacKay 2003). Although valid, both of these approaches are unsuitable in cases of correlated parameters in which the proper choice of direction can substantially accelerate mixing.
To address this issue, Tibbits et al. (2014) proposed to orthogonalise the parameter space using the sample covariance, thus getting rid of linear correlations between parameters. We will instead follow a different, perhaps more flexible, approach to construct an efficient slice sampler. Our aim is to utilise an ensemble of parallel chains/walkers that can exchange information about the covariance structure of the target distribution and thus by-pass the difficulties posed by correlations.
As hinted in the introduction, the ensemble of walkers collectively sample the target distribution and thus their positions encode information about the correlations between the parameters. One way to take advantage of this information is to use it to construct direction vectors along which slice sampling can take place. Many moves that generate direction vectors from the complementary ensemble are possible. zeus offers a collection of them, including some that utilise clustering algorithms and density estimation methods. As we will show in Section 3, such moves can help accelerate sampling in difficult cases such as strongly multimodal distributions. Any distribution of the complementary ensemble can be used as a valid proposal to generate such direction vectors and zeus offers a highly flexible interface for the user to define such a move or choose one (or a mixture) from the ones that are already implemented and tested.
Here is a list of the currently implemented moves in zeus: • Differential move: This is the default move used by zeus and shown in Figure 2. Using the differential move, Ensemble Slice Sampling updates the position of each walker in the ensemble by slice sampling along a direction defined by the difference between two uniformly selected walkers from the rest of the ensemble (i.e. the complementary ensemble).
• Gaussian move: The Gaussian move samples the direction vectors along which slice sampling is performed from a normal distribution that shares the same covariance structure as the complementary ensemble. This approach is very efficient in cases in which the target distribution is close to normal.
• Global move: The Global move utilises a Dirichlet Process Gaussian Mixture to fit the complementary ensemble and proposes directions along different peaks of the target distribution in cases of strong multi-modality.
• KDE move: The KDE move samples the direction vectors from a Gaussian Kernel Density Estimate of the complementary ensemble. This can be useful in cases of highly non-Gaussian target distributions.
• Random move: The Random move performs slice sampling along isotropic directions. This is equivalent of standard multivariate slice sampling and it is mostly offered for testing purposes as it cannot handle correlations efficiently.
For more information on how those moves work as well as a comparison of the Differential, Gaussian and Global moves we direct the interested reader to Karamanis & Beutler (2020). Unless stated otherwise the Differential move will be used for the following examples.
To parallelise this process and capitalise on the availability of multiple CPUs we randomly split the ensemble into two sets of walkers (i.e. active and passive sets) (Foreman-Mackey et al. 2013) and choose to update the positions of the active walkers along direction vectors defined by passive walkers. Then the passive become active and vice versa and the process is repeated. The ensemble splitting technique is required in order to parallelise the algorithm without violating detailed balance. Parallelisation is achieved in practise using either multiprocessing or MPI using the implemented ChainManager utility that can distribute both multiple ensembles and multiple chains in parallel computing environments at the same time. Heuristics to determine the number of required walkers per application are discussed in Section 4.

EMPIRICAL EVALUATION
For the empirical evaluation of zeus we use five toy examples that manifest significant aspects of real astronomical applications 1 (i.e. linear and non-linear correlations, multimodality, heavy tails, hard boundaries) and two real-world astronomical examples characteristic of modern astronomical analyses.

Toy examples
In order to understand the behaviour of zeus in various sampling scenarios it is important to study its performance in different toy examples that demonstrate different characteristics of common target distributions that arise in astronomical applications. For that reason we chose five such toy examples. The first one is a normal (Gaussian) It is important to note here that both emcee algorithms exhibit a peak at zero separation; zeus on the other hand does not due to its nonrejection nature. distribution which by definition is characterised only by linear correlation between its parameters. The second toy problem is the ring distribution, a characteristic example of strong non-linear correlations. The third example is a Gaussian mixture with two components. While the purpose of the first two examples is to study the behaviour of the algorithm in the presence of linear and non-linear correlations respectively, the goal of the third example is to demonstrate the ability of zeus to sample efficiently from multimodal target distributions. The fourth toy example investigates the effect that heavy tails have on the sampling efficiency and the fifth the shows the effects that hard boundaries have on sampling.
We compare zeus with two popular alternatives offered by emcee, namely affine-invariant ensemble sampling with the stretch move (emcee/AIES) and the differential evolution move (emcee/DEMC). The main goal of this analysis is to justify our choice of slice sampling as the basis of zeus instead of Metropolis updates through the use of simple yet instructive toy examples.
For all three toy examples discussed below we adopt the same analysis procedure, where we initialise the walkers by sampling from a normal distribution N (0, I) where I is the identity covariance matrix and we discarded 10 4 iterations as burn-in.
The main metric that we use to investigate the behaviour of the samplers in those toy examples and to compare their performance is the distribution of steps performed by the walkers. As a step we define the distance spanned in parameter space by a single walker in a single iteration. This is a fundamental measure of the efficiency of an MCMC method and it is directly related to the expected squared jump distance (ESJD; Pasarica & Gelman 2010) given by: where are the chain samples, 1 is the first-order autocorrelation, and Var ( ) ( ) is a function of the stationary distribution only. Assuming that the higher-order autocorrelations 2 , 3 , . . . are monotonically decreasing with respect to 1 , then maximising the ESJD leads to minimisation of the autocorrelation between chain elements and thus maximisation of the sampling efficiency. In other words, the further away (i.e. the greater the ESJD) the walkers jump per iteration, the higher the sampling efficiency of the method. A benefit of using ESJD instead of the autocorrelation time as a metric is that the former, as an expectation value, is more accurate when computed using short chains. In order to account for the different computational cost (i.e. different number of model evaluations per iteration) between zeus and emcee we thinned the chains of the latter method according to the average number of model evaluations of zeus. This allowed us to compare the distribution of steps of the three samplers as shown in

The correlated normal distribution
Starting with the normal target distribution it is important to note here that all three of the methods used in the comparison are affineinvariant 2 , meaning that their performance is immune to any linear correlations between the parameters. Since the normal distribution incorporates, by construction, only linear correlations (i.e. the 2D marginal distribution contours look like ellipses), it is the perfect testing ground to assess the effect that high dimensionality has on the three methods independently of other complications. For our example we used a zero-mean normal distribution with a covariance matrix in which the diagonal elements are set to 1 and the offdiagonal ones are equal to 0.95. We then proceed by sampling the aforementioned distribution in 10, 25 and 50 dimensions. Based on Figure 3 one can see that the walkers of emcee/AIES dissolve into an inefficient random walk characterised by low step size and high autocorrelation time as the number of parameters increases. zeus and emcee/DEMC are not so severely affected by the high number of parameters exhibiting a substantially lower autocorrelation. Let us now try to explain this difference in behaviour by looking into the distribution of the steps of the walkers in Figure 4. One thing to notice here is that the distribution of the steps of zeus's walkers extends significantly further away than those of emcee/AIES and emcee/DEMC. This should come as no surprise since the construction of the approximate slice allows for larger steps than Metropolis updates as shown in Table 1. This is because when a proposal is rejected in slice sampling the approximate slice shrinks and another sample is proposed instead. In this way zeus's walkers always move and the chance of staying fixed is zero -unlike MH-based updates in which frequent rejection of samples is a necessity. This aforementioned procedure leads to greater steps in parameter space. The difference between emcee/AIES and emcee/DEMC is attributed to the fact that DEMC uses a proposal scale 3 = 2.38/ √ that guar-2 Differential evolution Metropolis is only approximately affine-invariant due to the jitter that it is often added to its proposal. This however has a negligible effect. 3 The proposal scale is similar to used in ESS in the sense that its value determines the length scale of the proposed jumps in parameter space. A high value would lead to large steps that are often rejected and a low value would lead to small steps that are often accepted but do not carry the walkers far. For such methods, a balance must me found. antees a constant acceptance rate accounting for the number of dimensions . This proposal scale is however optimal only in the case of a normal target distribution such as the one that we are studying here and there is no guarantee that it would return acceptable results in non-Gaussian distributions. For the case of emcee/AIES, the relevant proposal scale is allowed to vary in the range between 1/ and where = 2 is often taken as the typical value. It is clear that in the latter case does not possess the desired scaling ∝ 1/ √ and thus, although the method generates proposals in the right overall direction, most of the samples do not reside in the typical set (Speagle 2019). In other words, the lack of proper scaling of the proposal scale with the number of dimensions leads to emcee/AIES "overshooting" the typical set where most of the posterior mass is located.
We can also draw some useful insights about the sampling efficiency of those samplers and its scaling with the number of dimensions by estimating the integrated autocorrelation time of the chains. Given the autocorrelation time we can also estimate the effective sample size as the percentage of effectively independent samples in a chain. Dividing the effective sample size with the computational cost of each method we can then estimate the sampling efficiency. The results of such a comparison are shown in Figure 5. We immediately notice here that the autocorrelation times of zeus and emcee/DEMC scale linearly with the number of dimensions, whereas the autocorrelation time of emcee/AIES scales exponentially. The computational cost of zeus per iteration per walker, although somewhat higher than that of emcee, does not vary with the number of dimensions. This means that in high dimensions, zeus dominates over emcee/AIES in terms of sampling efficiency. The above discussion allows us to clearly state a crucial distinction between the three methods, which is their response to the curse of dimensionality. As the number of dimensions increases, the probability mass of a distribution is concentrated into a thin shell within the tails of the distribution (i.e. the typical set). To account for this and maintain its efficiency, a sampling method has to adjust its proposal scale -otherwise the proposals will not be located in the typical set and thus they will not be accepted. The three methods that we mentioned so far deal with this in different ways. emcee/AIES's proposal scale is not adjusted and thus its proposals become increasingly inefficient in high dimensions. emcee/DEMC's proposal scale is adjusted based on the theoretical expectation for the case of the normal target distribution. Although both emcee methods perform well in this example, their sub-optimal scaling will degrade their performance in non-Gaussian target distributions as we will demonstrate in the next toy example. Finally, zeus's proposal scale is continuously adapted, as the slice expands and contracts in every iteration, thus guaranteeing optimal scaling. Huĳser et al. (2017) found that the suboptimal scaling of emcee/AIES with the number of dimensions can introduce biases into the expectation values derived from the chains in high dimensions that are hard to diagnose. The locally adaptive nature of zeus allows it to avoid this problem by adjusting its proposals accordingly.
Another kind of analysis we can perform is to use the highly correlated 25-dimensional normal distribution as the target distribution and estimate the convergence rate of the three samplers. Although simple, the normal distribution is a valid approximation of many realistic astronomical posterior distributions and as such we expect the results presented in this paragraph to be applicable to a wide range of other distributions that resemble the normal distribution to some extent. We acknowledge however that the no free lunch theorem also applies to this case, and there are bound to be cases in which the results would be qualitatively different. That being said, we initialised the walkers from a compact normal distribution (i.e. standard deviation equal to 10 −4 times that of the target distribution) centred around a point along the first axis of the parameter space at a distance of 100 standard deviations from the mode. We then measured the number of model evaluations required until the samplers have converged to the target distribution. The results for varying number of walkers are presented in Figure 6. It is important to note here that both emcee algorithms exhibit a peak at zero separation; zeus on the other hand does not. The existence of the zero-peak in emcee is due to the high number of rejected proposals (i.e. low acceptance rate). 6 in which the computational cost until convergence (i.e. number of model evaluations) for a single walker diminishes and then reaches a plateau as the number of walkers is increased. We notice however that, at the level of a single walker, the computational cost of emcee/AIES is significantly higher compared to that of either zeus or emcee/DEMC. This is due to the way that different samplers choose the directions along which walkers move. In particular, both zeus and emcee/DEMC define a direction vector as the difference between two walkers from the complementary ensemble, thus two walkers are required to define a direction. On the other hand, emcee/AIES requires only a single walker from the complementary ensemble as the direction is defined by the difference between the updated walker and the complementary one. This stark contrast between the way those samplers choose their direction vectors lies at the heart of the difference in the computational cost of emcee/AIES as compared to zeus and emcee/DEMC in the limit of low number of walkers. In order to dive a little deeper into this, we can compute the exact number of possible directions for all three methods. Since emcee/AIES requires only a single walker from the complementary ensemble the number of available directions is equal to the size of the complementary ensemble. On the other hand, zeus's and emcee/DEMC's requirement for a pair of walkers means that the number of available directions is equal to 2 , meaning the 2-combination from a set of walkers that comprise the complementary ensemble. Clearly, as shown in Figure 7, the latter increases faster with the size of the complementary ensemble, thus explaining the larger variety of possible directions available in the case of zeus and emcee/DEMC compared to emcee/AIES.
The discussion so far was about the computational cost of convergence in terms of the number of model evaluations for a single walker. Of course, the ensemble of walkers consists by definition of more than a single walker. Therefore, in order to compute the total number of model evaluations required until the ensemble converges we need to multiply the results of the single walker with the total number of walkers. Those results are presented in the right panel of Figure 6. From this plot we can see that both zeus and emcee/DEMC converge faster when the number of walkers is close to its minimum value i.e. 2 × . emcee/AIES on the other hand prefers a higher number of walkers (i.e. 32 × ) in order to overcome the sparsity of available directions in the limit of low number of walkers. This, however, means that even if we choose the optimal number of walkers for emcee/AIES it would still converge slower than either zeus or emcee/DEMC. Furthermore, we cannot know a priori the optimal number of walkers for emcee/AIES unlike for zeus and emcee/DEMC in which the optimal size of the ensemble is close to 2 × . Finally, the faster convergence of zeus compared to It is important to note here that both emcee algorithms exhibit a peak at zero separation; zeus on the other hand does not due to its non-rejection basis. emcee/DEMC can be attributed to the local adaptation that the former performs by extending the length of the slice and thus allowing larger steps in parameter space.

The ring distribution
The ring distribution defined as where = 2, = 1 and is the total number of parameters; this is an artificial target distribution that exhibits strong non-linear correlations between its parameters. This aspect of the ring distribution allows us to demonstrate the locally adaptive nature of zeus. Whereas emcee/AIES and emcee/DEMC use a single global proposal scale for all regions of the parameter space, zeus has the ability to adjust its proposal scale locally by expanding the slice appropriately. As expected, this will allow zeus to sample efficiently even in cases in which strong non-linear correlations are present. Looking at Figure  8 one can see that zeus manages to generate multiple samples efficiently even in high dimensions. On the other hand, emcee/AIES and emcee/DEMC do not efficiently produce valid proposals: for emcee/AIES this leads to an inefficient random walk, characterised by small steps; for emcee/DEMC the acceptance rate almost vanishes beyond = 2. The expected squared jump distance of each method for the case of = 25 is shown in Table 1. It is important to note here that out of the three samplers only zeus manages to converge in all three cases (i.e. in 2, 10 and 25 dimensions). emcee/AIES and emcee/DEMC on the other hand converge successfully only in 2 dimensions.
To explain this result one only has to look at the distribution of walker steps of the different methods at Figure 9. zeus's steps extend to large distances in parameter space whereas most of emcee/AIES's and emcee/DEMC's steps are rejected (i.e. shown as zero in the histogram). We can see that emcee/DEMC manages to perform some long distance steps but those are few and there is almost nothing in between. It is clear from this and the previous toy examples that the = 2.38/ √ scaling of emcee/DEMC's scale factor does not generalise well beyond the Gaussian case.

The two-component Gaussian mixture distribution
One other important aspect of astronomical posterior distributions is the fact that many of them exhibit multiple peaks. Multimodality can arise either from non-linear models or sparse and uninformative data. In either case, multimodal target distributions present a formidable challenge for most MCMC methods. Perhaps the simplest example of such a distribution is the two-component Gaussian mixture. In this example we will position the two, equal-mass, components at −0.5 and +0.5 respectively with standard deviation of 0.1. Sampling from multimodal distributions requires two types of proposals, local proposals that sample different modes individually and global proposals that transfer walkers from one mode to the other. For this reason we will make use of zeus's GlobalMove that uses a Dirichlet Process Gaussian Mixture model of the ensemble to efficiently propose between-mode and within-mode steps.
As seen in Figure 10, zeus's walkers manage to move from one mode to the other frequently enough for mixing to be efficient even in the = 25 case. Out of emcee/AIES and emcee/DEMC, only the latter proposes valid steps from one mode to the other in the = 2 case. As for the = 25 case, one can see in Figure 11 that zeus's walkers perform numerous jumps whereas emcee's walkers are unable to do so. The ability of the walkers to jump from mode to mode is of paramount importance if we want to sample correctly from the target distribution. Lack of such proposals will lead to an improper probability mass ratio between the two modes and thus biased inference. The expected squared jump distance of each method for the case of = 25 is shown in Table 1.
Clustering-based proposals have also been applied to MH-type ensemble MCMC methods but as shown in Karamanis & Beutler (2020), they fail to generate valid proposals in problems with moderate number of dimensions. The reason is, as discussed in Section 3, that MH has to propose a valid point in the other mode. In other words, whereas Ensemble Slice Sampling only needs to determine the direction of the other mode relative to the chosen walker correctly, MH needs to guess both the direction and the distance, a task that rapidly becomes very hard as the number of dimensions rises.

The Student's -distribution
The fourth toy example tests the case in which the target distribution is characterised by heavy-tails. In order to demonstrate zeus's ability to sample efficiency is such cases we chose to use the multivariate Student's -distribution with 2 degrees of freedom. The aforementioned density exhibits heavier tails than a normal distribution which means that it is more likely to produce samples that are far away from the mean. The -distribution arises when estimating the mean of a normally distributed sample with unknown standard deviation and small size. The probability density function of a -dimensional Student's -distribution with degrees of freedom is given by: where is the × positive semi-definite shape matrix and is the mean vector. We sampled the above distribution using the three samplers in 2, 10 and 25 dimensions respectively as shown in Figure 12. The diagonal elements of shape matrix were set to 1 and the offdiagonal elements to 0.95. The mean vector was set to 0. All three samplers managed to sample efficiently in 2, 10 and 25 dimensions as shown in Figure 12 and Table 1. Overall, zeus was the most efficient method with emcee/AIES being second and emcee/DEMC last. One can see from Figure 13 that the distributions of steps of zeus and emcee/AIES are very similar whereas that of emcee/DEMC is substantially shorter. Unlike the previous toy examples in which the proposal strategy of emcee/AIES was causing it to overshoot the bulk of posterior mass, in the case of the heavy-tailed -distribution more proposals are accepted. On the other hand, emcee/DEMC's proposals which are optimised for Gaussian targets are more conservative in the case of the -distribution and they do not extend far away. As also demonstrated in the previous toy examples, the locally adaptive nature of zeus allows it to perform efficient proposals that span large distances in parameter space.

The truncated normal distribution
The fifth and final toy example tests the case in which the target distribution is bounded from below or above. We chose to employ a truncated normal distribution similar to the one used in the first toy example, with the additional constraint being that > 0. This effectively introduces a hard boundary along all dimensions. One of the reasons that we study this distribution is to assess the bias introduced by the presence of the hard boundary.
We sampled the above distribution using the three samplers in 2, 10 and 25 dimensions respectively as shown in Figure 14. The diagonal elements of the covariance matrix were set to 1 and the off-diagonal to 0.95. The mean vector was set to 0. All three samplers managed to sample efficiently in 2, 10 and 25 dimensions as shown in Figure 14 and Table 1. Overall, zeus was the most efficient method with emcee/DEMC being second and emcee/AIES last. One can see from Figure 15 that the distributions of steps of zeus and emcee/AIES are very similar whereas that of emcee/AIES is slightly shorter. As shown in the right panels of Figure 14 zeus exhibits the least amount of bias compared to emcee/AIES and emcee/DEMC. In practical astronomical examples however, only one or two parameters would usually be bounded (e.g. the neutrino mass in galaxy clustering analyses) and thus unbiased sampling would be easier to perform by either of the three samplers.

Real astronomical analyses
The previous section employs toy examples in order to exhibit various scenarios that might emerge during sampling, and shows how zeus is better equipped to handle them. To demonstrate the efficiency of zeus compared to other samplers in realistic target distributions, we chose two common astronomical inference problems as the testing ground. Those are the cases of baryon acoustic oscillation (BAO) parameter inference and exoplanet parameter estimation. We used the same three samplers in our comparison, namely emcee with AIES and DEMC, and of course zeus. We performed three distinct tests: • The first test was to estimate the efficiency for each sampler, defined as the number of independent samples produced per loglikelihood evaluation. To this end, we ran the MCMC procedure 5 times for each sampler and computed the mean efficiency using the estimated autocorrelation time of the chains. The autocorrelation time was estimated using the method presented in Karamanis & Beutler (2020).
• The second test relates to the convergence rate of the three algorithms. As a measure of convergence rate, we adopt the inverse of the number of iterations required until all the convergence criteria specified below are met. In order to estimate the mean convergence rate we ran the sampling procedure 40 times for each sampler initialising the walkers close to the Maximum a Posteriori (MAP) estimate.
• Finally, we tested the sensitivity of the samplers to the initial conditions by running 40 realisations with the walkers initialised from a small sphere (of radius 10 −4 ) around a randomly chosen point in the prior volume, counting how many of those attempts led to converged chains before a predetermined number of likelihood evaluations.
To determine whether a chain has converged we used four different metrics: the Gelman-Rubin split-statistic (Gelman et al. , 2013 using four independent ensembles of walkers; the Geweke test (Geweke 1992); a minimum length of the chain as a multiple of the integrated autocorrelation time (IAT); as well as an upper bound on the rate of change of the IAT. Only the second half of the chains was used to evaluate the aforementioned criteria. The number of walkers used in both examples was close to the minimum value of 2 × as specified below. As we will discuss in Section 4 this often leads to faster convergence.

Cosmological inference
The particular inference problem that we consider here is that of the anisotropic BAO parameter inference using estimates of the galaxy power spectrum. The data we used come from the 12th data release (DR12) of the high-redshift North Galactic Cap (NGC) sample as observed by the Sloan Digital Sky Survey (SDSS) (Eisenstein et al. 2011) Baryon Oscillation Spectroscopic Survey (BOSS) (Dawson et al. 2013). Our analysis follows closely that of Beutler et al. (2017) with the difference that we chose not to fix any parameters and fit the hexadecapole multipole of the power spectrum as well as the monopole and quadrupole. Those choices were made solely to render the problem more challenging. Indeed the inclusion of the hexadecapole does not contribute any additional constraining power for the data that we used. However, such extended models will prove useful when analysing data from larger galaxy surveys such as DESI (DESI Collaboration et al. 2016). In terms of Bayesian inference, the problem has 22 free parameters. The results of our analysis are consistent with those of Beutler et al. (2017). We used weakly informative flat (uniform) priors for all parameters except for the two scaling parameters, and ⊥ for which we used normal (Gaussian) priors. We used 50 walkers in total.
In terms of efficiency, zeus generates at least 5 independent samples for each one generated by emcee/DEMC and at least 9 for each one generated by emcee/AIES factoring in the different computational cost of the methods. As for the convergence rate, zeus converges more than 3 times faster than either emcee variant. Finally, we found that zeus is less sensitive to the initialisation than either of the other two methods. In particular, out of the 40 tests conducted with different initialisation, zeus converged 36 times, emcee/DEMC 14 times and emcee/AIES 7 times prior to the predetermined maximum number of likelihood evaluations (i.e. 5 × 10 6 in this case). The aforementioned results are presented in detail in Table 2. The 1-D and 2-D marginal posterior distributions are shown in Figure 16 demonstrating the agreement between the three methods 4 . 4 No upper limit on the number of likelihood evaluations or iterations was used for this run and convergence was diagnosed using all the metrics that we introduced.

Exoplanet inference
Another common application of MCMC methods in astronomy is the problem of exoplanet parameter inference through modelling of Keplerian orbits and radial velocity time series data. In this section we demonstrate the performance of zeus using a two-planet model with 14 free parameters and real data from the K2-24 (EPIC-203771098) extrasolar system (Petigura et al. 2016) that is known to host two exoplanets. We used the popular Python package RadVel (Fulton et al. 2018) for the Keplerian modelling of the planetary orbits. The results of our analysis are consistent with published constraints for the aforementioned extrasolar system (Petigura et al. 2016). We used 30 walkers in total for sampling.
We performed the same suite of tests as in the cosmological inference case. In terms of efficiency, zeus generates more than 7 independent samples per each one generated by emcee/DEMC and more than 29 independent samples per each one generated by emcee/AIES. As for the convergence rate, zeus converges 7.5 times faster than emcee/AIES and 3.5 faster than emcee/DEMC on average. Finally, we found again that zeus is less sensitive to the specific initialisation of the walkers. In particular, out of the 40 tests conducted with different initialisation, zeus converged 38 times, emcee/DEMC 29 times and emcee/AIES 23 times prior to the predetermined maximum number of likelihood evaluations (i.e. 5 × 10 3 in this case). Detailed results about the values of the used metrics are shown in Table 2. The 1-D and 2-D marginal posterior distributions are shown in Figure 17, demonstrating the agreement between the three methods.

DISCUSSION
Following the analysis we conducted in Section 3 using the normal distribution there are two important questions that need to be answered about the initialisation of the walkers. First, how many walkers are necessary and, second, how to choose the initial positions of the walkers. Although there are many ways of answering those questions and there is no consistent solution that works for all target distributions, we will try to provide some general rules and heuristics to help ease the task of choosing the number and initial positions of the walkers for most cases.
Let us first discuss the effect of the number of the walkers on the general performance of zeus. Naively, one might expect that the minimum number of walkers should be + 1, where is the number of dimensions. However, the ensemble splitting technique, which was introduced in Section 3 to render the algorithm parallelisable, requires at least 2 × walkers in order to produce 2 linearly independent samples. If a smaller number is chosen then the walkers can be trapped in a lower-dimensional hyper-plane of the parameter space, being unable to sample properly and leading to erroneous results. Although there is no upper bound on the number of walkers, we recommend to use between two to four times the number of dimensions. The reason is that increasing the number of dimensions can increase the cost of the burn-in period as we explained in detail in Section 3. Ideally, one wants to use the minimum number (or close to that) of walkers until the burn-in period is over and then increase the number of walkers to rapidly produce a great number of independent samples. It is also worth noting that in cases in which either non-linear correlations or multiple modes are present it is recommended to use more walkers (e.g. 4-8 times the number of parameters for a bimodal target distribution).
As for the initialisation of the walkers, there are many ways to choose their starting positions ranging from prior sampling to more localised initial positions. Empirical tests indicate that the latter often outperforms the former (i.e. leads to shorter burn-in periods). That is not surprising since the total probability of a prior-sampled initialisation can be very small when the number of parameters is high. In particular we found that initialising the walkers from a tight region in parameter space (i.e. normal distribution with small variance) consistently leads to good performance. For low to moderate dimensional problems initialising the walkers from a tight ball around the Maximum A Posteriori (MAP) estimate can substantially reduce the burn-in period (Foreman-Mackey et al. 2013).
Finally, while emcee/AIES and emcee/DEMC can sample acceptably from most target distributions with 20, the efficient scaling of zeus with the number of parameters allows us to extend this range and efficiently test more complicated models (Karamanis & Beutler 2020). Like most gradient-free methods, zeus will fail to sample efficiently in very high dimensional problems in which = O (10 2 ). In such cases, more sophisticated algorithms (e.g. tempering, block updating, Hamiltonian dynamics etc.) need to be used (Robert et al. 2018).

CONCLUSIONS
The aim of this project was to develop a tool that could facilitate Bayesian parameter inference in computationally demanding astronomical analyses and tackle the challenges posed by the complexity of the models and data that are often used by astronomers. To this end, we introduced zeus, a parallel, general-purpose and gradientfree Python implementation of Ensemble Slice Sampling.
After introducing the method in Section 2, we thoroughly demonstrated its performance compared to two popular alternatives (i.e. emcee with affine-invariant ensemble sampling and differential evolution Metropolis) using a variety of artificial and realistic target distributions in Section 3. The artificial toy examples helped to shed light on the general behaviour of the samplers in target distributions characterised by linear and non-linear correlations as well as multimodal densities. When compared to emcee/AIES and emcee/DEMC in the problems of Baryon Acoustic Oscillation parameter inference and exoplanet radial velocity fitting, zeus consistently converges faster (i.e. its burn-in is shorter by a factor of at least 3), it is less sensitive to the initialisation of the walkers and generates substantially more independent samples per likelihood evaluation (i.e. approximately ×9 and ×29 speed-up compared to emcee/AIES in the cosmological and exoplanet examples, respectively).
We have shown that zeus performs similarly or better than existing MCMC methods in a range of problems. We hope that zeus will prove useful to the astronomical and cosmological community by complementing existing approaches and facilitating the study of novel models and data over the coming years. zeus is publicly available at https://github.com/minaskar/zeus with detailed documentation and examples that can be found at https: //zeus-mcmc.readthedocs.io.