First-order Methods for Geodesically Convex Optimization

Hongyi Zhang, Suvrit Sra

Introduction

Convex optimization is fundamental to numerous areas including machine learning. Convexity often helps guarantee polynomial runtimes and enables robust, more stable numerical methods. But almost invariably, the use of convexity in machine learning is limited to vector spaces, even though convexity per se is not limited to vector spaces. Most notably, it generalizes to geodesically convex metric spaces (Gromov, 1978; Bridson and Haefliger, 1999; Burago et al., 2001), through which it offers a much richer setting for developing mathematical models amenable to global optimization.

Our broader aim is to increase awareness about g-convexity (see Definition 2.2); while our specific focus in this paper is on contributing to the understanding of geodesically convex (g-convex) optimization. In particular, we study first-order algorithms for smooth and nonsmooth g-convex optimization, for which we prove iteration complexity upper bounds. Except for a fundamental lemma that applies to general g-convex metric spaces, we limit our discussion to Hadamard manifolds (Riemannian manifolds with global nonpositive curvature), as they offer the most convenient grounds for generalization while also being relevant to numerous applications (see e.g., Section 1.1).

Specifically, we study optimization problems of the form

Although Riemannian geometry provides tools that enable generalization of Euclidean algorithms (Udriste, 1994; Absil et al., 2009), to obtain iteration complexity bounds we must overcome some fundamental geometric hurdles. We introduce key results that overcome some of these hurdles, and pave the way to analyzing first-order g-convex optimization algorithms.

We recollect below a few items of related work and some examples relevant to machine learning, where g-convexity and more generally Riemannian optimization play an important role.

Standard references on Riemannian optimization are (Udriste, 1994; Absil et al., 2009), who primarily consider problems on manifolds without necessarily having access to g-convexity. Consequently, their analysis is limited to asymptotic convergence (except for (Theorem 4.2, Udriste, 1994) that proves linear convergence for functions with positive-definite and bounded Riemannian Hessians). The recent monograph (Bacák, 2014) is devoted to g-convexity and g-convex optimization on geodesic metric spaces, though without any attention to global complexity analysis. Bacák (2014) also details a noteworthy application: averaging trees in the geodesic metric space of phylogenetic trees (Billera et al., 2001).

At a more familiar level, implicitly the topic of “geometric programming” (Boyd et al., 2007) may be viewed as a special case of g-convex optimization (Sra and Hosseini, 2015). For instance, computing stationary states of Markov chains (e.g., while computing PageRank) may be viewed as g-convex optimization problems by placing suitable geometry on the positive orthant; this idea has a fascinating extension to nonlinear iterations on convex cones (in Banach spaces) endowed with the structure of a geodesic metric space (Lemmens and Nussbaum, 2012).

Perhaps the most important example of such metric spaces is the set of positive definite matrices viewed as a Riemannian or Finsler manifold; a careful study of this setup was undertaken by Sra and Hosseini (2015). They also highlighted applications to maximum likelihood estimation for certain non-Gaussian (heavy- or light-tailed) distributions, resulting in various g-convex and nonconvex likelihood problems; see also (Wiesel, 2012; Zhang et al., 2013). However, none of these three works presents a global convergence rate analysis for their algorithms.

There exist several nonconvex problems where Riemannian optimization has proved quite useful, e.g., low-rank matrix and tensor factorization (Vandereycken, 2013; Ishteva et al., 2011; Mishra et al., 2013); dictionary learning (Sun et al., 2015; Harandi et al., 2012); optimization under orthogonality constraints (Edelman et al., 1998; Moakher, 2002; Shen et al., 2009; Liu et al., 2015); and Gaussian mixture models (Hosseini and Sra, 2015), for which g-convexity helps accelerate manifold optimization to greatly outperform the Expectation Maximization (EM) algorithm.

2 Contributions

We summarize the main contributions of this paper below.

We develop a new inequality (Lemma 3.1) useful for analyzing the behavior of optimization algorithms for functions in Alexandrov space with curvature bounded below, which can be applied to (not necessarily g-convex) optimization problems on Riemannian manifolds and beyond.

For g-convex optimization problems on Hadamard manifold (Riemannian manifold with nonpositive sectional curvature), we prove iteration complexity upper bounds for several existing algorithms (Table 1). For the special case of smooth geodesically strongly convex optimization, a prior linear convergence result that uses line-search is known (Udriste, 1994); our results do not require line search. Moreover, as far as we are aware, ours are the first global complexity results for general g-convex optimization.

Background

Before we describe the algorithms and analyze their properties, we would like to introduce some concepts in metric geometry and Riemannian geometry that generalize concepts in Euclidean space.

The properties of geodesic triangles will be central to our analysis of optimization algorithms. A geodesic triangle △pqr\triangle pqr with vertices p,q,r∈Xp,q,r\in X consists of three geodesics pq‾,qr‾,rp‾\overline{pq},\overline{qr},\overline{rp}. Given △pqr∈X\triangle pqr\in X, a comparison triangle △pˉqˉrˉ\triangle\bar{p}\bar{q}\bar{r} in kk-plane is a corresponding triangle with the same side lengths in two-dimensional space of constant Gaussian curvature kk. A length space with curvature bound is called an Alexandrov space. In particular, we have the following important definition:

Let kk be a real number. A length space XX is a space of curvature ≥k\geq k if every point x∈Xx\in X has a neighborhood UU such that for any triangle △abc\triangle abc contained in UU and any point d∈ac‾d\in\overline{ac} the inequality ∣bd∣≥∣bˉdˉ∣|bd|\geq|\bar{b}\bar{d}| holds, where △aˉbˉcˉ\triangle\bar{a}\bar{b}\bar{c} is a comparison triangle in the kk-plane and dˉ∈aˉcˉ‾\bar{d}\in\overline{\bar{a}\bar{c}} is the point such that ∣aˉdˉ∣=∣ad∣|\bar{a}\bar{d}|=|ad|.

The notion of angle is defined in the following sense. Let γ:→X\gamma:\to X and η:→X\eta:\to X be two geodesics in (X,d)(X,d) with γ0=η0\gamma_{0}=\eta_{0}, we define the angle between γ\gamma and η\eta as α(γ,η):=lim⁡sup⁡s,t→0+∡γˉsγˉ0ηˉt\alpha(\gamma,\eta):=\lim\sup_{s,t\to 0_{+}}\measuredangle\bar{\gamma}_{s}\bar{\gamma}_{0}\bar{\eta}_{t} where ∡γˉsγˉ0ηˉt\measuredangle\bar{\gamma}_{s}\bar{\gamma}_{0}\bar{\eta}_{t} is the angle at γˉ0\bar{\gamma}_{0} of the corresponding triangle △γˉsγˉ0ηˉt\triangle\bar{\gamma}_{s}\bar{\gamma}_{0}\bar{\eta}_{t}. We use Toponogov’s theorem to relate the angles and lengths of any geodesic triangle in a geodesic space to those of a comparison triangle in a space of constant curvature (Burago et al., 1992, 2001).

2 Riemannian Geometry

As tangent vectors at two different points x,y∈Mx,y\in\mathcal{M} lie in different tangent spaces, we cannot compare them directly. To meaningfully compare vectors in different tangent spaces, one needs to define a way to move a tangent vector along the geodesics, while ‘preserving’ its length and orientation. We thus need to use an inner product structure on tangent spaces, which is called a Riemannian metric. A Riemannian manifold (M,g)(\mathcal{M},g) is a real smooth manifold equipped with an inner product gxg_{x} on the tangent space TxMT_{x}\mathcal{M} of every point xx, such that if u,vu,v are two vector fields on M\mathcal{M} then x↦⟨u,v⟩x:=gx(u,v)x\mapsto\langle u,v\rangle_{x}:=g_{x}(u,v) is a smooth function. On a Riemannian manifold, the notion of parallel transport (parallel displacement) provides a sensible way to transport a vector along a geodesic. Intuitively, a tangent vector v∈TxMv\in T_{x}\mathcal{M} at xx of a geodesic γ\gamma is still a tangent vector Γ(γ)xyv\Gamma(\gamma)_{x}^{y}v of γ\gamma after being transported to a point yy along γ\gamma. Furthermore, parallel transport preserves inner products, i.e. ⟨u,v⟩x=⟨Γ(γ)xyu,Γ(γ)xyv⟩y\langle u,v\rangle_{x}=\langle\Gamma(\gamma)_{x}^{y}u,\Gamma(\gamma)_{x}^{y}v\rangle_{y}.

The curvature of a Riemannian manifold is characterized by its Riemannian metric tensor at each point. For worst-case analysis, it is sufficient to consider geodesic triangles of any two-dimensional subspace. Sectional curvature is the Gauss curvature of a two dimensional subspace of a Riemannian manifold, which characterizes the metric space property within that subspace. A subspace with positive, zero or negative sectional curvature is locally isometric to a two dimensional sphere, a Euclidean plane, or a hyperbolic plane with the same Gauss curvature.

3 Function Classes on a Riemannian Manifold

We first define some key terms. Throughout the paper, we assume that the function ff is defined on a Riemannian manifold M\mathcal{M}, unless stated otherwise.

It can be shown that an equivalent definition is that for any x,y∈Mx,y\in\mathcal{M},

where gxg_{x} is a subgradient of ff at xx, or the gradient if ff is differentiable, and ⟨⋅,⋅⟩x\langle\cdot,\cdot\rangle_{x} denotes the inner product in the tangent space of xx induced by the Riemannian metric. In the rest of the paper we will omit the index of tangent space when it is clear from the context.

where Γyx\Gamma_{y}^{x} is the parallel transport from yy to xx.

Observe that compared to the Euclidean setup, the above definition requires a parallel transport operation to “transport” gyg_{y} to gxg_{x}. It can be proved that if ff is LgL_{g}-smooth, then for any x,y∈Mx,y\in\mathcal{M},

Convergence Rates of First-order Methods

General subgradient / gradient algorithms on Riemannian manifolds take the form

where ss is the iterate index, gsg_{s} is a subgradient of the objective function, and ηs\eta_{s} is a step-size. For brevity, we will use the word ‘gradient’ to refer to both subgradient and gradient, deterministic or stochastic; the meaning should be apparent from the context.

While it is easy to translate first-order optimization algorithms from Euclidean space to Riemannian manifolds, and similarly to prove asymptotic convergence rates (since locally Riemannian manifolds resemble Euclidean space), it is much harder to carry out non-asymptotic analysis, at least due to the following two difficulties:

Non-Euclidean trigonometry is difficult to use. Trigonometric geometry in nonlinear spaces is fundamentally different from Euclidean space. In particular, for analyzing optimization algorithms, the law of cosines in Euclidean space

where a,b,ca,b,c are the sides of a Euclidean triangle with AA the angle between sides bb and cc, is an essential tool for bounding the squared distance between the iterates and the minimizer(s). Indeed, consider the Euclidean update xs+1=xs−ηsgsx_{s+1}=x_{s}-\eta_{s}g_{s}. Applying (3) to the triangle △xsxxs+1\triangle x_{s}xx_{s+1}, with a=xxs+1‾a=\overline{xx_{s+1}}, b=xsxs+1‾b=\overline{x_{s}x_{s+1}}, c=xxs‾c=\overline{xx_{s}}, and A=∡xxsxs+1A=\measuredangle xx_{s}x_{s+1}, we get the frequently used formula

However, this nice equality does not exist for nonlinear spaces.

Linearization does not work. Another key technique used in bounding squared distances is inspired by the proximal algorithms. Here, gradient-like updates are seen as proximal steps for minimizing a series of linearizations of the objective function. Specifically, let ψ(x;xs)=f(xs)+⟨gs,x−xs⟩\psi(x;x_{s})=f(x_{s})+\langle g_{s},x-x_{s}\rangle be the linearization of the convex function ff, and let gs∈∂f(xs)g_{s}\in\partial f(x_{s}). Then, xs+1=xs−ηsgsx_{s+1}=x_{s}-\eta_{s}g_{s} is the unique solution to the following minimization problem

Since ψ(x;xs)\psi(x;x_{s}) is convex, we thus have (see e.g. Tseng (2009)) the recursively useful bound

But in nonlinear space there is no trivial analogy of a linear function. For example, for any given y∈My\in\mathcal{M} and gy∈TyMg_{y}\in T_{y}\mathcal{M}, the function

is geodesically both star-concave and star-convex in yy, but neither convex nor concave in general. Thus a nonlinear analogue of the above result does not hold.

We address the first difficulty by developing an easy-to-use trigonometric distance bound for Alexandrov space with curvature bounded below. When specialized to Hadamard manifolds, our result reduces to the analysis in (Bonnabel, 2013), which in turn relies on (Cordero-Erausquin et al., 2001, Lemma 3.12). However, unlike (Cordero-Erausquin et al., 2001), our proof assumes no manifold structure on the geodesic space of interest, and is fundamentally different in techniques.

As noted above, a main hurdle in analyzing non-asymptotic convergence of first-order methods in geodesic spaces is that the Euclidean law of cosines does not hold any more. For general nonlinear spaces, there are no corresponding analytical expressions. Even for the (hyperbolic) space of constant negative curvature −1-1, perhaps the simplest and most studied nonlinear space, the law of cosines is replaced by the hyperbolic law of cosines:

which is not amendable to the standard techniques of convergence rate analysis. With the goal of developing analysis for nonlinear space optimization algorithms, our first contribution is the following trigonometric distance bound for Alexandrov space with curvature bounded below. Owing to its fundamental nature, we believe that this lemma may be of broader interest too.

If a,b,ca,b,c are the sides (i.e., side lengths) of a geodesic triangle in an Alexandrov space with curvature lower bounded by κ\kappa, and AA is the angle between sides bb and cc, then

sketch. The complete proof contains technical details that digress from the main focus of this paper, so we leave them in the appendix. Below we sketch the main steps.

Our first observation is that by the famous Toponogov’s theorem (Burago et al., 1992, 2001), we can upper bound the side lengths of a geodesic triangle in an Alexandrov space with curvature bounded below by the side lengths of a comparison triangle in the hyperbolic plane, which satisfies (cf. (4)):

Second, we observe that it suffices to study κ=−1\kappa=-1, which corresponds to (4), since Eqn. (6) can be seen as Eqn. (4) with side lengths a=∣κ∣a′,b=∣κ∣b′,c=∣κ∣c′a=\sqrt{|\kappa|}a^{\prime},b=\sqrt{|\kappa|}b^{\prime},c=\sqrt{|\kappa|}c^{\prime} (see Lemma A.7).

Finally, we observe that in (4), ∂2∂b2cosh⁡(a)=cosh⁡(a)\frac{\partial^{2}}{\partial b^{2}}\cosh(a)=\cosh(a). Letting g(b,c,A):=cosh⁡(rhs(b,c,A))g(b,c,A):=\cosh(\sqrt{\text{rhs}(b,c,A)}), where rhs(b,c,A)\text{rhs}(b,c,A) is the right hand side of (5), we then see that it is sufficient to prove the following:

cosh⁡(a)\cosh(a) and g(b,c,A)g(b,c,A) are equal at b=0b=0.

the first partial derivatives of cosh⁡(a)\cosh(a) and g(b,c,A)g(b,c,A) w.r.t. bb agree at b=0b=0.

∂2∂b2g(b,c,A)≥g(b,c,A)\frac{\partial^{2}}{\partial b^{2}}g(b,c,A)\geq g(b,c,A) for b,c≥0b,c\geq 0 (Lemma A.1).

These three steps, if true, lead to the proof of cosh⁡(a)≤g(b,c,A)\cosh(a)\leq g(b,c,A) for b,c≥0b,c\geq 0, thus proving a special case of Lemma 3.1 for space with constant sectional curvature −1-1 as shown in Lemma A.3, A.5. Combing this special case with our first two observations concludes the proof of the lemma.

Inequality (5) provides an upper bound on the side lengths of a geodesic triangle in an Alexandrov space with curvature bounded below. Some examples of such spaces are Riemannian manifolds, including hyperbolic space, Euclidean space, sphere, orthogonal groups, and compact sets on a PSD manifold. However, our derivation does not rely on any manifold structure, thus it also applies to certain cones and convex hypersurfaces (Burago et al., 2001).

In the sequel, we use the notation ζ(κ,c)≜∣κ∣ctanh⁡(∣κ∣c)\zeta(\kappa,c)\triangleq\frac{\sqrt{|\kappa|}c}{\tanh(\sqrt{|\kappa|}c)} for the curvature dependent quantity from inequality (5). From Lemma 3.1 it is straightforward to prove the following corollary, which characterizes an important relation between two consecutive updates of an iterative optimization algorithm on Riemannian manfiold with curvature bounded below.

It is instructive to compare (7) with its Euclidean counterpart (for which actually ζ=1\zeta=1):

Corollary 3.4 furnishes the missing tool for analyzing non-asymptotic convergence rates of manifold optimization algorithms. We now move to the analysis of several such first-order algorithms.

2 Convergence Rate Analysis

The following two theorems show that both deterministic and stochastic subgradient methods achieve a curvature-dependent O(1/t)O(1/\sqrt{t}) rate of convergence for g-convex on Hadamard manifolds.

Summing over ss from 11 to tt and dividing by tt, we obtain

Plugging in d(x1,x∗)≤Dd(x_{1},x^{*})\leq D and η=DLfζ(κ,D)t\eta=\frac{D}{L_{f}\sqrt{\zeta(\kappa,D)t}} we further obtain

It remains to show that f(x‾t)≤1t∑s=1tf(xs)f(\overline{x}_{t})\leq\frac{1}{t}\sum_{s=1}^{t}f(x_{s}), which can be proved by an easy induction.

We note that Theorem 3.6 and our following results are all generalizations of known results in Euclidean space. Indeed, setting curvature κ=0\kappa=0 we can recover the Euclidean convergence rates (in some cases up to a difference in small constant factors). However, for Hadamard manifolds κ<0\kappa<0 and the theorem implies that the algorithms may converge more slowly. Also worth noting is that we must be careful in how we obtain the “average” iterate x‾t\overline{x}_{t} on the manifold.

The proof structure is very similar, except that for each equation we take expectation with respect to the sequence {xs}s=1t\{x_{s}\}_{s=1}^{t}. Since ff is g-convex, we have

Now arguing as in Theorem 3.6 the proof follows.

Strongly convex nonsmooth functions.

The following two theorems show that both subgradient method and stochastic subgradient method achieve a curvature dependent O(1/t)O(1/t) rate of convergence for g-strongly convex functions on Hadamard manifolds.

If ff is geodesically μ\mu-strongly convex and LfL_{f}-Lipschitz, and the sectional curvature of the manifold is lower bounded by κ≤0\kappa\leq 0, then the subgradient method with ηs=2μ(s+1)\eta_{s}=\frac{2}{\mu(s+1)} satisfies

Since ff is geodesically μ\mu-strongly convex, we have

which combined with Corollary 3.4 and LfL_{f}-Lipschitz condition yields

Multiply (12) by ss and sum over ss from 11 to tt; then divide the result by t(t+1)2\frac{t(t+1)}{2} to obtain

The final step is to show f(x‾t)≤2t(t+1)∑s=1tsf(xs)f(\overline{x}_{t})\leq\frac{2}{t(t+1)}\sum_{s=1}^{t}sf(x_{s}), which again follows by an easy induction.

The proof structure is very similar to the previous theorem, except that now we take expectations over the sequence {xs}s=1t\{x_{s}\}_{s=1}^{t}. We omit the details for brevity.

Theorems 3.10 and 3.12 are generalizations of their Euclidean counterparts (Lacoste-Julien et al., 2012), and follow the same proof structures. Our upper bounds depend linearly on ζ(κ,D)\zeta(\kappa,D), which implies that with κ<0\kappa<0 the algorithms may converge more slowly. However, note that for strongly convex problems, the distances from iterates to the minimizer are shrinking, thus the inequality (11) (or its stochastic version) may be too pessimistic, and better dependency on κ\kappa may be obtained with a more refined analysis. We leave this as an open problem for the future.

Smooth convex optimization.

The following two theorems show that gradient descent algorithm achieves a curvature dependent O(1/t)O(1/t) rate of convergence, whereas stochastic gradient achieves a curvature dependent O(1/t+1/t)O(1/t+\sqrt{1/t}) rate for smooth g-convex functions on Hadamard manifolds.

For simplicity we denote Δs=f(xs)−f(x∗)\Delta_{s}=f(x_{s})-f(x^{*}). First observe that with η=1Lg\eta=\frac{1}{L_{g}} the algorithm is a descent method. Indeed, we have

On the other hand, by the convexity of ff and Corollary 3.4 we obtain

Multiplying (14) by ζ(κ,D)\zeta(\kappa,D) and adding to (15), we get

Now summing over ss from 11 to t−1t-1, a brief manipulation shows that

Since for s≤ts\leq t we proved Δt≤Δs\Delta_{t}\leq\Delta_{s}, and by assumption Δ1≤LgD22\Delta_{1}\leq\frac{L_{g}D^{2}}{2}, for t>1t>1 we get

As before we write Δs=f(xs)−f(x∗)\Delta_{s}=f(x_{s})-f(x^{*}). First we observe that

Taking expectation, and letting η=1Lg+1/α\eta=\frac{1}{L_{g}+1/\alpha}, we obtain

On the other hand, using convexity of ff and Corollary 3.4 we get

Multiply (21) by ζ(κ,D)\zeta(\kappa,D) and add to (22), we get

Summing over ss from 11 to t−1t-1 and simplifying, we obtain

Now set α=Dσζ(κ,D)t\alpha=\frac{D}{\sigma\sqrt{\zeta(\kappa,D)t}}, and note that Δ1≤LgD22\Delta_{1}\leq\frac{L_{g}D^{2}}{2}; thus, for t>1t>1 we get

Finally, due to g-convexity of ff it is easy to verify by induction that

Smooth and strongly convex functions.

Next we prove that gradient descent achieves a curvature dependent linear rate of convergence for geodesically strongly convex and smooth functions on Hadamard manifolds.

As before we use Δs=f(xs)−f(x∗)\Delta_{s}=f(x_{s})-f(x^{*}). Observe that with η=1Lg\eta=\frac{1}{L_{g}} we have descent:

On the other hand, by the strong convexity of ff and Corollary 3.4 we obtain the bounds

Multiply (24) by ζ(κ,D)\zeta(\kappa,D) and add to (25) to obtain

Let ϵ=min⁡{1ζ(κ,D),μLg}\epsilon=\min\{\frac{1}{\zeta(\kappa,D)},\frac{\mu}{L_{g}}\}, multiply (27) by (1−ϵ)−(s−1)(1-\epsilon)^{-(s-1)} and sum over ss from 11 to t−1t-1, we get

Observe that since Δ1≤LgD22\Delta_{1}\leq\frac{L_{g}D^{2}}{2}, it follows that Δt≤(1−ϵ)t−2LgD22\Delta_{t}\leq\frac{(1-\epsilon)^{t-2}L_{g}D^{2}}{2}, as desired.

It must be emphasized that the proofs of Theorems 3.14, 3.16, and 3.18 contain some additional difficulties beyond their Euclidean counterparts. In particular, the term Δs\Delta_{s} does not cancel nicely due to the presence of the curvature term ζ(κ,D)\zeta(\kappa,D), which necessitates use of a different Lyapunov function to ensure convergence. Consequently, the stochastic gradient algorithm in Theorem 3.16 requires some unusual looking averaging scheme. In Theorem 3.18, since the distance between iterates and the minimizer is shrinking, better dependency on κ\kappa may also be possible if one replaces ζ(κ,D)\zeta(\kappa,D) by a tighter constant.

Experiments

To empirically validate our results, we compare the performance of a stochastic gradient algorithm with a full gradient descent algorithm on the matrix Karcher mean problem. Averaging PSD matrices have applications in averaging data of anisotropic symmetric positive-definite tensors, such as in diffusion tensor imaging (Pennec et al., 2006; Fletcher and Joshi, 2007) and elasticity theory (Cowin and Yang, 1997). The computation and properties of various notions of geometric means have been studied by many (e.g. Moakher (2005); Bini and Iannazzo (2013); Sra and Hosseini (2015)). Specifically, the Karcher mean of a set of NN symmetric positive definite matrices {Ai}i=1N\{A_{i}\}_{i=1}^{N} is defined as the PSD matrix that minimizes the sum of squared distance induced by the Riemannian metric:

is known to be nonconvex in Euclidean space but geometrically 2N2N-strongly convex , enabling the use of geometrically convex optimization algorithms. The full gradient update step is

where each index i(s)i(s) is drawn uniformly at random from {1,…,N}\{1,\dots,N\}. The step-sizes ηs\eta_{s} for gradient descent and stochastic gradient method have to be chosen according to the smoothness constant or the strongly-convex constant of the loss function. Unfortunately, unlike the Euclidean square loss, there is no cheap way to compute the smoothness constant exactly. In (Bini and Iannazzo, 2013) the authors proposed an adaptive procedure to estimate the optimal step-size. Empirically, however, we observe that an LgL_{g} estimate of 5N5N always guarantees convergence. We compare the performance of three algorithms that can be applied to this problem:

Gradient descent (GD) with ηs=15N\eta_{s}=\frac{1}{5N} set according to the estimate of the smoothness constant (Theorem 3.18).

Stochastic gradient method for smooth functions (SGD-sm) with ηs=1N(s+1)\eta_{s}=\frac{1}{N(s+1)} set according to the estimates of the smoothness constant, domain diameter and gradient variance (Theorem 3.16).

Stochastic subgradient method for strongly convex functions (SGD-st) with ηs=1N(s+1)\eta_{s}=\frac{1}{N(s+1)} set according to the 2N2N-strong convexity of the loss function (Theorem 3.12).

Our data are 100×100100\times 100 random PSD matrices generated using the Matrix Mean Toolbox (Bini and Iannazzo, 2013). All matrices are explicitly normalized so that their norms all equal 11. We compare the algorithms on four datasets with N∈{102,103}N\in\{10^{2},10^{3}\} matrices to average and the condition number QQ of each matrix being either 10210^{2} or 10810^{8}. For all experiments we initialize XX using the arithmetic mean of the dataset. Figure 2 shows f(X)−f(X∗)f(X)-f(X^{*}) as a function of number of passes through the dataset. We observe that the full gradient algorithm with fixed step-size achieves linear convergence, whereas the stochastic gradient algorithms have a sublinear convergence rate, but is much faster during the initial steps.

Discussion

In this paper, we make contributions to the understanding of geodesically convex optimization on Hadamard manifolds. Our contributions are twofold: first, we develop a user-friendly trigonometric distance bound for Alexandrov space with curvature bounded below, which includes several commonly known Riemannian manifolds as special cases; second, we prove iteration complexity upper bounds for several first-order algorithms on Hadamard manifolds, which are the first such analyses up to the best of our knowledge. We believe that our analysis is a small step, yet in the right direction, towards understanding and realizing the power of optimization in nonlinear spaces.

Many questions are not yet answered. We summarize some important ones in the following:

A long-time question is whether the famous Nesterov’s accelerated gradient descent algorithms have nonlinear space counterparts. The analysis of Nesterov’s algorithms typically relies on a proximal gradient projection interpretation. In nonlinear space, we have not been able to find an analogy to such a projection. Further study is needed to see if similar analysis can be developed, or a different approach is required, or Nesterov’s algorithms have no nonlinear space counterparts.

Another interesting direction is variance reduced stochastic gradient methods for geodesically convex functions. For smooth and convex optimization in Euclidean space, these methods have recently drawn great interests and enjoyed remarkable empirical success. We hypothesize that similar algorithms can achieve faster convergence over naive incremental gradient methods on Hadamard manifolds.

Finally, since in applications it is often favorable to replace exponential mapping with computationally cheap retractions, it is important to understand the effect of this approximation on convergence rate. Analyzing this effect is of both theoretical and practical interests.

References

Appendix A Proof of Lemma 1

If c=0c=0, g(b,c)=cosh⁡(b)=∂2∂b2g(b,c)g(b,c)=\cosh(b)=\frac{\partial^{2}}{\partial b^{2}}g(b,c). Now we focus on the case when c>0c>0. If c>0c>0, Let u=(1+x)b2+c2−2bccos⁡(A)u=\sqrt{(1+x)b^{2}+c^{2}-2bc\cos(A)} where x=x(c)x=x(c). We have

Since g(b,c)=cosh⁡(u)>0g(b,c)=\cosh(u)>0, it suffices to prove

Solving for h1′(u)=0h^{\prime}_{1}(u)=0, we get u=0u=0. Since lim⁡u→0+h1(u)=0\lim_{u\to 0_{+}}h_{1}(u)=0 and h1(u)>0,∀u>0h_{1}(u)>0,\forall u>0, h1(u)h_{1}(u) is monotonically increasing on u>0u>0. Thus h1(u)≥h1(umin⁡),∀u>0h_{1}(u)\geq h_{1}(u_{\min}),\forall u>0. Note that c2x(x+sin⁡2A)=1+xxumin⁡2\frac{c^{2}}{x}(x+\sin^{2}A)=\frac{1+x}{x}u^{2}_{\min}, thus it suffices to prove

Now fix cc and notice that tanh⁡(umin⁡)umin⁡\frac{\tanh(u_{\min})}{u_{\min}} as a function of sin⁡2A\sin^{2}A is monotonically decreasing. Therefore its minimum is obtained at sin⁡2A=1\sin^{2}A=1, where umin⁡2=u∗2=c2u^{2}_{\min}=u^{2}_{*}=c^{2}, i.e. u∗=cu_{*}=c. So it only remains to show

Suppose h(x)h(x) is twice differentiable on [r,+∞)[r,+\infty) with three further assumptions:

h′′(x)≤h(x),∀x∈[r,+∞)h^{\prime\prime}(x)\leq h(x),\forall x\in[r,+\infty),

then h(x)≤0,∀x∈[r,+∞)h(x)\leq 0,\forall x\in[r,+\infty)

It suffices to prove h′(x)≤0,∀x∈[r,+∞)h^{\prime}(x)\leq 0,\forall x\in[r,+\infty). We prove this claim by contradiction.

Suppose the claim doesn’t hold, then there exist some t>s≥rt>s\geq r so that h′(x)≤0h^{\prime}(x)\leq 0 for any xx in [r,s][r,s], h′(s)=0h^{\prime}(s)=0 and h′(x)>0h^{\prime}(x)>0 is monotonically increasing in (s,t](s,t]. It follows that for any x∈[s,t]x\in[s,t] we have

which leads to a contradiction with our assumption h′(t)>0h^{\prime}(t)>0.

If a,b,ca,b,c are the sides of a (geodesic) triangle in a hyperbolic space of constant curvature −1-1, and AA is the angle between bb and cc, then

For a fixed but arbitrary c≥0c\geq 0, define hc(x)=f(x,c)−g(x,c)h_{c}(x)=f(x,c)-g(x,c). By Lemma A.1 it is easy to verify that hc(x)h_{c}(x) satisfies the assumptions of Lemma A.3. Apply Lemma A.3 to hch_{c} with r=0r=0 to show hc≤0h_{c}\leq 0 in [0,+∞)[0,+\infty). Therefore f(b,c)≤g(b,c)f(b,c)\leq g(b,c) for any b,c≥0b,c\geq 0. Finally use the fact that cosh⁡(x)\cosh(x) is monotonically increasing on [0,+∞)[0,+\infty).

If a,b,ca,b,c are the sides of a (geodesic) triangle in a hyperbolic space of constant curvature κ\kappa, and AA is the angle between bb and cc, then

For hyperbolic space of constant curvature κ<0\kappa<0, the law of cosines is

which corresponds to the law of cosines of a geodesic triangle in hyperbolic space of curvature −1-1 with side lengths ∣κ∣a,∣κ∣b,∣κ∣c\sqrt{|\kappa|}a,\sqrt{|\kappa|}b,\sqrt{|\kappa|}c. Applying Lemma A.5 we thus get