Riemannian stochastic variance reduced gradient algorithm with retraction and vector transport

TL;DR

This paper introduces a Riemannian stochastic variance reduced gradient algorithm that effectively handles manifold optimization challenges, with proven convergence and superior performance on several geometric problems.

Contribution

It extends Euclidean SVRG to Riemannian manifolds using retraction and vector transport, providing convergence analysis and practical applications.

Findings

Outperforms Riemannian stochastic gradient descent in experiments

Provides global convergence and local rate analysis

Successfully applied to Riemannian centroid, PCA, and matrix completion

Abstract

In recent years, stochastic variance reduction algorithms have attracted considerable attention for minimizing the average of a large but finite number of loss functions. This paper proposes a novel Riemannian extension of the Euclidean stochastic variance reduced gradient (R-SVRG) algorithm to a manifold search space. The key challenges of averaging, adding, and subtracting multiple gradients are addressed with retraction and vector transport. For the proposed algorithm, we present a global convergence analysis with a decaying step size as well as a local convergence rate analysis with a fixed step size under some natural assumptions. In addition, the proposed algorithm is applied to the computation problem of the Riemannian centroid on the symmetric positive definite (SPD) manifold as well as the principal component analysis and low-rank matrix completion problems on the Grassmann manifold. The results show that the proposed algorithm outperforms the standard Riemannian stochastic gradient descent algorithm in each case.

Equations

$$\underset{w\in \mathcal{M}}{\min }f(w):=\frac{1}{N}\sum _{n=1}^N{f}_n(w),$$

$$w_{t+1}={\mathrm{Exp}}_{w_t}\left(s_t{\xi }_{w_t}\right),$$

$$\mathrm{D}f(w)[{\xi }_w]$$

$$w_{t+1}$$

$$v_t^s$$

$${\xi }_t^s$$

$${\mathbb{E}}_{i_t^s}\left[{\xi }_t^s\right]$$

$$\sum {\left({\alpha }_t^s\right)}^2<\infty \mathrm{and}\sum {\alpha }_t^s=\infty ,$$

$$\sum _{n=0}^{\infty }\left|\mathbb{E}\left[X_{n+1}-X_n\mid {\mathcal{F}}_n\right]\right|<\infty a.s.,\mathrm{and}{X}_n\mathrm{converges}a.s.$$

$$\Vert {\xi }_{t+1}^s{\Vert }_{w_t^s}$$

$${\Vert -{\alpha }_t^s{\xi }_{t+1}^s\Vert }_{w_t^s}\le {\alpha }_t^s{C}^{\prime }<I.$$

$$\frac{d^2}{d{\tau }^2}g(\tau ;w,\xi )\le 2k_1,\tau \in [0,1],\xi \in \mathbb{B}(0_w,I),w\in K,$$

$$g(\alpha ;w,\xi )=g(0;w,\xi )+g^{\prime }(0;w,\xi )\alpha +{\alpha }^2{\int }_0^1(1-\tau )g^{\prime \prime }(\alpha \tau ;w,\xi )d\tau$$

$$f\left(w_{t+1}^s\right)-f\left(w_t^s\right)=$$

$$=$$

$$\le$$

$${\mathcal{F}}_t^s=\left\{i_1^1,i_2^1,\dots ,i_{m_1}^1,i_1^2,i_2^2,\dots ,i_{m_2}^2,\dots ,i_1^{s-1},i_2^{s-1},\dots ,i_{m_{s-1}}^{s-1},i_1^s,i_2^s,\dots ,i_{t-1}^s\right\}.$$

$$\mathbb{E}\left[{\xi }_{t+1}^s\mid {\mathcal{F}}_{t+1}^s\right]$$

$$=$$

$$\mathbb{E}\left[{⟨{\xi }_{t+1}^s,\mathrm{grad}f\left(w_t^s\right)⟩}_{w_t^s}\mid {\mathcal{F}}_{t+1}^s\right]={\Vert \mathrm{grad}f\left(w_t^s\right)\Vert }_{w_t^s}^2,$$

$$\mathbb{E}\left[f\left(w_{t+1}^s\right)-f\left(w_t^s\right)\mid {\mathcal{F}}_{t+1}^s\right]$$

$$\sum _{s\ge s_0,t}{\alpha }_t^s{\Vert \mathrm{grad}f\left(w_t^s\right)\Vert }_{w_t^s}^2$$

$$\sum _{s\ge s_0,t}\mathbb{E}\left[\mathbb{E}{\left[f\left(w_{t+1}^s\right)-f\left(w_t^s\right)\mid {\mathcal{F}}_{t+1}^s\right]}^{+}\right]<\infty$$

$$\left|-\sum _{s\ge s_0,t}\mathbb{E}\left[f\left(w_{t+1}^s\right)-f\left(w_t^s\right)\mid {\mathcal{F}}_{t+1}^s\right]\right|\le \sum _{s\ge s_0,t}\left|\mathbb{E}\left[f\left(w_{t+1}^s\right)-f\left(w_t^s\right)\mid {\mathcal{F}}_{t+1}^s\right]\right|<\infty ,$$

$$p_{t+1}^s-p_t^s\le -2{\alpha }_t^s{⟨\mathrm{grad}f\left(w_t^s\right),\mathrm{Hess}f\left(w_t^s\right)\left[{\xi }_{t+1}^s\right]⟩}_{w_t^s}+{\left({\alpha }_t^s\right)}^2{\Vert {\xi }_{t+1}^s\Vert }_{w_t^s}^2{k}_2.$$

$$\mathbb{E}\left[p_{t+1}^s-p_t^s\mid {\mathcal{F}}_{t+1}^s\right]\le 2{\alpha }_t^s{\Vert \mathrm{grad}f\left(w_t^s\right)\Vert }_{w_t^s}^2{k}_3+{\left({\alpha }_t^s\right)}^2{C}^{\prime 2}{k}_2.$$

$$\Vert {\mathcal{T}}_{\eta }-{\mathcal{T}}_{R_{\eta }}\Vert \le c_0\Vert \eta \Vert ,\Vert {\mathcal{T}}_{\eta }^{-1}-{\mathcal{T}}_{R_{\eta }}^{-1}\Vert \le c_0\Vert \eta \Vert ,$$

$${\mathcal{T}}_{R_{\eta }}(\xi )=\mathrm{D}{R}_w(\eta )[\xi ],\eta ,\xi \in T_w\mathcal{M},w\in \mathcal{M}.$$

$${\Vert P_{\gamma }^{w\leftarrow z}(\mathrm{grad}f(z))-\mathrm{grad}f(w)\Vert }_w$$

$${\Vert {\mathcal{T}}_{\eta }(\xi )-P_{\eta }(\xi )\Vert }_z\le \theta \Vert \xi {\Vert }_w\Vert \eta {\Vert }_w,$$

Full text

Published in SIAM Journal on Optimization: https://epubs.siam.org/doi/abs/10.1137/17M1116787.

Keywords: Riemannian optimization, stochastic variance reduced gradient, retraction, vector transport, Riemannian centroid, principal component analysis, matrix completion

1 Introduction

A general loss minimization problem is defined as $\min_{w}f(w)$, where $f(w):=\frac{1}{N}\sum_{n=1}^{N}f_{n}(w)$, $w$ is the model variable, $N$ is the number of samples, and $f_{n}(w)$ is the loss incurred on the $n$-th sample. The full gradient descent (GD) algorithm requires the evaluation of $N$ derivatives, i.e., $\sum_{n=1}^{N}\nabla f_{n}(w)$, per iteration, which is computationally expensive when $N$ is extremely large. A well-known alternative uses only one derivative $\nabla f_{n}(w)$ per iteration for the $n$-th sample, and it forms the basis of the stochastic gradient descent (SGD) algorithm. When a relatively large step size is used in SGD, the training loss first decreases rapidly but results in large fluctuations around the solution. Conversely, when a small step size is used, a large number of iterations are required for SGD to converge. To circumvent this problem, SGD starts with a relatively large step size and gradually decreases it.

Recently, variance reduction techniques have been proposed to accelerate SGD convergence [7, 15, 21, 27, 29, 30, 32]. The stochastic variance reduced gradient (SVRG) algorithm is a popular technique with excellent convergence properties [15]. For smooth and strongly convex functions, SVRG has convergence rates similar to those of the stochastic dual coordinate ascent algorithm [30] and the stochastic average gradient (SAG) algorithm [27]. Garber and Hazan [9] analyzed the convergence rate of SVRG when $f$ is a convex function that is the sum of nonconvex (but smooth) terms, and they applied their result to the principal component analysis (PCA) problem. Shalev–Shwartz [28] also obtained similar results. Allen–Zhu and Yuan [3] further studied the same case with better convergence rates. Shamir [31] specifically studied the convergence properties of the variance reduction PCA algorithm. More recently, Allen-Zhu and Hazan [2] and Reddi et al. [26] independently proposed variance reduction methods for faster nonconvex optimization. However, it should be noted that all these cases assume a Euclidean search space.

In this paper, we handle problems in which the variables have a manifold structure:

[TABLE]

where $\mathcal{M}$ is a Riemannian manifold and $f_{n},n=1,2,\dots,N$ are real-valued functions on $\mathcal{M}$. These problems include, e.g., the low-rank matrix completion problem [23], the Riemannian centroid computation problem, and the PCA problem. In all these problems, optimization on Riemannian manifolds has shown state-of-the-art performance. The Riemannian framework exploits the geometry of the search space, which is characterized by the constraints of the optimization problem. Numerous efficient optimization algorithms have been developed [1]. Specifically, the problem $\min_{w\in\mathcal{M}}f(w)$, where $\mathcal{M}$ is a Riemannian manifold, is solved as an unconstrained optimization problem defined over the Riemannian manifold search space. Furthermore, although these algorithms mainly address batch-based approaches, Bonnabel [5] proposed a Riemannian stochastic gradient descent (R-SGD) algorithm that extends SGD from Euclidean space to Riemannian manifolds. Recently, more advanced stochastic optimization algorithms have also been generalized to Riemannian manifolds, including R-SQN-VR [18] and R-SRG [17].

Building upon the work of Bonnabel [5], we propose an extension of the stochastic variance reduced gradient algorithm to a Riemannian manifold search space (R-SVRG) and novel analyses. This extension is nontrivial and requires particular consideration for handling the averaging, addition, and subtraction of multiple gradients at different points on the manifold $\mathcal{M}$. Toward this end, this paper specifically leverages the notions of retraction and vector transport. The algorithm and convergence analysis presented in this paper are generalized, which in itself is a challenging problem, in the retraction and vector transport case, as well as in the exponential mapping and parallel translation case, allowing for extremely efficient implementation and making distinct contributions compared with an existing approach [34] that relies only on the exponential mapping and parallel translation case.

It should be mentioned that a recent study [34] by Zhang et al., which appeared simultaneously with our technical report [16], has also proposed R-SVRG on manifolds. The main difference between our work and [34] is that we provide convergence analyses for the algorithm with retraction and vector transport, whereas [34] deals with a special case in which exponential mapping and parallel translation are used as retraction and vector transport, respectively. There are additional differences. Our convergence analysis handles global and local convergence analyses separately, as in the typical analyses of batch algorithms on Riemannian manifolds [1]. Another difference is that our assumptions for the local convergence rate analysis are imposed only in a local neighborhood around a minimum, which is milder and more natural than the assumptions in [34], which assumes Lipschitz smoothness in the entire space. In other words, our global convergence analysis is not for a convergence rate or complexity, but it is an asymptotic convergence analysis. Here, according to classical usage in nonlinear programming, we use the term global convergence for convergence to a critical point from any initial point. On the other hand, our local convergence analysis is for a strongly convex function near the optimum, and the class of such functions includes many nonconvex functions that do not necessarily have global strong convexity. Consequently, our analysis should be applicable to different types of manifolds. For example, the parallel translation on the Stiefel manifold, which is an extremely important manifold in practice, is not available in a closed form. We can use a vector transport based on the orthogonal projection to the tangent space of the Stiefel manifold as an efficient implementation. Therefore, compared to the results of [34], our convergence results with retraction and vector transport enable us to deal with a wider variety of manifolds.

We emphasize that although we derive a local convergence rate with retraction and vector transport in this paper, it can be immediately used to derive a global iteration complexity if we assume, e.g., strong convexity of $f$ globally on the search space. Here, local convergence rate and global iteration complexity imply the rate at which the iterates approach a critical point in a sufficiently small neighborhood and the iteration complexity of obtaining a critical point from an arbitrary initial point, respectively. The result thus obtained can be regarded as a generalization of the global iteration complexity with strongly convex $f$ obtained in [34], where exponential mapping and parallel translation are used.

The remainder of this paper is organized as follows. Section 2 discusses Riemannian optimization theory, including the background on Riemannian geometry and some geometric tools used for optimization on Riemannian manifolds. Section 3 provides a detailed description of R-SVRG. Sections 4 and 5 present the global convergence analysis and local convergence rate analysis of the proposed R-SVRG, respectively. Section 6 highlights the superior performance of R-SVRG through numerical comparisons with R-SGD on three problems.

The proposed R-SVRG is implemented in the MATLAB toolbox Manopt [6]. The MATLAB code for all the proposed algorithms is available at https://github.com/hiroyuki-kasai/RSOpt.

2 Riemannian optimization

Optimization on Riemannian manifolds, or Riemannian optimization, seeks a critical point of a given real-valued function, called the objective or cost function, defined on a smooth Riemannian manifold $\mathcal{M}$. One of the advantages of using Riemannian geometry tools is that the intrinsic properties of the manifold enable us to handle constrained optimization problems as unconstrained optimization problems. This section introduces optimization on manifolds by summarizing [1]. Readers may refer to the various references therein as well as in [22, 24] for further details.

Let $f\colon\mathcal{M}\to\mathbb{R}$ be a smooth real-valued function on manifold $\mathcal{M}$. In optimization, we compute a minimum of $f$; typical methods for solving this minimization problem are iterative algorithms on manifold $\mathcal{M}$. In an iterative algorithm based on line search, with a given starting point $w_{0}\in\mathcal{M}$, we generate a sequence $\{w_{t}\}_{t\geq 0}$ on $\mathcal{M}$ that converges to $w^{*}$ whenever $w_{0}$ is in a neighborhood of $w^{*}$. In an iterative optimization algorithm, we compute a search direction and then move in the search direction. More specifically, an iteration on manifold $\mathcal{M}$ is performed by following geodesics emanating from $w_{t}$ and tangent to $\xi_{w_{t}}$ at $w_{t}$. The notion of geodesics on Riemannian manifolds is a generalized concept of straight lines in Euclidean space. For any tangent vector $\xi\in T_{w}\mathcal{M}$ at $w\in\mathcal{M}$, there exists an interval $I$ about [math] and a unique geodesic $\gamma_{e}(\cdot;w,\xi)\colon I\to\mathcal{M}$, such that $\gamma_{e}(0;w,\xi)=w$ and $\dot{\gamma_{e}}(0;w,\xi)=\xi$. The exponential mapping ${\rm Exp}_{w}\colon T_{w}\mathcal{M}\to\mathcal{M}$ at $w\in\mathcal{M}$ is defined by geodesics emanating from $w$ as ${\rm Exp}_{w}\xi=\gamma_{e}(1;w,\xi)$ for $\xi\in T_{w}\mathcal{M}$. If $\mathcal{M}$ is a complete manifold, the exponential mapping is defined for all vectors $\xi\in T_{w}\mathcal{M}$ [1, Section 5.4]. We can thus obtain an update formula using the exponential mapping

[TABLE]

where the search direction $\xi_{w_{t}}$ is in the tangent space $T_{w_{t}}\mathcal{M}$ of $\mathcal{M}$ at $w_{t}$, the scalar $s_{t}>0$ is the step size, and ${\rm Exp}_{w_{t}}(\cdot)$ is the exponential mapping, which induces a line search algorithm along the geodesics. In addition, given two points $w$ and $z$ on $\mathcal{M}$, the logarithm mapping, or simply log mapping, which is the inverse of the exponential mapping, maps $z$ to a vector $\xi\in T_{w}\mathcal{M}$ on the tangent space at $w$. The log mapping satisfies ${\rm dist}(w,z)=\|{\rm Log}_{w}(z)\|_{w}$, where ${\rm dist}\colon\mathcal{M}\times\mathcal{M}\to\mathbb{R}$ is the shortest distance between two points on $\mathcal{M}$ and $\|\cdot\|_{w}$ is the norm in $T_{w}\mathcal{M}$ defined through a Riemannian metric (see below).

The steepest descent or full gradient descent method for minimizing $f$ on $\mathcal{M}$ is an iterative algorithm obtained when $-{\rm grad}f(w_{t})$ is used as the search direction $\xi_{w_{t}}$. ${\rm grad}f(w_{t})$ is the Riemannian gradient of $f$ at $w_{t}$, which is computed according to the chosen metric $g$ at $w_{t}\in\mathcal{M}$. Collecting each metric $g_{w}\colon T_{w}\mathcal{M}\times T_{w}\mathcal{M}\to\mathbb{R}$ over $w\in\mathcal{M}$ gives a family called a Riemannian metric on $\mathcal{M}$. $g_{w}(\xi_{w},\zeta_{w})$ is an inner product of elements $\xi_{w}$ and $\zeta_{w}$ in the tangent space $T_{w}\mathcal{M}$ at $w$. Here, we use the notation $\langle\cdot,\cdot\rangle_{w}$ instead of $g(\cdot,\cdot)_{w}$ for simplicity. The gradient ${\rm grad}f(w)$ is defined as the unique element of $T_{w}\mathcal{M}$ that satisfies

[TABLE]

where ${\rm D}f(w)\colon T_{w}\mathcal{M}\to\mathbb{R}$ is the derivative of $f$ at $w$.

In searching for the next point along a geodesic, we need to compute tangent vectors obtained by the exponential mapping, which are expensive to compute in general. There are some Riemannian manifolds for which a closed form for the exponential mapping is not available. Alternatively, we can use curves other than geodesics as long as a starting point and its tangent vector at the initial time are the same as those of the geodesics. A more general update formula is then written as

[TABLE]

where $R_{w_{t}}$ is a retraction, which is any map $R_{w}\colon T_{w}\mathcal{M}\to\mathcal{M}$ that locally approximates the exponential mapping, up to the first order, on the manifold. The definition of retraction is as follows [1, Definition 4.1.1]:

Definition 2.1**.**

$R\colon T\mathcal{M}\to\mathcal{M}$* is called a retraction on $\mathcal{M}$ if the restriction $R_{w}\colon T_{w}\mathcal{M}\to\mathcal{M}$ to $T_{w}\mathcal{M}$ for any $w\in\mathcal{M}$ satisfies all of the following:*

$R_{w}(0_{w})=w$, where $0_{w}$ is the zero vector in $T_{w}\mathcal{M}$; 2. 2.

${\rm D}R_{w}(0_{w})[\xi]=\xi$* for any $\xi\in T_{w}\mathcal{M}$.*

Retractions include the exponential mapping as a special case. An advantage of using retractions is that the computational cost can be reduced compared to exponential mapping. It is worth noting that the convergence properties for the exponential mapping usually hold for retractions as well.

In the R-SVRG proposed in Section 3, we need to add tangent vectors that are in different tangent spaces, say $\tilde{w}$ and $w$ on $\mathcal{M}$. A mathematically natural way to do so is to use the parallel translation operator. Parallel translation $P_{\gamma}$ transports a vector field $\xi$ on the geodesic curve $\gamma$ that satisfies $P_{\gamma}^{a\leftarrow a}(\xi(a))=\xi(a)$ and $\frac{{\rm D}}{{\rm d}t}(P_{\gamma}^{t\leftarrow a}\xi(a))=0$ [1, Section 5.4], where $P_{\gamma}^{b\leftarrow a}$ is the parallel translation operator sending $\xi(a)$ to $\xi(b)$. However, parallel translation is sometimes computationally expensive, and no explicit formula is available for some manifolds, such as the Stiefel manifold. A vector transport on $\mathcal{M}$, which is a map $\mathcal{T}\colon T\mathcal{M}\oplus T\mathcal{M}\to T\mathcal{M}$, is used as an alternative. The definition of vector transport is as follows [1, Definition 8.1.1]:

Definition 2.2**.**

$\mathcal{T}\colon T\mathcal{M}\oplus T\mathcal{M}\to T\mathcal{M}$* is called a vector transport on $\mathcal{M}$ if it satisfies all of the following:*

$\mathcal{T}$* has an associated retraction $R$, i.e., for $w\in\mathcal{M}$ and $\xi,\eta\in T_{w}\mathcal{M}$, $\mathcal{T}_{\eta}(\xi)$ is a tangent vector at $R_{w}(\xi)$;* 2. 2.

$\mathcal{T}_{0_{w}}(\xi)=\xi$, where $\xi\in T_{w}\mathcal{M}$ and $w\in\mathcal{M}$; 3. 3.

$\mathcal{T}_{\eta}(a\xi+b\zeta)=a\mathcal{T}_{\eta}(\xi)+b\mathcal{T}_{\eta}(\zeta)$, where $a,b\in\mathbb{R}$, $\eta,\xi,\zeta\in T_{w}\mathcal{M}$, and $w\in\mathcal{M}$.

With a vector transport $\mathcal{T}$, $\mathcal{T}_{\eta}(\xi)$ can be regarded as a transported vector of $\xi$ along $\eta$. Parallel translation is an example of vector transport. In the following, we use the notations $P_{\eta}$ and $P_{\gamma}^{z\leftarrow w}$ interchangeably, where $\gamma$ is a curve connecting $w$ and $z$ on $\mathcal{M}$ defined by retraction $R$ as $\gamma(\tau):=R_{w}(\tau\eta)$ with $z=R_{w}(\eta)$. We also omit the subscript $\gamma$ when there is no confusion about the curve along which we transport a vector.

3 Riemannian stochastic variance reduced gradient algorithm

After a brief explanation of the variance reduced gradient variants in Euclidean space, we describe the proposed Riemannian stochastic variance reduced gradient algorithm on Riemannian manifolds.

3.1 Variance reduced gradient variants in Euclidean space

The SGD update in Euclidean space is $w_{t+1}=w_{t}-\alpha v_{t}$, where $v_{t}$ is a randomly selected vector called the stochastic gradient and $\alpha$ is the step size. SGD assumes an unbiased estimator of the full gradient as $\mathbb{E}_{n}[\nabla f_{n}(w_{t})]=\nabla f(w_{t})$. Many recent variants of the variance reduced gradient of SGD attempt to reduce its variance $\mathbb{E}[\|v_{t}-\nabla f(w_{t})\|^{2}]$ as $t$ increases to achieve better convergence [7, 15, 21, 27, 29, 30, 32]. SVRG, proposed in [15], introduces an explicit variance reduction strategy with double loops, where the $s$-th outer loop, called the $s$-th epoch, has $m_{s}$ inner iterations. SVRG first keeps $\tilde{w}^{s-1}=w_{m_{s-1}}^{s-1}$ or $\tilde{w}^{s-1}=w_{t}^{s-1}$ for randomly chosen $t\in\{1,2,\ldots,m_{s-1}\}$ at the end of the $(s\!\!-\!\!1)$-th epoch and also sets the initial value of the $s$-th epoch to $w_{0}^{s}=\tilde{w}^{s-1}$. It then computes a full gradient $\nabla f(\tilde{w}^{s-1})$. Subsequently, denoting the selected random index $i\in\{1,2,\ldots,N\}$ by $i_{t}^{s}$, SVRG randomly picks the $i_{t}^{s}$-th sample for each $t\geq 1$ at $s\geq 1$ and computes the modified stochastic gradient $v_{t}^{s}$ as

[TABLE]

It should be noted that SVRG can be regarded as a special case of S2GD (semi-stochastic gradient descent), which differs in terms of the number of inner loop iterations chosen [19].

3.2 Proposed Riemannian extension of SVRG (R-SVRG)

We propose a Riemannian extension of SVRG on a Riemannian manifold $\mathcal{M}$, called R-SVRG. Here, we denote the Riemannian stochastic gradient for the $i_{t}^{s}$-th sample as ${\rm grad}f_{i_{t}^{s}}$ and the modified Riemannian stochastic gradient as $\xi_{t}^{s}$ instead of $v_{t}^{s}$, in order to indicate the differences from the Euclidean case.

R-SVRG reduces the variance of the Riemannian stochastic gradient analogously to the SVRG algorithm in the Euclidean case. More specifically, R-SVRG keeps $\tilde{{w}}^{s-1}\in\mathcal{M}$ after $m_{s-1}$ stochastic update steps of the $(s\!\!-\!\!1)$-th epoch, and computes the full Riemannian gradient ${\rm grad}f(\tilde{{w}}^{s-1})=\frac{1}{N}\sum_{i=1}^{N}{\rm grad}f_{i}(\tilde{{w}}^{s-1})$ only for this stored $\tilde{{w}}^{s-1}$. The algorithm also computes the Riemannian stochastic gradient ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$ that corresponds to the $i_{t}^{s}$-th sample. Picking the $i_{t}^{s}$-th sample for each $t$-th inner iteration of the $s$-th epoch at ${w}_{t-1}^{s}$, we calculate $\xi_{t}^{s}$ in the same way as $v_{t}^{s}$ in (2), i.e., by modifying the stochastic gradient ${\rm grad}f_{i_{t}^{s}}({w}_{t-1}^{s})$ using both ${\rm grad}f(\tilde{{w}}^{s-1})$ and ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$. Translating the right-hand side of (2) to manifold $\mathcal{M}$ involves the sum of ${\rm grad}f_{i_{t}^{s}}({w}_{t-1}^{s})$, $-{\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$, and ${\rm grad}f(\tilde{{w}}^{s-1})$, which belong to two separate tangent spaces $T_{\scriptsize{w}_{t-1}^{s}}\mathcal{M}$ and $T_{\scriptsize\tilde{{w}}^{s-1}}\mathcal{M}$. This operation requires particular attention on a manifold, and a vector transport enables us to handle multiple elements on two separate tangent spaces flexibly. More specifically, ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$ and ${\rm grad}f(\tilde{{w}}^{s-1})$ are first transported to $T_{\scriptsize{w}_{t-1}^{s}}\mathcal{M}$ at the current point, ${w}_{t-1}^{s}$; then, they can be added to ${\rm grad}f_{i_{t}^{s}}({w}_{t-1}^{s})$ on $T_{\scriptsize{w}_{t-1}^{s}}\mathcal{M}$. Consequently, the modified Riemannian stochastic gradient $\xi_{t}^{s}$ at the $t$-th inner iteration of the $s$-th epoch is set to

[TABLE]

where $\mathcal{T}$ is a vector transport operator on $\mathcal{M}$ and $\tilde{\eta}_{t-1}^{s}\in T_{\tilde{w}^{s-1}}\mathcal{M}$ satisfies $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t-1}^{s})=w_{t-1}^{s}$. Specifically, we need to calculate the tangent vector from $\tilde{{w}}^{s-1}$ to ${w}_{t-1}^{s}$, which is given by the inverse of the retraction, i.e., $R^{-1}$, if available. Consequently, the final update rule of R-SVRG is defined as ${w}_{t}^{s}=R_{\scriptsize{w}_{t-1}^{s}}(-\alpha_{t-1}^{s}\xi_{t}^{s})$, where $\alpha_{t-1}^{s}>0$ is the step size at $w_{t-1}^{s}$. As shown in our local convergence analysis (Section 5.2), $\alpha_{t-1}^{s}$ can be fixed once the iterate becomes sufficiently close to a solution.

Let $\mathbb{E}_{i_{t}^{s}}[\cdot]$ be the expectation with respect to the random choice of $i_{t}^{s}$, conditioned on all the randomness introduced up to the $t$-th iteration of the inner loop of the $s$-th epoch. Conditioned on ${w}_{t-1}^{s}$, we take the expectation with respect to $i_{t}^{s}$ and obtain

[TABLE]

The theoretical convergence analysis of the Euclidean SVRG algorithm assumes that the initial vector ${w}_{0}^{s}$ of the $s$-th epoch is set to the average (or a random) vector of the $(s\!\!-\!\!1)$-th epoch [15, Figure 1]. However, the set of the last vectors in the $(s\!\!-\!\!1)$-th epoch, i.e., ${w}_{m_{s-1}}^{s-1}$, shows superior performance in the Euclidean SVRG algorithm. Therefore, for our local convergence rate analysis in Theorem 5.1, this study also uses, as option I, the mean value of $\tilde{{w}}^{s}=g_{m_{s}}({w}_{1}^{s},w_{2}^{s},\dots,{w}_{m_{s}}^{s})$ as $\tilde{{w}}^{s}$, where $g_{n}({w}_{1},w_{2},\ldots,{w}_{n})$ is the Riemannian centroid on the manifold. Alternatively, we can also simply choose $\tilde{{w}}^{s}={w}_{t}^{s}$ for $t\in\{1,2,\ldots,m_{s}\}$ at random. In addition, as option II, we can use the last iterate in the $s$-th epoch, i.e., $\tilde{{w}}^{s}={w}_{m_{s}}^{s}$. Note that option II is always of practical use because no additional computation is needed to obtain $\tilde{w}^{s}$. However, if we can compute the Riemannian centroid of $w_{1}^{s},w_{2}^{s},\dots,w_{m_{s}}^{s}$ relatively cheaply, option I may attain a better linear convergence rate. In the global convergence analysis in Section 4, we use option II, whereas both options are analyzed in the local convergence analysis in Section 5. The overall algorithm is summarized in Algorithm 1.

In addition, variants of the variance reduced SGD initially require a full gradient calculation at every epoch. This initially results in great overhead compared to the ordinary SGD algorithm and eventually induces a cold-start property in these variants. To avoid this, the use of the standard SGD update has been proposed only for the first epoch in Euclidean space [19]. We also adopt this simple modification of R-SVRG, denoted as R-SVRG+; we do not analyze this extension but leave it as an open problem.

As mentioned earlier, each iteration of R-SVRG has double loops to reduce the variance of the Riemannian stochastic gradient ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$. The $s$-th epoch, i.e., the outer loop, requires $N+2m_{s}$ gradient evaluations, where $N$ is for the full gradient ${\rm grad}f(\tilde{{w}}^{s-1})$ at the beginning of each $s$-th epoch and $2m_{s}$ is for the inner iterations, since each inner step needs two gradient evaluations, i.e., ${\rm grad}f_{i_{t}^{s}}({w}_{t-1}^{s})$ and ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$. However, if ${\rm grad}f_{i_{t}^{s}}(\tilde{{w}}^{s-1})$ for each sample is stored at the beginning of the $s$-th epoch, as in SAG, the evaluations over all the inner loops result in $m_{s}$. Finally, the $s$-th epoch requires $N+m_{s}$ evaluations. It is natural to choose an $m_{s}$ of the same order as $N$ but slightly larger (e.g., $m_{s}=5N$ for nonconvex problems has been suggested in [15]).

4 Global convergence analysis

In this section, we present a global convergence analysis of Algorithm 1 for Problem (1) after introducing some assumptions. Throughout this section, we let $R$ and $\mathcal{T}$ denote a retraction and vector transport used in Algorithm 1, respectively, and we make the following assumptions.

Assumption 4.1**.**

The retraction $R$ is such that $R_{w}\colon T_{w}\mathcal{M}\to\mathcal{M}$ for any $w\in\mathcal{M}$ is of class $C^{2}$. The vector transport $\mathcal{T}$ is continuous and isometric on $\mathcal{M}$, i.e., for any $w\in\mathcal{M}$ and $\xi,\eta\in T_{w}\mathcal{M}$, $\|\mathcal{T}_{\eta}(\xi)\|_{R_{w}(\eta)}=\|\xi\|_{w}$.

Note that we can construct an isometric vector transport as in [11, 12] such that Assumption 4.1 holds.

Assumption 4.2**.**

The objective function $f$ is thrice continuously differentiable and its components $f_{1},f_{2},\dots,f_{N}$ are twice continuously differentiable.

Assumption 4.3**.**

For a sequence $\{w_{t}^{s}\}$ generated by Algorithm 1, there exists a compact and connected set $K\subset\mathcal{M}$ such that $w_{t}^{s}\in K$ for all $s,t\geq 0$. In addition, for each $s\geq 1$ and $t\geq 0$, there exists $\tilde{\eta}_{t-1}^{s}\in T_{\tilde{w}^{s-1}}\mathcal{M}$ such that $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t-1}^{s})=w_{t-1}^{s}$. Furthermore, there exists $I>0$ such that for any $z\in K$, $R_{z}(\cdot)$ is defined in a ball $\mathbb{B}(0_{z},I)\subset T_{z}\mathcal{M}$, which is centered at the origin $0_{z}$ in $T_{z}\mathcal{M}$ with radius $I$.

Note that the existence of $\tilde{\eta}_{t-1}^{s}$ in Assumption 4.3 is guaranteed if $R_{z}(\cdot)$ for any $z\in K$ is a diffeomorphism in a ball $\mathbb{B}(0_{z},\rho)$ that satisfies $K\subset R_{z}(\mathbb{B}(0_{z},\rho))$. This assumption is a weakened version of the assumptions in some existing studies for special cases. For example, in [34], where the exponential mapping ${\rm Exp}$ is used as a retraction, the inverse of ${\rm Exp}$ is assumed to exist at all points $\{w_{t}^{s}\}$ in their global convergence analysis. Here, we note that if $I$ in Assumption 4.3 is too small, then $R_{\tilde{w}^{s-1}}^{-1}$ may not be defined in a neighborhood of $w_{t-1}^{s}$. We just suppose that there exists $\tilde{\eta}_{t-1}^{s}$ that is mapped to $w_{t-1}^{s}$ by $R_{\tilde{w}^{s-1}}$. Further analysis with respect to a specific manifold in question may be required to specify the value of $I$.

In the global convergence analysis, we also assume that the sequence of step sizes $\{\alpha_{t}^{s}\}_{t\geq 0,s\geq 1}$ satisfies the usual condition in stochastic approximation as follows:

Assumption 4.4**.**

The sequence $\{\alpha_{t}^{s}\}$ of step sizes satisfies

[TABLE]

where $\sum$ denotes $\sum_{s=1}^{\infty}\sum_{t=0}^{m_{s}}$.

This condition is satisfied, e.g., if $\{m_{s}\}$ is upper bounded and $\alpha_{t}^{s}=\alpha_{0}(1+\alpha_{0}\lambda\lfloor k/m_{s}\rfloor)^{-1}$ with positive constants $\alpha_{0}$ and $\lambda$, where $k$ is the total iteration number depending on $s$ and $t$ and $\lfloor\cdot\rfloor$ denotes the floor function. We also note the following proposition introduced in [8]:

Proposition 4.1** ([8]).**

Let $(X_{n})_{n\in\mathbb{N}}$ be a nonnegative stochastic process with bounded positive variations, i.e., such that $\sum_{n=0}^{\infty}\mathbb{E}[\mathbb{E}[X_{n+1}-X_{n}\mid\mathcal{F}_{n}]^{+}]<\infty$, where $X^{+}$ denotes the quantity $\max\{X,0\}$ for a random variable $X$ and $\mathcal{F}_{n}$ is the increasing sequence of $\sigma$-algebras generated by the variables just before time $n$. Then, the process is a quasi-martingale, i.e.,

[TABLE]

Now, we give the a.s. convergence of the proposed algorithm under the assumption that the generated sequence is in a compact set.

Theorem 4.1**.**

Suppose Assumptions 4.1–4.3 and consider Algorithm 1 with option II and step sizes $\{\alpha_{t}^{s}\}$ satisfying Assumption 4.4 on a Riemannian manifold $\mathcal{M}$. If $f\geq 0$, then $\{f(w_{t}^{s})\}$ converges a.s. and ${\rm grad}f(w_{t}^{s})\to 0$ a.s.

Proof.

The claim is proved similarly to the proof of the standard Riemannian SGD (see [5]). Since $K$ is compact, all continuous functions on $K$ can be bounded. Therefore, there exists $C>0$ such that for all $w\in K$ and $n\in\{1,2,\ldots,N\}$, we have $\|{\rm grad}f(w)\|_{w}\leq C$ and $\|{\rm grad}f_{n}(w)\|_{w}\leq C$. We use $C^{\prime}:=3C$ in the following. Moreover, as $\alpha_{t}^{s}\to 0$, there exists $s_{0}$ such that for $s\geq s_{0}$, we have $\alpha_{t}^{s}<1$ and $\alpha_{t}^{s}C^{\prime}<I$. Now, suppose that $s\geq s_{0}$. Let $\tilde{\eta}_{t}^{s}$ satisfy $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t}^{s})=w_{t}^{s}$. The existence of such $\tilde{\eta}_{t}^{s}$ is guaranteed from Assumption 4.3. It follows from the triangle inequality that

[TABLE]

where we have used the assumption that $\mathcal{T}$ is an isometry. Therefore, we have

[TABLE]

Hence, it follows from Assumption 4.3 that there exists a curve $R_{w_{t}^{s}}(-\tau\alpha_{t}^{s}\xi_{t+1}^{s})_{0\leq\tau\leq 1}$ linking $w_{t}^{s}$ and $w_{t+1}^{s}$. Defining $g(\tau;w,\xi):=(f\circ R_{w})(-\tau\xi)$, there exists a constant $k_{1}$ such that

[TABLE]

since $[0,1]\times\{(w,\xi)\mid\xi\in\mathbb{B}(0_{w},I),w\in K\}$ is compact. A trivial relation

[TABLE]

together with $\alpha_{t}^{s}<1$ then implies that

[TABLE]

Let $\mathcal{F}_{t}^{s}$ be the increasing sequence of $\sigma$-algebras defined by

[TABLE]

Since $w_{t}^{s}$ is computed from $i_{1}^{1},i_{2}^{1},\ldots,i_{t}^{s}$, it is measurable in $\mathcal{F}_{t+1}^{s}$. As $i_{t+1}^{s}$ is independent of $\mathcal{F}_{t+1}^{s}$, we have

[TABLE]

from the linearity of $\mathcal{T}_{\tilde{\eta}_{t}^{s}}$. Therefore, it holds that

[TABLE]

which yields

[TABLE]

We reindex the sequence $\{w_{0}^{1},w_{1}^{1},\dots,w_{m_{1}}^{1}(=w_{0}^{2}),w_{1}^{2},w_{2}^{2},\dots,w_{m_{2}}^{2},\dots,w_{t}^{s},\dots\}$

as $\{w_{1},w_{2},\dots,w_{k},\dots\}$. We also similarly reindex $\{\alpha_{t}^{s}\}$. As $f(w_{k})\geq 0$, (5) proves that $\{f(w_{k})+\sum_{l=k}^{\infty}(\alpha_{l})^{2}k_{1}\}$ is a nonnegative supermartingale. Hence, it converges a.s. Returning to the original indexing, this implies that $\{f(w_{t}^{s})\}$ converges a.s. Moreover, summing the inequalities, we have

[TABLE]

Here, we prove that the right-hand side is bounded. By doing so, we can conclude the convergence of the left term.

Summing (5) over $s$ and $t$, we have

[TABLE]

by Assumption (4). It then follows from Proposition 4.1 that

[TABLE]

which together with (6) implies that $\sum_{s\geq s_{0},t}\alpha_{t}^{s}\|{\rm grad}f(w_{t}^{s})\|_{w_{t}^{s}}^{2}$ converges a.s. If we further prove that $\{\|{\rm grad}f(w_{t}^{s})\|_{w_{t}^{s}}\}$ converges a.s., it can only converge to 0 a.s. because of Eq. (4).

As in the discussion in [5], we consider the nonnegative process $p_{t}^{s}=\|{\rm grad}f(w_{t}^{s})\|_{w_{t}^{s}}^{2}$. Bounding the largest eigenvalue of the Hessian of $\|{\rm grad}f(w)\|_{w}^{2}$ from above by $k_{2}$ on the compact set $K$ along the curve defined by the retraction $R$ linking $w_{t}^{s}$ and $w_{t+1}^{s}$, a Taylor expansion yields

[TABLE]

Let $k_{3}$ be a constant such that $-k_{3}$ is a lower bound of the minimum eigenvalue of the Hessian of $f$ on $K$. Then, we have

[TABLE]

We have proved the fact that the infinite series of the terms of the right-hand side is finite, implying that $\{p_{t}^{s}\}$ is a quasi-martingale, which further implies the a.s. convergence of $\{p_{t}^{s}\}$ to [math], as claimed. ∎

Theorem 4.1 includes a global convergence of the algorithm with exponential mapping and parallel translation as a special case. In this case, a sufficient condition for the assumptions can be easily written by the notions of injectivity radius as the following corollary.

Corollary 4.1**.**

Suppose Assumption 4.1 and consider Algorithm 1 with option II and step sizes $\{\alpha_{t}^{s}\}$ satisfying Assumption 4.4 on a Riemannian manifold $\mathcal{M}$, where exponential mapping and parallel translation are used as retraction and vector transport, respectively. Assume that $\mathcal{M}$ is connected and has an injectivity radius uniformly bounded from below by a positive constant. Assume also that there exists a compact and connected set $K\subset\mathcal{M}$ such that $w_{t}^{s}\in K$ for all $s,t\geq 0$. If $f\geq 0$, then $\{f(w_{t}^{s})\}$ converges a.s. and ${\rm grad}f(w_{t}^{s})\to 0$ a.s.

5 Local convergence rate analysis

In this section, we show the local convergence rate analysis of the R-SVRG algorithm; we analyze the convergence of any sequences generated by Algorithm 1 that are contained in a sufficiently small neighborhood of a local minimum point of the objective function. Hence, we can assume that the objective function is strongly convex in such a neighborhood. We first give formal expressions of these assumptions and then analyze the local convergence rate of the algorithm.

5.1 Assumptions and existing lemmas

We again suppose Assumption 4.1 and make the following assumption, which is a weaker version of Assumption 4.2.

Assumption 5.1**.**

The objective function $f$ and its components $f_{1},f_{2},\dots,f_{N}$ are twice continuously differentiable.

Let $w^{*}$ be a critical point of $f$. As discussed in [12], for a positive constant $\rho$, a $\rho$-totally retractive neighborhood $\Omega$ of $w\in\mathcal{M}$ is a neighborhood such that for all $z\in\Omega$, $\Omega\subset R_{z}(\mathbb{B}(0_{z},\rho))$, and $R_{z}(\cdot)$ is a diffeomorphism on $\mathbb{B}(0_{z},\rho)$, which is the ball in $T_{z}\mathcal{M}$ with center $0_{z}$ and radius $\rho$, where $0_{z}$ is the zero vector in $T_{z}\mathcal{M}$. The concept of a totally retractive neighborhood is analogous to that of a totally normal neighborhood for exponential mapping. Note that for each point, a $\rho$-totally retractive neighborhood for sufficiently small $\rho>0$ is proved to exist [12]. On the other hand, it should be noted that a concrete bound of the radii of locally retractive neighborhoods is difficult to specify. In the local convergence analysis, we assume the following.

Assumption 5.2**.**

The sequence $\{w_{t}^{s}\}$ generated by Algorithm 1 continuously remains in totally retractive neighborhood $\Omega$ of critical point $w^{*}$ of $f$, i.e., $R_{w_{t}^{s}}(-\alpha\xi_{t+1}^{s})\in\Omega$ for all $s,t\geq 0$ and for all $\alpha\in[0,\alpha_{t}^{s}]$. In addition, the radius of $\Omega$ is sufficiently small.

We note that the assumption on the radius of $\Omega$ in Assumption 5.2 is for Assumption 5.4 and for Lemma 5.7.

Assumption 5.3**.**

There exists a constant $c_{0}$ such that vector transport $\mathcal{T}$ satisfies the following conditions:

[TABLE]

where $\|\cdot\|$ denotes the induced (operator) norm of the Riemannian metric and $\mathcal{T}_{R}$ denotes the differentiated retraction, i.e.,

[TABLE]

Assumption 5.3 states that the difference between $\mathcal{T}$ and $\mathcal{T}_{R}$ is small if the tangent vector along which a vector is transported is close to [math]. See [12] for further details.

Furthermore, we assume the following.

Assumption 5.4**.**

$f$* is strongly retraction-convex with respect to $R$ in $\Omega$; i.e., there exist two constants $0<a_{0}<a_{1}$ such that $a_{0}\leq\frac{d^{2}}{d\alpha^{2}}f(R_{w}(\alpha\eta))\leq a_{1}$ for all $w\in\Omega$, $\eta\in T_{w}\mathcal{M}$ with $\|\eta\|_{w}=1$, and for all $\alpha$ satisfying $R_{w}(\tau\eta)\in\Omega$ for all $\tau\in[0,\alpha]$. Furthermore, $f$ is strongly geodesically convex in $\Omega$, i.e., $f$ is strongly retraction-convex with respect to ${\rm Exp}$.*

Letting $\Omega$ be smaller if necessary, we can guarantee the last assumption from the other assumptions using Lemma 3.1 in [12].

The next assumption states how close the modified stochastic gradient by retraction $R$ is to that obtained by the exponential mapping.

Assumption 5.5**.**

Let $\{\xi_{t}^{s}\}$ be generated by Algorithm 1 with a fixed step size of $\alpha_{t}^{s}:=\alpha$. ${\rm Exp}_{w}^{-1}$ for any $w\in\Omega$ is defined in $\Omega$ and there exists constant $c_{R}>0$ such that $\|{\rm Exp}_{w_{t-1}^{s}}^{-1}(w_{t}^{s})-(-\alpha\xi_{t}^{s})\|_{w_{t-1}^{s}}\leq c_{R}\|\alpha\xi_{t}^{s}\|_{w_{t-1}^{s}}^{2}$.

If $R_{w_{t-1}^{s}}^{-1}$ is defined in a ball whose image by $R_{w_{t-1}^{s}}$ contains $\Omega$, then the above assumption translates into $\|{\rm Exp}_{w_{t-1}^{s}}^{-1}(w_{t}^{s})-R_{w_{t-1}^{s}}^{-1}(w_{t}^{s})\|_{w_{t-1}^{s}}\leq c_{R}\|R_{w_{t-1}^{s}}^{-1}(w_{t}^{s})\|_{w_{t-1}^{s}}^{2}$, which is natural when we assume sufficient smoothness of ${\rm Exp}$ and $R$ because $R$ coincides with ${\rm Exp}$ up to the first order.

In the remainder of this section, we introduce some existing lemmas to evaluate the differences between using retraction and vector transport and using exponential mapping and parallel translation, and the effects of the curvature of the manifold in question.

Lemma 5.1** (In the proof of Lemma 3.9 in [12]).**

Under Assumptions 5.1 and 5.2, there exists a constant $\beta>0$ such that

[TABLE]

where $w$ and $z$ are in $\Omega$ and $\gamma$ is a curve $\gamma(\tau):=R_{z}(\tau\eta)$ for an arbitrary $\eta\in T_{z}\mathcal{M}$ defined by retraction $R$ on $\mathcal{M}$. $P_{\gamma}^{w\leftarrow z}(\cdot)$ is a parallel translation operator along curve $\gamma$ from $z$ to $w$.

Note that curve $\gamma$ in this lemma is not necessarily the geodesic on $\mathcal{M}$. Relation (7) is a generalization of the Lipschitz continuity condition.

Lemma 5.2** (Lemma 3.5 in [12]).**

Let $\mathcal{T}\in C^{0}$ be a vector transport associated with the same retraction $R$ as that of the parallel transport $P\in C^{\infty}$. Under Assumption 5.3, for any $\bar{w}\in\mathcal{M}$, there exists a constant $\theta>0$ and neighborhood $\mathcal{U}$ of $\bar{w}$ such that for all $w,z\in\mathcal{U}$,

[TABLE]

where $\xi,\eta\in T_{w}\mathcal{M}$ and $R_{w}(\eta)=z$.

We can derive the following lemma from the Taylor expansion as in the proof of Lemma 3.2 in [12].

Lemma 5.3** (In the proof of Lemma 3.2 in [12]).**

Under Assumptions 5.2–5.4, there exists a positive real number $\sigma$ such that

[TABLE]

Proof.

From the assumptions, there exists $\sigma>0$ such that $\frac{d^{2}}{d\alpha^{2}}f({\rm Exp}_{w}(\alpha\eta))\geq\sigma$ for all $w\in\Omega$, $\eta\in T_{w}\mathcal{M}$ with $\|\eta\|_{w}=1$, and for all $\alpha$ such that ${\rm Exp}_{w}(\tau\eta)\in\Omega$ for all $\tau\in[0,\alpha]$. From Taylor’s theorem, we can conclude that this $\sigma$ satisfies the claim. ∎

If we replace ${\rm Exp}$ in Lemma 5.3 with retraction $R$ on $\mathcal{M}$, we can obtain a similar result for $R$, which is shown in Lemma 3.2 in [12] where the constant $a_{0}$ corresponds to $\sigma$ in our Lemma 5.3. However, Lemma 5.3 for ${\rm Exp}$ is sufficient in the following discussion.

From Lemma 3 in [11], we have the following:

Lemma 5.4** (Lemma 3 in [11]).**

Let $\mathcal{M}$ be a Riemannian manifold endowed with retraction $R$ and let $\bar{w}\in\mathcal{M}$. Then, there exist $\mu>0$, $\nu>0$, and $\delta_{\mu,\nu}>0$ such that for all $w$ in a sufficiently small neighborhood of $\bar{w}$ and all $\xi\in T_{w}\mathcal{M}$ with $\|\xi\|_{w}\leq\delta_{\mu,\nu}$, the inequalities

[TABLE]

hold.

Since we have ${\rm dist}(w,{\rm Exp}_{w}(\xi))=\|\xi\|_{w}$, (9) is equivalent to ${\rm dist}(w,{\rm Exp}_{w}(\xi))\leq\mu{\rm dist}(w,R_{w}(\xi))$ and ${\rm dist}(w,R_{w}(\xi))\leq\nu{\rm dist}(w,{\rm Exp}_{w}(\xi))$, which give relations between the exponential mapping and a general retraction.

Now, we introduce Lemma 6 in [35] to evaluate the distance between $w_{t}^{s}$ and $w^{*}$ using the smoothness of our objective function. In the following, an Alexandrov space is defined as a length space whose curvature is bounded.

Lemma 5.5** (Lemma 6 in [35]).**

If $a$, $b$, and $c$ are the side lengths of a geodesic triangle in an Alexandrov space with curvature lower-bounded by $\kappa$, and $A$ is the angle between sides $b$ and $c$, then

[TABLE]

5.2 Local convergence rate analysis with retraction and vector transport

We now demonstrate the local convergence properties of the R-SVRG algorithm (i.e., local convergence to local minimizers) and its convergence rate. The main theorem (Theorem 5.1) follows three lemmas, the proofs of which are provided in Appendix A.

The following lemma shows a property of the Riemannian centroid on a general Riemannian manifold.

Lemma 5.6**.**

Let $w_{1},w_{2},\dots,w_{m}$ be points on a Riemannian manifold $\mathcal{M}$ and let $w$ be the Riemannian centroid of the $m$ points. For an arbitrary point $p$ on $\mathcal{M}$, we have

[TABLE]

Then, we state that the norms of modified stochastic gradients in R-SVRG are sufficiently small under Assumptions 4.1 and 5.1–5.3.

Lemma 5.7**.**

Under Assumptions 4.1 and 5.1–5.3, the norm of $\xi_{t}^{s}$ computed by (3) is sufficiently small; i.e., for any $\varepsilon>0$, there exists $r_{0}>0$ such that $\|\xi_{t}^{s}\|_{w_{t-1}^{s}}\leq\varepsilon$ for $r$ satisfying $r<r_{0}$, where $r$ is the radius of $\Omega$.

We now provide the upper bound of the variance of $\xi_{t}^{s}$ as follows.

Lemma 5.8**.**

Suppose Assumptions 4.1 and 5.1–5.3, which guarantee Lemmas 5.1, 5.2, and 5.4 for $\bar{w}=w^{*}$. Let $\beta>0$ be a constant such that

[TABLE]

The existence of such a $\beta$ is guaranteed by Lemma 5.1. The upper bound of the variance of $\xi_{t}^{s}$ is given by

[TABLE]

where the constant $\theta$ corresponds to that in Lemma 5.2, $C$ is the upper bound of $\|{\rm grad}f_{n}(w)\|,\ n=1,2,\dots,N$ for $w\in\Omega$, and $\mu>0$ appears in (9).

We proceed to the main theorem and prove it using Lemmas 5.6–5.8.

Theorem 5.1**.**

Let $\mathcal{M}$ be a Riemannian manifold whose curvature is lower-bounded by $\kappa$, $w^{*}\in\mathcal{M}$ be a nondegenerate local minimizer of $f$ (i.e., ${\rm grad}f(w^{*})=0$ and the Hessian ${\rm Hess}f(w^{*})$ of $f$ at $w^{*}$ is positive definite), $D$ be the diameter of the compact set $\Omega$, and $\zeta:=\sqrt{|\kappa|}D/\tanh(\sqrt{|\kappa|}D)$ if $\kappa<0$ and $\zeta:=1$ if $\kappa\geq 0$. Suppose Assumptions 4.1 and 5.1–5.5. Let $c_{R}$ be a constant in Assumption 5.5, let $\beta$, $\mu$, $\theta$, and $C$ be the same as in Lemma 5.8, and let $\nu$ be a constant in Lemma 5.4. Suppose that $\sigma>0$ is a constant in Lemma 5.3 satisfying (8). Let $\alpha$ be a positive number satisfying $0<\alpha(\sigma-28(\zeta\nu+2c_{R}D)(\beta^{2}+\mu^{2}C^{2}\theta^{2})\alpha)<1$. It then follows that for any sequence $\{\tilde{w}^{s}\}$ generated by Algorithm 1 with a fixed step size $\alpha_{t}^{s}:=\alpha$ and $m_{s}:=m$ converging to $w^{*}$, there exists $K>0$ such that for all $s>K$,

[TABLE]

where

[TABLE]

Before proving the theorem, we summarize the above theorem as the following corollary:

Corollary 5.1**.**

Let $\mathcal{M}$ be a Riemannian manifold whose curvature is lower-bounded and $w^{*}\in\mathcal{M}$ be a nondegenerate local minimizer of $f$. Suppose Assumptions 4.1 and 5.1–5.4, and let $\alpha$ be a positive number. If $\alpha$ is sufficiently small, for any sequence $\{\tilde{w}^{s}\}$ generated by Algorithm 1 with a fixed step size $\alpha_{t}^{s}:=\alpha$ and $m_{s}:=m$ converging to $w^{*}$, there exists $K>0$ such that for all $s>K$,

[TABLE]

where $\delta$ is a constant in $(0,1)$.

Proof of Theorem 5.14.

Since the function $x/\tanh(x)$ on $x$ monotonically increases in $[0,\infty)$, we have, from Lemma 5.5,

[TABLE]

for any geodesic triangle in $\Omega$ with side lengths $a$, $b$, and $c$, since $\sqrt{|\kappa|}c\leq\sqrt{|\kappa|}D$. Then, conditioned on $w_{t-1}^{s}$, the expectation of the distance between $w_{t}^{s}$ and $w^{*}$ with respect to the random choice of $i_{t}^{s}$ is evaluated by considering the geodesic triangle with $w_{t-1}^{s}$, $w^{*}$, and $w_{t}^{s}$ in $\Omega$ as

[TABLE]

Here, we have the following from Assumption 5.5:

[TABLE]

where the relation $\|{\rm Exp}_{w_{t-1}^{s}}^{-1}(w^{*})\|_{w_{t-1}^{s}}={\rm dist}(w_{t-1}^{s},w_{*})\leq D$ is incorporated. It follows that

[TABLE]

where the last inequality follows from $\mathbb{E}_{i_{t}^{s}}[\xi_{t}^{s}]={\rm grad}f(w_{t-1}^{s})$. Lemma 5.3, together with the relation $f(w^{*})\leq f(w_{t-1}^{s})$, yields

[TABLE]

We thus obtain, by Lemma 5.8,

[TABLE]

It follows for the unconditional expectation operator $\mathbb{E}$ that

[TABLE]

Summing (13) over $t=1,2,\ldots,m$ of the inner loop on the $s$-th epoch, we have

[TABLE]

Hence, with option II, where we compute $\tilde{w}^{s}$ as $\tilde{w}^{s}=w_{m}^{s}$, the facts $w_{0}^{s}=\tilde{w}^{s-1}$ and $\alpha(28(\zeta\nu+2c_{R}D)\alpha(\beta^{2}+\mu^{2}C^{2}\theta^{2})-\sigma)<0$ imply that

[TABLE]

On the other hand, if we use option I, where $\tilde{w}^{s}=g_{m}(w_{1}^{s},w_{2}^{s},\dots,w_{m}^{s})$, rearranging (14) yields

[TABLE]

Using $\tilde{w}^{s}=g_{m}(w_{1}^{s},w_{2}^{s},\ldots,w_{m}^{s})$ and Lemma 5.6, we obtain

[TABLE]

This completes the proof. ∎

In the above theorem, we note that, from the definitions of $\beta$ and $\sigma$, $\beta$ can be chosen to be arbitrarily large and $\sigma$ can be chosen to be arbitrarily small. Therefore, $\alpha=\sigma/56(\zeta\nu+2c_{R}D)(\beta^{2}+\mu^{2}C^{2}\theta^{2})$, e.g., satisfies $0<\alpha(\sigma-28(\zeta\nu+2c_{R}D)(\beta^{2}+\mu^{2}C^{2}\theta^{2})\alpha)<1$ for sufficiently large $\beta$ and small $\sigma$.

In fact, $\alpha$ satisfying the condition always exists for any values of $\beta$ and $\sigma$. Let $\beta^{\prime}:=28(\zeta\nu+2c_{R}D)(\beta^{2}+\mu^{2}C^{2}\theta^{2})$. We can analyze the inequality $0<\alpha(\sigma-28(\zeta\nu+2c_{R}D)(\beta^{2}+\mu^{2}C^{2}\theta^{2})\alpha)<1$, which is expressed as $0<\alpha(\sigma-\beta^{\prime}\alpha)<1$, with the condition $\alpha>0$ more specifically as

[TABLE]

Furthermore, we can show that the coefficient $\delta$ on the right-hand side of (11), which can be written as $4(7+4m\beta^{\prime}\alpha^{2})/7m\alpha(\sigma-\beta^{\prime}\alpha)=:r_{1}(\alpha)$ for Option I and $1-\sigma\alpha+(1+4m/7)\beta^{\prime}\alpha^{2}=:r_{2}(\alpha)$ for Option II, is less than $1$ when $m$ is sufficiently large and $\alpha$ is appropriately chosen. If $\alpha$ is fixed, $r_{1}(\alpha)\to 16\beta^{\prime}\alpha/7(\sigma-\beta^{\prime}\alpha)$ as $m\to\infty$, which is not necessarily less than $1$, and $r_{2}(\alpha)\to\infty$ as $m\to\infty$. Thus, we again need to specifically analyze an appropriate value of $\alpha$, which should depend on $m$.

For Option I, by calculating the derivative $r_{1}^{\prime}(\alpha)$ of $r_{1}(\alpha)$ on $\alpha$, we can show that $r_{1}(\alpha)$ takes the minimum value at

[TABLE]

which satisfies (15) when $m$ is sufficiently large, since $\lim_{m\to\infty}\alpha_{*}=0$. Note that we have $r_{1}^{\prime}(\alpha_{*})=0$, which yields $4m\beta^{\prime}\sigma\alpha_{*}^{2}+14\beta^{\prime}\alpha_{*}-7\sigma=0$. This relation gives the minimum value $r_{1}(\alpha_{*})$ as

[TABLE]

where we have used the facts that $\lim_{m\to\infty}\alpha_{*}=0$ and $\lim_{m\to\infty}m\alpha_{*}=\infty$. A simpler choice of $\alpha=1/\sqrt{m}$ also makes $r_{1}(\alpha)$ less than $1$ if $m$ is sufficiently large since

[TABLE]

Although $\alpha=1/\sqrt{m}$ does not achieve the best rate attained by $\alpha=\alpha_{*}$, this choice is practical because we do not know the exact values of $\sigma$, $\beta$, or $\alpha_{*}$ in general.

A similar discussion can be applied to the case of Option II. In this case, it is clear that a sufficiently small $\alpha$ satisfies (15) and $r_{2}(\alpha)<1$. Furthermore, if $\sigma^{2}-4\beta^{\prime}\leq 0$, $\alpha_{*}=\sigma/2(1+4m/7)\beta^{\prime}$ satisfies (15) and attains the minimum value of $r_{2}$ as

[TABLE]

Furthermore, this $\alpha_{*}$ satisfies (15) if $\sigma^{2}-4\beta^{\prime}>0$ and $m$ is sufficiently large.

We have thus shown that a local linear convergence rate is achieved under an appropriate fixed step size if $m$ is sufficiently large, which is the same as standard SVRG in Euclidean space (for nonconvex problems). We can also analyze the rate with decaying step sizes $\alpha_{0}^{s}>\alpha_{1}^{s}>\dots>\alpha_{m}^{s}$ (at the $s$-th epoch) as

[TABLE]

for Option I. This is larger than

[TABLE]

which is the coefficient in (11) with the fixed step size $\alpha=\alpha_{0}^{s}$. A similar discussion can be applied to the case of Option II. Consequently, using decaying step sizes also yields a local convergence, but at a worse rate than with a fixed step size. Both the above guarantees are quite similar to those available for batch gradient algorithms on manifolds. This setup, i.e., hybrid step sizes, follows our two convergence analyses, which first require decaying step sizes to approach a neighborhood of a local minimum and use a fixed step size to achieve faster linear convergence near the solution. As mentioned earlier, we guarantee global convergence and local linear convergence even if we use decaying step sizes from beginning to end. Therefore, this method is an improved version of decaying step size. Theoretical analysis of the switching between decaying and fixed step sizes is left for future work.

We make one more remark. Although this subsection describes the local convergence analysis with an objective function that is strictly retraction-convex in a sufficiently small neighborhood of a local minimizer, it is also applicable for a global convergence analysis if we assume that the function is globally strictly convex, as in other studies. As we have already discussed, the rate in Theorem 5.1 can lead to discussions on global iteration complexity. The key aspect of our local convergence analysis is that we do not have to assume global convexity to attain the local linear convergence rate.

5.3 Local convergence rate analysis with exponential mapping and parallel translation

In this subsection, we present a local convergence rate analysis of the R-SVRG algorithm with exponential mapping and parallel translation along the geodesics. This is a special case of the previous subsection where exponential mapping and parallel translation are chosen as retraction and vector transport, respectively. However, in this particular case, we can obtain a stricter rate than in a general case. Since the results are obtained by a similar discussion, we give a sketch of the proofs.

We obtain the following result as a corollary of the proof of Lemma 5.8, with $R={\rm Exp}$ and $\mathcal{T}=P$.

Corollary 5.2**.**

Suppose the conditions in Lemma 5.8 and consider Algorithm 1 with $R={\rm Exp}$ and $\mathcal{T}=P$, i.e., the exponential mapping and parallel translation case. Let $\beta$ be a constant in Lemma 5.8. Then, the upper bound of $\mathbb{E}_{i_{t}^{s}}[\|\xi_{t}^{s}\|_{w_{t-1}^{s}}^{2}]$ is given by

[TABLE]

Proof.

Putting $R={\rm Exp}$ and $\mathcal{T}=P$ in the middle of the proof of Lemma 5.8, which is in Appendix A, we obtain

[TABLE]

In a similar manner to Lemma 5.8, we have

[TABLE]

This completes the proof. ∎

Corollary 5.3**.**

Suppose the conditions in Theorem 5.1, except that a positive number $\alpha$ satisfies $0<\alpha(\sigma-14\zeta\beta^{2}\alpha)<1$, and consider Algorithm 1 with a fixed step size $\alpha_{t}^{s}:=\alpha$ for the exponential mapping and parallel translation case. For any sequence $\{\tilde{{w}}^{s}\}$ generated by the algorithm, there exists $K>0$ such that for all $s>K$,

[TABLE]

where

[TABLE]

Proof.

Note that constants in Theorem 5.1 are $c_{R}=\theta=0$ and $\mu=\nu=1$ in this case. By using (16) instead of (10), we obtain

[TABLE]

instead of (12). Summing over $t=1,2,\ldots,m$ of the inner loop on the $s$-th epoch, we have

[TABLE]

A similar discussion as in the proof of Theorem 5.1 yields the claimed convergence rates. ∎

6 Numerical comparisons

This section compares the performance of R-SVRG(+) (with option II) with that of the Riemannian extension of SGD, i.e., R-SGD, where the Riemannian stochastic gradient algorithm uses ${\rm grad}f_{i_{t}^{s}}(w_{t-1}^{s})$ instead of $\xi_{t}^{s}$ in (3). We also make a comparison with R-SD, which is the Riemannian steepest descent algorithm with backtracking line search [1, Chapters 4]. We consider both fixed step size and decaying step size sequences. The decaying step size sequence uses the decay $\alpha_{k}=\alpha_{0}(1+\alpha_{0}\lambda\lfloor k/m_{s}\rfloor)^{-1}$, where $k$ is the number of iterations. We select some values of $\alpha_{0}$ and consider three values of $\lambda$ ($10^{-1},10^{-2},\text{ and }10^{-3}$). In addition, since the global convergence analysis needs a decaying step size condition and the local convergence rate analysis holds for a fixed step size (Sections 4 and 5), we consider a hybrid step size sequence that follows the decaying step size before the $s_{\mathrm{TH}}$ epoch and subsequently switches to a fixed step size. All the experiments use $m_{s}=5N$ by following [15]. In all the figures, the $x$-axis represents the computational cost measured by the number of gradient computations divided by $N$. The algorithms are initialized randomly and are terminated when either the stochastic gradient norm is below $10^{-8}$ or the number of iterations exceeds a predefined threshold. It should be noted that all the results except those of R-SD are the best-tuned results. All the simulations were performed in MATLAB on a 2.6 GHz Intel Core i7 machine with 16 GB RAM. Hereinafter, this paper addresses three problems on the SPD manifold and the Grassmann manifold. In all the problems, full gradient methods, e.g., the steepest descent algorithm, become prohibitively computationally expensive when $N$ is extremely large. The stochastic gradient approach is a promising way to achieve scalability.

6.1 Problem on the SPD manifold and simulation results

We first consider the Riemannian centroid problem on the SPD manifold.

SPD manifold $\mathcal{S}_{++}^{d}$ and optimization tools. We designate the space of $d\times d$ SPD matrices as the SPD manifold, $\mathcal{S}_{++}^{d}$. If we endow $\mathcal{S}_{++}^{d}$ with the affine-invariant Riemannian metric (AIRM) [25] defined by $\langle\xi_{\scriptsize\mbox{\bf X}},\eta_{\scriptsize\mbox{\bf X}}\rangle_{\scriptsize\mbox{\bf X}}=\rm{trace}(\xi_{\scriptsize\mbox{\bf X}}\mbox{\bf X}^{-1}\eta_{\scriptsize\mbox{\bf X}}\mbox{\bf X}^{-1})$ for $\xi_{\scriptsize\mbox{\bf X}},\eta_{\scriptsize\mbox{\bf X}}\in T_{\scriptsize\mbox{\bf X}}\mathcal{S}_{++}^{d}$ at $\mbox{\bf X}\in\mathcal{S}_{++}^{d}$, the SPD manifold $\mathcal{S}_{++}^{d}$ forms a Riemannian manifold. The exponential mapping is written as

[TABLE]

for any $\xi_{\scriptsize\mbox{\bf X}}$ and X. The parallel translation of $\xi_{\scriptsize\mbox{\bf X}}$ along $\eta_{\scriptsize\mbox{\bf X}}$ on $\mathcal{S}_{++}^{d}$ is given by $P_{\eta_{\tiny\mbox{\bf X}}}(\xi_{\scriptsize\mbox{\bf X}})=\mbox{\bf X}^{1/2}\mbox{\bf Y}\mbox{\bf X}^{-1/2}\xi_{\scriptsize\mbox{\bf X}}\mbox{\bf X}^{-1/2}\mbox{\bf Y}\mbox{\bf X}^{1/2}$, where $\mbox{\bf Y}=\exp(\mbox{\bf X}^{-1/2}\eta_{\scriptsize\mbox{\bf X}}\mbox{\bf X}^{-1/2}/2)$. The logarithm map of Y at X is described as

[TABLE]

The exponential mapping and the parallel translation above are computationally expensive. Therefore, an efficient retraction is proposed as [14];

[TABLE]

This maps $\xi_{\scriptsize\mbox{\bf X}}$ onto $\mathcal{S}_{++}^{d}$ for all $\xi_{\scriptsize\mbox{\bf X}}\in T_{\scriptsize\mbox{\bf X}}\mathcal{S}_{++}^{d}$. Huang et al. proposed an efficient isometric vector transport [11, 33] defined as

[TABLE]

where $\mbox{\bf Y}=R_{\scriptsize\mbox{\bf X}}(\xi_{\scriptsize\mbox{\bf X}})$ and $a^{\flat}$ denotes the flat of $a\in T_{w}\mathcal{M}$, i.e., $a^{\flat}\colon T_{w}\mathcal{M}\to\mathbb{R}\colon v\mapsto\langle a,v\rangle_{w}$. $B_{\scriptsize\mbox{\bf X}}$ and $B_{\scriptsize\mbox{\bf Y}}$ are the orthonormal bases of $T_{\scriptsize\mbox{\bf X}}\mathcal{S}_{++}^{d}$ and $T_{\scriptsize\mbox{\bf Y}}\mathcal{S}_{++}^{d}$, respectively, where the basis is calculated based on the Cholesky decomposition. Consequently, the implementation of our algorithm for this particular problem uses the retraction (17) and the vector transport (18), which satisfy the requirements in the convergence analyses in Sections 4 and 5.

Riemannian centroid problem. We first evaluate the proposed algorithm in the Riemannian centroid problem on $\mathcal{S}_{++}^{d}$, which is frequently used for computer vision problems, such as visual object categorization and pose categorization [13]. Given $N$ points on $\mathcal{S}_{++}^{d}$ with matrix representations $\mbox{\bf X}_{1},\mbox{\bf X}_{2},\dots,\mbox{\bf X}_{N}$, the Riemannian centroid is derived from the solution to the problem

[TABLE]

where ${\rm dist}(\mbox{\bf A},\mbox{\bf B})=\|\log(\mbox{\bf A}^{-1/2}\mbox{\bf B}\mbox{\bf A}^{-1/2})\|_{F}$ represents the distance along the corresponding geodesic between the two points $\mbox{\bf A},\mbox{\bf B}\in\mathcal{S}_{++}^{d}$ with respect to the AIRM. The gradient of the loss function is computed as $\frac{1}{N}\sum_{n=1}^{N}-{\rm log}(\mbox{\bf X}_{n}\mbox{\bf C}^{-1})\mbox{\bf C}$.

Figures 1(a) and (b) show the results of the optimality gap and the norm of the gradient, respectively, where $N=1000$ and $d=3$. The choices of $\alpha_{0}$ are $\{10^{-3},2\times 10^{-3},4\times 10^{-3},6\times 10^{-3},8\times 10^{-3},10^{-2},2\times 10^{-2},\ldots,10^{-1}\}$. $s_{\mathrm{TH}}$ and the batch size are fixed to $3$ and $1$, respectively. The maximum number of iterations is $10$ for R-SVRG(+) and $60$ for the others. The optimality gap evaluates the performance against the minimum loss, which is obtained by R-SD with high precision in advance. From the figures, R-SVRG(+) outperforms R-SGD in terms of the gradient counts and exhibits much faster convergence than R-SD as expected.

6.2 Problems on the Grassmann manifold and simulation results

We focus on two popular problems on the Grassmann manifold: PCA and low-rank matrix completion problems.

Grassmann manifold and optimization tools. An element on the Grassmann manifold is represented by a $d\times r$ orthogonal matrix U with orthonormal columns, i.e., $\mbox{\bf U}^{T}\mbox{\bf U}=\mbox{\bf I}$. Two orthogonal matrices represent the same element on the Grassmann manifold if they are related by right multiplication of an $r\times r$ orthogonal matrix $\mbox{\bf O}\in\mathcal{O}(r)$, where $\mathcal{O}(r)$ is the orthogonal group of order $r$. Equivalently, an element of the Grassmann manifold is identified with a set of $d\times r$ orthogonal matrices $[\mbox{\bf U}]:=\{\mbox{\bf U}\mbox{\bf O}\mid\mbox{\bf O}\in\mathcal{O}(r)\}$. Thus, ${\rm Gr}(r,d):={\rm St}(r,d)/\mathcal{O}(r)$, where ${\rm St}(r,d)$ is the Stiefel manifold, which is the set of matrices of size $d\times r$ with orthonormal columns. The Grassmann manifold has the structure of a Riemannian quotient manifold [1, Section 3.4]. The exponential mapping for the Grassmann manifold from $\mbox{\bf U}(0):=\mbox{\bf U}\in{\rm Gr}(r,d)$ in the direction of $\xi\in T_{\scriptsize\mbox{\bf U}(0)}{\rm Gr}(r,d)$ is given in a closed form as [1, Section 5.4]

[TABLE]

where $\xi=\mbox{\bf W}\Sigma\mbox{\bf V}^{T}$ is the rank-$r$ singular value decomposition of $\xi$. The $\cos(\cdot)$ and $\sin(\cdot)$ operations are only on the diagonal entries. The parallel translation of $\zeta\in T_{\scriptsize\mbox{\bf U}(0)}{\rm Gr}(r,d)$ on the Grassmann manifold along $\gamma(t)$ with $\dot{\gamma}(0)=\mbox{\bf W}\Sigma\mbox{\bf V}^{T}$ is given in a closed-form by

[TABLE]

The logarithm map of $\mbox{\bf U}(t)$ at $\mbox{\bf U}(0)$ on the Grassmann manifold is given by

[TABLE]

where the rank-$r$ singular value decomposition of $(\mbox{\bf U}(t)-\mbox{\bf U}(0)\mbox{\bf U}(0)^{T}\mbox{\bf U}(t))\allowbreak(\mbox{\bf U}(0)^{T}\mbox{\bf U}(t))^{-1}$ is $\mbox{\bf W}\Sigma\mbox{\bf V}^{T}$. It should be noted that this experiment evaluates the projection-based vector transport and the QR-decomposition-based retraction, which do not satisfy the conditions in Sections 4 and 5, but are computationally efficient. The intention here is to show that our algorithm performs well empirically without using the specific vector transport.

The PCA problem. Given an orthonormal matrix projector $\mbox{\bf U}\in{\rm St}(r,d)$, the PCA problem is to minimize the sum of the squared residual errors between the projected data points and the original data, as

[TABLE]

where \mbox{\boldmathx}_{n} is a data vector of size $d\times 1$. Problem (21) is equivalent to maximizing the function \frac{1}{N}\sum_{n=1}^{N}\mbox{\boldmathx}_{n}^{T}\mbox{\bf U}\mbox{\bf U}^{T}\mbox{\boldmathx}_{n}. Here, the critical points in the space ${\rm St}(r,d)$ are not isolated because the cost function remains unchanged under the group action $\mbox{\bf U}\mapsto\mbox{\bf UO}$ for any orthogonal matrices O of size $r\times r$. Consequently, Problem (21) is reformulated as an optimization problem on the Grassmann manifold ${\rm Gr}(r,d)$.

Figures 2(a) and (b) show the optimality gap and gradient norm, respectively, where $N=10000$, $d=20$, and $r=5$. The choices of $\alpha_{0}$ are $\{10^{-3},2\times 10^{-3},\ldots,10^{-2}\}$. The minimum loss for the optimality gap evaluation is obtained by the MATLAB function pca. $s_{\mathrm{TH}}$ and the batch size are fixed to 5 and 10, respectively. The maximum number of iterations is $16$ for R-SVRG(+) and $100$ for the others. From Figure 2(a), we can observe that between R-SVRG and R-SVRG+, the latter shows superior performance for all the step size sequences. In Figure 2(b), while the gradient norm of SGD remains at higher values, those of R-SVRG and R-SVRG+ converge to lower values in all the cases.

Low-rank matrix completion. The matrix completion problem amounts to completing an incomplete matrix X, say of size $d\times N$, from a small number of entries by assuming a low-rank model for the matrix. If $\Omega$ is the set of indices for which we know the entries in X, the rank-$r$ matrix completion problem amounts to solving the problem

[TABLE]

where the operator $\mathcal{P}_{\Omega}$ acts as $\mathcal{P}_{\Omega}(\mbox{\bf X}_{ij})=\mbox{\bf X}_{ij}$ if $(i,j)\in\Omega$ and $\mathcal{P}_{\Omega}(\mbox{\bf X}_{ij})=0$ otherwise. This is called the orthogonal sampling operator and is a mathematically convenient way to represent the subset of known entries. Partitioning \mbox{\bf X}=[\mbox{\boldmathx}_{1},\mbox{\boldmathx}_{2},\ldots,\mbox{\boldmathx}_{n}], Problem (22) is equivalent to the problem

[TABLE]

where \mbox{\boldmathx}_{n}\in\mathbb{R}^{d} and the operator $\mathcal{P}_{\Omega_{n}}$ is the sampling operator for the $n$-th column. Given U, \mbox{\boldmatha}_{n} in (23) admits a closed form solution. Consequently, Problem (23) depends only on the column space of U and is on the Grassmann manifold [4].

The proposed algorithms are also compared with Grouse [4], a state-of-the-art stochastic descent algorithm on the Grassmann manifold. We first consider a synthetic dataset with $N=5000$ and $d=500$ with rank $r=5$. The algorithms are initialized randomly as suggested in [20]. The ten choices of $\alpha_{0}$ are $\{10^{-3},2\times 10^{-3},\ldots,10^{-2}\}$ for R-SGD and R-SVRG(+), and $\{0.1,0.2,\ldots,1.0\}$ for Grouse. This instance considers the loss on a test set $\Gamma$, which differs from training set $\Omega$. We also impose an exponential decay of the singular values. The condition number (CN) of a matrix is the ratio of the largest to the smallest singular values of the matrix. The over-sampling ratio (OS) determines the number of entries that are known. This instance uses CN=$5$ and OS=$5$. An OS of $5$ implies that $5(N+d-r)r$ randomly and uniformly selected entries are known a priori among the total $Nd$ entries. Figure 3(a) shows the results of loss on the test set $\Gamma$. These results show the superior performance of our proposed algorithms.

Next, we consider Jester dataset 1 [10], consisting of ratings of $100$ jokes by $24983$ users. Each rating is a real number from $-10$ to $10$. We randomly extract two ratings per user as the training set $\Omega$ and test set $\Gamma$. The algorithms are run by fixing the rank to $r=5$ with random initialization. $\alpha_{0}$ is chosen from $\{10^{-6},2\times 10^{-6},\ldots,10^{-5}\}$ for SGD and SVRG(+) and $\{10^{-3},2\times 10^{-3},\ldots,10^{-2}\}$ for Grouse. Figure 3(b) shows the superior performance of R-SVRG(+) on the test set of the Jester dataset.

As a final test, we compare the algorithms on the MovieLens-1M dataset downloaded from http://grouplens.org/datasets/movielens/. The dataset has a million ratings corresponding to $6040$ users and $3952$ movies. $\alpha_{0}$ is chosen from $\{10^{-5},2\times 10^{-5},\ldots,10^{-4}\}$. Figure 3(c) shows the results on the test set for all the algorithms except Grouse, which faces issues with convergence on this data set. R-SVRG(+) shows much faster convergence than the others, and R-SVRG is better than R-SVRG+ in terms of the final test loss for all step size algorithms.

7 Conclusion

We proposed a Riemannian stochastic variance reduced gradient (R-SVRG) algorithm with retraction and vector transport, which includes the algorithm with exponential mapping and parallel translation as a special case. The proposed algorithm stems from the variance reduced gradient algorithm in Euclidean space, but it has been extended herein to Riemannian manifolds. The main challenges of averaging, adding, and subtracting multiple gradients on a Riemannian manifold were handled with a vector transport. We proved that R-SVRG generates globally convergent sequences with a decaying step size and is locally linearly convergent with a fixed step size under some natural assumptions. Numerical comparisons of problems on the SPD manifold and the Grassmann manifold indicated the superior performance of R-SVRG on various benchmarks.

Appendix A Proofs of lemmas

In this section, we present complete proofs of Lemmas 5.6–5.8.

Proof of Lemma 5.6.

From the triangle inequality and $(a+b)^{2}\leq 2a^{2}+2b^{2}$ for real numbers $a$ and $b$, we have for $i=1,2,\ldots,m$,

[TABLE]

Since $w$ is the Riemannian centroid of $w_{1},w_{2},\ldots,w_{m}$, it holds that

[TABLE]

It then follows that

[TABLE]

This completes the proof. ∎

Proof of Lemma 5.7.

We first show that $\|{\rm grad}f_{n}\|$ for any $n\in\{1,2,\dots,N\}$ is upper-bounded in $\Omega$ when $\Omega$ is sufficiently small. Since $\Omega$ is sufficiently small, it is contained in a set $U\subset\mathcal{M}$ diffeomorphic to an open set $U^{\prime}\subset\mathbb{R}^{\dim\mathcal{M}}$ by a chart $\phi\colon U\to U^{\prime}$. Consider a sufficiently small closed ball $B$ in $U^{\prime}$ centered at $\phi(w^{*})$ such that $\phi^{-1}(B)\subset U$. Then, $w^{*}$ is in $\phi^{-1}(B)$. Note that $B$ is compact and $\phi$ is a diffeomorphism. Hence, $\phi^{-1}(B)$ is also compact. Replacing $\Omega$ with this $\phi^{-1}(B)$ if necessary, we can assume that $\Omega$ is compact. Therefore, $\|{\rm grad}f_{n}\|$ is upper bounded in $\Omega$. In this proof, let $C$ denote a constant such that $\|{\rm grad}f_{n}(z)\|_{z}\leq C$ for all $n\in\{1,2,\dots,N\}$ and $z\in\Omega$.

Let $\tilde{\eta}_{t-1}^{s}$ satisfy $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t-1}^{s})=w_{t-1}^{s}$. The definition of $\xi_{t}^{s}$, Eq. (3), and the triangle inequality yield

[TABLE]

where $\gamma_{t-1}^{s}(\tau)=R_{\tilde{w}^{s-1}}(\tau\tilde{\eta}_{t-1}^{s})$ and $P$ is the parallel translation operator along curve $\gamma_{t-1}^{s}$. The three terms on the right-hand side above are bounded as follows. For the first term, it follows from Lemma 5.1 that

[TABLE]

where $\beta$ is in the lemma. The second term can be evaluated from Lemmas 5.2 and 5.4 as

[TABLE]

where $\theta$ and $\mu$ are in the lemmas. From Lemma 5.1, we have

[TABLE]

for the third term, where $P_{*}$ is the parallel translation along curve $\tilde{\gamma}_{*}^{s-1}$ defined by $\tilde{\gamma}_{*}^{s-1}(\tau)=R_{w^{*}}(\tau\tilde{\eta}_{*}^{s-1})$ with $\tilde{\eta}_{*}^{s-1}$ satisfying $R_{w^{*}}(\tilde{\eta}_{*}^{s-1})=\tilde{w}^{s-1}$. Therefore, we have

[TABLE]

if we choose a sufficiently small $r$ such that $r<\varepsilon/(3\beta+2\theta C\mu)$. ∎

Proof of Lemma 5.8.

Let $\eta_{t-1}^{*s}\in T_{w^{*}}\mathcal{M}$ and $\tilde{\eta}_{t-1}^{s}\in T_{\tilde{w}^{s-1}}\mathcal{M}$ satisfy $R_{w^{*}}(\eta_{t-1}^{*s})=w_{t-1}^{s}$ and $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t-1}^{s})=w_{t-1}^{s}$, respectively. Let $P^{w\leftarrow z}$ be a parallel translation along the curve $R_{z}(\tau\eta)$, where $R_{z}(\eta)=w$. By Lemmas 5.1 and 5.2, the upper bound of $\mathbb{E}_{i_{t}^{s}}[\|\xi_{t}^{s}\|_{w_{t-1}^{s}}^{2}]$ in terms of the distance of $w_{t-1}^{s}$ and $\tilde{w}^{s-1}$ from $w^{*}$ is computed as

[TABLE]

where the relation $\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$ for vectors $a$ and $b$ in a norm space and the triangle inequality are used repeatedly. Note also that $\mathbb{E}_{i_{t}^{s}}[{\rm grad}f_{i_{t}^{s}}(\tilde{w}^{s-1})]={\rm grad}f(\tilde{w}^{s-1})$ and ${\rm grad}f(w^{*})=0$ and that $\mathbb{E}_{i_{t}^{s}}$ is a linear operator. Furthermore, we have evaluated the value $\mathbb{E}_{i_{t}^{s}}[\|{\rm grad}f_{i_{t}^{s}}(w_{t-1}^{s})-\mathcal{T}_{\eta_{t-1}^{*s}}({\rm grad}f_{i_{t}^{s}}(w^{*}))\|_{w_{t-1}^{s}}^{2}$ and again used the obtained relation in the third inequality. ∎