orangesize= tiny, tickmarkheight=0.25em ,shadow tealsize= tiny,
tickmarkheight=0.25em ,shadow, color=cyan

# An Ordering of Divergences for Variational Inference
with Factorized Gaussian Approximations

Charles C. Margossian1, Loucas Pillaud-Vivien2, and Lawrence K. Saul1
1Center for Computational Mathematics, Flatiron Institute
2CERMICS Laboratory, Ecole des Ponts ParisTech

###### Abstract

Given an intractable distribution $p$, the problem of variational inference
(VI) is to compute the best approximation $q$ from some more tractable family
$\mathcal{Q}$. Most commonly the approximation is found by minimizing a
Kullback-Leibler (KL) divergence. However, there exist other valid choices of
divergences, and when $\mathcal{Q}$ does not contain $p$, each divergence
champions a different solution. We analyze how the choice of divergence
affects the outcome of VI when a Gaussian with a dense covariance matrix is
approximated by a Gaussian with a diagonal covariance matrix. In this setting
we show that different divergences can be ordered by the amount that their
variational approximations misestimate various measures of uncertainty, such
as the variance, precision, and entropy. We also derive an impossibility
theorem showing that no two of these measures can be simultaneously matched by
a factorized approximation; hence, the choice of divergence informs which
measure, if any, is correctly estimated. Our analysis covers the KL
divergence, the Rényi divergences, and a score-based divergence that compares
$\nabla\log p$ and $\nabla\log q$. We empirically evaluate whether these
orderings hold when VI is used to approximate non-Gaussian distributions.

Keywords: variational inference, Kullback-Leibler divergence, Rényi
divergence, Fisher divergence, score matching

##  1 Introduction

In many statistical analyses, it is necessary to approximate an intractable
distribution, $p$, by a more tractable distribution, $q$. This problem can be
posed as an optimization, one whose goal is to discover the best-matching
distribution $q$ within some larger, parameterized family $\mathcal{Q}$ of
approximating distributions. Optimizations of this sort are the subject of a
large literature on variational inference (VI; Jordan et al., 1999; Wainwright
and Jordan, 2008; Blei et al., 2017). VI is not the only way to approximate an
intractable distribution, but with a finite computational budget, it often
achieves a smaller error than Markov chain Monte Carlo methods, even if the
latter possess better asymptotic properties (Robert and Casella, 2004). In
this sense, VI is well suited to the challenges posed by large data sets and
high-dimensional models. It is now widely applied to problems in Bayesian
inference and machine learning.

Despite its widespread adoption, many aspects of VI remain poorly or
incompletely understood. At the heart of VI is the question of what
constitutes a useful approximation of $p$. Suppose we wish to approximate a
target distribution $p$ over latent variables ${\bf z}=z_{1:n}\in\mathcal{Z}$,
and let $\mathcal{Q}$ be the family of approximating distributions. VI
attempts to solve an optimization of the form

$q^{*}=\underset{q\in\mathcal{Q}}{\text{argmin}}\ D(q,p),$ (1)

where $D$ is a divergence satisfying (i) $D(q,p)\geq 0$ for all
$q\in\mathcal{Q}$ and (ii) $D(q,p)=0$ if and only if $p=q$.

In principle any proper divergence $D(q,p)$ may be used for VI in eq. (1). By
far the most popular choice is the Kullback-Leibler divergence,
$\text{KL}(q||p)$, which measures the disagreement between $\log p$ and $\log
q$, but it is not the only one that can be considered. Other valid choices
include the “forward” $\text{KL}(p||q)$, which arises in expectation
propagation (Vehtari et al., 2020), the Rényi or $\alpha$-divergences
$\text{R}_{\alpha}(p||q)$ (Minka, 2005; Li and Turner, 2016), which
interpolate between $\text{KL}(q||p)$ and $\text{KL}(p||q)$ when
$\alpha\\!\in\\!(0,1)$, and score-based (or Fisher) divergences (Courtade,
2016; Cai et al., 2024), which measure the discrepancy between $\nabla\log p$
and $\nabla\log q$. Table 1 provides an overview of the divergences we analyze
in this paper.

Divergence |  Notation | Definition
---|---|---
(Reverse) Kullback-Leibler | $\text{KL}(q||p)$ | $\mathbb{E}_{q}\\!\left[\log\frac{q}{p}\right]$
(Forward) Kullback-Lebiler | $\text{KL}(p||q)$ | $\mathbb{E}_{p}\\!\left[\log\frac{p}{q}\right]$
Rényi of order $\alpha\in(0,1)$ | $\text{R}_{\alpha}(p||q)$ | $\frac{1}{\alpha(\alpha-1)}\,\mathbb{E}_{q}\\!\left[\frac{p^{\alpha}}{q^{\alpha}}-1\right]$
(Reverse) Score-based | $\text{S}(q||p)$ | $\mathbb{E}_{q}\\!\left[\big{\|}\nabla\log\frac{q}{p}\big{\|}^{2}_{\text{Cov}(q)}\right]$
(Forward) Score-based | $\text{S}(p||q)$ | $\mathbb{E}_{p}\\!\left[\big{\|}\nabla\log\frac{q}{p}\big{\|}^{2}_{\text{Cov}(p)}\right]$
Table 1: Divergences that we consider in this paper for VI. The score-based
divergences (Cai et al., 2024) are defined here for the special case that the
base distribution is Gaussian.

All of these divergences vanish if and only if $q\\!=\\!p$, but when $p\notin
Q$, as is typically the case in VI, the solution of eq. (1) can vary
dramatically from one divergence to the other. We illustrate this when $p$ is
a two-dimensional Gaussian with correlation $\varepsilon=0.75$ and
$\mathcal{Q}$ is the family of Gaussians with diagonal covariance matrices;
see Fig. 1. The goal of this paper is to understand which divergence works
best for a given application.

a. $\mathcal{Q}$$p$ b.
Divergence$\color[rgb]{0.69921875,0.09375,0.16796875}\definecolor[named]{pgfstrokecolor}{rgb}{0.69921875,0.09375,0.16796875}{\text{S}(p||q)}$$\color[rgb]{0.9375,0.54296875,0.3828125}\definecolor[named]{pgfstrokecolor}{rgb}{0.9375,0.54296875,0.3828125}{\text{KL}(p||q)}$$\color[rgb]{0.9921875,0.859375,0.78125}\definecolor[named]{pgfstrokecolor}{rgb}{0.9921875,0.859375,0.78125}{\text{R}_{\alpha}(q||p)}$$\color[rgb]{0.8203125,0.8984375,0.94140625}\definecolor[named]{pgfstrokecolor}{rgb}{0.8203125,0.8984375,0.94140625}{\text{KL}(q||p)}$$\color[rgb]{0.40234375,0.6640625,0.8125}\definecolor[named]{pgfstrokecolor}{rgb}{0.40234375,0.6640625,0.8125}{\text{S}(q||p)}$MatchEntropyVariancePrecisionTarget
$p$

Figure 1:  (a) The optimal projection of $p$ onto $\mathcal{Q}$ depends on the
divergence used to measure disagreement between $p$ and
$q\\!\in\\!\mathcal{Q}$. (b) Optimal factorized approximations $q$ to a two-
dimensional Gaussian $p$ with off-diagonal correlation, as measured by
different divergences. The marginal variances are matched by minimizing
$\text{KL}(p||q)$, the marginal precisions by minimizing $\text{KL}(q||p)$,
and the entropy (in this example) by minimizing $\text{R}_{\alpha}(p||q)$ with
$\alpha\\!=\\!0.5$.

For VI to be useful, these divergences must align with inferential goals.
These goals may vary from one application to the next. In Bayesian statistics,
the aim is often to calculate the posterior mean, variance, and quantiles of
interpretable univariate quantities (Gelman et al., 2013); in certain
applications, it may also or alternatively be of interest to calculate the
posterior precision (or inverse covariance) matrix (Bernardo and Smith, 2000).
In statistical physics, the aim is generally to compute the free energy
(Parisi, 1988), which is itself equal to the difference of the internal energy
and the entropy. In this paper, we consider three distinct measures of
uncertainty—(i) the marginal variances, (ii) the marginal precisions, and
(iii) the multivariate entropy—and for each of these measures, we analyze how
the fidelity of VI depends on the divergence that it attempts to minimize.

Our work is in line with several efforts to evaluate VI through the lens of
downstream inferential tasks. The mean and variances of $p$ are of most
interest for Bayesian statistics, but Huggins et al. (2020) caution that a low
$\text{KL}(q||p)$ or $R_{\alpha}(p||q)$ divergence does not guarantee that $q$
accurately estimates these quantities. Instead they recommend to estimate the
Wasserstein distance, which can be used to bound the difference in moments of
$q$ and $p$; see also Biswas and Mackey (2024). In latent variable models, VI
is used to bound a marginal likelihood which is then maximized with respect to
the model parameters (e.g., the weights of a neural network (Tomczak, 2022)).
Several works examine the tightness of this bound (e.g Li and Turner, 2016;
Dieng et al., 2017; Daudel et al., 2023) and the statistical properties of the
resulting maximum likelihood estimator (Wang and Blei, 2018). VI is also used
in tandem with other inference schemes—for example, to initialize Markov chain
Monte Carlo (Zhang et al., 2022) or to estimate a proposal distribution for
importance sampling. For the latter, the quality of VI can be assessed from
the distribution of importance weights (Yao et al., 2018; Vehtari et al.,
2024). In sum, there are many different ways to evaluate VI, and one lesson of
our paper is that different inferential goals may be incompatible—may, in
fact, compete against one another.

In most applications of VI, the variational family does not include the target
distribution (i.e., $p\not\in\mathcal{Q}$), and the resulting inference must
be compromised in some way. The simplest variational families are those that
only contain factorized distributions,

$q({\bf z})=\prod_{i=1}^{n}q_{i}(z_{i}).$ (2)

Variational families of this form originated in the mean-field approximations
for physical models of ferromagnetism (Parisi, 1988). In practice we do not
expect $p$ to factorize; the assumption is made out of computational
convenience. When ${\bf z}\in\mathbb{R}^{n}$, a common choice for $q_{i}$ is
that of a Gaussian density (e.g Kucukelbir et al., 2017). The resulting
algorithm is called factorized Gaussian VI (FG-VI). By its very nature, FG-VI
cannot estimate the correlations between the different components of ${\bf
z}$. Still, at least in theory, it is possible to construct a factorized
Gaussian that matches the variance, precision, or entropy of $p$ and, in this
sense, construct an approximation that (however misleading it may be in other
respects) achieves its primary inferential goal to estimate one of these
quantities.

In this paper, we analyze how the choice of divergence in eq. (1) affects the
ability of FG-VI to meet these goals. We focus on the case where $p$ is a
Gaussian distribution with a non-diagonal covariance matrix. This setting,
though simple, is already rich enough to provide several counter-intuitive
results, and it reveals in a stark way that different inferential goals are
best served by different divergences. Our analysis is a natural first step in
the study of the inherent trade-offs of factorized VI. There has been previous
work (MacKay, 2003; Turner and Sahani, 2011; Margossian and Saul, 2023) in the
setting where $p$ is Gaussian, but here we provide a more thorough treatment,
formalizing existing intuitions and extending many earlier results. Our
contributions are as follows:

  * •

In Section 2, we present an impossibility theorem: for any choice of the
divergence in eq. (1), the resulting approximation $q$ from FG-VI cannot
simultaneously match any two of the variance, precision, and entropy of the
target distribution $p$, and hence the choice of divergence determines which
(if any) of these measures of uncertainty is correctly estimated.

  * •

In Section 3, we study the optimizations that arise in FG-VI when $p$ is
Gaussian, and we analyze the solutions to these optimizations for all of the
divergences that appear in Table 1. It has been shown previously that FG-VI
matches the marginal precisions of $p$ when minimizing $\text{KL}(q||p)$ and
the marginal variances when minimizing $\text{KL}(p||q)$ (MacKay, 2003; Turner
and Sahani, 2011); we review these results in light of the impossibility
theorem of Section 2. Later, with the additional machinery of Section 4, we
show that there always exists a unique $\alpha\\!\in\\!(0,1)$, determined by
the covariance of $p$, such that FG-VI matches the entropy of $p$ when
minimizing the Rényi divergence $\text{R}_{\alpha}(p||q)$. In this section, we
also show that FG-VI can wildly misestimate the marginal variances when
minimizing the score-based divergences; in particular, we derive conditions
under which it estimates a vanishing variance when minimizing $\text{S}(q||p)$
and an infinite variance when minimizing $\text{S}(p||q)$. We refer to such
cases of FG-VI as instances of variational collapse, where the variational
approximation is not well-defined because
$\arg\\!\inf_{q\in\mathcal{Q}}D(q,p)\not\in\mathcal{Q}$.

  * •

Our main result, presented in Section 4, is an ordering of all the divergences
in Table 1 based on the estimates they yield of the marginal variances. This
ordering also induces corresponding orderings based on the marginal precisions
and the entropy.

  * •

In Section 5, we experiment with FG-VI on several data sets and investigate
numerically whether these orderings hold for target distributions that are not
Gaussian.

##  2 Impossibility theorem for FG-VI

First we explore various trade-offs when one normal distribution, $p$, whose
covariance matrix has off-diagonal elements, is approximated by another normal
distribution, $q$, whose covariance matrix is diagonal. Let

$\displaystyle p$
$\displaystyle=\text{normal}(\boldsymbol{\mu},\boldsymbol{\Sigma}),$ (3)
$\displaystyle q$
$\displaystyle=\text{normal}(\boldsymbol{\nu},\boldsymbol{\Psi}),$ (4)

where $\boldsymbol{\Psi}$ is diagonal but $\boldsymbol{\Sigma}$ is not; here
we have in mind the setting of FG-VI where $p$ is the target distribution and
$q$ is the variational approximation. Nevertheless we emphasize that the
results in this section are not tied to any particular choice of the
divergence $D(q,p)$ in eq. (1). Rather, these trade-offs are inherent to the
use of a factorized approximation, and they arise regardless of how $q$ is
chosen.

The main result of this section is an impossibility theorem: the theorem
considers three measures of $p$’s uncertainty—its marginal variances, its
marginal precisions, and its entropy—and states that if $p$ has a non-diagonal
covariance matrix, then no factorized approximation $q$ can correctly estimate
any two of these measures at the same time. To state the theorem, the
following notation is useful. We denote the entropy of a distribution by

$\mathcal{H}(p)=-\mathbb{E}_{p}[\log p],$ (5)

and we remind the reader that the precision matrix is given by the inverse of
the covariance matrix. Also, for any square matrix, we use
$\text{diag}(\cdot)$ to denote the vector formed from its diagonal elements,
and we use vector inequalities, such as
$\text{diag}(\boldsymbol{\Psi})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma})$, to
denote that $\Psi_{ii}\\!\leq\\!\Sigma_{ii}$ for all $i$. With this notation,
we have the following result.

######  Theorem 2.1 (Impossibility theorem for FG-VI).

Let $p$ and $q$ be normal distributions over $\mathbb{R}^{N}$ with covariances
$\boldsymbol{\Sigma}$ and $\boldsymbol{\Psi}$, respectively, where
$\boldsymbol{\Psi}$ is diagonal but $\boldsymbol{\Sigma}$ is not. Then 1.
(Variance matching)
If
$\text{diag}({\boldsymbol{\Psi}})\\!=\\!\text{diag}({\boldsymbol{\Sigma}})$,
then $\mathcal{H}(q)\\!>\\!\mathcal{H}(p)$ and
$\text{diag}(\boldsymbol{\Psi}^{-1})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma}^{-1})$,
and this last inequality is strict for at least one component along the
diagonal. 2. (Precision matching)
If
$\text{diag}({\boldsymbol{\Psi}}^{-1})\\!=\\!\text{diag}({\boldsymbol{\Sigma}^{-1}})$,
then $\mathcal{H}(q)\\!<\\!\mathcal{H}(p)$ and
$\text{diag}(\boldsymbol{\Psi})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma})$,
and this last inequality is strict for at least one component along the
diagonal. 3. (Entropy matching)
If $\mathcal{H}(q)\\!=\\!\mathcal{H}(p)$, then $\Psi_{ii}\\!<\\!\Sigma_{ii}$
and $\Psi^{-1}_{jj}\\!<\\!\Sigma^{-1}_{jj}$ for at least one $i$ and $j$.

The theorem states that when $q$ matches $p$ in any one of the above ways, it
cannot also match $p$ in any of the others (except in the trivial case that
$\boldsymbol{\Sigma}$ is diagonal). We prove the theorem with three lemmas,
each of which examines the trade-off in FG-VI between two different measures
of uncertainty.

###  2.1 Variance-entropy tradeoff

The trade-off between variance and entropy has been previously studied in the
classical setting of FG-VI where $q$ is chosen to minimize $\text{KL}(q||p)$.
There it was noted that the entropy of $q$ is governed by a shrinkage-
delinkage trade-off (Margossian and Saul, 2023): the entropy of $q$ is
decreased when $q$ underestimates the marginal variances of $p$ (the shrinkage
effect), but it is increased when $q$ ignores the correlations in $p$ between
off-diagonal components (the delinkage effect). We use similar tools to prove
the following result.

######  Lemma 2.2 (Variance-entropy trade-off).

If $p$ in eq. (3) has a non-diagonal covariance, then $q$ in eq. (4) with a
diagonal covariance cannot correctly match both its marginal variances and
entropy: i.e., if
$\text{diag}(\boldsymbol{\Psi})\\!=\\!\text{diag}(\boldsymbol{\Sigma})$, then
${\cal H}(q)\\!>\\!{\cal H}(p)$, and if ${\cal H}(q)\\!=\\!{\cal H}(p)$, then
$\Psi_{ii}\\!<\\!\Sigma_{ii}$ for some $i$.

###### Proof.

We begin by defining two auxiliary matrices whose elements are dimensionless
(unlike those of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Psi}$). The first of
these is the correlation matrix ${\bf C}\in\mathbb{R}^{n\times n}$ of the
target distribution $p$, with elements

$C_{ij}=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}},$ (6)

which is positive-definite and has unit diagonal entries, $C_{ii}\\!=\\!1$.
The second is the diagonal matrix ${\bf R}\in\mathbb{R}^{n\times n}$ that
records the ratios of marginal variances in $q$ and $p$; i.e.,

$R_{ii}=\frac{\Psi_{ii}}{\Sigma_{ii}}.$ (7)

Next we consider the difference in entropy of $p$ and $q$. From basic
properties of Gaussian distributions (Cover and Thomas, 2006) and the
definitions in eqs. (6–7), it follows that

$\displaystyle\mathcal{H}(p)-\mathcal{H}(q)$
$\displaystyle=\tfrac{1}{2}\log\Big{[}|\boldsymbol{\Sigma}|(2\pi
e)^{n}\Big{]}-\tfrac{1}{2}\log\Big{[}|\boldsymbol{\Psi}|(2\pi e)^{n}\Big{]},$
$\displaystyle=\tfrac{1}{2}\log|\boldsymbol{\Sigma}|-\tfrac{1}{2}\log|\boldsymbol{\Psi}|,$
$\displaystyle=\tfrac{1}{2}\log\left|\boldsymbol{\Psi}^{-\frac{1}{2}}\boldsymbol{\Sigma}\boldsymbol{\Psi}^{-\frac{1}{2}}\right|,$
$\displaystyle=\tfrac{1}{2}\log\left|{\bf R}^{-\frac{1}{2}}{\bf C}{\bf
R}^{-\frac{1}{2}}\right|,$ (8) $\displaystyle=\tfrac{1}{2}\log|{\bf
R}|^{-1}-\tfrac{1}{2}\log|{\bf C}|^{-1}.$ (9)

The first determinant in eq. (9) is easy to compute because the matrix ${\bf
R}$, as defined in eq. (7), is purely diagonal. Notably, we can bound the
second determinant in eq. (9) by the Hadamard inequality (Horn and Johnson,
2012, Theorem 7.8.1), which states that

$|{\bf C}|<\prod_{i=1}^{n}C_{ii}=1.$ (10)

We note that the inequality in eq. (10) is strict under the assumption of the
lemma that $\boldsymbol{\Sigma}$ (and hence also ${\bf C}$) is not diagonal.
Combining the last two results, we see that

$\mathcal{H}(p)-\mathcal{H}(q)\ <\
\frac{1}{2}\sum_{i}\log\left(\frac{\Sigma_{ii}}{\Psi_{ii}}\right).$ (11)

Now suppose that
$\text{diag}(\boldsymbol{\Psi})\\!=\\!\text{diag}(\boldsymbol{\Sigma})$; then
from eq. (11) we see that $\mathcal{H}(p)\\!<\\!\mathcal{H}(q)$, proving the
first claim of the lemma. Conversely, suppose that
$\mathcal{H}(p)\\!=\\!\mathcal{H}(q)$; then from eq. (11) we see that
$\sum_{i}\log\Sigma_{ii}\\!>\\!\sum_{i}\log\Psi_{ii}$, and this can only be
true if $\Sigma_{ii}\\!>\\!\Psi_{ii}$ for some $i$. This proves the second
claim of the lemma. ∎

###  2.2 Precision-entropy trade-off

Next we derive a similar trade-off between precision and entropy. This trade-
off is captured by the following lemma.

######  Lemma 2.3 (Precision-entropy trade-off).

If $p$ in eq. (3) has a non-diagonal covariance, then $q$ in eq. (4) with a
diagonal covariance cannot correctly match both its marginal precisions and
entropy: i.e., if
$\mathrm{diag}(\boldsymbol{\Psi}^{-1})\\!=\\!\mathrm{diag}(\boldsymbol{\Sigma}^{-1})$,
then ${\cal H}(q)\\!<\\!{\cal H}(p)$, and if ${\cal H}(q)\\!=\\!{\cal H}(p)$,
then $\Psi^{-1}_{ii}\\!<\\!\Sigma^{-1}_{ii}$ for some $i$.

###### Proof.

The proof follows the same structure as the proof of the preceding lemma. We
define auxiliary matrices $\tilde{{\bf C}}$ and $\tilde{{\bf R}}$ with the
dimensionless elements

$\tilde{C}_{ij}=\frac{\Sigma^{-1}_{ij}}{\sqrt{\Sigma^{-1}_{ii}\Sigma^{-1}_{jj}}},\quad\quad\tilde{R}_{ii}=\frac{\Psi_{ii}^{-1}}{\Sigma^{-1}_{ii}},$
(12)

where $\tilde{{\bf C}}$ has all ones along the diagonal and $\tilde{{\bf R}}$
has all zeroes off the diagonal. Again we examine the difference in entropy
between $p$ and $q$. Proceeding from eq. (8), we find that

$\mathcal{H}(p)-\mathcal{H}(q)=\tfrac{1}{2}\log\left|{\bf
R}^{-\frac{1}{2}}{\bf C}{\bf
R}^{-\frac{1}{2}}\right|=\tfrac{1}{2}\log\left|\tilde{{\bf
R}}^{\frac{1}{2}}\tilde{{\bf C}}^{-1}\tilde{{\bf
R}}^{\frac{1}{2}}\right|=\tfrac{1}{2}\log|\tilde{{\bf
R}}|-\tfrac{1}{2}\log|\tilde{{\bf C}}|.$ (13)

As before, the difference in entropy is expressed in terms of two
determinants. The determinant of $\tilde{{\bf R}}$ is easy to compute because
$\tilde{{\bf R}}$, as defined in eq. (12), is a diagonal matrix. Again we
deduce from the Hadamard inequality that $|\tilde{{\bf
C}}|<\prod_{i}\tilde{{\bf C}}_{ii}=1$. Thus we see that

$\mathcal{H}(p)-\mathcal{H}(q)\ >\
\frac{1}{2}\sum_{i}\log\left(\frac{\Psi_{ii}^{-1}}{\Sigma^{-1}_{ii}}\right).$
(14)

Note that the direction of inequality in eq. (14) is reversed from the
direction of inequality in eq. (11). Now suppose that
$\text{diag}(\boldsymbol{\Psi}^{-1})\\!=\\!\text{diag}(\boldsymbol{\Sigma}^{-1})$;
then from eq. (14) we see that $\mathcal{H}(p)\\!>\\!\mathcal{H}(q)$, proving
the first claim of the lemma. Conversely, suppose that
$\mathcal{H}(p)\\!=\\!\mathcal{H}(q)$; then from eq. (14) we see that
$\sum_{i}\log\Sigma^{-1}_{ii}\\!>\\!\sum_{i}\log\Psi^{-1}_{ii}$, and this can
only be true if $\Sigma^{-1}_{ii}\\!>\\!\Psi^{-1}_{ii}$ for some $i$. This
proves the second claim of the lemma. ∎

###  2.3 Variance-precision trade-off

Finally, we derive a similar trade-off between variance and precision. This
trade-off is captured by the following lemma.

######  Lemma 2.4 (Variance-precision trade-off).

If $p$ in eq. (3) has a non-diagonal covariance, then $q$ in eq. (4) with a
diagonal covariance cannot correctly match both its marginal variances and
precisions: i.e., if
$\text{diag}(\boldsymbol{\Psi})\\!=\\!\text{diag}(\boldsymbol{\Sigma})$, then
$\text{diag}(\boldsymbol{\Psi}^{-1})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma}^{-1})$
and this inequality is strict for at least one component along the diagonal;
likewise, if
$\text{diag}(\boldsymbol{\Psi}^{-1})\\!=\\!\text{diag}(\boldsymbol{\Sigma}^{-1})$,
then
$\text{diag}(\boldsymbol{\Psi})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma})$ and
this inequalitiy is strict for at least one component along the diagonal.

###### Proof.

Let ${\bf C}$ denote the correlation matrix in eq. (6), and let
$\lambda_{1},\lambda_{2},\ldots,\lambda_{n}$ denote its eigenvalues. Also, let
$\mathbf{e}_{i}$ denote the unit vector along the $i^{\textrm{th}}$ axis.
Since $C_{ii}\\!=\\!1$, we have

$\displaystyle C^{-1}_{ii}$ $\displaystyle=\mathbf{e}_{i}^{\top}({\bf
C}^{-1}\\!+{\bf C})\,\mathbf{e}_{i}-1,$ (15)
$\displaystyle\geq\min_{\|\mathbf{e}\|=1}\left[\mathbf{e}^{\top}({\bf
C}^{-1}\\!+{\bf C})\,\mathbf{e}-1\right],$ (16)
$\displaystyle=\min_{j}\left(\lambda_{j}^{-1}+\lambda_{j}-1\right),$ (17)
$\displaystyle\geq\min_{\lambda>0}\left(\lambda^{-1}+\lambda-1\right),$ (18)
$\displaystyle=1,$ (19)

where in eq. (17) we have exploited the properties of Rayleigh quotients.
Next, starting from eq. (6), we can verify by direct calculation that
$C^{-1}_{ij}=\Sigma^{-1}_{ij}\sqrt{\Sigma_{ii}\Sigma_{jj}}$, and in particular
along the diagonal (when $i\\!=\\!j$), we have

$C^{-1}_{ii}=\Sigma^{-1}_{ii}\Sigma_{ii}.$ (20)

Now suppose that the variances of $p$ and $q$ match; then combining eq. (20)
and the previous inequality, we see that

$1\leq
C^{-1}_{ii}=\Sigma^{-1}_{ii}\Sigma_{ii}=\Sigma^{-1}_{ii}\Psi_{ii}=\frac{\Sigma^{-1}_{ii}}{\Psi^{-1}_{ii}},$
(21)

and we conclude (since $i$ was arbitrary) that
$\text{diag}(\boldsymbol{\Psi}^{-1})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma}^{-1})$.
Conversely, suppose that the precisions of $p$ and $q$ match; then by the same
reasoning we see that

$1\leq
C^{-1}_{ii}=\Sigma^{-1}_{ii}\Sigma_{ii}=\Psi^{-1}_{ii}\Sigma_{ii}=\frac{\Sigma_{ii}}{\Psi_{ii}},$
(22)

and again we conclude (since $i$ was arbitrary) that
$\text{diag}(\boldsymbol{\Psi})\\!\leq\\!\text{diag}(\boldsymbol{\Sigma})$.
Finally, we recall that by assumption $\boldsymbol{\Sigma}$ (and hence ${\bf
C}$) is not diagonal; furthermore, if $C_{ij}\neq 0$, then the unit vector
$\mathbf{e}_{j}$ is not an eigenvector of ${\bf C}$. In this case, eq. (16)
becomes a strict inequality, and for the $j^{\text{th}}$ component along the
diagonal we see that $C_{jj}^{-1}\\!>\\!1$. But then eqs. (21–22) also become
strict inequalities, showing that $\Psi^{-1}_{jj}\\!<\\!\Sigma^{-1}_{jj}$ when
the variances match and $\Psi_{jj}\\!<\\!\Sigma_{jj}$ when the precisions
match. This proves the lemma. ∎

######  Remark 2.5.

The proof of Lemma 2.4 does not use the fact that $p$ and $q$ are Gaussian,
and so the result extends to any factorized distribution $q$ and non-
factorized distribution $p$ with finite second moments.

From the preceding lemma, we can also extract a more general result of linear
algebra (Horn and Johnson, 2012, Theorem 7.17) that will prove especially
useful in Section 4. Note that for any covariance matrix $\boldsymbol{\Sigma}$
and correlation matrix ${\bf C}$, we have shown above that
$\text{diag}(\mathbf{C}^{-1})\\!\geq\\!1$; also, this inequality is strict for
at least one component along the diagonal when $\boldsymbol{\Sigma}$ is not
itself diagonal. Since $C^{-1}_{ii}\\!=\\!\Sigma_{ii}\Sigma^{-1}_{ii}$, we
have therefore shown the following.

######  Proposition 2.6.

Let $\mathbf{A}\in\mathbb{R}^{n\times n}$. If
$\mathbf{A}\\!\succ\\!\mathbf{0}$, then $\displaystyle A_{ii}A^{-1}_{ii}$
$\displaystyle\geq 1\quad\mbox{for all $i$},$ (23) $\displaystyle
A_{jj}A^{-1}_{jj}$ $\displaystyle>1\quad\mbox{for some $j$};$ (24) in
particular, the inequality is strict along the $j^{\text{th}}$ component
whenever $\sum_{k\neq j}|A_{jk}|>0$.

To conclude the section, we note the following: we obtain the Impossibility
Theorem for FG-VI merely by piecing together the results of the preceding
lemmas. In particular, Lemma 2.2 proves the first claims of Theorem 2.1 under
the conditions for variance and entropy matching, Lemma 2.4 proves the second
claims of Theorem 2.1 under the conditions for variance and precision
matching, and Lemma 2.3 proves the remaining claims.

##  3 Divergences for FG-VI

In this section we derive the solutions for FG-VI that are obtained from eq.
(1) by minimizing the divergences in Table 1. We also highlight various
properties of these solutions. As in the previous section, we assume
everywhere that $p$ and $q$ are given by eqs. (3–4), where $\boldsymbol{\Psi}$
is diagonal but $\boldsymbol{\Sigma}$ is not.

###  3.1 KL divergence

Most VI is based on minimizing the reverse KL divergence in Table 1, and the
implications of this choice for FG-VI have been extensively studied. For
completeness, we review certain key results in the notation of this paper. For
$p$ and $q$ in eqs. (3–4), the reverse KL divergence is given by (Cover and
Thomas, 2006)

$\text{KL}(q||p)=\frac{1}{2}\left[(\boldsymbol{\nu}\\!-\\!\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\nu}\\!-\\!\boldsymbol{\mu})+\log\frac{|\boldsymbol{\Sigma}|}{|\boldsymbol{\Psi}|}+\textrm{Tr}(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Psi})-n\right].$
(25)

With respect to this divergence, the best factorized approximation $q$ is
obtained by choosing the parameters $\boldsymbol{\nu}$ and $\boldsymbol{\Psi}$
that minimize eq. (25). The variational mean appears only in the first term of
eq. (25), which is minimized by setting $\boldsymbol{\nu}=\boldsymbol{\mu}$,
and the remaining terms are minimized by setting

$\textrm{diag}(\boldsymbol{\Psi}^{-1})=\textrm{diag}(\boldsymbol{\Sigma}^{-1}).$
(26)

Further details of these calculations can be found in the references (Turner
and Sahani, 2011; Margossian and Saul, 2023).

It has been widely observed that VI based on the reverse KL divergence tends
to underestimate uncertainty (MacKay, 2003; Minka, 2005; Turner and Sahani,
2011; Blei et al., 2017; Giordano et al., 2018). We can now refine and qualify
this statement in the setting of FG-VI where both $p$ and $q$ are Gaussian. In
this case we have the following:

######  Proposition 3.1 (Precision matching).

If $q$ in eq. (4) has a diagonal covariance but $p$ in eq. (3) does not, then
the reverse KL divergence in eq. (25) is minimized when $q$ matches the mean
and marginal precisions of $p$ but strictly underestimates the joint entropy
and underestimates all of the marginal variances (at least one of them
strictly). The proof is immediate: the matching of precisions follows from
eqs. (26), and the underestimation of variances and entropy follows from the
Impossibility Theorem in Section 2. This is the first of several results
showing that the “uncertainty deficit” of FG-VI depends on the measure that is
used to quantify uncertainty. Indeed, when $p$ is Gaussian, the marginal
precisions are correctly estimated by minimizing the reverse KL divergence;
moreover, it has been shown that when $\boldsymbol{\Sigma}$ has constant off-
diagonal entries, the entropy gap between $p$ and $q$ only grows as
$\mathcal{O}(\log n)$, while the entropy itself grows as $\mathcal{O}(n)$,
meaning that the fractional entropy gap tends to 0 as $n\\!\to\\!\infty$
(Margossian and Saul, 2023, Theorem 3.6). A similar result is obtained in the
thermodynamic limit of Ising models with constant, long-range interactions,
where the mean-field approximation has its roots (Parisi, 1988).

Next we consider the forward KL divergence in Table 1. This divergence is not
generally minimized for VI because it involves an expectation, namely
$\mathbb{E}_{p}[\log\frac{p}{q}]$ with respect to the target distribution $p$.
When $p$ is not tractable, it is not possible to compute this expectation
analytically, and it may be too expensive to estimate it by sampling. However,
it is sometimes possible, as is done in expectation propagation (Vehtari et
al., 2020), to match moments between $q$ and an intermediate distribution
$q_{\backslash i}$, which equals $q$ save for its $i^{\text{th}}$ factor,
which is replaced by the $i^{\text{th}}$ factor of $p$. This is equivalent to
minimizing $\text{KL}(q_{\backslash i}||q)$. It is not true, however, that we
minimize $\text{KL}(p||q)$ by iterating through all the factors in this way.

It is illuminating to contrast the properties of the reverse and forward KL
divergences for FG-VI. For $p$ and $q$ in eqs. (3–4), the forward KL
divergence can be obtained from eq. (25) simply by swapping the parameters of
these distributions. In this way we see that

$\textrm{KL}(p||q)=\frac{1}{2}\left[(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu})^{\top}\boldsymbol{\Psi}^{-1}(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu})+\log\frac{|\boldsymbol{\Psi}|}{|\boldsymbol{\Sigma}|}+\textrm{Tr}(\boldsymbol{\Psi}^{-1}\boldsymbol{\Sigma})-n\right].$
(27)

This divergence is similarly minimized by setting
$\boldsymbol{\nu}=\boldsymbol{\mu}$, but now we obtain a different solution
when the remaining terms are minimized with respect to $\boldsymbol{\Psi}$. In
particular, the best factorized approximation minimizing eq. (27) is obtained
when

$\textrm{diag}(\boldsymbol{\Psi})=\textrm{diag}(\boldsymbol{\Sigma});$ (28)

namely, the forward KL divergence is minimized by matching the marginal
distributions of $p$ and $q$; a similar observation is made in MacKay (2003).
By combining this observation with the Impossibility Theorem, we obtain a
counterpart to our previous result:

######  Proposition 3.2 (Variance matching).

If $q$ in eq. (4) has a diagonal covariance but $p$ in eq. (3) does not, then
the forward KL divergence in eq. (27) is minimized when $q$ matches the mean
and marginal variances of $p$ but strictly overestimates the joint entropy and
underestimates all of the marginal precisions (at least one of them strictly).

###  3.2 Rényi divergence

Next we consider the Rényi divergence of order $\alpha$ (or the
$\alpha$-divergence) in Table 1,

$\text{R}_{\alpha}(p||q)=\frac{1}{\alpha(\alpha-1)}\int\
\left(\frac{p^{\alpha}({\bf z})}{q^{\alpha}({\bf z})}-1\right)\,q({\bf
z})\,\text{d}{\bf z},$ (29)

The $\alpha$-divergence has also been studied in the context of VI (Li and
Turner, 2016). It can be defined in various ways, but these various
definitions have in common that they all yield the same underlying
optimization for VI. We use the above definition (Cichocki and Amari, 2010) to
recover the KL divergences in the previous sections as limiting cases:

$\displaystyle\lim_{\alpha\to 0}\text{R}_{\alpha}(p||q)$
$\displaystyle=\text{KL}(q||p),$ (30) $\displaystyle\lim_{\alpha\to
1}\text{R}_{\alpha}(p||q)$ $\displaystyle=\text{KL}(p||q).$ (31)

Thus for $\alpha\\!\in\\!(0,1)$, the Rényi divergences provide a one-parameter
family of divergences interpolating between $\text{KL}(q||p)$ and
$\text{KL}(p||q)$. Likewise, for $\alpha\\!=\\!2$, eq. (29) recovers the
$\chi^{2}$-divergence, which can also be minimized for approximate inference
(Dieng et al., 2017).

Next we consider FG-VI based on minimizing the $\alpha$-divergence in eq.
(29). We restrict our attention to the case $\alpha\in(0,1)$, and as
shorthand, we define the matrices

$\displaystyle\boldsymbol{\Sigma}_{\alpha}$
$\displaystyle=\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma},$
(32) $\displaystyle\boldsymbol{\Phi}_{\alpha}$
$\displaystyle=\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1},$
(33)

which, for $\alpha\in(0,1)$, are convex combinations of the covariance and
precision matrices of $p$ and $q$ in eqs. (3-4). Note that these matrices are
also related by the identity
$\boldsymbol{\Sigma}_{\alpha}=\boldsymbol{\Sigma}\,\boldsymbol{\Phi}_{\alpha}\boldsymbol{\Psi}$.
In terms of these matrices, the $\alpha$-divergence in eq. (29) is given by

$\displaystyle\text{R}_{\alpha}(p||q)$
$\displaystyle=\frac{1}{\alpha(\alpha\\!-\\!1)}\left[e^{-\frac{\alpha(1-\alpha)}{2}(\boldsymbol{\mu}-\boldsymbol{\nu})^{\top}\boldsymbol{\Sigma}_{\alpha}^{-1}(\boldsymbol{\mu}-\boldsymbol{\nu})}\,|\boldsymbol{\Sigma}_{\alpha}|^{-\frac{1}{2}}|\boldsymbol{\Psi}|^{\frac{\alpha}{2}}|\boldsymbol{\Sigma}|^{\frac{1-\alpha}{2}}-1\right],$
(34)

and it is only defined for values of $\alpha$ such that
$\boldsymbol{\Phi}_{\alpha}\\!\succ\\!0$ (Burbea, 1984; Liese and Vajda, 1987;
Hobza et al., 2009; Gil et al., 2013). (If $\boldsymbol{\Phi}_{\alpha}$ is not
positive definite, then the integral in eq. (29) does not exist, and the
$\alpha$-divergence is undefined; when $\alpha\in(0,1)$, it is always the case
that $\boldsymbol{\Phi}_{\alpha}\succ 0$ by convexity of the positive-definite
matrix cone.) With this notation, we have the following result.

######  Proposition 3.3 (Rényi fixed-point equations).

Let $\alpha\\!\in\\!(0,1)$, and let $p$ and $q$ be given by eqs. (3–4) where
$\boldsymbol{\Psi}$ is diagonal. Then the Rényi divergence in eq. (29) is
minimized when $\boldsymbol{\nu}=\boldsymbol{\mu}$ and the estimated
precisions from $\boldsymbol{\Psi}$ satisfy the fixed-point equation ${\rm
diag}\left(\boldsymbol{\Psi}^{-1}\right)={\rm
diag}\left(\boldsymbol{\Sigma}_{\alpha}^{-1}\right),$ (35) or equivalently
when the estimated variances from $\boldsymbol{\Psi}$ satisfy the fixed-point
equation ${\rm diag}(\boldsymbol{\Psi})={\rm
diag}\left(\boldsymbol{\Phi}_{\alpha}^{-1}\right).$ (36)

###### Proof.

When $\alpha\\!\in\\!(0,1)$, the prefactor $1/\alpha(\alpha-1)$ in eq. (34) is
negative, and the expression as a whole is minimized by maximizing the first
term in brackets. Equivalently, upon taking logarithms, we see that

$\operatorname*{arg\hskip
0.43057ptmin}_{q\in\mathcal{Q}}\Big{[}R_{\alpha}(p||q)\Big{]}=\operatorname*{arg\hskip
0.43057ptmin}_{\boldsymbol{\nu},\boldsymbol{\Psi}}\left[\alpha(1\\!-\\!\alpha)(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu})^{\top}\boldsymbol{\Sigma}_{\alpha}^{-1}(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu})+\log|\boldsymbol{\Sigma}_{\alpha}|-\alpha\log|\boldsymbol{\Psi}|\right].$
(37)

The first term on the right is minimized (and zeroed) by setting
$\boldsymbol{\nu}\\!=\\!\boldsymbol{\mu}$, thus matching the means of $q$ and
$p$. Eliminating this term, and substituting for
$\boldsymbol{\Sigma}_{\alpha}$ from eq. (32), we see that

$\operatorname*{arg\hskip
0.43057ptmin}_{q\in\mathcal{Q}}\Big{[}R_{\alpha}(p||q)\Big{]}=\operatorname*{arg\hskip
0.43057ptmin}_{\boldsymbol{\Psi}}\Big{[}\log|\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}|-\alpha\log|\boldsymbol{\Psi}|\Big{]}.$
(38)

The expression on the right of eq. (38) diverges if any diagonal element
$\Psi_{ii}$ vanishes or goes to infinity, so the minimum must exist at some
$\boldsymbol{\Psi}\\!\succ\\!\mathbf{0}$. Solving for the minimum, we find

$\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial\Psi_{ii}}[\log|\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}|-\alpha\log|\boldsymbol{\Psi}|\Big{]},$
(39)
$\displaystyle=\alpha[\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}]^{-1}_{ii}-\alpha\Psi^{-1}_{ii},$
(40)
$\displaystyle=\alpha\left[\boldsymbol{\Sigma}_{\alpha}^{-1}-\boldsymbol{\Psi}^{-1}\right]_{ii},$
(41)

thus proving eq. 35. Finally, we observe that as $\alpha\\!\rightarrow\\!0$,
we recover the precision-matching solution obtained when minimizing
$\text{KL}(q||p)$.

To derive the fixed-point equation in eq. (36), we rewrite eq. (38) as

$\displaystyle\operatorname*{arg\hskip
0.43057ptmin}_{q\in\mathcal{Q}}\Big{[}R_{\alpha}(p||q)\Big{]}$
$\displaystyle=\operatorname*{arg\hskip
0.43057ptmin}_{\boldsymbol{\Psi}}\Big{[}\log|\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}|-\alpha\log|\boldsymbol{\Psi}|+\log|\boldsymbol{\Psi}|-\log|\boldsymbol{\Psi}|\Big{]},$
(42) $\displaystyle=\operatorname*{arg\hskip
0.43057ptmin}_{\boldsymbol{\Psi}}\Big{[}\log\left|\alpha+(1\\!-\\!\alpha)\boldsymbol{\Sigma}\boldsymbol{\Psi}^{-1}\right|-(1\\!-\\!\alpha)\log\left|\boldsymbol{\Psi}^{-1}\right|\Big{]},$
(43) $\displaystyle=\operatorname*{arg\hskip
0.43057ptmin}_{\boldsymbol{\Psi}^{-1}}\Big{[}\log\left|\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}\right|-(1\\!-\\!\alpha)\log\left|\boldsymbol{\Psi}^{-1}\right|\Big{]},$
(44)

where in the last step we have exploited the fact that additive constants,
such as $\log\left|\boldsymbol{\Sigma}^{-1}\right|$, do not change the
location of the minimum with respect to $\boldsymbol{\Psi}$ (or equivalently,
with respect to $\boldsymbol{\Psi}^{-1}$). Solving for the minimum, we find
that

$0=\frac{\partial}{\partial\Psi^{-1}_{ii}}\Big{[}\log\left|\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}\right|-(1\\!-\\!\alpha)\log\left|\boldsymbol{\Psi}^{-1}\right|\Big{]}=(1\\!-\\!\alpha)\left[\boldsymbol{\Phi}_{\alpha}^{-1}-\boldsymbol{\Psi}\right]_{ii},$
(45)

which is equivalent to the claim that
$\textrm{diag}(\boldsymbol{\Psi})=\textrm{diag}\left(\boldsymbol{\Phi}^{-1}_{\alpha}\right)$
in eq. (36). Finally, we observe that, as $\alpha\\!\rightarrow\\!1$, this
solution recovers the variance-matching solution obtained by minimizing
$\text{KL}(p||q)$.

∎

Proposition 3.3 provides two (equivalent) fixed-point equations for the
variational parameter $\boldsymbol{\Psi}$, and it will prove useful to work
with both in Section 4. In componentwise notation, these fixed-point equations
take the form

$\displaystyle\Psi^{-1}_{ii}$
$\displaystyle=\big{[}\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}\big{]}^{-1}_{ii},$
(46) $\displaystyle\Psi_{ii}$
$\displaystyle=\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}\right]^{-1}_{ii},$
(47)

and again we see that eq. (46) reduces to the precision-matching solution
$\Psi^{-1}_{ii}=\Sigma^{-1}_{ii}$ as $\alpha\rightarrow 0$ and eq. (47)
reduces to the variance-matching solution $\Psi_{ii}=\Sigma_{ii}$ as
$\alpha\rightarrow 1$.

###  3.3 Score-based divergence

Next we consider the score-based divergences presented in Table 1. These
divergences are similar to the relative Fisher information (Courtade, 2016,
eq. (7)) except that they use a weighted norm (Cai et al., 2024, Appendix A)
to measure the difference between the gradients of $\log p$ and $\log q$. When
$p$ and $q$ are Gaussian, these score-based divergences are given by

$\displaystyle S(q||p)$ $\displaystyle=\int\text{d}{\bf z}\,q({\bf
z})\,\big{\|}\nabla\log q({\bf z})-\nabla\log p({\bf
z})\big{\|}^{2}_{\text{Cov}(q)},$ (48) $\displaystyle S(p||q)$
$\displaystyle=\int\text{d}{\bf z}\,q({\bf z})\,\big{\|}\nabla\log q({\bf
z})-\nabla\log p({\bf z})\big{\|}^{2}_{\text{Cov}(p)},$ (49)

where as shorthand we use
$\|\mathbf{v}\|^{2}_{\mathbf{A}}=\mathbf{v}^{\top}\mathbf{A}\mathbf{v}$ to
denote the weighted norm for any $\mathbf{A}\\!\succ\\!\mathbf{0}$. The score-
based divergences in eqs. (48-49) are dimensionless measures of the
discrepancy between $p$ and $q$, and in particular, like the KL and Rényi
divergences in the previous sections, they are invariant to affine
reparameterizations of the support of these distributions (Cai et al., 2024,
Theorem A.4). For $p$ and $q$ in eqs. (3–4), these divergences are given by

$\displaystyle S(q||p)$
$\displaystyle=\text{tr}\left[\left(\mathbf{I}\\!-\\!\boldsymbol{\Psi}\boldsymbol{\Sigma}^{-1}\right)^{2}\right]+(\boldsymbol{\nu}\\!-\\!\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{\Psi}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\nu}\\!-\\!\boldsymbol{\mu}),$
(50) $\displaystyle S(p||q)$
$\displaystyle=\text{tr}\left[\left(\mathbf{I}\\!-\\!\boldsymbol{\Sigma}\boldsymbol{\Psi}^{-1}\right)^{2}\right]+(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu})^{\top}\boldsymbol{\Psi}^{-1}\boldsymbol{\Sigma}\boldsymbol{\Psi}^{-1}(\boldsymbol{\mu}\\!-\\!\boldsymbol{\nu}),$
(51)

as also shown in Cai et al. (2024). The rest of this section examines the
solutions that minimize these expressions with respect to the variational mean
$\boldsymbol{\nu}$ and diagonal covariance $\boldsymbol{\Psi}$.

First we consider the minimization of the reverse score-based divergence,
$S(q||p)$. This minimization can be formulated as a nonnegative quadratic
program (NQP)—that is, a convex instance of quadratic programming with
nonnegativity constraints.

######  Proposition 3.4 (NQP for minimizing $S(q||p)$).

Let $p$ and $q$ be given by eqs. (3–4) where $\boldsymbol{\Psi}$ is diagonal.
Then the reverse score-based divergence in eq. (50) is minimized by setting
$\boldsymbol{\nu}\\!=\\!\boldsymbol{\mu}$ and solving the quadratic program
$\min_{{\bf s}\geq 0}\left[\tfrac{1}{2}{\bf s}^{\top}{\bf H}\,{\bf
s}-\mathbf{1\\!}^{\top}\\!{\bf s}\right],$ (52) where ${\bf s}$ is constrained
to lie in the nonnegative orthant, $\mathbf{1}$ is the vector of all ones,
${\bf H}$ is the positive-definite matrix with elements
$H_{ij}=\frac{\left(\Sigma^{-1}_{ij}\right)^{2}}{\Sigma^{-1}_{ii}\Sigma^{-1}_{jj}},$
(53) and $\boldsymbol{\Psi}$ is obtained from the solution of eq. (52) by
identifying $s_{i}=\Psi_{ii}\Sigma^{-1}_{ii}$.

###### Proof.

It is clear that eq. (50) is minimized by setting
$\boldsymbol{\nu}\\!=\\!\boldsymbol{\mu}$ in the rightmost term. For the
remaining term, we note that

$\tfrac{1}{2}\,\text{tr}\left[\left(\mathbf{I}\\!-\\!\boldsymbol{\Psi}\boldsymbol{\Sigma}^{-1}\right)^{2}\right]\
=\
\tfrac{n}{2}\,-\,\sum_{i}\Psi_{ii}\Sigma_{ii}^{-1}\,+\,\tfrac{1}{2}\sum_{ij}\Psi_{ii}\Psi_{jj}\left(\Sigma^{-1}_{ij}\right)^{2},$
(54)

and the NQP in eq. (52) is obtained by defining
$s_{i}=\Psi_{ii}\Sigma_{ii}^{-1}$. Finally, we observe that the matrix
$\mathbf{H}$ in eq. (53) is always positive-definite. Indeed, $\mathbf{H}$ is
the Hadamard (i.e., component-wise) square of the correlation matrix built
from the precision matrix $\mathbf{\Sigma}^{-1}$. Positive-definiteness of
$\mathbf{H}$ follows then from the Schur product theorem that the Hadamard
product of two positive-definite matrices is also positive-definite. ∎

Given the symmetry of eqs. (50-51), it is not surprising that we obtain an
analogous optimization when minimizing the forward score-based divergence,
$S(p||q)$, in terms of the variational parameters $\boldsymbol{\nu}$ and
$\boldsymbol{\Psi}$. In particular, we have the following result.

######  Proposition 3.5 (NQP for minimizing $S(q||p)$).

Let $p$ and $q$ be given by eqs. (3–4) where $\boldsymbol{\Psi}$ is diagonal.
Then the forward score-based divergence in eq. (50) is minimized by setting
$\boldsymbol{\nu}\\!=\\!\boldsymbol{\mu}$ and solving the quadratic program
$\min_{{\bf t}\geq 0}\left[\tfrac{1}{2}{\bf t}^{\top}\mathbf{J}\,{\bf
t}-\mathbf{1\\!}^{\top}\\!{\bf t}\right],$ (55) where ${\bf t}$ is constrained
to lie in the nonnegative orthant, $\mathbf{J}$ is the positive-definite
matrix with elements
$J_{ij}=\frac{\left(\Sigma_{ij}\right)^{2}}{\Sigma_{ii}\Sigma_{jj}},$ (56) and
$\boldsymbol{\Psi}$ is obtained from the solution of eq. (55) by identifying
$t_{i}=\Psi^{-1}_{ii}\Sigma_{ii}$.

###### Proof.

The proof follows the same steps as in the previous proposition, but with the
matrices $\boldsymbol{\Psi}$ and $\boldsymbol{\Sigma}$ playing the reverse
roles as they did in eq. (54). ∎

We conclude this section by noting a peculiar property of the score-based
divergences in Table 1. It is possible for the solutions of the NQPs in
Propositions 3.4 and 3.5 to lie on the boundary of the nonnegative orthant.
That is, there may exist some $i^{\text{th}}$ component along the diagonal of
$\boldsymbol{\Psi}$ such that $S(q||p)$ is minimized by setting
$\Psi_{ii}\\!=\\!0$ (in Proposition 3.4) or such that $S(p||q)$ is minimized
by setting $\Psi_{ii}^{-1}\\!=\\!0$ (in Proposition 3.5). Strictly speaking,
solutions of this form—with zero or infinite marginal variances—do not define
proper distributions that lie in the family $\mathcal{Q}$ of factorized
approximations. Such solutions never arise with FG-VI based on the KL or Rényi
divergences for $\alpha\in\\{0,1\\}$. This phenomenon gives rise to the
following definition.

######  Definition 3.6.

Let $\mathcal{Q}$ denote the family of factorized Gaussian approximations in
eq. (2) that are proper distributions (i.e., with
$0\\!<\\!\mathcal{H}(q)\\!<\\!\infty$). For some target distribution $p$ and
divergence $D$, we say that FG-VI undergoes variational collapse when
$\arg\\!\inf_{q\in\mathcal{Q}}D(q,p)\not\in\mathcal{Q}$.

Fig. 2 illustrates this collapse when FG-VI is used to approximate a
multivariate Gaussian distribution in three dimensions with varying degrees of
correlation; the amount of correlation is determined by the magnitude of the
off-diagonal elements $C_{12}$, $C_{23}$, and $C_{13}$ in eq. (6). These off-
diagonal elements are constrained by the fact that the correlation matrix as a
whole must be positive definite. Each panel in the figure visualizes a two-
dimensional slice of the three-dimensional (convex) set of positive-definite
correlation matrices for a fixed value of $C_{23}$. The blue regions of each
slice indicate the areas where the NQPs in Propositions 3.4 and 3.5 yield
finite and nonzero estimates of the variances; the red regions indicate areas
of variational collapse. (It can be shown, in three dimensions, that the
forward and reverse score-based divergences share the same areas of
variational collapse.) Interestingly, the leftmost panel shows that
variational collapse does not occur in two dimensions or arise purely from
pairwise correlations. On the other hand, the remaining panels show that
variational collapse occurs in three dimensions whenever the correlation
matrix is dense and has at least one off-diagonal element of sufficiently
large magnitude.

Figure 2: When FG-VI is based on minimizing the score-based divergences in
Table 1, it may estimate zero or infinite values for the marginal variances.
The red areas indicate these occurrences of variational collapse when FG-VI
with a score-based divergence is used to approximate a three-dimensional
Gaussian with a non-diagonal correlation matrix $\mathbf{C}$.

##  4 Ordering of divergences for FG-VI

We have seen that different divergences $D(q,p)$ yield different solutions to
the problem of variational inference in eq. (1); we have also seen, in turn,
that these different solutions yield different estimates of the marginal
variances, precisions, and entropy. When it is possible to order these
different estimates of uncertainty, it is also natural to define a
corresponding ordering of the divergences from which they were derived. In
this section, we derive such an ordering over all the divergences listed in
Table 1.

We begin with the following definition, which bases this ordering on the
marginal variances that are estimated by FG-VI.

######  Definition 4.1.

Let $\mathcal{P}$ and $\mathcal{Q}$ denote the families of multivariate
Gaussian distributions in eqs. (3–4), where $\boldsymbol{\Psi}$ is diagonal
but $\boldsymbol{\Sigma}$ is not. Consider two divergences $D_{a}$ and
$D_{b}$, and for any $p\\!\in\\!\mathcal{P}$, let $\boldsymbol{\Psi}(a)$ and
$\boldsymbol{\Psi}(b)$ denote the covariances of the factorized approximations
obtained, respectively, by minimizing $D_{a}(q,p)$ and $D_{b}(q,p)$ over all
$q\\!\in\\!\mathcal{Q}$. We say $D_{a}$ dominates $D_{b}$, written $D_{a}\succ
D_{b}$, if for any $p\in\mathcal{P}$ we have that $\displaystyle\forall
i,\quad\Psi_{ii}(a)$ $\displaystyle\geq\Psi_{ii}(b),$ (57)
$\displaystyle\exists j,\\!\quad\Psi_{jj}(a)$ $\displaystyle>\Psi_{jj}(b).$
(58)

From the above definition, we see intuitively that $D_{a}$ dominates $D_{b}$
in the following sense: it yields no smaller estimates of the marginal
variances—and at least one larger estimate—when $D_{a}$ and $D_{b}$ are used
for FG-VI with a Gaussian target distribution.

The above definition was based on the estimates of marginal variances by FG-
VI, and we shall use it throughout the rest of the paper. But it is not the
only possible definition. Because FG-VI estimates a strictly diagonal
covariance matrix, the ordering in Definition 4.1 also immediately implies an
ordering of the marginal precisions and entropies.

######  Proposition 4.2.

Suppose $D_{a}\succ D_{b}$. Then for any $p\\!\in\\!\mathcal{P}$: (i) the
marginal precisions are ordered, with
$\mathrm{diag}(\boldsymbol{\Psi}^{-1}(a))\leq\mathrm{diag}(\boldsymbol{\Psi}^{-1}(b))$,
and this inequality is strict for at least one component along the diagonal.
(ii) the entropies are ordered, with $\mathcal{H}(q_{a})<\mathcal{H}(q_{b})$,
where $q_{a}$ and $q_{b}$ are the solutions obtained by minimizing
$D_{a}(q,p)$ and $D_{b}(q,p)$ over all $q\in\mathcal{Q}$. These implications
may be useful for applications which require the estimation of precisions
and/or entropies rather than variances.

Figure 3:  Comparison of variances, precisions, and entropy estimated by FG-VI
with different divergences. Here FG-VI was used to approximate a 2-dimensional
Gaussian with a non-diagonal correlation matrix $\mathbf{C}$. The Rényi
divergence was computed for $\alpha\\!=\\!0.5$. The ordering of the curves
matches the predictions of Theorem 4.3 and Proposition 4.2.

With Definition 4.1 in hand, we are now able to state the main result of this
paper; it is an ordering over all the divergences listed in Table 1.

######  Theorem 4.3 (Ordering of divergences).

For any $(\alpha_{1},\alpha_{2})$ satisfying
$0\\!<\\!\alpha_{1}\\!<\\!\alpha_{2}\\!<\\!1$,
$\text{S}(q||p)\prec\text{KL}(q||p)\prec R_{\alpha_{1}}(p||q)\prec
R_{\alpha_{2}}(p||q)\prec\text{KL}(p||q)\prec\text{S}(p||q).$ (59) Figure 3
illustrates the consequences of this ordering when FG-VI is used to
approximate a 2-dimensional Gaussian with varying degrees of correlation. The
amount of correlation is measured by the determinant of the correlation
matrix, $|\mathbf{C}|$; as this determinant decreases from one to zero—from no
correlation to perfect correlation—the target Gaussian more starkly violates
the assumption of factorizability in FG-VI. From left to right, the panels
plot the estimates of variance, precision, and entropy that are obtained from
FG-VI with different divergences, and we see that for all values of
$|\mathbf{C}|$, these estimates are ordered exactly as predicted by Theorem
4.3 and Proposition 4.2. The rest of this section is devoted to proving this
theorem.

In fact, we have already proven one of these orderings. Recall that when FG-VI
is based on minimizing $\text{KL}(q||p)$, it underestimates the marginal
variances, at least one of them strictly; this is the content of Lemma 2.4 and
Proposition 3.1. On the other hand, when FG-VI is based on minimizing
$\text{KL}(p||q)$, it correctly estimates the marginal variances (MacKay,
2003); this is the content of Proposition 3.2. Thus we have already shown
$\text{KL}(q||p)\\!\prec\\!\text{KL}(p||q)$.

###  4.1 Ordering of score-based divergences

In what follows we prove the two outermost orderings in Theorem 4.3;
specifically, we show that $S(q||p)\prec\text{KL}(q||p)$ and
$\text{KL}(p||q)\prec S(p||q)$. At a high level, these proofs are obtained by
analyzing the Karush-Kuhn-Tucker (KKT) conditions of the quadratic programs in
Propositions 3.4 and 3.5.

######  Proof of $S(q||p)\prec\text{KL}(q||p)$.

Consider the diagonal covariance matrix $\boldsymbol{\Psi}$ that FG-VI
estimates by minimizing $S(q||p)$ on the left side of this ordering. This
matrix is obtained from the solution of the NQP in eq. 52 for the unknown
nonnegative variables $s_{i}=\Psi_{ii}\Sigma^{-1}_{ii}$. The solution to this
NQP must satisfy the KKT conditions: namely, for each component of
$\mathbf{s}$, we have either that (i) $s_{i}\\!=\\!0$ and
$\mathbf{Hs}_{i}>\\!1$ or (ii) $s_{i}\\!>\\!0$ and $\mathbf{Hs}_{i}\\!=\\!1$.
We examine each of these cases in turn:

$s_{i}\\!=\\!0$ |  This is a case of variational collapse where by minimizing $S(q||p)$ we estimate that $\Psi_{ii}=0$.
---|---
$s_{i}\\!>\\!0$ |  In this case we observe that $(\mathbf{Hs})_{i}\\!=\\!1$ from the KKT condition and also that $H_{ij}\geq 0$ and $H_{ii}=1$ from eq. (53). Thus we see that $s_{i}=H_{ii}s_{i}=1-\sum_{j\neq i}H_{ij}s_{j}\leq 1,$ (60) or equivalently, that $\Psi_{ii}\Sigma^{-1}_{ii}\leq 1$.

We have thus shown that $\Psi_{ii}\Sigma_{ii}^{-1}\\!\leq\\!1$ for all of the
marginal variances estimated by minimizing the score-based divergence
$S(q||p)$. In fact, we claim that this inequality must be strict for at least
one component. The claim is trivially true if any $s_{i}\\!=\\!0$, since this
implies $\Psi_{ii}\Sigma_{ii}^{-1}\\!=\\!0$. Suppose to the contrary that
every element of $\mathbf{s}$ is strictly positive. Since by assumption
$\boldsymbol{\Sigma}$ is not diagonal, it must also be true that
$\boldsymbol{\Sigma}^{-1}$ is not diagonal, and therefore from eq. (53) we can
pick some $i$ and $k$ such that $H_{ik}\\!>\\!0$. But then, continuing from
eq. (60), we see that

$s_{i}=1-\sum_{j\neq i}H_{ij}s_{j}\leq 1-H_{ik}s_{k}<1,$ (61)

or equivalently that $\Psi_{ii}\Sigma_{ii}^{-1}\\!<\\!1$, thus proving the
claim. In sum we have shown that $\Psi_{ii}\\!\leq\\!1/\Sigma^{-1}_{ii}$ for
all of the marginal variances estimated by minimizing $S(q||p)$, and also that
this inequality is strict for at least one component. Finally, we recall from
eq. (26) that $1/\Sigma^{-1}_{ii}$ is the marginal variance estimated by
minimizing $\text{KL}(q||p)$. It follows from Definition 4.1 that
$S(q||p)\prec\text{KL}(q||p)$. ∎

######  Proof of $\text{KL}(p||q)\prec S(q||p)$.

Consider the diagonal covariance matrix $\boldsymbol{\Psi}$ that FG-VI
estimates by minimizing $S(p||q)$. This matrix is obtained from the solution
of the NQP in eq. 55 for the unknown nonnegative variables
$t_{i}=\Psi^{-1}_{ii}\Sigma_{ii}$. The rest of the proof follows exactly the
same steps as the previous one, but with $t_{i}$ playing the role of $s_{i}$.
In this way we can show that $\Psi^{-1}_{ii}\Sigma_{ii}\\!\leq\\!1$ (or
equivalently $\Sigma_{ii}\\!\leq\\!\Psi_{ii}$) for all of the marginal
variances estimated by minimizing $S(p||q)$, and also that this inequality is
strict for at least one component. Finally, we recall from eq. (28) that
$\Sigma_{ii}$ is the marginal variance estimated by minimizing
$\text{KL}(p||q)$, and it follows from Definition 4.1 that
$\text{KL}(p||q)\prec S(p||q)$. ∎

###  4.2 Ordering of KL divergences

Next we prove the two intermediate orderings in Theorem 4.3; specifically, for
any $\alpha\\!\in\\!(0,1)$, we show that
$\text{KL}(q||p)\prec\text{R}_{\alpha}(p||q)$ and
$\text{R}_{\alpha}(p||q)\prec\text{KL}(p||q)$. At a high level, these proofs
are obtained by applying the inequalities in Proposition 2.6 to the fixed-
point equations for the marginal variances in eqs. (46–47); recall that these
were the marginal variances from FG-VI estimated by minimizing the Rényi
divergence of order $\alpha$.

######  Proof of $\text{KL}(q||p)\prec\text{R}_{\alpha}(p||q)$ for
$\alpha\\!\in\\!(0,1)$.

Let $\boldsymbol{\Psi}$ denote the diagonal covariance matrix estimated by
minimizing $\text{R}_{\alpha}(p||q)$, and consider the matrix
$[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}]^{-1}$ that
appears on the right side of its fixed-point equation in eq. (47). This matrix
is positive-definite and non-diagonal for all $\alpha\\!\in\\!(0,1)$, and thus
we can apply the inequality in eq. (23) of Proposition 2.6 to its diagonal
elements. In this way, we find:

$\displaystyle\alpha\Psi_{ii}\Sigma^{-1}_{ii}$
$\displaystyle=\Psi_{ii}\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}\right]_{ii}\,-\,(1\\!-\\!\alpha),$
(62)
$\displaystyle\geq\frac{\Psi_{ii}}{\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}\right]^{-1}_{ii}}\,-\,(1\\!-\\!\alpha),$
(63) $\displaystyle=(\Psi_{ii}/\Psi_{ii})\,-\,(1\\!-\\!\alpha),$ (64)
$\displaystyle=\alpha.$ (65)

Dividing both sides by $\alpha\Sigma^{-1}_{ii}$, we see that $\Psi_{ii}\geq
1/\Sigma_{ii}^{-1}$ for all $i$. Next we appeal to the strict inequality of
Proposition 2.6 and conclude that $\Psi_{jj}>1/\Sigma_{jj}^{-1}$ for some $j$.
On the left and right sides of these inequalities appear, respectively, the
marginal variances estimated by minimizing $R_{\alpha}(p,q)$ and
$\text{KL}(q,p)$. Thus we have shown
$\text{R}_{\alpha}(p||q)\succ\text{KL}(q,p)$. ∎

######  Proof of $\text{R}_{\alpha}(p||q)\prec\text{KL}(p||q)$ for
$\alpha\\!\in\\!(0,1)$.

The proof follows the same steps as the previous one, but now we combine the
inequalities of Proposition 2.6 with the fixed-point equation in eq. (46). In
this way we find:

$(1\\!-\\!\alpha)\Psi^{-1}_{ii}\Sigma_{ii}\,=\,\Psi^{-1}_{ii}\left[\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}\right]_{ii}\\!-\\!\alpha\,\geq\,\frac{\Psi^{-1}_{ii}}{\left[\alpha\boldsymbol{\Psi}+(1\\!-\\!\alpha)\boldsymbol{\Sigma}\right]^{-1}_{ii}}\\!-\\!\alpha\,=\,1\\!-\\!\alpha.$
(66)

Dividing both sides by $(1\\!-\\!\alpha)\Psi^{-1}_{ii}$, we see that
$\Sigma_{ii}\geq\Psi_{ii}$ for all $i$, and we know from the strict inequality
in eq. (24) of Proposition 2.6 that $\Sigma_{jj}>\Psi_{jj}$ for some $j$. On
the left and right sides of these inequalities appear, respectively, the
marginal variances estimated by minimizing $\text{KL}(p||q)$ and
$\text{R}_{\alpha}(p||q)$. Thus we have shown
$\text{KL}(p||q)\succ\text{R}_{\alpha}(p||q)$. ∎

###  4.3 Ordering of Rényi divergences

Finally we prove the ordering of Rényi divergences in Theorem 4.3. This proof
is more technical, relying on a detailed analysis of the fixed-point equation
in eq. (47).

######  Proof of $R_{\alpha_{1}}(p||q)\prec R_{\alpha_{2}}(p||q)$ for
$0\\!<\alpha_{1}\\!<\\!\alpha_{2}\\!<\\!1$.

Let $\boldsymbol{\Psi}(\alpha)$ denote the diagonal covariance matrix that FG-
VI estimates by minimizing $\text{R}_{\alpha}(p||q)$; here, we regard this
matrix as a function over the interval $(0,1)$. The function is continuous,
and we shall prove the ordering by showing that the estimated variance
$\Psi_{ii}(\alpha)$ is strictly increasing whenever
$\Sigma_{ii}\Sigma_{ii}^{-1}\\!>\\!1$.

Let us review what we have already shown. For $\alpha\\!\in\\!(0,1)$, we know
from the orderings of Section 4.2 that $\Psi_{ii}(\alpha)$ is sandwiched
between its limiting values at zero and one:

$\left(\Sigma_{ii}^{-1}\right)^{-1}\,=\,\lim_{\alpha\rightarrow
0^{+}}\Psi_{ii}(\alpha)\,\leq\,\Psi_{ii}(\alpha)\,\leq\,\lim_{\alpha\rightarrow
1^{-}}\Psi_{ii}(\alpha)\,=\,\Sigma_{ii}.$ (67)

We also know that $\Sigma_{ii}\Sigma^{-1}_{ii}>1$ whenever $\sum_{j\neq
i}|\Sigma_{ij}|>0$, and that in this case, the above inequalities are strict.
A useful picture of this situation is shown in the left panel of Fig. 4.
Consider any component $\Psi_{ii}(\alpha)$ of the estimated variances that is
not constant on the unit interval. (As shown previously, there must be at
least one such component.) In this case, there are only two possibilities:
either $\Psi_{ii}(\alpha)$ is strictly increasing, or it is not. If the former
is always true, then it follows that $R_{\alpha_{1}}(p||q)\prec
R_{\alpha_{2}}(p||q)$ whenever $0\\!<\\!\alpha_{1}\\!<\\!\alpha_{2}\\!<\\!1$.
We shall prove by contradiction that this is indeed the case.

Assume that there are one or more diagonal components, $\Psi_{ii}(\alpha)$,
that are neither constant nor strictly increasing over the unit interval;
also, let $\mathcal{I}$ denote the set that contains the indices of these
components. Since each such component must have at least one stationary point,
we can also write:

$\mathcal{I}=\Big{\\{}i\,\Big{|}\,\Sigma_{ii}\Sigma_{ii}^{-1}\\!>\\!1\
\mbox{and}\ \Psi_{ii}^{\prime}(\alpha)=0\ \mbox{for some}\
\alpha\in(0,1)\Big{\\}}.$ (68)

For each such component, there must also exist some minimal point
$\tau_{i}\in(0,1)$ where its derivative vanishes. We define
$\tau=\min_{i\in\mathcal{I}}\tau_{i}$ and
$j=\arg\\!\min_{i\in\mathcal{I}}\tau_{i}$, so that by definition

$\displaystyle\Psi_{ii}^{\prime}(\tau)$ $\displaystyle\geq 0\quad\mbox{for all
$i$},$ (69) $\displaystyle\Psi_{jj}^{\prime}(\tau)$ $\displaystyle=0.$ (70)

To obtain the desired contradiction, we shall prove that there exists no point
$\tau$ (as imagined in Fig. 4) with these properties.

$\alpha$$\Psi_{jj}(\alpha)$$1/\Sigma^{-1}_{jj}$$\Sigma_{jj}$10${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\tau}$

$\alpha$$f(\alpha)$${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}f(0)=f(\tau)=\Psi_{jj}(\tau)}$0${\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}\tau}$1

Figure 4:  (Left) Either $\Psi_{jj}(\alpha)$ is strictly increasing over
$\alpha\\!\in\\!(0,1)$, or it is not, with some minimal point $\tau$ of
vanishing derivative. We prove the former by showing that no such point $\tau$
exists. (Right) The proof is based on properties of the function $f(\alpha)$
in eq. 75. The function is convex; it also satisfies
$f(0)\\!=\\!f(\tau)\\!=\\!\Psi_{jj}(\tau)$ and
$f^{\prime}(0)f^{\prime}(\tau)<0$.

We start by rewriting the fixed-point equation in LABEL:prop:renyi-variance-
matching and eq. (47) in a slightly different form,

$\Psi_{ii}(\alpha)=[\boldsymbol{\Phi}(\alpha)^{-1}]_{ii}=\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\Psi^{-1}(\alpha)\right]_{ii}^{-1},$
(71)

where on the right side we have explicitly indicated all sources of dependence
on $\alpha$. Next, we differentiate both sides of eq. (71) for $i\\!=\\!j$ and
at $\alpha\\!=\\!\tau$:

$\Psi_{jj}^{\prime}(\tau)=\mathbf{e}_{j}^{\top}\\!\left\\{\frac{d}{d\alpha}\Big{[}\tau\boldsymbol{\Sigma}^{-1}\\!+(1\\!-\\!\tau)\boldsymbol{\Psi}^{-1}(\alpha)\Big{]}^{-1}\\!\\!+\frac{d}{d\alpha}\Big{[}\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}(\tau)\Big{]}^{-1}\right\\}\bigg{|}_{\alpha=\tau}\\!\\!\\!\mathbf{e}_{j},$
(72)

where on the right side we have separated out the different sources of
dependence on $\alpha$. Evaluating the first term on the right side, we find
that

$\mathbf{e}_{j}^{\top}\\!\left\\{\frac{d}{d\alpha}\Big{[}\tau\boldsymbol{\Sigma}^{-1}\\!+(1\\!-\\!\tau)\boldsymbol{\Psi}^{-1}(\alpha)\Big{]}^{-1}\right\\}\bigg{|}_{\alpha=\tau}\\!\\!\\!\\!\mathbf{e}_{j}=\mathbf{e}_{j}^{\top}\boldsymbol{\Phi}^{-1}(\tau)\boldsymbol{\Psi}^{-1}(\tau)\boldsymbol{\Psi}^{\prime}(\tau)\boldsymbol{\Psi}^{-1}(\tau)\boldsymbol{\Phi}^{-1}(\tau)\,\mathbf{e}_{j},$
(73)

and now we recognize that this term is nonnegative: of the matrices on the
right side of eq. (73), note that $\boldsymbol{\Phi}(\alpha)$ and
$\boldsymbol{\Psi}(\alpha)$ are positive definite for all
$\alpha\\!\in\\!(0,1)$, and $\boldsymbol{\Psi}^{\prime}(\tau)$ is positive
semidefinite by virtue of the defining properties of $\tau$ in eqs. (69–70).
Dropping this first term from the right side of eq. (72), we obtain the
inequality

$0\geq\frac{d}{d\alpha}\left\\{\mathbf{e}_{j}^{\top}\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}(\tau)\right]^{-1}\mathbf{e}_{j}\right\\}\Big{|}_{\alpha=\tau},$
(74)

where we have exploited that the left side of eq. (72) vanishes by the second
defining property of $\tau$ in eq. (70).

Next we show that eq. (74) contains a contradiction. To do so, we examine the
properties of the function defined on $[0,1]$ by

$f(\alpha)=\mathbf{e}_{j}^{\top}\left[\alpha\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\alpha)\boldsymbol{\Psi}^{-1}(\tau)\right]^{-1}\mathbf{e}_{j}.$
(75)

This is, of course, exactly the function that is differentiated in eq. (74),
whose right side is equal to $f^{\prime}(\tau)$. This function possesses
several key properties. First, we note that

$f(0)=\mathbf{e}_{j}^{\top}\left[\boldsymbol{\Psi}^{-1}(\tau)\right]^{-1}\mathbf{e}_{j}=\Psi_{jj}(\tau)=\mathbf{e}_{j}^{\top}\left[\tau\boldsymbol{\Sigma}^{-1}+(1\\!-\\!\tau)\boldsymbol{\Psi}^{-1}(\tau)\right]^{-1}\mathbf{e}_{j}=f(\tau).$
(76)

Second, we note that $f$ is convex on the unit interval; here, we are
observing the convexity of the matrix inverse (Nordström, 2011). Third, we
note that

$f^{\prime}(0)=-\mathbf{e}_{j}^{\top}\boldsymbol{\Psi}(\tau)\left[\boldsymbol{\Sigma}^{-1}\\!-\\!\boldsymbol{\Psi}^{-1}(\tau)\right]\boldsymbol{\Psi}(\tau)\,\mathbf{e}_{j}=\Psi_{jj}(\tau)\left[\boldsymbol{\Sigma}^{-1}\\!-\\!\boldsymbol{\Psi}^{-1}(\tau)\right]_{jj}<0.$
(77)

Finally, we note that the sign of $f^{\prime}(\tau)$ must oppose the sign of
$f^{\prime}(0)$; as shown in the right panel of Fig. (4), this follows from
the preceding properties that $f(0)\\!=\\!f(\tau)$, that the function $f$ is
convex on $[0,1]$, and that $f^{\prime}(0)\\!<\\!0$. Thus we have shown that
$f^{\prime}(\tau)\\!>\\!0$. But this is directly in contradiction with eq.
(74), which states that $0\geq f^{\prime}(\tau)$, and we are forced to
conclude that no such point $\tau$ exists. This completes the proof. ∎

Figure 5: When $p$ is Gaussian over $\mathbb{R}^{n}$, there always exists a
unique $\alpha\in(0,1)$ such that the factorized approximation $q$ minimizing
$\text{R}_{\alpha}(p||q)$ matches the entropy of $p$. However, the entropy-
matching value of $\alpha$ depends on the dimension and covariance structure
of $p$. Above we vary the dimension $n$ and constant correlation
$\varepsilon$.

From the above, we can also prove the existence of an entropy-matching
solution for FG-VI analogous to the precision-matching and variance-matching
solutions of Section 3.1.

######  Corollary 4.4.

Let $q_{\alpha}=\arg\\!\min_{q\in\mathcal{Q}}\text{R}_{\alpha}(p||q)$. There
exists a unique value $\alpha\in(0,1)$ such that
$\mathcal{H}(q_{\alpha})=\mathcal{H}(p)$.

###### Proof.

Let $\boldsymbol{\Psi}(\alpha)$ denote the covariance of $q_{\alpha}$ for
$\alpha\\!\in\\!(0,1)$. Since $\Psi_{ii}(\alpha)$ is continuous with respect
to $\alpha$, so is
$\log|\boldsymbol{\Psi}(\alpha)|=\sum_{i}\log\Psi_{ii}(\alpha)$, and hence so
is the entropy $\mathcal{H}(q_{\alpha})$. Moreover, from the ordering of Rényi
divergences, we see that $\log|\Psi(\alpha)|$ and $\mathcal{H}(q_{\alpha})$
are not only continuous, but strictly increasing over the unit interval
$\alpha\\!\in\\!(0,1)$. In addition, we have

$1/\Sigma^{-1}_{ii}=\Psi_{ii}(0)\leq\Psi_{ii}(\alpha)\leq\Psi_{ii}(1)=\Sigma_{ii},\
\forall i.$ (78)

with $q_{\alpha}$ interpolating smoothly between the precision-matching and
the variance-matching approximations of FG-VI. From the impossibility result
in Theorem 2.1, the precision-matching approximation (at $\alpha\\!=\\!0$)
underestimates the entropy of $p$, while the variance-matching approximation
(at $\alpha\\!=\\!1$) overestimates it. From continuity and strict
monotonicity, it follows that there exists a unique $\alpha\\!\in\\!(0,1)$
such that $\mathcal{H}(q_{\alpha})=\mathcal{H}(p)$. ∎

The above proves the existence an entropy-matching solution, but it is not a
constructive proof. There does not exist a single value of $\alpha$ for which
$\mathcal{H}(q_{\alpha})=\mathcal{H}(p)$, even when $p$ is Gaussian. We show
this numerically for the case where the target covariance
$\boldsymbol{\Sigma}$ has unit diagonal elements and constant off-diagonal
elements $\varepsilon\\!\in\\!(0,1)$. As seen in Figure 5, the entropy-
matching value of $\alpha$ (obtained by a grid search) changes as we vary the
amount of correlation and the dimension of the problem.

##  5 Numerical experiments

Does the ordering of divergences hold when FG-VI is applied to non-Gaussian
targets? We study this question empirically on a range of models; in most of
these models, the target $p$ is only known up to a normalizing constant.

When $p$ is not Gaussian, some divergences in Table 1 are harder to minimize
numerically than others. We use the following approaches. We minimize
$\text{KL}(q||p)$ by (equivalently) maximizing the evidence lower bound
(ELBO), and we estimate the ELBO via Monte Carlo with draws from $q$, a
procedure at the core of modern implementations of “black-box” VI (Kucukelbir
et al., 2017). We use a similar procedure to optimize the $\alpha$-divergence;
however, this procedure is more fraught111For more discussion on implementing
VI with $\alpha$-divergences, see e.g., Hernandez-Lobato et al. (2016); Li and
Turner (2016); Dieng et al. (2017); Daudel et al. (2021, 2023). as Monte Carlo
estimators of the $\alpha$-divergence and its gradients suffer from large
variance, especially in large dimensions and as $\alpha$ approaches 1 (Geffner
and Domke, 2021). We minimize $\text{S}(q||p)$ by adapting a recent batch-and-
match (BaM) method inspired by proximal point algorithms (Cai et al., 2024);
BaM was originally developed for Gaussian variational families with dense
covariance matrices, but here we experiment with BaM updates for FG-VI
(derived in Appendix A.1). We minimize $\text{KL}(p||q)$ by estimating the
means and marginal variances of $p$ from long runs of Markov chain Monte
Carlo; the factorized approximation minimizing $\text{KL}(q||p)$ is the one
that matches these moments (MacKay, 2003). We do not empirically evaluate the
forward score-based divergence for FG-VI as we do not have a reliable method
to minimize $\text{S}(p||q)$.

Model | $n$ | Details | Variances | Entropy
---|---|---|---|---
Gaussian | 10 | Full-rank covariance matrix. | ✓ | ✓
Rosenbrock | 2 | High curvature. | ✓ | ✓
Eight Schools | 10 | Small hierarchical linear model. | ✗ | ✓
German Credit | 25 | Logistic regression. | ✓ | ✓
Radon Effect | 91 | Large hierarchical linear model. | ✗ | ✓
Stochastic Volatility | 103 | Stochastic volatility model. | ✗ | ✓
Table 2:  Target distributions $p$ with dimension $n$. The two right-most
columns indicate whether the divergences are ordered according to the marginal
variances and the entropy.

We experiment with the above methods on a diverse collection of target
distributions from the inference gym, a curated library of models for the
study of inference algorithms (Sountsov et al., 2020). We describe these
target distributions briefly in Table 2 and provide more details in Appendix
A.2. We implement each model in Python and use JAX to calculate derivatives
(Bradbury et al., 2018). We optimize the $\alpha$-divergences for
$\alpha\\!=\\!0.1$ and $\alpha\\!=\\!0.5$, but only on low-dimensional targets
($D\\!<\\!10$) where the optimizers are numerically stable. We run each VI
algorithm for 1,000 iterations, and at each iteration, we use $B=10,000$ draws
from $q$ to estimate the divergence and its gradient. This large number of
draws is overkill for many problems, but it helps in practice to stabilize the
optimizations of $\alpha$-divergences. In addition, we parallelize many
calculations on a GPU, keeping run times short. We provide code to run these
experiments at https://github.com/charlesm93/VI-ordering.

Figure 6:  Variance when targeting the posterior distribution of the Eight
Schools model.

The ordering of divergences that we proved for Gaussian targets is not
guaranteed to hold for non-Gaussian targets, and empirically we observe that
it is more likely to fail in high-dimensional problems; see the penultimate
column of Table 2. Intuitively, this is to be expected because the ordering in
Theorem 4.3 is voided by a violation in any component of the marginal
variance. This is illustrated in Fig. 6 for the Eight Schools models, where
some marginal variances are ordered in the same way as predicted for Gaussian
targets, but most are not. The true marginal variances in this model were
estimated by long runs of Markov chain Monte Carlo.

When $p$ is not Gaussian, an alternative ordering of divergences for FG-VI is
suggested by the estimates they yield of the joint entropy. As shown in the
final column of Table 2 and in Fig. 7, the ordering in Theorem 4.3 does
correspond—across all of the non-Gaussian targets in our experiments—to the
ordering of entropies estimated by FG-VI. In these experiments, however, we
are not able to compare the entropies estimated by FG-VI to the actual
entropies. To compute the entropy of non-Gaussian distributions, it is
necessary to estimate their normalizing constants; this can be done by
invoking a Gaussian-like approximation, as in bridge sampling (Meng and
Schilling, 2002; Gronau et al., 2020), however we do not wish to compare a
theory for Gaussian targets to an empirical benchmark which relies on a
Gaussian approximation.

Figure 7:  Ordering of entropy across targets. The zero on each vertical axis
corresponds to the average of the entropies estimated by the approximations in
each panel.

##  6 Discussion

Divergences over probability spaces are ubiquitous in the statistics and
machine learning literature, but it remains challenging to reason about them.
This challenge is salient in VI, where the divergence does not serve merely as
an object for theoretical study but also bears directly on the algorithmic
implementation of the method. Our work reveals how the choice of divergence
can align (or misalign) with intelligible inferential goals, and it is in line
with a rich and recent body of work on assessing the quality of VI (e.g, Yao
et al., 2018; Wang and Blei, 2018; Huggins et al., 2020; Dhaka et al., 2021;
Biswas and Mackey, 2024). A distinguishing feature of our paper is that it
analyzes multiple and competing inferential goals—goals which cannot be
simultaneously met and which are each best served by a different divergence.
Also, our work focuses on a setting that is especially conducive to analysis,
one in which the inherent trade-offs of factorized VI can be completely
characterized, as opposed to more general treatments of the problem. The
intuitions developed from this analysis suggest useful directions,
conjectures, and insights for studying broader classes of distributions.

There exists previous work comparing divergences, notably in the context of
latent variable models, where VI is used to approximate a marginal likelihood.
In this context, the KL and Rényi divergences have been ordered based on the
(lower or upper) bound they provide on the marginal likelihood (Li and Turner,
2016; Dieng et al., 2017). (Interestingly, this ordering agrees with the one
provided in this paper, even though it is based on a different criterion.)
This type of ordering cannot be extended to score-based divergences, however,
as they do not provide bounds on the marginal likelihood.

We have seen that the results of VI can be improved, for certain inferential
goals, by choosing a particular divergence. But this is not the only way to
improve the results of VI; another is to use a richer family $\mathcal{Q}$ of
variational approximations (e.g, Hoffman and Blei, 2015; Johnson et al., 2016;
Dhaka et al., 2021). The deficits of FG-VI may be addressed, for instance, by
Gaussian variational approximations with full covariance matrices (Opper and
Archambeau, 2009; Cai et al., 2024). In high dimensions, Dhaka et al. (2021)
argue that the results of VI can be improved more easily by optimizing
$\text{KL}(q||p)$ over a richer variational family than by optimizing certain
alternative divergences. Nevertheless, even as the variational family becomes
richer, it is not likely to contain the target distribution $p$, and therefore
the approximating distribution $q$ must in some way be compromised. For
example, if $\mathcal{Q}$ is a family of symmetric distributions and the
target $p$ is asymmetric, then an approximation may only be able to match the
mode or the mean of the target. It stands to reason that further impossibility
theorems, in the spirit of Theorem 2.1, can be derived to elucidate these
compromises.

## Aknowledgments

We thank David Blei, Diana Cai, Robert Gower, Chirag Modi, and Yuling Yao for
helpful discussions. We are also grateful to anonymous reviewers from the
conference on Uncertainty in Artificial Intelligence for their comments on a
previous paper (Margossian and Saul, 2023), which provided some of the
motivation for this manuscript.

## References

  * Bernardo and Smith (2000) Bernardo, J. M. and A. F. M. Smith (2000).  Bayesian Theory.  Wiley.
  * Biswas and Mackey (2024) Biswas, N. and L. Mackey (2024).  Bounding Wasserstein distance with couplings.  Journal of the American Statistical Association.  DOI: 10.1080/01621459.2023.2287773.
  * Blei et al. (2017) Blei, D. M., A. Kucukelbir, and J. D. McAuliffe (2017).  Variational inference: A review for statisticians.  Journal of the American Statistical Association 112, 859–877.
  * Bradbury et al. (2018) Bradbury, J., R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang (2018).  JAX: composable transformations of Python+NumPy programs.
  * Burbea (1984) Burbea, J. (1984).  The convexity with respect to Gaussian distributions of divergences of order $\alpha$.  Utilitas Mathematica 26, 171–192.
  * Cai et al. (2024) Cai, D., C. Modi, L. Pillaud-Vivien, C. C. Margossian, R. M. Gower, D. M. Blei, and L. K. Saul (2024).  Batch and match: black-box variational inference with a score-based divergence.  arXiv:2402.14758.
  * Cichocki and Amari (2010) Cichocki, A. and S.-i. Amari (2010).  Families of alpha- beta- and gamma-divergences: Flexible and robust measures of similarities.  Entropy 12, 1532–1568.
  * Courtade (2016) Courtade, T. A. (2016).  Links between the logarithmic Sobolev inequality and the convolution inequalities for entropy and Fisher information.  arXiv:1608.05431.
  * Cover and Thomas (2006) Cover, T. M. and J. A. Thomas (2006).  Elements of Information Theory.  John Wiley & Sons, Inc.
  * Daudel et al. (2023) Daudel, K., J. Benton, and A. Doucet (2023).  Alpha-divergence variational inference meets importance weighted auto-encoders: Methodology and asymptotics.  Journal of Machine Learning Research 24, 1–83.
  * Daudel et al. (2021) Daudel, K., R. Douc, and F. Roueff (2021).  Infinite-dimensional gradient-based descent for alpha-divergence minimisation.  The Annals of Statistics 49, 2250–2270.
  * Dhaka et al. (2021) Dhaka, A. K., A. Catalina, M. Welandawe, M. R. Andersen, J. H. Huggins, and A. Vehtari (2021).  Challenges and opportunities in high-dimensional variational inference.  In Advances in Neural Information Processing Systems 34, pp. 7787–7798.
  * Dieng et al. (2017) Dieng, A. B., D. Tran, R. Ranganath, J. Paisly, and D. M. Blei (2017).  Variational inference via $\chi$ upper bound minimization.  In Advances in Neural Information Processing Systems 30, pp. 2732–2741.
  * Dua and Graff (2017) Dua, D. and C. Graff (2017).  UCL machine learning repository.
  * Geffner and Domke (2021) Geffner, T. and J. Domke (2021).  On the difficulty of unbiased alpha divergence minimization.  In Proceedings of the 38th International Conference on Machine Learning, pp. 3650–3659.
  * Gelman et al. (2013) Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2013).  Bayesian Data Analysis.  Chapman & Hall/CRC Texts in Statistical Science.
  * Gelman and Hill (2007) Gelman, A. and J. Hill (2007).  Data Analysis Using Regression and Multilevel-Hierarchical Models.  Cambridge University Press.
  * Gil et al. (2013) Gil, M., F. Alajaji, and T. Linder (2013).  Rényi divergence measures for commonly used univariate continuous distributions.  Information Sciences 249, 124–131.
  * Giordano et al. (2018) Giordano, R., T. Broderick, and M. I. Jordan (2018).  Covariances, robustness, and variational bayes.  Journal of Machine Learning Research 19, 1–49.
  * Gronau et al. (2020) Gronau, Q. F., H. Singmann, and E.-J. Wagenmakers (2020).  bridgesampling: An R package for estimating normalizing constants.  Journal of Statistical Software 92.
  * Hernandez-Lobato et al. (2016) Hernandez-Lobato, J., Y. Li, M. Rowland, T. Bui, D. Hernández-Lobato, and R. Turner (2016).  Black-box alpha divergence minimization.  In Proceedings of the 33rd International Conference on Machine Learning, pp. 1511–1520.
  * Hobza et al. (2009) Hobza, T., D. Morales, and L. Pardo (2009).  Statistical Methodology 6, 424–436.
  * Hoffman and Blei (2015) Hoffman, M. D. and D. M. Blei (2015).  Stochastic structured variational inference.  In Proceedings of the 18th International Conference on Artificial Intelligence and Statistics, pp. 361–369.
  * Horn and Johnson (2012) Horn, R. A. and C. R. Johnson (2012).  Matrix analysis.  Cambridge University Press.
  * Huggins et al. (2020) Huggins, J. H., M. Kasprzak, T. Campbell, and T. Broderick (2020).  Validated variational inference via practical posterior error bounds.  In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics, pp. 1792–1802.
  * Johnson et al. (2016) Johnson, M. J., D. K. Duvenaud, A. Wiltschko, R. P. Adams, and S. R. Datta (2016).  Composing graphical models with neural networks for structured representations and fast inference.  In Advances in neural information processing systems 29, pp. 2946–2954.
  * Jordan et al. (1999) Jordan, M. I., Z. Ghahramani, T. S. Jaakkola, and L. K. Saul (1999).  An introduction to variational methods for graphical models.  Machine Learning 37, 183–233.
  * Kim et al. (1998) Kim, S., N. Shepard, and S. Chib (1998).  Stochastic volatility: likelihood inference and comparison with arch models.  The Review of Economic Studies 65, 361–393.
  * Kucukelbir et al. (2017) Kucukelbir, A., D. Tran, R. Ranganath, A. Gelman, and D. Blei (2017).  Automatic differentiation variational inference.  Journal of Machine Learning Research 18, 1–45.
  * Li and Turner (2016) Li, Y. and R. E. Turner (2016).  Rényi divergence variational inference.  In Advances in Neural Information Processing Systems 29, pp. 1073–1081.
  * Liese and Vajda (1987) Liese, F. and I. Vajda (1987).  Convex Statistical Distances.  Leipzig: Teubner.
  * MacKay (2003) MacKay, D. J. (2003).  Information theory, inference, and learning algorithms.
  * Margossian and Saul (2023) Margossian, C. C. and L. K. Saul (2023).  The shrinkage-delinkage trade-off: An analysis of factorized gaussian approximations for variational inference.  In Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence, pp. 1358–1367.
  * Meng and Schilling (2002) Meng, X. and S. Schilling (2002).  Warp bridge sampling.  Journal of Computational and Graphical Statistics 11, 552–586.
  * Minka (2005) Minka, T. (2005).  Divergence measures and message passing.  Technical Report MSR-TR-2005-173.
  * Nordström (2011) Nordström, K. (2011).  Convexity of the inverse and Moore-Penrose inverse.  Linear Algebra and Its Applications 434, 1489–1512.
  * Opper and Archambeau (2009) Opper, M. and C. Archambeau (2009).  The variational gaussian approximation revisited.  Neural Computation 21(3), 786–792.
  * Parisi (1988) Parisi, G. (1988).  Statistical Field Theory.  Addison-Wesley.
  * Robert and Casella (2004) Robert, C. P. and G. Casella (2004).  Monte Carlo Statistical Methods.  Springer.
  * Rosenbrock (1960) Rosenbrock, H. H. (1960).  An automatic method for finding the greatest or least value of a function.  Computer Journal 3, 175–184.
  * Rubin (1981) Rubin, D. B. (1981).  Estimation in parallelized randomized experiments.  Journal of Educational Statistics 6, 377–400.
  * Sountsov et al. (2020) Sountsov, P., A. Radul, and contributors (2020).  Inference gym.
  * Tomczak (2022) Tomczak, J. M. (2022).  Deep generative modeling.  Springer.
  * Turner and Sahani (2011) Turner, R. E. and M. Sahani (2011).  Two problems with variational expectation maximisation for time-series models.  In D. Barber, A. T. Cemgil, and S. Chiappa (Eds.), Bayesian Time series models, Chapter 5, pp. 109–130. Cambridge University Press.
  * Vehtari et al. (2020) Vehtari, A., A. Gelman, T. Sivula, P. Jylanki, D. Tran, S. Sahai, P. Blomstedt, J. P. Cunningham, Schiminovich, and C. P. Robert (2020).  Expectation propagation as a way of life: A framework for bayesian inference on partitioned data.  Journal of Machine Learning 21, 1–53.
  * Vehtari et al. (2024) Vehtari, A., D. Simpson, A. Gelman, Y. Yao, and J. Gabry (2024).  Pareto smoothed importance sampling.  Journal of Machine Learning Research.  To appear.
  * Wainwright and Jordan (2008) Wainwright, M. J. and M. I. Jordan (2008).  Foundations and Trends in Machine Learning 1(1–2), 1–305.
  * Wang and Blei (2018) Wang, Y. and D. M. Blei (2018).  Frequentist consistency of variational bayes.  Journal of the American Statistical Association 114, 1147–1161.
  * Yao et al. (2018) Yao, Y., A. Vehtari, D. Simpson, and A. Gelman (2018).  Yes, but did it work?: Evaluating variational inference.  In Proceedings of the 35th International Conference on Machine Learning, pp. 5577–5586.
  * Zhang et al. (2022) Zhang, L., B. Carpenter, A. Gelman, and A. Vehtari (2022).  Pathfinder: Parallel quasi-Newton variational inference.  Journal of Machine Learning Research 23(306), 1–49.

##  Appendix A Additional details for experiments

###  A.1 Batch and Match algorithm for FG-VI

Batch and Match (BaM) (Cai et al., 2024, BaM) is an iterative algorithm for VI
that attempts to estimate and minimize the reverse score-based divergence,
$\text{S}(q||p)$. At each iteration, BaM minimizes a regularized objective
function

$\mathcal{L}^{\text{BaM}}(q)=\widehat{\mathcal{D}}_{q_{t}}(q||p)+\frac{2}{\lambda_{t}}\text{KL}(q||q_{t}),$
(79)

where

  * •

$q_{t}\in\mathcal{Q}$ is a current approximation,

  * •

$\hat{\mathcal{D}}_{q_{t}}(q||p)$ is a Monte Carlo estimator of the score-
based divergence

$\hat{\mathcal{D}}_{q_{t}}(q||p)=\frac{1}{B}\sum_{b=1}^{B}||\nabla\log
q\big{(}{\bf z}^{(b)}\big{)}-\nabla\log p\big{(}{\bf
z}^{(b)}\big{)}||^{2}_{\boldsymbol{\Psi}},$ (80)

using draws ${\bf z}^{(b)}\sim q_{t}$.

  * •

$\lambda_{t}>0$ is the learning rate, or step size, at the $t^{\text{th}}$
iteration.

Cai et al. (2024, Appendix C) derive the BaM updates which minimize eq. 79
when $\mathcal{Q}$ is the family of Gaussian distributions with dense
covariance matrices. In this appendix, we derive the analogous updates to
minimize eq. 79 when $\mathcal{Q}$ is the family of Gaussian distributions
with diagonal covariance matrices. These were the updates used for the
experiments in Section 5 of the paper.

We follow closely the derivation in Cai et al. (2024), noting only the
essential differences for the case of FG-VI. We use ${\bf g}^{(b)}=\nabla\log
p\big{(}{\bf z}^{(b)}\big{)}$ as shorthand for the score at the
$b^{\text{th}}$ sample. As before, the following averages need to be computed
at each iteration:

$\bar{\bf z}=\frac{1}{B}\sum_{b=1}^{B}{\bf z}^{(b)},\quad\bar{\bf
g}=\frac{1}{B}\sum_{b=1}^{B}{\bf g}^{(b)}.$ (81)

For FG-VI, we instead compute diagonal matrices ${\bf C}$ and
$\boldsymbol{\Gamma}$ with nonnegative elements

$C_{ii}=\frac{1}{B}\sum_{b=1}^{B}\left(z_{i}^{(b)}-\bar{z}_{i}\right)^{2},\ \
\Gamma_{ii}=\frac{1}{B}\sum_{b=1}^{B}\left(g^{(b)}_{i}-\bar{g}_{i}\right)^{2}.$
(82)

Let $q_{t+1}$ denote the minimizer of eq. 79, with mean
$\boldsymbol{\nu}_{t+1}$ and diagonal covariance matrix
$\boldsymbol{\Psi}^{t+1}$. The update for the mean takes the same form as in
the original formulation of BaM:

$\boldsymbol{\nu}^{t+1}=\frac{\lambda_{1}}{1+\lambda_{t}}\left(\bar{\bf
z}+\boldsymbol{\Psi}^{t+1}\bar{\bf
g}\right)+\frac{1}{1+\lambda_{t}}\boldsymbol{\nu}^{t}.$ (83)

To obtain the covariance update, we substitute this result into eq. 79 and
differentiate with respect to the diagonal elements of
$\boldsymbol{\Psi}^{t+1}$. For FG-VI, we obtain a simple quadratic equation
for each updated variance:

$\left(\Gamma_{ii}+\frac{1}{1+\lambda_{t}}\bar{g}_{i}^{2}\right)\left(\Psi_{ii}^{t+1}\right)^{2}+\frac{1}{\lambda_{t}}\Psi^{t+1}_{ii}-\left(C_{ii}+\frac{1}{\lambda_{t}}\Psi^{t}_{ii}+\frac{(\nu^{t}_{i}-\bar{z}_{i})^{2}}{1+\lambda_{t}}\right)=0.$
(84)

Eq. (84) always admits a positive root, and for FG-VI, this positive root is
the BaM update for the $i^{\text{th}}$ diagonal element of
$\boldsymbol{\Psi}^{t+1}$.

At a high level, BaM works the same way with diagonal covariance matrices as
it does with full covariance matrices. Each iteration involves a batch step
which draws the $B$ samples from $q_{t}$ and a match step that performs the
calculations in eqs. (81-84). We find a constant learning rate
$\lambda_{t}\\!=\\!1$ to work well for the experiments in Section 5.

###  A.2 Models from the inference gym

In this appendix, we provide details about the models and data sets from the
inference gym used in Section 5.

#### Rosenbrock Distribution.

(n = 2) A transformation of a normal distribution, with non-trivial
correlations (Rosenbrock, 1960). The contour plots of the density function
have the shape of a crescent moon. The joint distribution is given by:

$\displaystyle p(z_{1})$ $\displaystyle=\text{normal}(0,10)$ (85)
$\displaystyle p(z_{2}|z_{1})$
$\displaystyle=\text{normal}(0.03(z_{1}^{2}-100),1).$ (86)

We use FG-VI to approximate $p(z_{1},z_{2})$.

#### Eight Schools.

($n\\!=\\!10$) A Bayesian hierarchical model of the effects of a test
preparation program across 8 schools (Rubin, 1981):

$\displaystyle p(\mu)$ $\displaystyle=\text{normal}(5,3^{2})$ (87)
$\displaystyle p(\tau)$ $\displaystyle=\text{normal}^{+}(0,10)$ (88)
$\displaystyle p(\theta_{i}|\mu,\sigma)$
$\displaystyle=\text{normal}(\mu,\tau^{2})$ (89) $\displaystyle
p(y_{i}|\theta_{i})$ $\displaystyle=\text{normal}(\theta_{i},\sigma_{i}^{2}),$
(90)

where $\text{normal}^{+}$ is a normal distribution truncated at 0. The
observations are ${\bf x}=(y_{1:8},\sigma_{1:8})$ and the latent variables are
${\bf z}=(\mu,\tau,\theta_{1:8})$. We use VI to approximate the posterior
distribution $p(\mathbf{z}|\mathbf{x})$.

#### German Credit.

($n\\!=\\!25$) A logistic regression applied to a credit data set (Dua and
Graff, 2017), with covariates ${\bf x}=x_{1:24}$, observations ${\bf
y}=y_{1:24}$, coefficients ${\bf z}=z_{1:24}$ and an intercept $z_{0}$:

$\displaystyle p({\bf z})$ $\displaystyle=\text{normal}(0,{\bf I}_{n}),$ (91)
$\displaystyle p(y_{i}|{\bf z},{\bf x})$
$\displaystyle=\text{Bernoulli}\left(\frac{1}{e^{-{\bf z}^{T}{\bf
x}_{i}-z_{0}}}\right)$ (92)

We use FG-VI to approximate $p(\mathbf{z}|y_{i},\mathbf{x})$.

#### Radon Effect.

($n\\!=\\!91$) A hierarchical model to measure the Radon level in households
(Gelman and Hill, 2007). Here, we restrict ourselves to data from Minnesota.
For each household, we have three covariates, ${\bf x}_{i}=x_{i,1:3}$, with
corresponding coefficients ${\bf w}=w_{1:3}$, and a county level covariate,
$\theta_{j[i]}$, where $j[i]$ denotes the county of the $i^{\text{th}}$
household. The model is:

$\displaystyle p(\mu)$ $\displaystyle=\text{normal}(0,1)$ (93) $\displaystyle
p(\tau)$ $\displaystyle=\text{Uniform}(0,100)$ (94) $\displaystyle
p(\theta_{j}|\mu,\tau)$ $\displaystyle=\text{normal}(\mu,\tau^{2})$ (95)
$\displaystyle p(w_{k})$ $\displaystyle=\text{normal}(0,1)$ (96)
$\displaystyle p(\sigma)$ $\displaystyle=\text{Uniform}(0,100)$ (97)
$\displaystyle p(\log y_{i})$ $\displaystyle=\text{normal}({\bf w}^{T}{\bf
x}_{i}+\theta_{j[i]},\sigma^{2}).$ (98)

We use FG-VI to approximate
$p(\mu,\tau,\mathbf{\theta},\mathbf{w},\sigma|\mathbf{x},y)$.

#### Stochastic volatility.

(n=103) A time series model with 100 observations. The original model by Kim
et al. (1998) is given by:

$\displaystyle p(\sigma)$ $\displaystyle=\text{Cauchy}^{+}(0,2)$ (99)
$\displaystyle p(\mu)$ $\displaystyle=\text{exponential}(1)$ (100)
$\displaystyle p((\phi+1)/2)$ $\displaystyle=\text{Beta}(20,1.5)$ (101)
$\displaystyle p(z_{i})$ $\displaystyle=\text{normal}(0,1)$ (102)
$\displaystyle h_{1}$ $\displaystyle=\mu+\sigma z_{1}/\sqrt{1-\phi^{2}}$ (103)
$\displaystyle h_{i>1}$ $\displaystyle=\mu+\sigma z_{i}+\phi(h_{i-1}-\mu)$
(104) $\displaystyle p(y_{i}\mid h_{i})$
$\displaystyle=\text{normal}(0,\exp(h_{i}/2)).$ (105)

The original model was developed for a time series of 3000 observations, but
we only work with the first 100 observations. For these observations, we use
FG-VI to model $p(\sigma,\mu,\phi,\mathbf{z}|\mathbf{y})$.