Provable convergence guarantees for black-box variational inference
Authors: Justin Domke
Robert M. Gower
Guillaume Garrigos
Venue: NeurIPS 2023 poster (2023)
Keywords: optimization
variational inference
TL;DR: Black-box variational inference with projected or proximal SGD converges at a 1/√T rate for log-smooth and log-concave targets, and 1/T if strongly log-concave.

Abstract: Black-box variational inference is widely used in situations where there is no proof that its stochastic optimization succeeds. We suggest this is due to a theoretical gap in existing stochastic optimization proofs—namely the challenge of gradient estimators with unusual noise bounds, and a composite non-smooth objective. For dense Gaussian variational families, we observe that existing gradient estimators based on reparameterization satisfy a quadratic noise bound and give novel convergence guarantees for proximal and projected stochastic gradient descent using this bound. This provides rigorous guarantees that methods similar to those used in practice converge on realistic inference problems.

Provable convergence guarantees for black-box variational inference
Justin Domke
Guillaume Garrigos
University of Massachusetts Amherst domke@cs. umass.
Université Paris Cité and Sorbonne Université, CNRS Laboratoire de Probabilités, Statistique et Modélisation
F-75013 Paris, France
petriggs@@psm.par..ee
Robert Gower
Center for Computational Mathematics
Flatiron Institute, New York
igower@flatironistiii org
Abstract
Black-box variational inference is widely used in situations where there is no proof that its stochastic optimization succeeds. We suggest this due to a theoretical gap in existing stochastic optimization proofs-named the challenge of gradient estimators with unusual noise bounds, and composite non-smooth objective. For dense Gaussian variational families, observe that existing gradient estimators based on reparameterization satisfy a quadratic noise bound and give novel con- vergence guarantees for proximal and projected stochastic gradient descent using this bound. This provides rigorous guarantees that methods similar to those used practice converge realistic inference problems.
Introduction
Variational inference tries approximate complex target distribution with a distribution from simpler family [1, 2, 3, 4]. If p(z, æ) is a target distribution with latent variables E g and observed data E _ and lw is variational family with parameters W, the goal minimize the negative evidence lower bound (ELBO)
((w) =
E logp(z, x:)
+ log qm (z).
2~-q
zruqun
((1))
k(w)
For subsequent discussion, have decomposed the objective into the free energy and the negative entropy h. Minimizing of is equivalent to minimizing the KL-divergence between fu to p(-[z), because KL(g.. [|p((==)) fww + log p(x). Recent research has focused "black box" variational inference, where the target distribution p is sufficiently complex that one can only access through evaluating probabilities (or gradients) at chosen points [5, 6, 7, 8, 9, 10, 11]. Crucially, one can still get stochastic estimates of the variational objective / and of its gradient, and use these to optimize.
Still, variational inference can sometimes be unreliable [12, 13, 14, 15], and some basic questions remain unanswered. Most notably: does stochastic optimization of f converge minimum of the objective? There has been various progress towards answering this question. One line of research seeks to determine if the variational objective has favorable structural properties, such as smoothness or (strong) convexity [13, 16, 17] (Sec. 2.1). Another line seeks to control the "noise" (variance or expected squared norm) of different gradient estimators [18, 19, 20, 21] (Sec. 2.2). However, few
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

full convergence guarantees are known. That is, there are few known cases where applying given stochastic optimization algorithm given target distribution known to converge given rate.
We identify two fundamental barriers preventing from analysing this VI problem as a standard stochastic optimization problem. First, the gradient noise depends the parameters in non- standard way (Sec. 2.3). This adds great technical complexity and renders many traditional stochastic optimization proofs inapplicable. Second, stochastic optimization theory typically requires that the (exact) objective function is Lipschitz continuous or Lipschitz smooth. But in our VI setting, under some fairly benign assumptions, the ELBO is neither Lipschitz continuous nor Lipschitz smooth.
We obtain non-asymptotic convergence guarantees for this problem, under simple assumptions.
Central contributions (Informal). Suppose that the target model log p(., ) is concave and Lipschitz-smooth, and that gui is in a Gaussian variational family parameterized by the mean and factor of the covariance matrix (Eq. 2). Consider minimizing the negative ELBO using either one of the two following algorithms:
a proximal stochastic gradient method, with the proximal step applied to h and the gradient step applied to I, estimating Vi(w) with standard reparameterization gradient estimator (Eq. 5), and using a triangular covariance factor;
a projected stochastic gradient method, with the gradient applied to f 1+ h, estimating V (w)) using either of two common gradient estimators (Eq. 7 or Eq. 9), with the
projection done over symmetric and non-degenerate (Eq. 3) covariance factors
Then, both algorithms converge with 1/vT complexity rate (Cor. 12), or 1/T if we further assume that log p(,,)) is strongly concave (Cor. 13).
We also give new bound on the noise of the "sticking the landing" gradient estimator, which leads to faster convergence when the target distribution p is closer to Gaussian, up to exponentially fast convergence when is exactly Gaussian (Cor. 14).
This is achieved through series of steps, that we summarize below.
1. We analyze the structural properties of the problem. Existing results show that with a Gaussian variational family, if - log x)) is (strongly) convex Lipschitz smooth, then so is the free energy This is for instance known to be the case for some generalized linear models, and we give proof of convexity and smoothness for some hierarchical models including hierarchical logistic regression (see Appendix 7.3). The remaining component of the ELBO, the neg- neg-entropy h, is convex when restricted to appropriate set. It is not smooth, but it was recently proved to be smooth over a certain non-degeneracy set.
2. We study the noise of three common gradient estimators. They do not satisfy usual noise bounds, but we show that they all satisfy new quadratic bound (Definition 5). For the sticking- the- landing estimator, our bound formalizes the longstanding intuition that should have lower noise when the variational approximation is strong (Thm. 4).
3. We identify and solve the key optimization challenges posed by the above issues via new convergence results for the proximal and projected stochastic gradient methods, when applied to objectives that are smooth (but not uniformly smooth) and with gradient estimators satisfying quadratic bound.
1.1 Related work
Recently, Xu and Campbell [22] analyzed projected- SSG (stochastic gradient descent) for Gaussian VI using the gradient estimator we will later call gent (7), with a particular rescaling. They show that, asymptotically in the number of observed data, their method converges locally with a rate of 1 1/VT, under mild assumptions. Our results are less general in requiring convexity, but apply with other gradient estimators, and give a faster 1 1/T rate with strong convexity.
asympotit,,,,
Lambert et al. [23] introduce VI-like SGD algorithm, derived from discretisation of a Wasserstein gradient flow, and show converges at a 1/" rate for the 2- 2-Wasserstein distance when the log posterior is smooth and strongly concave. This line continued by Diao et al. [24], who propose

---

proximal-SGi method based on the decomposition f f= 1 h. They obtain a 1/VT rate when the log posterior is smooth and convex, and a 1/T rate when it is smooth and strongly convex. Unlike typical black-box VI algorithms used in practice, these algorithms require computing the Hessian of the posterior. We analyze more straighforward applications of SGD to the VI problem using standard gradient estimators. Under the same assumptions, our algorithms have the same rates for KL-divergence, which imply the same rates for 2- 2-Wasserstein distance by of Pinsker's s inequality and that the total-variation norm upper-bounds Wasserstein distance [25, Remark 8.2].
In concurrent work, Kim al. [26], consider proximal -SD method similar to our approach in Sec. 6. They obtain a 1/ll convergence rate, similar to what Cor. 12 gives when the log posterior strongly concave. They also consider alternative parametrizations of the scale parameters that render the ELBO globally smooth, and obtain nonconvex result: Under a relaxed version of the inequality, they can guarantee 1/7 rate.
2 Properties of Variational Inference (VI) problems
Traditionally, the ELBO was optimized using message passing methods, which essentially assume that p and q are simple enough that coorcccoee updates are possible [2, 3, 4]. However, in the last decade a series of papers developed algorithms based on a "black-box" model where is assumed to be complex enough that one can only evaluate log p (or its gradient) at selected points z. The key observation is that even if is quite complex, it is still possible to compute stochastic estimates of the gradient of the BLBOO which can be deployed in stochastic optimization algorithm [5, 6, 8, 9, 10].
This paper seeks rigorous convergence guarantees for this problem. We study the setting where the variational family is the set of (dense) multivariate Gaussian distributions, parameterized in terms of w = (m, C), where m is the mean and C' is a factor of the covariance, i.e.
qu (z) = N (z/m, CC").
We optimize over E W, where W = [(m, C) C' > 0} and will further require C' to be either symmetric or triangular. This does not really change the problem, since for a given covariance matrix Σ there always exists symmetric (or lower triangular) positive definite factor C' such that CCT E.
Now, is possible to solve this optimization problem? Without further assumptions, it is unlikely any guarantee is possible, because it is easy to encode NP-hard problems into this VI framework [27, 28). We discuss below the assumptions we will make to be able to solve the problem.
2.1 Structural properties and assumptions
The properties of the free energy / depend on the properties of the target distribution p. necessary to assume that p is somehow "nice" to ensure the problem can be solved. Titsias and Greeilla [17, Proposition 1] showed that if - log p(., xx) is convex, then / is convex too. Challis and Barber [16, Sec 3.2] showed that if the likelihood p(-[z) is convex and the prior p(z) is Gaussian, then / is strongly-convex. Domke [13, Theorem 9] showed that if - logp(', x)) is ju-strongly convex, then 7 is j-strongly convex as well, and that the constant is sharp. Similarly, Domke [13, Theorem 1] showed that if log p(.,)) 2) is M smooth, then / also -smoohhh, and that the constant is sharp.
In this paper we make two assumptions about the target distribution p: the negative fog-probabiliiyy - -log(::;>) must be convex (or strongly convex), and Lipschitz smooth. Section 7.3 (Appendix) gives some example models where these assumptions are satisfied. For example, if the model is Bayesian linear regression, or logistic regression with a Gaussian prior, then - log æ) is smooth and strongly convex. In addition, if the target is hierarchical logistic regression model, then - llo p(o x) is smooth and convex.
(mm.C"" C*) of (w) can look like. Intuitively, if the target log p(., x) is justrrongly concave, then we would expect the target distribution to be "peaky". This means that the optimal distribution would be close to a delta function centered at the MAP
Assumptions on also impact what a minimizer w' = (m",
solution: mn* would be close to some maximum would not be too large. This intuition can be formalized: in this context we have |C*1|2 + I|m mi S d/p [13, Theorem 10). Similarly, if log p(,, æ)) is M-Lipschitz smooth, then we expect that the target is not too concentrated, so we might expect that the optimal covariance cannot not be too
((.p((,))) ar) noted m, and the covariance factor C*
3

---

small. Formally, it can be shown that the singular values of the covariance factor C* are greater than 1 1/VM (13, Theorem 7].
The properties of the neg-entropy h are inherited from the choice of the variational family and do not depend on p. Since we consider a Gaussian variational family, h is known in closed-form. To avoid some technical issues, h(w) is equal to - In det C' up to a constant if C is positive definite (see Appendix 7.1), and to too otherwise So inherits the properties of the negative log determinant, meaning is proper closed convex function whenever C' is symmetric or triangular (see Appendix 7.2). From an optimization perspective, it is natural use a gradient-based algorithm. But unfortunately h is not Lipschitz smooth on its domain, because its gradient can change arbitrarily quickly when the singular values of C' become small. However, can be shown [13, Lemma 12] that is M-smooth over the following set of non-singular parameters (which contains the solution w", see the previous paragraph) given by
will restrict / so that the covariance factor is positive definite, so
WM w (m, C): :CC > 0 and miin (C) minn
(3)
VM
Instead of computing gradients, an optimizer might want to use the proximal operator of h, which also be computed in closed form (see Sec. 4).
2.2 Gradient estimators
This paper considers gradient estimators based on the "path" method or "reparameterizationn. These assume some base distribution is known, together with some deterministic transformation Tw such that the distribution of Tw (u) is equal to Iw when u ~ s. Then, we can write
V E Φ(z) 2roque
E V (U). urve
(4)
In the case of multivariate Gaussians, the most common choice is to use -N(,,) and (u,
Cu + m, which we will consider in the rest of the paper.
We will consider three different gradient estimators based on the path-type strategy (Eq. 4). We will provide new noise bounds for each of them (all the proofs are Appendix 8). Our analysis is mostly based on the following general result [20, Thm. 3].
Theorem 1. Let (u) = Cu + foror w = (m, C). Let do IRd - R be M-smooth, suppose that ø is stationary m, and define
(m, 0). Then
(T= S (d + 1)M² || pmm-mli- mll (d+3)M2 | s (d (d+3) || wl12 (
E
Furthermore, the first inequality cannot improved.
Intuitively, this result says that the noise of gradient estimator is lower when the scale parameter C is small and the mean parameter m close to a stationary point.
The first estimator we consider is Vi(w) only. It is given by taking Φ log p into Eq. 4, i.e.
geuergy (u) = V log (o)) u ~~,
(5)
The next result gives bound the noise of this estimator, which is direct consequence of Theorem The second line uses Young's inequality to bound the noise in terms of (i) the distance of w from w" and (ii) a constant determined by the distance of w" from some fixed parameters w.
Theorem 2. Suppose that log x)) is M M-smooth and has maximum (or stationary point) at m. and define w (m, 0). Then, for every w and every solution w of the VI problem.
E genergy (u)l
<
(d 3)M² || u1/,
(6)
2(d + 3)M w= 12 2(d + 3)M
Second, we consider an estimator of the gradient of the full objective / ++ h. It obtained by simply taking the above estimator and adding the true (known) gradient of the neg-entropy, i.e.
Bent (u) := genergy (u) + Vh(w), u -~ u~

---

The noise of geut can be bounded since it only differs from genergy by the deterministic quantity Vh(w), and the fact that provided w E Wi,-the singular values of w cannot be too small and so Vh(w) cannot be too large. This is formalized in the following theorem, where again, the second line relaxes the result into term based on the distance of from w² plus constants. (Another type of noise bound [29] for gent was obtained by Kim et al. [21] for diagonal Gaussians and variety of parameterizations the covariance.)
Theorem 3. Suppose that log (()))) is M-smooth, that is maximal at m, and define w (m² 0).
Then, for every L > for every w E WE and every solution wt of the VI problem,
E ||gett (u)l2
2(d 3) || 111 + 2d..
(8)
4(d + 3) M² |/ w" 112 + 4(d 3)M ww? 2dT.
While gent used the exact gradient of h. it may be beneficial to use a stochastic estimator / instead. To derive our third estimator, write the gradient of the ELBO as'
Vi(w)) + Vh(w)= V w
E zm
log + log qu (z)] "=="
where the parameters u serve as a way to "hold w constant under differentiation". This leads to our third estimator, called the "sticking the landing" (STL) gradient estimator
(u) Benergy (u) + V. log qu (Tu
W(0,
(9)
Intuitively, we expect that ESTI. will tend to have lower noise than gent when the posterior is well- approximated by the variational distribution. The reason is that, as observed by Roeder al. [30], if (z) were a perfect approximation of p(z[z), then the two terms in Eq. 9 would exactly cancel (for every u) and so the estimator would have zero variance. Below we formalize this intuition in what we believe is the first noise bound on gsTL.
Theorem 4. Suppose that log log x) is M-smooth. Consider the residual r(2) := log x) log 9m. (2) for any solution w of the VI problem, assume that it has stationary point mn, and define v = (m, 0). Then r is K Ksmooohh for some E [0, 2M), and for all w E WM.
E ESTTL
I 8(d 8(d+3)MP w" w01 2(d + 3) m| -0ll
(10)
< A(d+3)( 2M2 || w 112 + 4(d + 32 Iw wl
Moreover, if p(12)) is Gaussian then K 0.
When the target is Gaussian then K = 0 in Eq. 10, meaning that when w = w., the STL estimator has no variance. This is distinguishing feature of gsTL, opposed to genergy or gent-
2.3 Challenges for optimization
Three major issues bar applying existing analysis of stochastic optimization methods to our VI setting: 1) non-smooth composite objective 2) lack of uniform smoothness and 3) lack of uniformly bounded noise of the gradient estimator.
The first issue is due to the non con-smothhness of neg-entropy function h. This means that under the benign assumption that the target log p is smooth, the full objective + cannot be smooth, since a nonsmooth function plus a smooth function is always nonsmooth. This renders stochastic optimization proofs (e.g. those for the "ABC" conditions [29]) that do not use projections or proximal operators inapplicable.
One way to tackle this first issue would be use a non-smooth proof technique, but these rely on having a uniform gradient noise bound (our third issue). Alternatively, one can overcome the non-smoothness of the neg-entropy function by either using proximal operator, or projecting onto a set where h is smooth, namely WM. This is the strategy we pursue.
'One could also build estimator without holding constant in this way. However, in the case of Gaussian variational families, this turns out to be mathematically equivalent to using closed form entropy because log q (u)) = = ||ul| log(2 log |C], which has the same gradient (independent of u) as h(w)..
5

---

Table 1: Table of known black-box VI properties relevant to optimization. In some cases there could be multiple valid w" or which case these results hold for all simultaneously.
Description
Definition
Estimator for VI
Senergy V
Estimator for V + Vh
gent - V. logp(Tic + Vh(w)
Estimator for VI + Vh
BSTL Ww + V log (Tu
Constraint set
WM {(m, C) C 0, min (C) 2 1/vM)
Optimum
w" E argmin (ww) + huw)
Stationary point of
l = (z, 0) for any that is stationary point of log p((p))
Condition on log 2(z
Consequence
none
h(w) is convex when is symmetric or triangular
k(w) is M-smooth over WM
convex
((w) is convex
j--strongly convex
(ww) is j-strongly convex
||0- w 1 - /M
M-smooth
((w) is M-smooth
w² E War
Benergy gent, and ESTL are quadratically bounded
The second issue is that existing analyses for proximal/projectedd stochastic methods either rely on a uniform noise bound ([31, Cor. 3.6)) or uniform smoothness [32, 33, 34, 35, 36), both of which not known to be true for VI. By uniform smoothness, we refer to the assumption that p(T. (w), can only guarantee that log (u), x) is smooth in expectation, i.e. that ((w) is smooth. Several works [9, Cond. 1][37, Thm. 1][38, Sec. 4][39, Assumption A A1][40, Thm. 1][41, Assumption 3.2] assumed that the full objective / + h is smooth, but we are not aware of cases where this holds in practice for VI and in our parameterization (Eq. this cannot be true log p smooth.
M-smooth for every u, with M being independent of u. Instead,
The third issue is the lack a uniform noise bound. Most non-smooth convergence guarantees within stochastic optimization [42] assume that the noise of the gradient estimator is uniformly bounded by constant. But this does not appear to be true even under favorable assumptions-e.g. it is untrue if the target distribution is standard Gaussian. The best that can hope is that the gradient noise can be bounded by a quadratic that depends on the current parameters the estimator-even this may only be true when the parameters optimization contribution is to provide theoretical guarantees for stochastic algorithms under such specific noise regime, which we make precise in the next definition.
and-depending on in the set War- Thus, our main
Definition 5. Let & be differentiable function. We say that a random vector g is a quadratically bounded estimator for VΦ w with parameters (a, b, w²), it is unbiased E ig] Vo(w) and the expected squared norm is bounded by quadratic function of the distance of parameters w to w", iee. E | w w'1l b.
The noise bounds derived in Section 2.2 for the estimators Benergy- Bents and ESTL imply that they are uniformly quadratically bounded estimators. See Appendix 8.2 for a table with the corresponding constants and b.
3 Stochastic Optimization with quadratically bounded estimators
In this section we give new convergence guarantees for the Prox-SGD algorithm and the Proj-SGD algorithm with quadratically bounded gradident estimators. Because these may be of independent

---

interest, they are presented generically, without any reference to the VI setting. We specialize these results VI in Section 4. For both algorithms, we present results assuming the problem be strongly convex or just convex. All proofs for this section can be found in the Appendix (Sec. 9).
3.1 Stochastic Proximal Gradient Descent
Here we want to minimize a function which the sum of two terms / + h, were both and are proper closed convex functions, and / is smooth. For this we will use stochastic proximal gradient.
Definition 6 (Prox-SGD). Let w be fixed initial parameters and let 71 72.- be sequence of step sizes. The stochastic proximal gradient (Prox-SGD) method given by
w y (w² - -xgg")
where g¹ is gradient estimator for Vl(w"), and the proximal operator is defined
1 2y 1/ v112.
pwo = argmin h(v) +
Theorem 7. Let 1 be p-strrongy convex and M-smooth function, and let i = argmin(l). Ler h be proper closed convex function, and let w = argmin(l + h). Let (w")en be generated by the Prox-SGI algorithm, with a constant stepsize y E (0, min Za 41] Suppose that g is quadratically bounded estimator (Def. 5) for VI with parameters (a, b, wt) Then,
E w T+1 w² 10 (1 "1))" | -=0=="
w011
27
(b M ||0" - m11)) -w0l)
(11)
28 +1
Alternatively, we the decaying stepsize "/t min
then
2a μ (t 1)2
Bllw" - 11? 16|a 16|a/4212
1/000 - w0 112 + 22T
(b + M² | - -m1l))
(12)
T2
In both cases, T !) iterations sufficient to guarantee that E||w w011
112 S
The above theorem gives an anytime 1/l rate of convergence when the stepsizes are chosen based on the strong convexity and gradient noise constants μ and a. Note that do not need to know precisely those constants: We show in Appendix 9.3 that using any constant step size proportional to T log I' leads a log TT rate.
Theorem 8. Let I. be proper convex and M-smooth function. Let h be proper closed convex function, and let w" E argminn( + h). Let (w" teen be generated by the Prox-SGD algorithm, with constant stepsize 7 E (0, ). Suppose that ge a quadratically bounded estimator (Def. for V with parameters (a, b, w*). Then,
ww10 -w*1l2 E [f(w²) inf 7
7(a
(1 HT) + +b).
1 20y? and T del Σt @c+hw²
where A del
Σt a++1 In particular, if
VaT then
_22 2
E [[(W))] - inf < VaT
(2allw² - w w01 +b) VT > max
Thus, = O(E 2) iterations are sufficient to guarantee that EE [f(w²) inf f] Se
3.2 Stochastic Projected Gradient Descent
Here we want to minimize a function f over set of constraints W. where W is a nonempty closed convex set and f is a proper closed convex function which differentiable on W. To solve this problem, we will consider the stochastic projected gradient algorithm.
7

---

Definition 9 (Proj-SGD), Let w" be some fixed initial parameter. and let 71, "2, be a sequence of step sizes. The projected stochastic gradient (Proj-SGD) method is given by
w ++ projw (w' ->gg).
where g² is gradient estimator for ff("'), and the projection operator is defined as
projw(w) argmin ||0 11lo-ll2.-
w12.
VEW
Note that we do not require f to be smooth, meaning that this setting is not particular case of the one considered in Section 3.1. In fact, the arguments we use in the proofs for Proj-SGD are more closely related to stochastic subgradient methods than stochastic gradient methods.
Theorem 10. Let W be a nonempty closed convex set. Let f be u-strongly convex function,
differentiable on W. Let w" argminn (f). Let (w")Le be generated by the Proj-SGD algorithm,
with constant stepsize 7 (0, min 20 ]]]. Suppose that g is a quadratically bounded estimator (Def. 5) for f with parameters (a, b, w²) Then,
E B||w²
w' w'll (1
p9) 2)
1| 100'-w"l'
27b
(13)
P
2 28 +
Alternatively, if we use the decaying stepsize "/ min
then
2a
11 (t 1)3
32a
16b
E [ w" [|3]
mml
w²1 +
(14)
2
12T
In both cases, T = O(c iterations are sufficient to guarantee that B||wT - wn Sc.
Note that Eq. 13 is sum of two terms: one that decays exponentially in T and one that decreases only when the stepsize 7 is sufficiently small. This has an important consequence: if one uses the gradient estimator gsTL and the target distribution exactly a Gaussian, then b = (Appendix 8 8.2), meaning that the algorithm will converge at an exponential rate. This is similar to many results in the stochastic optimization literature showing faster rates hold when interpolation holds [35].
Theorem 11. Let W be a nonempty closed convex Let be a convex function, differentiable on W. Let w" E ( Let (w*) )ee be generated by the Proj-SGD algorithm, with a constant stepsize 7 E (0, +o0) Suppose that g" is quadratically bounded estimator (Def. 5) for V ʃ at w with constant parameters (a, b, w*). Then,
w0
[(
inf F W
S 2
a
+
0r
1+ a and ut def ET-1 0++1 w²
where 0 def
Finally if
/2 Vat
and T 2 then
ΣT- g++1
2a
E [f(("))] inf inf
Iw00 - w²1
(15)
V2aT
Thus, T = 2) iterations are sufficient to guarantee that TE [[(=²)- inf s] Se.
4 Solving VI with provable guarantees
We now specialize the optimization results of the previous section to our VI setting. We aim to minimize the negative ELBO, f REd x V - IR U [+oo}, which decomposes f = h where
((w) = logp(z, x) and h(w) =
log qu (z) if C > 0, +oo otherwise,
(16)
zr
zrogu
and where V is the vector space of matrices in which the covariance factors C' belong. Under the assumption that the log-target p(.,)))) M M-smooth and concave, we propose two strategies to minimize this objective. All the proofs for this section can be found in Appendix 10.

---

1. Apply the Prox-SGD algorithm to the sum ++ h, using proximal step with respect to h, and stochastic gradient step with respect to L. The gradient of will be estimated with Benergy (5), which is globally quadratically bounded. The prox of h admits a closed form formula if we choose that V is the space of lower triangular matrices.
2. Apply the Proj-SGD algorithm to minimize ʃ over Was. The gradient of f will be estimated using either gent (7) or ESTL (9), both of which are quadratically bounded on WM. The projection onto WM admits closed form formula if choose that the space of symmetric matrices.
Corollary 12 (Prox-SGD for VI). Consider the VI problem where qu a multivariate Gaussian distribution (Eq. 2) with parameters
(m, C) R2d x T and assume that this problem admits a solution w". Suppose that log p(,,,) ) is M-smooth and concave (resp. ju-strongly concave). Generate sequence w' by using the Prox-SGI algorithm (Def. 6) applied to I and h (Eq. 16), using genergy (5) as an estimator of Vl. Let the stepsizes be constant and equal to 1/( denergy')) (resp. be decaying in Theorem 7 with Genergy = 2(d + 3) Then, for cerlain average w² of the iterates, we have for 2 2 thht
E [s(w")) inf s] = 0(1/V)) (resp. E [lIW" - w = O(1TT).
Corollary 13 (Proj-SGI for VI). Consider the VI problem where is a multivariate Gaussian distribution (Eq. with parameters w (m, C) E Rd Sa, and assume that this problem admits solution w². Suppose that log p(:, x) is M-smooth and concave (resp. pu-strongly concave). Generate a sequence w by using the Proj-SGD algorithm (Def. 9) applied to the function f 1+hh (Eq. 16) and the constraint Was (Eq. 3), using gent (7) or gSTL (9) as an estimator of V f. Let the stepsizes "it be constant and equal to / /2/(aT)) (resp. be decaying as in Theorem 10) with Gent = 4(d + 3) M2
or a = aSTT = 24(d 33 Then, for certain average w" of the iterates, we have T 2 2 that
E [s(w") inf s] O(1/VT) (resp. E [llw"-w""
O(ITT)-
Corollary 14 (Proj-SGD for VI Gaussian target). Consider the setting of Corollary 13, in the scenario that log p(,,)) is jesttronlyy concave, that we use the gSTL estimator, and that we take a constant stepsize Tt = E (0, min 2AgTL [33. Assume further that p(-1x) is Gaussian. Then,
E [|lWT w 12 < (1 17) T 1| w w'11
Let us now discuss the practical implementation of these two algorithms (more details and an explicit implementation of the methods are given in Appendix 10.2). These require computing either the proximal operator of the neg-entropy h, or the projection onto the set of non-degenerate covariance factors WM- These be computed [13] for every (m, C) as:
prox, (w) = (m, C + AC), where AC' is diagonal with AC, = 1 + Cis),
(w) DDUT, where C has the SVD C UDUT and D.s max(maa 1/ M).
Our theory for the Prox-SGD algorithm formally assumes that the covariance factor C lives in the vector space of lower-triangular matrices. This can be implemented by letting C Raxd and "clamping" the apper-triangulat entries to zero throughout computation. Then the gradient genergy (5) can be computed using automatic differentiation, and the proximal operator can be computed as the above equation. Or, one may choose w (m, c) where E [dx(d+1)/2 is the lower-triangular entries of the matrix so C = tril(c). Gradients with respect to w = (m, c) can again be estimated using automatic differentiation (including the tril operator), and the proximal operator can be computed by forming C, using the above formula, and then extracting the lower- triangular entries.
Similarly, theory for the Proj-SGD algorithm formally assumes the covariance factor C lives in the vector space of symmetric matrices. This can be implemented by letting E gdxd com- puting the gradient estimator gent (7) or ESTL (9) over the reals using automatic differentiation, and "symmetrizing" " C' throughout computation. That is, set C' + + + after applying each gradient update, right before projection. Or, one may choose = (m, c) where c E (dx(141)/2 is the trian-ur entries of the matrix, so C' symm(e). where symm is function similar to tril except creating symmetric matrices. The gradient estimator gent (7) or gsTL (9) can be computed

---

using automatic differentiation (including the symm operator). In this case, to project one would first form C. then use the above formula, and then extract the lower- trianntri matrices.
In terms of cost, computing log (z) or (wo(t) needs ((d*) operations. Computing the prox of requires d operations, but in contrast projecting onto WM requires diagonalizing a symmetric matrix, which takes ((d*) operations. We note that this is not necessarily an issue, because 1) eigenvalue decomposition has excellent computational constants in moderate dimensions; 2) computing the target log p( 2, .) may itself require (dd)) or more operations; 3) it is common to average a "minibatch" of gradient estimates, meaning each eigenvalue decomposition computation is amortized over many ELBO evaluations. Still, in high dimensions, when log p(z, æ) is inexpensive, and small minibatches are used, computing such decomposition in each iteration could become a computational bottleneck. In this scenario, the Proj-SGD algorithm is clearly cheaper, with its ((d) cost for the proximal step.
5 Discussion
While this paper has focused on convergence with dense Gaussian variational distributions, the results in Sec. 3 apply more generally. It is conceivable that the necessary smoothness, convexity, and noise bounds could be established for other variational families, in which case these optimization results would provide "plug in" guarantees. We suspect this might be fairly easy to do with, c.g., elliptical distributions location-scale families, using similar techniques as done with Gaussians. But may be difficult to do so for broader variational families.
If log p is strongly convex and smooth, another proof strategy is possible: Smoothness guarantees that w" E Was and strong convexity guarantees that w' E {w || - w S d/u)). So one could do projected SGD, projecting onto the intersection of these two sets. The negative ELBO would be strongly convex and smooth over that set, and since || wil2 is bounded, gradient noise is upper- bounded by a constant, meaning classical projected SGD convergence rates that assume smoothness, strong convexity, and uniform gradient noise apply. However, projecting onto the intersection of those sets poses difficulty and this uniform noise bound may be loose in practice.
Several variants of the gradient estimators we consider have been proposed to reduce noise, often based control variates sampling from different base distribution [30, 37, 43, 44, 45, 46, 47, 48, 49, 50J. It would be interesting to establish gradient noise bounds for these estimators.
Various VI variants have been proposed that use natural gradient descent. It would also be interesting to seek convergence guarantees for these algorithms.
VI can be done with "score" or "reinforce" gradient estimators [6, 8), which use the identity
E ((z)) V w zm
((z) log (he
in place of how we used Eq. 4 for the "path" estimators we considered. While path estimators often seem have lower variance, this not always true: Take an unnormalized log- posterrir ((z) where z is scalar. Now, define ø'(z) = <(z) + ve sin(z/e). where e is very small. Then, ó and d' represent almost the same posterior, and score estimators for them will have almost the same variance. Yet, the derivative of the added term is cos(2/c)/u VG, so path estimator for ø' will have much higher variance than one for Φ. This underlines why smoothness is essential for our guarantees excludes posteriors like Φ'. Future work further unravel when score estimators perform better than path estimators.
One attractive feature of VI is that data subsampling can often be used when computing gradients, decreasing the cost of each iteration. While we have not explicitly addressed this, our proofs can be casily generalized, at least for target distributions of the form p(2,x) = p(2) II p(x,, (z). The only issue is to bound the noise of the subsampled estimator. While our gradient variance guarantees use Theorem I from [20], more general result [20, Theorem 6] considers data subsampling. Very roughly speaking, with uniform data subsampling of 1 datum of time, the gradient noise bounds in Thms. 2 and 3 would increase by a factor between 1 (no increase) and the number of data, depending on how correlated the data are-less correlation leads to a larger increase. (More precisely, the second term in Thm. 3 would not increase.) These increases would manifest as larger constants and b in the quadratic bounds, after which exactly the same results hold. However, is not obvious if the bound for ESTL in Thm. 4 can be generalized in this way.
10

---

Acknowledgements
We thank Kyurae Kim and an anonymous reviewer for pointing out an error an earlier version of Lemma 32. This work was partially completed while the first two authors were visiting the Flatiron Institute and was supported in part by the National Science Foundation under Grant No. 2045900.
References
[1] Michael Jordan and Zoubin Ghahramani. An introduction variational methods for graphical models, Machine Learning, 37(2):183-23,, 1999.
[2] Zoubin Ghahramani and Matthew Beal. Propagation Algorithms for Variational Bayesian Learning. In NeurIPS, 2001.
[3] John Winn and Christopher M Bishop. Variational Message Passing. Journal of Machine Learning Research, 6:661-694, 2005.
[4] David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational Inference: A Review for Statisticians. Journal of the American Statistical Association, 1126118)859-877, 2017.
[5] Tim Salimans and David A. Knowles. Fixed-Form Variational Posterior Approximation through Stochastic Linear Regression. Bayesian Anul., 8(4):837-882, 2013.
[6] David Wingate and Theophane Weber. Automated Variational Inference in Probabilistic Pro- gramming, arXiv: 101.111111011 1299 [cs, stat], 2013.
[7] Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley, Stochastic Variational Inference. Journal of Machine Learning Research, 14:1303-1347, 2013.
[8] Rajesh Ranganath, Sean Gerrish, and David M. Blei. Black Box Variational Inference. In AISTATS, 2014.
[9] Jeffrey Regier, Michael I Jordan, and Jon McAuliffe. Fast Black-box Variational Inference through Stochastic Trust-Region Optimization. In NeurIPS, page 10, 2017.
[10] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M. Blei. Auto- matic Differentiation Variational Inference. Journal of Machine Learning Research, 18(14): 1-45, 2017.
[11] C. Zhang, J. Bütepage. H. Kjellström, and Mandt. Advances Variational Inference. IEEE Transactions Pattern Analysis and Machine Intelligence, (888820088-203 2019.
[12] Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Yes, but Did It Work?: Evaluating Variational Inference. In ICML, 2018.
[13] Justin Domke. Provable Smoothness Guarantees for Black-Bo Box Variational Inference. In ICML, 2020.
[14] Akash Kumar Dhaka, Alejandro Catalina, Manushi Welandawe, Michael Riss Andersen, Jonathan Huggins, and Aki Vehtari. Challenges and opportunities high- nigh-dimensional varia- tional inference. In NeurIPS, 2021.
[15] Manushi Welandawe, Michael Riis Andersen, Aki Vehtari, and Jonathan H. Huggins. Robust, automated, and accurate black-box variational inference, 2022.
[16] Edward Challis and David Barber. Gaussian Kullback-Leibler Approximate Inference. Journal of Machine Learning Research, 14:2239-2286, 2013.
[17] Michalis Titsias and Miguel I drzarl-Gredill. Doubly Stochastic Variational Bayes for non- Conjugate Inference. In /CML, 2014.
[18] Kai Fan, Ziteng Wang, Jeff Beck, James Kwok, and Katherine Heller. Fast Second-Order Stochastic Backpropagation for Variational Inference. In NeurlPS, 2015.
11

---

[19] Ming Xu, Matias Quiroz, Robert Kohn, and Scott A. Sisson. Variance reduction properties of the reparameteriizttion trick. In AISTATS, 2019.
[20] Justin Domke. Provable Gradient Variance Guarantees for Black-Box Variational Inference. In NeurlPS, 2019.
[21] Kyurae Kim, Kaiwen Wu, Jisu Oh, and Jacob R. Gardner. Practical and matching gradient variance bounds for black-box variational inference. In ICML,, 2023.
[22] Zuheng Xu and Trevor Campbell. The computational asymptoties of gaussian variational inference and the laplace approximation. Stat Comput, (32), 2023.
[23] Mare Lambert, Sinho Chewi, Francis Bach, Silvère Bonnabel, and Philippe Rigollet. Variational inference via wasserstein gradient flows. In Advances in Neural Information Processing Systems, volume 35, pages 14434 14447. Curran Associates, Inc., 2022.
[24] Michael Ziyang Diao, Krishna Balasubramanian, Sinho Chewi, and Adil Salim. Forward- backward Gaussian variational inference via JKO in the bures- Srres-Wasersstenn space. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 7960-7991. PMLR, 23-29 Jul 2023.
[25] Gabriel Peyre and Marco Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 1 11(5-6):35 607, 2019.
[26] Kyurae Kim, Jisu Oh, Kaiwen Wu, Yi-An Ma, and Jacob R. Gardner. On the convergence and scale parameterizations of black-box variational inference. In NeurIPS, 2023.
[27] Gregory Cooper. The computational complexity of probabilistic inference using bayesian belief networks. Artificial Intelligence, (2(2):393-403 1990.
[28) Paul Dagum and Michael Luby. networks is np-hard. Artificial Intelligence, 60(1):141-153, 1993. ISSN 0004-3702.
Approximating probabilistic inference in bayesian belief
[29] Ahmed Khaled and Peter Richtárik. Better theory for sgd in the nonconvex world. Transactions of Machine Learning Research, 2023.
[30] Geoffrey Roeder, Yuhuai Wu, and David K Duvenaud. Sticking the Landing: Simple, Lower- Variance Gradient Estimators for Variational Inference. In NeurlPS, 2017.
[31] Damek Davis and Dmitriy Drusvyatskiy. Stochastic model-based minimization of weakly convex functions. SIAM Journal on Optimization, 29(1):007-239, 2019. ISSN 1052-6234.
[32] Eric Moulines and Francis Bach. Non-asymptotic analysis of stochastic approximation algo- rithms for machine learning. In Shawe-Taylor, R. Zemel, P. Bartlett Pereira, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
[33] Deanna Needell, Nathan Srebro, and Rachel Ward. Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm. Mathematical Programming, 155(1): 549-573, 2015. ISSN 1436-4646.
[34] Ahmed Khaled, Othmane Sebbouh, Nicolas Loizou, Robert M. Gower, and Peter Richtárik. Unified analysis of stochastic gradient methods for composite convex and smooth optimization, 2020.
[35] Guillaume Garrigos and Robert M. Gower. Handbook of convergence theorems for (stochastic) gradient methods, 2023,
[36] Eduard Gorbunov, Filip Hanzely, and Peter Richtarik. A unified theory of sgd: Variance reduction, sampling, quantization and coordinate descent. Proceedings of the Twenty Third International Conference Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 680-690. PMLR, 26-28 Aug 2020.
[37] Alexander Buchholz, Florian Wenzel, and Stephan Mandt. Quasi-Monte Carlo Variational Inference. In JCML, 2018
12

---

[38] Mohammad E. Khan, Pierre Baqué, François Fleuret, and Pascal Fua. Kullback- proximal variational inference. In NeurlPS, pages 3402-3410, 2015.
[39] Mohammad Emtiyaz Khan, Reza Babanezhad, Wu Lin, Mark Schmidt, and Masashi Sugiyama. Faster Stochastic Variational Inference using Proximal-Gradien Methods with General Diver- gence Functions. In 2016.
[40] K. Fan, Y. Zhang, R. Henao, and K. Heller. Triply Stochastic Variational Inference for Non- linear Beta Process Factor Analysis. In 2016 IEEE 16th International Conference Data Mining (ICDM), pages 121-130, 2016.
[41] Pierre Alquier and James Ridgway. Concentration of tempered posteriors and of their variational approximations. arXiv: 09293 [cs, math, stat], 2017.
[42] Guanghui. Lan. First-order and Stochastic Optimization Methods for Machine Learning. Springer Series in the Data Sciences. Springer International Publishing, Cham, 1st ed. 2020. edition, 2020. ISBN 3-030-39568-5.
[43] Michalis K. Titsias and Miguel Lázaro-Grediln Local Expectation Gradients for Black Box Variational Inference. In NeurIPS, 2015.
[44] Francisco R. Ruiz, Michalis K. Titsias, and David M. Blei. The Generalized Reparameteriza- tion Gradient, In NeurlPS, 2016.
[45] Francisco J. R. Ruiz, Michalis K. Titsias, and David M. Blei. Overdispersed Black-Box Variational Inference. 1603.00603.0/ [star], 2016.
[46] Andrew Miller, Nicholas J. Foti, Alexander D' Amour, and Ryan P. Adams. Reducing Reparameterization Gradient Variance. arXiv: [stat], 2017.
[47] Linda S. L. Tan and David J. Nott. Gaussian variational approximation with sparse precision matrices. Stat Comput, 28(2):259-22) 275, 2018.
[48] Tomas Geffner and Justin Domke. Using Large Ensembles of Control Variates for Variational Inference. In NeurIPS, 2018
[49] Tomas Geffner and Justin Domke. Approximation based variance reduction for reparameteriza- tion gradients. In NeurlPS, 2020.
[50] Masahiro Fujisawa and Issei Sato. Multilevel monte carlo variational inference. Journal of Machine Learning Research, 2021.
[51] T. M. Cover and Joy A. Thomas. Elements of Information Theory. Wiley-Interscience, Hoboken, N.J, 2nd ed edition, 2006.
[52] Amir Beck. First-Order Methods in Optimization. SIAM-Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2017. ISBN 161 1974948.
[53] Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtárik. SGD: General analysis and improved rates. Proceedings of the 36th International Conference on Machine Learning, spotlight paper, volume 97, pages 5200-5209. PMLR, 09-15 Jun 2019.
[54] Sebastian U Stich. arXiv: 2019.
Unified optimal analysis of the (stochastic) gradient method.
[55] Amir Beck and Mare Teboulle. maudient-based algorithms with applications to signal recovery. Convex optimization in signal processing and communications, pages 42-88, 2009.
13

---

Contents
1
Introduction
1.1
Related work
2
Properties of Variational Inference (VT) problems
2.1
Structural properties and assumptions
2.2
Gradient estimators
2.3
Challenges for optimization...
3
Stochastic Optimization with quadratically bounded estimators
3.1
Stochastic Proximal Gradient Descent
3.2
Stochastic Projected Gradient Descent
Solving VI with provable guarantees
5
Discussion
10
6 Notations
15
7 VI problems and their structural properties
15
7.1 Modeling the problem
15
7.2 Smoothness and convexity for VI problems
17
7.3 A case study linear models
19
Estimators and variance bounds
22
8.1 Variance bounds for the estimators
22
8.2
Quadratically bounded estimators
25
9
Optimization proofs
25
9.1
Anytime Convergence Theorem for Strongly Convex
25
9.2
Complexity Lemma for Convex
26
9.3 Proximal gradient descent with strong convexity and smoothness
27
9.4 Proximal gradient descent with convexity and smoothness
29
9.5 Projected gradient descent with strong convexity
32
9.6 Projected gradient descent with convexity
33
10 Solving VI with stochastic optimization
35
10.1
Convergence of the algorithms
35
10.2 Computations for the algorithms
36
10.3 Explicit implementation of the algorithms
37
14

---

Notations
We note Md the space of d x matrices, S the subspace of symmetric matrices, and T the subspace of lower triangular matrices. We use Wa == {(m, C) E W Omin (C) 2 i to denote the set of parameters with positive definite covariance matrix whose singular values are greater or equal to VM
object
description
observed variables
2 E Rd
latent variables
p(2,x) x))
target distribution
= (m, C)
variational parameters (m E Rd and C' E MM, S4 or T
Qu (z) (1/mm, CCT)
Gaussian variational distribution
((w) Exng. logp(z,z))
free energy
h(w) Ex~qa log g(z))
negative entropy
f(w) (w) = l(w)= + h(w)
negative ELBO
w' argmin ((w) + h(w)
optimal parameters
w (z. 0)
MAP parameters (a E argmax p(=,x))
V[g] tr Clg).
variance of vector-valued random variables
7
VI problems and their structural properties
This section contains the proofs all the claims made Section 2.
7.1 Modeling the VI problem
We recall and detail here the setup of our problem. Given a target distribution and observed data a, want to solve the following Variational Inference problem
min KL(ghellp(-1a NEW
(VI)
where qu is a multivariate Gaussian variational family, whose density may written as
exp ( 1- (mm))=))))) w = (m, C).
qu (2) dct(CT))
(17)
We make the assumption (VT) has non-degeee solution, i.e. solution for which the covariance matrix is invertible.
We will impose that our parameters have the form (m, C) E Rd V, where V is some vector subspace of Md, the space of d x matrices. On top of that, we will also want the covariance factor be positive definite. In other words, we will impose that our parameters taken from the
W = {m,, C) E Rd x V C 0}.
For the sake of simplicity, the reader can assume for most of the paper that V = -Md But when coming to analyze our algorithms, we will see that we should take either V Td (the subspace of d x d lower triangular matrices) or = Sdd (the subspace of d x symmetric matrices). This choice dictated by two reasons: guarantees that our problem remains convex (see Lemma 19), and lowers the computational cost of our algorithms (see Section 10.2). Fortunately, requiring C' to be symmetric (or triangular) and positive definite does not change our problem, as stated in the next Lemma.
Lemma 15 (Equivalent problems for triangular/symmetric covariance factors). Suppose that w Rd x ddd is such that qu minimizes (1(g (p(-=))))))))))))))))))))))))))) over Rd x Md, Then there exists w" = =mm", CC)) such that qh and the have the same distribution, with C* being positive definite and symmetric (or lower triangular).
15

---

Proof. Write w = (m, C)2, and take m = m. Let Σ CC Since we assumed that (VI) has non-degenerated solution, we can assume without loss of generality that Σ invertible. Therefore we can apply the Cholesky decomposition and write Σ (*(C") where C* is lower triangular and positive definite. We can also define C (E)1/² to obtain symmetric and positive definite factor.
1/2
Once restricted to W, our problem (VI) becomes equivalent to l(w) + h(w),
min
(18)
(= CV
where
((w) Er
log logp(z, x)) is the free energy
h(w) { -log log det C' if w E W otherwise
+00
This equivalence is mostly standard, but we state it formally for the sake of completeness. Lemma 16 (VI reformulated). Problems (VI) and (18) are equivalent.
Proof. Let (m, C) E W. and let us show that KL(gw ll(((II(()))))))))))()(()))))l ((w) h(w). By definition of the Kullback-Lieller divergence, and using p(z|z) = 2222)//p(22)) we can write
E log (z)(z) = E log zroqu
9m (z)p((
KL(Gw ||((|=)))
=
p(z(x)
z~ww
p(z,
log qu (z) - E log logp(z,, x) + E logp(x)
zr-qs
z~-w
zm
E zrugu
log Zm (z) + l(w) + log p(a).
If we note (X) the entropy of a random variable, we can see that Ez-g, log (z) -z where Z 9u= We can rewrite via an affine transform as Cu + m with u NO, I). so we can use [51, Section 8.6] to write (see also the arguments in [13])
H(z) (u) + log det Cl.
Here (u). the entropy of is a constant equal to (d log(2m). Moreover since W, know that C' > 0, which implies that log det C log det C. Because log p(æ) and H (u) are both constants, deduce that minimizing (L(gwwllp(-laz) over W is equivalent to minimizing / + h. where h(w) = - log det C w W and +oo otherwise.
We emphasize that working with triangular symmetric representation for C is not restrictive, and that doing computations this setting pretty straightforward, as we illustrate next.
Lemma 17 (Computations with trianuularsymmmericc covariance factors). Let V be vector subspace
of M Let & Rd x dd + R and let Φv Rd x V - RR be the restriction of o to Rdd V. Then:
1. If á is convex then dy is convex.
2. If ó is differentiable then Φv is differemtiabl, with Vóv(w) projoa v ( (Vo(w)) In particular,
|| (w:ll < [[]((
3. If 6 is M-smooth over Q. C Rdd x Md, then dy is M-smooth over n (Rd x V).
4. (C) (C + )/2, and p tCi)(C), where tril(C) sets to zero all the entries of above the diagonal.
Proof.
1. This is immediate due to the fact that V is convex, since it is vector subspace of M
16

---

2. <: R x V + +R x M ((m, C) (m, C) be the canonical injection. Then we clearly have Φv o oo t. Moreover,
linear application whose adjoint is x* is exactly the orthogonal projection pr So the computation of the gradient follows after applying the chain rule (wv = [V() = projR'x v (V)(w)) As for the norm, it suffices to observe that the orthogonal projection is linear operator with norm less equal than 1.
3. Using that the orthogonal projection is contractive map gives
[Vvv(w) - Vvv(w')l = || x - p
[Vo(w) - Vo(w")ll - M|W w'll
4. This standard linear algebra result.
Smoothness and convexity for VI problems
Now we state the main smoothness and convex properties for VI problems. We will see in our analysis that the smoothness properties often revolve around the following subset of W
1
WM {(m, C) E W T (C) 2
VM
where PCin(C) refers to the smallest singular value of C. Lemma 18 (Smoothness for VI). Assume that lo, x) is M- M-smooth. Then
1. is M-smooth.
2. h is differentiable over W, with Vh(w) (0, (0,-proiv(CC) (O.-projv(C-T)) for all w e W.
3. h is M-smooth over WM {m, C) e W Jmin (C) 2 VMd
1}.
Proof.
1. Combine [13, Theorem 1] with Lemma 17.3.
2. Combine the fact that the gradient of log det at a matrix X is -x- with Lemma 17.2.
3. Combine [13, Lemma 12] with Lemma 17.3.
Lemma 19 (Convexity for VI). Assume that - log x(),)))))))))).)))(. is convex (resp. p-stronggy convex). Then
1. is convex (resp. p-strongly convex).
2. V = then is closed convex, W is convex, and WM convex and closed. Moreover,
W {m, C) E Idd Sd Amin (C) > 0} and WM (m, C*) E Rdd x Sd Amiu ( VM
3. VfV - then h closed convex and is convex. Moreover
W (m, () Ed Td C has positive diagonal }.
4. Vf Md then h can be not convex.
Proof.
1. Combine Titsias and Léraill [17, Proposition 1] (resp. Domke [13, Theorem 9]) with Lemma 17.1.
17

---

2. Convexity of h follows from [52, Example 7.13 and Theorem 7.17]. Convexity of W comes from the fact that is essentially the of symmetric positive definite matrices (which is equivalent to Amin > 0). Let now prove that WM is convex and closed. First we notice that, because of symmetry, we can rewrite Was as
1 VM
WM {(n, C) TR x St Amin Ammn(())
This set is closed because Amia is a continuous function. To see that WAs is convex we are going use the fact that Amin a concave function. Take C1, C2 E Was. a e [0, 1] and write
(((1-a)CC Amin ((1 a) + aC2) = min (((1 - a)C aC2)2,,2)
||||=]
min (1 - a](C] a(C2z,2)
||]|=]
min (1 a)(C12,2) min aC 2)
||||=
||||=
1
- a.) (Ci) + a: (C2) N VM
3. Here W = {(m, C) E Rd Tá C- 0}. Using Sylvester's criterion, is easy to see that lower triangular matrix is positive definite if and only if its diagonal has positive coefficients. In other words, W (mm, C)) E Rd T : Ca > 0}. which is clearly convex. It remains to verify that h is closed convex even W not closed. For w (m,C) E RRd Td, we can use the triangular structure of C to write
d
h(w) = - log det C' + 6w(w) = Σ - log Ci + tc0,
where used the notation 8. for the indicator function of set A, which is equal A and too outside. Since - log(t) + 6o) (t) is a closed convex function, we see that h is separable of closed convex functions, and thus is itself closed and convex.
4. Let 2, let be the identity matrix, let R (-16) be rotation matrix, and define Φ : R + R with v+1= h(T + LR). Observe that v is well defined because IR > 0 for all t R. On the one hand. if h was convex, then & would be convex too, as it is composition of convex function with the affine function I + tR. On the other hand, we can compute explicitly that (t) - log(1 + t²), which is not convex.
Lemma 20 (Solutions of VI). Assume that log p(-, 22 is M-smooth, and that V M sa or T Then argmin(f) C WM == {(m, C.) E W: Jmin (C) 2 VM
Proof. u = (m",C** E argmin(f)) f), then w² is in the domain of h, which is W by definition of h. So it remains prove that amin (C*) 2 1/ /M. When V = Md, this is proved in [13, Theorem 7]. Suppose now that V = Sá, meaning that w* is minimizer of f over W C Rd x Sd. We see that w² must also be minimizer of f over Rd Md, because if there was better solution u (m, C) with č E Md, could consider (ČCT)1/2 S which would itself be better solution than C", which contradiction. So
apply [13, Theorem 7] to w" to conclude that m (C*) > 1/VM. If V Td, we can use the same argument, using this time the Cholesky decomposition of ČÒT
Lemma 21. Let (m, C) with C' invertible. Then log is omin -
Proof. According to the definition of qu (see (17)), the Hessian of log qw is (loggg)) -(CCT)- and so V (log (wwlllggwwweeeee (z)]l Tmax ((CC )-1) Jmin (C)-2.
18

---

7.3 A case study
linear models
Table Table of models
Model
Smoothness constant
Strong convexity constant
Convex
Bayesian linear
1 +
x² 0 max (A)
1 o² m (m))
yes
regression (2 = !)
Logistic regression
m (A)²
yes
(I=
Heirarchical
exists
n/a
yes
logistic regression
Definition 22. twice differentiable function f is M-smooth if -MT 3 V25 MI.
Definition 23. A twice differentiable function f is p-strongly convex V2 2 ul.
Theorem 24. Take generic ii.dd linear model defined as p(z, = p(2) II p(zn|z, an), where p(z) = =[ Σ) and
p(znlzz, an) = exp (-o (z an Rn))
for some function 6. Suppose that the second derivative of & with respect to its first (scalar) argument bounded by muin S 6² S maax Let
Amin E + Omnn A AAT)
8
Xmax (Σ-1 0 O AAT).
where A a matrix with on the n-th column. Then
1. Mogg(zze0 2) M-smooth over z for M max (1u) 18)).
2. If 2 then log pezpp(z,1 x) is convex over z
3. If > 0 then - plz, x) is p-strongly convex over
Proof. It is easy show that
V20(2 an an. xn) (2 annnsn)aann
from which it follows that
? log = =--1 x))
Σ
an
a
Now, take some arbitrary vector We have that u (oggoog(2)))))
-- Σ bn TTpll
uT T Omin Σ 1oTTa amal 11
TT
Σ-1 + Amin Σ
Eanan
V.
uT (Σ-1)
+ Omnn AA )
U.
Similarly, we have that
v (oVioog(zzz))
vT -- Σ On 1100 Tanl
s u I 1 + mmax Σl |1oTa. anll 1|2
TT
Σ 1 + Omax
Eanam
TT (Σ-1 AA v.
19

---

Putting both the above results together, we have that
Σ + 0 min AA 1 V logg((,, =-ー + 0max
This implies that all eigenvalues of - V log are in the interval [e, (). All the claimed properties follow from this:
For any eigenvalue À of V ((), we know that [A( logg2,)((,, - M = max (11) Thus, (og,(2, M-smooth over 2.
μ 2 this means all eigenvalues of -V2 log are positive, which implies that log p2,,(2 is convex.
If y > 0, then all eigenvalues of V2 log(z, x) are bounded below by a, which impiles that (ogg((z,)) is j-strongly convex.
Corollary 25. Take a linear regression model with inputs on E Rdd and outputs In, IR Let
p(z) = (2) Σ) and p(æn|z) = N (x| _ an o²) for some fixed a. Then - log (22,00 x) is p
strongly convex for μ Amin (2-1 + [AAT) and M- -mmoooh for M =mmax22-1 {AA
Proof. In this case we have that
@ On ani In) - logp(#nl22)
20 (on - 2 =Ta,) =Tan)² = log
Or, more abstractly,
Φ(a,8) =
(a - -3)²
log (og(2To²))
20
2
1
)((a,))
(a 3)
0²
1
&" ("(a, B)
d²
Thus we have that 0min = 0 max zz. It follows from Theorem 24 that - -logp(2, x) is strongyy
convex for μ Amin (E- + / AAT) and M-smooth for M Amax (E-1 1114440 (AAA)
In the particular case where Σ I, we get that Amiin (I AAT) =1 a3 o min (A)² and M Amax AAAT
_ max
Corollary 26. Take logistic regression model with imputs On E Rd and oulputs in E f1,1,+1 Let p(z) = N[2200, [1000) E) and p(zn 1|2) = Sigmoid(zn on) = 1/(1+ + exp(-2n2
22, an,))) Then
loog p(2,x) p-strongly convex for μ Amiu (Σ-1) and M-smooth for M Amax ).AT).
Proof. In this case we have that
66z an, = logp(#n(z)
= log (1 + exp(- x,, Tan))
Or, more abstractly,
6(a,3) log (1 + exp( exx((-a33)
@ ((a,3)
+ exp((.)
8² e
("(a,
(1 ex
20

---

The second derivative is non-negative and goes to zero if a - -00 or a - +oc. It is maximized a in which case ø" 3.3 But of course, in this model, the second imput is In (-1,+1}. meaning / always. Thus we have that imii 0 and mmux
follows from
Theorem 24 that logllgpp2, p(2, a) is p-strongly convex for μ Amin (I-1) and M-smooth for M Amax (2- + )
In the particular case where = I, we get that p and M Amax + (AAT) = m (A)?.
Theorem 27. Take a heirarchical model of the form p(0, , 2, 2, p(0) II: p(z; 10)p((;00 ((2:1)0))111 21). where all log-densities are twice-differentiably Then, log p(0, 2, is smooth joinly over (0, 2) if and only if the spectral norm of the second derivative matrices Voo p(, 2, 2),, VesT log p(0, a 24 Tz), and
V = log Rig T4) are all bounded uniformly over (0,z). 2).
Proof. The following conditions are all equivalent:
1. p0g, x) smooth as a function of (0,
2. V lo 1(0, 2, x)) 1|1 is bounded above (uniformly over (0,z), where V denotes the Hessian with respect to (0, =) and |l-112 denotes the spectral norm.)
3. log 2, || is bounded above (uniformly over (0, z), where V denotes the Hessian with respect to (0, 2) and ||-|| denotes block infinity norm on top of spectral norm inside each block.
4. V logp(0, 1| above (uniformly over (0, 2))
Va., logp(0, and
V2 [.]]
logp(0,
are all bounded ||
5. Voe logp(0, above (uniformly over 2))
V2 .hg log (0g 20, and
V2 #ta log,(0,,, 2,) all bounded
6. G0T
dez logp(0, 2, and
V2 logp(0, 2, are all bounded
above (uniformly over (0, 2))
7. 2 TTT log 122 V log p(0, 2,, (0, above (uniformly over 2))
and V2 V2..T log 24, 2., are all bounded
Corollary 28. Take a hierarchical logistic regression model defined by p(0,2,x) pl. 11, p(z))lz)) where (0, zi} E Rdd Inj [[1,++],, and
0 (0,
24 ~ (0,A))
Zij p(zzg(zz)
p(z;j
1|2,)
Sigmoid(Z:j2 o:j),
where Sigmoid(b) = 14exp(- 6) Then, log p(0, x.) is convex and smooth with respect to (8, 2).
Proof. It's immediate that - log p is convex with respect to 20, We can show that this function is smooth, by verifying the three conditions of Theorem 27.
2 G0r log pogp(0,,.2) = -
is constant.
V2 GzтT logp(0, 2;, is also constant.
21

---

p(2z(0) is a Gaussian with constant covariance, so is smooth with constant independent of B. And we can write that log (2;) log (1 + exp (- axz)): This function is also smooth over 2; with a constant that doesn't depend on 0. Thus log p 2222 (0) + Σ, log p(zg)(z)) is smooth over 2, with a smoothness constant that is independent of e.
8 Estimators and variance bounds
We remember that unless specified otherwise, we consider that w = (m, C)) E IR. x V, where V is vector subspace of Md.
8.1 Variance bounds for the estimators
IRd - R be M-smooth, suppose that Theorem 1. Let (w) Cu + for = (m. C). Let Φ :R. stationary m, and define t (m, 0). Then
E w=N((0.1)
(T S (d + 1)M² I ml12 (d+ | (d+3) || tw-wl12
Furthermore, the first inequality cannot be improved.
Proof. This bound is proven by [20, Thm. 3] in the case that V = -Md. We conclude that this bound holds for every subspace V of M by using the inequality in Lemma 17.2.
Theorem 2. Suppose that log ax) is M-smooth and has a maximum (or stationary point) at m. and define w (m, 0). Then, for every w and every solution w" of the VI problem.
(d+ 3)M |/ 011 E gemergy (u)l s
(6)
2(d + 3)M w 11 2(d 3) -0112-
Proof. is direct consequence of Theorem and Young's inequality
genergy 112
(d 3) ||l 0
= (d (d+ ||w- - =" + + - -
< 2(d (dd+3)M² || w w01 2(d + 3)M² I| - -
Theorem 3. Suppose that log p(-;x)) is M-smooth, that is maximal at m. and define t = (m, 0). Then, for every L > 0, for every E Wt. and every solution w" of the VI problem,
Elgeet (u)l1
2(d + 3) 1 + 2dL
212
(8)
W 4(d + 3) u w1 + 4(d + 333M2 - - + 2dL.
Proof. The difference between gent and genergy is the addition of the constant vector Vh(w). Thus, we have that
E gent
2 = E genergy + Vh(w)l
2 Beuergy 1| + 2E [Vh(w)
M 2(d + 3 || m11 [[Vh(w))
remains bound the final term. We can this using the closed-fomm expression for the entropy gradient (see Lemma 18.2), plus the assumption that w E Wz C W, to see that
22

---

VVh(w)]2
= Iprojj(
S ||o I
E (c
E (C)-
Σ VL
= dL.
We conclude with Young's ss inequality
E geng
2(d 3)M² ||--
----wly u11 dL
2(d + 3) | - "² + w" w112 + dL
4(d 3) || ||w- w + 4(d + 3) 1| - 201 dL.
For the next Theorem, we will need the follwing technical Lemma Lemma 29. Lef NO0 I). A E M and b E Rd. Then
E Au 3ll + ||u|] (d + A|
Proof. Develop the squares to write
|||33) Al
Aull 2(A b, u) + ||61 || 2 + 6, ||b||=
Since IEu 0, immediately see that B(A 6, u) 0. Because of symmetry, the third-order moment Euilu| 1 is also zero, so we obtain B(A b. u||u||? = O. We also have the second-order moment E||u|| d, which means that |||b||||u|| d|bb/22 For the other terms, we use the identities
Euu I, Euu Tuu = (d+ - 2)T
from [20, Lemma 9] allows us write
E 2ul E Etr(A Auu 7) tr A) ||A|) E |u
1| = tr(AT Auu T (d+ |
The conclusion follows after gathering all those identities.
Theorem 4. Suppose that log pggp(..) is M-smooth. Consider the residual r(z) == log p(z, x)
log qu- (z) for any solution w" of the VI problem, assume that it has stationary point m, and define v = (m, 0). Then r is K K-smoothfo for some K E [0, 2M), and for all w E WM.
E E|gSTLI
I 8(d+3)M2
w01 2(d + 3² ||
(10)
4(d (K² 2MM || w0 1 + 4(d 3) I ui
Moreover, if p(1()))) is Gaussian then K = -.
Proof. Let E War be fixed. For the duration the proof, we take to be a second copy of w (held constant under differentiation with respect w). Now, we can rearrange the estimator in Eq. by making appear qw- (Tw (u)), namely
(u) = V., log go (u)) V log Jn (Tm + V... log Lu (To) Tw (Tw(
()),
p(Ti. (u), x))
23

---

Introducing ((2) == log (2) - log fw- (z), and recalling r(2) = logp(z, x) - log tw- (2), we have
gSTL = Vo - V (TT
We assume that log p( is M smooth, and that V = -M, T or S which implies that w E War (see Lemma 20). This in turn guarantees that log qu- is also M-smooth (see Lemma 21). Thus, is the sum two M -smothh functions, so is K Kssmooth with K [0, 2MM.. Using Young's 's inequality, together with Theorem applied to r, obtain
(u)l
2E T (u)) + 2E, VV (TT.
2(d + 3)
0 w|| + 2 (
(19)
where (m,, 0). with mn being a stationary point of Now remains control the last term of inequality (19). Apply the chain rule onto ofT (u)) to write (with the help of Lemma 37)
oTo (V. (T (u)), projy (V (V=o(T(vv)uu))
This gives us directly that
|V (V)((T
= V (T= (1 projy (1
S (T) [V (V) '
= (T |ull?)
(20)
where we used that || T1 ||> ||| /1l2 for any vectors v and u. Computing V.ø amounts to computing the gradient the Gaussian density lw- From definition (see Eq. 17) we have:
V log qw (z) - -CCCT - m) - (2 m).
Recalling the definition of = log 1o:00-lla -log qm and writing u = (m, C), w² (m,, E. C.C. we obtain
V=0(2) (E,1 Σ + ---mm Σ-1
In other words, we have
V. v.o(Tu(u)) = (2, - 1) (Cu + m) + + m - '.' m. = Au + b,
(21)
where := = (E,1_ E D and = Σ m m. We can now combine (20) and (21) and use Lemma 29 to write
E. V o
Eu 61 ||u||²) (d+ 1 + (d +3)1|
S (d+ (1b12+11A|13)
(22)
Let introduce the function ::: Rd Md RU(+ defined by s(w) == Fe- log que (z)
log qu- (z) if C' is invertible, +o0 otherwise. This function is nothing but the Kullback-Lieble divergence between the two Gaussian distributions qu and qu-, and can be computed in closed-form on its domain:
x (w)
2 (log de.( - log det d tr + (E- (m. - -m), (m. - -mm))
where we note before w (m, C)2 w" = (m,, E= CCT, E, = =C.CT. This allows us to compute its gradient at w E WM:
V Vmni(w) Von(w)
Σ, m))
-1 log + Vo
= - --
We see that Vmn(ww= and Vor(w) A, which means that ||b||+ (22) and the fact that V (i = 0 (because w" is minimizer of k). we deduce that
|| = [ From From
E. V o (d V x(w) - V (wn(w") 21.
24

---

Now we remind that w" E War which implies that log 9u- is M-smooth. This means that Exoqu log qw- is M M-smooth as well, according Domke [13, Theorem 1]. On the other hand, we know (see the proof of Lemma 16) that 1Zz-q loggu= qu. h(w) which is M-smooth on War (see Lemma 18.3). All this implies that K. 2M-smooth on WM. from which we conclude
Bu V or (T. (u))i² (d+3))lV ( V=n(w=))l= < (d + 3344 w" w*|l²
The Theorem's main incquality (10) follows after plugging the above inequality into (19). The second inequality follows after using Young's 's incquality
E 1
s 2(d K² || vil12 + 8(d+ 3) || w 1
(dd+3) ||l =² ---= -----1i
m112 + 3(d+3)M) || w0 1012
S 4(d+ 3) 2 ||0 - w + 1 -0l2 + 8(d + M2 w w 1
w01 4(d+3) 3) (K² + 2M22 1| + 4(d 3) I|0" 10
To conclude the proof, it remains check that K 0 whenever p(-[x) is Gaussian. In that case, because our main objective function f(w) the divergence between qu and p(-[z), is clear that is minimized whenever Sw p(-1x). Therefore, without loss of generality, we can assume that p(-[z) = qw- This implies that the residual r is constant (r(2) = logp(z)) p(z)), and so is O-smooth.
8.2
Quadratically bounded estimators
Theorem 30. The estimators obey E lg S a w²12 + b with the the following constants:
Estimator
valid
genergy
2)
2(d 32 2 [w- 00/2
any
gent
4(d+3)M²
4(d 3) |
011
dL
Wr.
ESTL
4(d + 3) (K2 2M²
4 + K2 I|w"" m11
WM
Proof. Combine the results of Theorems 2, 3 and 4.
9
Optimization proofs
This section contains all the proofs for the Theorems stated in Section 3.
9.1 Anytime Convergence Theorem for Strongly Convex
Theorem 31 ( Thhoremm 3.2 in [53]), Suppose
given sequence of wt iterates such that for step size Y < 22 and given constants μ > O and o² > 0 we have that
E
[|lwtt
<
(1 - ME) [Ile* w= 112 2022?
(23)
By switching to a decaying stepsize according to
fort <
2C
7/-
1 2t
for
4 (t + 1)2
where t* = 4[C/H,, we have that
16|_//u 12
a²
I| T+1
201
<
wl10
w" 112
(24)
(T' 1)2
T +
25

---

Proof. We divide the proof for steps t that are great or smaller than t². For t S t², we have by unrolling (23) starting from t" ........ 1, with i't = 7 gives
(1 (1-20 (1- 1| w w²12+ + 20
E |1w0" w" 112]
(25)
Now consider t 2 t" for which we have that "t = 2t-1 which when inserted into (23) gives (t-1)z
t²
+ 0² 2 (2t 1)²
E [1| [ w0
(t+ 1) [llo w w0]
2² (1 + 1)a
Multiplying through by (t + 1)² and re-arranging gives
2 (2t + 1)² 422 (t + 1)2
(t + 1) [1 w [13]] x
E I||w* - 11]
+ 0²
8
E [||w* w² 11]
2
where we used that gives
2+++1
< 2. Summing up over = ......, T and using telescopic cancellation I1
w" w11 I (t")²E Il w² w1] (T (T") 802
(T 1) w T+I
412
Now using that for & S tt we have that 7 and (25) holds, thus
S
(T 1) [1l - --=-1==11"" 1=1
(1))² (1 11-7/4 711)" "I10 w011 + 2022
+ (T - (*))
11)
11
(17)² Ilwe Ilwo -w1= 1=11++200 +2o22 (1)
80²
(²)
11
122
where we used that (1 (1-7p)= S 1. Substituting 7 22, 1" 4 4..// [p] and multiplying by (T' + 1)2 gives
[1lu - =11
16 16
11
112 +
a²
(8(T - -²) (1²)2 1
(T + 1)2
1T + 1)2 a2 8 21")
16 16[ c/p, 12
11
+ μ² (T 1)2
(T + 1)2
g2
(T 1)2 |!w0
w 11
112 +1
which concludes the proof of (24).
9.2 Complexity Lemma for Convex
Lemma 32. Suppose we are given a sequence :.. and constants a, b, B such that
E [[] - inf f S a o w"
+ b
0T
holds for
== + By
Ify = A for A 2 i and T 2 2, then
A - w0112 +b) b).
[f(xT) inf f S
(2all00
VT
Thus, T 2) iterations sufficient to guarantee that E [f(w]]] inf <
26

---

Proof. First we show that I-6r
2 is equivalent to
log
T
(26)
log + BA2/TT)
Indeed this follows since
1
2 2
IT
0"
log 2
T 2
log 1/0
log 2
T 2 log(1 1 + By
log
T log(1 + BA²/T)
Now note that (og(1+x) 2 14x this gives that
for 2 0, which implies that
S 114 ! for x 2 0. Applying
log 1+22
1
T
log(1 (1 + BA2/T) S1
BA?
So to guarantee 1-8r $ 2, from (26) it is sufficient enforce that
T 2 1+ BA²
T
log(2)
I+
BA²
Assuming A2 22/B, this last condition holds if
1
T 121+ BA² 1
(27)
Since also impose that BA? 2 2 we have that BA²-1 < 1 thus for (27) to hold it suffices that
1
A 2 2. Substituting V and the relaxation I-0r S into the original bound gives the claimed result.
9.3 Proximal gradient descent with strong convexity and smoothness
We start by recalling some standard properties of the proximal operator.
Lemma 33, Let > and w E argminn l(w) + h(w). Assume h(w) is convex and that I(w) is continuously differentiable at w". Then
1. For all w, w' and 7 > 0, ||prox;; (w) - prox, (u')ll2 I w|w - w'll2
2. w" -prox,, (w* VV((w*))
3. = prOx. (w) 4> W1-2 E @h((v)
Theorem 7. Let 1 be ju-strongly convex and M-smooth function, and let i argmin(l). Ler h be proper closed convex function, and let w' argmin( + h). Let (w") be generated by the Prox-SGD algorithm, with constant stepsize E (0, min Za 41]. Suppose that g is quadratically bounded estimator (Def. 5) for VI with parameters (a, b, w Then,
E 1 T+1 w" 1 < (1 - -74)" w w/ +
2y (b M 1 - 0016)) -00())
(1)
27

---

2t +
Alternativel, if we use the decaying stepsize " min
then
2a (t + 1)2
Bllw" -w" w=1 16[a/p²13
1
w' + N²T
8
(b + M² 110" - -0lā)
(12)
T
In both cases, T (EE iterations are sufficient to guarantee that Bllw w w011 112 Se.
Proof. We start by using the non-expansiveness and fixed-point properties of the proximal operator write
1 (
w012
|prox (w" - -9'' - proxyn (w² - >V((u"))
v |wo w² + (V((ww) 9²) g')
||0 - w011 + [|V[((") - 27 (w² - w², vI(w")- - g¹)
Above the first line uses the fixed point property of the proximal operator that w" prox,A (w* VV((w), while the second line uses the fact that the proximal operator is non-expansive. Now, we apply our assumption on the gradient estimator write
[V((")) - g e²ll w
z 2 [|V[(w>)] w0 + 2
|19" 112 wo
2 2|[V((w")!l )1{ + 2 (a | - w 112 + b)
< 2M | - !1l + 2 (a [a|lw'= - w0 112 + b)
For the last term, we have by strong convexity and the fact that gf is an unbiased estimator that E Kww' w", vl(w") g²) [e]] - (w" w', Vl(w=)) (((w))-Vi ((w""))V(ww')
I -H -m|lw' w 1l
where the second line follows from the fact that 7 is a strongly convex. Putting the pieces together, we get that
|| w: 112]] S E 100 - w= 112 27, 2-222 N2 1100" - -il12 + 2y (a (allue w0 10(2+1) -27 Ilwe" w-
(1 27/4 + 27 a) E 1|10" w01l + 2?² (o + M² I|u" 00l22)
Now the following conditions are equivalent:
27μ + 272a a
< mp
27μ1 27²a
7/
2/0 - 2y
2/1 /
2
/ Z 2710
|
7 2a
This means that
E
[1l
w0 [2] 5 (1 7/1) E
| w + 27² (b (b+M M2
1/ "1u"--ij -w(()
Now since E 1
1
5t - 1-74 2 we can apply the above recursively, which gives that
- -1 27² (b+ M²2 Il10 = Ea-
E
| ++1
w01
(1 710 | -
k=
But can bound this geometric sum by
yu)* 11--1-10
1- (1 7p)²
1
7/4
7112"
28

---

leading to the result that
2 w= 112] S 7/0)" 11100
w==12 + 27 (b M 1| .0(l)
B
To prove the anytime result, we can now apply Theorem 31. To do so, we need to map the notation we use here to the notation used in Theorem 31. The mapping we need is L
and o² je?
(b+ M Ilw" por--uj -wl2) which when inserting into Theorem 31 gives the result.
Corollary 34. Under the conditions of Theorem 7, i = AlogT for S A, and T is large enough
that
20A
log
2A lologT w 1125 T ||00 -w" ==1
(b M2 I10 -wll) = -2 log
E ||
T+1
μT
Proof. The previous theorem only be applied when y 20 i.c when Alog I
2u or
20A
< T Supposing that is true, observe that (1 (1- exp((-x'T) from which it follows that 1
log
(1-744² < (-x
THa S T
9.4 Proximal gradient descent with convexity and smoothness
Theorem 8. Let I be proper convex and M Mssmohth function. Let h be proper closed convex function, and let w E argmin(l + h). Let (w")rew Jeen be generated by the Prox-SGD algorithm, with constant stepsize 7 E (0, 33 Suppose that g quadratically bounded estimator (Def. f) VI with parameters (a, b, w Then,
(1 (1-01) 0T)) -0 + b)) -b).
11 110--"0000 w0112
E [f(w²) 1
inf f] S 7
1 20r? and in def ΣT=1 wcti
where 0
1== ac+1 In particullr, if
then
VaT"
I22 2 . Thus, T = O(E 2) iterations are sufficient to guarantee that E [f(w") inf f Se.
E [f(w")] - f I VaT (2allw² 12=12=+0) -b) VT 2 max
Proof. Let us start by looking at || - w01 11w0* - =='1 Expanding the squares, we have that
1
1 I1w0 w0112
-1 1 r0112 P-
w² (u' u0+1 1+1 w
2%1 |11 w²
2%
27/
7
Since r+1 prox (w' ng"), we know from Lemma 33 that
(w0
(w" Ne9')
++1
Oh(w"+1)
7t
and therefore that
w² #11
g' @1(w
7
where (w) is the subdifferential of h. Consequently there exists subgradient b++ ah(w++1) such that
1
1 - -11 _ I|o6 w"
(28)
27:
2%,
-1
|
20 - (g' + bi+1 w w²)
27/6
-1
1
w01 2 - (g) - ve(w"), 111
w - w²) -VV((V() + b+1 w w²).
2%
29

---

We decompose the last term of (28) as
(Ve(w")+bb* w w²) -bbtt1 with pott-="") (Ve(w"), a/+t-w") + (Vl(w"), w² w²))
(29)
For the first term in the above we can use that b E Oh(w is subgradient together with the defining property of subgradient to write
+b w1+ w²) (b++1, w²= - I h(ww) h(w
(30)
On the second term we can use the fact that I is M-smooth to write
M || w'1l + E(ww)) f(w
(VV((""), w - w²) -
(31)
2
On the last term we can use the convexity of f to write
(Vl(w"), w - w²) S f(ww) - l(w').
(32)
By inserting (30), (31), (32) into (29) gives
M
- -(Ve(w') + ++1 w w²) & h(w") + f(w") - h(w*+1) - -
+
1| w011
2
M
= =s(w") - s(w + Ilw fu'tt=w011² w0
2
Now using the above in (28), and our assumption that 7MM s 1, we obtain
1
1
w²
2%, I|w w²11
2%
$ -1
|1 - w/112++ +
1w'tl-ui - (f(w - inf f) - (g' ve(w"), -²")
2/1
2
< ((((wtt)) - -inf f) (g' - Vl(w"), ("t""),+++1_ w")
(33)
We now have to control the last term of (33) in expectation. To shorten our notation we temporarily introduce the operators
T(w) def
7t V((w),
[(w) def r1g"
(34)
Notice in particular that w = prox, (w*)) (w))). So the the last term of (33) can be decomposed as (g' V((ww)) w0t1 1-+1--ir -w") -(g² - Ve(w"), prox (f(("))) prox (T(w²)))
-(9" - ve(w"), prox,, (r (T(w")) w²)
(35)
We observe that the last term is, in expectation, is equal to zero. This is because prox,, (T (w=")) - is deterministic when conditioned on wf. Since we will later on take expectations, we drop this term now and keep on going. As for the first term, using the nonexpansiveness of the proximal operator (Lemma 33), we have that
Vl(w"), ) ,r (Tw")) prox
(T(w")) x
S ||g - Ve(w - T(w*)l
nlll - ve(w
Using the above two bounds in (35) we have proved that (after taking expectation)
-Bíg' - Vl(w"), w w²) < T ve(w")ll V [g"]
Injecting the above inequality into (33) and multiplying through by 7e, we obtain Fl w0 1|2 TE [llw*
w" 112] < - [f( inf ʃ| + 2 [o]]
(36)
From now assume that the stepsize sequence is constant Yee 7. Now, applying our assumption on the gradient estimator, we have
v [ < Ellg'' < aE w' w1|12+
30

---

Inject this inequality into (36), and reordering the terms gives
1 R [|lutt] - -w1 + 7 (B ||me - w00 -001 +1) -1(])])))))))
ylR [
inf sl
E [|lo' w w²113]
2
2
1 1E [1/00t1 -=1 + 726
2 (1 + 27(a2 a) E [[]lw* - "" w²1
-
2
CEE [llee* w=w1 ==11²] - 2 [llw w²113] +
(37)
where we introduced the constant C def introducing sequence ax+1 > 0 of weights. Multiplying the above inequality by some ax+1 0, and summing over O, O....,T T 1 gives
1 + 2a72. We will now do a weighted telescoping by
T-1
T-1 |110 w=11] +2 2
y Σ B [[w inf /| S Σ a141
(CHEllw²=
- w01 "I
Σ
act1
To be able to use a telescopic sum in the right-hand term, we need that ar+1 C = ar (an idea we
def J. which allows us to write borrowed from [54]). ). This holds by choosing Ox, Ct.
T-1
T- w w=1l3] + 72
o++1 E [ - in v
Ctt
[llw*
w w1]
0t++1
E
[|/
Σ
0++1
2
2
t=0
t=0
t=
T-1
a0
w100 w² - aT
E
[|lw" w
w=1
+ 22
04+
2
2
0=0
a0
2 1100 w0 3 + 2
01+1-
Let us now introduce w" def Elsr aerr ++1 Ei ax41
using that ap 1, and by using the convexity of
with Jensen's inequality, obtain
Ei=o a (w
T
E [[(w"]]
E
inf f
Et= ar+1
T
a+++ E [f(w++1]]
0 0
inf f
=i=0 a++1
T 1
(=0 Cue+ E [[(
- inf /
Σ/ 0141
Σ
1= a+++ EE [f( inf fl
T-1
ΣT ac+1
w0ll 112 + 76 EtTo a111
2
<
Σ ET-1 ac+1
1
mmll - =u' u01 + yb
T-1
27 21=0 a1+1
||w0 w w11 yb
27 Et=
a1+1
1|100 w 012 + 7b.
27 2tmo
0x+1
To finish the proof, remains to compute the geometric series
Σt a:+1 which is given by
T-1
- -
11-1111100
11-1/C" 1 1/CT
-
C
31

---

Using the above and substituting for C' = 1 + 2ay2 and 0 = 1/C gives
C -
E [[f(m")] inf] f
27 (1 1/C²) 1m10 w +
(1- (1-07) ) 1w00 w=11 +7b
Finally by applying Lemma 32 with where B = 2a, A = v
1 and thus
A
Ja
VT
VoT'
Note that with this choice of 7 have that
M
+ T>
VaT
M
Thus for T 2 2 the result of Lemma 32 gives
E [[s(w")] inf f s VavT
-all==-w'112-
+
vT
9.5 Projected gradient descent with strong convexity
Theorem 10. Let W be nonempty closed convex set. Let f be u-strongly convex function, differentiable on W. Let w" = argmin Let (w" twe)en be generated by the Proj-SGD algorithm, with constant stepsize 7 E (0, min{ # ]]]. Suppose that g¹ a quadratically bounded estimator (Def. for f with parameters (a, b, w*). Then,
E B||w²
w' w'1
1| w w²|l= +
2yb
(13)
P
11 21
(t + 1)2 } then
Alternativel, we use the decaying stepsize "/t = min
2a
32a
16b
E
[
w 1²]
w" w²1l 1²T
(14)
p2T2
In both cases, T = 2) iterations are sufficient to guarantee that El| - w w12 S..
S
Proof. Plugging in one step of Proj-SGD we have that
|| | w² = |/roo - Y:gg)) projw
< ww =" ?
where we used that projections are non-expansive. Taking expectation conditioned on w', and expanding squares we have that
E || w/ wt]
||="1|/ w( < I|w" w/ w011 2T. (VS(w'', w - -w") + 7/ E
< (1 - wrx) l| - w01 - 2 (J(")) f(w*)) )) + arillu² - w + bn?
where we used E [g"] = fw") and definition 5, and that f is j-strongly convex. Taking full expectation and re-arranging gives
T. Flll wr|l + br?- 2% Bif(ww)) f(w)) (1 re(p ary:)) EE I|w0 w²1
Now using that 7t I #h we have that
1 - 1-7e - -7x) < 7eff
2
and thus
2%, Blf(w') f(w=)) s x11- 147
El|u0* 1oo" w0 112
1 r+1 w'll + +br (38)
2
32

---

Re-arranging we have that
E w ++ w²l
-
107
E Iw' - w¹12 2n (l[(w")))))))))) - f(w")( + by
2
1
E Iluet w²1 + b
(39)
2
Using 7t = y constant, and since 7 S
1
we have by unrolling the recurrence that
E Bll
w0
10= 5
(1
M4
| w w²1
(1
by
2
2
k=0
Since
Σ
)
7/1
2
72 27² 11 - (1 -1/4/22
27
27
7p
7/
1
k=0
we have that
E (+1 -
B
w01 -
1 (11
1+1
w
w=1
w²
2yb
2
1
which concludes the proof of (13).
To prove (14) we will apply Theorem 31, for which we will switch our notation to match that of a
Theorem 31. That is, substituting L a o² 2 and y into (39) gives
E Bl
w 11 s (1 (im) I w 20
(40)
which now holds for & 22 and fits exactly the format of (23), excluding the additional hat μ. Thus we know from following verbatim the proof that continues after (23) that for a stepsize schedule
of
1
for t < t"
2C
7/-
1
21 +
t"
j (t + 1)2
t = C[C/ we have that
C/]
0² 1:2 T
E
[[l
11
<
(T + 1)2 wmll w011
After switching back notation -=
0² and ت gives (14).
9.6 Projected gradient descent with convexity
Theorem 11. Let W be nonempty closed convex set. Let f be a convex function, differentiable W. Let w" argminw (f). ). Let (w*) lee be generated by the Proj-SGD algorithmm, with a constant stepsize 7 E (0, +o0) Suppose that g" is quadratically bounded estimator (Def. 5) for V f at w with constant parameters (a, b, w²). Then.
1 w"
112 +
[((m²)
inf f W
?
S
a
2
1 - AT
1+ 1 and ao def IT- +++1 wt
where def
Finally if
-2 VaT
and T 2 then
ΣT' g++1
2a
E [f(w")] inf inf f
1100
w0l
(15)
V2aT
Thus, T
22 iterations sufficient to guarantee that E [f(w²) inf f] E.
33

---

Proof. Plugging in one step of the algorithm, we have that
w'|
|| ||proi Teg") proiw(w"))
$ 1|w w" - ng' 11
where we used that projections are non-expansive. Taking expectation conditioned on w', and expanding squares we have that
E | w"
S 1| - w²' 1| 11-11-220 -27 (VS(w"), w' - w²) + gp E || wo < || w² 1||- - f(w")) + a - w0 + bn?
where we used that gو is an unbiased estimator with bounded squared norm. Taking full expectation and re-arranging gives
2t E [f(w²) f(w")] s (1 an;2) E [|ret w² w01
[1] -E E ||
w0]]]]
+ by?
From now on telescopic sum. To make this happen, we multiply through by 0+4+1 and impose that 0++1 I C = 0t.
constant stepsize no y. Let C (1 + ar Next we would like to setup
def
This holds with at C.. Multiplying through by 0z41 and summing from / = O,..., T - 1 we have that
Y Σ a+t [f((") - f(ww)] 1 0=
0441 CE ||w0" w 1|]
a++1 E || w1((()) + Σ by²
0
1=
T
04 ||ul fuutww" w11] a+t1 [luttt-wr119) w((] Σ att1
00 00 I w¹11 orr ||10" - -w---w" w²1 + T
=
Using that og = 1, dropping the negative term -or E ||100 11110"-=0110 111 dividing through by
27 DT- 0441 and using Jensen's inequality we have that
E [[f(w') f(ww)] 5
0=0 a+11
E [f(w0) - f(w))]
Eho ax+1
w w0 2
+ by
27 Σk=0
ax+1
2
Finally using that
T-1
aK+1 =
T-1
--(++)) (k+1) =
C
T-1
Σ
C
C
11
1
1
C-1
1
(1+07237
ay
0=
k=0
k=0
gives that
ay²
1w0
w
by
E [](w') [(w")]
24
(1+ay2j7
ary
1/w00
w
by
2
(Itay1)
w w
112
7
+ b
2
[1+0797
| 0 w0
112
+ b
2
4T
34

---

Finally we can apply Lemma 32 with where B A = /2/B= v2/va. 0
1tay and
= A/vT - v2/(VavT) which gives
E [[(())) inf f 2 T
(2all - -0011 + +b)
(2all00 - w² w²1 +bb) =
V2 VavT
10
Solving VI with stochastic optimization
10.1
Convergence of the algorithms
Corollary 12 (Prox-SGD for VI). Consider the VI problem where qu a multivariate Gaussian distribution (Eq. 2) with parameters
(m, C) Rdd x T and assume that this problem admits solution w". Suppose that log log x)) is M-smooth and concave (resp. ju-strongly concave). Generate sequence w' by using the Prox-SGD algorithm (Def. 6) applied to I and h (Eq. 16), using genergy (5) as an estimator of Vl. Let the stepsizes be constant and equal to 1/( denorgyy)) (resp. be decaying in Theorem 7 with Genergy = 2(d + 1+ ). Then, for certain average w of the iterates, we have for 2 2 that
E [[f(w²) inf f] O(L/VT) (resp. E [1| w0 12 O(1/T)).
Proof. Our assumption on the log target imply that / is smooth and convex (resp. p-strongly convex) according to Lemmas 18.1 and 19.1. Furthermore, since log ((. is smooth, know from Theorem that the gradient estimator genergy quadratically bounded at every w with respect to constant parameters (a, b, w",, with a 2(d + 3)M Then, because we work with triangular scale parameters, we know from Lemma 19.3 that h is closed convex. Finally T since
/a, 2) 2
M
M
1.
a
2(d+3)M2
2)
We have now verified the hypotheses of Theorem 8 (resp. Theorem 7), which gives the desired bounds.
Corollary 13 (Proj-SGh for VI). Consider the VI problem where qu multivariate Gaussian distribution (Eq. 2) with parameters w (m, C) E Rd x Sa, and assume that this problem admits solution w". Suppose that log p(-xx) is M-smooth and concave (resp. p-strongly concave). Generate sequence w by using the Proj-SGD algorithm (Def. 9) applied to the function f 1+ (F4. 16) and the constraint Was (Eq. 3), using gent (7) or ESTL (9) as an estimator of V. Let the stepsizes 7/t be constant and equal to /2/(aTT) (resp. be decaying as in Theorem 10) with dent = 4(d + 3)M²
or a = aSTL= = 24(d 3) Then, for a certain average wf of the iterates, we have 2 2 that
E [f(w") inf fl O(I/VT) (resp. E [1|W] w1 O(LTT)).
Proof. Our assumption on the log target imply that is smooth and convex (resp. ju-strongly convex) according to Lemmas 18.1 and 19.1. Regarding the entropy, since pogp((.)) is M-smooth, we know from Lemma 18.3 that h is M-smooth on WM. We also know that is closed convex, since we consider symmetric scale parameters (see Lemma 19.2). Furthermore, the smoothness of log p(-,2) implies that the gradient estimator gent is quadratically bounded at every w² E WWM. with respect constant parameters (a, b, w*). This is stated in Theorem 30, where we can see that for gent we can take Qent = 4(d + 3) M and for gsTL we can take aSTL 2(d + (2K + (2M))². Since we Finally, we know from I Lemma 20 that the minimizers of belong to WM. meaning that argmin([)) argmin (f). We have now verified the hypotheses of Theorem 11 (resp. Theorem 10), which gives us the desired bounds.
know from Theorem that K 2M, we can consider that OSTI. S 24(d + 3)M2.
D
35

---

Corollary 14 (Proj-SGD for VI Gaussian target), Consider the setting of Corollary 13, in the scenario that log 2() p-strongly concave, that we use the ESTI. estimator, and that we take constant stepsize " = 7 E (0, min Zagri. 33. Assume further that p(-[z) Gaussian. Then,
1 T 1 - w=
[[[W" w
- (1-27
Proof. We use the same arguments as in the proof of Corollary 13, but this time we use a constant stepsize, which gives the bound (13) from Theorem 10, where b 4(d I w/| according to Theorem 4. All we need conclude is to verify that b = O. This comes from the
assumption that the target is Gaussian, from which we deduce that K 0 (see again Theorem 4). D
10.2 Computations for the algorithms
Lemma 35 (Prox of the entropy for triangular matrices). Assume that V T and let w (m, C) Id Tá. Then prox, (w) (m, C), where
[C +4g+ Ci,
if if H
Cig
Proof. This is essentially proved in [13, Theorem 13]. We provide proof here for the sake of completeness, also because our entropy here enforces positive definiteness. So we need to compute
argmin =ERd,XER X
log det X + 27 || ml|
1 - -
-
Because the objective here is separated in 2 and X, we can optimizae those variables separately. is clear that the optimal will be m. For the matrix component, it remains to solve
argmin log det X XeT
27 ||x
CIP.
Now remember that for a triangular matrix, being positive definite is equivalent to have positive diagonal elements (Lemma 19.3). Because X is lower triangular matrix, have log det X Σd - - log X5 and Xiş when <.. So our problem becomes
argmin Σ (- log X.:) + Σ 1 (X,, - C.y)- XETT,Xwca i=1
27
123
This is gain a separated optimization problem with respect to the variables Xij. So we deduce that when 2 j must take Cij Cij. and when j need take
C2 = argmin - log + (x c
(41)
27
#>0
This is the proximal operator of 11-logg logg which can be classically computed [55, Example 6.9], and gives the desired result.
Lemma 36 (Projection onto WM for symmetric matrices). Assume that V = Sd, and let (m, C) TRd x Sd. Let C UDUT be the eigenvalue decomposition of C, with D diagonal and U orthogonal. Then proiws, (w) = (m, U DÚT), where δ the diagonal matrix defined by
DR Dii
VM
Proof. Consider g = 6war to be the indicator function of Was which is equal to 0 on WMr and +oo elsewhere. Then projwar prox p. and we can compute it with [55, Theorem 7.18). For this, we need to write gCC) F(X(C)), where X(C') is the set of eigenvalues of C, and F is the indicator function of K = {a E R z; 2 1/VM). In that way, we see that g symmetric spectral function in the sense of [55, Definition 7.12]. and the conclusion follows.
36

---

Lemma 37 (Adjoint of the distribution transformation). Let u E R be fixed. Let T: R x V - R be such that T T'm. C) CutCun m. Then T is linear operator such that its adjoint T R R XV verifies T (z) (z, projy(zu)))
Proof. The fact that is linear is immediate. Let us consider now the extension T Rd M - Rd such that Tm,, C*) Cuu + If define the canonical injection : R2 x V - R M clearly have T = t. This means that = o T = projge T It remains to prove that (=) = (z, zu ). For this, it suffices to use the properties of the trace to write that
((2, zu 2, (m, ())
(z, m) + (zu ,C) = m2 + tr(uz C)
(z, m) + tr(=TCu) (z, m) + (z, Cu)
(2, Tm,, C)).
10.3 Explicit implementation of the algorithms
Lemma 38 (Explicit computation of the estimators). Let E IR4, and let w (m, C) E R2 x V where 0. Define
V= "'loop(Cuu m,
Then
1. Benergy (u) = (Λ., (rojy(TuT]]
2. gent (u) (T,, proiv ( C--T),,
3. ESTL (u) = (IT C-Ju, projy (Tu - C-TTu
Proof. Let W (m, C) and let T the linear operator defined by T(m, C) = Cu + m, as in Lemma 37. As shorthand we will note pz := p(;,).
1. From its definition in (5), the energy estimator is equal to
genergy (u) = V m (log o Pa. o (w))
Applying the chain rules, we obtain
genergy (u) V & (log P2) (Tw) T n.
The conclusion follows from the expression of T obtained in Lemma 37.
2. From its definition in (7), the entropy estimator gent (u) is equal to geuergy (u) Vh(w), where (we use the fact that we assumed that w W) we know from Lemma 18.2 that Vh(w) = The conclusion follows from the previous item and the lincarity of projy.
-00, (Cojv(C-T))
3. From its definition in (9), the STL estimator equal to
gSTL. (u) genergy (u) V o (w:)
where (x2) log (2), with v being a copy of w. From the definition of Gaussian density (17) we have that
<=((2) = -0 -C -C--TC-) mm),
so we deduce that
(w) --C
TC-
m) -T
The conclusion follows after combining the first item with Lemma 37.
Lemma 39 (Pscudocode of the algorithms).
37

---

1. Algorithm 1 is equivalent to the Prox-SGD algorithu (Def. 6) applied to and (Eq. 16), using genergy (5) as an estimator of VL and using lower triangular covariance factors.
2. Algorithm 2 is equivalent the Proj-SGD algorithm (Def. 9) applied to the function 1 (Eq. 16) and the constraint Wyr (Eq. 3), and using gent (7) as an estimator of f.
3. Algorithm is equivalent the Proj-SGD algorithm (Def. 9) applied to the function f
(Eq. 16) and the constraint WM (Eq. 3), and using gSTL (9) as estimator of f.
Proof. In this proof, we fix an iteration E N, consider ut N I) and we note w (m+, Ct).
1. From (Def. 6) we know that w = prox, y (w" - ~Yteenergy (u²)). According to Lemma 38, we have
w - 7/genorgy (1t) = (m,, Ct) - - projT= )),
where
Tt - V log p( (C me, 2)). Moreover, we know from Lemma 35 what prox is:
[2 + VC +4%)
if i j.
pr C? (m, Ry (C)), where Ry.
Cip
First, that with respect to the variable the proximal operator is the identity, meaning that we indeed have
M4+1 my T=
Second, we see that with respect to the C variable we have
C++1= R, (C, "/ projoe )).
It is immediate to see by induction that since we initialize with C. E Td, then C E T well. Therefore we can use the linearity of projy to write that
C++1 = Rr, (proj - -)eteu)))
The conclusion follows after observing that
C(i
ifi > j. if = j.
(Ry, o projT.) (C) =
1
((1,)) VC(, +4))
2
if < j.
2. From (Def. 9) we know that w*+1 projww (w* Mgent (ue)). According Lemma 38, have -² "/gee (wx) = (mi, C.) - proois((Te)) CTT))
Morcover, we know from Lemma 36 what projww
projwz (m, C) (m, R())))
where R(C) = UDUT, with UDU being SVD decomposition C, and D(i max i); 1/13 First, we see that with respect to the m variable the projection is the identity, meaning that indeed have
m441 = m, 7e.
Second, we see that with respect to the C variable we have
C++1 R(C, - 7 p CT )).
It is immediate to see by induction that since we initialize with Co Sd, then C, E S as well. Therefore we can use the linearity of projsa to write that
C++1 R (prois= (C. - - - CT ()))
The conclusion follows after setting Ō.+1 = =C, - - - CT TTT), and using the fact that (C) (C+ CT
3. Use the same reasoning as in the previous item, with the appropriate computation for the STL estimator from Lemma 38.
38

---

Algorithm I Prox-SGD with energy estimator and triangular factors
Let mg E IRd, Co E T such that Co > 0, (TH)EEN C (0, +00)
for i g N do
Sample dh ~ (),
Compute गH - log (CC, m)
Set m++1 = My Tt
Set C4+1 = Ci - 7
Co+1(1,
ifi> j ifi=j
Set Ce+1 via (t,5( =
(C+i((,0)) + V Ci+i 47H
2
0
ifi
Algorithm Proj-SGD with entropy estimator and symmetric factors
Let mg E R2 Co Sd such that Co O, (Te)xeN C +20)
for e N do
Sample u, No,
Compute
- V log (CC + m4, x)
Set m++1 my 7t
Set C++1 C. - 7/1 (Tp - 1)
Compute U241 D0+1 U,T a singular value decomposition (Č(+1 + 112/
Set C++1 Ul+1 D2+1U 4+1 where i max (i, M 1/23
Algorithm 3 Proj-SGD with STL estimator and symmetrie factors
Let mg E RRd Co E Sd such that Co O, (Te)ee C (0, +20)
for E N do
Sample Wt (), - log [CC, + m4 2)
Compute
Set m++1 = my - 7m( C, uz)
Compute -C.-2i (,,4, - ---14
Set Ōc+1 = C. C.-
441 singular value decomposition of (Č4+1 + 2T22
Set C4+1 = Un+1 Dz41 UTT where (i+)(i, = max(Di41ii.)) 1/23
39


Reviews:
[{"review_id": "peUss9lUB0", "reviewer": "Reviewer_TUBr", "summary": "This paper analyzed the convergence of Black-box VI, which has been widely used in variational inference in recent years. Under the assumption that the target joint distribution is convex or strongly convex and the variational posterior distribution is Gaussian, the convergence of the variational parameter obtained by Black-box VI to the ELBO optimal solution is first presented when using Prox-GD and Projected Gradient Descent. An important aspect of the analysis of Black-box VI is that the gradient is stochastic.\nTo deal with this, the authors introduced a class of quadratically bounded estimators as estimators of the gradient and combined them appropriately with the existing Prox-GD analysis.\n\n", "questions": "I would like to know more about the novelty of this study in the context of Stochastic Proximal Gradient Descent in terms of technical aspects and assumptions of the proof.", "limitations": "Yes, the limitations of the statements are clearly explained.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Although the proposed analysis makes a number of assumptions, the results are significant in that they include an important Bayesian model that is used in practice.\n\nFor the first time, an analysis of the convergence of Black-box VI is presented, albeit with a rather limited problem, such as the assumption of variational posterior distribution as Gaussian distribution.", "weaknesses": "It was very difficult to understand how novel and important the proposed analysis is from a single reading. In existing studies, each component of the presented analysis is already well known, but I am not sure how important and novel the presented  VI analysis is as the combination of those well-known results.\n\nBy assumption, since the objective function is convex (or strongly convex), the convergence using Prox SGD or other methods seems apparent, the difficulty of the problem stems from the noisy version of the gradient, as discussed in  Sec 2.3. On the other hand, as there have been many studies discussing noisy versions of prox GD, such as Stochastic Proximal Gradient Descent[1], so I am not sure how novel this analysis is.\n\nIn the proofs of Theorem 6 and 7, I understood that the novelty is mainly the technique to treat the gradient estimator that satisfies the quadratic bounded property shown in Definition 4. However, I still did not understand how important it is compared to existing studies such as the Stochastic Proximal Gradient Descent study.\n\n[1]Convergence of Stochastic Proximal Gradient Algorithm, L. Rosasco et al.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review.\n\n### Novelty of the analysis\n\nGiven the current state of knowledge, the convergence of Prox/Proj-SGD is not apparent, because the convexity of the objective function alone is not enough. A second essential hypothesis must be satisfied by the estimator of the gradient: typically, results require a noise bound that depends on the suboptimality gap, gradient norms and constants. This is for instance the case in the paper [1] suggested by the reviewer, and we also refer to [KR] for a good discussion on noise bounds that result from standard assumptions. Those standard noise bounds can be verified under the common assumption that the stochastic functions are (uniformly) smooth, or Lipschitz. But our problem is that our gradient estimators do not satisfy such a typical noise bound, but instead the quadratic error bound presented in the paper. There is indeed a significant literature on the analysis of Prox/Proj-SGD, but we carried out a thorough bibliographic review and consulted with experts to find that no existing analysis could accommodate such quadratic noise bound while at the same time handling a composite non-smooth objective. The goal of Section 2.3 in our paper is to explain to the reader exactly what in our paper is novel. We will try to explain this more clearly in the revision. \n\n[KR] Better Theory for SGD in the Nonconvex World, Ahmed Khaled & Peter Richtarik, TMLR 2022\n\n### Is the main novelty in treatment of the gradient estimators?\n\nBroadly speaking, yes. Our proof technique borrows several ingredients from previous analyses of SGD, and the way we handle convexity and smoothness is rather classical. The main difficulty we faced is how to handle the quadratic noise bound together with composite non-smooth objective.. For this we needed a relatively new weighted telescoping technique, the application of which involves substantial technical difficulties (see proofs on pages 21-30 of appendix).\n"}, {"review_id": "x2xGj6Iy1F", "reviewer": "Reviewer_6P4Y", "summary": "The paper addresses the lack of provable convergence guarantees for black-box variational inference (VI) and proposes convergence guarantees for two stochastic optimization algorithms applied to Gaussian variational families. \n\nThe authors identify challenges in analyzing VI as a standard stochastic optimization problem, including the non-smoothness of the objective function, lack of uniform smoothness, and lack of uniformly bounded noise for gradient estimators. \n\nThey provide theoretical results and noise bounds for two gradient estimators and propose proximal and projected stochastic gradient descent algorithms.\n\nI did not find technical flaws.", "questions": "See weakness.", "limitations": "None.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper addresses an important problem by providing provable convergence guarantees for black-box VI, which is widely used but lacks theoretical guarantees.\n\n2. The authors clearly articulate the challenges in analyzing VI as a standard stochastic optimization problem, such as the non-smoothness of the objective function and the lack of uniform smoothness and noise bounds. The paper provides rigorous theoretical results and noise bounds for gradient estimators used in black-box VI.\n\n3. The paper is well written. It provides a comprehensive and coherent review of the literature and clearly presents the technical contributions and proof ideas. I'm not a theory guy but I enjoy reading the paper.", "weaknesses": "1. I wonder if is it possible to perform some simple (even toy) examples to empirically check the results and illustrate the theory in a better way. \n\n2. The dense Gaussian family is still limited (at least in my opinion). It is better to discuss whether the techniques developed in the paper can be applied to more general cases (e.g., black-box variational inference with a neural network model) or at least discuss the challenges. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review. We will respond to two points.\n\n### Examples\n\nWe ask for some consideration for the constraints imposed by a 9-page limit. Given that this is a theoretical paper whose goal is to provide guarantees for algorithms already commonly used in practice, we prefer to focus entirely on theory.\n\n### Dense Gaussians\n\nWe agree this is an important issue that should be discussed more in the paper. We emphasize that our new optimization results in section 3 apply more broadly than just Gaussian VI (and indeed apply more broadly than VI)—they apply to any composite objective consisting of a smooth part and a non-smooth part where (1) Either (1a) The non-smooth part is proximalable or (1b) can be projected onto a constraint set where it is smooth and (2) the gradient noise obeys a quadratic bound. It is conceivable that such guarantees could be shown for more general variational families, and if so our convergence results would provide “plug in” optimization guarantees. We suspect it might be fairly easy to establish such guarantees for, e.g., elliptical distributions or location-scale families, using similar proof strategies as for Gaussians. But doing it for much broader classes would definitely be a challenge—the properties (1) and (2) we use above were developed over several years in a sequence of papers. We propose to add a short discussion of this topic to section 5 of the paper.\n"}, {"review_id": "e06OEzY5c5", "reviewer": "Reviewer_biC8", "summary": "The paper proves 1/sqrt(T) respectively 1/T convergence rates for black-box variational inference methods, when implemented with a proximal stochastic gradient method.  Such rates were not available in the literature until now, due to difficulties in bounding the gradient noise. The main contribution of the paper is to prove new bounds on the gradient noise and incorporate them into the proof techniques for proximal SGD. ", "questions": "1) Which algorithm would be preferable in practice, prox-SGD or proj-SGD? From the theory side, it seems they have similar convergence rates, but maybe there are other things to take into consideration?\n", "limitations": "All limitations are addressed.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well-written and clearly states its contributions.  While the topic seems rather narrow, the technical part of the paper appears to be a novel and non-trivial result. ", "weaknesses": "1) The paper solves a quite \"narrow\" problem, that is fixing some gaps in existing convergence proofs for variational inference. It is unclear whether this work will be of interest to the wider ICML community, as it is a niche topic, and the theory does not seem to directly suggest any practical improvements or changes. \n\n2) BBVI is sometimes less preferable in practice than natural-gradient algorithms. These methods use the KL / Fisher-geometry (https://arxiv.org/abs/2107.04562). This line of works could be mentioned in the introduction.\n\n3) There has been some recent works\nhttps://arxiv.org/abs/2205.15902\nhttps://proceedings.mlr.press/v202/diao23a.html\nwhich prove convergence for BBVI-like algorithms by interpreting them as Wasserstein gradient flow. Perhaps these works could be mentioned. \n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review.\n\n### Natural gradient algorithms\n\nWe agree this should be discussed. We will also mention that extending theory to address these algorithms remains an open problem.\n\n### Related work by Diao et al.\n\nWe agree we should mention these recent works. (We note that they make use of Hessians which can be computationally challenging and is distinct from most practical BBVI algorithms.)\n\n### Prefer prox-SGD or proj-SGD\n\nThe reviewer is correct when saying that from a theory side, prox-SGD and proj-SGD have similar convergence rates with respect to T. If we want to compare them further, we need to look closely at the multiplicative constants in the rates, which depend essentially on the constants (a,b) in the quadratic bound. Since prox-SGD needs the “energy” estimator and proj-SGD needs the “entropy” estimator, the question could be shifted to: which estimator has the best constants (for instance the lowest variance). From that perspective, the two estimators we present have roughly the same constants.\n\nAnother consideration in practice is the complexity of the projection vs proximal steps. In many cases, these will both be very cheap but in certain high-dimensional settings, the complexity of the proximal step could be higher (Θ(d^3) vs Θ(d^2), see lines 273-281). \n\nThis being said, the story could be completely different when considering *other* estimators. Intuitively, we could imagine that an estimator will be able to have lower variance if it is a stochastic estimator of the full objective function instead of just the energy function. That scenario would give an advantage to the proj-SGD algorithm. We have just now been able to prove exactly this for the “Sticking-The-Landing (STL)” estimator. We showed that it is quadratically bounded, and further that its constant “b” is directly related to how close to a Gaussian the target is. This formalizes the idea that the easier the problem is (here approximating an almost Gaussian with a Gaussian), the lower the variance of the estimator is. In the extreme case that the target is exactly a Gaussian, this variance “b” is exactly zero, meaning that the rates drastically improve from 1/T to exponential rates. In this setting, proj-SGD would be a better choice.\nWe propose to add this new small result in the paper, and to develop the discussion on the cost at the end of Section 4. We will underline the fact that the variance of the chosen estimator has a strong impact on the complexity (a well-known fact in stochastic optimization) and that some estimators can have a significantly smaller variance (like STL for easy problems)."}, {"review_id": "bUeBuKxnWF", "reviewer": "Reviewer_gXey", "summary": "This paper offers a convergence proof for the stochastic optimization problem inherent in full-rank Gaussian variational inference when the log-density of the target is concave. The primary challenge of the convergence proof lies in managing the non-smoothness present in the entropy term of Gaussian VI. This issue is addressed in the current work by considering proximal and projected stochastic gradient descent. To evaluate the validity of the imposed assumptions, case studies involving Bayesian (generalized) linear regression problems are provided.", "questions": "1. Compare to  [Xu&Campbell, 2022]: As I mentioned earlier, this is a very relavant past work on this topic. I wonder specifically, if we ignore all the asymptotic stuff and assumes $\\log p$ is strongly-log-concave, does [Xu&Campbell, 2022] achieve similar convergence rate ($O(1/T)$) as the present paper? \n\n2. I wonder if the proximal operator (line 264) has been noted in past VI literature? Or it is proposed by the author.\n\n", "limitations": "See Weakness", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "I appreciate the clarity of the writing of the present work; it clearly descibes the scope of the problem studied and states the main challenge of showing the convergence. The layout of each section forms a smooth flow of the proof strategies, in addition to a nice organization of the mathematical proofs, which makes the reading enjoyable.\n\nThe case study provided in appendix (sec 7.3) is helpful for reasoning imposed assaumptions. \n\nThe proximal operator introduced for optimizing the covariance matrix is novel to me. Although it's a standard techniques in optimization literature---splitting non-strongly smooth (but closed convex) part of the objective function using proximal operator, it's the first time I see in Gaussian VI literature. (This could of course be due to my lack of knowledge in this field) Additionally, the weighted telescope summation used when $-logp$ is only convex(rather than strongly convex) is very interesting, which is a setting that rarely considered in literature as far as I'm aware.   \n\nEven though the paper limits the scope to dense Gaussian VI problem, I think the present techniques can be applied to general location-scale variational families (as long as the variance quadratic bounds still holds).\n\n  ", "weaknesses": "In my view, the principal limitation of this work is the somewhat restricted scope of the Variational Inference (VI) problem it examines. Specifically, the variational family under consideration is: (1.) Gaussian (albeit full-rank), and (2.) the target distribution is log-concave or strongly-log-concave. Additionally, (3.) data subsampling on $\\log p$ is not taken into account. The proof techniques necessitated by this setting, in my opinion, are fairly standard within the stochastic optimization literature. While I acknowledge that any relaxation of condition (2) would likely preclude anything beyond convergence to some stationary point, and that the quadratic variance bound is heavily reliant on condition (1) (which is probably not a major concern for the broader community in regards to more complex VI families), offering a convergence result when data subsampling is implemented could greatly increase the impact of this work. I would be inclined to raise the score to 8 if data subsampling was considered, or if existing proof techniques could easily address it (though I doubt the current approach is effective in this case, please correct me if I'm mistaken).\n\nRegarding originality, even though I concur that explicit results may be lacking for this specific problem setting, I don't think this is the **first** optimization theory guarantee for full-rank Gaussian VI. For instance, [Xu & Campbell, 2022] studies the convergence of full-rank Gaussian VI without the assumption of a log-concave target (even though they utilize posterior asymptotics to somewhat reduce the underlying optimization problem on a strongly-log-concave target, they also employ a special scaling operator (Eq(8)) to handle the non-Lip-smoothness of the entropy term $\\log \\text{det}$). I recommend that the authors integrate an optimization analysis of the scaled stochastic gradient descent proposed in [Xu & Campbell, 2022], and carefully compare the results to those provided in [Xu & Campbell, 2022] within a strongly convex setting at least (ignoring the data asymptotics and merely assuming that $\\log p$ is strongly-log-concave should align their work with the current setting).\n\nFinally,  I encourage the authors to provide additional case studies, perhaps involving more Bayesian GLMs or even Bayesian sparse regression with a Horseshoe prior. From what I understand, most of these models yield a log-concave target, but the application of the variance quadratic bound is unclear to me. If negative results were to appear, for example, some models not having a quadratically bounded gradient noise, these examples would still significantly benefit the community.\n\n\nAt the end, I hope to stress that despite the aforementioned weakness/limitations, I think this is a really good work for this niche research direction.  \n\nReference:\n [Xu&Campbell, 2022]: The computational asymptotics of Gaussian variational inference and the Laplace approximation, Statistics and Computing", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review.\n\n### Convergence result when data subsampling \n\nOur proofs can indeed address subsampling with minor technical difficulty, namely bounding the variance of a slightly changed gradient estimator. Note that the main optimization results only depend on (1) the structural properties of the problem (which are unchanged with subsampling) and (2) the constants (a,b) bounding the gradient variance. Our gradient variance guarantees in Theorems 2 and 3 are all based on Theorem 1 from [20]. However, that reference also provides a more general version of Theorem 1 (Theorem 6 in [20]) which considers data subsampling. Very roughly speaking, with uniform data subsampling of 1 datum at a time, the expected squared norm of the gradient increases by a factor of between $1$ (no increase) and $ndata$, depending on how correlated the data are. (The less correlation, the larger the increase.) So, with subsampling, both Theorem 2 and Theorem 3 would increase by a factor of $1$ to $ndata$. (To be more precise, the second term in theorem 3 would be unchanged, though this would typically be dominated by the first term.) Those increases would then manifest as an increase of between $1$ and $ndata$ for both a and b in the quadratic noise bounds. After that point, exactly the same convergence results hold. We will add a discussion of this.\n\n### Compare to [Xu&Campbell]\n\nThank you for pointing out this paper which we were not aware of. It does indeed provide a convergence guarantee for an algorithm similar to proximal-SGD using the g_{energy} estimator, and it was an oversight not to cite and discuss it. The principal differences with our work is (1) their guarantee is asymptotic (holds in the limit of large T) (2) their guarantee is local (holds when started close to the solution) (3) our analysis also considers the g_{ent} estimator, and (4) our analysis considers the convex and strongly-convex cases, whereas theirs considers the non-convex case . Their analysis does not consider the strongly convex (equivalently when p is strongly-log-concave) assumption so only gives a O(1/√T) rate [Thm 3 in their paper] rather than O(1/T). We will discuss this paper and revise our claims appropriately.\n\n### Prox operator\n\nThe proximal operator has indeed been used in the past in VI [13] which we should note more clearly."}, {"review_id": "ZHwpFdu5Vq", "reviewer": "Reviewer_vBBM", "summary": "This paper offers the first convergence results for black-box VI, a widely used and popular framework for Bayesian problems. Under assumptions on the log model, $\\log p$, and given a Gaussian variational family of distributions, convergence rates are established by utilizing recent advances in the field. The assumptions are motivated by giving practical examples of when they are satisfied. ", "questions": "* In Sec. 2.1 it stated that $h$ is $M$-smooth over the set $\\mathcal{W}_M$. Based on this, do you make any assumptions on $q_w$ other than it being Gaussian? I was expecting some new assumption given the title of the section, but maybe this is a \"structural property\"?\n* Out of curiosity, could you expand on the constructed problem where the variance of the path-based estimator is \"vastly\" increased and the score-based estimator is not? Where is the noise added? To $\\phi$?  ", "limitations": "The assumptions are clearly highlighted and discussed in the paper. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "This work seems to lay the missing puzzle piece in a series of works towards establishing convergence rates for black-box VI. As such, it is clearly original, and it is significant as black-box VI is widely used in practice and is a popular research topic. Although the results build on seemingly strong assumptions on the log model, $\\log p$, the significance of the results is emphasized by exemplifying broad classes of problem settings where the assumptions are satisfied, which deserves credit. I believe this paper will be of interest to the NeurIPS community.\n\nThe paper follows a clear line of argument and builds up to its results in a pedagogical manner, giving a convincing account of the existing, related works.", "weaknesses": "**W1**: $C$ is a matrix (line 63), but in line 120 I read it as $C = 0$, by inspecting $\\bar{w} = (\\bar{m}, 0)$. It is not clear to me what this means. Are all elements in the covariance matrix zero (since $\\Sigma = C C^T$)?\n\n**W2**: In Sec. 2.1 I read it as, if $\\log p$ is $M$-smooth then we can expect that the optimal covariance parameter for $q_w$ should not be too small. This makes sense to me. In fact they should be at least $1 / \\sqrt{M}$ as shown in Domke's previous work. However in Theorem 2, where the log-target is assumed to be $M$-smooth, it says that $\\bar{w}=(\\bar{m}, 0)$ is the \"maximum of $\\log p$\". I am confused about this. Doesn't $\\bar{w}$ here imply that $\\log p$ is not smooth, and that the optimal parameters $w^* = \\bar{w}$, i.e. $C^* = 0 < 1 / \\sqrt{M}$?\n\nIf I am helped clear up this confusion, I am willing to raise my score.\n\n**W3:** I am not sure how the discussion about the score-type estimator in Eq. 4 contributes to the paper. The results in paper are based on the reparameterization trick (i.e. the path-type estimator in Eq. 5), no? What is the contribution of this discussion to the results in the paper?\n\n**Minor issues**\n\nIn line 16, $h$ is called the entropy, but as it is defined in Eq. 1 and 15 it is the neg-entropy. \n\nSometimes (e.g. in the Central contributions box and on line 93) $q_w$ is referred to a family of distributions and sometimes (line 243) it is used a density. In VI the \"family\" is often meant as the set of distributions in which we want to find the optimal approximation, no? So can it really be a family and a distribution simultaneously? Perhaps it is more clear if $\\mathcal{Q}_w$ is used to notate the family of Gaussian distributions parameterized by $w$?\n\nTypo: \n* line 35: objective is misspelled\n* line 141: \"that bar\" seems to be a typo", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### C is a matrix but C=0\n\nYes, in this context, of writing something like w=(m,0) then 0 means a matrix of zeros. If φ(z)=-log p(z,x), then the gradient noise is bounded in terms of how far the parameters w=(m,C) are from representing a delta function centered at the MAP solution. (Note this is different from the optimal solution w*.)\n\n### $\\bar{w}$ is the maximum of $\\log p$, doesn't this imply log p is not smooth?\n\nIt’s important to note here that $\\bar{w}$ is not the optimal solution to the problem, which we denote by $w^*$. We do not say that $\\bar{w}$ is the maximum of $\\log p$ (which is not true) but rather that $\\bar{w} = (\\bar m, 0)$, where $\\bar m$ is the maximum of $\\log p$. Intuitively, we can see $\\bar{w}$ as parameters that represent a delta function (because it has zero covariance) centered at the MAP solution. Note that this is just an intuition and we never make use of this interpretation. We propose to rewrite the Theorem to make this more clear.\n\n### Discussion of score-type estimators\n\nWe appreciate this point. Arguably this discussion is unnecessary and distracting at the point it currently happens. Our reason for including it was that we wanted to mention that there are situations where score-function type estimators have lower variance than reparameterization, and so these arguably deserve further investigation. We believe this comment would make more sense in the discussion section.\n\n### Any other assumptions needed for h to be M-smooth?\n\nThe assumption that h is M-smooth (over $W_M$) requires no other assumptions : it is indeed a structural property. We will clarify how this is discussed.\n\n### Expand on the problem with variance of path-based is increased and score-based is not?\n\nIn terms of how the variance of the path-based estimator could be increased, imagine some problem with an unnormalized posterior φ(z) which is equal (up to a constant) to -log p(z,x). Suppose for simplicity that z is just a scalar. Now imagine changing to a new unnormalized posterior φ’(z) = φ(z) + √ε sin(z/ε) for some very small value ε. The difference of φ and φ’ is trivial (they only vary by a tiny amount). However, the derivative of the added term is cos(z/ε) / √ε. When ε is very small, this derivative will be huge. This will vastly increase the noise of a reparameterization estimator but have little effect on a score-based estimator. (In a sense, this is why the assumption that log p(z,x) is M-smooth is so important—it prohibits posteriors that do things like this!)\n\n### errors / typos\n\nThank you for pointing these out. We will correct them. We agree that the way we used the terminology of “density” vs “family” is confusing. We will fix this by saying that q_w for a specific w is a Gaussian distribution, whereas the set {q_w | w ∈ W} (or Q_w)  is the family."}]

Author rebuttal:

Meta-review: This paper proved $1/\sqrt{T}$ and $1/T$ convergence rate for black-box VI with convex and strongly convex model log likelihood, respectively, under Gaussian variational distribution assumption and some Lipschitz conditions for the model. The main techniques are to provide control for the noise in path-based BBVI gradient estimators as well as to use theoretical results in stochastic proximal gradient descent in some smart way. 

Reviewers commend that the paper, although with assumptions that narrow the scope, makes a significant progress towards addressing the important theoretical challenge of convergence guarantees for BBVI. They provided comments regarding discussions with existing theoretical results analysing Gaussian VI as well as another very recent work [Xu & Campbell, 2022] related to BBVI theoretical studies. One reviewer also pointed out a bug in one of the lemma, which is then fixed by the authors without impacting on the main theoretical conclusions.

One reviewer mentioned potentially adding experiments to verify the theory. I personally think that experiments are nice to have, but given the significant contribution on the theory front, having no experiment shouldn't be a reason for rejection in this case.

I suggest the authors to include discussions regarding other theoretical papers for Gaussian VI and [Xu & Campbell 2022] in revision.
Decision: Accept (poster)