# The Shrinkage-Delinkage Trade-off: An Analysis of
Factorized Gaussian Approximations for Variational Inference

Charles C. Margossian Center for Computational Mathematics
Flatiron Institute
New York, NY, USA  Lawrence K. Saul  Center for Computational Mathematics
Flatiron Institute
New York, NY, USA

###### Abstract

When factorized approximations are used for variational inference (VI), they
tend to understimate the uncertainty—as measured in various ways—of the
distributions they are meant to approximate. We consider two popular ways to
measure the uncertainty deficit of VI: (i) the degree to which it
underestimates the componentwise variance, and (ii) the degree to which it
underestimates the entropy. To better understand these effects, and the
relationship between them, we examine an informative setting where they can be
explicitly (and elegantly) analyzed: the approximation of a Gaussian, $p$,
with a dense covariance matrix, by a Gaussian, $q$, with a diagonal covariance
matrix. We prove that $q$ always underestimates both the componentwise
variance and the entropy of $p$, though not necessarily to the same degree.
Moreover we demonstrate that the entropy of $q$ is determined by the trade-off
of two competing forces: it is decreased by the shrinkage of its componentwise
variances (our first measure of uncertainty) but it is increased by the
factorized approximation which delinks the nodes in the graphical model of
$p$. We study various manifestations of this trade-off, notably one where, as
the dimension of the problem grows, the per-component entropy gap between $p$
and $q$ becomes vanishingly small even though $q$ underestimates every
componentwise variance by a constant multiplicative factor. We also use the
shrinkage-delinkage trade-off to bound the entropy gap in terms of the problem
dimension and the condition number of the correlation matrix of $p$. Finally
we present empirical results on both Gaussian and non-Gaussian targets, the
former to validate our analysis and the latter to explore its limitations.

##  1 Introduction

Variational inference (VI) is a popular methodology for approximate Bayesian
inference [Jordan et al., 1999, Wainwright and Jordan, 2008, Blei et al.,
2017]. Given a target distribution, $p$, VI searches for a tractable
distribution, $q\in\mathcal{Q}$, that minimizes the Kullback-Leibler (KL)
divergence to $p$. A common choice for $\mathcal{Q}$ is to use a family of
factorized distributions. The KL-divergence can then be optimized in a
scalable manner for high-dimensional distributions [Wainwright and Jordan,
2008], which is crucial, for instance, to train models such as variational
auto-encoders over large data sets [Kingma and Welling, 2013].

Factorized VI has its roots in the mean-field approximations to certain Gibbs
distributions from statistical physics [Parisi, 1988, MacKay, 2003]. In this
approach, the approximating distribution is modeled as

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

In most applications, the target distribution $p({\bf z})$ does not factorize.
By its very nature, factorized VI cannot estimate the correlations between
different elements of ${\bf z}$. A more subtle shortcoming of factorized VI is
that it also fails to correctly estimate the marginal distributions,
$p(z_{i})$. This failure typically manifests as an approximation $q$ with an
uncertainty deficit relative to $p$, a phenomenon which has been studied both
empirically and theoretically [e.g MacKay, 2003, Wang and Titterington, 2005,
Bishop, 2006, Turner and Sahani, 2011, Blei et al., 2017, Giordano et al.,
2018]. There exists several measures of uncertainty, and we focus on two: (i)
the componentwise variance and (ii) the entropy. The componentwise variance
plays a crucial role in Bayesian modeling, especially when estimating the
posterior distribution over interpretable variables. Meanwhile the entropy
provides a multivariate notion of uncertainty and, in statistical physics, can
be linked to the free energy, a quantity of interest for many problems.

Intuitively, we expect factorized VI to shrink the variance of $q$ to minimize
its overlap with the tails of $p$. It is less clear how it should affect the
entropy of $q$: on the one hand, this entropy is decreased by any shrinkage in
the variance, but it is increased by the factorized approximation, which
delinks the nodes in the full-covariance graphical model of $p$. Hence entropy
is driven by a trade-off between two competing forces. We call this the
shrinkage-delinkage trade-off. This trade-off hints that the adequacy of
factorized VI may depend on the way we elect to measure its uncertainty
deficit.

The goal of this paper is to understand the uncertainty deficit of factorized
VI in the most informative setting where it can be rigorously analyzed. To
this end, we study the special case where $p$ is a Gaussian distribution over
$\mathbb{R}^{n}$ with a full covariance matrix and $q$ is a Gaussian
distribution over $\mathbb{R}^{n}$ with a diagonal covariance matrix. This
choice of $q$ is natural when $p$ is a multivariate distribution over ${\bf
z}\in\mathbb{R}^{n}$, and leads to factorized Gaussian variational inference
(FG-VI)—a popular method among practionners due notably to “black box”
implementations such as automatic differentiation variational inference (ADVI)
[Kucukelbir et al., 2017]. Our paper expands on previous analyses of FG-VI
[Bishop, 2006, Turner and Sahani, 2011] in many ways, but perhaps most
significantly by identifying—and elucidating—the shrinkage-delinkage trade-off
of factorized VI, which in the considered setting can be written explicitly.

Our analysis is grounded in two fundamental inequalities. First we show that
if $p$ is multivariate Gaussian, and if $q$ is the distribution (optimally)
estimated by FG-VI, then

$\text{Var}_{q}(z_{i})\leq\text{Var}_{p}(z_{i}).$ (2)

Second, under the same assumptions, we show that

${\cal H}(q)\leq{\cal H}(p),$ (3)

where ${\cal H}(\cdot)$ denotes the entropy. This second inequality, relating
the entropies of $p$ and $q$, formalizes an observation [MacKay, 2003, Bishop,
2006] that $q$ tends to be more “compact” than $p$. Our proofs of these
inequalities hold generally for Gaussian distributions over $\mathbb{R}^{n}$;
to the best of our knowledge, they are more direct and more general than
previous demonstrations While both inequalities reveal an uncertainty deficit,
we will see that the two notions of uncertainty are not equivalent. Indeed, we
provide one example where $q$ underestimates each componentwise variance by a
constant multiplicative factor, but the per-component entropy gap between $p$
and $q$ can be arbitrarily small. This discrepancy arises because the entropy
gap in FG-VI is in fact equal to the KL divergence minimized by FG-VI when it
targets a multivariate Gaussian. But, as we will see, this choice of objective
function can harm the estimation of marginal variances.

The inequalities in eq. (2–3) anchor our subsequent analysis. As shown in
Figure 1, the amount of shrinkage in FG-VI depends in general on the number of
components of ${\bf z}\in\mathbb{R}^{n}$ as well as the degree of correlation
between these components. With this motivation, we derive an upper bound on
the entropy gap in eq. (3) in terms of the problem dimensionality, $n$, and
the condition number of the true correlation matrix.

Finally we examine the relevance for some of our findings when FG-VI is
applied to non-Gaussian target distributions. For these experiments, we draw
on several examples from the Bayesian literature. We find that, while the
variance shrinkage (2) does not hold systematically, it holds on average in
the considered examples. We do not have a reliable method to empirically
estimate the entropy, but make an argument that eq. (3) may hold in the
studied examples.

Our results build on those of many previous studies. MacKay [2003], Bishop
[2006], Turner and Sahani [2011], and Blei et al. [2017] all use a two-
dimensional Gaussian to illustrate that the approximations from VI are more
“compact” than the distributions they target. We formalize this observation in
the general $n$-dimensional setting, while highlighting the difference between
componentwise variance and entropy as measures of uncertainty—a difference
that becomes more critical in high-dimensional settings. Experiments on non-
Gaussian models also suggest a more nuanced picture, showing for instance that
FG-VI does not always underestimate every componentwise variance, though in
the studied examples variance shrinkage holds on average. Previous studies
have also examined other measure of uncertainty, such as the frequentist
intervals obtained by variational Bayes estimators [e.g. Wang and
Titterington, 2005]. Finally, many have been motivated by the uncertainty
deficit of factorized VI to develop new methods for inference. These include
post-hoc corrections of variational approximations [Giordano et al., 2018] or,
for certain models, careful decompositions of $p$ using conditional
distributions to justify the assumption of factorization [Agrawal and Domke,
2021].

The code for all results and figures is available on GitHub.

Figure 1: FG-VI’s approximation of a multivariate Gaussian whose correlation
matrix has constant off-diagonal terms. FG-VI’s variance shrinkage grows with
both increasing dimensionality ($n$) and correlation ($\varepsilon$). For
$n=64$, the distributions are projected onto their first two coordinates.
Despite what the picture suggests, the entropy gap between the approximation
and the target is actually quite small (Section 3.2). In this sense, the
lower-dimensional projection is misleading.

##  2 Preliminaries

We analyze FG-VI in the setting where $p(\mathbf{z})$ is multivariate Gaussian
with mean $\boldsymbol{\mu}\in\mathbb{R}^{n}$ and covariance
$\boldsymbol{\Sigma}\in\mathbb{R}^{n\times n}$. In this setting FG-VI has a
particularly simple solution. (An earlier statement of this solution can be
found in Turner and Sahani [2011].)

######  Proposition 2.1.

Let $q(\mathbf{z})$ be multivariate Gaussian with mean $\boldsymbol{\nu}$ and
diagonal covariance $\boldsymbol{\Psi}$. Then the variational parameters
minimizing $\text{KL}(q||p)$ are given by $\boldsymbol{\nu}=\boldsymbol{\mu}$
and

$\Psi_{ii}=\frac{1}{\Sigma^{-1}_{ii}},$ (4)

where the denominator $\Sigma_{ii}^{-1}$ denotes a diagonal element of the
matrix inverse $\boldsymbol{\Sigma}^{-1}$.

###### Proof.

The variational parameters $\boldsymbol{\nu}$ and $\boldsymbol{\Psi}$ are
estimated by minimizing the KL-divergence

$\text{KL}(q||p)=\mathbb{E}_{q}[\log q({\bf z})]-\mathbb{E}_{q}[\log p({\bf
z})],$ (5)

where each expectation is taken with respect to the measure $q$. Note that
only the second term in eq. (5) depends on the variational mean
$\boldsymbol{\nu}$, and it is given by

$-\mathbb{E}_{q}[\log p({\bf
z})]=\tfrac{1}{2}(\boldsymbol{\nu}\\!-\\!\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\nu}\\!-\boldsymbol{\mu})\,+\,\ldots$
(6)

where the ellipses indicate terms that do not depend on $\boldsymbol{\nu}$. By
minimizing this expression, it follows at once that
$\boldsymbol{\nu}=\boldsymbol{\mu}$. With this substitution, eq. (5)
simplifies to

$\text{KL}(q||p)=\tfrac{1}{2}\left[\text{trace}\big{(}\boldsymbol{\Psi}\boldsymbol{\Sigma}^{-1}\big{)}-\log\big{|}\boldsymbol{\Psi}\boldsymbol{\Sigma}^{-1}\big{|}-n\right],$
(7)

and the result in eq. (4) follows by minimizing the above expression with
respect to the diagonal elements of $\boldsymbol{\Psi}$. ∎

In sections 2, 3, and 4 of the paper, we assume that $q$ is the factorized
Gaussian distribution whose variances are given by eq. (4). We emphasize in
general that $\Psi_{ii}\neq\Sigma_{ii}$. However, it is true that
$\boldsymbol{\Psi}=\boldsymbol{\Sigma}$ when $\boldsymbol{\Sigma}$ is
diagonal.

Many of our results will not be expressed directly in terms of
$\boldsymbol{\Sigma}$ and $\boldsymbol{\Psi}$, but in terms of two related
(but dimensionless) matrices. The first is the correlation matrix $\mathbf{C}$
with elements

$C_{ij}=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}.$ (8)

Note that $C_{ii}=1$, a simple fact that we will often exploit, and also that
$\mathbf{C}$ reduces to the identity matrix when $\boldsymbol{\Sigma}$ is
diagonal. At the other extreme, we may consider the case where all the off-
diagonal elements of $\mathbf{C}$ are equal to some constant
$\varepsilon\\!>\\!0$. This is explored visually in Figure. (1). In appendix A
we show that $\Psi_{ii}\rightarrow 0$ as $\varepsilon\rightarrow 1$ for fixed
$n$, and that $\Psi_{ii}\rightarrow(1\\!-\\!\varepsilon)\Sigma_{ii}$ as
$n\rightarrow\infty$ for fixed $\varepsilon$. Note that FG-VI underestimates
the variance in both limits.

In addition to the correlation matrix, we also define the diagonal shrinkage
matrix $\mathbf{S}$ with dimensionless entries

$S_{ii}=\frac{\Sigma_{ii}}{\Psi_{ii}}=\Sigma_{ii}\,\Sigma^{-1}_{ii}.$ (9)

We will use the matrices $\mathbf{C}$ and $\mathbf{S}$ to analyze how FG-VI
underestimates the uncertainty of $p$. The uncertainty in axis-aligned
directions is measured by the variances $\Sigma_{ii}$, but a multivariate
measure of uncertainty is provided by the entropy

$\mathcal{H}(p)=-\mathbb{E}_{p}\log p({\bf z}).$ (10)

To what extent does FG-VI underestimate this entropy? As shown next, the
answer is very naturally expressed in terms of the correlation matrix
$\mathbf{C}$ and the shrinkage matrix $\mathbf{S}$.

######  Proposition 2.2.

Let $p$ and $q$ be defined as above. Then their difference in entropy is given
by

$\mathcal{H}(p)-\mathcal{H}(q)=\tfrac{1}{2}\log|\mathbf{S}|-\tfrac{1}{2}\log|\mathbf{C}|^{-1}.$
(11)

###### Proof.

A standard calculation for multivariate Gaussian distributions [Cover and
Thomas, 2006] gives $\mathcal{H}(p)=\frac{1}{2}\log|\Sigma|(2\pi e)^{n}$, and
an analogous result holds for $\mathcal{H}(q)$. Let
$\Delta\mathcal{H}=\mathcal{H}(p)-\mathcal{H}(q)$. Then we see that

$\Delta\mathcal{H}=\tfrac{1}{2}\log|\mathbf{\Sigma}|-\tfrac{1}{2}|\log{\mathbf{\Psi}}|=\tfrac{1}{2}\log\left|\mathbf{\Psi}^{-\frac{1}{2}}\mathbf{\Sigma}\mathbf{\Psi}^{-\frac{1}{2}}\right|.$
(12)

Now from the definitions in eqs. (8–9), it can be verified by direct
substitution that
$\boldsymbol{\Psi}^{-\frac{1}{2}}\boldsymbol{\Sigma}\,\boldsymbol{\Psi}^{-\frac{1}{2}}=\mathbf{S}^{\frac{1}{2}}\mathbf{C}\,\mathbf{S}^{\frac{1}{2}}$.
It follows from the basic properties of determinants that

$\Delta\mathcal{H}=\tfrac{1}{2}\log\big{|}\mathbf{S}^{\frac{1}{2}}\mathbf{C}\mathbf{S}^{\frac{1}{2}}\big{|}=\tfrac{1}{2}\log|\mathbf{S}|-\tfrac{1}{2}\log|\mathbf{C}|^{-1}.$

∎

##  3 Shrinkage-Delinkage Trade-off

In this section we prove that FG-VI systematically underestimates the variance
and entropy of a multivariate Gaussian distribution. We will see, however,
that a large shrinkage in all componentwise variances does not imply a
correspondingly large shrinkage in the entropy.

###  3.1 Fundamental inequalities

######  Theorem 3.1 (Variance shrinkage).

The solution for FG-VI in eq. (4) underestimates the variance; that is,

$\Psi_{ii}\leq\Sigma_{ii},$ (13)

and the inequality is strict for some component of the variance (i.e.,
$\Psi_{ii}<\Sigma_{ii}$) if $\mathbf{\Sigma}$ is not purely diagonal.

###### Proof.

Let $\mathbf{C}$ denote the correlation matrix in eq. (8). It can be verified
by direct calculation that

$C^{-1}_{ij}=\Sigma^{-1}_{ij}\sqrt{\Sigma_{ii}\Sigma_{jj}}.$ (14)

As further notation, let $\lambda_{1},\ldots,\lambda_{n}$ denote the
eigenvalues of $\mathbf{C}$, and let $\mathbf{e}_{i}$ denote the unit vector
along the $i^{\rm th}$ axis. Then from the solution in eq. (4), it follows
that

$\displaystyle\frac{\Sigma_{ii}}{\Psi_{ii}}$ $\displaystyle=$ $\displaystyle
C^{-1}_{ii},$ (15) $\displaystyle=$ $\displaystyle C^{-1}_{ii}+C_{ii}-1,$ (16)
$\displaystyle=$
$\displaystyle\mathbf{e}_{i}^{\top}(\mathbf{C}^{-1}\\!+\mathbf{C}\,)\,\mathbf{e}_{i}-1,$
(17) $\displaystyle\geq$
$\displaystyle\min_{\|\mathbf{e}\|=1}[\mathbf{e}^{\top}(\mathbf{C}^{-1}\\!+\mathbf{C}\,)\,\mathbf{e}-1],$
(18) $\displaystyle=$
$\displaystyle\min_{i}\,(\lambda_{i}^{-1}+\lambda_{i}-1),$ (19)
$\displaystyle\geq$
$\displaystyle\min_{\lambda>0}\,(\lambda^{-1}+\lambda-1)=1,$ (20)

where in the last step we have used the fact that the correlation matrix
$\mathbf{C}$ has strictly positive eigenvalues. This proves eq. (13). Now
suppose that $\mathbf{\Sigma}$ is not purely diagonal. Then $\mathbf{C}$ is
also not diagonal; hence there must be some unit vector $\mathbf{e}_{i}$ that
is not an eigenvector of $\mathbf{C}$. In this case the inequality in eq. (18)
is strict, showing that $\Sigma_{ii}>\Psi_{ii}$. ∎

Next we examine the difference in entropy given by eq. (11). First we show
that this difference is determined by the trade-off of competing entropic
forces.

######  Theorem 3.2 (The Shrinkage-Delinkage Tradeoff).

Consider the entropy difference in eq. (11) from FG-VI:

$\mathcal{H}(p)-\mathcal{H}(q)=\tfrac{1}{2}\log|\mathbf{S}|-\tfrac{1}{2}\log|\mathbf{C}|^{-1}.$

Both terms in this difference are nonnegative: that is,

$\displaystyle\log|\mathbf{S}|\ $ $\displaystyle\geq\ 0,$ (21)
$\displaystyle\log\frac{1}{|\mathbf{C}|}$ $\displaystyle\geq\ 0.$ (22)

Before proving the theorem we consider the meaning of these inequalities.
Conceptually, the first inequality shows that any shrinkage of variances (from
Theorem 3.1) reduces the entropy of $q$ and thus contributes to a larger
difference in eq. (11). The second inequality shows that the factorization of
$q$ acts as a counterbalance to this effect: the entropy of $p$ is necessarily
reduced by the presence of correlations, but such correlations cannot be
modeled by $q$. Thus the factorization of $q$ must (to some extent) oppose the
entropy difference in eq. (11), and the net difference is determined by the
trade-off of these forces. Visually the factorization of $q$ is represented by
the delinkage of nodes in the full-covariance graphical model for $p$. This is
the essence of the shrinkage-delinkage tradeoff for FG-VI.

###### Proof.

The bound on $\log|\mathbf{S}|$ in eq. (50) follows at once from Theorem 3.1:

$\log|\mathbf{S}|=\sum_{i=1}^{n}\log\frac{\Sigma_{ii}}{\Psi_{ii}}\geq 0.$ (23)

As before, let $\lambda_{1},\ldots,\lambda_{n}$ denote the eigenvalues of
$\mathbf{C}$ so that $\log|\mathbf{C}|=\sum_{i}\log\lambda_{i}$. From Jensen’s
inequality, we have:

$\sum_{i=1}^{n}\log\lambda_{i}\leq
n\log\bigg{[}\tfrac{1}{n}\sum_{i=1}^{n}\lambda_{i}\bigg{]}=n\log\tfrac{1}{n}\,\text{trace}(\mathbf{C})=0,$
(24)

which proves eq. (52). ∎

To prove that $q$ underestimates the entropy of $p$, we need the following
result which is important in its own right.

######  Proposition 3.3.

The entropy gap between $p$ and $q$ is equal to the KL divergence minimized by
FG-VI:

$\mathcal{H}(p)-\mathcal{H}(q)={\rm KL}(q||p).$ (25)

###### Proof.

The identity follows by substituting the solution from eq. (4) into the KL
divergence in eq. (5). This yields the entropy gap, namely ${\rm
KL}(q,p)=\frac{1}{2}\log|\boldsymbol{\Sigma}|-\frac{1}{2}\log|\boldsymbol{\Psi}|$,
computed in eq. (12). ∎

It follows that FG-VI is minimizing the entropy gap between $p$ and $q$ when
it targets a multivariate Gaussian. As suggested by the trade-off in Theorem
3.2, however, the entropy gap can be minimized despite a large shrinkage in
componentwise variances.

The nonnegativity of the KL divergence in eq. (25) also leads to the other
fundamental inequality of this section.

######  Theorem 3.4 (Entropy gap).

The solution for FG-VI in eq. (4) underestimates the entropy; that is,

$\mathcal{H}(q)\leq\mathcal{H}(p),$ (26)

and this inequality is strict if $\mathbf{\Sigma}$ is not purely diagonal.

An immediate implication of this theorem is that the shrinkage term in eq.
(50) dominates the shrinkage-delinkage trade-off in Theorem 3.2.

######  Remark 3.5.

We see also from Theorem 3.4 that $|\mathbf{\Psi}|\leq|\mathbf{\Sigma}|$. This
inequality can be viewed as a multivariate analog of the result, in Theorem
3.1, that $\Psi_{ii}\leq\Sigma_{ii}$.

###  3.2 Demonstration of the trade-off

Figure 2 illustrates the shrinkage-delinkage trade-off in FG-VI, and how it is
resolved, for multivariate Gaussian distributions with two types of covariance
matrices:

  * •

Squared exponential kernel: This type of covariance matrix arises in models
involving Gaussian processes [e.g. Rasmussen and Williams, 2006, Chapter 2].
For this example we sampled a random input ${\bf
x}\sim\text{uniform}(0,200)^{n}$ and set the covariance matrix via the kernel
function $\Sigma_{ij}=\exp(-(x_{i}\\!-\\!x_{j})^{2}/\rho^{2})$. We use the
hyperparameter $\rho\\!>\\!0$ to vary the degree of correlation.

  * •

Constant off-diagonal: The posterior distributions of models with exchangeable
data [Gelman et al., 2013, chapter 5] can generate such covariance matrices,
or at least covariance matrices whose subblocks have the described structure.
For this example we set $\Sigma_{ii}\\!=\\!1$ along the diagonal and
$\Sigma_{ij}\\!=\\!\varepsilon$ for all $i\\!\neq\\!j$, and we used the
hyperparameter $\varepsilon\\!>\\!0$ to vary the degree of correlation.

Figure 2: Shrinkage-delinkage trade-off in FG-VI when the Gaussian target over
$\mathbb{R}^{n}$ has a squared-exponential-kernel covariance matrix (top) or a
covariance matrix with constant off-diagonal terms (bottom). Here
$n\\!=\\!10$.

Figure 2 plots the opposing contributions from the shrinkage and delinkage
terms in eq. (11) using solid and dashed lines. In each panel, the difference
between these curves reveals the degree to which FG-VI underestimates the
entropy of the multivariate Gaussian distribution it is being used to
approximate. It can also be seen that FG-VI manages the shrinkage-delinkage
trade-off differently for different types of covariance matrices. While in the
squared exponential kernel case the entropy gap is large, it is smaller when
the covariance matrix has constant off-diagonal terms: there, the shrinkage
and delinkage terms in eq. (11) are almost perfectly balanced.

This last finding may come as a surprise in light of earlier results, shown in
Figure 1, where the variational approximation is clearly too “compact.” But
the two-dimensional projections in Figure 1 are misleading. In higher
dimensions, the approximating sphere of FG-VI gains more in volume than its
target ellipse; this discrepancy arises because each added component is
independent for $q$ but strongly correlated for $p$. The overall effect is
that the opposing terms in eq. (11) are nearly balanced. Hence even when FG-VI
hardly underestimates the (per-component) entropy, it may still grossly
underestimate the componentwise variance. This contrast becomes more acute in
the asymptotic limit of $n$.

######  Theorem 3.6.

Suppose $\boldsymbol{\Sigma}$ has constant off-diagonal terms,
$\varepsilon>0$. Then the per-component entropy gap vanishes in the limit
$n\rightarrow\infty$, whereas every componentwise variance shrinks by a
constant factor:

$\displaystyle\underset{n\to\infty}{\lim}\
\tfrac{1}{n}\left(\mathcal{H}(q)-\mathcal{H}(q)\right)$ $\displaystyle=0.$
(27) $\displaystyle\underset{n\to\infty}{\lim}\ (\Psi_{ii}/\Sigma_{ii})$
$\displaystyle=1-\varepsilon.$ (28)

The proof is given in Appendix A. The theorem also shows that the average of
the diagonal elements in the shrinkage matrix also converges to a constant
factor:

$\underset{n\to\infty}{\lim}\ \tfrac{1}{n}\mathrm{trace}({\bf
S})=(1-\varepsilon)^{-1}.$ (29)

This example highlights the roots of FG-VI in mean-field approximations from
statistical physics [Parisi, 1988]. As is well known, the mean-field
approximation for the free energy becomes exact in the limit
$n\\!\rightarrow\\!\infty$ for certain spin systems with infinite-range
interactions. The infinite-range interactions in these systems are analogous,
for the Gaussian models we study here, to the assumption of constant off-
diagonal terms in the covariance matrix [e.g Mukherjee et al., 2018,
Margossian and Mukherjee, 2021]. An important takeaway is that factorized
approximations can work well to estimate the entropy—to wit, minimizing the
KL-divergence with FG-VI on a Gaussian target is equivalent to minimizing the
entropy gap—but still fail to accurately compute the componentwise variances.
This can become an important limitation as we apply VI beyond problems in
statistical physics and more broadly to Bayesian modeling, where estimation of
the variances is critical.

##  4 Bounds on ${\bf log}|\mathbf{S}|$ and ${\bf log}|\mathbf{C}|$

In the last section, we saw that the shrinkage-delinkage trade-off played out
differently for different types of covariance matrices; we also proved certain
asymptotic results that depended on the detailed structure of the covariance
matrix (e.g., constant off-diagonal). In this section, we derive more general
bounds on the terms in this trade-off that depend only the problem
dimensionality, $n$, and the condition number, $R$, of the correlation matrix,
${\bf C}$.

###  4.1 Optimizations for upper bounds

Consider the space of all correlation matrices with condition number $R$. We
denote this space by the set

$\mathcal{C}_{R}=\\{\mathbf{C}\in\mathcal{S}^{n}_{+}\,|\,C_{ii}\\!=\\!1\
\forall i,\,\lambda_{\rm max}(\mathbf{C})\\!=\\!R\lambda_{\rm
min}(\mathbf{C})\\}.$ (30)

The set contains the intersection of those $n\\!\times\\!n$ matrices that are
positive semidefinite (i.e., lying in the cone $\mathcal{S}_{+}^{n}$), whose
diagonal elements are equal to unity, and whose largest eigenvalue is $R$
times larger than its smallest one.

If the condition number of the correlation matrix $\mathbf{C}$ is known to be
$R$, then we can (in principle) compute the following upper bounds on the
terms in eq. (11):

$\displaystyle\log|\mathbf{S}|$
$\displaystyle\leq\max_{\mathbf{C}\in\mathcal{C}_{R}}\left[\sum_{i=1}^{n}\log
C_{ii}^{-1}\right],$ (31) $\displaystyle\log|\mathbf{C}|$
$\displaystyle\leq\max_{\mathbf{C}\in\mathcal{C}_{R}}\left[\sum_{i=1}^{n}\log\lambda_{i}(\mathbf{C})\right].$
(32)

In eq. (31), we have used the fact from eq. (15) that
$S_{ii}\\!=\\!C^{-1}_{ii}$, while in eq. (32), we have written the determinant
of a matrix as the product of its eigenvalues.

In practice, however, it is difficult to perform the optimizations over the
set $\mathcal{C}_{R}$ in eq. (31-32). Instead we consider a more tractable
relaxation; the essential idea is to optimize over a larger set of matrices,
one that is characterized only in terms of constraints on its eigenvalues. We
denote this constrained set of eigenvalues by

$\Lambda_{R}=\left\\{\boldsymbol{\lambda}\\!\in\\!\mathbb{R}^{n}_{+}\,\bigg{|}\,\lambda_{1}\\!\geq\\!\ldots\geq\\!\lambda_{n}=R\lambda_{1},\,\sum_{i=1}^{n}\lambda_{i}\\!=\\!n\right\\}\\!.$
(33)

Note that the set ${\cal C}_{R}$ of correlation matrices is contained strictly
within the set of matrices with eigenvalues in $\Lambda_{R}$. In particular, a
matrix in ${\cal C}_{R}$ is constrained to have ones along its diagonal, while
a matrix with eigenvalues in $\Lambda_{R}$ is only constrained to have a trace
equal to $n$. With the above relaxation, we obtain the following upper bounds
on the terms $\log|\mathbf{S}|$ and $\log|\mathbf{C}|$ in eq. (11).

######  Proposition 4.1.

Suppose that the correlation matrix $\mathbf{C}$ has condition number $R$.
Then

$\displaystyle\log|\mathbf{S}|$ $\displaystyle\leq
n\log\frac{1}{n}\\!\left[\max_{\boldsymbol{\lambda}\in\Lambda_{R}}\sum_{i=1}^{n}\lambda_{i}^{-1}\right]$
(34) $\displaystyle\log|\mathbf{C}|$
$\displaystyle\leq\max_{\boldsymbol{\lambda}\in\Lambda_{R}}\left[\sum_{i=1}^{n}\log\lambda_{i}\right].$
(35)

###### Proof.

The second bound is immediate from eq. (32) and the relaxation in eq. (33).
For the first bound, recall that $S_{ii}=C^{-1}_{ii}$, and note from Jensen’s
equality that

$\tfrac{1}{n}\mbox{$\sum_{i}$}\log
C_{ii}^{-1}\leq\log\tfrac{1}{n}\mbox{$\sum_{i}$}\,C^{-1}_{ii}=\log\left[\tfrac{1}{n}\mbox{$\sum_{i}$}\,\lambda_{i}^{-1}(\mathbf{C})\right].$

The bound in eq. (34) follows from the above in concert with the relaxion in
eq. (33). ∎

###  4.2 Solutions from symmetry

The optimizations over $\Lambda_{R}$ in eqs. (34–35) have a great deal of
structure that we can exploit to compute their solutions. We analyze each of
these optimizations in turn.

######  Lemma 4.2.

Let $\boldsymbol{\lambda}\in\Lambda_{R}$ be the solution that maximizes the
right side of eq. (34). Then at most one $\lambda_{i}$ is not equal to either
$\lambda_{1}$ or $\lambda_{n}$.

###### Proof.

We prove the lemma by contradiction. Suppose there exists a solution with
intermediate elements $\lambda_{i}$ and $\lambda_{j}$ that satisfy
$\lambda_{1}\\!>\\!\lambda_{i}\\!>\\!\lambda_{j}\\!>\\!\lambda_{n}$. Consider
the effect on this solution of a perturbation that adds some small amount
$\delta\\!>\\!0$ to $\lambda_{i}$ and subtracts the same amount from
$\lambda_{j}$. Note that for sufficiently small $\delta$, this perturbation
will not leave the set $\Lambda_{R}$; however, it will expand the separation
of $\lambda_{i}$ from $\lambda_{j}$. As a result the objective
$\sum_{i}\lambda_{i}^{-1}$ has a gain

$f(\delta)=\frac{1}{\lambda_{i}+\delta}-\frac{1}{\lambda_{i}}+\frac{1}{\lambda_{j}-\delta}-\frac{1}{\lambda_{j}}.$
(36)

Next we evaluate the derivative $f^{\prime}(\delta)$ at $\delta=0$; doing so
we find $f^{\prime}(0)=\lambda_{j}^{-2}\\!-\\!\lambda_{i}^{-2}>0$. But this
yields a contradiction, because any solution must be maximal, and hence
stationary (i.e., $f^{\prime}(0)\\!=\\!0$), with respect to small
perturbations. ∎

The above lemma greatly restricts the form of the solutions that we must
consider for the optimization in eq. (34). The next lemma does the same for
the optimization in eq. (35).

######  Lemma 4.3.

Let $\boldsymbol{\lambda}\in\Lambda_{R}$ be the solution that maximizes the
right side of eq. (35). Then $\lambda_{i}\\!=\\!\lambda_{j}$ whenever
$1\\!<\\!i\\!<\\!j\\!<\\!n$.

###### Proof.

We prove this lemma in similar fashion. Suppose there exists a solution with
intermediate elements $\lambda_{i}$ and $\lambda_{j}$ that satisfy
$\lambda_{1}\\!\geq\\!\lambda_{i}\\!>\\!\lambda_{j}\\!\geq\\!\lambda_{n}$.
Consider the effect on this solution of a perturbation that adds some small
amount $\delta\\!>\\!0$ to $\lambda_{j}$ and subtracts the same amount from
$\lambda_{i}$. Again, for sufficiently small $\delta$, this perturbation will
not leave the set $\Lambda_{R}$; however, it will diminish the separation of
$\lambda_{i}$ from $\lambda_{j}$. As a result the objective
$\sum_{i}\log\lambda_{i}$ has a gain

$g(\delta)=\log(\lambda_{i}-\delta)+\log(\lambda_{j}+\delta)-\log(\lambda_{i})-\log(\lambda_{j}).$
(37)

Evaluating the derivative, we find
$g^{\prime}(0)=\lambda_{j}^{-1}\\!-\\!\lambda_{i}^{-}>0$. As before this
yields a contradiction, because any solution must be maximal, and hence
stationary (i.e., $g^{\prime}(0)\\!=\\!0$), with respect to small
perturbations. ∎

Figure 3: Upper bounds on the competing terms in the shrinkage-delinkage
tradeoff, eq. (11), for a correlation matrix with condition number $R$. The
actual values for these terms are shown for the correlation matrices
introduced in section 3.2. Note that the entropy gap is minimal at
$n\\!=\\!100$ for the case of a constant off-diagonal correlation matrix.

Now let us consider, at a high level, how these lemmas simplify the
optimizations in eqs. (34–35). The lemmas show that for each optimization,
there exist three elements—the maximum element $\lambda_{1}$, the minimum
element $\lambda_{n}$, and some intermediate element $\lambda_{k}$ for
$1\\!<\\!k\\!<\\!n$—from which the remaining $n\\!-\\!3$ elements of the
solution can be deduced by symmetry. Note also that any solution in
$\Lambda_{R}$ must satisfy the two additional constraints that
$\sum_{i}\lambda_{i}\\!=\\!n$ and $\lambda_{1}\\!=\\!R\lambda_{n}$. In
Appendix B, we show that by exploiting these symmetries and constraints in
concert, we can reduce the optimizations in eqs. (34–35) to a sequence of one-
dimensional problems for which we have closed-form solutions.

Figure 3 plots the bounds on the shrinkage and delinkage terms as a function
of the condition number, $R$, for problems with dimensionalities $n\\!=\\!10$
(left) and $n\\!=\\!100$ (right). These bounds provide envelopes between which
the actual values of the competing terms in eq. (11) must lie. To illustrate
this, the figure also shows the corresponding values of these terms for the
squared-exponential-kernel and constant-off-diagonal covariance matrices
introduced in the previous section. (Notice that in this figure, unlike Figure
2, these values are plotted against the condition number of the correlation
matrix rather than the hyperparameters $\rho$ or $\varepsilon$.)

Using similar methods, it also possible to upper-bound the entropy gap and the
trace of the shrinkage matrix (which reflects the average shrinkage in
componentwise variance). The derivations of these additional bounds are
relegated to Appendices C and D.

##  5 Non-Gaussian models

Do our results extend in any way when FG-VI is applied to non-Gaussian models?
In this section we suppose that $p$ is a non-Gaussian target with covariance
$\boldsymbol{\Sigma}$. Our previous analysis of FG-VI targets was based on the
variance estimator,

$(\Psi_{G})_{ii}:=1/(\Sigma^{-1})_{ii},$

and the corresponding shrinkage matrix with diagonal elements
$(S_{G})_{ii}=\Sigma_{ii}/(\Psi_{G})_{ii}$. But neither $\mathbf{\Psi}_{G}$
nor $\mathbf{S}_{G}$ will be returned by FG-VI when it is applied to a non-
Gaussian target with covariance $\boldsymbol{\Sigma}$.

To explore these issues, we applied ADVI [Kucukelbir et al., 2017] with a
factorized Gaussian approximation to study the posterior distributions in
several Bayesian models as well as one “adversarial” target (Table 1). These
test targets represent a diversity of applications. The GLM and 8-schools
models are taken from the model data base PosteriorDB [Magnusson et al.,
2022], while the disease map Gaussian process model and sparse kernel
interaction model [Agrawal et al., 2019] are studied by [Margossian et al.,
2020]. We also included a mixture of well-separated spherical Gaussians; for
this target, the approximation by FG-VI collapses to one of the modes, so that
FG-VI can underestimate the componentwise variances by an arbitrarily large
amount (e.g., if the modes are widely separated). Note that in all cases,
before applying ADVI, we transformed any constrained (e.g., nonnegative)
variables of the target distribution to be unconstrained variables over
$\mathbb{R}$.

Call | $n$ | Description | $\frac{1}{n}\log|\boldsymbol{\Sigma}/\boldsymbol{\Psi}|$
---|---|---|---
glm_binomial | 3 | General linear model with a binomial likelihood. | 0.291
8schools_nc | 10 | Hierarchical model with a non-centered parameterization. | 0.011
8schools_pool | 9 | Same model but with a small, fixed population variance value to enforce | 0.339
|  | strong partial pooling and create a high posterior correlation. |
disease_map | 102 | Gaussian process model with Poisson likelihood. Applied to disease | 0.066
|  | map of Finland using 100 randomly sampled counties (out of 911). |
SKIM | 305 | Sparse kernel interaction model, applied to a Prostate cancer microarray | 0.033
|  | data set on a subset of 200 SNPS. |
Mixture | 2 | Mixture of well-separated Gaussians with spherical covariance matrices. | 3.051
Table 1: Non-Gaussian targets for numerical experiments.

We estimated the posterior covariance using long runs of Markov chain Monte
Carlo, specifically 16,000 draws using the software Stan [Carpenter et al.,
2017], except in the mixture example, where the covariance was calculated
analytically. We then estimated (i) the shrinkage matrix, $\mathbf{S}$, when
targeting the posterior and (ii) the shrinkage matrix, $\mathbf{S}_{G}$, when
targeting a Gaussian with the same covariance as the posterior. For non-
Gaussian posteriors, we observe that FG-VI does not always underestimate the
componentwise variance; see Figure 4. On the other hand, for all models in
Table 1, we see that $\frac{1}{n}\text{trace}(\mathbf{S})\\!>\\!1$, meaning
that the componentwise variances are underestimated on average (Figure 5). In
addition, for the Bayesian models, we observe that
$\text{trace}(\mathbf{S})\\!\approx\\!\text{trace}(\mathbf{S}_{G})$. The
mixture target, however, provides a counter-example where
$\text{trace}(\mathbf{S})\gg\text{trace}(\mathbf{S}_{G})$.

Figure 4: Shrinkage matrix for FG-VI when targeting the posterior distribution
of 8schools_nc versus targeting a Gaussian with the same covariance matrix.
For the non-Gaussian target, we may have $S_{ii}\\!<\\!1$. Figure 5: Trace of
shrinkage matrix for various models when targeting the true posterior versus
targeting a Gaussian with the same covariance matrix.

In addition to the shrinkage in componentwise variances, we would have also
liked to evaluate the entropy gap in these models. It is easy to obtain an
upper bound on this gap by observing that the Gaussian maximizes the entropy
among all continuous distributions with a given covariance
$\boldsymbol{\Sigma}$ [Cover and Thomas, 2006]. Thus we have

$\mathcal{H}(p)-\mathcal{H}(q)\leq\tfrac{1}{2}(\log|\boldsymbol{\Sigma}|-\log|\boldsymbol{\Psi}|).$
(38)

We observed this upper bound to be positive for all the models in Table 1. But
we also know that FG-VI can overestimate the entropy in non-Gaussian models:
Turner and Sahani [2011] demonstrated this for a mixture of largely
overlapping 1-dimensional Gaussians. It is an open question to understand the
general conditions under which FG-VI underestimates the entropy. We next note
that our upper-bound on the entropy gap does not immediately apply to (38),
because it is unclear how $\log|\boldsymbol{\Psi}|$ and
$\log|\boldsymbol{\Psi}_{G}|$ compare. To evaluate the entropy gap empirically
in a non-Gaussian model—as would be required to investigate this question
further—it is necessary to estimate the normalizing constant of the posterior.
Candidate methods for this, such as bridge sampling [e.g Gronau et al., 2020,
Meng and Schilling, 2002], rely on a proposal distribution which (roughly)
approximates the target. Typically, a Gaussian-like approximation is used for
these proposals, but this is precisely the assumption we want to relax. In
other words, we do not wish to compare a theory for Gaussian targets to an
empirical benchmark which relies on a Gaussian approximation. We leave this
issue to future work.

##  6 Discussion

In this paper we have shown that FG-VI underestimates the componentwise
variance and joint entropy of a multivariate Gaussian distribution.
Furthermore we expressed the entropy gap as a trade-off between two competing
terms and observed that it was equal to the KL divergence minimized by FG-VI.
Our analysis helps to understand why FG-VI can greatly underestimate the
componentwise variances even when it effectively minimizes the entropy gap and
KL divergence. Our results also suggest that better estimates of variance may
be obtained by changing the objective function or using a different family of
approximations.

An open question is whether a shrinkage-delinkage tradeoff operates when FG-VI
is applied to non-Gaussian targets. For such targets, we have produced a
counter-example showing that FG-VI can overestimate a particular componentwise
variance. On the other hand, we have observed that these variances are
underestimated on average, and moreover that the shrinkage term, $\log|{\bf
S}|$, remains positive. It requires further investigation to make more general
statements about the entropy gap. Finally, it would also be interesting to
extend our analysis beyond FG-VI—for instance to approximations based on a
diagonal plus low-rank covariance matrix, rather than a strictly diagonal one
[e.g Zhang et al., 2022].

##  7 Acknowledgments

We thank David Blei, Bob Carpenter, Justin Domke and Chirag Modi for feedback
on this manuscript.

## References

  * Agrawal and Domke [2021] Abhinav Agrawal and Justin Domke.  Amortized variational inference in simple hierarchical models.  _Advances in Neural Information Processing Systems (NeurIPS 2021)_ , 35, 2021.
  * Agrawal et al. [2019] Raj Agrawal, Jonathan H Huggins, Brians Trippe, and Tamara Broderick.  The Kernel interaction trick: Fast Bayesian discovery of pairwise interactions in high dimensions.  _Proceedings of the 36th International Conference on Machine Learning_ , 97, April 2019.
  * Bishop [2006] Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Blei et al. [2017] David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American Statistical Association_ , 112, 2017.  10.1080/01621459.2017.1285773.  URL https://arxiv.org/abs/1601.00670.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matt Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus A. Brubaker, Jiqiang Guo, Peter Li, and Allen Riddel.  Stan: A probabilistic programming language.  _Journal of Statistical Software_ , 76:1 –32, 2017.  10.18637/jss.v076.i01.
  * Cover and Thomas [2006] Thomas M. Cover and Joy A. Thomas.  _Elements of Information Theory_.  John Wiley & Sons, Inc, 2006.
  * Gelman et al. [2013] Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Ark Vehtari, and Donald B. Rubin.  _Bayesian Data Analysis_.  Chapman & Hall/CRC Texts in Statistical Science, 2013.
  * Giordano et al. [2018] Rian Giordano, Tamara Broderick, and Michael I. Jordan.  Covariances, robustness, and variational bayes.  _Journal of Machine Learning Research_ , 19:1 – 49, 2018\.
  * Gronau et al. [2020] Quantin F. Gronau, Henrik Singmann, and Eric-Jan Wagenmakers.  bridgesampling: An R package for estimating normalizing constants.  _Journal of Statistical Software_ , 92, 2020.
  * Jordan et al. [1999] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul.  An introduction to variational methods for graphical models.  _Machine Learning_ , 37:183 – 233, 1999.
  * Kingma and Welling [2013] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv:1312.6114_ , 2013.
  * Kucukelbir et al. [2017] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David Blei.  Automatic differentiation variational inference.  _Journal of machine learning research_ , 18:1 – 45, 2017\.
  * MacKay [2003] David J.C. MacKay.  _Information theory, inference, and learning algorithms_.  2003\.
  * Magnusson et al. [2022] M$\mathring{\text{a}}$ns Magnusson, Paul-Christian Bürkner, and Aki Vehtari.  posteriordb: a set of posteriors for bayesian inference and probabilistic programming, 2022.
  * Margossian and Mukherjee [2021] Charles C Margossian and Sumit Mukherjee.  Simulating ising and potts models at critical and cold temperatures using auxiliary gaussian variables.  _arXiv:2110.10801_ , 2021.
  * Margossian et al. [2020] Charles C. Margossian, Aki Vehtari, Daniel Simpson, and Raj Agrawal.  Hamiltonian Monte Carlo using an adjoint-differentiated Laplace approximation: Bayesian inference for latent gaussian models and beyond.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 33, October 2020.
  * Meng and Schilling [2002] Xl Meng and S Schilling.  Warp bridge sampling.  _Journal of Computational and Graphical Statistics_ , 11:552 – 586, 2002.  doi:10.1198/106186002457.
  * Mukherjee et al. [2018] Rajarshi Mukherjee, Sumit Mukherjee, and Ming Yuan.  Global testing against sparse alternatives under Ising models.  _Annals of Statistics_ , 46, 2018.  doi: 10.1214/17-AOS1612.
  * Parisi [1988] Georgio Parisi.  _Statistical Field Theory_.  Addison-Wesley, 1988.
  * Rasmussen and Williams [2006] Carl E. Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Turner and Sahani [2011] Richard E. Turner and Maneesh Sahani.  Two problems with variational expectation maximisation for time-series models.  In David Barber, A. Taylan Cemgil, and Silvia Chiappa, editors, _Bayesian Time series models_ , chapter 5, pages 109–130. Cambridge University Press, 2011.
  * Wainwright and Jordan [2008] Martin J. Wainwright and Michael I. Jordan.  _Foundations and Trends in Machine Learning_ , 1:1 – 305, 2008.
  * Wang and Titterington [2005] Bo Wang and Donald M. Titterington.  Inadequacy of interval estimates corresponding to variational bayesian approximations.  _Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics_ , R5:373 – 380, 2005.
  * Zhang et al. [2022] Lu Zhang, Bob Carpenter, Aki Vehtari, and Andrew Gelman.  Pathfinder: Parallel quasi-newton variational inference.  _Journal of Machine Learning Research_ , 23:1 – 49, 2022\.

## Appendix

In appendix A, we provide further details on FG-VI when the covariance matrix
has constant off-diagonal terms. Next appendix B provides an algorithm to
efficiently compute bounds on the shrinkage and delinkage terms, using
techniques developped in section 4. In appendix C, we show how to derive
bounds on the trace of the shrinkage matrix. In appendix D, we obtain an upper
bound on the KL divergence—equivalently the entropy gap—between $q$ and $p$.

##  Appendix A FG-VI for covariance with constant off-diagonal terms

Let $p$ be a multivariate Gaussian distribution over $\mathbb{R}^{n}$ with
covariance matrix $\boldsymbol{\Sigma}$, and let $q$ be the solution of FG-VI
with diagonal covariance matrix $\boldsymbol{\Psi}$, where

$\Psi_{ii}=\frac{1}{\Sigma^{-1}_{ii}}$ (39)

as in eq. (4). The elements of the correlation matrix are related to those of
the covariance matrix by

$C_{ij}=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}.$ (40)

In this section we assume that the correlation matrix has constant off-
diagonal terms, and we use $\varepsilon\in[0,1)$ to denote the value of these
terms. Note that we require $\varepsilon\geq 0$ since three or more random
variables cannot all be mutually anti-correlated. Also we require
$\varepsilon\\!<\\!1$ since otherwise $\mathbf{C}$ (and hence
$\boldsymbol{\Sigma}$) would not be positive-definite.

We now prove Theorem 3.5, broken up into a statement about the estimated
variance (Proposition A.1) and a statement about the entropy gap (Proposition
A.2).

######  Proposition A.1.

If the correlation matrix in eq. (40) has constant off-diagonal terms, then
the solution for FG-VI in eq. (39) obeys the following limits:

$\displaystyle\underset{n\to\infty}{\lim}\Psi_{ii}$ $\displaystyle=$
$\displaystyle(1\\!-\\!\varepsilon)\Sigma_{ii}$ (41)
$\displaystyle\underset{\varepsilon\to 1}{\lim}\Psi_{ii}$ $\displaystyle=$
$\displaystyle 0$ (42)
$\displaystyle\underset{n\to\infty}{\lim}\frac{1}{n}\mathrm{trace}({\bf S})$
$\displaystyle=$ $\displaystyle\frac{1}{1-\varepsilon}$ (43)

where $\varepsilon\in[0,1)$ denotes the value of $C_{ij}$ for $i\neq j$.

###### Proof.

Let $\mathbf{1}\in\mathbb{R}^{n}$ denote the vector of all ones. Then the
correlation matrix can be written as

$\mathbf{C}=(1-\varepsilon)\mathbf{I}+\varepsilon{\bf 1}{\bf 1}^{\top}.$ (44)

One can verify by direct substitution that the inverse correlation matrix has
elements

$\mathbf{C}^{-1}=\tfrac{1}{1-\varepsilon}\left[{\bf
I}-\tfrac{\varepsilon}{1+(n-1)\varepsilon}{\bf 1}{\bf 1}^{\top}\right].$ (45)

Recall from eq. (15) that $\Psi_{ii}=\Sigma_{ii}/C_{ii}^{-1}$. With some
algebra, it follows from eq. (45) that

$\Psi_{ii}=\left[\frac{(1-\varepsilon)(1+(n\\!-\\!1)\varepsilon)}{1+(n\\!-\\!2)\varepsilon}\right]\Sigma_{ii}.$
(46)

It is straightforward to take the limits of eq. (46) as $n\rightarrow\infty$
or $\varepsilon\rightarrow 1$, and these limits yield the results of the
proposition.

Finally, recalling that $\text{trace}({\bf
S})=\sum_{i=1}^{n}\Sigma_{ii}/\Psi_{ii}$, we immediately get

$\underset{n\to\infty}{\lim}\frac{1}{n}\text{trace}({\bf
S})=\frac{1}{1-\varepsilon}.$ (47)

∎

######  Proposition A.2.

If the correlation matrix in eq. (40) has constant off-diagonal terms, then
the per-component entropy gap from FG-VI vanishes in the limit
$n\rightarrow\infty$; that is,

$\underset{n\to\infty}{\lim}\
\frac{1}{n}\left[\mathcal{H}(p)-\mathcal{H}(q)\right]=0.$ (48)

###### Proof.

Recall from Theorem 3.2 of the main paper that the entropy gap for FG-VI is
given by

$\mathcal{H}(p)-\mathcal{H}(q)=\tfrac{1}{2}\log|\mathbf{S}|+\tfrac{1}{2}\log|\mathbf{C}|,$
(49)

where $\mathbf{S}$ is the diagonal shrinkage matrix with elements
$S_{ii}=\Sigma_{ii}/\Psi_{ii}$. We consider each term on the right side of
this equation in turn. It follows at once from eq. (46) that

$\log|\mathbf{S}|=n\left[\log\frac{1+(n\\!-\\!2)\varepsilon}{(1\\!-\\!\varepsilon)(1+(n\\!-\\!1)\varepsilon)}\right],$
(50)

where $\varepsilon\\!>\\!0$ denotes the amount of off-diagonal correlation.
Next we show how to evaluate $\log|\mathbf{C}|$. From eq. (44), we rewrite the
correlation matrix as

$\mathbf{C}=(1\\!-\\!\varepsilon){\bf
I}+n\varepsilon\left(\tfrac{1}{\sqrt{n}}{\bf
1}\right)\left(\tfrac{1}{\sqrt{n}}{\bf 1}\right)^{T}.$ (51)

Note that the second term on the right side of eq. (51) is a rank-one matrix
whose one nonzero eigenvalue is equal to $n\varepsilon$. By adding the first
term—which is a multiple of the identity matrix—we obtain a new matrix whose
eigenvalues are shifted by a uniform amount. It follows that this new matrix
(namely, $\mathbf{C}$) has $n\\!-\\!1$ eigenvalues at $1\\!-\\!\varepsilon$
and one eigenvalue at $1+(n\\!-\\!1)\varepsilon$, so that

$\log|\mathbf{C}|=(n\\!-\\!1)\log(1\\!-\\!\varepsilon)+\log(1+(n\\!-\\!1)\varepsilon).$
(52)

The entropy gap is related to the sum of $\log|\mathbf{S}|$ and
$\log|\mathbf{C}|$ by eq. (49). Adding the results in eq. (50) and eq. (52),
we find that

$\displaystyle\log|\mathbf{S}|+\log|\mathbf{C}|$
$\displaystyle=-\log(1\\!-\\!\varepsilon)+\log(1+(n\\!-\\!1)\varepsilon)$
$\displaystyle\mbox{\hskip 8.61108pt}\hskip 12.91663pt+\
n\log\left[\tfrac{1+(n-2)\varepsilon}{1+(n-1)\varepsilon}\right].$ (53)

Note that the first term on the right side is $\mathcal{O}(1)$, the second
term is $\mathcal{O}(\log n)$, and the third term can be written as

$n\log\left[\tfrac{1+(n-2)\varepsilon}{1+(n-1)\varepsilon}\right]=n\log\left[1-\tfrac{\varepsilon}{1+(n-1)\varepsilon}\right].$
(54)

For large $n$, the log term in this equation is $\mathcal{O}(\frac{1}{n})$ so
that the entire expression is $\mathcal{O}(1)$. From eq.(53), it therefore
follows that the entropy gap in eq. (49) is $\mathcal{O}(\log n)$. Dividing by
$n$, we see that the per-component entropy gap vanishes in the limit
$n\rightarrow\infty$, thereby completing the proof. ∎

From (53), we also see that the entropy gap becomes infinite as
$\varepsilon\to 1$ (for fixed $n$). Additional limits can be considered with
respect to both $\varepsilon$ and $n$, but we do not pursue those here.

##  Appendix B Solutions for the bounds on the shrinkage and delinkage terms

By exploiting the symmetries that we proved in Section 4.2, we can efficiently
compute bounds on the terms $\log|\mathbf{S}|$ and $\log|\mathbf{C}|$; these
are the terms that arise, respectively, from the effects of shrinkage and
delinkage.

###  B.1 Upper bound on $\log|{\bf S}|$

First we show how to compute the upper bound on $\log|\mathbf{S}|$. Recall
that to do so, we must solve the optimization problem

$\max_{\boldsymbol{\lambda}\in\Lambda_{R}}\sum_{i=1}^{n}\lambda_{i}^{-1}.$
(55)

From Lemma 4.2, we know that all the elements of the solution assume the edge
values of $\lambda_{1}$ or $\lambda_{n}$ save for at most one which we denote
$\lambda_{k}$. At a high level, we solve the optimization by exhaustively
computing the optimal solution for each candidate value of
$k\in\\{1,\ldots,n\\}$, then choosing the particular value of $k$ whose
solution maximizes the overall objective function.

It remains only to show how to compute the solution for a particular candidate
value of $k$. Recall the constraints that $\sum_{i=1}^{n}\lambda_{i}=n$ and
$\lambda_{1}=R\lambda_{n}$. It follows that

$\lambda_{k}=n-\left[(k-1)R+n-k\right]\lambda_{n}.$ (56)

Using the constraints to eliminate $\lambda_{1}$ and $\lambda_{k}$, we can
write the objective function entirely in terms of $\lambda_{n}$. In this way
we find

$\sum_{i=1}^{n}\lambda_{i}^{-1}=\frac{1}{n-\left[R(k-1)+n-k\right]\lambda_{n}}+\frac{(k-1)}{R\lambda_{n}}+\frac{n-k}{\lambda_{n}}.$
(57)

Crucially, we also need to enforce the boundary conditions
$\lambda_{n}\leq\lambda_{k}\leq\lambda_{1}$, or equivalently

$\frac{n}{Rk+n-k}\leq\lambda_{n}\leq\frac{n}{R(k-1)+n-k+1}.$ (58)

Note that the simplified objective in (57) for fixed $k$ is convex in
$\lambda_{n}$; hence the maximizer must lie at one of the boundary values in
(58). By computing the objective for each boundary value of $\lambda_{k}$, we
find the optimal solution for this candidate value of $k$. Finally, we obtain
the overall solution to eq. (55) by considering all $n$ candidate values of
$k$ and choosing the best one.

###  B.2 Upper Bound on $\log|{\bf C}|$

Next we show how to compute the upper bound on $\log|\mathbf{C}|$. Recall that
to do so, we must solve the optimization problem

$\max_{\boldsymbol{\lambda}\in\Lambda_{R}}\left[\sum_{i=1}^{n}\log\lambda_{i}\right].$
(59)

From Lemma 4.3, we know that all eigenvalues other than $\lambda_{1}$ and
$\lambda_{n}$ must have the same value; we denote this value by $\lambda_{k}$.
From the constraint $\sum_{i=1}^{n}\lambda_{i}=n$, it follows that

$\lambda_{k}=\frac{n-(1+R)\lambda_{n}}{n-2}.$ (60)

Again, using the constraints to eliminate $\lambda_{1}$ and $\lambda_{k}$, we
can write the objective function entirely in terms of $\lambda_{n}$. In this
way we find

$\sum_{i=1}^{n}\log\lambda_{i}=(n\\!-\\!2)\log\frac{n-(1+R)\lambda_{n}}{n-2}+\log
R\lambda_{n}+\log\lambda_{n}.$ (61)

This objective is concave in $\lambda_{n}$, so we can locate the maximum by
setting its derivative with respect to $\lambda_{n}$ equal to zero. Some
straightforward algebra shows that this derivative vanishes when

$\lambda_{n}=\frac{2}{1+R}.$ (62)

Finally we need to check that this solution does not violate the boundary
conditions of the problem; in particular, we require that
$\lambda_{n}\geq\lambda_{k}\geq R\lambda_{n}$, or equivalently that

$\frac{n}{1+R(n-1)}\leq\lambda_{n}\leq\frac{n}{n-1+R}.$ (63)

These conditions are always satisfied for $n\geq 3$. Hence we obtain an
analytical solution for the upper bound on $\log|\mathbf{C}|$. Finally, note
that while the solution for $\lambda_{n}$ does not depend on $n$, the
optimized objective function does depend on $n$ through eq. (61).

Algorithm 1 provides an implementation of the above-described method.

1 Input: $R,n$

2

3 Function  _ObjF(_$\lambda_{n}$ , $k$_)_:

4 return
$\left(n-k+\frac{k-1}{R}\right)\frac{1}{\lambda_{n}}+\frac{1}{n-[R(k-1)+n-k]\lambda_{n}}$

5

6

7

8for  _$k$ in $\\{2,\cdots,n-1\\}$_ do

9 $\lambda_{a}\leftarrow\frac{n}{Rk+n-k}$

10 $\lambda_{b}\leftarrow\frac{n}{R(k-1)+n-k+1}$

11
$F_{k}\leftarrow\text{max}(\text{{ObjF}}\;(\lambda_{a}),\text{{ObjF}}(\lambda_{b}))$

12 if ($k=1$) $F\leftarrow F_{k}$

13 else $F\leftarrow\text{max}(F,F_{k})$

14

15 end for

16$U_{s}\leftarrow n\log(F/n)$

17

18$\lambda_{n}\leftarrow\frac{2}{1+R}$

19
$U_{c}\leftarrow\log\frac{1}{\lambda_{n}}+\log\frac{1}{R\lambda_{n}}+(n\\!-\\!2)\log\frac{n-2}{n-(1+R)\lambda_{n}}$

20

21

22Return: $U_{s}$, $U_{c}$

Algorithm 1 Upper bounds on $\log|\mathbf{S}|$ and $\log|\mathbf{C}|^{-1}$

##  Appendix C Bounds on the average variance shrinkage

We can also derive bounds on the average shrinkage in componentwise variance
in terms of the problem dimensionality, $n$, and the condition number, $R$, of
the correlation matrix. The average in this case is performed over the
different components of $\mathbf{z}$. Recall that the shrinkage in each
componentwise variance is given by $S_{ii}=\Sigma_{ii}/\Psi_{ii}$. Hence we
can also express this bound in terms of the trace of the shrinkage matrix,
$\text{trace}({\bf S})$.

######  Proposition C.1.

Suppose that the correlation matrix $\mathbf{C}$ has condition number $R$.
Then the solution for FG-VI in section 2 satisfies

$\min_{\lambda\in\Lambda_{R}}\sum_{i=1}^{n}\lambda_{i}^{-1}\leq{\rm
trace}(\mathbf{S})\leq\max_{\lambda\in\Lambda_{R}}\sum_{i=1}^{n}\lambda_{i}^{-1},$
(64)

where $\Lambda_{R}$ is the set defined in section 4.

###### Proof.

We showed in the proof of Theorem 3.1 that

$S_{ii}=\frac{\Sigma_{ii}}{\Psi_{ii}}=C_{ii}^{-1}.$ (65)

It follows that ${\rm trace}(\mathbf{S})={\rm
trace}(\mathbf{C}^{-1})=\sum_{i}\lambda_{i}^{-1}$, where
$\lambda_{1},\ldots,\lambda_{n}$ are the eigenvalues of $\mathbf{C}$. The
bound then follows from the relaxtion from the set $\mathcal{C}_{R}$ to the
set $\Lambda_{R}$ in section 4. ∎

###  C.1 Lower bound on $\text{trace}({\bf S})$

The optimization implied by eq. (64) is convex, since both the set
$\Lambda_{R}$ and the objective function $\sum_{i}\lambda_{i}^{-1}$ are
convex. In fact, this bound can be evaluated in closed form by using similar
methods as in section 4.2.

######  Lemma C.2.

Let $\boldsymbol{\lambda}\in\Lambda_{R}$ be the solution that minimizes the
left side of eq. (64). Then $\lambda_{i}\\!=\\!\lambda_{j}$ whenever
$1\\!<\\!i\\!<\\!j\\!<\\!n$.

###### Proof.

This proof follows the same argument as the proof (by contradiction) for Lemma
4.3. Suppose there exists a solution with intermediate elements $\lambda_{i}$
and $\lambda_{j}$ that satisfy
$\lambda_{1}\\!\geq\\!\lambda_{i}\\!>\\!\lambda_{j}\\!\geq\\!\lambda_{n}$.
Consider the effect on this solution of a perturbation that adds some small
amount $\delta\\!>\\!0$ to $\lambda_{j}$ and subtracts the same amount from
$\lambda_{i}$. For sufficiently small $\delta$, this perturbation will not
leave the set $\Lambda_{R}$; however, it will diminish the separation of
$\lambda_{i}$ from $\lambda_{j}$. As a result the objective
$\sum_{k}(1/\lambda_{k})$ experiences a change

$g(\delta)=\frac{1}{\lambda_{i}-\delta}-\frac{1}{\lambda_{i}}+\frac{1}{\lambda_{j}+\delta}-\frac{1}{\lambda_{j}}.$
(66)

Evaluating the derivative, we find
$g^{\prime}(0)=\lambda_{i}^{-2}\\!-\\!\lambda_{j}^{-2}<0$, so that the
objective function is decreased for some $\delta>0$. As before this yields a
contradiction, because any solution must be maximal, and hence stationary
(i.e., $g^{\prime}(0)\\!=\\!0$), with respect to small perturbations. ∎

With the above lemma, the $n$-dimensional optimization over $\Lambda_{R}$ can
be reduced to a one-dimensional optimization that can be solved in closed
form. The methods are identical to those in the previous appendix.

First we rewrite the constraint, $\lambda_{n}\leq\lambda_{k}\leq
R\lambda_{n}$, as

$\frac{n}{R(n-1)+1}\leq\lambda_{n}\leq\frac{n}{R+n-1}.$ (67)

Since the minimization problem is convex, a minima can be found at a
stationary point of the objective function

$\sum_{i=1}^{n}\frac{1}{\lambda_{i}}=\frac{(n-2)^{2}}{n-(1+R)\lambda_{n}}+\frac{1}{\lambda_{n}}+\frac{1}{R\lambda_{n}},$
(68)

which now only depends on $\lambda_{n}$. Differentiating and setting to 0, we
obtain the root-finding problem,

$\left[R(n-2)^{2}-(1+R)^{2}\right]\lambda^{2}_{n}+2n(1+R)\lambda_{n}-n^{2}=0,$
(69)

which can be solved exactly. It remains to check whether the roots violate the
constraints in (67), and pick the non-offending root which maximizes the
objective in (68). If both roots violate the constraints then, by convexity of
the problem, the solution must lie at one of the boundary terms in (67).

###  C.2 Upper bound on trace(S)

A similar approach gives us an upper bound on $\text{trace}({\bf S})$. In
fact, we have already solved the problem of maximizing the right side of (64)
when upper-bounding $\log|{\bf S}|$. It remains to apply the same strategy.

##  Appendix D Tighter upper bound on entropy gap

In Proposition 4.1 we derived separate upper bounds on the individual terms
$\log|\mathbf{S}|$ and $\log|\mathbf{C}|$. One upper bound on the entropy gap
(or equivalently, on ${\rm KL}(q,p)$) is obtained simply by adding these
separate bounds. However, a tighter upper bound is obtained by replacing the
separate optimizations in Proposition 4.1 by a single joint optimization:

${\rm KL}(q,p)\ \leq\
\frac{1}{2}\max_{\boldsymbol{\lambda}\in\Lambda_{R}}\left[n\log\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}^{-1}+\sum_{i=1}^{n}\log\lambda_{i}\right].$
(70)

In this appendix we sketch how to solve this optimization and evaluate this
bound in closed form. The first step is to make the change of variables,

$\omega_{i}=\frac{\lambda_{i}^{-1}}{\sum_{j=1}^{n}\lambda_{j}^{-1}},$ (71)

and to translate the domain of optimization accordingly. Under this change of
variables, the original domain $\Lambda_{R}$ in section 4 is mapped onto the
set

$\Omega_{R}=\left\\{\boldsymbol{\omega}\in\mathbb{R}_{+}^{n}\,|\,\omega_{n}\geq\ldots\geq\omega_{1}=\frac{1}{R}\omega_{n},\sum_{i=1}^{n}\omega_{i}=1\right\\}.$
(72)

Likewise, a little algebra shows that the optimization in eq. (70) is
equivalent to the following:

${\rm KL}(q,p)\ \leq\
\frac{1}{2}\max_{\boldsymbol{\omega}\in\Omega_{R}}\left[\sum_{i=1}^{n}\log\frac{1}{\omega_{i}}-n\log
n\right].$ (73)

Now we can make a similar argument as in the proof of Lemma 4.2 to simplify
this optimization.

######  Lemma D.1.

Let $\boldsymbol{\omega}\in\Omega_{R}$ be the solution that maximizes the
right side of eq. (73). Then at most one $\omega_{i}$ is not equal to either
$\omega_{1}$ or $\omega_{n}$.

###### Proof.

We prove the lemma by contradiction. Suppose there exists a solution with
intermediate elements $\omega_{i}$ and $\omega_{j}$ that satisfy
$\omega_{n}\\!>\\!\omega_{i}\\!>\\!\omega_{j}\\!>\\!\omega_{1}$. Consider the
effect on this solution of a perturbation that adds some small amount
$\delta\\!>\\!0$ to $\omega_{i}$ and subtracts the same amount from
$\omega_{j}$. Note that for sufficiently small $\delta$, this perturbation
will not leave the set $\Omega_{R}$; however, it will expand the separation of
$\omega_{i}$ from $\omega_{j}$. As a result the objective in eq. (73) changes
by an amount

$f(\delta)=\frac{1}{2}\left[\log\frac{1}{\omega_{i}\\!+\\!\delta}-\log\frac{1}{\omega_{i}}+\log\frac{1}{\omega_{j}\\!-\\!\delta}-\log\frac{1}{\omega_{j}}\right].$
(74)

Next we evaluate the derivative $f^{\prime}(\delta)$ at $\delta=0$; doing so
we find $f^{\prime}(0)=\omega_{j}^{-1}\\!-\\!\omega_{i}^{-1}>0$, so that the
objective is increased for some $\delta>0$. But this yields a contradiction,
because any solution must be maximal, and hence stationary (i.e.,
$f^{\prime}(0)\\!=\\!0$), with respect to small perturbations. ∎

With the above lemma, we can reduce the $n$-dimensional optimization over
$\Omega_{R}$ to a one-dimensional optimization that can be solved in closed
form; the methods are identical to those in the previous appendix.

In details, let $\omega_{k}$ be the one variable which (potentially) does not
go to $\omega_{1}$ or $\omega_{n}$. Given $\sum_{i=1}^{n}\omega_{i}=1$,

$\omega_{k}=1-(k-1+R(n-k))\omega_{1}.$ (75)

The objective is then

$\displaystyle\sum_{i=1}^{n}\log\frac{1}{\omega_{i}}=$
$\displaystyle-(k-1)\log\omega_{1}-(n-k)\log R\omega_{1}$ (76)
$\displaystyle-\log(1-[k-1+R(n-k)\omega_{1}]).\ $

Since we are trying to maximize a convex function, the solution does not lie
at a stationary point, rather at a boundary set by the constraint,
$\omega_{1}\leq\omega_{k}\leq\omega_{n}$, or equivalently

$\frac{1}{k-1+R(n-k+1}\leq\omega_{1}\leq\frac{1}{k+R(n-k)}.$ (77)

It remains to test each candidate boundary for each choice of $k$ to obtain a
maximizer.