A Unified Perspective on Natural Gradient Variational Inference with Gaussian Mixture Models
Authors: Oleg Arenz
Philipp Dahlinger
Zihan Ye
Michael Volpp
Gerhard Neumann
Venue: Accepted by TMLR ()

Abstract: Variational inference with Gaussian mixture models (GMMs) enables learning of highly tractable yet multi-modal approximations of intractable target distributions with up to a few hundred dimensions. The two currently most effective methods for GMM-based variational inference, VIPS and iBayes-GMM, both employ independent natural gradient updates for the individual components and their weights. We show for the first time, that their derived updates are equivalent, although their practical implementations and theoretical guarantees differ. We identify several design choices that distinguish both approaches, namely with respect to sample selection, natural gradient estimation, stepsize adaptation, and whether trust regions are enforced or the number of components adapted. We argue that for both approaches, the quality of the learned approximations can heavily suffer from the respective design choices: By updating the individual components using samples from the mixture model, iBayes-GMM often fails to produce meaningful updates to low-weight components, and by using a zero-order method for estimating the natural gradient, VIPS scales badly to higher-dimensional problems. Furthermore, we show that information-geometric trust-regions (used by VIPS) are effective even when using first-order natural gradient estimates, and often outperform the improved Bayesian learning rule (iBLR) update used by iBayes-GMM. We systematically evaluate the effects of design choices and show that a hybrid approach significantly outperforms both prior works. Along with this work, we publish our highly modular and efficient implementation for natural gradient variational inference with Gaussian mixture models, which supports $432$ different combinations of design choices, facilitates the reproduction of all our experiments, and may prove valuable for the practitioner.

Published in Transactions on Machine Learning Research (07/2023)
A Unified Perspective on Natural Gradient Variational Inference with Gaussian Mixture Models
Oleg Arenz
oleg. datg.arren@tttuddmmsttdd.dd de
Intelligent Autonomous Systems
Technical University of Dormstadt
Philipp Dahlinger
philiipp.dahligeer Skit edu
Autonomous Learning Robots Karlsruhe Institute of Technology
Zihan Ye
zihan. wbwn.ye@lu.darmsaldd. de
Artificial Intelligence Machine Learning Hessian Center for AT (hessian.Al) Technical University of Darmstadt
Michael Volpp
michael, michael.volpp Skit gerhard. neumann Bkit, edu
Skit. edu
Gerhard Neumann
Autonomous Learning Robots
Karlsruhe Institute of Technology
Reviewed on OpenReview: https: hptp//veeerrr.e.... forum? id= tLBjaX4
Abstract
Variational inference with Gaussian mixture models (GMMs) enables learning of highly tractable yet multi-modal approximations of intractable target distributions with up to a few hundred dimensions. The two currently most effective methods for GMM-based variational inference, VIPS and iBayes-GMM, both employ independent natural gradient updates for the individual components and their weights. We show for the first time, that their derived updates are equivalent, although their practical implementations and thcorctical guaranteos differ. We identify several design choices that distinguish both approaches, namely with respect to sample selection, natural gradient estimation, stepsize adaptation, and whether trust regions are enforced or the number of components adapted. We argue that for both approaches, the quality of the learned approximations can heavily suffer from the respective design choices: By updating the individual components using samples from the mixture model, iBayes-GMi GMM often fails to produce meaningful updates to low-weight components, and by using a zero-order method for estimating the natural gradient, VIPS scales badly to digher-dimenssinall problems. Furthermore, we show that information geomettrrr trust- regions (used by VIPS) are effective even when using first-order natural gradient estimates, and often outperform the improved Bayesian learning rule (iBLR) update used by iBayes- GMM. We systematically evaluate the effects of design choices and show that a hybrid approach significantly outperforms both prior works. Along with this work, we publish our highly modular and efficient implementation for natural gradient variational inference with Gaussian mixture models, which supports 432 different combinations of design choices, facil- itates the reproduction of all our experiments, and may prove valuable for the practitioner.

---

Published in Transactions on Machine Learning Research (07/ (07/2023)
1 Introduction
Many problems in machine learning involve inference from intractable distributions p(x), that might further only be known up to a normalizing constant Z, that is, p(x) 5x) For example when learning latent vari- able models, p(x) corresponds to the intractable distribution over the latent variable (Volpp et al., 2023); and in maximum entropy reinforcement learning, p(x) corresponds to the exponentiated return, p(x) = exp R(z) of trajectory x (Ziebart, 2010). Bayesian inference is another example, where the intractable, unnormalized target distribution p(x) corresponds to the product of prior and likelihood. However, whereas in Bayesian inference, there is particular interest in scaling to high-dimensional problems and large datasets, we stress that our work considers problems of moderate dimensionality of up to a few hundred of dimensions, where modeling the full covariance matrix of Gaussian distríbutions is still tractable. Important applications can be found, for example, in robotics where p(x) could be multimodal distribution over joint-configurations that reach the desired pose (Pignal et. al., 2020) or over collision-free molions that reach a given goal (Ewerlon et al., 2020)-, or in non- con-amortized variational inference for latent variable models (Volpp et al., 2023).
Variational inference (VT) aims to approximate the intractable target distribution p(x) by mcans of a tractable, parametric model [x)), with parameters 0. Variational inference is typically framed as the problem of maximizing the evidence lower bound objective (ELBO), which is equivalent to minimizing the reverse Kullback-Leible (KL) divergence (Kullback and Leibler, 1951), KL(gollp), between approximation Go and larget distribution p. The reverse KL divergence is principled choice for the optimization prob- lem, as it directly corresponds to the average amount of information (measured in nats), that samples from the approximation (x( contain for discriminating between the approximation and the target distribu- tion. In this work, we focus on a particular choice of variational distribution - Gaussian mixture model (x(x)= = E. 2o (0(0))0(x()) with weights (o( and Gaussian components (x The ELBO is then given by
((O) = Σ (o) / logū(x)dx + (gg)),
(1)
where (go) = problem of maximizing Equation 1 with respect to the parameters 0, which correspond to the weights g(o) and the mean and covariance matrix for each Gaussian component q(x(o).
- - fx go (x) log(1o (x)) denotes the entropy of the GMM. Hence, this work studies the
Gaussian mixture models are simple yet powerful choice for a model family since they can approximate arbitrary distributions when assuming sufficiently high number of components. Compared to more complex models, such as normalizing flows (Kobyzev et al., 2020), they are more interpretable and tractable since not only sampling and evaluating GMMs is cheap, but also marginalizations and certain expectations (of linear or quadratic functions) can be computed in closed form. Furthermore, the simplicity of GMMs allows for sample-efficiet learning algorithus, that is, algorithus that require relatively few evaluations of the target distribution for learning the model parameters O. GMM-based variational inference is, in particular, relevant for problem settings with up to a few hundred dimensions (Volpp et al., 2023; Ewerton et al., 2020; Pignat et al., 2020), where learning and storing of full covariance matrices is still tractable.
Arguably the two most effective algorithms for maximizing Equation 1, both apply independent natural gradient (NG) updates on each component as well as on the categorical distribution over weights (Arenz et al., 2018; Lin et al., 2019a). Yet, both algorithms were derived from a different perspective, have different theoretical guarantees, and even different objectives for the independent updates. Namely, iBayes- (Lin et al., 2019a; 2020) uses the original GMM objective (Eq. 1) for each independent update to perform natural gradient descent also with respect to the full mixture model, whereas VIPS (Arenz et al., 2018; 2020) uses lower bound for expectation- eapestationmaxiiizaatie maximization procedure, which yields independent objective functions for each component and the mixture weights Their approach can be shown to converge, even when the M-Step does not consist of single natural gradient updates. However, it was not yet proven, that their proposed procedure, which does use single NG steps, also performs natural gradient descent on the full mixture.
In this work, we will further explore the previous works iBayes-GMM (Lin et al., 2019a; 2020) and VIPS (Arenz et al., 2018; 2020), and use our findings derive a generalized method that outperforms both of them. In particular, we will present the following contributions.
2

---

Published in Transactions on Machine Learning Research (07/2023)
In Section 2, will review and compare the derived update equations of IBayesGMM (Section 2.1) and VTPS (Section 2.2) without considerations of implementation details and design choices. We will prove in Section 2.3 that these update equations are actually equivalent and, therefore, the methods only differ due to design choices not mandated by the derivations. By connecting these two previously separated lines of research, we improve our theoretical understanding of GMM-based VI.
In Section 3 we will present a general framework for learning GMM approximation based on the up- date oquations discussed in Section 2. Our framework uses seven modules to independently select design choices, for example, regarding how samples are selected, how natural gradients are estimated or how the learning rate or the number of components is adapted. For each design choice, we will review and compare the different options that have been used in prior works, and we will discuss potential limita- tions. For example, VIPS uses an inefficient zero-order method for estimating natural gradients, whereas IBayesGMM updates the individual components based on samples from current GMM approximation, which can prevent component with low weight from receiving meaningful updates.
In Section 4, we evaluate several design choices based highly modular implementation of our gen- eralized framework. We propose a novel combination of design choices and show that it significantly outperforms both prior methods. In particular, we combine KL-constrained trust regions, which have been popularized in the gradient- reinforcement learning setting (Peters et al., 2010; Schulman et al., 2015; Outo et al., 2021), with gradient-based estimates of the NG (Lin et al., 2019b), use samples from each component and adapt the number of components. Test problems are used from both prior works.
We release the open-source implementation of our generalized framework for GMM-based VI. Our imple- mentation allows each design choice to be set independently and outperforms the reference implementa- tions of iBayes-GMii GMM and VIPS when using their respective design choices. A separate reproducibility package contains the scripts we used for starting each experiment, including hyperparameter optimization.
2 Two Derivations for the Same Updates
Maximizing the ELBO in Eq. I with respect to the whole mixture model can be roduced to individual updates of the components and the weight distribution (Arenz et al., 2018; Lin et al., 2019a). VIPS (Arenz et al., 2018) achieves this reduction using a lower-bound objective for an expectation-maxiiizai maximization procedure. In contrast, Lin et al. (2019a) investigated the natural gradient of Eq. I with respect to the GMM parameters and showed that it can be estimated independently for the paramcters of the individual components and the parameters of the weight distribution. While both decompositions can be applied to larger classes of latent variable models, in the following we restrict our discussion to Gaussian mixture models. We will briefly review the different derivations in Section 2.1 and 2.2 before presenting our main thcoretical result in Section 2.3, namely, that the derived updates are equivalent. Throughout this section, deliberately avoid the discussion of design choices and implementation details that are not mandated by the derivations. Instead we will discuss the different options which may significantly affect the performance in Section 3.
2.1 iBayesGMM: Independently Computing the Natural Gradient
The ELBO objective (Eq. 1) could be straightforwarddy optimized using (vanilla) gradient descent, using the reparameterization trick (Kingma and Welling, 2014; Rezende et al., 2014) for obtaining the gradient with respect to the parameters of each component and the weights. However, compared to gradient de- scent, natural gradient descent (Amari, 1998) has been shown to be much more efficient for variational inference (Khan and Nielsen, 2018). Whereas gradient descent performs steepest descent subject to the con- straint of (infinitesimal) small (in the Euclidean sense) changes to the parameters, natural gradient descent performs steepest descent subject lo small changes to the underlying distribution (with respect to the Fisher information metric). The natural gradient can be obtained from the vanilla gradient by preconditioning it with the inverse Fisher information matrix (FIM), but explicitly computing the FIM is expensive. Instead, Khan and Nielsen (2018) have shown that the natural gradient with respect to the natural parameters of exponential family distribution (such as Gaussian) is given by the vanilla gradient with respect to the
3

---

Published in Transactions on Machine Learning Research (07/2023)
expectation parameters (please refer to Appendix D for additional background). However, GMMs do not belong to the exponential family and are, thus, not directly amenable to this procedure.
To derive efficient natural gradient updates for broader class of models, Lin et al. (2019a) considered latent variable models, such as GMMs, where the marginal distribution of the latent variable q(o) and the condi- tional distribution of the observed variable g(xlo) are both minimal exponential family distributions. They showed that for such minimal conditionally exponential family (MCEF) distributions, the Fisher informa- tion matrix of the joint distribution g(2, o) is block-diagonal, which in turn justifies computing the natural gradients of the individual components and the weight distribution independently.
2.2 VIPS: Independently Maximizing Lower Bound
Whereas Lin et al. (2019a) showed that single natural gradient update can be performed for every com- ponent independently and that such procedure also performs natural gradient descent the whole mixture model, Arenz et al. (2018) proposed a method for GMM-based variational inference that is not derived for natural gradient descent and that allows for independent optimizations (going beyond single step updates).
For understanding how VIPS (Arenz et al., 2018) decomposed the ELBO objective (Eq. 1) into independent objectives, is helpful to first understand which terms prevent independent optimization in the first place. In particular, note that the first term of the ELBO given in Equation 1 already decomposes into indepen- dent objectives for cach component, and only the second term (the model entropy) prevents independent optimization. More specifically, when using Bayes' theorem to write the probability density of the GMM in terms of the marginals and conditionals, the interdependence between the different components can be narrowed down to the rrogosresposssiiiiiis (o(x)x) within the entropy,
H(go) = - - f (x()) (log
go (0(0)a0(x10) )ax (0
which is the only term in Equation 1 that creates a mutual dependence between different components. Hence, Arenz al. (2018) introduced an auxiliary distribution q(o|x), to derive the lower bound
1) = Ego(o) I/ (o(x|o) (logū(x) (10g())))))F( + log (((x) dx + H((o(x)o)) H(o(o))
(2)
=](() [ [
As an expected KL is non-negative for any two distributions and equals zero, if both distributions are equal, J(q. 0) is a lower bound of J(0) that is tight when ī(o|x) ()(o))). Hence, VIPS uses procedure similar to expectation maximization: In the E-step of iteration the auxiliary distribution is set to the current model, that is, ī(olx) == aatn (o|x). Then, in the M-step, the model parameters are optimized with respect to the corresponding lower bound [ggoo, 0). As the lower bound is tight before each M-step, i., 3[ (()) = J(O()), it also improves the original objective. As the auxiliary distribution ?(o|x) is independent of the model parameters (during the M-step), the weights (o(o) and each component @o (x(o) can be optimized independently, resulting in a component-wise loss function
Ja(ğ. 0) = (x(())) (log((x) + log ((olx) log go (x)o) dx.
While these derivations would justify any procedure that improves the model with respect to the lower bound objective during the M-step, VIPS, performs single natural gradient steps, closely connecting it to Bayc Indeed, we will now show that the updates of both mothods are equivalent, except for design choices that are not mandated by the derivations and can be interchanged arbitrarily.
2.3 The Equivalence of Both Updates
Both previously described methods iteratively improve the GMM by applying at each iteration a single NG step independently to each component and the weight distribution. Although the NG is computed with
4

---

Published in Transactions on Machine Learning Research (07/2023)
respect to different objective functions, namely the original ELBO (Eq. 1) and the lower bound (Eq. 2), we clarily in Theorem 2.1 that both natural gradients, and, thus, the derived updates, are indeed the same.
Theorem 2.1. Directly after the E-step (which sets the auriliary distribution to the current model with porameters o that is ((o[x)= = agt (olx)), the natural gradient (denoted by V) of the lower bound objective (Eq. 2) matches the natural gradient of the ELBO (Eq. 1), VoJ(O) (00
(gosss 0=0
Proof. See Appendix E.
This deep connection between both methods was previously not understood and has several implications. First of all, shows that VIPS actually performs natural gradient descent on the whole mixture model (not only on the independent components) and with respect to the original ELBO objective (not only to the lower bound). On the other hand, it may be interesting to exploit the milder convergence guarantees of VIPS (even bad NG estimates may lead to convergences long perform more robust natural gradient descent. Increasing the robustness by monitoring the lower bound may in particular become relevant when it is difficult to approximate the natural gradient sufficiently well, for example in amortized variational inference, where we want to learn a neural network that predicts the mean and covariances based on a given input. Finally, our new insights raise the question, why the reference implementations of the different methods perform quite differently, despite using the same underlying update equations. In the remainder of this work, we will shed some light on the latter question, by analyzing the different design choices and evaluating their effects on the quality of the learned approximations.
improve on the lower bound objective) to
3 A Modular and General Framework
We have shown that both derivations suggest exactly the same updates (although the decision to perform single natural gradient steps is voluntary for the EM-based derivation). However, for realizing these up- dates, several design choices can be made, for example, how to estimate the natural gradient when updating a Gaussian component. By closely comparing both previous methods and their respective implementations, we identified seven design choices that distinguish them. Hence, we can unify both methods in a common framework (Algorithm 1) with seven modules that can be implemented depending (1) The SampleSelector selects the samples that are used during the current iteration for estimating the natural gradients with respect to the parameters of the components and the weight distribution, (2) the SamppoeenSlepsiieAAuulloe
the design choice:
module chooses the stepsizes for the next component updates, (3) the NgEs- timator estimates the natural gradient for the component update based on the selected samples, (4) the Ng ByBaseaUpdlater performs the component updates based on the estimated natural gradients and the step- sizes, (5) the WeighttStepsizeeAdaatttti module chooses the stepsize for the next update of the weights, (6) the WeighlUpdater updates the weight distribution along the natural gradient based on the chosen stepsize, and (7) the ComponentAdaptatuinn module decides whether to add delete components. In the remainder of this section, we will discuss and compare the different options that have been employed by prior work.
A Algorithm Natural Gradient GMM Variational Inference repeat
samples, largedDensiiiessA
SampleSleetorseeeeeSSnpleeeeee SampleDatabase)
compStepsizes <- ComppreeeSSeeesizeetttpptttooo.u wttp::/ioouuuutdttttesssseeeeeee (mum (nmmCompUpdtteettttoo compStepsizes)
naturalGradients + NCFFstinator.gee/NSssssssaa/ aereeeDDDensitessnnnnaadada
GMM)
GMM + wwtpe//hheeelie....ooeee/aaaaaaa
stepsizes, GMM)
weight.St.essiee + Wwg.eegggieegiiiilaaaaaaaoooooia
GMM + Weiienntupplaaee...uuuaeeeeeeii htt.://wwwwtttpppssssoemeeaaaooo (nmmteratioons, weight.Stepsize)
samples, tangetDDnsitties, GMM)
GMM + C-t.0-0-aaaa.......aaaaaaaaaaaaa
GMM, SampleDatabase)
until stoppingCriteria
5

---

Published in Transactions on Machine Learning Research (07/ 2003)
3.1 Sample Selection
For estimating the natural gradients of the components and the weight distribution, the target distribution p(x) needs to be evaluated on samples from the respective distribution. However, both GMyyes-GMM and VIPS use importance sampling to use the same samples for updating all distributions at a given iteration.
iBayesGMM obtains the samples for a given ileration by sampling from the current. GMM, using standard importance sampling to share samples between different components.
VIPS samples from the individual components, rather than the GMM and also stores previous samples (and their function evaluations) in a database. Using two additional hyperparameters, Nreused and deee: VIPS starts by obtaining the Preused newest samples from the databasc. Then, the effective sample size net(()) = (Ei (o))) is computed using the self-normalized importance weights wft and Oldes - (ott(() new samples are drawn from each component.
In our implementation both oplions con make use of somples from previous ileralions by selling the by- perporameter neressed larger than zero (Arenz et al., 2020). The first option, which generalizes the procedure used by iBaaee-GGMi GMM (Lin et al., 2020) computes the effective sample size Theff on the GMM and draws max(0, N des nefy) new samples from the GMM. The second option computes neg(o) for each component and draws max (0, Rides - neg(()) new samples from each component, matching the procedure used by VIPS.
3.2 Natural Gradient Estimation
For estimating the natural gradient for the individual component updates, VIPS only uses black-box evalu- ations of the unnormalized target distribution log ((x;) on samples x, In contrast, the NG estimate used by iBayes-GMM (Lin et al., 2020), which is based on Stein's Lemma (Stein, 1981), uses first-order information Vx, log p(x,) which is less general but typically more sample efficient.
VIPS uses the policy search method MORE (Abdolmaleki et al., 2015), which is based on compatible function approximation (Pajarinen et al., 2019; Peters and Schaal, 2008; Sutton et al., 1999). Namely, as shown by Peters and Schaal (2008), an unbiased estimate of the natural gradient for an objective of the form ER (R()] is given by the weights w of a compatible function approximator (Sutton et al., 1999), R(x)= = w Φ(x), that is fitted using ordinary least squares to approximate R, based on data set X ((x, R(x)),,, with samples x, that are obtained by sampling from the distribution ma. A function approximator is compatible to the distribution тe, if the basis functions ((x) are given by the gradient of the log distribution, Φ(x) Vo log TO (x)). For updating Gaussian component g(x(o) parameterized with its natural parameters, Vo = [E Mo: -11-1.
1), the compatible function approximator can be written as
n(x) = x T Rx xxx+xx+r,
where the matrix R, the vector r and the scalar r, are the linear parameters that are learned using least squares. Here, the constant offset r can be discarded and R and r directly correspond to the natural gradients, that could be used to update the natural parameters,
1 -. 0 -
---1 old + +B.R.
Σ So'Mo = E Holdd + Bar.
Ho
(4)
2
The least squares targets basod on the ELBO (Eq. 1), (x) = p(x) log go (x), and the least squares targets based on the lower bound (Eq. 2), Rx, = log p(x) +log q(o) x) - log q(x(o) only differ by constant offset log go(o) (o) that would be consumed by the scalar and not affect the NG estimates R and r.
For compatible function approximation, the data set X' should be obtained by sampling the respective component. However, VIPS uses importance weighting to make use of samples from previous iterations and from different components, by performing weighted least squares, weighting each data point by its importance weight. The importance weights wA can be computed as ((x]o)
2(x;) where 2 is the actual distribution that was used for sampling the data set. However, VIPS uses self-normalized importance weights with = (E, - which yield lowervarrincce, biased, but asymptotically unbiased estimates. Furthermore, VIPS applics ridge regularization during weighted least squares, which further introduces a bias
6

---

Published in Transactions on Machine Learning Research (07/2023)
iBayesGMM exploits that for exponential-family distributions, the natural gradient with respect to the natural parameters (for Gaussians = - - - corresponds to the vanilla gradient with respect to the expectation parameters m (for Gaussians m = {u, Σ + ph }) (Khan and Nielsen, 2018). Using the chain rulc, the gradient with respect to the expectation paramcters can be expressed terms of the gradient with respect. lo the covariance Σ and mean l (Khan and Lin, 2017, Appendix B.1). Thus, the NG step for a Gaussian component q(x(o) with stepsize 3. and objective J(x) can be computed as
1 Σ - - 2
I- o,old + VV..J E no = - ooood/ooool + Bo (-2
V Ma,old + Vnr).
(5)
as shown by Khan et al. (2018, Appendix C). As the objective corresponds to an expected value, i.e., J = (x(x R(x) with R(x) = log f(x) expected gradient and Hessian (Opper and Archambeau, 2009),
[(x) the gradient with respect to mean and covariance are given by the
- (xt [xxx((x)) V = 2
V = Bo(x(o) Vx
(6)
Hence, using Monto-Carlo to estimate the gradients with respect to mean and covariance (Eq. 6), we can obtain unbiased estimates of the natural gradient (Eq. 5). However, evaluating the Hessian Vxx R(x) is computationally expensive, and therefore, Lin et al. (2019b) proposed an unbiased first-order estimate of the expected Hessian based on Stein's Lemma (Stcin, 1981) given by (Lin et al., 2019b, Lemma 11)
(x(x(o) V
Xxx Xxx) Eg(x)o [2, (x Mo) (VxR(x))
(7)
Similar to VIPS, (Bayes-GMM uses importance sampling to perform the Monte- Carlo estimates (Eq. 6 right, Eq. 7) based on samples from different components. However, in contrast to VIPS, iBayes-GMM uses standard (unbiased) importance weighting, rather than self- self-normalized importance weighting, which in our implementation, can be selected using a hyperparameter that available for both options.
Our implementation supports both options, where using standard importance weighting or self-normalized imporlance weighting can be chosen using a hyperparameler that is available to both oplions.
3.3 Natural Gradient based Component Updates
For performing the natural gradient update, we identified three different options in the related literature. Lin et al. (2019a) directly apply the natural gradient update (Eq. 4) based on the natural gradicnts R, and the stepsize Bo, which we assume given. However, this update may lead to indefinite covariance matrices, and therefore iBayes-GMM (Lin et al., 2020) uses the improved Bayesian learning rule (tBLR), which modifies the update direction (no longer performing natural gradient descent) in a way that ensures that the covariance matrix remains positive definite. Thirdly, VIPS indirectly control the stepsize by choosing an upper bound the KT divergence between the updated and the old component, [[[[(x optimization problem to find the largest & that satisfies the trust region constraint. Directly performing the natural gradient step is straightforward, and therefore, now only discuss the two latter options.
and solves a convex
iBayesGMM applies the improved Bayesian Learning Rule (iBLR), which approximates Riemannian gra- dient descent (Lin et al., 2020).
- - - E-
Σ
B.REe, 2 oolid + + (R-.RRREE R),
Σ = - Eooold
I'''D (Eoda + Bor)
(8)
2
The iBLR update (Eq. 8) differs from the NG update (Eq. 4) due to the additional terms - R.R.E. R and Σ- Doold Using this update, the new precision matrix can be shown to be an average of a positive definite
and a positive semidefinite term, Σ-¹ E + UTU), where U - 2.3L and L = ooold
VIPS updates the components using natural gradient descent (Eq. 4). However, rather than controlling the change of the component using the stepsize So, it upper- upper-bounds bounds the KL divergence KL (gx(x)))l((x(())) < c, by solving for every component update an optimization problem that finds the largest stepsize Ba, such that
7

---

Published in Transactions on Machine Learning Research (07/)023)
the updated component will have a positive definite covariance matrix and KL divergence to the original component that is smaller than a bound 6o, which is assumed given instead of the stepsize 8o (Abdolmaleki et al., 2015). The additional optimization problems for finding the stepsize add little computational overhead since the KL divergence is convex in the stepsize and can be efficiently computed for Gaussian distributions.
Our implementation supports three options: (1) directly applying the NG update (Eq. 4), (2) the iBLR update, and the (3) trust-region update of VIPS. However, VIPS uses an external L-BFGS-G optimizer to find the largest stepsize that respects the trust-region, which was difficult to efficiently integrate into our Tensorflow (Abadi et al,, 2015) implementation. Instead we implemented a simple bisection method so solve this convex oplimization problem, see Appendix I for details.
3.4 Weight Update
The natural gradient update for the categorical distribution is performed by updating the log weights in the direction of the expected reward of the corresponding component, that is,
q()) x q())oia exp (BR(O)).
(9)
The NG update is presented differently by Lin et al. (2020), but clarify in Appendix F that both updates are indeed equivalent'. However, we identified two different options to perform the weight update.
where R(o) is a Monte-Carle estimate of R(o) Ba(ojo) [log p(x) {(x)
IBayesGMM directly applies the stepsize g.
VIPS originally optimized the stepsize with respect to a desired KL bound KL (g(o)) (ota)))))))))) S e, analogously to the component update described in Section 3.3 (Arenz et al., 2018). Although in subsequent work, Arenz et al. (2020) found that fixed stepsize 3 = 1, corresponding a greedy step to the optimum of their lower bound objective, performs similarly well.
Our implementation supports both options, directly applying the stepsize and optimizing the stepsize with respect lo trust-region. Directly opplying the slepsize in combination with a AespisSStesizeAaapttioo that always selects fixed stepsize 3 = corresponds to the greedy procedure (Arena et al., 2020).
3.5 Weight and Component Stepsize Adaptation
While we use separate modules for Weightsteesiize and ComporentStpsiiee.dppttttt. we consider the same adaplation schemes in both modules.
iBayesGMM used fixed stepsizes in their implementation but used different stepsizes for weights and for the components. Often, a weight stepsize 8 = was used to maintain a uniform distribution over the different components. We believe that this concession was probably required due to their design choice, of using samples from the GMM for updating the individual components, which may prevent components with low weight to obtain meaningful updates. In prior work, Khan and Lin (2017) discussed backtracking method to adapt the stepsize for the component updates for ensuring positive definiteness, but they found it slow in practice and did not test it in the GMM setting, and Khan et al. (2018) used a decaying stepsize.
VIPS originally used fixed trust-regions for the weight and component update (Arenz et al., 2018). Later a fixed stepsize with updates, where the bounds were chosen independently for each component, by increasing the bound if the reward of the corresponding component improved during the last update, and by decrcasing it otherwise.
for the weights was combined with adaptive trust-regions for the component
Our implementation supports three options for each stepsize scheduler: (1) keeping the stepsize constant (2) using a decaying stepsize (Khan et al. 2018) and (3) adapting the stepsize based on the improvement of the last updote ( (Arenz et. al., 2020).
The respective implementations still differ slightly because VIPS uses nolf-normalized importance weighting, rather than standard importaane--eeighinngg which we control in our implementation using a hyperparameter.
8

---

Published in Transactions on Machine Learning Research (07 (07/2023)
MODULE
OPTIONS
NgEstimator
More
Z
Stein
ComponentAAaaptttion
Fixed
E
Vips
SampleSelector
Lin et al.
P
Vips
M
DyBassdCCmmponennUppaae
Direct
I
iBLR
Y
Trust-Region T
Fixed
F
Decaying
D
Adaptive
R
Weight WeighiSleesizeeAkäalltiiin
Direct
U
Trust-Region
O
Fixed
X
Decaying
G
Adaptive
N
Table 1: We assign unique letter to every option such that every combination of options can be specified with a 7-letter word (one letter per module).
3.6 Component Adaptation
VIPS dynamically adapts the number of components by deleting components at poor local optima that do not contribute to the approximation, and by adding new components promising regions. The initial mean for the new component is chosen based on a single-sample estimate of the initial reward R(Onew) for the new component, evaluated for samples from the database.
iBayes-GMM uses a fixed number of components K that is specified by hyperparameter.
Our implementation supports (1) a dummy option that does not do anything and (2) adapting the number of components dynamically (Arenz el al., 2020).
4 Experiments
To evaluate the effect the individual design choices, we implemented all previously discussed options in a common framework, such that they can be chosen independently. Prior to the experiments we ensured that our implementation performs comparably to the MATLAB implementation of iBayes-GMM (Lin et al., 2020) and to the C++-implemennatioon of VIPS (Arenz et al., 2020) when using the respective design choices. As shown in Appendix G, when comparing our implementation with the reference implementation on their target distributions, we always learn similar or better approximations, Additional details on the implementation are described in Appendix I. An overview of the available options shown in Table 1. Available hyperparamoters for each option are shown Appendix H. Our framework allows for 24 3 = 432 different combinations of options. We assign a unique letter each option (see Table 1), such that we can refer each of the combinations with 7-letter codeword. For example, SEPYFUX refers to iBayes-GMM (Lin et al., 2020) and ZAMTRUX to VIPS (Arenz et al., 2020). However, evaluating each combination is computationally prohibitive. Instead, we designed the following sequence of experiments.
In the first group of experiments, we evaluate the stability of the component updates, using fixed weights, and without adapting the number of components. In the second group of experiments, we use the most promising design choices for the component update that we identified in the first group of experiments, and ablate over the design choices that affect the update of the weights and the adaptation of the number of components. Based the results of the previous group of experiments, we then identify promising combinations and evaluate them the full test suite. For comparison with prior work, we also evaluate the design choices used by iBayes-GMM (Lin et al., 2020) and VIPS (Arenz et al., 2020).
4.1 Experiment 1: Component Update Stability
For evaluating the effects of the design choices on the component update stability, we do not change the number of components (Option E) and do not update the (uniformly initialized) weights (X and U, with initial weight-stepsize 3 = 0). Furthermore, we always sample from the GMM (P), since the two options for the SampleSelector hardly differ for uniform GMMs.. We evaluate the different options for the NgEstimator, ComponenuSteppizzeAttuau Adoplation and NgBosedCoompnneen Updater resulting in 18 different combinations.

---

Published in Transactions on Machine Learning Research (07/2023)
Design Choice
BreastCancer
Breast CrecerMM
Wine
MORE (2)
78.88
152.84
N/A
Component NG Estimation
Stein (S)
78.46
68.22
1431.09
Direct (1)
78.95
70.96
2263.78
Component Update
iBLR (Y)
78.69
70.37
1431.09
Trust-Region (T)
78.46
68.22
1443.40
Fixed (F)
78.46
69.01
1450.11
Component Stepsize Adaptation
Decaying (D)
78.47
70.69
1468.50
Adaptive (R)
78.46
68.22
1431.09
We show optimistic estimates of the best performance (negated ELBO) that we can achieve with every option when optimizing a uniform GMM with a fixed number components, by reporting the best performance achieved during Bayesian optimization among all runs that used the respective option.
Table 2:
We evaluate each candidate three target distributions: In the BreastCancer experiment Arenz et al. (2018) perform Bayesian logistic regression using the full "Breast Cancer" datatset (Lichman, 2013). Further, we introduce the variant. BreastCanaerMMB which uses minibatches of size 64 to evaluate the effect of stochastic target distributions. Thirdly, we evaluate the candidates in a more challenging Bayesian neural network experiment on the "Wine" dataset (Lichman, 2013), where we use a batch size of 128. We use two hidden layer with width 8, such that the 177-dimensional posterior distribution is still amenable to GMMs with full covariance matrices. Please refer to Appendix for details on the implementation for all experiments.
We tuned the hyperparameters of each of the 18 candidates separately for each test problem using Bayesian optimization. We granted each candidate exclusive access to compute node with 96 cores for two days per experiment. The results of the hyperparameter search are summarized in Table 2, where we report the best ELBO achieved for every design choice (maximizing over all tested hyperparameters and all combinations of the remaining options). On Wine, which uses almost three times as many parameters as have been tested by Arenz et al. (2020), we could not obtain reasonable performance when using MORE, as this zero-order method requires significantly more samples, resulting in very slow optimization.
The values in Table 2 give optimistic estimate of the best performance that we expect for given design choice, however, it is in general not guaranteed that this performance can be robustly achieved over different seeds. Hence, for the most promising candidates, SEPTRUX and SEPYRUX, we used the best hyper- parameters from the Bayesian optimization and evaluated the performance over ten seeds. The mean of the final performance and its 99.7% standard error are shown in Table 3. We also show for every experiment a second metric, where we use the maxinum mean discrepancy (Gretton et al., 2012) (MMD) for BreastCancer and BreastCancerMM (comparing samples from the GMM with baseline samples (Arenz et al., 2020)), and the mean squared error of the Bayesian inference predictions for Wine. Unfortunately, the hyperparameters for SEPYRUX BreastCancer and Cancer led to unstable optinization for some seeds, despite their good performance during hyperparameter search. We expect that using slightly more conservative hy- perparameters, SEPYRUX could reliably achieve a performance that is only slightly worse than the optimistic value provided in Table 3 for BreastCancer (78.69). However, the performance of SEPTRUX is even in expec- (ation over multiple seeds already better than this optimistic estimate. Furthermore, also on WINE, where SEPYRUX did not suffer from instabilities, the trust-region updates achieved better performance, albeit not statistically significant. Hence, we decided to use trust-region updates for the second group of experiments. We will further test the iBLR update in main experiments for more thorough comparison.
4.2 Experiment 2: Weight Update and Exploration
According to the first experiment, first-order estimates using Stein's Lemma, adaplive component stepsizes and trust-region updates are the most effective and stable options for the component updates, and therefore, we fix the corresponding design choices (S, T, R.) for our second group of experiments. The modules that we SampleSelector, WeightUpdater
did not evaluate in the first group of experiments (ComponentAAdataatin and WeighlSSepssizeAdapllaaa are in particular relevant for discovering and approximating multiple modes
10

---

Published in Transactions on Machine Learning Research (07/2023)
BreastCancer
Breast CancerrB
Wine
Candidate
-ELBO
-ELBO
MMD
MMD
-ELBO
MSE
(full batch)
1042.28
0.004
1130.81
0.033
1462.91
0.481
SEPYRUX 26699.32
#0.002
7799.66
#0.014
335.70
#0.02
78.53
0.002
81.21
0.002
1444.01
0.478
SEPTRUX
=0.02
+00000
#1.13
10.000
+30.78
#E0.221
We evaluated the best hyperparameters for the most promising candidates of our experiments for Group on 10 different seeds with respect to ELBO and secondary metric (maximum mean discrepancy mean squared error). SEPYRUX, which uses the iBLR update for the component updates did not achieve stable results the breast cancer experiments. SEPTRUX which uses trust-region updates for the components oulperformed SEPYRUX in all experiments, although in WINE the advantage is not statistically significant.
Table 3:
Design Choice
GMM20
STM20
PlanarRobot4
Non-Adaptive (E)
0.00
0.08
12.33
Component Adaptation
Adaptive (A)
0.00
0.05
12.15
From Mixture (P)
0.11
1.52
12.77
Sample Selection
From Components (M)
0.00
0.05
12.15
Natural Gradient Descent (U)
0.00
0.06
12.30
Weight Update
Trust-Regite (0)
0.00
0.05
12.15
Fixed (X)
0.00
0.06
12.15
Weight Stepsize Adaptation
Decaying (G)
0.00
0.05
12.33
Improvement-Based (N)
0.00
0.06
12.22
Table 4: We show optimistic estimates of the best performance (negated ELBO) that we can achieve with
every option (updating the components using the design choices identified in the first experiment). We report the best performance achieved during Bayesian optimization among all runs that used the respective option.
of the target distribution. Hence, we focus on multi-modal target distributions for the second set of experi- ments. All test problems for these experiments were taken from prior work. Namely, we chose GMM20 and PlanarRobot, from Arenz et al. (2020) and STM20 from Lin et al. (2020). For GMM20 and STM20 the target distribution is given by an unknown mixture of 20-dimensional Gaussians and Student- respec- tively. In PlanarRobot, we want to approximate distribution over joint configurations of 10-link planar robol, such that it approximately reaches any of four possible goal positions. For Bayesian optimization of the hyperparameters, we grant each of the 24 candidates exclusive access to our compute node for one day per test problem. Optimistic estimates of the best performance for each design choices are shown Table Based on our experiments sampling according to the mixture weights (Option P) seems to be clearly inferior to sampling from the components (Option M), it the only option that. was not able to solve the GMM20 experiment. Furthermore, adapting the number of components (Option A) and using trust-region updates for the weight update (Option O) seem beneficial for multimodal target distributions. However, for WeighhtStepsizeeAuapttte also a fixed stepsize (when used as trust-region) achieved good performance.
4.3 Experiment 3: Evaluating the Promising Candidates
For main experiment focus candidates that use first order NG estimates, adaptive number of components and adaptive component-stepsizes and sample from the individual components (Options S, A, M, R), and aim to better compare trust-region updates with the iBLR update (Option T vs. Y). For updating the weights, we evaluate fixed and adaptive trust-regions, and fixed direct NG updates (OX vs. ON vs. UX), resulting in candidates: SAMTRUX, SAMTROX, SAMTRON, SAMYRUX, SAMYROX and SAMYRON. Furthermore, compare with prior work, we also evaluate ZAMTRUX (Arenz et al., 2020) and SEPYPUX (Lin et al., 2020) well as the variant SEPYRUX, combining iBLR updates with adaptive stepsizes.
11

---

Published in Transactions on Machine Learning Research (07 (07/2023)
GermanCrediLMM
BreatCancer.BB
GermanCredit
PlanarRobol
BreastCancee
GMM100
GMM20
STM20
STM300
TALOS
WIND
Candidate
78.00
78.41
585.10
585.12
11.47
-0.00
0.01
0.00
14.96
1423.12
224.43
SAMTRON
+0.02
+0.03
+0.00
+0000
+0.04
+0.00
+0.03
+0.00
+0.48
+31.55
+0.16
78.61
79.96
585.11
585.12
12.93
0.00
0.16
0.03
22.50
1433.87
244.00
SAMYRON
+0.05
0.22
#0.01
+0.00
#0.13
+0.00
0.00
#0.01
+0.15
330.49
+0.23
iBayes-GMar GMM
79.78
80.13
585.12
585.12
17.26
0.43
4.54
0.46
26.87
3279.32
16.64
(SEPYFUX)
+0.40
0.66
+0.01
+0.00
+2.13
10.15
+0.52
#0.06
+0.45
11425.12
2526
VTPS
78.14
83.65
585.10
585.35
11.48
-0.00
1.21
0.54
163.3338
23.69
N/A
(ZAMTRUX)
+0.01
0.99
+0.00
+0.03
+0.04
+0000
+0.15
+0017
881.755
#0.16
Table 5: Along with the negated ELBO,, we show the 30 confidence intervals based the standard er- ror of its mean using ten different seeds. The proposed candidate clearly outperforms the prior methods VIPS (Arenz et al., 2020) and IBAYES-GMM (Lin et al., 2020). First-order natural gradient estimates with trust-region constraints (T) seem preferable over the iBLR update (Y). We observed instabilities for SEPYFUX on Planar Robot and TALOS and, thus, removed bad outliers when computing the reported values.
We evaluate these candidates on the full test suite, which uses the following test problems top of the previously discussed ones: GermanCredit (Arenz et al., 2020) and GermanCreditMB are similar to the BreastCancer experiment, but use the 25-dimensional GermanCreddi dataset (Lichman, 2013); GMM100 and STM300 are higher-dimensional variants of the GMM20 and STM20 experiments; and, finally TALOS is a new experiment that we introduced to evaluate the different options on test problem that was not used during development. Somewhat related to the PlanarRobot4 experiment, in the TALOS experiment we aim to learn distribution over joint configurations to reach goal position with the robot's endeffector. However, the TALOS experiment uses the kinematics of actual robot the humanoid TALOS from PAL Robotics (Stasse et al., 2017). Furthermore, the target distribution is based on the work by Pignat et al. (2020) and is more complex by using a product of expert that penalize pose errors for each foot and hand as well as unstable configurations (that occur when the robot's center mass projected on the ground is outside of the support polygon spanned by the feet), and that incorporate a prior distribution over the joint angles. The target distribution 34 dimensional (7 joint configurations for each leg, 6 joint angles for each and 6 additional parameters that specify the pose of the torso with respect to fixed reference frame. On TALOS, prior to the experiments of this group, we only performed few experiments using SEPIFOX to optimize single Gaussian component, to ensure that the target distribution is correctly implemented.
For selecting the hyperparameters, we perform for each candidate and each experiment a small grid search, where we make use of the results from the previous experiments to select suitable ranges. Extensive hyperpa- rameter search (e.g. using Bayesian optimization) would unfairly benefit options with more hyperparameters. Table compares the final performance (negated ELBO) of the best-performing candidate (SAMtRON) with the prior methods VIPS (Arenz et al., 2020) and iBayes-G.M GMM (Lin et al., 2020). We also show SAMYRON for comparing trust-region natural-gradient updates with the iBLR update. According to these experiments, we can clearly improve upon ZAMTRUX by using Stein's 's Lemma for estimating the natural gradient (in partic- ular for higher-dimenssinall problems)) and upon SEPIFUX, by sampling from the components and adapting their number during optimization. Interestingly, trust-region constraints seem to be beneficial also when using first-order estimates of the natural gradient, and showed slight but consistent advantage compared to the iBLR update (Lin et al., 2020) in our experiments. Full results of our main experiments are shown in A Appendix K, where we show the performance of all tested candidates, also with respect to secondary metrics.
12

---

Published in Transactions on Machine Learning Research (07 (022))
5 Conclusion
Although VIPS and iBayes- GMM are derived from different perspectives where the derivations for Bayes- GMM are less general (by requiring single NG steps for the component update) but enjoy stronger guaranlees the whole mixture model)-, we showed that both algorithms only differ in design choices and could have been derived from the other perspective, respectively. This unification of both perspective shows that we can derive approximate natural gradient descent algorithms also for mixtures of non-Gaussian components- where the approximation errors of the natural gradient are potentially much larger without having give up on convergence guarantees. Furthermore, our results are of high relevance for the practitioner, both due our extensive study on the effects of the individual design choices- which shows that both prior works can be improved by using combination of their design choices-and by relcasing our modular framework for natural gradient GMM-based variational inference, which is well-documented and easy to use and outperfoms the reference implementations by Arenz et al. (2020) and Lin et al. (2020) when using the respective design choices. Typical trade-offs for the practitioner are discussed in Appendix M. The limitations of this study and the potential for negative societal impact are discussed Appendix A and B.
(by proving natural gradient descent
Acknowledgements
This research supported by "The Adaptive Mind", funded by the Excellence Program of the Hessian Ministry of Higher Education, Science, Research and Art.
The authors gratefully acknowledge the computing time provided to them on the high-performance computer Lichtenberg at the NHR Centers NHRACES at TU Darmstadt. This funded by the Federal Ministry of Education and Research, and the state governments participating on the basis of the resolutions of the GWK for national high performance computing at universities.
The authors lednowledge support by the state of Haden-Württemilee through bwHPC.
This work was performed on the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden-Würteenn Württemberg and by the Federal Ministry of Education and Research.
References
Michael Volpp, Philipp Dahlinger, Philipp Becker, Christian Daniel, and Gerhard Neumann. Accurate bayesian meta-learning by accurate task posterior inference. In International Conference on Learning Representations, 2023. URL http:://ppppaaeeee..eeeeerooorr. -ksarrd0
Brian D. Zicbart. Modeling purposeful adaptive behavior with the principle of mazimum causal entropy. PhD thesis, Carnegie Mellon University, 2010.
Emnmanuel Pignat, Teguh Lembono, and Sylvain Calinon. Variational inference with mixture model ap- proximation for applications in robotics. In International Conference on Robotics and Automation, pages 8395-3401. IEEE, 2020.
Marco Ewerton, Oleg Arenz, and Jan Peters. Assisted teleoperation in changing environments with mixture of virtual guides. Advanced Robotics, 34(18):115--1170, 2020. doi: 10. 10.10800/11 00000066666...... 012.6966..222.177888. 1780.1785326.
Solomon Kullback and Richard A Leibler. On information and sufficiency. slalistics, 22(1):79 86, 1951.
The annals of mathematical
Ivan Kobyzev, Simon Prince, and Marcus Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions Pattern Analysis and Machine Intelligence, 2020.
Oleg Arenz, Mingjun Zhong, and Gerhard Neumann. Efficient gradient-free variational inference using policy search. Tn International Conference Machine Learning, 2018.
13

---

Published in Transactions on Machine Learning Research (07/2023)
Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt. Fast and simple natural-gradient variational inference with mixture of exponential-family approximations. In International Conference on Machine Learning, pages 3992 4002. PMLR, 2019a.
Wu Lin, Mark Schmidt, and Mohammad Emtiyaz Khan. Handling the positive-defiitte constraint in the Bayesian learning rule. In Hal Daumé III and Aarti Singh, editors, International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 6116-6126. PMLR, 13- 18 Jul 2020.
Oleg Arenz, Mingjun Zhong, and Gerhard Neumann. Trust-region variational inference with gaussian mixture models. Journal of Machine Learning Research, 21(163):1-60 2020. URL httt::/jmmzz org httg/apprr//21// 29-524.05 html. htmm2..mmm.
Jan Peters, Katharina Mulling, and Yasemin Altun. Relative entropy policy search. In Conference on Artificial Intelligence. AAAI, 2010.
John Schulman, Sergey Levine, Pieter Abbeel, Michael Jordan, and Philipp Moritz. Trust region policy optimization. In International Conference on Machine Learning, pages 1889 1897. PMLR, 2015.
Fabian Otto, Philipp Becker, Vien Anh Ngo, Hanna Carolin Maria Ziesche, and Gerhard Neumann. Dif- ferentiable trust region layers for deep reinforcement learning. In International Conference on Learning Representations, 2021. URL https://openrevier: wttp:eerrreeiee.eeeoooooooreiiie hetp/oorumiddgg//--Aaa.a A01Vn.
Wu Lin, Mohammad Emtiyaz Khan, and Mark Schmidt. Stein's lemma for the reparameterixation trick with exponential family mixtures. arXiv preprint arXi::1910.1383 2019b.
Diederik Kingma and Max Welling. Auto-encoding variational bayes. Tn International Conference on Learn- ing Representations (ICLR), 2014.
Danicl Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate in- ference in deep generative models. In International Conference on Machine Learning, pages 1278 1286, 2014.
Shun-Ichi Amari. Natural gradient works efficiently in learning. Neural compulation, 10(2):251-22) 276, 1998.
Mohammad Emtiyaz Khan and Didrik Nielsen. Fast yet simple natural-gradient descent for variational inference in complex models. In International Symposium on Information Theory and Its Applications (ISITA), pages 31 35. IEEE, 2018.
Charles M. Stein. Estimation of the mcan of a multivariate normal distribution. The annals of Statistics, pages 1135 1151, 1981
Abbas Abdolmaleki, Rudolf Lioutikov, Nuno Lau, Luis Paulo Reis, Jan Peters, and Gerhard Neumann. Model-based relative entropy stochastic search. In Advances in Neural Information Processing Systems, pages 153 154, 2015.
Joni Pajarinen, Hong Linh Thai, Riad Akrour, Jan Peters, and Gerhard Neumann. Compatible natural gradient policy search. Machine Learning, 108:1443 1466, 2019.
Jan Peters and Stefan Schaal. Natural actor-critic. Neurocomputing, (17-99)1180-11199 2008
Richard S Sutton, David McAllester, Satinder Singh, and Yishay Mansour. Policy gradient methods for reinforcement learning with function approximation. In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999. URL https://proceed:: .df..
Mohammad Khan and Wu Lin. Conjugate-computation variational inference: Converting variational infer- ence in non-conjugate models inferences in conjugate models. In Artificial Intelligence and Statistics, pages 878 887. PMIRR, 2017.
14

---

Published in Transactions on Machine Learning Research (07/2023)
Mohammad Khan, Didrik Nielsen, Voot Tangkaratt, Wu Lin, Yarin Gal, and Akash Srivastava. Fast and scalable Bayesian deep learning by reight-perturbation in Adam. In Jennifer Dy and Andreas Krause, editors, International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2611-2620. PMLR, 10-15 Jul 2018.
Manfred Opper and Cédric Archambeau. The variational gaussian approximation revisited. Newral compu- tation, (1()7778 792, 2009.
M. Abadi, M. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, T. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large-scale machine learning heteroge- ncous systems, 2015. URL https:/ https://www..r www.teensrffo.... org/. Software available from tensorflow.org.
M. Lichman. UCI machine learning repository, 2013. URL
wwt.:cic.......... uci. edu/ed
A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample test. Journal of Machine Learning Research (JMLR), 13:723 773, March 2012. ISSN 1532-4435.
O. Stasse, T. Flayols, R. Budhiraja, K. Giraud-Esclasse, J. Carpentier, J. Mirabel, A.
Del Prete, P. Souères, N. Mansard, F. Lamiraux, J.-P. Laumond, T.. Marchionni, H. Tome, and F. Ferro. Talos: A new humanoid research platform targeted for industrial applications. In 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids). pages 689 695, 2017. doi: 10.1109/HUU HHHHAAUMAISSII 82469947.
Andrew C. Miller, Nicholas J. Foti, A. D'Amour, and Ryan P. Adams. Variational boosting: Iteratively refining posterior approximations. In International Conference on Machine Learning, 2017.
Fangjian Guo, Xiangyu Wang, Kai Fan, Tamara Broderick, and David B. Dunson. Boosting variational Astat.ML/, 2016.
inference. arXiv:1611.0555000
Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(4):1303 1347, 2013. URL http://jmlr.. wrggpapprs/vvlle hoffman 13a. html.
Lucas Theis and Matthew Hoffman. A trust-region method for stochastic variational inference with appli- cations to streaming data. In International Conference on Machine Learning, pages 2503 2511. PMLR, 2015.
Jeffrey Regier, Michael Jordan, and Jon McAuliffe. Fast black-box variational inference through stochastic trust-region optimizaatinn In Advances in Neural Information Processing Systems, pages 2399 2408, 2017.
Mohammad Emtiyaz Khan, Pierre Baqué, François Fleuret, and Pascal Fua. Kullback-leible proximal variational inference. In Advances in Neural Information Processing Systems, 2015.
Mohammad Emtiyaz Khan, Reza Babanczhad, Wu Lin, Mark Schmidt, and Masashi Sugiyama. Faster stochastic variational inference using proximal-gradient methods with general divergence functions. In Conference Uncerteinty in Artificial Intelligence, pages 319-328, 2016.
Hugh Salimbeni, Stefanos Eleftheriadis, and James Hensman. Natural gradients in practice: Non-conjugate variational inference in gaussian process models. In International Conference on Artificial Intelligence and Statistics, pages 689 697. PMLR, 2018
Kamélia Daudel, Randal Douc, and François Roueff. Monotonic alpha-divergence minimisation for variational inference. Journal of Machine Learning Research, (4(2)):1-76, 2023.
Warren Morningstar, Sharad Vikram, Cusuh Ham, Andrew Gallagher, and Joshua Dillon. Automatic differentiation variational inference with mixtures. In Arindam Banerjee and Kenji Fukumizu, edi- tors, Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 3250-3258. PMLR, 13-15 Apr 2021. URL html.
attps://proceeet mlr. www..reeae
15

---

Published in Transactions on Machine Learning Research (07/2023)
S. J. Gershman, M. D. Hoffman, and D. M. Blei. Nonparametric variational inference. In International Conference on Machine Learning); 235-242 2012.
Wu Lin, Frank Nielsen, Khan Mohammad Emtiyaz, and Mark Schmidt. Tractable structured natural- gradient descent using local parameterizations. In Marina Meila and Tong Zhang, editors, International Conference Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 6680- .hmml.
6691. PMLR, 18 24 Jul 2021. URL httpss//ppl pttps://proceedingg.. mll wtt.eeediinn.......res//////1999
Philipp Bocker, Oleg Arenz, and Gerhard Neumann. Expectod information maximization: Using the i- projection for mixture density estimation. In International Conference on Learning Represenlations, 2019.
Geolfrey Roeder, Yuhuai Wu, and David K Duvenaud. Sticking the landing: Simple, lower-variance gra- dient estimators for variational inference. In Guyon, U. Von Luxburg, S. Bengio, Н. Wallach, R. Fer- gus, Vishwanathan, and R. Garnett, oditors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL 0F Paper pdf.
R.H. Byrd, Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. M Journal on Scientific Computing, 16(5):1190-1208, 1995.
16

---

Published in Transactions on Machine Learning Research (07/ 2203))
A Limitations
The scope of this work is narrow, focusing on two specific approaches for natural-gradient GMM-based variational inference. There are of course many other models that can be applied for variational inference, and, depending on the problem setting, some of these models are highly preferable over GMMs, for example, normalizing flows should likely be preferred for high-dimensional problem settings, such as (deep) Bayesian neural networks. However, for this work we assume that we indeed want to optimize a Gaussian mixture model, for example, because we require an interpretable model with smooth gradients (Ewerton el. al., 2020). Even in the field of GMM-based variational inference, alternative methods, based on boosting (Miller et al., 2017; Guo et al., 2016) or the reparamcterization trick are possible. By not using natural gradients, these methods can be applied more straightforwardly to sparse covariance parameterizations which can be beneficial for higher-dimensional problems. However, in the considered problem setting where we can learn GMMs with full covariance matrices, these methods are not competitive (Arenz et al., 2020) to the natural gradient based methods described in this work. Please refer to Appendix for a discussion of related work.
Regarding our empirical study, we want to point out several aspects that could lead to misinterpretations of the results. While we also report a secondary metric for each test problem in our main experiment (Table 8), please recall that the hyperparameter optimization was performed only with respect to the ELBO. We noticed on the planar robot experiment, that SAMTRUX could still achieve competitive ELBO when initializing the mixture model with fewer components than we used for our evaluation, while resulting in a significantly worse MMD. Hence, one should consider that for each method, the performance with respect to the secondary metric could potentially be better, if the hyperparameters were chosen correspondingly. Similarly, one should be careful when comparing the learning curves (which can be found in Appendix L), with respect to efficiency or stability, the hyperparameters have only been chosen based on the final ELBO of each run. We used this simple and hard criterion for selecting the hyperparametens, to make the experimental study more transparent and objective by removing human influence.
Human influence, of course, could not be completely avoided: For our main experiment, we only allowed for a coarse hyperparameter search over carefully chosen parameters. In contrast to the first two groups of experiments, where we used extensive Bayesian optimization to identify the best performance that we can expect for each design choice, in the main experiment we took into account that the computational budget for hyperparameter search (also ours) is limited. As a consequence, the results of our main experiment, are to some extend affected by our ability to propose good hyperparamcter ranges for the different design choices. However, we had good prior knowledge about suitable parameter ranges based on the previous groups of experiments and also based on the parameters reported by prior work, and furthermore the results of our main experiment are consistent with the performances we observed in the previous experiments (which used much less subjective Bayesian optimization), and, hence, we conclude that the effect of the human factor is rather small. To increase the transparency of our empirical study we release separate reproducibility package", which contains scripts for running each experiment, thereby documenting the exact conditions under which all our experiments have been started (including hyperparameter search).
B Potential for Negative Societal Impact
Machine learning methods can have significant impact on our daily lives. They can have positive impact on society by taking work off our hands, addressing challenges like the climate crisis, or helping to develop better medical treatments. But they can also have a negative impact on socicty by reinforcing projudices, discriminating against people, invading our privacy, wasting enormous amounts of energy, or exacerbating the imbalance of power and wealth. They can cause serious harm if we overestimate their capabilities, and they can be used maliciously, for example, to falsify data or carry out cyberattacks.
In terms of negative impact, this work is unlikely to contribute significantly to energy waste because we focus on structured representations with few parameters that can be learned efficiently. Furthermore, our work does not focus on processing or forgery of images, text, or speech, so it is unlikely to have a significant
negative impact on society due to privacy violations or misinformation. We are confident that we can use
htttp:://gtthub.ccm/lleeeeeee// retporddiiirddiilllit...
17

---

Published in Transactions on Machine Learning Research (07/2023)
our findings to help people in their daily lives and have net positive impact on society, although we of course can not foresee how our findings will be used by future work.
C Related Work
In addition to the aforementioned methods, several other methods have explored natural gradients and trust regions for variational inference. For mean-field approximations, Hoffman al. (2013) propose an efficient method for estimating the natural gradient from mini-batches, and Theis and Hoffman (2015) proposed related KL-constrained trust-region method. Regier et al. (2017) proposed sccond-order method for Gaussian variational inference based on Euclidean trust regions. Khan el al. (2015) introduced KI-based proximal point method for conjugate models and related it to natural gradient descent. For non-conjugate modols, they proposed local lincarizations. Khan et al. (2016) extended this approach to other divergences and stochastic gradients. Salimbeni et al. (2018) compute the natural gradient based on the Jacobian of the parameters of the Gaussian and its expectation parameters. The Jacobian can be computed using forward-mode differentiation, or by using inference, Daudel et al. (2023) explored the optimization of alternate divergences by presenting method for minimizing alpha divergences. Alpha divergences are a family of divergences parameterized by a scalar parameter that for a 1 also includes the KL divergence, which we consider in this work. However, the procedure proposed by Daudel et al. (2023) assumes 0 < a < 1 and is thus not applicable in our setting, Morningstar et al. (2021) recently proposed a method for amortized variational inference with GMMss. In amortized VT, the parameters of the GMM are not optimized directly. Instead neural network trained predict the GMM parameters based on given input. NPVI (Gershman et al., 2012) was presented noпparametric variational inference. However, can also be regarded as method for optimizing a very restricted type of GMM, namely, mixtures with uniform weights and Gaussian components with isotropic covariance matrices. Miller et al. (2017) and Guo et al. (2016) proposed boosting methods, which, however, can require unreasonably large mixture models, because components that have been added at previous iterations do not longer receive updates. Recently, Lin et al. (2021) investigated natural gradient GMM- based variational inference with a structured covariance matrix, which may be important for scaling these methods to higher dimensions.
movers-mmode differentiation twice. For GMM-based variational
D Background Material
We will now provide a brief introduction to natural gradient descent, Fisher information, and exponential family distributions.
D.1 Natural Gradient Descent
A possible way to update the parameters ( al iteration minimize an objective ((0) the (vanilla) gradient descent update,
g(++1) = g( (())).
where 3 is the stepsize. However, this update is not invariant to the chosen parameterization: Using a different parameterization (.)) J(O()), performing the gradient descent update based on VaJa (X) and changing the obtained parameters (++1) back to g(++1) may in general result an update g + ((++)
that is not along the gradient descent in the original paramoterization. To understand the reason for this dependency on the parameterization it helps to investigate the optimization problem that is solved by gradient descent. Namely- as be easily verified by the Tagrangian method- the gradient descent update maximizes first-order Taylor approximation of the objective subject to Euclidean trust region constraint,
arg min J(((1) + (Vo](0) ) ) 80 60
(10)
8.t.
T0 00 =e,
(11)
18

---

Published in Transactions on Machine Learning Research (07/2023)
where the stepsize B (2x) depends on the Lagrangian multiplier À and, therefore the trust region 6. However, measuring the distance of the update step using the Euclidean distance of the parameters is in general not meaningful as it depends on the chosen parameterization.
Instead, natural gradient descent (Amari, 1998) uses more general Riemannian distances d(60) 50 G(00))50, where G( is a matrix that may smoothly change as a function of 0, which also called the Riemannian metric tensor. The Riemannian metric tensor defines Riemannian manifold, a curved but locally flat manifold. Minimizing the first-order Taylor approximation subject to trust region constraint the Riemannian distance d(o0) produces update along the natural gradient
VeJ(O) (= -(0 (=
where the stepsize again dependent on the trust region. Hence, the natural gradient can be computed by preconditioning the vanilla gradient with the inverse Riemannian metric tensor.
When optimizing objective with respect to a probability distribution, scent based Riemannian distance that approximates the Kullback-Leibler divergence KL(pgc>] [L(ppco)|Ppocoo) (Poce)) between the old and the updated distribution. Such distance metric is typically much more appropriate than the Euclidean distance of the parameters, as it is invariant with respect to the chosen parameterization. The corresponding Riemannian metric tensor corresponds to the Fisher information matrix, which we will discuss in the next subsection.
use natural gradient de-
D.2 Fisher Information
The Fisher information matrix F(0) can be derived as the Riemannian distance metric that obtained from the second-order Taylor approximation of the KT. divergence at. 50 O,
pe (x)
[ vo(x) logpoxxxx
~ (x) log
dx -
2 500
po Voo log Pa (x)xx 60
po (x)
1 -000 Vepo(x)dx
/ po (x) Voo (og 80
-
T
po (x) Ven
80
x
[ (po(x)(Volog) (Vologro(x)) (ve Voloppo()))
_0
dx 50.
F())
Hence, the steepest descent update subject to trust-region constraint based the approximated KL divergence is given by the natural gradient
VoJ(O)
-Vo (=0()
We will now revisit exponential family distributions, family of distributions for which the natural gradient can be computed without explicitly computing the Fisher information matrix, mentioned in Section 3.2.
D.3 Exponential Family Distributions
Exponential family distributions distributions that can be written in the form
p(x; 0) 1i(x) exp (n) T(x) A(0)).
19

---

Published in Transactions on Machine Learning Research (07/2023)
with base measure h(x), natural parameters n/0), sufficient statistics T(x) and log partition function A(0). Furthermore, exponential family distributions can be also specified using the expectation parameters, that is, the expected value of the sufficient statistics m(0) = (=xx(x)T(x)dxx
Exponential family distributions include wide range of distributions, such as k-dimensional Gaussian distributions, where
h(x) = (2m) 4, (0))=
Σ ーー1
T(x) =
x A(O) = 2u'²0 μ + log Σ
xx T
and calegorical distributions, where
log PL PK
X x = 2
K-1 A(O) log (1
h(x) 1, m70)= log PK
T(x)
Σ P.):
(12)
[x K
and where T(x) is the one-hot encoding of the outcome x. For minimal exponential family distributions (where the natural parameters and sufficient statistics are linearly independent) the natural gradient with respect the natural paramoters coincides with the gradient with respect to the expectation paramcters, that is
F(0) 'V ((()) V mm(o)J.(m(O)
as shown by Khan and Nielsen (2018, c.f. Theorem 1).
E Proof of Theorem 2.1
The E-step of VIPS performed by sctting the auxiliary distribution to the current model. Using 0(*) denote the model parameters at iteration i, the auxiliary distribution at iteration i is given by
q(o)x) = gocn (o|x)
by construction. Following Becker et al. (2019) we express the auxiliary distribution in terms of ((o)
Got (o), q(x) Got (x), and q(x(o) gaon (x(o) using q(o|x) [o))(()))) objective as,
reformulate the lower bound
((x)
J(q. Egolo)] [/ (0(x)o) log p(x) ((x)
dx (L(go(X(0) - (L(g (())),
(13)
( (0(())))
log p(x)
# (o()))
dx J(0).
(14)
go (x)
We can see that the lower bound objective (Arenz et al., 2018; Becker et al., 2019) differs from the original ELBO objective by two additional KL divergences, and further, by computing the entropy based on the fixed auxiliary distribution q(x) Goin (x) rather than the model (x(x). Directly after the E-Stop, the lower bound is tight, J[goe,, 0(0) J(O))), because the auxiliary distributions equal to the current model. For proving Theorem 2.1 we need to prove that also the gradients match, that is,
VoJ(g 0) 0
VoJ(O) 0=00)
(15)
which would also imply that the natural gradients match, since the Fisher Information Matrix only depends on the current distribution go, not on the objective.
For proving the equivalence of the gradients, the following lemma will become handy:
Lemma E. The gradient of the Kullback- Leibler divergence between two equal distributions go (x) = π(x) is zero, that is,
4 (x) = =((x) = VoKL(goll)) 0=0-
20

---

Published in Transactions on Machine Learning Research (07 (022))
Proof.
VøKL(qellū) = Vo g0(x) log g0(x) dx = Vogo(x) loo g0(x) [(x)
dx + / Vogo(x)dx g(x)
x
Voqa(x): 0 dx + Vol -0
00 -=
As the gradient of the KL divergence between two equal distributions is zero, these terms do not affect the gradient of the objective function directly after the E-step. Hence, only need to show that
[/ go (x) log (x)
[[ go(x)) log ga (x)dx
V o
= Vo
0=00
A=
This equivalence be verified as follows:
V / (x) log (x) = / Vogo(x) log (x))xxx / Vogo(x)dx
(xyo(x) lx + Vol 0=0 / Voyo(x) log goco aoco (x)dx Vo
[/ (0(x) log yoro (x )x
Pretending for a moment that we are optimizing both objectives using the reparameterization trick, it is inspiring to relate the lower bound, the common practice of not backpropagating through the score function, yielding unbiased, but often lower-variance estimate of the ELBO gradient (Roeder et al., 2017).
F
Equivalence of the Weight Updates
We will now compare the weight update used by Lin et al. (2019a; 2020), with the weight updated used by Arenz al. (2018; 2020), and show their equivalence.
F.1 Weight Update of Lin et al. (2019a; 2020)
Lin et al. (2020) formulate the weight update in the natural paramoter space (sce Eq. 12),
= nold + 35
(16)
where each dimension 8, of the natural gradient is given by
8, (xlo) [logT(x) - log (x))
(x K) ((x) log 9o (x)].
Defining
R(()) = Ex,() (log - log (x)),
the i-th index of the natural parameters given by
n: Mold,t + (/(((() Bf((K)
21

---

Published in Transactions on Machine Learning Research (07 (022))
The new weight of component < K is given by
exp Nold, + +R((0) BR(K))
exp (7o)
((o;
Ei=1 exp1 m7
+1
k Σ
exp (Podd, + (B((k) - (B((K)
+
g(o)ald
exp (BR(O))
(())old exp (BR(0)
g(K) Jold
(k)ool exp (sf((t) (xp
(K))old exp (BR(k))
Hence, for all components (including K), we have
(0) x a(o)ata exp (B((O))
(17)
F.2 Weight Update (VIPS)
We start be expressing the component's reward R(o) in terms of R((), namely,
R(o)= = EE EE (x) o)
log + log [(/ + H(g(x)o))
E (x
o)
log + log ((0x) log ((x
(x)o)
(x)
log)(x) + (()) log ((z) R(()) + log ((0).
where we exploited that the auxiliary distribution is chosen according to the true responsibilities, ((0/2) g(olæ)). Expressing the weight update of VIPS (Arenz et al., 2018, Eq.8) in terms of M())
(o) x (ola(() 14 exp (R(o)) Tr
(ola(o)
2hc I+ow
exp (A(o) + log gota(o) 1+7
1 + Tlw A(O)
= (ola(o) exp
we can see that it exactly matches the update in Equation 17 for $ = I+m
G Comparisons with Reference Implementations
The target distributions BreastCancer, GermanCredit,, GMM20 and PlanarRobot4 were taken from VIPS (Arenz et al., 2020). The design choices of VIPS correspond to the codename ZAMTRUX our implementation. Table compares the final (negated) ELBO that we achieved our main experiments (cf. Table 5) with the performance of VIPS reported by Arenz al. (2020), that was obtained using their (C++) implerennattion. For all environments that have been tested in both works, the final ELBO performances published in this work are better than the results published in the original work, except for CermanCredit were, we could not measure any difference between both implementations.
We also compared our implementation with the Matlab implementation of Lin el al. (2020) on the target distributions that we took from their work (STM20 and STM300). However, our problem setting slightly differs from the original setting, as Lin al. (2020) use an expensive Hessian-based pre-training for initializing the GMM and compare the methods only with respect to their fine-tuning performance. We found that such pre-training is in general not necessary with implementation and directly use the respective algorithms
22

---

Published in Transactions on Machine Learning Research (07/2023)
Implementation
Breast Cancer
GermanCredit
GMM20
PlanarRobot4
VIPS//ZAMTRUX (theirs)
78.20 + 0.04
585.10 0.00
0.01 = 0.00
12.02 + 0.13
VIPS/ZAMTRUX (ours)
78.14 = 0.01
585.10 + 0.00
-0.00 = 0.00
1
Table 6: We compare the final (negated) ELBO achieved by both implementations. When using the same design choices as VIPS (Arenz et al. 2020), our implementation and hyperparameters led to better approx- imations their target distributions.
starting from the original initialization of Lin et al. (2020). Figure 1, which compares the learned model and the larget distributions based on the marginals on the STM20 experiment, demonstrates that even without pre-training, we can learn higher-quality approximations with our implementation (cf. Lin et al., 2020, Fig.3). Here, we used the SAMTRON design choices, which performed best in our experiments.
We ran the STM20 experiments on the reference implementation (with disabled pre-training) using our hyperparameters, and using the hyperparameters of Lin et al. (2020) and compare the learning curves with our SEPYFUX evaluation in Figure 2a. Our hyperparameters achieve better final ELBO even on the original implementation. Using the same hyperparameters, the learning curves of both implementations do not differ significantly, but our implementation performed slightly better. The STM300 experiment was not evaluated by Lin et al. (2020) in the first-order setting, but only when using Hessian information (using Eq. 6 left), which is often computationally prohibitive. Hence, we can only compare both implementations using our hyperparameters. The respective learning curves are shown in Figure 2b and very similar for both implementations.
Hyperparameters
We list the hyperparameters for each design choice in Table 7. Please refer to Appendix for a description of the different hyperparameters. The tested and eventually chosen hyperparameters for each experiment can be found in the reproducibility package; please refer to its README. file for links to all relevant config files.
Notes On Our Implementations
We implemented two different options for estimating the natural gradients for the component update in two separate classes. The MoreNgEstimmaoor uses weighted least squares to estimate the natural gradient using compatible function approximation; the SteinngEEstimator makes use of gradient information, using Stein's Lemma to estimate the natural gradient. For both options, a boolean hyperparameter is available to select whether self-normalizzed importance weighting, or standard importance weighting should be used to make use of samples from different distributions. Furthermore, a boolean hyperparameter can be used to disable importance sampling, only using samples from the respective component during the component update, which we enabled for all methods on WINE to reduce memory footprint. When using MORE, an additional hyperparameter can be used to select the initial ez-reuularizaton for linear least-squares; the regularizer is automatically adapted as described by Arenz et al. (2020).
We implemented three options for performing the natural gradient update of the components based on the estimated natural gradients and given stepsizes (or trust regions). The PiretNigaesCCompmeeette directly applies the natural gradient update given by Equation 4. If an updated component is no longer positive definite, the update is undone, hoping that the update would succeed in the next iteration (poten- tially with smaller stepsize). The ColtseddCoomoonnt(pppctttt the improved Bayesian learning rule (Eq. 8). We noticed that the update may still result in non-invertible covariance matrices (albeit much less frequently compared to the natural gradient update) due to numerical errors, in which case we also undo the respective updates. The solves an optimization problem to find stepsizes that result positive definite covariance matrices and updated distributions that respect the desired bound on the KL divergence, as discussed in Section 3.3. For solving the convex optimization problem, Arenz et al. (2018) used L-BFGS-B (Byrd et al., 1995)
updates the components according to
23

---

Published in Transactions on Machine Learning Research (07/2023)
0.3
025
0.15
020
0.2
0.10
0.10
0.05
0.05
0.25
0.15
0.20
0.25
0.15
0.10
0.20
0.10
0.05
0.05
0.05
0.00
0.00
0.00
0.00
0.20
0.25
0.25
0.20
0.15
0.20
0.10
0.10
0.30
0.05
0.05
0.05
0.00
0.00
0.00
0.00
0.4
0.25
0.15
0.20
0.20
0.2
0.1.0
0.35
0.15
0.20
0.10
0.1
0.05
0.05
0.05
0.20
0.20
0.15
0.2
0.30
0.10
0.1
0.1
0.0
0.0
0.00
0.00
-20
-10
100
-20
-10
100
-20
-10
10
-10
10
Figure 1: A representative plot of the 20 marginal distributions of the GMM learned with SAMTRON for the STM20 experiment is shown in red. The marginals of the Mixture of Student-T target distribution are shown in blue and hardly distinguishable.
24

---

Published in Transactions on Machine Learning Research (07 (07/202) 2023)
10a
reference
reference (our params)
reference (ref. params)
ours
ours (our params)
10a
10
0
1
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1e6
1e6
(a) STM20
(b) STM300
Figure 2: The learning curves plot the ELBO (in logarithmic scale) the number of samples for implementation and the reference implementation, where both use the SEPYFUX design choices. On STM20 our hyperparameters lead to better approximation than the hyperparameters used by Lin et al. (2020), and when using the same hyperparameters, our implementation performs slightly better. On STM300 both implementations perform similarly using our hyperparameters. Reference hyperparameters are not available for STM300.
Module
Design Choice
Hyperparameters
Component Adaptation
Fixed
Adaptive
adding iter.
deletion iter.
# prior samples
# DB samples
Component Stepsize
Fixed
initial stepsize
Decaying
initial stepsize
annealing exponent min stepsize
Adaptation
Adaptive
initial stepsize
max stepsize
NgBased Component
Direct iBLR Trust-Region
Updater
Stein
self-normalized IW self-normalized IW
samples use own samples
NgEstimator
MORE
lg fg-coefficient
Sample Selector
Comp.-Based
desired samples desired samples
reused samples reused samples
Mixture-Based
Weight Stepsize
Fixed
Suttial_stepsiie initial stepsize initial stepsize
Decaying
annealing exponent min stepsize
Adaptation
Adaptive
max stepsize
Weight Updater
Direct Trust-Region
self-normalized IW
self-normalized IW
Table 7: The table lists each hyperparameter that was tuned at any of the experiments.
25

---

Published in Transactions on Machine Learning Research (07/2023)
optimizer. However, it was difficult to efficiently integrate such optimizer into our Tensorflow (Abadi et al., 2015) implementation, and furthermore, the optimizer would sometimes fail for numerical reasons. Trstead, we implemented a simple bisection method, by iteratively refining an initially sufficiently large bracket on log,3, by evaluating the center of the bracket and using it as new upper value when the corresponding KL divergence is loo small or the covariance matrix not positive definite, or; otherwise, as new lower value. Although naive, this approach numerically robust by only requiring the sign of the gradient of the KL. divergence, converges to the optimum due to the convexity of the problem and can be efficiently compiled into compute graph.
For the SampleSelector module, we implemented two options. Both options can make use of samples from previous iterations (as discussed by Arenz et al. (2020)), by setting the hyperparameter Vreused larger than zero. The LinSempleSelector computes the effective sample size Nelf on the GMM and draws max(0, N ndee nett) new samples from the GMM. For nreused = 0 this approach corresponds to the procedure used by 1Bayes-GMM (Lin et al., 2020). The VipsSampleSelector computes (e)(o) for each component and draws max (0, fldes (etr(o)) new samples from each component, matching the procedure used by VIPS (Arenz et al., 2020).
We implemented two options for performing the weight update. The Direct Upiget directly updates the categorical distribution using Eq. 9, whereas the TrustRegion WeightUpplater uses our bracketing search to stay within a given trust region. Whether standard importance sampling or self- normalzeedd impor- fance sampling should be used for estimating the component rewards R(o) can be chosen with a hyperpa- rameter, that is available for both options and can be chosen independently to the respective hyperparameter of the NgEstimator.
For stepsize adaptation, we implemented three options for both, the ComponenllStepsiixxAdlilioo module and Fixed WeightStep-
and the WeightStepiizeAaaptattoi module. The Adapttatnni simply return a fixed stepsize, chosen as a hyperparameter. The DempyinnCooponentttler sizettaatttion and Decaying WeeightSeepsizeAdappee
return an exponentially decaying stepsize based on the number of times the respective distributton has been updated, as described by Khan et al. (2018). The and ImprovementBased AttgghtSSepsteeedptttttt increase the stepsize if the last update of the respective distribution improved its reward, and decrease it otherwise, as proposed by Arenz et al. (2020).
For the ComponentAdaptationn nothing, keeping the number of components fixed throughout optimization; the uses the procedure of VIPS (Arenz et al., 2020) to delete bad components, and to add new components in promising regions. Arenz et al. (2020) only considered samples from the sample database for initializing the mean of the new component, however, we had to disable the sample database for the high-dimensional prob- lems STM300 and WINE to reduce memory footprint. Hence, we introduced an additional hyperparameter to specify additional samples from the prior distribution (the same distribution that was used for drawing the mcans of the initial GMM), that should be drawn specifically for ComponentAdaptatiion
module, we implemented two options. The FizeeCCopponentAAdppttie does
J Test Problems
We performed several experiments on target distributions taken from related work. For additional details, please refer to the corresponding work. The BreastCancer, GermanCredill, PlanarRobol and GMM exper- iments were used by Arenz et al. (2020) and the Student-T experiments were used by Lin et al. (2020). For the most detailed (and fully accurate) specification of all target distributions, please refer to our code supplements httt:::///iituu..co//// wtt.:/ooo////////e//// ggmvi).
In the BreastCancer and GermanCredit experiments, we aim to approximate the postcrior distribu- lion of a logistic regression problem. The dimensions are 25 and 31, respectively. The data sets can be obtained from the UCI Machine Learning Repository (Lichman, 2013) and contain 1000 and 569 data points. Mimicking the setup of Arenz et al. (2020), the ELBO computed on the full data set during training and evaluation. For the minibatch variants GermanCreditMM2 and BreastCancerMB we also use batches from the full data set, such that the respective ELBOs are comparable with the
26

---

Published in Transactions on Machine Learning Research (07/2023)
original experiments. However, we report the full-batch ELBOs to remove unnecessary noise in the evaluation.
In the WINE experiment want to approximate the posterior distribution over the weights of neural network that predics the scalar wine quality based on cleven features using the WINE data set Lichman (2013). The network has two hidden layers of width 8, resulting in 177 parameters. The likelihood is given by the root mean squared error over a minibatch of size 128. We split the data set into a training and test set, and make sure that the training and test sets are deterministic given the seed. As we are prinarily interested in the ability of the different methods to approximate given target distribution (rather than testing how the approximations perform on the downstream task), report the training set ELBOs. However, Table 8 reports as secondary metric the mean squared error of the prediction using approximate Bayesian inference on the test set.
In the planar robot experiment, we aim to sample joint configurations of 10-link planar robol (all links have the same length) and aim to reach one of four goal positions. The target distribution is Gaussian in the endeffector configuration space (but non-Gaussian in the joint configuration space). A meao Gaussian prior on the joint angles additionally used to prevent non-smooth configuration.
In the GMM experiment we aim to approximate an unknown Gaussian mixture model with 10 components and a varying number of dimensions. Arenz et al. (2020) only considered 60 dimensions, but we increased the dimensionality to up to 100. The mean is sampled uniformly in the range [-50, 50) and the covariance matrices Σ = A A I are created by randomly sampling the elements of the square matrix A.
The mixture of Student-T experiment (STM) is similar to the GMM experiment but uses Student-TT components instead of Gaussians. We exactly follow Lin et al. (2020) by considering 20-dimensional mixture with 10 components with mean uniformly sampled in [-20, 20), and 300-dimensional mix- lure with 20 components sampled in [-25, 25]. We follow Lin et al. (2020) by initializing the components of the GMM by sampling the mean from a zero-mean Gaussian distribution with di- agonal covariance with standard deviation 100, and by initializing the diagonal covariance matrices with Emit = 300I. However, instead of pre-training with a second-order method, we directly start training from the initial GMM.
The TALOS experiment is based on the implementation by Pignat et al. (2020). The poses of both feet, as well as the positions of the left end-effector and the center-of-mass are computed for the given joint positions based on kinematic model of the robot. The target positions of the feet are given by [-0.02, 0.09, -0.] and |-0.02, -0.09, -0.] and their target orientation are given by identity rotation matrices R = Ig- The likelihood for cach foot, is given by a 12-dimensional Gaussian distribution that penalizes deviations from these goal-parameters using a standard deviation of 0.2 for the Cartesian positions and a standard deviation of 0.1 for each entry of the rotation matrix. The likelihood of the left endeffector position given by a Gaussian with mean [0.1, 1.] and diagonal standard deviations of 0.02. Violations of the inequality constraints that the joint angles should be within their limits, and the center-of-mass within the convex polygon spanned by the feet are penalized using the log-density of a Gaussian distribution placed on the violated bound (Pignat et al., 2020) using standard deviation of 0.01 for the center-of-mass and 0.05 for the joint limits.
K Full Table for Experiment 3
The complete table for Experiment 3, showing all tested candidates, can be found in Table 8.
L Learning Curves
The learning curves (ELBO over time) for our main experiments are shown in Figure 3.
27

---

Published in Transactions on Machine Learning Research (07 /023)
+7.8e1
+7.8e1
+5.8509e2
7.00
21.0
0.10
5.25
17.5
0.08
14.0
0.06
3.50
10.5
0.04
7.0
1.75
3.5
0.02
0.00
0.0
0.00
1000
2000
3000
1000
2000
3000
2000
4000
6000
(a) Breast Cancer
(b) Breast Cancer (MB)
(c) German Credit
+5.8509e2
1.00
1.00
1.20
0.75
0.75
0.96
0.72
0.50
0.50
0.48
0.25
0.25
0.24
0.00
0.00
0.00
2000
4000
6000
500
1000
1500
0
500
1000
1500
(d) German Credit (MB)
(e) GMM (20D)
(I) GMM (100D)
1.00
+1.4e1
+1.1e1
15.0
15.0
0.75
12.5
12.5
10.0
10.0
0.50
7.5
5.0
5.0
0.25
2.5
2.5
0.00
0.0
0
2000
4000
6000
20000
40000
60000
0
10000
20000
(g) STM (20D)
(h) STM (300D)
(1) Planar Robot (4 goals)
+1.35e3
-2.5e1
200
1.80
zamtrux
sepyrux
1.35
samyrux
samyron
100
0.90
samtron
samyrox
sepyfux samtrux
samtrox
50
0.45
0.00
0
20000
40000
60000
25000
50000
75000
(i) WINE
(k) TALOS
Figure 3: The learning curves for main experiment (Experiment 3) show the negated ELBO over time (in seconds). Shaded arcas show best and worst performance. Note that hyperparameters were selected with respect to final ELBO.
28

---

Published in Transactions on Machine Learning Research (07/2023)
SAMTRUX
SAMTROX
SAMTRON
SAMYRUX
SAMYROX
SAMYRON
SEFYFUX
SEPYRUX
ZAMTRUX
Experiment
Metric
-LBO
78.01
78.05
78.00
78.37
78.63
78.61
79.78
79.91
78.14
00.0
#0.02
10.02
+0.07
+0.07
+0.05
+0.40
+0.93
#0.01
BreastCancer
MMD
1.1e 3
1.1e 3
9.9c
1.2p
1.1e
1.3e 3
2.3e 3
2.8e
1.1e 3
:le
11
18c
Ele
12c
Lle
1:3e
13e
#1e
ELBO
79.66
79.71
78.41
79.68
79.94
79.96
80.13
79.38
83.65
BreastCancer
+0.32
+0.43
+0.03
10.19
0.223
+0.26
+0.69
25.25
+0.97
(minibatches)
MMD
2.9e
2.Ge
1.7e
2.2e
3
2.8e
2.9e
3
2.5e
1.9e
8.3e
#1e
+e
8
=9e
Se
tdl
ELBO
585
585.10
585.11
585.12 12
585.11
585.10
0.000
+0.00
#0.00
+0.00
0.000
+:0.0
+0.01
+0001
+0000
GermanCredit
MMD
5.9e
5.5e
5.7e
5.5e
5.7e
5.6e
6.0e
4 5
66
+5e
3: 5
:4ee
3:e 5
+:5e
+5e
77e 5
ELBO
585 12
585.19
585. 12
585.13
585. 13
585 .2
585. 685.12
585.12
585.35
GermanCredit (minibatches)
00.01
#0.01
10.00
+0.00
10.01
10.00
10.00
10.00
10.03
MMD
6.5e 4
8.1e 4
5.9c
6.0e
4
6.2c
4
5.4e 4
5.8e 4
5.7e 4
1.9e 3
#1
4
#:
#4e
5
=5e-
5
=7e
#Ge 5
#e
4
= 4
ELBO
11.47
11.47
11.47
13.16
12.98
12.93
17.26
16.35
11.48
0000
+0.04
#0.04
+0.20
0.00
+30.3
#:2.13
+0.65
#0.04
Planar Robot
MMD
1.6c 2
1.6c
2
1.6e 2
3.7e 2
2.6e
2.4e 2
4.5e
4.4e
1.6c 2
8e
2
#1e
3
7
2
# 3
+
#7e
= 4
ELBO
0.00
0.00
-0.00
0.00
0.01
0.00
0.43
0.28
-0.00
10.00
10.00
10.00
10.00
10.03
10.00
10.15
10.15
10.00
GMM20
Modes
10.00
10.00
10.00
10.00
9.90
10.00
7.90
8.80
10.00
10.00
10.00
10.00
10.00
10.28
10.00
11.16
11.11
10.00
ELBO
0.00
0.01
0.05
0.09
0.16
4.54
0.96
1.21
10.00
10.00
10.03
40.03
10.07
10.09
10.52
10.36
10.15
GMM100
Modes
10.00
9.90
10.00
9.70
9.40
0.00
3.00
+0.00
+0.00
+0.28
10.00
+0.43
10.46
+0.00
+0.42
-ELBO
0.12
0.02
0.00
0.04
0.05
0.03
0.46
0.46
0.54
10.07
10.04
10.00
+0.02
+0.03
10.01
+:006
+:007
10.17
STM20
Modes
8.90
9.80
10.00
10.00
9.90
10.00
9.20
9.20
6.10
10.66
10.38
10.00
10.00
10.28
10.00
10.83
10.71
10.99
-ELBO
15.00
15.24
14.96
22.41
22.28
22.50
26.69
26.87
N/A
10.38
10.36
10.48
10.23
10.13
10.15
10.39
10.45
STM3300
Modes
13.70
14.10
14.30
9.57
9.80
9.10
0.90
0.30
N/A
11.80
11.50
11.53 553
.1.73
11.58
10.99
1..89
10..3
EELBO
1435.62 15
1429.02
1423. 12
1432.72
1430.80
1433.87
3279.32
4100.92
16 503.38
+34.06
+:30.3
+:31.55
30..49
+33.66
39.449
+1425. 12
+ 12995.57
#813.75
WINE
MSE
0.48
0.48
0.47
0.48
0.47
0.48
0.67
0.76
0.80
+0.02
+0.02
+0.02
+0002
+0.02
0000
+0.14
#0.17
+0.13
-ELBO
24.32
24.30
24.43
24.13
23.91
24.00
16.64
19.00
23.69
#0.22
10.10
+0.16
#0.14
+0.13
+0.23
+5.26
+1.07
#0.16
TALOS
H(g)
16.81
16.88
16.82
-17.26
11.32
-17.25
25.03 1:5.46
22.34 +1.11
-16.91
0.00
10.07
10.10
+0.16
10.16
10.07
Table 8: The full table for our main experiment shows all tested candidates, well the secondary metrics.
29

---

Published in Transactions on Machine Learning Research (07/2023)
SAMTRUX
SAMTROX
SAMTRON
SAMYRUX
SAMYRON
SAMYRON
SEPYFUX
SEPYRUX
ZAMTRUX
Experiment
Metric
-ELBO
78.00
78.04
78.01
78.62
78.62
78.61
79.68
80.13
78.15
#0.02
+0.02
10.01
+0.08
=0.05
+0.06
+0.40
+0.79
+0002
Brenst
ACCURACY
98.72
98.79
98.75
98.79
98.70
98.70
98.77
98.75
98.70
10.11
10.12
10.09
10.09
1.008
10.11
10.13
1.009
100.1
-ELBO
79.54
79.54
78.40
79.70
79.90
79.95
79.83
79.39
85.88
Brnner
10.19
+0.24
+3.03
+0.20
00.22
+3.23
#0.70
+0.70
#3.41
(minibatches)
ACCURACY
98.73
98.79
98.66
98.72
98.77
98.79
98.80
98.75
98.89
+0.10
+0
+:008
+0.15
+0007
#0.14
+0.10
+0.05
+0.08
-ELBO
585.10
585.10
585.10
585.11
585.11
585.11
585.12
585.
585.10
+0.00
+:000
+:000
+:000
=0.00
+0000
#0.01
+:0.01
+00000
GermanCredit
ACCURACY
78.47
78.47
78.55
78.47
78.55
78.47
78.52
78.53
78.54
+0.13
+0.10
+0010
+0.10
+0.06
+0.10
+0112
+0.09
+0.11
-ELBO
585.12 12
585.14
585. .2
585.13
585.12 12
585.12
585 .2
585 11
585.32
GremtnCredit (minibatches)
10.01
+0.01
1:0.00
+:000
00.01
10.00
10.01
10.00
::0.02
ACCURACY
78.52
78 78.52
78.55
78.49
78.53
78.52
78.49
78.54
78. 78.46
+0.09
10.15
+0.10
10.09
+0.09
10.13
10.11
#0.10
10.14
Table 9: We reran the experiments on the logistic regression tasks where we also evaluatod to accuracy of the Bayesian inference predictions (based on 2000 samples from the learned model). These experiments the same hyperparameters and seeds that used for the experiments reported in Table 8. For these tasks, slightly better approximations of the posterior, did not result in significant differences in the prediction accuracy.
M Tips for the Practitioner
We do not want to make the impression that any variant, e.g. our proposed hybrid variant SAMTRON is to be unconditionally preferred. Instead, the diffcrent variants have their own advantages and disadvantages. this section, we would like discuss some tradeoffs that should be considered when applying our generalized framework to practical problem.
M.1 Exploitation vs. Exploration
We found that exploration a key concern, when aiming to learn a highly accurate approximation that covers many different modes of the target distribution. Several design choices are particular important for discovering different modes.
Sampling from the Components. Based our experience, sampling from the individual compo- nents (design choice M) is not only important for discovering different modes, but also for ensuring that the discovered modes are actually covered by components with appropriate weight. Sampling from the mixture (P) may lead to vicious circle, where components that don't approximate their respective mode very well, will have low weight and, thus, not receive sufficient samples improve the approximation and obtain larger weight. On the other hand, design choice M may end up sampling at locations that are not relevant for good approximation, which may lead to bad sample efficiency. Hence, when very fast learning is essential, P can be preferable.
Adapting the Number of Components. Adapting the mumber of components by adding new components promising regions very effective technique to improve exploration. Often, initial- izing a GMM with a single, high-entropy component and adding new components good sample locations (Arenz et al., 2020) is still able to discover multiple different modes of the target distribu-
30

---

Published in Transactions on Machine Learning Research (07/2023)
tion. However, found that initializing with many components is effective strategy to further improve exploration. As the component adaptation strategy also involves deleting components with low weights, it is even feasible to start with large initial GMMs and naturally reduce its size by deleting low weight components that no longer improve. Dynamically growing and shrinking the size of the GMM can not only help in exploration but may also improve sample efficiency, and therefore seems generally preferable compared to optimizing a fixed-size GMM.
Number of Samples. The hyperparameter of the SampleSelector that specifies the number of desired samples does not only trade off training stability and efficiency, but also affects exploration. Hence, may be beneficial to use relatively large sample sizes, even they are not needed for improving the natural gradient estimates.
M.2 Efficiency vs. Stability
The training stability another major concern. However, when using adaptive trust regions, we can often obtain efficient and stable updates with little hyperparameter tuning.
Sample Size and Learning Rate. The number of samples used for estimating the natural gradient and the step size (or trust region) of the natural gradient updale should be considered in combination, since better estimates allow for larger step sizes. We found that adaptive trust regions for the component updates (design choices and R. can ensure stable improvements for a wide range of different sample sizes, which results in good sample efficiency and stability with little hyperparameter tuning. However, too small step sizes can lead to computational overhead and too large step sizes can harm exploration. When evaluating the target distribution is cheap, recommend to use large sample sizes while upper-bounding the maximum trust-region to help exploration.
Trust Regions. While the performance of the iBLR update was often comparable to trust regions, lypically least slightly worse, and also the slep size adaplation seems work slightly better when controlling trust region rather than the learning rate. Hence, recommend trust regions for the component updates, although the iBLR update can also be worth trying. For the weight update, adaptive trust regions are preferably as they require little hyperparameter tuning, However, also directly controlling the stepsize, and even using fixed stepsize can achieve the same performance with good hyperparameters.
Natural Gradient Estimation. Exploiting first-order information by using Stein's lemma for estimating the natural gradient (design choice S) is often around one order of magnitude more efficient than use the zero-order method MORE. Furthermore, solving the weighted least-squares problem can become relatively slow for large number of dimensions. Hence, we recommend to always use Stein's lemma, when the target distribution is differentiable.
31


Reviews:
[{"review_id": "KZ0JIBzike", "reviewer": "Reviewer_qBLy", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The studied problem is important and the motivation is clear.\n- The provided new perspective appears technically sound\n- The comprehensive summary of GMM-based variational inference variants and provided codebase can be very useful for the practitioners\n\nWeaknesses:\n\n- It will be better if the authors can provide a summary of the pros and cons for those combinations and provide some suggestions about which method should be used for which scenarios. \n- The evaluation metric is mostly ELBO which sometimes may not align with the metrics that we care about in the end. It will be better to also include other metrics such as MSE or classification accuracy.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper shows the equivalence between two GMM-based variational inference methods, VIPS and iBayes-GMM. Based on this new insight, the authors further summarize the design components and design choices in GMM-based variational inference and test these combinations on Bayesian logistic regression on the Breast Cancer dataset, a small Bayesian neural network on the Wine dataset. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See Weaknesses.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "o6WkVrLsgJ", "reviewer": "Reviewer_fDiF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "STRENGTH:\n\nThis work is highly practical and has an impressive amount of empirical work. \n\nThe paper is also very well-written and has a clever way of presenting the material in such a way that the key idea & its most important execution would come across clearly even if the reader are not familiar with the prior work of VIPS and iBayes-GMM.\n\nThe key idea is also pretty intriguing to me: by forging a mapping between the high-level layout of the two approaches, one can then mix and match their specific approximation subroutines (originally devised to each component of the high-level map)\n\nWEAKNESS:\n\nDespite the strength above, I unfortunately am not quite convinced with the theoretical component of this paper, which might be both incomplete and flawed. Let me elaborate below:\n\nFirst, the established equivalence between VIPS & iBayes-GMM appears to rest of the claim that the gap between the vanilla ELBO and the relaxed version of VIPS would disappear after each expectation step in the EM scheme of VIPS. \n\nGiven that the gap disappears, the equivalence would obviously follow but the claim is assumed rather than being proven formally. This kind of makes the claimed equivalence vacuous.\n\nSecond, another key step in establishing this equivalence lies with the equation following Eq. (5) in which it is shown that the gradients of the entropy terms are equal. But, the derivation does not seem right to me. Please correct me if I miss any essential steps below:\n\nFrom the presented derivation, I get that:\n\n$\\frac{\\partial}{\\partial\\boldsymbol{\\theta}}\\Bigg[\\displaystyle\\int_{\\mathbf{x}} q_{\\boldsymbol{\\theta}}(\\mathbf{x})\\log q_{\\boldsymbol{\\theta}}(\\mathbf{x})\\mathrm{d}\\mathbf{x}\\Bigg] = \\displaystyle\\int_\\mathbf{x} \\Bigg[\\frac{\\partial q_{\\boldsymbol{\\theta}}(\\mathbf{x})}{\\partial\\boldsymbol{\\theta}}\\cdot \\log q_{\\boldsymbol{\\theta}}(\\mathbf{x})\\Bigg]\\mathrm{d}\\mathbf{x}$\n\nThen, using the same logic\n\n$\\frac{\\partial}{\\partial\\boldsymbol{\\theta}}\\Bigg[\\displaystyle\\int_{\\mathbf{x}} q_{\\boldsymbol{\\theta}}(\\mathbf{x})\\log \\tilde{q}(\\mathbf{x})\\mathrm{d}\\mathbf{x}\\Bigg] = \\displaystyle\\int_\\mathbf{x} \\Bigg[\\frac{\\partial q_{\\boldsymbol{\\theta}}(\\mathbf{x})}{\\partial\\boldsymbol{\\theta}}\\cdot \\log \\tilde{q}(\\mathbf{x})\\Bigg]\\mathrm{d}\\mathbf{x}$\n\nwhich would require $q_{\\boldsymbol{\\theta}}(\\mathbf{x}) = \\tilde{q}(\\mathbf{x})$. But it is not clear how it will be achieved.\n\n--\n\nFurther minor points:\n\nAlthough the presentation is made in a very clever way that (in my opinion) allows people to appreciate the core idea without needing exposure to previous work of VIPS and iBayes-GMM, I still think that for a journal, it is best that those prior work are covered comprehensively in separated sections to improve the educational value of this work.\n\nIt is also good to include detailed background on natural gradient, exponential family, Fisher Information Matrix and etc. A succinct review of the work of (Amari, 1998) would be nice to included (in the appendix if not in the main text)\n\nIn addition, it seems Eq. (3) might have missed $q_{\\boldsymbol{\\theta}}(o)$. Otherwise, how would it be optimized?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a unified perspective of variational inference for Gaussian Mixture Model (GMM) which combines two lines of independent work, including VIPS and iBayes-GMM.\n\nIn particular, it is shown that VIPS further relaxes the vanilla VI's ELBO of GMM with another lower-bound, which is interestingly optimized via an EM procedure. It is then claimed (by this paper) that the gap between the ELBO and its derived lower-bound in VIPS would zero out (with respect to the current estimate model parameter) after every expectation step (E-step).\n\nThus, at the beginning of each maximization step, with the above gap dismissed, it follows that per iteration, the objective of M-step in VIPS is the same with the objective in iBayes-GMM.\n\nWith this unified view, it follows naturally that the specific implementation subroutines (i.e., approximation to exact optimization or computation of intermediate quantities) used in each approach (i.e., VIPS & iBayes-GMM) can now be mixed and matched, resulting in a broader implementation scheme comprising 432 combinations of those implementation subroutines.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I request the following changes:\n\n1. Elaborate on how the gap between the ELBO & VIPS's relaxed version disappear. I believe this is key to establishing equivalence. Otherwise, if this is assumed away, the rest of the equivalence might be a bit too straightforward to make a contribution.\n\n2. Adding comprehensive coverage for VIPS, iBayes-GMM along with background on natural gradient & Fisher Information Matrix.\n\n3. Empirically, please demonstrate that the VIPS's ELBO gap is tight after each E-step as claimed.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2KqzUvaTub", "reviewer": "Reviewer_dU7d", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I appreciate the effort that the authors have put into providing a unified perspective on natural gradient variational inference with mixture models. However, I think the clarity of the work should be improved for it to be considered for publication.\n\nPlease find below a detailed review regarding the aspects that I think could be improved.\n\n### 1. Clarity\n \nI found really confusing the constant back and forth between \n- VIPS and iBayesGMM\n- existing work and contributions.\n\nFor example:\n- In Section 1: \n*“Namely, VIPS uses samples from the individual components to estimate the natural gradient using the policy search method MORE (Abdolmaleki et al., 2015), which is also applicable to non-differentiable target … outperforms both previous methods significantly\"*.\n\n- In Section 2:\n\n(a) In Section 2.2.1: Important aspects regarding iBayesGMM updates are deferred to Section 3.1.2, blurring the line between existing work and contributions. In my view Section 2.2.1 should be devoted to summarizing the key aspects of iBayesGMM, including theoretical guarantees, parameters updates and implementation details.\n\n(b) Similarly, I think that Section 2.2.2 should be devoted to summarizing the key aspects of VIPS. Please also refrain from using “us” when you are referring to previous existing work (e.g. “prevents us”).\n\n(c) Section 2.2.3 should be in a separate section to clearly indicate that it is a contribution. Please state your result in a proposition/theorem environment and write down formally the proof in a proof environment. Please also write down the proof that the gradient of the KL divergence between two equal distributions is zero. Can you explain the second inequality after “Furthermore, the gradients of the entropy terms are equal, since”? In addition, the authors should comment more on the implications of their result (theoretically and from an implementation perspective).\n\n - In Section 3: I found it difficult to differentiate between related work and the authors’ contribution here. Algorithm 1 appears to present a unified framework, yet most of Section 3 seems like related work?\nConsider using different names for each option instead of Option 1 and Option 2 everywhere (e.g. Option A1, A2, B1, B2, …)\n\nNB: It might help to have a table summarizing the differences between VIPS and iBayesGMM once these two algorithms have been introduced in lieu of the paragraph from Section 1 mentioned above (from both a theoretical and an implementation point of view).\n\n### 2. Related work\n- I think the authors should consider acknowledging recent works of Variational Inference using Mixture Models, such as:\n\n[1] Warren Morningstar, Sharad Vikram, Cusuh Ham, Andrew Gallagher, Joshua Dillon. Automatic Differentiation Variational Inference with Mixtures. Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, PMLR 130:3250-3258, 2021.\n\n[2] Kamélia Daudel, Randal Douc, and François Roueff. Monotonic alpha-divergence minimisation for variational inference. Journal of Machine Learning Research, 24(62):1–76, 2023.\n\n- I was also surprised not to see a reference to:\n\n[3] Samuel Gershman, Matt Hoffman, and David Blei. Nonparametric variational inference. In Proceedings of the 29 th International Conference on Machine Learning, Edinburgh, Scotland, UK, 2012.\n\n- *“Arguably the two most effective algorithms for GMM-based VI”*: Perhaps the authors could make a less assertive statement here, in light of new bodies of research for GMM-based VI such as [1,2].\n\n### 3. Minor comments\n- *“Their approach can be shown to converge, even when the M-Step does not consist of single natural gradient updates, however, it was not yet proven, that their proposed procedure, which does use single natural gradient steps, also performs natural gradient descent on the full mixture.”*\nConsider splitting this sentence in two.\n\n- *“two equal distribution is zero”*\nTypo: an s is missing\n\n- *“Pseudo-Code for our framework is show in Algorithm 1.”*\nTypo: shown\n\n- A parenthesis is missing at the end of eq. (4) and in the definition of R(o) after eq. (11)\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors are interested in two algorithms for GMM-based variational inference, namely VIPS and iBayes-GMM. They obtain that these algorithms have equivalent derived updates, despite having different implementations and theoretical guarantees. They then evaluate empirically the impact of different implementation choices.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Clarifying the manuscript according to the **1. Clarity** above is a critical adjustment for this work to be considered for publication. **2. Related work** would strengthen the work.", "broader_impact_concerns": "I do not see any concerns on the ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper has already undergone three revisions during the reviewing process, improving clarity and addressing reviewers' comments. My understanding is that most, if not all, of the reviewers' suggestions have been incorporated. Therefore, I don't see a need for a further revision, and I recommend acceptance as is.
Decision: Accept (poster)