Variational Bayes on Manifolds

TL;DR

This paper extends Variational Bayes to Riemannian manifolds, developing efficient algorithms that leverage geometric structures, with proven convergence and superior performance in numerical experiments.

Contribution

It introduces manifold-based VB algorithms that incorporate geometric and information-theoretic structures, broadening VB applicability beyond Euclidean spaces.

Findings

Algorithms are provably convergent with specific rates.

Manifold VB algorithms are stable and less sensitive to initialization.

Numerical experiments show improved performance over existing VB methods.

Abstract

Variational Bayes (VB) has become a widely-used tool for Bayesian inference in statistics and machine learning. Nonetheless, the development of the existing VB algorithms is so far generally restricted to the case where the variational parameter space is Euclidean, which hinders the potential broad application of VB methods. This paper extends the scope of VB to the case where the variational parameter space is a Riemannian manifold. We develop an efficient manifold-based VB algorithm that exploits both the geometric structure of the constraint parameter space and the information geometry of the manifold of VB approximating probability distributions. Our algorithm is provably convergent and achieves a convergence rate of order $\mathcal O(1/\sqrt{T})$ and $\mathcal O(1/T^{2-2\epsilon})$ for a non-convex evidence lower bound function and a strongly retraction-convex evidence lower bound function, respectively. We develop in particular two manifold VB algorithms, Manifold Gaussian VB and Manifold Neural Net VB, and demonstrate through numerical experiments that the proposed algorithms are stable, less sensitive to initialization and compares favourably to existing VB methods.

Tables

Table 1. Table 1: The numbers in brackets show the posterior variances. The manifold Wishart Variational Bayes stopped after 75 iterations, with S = 1000 𝑆 1000 S=1000 Monte Carlo samples used for estimating the score-function gradient of the lower bound. The gradients in ( 5.7 ) and ( 5.8 ) were used.

$v_{ij}$	True	VB	$v_{ij}$	True	VB	$v_{ij}$	True	VB
$v_{11}$	0.89 (0.03)	0.91 (0.04)	$v_{21}$	$-0.37\left(0.02\right)$	$-0.37\left(0.02\right)$	$v_{31}$	0.39 (0.02)	0.39 (0.02)
$v_{41}$	$-0.32\left(0.02\right)$	$-0.33\left(0.03\right)$	$v_{51}$	0.25 (0.02)	0.26 (0.02)	$v_{22}$	0.83 (0.03)	0.84 (0.03)
$v_{32}$	$-0.39\left(0.02\right)$	$-0.39\left(0.02\right)$	$v_{42}$	0.29 (0.02)	0.30 (0.02)	$v_{52}$	$-0.32\left(0.02\right)$	$-0.32\left(0.02\right)$
$v_{33}$	0.93 (0.04)	0.94 (0.04)	$v_{43}$	$-0.50\left(0.03\right)$	$-0.50\left(0.03\right)$	$v_{53}$	0.35 (0.02)	0.35 (0.02)
$v_{44}$	1.18 (0.06)	1.19 (0.07)	$v_{54}$	$-0.73\left(0.04\right)$	$-0.74\left(0.04\right)$	$v_{55}$	0.96 (0.04)	0.97 (0.04)

Equations

$$p(\theta |y)=\frac{p(\theta )p(y|\theta )}{p(y)},$$

$${\lambda }^{\ast }=\arg \underset{\lambda \in \mathcal{M}}{\min }\left\{\text{KL}(q_{\lambda }\Vert p(\cdot |y))=\int q_{\lambda }(\theta )\log \frac{q_{\lambda }(\theta )}{p(\theta |y)}d\theta \right\}.$$

$$\text{KL}(q_{\lambda }\Vert p(\cdot |y))=-\int q_{\lambda }(\theta )\log \frac{p(\theta )p(y|\theta )}{q_{\lambda }(\theta )}d\theta +\log p(y),$$

$$\mathcal{L}(\lambda )=\int q_{\lambda }(\theta )\log \frac{p(\theta )p(y|\theta )}{q_{\lambda }(\theta )}d\theta .$$

$$\mathcal{N}=\{q_{\lambda }(\theta ):\lambda \in \mathcal{M}\}$$

$$I_F(\lambda ):={\mathrm{Cov}}_{q_{\lambda }}({\nabla }_{\lambda }\log q_{\lambda }(\theta ))={\mathbb{E}}_{q_{\lambda }}[{\nabla }_{\lambda }\log q_{\lambda }(\theta )({\nabla }_{\lambda }\log q_{\lambda }(\theta ){)}^{\top }]$$

$$\text{KL}(q_{\lambda }||q_{\lambda +ϵ})$$

$$=-{\mathbb{E}}_{q_{\lambda }}({\nabla }_{\lambda }\log q_{\lambda }(\theta ){)}^{\top }ϵ+\frac{1}{2}{ϵ}^{\top }{I}_F(\lambda )ϵ$$

$$=\frac{1}{2}{ϵ}^{\top }{I}_F(\lambda )ϵ.$$

$$\mathcal{L}(\lambda +ϵ)\approx \mathcal{L}(\lambda )+{\nabla }_{\lambda }\mathcal{L}(\lambda {)}^{\top }ϵ.$$

$$\operatorname*{arg\,max}_{\epsilon:\epsilon^{\top}I_{F}(\lambda)\epsilon=l}\Big{\{}\nabla_{\lambda}\mathcal{L}(\lambda)^{\top}\epsilon\Big{\}}.$$

$$ϵ={\nabla }_{\lambda }^{\text{nat}}\mathcal{L}(\lambda ):=I_F^{-1}(\lambda ){\nabla }_{\lambda }\mathcal{L}(\lambda ).$$

$$\Sigma =B{B}^{\top }+D^2=B^{\prime }{B^{\prime }}^{\top }+D^2$$

$$[B]=\{BO:O{O}^{\top }=I_p\}.$$

$$\underset{\lambda \in \mathcal{M}}{\mathrm{argmax}}\mathcal{L}(\lambda ).$$

$$<{\zeta }_{\lambda },{\xi }_{\lambda }>={\zeta }_{\lambda }^{\top }{I}_F(\lambda ){\xi }_{\lambda }.$$

$$\underset{{\eta }_{\lambda }\in T_{\lambda }\overline{\mathcal{M}},\Vert {\eta }_{\lambda }\Vert =1}{\mathrm{argmax}}\text{D}\overline{\mathcal{L}}(\lambda )[{\eta }_{\lambda }],$$

$${\nabla }_{\lambda }^{\text{nat}}\overline{\mathcal{L}}(\lambda )=I_F^{-1}(\lambda ){\nabla }_{\lambda }\overline{\mathcal{L}}(\lambda ),\lambda \in \overline{\mathcal{M}}.$$

$${\nabla }_{\lambda }^{\text{nat}}\overline{\mathcal{L}}(\lambda )=I_F^{-1}(\lambda ){\nabla }_{\lambda }\overline{\mathcal{L}}(\lambda )=I_F^{-1}(\lambda ){\nabla }_{\lambda }\mathcal{L}(\lambda )={\nabla }_{\lambda }^{\text{nat}}\mathcal{L}(\lambda ).$$

$$<\text{grad}f(\lambda ),{\xi }_{\lambda }>=\text{D}f(\lambda )[{\xi }_{\lambda }],\forall {\xi }_{\lambda }\in T_{\lambda }\mathcal{M}.$$

$$<{\zeta }_{\lambda },{\xi }_{\lambda }>={\zeta }_{\lambda }^{\top }{I}_F(\lambda ){\xi }_{\lambda },$$

$$\text{grad}\overline{\mathcal{L}}(\lambda )={\text{Proj}}_{\lambda }\text{grad}\overline{\mathcal{L}}(\lambda )+{\text{Proj}}_{\lambda }^{\perp }\text{grad}\overline{\mathcal{L}}(\lambda )$$

$$\text{grad}\mathcal{L}(\lambda )={\text{Proj}}_{\lambda }\text{grad}\overline{\mathcal{L}}(\lambda ).$$

$$<\text{grad}\mathcal{L}(\lambda ),{\eta }_{\lambda }>$$

$$=<\text{grad}\overline{\mathcal{L}}(\lambda ),{\eta }_{\lambda }>$$

$$=\text{D}\overline{\mathcal{L}}(\lambda )[{\eta }_{\lambda }]=\text{D}\mathcal{L}(\lambda )[{\eta }_{\lambda }],\forall {\eta }_{\lambda }\in T_{\lambda }\mathcal{M}.$$

$$\mathcal{S}(p,n)=\{W\in \text{Mat}(n,p):W^{\top }W=I_p\},$$

$$<\zeta_{W},\xi_{W}>=\big{(}\text{vec}(\zeta_{W})\big{)}^{\top}I_{F}(W)\text{vec}(\xi_{W}),\;\;\zeta_{W},\xi_{W}\in T_{W}\overline{\mathcal{M}}\cong\text{Mat}(n,p),$$

$$\text{grad}\overline{\mathcal{L}}(W)=\text{vec}^{-1}\Big{(}I_{F}(W)^{-1}\nabla_{\text{vec}(W)}\mathcal{L}(W)\Big{)}\in\text{Mat}(n,p)$$

$$T_W\mathcal{M}=\{Z\in \text{Mat}(n,p):Z^{\top }W+W^{\top }Z=0_{p\times p}\}.$$

$$<\zeta_{W},\xi_{W}>_{Euc}=\text{trace}(\zeta_{W}^{\top}\xi_{W})=\big{(}\text{vec}(\zeta_{W})\big{)}^{\top}\text{vec}(\xi_{W}),\;\;\zeta_{W},\xi_{W}\in T_{W}\mathcal{M}.$$

$$\text{grad}\mathcal{L}(W)=(I_n-W{W}^{\top })G+W\text{skew}(W^{\top }G)$$

Full text

Keywords. Marginal likelihood, variational Bayes, natural gradient, stochastic approximation, Riemannian manifold.

1 Introduction

Increasingly complicated models in modern statistics and machine learning have called for more efficient Bayesian estimation methods. Of the Bayesian tools, Variational Bayes (VB) (Waterhouse et al.,, 1996; Jordan et al.,, 1999) stands out as one of the most versatile alternatives to conventional Monte Carlo methods for statistical inference in complicated models. VB approximates the posterior probability distribution by a member from a family of tractable distributions indexed by variational parameters $\lambda$ belonging to a parameter space $\mathcal{M}$. The best member is found by minimizing the Kullback-Leibler divergence from the candidate member to the posterior. VB methods have found their application in a wide range of problems including variational autoencoders (Kingma and Welling,, 2013), text analysis (Hoffman et al.,, 2013), Bayesian synthetic likelihood (Ong et al., 2018a, ), deep neural nets (Tran et al.,, 2019), to name but a few.

Most of the existing VB methods work with cases where the variational parameter space $\mathcal{M}$ is (a subset of) the Euclidean space $\mathbb{R}^{d}$. This paper considers the VB problem where $\mathcal{M}$ is a Riemannian manifold, which naturally arises in many modern applications. For example, in Gaussian VB where the VB approximating distribution is a multivariate Gaussian with mean $\mu$ and covariance $\Sigma$, $\lambda=(\mu,\Sigma)$ belongs to the product manifold $\mathcal{M}=\mathcal{M}_{1}\otimes\mathcal{M}_{2}$ where $\mathcal{M}_{1}$ is an Euclidean manifold and $\mathcal{M}_{2}$ is the manifold of symmetric and positive definite matrices. We develop manifold-based VB algorithms that cast Euclidean-based constrained VB problems as manifold-based unconstrained optimization problems under which the solution can be efficiently found by exploiting the geometric structure of the constraints. Optimization algorithms that work on manifolds often enjoy better numerical properties. See the monograph of Absil et al., (2009) for recent advances.

Many Euclidean-based VB methods employ (Euclidean) stochastic gradient decent (SGD) for solving the required optimization problem, and it is well-known that the natural gradient (Amari,, 1998) is of major importance in SGD. The natural gradient, a geometric object itself, takes into account the information geometry of the family of approximating distributions to help stabilize and speed up the updating procedure. For a comprehensive review and recent development of the natural gradient descent in Euclidean spaces, the reader is referred to Martens, (2020). Extending natural gradient decent for use in Riemannian stochastic gradient decent is a non-trivial task and of interest in many VB problems. This paper develops a mathematically formal framework for incorporating the natural gradient into manifold-based VB algorithms. The contributions of this paper are threefold:

• •

We develop a doubly geometry-informed VB algorithm that exploits both the geometric structure of the manifold constraints of the variational parameter space, and the information geometry of the manifold of the approximating family, which leads to a highly efficient VB algorithm for Bayesian inference in complicated models.

• •

The proposed manifold VB algorithm is provably convergent and achieves a convergence rate of order $\mathcal{O}(1/\sqrt{T})$ and $\mathcal{O}(1/T^{2-2\epsilon})$, with $\epsilon\in(0,1)$ and $T$ the number of iterations, for a non-convex lower bound function and a strongly retraction-convex lower bound function, respectively.

• •

We develop in detail a Manifold Gaussian VB algorithm and a Manifold Wishart VB algorithm, both can be used as a general estimation method for Bayesian inference. The numerical experiments demonstrate that these manifold VB algorithms work efficiently, are more stable and less sensitive to initialization as compared to some existing VB algorithms in the literature. We would like to emphasize that making VB more stable and less initialization-sensitive is of major importance in the current VB literature. We also apply our VB method to estimating a financial time series model and demonstrate its high accuracy in comparison with a “gold standard” Sequential Monte Carlo method.

The paper is organized as follows. Section 2 reviews the VB method on Euclidean spaces and sets up notations. Section 3 develops the manifold-based VB algorithm and Section 4 studies its convergence properties. Section 5 presents the Manifold Gaussian VB and the Manifold Wishart VB algorithms, and their applications. Section 6 concludes. The technical proofs are presented in the Appendix.

2 VB algorithms on Euclidean spaces

This section gives a brief overview of VB methods where the variational parameter $\lambda$ lies in (a subset of) a Euclidean space. It also gives the definition of the natural gradient, and the motivation for extending the Euclidean-based VB problem into manifolds.

Let $y$ be the data and $p(y|\theta)$ the likelihood function based on a postulated model, with $\theta$ the set of model parameters to be estimated. Let $p(\theta)$ be the prior. Bayesian inference requires computing expectations with respect to the posterior distribution whose density (with respect to some reference measure such as the Lebesgue measure) is

[TABLE]

where $p(y)=\int p(\theta)p(y|\theta)d\theta$, called the marginal likelihood. It is often difficult to compute such expectations, partly because the density $p(\theta|y)$ itself is intractable as the normalizing constant $p(y)$ is often unknown. For simple models, Bayesian inference is often performed using Markov Chain Monte Carlo (MCMC), which estimates expectations w.r.t. $p(\theta|y)$ by sampling from it. For models where $\theta$ is high dimensional or has a complicated structure, MCMC methods in their current development are either not applicable or very time consuming. In the latter case, VB is often an attractive alternative to MCMC. VB approximates the posterior $p(\theta|y)$ by a probability distribution with density $q_{\lambda}(\theta)$, $\lambda\in\mathcal{M}$ - the variational parameter space, belonging to some tractable family of distributions such as Gaussian. The best $\lambda$ is found by minimizing the Kullback-Leibler (KL) divergence from $q_{\lambda}(\theta)$ to $p(\theta|y)$

[TABLE]

It is easy to see that

[TABLE]

thus minimizing KL is equivalent to maximizing the lower bound on $\log\;p(y)$

[TABLE]

SGD techniques are often employed to solve this optimization problem. The VB approximating distribution $q_{\lambda}(\theta)$ with the optimized $\lambda$ is then used for Bayesian inference. See Tran et al., (2021) for an accessible tutorial introduction to VB.

Let

[TABLE]

be the set of VB approximating probability distributions parameterized by $\lambda$, and

[TABLE]

be the Fisher information matrix w.r.t. $q_{\lambda}$. By the Taylor expansion, we have that

[TABLE]

This shows that the local KL divergence around the point $q_{\lambda}\in\mathcal{N}$ is characterized by the Fisher matrix $I_{F}(\lambda)$. Formally, $\mathcal{N}$ can be made into a Riemannian manifold with the Riemannian metric induced by the Fisher information matrix (Rao,, 1945; Amari,, 1998).

Assume that the objective function $\mathcal{L}$ is smooth enough, then

[TABLE]

The steepest ascent direction $\epsilon$ for maximizing $\mathcal{L}(\lambda+\epsilon)$ among all the directions with a fixed length $\|\epsilon\|:=\epsilon^{\top}I_{F}(\lambda)\epsilon=l$ is

[TABLE]

By the method of Lagrangian multipliers, this steepest ascent is

[TABLE]

Amari, (1998) termed this the natural gradient and popularized it in machine learning. In the statistics literature, the steepest ascent in the form (2.5) has been used for a long time and is often known as Fisher’s scoring in the context of maximum likelihood estimation (see, e.g., Longford, (1987)). We adopt the term natural gradient in this paper. The efficiency of the natural gradient over the ordinary gradient has been well documented (Sato,, 2001; Hoffman et al.,, 2013; Tran et al.,, 2017). A remarkable property of the natural gradient is that is is invariant under parameterization (Martens,, 2020), i.e. it is coordinate-free and an intrinsic geometric object. This further motivates the use of natural gradient in optimization on manifolds.

Most of the VB methods and natural gradient descent are developed for cases where the variational parameter $\lambda$ lies in an unconstrained Euclidean space. In many situations, however, $\lambda$ belongs to a non-linear constrained space that forms a differential manifold. A popular example is Gaussian VB where the covariance matrix $\Sigma$ of size $d\times d$ is subject to the symmetric and positive definite constraint. Ong et al., 2018b avoid the difficulty of dealing with this constraint by using a factor decomposition $\Sigma=BB^{\top}+D^{2}$, where $B$ a full-rank matrix of size $d\times p$ with $p\leq d$, and $D$ a diagonal matrix. Such a decomposition is invariant under orthogonal transformations of $B$, i.e.

[TABLE]

for all $B^{\prime}=BO$ with $O$ an orthogonal matrix, i.e. $OO^{\top}=I_{p}$. That is, the variational parameter $B$ lies in a quotient manifold where each point in this manifold is an equivalence class

[TABLE]

This manifold structure is not considered in Ong et al., 2018b . Zhou et al., (2020) take into account this manifold structure and report some improvement over the plain VB methods. Another example is Wishart VB where the VB distribution $q_{\lambda}(\theta)$ is an inverse-Wishart distribution $IW(\nu,\Sigma)$. Here, the variational parameter $\Sigma$ lives in the manifold of the symmetric and positive definite matrices.

Related work

As we employ the SGD method for optimizing the lower bound $\mathcal{L}(\lambda)$, our paper is related to the recent development of SGD algorithms on Riemannian manifolds. Bonnabel, (2013) is one of the first to develop SGD where the cost function is defined on a Riemannian manifold. It is showed in his paper that under some suitable conditions the Reimannian SGD algorithm converges to a critical point of the cost function. In a recent paper, Kasai et al., (2019) propose an adaptive SGD on Riemannian manifolds, which uses different learning rates for different coordinates. Their method is proved to converge to a critical point of the cost function at a rate $\mathcal{O}(\log(T)/\sqrt{T})$. For a recent discussion of generalization of Euclidean adaptive SGD algorithms, such as Adam and Adagrad, to Riemannian manifolds, see Bécigneul and Ganea, (2018). The monograph of Absil et al., (2009) provides an excellent account of recent development on optimization on matrix manifolds. Companion user-friendly software such as Manopt (Boumal et al.,, 2014) has been developed to assist fast growing research in Riemannian optimization.

The natural gradient has been widely used in the machine learning literature; see, e.g., Sato, (2001), Hoffman et al., (2013), Khan and Lin, (2017) and Martens, (2020). However, most of the existing work only consider cases where the variational parameter $\lambda$ belongs to an Euclidean space. Zhou et al., (2020) is the only paper that we are aware of develops a VB method on manifolds. However, their paper only considers the Factor Gaussian VB for the particular quotient manifold in (2.6), and does not consider natural gradient. Our paper develops a general VB method for Riemannian manifolds that incorporates the natural gradient, and provides a careful convergence analysis.

3 VB on manifolds with the natural gradient

This section presents our proposed VB algorithm on manifolds. Recall that we are interested in a VB problem where the variational parameter $\lambda$ lies in a Riemannian manifold $\mathcal{M}$, i.e. we wish to solve the following optimization problem

[TABLE]

In order to incorporate the natural gradient into Riemannian SGD, we view the manifold $\mathcal{M}$ as embedded in a Riemannian manifold $\overline{\mathcal{M}}\subset\mathbb{R}^{d}$, whose Riemannian metric is defined by the Fisher information matrix $I_{F}(\lambda)$. Let $T_{\lambda}\overline{\mathcal{M}}$ be the tangent space to $\overline{\mathcal{M}}$ at $\lambda\in\overline{\mathcal{M}}$. The inner product between two tangent vectors $\zeta_{\lambda},\xi_{\lambda}\in T_{\lambda}\overline{\mathcal{M}}$ is defined as

[TABLE]

For VB on manifolds without using the natural gradient, this inner product is the usual Euclidean metric, i.e. $<\zeta_{\lambda},\xi_{\lambda}>=\zeta_{\lambda}^{\top}\xi_{\lambda}$. The metric in (3.1) is often referred to as the Fisher-Rao metric. Let $\overline{\mathcal{L}}$ be a differentiable function defined on $\overline{\mathcal{M}}$ such that its restriction on $\mathcal{M}$ is the lower bound $\mathcal{L}$. Similar to (2.4)-(2.5), it can be shown that the steepest ascent direction at $\lambda\in\overline{\mathcal{M}}$ for optimizing the objective function $\overline{\mathcal{L}}(\lambda)$, i.e. the direction of

[TABLE]

is the natural gradient

[TABLE]

Here, $\text{D}\overline{\mathcal{L}}(\lambda)[\eta_{\lambda}]$ denotes the directional derivative of $\overline{\mathcal{L}}$ at $\lambda$ in the direction of $\eta_{\lambda}$, and $\nabla_{\lambda}\overline{\mathcal{L}}(\lambda)$ is the usual Euclidean gradient vector of $\overline{\mathcal{L}}(\lambda)$. We note that, for $\lambda\in\mathcal{M}$,

[TABLE]

We recall that the Riemannian gradient of a smooth function $f(\lambda)$ on a Riemannian manifold $\mathcal{M}$, embedded in $\mathbb{R}^{d}$ and equipped with the Riemannian metric $<,>$, is the unique vector $\text{grad}f(\lambda)$ in the tangent space $T_{\lambda}\mathcal{M}$ at $\lambda\in\mathcal{M}$ such that

[TABLE]

The following lemma is important for the purpose of this paper. It shows that the natural gradient $\nabla_{\lambda}^{\text{nat}}\overline{\mathcal{L}}(\lambda)$ is a Riemannian gradient defined in the ambient manifold $\overline{\mathcal{M}}$, which leads to a formal framework for associating the natural gradient to the Riemannian gradient of the lower bound $\mathcal{L}$ defined in the manifold $\mathcal{M}$.

Lemma 3.1**.**

The natural gradient of the function $\overline{\mathcal{L}}$ on the Riemannian manifold $\overline{\mathcal{M}}$ with the Fisher-Rao metric (3.1) is the Riemannian gradient of $\overline{\mathcal{L}}$. In particular, the natural gradient at $\lambda$ belongs to the tangent space to $\overline{\mathcal{M}}$ at $\lambda$.

We now need to associate the Riemannian gradient $\nabla_{\lambda}^{\text{nat}}\overline{\mathcal{L}}(\lambda)$ to the Riemannian gradient of the lower bound $\mathcal{L}(\lambda)$ defined in $\mathcal{M}$; the latter is what we need for using Riemannian SGD to optimize $\mathcal{L}(\lambda)$. This is done in the two cases: $\mathcal{M}$ is a submanifold (Section 3.1) and $\mathcal{M}$ is a quotient manifold (Section 3.2).

3.1 Riemannian submanifolds

Suppose that $\mathcal{M}$ is a submanifold of $\overline{\mathcal{M}}$. In order to define the Riemannian gradient of the lower bound $\mathcal{L}$ defined on the manifold $\mathcal{M}$, we need to equip $\mathcal{M}$ with a Riemannian metric. In most cases, this metric is inherited from that of $\overline{\mathcal{M}}$ in a natural way. Since $T_{\lambda}\mathcal{M}$ is a subspace of $T_{\lambda}\overline{\mathcal{M}}$, the Riemannian metric of $\zeta_{\lambda},\xi_{\lambda}\in T_{\lambda}\mathcal{M}$ can be defined as

[TABLE]

with $\zeta_{\lambda},\xi_{\lambda}$ viewed as vectors in $T_{\lambda}\overline{\mathcal{M}}$. With this metric, we can define the orthogonal complement $(T_{\lambda}\mathcal{M})^{\bot}$ of $T_{\lambda}\mathcal{M}$ in $T_{\lambda}\overline{\mathcal{M}}$, i.e. $T_{\lambda}\overline{\mathcal{M}}=T_{\lambda}\mathcal{M}\oplus(T_{\lambda}\mathcal{M})^{\bot}$. Write

[TABLE]

where $\text{Proj}_{\lambda}$ and $\text{Proj}_{\lambda}^{\bot}$ denote the projections on $T_{\lambda}\mathcal{M}$ and $(T_{\lambda}\mathcal{M})^{\bot}$, respectively. Recall that $\text{grad}\overline{\mathcal{L}}(\lambda)=\nabla_{\lambda}^{\text{nat}}\overline{\mathcal{L}}(\lambda)=I_{F}^{-1}(\lambda)\nabla_{\lambda}{\mathcal{L}}(\lambda)$, $\lambda\in\mathcal{M}$. Then, the Riemannian gradient of $\mathcal{L}$ is the projection of $\text{grad}\overline{\mathcal{L}}$ on $T_{\lambda}\mathcal{M}$

[TABLE]

This is because

[TABLE]

In some cases, such as Gaussian VB in Section 5.1, $T_{\lambda}\mathcal{M}\cong T_{\lambda}\overline{\mathcal{M}}$, then $\text{grad}\mathcal{L}(\lambda)=\text{grad}\overline{\mathcal{L}}(\lambda)$, i.e. the natural gradient is the Riemannian gradient of $\mathcal{L}$. In some other cases, however, using the inherited metric might lead to a projection $\text{Proj}_{\lambda}$ that is cumbersome to compute. In such cases, one needs to use an alternative Riemannian metric on $\mathcal{M}$ such that the projection $\text{Proj}_{\lambda}$ is easy to compute. Below we give an example in the case of Stiefel manifold, which is a popular manifold in many applications.

Stiefel manifolds. Suppose that $\mathcal{M}$ is a Stiefel manifold $\mathcal{M}=\mathcal{S}(p,n)$ defined as

[TABLE]

where $\text{Mat}(n,p)$ is the set of real matrices of size $n\times p$. We can think of $\mathcal{M}$ as embedded in $\overline{\mathcal{M}}=\mathbb{R}^{d}$, $d=n\cdot p$, equipped with the Fisher-Rao metric

[TABLE]

where $\text{vec}(\cdot)$ denotes the vectorization operator, $I_{F}(W)=\text{cov}_{q_{W}}(\nabla_{\text{vec}(W)}\log q_{W}(\theta))$ is the Fisher matrix defined in (2.2) with $q_{W}(\theta)$ the variational distribution. The natural gradient of function $\overline{\mathcal{L}}$ at $W\in\mathcal{M}$ is

[TABLE]

where $\text{vec}^{-1}$ is the inverse of vec, sending a $np$-vector to the corresponding matrix in $\text{Mat}(n,p)$. It is easy to see that the tangent space of $\mathcal{M}$ at $W$ is

[TABLE]

If we equip $\mathcal{M}$ with the Riemannian metric defined in (3.5), the projection on $T_{W}\mathcal{M}$ is cumbersome to compute. We therefore opt to use the usual Euclidean metric

[TABLE]

The following lemma gives an expression for the Riemannian gradient of $\mathcal{L}$ defined on the Stiefel manifold, and is useful for many applications involving Stiefel manifolds. Similar results to Lemma 3.2 can be found in the literature (see., e.g., (Edelman et al.,, 1998)), however, here we state and prove the results specifically for the case $G$ is the natural gradient.

Lemma 3.2**.**

Let $\mathcal{L}$ be a function on the Stiefel manifold $\mathcal{M}$ equipped with the usual Euclidean metric. The Riemannian gradient of $\mathcal{L}$ at $W$ is

[TABLE]

with $G=\text{grad}\overline{\mathcal{L}}(W)$ given in (3.6), and $\text{skew}(A):=(A-A^{\top})/2$.

3.2 Quotient manifolds

This section derives the Riemannian gradient of $\mathcal{L}$ when $\mathcal{M}$ is a quotient manifold induced from the ambient manifold $\overline{\mathcal{M}}$. Suppose that $\overline{\mathcal{M}}\subset\mathbb{R}^{d}$ is a Riemannian manifold with the Riemannian metric $<\cdot,\cdot>_{\overline{\mathcal{M}}}$. Suppose that there is an equivalence relation on $\overline{\mathcal{M}}$ defined as

[TABLE]

and thus $\mathcal{L}(\lambda)=\mathcal{L}(\lambda^{\prime})$. This is the case of Gaussian VB with the covariance matrix $\Sigma$ having a factor decomposition. Define the equivalence class

[TABLE]

i.e., the class of all parameterizations $\lambda$ that represent the same distribution. Let

[TABLE]

and define the canonical projection

[TABLE]

Then we can endow the quotient set $\mathcal{M}$ with the topology induced from $\overline{\mathcal{M}}$ by the projection $\pi$. This makes $\mathcal{M}$ become a smooth manifold called the quotient manifold, see Absil et al., (2009). If we define $L:\mathcal{\mathcal{M}}\to\mathbb{R},\;[\lambda]\mapsto L([\lambda])=\mathcal{L}(\lambda)$, i.e. $\mathcal{L}=L\circ\pi$, then $L$ is the function defined on $\mathcal{\mathcal{M}}$ that we want to optimise. For optimisation on $\mathcal{\mathcal{M}}$, one needs to be able to represent numerically tangent vectors at each $[\lambda]\in\mathcal{M}$. Geometrical objects in $\overline{\mathcal{M}}$, such as points $\lambda$ and tangent vectors, are vectors in the usual sense, so they can be numerically represented in computer for numerical computation. However, geometrical objects in the quotient manifold ${\mathcal{M}}$ are abstract, much of reseach in quotient manifolds has been focused on how to represent these geometrical objects numerically. The key tool is the concept of horizontal lift; see, e.g. Absil et al., (2009); Kobayashi and Nomizu, (1969).

By the level set theorem (Tu,, 2011, Chapter 2), $\pi^{-1}([\lambda])$ is an embedded submanifold in $\overline{\mathcal{M}}$, hence, it admits a tangent space

[TABLE]

called the vertical space, which is a linear subspace of $T_{\lambda}\overline{\mathcal{M}}$. Let $\mathcal{H}_{\lambda}$ be the orthogonal complement of $\mathcal{V}_{\lambda}$ in $T_{\lambda}\overline{\mathcal{M}}$, called the horizontal space, i.e. $T_{\lambda}\overline{\mathcal{M}}=\mathcal{H}_{\lambda}\oplus\mathcal{V}_{\lambda}.$ The orthogonality here is w.r.t. the metric defined on $\overline{\mathcal{M}}$. For each tangent vector $\xi_{[\lambda]}$ at $[\lambda]\in\mathcal{M}$, there exists an unique vector $\overline{\xi}_{\lambda}$ in the horizontal space $\mathcal{H}_{\lambda}$ such that (Kobayashi and Nomizu,, 1969, Prop. 1.2)

[TABLE]

$\overline{\xi}_{\lambda}$ is called the horizontal lift of $\xi_{[\lambda]}$. Then, as $\mathcal{L}=L\circ\pi$,

[TABLE]

which shows that the directional derivative of $L$ in the direction of $\xi_{[\lambda]}$ is characterised by the directional derivative of $\mathcal{L}$ in the direction of the horizontal lift $\overline{\xi}_{\lambda}$. Intuitively, for the optimization purposes, we can ignore the vertical space and just focus on the horizontal space, as the objective function $\mathcal{L}$ doesn’t change along the vertical space. It’s worth noting that the property (3.9) does not depend on any particular choice $\lambda$ in $[\lambda]$.

Let $\overline{\text{grad}\mathcal{L}(\lambda)}$ be the Riemannian gradient of $\mathcal{L}$ at $\lambda\in\overline{\mathcal{M}}$. We have that, for all $\eta_{\lambda}\in\mathcal{V}_{\lambda}$

[TABLE]

as $\mathcal{L}(\lambda)$ doesn’t change along the vertical space, which shows that $\overline{\text{grad}\mathcal{L}(\lambda)}\in\mathcal{H}_{\lambda}$. Let $\text{grad}L([\lambda])$ be the tangent vector to $\mathcal{M}$ at $[\lambda]$ that has $\overline{\text{grad}\mathcal{L}(\lambda)}$ as its horizontal lift. Then, by equipping $\mathcal{M}$ with the inner product inherited from $\overline{\mathcal{M}}$,

[TABLE]

we have that $\text{grad}L([\lambda])$ is the Riemannian gradient of $L$ at $[\lambda]\in\mathcal{M}$. We note that (3.10) does not depend on the choice of $\lambda\in[\lambda]$. So, with the inherited inner product from $\overline{\mathcal{M}}$, the usual Riemannian gradient of $\mathcal{L}$ on $\overline{\mathcal{M}}$ is the horizontal lift of the Riemannian gradient of $L$ on $\mathcal{M}$. This remarkable property of quotient manifolds makes it convenient for numerical optimisation problems.

Remark 3.1**.**

Technically, in order for the inherited Riemannian metric on $\mathcal{M}$ to be well-defined, it is often required in the literature that $<\overline{\xi}_{\lambda},\overline{\eta}_{\lambda}>_{\overline{\mathcal{M}}}$ does not depend on $\lambda\in[\lambda]$. This condition is typically not satisfied when $\overline{\mathcal{M}}$ is equipped with the Fisher-Rao metric as considered in this paper. However, as we show above, the Riemannian gradient of $L$ is still well-defined without this requirement, as (3.10) holds for any $\lambda\in[\lambda]$.

3.3 Retraction

After deriving the Riemannian gradient, which is the steepest ascent direction of the lower bound function $\mathcal{L}$ at the current point on the manifold, we need to derive the *exponential map111In general, an exponential map at $\lambda$ is defined locally near $\lambda$. In order to define exponential map on the entire tangent space, the manifold needs to be complete, see, e.g., Tu, (2011). *, denoted by $\text{Exp}_{\lambda}(\xi_{\lambda})$, that projects a point on the tangent space back to the manifold. Exponential map is a standard concept in differential geometry. Intuitively, exponential maps are mappings that, given a point $\lambda$ on a manifold and a tangent vector $\xi_{\lambda}$ at $\lambda$, generalize the concept “$\lambda+\xi_{\lambda}$” in Euclidean spaces. $\text{Exp}_{\lambda}(\xi_{\lambda})$ is a point on the manifold that can be reached by leaving from $\lambda$ and moving in the direction $\xi_{\lambda}$ while remaining on the manifold. We refer to Absil et al., (2009) for a precise definition and examples. One major drawback of exponential maps is that their calculation is often cumbersome in practice. Retraction, the first order approximation of the exponential map, is often used instead. A retraction $R_{\lambda}:T_{\lambda}\mathcal{M}\to\mathcal{M}$ at $\lambda\in\mathcal{M}$ has the important property that it preserves gradients, i.e. the curve $\gamma_{\xi_{\lambda}}:t\mapsto R_{\lambda}(t\xi_{\lambda})$ satisfies $\text{D}\gamma_{\xi_{\lambda}}(0)[\xi_{\lambda}]=\xi_{\lambda}$ for every $\xi_{\lambda}\in T_{\lambda}\mathcal{M}$. See (Absil et al.,, 2009, Chapter 4) for a formal definition of retraction and Manton, (2002) for an interpretation of retraction from a local optimization perspective. Also see Figure 1 (Left) for a visualization.

Closed-form formulae for retractions on common manifolds are available in the literature, see, for example, Absil et al., (2009). For instance, a popular retraction on the Stiefel manifold is

[TABLE]

Here $\text{qf}(A)=Q$, where $A=QR$ is the QR decomposition of $A\in\text{Mat}(n,p)$, $Q\in\mathcal{S}(p,n)$ and $R\in\text{Mat}(n,p)$ is upper triangular. See Sato and Aihara, (2019) for an efficient computation of this retraction based on the Cholesky QR factorization. For quotient manifolds, a popular retraction is

[TABLE]

$\overline{\xi}_{\lambda}$ is the horizontal lift of $\xi_{[\lambda]}$, and $\pi$ is the canonical projection in (3.8).

3.4 Momentum

The momentum method, which uses a moving average of the gradient vectors at the previous iterates to accelerate convergence and also help reduce noise in the estimated gradient, is widely used in Euclidean-based stochastic gradient optimization. Extending the momentum method to manifolds requires parallel translation, a tool in differential geometry for moving tangent vectors from one tangent space to another, while still preserving the length and angle (to some fixed direction) of the original tangent vectors. Similar to exponential map, a parallel translation is often approximated by a vector transport which is much easier to compute; see (Absil et al.,, 2009, Chapter 8) for a formal definition. See Figure 1 (Right) for a visualization. Let $\Gamma_{\lambda_{t}\to\lambda_{t+1}}(\xi_{\lambda_{t}})$ denote the vector transport of tangent vector $\xi_{\lambda_{t}}\in T_{\lambda_{t}}\mathcal{M}$ to tangent space $T_{\lambda_{t+1}}\mathcal{M}$. A simple vector transport is the projection of $\xi_{\lambda_{t}}$ on $T_{\lambda_{t+1}}\mathcal{M}$, i.e. $\Gamma_{\lambda_{t}\to\lambda_{t+1}}(\xi_{\lambda_{t}})=\text{Proj}_{T_{\lambda_{t+1}}\mathcal{M}}(\xi_{\lambda_{t}})$. Roy and Harandi, (2017) is the first to use the momentum method in Riemannian SGD, but they do not provide any convergence analysis.

3.5 Manifold VB algorithm

The pseudo-code in Algorithm 1 summarizes our VB algorithm on manifolds. We use the “hat” notation in $\widehat{\nabla_{\lambda}^{\text{nat}}\mathcal{L}(\lambda)}$ to emphasize that the natural gradient is obtained from a noisy and unbiased estimator $\widehat{\nabla_{\lambda}\mathcal{L}(\lambda)}$ of the Euclidean gradient as we often don’t have access to the exact ${\nabla_{\lambda}\mathcal{L}(\lambda)}$.

That is, starting from an initial $\lambda_{0}$, the initial momentum gradient $Y_{0}$ can be found by projecting the natural gradient estimate $\widehat{\nabla_{\lambda}^{\text{nat}}\mathcal{L}(\lambda_{0})}$ on $T_{\lambda_{0}}\mathcal{M}$. The relevant projections were described in Sections 3.1 and 3.2. At a step $t$, from $\lambda_{t}$, we move on the tangent space $T_{\lambda_{t}}\mathcal{M}$ along the direction $Y_{t}$ to find the next iterate $\lambda_{t+t}$ by retraction, $\lambda_{t+1}=R_{\lambda_{t}}(\varepsilon Y_{t})$. Then, we calculate the natural gradient estimate $\widehat{\nabla_{\lambda}^{\text{nat}}\mathcal{L}(\lambda_{t+1})}$ at $\lambda_{t+1}$, whose projection on $T_{\lambda_{t+1}}\mathcal{M}$ gives the Riemannian gradient $\text{grad}\mathcal{L}(\lambda_{t+1})$. Finally, the new momentum gradient $Y_{t+1}$ is calculated from the vector transport of $Y_{t}$ (to $T_{\lambda_{t+1}}\mathcal{M}$) and $\text{grad}\mathcal{L}(\lambda_{t+1})$.

4 Convergence analysis

To be consistent with the standard notation in the optimization literature, and with an abuse of notation, let us define the cost function as $\mathcal{L}(\cdot):=-\mathcal{L}(\cdot)$. That is, our optimization problem is

[TABLE]

In this section, for notational simplicity, we will denote by $\nabla\mathcal{L}(\lambda)$ the Riemannian gradient of the cost function $\mathcal{L}$, and by $\widehat{\nabla\mathcal{L}}(\lambda)$ its unbiased estimator. Let $\{\lambda_{t},t\leq 0\}$ be the iterates from Algorithm 1, and $\mathcal{F}_{t}$ be the $\sigma$-field generated by $\{\lambda_{s},s\leq t\}$. Because of the unbiasedness, we can write

[TABLE]

with $\{\Delta M_{t}\}$ a martingale difference w.r.t. $\{\mathcal{F}_{t}\}$, i.e. $\mathbb{E}(\Delta M_{t+1}|\mathcal{F}_{t})=0$. For the purpose of convergence analysis, we write our manifold VB algorithm as follows

[TABLE]

The minus sign in $R_{\lambda_{t}}(-Y_{t})$ results from the change in the notation $\mathcal{L}(\cdot):=-\mathcal{L}(\cdot)$ above, and with the suitable choice of $\zeta$ and $\gamma$ we can recover Algorithm 1.

We next need some definitions.

Definition 4.1**.**

(Huang et al., 2015a*, **, Section 3.3)**

A neighborhood $\mathcal{S}\subset\mathcal{M}$ of $x^{*}$ is said to be totally retractive if there is $\delta>0$ such that for any $y\in\mathcal{S}$, $R_{y}(B(0_{y},\delta))\supset\mathcal{S}$ and $R_{y}$ is a diffeomorphism on $B(0_{y},\delta)$, where $B(0_{y},\delta)$ is the ball of radius $\delta$ in $T_{y}\mathcal{M}$ centered at the origin $0_{y}$.*

Definition 4.2**.**

(Huang et al., 2015b*, *, Definition 3.1)** For a function $f:\mathcal{M}\mapsto\mathbb{R}$ on a Riemannian manifold $\mathcal{M}$ with retraction $R$, define $m_{\lambda,\eta}(t)=f(R_{\lambda}(t\eta))$ for $\lambda\in\mathcal{M}$, $\eta\in T_{\lambda}\mathcal{M}$. The function f is retraction-convex with respect to the retraction $R$ in a set $\mathcal{S}$ if for all $\lambda\in\mathcal{M}$, $\eta\in T_{\lambda}\mathcal{M},|\eta|=1$, $m_{\lambda,\eta}(t)$ is convex for all t which satisfy $R_{\lambda}(s\eta)\in\mathcal{S}$ for all $s\in[0,t]$. Moreover, $f$ is strongly retraction-convex in $\mathcal{S}$ if $m_{\lambda,\eta}(t)$ is strongly convex, that is $\tilde{a}_{0}\leq\dfrac{d^{2}m_{\lambda,\eta}}{dt^{2}}(t)\leq\tilde{a}_{1}$ for some positive constants $\tilde{a}_{1},\tilde{a}_{0}$.

Remark 4.1**.**

If $m_{\lambda,\eta}(t)$ is strongly convex with $\tilde{a}_{0}\leq\dfrac{d^{2}m_{\lambda,\eta}}{dt^{2}}(t)\leq\tilde{a}_{1}$ then

[TABLE]

where $\lambda=R_{\lambda^{*}}(t\eta).$ If we assume that $f$ is strongly retraction-convex in $\mathcal{S}$ and for any $\tilde{\lambda}\in\mathcal{S}$, there exists $\tilde{\eta}$ such that $R_{\lambda^{*}}(\tilde{\eta})=\tilde{\lambda}$ and the derivative $DR_{\lambda}(\eta)$ is bounded, then the chain rule implies $|\dfrac{dm_{\lambda^{*},\eta}}{dt}(t)|\leq c|\nabla f(\lambda)|$ with $\lambda=R_{\lambda^{*}}(t\eta).$ As a result,

[TABLE]

In Theorem 1, we show the convergence of (4.2) under suitable conditions imposed on the objective function $\mathcal{L}$. It is worth emphasizing that the convergence analysis in this section is done in a general setting for Riemannian SGD with momentum rather than restricting on the VB problem in the previous sections. It can therefore be applied to more general settings.

Theorem 1**.**

• Assume that
• •

There exists a totally retractive neighborhood of $\lambda^{*}$, $\mathcal{S}$, such that $\lambda_{t}\in\mathcal{S}$ for any $t\geq 0$.

• •

$\nabla\mathcal{L}$* and $\Delta M_{t}$ are bounded such that $|\nabla\mathcal{L}|+|\Delta M|\leq b_{\mathcal{L}}$ almost surely for some constant $b_{\mathcal{L}}>0$.*

• •

$\nabla\mathcal{L}(\lambda)$* is $\tilde{L}$-Lipschitz with respect to retraction $R$, that is, $|\nabla\mathcal{L}(R_{\lambda}(\eta))-\nabla\mathcal{L}(\lambda)|\leq\tilde{L}|\eta|$ for $\lambda\in\mathcal{S}$, $\eta\in T_{\lambda}\mathcal{M}$.*

Consider the sequence $(\lambda_{t})_{t\in\mathbb{N}}$ obtained from (4.2) using $\gamma=\frac{1}{\sqrt{T}}$. The following holds true:

[TABLE]

Moreover, when the objective function $\mathcal{L}$ is strongly retraction-convex, for $\epsilon\in(0,1)$, by choosing $\gamma=\frac{1}{T^{\epsilon}}$, there exists a constant $C_{\epsilon}$ such that

[TABLE]

The proof can be found in the Appendix. We note that the first assumption in Theorem 1 is standard, see, for example Huang et al., 2015a . The condition $|\nabla\mathcal{L}|+|\Delta M|\leq b_{\mathcal{L}}$ is essential to make sure that $Y_{t}$ stays bounded. That is, it does not diverge to infinity, which is the key property of the algorithm. While the Lipschitz property of $\nabla\mathcal{L}(\lambda)$ guarantee that we can expand $\mathcal{L}$ using Taylor expansion. The assumption related to the martingale difference $\Delta M_{t}$ is justified by the fact that the estimator $\widehat{\nabla\mathcal{L}}$ is unbiased.

5 Applications

5.1 Manifold Gaussian VB

Gaussian VB (GVB) uses a multivariate Gaussian distribution $N(\mu,\Sigma)$ for the VB approximation $q_{\lambda}$, $\lambda=(\mu,\text{vec}(\Sigma))$. GVB has been extensively used in the literature, often with some simplifications imposed on $\Sigma$; e.g., $\Sigma$ is a diagonal matrix ${\rm diag}(\sigma_{1}^{2},...,\sigma_{d}^{2})$ or has a factor structure $\Sigma=BB^{\top}+D^{2}$. One of the reasons of imposing these simplifications is to deal with the symmetric and positive definiteness constraints on $\Sigma$. We do not impose any simplifications on $\Sigma$, and deal with these constraints by considering the VB optimization problem on the manifold $\mathcal{M}$ of symmetric and positive definite matrices $\mathcal{M}=\{\Sigma\in\text{Mat}(d,d):\Sigma=\Sigma^{\top},\;\;\Sigma>0\}$. We can think of $\mathcal{M}$ as embedded in $\overline{\mathcal{M}}=\mathbb{R}^{d^{2}}$.

From (2.1), the gradient of lower bound is (Tran et al.,, 2021)

[TABLE]

Here, we use the so-called score-function VB as in Tran et al., (2017), which does not require the gradient of the log-likelihood. We follow Tran et al., (2017) and use a control variate for the gradient of lower bound

[TABLE]

where $c$ is a vector selected to minimize the variance of the gradient estimate

[TABLE]

with $d_{\lambda}$ the size of $\lambda$, which can be estimated by sampling from $q_{\lambda}$.

Mardia and Marshall, (1984) show that the Fisher information matrix for the multivariate Gaussian distribution $N(\mu,\Sigma)$ is

[TABLE]

where $I_{F}(\Sigma)$ is an $d^{2}\times d^{2}$ matrix with entries

[TABLE]

One can derive that $I_{F}(\Sigma)\approx\Sigma^{-1}\otimes\Sigma^{-1}$, with $\otimes$ the Kronecker product. Therefore

[TABLE]

which gives a convenient form for obtaining an approximate natural gradient. The natural gradient w.r.t. $\mu$ and $\Sigma$ is approximated as

[TABLE]

As $T_{\Sigma}\mathcal{M}=\text{Mat}(d,d)\cong T_{\Sigma}\overline{\mathcal{M}}\cong\mathbb{R}^{d^{2}}$, the projection in (3.3) is the identity, hence $\nabla_{\Sigma}^{\text{nat}}\mathcal{L}(\lambda)$ is the Riemannian gradient of lower bound $\mathcal{L}$ w.r.t. $\Sigma$. The Manifold GVB algorithm is outlined in Algorithm 2.

One of the most popular retractions used for the manifold $\mathcal{M}$ of symmetric and positive definite matrices is (see the Manopt toolbox of Boumal et al., (2014))

[TABLE]

and vector transport

[TABLE]

The Matlab code implementing the Manifold GVB algorithm is made available online at **https://github.com/VBayesLab. **

Numerical experiments

We apply the Manifold GVB algorithm to fitting a logistic regression model using the German Credit dataset. This dataset, available on the UCI Machine Learning Repository https://archive.ics.uci.edu/ml/index.php, consists of observations on 1000 customers, each was already rated as being “good credit” (700 cases) or “bad credit” (300 cases). The covariates include credit history, education, employment status, etc. and lead to totally 25 predictors after using dummy variables to represent the categorical covariates.

A naive GVB implementation is to only update $\Sigma$ when its updated value satisfies the symmetric and positive definiteness constraint. This naive implementation didn’t work at all in this example. To see the usefulness of incorporating the natural gradient into the Manifold GVB, we compare Algorithm 2 with a version without using the natural gradient. As shown in Figure 2, using the natural gradient leads to a much faster and more stable convergence. Also, the Manifold GVB without the natural gradient requires a large number of samples used in estimating the gradient (5.1) (we used $S=10,000$), compared to $S=100$ for the Manifold GVB with the natural gradient. The CPU running time for the Manifold GVB algorithms with and without the natural gradient is 43 and 280 seconds, respectively.

Tran et al., (2019) develop a Gaussian VB algorithm where $\Sigma$ is factorized as $\Sigma=BB^{\top}+D^{2}$ with $B$ a vector, and term their algorithm NAGVAC. Figure 3 plots the lower bound estimates of the Manifold GVB and NAGVAC. Manifold GVB stopped after 921 iterations and NAGVAC stopped after 1280 iterations, and their CPU running times are 19 seconds and 12 seconds, respectively. As shown, the Manifold GVB algorithm converges quicker than NAGVAC and obtains a larger lower bound. We note, however, that NAGVAC is less computational demanding than Manifold GVB in high-dimensional settings such as deep neural networks.

To assess the training stability of the Manifold GVB and NAGVAC algorithms, we use the same initialization $\mu_{0}$ and $\Sigma_{0}$ for both algorithms and run each for 20 different replications. The standard deviations of the estimates of $\mu$ (across the different runs, then averaged over the 25 coordinates) for NAGVAC and Manifold GVB are 0.03 and 0.01, respectively. This demonstrates that the Manifold GVB algorithm is more stable than NAGVAC. To assess their sensitivity to the initialization, in each algorithm, we now use a random initialization but fix the random seed in the updating stage. The standard deviations of the estimates of $\mu$ (across the different runs, then averaged over the 25 coordinates) for NAGVAC and Manifold GVB are 0.0074 and 0.0009, respectively. This demonstrates that the Manifold GVB algorithm is less sensitive to the initialization than NAGVAC.

Application to financial time series data

This section applies the MGVB method to analyze a financial stock return data set, and compares MGVB to Sequential Monte Carlo (SMC). We consider the GARCH model of Bollerslev, (1986) for modelling the underlying volatility dynamics in the Standard & Poor 500 stock indices observed from 4 Jan 1988 to 26 Feb 2007 (1000 observations). Let $\{y_{t},t=1,...,n\}$ be the stock returns. We consider the following GARCH model

[TABLE]

The parameters $\theta$ are $w>0$, $\alpha>0$ and $\beta>0$, with the constraint $\alpha+\beta<1$ to ensure the stationarity. To impose this constraint, we parameterize $\alpha$ and $\beta$ as $\alpha=\psi_{1}(1-\psi_{2})$ and $\beta=\psi_{1}\psi_{2}$ with $0<\psi_{1},\psi_{2}<1$. We use an inverse Gamma prior IG$(1,1)$ for $w$ and an uniform prior $U(0,1)$ for $\psi_{1}$ and $\psi_{2}$. Finally, we use the following transformation

[TABLE]

and work with the unconstrained parameters $\widetilde{\theta}=(\theta_{w},\theta_{\psi_{1}},\theta_{\psi_{1}})$, but we will report the results in terms of the original parameters $\theta=(w,\alpha,\beta)$.

We compare the MGVB method to SMC. In this application with only three unknown parameters, SMC is applicable and can be considered as the “gold standard” as it produces an asymptotically exact approximation of the posterior $p(\widetilde{\theta}|y)$. We implement the likelihood annealing SMC method of Tran et al., (2014), which is a robust SMC sampling technique that first draws samples from an easily-generated distribution and then moves these samples via annealing distributions towards the posterior distribution through weighting, resampling and Markov moving. We run SMC with 10,000 particles and the annealing distributions are adaptively designed such that the effective sample size is always at least 80%. Figure 4 plots the posterior estimates for the model parameters $\theta$ by SMC and MGVB, which shows that the MGVB estimates are almost identical to that of SMC. The CPU running time for MGVB and SMC is 10.4 and 380.5 seconds, respectively.

5.2 Manifold Wishart Variational Bayes

Suppose that we are interested in approximating the posterior distribution of a covariance matrix $V$ of size $d$ by an inverse-Wishart distribution $q_{\lambda}(V)=\text{inverse-Wishart}(\nu_{q},\Sigma_{q})$. The density $q_{\lambda}$ is

[TABLE]

with

[TABLE]

the multivariate gamma function. The gradient of $\log q_{\lambda}(V)$ w.r.t. $\nu_{q}$ and $\Sigma_{q}$ are

[TABLE]

and

[TABLE]

where $\psi_{d}(\nu)\coloneqq\partial\log\Gamma_{d}(\nu)/\partial\nu$. From this, it is straightforward to estimate the gradients of the lower bound $\nabla_{\nu_{q}}\mathcal{L}(\lambda)$ and $\nabla_{\Sigma_{q}}\mathcal{L}(\lambda)$ using the score-function method. It is often much more efficient to use

[TABLE]

instead of $\nabla_{\Sigma_{q}}\mathcal{L}(\lambda)$, which is analogous to the use of the natural gradient in (5.4). The natural-gradient-like expression in (5.7) is often found very useful in practice; without it, it is difficult for the manifold Wishart VB algorithm that only uses the Euclidean $\nabla_{\Sigma_{q}}\mathcal{L}(\lambda)$ to converge. As $\Sigma_{q}$ lies on the manifold of symmetric and positive definite matrices, the retraction and vector transport in (5.5) and (5.6) can be used. To update $\nu_{q}$, it is recommended to use some adaptive learning method such as ADAM or AdaDelta, or an approximate natural gradient. If we set the correlation between $\nabla_{\nu_{q}}\log q_{\lambda}(V)$ and $\nabla_{\Sigma_{q}}\log q_{\lambda}(V)$ to be zero, then the natural gradient of the lower bound with respect to $\nu_{q}$ is approximated by

[TABLE]

where $\psi_{d}^{\prime}(\nu)\coloneqq\partial\psi_{d}(\nu)/\partial\nu$.

A numerical example

Data of size $n$ is generated from $\mathcal{N}_{d}(0,V)$, where the elements of the true covariance matrix $V$ are $v_{ij}=(-0.5)^{|i-j|}$. An inverse-Wishart($\nu_{0},S_{0}$) prior is used for $V$ with $\nu_{0}=d$ and $S_{0}=0.01I_{d}$. It is easy to see that the posterior distribution of $V$ is inverse-Wishart with the degree of freedom $\nu=n+d$ and scale matrix $S=S_{0}+\sum_{i}y_{i}y_{i}^{\top}$. We run the manifold Wishart VB algorithm with the inital value $\nu_{q}=n$ and $\Sigma_{q}=n\times S_{n}$, where $S_{n}$ is the sample covariance matrix of the $y_{i}$.

In the first simulation, we consider $n=50$ and $d=5$. The estimation result is summarized in Table 1. We consider the second simulation with $n=500$ and $d=50$. The lower bound and a plot of the true posterior means of the $v_{ij}$ v.s. their VB estimates are shown in Figure 5. These results show that the manifold Wishart VB algorithm appears to work effectively and efficiently in this example.

6 Conclusion

We proposed a manifold-based Variational Bayes algorithm that takes into account both information geometry and geometric structure of the constraint parameter space. The algorithm is provably convergent for the case when the objective function either non-convex or strongly retraction-convex. Our numerical experiments demonstrate that the proposed algorithm converges quickly, is stable and compares favourably to the existing VB methods. An interesting research direction is to develop a Manifold VB method for the Sylvester normalizing flows of Berg et al., (2018). We leave it as an interesting project for future research.

Acknowledgment

Tran would like to thank Prof. Pierre-Antoine Absil for useful discussions, and Nghia Nguyen for his help with the SMC implementation.

Appendix

Proof of Lemma 3.1.

By definition, the Riemannian gradient of $\overline{\mathcal{L}}$ at $\lambda\in\overline{\mathcal{M}}$, denoted by $\text{grad}\overline{\mathcal{L}}(\lambda)$, is the unique tangent vector in $T_{\lambda}\overline{\mathcal{M}}$ such that

[TABLE]

That is

[TABLE]

The natural gradient $\nabla_{\lambda}^{\text{nat}}\overline{\mathcal{L}}(\lambda)$ in (3.2) satisfies (6.1). Indeed, due to the symmetry of $I_{F}(\lambda)$, we have

[TABLE]

If $\zeta_{\lambda}\in T_{\lambda}\overline{\mathcal{M}}$ also satisfies (6.1), then

[TABLE]

which implies that $\zeta_{\lambda}=\text{grad}\overline{\mathcal{L}}(\lambda)$. ∎

Proof of Lemma 3.2.

Much of our proof is taken from Tagare, (2011). The idea of the proof is that we want to find a vector in $T_{W}\mathcal{M}$ that represents the action of the differential $\text{D}\mathcal{L}(W)$ via the natural gradient $G$. Let $W_{\bot}\in\text{Mat}(n,n-p)$ such that its columns together with the columns of $W$ form an orthonormal basis for $\mathbb{R}^{n}$. As $[W,W_{\bot}]$ is an orthogonal matrix, for any $U\in\text{Mat}(n,p)$, there exists a $C\in\text{Mat}(n,p)$ such that $U=[W,W_{\bot}]C$. Write $C=\begin{pmatrix}U_{W}\\ U_{W_{\bot}}\end{pmatrix}$ with $U_{W}$ the $p\times p$-matrix formed by the first $p$ rows of $C$, $U_{W_{\bot}}$ the $(n-p)\times p$ matrix formed by the last $n-p$ rows of $C$. That is, any matrix $U\in\text{Mat}(n,p)$ can be written as

[TABLE]

If $Z=WZ_{W}+W_{\bot}Z_{W_{\bot}}\in T_{W}\mathcal{M}$, then $Z^{\top}W+W^{\top}Z=0_{p\times p}$ implies that $Z_{W}+Z_{W}^{\top}=0$. So $T_{W}\mathcal{M}$ is a subset of the set

[TABLE]

It’s easy to check that this set is also a subset of $T_{W}\mathcal{M}$. We arrive at an alternative representation of the tangent space $T_{W}\mathcal{M}$

[TABLE]

We want to find a vector in $T_{W}\mathcal{M}$ that represents the action of the differential $\text{D}\mathcal{L}(W)$ on $T_{W}\mathcal{M}$, with $\text{D}\mathcal{L}(W)$ characterized by $G$; i.e. find $U=WU_{W}+W_{\bot}U_{W_{\bot}}\in T_{W}\mathcal{M}$, $U_{W}=-U_{W}^{\top}$ such that

[TABLE]

As the gradient $G=\text{grad}\overline{\mathcal{L}}(W)\in\text{Mat}(n,p)$, it can be written as $G=WG_{W}+W_{\bot}G_{W_{\bot}}$. Based on the natural gradient $G$,

[TABLE]

where we have used the fact that $G_{W}=\text{skew}(G_{W})+\text{sym}(G_{W})$ with $\text{sym}(G_{W})=(G_{W}+G_{W}^{\top})/2$, and that $\text{tr}(\text{sym}(G_{W})^{\top}Z_{W})=0$. We have that

[TABLE]

Comparing (Proof of Lemma 3.2.) and (6.5) gives

[TABLE]

As $G_{W}=W^{\top}G$ and $W_{\bot}G_{W_{\bot}}=G-WG_{W}=(I_{n}-WW^{\top})G$,

[TABLE]

From (6.3), $U$ is the Riemannian gradient of $\mathcal{L}$. This completes the proof. ∎

Proof of Theorem 1.

(i) First, we consider the case $\mathcal{L}$ is not assumed to be convex. Consider the update

[TABLE]

where $\zeta,\gamma\in(0,1)$, $\Delta M_{t+1}$ is a martingale difference. Since $|\nabla\mathcal{L}|+|\Delta M|\leq b_{\mathcal{L}}$ a.s, we have $|Y_{t}|\leq\frac{2\gamma b_{\mathcal{L}}}{1-\zeta}$. This can be proved as follows. First for the rest of this proof, we use $O(a)$ to denote vector/scalar with $|O(a)|\leq a$. Let $P_{x\to y}$ be a parallel translation. From Lemma 6 Huang et al., 2015a , there exists a constant $a_{1}>0$ such that

[TABLE]

which, together with $|P_{\lambda_{t}\to\lambda_{t+1}}(Y_{t})|=|Y_{t}|$ implies

[TABLE]

If we take norm of the second equation in (6.6), and using (6.8) we have

[TABLE]

With $\gamma<\frac{(1-\zeta)^{2}}{4a_{1}b_{\mathcal{L}}}$ and $|Y_{0}|\leq\frac{2\gamma b_{\mathcal{L}}}{1-\zeta}$ if $\zeta,\gamma$ are sufficiently small, using induction, we can show that $|Y_{t}|\leq\frac{2\gamma b_{\mathcal{L}}}{1-\zeta}$. Indeed assuming $|Y_{t}|\leq\frac{2\gamma b_{\mathcal{L}}}{1-\zeta}$, using (Proof of Theorem 1.), we have

[TABLE]

Next, applying the Taylor expansion to the function $t\mapsto\mathcal{L}(R_{\lambda}(t\eta))$ and using the Lipschitz continuity of $\nabla\mathcal{L}$ as well as the smoothness of $R_{\lambda}(\cdot)$ in $\mathcal{S}$, there exists a constant $L$ such that

[TABLE]

From (6.6), squaring both sides of the second equation we obtain

[TABLE]

where

[TABLE]

From (6.8), with $\gamma$ small enough so that $\frac{a_{1}\gamma b_{\mathcal{L}}}{1-\zeta}\leq 1$, we have

[TABLE]

Consider the second term without the constant in (6.11), we have

[TABLE]

Recall that $P_{x\to y}$ is an isometry, see, for example Absil et al., (2009). As a result we have

[TABLE]

In view of the fundamental theorem of calculus (see (Huang et al., 2015a, , Lemma 8)), we have $P_{\lambda_{t}\to\lambda_{t+1}}^{-1}\nabla\mathcal{L}(\lambda_{t+1})=\nabla\mathcal{L}(\lambda_{t})+O(a_{2})|Y_{t}|$ for some constant $a_{2}>0$. Then, for $\gamma\leq\frac{1-\zeta}{a_{2}}$, we have

[TABLE]

and

[TABLE]

For (II), since $\nabla\mathcal{L}$ is a.s bounded by $b_{\mathcal{L}}$, we derive from (6.8) that

[TABLE]

As a result of (I) and (II), we have

[TABLE]

Plugging (6.15) into (6.11),

[TABLE]

where $a_{3}=\frac{6a_{1}b_{\mathcal{L}}}{1-\zeta}+2(a_{1}b_{\mathcal{L}}+a_{2})\zeta$. Next, using the second equation in (6.6) and (6.15), we have

[TABLE]

Multiplying (6.16), (6.17) with $A$ (to be chosen later) and $-\frac{1}{1-\zeta}$ respectively and then adding to (6.10) we have

[TABLE]

Note that,

[TABLE]

(In the second inequality above, we have used the fact that $|<p,q>|\leq\frac{p^{2}}{2\epsilon}+\frac{\epsilon q^{2}}{2}$ with $\epsilon=(1-\zeta)$). Plugging back to (6.18), we obtain

[TABLE]

Now select $A=\frac{4}{(1-\zeta^{2})}\left(\frac{L}{2}{-\zeta}\frac{a_{1}b_{\mathcal{L}}+a_{2}}{1-\zeta}\right)$ and let $\gamma$ satisfy

[TABLE]

taking the expectation of both sides of (6.19), we have

[TABLE]

Take the sum, we have

[TABLE]

Then

[TABLE]

Hence

[TABLE]

If we choose $\gamma^{2}=\frac{1}{T}$, then

[TABLE]

(ii) Now we assume that $\mathcal{L}$ is strongly retraction convex with

[TABLE]

for some $\tilde{\mu}>0$ (see Remark 4.1). Let $\rho:=1-\frac{\gamma}{4\tilde{\mu}(1-\zeta)}>\frac{1+\zeta}{2}$ when $\gamma$ is small. Rewrite (6.17)

[TABLE]

Multiplying (6.16), (6.22) with $A$ (to be chosen later) and $-\frac{1}{\rho-\zeta}$ respectively and then adding to (6.10) we have

[TABLE]

We have

[TABLE]

Plugging this into (6.23), we have

[TABLE]

We deduce from (6.21) that

[TABLE]

On the other hand due to (6.12) we have

[TABLE]

Adding (6.24) and (Proof of Theorem 1.) and using (Proof of Theorem 1.) we have for

[TABLE]

that

[TABLE]

Since $\mathcal{L}$ is strongly retraction-convex, there exists $\overline{\mu}>0$ such that $|\nabla\mathcal{L}(\lambda)|^{2}\leq\tilde{L}^{2}|R^{-1}_{\lambda^{*}}(\lambda)|^{2}\leq\tilde{L}^{2}\overline{\mu}(\mathcal{L}(\lambda_{T})-\mathcal{L}(\lambda_{*}))$. As a result, if $A>\frac{2\tilde{L}^{2}\overline{\mu}}{(1-\zeta)^{2}}>\frac{\tilde{L}^{2}\overline{\mu}}{2(\rho-\zeta)^{2}}$ then we have from Cauchy’s inequality that

[TABLE]

Let $A$ and $\gamma$ be such that

[TABLE]

and

[TABLE]

we have

[TABLE]

As a result, recall the definition of $\rho$, we have

[TABLE]

where $C_{0}=\mathbb{E}V_{0}$.

With $\gamma=\frac{1}{T^{\varepsilon}}$, we have

[TABLE]

When $T$ is large then

[TABLE]

Thus for $C_{\varepsilon}=\max\{C_{0},C\}$,

[TABLE]

we have

[TABLE]

∎