s2fieldsofstudy:
Computer Science
Physics

text:
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.


Are Multilingual Models Effective in Code-Switching?

Multilingual language models have shown decent performance in multilingual and cross-lingual natural language understanding tasks. However, the power of these multilingual models in code-switching tasks has not been fully explored. In this paper, we study the effectiveness of multilingual language models to understand their capability and adaptability to the mixed-language setting by considering the inference speed, performance, and number of parameters to measure their practicality. We conduct experiments in three language pairs on named entity recognition and part-of-speech tagging and compare them with existing methods, such as using bilingual embeddings and multilingual meta-embeddings. Our findings suggest that pre-trained multilingual models do not necessarily guarantee high-quality representations on code-switching, while using meta-embeddings achieves similar results with significantly fewer parameters.

Introduction
Learning representation for code-switching has become a crucial area of research to support a greater variety of language speakers in natural language processing (NLP) applications, such as dialogue system and natural language understanding (NLU). Code-switching is a phenomenon in which a person speaks more than one language in a conversation, and its usage is prevalent in multilingual communities. Yet, despite the enormous number of studies in multilingual NLP, only very few focus on code-switching. Recently, contextualized language models, such as mBERT (Devlin et al., 2019) and XLM-R (Conneau et al., 2020) have achieved stateof-the-art results on monolingual and cross-lingual tasks in NLU benchmarks (Wang et al., 2018a;Hu et al., 2020;Wilie et al., 2020;. However, the effectiveness of these multilingual language models on code-switching tasks remains unknown. Several approaches have been explored in code-switching representation learning in NLU. Character-level representations have been utilized to address the out-of-vocabulary issue in codeswitched text (Winata et al., 2018c;Wang et al., 2018b), while external handcrafted resources such as gazetteers list are usually used to mitigate the low-resource issue in code-switching (Aguilar et al., 2017;Trivedi et al., 2018); however, this approach is very limited because it relies on the size of the dictionary and it is language-dependent. In another line of research, meta-embeddings have been used in code-switching by combining multiple word embeddings from different languages (Winata et al., 2019a,b). This method shows the effectiveness of mixing word representations in closely related languages to form languageagnostic representations, and is considered very effective in Spanish-English code-switched named entity recognition tasks, and significantly outperforming mBERT (Khanuja et al., 2020) with fewer parameters.
While more advanced multilingual language models (Conneau et al., 2020) than multilingual BERT (Devlin et al., 2019) have been proposed, their effectiveness is still unknown in codeswitching tasks. Thus, we investigate their effectiveness in the code-switching domain and compare them with the existing works. Here, we would like to answer the following research question, "Which models are effective in representing code-switching text, and why?." In this paper, we evaluate the representation quality of monolingual and bilingual word embeddings, multilingual meta-embeddings, and multilingual language models on five downstream tasks on named entity recognition (NER) and part-ofspeech tagging (POS) in Hindi-English, Spanish-English, and Modern Standard Arabic-Egyptian. We study the effectiveness of each model by considering three criteria: performance, speed, and the number of parameters that are essential for practical applications. Here, we set up the experimental setting to be as language-agnostic as possible; thus, it does not include any handcrafted features.
Our findings suggest that multilingual pretrained language models, such as XLM-R BASE , achieves similar or sometimes better results than the hierarchical meta-embeddings (HME) (Winata et al., 2019b) model on code-switching. On the other hand, the meta-embeddings use word and subword pre-trained embeddings that are trained using significantly less data than mBERT and XLM-R BASE and can achieve on par performance to theirs. Thus, we conjecture that the masked language model is not be the best training objective for representing code-switching text. Interestingly, we found that XLM-R LARGE can improve the performance by a great margin, but with a substantial cost in the training and inference time, with 13x more parameters than HME-Ensemble for only around a 2% improvement. The main contributions of our work are as follows: • We evaluate the performance of word embeddings, multilingual language models, and multilingual meta-embeddings on code-switched NLU tasks in three language pairs, Hindi-English (HIN-ENG), Spanish-English (SPA-ENG), and Modern Standard Arabic-Egyptian (MSA-EA), to measure their ability in representing code-switching text.
• We present a comprehensive study on the effectiveness of multilingual models on a variety of code-switched NLU tasks to analyze the practicality of each model in terms of performance, speed, and number of parameters.
• We further analyze the memory footprint re-quired by each model over different sequence lengths in a GPU. Thus, we are able to understand which model to choose in a practical scenario.

Representation Models
In this section, we describe multilingual models that we explore in the context of code-switching. Figure 1 shows the architectures for a word embeddings model, a multilingual language model, and the multilingual meta-embeddings (MME), and HME models.

FastText
In general, code-switching text contains a primary language the matrix language (ML)) as well as a secondary language (the embedded language (EL)).
To represent code-switching text, a straightforward idea is to train the model with the word embeddings of the ML and EL from FastText . Code-switching text has many noisy tokens and sometimes mixed words in the ML and EL that produce a "new word", which leads to a high number of out-of-vocabulary (OOV) tokens. To solve this issue, we utilize subword-level embeddings from FastText  to generate the representations for these OOV tokens. We conduct experiments on two variants of applying the word embeddings to the code-switching tasks: FastText (ML) and FastText (EL), which utilize the word embeddings of ML and EL, respectively.

MUSE
To leverage the information from the embeddings of both the ML and EL, we utilize MUSE (Lample et al., 2018) to align the embeddings space of the ML and EL so that we can inject the information of the EL embeddings into the ML embeddings, and vice versa. We perform alignment in two directions: (1) We align the ML embeddings to the vector space of the EL embeddings (denoted as MUSE (ML → EL)); (2) We conduct the alignment in the opposite direction, which aligns the EL embeddings to the vector space of the ML embeddings (denoted as MUSE (EL → ML)). After the embeddings alignment, we train the model with the aligned embeddings for the code-switching tasks.

Multilingual Pre-trained Models
Pre-trained on large-scale corpora across numerous languages, multilingual language models (Devlin et al., 2019;Conneau et al., 2020) possess the ability to produce aligned multilingual representations for semantically similar words and sentences, which brings them advantages to cope with codemixed multilingual text.

Multilingual BERT
Multilingual BERT (mBERT) (Devlin et al., 2019), a multilingual version of the BERT model, is pretrained on Wikipedia text across 104 languages with a model size of 110M parameters. It has been shown to possess a surprising multilingual ability and to outperform existing strong models on multiple zero-shot cross-lingual tasks (Pires et al., 2019;Wu and Dredze, 2019). Given its strengths in handling multilingual text, we leverage it for code-switching tasks.

XLM-RoBERTa
XLM-RoBERTa (XLM-R) (Conneau et al., 2020) is a multilingual language model that is pre-trained on 100 languages using more than two terabytes of filtered CommonCrawl data. Thanks to the largescale training corpora and enormous model size (XLM-R BASE and XLM-R LARGE have 270M and 550M parameters, respectively), XLM-R is shown to have a better multilingual ability than mBERT, and it can significantly outperform mBERT on a variety of cross-lingual benchmarks. Therefore, we also investigate the effectiveness of XLM-R for code-switching tasks.

Char2Subword
Char2Subword introduces a character-to-subword module to handle rare and unseen spellings by training an embedding lookup table (Aguilar et al., 2020b). This approach leverages transfer learning from an existing pre-trained language model, such as mBERT, and resumes the pre-training of the upper layers of the model. The method aims to increase the robustness of the model to various typography styles.

Multilingual Meta-Embeddings
The MME model (Winata et al., 2019a) is formed by combining multiple word embeddings from different languages. Let's define w to be a sequence of words with n elements, where w = [w 1 , . . . , w n ].
First, a list of word-level embedding layers is used E (w) i to map words w into embeddings x i . Then, the embeddings are combined using one out of the following three methods: concat, linear, and selfattention. We briefly discuss each method below.
Concat This method concatenates word embeddings by merging the dimensions of word representations into higher-dimensional embeddings. This is one of the simplest methods to join all embeddings without losing information, but it requires a larger activation memory than the linear method. (1) Linear This method sums all word embeddings into single word embeddings with equal weight without considering each embedding's importance. The method may cause a loss of information and may generate noisy representations. Also, though it is very efficient, it requires an additional layer to project all embeddings into a single-dimensional space if one embedding is larger than another.
Self-Attention This method generates a metarepresentation by taking the vector representation from multiple monolingual pre-trained embeddings in different subunits, such as word and subword. It applies a projection matrix W j to transform the dimensions from the original space x i,j ∈ R d to a new shared space x i,j ∈ R d . Then, it calculates attention weights α i,j ∈ R d with a non-linear scoring function φ (e.g., tanh) to take important information from each individual embedding x i,j . Then, MME is calculated by taking the weighted sum of the projected embeddings x i,j :

Hierarchical Meta-Embedings
The HME method combines word, subword, and character representations to create a mixture of embeddings (Winata et al., 2019b). It generates multilingual meta-embeddings of words and subwords, and then, concatenates them with character-level embeddings to generate final word representations. HME combines the word-level, subword-level, and character-level representations by concatenation, and randomly initializes the character embeddings. During the training, the character embeddings are trainable, while all subword and word embeddings remain fixed.

HME-Ensemble
The ensemble is a technique to improve the model's robustness from multiple predictions. In this case, we train the HME model multiple times and take the prediction of each model. Then, we compute the final prediction by majority voting to achieve a consensus. This method has shown to be very effective in improving the robustness of an unseen test set (Winata et al., 2019c). Interestingly, this method is very simple to implement and can be easily spawned in multiple machines, as in parallel processes.

Experiments
In this section, we describe the details of the datasets we use and how the models are trained.

Datasets
We evaluate our models on five downstream tasks in the LinCE Benchmark (Aguilar et al., 2020a). We choose three named entity recognition (NER) tasks, Hindi-English (HIN-ENG) , Spanish-English (SPA-ENG) (Aguilar et al., 2018) and Modern Standard Arabic (MSA-EA) (Aguilar et al., 2018), and two part-of-speech (POS) tagging tasks, Hindi-English (HIN-ENG) (Singh et al., 2018b) and Spanish-English (SPA-ENG) (Soto and Hirschberg, 2017). We apply Roman-to-Devanagari transliteration on the Hindi-English datasets since the multilingual models are trained with data using that form.

Experimental Setup
We describe our experimental details for each model.

Scratch
We train transformer-based models without any pretraining by following the mBERT model structure, and the parameters are randomly initialized, including the subword embeddings. We train transformer models with four and six layers with a hidden size of 768. This setting is important to measure the effectiveness of pre-trained multilingual models. We start the training with a learning rate of 1e-4 and an early stop of 10 epochs.

Word Embeddings
We use FastText embeddings  to train our transformer models. The model consists of a 4-layer transformer encoder with four heads and a hidden size of 200. We train a transformer followed by a Conditional Random Field (CRF) layer (Lafferty et al., 2001).
The model is trained by starting with a learning rate of 0.1 with a batch size of 32 and an early stop of 10 epochs. We also train our model with only ML and EL embeddings. We freeze all embeddings and only keep the classifier trainable. We leverage MUSE (Lample et al., 2018) to align the embeddings space between the ML and EL. MUSE mainly consists of two stages: adversarial training and a refinement procedure. For all alignment settings, we conduct the adversarial training using the SGD optimizer with a starting learning rate of 0.1, and then we perform the refinement procedure for five iterations using the Procrustes solution and CSLS (Lample et al., 2018). After the alignment, we train our model with the aligned word embeddings (MUSE (ML → EL) or MUSE (EL → ML)) on the code-switching tasks.

Pre-trained Multilingual Models
We use pre-trained models from Huggingface. 1 On top of each model, we put a fully-connected layer classifier. We train the model with a learning rate between [1e-5, 5e-5] with a decay of 0.1 and a batch size of 8. For large models, such as XLM-R LARGE and XLM-MLM LARGE , we freeze the embeddings layer to fit in a single GPU.

Multilingual Meta-Embeddings (MME)
We use pre-trained word embeddings to train our MME. Table 2 shows the embeddings used for each dataset. We freeze all embeddings and train a transformer classifier with the CRF. The transformer classifier consists of a hidden size of 200, a head of 4, and 4 layers. All models are trained with a learning rate of 0.1, an early stop of 10 epochs, and a batch size of 32. We follow the implementation from the code repository. 2 Table 2 shows the list of word embeddings used in MME.

Hierarchical Meta-Embeddings (HME)
We train our HME model using the same embeddings as MME and pre-trained subword embeddings from Heinzerling and Strube (2018). The subword embeddings for each language pair are shown in Table 3. We freeze all word embeddings and subword embeddings, and keep the character embeddings trainable.

Other Baselines
We compare the results with Char2subword and mBERT (cased) from Aguilar et al. (2020b). We also include the results of English BERT provided by the organizer of the LinCE public benchmark leaderboard (accessed on March 12nd, 2021). 3

LinCE Benchmark
We evaluate all the models on the LinCE benchmark, and the development set results are shown in Table 4. As expected, models without any pretraining (e.g., Scratch (4L)) perform significantly worse than other pre-trained models. Both Fast-Text and MME use pre-trained word embeddings, but MME achieves a consistently higher F1 score than FastText in both NER and POS tasks, demonstrating the importance of the contextualized selfattentive encoder. HME further improves on the F1 score of the MME models, suggesting that encoding hierarchical information from sub-word level, word level, and sentence level representations can improve code-switching task performance. Comparing HME with contextualized pre-trained mul-   tilingual models such as mBERT and XLM-R, we find that HME models are able to obtain competitive F1 scores while maintaining a 10x smaller model sizes. This result indicates that pre-trained multilingual word embeddings can achieve a good balance between performance and model size in code-switching tasks. Table 5 shows the models' performance in the LinCE test set. The results are highly correlated to the results of the development set. XLM-R LARGE achieves the best-averaged performance, with a 13x larger model size compared to the HME-Ensemble model.

Model Effectiveness and Efficiency
Performance vs. Model Size As shown in Figure 2, the Scratch models yield the worst average score, at around 60.93 points. With the smallest pre-trained embedding model, FastText, the model performance can improve by around 10 points compared to the Scratch models and they only have 10M parameters on average. On the other hand, the MME models, which have 31.6M parameters on average, achieve similar results to the mBERT models, with around 170M parameters. Interestingly, adding subwords and character embeddings to MME, such as in the HME models, further im-  proves the performance of the MME models and achieves a 81.60 average score, similar to that of the XLM-R BASE and XLM-MLM LARGE models, but with less than one-fifth the number of parameters, at around 42.25M. The Ensemble method adds further performance improvement of around 1% with an additional 2.5M parameters compared to the non-Ensemble counterparts.
Inference Time To compare the speed of different models, we use generated dummy data with various sequence lengths, [16,32,64,128,256,512,1024,2048,4096]. We measure each model's inference time and collect the statistics of each model at one particular sequence length by running the model 100 times. The experiment is performed on a single NVIDIA GTX1080Ti GPU. We do not include the pre-processing time in our analysis. Still, it is clear that the pre-processing time for meta-embeddings models is longer than for other models as pre-processing requires a tokenization step to be conducted for the input multiple times with different tokenizers. The sequence lengths are counted based on the input tokens of each model. We use words for the MME and HME models, and subwords for other models.
The results of the inference speed test are shown in Figure 3. Although all pre-trained contextualized language models yield a very high validation score, these models are also the slowest in terms of inference time. For shorter sequences, the HME model performs as fast as the mBERT and XLM-R BASE models, but it can retain the speed as the sequence length increases because of the smaller model dimension in every layer. The FastText, MME, and Scratch models yield a high throughput in shortsequence settings by processing more than 150  samples per second. For longer sequences, the same behavior occurs, with the throughput of the Scratch models reducing as the sequence length increases, even becoming lower than that of the HME model when the sequence length is greater than or equal to 256. Interestingly, for the FastText, MME, and HME models, the throughput remains steady when the sequence length is less than 1024, and it starts to decrease afterwards.

Memory Footprint
We record the memory footprint over different sequence lengths, and use the same setting for the FastText, MME, and HME models as in the inference time analysis. We record the size of each model on the GPU and the size of the activation after performing one forward operation to a single sample with a certain sequence length. The result of the memory footprint analysis for a sequence length of 512 is shown in Table  6. Based on the results, we can see that metaembedding models use a significantly smaller memory footprint to store the model and activation memory. For instance, the memory footprint of the HME

Related Work
Transfer Learning on Code-Switching Previous works on code-switching have mostly focused on combining pre-trained word embeddings with trainable character embeddings to represent noisy mixed-language text (Trivedi et al., 2018;Wang et al., 2018b;Winata et al., 2018c). Winata et al.
(2018a) presented a multi-task training framework to leverage part-of-speech information in a language model. Later, they introduced the MME in the code-switching domain by combining multiple word embeddings from different languages (Winata et al., 2019a). MME has since also been applied to Indian languages (Priyadharshini et al., 2020;Dowlagar and Mamidi, 2021). Meta-embeddings have been previously explored in various monolingual NLP tasks (Yin and Schütze, 2016;Muromägi et al., 2017;Coates and Bollegala, 2018;Kiela et al., 2018). Winata et al. (2019b) introduced hierarchical meta-embeddings by leveraging subwords and characters to improve the codeswitching text representation. Pratapa et al. (2018b) propose to train skip-gram embeddings from synthetic code-switched data generated by Pratapa  Gupta et al. (2020) proposed a generative-based model for augmenting codeswitching data from parallel data. Recently, Aguilar et al. (2020b) proposed the Char2Subword model, which builds representations from characters out of the subword vocabulary, and they used the module to replace subword embeddings that are robust to misspellings and inflection that are mainly found in a social media text. Khanuja et al. (2020) explored fine-tuning techniques to improve mBERT for code-switching tasks, while Winata et al. (2020) introduced a meta-learning-based model to leverage monolingual data effectively in code-switching speech and language models.

Conclusion
In this paper, we study multilingual language models' effectiveness so as to understand their capability and adaptability to the mixed-language setting. We conduct experiments on named entity recognition and part-of-speech tagging on various language pairs. We find that a pre-trained multilingual model does not necessarily guarantee highquality representations on code-switching, while the hierarchical meta-embeddings (HME) model achieve similar results to mBERT and XLM-R BASE but with significantly fewer parameters. Interestingly, we find that XLM-R LARGE has better performance by a great margin, but with a substantial cost in the training and inference time, using 13x more parameters than HME-Ensemble for only a 2% improvement. SAR Government, and School of Engineering Ph.D. Fellowship Award, the Hong Kong University of Science and Technology, and RDC 1718050-0 of EMOS.AI.