Fisher meets Feynman: score-based variational inference with a product of experts
Authors: Diana Cai
Robert M. Gower
David Blei
Lawrence K. Saul
Venue: NeurIPS 2025 spotlight ()
Keywords: black-box variational inference
score matching
product of experts
multivariate t-distribution
Feynman identity
Fisher divergence
expressive variational family
Primary area: probabilistic_methods
TL;DR: We develop a score-based variational inference method that learns a product-of-t-experts model via a Feynman-identity latent-variable formulation, reducing inference to a sequence of convex quadratic programs with provable convergence.

Abstract: We introduce a highly expressive yet distinctly tractable family for black-box variational inference (BBVI).  Each member of this family is a weighted product of experts (PoE), and each weighted expert in the product is proportional to a multivariate $t$-distribution.  These products of experts can model distributions with skew, heavy tails, and multiple modes, but to use them for BBVI, we must be able to sample from their densities.  We show how to do this by reformulating these products of experts as latent variable models with auxiliary Dirichlet random variables.  These Dirichlet variables emerge from a Feynman identity, originally developed for loop integrals in quantum field theory, that expresses the product of multiple fractions (or in our case, $t$-distributions) as an integral over the simplex. We leverage this simplicial latent space to draw weighted samples from these products of experts---samples which BBVI then uses to find the PoE that best approximates a target density. Given a collection of experts, we derive an iterative procedure to optimize the exponents that determine their geometric weighting in the PoE. At each iteration, this procedure minimizes a regularized Fisher divergence to match the scores of the variational and target densities at a batch of samples drawn from the current approximation.  This minimization reduces to a convex quadratic program, and we prove under general conditions that these updates converge exponentially fast to a near-optimal weighting of experts.  We conclude by evaluating this approach on a variety of synthetic and real-world target distributions.

Fisher meets Feynman: score-based variational inference with a product of experts
Diana Cail, Robert M. Gower', David M. Blei', Lawrence K. Saul
Flatiron Institute
Columbia University
Abstract
We introduce a highly expressive yet distinctly tractable family for black-box variational inference (BBVI). Each member of this family is weighted product experts (PoE), and each weighted expert the product is proportional to a multivariate (-distribution. These products of experts can model distributions with skew, heavy tails, and multiple modes, but to use them for BBVI, we must be able to sample from their densities. We show how to do this by reformulating these products experts as latent variable models with auxiliary Dirichlet random variables. These Dirichlet variables emerge from a Feynman identity, originally developed for loop integrals in quantum field theory, that expresses the product multiple fractions (or in our case, f-distributions) as integral over the simplex. We leverage this simplicial latent space to draw weighted samples from these products of experts experts-samples approximates a target density. Given a collection of experts, we derive an iterative procedure optimize the exponents that determine their geometric weighting in the PoF. At each iteration, this procedure minimizes a regularized Fisher divergence match the scores of the variational and target densities at batch of samples drawn from the current approximation. This minimization reduces to a convex quadratic program, and prove under general conditions that these updates converge exponentially fast to near-optimal weighting of experts. We conclude by evaluating this approach on a variety of synthetic and real- world target distributions.
samples which BBVI then uses find the PoE that best
Introduction
The goal of variational inference (VT) approximate intractable probability density p by the best-matching density y from some simpler parameterized family Q [4, 23, 53J. VI is typically used when it is difficult to draw samples from p. but often there exists a "black-box" way compute the gradient of log p (that is. the score) at any point in the domain of IR [28, 42]. Each such gradient considerably more than what is provided by mere
evaluation provides wealth of sample-and and for this reason a growing number of researchers have begun to investigate score-based methods for black-box VI (BBVI) [6, 7, 36, 37]. This direction of research is also motivated by the remarkable successes of basee methods for generative modeling [18, 19, 47-49].
Despite this allure, score-based BBVI still faces many challenges. There inherent trade-offs that arise between the expressivity of the variational family Q and its ease of use. For BBVI it must be tractable to evaluate and draw samples from each g Q: it must also be tractable optimize over Q and find its best approximation to the target These trade-off's must be managed by any practitioner of BBVI. At the same time, researchers need more than methods which are merely practical. BBVI replaces an intractable problem in inference by more tractable one in optimization, but still the challenge remains to prove theoretical guarantees for the optimizations that arise this framework [3, 11, 24, 57].
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In this paper we introduce a new family for bosee BBVI that navigates these trade-offs in an appealing fashion. The densities in this family are highly expressive and yet manageably tractable, and we are also able to provide certain theoretical guarantees for the optimizations required for score- based VI. Each density in this family is a product of experts (PoE) [16], and each (weighted) expert is proportional to a multivariate t-distribution [54] over RD_ In general, it can be challenging to work with products of experts, but for this family we show how to sample from and (in some cases) evaluate their censsttis-exaccio what is needed to use them for VI.
The full potential of this family is unlocked by a Feynman identity [12, 46], originally developed for loop integrals in quantum field theory, that expresses the product of K fractions as an integral over the simplex A KK-1 In particular, we show that the Feynman identity implies representation of the product of t densities as a continuous mixture of t-distributions. A consequence of this representation that a PoE in this family is more naturally equipped than a finite mixture model to approximate target densities with a continuum of modes that lie in convex set. We also leverage this representation to perform two important tasks for VB this variational f'amily, and second, how to estimate its normalizing constant. Notably, we transform these problems to samplers and integrals over the simplex A K -K1 as opposed to all of RRD.
how to draw samples from a PoE in
Using these techniques, we then introduce a score-based BBVI algorithm to find the member of this PoE with the minimum Fisher divergence to the target density. To do so, we first generate a large pool of experts (i.e., t-distributions with different modes and tails), and then we derive an iterative score-matching procedure to optimize the exponents that determine their geometric weighting in the PoE. In practice this procedure drives many exponents to zero, thus pruning irrelevant experts from the product. More specifically, each iteration updates the expert weights in the PoE by minimizing a regularized Fisher divergence. Score matching is particularly convenient for PoE models because the score is linear in the weights, and the optimization problem reduces to solving a sequence of convex quadratic programs, where each subproblem can be solved efficiently.
In addition, we analyze the convergence of the BBVI algorithm. First, we derive rates of convergence for the expert weights that are iteratively estimated by the algorithm with finite batch size. We prove that these weights converge exponentially quickly neighborhood of an optimal weighting of experts in the PoE, where the size of the neighborhood depends on the amount of misspecification of the variational family. We also demonstrate the benefits of this approach to BBVI empirically on a variety of synthetic and real-world target densities. In particular, we show that with this variational family, we can approximate diverse target densities, including ones that are skewed and heavy-tailed.
2
latent variable model for products of experts
Our goal perform score-based BBVI with a particular PoE variational family by minimizing the Fisher divergence between a variational density ((z) and a target density p(z) supported on RD:
((::p)
| log ((2) - V logg(=))1 ((=) dz.
Before tackling the broader problem of BBVI with products of experts, we begin by studying these models in their own right, showing how to enable the use of the PoE as variational family. Consider a (weighted) PoE with the density
K
(())
II (z) k=
(2)
Ca
where the exponents {ax K determine the geometric weighting of the nonnegative functions {qx 1K_ in the product, and C. is normalizing constant given by
K II (x)))² dz. 1=1
Ca /a=/In
(3)
We refer to the function qj as the kth expert in the PoE and to the exponentiated function of" as the kth weighted expert. We assume that the exponents (or weights) in the PoE are nonnegative (i.e., ag 2 0), and later we identify further constraints that ensure the integrability of the product in Eq. 3.

---

Product (w= =2)
Product (4= =5)
(μ =5)
0.3,
d: a;
0: 0.3, 0: 1.0
density
0.
Figure Product vs mixture of two t-distributions. The experts (gray dashed curves) all have the same scale of 1. The weights in the mixture (rightmost) are the normalized or values.
When the expert functions gA Gaussian, the product itself is Gaussian, leading closed-form normalizing constant and efficient sampling (see, e.g. Wu and Goodman [56]). However, more generally,
be difficult to work with products of experts, and in particular to sample from their densities or to estimate their normalizing constants. Because the Fisher divergence Eq. is generally intractable,
need samples from form empirical estimate of the divergence that can then be minimized. In addition, while the Fisher divergence does not rely on the normalizing constant of q, there many applications that require it; in this work, Kullback-Leible (KI.) divergence Section 4.
it compute the
In this section, develop particular family of products of experts, show how to reformulate the models in this family as latent variable models, and then (most critically) leverage the latent space in these models to draw samples from their densities. The ability to sample efficiently from these models will be at the heart of their use for BBVI in Section
2.1 Products of multivariate (-distributions
We now focus on the parameterized family of products of experts whose weighted experts are proportional to multivariate /-distributions over RDD. Specifically, we suppose that
(=(=)**
uR)) As(2 Ax(==1uz)
where μ4 ERD and A₂ 20 O. Recall that I-distribution is parameterized by its mean μ RD,, inverse scale matrix Q2> > 0, and degrees of freedom v > 0, and has the probability density function
D( 12/1/2
T(21m, = (TV) [1 (z- (u+D)/2
(5)
where I' is the gamma function and | denotes the determinant Ω. Thus the kth weighted expert Eq. is proportional to t-distribution with mean μk. inverse scale matrix Aş (when Aş 0), and degrees of freedom 2a D.
Note that do require all of the weighted experts in Eq. define proper densities; instead we allow some of the inverse scale matrices to be rank-deficient with |At/=0, Thus these products of experts have more flexibility than finite mixture models whose component densities must be normalizable. The PoE in Eqs. 2 and 4 will be normalizable if M > D. provided all inverse scale matrices Ad full rank. Consequently, we impose this linear sum constraint the expert weights. If not all inverse scale matrices are positive definite, then integrability must be verified.
We demonstrate the flexibility of this family in Figure 1, where the first two pancls show the product of two t-distributions that take the form 9(2) I 1 parameters are set to either wide range of behaviors. Finally, we plot the mixture of the same two experts as the middle panel, where the mixture weights are formed by normalizing the values. These two panels highlight the difference between products- whose "and" relationship yields larger density values when all experts large-and and mixtures-whan "or" relationship yields larger values when any expert is large.
(2- p )] -a where the location +3 or μ +5. With only two experts they already parameterize

---

of 2.2 Feynman parameterization for products of experts
We show how to rewrite the PoE in Eq. 2 in particularly revealing form. To do so, we rely on an identity that expresses a product of positive denominators as an integral over the sim- plex A - {w E [0, 1]K Ex wy = 1}. The simplest form of the identity (for two denominators) is straightforward to verify; it states that
A Jo (Atw+ + A2(1-w)) dw
for all A,, A2>02, The identity we generalizes the above to product of k denominators that geometrically weighted by exponents og- In this case, the integral over [0, 1] in Eq. 6 is replaced by an integral over the simplex AK- In particular, is also true that
T(Cx on) I.. ((ox) S.x. (Ex ux
06-1 II. w
dw,
A AK
This so-called Feynman parameterization is useful to simplify certain loop integrals in quantum field theory [12. 46]. Here we make telling observation about the terms in Eq. 7 that are independent of the denominators As; these terms are equal to the probability density function of the Dirichlet distribution. In particular, we can express the Feynman parameterization more compactly as
4 aL = Eiraa)((a))
(Ea
En
(8)
We now use the above identity to reinterpret the PoE introduced in the previous section. First we set A =1/x((2)) (2) in Eq. 8, where g& (2) is the expert in Eq. 4. Also, as shorthand, let ||>|| v Av denote the quadratic forms that appear in these experts. From these substitutions follows that
K
[ (1+ wwllz- paall--Mxxl3, 113.
Σn cx
II
go (z) Binaa)
(9)
We now show that the PoE can be re-expressed as a continuous mixture of multivariate t-distributions, all of which have the same number of degrees of freedom but differ in other respects; Appendix C.1 for a detailed derivation. To do so, we define
A(w) = "x
(10)
p(w)= = (wwww--E WgAApMa
(11)
which can be viewed as location and inverse scale parameters that are continuously indexed by W A K-1, Now, after expanding the quadratic form in the denominator of Eq. 9 and completing the square, we can use the definitions of A(w) and p(w) to rewrite Eq. 9 as
(1+ [1+EeuRllne (l)) + +1 ((=))()))))( -E ax
II (=) B.-DDiia)
(12)
Next we supplement Eq. 10 by defining another inverse scale matrix, ((w), that is continuously indexed by w A K-1 and absorbs the terms in Eq. 12 that are independent of 2. Let
v = 2 Σx k Ork - D, )²(ww = Σx wr I [())
(13)
(14)
((w) (A(ww/// (1=
(15)
We can now cast the expectation in Eq. 12 into a more revealing form, arriving the following representation of the PoE as a continuous mixture of t-distributions. This result is obtained by making the substitutions Eqs. 13 15 and appealing the form of the t-distribution Eq. 5.
Result 2.1 (PoE as continuous mixture of t-distributions). Consider the product II 9a (=))0x with t-distributed experts (Eq. 4), and let 2 Ex an -D. D. The product can be computed as
F Fw-DD((a) =
II k=1
(1+0
2+D
(1
(16)
Thus the PoE can be viewed as continuous mixture of t-distributions in Eq. 5 with v degrees freedom, location parameters (uw), and inverse scale matrices ((w), where w E K-1

---

This result suggests that this type of PoE may be well-suited for distributions with continuum of modes that lie convex set. Any such regions of high probability will be more naturally modeled by continuous mixture than finite one. As an aside, we note that the I-distribution can itself be written as (continuous) scale mixture of Gaussians [1]. By extension, the above result shows that product of t-distributed experts can be written even more generally as continuous location-and-scaa mixture of Gaussians. This correspondence provides further motivation for their use in VI.
2.3
Joint density and latent variable model
There are many useful results that follow from the Feynman parameterization, One such result arises when revisiting the computation of the PoE's normalizing constant Ca in Eq. 3. While C. is initially expressed as an integral over RD, our goal is to re-express as a (potentially simpler) integral over the simplex A K-1, Substituting Eq. 16 into Eq. 3, we find that
(w)) 4 / dz - (1=
C. = (iwDa(()))
(17)
where we have applied Fubini's theorem to move the integral over inside the expectation over w. can now perform the integral over z, as is given exactly by the normalizing constant Eq. 5 for a I-distribution with inverse scale matrix ((w) and degrees of freedom. In this way we obtain the following result.
Result 2.2 (Normalizing constant of PoE). Consider the PoE with the t-distributed experts in Eq. 4, and let == Σx -D The normalizing constant of this PoE can alternately be computed as
C. / 11 (2) dz Biraa)
[ (1+07(w))
(TE) r( T(OtD)
(18)
Next we use the Feynman parameterization to show how sample from q. In particular, construct a joint density g(w, 2) that yields the desired marginal ((2) in Eq. 2, and then we generate samples from the joint q(w, z). Combining Result 2.1 and Result 2.2, we find
C.(wr) Ca
q(w, 2)
Dir(w ( T())
(19)
where the new leading factor (w) in the numerator is used to account for the terms in Eq. that are not absorbed by the (normalized) Dirichlet and t-distributions t-distribution--nmelly
17
((5) (TT) r( P)
(w) = (1+o²)))
(20)
is then straightforward to verify that q(w, (00.2)dd 2)dw q(z), leading to natural auxiiiary-variable sampling procedure. Indeed, we can interpret this joint density as a latent variable model and use its marginal and conditional densities to draw samples from the PoE. Marginalizing over z in Eq. 19,
find that q(w)= = J [qw, 2) dz = C. recover the !-distribution for the conditional density q(2) w) upon dividing the joint in Eq. 19 by this result. In we have shown the following.
Cu(w) Dir(w a)), where Ca (w) is given by Eq. 20. Likewise, we
Result 2.3 (Latent variable model for PoE). Consider the PoE with the t-distributed experts in Eq. 4. We can draw samples from the PoE by sampling from the generative model
Wh (w(w)
Dir(w (a),
(21)
C.
2b Wb ~ T( [<(]) I(wx),
(22)
where Wy lies in the simplex, the de-depeddent terms on the right are given by Eqs. 10, 11, to and 20, and 2, is conditionally i-distributed given
5

---

Identify mode(s)
Refine mode(s)
Fit variational parameters
025
larga
sspert [iter
verieloral aporuamaton
(mode)
0.15
0.10
0.05
0.00
2
2
Figure 2: To select experts, we identify each mode and then add more experts to refine the fit.
The above result not yet a practical recipe for sampling from a PoE. The difficulty lies in the first step: it is not straightforward to sample from ((w) in Eq. 21 due to its leading dependence on Ca Ca(w). But we can circumvent this difficulty by using the Dirichlet distribution that appears in q(w) as proposal distribution for importance sampling.
Suppose we wish to estimate an expected value
We can draw batch of samples
Wh ~ Dir(w (x)
(23)
26 W/ T( ([w)), (w.))
(24)
from the Dirichlet and t-distributions in Eqs. 21 and 22 and weight these samples by Co (wx) in the calculation of the expected value. In this way we obtain the estimate
En Ca((w)) h(z))
E (x)))((2)]]
(25)
Es (un))
3 Score-based VI with products of experts
In this section we show how to approximate a target density by a PoE with t-distributed experts. To do so, we must specify the number of experts, the parameters of their t-distributions, and the exponents that determine their geometric weighting in the PoE. We discuss these problems in turn.
3.1 Selecting the experts
Of the many ways to select experts, we seek a simple heuristic that works in practice and avoids complicated, coupled optimization over the expert parameters (Mk, A&) and weights og in Eq. 4. Our basic strategy is to generate large, oversaturated pool of experts whose means 14 are concentrated near the modes of the target density. We shall see in later sections that poorly situated (and hence irrelevant) experts are efficiently pruned by the procedure for learning the weights aq.
Our strategy for selecting experts has three steps. The first is to locate the modes of the target density by hill-climbing in log p(2) from randomly chosen starts; if the target density a Bayesian posterior, then we can choose these starts by sampling from the corresponding prior. The second step is to place an expert each mode: the expert's mean μ4 and inverse scale parameter Aş to match the mode's 's location and curvature. The third step is to place additional experts nearby so that the weighted PoE can better model the shape (e.g., skew, kurtosis) of each mode: once new location μ4g is chosen, the corresponding inverse scale A₁ is set to the (negative) Hessian at l4 projected to the cone of positive semidefinite matrices. For further details of this step and a discussion of expert placing cost, see Appendix D.I. Figure 2 illustrates the intuition behind this strategy for target density with one mode, Overall we found this strategy to be quite effective in conjunction with a score-based procedure for optimizing the expert weights in the PoE. We describe this score-based procedure next.
3.2 Weighting the experts
Given experts qu in the form of Eq. 4, we now consider how to form the weighted PoE in Eq. 2 that best approximates the target density p. Specifically, for 4(2) 1I (2) we seek the weights (aw) that minimize the Fisher divergence between q and p, defined in Eq. 1. Since the PoE in Eq. 2 has support on DD, this divergence vanishes only when =pp

---

Though we cannot minimize Eq. 1 directly, in the spirit of BBVI, we can instead attempt to minimize an empirical estimate of the Fisher divergence, one that is based on drawing samples from q. A further simplification is achieved by decoupling the procedures for sampling from q and optimizing g. To do so, we iteratively minimize a closely related objective. In particular, let a E RK be the expert weights at the t iteration, and let q(z (a) denote the density of the PoE with expert weights a. Rather than minimizing Eq. 1 directly, at the th iteration we instead solve the simpler problem
(1+1) aargmin 6
{ V
Vlogp(z)]] ((2)a!") dz +
my 2/ ='') (26)
where 7/ > a learning rate and the domain & of the optimization constrains the expert weights to define normalizable PoE. Note that this update minimizes a biased estimate of the Fisher divergence; the estimate is biased because in the first term of the objective the expectation is performed with respect to (=) instead of (=]aa). But at the same time, the update attempts to compensate for this bias by penalizing solutions that move too far from one iteration to the next; this penalty is enforced by the regularizer in the second term of the objective. This the same intuition that is behind a recently proposed aatch-and-mattee algorithm for Gaussian BBVI [7].
The rest of this section fleshes out this iterative procedure and highlights its three main advantages. First, when is a weighted PoE with I-distributed experts, we can use samples to compute an empirical estimate of the objective in Eq. 26. Second, at each iteration, we can minimize this empirical estimate by solving strongly convex quadratic program. Third, this iterative procedure provably converges under fairly general conditions to a neighborhood of an optimally weighted PoE.
We now construct an empirical estimate (a(a) for the objective in Eq. 26 at the 4t iteration. Using the latent variable model in Eqs. 21 and 22, we generate a batch of B weighted samples {[w,, 2x) =))]m=
from the PoE with expert weights a (+), From these samples, we construct the empirical estimate
((()) =
Σ Th v loga(z)la) - Vlogg(>)))1 + 1l a01
(27)
where To X Caso (wy) is the importance weight of the 8th sample from Eq. 25. Note how by sampling the PoE with weights (1), we have decoupled these samples from the optimization over a in Eq. 26.
Next show that the empirical estimate (a) is minimized by solving a convex quadratic program in the expert weights a. First, we observe that for any PoE, as defined by ((=(a) x II qu (=)** in Eq. 2, the score is linear in the weights of its experts: namely, V log ((=(a)) Ex k Oa V log (z). We make this explicit for the I-distributed experts in Eq. 4 by writing
V Plogg(zs(a) = Qba,
(28)
where Qo is the Dx K matrix whose 4tth column is given by V log qu (=6) = 22x (z))A<(z))))))))()))))))))) pk)) From the linearity of the scores, it follows at once that El(a)) Eq. 27 quadratic in the expert weights a. In particular, we can rewrite Eq. 27 as
(a) = a. [ E.m. Qa ] a- OF V logp(=s) +
a + {,),,, (29) mt [1]
and from the above, also that (a) is strongly convex in for all ni E (0, oo). The expert weights are updated by minimizing this objective subject to a constraint that the newly weighted PoE is normalizable. To satisfy the parameter constraints of the PoE, we define the constraint set, such that a E 6, as
6 f a 20, az 20.....
as 2 0, Eeox 2 2 +
(30)
where E is a slack variable added ensure that the product of experts Il (2()))n integrable. Eq. 30 defines a convex set, and hence the overall optimization is convex; in particular, it is nonnegative least squares (NNLS) problem with linear constraints, for which there exist many efficient solvers [29. Ch. 23]. For our purposes, it is also interesting that problems in NNLS often yield sparse solutions where multiple constraints are active. In our setting, these are solutions in which irrelevant experts are assigned zero weights and do not contribute to the density of the PoE.
7

---

Algorithm Score-based VI for PoE: learning the weights o
1: Input: initial weights a (0) E E, target score function V log p. experts qN (z), learning rates 1/2 > 0.
2: for = 1,...., T: do
3: for b 1, :, B: do
4: Sample wa ~ Dir(w I ()) and 26 ~ [( p(ww), (((s)),())
5: Compute importance weight Tb & C, (ws) and target score gb = V logp(zs).
6: Compute matrix Q4 E DRX whose kth column is equal to V log (x()).
7: end for
8: Compute G. = E m0Q Q6 + 2 and h1 = = m 90 + m aa
d5
9: Update expert weights: a (1+1) argminoee [ Ga--ifa a).
end for
11:
Output: PoE weights (T) E &6 (Eq. 30)
3.3 Convergence theorem
We summarize the overall iterative procedure for learning expert weights in Algorithm 1. To prove convergence we make some basic assumptions about the Fisher divergence Eq. and its empirical estimate at each iteration of Algorithm 1. To state these assumptions, let
(a)a) = 1 E V log (z((aa) - V log p(z))
(31)
denote the first term in Eq. 27, and let V&, and H [.] denote the gradient and Hessian of this term with respect to the expert weights a. With these definttions, we can state the following theorem.
Suppose that 2(; p) in Eq. 1 is minimized by a unique a EC, and also that for all 2 0 there exists some 6'20 0 such that B 1//// V (ar))) < 8 and some >> 0 such that H [[ > AZ almost surely. Then for constant learning rates 1/+ E1) the expected error of the iterates satisfies
Theorem 3.1.
Ella! - a= (1+ ||a.0
- a²Il
(32)
This result ensures that the iterates a of Algorithm converge, in expectation, to a neighborhood around the optimal expert weights a", with the error decaying at a linear (geometric) rate. The proof relies a few key ideas, most notably that (i) the constrained least-squares problem in Eq. C.6 can be solved by projecting its unconstrained solution onto C, and (ii) this projection is with respect to an induced Mahalanobis norm, and it nonexpansive [2, Proposition 4.8]. See Appendix E for complete proof.
The bound in Eq. 32 separates two effects-the transient error, which shrinks exponentially fast with I and depends on the strong convexity parameter X, and the asymptotic error floor & which depends also on the misspecification parameter We briefly sketch the intuition behind these terms.
First, the error floor is proportional to the misspecification parameter ó. Let a denote the optimally =|V = O for all I; we see in
weighted PoE. On one hand, if P Q, then q" - p and this case that 8=0. On the other hand, if p g Q, then q" Zp. this case we expect that the stochastic gradients will have small norms at a² (and thus 8 will be small) whenever 2(4": itself is small. In the right panel of Figure D.2 we highlight on one experiment with an sinh-aresinh target, that for constant or decreasing step size schedules, the error floor does not go below 10-1 However, on the left panel of Figure D.2, we do show that this error floor improves as we increase the batch size. See Appendix D.2 and Appendix D.3 for details.
Second, the transient error depends the assumption that (a)a) is X-strongly convex for all t. Observe from Eq. C.6 that M[ 1 Eee TTQ Qo: i.e., the Hessian sum of B positive semidefinite matrices. Thus for larger batch sizes, the assumption of strong convexity is increasingly likely to be satisfied. We confirm this intuition in Figure E. 1, where we show that as we increase the batch size, the eigenvalues of the Hessian also increase and are bounded away from zero. We only encountered one indefinite Hessian in our experiments; this occurred when the number of experts (K = (K0) was much larger than the batch size (B= = 10).
B
8

---

Fisher
KL.
KL
Diwargen ganoo
Divergence
Ducrgence
400
Bold
Riee
Pdgl
A0M
Bold
Flowr
ADM
Bold
Pioee
PrE
ADM
Bold
(a) Gaussian mixture
(b) Product of experts
Diamond
(d)
Funnel
Figure 3: Synthetic 2D targets. Top: Gray contours represent the target, and blue points represent
samples from the fitted PoE. Bottom: The KL and Fisher divergences each method.
4 Experiments
In this section, we compare VI with this PoE family to Gaussian BBVI (with ADVI [28] and BaM [7)) and normalizing flow-based VI [43]. In Appendix F, we provide more details on these experiments and also results on several additional examples. We also present additional results for estimating PoE normalizing constants and evaluating sampling (Appendix F.I) and expert placement (Appendix D.2).
4.1 Synthetic 2D targets
We first consider several synthetic 2D target distributions: 1) Mixture of Gaussians, 2) Product of t-distributions, 3) Diamond, 4) Funnel. We report Monte Carlo estimates of two training-objective- agnostic metrics, KL(p: 9) and Ep log v log pl"), using 1000 samples from p. (The VI methods minimize the expectations with respect to g.) We visualize the PoE samples by first drawing weighted sample and then resampling the generated samples according to their normalized weights. For full details on each target, see Appendix
Figure 3 shows each target distribution's 's contours (gray curves) overlaid with samples from the fitted PoE (blue points). In these examples, the PoE achieves substantially lower values for both divergence metrics than the Gaussian approximation, due to its ability to better capture the tails of the target distributions. The normalizing flow also provides substantial improvement over the Gaussian fit. Notably, for heavy-tailed (product of experts and diamond) targets, the PoE outperforms the flow-based family as well, yielding lower divergence in these cases.
4.2 Sinh-arcsinh target
We now evaluate a higher-dimensions synthetic target distribution with skew and heavy tails. The sinh-aresinh distribution [21, 22] generalizes the Gaussian distribution with additional parameters e E R and r E R++ that control the skewness and the tail-weight. It transforms a Gaussian draw 2 (0, )) coordinate-wise via 2y
(a.)) + TT)) We construct the target density in D = 50 dimensions so that is positively skewed (e = 0.3) and heavy-tailed 077)..7).
In Figure 4, we show the KL and Fisher divergences between all approaches, and we found that the PoE had the lowest divergence for both metrics (left). We also plotted the divergences against the number of gradient evaluations for the iterative parameter fitting of the PoE and the normalizing flow. Note that this plot does not show the initial startup costs of each method (selecting the experts for the PoE and selecting learning rate for the flow), which are not the dominating cost. In addition, this plot only shows the first x 10d gradient evaluations. For both divergences, the values decrease slowly for the flow model; the final reported value in the bar graph is after 10" gradient evaluations.

---

50
K0.
10
Flow
Fisher divergence
Flww
Divergence
ADV
Baad
Plowr
Pot
Swi0f
10xlc
15x10
5x10*
10x20"
15xL0"
gradent evaluations
(a) All methods
(b) Divergence vs # of gradient evaluations
Figure 4: 50-dimensional sinh-arcsinh target with skew and heavy tails.
8
Negabbe LLH
2.0
Fannt chargern
8
Nogattoe LLM
Finner Coreereee
Fener chrrgenncs
1.5
Nagathe
ADM
Bad
Plae
PdB
Bal
Ploee
PdE
AD12
Bbdd
Fiaee
ADV1 Dald
Pioe
PaE
ADV Bas
Fiee
Puel
ADV
Bald
Floe
(a) garchi1
(b) 8schools noncoleeenneeeeeed
gp-regr
Figure 5: posteriordb targets that highlight a range of non-Gaussian posterior properties.
Posterior inference problems
Next we study several posterior inference test problems from posteriordb [32, 33]. each target, we use reference samples computed using Stan (drawn via Hamiltonian Monte Carlo) [9]. Because we do not have access to the normalized target density p(z), instead of reporting KL q). we report the negative expected log likelihood Reg-LLH(p:: q) = [log(()]) We again also consider the Fisher divergence F(p q), which does not require normalizing constants. All expectations are estimated using the reference samples from p.
In Figure 5, we highlight three particular examples from posteriords First, we consider garch11: this model has (half) uniform priors that skew the posterior when the support is transformed to R... The skew not modeled by Gaussian VI, but it present the PoE and flow-based approximations. Next, we consider dechoolssnonneenter introduced by a Half-Cauchy 5) prior. These properties are best modeled by the variational PoE. Finally, we study tailed example in gp-regr, where Gaussian BBVI performs well due the near symmetry of the posterior. Nevertheless, despite its heavy tails, the product of t-distributions still outperforms both Gaussian BBVI and the normalizing flow (with Gaussian base distribution) this example. This example illustrates how the expert weights a can optimized to match different tails via the value of Σ, a% D in the mixture representation of Result 2.1.
which exhibits skew and a heavy-tailed component
Discussion of contributions, limitations, and future work
We have shown that VI with products of t-experts can model a variety of target densities. A key technical insight, via Feynman identity, was represent each PoE continuous mixture indexed by a Dirichlet variable. We then used this representation to sample from the PoE, a core requirement for BBVI. To optimize the expert weights, sequence of convex quadratic programs, and we proved that its iterates converge exponentially to neighborhood that depends the degree of misspecification of the variational family.
developed score-based VI algorithm that solves
While our algorithm learns the expert weights, it is limited by relying on a fixed collection of experts selected upfront. In practice, the quality of the approximation depends heavily how these experts are chosen. While our current heuristic overspecifies the number of experts and prunes away the ineffectual ones, further gains could be achieved with more refined approach- for instance, boosting style approach [15, 35] that sequentially adds experts based on (say) score mismatch.
Several promising directions remain. First, the Feynman identity may have broader implications: for example, by providing semi-analytic procedure estimate the normalizing constant, it may open a new door to generative models with products of (-experts [26, 54]. Second, the PoE construction may be useful in certain sampling methods as an initialization proposal scheme. Finally, as mentioned above, by choosing experts carefully can hope to accelerate every aspect of our approach.
10

---

Acknowledgements
D.B. was supported by NSF IIS-2127869, NSF DMS- 2311 108, ONR N000142412243, and the Simons Foundation.
References
[1] D. F. Andrews and C. L. Mallows. Scale mixtures of normal distributions. Journal of the Royal Statistical Society: Series B (Methodological), 36(1):99-102 -029 1974.
[2] H. H. Bauschke, P. L. Combettes, H. H. Bauschke, and P. L. Combettes. Convex analysis and monotone operator theory in Hilbert spaces. Springer, 2017.
[3] K. Bhatia, N. L. Kuang, Y.-.. Ma, and Y. Wang. Statistical and computational trade-offs in vari- ational inference: case study inferential model selection. arXiv preprint 1rXiv:2207.11208, 2022.
[4] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859-8777 2017.
[5] A. Cabezas, A. Corenflos, Lao, R. Louf, A. Carnec, Chaudhari, R. Cohn-Gordon, Coul- lon, W. Deng, S. Duffield, et al. BlackJAX: composable Bayesian inference in JAX. arXiv preprint arXiv:2402.1073 2024.
[6] D. Cai, C. Modi, C. Margossian, R. Gower, D. Blei, and L. Saul. EigenVI: score-based variational inference with orthogonal function expansions. In Advances in Neural Information Processing Systems, 2024.
[7] D. Cai, C. Modi, L. Pillaud-Vivien, C. Margossian, R. Gower, D. Blei, and L. Saul. Batch and match: blxc-boox variational inference with a score-based divergence. In International Conference on Machine Learning, 2024.
[8] T. Campbell and X. Li. Universal boosting variational inference. Advances in Neural Information Processing Systems, 2019.
[9] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell. Stan: A probabilistic programming language. Journal of Statistical Software, 76(1):1-32, 2017.
[10] L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real NVP. In International Conference on Learning Representations, 2017.
[11] J. Domke, G. Garrigos, and R. Gower. Provable convergence guarantees for black-box varia- tional inference. Advances in Neural Information Processing Systems, 2023.
[12] R. P. Feynman. Space-time approach to quantum electrodynamics. Physical Review, 76(6): 769-789, 1949.
[13] G. Garrigos and R. M. Gower. Handbook of convergence theorems for (stochastic) gradient methods, 2023.
[14] S. Gershman, M. Hoffman, and D. Blei. Nonparametric variational inference. In International Conference on Machine Learning, 2012.
[15] F. Guo, X. Wang, K. Fan, T. Broderick, and D. B. Dunson. Boosting variational inference. arXiv preprint arXiv: 161 2016.
[16] G. E. Hinton. Products of experts. In International Conference on Artificial Neural Networks, 1999.
[17] G. E. Hinton. Training products experts by minimizing contrastive divergence. Neural Computation, 14(8):1771-1800,
2002. 2002.

---

[18] J. Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 2020.
[19] A. Hyvärinen. Estimation of normallzdd statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
[20] P. Jaini, I. Kobyzev, Y. Yu, and M. Brubaker. Tails of Lipschitz triangular flows. In International Conference on Machine Learning, 2020.
[21] C. Jones and A. Pewsey. Sinh-aresinh distributions. Biometrika, 96(4):761-780, 2009.
[22] C. Jones and A. Pewsey. The sinh-aresinh normal distribution. Significance, 16(2):6-7, 2019.
[23] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul. An introduction variational
methods for graphical models. Machine Learning, 37(2):183-233, 1999.
[24] K. Kim, J. Oh, K. Wu, Y. Ma, and J. R. Gardner. On the convergence of black-box variational inference. Advances in Neural Information Processing Systems, 2023.
[25] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
[26] D. P. Kingma and M. Welling. Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
[27] D. P. Kingma, T. Salimans, R. Jozefowiez, X. Chen, I. Sutskever, and M. Welling. Improved vari- ational inference with inverse autoregressive flow. Advances in Neural Information Processing Systems, 2016.
[28] A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei. Automatic differentiation variational inference. Journal of Machine Learning Research, 18(14):1-44 2017.
[29] C. L. Lawson and R. J. Hanson. Solving least squares problems. SIAM, 1995.
[30] Liang, M. Mahoney, and L. Hodgkinson. Fat-tailed variational inference with anisotropic tail adaptive flows. In International Conference Machine Learning, 2022.
[31] F. Locatello, G. Dresdner, R. Khanna, I. Valera, and G. Rätsch. Boosting black box variational inference. Advances in Neural Information Processing Systems, 2018.
[32] M. Magnusson, J. Torgander, P.-C. Bürkner, L. Zhang, B. Carpenter, and A. Vehtari. posteriordb: Testing, benchmarking and developing Bayesian inference algorithms. 2025.
[33] M. Magnusson, P. Bürkner, and A. Vehtari. posteriordb: set of posterions for Bayesian inference and probabilistic programming. https://githut com/ /ttp///ttaadd////seerrottdtdddd 2022.
[34] L. Martino, V. Elvira, and F. Louzada. Effective sample size for importance sampling based on discrepancy measures. Signal Processing, 131:386-401, 2017.
[35] A. C. Miller, N. J. Foti, and R. Adams. Variational boosting: Iteratively refining posterior approximations. In International Conference on Machine Learning. PMLR, 2017.
[36] C. Modi, C. Margossian, Y. Yao, R. Gower, D. Blei, and L. Saul. Variational inference with Gaussian score matching. Advances in Neural Information Processing Systems, 2023.
[37] C. Modi, D. Cai, and L. K. Saul. Batch, match, and patch: low-rank approximations for score-based variational inference. In International Conference on Artificial Intelligence and Statistics, 2025.
[38] W. Morningstar, S. Vikram, C. Ham, Gallagher, and J. Dillon. Automatic differentiation variational inference with mixtures. In International Conference on Artificial Intelligence and Statistics, 2021.
[39] E. Moulines and R. Bach. aon-asymptotic analysis of stochastic approximation algorithms for machine learning. In Advances in Neural Information Processing Systems, 2011.
12

---

[40] D. Needell, N. Srebro, and R. Ward. Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm. Mathematical Programming, Series A, 155()):54--573 2016.
[41] F. Orabona. A modern introduction to online learning. arXiv preprint
2025.
[42] R. Ranganath, S. Gerrish, and D. Blei. Black box variational inference. In Artificiaa Intelligence and Statistics. PMLR, 2014.
[43] D. Rezende and S. Mohamed. Variational inference with normalizing flows. In International Conference Machine Learning, 2015.
[44] E. Roualdes, B. Ward, S. Axen, and B. Carpenter. BridgeStan: Efficient in-memory access to Stan programs through Python, Julia, and R. https:///iithb: httpehuuu.cmmmmuueeeeeeee http:////oooaeeeeerereeeeeeeee 2023.
[45] Y. Shao, S. Yu, and Feng. Nonparametric automatic differentiation variational inference with spline approximation. In International Conference Artificial Intelligence and Statisticy, 2024.
[46] V. A. Smirnov. Evaluating Feynman integrals. Springer Tracts in Modern Physics, 211:1-244, 2004.
[47] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning, 2015.
[48] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 2019.
[49] Y. Song, S. Garg, J. Shi, and S. Ermon. Sliced score matching; A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence, 2020.
[50] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd. OSQP: An operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4):637-672, 2020.
[51] M. K. Titsias and F. Ruiz. Unbiased implicit variational inference. In International Conference on Artificial Intelligence and Statisties, 2019.
[52] A. Uppal, K. Stensho-Smidd, W. Boomsma, and Frellsen. Implicit variational inference for high-dimensional posteriors. Advances in Neural Information Processing Systems, 2023.
[53] M. J. Wainwright, M. I. Jordan, al. Graphical models, exponential families, and variational inference. Foundations and Trends@ in Machine Learning, 1(1-2):1-305, 2008.
[54] M. Welling, S. Osindero, and G. E. Hinton. Learning sparse topographic representations with products of student-t distributions. Advances in Neural Information Processing Systems, 2002.
[55] K. W. k. Wong, M. Gabrié, and D. Foreman-Mackey. flowMC: Normalizing flow enhanced sampling package for probabilistic inference in JAX. J. Open Source Softw., 8(83):5021, 2023.
[56] M. Wu and N. Goodman. Multimodal generative models for scalable scakly-supervised learning. In Advances in Neural Information Processing Systems, 2018.
[57] Z. Xu and T. Campbell. The computational asymptotics of Gaussian variational inference and the Laplace approximation. Statistics and Computing, 32(4):63, 2022.
[58] Z. Xu, N. Chen, and T. Campbell. MixFlows: principled variational inference via mixed flows. In International Conference on Machine Learning, 2023.
[59] M. Yin and M. Zhou. Semi-implicit variational inference. In International Conference on Machine Learning, 2018
[60] D. Zoltowski, D. Cai, and P. Adams. Slice sampling reparameterization gradients. Advances in Neural Information Processing Systems, 2021.
13

---

Appendix
A Overview
In this supplementary material we provide additional details several aspects of the paper. We start by providing a discussion of the broader VI literature in Appendix B, focusing on VI methods that employ expressive variational families. Then, in Appendix C provide additional derivations and details about the Feynman parameterization. In particular, we provide a full derivation of the continuous mixture representation of the PoE presented in Result 2.1, along with an expanded discussion of the normalizing constant. Next, in Appendix D we present one heuristic used for selecting the experts, along with an empirical study of the heuristic. In Appendix E, we provide the proof for the convergence theorem, which is supplemented with an empirical investigation. Finally, in Appendix F, we present additional empirical results. In particular, we first study the properties of the PoE family, including the estimation of the normalizing constant and the sampling procedure. We then provide details about the setup of the experiments Section 4, including more details about the tuning of the baselines and the implementation of the PoE score matching algorithm. Next, we provide additional plots and analysis for the experiments Section 4.
B Related work
Mixture modeling. There have been many efforts to extend VI beyond families of factorized or multivariate Gaussian distributions. One natural extension is to choose the approximating distribution from a parameterized family of mixture models (i.e., a mixture of experts), and in this line of work, it has been most common study Gaussian mixture models. Gershman et al. [14] introduced an infinite Gaussian mixture model that draws inspiration from kernel density estimation. In a different approach, Morningstar al. [38] integrated mixture families into ADVI using stratified sampling and optimizing a tighter evidence lower bound via gradient descent. More recently, Xu et al. [58] combined the variational families of mixture models and normalizing flows, establishing MCMC-like convergence guarantees.
Boosting VI approaches. Mixture models have also been used in boosting strategies for VI. In this approach, the variational approximation is iteratively improved by adding new components to the mixture. Guo et al. [15] and Miller et al. [35] proposed greedy algorithms that incrementally add Gaussian components to reduce the KL divergence. Locatello et al. [31] later reinterpreted this boosting procedure as Frank-Wolfe algorithm, deriving explicit convergence guarantees. Campbell and Li [8] replaced the KL divergence with the Hellinger distance; with this objective, the optimization of mixture weights reduces
problem in nonnegative least squares that is similar in flavor the optimization of PoE weights in our work.
Basis-set variational families. A related line of work has explored variational approximations that are parameterized by linear combinations of basis functions. For example, in Eigen [6] the variational approximation is parameterized as the square of a linear combination of orthogonal functions, and the best approximation is found by minimizing a Fisher divergence- an optimization that reduces in turn to minimum-eigs eigeevaale problem. However, EigenVI scales exponentially with the latent dimension, limiting its applicability. Shao et al. [45] used spline-basse approximations within an ADVI framework, increasing the flexibility of the variational approximation.
Normalizing hosebased families. Another line of work has explored the use of normalizing flows [27, 43), where the variational approximation is parameterized by an invertible neural network and found by optimizing the ELBO using the reparamcterization trick. Jaini et al. [20] and Liang et al. [30] introduced tai--adaptive components into flow-based models in order to better approximate non-Gaussian posteriors with heavy tails. Normalizing flows with these architectures have enjoyed empirical success, but they can be difficult to optimize and generally lack convergence guarantees when they are used for VI.
Implicit variational inference. Implicit VI [51, 52, 59] defines variational families using con- tinuous mixture of a kernel. The product of experts we consider can be viewed as fitting in the framework of implicit VI with particular kernel (l-distribution) and mixing measure (Dirichlet) due to Result 2.1; the PoE here is particular function of the auxiliary variables. On the other hand, the typical implicit VI setting specifies a neural network with a mixing distribution over the parameters
14

---

of that neural network. In this special case of continuously parameterized models, we are able to identify tools that lead both an expressive family of densities but also tractable normalizing constant via the Feyman representation.
Products of experts. There has been relatively less work products of experts for VI. In part they have been underexplored because it is generally intractable to compute their normalizing constants, and in one way or another these normalizing constants enter into the optimizations for ELBO-basel BBVI. Products of experts (including t-distributed experts [54]) have been used for generative modeling [16, 17]. but in these studies the learning was finessed by computing gradients of the contrastive divergence instead of the KL divergence. Products of experts can also be viewed as special case of energy-base based models. The latter have been used for VI [60], but in this context they are typically paired with Markov chain Monte Carlo (MCMC) methods.
C Latent variable model for products of experts: additional details
C.1 Feynman parameterization of product of t-experts (Result 2.1)
We now provide a full derivation of the Feynman parameterization of the PoE with (-experts. In this case, the denominators have the form
A == 1 + (z- μx) Ax - pk).
(C.1)
Now we substitute the denominators A₁ into the Feynman parameterization in Eq. 7:
II 0=)=)===
II
(C.2)
(2 -4)) An(2 1x) ak
=Ap
r ook) I T(Oz)
I
ak Wk
dw.
(C.3)
EE we[1 + (2 ux) T Ax(2 pkx )2)
JAK1
Now we will rewrite the denominator Eq. C.3 convenient form by expanding the terms and completing the square. In particular, the sum in the denominator can be written
Σ wik [1 + (z (x-1)) AA [[]=)=)]]
(C.4)
Σ Wk + Σ wa(= mm)) An(2 ()x)
(C.5)
=1 + Σ (.
2z
Σ Wk- ARJk))
+ Σ wopt Agtis
(C.6)
where in the last line we have isolated the terms that are quadratic and linear in z. Our next step is to complete the square, noting that
A To = ||2 - A -bl2 -T A
(C.7)
for any invertible, symmetric matrix A and vector Next we identify A Σg w& and b = Σn WxAsst from the quadratic and linear terms in Eq. C.6. Completing the square in this way, we can write
Σ wo [1 (2-1E) TA [k)] = + (o-m(w))1133 + +w?(w),
(C.8)
where we have defined
A(w) Σ AkAt
(C.9)
m(w) = -A-1(w) WeAstE
(C.10)
dww) = =n((w] (u(w)T((w))) A(())((w)
+ Σ
(C.1I)
k
15

---

Next we observe that the expression for (²(w) in Eq. C.11 1.11 can be written more simply as
d²(w) = Σ wk Il (uw)
(C.12)
so that Eq. C.8 reproduces the earlier result for this denominator given by Eq. 12.
Thus, after completing the square and rearranging terms, the product of experts can be written as:
II (=) Bi(a) km
(C.13)
Ex wikk (z u4)) As(2 - p* [x]])= az
1
Dir(a)
(C.14)
[[1 + 0 + u(w))) A(w)(===<(())) =m
Diri(()))
(1 (1+ at
Alw) 1 1 (z p(w))T
(1+0 (z u(w))
(C.15)
1 + (1+ au
1
Dit(a)
(C.16)
Σ ap
1 1 ||
a (
ll
The form in Eq. C.15 is useful because that is proportional to a multivariate t-distribution with location parameter μ(w), inverse scale matrix
identify the quantity inside the integral function
(ww)= = (w)) A( Aww), and degrees of freedom == >x ar- D.
Thus, have shown that the Feynman parameterization allows us to express the product of t-experts as a continuous mixture:
II gh (2) = Din(a) (1 (w) E. an [1 J || ((w))])(w
E. a+
(C.17)
C.2 Computing the normalizing constant (Result 2.2)
The definition of normalizing constant and Eq. C.17 gives
C. = / II (x) = k=1
K
/ (EunDiittoo
1 [|(2 --((u)))) (((w))
-e -- 0k
dz.
(I o² (w))
(C.18)
Because the integrand neganiga,,ve, we can apply Fubini's (Tonelli's 's) theorem, to interchange the integral with the expectation:
[ [1+1 ++11l(= 11( ((w)ll E. ax
C. = EwwDDi
dz
(C.19)
(1 (2(w))2+
Recall that the function in the inner integral is proportional to multivariate t-distribution with location parameter μ(w), inverse scale matrix ?(w) [-o²(w) A(w), and degrees of freedom
// = 2 E. ag - D. Thus, we can identify the integral as being equal to the inverse normalizing constant of this multivariate i-distribution multiplied by (1 + =²())) Σeor:
- ( - P( (1+ (1)))(w))- E. oz = T(Ex an - 2)
-2) (ExaH)
+ ( - -- [xxa)|((w)]+
where we simplified terms after expanding ((w). We also note that E. ausign(| needed to ensure integrability.
16

---

After substituting the definition of v, the normalizing constant can be written
/ II k=1
K
Co
(=
[A(cc)] (1 + o (w))
(C.20)
Thus, this formulation of the normalizing constant turns the problem of integral over TRD into integral over the K - simplex.
D Score-based variational inference
In this section, we first describe heuristic for placing the experts and then we perform an empirical study for this heuristic on an increasingly large number of experts. Then we perform an empirical study for the learning parameters of the variational inference algorithm.
D.1 Heuristies for placing the experts
In this section, we provide one heuristic to address the main intuitive goal outlined in Section 3.1. Here, the goal is to select large number of experts by first locating the mode(s) and then placing large number of experts nearby to help shape the mode(s). The irrelevant experts are pruned away during the fitting the weights via score matching.
will now assume there is one mode, but this discussion can be extended to multiple modes. First, locate the mode 2" by hill climbing with the objective log p(2), where p may be unnormalized. At the mode, we compute or estimate the curvature at 2", which is given by H²2") = this matrix arises from matching the second-orde Taylor expansions of log p and log 9u- The expert placed at this mode
1 V2 log p(=");
2
#i 2", A = -
(D.1)
(D.2)
Now we refine the mode by placing additional experts using the unnormalized target p(z). Using low discrepancy sequence, generate M candidate points from the hypercube
#1 / diag(AA))),u + 8V diag( A
(D.3)
where s 0 scaling parameter. Then we resample the candidate points by sampling N points with replacement according the weights: for a point z; in the candidate set, we compute
w; x exp(3 log p(2))
(D.4)
where 8 E (0, 1) is tempering parameter.
Finally, we refine the candidate set down to the final set of expert locations: beginning with the mode f11. we greedily construct the set of experts by iterating through the candidate set and adding an element μ4 the set of locations if [|1 - mall < т, where T
threshold parameter. The expert inverse scales are then set to Ag [ H ) where we project [- H (px) to the cone of positive semidefinite matrices.
This heuristic serves as a proof of concept and represents one possible way to place experts; future work will focus on more automated expert-selection strategies.
D.2 Empirical study: expert-placing heuristics
In this section, we study the heuristic for algorithm placing on the sinh-arcsinh target described Sec- tion 4 with D 5, = 0.3. and & 0.7. In all experiments, used sensorflow www.rroaooooooobbliieeeee to generate a randomized Halton sequence, which was then transformed appropriately to obtain the desired hypercube.
We began by placing M 50,000 candidate points, where we set 15, 30 0.5, r 6. When choosing the locations, greedily chose locations until hitting the desired number of experts
17

---

# of aclive experts
Fisher
Divergence
10
Ke10
Ke25
K-50
K-75
K=100
Ke10
Ke25
K=50
K=75
K=100
Figure D.I: Study of expert-placing heuristic. Left: KL and Fisher divergences for an increasing number of experts K. Right: The number of active experts after reweighting same sets of experts.
B=10
B=100
= 100
B-1000
= 1/fe 1]
Fither divergence
B 8-00000
Fisher divergance
Q- 1/1 11?
10"
10"
75
100
125
100
175
200
25
75
100
125
150
176
200
Iteration
Iteration
Figure D.2: Study of the learning parameters in score-based VI for the PoE family. Left: varying the number of samples B while fixing the learning rate. Right: varying the learning rate schedule while fixing the batch size.
K = 10, 20,50, 75, 100. We ran the variational inference algorithm for = 200 iterations with B 20,000 and constant learning rate l7. = 1.
In Figure D.I, we report the metrics KL q) and F(g:p). We observe that the KL and Fisher divergences decrease as the number of experts increase. We also report the number of active experts, ie., the experts with non-zero weight. As the number of experts increases, the number of active experts increases slowly; in particular for K = 50, 75, 100, there were about 20-25 active experts. Thus, even with an overspecified number of experts, the method can learn a sparse set of experts.
Computational cost of using experts. As showed in this section, using experts will result in an improved posterior accuracy. But will also increase the computational cost of solving for the weights. Once the experts fixed, determining the PoE parameters comes a relatively small cost, since we need only solve a quadratic program to determine a, which can be done in matter of seconds, solve is O(K3), which is not prohibitive for K < 1000. For the problems studied in this work, we did not require this many experts. For very large K, the cost will become prohibitive, and this issue provides motivation to work on better methods for expert selection in higher dimensions. In addition, the cost of computing scores is linear the dimensionality of the latent variable z, and in many real-world problems, evaluating the unnormalizing target and its score is the dominating bottleneck.
in high dimensions. The computational cost of the constrained least squares
D.3 Empirical study: learning parameters
We continue studying the same target example as Appendix D.2. We used the same expert selection parameters, fixing K 50. Now, we vary the learning parameters: the number of samples B, and the learning rate schedule л)- For the number of samples, we considered B 10, 100, 1000, 10,000; here we fixed Л)t 1. For the learning rate schedule, we fixed B = 10 and considered constant learning rates with = 100 and varying learning rates with 7/t 1/(1 + and 1/(L+ + 1)
In all experiments, initialized a
1. Here we observe that the larger batch sizes had faster and
18

---

more stable convergence. We additionally observe that the algorithm was not sensitive to the choice of 77 considered.
E Proofs: convergence of score-based variational inference for the PoE family
E.1 Statement of convergence theorem
For completeness and ease of reference, restate the main definitions needed to state the theorem. In what follows, we assume that ||-I| denotes the L2 norm, and |l-Il denotes the Mahalanobis norm matrix M.
The empirical Fisher divergence used the objective is defined
(a) = + log (zz logp(z;))1
(E.1)
and is the first term in the variational inference objective function (Fqq. 27). Let V V and H denote the gradient and Hessian of this term with respect to the expert weights a. With these definitions, we can state the following theorem.
Theorem E.1. [Statement of Theorem 3.1] Suppose that O(g: p) in Eq. is minimized by a unique a" E 6, and also that for all 2 there exists some 2 such that Ell V (a))) - 5 and some À >>0 0 such that H[ > A/ almost surely. Then for constant learning rates 7) = v,. the expected error of the iterates satisfies
B3 Elo0 -
a²1 (1=xx)
a²1
(E.2)
E.2 Proof of Theorem 3.1
Proof. The proof consists of several steps.
Step 1. We begin by showing that minimizing the VI objective subject to the constraint a E E is equivalent to projecting an unconstrained solution into E with respect to Mahalanobis norm. To show this, we first recall that the objective can be written as
()) + + l|a a = a G.a + a The
(E.3)
where
G.
1
Σ
Q Q+
11, J5
h4 = -2 m logr(zz) Lo(0)
(E.4)
T
Let (t+++) denote the unconstrained minimizer of this objective, which has the closed-form solution
ā (++1) -
- h hr-
(E.5)
Completing the square of the RHS of Eq. E.3, we have
a T G,a + a T h4 ((+1)(+1)) Gi(a (((+))) (ol-+)))TG,(a------ Thus, we conclude that minimizing the VI objective is equivalent projection onto E of the unconstrained solution in the Mahalanobis norm w.r.t. G;:
G.
(E.6)
arg min (a) + ||a 10 arg min ||ax ā (1+1) ag
12.
(E.7)
ae
Step 2. Our next step is show that bounds the error of the unconstrained minimizer translate bounds on the error of the minimizer in 66 To this end, we define the error and unconstrained error, respectively,
e.1 = a (t+1) -a م④ع == (1+1) - a²
(E.8)
(E.9)
19

---

Next we exploit the nonexpansive property of the projection in the Mahalanobis norm onto the closed convex set 6 [2, Proposition 4.8]. This property implies that the norm of the error ct+1 related to the weighted norm unconstrained enror 1 In particular, noting that 0++1 is the projection of ū and a² the projection of ax itself (since a* E 6), we have that
lec+alla. [|a a ll.. < [|a((+1) -a" a"Ilè, 112+1 Ila.
(E.10)
Thus, it suffices to bound the error of the unconstrained minimizer.
Step 3. Our next step is relate the unconstrained minimizer ((+1) the gradient V at the optimizer o² of p): By doing so, we will make explicit the relationship between the unconstrained solution (1+1) and 2 (a("", which provides notion of misspecification in terms of the scores. Define the Hessian of ((x) to be
H. [[)
T Qo.
(E.11)
Let & denote one-half of the gradient of 9; at a², which given by
B
5 = (V((a") = v
1
V log (0g(((bl)) log pog
B
b=
B
TeQu (V 1lgg((z - logp(zz))
B
Σ
T (Qea" V logp(z))).
(E.12)
Re-arranging Eq. E.12 and noting the definition of He gives
B
B
Σ
Q V logp(z)) = B
Σ
Q Qra" - 6, = H1a" - dc.
(E.13)
B
b=1
Using Eq. E.13 and the definition of h in Eq. E.4, we can now re-express the unconstrained solution in Eq. E.5 function of &c:
(1(+1) = - -2
Q V Vlog((>>) + 1 a()
(E.14)
d=1
= H1a" + La - 6c
(E.15)
m
Step 4. Our next step is to derive a recursion for the error in a Mahalanobis norm. We start by analyzing the error of the unconstrained minimizer. satisfies the following recursion w.r.. e,;
ē ((t+1) -a¹
(E.16)
G 11i" + (1( -6 -a"
(E.17)
n
I - 11 + 1. G-16
(F.18)
a" 1 "G-la" + m
G (t)) G.... 8 a"
(E.19)
nt G. (a aa) G. 5
(F.20)
G. G 5z),
(E.21)
where in the third line we expanded the terms and substituted H, = G, I. Now we can translate
7/t
this recursion on the unconstrained error into for the actual error. To do so, appeal to the
20

---

nonexpansive property of the projection onto 8, From Eq. E.21, we find in this way that
||e+11 lo. - |ēt+1l
ll..
(nonexpansiveness)
(E.22)
|
Gi
1 - 6z)
(using Eq. E.21)
(E.23)
Ga
Ler _t
(E.24)
GiT
Step 5. Our next slep is to translate the bounds the error in the Mahalanobis norm into bounds on the error in the Euclidean norm. In particular, the former can be bounded above and below by using the Euclidean norm and multiplying against cigenvalues of the matrices in the quadratic forms:
minn (G+) |let 21 -
||e
Le, de I
Avn (G+]
N/ Et
d
(E.25)
Consequently, we obtain the following bound on the error 6x+1:
|ee+1 1 S (G.
l Et - 81
(E.26)
2
l (et 71:0a)
(E.27)
Nama (G.)²
|le, n,.6112
(E.28)
Aua (G,)² mi
Now, we note that Amia (Gi) Amia (H₂) + ne
and so Ain (G₁) "m (ania(H)))) (H₂) + 1. Making this
substitution and taking square roots, have
|ec+1ll2 S
m/e [C.) ||e, n76eel2
(E.29)
Htmedam(HH) (H.) |le, n1 112
(E.30)
Tptt llet 3/6el12
(using H1 - AI)
(E.31)
S Tma (1leell2 + (112)
(triangle inequality)
(E.32)
T (lleell2 + (1
(using ?) > 0)
(F.33)
Step 6. Our next step is translate this (per-iteration) recursion for the error in particular run of the algorithm into bounds on the expected error after t iterations. Because each update the iterates (1+1)) depends on the samples {(zb; =b)) p= q(> 0z),, the iterates are a stochastic process. When we take the expectation, we are taking expectation with respect to the filtration defined by {1200
We proceed by taking expectations on both sides of the recursion in Eq. E.33 and using our assumption ||| < 6. In this way find
that LB (a*)))
E leetill2 S -- I+mx
(E.34)
Unrolling this recurrence, obtain the bound
E lleell2 S Pe ||eoll2 + nc.
(E.35)
where
Pt:
I+
(B.36)
j=1
m +
4/ = Σ
(E.37)
1tt
na
Step 7. Our final step is to examine the above bound for a constant step size = 17. Define the quantity
x)n,A)==
ItπX <1
(E.38)
21

---

Relative min eigenvalue of H[Dr)
Relative algenvalue HLDd
1020
100
lleration
Iteration
Figure E.I: Relative minimum eigenvalues with K - 20 (left) and K 100 (right) experts.
In this case have that
Pr (1n; A)
(E.39)
1/ = ΣID ppnn, 170
(E.40)
77 Ej (011, \)
(change of variables = t- -
(E.41)
no (n, A) ((n, a)
(sum of geometric series)
(E.42)
- =pnn, 2)
no 1+nd
[p(1,, A) -s(n, 2)
(E.43)
À p(n, 2)
(E,44)
= [1 - pr]
(E.45)
W. 6 A
(E.46)
and the theorem is proved by substituting the results in Eq. E.39 and Eq. E,46 into Eq. E.35.
E.3 Empirical verification of strong convexity constant À
We ran 100 iterations of Algorithm on the sinh-arcsinh target described in Section 4 with D = Ed 0.3, and d₂ 0.7. We ran the algorithm for increasing number of samples B, In all iterations 0. For each iteration, computed the relative minimum eigenvalue of H[D.]:
Tmin Amin/ Amax:
(E.47)
where higher values indicate
well-behaved positive definite matrix. In Figure E. plot the relative eigenvalues for K = 20 and K = 100. In both cases, we found that larger batch sizes tended to lead more well-behaved matrices. For K 100 and B 10, the matrix was effectively singular, and so we omit that case from the plot.
E.4 Relationship with broader literature
Finally, we highlight the relationship of this convergence result with the broader convergence literature. Theorem 3.1 shares similarities with the convergence behavior of stochastic gradient descent (SGD) applied to strongly convex functions, as established in prior works [13, 39, 40J. However, traditional SGD analyses typically assume fixed objective function throughout the optimization process. In contrast, our setting involves an objective function that changes at each iteration, rendering standard SGD theory inapplicable.
The online learning framework accommodates such scenarios by allowing the objective function to vary over time, even in an adversarial manner [41]. Within this framework, convergence results often focus on bounding the regret, defined as the cumulative difference between the loss incurred by the algorithm at each iteration and the loss of comparator chosen in hindsight. Notably, even in
22

---

9M
0.00
Ressive altecte sample 8ize
0.18
0.12
0.0
3.00
Number of samples
Number samples
(a) Example 1. Skew and heavy tails.
0.00
Snlabrre eftedive sampl 525
030
024
040
0.1
0.12
-2
0.00
Number samples
Number samples
(b) Example 2. High anisotropy.
100
Careblent
0212
25
0.112
25
0.0e
0.04
Eiter
-000
0.008
17
12'
12
-10
samppis
Namber samppis.
(c) Example Many sparsely weighted experts.
Figure F.1: Empirical study for the normalizing constant estimate (second column) and the sampling procedure (third column) for several PoE densities (first column). For each metric, we compute the mean value with respect to 10 random seeds. The shaded region represents the standard error taken with respect to 10 random seeds.
the strongly convex setting, these results provide bounds on regret rather than on the distance the optimal solution.
In contrast, Theorem 3.1 establishes bound the expected & distance between the iterates and the optimal solution a", offering a different perspective on convergence in the presence of time-varying objective functions.
F Experiments: additional studies
In Appendix D.2 and Appendix D.3, we provided empirical studies regarding the expert selection heuristic and the VI algorithm reweight the experts, and in Appendix E.3 provided empirical study to supplement the theoretical convergence results.
In this section, we begin by providing empirical study on the normalizing constant and sampling procedure. In addition, we expand on the details and also provide additional results related the evaluation of the score-based VI method for the PoE family.
23

---

F1 Empirical study: normalizing constant and sampling procedure
In this section, we perform a study for the estimate of the normalizing constant using the Feynman parameterization. We also study the performance of the importance sampling procedure. We considered the following three examples.
Example 1: skew and heavy tails. Consider a product of experts in D = with locations
A1 = [-1,-
-1]" 122 = [0, 0] 13 [1, 1]"
inverse scale matrices
[/3 1/2
1/2 1 .A 1/3
A1
0 1/3 A2
1
and weights a = [1, 11,1,.211"
Example 2: high anisotropy. In our second example, we place 2 experts with locations
M1 [0.0]"
122 [0, O)T
inverse scale matrices
0 1/500
1/500 0
A1
,A
0
1
and weights 2.0].
Example 3: many sparsely weighted experts. Here we considered K 100 experts, but with the majority of a weights on 30 "active" experts. We generated the expert locations M4 unif[-5, 5]. For the expert inverse scale matrices, we randomly generated two eigenvalues Ay. uniformly from [0.05, 2]. and random rotation Rt and then formed the matrix A& R diag(A>>)?? The weight for the a vector was constructed so that the first 30 received weight 0.25, and the remaining received weight 10 12; then the entries were randomly shuffled.
For all three of the above examples, we computed a Monte Carlo estimate of the normalizing constant:
π rI() C. ~ D(
B
Σ
[[(w))] b=1
(1 (w))))
idd Wh ~ Dir(a).
(F.1)
As "ground truth," we computed the estimate of the normalizing constant using 5 million sam- ples. Then, we estimated the normalizing constant using an increasing number of samples B = 10, 100, 500, 1000, 5000, 10,000, 50,000, 100,000, 500,000, 1, 1,000,000. In Figure F.1, we plot the mean absolute error between the B samples and the "true" value. where the mean is taken over 10 random seeds.
For the sampling study, we ran the importance sampling procedure for the same increasing numbers of samples as above. The sampling procedure outputs normalized weights {m) B b=1- Using these
weights, we then calculated (an estimate of) the effective sample size (ESS) [34], which provides measure of how much efficiency in sampling is lost to weight variability:
{ m2
B
ESS
(F.2)
Σb=1 T
In Figure F.1, plot the relative ESS, which divides the ESS by B. The relative ESS measures the efficiency of the sampling procedure, where values closer to are better.
Here we find that in the first two examples, the sampling procedure has high relative ESS, with values above 0.8. For the third example, there is a decrease in efficiency, with values roughly between 0.3-0.5. While this procedure can still be effective for variational inference, it indicates that more samples may need to be taken.
The final result also suggests several opportunities for improvement. First, the sampler could be specialized to take advantage of the sparsity of the weights. These results also suggest that more compact set of experts may lead to more effective variational inference procedures. This motivates developing more advanced approaches for expert selection than the heuristic considered here, which places down a large number of experts and then reweights their importance. However, as we observe in the later VI experiments, given enough samples, the VI fits are still able to provide accurate posterior approximations, even when using large number of experts.
24

---

10"
10
Divergnnca
Floter
Divegencc
Fisher
Divergence
Fisher
16"
10
4012
Badd
Ploer
PdE
ADW
Batd
AD4
PdE.
Rosenbrock,
(b) Sinh-aresing
(c) sunh-aresinh 10
Figure F.2: Additional synthetic target examples
F.2 Target distributions in Section 4
1) Gaussian mixture target: We constructed mixture of Gaussian distributions with weights
0.3, 0.3, 0.3. The means were set μn 1.00, 0.0) M1 [1.0, 0.0) .4g [0.0, 1.0] and the covariances were set to 2 Σ 0.0 0.5, and 2g [.3 1.8.]
2) PoE target: We constructed multimodal PoE target with locations f41 [-2, - and
#2 2, 2] The inverse scale matrices were set to A weights were all set to I.
1.0 0.5 and Ag [ 0.1
0.2 1.0
1.0. The
3) Diamond target: We constructed a PoE target in D=2 with μ1 42= A1
[, 100 9
A2 = = 100 ]_00]__.____ and a= [12,, 1.2]T This target density forms a diamond- shape.
4) Funnel target: The model is 21 22, ... 2D cx)) where set
a² 1.1.
F3
Additional synthetic target results
Additional synthetic target results are in Figure F.2. The Rosenbrock target used was
log ((2) [(1 2 2(22 - =]
(F:3)
The normalizing constant was estimated via quadrature and samples from were generated via a blackjax [5] implementation HMC. The sinh-aresinh targets used were all positively skewed (e- 0.3) and heavy-tailed (77=0.7).
F4 Experimental Setup for Section 4
Computing resources. All experiments were run on CPU, We used Linux workstation with 32-core Intel processor and with 503 GB of memory.
Tuning parameters for baseline methods. The normalizing flow variational family based on a realNVP [10]. We followed the implementation details of the package FlowMC [55]. The base distribution was standard multivariate normal. In all lower-dimensional synthetic experiments (D < 50), we fixed the number of layers to 8 and the number of hidden units to 128. For the 50-dimensional sinh-aresinh target and the posterior-db experiments, we set the number of layers to 16 and the number of hidden units to 128. The flow weights in every coupling layer were initialized randomly using a draw from normal centered
with scale=le-4.
For the Gaussian variational families, the initial values were by default set to a standard Gaussian. For the stochastic gradient- trasient-based methods (ADVI, Flow), optimized the ELBO using Adam [25]. In order to choose the learning rate, we conducted pilot studies where we ran each method for at least 2000 iterations, using grid of learning rate values [0.001, 0.005, 0.01, 0.02, 0.03). The learning rate chosen was the one resulting in the lowest loss value. Some examples are shown in Figure F.3. For BaM, the learning rate schedule was fixed to _BD in all experiments; due the highly Gaussann target problems considered, the performance of this method could likely be improved with additional tuning. For the lower-dimensionall targets (D) < 50), all three methods used the same fixed batch size B 64 and number of iterations = 5000. For the 50-dimensional target, all three methods used batch size of B 128 and number of iterations T 10,000. Thus, allocated the same budget
(1+1]
25

---

0331
1235
C:33
12ze
0000
250
750
1993
1990
152
270
2000
4000
Tonatio
1-1000 0.001
1- 0.005
++001
127
0.50
10-
18
16' 011072
(a) Sinh-arcsinh target example
PoE target example
Figure F3: Examples of pilot studies conducted for hyperparameter grid search. Top: ADVI. Bottom:
Normalizing flow.
Table F.1: Summary of PoE settings for synthetic targets.
Target
K
Active
M
Gaussian mixture
143
50,000
28
0.4
Product of experts
157
22
50,000
28
0.5
Diamond
60
47
50,000
28
0.1
Funnel
50
17
50,000
15
0.5
Sinh-arcsinh (D=5)
100
27
50,000
15
0.5
Sinh-arcsinh (D=10)
480
243
1,000,000
28
0.5
Sinh-arcsinh (D=50)
400
184
3,500,000
28
0.4
Rosenbrock
90
26
3,500,000
20
0.5
in terms of number of gradient evaluations of log (we note that in many real-world problems, evaluating the score of the target dominates the computational cost).
Details PoE VI procedure. The PoE-VI code implemented using Python. To solve the quadratic program in Algorithm 1, we used the qpsolvers package with the solver jaxopt_osqp. which JAX wrapper around the OSQP solver [50]. In order to enforce strict inequality con- straints for integrability, we added small slack parameter inverse scale matrices chosen be positive definite, sufficient condition for integrability is 2 02 sign( (el) > D, but we found that it was often overly strong. Instead we used the constraint 2 > D and verified that the solution resulted an integrable density. By default, we used a batch size of B 10,000 and constant learning rate of 100; the effects of varying these parameters were investigated in earlier experiments. For all experiments, we used the default sampling procedure described in Eq. 25, For the normalizing constant, 500,000 samples were used to estimate the nor- malizing constant, which was used only when computing the estimate of the forward KL divergence. Based on the learning results in Appendix D.2, we set the learning rate to d 1, and ran the algorithm for T 20 iterations. In all experiments,
le-12. We note that when not all
initialized a (0) 1.
In Table F, we summarize the expert selection details used on each of the synthetic targets, along with the number of experts selected and number of active experts in the final PoE.
26

---

posterior-db targets. For all posterior-db experiments, we use the BridgeStan library [44] to first transform each target distribution have support in E. before fitting the VI procedures.
posterior-db: garch11. This model is discrete-time time volatility model with given constant d1. The model uses improper uniform priors on the parameters μ and Ctg > 0, resulting in:
/ uniform
(F.4)
Ogg half-uniform
(F.5)
01 uniform 1)
(F.6)
81 uniform (0, a1).
(F.7)
For each time step = .....,, T::
of = ao + ai =
310
(F.8)
~ (0, 31 = μ + ae E
(F.9)
(F.10)
The PoE method used K = 100 experts, with 24 experts active in the final variational approximation. Experts were chosen with M = 80,000, 12, and 3 = 0.4.
posterior-db: eight-schools. This is classic hierarchical model where each school has its own effect 0y. this model, the global hyperpriors are
~ ((0
(F.11)
T half Salf--Cacchy(O.3 5)
(F.12)
For each school j 1,...,&, 8, the parameters and observations generated as:
πJ (0,
(F.13)
0j ~ μ1 + Th
(F.14)
y/ ~ NV(0,
(F.15)
The PoE method used K 50 experts, with 24 experts active in the final variational approximation. Experts were chosen with M 300,000, s = 14, and = 0.45.
posterior-db: gp-regr. This example is a Gaussian process regression model with mean 0 and covariance function
K K(x,, a² exp (- (x; 104-20 xz2)
(F.16)
2p²
Given N data points, let K denote the N x N gram matrix. The generative model is:
p ~ gamma 4)
(F.17)
a ~ truncated-normal (0, 2²)
(F.18)
truncated-normal (0, 1)
(F.19)
f ~ K)
(F:20)
y ~ [N([ffo')
(F.21)
The PoE method used K 70 experts, with experts active the final variational approximation. Experts were chosen with M = 10,000, 10, and 3 = 0.6.
27

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Yes, we provide two primary contributions: one is introduce new family for variational inference and study its properties, and the other is to propose, analyze, and evaluate an algorithm fit the best variational distribution within this family.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discuss limitations in Section 5.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-specfiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
28

---

Answer: [Yes]
Justification: Yes, we discuss assumptions and provide a sketch of the proof in the main paper, and the full details are in the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We provide details on the experiments in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
29

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [No]
Justification:
We will upload Python implementation online.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cr cc/
sanlliicggiieee/oeeeuutiiii sos.oddeessmmmisioooooooiorc for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: //nips.coo cttps/iccppbl//////u guides/ /CodeSubmi ssi agionyllccy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We provide this information in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification:
The error bars that are included refer to standard error.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
30

---

The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative emrr rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We discuss compute in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: http///eerrr...eeeee
Answer: [Yes]
Justification:
We have reviewed the ethics guidelines.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer:
[Yes]
Justification: have a section on broader impacts in the Appendix.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
31

---

The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: The paper does not pose such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
The data and models used are provided in posteriordb, which we cite in the paper.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
32

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets released.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
We did not use crowdsourcing or human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question:
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
We did not use crowdsourcing or human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
33

---

Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs did not affect the core methodology of this work and were not used in non-standard components.
Guidelines:
The answer NA
that the method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: : htt///eerrss..eee httpe//ris..cc/oooeeeeeeeeeeeeee httpc////o.eeeeee///////////1/// for what should or should be described.
34


Reviews:
[{"review_id": "3T60aNEe5k", "reviewer": "Reviewer_qHA4", "summary": "This work introduces a new family of variational approximations where each member is a weighted product of experts and each expert is proportional to a multivariate t-distribution. Sampling from these distributions is non-trivial and this work presents a Feynman reparameterisation to tackle this issue. This reparameterisation reformulates these PoE distributions into latent variable models with auxiliary Dirichlet random variables. They then propose a Fisher Divergence based objective to fit this distribution to a given target distribution and prove that the corresponding updates converge to the optimal weights for the experts. They appraise this framework's performance by fitting synthetic datasets and datasets from posteriordb in Stan.", "questions": "For most of my primary concerns please refer to the weaknesses section above.\n\nTwo more questions I would like to ask -\n1. What happens when you initialise the mixture components of a GMM at the right spots and only optimise the mixing proportion parameter? In terms of divergences, does the GMM not capture the target distribution accurately? How much better would a PoE be compared to this? When a GMM with a large number of components is initialised for the same target, does standard VI not infer a sparser GMM by driving down mixing proportion values for some of the components?\n2.  From `Figure D.1` it seems that you have to always start with a very high number of experts, off which some are found to be redundant. Does this not lead to an increase in computational cost? Is there any way to be more efficient with this initialisation?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n1. This work, to the best of my novel appears to be novel and original. The work seems inspired and geared towards solving present-day problems.\n2. The theoretical foundations are solidly developed and clearly presented. \n\n**Weaknesses**\n1. Though this work motivates the development of this novel variational approximation from the perspective of having expressive and practical approximate posteriors, I do not see strong supporting evidence in the work to conclude either of these. Firstly, a limited evaluation has been conducted here and without more exhaustive testing on practical applications and uncertainty estimation benchmarks, it's hard to ascertain the quality of the posterior inferred.\n2. The proposed approximation has very elegant foundations but, in the authors own words, optimising the parameters is `complicated` and without a strategy to place the experts at the right spots this approach cannot infer the right expert parameters. In this regard I would also like the authors to comment on the case where experts are placed at completely random locations. This seemingly non-negotiable requirement to correctly initialise the experts, shrinks the applicability of the proposed variational approximation.\n3. According to `Result 2.1` text in `l 113- l 115 `, this PoE can be viewed as a continuous mixture of t-distributions. In this regard I would like the authors to comment on relevant approaches in Implicit variational inference [1][2] where continuous mixtures of Gaussians have been proposed as expressive variational approximations. How does the proposed approach compare to such formulations?\n\n\n**Minor**\n1. `l 25` >Each of such gradient evaluation..... . I think this line is a little vague. In my opinion, it really depends on the downstream task or purpose which determines whether a sample or a gradient is more useful. Kindly consider rephrasing this line or adding more information.\n2.  It's hard to draw conclusions from the figures in the appendix, namely figure E.1 and F.2. In F.2 the legend is not readable.\n\n[1] Unbiased Implicit Variational Inference, Titsias et. al. 2019\n[2] Implicit Variational Inference for high dimensional posteriors, Uppal et. al. 2023", "was_revised": true, "final_justification": "I've upgraded my rating from 3 to 4, please check our discussion below.", "rebuttal": "Thank you for your time and your feedback on our paper. We are glad that you appreciated the “elegant foundations” and novelty of the Feynman representation of the product of experts, the clarity of the writing, and the theoretical grounding of this work, such as the auxiliary variable formulation and the analysis of the algorithm’s convergence.  \n\nBelow, we respond to your feedback and questions on the paper.\n\n**1. Expressiveness of the variational family**  \n\n> “Though this work motivates the development of this novel variational approximation from the perspective of having expressive and practical approximate posteriors, I do not see strong supporting evidence in the work to conclude either of these.”\n\nWe show that the family is expressive in several ways. First, the Feynman representation shows that the product of experts can be represented by a continuous location-scale mixture of t-distributions; indeed with just two experts, Figure 1 shows a variety of behaviors can be captured. Second, we show how the expressivity of the family impacts the convergence of the algorithm: in particular, the iterates converge to a neighborhood of the optimal weighting, where the radius of the neighborhood is determined by the amount of misspecification of the variational family (due to e.g., expert selection). Finally, we show how the PoE family can be used to model a variety of target distributions – we have added several new benchmarks in the rebuttal (see discussion below in #2).\n\n**2. Empirical results on posterior benchmarks**\n\n> “A limited evaluation has been conducted here and without more exhaustive testing on practical applications and uncertainty estimation benchmarks, it's hard to ascertain the quality of the posterior inferred.”\n\nThanks for this feedback. We have now run the procedure on _several new target distributions_, including standard synthetic benchmarks and real-world data benchmarks from posteriordb. We believe that this addition adds to the diversity of targets studied in the empirical results. In particular, we have results covering the following targets:\n\n1. Cross shape\n2. Funnel\n3. Skewed and heavy-tailed\n4. Banana / Rosenbrock **(new)**\n5. Multimodal, well-separated modes **(new)**\n6. Multimodal, modes from mixing **(new)**\n7. PosteriorDB - garch11: heavy-tailed + skewed\n8. PosteriorDB - 8schools: funnel\n9. PosteriorDB - gp-regr: light tailed and symmetric \n10. PosteriorDB - sesame: **(new)**\n\nWe believe this provides a range of target distributions that often arise in hierarchical modeling settings. The initial results for the new targets can be found in our response to Reviewer TPhU, and we plan to include additional posteriorDB models in the final version. \n\n**3. Expert placement at random locations**\n\n> “The proposed approximation has very elegant foundations but, in the authors own words, optimising the parameters is complicated and without a strategy to place the experts at the right spots this approach cannot infer the right expert parameters. In this regard I would also like the authors to comment on the case where experts are placed at completely random locations. This seemingly non-negotiable requirement to correctly initialise the experts, shrinks the applicability of the proposed variational approximation.”\n\nWe understand that you are concerned about the need to place experts at the “right” locations. But it is not the case that there is “ground truth” expert placement and the model fails unless it recovers that ground truth.\n\nA concrete example is when the target is itself a product of t-experts – our expert placing procedure has no prior knowledge of the true locations, nor does it try to “infer” these locations. Instead, we place down a lot of experts and use the flexibility of reweighting to match the scores. We show that our procedure can still recover the target distribution without knowing or finding the true locations. \n\n**4. Relationship with implicit VI approaches**\n\n> “According to Result 2.1 text in l 113- l 115 , this PoE can be viewed as a continuous mixture of t-distributions. In this regard I would like the authors to comment on relevant approaches in Implicit variational inference [1][2] where continuous mixtures of Gaussians have been proposed as expressive variational approximations. How does the proposed approach compare to such formulations?”\n\nThe implicit VI setting is an interesting one, and thanks for pointing us to these references. We will add a discussion of this setting and include additional relevant references to the updated paper.\n\nVia the Feynman representation and Result 2.1, the product of experts we consider can be viewed from the framework of implicit VI with a different kernel (t-distribution) and mixing measure (Dirichlet). In particular, the product of experts has a particular way of combining the auxiliary variables via Equations 15 and 9–10, whereas the typical implicit VI setting specifies a neural network with a mixing distribution over the parameters of that neural network. \n\nIn this special case of continuously parameterized models, we are able to identify tools that lead to both an expressive family of densities but also a tractable normalizing constant via the Feyman representation. Our goal in this work has been to understand the potential and limitations of this special case.\n\n**5. Comparison with mixtures of Gaussians**\n\n> “What happens when you initialise the mixture components of a GMM at the right spots and only optimise the mixing proportion parameter? In terms of divergences, does the GMM not capture the target distribution accurately? How much better would a PoE be compared to this? When a GMM with a large number of components is initialised for the same target, does standard VI not infer a sparser GMM by driving down mixing proportion values for some of the components?”\n\nThe weights in a mixture of experts play a different role than a product of experts. A useful way to distinguish between the two constructions is that mixture weights act like an “or,” whereas product weights act like an “and”. While the latter may seem restrictive, we also don’t require the experts themselves to be normalized, unlike in the mixtures case.  \n\nThus, we expect the two models to excel in different types of applications. To give an extreme and pedagogical example, if the target distribution is a student’s t distribution, we need an infinite mixture of Gaussians to model this. \n\nIn addition, we note there are differences in _optimization_: products of experts are more amenable to score matching, as the objective is linear in the weights, leading to a constrained, nonnegative least squares (NNLS) problem. On the other hand, mixtures do not have efficient solutions for score-based divergences (or even the KL divergence), and are often optimized with gradient descent. \n\nThese differences are important because the NNLS objective used in the PoE setting tends to learn _sparse_ solutions, so the parameters are driven down to 0. It’s unclear whether optimizing the GMM weights (e.g., using \"standard VI\" or ELBO optimization) would lead to a sparse solution. \n\nWe will add a discussion to the revised paper highlighting the conceptual differences between these models.\n\n**5. Computational cost of using many experts**\n\n> “From Figure D.1 it seems that you have to always start with a very high number of experts, off which some are found to be redundant. Does this not lead to an increase in computational cost? Is there any way to be more efficient with this initialisation?”\n\nUsing more experts increases the cost of determining the PoE parameter alpha, but this comes at a _relatively small cost_, since finding alpha boils down to solving a quadratic program, which can be done in a matter of seconds, even in high dimensions.\n\nWe performed timings for solving a QP with linear inequality constraints for different matrix sizes $K$ (note that in this work, we only solve QPs up to size $K=100$):\n\n| K (matrix size) | Mean time (s) |\n|---|-------------|\n| 10                      | 0.0007 |\n| 50                      | 0.0013 |\n| 100                    | 0.0043 |\n| 500                    | 0.0900 |\n| 1000                  | 0.3780 |\n\nHowever, your question about other ways of initializing is an interesting one. We believe there are many ways this can be done, and for this paper we decided to focus on the simplest strategy. For future work, we have been actively investigating other strategies. See our response to Review xP2U for some concrete examples.\n\n**6. Additional clarity points**\n\n> “l 25 >Each of such gradient evaluation..... . I think this line is a little vague. In my opinion, it really depends on the downstream task or purpose which determines whether a sample or a gradient is more useful. Kindly consider rephrasing this line or adding more information.”\n\nThanks for this feedback. Due to space limitations, we cut a contextual footnote related to this after the phrase “considerably more than what is provided”. \n\nThe footnote said: “Here is a simple illustration of this point: suppose that $p$ is univariate Gaussian. Then one can determine the mean and variance of $p$ from the derivatives of $\\log p$ at any two points in $\\mathbb{R}^1$, whereas from $n$ samples (but no scores) one can only estimate these moments to $O(1/\\sqrt{n})$ accuracy.”\n\nWe are happy to add this back to the final revision of the paper.\n\n> “It's hard to draw conclusions from the figures in the appendix, namely figure E.1 and F.2. In F.2 the legend is not readable.”\n\nWe will increase the sizes of the figures and the axis fonts."}, {"review_id": "jcBBYLmS8H", "reviewer": "Reviewer_xP2U", "summary": "This paper explores a new family of models for black box variational inference, where the models are a weighted product of experts, with each expert being a multivariate t-distribution. Noting a relationship to a Feynman identity of integrals from physics, the authors is how how this model can be written as a latent variable model, where the latent variable is Dirichlet distributed. This enable the authors to define a sampling strategy and a means to compute expectation values. From the standpoint of VI, the authors introduce a heuristic scheme for selecting experts, but also show how using a modifier Fisher divergence for the optimization of the weights of experts can lead to a tractable Quadratic program. A convergence theorem is also proven. On several synthetic experiments, the model is seen to work well.", "questions": "Could you move more information for both the heuristics for placing experts, and the details of the experiments into the main body? these two aspects seem under described in the text, and the paper would benefit from more details in the main body.\n\nDo you have a possible path to move away from heuristics for selecting experts, and instead optimize their parameters?\n\nDo you a have a notion of how big of a bias is introduced by the empirical estimate for the loss? Do you have any experiments that show if this significantly affects the optimization?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "This is an interesting paper, introducing a new model and relying on mathematical tools from physics. The paper is clearly written, and includes description of the model, the sampling, the inference strategy, and the convergence.\n\nOne weakness is the method for selection of experts. This is essentially a heuristic method, that requires hill climbing and hand selecting means, inverse scale, and numbers of experts. Without a more automated and optimizable method for selecting experts and estimating their parameters, is is not clear how broadly useful this method can be, especially in anything but low dimensions of the data. Nonetheless, this work provides a good first step towards exploring this family of models.\n\nThe documentation of the experiments is also difficult to assess in the main text, although significantly more useful information is made available in the appendix.", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your constructive feedback on the paper. We are happy to hear that you found the paper interesting, clearly written, and that you appreciated the overall approach. We respond to your comments and questions below.\n\n**1. Expert selection heuristics**\n> “One weakness is the method for selection of experts. This is essentially a heuristic method, that requires hill climbing and hand selecting means, inverse scale, and numbers of experts. Without a more automated and optimizable method for selecting experts and estimating their parameters, is is not clear how broadly useful this method can be, especially in anything but low dimensions of the data. Nonetheless, this work provides a good first step towards exploring this family of models.”\n\nWith regard to the selection of experts, there is a _tradeoff_: if the experts are placed randomly, then we may require very many, but the cost of each added expert is very small; on the other hand, with very carefully selected experts, we may only require a small number, but we will require a great deal more computation to place them effectively.  \n\nWe have explored the first regime in which the expert placement is extremely simple and cheap, but as suggested by the reviewers, we are also looking at other places on this spectrum, where we consider more sophisticated and expensive ways to choose/adapt a smaller pool of experts (see #3 below for some strategies). \n\n**2. Requested additional details in main text**  \n\n> “Could you move more information for both the heuristics for placing experts, and the details of the experiments into the main body? these two aspects seem under described in the text, and the paper would benefit from more details in the main body.”\n\nYes, due to space, we moved many results and details to the experiments. With the extra page for the final revision, we will be able to better describe the expert placing and move more details of the experiments into the main body (e.g., details currently in Appendix F). \n\nIn the main paper, we will revise the experiments section to provide a more complete overview of what to check out in the appendix, so the reader knows where to look, along with a _summary of the additional results in the appendix_, e.g.,: \n* empirical studies on expert placement and varying the number of experts, \n* understanding sensitivity to varying learning hyperparameters, \n* empirical verification of convergence theory, and \n* simulation studies on the normalizing constant estimation procedure. \n\n**3. Approaches beyond random expert placement** \n\n> “Do you have a possible path to move away from heuristics for selecting experts, and instead optimize their parameters?”\n\nYes, there are many interesting directions that we are pursuing but did not make it into the paper for lack of space or being out of scope; here we chose to focus on other elements of the paper (e.g., studying the PoE family and theoretically characterizing convergence) and to use the simplest approach for expert selection. Below, we describe three approaches of increasing complexity.\n\nOne idea is to use a score-based metric for finding where to place experts (in contrast to what we are currently doing, which involves placing a lot of experts randomly and pruning): that is, we check where the scores between the target and variational approximation are a poor fit, and place new experts to patch up the scores. \n\nAnother approach is to select experts in a greedy “boosting” way (e.g., [1]) – that is, we iteratively repeat: 1) adding a new expert according to regions where there is large score mismatch, and then 2) running the described reweighting procedure to refit the weights of the experts. \n\nA final approach is that this family can be directly optimized using gradient descent. In addition, this family could be used in an amortized inference setting by letting the expert parameters and weights be a function of a shared set of parameters (over a data set). While we expect such approaches to utilize more compute, it could be interesting in settings where the user has a large computational budget for training, and would like posterior inference to be fast at test time. \n\nWe are actively investigating all of these directions.\n\n[1] Miller et al. Variational boosting: iteratively refining posterior approximations. ICML, 2017.\n\n**4. Bias from empirical loss estimate**  \n\n> “Do you a have a notion of how big of a bias is introduced by the empirical estimate for the loss? Do you have any experiments that show if this significantly affects the optimization?”\n\nThe bias introduced due to sampling from $q(z| \\alpha^{(t)})$ instead of $q(z| \\alpha)$ does not prevent the iterates from converging to the optimal parameters. In particular, we prove in Theorem 3.1 that the iterates $\\alpha^{(t)}$ obtained from solving Equation 27 converge to the optimal parameters (under no misspecification; otherwise, they converge to a neighborhood of the optimal parameters). We will add a sentence in this section clarifying this point. \n\nFinally, we also show experiments investigating the convergence under different values of $\\eta_t,$ which controls the step size, and we find that the optimization is not very sensitive to this parameter."}, {"review_id": "vG4IwgkcLt", "reviewer": "Reviewer_TPhU", "summary": "This work presents a novel score-based BBVI procedure based on a Product of Experts (PoE) that models skewed and/or heavy tailed distributions more accurately than state-of-the-art methods.\n\nTo that end, the authors introduce a new family of PoEs, each proportional to a multivariate t-distribution. Notably, they show how to draw samples from this (complex) family and how to estimate its normalizing constant, by transforming the problem via the Feynman identity (Section 2).\n\nWith this PoE (plus importance sampling tricks), the authors derive an iterative score-matching procedure to optimize the PoE's geometric weighting of members, via minimization of a regularized Fisher divergence (Section 3). This minimization is proven to be a convex quadratic problem, hence resulting in rates of convergence to the neighborhood of an optimal weighting of experts.\n\nThe empirical experiments in synthetic and a Bayesian posterior benchmark of Section 4 illustrate how BBVI with this PoE better approximates skewed and/or heavy tailed densities.", "questions": "- Can the authors elaborate on the benefits (and if any, drawbacks) of the proposed PoE based BBVI technique beyond the skewed and heavy-tailed distributions presented in Section 4?\n    - Namely, what is the performance of PoE-BBVI in other posteriors?\n    - Is it better/comparable/worse than alternatives?\n\n- Can the authors elaborate on the computational cost/complexity of computing the PoE within the BBVI framework?\n    - what is the tradeoff between compute time and posterior accuracy?\n    - can these metrics be reported for the evaluations in Section 4?\n\n- Even though it is listed as a limitation, more details on the initialization procedure would be appreciated:\n    - How sensitive are the presented results to a good initialization of the fixed collection of experts selected upfront?\n    - What is the number of experts used in the experiments, and how do results change with respect to varying it?\n    - What's the computational complexity of such initialization procedure?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 4, "strengths_and_weaknesses": "**The contributions of this work are two-fold**:\n- The use of the Feynman identity to demonstrate how to draw samples, compute the marginal and estimate quantities of interest from a PoE with experts proportional to multivariate t-distributions (efficiently over the simplex, instead of over the reals)\n- The use of such PoE within BBVI to derive an iterative score-matching procedure for optimization of the geometric weighting of experts for successful approximation to skewed and heavy-tailed distributions\n    \nAdditionally, I would highlight the theoretical analysis of the weighting of experts procedure, resulting in convergence guarantees.\n\nAs the **main weaknesses** of this work, I would identify:\n- The heuristic proposed to select the experts (and the lack of evaluation of its computational cost and impact on performance)\n- The (relatively) complex procedure to attain gains (only?) on skewed and heavy-tailed posteriors\n    - See my questions below on the method's tradeoffs, and whether results in other posteriors are comparable to SOTA", "was_revised": true, "final_justification": "I have increased my score towards acceptance, as the authors' rebuttal has been thorough and informative, clarifying many of the minor concerns raised, as well as by providing new empirical evidence of the benefits of the proposed method.", "rebuttal": "We appreciate your feedback and engagement with the paper. We are glad to hear you found the paper clearly written and highly original. Thanks also for highlighting our _theoretical analysis of the algorithm_ – due to the changing objective function, we could not apply traditional SGD theory for strongly convex functions, and needed to develop a specialized analysis of the algorithm. \n\nBelow, we answer your questions and comment on feedback in the review.\n\n**1. Performance beyond skewed / heavy-tailed posteriors**  \n\n> “Can the authors elaborate on the benefits (and if any, drawbacks) of the proposed PoE based BBVI technique beyond the skewed and heavy-tailed distributions presented in Section 4?\nNamely, what is the performance of PoE-BBVI in other posteriors?\nIs it better/comparable/worse than alternatives?”\n\nFirst, we would like to emphasize that we evaluated more than _just_ skewed and heavy-tailed posteriors – these posteriors were emphasized because some standard baselines do not work well in these settings (e.g., Gaussians cannot model skew or heavy tails, and the standard flow-based models are not designed for heavy tails). \n\nIn our experiments, not only did we study skew and heavy tails, we also studied examples that are _light tailed_ (e.g., garch, funnel, 8-schools, gp-regr) and have _complex geometries_ (e.g., the funnel, the cross, and 8-schools). Furthermore, we found that our PoE method achieved the best resulting Fisher divergence on the posteriordb targets in Figure 4. \n\nThat said, we appreciate your curiosity in understanding the performance in a broader class of models. We have added a few other examples that showcase different target distributions. The _summary of target distributions_ is below, where we have indicated which ones are new:\n1. Cross shape\n2. Funnel\n3. Skewed and heavy-tailed\n4. Banana / Rosenbrock **(new)**\n5. Multimodal, well-separated modes **(new)**\n6. Multimodal, modes from mixing distributions **(new)**\n7. PosteriorDB - garch11: heavy-tailed + skewed\n8. PosteriorDB - 8schools: funnel\n9. PosteriorDB - gp-regr: light tailed and symmetric \n10. PosteriorDB - sesame: **(new)**\n\nWe believe this provides a diverse set of target distributions that often arise in hierarchical modeling settings. For the final revision, we plan to include additional posteriorDB models. \n\n_**Additional results**_\n\n**Rosenbrock function (banana)**\n\n| Approach |   KL  | Fisher |\n|----------|-------|--------|\n| BaM      | 0.89  | 60.66  |\n| ADVI     | 0.60  | 36.71  |\n| Flow     | 0.009 | 0.28  |\n| PoE      | 0.002 | 2.29  |\n\n\n\n**Multimodal example (well-separated modes)**\n\n| Approach |  KL  | Fisher |\n|----------|------|--------|\n| BaM      | 1.59 |  1.14  |\n| ADVI     | 0.93 |  0.64  |\n| Flow     | 0.05 |  0.12  |\n| PoE      | 0.08 |  0.07  |\n\n\n**Multimodal example (modes from mixing distributions)**\n\n| Approach |  KL  | Fisher |\n|----------|------|--------|\n| BaM      | 0.07 |  0.55  |\n| ADVI     | 0.05 |  0.49  |\n| Flow     | 0.03 |  0.43  |\n| PoE      | 0.01 |  0.02  |\n\n\n**PosteriorDB: sesame**\n\n| Approach |  NLL  | Fisher |\n|----------|-------|--------|\n| BaM      | -5.52 |  27.56 |\n| ADVI     | -5.48 | 259.78 |\n| Flow     | -5.45 | 205.76 |\n| PoE      | -5.54 |   2.17 |\n\n\n**2. Computational cost and accuracy trade-off of the VI algorithm**  \n\n> “Can the authors elaborate on the computational cost/complexity of computing the PoE within the BBVI framework?\nwhat is the tradeoff between compute time and posterior accuracy?\ncan these metrics be reported for the evaluations in Section 4?”\n\nHaving more experts will result in an improved posterior accuracy, as we show in the appendix in Figure D.1, but it will also increase the computational cost of solving for the weights. Once the experts are fixed, determining the PoE parameters $\\alpha$ comes at a _relatively small cost_, since we need only solve a quadratic program to determine $\\alpha$, which can be done in a matter of seconds, even in high dimensions. \n\nThe timings of solving a QP with linear inequality constraints are below for different matrix dimensions $K$ (note that in this work, we only solve QPs up to size $K=100$):\n| $K$ | Mean time (s) |\n|---|-------------|\n| 10                      | 0.0007 |\n| 50                      | 0.0013 |\n| 100                    | 0.0043 |\n| 500                    | 0.0900 |\n| 1000                  | 0.3780 |\n\n\nThe computational cost of the constrained least squares solve is $O(K^3)$; we note that this cost is not prohibitive for $K < 1000$, and for the problems in this paper, we did not require even this many randomly placed experts. For very large $K$, the cost will become prohibitive, and this provides motivation to work on better methods for expert selection in higher dimensions. \n\nIn addition, the cost of computing scores is linear in the dimensionality of the latent variable $z$, and in many real-world problems, evaluating the unnormalizing target and its score is the dominating bottleneck. \n\nFinally, we are also happy to report additional metrics for the computational cost. \nWe will revise this section of the paper to discuss the computational cost in more depth.\n\n**3. Expert placement: sensitivity, number of experts, and cost**  \n\n\n > “Even though it is listed as a limitation, more details on the initialization procedure would be appreciated:\nHow sensitive are the presented results to a good initialization of the fixed collection of experts selected upfront?\nWhat is the number of experts used in the experiments, and how do results change with respect to varying it?\nWhat's the computational complexity of such initialization procedure?”\n\n**Sensitivity:** We place experts randomly on a hypercube around each mode, where the bounds are set according to the scale of the target. Thus, if the bounds of the hypercube are set incorrectly, then the tails of the target will not be modeled as accurately. We will add a discussion of this point to the paper. \n\n**Varying the number of experts:** In Appendix F.2, we describe the number of experts used (including the final number of “active” experts).\n\n**How to do results change with the number of experts:** We provide an empirical study in Appendix D.2. We show in Figure 2 the KL and Fisher divergences computed for K=10, 25, 50, 75, and 100. Indeed, the divergences decrease as we add more experts. In addition, we report the number of active experts in each case. \n\n**Computational cost:** With regard to the selection of experts, there is a _tradeoff_: if the experts are placed randomly, then we may require many, but the cost of each added expert is very small; on the other hand, with very carefully selected experts, we may only require a small number, but we will require a great deal more computation to place them effectively. \n\nWe have explored the first regime in which the expert placement is extremely simple and cheap. The main bottleneck of the VI procedure is in the estimation of the weights, not in the expert procedure. \n\nWe will revise the current manuscript to make this point about the cost of expert placement more explicit."}, {"review_id": "tJi7pIWUj4", "reviewer": "Reviewer_UfQj", "summary": "This paper addresses the variational inference (VI) where the approximate variational distribution is assumed to be a weighted product of experts, and each expert is a multivariate t-distribution. To optimize the expert weights, this paper first use Feynman identity to reformulate the product of t-distributions as a latent variable model where the latent variable (denoted as w) lies in the simplex and the observed variable (denoted as z) is conditionally t-distributed given w. Then propose to minimize the Fisher divergence between the PoE and the target density. Some experiments on synthetic and real-world datasets are conducted to demonstrate the performance of the proposed method.", "questions": "See the weaknesses.", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\nOverall, I enjoy reading this paper.\n1. The idea of using the Feynman identity to reformulate the PoE into a continuous mixture of t-distributions is interesting. In this formulation, the mixture weights w are Dirichlet-dsitributed, and the variable z follows t-distribution given w. This makes sampling from the PoE more tractable: one can sample w from a Dirichlet and then sample z from conditional t-distribution.\n\n2. The Fisher div between the PoE and target is optimized using \"batch-and-match\" strategy. Since this requires sampling from PoE for each iteration, the Feynman identity-based reformulation provides a practical advantages by facilitating this sampling process.\n\n3. Experimental results demonstrate that the proposed PoE model captures the heavy tails of the target distributions more effectively than Gaussian approximations, benefiting from the properties of t-distributions.\n\nWeaknesses:\n1. The proposed algorithm only learns the expert weights (\\alpha ), while the set of expert distributions (t-distributions) must be predefined. The parameters of these experts remain fixed during training, which could limit the model’s ability to approximate complex target distributions.\n\n2. The Feynman identity-based reformulation uses the fact that the density of t-distribution is proportional to the inverse of a quadratic from. This property is essential for applying  The Feynman identity, and may not generalize easily to other distribution families. This might limit the broader applicability of the proposed approach.", "was_revised": true, "final_justification": "Considering rebuttal and discussions with authors, I am happy to increase my score.", "rebuttal": "Thank you for your time in reading our paper. We are happy to hear that you enjoyed reading this paper and appreciated the novelty of the Feynman identity approach and the tractability it adds to this PoE framework. \n\nWe respond to your comments and feedback below.\n\n**1. Fixed experts and flexibility** \n\n> “The proposed algorithm only learns the expert weights (\\alpha ), while the set of expert distributions (t-distributions) must be predefined. The parameters of these experts remain fixed during training, which could limit the model’s ability to approximate complex target distributions.”\n\nThere are many types of approaches to forming a flexible density estimator. One way is the “basis set” view, where you attempt to learn the coefficients of a basis expansion. We deliberately adopt a view that is reminiscent of the basis-expansion perspective, where we begin in a large, expressive set of experts that is then reweighted by minimizing the score divergence. Indeed, Figure 1 (leftmost and rightmost panels) shows that with just two fixed experts, changing the weights can result in very different shaped densities.  Furthermore, we show in our experiments that the variational approximation can approximate a variety of target distributions, including challenging geometries such as banana, funnel, cross, and multimodal target distributions (see also new results added, see our response to Review TPhU). \n\nFinally, we note that placing experts randomly and reweighting is one of the simplest first approaches to optimizing this family; there are other more sophisticated approaches that we are actively investigating for future work (see our response to Review xP2U). \n\n**2. Feynman identity and broader applicability**\n> “The Feynman identity-based reformulation uses the fact that the density of t-distribution is proportional to the inverse of a quadratic from. This property is essential for applying The Feynman identity, and may not generalize easily to other distribution families. This might limit the broader applicability of the proposed approach.”\n\nWe note that the Feynman representation of the product of t-experts leads to an already fairly rich variational family: we can model distributions with different tail weights, and we note that the experts themselves don’t need to be normalized. It is an interesting question what other family of functions remain “closed” under convex combinations, are easy to integrate, and lead to an expressive family of densities. We chose to focus on the t-distribution case given its use in the literature [1], and were able to show in Result 2.1 that this already produces an expressive variational family. \n\n[1] Welling et al. Learning sparse topographic representations with products of student-t distributions. NeurIPS, 2002."}]

Author rebuttal:

Meta-review: This paper introduces a neat trick: by leveraging the Feynman parameterization (that allows to re-write a product of denominators into an expectation w.r.t. a Dirichlet helper variable), the authors re-write a product of t-distributions (PoE) into an infinite Dirichlet-mixture over t-distributions. This trick is elegant and new, and shows promising results on a few inference tasks. The paper is very clearly written. All reviewers have voted for accepting the paper and the AC has also read large parts of the paper and fully agrees with the reviewers.

I recommend accepting the paper and assign it a spotlight (disclaimer to the authors: this is the AC's preliminary recommendation which might be changed by the SAC or the PC), as it introduces interesting, theoretical sound and strong method to treat PoEs.

Reasons why not give it a full oral (this was discussed among the reviewers) are that the trick is perhaps a bit niche (not clear whether it can lead to a broader theory) and the experimental evaluation which could have been a bit more ample.

Minor, Typo in the appendix: I think in C.10 it should be \Lambda instead of \Lambda^{-1}. The exposition in the appendix, leading to C.17, is also clearer than in the main paper. The authors should perhaps adopt it into the main paper in the revised version.
Decision: Accept (spotlight)