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 with vertices consists of three geodesics . Given , a comparison triangle in -plane is a corresponding triangle with the same side lengths in two-dimensional space of constant Gaussian curvature . A length space with curvature bound is called an Alexandrov space. In particular, we have the following important definition:
Let be a real number. A length space is a space of curvature if every point has a neighborhood such that for any triangle contained in and any point the inequality holds, where is a comparison triangle in the -plane and is the point such that .
The notion of angle is defined in the following sense. Let and be two geodesics in with , we define the angle between and as where is the angle at of the corresponding triangle . 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 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 is a real smooth manifold equipped with an inner product on the tangent space of every point , such that if are two vector fields on then 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 at of a geodesic is still a tangent vector of after being transported to a point along . Furthermore, parallel transport preserves inner products, i.e. .
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 is defined on a Riemannian manifold , unless stated otherwise.
It can be shown that an equivalent definition is that for any ,
where is a subgradient of at , or the gradient if is differentiable, and denotes the inner product in the tangent space of 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 is the parallel transport from to .
Observe that compared to the Euclidean setup, the above definition requires a parallel transport operation to “transport” to . It can be proved that if is -smooth, then for any ,
Convergence Rates of First-order Methods
General subgradient / gradient algorithms on Riemannian manifolds take the form
where is the iterate index, is a subgradient of the objective function, and 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 are the sides of a Euclidean triangle with the angle between sides and , is an essential tool for bounding the squared distance between the iterates and the minimizer(s). Indeed, consider the Euclidean update . Applying (3) to the triangle , with , , , and , 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 be the linearization of the convex function , and let . Then, is the unique solution to the following minimization problem
Since 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 and , the function
is geodesically both star-concave and star-convex in , 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 , 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 are the sides (i.e., side lengths) of a geodesic triangle in an Alexandrov space with curvature lower bounded by , and is the angle between sides and , 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 , which corresponds to (4), since Eqn. (6) can be seen as Eqn. (4) with side lengths (see Lemma A.7).
Finally, we observe that in (4), . Letting , where is the right hand side of (5), we then see that it is sufficient to prove the following:
and are equal at .
the first partial derivatives of and w.r.t. agree at .
for (Lemma A.1).
These three steps, if true, lead to the proof of for , thus proving a special case of Lemma 3.1 for space with constant sectional curvature 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 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 ):
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 rate of convergence for g-convex on Hadamard manifolds.
Summing over from to and dividing by , we obtain
Plugging in and we further obtain
It remains to show that , 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 we can recover the Euclidean convergence rates (in some cases up to a difference in small constant factors). However, for Hadamard manifolds 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 on the manifold.
The proof structure is very similar, except that for each equation we take expectation with respect to the sequence . Since 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 rate of convergence for g-strongly convex functions on Hadamard manifolds.
If is geodesically -strongly convex and -Lipschitz, and the sectional curvature of the manifold is lower bounded by , then the subgradient method with satisfies
Since is geodesically -strongly convex, we have
which combined with Corollary 3.4 and -Lipschitz condition yields
Multiply (12) by and sum over from to ; then divide the result by to obtain
The final step is to show , 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 . 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 , which implies that with 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 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 rate of convergence, whereas stochastic gradient achieves a curvature dependent rate for smooth g-convex functions on Hadamard manifolds.
For simplicity we denote . First observe that with the algorithm is a descent method. Indeed, we have
On the other hand, by the convexity of and Corollary 3.4 we obtain
Multiplying (14) by and adding to (15), we get
Now summing over from to , a brief manipulation shows that
Since for we proved , and by assumption , for we get
As before we write . First we observe that
Taking expectation, and letting , we obtain
On the other hand, using convexity of and Corollary 3.4 we get
Multiply (21) by and add to (22), we get
Summing over from to and simplifying, we obtain
Now set , and note that ; thus, for we get
Finally, due to g-convexity of 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 . Observe that with we have descent:
On the other hand, by the strong convexity of and Corollary 3.4 we obtain the bounds
Multiply (24) by and add to (25) to obtain
Let , multiply (27) by and sum over from to , we get
Observe that since , it follows that , 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 does not cancel nicely due to the presence of the curvature term , 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 may also be possible if one replaces 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 symmetric positive definite matrices 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 -strongly convex , enabling the use of geometrically convex optimization algorithms. The full gradient update step is
where each index is drawn uniformly at random from . The step-sizes 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 estimate of always guarantees convergence. We compare the performance of three algorithms that can be applied to this problem:
Gradient descent (GD) with set according to the estimate of the smoothness constant (Theorem 3.18).
Stochastic gradient method for smooth functions (SGD-sm) with 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 set according to the -strong convexity of the loss function (Theorem 3.12).
Our data are random PSD matrices generated using the Matrix Mean Toolbox (Bini and Iannazzo, 2013). All matrices are explicitly normalized so that their norms all equal . We compare the algorithms on four datasets with matrices to average and the condition number of each matrix being either or . For all experiments we initialize using the arithmetic mean of the dataset. Figure 2 shows 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 , . Now we focus on the case when . If , Let where . We have
Since , it suffices to prove
Solving for , we get . Since and , is monotonically increasing on . Thus . Note that , thus it suffices to prove
Now fix and notice that as a function of is monotonically decreasing. Therefore its minimum is obtained at , where , i.e. . So it only remains to show
Suppose is twice differentiable on with three further assumptions:
,
then
It suffices to prove . We prove this claim by contradiction.
Suppose the claim doesn’t hold, then there exist some so that for any in , and is monotonically increasing in . It follows that for any we have
which leads to a contradiction with our assumption .
If are the sides of a (geodesic) triangle in a hyperbolic space of constant curvature , and is the angle between and , then
For a fixed but arbitrary , define . By Lemma A.1 it is easy to verify that satisfies the assumptions of Lemma A.3. Apply Lemma A.3 to with to show in . Therefore for any . Finally use the fact that is monotonically increasing on .
If are the sides of a (geodesic) triangle in a hyperbolic space of constant curvature , and is the angle between and , then
For hyperbolic space of constant curvature , the law of cosines is
which corresponds to the law of cosines of a geodesic triangle in hyperbolic space of curvature with side lengths . Applying Lemma A.5 we thus get