remarkRemark AssumptionAssumption Riemannian Acceleration with
PreconditioningNian Shao and Wenbin Chen

# Riemannian Acceleration with Preconditioning for symmetric eigenvalue
problems ††thanks: . Chen is supported by the National Natural Science
Foundation of China (NSFC) 12241101 and 12071090.

Nian Shao Corresponding author. Institute of Mathematics, EPF Lausanne, 1015
Lausanne, Switzerland (nian.shao@epfl.ch). Part of this work was performed
when Shao was at School of Mathematical Sciences, Fudan University. Wenbin
Chen School of Mathematical Sciences and Shanghai Key Laboratory for
Contemporary Applied Mathematics, Fudan University, Shanghai, 200433, People’s
Republic of China (wbchen@fudan.edu.cn).

###### Abstract

Nesterov accelerated gradient method is popular in convex optimization due to
its acceleration property. Zhang and Sra [COLT, 75 (2018), pp. 1703-1723]
generalized it into Riemannian manifold, and proposed a RAGD method for
geodesically convex problems with acceleration. Similar to Euclidean space,
preconditioning techniques are also powerful in Riemannian manifold when
calculating large scale problems for fast convergence. In this paper, we
propose a Riemannian Acceleration with Preconditioning (RAP) for symmetric
eigenvalue problems, which is one of the most important geodesically convex
optimization problem on Riemannian manifold, and obtain the acceleration.
Firstly, the preconditioning for symmetric eigenvalue problems from the
Riemannian manifold viewpoint is discussed. In order to obtain the local
geodesic convexity, we develop the leading angle to measure the quality of the
preconditioner for symmetric eigenvalue problems. A new Riemannian
acceleration, called Locally Optimal Riemannian Accelerated Gradient (LORAG)
method, is proposed to overcome the local geodesic convexity for symmetric
eigenvalue problems. With similar techniques for RAGD and analysis of local
convex optimization in Euclidean space, we analyze the convergence of LORAG.
Incorporating the local geodesic convexity of symmetric eigenvalue problems
under preconditioning with the LORAG, we propose the Riemannian Acceleration
with Preconditioning (RAP) and prove its acceleration. Additionally, when the
Schwarz preconditioner, especially the overlapping or non-overlapping domain
decomposition method, is applied for elliptic eigenvalue problems, we also
obtain the rate of convergence as $1-C\kappa^{-1/2}$, where $C$ is a constant
independent of the mesh sizes and the eigenvalue gap,
$\kappa=\kappa_{\nu}\lambda_{2}/(\lambda_{2}-\lambda_{1})$, $\kappa_{\nu}$ is
the parameter from the stable decomposition, $\lambda_{1}$ and $\lambda_{2}$
are the smallest two eigenvalues of the elliptic operator. Numerical results
show the power of Riemannian acceleration and preconditioning.

###### keywords:

Riemannian acceleration, preconditioning, symmetric eigenvalue problems.

15A18, 65F08, 65F15, 65N25, 90C25.

##  1 Introduction

In this paper, we consider the smallest eigenvalue problem

(1) $Au_{1}=Mu_{1}\lambda_{1},$

where $A$ and $M$ are both symmetric positive definite matrices, $\lambda_{1}$
is the smallest eigenvalue and $u_{1}$ is the corresponding $M$-unit
eigenvector. Assume $\lambda_{1}$ is simple. Let

(2) $f(x)=\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Mx}$

be the Rayleigh quotient, it is well-known that the eigensystem can be solved
by

(3) $u_{1}=\mathop{\arg\min}_{\mathbf{0}\neq
x\in\mathbb{R}^{n}}f(x)\quad\text{and}\quad\lambda_{1}=f(u_{1}).$

From the Riemannian manifold viewpoint, a classical approach [1] for solving
Eq. 3 is

(4) $\displaystyle\min$ $\displaystyle x^{\mathsf{T}}Ax$
$\displaystyle\mathrm{s.t.}$ $\displaystyle x^{\mathsf{T}}Mx=1.$

#### Symmetric eigenvalue problems as convex optimization

The symmetric eigenvalue problem is well-studied in numerical analysis and
scientific computing [8, 43, 53, 44]. The convexity of symmetric eigenvalue
problems is a very important issue for understanding the convergence of
eigensolvers. Actually, a weak version of the convexity, i.e., Polyak-
Lojasiewicz condition, is a textbook-level result in [15, Theorem 5.5]. There
are also some results about the geodesic convexity of the symmetric eigenvalue
problem since it can be viewed as an optimization problem with orthogonality
constraints [1, 19, 20, 14], such as the systematical study of the geodesic
convexity for symmetric eigenvalue problems by Alimisis and Vandereycken [5].
Another convexity structure about symmetric eigenvalue problems is the
implicit convexity proposed by Shao, Chen and Bai in [46]. Besides the
implicit convexity, a symmetric Eigensolver based on Preconditioning and
Implicit Convexity (EPIC) is also proposed.

#### Preconditioned eigensolver

Preconditioned eigensolvers have been studied for more than forty years. For
early developments of this topic, we refer readers to D’yakonov’s book [18]
and Knyazev’s survey paper [23]. One of the most popular preconditioned
eigensolvers is the Locally Optimal Block Preconditioned Conjugate Gradient
(LOBPCG) method [24], but its convergence analysis is still a challenging
problem now. Previous studies about the theoretical analysis of the LOBPCG
method are mainly from two viewpoints. One is the sharp estimation of the
convergence rate for the Preconditioned Steepest Descent (PSD) method [25, 37,
36, 6, 26]. These work focuses on the analysis of the preconditioner, but the
momentum term is not involved. The other one viewpoint is the understanding of
the momentum term [39, 40]. Ovtchinnikov started from the Jacobi orthogonal
complement correction operator, and gave an asymptotic convergence analysis of
LOBPCG. But the acceleration property can not be exploited from his asymptotic
rate of convergence.

#### Acceleration on Riemannian optimization

The gradient method, which dates back to Euler and Lagrange, is almost the
earliest method for convex optimization problems. Forty years ago, Nesterov
[34] proposed an Accelerated Gradient (NAG) method, which attracts lots of
practical and theoretical interests, such as [9, 10, 35, 41]. In 2014, a
dynamical system analogy for the NAG method was proposed by Su, Boyd and
Candes in [49]. With such a new analogy, the NAG method was further studied in
the past ten years [48, 47, 32, 33, 29, 42]. For more results about Nesterov
method, we refer to the recent education paper [51]. When it comes to the
Riemannian optimization [1], there are relative few results. The earliest
research about accelerated method on Riemannian manifold was [28], and the
first computationally tractable Riemannian acceleration method was proposed by
Zhang and Sra in [54]. They analyzed the convergence by generalizing the
estimate sequence technique from the traditional analysis of NAG flow.
Recently, lots of Riemannian accelerated gradient methods were proposed for
the Riemannian optimization problems, such as [2, 3, 17, 21, 4, 30].

#### Contributions

In section 2, we study the geodesic convexity of preconditioned symmetric
eigenvalue problems. A new framework of preconditioning on Riemannian manifold
for symmetric eigenvalue problems is proposed. An interesting property is that
our new preconditioning framework on manifold is similar to the traditional
preconditioned eigensolver in Euclidean space, which is different from the
Riemannian preconditioning by Mishra and Sepulchre in [31]. For better
illustration of the preconditioning for eigenvalue problems, some new
measurements, such as the leading angle, are proposed. And as the final result
of this section, we give the estimation of the local geodesic convexity and
Lipschitz smoothness of symmetric eigenvalue problems under preconditioning.
In section 3, we propose a Locally Optimal Riemannian Accelerated Gradient
(LORAG) method for general Riemannian convex optimization problems with only
local geodesic convexity and Lipschitz smoothness, and prove the acceleration.
The Riemannian Acceleration with Preconditioning (RAP) is proposed by applying
the LORAG method for the preconditioned eigensystems, and the practical
implementations are also discussed. For the theoretical aspect, when the
initial vector $x_{0}$ is sufficiently close to the eigenvector $u_{1}$, we
obtain the asymptotic rate of convergence as

$f(x_{m})-\lambda_{1}\leq
2\Bigl{(}1-\frac{1}{2\sqrt{\kappa}}\Bigr{)}^{m}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)},$

where $f(x)$ is the Rayleigh quotient, $\lambda_{1}$ is the smallest
eigenvalue,

$\kappa\approx\frac{\nu_{\max}}{\nu_{\min}}\frac{1-\lambda_{1}/\lambda_{n}}{1-\lambda_{1}/\lambda_{2}},$

$\nu_{\min}$ and $\nu_{\max}$ are the smallest and largest eigenvalues of the
matrix pencil $(A,B)$, respectively, and $B$ is the preconditioner for $A$,
$\lambda_{2}$ and $\lambda_{n}$ are the second smallest and the largest
eigenvalues of $(A,M)$, respectively. Compared with previous results on
preconditioned eigensolvers, the acceleration property of RAP is proved. For
better understanding our theoretical results, we consider a specific example,
i.e., solving elliptic eigenvalue problems with Schwarz preconditioners, in
section 4. When the overlapping or non-overlapping domain decomposition method
is applied, with sufficiently small coarse mesh size $H$, the asymptotic rate
of convergence as

$f(x_{m})-\lambda_{1}\leq
2\biggl{(}1-C_{\text{RAP}}\Bigl{(}\frac{\lambda_{2}-\lambda_{1}}{\kappa_{\nu}\lambda_{2}}\Bigr{)}^{1/2}\biggr{)}^{m}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)},$

is proved, where $C_{\text{RAP}}$ is a constant independent of the coarse and
fine mesh sizes and the eigenvalues gaps. Compared with the rate of
convergence for the PSD method as

$f(x_{m})-\lambda_{1}\leq\Bigl{(}1-C_{\text{PSD}}\frac{\lambda_{2}-\lambda_{1}}{\kappa_{\nu}\lambda_{2}}\Bigr{)}\bigl{(}f(x_{m-1})-\lambda_{1}\bigr{)},$

the acceleration is clear. Numerical experiments are presented in section 5,
which illustrate the power of Riemannian acceleration and preconditioning.

#### Notations

In order to avoid the proliferation of constants, let $\alpha\lesssim\beta$
represent the statement that $\alpha\leq c\beta$, where the constant $c$ is
positive and independent of the variables in $\alpha$ and $\beta$. The
notation $\alpha\gtrsim\beta$ means $\beta\lesssim\alpha$, and
$\alpha\approx\beta$ means that $\alpha\lesssim\beta$ and
$\beta\lesssim\alpha$. Moreover, for a symmetric positive definite matrix $A$,
the notations $(\cdot,\cdot)_{A}$ and $\lVert\cdot\rVert_{A}$ represent the
corresponding inner product and norm, respectively. For standard Euclidean
inner product and norm, the subscripts are omitted. For a symmetric matrix
pencil, denote $\lambda_{\max}(A,B)$ and $\lambda_{\min}(A,B)$ as the largest
and smallest generalized eigenvalues of $(A,B)$, i.e., $Ax=\lambda Bx$,
respectively.

##  2 Geodesic convexity of preconditioned symmetric eigensystems

###  2.1 New framework for preconditioned eigensystems from Riemannian
manifold viewpoint

Consider the following Riemannian optimization problem:

(5) $\displaystyle\min$ $\displaystyle f(x)$ $\displaystyle\mathrm{s.t.}$
$\displaystyle x^{\mathsf{T}}Bx=1,$

where $B$ is a symmetric positive definite preconditioner for $A$. The
following lemma explains the connection between Eq. 5 and preconditioning.

######  Lemma 2.1.

Consider the following preconditioned gradient flow:

(6) $\left\\{\begin{aligned}
\frac{\,\mathrm{d}x}{\,\mathrm{d}t}&=-B^{-1}\nabla f(x),\\\
x(0)&=x_{0},\end{aligned}\right.$

where $f(x)$ is the Rayleigh quotient defined in Eq. 2. Suppose $\lVert
x_{0}\rVert_{B}=1$, then the gradient flow $x(t)$ from Eq. 6 lies on the
$B$-sphere for all $t>0$, i.e., $\lVert x(t)\rVert_{B}=1$. When
$f(x_{0})<\lambda_{2}$, where $\lambda_{2}$ is the second smallest eigenvalue
of $(A,M)$, then

$f\bigl{(}x(t)\bigr{)}-\lambda_{1}\leq
e^{-Ct}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)},$

where

$C=\frac{4\lambda_{1}\bigl{(}\lambda_{2}-f(x_{0})\bigr{)}}{\kappa_{\nu}\lambda_{2}}>0\quad\text{and}\quad\kappa_{\nu}=\frac{\lambda_{\max}(A,B)}{\lambda_{\min}(A,B)}.$

######  Proof 2.2.

From the definition of the Rayleigh quotient $f(x)$, we know

(7) $\nabla f(x)=\frac{2}{x^{\mathsf{T}}Mx}\bigl{(}Ax-f(x)Mx\bigr{)}.$

Then we have the orthogonality of $\nabla f(x)$ and $x$ in Euclidean space:

$\langle\nabla f(x),x\rangle=0,$

and by (6),

$\frac{\,\mathrm{d}\lVert
x\rVert_{B}^{2}}{\,\mathrm{d}t}=2\Bigl{\langle}\frac{\,\mathrm{d}x}{\,\mathrm{d}t},x\Bigr{\rangle}_{B}=-2\langle\nabla
f(x),x\rangle=0,$

which is the conservation law of $\lVert x(t)\rVert_{B}$. According to the
energy declination of the gradient system as

$\frac{\,\mathrm{d}f}{\,\mathrm{d}t}=\Bigl{\langle}\nabla
f(x),\frac{\,\mathrm{d}x}{\,\mathrm{d}t}\Bigr{\rangle}=-\lVert\nabla
f(x)\rVert_{B^{-1}}^{2}\leq 0,$

we know that $f\bigl{(}x(t)\bigr{)}$ decreases monotonically. Let $\rho=f(x)$
and $r=Ax-\rho Mx$. By the Courant-Fischer minimax theorem, we have

$\frac{\lVert r\rVert_{B^{-1}}^{2}}{\lVert
r\rVert_{A^{-1}}^{2}}\geq\lambda_{\min}(B^{-1},A^{-1})=\lambda_{\min}(A,B)\quad\text{and}\quad\frac{\lVert
x\rVert_{M}^{2}}{\lVert x\rVert_{B}^{2}}=\frac{\lVert
x\rVert_{A}^{2}}{\rho\lVert
x\rVert_{B}^{2}}\leq\frac{\lambda_{\max}(A,B)}{\rho}.$

Using the expression of the $\nabla f$ in Eq. 7, we obtain

$\lVert\nabla f(x)\rVert_{B^{-1}}^{2}=\frac{4\lVert
r\rVert_{B^{-1}}^{2}}{\lVert
x\rVert_{M}^{4}}\geq\frac{4\lambda_{\min}(A,B)\rho\lVert
r\rVert_{A^{-1}}^{2}}{\lambda_{\max}(A,B)\lVert x\rVert_{B}^{2}\lVert
x\rVert_{M}^{2}}=\frac{4\rho\lVert r\rVert_{A^{-1}}^{2}}{\kappa_{\nu}\lVert
x\rVert_{M}^{2}}.$

According to [18, Lemma 3 in Chapter 9.3.1],

$\frac{\lVert r\rVert_{A^{-1}}^{2}}{\lVert
x\rVert_{M}^{2}}\geq\frac{\rho(\rho-\lambda_{1})(\lambda_{2}-\rho)}{\lambda_{1}\lambda_{2}}.$

Thus, combining these relationships above, we know

$\frac{\,\mathrm{d}f}{\,\mathrm{d}t}\leq-\frac{4\rho^{2}(\rho-\lambda_{1})(\lambda_{2}-\rho)}{\kappa_{\nu}\lambda_{1}\lambda_{2}}\leq-\frac{4\lambda_{1}\bigl{(}\lambda_{2}-f(x_{0})\bigr{)}}{\kappa_{\nu}\lambda_{2}}\bigl{(}f(x)-\lambda_{1}\bigr{)}.$

Taking integral we obtain

$f\bigl{(}x(t)\bigr{)}-\lambda_{1}\leq
e^{-Ct}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)},\quad\text{where
}C=\frac{4\lambda_{1}\bigl{(}\lambda_{2}-f(x_{0})\bigr{)}}{\kappa_{\nu}\lambda_{2}}>0,$

which finishes the proof.

Lemma 2.1 is intuitive, and it illustrates how the preconditioner influences
the gradient flow. Due to the homogeneity of the Rayleigh quotient, i.e.,
$f(x)=f(\alpha x)$ for all scalar $\alpha\neq 0$, the minimum of $f(x)$ in any
$B$-sphere is always $\lambda_{1}$. However, as discussed in [23], a good
choice of $B$ may accelerate the convergence. Here we give an interpretation
of this phenomenon: When a good preconditioner is applied, the new trajectory
may lead to fast convergence.

On the other hand, such a phenomenon also provides a connection between
preconditioning and Riemannian optimization: The iterates should be restricted
on a $B$-sphere. Mathematically, the preconditioned eigensystem should be
written as Eq. 5. Compared with Eq. 4, the optimization problem Eq. 5
sacrifices some computational simplicity from cancelling the denominator of
Rayleigh quotient with the manifold constraints, but it enjoys the fast
convergence by computing on a better sphere.

In the rest of this section, we develop some measurements about the
preconditioned eigensystem Eq. 5 and show its local geodesic convexity and
Lipschitz smoothness. In the next section, we propose a new algorithm to
overcome the local convexity, and prove the acceleration property.

###  2.2 Some measurements about preconditioned eigensystems

First, we list some useful results for further analysis.

######  Lemma 2.3.

Suppose $x\in\mathbb{R}^{n}$ with
$\lambda_{1}\leq\rho<(\lambda_{1}+\lambda_{2})/2$, where $\rho=f(x)$, then:

  * •

For any $x_{1}$ and $x_{2}$ satisfying $\lambda_{1}\leq
f(x_{1}),\,f(x_{2})\leq\rho$,

(8) $\frac{\lvert x_{1}^{\mathsf{T}}Mx_{2}\rvert}{\lVert x_{1}\rVert_{M}\lVert
x_{2}\rVert_{M}}\geq 1-\frac{2(\rho-\lambda_{1})}{\lambda_{2}-\lambda_{1}}>0.$

  * •

Let $u_{1}$ be the smallest eigenvector of $(A,M)$ corresponding to
$\lambda_{1}$,

(9) $\frac{(x^{\mathsf{T}}Mu_{1})^{2}}{\lVert x\rVert_{M}^{2}\lVert
u_{1}\rVert_{M}^{2}}\geq\frac{\lambda_{2}-\rho}{\lambda_{2}-\lambda_{1}}\quad\text{and}\quad\frac{(x^{\mathsf{T}}Au_{1})^{2}}{\lVert
x\rVert_{A}^{2}\lVert
u_{1}\rVert_{A}^{2}}\geq\frac{\rho^{-1}-\lambda_{2}^{-1}}{\lambda_{1}^{-1}-\lambda_{2}^{-1}}.$

  * •

For any $v$ satisfying $v^{\mathsf{T}}Mx=0$,

(10) $\lambda_{1}+\lambda_{2}-\rho\leq f(v)\leq\lambda_{n}.$

  * •

Let $r=Ax-\rho Mx$, then

(11)
$\frac{r^{\mathsf{T}}A^{-1}r}{x^{\mathsf{T}}Mx}\leq\frac{\rho(\rho-\lambda_{1})}{\lambda_{1}}.$

######  Proof 2.4.

Let $u_{i}$ be the $M$-orthogonal eigenvector corresponding to the eigenvalue
$\lambda_{i}$ for $i=1,\dotsc,n$. For Eq. 8, assume
$x_{j}=\sum_{i=1}^{n}c_{i,j}u_{i}$, where $j=1$ and $2$, we have

$\rho\geq f(x_{j})=\sum_{i=1}^{n}c_{i,j}^{2}\lambda_{i}\geq
c_{1,j}^{2}\lambda_{1}+(1-c_{1,j}^{2})\lambda_{2}.$

Combining it with $\rho<(\lambda_{1}+\lambda_{2})/2$, we know

$c_{1,j}^{2}\geq\frac{\lambda_{2}-\rho}{\lambda_{2}-\lambda_{1}}>\frac{1}{2}.$

By the Cauchy-Schwarz inequality,

$\displaystyle\lvert x_{1}^{\mathsf{T}}Mx_{2}\rvert$
$\displaystyle=\Bigl{\lvert}\sum_{i=1}^{n}c_{i,1}c_{i,2}\Bigr{\rvert}\geq\lvert
c_{1,1}c_{1,2}\rvert-\frac{1}{2}\sum_{i=2}^{n}(c_{i,1}^{2}+c_{i,2}^{2})=\frac{(\lvert
c_{1,1}\rvert+\lvert c_{1,2}\rvert)^{2}}{2}-1$ $\displaystyle\geq
1-\frac{2(\rho-\lambda_{1})}{\lambda_{2}-\lambda_{1}}>0,$

which finishes the proof of Eq. 8. Since $f(u_{1})=\lambda_{1}$, Eq. 9 are
direct results of Eq. 8. The estimation Eq. 10 has been proved in [38, Lemma
3.1]. For the last estimation Eq. 11, without loss of generality, suppose
$x=\sum_{i=1}^{n}c_{i}u_{i}$, where $\sum_{i=1}^{n}c_{i}^{2}=1$, then

$\displaystyle\frac{r^{\mathsf{T}}A^{-1}r}{x^{\mathsf{T}}Mx}$
$\displaystyle=\sum_{i=1}^{n}\frac{c_{i}^{2}(\lambda_{i}-\rho)^{2}}{\lambda_{i}}\leq\frac{c_{1}^{2}(\rho-\lambda_{1})^{2}}{\lambda_{1}}+\sum_{i=2}^{n}c_{i}^{2}(\lambda_{i}-\rho)$
$\displaystyle=\frac{c_{1}^{2}\rho(\rho-\lambda_{1})}{\lambda_{1}}\leq\frac{\rho(\rho-\lambda_{1})}{\lambda_{1}},$

where the equation is based on the fact $\sum_{i=1}^{n}c_{i}^{2}=1$ and
$\sum_{i=1}^{n}c_{i}^{2}\lambda_{i}=\rho$.

A general measurement for the preconditioner $B$ is the ratio
$\kappa_{\nu}=\nu_{\max}/\nu_{\min}$, where

(12)
$\nu_{\min}=\min_{x\in\mathbb{R}^{n}}\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Bx}\quad\text{and}\quad\nu_{\max}=\max_{x\in\mathbb{R}^{n}}\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Bx}.$

Let $x=B^{-1}Ay$, then $Ax=Bx\lambda$ can be written as $AB^{-1}Ay=Ay\lambda$,
which means matrix pencils $(A,B)$ and $(AB^{-1}A,A)$ share same eigenvalues.
By the Courant-Fischer minimax theorem, it can be proved that

(13)
$\nu_{\min}\leq\frac{x^{\mathsf{T}}AB^{-1}Ax}{x^{\mathsf{T}}Ax}\leq\nu_{\max}.$

Besides $\kappa_{\nu}$, we introduce a new measurement $\vartheta$, which is
called as leading angle, for the preconditioned eigenvalue problem.

######  Definition 2.5 (Leading Angle).

Suppose $(A,M)$ is a symmetric positive definite matrix pencil with
eigenvalues $0<\lambda_{1}<\lambda_{2}\leq\dotsb\leq\lambda_{n}$, and $B$ is a
symmetric positive definite matrix. Let $f(x)$ be the Rayleigh quotient
defined in Eq. 2 and $\lambda_{1}\leq\rho_{*}\leq\lambda_{2}$ be a scalar.
Define the leading angle $\vartheta$ as

(14) $\vartheta(B,\rho_{*};f)=\inf_{\begin{subarray}{c}v^{\mathsf{T}}Bx=0\\\
f(x)\leq\rho_{*}\end{subarray}}\arccos\Bigl{(}\frac{\lvert
v^{\mathsf{T}}Ax\rvert}{\lVert v\rVert_{A}\lVert x\rVert_{A}}\Bigr{)},$

and we denote $\vartheta=\vartheta(B,\rho_{*};f)$ for the simplicity of
notations.

Different from $\kappa_{\nu}$, $\vartheta$ focuses more on the approximation
property in the neighborhood of $u_{1}$. Definitely, $\vartheta$ can be
bounded by $\kappa_{\nu}$ as the following lemma.

######  Lemma 2.6.

With notations in Definition 2.5, then

$\cos\vartheta\leq\sqrt{\kappa_{\nu}-1}.$

######  Proof 2.7.

By Eq. 13, for any $z\in\mathbb{R}^{n}$,

(15) $\min_{\alpha\in\mathbb{R}}\lVert B^{-1}(A-\alpha
B)z\rVert_{B}=\Bigl{(}z^{\mathsf{T}}AB^{-1}Az-\frac{(z^{\mathsf{T}}Az)^{2}}{z^{\mathsf{T}}Bz}\Bigr{)}^{1/2}\leq\sqrt{\nu_{\max}-\nu_{\min}}\lVert
z\rVert_{A}.$

Then using the fact $x^{\mathsf{T}}Bv=0$ and the Cauchy-Schwarz inequality, we
have

(16) $\displaystyle\lvert x^{\mathsf{T}}Av\rvert$
$\displaystyle=\min_{\alpha\in\mathbb{R}}\bigl{\lvert}\alpha
x^{\mathsf{T}}Bv+x^{\mathsf{T}}(A-\alpha
B)v\bigr{\rvert}=\min_{\alpha\in\mathbb{R}}\bigl{\lvert}x^{\mathsf{T}}(A-\alpha
B)v\bigr{\rvert}$ $\displaystyle\leq\lVert
v\rVert_{B}\min_{\alpha\in\mathbb{R}}\lVert B^{-1}(A-\alpha
B)x\rVert_{B}\leq\Bigl{(}\frac{\nu_{\max}-\nu_{\min}}{\nu_{\min}}\Bigr{)}^{1/2}\lVert
v\rVert_{A}\lVert x\rVert_{A}$ $\displaystyle=\sqrt{\kappa_{\nu}-1}\lVert
v\rVert_{A}\lVert x\rVert_{A},$

which finishes the proof.

Here we briefly explain why we need a new measurement $\vartheta$ instead of
$\kappa_{\nu}$. Let us consider an extreme example, suppose $u_{1}$ is also an
eigenvector of $(A,B)$, let $\vartheta=\vartheta(B,\lambda_{1};f)$, then
$\cos\vartheta=0$ for any $B$. But the bound for $\cos\vartheta$ in Lemma 2.6
is much weaker due to the fact $u_{1}$ is an eigenvector of $(A,B)$ only
provides information about the approximation of $B$ and $A$ in one dimension
$\mathrm{span}\\{u_{1}\\}$ rather the global information, which is measured by
$\kappa_{\nu}$. However, in the smallest eigenvalue problem, what we are
really interested in is the information around $u_{1}$. So we use $\vartheta$
as a measurement for the local property preconditioner $B$.

The leading angle defined in Definition 2.5 is a bit abstract, here we give a
decomposition for it.

######  Theorem 2.8.

Suppose $\lambda_{1}\leq\rho_{*}<(\lambda_{1}+\lambda_{2})/2$, the leading
angle $\vartheta$ defined in Eq. 14 satisfies

(17) $\cos\vartheta\leq\frac{\lVert
B^{-1}Au_{1}-\nu_{1}u_{1}\rVert_{A}}{\nu_{\min}}+\Bigl{(}\frac{\nu_{1}}{\nu_{\min}}+\sqrt{\kappa_{\nu}}\Bigr{)}\Bigl{(}\frac{\lambda_{1}^{-1}-\rho_{*}^{-1}}{\lambda_{1}^{-1}-\lambda_{2}^{-1}}\Bigr{)}^{1/2},$

where $(\lambda_{1},u_{1})$ is the smallest eigenpair of $(A,M)$, $\nu_{\min}$
and $\nu_{\max}$ are the smallest and largest eigenvalue of $(A,B)$, and

(18) $\nu_{1}=\mathop{\arg\min}_{\alpha\in\mathbb{R}}\lVert
B^{-1}Au_{1}-\alpha
u_{1}\rVert_{A}=\frac{u_{1}^{\mathsf{T}}AB^{-1}Au_{1}}{u_{1}^{\mathsf{T}}Au_{1}}.$

######  Proof 2.9.

Without loss of generality, assume $\lVert x\rVert_{M}=1$. Let
$\alpha_{1}=x^{\mathsf{T}}Mu_{1}$. For any $x$ satisfying $f(x)\leq\rho_{*}$
and $v^{\mathsf{T}}Bx=0$, using the Cauchy-Schwarz inequality, we have

(19) $\frac{v^{\mathsf{T}}Ax}{\lVert v\rVert_{A}\lVert
x\rVert_{A}}=\frac{v^{\mathsf{T}}B(B^{-1}Ax-\nu_{1}x)}{\lVert
v\rVert_{A}\lVert x\rVert_{A}}\leq\frac{\lVert v\rVert_{B}\lVert
B^{-1}Ax-\nu_{1}x\rVert_{B}}{\lVert v\rVert_{A}\lVert x\rVert_{A}}.$

Let $u_{1}$ be the $M$-unit eigenvector of $(A,M)$ corresponding to
$\lambda_{1}$, by Eq. 13,

(20) $\displaystyle\lVert B^{-1}Ax-\nu_{1}x\rVert_{B}$ $\displaystyle\leq$
$\displaystyle\lVert
B^{-1}A(x-\alpha_{1}u_{1})\rVert_{B}+\nu_{1}\lVert\alpha_{1}u_{1}-x\rVert_{B}+\lvert\alpha_{1}\rvert\lVert
B^{-1}Au_{1}-\nu_{1}u_{1}\rVert_{B}$ $\displaystyle\leq$
$\displaystyle\frac{\lvert\alpha_{1}\rvert}{\sqrt{\nu_{\min}}}\lVert
B^{-1}Au_{1}-\nu_{1}u_{1}\rVert_{A}+\Bigl{(}\frac{\nu_{1}}{\sqrt{\nu_{\min}}}+\sqrt{\nu_{\max}}\Bigr{)}\lVert\alpha_{1}u_{1}-x\rVert_{A}.$

According to Eq. 9, we know

(21) $\frac{\lVert\alpha_{1}u_{1}-x\rVert_{A}^{2}}{\lVert
x\rVert_{A}^{2}}=1-\frac{\alpha_{1}^{2}\lambda_{1}}{\lVert
x\rVert_{A}^{2}}\leq
1-\frac{\alpha_{1}^{2}\lambda_{1}}{\rho_{*}}\leq\frac{\lambda_{2}(\rho_{*}-\lambda_{1})}{\rho_{*}(\lambda_{2}-\lambda_{1})}=\frac{\lambda_{1}^{-1}-\rho_{*}^{-1}}{\lambda_{1}^{-1}-\lambda_{2}^{-1}}.$

Then the lemma is proved by $\lvert\alpha_{1}\rvert\leq 1$, Eq. 13 and Eqs.
19, 20 and 21.

######  Corollary 2.10.

Following the notations in Theorem 2.8, then

(22) $\cos\vartheta\lesssim\sqrt{\kappa_{\nu}}\varepsilon+\Bigl{(}\frac{\lVert
B^{-1}Au_{1}-\nu_{1}u_{1}\rVert_{A}}{\nu_{\min}}+\frac{\nu_{1}}{\nu_{\min}}\varepsilon\Bigr{)},$

where

(23)
$\varepsilon=\Bigl{(}\frac{\rho_{*}-\lambda_{1}}{\lambda_{2}-\lambda_{1}}\Bigr{)}^{1/2}.$

######  Remark 2.11.

Theorems 2.8 and 2.10 give a clearer illustration about the leading angle.

  1. 1.

The term $\sqrt{\kappa_{\nu}}\varepsilon$ in Eq. 22 measure the global
property of the preconditioner $B$. Compared with Lemma 2.6, the factor
$\varepsilon$ reduces the effect of the ratio $\kappa_{\nu}$ when
$\rho_{*}-\lambda_{1}$ is close to zero.

  2. 2.

For the other term, note that $\nu_{1}$ only depends on the behavior of
$B^{-1}$ in one dimension, i.e., $\mathrm{span}\\{Au_{1}\\}$. As $u_{1}$ is
the eigenvector of $(A,M)$, when $B$ approximate $A$ well in
$\mathrm{span}\\{Au_{1}\\}$, $\nu_{1}$ may be bounded independent of
$\kappa_{\nu}$.

The following theorem introduces a new parameter $\varrho$, which is a
preconditioned generalization of Eq. 10.

######  Theorem 2.12.

Using the notations in Definition 2.5, let

(24) $\varrho(B,\rho_{*};f)=\inf_{\begin{subarray}{c}v^{\mathsf{T}}Bx=0\\\
f(x)\leq\rho_{*}\end{subarray}}f(v).$

Denote $\varrho=\varrho(B,\rho_{*};f),$ then

(25)
$\varrho^{-1}\leq\lambda_{2}^{-1}+2(\lambda_{1}^{-1}-\lambda_{2}^{-1})\cos^{2}\vartheta+2(\lambda_{1}^{-1}-\rho_{*}^{-1}).$

######  Proof 2.13.

Without loss of generality, assume $\lVert x\rVert_{A}=\lVert v\rVert_{A}=1$.
Let $\\{\widetilde{u}_{i}\\}$ be $A$-unit eigenvectors of $(A,M)$, i.e.,
$\widetilde{u}_{i}=u_{i}/\sqrt{\lambda_{i}}$ and
$\lVert\widetilde{u}_{i}\rVert_{A}=1$. Suppose $x$ and $v$ have the
eigendecompositions

$x=\sum_{i=1}^{n}c_{i}\widetilde{u}_{i}\quad\text{and}\quad
v=\sum_{i=1}^{n}\widetilde{c}_{i}\widetilde{u}_{i},\quad\text{where
}\sum_{i=1}^{n}c_{i}^{2}=\sum_{i=1}^{n}\widetilde{c}_{i}^{2}=1.$

According to Definition 2.5, we have
$\lvert\sum_{i=1}^{n}c_{i}\widetilde{c}_{i}\rvert=\lvert
x^{\mathsf{T}}Av\rvert\leq\cos\vartheta$. Let $\lvert c_{1}\rvert=\cos\theta$
and $\lvert\widetilde{c}_{1}\rvert=\cos\widetilde{\theta}$, where $\theta$ and
$\widetilde{\theta}$ are both in $[0,\pi/2]$, we know that

$\cos\theta\cos\widetilde{\theta}-\cos\vartheta\leq\Bigl{\lvert}\sum_{i=2}^{n}c_{i}\widetilde{c}_{i}\Bigr{\rvert}\leq\Bigl{(}\sum_{i=2}^{n}c_{i}^{2}\Bigr{)}^{1/2}\Bigl{(}\sum_{i=2}^{n}\widetilde{c}_{i}^{2}\Bigr{)}^{1/2}=\sin\theta\sin\widetilde{\theta}.$

Rearranging this inequality we have

$\cos\vartheta\geq\cos\theta\cos\widetilde{\theta}-\sin\theta\sin\widetilde{\theta}=\cos(\theta+\widetilde{\theta}),$

which implies $\vartheta\leq\widetilde{\theta}+\theta$. Then
$\lvert\widetilde{c}_{1}\rvert$ is bounded by

$\lvert\widetilde{c}_{1}\rvert=\cos\widetilde{\theta}\leq\cos(\vartheta-\theta)\leq\lvert\cos\vartheta\rvert+\lvert\sin\theta\rvert=\cos\vartheta+\sqrt{1-c_{1}^{2}}.$

According to Eq. 9,

$c_{1}^{2}=(x^{\mathsf{T}}A\widetilde{u}_{1})^{2}=\frac{(x^{\mathsf{T}}A\widetilde{u}_{1})^{2}}{\lVert
x\rVert_{A}^{2}\lVert\widetilde{u}_{1}\rVert_{A}^{2}}=\frac{(x^{\mathsf{T}}Au_{1})^{2}}{\lVert
x\rVert_{A}^{2}\lVert
u_{1}\rVert_{A}^{2}}\geq\frac{\lambda_{1}(\lambda_{2}-\rho_{*})}{\rho_{*}(\lambda_{2}-\lambda_{1})}.$

Combining these two inequalities above, we know

$\widetilde{c}_{1}^{2}\leq 2\cos^{2}\vartheta+2(1-c_{1}^{2})\leq
2\cos^{2}\vartheta+2\Bigl{(}\frac{\lambda_{1}^{-1}-\rho_{*}^{-1}}{\lambda_{1}^{-1}-\lambda_{2}^{-1}}\Bigr{)}.$

By combining this inequality with

$\frac{1}{f(v)}=\frac{v^{\mathsf{T}}Mv}{v^{\mathsf{T}}Av}=\sum_{i=1}^{n}\frac{\widetilde{c}_{i}^{2}}{\lambda_{i}}\leq\frac{\widetilde{c}_{1}^{2}}{\lambda_{1}}+\frac{1-\widetilde{c}_{1}^{2}}{\lambda_{2}}=\lambda_{2}^{-1}+(\lambda_{1}^{-1}-\lambda_{2}^{-1})\widetilde{c}_{1}^{2},$

we prove the theorem.

######  Corollary 2.14.

Suppose $\lambda_{1}\leq\rho_{*}<(\lambda_{1}+\lambda_{2})/2$. For any
symmetric positive definite $B$,

$\varrho(B,\lambda_{1};f)>\lambda_{1}.$

For any $\lambda_{1}\leq\rho_{*}^{(1)}\leq\rho_{*}^{(2)}<\lambda_{2}$,

$\vartheta(B,\rho_{*}^{(1)};f)\geq\vartheta(B,\rho_{*}^{(2)};f)\quad\text{and}\quad\varrho(B,\rho_{*}^{(1)};f)\geq\varrho(B,\rho_{*}^{(2)};f).$

With $\vartheta$, $\varepsilon$ and $\varrho$ defined in Eqs. 24, 14 and 23,
respectively, then

(26)
$\frac{\lambda_{2}-\varrho}{\lambda_{2}-\lambda_{1}}\lesssim\cos^{2}\vartheta+\varepsilon^{2}.$

The relationship $\varrho>\rho_{*}$ holds when

(27)
$\cos^{2}\vartheta\leq\frac{3\rho_{*}^{-1}-2\lambda_{1}^{-1}-\lambda_{2}^{-1}}{2(\lambda_{1}^{-1}-\lambda_{2}^{-1})}.$

######  Proof 2.15.

The first two items come from the definition of $\varrho$ and $\vartheta$
directly. The third and fourth items Eqs. 26 and 27 are from the result in
Theorem 2.12.

######  Remark 2.16.

When some specific preconditioners are applied, such as $B=M$ or $B=A$, by
using Eq. 10,

$\varrho(M,\rho_{*};f)\geq\lambda_{1}+\lambda_{2}-\rho_{*}\quad\text{and}\quad\frac{1}{\varrho(A,\rho_{*};f)}\leq\frac{1}{\lambda_{1}}+\frac{1}{\lambda_{2}}-\frac{1}{\rho_{*}}.$

Comparing with the result in Theorem 2.12, we know it is sharp up to a
constant $2$.

######  Theorem 2.17.

Suppose $\lambda_{1}\leq\rho_{*}<(\lambda_{1}+\lambda_{2})/2$, let

$\sigma=\frac{u_{1}^{\mathsf{T}}Au_{1}}{u_{1}^{\mathsf{T}}Bu_{1}}\quad\text{and}\quad\varsigma(B,\rho_{*};f)\equiv\sup_{f(x)\leq\rho_{*}}\Bigl{\lvert}\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Bx}-\sigma\Bigr{\rvert}$

then

$\varsigma(B,\rho_{*};f)\lesssim\frac{\sigma^{2}\varepsilon_{*}^{2}}{\nu_{\min}},$

where
$\varepsilon_{*}=\lambda_{2}\cos\vartheta/(\lambda_{2}-\lambda_{1})+\varepsilon$
with $\vartheta$ and $\varepsilon$ defined in Eqs. 14 and 23 respectively.

######  Proof 2.18.

For any $x$ satisfying $f(x)\leq\rho_{*}$, we first claim that there exist
scalars $\alpha\neq 0$ and $\theta\in[0,\pi/2]$ such that

(28) $\alpha x=\widetilde{u}_{1}\cos\theta+v\sin\theta,$

where $\widetilde{u}_{1}$ is the $A$-unit eigenvector of $(A,M)$ corresponding
to $\lambda_{1}$, $\widetilde{u}_{1}^{\mathsf{T}}Bv=0$,
$\lVert\widetilde{u}_{1}\rVert_{A}=\lVert v\rVert_{A}=1$ and
$0\leq\theta\leq\pi/2$. Consider the $B$-orthogonal decomposition of $x$ as

$x=\widetilde{u}_{1}(\widetilde{u}_{1}^{\mathsf{T}}Bx)+(I-\widetilde{u}_{1}\widetilde{u}_{1}^{\mathsf{T}}B)x,$

where $I$ is the identity operator. Let

$\displaystyle v$
$\displaystyle=\frac{(I-\widetilde{u}_{1}\widetilde{u}_{1}^{\mathsf{T}}B)x}{\lVert(I-\widetilde{u}_{1}\widetilde{u}_{1}^{\mathsf{T}}B)x\rVert_{A}},\quad\theta=\arctan\Bigl{(}\frac{\lVert(I-\widetilde{u}_{1}\widetilde{u}_{1}^{\mathsf{T}}B)x\rVert_{A}}{\widetilde{u}_{1}^{\mathsf{T}}Bx}\Bigr{)},$
$\displaystyle\text{and}\quad\alpha=\bigl{(}(\widetilde{u}_{1}^{\mathsf{T}}Bx)^{2}+\lVert(I-\widetilde{u}_{1}\widetilde{u}_{1}^{\mathsf{T}}B)x\rVert_{A}^{2}\bigr{)}^{-1/2},$

It is easy to verify that the decomposition Eq. 28 holds. Since
$\varsigma(B,\rho_{*};f)$ and the Rayleigh quotient $f(x)$ are homogeneous of
$x$, we can take $\alpha=1$ in Eq. 28 as

(29) $x=\widetilde{u}_{1}\cos\theta+v\sin\theta.$

With $\varrho$ defined in Eq. 24, expanding $x^{\mathsf{T}}Mx$ with Eq. 29 we
obtain

$\frac{\cos^{2}\theta}{\lambda_{1}}+\frac{2\cos\theta\sin\theta\cos\vartheta}{\lambda_{1}}+\frac{\sin^{2}\theta}{\varrho}\geq
x^{\mathsf{T}}Mx\geq\frac{\lVert
x\rVert_{A}^{2}}{\rho_{*}}\geq\frac{1}{\rho_{*}}(1-2\cos\theta\sin\theta\cos\vartheta),$

where the first and last inequalities comes from the definition of $\vartheta$
and $\varrho$, and the second inequality comes from $f(x)\leq\rho_{*}$.
Rearranging the inequalities above with trigonometric formulas, we have

$\rho_{*}^{-1}-\frac{\lambda_{1}^{-1}+\varrho^{-1}}{2}-\frac{(\lambda_{1}^{-1}-\varrho^{-1})\cos
2\theta}{2}\leq(\lambda_{1}^{-1}+\rho_{*}^{-1})\cos\vartheta\sin 2\theta.$

Using the trigonometric formulas, the inequality above and Eq. 26 in order, we
know111Throughout this paper, we assume $\lambda_{2}/\lambda_{1}\lesssim 1$ to
simplify the parameters.

$\displaystyle(\lambda_{1}^{-1}-\lambda_{2}^{-1})\sin^{2}\theta=2\lambda_{1}^{-1}(1-\cos
2\theta)-\lambda_{2}^{-1}\sin^{2}\theta$ $\displaystyle=$
$\displaystyle\rho_{*}^{-1}-\frac{\lambda_{1}^{-1}+\varrho^{-1}}{2}-\frac{(\lambda_{1}^{-1}-\varrho^{-1})\cos
2\theta}{2}+(\lambda_{1}^{-1}-\rho_{*}^{-1})+(\varrho^{-1}-\lambda_{2}^{-1})\sin^{2}\theta$
$\displaystyle\leq$
$\displaystyle(\lambda_{1}^{-1}+\rho_{*}^{-1})\cos\vartheta\sin
2\theta+(\lambda_{1}^{-1}-\rho_{*}^{-1})+(\varrho^{-1}-\lambda_{2}^{-1})\sin^{2}\theta$
$\displaystyle\lesssim$
$\displaystyle(\lambda_{1}^{-1}-\lambda_{2}^{-1})\varepsilon_{*}\sin\theta+(\lambda_{1}^{-1}-\lambda_{2}^{-1})\varepsilon_{*}^{2}.$

Solving the quadratic inequality we obtain
$\sin\theta\lesssim(\lambda_{2}/\lambda_{1})\varepsilon_{*}$. By
$v^{\mathsf{T}}B\widetilde{u}_{1}=0$, we have

$\displaystyle\lvert
x^{\mathsf{T}}Bx-\widetilde{u}_{1}^{\mathsf{T}}B\widetilde{u}_{1}\rvert=\lvert
v^{\mathsf{T}}Bv-\widetilde{u}_{1}^{\mathsf{T}}B\widetilde{u}_{1}\rvert\sin^{2}\theta\lesssim\nu_{\min}^{-1}\varepsilon_{*}^{2},$

where we use facts Eq. 29, $\lVert
v\rVert_{A}=\lVert\widetilde{u}_{1}\rVert_{A}=1$ and
$\nu_{\min}^{-1}=\lambda_{\max}(B,A)$. Similarly,

$\lvert
x^{\mathsf{T}}Ax-\widetilde{u}_{1}^{\mathsf{T}}A\widetilde{u}_{1}\rvert=\lvert
x^{\mathsf{T}}Ax-1\rvert\leq
2\sin\theta\cos\theta\cos\vartheta\lesssim\varepsilon_{*}^{2}.$

Combining these two inequalities above, we obtain

$\Bigl{\lvert}\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Bx}-\frac{{u}_{1}^{\mathsf{T}}A{u}_{1}}{{u}_{1}^{\mathsf{T}}B{u}_{1}}\Bigr{\rvert}=\Bigl{\lvert}\frac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Bx}-\frac{\widetilde{u}_{1}^{\mathsf{T}}A\widetilde{u}_{1}}{\widetilde{u}_{1}^{\mathsf{T}}B\widetilde{u}_{1}}\Bigr{\rvert}\lesssim\frac{\sigma^{2}\varepsilon_{*}^{2}}{\nu_{\min}},$

which finishes the proof.

######  Remark 2.19.

By the Cauchy-Schwarz inequality, it can be proved that $\sigma\leq\nu_{1}$,
where $\nu_{1}$ is defined in Eq. 18, and

$\frac{\sigma+\varsigma}{\sigma}-1\leq\frac{\sigma}{\sigma-\varsigma}-1\leq\frac{\sigma+\varsigma}{\sigma-\varsigma}-1\lesssim\frac{\sigma\varepsilon_{*}^{2}}{\nu_{\min}}\lesssim\frac{\sigma}{\nu_{\min}}\biggl{(}\Bigl{(}\frac{\lambda_{2}}{\lambda_{2}-\lambda_{1}}\Bigr{)}^{2}\cos^{2}\vartheta+\varepsilon^{2}\biggr{)}.$

###  2.3 Instruments on Riemannian manifold

In this part, we briefly review some concepts in Riemannian manifold related
to our analysis, and give some examples on sphere. For more details, we refer
to the literatures about Riemannian manifold in [27, 22] and about Riemannian
optimization in [11, 1].

######  Proposition 2.20.

Let $\mathcal{S}^{n\\!-\\!1}$ be the unit sphere in $\mathbb{R}^{n}$, with the
Riemannian metric induced by the standard Euclidean inner product, the
following properties hold.

  * •

For any $x\in\mathcal{S}^{n\\!-\\!1}$, the tangent space
$T_{x}\mathcal{S}^{n\\!-\\!1}$ is

$T_{x}\mathcal{S}^{n\\!-\\!1}\equiv\\{v\in\mathbb{R}^{n}\mid
x^{\mathsf{T}}v=0\\}.$

And the projection of $\xi\in\mathbb{R}^{n}$ to $T_{x}\mathcal{S}^{n\\!-\\!1}$
is

$P_{x}\xi=(I-xx^{\mathsf{T}})\xi.$

  * •

The geodesic is a curve $\gamma\colon[0,1]\mapsto\mathcal{M}$ of constant
speed and locally minimum length, which plays the role of straight line in
Riemannian manifold. For any $x\in\mathcal{S}^{n\\!-\\!1}$ and $v\in
T_{x}\mathcal{S}^{n\\!-\\!1}$, the geodesic $\gamma$ with initial value $x$
and speed $v$ is

$\gamma(t)=\cos(\lVert v\rVert t)x+\sin(\lVert v\rVert t)\frac{v}{\lVert
v\rVert}.$

  * •

Let $\exp_{x}\colon T_{x}\mathcal{M}\mapsto\mathcal{M}$ be the exponential
map, which is defined as $\exp_{x}(v)=\gamma(1)$, where $\gamma$ is the unique
geodesic satisfying $\gamma(0)=x$ and $\gamma^{\prime}(0)=v$. The exponential
map on sphere is

$\exp_{x}(v)=\gamma(1)=\cos(\lVert v\rVert)x+\sin(\lVert
v\rVert)\frac{v}{\lVert v\rVert}.$

  * •

For a local neighborhood $\mathcal{U}$ of $x$, the exponential map $\exp_{x}$
is a diffeomorphism in $\mathcal{U}$, and its inverse, which is called
logarithmic map, can be defined as
$\log_{x}=\exp_{x}^{-1}\colon\mathcal{U}\mapsto T_{x}\mathcal{M}$. For $x\neq
y\in\mathcal{S}^{n\\!-\\!1}$, the logarithmic map is

$\log_{x}(y)=\arccos(x^{\mathsf{T}}y)\frac{P_{x}y}{\lVert P_{x}y\rVert}.$

Moreover, the norm $\lVert\log_{x}(y)\rVert$ is the distance between $x$ and
$y$ on $\mathcal{M}$.

  * •

For any $x$ and $y$ in $\mathcal{M}$, let $\Gamma_{x}^{y}\colon
T_{x}\mathcal{M}\mapsto T_{y}\mathcal{M}$ be the parallel transport, which
transports tangent vector from $T_{x}\mathcal{M}$ to $T_{y}\mathcal{M}$ along
the geodesic $\gamma$, where $\gamma(0)=x$ and $\gamma(1)=y$. On sphere, let
$v=\log_{x}(y)$, then

$\Gamma_{x}^{y}(u)=\Bigl{(}I+\bigl{(}\cos(t\lVert
v\rVert)-1\bigr{)}\frac{vv^{\mathsf{T}}}{\lVert v\rVert^{2}}-\sin(t\lVert
v\rVert)\frac{xv^{\mathsf{T}}}{\lVert v\rVert}\Bigr{)}u.$

Apart from these basic properties in Proposition 2.20, we also need to handle
the distortion of Riemannian manifold. Suppose $U$ is a $2$-dimensional plane
of $T_{x}\mathcal{M}$, then the sectional curvature of $\mathcal{M}$ at $x$
associated with the plane $U$, which is denoted as $\kappa(x,U)$, is the
Gaussian curvature of $U$ at $x$. For the Riemannian manifold with bounded
sectional curvature, Zhang and Sra also gave an estimation in [54, Theorem 10]
as follows.

######  Proposition 2.21.

For any three points $x$, $y$ and $z\in\mathcal{X}$, where $\mathcal{X}$ is a
geodesically convex set with sectional curvature bounded in $[-K,K]$. Let

$\eta=\sqrt{1+2K\lVert\log_{x}(y)\rVert^{2}},$

then we have

$\displaystyle\lVert\log_{y}(x)-\log_{y}(z)\rVert\leq\eta\lVert\log_{x}(z)\rVert\leq\eta^{2}\lVert\log_{y}(x)-\log_{y}(z)\rVert.$

Since we are studying the convex optimization problem, some background the
convexity and gradients are also needed.

######  Definition 2.22 (Geodesically convex set).

Let $\mathcal{X}\subset\mathcal{M}$. If for any $x,\,y\in\mathcal{X}$, there
is a geodesic $\gamma$ with $\gamma(0)=x$ and $\gamma(1)=y$ such that
$\gamma(t)\in\mathcal{X}$ for any $0\leq t\leq 1$, then $\mathcal{X}$ is
called geodesically convex.

######  Definition 2.23 (Geodesically convex function).

Let $f\colon\mathcal{X}\subset\mathcal{M}$, where $\mathcal{X}$ is a
geodesically convex set. If for any $x,\,y\in\mathcal{X}$, there is a geodesic
$\gamma$ with $\gamma(0)=x$ and $\gamma(1)=y$ such that

$f\bigl{(}\gamma(t)\bigr{)}\leq(1-t)f(x)+tf(y)\quad\forall\,0\leq t\leq 1,$

then $f$ is called geodesically convex.

######  Definition 2.24 (Gradient).

Let $f\colon\mathcal{M}\mapsto\mathbb{R}$ be a differentiable function, then
the gradient of $f$ at $x\in\mathcal{M}$, denoted by $\mathrm{grad}f(x)$ is
defined as the unique element of $T_{x}\mathcal{M}$ that satisfies

$\langle\mathrm{grad}f(x),\xi\rangle=\mathrm{D}f(x)[\xi]\quad\forall\,\xi\in
T_{x}\mathcal{M},$

where $\mathrm{D}f(x)[\xi]$ is the directional derivative of $f$ along $\xi$.

######  Proposition 2.25.

If $\log_{x}$ is well-defined for any $x\in\mathcal{X}$ and $f$ is a
differentiable function, an equivalent definition for a geodesically convex
function is that: there exists a constant $\mu>0$ such that

$f(y)\geq
f(x)+\langle\mathrm{grad}f(x),\log_{x}(y)\rangle+\frac{\mu}{2}\lVert\log_{x}(y)\rVert^{2}\quad\forall\,x,\,y\in\mathcal{X}.$

In this case, we call $f$ as $\mu$-geodesically convex.

######  Definition 2.26 (Lipschitz smooth function).

Let $f\colon\mathcal{X}\mapsto\mathbb{R}$ be a differentiable function, assume
$\log_{x}$ is well-defined for any $x\in\mathcal{X}$. If there exists a
constant $L>0$ such that

$f(y)\leq
f(x)+\langle\mathrm{grad}f(x),\log_{x}(y)\rangle+\frac{L}{2}\lVert\log_{x}(y)\rVert^{2}\quad\forall\,x,\,y\in\mathcal{X},$

then $f$ is called $L$-Lipschitz smooth.

######  Proposition 2.27.

If $f$ is $L$-Lipschitz smooth on $\mathcal{X}$, then

$\lVert\mathrm{grad}f(x)-\Gamma_{y}^{x}\mathrm{grad}f(y)\rVert\leq
L\lVert\log_{x}(y)\rVert\quad\forall\,x,\,y\in\mathcal{X}.$

######  Proposition 2.28 (Second order condition).

Let $f\colon\mathcal{X}\mapsto\mathbb{R}$ be a smooth function, assume
$\log_{x}$ is well-defined for any $x\in\mathcal{X}$. If there exist constants
$0<\mu<L$ such that for any $v\in T_{x}\mathcal{X}$,

$\mu\lVert v\rVert^{2}\leq\langle v,\nabla^{2}f(x)[v]\rangle\leq L\lVert
v\rVert^{2}$

holds, where $\nabla^{2}f(x)$ is the Hessian of $f$ at $x$ defined as

$\nabla^{2}f(x)[v]=\nabla_{v}\mathrm{grad}f(x),$

then $f$ is $\mu$-geodesically convex and $L$-Lipschitz smooth.

###  2.4 Geodesic convexity on $M$-sphere

Let us first consider the geodesic convexity of Eq. 4. For simplicity of our
discussion, we use some transformation to transform $M$-sphere to a standard
sphere. Let $z=M^{1/2}x$, then we know that

$z^{\mathsf{T}}z=x^{\mathsf{T}}Mx=1\quad\text{and}\quad
z^{\mathsf{T}}A_{M}z=x^{\mathsf{T}}Ax,$

where $A_{M}=M^{-1/2}AM^{-1/2}$. In this part, we consider

(30) $\min_{z^{\mathsf{T}}z=1}\quad f_{M}(z)$

instead of Eq. 3, where $f_{M}(z)\equiv z^{\mathsf{T}}A_{M}z$. A direct result
from the transformation is that $\\{\lambda_{i}\\}_{i=1}^{n}$ are eigenvalues
of $A_{M}$, and

(31) $u_{M}=\frac{M^{1/2}u_{1}}{\lVert u_{1}\rVert_{M}}$

is the unit eigenvector corresponding to $\lambda_{1}$.

Now we need to define a geodesically convex set $\mathcal{X}_{M}$ for
computation since most iterative eigensolvers are locally convergent.

######  Lemma 2.29.

Suppose $\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$, then
the set

$\mathcal{X}_{M}=\\{x\in\mathcal{S}^{n\\!-\\!1},\,x^{\mathsf{T}}u_{M}>0\mid
x^{\mathsf{T}}A_{M}x\leq\rho_{\mathcal{X}}\\}$

is geodesically convex, where $u_{M}$ is defined in Eq. 31.

######  Proof 2.30.

Since the hemisphere
$\mathcal{S}^{n\\!-\\!1}_{+}=\\{x\in\mathcal{S}^{n\\!-\\!1},\,x^{\mathsf{T}}u_{M}>0\\}$
is geodesically convex, it is sufficient to prove for any $x$ and $y$ in
$\mathcal{X}_{M}$ and the geodesic $\gamma$ satisfying $\gamma(0)=x$ and
$\gamma(1)=y$, the inequality

(32) $\lVert\gamma(t)\rVert_{A_{M}}^{2}\leq\rho_{\mathcal{X}}$

holds for all $0\leq t\leq 1$. Note that for any
$x\in\mathcal{S}^{n\\!-\\!1}$, if $x$ is orthogonal to $u_{M}$, then
$x^{\mathsf{T}}A_{M}x\geq\lambda_{2}>\rho_{\mathcal{X}}$, which means
$\mathcal{X}_{M}$ is a closed subset of the hemisphere
$\mathcal{S}^{n\\!-\\!1}_{+}$. If Eq. 32 does not hold for some $x$ and $y$ in
$\mathcal{X}_{M}$, since $\mathcal{X}_{M}$ is closed, we may assume

(33)
$x^{\mathsf{T}}A_{M}x=\rho_{\mathcal{X}},\quad\gamma(0)=x,\quad\gamma(1)=y,\quad\lVert\gamma(t)\rVert_{A_{M}}^{2}>\rho_{\mathcal{X}}\quad\text{for
all }0<t<1.$

Let $v=\log_{x}(y)$ be the logarithmic map and

$\phi(t)=\lVert\gamma(t)\rVert_{A_{M}}^{2}-\rho_{\mathcal{X}}=\Bigl{(}\frac{v^{\mathsf{T}}A_{M}v}{v^{\mathsf{T}}v}-\rho_{\mathcal{X}}\Bigr{)}\sin^{2}(\lVert
v\rVert t)+\frac{x^{\mathsf{T}}A_{M}v}{\lVert v\rVert}\sin(2\lVert v\rVert
t).$

According to Eq. 33, we have

(34) $\phi^{\prime}(0)=2x^{\mathsf{T}}A_{M}v\geq 0.$

Using Eqs. 10 and 8, when
$\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$, we obtain

(35)
$\frac{v^{\mathsf{T}}A_{M}v}{v^{\mathsf{T}}v}\geq\lambda_{1}+\lambda_{2}-\rho_{\mathcal{X}}\geq\rho_{\mathcal{X}}\quad\text{and}\quad\cos(\lVert
v\rVert)=x^{\mathsf{T}}y\geq
1-\frac{2(\rho_{\mathcal{X}}-\lambda_{1})}{\lambda_{2}-\lambda_{1}}>0.$

Since $\lVert v\rVert$ is the length of geodesic, which is a minor arc, when
$x$ and $y$ are different points, we know $0<\lVert v\rVert<\pi/2$. Combining
Eqs. 34 and 35, we have

$\displaystyle 0$
$\displaystyle\geq\phi(1)=\Bigl{(}\frac{v^{\mathsf{T}}A_{M}v}{v^{\mathsf{T}}v}-\rho_{\mathcal{X}}\Bigr{)}\sin^{2}(\lVert
v\rVert)+\frac{x^{\mathsf{T}}A_{M}v}{\lVert v\rVert}\sin(2\lVert v\rVert)\geq
0,$

where the equalities hold when

$\frac{v^{\mathsf{T}}A_{M}v}{v^{\mathsf{T}}v}=\rho_{\mathcal{X}}\quad\text{and}\quad
x^{\mathsf{T}}A_{M}v=0.$

In this situation, $\phi(t)=0$ for all $t\in[0,1]$, a contradiction to Eq. 33.

Here is the geodesic convexity and Lipschitz smoothness of Eq. 30.

######  Theorem 2.31.

Suppose $\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$, let
$f_{M}(x)=x^{\mathsf{T}}A_{M}x$, then for any $x\in\mathcal{X}_{M}$ and $v\in
T_{x}\mathcal{S}^{n\\!-\\!1}$, the following inequality holds:

$2(\lambda_{1}+\lambda_{2}-2\rho_{\mathcal{X}})\lVert v\rVert^{2}\leq\langle
v,\nabla^{2}f_{M}(x)[v]\rangle\leq 2(\lambda_{n}-\lambda_{1})\lVert
v\rVert^{2}.$

These inequalities imply $f_{M}$ is
$2(\lambda_{1}+\lambda_{2}-2\rho_{\mathcal{X}})$-geodesically convex and
$2(\lambda_{n}-\lambda_{1})$-Lipschitz smooth.

######  Proof 2.32.

Without loss of generality, assume $\lVert v\rVert=1$, we know

$\langle
v,\nabla^{2}f_{M}(x)[v]\rangle=2v^{\mathsf{T}}\bigl{(}A_{M}-f_{M}(x)I\bigr{)}v=2f_{M}(v)-2f_{M}(x).$

Then the lemma is proved by Eq. 10.

######  Remark 2.33.

The first order approximation for the condition number of Eq. 30 is

$\kappa=\frac{L}{\mu}=\frac{\lambda_{n}}{\lambda_{2}}\frac{1-\lambda_{1}/\lambda_{n}}{1-\lambda_{1}/\lambda_{2}}+\mathcal{O}\Bigl{(}\frac{(\lambda_{n}-\lambda_{1})(\rho_{\mathcal{X}}-\lambda_{1})}{(\lambda_{2}-\lambda_{1})^{2}}\Bigr{)}.$

###  2.5 Geodesic convexity on $B$-sphere

Let us turn to the geodesic convexity of Eq. 5. First, some transformations
are also needed. Let

$A_{B}=B^{-1/2}AB^{-1/2}\quad\text{and}\quad M_{B}=B^{-1/2}MB^{-1/2},$

the optimization problem Eq. 5 becomes

(36)
$\min_{x\in\mathcal{S}^{n\\!-\\!1}}f_{B}(x)\equiv\frac{x^{\mathsf{T}}A_{B}x}{x^{\mathsf{T}}M_{B}x},$

and the eigenvalues of $(A_{B},M_{B})$ are $\\{\lambda_{i}\\}_{i=1}^{n}$, and

$u_{B}=\frac{B^{1/2}u_{1}}{\lVert u_{1}\rVert_{B}}$

is the eigenvector corresponding to $\lambda_{1}$. Similar to Lemma 2.29, the
convex set $\mathcal{X}_{B}$ needs to be defined first.

######  Lemma 2.34.

Suppose $\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$, then
the set

$\mathcal{X}_{B}=\\{x\in\mathcal{S}^{n\\!-\\!1},\,x^{\mathsf{T}}u_{B}>0\mid\frac{x^{\mathsf{T}}A_{B}x}{x^{\mathsf{T}}M_{B}x}\leq\rho_{\mathcal{X}}\\}$

is geodesically convex.

######  Proof 2.35.

According to Lemma 2.29, we know

(37) $\displaystyle\mathcal{X}$
$\displaystyle=\\{x^{\mathsf{T}}M_{B}x=1,\,x^{\mathsf{T}}u_{B}>0\mid\frac{x^{\mathsf{T}}A_{B}x}{x^{\mathsf{T}}M_{B}x}\leq\rho_{\mathcal{X}}\\}$
$\displaystyle=\\{x^{\mathsf{T}}M_{B}x=1,\,x^{\mathsf{T}}u_{B}>0\mid
x^{\mathsf{T}}A_{B}x\leq\rho_{\mathcal{X}}\\},$

is geodesically convex. Let $\psi\colon\mathcal{X}_{B}\mapsto\mathcal{X}$ be
the scaling transformation defined as

$\psi(x)=\frac{x}{\lVert x\rVert_{M_{B}}},$

then for any $x_{1}$ and $x_{2}$ in $\mathcal{X}_{B}$, and geodesic
$\gamma_{B}(t)\in\mathcal{X}_{B}$ satisfying $\gamma_{B}(0)=x_{1}$ and
$\gamma_{B}(1)=x_{2}$, there exists a geodesic $\gamma(t)\in\mathcal{X}$ such
that $\gamma(0)=\psi(x_{1})$ and $\gamma(1)=\psi(x_{2})$. It can be proved
that for any $t\in[0,1]$, there exists unique $t_{B}\in[0,1]$ such that
$\psi\bigl{(}\gamma(t)\bigr{)}=\gamma_{B}(t)$. Thus, $\mathcal{X}$ is
geodesically convex leads to $\mathcal{X}_{B}$ is geodesically convex.

Now we can give the convexity of $f_{B}$ on $\mathcal{X}_{B}$ like Theorem
2.31.

######  Theorem 2.36.

Suppose $\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$, let
$f_{B}(x)$ be the Rayleigh quotient defined in Eq. 36, denote

$\varrho_{B}=\varrho(I,\rho_{\mathcal{X}};f_{B}),\quad\vartheta_{B}=\vartheta(I,\rho_{\mathcal{X}};f_{B})\quad\text{and}\quad\varsigma_{B}=\varsigma(I,\rho_{\mathcal{X}};f_{B}),$

under the condition $\varrho_{B}>\rho_{\mathcal{X}}$, the following inequality
holds:

(38) $\mu_{B}\lVert v\rVert^{2}\leq\langle v,\nabla^{2}f_{B}(x)[v]\rangle\leq
L_{B}\lVert v\rVert^{2},$

where

(39) $\displaystyle\mu_{B}$
$\displaystyle=\frac{2\nu_{\min}\lambda_{1}}{\sigma+\varsigma_{B}}\Bigl{(}1-\frac{\rho_{\mathcal{X}}}{\varrho_{B}}\Bigr{)}-C_{\mathcal{X}},\quad
L_{B}=\frac{2\nu_{\max}\rho_{\mathcal{X}}}{\sigma-\varsigma_{B}}\Bigl{(}1-\frac{\lambda_{1}}{\lambda_{n}}\Bigr{)}+C_{\mathcal{X}},$
$\displaystyle C_{\mathcal{X}}$
$\displaystyle=8\kappa_{\nu}\rho_{\mathcal{X}}\biggl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{1}}+\Bigl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{1}}\Bigr{)}^{1/2}\cos\vartheta_{B}\biggr{)},\quad\sigma=\frac{u_{1}^{\mathsf{T}}Au_{1}}{u_{1}^{\mathsf{T}}Bu_{1}},$

$\nu_{\min}$ and $\nu_{\max}$ are the smallest and largest eigenvalues of
$A_{B}$, respectively.

######  Proof 2.37.

Without loss of generality, assume $\lVert v\rVert=1$, then

$\langle v,\nabla^{2}f_{B}(x)[v]\rangle=J_{1}-J_{2},$

where

$J_{1}=\frac{2\lVert v\rVert_{A_{B}}^{2}f_{B}(x)}{\lVert
x\rVert_{A_{B}}^{2}}\Bigl{(}1-\frac{f_{B}(x)}{f_{B}(v)}\Bigr{)}\text{ and
}J_{2}=\frac{8(v^{\mathsf{T}}M_{B}x)}{\lVert
x\rVert_{M_{B}}^{4}}\Bigl{(}v^{\mathsf{T}}\bigl{(}A_{B}x-f_{B}(x)M_{B}x\bigr{)}\Bigr{)}.$

Under the condition $\varrho_{B}>\rho_{\mathcal{X}}$, the first term $J_{1}$
can be bounded with

(40)
$\frac{2\nu_{\min}\lambda_{1}}{\sigma+\varsigma_{B}}\Bigl{(}1-\frac{\rho_{\mathcal{X}}}{\varrho_{B}}\Bigr{)}\leq
J_{1}\leq\frac{2\nu_{\max}\rho_{\mathcal{X}}}{\sigma-\varsigma_{B}}\Bigl{(}1-\frac{\lambda_{1}}{\lambda_{n}}\Bigr{)},$

by $\lambda_{1}\leq f_{B}(x)\leq\rho_{\mathcal{X}}$, $\varrho_{B}\leq
f_{B}(v)\leq\lambda_{n}$, $\nu_{\min}\leq\lVert
v\rVert_{A_{B}}^{2}\leq\nu_{\max}$ and Theorem 2.17. For the second term
$J_{2}$, by the Cauchy-Schwarz inequality, we have

$\bigl{\lvert}v^{\mathsf{T}}\bigl{(}A_{B}x-f_{B}(x)M_{B}x\bigr{)}\bigr{\rvert}\leq\lVert
v\rVert_{A_{B}}\bigl{\lVert}\bigl{(}A_{B}x-f_{B}(x)M_{B}x\bigr{)}\bigr{\rVert}_{A_{B}^{-1}}.$

According to the definition of $\vartheta_{B}$,

$\displaystyle\frac{\lvert v^{\mathsf{T}}M_{B}x\rvert}{\lVert
x\rVert_{M_{B}}^{2}}$ $\displaystyle=\frac{\lvert
v^{\mathsf{T}}f_{B}(x)M_{B}x\rvert}{\lVert
x\rVert_{A_{B}}^{2}}\leq\frac{\bigl{\lvert}v^{\mathsf{T}}\bigl{(}A_{B}x-f_{B}(x)M_{B}x\bigr{)}\bigr{\rvert}+\lvert
v^{\mathsf{T}}A_{B}x\rvert}{\lVert x\rVert_{A_{B}}^{2}}$
$\displaystyle\leq\frac{\lVert
v\rVert_{A_{B}}\bigl{\lVert}\bigl{(}A_{B}x-f_{B}(x)M_{B}x\bigr{)}\bigr{\rVert}_{A_{B}^{-1}}+\lVert
v\rVert_{A_{B}}\lVert x\rVert_{A_{B}}\cos\vartheta_{B}}{\lVert
x\rVert_{A_{B}}^{2}}$

Combining these two inequalities above with $f_{B}(x)\leq\rho_{\mathcal{X}}$
and

$\nu_{\min}\leq\lVert x\rVert_{A_{B}}^{2},\lVert
v\rVert_{A_{B}}^{2}\leq\nu_{\max},$

we have

$\displaystyle\lvert J_{2}\rvert$ $\displaystyle\leq\frac{8\lVert
v\rVert_{A_{B}}^{2}}{\lVert x\rVert_{A_{B}}^{2}}\Biggl{(}\frac{\lVert
A_{B}x-f_{B}(x)M_{B}x\rVert_{A_{B}^{-1}}^{2}}{\lVert
x\rVert_{M_{B}}^{2}}+\frac{\lVert
A_{B}x-f_{B}(x)M_{B}x\rVert_{A_{B}^{-1}}\lVert
x\rVert_{A_{B}}\cos\vartheta_{B}}{\lVert x\rVert_{M_{B}}^{2}}\Biggr{)}$
$\displaystyle\leq\frac{8\nu_{\max}}{\nu_{\min}}\Biggl{(}\frac{\bigl{\lVert}A_{B}x-f_{B}(x)M_{B}x\bigr{\rVert}_{A_{B}^{-1}}^{2}}{\lVert
x\rVert_{M_{B}}^{2}}+\frac{\bigl{\lVert}A_{B}x-f_{B}(x)M_{B}x\bigr{\rVert}_{A_{B}^{-1}}}{\lVert
x\rVert_{M_{B}}}\sqrt{\rho_{\mathcal{X}}}\cos\vartheta_{B}\Biggr{)}.$

By Eq. 11, we know that

$\frac{\bigl{\lVert}A_{B}x-f_{B}(x)M_{B}x\bigr{\rVert}_{A_{B}^{-1}}}{\lVert
x\rVert_{M_{B}}}\leq\Bigl{(}\frac{\rho_{\mathcal{X}}(\rho_{\mathcal{X}}-\lambda_{1})}{\lambda_{1}}\Bigr{)}^{1/2}.$

Combining these two inequalities above with
$\kappa_{\nu}=\nu_{\max}/\nu_{\min}$, we have

(41) $\lvert J_{2}\rvert\leq
8\kappa_{\nu}\rho_{\mathcal{X}}\biggl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{1}}+\Bigl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{1}}\Bigr{)}^{1/2}\cos\vartheta_{B}\biggr{)}.$

Then the lemma is proved by Eqs. 40 and 41.

######  Corollary 2.38.

The condition $\varrho_{B}>\rho_{\mathcal{X}}$ holds when $\varepsilon$
defined in Eq. 23 is sufficiently small due to Eqs. 27 and 26. Moreover, we
have

$\displaystyle\mu_{B}$
$\displaystyle=\frac{2\nu_{\min}\lambda_{1}}{\sigma}\Bigl{(}1-\frac{\lambda_{1}}{\lambda_{2}}\Bigr{)}+\mathcal{O}\biggl{(}\kappa_{\nu}\Bigl{(}(\lambda_{2}-\lambda_{1})\varepsilon^{2}+\frac{\lambda_{2}^{2}}{\lambda_{2}-\lambda_{1}}\cos^{2}\vartheta_{B}\Bigr{)}\biggr{)},$
$\displaystyle L_{B}$
$\displaystyle=\frac{2\nu_{\max}\lambda_{1}}{\sigma}\Bigl{(}1-\frac{\lambda_{1}}{\lambda_{n}}\Bigr{)}+\mathcal{O}\biggl{(}\kappa_{\nu}\lambda_{1}\Bigl{(}\varepsilon^{2}+\frac{\lambda_{2}^{2}}{(\lambda_{2}-\lambda_{1})^{2}}\cos^{2}\vartheta_{B}\Bigr{)}\biggr{)},$

which means $f_{B}$ is $\mu_{B}$-geodesically convex and $L_{B}$-Lipschitz
smooth. Furthermore, the condition number of Eq. 36 is

$\kappa_{B}=\frac{L_{B}}{\mu_{B}}=\frac{\nu_{\max}}{\nu_{\min}}\frac{1-\lambda_{1}/\lambda_{n}}{1-\lambda_{1}/\lambda_{2}}+\mathcal{O}\biggl{(}\frac{\sigma}{\nu_{\min}}\frac{\kappa_{\nu}^{2}\lambda_{2}}{\lambda_{2}-\lambda_{1}}\Bigl{(}\varepsilon^{2}+\frac{\lambda_{2}^{2}}{(\lambda_{2}-\lambda_{1})^{2}}\cos^{2}\vartheta_{B}\Bigr{)}\biggr{)}.$

Compared with the condition number of Eq. 30 in Remark 2.33, i.e.,

$\kappa=\frac{\lambda_{n}}{\lambda_{2}}\frac{1-\lambda_{1}/\lambda_{n}}{1-\lambda_{1}/\lambda_{2}}+\mathcal{O}\Bigl{(}\frac{\lambda_{n}-\lambda_{1}}{\lambda_{2}-\lambda_{1}}\varepsilon^{2}\Bigr{)},$

when $B$ is a good preconditioner for $A$, i.e.,
$\kappa_{\nu}=\nu_{\max}/\nu_{\min}$ and $\cos\vartheta_{B}$ are small, the
condition number $\kappa_{B}$ is significantly improved.

##  3 Locally optimal Riemannian accelerated gradient method

In this section, we consider a general convex Riemannian optimization problem

(42) $x_{*}=\mathop{\arg\min}_{x\subset\mathcal{M}}f(x),$

where $\mathcal{M}$ is a Riemannian manifold and
$f\colon\mathcal{M}\mapsto\mathbb{R}$ is a smooth function. In order to save
notations, we reuse the notations in the previous section. Here we list our
basic assumption about Eq. 42 in this section. Suppose
$\mathcal{X}\subset\mathcal{M}$ is a geodesically convex region, where the
sectional curvature of $\mathcal{M}$ is bounded in $[-K,K]$. Assume
$x_{*}\in\mathcal{X}$ and $f$ is $\mu$-geodesically convex and $L$-Lipschitz
smooth on $\mathcal{X}$. Moreover, let $\mathcal{B}_{R}$ be the closed ball on
$\mathcal{M}$ defined as

$\mathcal{B}_{R}\equiv\\{x\in\mathcal{M}\mid\lVert\log_{x_{*}}(x)\rVert\leq
R\\},$

assume for any $R>0$, if $\mathcal{B}_{R}\subset\mathcal{X}$, then
$\mathcal{B}_{R}$ is geodesically convex.

######  Remark 3.1.

Different from most literatures about Riemannian acceleration, only local
geodesic convexity is assumed in Eq. 42, which is coincided with our previous
results about preconditioned eigensystem in Theorem 2.36.

###  3.1 Estimate sequence on Riemannian manifold

The estimate sequence, which was first proposed by Nesterov in [34], is a
traditional and popular technique for the analysis of the momentum term. Zhang
and Sra generalized it into the Riemannian manifold, called weak estimate
sequence, as follows.

######  Definition 3.2 (Weak estimate sequence, [54, Definition 1]).

Let $f\colon\mathcal{M}\mapsto\mathbb{R}$, if there exists $0<\alpha<1$ and a
sequence $\phi_{m}\colon\mathcal{M}\mapsto\mathbb{R}$ such that

$\phi_{m}(x_{*})-f(x_{*})\leq(1-\alpha)^{m}\bigl{(}\phi_{0}(x_{*})-f(x_{*})\bigr{)},$

where $x_{*}$ is the solution of Eq. 42, then $(\alpha,\\{\phi_{m}\\})$ is
called a weak estimate sequence.

With the weak estimate sequence, Zhang and Sra proposed a Riemannian
Accelerated Gradient Descent (RAGD) method, i.e., [54, Algorithm 2], and
analyzed its convergence. Here we slightly generalize it as follows. For
initial point $v_{0}=x_{0}\in\mathcal{X}$, step size $h\leq 1/L$ and shrinkage
parameter $\beta>0$, let

$\alpha=\frac{\sqrt{\beta^{2}+4(1+\beta)\mu
h}-\beta}{2},\quad\gamma=\frac{\alpha\mu}{\alpha+\beta}\quad\text{and}\quad\overline{\gamma}=(1+\beta)\gamma.$

For $m=0,1,\dotsc$, consider the following recurrence

(43) $\left\\{\begin{aligned}
&y_{m}=\exp_{x_{m}}\Bigl{(}\dfrac{\alpha}{1+\alpha+\beta}\log_{x_{m}}(v_{m})\Bigr{)},\\\
&v_{m+1}=\exp_{y_{m}}\Bigl{(}\frac{1-\alpha}{1+\beta}\log_{y_{m}}(v_{m})-\frac{\alpha}{\overline{\gamma}}\mathrm{grad}f(y_{m})\Bigr{)},\\\
&\text{update $x_{m+1}$ satisfying $f(x_{m+1})\leq
f(y_{m})-\frac{1}{2L}\lVert\mathrm{grad}f(y_{m})\rVert^{2}$.}\end{aligned}\right.$

In RAGD of [54], the iterates $x_{m+1}$ are constructed as

$x_{m+1}=\exp_{y_{m}}(-L^{-1}\mathrm{grad}f(y_{m})),$

which satisfies the condition for $x_{m+1}$ in Eq. 43 when
$x_{m+1}\in\mathcal{X}$. Since in the estimate sequence approach, the only
requirement for $x_{m+1}$ is $f(x_{m+1})$ should be small enough, thus
theoretical results in [54] can be generalized into Eq. 43. Here we list the
core theoretical result about the convergence of it, whose proof is very
similar to the analysis of RAGD, i.e., [54, Lemmas 2–6].

######  Proposition 3.3.

For the iterates from Eq. 43 with any shrinkage parameter $\beta>0$ and
initial value $\phi_{0}^{*}=f(x_{0})$ and
$\phi_{0}(x)=\phi_{0}^{*}+\frac{\gamma}{2}\lVert\log_{x_{0}}(x)\rVert^{2}$,
let

$\phi_{i+1}(x)=\phi_{i+1}^{*}+\frac{\overline{\gamma}}{2}\lVert\log_{y_{i}}(x)-\log_{y_{i}}(v_{i+1})\rVert^{2},$

for $i\geq 0$, where

$\displaystyle\phi_{i+1}^{*}=$ $\displaystyle(1-\alpha)\phi_{i}^{*}+\alpha
f(y_{i})-\frac{\alpha^{2}}{2\overline{\gamma}}\lVert\mathrm{grad}f(y_{i})\rVert^{2}$
$\displaystyle+\frac{\alpha(1-\alpha)}{1+\beta}\Bigl{(}\frac{\mu}{2}\lVert\log_{y_{i}}(v_{i})\rVert^{2}+\langle\mathrm{grad}f(y_{i}),\log_{y_{i}}(v_{i})\rangle\Bigr{)}.$

If $\\{x_{i}\\}_{i=0}^{m+1}$, $\\{y_{i}\\}_{i=0}^{m}$ and
$\\{v_{i}\\}_{i=0}^{m+1}$ are all in $\mathcal{X}$ and

(44)
$\lVert\log_{y_{i+1}}(x_{*})-\log_{y_{i+1}}(v_{i+1})\rVert^{2}\leq(1+\beta)\lVert\log_{y_{i}}(x_{*})-\log_{y_{i}}(v_{i+1})\rVert^{2}$

holds for $0\leq i\leq m-1$, then the following relationships hold:

  1. 1.

$(\alpha,\\{\phi_{i}\\}_{i=0}^{m+1})$ is a weak estimate sequence.

  2. 2.

$f(x_{i})\leq\phi_{i}^{*}$ holds for all $0\leq i\leq m+1$.

  3. 3.

$\phi_{i+1}(x_{*})-f(x_{*})\leq(1-\alpha)\bigl{(}\phi_{i}(x_{*})-f(x_{*})\bigr{)}$
holds for all $0\leq i\leq m$.

######  Corollary 3.4 ([54, Theorem 7]).

The iterates $\\{x_{i}\\}_{i=0}^{m+1}$ from Eq. 43 satisfy

$\displaystyle f(x_{m+1})-f(x_{*})$
$\displaystyle\leq\phi_{m+1}(x_{*})-f(x_{*})\leq(1-\alpha)^{m+1}\bigl{(}\phi_{0}(x_{*})-f(x_{*})\bigr{)}$
$\displaystyle=(1-\alpha)^{m+1}\Bigl{(}f(x_{0})-f(x_{*})+\frac{\gamma}{2}\lVert\log_{x_{0}}(x_{*})\rVert^{2}\Bigr{)}.$

###  3.2 Algorithm and convergence analysis

In this part, we propose a Locally Optimal Riemannian Accelerated Gradient
(LORAG) method, whose pseudocode is given in Algorithm 1, by combining the
RAGD method with the locally optimal technique from LOBPCG [24].

Input: Initial point $x_{0}$, parameters $\beta$, $\mu$ and $L$.

Set $0<\alpha\leq\frac{\sqrt{\beta^{2}+4(1+\beta)\mu/L}-\beta}{2}$,
$\gamma=\frac{\alpha\mu}{\alpha+\beta}$, $\overline{\gamma}=(1+\beta)\gamma$;

Set $v_{0}=x_{0}$ and $m=0$;

repeat

Choose
$y_{m}=\exp_{x_{m}}\Bigl{(}\dfrac{\alpha}{1+\alpha+\beta}\log_{x_{m}}(v_{m})\Bigr{)}$;

Find $x_{m+1}$ such that
$x_{m+1}=\mathop{\arg\min}\limits_{g\in\mathrm{span}\\{\mathrm{grad}f(y_{m}),\log_{y_{m}}(x_{m})\\}}f\bigl{(}\exp_{y_{m}}(g)\bigr{)}$;

Set
$v_{m+1}=\exp_{y_{m}}\Bigl{(}\frac{1-\alpha}{1+\beta}\log_{y_{m}}(v_{m})-\frac{\alpha}{\overline{\gamma}}\mathrm{grad}f(y_{m})\Bigr{)}$;

$m=m+1$;

until  _Convergence_ ;

Algorithm 1 Locally Optimal Riemannian Accelerated Gradient method

Compared with the RAGD method or Eq. 43, one advantage of the LORAG method is
the monotonicity of function values $f(x_{m})$, which is claimed by the
following lemma.

######  Lemma 3.5.

For any $\\{x_{m}\\}$ from Algorithm 1, $f(x_{m+1})\leq f(x_{m})\leq
f(x_{0})$.

######  Proof 3.6.

By the locally optimal step in Algorithm 1, we have

$f(x_{m+1})\leq
f\Bigl{(}\exp_{y_{m}}\bigl{(}\log_{y_{m}}(x_{m})\bigr{)}\Bigr{)}=f(x_{m})\leq\dotsb\leq
f(x_{0}),$

which finishes the proof.

Besides the monotonicity, we are also interested in the asymptotic convergence
of the LORAG method. Different from most existing Riemannian Nestrov methods,
which are based on the assumption as all iterates lying in $\mathcal{X}$,
where $f$ is geodesically convex and Lipschitz smooth, we discuss how to
obtain the convergence without the prior assumption about iterates lying in
$\mathcal{X}$, which is very general in practical problems, such as eigenvalue
problems in the last section. Actually, similar theory and algorithm for
Euclidean space has been discussed in [46]. From the Riemannian manifold
viewpoint, the framework is close, but technical details are harder.

######  Theorem 3.7.

Under Eq. 42, if there exist $0<R_{1}<R_{2}$ such that

(45) $\\{x\in\mathcal{M}\mid f(x)\leq
f(x_{0})\\}\subset\mathcal{B}_{R_{1}}\subset\mathcal{B}_{R_{2}}\subset\mathcal{X},$

and

(46) $(1+KR_{1}^{2})^{2}\leq 1+\beta,$

where

$\displaystyle R_{1}$
$\displaystyle=\Bigl{(}\frac{2E_{0}}{\gamma}\Bigr{)}^{1/2},\quad
R_{2}=\max\Bigl{\\{}2R_{1},\Bigl{(}1+\frac{\alpha
L}{(1+\beta)\gamma}\Bigr{)}R_{1}\Bigr{\\}}\quad\text{and}$ $\displaystyle
E_{0}$
$\displaystyle=\phi_{0}(x_{*})-f(x_{*})=f(x_{0})-f(x_{*})+\frac{\gamma}{2}\lVert\log_{x_{0}}(x_{*})\rVert^{2},$

then the iterates $\\{x_{m}\\}$ from Algorithm 1 lie in
$\mathcal{B}_{R_{1}}\subset\mathcal{X}$.

######  Proof 3.8.

By Lemma 3.5, we know $f(x_{m})\leq f(x_{0})$ for all $m\geq 0$, then
$x_{m}\in\mathcal{B}_{R_{1}}$ comes directly from Eq. 45. Once
$v_{m}\in\mathcal{B}_{R_{1}}$ is proved, since $\mathcal{B}_{R_{1}}$ is
geodesically convex, we can obtain $y_{m}\in\mathcal{B}_{R_{1}}$ from the fact
$y_{m}$ is in the geodesic from $x_{m}$ to $v_{m}$ directly.

Now, we only need to prove $v_{m}\in\mathcal{B}_{R_{1}}$ for all $m\geq 0$.
The proof is based on the mathematical induction. The case $m=0$ is straight
from $v_{0}=x_{0}\in\mathcal{B}_{R_{1}}$. Assume $v_{i}\in\mathcal{B}_{R_{1}}$
for all $0\leq i\leq m$.

First, we need to prove

(47) $f(x_{i+1})\leq
f(y_{i})-\frac{1}{2L}\lVert\mathrm{grad}f(y_{i})\rVert^{2}$

holds for all $0\leq i\leq m$. Since
$y_{i}\in\mathcal{B}_{R_{1}}\subset\mathcal{X}$, we know that

$\lVert\mathrm{grad}f(y_{i})\rVert=\lVert\mathrm{grad}f(y_{i})-\Gamma_{x_{*}}^{y_{i}}\mathrm{grad}f(x_{*})\rVert\leq
L\lVert\log_{x_{*}}(y_{i})\rVert\leq LR_{1},$

where $\Gamma_{x_{*}}^{y_{i}}$ is the parallel transport. By the triangular
inequality, let

$\overline{y}_{i}=\exp_{y_{i}}\Bigl{(}-\frac{1}{L}\mathrm{grad}f(y_{i})\Bigr{)},$

we know

$\lVert\log_{x_{*}}(\overline{y}_{i})\rVert\leq\lVert\log_{x_{*}}(y_{i})\rVert+\frac{1}{L}\lVert\mathrm{grad}f(y_{i})\rVert\leq
2R_{1}\leq R_{2}.$

Using the locally optimal step in Algorithm 1 and the $L$-Lipschitz
smoothness, we obtain Eq. 47.

Then we will show

(48) $v_{m+1}\in\mathcal{B}_{R_{2}}\subset\mathcal{X}.$

From the recurrence of $v_{m+1}$ in Algorithm 1 and the triangular inequality,
we know

$\lVert\log_{x_{*}}(v_{m+1})\rVert\leq\lVert\log_{x_{*}}(\overline{v}_{m})\rVert+\frac{\alpha}{\overline{\gamma}}\lVert\mathrm{grad}f(y_{m})\rVert,$

where

$\overline{v}_{m}=\exp_{y_{m}}\Bigl{(}\frac{1-\alpha}{1+\beta}\log_{y_{m}}(v_{m})\Bigr{)}$

is on the geodesic from $y_{m}$ to $v_{m}$. Similar the analysis of
$\overline{y}_{i}$, we know that

$\lVert\log_{x_{*}}(v_{m+1})\rVert\leq R_{1}+\frac{\alpha
L}{\overline{\gamma}}R_{1}=\Bigl{(}1+\frac{\alpha
L}{(1+\beta)\gamma}\Bigr{)}R_{1}\leq R_{2},$

which implies Eq. 48.

Now we can use Propositions 3.3 and 2.21 to obtain

$\phi_{m+1}(x_{*})-f(x_{*})\leq(1-\alpha)^{m+1}E_{0}\leq E_{0}.$

Note that $\phi_{m+1}^{*}\geq f(x_{m+1})\geq f(x_{*})$ and

$\phi_{m+1}(x_{*})-f(x_{*})=\phi_{m+1}^{*}-f(x_{*})+\frac{\overline{\gamma}}{2}\lVert\log_{y_{m}}(x)-\log_{y_{m}}(v_{m+1})\rVert^{2},$

we have

$\frac{\overline{\gamma}}{2}\lVert\log_{y_{m}}(x)-\log_{y_{m}}(v_{m+1})\rVert^{2}\leq\phi_{m+1}(x_{*})-f(x_{*})\leq
E_{0}.$

Using Propositions 2.21 and 46 we know

$\lVert\log_{x_{*}}(v_{m+1})\rVert^{2}\leq(1+KR_{1}^{2})\lVert\log_{y_{m}}(x)-\log_{y_{m}}(v_{m+1})\rVert^{2}\leq\frac{2(1+KR_{1}^{2})E_{0}}{(1+\beta)\gamma}\leq
R_{1}^{2},$

which finishes the proof.

Combining Theorems 3.7 and 3.3, we can obtain the convergence of Algorithm 1
as follows.

######  Corollary 3.9.

With notations and conditions Eqs. 45 and 46 in Theorem 3.7, the iterates from
Algorithm 1 satisfy

(49) $\displaystyle f(x_{m})\leq f(x_{m-1})\leq\dotsb\leq f(x_{0}),$ (50)
$\displaystyle f(x_{m})-f(x_{*})\leq(1-\alpha)^{m}E_{0}.$

######  Proof 3.10.

The monotonicity comes from Lemma 3.5. For the asymptotic convergence, from
Theorem 3.7 we know all iterates $(y_{m},x_{m+1},v_{m+1})$ are all in
$\mathcal{B}_{R_{1}}\subset\mathcal{X}$, then the convergence is obtained from
Proposition 3.3.

######  Proposition 3.11.

Denote $\kappa=L/\mu$, assume $\kappa\geq 9$, let

$\displaystyle\beta=\frac{3}{2\sqrt{\kappa}-4},\quad\alpha=\frac{1}{2\sqrt{\kappa}},\quad\gamma=\frac{(\sqrt{\kappa}-2)\mu}{2(2\sqrt{\kappa}-1)}$
$\displaystyle
R_{1}=\Bigl{(}\frac{2E_{0}}{\gamma}\Bigr{)}^{1/2}\quad\text{and}\quad
R_{2}=(1+2\sqrt{\kappa})R_{1},$

where
$E_{0}=f(x_{0})-f(x_{*})+\frac{\gamma}{2}\lVert\log_{x_{*}}(x_{0})\rVert^{2}$.
When $x_{0}$ is sufficiently close to $x_{*}$, i.e.,

$(1+KR_{1}^{2})^{2}\leq
1+\beta\quad\text{and}\quad\mathcal{B}(R_{2})\subset\mathcal{X},$

then all conditions in Eqs. 46 and 45 holds, and the asymptotic rate of
convergence is

(51)
$f(x_{m})-f(x_{*})\leq\Bigl{(}1-\frac{1}{2\sqrt{\kappa}}\Bigr{)}^{m}E_{0}.$

######  Remark 3.12.

The coefficients in Proposition 3.11 are not essential. Actually, when $\beta$
is sufficiently small, $\alpha$ can be any positive scalar smaller than
$\kappa^{-1/2}$.

######  Proposition 3.13.

Under conditions in Proposition 3.11, we have

$R_{1}^{2}\leq\frac{22}{\mu}\bigl{(}f(x_{0})-f(x_{*})\bigr{)}\quad\text{and}\quad
R_{2}^{2}\leq\frac{22}{\mu}(1+2\sqrt{\kappa})^{2}\bigl{(}f(x_{0})-f(x_{*})\bigr{)}.$

######  Proof 3.14.

For $R_{1}^{2}$, by the geodesic convexity of $f$, we have

$E_{0}=f(x_{0})-f(x_{*})+\frac{\gamma}{2}\lVert
x_{0}-x_{*}\rVert^{2}\leq\frac{\mu+\gamma}{\mu}\bigl{(}f(x_{0})-f(x_{*})\bigr{)}.$

Note that $0<\alpha\leq 1/6$, we know $4<\mu/\gamma\leq 10$ and

$R_{1}^{2}\leq\frac{2E_{0}}{\gamma}=\frac{2(\mu+\gamma)}{\gamma\mu}\bigl{(}f(x_{0})-f(x_{*})\bigr{)}\leq\frac{22}{\mu}\bigl{(}f(x_{0})-f(x_{*})\bigr{)}.$

The estimation for $R_{2}$ comes from $R_{2}=(1+2\sqrt{\kappa})R_{1}$.

###  3.3 Solving symmetric eigenvalue problem by preconditioned LORAG method

Now let us turn back to the symmetric eigenvalue problem. As discussed in
Theorem 2.36, under some mild conditions, the symmetric eigenvalue problem is
geodesically convex and Lipschitz smooth on the manifold, which means we can
solve it by Algorithm 1.

One practical difficulty in the computation of Algorithm 1 is the
multiplication of the preconditioner $B$. For example, when computing $y_{m}$,
we need to compute

$\displaystyle\theta_{m}$
$\displaystyle=\dfrac{\alpha}{\alpha+\beta+1}\arccos(x_{m}^{\mathsf{T}}Bv_{m}),$
$\displaystyle w_{m}$ $\displaystyle=v_{m}-(v_{m}^{\mathsf{T}}Bx_{m})x_{m},$
$\displaystyle y_{m}$
$\displaystyle=\cos(\theta_{m})x_{m}+\sin(\theta_{m})\dfrac{w_{m}}{\lVert
w_{m}\rVert_{B}}.$

In most situation, the preconditioner $B$ is only available for solving linear
system rather than multiplication. In order to overcome it, we use some “co-
iterates”:

(52)
$\widehat{x}_{m}=Bx_{m},\quad\widehat{v}_{m}=Bv_{m},\quad\widehat{y}_{m}=By_{m}.$

With Eq. 52, we can compute the $B$ inner-product and norm by a standard inner
product of one iterate and one co-iterate. Moreover, since most iterations in
Algorithm 1 is linear, the co-iterates can be obtained without $B$ or
$B^{-1}$. With the co-iterates in Eq. 52, we propose the Riemannian
Acceleration with Preconditioning (RAP) for symmetric eigenvalues problems as
Algorithm 2.

Input: Initial $\widehat{x}_{0}$, matrices $(A,M)$, preconditioner $B$,
parameters $(\mu,L)$.

Set $\kappa=\dfrac{L}{\mu}$, $\beta=\dfrac{3}{2\sqrt{\kappa}-4}$,
$\alpha=\dfrac{\sqrt{\beta^{2}+4(1+\beta)\kappa^{-1}}-\beta}{2}$,
$\gamma=\dfrac{\alpha\mu}{\alpha+\beta}$;

Compute $x_{0}=B^{-1}\widehat{x}_{0}$;

Normalize $x_{0}$ and $\widehat{x}_{0}$ with
$x_{0}=\dfrac{x_{0}}{\sqrt{x_{0}^{\mathsf{T}}\widehat{x}_{0}}}$ and
$\widehat{x}_{0}=\dfrac{\widehat{x}_{0}}{\sqrt{x_{0}^{\mathsf{T}}\widehat{x}_{0}}}$;

Set $\widehat{v}_{0}=\widehat{x}_{0}$, $v_{0}=x_{0}$ and $m=0$;

repeat

Update $y_{m}$ and $\widehat{y}_{m}$:
$\displaystyle\theta_{m}=\dfrac{\alpha}{\alpha+\beta+1}\arccos(x_{m}^{\mathsf{T}}\widehat{v}_{m}),$
$\displaystyle
w_{m}=v_{m}-(x_{m}^{\mathsf{T}}\widehat{v}_{m})x_{m}\quad\text{and}\quad\widehat{w}_{m}=\widehat{v}_{m}-(x_{m}^{\mathsf{T}}\widehat{v}_{m})\widehat{x}_{m},$
$\displaystyle\eta_{m}^{w}=\sqrt{w_{m}^{\mathsf{T}}\widehat{w}_{m}},\quad
w_{m}=\dfrac{w_{m}}{\eta_{m}^{w}}\quad\text{and}\quad\widehat{w}_{m}=\dfrac{\widehat{w}_{m}}{\eta_{m}^{w}},$
$\displaystyle
y_{m}=x_{m}\cos\theta_{m}+w_{m}\sin\theta_{m}\quad\text{and}\quad\widehat{y}_{m}=\widehat{x}_{m}\cos\theta_{m}+\widehat{w}_{m}\sin\theta_{m}.$

Update $v_{m+1}$ and $\widehat{v}_{m+1}$:  $\displaystyle\sigma_{m}$
$\displaystyle=\dfrac{y_{m}^{\mathsf{T}}Ay_{m}}{y_{m}^{\mathsf{T}}My_{m}},$
$\displaystyle p_{m}$
$\displaystyle=v_{m}-(y_{m}^{\mathsf{T}}\widehat{v}_{m})y_{m}\quad\text{and}\quad\widehat{p}_{m}=\widehat{v}_{m}-(y_{m}^{\mathsf{T}}\widehat{v}_{m})\widehat{y}_{m},$
$\displaystyle\eta_{m}^{p}$
$\displaystyle=\sqrt{p_{m}^{\mathsf{T}}\widehat{p}_{m}},\quad
p_{m}=\dfrac{p_{m}}{\eta_{m}^{p}}\quad\text{and}\quad\widehat{p}_{m}=\dfrac{\widehat{p}_{m}}{\eta_{m}^{p}},$
$\displaystyle\widehat{g}_{m}$
$\displaystyle=\dfrac{2}{y_{m}^{\mathsf{T}}My_{m}}(Ay_{m}-\sigma_{m}My_{m})\quad\text{and}\quad
g_{m}=B^{-1}\widehat{g}_{m},$ $\displaystyle q_{m}$
$\displaystyle=\dfrac{(1-\alpha)\theta_{m}}{\alpha}p_{m}-\dfrac{\alpha}{(1+\beta)\gamma}g_{m},$
$\displaystyle\widehat{q}_{m}$
$\displaystyle=\dfrac{(1-\alpha)\theta_{m}}{\alpha}\widehat{p}_{m}-\dfrac{\alpha}{(1+\beta)\gamma}\widehat{g}_{m},$
$\displaystyle\eta_{m}^{q}$
$\displaystyle=\sqrt{q_{m}^{\mathsf{T}}\widehat{q}_{m}},\quad
q_{m}=\frac{q_{m}}{\eta_{m}^{q}}\quad\text{and}\quad\widehat{q}_{m}=\frac{\widehat{q}_{m}}{\eta_{m}^{q}},$
$\displaystyle v_{m+1}$
$\displaystyle=y_{m}\cos\eta_{m}^{q}+q_{m}\sin\eta_{m}^{q}$
$\displaystyle\widehat{v}_{m+1}$
$\displaystyle=\widehat{y}_{m}\cos\eta_{m}^{q}+\widehat{q}_{m}\sin\eta_{m}^{q}.$

Update $x_{m+1}$ and $\widehat{x}_{m+1}$:  $\displaystyle
x_{m+1}=\mathop{\arg\min}\limits_{x\in\mathrm{span}\\{x_{m},y_{m},g_{m}\\}}\dfrac{x^{\mathsf{T}}Ax}{x^{\mathsf{T}}Mx}=c_{m}^{(1)}x_{m}+c_{m}^{(2)}y_{m}+c_{m}^{(3)}g_{m},$
$\displaystyle\widehat{x}_{m+1}=c_{m}^{(1)}\widehat{x}_{m}+c_{m}^{(2)}\widehat{y}_{m}+c_{m}^{(3)}\widehat{g}_{m},$
$\displaystyle\eta_{m+1}^{x}=\sqrt{x_{m+1}^{\mathsf{T}}\widehat{x}_{m+1}},\quad
x_{m+1}=\frac{x_{m+1}}{\eta_{m+1}^{x}}\quad\text{and}\quad\widehat{x}_{m+1}=\frac{\widehat{x}_{m+1}}{\eta_{m+1}^{x}}.$
$m=m+1$;

until  _Convergence_ ;

Algorithm 2 Riemannian Acceleration with Preconditioning (RAP)

Note that the convexity radius of the sphere is $\pi/2$ [22, Page 84], i.e.,
any geodesic ball in $\mathcal{X}_{B}$ is convex, combining it with Theorems
3.7, 3.11 and 2.36, we can obtain the convergence of Algorithm 2.

######  Theorem 3.15.

Let $f(x)$ be the Rayleigh quotient defined in Eq. 2, suppose
$\lambda_{1}\leq\rho_{\mathcal{X}}<(\lambda_{1}+\lambda_{2})/2$ and the
parameters $\mu$ and $L$ defined in Eq. 38 satisfy $L\geq 9\mu>0$. When the
initial point $x_{0}$ satisfies

(53)
$f(x_{0})-\lambda_{1}\leq\frac{(\rho_{\mathcal{X}}-\lambda_{1})}{22(1+2\sqrt{\kappa})^{2}\kappa},$

where $\kappa=L/\mu$, then Algorithm 2 is monotonic convergence. The
asymptotic rate of convergence is

(54) $f(x_{m})-\lambda_{1}\leq
2\Bigl{(}1-\frac{1}{2\sqrt{\kappa}}\Bigr{)}^{m}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)}.$

######  Proof 3.16.

It is sufficient to verify conditions in Proposition 3.11. Using the
$L$-Lipschitz smooth of $f$, we know when

$\frac{LR_{2}^{2}}{2}\leq\rho_{\mathcal{X}}-\lambda_{1},$

it can be proved that $\mathcal{B}_{R_{2}}\subset\mathcal{X}$. Then, according
to Proposition 3.13, it is sufficient to show

(55) $\displaystyle f(x_{0})-f(x_{*})$
$\displaystyle\leq\frac{\mu}{22}\bigl{(}\sqrt{1+\beta}-1\bigr{)}$ (56)
$\displaystyle f(x_{0})-f(x_{*})$
$\displaystyle\leq\frac{(\rho_{\mathcal{X}}-\lambda_{1})}{22(1+2\sqrt{\kappa})^{2}\kappa}.$

Note that Eq. 56 is same as the condition Eq. 53, we need to prove Eq. 53
implies Eq. 55. Since $0<\beta\leq 3/2$, we obtain that

$\sqrt{1+\beta}-1\geq\frac{\beta}{3}=\frac{1}{2\sqrt{\kappa}-4}\geq\frac{1}{2\sqrt{\kappa}}.$

Combining it with $\kappa=L/\mu$, it is sufficient to show

$\frac{\mu(\rho_{\mathcal{X}}-\lambda_{1})}{(1+2\sqrt{\kappa})^{2}L}\leq\frac{\mu}{2\sqrt{\kappa}},$

which is equivalent to

$L\geq\frac{2\sqrt{\kappa}}{(1+2\sqrt{\kappa})^{2}}(\rho_{\mathcal{X}}-\lambda_{1}).$

According to Eq. 39, we have

$L=\frac{2\nu_{\max}\rho_{\mathcal{X}}}{\sigma-\varsigma}\Bigl{(}1-\frac{\lambda_{1}}{\lambda_{n}}\Bigr{)}+C_{\mathcal{X}}\geq\frac{2\nu_{\max}\rho_{\mathcal{X}}}{\sigma}\Bigl{(}1-\frac{\lambda_{1}}{\rho_{\mathcal{X}}}\Bigr{)}\geq
2(\rho_{\mathcal{X}}-\lambda_{1}).$

Then the theorem is proved by $\kappa\geq 9$.

##  4 Second order symmetric elliptic eigenvalue problems with Schwarz
preconditioners

In this part, we turn to the eigensolver based on Schwarz preconditioners [50,
Section 2], and give an illustration on why $\vartheta$ is a better
measurement for preconditioned eigenvalue problems. Recently, Schwarz
preconditioners, especially the Domain Decomposition Method (DDM), becomes a
popular technique for elliptic eigenvalue problems [52, 45, 13], and in our
previous work [13], the acceleration by involving the momentum is observed in
the numerical examples. In this part, we give a theoretical explanation about
this phenomenon.

Since our discussion in this part stands on PDE and FEM viewpoints, where we
are more interested in the relationship with mesh sizes $h$ and $H$, the
notation “$\lesssim$” also means independent of $h$ and $H$.

###  4.1 Theory of Schwarz methods and domain decomposition method

Suppose $\Omega\subset\mathbb{R}^{d}$ for $d=2$ or $3$ is a convex polygonal
domain, the matrix $\\{a_{ij}(x)\\}_{i,\,j=1}^{d}$ is symmetric and uniformly
positive definite and $a_{ij}(x)\in C^{0,1}(\overline{\Omega})$ for
$i,j=1,\dotsc,d$. Let $V_{0}\subset V\subset H_{0}^{1}(\Omega)$ be the
continuous and piecewise linear element spaces based on quasi-uniform
triangular partitions $\mathcal{T}_{H}$ and $\mathcal{T}_{h}$, where
$\mathcal{T}_{h}$ is refined by $\mathcal{T}_{H}$, and $0<h<H<1$ are mesh
sizes of $\mathcal{T}_{H}$ and $\mathcal{T}_{h}$, respectively. We consider
the following second order symmetric elliptic eigenvalue problem

(57) $\displaystyle a(u_{1},v)\equiv$
$\displaystyle\sum_{i,j=1}^{d}\int_{\Omega}a_{ij}(x)\frac{\partial
u_{1}}{\partial x_{i}}\frac{\partial v}{\partial
x_{j}}\,\mathrm{d}x=\lambda_{1}\int_{\Omega}u_{1}v\,\mathrm{d}x\quad\forall\,v\in
V,$ $\displaystyle\int_{\Omega}u_{1}^{2}\,\mathrm{d}x=1\quad\text{and}\quad
u_{1}\in V.$

Consider a family of subspaces of $V$ as $\\{V_{i}\\}_{i=1}^{N}$, assume $V$
admits the following decomposition

(58) $V=R_{0}^{\mathsf{T}}V_{0}+\sum_{i=1}^{N}R_{i}^{\mathsf{T}}V_{i},$

where $R_{i}^{\mathsf{T}}\colon V_{i}\mapsto V$ are the interpolation
operators for $i=0,\dotsc,d$. The additive Schwarz preconditioner can be
defined as

(59) $B^{-1}=\sum_{i=0}^{N}R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i},$

where $B_{i}^{-1}$ is the solver on $V_{i}$. Here are some assumptions about
the additive Schwarz preconditioner. [[50, Assumptions 2.2–2.4]]

  * •

Stable decomposition: There exists a constant $C_{s}$ such that every $v\in V$
admits a decomposition $v=\sum_{i=0}^{N}R_{i}^{\mathsf{T}}v_{i}$ that
satisfies $\sum_{i=1}^{N}\lVert v_{i}\rVert_{B_{i}}^{2}\leq C_{s}^{2}\lVert
v\rVert_{A}^{2}$.

  * •

Strengthened Cauchy-Schwarz inequalities: There exist constants
$0\leq\epsilon_{ij}\leq 1$, $1\leq i,j\leq N$, such that
$(R_{i}^{\mathsf{T}}v_{i},R_{j}^{\mathsf{T}}v_{j})_{A}\leq\epsilon_{ij}\lVert
R_{i}^{\mathsf{T}}v_{i}\rVert_{A}\lVert R_{j}^{\mathsf{T}}v_{j}\rVert_{A}$ for
$v_{i}\in V_{i}$ and $v_{j}\in V_{j}$.

  * •

Finite covering: Assume there are at most $N_{0}$ nonzero elements at each
column of $\mathcal{E}=\\{\epsilon_{ij}\\}_{i,j=1}^{N}$, where $N_{0}$ is
independent of $H$ and $h$. Then the spectral radius of $\mathcal{E}$ is
bounded by $N_{0}$.

  * •

Local stability: There exists $\omega>0$, such that $\lVert
R_{i}^{\mathsf{T}}v_{i}\rVert_{A}^{2}\leq\omega\lVert v_{i}\rVert_{B_{i}}^{2}$
holds for all $v_{i}\in V_{i}$ and $0\leq i\leq N$.

######  Proposition 4.1 ([50, Theorem 2.7]).

Under Eq. 59, the smallest and largest eigenvalues of $A,B$, i.e.,
$\nu_{\min}$ and $\nu_{\max}$, and the condition number $\kappa_{\nu}$ satisfy

(60) $\nu_{\min}\geq 1,\quad\nu_{\max}\leq
C_{s}^{2}\omega(N_{0}+1),\quad\kappa_{\nu}\leq C_{s}^{2}\omega(N_{0}+1).$

Specifically, in DDM, $N_{0}\approx 1$, $\omega\approx 1$, and $C_{s}\lesssim
1$ for the two-level overlapping method [16] and $C_{s}\lesssim 1+\ln(H/h)$
for the Bramble-Pasciak-Schatz (BPS) method [12].

Besides traditional assumptions for additive Schwarz theory, we introduce some
parameters, and give their estimations in DDM, whose details can be found in
[50, 7].

  * •

Let $(\lambda_{0},u_{0})$ be the smallest eigenpair of $B_{0}$, i.e.,

$(B_{0}u_{0},v_{0})=\lambda_{0}(u_{0},v_{0})\quad\forall\,v_{0}\in V_{0},$

suppose $\lVert u_{0}\rVert=1$, then $C_{c}$ is defined as

(61) $C_{c}\equiv\min_{\alpha\in\mathbb{R}}\lVert u_{1}-\alpha
R_{0}^{\mathsf{T}}u_{0}\rVert_{A}.$

In DDM, $C_{c}\lesssim H$ and $0\leq\lambda_{0}-\lambda_{1}\lesssim H^{2}$.

  * •

The parameter $C_{P}$ is defined as

(62) $C_{P}\equiv\max_{1\leq i\leq N}\lambda_{\max}(B_{i}^{-1})=\max_{1\leq
i\leq N}\lambda_{\min}^{-1}(B_{i}).$

In DDM, $C_{P}\lesssim H^{2}$.

  * •

The parameter $C_{a}$ is defined as

(63) $C_{a}=\max_{v\in V}\frac{\sum_{i=1}^{N}\lVert R_{i}v\rVert^{2}}{\lVert
v\rVert^{2}}.$

In DDM, $C_{a}\lesssim 1$.

At the end of this subsection, we summarize the parameters defined above in
Table 1.

Table 1: Bounds for parameters about DDM. $C_{s}$ | $N_{0}$ | $\omega$ | $C_{c}$ | $\lambda_{0}$ | $C_{P}$ | $C_{a}$
---|---|---|---|---|---|---
$\mathcal{O}(1)$ or $\mathcal{O}(1+\ln(H/h))$ | $\mathcal{O}(1)$ | $\mathcal{O}(1)$ | $\mathcal{O}(H)$ | $\lambda_{1}+\mathcal{O}(H^{2})$ | $\mathcal{O}(H^{2})$ | $\mathcal{O}(1)$

###  4.2 Analysis for leading angle

Now let use consider the leading angle defined in Definition 2.5 in the
Schwarz framework. In order to use Eq. 17, we first consider the global term
$R_{0}^{\mathsf{T}}B_{0}^{-1}R_{0}u_{1}-\alpha u_{1}$ and the local term
$\sum_{i=1}^{N}R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i}u_{1}$, respectively.

######  Lemma 4.2.

Let $u_{1}$ be the eigenvector of Eq. 57, under Eq. 59,

(64) $\lVert
R_{0}^{\mathsf{T}}B_{0}^{-1}R_{0}u_{1}-\lambda_{0}^{-1}u_{1}\rVert_{A}\leq\frac{\omega\lambda_{1}^{-1/2}+1}{\lambda_{0}}\,C_{c}\lesssim
H.$

######  Proof 4.3.

For any $\alpha\in\mathbb{R}$, by Eq. 59

$\displaystyle\lVert
R_{0}^{\mathsf{T}}B_{0}^{-1}R_{0}u_{1}-\lambda_{0}^{-1}u_{1}\rVert_{A}$
$\displaystyle\leq\lVert R_{0}^{\mathsf{T}}B_{0}^{-1}(R_{0}u_{1}-\alpha
u_{0})\rVert_{A}+\lambda_{0}^{-1}\lVert\alpha
R_{0}^{\mathsf{T}}u_{0}-u_{1}\rVert_{A}$ $\displaystyle\leq\omega\lVert
R_{0}u_{1}-\alpha u_{0}\rVert_{B_{0}^{-1}}+\lambda_{0}^{-1}\lVert\alpha
R_{0}^{\mathsf{T}}u_{0}-u_{1}\rVert_{A}$
$\displaystyle\leq\lambda_{0}^{-1}\bigl{(}\omega\lVert\alpha
u_{0}-R_{0}u_{1}\rVert+\lVert\alpha
R_{0}^{\mathsf{T}}u_{0}-u_{1}\rVert_{A}\bigr{)}.$

Note that for any $v_{0}$ and $w_{0}$ in $V_{0}$,

$(v_{0},w_{0})=(R_{0}^{\mathsf{T}}v_{0},R_{0}^{\mathsf{T}}w_{0})=(v_{0},R_{0}R_{0}^{\mathsf{T}}w_{0}),$

we know $R_{0}R_{0}^{\mathsf{T}}$ is the identity operator on $V_{0}$. For any
$v$ in $V$, since

$\lVert R_{0}v\rVert=\sup_{v_{0}\in V_{0}}\frac{(R_{0}v,v_{0})}{\lVert
v_{0}\rVert}=\sup_{v_{0}\in V_{0}}\frac{(v,R_{0}^{\mathsf{T}}v_{0})}{\lVert
v_{0}\rVert}\leq\lVert v\rVert\sup_{v_{0}\in V_{0}}\frac{\lVert
R_{0}^{\mathsf{T}}v_{0}\rVert}{\lVert v_{0}\rVert}=\lVert v\rVert,$

we have

$\lVert\alpha u_{0}-R_{0}u_{1}\rVert=\lVert\alpha
R_{0}R_{0}^{\mathsf{T}}u_{0}-R_{0}u_{1}\rVert\leq\lVert\alpha
R_{0}^{\mathsf{T}}u_{0}-u_{1}\rVert\leq\lambda_{1}^{-1/2}\lVert\alpha
R_{0}^{\mathsf{T}}u_{0}-u_{1}\rVert_{A}.$

Then the lemma is proved by Eq. 61.

######  Lemma 4.4.

Under Eq. 59, let

(65) $B_{L}^{-1}=\sum_{i=1}^{N}R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i},$

then for any $u$ with $\lVert u\rVert=1$,

(66) $\lVert B_{L}^{-1}u\rVert_{A}^{2}\leq\omega N_{0}C_{a}C_{P}\lesssim
H^{2}.$

######  Proof 4.5.

By the Cauchy-Schwarz inequality and Eqs. 59 and 62,

(67) $\lVert
R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i}u\rVert_{A}^{2}\leq\omega(R_{i}u,B_{i}^{-1}R_{i}u)\leq\omega
C_{P}\lVert R_{i}u\rVert^{2}$

Combining it with Eqs. 59 and 63, we have

$\lVert B_{L}^{-1}u\rVert_{A}^{2}\leq N_{0}\sum_{i=1}^{N}\lVert
R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i}u\rVert_{A}^{2}\leq\omega
N_{0}C_{P}\sum_{i=1}^{N}\lVert R_{i}u\rVert^{2}\leq\omega N_{0}C_{a}C_{P},$

which finishes the proof.

######  Lemma 4.6.

Let $u_{1}$ be the eigenvector of Eq. 57, under Eq. 59,

$\nu_{1}=\frac{\lVert B^{-1}Au_{1}\rVert_{B}^{2}}{\lVert
u_{1}\rVert_{A}^{2}}\leq\lambda_{0}^{-1}\lambda_{1}(\omega\lambda_{1}^{-1/2}+1)C_{c}+\lambda_{1}(\omega
N_{0}C_{a}C_{P})^{1/2}+\lambda_{0}^{-1}\lambda_{1}\lesssim 1.$

######  Proof 4.7.

Without loss of generality, assume $\lVert u_{1}\rVert_{A}=1$. By Lemmas 4.2
and 4.4,

$\displaystyle\lVert B^{-1}Au_{1}\rVert_{B}^{2}$
$\displaystyle=(u_{1},B^{-1}Au_{1})_{A}\leq\lVert
B^{-1}Au_{1}\rVert_{A}\leq\lambda_{1}\lVert
B^{-1}u_{1}-\lambda_{0}^{-1}u_{1}\rVert_{A}+\lambda_{0}^{-1}\lambda_{1}$
$\displaystyle\leq\lambda_{1}\lVert
R_{0}^{\mathsf{T}}B_{0}^{-1}R_{0}u_{1}-\lambda_{0}^{-1}u_{1}\rVert_{A}+\lambda_{1}\lVert
B_{L}^{-1}u_{1}\rVert_{A}^{2}+\lambda_{0}^{-1}\lambda_{1}$
$\displaystyle\leq\lambda_{0}^{-1}\lambda_{1}(\omega\lambda_{1}^{-1/2}+1)C_{c}+\lambda_{1}(\omega
N_{0}C_{a}C_{P})^{1/2}+\lambda_{0}^{-1}\lambda_{1},$

which finishes the proof.

Combining Lemmas 4.2, 4.4 and 4.6, we obtain the bound for $\vartheta$ defined
Definition 2.5.

######  Theorem 4.8.

With Eq. 59 and parameters for DDM in Table 1, the leading angle $\vartheta$
in Eq. 17 can be bounded by

$\cos\vartheta\lesssim\sqrt{\kappa_{\nu}}\Bigl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{2}-\lambda_{1}}\Bigr{)}^{1/2}+H,$

where the condition number $\kappa_{\nu}$ can be bounded by

(68) $\kappa_{\nu}\lesssim\left\\{\begin{aligned} &1\quad&\text{for two-level
overlapping method,}\\\ &\bigl{(}1+\ln(H/h)\bigr{)}^{2}\quad&\text{for BPS
method.}\end{aligned}\right.$

###  4.3 Analysis for asymptotic convergence rate

In Algorithm 2, the initial point $x_{0}$ is computed as
$x_{0}=B^{-1}\widehat{x}_{0}$, where $\widehat{x}_{0}$ is the input. The
following lemma gives an estimation about $f(x_{0})$.

######  Lemma 4.9.

Assume the $\widehat{x}_{0}$ is the eigenvector calculated on the coarse space
$V_{0}$, i.e., $\widehat{x}_{0}=R_{0}^{\mathsf{T}}u_{0}$, let
$x_{0}=B^{-1}\widehat{x}_{0}$, then,

$f(x_{0})-\lambda_{1}\lesssim H^{2}.$

######  Proof 4.10.

By the estimation of eigenvector, we know

$f(\widehat{x}_{0})-\lambda_{1}=\lambda_{0}-\lambda_{1}\lesssim H^{2}.$

According to $x_{0}=B^{-1}\widehat{x}_{0}$, we have

$x_{0}=\lambda_{0}^{-1}\widehat{x}_{0}+B_{L}^{-1}\widehat{x}_{0},$

where $B_{L}$ is defined in Lemma 4.4. Then, by Lemma 4.4,

(69)
$\rho_{0}=f(x_{0})=\frac{\lVert\lambda_{0}^{-1}\widehat{x}_{0}+B_{L}^{-1}\widehat{x}_{0}\rVert_{A}^{2}}{\lVert\lambda_{0}^{-1}\widehat{x}_{0}+B_{L}^{-1}\widehat{x}_{0}\rVert^{2}}\lesssim\frac{\lambda_{0}+\lambda_{0}(\widehat{x}_{0},B_{L}^{-1}\widehat{x}_{0})_{A}+H^{2}}{1+\lambda_{0}(\widehat{x}_{0},B_{L}^{-1}\widehat{x}_{0})+\lVert
B_{L}^{-1}\widehat{x}_{0}\rVert^{2}}.$

Due to Eq. 59, we know

(70)
$(\widehat{x}_{0},B_{L}^{-1}\widehat{x}_{0})=\sum_{i=1}^{N}(R_{i}\widehat{x}_{0},B_{i}^{-1}R_{i}\widehat{x}_{0})\geq\frac{1}{\omega}\sum_{i=1}^{N}\lVert
R_{i}^{\mathsf{T}}B_{i}^{-1}R_{i}\widehat{x}_{0}\rVert_{A}^{2}\geq 0.$

Thus, the denominator in Eq. 69 is greater than or equal to $1$. For the
numerator, let

$r_{0}=\alpha_{0}\widehat{x}_{0}-u_{1}=\alpha_{0}R_{0}^{\mathsf{T}}u_{0}-u_{1},$

where

$\alpha_{0}=\mathop{\arg\min}_{\alpha\in\mathbb{R}}\lVert u_{1}-\alpha
R_{0}^{\mathsf{T}}u_{0}\rVert_{A}=\frac{(R_{0}^{\mathsf{T}}u_{0},u_{1})_{A}}{\lambda_{0}}\approx
1.$

Using the Cauchy-Schwarz inequality and Eq. 61, we have

$\displaystyle\alpha_{0}^{2}(\widehat{x}_{0},B_{L}^{-1}\widehat{x}_{0})_{A}$
$\displaystyle=\alpha_{0}(r_{0},B_{L}^{-1}r_{0})_{A}+2(u_{1},B_{L}^{-1}r_{0})_{A}+(u_{1},B_{L}^{-1}u_{1})_{A}$
$\displaystyle\leq\alpha_{0}\lVert r_{0}\rVert_{A}\lVert
B_{L}^{-1}r_{0}\rVert_{A}+2\sqrt{\lambda_{1}}\lVert
B_{L}^{-1}r_{0}\rVert_{A}+(u_{1},B_{L}^{-1}u_{1})_{A}$
$\displaystyle\lesssim\lVert
B_{L}^{-1}r_{0}\rVert_{A}+(u_{1},B_{L}^{-1}u_{1})_{A}.$

By the Poincaré inequality, Lemmas 4.4 and 61,

$\lVert B_{L}^{-1}r_{0}\rVert_{A}\lesssim H\lVert r_{0}\rVert\lesssim H\lVert
r_{0}\rVert_{A}\lesssim H^{2}.$

For the other term, similar to Eq. 67,

$(u_{1},B_{L}^{-1}u_{1})_{A}=\lambda_{1}(u_{1},B_{L}^{-1}u_{1})=\lambda_{1}\sum_{i=1}^{N}(R_{i}u_{1},B_{i}^{-1}R_{i}u_{1})\lesssim
H^{2}\sum_{i=1}^{B}\lVert R_{i}u_{1}\rVert^{2}\lesssim H^{2}.$

Combining these three inequalities above, we have

(71) $(\widehat{x}_{0},B_{L}^{-1}\widehat{x}_{0})_{A}\lesssim H^{2}.$

Combining Eqs. 69, 70 and 71 with $\lambda_{0}-\lambda_{1}\lesssim H^{2}$, we
obtain

$\rho_{0}-\lambda_{1}=f(x_{0})-\lambda_{1}\lesssim\lambda_{0}+H^{2}-\lambda_{1}\lesssim
H^{2},$

which finishes the proof.

######  Theorem 4.11.

Assume the initial point $\widehat{x}_{0}$ is the eigenvector calculated on
the coarse space $V_{0}$, i.e., $\widehat{x}_{0}=R_{0}^{\mathsf{T}}u_{0}$,
when $H$ is sufficiently small and

(72)
$H\kappa_{\nu}^{4}\Bigl{(}\frac{\lambda_{2}}{\lambda_{2}-\lambda_{1}}\Bigr{)}^{3}\lesssim
1\quad\text{and}\quad
H\frac{\kappa_{\nu}^{2}\lambda_{2}}{\lambda_{2}-\lambda_{1}}\leq\frac{1}{2},$

then Algorithm 2 is monotonic convergence, and the asymptotic convergence rate
is

$f(x_{m})-\lambda_{1}\leq
2\biggl{(}1-C\Bigl{(}\frac{\lambda_{2}-\lambda_{1}}{\kappa_{\nu}\lambda_{2}}\Bigr{)}^{1/2}\biggr{)}^{m}\bigl{(}f(x_{0})-\lambda_{1}\bigr{)},$

where $C$ is a constant independent of $h$, $H$ and
$(\lambda_{2}-\lambda_{1})/\lambda_{2}$.

######  Proof 4.12.

First, let us pick a $\rho_{\mathcal{X}}$ as

$\rho_{\mathcal{X}}=\lambda_{1}+H\kappa_{\nu}^{2}\lambda_{2}\leq\frac{\lambda_{1}+\lambda_{2}}{2}.$

By Corollaries 2.38 and 72, we know that

$\displaystyle\kappa_{B}$
$\displaystyle=\frac{L_{B}}{\mu_{B}}\lesssim\frac{\kappa_{\nu}\lambda_{2}}{\lambda_{2}-\lambda_{1}}+\frac{\kappa_{\nu}^{2}\lambda_{2}^{3}}{(\lambda_{2}-\lambda_{1})^{3}}\Bigl{(}\frac{\rho_{\mathcal{X}}-\lambda_{1}}{\lambda_{2}-\lambda_{1}}\kappa_{\nu}+H^{2}\Bigr{)}$
$\displaystyle\lesssim\frac{\kappa_{\nu}\lambda_{2}}{\lambda_{2}-\lambda_{1}}+H\kappa_{\nu}^{5}\Bigl{(}\frac{\lambda_{2}}{\lambda_{2}-\lambda_{1}}\Bigr{)}^{4}\lesssim\frac{\kappa_{\nu}\lambda_{2}}{\lambda_{2}-\lambda_{1}}.$

In this case, by Lemma 4.9,

$f(x_{0})-\lambda_{1}\lesssim
H^{2}=\frac{(\rho_{\mathcal{X}}-\lambda_{1})H}{\kappa_{\nu}^{2}\lambda_{2}}\lesssim\frac{H\lambda_{2}}{(\lambda_{2}-\lambda_{1})^{2}}\frac{(\rho_{\mathcal{X}}-\lambda_{1})}{11\sqrt{10}(3+2\sqrt{\kappa_{B}})^{2}\kappa_{B}},$

which means when $H$ is small enough, all assumptions in Theorem 3.15 hold.
Then the theorem is obtained by Theorem 3.15.

######  Remark 4.13.

When applying the two-level overlapping or BPS domain decomposition
preconditioner, the requirement Eq. 72 is a mild condition due to Eq. 68.

######  Remark 4.14.

Compared to the Preconditioned Steepest Descent (PSD) method, the (asymptotic)
rate of convergence is improved

$\text{from}\quad
1-C_{\text{PSD}}\frac{\lambda_{2}-\lambda_{1}}{\kappa_{\nu}\lambda_{2}}\quad\text{to}\quad
1-C_{\text{RAP}}\Bigl{(}\frac{\lambda_{2}-\lambda_{1}}{\kappa_{\nu}\lambda_{2}}\Bigr{)}^{1/2}.$

Such an improvement significantly accelerates the convergence when
$\kappa_{\nu}$ is large or the spectral gap $\lambda_{2}-\lambda_{1}$ is
small.

##  5 Numerical experiments

In this section, we compare

  * •

RAP: Riemannian acceleration with preconditioning,

  * •

PSD: preconditioned steepest descent,

  * •

RA: Riemannian acceleration,

  * •

SD: steepest descent,

to observe the acceleration phenomenon. The 2D Laplacian eigenvalue problem

$-\Delta u=\lambda u$

on the unit square $\Omega=(0,1)^{2}$ with the Dirichlet boundary condition
$u=0$ on $\partial\Omega$ is computed. Our experiments are performed in Matlab
R2020b, and the stopping criteria is $\rho-\lambda^{h}\leq
10^{-10}\lambda^{h}$, where $\rho$ is the computed Rayleigh quotient, and
$\lambda^{h}$ is the smallest eigenvalue of the discrete Laplacian operator,
or the iteration number reaches $20000$.

The standard finite element method and the two-level overlapping domain
decomposition method [50] are applied for the discretization and as the
preconditioner. The coarse mesh size is set as $H=2^{-2}$, the overlapping
ratio is set as $1/2$ and the fine mesh sizes range from $2^{-3}$ to
$2^{-10}$. Here are some observations from the iteration numbers in Table 2:

  * •

Without momentum term and preconditioning, SD converges very slow. Even for
medium size problem, i.e., $h=2^{-7}$, SD can not converge in $20000$
iterations.

  * •

The momentum can accelerate the convergence, but for large scale problems, RA
still needs numerous iterations.

  * •

The iteration numbers required by RA and SD are $\mathcal{O}(\kappa^{1/2})$
and $\mathcal{O}(\kappa)$, respectively, since the condition number of
eigensystem without preconditioning is about $\kappa\approx h^{-2}$.

  * •

The preconditioning technique is very powerful in computing large scale
problems. When the two-level overlapping domain decomposition preconditioner
is applied, both RAP and PSD converge in very few iterations.

  * •

Compared with the PSD method, RAP reaches the convergence with only about half
iterations. Since the computational cost of one iteration for both methods are
close, i.e., solving one preconditioned linear system, the RAP is effective.

Table 2: Iteration numbers for RAP, PSD, RA and SD method. The “$\times$” means algorithm fails to converge in $20000$ iterations. $h$ | $2^{-3}$ | $2^{-4}$ | $2^{-5}$ | $2^{-6}$ | $2^{-7}$ | $2^{-8}$ | $2^{-9}$ | $2^{-10}$
---|---|---|---|---|---|---|---|---
RAP | $15$ | $15$ | $15$ | $15$ | $15$ | $15$ | $17$ | $16$
PSD | $25$ | $28$ | $30$ | $32$ | $34$ | $35$ | $35$ | $36$
RA | $59$ | $108$ | $221$ | $387$ | $821$ | $1692$ | $3285$ | $6444$
SD | $91$ | $413$ | $1732$ | $6984$ | $\times$ | $\times$ | $\times$ | $\times$

##  Appendix A Matlab code for Riemannian Acceleration with Preconditioning
(RAP)

Code 1: Matlab code for Riemannian Acceleration with Preconditioning (RAP)

⬇

1ction err = RAP(A,M,B,x0,mu,L,iterMax,tol)

2rq = @(x) x'*A*x/(x'*M*x);

3kappa = sqrt(mu/L); beta = 3*kappa/(2-4*kappa);

4alpha = (sqrt(beta^2+4*(beta+1)*kappa^2)-beta)/2;

5gamma = alpha*mu/(alpha+beta); c1 = alpha/(alpha+beta+1);

6c2 = (1-alpha)/alpha; c3 = alpha/(1+beta)/gamma;

7x = B(x0); eta = sqrt(x'*x0); x = x/eta; x0 = x0/eta; v = x; v0 = x0;

8lambda = zeros(1,iterMax); lambda(1) = rq(x);

9for iter = 2:iterMax

10 % update y

11 eta = (x'*v0+v'*x0)/2; w = v-x*eta;

12 if norm(w) < eps % x=v=y

13 sigma = rq(x); g0 = (A*x-sigma*(M*x))*2/(x'*M*x); g = B(g0);

14 eta = c3*sqrt(g'*g0); g = g/eta; g0 = g0/eta;

15 v = v*cos(eta)+g*sin(eta); v0 = v0*cos(eta)+g0*sin(eta);

16 [U,R] = qr([x,g],0); [c,$\neg$] = eigs(U'*A*U,U'*M*U,1,'smallestabs');

17 x = U*c; c = R\c; x0 = x0*c(1)+g0*c(2);

18 else

19 w0 = v0-x0*eta; eta = sqrt(w'*w0); w = w/eta; w0 = w0/eta;

20 y = x*cos(theta)+w*sin(theta); y0 = x0*cos(theta)+w0*sin(theta);

21 % update v

22 eta = (y'*v0+v'*y0)/2; p = v-y*eta; p0 = v0-y0*eta;

23 eta = sqrt(p'*p0); p = p/eta; p0 = p0/eta;

24 sigma = rq(y); g0 = (A*y-sigma*(M*y))*2/(y'*M*y); g = B(g0);

25 q = c2*theta*p-c3*g; q0 = c2*theta*p0-c3*g0;

26 eta = sqrt(q'*q0); q = q/eta; q0 = q0/eta;

27 v = y*cos(eta)+q*sin(eta); v0 = y0*cos(eta)+q0*sin(eta);

28 % update x

29 [U,R] = qr([x,y,g],0);

30 [c,$\neg$] = eigs(U'*A*U,U'*M*U,1,'smallestabs');

31 x = U*c; c = R\c; x0 = x0*c(1)+y0*c(2)+g0*c(3);

32 end

33 eta = sqrt(x'*x0); x = x/eta; x0 = x0/eta; lambda(iter) = rq(x);

34end

35end

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds, Princeton University Press, Princeton, 2009, https://doi.org/10.1515/9781400830244.
  * [2] K. Ahn and S. Sra, From Nesterov’s estimate sequence to Riemannian acceleration, in Conference on Learning Theory, vol. 125, 2020, pp. 84–118, https://proceedings.mlr.press/v125/ahn20a.html.
  * [3] F. Alimisis, A. Orvieto, G. Becigneul, and A. Lucchi, A continuous-time perspective for modeling acceleration in Riemannian optimization, in International Conference on Artificial Intelligence and Statistics, vol. 108, 2020, pp. 1297–1307, https://proceedings.mlr.press/v108/alimisis20a.html.
  * [4] F. Alimisis, A. Orvieto, G. Becigneul, and A. Lucchi, Momentum improves optimization on Riemannian manifolds, in International Conference on Artificial Intelligence and Statistics, 2021, pp. 1351–1359, https://proceedings.mlr.press/v130/alimisis21a.html.
  * [5] F. Alimisis and B. Vandereycken, Geodesic convexity of the symmetric eigenvalue problem and convergence of Riemannian steepest descent, arXiv preprint arXiv:2209.03480, (2022), https://doi.org/10.48550/arXiv.2209.03480.
  * [6] M. E. Argentati, A. V. Knyazev, K. Neymeyr, E. E. Ovtchinnikov, and M. Zhou, Convergence theory for preconditioned eigenvalue solvers in a nutshell, Foundations of Computational Mathematics, 17 (2017), pp. 713–727, https://doi.org/10.1007/s10208-015-9297-1.
  * [7] I. Babuška and J. Osborn, Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems, Mathematics of Computation, 52 (1989), pp. 275–297, https://doi.org/10.1090/S0025-5718-1989-0962210-8.
  * [8] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, eds., Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia, 2000, https://doi.org/10.1137/1.9780898719581.
  * [9] A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences, 2 (2009), pp. 183–202, https://doi.org/10.1137/080716542.
  * [10] S. Becker, J. Bobin, and E. J. Candès, NESTA: A fast and accurate first-order method for sparse recovery, SIAM Journal on Imaging Sciences, 4 (2011), pp. 1–39, https://doi.org/10.1137/090756855.
  * [11] N. Boumal, An introduction to optimization on smooth manifolds, Cambridge University Press, Cambridge, 2023, https://doi.org/10.1017/9781009166164.
  * [12] J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring. I, Mathematics of Computation, 47 (1986), pp. 103–134, https://doi.org/10.1090/S0025-5718-1986-0842125-3.
  * [13] W. Chen, N. Shao, and X. Xu, A locally optimal preconditioned Newton-Schur method for symmetric elliptic eigenvalue problems, To appear in Mathematics of Computation, (2023), https://doi.org/10.1090/mcom/3860.
  * [14] X. Dai, Q. Wang, and A. Zhou, Gradient flow based Kohn–Sham density functional theory model, Multiscale Modeling & Simulation, 18 (2020), pp. 1621–1663, https://doi.org/10.1137/19M1276170.
  * [15] J. W. Demmel, Applied Numerical Linear Algebra, SIAM Philadelphia, 1997, https://doi.org/10.1137/1.9781611971446.
  * [16] M. Dryja, An additive Schwarz algorithm for two-and three-dimensional finite element elliptic problems, in Domain Decomposition Methods. Second International Symposium on Domain Decomposition Methods,, SIAM Philadelphia, 1989, pp. 168–172.
  * [17] V. Duruisseaux and M. Leok, A variational formulation of accelerated optimization on Riemannian manifolds, SIAM Journal on Mathematics of Data Science, 4 (2022), pp. 649–674, https://doi.org/10.1137/21M1395648.
  * [18] E. G. D’yakonov, Optimization in solving elliptic problems, CRC Press, Boca Raton, 1996, https://doi.org/10.1201/9781351075213.
  * [19] A. Edelman, T. A. Arias, and S. T. Smith, The geometry of algorithms with orthogonality constraints, SIAM journal on Matrix Analysis and Applications, 20 (1998), pp. 303–353, https://doi.org/10.1137/S0895479895290954.
  * [20] B. Jiang and Y.-H. Dai, A framework of constraint preserving update schemes for optimization on Stiefel manifold, Mathematical Programming, 153 (2015), pp. 535–575, https://doi.org/10.1007/s10107-014-0816-7.
  * [21] J. Kim and I. Yang, Accelerated gradient methods for geodesically convex optimization: Tractable algorithms and convergence analysis, in International Conference on Machine Learning, PMLR, 2022, pp. 11255–11282, https://proceedings.mlr.press/v162/kim22k.html.
  * [22] W. P. Klingenberg, Riemannian Geometry, De Gruyter, Berlin, New York, 1995, https://doi.org/doi:10.1515/9783110905120.
  * [23] A. V. Knyazev, Preconditioned eigensolvers - an oxymoron?, Electronic Transactions on Numerical Analysis, 7 (1998), pp. 104–123.
  * [24] A. V. Knyazev, Toward the optimal preconditioned eigensolver: locally optimal block preconditioned conjugate gradient method, SIAM Journal on Scientific Computing, 23 (2001), pp. 517–541, https://doi.org/10.1137/S1064827500366124.
  * [25] A. V. Knyazev and K. Neymeyr, A geometric theory for preconditioned inverse iteration III: A short and sharp convergence estimate for generalized eigenvalue problems, Linear Algebra and its Applications, 358 (2003), pp. 95–114, https://doi.org/10.1016/S0024-3795(01)00461-X.
  * [26] A. V. Knyazev and K. Neymeyr, Gradient flow approach to geometric convergence analysis of preconditioned eigensolvers, SIAM Journal on Matrix Analysis and Applications, 31 (2009), pp. 621–628, https://doi.org/10.1137/080727567.
  * [27] J. M. Lee, Introduction to Riemannian manifolds, vol. 176, Springer, Cham, 2018, https://doi.org/10.1007/978-3-319-91755-9.
  * [28] Y. Liu, F. Shang, J. Cheng, H. Cheng, and L. Jiao, Accelerated first-order methods for geodesically convex optimization on Riemannian manifolds, in Advances in Neural Information Processing Systems, vol. 30, 2017, pp. 4868–4877, https://proceedings.neurips.cc/paper/2017/file/6ef80bb237adf4b6f77d0700e1255907-Paper.pdf.
  * [29] H. Luo and L. Chen, From differential equation solvers to accelerated first-order methods for convex optimization, Mathematical Programming, (2021), pp. 1–47, https://doi.org/10.1007/s10107-021-01713-3.
  * [30] D. Martínez-Rubio, Global Riemannian acceleration in hyperbolic and spherical spaces, in International Conference on Algorithmic Learning Theory, 2022, pp. 768–826, https://proceedings.mlr.press/v167/martinez-rubio22a.html.
  * [31] B. Mishra and R. Sepulchre, Riemannian preconditioning, SIAM Journal on Optimization, 26 (2016), pp. 635–660, https://doi.org/10.1137/140970860.
  * [32] M. Muehlebach and M. Jordan, A dynamical systems perspective on Nesterov acceleration, in International Conference on Machine Learning, 2019, pp. 4656–4662, http://proceedings.mlr.press/v97/muehlebach19a/muehlebach19a.pdf.
  * [33] M. Muehlebach and M. Jordan, Optimization with momentum: Dynamical, control-theoretic, and symplectic perspectives., Journal of Machine Learning Research, 22 (2021), pp. 1–50, http://jmlr.org/papers/v22/20-207.html.
  * [34] Y. Nesterov, A method for solving the convex programming problem with convergence rate $\mathcal{O}(1/k^{2})$, Soviet Mathematics Doklady, 269 (1983), pp. 543–547.
  * [35] Y. Nesterov, Smooth minimization of non-smooth functions, Mathematical Programming, 103 (2005), pp. 127–152, https://doi.org/10.1007/s10107-004-0552-5.
  * [36] K. Neymeyr, A geometric convergence theory for the preconditioned steepest descent iteration, SIAM Journal on Numerical Analysis, 50 (2012), pp. 3188–3207, https://doi.org/10.1137/11084488X.
  * [37] K. Neymeyr, E. Ovtchinnikov, and M. Zhou, Convergence analysis of gradient iterations for the symmetric eigenvalue problem, SIAM Journal on Matrix Analysis and Applications, 32 (2011), pp. 443–456, https://doi.org/10.1137/100784928.
  * [38] Y. Notay, Combination of Jacobi–Davidson and conjugate gradients for the partial symmetric eigenproblem, Numerical Linear Algebra with Applications, 9 (2002), pp. 21–44, https://doi.org/10.1002/nla.246.
  * [39] E. E. Ovtchinnikov, Sharp convergence estimates for the preconditioned steepest descent method for Hermitian eigenvalue problems, SIAM Journal on Numerical Analysis, 43 (2006), pp. 2668–2689, https://doi.org/10.1137/040620643.
  * [40] E. E. Ovtchinnikov, Jacobi correction equation, line search, and conjugate gradients in Hermitian eigenvalue computation I: Computing an extreme eigenvalue, SIAM Journal on Numerical Analysis, 46 (2008), pp. 2567–2592, https://doi.org/10.1137/070688742.
  * [41] B. O’donoghue and E. Candes, Adaptive restart for accelerated gradient schemes, Foundations of Computational Mathematics, 15 (2015), pp. 715–732, https://doi.org/10.1007/s10208-013-9150-3.
  * [42] J.-H. Park, A. J. Salgado, and S. M. Wise, Preconditioned accelerated gradient descent methods for locally lipschitz smooth objectives with applications to the solution of nonlinear pdes, Journal of Scientific Computing, 89 (2021), pp. 1–37, https://doi.org/10.1007/s10915-021-01615-8.
  * [43] B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, Philadelphia, 1998, https://doi.org/10.1137/1.9781611971163.
  * [44] Y. Saad, Numerical methods for large eigenvalue problems: revised edition, SIAM, Philadelphia, 2011, https://doi.org/10.1137/1.9781611970739.
  * [45] N. Shao and W. Chen, Convergence analysis of Newton–Schur method for symmetric elliptic eigenvalue problem, SIAM Journal on Numerical Analysis, 61 (2023), pp. 315–342, https://doi.org/10.1137/21M1448847.
  * [46] N. Shao, W. Chen, and Z. Bai, A symmetric eigensolver based on preconditioning and implicit convexity (EPIC).  In preparation, 2023.
  * [47] B. Shi, S. S. Du, M. Jordan, and W. J. Su, Understanding the acceleration phenomenon via high-resolution differential equations, Mathematical Programming, (2021), pp. 1–70, https://doi.org/10.1007/s10107-021-01681-8.
  * [48] B. Shi, S. S. Du, W. Su, and M. Jordan, Acceleration via symplectic discretization of high-resolution differential equations, Advances in Neural Information Processing Systems, 32 (2019), https://proceedings.neurips.cc/paper/2019/file/a9986cb066812f440bc2bb6e3c13696c-Paper.pdf.
  * [49] W. Su, S. Boyd, and E. Candes, A differential equation for modeling Nesterov’s accelerated gradient method: Theory and insights, in Advances in Neural Information Processing Systems, vol. 27, 2014, https://proceedings.neurips.cc/paper/2014/file/f09696910bdd874a99cd74c8f05b5c44-Paper.pdf.
  * [50] A. Toselli and O. Widlund, Domain Decomposition Methods - Algorithms and Theory, Springer, Berlin Heidelberg, 2005, https://doi.org/10.1007/b137868.
  * [51] N. J. Walkington, Nesterov’s method for convex optimization, SIAM Review, 65 (2023), pp. 539–562, https://doi.org/10.1137/21M1390037, https://doi.org/10.1137/21M1390037.
  * [52] W. Wang and X. Xu, A two-level overlapping hybrid domain decomposition method for eigenvalue problems, SIAM Journal on Numerical Analysis, 56 (2018), pp. 344–368, https://doi.org/10.1137/16M1088302.
  * [53] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965, https://doi.org/10.1017/S0013091500012104.
  * [54] H. Zhang and S. Sra, An estimate sequence for geodesically convex optimization, in Conference On Learning Theory, vol. 75, 2018, pp. 1703–1723, https://proceedings.mlr.press/v75/zhang18a.html.