Accelerated Stochastic Quasi-Newton Optimization on Riemann Manifolds

TL;DR

This paper introduces an accelerated stochastic quasi-Newton method on Riemannian manifolds that achieves fast convergence without line searches, supported by new theoretical proofs and empirical validation on matrix and eigenvalue problems.

Contribution

It presents a novel L-BFGS algorithm on Riemannian manifolds with variance reduction, providing convergence proofs without curvature assumptions and demonstrating effectiveness in practical experiments.

Findings

Fast convergence on synthetic experiments

Effective computation of Karcher means and eigenvalues

Outperforms VR-PCA and Riemannian SVRG in tests

Abstract

We propose an L-BFGS optimization algorithm on Riemannian manifolds using minibatched stochastic variance reduction techniques for fast convergence with constant step sizes, without resorting to linesearch methods designed to satisfy Wolfe conditions. We provide a new convergence proof for strongly convex functions without using curvature conditions on the manifold, as well as a convergence discussion for nonconvex functions. We discuss a couple of ways to obtain the correction pairs used to calculate the product of the gradient with the inverse Hessian, and empirically demonstrate their use in synthetic experiments on computation of Karcher means for symmetric positive definite matrices and leading eigenvalues of large scale data matrices. We compare our method to VR-PCA for the latter experiment, along with Riemannian SVRG for both cases, and show strong convergence results for a range of datasets.

Equations

$$f(\gamma (t))\le tf(x)+(1-t)f(y)\forall t\in [0,1].$$

$$f(y)\ge f(x)+g_x({\nabla }_x,{\text{Exp}}_x^{-1}(y))+\frac{S}{2}\Vert {\text{Exp}}_x^{-1}(y){\Vert }^2.$$

$$\Vert \nabla (x)-{\Gamma }_{\gamma }\nabla (y)\Vert \le L\Vert {\text{Exp}}_x^{-1}(y)\Vert ,$$

$$f(y)\le f(x)+g_x({\nabla }_x,{\text{Exp}}_x^{-1}(y))+\frac{L}{2}\Vert {\text{Exp}}_x^{-1}(y){\Vert }^2.$$

$$\underset{x\in \mathcal{M}}{\min }f(x)=\frac{1}{N}\sum _{i=1}^N{f}_i(x),$$

$${\hat{B}}_r={\Gamma }_{\gamma }\left[{\hat{B}}_{r-1}-\frac{B_{r-1}(s_{r-1},\cdot ){\hat{B}}_{r-1}{s}_{r-1}}{B_{r-1}(s_{r-1},s_{r-1})}\right]{\Gamma }_{\gamma }^{-1},$$

$$H_r={\Gamma }_{\gamma }\left[G^{-1}{H}_{r-1}G+\frac{g_{x_{r-1}}(s_{r-1},\cdot )s_{r-1}}{y_{r-1}{s}_{r-1}}\right]{\Gamma }_{\gamma }^{-1},$$

$$\mathrm{tr}({\hat{B}}_r^{s+1})\le \mathrm{tr}({\hat{B}}_0^{s+1})+M\Lambda ,\det ({\hat{B}}_r^{s+1})$$

$$\mathbb{E}\left[f(x^{T+1})-f(x^{\ast })\right]\le L{S}^{-1}{\beta }^T\mathbb{E}\left[f(x^0)-f(x^{\ast })\right],$$

$$\mathbb{E}\Vert \nabla f(x^T){\Vert }^2\le (Kϵ{)}^{-1}L{\eta }_2^{{\alpha }_1}{\zeta }^{{\alpha }_2}\left(f(x^0)-f(x^{\ast })\right).$$

$$\underset{\mathbf{W}⪰0}{\mathrm{argmin}}\left\{\sum _{i=1}^N\Vert \log \left({\mathbf{W}}^{-\text{nicefrac}12}{X}_i{\mathbf{W}}^{-\text{nicefrac}12}\right){\Vert }_F^2\right\}.$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^d:{\mathbf{z}}^T\mathbf{z}=1}{\min }-{\mathbf{z}}^T\left(\frac{1}{N}\sum _{i=1}^N{d}_i{d}_i^T\right)\mathbf{z},$$

$$\underset{x\in \mathcal{M}}{\min }f(x)=\frac{1}{N}\sum _{i=1}^N{f}_i(x),$$

$${\hat{B}}_r={\Gamma }_{\gamma }\left[{\hat{B}}_{r-1}-\frac{B_{r-1}(s_{r-1},\cdot ){\hat{B}}_{r-1}{s}_{r-1}}{B_{r-1}(s_{r-1},s_{r-1})}\right]{\Gamma }_{\gamma }^{-1},$$

$$H_r={\Gamma }_{\gamma }\left[G^{-1}{H}_{r-1}G+\frac{g_{x_{r-1}}(s_{r-1},\cdot )s_{r-1}}{y_{r-1}{s}_{r-1}}\right]{\Gamma }_{\gamma }^{-1},$$

$$\mathrm{tr}({\hat{B}}_r^{s+1})\le \mathrm{tr}({\hat{B}}_0^{s+1})+M\Lambda ,\det ({\hat{B}}_r^{s+1})$$

$$G_r(\cdot ,\cdot )=\underset{0}{\overset{1}{\int }}{D}^2[f(t{z}_r)](\cdot ,\cdot )dt,$$

$$\frac{y_r{z}_r}{\Vert z_r{\Vert }^2}=\frac{G_r(s_r,s_r)}{\Vert z_r{\Vert }^2}\ge \lambda ,\frac{\Vert y_r{\Vert }^2}{y_r{z}_r}\le \Lambda .$$

$$\mathrm{tr}({\hat{B}}_r)$$

$$\le \mathrm{tr}({\Gamma }_{\gamma }{\hat{B}}_{r-1}{\Gamma }_{\gamma }^{-1})+\frac{\Vert {\Gamma }_{\gamma }{\hat{y}}_{r-1}{\Vert }^2}{y_{r-1}{s}_{r-1}}$$

$$\le \mathrm{tr}(B_0)+M\Lambda .$$

$$\det ({\hat{B}}_r)=$$

$$=$$

$$\ge$$

$$\det ({\hat{B}}_r)\ge \det (B_0)\frac{{\lambda }^M}{(\mathrm{tr}(B_0)+M\Lambda {)}^M}.$$

$$f(x_{t+1}^{s+1})\le f(x_t^{s+1})+g(\nabla f(x_t^{s+1}),{\text{Exp}}_{x_t^{s+1}}^{-1}(x_{t+1}^{s+1}))+\frac{L}{2}\Vert {\text{Exp}}_{x_t^{s+1}}^{-1}(x_{t+1}^{s+1}){\Vert }^2.$$

$$\mathbb{E}\left[f(x^{T+1})-f(x^{\ast })\right]\le L{S}^{-1}{\beta }^T\mathbb{E}\left[f(x^0)-f(x^{\ast })\right],$$

$$f(x_{i+1}^{t+1})$$

$$=f(x_i^{t+1})-{\eta }_2\cdot g\left(\nabla f(x_i^{t+1}),H_r^{t+1}{\nu }_i^{t+1}\right)+\frac{L{\eta }_2^2}{2}\Vert H_r^{t+1}{\nu }_i^{t+1}{\Vert }^2,$$

$$\mathbb{E}f(x_{i+1}^{t+1})\le \mathbb{E}f(x_i^{t+1})-{\eta }_2\gamma \Vert \nabla f(x_i^{t+1}){\Vert }^2+{\eta }_2^2{L}^3{\Gamma }^2\left[2\Vert {\text{Exp}}_{x_i^{t+1}}^{-1}(x^{\ast }){\Vert }^2+3\Vert {\text{Exp}}_{x^t}^{-1}(x^{\ast }){\Vert }^2\right],$$

$$\mathbb{E}\Vert {\nu }_i^{t+1}{\Vert }^2\le 4L^2\Vert {\text{Exp}}_{x_i^{t+1}}^{-1}(x^{\ast }){\Vert }^2+6L^2\Vert {\text{Exp}}_{x^t}^{-1}(x^{\ast }){\Vert }^2.$$

$$\Vert \nabla f(x_i^{t+1}){\Vert }^2\ge 2\kappa \left(f(x_i^{t+1})-f(x^{\ast })\right)\ge S\kappa \Vert {\text{Exp}}_{x_i^{t+1}}^{-1}(x^{\ast }){\Vert }^2,$$

Full text

1 Introduction

Optimization algorithms are a mainstay in machine learning research, underpinning solvers for a wide swath of problems ranging from linear regression and SVMs to deep learning. Consequently, scaling such algorithms to large scale datasets while preserving theoretical guarantees is of paramount importance. An important challenge in this field is designing scalable algorithms for optimization problems in the presence of constraints on the search space, a situation all too often encountered in real life. One approach to handling such constrained optimization problems on vector spaces is to reformulate them as optimization tasks on a suitable Riemannian manifold, with the constraints acting as manifold parametrization. Often, the problems can be shown to possess desirable geometric properties like convexity with respect to distance-minimizing geodesics on the manifold, leading to provably efficient optimization algorithms [1, 2, 3, 4]. These ideas can then be combined with stochastic optimization techniques influenced by [5], to deal with large datasets with theoretical convergence guarantees. See [6, 7] for recent examples. For instance, we can consider the problem of computing leading eigenvectors in the PCA setting [8] with unit-norm constraints. Projection-based strategies are normally used for this kind of problems [9], but alternating between solving and projecting can be prohibitively expensive in high dimensions. However, the unit-norm constraint can be used to cast the eigenvector problem into an unconstrained optimization scenario on the unit sphere, which happens to be one of the most well-behaved Riemannian manifolds.

Once the problems have been cast onto manifolds, one would want fast optimization algorithms that potentially use stochastic minibatches to deal with very large datasets. Such algorithms operating in Euclidean space have been widely researched in the optimization literature, but their development for Riemannian manifolds has been limited so far. In particular, one should note that the convergence speed limitations of unconstrained stochastic algorithms in the Euclidean case apply to manifold optimization as well; for instance a straightforward port of stochastic gradient descent to Riemannian manifolds [2] attains the same sublinear convergence seen in Euclidean space. There has been extensive work in the Euclidean domain using variance-reduced gradients to address this issue, with the aim of improving convergence rates by explicitly reducing the variance of stochastic gradients with suitably spaced full-gradient evaluations [8, 10]. Another nice advantage of this technique is the removal of the need for decaying learning rates for proving convergence, thereby solving the sublinearity issue as well sidestepping the nontrivial task of selecting an appropriate decay rate for SGD-like algorithms in large-scale optimization scenarios. Researchers have begun porting these methods to the manifold optimization domain, with a stochastic first-order variance reduced technique [7] showing robust convergence guarantees for both convex and nonconvex problems on geodesically complete manifolds.

Another complementary approach to improving convergence rates is of course using second-order updates for the iterates. In the Euclidean setting, one can show quadratic convergence rates for convex problems using Newton iterations, but these tend to be prohibitively expensive in high-dimensional big-data settings due to the need to store and invert the Hessian matrix. This limitation has led to the development of quasi-Newton methods, most notably L-BFGS [11], which uses lower-order terms to approximate the inverse Hessian. The curvature information provided by the Hessian estimate allows superlinear convergence in ideal settings [12]. While widely used for small-to-medium scale problems, adoption of these methods for big data problems has been limited, since the second order updates can be prohibitively expensive to compute in these situations. However, most optimization algorithms in the literature that use stochastic minibatching techniques to deal with large datasets are modifications of first order gradient-descent [13, 14] with relatively slower convergence in practical situations. This has recently begun to be addressed, with researchers devising stochastic variants of the L-BFGS technique [15], with straightforward convergence analyses. This has also been combined with variance reduction techniques and shown to have a linear convergence rate for convex problems in Euclidean space [16]. Our work in this paper is in a similar vein: we study quasi-Newton L-BFGS updates with stochastic variance reduction techniques for optimization problems on Riemannian manifolds, and analyze their convergence behavior for convex and nonconvex functions.

Contributions: The main contributions of this work may be summarized as follows:

1. We propose a stochastic L-BFGS method for Riemannian manifolds using stochastic variance reduction techniques for the first-order gradient estimates, and analyze the convergence for both convex and nonconvex functions under standard assumptions.

2. Our proof for strongly convex functions is different from those of recently proposed stochastic L-BFGS algorithms using variance-reduced gradients in Euclidean space [16] due to different bounds on the stochastic gradients. We do not use sectional curvature bounds in our proof for the convex case, making it structurally different from that of Riemannian SVRG [7].

3. We show strong experimental results on Karcher mean computations and calculation of leading eigenvalues, with noticeably better performance than Riemannian SVRG and VR-PCA; the latter is one of the best performing Euclidean algorithms for the finding dominant eigenvalues, that also uses stochastic variance-reduced gradients.

2 Preliminaries

2.1 Riemannian geometry

We begin with a brief overview of the differential geometric concepts we use in this work. We consider $C^{\infty}$ (smooth) manifolds that are locally homeomorphic to open subsets of $\mathbb{R}^{D}$, in the sense that the neighborhood of each point can be assigned a system of coordinates of appropriate dimensionality. Formally, this is defined with the notion of a chart $c:U\rightarrow\mathbb{R}^{D}$ at each $x\in\mathcal{M}$, where $U\subset\mathcal{M}$ is an open subspace containing $x$. Smooth manifolds are ones with covering collections of differentiable ($C^{\infty}$) charts. A Riemannian metric $g(\cdot,\cdot)$ is a bilinear $C^{\infty}$ tensor field of type $0\choose 2$, that is also symmetric and positive definite. A manifold endowed with such a metric is called a Riemannian manifold. The tangent space $T_{x}\mathcal{M}$ at every $x\in\mathcal{M}$ is a vector space, with the Riemannian metric $g:T_{x}\mathcal{M}\times T_{x}\mathcal{M}\rightarrow\mathbb{R}$ as the attendant metric. $g$ then induces a norm for vectors in the tangent space , which we denote by $\|\cdot\|$.

Riemannian manifolds are endowed with the Levi-Civita connection, which induces the notion of parallel transport of vectors from one tangent space to another along a geodesic, in a metric preserving way. That is, we have an operator $\Gamma_{\gamma}:T_{x}\mathcal{M}\rightarrow T_{y}\mathcal{M}$ where, informally speaking, $\gamma$ joins $x\text{ and }y$, and for any $u,\nu\in T_{x}\mathcal{M}$, we have $g(u,\nu)=g(\Gamma(u),\Gamma(\nu))$ . The parallel transport can be shown to be an isometry.

For every smooth curve $\gamma:[0,1]\rightarrow\mathcal{M}$ lying in $\mathcal{M}$, we denote its velocity vector as $\dot{\gamma}(t)\in T_{x}\mathcal{M}$ for each $t\in[0,1]$, with the “speed” given by $\|\dot{\gamma}(t)\|$. The length of such a curve is usually measured as $L(\gamma)=\int\limits_{0}^{1}\|\dot{\gamma}(t)\|dt.$ Denoting the covariant derivative along $\gamma$ of some $\nu\in T_{x}\mathcal{M}$, with respect to the Riemannian (Levi-Civita) connection by $A\nu$, we call $A\dot{\gamma}$ the acceleration of the curve. Curves with constant velocities ($A\dot{\gamma}\equiv 0$) are called geodesics, and can be shown to generalize the notion of Euclidean straight lines. We assume that every pair $x,y\in\mathcal{M}$ can be connected by a geodesic $\gamma$ s.t. $\gamma(0)=x\text{ and }\gamma(1)=y$. Immediately we have the notion of “distance” between any $x,y\in\mathcal{M}$ as the minimum length of all geodesics connecting $x\text{ and }y$, assuming the manifolds are geodesically complete as mentioned above, in that every countable decreasing sequence of lengths of geodesics connecting a pair of points has a well-defined limit.

The geodesic induces a useful operator called the exponential map, defined as $\text{Exp}_{x}:T_{x}\mathcal{M}\rightarrow\mathcal{M}\text{ s.t. }\text{Exp}_{x}(\nu)=\gamma(1)\text{ where }\gamma(0)=x,\gamma(1)=y\text{ and }\dot{\gamma}(0)=\nu$. If there is a unique geodesic connecting $x\text{ and }y$, then the exponential map has an inverse, denoted by $\text{Exp}^{-1}_{x}(y)$. The length of this geodesic can therefore be seen to be $\|\text{Exp}^{-1}_{x}(y)\|$.

The derivative $D$ of a differentiable function is defined using the Riemannian connection by the following equivalence: $Df(x)\nu=\nu f$, where $\nu\in T_{x}\mathcal{M}$. Then, by the Riesz representation theorem, there exists a gradient $\nabla f(x)\in T_{x}\mathcal{M}$ s.t. $\forall\nu\in T_{x}\mathcal{M},\>Df(x)\nu=g_{x}(\nabla f(x),\nu)$. Similarly, we can denote the Hessian as follows: $D^{2}f(x)(\cdot,\cdot):T_{x}\mathcal{M}\times T_{x}\mathcal{M}\rightarrow\mathbb{R}$. We denote the mapping from $\nu\in T_{x}\mathcal{M}$ to the Riesz representation of $D^{2}f(x)(\nu,\cdot)$ by $\nabla^{2}f(x)$. One can consult standard textbooks on differential geometry [17, 18] for more details.

2.2 Convexity and Lipschitz smoothness on manifolds

Similar to [2, 3, 7], we define manifold (or geodesic) convexity concepts analogous to the Euclidean baselines, as follows : a set $U\subset\mathcal{M}$ is convex on the manifold if $\forall x,y\in U$ there exists a geodesic $\gamma$ connecting $x,y$ that completely lies in $U$, i.e. $\gamma(0)=x,\gamma(1)=y$. Then, a function can be defined as convex w.r.t. geodesics if $\forall x,y\in U\text{ where }\exists\gamma$ connecting $x,y$ on the manifold, we have:

[TABLE]

We can also define a notion of strong convexity as follows: a function $f$ is called $S-$strongly convex if for any $x,y\in U$ and (sub)gradient $\nabla_{x}$, we have

[TABLE]

We define Lipschitz smoothness of a function $f$ by imposing Lipschitz continuity on the gradients, as follows: $\forall x,y\in U,$

[TABLE]

where L is the smoothness parameter. Analogous to the Euclidean case, this property can also be formulated as::

[TABLE]

3 Stochastic Riemannian L-BFGS

In this section we present our stochastic variance-reduced L-BFGS algorithm on Riemannian manifolds and analyze the convergence behavior for convex and nonconvex differentiable functions on Riemannian manifolds. We assume these manifolds to be $L$-Lipschitz smooth, as defined above, with existence of unique distance-minimizing geodesics between every two points, i.e. our manifolds are geodesically complete; this allows us to have a well-defined inverse exponential map that encodes the distance between a pair of points on the manifold. For the convergence analysis, we also assume $f$ to have a unique minimum at $x^{*}\in U$, where $U$ is a compact convex subset of the manifold.

3.1 The Algorithm

The pseudocode is shown in Algorithm 1. We provide a brief discussion of the salient properties, and compare it to similar algorithms in the Euclidean domain, for example [15, 16], as well as those on Riemannian manifolds, for example [3]. To begin, note that $\nabla$ denotes the Riesz representation of the gradient $D$, as defined in $\S$2.1. We denote full gradients by $\nabla$ and stochastic gradients by $\tilde{\nabla}$. Similar to other stochastic algorithms with variance-reduction, we use two loops: each iteration of the inner loop corresponds to drawing one minibatch from the data and performing the stochastic gradient computations (Steps $10$, $11$), whereas each outer loop iteration corresponds to two passes over the full dataset, once to compute the full gradient (Step $6$) and the other to make multiple minibatch runs (Steps $8$ through $30$). Compared to the Euclidean setting, note that the computation of the variance-reduced gradient in Step $11$ involves an extra step: the gradients ($\nabla f(x)$-s) reside in the tangent spaces of the respective iterates, therefore we have to perform parallel transport to bring them to the same tangent space before performing linear combinations.

To avoid memory issues and complications arising from Hessian-vector computations in the Riemann setting, we chose to update the second correction variable $y_{r}$ using a simple difference of gradients approximation: $y_{r}=\tilde{\nabla}f(u_{r})-\Gamma_{\gamma}\tilde{\nabla}f(u_{r-1})$. We should note here that the parallel transport parametrization should be clear from the context; $\Gamma_{\gamma}$ here denotes transporting the vector $\tilde{\nabla}f(u_{r-1})\in T_{u_{r-1}}\mathcal{M}$ to $T_{u_{r}}\mathcal{M}$ along the connecting geodesic $\gamma$. We omit any relevant annotations from the transport symbol to prevent notational overload. We calculate the first correction pair $z_{r}$ in one of two ways: (a) as $z_{r}=\Gamma_{\gamma}\left(\eta_{2}\rho_{t-1}\right)$, or (b) as $z_{r}=\Gamma_{\gamma}\left(-\eta_{1}\nu_{\text{prev}}\right)$. We denote these by Option $\mathbf{1}$ and Option $\mathbf{2}$ respectively in Alg. 1. Note that in both cases, $\Gamma_{\gamma}$ denotes the parallel transport of the argument to the tangent space at $x_{t}^{s+1}$. In our experiments, we noticed faster convergence for the strongly convex centroid computation problem with Option $\mathbf{1}$, along with computation of the correction pairs every iteration and a low memory pool. For calculating dominating eigenvalues on the unit-radius sphere, Option $\mathbf{2}$ yielded better results. Once the corrections pairs $z_{r},y_{r}$ have been computed, we compute the descent step in Step 21 using the standard two-loop recursion formula given in [12], using the $M$ correction pairs stored in memory. Note that we use fixed stepsize in the update steps and in computing the correction pairs, and do not impose or perform calculations designed to have them satisfy Armijo or Wolfe conditions.

Compared to the Euclidean algorithms [15, 16], Alg.1 has some key differences: 1) we did not notice any significant advantage from using separate minibatches in Steps $10$ and $18$, therefore we use the same minibatch to compute the VR gradient and the correction elements $y_{r}$; 2) we do not keep a running average of the iterates for computing the correction element $z_{r}$ (Steps $15$ through $17$); 3) we use constant stepsizes throughout the whole process, in contrast to [15] that uses a decaying sequence. Note that, as seen in Step $24$, we use the first-order VR gradient to update the iterates for the first $2R$ iterations; this is because we calculate correction pairs every $R$ steps and evaluate the gradient-inverse Hessian product (Step $26$) once at least two pairs have been collected. Similar to [12], we drop the oldest pair to maintain the memory depth $M$. Compared to the algorithms in [3, 19], ours uses stochastic VR gradients, with all the attendant modifications and advantages, and does not use linesearch techniques to satisfy Wolfe conditions.

3.2 Analysis of convergence

In this section we provide the main convergence results of the algorithm. We analyze convergence for finite-sum empirical risk minimization problems of the following form:

[TABLE]

where the Riemannian manifold is denoted by $\mathcal{M}$. Note that the iterates are updated in Algorithm 1 by taking the exponential map of the descent step multiplied by the stepsize, with the descent step computed as the product of the inverse Hessian estimate and the stochastic variance-reduced gradient using the standard two-loop recursion formula. Thus, to bound the optimization error using the iterates, we will need bounds on both the stochastic gradients and the inverse Hessians. As mentioned in [7], the methods used to derive the former bounds for Euclidean algorithms cannot be ported directly to manifolds due to metric nonlinearities; see the proof of Proposition 1 for details. For the latter, we follow the standard template for L-BFGS algorithms in the literature [4, 12, 15]. To begin, we make the following assumptions:

Assumption 1. The function $f$ in (5) is strongly convex on the manifold, whereas the $f_{i}$s are individually convex.

Assumption 2. There exist $\lambda,\Lambda\in(0,\infty),\>\lambda<\Lambda$ s.t. $\lambda\|\nu\|_{x}^{2}\leq D^{2}f\leq\Lambda\|\nu\|_{x}^{2}\quad\forall\nu\in T_{x}\mathcal{M}$.

These two assumptions allow us to (a) guarantee that $f$ has a unique minimizer $x^{*}$ in the convex sublevel set $U$, and (b) derive bounds on the inverse Hessian updates using BFGS update formulae for the Hessian approximations. Similar to the Euclidean case, these can be written as follows:

[TABLE]

and by the Sherman-Morrison-Woodbury lemma, that of the inverse:

[TABLE]

where $G=I-\frac{g_{x_{r-1}}(s_{r-1},\cdot)\hat{y}_{r-1}}{y_{r-1}s_{r-1}}$, and $\hat{B}_{r}$ is the Lax-Milgram representation of the Hessian. Details on these constructs can be found in [4], in addition to [17, 18].

3.2.1 Trace and determinant bounds

To start off our convergence discussions for both convex and nonconvex cases, we derive bounds for the trace and determinants of the Hessian approximations, followed by those for their inverses. The techniques used to do so are straightforward ports of the Euclidean originals [12], with some minor modifications to account for differential geometric technicalities. Using the assumptions above, we can prove the following bounds [4]:

Lemma 1**.**

Let $B_{r}^{s+1}=\left(H_{r}^{s+1}\right)^{-1}$ be the approximation of the Hessian generated by Algorithm 1, and $\hat{B}_{r}^{s+1}$ and $\hat{H}_{r}^{s+1}$ be the corresponding Lax-Milgram representations. Let $M$, the memory parameter, be the number of correction pairs used to update the inverse Hessian approximation. Then, under Assumptions 1 and 2, we have:

[TABLE]

Also, $\gamma I\preceq\hat{H}_{r}^{s+1}\preceq\Gamma I$, for some $\Gamma\geq\gamma>0$.

From a notational perspective, recall that our notation for the parallel transport operator is $\Gamma_{\gamma}$, with the subscript denoting the geodesic. The symbols $\gamma$ and $\Gamma$ in Lemma 1 above are unrelated to these geometric concepts, merely being the derived bounds on the eigenvalues of inverse Hessian approximations. The proof is given in the supplementary for completeness.

3.2.2 Convergence result for strongly convex functions

Our convergence result for strongly convex functions on the manifold can be stated as follows:

Proposition 1**.**

Let the Assumptions 1 and 2 hold. Further, let the $f(\cdot)$ in (5) be S-strongly convex, and each of the $f_{i}$ be L-smooth, as defined earlier. Define the following constants: $p=\left[LS^{-1}+2\eta_{2}S^{-1}\left\{2\eta L^{3}\Gamma^{2}-S\kappa\gamma\right\}\right]$, and $q^{\prime}=6\eta^{2}L^{3}\Gamma^{2}S^{-1}$. Denote the global optimum by $x^{*}$. Then the iterate $x^{T+1}$ obtained after $T$ outer loop iterations will satisfy the following condition:

[TABLE]

where the constants are chosen to satisfy $\beta=\left(1-p\right)^{-1}\left(q^{\prime}+p^{T}(1-p-q^{\prime})\right)<1$ for linear convergence.

For proving this statement, we will use the $L$-smoothness (2) and $S$-strong convexity (1) conditions mentioned earlier. As in the Euclidean case [10, 16], we will also require a bound on the stochastic variance-reduced gradients. These can be bounded using triangle inequalities and $L$-smoothness on Riemannian manifolds, as shown in [7]. This alternative is necessary since the Euclidean bound first derived in [10], using the difference of the objective function at the iterates, cannot be ported directly to manifolds due to metric nonlinearities. Thus we take a different approach in our proof compared to the Euclidean case of [16], using the interpoint distances defined with the norms of inverse exponential maps. We do not use trigonometric distance inequalities [3, 6] for the convex case either, making the overall structure different from the proof of Riemannian SVRG as well. The details are deferred to the supplementary due to space limitations. However we do use the trigonometric inequality along with assumed lower bounds on sectional curvature for showing convergence for nonconvex functions, as described next.

3.2.3 Convergence for the nonconvex case

Here we provide a convergence result for nonconvex functions satisfying the following condition: $f(x^{t})-f(x^{*})\leq\kappa^{-1}\|\nabla f(x^{t})\|^{2}$, which automatically holds for strongly convex functions. We assume this to hold even if $f$ is nonconvex, since it allows us to show convergence of the iterates using $\|\nabla f(x^{t})\|^{2}$. Further, similar to [2, 7] we assume that the sectional curvature of the manifold is lower bounded by $c_{\delta}$. This allows us to derive a trigonometric inequality analogous to the Euclidean case, where the sides of the “triangle” are geodesics [6]. The details are given in the supplementary. Additionally, we assume that the eigenvalues of the inverse Hessian are bounded by $(\gamma,\Gamma)$ within some suitable region around an optimum. The main result of this section may be stated as follows:

Proposition 2**.**

Let the sectional curvature of the manifold be bounded below by $c_{\delta}$, and the $f_{i}$ be $L$-smooth. Let $x^{*}$ be an optimum of $f(\cdot)$ in (5). Assume the eigenvalues of the inverse Hessian estimates are bounded. Set $\eta_{2}=\mu_{0}/\left(\Gamma Ln^{\alpha_{1}}\eta^{\alpha_{2}}\right)$, $K=mT$, and $m=\lfloor n^{\nicefrac{{3\alpha_{1}}}{{2}}}/\left(3\mu_{0}\zeta^{1-2\alpha_{2}}\right)\rfloor$, where $\alpha_{1}\in(0,1]$ and $\alpha_{2}\in[0,2]$. Then, for suitable choices of the inverse Hessian bounds $\gamma,\Gamma$, we can find values for the constants $\mu_{0}>0$ and $\epsilon>0$ so that the following holds:

[TABLE]

$\zeta$ is defined as $\zeta=\left(\tanh\left(d\sqrt{|c_{\delta}|}\right)\right)^{-1}d\sqrt{|c_{\delta}|}$ if $c_{\delta}<0$, and [math] otherwise; $d$ is an upper bound on the diameter of the set $U$ mentioned earlier, containing an optimum $x^{*}$. The proof is inspired by similar results from both Euclidean [20] and Riemannian [7] analyses, and is given in the supplementary. One way to deal with negative curvature in Hessians in Euclidean space is by adding some suitable positive $\alpha$ to the diagonal, ensuring bounds on the eigenvalues. Investigation of such “damping” methods in the Riemannian context could be an interesting area of future work.

4 Experiments

4.1 Karcher mean computation for PD matrices

We begin with a synthetic experiment on learning the Karcher mean (centroid) [21] of positive definite matrices. For a collection of matrices $\left\{\mathbf{W}_{i}\right\}_{i=1}^{N}$, the optimization problem can be stated as follows:

[TABLE]

We compare our minibatched implementation of the Riemannian SVRG algorithm from [7], denoted as rSVRG, with the stochastic variance-reduced L-BFGS procedure from Algorithm 1, denoted rSV-LBFGS. We implemented both algorithms using the Manopt [22] and Mixest [23] toolkits. We generated three sets of random positive definite matrices, each of size $100\times 100$, with condition numbers $10,1e2,$ and $1e3$, and computed the ground truths using code from [24]. Matrix counts were $100$ for condition number $1e-2$, and $1000$ for the rest. Both algorithms used equal batchsizes of $50$ for the first and third datasets, and $5$ for the second, and were initialized identically. Both used learning rates satisfying their convergence conditions. In general we found rSV-LBFGS to perform better with frequent correction pair calculations and a low retention rate, ostensibly due to the strong convexity; therefore we used $R=1$, $M=2$ for all three datasets. As mentioned earlier, the $z_{r}$ correction pair was calculated using Option 1: $z_{r}=\Gamma_{\gamma}\left(\eta_{2}\rho_{t}\right)$ where $\rho_{t}$ is calculated using the two-loop recursion. We used standard retractions to approximate the exponential maps. The retraction formulae for both symmetric PD and sphere manifolds used in the next section are given in the supplementary.

We calculated the error of iterate $\mathbf{W}$ as $\|\mathbf{W}-\mathbf{W}^{*}\|_{F}^{2}$, with $\mathbf{W}^{*}$ being the ground truth. The log errors are plotted vs the number of data passes in Fig.1. Comparing convergence speed in terms of # data passes is often the preferred approach for benchmarking ML algorithms since it is an implementation-agnostic evaluation and focuses on the key bottleneck (I/O) for big data problems. Comparisons of rSVRG with Riemannian gradient descent methods, both batch and stochastic, can be found in [7]. From Fig.1, we find rSV-LBFGS to converge faster than rSVRG for all three datasets.

4.2 Leading eigenvalue computation

Next we conduct a synthetic experiment on calculating the leading eigenvalue of matrices. This of course is a common problem in machine learning, and as such is a unit-norm constrained nonconvex optimization problem in Euclidean space. It can be written as:

[TABLE]

where $D^{d\times N}$ is the data matrix, and $d_{i}$ are its columns. We can transform this into an unconstrained manifold optimization problem on the sphere defined by the norm constraint. To that end, we generated four sets of datapoints, for eigengaps $0.005$, $0.05$, $0.01$ and $0.1$, using the techniques described in [8]. Each dataset contains $100,000$ vectors of dimension $1000$. We used a minibatch size of $100$ for the two Riemannian algorithms. As before, learning rates for rSVRG were chosen according to the bounds in [7]. Selecting appropriate values for the four parameters in rSV-LBFGS (the first and second-order learning rates, $L$ and $M$) was a nontrivial task; after careful grid searches within the bounds defined by the convergence conditions, we chose $\eta_{1}=0.001,\eta_{2}=0.1,$ and $M=10$ for all four datasets. $L$ was set to $5$ for the dataset with eigengap $0.005$, and $10$ for the rest. The $z_{r}$ correction pair was calculated using Option 2: $z_{r}=\Gamma_{\gamma}\left(-\eta_{1}\nu_{\text{prev}}\right)$. We plot the performance of rSV-LBFGS, rSVRG and VR-PCA in Fig.2. Extensive comparisons of VR-PCA with other Euclidean algorithms have been conducted in [8]; we do not repeat them here. We computed the error of iterate $\mathbf{z}$ as $1-\left(Ne^{*}\right)^{-1}\|D^{T}\mathbf{z}\|_{2}^{2}$, $e^{*}$ being the ground truth obtained from Matlab’s $eigs$.

We see that the rSV-LBFGS method performs well on all four datasets, reaching errors of the order of $1e-15$ well before VR-PCA and rSVRG in the last three cases. The performance delta relative to VR-PCA is particularly noticeable in each of the four cases; we consider this to be a noteworthy result for fixed-stepsize algorithms on Riemannian manifolds.

5 Conclusion

We propose a novel L-BFGS algorithm on Riemannian manifolds with variance reduced stochastic gradients, and provide theoretical analyses for strongly convex functions on manifolds. We conduct experiments on computing Riemannian centroids for symmetric positive definite matrices, and calculation of leading eigenvalues, both using large scale datasets. Our algorithm outperforms other Riemannian optimization algorithms with fixed stepsizes in both cases, and performs noticeably better than one of the fastest stochastic algorithms in Euclidean space, VR-PCA, for the latter case.

6 Appendices

We present the convergence results from Propositions 1 and 2 in the main text in this section.

6.1 Analysis of convergence

We analyze convergence for finite-sum empirical risk minimization problems of the following form:

[TABLE]

where the Riemannian manifold is denoted by $\mathcal{M}$. Note that the iterates are updated in Algorithm 1 by taking the exponential map of the descent step multiplied by the stepsize, with the descent step computed as the product of the inverse Hessian estimate and the stochastic variance-reduced gradient using the standard two-loop recursion formula. Thus, to bound the optimization error using the iterates, we will need bounds on both the stochastic gradients and the inverse Hessians. As mentioned in [7], the methods used to derive the former bounds for Euclidean algorithms cannot be ported directly to manifolds due to metric nonlinearities; see the proof of Proposition 1 for details. For the latter, we follow the standard template for L-BFGS algorithms in the literature [4, 12, 15]. To begin, we make the following assumptions:

Assumption 1. The function $f$ in (5) is strongly convex on the manifold, whereas the $f_{i}$s are individually convex.

Assumption 2. There exist $\lambda,\Lambda\in(0,\infty),\>\lambda<\Lambda$ s.t. $\lambda\|\nu\|_{x}^{2}\leq D^{2}f\leq\Lambda\|\nu\|_{x}^{2}\quad\forall\nu\in T_{x}\mathcal{M}$.

These two assumptions allow us to (a) guarantee that $f$ has a unique minimizer $x^{*}$ in the convex sublevel set $U$, and (b) derive bounds on the inverse Hessian updates using BFGS update formulae for the Hessian approximations. Similar to the Euclidean case, these can be written as follows:

[TABLE]

and by the Sherman-Morrison-Woodbury lemma, that of the inverse:

[TABLE]

where $G=I-\frac{g_{x_{r-1}}(s_{r-1},\cdot)\hat{y}_{r-1}}{y_{r-1}s_{r-1}}$. The $\hat{B}_{r}$ is the Lax-Milgram representation of the Hessian [4].

6.1.1 Trace and determinant bounds

To start off our convergence discussions for both convex and nonconvex cases, we derive bounds for the trace and determinants of the Hessian approximations, followed by those for their inverses. The techniques used to do so are straightforward ports of the Euclidean originals [12], with some minor modifications to account for differential geometric technicalities. Using the assumptions above, we can prove the following bounds [4]:

Lemma 1**.**

Let $B_{r}^{s+1}=\left(H_{r}^{s+1}\right)^{-1}$ be the approximation of the Hessian generated by Algorithm 1, and $\hat{B}_{r}^{s+1}$ and $\hat{H}_{r}^{s+1}$ be the corresponding Lax-Milgram representations. Let $M$, the memory parameter, be the number of correction pairs used to update the inverse Hessian approximation. Then, under Assumptions 1 and 2, we have:

[TABLE]

Also, $\gamma I\preceq\hat{H}_{r}^{s+1}\preceq\Gamma I$, for some $\Gamma\geq\gamma>0$.

Proof.

For brevity of notation we temporarily drop the $(s+1)$ superscript. The proof for the Euclidean case [15, 25] can be generalized to the Riemannian scenario in a straightforward way, as follows. Define the average Hessian $G_{r}$ by

[TABLE]

such that $y_{r}=G_{r}(z_{r},\cdot)$. Then, it can be easily shown that $G_{r}$ satisfies the bounds in Assumption 2. Furthermore, we have the following useful inequalities

[TABLE]

Let $\hat{y}_{r}$ be the Riesz representation of $y_{r}$. Recall that parallel transport is an isometry along the unique geodesics, which implies invariance of the trace operator. Then using the L-BFGS update (6) and (7), we can bound the trace of the Lax-Milgram representation of the Hessian approximations as follows:

[TABLE]

This therefore proves boundedness of the largest eigenvalue of the $\hat{B}_{r}$ estimates.

Similarly, to get a lower bound for the minimum eigenvalue, we bound the determinant as follows:

[TABLE]

where we use $\lambda_{\text{max}}$ to denote the maximum eigenvalue of $B_{r-1}$, and use (7). Since $\lambda_{\text{max}}$ is bounded above by the trace of $(\hat{B}_{r})$, we can telescope the inequality above to get

[TABLE]

The bounds on the maximum and minimum eigenvalues of $B_{r}$ thus derived allows us to infer corresponding bounds for those of $H_{r}$ as well, since by definition $H_{r}=\hat{B}_{r}^{-1}$. ∎

6.1.2 Convergence results for the strongly convex case

Next we provide a brief overview of the bounds necessary to prove our convergence result. First, note the following bound implied by the Lipschitz continuity of the gradients:

[TABLE]

Note the update step fom line 11 of Algorithm 1: $x_{t+1}^{s+1}=\text{Exp}_{x_{t}^{s+1}}(-\eta H_{r}^{s+1}\nu_{t}^{s+1})$. We can replace the inverse exponential map in the inner product above by the quantity in the parentheses. In order to replace $H_{r}^{s+1}$ by the eigen-bounds from Lemma 1, we invoke the following result (Lemma 5.8 from [17]:

Lemma 2**.**

For any $D\in T_{x}\mathcal{M}$ and $c,t\in\mathbb{R}$, $\gamma_{cD}(t)=\gamma_{D}(ct)$,

where $\nu$ is the “speed” of the geodesic. This allows us to write $\text{Exp}_{x}(c\nu)=\gamma_{\nu}(c)=\gamma_{c\nu}(1)$. Recall that for Riemann geodesics we have $\|\dot{\gamma(t)}\|=\bar{s}\text{ for all }t\in[0,1]$, a constant.

Proposition 1**.**

Let the Assumptions 1 and 2 hold. Further, let the $f(\cdot)$ in (5) be S-strongly convex, and each of the $f_{i}$ be L-smooth, as defined earlier. Define the following constants: $p=\left[LS^{-1}+2\eta_{2}S^{-1}\left\{2\eta L^{3}\Gamma^{2}-S\kappa\gamma\right\}\right]$, and $q^{\prime}=6\eta^{2}L^{3}\Gamma^{2}S^{-1}$. Denote the global optimum by $x^{*}$. Then the iterate $x^{T+1}$ obtained after $T$ outer loop iterations will satisfy the following condition:

[TABLE]

where the constants are chosen to satisfy $\beta=\left(1-p\right)^{-1}\left(q^{\prime}+p^{T}(1-p-q^{\prime})\right)<1$ for linear convergence.

Proof.

From the $L$-smoothness condition (2), we have the following:

[TABLE]

where we have omitted subscripts from the metric. Taking expectations, and using the bounds on the inverse Hessian estimates derived in Lemma 1, we have the following:

[TABLE]

where we have used the following bound on the stochastic variance-reduced gradients derived in [7]:

[TABLE]

This can be derived using triangle inequalities and the $L$-smoothness assumption. Note that the bound is different from the Euclidean case [10], due to technicalities introduced by the Riemannian metric not being linear in general.

Now, recall the condition $f(x^{t})-f(x^{*})\leq 2\kappa)^{-1}\|\nabla f(x^{t})\|^{2}$, which follows from strong convexity. Using this, we can derive a bound on the $\nabla$ term in (8) as follows:

[TABLE]

where the second inequality follows from $S$-strong convexity (1), since $\nabla f(x^{*})=0$. Plugging this into (8), we have the following:

[TABLE]

Now, note that $S$-strong convexity allows us to write the following:

[TABLE]

where the last step follows from $L$-smoothness. Taking expectations of both sides, and using (9) for the first component on the right, we have

[TABLE]

Now, we denote $p=\left[\frac{L}{S}+\frac{2\eta_{2}}{S}\left\{2\eta L^{3}\Gamma^{2}-S\kappa\gamma\right\}\right]$, and $q^{\prime}=\frac{6\eta_{2}^{2}L^{3}\Gamma^{2}}{S}$. Then, taking expectations over the sigma algebra of all the random variables till minibatch $m$, and some algebra, it can be shown that:

[TABLE]

where $q=(1-p)^{-1}q^{\prime}$. Note that this provides a bound on the iterate at the end of the inner minibatch loop. Telescoping further, we have the bound

[TABLE]

where $\beta=\frac{q^{\prime}+p^{T}(1-p-q^{\prime})}{1-p}$. Then, using this result with a final appeal to the $L$-Lipschitz and $S$-strong convexity conditions, we have the bounds

[TABLE]

thereby completing the proof. ∎

6.1.3 Convergence results for the nonconvex case

We begin with the following inequality involving the side lengths of a geodesic “triangle” [3, 6]:

Lemma 3**.**

Let the sectional curvature of a Riemannian manifold be bounded below by $c_{\delta}$. Let $A$ be the angle between sides of length $b$ and $c$ in a triangle on the manifold, with the third side of length $a$, as usual. Then the following holds:

[TABLE]

The cosine is defined using inner products, as in the Euclidean case, and the distances using inverse exponential maps, as seen above. The following sequence of results and proofs are inspired by the basic structure of [20], with suitable modifications involving the inverse Hessian estimates from the L-BFGS updates.

Lemma 4**.**

Let the assumptions of Proposition 2 hold. Define the following functions:

[TABLE]

where $c_{i},c_{i+1},\beta,\eta_{2}>0$. Further, for $0\leq t\leq T-1$, define the Lyapunov function $R_{i}^{t+1}=\mathbb{E}\left[f(x_{i}^{t+1})+c_{i}\|\text{Exp}^{-1}_{x^{t}}(x_{i}^{t+1})\|^{2}\right]$. Then we have the following bound:

[TABLE]

Proof.

As with the proof of Proposition 1, we begin with the following bound derived from $L$-smoothness:

[TABLE]

where we have bound the bounds on the inverse Hessian derived earlier. Then using Lemma 3 above, we have:

[TABLE]

where we have used $g(a,b)\leq\frac{1}{2\beta}\|a\|^{2}+\frac{\beta}{2}\|b\|^{2}$. Note that we have used the norm of the inverse exponential maps as the side lengths in Lemma 3. Using these last two results, we can derive the following bound for the Lyapunov functions $R_{i+1}^{t+1}$:

[TABLE]

The norm of the stochastic variance reduced gradient can be bounded as follows [7, 20]:

[TABLE]

This allows us to bound the Lyapunov function above as:

[TABLE]

which completes the proof. ∎

Next we present a bound on $\|\nabla f(\cdot)\|^{2}$ using the $\delta_{i}$’s defined above (Thm 6 of [7]):

Lemma 5**.**

Let the conditions of Lemma 4 hold, and define the quantities therein. Let $\delta_{i}>0$ $\forall i\in[0,m]$, and $c_{m}=0$. Let $\delta_{\delta}=\min_{i}\delta_{i}$, and $K=mT$. Then if we randomly return one of the iterates $\left\{x_{i}^{t+1}\right\}_{i=1}^{m}$ as $x^{t+1}$, then:

[TABLE]

This result can be shown by telescoping the bound derived in the previous lemma for the Lyapunov functions, using $c_{m}=0$.

Proposition 2**.**

Let the sectional curvature of the manifold be bounded below by $c_{\delta}$, and the $f_{i}$ be $L$-smooth. Let $x^{*}$ be an optimum of $f(\cdot)$ in (5). Assume the eigenvalues of the inverse Hessian estimates are bounded. Set $\eta_{2}=\mu_{0}/\left(\Gamma Ln^{\alpha_{1}}\eta^{\alpha_{2}}\right)$, $K=mT$, and $m=\lfloor n^{\nicefrac{{3\alpha_{1}}}{{2}}}/\left(3\mu_{0}\zeta^{1-2\alpha_{2}}\right)\rfloor$, where $\alpha_{1}\in(0,1]$ and $\alpha_{2}\in[0,2]$. Then, for suitable choices of the inverse Hessian bounds $\gamma,\Gamma$, we can find values for the constants $\mu_{0}>0$ and $\epsilon>0$ so that the following holds:

[TABLE]

Proof.

We define $\beta=L\zeta^{1-\alpha_{2}}/\left(n^{\nicefrac{{\alpha_{1}}}{{2}}}\Gamma\right)$. Also, as mentioned in the proposition, $\eta_{2}=\mu_{0}/\left(\Gamma Ln^{\alpha_{1}}\eta^{\alpha_{2}}\right)$, with appropriate $\alpha_{1}$, $\alpha_{2}$. Note that we need a bound for $\delta_{\delta}$ to plug in the denominator of the bound in Lemma 5 above. This quantity can be lower bounded as follow:

[TABLE]

Now we need to bound $c_{0}$. To that end, telescoping the $c_{i+1}$ function defined in Lemma 4 above with $c_{m}=0$, and denoting $\theta=\eta_{2}\beta\Gamma+2\zeta\eta_{2}^{2}L^{2}\Gamma^{2}$, we get the following:

[TABLE]

Using the definitions of $\eta_{2}$ and $\beta$ above, we note that $\theta<1/m$, implying $c_{0}\leq\frac{L\mu_{0}}{\zeta n^{\nicefrac{{\alpha_{1}}}{{2}}}}(e-1)$. Plugging this in the bound above, we posit that $\delta_{\delta}$ can be bounded below as follows:

[TABLE]

for some sufficiently small $\epsilon$, and suitable choices of the inverse Hessian bounds $\gamma,\Gamma$ and the rest of the parameters. Using this bound in the denominator of the right hand side of Lemma 5 above completes the proof. ∎

6.2 Retractions

We approximated the exponential maps with retractions from the Manopt [22] toolbox. We used the following formulae: $\mathbb{R}_{x}(\eta\rho)=x\cos\|\eta\rho\|_{F}+\frac{\eta\rho}{\|\eta\rho\|_{F}}\sin\|\eta\rho\|_{F}$ for the sphere manifold, and $\mathbb{R}_{x}(\eta\rho)=x\cdot M_{x}(x\setminus\eta\rho)$ for the manifold of symmetric PD matrices, where $M_{x}$ denotes the matrix exponential, and $\setminus$ is matrix division. Here $x\in\mathcal{M}$ is some point on the manifold, $\rho\in T_{x}\mathcal{M}$ is some descent step evaluated at $x$, and $\eta$ is the stepsize.