Exact manifold Gaussian Variational Bayes
Authors: Martin Magris
Mostafa Shabani
Alexandros Iosifidis
Venue: Submitted to ICLR 2023 (2023)
Keywords: variational inference
Bayes
Riemann
black box
deep learning

Abstract: We propose an optimization algorithm for Variational Inference (VI) in complex models. Our approach relies on natural gradient updates where the variational space is a Riemann manifold. We develop an efficient algorithm for Gaussian Variational Inference that implicitly satisfies the positive definite constraint on the variational covariance matrix. Our Exact manifold Gaussian Variational Bayes (EMGVB) provides exact but simple update rules and is straightforward to implement. Due to its black-box nature, EMGVB stands as a ready-to-use solution for VI in complex models. Over five datasets, we empirically validate our feasible approach on different statistical, econometric, and deep learning models, discussing its performance with respect to baseline methods.

EXACT MANIFOLD GAUSSIAN VARIATIONAL BAYES
Anonymous authors
Paper under double-blind review
ABSTRACT
We propose an optimization algorithm for Variational Inference (VI) in complex models. Our approach relies on natural gradient updates where the variational space is Riemann manifold. We develop an efficient algorithm for Gaussian Variational Inference that implicitly satisfies the positive definite constraint on the variational covariance matrix. Our Exact manifold Gaussian Variational Bayes (EMGVB) provides exact but simple update rules and is straightforward to imple- ment. Due to its black-box nature, EMGVB stands VI in complex models. Over five datasets, we empirically validate our feasible ap- proach on different statistical and econometric models, discussing its performance with respect to baseline methods.
ready-to-use solution for
INTRODUCTION
Although Bayesian principles are not new to Machine Learning (ML) (e.g. Mackay
1992
1995
Lampinen & Vehtari 2001), it has been only recently that feasible methods boosted a growing use of Bayesian methods within the field (e.g. Zhang et al. 2018 Trusheim et al. 2018 Osawa et al.) 2019 Khan et al. 2018b Khan Nielsen 2018). In the typical ML settings the applicability of sampling methods for the challenging computation of the posterior is prohibitive, however approx- imate methods such as Variational Inference (VT) have been proved suitable and successful (Saul
et al. 1996; Wainwright & Jordan| 2008 Hoffman et al., 2013) Blei et al.) (017). VI is generally
performed with Stochastic Gradient Descent (SGD) methods (Robbins
Monro 1951 Hoffman
ct al.
2013; Salimans & Knowles, 2014), boosted by the use of natural gradients (Hoffman et al.)
2013
Wierstra et al. [0144] Khan et al. (0188)), and the updates often take a simple form (Khan &
Nielsen 2018 Osawa et al. Magriss et al. 2022).
The majority of VI algorithms rely on the extensive use of models' gradients and the form of the variational posterior implies additional model-specific derivations that are not easy to adapt general, plug-and-play optimizer. Black box methods (Ranganath et al., 2014), are straightforward implement and versatile use as they avoid model-specific derivations by relying on stochastic sampling (Salimans & Knowles 2014 Paisley et al. | Kingma & Welling 2013). The increased variance in the gradient estimates as opposed e.g. methods relying on the Reparametrization Trick (RT) (Blundell et al. 2015 Xu et al. 2019) can be alleviated with variance reduction techniques (e.g
Magris al. 2022).
Furthermore, the majority of existing algorithms do not directly address parameters' constraints. Under the typical Gaussian variational assumption, granting positive-deffnitenes of the covariance
matrix is an acknowledged problem (e.g TTan et al. 2021a] Khan et al.) 2018b] Lin et al. 2020).
Only a few algorithms directly tackle the problem (Osawa et al. 2019 Lin et al. 2020). see Section A recent approximate approach based on manifold optimization provided by Tran et al. (202
3
On the theoretical results Khan & Lin (017)) : Khan et al. (018a) we develop an exact version Tran et al.) (2021a), resulting in algorithm that explicitly tackles the positive-deffiiteness
straint for the variational covariance matrix and resembles the readily-applicable natural-gres -attrel-gradient black-box framework of (Magris al. 2022]. For its implementation, we discuss recommendations and practicalities, show that EMGVB is of simple implementation, and demonstrate its feasibility in extensive experiments
four datasets, 12 models, and three competing VI optimizers.
In Section 2 we review the basis of VI. in Section we review the Manifold Gaussian Variational Bayes approach and other related works, in Section 4 we discuss proposed approach. Experi-

---

while Section concludes. Appendices A HBB complement the main discussion, Appendix reinforces and expands the experiments, Appendix provides proofs.
ments are found in Section 5
2 VARIATIONAL INFERENCE
Variational Inference (VI) stands convenient and feasible approximate method for Bayesian inference. Let denote the data, p(y(0) the likelihood of the data based on some model whose k-dimensional parameter is 0. Let p(0) be the prior distribution ence the posterior is retrieved via the Bayes theorem as p(Oly) p( p())p()) ())((y) As the marginal likelihood p(y) generally intractable, Bayesian inference is often difficult for complex models. Though the problem can be tackled with sampling methods, Monte Carlo techniques, although non- parametric and asymptotically exact may be slow, especially in nigh-dimensional applications (Sal-) imans et al. 2015).
0. In standard Bayesian infer-
VI approximates the true unknown posterior with a probability density within a tractable class of distributions Q, such as the exponential family. VI turns the Bayesian inference problem into that of finding the best variational distribution q" Q minimizing the Kullback-Leibhe from to p(Oly): 9º arg min DKL (g)l It can be shown that the KL minimization problem is equivalent to the maximization of the so-called Lower Bound (LB) on log p(y). (e.g. Tran al., 2021b). fixed-form variational Bayes, parametric form of the variational posterior is set. The optimization problem accounts for finding the optimal variational parameter parametrizing q = qg that maximizes the LB (C), that is:
(KL.) divergence
5* = arg max (()) = qc ((((0) (o) p(0)p(y10) (0)
p(0)p(y10)] (c(()
dB EEg
log
CEZ
where Ey means that the expectation is taken with respect to the distribution ạç, and Z is the param- eter space for Ç.
The maximization of the LB is generally tackled with a gradient-descent method such as SGD (Robbins & Monro 1951 ADAM Kingma Ba (2014), ADAGRAD Duchi al. (001I) The learning of the parameter based on standard gradient descent is however problematic as it ignores the information geometry the distribution qc, is not scale invariant, unstable, and very susceptible to the initial values (Wierstra et al.) [014). SGD implicitly relies on the Euclidean norm for cap- turing the dissimilarity between two distributions, which can be a poor and misleading measure of dissimilarity (Khan & Nielsen 2018). By using the KI. divergence in place of the Fuclidean norm, the SGD update results in the following natural gradient update:
= c + Be [V<[(S)]
where B, is possibly adaptive learning rate, and denotes the iteration. The above update results in improved steps towards the maximum of the LB when optimizing it for the variational parameter Ç- VcL(C) is obtained by rescaling the euclidean gradient VcL(C) by the inverse of the Fisher Information Matrix (FIM) I,, i.e. VcL(() C(()). For readability, we shall write L in place of L()).
The natural gradient
major issue following this approach is that Ç is unconstrained. Think of Gaussian variational posterior: under the above update, there is no guarantee that the covariance matrix is iteratively updated onto a symmetric and positive definite matrix. As discussed in the introduction, manifold optimization is an attractive possibility.
3 RELATED WORK
In Tran et al. (02 a d-dimensional Gaussian distribution W(, E), provides the fixed-form of the variational posterior 9ç (n, vec(E)). There are no restrictions μ yet the covari- ance matrix Σ constrained to the manifold M symmetric positive-defing matric es, M = { Rdxd Σ DT 0}. see e.g. (Abraham et al. 2012 Hu et al. 2020).
2

---

The exact form of the Fisher information matrix for the multivariate normal distribution is, e.g., provided in (Mardia & Marshall, 1984) and reads
Σ
1 tr
02 Jos
Ie
Σ
Σー
(I)
I<(E)
2
dow
with Ie (<(2),,, (.))...r. being the generic element of the de x dl matrix Ie (E). The MGVB method relies on the approximation (x) s (Σ-1 & ) which leads to a convenient approximate form of the natural gradients of the lower bound with respect to / and Σ,, respectively computed as A Vz((()) = vec vec ((I @ V (()) = [()E (2)
V)C(()= = VV (()) and
where @ denotes the Kronecker product. In virtue of the natural gradient definition, the first natural gradient is exact while the second approximate. Thus, Tran et al. (2021a) adopt the following
updates for the parameters of the variational posterior:
(3) where (x(() denotes suitable retraction for Σ on the manifold M. Momentum gradients can be used in place of plain natural ones. particular, the momentum gradient for the update of I relies on a vector transport granting that at each iteration the weighted gradient remains in the tangent space of the manifold M. Refer to Section 4.2 for more information on retraction and vector transport. Besides the relationship between EMGVB and MGVB already discussed, method of handling the positivity constraint diagonal covariance matrices is VOGN optimizer (Khan et al., 2018b Osawa et al. 2019). VOGN relates to the VON update (see Appendix B it indirectly updates and Σ from the Gaussian natural parameters updates. Following a non-Black Black-Boo approach, VOGN uses some theoretical results on the Gaussian distribution to recover an update for Σ that involves the Hessian of the likelihood. Such Hessian is estimated as the samples' mean squared gradient, granting the non-negativity of the diagonal covariance update. [Osawa et al. (209)) devise the computation of the approximate Hessian in a block-diagonal fashion within the layers of Deep-learning model.
= + BV (((( and Σ= ==
Lin et al. (2020) extend the above to handle the positive definiteness constraint by adding an ad- ditional term to the update rule for Σ, applicable to certain partitioned structures of the FIM. The retraction map in (Lin et al., (2020) is more general than (Tran et al. 2021a) and obtained through different Riemann metric, from which MGVB is retrieved as a special case. As opposed to EMGVB, the use of the RT in (Lin et al. 2020) requires model-specific computation or uuto-differentiation. See (Lin et al., 2021 for an extension on stochastic, non-convex problems. Lin et al. (2020) under- line that in Tran et al. (2021a) the chosen form of the retraction not well-justified as it specific for the SPD matrix manifold, whereas the natural gradient is computed for the Gaussian manifold. An extensive discussion on this point and its relationship with the EMGVB optimizer here proposed is found in Appendix D.3
Alternative methods that rely on unconstrained transformations (c.g. Cholesky factor) (e.g Tan 2021), the adaptive adjustment of the learning rate (e.g. Khan & Lin 2017) lie outside the manifold context here discussed. Among the methods that do not control for the positive definiteness constraint, the QBVI update (Magris et al., (022)) provides comparable black-bock method that, despite other black-bock VI algorithms, uses exact natural gradients updates obtained without the computation of the FIM.
4 EXACT MANIFOLD GAUSSIAN VB
Consider a variational Gaussian distribution yx with mean μ and -ositive-definite covariance matrix" Σ. Be X Σ Σ- and X2 - its natural parameters and define À = (Aı, vec()2)). The
'We present the MGVB optimizer by exactly following Tran al. (Tran et 2021a) there a typo as their I. ->- term reads
2021a
Lin et al. (2020) assert that in place of 2(2 @ Σ I))
3 Σ
which would lead to the actual natural gradient 2 E (e.g. Barfoot 2020). While their observation is valid, we argue that the omission of the constant is embedded in the approximation, as it is also omitted from the implementation codes for MGVB, where Vx is computed as 3V To clarify, V..L- = 2 is an exact relationship. while VzL >V not.
"This is the general case of practical relevance in applications, ruling out singular Gaussian distributions. For such peculiar distributions, 2 is singular, - does not exist, and neither does the density. Though this might he theoretically interesting to develop, the discussion is here out of scope. Assuming E to be positive- definite is restrictive and aligned with (Tran al. 202
3

---

corresponding mean expectation parameters (m1, vec(m2)) are given by m. gg [ and m2 En [00T] np - Σ. When required, in place of the somewhat vague notation L whose precise meaning is to be inferred from the context, we shall use ((m) to explicitly denote the lower
bound expressed in terms of the expectation parameter m, opposed (()) expressed in terms of X.
Proposition 1 For differentiable function L, and being a Gaussian distribution with mean and covariance matrix S,
VuC = LV
vz 22, -2V
where X2 = denotes the second natural parameter of qx-
The covariance matrix Σ is positive definite, its inverse exists and it as well symmetric and positive definite. Therefore E-1 lies within the manifold M and can be updated with suitable retraction
algorithm as for E in equation 3
Σー RE- (3VE-1 = - (15-2-))
(4) Opposed to the update in eq. which relies on the approximation (Z) s E-' for tacking positive-definit update of Σ, we target at updating -ー1 for which its natural gradient is available in am exact form, by primarily exploiting the duality between the gradients in the natural and expectation parameter space (Appendix D.I eq. 25) that circumvents the computation of the FIM.
For coherency with the literature VI for Bayesian deep learning (e.g. |Ranganath al. 2014 among many others), we specify the variational posterior in terms of the covariance matrix >, but update Σ-1 Yet nothing prevents specifying ςx in terms of its precision matrix Σ-1, as is often the case in Bayesian statistics textbooks, in which this case, the update 4 corresponds to an update for the actual variational precision parameter.
For updating y is reasonable adopt plain SGD-like step driven by the natural parameter V VV, as in (Tran et al. 2021a). We refer to the following update rules as Exact Manifold Gaussian Variational Bayes, or shortly EMGVB,
P++1 14, S.EVV n.C, and E = (-28
(5)
where the gradients w.r.t. L. intended as evaluated at the current value of the parameters, e.g- VrLC, = V= With respect to the MGVB update of Tran et al. (2021a), there are no approximations, e.g. regarding the FIM, yet the cost of updating Σ appears to be that of introducing additional inversion for retrieving E that is involved in the EMGVB update for μ. In the following Section, we show that with a certain gradient estimator such an inversion irrelevant. Furthermore, in Appendix B we point out that covariance matrix inversion is implicit in both MGVB and EMGVB due to the second-order form of the retraction and also show that the update for μ is optimal in the sense therein specified.
4.1 IMPLEMENTATION
We elaborate on how to evaluate the gradients v. and V,C. We follow the Black-box approach (Ranganath al. 2014) under which such gradients are approximated via Monte Carlo (MC) sam- pling and rely on function queries only. The implementation of the EMGVB updates does not require the model's gradient to be specified nor to be computed numerically, e.g. with backpropagation. By
use of the so-called log-derivative trick (see e.g. (Ranganath et al., (014)) it is possible to evaluate the gradients of the LB as an expectation with respect to the variational distribution. In particular, for generic differentiation variable Ç.
p p00)p(y10)
Vel(C) = EEy, cellog qe (0) (O),, where h0) = log
4c (0)
In the EMGVB context with q ~ ( Σ). < = (1,, vec(E)) and L(() = C(u, E). The gradient
of the L w.r.t. ( evaluated at < G. can be easily estimated using S samples from the variational posterior through the unbiased estimator
5 [eelog g (0.)) bc (o.) 0, ~ (N),
VcL(G)= = Vef(O)<=G. ~
(6)
S
4

---

where the h--function is evaluated in the current values of the parameters, i.. in (W/ vec (E,) ))
For a Gaussian distribution q ~ Mp, Σ) it can be shown that (e.g. Wierstra et al. 2014 Magris
et al. 2022):
V a log (0) Σ- u),
V log (() = 2 Σー ---000--11-----1- 1)) - -1) --)=
(8)
2
Equations 7 along with and Proposition
immediately lead to the feasible natural gradients
estimators:
1
[00 (.)h..(0.)..
(9)
a=1
V5.- - -V2
1 Σ
E (
- E, -
Ny (Ba)
(10)
S
As for the MGVB update, the EMGVB update applies exclusively to Gaussian variational posteriors, yet no constraints are imposed on the parametric form of p. When considering a Gaussian prior, the implementation of the EMGVB update take advantage of some analytical results leading MC estimators of reduced variance, namely implemented over the log-likelihood log p(yl0.) rather than the h-ffucciinn
In Appendix [.2] we show that, under a Gaussian prior specification, the above updates can be also implemented in terms of the model likelihood than in terms of the h-function. The general form of the EMGVB updates then writes:
PpG(G) = cp. + E5(0
1
[00 14) log f(0))]
(11)
s=1
V .f(G) ~e Cr. +
1
Σ
E (x
(0, HD) [0) log f(0.)
(12)
S
where, (i) if p is Gaussian Cx. = - + Eo Cu. = -5,50 (14 n0) and log f(0.)
=
log p(y)0,), (ii) if p is Gaussian not Cx. Cu. 0 and log f(0.) he. (0.).
and 0g ~ 9c. = (W,,, Σx), 8 = 1,..., S. log(y/0,) and h;. (0.) respectively denote the model likelihood and the h-function evaluated in 0,. Note that the latter depends on it involves the variational posterior, evaluated at the parameters at the value of the parameters for iteration t. denotes the prior. It is clear that under the Gaussian specification for p the MC estimator is of reduced variance, compared to the general one based on the h-function. Note that the log likeliooddood case does not involve additional inversion for retrieving in Cu. as anyway required in the second-order retraction (for both MGVB and EMGVB). This aspect is further discussed Appendix For Inverting 1 we suggest inverting the Cholesky factor --1 of E-1 and compute Σ as L.I." This takes advantage of the triangular form of - which can be inverted with back-substitution, which is A //3 flops cheaper than inverting -1 but still O (k**). L is furthermore used for generating the draws 0g as either 0g = + Ls or 0. = + I.- T with e, with ~ (0, I). As outlined in Appendix A we devise the use of control variates reduce the variance of the stochastic gradient estimators.
Though the lower bound is not directly involved in EMGVB updates, it can be naively estimated each iteration as
c. Σ
((()) log p(y log 9: (0)).
(13)
As discussed in Appendix A L, is needed for terminating the optimization routine, verifying anoma-
lies in the algorithm works (the LB should actually increase and converge) and comparing EMGVB with MGVB see Section 5
5

---

4.2 RETRACTION AND VECTOR TRANSPORT
Aligned with Tran et al. (02 we adopt the retraction method advanced in (Jeuris et al., 2012) for the manifold M of symmetric and positive definite matrices
(g) - + + 2 5 where E Tg-MMM
(14)
with & being the rescaled natural gradient V IC v.C In practice, whenever applicable, as e.g. in the retraction, for granting the symmetric from of a matrix (or gradient matrix) S, we compute as 1/2(5 + ST). Vector transport is as well easily implemented by
T3 -
EGE", where E (2))
2,)² & T5 -M.
(15)
We refer to the Manopt toolbox (Boumal et al., 2014) for the practical details of implementing the above two algorithms in a numerically stable fashion. This translates into the momentum gradients
Vic L++1 w WT, -
(ynen f.) (1 V (
(16)
Fimom L++1 w mmam L. + (1 -1--w w
(17)
where the weight 0 < w < is hyper-parameter.
The attentive reader may recognize the adoption of the form of retraction and parallel transport obtained from the SPD (matrix) manifold on the natural gradient obtained from the Gaussian manifold. This apparent inconsistency in mixing elements of different manifold structures dis- cussed in Appendix D.3 We show that, from a learning perspective, the discrepancy between the form of the SPD manifold Riemann gradient Σ V ent Vx IC -2 -2 -2 ---1 is absorbed in the learning rate In particular, our update rule can be derived within a fully consistent SPD manifold setting by updating (11, 22-
-V and the natural gradi-
In the above view, we can now further clarify that the wording "Exact" in EMGVB is twofold. (i) In Tran et al. (2021a) the natural gradient Σ is in place of the actual one 2ΣΣ- E, whose corresponding one for Σ is the one that EMGVB actually adopts. (ii) Even by the adoption the actual natural gradient -2 14 as of TTan et al. (2021a), are not well-justified: in Appendix [D.3] these are justified, and EMGVB is shown to be consistent approach. Note that EMGVB is exact in the sense of the above, yet still approximate in absolute terms due to the use of retraction. Retractions are approximate forms of the exponential map tracing back vectors the tangent space to the manifold, which is generally cumbersome transform to compute and impractical (e.g. Absil et al. 2009 Hu et al. 2020).
the use of the SPD retraction and vector transport forms
Algorithm| I summarizes the EMGVB update for the Gaussian prior- posterior case. Com- putational aspects are discussed in Appendix B.2
Algorithm EMGVB implementation
I: Set hyper- parameters: B, w S.
Set type of gradient estimator, i.e. function log (0=)
3: Set prior p(0), likelihood form p(y]0), and initial values μ4, E-1
4: t == 1, Stop = false
5: Generate: 0. ,,,, 1...S
N,,
6: Compute: 0N- = VV Or = - 22
b eqs.
11 12
7: my Iu. my-1 ~2-1
initialize momentum
8: while Stop true do
9: 11 = 1 Brr4,
EMGVB update for μ > EMGVB update for Σ
>
10: D --==--1- --!"" E-1 = Ry-1 (Bmm.-1)
11:
Generate: O, S.. 1...S S
12:
Compute: gN. -1:-1
as in line
13:
mn wm, + (1 w)g
cq.
my T + (mg 1) (1 w)
cq.
17
15:
Compute: L.
eq. > see Appendix
13
16:
t=t 1. Stop = Seut (C, P, tmax)
17: end while

---

4.3 FURTHER CONSIDERATIONS
Along with the choice of the gradient estimator and the use of momentum, there are other aspects of relevance in the implementation of EMGVB. Details are discussed in Appendix
4.4 ISOTROPIC PRIOR
For mid-sized large-scale problems, the prior is commonly specified as an isotropic Gaussian of mean 140, often 140 O, and covariance matrix E5- = rI, with T > 0 scalar precision parameter. The covariance matrix of the variational posterior can be either diagonal or not. Whether full co- variance specification (df d parameters)
provide additional degrees of freedom that can gauge models' predictive ability, a diagonal posterior (d parameters) can be practically and computation- ally convenient to adopt c.g. in large-sized problems. The diagonal-postene
assumption is largely
adopted in Bayesian inference and VI (e.g. Blundell et al., 2015; Ganguly & Earp 2021
Tran et al,
2021b) and Bayesian ML applications (e.g. Kingma & Welling 2013] Graves 2011
Khan et al.
2018b
Osawa et al. 2019), in Appendix we provide block-diagonal variant.
4.4.1 ISOTROPIC PRIOR AND DIAGONAL GAUSSIAN POSTERIOR
Assume a d-variate diagonal Gaussian variational specification, that is g ~ N(n, Σ) with diag(E) = 0², Eij = 0, for i,j = ,,...,d d and + j. In this case, Σ-1 diag (1//²)) where the division is intended element-wise, and VyL = diag( V =),,, is now d x 1 vector. Updating Σ-1 amounts updating o 2; the natural gradient retraction- bassdd update for 2 is now based on the equality
V 2C -V so that the general-case EMGVB update reads
0 R 1 (-29V
and M++1 Mt + o V
(18)
where denotes the element-wise product. The corresponding MC estimators for the gradients are
8
24 oZ ~ co² + 2
1a (0, - 0 lo
(19)
5
8
v F 2 Cp +
E1(0.
10. Mx) Po
(20)
a=1
where co? = = T, Cu. = To² @ (p), 0, ~ NNu, diag(o")). S, (0, m1)² is intended element-wise, and 1, = (1, (I,..., 1) E Rdd In the Gaussian case with general diagonal covariance matrix, retrieving 6² from the updated o a- is inexpensive as of = 1/o, -2 indicating that in this context the use of the h-function estimator never advisable.
4.4.2 ISOTROPIC PRIOR AND FULL GAUSSIAN POSTERIOR
Because of the full form of the covariance matrix, this case is rather analogous the general one. In particular, factors Cu, and ex, in eq. 29 are replaced by (i) CE, - T, Cu: TEi, mo) or (ii) CS, = 0, C₁. = 0, respectively under the Gaussian- case (log (0.) = log p(y]0x) and the general one (log fi (e.) h(0).. The MC estimators 7 and g apply: (i) leads to estimator of reduced variance, while (ii) is identical to the general case.
5 EXPERIMENTS
We validate and explore the empirical validity and feasibility of our suggested optimizer over four datasets and 12 models. These include logistic regression (Labor dataset), different volatility mod- els S&P 500 returns (Volatility dataset), and linear regression on Stock indexes (Istanbul data). Details on the datasets and models are summarized in Appendix [1 The main baseline for model comparison is the MGVB optimizer and (sequential) MCMC estimates representative of the true posterior. Additionally, we also include results related to the QB optimizer| Magris et (2022). In this section, report synthetic results on two tasks: logistic regression (classification) and volatility
7

---

modeling with the FIGARCH model (regression). Results the other datasets and models appear in Appendix C.4 Matlab codes are available at github. com/ blinedd
The bottom rows in Figures clearly show that results align with the sampling- based MCMC results and with the Maximum Likelihood (ML) estimates. Whereas marginal posterior approxi- mations rather close between EMGVB and MGVB, the top row in Figures T 5 show that the parameter learning is qualitatively different. The panels in Figure [ the LB optimization pro- cess across the iterations. In diagonal posterior setting, MGVB exact and aligns with EMGVB (middle panel), however for non-diagonal posteriors, EMGVB's lower bound shows an improved convergence rate on both the training and test sets (left and right panels respectively). Furthermore, we observe that the adoption of the h-function estimator has a minimal impact. From the point of view of standard performance measures, Figure 2 shows that compared to MGVB, at early itera- tions, EMGVB displays a steeper growth in model accuracy, precision, recall, and fl-score both on the training test and test set. Ultimately EMGVB and MGVB measures converge the same value, yet the exact nature of the EMGVB update leads to convergence after approximately 200 iterations on the training set as opposed to 500 for MGVB. A similar behavior is observed for the FIGARCH(1,11) model, the top of Figure 5
-360
Lowert bound
Lowee bound
-300
EMGNN FMGXI
-300
Lower
SCGVB
EAGGVE Tol
MGV
-145
MOVI Tick
-400
-400
200
400
600
800
1000
600
600
800
1000
400
600
800
1000
Dieration
Iteration
Iteeattion
30
A2
8g
MGVB
EGVB
2
MCMC
ML
0.5
1
1.5
-2
-1.5
-1
00.4
0.2
0
-0.5
Figure 1: Top row: lower bound optimization. Bottom row: variational posteriors (for four of the eight parameters).
Train performance
72%
Testperiormance mamnace
70%
EMGVE 30VVB
68%
64%
Accurney
66%
61%
Procision
Rerall
58%
62%
55%
200
400
600
800
1000
1200
200
400
600
800
1000 1200
Iteration
Iteration
Figure 2: EMGVB and MGVB performance on the Labour dataset across the iterations.
Tables I and 2 report such performance measures for the optimized variational posterior along with the value of the maximized lower bound. EMGVB is very close to the basclines and well-aligned with the MCMC and MI. estimates, which, along with the estimates (in Table S for the logistic regression), show that EMGVB converges towards the same LB maximum, with comparable pre- dictive ability with respect to the alternatives. It thus not surprising that the estimates, performance metrics, and value of the optimized LB are similar across the optimizers: they all converge to the same minimum but in a qualitatively different way. Also estimated variational covariance matrix for the full-covariance cases, closely replicates the from the MCMC chain (see tables 4 G). For the diagonal cases, MCMC and ML covariance matrices are not suitable for a direct comparison (see Appendix C.3 The sanity check in Figure furthermore shows that the learning of either the
8

---

mean and covariance variation parameters is smooth and steady without wigglings anomalies. As expected, the non-diagonal version leads to faster convergence while the use the h-function esti- mator slightly stabilizes the learning process. In Table 7 we also show that the impact the number of MC samples has on posterior means, likelihood, performance measures, and the optimized lower bound is minor for both the training and test phases.
Train
Test
[(0))
Accuracy
Precision
Recall
[((*)
Accuracy
Precision
Recall
fl
EMGVB
-356.642
0.713
0.712
0.703
0.708
-134,814
0.698
0.679
0.674
0.676
MGVB
-356.642
0.712
0.703
0.708
-134.814
0.698
0.679
0.674
0.676
QBVI
-35.6442
0.713
0.712
0.703
0.708
-134.804
0.698
0.679
0.674
0.676
MCMC
0.711
0.710
0.701
0.706
0.698
0.679
0.674
0.676
ML
0.709
0.708
0.699
0.704
0.709
0.708
0.699
0.704
Table I: Optimizers' performance for the Labor data on the train and test sets. See Appendix C.4 for extended results, including the use of the h-function estimator, diagonal and block block-diagonal covariance specifications.
-605
Tower bound
Lowee bomd
MSS (teii))
(toat)
-606
25
3.30
2009
ENCVB Train
EMGVB Test
ENGVE
MGVB Train
MGVE Twst
NCVU
-2000
-607
0
200
400
600
800
1000
1200
200
400
600
800
1000
1200
200
400
600
800
1000
1200
Dierationn
Iteraiion
Iraertionn
d
20
15.
10
2
MMGVE WCMC
WEV3
0
Figure 3: FIGARCH(1,1,1) model. Top row: lower bound optimization. Bottom variational marginals.
Train
C(0*)
MSE
L(0")
MSE
EMGVB
0.100
0.059
0.663
0.481
-2007.66
25.773
-604.74
3.988
MGVB
0.100
0.059
0.663
0.481
-2007.62
25.773
-604.72
3.988
QBVI
0.100
0.059
0.663
0.480
-2007.62
25.771
-604.73
3.988
MCMC
0.100
0.062
0.656
0.480
25.784
3.979
MI.
0.099
0.060
0.669
0.483
25.767
3.996
Table 2: Optimizers' estimates and performance for the FIGARCH(I,I,1) model on the Volatility dataset.
CONCLUSION
Within a Gaussian variational framework, we propose an algorithm based on manifold optimization guarantee the positive-defintt constraint on the covariance matrix, employing exact analytical solutions for the natural gradients. Our black-box optimizer results in ready-to-use solution for VI, scalable to structured covariance matrices, that can take advantage of control variates, momentum, and alternative forms of the stochastic gradient estimator. We show the feasibility of our solution on multitude of models. Our approach aligns with sampling methods and provides advantages over paaee-ofthee-ae baselines. Future research may investigate the applicability of our approach to a broader set of variational distributions, explore the advantages and limitations of the blxx framework, attempt at addressing the online inversion bottleneck of manifold-based VI.

---

REFERENCES
Ralph Abraham, Jerrold E Marsden, and Tudor Ratiu. Manifolds, tensor analysis, and applications, volume 75. Springer Science & Business Media, 2012.
P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. In Optimization Algorithms Matrix Manifolds. Princeton University Press, 2009.
Oguz Akbilgic, Hamparsum Bozdogan, and M Erdal Balaban. A novel hybrid rbf neural networks model as forecaster. Statistics and Computing, 24(3):365-375, 2014.
Timothy D Barfoot. Multivariate gaussian variational inference by natural gradient descent. arXiv preprint arXiv:2001.10025, 2020.
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisti- cians. Journal of the American statistical Association, 112(518):859 877, 2017.
Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight uncertainty in neural network. In International Conference Machine Learning, pp. 1613-1622, 2015.
N. Boumal, B. Mishra, P.-A. Absil, and Sepulchre. Manopt, a Matla toolbox for optimization on manifolds. Journal of Machine Learning Research, 15(42):1455-1455) 2014. URL https: //wwww w/www. manopt......... org
Fulvio Corsi. simple approximate long-memory model of realized volatility. Journal of Financial Econometrics, (22):1441999 2009.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research, 12(7), 2011.
Ankush Ganguly and Samuel W. F. Earp. An introduction to variational inference. arXiv preprint arXiv:2108. 13083, 2021.
Alex Graves. Practical variational inference for neural networks. Advances in neural information processing systems, 24, 2011.
Matthew D. Hoffman, David M. Blei, Chong Wang. and John Paisley. Stochastic variational infer- ence. Journal of Machine Learning Research, 2013.
Reshad Hosseini and Suvrit Sra. Matrix manifold optimization for gaussian mixtures. Advances in Neural Information Processing Systems, 28, 2015.
Jiang Hu, Xin Liu, Zai-Wen Wen, and Ya-Xiang Yuan. A brief introduction manifold optimiza- tion. Journal of the Operations Research Society of China, 6(2):199-249 2020.
Ben Jeuris, Raf Vandebril, and Bart Vandereycken. A survey and comparison of contemporary algorithms for computing the matrix geometric mean. Electronic Transactions on Numerical Analysis, 39(ARTICLE)3379-400 2012.
Mohammad Khan and Wu Lin. inference non-conjugate models to inferences conjugate models. In Artificial Intelligence and Statistics, pp. 878-887. Proceedings of Machine Learning Research. 2017.
variational inference: Converting variational
Mohammad Khan, Didrik Nielsen, Voot Tangkaratt, Wu Lin, Yarin Gal, and Akash Srivastava. Fast and scalable bayesian deep learning by weight-perturbatinn in adam. In International Conference on Machine Learning, pp. 2611-2620. Proceedings of Machine Learning Research, 2018a.
Mohammad Khan, Didrik Nielsen, Voot Tangkaratt, Wu Lin, Yarin Gal, and Akash Srivastava. Fast and scalable bayesian deep learning by pertutnnntt adam. In International Conference on Machine Learning, pp. 2611-26222 1-2620, 2018b.
Mohammad Emtiyaz Khan and Didrik Nielsen. Fast yet simple natural-gradien. descent for vari- ational inference in complex models. In 2018 International Symposium on Information Theory and Its Applications, pp. 31-35, 2018.
10

---

Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Diederik P. Kingma and Max Welling. arXiv:1312.61 2013.
Auto-encoding variational bayes.
arXiv preprint
Alp Kucukelbir, Rajesh Ranganath, Andrew Gelman, and David Blei. Automatic variational infer- ence stan. Advances in neural information processing systems, 28, 2015.
Jouko Lampinen and Aki Vehtari. Bayesian approach for neural networks-reviee and case studies. Neural networks, 14(3):257-273) 2001.
Erich L. Lehmann and George Casella. Theory of point estimation. Springer texts in statistics.
Springer-Verlag, New York, second edition, 1998. ISBN 9780387985022.
Wu Lin, Mark Schmidt, and Mohammad Emtiyaz Khan. Handling the positive-definite constraint in the bayesian learning rule. In International Conference Machine Learning, pp. 6116-6122 6126. PMLR, 2020.
Wu Lin, Frank Nielsen, Mohammad Emtiyaz Khan, and Mark Schmidt. Structured second-order methods via natural gradient descent. arXiv preprint wxXi::2107.1088/ 2021.
David J. C. Mackay. Probable networks and plausible predictions review of practical bayesian methods for supervised neural networks. Network: Computation In Neural Systems, 6:469-505, 1995.
David John Cameron Mackay. Bayesian methods for adaptive models. PhD thesis, California Institute of Technology, 1992.
Martin Magris, Mostafa Shabani, and Alexandros Iosifidis. Quasi black-box variational inference with natural gradients for bayesian learning. arXiv preprint arXiv:2205./1568 2022.
Kanti V Mardia and Roger Marshall. Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 71(1):135-146, 1984.
Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih. Monte carlo gradient esti- mation in machine learning. Journal of Machine Learning Research, 21(132):1-62, 2020.
Thomas Alvin Mroz. The sensitivity of empirical model of married women's 's hours of work to economic and statistical assumptions. PhD thesis, Stanford University, 1984.
Arkadij Semenović Nemirovskij and David Borisovich Yudin. Problem complexity and method efficiency in optimization. Wiley-Interscience, 1983.
Kazuki Osawa, Siddharth Swaroop, Mohammad Emtiyaz. Khan, Anirudh Jain, Runa Eschenhagen, Richard E. Turner, and Rio Yokota. Practical deep learning with bayesian principles. In Advances in Neural Information Processing Systems, volume 32, pp. 1-13, 2019.
John Paisley, David M. Blei, and Michael Jordan. Variational bayesian inference with stochastic search. arXiv preprint arXiv: 220664306. 2012.
Rajesh Ranganath, Sean Gerrish, and David M. Blei. Black box variational inference. Artificial intelligence and statistics, pp. 314-822, 2014.
Garvesh Raskutti and Sayan Mukherjee. The information geometry of mirror descent. IEEE Trans- actions on Information Theory, 61(3):1451-1457 2015.
Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathemati- cal statistics, pp. 400-407, 1951.
Tim Salimans and David A. Knowles. On using control variates with stochastic approxima- tion for variational bayes and its connection to stochastic linear regression. arXiv preprint arXiv: 1401. 1022, 2014.
11

---

Tim Salimans, Diederik Kingma, and Max Welling. Markov chain monte carlo and variational inference: Bridging the gap. In International conference on machine learning, pp. 1218-1226. PMLR, 2015.
Lawrence K. Saul, Tommi Jaakkola, and Michael Jordan. Mean field theory for sigmoid belief networks. Journal of artificial intelligence research, 4:61-76, 1996.
Linda SL Tan. Natural gradient updates for cholesky factor in gaussian and structured variational inference. arXiv preprint arXiv:2109.00375, 2021.
Timo Teräsvirta. An introduction to univariate garch models. In Handbook of financial time series, pp. 17-42. Springer, 2009.
Minh-Ngoe Tran, Dang H. Nguyen, and Duy Nguyen. Variational Bayes manifolds. Statistics and Computing, 31(6):71, September 2021a. ISSN 1573-1375. doi: 10.1007/ss 1222-011-10047-1
Minh-Ngoc Tran, Trong-Nghtn Nguyen, and Viet-Hung Dao. A practical tutorial on variational bayes. arXiv preprint arXiv:2103.01327, 2021b.
Felix Trusheim, Alexandru Condurache, and Alfred Mertins, Boosting black-box variational infer- ence by incorporating the natural gradient. In 2018 24th International Conference Pattern Recognition (ICPR), pp. 19-24. IEEE, 2018.
Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends@ in Machine Learning, 1(1-2):1-305, 2008.
Daan Wierstra, Tom Schaul, Tobias Glasmachers, Yi Sun, Jan Peters, and Jürgen Schmidhuber.
Natural evolution strategies. The Journal of Machine Learning Research, 15(1):949-980) 2014.
Ming Xu, Matias Quiroz, Robert Kohn, and Scott A Sisson. Variance reduction properties of the reparameterizatiio trick. The 22nd International Conference on Artificial Intelligence and Statistics, pp. 2711-2720 PMLR, 2019.
Guodong Zhang, Shengyang Sun, David Duvenaud, and Roger Grosse. Noisy natural gradient as variational inference. In International Conference on Machine Learning, pp. 5852-5861. PMLR, 2018.
12

---

A FURTHER CONSIDERATIONS ON EMGVB IMPLEMENTATION
A.1 VARIANCE REDUCTION
As EMGVB does not involve model gradients the use of the reparametrization trick (RT) (Blundell al. 2013) is not immediate. While eq. would generally hold, the form of the EMGVB gradient estimators under the RT would differ from eqs. 7 8 we develop EMGVB as general and ready- to-use solution for VI that does not require -pdel-specific derivations, yet one may certainly enable the RT within EMGVB. Though the use of the RT is quite popular in VI and ML as it empirically yields more accurate estimates of the gradient of the variational objective than alternative approaches (Mohamed et al., 2020), note that the variance of the RT' estimator can be higher than that of the score-function estimator and the path-wise RT estimator not necessarily preferable (Xu et al.] 2019 Mohamed al., 2020). Not less importantly, note that the use of the score estimator is broader as does not require log p(y(0) to be differentiable.
Control Variates (CV) stand simple and effective approach for reducing the variance of the MC gradient estimator, e.g. (Paisley et al. 2012). The CV estimator
V Vlloog ((0.) (lo)((())))))))) c).
S
is unbiased for the expected gradient, but of equal smaller variance that the naive MC one. For i = 1,2 the optimal 6j minimizing the variances of the CV estimator is
c" = C-C(vV (0)) log p(v(0), V log ((()))/Va((V{oo)) (())).
(21)
By enabling CVs, S can be tuned to balance the estimates' variance and computational performance.
In Table 7 we asses that for logistic regression values of as little as 10 appear satisfactory, yet the iterative computation of the log-likelihood is not prohibitive we suggest the adoption of a more generous value, e.g. S xe 100. Magris et al. (2022) furthermore shows that the denominator in 21 analytically tractable for Gaussian q. reducing the variance of estimated c* and thus improving the overall CVs' efficiency. model gradients are available one may use CV along with the RT to further enhance the efficiency of the expected gradient estimation.
A.2 LB SMOOTHING AND STOPPING CRITERION
The stochastic nature of the gradient estimator introduces some noise in the estimated LB å that can violate its expected non-decreasing behavior across the iterations. By setting window of size we rather consider the moving average on the LB, L. 1/w Σ(=1 Lr-it 1, whose variance is reduced and behavior stabilized. By keeping track of max C we terminate the the learning after max & did not improved for P iterations (patience parameter) or after a maximum number of iteration (fmax) reached (stopping criterion fexit function in Algorithm
A.3 CONSTRAINTS ON MODEL PARAMETERS
EMGVB assumes Gaussian variational posterior, that is the parameters are unbounded and de- fined over the entire real line. Assuming that a model parameter is required to lie on a support S, impose such constraint suffices to identify a feasible transform T R - S and apply the EMGVB update to the unconstrained parameter // T0). Certainly, by applying VI on v we require that the variational posterior assumption holds for & rather than 0. The actual distribution for 8 under a Gaussian variational &b can be computed (or approximated with a sampling method) as (T (0); μ., Σ)) (Jr (f)),, with Jr 1 the Jacobin of the inverse transform (Kucukelbir et al.) 2015). Example. For the GARCH(1,1) model (see Section 5) the intercept the autore- gressive coefficient of the lag-one squared return a and moving-average coefficient 3 of the lag-one conditional variances need to satisfy the stationarity conditions + 3 < I and w > 0, > 0, 2 0. Such conditions are unfeasible under Gaussian variational approximation: we estimate the uncon- strained parameters thun Va, Vig where w T(t)), T(Pa) T( T(ws)). 3 T(*a)T(Wpa) with T(x) = exp(x)/ exp(x)) for x real, on which Gaussian's apply.
assumptions
13

---

A.4 GRADIENT CLIPPING
Especially for low values of S, and even more if variance control method is not adopted, the stochastic gradient estimate may be poor and the offset from its actual value large. This may result in updates whose magnitude is too big either iterations, and with poor initial values, this issue may c.g. cause complex roots in eq. 15 At each iteration t, control for the magnitude of the stochastic gradient 9. we rescale its Io-norm ll Il whenever it is larger than fixed threshold (max by replacing Ôt with m lldell which preserves its norm. Gradient clipping
positive or negative direction. Especially at early
be either applied to the gradients Vp, Vy or to the natural gradients ΣV μ. -2 and in any case before obtaining momentum gradients. We suggest applying gradient clipping readily to V, Vx to promptly mitigate the impact that far-from-the-mean estimates may have on successive computations.
A.5 ADAPTIVE LEARNING RATE
is convenient to adopt an adaptive learning rate or scheduler for decreasing 3 after a certain number of iterations. Typical options are that of reducing / by certain factor (e.g. 0.2) every set number of iterations (e.g. 100), or decrease it after iteration l' e.g. by setting 3, min (3, 34 where l' is a fraction (e.g. 0.7) of the maximum number of iterations max allowed before the LB optimization stopped.
A.6 CLASSIFICATION VS. REGRESSION
We point out that the EMGVB framework is applicable to both regression and classification prob- lems. In generic DL classification problems, predictions are based on the class of maximum proba- bility which is computed by applying a softmax function at the last layer returning to the probability (:())) of certain class Cj for the i-th sample, i = 1,..., M. From these probabilities it is straight- forward to compute the model log- ligllikld as E Vi.come log Ps (Cmu). with yi,cow representing the one-key-hot encoding of the i-th sample, whose true class Caue- For regression, the paramet- ric form of log p(y is clearly different and model-specific (e.g. regression with normal errors as opposed to Poisson regression, with the latter being feasible as the use of the score estimator does not require the likelihood to be differentiable). Note that however additional parameters may enter into play besides the ones involved in the back-bone forward model: e.g. for regression with normal errors tackled with artificial neural network, the Gaussian-form likelihood involves the regression variance, which is an additional parameter over the network's 's ones, or for Student-t errors the degree of freedom parameter y (with the constraint 2). See the application in Appendix C.3
A.7 MEAN-FIELD VARIANT
Assume that for d-variate model the Gaussian variational posterior is factorized as
(c( = 9c, (0)) 9c, (0h) II 9c, (0%). 1=1
with h < k. I = d this corresponds to dull-diagonal case where each 0 is scalar and the the variational covariance matrix E of 93
covariance between ......,St is ignored. If matrix, and some of the Os are indeed vectors. In any case, the
corresponds to block-diagone expected gradients with respect to each block of parameters can be computed independently, given the scalars hGg (0) or log p(y(0), depending on whether the f--function estimator is used. For Gaussian prior, its covariance matrix can be diagonal, full, block-diagonal with a structure matching not that of S. Eqs. I1 T2 with the condition 29 can be used starting point to derive case- specific EMGVB variants based on the form of the prior covariance.
Algorithm 2 summarizes the case with isotropic Gaussian prior of zero-mean and variance т, using the gradient estimator based on the log- ligeliiihhood: μ4, Σ, (E,-) respectively denote the mean and covariance (precision) matrix of the i-th block of Σ. In this case, the block-wise natural gradients
14

---

estimated
8 [(0,
= 14
24+ 5
mb) 1(g)()))
s=1
v + diag(7 ( +
2 Ex (0,, 14) us)) -x))) log
log p(y p(y)
S
#=1
where 8, is a sample from the variational posterior. 6, can be obtained by concatenating marginal samples from each block, 82 ........,. 0,,l, with 8s, 9u. i I.....h
Algorithm 2 EMGVB for block-diagonal covariance matrix (prior with zero-mean and covariance
matrix r
1: Set hyper- paper-parameters: < 3, < S
2: Set the type of gradient estimator, i.c. function log (0x)
3: Set prior p00 0, r). likelihood p(y 0). and initial values μa, Σ
4: 1 Stop false
Generate: 0, = 10.y..... Osal. Os, ~ 9,, 1...S, 1.....
6: Compute: log p(y]0.)
7: for i= I,..., h do
Compute: 9r., Σ, Vnih C. 85;1 x
9: mp; = 9p:- my,- = gx,
10: end for
11: while Stop tru do
12:
L=
13:
for i = 1,..., do
14:
H4 M4 Bmp
Vold, 2 2 R2z1 (mm=r1)
16:
end for
17:
Generate: 0, 10,..... Bop On ~ Du, 1,.,,,;;= l....,
18:
Compute: logp(y)0.). log2(0.)
for l,...,
20:
Compute: Juc 9y,
21:
My, wm + (1 - w)
22:
-z TT-- -2,1 (z-1) + w)gs-1
23:
 = L+ log logp(y)(s)) log qm, F, (0..)
24:
1= 1 + 1. Stop far (C, bmxx)
25:
end for end while
26:
B OPTIMALITY AND EFFICIENCY
B.1 OPTIMALITY
Several authors (e.g. Khan & Lin 2017 Khan et al. 2018a) obtained update rules for VI by devel- oping over SGD-like updates for the natural parameters of the variational posterior. By updating the natural parameter À and exploiting its definition, it is relatively simple recover the update rules for and Σ. Indeed from the SGD update for the natural parameter Aet1 = A + (()) it follows
Mt+1 = Et+1 Σ, + 3 [AC(A))]
Ex+1 [(S]] 2,8 V.: + BIV C(,O]]
(22)
and
Σd1 = E Vz
(23)
By replacing Σ, [ with Σ, in the update for μ. Khan & Lin 2017 (2018a) obtain
Khan et al.
P++1 14 0 (V ((A:))
(24)
15

---

Eq. 23 does not apply to the EMGVB update the update for E- is carried out with retraction and does not result from an SGD update of the natural parameter X2 = - yet the form of equation 22 does apply. We refer eq. - as indirect update since derived from the natural parameter update. Note that, the update exploits E, resulting a one-step forward-looking rule. It is relevant investigate whether the update 22 is preferable the EMGVB update for μs. Intuitively one might expect that eq. [22] is preferable as it somewhat readily exploits the updated I value as soon it becomes available. The following theorem however proves that this is the case, a proof is provided in Appendix D.4
Theorem For the Gaussian distribution with parameters = (1,, vec()) the oplimization prob- lem
S++1 aang Vę(((s)) 1
/3 DKL(PcllPE.)
where (:--) denotes the inner product and V = (Vuf()), ve Mc==. convex with respect to Ç. The optimum update for μ is available in closed form and analogous that of the EMGBV update.
The objective in Theorem I is that of the mirror descent developed by Nemirovskij & Yudin (983). where non-Euclidean geometry is induced by considering penalized optimization ob- tained through proximity function such as the Bergman divergence, which equals the KL divergence for exponential- famill distributions. In this regard see (Raskutti & Mukherjee 2015).
Following Theorem I the update for μ in eq. 5 is optimal in terms of the above objective and perhaps counter-intuitively the indirect forward-looking update 22 is proved to provide non-optimal steps toward the maximization of the lower bound. The EMGVB update for is thus preferable over the alternative recovering an indirect update rule for jı starting from an SGD update on the natural parameter as e.g. in (Khan & Lin 2017 Khan al., 2018a): in the above terms of Theorem I the EMGVB update for y is the best one could take.
B.2 COMPUTATIONAL ASPECTS
In terms of computational complexity, the exact EMGVB implementation is at no additional cost. Actually, the cost of computing the natural gradient in EMGVB as -2 y L is cheaper than the one in MGVB, LV, (k*)) operations for each matrix multiplication. However, both MGVB and EMGVB share a cumbersome matrix inversion.
Going back to eqs. 12 and II it is noticeable that under the most general estimator based on the h-function, Σ not involved in any computation, neither in the gradient involved in the updated for Σ-1 nor in that for μ., suggesting that implicitly the EMGVB optimization routine does not require the inversion of E-1 The above point however ignores that the update for Σ-1 masked by the underlying retraction. For the retraction form in eq. 14 both E-1 and its inverse Σ are needed, thus implying matrix inversion at every iteration. That is, the covariance matrix inversion is implicit in MGVB and EMGVB methods, which both require - and Σ at every iteration (with little surprise, as the form of the retraction is a second-order approximation though neither eq. 7 nor involve E, the inversion of Σ-1 is
the exponential map).
With the h-function estimator still necessary, as Σ is required retraction. Similarly, the adoption the log-likelihood estimator under the Gaussian regime in eq. 29 is not computationally more expensive than the h-function case, as Σ, involved in Ca. is anyway required in retraction. As outlined in Section 4 Σ can be conveniently recovered from the Cholesky factor of --1, with a lesser number flips. Lastly, E-1 (E) is diagonal, the inversion is trivial and, when applicable, eq. 11 is preferred.
C EXPERIMENTS
C.1 RESULTS FOR THE LABOUR DATA
16

---

61
B2
Bg
81
96
MGVB
0.678
0.082
0.573
0.494
0.637
0.608
0,041
EMGVB
0.679
1.485
0,082
0,573
0.494
0,637
0.608
0,041
QBVI
0.679
1.485
0.082
0.573
0.494
0.637
0.608
0.041
MGVB
0.678
0,084
0.574
0.493
0.638
0.609
0.048
EMGVB
0.678
1.487
0,084
0.493
0.609
QB
0.678
1.487
0.084
0.574
0.493
0.638
0.609
0.048
MGVB
0.576
1.430
0.056
0.541
0.493
-0 -0.638
0.593
0.105
EMGVB
0.578
1.431
0.057
0.542
0.493
-0.637
0.593
0.104
QBV
0.579
-1.432
0.057
0.542
0.493
-0.637
0.593
0.103
MGVB" dagg
0.573
1.430
0.057
0.542
0.493
-0.635
0.593
0.109
EMGVB1 diag.
0.575
1.431
0.057
0.543
0.493
-0.635
0.593
0.107
QBVT diag.
0.576
1.431
0.057
0.543
0.493
0.635
0.593
0.107
MCMC
0.675
1.484
0.084
0.570
0.493
0.639
0.608
0.051
ML
0.679
1.476
00.055
0.571
0.485
0.625
0.399
0.045
Table 3: Parameters' estimates for the labour dataset. Top: posterior means, middle: variances, covariances (x (x10)). denotes the use the h-function gradient estimator, diag the use of diagonal variational posterior.
bottom:
MGVB
81
fg
84
96
80
1.639
-..488
0.697
0.163
0.262
0.192
2.709
81
1,634
0.298
1.425
0.451
0.020
0.084
0.118
82
-0.810
0.288
0.413
0.066
0.050
-0.064
-0.109
EMGVB
83
0.700
1.393
0.406
0.039
0.128
-0.189
00.97
BA
0.077
0.361
0.127
0.083
0.135
0.259
0.191
0.312
-0.551
0.008
0.068
-0.144
-1.322
-0.611
86
0.176
0.230
11200
0.028
By
2.809
0.048
0.103
0.381
0.218
0.633
0.073
ML
3g
81
Bg
81
A5
&
37
&
-1.578
-0.800
-0.684
0.067
0.318
0.194
22.00
-1.510
0.274
1.378
-0.422
-0.027
0.074
-0.884
B2
00.79
0.249
0.393
0.102
0.024
-0.058
-0.093
MCMC
-0.654
1.307
0.381
0.067
0.117
-0.164
00352
84
0.047
-0.414
0.111
0.072
-0.172
-0.265
-0.109
B8
0.334
0.049
0.006
0.089
-0.177
-1.286
00613
86
0.172
0.111
-0.036
-0.137
-0.266
- 11.822
-0.097
17
2.603
00.055
-0.097
-0.349
-0.119
00.90
-0.102
Table 4: Parameters' covariance matrices for the labour dataset. Entries are multiplied by 103.
0.5
350
300
250
.0.5
200
150
-1.5
50
EMGV3
BMGVB
En
-2
KMCV3L-/e fen
EMGVRR-/AK
EMCVEl-
-50
EMGVE
-2.5
-100
0
200
400
600
800
1000
1200
200
400
600
800
1000
1200
Berstion (8)
Iteration (0)
Figure 4: Parameter learning across the iterations under different variants of the EMGVB algorithm for the labour dataset.
17

---

Train
C(0*)
logplyl0*)
Accuracy
Precision
Recall
f1
EMGVB
355.642
-332.992
0.713
0.712
0.703
0.708
MGVB
335.6442
-332.991
0.713
0.712
0.703
0.708
QBVI
-356.642
332.992
0.713
0.712
0.703
0.708
EMGVR
-356.635
-332.990
0.711
0.710
0.701
0.706
MGVB
335.6335
-332.900
0.711
0.710
0.701
0.706
QBVT
336.635
332.990
0.711
0.710
0.701
0.706
EMG
358.423
333.124
0.711
0.710
0.701
0.706
MGVB
358.428
333. 139
0.711
0.710
0.701
0.706
QBV
358.419
333.
0.711
0.710
0.701
0.706
EMGVB!
358.421
333.127
0.709
0.708
0.700
0.704
MGVB dag
3358.427
-333.133
0.709
0.708
0.700
0.704
QBVI
-358.418
-333.124
0.709
0.708
0.700
0.704
MCMC
-332.990
0.711
0.710
0.701
0.706
ML
-332.983
0.709
0.708
0.699
0.704
Test
C(#*)
1og
Accuracy
Precision
Recall
fl
EMGVB
134.814
113.854
0.698
0.679
0.674
0.676
MGVB
134.814
113.854
0,698
0.679
0.674
0.676
QBVI
134.804
113.854
0.698
0.679
0.674
0.676
EMGVB
134.805
113.855
0.698
0.679
0.674
0.676
MGVB
134.805
113.855
0.698
0.679
0.674
0.676
QBVI
134.805
113.855
0.698
0.679
0.674
0.676
EMGVB-le
136.864
114.171
0.667
0.644
0.640
0.642
MGVBligg
136.877
-114.178
0.667
0.644
0.640
0.642
QBV
136.857
-114.167
0.667
0.644
0.640
0.642
EMGVB
136.876
-114.220
0.667
0.644
0.640
0.642
MGVB"
136.883
-114.228
0.667
0.644
0.640
0.642
QBVI 136.871
114.217
0.667
0.644
0.640
0.642
MCMC
-13.8838
0.698
0.679
0.674
0.676
ML
113.800
0.709
0.708
0.699
0.704
Table 5: Models' performance for the labour data on the train and test sets. h-function gradient estimator, diag the use of diagonal variational posterior.
denotes the use of the
6
Bg
Ag
8
3
By
EMGVB
4.242
5.465
0.601
1.479
1.220
1.931
1.879
4.477
MGVB
4.240
5.574
0.615
1.506
1.288
2.127
1.927
4.393
QBVI
4.255
5.484
0.603
1.484
1.224
1.938
1.886
4.493
EMGVB
4.068
5.415
0.594
1.464
1.290
2.043
1.968
4.385
MGVB
4.069
5.416
0.598
1.469
1.287
2.058
1.974
4.387
QBVI
4.068
5.415
0.594
1.464
1.290
2.043
1.968
4.385
EMGVBaug
0.909
3.356
0.255
0.895
0.988
0.955
0.994
1,286
MGVB
0.909
3.356
0.255
0.894
0.988
0.955
0.993
1,286
0.910
3.363
0.255
0.896
0.989
0.956
0.995
1,288
EMGVBidig
0.871
3.391
0.252
0.897
1.051
0.999
1.006
1.313
MGVB rdia
0.871
3.390
0.252
0.897
1.051
1.000
1.006
1.313
QBVI diag.
0.872
3.394
0,252
0.897
1.052
1.000
1.007
1.315
MCMC
3.967
5.364
0.589
1.406
1.274
2.071
1.966
4.248
MI.
4.079
5.436
0.589
1.457
1.276
2.059
1.960
4.393
Table 6: Variances of the parameters for the labour dataset. Entries multiplied 102.
denotes
the use of the h-function gradient estimator, diag the use of diagonal variational posterior. For discussion on the diagonal case see Appendix C.3
18

---

B2
A4
86
C(06)
0.687
1.498
0,086
0.576
0,500
0.605
0.037
430,216
20
6.984
0.682
-1.498
0.085
-0.571
0.499
-0.641
0.607
429.626
30
8.137
0.684
1.493
0.087
0.573
0.495
0.640
0.609
0.043
434.863
50
9.015
0.680
1.492
0.085
0.576
0.494
0.643
0.611
0.047
435.890
75
11.788
0.682
-1 -1.499
-0.85
-05575
0.492
-0.642
0.611
0.045
- 436.252
100
14.652
0.682
1.490
0.085
0.574
0.493
0.639
0.609
0.044
436.924
22.862
0.680
1.488
0.085
0.574
0.492
0.046
436.759
200
24.169
0.681
1,487
0.085
005574
0.492
-0.638
0.608
-37.036
300
31.708
0.680
1.488
0.085
0.575
0.493
0.638
0.609
0.046
436.529
Train
(((")
lly (**)
Accuracy
Precision
Recall
fl
10
356.686
332.991
0.711
0.710
0.701
0.706
20
-356.667
-332.994
0.713
0.712
0.703
0.708
30
356.653
-332.990
0.709
0.708
0.699
0.704
356.645
332.988
0.713
0.712
0.703
0.708
-356.642
-332
0.713
0.712
0.703
0.708
100
-35..440
-332.988
0.711
0.710
0.701
0.706
150
336.6440
332.990
0.711
0.710
0.701
0.706
200
-356.638
-332 990
0.713
0.712
0.703
0.708
300
-356.639
332.989
0.711
0.710
0.701
0.706
Test
(((")
logply
Accuracy
Precision
Recall
fl
10
-134.865
-113.910
0.698
0.679
0.674
0.676
20
-134.999
-113.856
0.698
0.679
0.674
0.676
134.793
-113.827
0.698
0.679
0.674
0.676
50
-134.805
-113.838
0.698
0.679
0.674
0.676
75
-134.814
-113.854
0.698
0.679
0.674
0.676
100
-13.8800
-113.876
0.698
0.679
0.674
0.676
150
-134.791
-113.845
0.698
0.679
0.674
0.676
200
-134.788
-113.840
0.698
0.679
0.674
0.676
300
133.8332
-113.842
0.679
0.674
Table 7: Estimated parameters and performance measures on the labour dataset for EMGVB (full- posterior) for different sizes of the number of MC draws for the estimation of the stochastic gradients S. t refers to the run-timg on-time per iteration (in milliseconds), L(0)) to the LB evaluated at the initial parameters. For each S, common random seed used.
19

---

C.2 VOLATILITY MODELS
Our second set of experiments involves the estimation of several GARCH-family volatility models. The models in Table [ differ for the number of estimated parameters, the form of the likelihood function (which can be quite complex as for the FIGARCH models), and constraints imposed on the parameters. Besides the GARCH-type models, we include the well-known linear HAR model for realized volatility (Corsi 2009). We performed preliminary study for retaining only relevant mod- els, e.g. we observed that for a GARCH(1,022 A3 is not significant, so we trained a GARCH( (1,01), or that the autoregressive coefficient of the squared innovations is always significant only at lag one, so we did not consider further lags for a. For a, 3, we restricted the search up to lag 2. Except for HAR's parameter $3, all the parameters of all the models are statistically significant under standard ML at 5%. Note that the aim this experiment is that of applying VI and EMGVB to the above class of models, not to discuss their empirical performance or forecasting ability. For the reader unfamiliar with the above (standard) models, discussion and notation we refer e.g. to the accessible introduction of Teräsvirta 2009)..
As for the Labour data, we report the values of the smoothed lower bound computed at the optimized parameter L((*), the model's log-likelihood in the estimates posterior parameter p(y(0*) and the MSE between the fitted values and squared daily returns, used as a volatility proxy. Details on the data and hyperparameters
provided in Table
Figures 6 provide sample illustrations of the lower bound maximization for the GJR(I,1,1) model (perhaps the most used and effective in applications beyond the standard GARCH(1,0.D) and the FI- GARCH(1,111) model, the most complex one among our selection due to the form of the likelihood, constraints, and econometric interpretation. In general, beyond figure 6 witness slight but consistent improved convergence of the lower bound towards its maximum the train set for EMGVB with respect MGVB, the convergence of the I.B similar level the test sets and MSE that eventually converge to rather similar values but that in some cases can be quite different at early it- erations (which is expected but irrelevant in applications as 6* the measures rather analogous. These observations are quantitatively supported by the results in g where all the optimizers lead to rather similar estimates and statistics.
A visual inspection of the marginal densities e.g. in Figures and reveals that in general both EMGVB and MGVB perform quite well compared to MCMC sampling and that the variational Gaussian assumption is quite feasible for all the volatility models. Note that the skew observed e.g. Figure 7 for the w parameters and the non-standard form of e.g. v/ for the FIGARCH models due to the parameter transformation: VI is applied on the unconstrained parameters (s/w, 1%) and such variational Gaussians are back transformedd on the original constrained parameter space where the distributions are generally no longer Gaussian (Figure opposed to Figure 5
4/
2.5
0.55
1.5
0.1
0.4
EMGVB
NCMC
0.5
MGVB
-2
-2
Figure 5: FIGARCH(L,I,) model. Variational and MCMC marginals for the unconstrained param- eters, complement Figure i the main text.
20

---

.604.5
27.3
-1989
bocndd
boand
-504.55
MSK. (taiin)
27.25
1989.5
-4.49
-1960
1190.5
EMGVS Trann MG Train
604.65
ENIGVE Test MGVB Test
27.15
ENGV8 MGVS
4.47
-191
604.7
27.1
4.46
Iteration
Iteration
Iteeatiion
Figure 6: GJR(1,1,1) model. lower bound and mean squared error in the train and test set. Middle and bottom rows.
V/
1.5-
0.5
0.2
-2
13
12
-1
-5
0
10
20
10
20
60
15.
40
10-
BINGVR
20
NGVB
NCMC
0
0.05
0.1
0.1
0.2
0.2
0.4
0.6
0.6
0.8
Figure 7: GJR(1,1,1) model. MCMC and variational marginals in for the unconstrained parameters (We, thay yyy, 122) (top row) and constrained parameters (w,a, 79) (bottom row).
21

---

Train
Test
@
a
81
B2
(()")
p(y
MSE
[()*)
p(a19*)
MSE
ARCH
EMGVB
0.519
0.640
2247.14
-2241.03
29.818
-651.11
-645.75
5.447
MGVB
0,319
0.640
-2247 .14
-2241.03
29.818
-651.11
-645.75
5.447
-2247.14
-2241.03
MCMC
0.519
0.628
-2241.01
29.694
-645.66
5.391
ML.
0.519
0.640
-2241.04
29.820
-65.73
5.447
GARCH(L.0.1)
EMGVB
0.043
0.230
0.737
-2012.40
-2002.56
25.691
-607.70
-99.668
3.953
MGVB
0.043
0.230
0.737
-2012.40
-2002.56
25.693
-607.69
-398.67
3.954
0.043
0.230
0.737
-2012.41
-2002.56
25,691
-607.70
-988668
3,953
MCMC
0.042
0.226
0.738
2002..52
25.683
-99.99
3.944
ML
0.043
0.231
0.738
-2022.67
25.723
-398.43
3.957
GJR(1,1,1)
EMGVB
0.044
0.108
-1988.7 71
-1976.08
27.214
-64.556
-92.91
4.492
MGVB
0.044
0.108
0.293
0.721
11888.70
-1976.09
27.243
-604.56
-992.90
4,497
QBVI
0.044
0.108
0.293
0.721
-1988.77 70
-1976.08
27.238
-604.59
-392.93
4.497
MCMC
0.042
0.108
0.289
0.722
-1976.03
27.143
-53.229
4,476
ML
0.043
0.109
0,294
0.723
-1976.24
27,386
-592,58
4,515
G..().
EMGVB
0.045
0.116
0.323
0.049
-1988.
-1976.12
27.816
-604.33
-593.68
MGVB
0.045
0.115
0.322
0.653
0.050
11888.8
-1976.10
27.767
-64.33
-93.665
4,628
QBVI
0.045
0.116
0,323
0.653
0.050
1988..39
-1976.11
27.797
-604.35
-593.66
4.635
MCMC
0.044
0.110
0.672
0.041
-1976.01
27.267
-593338
4,515
MI.
0.044
0.114
0.321
0.670
0.036
11976.12
27.778
-593.63
4.629
EGARCH(
EMGVB
-0.003
0.414
0.929
20133.15
-2117.67
26.8411
-613.73
-599.43
3.931
MGVB
-0.003
0.416
12033.16
-2017.69
26.837
-613.68
-99.44
3.933
QBVI
-0.003
0.416
0.929
20133.16
-0117.69
26.837
-63.667
599.44
3,933
MCMC
-0.003
0.406
0.932
-2017.64
26.858
599.37
3,923
ML
-0.003
0,415
0.929
-2017.68
26.819
-599.44
3,932
EGARCH(1.1.))
EMGVB
-0.015
0.350
-0.171
0.930
-1995 54
-1974.97
163
-608.06
-589.08
4.128
MGVB
-0.015
-1995.69
-1975.21
26.231
-608.0
QBVI
-0.015
0.364
-0.174
0.926
-1995.69 69
-1975.21
26.234
-68.09
589.57
4.170
MCMC
-0.015
0.340
-0.169
0.932
-1974.92
26. 165
-58..65
4,094
ML
-0.015
0.350
-0.172
-1974.98
26.208
-89.09
4.129
EGARCH
EMGVB
-0.015
0.355
-0.173
0.917
0.011
-1998 67
-1975.03
26.235
-610.85
-589.17
4.143
MGVB
-0.016
0.395
-0.185
0.859
0.061
-1999.55
-1976.00
26.572
-611.10
-500.33
4.291
QBVI
-0.016
0.395
-0.186
0.855
0.065
-1999.55
-1976.01
26.589
-611.08
-50030
4,293
MCMC
-0.015
0.335
-0.166
0.955
-0.022
-1974.90
26.
588.63
4.078
ML
-0.015
0.351
-0.171
0.925
0.004
-1974.99
26.214
-589.07
4.128
FIGARCHOO
EMGVB
0.100
0.648
0.418
-2077 .49
-2000.06
25.673
-604.27
-598.35
3.978
MGVB
-0007.4
-60.227
3.97%
QBVI
0.100
0.649
0.419
-2007.50 .50
-2000.06
25.675
-604.26
-598,36
3.978
MCMC
0.102
0.657
0.432
2000.02
25.722
-598.13
3.971
ML
0.100
0.658
0.424
-2000.13
25.670
-98.666
3.989
FIGARCH(1111)
EMGVB
0.100
0.059
0.663
0.481
-2007.6 .62
-1999.69
25.773
-604.74
-598.44
3,988
MGVB
0.100
0.059
0.663
0.481
-2007.0
-1999.69
25.773
-604.72
-598.44
3.988
QBVI
0.100
0.059
0.663
0.480
-2007 .266
11999.69
25.771
-604.73
-598.44
3.988
MCMC
0.100
0.062
0.656
0.480
1999..77
25.784
-988.8
3.979
ML
0.099
0.060
0.669
0.483
-1999.74
25.767
-598.65
3.996
Ao
HAR
EMGVB
1.078
0.488
0.420
-0.012
-5082.83
-5660.65
14.179
-1240.01
-1220.23
18.857
MGVB
1.063
0.488
0.421
-0.011
5083.02
-0660.65
24.179
-1240.30
-1220.26
18.860
QBVI
1.788
0.488
0.420
-0.012
-082.83
-5060.65
24.179
-1240.00
-1220.23
18.857
EMGVB
1.085
0.491
0.421
-0.016
5085 54
-5660.68
24.179
-1242.81
-1220.40
18.867
MGVB
1.067
0.493
0.420
-0.016
-085.61
-5060.70
24,180
-1242.97
-1220.42
18.869
QBVI"
1.084
0.490
0.421
-0.015
50K3.54
-0660.67
24.179
-1242.81
-1220.40
18.866
MCMC
1.067
0.489
0.420
-0.011
-5060.64
24.179
-1220.04
18.855
MI.
1.079
0.488
0.420
-0.012
-060063
24.179
-1220.03
18.857
Table 8: Parameter estimates (on the actual constrained parameter space) and statistics models' performance on the train and test set. denotes the use of the h-function gradient estimator.
22

---

C.3 ISTANBUL DATASET: BLOCK-DIAGONAL COVARIANCE
In this section, we apply EMGVB under different assumptions for the structure of the variational covariance matrix. We use the Istanbul stock exchange dataset of Af et al. (2014), (details provided in ind Table To demonstrate the feasibility of the block-diagonal estimation under the mean-field framework outlined in Appendix A.7 100 index (ISE):
we model the Istanbul stock exchange national
ISE, = 3o + S.SP + 92NIK; + ByBOV, + 9,DAX, + 3EFTSE, + SgEU, + 3,EM, + Ex
with e ~ o0) and the covariates respectively correspond to the S&P 500 index, Japanese Nikkei index, Brazilian Bovespa index, German DAX index, UK FTSE index, MSCI European index, and MSCI emerging market index. We estimate the coefficients So...., 87 and the trans- formed parameter vi log(ar), from which (standard error of the disturbances) is computed as exp(1>=)+ Var(t>>) with Var(w>) read from the variation posterior covariance matrix, while for MI. regression corresponds to the residuals' root mean squared error.
We consider the following structures for the variational posterior: (i) full covariance matrix (Full), (ii) diagonal covariance matrix (Diagonal), (iii) block-diagonal structure with two blocks of sizes 8 8 and 1 x 1 (Block I) and, (iv) block diagonal structure with blocks of sizes 1 x 3 3, 2 x 2, 2 x 2 and x 1 (Block 2). Case (iii), models the covariance between the actual regressors but ignores their covariance with the regression standard error. Case (iv) groups in the 3 x block the indices traded in non- -ur-Eun stock exchanges, and in the remaining 2 x 2 blocks the indices referring to European exchanges and the two MSCI indexes. Furthermore, covariances between the intercept and regression's 's standard error with all the other variables are set to zero. The purpose of this application is that of providing an example for Algorithm _ and the discussion in Appendix A.7 rather than providing an effective forecasting model supported by solid econometric rationale. Yet, structures (iii) and (iv) correspond to quite intuitive grouping of the variables involved in our regression problem, motivating the choice of the dataset.
Tables and 10 summarize the estimation results. Table 9 shows that the impact of the different structures of the covariance matrix is somewhat marginal in terms of the performance measures, with respect to each other and with respect to the MI. estimates. As for the logistic regression exam- ple, we observe in the most constrained cases (ii) and (iii) certain the estimates of certain posterior means slightly deviate from the other cases indicating that the algorithm perhaps terminates at a different local maximum. Regarding variational covariances reported in Table TO there is remark- able accordance between the covariance structures (i), (iii) and ML, while for the diagonal structure (ii) and block-diagonal structure (iv) the covariances are misaligned with the ML and full-diagonal ones, further suggesting that the algorithms convergence at different maxima of the lower bound.
From a theoretical perspective, if is the covariance matrix of the joint distribution of the eight variates (case (i)), by the properties of the Gaussian distribution, blockse.gg in case (iv) and diagonal entries should match the corresponding elements in Σ. It is however not surprising to observe that the elements in the sub-matrices e.g. in cases (ii) and (iv) deviate from those Σ. Indeed, the results refer to independent optimizations of alternative models (over the same dataset) that are not granted to converge at the same maximum (and thus distribution). Across the covariance structures (i) to (iv) the optimal variational parameters correspond to different multivariate distributions, that independently maximize the lower bound, and that are not constrained to be related to each other. This is indeed confirmed by the differences in the maximized Lower bound Lo in Table 9 and in the different levels at which the curves in Figure are observed converge. Thus the blocks in the covariance matrix under case (iv) do not match the entries in Σ. In this light, the ML estimates' variances in the third panel of Table 10 can be compared those of case (i), but are misleading for the other cases, as the covariance matrix of the asymptotic (Gaussian) distribution of the ML. estimator is implicitly full.
23

---

8o
Bi
A2
83
B4
Ag
[
a
Full
0.001
0.098
0.079
-0.272
-0.168
-0.354
1.165
0.943
0.014
Diagonal
0.001
0.054
0.107
-0.211
0.135
-0.008
0.519
0.888
0.014
Block I
0.001
0.098
0.079
-0.272
0..16
00.35
1.164
0.944
0.014
Block 2
0.001
0.066
0.127
00.205
0.244
0.225
0.188
0.844
0.014
ML
0.001
0.099
0.078
-0.273
00.144
0.363
1.179
0.946
0.014
Train
Test
[((0)
p(y(0*)
MSE
p(y 0*)
MSE
Full
1186.079
1223.699
0.194
316.849
0.041
Diagonal
1177.810
1219.055
0.198
317.042
0.041
Block
1186.082
1223.701
0.194
316.848
0.041
Block 1173.852
1213.509
0.204
316.696
0.041
ML
1223.728
0.194
316.871
0.042
Table Posterior means, ML estimates and performance measures on the train and test set. a refers to the disturbances' standard deviation. Block corresponds to case (iii) and Block 2 to case (iv).
8
82
8g
8,
&
86
8y
8
5g
0.001
0.003
-0.001
-0.002
-0.001
0.006
-0.009
0.000
31
0.001
-0.046
-2,970
-1.858
-0.789
-0.186
1.603
0.009
32
0.003
-0.032
1,094
0.691
44.137
-0.001
-2.992
1.094
0.897
0.596
-0.949
-4,993
-0.017
ML
34
-0.02
-1.906
0.297
0.961
3.687
-20.246
-0.338
-0.032
3g
0.002
-0.729
0.679
0.622
3.836
-28.377
-1.364
0.032
56
0.007
-0.170
-0.393
-1.065
-20.477
-28.813
-3.426
-0.004
37
0.009
1.605
-4.148
-4,700
-0.563
-1,386
-3.324
0.017
Block
81
82
6g
8g
&
86
87
3g
8g
0.001
0.003
-0.001
-0.002
-0.002
0.006
-0.009
31
0.015
-2.953
-1.905
-0.743
-0.181
1.547
32
1,071
0.652
Block
(5
-2.981
-0.364
0.953
0.615
-1,092
-4,669
34
3.858
-20.660
-0.632
3g
-8.490
-28,409
-1.408
96
-3.195
37
-4,325
B1
8a
81
3g
3x
5.610
3.174
3.324
16.804
27.424
322868
15.323
1.170
Diagonal
0.001
1.910
1.779
1.860
2.152
2,662
2.257
3.406
1.211
Block I
0.000
5.599
3.199
5.361
16.996
27.320
53.149
15.278
1.184
Block 2
0.001
4.434
2.009
3.869
8.517
11.167
4.957
7,798
1.170
ML
0.000
5.609
3.171
5.361
16.893
27.676
33.514
15.283
Table 10: Covariance matrices. Top table: covariance matrix of the estimated coefficients for the full variational posterior and covariances of the ML estimate. Second and third table: block covariance matrices where Block corresponds to case (iii) and Block 2 to case (iv). Bottom table: variances of the full and diagonal posteriors along with ML variances. All the entries across the tables multiplied by 1
24

---

1200
Lower bound
1000
800
Full
Diagonal
Block 1 (1B)
Block 2 (iv)
600
200
400
600
800
1000
1200
Iteration
Figure 8: Lower bound optimization for the Istanbul data under different variational covariance structures.
C.4 DATASETS AND HYPERPARAMETER
Dutuset
Model
Number of Purameters
Number of sumples
Sumples train set
Samples in test set
Period
Labour Volatility
Logistic regression
753
564 (75%)
189 (25%%)
ARCH
2112
1689 (80%)
423 (20%)
3-Jan-2014 009-Jun-2022
GARCHO(10011
2112
1689 (80%))
(20%)
GJR(L.I.I)
2112
1689 (80%)
423 (20%)
GJR(1.1,2)
2112
1689 (80%)
423 (20%)
EGARCH(1O..)
2112
1689 (80%)
423 (20%)
EGARCH(I,LI))
2112
1689 (80%)
423 (20%)
EGARCH(.112)
2112
1689 (80%)
423 (20%)
FIGARCH(,,LD
2112
1689 (80%)
423 (20%))
FIGARCH(L.22
2112
1689 (80%)
423 (20%)
HAR Linear regr.)
4+1
2102
1681 (80%)
421 (20%)
1iee-2014//8-Jumm
Istanbul
Linear regression
8+1
536
428 (80%)
108 (20%S)
5-Jam-2009 22--Feb-2022
EMGVB optimizer
Initial values
Prior
Experiment
8
Grad clip
Grad clip inil.
10
Las
Pi
Hic
Eg
Labour data
0.01
3000
1000
0.4
1200
1000
500
75
(/ E1)
0.05
ARCH GARCH- GIR
0.01
1000
1000
0.4
1200
1000
500
150
ML.
0.05
EGARCH
0.01
1000
1000
0,4
3000
2500
500
ML
0.05
FIGARCH
0.01
1000
1000
0.4
1200
1000
500
150
ML
0.05
HAR Istanbul data
0.001
50000
500
0.4
1200
1000
500
ML
0.01
0.07
50000
500
0.4
1200
1000
500
100
~ N/
E,)
0.01
Table 11: Details on the datasets and corresponding models, as well as the hyperparameters used in the experiments.
Table 11 summarized some information about the datasets and the setup used across the experiments. For the experiments on the Labour and Volatility datasets, the same set of hyperparameters applies to EMGVB, MGVB (and QBVI). While the Labou and Istanbu datasets are readily available, the volatility dataset extracted from the Oxford-Man Institute realized volatility library. We use daily close-to-che to-coee demeaned returns for the MARCH-famile models and 5-minute sub-sampled daily measures of realized volatilities (further annualized) for the HAR model.
Publicly available at Ikey2stats wttp://ttst..oomtataaaeteeeeeeee sesee/viiwww/e See data is also adopted in VI applications e.g. by (Tran et al. 202lb Magris et al.
Mroz 1984 for details. The
2022
PPublicyy available the UCI repository, arch ve. iss vet edu /ataaaeta /istttaiissttttttt
stock+exchanne See (Akbilgic et al. 2014) for details.
Publicly available real 1 zed. oxford man. ox..o.....ukk
25

---

D PROOFS
D.1 PROOF OF PROPOSITION
Preliminaries:
Noticing that for exponential family distributions Ze- Vam and using the chain
rule, Vx((A) "mmVmC(m I vn (()), implying that the natural gradient
=
Vxc(A) - = (m)(m)
(25)
can be easily computed as the euclidean gradient with respect to the expectation parameters, without requiring the inverse FIM (Khan Lin 2017 Khan al. 2018a). Khan et al. (2018a) derive
V m, ((m) (mC(m) 2 [(m))
and
mmz((m) VyL(m),
(26)
which allows expressing the euclidean gradients with the respect the expectation parameters as euclidean gradients with respect to and Σ, thus providing an exact relationship between the natural gradients of the LB and its euclidean gradients with respect to the common (n, Σ) parametrization for yy- Note that the above (and the following) applies to Gaussian distributions only.
The first natural gradient in Proposition Lis trivial as follows from the definition of natural gradient and the Gaussian FIM in eq.
Proof:
If g = É(A) is smooth reparametrization of the variational density,
Ie = go V log qe (0) = JI
with J = VeA being the Jacobian matrix (Lehmann & Casella 1998). If in addition & is an invertible function of X, then // is itself invertible. Therefore for Σ-1 2222, the above implics that
T (Vz 122) 122)-1 Ixd (V 112)
with (Vx=-d22) = Vxa E
Thus for the natural gradient Vz-
0d2 ac 0I-I0N
Vz L = Iz--1 V IC -
-I-- (Vg-1 12)
(Vg-1122 I/' (V (Vg-1d)))V
(Vx [
-22
From eqs. 25 and 26 V C m. L VsL, so that
V - =
which proves the proposition.
D.2 GENERAL FORM OF THE EMGVB UPDATE
For prior p N(o, En) and a variational posterior y Nu, Nu,E), E), by rewriting the LB
Eyg. (0)) = Egg log p(0) log (0) log p(y)()))) Eg (() log (0))+ [log
Egg
decompose VcL. E. (()) log 9. (0) + V F. Plg 00). As in
apply the
log-derivative trick on the last term and write F p(y]0)] = Eg V [log de (0)) log (0)]
On the other hand, it is easy to show that up to a constant that does not depend on j and
Eyg (log p(0) log (0)) log Eol
log [Σ]
+ 7d
1
tr
(00
(u nO)
"no (μ no)-
2
26

---

so that
1 Dー1
1
V Vr:B.ce [log((0) log 2 (0)]
-0
2 2
ne [log (0) log a< (0)]
-o (1r pa).
For the natural gradients, have
V s: Byg p(0) log y5 (g() 22 By ((0 log u5 [0)] ---1 En
4 Egg [lo) log 2c (0)) = EV VV,, Be Dog p(0) log 2< (0)] = -EN (1= (1)
while the feasible naive estimators for V E.g ply and 5 5.-1Eg. ggg log p(y]0)] turn analogous to the right-hand sides of eqs. 9 with he replaced with log p ()). This leads to the general form of the EMGVB update, based either on the h--fuccion gradient estimator (generally applicable) or the above decomposition (applicable under Gaussian prior):
1 D51
pp((G) cp
10. - -14) log [(0s)]
(27)
5
g
((G 2 Cs., +
Σ
[ Ej - (0, ( -+)(0, "=)==,")) log s(0.)
(28)
5
mm1
where
Cr. = = + Σn-
Cp. -EEo (me mm)
if is Gaussian
(0,) = log p(y
(29)
Ch. =0
Cu. = 0 log f0.) = (0.)
if is Gaussian or not
D.3 JUSTIFICATION FOR THE EMGVB UPDATE
For any positive definite matrix S the Riemann gradient VsC, for a differentiable function is SV LS (Hosseini & Sra 2015). This is the form of the Riemann gradient obtained from the SPD (Symmetric and Positive Definite) matrix manifold, for which the following retraction and parallel transport equations for the SPD (matrix) manifold apply:
Rs(E) S ss where TsM.
(30)
T5, FÇET, where E (8++1 E EETSM
(31)
with & being the rescaled Riemann gradient & BSV obtained from the SPD manifold. > 0 rescales the tangent vector and is arbitrary. From algorithmic perspective, B is inter- preted as learning rate, driving the magnitude of the gradient component in the retraction. For the precision matrix Σ-1
V - -
(32)
On the other hand, for the natural gradient
v - -V 1-
where the first equality comes from Proposition Jadd the second is casy to prove with simple matrix algebra and is furthermore analogous to the form o In this regard, more be found in (Barfoot 2020).
The natural gradient is obtained from the Gaussian manifold, that applying the above SPD man- ifold retraction and parallel transport equations is technically incorrect. The concept of retraction R.)
is general and indeed eq. 4 denotes generic retraction function R5 (-), whose specific form is
27

---

specified in Section 4.2 and coincides with eq. 30 The of eq. 30 which is specific for the SPD manifold, with the natural gradient obtained from the Gaussian manifold, appears incorrect. The gradient V IL = -V is not a Riemannian gradient for the SPD manifold but a natural gradient obtained from the Gaussian manifold and equal to 2vxL. The actual exponential map and retraction for updating Σ- based Vv need to be separately worked out for the Gaussian manifold by solving a system of ordinary differential equations, whose curbstone solutions do not coincide with those obtained within the SPD manifold.
In this light, the use of 30 for E and the natural gradient in Tran et al. (022a)) is not well-justified (Lin et al. 2020). Indeed, their approach can be thought of inexact in two ways.
(i) The natural
gradient should read as 22V SSS in place of >V C (eq. 5,4 in (Tran et al. (021a)): the form of retraction and parallel transport therein applied is consistent with the adopted form of the natural gradient (S) which is indeed a Riemann gradient for the SPD manifold, however, the actual natural gradient 2 3 (ii) Perhaps the form g in place of LS is a typo, thus the application of the applied to a Riemann gradient obtained from the SPD manifold but to the natural gradient obtained from the Gaussian manifold, as also discussed in (Lin et al. 2020).
form of the retraction and parallel transport is inexact as it is not
In this view, our discussion in Section 4 is subject to the same inexact setting as the above case (ii) -a correct form for the natural gradient but mixing manifold structures-, under which the adoption of 100 and 3 is, as in Tran et al. (2021a) not justified (though working in practice).
We now show that the forms of retraction and parallel transport in Section 4.2 are justified and arise from the consistent use of the SPD manifold for updating 2Σ-1, from which the update for E-1 follows, and corresponds the retraction from in Section 4.2 Later, we show that parallel transport in eq. 15 also applies.
Consider of updating (n, 2E-1 in place of E-) According eq. 32 the Riemann gradient of 2 is for the SPD manifold:
1 OC VEr- = (2Σ (2E (2 =(2E-1) 2 205-1 (2
2Σ 2-- VVi-i
=V--IC
Thus 2-- = is the Riemann gradient w.r.t. 2Σ-1 obtained from the SPD) manifold. 2Σ can be updated by the retraction in eq. 30 with the Riemann gradient V2 1.. The update is now legit and justified, as it consistently adopts the SPD manifold Riemann gradient and the SPD manifold form of retraction:
2Σ < Ros- (BV = -Rox- (5V=-1))
that is,
2 = 2 8 (5
+ 2 < [
(33)
As 3 > 0 is arbitrary rewrite the above in terms of an arbitrary G¹ 28 (i.e. simple reparametrization of the hyperparameter),
2 = 8 [--
1
(34)
The update for E-1 follows,
8¹
1
4 v [ [5 []
E x, +
Vz L
+
2
2 4
V c
1 3
2
[ 2 [
+
+
2
2
B
A c
1
32
V
(35)
2
Rx-1 ()
Rs (-228v2L)
28

---

where 8' B/2.. From the above, the main equalities are
Ras (5V = Rg-1 (8s-(0VV--1 (DVz-12) Rr-11-22000 R2--1-228Vv.),
whose interpretation is as follows. The retraction 2M-1 with the Riemann gradient 2 2 leads to an update for 2-1 (eq. 35) which the same update obtained with the retraction in 30 on Σ-1 with g-g-t This last gradient analogous to 22 (Proposition and the corre- sponding retraction for is that presented in Section 4.2
Blindly updating E-1 with Rg-1 (--33- is by itself inexact as it involves the natural gradient from the Gaussian manifold, however, this is equivalent to the update for Σ-1 that one obtains by in updating 2 with the consistent retraction for the SPD manifold R2 (3V2x-1) From the above we also have that, logically, the arbitrary stepsize 3 for updating Σ-1 is half of that used for updating 2Σ-1, 23.
similar argument holds for vector transport. Consider the vector transport for 2M-1,
225 -22,1 2V I)) = 2225V22-11))
with F= = (2 (2 -17 -1 = - Ez]
1/2 The vector transport for E-1 is then
T3 -2 (V =L)
T 22, +22,, (V2r L) E(SV 25 1)) FT
(36)
=E((Vs-11) ET
Now note that,
E ( () = T5-1 - (5Vz-1L) = -T - (-2
(37)
Ex
The vector transport in the form of eq. 31 is consistently applied to 2Σ-1 based on the Riemann SPD manifold gradient Vaz-1 L. from which the vector transport for Σ-1 follows (eq. 36). which equals to the vector transport in eq. 31 applied to E-1 based on the rescaled natural
as in Section 4.2
gradient SV 52, eq. 37
D.4 PROOF OF THBOREM 1
For Gaussian distribution distribution with parameter = (n, vec(E)), be VcL(G)
(C,, (ee(V>;((()) where V.L(G:) is the derivative of L({) with respect to æ evaluated at
S G. Furthermore, adopt the following short-hand notation 95 q(u, E) and qç. := q(u, Et).
By using the well-known form of the KL divergence between two multivariate Gaussian distribu- tions, the optimization problem can be written as:
1 3 DKl(9<llg.) = P T V C(S.) + vec(E) vec(V>[(()))
(C. vcF(G))
1
Σ
d 1r+t((E77)) + (m ux) ma)
2 log
J
V (C.)
+ x((EV(((G))
1
Σ, d tr( T + (μ (p--))) p1) ( mm)
log
2,3
Note that V D(G) and VxL(() are now constants and that the Hessian of the above equation amounts to the Hessian of the KL divergence. Furthermore the FIM of &c equals the Hessian of the function 4. + DKL (ge llgc,) evaluated at G, = Ç, that of the above KL divergence:
I V DKL lc..c.
The multivariate Gaussian distribution is a full-rank exponential family for which the negative Hes- sian of the likeliield (FIM) is the covariance matrix of the sufficient statisties. The FIM is
29

---

convex combination of such positive semi-definite matrices, so it positive-definite Thus the above expression is convex with respect to Çt-
The objective is optimized by setting its derivatives with respect to Ç to zero:
vxl(Gr)
=x (1 - 14)) = 0,
8
1 [n - Σー 0, VxL(G))
2
From which it follows that the optimal updates read
μ= μ4 + VE, (. Σ - - 2
The update for j is analogous to the EMGVB update, thus optimal in terms of Theorem ! On the other hand, we have that the optimal update for Σ-1 is the above one, which differs from the EMGVB update. Indeed EMGVB accounts for Dositive-definiteesss constraint on Σ-1, not in the hypotheses of the Theorem. Note however that the EMGVB update for Σ-1 corresponds to the above one when the second-order term the retraction (which indeed accounts for the positive definiteness of 1 is ignored.
30


Reviews:
[{"review_id": "tBi94K7HJRV", "reviewer": "Reviewer_2pY7", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The submission proposed a novel natural gradient method for Gaussian variational inference to guaranttee the positive definite constraint on the variational covariance matrix. The propsed method does not need to compute the inverse of FIM explicitly, so it provides exact update rules. Moreover, the proposed method is validated on different models and compared with baseline methods. ", "summary_of_the_review": "The submission proposed a novel natural gradient method for Gaussian variational inference to guaranttee the positive definite constraint on the variational covariance matrix. To avoid the approximation on the computation of the inverse of FIM, the proposed method targets to update the inverse of covariance matrix whose natural gradietn is available in an exact form. Moreover, to facilitate the compuation of gradient of LB, the method also utilizes the log-derivative trick. The method is interesting in my eyes and I am inclined to acceptance, but I would also like to see the opinions of other reviewers as I am not an expert in natural gradient. ", "strength_and_weaknesses": "Stength: The submission provided a natural gradient method for Gaussian variational inference to guaranttee the positive definite constraint on covariance matrix. \n\nWeakness: I did not find any weakness in the submission. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "rW_VCntg2L", "reviewer": "Reviewer_VZUh", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors provide a variational inference optimizer on the manifold of positive-definite matrices with all its implementation details and several experiments with results that are, in a wide set of experiments, surprisingly analogous to sampling methods.  ", "summary_of_the_review": "I would first like to say that I am not an expert on the variational inference field but of the Bayesian optimization field, so some of my arguments here may be invalid or obsolete. The paper presents a new VI optimizer whose results seem to be the same as some sampling methods, thing that I find very suprising and that if it is validated by VI/MCMC experts it would be very great! The explanation, clarity and quality of the method is also great, so I believe that this can be a great paper.\n\nSome minor things that I can say about the paper to make it better:\n\nIn the introduction section, I miss a paragraph that describes the organization of the paper.\nIn the second section, we find the claim that \"Through the problem can be tackled with sampling methods, in high-dimensional apps such as DL, MC techniques are challenging and unfeasible\" without references. This may be controversial. In particular, some Bayesian DL methodologies use MCMC sampling and the same happens with Deep Gaussian processes, where I think that this can also be applied successfully BTW (this may be included as a future research line). This needs to be better supported by arguments as it may \"bother\" some researchers :-)\nProvide a reference of the LB maximizers ADAM and company as well as to the beautiful topic of the information geometry of distributions.\nThe PSD manifold may be better defined as the whole paper depend on this concept, at least the PSD matrix, it may help some novel readers.\nIMHO, the related work section must be number 2 to give more cohesion to the paper.\nClarify whether the product operator after equation 6 is a Kronecker product. \n\"For coherency with the literature available on variational inference and BDL....\" provide here a reference. \nProvide the complexity of algorithm 1.\nIn the experiments section: how many repetitions? And the most important thing: provide a link to the experiments.\n\nThanks for your work and to the other reviewers please validate the method section.", "strength_and_weaknesses": "Strength: Rigurous, state-of-the-art results, style.\nWeaknesses: I find the experimental section lack some details in the main manuscript. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "guGZ1tyduz", "reviewer": "Reviewer_8VC5", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this work, the authors consider an \"exact\" version of Tran et al 2021a in gradient-free settings via the reinforce gradient estimation.\nSeveral standard manifold techniques such as retraction and vector transport are included in the proposed method.\nThe authors claim that Eq (2) is the natural-gradient w.r.t. $\\color{red} \\Sigma$ is not exact and propose an exact natural-gradient w.r.t. $\\color{red}  \\Sigma^{-1}$ in Proposition 1.\nUnfortunately, there are many incorrect and misleading statements. In fact, both natural-gradient (Eq (2)) w.r.t. $\\color{red} \\Sigma$  and natural-gradient (Proposition 1) w.r.t. $\\color{red} \\Sigma^{-1}$  are exact/unbiased (see pointer 2.1 in the weakness section). Thus, there is no such \"exact\" version of  Tran et al 2021a. \n\n", "summary_of_the_review": "The main claim of this work: the authors claim that Eq (2) (Tran et al 2021a) is an inexact/biased  natural-gradient w.r.t. $\\color{red} \\Sigma$, and propose an exact/unbiased natural-gradient w.r.t. $\\color{red}  \\Sigma^{-1}$ in Proposition 1.\n\n1. Unfortunately, there are many incorrect and misleading statements.  Indeed, Eq (2) (Tran et al 2021a) is the exact/unbiased natural-gradient (obtained from the Gaussian manifold) w.r.t. $\\color{red} \\Sigma$. The main claim is incorrect.\n2. If the authors want to give an \"exact\" version of Tran et al 2021a, **the authors should give an exact version of natural-gradient w.r.t.  $\\color{red} \\Sigma$ instead of the one w.r.t. $\\color{red}  \\Sigma^{-1}$**. Note that natural-gradient decent depends on the choice of parametrization. Although a natural gradient can be non-linearly invariant,  a linear update in the natural gradient (A.K.A., natural-gradient decent) is only linearly invariant (see https://informationgeometryml.github.io/posts/2021/11/Geomopt04//#natural-gradient-descent-is-linearly-invariant )\nThus, the performance of natural-gradient descent depends on the choice of parametrization (e.g., the choice to update either $\\color{red} \\Sigma$ or $\\color{red} \\Sigma^{-1}$). In other words, the improvement in the experiments is mainly due to (1) adding the momentum and (2) the change of parametrization.  \n\n3. The title, **Exact** manifold **Gaussian** Variational Bayes, is misleading/confusing. It is unclear why the authors use the retraction map and the transport map for the **matrix manifold** instead of the ones for the Gaussian manifold. Moreover, the retraction map is **inexact** even in the full positive definite matrix manifold cases (see Table 1 of [1] for the exact retraction/exponential map).\n\n ", "strength_and_weaknesses": "* Strength\n1. Study the method of Tran et al 2021a for black-box variational inference  \n2. introduce a momentum via the vector transport \n\n* Weakness \n There are many incorrect and misleading statements.\n---\n1. Ad-hoc usage of the manifold structures. There are two **distinct** types of manifold structures.\nTechnically speaking, the authors should stick to using one manifold structure at a time. \n\n*  **Matrix Manifold**: Riemannian manifold structure for the positive definite covariance/precision matrix  (the authors use the retraction and the vector transport obtained from the matrix manifold)\n*  **Gaussian Manifold**: Riemannian manifold structure for the Gaussian distribution    (the authors use natural/Riemannian gradients  obtained from the Gaussian manifold) \n\nEssentially, the authors compute natural/Riemannian gradients derived from the Gaussian manifold while using the retraction and the vector transport derived from the positive definite matrix manifold.\nAs discussed by Lin et al 2020., the standard Riemannian metric for the positive definite matrix manifold is **different** from the standard Riemannian metric  (A.K.A., the Fisher information matrix) for the Gaussian manifold.\n\nFor any positive definite matrix $P$, the Riemannian gradient (obtained from the matrix manifold) is always $P \\nabla_P (\\mathcal{L}) P$ as shown in Table 1 of [1].\nAccording to the matrix calculus, we have this identity: $P \\nabla_P (\\mathcal{L})   P= {\\color{red}-  } \\nabla_{P^{-1}} (\\mathcal{L}) $\n\n* difference 1: For the precision matrix $\\Sigma^{-1}$ (when $P=\\Sigma^{-1}$),  the Riemannian gradient (obtained from the matrix manifold) is $ {\\color{red}-  } \\nabla_{\\Sigma} (\\mathcal{L})$, while the natural gradient (obtained from the Gaussian manifold) w.r.t.  $\\Sigma^{-1}$  is $ {\\color{red}- 2 } \\nabla_{\\Sigma} (\\mathcal{L})$. See proposition 1 in the main text, where the authors claim this is **exact**/unbiased.\n\n* difference 2: Similarly,  for the covariance matrix $\\Sigma$ (when $P=\\Sigma$),  the Riemannian gradient (obtained from the  matrix manifold) is $ {\\color{red}-  } \\nabla_{\\Sigma^{-1}} (\\mathcal{L})$, while the natural gradient (obtained from the Gaussian manifold) w.r.t.  $\\Sigma$  is $ {\\color{red}- 2 } \\nabla_{\\Sigma^{-1}} (\\mathcal{L})=2 \\Sigma \\nabla_{\\Sigma} (\\mathcal{L})   \\Sigma $. See Eq (2) in the main text, where the authors claim this is **inexact**/biased.\n\n* difference 3: The retraction map and the vector transport map for the Gaussian manifold are different from the ones for the matrix manifold.  Why do the authors want to use the two maps (in Sec 3.2) obtained from the matrix manifold, instead? Any justification?\n\n---\n2. The natural gradient   (obtained from the Gaussian manifold)  for $\\Sigma$ in Eq (2) is indeed exact/unbiased. Thus, the main claim of this submission is incorrect.\n\nThere are some caveats when it comes to the Fisher matrix computation as discussed in this blog, https://informationgeometryml.github.io/posts/2021/09/Geomopt01/#caveats-of-the-fisher-matrix-computation\n\n2.1 Technically speaking (as discussed in [2]),  the Fisher information matrix $\\mathcal{I}(\\Sigma)$ for $\\Sigma$ is singular if we consider $\\Sigma$ to be a $d$-by-$d$ matrix. Note that by the definition of Gaussian, we have to **explicitly enforce the symmetric constraint** in $\\Sigma$ as ${\\color{red} \\text{highlighted in red}}$ to compute the Fisher information matrix:  $\\mathcal{I}(\\Sigma) = - E_{p(z)}[ \\nabla_\\Sigma^2  \\log p(z|\\mu,\\Sigma)] $, where $\\log p(z|\\mu,\\Sigma)=-\\frac{1}{2} [d\\log (2\\pi) + \\log \\mathrm{det}({ \\color{red} \\frac{\\Sigma+\\Sigma^T}{2} } ) + (z-\\mu)^T ( { \\color{red} \\frac{\\Sigma+\\Sigma^T}{2} }  )^{-1} (z-\\mu) ] $.\nIn other words,  **the Fisher information matrix for $\\mathrm{vec}(\\Sigma)$ is singular**.  I encourage the authors to work on a 2-dimensional Gaussian case for verification. In this singular case, the **Moore-Penrose inverse** of the Fisher information matrix  (where the authors refer to this inverse as **an inexact approximation**) should be used as discussed in [5]. \n\n2.2 A more involved but equivalent way is to consider **the Fisher information matrix for ${\\color{red}\\mathrm{vech}}(\\Sigma)$**, where ${\\color{red}\\mathrm{vech}}(\\Sigma)$ vectorizes symmetric matrix  $\\Sigma$ to a $\\frac{d(d+1)}{2}$-by-$1$  vector obtained by only the lower triangular part. **The Fisher information matrix w.r.t.  **${\\color{red}\\mathrm{vech}}(\\Sigma)$** is non-singular**.\n \nThe authors can also find related discussions in [3].\nIn summary,  in many existing works ([3-4]) ,  the (**exact**) natural gradient  (obtained from the Gaussian manifold)  w.r.t. $\\Sigma$ is exactly Eq 2, which can be computed by either using the the Moore-Penrose inverse for $\\mathrm{vec}(\\Sigma)$ or using ${\\color{red}\\mathrm{vech}^{-1}}(\\mathbf{a})$,  where $\\mathbf{a}$ is a natural gradient w.r.t.  ${\\color{red}\\mathrm{vech}}(\\Sigma)$,\n\nNote that: **the Fisher information matrix for $\\mathrm{vec}(\\Sigma^{-1})$  is indeed also singular** while the one for $\\mathrm{vech}(\\Sigma^{-1})$ is not.\n\n---\n3. The retraction and vector transport in Sec 3.2 are not exact and only for the full-rank and diagonal matrix manifold instead of the Gaussian manifold.\n\nThese two maps are **only valid** for full-rank/dense and diagonal positive definitive matrices.\nIt is an open problem to find these two maps for (1) the Gaussian manifold and (2) a structured positive definitive matrix manifold (e.g., rank-1 plus diagonal positive definitive matrices).\n\nThe authors should note that the retraction and vector transport are often **approximations** of the (exact) manifold exponential map and the (exact) parallel transport map. \n The exponential map and the parallel transport map are defined by solving a system of matrix ODEs, which is hard to obtain a closed-form expression both in the Gaussian case and structured matrix case.\n\n---\nReferences\n* [1] Hosseini, Reshad, and Suvrit Sra. \"Matrix manifold optimization for Gaussian mixtures.\" Advances in Neural Information Processing Systems 28 (2015).\n* [2] Lin, W., Nielsen, F., Khan, M. E., & Schmidt, M., Introduction to Natural-gradient Descent,  https://informationgeometryml.github.io/posts/2021/09/Geomopt01/#dimensionality-of-a-manifold, https://informationgeometryml.github.io/posts/2021/10/Geomopt02/#definition-of-natural-gradients, and https://informationgeometryml.github.io/posts/2021/12/Geomopt05/#efficient-ngd-for-multivariate-gaussians\n* [3] Barfoot, Timothy D. \"Multivariate Gaussian variational inference by natural gradient descent.\" arXiv preprint arXiv:2001.10025 (2020).\n* [4] Zhang, Guodong, et al. \"Noisy natural gradient as variational inference.\" International Conference on Machine Learning. PMLR, 2018.\n* [5] van Oostrum, Jesse, Johannes Müller, and Nihat Ay. \"Invariance properties of the natural gradient in overparametrised systems.\" Information Geometry (2022): 1-17.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "cOL4aYVMBzr", "reviewer": "Reviewer_y1Je", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors proposed a new variational inference algorithm for settings where the parameters of the approximate posterior form a Riemannian manifold. Unlike previous methods that relied on various forms approximations, the proposed method is based on an exact formula for the gradient.", "summary_of_the_review": "I think this paper is currently below the threshold of acceptance, but could be over the threshold after an improved presentation and a more thorough exploration of experimental results.", "strength_and_weaknesses": "STRENGTHS:\n* The paper provides a novel estimator for the gradient in a broad class of models fit with variational inference.\n* The derivations and the proposed algorithms are novel, technically sound, and non-trivial. The level of technical depth required to derive the estimator is significant.\n\nWEAKNESSES:\n* The experiments do not seem to convey a strong sense of the improvements offered by the new methods. The performance improvements seem marginal at first glance.\n* The writing quality and the overall structure of the paper is acceptable, but not quite at the level of polish and organization seen in a typical ICLR accepted paper.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject