Optimization Techniques on Riemannian Manifolds

Introduction

The preponderance of optimization techniques address problems posed on Euclidean spaces. Indeed, several fundamental algorithms have arisen from the desire to compute the minimum of quadratic forms on Euclidean space. However, many optimization problems are posed on non-Euclidean spaces. For example, finding the largest eigenvalue of a symmetric matrix may be posed as the maximization of Rayleigh’s quotient defined on the sphere. Optimization problems subject to nonlinear differentiable equality constraints on Euclidean space also lie within this category. Many optimization problems share with these examples the structure of a differentiable manifold endowed with a Riemannian metric. This is the subject of this paper: the extremization of functions defined on Riemannian manifolds.

The minimization of functions on a Riemannian manifold is, at least locally, equivalent to the smoothly constrained optimization problem on a Euclidean space, because every $C^{\infty}$ Riemannian manifold can be isometrically imbedded in some Euclidean space [47, Vol. V]. However, the dimension of the Euclidean space may be larger than the dimension of the manifold; practical and aesthetic considerations suggest that one try to exploit the intrinsic structure of the manifold. Elements of this spirit may be found throughout the field of numerical methods, such as the emphasis on unitary (norm preserving) transformations in numerical linear algebra , or the use of feasible direction methods .

An intrinsic approach leads one from the extrinsic idea of vector addition to the exponential map and parallel translation, from minimization along lines to minimization along geodesics, and from partial differentiation to covariant differentiation. The computation of geodesics, parallel translation, and covariant derivatives can be quite expensive. For an $n$-dimensional manifold, the computation of geodesics and parallel translation requires the solution of a system of $2n$ nonlinear and $n$ linear ordinary differential equations. Nevertheless, many optimization problems are posed on manifolds that have an underlying algebraic structure that may be exploited to greatly reduce the complexity of these computations. For example, on a real compact semisimple Lie group endowed with its natural Riemannian metric, geodesics and parallel translation may be computed via matrix exponentiation . Several algorithms are available to perform this computation . This algebraic structure may be found in the problems posed by Brockett , Bloch et al. , Smith , Faybusovich , Lagarias , Chu et al. , Perkins et al. , and Helmke . This approach is also applicable if the manifold can be identified with a symmetric space or, excepting parallel translation, a reductive homogeneous space . Perhaps the simplest nontrivial example is the sphere, where geodesics and parallel translation can be computed at low cost with trigonometric functions and vector addition. Furthermore, Brown and Bartholomew-Biggs show that in some cases function minimization by following the solution of a system of ordinary differential equations can be implemented such that it is competitive with conventional techniques.

The outline of the paper is as follows. In Section 2, the optimization problem is posed and conventions to be held throughout the paper are established. The method of steepest descent on a Riemannian manifold is described in Section 3. To fix ideas, a proof of linear convergence is given. The examples of Rayleigh’s quotient on the sphere and the function $\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$ on the special orthogonal group are presented. In Section 4, Newton’s method on a Riemannian manifold is derived. As in Euclidean space, this algorithm may be used to compute the extrema of differentiable functions. It is proved that this method converges quadratically. The example of Rayleigh’s quotient is continued, and it is shown that Newton’s method applied to this function converges cubically, and is approximated by the Rayleigh quotient iteration. The example considering $\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$ is continued. In a related example, it is shown that Newton’s method applied to the sum of the squares of the off-diagonal elements of a symmetric matrix converges cubically. This provides an example of a cubically convergent Jacobi-like method. The conjugate gradient method is presented in Section 5 with a proof of superlinear convergence. This technique is shown to provide an effective algorithm for computing the extreme eigenvalues of a symmetric matrix. The conjugate gradient method is applied to the function $\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$.

Preliminaries

This paper is concerned with the following problem.

Let $M$ be a complete Riemannian manifold, and $f$ a $C^{\infty}$ function on $M$. Compute

There are many well-known algorithms for solving this problem in the case where $M$ is a Euclidean space. This paper generalizes several of these algorithms to the case of complete Riemannian manifolds by replacing the Euclidean notions of straight lines and ordinary differentiation with geodesics and covariant differentiation. These concepts are reviewed in the following paragraphs. We follow Helgason’s and Spivak’s treatments of covariant differentiation, the exponential map, and parallel translation. Details may be found in these references.

Let $M$ be a complete $n$-dimensional Riemannian manifold with Riemannian structure $g$ and corresponding Levi-Civita connection $\nabla$. Denote the tangent plane at $p$ in $M$ by $T_{p}$ or $T_{p}M$. For every $p$ in $M$, the Riemannian structure $g$ provides an inner product on $T_{p}$ given by the nondegenerate symmetric bilinear form $g_{p}\colon T_{p}\times T_{p}\to{\bf R}$. The notation $\langle X,Y\rangle=g_{p}(X,Y)$ and $\|X\|=g_{p}(X,X)^{1/2}$, where $X$, $Y\in T_{p}$, is often used. The distance between two points $p$ and $q$ in $M$ is denoted by $d(p,q)$. The gradient of a real-valued $C^{\infty}$ function $f$ on $M$ at $p$, denoted by $(\mathop{\rm grad}\nolimits{\!f})_{p}$, is the unique vector in $T_{p}$ such that $df_{p}(X)=\langle(\mathop{\rm grad}\nolimits{\!f})_{p},X\rangle$ for all $X$ in $T_{p}$.

Denote the set of $C^{\infty}$ functions on $M$ by $C^{\infty}(M)$ and the set of $C^{\infty}$ vector fields on $M$ by ${X}(M)$. An affine connection on $M$ is a function $\nabla$ which assigns to each vector field $X\in{X}(M)$ an ${\bf R}$-linear map $\nabla_{\!X}\colon{X}(M)\to{X}(M)$ which satisfies

for all $f\!$, $g\in C^{\infty}(M)$, $X$, $Y\in{X}(M)$. The map $\nabla_{\!X}$ may be applied to tensors of arbitrary type. Let $\nabla$ be an affine connection on $M$ and $X\in{X}(M)$. Then there exists a unique ${\bf R}$-linear map $A\mapsto\nabla_{\!X}A$ of $C^{\infty}$ tensor fields into $C^{\infty}$ tensor fields which satisfies

where $f\in C^{\infty}(M)$, $Y\in{X}(M)$, and $A$, $B$ are $C^{\infty}$ tensor fields. If $A$ is of type $(k,l)$, then $\nabla_{\!X}A$, called the covariant derivative of $A$ along $X$, is of type $(k,l)$, and ${\nabla\!A}\colon X\mapsto\nabla_{\!X}A$, called the covariant differential of $A$, is of type $(k,l+1)$.

Let $M$ be a differentiable manifold with affine connection $\nabla$. Let $\gamma\colon I\to M$ be a smooth curve with tangent vectors $X(t)=\dot{\gamma}(t)$, where $I\subset{\bf R}$ is an open interval. The curve $\gamma$ is called a geodesic if $\nabla_{\!X}X=0$ for all $t\in I$. Let $Y(t)\in T_{\gamma(t)}$ ($t\in I$) be a smooth family of tangent vectors defined along $\gamma$. The family $Y(t)$ is said to be parallel along $\gamma$ if $\nabla_{\!X}Y=0$ for all $t\in I$.

For every $p$ in $M$ and $X\neq 0$ in $T_{p}$, there exists a unique geodesic t\mapsto\gamma_{\lower 1.0pt\hbox{\scriptstyle X}}(t) such that \gamma_{\lower 1.0pt\hbox{\scriptstyle X}}(0)=p and \dot{\gamma}_{\lower 1.0pt\hbox{\scriptstyle X}}(0)=X. We define the exponential map $\exp_{p}\colon T_{p}\to M$ by \exp_{p}(X)=\gamma_{\lower 1.0pt\hbox{\scriptstyle X}}(1) for all $X\in T_{p}$ such that $1$ is in the domain of \gamma_{\lower 1.0pt\hbox{\scriptstyle X}}. Oftentimes the map $\exp_{p}$ will be denoted by “$\exp$” when the choice of tangent plane is clear, and \gamma_{\lower 1.0pt\hbox{\scriptstyle X}}(t) will be denoted by $\exp tX$. A neighborhood $N_{p}$ of $p$ in $M$ is a normal neighborhood if $N_{p}=\exp N_{0}$, where $N_{0}$ is a star-shaped neighborhood of the origin in $T_{p}$ and $\exp$ maps $N_{0}$ diffeomorphically onto $N_{p}$. Normal neighborhoods always exist.

Given a curve $\gamma\colon I\to M$ such that $\gamma(0)=p$, for each $Y\in T_{p}$ there exists a unique family $Y(t)\in T_{\gamma(t)}$ ($t\in I$) of tangent vectors parallel along $\gamma$ such that $Y(0)=Y$. If $\gamma$ joins the points $p$ and $\gamma(\alpha)=q$, the parallelism along $\gamma$ induces an isomorphism $\tau_{pq}\colon T_{p}\to T_{q}$ defined by $\tau_{pq}Y=Y(\alpha)$.

Given a Riemannian structure $g$ on $M$, there exists a unique affine connection $\nabla$ on $M$, called the Levi-Civita connection, which for all $X$, $Y\in{X}(M)$ satisfies

Length minimizing curves on $M$ are geodesics of the Levi-Civita connection. We shall use this connection throughout the paper.

Unless otherwise specified, all manifolds, vector fields, and functions are assumed to be smooth. When considering a function $f$ to be minimized, the assumption that $f$ is differentiable of class $C^{\infty}$ can be relaxed throughout the paper, but $f$ must be continuously differentiable at least beyond the derivatives that appear. As the results of this paper are local ones, the assumption that $M$ be complete may also be relaxed in certain instances.

We will use the the following definitions to compare the convergence rates of various algorithms.

Let $\{p_{i}\}$ be a Cauchy sequence in $M$ that converges to ${\hat{p}}$. (i) The sequence $\{p_{i}\}$ is said to converge (at least) linearly if there exists an integer $N$ and a constant $\theta\in[0,1)$ such that $d(p_{i+1},{\hat{p}})\leq\theta d(p_{i},{\hat{p}})$ for all $i\geq N$. (ii) The sequence $\{p_{i}\}$ is said to converge (at least) quadratically if there exists an integer $N$ and a constant $\theta\geq 0$ such that $d(p_{i+1},{\hat{p}})\leq\theta d^{2}(p_{i},{\hat{p}})$ for all $i\geq N$. (iii) The sequence $\{p_{i}\}$ is said to converge (at least) cubically if there exists an integer $N$ and a constant $\theta\geq 0$ such that $d(p_{i+1},{\hat{p}})\leq\theta d^{3}(p_{i},{\hat{p}})$ for all $i\geq{N}$. (iv) The sequence $\{p_{i}\}$ is said to converge superlinearly if it converges faster than any sequence that converges linearly.

Steepest descent on Riemannian manifolds

The method of steepest descent on a Riemannian manifold is conceptually identical to the method of steepest descent on Euclidean space. Each iteration involves a gradient computation and minimization along the geodesic determined by the gradient. Fletcher , Botsaris , and Luenberger describe this algorithm in Euclidean space. Gill and Murray and Sargent apply this technique in the presence of constraints. In this section we restate the method of steepest descent described in the literature and provide an alternative formalism that will be useful in the development of Newton’s method and the conjugate gradient method on Riemannian manifolds.

Let $M$ be a complete Riemannian manifold with Riemannian structure $g$ and Levi-Civita connection $\nabla$, and let $f\in C^{\infty}(M)$.

Select $p_{0}\in M$, compute $G_{0}=-(\mathop{\rm grad}\nolimits{\!f})_{p_{0}}$, and set $i=0$.

It is easy to verify that $\langle G_{i+1},\tau G_{i}\rangle=0$, for $i\geq 0$, where $\tau$ is the parallelism with respect to the geodesic from $p_{i}$ to $p_{i+1}$. By assumption, the function $\lambda\mapsto f(\exp\lambda G_{i})$ is minimized at $\lambda_{i}$. Therefore, we have $0={(d/dt)|_{t=0}}\penalty 0{f(\exp(\lambda_{i}+t)G_{i})}=df_{p_{i+1}}(\tau G_{i})=\langle(\mathop{\rm grad}\nolimits{\!f})_{p_{i+1}},\tau G_{i}\rangle$. Thus the method of steepest descent on a Riemannian manifold has the same deficiency as its counterpart on a Euclidean space, i.e., it makes a ninety degree turn at every step.

The convergence of Algorithm 3.1 is linear. To prove this fact, we will make use of a standard theorem of the calculus, expressed in differential geometric language. The covariant derivative $\nabla_{\!X}{\!f}$ of $f$ along $X$ is defined to be $X{\!f}$. For $k=1$, $2$, …, define $\nabla_{\!X}^{k}{\!f}=\nabla_{\!X}\circ\cdots\circ\nabla_{\!X}{\!f}$ ($k$ times), and let $\nabla_{\!X}^{0}{\!f}=f$.

The following special cases of Remark 3.2 will be particularly useful. When $n=2$, Eq. (1) yields

The convergence proofs require a characterization of the second order terms of $f$ near a critical point. Consider the second covariant differential $\nabla\nabla{\!f}=\nabla^{2}{\!f}$ of a smooth function $f\colon M\to{\bf R}$. If $(U,x^{1},\ldots,x^{n})$ is a coordinate chart on $M$, then at $p\in U$ this $(0,2)$ tensor takes the form

Let $M$ be a complete Riemannian manifold with Riemannian structure $g$ and Levi-Civita connection $\nabla$. Let $f\in C^{\infty}(M)$ have a nondegenerate critical point at ${\hat{p}}$ such that the Hessian $(d^{2}{\!f})_{\hat{p}}$ is positive definite. Let $p_{i}$ be a sequence of points in $M$ converging to ${\hat{p}}$ and $H_{i}\in T_{p_{i}}$ a sequence of tangent vectors such that

where $\lambda_{i}$ is chosen such that $f(\exp\lambda_{i}H_{i})\leq f(\exp\lambda H_{i})$ for all $\lambda\geq 0$. Then there exists a constant $E$ and a $\theta\in[0,1)$ such that for all $i=0$, $1$, …,

Proof.The proof is a generalization of the one given in Polak [37, p. 242ff] for the method of steepest descent on Euclidean space.

The existence of a convergent sequence is guaranteed by the smoothness of $f$. If $p_{j}={\hat{p}}$ for some integer $j$, the assertion becomes trivial; assume otherwise. By the smoothness of $f\!$, there exists an open neighborhood $U$ of ${\hat{p}}$ such that $(\nabla^{2}{\!f})_{p}$ is positive definite for all $p\in U$. Therefore, there exist constants $k>0$ and $K\geq k>0$ such that for all $X\in T_{p}$ and all $p\in U$,

Define $X_{i}\in T_{{\hat{p}}}$ by the relations $\exp X_{i}=p_{i}$, $i=0$, $1$, … By assumption, $df_{\hat{p}}=0$ and from Eq. (2), we have

Combining this equality with the inequalities of (5) yields

Next, use (6) with Schwarz’s inequality and the first inequality of (7) to obtain

Define the function ${\mit\Delta}\colon T_{p}\times{\bf R}\to{\bf R}$ by the equation ${\mit\Delta}(X,\lambda)=f(\exp_{p}\lambda X)-f(p)$. By Eq. (2), the second order Taylor formula, we have

Using assumption (ii) of the theorem along with (5) we establish for $\lambda\geq 0$

We may now compute an upper bound for the rate of linear convergence $\theta$. By assumption (i) of the theorem, $\lambda$ must be chosen to minimize the right hand side of (9). This corresponds to choosing \lambda=c\|(\mathop{\rm grad}\nolimits{\!f})_{p_{i}}\|\big{/}K\|H_{i}\|. A computation reveals that

Applying (7) and (8) to this inequality and rearranging terms yields

where \theta=\bigl{(}1-(ck/K)^{2}\bigr{)}. By assumption, $c\in(0,1]$ and $0<k\leq K$, therefore $\theta\in[0,1)$. (Note that Schwarz’s inequality bounds $c$ below unity.) From (10) it is seen that \bigl{(}f(p_{i})-f({\hat{p}})\bigr{)}\leq E\theta^{i} where E=\bigl{(}f(p_{0})-f({\hat{p}})\bigr{)}. From (7) we conclude that for $i=0$, $1$, …,

If Algorithm 3.1 converges to a local minimum, it converges linearly.

The choice $H_{i}=-(\mathop{\rm grad}\nolimits{\!f})_{p_{i}}$ yields $c=1$ in the second assumption the Theorem 3.3, which establishes the corollary.

Let $S^{n-1}$ be the imbedded sphere in ${\bf R}^{n}$, i.e., $S^{n-1}=\{\,x\in{\bf R}^{n}:x^{\scriptscriptstyle\rm T}x=1\,\}$, where $x^{\scriptscriptstyle\rm T}y$ denotes the standard inner product on ${\bf R}^{n}$, which induces a metric on $S^{n-1}$. Geodesics on the sphere are great circles and parallel translation along geodesics is equivalent to rotating the tangent plane along the great circle. Let $x\in S^{n-1}$ and $h\in T_{x}$ have unit length, and $v\in T_{x}$ be any tangent vector. Then

where $\tau$ is the parallelism along the geodesic $t\mapsto\exp th$. Let $Q$ be an $n$-by-$n$ positive definite symmetric matrix with distinct eigenvalues and define $\rho\colon S^{n-1}\to{\bf R}$ by $\rho(x)=x^{\scriptscriptstyle\rm T}Qx$. A computation shows that

The function $\rho$ has a unique minimum and maximum point at the eigenvectors corresponding to the smallest and largest eigenvalues of $Q$, respectively. Because $S^{n-1}$ is geodesically complete, the method of steepest descent in the opposite direction of the gradient converges to the eigenvector corresponding to the smallest eigenvalue of $Q$; likewise for the eigenvector corresponding to the largest eigenvalue. Chu considers the continuous limit of this problem. A computation shows that $\rho(x)$ is maximized along the geodesic $\exp_{x}th$ ($\|h\|=1$) when $a\cos 2t-b\sin 2t=0$, where $a=2x^{\scriptscriptstyle\rm T}Qh$ and $b=\rho(x)-\rho(h)$. Thus $\cos t$ and $\sin t$ may be computed with simple algebraic functions of $a$ and $b$ (which appear below in Algorithm 5.5). The results of a numerical experiment demonstrating the convergence of the method of steepest descent applied to maximizing Rayleigh’s quotient on $S^{20}$ are shown in Figure 1 on page 1.

Consider the function $f(\Theta)=\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$ on the special orthogonal group $\mathop{\it SO}\nolimits(n)$, where $Q$ is a real symmetric matrix with distinct eigenvalues and $N$ is a real diagonal matrix with distinct diagonal elements. It will be convenient to identify tangent vectors in $T_{\Theta}$ with tangent vectors in $T_{I}\cong\mathop{so}\nolimits(n)$, the tangent plane at the identity, via left translation. The gradient of $f$ (with respect to the negative Killing form of $\mathop{so}\nolimits(n)$, scaled by $1/(n-2)$) at $\Theta\in\mathop{\it SO}\nolimits(n)$ is $\Theta[H,N]$, where $H=\mathop{\rm Ad}\nolimits_{\Theta^{\scriptscriptstyle\rm T}}(Q)=\Theta^{\scriptscriptstyle\rm T}Q\Theta$. The group $\mathop{\it SO}\nolimits(n)$ acts on the set of symmetric matrices by conjugation; the orbit of $Q$ under the action of $\mathop{\it SO}\nolimits(n)$ is an isospectral submanifold of the symmetric matrices. We seek a ${\hat{\Theta}}$ such that $f({\hat{\Theta}})$ is maximized. This point corresponds to a diagonal matrix whose diagonal entries are ordered similarly to those of $N$. A related example is found in Smith , who considers the homogeneous space of matrices with fixed singular values, and in Chu .

The Levi-Civita connection on $\mathop{\it SO}\nolimits(n)$ is bi-invariant and invariant with respect to inversion; therefore, geodesics and parallel translation may be computed via matrix exponentiation of elements in $\mathop{so}\nolimits(n)$ and left (or right) translation [25, Ch. II, Ex. 6]. The geodesic emanating from the identity in $\mathop{\it SO}\nolimits(n)$ in direction $X\in\mathop{so}\nolimits(n)$ is given by the formula $\exp_{I}tX=e^{tX}$, where the right hand side denotes regular matrix exponentiation. The expense of geodesic minimization may be avoided if instead one uses Brockett’s estimate for the step size. Given $\Omega\in\mathop{so}\nolimits(n)$, we wish to find $t>0$ such that $\phi(t)=\mathop{\rm tr}\nolimits\mathop{\rm Ad}\nolimits_{e^{-t\Omega}}(H)N$ is minimized. Differentiating $\phi$ twice shows that $\phi^{\prime}(t)=-\mathop{\rm tr}\nolimits\mathop{\rm Ad}\nolimits_{e^{-t\Omega}}(\mathop{\rm ad}\nolimits_{\Omega}H)N$ and $\phi^{\prime\prime}(t)=-\mathop{\rm tr}\nolimits\mathop{\rm Ad}\nolimits_{e^{-t\Omega}}(\mathop{\rm ad}\nolimits_{\Omega}H)\mathop{\rm ad}\nolimits_{\Omega}N$, where $\mathop{\rm ad}\nolimits_{\Omega}A=[\Omega,A]$. Hence, $\phi^{\prime}(0)=2\mathop{\rm tr}\nolimits H\Omega N$ and, by Schwarz’s inequality and the fact that $\mathop{\rm Ad}\nolimits$ is an isometry, $|\phi^{\prime\prime}(t)|\leq\|\mathop{\rm ad}\nolimits_{\Omega}H\|\;\|\mathop{\rm ad}\nolimits_{\Omega}N\|$. We conclude that if $\phi^{\prime}(0)>0$, then $\phi^{\prime}$ is nonnegative on the interval

which provides an estimate for the step size of Step 1 in Algorithm 3.1. The results of a numerical experiment demonstrating the convergence of the method of steepest descent (ascent) in $\mathop{\it SO}\nolimits(20)$ using this estimate are shown in Figure 2.

Newton’s method on Riemannian manifolds

As in the optimization of functions on Euclidean space, quadratic convergence can be obtained if the second order terms of the Taylor expansion are used appropriately. In this section we present Newton’s algorithm on Riemannian manifolds, prove that its convergence is quadratic, and provide examples. Whereas the convergence proof for the method of steepest descent relies upon the Taylor expansion of the function $f\!$, the convergence proof for Newton’s method will rely upon the Taylor expansion of the one-form $df$. Note that Newton’s method has a counterpart in the theory of constrained optimization, as described by, e.g., Fletcher , Bertsekas , or Dunn . The Newton method presented in this section has only local convergence properties. There is a theory of global Newton methods on Euclidean space and computational complexity; see the work of Hirsch and Smale , Smale , and Shub and Smale .

Let $M$ be an $n$-dimensional Riemannian manifold with Riemannian structure $g$ and Levi-Civita connection $\nabla$, let $\mu$ be a $C^{\infty}$ one-form on $M$, and let $p$ in $M$ be such that the bilinear form $({\nabla\!\mu})_{p}\colon T_{p}\times T_{p}\to{\bf R}$ is nondegenerate. Then, by abuse of notation, we have the pair of isomorphisms

with the forward map defined by $X\mapsto(\nabla_{\!X}\mu)_{p}=({\nabla\!\mu})_{p}(\mathchar 513\relax,X)$, which is nonsingular. The notation $({\nabla\!\mu})_{p}$ will henceforth be used for both the bilinear form defined by the covariant differential of $\mu$ evaluated at $p$ and the homomorphism from $T_{p}$ to $T_{p}^{*}$ induced by this bilinear form. In case of an isomorphism, the inverse can be used to compute a point in $M$ where $\mu$ vanishes, if such a point exists. The case $\mu=df$ will be of particular interest, in which case ${\nabla\!\mu}=\nabla^{2}{\!f}$. Before expounding on these ideas, we make the following remarks.

This remark can be generalized in the following way.

The remark follows by applying Remark 4.1 and the Taylor’s theorem of real analysis to the function $\lambda\mapsto(\tau_{\lambda}^{-1}\mu_{p_{\lambda}})(A)$ for any $A$ in $T_{p}$.

Remarks 4.1 and 4.2 can be generalized to $C^{\infty}$ tensor fields, but we will only require Remark 4.2 for case $n=2$ to make the following observation.

Let $\mu$ be a one-form on $M$ such that for some ${\hat{p}}$ in $M$, $\mu_{\hat{p}}=0$. Given any $p$ in a normal neighborhood of ${\hat{p}}$, we wish to find $X$ in $T_{p}$ such that $\exp_{p}X={\hat{p}}$. Consider the Taylor expansion of $\mu$ about $p$, and let $\tau$ be the parallel translation along the unique geodesic joining $p$ to ${\hat{p}}$. We have by our assumption that $\mu$ vanishes at ${\hat{p}}$, and from Eq. (14) for $n=2$,

If the bilinear form $({\nabla\!\mu})_{p}$ is nondegenerate, the tangent vector $X$ may be approximated by discarding the higher order terms and solving the resulting linear equation

This approximation is the basis of the following algorithm.

Let $M$ be a complete Riemannian manifold with Riemannian structure $g$ and Levi-Civita connection $\nabla$, and let $\mu$ be a $C^{\infty}$ one-form on $M$.

Select $p_{0}\in M$ such that $({\nabla\!\mu})_{p_{0}}$ is nondegenerate, and set $i=0$.

(assume that $({\nabla\!\mu})_{p_{i}}$ is nondegenerate), increment $i$, and repeat.

It can be shown that if $p_{0}$ is chosen suitably close (within the so-called domain of attraction) to a point ${\hat{p}}$ in $M$ such that $\mu_{\hat{p}}=0$ and $({\nabla\!\mu})_{\hat{p}}$ is nondegenerate, then Algorithm 4.3 converges quadratically to ${\hat{p}}$. The following theorem holds for general one-forms; we will consider the case where $\mu$ is exact.

Let $f\in C^{\infty}(M)$ have a nondegenerate critical point at ${\hat{p}}$. Then there exists a neighborhood $U$ of ${\hat{p}}$ such that for any $p_{0}\in U$, the iterates of Algorithm 4.3 for $\mu=df$ are well defined and converge quadratically to ${\hat{p}}$.

The proof of this theorem is a generalization of the corresponding proof for Euclidean spaces, with an extra term containing the Riemannian curvature tensor (which of course vanishes in the latter case).

Proof.If $p_{j}={\hat{p}}$ for some integer $j$, the assertion becomes trivial; assume otherwise. Define $X_{i}\in T_{p_{i}}$ by the relations ${\hat{p}}=\exp X_{i}$, $i=0$, $1$, …, so that $d(p_{i},{\hat{p}})=\|X_{i}\|$ (n.b. this convention is opposite that used in the proof of Theorem 3.3). Consider the geodesic triangle with vertices $p_{i}$, $p_{i+1}$, and ${\hat{p}}$, and sides $\exp tX_{i}$ from $p_{i}$ to ${\hat{p}}$, $\exp tH_{i}$ from $p_{i}$ to $p_{i+1}$, and $\exp tX_{i+1}$ from $p_{i+1}$ to ${\hat{p}}$, for $t\in$. Let $\tau$ be the parallelism with respect to the side $\exp tH_{i}$ between $p_{i}$ and $p_{i+1}$. There exists a unique tangent vector ${\mit\Xi}_{i}$ in $T_{p_{i}}$ defined by the equation

(${\mit\Xi}_{i}$ may be interpreted as the amount by which vector addition fails). If we use the definition $H_{i}=-(\nabla^{2}{\!f})_{p_{i}}^{-1}df_{p_{i}}$ of Algorithm 4.3, apply the isomorphism $(\nabla^{2}{\!f})_{p_{i}}\colon T_{p_{i}}\to T_{p_{i}}^{*}$ to both sides of Eq. (15), we obtain the equation

By Taylor’s theorem, there exists an $\alpha\in$ such that

By the smoothness of $f$ and $g$, there exists an $\epsilon>0$ and constants $\delta^{\prime}$, $\delta^{\prime\prime}$, $\delta^{\prime\prime\prime}$, all greater than zero, such that whenever $p$ is in the convex normal ball ${B_{\epsilon}({\hat{p}})}$,

where the induced norm on $T_{p}^{*}$ is used in all three cases. Taking the norm of both sides of Eq. (18), applying the triangle inequality to the right hand side, and using the fact that parallel translation is an isometry, we obtain the inequality

The length of ${\mit\Xi}_{i}$ can be bounded by a cubic expression in $d(p_{i},{\hat{p}})$ by considering the distance between the points $\exp(H_{i}+\tau^{-1}X_{i+1})$ and $\exp X_{i+1}={\hat{p}}$. Given $p\in M$, $\epsilon>0$ small enough, let $a$, $v\in T_{p}$ be such that $\|a\|+\|v\|\leq\epsilon$, and let $\tau$ be the parallel translation with respect to the geodesic from $p$ to $q=\exp_{p}a$. Karcher [29, App. C2.2] shows that

where $K$ is the sectional curvature of $M$ along any section in the tangent plane at any point near $p$.

If $(\nabla^{2}{\!f})_{\hat{p}}$ is positive (negative) definite and Algorithm 4.3 converges to ${\hat{p}}$, then Algorithm 4.3 converges quadratically to a local minimum (maximum) of $f$.

Let $S^{n-1}$ and $\rho(x)=x^{\scriptscriptstyle\rm T}Qx$ be as in Example 3.5. It will be convenient to work with the coordinates $x^{1}$, …, $x^{n}$ of the ambient space ${\bf R}^{n}$, treat the tangent plane $T_{x}S^{n-1}$ as a vector subspace of ${\bf R}^{n}$, and make the identification $T_{x}S^{n-1}\cong T_{x}^{*}S^{n-1}$ via the metric. In this coordinate system, geodesics on the sphere obey the second order differential equation $\ddot{x}^{k}+x^{k}=0$, $k=1$, …, $n$. Thus the Christoffel symbols are given by $\Gamma_{ij}^{k}=\delta_{ij}x^{k}$, where $\delta_{ij}$ is the Kronecker delta. The $ij$th component of the second covariant differential of $\rho$ at $x$ in $S^{n-1}$ is given by (cf. Eq. (4))

Let $u$ be a tangent vector in $T_{x}S^{n-1}$. A linear operator $A\colon{\bf R}^{n}\to{\bf R}^{n}$ defines a linear operator on the tangent plane $T_{x}S^{n-1}$ for each $x$ in $S^{n-1}$ such that

If $A$ is invertible as an endomorphism of the ambient space ${\bf R}^{n}$, the solution to the linear equation $A\mathchar 513\relax u=v$ for $u$, $v$ in $T_{x}S^{n-1}$ is

For Newton’s method, the direction $H_{i}$ in $T_{x}S^{n-1}$ is the solution of the equation

Combining Eqs. (12), (21), and (22), we obtain

where \alpha_{i}=1\big{/}x_{i}^{\scriptscriptstyle\rm T}(Q-\rho(x_{i})I)^{-1}x_{i}. This gives rise to the following algorithm for computing eigenvectors of the symmetric matrix $Q$.

Let $Q$ be a real symmetric $n$-by-$n$ matrix.

Select $x_{0}$ in ${\bf R}^{n}$ such that $x_{0}^{\scriptscriptstyle\rm T}x_{0}=1$, and set $i=0$.

and set \alpha_{i}=1\big{/}x_{i}^{\scriptscriptstyle\rm T}y_{i}.

The quadratic convergence guaranteed by Theorem 4.4 is in fact too conservative for Algorithm 4.7. As evidenced by Figure 1, Algorithm 4.7 converges cubically.

If $\lambda$ is a distinct eigenvalue of the symmetric matrix $Q$, and Algorithm 4.7 converges to the corresponding eigenvector ${\hat{x}}$, then it converges cubically.

Proof 1.In the coordinates $x^{1}$, …, $x^{n}$ of the ambient space ${\bf R}^{n}$, the $ijk$th component of the third covariant differential of $\rho$ at ${\hat{x}}$ is $-2\lambda{\hat{x}}^{k}\delta_{ij}$. Let $X\in T_{\hat{x}}S^{n-1}$. Then $(\nabla^{3}\rho)_{\hat{x}}(\mathchar 513\relax,X,X)=0$ and the second order terms on the right hand side of Eq. (18) vanish at the critical point. The proposition follows from the smoothness of $\rho$.

Proof 2.The proof follows Parlett’s [35, p. 72ff] proof of cubic convergence for the Rayleigh quotient iteration. Assume that for all $i$, $x_{i}\neq{\hat{x}}$, and denote $\rho(x_{i})$ by $\rho_{i}$. For all $i$, there is an angle $\psi_{i}$ and a unit length vector $u_{i}$ defined by the equation $x_{i}={\hat{x}}\cos\psi_{i}+u_{i}\sin\psi_{i}$, such that ${\hat{x}}^{\scriptscriptstyle\rm T}u_{i}=0$. By Algorithm 4.7

where $\beta_{i}=\cos\theta_{i}-\sin\theta_{i}/\theta_{i}$. Therefore,

The following equalities and low order approximations in terms of the small quantities $\lambda-\rho_{i}$, $\theta_{i}$, and $\psi_{i}$ are straightforward to establish: ${\lambda-\rho_{i}}={(\lambda-\rho(u_{i}))}\sin^{2}\psi_{i}$, $\theta_{i}^{2}=\cos^{2}\psi_{i}\sin^{2}\psi_{i}+{\rm h.o.t.}$, $\alpha_{i}={(\lambda-\rho_{i})}+{\rm h.o.t.}$, and $\beta_{i}=-\theta_{i}^{2}/3+{\rm h.o.t.}$ Thus, the denominator of the large fraction in Eq. (23) is of order unity and the numerator is of order $\sin^{2}\psi_{i}$. Therefore, we have

If Algorithm 4.7 is simplified by replacing Step 2 with

then we obtain the Rayleigh quotient iteration. These two algorithms differ by the method in which they use the vector $y_{i}=(Q-\rho(x_{i})I)^{-1}x_{i}$ to compute the next iterate on the sphere. Algorithm 4.7 computes the point $H_{i}$ in $T_{x_{i}}S^{n-1}$ where $y_{i}$ intersects this tangent plane, then computes $x_{i+1}$ via the exponential map of this vector (which “rolls” the tangent vector $H_{i}$ onto the sphere). The Rayleigh quotient iteration computes the intersection of $y_{i}$ with the sphere itself and takes this intersection to be $x_{i+1}$. The latter approach approximates Algorithm 4.7 up to quadratic terms when $x_{i}$ is close to an eigenvector. Algorithm 4.7 is more expensive to compute than—though of the same order as—the Rayleigh quotient iteration; thus, the RQI is seen to be an efficient approximation of Newton’s method.

If the exponential map is replaced by the chart $v\in T_{x}\mapsto(x+v)/\|x+v\|\in S^{n-1}$, Shub shows that a corresponding version of Newton’s method is equivalent to the RQI.

Let $\Theta$, $Q$, $H=\mathop{\rm Ad}\nolimits_{\Theta^{\scriptscriptstyle\rm T}}(Q)$, and $\Omega$ be as in Example 3.6. The second covariant differential of $f(\Theta)=\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$ may be computed either by polarization of the second order term of $\mathop{\rm tr}\nolimits\mathop{\rm Ad}\nolimits_{e^{-t\Omega}}(H)N$, or by covariant differentiation of the differential $df_{\Theta}=-\mathop{\rm tr}\nolimits[H,N]\Theta^{\scriptscriptstyle\rm T}(\mathchar 513\relax)$:

where $X$, $Y\in\mathop{so}\nolimits(n)$. To compute the direction $\Theta X\in T_{\Theta}$, $X\in\mathop{so}\nolimits(n)$, for Newton’s method, we must solve the equation $(\nabla^{2}{\!f})_{\Theta}(\Theta\mathchar 513\relax,\Theta X)=df_{\Theta}$, which yields the linear equation

The linear operator $L_{\Theta}\colon\mathop{so}\nolimits(n)\to\mathop{so}\nolimits(n)$ is self-adjoint for all $\Theta$ and, in a neighborhood of the maximum, negative definite. Therefore, standard iterative techniques in the vector space $\mathop{so}\nolimits(n)$, such as the classical conjugate gradient method, may be used to solve this equation near the maximum. The results of a numerical experiment demonstrating the convergence of Newton’s method in $\mathop{\it SO}\nolimits(20)$ are shown in Figure 2. As can be seen, Newton’s method converged within round-off error in two iterations.

If Newton’s method applied to the function $f(\Theta)=\mathop{\rm tr}\nolimits\Theta^{\scriptscriptstyle\rm T}Q\Theta N$ converges to the point ${\hat{\Theta}}$ such that $\mathop{\rm Ad}\nolimits_{{\hat{\Theta}}^{\scriptscriptstyle\rm T}}(Q)=H_{\infty}=\alpha N$, $\alpha\in{\bf R}$, then it converges cubically.

Proof.By covariant differentiation of $\nabla^{2}{\!f}\!$, the third covariant differential of $f$ at $\Theta$ evaluated at the tangent vectors $\Theta X$, $\Theta Y$, $\Theta Z\in T_{\Theta}$, $X$, $Y$, $Z\in\mathop{so}\nolimits(n)$, is

If $H=\alpha N$, $\alpha\in{\bf R}$, then $(\nabla^{3}{\!f})_{\Theta}(\mathchar 513\relax,\Theta X,\Theta X)=0$. Therefore, the second order terms on the right hand side of Eq. (18) vanish at the critical point. The remark follows from the smoothness of $f$.

This remark illuminates how rapid convergence of Newton’s method applied to the function $f$ can be achieved in some instances. If $E_{ij}\in\mathop{so}\nolimits(n)$ is a matrix with entry $+1$ at element $(i,j)$, $-1$ at element $(j,i)$, and zero elsewhere, $X=\sum_{i<j}x^{ij}E_{ij}$, $H=\mathop{\rm diag}\nolimits(h_{1},\ldots,h_{n})$, and $N=\mathop{\rm diag}\nolimits(\nu_{1},\ldots,\nu_{n})$, then

If the $h_{i}$ are close to $\alpha\nu_{i}$, $\alpha\in{\bf R}$, for all $i$, then $(\nabla^{3}{\!f})_{\Theta}(\mathchar 513\relax,\Theta X,\Theta X)$ may be small, yielding a fast rate of quadratic convergence.

Let $\pi$ be the projection of a square matrix onto its diagonal, and let $Q$ be as above. Consider the maximization of the function $f(\Theta)=\mathop{\rm tr}\nolimits H\pi(H)$, $H=\mathop{\rm Ad}\nolimits_{\Theta^{\scriptscriptstyle\rm T}}(Q)$, on the special orthogonal group. This is equivalent to minimizing the sum of the squares of the off-diagonal elements of $H$ (Golub and Van Loan derive the classical Jacobi method). The gradient of this function at $\Theta$ is $2\Theta[H,\pi(H)]$ . By repeated covariant differentiation of $f\!$, we find

where $I$ is the identity matrix and $X$, $Y$, $Z\in\mathop{so}\nolimits(n)$. It is easily shown that if $[H,\pi(H)]=0$, i.e., if $H$ is diagonal, then $(\nabla^{3}{\!f})_{\Theta}(\mathchar 513\relax,\Theta X,\Theta X)=0$ (n.b. $\pi(\mathop{\rm ad}\nolimits_{X}H)=0$). Therefore, by the same argument as the proof of Remark 4.11, Newton’s method applied to the function $\mathop{\rm tr}\nolimits H\pi(H)$ converges cubically.

Conjugate gradient on Riemannian manifolds

The method of steepest descent provides an optimization technique which is relatively inexpensive per iteration, but converges relatively slowly. Each step requires the computation of a geodesic and a gradient direction. Newton’s method provides a technique which is more costly both in terms of computational complexity and memory requirements, but converges relatively rapidly. Each step requires the computation of a geodesic, a gradient, a second covariant differential, and its inverse. In this section we describe the conjugate gradient method, which has the dual advantages of algorithmic simplicity and superlinear convergence.

Hestenes and Stiefel first used conjugate gradient methods to compute the solutions of linear equations, or, equivalently, to compute the minimum of a quadratic form on ${\bf R}^{n}$. This approach can be modified to yield effective algorithms to compute the minima of nonquadratic functions on ${\bf R}^{n}$. In particular, Fletcher and Reeves and Polak and Ribière provide algorithms based upon the assumption that the second order Taylor expansion of the function to be minimized sufficiently approximates this function near the minimum. In addition, Davidon, Fletcher, and Reeves developed the variable metric methods , but these will not be discussed here. One noteworthy feature of conjugate gradient algorithms on ${\bf R}^{n}$ is that when the function in question is quadratic, they compute its minimum in no more than $n$ steps.

The conjugate gradient method on Euclidean space is uncomplicated. Given a function $f\colon{\bf R}^{n}\to{\bf R}$ with continuous second derivatives and a local minimum at ${\hat{x}}$, and an initial point $x_{0}\in{\bf R}^{n}$, the algorithm is initialized by computing the (negative) gradient direction $G_{0}=H_{0}=-(\mathop{\rm grad}\nolimits{\!f})_{x_{0}}$. The recursive part of the algorithm involves (i) a line minimization of $f$ along the affine space $x_{i}+tH_{i}$, $t\in{\bf R}$, where the minimum occurs at, say, $t=\lambda_{i}$, (ii) computation of the step $x_{i+1}=x_{i}+\lambda_{i}H_{i}$, (iii) computation of the (negative) gradient $G_{i+1}=-(\mathop{\rm grad}\nolimits{\!f})_{x_{i+1}}$, and (iv) computation of the next direction for line minimization,

where $\gamma_{i}$ is chosen such that $H_{i}$ and $H_{i+1}$ conjugate with respect to the Hessian matrix of $f$ at ${\hat{x}}$. When $f$ is a quadratic form represented by the symmetric positive definite matrix $Q$, the conjugacy condition becomes $H_{i}^{\scriptscriptstyle\rm T}QH_{i+1}=0$; therefore, $\gamma_{i}=-H_{i}^{\scriptscriptstyle\rm T}QG_{i+1}/H_{i}^{\scriptscriptstyle\rm T}QH_{i}$. It can be shown in this case that the sequence of vectors $G_{i}$ are all mutually orthogonal and the sequence of vectors $H_{i}$ are all mutually conjugate with respect to $Q$. Using these facts, the computation of $\gamma_{i}$ may be simplified with the observation that $\gamma_{i}=\|G_{i+1}\|^{2}/\|G_{i}\|^{2}$ (Fletcher-Reeves) or $\gamma_{i}=(G_{i+1}-G_{i})^{\scriptscriptstyle\rm T}G_{i+1}/\|G_{i}\|^{2}$ (Polak-Ribière). When $f$ is not quadratic, it is assumed that its second order Taylor expansion sufficiently approximates $f$ in a neighborhood of the minimum, and the $\gamma_{i}$ are chosen so that $H_{i}$ and $H_{i+1}$ are conjugate with respect to the matrix $(\partial^{2}{\!f}/\partial x^{i}\partial x^{j})(x_{i+1})$ of second partial derivatives of $f$ at $x_{i+1}$. It may be desirable to “reset” the algorithm by setting $H_{i+1}=G_{i+1}$ every $r$th step (frequently, $r=n$) because the conjugate gradient method does not, in general, converge in $n$ steps if the function $f$ is nonquadratic. However, if $f$ is closely approximated by a quadratic function, the reset strategy may be expected to converge rapidly, whereas the unmodified algorithm may not be.

Many of these ideas have straightforward generalizations in the geometry of Riemannian manifolds; several of them have already appeared. We need only make the following definition.

Given a tensor field $\omega$ of type $(0,2)$ on $M$ such that for $p$ in $M$, $\omega_{p}\colon T_{p}\times T_{p}\to{\bf R}$ is a symmetric bilinear form, the tangent vectors $X$ and $Y$ in $T_{p}$ are said to be $\omega_{p}$-conjugate or conjugate with respect to $\omega_{p}$ if $\omega_{p}(X,Y)=0$.

An outline of the conjugate gradient method on Riemannian manifolds may now be given. Let $M$ be an $n$-dimensional Riemannian manifold with Riemannian structure $g$ and Levi-Civita connection $\nabla$, and let $f\in C^{\infty}(M)$ have a local minimum at ${\hat{p}}$. As in the conjugate gradient method on Euclidean space, choose an initial point $p_{0}$ in $M$ and compute the (negative) gradient directions $G_{0}=H_{0}=-(\mathop{\rm grad}\nolimits{\!f})_{p_{0}}$ in $T_{p_{0}}$. The recursive part of the algorithm involves minimizing $f$ along the geodesic $t\mapsto\exp_{p_{i}}tH_{i}$, $t\in{\bf R}$, making a step along the geodesic to the minimum point $p_{i+1}=\exp\lambda_{i}H_{i}$, computing $G_{i+1}=-(\mathop{\rm grad}\nolimits{\!f})_{p_{i+1}}$, and computing the next direction in $T_{p_{i+1}}$ for geodesic minimization. This direction is given by the formula

where $\tau$ is the parallel translation with respect to the geodesic step from $p_{i}$ to $p_{i+1}$, and $\gamma_{i}$ is chosen such that $\tau H_{i}$ and $H_{i+1}$ are $(\nabla^{2}{\!f})_{p_{i+1}}$-conjugate, i.e.,

Eq. (26) is, in general, expensive to use because the second covariant differential of $f$ appears. However, we can use the Taylor expansion of $df$ about $p_{i+1}$ to compute an efficient approximation of $\gamma_{i}$. By the fact that $p_{i}=\exp_{p_{i+1}}(-\lambda_{i}\tau H_{i})$ and by Eq. (14), we have

Therefore, the numerator of the right hand side of Eq. (26) multiplied by the step size $\lambda_{i}$ can be approximated by the equation