The Intrinsic Riemannian Proximal Gradient Method for Nonconvex Optimization
June 11, 2025
Abstract
We consider the proximal gradient method on Riemannian manifolds for functions that are possibly not geodesically convex. Starting from the forward-backward-splitting, we define an intrinsic variant of the proximal gradient method that uses proximal maps defined on the manifold and therefore does not require or work in any embedding. We investigate its convergence properties and illustrate its numerical performance, particularly for nonconvex or nonembedded problems that are hence out of reach for other methods.
1 Introduction↩︎
Nonsmooth optimization on Riemannian manifolds has recently gained considerable interest. Starting with inaugural algorithms like the subgradient method FerreiraOliveira:1998:1? and the proximal point algorithm FerreiraOliveira:2002:1?, introducing more properties of either the objective or the manifold has led to the book Bacak:2014:2?. For methods involving a bundle or approaches that involve a bundle of subgradients, see HoseiniMonjeziNobakhtianPouryayevali:2021:1?, BergmannHerzogJasa:2024:1?. Especially the proximal bundle method from HoseiniMonjeziNobakhtianPouryayevali:2021:1? is one of the first methods to investigate the non-convex case.
A prominent approach to nonsmooth optimization are splitting methods, where the objective function can be split. For a Riemannian manifold \(\cM\), we consider the problem \[\label{eq:splitting} \argmin_{p \in \cM} f(p) , \qquad \text{where} \quad f(p) = g(p) + h(p) .\tag{1}\] Some, like the cyclic proximal point algorithm (CPPA) Bacak:2014:1? involve an arbitrary number of summands, while others also work for various splittings, , convex-concave-splittings or difference of convex (DC) problems [1]–[3]. In the case with two functions, as in 1 , the goal is to exploit properties of \(g\) and \(h\) that allow for an efficient minimization. For example when both have “cheap” proximal maps, the Douglas-Rachford algorithm BergmannPerschSteidl:2016:1? can be used also on Hadamard manifolds. Another prominent algorithm is the so-called Alternating Direction Method of Multipliers (ADMM) [4], for which variants that work on manifolds were, , considered in [5], [6]. Finally, together with an investigation of Fenchel duality, [7], SilvaLouzeiroBergmannHerzog:2022:1? allows for the so-called Chambolle-Pock algorithm BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez:2021:1? to be defined on manifolds, also involving proximal maps. A prominent algorithm that is based on splitting is the proximal-gradient method, see Beck:2017:1? for a comprehensive overview in the Euclidean case. In this case, \(g\) is assumed to be (\(L_g\)-)smooth and \(h\) to be convex and potentially nonsmooth. If \(h\) has a proximal map that is fast to evaluate, maybe even in closed form, the algorithm is efficient, and a convergence rate can be derived.
To the best of our knowledge, the authors in [8] first introduced a partially convex proximal-gradient version for the Stiefel manifold, phrasing the nonsmooth part as a constrained problem in the embedding, and requiring the nonsmooth function to be the restriction of a convex function in the embedding. This was generalized to general Riemannian manifolds in [9]. In this work, the authors also phrase the nonsmooth subproblem as a problem in the tangent space at the current iterate, but forgo convexity assumptions in their general analysis. The same authors generalize their framework to an inexact version in [10] and still keep their convergence analysis from the previous paper. The linear convergence of [8] and [9] under strong convexity has been studied in [11].
The authors in [12] rephrase the nonsmooth subproblem as a trust-region problem. All these works are based on a “generalized notion of the proximal map” [9], and solving the corresponding subproblem either in the embedding or the tangent space of the current iterate. Most of these works assume \(g\) and \(h\) fulfil some convexity property, either along geodesics or along curves generated by retractions when using a retraction-based generalization of the proximal subproblem. In the work of [9], a “retraction smoothness” of \(g\) is required, but \(h\) can be nonsmooth and nonconvex for their convergence investigations. However, their algorithm employs the embedding in a way that \(h\) is required to have certain properties thereon in order to solve their approximation to the proximal map. In [13] the authors do employ a proximal map, but restrict their investigation to Hadamard manifolds and employing the Kurdyka-Łojasiewicz property on manifolds. In a related work, [14] consider gradient schemes that combine approximate solutions to proximal maps defined on the manifold, leading to a globally accelerated gradient descent algorithm on Hadamard manifolds.
1.1 Contributions↩︎
In this work, we introduce a formulation of the Riemannian proximal gradient method that actually employs proximal maps for the nonsmooth step. This relates the new proximal gradient method directly back to both the PPA FerreiraOliveira:2002:1? as well as gradient descent. This new algorithm can still be generalized to a corresponding retraction-based proximal map, which, similarly to [9], can be understood as an approximate variant. This new formulation allows to investigate several properties from Beck:2017:1? on manifolds, for the case where neither the smooth part \(g\) nor the nonsmooth part \(h\) are convex. The aim is to derive a convergence theory that is in accordance with [9] also for our intrinsic variant, that is, it is also based on intrinsic properties that \(g\) and \(h\) have as functions defined on the manifold. A major difference to [9] is that the method works intrinsically, entirely on the manifold, without solving a subproblem on a tangent space or requiring an embedding. This especially allows using Riemannian proximal maps defined intrinsically on the manifold, such as the one introduced in, , FerreiraOliveira:2002:1?, instead of Euclidean-flavored approximations thereof. Furthermore, this approach can directly be applied also to manifolds that are not submanifolds of a Euclidean space. We show two main results:
In 2, we show sublinear convergence to \(\varepsilon\)-stationary points, with rate \(\cO(1/\sqrt k)\).
In 3, we show that any accumulation point \(p^\ast\) of the sequence of iterates is a stationary point of \(f\), , \(0_{p^\ast} \in \partial f(p^\ast)\).
Crucially, our analysis does not make any assumption on the convexity of either of the summands \(g, h\) that the objective is split into, whereas a more standard approach usually involves assuming the nonsmooth part \(h\) to be convex. As a result, our framework captures manifolds of any bounded curvature, be it positive or negative. Finally, we validate the performance of our method on realistic numerical examples. We show that in many cases, the proximal map of \(h\) either exists in closed-form or is efficiently computable intrinsically on the manifold.
1.2 Organization↩︎
The remainder of this paper is organized as follows. After introducing the necessary notation and some preliminary properties in 2, we derive the new formulation of the nonconvex proximal gradient method in 3. In 4, we discuss convergence properties of our proposed algorithm. We introduce a generalization of the algorithm that uses retractions and employs a retraction based distance and proximal map in 5. 6 illustrates the performance of the algorithm and its convergence, and 7 concludes the paper.
2 Preliminaries↩︎
In this section, we recall some notions from Riemannian geometry. For more details, the reader may wish to consult DoCarmo:1992:1?, Boumal:2023:1?. Let \(\cM\) be a smooth manifold with Riemannian metric \(\riemannian{\cdot}{\cdot}\) and let \(\covariantDerivativeSymbol\) be its Levi-Civita connection. Given \(p, q \in \cM\), we denote by \(C_{pq}\) the set of all piecewise smooth curves \(\gamma \colon[0,1] \to \cM\) from \(p = \gamma(0)\) to \(q = \gamma(1)\). The arc length of \(\gamma \in C_{pq}\), \(L \colon C_{pq} \to \bbR_{\ge 0}\) is defined as \({ L(\gamma) = \int_0^1 \riemanniannorm{\dot{\gamma}}[\gamma(t)] \, \d t , }\) where \(\riemanniannorm{\cdot}[p]\) is the norm induced by the Riemannian metric in the tangent space \(\tangentSpace{p}\) at the point \(p\). The inner product given by the Riemannian metric at \(p\) is \({\riemannian{\cdot}{\cdot}[p] \colon \tangentSpace{p} \times \tangentSpace{p} \to \bbR}\). The explicit reference to the base point will be omitted whenever the base point is clear from context. The sectional curvature at \(p \in \cM\) is denoted by \(K_p(X_p,Y_p)\) for tangent vectors \(X_p, Y_p \in \tangentSpace{p}\).
The Riemannian distance between two points, denoted with \(\dist \colon \cM \times \cM \to \bbR_{\ge 0} \cup \{+\infty\}\), is given by \(\dist(p,q) = \inf_{\gamma \in C_{pq}} L(\gamma) .\) A geodesic arc from \(p \in \cM\) to \(q \in \cM\) is a smooth curve \(\gamma \in C_{pq}\) that is parallel along itself, , \(\covariantDerivative{\dot{\gamma}(t)}[\dot{\gamma}(t)] = 0\) for all \(t \in (0,1)\). A finite-length geodesic arc is minimal if its arc length coincides with the Riemannian distance between its extreme points. Given two points \(p, q \in \cM\), we denote with \(\geodesic{p}{q}\) a minimal geodesic arc that connects \(p = \geodesic{p}{q}(0)\) to \(q = \geodesic{p}{q}(1)\). Given a point \(p\in \cM\) and a tangent vector \(X_p \in \tangentSpace{p}\), a geodesic can also be specified as the curve \(\geodesic{p}{X_p}\) such that \(\geodesic{p}{X_p}(0)=p\) and \(\geodesic{p}{X_p}'(0) = X_p\). The exponential map at \(p \in \cM\) maps a tangent vector \(X_p \in \tangentSpace{p}\) to the endpoint of the corresponding geodesic, \(\exponential{p}(X_p) = \geodesic{p}{X_p}(1)\). If it exists, the inverse map is called the logarithmic map \(\logarithm{p}\colon \cM \rightarrow \tangentSpace{p}\). We call a set uniquely geodesic if every two points in that space are connected by one and only one geodesic in that set. In this case the logarithmic map exists everywhere. A set \(\cU \subseteq \cM\) is said to be (uniquely) geodesically convex if for any \(p,q \in \cU\), there exists a (unique) minimal geodesic arc \(\geodesic{p}{q}\) which lies entirely in \(\cU\). Observe that if a manifold \(\cM\) has positive sectional curvature, the set \(\cU\) may not be allowed to have arbitrarily large diameter. For example, if \(\kmin > 0\), for \(\cU\) to be uniquely geodesically convex, its diameter must be bounded by \(\diam(\cU) < \pi/\sqrt{\kmin}\). See, , [15], WangLiYao:2016:1? or BridsonHaefliger:1999:1?.
We denote the parallel transport from \(p\) to \(q\) along \(\geodesic{p}{q}\) with respect to the connection \(\covariantDerivativeSymbol\) with \(\parallelTransport{p}{q} \colon \tangentSpace{p} \to \tangentSpace{q}\), defined by \(\parallelTransport{p}{q} Y_p = X_q\) for \(Y_p \in \tangentSpace{p} ,\) where \(X\) is the unique smooth vector field along \(\geodesic{p}{q}\) which satisfies \(X_p = Y_p \text{ and } \covariantDerivative{\dot{\gamma}}{X} = 0 .\) Notice that \(\parallelTransport{p}{q}\) depends on the choice of \(\geodesic{p}{q}\).
Given a differentiable function \(f \colon \cM \to \bbR\), we denote the Riemannian gradient of \(f\) at \(p \in \cM\) by \(\grad f(p) \in \tangentSpace{p}\). See Boumal:2023:1? for more details. Given a uniquely geodesic set \(\cU \subseteq \cM\) and a real number \(L_f > 0\), the function \(f\) is said to be \(L_f\)-smooth over \(\cU\) if, for any \(p, q \in \cU\), \[\riemanniannorm[auto] {\grad f(p) - \parallelTransport{q}{p} \grad f(q)} \le L_f \dist(p, q) ;\] see, , ZhangSra:2016:1?, Boumal:2023:1?. Furthermore, it can be shown that, if \(f\) is \(L_f\)-smooth over \(\cU\), then \[f(q) \le f(p) + \inner{\grad f(p)}{\logarithm{p}(q)} + \frac{L_f}{2} \dist^2(p, q) ,\] for any \(p, q \in \cU\); Boumal:2023:1?. Given a real number \(L \ge 0\), a function \(f\) is \(L\)-Lipschitz if, for any \(p, q \in \cM\), \[\abs{f(p) - f(q)} \le L \dist(p, q) .\] Furthermore, \(f\) is locally Lipschitz if, for every \(p \in \cM\), there exists a neighborhood \(\cU\) of \(p\) such that \(\restr{f}{\cU}\) is Lipschitz. Let \(f \colon \cM \to \bbR\) be a locally Lipschitz function, and let \(\hat{f} \coloneq f \circ \expOp \colon \tangentBundle \to \bbR\). Then, the generalized Clarke directional derivative of \(f\) in the direction \(X_p \in \tangentSpace{p}\) at \(p \in \cM\) is defined as \[\label{eq:generalized-clarke-directional-derivative} f^\circ(p, X_p) \coloneq \hat{f}^\circ(0_p, X_p) \coloneq \limsup_{\substack{Y_p \to 0_p \\ t \downarrow 0}} \frac{ f \circ \exponential{p}(Y_p + tX_p) - f \circ \exponential{p}(0_p) }{t} .\tag{2}\] The generalized Clarke subdifferential of \(f\) at \(p \in \cM\) is then defined as \[\label{eq:generalized-clarke-subdifferential} \partial f(p) \coloneq \setDef[auto]{X_p \in \tangentSpace{p}}{ f^\circ(p, Y_p) \ge \riemannian{X_p}{Y_p} \; \text{for all } Y_p \in \tangentSpace{p} } .\tag{3}\] Let \(f\) be a locally Lipschitz function on \(\cM\) and \(p^\ast \in \cM\) a local minimizer of \(f\), then \[0_{p^\ast} \in \partial f(p^\ast) ,\] see, , YangZhangSong:2014:1?. If now \(f = g + h\) with \(g\) of class \(\cC^1\) and \(h\) locally Lipschitz, then \[\begin{gather} f^\circ(p, X_p) = \limsup_{\substack{Y_p \to 0_p \\ t \downarrow 0}} \frac{1}{t} \bigl( g \circ \exponential{p}(Y_p + tX_p) - g \circ \exponential{p}(0_p) \\ + h \circ \exponential{p}(Y_p + tX_p) - h \circ \exponential{p}(0_p) \bigr) , \end{gather}\] and since \(\limsup(a_n + b_n) \le \limsup a_n + \limsup b_n\), we get \[\label{eq:directional-derivative-of-sum-inequality-1} f^\circ(p, X_p) = (g + h)^\circ(p, X_p) \le g'(p, X_p) + h^\circ(p, X_p) ,\tag{4}\] where \(g'(p, X_p) = \riemannian{\grad g(p)}{X_p}\). Since \(h^\circ(p, X_p) = ((h + g) - g)^\circ(p, X_p)\), applying 4 , we get \[\label{eq:directional-derivative-of-sum-inequality-2} h^\circ(p, X_p) + g'(p, X_p) \le (h + g)^\circ(p, X_p) .\tag{5}\] 4 and 5 now imply \[\label{eq:directional-derivative-of-sum-equality} (g + h)^\circ(p, X_p) = g'(p, X_p) + h^\circ(p, X_p) .\tag{6}\] This and the definition of the Clarke subdifferential now give \[\label{eq:clarke-subdifferential-of-sum} \partial (g + h) (p) = \grad g(p) + \partial h(p) \quad \text{for all } p \in \cM .\tag{7}\] This is a non-Hadamard version of BentoFerreiraOliveira:2015:1?. Therefore, \(0 \in \partial(g+h)(p) = \grad g(p) + \partial h(p)\) is equivalent to the stationarity condition \[\label{eq:stationarity-condition} - \grad g(p) \in \partial h(p) .\tag{8}\]
Given a function \(f\) and a real number \(\lambda > 0\), we introduce the proximal map as \[\prox{\lambda f}(p) \coloneq \argmin_{q \in \cM} f(q) + \frac{1}{2 \lambda} \dist^2(p, q) ,\] if the minimizer exists. Note that, as opposed to other works, the minimum is taken on the manifold.
In the analysis of the proximal gradient method, we will also need the following curvature-dependent quantities. Given real numbers \(\kappa_1, \kappa_2\), which will typically be lower and upper bounds to the sectional curvature of \(\cM\), and \(s \in \bbR\), we define the following quantities, see also LezcanoCasado:2020:1?: \[\begin{align} \tag{9} \newtarget{def:zeta-first}{ \zeta_{1, \kappa_1}(s) } & \coloneq \begin{cases} 1 & \text{if } \kappa_1 \ge 0 , \\ \sqrt{- \kappa_1} \, s \coth(\sqrt{- \kappa_1} \, s) \quad & \text{if } \kappa_1 < 0 , \end{cases} \\ \tag{10} \newtarget{def:zeta-second}{ \zeta_{2, \kappa_2}(s) } & \coloneq \begin{cases} 1 & \text{if } \kappa_2 \le 0 , \\ \sqrt{\kappa_2} \, s \cot(\sqrt{\kappa_2} \, s) \quad & \text{if } \kappa_2 > 0 , \end{cases}\\ \tag{11} \newtarget{def:sigma-curvature}{ \sigma_{\kappa_1, \kappa_2}(s) } &\coloneq \max \paren[big]\{\}{ \zeta_{1, \kappa_1}(s) , \; \abs{\zeta_{2, \kappa_2}(s)} } ,\\ \tag{12} \newtarget{def:generalized-sine-positive-zeros}{ \pikappa } &\coloneq \begin{cases} \infty & \text{if } \kappa \le 0 , \\ \frac{\pi}{\sqrt{\kappa}} & \text{if } \kappa > 0 . \end{cases} \end{align}\]
3 The Riemannian Proximal Gradient Method↩︎
To derive the proximal gradient method, the function \(g\) is replaced by a proximal regularization of its linearization; see BeckTeboulle:2009:1? for the derivation in the Euclidean case. For \(\lambda > 0\) and a point \(p\) we compute a new candidate as \[\argmin_{q\in\cM} g(p) + \riemannian{\grad g(p)}{\logarithm{p}q} + \frac{1}{2\lambda} \dist^2(p,q) + h(q) .\] We introduce \(r(p) \coloneq \frac{1}{2\lambda} \dist^2(p,q)\) and Taylor expand this function ( Boumal:2023:1?): \[\begin{align} r(\exponential{p}(X_p)) & = r(p) + \riemannian{\grad r(p)}{X_p} + \cO(\riemanniannorm{X_p}^2) , \\ & = r(p) - \frac{1}{\lambda} \riemannian{\logarithm{p}q}{X_p} + \cO(\riemanniannorm{X_p}^2 ) , \end{align}\] where \(X_p \in \tangentSpace{p}\). For the case \(X_p = -\lambda\grad g(p)\) this reads \[\frac{1}{2\lambda} \dist^2( \exponential{p}(-\lambda\grad g(p)) , q ) = \frac{1}{2\lambda} \dist^2(p, q) + \riemannian{\logarithm{p}q}{\grad g(p)} + \cO(\lambda^2) ,\] using this as a second quadratic approximation we obtain the new candidate as \[\argmin_{q\in\cM} g(p) + \frac{1}{2\lambda} \dist^2( \exponential{p}(-\lambda\grad g(p)) , q ) + h(q) .\] Since \(g(p)\) does not change the minimizer, the resulting expression is equal to the proximal map \[\prox[big]{\lambda h}( \exponential{p}(-\lambda\grad g(p)) ) = \argmin_{q\in\cM} \frac{1}{2\lambda} \dist^2( \exponential{p}(-\lambda\grad g(p)) , q ) + h(q) .\]
In this article, we consider the general case where the function \(h\) or the Riemannian distance are not necessarily geodesically convex. Then, the proximal map \(\prox{\lambda h}\) might not be single-valued, and the iteration might not be well-defined. In this case, inspired by [9], we only demand finding a stationary point \(p^\ast\) of the function \[\label{eq:nonconvex-subproblem-objective} H(q; p, \lambda) \coloneq h(q) + \frac{1}{2 \lambda} \dist^2 \left( q , \exponential[big]{p}( -\lambda \grad g(p)) \right) ,\tag{13}\] \(H(p^\ast; p, \lambda) \le H(p; p, \lambda)\). We state the nonconvex Riemannian Proximal Gradient method (NCRPG) in 1.
We point out that, in our experiments, the proximal subproblem is either solvable in closed-form or can be efficiently solved iteratively on the manifold. See 6 for more details.
4 Convergence↩︎
In this section, we analyze the convergence behavior of the nonconvex Riemannian proximal gradient (NCRPG) method from 1. We will assume the curvature of the manifold to be bounded, and no convexity assumptions are made for the functions \(g\) and \(h\).
The following assumptions are made in this section:
Assumption 1.
There exist two numbers \(\kmin, \kmax \in \bbR\) such that the sectional curvature satisfies \(\newtarget{def:k-min}{\kmin} \le K_p(X_p,Y_p) \le \newtarget{def:k-max}{\kmax}\) for all \(p \in \cM\) and all linearly independent unit vectors \(X_p , Y_p\in \tangentSpace{p}\).
For every starting point \(\sequence{p}{0} \in \cM\) the sublevel set \[\newtarget{def:level-set}{ \startlevelset } \coloneq \setDef[auto]{p \in \cM}{ f(p) \le f(\sequence{p}{0}) }\] is compact in \(\cM\).
\(g \colon \cM \to \bbR\) is continuously differentiable and \(\newtarget{def:gradg-lipschitz-constant}{\lipgrad}\)-smooth on \(\startlevelset\).
\(h \colon \cM \to \bbR\) is locally Lipschitz continuous and lower bounded.
In each step we find a stationary point \(\sequence{p}{k+1} \in \cM\) of \[H(q; \sequence{p}{k}, \sequence{\lambda}{k} ) \coloneq h(q) + \frac{1}{2 \sequence{\lambda}{k}} \dist^2 \Bigl(q, \exponential[big]{\sequence{p}{k}}( -\sequence{\lambda}{k} \grad g(\sequence{p}{k}) ) \Bigr) ,\] \(H(\sequence{p}{k+1}; \sequence{p}{k}, \sequence{\lambda}{k} ) \le H(\sequence{p}{k}; \sequence{p}{k}, \sequence{\lambda}{k} )\).
Note that since our sublevel set is compact and our objective function \(f\) is continuous, it is lower bounded by \(f_{opt} \coloneq \min_{q \in \startlevelset} f(q)\). Furthermore, since \(h\) and \(\grad g\) are continuous we obtain bounds \(\hlowerbound, \hupperbound, \gradgupperbound \in \bbR\) such that \[\newtarget{def:h-lower-bound}{ \hlowerbound } \le h(q) \le \newtarget{def:h-upper-bound}{ \hupperbound } \text{ for all } q \in \startlevelset, \text{ and } \norm{\grad g(q)} \le \newtarget{def:gradg-upper-bound}{ \gradgupperbound } \text{ for all } q \in \startlevelset.\]
Let \(p \in \cM\) and \(\lambda > 0\). We introduce the iteration map \(\imapOp[auto]{g}{h}{\lambda} \colon \cM \to \cM\) that to each point \(p \in \cM\) assigns a stationary point \(\newtarget{def:iteration-map}{\imap[auto]{g}{h}{\lambda}{p}} \in \cM\) of the function \(H(\cdot ;p, \lambda)\) defined in 13 . Note that, if \(h\) is geodesically convex, then \[\label{eq:iteration-mapping} \imap[auto]{g}{h}{\lambda}{p} = \prox{\lambda h}\Bigl( \exponential[big]{p}(-\lambda \grad g(p)) \Bigr) \in \cM .\tag{14}\] Furthermore, we introduce the gradient map \(\gmapOp[auto]{g}{h}{\lambda} \colon \cM \to \tangentBundle\) defined by \[\label{eq:gradient-mapping} \newtarget{def:gradient-map}{ \gmap[auto]{g}{h}{\lambda}{p} } \coloneq - \frac{1}{\lambda} \logarithm[big]{p}( \imap[auto]{g}{h}{\lambda}{p} ) \in \tangentSpace{p} .\tag{15}\] When \(g, h\) are clear from context, we denote \(\newtarget{def:iteration-map}{\Imap[auto]{\lambda}{p}} \equiv \imap[auto]{g}{h}{\lambda}{p}\), and \(\newtarget{def:gradient-map}{\Gmap{\lambda}{p}} \equiv \gmap[auto]{g}{h}{\lambda}{p}\), respectively.
4.1 Sufficient decrease↩︎
This first step towards a convergence result is to show a sufficient decrease of the cost function. For this, we need the following technical lemma, that provides a curvature-dependent upper bound on the distance between the points \(p\) and \(\Imap{\lambda}{p}\) and ensures the strict positivity of \(\zeta_{2, \kmax}\) within the sublevel set \(\startlevelset\). Notably, in case \(\kmax >0\), this lemma is crucial as it ensures that the gradient step \({\gradientstep{p} = \exponential{p}(-\lambda \grad g(p))}\) remains within a bounded distance from \(p\) and from \(\Imap{\lambda}{p}\) that is small enough to ensure that the logarithmic map, and hence also the gradient mapping \(\Gmap{\lambda}{p}\), is defined whenever it is needed.
Lemma 1. Let \(p \in \startlevelset\), \(\delta > 0\), \(\hlowerbound\), and \(\hupperbound\) be the numbers in [remark:nonconvex-assumptions:boundedness-of-h-and-grad-g], \[\label{eq:lambda-epsilon} \newtarget{def:delta-stepsize}{ \lambdadelta } \coloneq \frac{ \sqrt{ 4 \, (\hupperbound - \hlowerbound)^2 + \frac{{\pikmax}^2}{(2 + \delta)^2} \, \gradgupperbound^2 } - 2 \, (\hupperbound - \hlowerbound) }{2 \, \gradgupperbound^2} ,\qquad{(1)}\] and \[\label{eq:zeta-delta} \newtarget{def:delta-curvature}{ \zetadelta } \coloneq \zeta_{2, \kmax}\left( \deltaradius \right) = \begin{cases} 1 & \text{if } \kmax \le 0 , \\ \frac{\pi}{2 + \delta} \cot \left(\frac{\pi}{2 + \delta} \right) \quad & \text{if } \kmax > 0 , \end{cases}\qquad{(2)}\] with \(\newtarget{def:delta-radius}{\deltaradius} \coloneq \frac{\pikmax}{2 + \delta}\). Then, for all \(\lambda \in (0, \lambdadelta]\) we have \(\Imap{\lambda}{p}, \gradientstep{p} \in \cB(p, \deltaradius)\), where \(\cB(p, \deltaradius)\) is the uniquely geodesically convex open ball centered at \(p\) with radius \(\deltaradius\) and \(\newtarget{def:gradient-step}{\gradientstep{p} = \exponential{p}(-\lambda \grad g(p))}\). Furthermore, for all \(q \in \cB(p, \deltaradius)\), we have \[\label{eq:zetasecond-dist-bound} \zetasecond{\dist(p, q)} \ge \zetadelta > 0 .\qquad{(3)}\]
Proof. Let \(\delta > 0\), and observe that \(\Imap{\lambda}{p}, \gradientstep{p} \in \cB(p, \deltaradius)\) is achieved, , if both \(\dist(p, \gradientstep{p})\) and \(\dist(\Imap{\lambda}{p}, \gradientstep{p})\) are upper bounded by \(\frac{\deltaradius}{2}\). We have \(\dist(p, \gradientstep{p}) = \lambda \riemanniannorm{\grad g(p)}\) and [assumption:nonconvex-assumptions:decreaseprox] in 1 implies \[\begin{align} h(\Imap{\lambda}{p}) + \frac{1}{2 \lambda} \dist^2(\Imap{\lambda}{p}, \gradientstep{p}) \le h(p) + \frac{1}{2 \lambda} \dist^2(p, \gradientstep{p}) = h(p) + \frac{1}{2 \lambda} \lambda^2 \riemanniannorm{\grad g(p)}^2 , \end{align}\] so that \[\dist^2(\Imap{\lambda}{p}, \gradientstep{p}) \le 2 \lambda \left( h(p) - h(\Imap{\lambda}{p}) \right) + \lambda^2 \riemanniannorm{\grad g(p)}^2 .\] By [assumption:nonconvex-assumptions:h] and [remark:nonconvex-assumptions:boundedness-of-h-and-grad-g], we further obtain \[\dist^2(\Imap{\lambda}{p}, \gradientstep{p}) \le 2 \lambda \left( \hupperbound - \hlowerbound \right) + \lambda^2 \gradgupperbound^2 .\] Imposing now that \[2 \lambda \left( \hupperbound - \hlowerbound \right) + \lambda^2 \gradgupperbound^2 \le \left( \frac{\deltaradius}{2} \right)^2 ,\] we obtain \(\lambda \in \left(0, \lambdadelta \right]\), with \(\lambdadelta\) is as in ?? . Also note that this implies \(\lambda \le \frac{\deltaradius}{2 \gradgupperbound}\), which in turn implies \(\dist(p, \gradientstep{p}) \le \frac{\deltaradius}{2}\). Putting everything together we obtain \(\Imap{\lambda}{p}, \gradientstep{p} \in \cB(p, \deltaradius)\). The monotonicity of \(\zetasecond{s}\) on \((0, \deltaradius]\) then implies ?? , with \(\zetadelta\) as in ?? . ◻
With the above Lemma, we are able to show sufficient decrease of the cost function.
Lemma 2. Let \(g\) and \(h\) satisfy 1 and \(\delta > 0\). Then, for any \(p \in \startlevelset\) and \({\lambda \in \left(0, \min\bigl\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\bigr\}\right)}\) the following inequality holds \[\label{eq:sufficient-decrease} f(p) - f(\Imap[auto]{\lambda}{p}) \ge \frac{ \lambda\zetadelta - \lambda^2\lipgrad }{2} \riemanniannorm{\Gmap{\lambda}{p}}^2 .\qquad{(4)}\]
Proof. Let \(p \in \startlevelset\). Since \(\Imap[auto]{\lambda}{p}\) is a stationary point of the function \(H(\cdot ;p, \lambda)\) defined in 1, [assumption:nonconvex-assumptions:decreaseprox] and \(H(\Imap[auto]{\lambda}{p}; p, \lambda) \le H(p; p,\lambda)\), we have \[\begin{align} f(p) & \ge g(\Imap[auto]{\lambda}{p}) - \riemannian{\grad g(p)}{\logarithm{p}\Imap[auto]{\lambda}{p}} - \frac{\lipgrad}{2} \dist^2 \left( p , \Imap[auto]{\lambda}{p} \right) + h(\Imap[auto]{\lambda}{p}) \\ & \quad + \frac{1}{2 \lambda} \dist^2 \left( \exponential{p}(-\lambda \grad g(p)) , \Imap[auto]{\lambda}{p} \right) - \frac{1}{2 \lambda} \dist^2 \left( p , \exponential{p}(-\lambda \grad g(p)) \right) , \end{align}\] where we also used [assumption:nonconvex-assumptions:g]. Let now \({\gradientstep{p} \coloneq \exponential{p}(-\lambda \grad g(p))}\). We get \[\begin{gather} f(p) \ge f(\Imap{\lambda}{p}) + \frac{1}{2 \lambda} \left[ \dist^2(\Imap{\lambda}{p}, \gradientstep{p}) - \dist^2(p, \gradientstep{p}) - \lipgrad \lambda \dist^2(p, \Imap{\lambda}{p}) \right. \\ \left. + 2 \riemannian{\logarithm{p}\gradientstep{p}}{\logarithm{p}\Imap{\lambda}{p}} \right] . \end{gather}\] We are now in a position to apply the law of cosines on manifolds; see, , LezcanoCasado:2020:1?, or [16] which in this case reads \[\dist^2(\Imap{\lambda}{p}, \gradientstep{p}) \ge \zetasecond{\dist(m, p)} \dist^2(p, \Imap{\lambda}{p}) + \dist^2(p, \gradientstep{p}) - 2 \riemannian{\logarithm{p}\gradientstep{p}}{\logarithm{p}\Imap{\lambda}{p}} ,\] where \(m\) is a point on the minimal geodesic arc connecting \(\Imap{\lambda}{p}\) and \(\gradientstep{p}\). Since \(\lambda \leq \lambdadelta\), we obtain from 1 that \(\Imap{\lambda}{p}, \gradientstep{p} \in \cB(p, \deltaradius)\) and because this ball is uniquely geodesically convex, \(m \in \cB(p, \deltaradius)\), and thus \(\zetasecond{\dist(m, p)} \ge \zetadelta >0\). This yields \[\begin{align} f(p) & \ge f(\Imap{\lambda}{p}) + \frac{ \frac{\zetadelta}{\lambda} - \lipgrad }{2} \dist^2(p, \Imap{\lambda}{p}) . \end{align}\] Choosing a stepsize \(\lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\) gives the claim. ◻
4.2 Stepsizes↩︎
In 2, it was necessary to choose a stepsize \(\lambda\). We now present several strategies for this choice. Notice that, whenever \(\cM\) has nonpositive curvature \(\kmax \le 0\), it holds that \(\zetadelta = 1\). This means that to guarantee sufficient decrease, the stepsize \(\lambda\) can be then chosen as \(\lambda \in \left(0, \frac{1}{\lipgrad} \right)\) in this case.
We will consider a constant stepsize and a backtracking procedure. For the constant stepsize strategy, we choose \(\sequence{\lambda}{k} \coloneq \bar \lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\). The backtracking procedure is described in 2.
This means we look for the smallest \(i_k \in \bbN\) such that \[\label{eq:armijo-condition} f(\sequence{p}{k}) - f( \Imap{s \eta^{i_k}}{\sequence{p}{k}} ) \ge \beta s \eta^{i_k} \riemanniannorm[auto]{ \Gmap{s\eta^{i_k}}{\sequence{p}{k}} }^2\tag{16}\] holds. Notice that this procedure is finite since for \[\label{eq:backtracking-lambda-upper-bound} \sequence{\lambda}{k} < \min \left\{ \lambdadelta , \frac{\zetadelta - 2 \beta}{\lipgrad} \right\} ,\tag{17}\] we have \(\frac{\zetadelta- \sequence{\lambda}{k} \lipgrad}{2} \ge \beta ,\) and thus by 2 \[\begin{align} f(\sequence{p}{k}) - f( \Imap{\sequence{\lambda}{k}}{\sequence{p}{k}} ) \ge \sequence \lambda k \frac{ \zetadelta - \sequence{\lambda}{k}\lipgrad}{2} \riemanniannorm[auto]{\Gmap{\sequence{\lambda}{k}}{\sequence{p}{k}}}^2 \ge \beta \sequence \lambda k \riemanniannorm[auto]{\Gmap{\sequence{\lambda}{k}}{\sequence{p}{k}}}^2 , \end{align}\] so the procedure must end as soon as \(\sequence{\lambda}{k}\) satisfies 17 .
We can use this observation to obtain a lower bound for \(\sequence \lambda k\). Let \(\sequence \lambda k\) be the stepsize obtained by the backtracking procedure. Then either \(\sequence \lambda k = s\), or \(\frac{\sequence \lambda k}{\eta}\) has not fulfilled [eq:backtracking-sufficient-decrease], meaning \[\label{eq:backtracking-lambda-lower-bound} \sequence{\lambda}{k} \ge \min \left\{ s, \eta\lambdadelta, \eta\frac{\zetadelta - 2 \beta}{\lipgrad} \right\} .\tag{18}\] Taking either of the stepsizes above allows us to obtain a sufficient decrease property.
Corollary 1. Suppose that 1 hold. Let \(\sequence(){p}{k}\) be the sequence generated by 1 for some starting point \(\sequence p 0\) and either a constant stepsize \(\sequence{\lambda}{k} \coloneq \bar \lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\) or the backtracking procedure with constants \((s, \beta, \eta)\) as described in 4.2. Then for any \(k \ge 0\), \[\label{eq:sufficient-decrease-with-constant} f(\sequence{p}{k}) - f(\sequence{p}{k+1}) \ge M \riemanniannorm[auto]{ \Gmap{\sequence \lambda k}{\sequence p k} }^2,\qquad{(5)}\] where \[\label{eq:sufficient-decrease-M-d} \begin{align} M & \coloneq \begin{cases} \frac{ \bar\lambda\zetadelta - \bar\lambda^2\lipgrad }{2} & \text{constant stepsize} , \\ \beta \min \left\{ s, \eta\lambdadelta, \eta\frac{\zetadelta - 2 \beta}{\lipgrad} \right\} & \text{backtracking} . \end{cases} \end{align}\qquad{(6)}\]
4.3 Stationary points are zeros of the gradient mapping↩︎
The following result gives a characterization of the stationary points of the optimization problem 1 in terms of the gradient mapping.
Theorem 1. Suppose that 1 hold. Then \(\gmap[auto]{g}{h}{\lambda}{p^\ast} = 0\) if and only if \(p^\ast \in \startlevelset\) is a stationary point of 1 .
Proof. Observe that \(\gmap{g}{h}{\lambda}{p^\ast} = 0 \in \tangentSpace{p^\ast}\) if and only if \[p^\ast = \imap{g}{h}{\lambda}{p^\ast} .\] This is equivalent to \(H(p^\ast; p^\ast, \lambda) = H(\imap{g}{h}{\lambda}{p^\ast}; p^\ast, \lambda)\) with \[0 \in \partial H(p^\ast; p^\ast, \lambda) = \partial h(p^\ast) - \frac{1}{\lambda} \logarithm[big]{ p^\ast } ( \exponential{p^\ast}(-\lambda \grad g(p^\ast)) ) ,\] which is equivalent to \[\frac{1}{\lambda} \logarithm[big]{p^\ast}(\exponential{p^\ast}(-\lambda \grad g(p^\ast))) \in \partial h(p^\ast) ,\] , if and only if \[- \grad g(p^\ast) \in \partial h(p^\ast) ,\] which is the condition for stationarity given in 8 . ◻
Lemma 3. Suppose that 1 hold. Let \(\sequence(){p}{k}\) be the sequence generated by 1 for some starting point \(\sequence p 0\) and either a constant stepsize \(\sequence{\lambda}{k} \coloneq \bar \lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\) or the backtracking procedure with constants \((s, \beta, \eta)\). Then,
the sequence \(\left(f(\sequence{p}{k})\right)\) is nonincreasing, and \(f(\sequence{p}{k+1}) < f(\sequence{p}{k})\) if and only if \(\sequence{p}{k}\) is not a stationary point of 1 ; in particular \(\sequence(){p}{k} \subseteq \startlevelset\).
\(\Gmap{\sequence \lambda k}{\sequence p k} \to 0\) as \(k \to \infty\);
\(\min_{n = 0, \ldots, k} \riemanniannorm{\Gmap{\sequence \lambda n}{\sequence{p}{n}}} \le \sqrt{\frac{f(\sequence p 0)- f_{opt}}{M(k+1)}}\), where \(M\) is given in ?? and \(f_{opt}\) is the optimal solution of \(f\) in the sublevel set \(\startlevelset\) as described in [remark:nonconvex-assumptions:boundedness-of-h-and-grad-g].
Proof.
If \(\sequence{p}{k}\) is not a stationary point of 1 , then \(\Gmap{\sequence \lambda k}{\sequence{p}{k}} \neq 0\) by 1. Thus, by ?? we have \(f(\sequence{p}{k+1}) < f(\sequence{p}{k})\). Vice versa, if \(\sequence{p}{k}\) is a stationary point of 1 , then \(\Gmap{\sequence \lambda k}{\sequence{p}{k}} = 0\), which implies \(\sequence{p}{k+1} = \sequence{p}{k}\), hence \(f(\sequence{p}{k}) = f(\sequence{p}{k+1})\).
The sequence \(\left(f(\sequence{p}{k})\right)\) converges because it is nonincreasing and bounded from below. This implies \(f(\sequence{p}{k}) - f(\sequence{p}{k+1}) \to 0\) as \(k \to \infty\), which combined with 2 yields \(\Gmap{\sequence \lambda k}{\sequence{p}{k}} \to 0\).
Summing ?? over \(n=0, \ldots, k\) yields \[\begin{align} f(\sequence{p}{0}) - f_{opt} & \ge f(\sequence{p}{0}) - f(\sequence{p}{k+1}) \ge M \sum_{n=0}^k \riemanniannorm[auto]{ \Gmap{\sequence \lambda n}{\sequence p n} }^2 \\ & \ge M(k+1) \min_{n = 0, \ldots, k} \riemanniannorm[auto]{ \Gmap{\sequence \lambda n}{\sequence{p}{n}} }^2 . \qedhere \end{align}\]
◻
4.4 Sublinear convergence to \(\varepsilon\)-stationary points↩︎
The next lemma shows that the subdifferential of the cost function at each iteration is bounded by the norm of gradient mapping.
Lemma 4. Let \(g\) and \(h\) satisfy 1. Let \(\sequence(){p}{k}\) be the sequence generated by 1 for some starting point \(\sequence p 0\) and either a constant stepsize \(\sequence{\lambda}{k} \coloneq \bar \lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\) or the backtracking procedure with constants \((s, \beta, \eta)\). Then, for every \(k \in \bbN\), there exists a tangent vector \(X_{\sequence{p}{k}} \in \partial f\left(\sequence{p}{k}\right)\) such that \[\riemanniannorm{X_{\sequence{p}{k}}} \le ( \lipgrad + \sigmakmax ) \sequence \lambda k \, \riemanniannorm[auto]{ \Gmap{\sequence \lambda k}{\sequence{p}{k}} } ,\] where \[\label{eq:stationary-sigma-definition} \newtarget{def:curvature-upper-bound}{ \sigmakmax } \coloneq \begin{cases} \sigmacurvature{ 2 \deltaradius } & \text{if } \kmax > 0 , \\ \sigmacurvature{ t \, \gradgupperbound + \diam(\startlevelset) } & \text{if } \kmax \le 0 , \end{cases}\qquad{(7)}\] with \(\deltaradius\) as defined in 1. Here, \(t\) is chosen as \(\lambda\) if a constant stepsize is used, and as the initial guess \(s\) in the backtracked case.
Proof. Let \(\gradientstep{\sequence{p}{k}} \coloneq \exponential{\sequence{p}{k}}({-\sequence{\lambda}{k} \grad g(\sequence{p}{k})})\), and \(\sequence{Y}{k+1} \coloneq \grad g(\sequence{p}{k+1})\), \(\sequence{Z}{k+1} \coloneq \logarithm{\sequence{p}{k+1}}(\gradientstep{\sequence{p}{k}})\). Since \(\sequence{p}{k+1}\) is stationary for \(h(p) + \frac{1}{2 \sequence{\lambda}{k}} \dist^2( p , \gradientstep{\sequence{p}{k}} ),\) we have \[\begin{align} 0_{\sequence{p}{k+1}} \in \partial h(\sequence{p}{k+1}) - \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} = \partial f( \sequence{p}{k+1} ) - \sequence{Y}{k+1} - \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} , \end{align}\] where we used the relation \(\partial f(p) = \partial h(p) + \grad g(p)\); see 7 . Rearranging the terms we get \[\sequence{Y}{k+1} + \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} \in \partial f(\sequence{p}{k+1}) .\] We estimate the norm as follows: \[\label{eq:eps95stat95first95part} \begin{align} & \riemanniannorm[auto]{ \sequence{Y}{k+1} + \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} } \\ & \qquad \le \riemanniannorm[auto]{ \sequence{Y}{k+1} - \parallelTransport{\sequence{p}{k}}{\sequence{p}{k+1}} \sequence{Y}{k} } + \riemanniannorm[auto]{ \parallelTransport{\sequence{p}{k}}{\sequence{p}{k+1}} \sequence{Y}{k} + \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} } \\ & \qquad \le \lipgrad \dist(\sequence{p}{k}, \sequence{p}{k+1}) + \frac{1}{\sequence{\lambda}{k}} \riemanniannorm[auto]{ \parallelTransport{\sequence{p}{k}}{\sequence{p}{k+1}} \sequence{W}{k} - \sequence{Z}{k+1} } , \end{align}\tag{19}\] where the second inequality is due to the \(\lipgrad\)-smoothness of the function \(g\), and we also used \({\sequence{W}{k} \coloneq \logarithm{\sequence{p}{k}}(\gradientstep{\sequence{p}{k}}) = - \sequence{\lambda}{k}\sequence{Y}{k}}\). It holds that \[\riemanniannorm[auto]{ \parallelTransport{\sequence{p}{k}}{\sequence{p}{k+1}} \sequence{W}{k} - \sequence{Z}{k+1} } \le \sigmacurvature{\dist(m, \gradientstep{\sequence{p}{k}})} \dist(\sequence{p}{k+1}, \sequence{p}{k}) ,\] where \(m\) is a point on the geodesic arc that connects \(\sequence{p}{k}\) and \(\sequence{p}{k+1}\). This is an implication of AlimisisOrvietoBecigneulLucchi:2021:1? with the notation adapted to our setting. Using the triangle inequality together with the fact that \(m\) lies on such geodesic, we have \(\dist(m, \gradientstep{\sequence{p}{k}}) \le \dist(\sequence{p}{k}, \gradientstep{\sequence{p}{k}}) + \dist(\sequence{p}{k}, \sequence{p}{k+1})\). Thus by 1, \(\dist(m, \gradientstep{\sequence{p}{k}})\) is bounded by \(2 \deltaradius < \pikmax\) for all \(k \in \bbN\), which in turn implies \(\sigmacurvature{\dist(m, \gradientstep{\sequence{p}{k}})} \le \sigmacurvature{2 \deltaradius}\) for all \(k \in \bbN\) if \(\kmax > 0\). If instead \(\kmax \le 0\), \(\pikmax\) and \(\deltaradius\) are unbounded, and we thus need to argue with the boundedness of the sublevel sets. We get \[\dist(m , \gradientstep{\sequence{p}{k}}) \le \dist(\sequence{p}{k}, \gradientstep{\sequence{p}{k}}) + \dist(\sequence{p}{k}, \sequence{p}{k+1}) \le t \gradgupperbound + \diam(\startlevelset) ,\] where \(t\) is chosen as \(\lambda\) if a constant stepsize is used, or as the initial guess \(s\) in the backtracked case, and \(\gradgupperbound\) is the bound on \(\grad g\) on \(\startlevelset\); see [remark:nonconvex-assumptions:boundedness-of-h-and-grad-g]. This implies \({\sigmacurvature{\dist(m, \gradientstep{\sequence{p}{k}})} \le \sigmacurvature{t \gradgupperbound + \diam(\startlevelset)}}\) for all \(k \in \bbN\), which is finite by [assumption:nonconvex-assumptions:compact95levelset]. Inserting ?? into 19 yields \[\label{eq:eps-stat-second-part} \riemanniannorm[auto]{ \sequence{Y}{k+1} + \frac{1}{\sequence{\lambda}{k}} \sequence{Z}{k+1} } \le ( \lipgrad + \sigmakmax ) \dist(\sequence{p}{k+1}, \sequence{p}{k}) ,\tag{20}\] The claim follows after substituting \(\dist(\sequence{p}{k+1}, \sequence{p}{k}) = \sequence \lambda k \riemanniannorm[auto]{\Gmap{\sequence \lambda k}{\sequence{p}{k}}}\). ◻
The following is a complexity result for the convergence of the iterates to stationary points.
Theorem 2. Let \(g\) and \(h\) satisfy 1 and \(\varepsilon > 0\). Then 1 finds a point \(\sequence{p}{k}\) such that \(\riemanniannorm[auto]{X_{\sequence{p}{k}}} \le \varepsilon\), with \(X_{\sequence{p}{k}} \in \partial f(\sequence{p}{k})\), in at most \[\label{eq:complexity-rate} k = \frac{ ( \lipgrad + \sigmakmax )^2 (f(\sequence p 0) - f_{opt}) t^2 }{ \varepsilon^2 M } - 1\qquad{(8)}\] iterations, where \(t\) is chosen as \(\lambda\) in the case of a constant stepsize and the initial guess \(s\) for the backtracking procedure.
Proof. Choose \(n \ge \frac{ ( \lipgrad + \sigmakmax )^2 (f(\sequence p 0) - f_{opt}) t^2 }{ \varepsilon^2 M } - 1\). Notice that our choice of \(t\) ensures \(\sequence \lambda k \le t\) holds for all \(k\). Then, by 3, there exists a \(k \le n\) such that \[\riemanniannorm{\Gmap{\sequence \lambda k}{\sequence{p}{k}}} \le \sqrt{ \frac{f(\sequence p 0)- f_{opt}}{M(n + 1)} } \le \frac{ \varepsilon }{ ( \lipgrad + \sigmakmax ) t } .\] Combined with 4, we obtain an \(X_{\sequence{p}{k}} \in \partial f\big(\sequence{p}{k}\big)\) such that \[\riemanniannorm[auto]{ X_{\sequence{p}{k}} } \le ( \lipgrad + \sigmakmax ) \sequence \lambda k \riemanniannorm{\Gmap{\sequence \lambda k}{\sequence{p}{k}}} \le \frac{ \sequence \lambda k }{ t } \varepsilon \le \varepsilon .\qedhere\] ◻
4.5 Global convergence to stationary points↩︎
We can now give a global convergence result for 1, inspired by [10].
Theorem 3. Let \(g\) and \(h\) satisfy 1. Let \(\sequence(){p}{k}\) be the sequence generated by 1 for some starting point \(\sequence{p}{0}\) and either a constant stepsize \(\sequence{\lambda}{k} \coloneq \bar \lambda \in \left(0, \min\left\{\lambdadelta, \frac{\zetadelta}{\lipgrad}\right\}\right)\) or the backtracking procedure with constants \((s, \beta, \eta)\). Then, the sequence \(\sequence(){p}{n}\) has at least one accumulation point, and any accumulation point is a stationary point of \(f\).
Proof. By the compactness of the sublevel set \(\startlevelset\) of \(f\) and 1, the sequence \(\sequence(){p}{k}\) has at least one accumulation point. Let now \(\sequence(){p}{n_j}\) be a subsequence converging to a point \(p^\ast \in \startlevelset\). Then, by 4 we find a sequence of tangent vectors \(X_{\sequence{p}{n_j}} \in \partial f\big(\sequence{p}{n_j}\big)\) such that \[\lim_{j \to \infty} \riemanniannorm[auto]{X_{\sequence{p}{n_j}}} \le \lim_{j \to \infty} ( \lipgrad + \sigmakmax ) \sequence \lambda {n_j} \, \riemanniannorm[auto]{ \Gmap{\sequence \lambda {n_j}}{\sequence{p}{n_j}} } = 0,\] where the last step follows from 3. 1 implies the continuity of \(f\). Thus, \(f\left(\sequence{p}{n_j}\right) \to f(p^\ast)\), and since we found a (sub)sequence \(X_{\sequence{p}{n_j}} \in \partial f\left(\sequence{p}{n_j}\right)\) such that \(X_{\sequence{p}{n_j}} \to 0_{p^\ast}\), we can conclude that \(0_{p^\ast} \in \partial f(p^\ast)\) by [17], implying the stationarity of \(p^\ast\). ◻
5 A Riemannian proximal gradient method using retractions↩︎
Often, we do not have access to a closed form of the exponential map \(\exponential{p}{X}\) or its inverse, the logarithmic map. In other cases, computing these maps is simply too costly. For first order convergence results, it is often enough to consider first order approximations of the exponential map, so called retractions, and their inverses. In [9], the authors have shown convergence of their algorithm also for such retractions. This alters the proximal problem on the tangent space but it does not directly affect the algorithm itself. In this section, we will discuss the use of retractions in the NCRPG method and conclude that convergence theory remains largely an open problem.
A retraction on a Riemannian manifold \(\cM\) is a first-order approximation of the exponential map \(\retractionSymbol \colon \tangentBundle \to \cM , (p, X_p) \mapsto \retract{p}(X_p) ,\) such that each curve \(c(t) = \retract{p}(tX_p)\) satisfies \(c(0) = p\) and \({c'(0) = X_p}\). Similarly, the inverse retraction \(\inverseRetract{p}(q)\) is a first order approximation of the logarithmic map \(\logarithm{p} q\). Note that analogously to the logarithmic map, both the retraction and the inverse retraction might only be locally defined around \(p\). We can use these approximations to consider both an approximation of the Riemannian distance \(\dist\) and the proximal map \(\prox{\lambda f}\). We denote \[\distR(p,q) \coloneq \riemanniannorm{\inverseRetract{p}(q)}[p] \qquad \text{and} \qquad \proxR{\lambda f}(p) \coloneq \argmin_{q\in\cM} \frac{1}{2\lambda}\distR(q,p)^2 + f(q) ,\] for a given retraction \(\retractionSymbol\), respectively. We note that the retraction distance \(\distR\) is not symmetric but both approximations \(\distR(p,q) \approx \dist(p,q) \approx \distR(q,p)\) are equally valid. We have chosen the latter in the retraction-based proximal map, as it leads to nicer expressions.
We can then replace the minimization of 13 by a retraction-based proximal map. \[\begin{align} \label{eq:retraction-prox} \proxR[big]{\lambda h}[\retractionSymbol]( \mathrm{R}_{p}(-\lambda\grad g(p)) ) &= \argmin_{q\in\cM} \frac{1}{2\lambda} \Bigl( \distR[big]( q,\mathrm{R}_{p}(-\lambda\grad g(p)) )[\retractionSymbol]\Bigr)^2 + h(q) \\ &= \argmin_{q\in\cM} \frac{1}{2\lambda} \riemanniannorm[auto]{ \inverseRetract{q}\Bigl(\mathrm{R}_{p}\bigl(-\lambda\grad g(p)\bigr)\Bigr) }^2 + h(q) . \end{align}\tag{21}\] Or, as above, by finding a stationary point \(p^\ast\) of the objective function \[\label{equation:retraction-objective} H^\retractionSymbol(q; p, \lambda ) \coloneq h(q) + \frac{1}{2 \lambda} \left( \distR[Big]( q, \retract{p}\left( -\lambda \grad g(p)\right) )[\retractionSymbol] \right)^2,\tag{22}\] which fulfills \(H^\retractionSymbol(p^\ast; p, \lambda) \leq H^\retractionSymbol(p; p, \lambda).\) This is again different from [9] in that the problem is posed on the manifold as opposed to a tangent space. We then summarize the retraction-based algorithm in 3.
The main problem with obtaining a convergence theory of this method - at least using the considerations above - is that 2 only holds for the exact Riemannian proximal gradient method, where steps are taken along geodesics. Generalizing this result to approximations thereof, where the steps are taken using retractions, will add one of the following difficulties:
Either it would be necessary to generalize the law of cosines for (some special) retractions. So far, this seems out of reach. Alternatively, it could be possible to rewrite the update steps in terms of geodesic distances and then apply the standard law of cosines. One way to do this requires 23 below. We keep the modified proximal objective function \(H^\retractionSymbol\) but the \(\lipgrad\) smoothness in [assumption:nonconvex-assumptions:g] is still stated in terms of logarithm and geodesic distance. That is, \[g(p) \ge g(q) - \riemannian{\grad g(p)}{\logarithm{p}q} - \frac{\lipgrad}{2} \dist^2 \left( p , q \right).\] Combining these, we now get \[\begin{align} f(p) & \ge f(\Imap[auto]{\lambda}{p}) + \frac{1}{\lambda} \riemannian{ \logarithm[big]{p}({\exponential{p}(- \lambda \grad g(p))}) }{ \logarithm{p} \Imap[auto]{\lambda}{p} } - \frac{\lipgrad}{2} \dist^2 \left( p , \Imap[auto]{\lambda}{p} \right) \\ & \quad + \frac{1}{2 \lambda} \distR[auto]( \Imap[auto]{\lambda}{p} , \retract{p}(-\lambda \grad g(p)) )[\retractionSymbol]^2 - \frac{1}{2 \lambda} \distR[auto]\left( p , \retract{p}(-\lambda \grad g(p))\right)^2 \end{align}\] Define now \(q \coloneq \Imap[auto]{\lambda}{p}\), and \(\gradientstep{p} \coloneq \exponential{p}(-\lambda \grad g(p))\), and get \[\begin{align} f(p) & \ge f(q) + \frac{1}{\lambda} \riemannian{\logarithm{p}(\gradientstep{p})} {\logarithm{p}(q)} - \frac{\lipgrad}{2} \dist^2 \left( p , q \right) \\ & \quad + \frac{1}{2 \lambda} \riemanniannorm[auto]{ \inverseRetract{q}{\retract{p}(-\lambda \grad g(p))} }^2 - \frac{1}{2 \lambda} \dist^2\left( p , \gradientstep{p}\right) \end{align}\]
At this point, if there is a constant \(c \in \bbR\) such that for all \(p, q \in \cM\) and \(v \in \tangentSpace{p}\), \[\label{eq:retraction-property} \riemanniannorm[auto]{ \inverseRetract{q}{\retract{p}(v)} }^2 \ge \riemanniannorm[auto]{ \logarithm{q}{\exponential{p}(v)} }^2 - c \dist^2 \left(p,q \right)\tag{23}\] holds, we obtain \[\begin{align} f(p) & \ge f(q) + \frac{1}{\lambda} \riemannian{\logarithm{p}{\gradientstep{p}}} {\logarithm{p}{q}} - \frac{\tilde{L}_g}{2} \dist^2 \left( p , q \right) + \frac{1}{2 \lambda} \dist^2 \left( \gradientstep{p} , q \right) - \frac{1}{2 \lambda} \dist^2\left( p , \gradientstep{p}\right), \end{align}\] where \(\tilde{L}_g = \lipgrad + \frac{c}{\lambda}\). The rest of the proof would now be the same as before.
The requirement in 23 can be seen as a first-order approximation property of the exponential (and logarithmic map) in the first, manifold argument, instead of the second, tangent space argument. Whether we can find retractions that satisfy this requirement, or if that is even always possible, is out of the scope of this work and we will leave it, as well as the whole convergence theory of 3, as an open problem.
6 Numerical Examples↩︎
In this section, we present numerical experiments to evaluate the performance of the Nonconvex Riemannian Proximal Gradient (NCRPG) method described in 1. We consider three applications on manifolds that are non-Hadamard (, they do not exhibit nonpositive curvature everywhere). The NCRPG method has been implemented in within . The implementation is available with the level of generality provided by retractions as well, as discussed in 5. The numerical experiments were performed on a MacBook Pro M1, 16 GB RAM running macOS 26.5. The experiments are available as notebooks in the package [18].
6.1 Sparse PCA↩︎
The first example is Sparse PCA. Here, we consider the following optimization problem: \[\min_{X \in \text{OB}(n,r)} f(X) = \frac{1}{2} \| X^\top A^\top A X - D^2 \|^2 + \mu \norm{X}_1 ,\] where \(A \in \bbR^{m\times n}\) is a given data matrix, \(D\) contains the \(r\) leading singular values of \(A\) on the diagonal, \(\mu \ge 0\) is a weighting parameter, and \(\text{OB}(n, r) = \bigotimes_{i = 1}^r \mathbb{S}^{n-1}\) is the Oblique manifold of matrices with normalized columns. Hence, the goal to diagonalize \(A^\top A\) as much as possible using normalized components that are forced to be sparse by the \(\ell_1\)-penalty. This variation of the problem was introduced in [19]. We note that there are alternative formulations of the Sparse PCA problem on the Stiefel manifold that we do not consider here.
6.1.0.1 Soft Thresholding on the Sphere
The main advantage of our algorithm in this case is that the prox-problem can be efficiently solved using a fixed-point iteration, as we will show in the following. Since the Oblique manifold is a product of spheres, the prox-problem decomposes and it is sufficient to solve it on the sphere. Thus, we aim to solve \[\label{eq:soft-thresholding-on-the-sphere:problem} \argmin_{y \in \bbS^{n-1}} \phi(y) \coloneq \frac{1}{2} \dist(x, y)^2 + \xi \norm{y}_1 ,\tag{24}\] where \(x \in \bbS^{n-1}\), \(\xi \ge 0\), and \(\dist(x, y) = \arccos(\langle x, y\rangle)\) is the Riemannian distance.
We will now show that the minimum \(\bar y = \argmin_{y \in \bbS^{n-1}} \phi(y)\) is attained at \(\bar y = p_x(\bar t)\) and \(\bar t\) is a fixed point \(\bar t = \sigma(p_x(\bar t))\) where \[\label{eq:def-sigma} \sigma(y) \coloneq \xi \frac{ \sqrt {1 - \langle x ,{y}\rangle ^2} }{ \arccos \left( \langle x, y\rangle \right) } , \qquad \text{and} \qquad p_x(t) \coloneq \frac{\prox{ t \norm{\cdot}_1}^{\bbR^{n+1}}(x)}{\norm{\prox{ t \norm{\cdot}_1}^{\bbR^{n+1}}(x)}_2}\tag{25}\] is the normalized Euclidean proximal operator. Then, \(\bar t\), and with that also \(\bar y\) can be computed in linear time using Banach’s fixed point theorem.
Notice that, since \(\phi\) is continuous and \(\bbS^{n-1}\) is compact, such a \(\bar y\) exists, and \(0_{\bar y} \in \partial \phi(\bar y)\). Hence, \[\label{eq:soft-thresholding-on-the-sphere:zero-subgrad} \begin{align} 0_{\bar y} & = - \logarithm{\bar y} x + \xi \proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(z) \\ & = -\frac{ \arccos \left(\langle x, {\bar y}\rangle \right) }{ \Vert\proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(x)\Vert_2 } \proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(x) + \xi \proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(z) , \end{align}\tag{26}\] for some \(z \in \partial \norm{\bar y}_1\). With the projection onto the tangent space at \(\bar y\) given by \[\proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(x) = x - \langle x, \bar y \rangle \bar y ,\] we see that 26 simplifies to \(0_{\bar y} = \proj{\tangentSpace{\bar y}[\bbS^{n-1}]}(x-\sigma(y) z)\), which is equivalent to \[\label{eq:soft-thresholding-on-the-sphere:zero-subgrad-simplified} x - \sigma(\bar y) z = \langle{x - \sigma(\bar y) z},{\bar y}\rangle \bar y .\tag{27}\] The solution of the proximal operator of the function \(t \norm{\cdot}_1\) in Euclidean space \(\bbR^{n+1}\) for a point \(x\in \bbR^{n+1}\) is given elementwise by \[\left( \prox{t \norm{\cdot}_1}^{\bbR^{n+1}}(x) \right)_i = \sgn(x_i) \max \left\{ 0 , \abs[auto]{x_i} - t \right\} , \qquad i = 1 , \dots , n+1 ,\] see, , Beck:2017:1?.
The following theorem shows that, given \(x \in \bbS^{n-1}\), the solution of 24 can be obtained by applying the normalized Euclidean proximal operator of \(t \norm{\cdot}_1\) to \(x\), for a particular choice of \(t\).
Theorem 4. Let \(x\in \bbS^{n-1}\) and let \(\bar y\) be a solution of the optimization problem 24 with \(0 < \xi < \norm{x}_\infty\). Then \(\bar y = p_x(\bar t)\) with \(\bar t\) being a fixed point of \(\sigma(p_x(t))\).
Proof. Let \(\bar y\) be a solution of 24 and \(i \in \{1,\ldots,n\}\) be a fixed index. We first show that, if \(x_i > 0\) (\(x_i < 0\)), then \(\bar y_i \ge 0\) (\(\bar y_i \le 0\)): Assume \(x_i > 0\) and \(\bar y_i < 0\). Define \(\hat{y}\) as \(\hat{y}_i \coloneq - \bar y_i\) and \(\hat{y}_j \coloneq \bar y_j\) for \(j \ne i\). Then \(\hat{y} \in \bbS^{n-1}\) and \(\norm {\hat{y}}_1 = \norm {\bar y}_1\). But \(\langle x, {\hat{y}}\rangle = \langle{x},{\bar y}\rangle - 2 x_i \bar y_i > \langle {x}, {\bar y}\rangle\), hence \(\arccos \left(\langle{x},{\hat{y}}\rangle \right) < \arccos \left(\langle{x},{\bar y} \rangle\right)\), implying \(\phi(\hat{y}) < \phi({\bar y})\), contradicting the optimality of \(\bar y\). The same argument holds for \(x_i < 0\) and \(\bar y_i > 0\).
Let \(z \in \partial \norm {\bar y}_1\) be such that 27 holds. Notice that, since \(\langle{\bar y},{x}\rangle \in [0, 1]\), we have \(\sigma(\bar y) \in \left[\frac{2\xi}{\pi}, \xi\right]\). W.l.o.g. we can assume that \(x_1 = \norm{x}_\infty\). Since \(\xi < x_1\) and \(\sigma(\bar y) \le \xi\), we have \(x_1 - \sigma(\bar y) z_1 > 0\). Hence, by 27 , \(\bar y_1 \ne 0\). Combined with \(\bar y_1 \ge 0\) we have \[\langle{ x - \sigma(\bar y) z },{\bar y} \rangle = \frac{ x_1 - \sigma(\bar y) z_1 }{\bar y_1} > 0 .\] Let now \(i \in \{2, \ldots, n\}\). We first consider the case \(\abs[auto]{x_i} \le \sigma(\bar y)\). Assume that \(\bar y_i \ne 0\), then \({\sgn(x_i - \sigma(\bar y) z_i) \in \{-\sgn(\bar y_i), 0\}}\). But then by 27 \[\sgn(\bar y_i) = \sgn( \langle{ x - \sigma(\bar y) z },{\bar y} \rangle ) \sgn(\bar y_i) = \sgn( x_i - \sigma(\bar y) z_i ) \in \{ -\sgn(\bar y_i) , 0 \} ,\] which is a contradiction. Thus \(\bar y_i = 0\) and \(z_i = \frac{x_i}{\sigma(\bar y)}\) must hold.
If \(\abs[auto]{x_i} > \sigma(\bar y)\) we have \(\sgn(x_i - \sigma(\bar y) z_i) = \sgn(x_i)\), implying \(\bar y_i \ne 0\) by 27 . In summary, we obtain \[x_i - \sigma(\bar y) z_i = \begin{cases} 0 & \text{ if } \abs[auto]{x_i} \le \sigma(\bar y) , \\ x_i - \sigma(\bar y) \sgn(x_i) & \text{ otherwise} , \end{cases}\] , \(x - \sigma(\bar y) z = \prox{\sigma(\bar y) \norm{\cdot}_1}^{\bbR^{n+1}}(x)\). Now, since \(\norm{\bar y}_2 = 1\) and \(\langle{x - \sigma(\bar y) z},{\bar y}\rangle > 0\), we have \[\bar y = \frac{ x - \sigma(\bar y) z }{ \norm{ x - \sigma(\bar y) z }_2 } = p_x(\sigma(\bar y)) .\] In particular this implies that \(\bar t = \sigma (\bar y)\) is a fixed point of \(\sigma(p_x(t))\). ◻
With the above theorem, we can conclude that we are looking for a fixed point of \(\sigma(p_x(t))\). We will now show that this fixed point can be computed using Banach’s fixed-point theorem.
Lemma 5. Let \(x \in \bbS^{n-1}\) and \[\label{eq:mu-upper-bound} 0 < \xi < \min\left\{ \frac{\pi}{8 \sqrt{\norm x_0}} \left( \sqrt{ 16 \sqrt{\norm x_0} \norm{x}_\infty + \pi^2 } - \pi \right) , \norm{x}_\infty \right\} .\qquad{(9)}\] Then \(\sigma(p_x(t)) \colon (0, \xi] \to (0,\xi]\) has a unique fixed point \(\bar t\), and for every \(t_0 \in (0, \xi]\) the sequence \({t_{k+1} \coloneq \sigma(p_x(t_k))}\) converges linearly to \(\bar t\).
Proof. Notice that, by our choice of \(\xi\), \(p_x(t)\) is well-defined for all \(t \in (0, \xi]\) and \(\sigma(p_x(t)) \in (0, \xi]\). We now show that \(p_x(t)\) is a contraction mapping. Define \(\psi(s) \coloneq \xi \frac{\sqrt {1- s^2}}{\arccos(s)}\) for \(s \in [0, 1]\). We then have \[\sup_{s \in [0, 1)} \psi'(s) = \sup_{s \in [0, 1)} \xi \frac{ 1 - \frac{s \arccos(s)}{\sqrt{1 - s^2}} }{\arccos(s)^2} = \frac{4 \xi}{\pi^2} .\] Thus, for any \(t_1, t_2 \in (0, \xi]\), we obtain \[\begin{align} \abs[auto]{ \sigma(p_x(t_1)) - \sigma(p_x(t_2)) } & = \abs[auto]{ \psi(\langle{p_x(t_1)},{x}\rangle) - \psi(\langle{p_x(t_2)},{x}\rangle) } \\ & \le \frac{4 \xi}{\pi^2} \abs[auto]{ \langle{p_x(t_1)},{x}\rangle - \langle{p_x(t_2)},{x}\rangle } . \end{align}\] Since \(p_x(t), x \in \bbS^{n-1}\), we have \[\langle{p_x(t)},{x}\rangle = 1 - \frac{1}{2} \norm{p_x(t) - x}_2^2 .\] Hence, for any \(t_1, t_2 \in (0, \xi]\), we rewrite \[\begin{align} \abs[auto]{ \langle{p_x(t_1)},{x}\rangle - \langle{p_x(t_2)},{x}\rangle } & = \abs[auto]{ \frac{1}{2} \norm{p_x(t_1) - x}_2^2 - \frac{1}{2} \norm{p_x(t_2) - x}_2^2 } . \end{align}\] W.l.o.g., \(\frac{1}{2} \norm{p_x(t_1) - x}_2^2 \ge \frac{1}{2} \norm{p_x(t_2) - x}_2^2\). We now use that \(p_x(t_1)\) is a solution of \[\argmin_{y\in \bbS^{n-1}} \quad \frac{1}{2} \norm{x - y}_2^2 + t_1\norm y_1 .\] A proof can be found in [10]. Therefore, we have \[\begin{align} \abs[auto]{ \langle{p_x(t_1)},{x}\rangle - \langle{p_x(t_2)},{x}\rangle } & = \frac{1}{2} \norm{p_x(t_1) - x}_2^2 - \frac{1}{2} \norm{p_x(t_2) - x}_2^2 \\ & \le t_1 \left( \norm{p_x(t_2)}_1 - \norm{p_x(t_1)}_1 \right) \\ & \le t_1 \norm{p_x(t_2) - p_x(t_1)}_2 , \end{align}\] where we used the Lipschitz continuity of the \(\ell_1\)-norm with constant \(1\) for the last step.
Since \(\norm{\prox{t \norm{\cdot}_1}^{\bbR^{n+1}}(x)}_2 \ge \norm x_\infty - \xi\) for \(t \le \xi\), we get \[\begin{align} \norm{p_x(t_2) - p_x(t_1)}_2 & \le \frac{1}{\norm x_\infty - \xi} \norm{ \prox{ t_1 \norm{\cdot}_1 }^{\bbR^{n+1}}(x) - \prox{ t_2 \norm{\cdot}_1 }^{\bbR^{n+1}}(x) } \\ & = \frac{1}{\norm x_\infty - \xi} \norm{ \prox{ t_1 \norm{\cdot}_1 }^{\bbR^{n+1}}(x) - \prox{ t_1 \norm{\cdot}_1 }^{\bbR^{n+1}}( x + \bar x (t_1 - t_2) ) }_2 \\ & \le \frac{1}{\norm x_\infty - \xi} \norm{ x - ( x + \bar x (t_1 - t_2) ) }_2 \\ & = \frac{\sqrt{\norm x_0}}{\norm x_\infty - \xi} (t_1 - t_2) , \end{align}\] where \(\bar x\) is defined as \(\bar x_i \coloneq \sgn(x_i)\) and the second to last step uses the non-expansivity of the proximal operator of the \(\ell_1\)-norm, see, , Beck:2017:1?. All in all we get \[\begin{align} \abs[auto]{ \sigma(p_x(t_1)) - \sigma(p_x(t_2)) } & \le \frac{ 4 \xi t_1 \sqrt{\norm x_0} }{(\norm x_\infty - \xi) \pi^2}(t_1-t_2) \\ & \le \frac{ 4 \xi^2 \sqrt{\norm x_0} }{(\norm x_\infty - \xi) \pi^2}(t_1-t_2) . \end{align}\] The upper bound on \(\xi\) in ?? ensures \(\frac{4 \xi^2 \sqrt{\norm x_0}}{(\norm x_\infty - \xi) \pi^2} < 1\). Applying the Banach fixed-point theorem yields the stated result. ◻
5 implies that, if additionally \[\xi < \frac{\pi}{8 \sqrt{\norm x_0}} \left( \sqrt{ 16 \sqrt{\norm x_0} \norm{x}_\infty + \pi^2 } - \pi \right)\] holds, the solution of 24 is unique and can be computed in linear time. Since our method solves the proximal mapping with the parameter \(\xi = \sequence \lambda k \mu\), this bound is never exceeded in practice.
6.1.0.2 Implementation Details
1 We tested the NCRPG method using a data matrix \(A \in \mathbb{R}^{m \times n}\), where \(m = 20\), that was generated with normally distributed random entries. Then, the columns were centered and normalized by their standard deviation.
We ran our method as in 1 using the above fixed point iteration for the computation of the prox problem. As a comparison, we tested the Riemannian Proximal Gradient method (RPG) from [9] and the Manifold Proximal Gradient method (ManPG) [8]. The RPG and ManPG solve a different proximal problem stated over the tangent space. At each iteration, the RPG aims to find a stationary point \(\eta^*_{\sequence{x}{k}} \in \tangentSpace{\sequence{x}{k}}[\text{OB}(n,r)]\) of the function \[\ell_{\sequence{x}{k}}(\eta) = \riemannian[auto] {\grad g\left(\sequence{x}{k}\right)} {\eta} + \frac{1}{2\lambda} \norm{\eta}^2 + h(\retract{\sequence{x}{k}}(\eta)) ,\] s.t. \(\ell_{\sequence{x}{k}}(\eta^*_{\sequence{x}{k}}) \leq \ell_{\sequence{x}{k}}(0)\). Although the problem is stated over a linear space, it remains hard to solve in most cases. However, for the sparse PCA problem stated over the oblique manifold, an explicit iterative algorithm is given in [9] when the exponential map is used as the retraction. ManPG computes a minimizer \(\eta^*_{\sequence{x}{k}} \in \tangentSpace{\sequence{x}{k}}[\text{OB}(n,r)]\) of \(\ell_{\sequence{x}{k}}\) at each iteration but relaxes the retraction to the addition in the embedding space. Additionally, a backtracking procedure with contraction factor \(\beta = 0.5\) is applied to compute the next iterate. As suggested in [8], the subproblem is solved using a regularized semi-smooth Newton method (SSN).
We use \(\lambda = \frac{1}{2 \norm{A}^2_F}\) as a constant stepsize for NCRPG, RPG and ManPG and an initial guess \(s = \frac{5}{\norm{A}^2_F}\) and parameters \(\beta = \eta = 0.5\) in 2 for NCRPG with backtracking. Each run was initialized with a random starting point using the uniform distribution on the oblique manifold. We stopped NCRPG when the norm of the gradient mapping fell below \(10^{-4}\) with the fixed point iteration described above terminating when either \(\vert \sequence{t}{k+1} - \sequence{t}{k} \vert < 10^{-10}\) or after \(10\) iterations (usually after only 2-3 iterations). Similarly, RPG and ManPG were stopped when \(\norm{\frac{1}{\lambda} \eta^*_{\sequence{x}{k}}} < 10^{-4}\). The SSN subsolver of ManPG was terminated when either the residuum fell below \(10^{-12}\) or the number of iterations exceeded \(10\). Additionally, all algorithms were terminated after exceeding \(100\,000\) iterations.
6.1.0.3 Discussion
In 4, we can see that both the norm of the gradient mapping and the error of the function value converge as expected for all four methods. In fact, the local convergence rate seems to be linear. Our method with constant stepsize is a little faster than the RPG with the same stepsize, whereas ManPG took significantly longer. We note that our method also ran with larger stepsizes and was even faster then, but we chose the same stepsize for all methods. The backtracking procedure significantly reduces the number of iterations and it was also faster in terms of CPU time. See 5 for a summary. We note that all methods converge to solutions with a relative difference in the cost function of less than \(0.5\%\). In particular, the solutions obtained by RPG and NCRPG with a constant step size have a relative difference of at most \(10^{-10}\) and yield solutions with identical sparsity. For larger penalty parameter \(\mu\), sparsity is enhanced, as expected. This comes at the cost of either orthogonality or similarity to the exact principal components.
6.2 Experiments on the Real Grassmannian↩︎
The next experiment on the real Grassmannian is particular because here, we do not assume that the manifold is embedded. While [9] also do not need embeddedness in their theory, they at least require that the total space of a quotient manifold is embedded (see [9]), which, in the case of the Grassmannian, would be the Stiefel manifold. Since they do not provide details on how their method would be implemented for the Grassmann manifold, we cannot compare our algorithm to the RPG method in this case.
We aim to solve the optimization problem in 1 with \[g(p) = \frac{1}{2N} \sum_{j=1}^N \dist^2(p, q_j), \quad h(p) = \tau \dist(p, \bar{q}),\] where \(\{q_1, \ldots, q_N\} \subset \cM\) is a set of \(N=1000\) randomly generated data points, \(\bar{q} \in \cM\) is another fixed point, \(\tau \coloneq \frac{1}{2}\) is a weighting parameter, and \(\cM = \grassmann{n}{r}\) is the Grassmann manifold of \(r\)-dimensional subspaces in \(\bbR^n\). The function \(g\) represents a generalization of the sum of squared distances to data points, whose Riemannian gradient is \(\grad g(p) = -\frac{1}{N} \sum_{j=1}^N \logarithm{p} q_j\). The function \(h\) acts as a regularization term encouraging proximity to the reference point \(\bar{q}\), and its proximal mapping is given by \(\prox_{\lambda h}(p) = \geodesic<l>{p}{\bar q}( \min\{ \lambda , \dist(p, \bar q) \} ) ,\) where \(\geodesic<l>{p}{\bar q}\) is the unit-speed geodesic that starts from \(p\) in the direction of \(\bar q\). See, , WeinmannDemaretStorath:2014:1? and, , [20] for computational aspects and properties of the Grassmann manifold.
6.2.0.1 Implementation Details
2 The data points were generated by first choosing an anchor point \(q_0 \in \cM\) uniformly at random. Then, we generated \(N+2\) normally distributed random tangent vectors to \(q_0\) with a standard deviation of \(\sigma = 1\), and mapped to \(\cM\) using \(\exponential{q_0}\). Each tangent vector was first scaled down to be of length strictly less than \(\frac{\pikmax}{2}\). This ensures that the data points thus generated are contained in a uniquely geodesically convex set \(\cU \subset \cM\). Finally, we selected the first \(N\) points as the data points \(\{q_1, \ldots, q_N\}\) for the function \(g\), and the last two points as the reference point \(\bar{q}\) for \(h\), and the initial point \(p_0\) of the optimization respectively. Since \(\cU\) is a uniquely geodesically convex set containing all data points \(q_j\), the function \(g\) is \(\zeta_{1, \kmin}(D)\)-smooth over \(\cU\), where \(D = \diam(\cU)\), see, , [16]. Hence, because \(\kmax = 2\), and \(\kmin = 0\) on the Grassmannian, we can use the constant \(\lipgrad = \zeta_{1, \kmin}(D) = 1\).
max width=
We compared 1 with both a constant and a backtracked stepsize. The constant stepsize was estimated on \(\cU\) in accordance with the prescription in 4.2 with a choice of \(\delta = 0.01\) for \(\lambdadelta\) and \(\zetadelta\). The backtracking procedure was implemented as described in 2 with an initial guess \(s = 1\) and parameters \(\beta = \frac{\zetadelta}{4}\), \(\eta = 0.5\). The results are summarized in ¿tbl:table:NCRPG-CPPA-comparisons-time-iterations?. Both methods were terminated when the gradient mapping norm dropped below \(10^{-7}\). For every experiment, we report the manifold dimension and compare number of iterations and runtime of NCRPG with a constant stepsize and with backtracking.
6.2.0.2 Discussion
The results show that the algorithm clearly works in the nonembedded case. It can be seen that NCRPG with a backtracked stepsize outperforms the constant stepsize variant in all cases.
6.3 Row-sparse Low-rank Matrix recovery↩︎
A common problem in signal recovery applications such as hyperspectral imaging [21] or blind deconvolution [22] is the task of recovering a signal \(\bar X \in \bbR^{M \times N}\) from a set of \(m\) observations that is simultaneously of low rank and admits a row-sparse structure.
Row-sparsity can be achieved by adding a \(\ell_{1,2}\) penalty (the sum of the \(\ell_2\) norm of the rows) to the cost function. Instead of minimizing the nuclear norm, as is often done, we assume the rank \(r\) of the solution to be known. We get an optimization problem over the manifold of rank \(r\) matrices \(\cM_r\) [22]: \[\argmin_{X \in \cM_r} \, \norm{\boldsymbol{A} (X) - y}_2^2 + \mu \norm{X}_{1,2} ,\] where \(\boldsymbol{A} \colon \mathbb{R}^{M \times N} \rightarrow \mathbb{R}^m\) is a linear measurement operator that is often assumed to fulfill the restricted isometry property (RIP). The goal is to recover the (low-rank and row-sparse) signal \(\bar X\) from as few measurements \(y \coloneq \boldsymbol{A} (\bar X)\) as possible.
A convenient way to express a fixed rank matrix \(X \in \cM_r\) is in the form of its compact SVD, , \(X = U \Sigma V^\top\), where \(U^\top U = V^\top V = I_r\) and \(\Sigma\) is the diagonal matrix containing the singular values of \(X\). Using this notation, the tangent space at \(X\) can be expressed as Boumal:2023:1? \[\begin{align} \tangentSpace{X}[\cM_r] = \setDef{ U M V^\top + U_p V^\top + U V_p^\top }{ & M \in \bbR^{r \times r}, U_p \in \bbR^{M \times r}, V_p \in \bbR^{N \times r} \\ & \text{s.t. } U^\top U_p = 0, V^\top V_p = 0 } . \end{align}\] The problem of low-rank and row-sparse matrix recovery does not fulfill all assumptions in 1: The sectional curvature is bounded only locally and the sublevel sets are not compact since their intersection with the manifold is not closed. However, locally, this should not be a problem, which is why we will test only local convergence.
Furthermore, to our knowledge, there is no closed form for geodesics on the fixed-rank manifold. Consequently, the exponential and logarithmic maps cannot be used for optimization over \(\cM_r\). In our experiment, we will therefore use retractions (see 5) instead. As explained there, the convergence of a retraction based NCRPG remains an open problem. A detailed analysis of different retraction and inverse retraction methods, including their computation time and error compared to the exponential map, can be found in AbsilOseledets:2014:1?. In the experiment, we therefore use 3.
We will thus use the orthographic retraction which, for a given tangent vector expressed as \({H = UMV^\top + U_pV^\top + UV_p^\top}\), is given by \[\retract{X}(H) = \left( U(\Sigma + M) + U_p \right) \left( \Sigma + M \right)^{-1} \left( \left( \Sigma + M \right) V^\top + V_p^\top \right) .\] For \(Y \in \cM_r\), the corresponding inverse orthographic retraction is simply: \[\inverseRetract{X}(Y) = \proj{\tangentSpace{X}[\cM_r]}(Y - X) ,\] with the projection onto the tangent space at \(X\) being given by \[\proj{\tangentSpace{X}[\cM_r]}(Z) = (I_M - U U^\top) Z V V^\top + U U^\top Z .\]
We derive a closed-form solution for the proximal subproblem 22 with \(h = \mu \lVert \cdot \rVert_{1,2}\) for \(\lambda > 0\) at \(X\) using the retraction distance, i.e. finding a stationary point of \[H^R(Y;X,\lambda) = \Bigl\| \proj{\tangentSpace{Y}[\cM_r]}\bigl( \mathrm{R}_{X}(-\lambda\grad g(X)) - Y \bigr) \Bigr\|_F^2 + \mu \lVert Y \rVert_{1,2}\] In particular, 3 can be viewed as an generalization of the Riemannian proximal gradient method proposed in [22] for this specific problem. In that method, a proximal update using the Euclidean distance is solved in each iteration \(k\): \[\begin{align} \label{equation:prox-update-fixed-rank-eucl-dist} &\proxR[big]{\lambda \mu(X) \lVert \cdot \rVert_{1,2}}[\retractionSymbol]( \mathrm{R}_{X}(-\lambda\grad g(X)) ) \\ & \qquad\qquad = \underset{Y \in \mathcal{M}_r}{\mathrm{argmin}} \; \frac{1}{2\lambda} \lVert Y - \mathrm{R}_{X}(-\lambda\grad g(X)) \rVert_{\mathrm{F}}^{2} + \mu(X) \lVert Y \rVert_{1,2}, \end{align}\tag{28}\] where \(\mu(X) = \tau^k \max_r (\Vert \mathrm{R}_{X}(-\lambda\grad g(X))_i \Vert \colon i = 1, \ldots, M)\), the norm of the \(r\)-th largest row of the current iterate, with \(\tau \in (0, 1)\), is generally suggested. We will see that \(\proxR[big]{\lambda \mu \lVert \cdot \rVert_{1,2}}[\retractionSymbol]( \mathrm{R}_{X}(-\lambda\grad g(X)))\) is also a viable update step for our method, and thus, the two methods coincide when \(\mu(X) = \mu\) is chosen. For \(Z \in \cM_r\) , let \(D^{\lambda \mu}[Z] \in \mathbb{R}^{M \times M}\) be the row-wise soft thresholding operator with \[D^{\lambda\mu}[Z]_{ij} = \begin{cases} \frac{\lVert Z[i,:] \rVert_{2} - \lambda\mu}{\lVert Z[i,:] \rVert_{2}}, & \qquad \text{if} \;i = j \text{ and } \lVert Z[i,:] \rVert_{2} > \lambda\mu, \\ 0, &\qquad \text{else}. \end{cases}\] From [22] we know that \(D^{\lambda\mu}[Z]Z = \proxR[big]{\lambda \mu \lVert \cdot \rVert_{1,2}}[\retractionSymbol](Z)\). Therefore, it holds that \[0 \in \mu \partial \lVert D^{\lambda \mu}[Z]Z \rVert_{1,2} + \frac{1}{\lambda} \proj{\tangentSpace{D^{\lambda \mu}[Z]Z}[\cM_r]} ( D^{\lambda \mu}[Z]Z - Z ) = \partial H^R(Y;X,\lambda) ,\] where the equality is due to \(Z = U\Sigma V^\top\) and \(D^{\lambda \mu}[Z]Z = D^{\lambda \mu}[Z]U\Sigma V^\top\) lying in the tangent space of \(D^{\lambda \mu}[Z]Z\) (because they have the same row-space \(V\)). To see that this holds, one can compute the Euclidean differential of \(f(Y) = \frac{1}{2} \left\| \proj{\tangentSpace{Y}[\cM_r]}\left(Y-Z \right) \right\|_F^2\) as follows. \[\begin{align} \d f(Y)[\xi] &= \langle \proj{\tangentSpace{Y}[\cM_r]}\left(Y-Z\right) , \proj{\tangentSpace{Y}[\cM_r]} \xi + \d (Y \mapsto \proj{\tangentSpace{Y}[\cM_r]}) \left(Y-Z\right) [\xi] \rangle \\ &= \langle \proj{\tangentSpace{Y}[\cM_r]}\left(Y-Z\right) , \xi \rangle \\ & \quad + \langle Y-Z , \proj{\tangentSpace{Y}[\cM_r]} \d (Y \mapsto \proj{\tangentSpace{Y}[\cM_r]}) \circ \proj{\tangentSpace{Y}[\cM_r]} \left(Y-Z\right) [\xi] \rangle \\ &= \langle \proj{\tangentSpace{Y}[\cM_r]}\left(Y-Z\right) , \xi \rangle , \end{align}\] where the second step uses \(Y-Z = \proj{\tangentSpace{Y}[\cM_r]}\left(Y-Z\right)\) and that the projection is self-adjoint. The last step follows from the fact that for an operator \(P\) with \(P=P^2\) the identity \(P (DP) P = 0\) holds. Hence we get the closed form update step \[\sequence{X}{k+1} = D^{\lambda \mu}[\mathrm{R}_{\sequence{X}{k}}(-\lambda\grad g(\sequence{X}{k}))] \mathrm{R}_{\sequence{X}{k}}(-\lambda\grad g(\sequence{X}{k})),\] , the prox step is simply a row-wise soft thresholding of the gradient step.
6.3.0.1 Implementation Details
3 We tested the local convergence behavior of NCRPG versus the Riemannian Alternating Direction Method of Multipliers (RADMM) [6] for the problem described in 6.3. This algorithm replaces the first prox-step (of the smooth part \(g\)) in the Euclidean ADMM Beck:2017:1? with a gradient step on the manifold, whereas the second prox-step (of the nonsmooth part \(h\)) remains in the embedding space. As above, we used our closed form of the prox-step for \(h\).
The measurement operator \(\boldsymbol{A} \in \mathbb{R}^{m \times MN}\) was taken as a matrix with normally distributed entries in \(\mathcal{N}(0, 1/\sqrt{m})\), such that \(\boldsymbol{A}(X) = \boldsymbol{A} \mathrm{vec}(X) = y\). We note that this choice fulfills the RIP condition with very high probability.
To show local convergence of our method, we first generated a random ground truth matrix \({X \in \bbR^{M\times N}}\) of rank \(r\) with \(s\) nonzero rows, which we perturb with a random matrix \(Y \sim \mathcal{N}(0, 1/\sqrt{m})\) and truncate the rank. The stepsize for NCRPG and RADMM with fixed stepsize was set to \(\lambda = 0.25\). Backtracking was done with an initial guess of \(s = 0.5\) and parameters \(\beta = 0.1\), \(\eta = 0.5\) for 2. The RADMM was implemented in and the parameters \(\rho = 0.1, \beta = 10^{-7}\) were obtained by a grid search. We stopped NCRPG when the norm of the gradient mapping with respect to the orthographic retraction fell below \(10^{-7}\) and the RADMM analogously when \(\frac{1}{\eta} \distR(\sequence {X}{k}, \sequence {X}{k-1}) < 10^{-7}\).
max width=
6.3.0.2 Discussion
Despite the problem not fulfilling all assumptions, we observe (local) convergence in all cases, see ¿tbl:table:NCRPG-RADMM-comparisons-fixed-rank?. The backtracking procedure is slower than NCRPG with a constant stepsize, because more costly evaluations of the objective are required. All algorithms find the row-sparsity pattern of the ground truth, but the RADMM is less confident: the average norm of those rows is higher.
We note that the number of measurements \(m\) is chosen such that the problem is underdetermined if either of the constraints (low rank, row-sparse) is considered alone. Our results show that we recover a unique solution, indicating that, at least locally, the problem is well-posed.
7 Conclusion & Outlook↩︎
In this article, we introduced a nonconvex Riemannian proximal gradient method that works intrinsically on a Riemannian manifold. We established state-of-the-art convergence results with only mild assumptions. Convergence theory for an analogous method that uses retractions as an approximation of the exponential map remains an open problem. In the experiments, we have demonstrated the efficiency of our method in many nonconvex and intrinsic cases that have so far been out of reach for nonsmooth Riemannian optimization.
In an upcoming work [23], we will consider the same method under stronger convexity assumptions, both on the functions \(g\) and \(h\), as well as on the manifolds, considering in particular Hadamard manifolds. We will show that under these stronger assumptions, we can obtain improved convergence rates and even show convergence to stationary points.
Acknowledgements↩︎
We would like to thank Lukas Klingbiel for his help with the low-rank row-sparse matrix recovery problem. P.J. was funded by the DFG – Projektnummer 448293816.
References↩︎
The code is available at https://juliamanifolds.github.io/ManoptExamples.jl/stable/examples/NCRPG-Sparse-PCA/↩︎
The code is available at https://juliamanifolds.github.io/ManoptExamples.jl/stable/examples/NCRPG-Grassmann/↩︎
The code is available at https://juliamanifolds.github.io/ManoptExamples.jl/stable/examples/NCRPG-Row-Sparse-Low-Rank/↩︎