Manifold Optimization Assisted Gaussian Variational Approximation

TL;DR

This paper introduces manifold optimization techniques to improve Gaussian variational approximation in Bayesian inference, addressing covariance structure constraints and enhancing optimization stability and efficiency.

Contribution

It proposes novel Riemannian stochastic gradient descent schemes on Stiefel and Grassmann manifolds for Gaussian variational approximation, overcoming covariance constraint issues.

Findings

Achieves competitive accuracy in high-dimensional Bayesian tasks.

Demonstrates faster convergence with the proposed methods.

Maintains stability and robustness during optimization.

Abstract

Gaussian variational approximation is a popular methodology to approximate posterior distributions in Bayesian inference especially in high dimensional and large data settings. To control the computational cost while being able to capture the correlations among the variables, the low rank plus diagonal structure was introduced in the previous literature for the Gaussian covariance matrix. For a specific Bayesian learning task, the uniqueness of the solution is usually ensured by imposing stringent constraints on the parameterized covariance matrix, which could break down during the optimization process. In this paper, we consider two special covariance structures by applying the Stiefel manifold and Grassmann manifold constraints, to address the optimization difficulty in such factorization architectures. To speed up the updating process with minimum hyperparameter-tuning efforts, we design two new schemes of Riemannian stochastic gradient descent methods and compare them with other existing methods of optimizing on manifolds. In addition to fixing the identification issue, results from both simulation and empirical experiments prove the ability of the proposed methods of obtaining competitive accuracy and comparable converge speed in both high-dimensional and large-scale learning tasks.

Figures

Tables

Table 1. Table 1: Riemannian operations for the Stiefel manifold structure.

Manifolds	Stiefel $𝒮(m,p)$
Tangent Spaces	$T_B𝒮(m,p)=\{U\in {ℝ}^{m\times p}:\text{sym}\left(B^{\top }U\right)=0\}$
Projection ${\pi }_B\left(Z\right)$	$Z-B\text{sym}\left(B^{\top }Z\right)$
Retraction $R_B\left(U\right)$	$\left(B+U\right){\left(I+U^{\top }U\right)}^{-\frac{1}{2}}$
Parrallel Transport ${\Gamma }_{B_1\to B_2}\left(U\right)$	${\pi }_{B_2}\left(U\right)$

Table 2. Table 2: Riemannian operations for the Grassmann manifold structure.

Manifolds	Grassmann $𝒢(m,p)$
Tangent Spaces	$T_B𝒢(m,p)=\{U\in {ℝ}^{m\times p}:,B^{\top }U=0\}$
Projection ${\pi }_B\left(Z\right)$	$\left(I-B{B}^{\top }\right)Z$
Retraction $R_B\left(U\right)$	$\text{polar}\left(B+U\right)$
Parrallel Transport ${\Gamma }_{B_1\to B_2}\left(U\right)$	${\pi }_{B_2}\left(U\right)$

Table 3. Table 3: The average training error and test error rates for the ionosphere data with five-fold cross-validation and factor p = 4 𝑝 4 p=4 .

Method	Training Error(%)	Test Error(%)	Time(sec)
VAFC ($p=0$)	0.64	7.65	12.51
VAFC	0.36	8.22	23.04
SLANG	7.55	14.20	190.93
S-RGD-Basic	0.36	7.65	25.02
S-cRGD-M	0.43	7.65	24.71
S-RGD-RMSProp	0.43	7.09	25.29
S-RGD-AdaDelta	0.43	7.37	27.27
G-RGD-Basic	1.00	7.38	24.78
G-cRGD-M	0.71	7.38	24.75
G-RGD-RMSProp	0.71	7.09	24.43
G-RGD-AdaDelta	0.64	7.09	25.59

Table 4. Table 4: Average training error, test error and running time for the Leukemia Cancer data with p = 4 𝑝 4 p=4 and iteration = 5000 absent 5000 =5000 .

Method	Training Error	Test Error	Time(sec)
VAFC ($p=0$)	0/38	7/34	64.08
VAFC	0/38	3/34	8021.71
SLANG	-	-	-
S-RGD-Basic	0/38	2/34	10013.90
S-cRGD-M	0/38	2/34	10031.65
S-RGD-RMSProp	0/38	2/34	10011.27
S-RGD-AdaDelta	0/38	2/34	10005.28
G-RGD-Basic	0/38	2/34	6482.94
G-cRGD-M	0/38	2/34	6438.22
G-RGD-RMSProp	0/38	3/34	6590.98
G-RGD-AdaDelta	0/38	2/34	6202.85

Equations

$$\text{KL}(q_{\lambda }(\mathbf{\theta })\Vert p(\mathbf{\theta }|\mathbf{x}))$$

$$\mathcal{L}(\lambda )=$$

$${\lambda }^{(t+1)}={\lambda }^{(t)}+{\alpha }_t\widehat{{\nabla }_{\lambda }\mathcal{L}}$$

$$\mathbf{\Sigma }={\mathbf{BB}}^{\top }+{\mathbf{D}}_2^2,$$

$$\mathcal{L}(\lambda )={\mathbb{E}}_q[\log h(\mathbf{\theta })]+\frac{1}{2}\log |{\mathbf{BB}}^{\top }+{\mathbf{D}}_2^2|+\frac{m}{2}(\log (2\pi )+1).$$

$$\mathcal{L}(\lambda )={\mathbb{E}}_q[\log h(\mathbf{\theta })]+\frac{1}{2}\log |{\mathbf{BB}}^{\top }+{\mathbf{D}}_2^2|.$$

$$\mathbf{\theta }=\mathbf{\mu }+\mathbf{Bz}+{\mathbf{d}}_2\circ \mathbf{ϵ}=\mathbf{\mu }+\mathbf{Bz}+{\mathbf{D}}_2\mathbf{ϵ},$$

$$\mathcal{L}(\lambda )={\mathbb{E}}_f[\log h(\mathbf{\mu }+\mathbf{Bz}+{\mathbf{d}}_2\circ \mathbf{ϵ})]+\frac{1}{2}\log |{\mathbf{BB}}^{\top }+{\mathbf{D}}_2^2|.$$

$$|\mathbf{BQ}(\mathbf{BQ}{)}^{\top }+{\mathbf{D}}_2^2|$$

$${\mathbb{E}}_f[\log h(\mathbf{\mu }+\mathbf{B}(\mathbf{Qz})+{\mathbf{d}}_2\circ \mathbf{ϵ})]$$

$$\mathcal{L}(\lambda (\mathbf{\mu },\mathbf{B},{\mathbf{D}}_2))=\mathcal{L}(\lambda (\mathbf{\mu },\mathbf{BQ},{\mathbf{D}}_2)).$$

$$\mathcal{S}(p,m)=\{\mathbf{B}:{\mathbf{B}}^{\top }\mathbf{B}={\mathbb{I}}_p\}.$$

$$\mathbf{\Sigma }={\mathbf{B}{\mathbf{D}}_1^2\mathbf{B}}^{\top }+{\mathbf{D}}_2^2,$$

$$\mathbf{\theta }=\mathbf{\mu }+\mathbf{B}{\mathbf{D}}_1\mathbf{z}+{\mathbf{d}}_2\circ \mathbf{ϵ}$$

$$\mathcal{L}(\lambda )=$$

$$[\mathbf{B}]=\{\mathbf{BQ}:\mathbf{B}\in \mathcal{S}(p,m)\text{ and }\mathbf{Q}\in \mathcal{O}(p)\},$$

$$\mathcal{L}(\lambda )=$$

$$x_{k+1}=x_k+{\alpha }_k{\xi }_k,$$

$$x_{k+1}={\mathrm{R}}_{x_k}({\alpha }_k{\xi }_k),$$

$${\xi }_k=-\text{grad}\mathcal{L}(x_k)+{\alpha }_k{\Gamma }_{x_{k-1}\to x_k}({\xi }_{k-1}),$$

$${\nabla }_{\mathbf{\mu }}\mathcal{L}={\mathbb{E}}_f[{\nabla }_{\mathbf{\theta }}\log h(\mathbf{\mu }+\mathbf{B}{\mathbf{D}}_1\mathbf{z}+{\mathbf{d}}_2\circ \mathbf{ϵ})];$$

$${\nabla }_{\mathbf{B}}\mathcal{L}={\mathbb{E}}_f[{\nabla }_{\mathbf{\theta }}\log h(\mathbf{\mu }+\mathbf{B}{\mathbf{D}}_1\mathbf{z}+{\mathbf{d}}_2\circ \mathbf{ϵ})(\mathbf{z}\circ {\mathbf{d}}_1{)}^{\top }]+({\mathbf{B}{\mathbf{D}}_1^2\mathbf{B}}^{\top }+{\mathbf{D}}_2^2{)}^{-1}{\mathbf{BD}}_1^2;$$

$${\nabla }_{{\mathbf{d}}_1}\mathcal{L}={\mathbb{E}}_f[(B^{\top }{\nabla }_{\mathbf{\theta }}\log h(\mathbf{\mu }+\mathbf{B}{\mathbf{D}}_1\mathbf{z}+{\mathbf{d}}_2\circ \mathbf{ϵ}))\circ \mathbf{z}]+\text{diag}(B^{\top }({\mathbf{B}{\mathbf{D}}_1^2\mathbf{B}}^{\top }+{\mathbf{D}}_2^2{)}^{-1}\mathbf{B})\circ {\mathbf{d}}_1;$$

$${\nabla }_{{\mathbf{d}}_2}\mathcal{L}={\mathbb{E}}_f[\text{diag}({\nabla }_{\mathbf{\theta }}\log h(\mathbf{\mu }+\mathbf{B}{\mathbf{D}}_1\mathbf{z}+{\mathbf{d}}_2\circ \mathbf{ϵ}){\mathbf{ϵ}}^{\top })]+\text{diag}(({\mathbf{B}{\mathbf{D}}_1^2\mathbf{B}}^{\top }+{\mathbf{D}}_2^2{)}^{-1})\circ {\mathbf{d}}_2,$$

$$({\mathbf{B}{\mathbf{D}}_1^2\mathbf{B}}^{\top }+{\mathbf{D}}_2^2{)}^{-1}={\mathbf{D}}_2^{-2}-{\mathbf{D}}_2^{-2}{\mathbf{BD}}_1({\mathbb{I}}_p+{{\mathbf{D}}_1{\mathbf{B}}^{\top }\mathbf{D}}_2^{-2}{\mathbf{BD}}_1{)}^{-1}{{\mathbf{D}}_1{\mathbf{B}}^{\top }\mathbf{D}}_2^{-2}.$$

$${\nabla }_{\mathbf{\mu }}\mathcal{L}$$

$${\nabla }_{\mathbf{B}}\mathcal{L}$$

$${\nabla }_{{\mathbf{d}}_2}\mathcal{L}$$

$$\text{grad}\mathcal{L}(\mathbf{B})={\pi }_{\mathbf{B}}({\nabla }_{\mathbf{B}}\mathcal{L}(\mathbf{B})).$$

$${\mathbf{B}}^{(t+1)}={\mathrm{R}}_{{\mathbf{B}}^{(t)}}\left(-\alpha {\text{grad}\mathcal{L}(\mathbf{B})|}_{{\mathbf{B}}^{(t)}}\right),$$

$${\mathbf{m}}^{(t+1)}=$$

$${\mathbf{B}}^{(t+1)}=$$

$${\mathbf{B}}^{(t+1)}={\mathrm{R}}_{{\mathbf{B}}^{(t)}}\left(-\alpha {\pi }_{{\mathbf{B}}^{(t)}}\left(\frac{\nabla \mathcal{L}({\mathbf{B}}^{(t)})}{\text{sgn}(E(g_{\mathbf{B}}^2{)}^{(t+1)})\odot \sqrt{|E(g_{\mathbf{B}}^2{)}^{(t+1)}|}+ϵ}\right)\right)$$

$$E(g_{\mathbf{B}}^2{)}^{(t+1)}=$$

Full text

Keywords: Bayesian Variational Auto-Encoder; Riemannian Manifolds; Stiefel Manifolds; Grassmann Manifolds; Riemannian Stochastic Gradient Method

1 Introduction

Variational inference is a well-established approach to approximate intractable posterior distributions in Bayesian statistics (Fox and Roberts,, 2012; Blei et al.,, 2017; Zhang et al., 2018a, ). Compared to sampling-based Monte Carlo algorithms (Gamerman and Lopes,, 2006), this optimization-based method is computationally efficient and scalable to large datasets. Variational inference turns the Bayesian inference problem into an optimization problem that approximates the posterior distribution by the closest member from a family of tractable distributions such as Gaussians. The optimal variational approximation is often solved by stochastic gradient descent algorithms (Hoffman et al.,, 2013), which are able to handle massive datasets.

Successful estimation of the posterior relies on determining an appropriate variational family. Typically, the multivariate Gaussian family is often used as it provides notable tractability and simplicity (Opper and Archambeau,, 2009; Ormerod and Wand,, 2012). However, estimating a full covariance matrix in Gaussian variational approximation is generally challenging, since the number of matrix elements grows quadratically with the number of latent variables and the computational cost quickly becomes too expensive. Consequently, special structures are applied to the covariance matrix for scalability. For example, the mean-field variational family assumes a diagonal structure for the covariance matrix (Xing et al.,, 2002; Gershman et al.,, 2012). However, this assumption of independence is inappropriate to many complicated data structures and tends to underestimate the variance in practice (Blei et al.,, 2017). An alternative is to factorize the covariance matrix using a low-dimensional representation. Ideally, such low-dimensional representation should be simple enough to maintain a tractable approximation, but expressive enough to capture well the posterior correlations. Various low-dimensional representation structures have been considered so far, where one common choice is the low-rank plus diagonal (LR+D) decomposition, see Nickisch and Seeger, (2009); Seeger, (2010); Guo et al., (2016); Miller et al., (2016). In particular, Ong et al., (2018) proposed a Gaussian variational approximation method with the factorized covariance structure (VAFC) via the LR+D representation. Much of our paper will build on the VAFC method.

Originated from the factor analysis model in statistics, the LR+D decomposition requires restrictions to guarantee global identification of the factors to be estimated. The topic has been studied for a long time, see Bekker, (1986); Grayson and Marsh, (1994); Grzebyk et al., (2004), to mention but a few. For example, the restrictions in Bai and Ng, (2013) require the underlying factor matrix to be orthogonal lower trapezoidal. In the Euclidean space, imposing extra constraints can result in a non-convex optimization problem because of the changes in the feasible space. Ong et al., (2018) retained Euclidean optimization by ignoring the full-rank restriction on the factor matrix. This might lead to an inefficient optimization procedure and make it likely trapped into a bad local mode because the factor matrix is not uniquely determined.

Evidence shows that a constrained optimization problem with non-Euclidean nature can be transformed into a geodesically convex problem in the Riemannian space (Edelman et al.,, 1998; Boumal,, 2014; Cherian and Sra,, 2016; Ferreira et al.,, 2018). Various manifolds have been applied for geometry-aware learning tasks (Huo et al.,, 2007; Fiori,, 2011; Cunningham and Ghahramani,, 2015). For instance, Grassmann manifolds are well-explored in learning low-dimensional representation (Ngo and Saad,, 2012; Dong et al.,, 2014; Zhang et al., 2018c, ), and Stiefel manifolds are applied widely in fields like pattern recognition (Wright et al.,, 2009; Browne and Mcnicholas,, 2014) and dimension reduction (Suzuki and Sugiyama,, 2013). Other structures push forward low-rank approximations while guaranteeing high accuracy and computational feasibility (Vandereycken and Vandewalle,, 2010; Sato and Iwai,, 2013; Zhou et al.,, 2015).

Optimization on a Riemannian manifold has attracted great attention in the last decade, especially for large-scale learning tasks (Absil et al.,, 2008; Bonnabel,, 2013; Sra and Hosseini,, 2015). Early research on stochastic gradient optimization on manifolds uses a fixed learning rate (Bonnabel,, 2013; Zhang et al.,, 2016; Zhang and Sra,, 2016; Ferreira et al.,, 2018). Some recent works have generalized the Euclidean-based adaptive learning schemes to the geometry-aware learning. Roy et al., (2018) proposed coordinate-wise adaptive operations with momentum based on RMSProp. Bécigneul and Ganea, (2018) generalized the popular Euclidean adaptive learning rate schemes, AdaGrad and Adam, to adaptive learning on manifolds. Kasai et al., (2019) adapted separate weight matrices corresponding to the row and column subspaces.

It is a natural idea to treat special structures in variational parameters, such as symmetric positive definiteness for covariance matrix parameters or orthogonality constraints, as non-Euclidean spaces and apply Riemannian optimization. Existing methods usually deal with such constraints in an ad hoc manner that is specifically designed for the constraints under consideration (van den Berg et al.,, 2018). It is worth mentioning that the natural gradient (Honkela et al.,, 2007; Mishkin et al.,, 2018; Zhang et al., 2018b, ; Tran et al.,, 2019) has been proposed as a speed-up gradient descent algorithm. The natural gradient, although operating on the Euclidean space, is a geometric object that is the steepest descent direction of the cost function while taking into account the geometry structure of the variational family.

In this paper, we develop stochastic Gaussian variational approximation methods that guarantee the unique optimal variational parameters in the covariance matrix factorization. We do so by imposing constraints on the variational parameters and exploit their geometry to perform manifold-assisted optimization. More specifically, we factorize the covariance matrix using an LR+D structure, and ensure the uniqueness of that structure by introducing two manifold constraints. The main contributions of this work are:

• •

We combine the idea of Riemannian optimization with Euclidean adaptive stochastic gradients techniques, and propose two adaptive learning algorithms, namely RGD-RMSProp and RGD-AdaDelta, for stochastic optimization on manifolds. The two new optimizers enable computational speed-up in large-scale variational inference learning.

• •

We impose the Stiefel and Grassmann manifold constraints on the parameters in the LR+D structure for the Gaussian variational approximation. The optimization on the Grassmannian equivalence classes solves the non-identifiability issue arising from the LR+D structure. It is possible to generalize the method to other constrained variational families such as the Sylvester normalizing flow of van den Berg et al., (2018).

• •

We implement our methods on both low and high-dimensional cases to empirically test their performance. Our methods outperform baseline variational inference methods in terms of approximation accuracy and computational speed.

The paper is organized as follows. Section 2 reviews the VAFC method, Stiefel and Grassmann manifolds, and Riemannian gradient descent. Section 3 details our Riemannian stochastic gradient schemes for Gaussian variational approximation. Section 4 empirically examines the performance of our methods. Section 5 concludes the paper with potential future extensions.

2 Preliminary

This section introduces the key concepts and sets up the notation. We start from Section 2.1 to formulate the objective function in variational inference. The low rank plus diagonal (LR+D) factorization on the covariance matrix is discussed and its non-identifiability is identified. Two related manifolds, the Stiefel and Grassmann manifolds, are presented in Section 2.2 and 2.3, respectively. Section 2.4 describes the general Riemannian stochastic gradient decent method.

2.1 Gaussian Variational Approximation

Consider the Gaussian variational approximation problem for Bayesian inference. Let $q_{\lambda}(\boldsymbol{\theta})$ denote a member of the variational family where $\lambda$ denotes the variational parameters. For example, if the variational family is multivariate Gaussians, then $\lambda$ consists of a mean vector and a covariance matrix. The objective is to find the best approximation $q_{\lambda^{*}}(\boldsymbol{\theta})$ to the posterior $p(\boldsymbol{\theta}|\mathbf{x})$ where the optimal variational parameter $\lambda^{*}$ is found by minimizing the Kullback-Leibler divergence between $q_{\lambda}(\boldsymbol{\theta})$ and $p(\boldsymbol{\theta}|\mathbf{x})$

[TABLE]

Due to the intractable evidence term $\log p(\mathbf{x})$, computing the KL-divergence is generally impossible. Fortunately, as this evidence term does not involve $q_{\lambda}(\boldsymbol{\theta})$, minimizing the KL-divergence is equivalent to maximizing the variational evidence lower bound (ELBO). Formally, we rewrite this ELBO objective function as

[TABLE]

where $h(\boldsymbol{\theta})=p(\boldsymbol{\theta})p(\mathbf{x}|\boldsymbol{\theta})$ with the prior $p(\boldsymbol{\theta})$ and the model likelihood $p(\mathbf{x}|\boldsymbol{\theta})$. We denote by $m$ the dimension of $\lambda$, which can be considerably large. To find the optimal variational parameters from (1), we apply the stochastic gradient descent method by repeating

[TABLE]

until the sequence converges. Here $\alpha_{t}$ is the learning step size at the $t$-th iteration, and $\widehat{\nabla_{\lambda}\mathcal{L}}$ is an unbiased approximation of the exact gradient $\nabla_{\lambda}\mathcal{L}\left(\lambda^{(t)}\right)$.

In the context of Gaussian variational approximation, the variational distribution to be estimated is $q_{\lambda}(\boldsymbol{\theta})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, $\lambda=(\boldsymbol{\mu},\boldsymbol{\Sigma})$. To preserve some off-diagonal covariance structure with a small computational cost, Ong et al., (2018) factorize the full matrix using the LR+D representation, i.e.,

[TABLE]

where $\mathbf{B}$ is an $m\times p$ full rank matrix with $p\leq m$ and $\mathbf{D}_{2}$ is a diagonal matrix with non-zero diagonal elements in vector $\mathbf{d}_{2}=(d_{21},...,d_{2m})^{\top}$. The reason for using the subscript in $\mathbf{D}_{2}$ will become clear shortly. We use the number of factors $p$ to control the complexity of $q_{\lambda}(\boldsymbol{\theta})$. With a small $p$ value, both the storage and computational cost are reduced on approximating the dense covariance matrix $\boldsymbol{\Sigma}$, so that it is possible to solve large-scale problems. Specifically, when $p=0$, the model degenerates to the mean-field variational structure. With the factorized form (2), the variational parameter $\lambda$ collects $\{\boldsymbol{\mu},\mathbf{B},\mathbf{d}_{2}\}$. We rewrite (1) as

[TABLE]

We ignore the constant term $\frac{m}{2}(\log(2\pi)+1)$ since it is independent from $\lambda$, the parameters to be estimated. The final objective function is

[TABLE]

To obtain efficient stochastic gradient estimates, we use the reparameterization trick (Kingma and Welling,, 2014). Let $\mathbf{z}\in\mathbb{R}^{p}$ and $\boldsymbol{\epsilon}\in\mathbb{R}^{m}$ be two standard Gaussian variables, i.e., $(\mathbf{z},\boldsymbol{\epsilon})\sim f(\mathbf{z},\boldsymbol{\epsilon})=\mathcal{N}(\mathbf{0},\mathbb{I})$. The random variable $\boldsymbol{\theta}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ can be reparameterized as

[TABLE]

where $\circ$ indicates the Hadamard product. The objective function (3) then becomes

[TABLE]

The advantage of (4) over (3) is that, when the expectation is approximated by a sampling method, the samples $(\mathbf{z},\boldsymbol{\epsilon})$ are independent of $\lambda$, which can be explained as a control variate scheme (Paisley et al.,, 2012; Xu et al.,, 2019).

Such LR+D factorization for Gaussian variational approximation has many attractive properties and enjoys competitive performance (Ong et al.,, 2018). However, this VAFC method of Ong et al., (2018) suffers from the non-identifiability issue. Consider any orthogonal matrix $\mathbf{Q}$ of size $p\times p$, it is easy to see that

[TABLE]

as $\mathbf{Q}\mathbf{z}$ still follows the standard Gaussian distribution. Then for any orthogonal matrix $\mathbf{Q}$,

[TABLE]

Hence, if $\lambda(\boldsymbol{\mu},\mathbf{B},\mathbf{D}_{2})$ is the variational solution, so is $\lambda(\boldsymbol{\mu},\mathbf{B}\mathbf{Q},\mathbf{D}_{2})$. This non-identifiability may bring in certain optimization difficulties. While there are potential practical advantages that determining a small number of factors $p$ helps reduce the probability of occurring the identification problem, the method, from a purely statistical point of view, still dissatisfies the non-ignorable arbitrariness (Bartholomew et al.,, 2011).

To avoid the potential non-identifiability issue, restrictions on the factor matrix $\mathbf{B}$ are often imposed. For example, $\mathbf{B}$ can be a column full-rank lower triangular matrix with positive diagonal elements (Ong et al.,, 2018; Tran et al.,, 2018). Unfortunately, without careful treatment, the updating process might fail to guarantee such constraints on $\mathbf{B}$. Consider a stochastic gradient process with a fixed learning rate $\alpha$ = 0.01. Suppose we have $\mathbf{B}^{(t)}$ = $\begin{bmatrix}1&0\\ 2&1\end{bmatrix}$ and a calculated gradient $\frac{\partial\mathcal{L}(\mathbf{B})}{\partial\mathbf{B}}$ = $\begin{bmatrix}-2&0\\ -50&100\end{bmatrix}$. This gives the updated $\mathbf{B}^{(t+1)}=\begin{bmatrix}1.02&0\\ 2.5&0\end{bmatrix}$, which is a non-feasible solution to the optimization problem because the $\mathbf{B}^{(t+1)}$ is rank-deficient.

To this end, we propose an optimization framework on manifolds to guarantee a feasible solution that respects the necessary restrictions. We do so by utilizing the concepts of Stiefel manifold and Grassmann manifold that we review next.

Remark 1*.*

We do not consider the structure in Journée et al., (2010) as their manifold is not complete. The geometry is only suggested for Fixed-rank situations, where the rank of $\mathbf{BB}^{\top}$ equals exactly to some $k$ ($k\leq m$). The method also failed in our empirical experiments. Evidence is provided in the first row of Figure 3 where the algorithm does not converge properly.

2.2 Stiefel Manifold Constraint

Definition 1** (The Stiefel Manifold (Absil et al.,, 2008)).**

Let $p\leq m$, the Stiefel manifold $\mathcal{S}(p,m)$ is the set of $m\times p$-dimensional matrices consisting of orthonormal columns. That is

[TABLE]

A Stiefel manifold $\mathcal{S}(p,m)$ consists of all $p$ orthogonormal basis vectors that are rigidly connected to each other and represents all the equivalent results, i.e., a Stiefel manifold point $\mathbf{B}\in\mathcal{S}(p,m)$ includes all the equivalent $\mathbf{BQ}$ by a set of orthonormal bases. Formally, the covariance matrix in the Gaussian variational approximation $q_{\lambda}(\boldsymbol{\theta})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ is factorized as

[TABLE]

where $\mathbf{D}_{1}$ and $\mathbf{D}_{2}$ are diagonal matrices of size $p\times p$ and $m\times m$. We introduce $\mathbf{D}_{1}$ to solve the non-identifiability issue for the parameter $\mathbf{B}$. For the new covariance in (6), its corresponding Gaussian $q_{\lambda}(\boldsymbol{\theta})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ can be found by a new deterministic transformation

[TABLE]

where $(\mathbf{z},\boldsymbol{\epsilon})$ follows the standard $(p+m)$-dimensional Gaussian distribution $f(\mathbf{z},\boldsymbol{\epsilon})$.

Remark 2*.*

We adopt the scaling factor $\mathbf{D}_{1}$ in (6) to compensate for the loss of scale from the orthogonality requirement. Also, the factor effectively prevents the non-identifiability of $\mathbf{B}$.

Applying the reparameterization trick in (7) will update the objective function (3) to

[TABLE]

where the expectation $\mathbb{E}_{f}$ is with respect to the standard Gaussian distribution $f(\mathbf{z},\boldsymbol{\epsilon})=\mathcal{N}(0,\mathbb{I})$.

The Stiefel manifold constraint guarantees a full column rank solution with an orthonormal $\mathbf{B}$. The stochastic gradient estimates on $\mathbf{B}$ will be conducted with its manifold geometry.

2.3 Grassmann Manifold Constraint

Definition 2** (The Grassmann Manifold (Absil et al.,, 2008)).**

The Grassmann manifold, denoted by $\mathcal{G}(p,m)$, consists of all the $p$-dimensional subspaces in $\mathbb{R}^{m}$ ($p\leq m$). Any Grassmann point on $\mathcal{G}(p,m)$ can be represented by an $m\times p$ matrix $\mathbf{B}$ with orthonormal columns, that is, $\mathbf{B}^{\top}\mathbf{B}=\mathbb{I}_{p}$.

A Grassmann point on $\mathcal{G}(p,m)$ is the subspace spanned by the columns of the factor matrix $\mathbf{B}$. It can be considered as an equivalence class over Stiefel manifold $\mathcal{S}(p,m)$, i.e., for a representative $\mathbf{B}\in\mathcal{S}(p,m)$, it has an equivalent class

[TABLE]

where $\mathcal{O}(p)$ is the orthogonal group of order $p$ (Absil et al.,, 2008).

The property in (5) shows that the ELBO $\mathcal{L}(\lambda)$ in (4) is a well defined function on $\mathbb{R}^{m}\otimes\mathcal{G}(p,m)\otimes\mathbb{R}^{m}$ where $\boldsymbol{\mu}\in\mathbb{R}^{m}$, $\mathbf{B}\in\mathcal{G}(p,m)$ and $\mathbf{d}_{2}\in\mathbb{R}^{m}$. We then rewrite (3) as

[TABLE]

The equation is arranged similarly to (4), except $\mathbf{B}$ is a Grassmann point, i.e., a representative Stiefel matrix of its equivalence class. The optimal solution $\mathbf{B^{*}}$ from (9) is a Grassmann point with all equivalence solutions $[\mathbf{B^{*}}]$, and it is free from the identification problem. In contrast, (4) searches in the Euclidean space for a single solution $\mathbf{B}$ (among many) instead of its unique equivalence class solution $[\mathbf{B}]$.

Due to the manifold constraints, the conventional gradient update rules working on the Euclidean space fail to guarantee convergence for the objective functions in (8) and (9). Instead, we will apply the Riemannian gradient descent rules that we review next.

2.4 Riemannian Gradient Descent

We have formulated two new objective functions with the Stiefel and the Grassmann manifold constraints respectively. This section reviews how an optimization problem on manifolds is solved by standard Riemannian optimization algorithms. We start from the basic line-search method that involves a one-step gradient descent. The notion of retraction mapping will be introduced to enable updating on a manifold. We then review the gradient momentum methods that help accelerate convergence. By vector transport, we accumulate tangent vectors from different vector spaces onto the same space. Further, we show how a point on a manifold is connected to a tangent space by projection.

For a differentiable function $\mathcal{L}(x)$ in an Euclidean space, the line-search method defines

[TABLE]

where $\xi_{k}\in\mathbb{R}^{m}$ is the search direction and $\alpha_{k}\in\mathbb{R}_{+}$ is the step size. The gradient, however, relies on the vector space structure of $\mathbb{R}^{m}$, and this dependency is not fundamental for an objective function defined on nonlinear manifolds. Instead, for objective functions defined on a manifold $\mathcal{M}$, the general update rule is

[TABLE]

where $\xi_{k}$ is a tangent vector to $\mathcal{M}$ at $x_{k}$. We define $\mathrm{R}_{x}$ as a retraction at ${x}$ to map $\xi_{k}$ smoothly from the tangent space $T_{x}\mathcal{M}$ onto $\mathcal{M}$. The retraction mapping guarantees $x$ is continuously updated along the steepest descent (Absil et al.,, 2008), i.e., its moving direction is on a tangent vector, while the point is staying on the manifold. Figure. 1 visualizes the process of a retraction mapping to preserve the gradient at $x$.

In the line-search method, the tangent vector $\xi$ is the steepest-descent direction of the objective function, i.e., $\xi_{k}=-\text{grad}\mathcal{L}(x_{k})$ the Riemannian gradient of $\mathcal{L}$ at $x_{k}$. Many modern gradient descent methods include momentum to dampen oscillations by accumulating the stochastic gradient descents. In Euclidean space, the search direction is produced by summing up $-\text{grad}\mathcal{L}(x_{k})$ and $\xi_{k-1}$. On a nonlinear manifold, the two tangent vectors with respect to $x_{k}$ and $x_{k-1}$ do not belong to the same subspace, thus they cannot be combined arbitrarily. This issue can be fixed by applying a vector transport to bring $\xi_{k-1}$ to the tangent space at $x_{k}$. The search direction is then

[TABLE]

where $\Gamma_{x_{k-1}\rightarrow x_{k}}$ is the parallel transport that maps a vector $\xi_{k-1}\in T_{x_{k-1}}\mathcal{M}$ to a vector in $T_{x_{k}}\mathcal{M}$.

For manifolds that can be embedded in the ambient Euclidean space, their Riemannian gradient is found by projecting the Euclidean gradient in the ambient space of the manifold onto the tangent space at the point $x$. The way to implement the projection operator depends on manifolds. Table 1 and Table 2 summarize the mapping formulas of Stiefel and Grassmann manifolds.

Note that we define the symmetric operation of $\mathbf{B}^{\top}\mathbf{U}$ as $\text{sym}(\mathbf{B}^{\top}\mathbf{U})=\frac{1}{2}(\mathbf{B}^{\top}\mathbf{U}+\mathbf{U}^{\top}\mathbf{B})$. The retraction $\mathrm{R}_{\mathbf{B}}(\mathbf{U})$ on the Grassmann manifold can be implemented as the polarization $\text{polar}(\mathbf{B}+\mathbf{U})$. When $\mathbf{U}\neq 0$, the matrix $\mathbf{B}+\mathbf{U}$ is not orthogonal. To retain a representative on Grassmann manifold, we use the polar decomposition to avoid ill-conditioning.

In a nutshell, the Riemannian gradient descent algorithm is applied to update a point $x$ on a manifold. The search direction is found by the Riemannian gradient over the objective function, which projects the Euclidean gradient onto the tangent space. For vectors from different tangent spaces, parallel transport is adopted to accumulate them. Finally, the searching variable is mapped back to the manifold space by the retraction operator.

3 Approximation with Riemannian Operations

We have reformulated the objective functions in (8) and (9) in Section 2.2 and 2.3, respectively. The task is to find the optimal variational parameters $\lambda^{*}=\{\boldsymbol{\mu}^{*},\mathbf{B}^{*},\mathbf{D}^{*}_{1},\mathbf{D}^{*}_{2}\}$ that maximize the objective function $\mathcal{L}(\lambda)$. For the unconstrained parameters $\{\boldsymbol{\mu},\mathbf{D}_{1},\mathbf{D}_{2}\}$, we apply the stochastic gradient descent (SGD) update rules in the Euclidean space, see Section 3.1. In Section 3.2, we introduce the Riemannian (stochastic) gradient descent (RGD) update rules for the constrained parameter $\mathbf{B}$. We consider four different variants, where the first two are conventional methods for Riemannian optimization, and the last two adaptive update rules are first developed in this paper.

3.1 Euclidean Gradient of the Log-Likelihood Function

We first calculate the Euclidean partial derivative of the objective function with respect to each parameter. The unconstrained parameters can be updated directly with the SGD rules, while the constrained parameters require these derivatives as the input to implement the RGD rules. In the case of Stiefel constraint in (8), we calculate the derivatives as follow

[TABLE]

where $\nabla_{\boldsymbol{\theta}}\log h(\boldsymbol{\mu}+\mathbf{BD_{1}z}+\mathbf{d}_{2}\circ\boldsymbol{\epsilon})=\nabla_{\boldsymbol{\theta}}\log h(\boldsymbol{\theta})$, which is easy to calculate, e.g., by automatic differentiation. For high-dimensional problems, it is challenging to compute the inverse matrix $(\mathbf{BD^{2}_{1}B}^{\top}+\mathbf{D}^{2}_{2})^{-1}$ efficiently. Instead, we follow Ong et al., (2018) and apply the Woodbury identity (Woodbury,, 1950) to convert the inverse of an $m\times m$ matrix to an inverse of a $p\times p$ matrix as follows

[TABLE]

Similarly, for (9) with the Grassmann manifold constraint, we have

[TABLE]

3.2 Riemannian Stochastic Gradient Descent

We now describe four gradient-based updating rules under the RGD framework. We start with the vanilla RGD method (Bonnabel,, 2013) that uses a fixed learning rate. To accelerate convergence, we consider the second RGD update rule with momentum (Roy and Harandi,, 2017). We propose two adaptive learning rules for RGD, constrained RMSProp rule and the simulated AdaDelta rule. Both methods use a new learning rate for every parameter at every time step, and this learning rate is determined automatically during the update process. These two adaptive learning rules for Riemannian stochastic optimization, to the best of our knowledge, are first developed in this paper. In the following, we will use the notation $\mathcal{L}(\mathbf{B})$ for $\mathcal{L}(\lambda)$ where we consider the optimization over the manifold variable $\mathbf{B}$.

3.2.1 Basic RGD method with Fixed Learning Rate

The first update scheme follows the procedure in Bonnabel, (2013) and extends the vanilla stochastic gradient descent algorithm to a Riemannian manifold. As reviewed in Section 2.4, both Stiefel and Grassmann manifolds can be embedded in the ambient Euclidean space. The Riemannian gradient of the objective function at a point $\mathbf{B}$ is calculated by projecting its Euclidean gradient $\nabla_{\mathbf{B}}\mathcal{L}(\mathbf{B})$ from the ambient space to the tangent space, i.e.,

[TABLE]

The explicit projection mapping function $\pi_{\mathbf{B}}$ of the two manifold constraints are defined in Table 1 and Table 2. Once the Riemannian gradient is prepared, the RGD update rule is applied

[TABLE]

where $\mathrm{R}_{\mathbf{B}}(\cdot)$ is the retraction operator that maps the point on tangent space onto the manifold, and $\alpha$ is a fixed learning rate. We refer to this basic RGD method as RGD-Basic.

3.2.2 The Constrained SGD with Momentum Rule

The basic RGD updating rule in (18) provides an intuitive way of optimization over manifolds. However, that line-search method over the steepest-descent direction could be inefficient around local optima where the cost surface is usually steeper in one dimension than in another. Based on the momentum idea from the SGD literature, , Roy and Harandi, (2017) proposed Constrained RGD with Momentum (cRGD-M). Using their method, the manifold-constrained variable $\mathbf{B}$ is updated by

[TABLE]

where $0<\beta$ is the momentum constant and $\Gamma_{\mathbf{B}_{1}\rightarrow\mathbf{B}_{2}}(\mathbf{U})$ is the vector transport which maps the tangent vector $\mathbf{U}$ from the tangent space at $\mathbf{B}_{1}$ to that at $\mathbf{B}_{2}$. The implementation of this vector transport can be found in Table 1 and Table 2.

3.2.3 The Modified Constrained RMSProp Rule

Tuning a proper learning rate can be tricky especially for large-scale learning tasks. Also, the same learning rate might not apply to all parameter updates during the training process. Roy et al., (2018) considered the parameter-wise adaptive learning rate for the sparse data structure so that the optimizer gains a rapid convergence speed. Their method is a constrained RMSProp rule that ensures $\mathbf{B}^{\top}\mathbf{B}=\mathbb{I}$. However, this update rule breaks when the combination of a parallel transport and a projection on the tangent space becomes negative. Instead, we modify the update rule in Roy et al., (2018) and propose RGD-RMSProp as follows. We include the absolute operator over the element $E(g_{\mathbf{B}}^{2})$ and compensate its sign after the square root operation.

[TABLE]

where

[TABLE]

with $\odot$ as the element-wise product of two matrices and $\nabla\mathcal{L}$ as the Euclidean gradient of the objective function. For the empirical experiments, we suggest $\alpha=0.05,\beta=0.95$ and $\epsilon=10^{-6}$.

Remark 3*.*

Our RGD-RMSProp can be regarded as the basic RGD with a constrained normalized Euclidean gradient. The method differs from the Euclidean RMSProp (Hinton et al.,, 2016) and requires sophisticated mappings to accelerate multiple tangent vectors.

3.2.4 The Simulated AdaDelta Rule

Compared to the basic RGD, the RGD-RMSProp includes the adaptive learning rate scheme that effectively speeds up convergence. However, (21) has two limitations. First, selecting the learning rate $\alpha$ can be challenging. Second, when updating the parameter $\mathbf{B}$, the hypothetical units are related to the gradients $\nabla\mathcal{L}(\mathbf{B})$ rather than the parameters (Zeiler,, 2012). To fix the two issues, we propose using element-wise division in (22). The update rule is arranged as follows

[TABLE]

At each step, the learning rate is self-updated with the previous performance of the gradient. The update rules above could be considered as the Riemannian version of its Euclidean counterpart (Zeiler,, 2012). We refer to this method as RGD-AdaDelta.

The overall algorithm for the Manifold assisted Gaussian Variational Approximation algorithm is summarized in Algorithm 1. For notational simplicity, in the following, we use S for Stiefel and G for Grassmann.

4 Experiment

In this section, we investigate the performance of our proposed Gaussian variational approximation methods. There are eight of them: the Gaussian variational approximation with Stiefel manifold in Section 2.2 and the Gaussian variational approximation with Grassmann manifold in Section 2.3, each is estimated with four manifold learning rules described in Section 3.2. We start from Section 4.1 with a toy experiment where the data is generated from a non-linear regression model with outliers. We evaluate all eight proposed methods’ sensitivity to the different numbers of factors and report the fitting performance as well as the convergence speed. We also show the infeasibility of the naive low-rank representation structure mentioned in Section 2.1. In Section 4.2 and 4.3, we demonstrate the methods’ prediction power with low-dimensional ($m<n$) and high-dimensional ($m\gg n$) datasets, respectively. The accuracy is compared against three baseline methods, including VAFC (Ong et al.,, 2018), SLANG (Mishkin et al.,, 2018) and the conventional mean-field variational inference method. SLANG is a variational inference algorithm that uses the natural gradient and the mean-field variational inference method is a special case of our framework with the factor number $p=0$. All experiments are implemented in Matlab on a laptop with an Intel Core 64-bit i5-6600 3.3GHz CPU and 8G RAM.

4.1 Simulation Study

We follow Ong et al., (2018) and generate data to investigate our methods’ performance. We are specifically interested in their convergence speed, fitting performance, and their sensitivity to different numbers of factors $p$.

Consider a regression model

[TABLE]

where $g(\mathbf{x})=\sin(4\pi\mathbf{x})$ is the underlying function to be approximated. We generate the $\mathbf{x}_{i}$ from $U(0,1)$ and $\epsilon_{i}$ from the mixture

[TABLE]

We set $\sigma^{2}=1$ and $a=5$. The parameter $a$ controls the occurrence of extreme values. The methods’ degree of robustness can then be justified by setting different $a$s. A visualization of the underlying function and the generated 250 samples is given in Figure 2.

To approximate the underlying function $g$, we consider the penalized B-spline (P-spline) regression model (Lang and Brezger,, 2004). Denote by $b(\mathbf{x})$ a basis of P-splines with $k$ basis functions $b(\mathbf{x})=(b_{1}(\mathbf{x}),b_{2}(\mathbf{x}),...,b_{k}(\mathbf{x}))^{\top}$, and $\boldsymbol{\beta}=(\beta_{1},\beta_{2},...,\beta_{k})^{\top}$ the corresponding $k$ coefficients. We use $k=25$ in this example. The B-spline regression model is

[TABLE]

The model has a total number of $k+3$ random variables to infer. Aside from the aforementioned $k$ coefficients $\boldsymbol{\beta}$, there are three additional parameters $(\tau^{2},\psi,\sigma^{2})$ for describing the density of noise $\boldsymbol{\epsilon}$. Below we give full details to construct them under the Bayesian learning framework.

Ong et al., (2018) use the likelihood

[TABLE]

and the following prior for $\beta$

[TABLE]

with $P(\psi)^{-1}$ the matrix formed by the AR(1) model with the persistence coefficient $\psi$. The prior of the hyperparameters $\tau^{2},\psi$ is

[TABLE]

where $p(1/\tau^{2})$ is the Weibull density

[TABLE]

Finally, the prior for $\sigma^{2}$ is an inverse Gamma distribution with shape parameter $\alpha_{0}$ and scale parameter $\beta_{0}$

[TABLE]

The set of model parameters is $\boldsymbol{\theta}=\{\boldsymbol{\beta},\tau^{2},\psi,\sigma^{2}\}$. The Gaussian variational approximation methods for the posterior of $\boldsymbol{\theta}$ are obtained by the procedures outlined in Algorithm 1. We are interested in their convergence speed, accuracy, and their sensitivity to different numbers of factors $p$. The factorization of $\boldsymbol{\Sigma}$ follows (6) for Stiefel manifold constraints, and (2) for Grassmann manifold constraints. For the optimizers, we specify the learning rates $\alpha=0.05,\beta=0.95$ and $\epsilon=10^{-6}$.

The lower bounds of the proposed methods are visualized in Figure 3 with $p=0,3,6$ and $10$ factors. We only plot the first 1000 iterations for the fixed-rank and Grassmann methods and 2000 iterations for the methods with Stiefel manifold, as it takes more iterations to converge. All values are averaged from 10 repetitions.

The first row shows the results from the Fixed-rank geometry methods. Their failure to converge properly supports the claims in Remark 1. In the second row, the four Stiefel manifold-based methods share a similar convergence speed. The selection of factor $p$ does not dramatically influence the convergence speed, as the lower bounds for the different $p$ choices do not significantly differ from each other. Also, including extra covariance information does not benefit the model approximation as expected. As displayed in the subplots, not all the lower bound lines are above the magenta dashed line (factor $p=0$). In other words, more complex factor structures do not always outperform the naive independent covariance structure when implementing the Stiefel manifold method. The third row shows the lower bounds from the Grassmann manifold method, which tells a different story. The lines are generally smoother and stabler than the ones in the second row, and attain the maximum in the first 1000 iterations. In addition, the approximation process of the Grassmann methods can benefit from the covariance structure. For all four methods, selecting a larger factor of $p$ always results in a better lower bound curve. The estimation results of the proposed Gaussian variational approximation methods are evaluated in Figure. 4.

All 32 regression models (with 4 factor values and 8 methods) capture the underlying function $g(x)$ reasonably well. The results with $p>0$ factors are generally less sensitive to the extreme values than the baseline mean-field methods with a diagonal covariance ($p=0$). Since the samples are distributed less dense at $[0,0.2]$ than at $[0.5,0.7]$, the estimated lines are slimmer than the ground truth (green line) in the former region but wider in the latter region. Also, within $[0.5,0.7]$ there are more observations above the green line with considerably large values, see Figure 2. Consequently, higher peaks are approximated than the true function. The four Gaussian variational approximation methods with the Grassmann manifold constraint are more robust to extreme values than the ones with the Stiefel manifold. The fitted lines around stationary points are more similar to the ground truth. Also, the factorized methods outperform the baseline mean-field method, especially when working around the outliers. In terms of the four optimizers, they obtain similar performance with Stiefel and Grassmann geometry. We will leave further comparison and discussion in Section 4.2 and 4.3 with complex data structure.

Based on what we have observed from these toy experiments, we conclude that both the Stiefel and Grassmann manifolds provide accurate results at a considerably fast speed, while the conventional Fixed-rank geometry fails to handle the constrained optimization. In addition, the Grassmann methods outperform the Stiefel methods in terms of both convergence rates and approximation accuracy.

In the next two sections, we will consider empirical datasets including both low-dimensional and high-dimensional settings to examine our Gaussian variational approximation methods with more complex structures.

4.2 Low Dimensional Predictive Inference

We first evaluate our methods on the logistic regression for a low dimensional binary classification problem. We select the ionosphere data from the UCI Machine Learning Repository (Dua and Graff,, 2017) with 351 samples. The task is to classify 225 positive signals that could identify an ionosphere structure and 126 negative instances that are insensitive to the free electrons. We follow the same data processing operations as Ong et al., (2018) and expand the attributes from 34 to 111. Given $\widehat{\mathbf{x}}_{i}\in\mathbb{R}^{q}$, the likelihood of the response $y_{i}\in\{-1,1\}$ is

[TABLE]

where $\mathbf{x}_{i}=[1,\widehat{\mathbf{x}}^{\top}_{i}]^{\top}$ and $\boldsymbol{\beta}$ denotes the coefficient vector. We use the Gaussian prior for $\boldsymbol{\beta}$. The prediction performance of five-fold cross-validation is recorded in Table 3. The best results are highlighted in terms of the training error, test error and training time.

The overall accuracy of our proposed methods improves slightly from the two benchmarks VAFC (Ong et al.,, 2018) and mean-field models, but outperforms SLANG (Mishkin et al.,, 2018) significantly. All eight proposed methods have a lower test error, where around 5 samples are misclassified among 70 test instances. In particular, the Stiefel methods obtain both lower training error and test error, while the Grassmann methods further cut down the test error with a higher training error. Therefore, we conclude that the imposed manifold constraints, especially the Grassmann manifold constraint, make the methods less likely to overfit the data. The computational times of all eight methods are similar and similar to that of VAFC, and much less than that of SLANG. In terms of the optimizers, the two new methods (RGD-RMSProp and RGD-AdaDelta) outperforms the two conventional update rules (RGD-Basic and cRGD-M) in both Stiefel and Grassmann settings. In particular, the RGD-RMSProp methods obtain the lowest test error with a fast speed.

Figure 5 compares the estimation of lower bounds for the 8 proposed methods with $p=4$ factors over $5000$ iterations. Similar to the toy example, the shape of each curve depends on the manifold constraints. All the 8 models shrink their lower bounds quickly in the first $2000$ iterations. However, it takes longer for Stiefel methods to flatten the curve, and the oscillation is more significant than the Grassmann methods along the update process. The reason is that the Stiefel methods have an extra scaling factor $\mathbf{D}_{1}$ to update while the Grassmann methods only take care of $\mathbf{B}$ and $\mathbf{D}_{2}$. For different optimizers on the same manifold constraints, the latter two optimizers (RGD-RMSProp and RGD-AdaDelta) update slightly faster than the former two (RGD-Basic and cRGD-M), but the advantage does not distinguish clearly.

4.3 High Dimensional Predictive Inference

In the last experiment, we consider the binary Leukemia Cancer dataset with $7120$ predictors available at http://www.ntu.edu.sg/home/elhchen/data.htm. This classification problem contains 38 training instances and 34 testing instances. We evaluate the methods’ scalability with $m\gg n$. The Horseshoe prior (Carvalho et al.,, 2010) is applied due to the sparse data structure. The predictive accuracy is reported in Table 4. The performance of the 8 methods is compared against the mean-field variational approximation and VAFC. We don’t report the results from SLANG as it fails to converge after 96 hours.

As the sample size is relatively small, the proposed methods achieve similar performance in terms of predictive accuracy. Still, they all beat the two benchmark models. Specifically, the two negative samples are identically misclassified among all methods, which are wrongly predicted as positive. One of the possible reasons is that the training dataset is unbalanced with 27 positive samples against 11 negative samples. In terms of the operating speed, the conventional mean-field method stands out with its trivial structure. However, the fast computation sacrifices model accuracy, where the test error is significantly higher than the other methods. For the rest of the models, the Grassmann models are faster than other structures. In terms of the update methods, the RGD-AdaDelta scheme is faster than other optimization methods.

5 Discussion and Conclusion

This paper proposes two new manifolds-assisted optimization methods to resolve the identification issue of matrix factorization in the Gaussian variational approximation. We also develop two adaptive Riemannian gradient descent schemes, RGD-RMSProp and RGD-AdaDelta, to accelerate the optimization process. The proposed methods are assessed in both low- and high-dimensional situations. The experiments demonstrate the better performance of the proposed methods than the benchmark mean-field, VAFC, and SLANG models in terms of prediction accuracy and convergence speed.

It should be aware that the overhead cost of optimization on manifolds is usually higher than its Euclidean counterpart. This is because of the extra calculation of the retraction map and vector transport. The current research restricts the distribution of the latent variables to the multivariate Gaussian. Extension to other variational distributions, such as the exponential family, would be a potential research direction.