# Geodesic Convexity of the Symmetric Eigenvalue Problem and Convergence of
Riemannian Steepest Descent

Foivos Alimisis Bart Vandereycken  Department of Mathematics, University of
Geneva

###### Abstract

We study the convergence of the Riemannian steepest descent algorithm on the
Grassmann manifold for minimizing the block version of the Rayleigh quotient
of a symmetric and positive semi-definite matrix. Even though this problem is
non-convex in the Euclidean sense and only very locally convex in the
Riemannian sense, we discover a structure for this problem that is similar to
geodesic strong convexity, namely, weak-strong convexity. This allows us to
apply similar arguments from convex optimization when studying the convergence
of the steepest descent algorithm but with initialization conditions that do
not depend on the eigengap $\delta$. When $\delta>0$, we prove exponential
convergence rates, while otherwise the convergence is algebraic. Additionally,
we prove that this problem is geodesically convex in a neighbourhood of the
global minimizer of radius $\mathcal{O}(\sqrt{\delta})$.

##  1 Introduction

We consider the problem of computing the top $k$ eigenvectors of a symmetric
and positive semi-definite matrix $A\in\mathbb{R}^{n\times n}$, which has many
applications in numerical linear algebra (low rank approximation), statistics
(principal component analysis) and signal processing.

We denote by $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ the
eigenvalues of $A$ counted with multiplicity and by
$\delta:=\lambda_{k}-\lambda_{k+1}$ the eigengap for some $k$ between $1$ and
$n-1$. We also denote
$\Lambda_{\alpha}=\textnormal{diag}(\lambda_{1},\ldots,\lambda_{k})$ and
$\Lambda_{\beta}=\textnormal{diag}(\lambda_{k+1},\ldots,\lambda_{n})$.

A set of $k$ leading eigenvectors of $A$ can be found by minimizing the
function

$f(X)=-\operatorname{Tr}(X^{T}AX)$

over the set of $n\times k$ matrices with orthonormal columns. Indeed, from
the Ky-Fan theorem we know that

$\min\\{f(X)\colon X\in\mathbb{R}^{n\times
k},X^{T}X=I_{k}\\}=-(\lambda_{1}+\cdots+\lambda_{k})=-\operatorname{Tr}(\Lambda_{\alpha})=:f^{*}.$
(1)

Since $A$ is symmetric, we can define the matrix
$V_{\alpha}=\begin{bmatrix}v_{1}\hskip 5.69054pt\cdots\hskip
5.69054ptv_{k}\end{bmatrix}$ such that $V_{\alpha}^{T}V_{\alpha}=I_{k}$ and
with $v_{i}\in\mathbb{R}^{n}$ a unit-norm eigenvector corresponding to
$\lambda_{i}$. If the eigengap $\delta$ is strictly positive, then
$\operatorname{span}(V_{\alpha})$ is unique; otherwise, we can choose any
$v_{k}$ from a subspace with dimension equal to the multiplicity of
$\lambda_{k}$. It is readily seen that
$f(V_{\alpha})=-(\lambda_{1}+\cdots+\lambda_{k})$. In fact, all minimizers of
(1) are of the form $V_{\alpha}Q$ with $Q$ a $k\times k$ orthogonal matrix. We
also define $V_{\beta}=\begin{bmatrix}v_{k+1}\hskip 5.69054pt\cdots\hskip
5.69054ptv_{n}\end{bmatrix}$ that contains the eigenvectors corresponding to
the eigenvalues $\lambda_{k+1},\ldots,\lambda_{n}$. Its columns span the
orthogonal complement of $\operatorname{span}(V_{\alpha})$ in $\mathbb{R}^{n}$
and thus $V_{\beta}^{T}V_{\beta}=I_{n-k}$ and
$V_{\alpha}^{T}V_{\beta}={0}_{k\times(n-k)}$.

Since $\operatorname{span}(V_{\alpha})=\operatorname{span}(V_{\alpha}Q)$, it
is more natural to consider this problem as a minimization problem on the
Grassmann manifold $\operatorname{Gr}(n,k)$, the set of $k$-dimensional
subspaces in $\mathbb{R}^{n}$. Let us therefore redefine the objective
function as

$f(\mathcal{X})=-\operatorname{Tr}(X^{T}AX)\text{ where
$\mathcal{X}=\operatorname{span}(X)$ for $X\in\mathbb{R}^{n\times k}$ s.t.
$X^{T}X=I_{k}$}.$ (2)

This cost function can be seen as a block version of the standard Rayleigh
quotient. An immediate benefit is that, if $\delta>0$, the minimizer of (2) is
isolated since it is the subspace
$\mathcal{V_{\alpha}}=\operatorname{span}(V_{\alpha})$.

To minimize $f$ on $\operatorname{Gr}(n,k)$, we shall use the Riemannian
steepest descent method (RSD) along geodesics in $\operatorname{Gr}(n,k)$.
Quite remarkably, for $\operatorname{Gr}(n,k)$ these geodesics can be
implemented efficiently in closed form.

For analyzing the convergence properties of steepest descent on
$\operatorname{Gr}(n,k)$, we extend results of the recent work [4], where it
is shown that the Rayleigh quotient on the sphere enjoys favourable geodesic
convexity-like properties, namely, _weak-quasi-convexity_ and _quadratic
growth_. In this work, we show that these convexity-like properties continue
to hold in the more general case of the block Rayleigh quotient function
$f\colon\operatorname{Gr}(n,k)\rightarrow\mathbb{R}$. These results are of
general interest, but also sufficient to prove a local convergence rate for
steepest descent for minimizing $f$ when started from an initial point outside
the region of local convexity. For the latter, a crucial help is provided by
the fact that the Grassmann manifold is positively curved.

In particular, assuming a _strictly positive eigengap_ $\delta$ between
$\lambda_{k}$ and $\lambda_{k+1}$, we prove an exponential convergence rate to
the subspace spanned by the $k$ leading eigenvectors, similar to the
convergence of power method and subspace iteration (Theorem 11). If we do not
assume any knowledge regarding the eigengap, then we can still prove a sub-
exponential (polynomial) convergence rate of the function values to the global
minimum (Theorem 13), but we cannot directly study the convergence to a global
minimizer. This is in line with previous work but our analysis does not use
standard notions of geodesic convexity and allows for an initial guess further
from the global minimizer. In Appendix B we present related convergence
results for steepest descent with a more tractable step size but at the
expense of needing a slightly better initialization.

##  2 Related work

Over the last few years, different aspects of the convexity of eigenvalue
problems have received quite some attention. In [26], the authors prove (in
Theorem 4) that the Rayleigh quotient is geodesically gradient dominated in
the sphere, that is, it satisfies a spherical version of the
Polyak–Łojasiewicz inequality. In [4], it is shown that this result of [26]
can be strengthened to a geodesic weak-quasi-convexity and quadratic growth
property, which imply gradient dominance when combined. Finally, the recent
paper [3] examines (among other contributions) the convexity structure of the
same block version of the symmetric eigenvalue problem on the Grassmann
manifold that we introduced above. Unfortunately, the characterization of the
geodesic convexity region independently of the eigengap $\delta$ (Corollary 5
in [3]) is wrong (see our Appendix A). As we will prove in Theorem 18, the
geodesic convexity region of $f$ (and the one of the equivalent cost function
used in [3]) needs to depend on the eigengap, as appears also in [16, Lemma 7]
in the case of the sphere ($k=1$).

To the best of our knowledge, the current work is the first that deals with
the convergence of the steepest descent algorithm for the multiple eigenvalue-
eigenvector problem on the Grassmann manifold. The work [4] proves exponential
convergence of steepest descent only in the case of $k=1$, that is, for the
leading eigenvector. In this paper, we take a reasonable but highly non-
trivial step forward by extending this analysis for general $k$, that is, for
a block of $k$ leading eigenvectors.

The standard algorithm for computing the leading eigenspace of dimension $k$
is subspace iteration (or power method when $k=1$).111Krylov methods are
arguably the most popular algorithms but they do not iterate on a subspace
directly and are typically started from a single vector. In particular, they
cannot easily improve a given approximation of a subspace for large $k>1$.
However, there are reasons to believe that, in certain cases, Riemannian
steepest descent (and its accelerated version with non-linear conjugate
gradients) should be preferred, especially in noisy settings [4] or in
electronic structure calculations where the leading eigenspace of many varying
matrices $A$ needs to be computed.222Personal communication by Yousef Saad. In
particular, [4] presents strong experimental evidence that steepest descent is
more robust to perturbations of the matrix-vector products than subspace
iteration close to the optimum. While subspace iteration still behaves better
at the start of the iteration, it asymptotically fails to converge to an
approximation of the leading subspace that is as good as the one estimated by
Riemannian steepest descent. While [4] dealt with a noisy situation due to
calculations in a distributed setting with limited communication, exactly the
same effect can be observed when we inject the matrix-vector products with
Gaussian noise. Thus, we expect steepest descent to perform better than
subspace iteration close to the optimum in any stochastic regime [12].

Regarding worst-case theoretical guarantees, the strongest convergence result
for subspace iteration in the presence of a strictly positive eigengap
$\delta$ is in terms of the largest principal angle between the iterates and
the optimum [11], that is, the $\ell_{\infty}$-norm of the vector of principal
angles. In contrast, our convergence result for steepest descent for
$\delta>0$ (Theorem 11) is in terms of the $\ell_{2}$-norm of the same vector
of angles, which is in general stronger. When $\delta=0$, it is known from
[21, 17] that the largest eigenvalue ($k=1$) can still be efficiently
estimated. We extend this result for $k>1$ and prove a convergence rate of
steepest descent for the function values $f$ (Theorem 13), relying only on
weak-quasi-convexity (and thus using a different argument from [21, 17]).

##  3 Geometry of the Grassmann manifold and block Rayleigh quotient

We present here a brief introduction into the geometry of the Grassmann
manifold. The content is not new and for more details, we refer to [2, 7, 10].

The $(n,k)$-Grassmann manifold is defined as the set of all $k$-dimensional
subspaces of $\mathbb{R}^{n}$:

$\operatorname{Gr}(n,k)=\\{\mathcal{X}\subseteq\mathbb{R}^{n}\colon\mathcal{X}\hskip
2.84526pt\text{is a subspace and}\dim(\mathcal{X})=k\\}.$

Any element $\mathcal{X}$ of $\operatorname{Gr}(n,k)$ can be represented by a
matrix $X\in\mathbb{R}^{n\times k}$ that satisfies
$\mathcal{X}=\operatorname{span}(X)$. Such a representative is not unique
since $Y=XQ$ for some invertible matrix $Q\in\mathbb{R}^{k\times k}$ satisfies
$\operatorname{span}(Y)=\operatorname{span}(X)$. Without loss of generality,
_we will therefore always take matrix representatives $X$ of subspaces
$\mathcal{X}$ that have orthonormal columns throughout the paper_. With some
care, the non-uniqueness of the representatives is not a problem.333This can
be made very precise by describing $\operatorname{Gr}(n,k)$ as the quotient of
the Stiefel manifold with the orthogonal group. The elegant theory of this
quotient manifold is worked out in [2]. For example, our objective function
(2) is invariant to $Q$.

#### Riemannian structure.

The set $\operatorname{Gr}(n,k)$ admits the structure of a differential
manifold with tangent spaces

$T_{\mathcal{X}}\operatorname{Gr}(n,k)=\\{G\in\mathbb{R}^{n\times k}\colon
X^{T}G=0\\},$ (3)

where $\mathcal{X}=\operatorname{span}(X)$. Since $X^{T}G=0$ if and only if
$X^{T}(GQ)=0$, for any invertible matrix $Q\in\mathbb{R}^{k\times k}$, this
description of the tangent space does not depend on the representative $X$.
However, a specific tangent vector $G$ will depend on the chosen $X$. With
slight abuse of notation,444Using the quotient manifold theory, one would use
horizontal lifts. the above definition should therefore be interpreted as:
given a fixed $X$, we define tangent vectors $G_{1},G_{2},\ldots$ of
$\operatorname{Gr}(n,k)$ at $\mathcal{X}=\operatorname{span}(X)$.

This subtlety is important, for example, when defining an inner product on
$T_{\mathcal{X}}\operatorname{Gr}(n,k)$:

$\langle G_{1},G_{2}\rangle_{\mathcal{X}}=\operatorname{Tr}(G^{T}_{1}G_{2})\
\text{\ with\ $G_{1},G_{2}\in T_{\mathcal{X}}\operatorname{Gr}(n,k)$ }.$

Here, $G_{1}$ and $G_{2}$ are tangent vectors of the same representative $X$.
Observe that the inner product is invariant to the choice of orthonormal
representative: If $\bar{G}_{1}=G_{1}Q$ and $\bar{G}_{2}=G_{2}Q$ with
orthogonal $Q$, then we have

$\langle\bar{G}_{1},\bar{G}_{2}\rangle_{\mathcal{X}}=\operatorname{Tr}(\bar{G}^{T}_{1}\bar{G}_{2})=\operatorname{Tr}(Q^{T}G_{1}^{T}G_{2}Q)=\operatorname{Tr}(G_{1}^{T}G_{2}QQ^{T})=\operatorname{Tr}(G_{1}^{T}G_{2}).$

It is easy to see that the norm induced by this inner product in any tangent
space is the Frobenius norm, which we will denote throughout the paper as
$\|\cdot\|:=\|\cdot\|_{F}$.

#### Exponential map.

Given the Riemannian structure of $\operatorname{Gr}(n,k)$, we can compute the
exponential map at a point $\mathcal{X}$ as [1, Thm. 3.6]

$\displaystyle\operatorname{Exp}_{\mathcal{X}}:T_{\mathcal{X}}\operatorname{Gr}(n,k)$
$\displaystyle\rightarrow\operatorname{Gr}(n,k)$ (4) $\displaystyle G$
$\displaystyle\mapsto\operatorname{span}(\,XV\cos(\Sigma)+U\sin(\Sigma)\,),$

where $U\Sigma V^{T}$ is the _compact_ SVD of $G$ such that $\Sigma$ and $V$
are square matrices.

The exponential map is invertible in the domain [7, Prop. 5.1]

$\left\\{G\in
T_{\mathcal{X}}\operatorname{Gr}(n,k)\colon\|G\|_{2}<\frac{\pi}{2}\right\\},$
(5)

where $\|G\|_{2}$ is the spectral norm of $G$. The inverse of the exponential
map restricted to this domain is the logarithmic map, denoted by
$\operatorname{Log}$. Given two subspaces
$\mathcal{X},\mathcal{Y}\in\operatorname{Gr}(n,k)$, we have

$\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})=U\operatorname{atan}(\widehat{\Sigma})\,V^{T},$
(6)

where $U\widehat{\Sigma}V^{T}=(I-XX^{T})Y(X^{T}Y)^{-1}$ is again a compact
SVD. This is well-defined if $X^{T}Y$ is invertible, which is guaranteed if
all principal angles between $\mathcal{X}$ and $\mathcal{Y}$ are strictly less
than $\pi/2$ (see below). By taking
$G=\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})$, we see that
$\Sigma=\operatorname{atan}(\widehat{\Sigma})$.

#### Principal angles.

The Riemannian structure of the Grassmann manifold can be conveniently
described by the notion of the principal angles between subspaces. Given two
subspaces $\mathcal{X},\mathcal{Y}\in\operatorname{Gr}(n,k)$, the principal
angles between them are $0\leq\theta_{1}\leq\cdots\leq\theta_{k}\leq\pi/2$
obtained from the SVD

$Y^{T}X=U_{1}\cos\theta\ V_{1}^{T}$ (7)

where $U_{1}\in\mathbb{R}^{k\times k},V_{1}\in\mathbb{R}^{k\times k}$ are
orthogonal and the diagonal matrix
$\cos\theta=\operatorname{diag}(\cos\theta_{1},...,\cos\theta_{k})$.

We can express the Riemannian logarithm using principal angles and the
intrinsic distance induced by the Riemannian inner product discussed above is

$\operatorname{dist}(\mathcal{X},\mathcal{Y)}=\|\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})\|=\|\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})\|=\sqrt{\theta_{1}^{2}+...+\theta_{k}^{2}}=\|\theta\|_{2},$
(8)

where $\theta=(\theta_{1},\ldots,\theta_{k})^{T}$.

If $X\in\mathbb{R}^{n\times k}$ is an arbitrary matrix with orthonormal
columns, then, generically, these columns will not be exactly orthogonal to
the $k$ leading eigenvectors $v_{1},\ldots,v_{k}$ of $A$. Thus, we have with
probability one that the principal angles between $\mathcal{X}$ and the space
of $k$ leading eigenvectors satisfy
$0\leq\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$.

#### Curvature.

We can compute exactly the sectional curvatures in $\operatorname{Gr}(n,k)$,
but for our purposes we only need that they are everywhere non-negative [25,
7]. This means that the geodesics on the Grassmann manifold spread more slowly
than in Euclidean space. This is consequence of the famous Toponogov’s theorem
that we state here in the form of the following technical lemma, which will be
important in our convergence analysis.

######  Lemma 1.

Let $\mathcal{X},\mathcal{Y},\mathcal{Z}\in\operatorname{Gr}(n,k)$, such that

$\max\\{\textnormal{dist}(\mathcal{X},\mathcal{Z}),\textnormal{dist}(\mathcal{Y},\mathcal{Z})\\}<\frac{\pi}{2}.$

Then

$\textnormal{dist}(\mathcal{X},\mathcal{Y})\leq\|\log_{\mathcal{Z}}(\mathcal{X})-\log_{\mathcal{Z}}(\mathcal{Y})\|.$

######  Lemma 2.

(Law of cosines) Let $\mathcal{X},\mathcal{Y},\mathcal{Z}$ as in Lemma 1. Then

$\operatorname{dist}^{2}(\mathcal{X},\mathcal{Y})\leq\operatorname{dist}^{2}(\mathcal{Z},\mathcal{X})+\operatorname{dist}^{2}(\mathcal{Z},\mathcal{Y})-2\langle\operatorname{Log}_{\mathcal{Z}}(\mathcal{X}),\operatorname{Log}_{\mathcal{Z}}(\mathcal{Y})\rangle.$

###### Proof.

Apply Lemma 1 and expand
$\|\log_{\mathcal{Z}}(\mathcal{X})-\log_{\mathcal{Z}}(\mathcal{Y})\|^{2}$. ∎

#### Block Rayleigh quotient.

Our objective function for minimization is the block version of the Rayleigh
quotient:

$f(\mathcal{X})=-\operatorname{Tr}(X^{T}AX)\text{ where
$\mathcal{X}=\operatorname{span}(X)\in\operatorname{Gr}(n,k)$ s.t.
$X^{T}X=I_{k}$}.$

This function has
$\mathcal{V}_{\alpha}=\operatorname{span}(\begin{bmatrix}v_{1}\hskip
5.69054pt\cdots\hskip 5.69054ptv_{k}\end{bmatrix})$ as global minimizer. This
minimizer is unique on $\operatorname{Gr}(n,k)$ if and only if $\delta>0$.

Given any differentiable function
$f\colon\operatorname{Gr}(n,k)\rightarrow\mathbb{R}$, we can define its
Riemannian gradient as the vector field that satisfies

$df(\mathcal{X})(\mathcal{G})=\langle\operatorname{grad}f(\mathcal{X}),\mathcal{G}\rangle_{\mathcal{X}}.$

For a given representative $X$ of $\mathcal{X}$, the Riemannian gradient of
the block Rayleigh quotient satisfies

$\operatorname{grad}f(\mathcal{X})=-2(I-XX^{T})AX.$

Using the notions of the Riemannian gradient and Levi-Civita connection, we
can define also a Riemannian notion of Hessian. For the block Rayleigh
quotient $f$, the Riemannian Hessian $\operatorname{Hess}f$ evaluated as
bilinear form satisfies

$\operatorname{Hess}f(\mathcal{X})[{G},{G}]=2\langle G,GX^{T}AX-AG\rangle,$
(9)

for $G\in T_{\mathcal{X}}\operatorname{Gr}(n,k)$; see [10, §4.4] or [2,
§6.4.2].

##  4 Convexity-like properties of the block Rayleigh quotient

We now prove the new analytic properties of the block Rayleigh quotient
$f(\mathcal{X})=-\operatorname{Tr}(X^{T}AX)$. These are important in their own
right but will also be used later for the convergence of the Riemannian
steepest descent method.

###  4.1 Smoothness

A $C^{2}$ function defined on the Grassmann manifold is called $\gamma$-smooth
if the maximum eigenvalue of its Riemannian Hessian is everywhere upper
bounded by $\gamma$. By the second-order Taylor expansion of $f$, it is easy
to see that such a function then satisfies (see, e.g., [2, Thm. 7.1.2])

$f(\mathcal{X})\leq
f(\mathcal{Y})+\langle\operatorname{grad}f(\mathcal{Y}),\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})\rangle+\frac{\gamma}{2}\operatorname{dist}^{2}(\mathcal{X},\mathcal{Y}),$
(10)

for any $\mathcal{X},\mathcal{Y}\in\operatorname{Gr}(n,k)$ with
$\textnormal{dist}(\mathcal{X},\mathcal{Y})<\frac{\pi}{2}$.

As in the introduction, denote the global minimum of $f$ by $f^{*}$ which is
attained at $V_{\alpha}\in\operatorname{Gr}(n,k)$. The previous inequality
also leads to

$f(\mathcal{X})-f^{*}\geq\frac{1}{2\gamma}\|\operatorname{grad}f(\mathcal{X})\|^{2},$
(11)

for any $\mathcal{X}\in\operatorname{Gr}(n,k)$ with
$\textnormal{dist}(\mathcal{X},\mathcal{V}_{\alpha})<\frac{\pi}{2}$. We
present a proof below although the result is well known.

###### Proof of (11).

Since $f^{*}$ is a global minimum of $f$, we have from (10) that

$f^{*}\leq f(\mathcal{X})\leq
f(\mathcal{Y})+\langle\textnormal{grad}f(\mathcal{Y}),\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})\rangle+\frac{\gamma}{2}\|\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})\|^{2},$

for any $\mathcal{X},\mathcal{Y}\in\operatorname{Gr}(n,k)$ with
$\textnormal{dist}(\mathcal{X},\mathcal{Y})<\frac{\pi}{2}$.

We set
$\mathcal{X}:=\operatorname{Exp}_{\mathcal{Y}}\left(-\frac{1}{\gamma}\textnormal{grad}f(\mathcal{Y})\right)$.
By the mean value theorem in $\operatorname{grad}f$, we have
$\tfrac{1}{\gamma}\|\operatorname{grad}f(\mathcal{Y})\|\leq\textnormal{dist}(\mathcal{Y},\mathcal{V}_{\alpha})$.
Indeed, consider the geodesic connecting $\mathcal{Y}$ and
$\mathcal{V}_{\alpha}$ and apply Lemma 4 from [5] for the vector field
$A=\textnormal{grad}f$. Taking norms of both sides and using that
$\textnormal{grad}f(\mathcal{V}_{\alpha})=0$, the norm of the integral is less
or equal than the integral of the norms. Since parallel transport preserves
the norm and the covariant derivative of the gradient is the Riemannian
Hessian, we obtain the claimed inequality. If we now assume that
$\textnormal{dist}(\mathcal{Y},\mathcal{V}_{\alpha})<\frac{\pi}{2}$, then
$-\tfrac{1}{\gamma}\textnormal{grad}f(\mathcal{Y})$ satisfies the condition of
the domain in (5) and $\operatorname{Log}$ is well defined:
$\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})=-\frac{1}{\gamma}\textnormal{grad}f(\mathcal{Y})$.
In that case, we also have

$\textnormal{dist}(\mathcal{X},\mathcal{Y})=\|\operatorname{Log}_{\mathcal{Y}}(\mathcal{X})\|<\frac{\pi}{2}$
(12)

and the right hand side of the initial inequality becomes

$\displaystyle f^{*}\leq
f(\mathcal{Y})-\frac{1}{\gamma}\|\textnormal{grad}f(\mathcal{Y})\|^{2}+\frac{1}{2\gamma}\|\textnormal{grad}f(\mathcal{Y})\|^{2}=f(\mathcal{Y})-\frac{1}{2\gamma}\|\textnormal{grad}f(\mathcal{Y})\|^{2},$

for any $\mathcal{Y}\in\operatorname{Gr}(n,k)$ with
$\textnormal{dist}(\mathcal{Y},\mathcal{V}_{\alpha})<\frac{\pi}{2}$.
Rearranging the last inequality and substituting $\mathcal{Y}=\mathcal{X}$, we
get the desired result. ∎

We start our analysis by showing that the largest eigenvalue of the Riemannian
Hessian of the block Rayleigh quotient $f$ is indeed upper bounded on the
Grassmann manifold. Thus the properties from above hold for the stated
$\gamma$.

######  Proposition 3 (Smoothness).

The eigenvalues of the Riemannian Hessian of $f$ on $\operatorname{Gr}(n,k)$
are upper bounded by $\gamma:=2(\lambda_{1}-\lambda_{n})$.

###### Proof.

Let $G$ be a tangent vector of $\operatorname{Gr}(n,k)$ at $X$. Then the
Riemannian Hessian satisfies (see (9))

$\tfrac{1}{2}\operatorname{Hess}f(\mathcal{X})[G,G]=\operatorname{Tr}(G^{T}GX^{T}AX)-\operatorname{Tr}(AGG^{T}).$

Since $A,X^{T}AX,GG^{T}$, and $G^{T}G$ are all symmetric and positive definite
matrices, standard trace inequality (see, e.g, [15, Thm. 4.3.53]) gives

$\operatorname{Hess}f(\mathcal{X})[G,G]\leq
2(\lambda_{\max}(X^{T}AX)-\lambda_{\min}(A))\|G\|^{2}.$

Since $X$ has orthonormal columns,
$\lambda_{\max}(X^{T}AX)\leq\lambda_{\max}(A)$; see, e.g., [15, Cor. 4.3.37].
The proof is now complete with the definition of $\lambda_{1}$ and
$\lambda_{n}$. ∎

The result in Prop. 3 is tight: Choosing $X=V_{\alpha}$ and
$G=v_{n}e_{1}^{T}$, it is readily verified that the upper bound is attained.

We will also need the following upper bound of the spectral norm of the
Riemannian gradient, independently of $\mathcal{X}$.

######  Lemma 4.

For all $\mathcal{X}\in\operatorname{Gr}(n,k)$, the Riemanian gradient of $f$
satisfies $\|\operatorname{grad}f(\mathcal{X})\|_{2}\leq\frac{\gamma}{2}.$

###### Proof.

Since $X$ has orthonormal columns, we can complete it to the orthogonal matrix
$Q=\begin{bmatrix}X&X_{\perp}\end{bmatrix}$. Hence,
$\|\operatorname{grad}f(\mathcal{X})\|_{2}=\|2(I-XX^{T})AX\|_{2}=2\|X_{\perp}^{T}AX\|_{2}$.
The result now follows directly from [18, Thm. 2] since $A$ is real symmetric
and the definition of $\gamma=2(\lambda_{1}-\lambda_{n})$. ∎

###  4.2 Weak-quasi-convexity and quadratic growth

We now turn our interest in the convexity properties of the block Rayleight
quotient function. We start by proving a property which is known in the
literature as _quadratic growth_.

######  Proposition 5 (Quadratic growth).

Let $0\leq\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$ be the principal angles
between the subspaces $\mathcal{X}$ and $\mathcal{V}_{\alpha}$. The function
$f$ satisfies $f(\mathcal{X})-f^{*}\geq
c_{Q}\,\delta\,\operatorname{dist}^{2}(\mathcal{X},\mathcal{V_{\alpha}})$
where $c_{Q}=4/\pi^{2}>0.4$.

###### Proof.

The spectral decomposition of
$A=V_{\alpha}\Lambda_{\alpha}V_{\alpha}^{T}+V_{\beta}\Lambda_{\beta}V_{\beta}^{T}$
implies

$X^{T}AX=X^{T}V_{\alpha}\Lambda_{\alpha}V_{\alpha}^{T}X+X^{T}V_{\beta}\Lambda_{\beta}V_{\beta}^{T}X.$
(13)

Since $f(\mathcal{X})=-\operatorname{Tr}(X^{T}AX)$, we have

$\displaystyle f(\mathcal{X})-f^{*}$
$\displaystyle=\operatorname{Tr}(\Lambda_{\alpha})-\operatorname{Tr}(X^{T}V_{\alpha}\Lambda_{\alpha}V_{\alpha}^{T}X)-\operatorname{Tr}(X^{T}V_{\beta}\Lambda_{\beta}V_{\beta}^{T}X)$
$\displaystyle=\operatorname{Tr}(\Lambda_{\alpha})-\operatorname{Tr}(\Lambda_{\alpha}V_{\alpha}^{T}XX^{T}V_{\alpha})-\operatorname{Tr}(\Lambda_{\beta}V_{\beta}^{T}XX^{T}V_{\beta})$
$\displaystyle=\operatorname{Tr}(\Lambda_{\alpha}(I_{k}-V_{\alpha}^{T}XX^{T}V_{\alpha}))-\operatorname{Tr}(\Lambda_{\beta}V_{\beta}^{T}XX^{T}V_{\beta}).$

From the definition (7) of the principal angles between $X$ and $V_{\alpha}$,
we recall that

$V_{\alpha}^{T}X=U_{1}\cos\theta\,V_{1}^{T},$ (14)

where $\cos\theta=\operatorname{diag}(\cos\theta_{1},\ldots,\cos\theta_{k})$
is a diagonal matrix and $U_{1},V_{1}$ are orthogonal matrices. Plugging this
equality in, we get that the $j$th eigenvalue of the matrix
$I_{k}-V_{\alpha}^{T}XX^{T}V_{\alpha}$ is equal to
$1-\cos^{2}\theta_{j}=\sin^{2}\theta_{j}\geq 0$. Thus, by standard trace
inequality for symmetric and positive definite matrices (see, e.g., [15, Thm.
4.3.53]), the first summand above satisfies

$\operatorname{Tr}(\Lambda_{\alpha}(I_{k}-V_{\alpha}^{T}XX^{T}V_{\alpha}))\geq\lambda_{k}\sum_{j=1}^{k}\sin^{2}\theta_{j}.$

The matrix $V_{\beta}^{T}XX^{T}V_{\beta}$ has the same non-zero eigenvalues
with the same multiplicity as the matrix

$X^{T}V_{\beta}V_{\beta}^{T}X=I_{k}-V_{1}\cos^{2}\theta\,V_{1}^{T}=V_{1}\sin^{2}\theta\,V_{1}^{T}$

where we used $V_{\beta}V_{\beta}^{T}=I_{n}-V_{\alpha}V_{\alpha}^{T}$ and the
SVD of $V_{\alpha}^{T}X$. Thus the $j$th eigenvalue of
$V_{\beta}^{T}XX^{T}V_{\beta}$ is $\sin^{2}\theta_{j}\geq 0$. By trace
inequality again, the second summand therefore satisfies

$\operatorname{Tr}(\Lambda_{\beta}V_{\beta}^{T}XX^{T}V_{\beta})\leq\lambda_{k+1}\sum_{j=1}^{k}\sin^{2}\theta_{j}.$

Putting both bounds together, we get

$f(\mathcal{X})-f^{*}\geq(\lambda_{k}-\lambda_{k+1})\sum_{j=1}^{k}\sin^{2}\theta_{j}\geq\delta\sum_{j=1}^{k}\frac{4}{\pi^{2}}\theta_{j}^{2}$

and the proof is complete by the definition (8) of $\operatorname{dist}$. ∎

We say that $f$ is geodesically convex if for all $\mathcal{X}$ and
$\mathcal{Y}$ in a suitable region it holds

$f(\mathcal{X})-f(\mathcal{Y})\leq\langle\textnormal{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})\rangle.$

This generalizes the classical convexity of differentiable functions on
$\mathbb{R}^{n}$ to manifolds by taking the logarithmic map instead of the
difference $\mathcal{X}-\mathcal{Y}$.

In Appendix A, we prove that our objective function $f$ is only geodesically
convex in a small neighbourhood of size $\mathcal{O}(\sqrt{\delta})$ around
the minimizer $\mathcal{V}_{\alpha}$. Fortunately, our key result of this
section shows that $f$ satisfies a much weaker notion of geodesic convexity,
known in the literature as _weak-quasi-convexity_ , that does not depend on
the eigengap $\delta$.

We first need the following lemma which is the general version of the CS
decomposition but applied to our setting of square blocks.

######  Lemma 6.

Let $X,Y\in\mathbb{R}^{n\times k}$ be such that $X^{T}X=Y^{T}Y=I_{k}$ with
$k<n$. Choose $X_{\perp},Y_{\perp}\in\mathbb{R}^{n\times(n-k)}$ such that
$X_{\perp}^{T}X_{\perp}=Y_{\perp}^{T}Y_{\perp}=I_{n-k}$ and
$\operatorname{span}(X_{\perp})=\operatorname{span}(X)^{\perp}$,
$\operatorname{span}(Y_{\perp})=\operatorname{span}(Y)^{\perp}$. Then there
exists $0\leq r,s\leq k$ such that

$\displaystyle Y^{T}X$ $\displaystyle=U_{1}\begin{bmatrix}I_{r}\\\ &C_{s}\\\
&&O_{p\times p}\end{bmatrix}V_{1}^{T},$ $\displaystyle Y^{T}X_{\perp}$
$\displaystyle=U_{1}\begin{bmatrix}O_{r\times m}\\\ &S_{s}\\\
&&I_{p}\end{bmatrix}V_{2}^{T}$ $\displaystyle Y_{\perp}^{T}X$
$\displaystyle=U_{2}\begin{bmatrix}O_{m\times r}\\\ &S_{s}\\\
&&I_{p}\end{bmatrix}V_{1}^{T},$ $\displaystyle Y_{\perp}^{T}X_{\perp}$
$\displaystyle=U_{2}\begin{bmatrix}-I_{m}\\\ &-C_{s}\\\ &&O_{p\times
p}\end{bmatrix}V_{2}^{T}$

with $p=k-r-s$ and $m=n-2k+r$, and we have

  * •

orthogonal matrices $U_{1},V_{1}$ of size $k$ and $U_{2},V_{2}$ of size $n-k$;

  * •

identity matrices $I_{q}$ of size $q$;

  * •

zero matrices $O_{q\times t}$ of size $q\times t$;

  * •

diagonal matrices $C_{s}=\operatorname{diag}(\alpha_{1},\ldots,\alpha_{s})$
and $S_{s}=\operatorname{diag}(\beta_{1},\ldots,\beta_{s})$ such that
$1>\alpha_{1}\geq\cdots\geq\alpha_{s}>0$,
$0<\beta_{1}\leq\cdots\leq\beta_{s}<1$ and $C_{s}^{2}+S_{s}^{2}=I_{s}$.

###### Proof.

Since $\begin{bmatrix}X&X_{\perp}\end{bmatrix}$ and
$\begin{bmatrix}Y&Y_{\perp}\end{bmatrix}$ are orthogonal, the result follows
directly from the CS decomposition of the orthogonal matrix
$P=\begin{bmatrix}Y&Y_{\perp}\end{bmatrix}^{T}\begin{bmatrix}X&X_{\perp}\end{bmatrix}$;
see the Theorem of §4 in [22]. ∎

Observe that the matrix $\operatorname{diag}(I_{r},C_{s},O_{p\times p})$ in
this lemma corresponds to the matrix $\cos(\theta)$ in (7) with $\theta$ the
vector of principal angles $0\leq\theta_{1}\leq\cdots\leq\theta_{k}\leq\pi/2$
between $\operatorname{span}(X)$ and $\operatorname{span}(Y)$. However, the
lemma explicitly splits off the angles that are zero and $\pi/2$ so that it
can formulate the related decompositions for $Y^{T}X_{\perp},Y_{\perp}^{T}X,$
and $Y_{\perp}^{T}X_{\perp}$ with $C_{s}$ and $S_{s}$.

We are now ready to state our weak quasi-convexity result. In the statement of
the proposition below (and throughout this paper), we use the convention that
$\frac{0}{\tan 0}=1$.

######  Proposition 7 (Weak-quasi-convexity).

Let $0\leq\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$ be the principal angles
between the subspaces $\mathcal{X}$ and $\mathcal{V}_{\alpha}$. Then, $f$
satisfies
$2a(\mathcal{X})\,(f(\mathcal{X})-f^{*})\leq\langle\textnormal{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V_{\alpha})}\rangle$
with $a(\mathcal{X}):=\theta_{k}/\tan\theta_{k}$.

###### Proof.

Take $X$ and $V_{\alpha}$ matrices with orthonormal columns such that
$\mathcal{X}=\operatorname{span}(X)$ and
$\mathcal{V}_{\alpha}=\operatorname{span}(V_{\alpha})$. Since
$\theta_{k}<\pi/2$, we know that $p=0$ in Lemma 6 and thus $s=k-r$ with $r$
the number of principal angles that are equal to zero. Choosing a matrix
$X_{\perp}$ with orthonormal columns such that
$\operatorname{span}(X_{\perp})=\operatorname{span}(X)^{\perp}$, we therefore
get from Lemma 6 that there exist orthogonal matrices $U_{1},V_{1}$ of size
$k$ and $V_{2}$ of size $n-k$ such that

$\displaystyle V_{\alpha}^{T}X$ $\displaystyle=U_{1}\begin{bmatrix}I_{r}\\\
&C_{k-r}\end{bmatrix}V_{1}^{T},$ $\displaystyle V_{\alpha}^{T}X_{\perp}$
$\displaystyle=U_{1}\begin{bmatrix}O_{r\times m}\\\
&S_{k-r}\end{bmatrix}V_{2}^{T}.$ (15)

Comparing with (7), we deduce that
$C_{k-r}=\operatorname{diag}(\cos\theta_{r+1},\ldots,\cos\theta_{k})$ and
$S_{k-r}=\operatorname{diag}(\sin\theta_{r+1},\ldots,\sin\theta_{k})$ since
$C_{k-r}^{2}+S_{k-r}^{2}=I$.

We recall from (6) that

$\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})=U\operatorname{atan}(\Sigma)V^{T},$
(16)

where $U\Sigma V^{T}=(I_{n}-XX^{T})V_{\alpha}(X^{T}V_{\alpha})^{-1}=:M$ is a
compact SVD (without the requirement that the diagonal of $\Sigma$ is non-
increasing). Using $X_{\perp}$ from above, we can also write
$M=X_{\perp}X_{\perp}^{T}V_{\alpha}(X^{T}V_{\alpha})^{-1}$. Substituting (15)
and using that $U_{1}$ and $V_{1}$ are orthogonal gives

$M=X_{\perp}V_{2}\begin{bmatrix}O_{m\times r\textbf{}}\\\
&S_{k-r}C_{k-r}^{-1}\end{bmatrix}V_{1}^{T}=X_{\perp}\tilde{V}_{2}\begin{bmatrix}O_{r\times
r}\\\ &S_{k-r}C_{k-r}^{-1}\end{bmatrix}V_{1}^{T},$

where $\tilde{V}_{2}\in\mathbb{R}^{(n-k)\times k}$ contains the last $k$
columns of $V_{2}$ in order. Since $\theta_{1}=\cdots=\theta_{r}=0$, we can
therefore formulate the compact SVD of $M$ using the vector $\theta$ of all
principal angles as follows:

$M=U\Sigma V^{T}\quad\text{with $U=X_{\perp}\tilde{V}_{2},\
\Sigma=\tan(\theta),\ V=V_{1}$}.$

Hence from (16) we get directly that

$\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})=X_{\perp}\tilde{V}_{2}\,\theta\,V_{1}^{T},$
(17)

where $\theta$ is a diagonal matrix.

We now claim that (17) also satisfies

$\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})=X_{\perp}X_{\perp}^{T}V_{\alpha}U_{1}\frac{\theta}{\sin\theta}V_{1}^{T},$
(18)

where $\frac{\theta}{\sin\theta}$ is a diagonal matrix for which
$\frac{0}{\sin 0}=1$. Indeed, recalling that $\theta_{1}=\cdots=\theta_{r}=0$
and using the identities

$X_{\perp}^{T}V_{\alpha}=\tilde{V}_{2}\begin{bmatrix}O_{r\times r}\\\
&S_{k-r}\end{bmatrix}U_{1}^{T},\quad\frac{\theta}{\sin\theta}=\begin{bmatrix}I_{r}\\\
&S_{k-r}^{-1}\end{bmatrix}\begin{bmatrix}I_{r}\\\ &T_{k-r}\end{bmatrix}$

where $T_{k-r}=\operatorname{diag}(\theta_{r+1},\ldots,\theta_{k})$, we obtain

rhs of (18) $\displaystyle=X_{\perp}\tilde{V}_{2}\begin{bmatrix}O_{r\times
r}\\\ &S_{k-r}\end{bmatrix}\begin{bmatrix}I_{r}\\\
&S_{k-r}^{-1}\end{bmatrix}\begin{bmatrix}I_{r}\\\
&T_{k-r}\end{bmatrix}\,V_{1}^{T}$
$\displaystyle=X_{\perp}\tilde{V}_{2}\begin{bmatrix}O_{r\times r}\\\
&T_{k-r}\end{bmatrix}\,V_{1}^{T}=X_{\perp}\tilde{V}_{2}\,\theta\,V_{1}^{T}=\textrm{rhs
of~{}\eqref{eq:Log_with_V2}}.$

Next, we work out

$s:=\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V_{\alpha})}\rangle.$

Since $\operatorname{grad}f(\mathcal{X})$ and
$\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})$, respectively, give
tangent vectors for the same representative $X$ of $\mathcal{X}$, the inner
product above is the trace of the corresponding matrix representations. Using
(18) with $I-XX^{T}=X_{\perp}X_{\perp}^{T}$, we therefore get

$\displaystyle s$
$\displaystyle=2\Big{\langle}(I-XX^{T})AX,(I-XX^{T})V_{\alpha}U_{1}\frac{\theta}{\sin(\theta)}V_{1}^{T}\Big{\rangle}$
$\displaystyle=2\operatorname{Tr}\Big{(}\frac{\theta}{\sin(\theta)}U_{1}^{T}V_{\alpha}^{T}(I-XX^{T})AXV_{1}\Big{)}.$

Since $AV_{\alpha}=V_{\alpha}\Lambda_{\alpha}$, we can simplify

$V_{\alpha}^{T}(I-XX^{T})AX=\Lambda_{\alpha}V_{\alpha}^{T}X-V_{\alpha}^{T}XX^{T}AX.$
(19)

Substituting in the expression above and using that
$V_{\alpha}^{T}X=U_{1}\cos\theta\,V_{1}^{T}$, we get

$\displaystyle\tfrac{1}{2}s$
$\displaystyle=\operatorname{Tr}\Big{(}\frac{\theta}{\sin(\theta)}U_{1}^{T}\Lambda_{\alpha}U_{1}\cos(\theta)\Big{)}-\operatorname{Tr}\Big{(}\frac{\theta}{\sin(\theta)}\cos(\theta)V_{1}^{T}X^{T}AXV_{1}\Big{)}$
$\displaystyle=\operatorname{Tr}\Big{(}\frac{\theta}{\tan(\theta)}\Big{(}U_{1}^{T}\Lambda_{\alpha}U_{1}-V_{1}^{T}X^{T}AXV_{1}\Big{)}\Big{)},$

with the convention $\frac{0}{\tan 0}=1$.

Denote the symmetric matrix

$S:=U_{1}^{T}\Lambda_{\alpha}U_{1}-V_{1}^{T}X^{T}AXV_{1}.$ (20)

We show below that all diagonal entries $S_{11},\ldots,S_{kk}$ of $S$ are
positive. Hence, by diagonality of the matrix $\tfrac{\theta}{\tan(\theta)}$,
we obtain

$\displaystyle\tfrac{1}{2}s$
$\displaystyle=\sum_{j}\frac{\theta_{j}}{\tan\theta_{j}}\,S_{jj}\geq\min_{j}\frac{\theta_{j}}{\tan\theta_{j}}\,\operatorname{Tr}(S)=\frac{\theta_{k}}{\tan\theta_{k}}\,\Big{[}\operatorname{Tr}(\Lambda_{\alpha})-\operatorname{Tr}(X^{T}AX)\Big{]}$

since $U_{1}$ and $V_{1}$ are orthogonal matrices. We recover the desired
result after substituting $f(\mathcal{X})=-\operatorname{Tr}(X^{T}AX)$ and
$f_{*}=-\operatorname{Tr}(V_{\alpha}^{T}AV_{\alpha})=-\operatorname{Tr}(\Lambda_{\alpha})$,

It remains to show that $S_{jj}\geq 0$ for $j=1,\ldots,k$. Since
$\operatorname{span}(V_{\beta})=\operatorname{span}(V_{\alpha})^{\perp}$,
Lemma 6 gives us in addition to (15) also

$V_{\beta}^{T}X=U_{2}\begin{bmatrix}O_{m\times r}\\\
&S_{r}\end{bmatrix}V_{1}^{T}=\tilde{U}_{2}\sin\theta\,V_{1}^{T},$ (21)

where $\tilde{U}_{2}\in\mathbb{R}^{(n-k)\times k}$ contains the last $k$
columns of the orthogonal matrix $U_{2}$ in order. A short calculation using
(13) then shows that (20) satisfies

$S=U_{1}^{T}\Lambda_{\alpha}U_{1}-\cos\theta\,U_{1}^{T}\Lambda_{\alpha}U_{1}\cos\theta-\sin\theta\,\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2}\sin\theta$

with diagonal elements

$S_{jj}=\sin^{2}\theta_{j}\,(U_{1}^{T}\Lambda_{\alpha}U_{1}-\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2})_{jj}.$

Since $U_{1}$ and $\tilde{U}_{2}$ have orthonormal columns, we obtain

$\lambda_{\min}(U_{1}^{T}\Lambda_{\alpha}U_{1})\geq\lambda_{\min}(\Lambda_{\alpha})=\lambda_{k},\quad\lambda_{\max}(\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2})\leq\lambda_{\max}(\Lambda_{\beta})=\lambda_{k+1},$

from which we get with Weyl’s inequality that

$\lambda_{\min}(U_{1}^{T}\Lambda_{\alpha}U_{1}-\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2})\geq\lambda_{\min}(U_{1}^{T}\Lambda_{\alpha}U_{1})-\lambda_{\max}(\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2})\geq\lambda_{k}-\lambda_{k+1}\geq
0.$

Hence, the matrix

$U_{1}^{T}\Lambda_{\alpha}U_{1}-\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2}$
(22)

is symmetric and positive semi-definite. Its diagonal entries, and thus also
$S_{jj}$, are therefore positive. ∎

We now arrive at a useful property of $f$ that will later allow us to analyze
the convergence of Riemannian steepest descent. It is a _weaker version of
strong geodesic convexity_ and can be proved easily using quadratic growth and
weak-quasi-convexity.

######  Theorem 8 (Weak-strong convexity).

Let $0\leq\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$ be the principal angles
between the subspaces $\mathcal{X}$ and $\mathcal{V}_{\alpha}$. Then, $f$
satisfies
$f(\mathcal{X})-f^{*}\leq\frac{1}{a(\mathcal{X})}\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\rangle-
c_{Q}\delta\,\operatorname{dist}^{2}(\mathcal{X},\mathcal{V}_{\alpha})$ with
$a(\mathcal{X})=\theta_{k}/\tan\theta_{k}>0$, $c_{Q}=4/\pi^{2}>0.4$, and
$\delta=\lambda_{k}-\lambda_{k+1}\geq 0$.

###### Proof.

Combining Propositions 5 and 7 leads to

$c_{Q}\delta\,\operatorname{dist}^{2}(\mathcal{X},\mathcal{V}_{\alpha})\leq
f(\mathcal{X})-f^{*}\leq\frac{1}{2a(\mathcal{X})}\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\rangle.$

At the same time, Proposition 7 also implies

$\displaystyle f(\mathcal{X})-f^{*}$
$\displaystyle\leq\frac{1}{2a(\mathcal{X})}\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\rangle-
c_{Q}\delta\,\operatorname{dist}^{2}(\mathcal{X},\mathcal{V}_{\alpha})$
$\displaystyle\qquad+c_{Q}\delta\,\operatorname{dist}^{2}(\mathcal{X},\mathcal{V}_{\alpha}).$

Substituting this inequality in the first one gives the desired result. ∎

######  Remark.

Theorem 8 is also valid when the eigengap $\delta=0$. In that case,
$\mathcal{V}_{\alpha}$ is any subspace spanned by $k$ leading eigenvectors of
$A$ and the theorem reduces to Proposition 7.

While not needed for our convergence proof, the next result is of independent
interest and shows that $f$ is gradient dominated in the Riemannian sense when
the eigengap $\delta$ is strictly positive. This property is the _Riemannian
version of the Polyak–Łojasiewicz inequality_ and generalizes a result by [26]
for the Rayleigh quotient on the sphere.

######  Proposition 9 (Gradient dominance).

The function $f$ satisfies $\|\textnormal{grad}f(\mathcal{X})\|^{2}\geq
4\,c_{Q}\,\delta\,a^{2}(\mathcal{X})(f(\mathcal{X})-f^{*})$ for all subspaces
$\mathcal{X}$ that have a largest principal angle $<\pi/2$ with
$\mathcal{V}_{\alpha}$.

###### Proof.

We assume that $\delta>0$ since otherwise the statement is trivially true. By
Theorem 8, we have

$f(\mathcal{X})-f^{*}\leq\frac{1}{a(\mathcal{X})}\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\rangle-
c_{Q}\delta\operatorname{dist}^{2}(\mathcal{X},\mathcal{V}_{\alpha}).$

Since $\langle G_{1},G_{2}\rangle\leq(\|G_{1}\|^{2}+\|G_{2}\|^{2})/2$ for all
matrices $G_{1},G_{2}$, we can write for any $\rho>0$ that

$\langle\operatorname{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\rangle\leq\frac{\rho}{2}\|\textnormal{grad}f(\mathcal{X})\|^{2}+\frac{1}{2\rho}\|\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\|^{2}.$

Using that
$\operatorname{dist}(\mathcal{X},\mathcal{V}_{\alpha})=\|\operatorname{Log}_{\mathcal{X}}(\mathcal{V}_{\alpha})\|$
and choosing $\rho=1/(2c_{Q}\delta a(\mathcal{X}))$, we get the desired
result. ∎

##  5 Convergence of Riemannian steepest descent

We now have everything in place to prove the convergence of the Riemannian
steepest descent (RSD) method on the Grassmann manifold for minimizing $f$.
Starting from a subspace $\mathcal{X}_{0}\in\operatorname{Gr}(n,k)$, we
iterate

$\mathcal{X}_{t+1}=\operatorname{Exp}_{\mathcal{X}_{t}}(-\eta_{t}\,\operatorname{grad}f(\mathcal{X}_{t})).$
(23)

Here, $\eta_{t}>0$ is a step size that may depend on the iteration $t$ and
will be carefully chosen depending on the specific case.

We start by a general result which shows that the distance to the optimal
subspace contracts after one step of steepest descent. The step size depends
on the smoothness and weak-quasi-convexity constants of $f$ from Propositions
3 and 7. This is crucial since the constant $a(\mathcal{X})$ depends on the
biggest principal angle between $\mathcal{X}$ and $\mathcal{V}_{\alpha}$ and
bounding the evolution of distances of the iterates to the minimizer will help
us also bounding this constant.555The analysis of [16] is wrong with respect
to this issue as discussed in detail in [4]. An alternative contraction
property with a more tractable step size is presented in Proposition 20 of
Appendix B.

######  Lemma 10 (Contraction of RSD).

Let $\mathcal{X}_{t}$ and $\mathcal{V}_{\alpha}$ have principal angles
$0\leq\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$. Then, iteration (23) with
$0\leq\eta_{t}\leq\frac{a(\mathcal{X}_{t})}{\gamma}$ satisfies
$\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\big{(}1-2c_{Q}\delta
a(\mathcal{X}_{t})\,\eta_{t}\big{)}\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$

Observe that $\gamma=0$ implies $A=\lambda_{1}I$ and any subspace
$\mathcal{X}$ of dimension $k$ will be an eigenspace of $A$ with
$\operatorname{dist}(\mathcal{X},\mathcal{V}_{\alpha})=0$. We will therefore
not explicitly prove this lemma and all forthcoming convergence results for
$\gamma=0$ since the statements will be trivially true.

###### Proof of Lemma 10.

By the assumption on the principal angles, we get that
$0<a(\mathcal{X}_{t})=\theta_{k}/\tan\theta_{k}\leq 1$. The hypothesis on
$\eta_{t}$ and Lemma 4 then gives

$\eta_{t}\|\operatorname{grad}f(\mathcal{X}_{t})\|_{2}\leq\frac{a(\mathcal{X}_{t})}{\gamma}\|\operatorname{grad}f(\mathcal{X}_{t})\|_{2}\leq\frac{1}{2}<\frac{\pi}{2}.$

By (5), this guarantees that the geodesic
$\tau\mapsto\operatorname{Exp}(-\tau\eta_{t}\,\operatorname{grad}f(\mathcal{X}_{t}))$
lies within the injectivity domain at ${\mathcal{X}_{t}}$ for $\tau\in[0,1]$.
Hence, $\operatorname{Exp}$ is bijective along this geodesic and thus
$\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1})=-\eta_{t}\,\operatorname{grad}f(\mathcal{X}_{t})$.
We can thus apply Lemma 1 to obtain

$\displaystyle\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})$
$\displaystyle\leq\|-\eta_{t}\operatorname{grad}f(\mathcal{X}_{t})-\textnormal{log}_{\mathcal{X}_{t}}(\mathcal{V}_{\alpha})\|^{2}$
$\displaystyle=\eta_{t}^{2}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}+\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})+2\eta_{t}\,\sigma\
$ (24)

with

$\sigma:=\langle\operatorname{grad}f(\mathcal{X}_{t}),\log_{\mathcal{X}_{t}}(\mathcal{V}_{\alpha})\rangle.$

Theorem 8 and (11) together with Proposition 3 give

$\displaystyle\frac{\sigma}{a(\mathcal{X}_{t})}$ $\displaystyle\leq
f^{*}-f(\mathcal{X}_{t})-c_{Q}\delta\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$
$\displaystyle\leq-\frac{1}{2\gamma}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}-c_{Q}\delta\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha}).$

Multiplying by $2a(\mathcal{X}_{t})\,\eta_{t}$ and using $\eta_{t}\leq
a(\mathcal{X}_{t})/\gamma$, we get

$\displaystyle 2\eta_{t}\,\sigma$
$\displaystyle\leq-\frac{a(\mathcal{X}_{t})\,\eta_{t}}{\gamma}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}-2c_{Q}\delta
a(\mathcal{X}_{t})\,\eta_{t}\,\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$
$\displaystyle\leq-\eta_{t}^{2}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}-2c_{Q}\delta
a(\mathcal{X}_{t})\,\eta_{t}\,\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha}).$

Substituting into (5), we obtain the first statement of the lemma. ∎

######  Remark.

When $\delta=0$, Lemma 10 still holds _for any subspace $\mathcal{V}_{\alpha}$
spanned by $k$ leading eigenvectors of $A$_. In that case, the lemma only
guarantees that the distance between the iterates of steepest descent and this
$\mathcal{V}_{\alpha}$ does not increase.

###  5.1 Linear convergence rate under positive eigengap

When $\delta>0$, we can extend Lemma 10 to a linear convergence rate of
distances to the minimizer:

######  Theorem 11.

If $\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$ then the
iterates $\mathcal{X}_{t}$ of Riemannian steepest descent (23) with step size
$\eta_{t}$ such that
$0<\eta\leq\eta_{t}\leq\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))/\gamma$
satisfy
$\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\left(1-c_{Q}\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\,\delta\,\eta\right)^{t}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha}).$

###### Proof.

We first claim that
$\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$
for all $t\geq 0$. This would then also imply that
$\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha})<\pi/2$ for all $t\geq 0$
since

$\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\sqrt{\sum_{i=1}^{k}\theta_{i}(\mathcal{X}_{t},\mathcal{V}_{\alpha})^{2}}=\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha}).$

For $t=0$, we have $\theta_{k}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$ by
hypothesis on $\mathcal{X}_{0}$ and thus

$a(\mathcal{X}_{0})=\frac{\theta_{k}(\mathcal{X}_{0},\mathcal{V}_{\alpha})}{\tan(\theta_{k}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))}\geq\cos(\theta_{k}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\geq\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})).$

Since by construction
$\eta_{0}\leq\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))/\gamma$
, this implies that $\eta_{0}\leq a(\mathcal{X}_{0})/\gamma$ and Lemma 10
guarantees that
$\operatorname{dist}(\mathcal{X}_{1},\mathcal{V}_{\alpha})\leq\textnormal{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$.
In particular, we also have
$\theta_{k}(\mathcal{X}_{1},\mathcal{V}_{\alpha})<\pi/2$.

Next, assume that

$\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}),$

which implies $\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha})<\pi/2$. Then
by a similar argument like above, we have

$a(\mathcal{X}_{t})\geq\cos(\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha}))\geq\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})).$
(25)

By hypothesis on $\eta_{t}$, we observe

$\eta_{t}\leq\frac{\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))}{\gamma}\leq\frac{\cos(\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha}))}{\gamma}\leq\frac{a(\mathcal{X}_{t})}{\gamma}.$

Applying Lemma 10 once again with the induction hypothesis proves the claim:

$\operatorname{dist}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\operatorname{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}).$

The main statement of the theorem now follows easily: Since $\eta_{t}\leq
a(\mathcal{X}_{t})/\gamma$ and
$\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha})<\pi/2$ for all $t\geq 0$,
Lemma 10 gives

$\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\left(1-c_{Q}a(\mathcal{X}_{t})\delta\eta_{t}\right)\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$

Combining with (25) and $\eta_{t}\geq\eta$ shows the desired result by
induction. ∎

If the eigengap $\delta$ is strictly positive, then Theorem 11 gives an
exponential convergence rate towards the optimum $\mathcal{V}_{\alpha}$. If
$\delta=0$, then Theorem 11 _does not provide a convergence rate_ but rather
implies that the intrinsic distances of the iterates to the optimum do not
increase.

From Theorem 11 we get immediately the following iteration complexity.

######  Corollary 12.

Let Riemannian steepest descent be started from a subspace $\mathcal{X}_{0}$
that satisfies
$\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$. Then after
at most

$T\leq
2\frac{\log(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))-\log(\varepsilon)}{\log(1-0.4\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\delta\eta)}+1=\mathcal{O}\left(\frac{\log(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))-\log(\varepsilon)}{\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\delta\eta}\right)$

many iterations, $\mathcal{X}_{T}$ will satisfy
$\operatorname{dist}(\mathcal{X}_{T},\mathcal{V}_{\alpha})\leq\varepsilon$.
With the maximal step size allowed in Theorem 11, we get

$T=\mathcal{O}\left(\frac{\lambda_{1}-\lambda_{n}}{\delta}\frac{1}{\cos^{2}(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))}\log\left(\frac{\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})}{\varepsilon}\right)\right).$

As expected, $T$ depends inversely proportional on the eigengap $\delta$ and
proportional to the spread of the eigenvalues. In addition, we also have an
extra term
$1/\cos^{2}(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))$ that
depends on the initial distance
$\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$, which is due to
the weak-quasi-convexity property of $f$. This is a conservative
overestimation, since this quantity improves as the iterates get closer to the
optimum.

######  Remark.

If $\delta>0$, the exponential convergence rate is in terms of the intrinsic
distance on the Grassmann manifold, that is, the $\ell_{2}$ norm of the
principal angles. Standard convergence results for subspace iteration are
stated for the biggest principal angle, that is, the $\ell_{\infty}$ norm.
This is weaker than the intrinsic distance. For subspace iteration with
projection, the convergence result from [24, Thm. 5.2] shows that all
principal angles $\theta_{i}$ converge to zero and eventually gives
convergence of the $\ell_{4}$ norm of the principal angles. This is also
weaker than the intrinsic distance.

###  5.2 Convergence of function values without an eigengap assumption

When $\delta=0$, Theorem 11 still holds, but does not provide a rate of
convergence as discussed above. Instead, we can prove the following result:

######  Theorem 13.

If the distance $\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$ of
the initial subspace $\mathcal{X}_{0}$ to the minimizer satisfies
$\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$ for a
subspace $\mathcal{V}_{\alpha}$ that is spanned by any $k$ leading
eigenvectors of $A$, then the iterates $\mathcal{X}_{t}$ of Riemannian
steepest descent (23) with fixed step size
$\eta\leq\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))/\gamma$
satisfy
$f(\mathcal{X}_{t})-f^{*}\leq\frac{2\gamma+\frac{1}{\eta}}{\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))t+1}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})=\mathcal{O}\left(\frac{1}{t}\right).$

###### Proof.

Since we satisfy all the hypotheses of Theorem 11, we know that for all $t\geq
0$ it holds
$\textnormal{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\textnormal{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$
and thus also that $\mathcal{X}_{t}$ is in the injectivity domain of
$\operatorname{Exp}$ at $\mathcal{V}_{\alpha}$. In addition, its proof states
in (25) that

$a(\mathcal{X}_{t})\geq
C_{0}:=\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))>0,$

which implies that the function $f$ is weakly-quasi-convex at every
$\mathcal{X}_{t}$ with constant $2C_{0}$. Hence

$2C_{0}\Delta_{t}\leq\langle\operatorname{grad}f(\mathcal{X}_{t}),-\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{V}_{\alpha})\rangle,$
(26)

where we defined

$\Delta_{t}:=f(\mathcal{X}_{t})-f^{*}.$

Similar to the proof of Theorem 11, by the hypothesis on the step size
$\eta_{t}$, Lemma 10 shows that $\mathcal{X}_{t+1}$ is in the injectivity
domain of $\operatorname{Exp}$ at $\mathcal{X}_{t}$. Hence, by the definition
of Riemannian steepest descent, we have

$\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1})=-\eta\operatorname{grad}f(\mathcal{X}_{t}).$
(27)

In addition, the smoothness property (10) of $f$ gives

$\Delta_{t+1}-\Delta_{t}\leq\langle\operatorname{grad}f(\mathcal{X}_{t}),\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1})\rangle+\frac{\gamma}{2}\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{X}_{t+1}).$

Substituting (27), we obtain

$\Delta_{t+1}-\Delta_{t}\leq\left(-\eta+\frac{\gamma}{2}\eta^{2}\right)\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}\leq
0$ (28)

since $\eta\leq C_{0}/\gamma$ with
$0<C_{0}:=\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\leq
1$ and $\gamma>0$.

Since $\operatorname{Gr}(n,k)$ has nonnegative sectional curvature, Lemma 2
implies

$\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{X}_{t+1})+\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-2\langle\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1}),\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{V}_{\alpha})\rangle.$

Substituting (27) into the above and rearranging terms gives

$2\eta\langle\operatorname{grad}f(\mathcal{X}_{t}),-\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{V}_{\alpha})\rangle\leq\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})+\eta^{2}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}.$

Combining with (26), we get

$\Delta_{t}\leq\frac{1}{4C_{0}\eta}(\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha}))+\frac{\eta}{4C_{0}}\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}.$
(29)

Now multiplying (28) by $\frac{1}{C_{0}}$ and summing with (29) gives

$\frac{1}{C_{0}}\Delta_{t+1}-\left(\frac{1}{C_{0}}-1\right)\Delta_{t}\leq\frac{1}{4C_{0}\eta}(\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha}))\\\
+\frac{1}{C_{0}}\left(-\eta+\frac{\gamma}{2}\eta^{2}+\frac{\eta}{4}\right)\|\operatorname{grad}f(\mathcal{X}_{t})\|^{2}.$
(30)

By assumption $\eta\leq C_{0}/\gamma$ where
$0<C_{0}:=\cos(\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha}))\leq
1$ and $\gamma>0$. Since

$\frac{\eta}{C_{0}}\left(-1+\frac{\gamma}{2}\eta+\frac{1}{4}\right)\leq\frac{1}{\gamma}\left(\frac{C_{0}}{2}-\frac{3}{4}\right)\leq-\frac{1}{4\gamma}<0.$

the inequality (30) can be simplified to

$\frac{1}{C_{0}}\Delta_{t+1}-\left(\frac{1}{C_{0}}-1\right)\Delta_{t}\leq\frac{1}{4C_{0}\eta}(\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-\operatorname{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})).$

Summing from $0$ to $t-1$ gives

$\frac{1}{C_{0}}\Delta_{t}+\sum_{s=1}^{t-1}\Delta_{s}-\left(\frac{1}{C_{0}}-1\right)\Delta_{0}\leq\frac{1}{4C_{0}\eta}\left(\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})-\operatorname{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\right).$

From the smoothness property (10) at the critical point $\mathcal{V}_{\alpha}$
of $f$, we get

$\Delta_{0}\leq\frac{\gamma}{2}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha}).$

Combining these two inequalities then leads to

$\displaystyle\frac{1}{C_{0}}\Delta_{t}+\sum_{s=0}^{t-1}\Delta_{s}$
$\displaystyle\leq\frac{1}{C_{0}}\Delta_{0}+\frac{1}{4C_{0}\eta}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$
$\displaystyle=\frac{1}{2C_{0}}\left(\gamma+\frac{1}{2\eta}\right)\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha}).$

Since (28) holds for all $t\geq 0$, it also implies $\Delta_{t}\leq\Delta_{s}$
for all $1\leq s\leq t$. Substituting

$t\Delta_{t}\leq\sum_{s=0}^{t-1}\Delta_{s}.$

into the inequality from above,

$\Delta_{t}\leq\frac{1}{2C_{0}}\frac{\gamma+\frac{1}{2\eta}}{\frac{1}{C_{0}}+t}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})=\frac{\gamma+\frac{1}{2\eta}}{2(C_{0}t+1)}\operatorname{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha}),$

we obtain the desired result. ∎

######  Remark.

This type of result is standard for functions that are geodesically convex
(see, e.g. [27]). Our objective function does not satisfy this property, but
we can still have a similar upper bound on the iteration complexity for
convergence in function value. We note that this does not imply convergence of
the iterates to a specific $k$-dimensional subspace, but only convergence of a
subsequence of the sequence of the iterates.

###  5.3 Sufficiently small step sizes

The convergence results in Theorems 11 and 13 require that the initial
subspace $\mathcal{X}_{0}$ lies within a distance strictly less than $\pi/2$
from a global minimizer $\mathcal{V}_{\alpha}$. While this condition is
independent from the eigengap (unlike results that rely on standard convexity,
see appendix), it is also not fully satisfactory: it is hard to verify in
practice, and it is unnecessarily severe in numerical experiments. In fact,
this condition is only used to obtain a uniform lower bound on the weak-quasi-
convexity constant
$a(\mathcal{X}_{t})=\theta_{k}^{(t)}/\tan(\theta_{k}^{(t)})$ with
$\theta_{k}^{(t)}$ the largest principal angle between $\mathcal{X}_{t}$ and
$\mathcal{V}_{\alpha}$. Since the Riemannian distance is the $\ell_{2}$ norm
of the principal angles, a contraction in this distance leads automatically to
$\theta_{k}^{(t)}<\pi/2$ if $\theta_{k}^{(0)}<\pi/2$. If one could guarantee
by some other reasoning that $\theta_{k}^{(t)}$ does not increase after one
step, the condition
$\operatorname{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})<\pi/2$ would not be
needed.

We now show that for sufficiently small step sizes $\eta_{t}$, the largest
principal angle $\theta_{k}^{(t)}$ between $\mathcal{X}_{t}$ and
$\mathcal{V}_{\alpha}$ does indeed not increase after each iteration of
Riemannian steepest descent regardless of the initial subspace
$\mathcal{X}_{0}$. While it does not explain what we observe in numerical
experiments where large steps can be taken, it is a first result in explaining
why we can initialize the iteration at a random initial subspace
$\mathcal{X}_{0}$.

######  Proposition 14.

Riemannian steepest descent started from a subspace $\mathcal{X}_{t}$ returns
a subspace $\mathcal{X}_{t+1}$ such that
$\theta_{k}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha}),$
for all step sizes $0\leq\eta\leq\bar{\eta}$ where $\bar{\eta}>0$ is
sufficiently small.

For the proof of this proposition, we will need the derivatives of certain
singular values. While this is well known for isolated singular values, it is
possible to generalize to higher multiplicities as well by relaxing the
ordering and sign of singular values [9]. For a concrete formula, we use the
following result from Lemma A.5 in [19].

######  Lemma 15.

Let $\sigma_{1}\geq\cdots\geq\sigma_{n}$ be the singular values of
$A\in\mathbb{R}^{n\times n}$ with $u_{1},\ldots,u_{n}$ and
$v_{1},\ldots,v_{n}$ the associated left and right orthonormal singular
vectors. Suppose that $\sigma_{j}$ has multiplicity $m$, that is,

$\sigma_{j_{0}-1}>\sigma_{j_{0}}=\cdots=\sigma_{j}=\cdots=\sigma_{j_{0}+m-1}>\sigma_{j_{0}+m}.$

Then, the $j$th singular value of $A+\eta B$ satisfies

$\sigma_{j}(A+\eta
B)=\sigma_{j}+\eta\,\lambda_{j-j_{0}+1}+\mathcal{O}(\eta^{2}),\quad\eta\to
0^{+},$

where $\lambda_{j}$ is the $j$th largest eigenvalue of
$\tfrac{1}{2}(U^{T}BV+V^{T}B^{T}U)$ with

$U=\begin{bmatrix}u_{j_{0}}&\cdots&u_{j_{0}+m-1}\end{bmatrix}\quad\text{and}\quad
V=\begin{bmatrix}v_{j_{0}}&\cdots&v_{j_{0}+m-1}\end{bmatrix}.$

###### Proof of Proposition 14.

For ease of notation, let $X:=X_{t}$ and $X_{+}:=X_{t+1}$ such that
$\mathcal{X}_{t}=\operatorname{span}(X)$ and
$\mathcal{X}_{t+1}=\operatorname{span}(X_{+})$. By definition of the
exponential map on Grassmann, the next iterate of the Riemannian SD iteration
(23) with step $\eta$ satisfies

$X_{+}=XV\cos(\eta\Sigma)V^{T}+U\sin(\eta\Sigma)V^{T}$

with

$U\Sigma V^{T}=-\operatorname{grad}f(\mathcal{X}_{t}).$

Since $V$ is orthogonal, we can write

$U\sin(\eta\Sigma)V^{T}=U(\eta\Sigma)V^{T}V\left(\frac{\sin(\eta\Sigma)}{\eta\Sigma}\right)V^{T}=-\eta\operatorname{grad}f(\mathcal{X}_{t})V\left(\frac{\sin(\eta\Sigma)}{\eta\Sigma}\right)V^{T}$

where $1/\Sigma:=\Sigma^{-1}$ and $\frac{\sin 0}{0}=1$. Taking Taylor
expansions of $\sin$ and $\cos$,

$\displaystyle V\cos(\eta\Sigma)V^{T}$
$\displaystyle=V\left(I-\mathcal{O}(\eta^{2})\right)V^{T}=I-\mathcal{O}(\eta^{2})$
$\displaystyle V\frac{\sin(\eta\Sigma)}{\eta\Sigma}V^{T}$
$\displaystyle=V\left(I-\mathcal{O}(\eta^{2})\right)V^{T}=I-\mathcal{O}(\eta^{2}),$

we obtain

$\displaystyle{V}_{\alpha}^{T}X_{+}$
$\displaystyle={V}_{\alpha}^{T}X(I-\mathcal{O}(\eta^{2}))+{V}_{\alpha}^{T}(-\eta\operatorname{grad}f(\mathcal{X}))(I-\mathcal{O}(\eta^{2}))$
$\displaystyle={V}_{\alpha}^{T}(X-\eta\operatorname{grad}f(\mathcal{X}_{t}))(I-\mathcal{O}(\eta^{2}))$
(31)

since $\|V_{\alpha}\|_{2}=\|X\|_{2}=1$.

Let now $\theta$ be the vector of $k$ principal angles between
$\mathcal{X}_{t}$ and $\mathcal{V}_{\alpha}$. As in (14) and (21), we
therefore have the SVDs

$V_{\alpha}^{T}X=U_{1}\cos\theta\,V_{1}^{T}\qquad\text{and}\qquad
V_{\beta}^{T}X=\tilde{U}_{2}\sin\theta\,V_{1}^{T},$ (32)

where $U_{1},V_{1}\in\mathbb{R}^{k\times k}$ and
$\tilde{U}_{2}\in\mathbb{R}^{(n-k)\times k}$ have orthonormal columns. Next,
we write (5.3) in terms of

$M:=\sin^{2}\theta\,U_{1}^{T}\Lambda_{\alpha}U_{1}\cos\theta-\cos\theta\sin\theta\,\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2}\sin\theta.$

Since $\operatorname{grad}f(\mathcal{X}_{t})=-2(I-XX^{T})AX$, the identity
(19) gives

${V}_{\alpha}^{T}(X-\eta\operatorname{grad}f(\mathcal{X}_{t}))=V_{\alpha}^{T}X+2\eta\Lambda_{\alpha}V_{\alpha}^{T}X-2\eta
V_{\alpha}^{T}XX^{T}AX.$

After substituting (13) and (32), a short calculation using
$\cos^{2}\theta=I-\sin^{2}\theta$ and the orthogonality of $U_{1}$ and $V_{1}$
then shows

$V_{\alpha}^{T}(X-\eta\operatorname{grad}f(\mathcal{X}_{t}))=U_{1}(\cos\theta+2\eta
M)V_{1}^{T}.$

Relating back to (5.3), we thus obtain

$\displaystyle V_{\alpha}^{T}{X}_{+}$ $\displaystyle=U_{1}(\cos\theta+2\eta
M)V_{1}^{T}(I-\mathcal{O}(\eta^{2}))$ $\displaystyle=U_{1}(\cos\theta+2\eta
M)(I-V_{1}^{T}\mathcal{O}(\eta^{2})V_{1})V_{1}^{T}$
$\displaystyle=U_{1}(\cos\theta+2\eta M-\mathcal{O}(\eta^{2}))V_{1}^{T}.$

The singular values of $V_{\alpha}^{T}{X}_{+}$ are therefore the same as the
singular values of the matrix $\cos\theta+2\eta M+\mathcal{O}(\eta^{2})$.

By Weyl’s inequality (see, e.g., [15, Cor. 7.3.5]), each singular value of
$\cos\theta+2\eta M+\mathcal{O}(\eta^{2})$ is $\mathcal{O}(\eta^{2})$ close to
some singular value of $\cos\theta+2\eta M$. Let $1\leq j\leq k$. Denote the
$j$th singular of $\cos\theta+2\eta M$ by $\sigma_{j}(\eta)$ to which we will
apply Lemma 15. Since $\cos\theta$ is a diagonal matrix with decreasing
diagonal, its $j$th singular value equals $\cos\theta_{j}$ and its associated
left/right singular vector is the $j$th canonical vector $e_{j}$. Denoting

$E=\begin{bmatrix}e_{j_{0}}&\cdots&e_{j_{0}+m-1}\end{bmatrix},$

observe that $\cos\theta\,E=\cos\theta_{j_{0}}\,E$ (here, $\cos\theta$ is a
diagonal matrix and $\cos\theta_{j_{0}}$ is a scalar) and likewise for
$\sin\theta\,E$. We thus get

$E^{T}ME=\sin^{2}\theta_{j_{0}}\cos\theta_{j_{0}}(U_{1}^{T}\Lambda_{\alpha}U_{1}-\tilde{U}_{2}^{T}\Lambda_{\beta}\tilde{U}_{2}).$

In the proof of Proposition 7, we showed that the matrix in brackets above is
symmetric and positive semi-definite (see (22)). Since
$0\leq\theta_{j_{0}}\leq\pi/2$, the eigenvalues of $E^{T}ME$ are therefore all
non-negative. Lemma 15 thus gives that $\sigma_{j}(\eta)\geq\sigma_{j}$ for
sufficiently small and positive $\eta$. Since the singular values of
$V_{\alpha}^{T}X_{+}$ are the cosines of the principal angles between
$\mathcal{V}_{\alpha}$ and $\mathcal{X}_{t+1}$ with step size $\eta\geq 0$, we
conclude that there exists $\bar{\eta}>0$ such that for all
$\eta\in[0,\bar{\eta}]$ it holds

$\theta_{j}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})\leq\theta_{j}(\mathcal{X}_{t},\mathcal{V}_{\alpha}).$

Since $j$ was arbitrary, this finishes the proof. ∎

##  6 Conclusion and future work

We provided the first systematic study of Riemannian steepest descent on the
Grassmann manifold for computing a subspace spanned by $k$ leading
eigenvectors of a symmetric and positive semi-definite matrix $A$.

Our main idea was to exploit a convexity-like structure of the block Rayleigh
quotient, which can be of much more general interest than for only analyzing
steepest descent. One example is line search methods, which have usually
favourable properties compared to vanilla steepest descent. Also, weakly-
quasi-convex functions have been proven to admit accelerated algorithms [20],
while accelerated or almost accelerated Riemannian algorithms have been
developed in [28, 5, 6]. It would naturally be interesting to examine whether
a provable accelerated method can be developed for the block Rayleigh quotient
on the Grassmann manifold. This would hopefully reduce the dependence of the
iteration complexity on the eigengap $\delta$ from $\mathcal{O}(1/\delta)$ to
$\mathcal{O}(1/\sqrt{\delta})$.

Another interesting direction is to extend the analysis of [4] from the
computation of just one leading eigenvector to computation of a whole
subspace, using the generalized machinery developed in this work, or develop a
noisy version of steepest descent and compare with noisy power method [12].

#### Acknowledgements

This work was supported by the SNSF under research project 192363.

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre.  Riemannian geometry of Grassmann manifolds with a view on algorithmic computation.  Acta Applicandae Mathematicae, 80(2):199–220, 2004.
  * [2] P.-A. Absil, R. Mahony, and R. Sepulchre.  Optimization Algorithms on Matrix Manifolds.  Princeton University Press, Princeton, NJ, 2008.
  * [3] Kwangjun Ahn and Felipe Suarez.  Riemannian perspective on matrix factorization.  arXiv preprint arXiv:2102.00937, 2021.
  * [4] Foivos Alimisis, Peter Davies, Bart Vandereycken, and Dan Alistarh.  Distributed principal component analysis with limited communication.  Advances in Neural Information Processing Systems, 34, 2021.
  * [5] Foivos Alimisis, Antonio Orvieto, Gary Bécigneul, and Aurelien Lucchi.  A continuous-time perspective for modeling acceleration in riemannian optimization.  In International Conference on Artificial Intelligence and Statistics, pages 1297–1307. PMLR, 2020.
  * [6] Foivos Alimisis, Antonio Orvieto, Gary Becigneul, and Aurelien Lucchi.  Momentum improves optimization on riemannian manifolds.  In International Conference on Artificial Intelligence and Statistics, pages 1351–1359. PMLR, 2021.
  * [7] Thomas Bendokat, Ralf Zimmermann, and P.-A. Absil.  A Grassmann manifold handbook: Basic geometry and computational aspects.  arXiv:2011.13699 [cs, math], December 2020.
  * [8] Nicolas Boumal.  An introduction to optimization on smooth manifolds.  To appear with Cambridge University Press, Apr 2022.
  * [9] A. Bunse-Gerstner, R. Byers, V. Mehrmann, and N. K. Nichols.  Numerical computation of an analytic singular value decomposition of a matrix valued function.  Numerische Mathematik, 60(1):1–39, 1991.
  * [10] A. Edelman, T. A. Arias, and S. T. Smith.  The geometry of algorithms with orthogonality constraints.  SIAM Journal on Matrix Analysis and Applications, 20(2):303–353, 1999.
  * [11] Gene H Golub and Charles F Van Loan.  Matrix computations.  JHU press, 2013.
  * [12] Moritz Hardt and Eric Price.  The noisy power method: A meta algorithm with applications.  Advances in neural information processing systems, 27, 2014.
  * [13] N. J. Higham and S. Cheng.  Modifying the inertia of matrices arising in optimization.  Lin. Alg. Appl., 275–276:261–279, 1998.
  * [14] Roger Horn and Charles R. Johnson.  Topics in Matrix Analysis.  Cambridge University Press, 1991.
  * [15] Roger A. Horn and Charles R. Johnson.  Matrix Analysis.  Cambridge University Press, Cambridge ; New York, 2nd ed edition, 2012\.
  * [16] Long-Kai Huang and Sinno Pan.  Communication-efficient distributed PCA by Riemannian optimization.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 4465–4474. PMLR, 13–18 Jul 2020\.
  * [17] J. Kuczynski and H. Wozniakowski.  Estimating the largest eigenvalue by the power and lanczos algorithms with a random start, 1992.
  * [18] Chi-Kwong Li and Roy Mathias.  Inequalities on the Singular Values of an Off-Diagonal Block of a Hermitian Matrix.  Journal of Inequalities and Applications, 3(2):137–142, 1999.
  * [19] Ross A. Lippert.  Fixing two eigenvalues by a minimal perturbation.  Linear Algebra and its Applications, 406:177–200, September 2005\.
  * [20] Yurii Nesterov, Alexander Gasnikov, Sergey Guminov, and Pavel Dvurechensky.  Primal–dual accelerated gradient methods with small-dimensional relaxation oracle.  Optimization Methods and Software, pages 1–38, 2020.
  * [21] Dianne P O’Leary, GW Stewart, and James S Vandergraft.  Estimating the largest eigenvalue of a positive definite matrix.  Mathematics of Computation, 33(148):1289–1292, 1979.
  * [22] C. C. Paige and M. A. Saunders.  Towards a Generalized Singular Value Decomposition.  SIAM Journal on Numerical Analysis, 18(3):398–405, June 1981.
  * [23] Li Qiu, Yanxia Zhang, and Chi-Kwong Li.  Unitarily Invariant Metrics on the Grassmann Space.  SIAM Journal on Matrix Analysis and Applications, 27(2):507–531, January 2005.
  * [24] Y. Saad.  Numerical Methods for Large Eigenvalue Problems.  SIAM, 2nd edition edition, 2011.
  * [25] Yung-Chow Wong.  Sectional curvatures of Grassmann manifolds.  Proceedings of the National Academy of Sciences, 60(1):75–79, May 1968.
  * [26] Hongyi Zhang, Sashank J Reddi, and Suvrit Sra.  Riemannian svrg: Fast stochastic optimization on riemannian manifolds.  Advances in Neural Information Processing Systems, 29, 2016.
  * [27] Hongyi Zhang and Suvrit Sra.  First-order Methods for Geodesically Convex Optimization.  arXiv:1602.06053 [cs, math, stat], February 2016.
  * [28] Hongyi Zhang and Suvrit Sra.  Towards riemannian accelerated gradient methods.  arXiv preprint arXiv:1806.02812, 2018.

##  Appendix A Geodesic convexity

Let $\delta>0$ and thus $\mathcal{V}_{\alpha}$ is the unique minimizer of $f$.
Define the following neighbourhood of $\mathcal{V}_{\alpha}$ in
$\operatorname{Gr}(n,k)$:

$N_{*}(\varphi)=\\{\mathcal{X}\in\operatorname{Gr}(n,k)\colon\theta_{k}(\mathcal{X},\mathcal{V}_{\alpha})<\varphi\\}\qquad\text{with
$\varphi\in[0,\pi/4]$}.$ (33)

Here, $\theta_{k}(\mathcal{X},\mathcal{V}_{\alpha})$ denotes the largest
principle angle between $\mathcal{X}$ and $\mathcal{V}_{\alpha}$. Since
$\theta_{k}$ is a metric on $\operatorname{Gr}(n,k)$ (see [23]), any two
subspaces $\mathcal{X},\mathcal{Y}\in N_{*}(\varphi)$ will satisfy
$\theta_{k}(\mathcal{X},\mathcal{Y})<\pi/2$ by triangle inequality. They thus
have a unique connecting geodesic. It is shown in [3, Lemma 2] that for any
fixed $\varphi\in[0,\pi/4]$ this geodesic remains in $N_{*}(\varphi)$. Each
set $N_{*}(\varphi)$ is thus an open totally geodesically convex set as
defined in, e.g., [8, Def. 11.16].

One of the main results in [3], namely Cor. 4, states that $f$ is geodesically
convex on $N_{*}(\pi/4)$. This is unfortunately wrong and we present a small
counterexample.

#### Counterexample for Cor. 4 in [3].

Here we use the notation of [3]. The reader is encouraged to take a look there
for notational purposes.

Take $c:=\cos(\pi/4)=\sqrt{2}/2$ and $0\leq\varepsilon<1$. Define the matrices

$X_{p}:=\begin{pmatrix}1&0\\\ 0&1\\\ 0&0\\\ 0&0\end{pmatrix},\quad
U_{p}:=\begin{pmatrix}c&0\\\ 0&c\\\ c&0\\\ 0&c\end{pmatrix},\quad
M:=U_{p}\begin{pmatrix}1&0\\\ 0&\varepsilon\end{pmatrix}.$

These matrices satisfy the conditions posed in [3]:

  * •

Principal alignment: $X_{p}^{T}U_{p}=\begin{pmatrix}c&0\\\ 0&c\end{pmatrix}$.

  * •

Principal angles between $X_{p}$ and $U_{p}$ are in $[0,\pi/4]$.

  * •

$U=U_{p}$ since $Q=I$.

Now consider the following tangent vector of unit Frobenius norm:

$\Delta=\begin{pmatrix}0&0\\\ 0&0\\\ 0&1\\\ 0&0\end{pmatrix}.$

It is clearly a tangent vector of $[X_{p}]$ since $X_{p}^{T}\Delta=0$. The
Hessian of $f_{full}$ at $[X_{p}]$ in the direction of $\Delta$ satisfies (see
equation (4.2) in [3])

$\textnormal{Hess}f_{full}([X_{p}])[\Delta,\Delta]=-2\operatorname{Tr}(M^{T}\Delta\Delta^{T}(I-X_{p}X_{p}^{T})M)+\|(\Delta
X_{p}^{T}+X_{p}\Delta^{T})M\|_{F}^{2}.$

Simple calculation shows that

$\textnormal{Hess}f_{full}([X_{p}])[\Delta,\Delta]=-2c^{2}+(1+\varepsilon^{2})c^{2}.$

Hence for $\varepsilon<1$, we have
$\textnormal{Hess}f_{full}([X_{p}])[\Delta,\Delta]<0$ and the $f_{full}$ is
non-convex which is in contrast with Corollary 4. ∎

Instead, our Theorem 18 guarantees convexity when $\varphi$ depends on the
spectral gap. Since $f$ is smooth, the function is geodesically convex on
$N_{*}(\varphi)$ if and only if its Riemannian Hessian is positive definite on
$N_{*}(\varphi)$; see, e.g., [8, Thm. 11.23]. We will therefore compute the
eigenvalues of $\operatorname{Hess}f$ based on its matrix representation. This
requires us to first vectorize the tangent space.

From (3), a matrix $G$ is a tangent vector if and only if $G^{T}X=0$. Hence,
taking $X_{\perp}\in\mathbb{R}^{n\times(n-k)}$ orthonormal such that
$\mathcal{X}^{\perp}=\operatorname{span}(X_{\perp})$, we have the equivalent
definition

$T_{X}\operatorname{Gr}(n,k)=\\{X_{\perp}M\colon M\in\mathbb{R}^{(n-k)\times
k}\\}.$

The matrix $M$ above can be seen as the coordinates of $G=X_{\perp}M$ in the
basis $X_{\perp}$. More specifically, by using the linear isomorphism
$\operatorname{vec}\colon\mathbb{R}^{n\times k}\to\mathbb{R}^{nk}$ that stacks
all columns of a matrix under each other, we can define the tangent vectors of
$\operatorname{Gr}(n,k)$ as standard (column) vectors in the following way:

$\operatorname{vec}(G)=\operatorname{vec}(X_{\perp}M)=(I_{k}\otimes
X_{\perp})\operatorname{vec}(M).$

Here, the Kronecker product $\otimes$ appears due to [14, Lemma 4.3.1]. By
well-known properties of $\otimes$ (see, e.g., [14, Chap. 4.2]), the matrix
$I_{k}\otimes X_{\perp}$ has orthonormal columns. We have thus obtained an
orthonormal basis for the (vectorized) tangent space. With this setup, we can
now construct the Hessian.

######  Lemma 16.

Let $I_{k}\otimes X_{\perp}$ be the orthonormal basis for the vectorization of
$T_{\mathcal{X}}\operatorname{Gr}(n,k)$. Then the Riemannian Hessian of $f$ at
$\mathcal{X}$ in that basis has the symmetric matrix representation

$H_{X}=2(X^{T}AX\otimes I_{n-k}-I_{k}\otimes X_{\perp}^{T}AX_{\perp}).$ (34)

Furthermore, with $1\leq i\leq k$ and $1\leq j\leq n-k$ its $k(n-k)$
eigenvalues satisfy

$\lambda_{i,j}(H_{X})=2(\lambda_{i}(X^{T}AX)-\lambda_{j}(X_{\perp}^{T}AX_{\perp})).$

###### Proof.

Since $\operatorname{vec}$ is a linear isomorphism, the symmetric matrix
$H_{X}$ satisfies

$\operatorname{Hess}f(X)[X_{\perp}M,X_{\perp}M]=\langle\operatorname{vec}(M),H_{X}\operatorname{vec}(M)\rangle,\qquad\forall
M\in\mathbb{R}^{n\times(n-k)},$

where $\langle\cdot,\cdot\rangle$ is the Euclidean inner product. Define
$m=\operatorname{vec}(M)$. Plugging in the formula (9) for
$\operatorname{Hess}f$, we calculate

$\displaystyle\operatorname{Hess}f(X)[X_{\perp}M,X_{\perp}M]$
$\displaystyle=2\langle X_{\perp}M,X_{\perp}MX^{T}AX-AX_{\perp}M\rangle$
$\displaystyle=2\langle(I\otimes X_{\perp})m,(X^{T}AX\otimes
X_{\perp})m-(I\otimes AX_{\perp})m\rangle$ $\displaystyle=2\langle m,(I\otimes
X_{\perp})^{T}(X^{T}AX\otimes X_{\perp}-I\otimes AX_{\perp})m\rangle$
$\displaystyle=2\langle m,(X^{T}AX\otimes I-I\otimes
X_{\perp}^{T}AX_{\perp})m\rangle$

Here, we used typical calculus rules for the Kronecker product (see, e.g.,
[14, Chap. 4.2]). We recognize the matrix $H_{X}$ directly.

The eigenvalues of (34) can be directly obtained using [14, Thm. 4.4.5]. ∎

Taking $X=V_{\alpha}$ and $X_{\perp}=V_{\beta}$, Lemma 16 shows immediately
that the minimal eigenvalue of $\operatorname{Hess}f(\mathcal{V}_{\alpha})$ is
equal to $2\delta=2(\lambda_{k}-\lambda_{k+1})$. Since $\delta>0$,
$\operatorname{Hess}f$ will remain strictly positive definite in a
neighbourhood of $\mathcal{V}_{\alpha}$ by continuity. To quantify this
neighbourhood, we will connect $\mathcal{V}_{\alpha}$ to an arbitrary
$\mathcal{X}$ using a geodesic and see how this influences the bounds of Lemma
16. This also requires connecting $\mathcal{V}_{\beta}$ to
$\mathcal{X}^{\perp}$. The next lemma shows that both geodesics are closely
related. Recall that $\sin(t\theta)$ and $\cos(t\theta)$ denote diagonal
matrices of size $k\times k$. For convenience, we will denote by $O$ a zero
matrix whose dimensions are clear from the context and is not always square.

######  Lemma 17.

Let $X,Y\in\mathbb{R}^{n\times k}$ be such that $X^{T}X=Y^{T}Y=I_{k}$ with
$k\leq n/2$. Denote the principal angles between $\operatorname{span}(X)$ and
$\operatorname{span}(Y)$ by $\theta_{1}\leq\cdots\leq\theta_{k}$ and assume
that $\theta_{k}<\pi/2$. Choose
$X_{\perp},Y_{\perp}\in\mathbb{R}^{n\times(n-k)}$ such that
$X_{\perp}^{T}X_{\perp}=Y_{\perp}^{T}Y_{\perp}=I_{n-k}$ and
$\operatorname{span}(X_{\perp})=\operatorname{span}(X)^{\perp}$,
$\operatorname{span}(Y_{\perp})=\operatorname{span}(Y)^{\perp}$. Define the
curves

$\displaystyle\gamma(t)$ $\displaystyle\colon[0,1]\to\mathbb{R}^{n\times k},$
$\displaystyle t$ $\displaystyle\mapsto
XV_{1}\cos(t\theta)+X_{\perp}V_{2}\begin{bmatrix}O\\\
\sin(t\theta)\end{bmatrix},$ $\displaystyle\gamma_{\perp}(t)$
$\displaystyle\colon[0,1]\to\mathbb{R}^{n\times(n-k)},$ $\displaystyle t$
$\displaystyle\mapsto X_{\perp}V_{2}\begin{bmatrix}I\\\
&\cos(t\theta)\end{bmatrix}-XV_{1}\begin{bmatrix}O&\sin(t\theta)\end{bmatrix},$

where the orthogonal matrices $V_{1},V_{2}$ are the same as in Lemma 6. Then
$\operatorname{span}(\gamma(t))$ is the connecting geodesic on
$\operatorname{Gr}(n,k)$ from $\operatorname{span}(X)$ to
$\operatorname{span}(Y)$. Likewise, $\operatorname{span}(\gamma_{\perp}(t))$
is a connecting geodesic on $\operatorname{Gr}(n,n-k)$ from
$\operatorname{span}(X_{\perp})$ to $\operatorname{span}(Y_{\perp})$.
Furthermore, $\gamma(t)$ and $\gamma_{\perp}(t)$ are orthonormal matrices for
all $t$.

###### Proof.

Assume $\theta_{1}=\cdots=\theta_{r}=0$, where $r=0$ means that
$\theta_{1}>0$. Like in the proof of Prop. 7, the CS decomposition of $X$ and
$Y$ from Lemma 6 can be written in terms of their principal angles
$\theta_{1},\ldots,\theta_{k}$. Since $\theta_{k}<\pi/2$ and $n\leq k/2$, this
gives after dividing certain block matrices the relations

$\displaystyle Y^{T}X$ $\displaystyle=U_{1}\,\cos(\theta)\,V_{1}^{T},$
$\displaystyle Y^{T}X_{\perp}$
$\displaystyle=U_{1}\begin{bmatrix}O_{k\times(n-2k)}&\sin(\theta)\end{bmatrix}V_{2}^{T}$
$\displaystyle Y_{\perp}^{T}X$
$\displaystyle=U_{2}\begin{bmatrix}O_{(n-2k)\times k}\\\
\sin(\theta)\end{bmatrix}V_{1}^{T},$ $\displaystyle Y_{\perp}^{T}X_{\perp}$
$\displaystyle=U_{2}\begin{bmatrix}-I_{n-2k}\\\
&-\cos(\theta)\end{bmatrix}V_{2}^{T},$

where $U_{1},V_{1}$ and $U_{2},V_{2}$ are orthogonal matrices of size $k\times
k$ and $(n-k)\times(n-k)$, resp.

Denote $\mathcal{X}=\operatorname{span}(X)$ and
$\mathcal{Y}=\operatorname{span}(Y)$. By definition, the connecting geodesic
$\gamma(t)$ is determined by the tangent vector
$\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})$, which can be computed from
(6). To this end, we first need the compact SVD of
$M:=X_{\perp}X_{\perp}^{T}Y(X^{T}Y)^{-1}$. Substituting the results from
above, we get (cfr. (17))

$M=X_{\perp}V_{2}\begin{bmatrix}O_{(n-2k)\times k}\\\
\sin(\theta)\end{bmatrix}U_{1}^{T}U_{1}\,(\cos(\theta))^{-1}\,V_{1}^{T}=X_{\perp}V_{2}\begin{bmatrix}O_{(n-2k)\times
k}\\\ I_{k}\end{bmatrix}\,\tan(\theta)\,V_{1}^{T}.$

Observe that this is a compact SVD. Applying (6), we therefore get

$G:=\operatorname{Log}_{\mathcal{X}}(\mathcal{Y})=U\Sigma V^{T}\quad\text{with
$U=X_{\perp}V_{2}\begin{bmatrix}O\\\ I_{k}\end{bmatrix},\ \Sigma=\theta,\
V=V_{1}$}$

and from (4), the connecting geodesic satisfies

$\operatorname{Exp}_{\mathcal{X}}(tG)=\operatorname{span}(\,XV_{1}\cos(t\theta)+X_{\perp}V_{2}\begin{bmatrix}O\\\
I_{k}\end{bmatrix}\sin(t\theta)\,).$

We have proven the stated formula for $\gamma(t)$. Verifying that
$\gamma(t)^{T}\gamma(t)=I_{k}$ follows from a simple calculation that uses
$\cos^{2}(t\theta)+\sin^{2}(t\theta)=I_{k}$.

Denote $\mathcal{X}^{\perp}=\operatorname{span}(X_{\perp})$ and
$\mathcal{Y}^{\perp}=\operatorname{span}(Y_{\perp})$. To prove
$\gamma_{\perp}(t)$, we proceed similarly by computing
$G^{\perp}:=\operatorname{Log}_{\mathcal{X}^{\perp}}(\mathcal{Y}^{\perp})$,
which requires now the SVD of
$M^{\perp}:=XX^{T}Y_{\perp}(X_{\perp}^{T}Y_{\perp})^{-1}$. Again substituting
the results from the CS decomposition, we get

$\displaystyle M^{\perp}$
$\displaystyle=XV_{1}\begin{bmatrix}O_{k\times(n-2k)}&\sin(\theta)\end{bmatrix}U_{2}^{T}U_{2}\begin{bmatrix}-I_{n-2k}\\\
&-\cos(\theta)\end{bmatrix}^{-1}V_{2}^{T}$
$\displaystyle=XV_{1}\begin{bmatrix}O_{k\times(n-2k)}&-\tan(\theta)\end{bmatrix}V_{2}^{T}$

Since (6) requires a compact SVD with a _square_ $\Sigma$, we rewrite this as

$M^{\perp}=\begin{bmatrix}\widetilde{X}&XV_{1}\end{bmatrix}\begin{bmatrix}O_{(n-2k)\times(n-2k)}\\\
&-\tan(\theta)\end{bmatrix}V_{2}^{T}$

where $\widetilde{X}$ contains $n-2k$ columns that are orthonormal to $X$ (the
final result will not depend on $\widetilde{X}$). Let
$\theta^{\perp}_{1}\leq\cdots\leq\theta^{\perp}_{n-k}$ denote the principal
angles between $\mathcal{X}^{\perp}$ and $\mathcal{Y}^{\perp}$. Up to zero
angles, they are the same as those between $\mathcal{X}$ and $\mathcal{Y}$.
Since $k\leq n/2$, we thus have

$\theta^{\perp}_{1}=\cdots=\theta^{\perp}_{n-2k}=0,\
\theta^{\perp}_{n-2k+1}=\theta_{1},\ldots,\theta^{\perp}_{n-k}=\theta_{k}.$

Applying (6) with these principal angles, we obtain

$G^{\perp}:=\operatorname{Log}_{\mathcal{X}^{\perp}}(\mathcal{Y^{\perp}})=U\Sigma
V^{T}\quad\text{with $U=-\begin{bmatrix}\widetilde{X}&XV_{1}\end{bmatrix},\
\Sigma=\theta^{\perp},\ V=V_{2}$}.$

From (4), the corresponding geodesic satisfies

$\displaystyle\operatorname{Exp}_{\mathcal{X}^{\perp}}(tG^{\perp})$
$\displaystyle=\operatorname{span}(\,X_{\perp}V_{2}\cos(t\theta^{\perp})-\begin{bmatrix}\widetilde{X}&XV_{1}\end{bmatrix}\sin(t\theta^{\perp})\,)$
$\displaystyle=\operatorname{span}(\,X_{\perp}V_{2}\begin{bmatrix}I_{n-2k}\\\
&\cos(t\theta)\end{bmatrix}-\begin{bmatrix}O_{n\times(n-2k)}&XV_{1}\sin(t\theta)\end{bmatrix}\,).$

Rewriting the block matrix, we have proven $\gamma_{\perp}(t)$. Its
orthonormality is again a straightforward verification. ∎

With the previous lemma, we can now investigate the Riemannian Hessian of $f$
near $\mathcal{V}_{\alpha}$ when it is given in the matrix form $H_{X}$ of
Lemma 16. Let $\mathcal{X}=\operatorname{span}(X)\in\operatorname{Gr}(n,k)$
with orthonormal $X$. Its principal angles with $\mathcal{V}_{\alpha}$ are
$\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$. Use the substitutions $X\mapsto
V_{\alpha},Y\mapsto X$ and $X_{\perp}\mapsto V_{\beta},Y_{\perp}\mapsto
X_{\perp}$ in Lemma 17 to define the geodesics $\gamma(t)$ and
$\gamma_{\perp}(t)$ that connect $\mathcal{V}_{\alpha}$ to $\mathcal{X}$, and
$\mathcal{V}_{\beta}$ to $\mathcal{X}^{\perp}$, resp. Denoting

$C:=\cos(\theta),\ S:=\sin(\theta),\ \widetilde{C}:=\begin{bmatrix}I\\\
&C\end{bmatrix},\ \widetilde{S}:=\begin{bmatrix}O\\\ S\end{bmatrix},$

we get the following expressions for the geodesics:

$\gamma(t)=V_{\alpha}V_{1}C+V_{\beta}V_{2}\widetilde{S},\quad\gamma_{\perp}(t)=V_{\beta}V_{2}\widetilde{C}-V_{\alpha}V_{1}\widetilde{S}^{T}.$

Recall that $H_{X}$ is defined using $X^{T}AX$ and $X_{\perp}^{T}AX_{\perp}$.
Since $\gamma(1)=XQ_{1}$ and $\gamma_{\perp}(1)=X^{\perp}Q_{2}$ for some
orthogonal matrices $Q_{1},Q_{2}$, we can write with
$A=V_{\alpha}\Lambda_{\alpha}V_{\alpha}^{T}+V_{\beta}\Lambda_{\beta}V_{\beta}^{T}$
that

$\displaystyle Q_{1}^{T}X^{T}AXQ_{1}$ $\displaystyle=\gamma(1)^{T}A\gamma(1)$
(35)
$\displaystyle=C\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,C+\widetilde{S}^{T}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{S}$
$\displaystyle Q_{2}^{T}X_{\perp}^{T}AX_{\perp}Q_{2}$
$\displaystyle=\gamma_{\perp}(1)^{T}A\gamma_{\perp}(1)$
$\displaystyle=\widetilde{C}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{C}+\widetilde{S}\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,\widetilde{S}^{T}.$

Here we used simplifications like
$V_{\beta}^{T}AV_{\alpha}=V_{\beta}^{T}V_{\alpha}\Lambda_{\alpha}=0$.

A simple bounding of the eigenvalues of the difference of these matrices
results in the main result.

######  Theorem 18.

Let $k\leq n/2$. Define the neighbourhood

$B_{*}=\left\\{\mathcal{X}\in\operatorname{Gr}(n,k)\colon\sin^{2}(\theta_{k}(\mathcal{X},\mathcal{V}_{\alpha}))\leq\frac{\delta}{\lambda_{1}+\lambda_{k}}\right\\},$

then $f$ is geodesically convex on $B_{*}$.

###### Proof.

Our aim is to show that $\lambda_{i,j}(H_{X})$ remains positive given the
bound on $\theta_{k}$. From Lemma 16, we see that

$\lambda_{\min}(H_{X})\geq
0\quad\iff\quad\lambda_{\min}(X^{T}AX)\geq\lambda_{\max}(X_{\perp}^{T}AX_{\perp}).$
(36)

Since $Q_{1},Q_{2}$ are orthogonal in (35), it suffices to find a lower and
upper bound of, resp.,

$\displaystyle\lambda_{\min}(X^{T}AX)$
$\displaystyle=\lambda_{\min}(C\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,C+\widetilde{S}^{T}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{S})$
$\displaystyle\lambda_{\max}(X_{\perp}^{T}AX_{\perp})$
$\displaystyle=\lambda_{\max}(\widetilde{C}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{C}+\widetilde{S}\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,\widetilde{S}^{T}).$

Standard eigenvalue inequalities for symmetric matrices (see, e.g., [15, Cor.
4.3.15]) give

$\displaystyle\lambda_{\min}(X^{T}AX)$
$\displaystyle\geq\lambda_{\min}(C\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,C)+\lambda_{\min}(\widetilde{S}^{T}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{S})$
$\displaystyle\lambda_{\max}(X_{\perp}^{T}AX_{\perp})$
$\displaystyle\leq\lambda_{\max}(\widetilde{C}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{C})+\lambda_{\max}(\widetilde{S}\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,\widetilde{S}^{T}).$

Recall that $\lambda_{1}\geq\cdots\geq\lambda_{n}$ are the eigenvalues of $A$.
Since $\widetilde{S}$ is a tall rectangular matrix, we apply the generalized
version of Ostrowski’s theorem from [13, Thm. 3.2] to each term
above666Observe that the cited theorem orders the eigenvalues inversely to the
convention used in this paper. and obtain

$\displaystyle\lambda_{\min}(C\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,C)$
$\displaystyle\geq\lambda_{\min}(C^{2})\lambda_{\min}(\Lambda_{\alpha})=\cos^{2}(\theta_{k})\lambda_{k}$
$\displaystyle\lambda_{\min}(\widetilde{S}^{T}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{S})$
$\displaystyle\geq\lambda_{\min}(\tilde{S}^{T}\tilde{S})\lambda_{\min}(\Lambda_{\beta})=\sin^{2}(\theta_{1})\lambda_{n},$

since the matrices $V_{1},V_{2}$ are orthogonal and
$\theta_{1}\leq\cdots\leq\theta_{k}<\pi/2$. Adding this gives the lower bound

$\lambda_{\min}(X^{T}AX)\geq\cos^{2}(\theta_{k})\lambda_{k}+\sin^{2}(\theta_{1})\lambda_{n}\geq\cos^{2}(\theta_{k})\lambda_{k}.$
(37)

Likewise, using the block structure of $\widetilde{S}$, we get

$\displaystyle\lambda_{\max}(\widetilde{C}\,(V_{2}^{T}\Lambda_{\beta}V_{2})\,\widetilde{C})$
$\displaystyle\leq\lambda_{\max}(C^{2})\lambda_{\max}(\Lambda_{\beta})=\cos^{2}(\theta_{1})\lambda_{k+1}$
$\displaystyle\lambda_{\max}(\widetilde{S}\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,\widetilde{S}^{T})$
$\displaystyle=\lambda_{\max}(S\,(V_{1}^{T}\Lambda_{\alpha}V_{1})\,S)$
$\displaystyle\leq\lambda_{\max}(S^{2})\lambda_{\max}(\Lambda_{\alpha})=\sin^{2}(\theta_{k})\lambda_{1}$

and thus

$\lambda_{\max}(X_{\perp}^{T}AX_{\perp})\leq\cos^{2}(\theta_{1})\lambda_{k+1}+\sin^{2}(\theta_{k})\lambda_{1}\leq\lambda_{k+1}+\sin^{2}(\theta_{k})\lambda_{1}.$
(38)

The condition (36) is thus satisfied when

$\cos^{2}(\theta_{k})\lambda_{k}=\lambda_{k}-\sin^{2}(\theta_{k})\lambda_{k}\geq\lambda_{k+1}+\sin^{2}(\theta_{k})\lambda_{1},$

which reduces to the bound on $\theta_{k}$ in the statement of the theorem.

It remains to show that $B_{*}$ is an open totally geodesically convex set.
Since $\lambda_{1}\geq\lambda_{k}\geq\lambda_{k+1}\geq 0$, we get

$\frac{\lambda_{k}-\lambda_{k+1}}{\lambda_{1}+\lambda_{k}}\leq\frac{\lambda_{k}}{2\lambda_{k}}=\frac{1}{2}.$

Hence, $B_{*}=N_{*}(\varphi)$ with $\varphi\leq\pi/4$ since
$\sin^{2}(\pi/4)=1/2$. ∎

If $k=1$, the proof above can be simplified.

######  Corollary 19.

Let $k=1$ and define the neighbourhood

$B_{*}=\left\\{\mathcal{X}\in\operatorname{Gr}(n,1)\colon\sin^{2}(\theta_{1}(\mathcal{X},\mathcal{V}_{\alpha}))\leq\frac{\delta}{\delta+\lambda_{1}-\lambda_{n}}\right\\}.$

Then $f$ is geodesically convex on $B_{*}$.

###### Proof.

Since $k=1$, there is no need to simplify the bounds (37) and (38) as was done
above. This gives that $f$ is convex as long as

$\cos^{2}(\theta_{1})\lambda_{1}+\sin^{2}(\theta_{1})\lambda_{n}\geq\cos^{2}(\theta_{1})\lambda_{2}+\sin^{2}(\theta_{1})\lambda_{1}.$

Rewriting leads directly to the stated condition on $\sin^{2}(\theta_{1})$. ∎

Remark that optimizing $f$ on $\operatorname{Gr}(n,1)$ is equivalent to

$\min_{x\in\mathbb{R}^{n}}-x^{T}Ax\qquad\text{s.t.}\quad\|x\|=1,$ (39)

which is the minimization of the Rayleigh quotient problem on the unit sphere
$S^{n-1}=\\{x\in\mathbb{R}^{n}\colon x^{T}x=1\\}$. Cor. 19 can therefore also
be phrased in terms of a geodesically convex region for this problem. Denoting
a unit norm top eigenvector of $A$ by $v_{1}$ and using that
$\sin^{2}\theta_{1}=1-\cos^{2}\theta_{1}$, we get that (39) is geodesically
convex on

$\hat{B}_{*}=\left\\{x\in S^{n-1}\colon(x^{T}v_{1})^{2}\geq
1-\frac{\delta}{\delta+\lambda_{1}-\lambda_{n}}\right\\}.$

This result can now be directly compared to [16, Lemma 7] where the
corresponding region is defined as $(x^{T}v_{1})^{2}\geq
1-\frac{\delta}{\delta+\lambda_{1}}$. This is a stricter condition and our
result is therefore a small improvement.

##  Appendix B Convergence of steepest descent with step $\frac{1}{\gamma}$

We now prove convergence of steepest descent with a more tractable choice of
step size compared to the analysis of the main paper, where it depended on the
weak convexity constant $a(\mathcal{X})$. However, this requires a slightly
better initialization at most $\frac{\pi}{2\sqrt{2}}$ away from the minimizer.

###  B.1 Maximum extent of the iterates

We first prove that, while steepest descent with step size at most
$\frac{1}{\gamma}$ does not guarantee contraction in Riemannian distance to
the global minimizer, the distance after $t$ steps is always at most a
constant factor (independent from $t$) within the initial distance.

######  Proposition 20.

Consider steepest descent applied to $f$ with step-size
$\eta\leq\frac{1}{\gamma}$. If the iterates $\mathcal{X}_{t}$ satisfy
$\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha})<\frac{\pi}{2}$, then they
also satisfy

$\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq
2\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha}).$

###### Proof.

Consider the discrete Lyapunov function

$\mathcal{E}(t)=\frac{1}{\gamma}(f(\mathcal{X}_{t})-f^{*})+\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha}).$

Then

$\mathcal{E}(t+1)-\mathcal{E}(t)=\frac{1}{\gamma}(f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t}))+\frac{1}{2}(\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})-\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})).$

Recall that
$\mathcal{X}_{t+1}=\operatorname{Exp}_{\mathcal{X}_{t}}(-\eta\operatorname{grad}f(\mathcal{X}_{t}))$.
By $\gamma$-smoothness of $f$ (cfr. (10)), we have

$\displaystyle f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t})$
$\displaystyle\leq\langle\textnormal{grad}f(\mathcal{X}_{t}),\textnormal{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1})\rangle+\frac{\gamma}{2}\textnormal{dist}(\mathcal{X}_{t},\mathcal{X}_{t+1})^{2}$
$\displaystyle=\left(-\eta+\frac{\gamma}{2}\eta^{2}\right)\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}.$
(40)

We also know by Proposition 7 that

$\langle\textnormal{grad}f(\mathcal{X}),-\operatorname{Log}_{\mathcal{X}}(\mathcal{V_{\alpha})}\rangle\geq
0,$

for any $\mathcal{X}$ with
$\theta_{k}(\mathcal{X},\mathcal{V}_{\alpha})<\pi/2$.

By the fact that the sectional curvatures of the Grassmann manifold are non-
negative, we have

$\displaystyle\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})$
$\displaystyle\leq\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})+\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{X}_{t})-2\langle\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1}),\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{V_{\alpha})}\rangle$
$\displaystyle=\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})+\eta^{2}\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}+2\eta\langle\textnormal{grad}f(\mathcal{X}_{t}),\operatorname{Log}_{\mathcal{X}_{t}}(\mathcal{V_{\alpha})}\rangle$
$\displaystyle\leq\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})+\eta^{2}\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}.$

From $\eta\leq\frac{1}{\gamma}$, we therefore get

$\displaystyle\mathcal{E}(t+1)-\mathcal{E}(t)$
$\displaystyle\leq\left(-\frac{\eta}{\gamma}+\frac{\eta^{2}}{2}\right)\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}+\frac{\eta^{2}}{2}\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}$
$\displaystyle\leq\left(-\frac{\eta}{\gamma}+\eta^{2}\right)\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}\leq
0.$

Since $\mathcal{E}(t)$ does not increase, we have

$\displaystyle\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$
$\displaystyle\leq\mathcal{E}(t)\leq\mathcal{E}(0)=\frac{1}{\gamma}(f(\mathcal{X}_{0})-f^{*})+\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$
$\displaystyle\leq\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})+\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})=\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})$

and the desired result follows. ∎

###  B.2 Convergence under positive eigengap

When $\delta>0$, we can use gradient dominance to prove convergence of
steepest descent to the (unique) minimizer in terms of function values.

######  Proposition 21.

Steepest descent with step-size $\eta=\frac{1}{\gamma}$ initialized at
$\mathcal{X}_{0}$ such that

$\textnormal{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})\leq\frac{\pi}{4}$

satisfies

$f(\mathcal{X}_{t})-f^{*}\leq\left(1-0.32c_{Q}\frac{\delta}{\gamma}\right)^{t}(f(\mathcal{X}_{0})-f^{*}).$

###### Proof.

By the previous result and an induction argument to guarantee that the biggest
angle between $\mathcal{X}_{t}$ and $\mathcal{V}_{\alpha}$ stays strictly less
than $\pi/2$, we can bound the quantities $a(\mathcal{X}_{t})$ uniformly from
below:
Indeed, since
$\textnormal{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\sqrt{2}\cdot\textnormal{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})\leq\frac{\sqrt{2}\pi}{4}$,
we have

$a(\mathcal{X}_{t})\geq\cos(\theta_{k}(\mathcal{X}_{t},\mathcal{V}_{\alpha}))\geq\cos(\textnormal{dist}(\mathcal{X}_{t},\mathcal{V}_{\alpha}))\geq\cos\left(\frac{\sqrt{2}\pi}{4}\right)\geq
0.4.$

Since the step size $\eta=\frac{1}{\gamma}$, the bound (40) implies

$f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t})\leq-\frac{\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}}{2\gamma}$

Applying gradient dominance (Proposition 9), we therefore obtain

$f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t})\leq-\frac{2c_{Q}\delta
a^{2}(\mathcal{X}_{t})}{\gamma}(f(\mathcal{X}_{t})-f^{*})$

and thus

$f(\mathcal{X}_{t+1})-f^{*}\leq\left(1-2c_{Q}a^{2}(\mathcal{X}_{t})\frac{\delta}{\gamma}\right)(f(\mathcal{X}_{t})-f^{*})\leq\left(1-0.32c_{Q}\frac{\delta}{\gamma}\right)(f(\mathcal{X}_{t})-f^{*}).$

By induction the desired result follows. ∎

We now state the iteration complexity of steepest descent algorithm:

######  Theorem 22.

Steepest descent with step-size $\frac{1}{\gamma}$ starting from a subspace
$\mathcal{X}_{0}$ with Riemannian distance at most $\frac{\pi}{4}$ from
$\mathcal{V}_{\alpha}$ computes an estimate $\mathcal{X}_{T}$ of
$\mathcal{V}_{\alpha}$ such that
$\textnormal{dist}(\mathcal{X}_{T},\mathcal{V}_{\alpha})\leq\epsilon$ in at
most

$T=\mathcal{O}\left(\frac{\gamma}{\delta}\log\frac{f(\mathcal{X}_{0})-f^{*}}{\delta\epsilon}\right).$

###### Proof.

For $\textnormal{dist}(\mathcal{X}_{T},\mathcal{V}_{\alpha})<\epsilon$, it
suffices to have

$f(\mathcal{X}_{T})-f^{*}\leq c_{Q}\epsilon^{2}\delta$

by quadratic growth of $f$ in Proposition 5. Using $(1-c)^{T}\leq\exp(-cT)$
for all $T\geq 0$ and $0\leq c\leq 1$, the previous result gives that it
suffices to choose $T$ as the smallest integer such that

$f(\mathcal{X}_{T})-f^{*}\leq\exp\left(-0.32c_{Q}\frac{\delta}{\gamma}T\right)(f(\mathcal{X}_{0})-f^{*})\leq
c_{Q}\epsilon^{2}\delta.$

Solving for $T$ and substituting $c_{Q}=4/\pi^{2}$, we get the required
result. ∎

######  Remark.

The step size in the above theorems satisfies $\eta=\frac{1}{\gamma}$, which
may seem unrealistic. Since an overestimation for $\gamma$ is always a valid
(but less tight) smoothness constant, the previous theorems can also be
phrased for a step size $\eta\leq\frac{1}{\gamma}$.

###  B.3 Gap-less result

We also prove a convergence result for the function values when $\delta$ is
unknown and can be, in particular, equal to $0$.

######  Theorem 23.

Steepest descent with step-size $\eta=\frac{1}{\gamma}$ initialized at
$\mathcal{X}_{0}$ such that

$\textnormal{dist}(\mathcal{X}_{0},\mathcal{V}_{\alpha})\leq\frac{\pi}{4}$

satisfies

$f(\mathcal{X}_{t})-f^{*}\leq\frac{f(\mathcal{X}_{0})-f^{*}+\frac{\gamma}{2}\textnormal{dist}^{2}(\mathcal{X}_{0},\mathcal{V}_{\alpha})}{0.4t+1}=\mathcal{O}\left(\frac{1}{t}\right).$

###### Proof.

By Proposition 20, we have that
$\textnormal{dist}(X_{t},\mathcal{V}_{\alpha})\leq\frac{\sqrt{2}\pi}{4}$ and
$f$ satisfies the weak-quasi-convexity inequality at any iterate
$\mathcal{X}_{t}$ of steepest descent with constant $C_{0}:=0.4$.

Consider the discrete Lyapunov function

$\mathcal{E}(t)=\frac{C_{0}t+1}{\gamma}(f(\mathcal{X}_{t})-f^{*})+\frac{1}{2}\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$

We have that

$\displaystyle\mathcal{E}(t+1)-\mathcal{E}(t)=$
$\displaystyle\frac{C_{0}t+C_{0}+1}{\gamma}(f(\mathcal{X}_{t+1})-f^{*})-\frac{C_{0}t+1}{\gamma}(f(\mathcal{X}_{t})-f^{*})$
$\displaystyle+\frac{1}{2}(\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})-\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})).$

Now we have to estimate a bound for
$\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})-\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})$.
By $\gamma$-smoothness of $f$ and denoting
$\Delta_{t}=f(\mathcal{X}_{t})-f^{*}$ we have

$\Delta_{t+1}-\Delta_{t}\leq\langle\textnormal{grad}f(\mathcal{X}_{t}),\textnormal{Log}_{\mathcal{X}_{t}}(\mathcal{X}_{t+1})\rangle+\frac{\gamma}{2}\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{X}_{t+1})=-\frac{\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}}{2\gamma}$

By $C_{0}$-weak-strong-convexity of $f$ and the fact that the Grassmann
manifold is of positive curvature, we have

$C_{0}\Delta_{t}\leq\frac{\gamma}{2}(\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})-\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha}))+\frac{\|\textnormal{grad}f(\mathcal{X}_{t})\|^{2}}{2\gamma}$

Summing this to the previous inequality, we get

$\textnormal{dist}^{2}(\mathcal{X}_{t+1},\mathcal{V}_{\alpha})-\textnormal{dist}^{2}(\mathcal{X}_{t},\mathcal{V}_{\alpha})\leq\frac{2}{\gamma}((1-C_{0})(f(\mathcal{X}_{t})-f(\mathcal{X}_{t+1}))-C_{0}(f(\mathcal{X}_{t+1})-f^{*})).$

Thus

$\displaystyle\mathcal{E}(t+1)-\mathcal{E}(t)$
$\displaystyle\leq\frac{C_{0}t+1}{\gamma}(f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t}))+\frac{C_{0}}{\gamma}(f(\mathcal{X}_{t+1})-f^{*})$
$\displaystyle+\frac{1-C_{0}}{\gamma}(f(\mathcal{X}_{t})-f(\mathcal{X}_{t+1}))-\frac{C_{0}}{\gamma}(f(\mathcal{X}_{t+1})-f^{*})$
$\displaystyle=\frac{C_{0}t+C_{0}}{\gamma}(f(\mathcal{X}_{t+1})-f(\mathcal{X}_{t}))\leq
0.$

Thus $\mathcal{E}(t)\leq\mathcal{E}(0)$ and the result follows. ∎