Fearless Stochasticity in Expectation Propagation

Authors: jonathan  so
richard e turner
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 17:10:33

Abstract: Expectation propagation (EP) is a family of algorithms for performing approximate inference in probabilistic models. The updates of EP involve the evaluation of moments—expectations of certain functions—which can be estimated from Monte Carlo (MC) samples. However, the updates are not robust to MC noise when performed naively, and various prior works have attempted to address this issue in different ways. In this work, we provide a novel perspective on the moment-matching updates of EP; namely, that they perform natural-gradient-based optimisation of a variational objective. We use this insight to motivate two new EP variants, with updates that are particularly well-suited to MC estimation. They remain stable and are most sample-efficient when estimated with just a single sample. These new variants combine the benefits of their predecessors and address key weaknesses. In particular, they are easier to tune, offer an improved speed-accuracy trade-off, and do not rely on the use of debiasing estimators. We demonstrate their efficacy on a variety of probabilistic inference tasks.

Main paper: 
 introduction: 
 Expectation propagation (EP) [37, 40] is a family of algorithms that is primarily used for performing approximate inference in probabilistic models [9, 10, 12, 13, 16, 18, 20, 19, 21, 22, 23, 24, 28, 29, 30, 32, 35, 36, 38, 39, 41, 44, 47, 48, 50, 52, 53] , although it can be used more generally for approximating certain kinds of functions and their integrals [14] .EP involves the evaluation of moments-expectations of certain functions-under distributions that are derived from the model of interest. EP is usually applied to models for which these moments have convenient closed-form expressions or can be accurately estimated using deterministic methods. Moments can also be estimated using Monte Carlo (MC) samples, significantly expanding the set of models EP can be applied to. However, the updates of EP are not robust to MC noise when performed naively, and various prior works have attempted to address this issue in different ways [17, 48, 52] .In this work we provide a novel perspective on the moment-matching updates of EP; namely, that they perform natural-gradient-based optimization of a variational objective (Section 3). We use this insight to motivate two new EP variants, EP-η (Section 3.2) and EP-µ (Section 3.3), with updates that are particularly well-suited to MC estimation, remaining stable and being most sample-efficient when estimated with just a single sample. These new variants combine the benefits of their predecessors and address key weaknesses. In particular, they are easier to tune, offer an improved speed-accuracy trade-off, and do not rely on the use of debiasing estimators. We demonstrate their efficacy on a variety of probabilistic inference tasks (Section 4).Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes]Justification: The novel perspective on EP is given in Section 3, with proof in Appendix D. The new EP variants are also introduced in Section 3, and we demonstrate the stated advantages in the results of Section 4, and Appendices K and L. The abstract and introduction also make clear that our proposed methods are primarily of use when the updates of EP are to be estimated with MC samples.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
background: 
 In this section, we first introduce the problem setting. We then give an overview of EP, followed by a discussion of issues related to sampled estimation of EP updates.Let F be the tractable, minimal exponential family of distributions (see Appendix A), defined by the statistic function s(.) with respect to base measure ν(.). Let Ω, M, A(.) denote the natural domain, mean domain and log-partition function of F, respectively, and A * (.) the convex dual of A(.).Let p 0 be the member of F with natural parameter η 0 , so that p 0 (z) = exp(η ⊤ 0 s(z) -A(η 0 )) 1 , and assume that a distribution of interest p * , the target distribution, has a density of the form p * (z) ∝ p 0 (z) m i=1 exp ℓ i (z) .(EQUATIONIn Bayesian inference settings, p * would be the posterior distribution over parameters z given some observed data D, where p 0 may correspond to a prior distribution, and {ℓ i (z)} i to log-likelihood terms given some partition of D. 2 However, we consider the more general setting in which the target distribution is simply a product of factors. Note that we can assume form (1) without loss of generality if we allow p 0 to be improper. The inference problem typically amounts to computing quantities derived from the normalised density p * (z), such as samples, summary statistics, or expectations of given functions. In some special cases, these quantities can be computed exactly, but this is not feasible in general, and approximations must be employed.Given a target distribution density of form (1) , EP aims to find an approximation p ∈ F such that p(z) ∝ p 0 (z)EQUATIONBy assumption F is tractable, and so provided that (η 0 + i λ i ) ∈ Ω, p is a tractable member of F. Each factor exp(λ ⊤ i s(z)) is known as a site potential, and can be roughly interpreted as a F-approximation to the i-th target factor, exp(ℓ i (z)). 3 λ i is known as the i-th site parameter.Variational problem EP, and several of its variants [17, 24, 37] , can be viewed as solving a variational problem which we now introduce following the exposition of Hasenclever et al. [17] .Let the i-th locally extended family, denoted F i , be the exponential family defined by the statistic function s i (z) = (s(z), ℓ i (z)) with respect to base measure ν(.). Let Ω i , M i and A i (.) denote the natural domain, mean domain, and log partition function of F i , respectively. Note that a member of F i roughly corresponds to a distribution whose density is the (normalised) product of a member of F with the i-th target factor raised to a power. Unlike F, we do not assume F i is minimal. EQUATIONAt a solution to (3), the EP approximation is given by (2) . The hyperparameters {β i } i control the characteristics of the approximation, and correspond to the power parameters of power EP.EP updates EP [37, 40] , power EP [34] and double-loop EP [24] , can all be viewed as alternating between some number (≥ 1) of inner updates to decrease L with respect to {λ i } i , with an outer update to increase L with respect to θ [17, 25] . EP is not typically presented in this way, but by doing so we will be able to present a unified algorithm that succinctly illustrates the relationship between the different variants. We show equivalence with the conventional presentation of EP in Appendix B. The inner and outer updates are given by Inner update:EQUATIONOuter update:EQUATIONwhere p i ∈ F i denotes the i-th tilted distribution, with densityEQUATIONThe hyperparameter α controls the level of damping, which is used to aid convergence. An undamped inner update (with α = 1) can be seen as performing moment matching between p and p i . The expectation in (4) is most often computed analytically, or estimated using deterministic numerical methods. It can also be estimated by sampling from p i , however, we will later show that the resulting stochasticity can lead to a biased and unstable procedure. The inner updates can be performed either serially or in parallel (over i). In this work we assume they are applied in parallel, but the ideas presented can easily be extended to the serial case. Heskes and Zoeter [25] showed that (4) follows a decrease direction in L with respect to λ i , and so is guaranteed to decrease L when α is small enough. See Appendix C for a derivation of the inner update. When {λ i } i are at a partial minimum of L (for fixed θ), the outer update performs an exact partial maximisation of the primal form of the variational problem -see Hasenclever et al. [17] for details. All of the EP variants considered in this paper differ only in their handling of the inner minimisation, and this is our focus.Unified EP algorithm The double-loop EP algorithm of Heskes and Zoeter [24] repeats the inner update to convergence before performing each outer update, which ensures convergence of the overall procedure. The usual presentation of EP combines (4) and ( 5) into a single update in λ i (see Appendix B), however, Jylänki et al. [30] observed that it can also be viewed as performing double-loop EP with just a single inner update per outer update (which is not guaranteed to converge in general). By taking this view, we are we are able to present EP, power EP, and their double-loop counterparts as a single algorithm, presented in Algorithm 1. We do so primarily to illustrate how these variants are related to one another, and to the new variants of Section 3. Algorithm 1 EP (β i =1,n inner =1), power EP (β i ̸ =1,n inner =1), and their double-loop variants (n inner >1)Require: F, η 0 , {β i } i , {ℓ i (z)} i , {λ i } i , n inner , αwhile not converged do θ ← η 0 + j λ j for 1 to n inner do for i = 1 to m in parallel doStochastic estimation of this expectation can lead to biased and unstable updates.λ i ← λ i -α(η 0 + j λ j -∇A * (E pi(z) [s(z)])) return {λ i } i Note that dependence on {β i } i comes through the definition of p i (.). We can think of "exact" double-loop EP as corresponding to n inner = ∞, with the inner loop exiting once some convergence criterion has been satisfied. In practice, a truncated inner loop is often used by setting n inner to some small number [17, 30] . Upon convergence of Algorithm 1, the approximation p(z) ≈ p * (z) is given by (2) . Going forwards we will refer to the family of algorithms encompassed by Algorithm 1 simply as EP.Stochastic moment estimation EP is typically applied to models for which the updates either have closed-form expressions, or can be accurately estimated using deterministic numerical methods. However, as update (4) only depends on the target distribution through expectations under the tilted distributions, this suggests that updates could also be performed using sampled estimates of those expectations. By using MC methods to estimate the tilted distribution moments, EP can be used in a black-box manner, dramatically expanding the set of models it can be applied to. Instead of performing a single large sampling task-as would be required by applying MC methods directly-EP can instead solve several simpler ones, gaining significant computational advantages [17, 48, 52] . Unfortunately, when performed naively, update (4) is not robust to MC noise. This is because the moment estimates are converted to the natural (site) parameter space by mapping through ∇A * (.), which is not linear in general, and so MC noise in the estimates leads to biased updates of λ i .Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: Limitations are discussed in Section 5, and computational efficiency is also discussed separately in Appendix I.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
fearlessly stochastic ep algorithms: 
 In this section we show an equivalence between the moment-matching updates of EP and natural gradient descent (NGD). We use this view to motivate two new EP variants which have several advantages when updates are estimated using samples. We conclude with a review of related work.Several EP variants-EP [40, 37] , power EP, double loop EP [34] , stochastic natural gradient EP [17] , and the new ones to follow-differ only in how they perform the inner optimisation in (3) with respect to {λ i } i . The optimisation can be viewed as one with respect to the m + 1 distributions p, p 1 , . . . , p m , which are jointly parameterised by {λ i } i [17] . Natural gradient descent (NGD) [1] , which involves preconditioning a gradient with the inverse Fisher information matrix (FIM) of a distribution, is an effective method for optimising parameters of probability distributions. It would therefore seem desirable to apply NGD to the inner optimisation, yet doing so is not straightforward. While the FIM can be naturally extended to the product of statistical manifolds that is the solution space, computing and inverting it will not generally be tractable, making (the natural extension of) NGD in this space infeasible.(t) ∈ F with natural parameter η i (t) (λ i ) = η 0 + j λ j (t) + α -1 (λ i -λ i (t)). The t superscript on the site parameters indicates that they are fixed for the current iteration, and so the distribution is fully parameterised by λ i . Note that when λ j = λ j (t) ∀ j, we have pi (t) = p. Proposition 1, which we prove in Appendix D, states that the moment-matching updates of EP can be viewed as performing NGD in L with respect to mean parameters of pi (t) .We will make use of the following properties: for an exponential family with log partition function A(.), the gradient ∇A(.) provides the map from natural to mean parameters, and this mapping is one-to-one for minimal families, with the inverse given by ∇A * (.). See Appendix A for details. Proposition 1. For α > 0, the moment-matching update of EP (4) is equivalent to performing an NGD step in L with respect to the mean parameters of pi (t) with step size α -1 . That is, forµ i = E pi (t) (z) [s(z)],and Fi (t) (µ i ) the FIM of pi (t) with respect to µ i , we haveEQUATIONNote that the right hand side of (7) just maps update (4) to mean parameters of pi (t) . 4 We discussed in Section 2 that noise in the moment estimates results in bias in the update to λ i of EP, due to the noise being passed through the nonlinear map ∇A * (.). This alone would not be a problem if the updates followed unbiased descent direction estimates in some other fixed parameterisation. Let η i = ∇A * (µ i ), then, from the definition of L and µ i , we haveEQUATIONand so when E pi(z) [s(z)] is estimated with unbiased samples, the NGD update (7) in µ i is also unbiased. However, a sequence of updates introduces bias. This is because the parameterisation changes from one update to the next, with the mean parameters of pi (t) mapped to those of pi (t+1) byEQUATIONwhere (η i (t) ) -1 (.) is the inverse of η i (t) (.), and η i (t+1) (.) is defined using the site parameters after the update at time t. This map is not linear in general, and so an unbiased, noisy update of µ i at one time step, results in bias when mapped to the next. This bias can make the EP updates unstable in the presence of MC noise, leading previous work to use methods, specific to F, for obtaining approximately unbiased natural parameter estimates from samples. However, even with these debiasing methods, a relatively large number of samples is still needed in practice [48, 52] . Here we have assumed parallel application of (4), but the bias at site i is induced whenever the parameters for other sites change from one update to the next, which occurs in both parallel and serial settings.We now use the natural gradient interpretation of the EP updates to motivate two new variants that are far more robust to MC noise. In particular, they are stable, and most sample-efficient, when updates are estimated with just a single sample. Furthermore, unlike methods that rely on deibasing estimators, they do not require sample thinning (see Section 3.4) . This largely eliminates two hyperparameters for the practitioner. While practical considerations, such as ensuring efficient use of parallel hardware, may justify using more than one sample per update, this can in principle be tuned a priori, much like the batch size hyperparameter of stochastic gradient descent. [37] . EP (naive) uses maximum likelihood estimation for the updates, and EP (debiased) uses the estimator of Xu et al. [52] .Step size corresponds to α for EP, and ϵ for EP-µ and EP-η. Only EP-η and EP-µ can perform 1-sample updates, hence the other traces are not visible. The left panel shows the expected decrease in L after 100/n samp steps. Performing e.g. 100× 1-sample steps, or 10× 10-sample steps, achieves a much larger decrease in L than a single 100-sample step. The right panel shows the magnitude of the bias in λ i after a single parallel update, averaged over all sites and dimensions. The bias of EP-µ shrinks far faster as the step size decreases than that of EP. EP-η is always unbiased and so is not visible.We showed in the previous section that bias is introduced into the sequence of NGD updates due to a nonlinear map from one parameterisation to the next. If the map were affine, the sequence would remain unbiased. This can be achieved by performing NGD with respect to η i , the natural parameters of pi (t) . Then, the map from one parameterisation to the next (by equating λ i ) is given by η i (t+1) • (η i (t) ) -1 , which is linear. The NGD direction in L with respect to pi (t) and η i is given byEQUATIONSee Appendix E for a derivation. Note that α plays a similar role as the NGD step size in this parameterisation, but it also affects pi (t) , and so to obtain a more faithful NGD interpretation we fix α = 1 and introduce an explicit NGD step size ϵ to decouple from the effect on pi (t) . The resulting update can be expressed directly as an update in λ i , by applying (η i (t) ) -1 to the updated η i , givingEQUATIONWe can use automatic differentiation to efficiently compute (11) , by recognising the second term as a Jacobian-vector product (JVP) with respect to η i (µ i ) = ∇A * (µ i ). 5 In summary, by performing NGD with respect to natural parameters of pi (t) , instead of mean parameters, the resulting sequence of updates is unbiased in λ i . We call the resulting procedure-which is given in Algorithm 2-EP-η, to emphasise that we have simply changed the NGD parameterisation of EP. The unbiased updates of EP-η allow it to be more sample-efficient than EP by using a smaller number of samples per iteration. The left panel of Figure 1 demonstrates this effect in a stochastic version of the clutter problem [37] .The bias introduced by the updates of EP cannot be mitigated by reducing α (increasing the damping) without also proportionally sacrificing the amount of progress made by each update. More concretely, let the bias in dimension d of µ i (t+1) , after an update at time t, be defined as E[µ i (t+1) -μi (t+1) ] d , where μi (t+1) is the value of µ i (t+1) after a noise-free update, and expectation is taken over the sampling distributions of all parallel updates at time t. Then, Proposition 2 below, which we prove in Appendix G, summarises the effect of decreasing α on EP. Proposition 2. After update (4) is executed in parallel over i, as α → 0 + , both the expected decrease in L, and the bias E[µ i (t+1) -μi (t+1) ] d , are O(α) for all d.Proposition 1 at the beginning of this section states that update (4) can be viewed as performing an NGD step with respect to µ i with a step size of α -1 . However, changing α also has the effect of changing the definition of pi (t) . It is then natural to wonder what happens if we fix α = 1 and introduce an explicit step size for NGD, ϵ, resulting in the updateEQUATIONSee Appendix F for a derivation. It turns out that in doing so, when we decrease ϵ, the bias shrinks far faster than the expected decrease in L. Proposition 3. After update (12) is executed in parallel over i, as ϵ → 0 + , the expected decrease in L is O(ϵ), and the biasE µ i (t+1) -μi (t+1) d is O(ϵ 2 ), for all d.Proposition 3, which we prove in Appendix G, tells us that by using update (12) , we can reduce the bias to arbitrarily small levels while still making progress in decreasing L. Update (12) can also be expressed directly as an update in λ i by applying (η (t) ) -1 • ∇A * , givingEQUATIONwhich has the simple interpretation of performing EP, but with damping of the moments instead of the site (natural) parameters. In summary, if we perform NGD with respect to the same mean parameterisation as EP, but treat the NGD step size as a free parameter, ϵ, we can obtain updates that are approximately unbiased while still making progress. We call this variant EP-µ, again to indicate that we are simply performing the NGD update of EP in the mean parameter space. 6 The resulting procedure is also given by Algorithm 2.Algorithm 2 EP-η and EP-µ (differences with Algorithm 1 are highlighted in green)Require: F, η 0 , {β i } i , {ℓ i (z)} i , {λ i } i , n inner , ϵ while not converged do θ ← η 0 + j λ j for 1 to n inner do for i = 1 to m in parallel do update λ i using (11) for EP-η, or (13) for EP-µ return {λ i } iThe computational cost of EP-µ is lower than that of EP-η because it does not require JVPs through ∇A * (.). 7 The drawback is that the updates of EP-µ do still retain some bias, however, we find that the bias of EP-µ typically has negligible impact on its performance relative to EP-η. This is evident in the clutter problem of Minka [37] , as demonstrated in the left panel of Figure 1 , as well as in the larger scale experiments of Section 4. The right panel of Figure 1 illustrates how quickly the bias in λ i shrinks as the step size of EP-µ is decreased. Note that the results of this subsection are stated in terms of bias in µ i , but it is straightforward to show that equivalent results also hold for λ i using Taylor series arguments.The link between natural gradients and moment matching in an exponential family is well known [6] , but the connection with EP shown here -which relies on identification of the distribution pi(t)and parameterisation η i (t) -is new, to the best of our knowledge. Bui et al. [11] showed that the updates of (power) EP are equivalent to performing NGD in a local variational free energy objective when taking the limit of the power parameter as β i → ∞, and Wilkinson et al. [51] showed that this extends to the (global) variational free energy objective for models with a certain structure. Our result is different in that it relates to NGD of the EP variational objective, and applies for all values of β i .EP was combined with Markov chain Monte Carlo (MCMC) moment estimates by Xu et al. [52] in a method called sampling via moment sharing (SMS), and later by Vehtari et al. [48] in the context of hierarchical models. In both works, when F was multivariate normal (MVN) -arguably the most useful case-the authors used an estimator for the updates that is unbiased when pi is also MVN. Although this estimator is only approximately unbiased in general, it can help to stabilise the updates significantly. Even so, it is necessary to use a relatively large number of samples for the updates, with several hundred or more being typical. Using many samples per update is inefficient, as the update direction can change from one iteration to the next. Furthermore, such estimators typically rely on a known number of independent samples. Samples drawn using MCMC methods are generally autocorrelated, necessitating the use of sample thinning before the estimators can be applied. This adds another hyperparameter to the procedure-the thinning ratio-and it is also inefficient due to the discarding of samples. The practitioner is required to choose the number of samples, the thinning ratio, and the amount of damping, all of which affect the accuracy, stability, and computational efficiency of the procedure. There is no easy way to choose these a priori, forcing the practitioner to either set them conservatively (favouring accuracy and stability), or to find appropriate settings by trial and error, both of which are likely to expend time and computation unnecessarily.The stochastic natural gradient EP (SNEP) method of Hasenclever et al. [17] , which is closely related to our work, also optimises the EP variational objective using a form of NGD. The SNEP updates are unbiased in the presence of MC noise, allowing them to be performed with as little as one sample, without relying on F-specific debiasing estimators, or the sample thinning that would entail. In SNEP, NGD is performed with respect to mean parameters of the site potentials, which are treated as bona fide distributions in F. In contrast, we showed that EP can already be viewed as performing NGD, but with respect to the distributions {p i (t) } i , and our new variants, EP-η and EP-µ, are able to gain the same advantages as SNEP, but using the same distributions for NGD as EP. Hasenclever et al. [17] showed that in some settings, SNEP can obtain accurate point estimates fairly rapidly. However, we find that it typically converges far slower than both EP and our new variants, consistent with findings in Vehtari et al. [48] . We argue that this is because the site potentials (when considered as distributions) bear little resemblance with the distributions that are ultimately being optimised, and so their geometry is largely irrelevant for the optimisation of L. In contrast, the geometry of p i (t) is closely related to that of L, as we show in Appendix H.We note that Xu et al. [52] and Hasenclever et al. [17] also proposed methods for performing updates in an asynchronous fashion, significantly reducing the frequency and cost of communication between nodes in distributed settings. In principle, these methods could be combined with those presented in this paper, but we do not consider them further here.some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
evaluation: 
 In this section we demonstrate the efficacy of EP-η and EP-µ on a variety of probabilistic inference tasks. In each experiment, the task was to perform approximate Bayesian inference of unobserved parameters z, given some observed data D. All of the models in these experiments followed the same general structure, consisting of a minimal exponential family prior over z, p 0 , and a partition of the data {D i } m i , where each block D i depends on both z and an additional vector of local latent variables w i that also depend on z. That is, the joint density has the formEQUATIONwhere p 0 ∈ F. This structure is shown graphically in Appendix J. To perform approximate inference of z, we first define ℓ i (z) as the log likelihood of D i given z, with w i marginalised out. That is,EQUATIONThen, given p 0 (z) and {ℓ i (z)} i , we define p * (z) as (1), and proceed to find an approximation p(z) ≈ p * (z) using the methods described in earlier sections. Note that sampling from the tilted distribution p i (z) requires jointly sampling over z and w i and taking the marginal. EP is a particularly appealing framework for performing inference in this setting, as the dimensionality of the sampled distributions is constant with respect to m, mitigating the curse of dimensionality experienced by conventional MCMC approaches [48] . In our experiments we used NUTS [27] to perform the sampling, consistent with prior work [48, 52] . In each experiment we compared EP-η and EP-µ with EP, in the manner used by Xu et al. [52] and Vehtari et al. [48] . When F was MVN we also compared to SNEP [17] . For the models with normal-inverse Wishart (NIW) F, we were unable to find an initialisation for SNEP that would allow us to perform a meaningful comparison. We discuss this point further in Appendix J.To evaluate the performance of the different variants, we monitored KL divergence to an estimate of the optimum, obtained by running EP with a large number of samples. We used 500 different hyperparameter settings for each variant, chosen using random search, and repeated each run 5 times using different random seeds for NUTS. All runs of EP-η and EP-µ were performed with n samp = 1 and n inner = 1. The results of our evaluation are summarised by the Pareto frontiers in Figure 2 .We see that EP-η and EP-µ can typically reach a given level of accuracy (distance from an estimate of the optimum) faster than EP, often significantly so. We stress that the frontiers for EP-η and EP-µ are traced out with n samp and n inner each fixed to 1, with no sample thinning, varying only the step size hyperparameter ϵ, with higher accuracy/cost regions of the frontier corresponding to lower values of ϵ. In contrast, to trace out the frontier of EP we must jointly vary α, n samp , and the thinning ratio. See Appendix K for examples of this effect. The performances of EP-η and EP-µ are very similar. We found SNEP to be significantly slower than the other methods, consistent with findings by Vehtari et al. [48] . Pareto frontiers with respect to wall-clock time can be found in Appendix L. Code for these experiments can be found at https://github.com/cambridge-mlg/fearless-ep. We now provide a brief overview of individual experiments, with further details given in Appendix J.Hierarchical logistic regression with MVN prior Hierarchical logistic regression (HLR) is used to perform binary classification in a number of groups when the extent to which data should be pooled across groups is unknown [15] . We performed approximate Bayesian inference in a HLR model using a MVN prior over the group-level parameters. We applied this model to synthetic data, generated using the same procedure as Vehtari et al. [48] , which was designed to be challenging for EP. In this experiment z ∈ R 8 , w i ∈ R 4 ∀ i, and each of the m = 16 groups had n = 20 observations. In the upper-left panel of Figure 2 we see that EP-η and EP-µ typically reach a given level of accuracy faster than EP and SNEP.We also performed approximate Bayesian inference in a similar HLR model, but using a normal-inverse Wishart (NIW) prior over group-level parameters, allowing for correlation of regression coefficients within groups. We applied this model to political survey data, where each of m = 50 groups corresponded to responses from a particular US state, with 7 predictor variables corresponding to characteristics of a given survey respondent, so that w i ∈ R 7 ∀ i. There were n = 97 survey respondents per state. The data, taken from the 2018 Cooperative Congressional Election Study, related to support for allowing employers to decline coverage of abortions in insurance plans [33, 43] . The upper-right panel of Figure 2 shows that EP-η and EP-µ consistently reach a given level of accuracy significantly faster than EP.Cosmic radiation model A hierarchical Bayesian model was used by Vehtari et al. [48] to capture the nonlinear relationship between diffuse galactic far ultraviolet radiation and 100-µm infrared emission in various sectors of the observable universe, using data obtained from the Galaxy Evolution Explorer telescope. In this model each w i ∈ R 9 parameterised a nonlinear regression model using data obtained from one of m sections of the observable universe. The m regression problems were related through hyperparameters z ∈ R 18 , which parameterised the section-level parameter densities. The prior, p 0 ∈ F, was MVN. The specifics of the nonlinear regression model are quite involved and we refer the reader to Vehtari et al. [48] for details. We were unable to obtain the dataset, and so we generated synthetic data using parameters that were tuned by hand to try and match the qualitative properties of the original data -see Appendix N for examples. We used a reduced number of m = 36 sites and n = 200 observations per site to allow us to perform a comprehensive hyperparameter search. The lower-left panel of Figure 2 shows once again that EP-η and EP-µ consistently reach a given level of accuracy significantly faster than EP and SNEP.Neural response model A common task in neuroscience is to model the firing rates of neurons under various conditions. We performed inference in a hierarchical Bayesian neural response model, using recordings of V1 complex cells in an anesthesised adult cat [4] . In this dataset, 10 neurons in a specific area of cat V1 were simultaneously recorded under the presentation of 18 different visual stimuli, each repeated 8 times, for a total of 144 trials. We modelled the observed spike counts of the 10 neurons in each trial as Poisson, with latent (log) intensities being MVN, with mean and covariance (collectively z ∈ R 65 ) drawn from a NIW prior p 0 ∈ F. Inference of z amounted to inferring the means, and inter-neuron covariances, of latent log firing rates across stimuli. We grouped the data into m = 8 batches of n = 18 trials each, so that the latent log firing rates for all trials in batch i were jointly captured by w i ∈ R 180 . The lower-right panel of Figure 2 shows that EP-η and EP-µ either match or beat the speed of EP for reaching any given level of accuracy. For high levels of accuracy the performances of the methods appear to converge with one another, but note that the asymptote is not exactly zero because the optimal solution was estimated by running EP with a large (but finite) number of samples. 
limitations: 
 In this work we have largely assumed that the drawing of MC samples is the dominant cost. In practice, other costs become relevant, which may shift the balance in favour of using more samples to estimate updates. We discuss this in Appendix I, along with strategies for reducing computational overheads. We note, however, that wall-clock time results for our evaluation experiments (in Appendix L) are in line with the results of Section 4, indicating that the drawing of MC samples was indeed the dominant cost.When EP and its variants-including those introduced here-are combined with MC methods, the underlying samplers typically have hyperparameters of their own, which are often adapted using so-called warmup phases. While EP-η and EP-µ significantly reduce the complexity of tuning hyperparameters specific to EP, they do not help with those of the underlying samplers. This means that the complexity of hyperparameter tuning (for all EP variants) is still greater than that for direct MC methods. We used comprehensive hyperparameter searches in our experiments in order to perform a meaningful comparison with baseline methods, necessarily limiting the scale of problems we could tackle.Our focus has been on developing improved methods for performing EP when the updates must be estimated with noise, with the motivation for this kind of setting having been discussed by prior work [52, 17, 48] . EP and its variants are agnostic to the choice of sampling kernel, and as such we have not discussed issues relating to the performance of the underlying samplers. In the left panel of Figure 3 however, we highlight the relative inefficiency of using EP-η with NUTS as the underlying sampler when compared to conjugate-computation variational inference (CVI), an efficient variational inference method [31] . While EP-η is able to obtain much more faithful posterior approximations than CVI, it is significantly slower when measured in wall-clock time, despite the gains made over its predecessors. However, this is due to the relative cost of drawing independent samples with NUTS, as the sample efficiency of EP-η is roughly equivalent to that of CVI when an (illustrative) "oracle" sampling kernel is used, as demonstrated in the left-middle panel of Figure 3 . Further discussion and details of the experiments behind Figure 3 can be found in Appendix M. 
future work: 
 In Section 5, we demonstrated that the efficiency of EP variants is constrained by that of of the underlying samplers. There are several ways in which their performance could be improved by using knowledge of p to guide sampling from p i . One such improvement would be to use p to set MCMC hyperparameters directly, obviating the need to adapt them during warm-up phases, e.g. when p ∈ F is MVN, its precision matrix could be used directly as the mass matrix in Hamiltonian Monte Carlo schemes. Another potential improvement is that samples from p could be used to re-initialise MCMC chains at each iteration, providing approximate independence between updates for a small additional cost. It would also be worthwhile considering how our methods can be efficiently implemented on modern compute hardware to better take advantage of the inherently parallel nature of EP, e.g. practical performance could be improved by using sampling kernels that are themselves designed for efficient parallel execution [26] . Performance may be further improved in distributed settings by using asynchronous extensions of EP in order to minimise communication-related overheads and delays [17, 52] . We leave investigation of these ideas to future work.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: The experimental settings are summarised in Section 4, with full details, including hyperparameter search spaces, given in Appendix J. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 In this work, we used a novel interpretation of the moment-matching updates of EP to motivate two new EP variants that are far more robust to MC noise. We demonstrated that these variants can offer an improved speed-accuracy trade-off compared to their predecessors, and are easier to tune.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: The main results are summarised by the Pareto frontiers in Figures 2 and 7 . The figures display error bars, corresponding to the minimum and maximum value observed for chosen hyperparameter settings over 5 random seeds, which were used to seed the MCMC samplers, as we state in Section 4.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: Appendix J gives details of the hardware used to run the experiments, as well as the number of runs carried out. This, combined with the wall-clock time results in L, give an estimate of the compute resources needed to recreate our experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 • The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: We have made our experiment code public, but do not consider it one of the paper's main contributions.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: No crowdsourced experiments or research involving human subjects were conducted as part of this work.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? 
 
 Appendix: In this appendix, we give a very brief overview of exponential families of distributions. This is largely a condensed summary of relevant material from Wainwright and Jordan [49] , and we refer the reader to the original work for a far more detailed treatment.The exponential family of distributions F, defined by the d-dimensional, vector-valued sufficient statistic function s(.), and base measure ν(.), has density functionEQUATIONtaken with respect to ν(.), where A(η) = log exp(η ⊤ s(z))dν(z) is known as the log-partition function. η ∈ Ω are known as the natural parameters, and are used to index a specific member of F. Ω is the set of normalisable natural parameters of F, given byEQUATIONand is known as the natural domain of F. Ω is a convex set, and when it is open, the family F is said to be regular -we shall only consider regular families in this work. When the components of s(.) are linearly independent, F is said to be minimal. For minimal families, each distribution in the family is associated with a unique natural parameter vector η. An exponential family that is not minimal is said to be overcomplete, in which case, each distribution is associated with an entire affine subspace of Ω.The expected sufficient statistics of a distribution with density p(.) with respect to base measure ν(.), are given by E p(z) [s(z)]. Let M be defined as the set of expected sufficient statistics that can be attained by any density p(.) with respect to base measure ν(.) -that is,EQUATIONAll vectors in the interior of M, M • , are realisable by a member of F, and so µ ∈ M • provides an alternative parameterisation of F, in which µ are known as the mean parameters, and M • is known as the mean domain,.For minimal families, A(.) is strictly convex on Ω. The convex dual of A(.) is defined asEQUATIONand on M • , A * (µ) is equal to the negative entropy of the member of F with mean parameter µ.The mean parameters of the member of F with natural parameter η can be obtained by µ = ∇A(η).For minimal families, this mapping is one-to-one, and the reverse map is given by η = ∇A * (µ). For this reason, ∇A(.) and ∇A * (.) are sometimes referred to as the forward mapping and backward mapping respectively. In this work, we say that a family F is tractable if both the forward and backward mappings can be evaluated efficiently.The Fisher information matrix (FIM) of an exponential family distribution, with respect to natural parameters η, is given by, F (η) = ∇ 2 A(η). Furthermore, in minimal families, the FIM with respect to mean parameters µ is given by F (µ) = ∇ 2 A * (µ). Using the forward and backward mappings, we also haveEQUATIONEQUATIONand from the inverse function theorem,EQUATIONEQUATIONIn this appendix we give a more conventional presentation of the EP updates, and show equivalence with the atypical presentation of Section 2. We will assume here that β i = 1 ∀ i, but similar reasoning can be used to show equivalence in the power EP case, where β i ̸ = 1. Note that the update of power EP is often presented as having scale proportional to β i , with a separate damping parameter. We subsume both into α, with α = 1 corresponding to the default damping suggested by Minka [34] .Given a target distribution with the formEQUATIONEP attempts to find an approximation p ∈ F such thatEQUATION)The i-th site parameter λ i can loosely be interpreted as the natural parameters of a F-approximation to the i-th target factor, exp(ℓ i (z)). EP is also often used to approximate the normalising constant of p * (z) -we have omitted details here, but see e.g. Bishop [3] .In order to update λ i , EP performs a KL-divergence minimisation between p(z) and a local approximation to the target distribution, consisting of just one "real" factor and (m -1) approximate ones. Specifically,EQUATIONwhere pi (z) ∝ exp ℓ i (z)EQUATIONand we have used the notation p(z; λ i ) to make the dependence on λ i explicit. It can be easily shown that the solution to ( 26) is found by moment matching. That is, the KL-divergence is minimised whenEQUATIONFor a minimal exponential family F, there is a unique member of F with mean parameters E pi(z) [s(z)], and its natural parameters are given by ∇A * (E pi(z) [s(z)]). By equating this with the parameters of p(z), we haveEQUATIONUpdate (29) is the "standard" update of EP, and can be performed either serially or in parallel (over i). Often a damping factor α is used, givingEQUATIONIn Section 2 we presented EP and several variants as performing some number of "inner" updates, to decrease L with respect to λ i , with an "outer" update with respect to θ. These updates are given by Inner update:EQUATIONOuter update:EQUATIONImmediately after an outer update we have that θ = η 0 + j λ j , and so for β i = 1 we haveEQUATIONWe stated in Section 2 that (non-double-loop) EP can be viewed as performing a single inner update per outer update, and so by rolling (31) and (32) into a single update for λ i , we haveEQUATIONwhich is identical to (30) .Begin by taking gradients of (3), with respect to λ i , givingEQUATIONAt a minimum of L with respect to λ i , both sides of ( 35) must be equal to zero. We can use this to define a fixed-point condition with respect to λ i , and its updated valueEQUATIONwhere p i ∈ F i denotes the i-th tilted distribution, defined by (6) . Applying ∇A * (.), the inverse of ∇A(.), to both sides of (36) , and rearranging terms, we recover the moment-matching update of EPEQUATION)Often a damping parameter is used to aid convergence, giving the more general updateEQUATIONwhere the level of damping is given by (1 -α). It was shown by Heskes and Zoeter [25] that (4) follows a decrease direction in L with respect to λ i , and so is guaranteed to decrease L when α is small enough.Let pi (t) be the member of F with natural parameter η i (t) (λ i ) = η 0 + j λ j (t) + α -1 (λ i -λ i (t) ). Proposition 1 states that that the moment-matching updates of EP can be viewed as performing NGD in L with respect to the mean parameters, µ i , of pi (t) . We restate the proposition below, and give a proof. Note that the right hand side of (7) maps the EP update (4) onto mean parameters of pi (t) . The left hand side of ( 7) is simply a NGD update in µ i . Proposition 1. For α > 0, the moment-matching update of EP (4) is equivalent to performing an NGD step in L with respect to the mean parameters of pi (t) with step size α -1 . That is, forµ i = E pi (t) (z) [s(z)],and Fi (t) (µ i ) the FIM of pi (t) with respect to µ i , we haveEQUATIONProof. Taking the left hand side of (7), we haveEQUATIONwhere the penultimate equality follows from the definitionµ i = E pi (t) (z) [s(z)],and the fact that λ j (t) = λ j ∀ j before the update. All that remains is to show that the right hand side of ( 7) is equalEQUATIONUsing the relation η i = η i (t) (λ i ), and thereforeEQUATIONFor α = 1, we haveEQUATIONLet the bias in dimension d of µ i (t+1) , after an update at time t, be defined as E[µ i (t+1) -μi (t+1) ] d , where μi (t+1) is the value of µ i (t+1) after a noise-free update, and expectation is taken over the sampling distributions of all parallel updates at time t. Propositions 2 and 3, which we restate and prove below, summarise the effect of decreasing α and ϵ on EP and EP-µ, respectively. Proposition 2. After update (4) is executed in parallel over i, as α → 0 + , both the expected decrease in L, and the bias E[µ i (t+1) -μi (t+1) ] d , are O(α) for all d.Proof. Let η i (t) and µ i (t) be the pre-update natural and mean parameters of pi (t) respectively. Furthermore, let µEQUATIONbe the mean parameters of pi (t) after an EP inner update (4) but before being mapped to mean parameters of pi (t+1) through map (9) . ξ i is some zero-mean noise, e.g. due to MC variation. The second equality follows from the results of Appendix D. Similarly, let η ′ i(t) = ∇A * (µ ′ i (t)) be the corresponding natural parameters. Now apply map (9) to convert to mean parameters of pi (t+1) ,EQUATIONwhere the last equality follows from the definition of η i (t) (.), and by observing thatEQUATIONfor all j. Let μ′ i (t) be the mean parameters of pi (t) after a noise-free update, but before mapping to mean parameters of pi (t+1) , so that µ ′ i (t) = μ′ i (t) + ξ i . Substituting this in, we haveEQUATION)Taking a Taylor expansion around α = 0, we haveEQUATIONNow subtract the noise-free update μi (t+1) , and take expectations with respect to ξ j for j = 1, . . . , m,EQUATION)the first term on the right hand side does not have zero expectation in general-∇A * (.) is not generally affine-and so the bias is O(α) as α → 0 + .To see that the expected change in L is also O(α), take the Taylor expansion of L along the update direction-obtained by subtracting λ i from the right hand side of ( 4))-as a function of the step size around α = 0,EQUATIONwhich is clearly O(α) as α → 0 + . Proposition 3. After update (12) is executed in parallel over i, as ϵ → 0 + , the expected decrease in L is O(ϵ), and the biasE µ i (t+1) -μi (t+1) d is O(ϵ 2 ), for all d.Proof. Let η i (t) and µ i (t) be the pre-update natural and mean parameters of pi (t) respectively. Furthermore, let µEQUATIONbe the mean parameters of pi (t) after an EP-µ update (12) but before being mapped to mean parameters of pi (t) through (9) . ξ i is some zero-mean noise, e.g. due to MC variation. The second equality follows from the definition µ i (t) = E pi (t) (z) [s(z)], and the fact that pi (t) = p before an update. Now apply map (9) to convert to mean parameters of pi (t+1) ,EQUATIONwhere the last equality is obtained from the definition of η i (t) (.) for α = 1, and usingEQUATIONfor all j. Substituting in the definition of µ ′ j (t) , we haveEQUATION)Taking a Taylor expansion around ϵ = 0, we haveEQUATIONFinally, subtract the noise free μi (t+1) , and take expectations (over ξ j for j = 1, . . . , m) to obtain the biasEQUATIONBy assumption E[ξ j ] = 0 ∀ j, hence the first order term disappears, and the bias is O(ϵ 2 ).To see that the expected decrease in L is also O(ϵ), note that we are simply following the gradient in L with respect to µ, multiplied by the inverse Fisher and scaled by ϵ. Neither the reparameterisation (from λ i to µ i ) or the Fisher depend on ϵ, and so the change in L must be O(ϵ) as ϵ → 0 + by simple Taylor expansion arguments.In this paper we have shown that the updates of EP, EP-η, and EP-µ can be interpreted as performing NGD of L with respect to the distributions {p i (t) } i . The SNEP method of Hasenclever et al. [17] , in contrast, performs NGD of L with respect to the site potentials, treating them as bona fide distributions in F for the purposes of NGD -that is, NGD is performed with respect to the distributions with densities given by exp(λ ⊤ i s(z) -A(λ i )) for i = 1, . . . , m. In Section 3.4, we argued that the statistical manifolds of {p i (t) } i are far more relevant for the optimisation of L than those of the site potential pseudo-distributions. In this appendix we shall justify this claim.The loss function L, given by (3), has a unique optimum with respect to {λ i } i when ∇ λi L(θ, λ 1 , . . . , λ m ) = 0 ∀ i, which implies thatEQUATIONThe minimisation of L with respect to {λ i } i can therefore be viewed as a joint optimisation in the m + 1 distributions p, p 1 , . . . , p m , with the minimum found when their expected F-statistics match.Intuitively, it seems that optimisation of L with respect to λ i should be performed with consideration for the geometry of p and p i . By instead performing NGD with respect to the distributions {p i (t) } i , in the case of EP, EP-η and EP-µ, or with respect to the site potential psuedo-distributions, in the case of SNEP, these methods can be interpreted as using surrogate distributions for NGD [46] . It is possible to show that the Fisher information matrix (FIM) of pi (t) with respect to λ i is equal to that of p, up to a scalar constant, and furthermore, it can also be seen as a close approximation to that of p i . This motivates the use of NGD with respect to pi (t) for the optimisation of L. In contrast, the FIMs of the site-potential pseudo-distributions, used by SNEP, can bear little relation to those of p or p i . In particular, the difference will be greatest when the full approximation η 0 + j λ j is significantly different than λ i , and we believe this may be why Hasenclever et al. [17] found that increasing the number of sites reduced the convergence speed of SNEP.We can make these arguments more concrete by considering the curvature of L. Taking second derivatives of L with respect to λ i , we haveEQUATIONLet λEQUATIONwhere η i and µ i are the natural and mean parameters of pi (t) , respectively. The Hessian of the second term will not generally be available in closed form, or easily invertible. We can however find a convenient approximation. Note that after an outer update, at time t, we have θ = η 0 + j λ j (t) . Also we have that c exp λ ⊤ i s(z) ≈ exp(ℓ i (z)), at least near convergence, where c is some unknown scalar constant. Combining these, we haveEQUATION)Combining this with (58), we haveEQUATIONand if we use this curvature approximation to perform a quasi-Newton step in L with respect to λ i and step size ε, we haveEQUATIONBy letting ϵ = ε(1 + β -1 i ) -1 we have equivalence with the update of EP-η, given by (11) . We therefore have that EP-η is equivalent to performing a quasi-Newton step in L. Furthermore, the approximate curvature becomes exact at the optimum if ℓ i (z) is an affine function of s(z), e.g. if exp(ℓ i (z)) is an (unnormalised) member of F.In contrast, let us consider the implicit curvature approximation used by SNEP. Let γ i be mean parameters of the member of F with natural parameter λ i . To perform the inner minimisation of the variational objective, SNEP performs NGD in L with respect to γ i . The implicit curvature matrix used by the SNEP update is then the Fisher with respect to γ i , which is given by ∂λi ∂γi = ∇A * (γ i ). Now consider the actual curvature of L with respect to γ i ,EQUATIONwhere we have used λ i (γ i ) = ∇A * (γ i ) to denote the backwards map. Near the optimum we have ∂L ∂λi ≈ 0, and so the curvature will be approximately equal to the first term as we approach the optimum. We have already shown thatEQUATIONand so the curvature with respect to γ i isEQUATION)Compare this with the implicit curvature used by SNEP of ∂λi ∂γi , and it is clear that the two are proportional only when η 0 + j λ j = λ i . This suggests that if there are many sites (m is large), or if the prior is very informative, the curvature matrix implicity used by SNEP is likely to be significantly different from the true curvature, even near the optimum.Let c samp be the cost of drawing a single sample from one of the tilted distributions. Let c fwd be the cost of converting from natural to mean parameters in the base family F (the forward mapping). Similarly, let c bwd be the cost of converting from mean to natural parameters in F (the backward mapping). Recall that m is the number of sites. We now state the computational complexity of a round of parallel updates in each of the EP variants evaluated in this paper.EP EP draws n samp samples from each of the m sites. It then performs m backward mappings, to map from moments back to updated site parameters. The total cost is then mn samp c samp + mc bwd .EP-η EP-η draws n samp samples from each of the m sites. Each update also involves one forward mapping per iteration. In addition, each iteration also requires m JVPs through the backward mapping ∇A * (.). The primals of this JVP are the same for each site, so the linearisation only needs to be performed once per update, costing ≈ 2c bwd , but we have m tangents. The overall cost of an EP-η iteration is therefore ≈ mn samp c samp + c fwd + (2 + m)c bwd .EP-µ EP-µ draws n samp samples from each of the m sites. Each update involves one forward mapping, and m backward mappings per iteration, but notably does not require any JVPs. The overall cost of an EP-µ iteration is therefore mn samp c samp + c fwd + mc bwd .SNEP SNEP draws n samp samples from each of the m sites, and performs m backward mappings per iteration. The total cost is then mn samp c samp + mc bwd , which is the same as EP.In this work we largely assume that c samp is the dominant cost. This is often the case when the sampling is performed using MCMC. For example, in the NUTS sampler with default numpyro settings -as used in our evaluation -each sample can involve evaluating up to 1024 gradients of the tilted distribution log density.The other costs, c fwd and c bwd , depend on F. In the case of a MVN family with dense covariance, the foward and backward mappings are O(d 3 ), for d the dimensionality of z. When d is large, c fwd and c bwd become more significant, in which case the balance will shift in favour of using more samples per update (and taking larger steps). When d is very large, however, it may be necessary to choose a diagonal covariance base family instead, in which case the mappings can be performed in O(d) time. Also note that when d increases, c samp will also typically increase at a rate faster than d, even with optimal tuning [8] . In our evaluation we used just a single sample per update for EP-η and EP-µ, and they significantly outperformed the baselines when measure in both NUTS steps and wall-clock-time.We also note that the per-iteration cost of c fwd and c bwd could be reduced by exploiting the results of previous iterations. For example, in the case of a MVN with dense covariance, the cost of the forward and backward mappings are dominated by matrix inversions. This cost could be significantly reduced by using the result from the previous iteration to warm-start a Newton-style iterative inversion routine. This approach was used by Anil et al. [2] to reduce the cost of computing inverse matrix roots as part of a wider deep-learning optimisation scheme. We did not make use of such optimisations in our experiments.In this appendix, we give details about the experiments presented in Section 4. To evaluate the performance of the different variants, we monitored KL divergence to an estimate of the optimum, obtained by running EP to convergence with a large number of samples (n samp = 10 5 for the final iterations). We used 500 different hyperparameter settings for each variant, chosen using random search, and repeated each run using 5 different random seeds, which were used to seed the MCMC samplers. Hyperparameter settings that failed in any of the 5 runs were discarded. We used this setup for all experiments in Section 4.The x-axis in Figure 2 of Section 4 corresponds to the number of NUTS steps, or more specifically, the number of NUTS candidate evaluations. This number is hardware and implementation agnostic, and roughly corresponds to the total computational cost incurred up to that point. Wall-clock time results are given in Appendix L.Hyperparameter search spaces The step size for EP, α, was drawn log-uniformly in the range (10 -4 , 1), and n samp was drawn log-uniformly between [d + 2.5, 10000.5) and then rounded to the nearest integer, where d is the dimensionality of z. Note that the debiasing estimator used by Xu et al. [52] for MVN families is not defined for n samp ≤ d + 2. The thinning ratio was drawn uniformly from {1, 2, 3, 4}. The step size for EP-η and EP-µ, ϵ, was drawn log-uniformly in the range (10 -5 , 10 -2 ), and n samp was fixed to 1. The step size for SNEP was drawn log-uniformly in the range (10 -5 , 10 -2 ), with n samp and n inner (independently) drawn log-uniformly in the range [.5, 10.5) and then rounded to the nearest integer.All the methods considered can be used in a double-loop manner by taking n inner > 1 inner steps per outer update. We did not find this necessary in any of our experiments, and so we fixed n inner = 1 for EP, EP-η and EP-µ. For SNEP we included n inner in the search space (as described above) to ensure that setting it to 1 was not negatively impacting its performance.Hardware All experiments were executed on 76-core Dell PowerEdge C6520 servers, with 256GiB RAM, and dual Intel Xeon Platinum 8368Q (Ice Lake) 2.60GHz processors. Each individual run was assigned to a single core.Software Implementations were written in JAX [7] , with NUTS [27] used as the underlying sampler. We used the numpyro [42] implementation of NUTS with default settings. For experiments with a NIW base family, we performed mean-to-natural parameter conversions using the method of So [45] , with JAXopt [5] used to perform implicit differentiation through the iterative solve.We performed regular warm-up phases in order to adapt the samplers to the constantly evolving tilted distributions, consistent with prior work [52, 48] . The frequency and duration of these warm-up phases was also included in the hyperparameter search as follows. We drew the duration (number of samples) of each warm-up phase log-uniformly in the range [99.5, 1000.5) and then rounded to the nearest integer. We then drew the sampling-to-warm-up ratio log-uniformly in the range (1, 4) . This ratio then determined how frequently the warm-up was performed, which we rounded to the nearest positive integer number of updates. Note that the Pareto frontier plots in Figures 2 and 7 include the computation / time spent during warm-up phases.Site parameter initialisation For EP, EP-η and EP-µ, we initialised the site parameters to λ i = 0 ∀ i for all models, with the exception of the cosmic radiation model for reasons we discuss later. SNEP is not compatible with improper site potentials, and so for models where F was MVN we used λ i = (2m) -1 η 0 ∀ i, consistent with Vehtari et al. [48] . Unfortunately this initialisation strategy for SNEP is not always valid when F is NIW. This is because the site potentials in SNEP are restricted to being proper distributions, and this constraint was violated when using this initialisation in our experiments. We tried various other initialisation strategies, but none produced satisfactory results or allowed us to make a meaningful comparison, hence we omitted SNEP from the experiments with NIW F; namely, the political survey hierarchical logistic regression and neural response model experiments. We set β i = 1 ∀ i in all experiments, corresponding to the original (not power) EP of Minka [37] . All of the models in these experiments followed the same general structure, given byD w z i = 1, . . . , mEQUATIONwhere p 0 is a member of a tractable, minimal exponential family F. This is shown graphically in Figure 4 . We now provide details about individual experiments.In the hierarchical logistic regression (HLR) experiments, there were m groups (logistic regression problems) each with n covariate/response pairs, so that D i = {(x i,j , y i,j )}, x i,j ∈ R d , and y i,j ∈ {0, 1}, for i = 1 to m, and j = 1 to n. Each group had its own unobserved vector of regression coefficients w i ∈ R d , with density p(w i | z) parameterised by global parameters z.In the synthetic experiment, we had m = 16 groups, with d = 4 and n = 20. F was an MVN family, and z = (µ 1 , log σ 2 1 , ..., µ 4 , log σ 2 4 ) ∈ R8 corresponded to the means and log-variances for the independent normal density p(w i | z) = j N (w i,j | µ j , σ 2 j ). 8 The prior on z had density p 0 (z) = N (0, diag((4, 2, 4, 2, . . .))), where diag(.) constructs a diagonal matrix from its vector-valued argument. The dataset was generated using the same procedure as Vehtari et al. [48] .Political survey data with NIW prior In the political survey data experiment, F was a normalinverse-Wishart (NIW) family of distributions, and z = (µ, vech(Σ)) was used to parameterise the group-level coefficient density p(w i | z) = N (w i | µ, Σ). p 0 was a NIW prior on µ and Σ, with parameters µ 0 = 0, ν = 9, λ = 1, and Ψ = I (following Wikipedia notation). We used a log-Cholesky parameterisation of Σ for the purposes of sampling. The dataset consisted of binary responses to the statement "Allow employers to decline coverage of abortions in insurance plans (Support / Oppose)". We constructed m = 50 regression problems, corresponding to the 50 US states, and truncated the data so that there were exactly n = 97 responses for each state. We used 6 predictor variables from the dataset, corresponding to characteristics of a given survey participant. The predictors were binary variables conveying: age (3 groups), ethicity (2 groups), education (3 groups) and gender. We also included a state-level intercept to capture variation in the base response level between states, so that d = 7. This setup is based on one used by Lopez-Martin et al. [33] . The data is available at https://github.com/JuanLopezMartin/MRPCaseStudy.In this experiment, there were m nonlinear regression problems, each corresponding to a model of the relationship between diffuse galactic far ultraviolet radiation (FUV) and 100-µm infrared (i100) emission in a particular sector of the observable universe.The nonlinear regression model for sector i had parameters w i ∈ R 9 . The m regression problems were related through common hyperparameters z ∈ R 18 , which parameterised the section-level densities p(D i , w i | z) where D i is the observed data for sector i. F was the family of MVN distributions, and the prior on z ∈ R 18 had density p 0 (z) = N (0, 10I). The specifics of p(D i , w i | z) are quite involved and we refer the reader to Vehtari et al. [48] for further details.Vehtari et al. [48] applied this model to data obtained from the Galaxy Evolution Explorer telescope. We were unable to obtain the dataset, and so we generated synthetic data using hyperparameters that were tuned by hand to try and match the qualitative properties of the original data set (see Appendix N for examples). We used a reduced number of m = 36 sites and n = 200 observations per site to reduce computation, allowing us to perform a comprehensive hyperparameter search.Using the conventional site parameter initialisation of λ i = 0 ∀ i resulted in most runs of EP failing during early iterations. We found that this was resolved by initialising with the method used by Vehtari et al. [48] for SNEP, that is, λ i = (2m) -1 η 0 ∀ i, and so we used this initialisation for all methods. We note, however, that the performances of EP-η and EP-µ were largely unaffected by this change.In this experiment we performed inference in a hierarchical Bayesian neural response model, using recordings of V1 complex cells in an anaesthetised adult cat. 10 neurons in a specific area of cat V1 were simultaneously recorded under the presentation of 18 different visual stimuli, each repeated 8 times, for a total of 144 trials. Neural data were recorded by Tim Blanche in the laboratory of Nicholas Swindale, University of British Columbia, and downloaded from the NSF-funded CRCNS Data Sharing website https://crcns.org/data-sets/vc/pvc-3 [4] . z = (µ, vech(Σ)) ∈ R 65 was used to parameterise N (log r j ; µ, Σ), the density for latent log firing rates in trial j, for j = 1, . . . , 144. p 0 ∈ F was NIW with parameters µ 0 = 1, ν = 12, λ = 2.5, and Ψ = 1.25I (again following Wikipedia notation). The observed spike count for neuron k in trial j, x j,k ∈ N, was modelled as Poisson(c j,k , exp(r j,k )), for j = 1, . . . , 144, k = 1, . . . , 10. We again used a log-Cholesky parameterisation of Σ for sampling.We grouped trials together into m = 8 batches of n = 18 trials, so that w i = (log r 180(i-1)+1 , . . . , log r 180i ) ∈ R 180 was the concatenation of log firing rates for all 10 neurons across 18 trials, with D i = {(x j,1 , . . . , x j,10 )} the observed spike counts, for j = 180(i -1) + 1 to 180i.In this appendix, we examine the effect of varying hyperparameters of EP and EP-η on the synthetic hierarchical logistic regression experiment of Section 4. The results for EP-µ are similar to those of EP-η and so we do not consider it separately here. EP In the left panel of Figure 5 we show the effect of varying the number of samples used for estimating updates in EP -more accurate regions of the frontier generally require more samples. The middle panel, similarly, shows the effect of varying the step size, and the right panel shows the effect of varying the thinning ratio τ . Together these plots illustrate the difficulty of tuning EP in stochastic settings -tracing out the frontier is a three-dimensional problem. For example, to achieve the best accuracy within a compute budget of 10 7 steps, the practitioner would need to set 400 < n samp ≤ 500, 0.1 < α ≤ 0.3, and τ = 2, and any deviation from this would seemingly result in suboptimal accuracy. EP-η In the left panel of Figure 6 we show the effect of varying the number of samples used for estimating updates in EP-η -the frontier is traced out with n samp = 1. We note however that the difference between n samp = 1 and n samp = 10 is relatively small, and so it may be sensible to choose n samp > 1 to make efficient use of parallel hardware, or to minimise per-iteration overheads (see Appendix I). The middle panel, similarly, shows the effect of varying the step size, with more accurate regions of the frontier corresponding to a smaller step size. This also suggests that in practice it may make sense to set ϵ relatively large at first and then gradually decay it to improve the accuracy. Finally, the right panel shows the effect of varying n inner , the number of inner steps performed per outer update. n inner > 1 corresponds to using EP-η in "double-loop" mode. We did not find it necessary to increase n inner above one to obtain convergence in our experiments, but Figure 6 (right) demonstrates that doing so would have a relatively small impact on performance. Together these plots illustrate that hyperparameter tuning for EP-η is relatively straightforward. The frontier can be largely be traced out by varying ϵ, and is relatively insensitive to n samp and n inner . In Figure 7 , we show Pareto frontiers, with the y-axis showing the lowest KL divergence achieved by any hyperparameter setting, and the x-axis shows the cumulative number of seconds elapsed -that is, the wall-clock time equivalent of Figure 2 . These results are in broad agreement with those of Figure 2 , which suggests that the sampling cost does indeed dominate the computational overheads of EP-η and EP-µ in these experiments. The wall-clock time results for EP-η and EP-µ would likely be improved by using more than one sample per update, by making more efficient use of hardware resources and minimising per-iteration overheads.We note that these times are necessarily implementation and hardware dependent. We did not make particular efforts to optimise for wall-clock time. In Appendix I we discuss approaches for minimising per-iteration overheads. These would likely improve wall-clock time performance for all methods, but should disproportionately favour EP-η and EP-µ, due to their frequent updates and larger per-iteration overheads compared to EP and SNEP. Our focus in this work has been on developing improved methods for performing EP when the updates must be estimated with noise. A major motivation for this setting is that it potentially enables EP to be used as a so-called black-box inference method, significantly expanding the set of models it can be applied to. EP has several computational advantages over direct MCMC approaches, the prevailing dominant class of black-box mehods, as discussed by earlier works [52, 17, 48] . There is another popular class of black-box inference methods however; namely, those of variational inference (VI). 9In order to understand the trade-offs involved in using EP and its variants over VI, we performed experiments to explore the relative strengths and weaknesses of EP-η and conjugate-computation variational inference (CVI) [31] , an efficient VI method.Experimental setup We used EP-η and CVI to perform approximate Bayesian inference using the same hierarchical logistic regression model (with MVN prior) and synthetic dataset as Section 4. For our evaluation metrics, we used both forward and reverse KL divergences between the current (MVN) approximation, and a MVN "target" distribution. The target distribution was a MVN distribution fitted using 500,000 samples from the posterior, obtained by running NUTS for 1 million samples and then discarding the first half as a warm-up. For CVI, we used a structured MVN approximation that had the same conditional independence structure as the true posterior -that is, it modelled all pairwise dependencies within z, within each w i , and between z and each w i , but did not model direct dependence between w i , w j for i ̸ = j. We used 500 random hyperparameter settings for each method, using the search spaces described below. Unless stated otherwise (below), all other details were as described in Appendix J.Hyperparameter search spaces For EP-η we used the same hyperparameter search space as described in Section J. For CVI, the step size was drawn log-uniformly in the range (10 -5 , 1), and the number of samples used to estimate the update was drawn log-uniformly in the range [.5, 100.5) and then rounded to the nearest integer. The site parameters for CVI were initialised to be the parameters of a zero-mean MVN distribution with scaled identity covariance matrix, with the scale drawn log-uniformly from (10 -5 , 10 5 ).We compared the wall-clock time performance of EP-η and CVI, with the former using NUTS as the underlying sampler. For both methods, we chose the hyperparameters that gave the lowest reverse divergence between the approximation and the MCMC target distribution after 100 seconds. The results of this experiment are shown in the left panel of Figure 3 . We can see that EP-η converges at least one order of magnitude slower than CVI on this problem, when measured in wall-clock time. This difference is largely due to the relative cost of drawing samples for the two methods. While CVI also uses samples to estimate its updates, it uses samples from the current MVN approximation, which are cheap to generate using standard methods. EP-η on the other hand (as with other EP variants), requires samples from the tilted distributions to estimate the updates, and in this experiment we used NUTS to draw those samples. Each NUTS sample can involve many evaluations of the energy function gradient, and furthermore, because MCMC samples are generally autocorrelated, the cost of drawing approximately independent samples can be even higher. We used NUTS to be consistent with prior work, and because our focus in this work has primarily been the relative performance of different EP variants, we did not make particular efforts to improve the performance of the underlying samplers, beyond tuning their hyperparameters in a fairly standard manner. However, all of the EP variants considered in this paper, including EP-η, are agnostic to the choice of sampling kernel used, and we believe there is scope to significantly improve the efficiency of the underlying samplers, for reasons we briefly discuss in Section 6.In order to decouple the performance of EP-η from that of the underlying sampler, we repeated the experiment using an "oracle" sampling kernel for EP-η. This oracle kernel was simply NUTS but with a thinning ratio of 100, so that each sample was an approximately independent sample from the tilted distribution. We then compared EP-η and CVI using the same metrics as the previous experiment, but with "time" measured in samples. We chose the hyperparameter setting for each method that achieved the lowest reverse KL divergence after 1,500 samples. When measured on this basis we see that the convergence speeds of CVI and EP-η are roughly equivalent, demonstrating that the sample efficiency of the two methods are similar. This should not be too surprising, as CVI is (mathematically, but not algorithmically) equivalent to the limiting case of power EP (as β i → ∞) [11, 51] .Approximation quality Figure 3 demonstrates that EP-η is able to achieve a more faithful approximation of the target distribution when measured by either the forward or reverse KL. At first it may seem surprising that EP is able to obtain a lower reverse KL divergence than CVI, but there are two reasons for this apparent anomaly. First, with EP, we have a MVN approximation for z only, with the approximate posterior over local variables (w i for i = 1, . . . , m) represented only implicitly with samples. CVI on the other hand necessarily optimises a joint MVN over all variables, and the z marginal of the optimal VI posterior over all variables is not the same as the optimal VI posterior over z. The second reason is that our target distribution is essentially a moment-matched MVN approximation of the true posterior, which would naturally tend to favour the forward KL divergence that is (approximately) targeted by EP. Nevertheless, a qualitative assessment of the pairwise marginals, as seen in the two rightmost panels of Figure 3 , shows that the true posterior is more faithfully approximated by EP-η than CVI, with the latter significantly underestimating uncertainty in the inter-group variability parameters, {log σ i } i . The full set of pairwise posterior marginals can be seen in Figure 8 . Note that the CVI approximation appears more constrained than might be expected based on the marginal plots alone. This is due to the zero-avoiding nature of the reverse KL divergence -non-Gaussianity in one variable can constrain the marginal distributions of others, particularly when there is strong dependence. The local variables, {w i }, can have both significant non-Gaussian posterior structure, and strong dependence on {log σ i } i .Vehtari et al. [48] used a hierarchical Bayesian model to capture the nonlinear relationship between diffuse galactic far ultraviolet radiation (FUV) and 100-µm infrared emission (i100) in various sectors of the observable universe, using data obtained from the Galaxy Evolution Explorer telescope. We were unable to obtain the dataset, and so we generated synthetic data using hyperparameters that were tuned by hand to try and match the qualitative properties of the original data set. Example data generated using these hyperparameters is shown in Figure 9 . For comparison, see Figure 9 of Vehtari et al. [48] .Figure 9 : Synthetic data, generated using the cosmic radiation model of Vehtari et al. [48] . Each plot shows galactic far ultraviolet radiation (FUV) versus infrared radiation (i100) for a single sector of the observable universe.Justification: We have read the NeurIPS Code of Ethics and believe our work confirms with it in every respect.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).We use e.g. p to refer to a distribution, and p(.) or p(z) for its density, throughout.Going forward, we assume i and j to be taken from the index set {1, . . . , m} unless otherwise specified.This interpretation is not precise, since it is not necessarily the case that λi ∈ Ω.Note the apparent contradiction that by decreasing α (increasing the damping) we actually increase the NGD step size. This is resolved by observing that the definition of p(t) i also changes with α.This can be computed using either forward or reverse mode automatic differentiation, as ∂η i ∂µ i is symmetric.EP too is performing NGD in mean parameters, but using a "step size" that also affects the distribution pi.See Appendix I for a detailed discussion of the computational costs of EP-η and EP-µ.Note that in this section we use µ to denote the mean of a normal or multivariate normal distribution. This should not to be confused with the usage for exponential family mean parameters in the main text.Expectation propagation is also a variational inference method, but we follow convention here by using variational inference to refer specifically to variational optimisation of the free energy / evidence lower bound. 


Decision: This paper presents a new view of the moment-matching updates of expectation propagation (EP) as  natural gradient descent (NGD) of a variational objective. This view allows the authors to propose two new methods that address the deficiencies of previous approaches, in particular that of not being robust to Monte Carlo noise when estimating the corresponding expectations. 
 
Overall, the reviews acknowledge the novelty of the approach (which has been further clarified in the rebuttal) and have praised the paper for the “clever” trick of reducing the bias of the MC expectation via change of parameterization and for being technically very strong. Indeed, one of the reviewers mentions “I'm indeed impressed that EP can work with single-sample estimation in the inner loop”. The main issues raised by the reviewers concern (1) the clarity of the presentation of the EP algorithm, (2) description of the results and (3) relation to previous work (in particular wrt to the work of Hasenclever et al, 2017). I believe the authors have addressed these issues in their rebuttal successfully with some of the reviewers updating their score accordingly. 

With this, I recommend acceptance.

Reviews: review: paper_summary: paper_summary: This paper considers new inference algorithms for EP. By framing the moment-matching equations of EP as a natural gradient update of a variational objective, they propose two new algorithms that are better suited for reducing/removing the bias introduced when sampling is required.
strength_weakness: strength_weakness: 1)	In general I find the paper quite hard to follow and of course this is not helped by the form of the standard EP equations. For example the main point of the paper is to handle the bias introduced by using sampling to estimate certain quantities within the EP equations. However the source of this sampling is only explicitly mentioned on line 132 in text. It would be much clearer if this was explicitly in Alg 1, and would also make it more convincing that Alg. 1 actually encapsulated many EP style algorithms.

2)	Following above some terms are not defined until later on in the paper. For example the convex conjugate in eqn 4, and on line 100 only defined one page later on line 123. 

3)	The paper is very descriptive but does not explain results. For example, the experiments only describe the set up, but all results and any conclusions are pushed into the appendix. This makes it hard to actually assess the contributions.
questions: questions: 1)	What is c_samp?

2)	In Eqn 24 should there be a $beta_i$ scaling term?

3)	In fig 1 a) why do most of the curves have a `u’ shape? Additionally why do the dashed orange/green curves not?

4)	Does sequential update (instead of in parallel) affect the bias?

5)	Why do you fix alpha to 1 and introduce a new step size eps? This looks like the only difference between Eqn 10 and the natural gradient part of 11. 

6)	The link between natural gradient steps  and EP updates has been established previously (in Heskes and also in Hasenclever). Why is the view taken in paper a novel perspective compared to these previous works?

7)	How does this work relate to `Bayes-Newton Methods for Approximate Bayesian Inference with PSD Guarantees’  Wilkinson et al, 2021 ? Which consider EP style algorithms with a natural gradient framework.
limitations: limitations: yes
paper_summary: paper_summary: This work addresses the sensitivity of expectation propagation (EP) to the randomness of Monte Carlo (MC) estimates involved in its update steps. It tackles this issue by recasting the moment matching step in EP as natural gradient descent (NGD) in the mean space of an exponential family distribution. The author identifies that the instability of EP to MC noise is due to the nonlinearity of the mirror map defined by the log partition function of the exponential family. By cleverly moving the NGD from the mean space to the natural parameter space, this issue is bypassed.
The author also studies the influence of the stepsize on the accumulation of MC error, finding that decreasing the stepsize of the NGD updates helps to reduce the bias.
strength_weakness: strength_weakness: The technical side of this work is very strong in my perspective, and I don't observe many weakness on this regard. 
However, I think the novelty of this NGD perspective is overclaimed. To my knowledge, the relationship between moment matching in the exponential family and NGD is very well known (e.g., [1][2][3]). Once it is identified that $\tidle p_i$ is in the exponential family, the derivation of NGD (in either mean space Prop1 or natural parameter space) becomes quite straightforward. I hope the authors can elaborate on this point, and include discussion on some relavant literature arolund prop 1.

I believe the real contribution of this work lies in the rigorous understanding of the influence of the stepsize on the accumulation of Monte Carlo (MC) error and the identification of the nonlinearity of the mirror map that leads to this error accumulation, which is already impactful to me.




- [1]: Conjugate-Computation Variational Inference : Converting Variational Inference in Non-Conjugate Models to Inferences in Conjugate Models, 2018
- [2]: MCMC-driven learning, 2024
- [3]: Distributed Bayesian Learning with Stochastic Natural Gradient Expectation Propagation and the Posterior Server
questions: questions: - How do you estimate the KL divergence in the experiments? Maybe I'm missing something, but the target posterior is only known up to a constant, so unbiased KL est is not feasible?

- I'm curious to see the relative performance of this work to NGVI [1] (specified in the Weakness section).
limitations: limitations: The limitation is well discussed.
paper_summary: paper_summary: The authors introduce two methods for Expectation Propagation [EP] (which itself can be understood as a method for approximate Bayesian inference) that are robust to the stochastic noise introduced by using approximate expectations in the inner loop of EP. To do so the authors re-interpret the moment-matching update as a particular Natural Gradient Descent (NGD) step. While stochasticity remains due to the expectation, by moving to natural parameters MCMC estimates of the expectation do not suffer from bias-introducing non-linearities. By significantly reducing bias (to zero in the case of EP-$\eta$) the authors show that stable and computationally efficient EP is possible with single-sample moment-matching updates. In experiments the authors show that the method is easier to tune than alternatives and exhibits promising empirical performance.
strength_weakness: strength_weakness: To my mind the major weakness of this submission is that the exposition is rather dense and at times hard to follow. This may be somewhat inherent in the technical nature of the topic, but I believe the authors could do a better job of guiding the reader. In particular a lot of space is devoted to somewhat extraneous modeling details ("...between diffuse galactic far ultraviolet radiation and 100-μm infrared emission in various sectors of the observable universe, using data obtained from the Galaxy Evolution Explorer..."), and this space in the main text could be re-purposed to offer more details on the algorithm, discuss the implications of the experimental results in more detail, or give more color on the "unified" EP presentation.
questions: questions: - math typos: line 804
- NIW is undefined where first introduced
- i guess in sec. 2.1 you should emphasize more that the domain of $\theta_i$ is $\Omega$, since $\theta_i \notin \Theta_i$ as one might expect? similarly i guess in practice given the form of the outer update $\theta_i = \theta_j$ for all $i, j > 0$? maybe this should be emphasized? can you please clarify?
- i understand that this is a paper about EP, but wouldn't it be valuable to benchmark against variational inference in at least one case? whether the authors like it or not, this remains a relevant comparison in my opinion.
limitations: limitations: To my mind the major limitation is the lack of comparison to alternative methods for approximate inference like variational inference.

Comments: comment: I thank the authors for the detailed reponse to address my questions. I'm happy to keep my score.
Thank you for your response. I think further clarifying and being more specific with the distinction with SNEP (Hasenclever et al. [13]) would be beneficial to the paper. For example on line 219 you state that SNEP can be viewed as doing natural gradient descent 'but with distributions that are more closely matched with those being optimised'  however it is not clear to me what 'distributions that are more closely matched with those being optimised'. 

In light of your response and the other reviews, I will increase my score by one point.
I thank the authors for their careful response. I maintain my score as is and hope the other reviewers will join me in recommending publication.