Riemannian stochastic variance reduced gradient
Hiroyuki Sato\thanks{Department of Information and Computer Technology
Tokyo University of Science
Tokyo
Japan ({\tt hsato@rs.tus.ac.jp
Stochastic variance reduction algorithms have recently become popular for minimizing the average of a large but finite number of loss functions. In this paper, we propose a novel Riemannian extension of the Euclidean stochastic variance reduced gradient algorithm (R-SVRG) to a manifold search space. The key challenges of averaging, adding, and subtracting multiple gradients are addressed with retraction and vector transport. We present a global convergence analysis of the proposed algorithm with a decay step size and a local convergence rate analysis under a fixed step size under some natural assumptions. The proposed algorithm is applied to problems on the Grassmann manifold, such as principal component analysis, low-rank matrix completion, and computation of the Karcher mean of subspaces, and outperforms the standard Riemannian stochastic gradient descent algorithm in each case\footnote{This paper extends the earlier work \cite{Kasai_arXiv_2016 to include more general results..
R-SVRG_SIOPT_arXiv.tex
\documentclass[11pt,a4paper]{article}
\bibliographystyle{siam}
\usepackage{a4wide}
%\usepackage[truedimen,margin=25truemm]{geometry}

%\renewcommand{\theequation}{\thesection.\arabic{equation}}
%\renewcommand{\thefigure}{\thesection.\arabic{figure}}
%\renewcommand{\thetable}{\thesection.\arabic{table}}
%\renewcommand{\thealgorithm}{\thesection.\arabic{algorithm}}

% For figures
\usepackage{graphicx} % more modern
%\usepackage{epstopdf}

% For citations
%\usepackage{natbib}

% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}


% As of 2011, we use the hyperref package to produce hyperlinks in the
% resulting PDF.  If this breaks your system, please commend out the
% following usepackage line and replace \usepackage{icml2016} with
% \usepackage[nohyperref]{icml2016} above.
%\usepackage{hyperref}
\usepackage{url}
%\usepackage{breakurl}
%\PassOptionsToPackage{hyphens}{url}\usepackage{hyperref}


%%%%%%% HK %%%%%%%%%%%%%%%%%%%%%%%%

%% For math
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
% \usepackage{mathtools}
% \usepackage[raggedright]{titlesec} % avoid hyphen in sections
% \newtheorem{prop}{Proposition}\newtheorem{lemma}{Lemma}
\newtheorem{assumption}{Assumption}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{lemma}{Lemma}[section]

\usepackage{cite}

%% Vector, Matrix
\def\vec#1{\mbox{\boldmath $#1$}}
\def\mat#1{\mbox{\bf #1}}
\def\L#1{\mbox{$\ell_{#1}$}}		% for l1,l2-norm


%%
\newcommand{\gradf}{{\rm grad} f}
\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}


% Personal macros, pick colors according to your choice
%\newcommand{\changeBM}[1]{\textcolor{magenta}{#1}}
%\newcommand{\changeHK}[1]{\textcolor{red}{#1}}
%\newcommand{\changeHS}[1]{\textcolor{blue}{#1}}
\newcommand{\changeHKK}[1]{#1}
\newcommand{\changeHSS}[1]{{#1}}
\newcommand{\changeBMM}[1]{{#1}}

% Uncomment these when for the final submission 
\newcommand{\changeBM}[1]{#1}
\newcommand{\changeHK}[1]{#1}
%\newcommand{\changeHKK}[1]{#1}
\newcommand{\changeHS}[1]{#1}
%\newcommand{\changeHSS}[1]{#1}
%\newcommand{\changeBMM}[1]{#1}
\newcommand{\changeHSSS}[1]{#1}

\newcommand{\changehs}[1]{#1}
\newcommand{\changehss}[1]{#1}


\title{Riemannian stochastic variance reduced gradient}
\date{\today}


%\date{Originally Compiled on May 24th, 2016\\
%Compiled on \today}

\author{Hiroyuki Sato\thanks{Department of Information and Computer Technology, Tokyo University of Science, Tokyo, Japan ({\tt hsato@rs.tus.ac.jp}).} \and Hiroyuki Kasai\thanks{Graduate School of Informatics and Engineering, The University of Electro-Communications, Tokyo, Japan ({\tt kasai@is.uec.ac.jp}).} \and Bamdev Mishra\thanks{\changehss{Core Machine Learning Team, Amazon.com, Bangalore, India} ({\tt bamdevm@amazon.com}.)}}

\begin{document}

\maketitle



\begin{abstract}
Stochastic variance reduction algorithms have recently become popular for minimizing the average of a large but finite number of loss functions. In this paper, we propose a novel Riemannian extension of the Euclidean stochastic variance reduced gradient algorithm (R-SVRG) to a manifold search space. The key challenges of averaging, adding, and subtracting multiple gradients are addressed with retraction and vector transport. We present a global convergence analysis of the proposed algorithm with a decay step size and a local convergence rate analysis under a fixed step size under some natural assumptions. The proposed algorithm is applied to problems on the Grassmann manifold, such as principal component analysis, low-rank matrix completion, and computation of the Karcher mean of subspaces, and outperforms the standard Riemannian stochastic gradient descent algorithm in each case\footnote{This paper extends the earlier work \cite{Kasai_arXiv_2016} to include more general results.}.
\end{abstract}

\noindent {\bf Keywords:} 
Riemannian optimization, stochastic variance reduced gradient, retraction, vector transport, matrix completion

\section{Introduction}
\label{Sec:intro}

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\changeBM{,} and $f_n(w)$ is the loss incurred on \changehs{the} $n$-th sample. The {\it full gradient descent} (GD) algorithm requires evaluating $N$ derivatives, i.e., $\sum_{n=1}^N \nabla f_n(w)$, per iteration, which is computationally expensive when $N$ is very large. A popular alternative uses only one derivative $\nabla f_n(w)$ per iteration for \changehs{the} $n$-th sample, and forms the basis of the {\it stochastic gradient descent} (SGD) algorithm. When a relatively large step size is used in SGD, train loss first decreases rapidly, \changeBM{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. \changeBM{To circumvent this issue, SGD starts with a relatively large step size and gradually decreases it.}


\changeBM{Recently, {\it variance reduction} techniques have been proposed to} accelerate SGD convergence \cite{Johnson_NIPS_2013_s, Mairal_SIAMJOPT_2015,Roux_NIPS_2012_s,Shalev_arXiv_2012_s,Shalev_JMLR_2013_s,Defazio_NIPS_2014_s,Zhang_SIAMJO_2014_s}.
The stochastic \changeHK{v}ariance \changeHK{r}educed \changeHK{g}radient \changeHK{(SVRG)} is a popular technique with superior convergence properties \cite{Johnson_NIPS_2013_s}. \changeBM{For smooth and strongly convex functions,} SVRG \changeBM{has convergence rates similar} to those of stochastic dual coordinate ascent  \cite{Shalev_JMLR_2013_s} and stochastic average gradient (SAG) algorithms \cite{Roux_NIPS_2012_s}. Garber and Hazan \cite{Garber_arXiv_2015_s} analyzed the convergence rate for SVRG when $f$ is a convex function that is the sum of nonconvex (but smooth) terms and applied this result \changeBM{to} principal component analysis (PCA) problem. Shalev-Shwartz  \cite{Shalev_arXiv_2015_s} also found similar results. Allen-Zhu and Yuan \cite{Allen-Zhu-arXiv_2015_s} further studied the same case with better convergence rates. Shamir \cite{Shamir_arXiv_2015_s} specifically studied the convergence properties of the variance reduction PCA algorithm. Very recently, Allen-Zhu and Hazan \cite{Allen-Zhu-arXiv_2016_s} proposed a variance reduction method for faster nonconvex optimization. \changeBM{However, it should be noted that all these cases assume a Euclidean search space.}


In this paper, we handle problems in which \changeBM{the} variables have a manifold structure:
\begin{equation}
\label{Prob:general}
\min_{w \in \mathcal{M}} f(w):=\frac{1}{N}\sum_{n=1}^N f_n(w),
\end{equation}
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, for example, low-rank matrix completion problem \cite{Mishra_ICDC_2014_s}, Karcher mean computation problem, and PCA problem. In all these problems, optimization on {\it Riemannian manifolds} has shown state-of-the-art performance. The Riemannian framework exploits the geometry of the \changehss{search space which are characterized by constraints of the optimization problem. A lot of efficient optimization algorithms have been developed} \cite{Absil_OptAlgMatManifold_2008}. Specifically, the problem $\min_{w \in \mathcal{M}} f(w)$, where $\mathcal{M}$ is a Riemannian manifold, is solved as an \emph{unconstrained optimization problem} defined over the Riemannian manifold search space. Bonnabel \cite{Bonnabel_IEEETAC_2013_s} proposed a {\it Riemannian stochastic gradient descent} algorithm (R-SGD) that extends SGD from Euclidean space to Riemannian manifolds.

It should be mentioned that recent work \cite{Zhang_NIPS_2016}, which appeared simultaneously with our technical report \cite{Kasai_arXiv_2016}, also proposes R-SVRG on manifolds. The main difference \changehs{between} our work \changehs{and} \cite{Zhang_NIPS_2016} is that we provide convergence analyses for the algorithm with retraction and vector transport, \changehs{whereas} \cite{Zhang_NIPS_2016} deals with a special case \changehs{in which} exponential mapping and parallel translation are used as retraction and vector transport, respectively.
There are further differences; our convergence analysis handles global and local convergence analyses separately, as in typical analyses \changehs{of} batch algorithms on Riemannian manifolds \cite{Absil_OptAlgMatManifold_2008}. Another difference is that our assumptions for the local rate of convergence analysis are imposed only in a local neighborhood around a minimum, which is milder and more natural than the assumptions in \cite{Zhang_NIPS_2016}, which assumes Lipschitz smoothness in the entire space. Consequently, our analysis should be applicable to \changehs{a} wider \changehs{variety} of manifolds than \changehs{that of} \cite{Zhang_NIPS_2016}.

\changeBM{Building upon the work of Bonnabel \cite{Bonnabel_IEEETAC_2013_s} and \cite{Zhang_NIPS_2016}, we propose \changehs{an} extension of the stochastic variance reduction gradient algorithm to Riemannian manifold search space (R-SVRG) and novel analyses.}
This extension is nontrivial and requires particular consideration in handling the averaging, addition, and subtraction of multiple gradients at different points on the manifold $\mathcal{M}$. To this end, this paper specifically \changeHKK{leverages the notions of \changeHSSS{retraction and vector transport}.
The algorithm and convergence analysis in this paper are challengingly generalized in the retraction and vector transport case, as well as in the exponential mapping and parallel translation case,
allowing extremely efficient implementation and making distinct contributions compared with existing works \cite{Zhang_NIPS_2016} that rely only on the exponential mapping and parallel translation case.}

The paper is organized as follows. Section \ref{Sec:GrassmannAndProblems} discusses Riemannian optimization \changehs{theory, including background on Riemannian geometry and some geometric tools used in optimization on Riemannian manifolds}.
The detailed description of R-SVRG \changehs{is} given in Section \ref{Sec:R-SVRG}.
\changeHSSS{Sections \ref{Sec:GlobalAnalysis} and \ref{Sec:LocalAnalysis} present the global convergence analysis and the local convergence rate analysis of the proposed R-SVRG, respectively.} In Section \ref{Sec:NumericalComparison}, numerical comparisons with R-SGD on three problems suggest \changehs{the} superior performance of R-SVRG.

Our proposed R-SVRG is implemented in the Matlab toolbox Manopt \cite{Boumal_Manopt_2014_s}. The Matlab codes for the proposed algorithms are available at \url{https://bamdevmishra.com/codes/RSVRG/}.

%Additional numerical experiments are shown in Section \ref{AppenSec:AdditionalNumericalComparison}.
%As part of supplementary material, we also provide an implementation of R-SVRG compatible to Manopt (Matlab file). 


%
%
%
\section{\changeHKK{Riemannian optimization}}
\label{Sec:GrassmannAndProblems}
Optimization on manifolds, or {\it Riemannian \changeHS{optimization}}, seeks a global or local optimum of a \changehss{given} real-valued function\changehss{, called objective or cost function,} defined \changehss{on} a smooth \changehss{Riemannian} manifold $\mathcal{M}$. One of the advantages of using Riemannian geometry tools is that the intrinsic properties of the manifold allow constrained optimization problems to be handled as unconstrained optimization problems.
This section introduces optimization on manifolds by summarizing \cite{Absil_OptAlgMatManifold_2008}.
We refer to the many references therein, and in \cite{Mishra_SIAMJOPT_2016,Meyer_ICML_2011}, for further detail.

Let $f: \mathcal{M} \rightarrow \mathbb{R}$ be a smooth real-valued function on manifold $\mathcal{M}$.
\changehss{In optimization, we compute a minimum of $f$; usual methods for solving this minimization problem are {\it iterative algorithms} on manifold $\mathcal{M}$. In an iterative algorithm, with a given starting point $\changeHS{w_0} \in \mathcal{M}$, we generate a sequence $\{w_t\}_{t\geq0}$ on $\mathcal{M}$ that converges to $w^{*}$ whenever $\changeHS{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 concretely, an iteration on manifold $\mathcal{M}$ is performed by following geodesics \changehss{emanating} from $w_t$ and tangent to $\xi_{w_t}$ \changehss{at $w_t$}.
\changehss{Notion of} {\it Geodesics} on \changehss{Riemannian} manifolds \changehss{is a generalized concept} of straight lines in \changehss{the} Euclidean space. For \changehss{any tangent} vector in $\xi \in T_w \mathcal{M}$ at $w \in \mathcal{M}$, there exists an interval $I$ about $0$ and a unique geodesic $\gamma_e(t,w,\xi):I \rightarrow \mathcal{M}$, such that $\gamma_e(0)=w$ and $\dot{\gamma_e}(0) = \xi$. \changehss{The {\it exponential mapping} ${\rm Exp}_{w}: T_{w}\mathcal{M} \rightarrow \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, exponential mapping is defined for all vectors $\xi \in T_{w}\mathcal{M}$.
We can thus obtain an update formula using the exponential mapping:
%
\begin{eqnarray*}
	w_{t+1} & = & {\rm Exp}_{w_t}(s_t \xi_{w_t}),	
\end{eqnarray*}
%
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 {\it exponential mapping} \cite[Section~5.4]{Absil_OptAlgMatManifold_2008} that induces a line-search algorithm along geodesics.
Additionally, given two points $w$ and $z$ on $\mathcal{M}$, the \emph{logarithm mapping}, or simply {\it 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$.
\changehs{The log mapping satisfies ${\rm dist} (w,z) = \| {\rm Log}_{w} (z) \|_{w}$,
where ${\rm dist}: \mathcal{M} \times \mathcal{M} \to \mathbb{R}$ is the shortest distance between two points on $\mathcal{M}$.
}

\changehss{The {\it steepest descent} or {\it 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 \changehs{the} {\it Riemannian gradient} \changeHS{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: T_w\mathcal{M} \times T_w\mathcal{M} \rightarrow \mathbb{R}$ over $w \in \mathcal{M}$, the family is called a {\it Riemannian metric} on $\mathcal{M}$. $g_w(\xi_w, \zeta_w)$ is an inner product \changehss{elements} $\xi_w$ and $\zeta_w$ of the tangent space $T_{\changeHS{w}}\mathcal{M}$ at $w$.
\changeHSSS{Herein, 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
%
\begin{eqnarray*}
	 {\rm D} f(w)[\xi_w]& = & \changeHSSS{\langle}{\rm grad} f(w), \xi_w\changeHSSS{\rangle_w}, \ \ \forall \xi_w \in T_w \mathcal{M},
\end{eqnarray*}
where ${\rm D} f(w)[\xi_w]$ is the {\it Fr\'{e}chet derivative} of $f(w)$ in the direction $\xi_w$.

\changehss{In searching a 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 even some Riemannian manifolds for which a closed form for 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 is the same as those of geodesics. A more general update formula is then written as}
%
\begin{eqnarray*}
	w_{t+1} & = & R_{w_t}(s_t \xi_{w_t}),
\end{eqnarray*}
%
where $R_{w_t}$ is a {\it retraction}, which is any map $R_{w}: T_{w}\mathcal{M} \rightarrow \mathcal{M}$ that locally approximates the exponential mapping, up to first-order, on the manifold 
\cite[Definition~4.1.1]{Absil_OptAlgMatManifold_2008}. \changehss{Retractions include the exponential mapping as a special case. An advantage of using retractions is that the computational cost can be reduced compared with using the exponential mapping. It is worth noting that convergence properties for the exponential mapping usually holds for retractions as well.}

In the R-SVRG proposed in Section \ref{Sec:R-SVRG},
we need to add tangent vectors that are in different tangent spaces, say $\tilde{w}$ and $w$ on $\mathcal{M}$.
A mathematically natural way \changehs{to do so} is to use the parallel translation operator.
\changeHKK{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$ \cite[Section 5.4]{Absil_OptAlgMatManifold_2008}, 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 {\it vector transport} on $\mathcal{M}$, which is a map $\mathcal{T}:T\mathcal{M} \oplus T\mathcal{M} \to T \mathcal{M}$ \cite[Definition~8.1.1]{Absil_OptAlgMatManifold_2008}, is used as an alternative.
A vector transport $\mathcal{T}$ has an associated retraction $R$.
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)$, which can be regarded as a transported vector $\xi$ along $\eta$.
Parallel translation is an example of vector transport.
\changeHSSS{
In the following, we use the notations $P_{\eta}$, $P_{\gamma}^{b \leftarrow a}$, 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-a)/(b-a))\eta)$ with $z = R_w(\eta)$.
}

%
%
%
\section{Riemannian stochastic variance reduced gradient}
\label{Sec:R-SVRG}

After a brief explanation of the variance reduced gradient variants in Euclidean space, we describe the proposed Riemannian stochastic variance reduced gradient on \changeHKK{Riemannian} \changeHSSS{manifolds}.  

\subsection{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 {\it stochastic gradient} \changeBM{and $\alpha$ is the step size}. SGD assumes an {\it unbiased estimator} of the full gradient, as $\mathbb{E}_n[\nabla f_n(w_t)] = \nabla f (w_t)$. Many \changeBM{recent} variants of the variance reduced gradient of SGD attempt to reduce its variance $\mathbb{E}[\| v_t - \nabla f(w_t)\|^2]$ \changeBM{as} $t$ increases to achieve \changeBM{better} convergence \cite{Johnson_NIPS_2013_s, Mairal_SIAMJOPT_2015,Roux_NIPS_2012_s,Shalev_arXiv_2012_s,Shalev_JMLR_2013_s,Defazio_NIPS_2014_s,Zhang_SIAMJO_2014_s}. \changeBM{SVRG, proposed in \cite{Johnson_NIPS_2013_s},} introduces an explicit variance reduction strategy with double loops, where \changehs{the} $s$-th outer loop, called \changehs{the} $s$-th {\it 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,\ldots,m_{s-1}\}$ at the end of \changehs{the} $(s\!\!-\!\!1)$-th epoch, and also sets the initial value of \changehs{the} $s$-th epoch \changehs{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, \ldots, N\}$ by $i_t^s$, SVRG randomly picks \changehs{the} $i_t^s$-th sample for each $t \geq 1$ at $s \geq 1$ and computes the {\it modified stochastic gradient} $v_t^s$ as
%
\begin{eqnarray}
\label{Eq:E-SVRG}
v_t^s & = & \nabla f_{i_t^s} (w_{t-1}^{s}) -  \nabla f_{i_t^s} (\tilde{w}^{s-1})  + \nabla f(\tilde{w}^{s-1}). 
\end{eqnarray}
It should be noted that SVRG can be regarded as one special case of S2GD (semi-stochastic gradient descent), which differs in the number of inner loop iterations chosen \cite{Konecny_arXiv_2013}. 


%
\subsection{Proposed Riemannian extension of SVRG (R-SVRG)}
%
We propose a Riemannian extension of SVRG on a Riemannian manifold $\mathcal{M}$, \changehs{called R-SVRG}. Here, we denote the Riemannian stochastic gradient for \changehs{the} $i_t^s$-th sample as $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$ and the {\it modified Riemannian stochastic gradient} as $\xi_t^s$ instead of $v_t^s$, to show the differences from the Euclidean case.


R-SVRG reduces variance 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 \changehs{the} $(s\!\!-\!\!1)$-th epoch, and computes the full Riemannian gradient $\gradf (\tilde{{w}}^{s-1})=\frac{1}{N} \sum_{i=1}^{N} \gradf_i(\tilde{{w}}^{s-1})$ only for this stored $\tilde{{w}}^{s-1}$. The algorithm also computes the Riemannian stochastic gradient $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$ that corresponds to this $i_t^s$-th sample. Then, \changeHK{picking \changehs{the} $i_t^s$-th sample for each $t$-th inner iteration of \changehs{the} $s$-th epoch} at ${w}_{t-1}^{s}$, we calculate $\xi_t^s$ \changeBM{in} the same \changeBM{way} as $v_t^s$ in (\ref{Eq:E-SVRG}), i.e., by modifying the stochastic gradient $\gradf_{i_t^s}({w}_{t-1}^{s})$ using both $\gradf (\tilde{{w}}^{s-1})$ and $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$. \changeBM{Translating the right-hand side of (\ref{Eq:E-SVRG}) to manifold $\mathcal{M}$} involves the sum of \changeBM{$\gradf_{i_t^s}({w}_{t-1}^{s})$, $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$, and $\gradf(\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}$.} \changeBM{This operation requires particular attention on a manifold, and a vector transport \changehss{enables us to handle} multiple elements on two separated tangent spaces \changehss{flexibly}.} More concretely, $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$ and $\gradf(\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 $\gradf_{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 \changehs{the} $t$-th inner iteration of \changehs{the} $s$-th epoch is set \changehs{to}
%
\begin{eqnarray}
\label{Eq:R-SVRG-Grad-paralleltrans}
\xi_t^s & = & \gradf_{i_t^s}({w}_{t-1}^{s}) -  \mathcal{T}^{\scriptsize {w}_{t-1}^{s} \leftarrow \tilde{{w}}^{s-1}} \left(\gradf_{i_t^s}(\tilde{{w}}^{s-1}) - \gradf(\tilde{{w}}^{s-1})\right), 
\end{eqnarray}
%
where $\mathcal{T}_{\gamma}^{\scriptsize {w}_{t-1}^{s} \leftarrow \tilde{{w}}^{s-1}}(\cdot)$ represents a vector transport operator from $\tilde{{w}}^{s-1}$ to ${w}_{t-1}^{s}$ on $\mathcal{M}$.
Specifically, we need to calculate the tangent vector from $\tilde{{w}}^{s-1}$ to ${w}_{t-1}^{s}$,
which is given by \changeBM{the} inverse of the retraction\changeHSSS{, i.e., $R^{-1}$}.
Consequently, the final update rule of R-SVRG is defined as ${w}_{t}^{s} =  R_{\scriptsize {w}_{t-1}^{s}}( - \changeHSSS{\alpha_{t-1}^s} \xi_t^s)$\changeHSSS{, where $\alpha_{t-1}^s > 0$ is a step size at $w_{t-1}^s$}.
\changeHSSS{As shown in our local convergence analysis (Section \ref{AppenSec:LocalConvergenceAnalysis}), $\alpha_{t-1}^s$ can be fixed once the iterate approaches sufficiently close to a solution.}
It should be noted that the modified direction $\xi_t^s$ is also a Riemannian stochastic gradient of $f$ at ${w}_{t-1}^{s}$. 

\changehss{
Let $\mathbb{E}_{i_t^s}[\cdot]$ be the expectation with respect to the random choice of $i_t^s$, conditioned on all 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
\begin{eqnarray*}
\mathbb{E}_{i_t^s}[\xi_t^s] & = & 
\mathbb{E}_{i_t^s}[\gradf_{i_t^s}({w}_{t-1}^{s})] - 
\mathcal{T}^{\scriptsize {w}_{t-1}^{s} \leftarrow \tilde{{w}}^{s-1}} 
\left(
\mathbb{E}_{i_t^s}[\gradf_{i_t^s}(\tilde{{w}}^{s-1}) ] -
\gradf(\tilde{{w}}^{s-1})
\right) \nonumber\\
& = & \gradf({w}_{t-1}^{s}) -
\mathcal{T}^{\scriptsize {w}_{t-1}^{s} \leftarrow \tilde{{w}}^{s-1}} 
\left(
\gradf(\tilde{{w}}^{s-1})  -
\gradf(\tilde{{w}}^{s-1})
\right)  \nonumber\\
& = & \gradf({w}_{t-1}^{s}).
\end{eqnarray*}

The theoretical \changehs{convergence} analysis of the Euclidean SVRG algorithm assumes that the \changehs{initial} vector ${w}_0^s$ of \changehs{the} $s$-th epoch is set to the average (or a random) vector of the $(s\!\!-\!\!1)$-th epoch \cite[Figure~1]{Johnson_NIPS_2013_s}. 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 \ref{Thm:LocalConvergence}, this study also uses, as {\bf option I}, the mean value of $\tilde{{w}}^{s} = g_{m_s}({w}_1^s, \ldots {w}_{m_s}^s)$ as $\tilde{{w}}^s$, where $g_n({w}_1, \ldots, {w}_n)$ is the Karcher mean on the manifold. Alternatively, we can also simply choose $\tilde{{w}}^{s}={w}_t^s$ for $t \in \{1, \ldots, m_s\}$ at random. In addition, as {\bf option II}, we can use the last vector in the $(s\!\!-\!\!1)$-th epoch, i.e., $\tilde{{w}}^{s}={w}_{m_s}^s$.
\changeHSSS{In the global convergence analysis in Section \ref{Sec:GlobalAnalysis}, we use {\bf option II}.}
%by taking the mean to be the unique minimizer of the sum of the squares of the distances defined as $g_n({w}_1, \ldots, {w}_n)= {\rm arg\ min}_{\scriptsize {w} \in {\rm Gr}} \sum_{i}^{n}({\rm dist}({w}_i, {w}))^2$.
%
%
The overall algorithm, with a fixed step size, is summarized in Algorithm \ref{Alg:R-SVRG}. 


%\begin{algorithm}
%\caption{Algorithm for R-SVRG with a fixed step size.}
%\label{Alg:R-SVRG}
%\begin{algorithmic}[1]
%\REQUIRE{Update frequency $m_s>0$ and step size $\eta>0$.}
%\STATE{Initialize $\tilde{{w}}^0$.}
%\FOR{$s=1,2, \ldots$} 
%\STATE{Calculate the Riemannian full gradient $\gradf(\tilde{{w}}^{s-1})$}.
%\STATE{Store ${w}_0^s = \tilde{{w}}^{s-1}$.}
%	\FOR{$t=1,2, \ldots, m_s$} 
%	\STATE{Choose $i_t^s \in \{1, \ldots, N\}$ uniformly at random.}
%	\STATE{Calculate the tangent vector $\zeta$ from $\tilde{{w}}^{s-1}$ to ${w}_{t-1}^{s}$ by logarithm mapping in (\ref{Eq:logarithm_map}).}
%	\STATE{Calculate the modified Riemannian stochastic gradient $\xi_t^s$ \changeHK{in (\ref{Eq:R-SVRG-Grad-paralleltrans})} by parallel-translating $\gradf(\tilde{{w}}^{s-1})$ and $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$ \changeHS{along} $\zeta$ \changeHK{in (\ref{Eq:parallel_translation})} as\\
%	$\xi_t^s = \gradf_{i_t^s}({w}_{t-1}^{s}) -  P_{\gamma}^{{\scriptsize {w}}_{t-1}^{s} \leftarrow \tilde{{\scriptsize {w}}}^{s-1}} \left(\gradf_{i_t^s}(\tilde{{w}}^{s-1})  -  \gradf(\tilde{{w}}^{s-1})\right)$.}
%	\STATE{Update ${w}_t^s$ from ${w}_{t-1}^s$ as ${w}_t^s = {\rm Exp}_{\scriptsize {w}_{t-1}^s}\left(- \eta \xi_t^s \right)$ \changeHK{with the exponential mapping (\ref{Eq:exponential_map})}.}	
%	\ENDFOR
%	\STATE{{\bf option I}: $\tilde{{w}}^{s}=g_{m_s}({w}_1^s,\ldots,{w}_{m_s}^s)$ (or $\tilde{{w}}^{s}={w}_t^s$ for randomly chosen $t \in \{1, \ldots, m_s\}$).}	
%	\STATE{{\bf option II}: $\tilde{{w}}^{s}={w}^s_{m_s}$.}
%\ENDFOR
%\end{algorithmic}
%\end{algorithm}

\begin{algorithm}
\caption{\changeHSSS{R-SVRG for Problem \eqref{Prob:general}.}}
\label{Alg:R-SVRG}
\begin{algorithmic}[1]
\REQUIRE{Update frequency $m_s>0$ \changeHSSS{and sequence $\{\alpha_t^s\}$ of positive step sizes}.}
\STATE{Initialize $\tilde{{w}}^0$.}
\FOR{$s=1,2, \ldots$} 
\STATE{Calculate the full Riemannian gradient $\gradf(\tilde{{w}}^{s-1})$}.
\STATE{Store ${w}_0^s = \tilde{{w}}^{s-1}$.}
	\FOR{$t=1,2, \ldots, m_s$} 
	\STATE{Choose $i_t^s \in \{1, \ldots, N\}$ uniformly at random.}
	\STATE{Calculate the tangent vector $\zeta$ from $\tilde{{w}}^{s-1}$ to ${w}_{t-1}^{s}$ \changeHSS{by $\zeta = R^{-1}_{\tilde{{w}}^{s-1}}({w}_{t-1}^s)$}.}
	\STATE{Calculate the modified Riemannian stochastic gradient $\xi_t^s$ by transporting $\gradf(\tilde{{w}}^{s-1})$ and $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$ \changeHS{along} $\zeta$ as\\
	\changeHSS{$\xi_t^s = \gradf_{i_t^s}({w}_{t-1}^{s}) -  \mathcal{T}_{\zeta} \left(\gradf_{i_t^s}(\tilde{{w}}^{s-1})  -  \gradf(\tilde{{w}}^{s-1})\right)$}.}
	\STATE{Update ${w}_t^s$ from ${w}_{t-1}^s$ as \changeHSS{${w}_t^s = R_{\scriptsize {w}_{t-1}^s}\left(- \changeHSSS{\alpha_{t-1}^s} \xi_t^s \right)$} \changeHSS{with retraction $R$}.}	
	\ENDFOR
	\STATE{{\bf option I}: $\tilde{{w}}^{s}=g_{m_s}({w}_1^s,\ldots,{w}_{m_s}^s)$ (or $\tilde{{w}}^{s}={w}_t^s$ for randomly chosen $t \in \{1, \ldots, m_s\}$).}	
	\STATE{{\bf option II}: $\tilde{{w}}^{s}={w}^s_{m_s}$.}
\ENDFOR
\end{algorithmic}
\end{algorithm}

%
Additionally, variants of the variance reduced SGD \changehs{initially require a} full gradient calculation at every epoch.
This initially results in more overhead than the ordinal SGD algorithm and eventually causes a {\it cold-start} property \changehs{in these variants}. \changeBM{To avoid this, \cite{Konecny_arXiv_2013}, in Euclidean space, proposes using the standard SGD update only for \changehs{the} first epoch.}
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 modified stochastic gradient $\xi_t^s$. The $s$-th epoch, i.e., the outer loop, requires $N+2m_s$ gradient evaluations, where $N$ is for the full gradient $\gradf(\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., $\gradf_{i_t^s}({w}_{t-1}^{s})$ and $\gradf_{i_t^s}(\tilde{{w}}^{s-1})$. However, if $\gradf_{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 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 is suggested in \cite{Johnson_NIPS_2013_s}).

\section{\changeHSSS{Global convergence analysis}}
\label{Sec:GlobalAnalysis}
In this section, we provide a global \changeHSSS{convergence analysis} of \changeHSSS{Algorithm \ref{Alg:R-SVRG} for Problem \eqref{Prob:general}} after introducing some assumptions.
Throughout this section, we let $R$ and $\mathcal{T}$ denote a retraction and vector transport used in Algorithm \ref{Alg:R-SVRG}, respectively,
and suppose the following assumptions.
\begin{assumption}
\label{assump:T_isometry}
The vector transport $\mathcal{T}$ is isometric on $\mathcal{M}$, i.e., for any $w \in \mathcal{M}$ and $\xi_w, \eta_w \in T_w \mathcal{M}$, $\|\mathcal{T}_{\eta_w}(\xi_w)\|_{R_w(\eta_w)} = \|\xi_w\|_w$.
\end{assumption}
Note that we can construct an isometric vector transport as in \cite{huang2015broyden} so that Assumption \ref{assump:T_isometry} holds.

\begin{assumption}
\label{Appen_assump:C2}
The objective function \changeHSSS{$f$} \changeHSSS{and its components $f_1, f_2, \dots, f_N$} are twice continuously differentiable.
\end{assumption}

\changeHSSS{
As discussed in \cite{huang2015broyden}, 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)$, \changeHSS{which is the ball in $T_w \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.
\begin{assumption}
\label{assump:GlobalTRN}
For a sequence $\{w_t^s\}$ generated by Algorithm \ref{Alg:R-SVRG}, there exists a compact and connected set $K \subset \mathcal{M}$ such that $w_t^s \in K$ for all $s, t \ge 0$.
Additionally, for each $s\ge 1$, there exists a totally retractive neighborhood $\Omega_{s-1}$ of $\tilde{w}^{s-1}$ such that $w_t^s$ stays in $\Omega_{s-1}$ for any $t\ge 0$.
Furthermore, suppose that there exists $I>0$ such that $\inf_{s \ge 1} \{\sup_{z \in \Omega_{s-1}} \| R_{\tilde{w}^{s-1}}^{-1}(z)\|_{\tilde{w}^{s-1}}\} \ge I$.
\end{assumption}
}

\changeHSSS{In the global convergence analysis,} we also assume that the sequence of step sizes $(\alpha_t^s)_{t \geq 1, s \geq 1}$ satisfies the usual condition in stochastic approximation as follows:
\begin{assumption}
\label{assump:stepsizes}
The sequence $\{\changeHSSS{\alpha}_t^s\}$ of step sizes satisfies
\begin{eqnarray}
\label{Eq:StepsizeCondition}
\sum (\alpha_t^s)^2 \ <\  \infty {\ \ \ \rm and\ \ \ } \sum \alpha_t^s \ =\  \changeHSSS{\infty}.
\end{eqnarray}
\end{assumption}
\changeHSSS{This condition is satisfied, for example, if $\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$.}
We also note the following proposition \changeHSSS{introduced in \cite{Fisk_1965_s}}:
%
\begin{proposition}[\cite{Fisk_1965_s}]
\label{Prop:FiskPro}
Let $(X_n)_{n \in \mathbb{N}}$ be a nonnegative stochastic process with bounded positive variations, i.e., such that $\sum_{\changeHSSS{n=}0}^{\infty} \changeHSSS{\mathbb{E}[\mathbb{E}[X_{n+1} - X_n | \mathcal{F}_n]^{+}]} < \infty$, where \changeHSSS{$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., 
\begin{eqnarray*}
\sum_{n=0}^{\infty} |\mathbb{E}[ X_{n+1} - X_n| \mathcal{F}_n] | < \infty \ \ \ {\rm a.s.},\ {\rm and\ } X_n\ {\rm converges\ } {\rm a.s}. 
\end{eqnarray*}
\end{proposition}


Now, we give the a.s. convergence of the proposed algorithm under \changehss{the assumption that the generated sequenc is in a compact set}.

\begin{theorem}
\label{Thm:GlobalConvergence}
\changeHSSS{Suppose Assumptions \ref{assump:T_isometry}--\ref{assump:GlobalTRN} and} consider Algorithm \ref{Alg:R-SVRG} \changeHSSS{with {\bf option II} and} step sizes \changeHSSS{$\{\alpha_t^s\}$} satisfying Assumption \ref{assump:stepsizes} on a Riemannian manifold $\mathcal{M}$.
\changeHSSS{If $f\ge 0$, }then $\{f(w_t^s)\}$ converges a.s. and $\gradf (w_t^s) \rightarrow 0$ a.s.
\end{theorem}

\begin{proof}
The claim is proved similarly to the proof of the standard Riemannian SGD (see \cite{Bonnabel_IEEETAC_2013_s}).
Since $\changeHSSS{K}$ is compact, all continuous functions on $\changeHSSS{K}$ can be bounded.
Therefore, there exists $\changeHSSS{C} > 0$ such that for all $w \in \changeHSSS{K}$ and $n \in \{1,2,\ldots,N\}$, we have $\| \gradf (w) \|_{\changeHSSS{w}} \leq \changeHSSS{C}$ \changeHS{and $\| \gradf_{n}(w)\|_{\changeHSSS{w}} \le \changeHSSS{C}$}.
\changeHSSS{We use $C':= 3C$ in the following.}
Moreover, as $\alpha_t^s \rightarrow 0$, there exists $\changeHSSS{s_0}$ such that for $\changeHSSS{s} \geq \changeHSSS{s_0}$ we have $\alpha_t^s \changeHSSS{C'} < I$.
Suppose now that $\changeHSSS{s} \geq \changeHSSS{s_0}$.
\changeHSSS{
Let $\tilde{\eta}_{t}^s$ satisfy $R_{\tilde{w}^{s-1}}(\tilde{\eta}_{t}^s) = w_{t}^s$, i.e., $\tilde{\eta}_{t}^s=R^{-1}_{\tilde{w}^{s-1}}(w_{t}^s)$.
}
It follows from the triangle inequality that
\changeHSSS{
\begin{eqnarray*}
\|\xi_{t+1}^s \|_{w_{t}^s} & = & \| \gradf_{i_{t+1}^s}(w_{t}^{s}) -  \changeHSS{\mathcal{T}_{\tilde{\eta}_{t}^s}} \left(\gradf_{i_{t+1}^s}(\tilde{w}^{s-1}) \right) + \changeHSS{\mathcal{T}_{\tilde{\eta}_{t}^s}} \left(\gradf(\tilde{w}^{s-1})\right) \|_{w_{t}^s} \\
& \le & \changeHSSS{C'}/3 + \changeHSSS{C'}/3 + \changeHSSS{C'}/3 = \changeHSSS{C'},
\end{eqnarray*}
where we have used the assumption that $\mathcal{T}$ is an isometry.
Therefore, we have
\begin{equation*}
\|R_{w_t^s}^{-1}(w_{t+1}^s)\|_{w_t^s} = \| -\alpha_t^s \xi_{t+1}^s \|_{w_t^s} \le \alpha_t^s \changeHSSS{C'} < I.
\end{equation*}
Hence, it follows from the assumption that
}
there exists a curve $R_{w_t^s}(-\changeHSSS{\tau} \alpha_t^s \xi_{t+1}^s)_{0 \leq \changeHSSS{\tau} \leq 1}$ linking $w_t^s$ and $w_{t+1}^s$.
Thus, the Taylor formula implies that 
\changeHSS{
\begin{eqnarray*}
& & f(w_{t+1}^s) - f(w_{t}^s) \changeHSSS{=f(\changeHSS{R_{w_t^s}}(-\alpha_t^s \xi_{t+1}^s)) - f(w_t^s)} \\
& = & - \alpha_t^s 
\langle \xi_{t+1}^s, \gradf (w_t^s) 
\rangle_{w_t^s} + (\alpha_t^s)^2 \int_0^1 (1-\changeHSSS{\tau}) \langle \nabla^2 (f \circ R_{w_{t}^s}) (-\changeHSSS{\tau}\alpha_t^s \xi_{t+1}^s)[\xi_{t+1}^s], \xi_{t+1}^s\rangle_{w_t^s} \changeHSSS{d\tau}\\
& \leq & - \alpha_t^s 
\langle \xi_{t+1}^s, \gradf (w_t^s) 
\rangle_{w_t^s} + (\alpha_t^s)^2 
\| \xi_{t+1}^s \|_{w_t^s}^2 k_1,
\end{eqnarray*}
where $2k_1$ is an upper bound of the largest eigenvalues of the Euclidean Hessian $\nabla^2 (f \circ R_{w_t^s})$ of $f \circ R_{w_t^s}$ in the compact ball $\{ \xi \in T_{w_t^s} \mathcal{M} \,|\, \| \xi \|_{w_t^s} \le I\}$.}
Let $\mathcal{F}_t^s$ be the increasing sequence of $\sigma$-algebras \changehss{defined by}
%
\begin{eqnarray*}
\mathcal{F}_t^s=\{i_1^1, \ldots, i_{m_1}^1, \ldots, i_1^{s-1}, \ldots, i_{m_{s-1}}^{s-1}, i_1^{s}, \ldots, i_{t-1}^{s}\}.
\end{eqnarray*}
Since $w_t^s$ is computed from $i_1^1\changeHSSS{,} \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 
\begin{eqnarray*}
& & \mathbb{E}[\langle \xi_{t+1}^s, \gradf (w_t^s)\rangle_{\changeHSSS{w_t^s}} | \mathcal{F}_{t+1}^s] \\
& = & \mathbb{E}_{i_{t+1}^s}[\langle \xi_{t+1}^s, \gradf (w_t^s)\rangle_{\changeHSSS{w_t^s}}] \\
& = & \changeHSSS{\langle \mathbb{E}_{i_{t+1}^s}[\gradf_{i_{t+1}^s}(w_{t}^{s})], \gradf(w_t^s) \rangle_{\changeHSSS{w_t^s}}}
-\changeHSSS{\langle \changeHSSS{\mathcal{T}_{\tilde{\eta}_{\changeHSSS{t}}^s}}(\mathbb{E}_{\changeHSSS{i_{t+1}^s}}[\gradf_{i_{t+1}^s}(\tilde{w}^{s-1})]), \gradf(w_t^s) \rangle_{\changeHSSS{w_t^s}}}\\
&&
\changeHSSS{+\langle  \mathcal{T}_{\tilde{\eta}_{\changeHSSS{t}}^s}(\gradf(\tilde{w}^{s-1})), \gradf(w_t^s) \rangle_{\changeHSSS{w_t^s}}} \\
%
& = & \| \gradf (w_t^s)\|_{\changeHSSS{w_t^s}}^2,
\end{eqnarray*}
%
which yields that
\begin{eqnarray}
\label{Eq:ExpectationDiff}
\mathbb{E}[f(w_{t+1}^s) - f(w_{t}^s) | \mathcal{F}_{t+1}^s] & \leq & -\alpha_t^s \| \gradf (w_t^s)\|^2 + (\alpha_t^s)^2 \changeHSSS{C'}^2 k_1 \changeHSSS{\le (\alpha_t^s)^2 \changeHSSS{C'}^2 k_1},
\end{eqnarray}
as $\| \xi_{t+1}^s \|_{\changeHSSS{w_t^s}} \leq \changeHSSS{C'}$.

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_{\changeHSSS{k}}) \geq 0$, \changeHSSS{\eqref{Eq:ExpectationDiff}} proves that $\{f(w_{\changeHSSS{k}}) +\sum_{\changeHSSS{l=k}}^{\infty} (\alpha_{\changeHSSS{l}})^2 \changeHSSS{C'}^2 k_1\}$ is a nonnegative supermartingale. Hence, it converges a.s.
\changeHSSS{Returning to the original indexing, this implies} that $\changeHSSS{\{}f(w_{t}^s)\changeHSSS{\}}$ converges a.s. Moreover, summing the inequalities, we have
\begin{eqnarray}
\label{Eq:ExpectationDiff2}
\sum_{\changeHSSS{s \geq s_0},t} \alpha_t^s \| \gradf (w_t^s) \|^2 & \leq & - \sum_{\changeHSSS{s \geq s_0},t} \mathbb{E}[f(w_{t+1}^s) - f(w_{t}^s) | \mathcal{F}_{\changeHSSS{t+1}}^s] + \sum_{\changeHSSS{s \geq s_0},t}  (\alpha_t^s)^2 \changeHSSS{C'}^2 k_1.
\end{eqnarray}
Here, we prove that the right term is bounded\changehss{. By doing so, we can conclude the convergence of the left term}.

Summing (\ref{Eq:ExpectationDiff}) over $s$ and $t$, it is clear that $f(w_t^s)$ satisfies the assumption of Proposition \ref{Prop:FiskPro}; thus, $\{f(w_t^s)\}$ is a quasi-martingale, implying \changehs{that} $\sum_{t \geq t_0} \alpha_t^s \| \gradf (w_t^s) \|^2$ converges a.s. because of inequality (\ref{Eq:ExpectationDiff2}), where
$\left| - \sum_{\changeHSSS{s \geq s_0},t} \mathbb{E}[f(w_{t+1}^s) - f(w_{t}^s) | \mathcal{F}_{\changeHSSS{t+1}}^s]\right| < \infty$, by Proposition \ref{Prop:FiskPro}.
If $\{\|  \gradf (w_t^s) \|\}$ is proved to converge a.s., it can only converge to 0 a.s. because of condition (\ref{Eq:StepsizeCondition}).

Now consider the nonnegative process $p_t^s = \|  \gradf (w_t^s) \|^2$.
Bounding the second derivative of $\|  \gradf \|^2$ by $k_2$, along the \changeHSSS{curve defined by the retraction $R$} linking $w_t^s$ and $w_{t+1}^s$, a Taylor expansion yields 
\begin{eqnarray*}
p_{t+1}^s - p_t^s \leq -2 \alpha_t^s \langle \gradf (w_t^s), {\rm Hess} f(w_t^s)  [\xi_{t+1}^s] \rangle_{\changeHSSS{w_t^s}} + (\alpha_t^s)^2 \|  \xi_{t+1}^s \|_{\changeHSSS{w_t^s}}^2 k_2.
\end{eqnarray*}
\changehss{Let $k_3$ be a constant such that $-k_3$ is a lower bound of the Hessian of $f$ on the compact set. Then, we have}
\begin{eqnarray*}
\mathbb{E}(p_{t+1}^s - p_t^s | \mathcal{F}_{t+1}^s) \leq 2 \alpha_t^s \|  \gradf (w_t^s) \|_{\changeHSSS{w_t^s}}^2 k_3 + (\alpha_t^s)^2 \changeHSSS{C'}^2 k_2.
\end{eqnarray*}
\changehss{We have proved the fact that the sum of the right term is finite},
implying that $\changeHSSS{\{}p_t^s\changeHSSS{\}}$ is a quasi-martingale, thus further implying the a.s. convergence of $\changeHSSS{\{p_t^s\}}$ \changehss{to $0$}, \changeHSSS{as claimed}.
\end{proof}

\changeHSSS{
Theorem \ref{Thm:GlobalConvergence} contains 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 notions of injectivity radius.
}

\begin{corollary}
\label{Thm:GlobalConvergenceExp}
\changeHSSS{
Suppose Assumption \ref{assump:T_isometry} and consider Algorithm \ref{Alg:R-SVRG} with {\bf option II} and step sizes $\{\alpha_t^s\}$ satisfying Assumption \ref{assump:stepsizes} on a Riemannian manifold $\mathcal{M}$, where the exponential mapping and parallel translation are used as retraction and vector transport, respectively.
Assume that $\mathcal{M}$ is connected and has injectivity radius uniformly bounded from below by $I > 0$.
Assume also that there exists a compact set $K \subset \mathcal{M}$ such that $w_t^s \in K$ for all $s, t \ge 0$.
If $f\ge 0$, then $\{f(w_t^s)\}$ converges a.s. and $\gradf (w_t^s) \rightarrow 0$ a.s.
}
\end{corollary} 






































































%
\section{\changeHSSS{Local convergence rate analysis}}
\label{Sec:LocalAnalysis}
\changeHSSS{
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 \ref{Alg:R-SVRG} 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 convex in such a neighborhood.
We first give formal expressions of these assumptions and then analyze the local convergence rate of the algorithm.
}

\subsection{Assumptions and existing lemmas}
\changeHSSS{
We again suppose Assumptions \ref{assump:T_isometry} and \ref{Appen_assump:C2}.
That is, the objective function is $C^2$ and the vector transport is isometric.
}

\changeHSSS{
Let $w^*$ be a critical point.
As discussed in \cite{huang2015broyden}, we note that there exists a positive constant $\rho$ and $\rho$-totally retractive neighborhood $\Omega$ of $w^* \in \mathcal{M}$.
In the local convergence analysis, we assume the following.
\begin{assumption}
\label{assump:Omega}
The sequence $\{w_t^s\}$ generated by Algorithm \ref{Alg:R-SVRG} continuously remains in totally retractive neighborhood $\Omega$ of critical point $w^*$, i.e., $R_{w_t^s}(-\alpha\changeHSSS{\xi_{t+1}^s}) \in \Omega$ for all $s, t \ge 0$ and for all $\alpha\in [0,\alpha_t^s]$.
In addition, the radius of $\Omega$ is sufficiently small.
\end{assumption}
}

\begin{assumption}
\label{assump:T_Lipschitz}
There exists a constant $c_0$ such that vector transport $\mathcal{T}$ satisfies the following conditions:
\begin{equation*}
\| \mathcal{T}_\eta-\mathcal{T}_{R_\eta}\| \le c_0 \| \eta\|, \qquad \| \mathcal{T}^{-1}_\eta-\mathcal{T}^{-1}_{R_\eta}\| \le c_0 \| \eta\|,
\end{equation*}
where $\mathcal{T}_R$ denotes the differentiated retraction, i.e.,
\begin{equation*}
\mathcal{T}_{R_{\eta_w}}(\xi_w) = {\rm D} R(\eta_w) [\xi_w], \qquad \eta_w, \xi_w \in T_w \mathcal{M}, \quad w \in \mathcal{M}.
\end{equation*}
\end{assumption}
Note that it follows from Taylor's theorem and the fact \changehs{that} $\mathcal{T}_{0} = \mathcal{T}_{R_0}$ that Assumption \ref{assump:T_Lipschitz} holds if $\mathcal{T}$ and $\mathcal{T}_R$ are $C^1$.

Also, we assume the following.
\begin{assumption}
\label{Assump:ret_convex}
$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 \le \frac{d^2}{d\alpha^2} f(R_w(\alpha\eta)) \le a_1$ for all $w \in \Omega$, $\eta \in T_w \mathcal{M}$ with $\| \eta \|_{\changeHSSS{w}} = 1$, and for all $\alpha$ \changeHSSS{satisfying} $R_w(\tau \eta) \in \Omega$ for all $\tau \in [0,\alpha]$.
\changeHSSS{
Furthermore, $f$ is strongly convex in $\Omega$, i.e., the above claim is true even if $R$ is replaced with the exponential mapping ${\rm Exp}$.
}
\end{assumption}
Letting $\Omega$ be smaller if necessary, we can guarantee the last assumption from other assumptions using Lemma 3.1 in \cite{huang2015broyden}.


\changeHSSS{
In the rest of this section, we introduce some existing lemmas to evaluate the differences of using retraction and vector transport instead of exponential mapping and parallel translation, and the effects of the curvature of the manifold in question.
}
\begin{lemma}[In the proof of Lemma 3.9 in \cite{huang2015broyden}] \label{Lem:pseudo_Lipschitz}
Under Assumptions \ref{Appen_assump:C2} and \ref{assump:Omega}, there exists a constant $\beta>0$ such that
\begin{eqnarray}
     \label{pseudo_Lipschitz}
	\| P_{\changeHSSS{\gamma}}^{w \leftarrow z} ({\rm grad} f(z)) - {\rm grad} f(w) \|_{\changeHSSS{w}} & \leq & \beta {\rm dist}(z,w),
\end{eqnarray}
where $w$ and $z$ are in $\Omega$ and \changeHSSS{$\gamma$} is a curve $\gamma(\tau):=R_z(\tau\eta)$ for \changeHSSS{an arbitrary} $\eta \in T_z \mathcal{M}$ defined by retraction $R$ on $\mathcal{M}$.
$P_{\changeHSSS{\gamma}}^{w \leftarrow z}(\cdot)$ is a parallel translation operator along curve \changeHSSS{$\gamma$} from $z$ to $w$.
\end{lemma}
Note that curve $\gamma$ in this lemma is not necessarily the geodesic on $\mathcal{M}$.
Relation \eqref{pseudo_Lipschitz} is a generalization of the Lipschitz continuity condition.
\begin{lemma}[Lemma 3.5 in \cite{huang2015broyden}] \label{Lem:T_Lipschitz}
\changeHSSS{
Let $\mathcal{T}\in C^{0}$ be a vector transport associated with the same retraction $R$ as that of \changeHSSS{the parallel transport} $P\changeHSSS{\in C^{\infty}}$.
}
Under Assumption \ref{assump:T_Lipschitz}, for any $\bar{w}\in\mathcal{M}$, there exists a constant $\changeHSSS{\theta}>0$ and neighborhood $\mathcal{U}$ of $\bar{w}$ such that for all $w, z \in \mathcal{U}$,
\begin{equation*}
\|\mathcal{T}_\eta(\xi)-P_{\eta}(\xi)\|_{\changeHSSS{z}} \le \changeHSSS{\theta}\| \xi\|_{\changeHSSS{w}} \| \eta\|_{\changeHSSS{w}},
\end{equation*}
where $\xi, \eta \in T_w\mathcal{M}$ and $R_w(\eta)=z$.
\end{lemma}
We can derive the following lemma from the Taylor expansion \changeHSSS{as in the proof of Lemma 3.2 in \cite{huang2015broyden}}.
\begin{lemma}[\changeHSSS{In the proof of Lemma 3.2 in \cite{huang2015broyden}}]
\label{Lem:mean-value}
%{\rm HS: Here is my idea for the statement and proof of the Lemma which should be here.}\\
\changeHSSS{
Under Assumptions \ref{assump:Omega}--\ref{Assump:ret_convex}, there exists a positive real number $\sigma$ such that
\begin{eqnarray}
\label{Eq:mean-value}
f(z) \ge f(w) + \langle {\rm Exp}_{w}^{-1}( z), \gradf(w)\rangle_{w} + \frac{\sigma}{2} \|{\rm Exp}_{w}^{-1}(z)\|_{w}^2, \qquad w, z \in \Omega.
\end{eqnarray}
}
\end{lemma} 
\begin{proof}
\changeHSSS{
From the assumptions, there exists $\sigma>0$ such that $\frac{d^2}{dt^2}f({\rm Exp}_w(\alpha\eta)) \ge \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.
}
\end{proof}
If we replace ${\rm Exp}$ in Lemma \ref{Lem:mean-value} with retraction $R$ on $\mathcal{M}$, we can obtain a similar result for $R$,
\changeHSSS{which is shown in Lemma 3.2 in \cite{huang2015broyden} where the constant $a_0$ corresponds to $\sigma$ in our Lemma \ref{Lem:mean-value}}.
However, Lemma \ref{Lem:mean-value} for ${\rm Exp}$ is sufficient in the following discussion.

From Lemma 3 in \changeHSSS{\cite{huang2015riemannian}}, we have the following lemma:
\begin{lemma}[\changeHSSS{Lemma 3 in \cite{huang2015riemannian}}]
\label{lemma:retraction_dist}
Let $\mathcal{M}$ be a Riemannian manifold endowed with retraction $R$ and let $\bar{w} \in \mathcal{M}$.
Then, there exist $\changeHSSS{\mu} > 0$ and $\delta_{\changeHSSS{\mu}}\changeHSSS{>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 \le \delta_{\changeHSSS{\mu}}$, the inequality
\begin{equation}
\label{Eq:tau}
\|\xi\|_w \le \changeHSSS{\mu} {\rm dist} (w, R_w(\xi))
\end{equation}
holds.
\end{lemma}
\changeHSSS{
Since we have ${\rm dist}(w, {\rm Exp}_w(\xi)) = \| \xi \|_w$, Eq.~\eqref{Eq:tau} is equivalent to ${\rm dist}(w, {\rm Exp}_w(\xi)) \le \changeHSSS{\mu} {\rm dist} (w, R_w(\xi))$, which gives a relation between the exponential mapping and a general retraction.
}

Now, we introduce Lemma 6 in \cite{Zhang_JMLR_2016_s} to evaluate the distance between $w_{t}^s$ and $w^{*}$ using the smoothness of our objective function.
\changeHSSS{In the following, an Alexandrov space is defined as a length space whose curvature is bounded.}
\begin{lemma}[Lemma 6 in \cite{Zhang_JMLR_2016_s}]
\label{LemZhang}
If $a$, $b$, \changehs{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
\begin{equation*}
a^2 \le \frac{\sqrt{|\kappa|}c}{\tanh(\sqrt{|\kappa|}c)}b^2 + c^2 - 2bc\cos(A).
\end{equation*}
\end{lemma}















\subsection{Local convergence rate analysis with retraction and vector transport}
\label{AppenSec:LocalConvergenceAnalysis}

We now demonstrate the local convergence properties of the R-SVRG algorithm (i.e., local convergence to local minimizers) and its convergence rate. 

We first describe an assumption \changeHSSS{that generated points are sufficiently close to a critical point}.

We show a property of the Karcher mean on a general Riemannian manifold.
\begin{lemma}
\label{AppenLem:KarcherMeanDistance}
Let $w_1,\dots,w_m$ be points on a Riemannian manifold $\mathcal{M}$ and let $w$ be the Karcher mean of the $m$ points.
For an arbitrary point $p$ on $\mathcal{M}$, we have
\begin{eqnarray*}
({\rm dist}(p,w))^2\le \frac{4}{\changeHK{m}}\sum_{i=1}^m({\rm dist}(p,w_i))^2.
\end{eqnarray*}
\end{lemma}
\begin{proof}
From the triangle inequality \changeHK{and $(a + b)^2 \leq 2 a^2 + 2 b^2$ for real numbers $a$ and $b$}, we have for $i=1,2,\ldots,m$
\begin{equation*}
({\rm dist}(p,w))^2 \le \left({\rm dist}(p,w_i)+{\rm dist}(w_i,w)\right)^2 \le 2({\rm dist}(p,w_i))^2 + 2({\rm dist}(w_i,w))^2.
\end{equation*}
Since $w$ is the Karcher mean of $w_1,w_2,\ldots,w_m$, it holds that
\begin{equation*}
\sum_{i=1}^{\changeHK{m}} ({\rm dist}(w,w_i))^2 \le \sum_{i=1}^m ({\rm dist}(p,w_i))^2.
\end{equation*}
It then follows that
\begin{equation*}
\changeHK{m\ } \changeHSSS{(}{\rm dist}(p,w)\changeHSSS{)}^2 \le 2\sum_{i=1}^m ({\rm dist}(p,w_i))^2 + 2\sum_{i=1}^m ({\rm dist}(w_i,w))^2 \le 4 \sum_{i=1}^m ({\rm dist}(p,w_i))^2.
\end{equation*}
This completes the proof.
\end{proof}
%
%
%

Then, we prove that the norms of modified stochastic gradients in R-SVRG are sufficiently small under Assumption \ref{assump:Omega}.
\begin{lemma}
Under Assumptions \ref{assump:T_isometry}--\ref{Appen_assump:C2}, \ref{assump:Omega}, and \ref{assump:T_Lipschitz}, \changeHSSS{the norm of} $\xi_t^s$ computed by \eqref{Eq:R-SVRG-Grad-paralleltrans} is sufficiently small; i.e.,
for any $\varepsilon >0$, \changeHSSS{there exists $r_0>0$ such that $\|\xi_t^s\|_{w_{t-1}^s} \le \varepsilon$ for $r < r_0 $, where $r$ is the radius of $\Omega$.}
\end{lemma}

\begin{proof}
\changeHSSS{
We first show that $\|\gradf_n\|$ for any $n \in \{1,2,\dots, N\}$ is upper-bounded in $\Omega$ if $\Omega$ is sufficiently small.
Since $\Omega$ is sufficiently small, $\Omega$ is contained in a set $U \subset \mathcal{M}$ diffeomorphic to an open set $U' \subset \mathbb{R}^{\dim \mathcal{M}}$ by a chart $\phi: U \to U'$.
Consider a sufficiently small closed ball $B$ in $U'$ 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, $\|\gradf_n\|$ is upper bounded in $\Omega$.
In this proof, let $\changeHSSS{C}$ denote an upper bound, i.e., $\|\gradf_n(z)\|_z \le \changeHSSS{C}$ for all $n \in \{1,2,\dots,N\}$ and $z \in \Omega$.
}

Let $\tilde{\eta}_{t-1}^s = R^{-1}_{\tilde{w}^{s-1}}(w_{t-1}^s)$.
The definition of $\xi_t^s$, Eq.~\eqref{Eq:R-SVRG-Grad-paralleltrans}, and the triangle inequality yield
\begin{align*}
\|\xi_t^s\|_{w_{t-1}^s} = & \| \gradf_{i_t^s}(w_{t-1}^s) - \mathcal{T}_{\tilde{\eta}_{t-1}^s}(\gradf_{i_t^s}(\tilde{w}^{s-1})-\gradf(\tilde{w}^{s-1}))\|_{w_{t-1}^s}\\
\le & \| \gradf_{i_t^s}(w_{t-1}^s) - \mathcal{T}_{\tilde{\eta}_{t-1}^s}(\gradf_{i_t^s}(\tilde{w}^{s-1})) \|_{w_{t-1}^s} + \| \mathcal{T}_{\tilde{\eta}_{t-1}^s}(\gradf(\tilde{w}^{s-1}))\|_{w_{t-1}^s} \\
= & \| \gradf_{i_t^s}(w_{t-1}^s) - P_{\gamma_{t-1}^s}^{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) \|_{w_{t-1}^s}\\
& + \|P_{\gamma_{t-1}^s}^{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) - \mathcal{T}_{\tilde{\eta}_{t-1}^s}(\gradf_{i_t^s}(\tilde{w}^{s-1}))\|_{w_{t-1}^s} \\
& + \| \gradf(\tilde{w}^{s-1})\|_{\tilde{w}^{s-1}},
\end{align*}
where $\gamma_{t-1}^s(\tau) = R_{\tilde{w}^{s-1}}(\tau\tilde{\eta}_{t-1}^s)$ and $P_{\gamma_{t-1}^s}$ is the parallel translation operator along curve $\gamma_{t-1}^s$.
The three terms on the right-hand side above are bounded as
\begin{equation*}
\| \gradf_{i_t^s}(w_{t-1}^s) - P_{\gamma_{t-1}^s}^{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) \|_{w_{t-1}^s} \le \beta {\rm dist}(w_{t-1}^s, \tilde{w}^{s-1}) \le 2\beta r,
\end{equation*}
\begin{align*}
& \|P_{\gamma_{t-1}^s}^{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) - \mathcal{T}_{\tilde{\eta}_{t-1}^s}(\gradf_{i_t^s}(\tilde{w}^{s-1}))\|_{w_{t-1}^s} \\
\le & \changeHSSS{\theta}\|\gradf_{i_t^s}\changeHSSS{(\tilde{w}^{s-1})}\changeHSSS{\|_{\tilde{w}^{s-1}}} \| \tilde{\changeHSSS{\xi}}_{t\changeHSSS{-1}}^s\|_{\tilde{w}^{s-1}}
\le \changeHSSS{\theta}\changeHSSS{C}\changeHSSS{\mu}{\rm dist}(w_{t-1}^s, \tilde{w}^{s-1}) \le 2\changeHSSS{\theta}\changeHSSS{C}\changeHSSS{\mu} r,
\end{align*}
\begin{align*}
\| \gradf(\tilde{w}^{s-1})\|_{\tilde{w}^{s-1}}
= &\| \gradf(\tilde{w}^{s-1}) - P_{\tilde{\gamma}_{*}^{s-1}}^{\tilde{w}^{s-1} \leftarrow w^{\changeHSSS{*}}}(\gradf({w}^{\changeHSSS{*}})) \|_{\tilde{w}^{s-1}}\\
 \le & \beta {\rm dist}(\tilde{w}^{s-1}, w^{\changeHSSS{*}}) \le \beta r,
\end{align*}
using Lemmas \ref{Lem:pseudo_Lipschitz}, \ref{Lem:T_Lipschitz}, and \ref{Lem:pseudo_Lipschitz}, respectively, where $\beta$ and $\changeHSSS{\theta}$ are in the lemmas \changeHSSS{and $\tilde{\gamma}_*^{s-1}$ is a curve defined by $\tilde{\gamma}_*^{s-1}(\tau) = R_{w^*}(\tau R_{w^*}^{-1}(\tilde{w}^{s-1}))$}.
Therefore, we have
\begin{equation*}
\|\xi_t^s\|_{w_{t-1}^s} \le r(3\beta+2\changeHSSS{\theta}\changeHSSS{C}\changeHSSS{\mu}) \changeHSSS{< \varepsilon}
\end{equation*}
\changeHSSS{if we choose a sufficiently small $r$ such that $r<\varepsilon/(3\beta+2\changeHSSS{\theta}\changeHSSS{C}\changeHSSS{\mu})$}.
\end{proof}

We now derive the upper bound of the variance of $\xi_t^s$ as follows. 


%
%
%

\begin{lemma}
\label{AppenLem:UpperBoundVariance}
Suppose Assumptions \ref{assump:T_isometry}--\ref{Appen_assump:C2}, \ref{assump:Omega}, and \ref{assump:T_Lipschitz}, which guarantee Lemmas \ref{Lem:pseudo_Lipschitz}, \ref{Lem:T_Lipschitz}, and \ref{lemma:retraction_dist} \changeHSSS{for $\bar{w} = w^*$}.
\changeHSSS{
Let $\beta>0$ be a constant such that
\begin{eqnarray*}
	\| P_{\changeHSSS{\gamma}}^{w \leftarrow z} ({\rm grad} f_n(z)) - {\rm grad} f_n(w) \|_{\changeHSSS{w}} & \leq & \beta {\rm dist}(z,w),\qquad w, z \in \Omega,\ n = 1,2,\dots, N.
\end{eqnarray*}
The existence of such a $\beta$ is guaranteed by Lemma \ref{Lem:pseudo_Lipschitz}.
}
The upper bound of the variance of $\xi_t^s$ is given by
\begin{eqnarray}
\label{Append_Eq:UpperBoundVariance}
	\mathbb{E}_{i_t^s}\left[\| \xi_t^s \|_{\changeHSSS{w_{t-1}^s}}^2\right] &\leq &
4(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)(7({\rm dist}(w_{t-1}^s,w^*))^2 + 4({\rm dist}(\tilde{w}^{s-1},w^*))^2),
\end{eqnarray}
where the constant $\changeHSSS{\theta}$ corresponds to that in Lemma \ref{Lem:T_Lipschitz},
$C$ is the upper bound of $\|\gradf_n(w)\|,\ n=1, 2, \dots, N$ for $w \in \Omega$,
and $\changeHSSS{\mu} > 0$ appears in \eqref{Eq:tau}.
\end{lemma}

\begin{proof}
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$.
\changeHSSS{
By Lemma \ref{Lem:pseudo_Lipschitz} and \ref{Lem:T_Lipschitz},
}
the upper bound of the variance of $\xi_t^s$ in terms of the distance of $w_t^s$ and $\tilde{w}^{s-1}$ from $w^*$, \changehs{is}
\allowdisplaybreaks[1]
\begin{eqnarray*}
& & \mathbb{E}_{i_t^s}[\| \xi_t^s \|_{\changeHSSS{w_{t-1}^s}}^2]   \nonumber\\
&=&\mathbb{E}_{i_t^s}\left[\| \bigl( \gradf_{i_t^s}(w_{t-1}^{s}) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) \bigr)  \right. \nonumber\\
&& \left. + \bigl( \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) -  \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) + \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf(\tilde{w}^{s-1}) \right)\bigr) \|_{\changeHSSS{w_{t-1}^s}}^2 \right] \nonumber\\
%
& \leq & 2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*}))  \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&& + 2\mathbb{E}_{i_t^s}\left[ \| \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \mathcal{T}_{\eta_{t-1}^{*s}}\left(\gradf_{i_t^s}(w^{*})\right) - \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf(\tilde{w}^{s-1}) \right) \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
%
& = & 2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*}))  \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&&+ 2\mathbb{E}_{i_t^s}\left[ \| \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&& - 4 \left\langle \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf(\tilde{w}^{s-1}) \right), \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf(\tilde{w}^{s-1}) \right) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf(w^{*}))\right\rangle_{\changeHSSS{w_{t-1}^s}}  \nonumber\\
&&+ 2 \| \gradf(\tilde{w}^{s-1}) \|_{\changeHSSS{\tilde{w}^{s-1}}}^2 \nonumber\\
%
& = & 2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))\right.\nonumber\\
&& \left.\qquad \qquad \qquad \qquad \qquad + P^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*})) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&&+ 2\mathbb{E}_{i_t^s}\left[ \| \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \gradf_{i_t^s}(w_{t-1}^{s}) \right. \nonumber\\
&& \left. \qquad \qquad \qquad \qquad \qquad \qquad \qquad + \gradf_{i_t^s}(w_{t-1}^{s})- \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&& - 2 \| \gradf(\tilde{w}^{s-1}) \|_{\changeHSSS{w^{s-1}}}^2 \nonumber\\
%
& \leq & 4\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))  \|_{\changeHSSS{w_{t-1}^s}}^2 \right] \nonumber\\
&&+ 4\mathbb{E}_{i_t^s}\left[ \| P^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*})) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^{*})) \|_{\changeHSSS{w_{t-1}^s}}^2 \right]  \nonumber\\
&&+ 4\mathbb{E}_{i_t^s}\left[ \| \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \gradf_{i_t^s}(w_{t-1}^{s})\|_{\changeHSSS{w_{t-1}^s}}^2 \right]
\nonumber\\
&&
 + 4\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s} (w^{*})) \|_{\changeHSSS{w_{t-1}^s}}^2 \right] \nonumber\\
& \le & 
4\beta^2\changeHSSS{({\rm dist}(w^{s}_{t-1}, w^*))}^2 + 4\changeHSSS{\theta}^2\| \eta_{t-1}^{*s}\|_{\changeHSSS{w^*}}^2\|\gradf_{i_t^s}(w^{*})\|_{\changeHSSS{w^*}}^2\nonumber\\
&&
+ 4\mathbb{E}_{i_t^s}\left[ \| \mathcal{T}_{\tilde{\eta}_{t-1}^s} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - P^{\changeHSSS{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) \right.\nonumber\\
&& \qquad \qquad \qquad \qquad \qquad \qquad  \left.+ P^{\changeHSSS{w_{t-1}^s \leftarrow \tilde{w}^{s-1}}}(\gradf_{i_t^s}(\tilde{w}^{s-1})) -\gradf_{i_t^s}(w_{t-1}^{s})\|_{\changeHSSS{w_{t-1}^s}}^2 \right]
\nonumber\\
&&
 + 2(4\beta^2\changeHSSS{({\rm dist}(w^{s}_{t-1}, w^*))}^2 + 4\changeHSSS{\theta}^2\| \eta_{t-1}^{*s}\|_{\changeHSSS{w^*}}^2\|\gradf_{i_t^s}(w^{*})\|_{\changeHSSS{w^*}}^2) \nonumber
\\
& \leq & 12(\beta^2\changeHSSS{({\rm dist}(w^{s}_{t-1}, w^*))}^2+C^2\changeHSSS{\theta}^2\| \eta^{*s}_{t-1}\|_{\changeHSSS{w^*}}^2) \nonumber\\
&& + 8\changeHSSS{\theta}^2\changeHSSS{\mathbb{E}_{i_t^s}[\|\tilde{\eta}_{t-1}^s\|_{\changeHSSS{\tilde{w}^{s-1}}}^2]}\|\gradf_{i_t^s}(\tilde{w}^{s-1})\|_{\changeHSSS{\tilde{w}^{s-1}}}^2 + 8\beta^2\changeHSSS{({\rm dist}(\tilde{w}^{s-1}, w_{t-1}^s))}^2\nonumber\\
& \le &
4(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)(3({\rm dist}(w_{t-1}^s,w^*))^2 + 2({\rm dist}(\tilde{w}^{s-1},w_{t-1}^s))^2)
 \nonumber\\
%
& \le & 
4(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)(3({\rm dist}(w_{t-1}^s,w^*))^2 + 2({\rm dist}(\tilde{w}^{s-1},w^*)+{\rm dist}(w^*,w_{t-1}^s))^2)
 \nonumber\\
%
& \le & 
4(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)(7({\rm dist}(w_{t-1}^s,w^*))^2 + 4({\rm dist}(\tilde{w}^{s-1},w^*))^2),
\end{eqnarray*}
%
where the relation \changeHSSS{$\| 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}[\gradf_{i_t^s}(\tilde{w}^{s-1})]=\gradf(\tilde{w}^{s-1})$ and $\gradf(w^*)=0$ \changeHSSS{and that $\mathbb{E}_{i_t^s}$ is a linear operator}.
\changeHSSS{
Furthermore, we have evaluated the value $\mathbb{E}_{i_t^s}[\|\gradf_{i_t^s}(w_{t-1}^s) - \mathcal{T}_{\eta_{t-1}^{*s}}(\gradf_{i_t^s}(w^*))\|_{w_{t-1}^s}^2$ and again used the obtained relation in the third inequality.
}
\end{proof}



%
%


\begin{theorem}
\label{Thm:LocalConvergence}
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 \changeHS{the Hessian ${\rm Hess}f(w^*)$ of $f$ at \changeHSSS{$w^*$}} is positive definite), $D$ be the diameter of the compact set $\Omega$,
and $\zeta:=\changeHSSS{\sqrt{|\kappa|}D/\tanh(\sqrt{|\kappa|}D)}$.
Suppose Assumptions \ref{assump:T_isometry}--\ref{Appen_assump:C2} and \ref{assump:Omega}--\ref{Assump:ret_convex}.
\changeHSSS{
Let the constants $\beta$, $\mu$, $\theta$, and $C$ be the same as in Lemma \ref{AppenLem:UpperBoundVariance}.
Suppose that $\sigma>0$ is a constant in Lemma \ref{Lem:mean-value} satisfying \eqref{Eq:mean-value}.
Let $\alpha$} be a positive number satisfying \changeHSS{$0 < \alpha(\sigma-28\changehs{\zeta}\alpha(\beta^2+\changehs{\mu^2}C^2\changeHSSS{\theta}^2)) < 1$.}
It then follows that for any sequence $\{\tilde{w}^s\}$ generated by Algorithm \ref{Alg:R-SVRG} with a fixed step size $\alpha_t^s:=\alpha$ \changeHSSS{and $m_s:=m$} converging to $w^*$, there exists $K>0$ such that for all $s>K$,
\changeHSS{
\begin{eqnarray}
\label{Eq:LinearConvergence}
\mathbb{E}[({\rm dist}(\tilde{w}^s,w^*))^2]& \leq & 
\frac{4(1+ 16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2))}{\alpha m(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2))} 
\mathbb{E}[({\rm dist}(\tilde{w}^{s-1},w^*))^2].
\end{eqnarray}}
\end{theorem}

\begin{proof}
Since the function $x/\tanh(x)$ on $x$ monotonically increases in $[0, \infty)$, we have, \changeHSSS{from Lemma \ref{LemZhang}},
\begin{equation*}
a^2 \le \zeta b^2 + c^2 - 2bc \cos(A)
\end{equation*}
for any geodesic triangle in $\Omega$ with side lengths $a$, $b$, and $c$, \changeHSSS{since $\sqrt{|\kappa|}c \le \sqrt{|\kappa|}D$}.
Then, \changeHS{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 \changeHSSS{by considering the geodesic triangle with $w_{t-s}^s$, $w^*$, and $w_t^s$ in $\Omega$} as
%
\begin{eqnarray*}
&& \mathbb{E}_{i_t^s}\left[({\rm dist}(w_{t}^s, w^{*}))^2\right] \\
& \le & \mathbb{E}_{i_t^s}\left[\changeHSS{\zeta}({\rm dist}(w_{t-1}^s, w_t^s))^2+({\rm dist}(w_{t-1}^s, w^*))^2 - 2 \langle{\rm Exp}_{w_{t-1}^s}^{-1}(w_{t}^s), {\rm Exp}_{w_{t-1}^s}^{-1}(w^*)\rangle_{w_{t-1}^s}\right].
\end{eqnarray*}
It follows that
\begin{eqnarray*}
&& \mathbb{E}_{i_t^s}\left[({\rm dist}(w_{t}^s, w^{*}))^2 - ({\rm dist}(w_{t-1}^s, w^*))^2\right]\\
& \le & \mathbb{E}_{i_t^s}[\changeHSS{\zeta}({\rm dist}(w_{t-1}^s, w_{t}^s))^2 - 2 \langle -\alpha \xi_t^s, {\rm Exp}_{w_{t-1}^s}^{-1}(w^*)\rangle_{w_{t-1}^s}] \\
& = & \mathbb{E}_{i_t^s}[\changeHSSS{\zeta\| -\alpha \xi_t^s\|_{w_{t-1}^s}^2} + 2\alpha\langle\gradf(w_{t-1}^s), {\rm Exp}_{w_{t-1}^s}^{-1}(w^*)\rangle_{w_{t-1}^s}],
\end{eqnarray*}
where the last equality follows from $\mathbb{E}_{i_t^s}[\xi_t^s] = \gradf(w_{t-1}^{s})$.
Lemma \ref{Lem:mean-value}, together with the relation $f(w^*) \le f(w_{t-1}^s)$, yields that
\begin{eqnarray*}
\langle \gradf(w_{t-1}^s), {\rm Exp}_{w_{t-1}^s}^{-1}(w^*)\rangle_{w_{t-1}^s} \le -\frac{\sigma}{2} \|{\rm Exp}_{w_{t-1}^s}^{-1}(w^*)\|_{w_{t-1}^s}^2 \changeHSSS{=} -\frac{\sigma}{2} ({\rm dist}(w_{t-1}^s, w^*))^2.
\end{eqnarray*}
We thus obtain, by Lemma \ref{AppenLem:UpperBoundVariance},
\begin{align}
& \mathbb{E}_{\changeHSSS{i_t^s}}\left[({\rm dist}(w_{t}^s, w^{*}))^2 - ({\rm dist}(w_{t-1}^s, w^*))^2\right] \notag\\
 \le & \mathbb{E}_{\changeHSSS{i_t^s}}[\changeHSS{\zeta}\|\alpha\xi_t^s\|_{\changeHSSS{w_{t-1}^s}}^2  - \sigma\alpha({\rm dist}(w_{t-1}^s, w^*))^2] \notag\\
 \le & \mathbb{E}_{\changeHSSS{i_t^s}}[4\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)(7({\rm dist}(w_{t-1}^s, w^*))^2 \notag\\
& +  4({\rm dist}(\tilde{w}^{s-1}, w^*))^2) - \sigma\alpha({\rm dist}(w_{t-1}^s, w^*))^2] \notag\\
\label{Eq:forCorollary}
 \changeHK{=} & \alpha(28\changeHSS{\zeta}\alpha(\beta^2+\changeHSS{\changeHSSS{\mu}^2}C^2\changeHSSS{\theta}^2)-\sigma)({\rm dist}(w_{t-1}^s, w^*))^2 \\
 &+ 16\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)({\rm dist}(\tilde{w}^{s-1}, w^*))^2.\notag
\end{align}
It follows for the unconditional expectation operator $\mathbb{E}$ that
\begin{eqnarray*}
&& \mathbb{E}\left[({\rm dist}(w_{t}^s, w^{*}))^2\right] - \changeHSSS{\mathbb{E}}\left[({\rm dist}(w_{t-1}^s, w^*))^2\right]\\
& \le & \alpha(28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)-\sigma)\changeHSSS{\mathbb{E}}\left[({\rm dist}(w_{t-1}^s, w^*))^2\right] \\
&& + 16\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)\changeHSSS{\mathbb{E}}\left[({\rm dist}(\tilde{w}^{s-1}, w^*))^2\right].
\end{eqnarray*}
%
Summing over $t=1, \ldots, m$ of the inner loop on $s$-th epoch, we have 
\begin{align*}
&\mathbb{E}[({\rm dist}(w_{m}^s, w^{*}))^2] - \changeHSSS{\mathbb{E}}[({\rm dist}(w_{0}^s,w^{*}))^2] \\
\leq & \alpha(28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)-\sigma)
\sum_{t=1}^m \mathbb{E}[({\rm dist}(w_{t-1}^s,w^*))^2] \\
& + 16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2) \changeHSSS{\mathbb{E}}[({\rm dist}(\tilde{w}^{s-1},w^*))^2].
\end{align*}
Rearranging and using $w_0^s = \tilde{w}^{s-1}$, we obtain 
\begin{eqnarray*}
&& \alpha(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) \sum_{t=1}^{m} \mathbb{E}[({\rm dist}(w_{t}^s,w^*))^2] \\
& = & \changeHS{\alpha(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) \mathbb{E}\Biggl[\sum_{t=0}^{m-1} ({\rm dist}(w_{t}^s,w^*))^2 \\
&& \qquad \qquad \qquad \qquad \qquad \qquad \qquad + ({\rm dist}(w_{m}^s,w^*))^2 - ({\rm dist}(w_{0}^s,w^*))^2 \Biggr]}\\
& \leq & \mathbb{E}\left[\changeHS{({\rm dist}(w_{0}^s,w^*))^2 - ({\rm dist}(w_{m}^s,w^*))^2 +16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2) ({\rm dist}(w_{0}^s,w^*))^2  \right. \\
& & \left. -\alpha(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) ( ({\rm dist}(w_{0}^s,w^*))^2 - ({\rm dist}(w_{m}^s,w^*))^2)}\right]\\
& \leq & \changeHS{(1 - \alpha(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) + 16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2))\mathbb{E}[({\rm dist}(w_{\changeHSS{0}}^s,w^*))^2]}\\
& \leq & (1+ 16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) \mathbb{E}[({\rm dist}(\tilde{w}^{s-1},w^*))^2].
\end{eqnarray*}
%
Using $\tilde{w}^{\changeHK{s}} = g_m(\changeHK{w_1^s, \ldots,w_m^s})$ and Lemma \ref{AppenLem:KarcherMeanDistance}, we obtain
\begin{eqnarray*}
\mathbb{E}[({\rm dist}(\tilde{w}^s,w^*))^2]& \leq & 
\frac{4(1+ 16m\changeHSS{\zeta}\alpha^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2))}{\alpha m(\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2))} 
\mathbb{E}[({\rm dist}(\tilde{w}^{s-1},w^*))^2].
\end{eqnarray*}
This completes the proof.
\end{proof}
In the above theorem, we note that, from the definitions of $\beta$ and $\sigma$, $\beta$ can be chosen \changehs{to be} arbitrarily large and $\sigma$ to be arbitrarily small.
Therefore, \changeHSS{$\alpha = \sigma / 56\changeHSS{\zeta}(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)$, for example, satisfies $0 < \alpha (\sigma-28\changeHSS{\zeta}\alpha(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) < 1 $} for sufficiently large $\beta$ and small $\sigma$.

\changehs{
In fact, $\alpha$ satisfying the condition always exists for any values of $\beta$ and $\sigma$.
Let $\beta' := 28\zeta(\beta^2+\mu^2C^2\theta^2)$.
We can analyze the inequality $0<\alpha(\sigma-28\zeta\alpha(\beta^2+\mu^2C^2\theta^2))<1$, which is expressed as $0<\alpha(\sigma-\beta'\alpha)<1$, with the condition $\alpha>0$ more specifically as
\begin{equation}
\label{Eq:alpha}
\begin{cases}
\displaystyle
0 < \alpha < \frac{\sigma}{\beta'} \qquad & \text{if}\quad \sigma^2 - 4\beta' < 0,\\
\displaystyle
0 < \alpha < \frac{\sigma}{2\beta'},\ \frac{\sigma}{2\beta'} < \alpha < \frac{\sigma}{\beta'} \qquad & \text{if}\quad \sigma^2 - 4\beta' = 0,\\
\displaystyle
0 < \alpha < \frac{\sigma - \sqrt{\sigma^2-4\beta'}}{2\beta'},\ \frac{\sigma + \sqrt{\sigma^2-4\beta'}}{2\beta'} < \alpha < \frac{\sigma}{\beta'} \qquad & \text{if}\quad \sigma^2 - 4 \beta' >0.
\end{cases}
\end{equation}
}

\changehs{
Furthermore, we can show that the coefficient in the right-hand side of \eqref{Eq:LinearConvergence}, which can be written as $4(7+4m\beta'\alpha^2)/7m\alpha(\sigma-\beta'\alpha)$, is less than $1$ when $m$ is sufficiently large.
If $\alpha$ is fixed, $4(7+4m\beta'\alpha^2)/7m\alpha(\sigma-\beta'\alpha) \to 16\beta'/7(\sigma-\beta'\alpha)$ as $m \to \infty$, which is not necessarily less than $1$.
Thus, we again need to specifically analyze an appropriate value of $\alpha$, which should depend on $m$.
By calculating the derivative $r'(\alpha)$ of $r(\alpha):=4(7+4m\beta'\alpha^2)/7m\alpha(\sigma-\beta'\alpha)$ on $\alpha$,
we can show that $r(\alpha)$ takes the minimum value at
\begin{equation*}
\alpha = \frac{-7\beta' + \sqrt{49\beta'^2+28m\beta'\sigma^2}}{4m\beta'\sigma} =: \alpha_*,
\end{equation*}
which satisfies \eqref{Eq:alpha} when $m$ is sufficiently large, since $\lim_{m\to\infty}\alpha_* = 0$.
Note that we have $r'(\alpha_*) = 0$, which yields $4m\beta' \sigma\alpha_*^2 + 14\beta'\alpha_*-7\sigma = 0$.
This relation gives the minimum value $r(\alpha_*)$ as
\begin{equation*}
r(\alpha_*) = \frac{32(1-\sigma\beta'\alpha_*)}{2(7\beta'+2m\sigma^2)\alpha_*-7\sigma} \to 0 \quad (m \to \infty),
\end{equation*}
where we have used the facts that $\lim_{m\to \infty} \alpha_* = 0$ and $\lim_{m\to \infty} m\alpha_* =\infty$.
A more simple choice of $\alpha=1/\sqrt{m}$ also makes $r(\alpha)$ less than $1$ if $m$ is sufficiently large since
\begin{equation*}
r\left(\frac{1}{\sqrt{m}}\right) = \frac{4(7+4\beta')}{7(\sigma\sqrt{m}-\beta')}\to 0 \quad (m\to \infty).
\end{equation*}
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.
}

\changehs{
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 same as standard SVRG in the 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
$4(1+ 16m\changeHSS{\zeta}(\alpha_0^s)^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) / \alpha_m^s m(\sigma-28\changeHSS{\zeta}\alpha_0^s(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2))$.
This is larger than $4(1+ 16m\changeHSS{\zeta}(\alpha_0^s)^2(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2)) / \alpha_0^s m(\sigma-28\changeHSS{\zeta}\alpha_0^s(\beta^2+\changeHSSS{\mu}^2C^2\changeHSSS{\theta}^2))$, which is the coefficient in \eqref{Eq:LinearConvergence} with the fixed step size $\alpha=\alpha_0^s$.
Consequently, using decaying step sizes also yields a local convergence, but it gives a worse rate than with a fixed step size.
Both 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 firstly requires decaying step sizes to approach a neighborhood of a local minimum, and uses a fixed step size to achieve a faster linear convergence rate 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. However, analyzing the switching between decaying and fixed step sizes theoretically is left for future work.
}

\subsection{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 the exponential mapping and parallel translation are chosen as the 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.

\changeHSS{We obtain the following result as a corollary of the proof of Lemma \ref{AppenLem:UpperBoundVariance}, with $R={\rm Exp}$ and $\mathcal{T} = P$.
\begin{corollary} 
\label{Cor:UpperBoundVariance}
Consider Algorithm \ref{Alg:R-SVRG} with $R={\rm Exp}$ and $\mathcal{T} = P$, i.e., the exponential mapping and parallel translation case.
When each $\gradf_{n}$ is $\beta_0$-Lipschitz continuously differentiable, that is,
\begin{eqnarray*}
	\| P_{\changeHSSS{\gamma}}^{w \leftarrow z} ({\rm grad} f_n(z)) - {\rm grad} f_n(w) \|_{\changeHSSS{w}} & \leq & \beta_0 {\rm dist}(z,w),\qquad w, z \in \Omega,\ n = 1,2,\dots, N,
\end{eqnarray*}
where $\gamma$ is the geodesic connecting $z$ and $w$ and the upper bound of the variance of $\xi_t^s$ is given by
\begin{eqnarray}
\label{Eq:UpperBoundVariance}
	\mathbb{E}_{i_t^s}[\| \xi_t^s \|^2] &\leq &
\beta_0^2 (14({\rm dist}(w_{t-1}^s,w^*))^2 + 8{\rm dist}(\tilde{w}^{s-1},w^*))^2  ).
\end{eqnarray}
\end{corollary}
}
\changeHSS{
\begin{proof}
Putting $R={\rm Exp}$ and $\mathcal{T} = P$ in the middle of the proof of Lemma \ref{AppenLem:UpperBoundVariance}, 
we obtain
\begin{eqnarray*}
& & \hspace*{-1cm}\mathbb{E}_{i_t^s}[\| \xi_t^s \|^2]   \nonumber\\
%
& \le & 2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))  \|^2 \right]  \nonumber\\
&&+ 2\mathbb{E}_{i_t^s}\left[ \| P_{\gamma}^{w_{t-1}^s \leftarrow \tilde{w}_{s-1}} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*})) \|^2 \right]  \nonumber\\
&& - 2 \| P_{\gamma}^{w_{t-1}^s \leftarrow \tilde{w}_{s-1}} \left(\gradf(\tilde{w}^{s-1}) \right)\|^2.  \nonumber
\end{eqnarray*}
Similarly in Lemma \ref{AppenLem:UpperBoundVariance}, we have
\allowdisplaybreaks[1]
\begin{eqnarray*}
& & \mathbb{E}_{i_t^s}[\| \xi_t^s \|^2]   \nonumber\\
%& \leq & 2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))  \|^2 \right]  \nonumber\\
&&+ 2\mathbb{E}_{i_t^s}\left[ \| P_{\gamma}^{w_{t-1}^s \leftarrow \tilde{w}_{s-1}} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*})) \|^2 \right]  \nonumber\\
& \le & 
2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))  \|^2 \right]
 \nonumber\\
&&
+ 2\mathbb{E}_{i_t^s}\left[ \| P_{\gamma}^{w_{t-1}^s \leftarrow \tilde{w}_{s-1}} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \gradf_{i_t^s}(w_{t-1}^{s}) \right.\nonumber\\
&& \left. \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad + \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*})) \|^2 \right] \nonumber
\\
%
& \le & 
2\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s}(w^{*}))  \|^2 \right]
  \nonumber\\
&&
+ 4\mathbb{E}_{i_t^s}\left[ \| P_{\gamma}^{w_{t-1}^s \leftarrow \tilde{w}_{s-1}} \left(\gradf_{i_t^s}(\tilde{w}^{s-1}) \right) - \gradf_{i_t^s}(w_{t-1}^{s})\|^2 \right]
\nonumber\\
&&
 + 4\mathbb{E}_{i_t^s}\left[ \| \gradf_{i_t^s}(w_{t-1}^{s}) - P_{\gamma}^{w_{t-1}^{s} \leftarrow w^{*}}(\gradf_{i_t^s} (w^{*})) \|^2 \right] \nonumber
\\
%
& \leq & 
\beta_0^2 (6({\rm dist}(w_{t-1}^s,w^*))^2 + 4({\rm dist}(\tilde{w}^{s-1},w_{t-1}^s))^2  )
 \nonumber\\
%
& \le & 
\beta_0^2 (6({\rm dist}(w_{t-1}^s,w^*))^2 + 4({\rm dist}(\tilde{w}^{s-1},w^*)+{\rm dist}(w^*,w_{t-1}^s))^2  )
 \nonumber\\
%
& \le & 
\beta_0^2 (6({\rm dist}(w_{t-1}^s,w^*))^2 + 8({\rm dist}(\tilde{w}^{s-1},w^*))^2+8({\rm dist}(w^*,w_{t-1}^s))^2  )
 \nonumber\\
%
& = & 
\beta_0^2 (14({\rm dist}(w_{t-1}^s,w^*))^2 + 8({\rm dist}(\tilde{w}^{s-1},w^*))^2  ).
\end{eqnarray*}
%
This completes the proof.
\end{proof}
}

\begin{corollary}
Suppose the conditions in Theorem \ref{Thm:LocalConvergence} \changehs{with $\mu=\theta=0$, except that a positive number $\alpha$ satisfies $0<\alpha(\sigma-14\zeta\alpha\beta^2)<1$}.
Consider Algorithm \ref{Alg:R-SVRG} with a fixed step size $\alpha_t^s := \alpha$ and with the exponential mapping and parallel translation as retraction and vector transport, respectively.
For any sequence $\{\changeHK{\tilde{{w}}}^s\}$ generated by the algorithm, there exists $K>0$ such that for all $s>K$,
\begin{eqnarray*}
\mathbb{E}[({\rm dist}(\tilde{{w}}^s,{w}^*))^2]& \leq & 
\frac{4(1+ 8 m \zeta \alpha^2  \beta^2)}{ \alpha m (\sigma - 14 \zeta \alpha \beta^2)} 
\mathbb{E}[({\rm dist}(\tilde{{w}}^{s-1},{w}^*))^2].
\end{eqnarray*}
\end{corollary}

\begin{proof}
By using \eqref{Eq:UpperBoundVariance} instead of \eqref{Append_Eq:UpperBoundVariance}, we obtain
\begin{eqnarray*}
&& \mathbb{E}_{i_t^s}\left[({\rm dist}({w}_{t}^s, {w}^{*}))^2 - ({\rm dist}({w}_{t-1}^s, {w}^*))^2\right]\\
& \le & \alpha(14\zeta\alpha\beta^2-\sigma)\mathbb{E}_{i_{t}^s}\left[({\rm dist}({w}_{t-1}^s, {w}^*))^2 + 8\zeta\alpha^2\beta^2({\rm dist}(\tilde{{w}}^{s-1}, {w}^*))^2\right],
\end{eqnarray*}
instead of \eqref{Eq:forCorollary}.
%
Summing over $t=1, \ldots, m$ of the inner loop on the $s$-th epoch, we have 
\begin{equation*}
\begin{split}
&\mathbb{E}[({\rm dist}({w}_{m}^s, {w}^{*}))^2 - ({\rm dist}({w}_{0}^s,{w}^{*}))^2] \\
&\hspace*{1cm}\leq \alpha (14\alpha \beta^2 - \sigma) 
\sum_{t=1}^m \mathbb{E}[({\rm dist}({w}_{t-1}^s,{w}^*))^2]
 + 8 m \alpha^2  \beta^2 ({\rm dist}(\tilde{{w}}^{s-1},{w}^*))^2.
\end{split}
\end{equation*}
A similar discussion as in the proof of Theorem \ref{Thm:LocalConvergence} yields the claimed convergence rate.
\end{proof}




%
%
%
\section{Numerical comparisons}
\label{Sec:NumericalComparison}

This section compares the performance of R-SVRG(+) with the Riemannian extension of SGD, i.e., R-SGD, where the Riemannian stochastic gradient algorithm uses $\gradf_{i_t^s}(w_{t-1}^{s})$ instead of $\xi_t^s$ in (\ref{Eq:R-SVRG-Grad-paralleltrans}).
We also compare with R-SD, which is the Riemannian steepest descent algorithm with backtracking line search \cite[Chapters~4]{Absil_OptAlgMatManifold_2008}. \changeBM{We consider both \emph{fixed} step size and \changeHS{\emph{decaying}} step size sequences}. The \changehs{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 and $\lfloor \cdot \rfloor$ denotes the floor function. We select ten choices of $\alpha_0$ and consider three choices of $\lambda$ ($10^{-1}, 10^{-2}, and 10^{-3}$). \changeBM{In addition, since the global convergence analysis needs a \changehs{decaying} step size condition and the local convergence rate analysis holds for a fixed step size (\changeHSSS{Sections \ref{Sec:GlobalAnalysis} and \ref{Sec:LocalAnalysis}}),} we consider a {\it hybrid} step size sequence that follows the decay step size before the $s_{TH}$ epoch, and subsequently switches to a fixed step size. All experiments \changehs{herein} use $s_{TH}=5$. $m_s=5N$ is also fixed by following \cite{Johnson_NIPS_2013_s}\changeHK{, and  batch size is fixed to 10}. In all figures, the $x$-axis is the computational cost measured by the number of gradient computations divided by $N$. Algorithms are initialized randomly and are stopped when either the stochastic gradient norm is below $10^{-8}$ or the number of iterations exceeds $100$. It should be noted that all results except R-SD are the best-tuned results. All simulations are performed in Matlab on a 2.6 GHz Intel Core i7 machine with 16 GB RAM. 

%Additional numerical experiments are shown in Section \ref{AppenSec:AdditionalNumericalComparison}. 

\subsection{\changeHKK{Grassmann manifold and optimization tools}}
\changeHKK{This paper addresses problems on the Grassmann manifold, on which
an element is \changeHS{represented by a} $d \times r$ orthogonal matrix \mat{U} with \changeBM{orthonormal columns}, i.e., $\mat{U}^T\mat{U}=\mat{I}$. Two \changeBM{orthogonal} matrices \changeBM{represent the same element on} the Grassmann manifold if they are related by right multiplication of an $r\times r$ orthogonal matrix $\mat{O} \in \mathcal{O}(r)$. Equivalently, an element of the Grassmann manifold is identified with a set of $d \times r$ orthogonal matrices $[\mat{U}]: =\{\mat{U}\mat{O} :\mat{O} \in \mathcal{O}(r)\}$. That is, ${\rm Gr}(r,d) :={\rm St}(r,d)/ \mathcal{O}(r)$, where ${\rm St}(r,d)$ is the {\it Stiefel manifold} that is the set of matrices of size $d \times r$ with orthonormal columns. The Grassmann manifold has the structure of a Riemannian quotient manifold \cite[Section~3.4]{Absil_OptAlgMatManifold_2008}. 
%
The exponential mapping for the Grassmann manifold from $\mat{U}(0) := \mat{U} \in {\rm Gr}(r,d)$ in the direction of $\xi \in T_{\scriptsize \mat{U}(0)} {\rm Gr}(r,d)$ is given in a closed form as \cite[Section 5.4]{Absil_OptAlgMatManifold_2008}
\begin{eqnarray*}
\label{Eq:exponential_map}
\mat{U} (t) & = & [\mat{U}(0)  \mat{V}\ \  \mat{W}] 
	\left[
    		\begin{array}{c}
      		\cos t \Sigma  \\
      		\sin t \Sigma \\
    		\end{array}
	\right]
    \mat{V}^T,
\end{eqnarray*}
where $\xi=\mat{W} \Sigma \mat{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 \mat{U}(0)} {\rm Gr}(r,d)$ on the Grassmann manifold along $\gamma(t)$ with $\dot \gamma(0) = \mat{W} \Sigma \mat{V}^T$ is given in a closed form by
\begin{eqnarray*}
\label{Eq:parallel_translation}
\zeta(t) & = & \left( [\mat{U}(0) \mat{V}\ \  \mat{W}] 
	\left[
    		\begin{array}{c}
      		-\sin t \Sigma \\
      		\cos t \Sigma \\
    		\end{array}
	\right]
    \mat{W}^T + (\mat{I} - \mat{W}\mat{W}^T)
    \right) \zeta.
\end{eqnarray*}
%
The logarithm map of $\mat{U}(t)$ at $\mat{U}(0)$ on the Grassmann manifold is given by
\begin{eqnarray}
\label{Eq:logarithm_map}
\xi  &=& \log_{\scriptsize \mat{U}(0)}(\mat{U}(t)) \ = \ \mat{W} \arctan(\Sigma) \mat{V}^T,
\end{eqnarray}	
where the rank-$r$ singular value decomposition of $(\mat{U}(t) - \mat{U}(0) \mat{U}(0)^T  \mat{U}(t))\allowbreak(\mat{U}(0)^T \mat{U}(t))^{-1}$ is $\mat{W}\Sigma \mat{V}^T$.
}



\begin{figure}[t]
\begin{center}
	%\hspace*{-0.2cm}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{pca_train_error_N10000_d20_r5.eps}\\
		{\small (a) Train loss (enlarged).}
	\end{center} 
	\end{minipage}
	\hspace*{-0.1cm}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{pca_opt_gap_N10000_d20_r5.eps}\\
		{\small (b) Optimality gap. }
	\end{center} 
	\end{minipage}
	\hspace*{-0.1cm}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{pca_gnorm_N10000_d20_r5.eps}\\
		{\small (c) Norm of gradient.}
	\end{center} 
	\end{minipage}
\caption{Performance evaluations on PCA problem. }
\label{fig:PCA_results}
\end{center}		
\end{figure}



\subsection{\changeHKK{Problems on the Grassmann manifold and simulation results}}
\changeHKK{We} focus on three popular problems on the Grassmann manifold: PCA, low-rank matrix completion, and Karcher mean computation problems. In all these problems, full gradient methods, e.g., the steepest descent algorithm, become prohibitively computationally expensive when $N$ is very large;
the stochastic gradient approach is one promising way to achieve scalability.

{\bf The PCA problem.} Given \changeBM{an} orthonormal matrix projector $\mat{U} \in {\rm St}(r,d)$, the PCA problem is to minimize the sum of the squared residual errors between projected data points and the original data, as
\begin{equation}\label{Eq:PCA}
\begin{array}{ll}
{\displaystyle \min_{{\scriptsize \mat{U} \in {\rm St}(r,d)}}} & {\displaystyle \frac{1}{N}  \sum_{n=1}^N \| \vec{x}_n -  \mat{U}\mat{U}^T \vec{x}_n \|_2^2},
\end{array}
\end{equation}
where $\vec{x}_n$ is a data vector of size $d\times 1$. Problem \eqref{Eq:PCA} is equivalent to maximizing the function $\frac{1}{N} \sum_{n=1}^N \vec{x}_n^T\mat{U}\mat{U}^T\vec{x}_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 $\mat{U} \mapsto \mat{UO}$ for any orthogonal matrices $\mat{O}$ of size $r \times r$. Consequently, problem \eqref{Eq:PCA} is an optimization problem on the Grassmann manifold ${\rm Gr}(r,d)$.

\changeBM{Figures} \ref{fig:PCA_results}(a)--(c) show the train loss results, {\it optimality gap}, and gradient norm, respectively, where $N=10000$, $d=20$, and $r=5$. $\alpha_0$ is chosen as $\alpha_0 = 10^{-3},2\times10^{-3}, \ldots, 10^{-2}$. The optimality gap evaluates the performance against the minimum loss, which is obtained by the Matlab function {\tt pca}. Figure \ref{fig:PCA_results}(a) shows the enlarged results of the train loss, where all R-SVRG(+) algorithms yield better convergence properties. \changehs{Of all} the step size sequences of R-SVRG(+), the hybrid sequence performs best. Between R-SVRG and R-SVRG+, the latter shows superior performance for all step size \changeBM{sequences}. For the \changeHK{optimality} gap plots in Figure \ref{fig:PCA_results}(b), the results follow similar trends as those of \changehs{the} train loss plots. In Figure \ref{fig:PCA_results}(c), while the gradient norm of SGD \changehs{remains} at higher values, those of R-SVRG and R-SVRG+ converge to lower values in all cases.


\begin{figure}[t]
\begin{center}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{karchermean_train_error_N1000_d300_r5.eps}\\
		{\small (a) Train loss.}
	\end{center} 
	\end{minipage}
	\hspace*{-0.1cm}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{karchermean_train_error_enlarge_N1000_d300_r5.eps}\\
		{\small (b) Train loss (enlarged). }
	\end{center} 
	\end{minipage}
	\hspace*{-0.1cm}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{karchermean_gnorm_N1000_d300_r5.eps}\\
		{\small (c) Norm of gradient.}
	\end{center} 
	\end{minipage}
\caption{Performance evaluations on Karcher mean problem.}
\label{fig:KarcherMean_results}
\end{center}
\end{figure}




{\bf Computation of the Karcher mean of subspaces.} The Karcher mean was introduced as the notion of a \emph{mean} on Riemannian manifolds \cite{Karcher_1977_s}. It generalizes the notion of an ``average''  on the manifold. Given $N$ points on the Grassmann manifold with matrix representations $\mat{Q}_1,\ldots,\mat{Q}_N$, the Karcher mean is defined as the solution to the problem
\begin{equation} \label{Eq:KarcherMean}
\begin{array}{lll}
{\displaystyle \min_{{\scriptsize \mat{U} \in {\rm St}(r, d)}}} &{\displaystyle \frac{1}{2N} \sum_{n=1}^N ({\rm dist}(\mat{U}, \mat{Q}_n))^2},
\end{array}
\end{equation}
\changeBM{where ${\rm dist}(\cdot, \cdot)$ is the geodesic distance between the elements on the Grassmann manifold. The gradient of this loss function is $ \frac{1}{N} \sum_{n=1}^{N}  -{\rm Log}_{\scriptsize \mat{U}} (\mat{Q}_n)$, where ${\rm Log}$ is the log map defined in (\ref{Eq:logarithm_map}).} The Karcher mean on the Grassmann manifold ${\rm Gr}(r,d)$ is frequently used for computer vision problems, such as visual object and pose categorization \cite{Jayasumana_IEEETPAMI_2015_s}. Since recursive calculations of the Karcher mean are needed with each new visual image, the stochastic gradient algorithm becomes \changeBM{an appealing choice} for large datasets. 

Figures \ref{fig:KarcherMean_results}(a)--(c) show the results of the train loss, the enlarged train loss, and the norm of \changehs{the} gradient, respectively, where $N=1000$, $d=300$, and $r=5$. The ten choices of $\alpha_0$ are $\{0.1,0.2, \ldots, 1.0\}$. R-SVRG(+) outperforms R-SGD, and the final loss of R-SVRG(+) is less than that of R-SD. It should be noted that R-SVRG+, with fixed and \changehs{decaying} step sizes, decreases faster in the beginning, but eventually, R-SVRG converges to lower losses.




{\bf Low-rank matrix completion.} \changeBM{The matrix completion problem amounts to completing an incomplete matrix $\mat{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 $\mat{X}$, the rank-$r$ matrix completion problem amounts to solving the problem
\begin{equation}\label{Eq:MC_batch}
\begin{array}{ll}
{\displaystyle \min_{{\scriptsize \mat{U}} \in \mathbb{R}^{d \times r},\ {\scriptsize \mat{A}} \in \mathbb{R}^{r \times N}}} & \|\mathcal{P}_{\Omega}(\mat{UA}) - \mathcal{P}_{\Omega}(\mat X) \|_F^2,
\end{array}
\end{equation}
where the operator $\mathcal{P}_{\Omega}$ acts as $\mathcal{P}_{\Omega}(\mat{X}_{ij})=\mat{X}_{ij}$ if $(i,j) \in \Omega$ and $\mathcal{P}_{\Omega}(\mat{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 $\mat{X} = [\vec{x}_1, \vec{x}_2, \ldots, \vec{x}_n] $, problem (\ref{Eq:MC_batch}) is equivalent to the problem
\begin{equation}\label{Eq:MC}
\begin{array}{lll}
{\displaystyle \min_{{\scriptsize \mat{U}} \in \mathbb{R}^{d \times r},\ \vec{a}_n \in \mathbb{R}^{r}}} 
& 
{\displaystyle \frac{1}{N} \sum_{n=1}^N \|  \mathcal{P}_{\Omega_n}(\mat{U} \vec{a}_n) - \mathcal{P}_{\Omega_n}(\vec{x}_n) \|_2^2,
}
\end{array}
\end{equation}
where $\vec{x}_n \in \mathbb{R}^d$ and the operator $\mathcal{P}_{\Omega_n}$ \changehs{is} the sampling operator for the $n$-th column. Given \mat{U}, $\vec{a}_n$ in (\ref{Eq:MC}) admits a closed form solution. Consequently, problem (\ref{Eq:MC}) depends only on the column space of $\mat{U}$ and is on the Grassmann manifold \cite{Balzano_arXiv_2010_s}.}


\begin{figure}[t]
\begin{center}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{mc_test_MSE_N5000_d500_r5.eps}\\
		{\small (a) Test loss (synthetic).}
	\end{center} 
	\end{minipage}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]   
		{mc_jester_test_MSE_enlarge_N100_d24983_r5.eps}\\
		{\small (b) Test loss (Jester).}
	\end{center} 
	\end{minipage}
	\begin{minipage}[t]{0.32\textwidth}
	\begin{center}
		\includegraphics[width=\textwidth]{mc_movielens_test_MSE_enlarge_N3952_d6040_r5.eps}\\
		{\small (c) Test loss (MovieLens-1M).}
	\end{center} 
	\end{minipage}	
	% \hspace*{-0.1cm}
%	\begin{minipage}[t]{0.32\textwidth}
%	\begin{center}
%		\includegraphics[width=\textwidth]{results/mc/mc_gnorm_N5000_d500_r5.eps}\\
%		
%		{\small (b) Norm of gradient (synthetic). }
%		
%	\end{center} 
%	\end{minipage}
%	\hspace*{-0.1cm}
%	\begin{minipage}[t]{0.32\textwidth}
%	\begin{center}
%		\includegraphics[width=\textwidth]{results/mc_jester/mc_jester_train_MSE_enlarge_N100_d24983_r5.eps}\\
%		
%		{\small (c) Train loss (Jester).}
%		
%	\end{center} 
%	\end{minipage}
	% \vspace*{0.4cm}
	
	%\hspace*{-0.1cm}
%	\begin{minipage}[t]{0.32\textwidth}
%	\begin{center}
%		\includegraphics[width=\textwidth]   
%		{results/mc_jester/mc_jester_test_MSE_enlarge_N100_d24983_r5.eps}\\
%		{\small (d) Test loss (Jester).}
%	\end{center} 
%	\end{minipage}
%	\begin{minipage}[t]{0.32\textwidth}
%	\begin{center}
%		\includegraphics[width=\textwidth]{results/mc_movielens/mc_movielens_train_MSE_enlarge_2_N3952_d6040_r5.eps}\\
%		{\small (e) Train loss (MovieLens-1M).}
%	\end{center} 
%	\end{minipage}
	%\hspace*{-0.1cm}
%	\begin{minipage}[t]{0.32\textwidth}
%	\begin{center}
%		\includegraphics[width=\textwidth]{results/mc_movielens/mc_movielens_test_MSE_enlarge_N3952_d6040_r5.eps}\\
%		{\small (f) Test loss (MovieLens-1M).}
%	\end{center} 
%	\end{minipage}	
\caption{Performance evaluations on low-rank matrix completion problem.}
\label{fig:MC_results}
\end{center}	
\end{figure}


The proposed algorithms are also compared with Grouse \cite{Balzano_arXiv_2010_s}, 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$. Algorithms are initialized randomly as suggested in \cite{Kressner_BIT_2014_s}. The ten choices of $\alpha_0$ are $\{10^{-3},2\times10^{-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 of the total $Nd$ entries. \changeBMM{Figure \ref{fig:MC_results}(a) shows the results of loss on the test set $\Gamma$. \changehs{These} results show the superior performance of our proposed algorithms.}

%\changeBM{Figures} \ref{fig:MC_results}(a) and (b) show the results of loss on test set $\Gamma$ and the norm of \changehs{the} gradient, respectively. \changehs{These} results show the superior performance of our proposed algorithms.

Next, we consider Jester dataset 1 \cite{Goldberg_IR_2001_s}, 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\times10^{-6}, \ldots, 10^{-5}\}$ for SGD and SVRG(+) and $\{10^{-3},2\times10^{-3}, \ldots, 10^{-2}\}$ for Grouse. \changeBMM{Figure \ref{fig:MC_results}(b) shows the superior performance of R-SVRG(+) on the test set of the Jester dataset.}

%Figures \ref{fig:MC_results}(c) and (d) show the superior performance of R-SVRG(+) on both the training and test sets.

As a final test, we \changeBMM{also} compare the algorithms on the MovieLens-1M dataset, which is downloaded from \url{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\times10^{-5}, \ldots, 10^{-4}\}$. \changeBMM{Figure \ref{fig:MC_results}(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 others, and R-SVRG is better than R-SVRG+ in terms of the final test loss for all step size algorithms.}

%Figures \ref{fig:MC_results}(e) and (f) show the results on the training and test sets of all the algorithms except Grouse, which faces issues with convergence on this data set. R-SVRG(+) shows much faster convergence than others, and R-SVRG is better than R-SVRG+ in terms of the final test loss for all step size algorithms.







%
\section{Conclusion}
\changeBM{We have proposed a Riemannian stochastic variance reduced gradient algorithm (R-SVRG) 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 \changeBM{space, but is now extended to Riemannian manifolds}. The central difficulty of averaging, \changehs{adding}, and \changehs{subtracting} multiple gradients on a Riemannian manifold is handled with a vector transport. We proved that R-SVRG generates globally convergent sequences with a \changehs{decaying} step size condition and is locally linearly convergent with a fixed step size under some natural assumptions. Numerical comparisons on three popular problems on the Grassmann manifold suggested the superior performance of R-SVRG on various benchmarks.}


%  Supplementary material + references
%\clearpage
%\begin{comment}
\appendix

%\renewcommand\thefigure{\thesection.\arabic{figure}}  
%\renewcommand\thefigure{A.\arabic{figure}}  
%\setcounter{figure}{0} 

%\renewcommand\thetable{A.\arabic{table}}  
%\setcounter{table}{0} 

%\renewcommand{\theequation}{A.\arabic{equation}}
%\setcounter{equation}{0}


%\hrule height 1mm depth 0mm width 140mm
%\vspace{0.2cm}

%\begin{center}
%{\bf\LARGE Supplementary material}
%\end{center}
%\vspace{0.3cm}

%\hrule height 0.2mm depth 0mm width 140mm

%
%


%\clearpage
%\small
\bibliographystyle{unsrt}
\bibliography{stochastic_online_learning,matrix_tensor_completion}



\end{document}















R-SVRG_SIOPT_arXiv.bbl
\begin{thebibliography}{10}

\bibitem{Absil_OptAlgMatManifold_2008}
{\sc P.-A. Absil, R.~Mahony, and R.~Sepulchre}, {\em Optimization Algorithms on
  Matrix Manifolds}, Princeton University Press, 2008.

\bibitem{Allen-Zhu-arXiv_2016_s}
{\sc Z.~Allen-Zhu and E.~Hazan}, {\em Variance reduction for faster non-convex
  optimization}, tech. rep., arXiv preprint arXiv:1603.05643, 2016.

\bibitem{Allen-Zhu-arXiv_2015_s}
{\sc Z.~Allen-Zhu and Y.~Yan}, {\em Improved {SVRG} for non-strongly-convex or
  sum-of-non-convex objectives}, tech. rep., arXiv preprint arXiv:1506.01972,
  2015.

\bibitem{Balzano_arXiv_2010_s}
{\sc L.~Balzano, R.~Nowak, and B.~Recht}, {\em Online identification and
  tracking of subspaces from highly incomplete information}, in Allerton, 2010,
  pp.~704--711.

\bibitem{Bonnabel_IEEETAC_2013_s}
{\sc S.~Bonnabel}, {\em Stochastic gradient descent on {Riemannian} manifolds},
  IEEE Trans. on Automatic Control, 58 (2013), pp.~2217--2229.

\bibitem{Boumal_Manopt_2014_s}
{\sc N.~Boumal, B.~Mishra, P.-A. Absil, and R.~Sepulchre}, {\em Manopt: a
  {Matlab} toolbox for optimization on manifolds}, JMLR, 15 (2014),
  pp.~1455--1459.

\bibitem{Defazio_NIPS_2014_s}
{\sc A.~Defazio, F.~Bach, and S.~Lacoste-Julien}, {\em {SAGA}: A fast
  incremental gradient method with support for non-strongly convex composite
  objectives}, in NIPS, 2014, pp.~1646--1654.

\bibitem{Fisk_1965_s}
{\sc D.~L. Fisk}, {\em Quasi-martingales}, Trans. Amer. Math. Soc., 120 (1965).

\bibitem{Garber_arXiv_2015_s}
{\sc D.~Garber and E.~Hazan}, {\em Fast and simple {PCA} via convex
  optimization}, tech. rep., arXiv preprint arXiv:1509.05647, 2015.

\bibitem{Goldberg_IR_2001_s}
{\sc K.~Goldberg, T.~Roeder, D.~Gupta, and C.~Perkins}, {\em Eigentaste: a
  constant time collaborative filtering algorithm}, Inform. Retrieval, 4
  (2001), pp.~133--151.

\bibitem{huang2015riemannian}
{\sc W.~Huang, P.-A. Absil, and K.~A. Gallivan}, {\em A riemannian symmetric
  rank-one trust-region method}, Mathematical Programming, 150 (2015),
  pp.~179--216.

\bibitem{huang2015broyden}
{\sc W.~Huang, K.~A. Gallivan, and P.-A. Absil}, {\em A {B}royden class of
  quasi-{N}ewton methods for {R}iemannian optimization}, SIAM Journal on
  Optimization, 25 (2015), pp.~1660--1685.

\bibitem{Jayasumana_IEEETPAMI_2015_s}
{\sc S.~Jayasumana, R.~Hartley, M.~Salzmann, H.~Li, and M.~Harandi}, {\em
  Kernel methods on riemannian manifolds with gaussian rbf kernels}, IEEE
  Trans. Pattern Anal. Mach. Intell., 37 (2015), pp.~2464 -- 2477.

\bibitem{Johnson_NIPS_2013_s}
{\sc R.~Johnson and T.~Zhang}, {\em Accelerating stochastic gradient descent
  using predictive variance reduction}, in NIPS, 2013, pp.~315--323.

\bibitem{Karcher_1977_s}
{\sc H.~Karcher}, {\em Riemannian center of mass and mollifier smoothing},
  Comm. Pure Appl. Math., 30 (1977), pp.~509--541.

\bibitem{Kasai_arXiv_2016}
{\sc H.~Kasai, H.~Sato, and B.~Mishra}, {\em Riemannian stochastic variance
  reduced gradient on grassmann manifold}, arXiv preprint: arXiv:1605.07367,
  (2016).

\bibitem{Konecny_arXiv_2013}
{\sc J.~Kone\v{c}n\'y and P.~Richt\'arik}, {\em Semi-stochastic gradient
  descent methods}, tech. rep., arXiv preprint arXiv:1312.1666, 2013.

\bibitem{Kressner_BIT_2014_s}
{\sc D.~Kressner, M.~Steinlechner, and B.~Vandereycken}, {\em Low-rank tensor
  completion by {Riemannian} optimization}, BIT Numer. Math., 54 (2014),
  pp.~447--468.

\bibitem{Mairal_SIAMJOPT_2015}
{\sc J.~Mairal}, {\em Incremental majorization-minimization optimization with
  application to largescale machine learning}, SIAM J. Optim., 25 (2015),
  pp.~829--855.

\bibitem{Meyer_ICML_2011}
{\sc G.~Meyer, S.~Bonnabel, and R.~Sepulchre}, {\em Linear regression under
  fixed-rank constraints: A {R}iemannian approach}, in ICML, 2011.

\bibitem{Mishra_ICDC_2014_s}
{\sc B.~Mishra and R.~Sepulchre}, {\em {R3MC}: A {Riemannian} three-factor
  algorithm for low-rank matrix completion}, in IEEE CDC, 2014, pp.~1137--1142.

\bibitem{Mishra_SIAMJOPT_2016}
\leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em {Riemannian}
  preconditioning}, SIAM J. Optim., 635-660 (2016).

\bibitem{Roux_NIPS_2012_s}
{\sc N.~L. Roux, M.~Schmidt, and F.~R. Bach}, {\em A stochastic gradient method
  with an exponential convergence rate for finite training sets}, in NIPS,
  2012, pp.~2663--2671.

\bibitem{Shalev_arXiv_2015_s}
{\sc S.~Shalev-Shwartz}, {\em {SDCA} without duality}, tech. rep., arXiv
  preprint arXiv:1502.06177, 2015.

\bibitem{Shalev_arXiv_2012_s}
{\sc S.~Shalev-Shwartz and T.~Zhang}, {\em Proximal stochastic dual coordinate
  ascent}, tech. rep., arXiv preprint arXiv:1211.2717, 2012.

\bibitem{Shalev_JMLR_2013_s}
\leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em Stochastic dual
  coordinate ascent methods for regularized loss minimization}, JMRL, 14
  (2013), pp.~567--599.

\bibitem{Shamir_arXiv_2015_s}
{\sc O.~Shamir}, {\em Fast stochastic algorithms for {SVD} and {PCA}:
  Convergence properties and convexity}, tech. rep., arXiv preprint
  arXiv:1507.08788, 2015.

\bibitem{Zhang_NIPS_2016}
{\sc H.~Zhang, S.~J. Reddi, and S.~Sra}, {\em Fast stochastic optimization on
  {Riemannian} manifolds}, in Accepted for publication in NIPS, 2016.

\bibitem{Zhang_JMLR_2016_s}
{\sc H.~Zhang and S.~Sra}, {\em First-order methods for geodesically convex
  optimization}, in COLT, 2016.

\bibitem{Zhang_SIAMJO_2014_s}
{\sc Y.~Zhang and L.~Xiao}, {\em Stochastic primal-dual coordinate method for
  regularized empirical risk minimization}, SIAM J. Optim., 24 (2014),
  pp.~2057--2075.

\end{thebibliography}
