EigenVI: score-based variational inference with orthogonal function expansions
Authors: Diana Cai
Chirag Modi
Charles Margossian
Robert M. Gower
David Blei
Lawrence K. Saul
Venue: NeurIPS 2024 spotlight (2024)
Keywords: variational inference
black-box variational inference
Bayesian inference
probabilistic modeling
score-based divergence
score matching
non-Gaussian variational families
Primary area: probabilistic_methods
TL;DR: EigenVI uses score matching with flexible yet tractable variational families built using orthogonal function expansions and reduces to an eigenvalue problem.

Abstract: We develop EigenVI, an eigenvalue-based approach for black-box variational inference (BBVI). EigenVI constructs its variational approximations from orthogonal function expansions. For distributions over $\mathbb{R}^D$, the lowest order term in these expansions provides a Gaussian variational approximation, while higher-order terms provide a systematic way to model non-Gaussianity. These approximations are flexible enough to model complex distributions (multimodal, asymmetric), but they are simple enough that one can calculate their low-order moments and draw samples from them. EigenVI can also model other types of random variables (e.g., nonnegative, bounded) by constructing variational approximations from different families of orthogonal functions. Within these families, EigenVI computes the variational approximation that best matches the score function of the target distribution by minimizing a stochastic estimate of the Fisher divergence. Notably, this optimization reduces to solving a minimum eigenvalue problem, so that EigenVI effectively sidesteps the iterative gradient-based optimizations that are required for many other BBVI algorithms. (Gradient-based methods can be sensitive to learning rates, termination criteria, and other tunable hyperparameters.) We use EigenVI to approximate a variety of target distributions, including a benchmark suite of Bayesian models from posteriordb. On these distributions, we find that EigenVI is more accurate than existing methods for Gaussian BBVI.

Eigen VI: score-based variational inference with orthogonal function expansions
Diana Cai
Chirag Modi Flatiron Institute
Flatiron Institute
acal@flaatironhhstttuu....g
cmodi mmdd:fflaairoonniitute.ooom
Charles C. Margossian
Robert M. Gower
Flatiron Institute
Flatiron Institute
cnargpsssian@fiaairroonnttttttte
rgower@ffaairroiisssituee.. org
David M. Blei
Lawrence K. Saul
Columbia University
Flatiron Institute
david. ble elei@columbia.co
Isauloflatirrmmissiitte.
Abstract
We develop EigenVI, an eigenvalue-based approach for black-box variational infer- ence (BBVI). Eigen'II constructs its variational approximations from orthogonal function expansions. For distributions over D, the lowest order term in these expansions provides a Gaussian variational approximation, while higher-order terms provide systematic way to model non- Non-Gaussianity. These approximations are flexible enough model complex distributions (multimodal, asymmetric), but they are simple enough that one can calculate their low-order moments and draw samples from them. Eigen VI can also model other types of random variables (e.g.. nonnegative, bounded) by constructing variational approximations from different families of orthogonal functions. Within these families, Eigen computes the variational approximation that best matches the score function of the target distri- bution by minimizing a stochastic estimate of the Fisher divergence. Notably, this optimization reduces to solving minimum eigenvalue problem, so that Eigen VI effectively sidesteps the iterative mundient-based optimizations that are required for many other BBVI algorithms. (Gradient-baat based methods be sensitive learning rates, termination criteria, and other tunable hyperparameters.) We use EigenV VI approximate variety of target distributions, including benchmark suite of Bayesian models from posteriors torter.rrrbb..... On these distributions, we find that Eigen VI is more accurate than existing methods for Gaussian BBVI.
Introduction
Probabilistic modeling is cornerstone of modern data analysis, uncertainty quantification, and decision making. key challenge of probabilistic inference is computing target distribution of interest; for instance, in Bayesian modeling, the goal is to compute posterior distribution, which is often intractable. Variational inference (VI) [5, 21, 45] is popular method for scalable probabilistic inference that has worked across range of applications. The idea behind VI is to approximate the target distribution by the closest member of some tractable family.
One major focus of research is to develop black-box algorithms for variational inference [6, 15, 23, 27, 32, 37, 40, 44, 46]. Algorithms for black-box variational inference (BBVI) can be used to
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

approximate any target distribution that is differentiable and computable up to some multiplicative (normalizing) constant; as such. they are extremely flexible. These algorithms have been widely implemented in popular probabilistic programming languages, and they are part of the modern toolbox for practitioners in computational statistics and data analysis [1, 4, 7, 13, 43].
Traditionally, the variational approximations in BBVI are optimized by minimizing the Kullback- Leibler (KL) divergence between the variational family and the target (equivalently, maximizing the ELBO). This strategy is powerful and scalable, but it relies on stochastic gradient descent (SGD), which can be difficult to tune [10, 11, 51]. These difficulties can be acute even for Gaussian variational approximations [27, 40], particularly if these approximations employ full covariance matrices.
More recently, researchers have proposed algorithms for Gaussian BBVI that do not require the use of SGD [6, 37J. Instead of minimizing the KI. divergence, these methods aim to match the scores, or the gradients of the log densities, between the variational distribution and the target density. These methods exploit the special form of Gaussian distributions to derive closed-form proximal point updates for score-matching. These updates are as inexpensive SGD, but not as brittle. They show that score-based BBVI can be applied in an elegant way to Gaussian variational families.
In this paper, we show that score-based BBVI also yields simple, closed-form updates for much broader family of variational approximations. Specifically, we propose a new class of variational families constructed from orthogonal function expansions and inspired by solutions to the Schrödinger equation in quantum mechanics. These families are expressive enough to parameterize a wide range of target distributions; the same time, the distributions in these families are sufficiently tractable that one can calculate low-order moments and draw samples from them. In this paper, we mostly orthogonal function expansions to construct distributions supported on RR,; in this case, the lowest-order term in the expansion is sufficient to model Gaussian behavior, while higher-order terms account for increasing amounts of non-Gaussianity. More generally, we also show how different basis sets of orthogonal functions can be used to construct variational families over other spaces.
To optimize over variational family from this class, minimize an estimate of the Fisher divergence, which measures the scores of the variational distribution against those of the target distribution. We show that this optimization reduces to minimum eigenvalue problem, thus avoiding the need for gradient-based methods. For this reason, we call our approach EigenVI.
We study Vigen with variational family constructed from weighted Hermite polynomials. We first demonstrate the expressiveness of this family on variety of multimodal, asymmetric, and heavy- tailed distributions. We then use Eigen VI to approximate diverse collection of non-Gaussian target distributions from posteriordb [35], a benchmark suite of Bayesian hierarchical models. On these problems, Eigen VI provides more accurate posterior approximations than leading implementations of Gaussian BB VI based on KL minimization and soore-matching.
The organization of this paper is as follows. In Section 2 we introduce the variational families that arise from orthgonal function expansions. and we show how score-matching these families reduces to an eigenvalue problem. In Section 3 we review the literature related to EigenVI. In Section 4, we evaluate Eigen' VI on a variety of synthetic and real-data targets. Finally, in Section 5, we discuss limitations and future work.
2 Score-based variational inference with orthogonal function expansions
In this section we use orthogonal function expansions to develop new variational families for ap- proximate probabilistic inference. In Section 2.1, we review the basic properties of these expansions. In Section 2.2, we introduce a acore-based divergence for VI with these families; notably, for this divergence, the optimization for VI reduces to an eigenvalue problem. Finally in Section 2.3, we consider how to use these variational approximations for unstandardized distributions; in these settings we must carefully manage the trade-off between expressiveness and computational cost.
2.1 Orthogonal function expansions
Let 2 C D denote the support of the target distribution p. Suppose that there exists a complete set of orthonormal basis functions (=)) on this set. By complete, we mean that any sufficiently well- behaved function f Z IR can be approximated, to arbitrary accuracy, by particular weighted

---

Table I: Examples of orthogonal function expansions in dimension. The basis functions in the table are not normalized, but they can be rescaled so that their squares integrate to one.
support
orthogonal family
basis functions (x(-)
* E [-1.1]
Legendre polynomials
{1, 2, 322 -3,,
2 e St
Fourier basis
(1,00s 0, sin cos 20, sin .0.....
& E 0, 00)
weighted Laguerre polynomials
2,, -4 e +
E IR
weighted Hermite polynomials
sum of these basis functions, and by orthonormal, we mean that the basis functions satisfy
/ (x))
k
(1)
otherwise,
where the integral is over Z. Define the K -order variational family 2k to be the containing all distributions of the form
K where Σ
q(2)
Σ ag (x))(2))
a 1,
k=
and where a₁ all distributions that can be obtained by taking weighted sums of the first basis functions and then squaring the result.
R for k = 1,..., K are the parameters of the family QK. In words, Qk contains
Eq. involves a squaring operation, a sum-of-squares constraint, and a weighted sum, The squaring operation ensures that the density functions in Qx are nonnegative (i.e., with g(z) 20 0 for all 2), while the sum-of-squares constraint ensures that they are normalized:
K jule)d= [
andb(())
d2 E kAA==1
Σ
(3)
The weighted sum in Eq. 2 bears a superficial similarity to mixture model, but note that neither the basis functions (=)(z) the weights on in Eq. 2 are constrained to be nonnegative. Distributions of this form arise naturally in physics from the quantum -echttnica wave functions that satisfy Schrödinger's 's equation [16]. In that setting, though, is typical consider complex-vall -valued weights and basis functions, whereas here we only consider real-valued ones.
The simplest examples of orthogonal function expansions arise in dimension. For example, functions on the interval [-1, can be represented as weighted sums of Legendre polynomials, while functions the unit circle can be represented by Fourier series of sines and cosines; see Table 1. Distributions on unbounded intervals can also be represented in this way. On the real line, for example, may consider approximations of the form in Eq. where
Φk+1 (z) (V2rk)
(c 1=²
H (*))),
and Hz(2) are the probabilist' 's Hermite polynomials given by
de de
H²(z) = (-1)*
(5)
Note how the lowest-order basis function d1 (z) in this family gives rise (upon squaring) Gaussian distribution with zero mean and unit variance.
Figure shows how various multimodal distributions with dieeennnnio support can be approxi- mated by computing weighted sums of basis functions and squaring their result. We emphasize that the more basis functions in the sum, the better the approximation.
Orthogonal function expansions in one dimension are also important because their Cartesian products can be used to generate orthogonal function expansions in higher dimensions. For example, we can
3

---

K-8
larget
K=1
K=10
K=3
0.5
0.5
(a) Legendre polynomial expansion
(b) Fourier series expansion
(c) Hermite polynomial expansion
Figure I:
Target probability distributions (black dashed curves) on the interval [-1, 1] (left), the unit circle (middle), and the real line (right), and their approximations by orthogonal function expansions from different families and of different orders; see Eq. 2 and Table 1.
approximate distributions over (say) R by
2
Kg K2 ΣΣΣ
((21, 22; 23)
Bija dee
where
Σ
EBg. 1, (6) kgk
where Byz IR now parametrize the family. Note that there are total K₃ K3 parameters in the above expamsion, so that this method of Cartesian products does not scale well to high dimensions multiple basis functions random variables of mixed type: for example, from Table 1, we can create a variational family of distributions over IR x 1-1, 1]x [0, oo) from the Cartesian product of orthogonal function expansions involving Hermite, Legendre, and Laguerre polynomials.
used per dimension. Note that the same strategy
also be used for
As shown in Figure 1, the approximating distributions from T -order expansions can model the presence of multiple modes as well as many types asymmetry, and this expressiveness also extends to higher dimensions. Nevertheless, it remains tractable to sample from these distributions and even to calculate (analytically) their low-order moments, as we show in Appendices and
For concreteness, consider the distribution over R3 in Eq. 6. The marginal distribution ((21) is
((21)= = / q(2,, 22, 23) dzg d23 =
Σ 13'
Σ gk
BegkBung (1()))r((a),
and from this expression, moments such as E[21] and Var[21] can be calculated by evaluating integrals involving the elementary functions in Table 1. (In practice, these integrals are further simplified by recursion relations that relate basis functions of different orders; we demonstrate how to compute the first two moments for the normalized Hermite family in Eqs. B. and B.22.)
To generate samples (2(*),, each dimension is sampled as follows: we draw 2 ((2) by comput- ing the cumulative distribution function (CDF) of this marginal distribution and then numerically inverting this CDF. Finally, extending these idcas, we can calculate higher-order moments and obtain joint samples via the nested draws
20 4(2)), 10 ~ (1)) 23
g æ1, 22).
(8)
The overall complexity of these procedures scales no worse than quadratically in the number of basis functions in the expansion. These extensions are discussed further in Appendices A and B.
2.2 VI
In variational inference, we posit parameterized family of approximating distributions and then compute the particular approximation in this family that is closest to a target distribution of interest. Eq. 2 constructs a variational family Qx from the orthogonal functions (cx(z) KK whose variational parameters are the weights 1- We now derive Figen VI, a method to find q Qx that is close to the target distribution p(z).

---

We first define the measure of closeness that we will minimize. Eigen VI measures the quality of an approximate density by the Fisher divergence [18],
(g.
=
V log ((2) v log p(>))l g())d2,
where V q() and V p(z) are the score functions of the variational approximation and target, respectively. Suppose that and have the same support; then the Fisher divergence vanishes if and only if the scores of q and p are everywhere equal.
Though p is, by assumption, intractable to compute, in many applications it is possible to efficiently compute the score V logg p at any point 2 E Z. For example, Bayesian models the score of the larget posterior is equal to the gradient of the log joint. This observation is the main motivation for score-based methods in probabilistic modeling [6. 31, 37, 48].
Here we seek the E QK that minimizes ((q.p), But now a challenge arises: it is generally difficult to evaluate the integral for (gg. p) in Eq. 9, let alone to minimize it function of q. While it is possible to sample from the distribution q. it is not straightforward to simultaneously sample from q and optimize over the variational parameters {ax) k= in terms of which it is defined. Instead, we construct an unbiased estimator of ((g.p) p) by importance sampling, which also decouples the sampling distribution from the optimization. Let {2 2² zB) denote batch of B samples drawn from some proposal distribution on Z. From these samples we can form the unbiased estimator
B
q( π (2))
(g.g,) =
log g(2%) - logp(=")1²²
(10)
This estimator should be accurate for appropriately broad proposal distributions and for sufficiently large batch sizes. We can therefore attempt to minimize Eq. 10 in place of Eq. 9.
Now show that the minimization of Eq. 10 q E Qx simplifies to minimum eigenvalue problem for the weights {ax} 1
To obtain the eigenvalue problem, we substitute the orthogonal
function expansion Eq. 2 into Eq. 10 for the unbiased estimator of (q, p). As intermediate
step, we differentiate Eq. 2 to obtain the scores
V log 9(22) 2 2I, an (2)(
(1I)
Ee (x²)
Further substitution the scores provides the key result behind our approach: the unbiased estimator in Eq. 10 is simple quadratic form in the weights a == [a1....,,x'] ax] of the orthogonal function expansion,
(g. = a TMa
(12)
where the coefficients of the quadratic form are given by
B
M₁s Σ
T(>b) [2V))(2² Vlog"\Vo [2Von(=")
1
du logg"))
(13)
Note that the elements of the K x K symmetric matrix M capture all of the dependence the batch of samples {2ah B 1b=1. the scores of p and q at these samples, and the choice of the family of orthogonal functions. Next minimize the quadratic form in Eq. 12 subject to the sum-of-squares constraint E. o% = in Eq. 2. In this way we obtain the eigenvalue problem [8]
min V
(.(a.p)
min [ ||l||
Ma] =: Amin (M),
(14)
where Amin (M) is the minimal eigenvalue M, and the optimal weights are given (up to an arbitrary sign) by its corresponding eigenvector; see Appendix C for a proof. Eigen solves Eq. 14.
We note that the eigenvalue problem in EigenVI arises from the curious alignment of three particular maniccs-nammly... (i) the choice of variational family (based on orthogonal function expansions), (ii) the choice of divergence (based on scote-matching). and (iii) the choice of estimator for the divergence (based on importance sampling). The simplicity of this eigenvalue problem stands in contrast to the many heuristies of gradient-based pnniiizzatios-iivolviiii learning rates, terminating criteria, and perhaps other algorithmic hyperparameters that typically required for ELBO-based
5

---

BBVI [10, 11]. But EigenVT VI is also not entirely free of heuristics; to compute the estimator in Eq. 10 we must also specify the proposal distribution x and the number of samples B; see Appendix D for discussion.
The size of the eigenvalue problem in Eigen is equal to the number of basis functions K in the orthogonal function expansion of Eq. 2. The eigenvalue problem also generalizes to orthogonal function expansions that are formed from Cartesian products of one-dimensional families, but in this case, if multiple basis functions are used per dimension, then the overall basis size grows exponentially in the dimensionality. Thus, for example, the eigenvalue problem would be of size K2 K3 for the approximation in Eq. 6, as can be seen by "flattening" the tensor of weights 3 in Eq. 6 into the vector of weights a vec(.5) in Eq. 2. Finally, we note that EigenVI only needs compute the minimal eigenvector of M in Eq. 14, and therefore it can benefit from specialized routines that are much less expensive than a full diagonalization.
2.3 Eigen VI in RD: the Hermite family and standardization
We now discuss the specific case of EigenVI for 2 = TR with the Hermite-based variational family in Eq. 4. For this case, we propose transformation of the domain that serves to precondition or standardize the target distribution before applying EigenVI. While this standardization is not required to use Eigen VI, it helps reduce the number of basis functions needed approximate the target, leading to more computationally efficient procedure. It also suggests natural default choices for the proposal distribution π in Eq. 10.
Recall that the eigenvalue problem grows linearly in size with the number of basis functions. Before applying Eigen VI, our goal is therefore transform the domain in way that reduces the number of basis functions needed for good approximation. To meet this goal for distributions over RP, we observe that the lowest-order basis function of the Hermite family in Eq. yields (upon squaring) standard multivariate Gaussian, with zero mean and unit covariance. Intuitively, we might expect the approximation of Eigen to require fewer basis functions the statistics of the target distribution nearly match those of this lowest-order basis function. The goal of standardization is to achieve this match, to whatever extent possible, by a suitable transformation of the underlying domain. Having done so, Eigen VI in R2 can then be viewed as a systematic framework to model Gon-Gaassi effects via a small number of higher-order terms in its orthogonal function expansion.
Concretely, we consider a linear transformation of the domain:
Σ -
(15)
where ja and Σ are estimates of the mean and covariance obtained from some other algorithm (e.g.. a Laplace approximation, Gaussian variational inference, Monte Carlo, or domain domain-specific knowledge). We then apply the VI to fit a K Kth-order variational approximation (()) to the target distribution p(z) that is induced by this transformation; afterwards, we reverse the change-of-variables to obtain the final approximation to (.
((2) = ((3)
(16)
Figure 2 shows why it is more difficult to approximate distributions that are badly centered or poorly scaled. The left panel shows the effect of translating a standard Gaussian away from the origin and shrinking its variance; note how a comparable approximation to the uncentered Gaussian now requires a 16th-order expansion. On the other hand, after standardization, the target can be perfectly fitted by the base distribution in the orthogonal family of reweighted Hermite polynomials. The right panel shows the similar effect of translating the mixture distribution in Figure I (right panel); comparing these panels, we see that twice as many basis functions (K- 14 versus K = 7) are required to provide a comparable fit of the uncentered mixture.
Finally, we note another benefit of standardizing the target before fitting Eigen when the target has nearly zero mean and unit covariance, becomes simpler to identify natural choices for the proposal distribution Intuitively, in this case, we want proposal distribution that has the same mean but heavier tails than a standard Gaussian. In our experiments, we use two types of centered proposal Intributions-unifln and isotropic Gaussian- whose variances are greater than one.

---

K-1
K-1
K-4
K=9
K=16
K=14
0.5
(a) Gaussian target, mean and variance (b) Mixture target (translation of Figure lc)
Figure 2: Higher-order expansions may be required to approximate target distributions (black) that are not standardized. Left: approximation of a non-standardized Gaussian. Right: approximation of the mixture distribution in Figure I after translating its largest modes away from the origin.
3 Related work
Several recent works have considered BBVI methods based on score-matching. These methods take particularly simple form for Gaussian variational families [6, 37]. The Fisher divergence [18] has been previously studied as divergence for variational inference [47]. and Zhang [48] propose minimizing Fisher divergence for semi-implicit (non-Gaussian) variational families; the divergence minimized using gradient- based optimization. In another line of work, Zhang al. [50] consider variational families of energy- based models and derive closed-form solution to minimize the Fisher divergence in this setting.
More generally, there have many studies of VI with non-Gaussian variational families. One common extension is to consider families of mixture models [14, 17, 36]; these are typically optimized via ELBO maximization. BBVI algorithms have also been derived for more expressive variational
families of energy- based models [9, 22, 28, 29, 52, 53] and normalizing flows [3, 24, 25, 34. 39, 41]. However the performance of these models, especially the normalizing flows, is often sensitive to the hyperparameters of the flow architecture and optimizer, as well as the parameters of the base distribution [2, 11]. Other aspects of these variational approximations are also less straightforward; for example, one cannot compute their low-order moments, and one cannot casily evaluate or draw samples from the densities of energy- based models.
The variational approximation in Eigen is based on the idea of squaring a weighted sum of basis functions. Probability distributions of this form arise most famously in quantum mechanics [16]. This idea has also been used to model distributions machine learning, though not quite in the way proposed here. Novikov et al. [38] propose a tensor train-based model for density estimation, but they do not consider orthogonal basis sets. Similarly, Loconte et [33] obtain distributions by squaring mixture model with negative weights, and they study this model in conjunction with probabilistic circuits. By contrast in this work, we consider this idea the context of variational inference, and we focus specifically the use of orthogonal function expansions, which have many simplifying properties; additionally, the specific objective we optimize leads to a minimum eigenvalue problem.
4 Experiments
evaluate Eigen on synthetic targets and real data targets. these experiments we use orthogonal family induced by normalized Hermite polynomiuls (see Table 1), whose lowest-order expansion is Gaussian. Thus, this variational family can model non-Gaussian behavior with the higher-order functions in its basis. We first study 2D synthetic targets and use them to demonstrate the expressiveness these higher-order expansions. Next, we experiment with target distributions where systematically vary the tail heaviness and amount of skew. Finally, apply EigenVI to set of hierarchical Bayesian models from real-world applications and benchmark its performance against other Gaussian BBVI algorithms.
4.1 2D synthetic targets
We first demonstrate how higher-order expansions of the variational family yield more accurate approximations range of 2D non-Gaussian target distributions (Figure 3); see Appendix E.2
7

---

larget
Geuesian:
K=20,
-L-1122-02
tanget
KL-1.29
KL=0.29
KL=1195-00
larget
KL-2.72
K1.0112
KL 776-02
-2.36-07
Figure 3: 2D target functions (column I): 3-component Gaussian mixture distribution (row 1), funnel distribution (row 2), and a cross distribution (row 3). We report the KL(p: q) for the resulting optimal variational distributions obtained using score-based VI with a Gaussian variational family (column 2) and the EigenVI variational family (columns
where K K2
for details. We report estimate of KL(p; q) above each variational approximation. The Gaussian variational approximation is fit using batch and match VI [6]. which minimizes a score-based divergence. For EigenVI, the target distributions were standardized before fitting Eigen (we compare the costs of the methods in Figure E. E),, and the total number of basis functions is K = K K2.
4.2 Non-Caussianity: varying skew and tails in the sinh-arcsinh distribution
We now consider the sinh-aresinh normal distribution [19, 20], which is induced by transforming a multivariate Gaussian using parameters that control the amount of skew and the weight of the tails. We construct several targets (D = 2, 5) of increasing amounts of -ansaaassanntty in the skew or the tails of the distribution, and refer these targets as slight skew and tails,
skew and tails, and slight skew and heavier tails; see Appendix E.3 for details. In Figure 4a, we visualize the 2D targets and the Eigen'II fits along with their forward KLs. Before applying Eigen we standardize the target using mean and covariance estimated from batch and match VI [6]. In Figure 4b, we measure the EigenVI VI forward KL under varying numbers of samples B and across increasing numbers of basis functions, given by K = Ka. We also present forward match VI (BaM) and automatic differentiation VI (ADVI), which both use Gaussian variational familics and are run using the same budget in terms of number gradient evaluations. Next we consider similar targets with D 5, which are visualized in in Figure E.2, along with the resulting Eigen VI variational approximations. In Figure 4c, we observe greater differences the number of importance samples needed lead to good approximations, especially the number of basis functions increase.
resulting from batch and
4.3 Hierarchical modeling benchmarks from posteriordb
We now evaluate Eigen VI on set of hierarchical Bayesian models [7, 75,3,,22 42), which are summarized Table E.I. The goal is posterior inference: given data observations T1.N. the posterior of = is
p (1:N) x p(=)p(Z1 (x) =: p(z),
(17)
where p(z) the prior and p denotes the likelihood.
We compare VI to 1) automatic differentation VI (ADVI) [27], which maximizes the ELBO over a full-covariance Gaussian family (ADVI), 2) Gaussian score matching (GSM) [37], score-based
BBVI approach with full-covariance Gaussian family, and 3) batch and match VI (BaM) [6], which

---

night slow and take
more skew and balls
sight slmu and beader tide
11-1155-05 K1-8.31-04 Ki-2.00-45
(a) Example 2D targets (left) varying the skew or tail weight 7 components and their Eigen VI fits (right).
slight skew and tails
more skew and talls
slight scew and heavler tails
6.8
ADM
0:11
0.12
Ball
Farverd
2.19
EigenVil
5-1000
Forwerd
0.00
21
15 333
3n basis functions
basis functions
(b) D 2
slight
more
tails
0353
Eigeritf,
Eigenvi, B-20000
Elgervi, B-10000
Forward KL
Forward KL
0.75
BaM
0:53 ADMI
535 * of basis functions
126 # of basis functions
of basis functions
(c) D
Figure 4: Sinh-aresinh normal distribution synthetic target. Panel (a) shows the three targets consider 2D, and their resulting EigenVI fit. Panel (b) shows measures Kl-(p: q) for D = 2, and panel (c) shows KL(p: q) for the r-axis shows the number of basis functions, K =II Ks-
minimizes a regularized score-based divergence over a full-covariance Gaussian family. In these examples, standardize the target using either GSM or BaM before applying EigenVI.
In these models, we do not have access to the target distribution, p(z æ1:N), only the unnormalized target p. Thus, cannot evaluate an estimate of the forward KL. Instead, evaluate the fidelity of the fitted variational distributions, we compute the empirical Fisher divergence using reference samples from the posterior obtained via Hamiltonian Monte Carlo (HMC):
g
Σ
V log p(=") V log ( = >0>p(> p(2|<1/N)..
(18)
Note that this measure is not the objective that Eigen minimizes; it is analogous to the forward KL divergence, as the expectation is taken with respect to p. report the results Figure 5, computing the Fisher divergence for Eigen VI with increasing numbers of basis functions. We typically found that with more basis functions, the scores becomes closer to that of the target.
Finally, we provide qualitative comparison with real-NVP normalizing flows (NFs) [12], flexible variational family that fit by minimizing the reverse KL. We found that after tuning the batch-size and learning rate, NFs generally had a suitable fit. We visualize the posterior marginals for subset of dimensions from 8schools the top three rows, comparing Eigen VI, the NF, and BaM. Here, fit the tails of this target distribution. On the other hand,
we observe that the Gaussian struggles Eigen provides a competitive fit to the normalizing flow. Appendix E.4, we show full corner plot in Figure E.3 and marginals of the garch11 model in Figure E.4.
5 Discussion of limitations and future work
In this work, we introduced EIgen new approach for score-based variational inference based on orthgonal function expansions. The score-based objective for Eigen VI is minimized by solving an eigenvalue problem, and thus this framework provides an alternative to gradient-based methods for BBVL Importantly, many computations in Eigen VI can be parallelized with respect to the batch of

---

0.2
2.5
0.0
-2.5
2.5
-5
(a) EigenVI with normalized Hermite polynomial family
0.4
0.2
-25 0.0 25
-25 2.5
-2.5 0.0 2.5
-2.5 0.0 2.5
-5
(b) VI with normalizing flow family
0.4
0.2
-2.5 0.0
2.5
-2.5
0.0
2.5
-25
0.0
2.5
-2.5
2.5
-5
(c) Batch and match VI with full covariance Gaussian family
Riher divergence
ADVI-G
Risher divergence
302
Risher divergence
Rihee divergence
204
107
GSM
301
1.0
Bal
EigenVI
207
10a
10:4
100
10
10-
of basis functions
basis functions
of basis functions
# basis functions
(d) kidscore,
sesame, D
(f) gegr.
(g) logearn, D 4
Fisher divergence
10
Fisher divergence
104
Fisher divergence
Fisher divergence
107
x 1x
10
107
10
12
10a
101 "0o "10
of basis functions
of basis functions
of basis functions
basis functions
(h) garch11,
(i) arK-arK, D =
(i) logmesquite, = (k) 8-schools, D 10
Figure 5: Results on posteriordb models. Top three rows: marginal distributions of the even dimensions from 8-schools. Reference samples from HMC outlined in gray, and the VI samples are in green. Bottom two rows: evaluation of methods with the (forward) Fisher divergence. The rraxis shows the number of basis functions, K K Kg. Shaded regions represent standard errors computed with respect to random seeds.
samples, unlike in iterative methods. We applied Eigen VI to many synthetic and real-world targets, and these experiments show that Eigen VI provides principled way of improving upon Gaussian variational families.
Many future directions remain. First, the approach described in this paper relies importance sampling, and thus it may benefit from more sophisticated methods for adaptive importance sam- pling. Second. it may be useful to construct variational families from different orthogonal function expansions, Our empirical study focused on the family built from normalized Hermite polynomials. But this family may require very high-order expansion to model highly non-Gaussian targets, and such an expansion will be very expensive in high dimensions. Though this family was sufficient for many of the targets we simulated, others will be crucial for modeling highly non- oon-Gassiian targets. Another direction is to develop variational families whose orthogonal function expansions scale more favorably with the dimension, perhaps by incorporating low rank structure the target's covariance. Finally, it would be interesting to explore iterative versions of EigenVI in which each iteration solves a minimum eigenvalue problem on some subsample of data points. With such an approach, Figen could potentially be applied to very large-scale problems in Bayesian inference.
10

---

Acknowledgments and Disclosure of Funding
We thank Bob Carpenter and Yuling Yao for helpful discussions and anonymous reviewers for their time and feedback on the paper. The Flatiron Institute is division of the Simons Foundation. This work was supported in part by NSF IIS-2127869, NSF DMS-231 1108, NSF/DoD PHY-2229929, ONR N00014-17-1-2131, ONR N00014-15-1-2200 the Simons Foundation, and Open Philanthropy.
References
[1] O. Abril-Pla, V. Andreani, C. Carroll, L. Dong, C. J. Fonnesbeck, M. Kochurov, R. Kumar, J. Lao, C. C. Luhmann, O. A. Martin, et al. PyMC: a modern, and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 9:e 2023.
[2] A. Agrawal, D. R. Sheldon, and J. Domke. Advances in black-box VI: Normalizing flows. importance weighting, and optimization. Advances in Neural Information Processing Systems, 33, 2020.
[3] R. v. d. Berg, L. Hasenclever, J. M. Tomczak, and M. Welling. Sylvester normalizing flows for variational inference. Unceriainty in Artificial Intelligence, 2018.
[4] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh,
P. Szerlip, P. Horsfall, and N. D. Goodman. Pyro: Deep universal probabilistic programming. The Journal of Machine Learning Research, 20(1):973 978, 2019.
[5] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 18:55177 2017.
[6] D. Cai, C. Modi, L. Pillaud- Vivien, and match: black-box variational inference with score-based divergence. In International Conference Machine Learning, 2024.
C. Margossian, R. Gower, D. Blei, and L. Saul. Batch
[7] 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.
[8] R. Courant and D. Hilbert. Methoden der Mathematischen Physik, volume 1. Julius Springer, Berlin, 1924.
[9] B. Dai, H. Dai, A. Gretton, Song, D. Schuurmans, and He. Kernel exponential family estimation via doubly dual embedding. In International Conference on Artificial Intelligence and Statistics. PMLR, 2019.
[10] A. K. Dhaka, A. Catalina, M. R. Andersen, M. Magnusson, J. Huggins, and A. Vehtari. Robust, accurate stochastic optimization for variational inference. Advances in Neural Information Processing Systems, 33, 2020.
[11] A. K. Dhaka, A. Catalina, M. Welandawe, M. R. Andersen, J. Huggins, and A. Vehtari. Challenges and opportunities in high dimensional variational inference. Advances in Neural Information Processing Systems, 34, 2021.
[12] L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real NVP. International Conference on Learning Representations, 2017.
[13] H. Ge, K. Xu, and Z. Ghahramani. Turing: a language for flexible probabilistic inference. In International Conference on Artificial Intelligence and Statistics. PMLR, 2018.
[14] S. Gershman, M. Hoffman, and D. Blei. Nonparametric variational inference. International Conference on Machine Learning. 2012.
[15] R. Giordano, M. Ingram, and T. Broderick. Black box variational inference with deterministic objective: Faster, more accurate, and even more black box. Journal of Machine Learning Research, 25(18):1-3,, 2024.

---

[16] D. Griffiths and D. F. Schroeter. Introduction to Quantum Mechanics. Cambridge University Press, 2018.
[17] Guo, X. Wang, K. Fan, T. Broderick, and D. B. Dunson. Boosting variational inference. arXiv preprint 1.1.::56/1000559 2016.
[18] A. Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
[19] C. Jones and A. Pewsey. Sinh-aresinh distributions. Biometrika, 96(4):661-78) 2009.
[20] C. Jones and A. Pewsey. The sinh-arcsinh normal distribution. Significance, 16(2):677, 2019,.
[21] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul. An introduction to variational methods for graphical models. Machine Learning, 37:183-233, 1999.
[22] Kim and Y. Bengio. Deep directed generative models with energy- based probability estima- tion. arXiv preprint arXiv: 03439, 2016.
[23] D. P. Kingma and M. Welling, Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
[24] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling. Improved vari- ational inference with inverse autoregressive flow. Advances in Neural Information Processing Systems, 29, 2016.
[25] I. Kobyzev, S. J. Prince, and M. A. Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11): 3964-3979, 2020.
[26] J. Köhler, A. Krämer, and F. Noé. Smooth normalizing flows. Advances in Neural Information Processing Systems, 34, 2021.
[27] A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei. Automatic differentiation variational inference. Journal of Machine Learning Research, 2017.
[28] J. Lawson, G. Tucker, B. Dai, and R. Ranganath. Energy-inspired models: Learning with sampler-induced distributions. Advances in Neural Information Processing Systems, 32, 2019.
[29] Y. LeCun, Chopra, R. Hadsell, M. Ranzato, and Huang. tutorial on energy-based learning. Predicting Structured Data, 1(0), 2006.
[30] R. B. Lehoueq, D. C. Sorensen, and C. Yang. ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM, 1998. Available at http:/ /www caam rice. edu/ http//ofwwrre///e boftware/ARPACK/ ACK
[31] Q. Liu and D. Wang. Stein variational gradient descent: a general purpose Bayesian inference algorithm. Advances in Neural Information Processing Systems, 29, 2016
[32] F. Locatello, G. Dresdner, R. Khanna, I. Valera, and G. Rätsch. Boosting black box variational inference. Advances in Neural Information Processing Systems, 31, 2018.
[33] L. Loconte, A. M. Sladek, S. Mengel, M. Trapp, A. Solin, N. Gillis, and A. Vergari. Subtractive mixture models via squaring: Representation and learning. In International Conference Learning Representations, 2024.
[34] C. Louizos and M. Welling. Multiplicative normalizing flows for variational Bayesian neural networks. In International Conference Machine Learning. PMLR, 2017.
[35] M. Magnusson, P. Bürkner, and A. Vehtari. posteriordb: set of posteriors for Bayesian inference and probabilistic programming. https://githut: htm.:/tth.oooooeeeeeooommooottte 2022.
[36] A. C. Miller, J. Foti, and R. P. Adams. Variational boosting: Iteratively refining posterior approximations. International Conference on Machine Learning, pages 2420-2429. PMLR, 2017.
12

---

[37] 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, 36, 2023.
[38] G. S. Novikov, M. E. Panov, and I. V. Oseledets. Tensor-train density estimation. In Uncertainty in Artificial Intelligence, pages 1321-1333 PMLR, 2021.
[39] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan. Nor- malizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 2021.
[40] R. Ranganath, S. Gerrish, and D. Blei. Black box variational inference. In Artificial Intelligence and Statistics, pages 814-822, PMLR, 2014.
[41] D. Rezende and S. Mohamed. Variational inference with normalizing flows. In International Conference on Machine Learning. PMLR, 2015.
[42] E. Roualdes, B. Ward, Axen, and B. Carpenter. BridgeStan: Efficient in-memory access to Stan programs through Python, Julia, and R. 2023.
[43] Salvatier, T. V. Wiecki, and C. Fonnesbeck. Probabilistic programming in Python using PyMC3. PeerJ Computer Science, 2:e55, 2016.
[44] M. Titsias and M. Lázaro-Gredilla. Doubly stochastic variational Bayes for non-conjugate inference. In International Conference on Machine Learning. PMLR, 2014.
[45] M. Wainwright, M. I. Jordan, et al. Graphical models, exponential families, and variational inference. Foundations and Trends@ in Machine Learning, 1(1-2):1-38) 2008.
[46] X. Wang, Geffner, and J. Domke. Dual control variate for faster black-box variational inference. In International Conference on Artificial Intelligence and Statisties, 2024.
[47] Y. Yang, R. Martin, and H. Bondell. Variational approximations using Fisher divergence. arXiv preprint
2019.
[48] L. Yu and C. Zhang. Semi-implicit variational inference via score matching, In International Conference on Learning Representations, 2023.
[49] Zeghal, Lanusse, A. Boucaud, B. Remy, and Aubourg. Neural Posterior Estimation with Differentiable Simulators. In International Conference on Machine Learning Conference, 2022.
[50] C. Zhang, B. Shahbaba, and H. Zhao. Variational Hamiltonian Monte Carlo via score matching. Bayesian Analysis, 13(2):485, 2018
[51] L. Zhang, B. Carpenter, A. Gelman, and A. Vehtari. Pathfinder: Parallel quasi- variational inference. Journal of Machine Learning Research, 23(306):1-49, 2022.
[52] S. C. Zhu, Y. Wu, and D. Mumford. Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling. International Journal of Computer Vision, 27: 107-126, 1998.
[53] D. Zoltowski, D. Cai, and Adams. Slice sampling reparameterization gradients. Advances in Neural Information Processing Systems, 34:23532-253 24332222344, 2021.
13

---

A Sampling from orthogonal function expansions
In this appendix we show how to sample from density on DD constructed from Cartesian product of orthogonal function expansions. Specifically, we assume that the density is of the form
K Σ
K
((21, ,
2D)
Ok.kg..kk Okg (2) )oi (22) UND (ED)
(A.1)
11=1
kp=1
where (+)(.)) define family of orthonormal functions on Il and where the density is normalized by requiring that
Σ kikg...k
(A.2)
hg...kk
To draw samples from this density, we describe sequential procedure based on inverse transform sampling. In particular, we obtain sample IR D by the sequence of draws
21 ((1),,
(A.3)
22 ~ g 21):
(A.4)
2D
(A.5)
This basic strategy can also be used to sample from distributions whose domains are Cartesian products of different one-dimensional spaces.
In what follows, we first introduce "core primitive" density, and we show how to sample efficiently from its distribution. We then show how the sampling procedure in Eqs. A.3-A.5 reduces to sampling from this core primitive; key component of this procedure the property orthogonality, which helps facilitate the efficient computation of marginal distributions.
Core primitive
First we describe the core primitive that we will use for each the draws in Eqs. A.3-A.5. To begin, we observe the following: if S is any positive semidefinite matrix with trace(S)= then
K
p(E)
Σ SeeD)(E)
(A.6)
R.E=1
defines a normalized density over R. In particular, since z O, follows that p(E) 20 for all E and since trace(5)==1) 1, =1, it follows that
[ [P(t)
Sue [ df
p())dE= de
P(L)
Σ
Σ k.t=1
Suedg trace(S) = 1.
(A.7)
k.C=1
The core primitive that we need is an efficient procedure to sample from normalized density of this form. We will see later that all of the densities in Eqs. A.3-A.5 can be expressed in this form.
Inverse transform sampling
Since the density in Eq. A.6 is one-dimensional, can obtain the draw need by inverse transform sampling. In particular, let C(E) denote the cumulative distribution function (CDF) associated with Eq. A.6, which is given by
C(E) = J p(z) dz,
(A.8)
and let CE) denote the inverse CDF. Then at least in principle, we can draw a sample from by the two-step procedure
u ~ Uniform (0, 1],
(A.9)
E = (u).
(A.10)
14

---

Next we consider how to implement this procedure efficiently in practice, and in particular, how to calculate the definite integral for the CDF in Eq. A.8. As shorthand, define the doubly- indbe-i set of real-valued functions
[ )z<)) dz.
(E)
(A.11)
follows from orthogonality that 2x((+00) = őạt and from the Cauchy-Schwarre inequality that (()e(E) for all E IR. Our interest these functions stems from the observation that
K Σ
C(E) =
Suc = trace (E)),
(A.12)
k,t=1
so that if we have already computed the functions 4z (5). then we can use Eq. A.12 to compute the CDF whose inverse need in Eq. A.10. In practice, (E) for many values along the real line and then solve Eq. A. quickly by interpolation; that is, given u, find satisfying trace[S#(()) density p(E) in Eq. A.6.
numerical quadrature pre-compute
The result is unbiased sample drawn from the
Sequential sampling
Finally we show that each draw in Eqs. A.3-A.5 reduces to the problem described above. As in Section 2, we work the steps specifically for example in D =3 where we must draw the samples q(z1). 22 ~ ((22(21) and 23 q(2a (2.1)) 22). This example illustrates all the ideas needed for the general case but with minimum of indices.
Consider the joint distribution given by
K2 K
((21, 22, 23)
ΣΣ>
Bijk di (2)))s))
where
Σ
B = 1. (A.13) ksk
From this joint distribution, we can compute marginal distributions by integrating out subsets of variables, and each integration over E gives rise contraction of indices, as in Eq. 7, due to the property of orthogonality. In particular, expanding the square in Eq. A.13, we can write this joint distribution
K;
K₁
K₂
q(21, 22, 23)
=
Σ
Σ
Σ
Bizk Buyw (1 (2))) (z
k,k''
j.
(A.14)
and we can then contract the index k when integrating over 23, since f (za) (2))) = 8k In this way we find that the marginal distributions are
22),
Begbeera ()(2))
(A.15)
Ri
K
Kas
((21)
Σ
Bish (i(z))o((a))
B
(A.16)
=,²=1
k=1
Now note from the brackets in Eq. A.16 that this marginal distribution is already in the quadratic form of Eq. A.6 with coefficients
K2
Ks
Bi
(A.17)
From this first quadratic form, we can therefore use inverse transform sampling to obtain draw 21 q(11).
15

---

Next we consider how to sample from the conditional ( = ((21.12)2///))) (2)) Again, from the brackets in Eq. A.15, we see that this conditional distribution also in the quadratic form Eq. A.6 with coefficients
f/ Eir=
Ks Σk 1=1 BennBeru
(A.18)
((1))
From this second quadratic form, we can therefore use inverse transform sampling to obtain draw
22 ((22|=1). Finally, we consider how sample from g(23 [1,]=1122) 22) q(21; 22, 23)/q(22 22). From Eq. A.14, we see that this conditional distribution is also in the quadratic form of Eq. A.6 with coefficients
S(3) Li Ens ye (o (z))
(A.19)
(
From this third quadratic form, we can therefore use inverse transform sampling to obtain draw 23 (223 21,21,20 22). Finally, from the sums in Eq. A.19, we see that the overall cost of this procedure OKF or quadratic in the total number of basis functions.
B Calculation of moments
In this appendix we show how calculate the low-order moments of density constructed from the Cartesian product of orthogonal function expansions. In particular, we assume that the density is over RD and of the form
Ki Σ
q(21, 22, ..,............................. ZD)
Σ kp=1
ka. Phi (21) (22) OKD (ZD)
(B.1)
where (o)())) are orthogonal functions on R and where the coefficients are properly normalized so that the density integrates to one. For such density, we show that the calculation of first and second-order moments boils down to evaluating one-dimensional integrals the form
Fig
(z))
(B.2)
M/ =
(z) dz.
(B.3)
We also show how to evaluate these integrals specifically for the orthogonal family of weighted Hermite polynomials.
First we consider how calculate moments such as [,, where p E {1,2), and without loss of generality we focus calculating [x]) We start from the joint distribution in Eq. B.I and proceed by marginalizing over variables (E2, 23...., zn). Exploiting orthogonality, we find that
B,(==1 = f ((21, .., ................................ ..). z dz dz2 dap,
(B.4)
Ki
Kp Σ
al kz kp Φk, (21) ... (20)) )
of dzi da2 ... .eddp (B.5)
=
kp=1
K₂
KD
ak akikg..kk ag kDD Onn (z1) (21))))((1))))) (21) dzy.
(B.6)
L
kp=1
We can rewrite this expression more compactly as a quadratic form over integrals of the form in Eqs. B.2-B.3. To this end, we define the coefficients
K2
KP
Aij
Oskz...kp Oijka...pp
(B.7)
16

---

which simply encapsulate the bracketed term in Eq. B.6. Note that there are K of these coefficients, each of which can be computed in O (K2K3 ... D). With this shorthand, we can write
B [x[21]
Aijuij
(B.8)
B [x1
Asjij
(B.9)
f,=1
where Pij and Pij are the integrals defined in Eqs. B.2-B.3. Thus the problem has been reduced to weighted sum of one-dimensional integrals.
similar calculation gives the result we need for correlations. Again, without loss of generality, we focus calculating [2] Analogous to Eq. B.7, we define the tensor of coefficients
Ka
Kp
Byjkt
Oskky Ap Ajcks...DD
(B.10)
which arises from marginalizing over the variables (23, 21:-.., 2D)- There K2K3 of these coefficients, each of which can be computed in O(KgK4 ... KD). With this shorthand, we can write
K1
K₂
[212212]
=
Σ
BiskeliMMee
(B.11)
1,1=1
K,E=1
where pay is again the integral defined in Eq. B.2). Thus the problem has been reduced weighted sum of (the product of) one-dimensional integrals.
Finally, we show how evaluate the integrals in Eqs. B.2-B.3 for the specific case of orthogonal function expansions with weighted Hermite polynomials; similar computations apply in the case of Legendre polynomials. Recall in this case that
(x+1 = (Virk)) (c- Hx(2),
(B.12)
where H₂ (z(z) are the probabilist's Hermite polynomials given by d
H (-1)*
(B.13)
dz
To evaluate the integrals for this particular family, we can exploit the following recursions that are satisfied by Hermite polynomials:
Hx+1(=) = ==) H₁((),, H²(=)= = k
(B.14)
(B.15)
Eliminating the derivatives Hz((2) in Eqs. B.14-B.15, we see that x = (=) + kfff-116 We then substitute Eq. B. obtain a recursion for the orthogonal basis functions themselves:
(2) Viob41
1k-1óx-112
(B.16)
With the above recursion, we can now read off these integrals from the property orthogonality. For example, starting from Eq. B.2, we find that
pij / (z))))())) ddz
(B.17)
[ di(2) [T0s+1(=) + VJ- 10)- dz,
(B.18)
- + 5 51.-1Vi,
(B.19)
17

---

where diy is the Kronecker delta function. Next we consider the integral in Eq. B.3, which involves power of 2² in the integrand. In this case we can make repeated use of the recursion:
Vij
d1(2)o((2) dz,
(B.20)
Vios+1(2) + Vi-I [Vo5+1((s) + Vī-10j-1(2) d2,
(B.21)
àb Vij V ((-11)(i- 1)6-11) + (6-1-11/00+1 + 1 i(++1).
(B.22)
Note that the matrices in Eq. B.19 and Eq. B.22 can be computed for whatever size is required by the orthogonal basis function expansion in Eq. B. Once these matrices are computed, is simple matter of substitution' to compute the moments [g[21] [z],, and Ey(212) 22] from Eqs. B.8-B.9 and Eq. B. Finally, we can compute other low-order moments (such By[25) or Bg (322)) by an appropriate permutation of indices.
C Eigenvalue problem
In this appendix show in detail how the optimization for Eigen reduces to minimum eigenvalue problem. In particular we prove the following.
Lemma C.1. Let {0x(2)}== be an orthogonal function expansion, and let E Qk be the variational approximation parameterized by
K Σ
((=) =
(C.1)
where the weights satisfy E o 1, thus ensuring that the distribution is normalized. Suppose furthermore that g chosen to minimize the empirical estimate of the Fisher divergence given, as in eq. (10), by
B
q(
(g.p))
Σ
V logg(z")) - V logp(=")11²
Then the optimal variational approximation in this family can be computed by solving the minimum cigenvalue problem
[~(())) = min a TM =: Amin(M) gE
min
(C.2)
=1
where M is given in Eq. 13 and (a1..... Øx] E R.K The optimal weights a are given (up an arbitrary sign) by the corresponding eigenvector this minimal eigenvalue.
Proof. The scores of q in this variational family are given by
2 2 a (24
V log a(2²)
E. (20)
With further bookkeeping, one can also exploit the sparsity of μıj and Miy derive more efficient calculations of these moments.
18

---

Substituting the above into the empirical divergence, we find that
(g.p)
V log ((=))
V [logp(=")]1
B
En (x)))
2 Ix
k ae (on(z")
v Vlogp(z")
a2²))
Σ (x0x(2²))
2
2 Σ
Σ Camdu(Eh) & Vlogp(z")
(2²b)
1
Σ OKk [2Von(=²) - (x((2)) VolVioog))))) 1ogn(="11
(==))
6=1
= a Ma.
where M is given in (13) and a [o1,-... 0k) R Thus the optimal weights a are found by minimizing the quadratic form a Ma subject to the constraint a a 1. Equivalently, solution can be found by minimizing the Rayleigh quotient
U Me
argmin
(C.3)
v T
and setting w/|lo]l. It then follows from the Rayleigh- Ritz theorem [8] for symmetric matrices that a is the cigenvector corresponding to the minimal eigenvalue of M, and this proves the lemma.
D Practical considerations of Eigen
D.1 EIgenVI vs gradient-based BBVI
Recall that Figen VT has two hyperparameters: the number of basis functions K and the number of importance samples B. We note there is an important difference hetween these two hyperparameters and the learning rate in ADVI and other gradient- based methods. Here, as we use more basis functions and more samples, the resulting fit is a better approximation. So, we can increase the number of basis functions and importance samples until a budget is reached, or until the resulting variational approximation is a sufficient fit. On the other hand, tuning the learning rate gradient- based optimization is much more sensitive because it cannot be too large or too small. If it is too large, ADVI may diverge. If the learning rate too small, it may take too long to converge in which case may exceed computational budgets.
Another fundamental difference in setting the number of basis functions compared to the learning rate or batch size of gradient based optimization is that once we have evaluated the score of the target distribution for the samples, these same samples can be reused for solving the eigenvalue problem with any choice of the number of basis functions, as these tasks are independent. By contrast, in iterative BB VI, the optimization problem needs be re-solved for every choice of hyperparameters, and the samples from different runs cannot be mixed together.
Furthermore, solving the eigenvalue problem is fast, and scores can be computed in parallel. In our implementation, we use off-the-shelf eigenvalue solvers, such as ARPACK [30] or Julia's eigenvalue decomposition function, eigen. In many problems with complicated targets, the main cost comes from gradient evaluation and not the eigenvalue solver.
D.2 Choosing the number of samples
Intuitively, if the target p is in the variational family Q (i.e., it can be represented using an order- K expansion), then we should choose the number of samples to roughly equal the number of basis function K. If p is very different from Q, we need more samples, and in our experiments, we use a multiple of the number of basis functions (say of order 10). As discussed before, once we have evaluated set of scores, these can be reused to fit larger number of basis functions.
19

---

3-component
101
Forward KL
Fomwar
10-2
2.00
Forward
10-1
10-4
ful-6
full-6
10-4
full10
full-6
ful-18
full-10
Gaussian
ful--10
Gaussian
Gaussian
1014
30
Number score samples
Number score samples
Number score samples
105
10-4
10
Forward KL.
Forward RL
10"2
10-2
20-1
10t
10
Figure E.1: We compare the number of score evaluations wallclock FKL divergence for the target distributions in Figure 3: the Gaussian mixture (column 1), the funnel (column 2), and the cross (column 3) distributions. The K used for EigenVI is reported in each figure legend, where K K2. The black star denotes the number of gradient evaluations for the Gaussian method.
E
Additional experiments and details
E.1
Computational resources
The experiments were processor and with 503 GB memory. Experiments were run on CPU. In the sinh-arcsinh and posteriordb experiments, computations construct the matrix M were parallelized over 28 threads.
Linux workstation with 32-core Intel(R) Xeon(R) w5-3435X
E.2 2D synthetic targets
We considered the following synthetic 2D targets:
3-component Gaussian mixture:
0.57) + 0.3 -101,,,,,,,5, p(=) 0.4 _ Σ) + 0.3N( [1.1,11 .0.57)
2 0.1
where we define Σ-
0.1
Funnel distribution with 0² = 1.2:
p(z) = N(= 0, o (x,ex(2))
Cross distribution:
p(z) = (2 [0, 2TT, + (2
4V [| -2)
+
12, E2) [0, Ei),
1..15.9
1 0.15
1
where 2
and 3g
0
These experiments were conducted without standardization with a Gaussian VI estimate. The Eigen proposal distribution π used was uniform([-9) 9)) distribution.
In Figure E.2, we run VI for increasing numbers of importance samples and report the resulting forward KL divergence. The blue curves denote variational families with different K₁
20

---

K2 K values used, i... 3, 6, and 10 (resulting in a total number of basis functions of 3 6², and 10%). the bottom row of the plot, we also show wall clock timings (computed without parallelization) to show how the cost grows with the increase in the number of basis functions and importance samples. The horizontal dotted line denotes the result from batch and match VI, which fits a Gaussian via score matching; here a batch size of 16 was used and a learning rate of X 2+I
BD
The black star denotes the number of score evaluations used by the Gaussian VI method.
E.3 Sinh-aresinh targets
The sinh-aresinh normal distribution [19, 20] has parameters RD, RP S++;
is
induced by transforming a Gaussian Zo ~ (O, Σ) to Z = Se,r (Zo). where
Sagsa = sinh ( sinh (aah (za) (z) = s (EI), (.p.....,,sp.r(3 (zn))
84
(E.1)
Here sa controls the amount of skew in the dth dimension, and Ta controls the tail weight in that dimension. When =0=0 and =x=1 1 in all dimensions d, the distribution is Gaussian.
The sinh-arcsinh normal distribution has the following density:
D
p(238,T, = [2x) |I
II d=1
(1 + = - Coare(za) exp (-255 (
(E.2)
where we define the functions
Csaraa (2a) = (1+ 52 (z))
(E.3)
and
Sania( == sinh(Ta sinh -1 (za) - sa), (x) [S (z.)..... (ED))
(E.4)
We constructed 3 targets in 2 dimensions and 3 targets in dimensions, each with varying amounts of non-Gaussianity. The details each target are below. all experiments, EigenVI was applied with standardization, where a Gaussian was fit using batch and match VI with batch size of 16 and learning rate A ++1'
For all experiments, we used a proposal distribution π that was uniform on [-5,5]3.
2D sinh-aresinh normal experiment For D 2 (Figure 4b), we consider the slight skew and tails targer with parameters
0.2), [1.1, 1.1), the more skew and tails target with = [0.2, 0.5), = [1.1, and the slight skew and heavier tails with = [0.2, 0.2), = [1.4. 1.1]. Note that [0, 0), [1, 1] recovers the multivariate Gaussian. These three target are visualized in Figure 4a.
5D sinh-archsinh normal experiment We constructed three targets P1 (slight skew and tails), Pa (more skew and tails), and P3 (slight skew and heavier tails) each with
2.2
0.3
0
0
0.37
0.3
2.2
0
Σ
0
0
2.2
0.3
0
(E.S)
0
0.3
2.2
0.3
0
0
0
2.2
The skew and tail weight parameters used were:
[0., 0., 0.2.
1. 1.1]
82 = 00.0, 0.0, 0.6, 0.4, 0.5] To 11, 1., 1.1], and 83 = 0022,00.,,02,,,22,,23.23
[1.1, 1.4, 1.6). See Figure E.2 for a visualization of the marginals of each target distribu- tion. the second row, we show examples of resulting Figen VI fit (visualized using samples from ②)
from B 20,000 and K 10.
21

---

Slight skew and talls
More skew and talls
Slight skew and heavier tails
You 3.
8009 3882 2292 201
JAN 201 1898338300000
03-5.88-03
KL-4.8E.00
KL-3.11.02
3,299
Figure E.2: Targets (top) for the 5D sinh-av mach-arcsinh normal distribution example and Eigen VI fits (bottom) with the KL divergence in the figure title.
Table E.I: Summary of posteriordb models
Name
Dimension
Model description
kidscore
linear model with a Cauchy noise prior
seaame
3
linear model with uniform prior
gP_regr
3
Gaussian process regression with squared exponential kernel generalized autoregressive conditional heteroscedastic model
garch11
4
logearn
log-log linear model with multiple predictors autoregressive model for time series multiple predictors log-log model non-centered hicrarchical model for 8-schools
arK-arK
7
logmesquite
7
8-schools
10
E.4 Posteriordb experiments
We consider 8 real data targets from posteriordb, a suite of benchmark Bayesian models for real data problems. In Table 1, we summarize the models considered in the study. These target distributions are non Gaussian, typically with some skew or different tails. To access the log target probability and their gradients, we used the BridgeStan library [42], which by default transforms the target to be supported on RDD
For all experiments, we fixed the number of importance samples to be = 40,000; to construct the Eigen VI matrix M, the computations were parallelized over the samples. These experiments were repeated over 5 random seeds, and we report the mean and standard errors in Figure 5; for lower dimensions, there was little variation between runs.
The target distributions were standardized using a Gaussian fit from score matching before applying Eigen VI. In most cases, the proposal distribution æ chosen to be uniform over [-6, 6 For the models 8-schools, which has a longer tail, we used multivariate Gaussian proposal with zero mean and scaled diagonal covariance oT, with = 32.
For the Gaussian score matching (GSM) method [37], we chose a batch size 16 for experiments. We generally found the results were not too sensitive in comparison to other batch sizes of 4,8, and
22

---

EigenVi
BganWi
ADVI with normalizing izing flow
Batch and match VI
HMC
Figure E.3: Comparison of EigenVI, normalizing flow, and Gaussian score-based BBVI methods on 8schools.
23

---

EigenVi
INC BigeM
Normalizing flow
HMC
C 2
50 si
Gaussian score-based BBVI
25
45
359
Figure E.4: Comparison of Eigen normalizing flow, and Gaussian score-based BBVI methods on garchi1. Note that the Gaussian approximation over/underestinatts the tails, while the more expressive families fit the tails better.
24

---

32. For the batch and match (BaM) method [6], we chose batch size of 16. The learning rate was fixed at X BD which was a recommended schedule for Gaussann targets.
r+I
For all ELBO optimization methods (full covariance Gaussian family and normalizing flow fam- ily), we used Adam to optimize the ELBO, We performed grid search over the learning rate 0.01, 0.02, 0.05, 0.1 and batch size B 4, 8, 16, 32. For the normalizing flow model, used real NVP [12] with 8 layers and 32 neurons. We found empirically that the computational of the scores unreliable [26, 49]; hence we do not show their Fisher divergence in Figure 5.
In Figure E.3 and Figure E.4, we show the corner plots that compare an Eigen VI fit, a normalizing flow fit, and Gaussian fit (BaM). In each plot, we plot the samples from the variational distribution against samples from Hamiltonian Monte Carlo. We observe that the two more expressive families EigenVI VI and the normalizing flow able model the tails of the distribution better than the Gaussian fit.
Broader impacts
EigenVI adds the literature BBVI, which has been important line of work for developing automated, approximate Bayesian inference methods. In terms of positive societal impacts, Bayesian models are used throughout the sciences and engineering, and advances in fast and automated inference will aid in advances in these fields. In terms of negative societal impacts, advances in BBVI could be used train generative models with malicious unintended uses.
25

---

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: We describe propose Eigen VI, new eigenvalue-based approach for BBVI based on score matching.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
Justification: Limitations are discussed throughout the paper, and we summarize the main limitations again Section
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-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on 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?
26

---

Answer: [Yes]
We provide further derivations for moments, sampling, and the objective the appendix.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proof' 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 encouraged 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
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We have aimed to provide these details in the additional experiments section 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 novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either 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 reprodueibility may be tricky in some cases, in which authors are welcome to 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 be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
27

---

The authors should provide scripts to reproduce all experimental results for the new reproducible, they Question: 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?
Answer: [No]
Justification: We provide a Julia implementation of Eigen VI at https://githut. com/ di dicai/eigen// and a demonstration on several examples.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.c cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
/ httppsscccppbccccceeeeceeecc guwi/uidue des/ CoodeSubmis ssi aesionPolicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
proposed method and baselines. only subset of experiments should state which ones are omitted from the script and why.
At submission time, to 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 is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide additional experimental details 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 a level of detail that is necessary 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:
We report the standard errors as shaded regions around the mean over 5 random
seeds.
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, 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 with given experimental conditions).
28

---

The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 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 figures symmetric error bars that would yield results that are out of range (c.g. negative emror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables 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 experimental compute resources 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, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: : //t.r//eerrpp.eeeee htcpppbllcc/thhcccuteeeeeeeeeeee
Answer: [Yes]
Justification:
Yes we conform to the code of ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in 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:
We have a broader impacts section in the appendix.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or 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.
29

---

The conference expects that many papers will be foundational research and not tied 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 intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We do not believe the paper poses such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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. Lieenses for existing assets
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?
Question:
Answer: [Yes]
Justification: We cite the benchmark models used.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include 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 (c.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 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.
30

---

If this information is not available online, the authors are encouraged to reach out to the asset's '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:
We do not release assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
We do not perform research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing 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 in the main paper.
According the NeurIPS Code of Ethies, 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: There were no study participants used.
Guidelines:
The answer NA means that the paper does involve crowdsourcing research with human subjects.
Depending on 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 we 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 as the institution conducting the review.
31


Reviews:
[{"review_id": "dh9RII9ZUh", "reviewer": "Reviewer_XdZT", "summary": "The authors introduce a novel family of variational distributions based on orthogonal function expansions which is optimized by minimizing a Fisher divergence. They show that an unbiased importance sampling estimate of the divergence results in a quadratic form, which can be minimized by finding its smallest eigenvector.", "questions": "- The authors describe a procedure to sample from the variational approximation, so I am not sure why importance sampling needed to estimate the Fisher divergence? Is it purely for computational reasons or am I missing something?\n- How is the proposal distribution chosen that is used for importance sampling?", "limitations": "The authors have adequately addressed the limitations of their work.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The work is original and of high overall quality. The proposed variational family and training methodology is novel and of interest to the variational inference community. The empirical evaluation is fairly extensive and shows that the method works well for low dimensional problems. At the same time the authors are very forward with discussing the current limitations of their method which opens up interesting questions for future research.", "weaknesses": "As discussed by the authors, the work seems to be mostly limited by the fact that higher-order function expansions, which are needed to model non-Gaussianity, become exceedingly expensive in high dimensions as a result of the blow up of the size of the eigenvalue problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your interest in the paper. We will revise the final version of the paper to make the following points more clear.\n\n> 1. The authors describe a procedure to sample from the variational approximation, so I am not sure why importance sampling is needed to estimate the Fisher divergence? Is it purely for computational reasons or am I missing something?\n\nFirst, because we are optimizing the empirical Fisher divergence with respect to $q$, we need the samples to be independent of $q$ (we cannot simultaneously sample from and optimize over $q$). In addition, as we show in Appendix D, the importance weights are important and lead to the form of the eigenvalue problem. \n\nThe procedure to sample from $q$ is primarily used after the variational distribution has been fit to, for instance, approximate functionals of $p$.\n\nWe will revise the paragraph starting from Line 133 to explain these points more thoroughly. \n\n> 2. How is the proposal distribution chosen that is used for importance sampling?\n\nPlease see our general response above for intuition behind choosing the proposal distribution. See Appendix E for the exact proposal distribution used for each experiment. \n\nWe will add a paragraph in the Section 2.3 that details this intuition behind choosing the proposal distribution.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qYBPIC9YWw", "reviewer": "Reviewer_jsA6", "summary": "The paper proposes Eigen-VI, a black-box variational inference method that uses orthogonal function expansions to parameterize the variational distribution and uses Fisher divergence as an objective. The method does not require gradient-based optimization method and behaves well in a set of synthetic and real application experiments.", "questions": "* In the experiment part of this paper, why does the author only use the Hermite polynomials as the variational family? Is this because other variational family performs badly or have intractable issues?", "limitations": "See Weaknesses above.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "Strengths\n* The paper is well-written and easy to follow.\n* The proposed Eigen-VI is simple and easy to implement since the optimization problem can be solved by computing the minimum eigenvalue of a given matrix, and it also achieves comparable or better performances on hierarchical modeling benchmarks compared to other black-box variational inference methods.", "weaknesses": "Weaknesses\n* Compared to other BBVI methods, Eigen-WI restricts the variational families to be $K^{\\text {th}}$-order variational family(defined in eq. 2). When generalize to other variational families, gradient-free optimization in this paper may fail. Similarly, if variational distribution is set to be a linear combination of function basis, the optimization problem in other BBVI methods would also be much simpler. Therefore, the contribution of this paper may not be very significant.\n\n* Different proposal distributions in importance sampling may affect the property of the matrix $M$, and further affect the numerical efficiency or stability of computing $\\lambda_{\\text{min}}(M)$. Adding some analysis regarding this matter would enhance the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your feedback on the paper. We discuss several of the comments and questions below.\n\n> 1. Compared to other BBVI methods, Eigen-VI restricts the variational families to be 𝐾th-order variational family(defined in eq. 2). When generalize to other variational families, gradient-free optimization in this paper may fail. \n\nAs the reviewer points out, the EigenVI approach works for any orthogonal basis set. Importantly, this approach does not just work for a single variational family but a very large class of variational families. In particular, there are many orthogonal basis sets beyond the ones listed in the paper in Table 1. For instance, one may apply a procedure such as Gram-Schmidt to obtain such a set. Finally, new variational families can be formed by even mixing different basis sets. Thus, this approach opens the door to applying VI to many new types of families.  \n\n> 2. Similarly, if variational distribution is set to be a linear combination of function basis, the optimization problem in other BBVI methods would also be much simpler. \n\nClassical BBVI methods also have assumptions that restrict the variational families. For instance, it is common that the family supports reparameterization (which is not true here). \n\n> 3. Different proposal distributions in importance sampling may affect the property of the matrix M, and further affect the numerical efficiency or stability of computing λmin(M). Adding some analysis regarding this matter would enhance the paper.\n\nIn principle we agree that a poor proposal distribution might cause such numerical issues; however, in the current experiments, we do not observe this. Standardizing the distribution as we currently do it also helps here, as it provides a good proposal distribution which can easily be adjusted to have heavier tails if necessary. Furthermore, with enough samples, and since solving for minimum eigenvalues is typically so fast even for large matrices, we expect these issues to be less important. Experimenting with different proposal distributions is beyond the scope of this rebuttal, but we acknowledge the good point raised and will add a remark to this end in the revised manuscript on how, e.g., this could affect the spectral gap.\n\n> 4. In the experiment part of this paper, why does the author only use the Hermite polynomials as the variational family? Is this because other variational family performs badly or have intractable issues?\n\nThe parameters for most experiments considered lie on an unconstrained scale where Hermite polynomials form a natural basis; the Hermite family can also be seen as a natural extension of non-Gaussianity.\n\nThe EigenVI method itself applies to _any orthogonal basis set_; in Figure 1, we additionally show 1D target distributions and fitted distributions arising from Legendre polynomial and Fourier series expansions.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bfol71rQeP", "reviewer": "Reviewer_MKWn", "summary": "This paper proposes EigenVI, which uses orthonormal distribution functions as the basis of the variational distribution family $q$. Then, minimizing the difference between $q$ and the target distribution $p$ via minimizing the Fisher divergence (2-norm score distance) is turned into an eigen-decomposition problem. Authors show that EigenVI is very competitive in distribution approximation compared with alternatives.", "questions": "* Eq. 9, should be $\\mathrm{d} z$? Line 130, the score of the target is not equal to the gradient of the log joint, but up to a scaling?\n* Could the authors explain a bit more why centering the distribution is very important for better approximation?\n* Is there any guidance about how to choose a good proposal distribution? How do authors choose their proposal distributions used in the experiments? Is the choice of the proposal distribution very significant or quite robust?\n* Eq. 16 seems a bit weird, especially the term $L_z(x_{1:N})$. Besides, it is better to clearly state whether $x_{1:N}$ are independent samples or a sequence, and how it relates to a single $z$.\n* Why is forward KL used in the whole paper rather than the commonly used backward KL as in ELBO?\n* How about the actual running time comparison on synthetic and real-world experiments? I understand that theoretically comparison is hard because eigenVI's complexity is w.r.t. number of basis (order) but the stochastic algorithm's complexity is mainly w.r.t. # epochs. However, solving large eigen-decomposition problems in practice might be very time-consuming. So it is very helpful to have an empirical understanding of the running time of the algorithm.", "limitations": "/", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The presentation is clear. The math is solid.\n* The proposed EigenVI is very intuitive and reasonable.\n* The low-order basis and higher-order bases have their own interpretation of the distribution's feature, which is desirable.\n* Experiments are extensive, including synthetic and real-world data.", "weaknesses": "See Questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your feedback on the paper. In the revised manuscript, we will work on clarifying the following points.\n\n> 1. Should Equation 9 be dz?\n\nWe could alternatively write $q(z) dz$ here (with the appropriate adjustments of $q(z)$ and $p(z)$).\n\n> 2. Line 130: the score of the target is not equal to the gradient of the log joint but up to a scaling?\n\nThe log target is equal to the log joint + constant, and so the constant disappears when we take the gradient. Thus, $\\nabla \\log p(z) = \\nabla \\log \\rho(z)$, where $\\rho$ represents the unnormalized model (i.e., the joint in the context of Bayesian inference).\n\n> 3. On Equation 16 and the likelihood:\n\nThis equation describes a general set of Bayesian problems from a set of benchmark models, and so we don’t a priori make assumptions on the conditional independence structure of the likelihood. We are happy to write $p(x_{1:N} | z)$ instead.\n\n> 4. Could the authors explain a bit more why centering the distribution is very important for better approximation?\n\nThe standardization is used with the Hermite polynomial basis, which has a centered Gaussian as the base distribution. As we show in Figure 2, higher order expansions are needed to model distributions whose largest mode is shifted away from the origin. While this is fine in low dimensions, in order to handle higher dimensional targets, we apply the standardization technique as a way to reduce the order needed to model the target. \n\n> 5. Why is forward KL used rather than the commonly used backward KL as in ELBO?\n\nThe ELBO is often used when comparing VI algorithms that use the ELBO / reverse KL as the objective for all VI algorithms. Here we are comparing VI algorithms with _different_ objectives, and so we chose an evaluation metric that was agnostic to the algorithms’ objectives.\n\n> 6. How to choose a good proposal distribution? How were they chosen in experiments? \n\nPlease see our general response above for intuition behind choosing the proposal distribution. See Appendix E for the exact proposal distribution used for each experiment. \n\n> 7. How about the actual running time comparison on synthetic and real-world experiments? [...] Solving large eigen-decomposition problems in practice might be very time-consuming. \n\nWe do not need to solve the full eigendecomposition problem, but only need the smallest eigenvalue pair for a KxK symmetric matrix, which can be a much easier problem to solve. As we mention in the general rebuttal comment, for 3000 basis functions, which is around the largest number of basis functions considered in experiments, this takes under a second. In a qualitative context, in the rebuttal pdf, we report the time it takes to run iterative BBVI methods, and we report the time needed for the eigenvalue solve. For instance, our implementation of ADVI takes over 12 seconds to converge in this example.\n\nThus overall, we expect the overhead in terms of computational time due to solving the eigenvalue problem to be minimal.\n\nThere are two other steps in (Hermite) EigenVI - standardization and evaluating scores of importance samples. While the actual computational time for these depends on the efficiency of implementation, expect this to be competitive to iterative BBVI approaches for two reasons -- 1) Standardization is done using BBVI approaches themselves (like ADVI, BaM or GSM), and hence the cost of standardization is the same as the cost of doing VI inference using these. 2) The scores for the importance samples can be evaluated in a fully parallelized manner and hence, depending on the implementation, adds minimal overhead in terms of computational time.\n\nFinally, to present a quantitative comparison, we previously included a running time comparison for the synthetic 2D experiments, see Appendix E (Figure 6 shows both # of gradient evaluations and wallclock time without parallelization). Here we did not need to use standardization, and so we indicate when EigenVI and the iterative method for the Gaussian fit use the same number of gradient evaluations.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1Rc6kJuTbq", "reviewer": "Reviewer_35nb", "summary": "The authors propose EigenVI, a new method for black-box variational inference (BBVI).  The method uses the Fisher divergence (rather than the more typical reverse KL divergence) and a variational family based on orthogonal function expansions (which I've never seen used before).  Advantages of the approach include the ability to model complex distributions (multimodal, skewed, heavy-tailed, etc.) while also giving closed-form calculations for lower-order moments as well as easy posterior sampling.", "questions": "Clarifications:\n\n(1) pp.6: What is the proposed advantage of EigenVI over Zhang et al.'s energy-based approach with a closed-form solution?\n\n(2) pp.7, Figure 4:  If the normalized Hermite polynomials contain Gaussian distributions as a \"base\" (i.e. when the number of basis functions is K=1), then why does EigenVI still beat ADVI when K=1?  I expected them to perform identically in this setting. \n\n(3) pp.8: The authors state that they \"standardize the target using either GSM or BaM before applying EigenVI\".   Any reason that they pick those methods relative to others? Would they expect any arbitrary BBVI method (e.g. ADVI of Kucukelbir et al. 2017) to work just as well for this purpose?\n\n(4) pp.9:  Should \"this family is limited to target functions that are close to Gaussian\" read \"... that are close to mixtures of Gaussians\"?  For example, the cross-distribution (Figure 3, row 3) does not look close to Gaussian to me, although it does look close to a mixture of Gaussians. \n\nTypos:\n\n(1) Abstract: \"novel class variational approximations\" should read \"novel class of variational approximations\".", "limitations": "Yes, authors describe limitations well.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper is exceptionally strong in originality, clarity, and significance.\n\nFor originality,  I have not previously seen a paper in machine learning using orthogonal function expansions to define a variational family.  The authors state that the approach comes from  quantum mechanics, an atypical source domain for a machine learning paper.  The construction is interesting, and normalizes elegantly.  Given the importance of orthogonal bases throughout mathematics, it is likely that this approach to variational inference will connect with many fields of academic inquiry.  Thus, this approach has the potential to greatly enrich research in variational inference.\n\nFor clarity, the exposition is exceptionally clear nearly everywhere (although Sec 2.3 could use a couple iterations of rewriting to streamline the narrative).   In particular, orthogonal function expansions and score-based divergences were introduced in an exceedingly clear (and interesting) manner.  The high-level advantages of their method vs. others in the literature were also exceedingly clearly presented. \n\nFor significance, the approach is broadly useful (at least for low-dimensional posteriors); this is inherited from the fact that BBVI is constructed so as to be applicable to a wide set of applications without requiring any pen-and-paper calculations.   Moreover, by avoiding gradient-based inference in favor of solving an eigenvalue problem, their method would additionally remove the common headache of having to specify or learn a learning rate.", "weaknesses": "Some weaknesses:\n\n(1) As mentioned by the authors, the number of basis functions scales exponentially in the dimensionality of the posterior, and therefore so does inference time.   Given the restriction to low-dimensional posteriors, why use VI at all, rather than MCMC sampling?  It would have been nice to see this addressed in the exposition.  It would have also been nice to see an experimental comparison showing some advantage over using MCMC samplers.   For instance, do the authors expect a computational advantage of their approach over MCMC methods?  If so, it's hard to trust this since computation times are not given.\n\n(2) Following up on the preceding point, it is completely unclear how long it takes to run this method, since no computation times are given relative to competitors.  On pp.9, the authors say, \"In these examples, we standardize the target using either GSM or BaM before applying EigenVI; for this reason, we do not compare the costs of these methods to EigenVI.\"  I found this sentence quite puzzling.  Presumably the authors intend to argue that since EigenVI must take strictly longer than GSM or BaM,  there's not point in reporting the computation times.  If so, I don't see how the conclusion follows. A practitioner would want to know how MUCH additional compute is required to obtain any advantages.  Moreover, how does the overall (combined) computation time compare to other approaches, such as VI with a normalizing flow family or standard ADVI (Kucukelbir et al. 2017)?\n\n(3) The method seems fairly complex to run, since a competitor BBVI method must be run simply as a preprocessing step to standardize the distribution.  (They use either Gaussian Score Matching or Batch-and-Match VI, per pp. 8.)\n\n(4) A major selling point emphasized by the authors is that unlike gradient-based approaches, EigenVI does not require tunable hyper-parameters, such as learning rates or termination criteria.  However, EigenVI is clearly sensitive to (a) the number of basis functions and (b) the number of importance samples. Indeed, Figure 4 (c) demonstrates how inference can go catastrophically bad when the number of importance samples is not sufficiently large relative to the number of basis functions. This problem is reminiscent of the problem of a poorly chosen learning rate in standard ADVI (Kucukelbir et al., 2017).  Can the authors provide any guidance on how to choose a good number of basis functions for a given problem, and then, given that, how to choose a good number of importance samples?  A theoretical result would be nice, although admittedly it's unreasonable to expect that during the rebuttal period.  Do they at least have any practical insight to share?\n\n(5) I feel that the authors undersell their own method by neglecting to reinforce the primary selling points in the experiments section.  For example, the authors mention that their approach allows sampling, closed-form moment calculations, and modeling of distributions with mixed support. However, none of these things  appear explicitly in the experimental section.   Similarly, the authors emphasize that the normalizing flow competitor is \"often sensitive to the hyperparameters of the flow architecture, optimization algorithm, and parameters of the base distribution\" (pp.6).  Can the authors demonstrate that the posterior approximation deteriorates along these dimensions more easily than their own approach deteriorates as a function of number of importance samples and basis functions?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your feedback; in the revised manuscript, we will add or expand our discussion to address several of the points you bring up, as detailed below and in the main rebuttal comment.\n\n> 1. A practitioner would want to know how MUCH additional compute is required to obtain any advantages. \n\nWhile the actual overhead in terms of computational time depends on the efficiency of implementation, we provide qualitative reasons to justify why this is competitive to other approaches:\n\nStandardization can use BBVI approaches (like ADVI, BaM or GSM), and hence the cost of standardization is the same as the cost of doing VI inference using these. Following this, we need to evaluate the scores for importance samples which can be fully parallelized and hence, depending on the implementation, adds minimal overhead in terms of computational time. The final step is to solve the eigenvalue problem which is independent of the dimensions and is fast in the regimes we consider (see general comment and attached pdf for details).\n\n> 2. EigenVI is clearly sensitive to (a) the number of basis functions and (b) the number of importance samples. [...] This problem is reminiscent of the problem of a poorly chosen learning rate in standard ADVI \n\nWe will modify the text to clarify that we do have two hyper-parameters: the number of basis functions and the number of importance samples. But we do find there is an important difference between our two hyperparameters and the learning rate in ADVI and other gradient-based methods. As we use more basis functions and more samples, the resulting fitted $q$ is a better approximation. So, we can increase the number of basis functions and importance samples until a budget is reached, or until the resulting q is a good enough fit. Tuning the learning rate is much more sensitive because it cannot be too large or too small. If it is too large, ADVI may diverge. If it is too small, it may take too long to converge. \n\nAnother fundamental difference in setting the # of bases as compared to the learning rate or batch size of gradient based optimization is that once we have evaluated the score of the target distribution for the samples, these same samples can be reused for solving the eigenvalue problem with any choice of the number of basis functions, as these tasks are independent. By contrast, in iterative BBVI, the optimization problem needs to be solved from scratch for every choice of hyperparameters, and the samples from different runs cannot be mixed together. Furthermore, solving the eigenvalue problem is fast, and scores can be computed in parallel.\n\nFinally, for choosing the proposal distribution, in the case of the Hermite family, we have reasonable defaults because we are standardizing and hence can use e.g. an isotropic Gaussian with a larger variance to have stable importance weights. \n\n> 3. Can the authors provide any guidance on how to choose a good number of basis functions for a given problem, and then, given that, how to choose a good number of importance samples? [...] Do they at least have any practical insight to share?\n\nIn terms of practical insights, if the target p is in the variational family Q, then we need the number of samples = # of basis functions. If it is very different from Q, we need more samples, and we use a multiple of the number of basis functions (say of order 10). As discussed before, once we have evaluated a set of scores, these can be reused to fit a larger number of basis functions.\n\n> 4. Similarly, the authors emphasize that the normalizing flow competitor is \"often sensitive to the hyperparameters of the flow architecture, optimization algorithm, and parameters of the base distribution\"\n\nIn the rebuttal pdf Fig. 1, we focus on a particular model considered in the paper, and we show the sensitivity of posterior estimates to the learning rate and the batch size for the normalizing flow and standard ADVI. \n\n> 5. What is the proposed advantage of EigenVI over Zhang et al.'s energy-based approach?\n\nThe energy-based model is not normalized, and thus it cannot be evaluated or sampled from (without using MCMC). Indeed, their paper details an approach that uses HMC with the energy-based VI approach. \n\n> 6. If the normalized Hermite polynomials contain Gaussian distributions as a \"base\" (i.e. when the number of basis functions is K=1), then why does EigenVI still beat ADVI when K=1?\n\nIn this case, the Gaussian fit by BaM is used as the standardizing distribution. Thus, K=1 can beat ADVI if the standardizing distribution converges faster / to a better distribution.\n\n> 7. The authors state that they \"standardize the target using either GSM or BaM before applying EigenVI\". Any reason that they pick those methods relative to others? Would they expect any arbitrary BBVI method (e.g. ADVI of Kucukelbir et al. 2017) to work just as well for this purpose?\n\nWe discuss the motivation behind the standardization in Section 2.3; here you can pick your favorite Gaussian approximation to do this, e.g., a Laplace approximation or ADVI or directly estimating the target mean and covariance (we chose the score-based VI methods in the experiments because they require less tuning and often converge faster). In general, we expect a better estimate of the target mean and covariance to help due to the base distribution of the Hermite family being a standard Gaussian (and thus, lower-order expansions will be more accurate). \n\n> 8. Should \"this family is limited to target functions that are close to Gaussian\" read \"... that are close to mixtures of Gaussians\"? \n\nWe agree this sentence could be more clear, and we’ll modify it. What we meant was that to get to higher-dimensional spaces, we are limited to target functions that are close to Gaussian. \nFor lower dimensional models (e.g., the 2D targets), we can afford to take higher-order expansions, and are able to model highly non-Gaussian distributions; we generally found this to be true for up to 3 or 4 dimensions.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "11uOQw3rEg", "reviewer": "Reviewer_JZVK", "summary": "This work proposes a new variational family based on orthogonal function expansions from the quantum mechanics literature. Furthermore, the paper proposes to find the basis coefficients by solving a score matching problem, where the problem cleverly reduces to an eigenvalue problem. The overall approach, called EigenVI, is compared to recent score-matching-based Gaussian VI algorithms and BBVI.", "questions": "* In Line 103 to 105, it is stated that one can circumpass distributions with constrained supports. Given bijectors (ADVI), constrained supports are rarely a problem for VI algorithms these days. But dealing with *discrete support* is still a major challenge. Would there be a way to use the proposed method for problems with discrete support?\n* Eq (10): Why is importance sampling used here when one can sample directly from $q$?", "limitations": "The paper would benefit from a dedicated limitation section.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "* The idea of using orthogonal function expansions is new, although it reminds me of some other non-parametric/almost parametric variational families.\n* The fact that the score-matching problems reduce to an eigenvalue problem is beautiful in its own right.\n* The writing of the article is the best in my batch.", "weaknesses": "* While the authors are upright about the major limitation of the work, scaling with respect to dimensionality, this, unfortunately, strongly limits the practicality of the work. To me, for a VI algorithm to be useful, it should at least be scalable with respect to $d$, the dimensionality of the posterior, or $n$, the number of datapoints. Unfortunately, the proposed approach seems to achieve neither. In particular, the eigenvalue problem formulation is very cute but it isn't clear if subsampling results in an unbiased objective to be scalable with respect to $n$. (Please correct me if this is not the case!)\n* The evaluation is not adequate for a work proposing a new variational family. The two baseline variational families considered here are the full-rank Gaussian variational family and normalizing flows. There is a rather long history of developing less-parametric variational families in VI. For instance, boosting VI [1], copulas [2,3], splines [4], semi-implicit families [5], mixtures [6], and many more. These works should be the main focus of the evaluation, not Gaussian VI.\n* Unless I am not mistaken, the method suffers in the case the posterior is not *a priori* standardized. It doesn't seem like such standardization would be possible unless some preliminary inference has already been performed. This would defeat the point of using a VI algorithm. If one has already obtained preconditioners, why not just provide them to an MCMC algorithm instead?\n* The paper doesn't discuss how to solve the eigenvalue problem. Given that we only need to solve a specific eigenvalue problem, there might be clever ways to solve it in a scalable way, maybe by using randomized linear algebra methods, which could technically be interesting.\n\nOverall, I hope to see at least some evidence that the proposed methodology will eventually lead to a practical algorithm in the future. For instance, whether subsampling results in an unbiased/consistent algorithm. If this is the case, I would be happy to re-assess the paper. \n\n## Additional Comments\n* The Fisher divergence and the score-matching problem were introduced to the machine-learning community through the work of Hyvarinen [7], which isn't cited.\n* Related works: A discussion on the literature of non-parametric/\"less-parametric\" VI algorithms seems necessary. \n* Line 30 \"One thread of BBVI research focuses on Gaussian variational approximations. (ADVI, BBVI)\": I wouldn't quite agree with this statement since ADVI is applicable to any family where the reparameterization gradient is applicable. In fact, when normalizing flows are used with ELBO maximization, this is pretty much the same as ADVI. Furthermore, the original BBVI paper [8] requires conjugacy, meaning that one cannot simply use Gaussian variational families at all in a lot of cases.\n* A lot of the references are citing the arXiv version. Berg et al. 2018 was published at UAI; Dinh et al. 2016 was published at ICLR 2017; Giordano et al. 2023 was published in JMLR; and etc. Please thoroughly review the reference section.\n* Line 39, Section 2.1: Citing something when mentioning orthogonal functions would be better so that readers interested in the topic could have a look into more in-depth material.\n\n## References\n1. Campbell, T., & Li, X. (2019). Universal boosting variational inference. NeurIPS\n2. Smith, M. S., & Loaiza-Maya, R. (2023). Implicit copula variational inference. Journal of Computational and Graphical Statistics, 32(3), 769-781.\n3. Han, S., Liao, X., Dunson, D., & Carin, L. (2016, May). Variational Gaussian copula inference. In AISTATS.\n4. Shao, Y., Shan, N. Y., & Feng, T. (2024, April). Nonparametric Automatic Differentiation Variational Inference with Spline Approximation. In AISTATS.\n5. Yu, L., & Zhang, C. (2023). Semi-Implicit Variational Inference via Score Matching. arXiv preprint arXiv:2308.10014.\n6. Lambert, M., Chewi, S., Bach, F., Bonnabel, S., & Rigollet, P. (2022). Variational inference via Wasserstein gradient flows. NeurIPS.\n7. Hyvärinen, A., & Dayan, P. (2005). Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4).\n8. Ranganath, R., Gerrish, S., & Blei, D. (2014, April). Black box variational inference. In AISTATS.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your feedback on the paper. We discuss several of the points below.\n\n> 1. It isn't clear if subsampling results in an unbiased objective \n\nAny unbiased estimate of the gradient to our objective results in an unbiased estimate of the objective. And so EigenVI can scale to large n with subsampling.\n\nLet $p(z,x)$ be the target joint distribution between data $x$ and latent variables $z$.\n\nLet $g(z)$ be an unbiased estimate of the score $\\nabla \\log p(z,x)$, i.e., $\\mathbb{E}[g(z) \\mid z] = \\nabla \\log p(z,x)$.\n\nExpanding the Fisher divergence,\n\n$\\mathbb{E}[\\|\\nabla \\log q(z_b) - g(z_b)\\|^2 \\mid z_b] = \\|\\nabla \\log q(z_b)\\|^2 - 2 \\langle \\nabla \\log q(z), \\mathbb{E}[g(z_b) \\mid z_b] \\rangle + \\mathbb{E}[\\|g(z_b)\\|^2 \\mid z_b]$.\n\nUsing the unbiased property, this becomes\n\n$\\|\\nabla \\log q(z)\\|^2 - 2 \\langle \\nabla \\log q(z), \\nabla \\log p(z,x) \\rangle + \\mathbb{E}[\\|g(z_b)\\|^2 \\mid z_b]$.\n\nSo we get \n \n$\\|\\nabla \\log q(z) - \\nabla \\log p(z,x)\\|^2 + \\mathbb{E}[\\|g(z_b)\\|^2 \\mid z_b] - \\|\\nabla \\log p(z,x)\\|^2$.\n\nThus, $\\|\\nabla \\log q(z_b) - g(z_b)\\|^2$ is an unbiased estimate of $\\|\\nabla \\log q(z) - \\nabla \\log p(z,x)\\|^2$ up to constants that do not depend on $q$. Because of this, we can use $g(z_b)$ as a drop-in replacement within the Fisher divergence within our method since \n    \n$\\min_{q \\in \\mathcal{Q}} \\mathbb{E}\\_{q(z)} [|\\nabla \\log q(z) - \\nabla \\log p(z,x)|^2] = \\mathbb{E}_{q(z)}[|\\nabla \\log q(z) - g(z)|^2] + C.$\n\nTo build this unbiased estimator of the score, we can use data subsampling as follows: suppose data is sampled i.i.d such that \n\n$p(z,x) = \\prod_{i=1}^n p(z,x_i)$.\n\nThus, \n\n$\\nabla \\log p(z,x) = \\sum_{i=1}^n \\nabla \\log p(z,x_i)$.\n\nWe can use this to get unbiased estimates of the full score by sampling the data $x_i$. For instance, by sampling $x_i \\sim \\frac{1}{n}$ we have that \n\n$\\mathbb{E}[n \\nabla \\log p(z_b, x_i) \\mid z_b] = \\frac{1}{n} \\sum_{j=1}^n n \\nabla \\log p(z_b, x_i) = \\nabla \\log p(z, x)$.\n\nWe could also use any form of mini-batching over data to also build an unbiased estimate of the gradient.\n\n> 2. Scaling with dimension\n\nThere are many directions for future work to help scale to higher dimensions. For instance, in many situations, one might expect non-Gaussianity on a low-dimensional subspace, and that the remaining dimensions can be suitably modeled with a Gaussian. This family is a special case of the orthonormal Hermite family considered in the paper and so is still an eigenvalue problem. Another direction is a low-rank approximation of the tensor – this family could be optimized with gradient-based methods. Finally, one may consider refining an existing basis to develop better basis functions. \n\n> 3. Why not use VI fit as a preconditioner for MCMC?\n\nAs noted in above, MCMC struggles with targets with varying curvature (using a preconditioner does not help with multiscale distributions where curvature varies over different regions of the distribution, e.g., Neal's funnel). In the absence of pathological geometries, HMC can work very well in high-dimensions and so dimension itself doesn’t strike us as the relevant criterion to choose between VI or MCMC. Using VI with MCMC is an active research area, and this motivates developing new VI methods to benefit from their synergies.\n\n> 4. Comparing to non-parametric VI families\n\nWe thank the referee for pointing to additional literature on non-parametric BBVI families. We will refer to these in the related works in the revised manuscript. However given the breadth of this literature, individual comparison with each of these is beyond the scope of this work.\n\nThere are a few common elements to all these approaches, e.g., they all use gradient based optimization and hence are similarly sensitive as the baselines that we have considered. In addition, there are some individual challenges, e.g., mixtures are prone to mode collapse, and spline VI assumes a factorized posterior distribution.\n\nGiven these considerations, we consider our baselines of full rank Gaussian VI and normalizing flows (NFs) to be adequate. Together, these are the most widely used variational families for BBVI, they encompass the two extremes on the spectrum of parametric vs non-parametric families, and especially NFs being universal approximators should typically encompass other variational families suggested by the reviewer.\n\nFinally, note that we do not suggest that EigenVI results in better final approximation than NFs, but instead our focus is introducing a new variational family and approach to inference which does not rely on typical iterative gradient based optimization and associated hyperparameter tuning (which is common to all aforementioned methods).\n\n> 5. The paper doesn't discuss how to solve the eigenvalue problem. \n\nWithout further information on the scores of p, the resulting matrix in our eigenvalue problem is dense and unstructured, thus there is little we can leverage to design a specialized algorithm. In our experiments, we used off-the-shelf eigenvalue solvers. For typical problems, this is not a bottleneck (see the general response on the discussion for computational cost). We will add a discussion of this.\n\n> 6. References and related work \n\nWe will expand our discussion of related nonparametric families, and we will add references to additional literature on orthogonal families and score matching. \n\n> 7. Given bijectors, constrained supports are rarely a problem. Would there be a way to use the proposed method for problems with discrete support?\n\nTransforming supports to be real-valued can sometimes make the problem more challenging (and adds an additional level of tuning to the problem). If we know that the variables are bounded, it may be more natural to model them in the original space. Regarding discrete support, this is an interesting direction for future work. With an appropriate basis set and replacing the scores with log ratios, one would still get an eigenvalue problem.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for their time and engagement – we are pleased that they found the work to be highly original, well-written, and of interest to the variational inference community. 

In this work, we propose a novel variational family and an efficient algorithm to fit it via score matching. Here variational distributions are obtained by _squaring_ a linear combination of orthogonal basis functions – we demonstrate the flexibility of this family using several basis sets (e.g., Hermite, Fourier series, and Legendre polynomials) in Figure 1, but note that our approach applies to _any orthogonal basis set_, paving the way for many other possible families. We show that minimizing an estimate of the Fisher divergence (which matches the scores of q and p) is equivalent to solving for the minimum eigenvalue of some matrix. EigenVI avoids gradient-based optimization, and the scores and the matrix of interest can be evaluated and formed in parallel.

_**After thorough consideration of reviewers’ comments, we will expand our discussion of several points in the manuscript. Our response includes additional plots of hyperparameters and times for baseline methods.**_ 

## Score-based VI vs sampling & MCMC [Reviewers 35nb, JZVK]

First, a score-based approach has potential advantages over a purely sampling-based approach that does not exploit the availability of scores. These advantages are perhaps most easily understood by considering the simple problem of estimating the mean and variance of a one-dimensional Gaussian distribution, $p(z)$. By drawing $n$ samples from $p$, one can estimate the mean and variance to accuracy $O(1/\sqrt{n})$. However, one can determine the mean and variance _exactly_ from the scores $d\log(p)/dz$ at just two samples $z_1$ and $z_2$, provided that $z_1\neq z_2$. Most settings are not this contrived, but the broader point remains: the scores of a target distribution $p$ provide a great deal of information to constrain the search for the best variational approximation $q$ in some parameterized family of tractable distributions. Methods that do not exploit this information will generally be at a considerable disadvantage to those that do.

In addition, reviewers have asked specifically about MCMC vs EigenVI. We agree MCMC may be applied to several of the examples in the paper and a comprehensive comparison with MCMC strikes us as an important direction for future work. However, most off-the-shelf algorithms, such as Stan’s adaptive Hamiltonian Monte Carlo, do not work well on targets with a varying curvature (e.g. funnel) or well separated modes (e.g. cross), even with preconditioners, in low or high dimensions. In all cases, tuning MCMC—and in particular (i) determining the length of the warmup and sampling phase of the Markov chain, and (ii) deciding on a number of chains, are non-trivial tasks which severely impact the reported performance. Therefore, benchmarking against MCMC requires careful study. 

## Cost of EigenVI & solving the eigenvalue problem [Reviewers 35nb, MKWn, JZVK]

EigenVI needs to compute only the smallest eigenpair of a K x K symmetric matrix, where K is the # of basis functions, and not the full eigendecomposition. We use up to several thousand basis functions. In this regime, standard solvers for the minimum eigenvalue take under a second; e.g., ARPACK, a popular backend wrapped in Python, Julia, and Matlab that implements the restarted Lanczos method, takes 510 ms to find the smallest eigenvalue of a 3000 x 3000 randomly generated positive definite matrix. When the dimension K is very large, too large to store the matrix M, we can use iterative solvers for computing the smallest eigenvalue, e.g., gradient descent applied to the Rayleigh coefficient $\alpha^\top M \alpha/\alpha^\top\alpha $. Each iteration of gradient descent would require a matrix-vector product, and costs O(K^2). Furthermore, we do not need to explicitly form the matrix M, but need only implement the operation $\alpha \mapsto M \alpha,$ and thus never need to store the K x K matrix. 

Finally, we note in many problems with complicated targets, the main cost comes from gradient evaluation and not the eigenvalue solver. 

## Why use importance sampling if we can sample directly from q? [Reviewers XdZT, JZVK]

Because we are optimizing the empirical Fisher divergence with respect to q, we need the samples to be independent of q (we cannot simultaneously sample from and optimize over q, and it’s unclear how to apply reparameterization). In addition, the importance-sampled estimate leads to the quadratic form and thus the eigenvalue problem.

## Proposal dist. for EigenVI [Reviewers XdZT, jsA6, 35nb]

In the particular case of the Hermite variational family, we can first standardize the target distribution using the mean and covariance of a Gaussian fit. Thus, intuitively, we want a proposal that has heavier tails than a standard Gaussian. We found reasonable defaults to be centered distributions such as a uniform (if most of the target’s mass is between certain bounds) or multivariate Gaussian with long enough tails. 

## EigenVI does not require standardization, but it helps [Reviewers 35nb, JZVK]

If the variables are bounded, a uniform distribution is a reasonable base distribution, and no standardization is needed. 

For the Hermite family, we do not need standardization in the lower dimensions (no standardization is used for the 2D targets in the experiments). Standardization is used specifically with the Hermite family to help scale the method to larger dimensions, as the base distribution of the Hermite family is a standard Gaussian. (As we show in Fig. 2 for a 1D example, uncentered distributions may require more basis functions, so standardization helps to reduce the # of basis functions needed.) Given standardization with a Gaussian distribution, one can also view EigenVI as a post-processing step of a Gaussian BBVI (with any BBVI technique of choice--ADVI, BaM or GSM).

Meta-review: This paper had 5 reviewers instead of the normal 4 (or minimum of 3), and yet all reviewers were unanimous: this is an original paper of high quality. The idea is novel, the evaluation is solid, and the paper is honest about shortcomings. The topic, Variational Inference, is of wide-ranging interest in ML.

Therefore I'm pleased to recommend acceptance.

For the camera-ready version, reviewer JZVK strongly wants the authors to do a thorough review of non-parametric variational families. Please make any other changes requested by the reviewers that have not already been addressed in the revisions.
Decision: Accept (spotlight)