Keywords: Mirror descent, Riemannian optimization, Reparameterization, Non-asymptotic convergence analysis, Stiefel manifolds

MSCcodes: 65K05, 90C06, 90C30

1 Introduction↩︎

In this paper, we study the Riemannian optimization problem \[\min_{x\in\mathcal{M}} f(x),\] where \(\mathcal{M}\) is a Riemannian manifold, and \(f:\mathcal{M}\rightarrow\mathbb{R}\) is a smooth function on \(\mathcal{M}\). This problem can be viewed as constrained nonlinear optimization with a manifold constraint. The manifold structure enables us to design efficient algorithms and perform detailed analysis.

Riemannian optimization finds wide applications in practice. For instance, it can be applied to dictionary learning [2] and matrix decomposition [3], [4]. In deep learning, imposing orthogonal constraints on weight matrices can enhance training stability and improve generalization performance [5], [6]. Such constraints can be modeled by optimizing objective functions over the Stiefel manifold.

One line of research in this field focuses on generalizing unconstrained optimization algorithms to Riemannian manifolds. For example, Riemannian Gradient Descent [7], Riemannian Proximal Point [8], Riemannian Momentum [9], and Riemannian Stochastic Variance Reduction Gradient [10], [11] have been explored. However, Mirror Descent (MD), an important class of algorithms in Euclidean space, has been surprisingly unexplored in the existing literature. This paper takes the first step to fill this gap.

The MD algorithm was proposed by Nemirovski and Yudin [12] as a generalization of gradient descent and was later popularized by the widely used formulation of Beck and Teboulle [13]. Consider the Euclidean optimization problem \(\min_{x\in\Omega} f(x)\), where \(\Omega\subseteq\mathbb{R}^d\) is a closed convex feasible set. Let \(\psi:\Omega\to\mathbb{R}\) be a differentiable and strongly convex function, often called a potential function. The associated Bregman divergence [14] is defined as \[D_{\psi}(x,y)=\psi(x)-\psi(y)-\langle\nabla\psi(y),x-y\rangle.\] Given a step size \(\eta_t>0\), the MD update is \[x_{t+1}=\arg\min_{x\in\Omega}\left\{\left\langle \eta_t\nabla f(x_t),\,x-x_t\right\rangle+D_\psi(x,x_t)\right\},\] where the Bregman term can be interpreted as a geometry-adapted regularizer. When \(\Omega=\mathbb{R}^d\), the first-order optimality condition yields the implicit dual update \[\nabla\psi(x_{t+1})=\nabla\psi(x_t)-\eta_t\nabla f(x_t).\] Different choices of \(\psi\) recover algorithms suited to different geometries. For instance, taking \(\psi(x)=\tfrac12\|x\|_2^2\) reduces MD to the standard gradient descent update. When \(\Omega\) is the probability simplex, the entropy potential \(\psi(x)=\sum_{i=1}^d x_i\log x_i\) (with the convention \(0\log 0=0\)) leads to Entropy Mirror Descent, and \(D_\psi\) becomes the Kullback–Leibler divergence [15]. By adapting to the problem geometry with mild dependence on dimension, MD has found widespread applications, including image processing [13], policy optimization [16], [17], and neural network training [18], [19].

Subsequent work further broadened the study of MD, with particular emphasis on its extensions and connections to other methods. Nemirovski and Juditsky [15] considered using MD to solve stochastic optimization and stochastic saddle-point problems. Duchi et al. [20] generalized these results to composite functions. Recent research has focused on providing alternative explanations for MD. Gunasekar et al. [21] demonstrated that MD can be interpreted as the gradient flow on the manifold as the step sizes tend to zero. Amid and Warmuth [22] observed that MD becomes equivalent to gradient descent under different parameterizations. Li et al. [23] utilized this observation to explain the implicit regularization in over-parameterized models. Lei and Zhou [24] studied online mirror descent in a setting where data arrive sequentially.

The literature on Riemannian optimization is extensive. We briefly review several related results. Zhang and Sra [7] establish the first non-asymptotic convergence rates of geodesic gradient descent in both deterministic and stochastic settings. Later, Bento et al. [8] analyzed Riemannian subgradient and proximal point methods; however, their results rely on the assumption that the underlying manifold is Hadamard. More recently, Srinivasan and Wirth [25] extended this non-asymptotic convergence analysis to general Riemannian manifolds by constructing a suitable potential function. Alimisis et al. [9] studied momentum-based acceleration of first-order methods on manifolds. In the stochastic setting, Zhang and Sra [10] and Sato et al. [11] applied variance-reduction techniques and derived improved complexity bounds for finite-sum problems. Tripuraneni et al. [26] generalized Polyak–Ruppert averaging to the manifold setting and showed that the resulting averaged iterates enjoy improved convergence guarantees for stochastic gradient descent on manifolds. Recently, Shu et al. [27] introduced a quasilinearization framework to analyze proximal-based methods on Hadamard manifolds.

Our contributions can be summarized as follows:

• We generalize Euclidean Mirror Descent to optimization on Riemannian manifolds by developing a Riemannian Mirror Descent (RMD) framework via reparameterization. Building on the proposed framework, we further propose a stochastic variant of RMD (stochastic RMD) suitable for large-scale settings in which only stochastic gradient estimates are available.

• We establish non-asymptotic convergence guarantees for both RMD and stochastic RMD. Let \(T\) denote the total number of iterations. For geodesically convex objectives, we prove that the objective error decays at a sublinear rate of \(\mathcal{O}(1/T)\) in the deterministic setting and \(\mathcal{O}(1/T^{1/3})\) in the stochastic setting. For smooth nonconvex objectives, we show that the average squared gradient norm decays at a rate of \(\mathcal{O}(1/T)\) in the deterministic case and \(\mathcal{O}(1/T^{1/2})\) in the stochastic case. Our theoretical results fill the gap of [21], [22].

• As an application to optimization on the Stiefel manifold, we show that our RMD framework recovers the Curvilinear Gradient Descent (CGD) method of [1]. Moreover, specializing stochastic RMD to the Stiefel setting yields a stochastic extension of CGD, which we term Stochastic Curvilinear Gradient Descent (SCGD). The proposed SCGD is scalable and computationally efficient for large-scale manifold optimization problems and may be of independent interest. We further validate its effectiveness through numerical experiments.

Outline of this paper. In Section 2, we introduce the notation and assumptions used throughout the paper. In Section 3, we present the framework of RMD and Stochastic RMD, and provide several concrete instantiations. In Section 4, we analyze the convergence rates of our algorithms. In Section 5, we demonstrate the effectiveness of the proposed methods through numerical experiments. We conclude with a discussion and future research directions in Section 6.

2 Notation and Theoretical Background↩︎

We first recall some basic notions in Riemannian geometry.

Differentials and gradients are indispensable for optimization on manifolds.

For a differentiable map \(F:\mathcal{M}\to\mathcal{N}\) between manifolds, we denote by \(DF(x): \mathcal{T}_x\mathcal{M} \to \mathcal{T}_{F(x)}\mathcal{N}\) for its differential at \(x\in \mathcal{M}\), and write \(DF(x)[v]\in T_{F(x)}\mathcal{N}\) for its action on \(v\in \mathcal{T}_x\mathcal{M}\). We also denote by \(\nabla\) the Levi–Civita connection on \((\mathcal{M},g)\). In what follows, we will mainly use it in the form of the covariant derivative along a smooth curve.

In this paper, we assume that \((\mathcal{M},g)\) is connected and complete (in the Riemannian sense). Equivalently, with the distance function \(d_{\mathcal{M}}\) induced by Riemannian metric, \((\mathcal{M},d_{\mathcal{M}})\) is a complete metric space, or every geodesic can be extended for all time, i.e., the exponential map is defined on \(\mathcal{T}\mathcal{M}\). By the Hopf–Rinow theorem, for any \(x,y\in\mathcal{M}\), there exists a minimizing geodesic segment \(\gamma:[0,1]\to\mathcal{M}\) joining them, satisfying \(\gamma(0)=x, \gamma(1)=y\), and the length of \(\gamma\) is equal to \(d_{\mathcal{M}}( x,y)\). But this minimizing geodesic \(\gamma\) is not necessarily unique.

In general, since the geodesic equation \(\frac{d^{2} \gamma ( t)}{dt^{2}} =0\) is of second order, the exponential map does not have a closed form expression and is extremely expensive to compute. In this case, a retraction map is used to approximate it.

A typical example of retraction is the unit sphere \(\mathbb{S}^{d-1}\). For any point \(x\in \mathbb{S}^{d-1}\) and \(v\in \mathcal{T}_{x}\mathbb{S}^{d-1}\), a retraction is \(R_{x}(v) =\frac{x+v}{\|x+v\|}\), that is, the projection of \(x+v\). The following definition measures how close a retraction is to the exponential map, in view of Lemma 2.

With the notion of geodesics, we can define convexity on Riemannian manifolds.

Note that we do not require the minimizing geodesic segment connecting \(x\) and \(y\) to be unique. Since \((\mathcal{M},g)\) is assumed to be complete, by the Hopf–Rinow theorem any two points in \(\mathcal{M}\) can be joined by a minimizing geodesic segment; hence the whole manifold \(\mathcal{M}\) is geodesically convex. For an arbitrary subset \(B\subseteq \mathcal{M}\), the convex hull of \(B\) is defined as the intersection of all geodesically convex subsets \(A\subseteq \mathcal{M}\) that contain \(B\). By construction, this convex hull is geodesically convex.

Next, we introduce the convexity of a differentiable function. Recall that in Euclidean space, a function \(h\) is convex if for any \(x,y\), \[h( tx+( 1-t) y) \leqslant th( x) +( 1-t) h( y) ,\quad \forall t\in [ 0,1] .\] For function on the Riemannian manifold, we adopt the definition of Zhang and Sra [7].

If we replace the exponential map \(\mathrm{Exp}_{x}\) with retraction \(R_{x}\) in Definition 6, we get the definition of retraction convex function. If a function is both geodesically convex and differentiable, we have \[f\left(\mathrm{Exp}_{x}( v)\right) -f( x) \geqslant \langle \nabla f( x) ,v\rangle _{x} .\]

In optimization, a common assumption is the Lipschitz property of gradient.

We conclude by introducing some notation used throughout the paper. Given a matrix \(A=[a_{ij}]\) and index sets \(\mathcal{I},\mathcal{J}\), we denote by \(A(\mathcal{I},\mathcal{J})=[a_{ij}]_{i\in\mathcal{I},\,j\in\mathcal{J}}\) the submatrix of \(A\) with row indices \(\mathcal{I}\) and column indices \(\mathcal{J}\). We use the shorthand \([n]=\{1,2,\dots,n\}\). Finally, \(\nabla_\mathcal{M}\) means taking gradient on manifold \(\mathcal{M}\).

3 Riemannian Mirror Descent: Algorithms and Convergence Theorems↩︎

In this section, we formulate Riemannian Mirror Descent (RMD) and stochastic RMD on manifolds. We establish non-asymptotic convergence guarantees for both algorithms. We then present several concrete examples arising from our RMD framework. In particular, we obtain Stochastic Curvilinear Gradient Descent (SCGD) as a special case of particular interest.

3.1 Motivation and Algorithm Framework↩︎

In this subsection, the ambient space is \(\mathbb{R}^d\) equipped with its standard coordinates, and \(\nabla\) and \(\nabla^2\) denote the usual Euclidean gradient and Hessian, respectively. Our main motivation comes from the observation of Raskutti and Mukherjee [28]. They consider a different Riemannian structure on \(\mathbb{R}^d\). Given a strongly convex smooth function \(\psi\), for each point \(x \in \mathbb{R}^{d}\), we define \[g_{x}^{\psi}(v,u) = v^{\mathsf{T}} \nabla^{2} \psi(x)u, \quad \forall u,v \in \mathbb{R}^{d}.\] Notice that \(g_{x}^{\psi}(\cdot,\cdot)\) is an inner product on \(\mathbb{R}^{d}\) since \(\nabla^{2} \psi(x)\) is positive definite. Under the standard identification \(\mathcal{T}_x\mathbb{R}^d\simeq\mathbb{R}^d\), the matrix representation of \(g_{x}^{\psi}(\cdot,\cdot)\) is \(g^{\psi}(x)=\nabla^2\psi(x)\). We denote by \(\mathcal{M}_{\psi} = (\mathbb{R}^{d}, g^{\psi})\) the Riemannian manifold obtained by equipping \(\mathbb{R}^{d}\) with the Riemannian metric \(g^{\psi}\). Let \(\psi^*\) be the conjugate of \(\psi\), defined by \(\psi^{*}(y) = \sup_{x} \left\{\langle y, x\rangle - \psi(x)\right\}\). Then \(\psi^{*}\) is also strongly convex and smooth. Let \(\mathcal{M}_{\psi^{*}} = (\mathbb{R}^{d}, g^{\psi^{*}})\) be the manifold in the dual variable \(y\). Then the map \(x\mapsto \nabla\psi(x)\) is from \(\mathcal{M}_{\psi}\) to \(\mathcal{M}_{\psi^{*}}\), and \(\nabla \psi^{*}(\nabla \psi(x))=x\), \(\nabla^{2} \psi(x) \nabla^{2} \psi^*(\nabla\psi(x))=\mathrm{Id}_{\mathbb{R}^d}\).

The update in mirror descent (MD) can be rewritten as \[\begin{align} y_{t+1} & = \nabla \psi(x_{t}) - \eta_{t} \nabla f(x_{t}), \label{eq:MD1} \\ x_{t+1} & = \nabla \psi^{*}(y_{t+1}). \end{align}\tag{1}\] In the first step, \(\nabla\psi\) maps \(x_t\in\mathcal{M}_\psi\) to \(\nabla\psi(x_t)\in\mathcal{M}_{\psi^*}\), and the update subtracts \(\eta_t\nabla f(x_t)\) in the dual variable. Under our setting, \(\nabla\psi\) is invertible with inverse \((\nabla\psi)^{-1}=\nabla\psi^*\). The invertible map \(\nabla \psi\) pulls back \(f\) on \(\mathcal{M}_{\psi}\) to \(\tilde{f}=f\circ(\nabla \psi)^{-1}=f\circ\nabla\psi^{*}\) on \(\mathcal{M}_{\psi^{*}}\). We can calculate the Riemannian gradient of \(\tilde{f}\) at \(\nabla \psi(x_{t})\) on \(\mathcal{M}_{\psi^{*}}\). Note that \(y_{t}=\nabla \psi(x_{t})\), we have \[\begin{align} \nabla_{\mathcal{M}_{\psi^{*}}} \tilde{f}(y_{t}) &= (g^{\psi^{*}}(y_{t}))^{-1} \nabla \tilde{f} (y_{t})\\ &= (\nabla^{2} \psi^{*}(y_{t}))^{-1} \cdot \nabla\nabla\psi^{*}(y_{t}) \cdot \nabla f(\nabla\psi^{*}(y_{t})) \\ &= (\nabla^{2} \psi^{*}(y_{t}))^{-1} \cdot \nabla^{2} \psi^{*}(y_{t}) \cdot \nabla f(x_{t}) = \nabla f(x_{t}), \end{align}\] where the first line is the coordinate expression for Riemannian gradient and in the second line we use the chain rule. This computation means that \(\nabla f(x_{t})\) is equivalent to the gradient of \(\tilde{f}\) in \(\mathcal{M}_{\psi^*}\) at \(y_t\). As a result, the first step is a gradient step \(y_{t+1}=y_{t}-\eta_{t}\nabla_{\mathcal{M}_{\psi^{*}}} \tilde{f}(y_{t})\) in the space \(\mathcal{M}_{\psi^*}\). In the second step, the inverse of \(\nabla \psi\), i.e., \(\nabla \psi^{*}\) maps \(y_{t+1}\) back to the primal manifold \(\mathcal{M}_{\psi}\).

The MD update process can be viewed as gradient descent via reparameterization. More precisely, by changing of variables \(y=\nabla\psi(x)\), the algorithm maps the primal iterate \(x_t\in\mathcal{M}_\psi\) to the dual coordinate \(y_t\in\mathcal{M}_{\psi^*}\), performs gradient descent in the dual space \(\mathcal{M}_{\psi^{*}}\), and then maps back to \(\mathcal{M}_\psi\) via \(\nabla\psi^*\). Motivated by this idea, we propose the Riemannian Mirror Descent (RMD) algorithm (Algorithm 1) as an extension of Euclidean MD to general Riemannian manifolds.

In Algorithm 1, we replace the global map \(x\mapsto \nabla\psi(x)\) by a sequence of local diffeomorphisms \(\varphi_t\). At each iteration, we choose a local diffeomorphism \(\varphi_t:\hat{\mathcal{N}}_{r_t}(x_t)\to \mathcal{M}_t\) onto a reparameterized manifold \(\mathcal{M}_t\). Here \(\hat{\mathcal{N}}_{r}(x)=\{y\in\mathcal{M}:d_{\mathcal{M}}(x,y)< r\}\) is the open geodesic ball. For simplicity, when \(\mathrm{inj}(\mathcal{M})<+\infty\), we assume \(r_t< \mathrm{inj}(\mathcal{M})\) so that \(\hat{\mathcal{N}}_{r_t}(x)\) is diffeomorphic to open ball \(\{v\in \mathcal{T}_x \mathcal{M}:\|v\|<r\}\), and so is \(\mathcal{M}_t\). When \(\mathrm{inj}(\mathcal{M})=+\infty\), \(r_t\) is free to be finite or infinite. We then perform a gradient step on \(\mathcal{M}_t\) via retraction \(R^t\) and map the iterate back via \(\varphi_t^{-1}\). Updating \(y_{t+1}\) in Algorithm 1 corresponds to 1 . In essence, RMD amounts to carrying out gradient descent after a local reparameterization. One benefit is that the retraction map can be easier to compute in the reparameterized space.

Next, we present the convergence theorem for RMD in Algorithm 1. To establish these results, we introduce the following assumptions.

We are now in a position to state the convergence theorem.

A detailed proof is given in Section 4.

3.2 Examples of Riemannian Mirror Descent↩︎

Algorithm 1 establishes an abstract framework for RMD. By choosing different mirror maps \(\varphi_t\) and manifolds \(\mathcal{M}_t\), one can recover a variety of algorithms as special cases of this framework. In what follows, we present several concrete examples to illustrate our framework.

In this example, the mirror map is in fact an isometry. Fix \(x\in \mathcal{M}_{\psi }\) and \(u,v\in \mathcal{T}_{x}\mathcal{M}_{\psi }\). Let \(y=\nabla\psi(x)\) and denote \(\hat{u} = D(\nabla\psi)(x)[u]=\nabla^2\psi(x)\,u\) and \(\hat{v}=\nabla^2\psi(x)\,v\). Then \[g^{\psi^*}_{y}(\hat{u},\hat{v}) = \hat{u}^{\mathsf T}\nabla^2\psi^*(y)\hat{v} = u^{\mathsf T}\nabla^2\psi(x)\,\nabla^2\psi^*(\nabla\psi(x))\,\nabla^2\psi(x)\,v = u^{\mathsf T}\nabla^2\psi(x)\,v = g^{\psi}_{x}(u,v),\] where we used the identity \(\nabla^2\psi^*(\nabla\psi(x))=(\nabla^2\psi(x))^{-1}\). Therefore, using the conjugacy relation between \(\psi\) and \(\psi^*\), the mirror map \(\nabla\psi\) in classical MD is a Riemannian isometry between \(\mathcal{M}_{\psi}\) and \(\mathcal{M}_{\psi^*}\).

On a general Riemannian manifold, one cannot in general expect to construct a globally defined Riemannian isometry to serve as a mirror map. However, the tangent space provides a natural local parameterization of the manifold. By using the exponential map as a mirror map, Algorithm 1 recovers the Geodesic Gradient Descent algorithm.

Indeed, \(\mathrm{Exp}_{x_t}(0)=x_t\) and \(D\mathrm{Exp}_{x_t}(0)=\mathrm{Id}_{\mathcal{T}_{x_t}\mathcal{M}}\). Since \(\varphi_t(x_t)=\mathrm{Exp}_{x_t}^{-1}(x_t)=0\), we obtain \[y_{t+1} =R_{\varphi_t(x_t)} \left(-\eta_t\,D\varphi_t(x_t)[\nabla f(x_t)]\right) =-\eta_t \nabla f(x_t),\] and then \[x_{t+1}=\varphi_t^{-1}(y_{t+1}) =\mathrm{Exp}_{x_t} \left(-\eta_t \nabla f(x_t)\right),\] which is exactly Geodesic Gradient Descent.

Next, we discuss how to apply our framework to optimization over the Stiefel manifold.

From this definition, the Stiefel manifold with \(p=n\) is \(\mathrm{St}(n,n)=\{X\in\mathbb{R}^{n\times n}:X^{\mathsf T}X=I_n\}.\) Then \(I_n\in \mathrm{St}(n,n)\) and \[\mathcal{T}_{I_n}\mathrm{St}(n,n)=\mathrm{skew}(n),\] where \(\mathrm{skew}(n)=\{W\in\mathbb{R}^{n\times n}: W+W^{\mathsf T}=0\}\) denotes the set of \(n\times n\) skew-symmetric matrices. Therefore, the Cayley transform provides a natural local reparameterization of \(\mathrm{St}(n,n)\) near the identity. For any \(W\in \mathrm{skew}(n)\), the Cayley transform is defined by \[C(W)=(I_n-W)^{-1}(I_n+W)\in \mathrm{St}(n,n).\]Then \(C\) is a diffeomorphism between \(\mathrm{skew}(n)\) and the open subset \(\{X\in \mathrm{St}(n,n): X+I_n \text{ is invertible}\}\). On this subset, its inverse is \[\varphi(X)=C^{-1}(X)=(X-I_n)(X+I_n)^{-1}.\] Viewing \(\varphi\) as a smooth map into \(\mathrm{skew}(n)\), we have \[D\varphi(I_n)[W]=\tfrac12 W,\qquad \forall\, W\in \mathrm{skew}(n).\]Thus, when \(x_t=I_n\) and using the Euclidean retraction on \(\mathrm{skew}(n)\), the update in dual variable \(y\) is \[y_{t+1}=-\eta_t\,D\varphi(I_n)[\nabla f(I_n)] = -\tfrac{\eta_t}{2}\nabla f(I_n)\in \mathrm{skew}(n),\] and mapping back yields \[x_{t+1}= \varphi^{-1}(y_{t+1}) = C(y_{t+1}) = (I_n-y_{t+1})^{-1}(I_n+y_{t+1})\in \mathrm{St}(n,n).\] For a general base point \(x_t=X_0\in \mathrm{St}(n,n)\), since \(\mathrm{St}(n,n)\) is a homogeneous manifold, we use the translation \[O_{X_0}:\mathrm{St}(n,n)\to\mathrm{St}(n,n),\qquad O_{X_0}(X)=X X_0^{\mathsf T},\] which satisfies \(O_{X_0}(X_0)=I_n\). On a neighborhood of \(X_0\) where \(O_{X_0}(X)+I_n\) is invertible, we define \[\varphi(X)=C^{-1} \big(O_{X_0}(X)\big) =\big(XX_0^{\mathsf T}-I_n\big)\big(XX_0^{\mathsf T}+I_n\big)^{-1}.\] This provides a local reparameterization of \(\mathrm{St}(n,n)\) around \(X_0\) with \(\varphi(X_0)=\mathbf{0}\).

Since \(\mathcal{T}_{X_t}\mathrm{St}(n,n)=\{W X_t: W\in \mathrm{skew}(n)\}\), the Riemannian gradient can be written as \(\nabla f(X_t)=W_tX_t\) for some \(W_t\in \mathrm{skew}(n)\). The RMD update then in Algorithm 1 admits the explicit form \[\label{eq:rmd95sti} X_{t+1}=\big(I_n+\tfrac{\eta_t}{2}W_t\big)^{-1}\big(I_n-\tfrac{\eta_t}{2}W_t\big)X_t,\tag{2}\] which coincides with the Curvilinear Gradient Descent (CGD) update in [1]. For the general Stiefel manifold \(\mathrm{St}(p,n)\) with \(p<n\), we may extend any \(X\in \mathrm{St}(n,p)\) to an orthogonal matrix by choosing an orthonormal complement \(X_{\perp}\in\mathbb{R}^{n\times(n-p)}\) such that \(\bar X = [\,X,\;X_{\perp}\,]\in \mathrm{St}(n,n).\) Applying the same construction to \(\bar X\) and then taking its first \(p\) columns yields the iteration 2 on \(\mathrm{St}(n,p)\).

3.3 Stochastic Riemannian Mirror Descent↩︎

In this subsection, we introduce a stochastic variant of Riemannian Mirror Descent (RMD). Suppose that, instead of the exact Riemannian gradient \(\nabla f(x)\), we have access to an unbiased stochastic estimator \(\nabla f(x,\xi)\) satisfying \[\mathbb{E}\left[\nabla f(x,\xi)\right]=\nabla f(x),\] where \(\xi\) denotes the randomness in the oracle. In this setting, we can run Algorithm 1 by replacing the exact gradient with its stochastic estimate. In particular, the dual update becomes \[\label{dualupdate95in95RMBD} y_{t+1} = R^{t}{\varphi_t(x_t)}\left(-\eta_t D\varphi_t(x_t)\big[\nabla f(x_t,\xi_t)\big]\right),\tag{3}\] where \(\{\xi_t\}_{t=1}^T\) are independent sampling noises. We refer to this algorithm as Stochastic Riemannian Mirror Descent (SRMD). Analogous to Theorem 1, we establish the following convergence theorem for SRMD.

A detailed proof is given in Section 4.

As an illustration, we specialize SRMD to optimization over the Stiefel manifold \(\mathrm{St}(n,p)\) and obtain a randomized algorithm suitable for the regime where \(p\) is large. Wen and Yin [1] use the Sherman–Morrison–Woodbury (SMW) formula to accelerate the Cayley-type update 2 when \(2p \ll n\), exploiting the fact that the associated skew-symmetric matrix has rank at most \(2p\). When \(2p \geq n\), this low-rank reduction no longer yields a smaller system to invert (and hence provides little or no computational advantage); in this regime, one typically resorts to direct solvers or other implementations; see [1]. Motivated by SRMD, we propose the Stochastic Curvilinear Gradient Descent (SCGD) algorithm, which is amenable to parallel implementation and is effective when \(p\) is large.

We need the following lemma (Lemma 1 in [1]) to express the tangent-space gradient in terms of a skew-symmetric matrix.

The CGD update requires applying \((I_n+\eta_t W_t/2)^{-1}\). When \(2p\ll n\), \(W_t\) has rank at most \(2p\) and the SMW formula can substantially reduce the computational cost. When \(2p\geq n\), this reduction no longer yields a smaller system to invert, and computing \((I_n+\eta_t W_t/2)^{-1}\) can be expensive.

To address the large-\(p\) regime, we construct a block-diagonal approximation by randomly partitioning \([n]\) into \(K\) equal-sized subsets \(\mathcal{I}_1,\ldots,\mathcal{I}_K\) (with \(n=mK\)) and defining \(\hat{W}\) by keeping only the within-block entries: \[\hat{w}_{ij}= \begin{cases} w_{ij}, & i,j\in \mathcal{I}_k \text{ for some } k,\\ 0, & \text{otherwise}. \end{cases}\]For example, if \(\mathcal{I}_k = \{(k-1)m+1,\cdots,km\}\), then \[\hat{W} =\begin{pmatrix} W_{1} & O & O & O\\ O & W_{2} & O & O\\ O & O & \cdots & O\\ O & O & O & W_{K} \end{pmatrix} .\] For any \(i\neq j\), we have \(\mathbb{P}(i,j\text{ belong to the same block})=\frac{m-1}{n-1}\), and hence \[\mathbb{E}[\hat{W}] = \frac{n/K-1}{n-1} W.\] Therefore, \(\frac{n-1}{n/K-1}\hat{W}\) is an unbiased estimator of \(W\). Since \(\hat{W}\) is block diagonal, applying \((I+\eta_t \hat{W}/2)^{-1}\) reduces to inverting \(K\) independent \(m\times m\) systems \((I+\eta_t \hat{W}(\mathcal{I}_k,\mathcal{I}_k)/2)^{-1}\), which can be readily parallelized. This enables the algorithm to handle large-scale problems; see Algorithm 2 for details.

4 Convergence Analysis↩︎

This section is devoted to proving Theorems 1 and  2. To this end, we first establish the following lemma on the difference between retraction and the exponential map. Let \(\Phi_x(u) = \mathrm{Exp}_x^{-1}(R_x(u))\) be well-defined for a retraction \(R_x\).

For \(u\in B_r(x)\), consider the Taylor expansion of \(g( t) =\Phi _{x}( tu)\) at \(t=0\), \[\begin{align} g( 1) =\Phi _{x}( u) & =\Phi _{x}( 0) +D\Phi _{x}( 0)[ u] +\int _{0}^{1}( 1-t) D^{2} \Phi ( tu)[ u,u] dt\\ & =u+\int _{0}^{1}( 1-t) D^{2} \Phi ( tu)[ u,u] dt, \end{align}\] where we use \(R_{x}(0)=x\) and \(D\Phi_{x}(0)=\mathrm{Id}\).

By assumption, there holds \[\Vert \Phi _{x}( u) -u\Vert _{x} \leqslant L_{\Phi }\Vert u\Vert _{x}^{2}\int _{0}^{1} (1-t)dt=\frac{L_{\Phi }}{2}\Vert u\Vert _{x}^{2},\] which implies the desired result. \(\Box\) pt

Using Lemma 2, we can analyze the one-step progress of RMD. Recall that \[\hat{R}_{x}^t( u) =\varphi _{t}^{-1}( R_{\varphi _{t}( x)}^t( D\varphi _{t}( x)[u])).\] Then \(D\hat{R}_{x}^t( 0) =( D\varphi _{t}( x))^{-1} D\varphi _{t}( x) =\mathrm{Id}\), which implies that \(\hat{R}_x^t\) is a retraction.

Let \(d_{t} =-\mathrm{Exp}^{-1}_{x_t}(\hat{R}^t_{x_{t}}( -\eta _{t} \nabla f( x_{t})))\), by Assumption 3, \[\begin{align} f( x_{t+1}) -f( x_{t}) & \leqslant -\langle \nabla f( x_{t}) ,d_{t} \rangle _{x_{t}} +\frac{L_{f}}{2}\Vert d_{t}\Vert _{x_{t}}^{2} \notag\\ & =-\eta _{t}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} -\langle \nabla f( x_{t}) ,d_{t} -\eta _{t} \nabla f( x_{t}) \rangle _{x_{t}} +\frac{L_{f}}{2}\Vert d_{t}\Vert _{x_{t}}^{2}\notag \\ & \leqslant -\left( \eta _{t} -\frac{L_{\Phi } G\eta _{t}^{2}}{2}\right)\Vert \nabla f( x_{t})\Vert _{x_t}^{2} +\frac{L_{f}}{2}\Vert d_{t}\Vert _{x_{t}}^{2}.\label{eq:descent1} \end{align}\tag{4}\] Here we use Cauchy inequality and Assumption 4 to bound \(\langle \nabla f( x_{t}) ,d_{t} -\eta _{t} \nabla f( x_{t}) \rangle _{x_{t}}\) in the second line as \[\begin{align} \langle \nabla f( x_{t}) ,d_{t} -\eta _{t} \nabla f( x_{t}) \rangle _{x_{t}} &\leqslant \Vert \nabla f( x_{t})\Vert _{x_{t}}\Vert d_{t} -\eta _{t} \nabla f( x_{t})\Vert _{x_{t}} \\ &\leqslant G\Vert d_{t} -\eta _{t} \nabla f( x_{t})\Vert _{x_{t}} \\ &\leqslant \frac{L_{\Phi } G\eta _{t}^{2}}{2}\Vert \nabla f( x_{t})\Vert _{x_t}^{2}. \end{align}\] By Lemma 2, \[\begin{align} \Vert d_{t}\Vert _{x_{t}}^{2} &=\Vert d_{t} -\eta _{t} \nabla f( x_{t}) +\eta _{t} \nabla f( x_{t})\Vert _{x_{t}}^{2}\\ & \leqslant 2\Vert d_{t} -\eta _{t} \nabla f( x_{t})\Vert _{x_{t}}^{2} +2\Vert \eta _{t} \nabla f( x_{t})\Vert _{x_{t}}^{2}\\ & \leqslant \frac{L_{\Phi }^{2} \eta _{t}^{4}}{2}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{4} +2\eta _{t}^{2}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2}. \end{align}\] Together with (4 ) we obtain \[f( x_{t+1}) -f( x_{t}) \leqslant -\left( \eta _{t} -\frac{L_{\Phi } G\eta _{t}^{2}}{2} -\eta _{t}^{2} L_{f}\right)\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} +\frac{L_{f} L_{\Phi }^{2} \eta _{t}^{4}}{4}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{4} .\] We thus conclude the proof. \(\Box\) pt

Now we are ready to prove Theorem 1.

[Proof of Theorem 1] If \(\eta < \min\left\{1/\left( L_{\Phi } G/2+L_{f} +L_{f} L_{\Phi }^{2}\right) ,2/G, r/G\right\}\), we have \[\eta ^{4}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{4} \leqslant 4\eta ^{2}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} .\] By Lemma 3, \[f( x_{t+1}) -f( x_{t}) \leqslant -\left( \eta -L_{\Phi } G\eta ^{2} /2-L_{f} \eta ^{2} -L_{f} L_{\Phi }^{2} \eta ^{2}\right)\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} .\] Let \(C_{d} =\eta -L_{\Phi } G\eta ^{2} /2-L_{f} \eta ^{2} -L_{f} L_{\Phi }^{2} \eta ^{2} >0.\)

Nonconvex Case: Summing up over \(\displaystyle t\) and rearranging, we get \[\begin{align} \frac{1}{T}\sum _{t=1}^{T} \| \nabla f(x_{t} )\| _{x_{t}}^{2} & \leqslant \frac{f( x_{1}) -f( x_{T+1})}{C_{d} T} . \end{align}\]

Geodesically convex case: We construct a potential function to prove the convergence rate. Let \(A_{t+1} =A_{t} +a_{t} ,A_{0} =0,\;E_{t} =A_{t}\left( f( x_{t}) -f\left( x^{*}\right)\right),\) where \(a_{t}\) will be specified later. We have \[\begin{align} E_{t+1} -E_{t} & =A_{t+1}\left( f( x_{t+1}) -f\left( x^{*}\right)\right) -(A_{t+1}-a_t)\left( f( x_{t}) -f\left( x^{*}\right)\right)\notag \\ & =A_{t+1}( f( x_{t+1}) -f( x_{t})) + a_{t}\left( f( x_{t}) -f\left( x^{*}\right)\right) .\label{eq:potential95descent1} \end{align}\tag{5}\] By the definition of convexity, \(f( x_{t}) -f\left( x^{*}\right) \leqslant -\langle \nabla f(x_{t}),\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right) \rangle _{x_{t}}\). Combining with (5 ), we have \[\begin{align} E_{t+1} -E_{t} & \leqslant -C_{d} A_{t+1}\Vert \nabla f( x_{t})\Vert_{x_t} ^{2} -a_{t} \langle \nabla f( x_{t}), \mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right) \rangle _{x_{t}}\\ & =-C_{d} A_{t+1}\left(\Vert \nabla f( x_{t})\Vert_{x_t} ^{2} + \langle \nabla f( x_{t}) ,\frac{a_{t}}{C_{d} A_{t+1}}\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right) \rangle_{x_t} \right) . \end{align}\] Let \(H( v) =\Vert v\Vert _{x_{t}}^{2} + \langle v,\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right) /C\rangle _{x_{t}}\), where \(C\) is a constant. Since \(H\) is a convex function, it attains its minimal at stationary point \(v^* = -\frac{\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right)}{2C}\). \[H( v) \geqslant H\left(\frac{\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right)}{2C}\right) =-\frac{1}{4C^{2}}\left\Vert \mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right)\right\Vert _{x_{t}}^{2} .\] By this equation, \[E_{t+1} -E_{t} \leqslant \frac{a_{t}^{2}}{4C_{d} A_{t+1}}\left\Vert \mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right)\right\Vert _{x_{t}}^{2} .\label{eq:potential95descent2}\tag{6}\] Since we have \(\left\Vert \mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right)\right\Vert_{x_t} =d\left( x_{t} ,x^{*}\right) \leqslant \text{diam}( A)\), plug it into (6 ), \[E_{t+1} -E_{t} \leqslant \frac{a_{t}^{2}\text{diam}( A)^{2}}{4C_{d} A_{t+1}} .\] Summation with respect to \(t\), we get \[E_{T} \leqslant E_{0} +\frac{\text{diam}( A)^{2}}{4C_{d}}\sum _{t=0}^{T-1}\frac{a_{t}^{2}}{A_{t+1}} .\] Let \(A_{t} = t^{2} ,a_{t} =A_{t+1} -A_{t} =2t+1\), we have \[\begin{align} E_{T} &\leqslant E_{0} +\frac{\text{diam}( A)^{2}}{4C_{d}}\sum _{t=0}^{T-1}\frac{( 2t+1)^{2}}{( t+1)^{2}} \\ &< E_{0} +\frac{\text{diam}( A)^{2}}{4C_{d}}\sum _{t=0}^{T-1}\frac{( 2t+2)^{2}}{( t+1)^{2}} \\ &=E_{0} +\frac{\text{diam}( A)^{2} T}{C_{d}} . \end{align}\] By the definition of \(E\), \(E_{T} =A_{T}\left( f( x_{T}) -f\left( x^{*}\right)\right)\), \[f( x_{T}) -f\left( x^{*}\right) \leqslant \frac{E_{0}}{T^{2}} +\frac{\text{diam}( A)^{2}}{C_{d} T} =\mathcal{O}\left(\frac{1}{\eta T}\right),\] where the \(\eta\) factor in the denominator comes from \(C_d\). This concludes the proof of Theorem 1. \(\Box\) pt

To establish the convergence guarantee for stochastic RMD, i.e., Theorem 2, we first prove the following lemma. Its role is analogous to Lemma 3.

Let \(d_{t} =-\mathrm{Exp}^{}_{x_{t}}(\hat{R}^t_{x_{t}}( -\eta _{t} \nabla f( x_{t},\xi_t)))\), by Assumption 3, \[\begin{align} \mathbb{E}_{t} f( x_{t+1}) -f( x_{t}) & \leqslant -\mathbb{E}_{t}\langle \nabla f( x_{t}) ,d_{t} \rangle _{x_{t}} +\frac{L_{f}}{2}\mathbb{E}_{t}\Vert d_{t}\Vert _{x_{t}}^{2} \notag\\ & =-\eta _{t}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} -\mathbb{E}_{t} \langle \nabla f( x_{t}) ,d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t}) \rangle _{x_{t}} +\frac{L_{f}}{2}\mathbb{E}_{t}\Vert d_{t}\Vert _{x_{t}}^{2} \notag\\ & \leqslant -\eta _{t}\Vert \nabla f( x_{t})\Vert _{x_t}^{2} +\frac{L_{\Phi } G\eta _{t}^{2}}{2}\mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} +\frac{L_{f}}{2}\mathbb{E}_{t}\Vert d_{t}\Vert _{x_{t}}^{2}, \label{eq:stoc95d} \end{align}\tag{7}\] where in the last inequality we use \[\begin{align} \langle \nabla f( x_{t}) ,d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t}) \rangle _{x_{t}} & \leqslant \Vert \nabla f( x_{t})\Vert _{x_{t}}\Vert d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}\\ & \leqslant G\Vert d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}} \\ &\leqslant \frac{GL_{\Phi } \eta_t^2}{2}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} , \end{align}\] and Lemma 2, \[\begin{align} \Vert d_{t}\Vert _{x_{t}}^{2} & =\Vert d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t}) +\eta _{t} \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2}\\ & \leqslant 2\Vert d_{t} -\eta _{t} \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} +2\Vert \eta _{t} \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2}\\ & \leqslant \frac{L_{\Phi }^{2} \eta _{t}^{4}}{2}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{4} +2\eta _{t}^{2}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2}. \end{align}\] Together with (7 ), we get \[\begin{align} \mathbb{E}_{t}f( x_{t+1}) -f( x_{t}) \leqslant -\eta _{t}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} +\left(\frac{L_{\Phi } G\eta _{t}^{2}}{2} +\eta _{t}^{2} L_{f}\right)\mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} \\ \quad\quad\quad\quad\quad+\frac{L_{f} L_{\Phi }^{2} \eta _{t}^{4}}{4}\mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{4}, \end{align}\] which concludes the proof. \(\Box\) pt

Now we are ready to prove Theorem 2.

[Proof of Theorem 2] Let \(\displaystyle C_{sd} =\left( L_{\Phi } G+2L_{f} +2L_{f} L_{\Phi }^{2}\right)\), if \(\eta _{t} = \eta < \min\{1/C_{sd},2/G, r/G\}\), we have \[\begin{align} \mathbb{E}_t f( x_{t+1}) -f( x_{t}) & \leqslant -\eta\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} +\left(\frac{L_{\Phi } G\eta^{2}}{2} +\eta^{2} L_{f}\right)\mathbb{E}_t\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} \\ & \quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{L_{f} L_{\Phi }^{2} \eta^{4}}{4}\mathbb{E}_t\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{4}\\ & \leqslant -\eta\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} + C_{sd} \eta^{2}\mathbb{E}_t\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2}, \end{align}\] where we use \(\eta ^{4}\Vert \nabla f( x_{t},\xi_{t})\Vert _{x_{t}}^{4} \leqslant 4\eta ^{2}\Vert \nabla f( x_{t}, \xi_{t})\Vert _{x_{t}}^{2}\) since by assumption, \(\Vert \nabla f( x_{t}, \xi_{t})\Vert _{x_{t}}^{2} \leqslant G\). For the second term, we have \[\begin{align} \mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t})\Vert _{x_{t}}^{2} & =\mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t}) -\nabla f( x_{t}) +\nabla f( x_{t})\Vert _{x_{t}}^{2}\\ & \leqslant 2\mathbb{E}_{t}\Vert \nabla f( x_{t} ,\xi _{t}) -\nabla f( x_{t})\Vert _{x_{t}}^{2} + 2\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2}\\ & \leqslant 2\sigma ^{2} +2\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} . \end{align}\] Thus, we have \[\mathbb{E}_tf( x_{t+1}) -f( x_{t}) \leqslant -\left( \eta -C_{sd} \eta^{2}\right)\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2} +C_{sd} \eta ^{2} \sigma ^{2}.\] Let \(F_t = \mathbb{E}f(x_{t}) - f(x^*) - {t}C_{sd} \eta ^{2} \sigma ^{2}\), then \[F_{t+1} - F_t \leqslant -\left( \eta -C_{sd} \eta^{2}\right) \mathbb{E}\Vert \nabla f( x_{t})\Vert _{x_{t}}^{2}.\] Using the same argument as the proof of Theorem 1, we will prove the theorem in different cases.

Nonconvex Case: Summing up over \(t\) and rearranging leads to \[\begin{align} \frac{1}{T}\sum _{t=1}^{T}\mathbb{E} \| \nabla f(x_{t} )\| _{x_{t}}^{2} & \leqslant \frac{F_{1} -F_{T+1}}{\left( \eta -C_{sd} \eta ^{2}\right) T} \\ & \leqslant \frac{f( x_{1}) -f( x_{*})}{\left( \eta -C_{sd} \eta ^{2}\right) T} +\frac{C_{sd} \eta \sigma ^{2}}{( 1-C_{sd} \eta )}. \end{align}\]

Geodesically convex case: Let \(\displaystyle A_{t} =t^{2} ,\;a_{t} =A_{t} -A_{t-1} =2t+1\). We define \(\displaystyle E_{t} =A_{t}\left(\mathbb{E} f( x_{t}) -f\left( x^{*}\right)\right)\). Thus, \[\begin{align} E_{t+1} -E_{t} & =A_{t+1}\mathbb{E}[ f( x_{t+1}) -f( x_{t})] +a_{t}\mathbb{E}\left[ f( x_{t}) -f\left( x^{*}\right)\right] . \end{align}\] By Lemma 4 and convexity, we have \[\begin{align} E_{t+1} -E_{t} & \leqslant -A_{t+1}\left( \eta -C_{sd} \eta ^{2}\right)\mathbb{E} \| \nabla f( x_{t}) \|_{x_t}^{2} -a_{t} \langle \nabla f(x_{t} ),\mathrm{Exp}^{-1}_{x_{t}}\left( x^{*}\right) \rangle _{x_{t}} +A_{t+1} C_{sd} \eta ^{2} \sigma ^{2}\\ & \leqslant \frac{a_{t}^{2}}{4A_{t+1}\left( \eta -C_{sd} \eta ^{2}\right)}\mathbb{E} \| \nabla f( x_{t}) \|_{x_t}^{2} +A_{t+1} C_{sd} \eta ^{2} \sigma ^{2}\\ & \leqslant \frac{a_{t}^{2}\text{diam}( A)}{4A_{t+1}\left( \eta -C_{sd} \eta ^{2}\right)} +A_{t+1} C_{sd} \eta ^{2} \sigma ^{2} \end{align}\]where we use a similar argument as the proof of Theorem 1. Therefore, \[\begin{align} E_{T} & \leqslant E_{0} +\frac{\text{diam}( A)}{4\left( \eta -C_{sd} \eta ^{2}\right)}\sum _{t=0}^{T-1}\frac{( 2t+1)^{2}}{( t+1)^{2}} +\sum _{t=0}^{T-1}( t+1)^{2} C_{sd} \eta ^{2} \sigma ^{2}\\ & \leqslant O\left( E_{0} +\frac{\text{diam}( A) T}{\eta -C_{sd} \eta ^{2}} +C_{sd} \eta ^{2} \sigma T^{3}\right) . \end{align}\] Note that \(\displaystyle E_{T} =T^{2}\left(\mathbb{E} f( x_{T}) -f\left( x^{*}\right)\right)\), we get \[\begin{align} \mathbb{E} f( x_{T}) -f\left( x^{*}\right) & \leqslant O\left(\frac{1}{\eta T} +\eta ^{2} \sigma T\right). \end{align}\] We thus finish the proof of Theorem 2. \(\Box\) pt

5 Numerical Experiments↩︎

In this section, we do some experiments on Stiefel manifold optimization to verify our algorithm. We compare our algorithm with the classic algorithm [1], which is a special case of Algorithm 2. Since OptStiefel adapts non-monotonic line search strategy [30], which can greatly improve the performance of the algorithm, we adopt the same line search technique [30] to Algorithm 2. In [1], the authors also use the Barzilai-Borwein (BB) step size heuristic to accelerate their algorithm. However, it is not clear how to use BB step size in the stochastic setting. To highlight the improvement made by randomization, we do not use BB step size in our experiment. We set the stop criterion to be \[\| \nabla_\mathcal{M}f(x) \| \leqslant 10^{-5} \;\text{or} \quad \|x_t - x_{t-1}\|\leqslant 10^{-5} \;\text{or}\quad \|f(x_t) - f(x_{t-1})\| \leqslant 10^{-8}.\] All the experiments are tested on a machine with Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz 2.59 GHz and 16.0GB RAM².

5.1 Linear Eigenvalue Problem↩︎

Given a symmetric matrix \(A\), the linear eigenvalue problem is \[\begin{align} \max_{X\in\mathrm{St}(n,p)} \text{Tr}(X^\top A X) = \sum_{k=1}^p \lambda_k, \end{align}\] where \(\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_p\) are the \(p\)th largest eigenvalues of \(A\). In the experiments, we generate a random matrix \(N\) whose entries are i.i.d. standard normal random variables and set \(A=N^\top N\). The test results are reported in Table 1. We observe that while both algorithms produce solutions of similar quality, SCGD achieves a lower runtime than CGD. One reason is that SCGD only needs to solve several smaller linear systems rather than a large one in each iteration.

5.2 Orthogonal Procrustes Problem↩︎

The orthogonal procrustes problem [31] is \[\begin{align} \min_{X\in\mathrm{St}(n,p)} & \|AX-B\|_F^2. \end{align}\] In the experiments, we generate \(A=(A_{ij})\) with \(A_{ij}\overset{\text{i.i.d.}}{\sim}\mathrm{Unif}(0,1)\), and construct \(X^\star\in\mathrm{St}(n,p)\) by taking the \(Q\) factor of the QR decomposition of a Gaussian random matrix. We then set \(B=AX^\star\). Since the objective is nonnegative and \(X^\star\) attains value \(0\), the optimal value is \(0\). The results are reported in Table 2.

6 Conclusion↩︎

In this paper, we propose a new framework for mirror descent (MD) on Riemannian manifolds, based on the key observation that MD can be interpreted as an optimization method under a suitable reparameterization. Under mild assumptions, we obtain the first non-asymptotic convergence result of Riemannian Mirror Descent in both the deterministic and stochastic settings. Furthermore, we introduce the Stochastic Curvilinear Gradient Descent algorithm for large-scale Stiefel manifold optimization under the RMD framework, which we believe is of independent interest. Similar techniques may be applied to address optimization problems in other settings.

An interesting direction for future work is the explicit construction of problem-adapted reparameterizations. For instance, Li et al. [23] showed that, under any commuting parametrization, the gradient flow is equivalent to the continuous-time mirror flow. In our general Riemannian setting, where no additional global structure is imposed, the resulting reparameterizations are inherently local. It would be of interest to investigate whether stronger geometric assumptions, such as those satisfied by Hessian manifolds, permit global reparameterizations and lead to sharper guarantees and improved algorithms.

References↩︎

---

1. Authors are listed in alphabetical order and contributed equally to this work. The corresponding author is Lei Shi. The work of Lei Shi is supported by the National Natural Science Foundation of China [Grant No.12171093]. Email addresses: jxjiang20@fudan.edu.cn (J. Jiang), leishi@fudan.edu.cn (L. Shi), jiyuantan19@gmail.com (J. Tan).↩︎

2. The code and data set are available at https://github.com/JiyuanTan/RMD/tree/main.↩︎