EigenVI: score-based variational inference with orthogonal function expansions

Authors: diana  cai
chirag  modi
charles c margossian
robert m gower
david m blei
lawrence k saul
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 03:37:40

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.

Main paper: 
 introduction: 
 Probabilistic modeling is a cornerstone of modern data analysis, uncertainty quantification, and decision making. A key challenge of probabilistic inference is computing a target distribution of interest; for instance, in Bayesian modeling, the goal is to compute a posterior distribution, which is often intractable. Variational inference (VI) [5, 21, 45] is a popular method for scalable probabilistic inference that has worked across a 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 (NeurIPS 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, 37] . Instead of minimizing the KL 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 as 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 a 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; at 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 use orthogonal function expansions to construct distributions supported on R D ; 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 a variational family from this class, we 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 a minimum eigenvalue problem, thus avoiding the need for gradient-based methods. For this reason, we call our approach EigenVI.We study EigenVI with a variational family constructed from weighted Hermite polynomials. We first demonstrate the expressiveness of this family on a variety of multimodal, asymmetric, and heavytailed distributions. We then use EigenVI to approximate a diverse collection of non-Gaussian target distributions from posteriordb [35] , a benchmark suite of Bayesian hierarchical models. On these problems, EigenVI provides more accurate posterior approximations than leading implementations of Gaussian BBVI based on KL minimization and score-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 in these families reduces to an eigenvalue problem. In Section 3 we review the literature related to EigenVI. In Section 4, we evaluate EigenVI on a variety of synthetic and real-data targets. Finally, in Section 5, we discuss limitations and future work. 
score-based variational inference with orthogonal function expansions: 
 In this section we use orthogonal function expansions to develop new variational families for approximate probabilistic inference. In Section 2.1, we review the basic properties of these expansions. In Section 2.2, we introduce a score-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.Let Z ⊆ R D denote the support of the target distribution p. Suppose that there exists a complete set of orthonormal basis functions {ϕ k (z)} ∞ k=1 on this set. By complete, we mean that any sufficiently wellbehaved function f : Z → R can be approximated, to arbitrary accuracy, by a particular weighted (•) z ∈ [-1, 1] Legendre polynomials {1, z, 3z 2 -1, 5z 3 -3z, . . .} z = e iθ ∈ S 1 Fourier basis {1, cos θ, sin θ, cos 2θ, sin 2θ, . . .} z ∈ [0, ∞) weighted Laguerre polynomials e -z 2 {1, 1-z, z 2 -4z +2, . . .} z ∈ Rweighted Hermite polynomials e -z 2 4 {1, z, (z 2 -1), (z 3 -3z), . . .} sum of these basis functions, and by orthonormal, we mean that the basis functions satisfyEQUATIONwhere the integral is over Z. Define the K th -order variational family Q K to be the set containing all distributions of the formEQUATIONand where α k ∈ R for k = 1, . . . , K are the parameters of the family Q K . In words, Q K contains all distributions that can be obtained by taking weighted sums of the first K basis functions and then squaring the result.Eq. 2 involves a squaring operation, a sum-of-squares constraint, and a weighted sum. The squaring operation ensures that the density functions in Q K are nonnegative (i.e., with q(z) ≥ 0 for all z ∈ Z), while the sum-of-squares constraint ensures that they are normalized:EQUATIONThe weighted sum in Eq. 2 bears a superficial similarity to a mixture model, but note that neither the basis functions ϕ k (z) nor the weights α k in Eq. 2 are constrained to be nonnegative. Distributions of this form arise naturally in physics from the quantum-mechanical wave functions that satisfy Schrödinger's equation [16] . In that setting, though, it is typical to consider complex-valued weights and basis functions, whereas here we only consider real-valued ones.The simplest examples of orthogonal function expansions arise in one dimension. For example, functions on the interval [-1, 1] can be represented as weighted sums of Legendre polynomials, while functions on 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, we may consider approximations of the form in Eq. 2 whereEQUATIONand H k (z) are the probabilist's Hermite polynomials given byEQUATION)Note how the lowest-order basis function ϕ 1 (z) in this family gives rise (upon squaring) to a Gaussian distribution with zero mean and unit variance.Figure 1 shows how various multimodal distributions with one-dimensional support can be approximated 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 EQUATIONwhere β ijk ∈ R now parametrize the family. Note that there are a total K 1 K 2 K 3 parameters in the above expansion, so that this method of Cartesian products does not scale well to high dimensions if multiple basis functions are used per dimension. Note that the same strategy can also be used for random variables of mixed type: for example, from Table 1 , we can create a variational family of distributions over R×[-1, 1]×[0, ∞) from the Cartesian product of orthogonal function expansions involving Hermite, Legendre, and Laguerre polynomials.As shown in Figure 1 , the approximating distributions from K th -order expansions can model the presence of multiple modes as well as many types of 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 A and B.For concreteness, consider the distribution over R 3 in Eq. 6. The marginal distribution q(z 1 ) isEQUATIONand from this expression, moments such as E[z 1 ] and Var[z 1 ] can be calculated by evaluating integrals involving the elementary functions in To generate samples {z (t) }, each dimension is sampled as follows: we draw z (t)1 ∼ q(z 1 ) by computing the cumulative distribution function (CDF) of this marginal distribution and then numerically inverting this CDF. Finally, extending these ideas, we can calculate higher-order moments and obtain joint samples via the nested drawsEQUATIONThe 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.In variational inference, we posit a 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 Q K from the orthogonal functions {ϕ k (z)} K k=1 whose variational parameters are the weights {α k } K k=1 . We now derive EigenVI, a method to find q ∈ Q K that is close to the target distribution p(z).We first define the measure of closeness that we will minimize. EigenVI measures the quality of an approximate density by the Fisher divergence [18] ,EQUATIONwhere ∇ log q(z) and ∇ log p(z) are the score functions of the variational approximation and target, respectively. Suppose that q and p 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 ∇ log p at any point z ∈ Z. For example, in Bayesian models the score of the target 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 q ∈ Q K that minimizes D(q, p). But now a challenge arises: it is generally difficult to evaluate the integral for D(q, p) in Eq. 9, let alone to minimize it as a 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 {α k } K k=1 in terms of which it is defined. Instead, we construct an unbiased estimator of D(q, p) by importance sampling, which also decouples the sampling distribution from the optimization. Let {z 1 , z 2 , . . . z B } denote a batch of B samples drawn from some proposal distribution π on Z. From these samples we can form the unbiased estimatorEQUATION)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 we show that the minimization of Eq. 10 over q ∈ Q K simplifies to a minimum eigenvalue problem for the weights {α k } K k=1 . To obtain the eigenvalue problem, we substitute the orthogonal function expansion in Eq. 2 into Eq. 10 for the unbiased estimator of D(q, p). As an intermediate step, we differentiate Eq. 2 to obtain the scoresEQUATION)Further substitution of the scores provides the key result behind our approach: the unbiased estimator in Eq. 10 is a simple quadratic form in the weights α := [α 1 , . . . , α K ] ⊤ of the orthogonal function expansion,EQUATIONwhere the coefficients of the quadratic form are given byEQUATION)Note that the elements of the K ×K symmetric matrix M capture all of the dependence on the batch of samples {z b } B b=1 , the scores of p and q at these samples, and the choice of the family of orthogonal functions. Next we minimize the quadratic form in Eq. 12 subject to the sum-of-squares constraintk α 2 k = 1 in Eq. 2.In this way we obtain the eigenvalue problem [8] minEQUATIONwhere λ min (M ) is the minimal eigenvalue of M , and the optimal weights are given (up to an arbitrary sign) by its corresponding eigenvector; see Appendix C for a proof. EigenVI solves Eq. 14.We note that the eigenvalue problem in EigenVI arises from the curious alignment of three particular choices-namely, (i) the choice of variational family (based on orthogonal function expansions), (ii) the choice of divergence (based on score-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 heuristics of gradient-based optimizations-involving learning rates, terminating criteria, and perhaps other algorithmic hyperparameters-that are typically required for ELBO-based BBVI [10, 11] . But EigenVI is also not entirely free of heuristics; to compute the estimator in Eq. 10 we must also specify the proposal distribution π and the number of samples B; see Appendix D for a discussion.The size of the eigenvalue problem in EigenVI 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 K 1 K 2 K 3 for the approximation in Eq. 6, as can be seen by "flattening" the tensor of weights β in Eq. 6 into the vector of weights α = vec(β) in Eq. 2. Finally, we note that EigenVI only needs to 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.We now discuss the specific case of EigenVI for Z = R D with the Hermite-based variational family in Eq. 4. For this case, we propose a transformation of the domain that serves to precondition or standardize the target distribution before applying EigenVI. While this standardization is not required to use EigenVI, it helps to reduce the number of basis functions needed to approximate the target, leading to a 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 EigenVI, our goal is therefore to transform the domain in a way that reduces the number of basis functions needed for a good approximation. To meet this goal for distributions over R D , we observe that the lowest-order basis function of the Hermite family in Eq. 4 yields (upon squaring) a standard multivariate Gaussian, with zero mean and unit covariance. Intuitively, we might expect the approximation of EigenVI to require fewer basis functions if 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, EigenVI in R D can then be viewed as a systematic framework to model non-Gaussian effects via a small number of higher-order terms in its orthogonal function expansion.Concretely, we consider a linear transformation of the domain:EQUATIONwhere µ 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-specific knowledge).We then apply the EigenVI to fit a Kth-order variational approximation q(z) 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 p(z), i.e.,EQUATION)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 1 (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 EigenVI; when the target has nearly zero mean and unit covariance, it becomes simpler to identify natural choices for the proposal distribution π. Intuitively, in this case, we want a proposal distribution that has the same mean but heavier tails than a standard Gaussian. In our experiments, we use two types of centered proposal distributions-uniform and isotropic Gaussian-whose variances are greater than one. 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 1 after translating its largest modes away from the origin.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 in Section 5.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
related work: 
 Several recent works have considered BBVI methods based on score-matching. These methods take a particularly simple form for Gaussian variational families [6, 37] . The Fisher divergence [18] has been previously studied as a divergence for variational inference [47] . Yu and Zhang [48] propose minimizing a Fisher divergence for semi-implicit (non-Gaussian) variational families; the divergence is minimized using gradient-based optimization. In another line of work, Zhang et al. [50] consider variational families of energy-based models and derive a 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 easily evaluate or draw samples from the densities of energy-based models.The variational approximation in EigenVI 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 in 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 al. [33] obtain distributions by squaring a mixture model with negative weights, and they study this model in conjunction with probabilistic circuits. By contrast in this work, we consider this idea in the context of variational inference, and we focus specifically on the use of orthogonal function expansions, which have many simplifying properties; additionally, the specific objective we optimize leads to a minimum eigenvalue problem.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We evaluate EigenVI on 9 synthetic targets and 8 real data targets. In these experiments we use the orthogonal family induced by normalized Hermite polynomials (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 of these higher-order expansions. Next, we experiment with target distributions where we systematically vary the tail heaviness and amount of skew. Finally, we apply EigenVI to a set of hierarchical Bayesian models from real-world applications and benchmark its performance against other Gaussian BBVI algorithms.We first demonstrate how higher-order expansions of the variational family yield more accurate approximations on a range of 2D non-Gaussian target distributions (Figure 3 ); see Appendix E. Gaussian: KL=2.723 0 3 K = 6 2 , KL=0.12 3 0 3 K = 8 2 , KL=4.7E-02 3 0 3 K = 14 2 , KL=2.3E-02Figure 3 : 2D target functions (column 1): a 3-component Gaussian mixture distribution (row 1), a 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 3-5), whereK = K 1 K 2 .for details. We report an 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 not standardized before fitting EigenVI (we compare the costs of the methods in Figure E.1), and the total number of basis functions isK = K 1 K 2 .We now consider the sinh-arcsinh 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 non-Gaussianity in the skew or the tails of the distribution, and we refer to these targets as slight skew and tails, more 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 EigenVI fits along with their forward KLs. Before applying EigenVI, we standardize the target using a mean and covariance estimated from batch and match VI [6] . In Figure 4b , we measure the EigenVI forward KL under varying numbers of samples B and across increasing numbers of basis functions, given by K = D d=1 K d . We also present the forward KL resulting from batch and match VI (BaM) and automatic differentiation VI (ADVI), which both use Gaussian variational families 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 EigenVI variational approximations. In Figure 4c , we observe greater differences in the number of importance samples needed to lead to good approximations, especially as the number of basis functions increase.We now evaluate EigenVI on a set of hierarchical Bayesian models [7, 35, 42] , which are summarized in Table E .1. The goal is posterior inference: given data observations x 1:N , the posterior of z isEQUATIONwhere p(z) is the prior and p(x 1:N | z) denotes the likelihood.We compare EigenVI to 1) automatic differentation VI (ADVI) [27] , which maximizes the ELBO over a full-covariance Gaussian family (ADVI), 2) Gaussian score matching (GSM) [37] , a score-based BBVI approach with a full-covariance Gaussian family, and 3) batch and match VI (BaM) [6] , which In these models, we do not have access to the target distribution, p(z | x 1:N ), only the unnormalized target ρ. Thus, we cannot evaluate an estimate of the forward KL. Instead, to 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):EQUATIONNote that this measure is not the objective that EigenVI minimizes; it is analogous to the forward KL divergence, as the expectation is taken with respect to p. We report the results in Figure 5 , computing the Fisher divergence for EigenVI 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 a qualitative comparison with real-NVP normalizing flows (NFs) [12] , a flexible variational family that is 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 a subset of dimensions from 8schools in the top three rows, comparing EigenVI, the NF, and BaM. Here, we observe that the Gaussian struggles to fit the tails of this target distribution. On the other hand, EigenVI provides a competitive fit to the normalizing flow. In Appendix E. samples, unlike in iterative methods. We applied EigenVI to many synthetic and real-world targets, and these experiments show that EigenVI provides a principled way of improving upon Gaussian variational families.Many future directions remain. First, the approach described in this paper relies on importance sampling, and thus it may benefit from more sophisticated methods for adaptive importance sampling. 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 a 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-Gaussian targets.Another direction is to develop variational families whose orthogonal function expansions scale more favorably with the dimension, perhaps by incorporating low rank structure in 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, EigenVI could potentially be applied to very large-scale problems in Bayesian inference.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, we define the doubly-indexed set of real-valued functionsΦ kℓ (ξ) = ξ -∞ ϕ k (z)ϕ ℓ (z) dz. (A.11)It follows from orthogonality that Φ kl (+∞) = δ kl and from the Cauchy-Schwartz inequality that |Φ kℓ (ξ)| ≤ 1 for all ξ ∈ R. Our interest in these functions stems from the observation thatC(ξ) = K k,ℓ=1S kℓ Φ kl (ξ) = trace[SΦ(ξ)], (A.12) so that if we have already computed the functions Φ kℓ (ξ), then we can use Eq. A.12 to compute the CDF whose inverse we need in Eq. A. 10 . In practice, we can use numerical quadrature to pre-compute Φ kℓ (ξ) for many values along the real line and then solve Eq. A.10 quickly by interpolation; that is, given u, we find ξ satisfying trace[SΦ(ξ)] = u. The result is an unbiased sample drawn from the density ρ(ξ) in Eq. A.6.Finally we show that each draw in Eqs. A.3-A.5 reduces to the problem described above. As in Section 2.1, we work out the steps specifically for an example in D = 3, where we must draw the samples z 1 ∼ q(z 1 ), z 2 ∼ q(z 2 |z 1 ) and z 3 ∼ q(z 3 |z 1 , z 2 ). This example illustrates all the ideas needed for the general case but with a minimum of indices.Consider the joint distribution given byq(z 1 , z 2 , z 3 ) =   K1 i=1 K2 j=1 K3 k=1 β ijk ϕ i (z 1 )ϕ j (z 2 )ϕ k (z 3 )   2 where ijk β 2 ijk = 1. (A.13)From this joint distribution, we can compute marginal distributions by integrating out subsets of variables, and each integration over R gives rise to a 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 asq(z 1 , z 2 , z 3 ) = K3 k,k ′ =1   K1 i,i ′ =1 K2 j,j ′ =1 β ijk β i ′ j ′ k ′ ϕ i (z 1 )ϕ i ′ (z 1 )ϕ j (z 2 )ϕ j ′ (z 2 )   ϕ k (z 3 )ϕ k ′ (z 3 ),(A.14) and we can then contract the index k ′ when integrating over z 3 , since ϕ k (z 3 )ϕ k ′ (z 3 )dz 3 = δ kk ′ .In this way we find that the marginal distributions areq(z 1 , z 2 ) = K2 j,j ′ =1   K1 i,i ′ =1 K3 k=1 β ijk β i ′ j ′ k ϕ i (z 1 )ϕ i ′ (z 1 )   ϕ j (z 2 )ϕ j ′ (z 2 ), (A.15) q(z 1 ) = K1 i,i ′ =1   K2 j=1 K3 k=1 β ijk β i ′ jk   ϕ i (z 1 )ϕ i ′ (z 1 ). (A.16)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(1)ii ′ = K2 j=1 K3 k=1 β ijk β i ′ jk . (A.17)From this first quadratic form, we can therefore use inverse transform sampling to obtain a draw z 1 ∼ q(z 1 ).Next we consider how to sample from the conditional q(z 2 |z 1 ) = q(z 1 , z 2 )/q(z 1 ). Again, from the brackets in Eq. A.15, we see that this conditional distribution is also in the quadratic form of Eq. A.6 with coefficients(2)jj ′ = K1 i,i ′ =1 K3 k=1 β ijk β i ′ j ′ k ϕ i (z 1 )ϕ i ′ (z 1 ) q(z 1 ) . (A.18)From this second quadratic form, we can therefore use inverse transform sampling to obtain a draw z 2 ∼ q(z 2 |z 1 ). Finally, we consider how to sample from q(z 3 |z 1 , z 2 ) = q(z 1 , z 2 , z 3 )/q(z 1 , z 2 ). From Eq. A.14, we see that this conditional distribution is also in the quadratic form of Eq. A.6 with coefficients(3)kk ′ = K1 i,i ′ =1 K2 j,j ′ =1 β ijk β i ′ j ′ k ′ ϕ i (z 1 )ϕ i ′ (z 1 )ϕ j (z 2 )ϕ j ′ (z 2 ) q(z 1 , z 2 ) (A.19)From this third quadratic form, we can therefore use inverse transform sampling to obtain a draw z 3 ∼ q(z 3 |z 1 , z 2 ). Finally, from the sums in Eq. A. 19 , we see that the overall cost of this procedure isO(K 2 1 K 2 2 K 2 3), or quadratic in the total number of basis functions.In this appendix we show how to calculate the low-order moments of a density constructed from the Cartesian product of orthogonal function expansions. In particular, we assume that the density is over R D and of the formq(z 1 , z 2 , . . . , z D ) = K1 k1=1 • • • K D k D =1 α k1k2...k D ϕ k1 (z 1 )ϕ k2 (z 2 ) • • • ϕ k D (z D ) 2 , (B.1)where {ϕ k (•)} ∞ k=1 are orthogonal functions on R and where the coefficients are properly normalized so that the density integrates to one. For such a density, we show that the calculation of first and second-order moments boils down to evaluating one-dimensional integrals of the formµ ij = ∞ -∞ ϕ i (z)ϕ j (z) z dz, (B.2) ν ij = ∞ -∞ ϕ i (z)ϕ j (z) z 2 dz. (B.3)We also show how to evaluate these integrals specifically for the orthogonal family of weighted Hermite polynomials.First we consider how to calculate moments such as E q [z p d ], where p ∈ {1, 2}, and without loss of generality we focus on calculating E q [z p 1 ]. We start from the joint distribution in Eq. B.1 and proceed by marginalizing over the variables (z 2 , z 3 , . . . , z D ). Exploiting orthogonality, we find thatE q [z p 1 ] = q(z 1 , z 2 , . . . , z D ) z p 1 dz 1 dz 2 . . . dz D , (B.4) = K1 k1=1 • • • K D k D =1 α k1k2...k D ϕ k1 (z 1 )ϕ k2 (z 2 ) • • • ϕ k D (z D ) 2 z p 1 dz 1 dz 2 . . . dz D , (B.5) = K1 k1,k ′ 1 =1 K2 k2=1 • • • K D k D =1 α k1k2...k D α k ′ 1 k2...k D ϕ k1 (z 1 )ϕ k ′ 1 (z 1 ) z p 1 dz 1 . (B.6)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 coefficientsA ij = K2 k2=1 • • • K D k D =1 α ik2...k D α jk2...k D , (B.7)which simply encapsulate the bracketed term in Eq. B.6. Note that there are K 2 1 of these coefficients, each of which can be computed in O(K 2 K 3 . . . K D ). With this shorthand, we can writeE q [z 1 ] = K1 i,j=1 A ij µ ij , (B.8) E q [z 2 1 ] = K1 i,j=1 A ij ν ij , (B.9)where µ ij and ν ij are the integrals defined in Eqs. B.2-B.3. Thus the problem has been reduced to a weighted sum of one-dimensional integrals.A similar calculation gives the result we need for correlations. Again, without loss of generality, we focus on calculating E q [z 1 z 2 ]. Analogous to Eq. B.7, we define the tensor of coefficientsB ijkℓ = K3 k3=1 • • • K D k D =1 α ikk3...k D α jℓk3...k D , (B.10)which arises from marginalizing over the variables (z 3 , z 4 , . . . , z D ). There are K 2 1 K 2 2 of these coefficients, each of which can be computed in O(K 3 K 4 . . . K D ). With this shorthand, we can writeE q [z 1 z 2 ] = K1 i,j=1 K2 k,ℓ=1 B ijkℓ µ ij µ kℓ . (B.11)where µ ij is again the integral defined in Eq. B.2). Thus the problem has been reduced to a weighted sum of (the product of) one-dimensional integrals.Finally, we show how to 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ϕ k+1 (z) = √ 2πk! -1 2 e -1 2 z 2 1 2 H k (z), (B.12)where H k (z) are the probabilist's Hermite polynomials given byH k (z) = (-1) k e z 2 2 d k dz k e -z 2 2 . (B.13)To evaluate the integrals for this particular family, we can exploit the following recursions that are satisfied by Hermite polynomials:H k+1 (z) = zH k (z) -H ′ k (z), (B.14) H ′ k (z) = kH k-1 (z). (B.15)Eliminating the derivatives H ′ k (z) in Eqs. B.14-B.15, we see that zH k (z) = H k+1 (z) + kH k-1 (z). We can then substitute Eq. B.12 to obtain a recursion for the orthogonal basis functions themselves:zϕ k (z) = √ kϕ k+1 (z) + √ k-1ϕ k-1 (z). (B.16)With the above recursion, we can now read off these integrals from the property of orthogonality. For example, starting from Eq. B.2, we find thatµ ij = ∞ -∞ ϕ i (z)ϕ j (z) z dz, (B.17) = ∞ -∞ ϕ i (z) jϕ j+1 (z) + j -1ϕ j-1 (z) dz, (B.18) = δ i,j+1 j + δ i,j-1 √ i, (B.19)where δ ij is the Kronecker delta function. Next we consider the integral in Eq. B.3, which involves a power of z 2 in the integrand. In this case we can make repeated use of the recursion:ν ij = ∞ -∞ ϕ i (z)ϕ j (z) z 2 dz, (B.20) = ∞ -∞ √ iϕ i+1 (z) + √ i-1ϕ i-1 (z) jϕ j+1 (z) + j -1ϕ j-1 (z) dz, (B.21) = δ ij ij + (i-1)(j -1) + δ i-1,j+1 j(j +1) + δ j-1,i+1 i(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.1. Once these matrices are computed, it is a simple matter of substitution1 to compute the momentsE q [z 1 ], E q [z 2 1 ], and E q [z 1 z 2 ] from Eqs. B.8-B.9 and Eq. B.11. Finally, we can compute other low-order moments (such as E q [z 5 ] or E q [z 3 z 7 ]) by an appropriate permutation of indices.In this appendix we show in detail how the optimization for EigenVI reduces to a minimum eigenvalue problem. In particular we prove the following.Lemma C.1. Let {ϕ k (z)} ∞k=1 be an orthogonal function expansion, and let q ∈ Q K be the variational approximation parameterized byq(z) = K k=1 α k ϕ k (z) 2 , (C.1)where the weights satisfy K k=1 α 2 k = 1, thus ensuring that the distribution is normalized. Suppose furthermore that q is chosen to minimize the empirical estimate of the Fisher divergence given, as in eq. ( 10), byD π (q, p) = B b=1 q(z b ) π(z b ) ∇ log q(z b ) -∇ log p(z b ) 2 .Then the optimal variational approximation q in this family can be computed by solving the minimum eigenvalue problemmin q∈Q K D π (q, p) = min ∥α∥=1 α ⊤ M α =: λ min (M ), (C.2)where M is given in Eq. 13 and α = [α 1 , . . . , α K ] ∈ R K . The optimal weights α are given (up to an arbitrary sign) by the corresponding eigenvector of this minimal eigenvalue.Proof. The scores of q in this variational family are given by∇ log q(z b ) = 2 k α k ∇ϕ k (z b ) k α k ϕ k (z b ).Substituting the above into the empirical divergence, we find thatD π (q, p) = B b=1 q(z b ) π(z b ) ∇ log q(z b ) -∇ log p(z b ) 2 = B b=1 k α k ϕ k (z b ) 2 π(z b ) 2 k α k ∇ϕ k (z b ) k α k ϕ k (z b ) -∇ log p(z b ) 2 = B b=1 1 π(z b ) 2 k α k ∇ϕ k (z b ) - k α k ϕ k (z b ) ∇ log p(z b ) 2 = B b=1 1 π(z b ) k α k 2∇ϕ k (z b ) -ϕ k (z b )∇ log p(z b ) 2 = α ⊤ M α,where M is given in ( 13) and α = [α 1 , . . . , α K ] ∈ R K . Thus the optimal weights α are found by minimizing the quadratic form α ⊤ M α subject to the constraint α ⊤ α = 1. Equivalently, a solution can be found by minimizing the Rayleigh quotientargmin v v ⊤ M v v ⊤ v (C.3)and setting α = v/∥v∥. It then follows from the Rayleigh-Ritz theorem [8] for symmetric matrices that α is the eigenvector corresponding to the minimal eigenvalue of M , and this proves the lemma.Recall that EigenVI has two hyperparameters: the number of basis functions K and the number of importance samples B. We note there is an important difference between 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 in 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 is too small, it may take too long to converge in which case it may exceed computational budgets.Another fundamental difference in setting the number of basis functions 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 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.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 B 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 a set of scores, these can be reused to fit a larger number of basis functions. We compare the number of score evaluations wallclock vs 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, whereK = K 1 K 2 .The black star denotes the number of gradient evaluations for the Gaussian method.The experiments were run on a Linux workstation with a 32-core Intel(R) Xeon(R) w5-3435X processor and with 503 GB of memory. Experiments were run on CPU. In the sinh-arcsinh and posteriordb experiments, computations to construct the matrix M were parallelized over 28 threads.We considered the following synthetic 2D targets:• 3-component Gaussian mixture:p(z) = 0.4N (z | [-1, 1] ⊤ , Σ) + 0.3N (z | [1.1, 1.1] ⊤ , 0.5I) + 0.3N (z | [-1, -1] ⊤ , 0.5I),where we define Σ = 2 0.1 0.1 2 .• Funnel distribution with σ 2 = 1.2:p(z) = N (z 1 | 0, σ 2 )N (z 2 | 0, exp(z 1 /2)).• Cross distribution:p(z) = 1 4 N (z | [0, 2] ⊤ , Σ 1 ) + 1 4 N (z | [-2, 0] ⊤ , Σ 2 ) + 1 4 N (z | [2, 0] ⊤ , Σ 2 ) + 1 4 N (z | [0, -2] ⊤ , Σ 1 ),where Σ 1 = 0.15 0.9 0 0 1 and Σ 2 = 1 0 0 0.15 0.9 . These experiments were conducted without standardization with a Gaussian VI estimate. The EigenVI proposal distribution π used was a uniform([-9, 9]) distribution.In Figure E.2, we run EigenVI for increasing numbers of importance samples B and report the resulting forward KL divergence. The blue curves denote variational families with different K 1 = K 2 = K values used, i.e., 3, 6, and 10 (resulting in a total number of basis functions of 3 2 , 6 2 , and 10 2 ). In 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 λ t = BD t+1 . The black star denotes the number of score evaluations used by the Gaussian VI method.The sinh-arcsinh normal distribution [19, 20] has parameters s ∈ R D , τ ∈ R D + , Σ ∈ S ++ ; it is induced by transforming a Gaussian Z 0 ∼ N (0, Σ) to Z = S s,τ (Z 0 ), where S s,τ (z) := [S s1,τ1 (z 1 ), . . . , S s D ,τ D (z D )] ⊤ , S s d ,τ d (z d ) := sinh 1 τ d sinh -1 (z d ) + s d τ d . (E.1)Here s d controls the amount of skew in the dth dimension, and τ d controls the tail weight in that dimension. When s d = 0 and τ d = 1 in all dimensions d, the distribution is Gaussian.The sinh-arcsinh normal distribution has the following density:p(z; s, τ, Σ) = [(2π) D |Σ|] -1 2 D d=1 (1 + z 2 d ) -1 2 τ d C s d ,τ d (z d ) exp - 1 2 S s,τ (z) ⊤ Σ -1 S s,τ , (E.2)where we define the functions We constructed 3 targets in 2 dimensions and 3 targets in 5 dimensions, each with varying amounts of non-Gaussianity. The details of each target are below. In all experiments, EigenVI was applied with standardization, where a Gaussian was fit using batch and match VI with a batch size of 16 and a learning rate λ t = BD t+1 . For all experiments, we used a proposal distribution π that was uniform on [-5, 5] 2 . EQUATIONWe consider 8 real data targets from posteriordb, a suite of benchmark Bayesian models for real data problems. In Table E .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 R D .For all experiments, we fixed the number of importance samples to be B = 40,000; to construct the EigenVI 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 EigenVI. In most cases, the proposal distribution π was chosen to be uniform over [-6, 6] D . For the models 8-schools, which has a longer tail, we used a multivariate Gaussian proposal with zero mean and a scaled diagonal covariance σI, with σ = 3 2 .For the Gaussian score matching (GSM) method [37] , we chose a batch size of 16 for all experiments. We generally found the results were not too sensitive in comparison to other batch sizes of 4,8, and 32. For the batch and match (BaM) method [6] , we chose a batch size of 16. The learning rate was fixed at λ t = BD t+1 , which was a recommended schedule for non-Gaussian targets. For all ELBO optimization methods (full covariance Gaussian family and normalizing flow family), we used Adam to optimize the ELBO. We performed a 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, we used a real NVP [12] with 8 layers and 32 neurons. We found empirically that the computational of the scores was 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 EigenVI fit, a normalizing flow fit, and a 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 and the normalizing flow are able to model the tails of the distribution better than the Gaussian fit.EigenVI adds to the literature on BBVI, which has been an 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 to train generative models with malicious or unintended uses. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: 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 to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
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, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Justification: Yes we conform to the code of ethics. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 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.With further bookkeeping, one can also exploit the sparsity of µij and νij to derive more efficient calculations of these moments. 
 
 Appendix:  


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

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: 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?
- How is the proposal distribution chosen that is used for importance sampling?
limitations: limitations: The authors have adequately addressed the limitations of their work.
paper_summary: paper_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.
strength_weakness: strength_weakness: Weaknesses
* 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.

* 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.
questions: 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: limitations: See Weaknesses above.
paper_summary: paper_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.
strength_weakness: strength_weakness: See Questions.
questions: 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?
* Could the authors explain a bit more why centering the distribution is very important for better approximation?
* 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?
* 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$.
* Why is forward KL used in the whole paper rather than the commonly used backward KL as in ELBO?
* 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: limitations: /
paper_summary: paper_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.
strength_weakness: strength_weakness: Some weaknesses:

(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.

(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)?

(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.)

(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?

(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?
questions: questions: Clarifications:

(1) pp.6: What is the proposed advantage of EigenVI over Zhang et al.'s energy-based approach with a closed-form solution?

(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. 

(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?

(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. 

Typos:

(1) Abstract: "novel class variational approximations" should read "novel class of variational approximations".
limitations: limitations: Yes, authors describe limitations well.
paper_summary: paper_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.
strength_weakness: strength_weakness: * 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!)
* 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.
* 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?
* 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.

Overall, 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. 

## Additional Comments
* 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.
* Related works: A discussion on the literature of non-parametric/"less-parametric" VI algorithms seems necessary. 
* 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.
* 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.
* 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.

## References
1. Campbell, T., & Li, X. (2019). Universal boosting variational inference. NeurIPS
2. Smith, M. S., & Loaiza-Maya, R. (2023). Implicit copula variational inference. Journal of Computational and Graphical Statistics, 32(3), 769-781.
3. Han, S., Liao, X., Dunson, D., & Carin, L. (2016, May). Variational Gaussian copula inference. In AISTATS.
4. Shao, Y., Shan, N. Y., & Feng, T. (2024, April). Nonparametric Automatic Differentiation Variational Inference with Spline Approximation. In AISTATS.
5. Yu, L., & Zhang, C. (2023). Semi-Implicit Variational Inference via Score Matching. arXiv preprint arXiv:2308.10014.
6. Lambert, M., Chewi, S., Bach, F., Bonnabel, S., & Rigollet, P. (2022). Variational inference via Wasserstein gradient flows. NeurIPS.
7. Hyvärinen, A., & Dayan, P. (2005). Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4).
8. Ranganath, R., Gerrish, S., & Blei, D. (2014, April). Black box variational inference. In AISTATS.
questions: 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?
* Eq (10): Why is importance sampling used here when one can sample directly from $q$?
limitations: limitations: The paper would benefit from a dedicated limitation section.

Comments: title: Re: Rebuttal by Authors
Response
Rebuttal?
comment: Thank you authors for your detailed response.  After reading the rebuttal and other reviews, I continue to find the paper both technically solid and innovative, and will maintain my positive score.
Thanks for addressing my questions!

After reading all reviews and corresponding rebuttals, I will maintain my positive score. I believe this is a solid paper!
Thank you for the follow up &ndash; we now understand your question. 

Developing algorithms that support subsampling is an important line of future work, and we believe EigenVI develops some foundations that can be adapted to support subsampling. This may be achieved, for instance, with an iterative approach for updating variational parameters $\alpha_t$ that uses the score matching objective, where each iteration would involve solving for the minimum eigenvalue of a particular matrix&mdash;here each iteration could also take in scores from a subsampled batch of data points. However, a proper exploration of such a method warrants its own dedicated paper to explore a new algorithm (even before incorporating subsampling) and its properties, and we are planning to investigate these directions in future work. In the current work, we have maintained focus on introducing a novel variational family and non-iterative VI algorithm that we believe has its own merits.

We think the question of subsampling will be of broad interest to other readers, and we will add a discussion to this end in the discussion & future work section of the revised paper. We thank the reviewer for bringing up this interesting point.
I thank the author for the detailed response. Despite this, I think the weaknesses of this paper in my review still remain. Therefore, I would keep this score.
Thank you for the detailed response.

For the comment on unbiasedness, I should have phrased my original comment a little better. It is well known the score-matching objective itself is valid with subsampling. I was more concerned about the fact that it is unclear how to incorporate subsampling within the proposed algorithm in a consistent way.

Nevertheless, I am convinced that the proposed methodology is technically interesting and will raise my score. However, I am still concerned that there isn't enough evidence that it would be practically applicable to larger-scale real-world problems. Thus, the paper feels somewhat less complete in its current state. For instance:

> There 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.

Demonstrating this would have been much more convincing and made the work more complete.
Thanks for the responses. I will keep my positive score.
Dear authors (and reviewers),

There were some positive reviews, but also one negative and one borderline review.  At the end of the day, we don't simply average scores, so this paper is not guaranteed acceptance, though it certainly has a good chance.

Hence I'd advise the authors to write rebuttals, at least to the two most negative/neutral reviews.

Reviewers, please take a look at each others' reviews and see if that changes your mind. As a reminder, the OpenReview platform allows us to have a dialogue, so you can post small questions to the authors and authors can respond.

Best regards

Area Chair