# Convergence of Oja’s online principal component flow

Jian-Guo Liu  Department of Mathematics and Department of Physics, Duke
University, Durham, NC<EMAIL_ADDRESS>and  Zibu Liu  Department of
Mathematics, Duke University, Durham, NC<EMAIL_ADDRESS>
###### Abstract.

Online principal component analysis (PCA) has been an efficient tool in
practice to reduce dimension. However, convergence properties of the
corresponding ODE are still unknown, including global convergence, stable
manifolds, and convergence rate. In this paper, we focus on the stochastic
gradient ascent (SGA) method proposed by Oja. By regarding the corresponding
ODE as a Landau-Lifshitz-Gilbert (LLG) equation on the Stiefel manifold, we
proved global convergence of the ODE. Moreover, we developed a new technique
to determine stable manifolds. This technique analyzes the rank of the initial
datum. Using this technique, we derived the explicit expression of the stable
manifolds. As a consequence, exponential convergence to stable equilibrium
points was also proved. The success of this new technique should be attributed
to the semi-decoupling property of the SGA method: iteration of previous
components does not depend on that of later ones. As far as we know, our
result is the first complete one on the convergence of an online PCA flow,
providing global convergence, explicit characterization of stable manifolds,
and closed formula of exponential convergence depending on the spectrum gap.

###### Key words and phrases:

machine learning, dimensionality reduction, online principal component
analysis, gradient flow, stable manifold

##  1\. Introduction

The recursive method for principal component analysis (PCA) proposed by Oja in
[12, 13] pioneered the discussion of several online PCA algorithms and their
applications in the field of machine learning. Suppose that
$\mathbf{x}\in\mathbb{R}^{n}$ is a random variable with zero expectation and
distribution $\nu$. PCA aims to find a group of orthonormal vectors
$\mathbf{w_{1}},\ \mathbf{w_{2}},\ ,...,\mathbf{w_{p}}\in\mathbb{R}^{n},1\leq
p\leq n$ such that they maximize

(1.1) $\displaystyle
R_{p}:=\mathbb{E}\left[\sum_{i=1}^{p}(\mathbf{w_{i}}^{T}\mathbf{x})^{2}\right].$

Equivalently, $\mathbf{w_{1}},\ \mathbf{w_{2}},\ ...,\ \mathbf{w_{p}}$ are the
$p$ dominant eigenvectors of the covariance matrix

(1.2)
$\displaystyle\mathbf{A}:=\mathbb{E}\left[\mathbf{x}\mathbf{x}^{T}\right].$

If $\mathbf{A}$ is given, then PCA reduces to diagonalizing a symmetric matrix
which is a standard linear algebra problem. Nevertheless, $\mathbf{A}$ is
generally unknown in practice. Instead, only a series of independent samplings
of the random variable $\mathbf{x}$ are available. Moreover, due to the
limitation of storage, only several or a limited amount of samplings can be
addressed at the same time. Therefore, an online version of the PCA algorithm
is desired: it only requires several most recent samplings of $\mathbf{x}$ to
complete the iteration and efficiently converges to eigenvectors. The term
’online’ here means that only a limited amount of samplings of $\mathbf{x}$
are available in one iteration.

To solve this problem, Oja [12, 13] derived the following online PCA algorithm
which approximates dominant eigenvectors of $\mathbf{A}$ efficiently. Let
$\mathbf{W}=[\mathbf{w}_{1},\ \mathbf{w_{2}},...,\
\mathbf{w_{p}}]\in\mathbb{R}^{n\times p}$. Given a stream of data
$\\{\mathbf{x}(k)\\}$ which are independent samples of the random variable
$\mathbf{x}$, i.e., $\mathbf{x}(k)\sim\nu$ are independent, consider
$\mathbf{W}(k),k=1,2,...$ which are defined as:

(1.3) $\displaystyle\left\\{\begin{split}\mathbf{W}(0)&\in O(n\times p),\\\
\tilde{\mathbf{W}}(k)&=\mathbf{W}(k-1)+\eta_{k}\mathbf{x}(k)\mathbf{x}^{T}(k)\mathbf{W}(k-1),\
k\geq 1,\\\ \mathbf{W}(k)&=\tilde{\mathbf{W}}(k)\mathbf{S}(k),\ k\geq
1.\end{split}\right.$

Here $\eta_{k},k=1,2,...$ are real numbers representing learning rates and
$\mathbf{S}(k)\in\mathbb{R}^{p\times p}$ are matrices depending on
$\tilde{\mathbf{W}}(k)$ which orthonormalize $\tilde{\mathbf{W}}(k)$.
Therefore, $\mathbf{W}^{T}(k)\mathbf{W}(k)=\mathbf{I_{p}}$ holds for all $k$,
where $\mathbf{I_{p}}$ is the identity matrix of size $p$.

In (1.3), if $\mathbf{S}(k)$ conducts the Gram-Schmidt orthonormalization
(GSO) on the columns of $\tilde{\mathbf{W}}(k)$, then the stochastic gradient
ascent (SGA) algorithm is derived [12]. To reduce the computational complexity
of the SGA method, authors of [11, 13] developed a first order approximation
of the SGA scheme:

(1.4) $\displaystyle\mathbf{w_{j}}(k)$
$\displaystyle=\mathbf{w_{j}}(k-1)+\eta(k)\mathbf{x}^{T}(k)\mathbf{w_{j}}(k-1)[\mathbf{x}(k)-(\mathbf{x}^{T}(k)\mathbf{w_{j}}(k-1))\mathbf{w_{j}}(k-1)$
$\displaystyle-2\sum_{i=1}^{j-1}(\mathbf{x}^{T}(k)\mathbf{w_{i}}(k-1))\mathbf{w_{i}}(k-1)],\
j=1,\ 2,\ ...,\ p.$

Because higher-order terms are omitted, this scheme does not preserve
orthogonality. This approximated scheme is easy to apply in simulations and
appropriate for neural network implementation [12]. The reason is that
iteration of $\mathbf{w_{j}}$ only depends on $\mathbf{w_{i}},i=1,2,...,j-1$.
We will call this property the ’semi-decoupling’ property.

From now on, we assume $p=n$, i.e., we consider extracting all eigenvectors
instead of only the first several dominant ones. Due to the semi-decoupling
property, results for $p<n$ are actually a consequence of the one for $p=n$.
This will be explained in Section 5.4.

###  1.1. The corresponding differential equation for the SGA method

Formally omitting the higher-order terms and replacing
$\mathbf{x}(k)\mathbf{x}^{T}(k)$ by its average $\mathbf{A}$ in (1.4) as in
[12], one can derive the following differential equation of
$\mathbf{q_{i}}\in\mathbb{R}^{n},\ i=1,2,...,n$:

(1.5)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{q_{1}}}&=\mathbf{Aq_{1}}-\mathbf{(q_{1}\cdot
Aq_{1})q_{1}},\\\ \mathbf{\dot{q_{j}}}&=\mathbf{Aq_{j}}-\mathbf{(q_{j}\cdot
Aq_{j})q_{j}}-2\sum_{i=1}^{j-1}\mathbf{(q_{i}\cdot Aq_{j})q_{i}},\ j=2,\ 3,\
...,\ n.\\\ \mathbf{q_{i}}(0)&=\mathbf{q_{i,0}},\ i=1,\ 2,\ ...,\
n.\end{split}\right.$

Here
$\mathbf{Q_{0}}:=[\mathbf{q_{1,0}},\mathbf{q_{2,0}},...,\mathbf{q_{n,0}}]\in
O(n)$. One can rewrite (1.5) in matrices by denoting
$\mathbf{Q}(t)=[\mathbf{q_{1}}(t),\mathbf{q_{2}}(t),...,\mathbf{q_{n}}(t)]$,
which yields

(1.6)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{AQ}-\mathbf{Q}\sum_{j=1}^{n}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQ}\mathbf{E_{j}}-2\mathbf{Q}\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n).\end{split}\right.$

Here $\mathbf{E_{j}},j=1,2,...,n$ represent the matrix with 1 at $j$th row,
$j$th column, but 0 at other positions.

System (1.6) admits a unique solution by standard ODE theories. In fact, one
can prove that system (1.6) defines a flow on the Stiefel manifold $O(n)$:

######  Lemma 1.1.

Suppose that $\mathbf{Q}(t)$ is the solution to the system (1.6). Then for any
initial value $\mathbf{Q_{0}}\in O(n)$, we have $\mathbf{Q}(t)\in O(n)$, i.e.

(1.7)
$\displaystyle\mathbf{Q}(t)\mathbf{Q}^{T}(t)=\mathbf{Q}^{T}(t)\mathbf{Q}(t)=\mathbf{I_{n}}.$

Here $\mathbf{I_{n}}\in\mathbb{R}^{n\times n}$ is the identity matrix of size
$n$.

See proof of Lemma 1.1 in Section 2. In this paper, we will provide a complete
convergence analysis of (1.6), including global convergence, stable manifolds,
and convergence rate. In the following part of the introduction, we will first
review literature and then highlight our contribution to the research of
online PCA.

###  1.2. Previous results and unsolved problems

Another PCA flow developed by Oja and Brocket [4] also attracted much
attention:

(1.8)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{AQD}-\mathbf{QDQ}^{T}\mathbf{AQ},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n).\end{split}\right.$

where $\mathbf{D}\in\mathbb{R}^{n\times n}$ is a symmetric matrix. As (1.6),
(1.8) is also invariant on the Stiefel manifold. Global convergence to
equilibria of (1.8) was derived first in [7] under the assumption that
$\mathbf{D}$ is diagonal. Later on, if $\mathbf{D}$ is only symmetric, global
convergence was again verified in [16]. However, to our best knowledge,
similar result on the flow (1.6) is unknown. Only local convergence result
[13] and the global convergence of the first component $\mathbf{p_{1}}$ [10]
are available. The barrier of a global analysis of (1.6) is probably due to
its complicated matrix representation, unlike (1.8) which is clean.

Some other convergence properties of both (1.8) and (1.6) also remain unknown.
First, the stable manifolds of both (1.8) and (1.6) (if exists) are
undetermined [3]. Second, a global convergence rate was not computed for
either (1.8) or (1.6). Only local asymptotic analysis was conducted on both
(1.8) [7, 16] and (1.6) [11, 14]. In conclusion, only a few results were
derived for both (1.8) and (1.6).

###  1.3. Main results

First, we studied convergence property of (1.6). In analog to the Landau-
Lifshitz-Gilbert equation [8, 5, 6], we recast (1.6) in the following form on
a Riemannian manifold:

$\displaystyle\dfrac{\mathrm{d}\mathbf{m}(t)}{\mathrm{d}t}=\widehat{\dot{\mathbf{m}}}_{1}(\mathbf{m}(t),t)-\nabla_{g}\mathcal{E}(\mathbf{m}(t)).$

See (3.9) and Proposition 3.1. This equation consists of two parts: the
Hamiltonian part $\widehat{\dot{\mathbf{m}}}_{1}$ and the dissipative part
$-\nabla_{g}\mathcal{E}$, which is the variation of $-\mathcal{E}$ in some
sense. They are perpendicular under the Riemannian metric $g$. The critical
relationship between them is that the Hamiltonian conserved by the Hamiltonian
part is exactly the free energy $\mathcal{E}$ minimized by the dissipative
part.

To regard (1.6) as an LLG equation on the Stiefel manifold, we specified the
free energy $\mathcal{E}$ and the Riemmanian metric $g$. Instead of directly
using the Rayleigh quotient, we adopted a weighted version of it:

$\displaystyle\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})=\mathrm{tr}(\mathbf{N}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}).$

See (4.2). Here $\mathbf{N}$ is a given diagonal matrix with entries on the
diagonal line aligned in a descending order. The critical advantage of
utilizing $\mathbf{N}$ here is that it helps align the eigenbasis in order. In
fact, by the Wielandt-Hoffman inequality (see Lemma 4.1),
$\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ is maximized on $O(n)$ if and
only if $\mathbf{Q}$ is exactly the eigenbasis (up to sign) of $\mathbf{A}$
aligned in a descending order according to the eigenvalues of $\mathbf{A}$.
This is exactly the result one desires from conducting PCA. Therefore, we
selected $-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ as the free energy
(minimizing $-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ is equivalent to
maximizing $\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$).

Observing that $-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ is indeed a
Lyapunov function of (1.6) (see Lemma 4.2), we constructed the Reimannian
metric $\tilde{g}$ by the method in [1]. For two tangent vector fields
$\mathbf{X}$ and $\mathbf{Y}$ on the Stiefel manifold, $\tilde{g}$ is defined
as

$\displaystyle\langle\mathbf{X},\mathbf{Y}\rangle_{\tilde{g}}$
$\displaystyle:=\langle\mathbf{X_{0}},\mathbf{Y_{0}}\rangle_{g_{e}}-\dfrac{1}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\langle\mathcal{E}^{\prime},\
\mathbf{X_{1}}\rangle_{F}\langle\mathcal{E}^{\prime},\
\mathbf{Y_{1}}\rangle_{F}.$

Here $g_{e}$ is the Eucledian metric, $\mathbf{F}$ is defined in (2.9),
$\mathcal{E}^{\prime}$ is the derivative of the $\mathcal{E}$ w.r.t.
$\mathbf{Q}$, $\mathbf{X}=\mathbf{X_{0}}+\mathbf{X_{1}}$ and
$\mathbf{Y}=\mathbf{Y_{0}}+\mathbf{Y_{1}}$ are certain decompositions of
$\mathbf{X}$ and $\mathbf{Y}$ respectively, see (4.19). For details of
$\tilde{g}$, see (4.20) and Proposition 4.2.

Finally, we can reformulate (1.6) as an LLG equation on the Stiefel manifold:
in (3.9), let
$\widehat{\dot{\mathbf{m}}}_{1}=\mathbf{0},\mathcal{E}=-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N}),g=\tilde{g}$.
We emphasize that the Hamiltonian part $\widehat{\dot{\mathbf{m}}}_{1}$ is
degenerate in (1.6). Thus, we will explore usage of this Hamiltonian part in
future.

After the above reformulation, we are ready to study the convergence property
of (1.6). First, in Theorem 5.1, we proved the global convergence of (1.6):
for arbitrary initial datum in $O(n)$, (1.6) converges to an equilibrium of
the flow.

Moreover, we comprehensively characterized the stable manifolds of (1.6). To
our best knowledge, this is the first complete result on stable manifolds for
a PCA flow. To derive that, we developed a new technique based on rank
analysis on the initial datum $\mathbf{Q}_{0}$. For invertible
$\mathbf{Q}\in\mathbb{R}^{n\times n}$, we recursively defined
$\sigma_{i}(\mathbf{Q}),\ i=1,2,...,n$ which is a permutation of
$\\{1,2,...,n\\}$:

$\displaystyle\sigma_{1}(\mathbf{Q}):=\inf\limits_{1\leq k\leq n}\\{k:\
q_{1,k}\neq 0\\},$ $\displaystyle\sigma_{m}(\mathbf{Q}):=\inf\limits_{1\leq
k\leq n}\\{k:\ \mathrm{det}(\mathbf{Q}[1,2,3,...,m;\ \sigma(1),\ \sigma(2),\
...,\ k])\neq 0\\}.$

See (5.9). The geometric meaning of $\sigma_{m}(\mathbf{Q})$ is to determine
which component of $\mathbf{Q}$ converges to the $m$-th eigenvector of
$\mathbf{A}$.

By Theorem 5.2, the convergent point of $\mathbf{Q}(t)$ in (1.6) with initial
value $\mathbf{Q}(0)=\mathbf{Q_{0}}$ is given by

$\displaystyle\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}(t)\to\mathrm{sgn}(z_{m})\mathbf{e_{m}},\
m=1,2,...,n,$

here $z_{m}$ is a number defined in Lemma 5.3. Geometrically, $z_{m}$
determines the orientation of the orthogonal frame. One can see that the
convergence point of $\mathbf{Q}(t)$ as $t\to\infty$ only depends on
$\sigma_{m}(\mathbf{Q_{0}}),m=1,2,...,n$, so by Corollary 5.1, we have
$\sigma_{m}(\mathbf{Q}(t))=\sigma_{m}(\mathbf{Q_{0}}),m=1,2,...,n$ for all
$t>0$, thus rank analysis is only necessary for the initial datum.

As a corollary, we also determined the convergence rate to stable equilibria.
If $\mathbf{Q}(t)$ converges to the
$[\mathbf{e_{1}},\mathbf{e_{2}},...,\mathbf{e_{n}}]$ (up to sign), then by
Theorem 5.3,

$\|\mathbf{q_{i}}-\mathbf{e_{i}}\|\leq Ce^{-\nu_{i}t},\ i=1,2,...,n,$

where

$\displaystyle\nu_{k}$ $\displaystyle:=\min\\{\lambda_{1}-\lambda_{2},\
\lambda_{2}-\lambda_{3},\ ...,\ \lambda_{k}-\lambda_{k+1}\\}>0,\
k=1,2,...,n-1,$ $\displaystyle\nu_{n}$ $\displaystyle:=\nu_{n-1}$

and $\lambda_{1}>\lambda_{2}>...>\lambda_{n}$ are eigenvalues of $\mathbf{A}$.
See (5.29) for details.

Theorem 5.2 and 5.3 reveals the mechanism of the SGA method. If $\mathbf{Q}$
converges to a stable equilibrium, because $v_{k}\geq v_{k+1}$, convergence of
$\mathbf{q_{k}}$ is faster than that of $\mathbf{q_{k+1}}$ for all
$k=1,2,...,n-1$. Thus, alignment of $\mathbf{q_{1}}$ to $\mathbf{e_{1}}$ is
the first to complete, and due to orthogonality,
$\mathbf{q_{2}},...,\mathbf{q_{n}}$ will be forced into the orthogonal
complement of $\mathrm{span}\\{\mathbf{q_{1}}\\}$, which is very close to the
subspace of $\mathrm{span}\\{\mathbf{e_{2}},...,\mathbf{e_{n}}\\}$. Then,
alignment of $\mathbf{q_{2}}$ to $\mathbf{e_{2}}$ will be completed and
$\mathbf{q_{3}},...,\mathbf{q_{n}}$ will be forced into the orthogonal
complement of $\mathrm{span}\\{\mathbf{q_{1}},\mathbf{q_{2}}\\}$, which is
close to the subspace of
$\mathrm{span}\\{\mathbf{e_{3}},...,\mathbf{e_{n}}\\}$. Meanwhile, in this
process, alignment of $\mathbf{q_{1}}$ will not be influenced due to semi-
decoupling. Alignment of other components is then conducted in order. Finally,
alignment of $\mathbf{q_{n}}$ is completed and the solution converges. This is
the real mechanism of the SGA method in the case where $\mathbf{A}$ has all
single eigenvalues.

Finally, results for the case $p<n$ is discussed in Section 5.4. By the semi-
decoupling property, all results in the case of $p=n$ can be extended to this
case.

The rest of paper will be organized in the following pattern: in Section 2, we
clarify assumptions and notations. In Section 3, we introduced the Landau-
Lifshitz-Gilbert equation and generalize it into Riemannian manifolds. In
Section 4 and 5, we first regard (1.6) as a Landau-Lifshitz-Gilbert equation
on the Stiefel manifold. Then we derive global convergence of (1.6) by using
its dissipative property. Moreover, we develop other convergence properties,
including the stable manifolds and the convergence rate by exploiting the
semi-decoupling property. This critical property is reinterpreted as a
solution formula of (1.6) which involves the Cholesky decomposition. Omitted
proofs of lemmas and calculation can be found in Section 6.

##  2\. Premier

###  2.1. Reformulation of (1.6)

Before we introduce our main results in detail, we conduct preliminary
analysis of (1.6) and recast it into a more compact form.

Notice the following identity: for all $\mathbf{B}\in\mathbb{R}^{n\times n}$,

$\displaystyle\sum_{1\leq k\leq j\leq n}\mathbf{E_{k}BE_{j}}+\sum_{1\leq
j<k\leq n}\mathbf{E_{k}BE_{j}}=\sum_{1\leq j,k\leq
n}\mathbf{E_{k}BE_{j}}=\mathbf{B},$

we have

$\displaystyle\ \ \ \
\sum_{j=1}^{n}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQ}\mathbf{E_{j}}+2\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}$
$\displaystyle=\sum_{j=1}^{n}\sum_{k=1}^{j}\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}+\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}$
$\displaystyle=\mathbf{Q}^{T}\mathbf{AQ}+\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}-\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}},$

thus (1.6) can be reformulated as

(2.1)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{AQ}-\mathbf{Q}\mathbf{Q}^{T}\mathbf{AQ}+\mathbf{Q}\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}:=[\mathbf{q_{1,0}},\mathbf{q_{2,0}},...,\mathbf{q_{n,0}}]\in
O(n).\end{split}\right.$

Similarly, we can rewrite (1.4) by matrices. For
$\mathbf{\Lambda},\mathbf{Q}\in\mathbb{R}^{n\times n}$, define

(2.2) $\displaystyle\mathbf{\Sigma}(\mathbf{\Lambda},\mathbf{Q})$
$\displaystyle:=\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{\Lambda
QE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{\Lambda QE_{j}},$
$\displaystyle\mathbf{G}(\mathbf{\Lambda},\mathbf{Q})$
$\displaystyle:=\mathbf{\Lambda
Q}-\mathbf{QQ}^{T}\mathbf{\Lambda}\mathbf{Q}+\mathbf{Q}\mathbf{\Sigma}(\mathbf{\Lambda},\mathbf{Q}).$

Then (1.4) also reads as

(2.3)
$\displaystyle\left\\{\begin{split}\mathbf{A}(k)&=\mathbf{x}(k)\mathbf{x}^{T}(k),\\\
\mathbf{W}(k)&=\mathbf{W}(k-1)+\eta_{k}\mathbf{G}(\mathbf{A}(k),\mathbf{W}(k-1)).\end{split}\right.$

Remember that by Lemma 1.1, solutions to (1.6) satisfy
$\mathbf{Q}(t)\mathbf{Q}^{T}(t)=\mathbf{I_{n}},t>0$ if $\mathbf{Q}(0)\in
O(n)$. We prove this lemma here.

###### Proof of Lemma 1.1.

Let
$\mathbf{F}(\mathbf{Q})=\mathbf{Q}\displaystyle\sum_{j=1}^{n}\sum_{k=1}^{j-1}(\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}})$.
So $\mathbf{F}(\mathbf{Q})=\mathbf{Q}\mathbf{\Sigma({A,Q})}$. Because
$\mathbf{A}$ is symmetric, so $\mathbf{\Sigma({A,Q})}$ is a skew-symmetric
matrix. Direct computation implies that

(2.4)
$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}(\mathbf{Q}(t)\mathbf{Q}^{T}(t))$
$\displaystyle=\mathbf{A}\mathbf{Q}\mathbf{Q}^{T}+\mathbf{Q}\mathbf{Q}^{T}\mathbf{A}-2\mathbf{Q}\mathbf{Q}^{T}\mathbf{AQ}\mathbf{Q}^{T}+\mathbf{Q}\mathbf{\Sigma({A,Q})}\mathbf{Q}^{T}+\mathbf{Q}\mathbf{\Sigma({A,Q})}^{T}\mathbf{Q}^{T}.$

Because $\mathbf{\Sigma({A,Q})}$ is skew-symmetric, thus
$\mathbf{P}(t):=\mathbf{Q}(t)\mathbf{Q}^{T}(t)$ satisfies an algebraic Ricatti
eqaution [2]:

(2.5)
$\displaystyle\mathbf{\dot{P}}(t)=\mathbf{AP}+\mathbf{PA}-2\mathbf{PAP},\
\mathbf{P}(0)=\mathbf{I_{n}}.$

By uniqueness of solutions, we know that $\mathbf{P}(t)$ stays invariant, i.e.
$\mathbf{P}(t)=\mathbf{I_{n}}$ for any $t\geq 0$. So $\mathbf{Q}(t)\in O(n)$.
∎

######  Remark 2.1.

Notice that even though $\mathbf{Q}(0)$ does not stay on the Stiefel manifold,
$\mathbf{P}(t)=\mathbf{Q}(t)\mathbf{Q}^{T}(t)$ also satisfies (2.5) with
$\mathbf{P}(0)\neq\mathbf{I_{n}}$. Based on this observation, it is proved in
[15] that if rank$(\mathbf{Q}(0))=n$, then

(2.6)
$\displaystyle\lim\limits_{t\to\infty}\|\mathbf{Q}^{T}(t)\mathbf{Q}(t)-\mathbf{I_{n}}\|_{F}=0.$

This convergence is also exponential convergence. Main steps of this (2.6)
will be provided in Section 6.

So we can simplify (2.1) by cancelling the first two terms on the R.H.S.:

(2.7)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{Q}\sum_{j=1}^{n}\sum_{k=1}^{j-1}(\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}),\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n).\end{split}\right.$

Preserving the notation in the proof of Lemma 1.1 and notations in (2.2), we
define

(2.8)
$\displaystyle\mathbf{F(Q)}:=\mathbf{Q}\mathbf{\Sigma}(\mathbf{A},\mathbf{Q}).$

Thus one can rewrite (2.7) as

(2.9) $\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{F(Q)},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n).\end{split}\right.$

From now on, we will use (2.9) instead of (1.6) in all proofs.

###  2.2. Notations and assumptions

In this paper, we assume that $\nu$ has compact support, i.e., there exists a
constant $M>0$ such that

(2.10) $\displaystyle\mathbb{P}(\|\mathbf{x}\|_{2}\leq M)=1.$

In the following sections, we will adopt both the matrix representation and
the component-wise representation of (2.9), thus we clarify the notation here.
$\mathbf{q_{i}},\ i=1,2,...,n$ represent the column vectors of $\mathbf{Q}$ in
order, i.e.

(2.11) $\displaystyle\mathbf{Q}=[\mathbf{q_{1}},\ \mathbf{q_{2}},\ ...,\
\mathbf{q_{n}}],$

while $\tilde{\mathbf{q_{i}}},\ i=1,2,...,n$ represent the row vectors of
$\mathbf{Q}$ in order, i.e.

(2.12) $\displaystyle\mathbf{Q}^{T}=[\tilde{\mathbf{q_{1}}}^{T},\
\tilde{\mathbf{q_{2}}}^{T},\ ...,\ \tilde{\mathbf{q_{n}}}^{T}].$

For each entry, $q_{i,j},\ i,j=1,2,...,n$ represent the entries at $i$th row,
$j$th column of the matrix $\mathbf{Q}$, i.e.

(2.13) $\displaystyle\mathbf{q_{j}}=(q_{1,j},\ q_{2,j},\ ...,\ q_{n,j})^{T}.$

The canonical orthonormal basis in $\mathbb{R}^{n}$ is denoted as
$\mathbf{e_{j}},j=1,2,...,n$, which are written in column vectors, i.e.

(2.14) $\displaystyle\mathbf{I_{n}}=[\mathbf{e_{1}},\ \mathbf{e_{2}},\ ...,\
\mathbf{e_{n}}].$

Here $\mathbf{I_{n}}$ is the identity matrix of size $n$.

For $\mathbf{M},\mathbf{N}\in\mathbb{R}^{n\times n}$, $\|\mathbf{M}\|_{F}$
represents the Frobenius norm of $\mathbf{M}$ and $\langle\mathbf{M},\
\mathbf{N}\rangle_{F}$ represents the inner product in the Frobenius sense:

(2.15)
$\displaystyle\|\mathbf{M}\|=\sqrt{\mathrm{tr}(\mathbf{M}\mathbf{M}^{T})},\
\langle\mathbf{M},\
\mathbf{N}\rangle_{F}=\mathrm{tr}(\mathbf{M}^{T}\mathbf{N}).$

For $\mathbf{x}=(x_{1},x_{2},...,x_{n})\in\mathbb{R}^{n}$, $\|x\|_{2}$
represents the $\ell_{2}$ norm of $\mathbf{x}$, i.e.

(2.16) $\displaystyle\|\mathbf{x}\|_{2}=\sqrt{\sum_{j=1}^{n}|x_{j}|^{2}}$

For notation of submatrices, given $\mathbf{M}\in\mathbb{R}^{n\times n}$, row
indices $\\{a_{1},a_{2},...,a_{r}\\}$ and column indices
$\\{b_{1},b_{2},...,b_{c}\\}$, the submatrix that is formed from rows
$\\{a_{1},a_{2},...,a_{r}\\}$ and columns $\\{b_{1},b_{2},...,b_{c}\\}$ is
denoted as

(2.17) $\displaystyle\mathbf{M}[a_{1},a_{2},...,a_{r};\
b_{1},b_{2},...,b_{c}].$

The rank of a matrix $\mathbf{M}$ is denoted as $\mathrm{rank}(\mathbf{M})$.

Assume that the eigenvalues of $\mathbf{A}$ are all single, i.e. of
multiplicity one. Denote them as

(2.18) $\displaystyle\lambda_{1}>\lambda_{2}>...>\lambda_{n}>0$

in descending order. Without loss of generality, we assume that $\mathbf{A}$
is diagonal:

$\displaystyle\mathbf{A}=\mathrm{diag}\\{\lambda_{1},\lambda_{2},...,\lambda_{n}\\}.$

By default, omitted proofs of Lemmas and other important but complicated
computations are available in Section 6.

##  3\. The Landau-Lifshitz-Gilbert equation

###  3.1. The Landau-Lifshitz-Gilbert equation in $\mathbb{R}^{3}$

As a motivation, we first recall the Landau-Lifshitz-Gilbert equation [8, 5,
6] in this section. The magnetization $\mathbf{m}\in\mathbb{R}^{3}$ in a
ferromagnet is varying at each point while preserving the magnitude, which
equals to the saturation magnetization $m_{s}$. Landau and Lifshitz [8]
proposed the following differential equation of $\mathbf{m}$, which describes
the rotation of magnetization in reponse to external torque:

(3.1)
$\displaystyle\dfrac{\mathrm{d}\mathbf{m}}{\mathrm{d}t}=-\gamma\mathbf{m}\times\mathbf{h}-\lambda\mathbf{m}\times(\mathbf{m}\times\mathbf{h}).$

Here ’$\times$’ represents the cross product of two vectors in
$\mathbb{R}^{3}$, $\mathbf{h}\in\mathbb{R}^{3}$ is the effective magnetic
field applied to the magnetic moment, $\gamma$ is the electron gyromagnetic
ratio and $\lambda>0$ is a damping parameter which is related to $\gamma$ and
$m_{s}$. In 1955, Gilbert [5, 6] modified (3.1) by introducing parameters
characterizing the material property. The modified version was of the same
form as (3.1), but with $\gamma$ and $\lambda$ of different physical meanings.
Therefore, we refer to (3.1) as the Landau-Lifshitz-Gilbert (LLG) equation.

The effective magnetic field $\mathbf{h}$ is the negative derivative of a
magnetic energy density function w.r.t the magnetization $\mathbf{m}$, i.e.

(3.2) $\displaystyle\mathbf{h}=-\dfrac{\partial
F(\mathbf{m})}{\partial\mathbf{m}}.$

Versatile choices of the magnetic energy were considered under different
physical contexts. For instance, one can choose the exchange energy
$F(\mathbf{m})=\dfrac{1}{2}\int_{\Omega}|\nabla\mathbf{m}|^{2}\mathrm{d}\mathbf{x}$
to arrange the molecular magnetic field in order.

Now suppose that the energy density function $F(\mathbf{m})$ is given. A
direct calculation yields

(3.3) $\displaystyle\dfrac{\mathrm{d}F(\mathbf{m}(t))}{\mathrm{d}t}$
$\displaystyle=-\dfrac{\mathrm{d}\mathbf{m}}{\mathrm{d}t}\cdot\mathbf{h}$
$\displaystyle=\gamma(\mathbf{m}\times\mathbf{h})\cdot\mathbf{h}+\lambda\mathbf{h}\cdot(\mathbf{m}\times(\mathbf{m}\times\mathbf{h}))$
$\displaystyle=-\lambda|\mathbf{m}\times\mathbf{h}|^{2}\leq 0.$

Thus $F(\mathbf{m})$ is dissipated. As the first term is perpendicular to the
gradient direction, i.e. $\mathbf{m}\times\mathbf{h}\perp\mathbf{h}$, the
dissipation of the free energy is totally contributed by the second term, i.e.
$-\lambda\mathbf{m}\times(\mathbf{m}\times\mathbf{h})$. Therefore, we refer to
this term as the dissipative term.

Another observation is that

(3.4)
$\displaystyle\dfrac{\mathrm{d}\|\mathbf{m}\|^{2}}{\mathrm{d}t}=2\mathbf{m}\cdot(-\gamma\mathbf{m}\times\mathbf{h}-\lambda\mathbf{m}\times(\mathbf{m}\times\mathbf{h}))=0.$

Therefore, $\|\mathbf{m}\|$ is also preserved. So if
$\mathbf{m}(0)\in\mathbb{S}^{2}$, then $\mathbf{m}(t)\in\mathbb{S}^{2}$ for
all $t>0$.

If $\lambda=0$ in (3.1), then $F(\mathbf{m}(t))$ is conserved. If
$F(\mathbf{m})$ is selected as the Kinetic energy and $\mathbf{m}$ represents
the body angular velocity of a rigid body under free rotation, then the LLG
equation reduces to the Euler equation of a rigid body. Therefore, the first
term should be regarded as the Hamiltonian part of the LLG equation, which
conserves the free energy $F(\mathbf{m})$.

###  3.2. Generalization of the LLG equation: on Riemmanian manifolds

In the previous subsection, we observed that the critical structure of the LLG
equation is that: first, it possesses both a Hamiltonian part and a
dissipative part, and they are perpendicular to each other; second, the
Hamiltonian for the Hamiltonian part is exactly the free energy minimized by
the dissipative part, i.e. $F(\mathbf{m})$. Moreover, $\|\mathbf{m}\|$ is
preserved, so solutions of (3.1) is invariant on $\mathbb{S}^{2}$.

Therefore, we desire the LLG equation on a Riemannian manifold
$(\mathcal{M},g)$ should preserve the above properties:

  * •

It possesses the Hamiltonian-dissipation structure: the Hamiltonian term $H$
is perpendicular to the dissipative term $D$ on $(\mathcal{M},g)$; $D$ is the
variation of a free energy $\mathcal{E}$; $H$ preserves $\mathcal{E}$.

  * •

The solution $\mathbf{m}(t)\in\mathcal{M}$. For instance, the solution of
(3.1) stays on $\mathbb{S}^{2}$ and the solution of (2.9) satys on $O(n)$.

Suppose that $(\mathcal{M},g)$ is a Riemannian manifold (see Section 6 for
definition). Let $\mathcal{E}(\mathbf{m}):\ \mathcal{M}\to\mathbb{R}$ be a
smooth function on $\mathcal{M}$ and $\nabla_{g}\mathcal{E}$ be the gradient
of $\mathcal{E}$ w.r.t. the Riemannian metric $g$ (see Section 6 for
definition). Then by definition of the gradient, for any tangent vector field
$\dot{\mathbf{x}}$, we have

(3.5) $\displaystyle\langle\mathcal{E}^{\prime}(\mathbf{m}),\
\dot{\mathbf{x}}\rangle=\langle\nabla_{g}\mathcal{E},\
\dot{\mathbf{x}}\rangle_{g}.$

Here $\mathcal{E}^{\prime}(\mathbf{m})$ is the derivative of $\mathcal{E}$,
which is a cotangent vector field. The bracket $\langle\cdot,\cdot\rangle$ is
understood as the action of a cotangent vector field (in this case,
$\mathcal{E}^{\prime}$) on a tangent vector field.

The gradient flow of the free energy $\mathcal{E}$ on $(\mathcal{M},g)$ is the
following differential equation:

(3.6)
$\displaystyle\dfrac{\mathrm{d}\mathbf{m}(t)}{\mathrm{d}t}=-\nabla_{g}\mathcal{E}(\mathbf{m}).$

The term $-\nabla_{g}\mathcal{E}(\mathbf{m})$ will serve as the dissipative
term in the LLG equation on $\mathcal{M}$.

The Hamiltonian part is defined as following: for any time dependent tangent
vector field $\dot{\mathbf{m}}_{1}(\mathbf{m},t):\
\mathcal{M}\times[0,\infty)\to T\mathcal{M}$, define

(3.7)
$\displaystyle\widehat{\dot{\mathbf{m}}}_{1}:=\|\nabla_{g}\mathcal{E}(\mathbf{m})\|_{g}^{2}\cdot\dot{\mathbf{m}}_{1}-\langle\nabla_{g}\mathcal{E}(\mathbf{m}),\
\dot{\mathbf{m}}_{1}\rangle_{g}\cdot\nabla_{g}\mathcal{E}(\mathbf{m}).$

Then at each fixed time $t$, $\widehat{\dot{\mathbf{m}}}_{1}(\cdot,t)$ is
still a tangent vector field since it is the sum of two tangent vector fields.
Moreover, direct calculation yields

(3.8) $\displaystyle\langle\nabla_{g}\mathcal{E}(\mathbf{m}),\
\widehat{\dot{\mathbf{m}}}_{1}\rangle_{g}$
$\displaystyle=\|\nabla_{g}\mathcal{E}(\mathbf{m})\|_{g}^{2}\langle\nabla_{g}\mathcal{E}(\mathbf{m}),\
{\dot{\mathbf{m}}}_{1}\rangle_{g}-\|\nabla_{g}\mathcal{E}(\mathbf{m})\|_{g}^{2}\langle\nabla_{g}\mathcal{E}(\mathbf{m}),\
{\dot{\mathbf{m}}}_{1}\rangle_{g}$ $\displaystyle=0.$

Thus $\widehat{\dot{\mathbf{m}}}_{1}$ is perpendicular to the gradient
direction and preserves the free energy.

Upon (3.6) and (3.7), we are ready to develop the LLG equation on
$(\mathcal{M},g)$. For any time dependent tangent vector field
$\dot{\mathbf{m}}_{1}$, the following differential equation serves as an
analog of (3.1), i.e. the LLG equation on the Riemannian manifold
$(\mathcal{M},g)$:

(3.9)
$\displaystyle\dfrac{\mathrm{d}\mathbf{m}(t)}{\mathrm{d}t}=\widehat{\dot{\mathbf{m}}}_{1}(\mathbf{m}(t),t)-\nabla_{g}\mathcal{E}(\mathbf{m}(t)),$

where $\widehat{\dot{\mathbf{m}}}_{1}$ is defined in (3.7). Then we have the
following proposition for (3.9):

######  Proposition 3.1.

(the LLG equation on a Riemannian manifold) Suppose that $(\mathcal{M},g)$ is
a Riemannian manifold. Let $\mathcal{E}(\mathbf{m}):\
\mathcal{M}\to\mathbb{R}$ be a smooth function on $\mathcal{M}$ and
$\dot{\mathbf{m}}_{1}$ be a time dependent tangent vector field
$\dot{\mathbf{m}}_{1}(\mathbf{m},t):\ \mathcal{M}\times[0,\infty)\to
T\mathcal{M}$. Consider equation (3.9) with initial value
$\mathbf{m}(0)=\mathbf{m}_{0}\in\mathcal{M}$, i.e.

$\displaystyle\dfrac{\mathrm{d}\mathbf{m}(t)}{\mathrm{d}t}=\widehat{\dot{\mathbf{m}}}_{1}-\nabla_{g}\mathcal{E}(\mathbf{m}),\
\mathbf{m}(0)=\mathbf{m}_{0},$

where $\hat{\dot{\mathbf{m}}}_{1}$ is defined in (3.7). Then:

  1. (i)

($\mathcal{M}$ invariance) for any $t\geq 0$, $\mathbf{m}(t)\in\mathcal{M}$;

  2. (ii)

(energy dissipation) for any $t\geq 0$,

(3.10)
$\displaystyle\dfrac{\mathrm{d}\mathcal{E}(\mathbf{m}(t))}{\mathrm{d}t}=-\|\nabla_{g}\mathcal{E}(\mathbf{m}(t))\|_{g}^{2}\leq
0.$

###### Proof.

Because both $\nabla_{g}\mathcal{E}$ and $\dot{\mathbf{m}_{1}}$ are tangent
vector fields, so $(i)$ is straightforward. For $(ii)$, from definition (6.2)
(the definition of gradients on Riemannian manifolds), direct calculation
implies

$\displaystyle\dfrac{\mathrm{d}\mathcal{E}(\mathbf{m}(t))}{\mathrm{d}t}$
$\displaystyle=\langle\mathcal{E}^{\prime}(\mathbf{m}(t)),\
\widehat{\dot{\mathbf{m}}}_{1}-\nabla_{g}\mathcal{E}(\mathbf{m}(t))\rangle$
$\displaystyle=\langle\nabla_{g}\mathcal{E}(\mathbf{m}(t)),\widehat{\dot{\mathbf{m}}}_{1}\rangle_{g}-\langle\nabla_{g}\mathcal{E}(\mathbf{m}(t)),\nabla_{g}\mathcal{E}(\mathbf{m}(t))\rangle_{g}.$

By (3.8), we konw that
$\langle\nabla_{g}\mathcal{E}(\mathbf{m}(t)),\widehat{\dot{\mathbf{m}}}_{1}\rangle_{g}=0$.
Therefore,

$\displaystyle\dfrac{\mathrm{d}\mathcal{E}(\mathbf{m}(t))}{\mathrm{d}t}=-\langle\nabla_{g}\mathcal{E}(\mathbf{m}(t)),\nabla_{g}\mathcal{E}(\mathbf{m}(t))\rangle_{g}=-\|\nabla_{g}\mathcal{E}(\mathbf{m}(t))\|_{g}^{2}\leq
0.$

This proves (3.10). ∎

In summary, we generalized the LLG equation to any Riemmanian manifold. Based
on it, we are going to embed the problem of online PCA into the macroscopic
framework of the LLG equation on Riemannian manifolds.

##  4\. Online PCA: an LLG equation on the Stiefel manifold

In this section, we formulate the corresponding ODE of online PCA, i.e.
equation (2.9) as an LLG equation on the Stiefel manifold. Lemma 1.1 ensures
that solutions of (2.9) will stay on the Stiefel manifold if the initial value
$\mathbf{Q}_{0}\in O(n)$, so we only need to select a free energy
$\mathcal{E}(\mathbf{Q}):O(n)\to\mathbb{R}$ and assign an appropriate
Riemannian metric $g$ to the Stiefel manifold.

###  4.1. The weighted Rayleigh quotient and the Wielandt-Hoffman inequality

A natural choice of the free energy is the Rayleigh quotient [7]:

(4.1) $\displaystyle
R_{\mathbf{A}}(\mathbf{Q})=\mathrm{tr}(\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}),$

since maximizing the Rayleigh Quotient is equivalent to exploring the
principle eigenspaces. However, if one desires to not only detect the
principle eigenspaces, but also to distinguish each principle eigenvector from
each other, certain modification is necessary on the Rayleigh quotient to
achieve this aim. In fact, as long as $\mathbf{Q}\in O(n)$, the Rayleigh
quotient is constant and equals to $\mathrm{tr}(\mathbf{A})$:

$\displaystyle
R_{\mathbf{A}}(\mathbf{Q})=\mathrm{tr}(\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})=\mathrm{tr}(\mathbf{Q}\mathbf{Q}^{T}\mathbf{A})=\mathrm{tr}(\mathbf{A}).$

Therefore, the Rayleigh quotient can not be selected as the free energy
directly.

To distinguish different eigenvectors, one way is to assign different weights
to different components [7]. Let $\mathbf{N}\in\mathbb{R}^{n\times n}$ be a
symmetric matrix. Consider

(4.2)
$\displaystyle\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})=\mathrm{tr}(\mathbf{N}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}).$

As long as $\mathbf{N}$ is given, this is equivalent to selecting the
diagonalization of $\mathbf{N}$. Suppose that $\mathbf{P}\in O(n)$
diagonalizes $\mathbf{N}$, i.e.,
$\mathbf{N}=\mathbf{P}\mathbf{D}\mathbf{P}^{T}$ where $\mathbf{D}$ is
diagonal. Then

$\displaystyle\max\limits_{\mathbf{Q}\in
O(n)}\mathrm{tr(\mathbf{N}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})}$
$\displaystyle=\max\limits_{\mathbf{Q}\in
O(n)}\mathrm{tr(\mathbf{P}\mathbf{D}\mathbf{P}^{T}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})}=\max\limits_{\mathbf{Q}\in
O(n)}\mathrm{tr(\mathbf{D}\mathbf{P}^{T}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}\mathbf{P})}$
$\displaystyle=\max\limits_{\mathbf{QP}\in
O(n)}\mathrm{tr(\mathbf{D}\mathbf{P}^{T}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}\mathbf{P})}=\max\limits_{\mathbf{Q}\in
O(n)}\mathrm{tr(\mathbf{D}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q})}.$

So without loss of generality, we assume that $\mathbf{N}$ is diagonal.

To justify that maximizing (4.2) does provide the eigenspace decomposition, we
need the Wielandt-Hoffman inequality:

######  Lemma 4.1.

(the Wielandt-Hoffman inequality) Suppose that
$\mathbf{M},\mathbf{N}\in\mathbb{R}^{n\times n}$ are symmetric matrices. Let
$\lambda_{1,i},i=1,2,...,n$ and $\lambda_{2,i},i=1,2,...,n$ be eigenvalues of
$\mathbf{M}$ and $\mathbf{N}$ respectively in descending order, i.e.

$\displaystyle\lambda_{1,1}$
$\displaystyle\geq\lambda_{1,2}\geq...\geq\lambda_{1,2},$
$\displaystyle\lambda_{2,1}$
$\displaystyle\geq\lambda_{2,2}\geq...\geq\lambda_{2,n}.$

Then

(4.3)
$\displaystyle\|\mathbf{M}-\mathbf{N}\|_{F}^{2}\geq\sum_{i=1}^{n}(\lambda_{1,i}-\lambda_{2,i})^{2}.$

The equality holds if and only if $\mathbf{M}$ and $\mathbf{N}$ can be
simultaneously diagonalized with diagonal entries aligned in descending order,
i.e., there exists $\mathbf{P}\in O(n)$ such that

(4.4)
$\displaystyle\mathbf{P}^{T}\mathbf{M}\mathbf{P}=\mathrm{diag}\\{\lambda_{1,1},\lambda_{1,2},...\lambda_{1,n}\\},\
\mathbf{P}^{T}\mathbf{N}\mathbf{P}=\mathrm{diag}\\{\lambda_{2,1},\lambda_{2,2},...\lambda_{2,n}\\}$

See Section 6 for the proof of the above inequality.

By Lemma 4.1, we can prove that maximizing the weighted Rayleigh quotient is
exactly decomposing $\mathbb{R}^{n}$ into eigenspaces of $\mathbf{A}$ in
order:

######  Proposition 4.1.

Suppose that $\mathbf{N}=\mathrm{diag}\\{\mu_{1},\mu_{2},...,\mu_{n}\\}$ with
$\mu_{1}>\mu_{2}>\mu_{3}>...>\mu_{n}$. Consider the weighted Rayleigh quotient
defined in (4.2). Then

(4.5) $\displaystyle\max\limits_{\mathbf{Q}\in
O(n)}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})=\sum_{i=1}^{n}\mu_{i}\lambda_{i}.$

Here $\lambda_{i}$ are eigenvalues of $\mathbf{A}$ defined in (2.18). The
maximum is attained if and only if

(4.6)
$\displaystyle\mathbf{Q}^{T}\mathbf{AQ}=\mathrm{diag}\\{\lambda_{1},\lambda_{2},...,\lambda_{n}\\}.$

###### Proof.

Using the fact

$\displaystyle\|\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}-\mathbf{N}\|^{2}_{F}=\|\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}\|^{2}+\|\mathbf{N}\|^{2}-2\mathrm{tr}(\mathbf{N}\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}),$

we derive

$\displaystyle\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$
$\displaystyle=\mathrm{tr}(\mathbf{N}\mathbf{Q}^{T}\mathbf{AQ})$
$\displaystyle=\dfrac{\|\mathbf{Q}^{T}\mathbf{AQ}\|_{F}^{2}+\|\mathbf{N}\|_{F}^{2}-\|\mathbf{Q}^{T}\mathbf{AQ}-\mathbf{N}\|_{F}^{2}}{2}$
$\displaystyle=\dfrac{\|\mathbf{A}\|_{F}^{2}+\|\mathbf{N}\|_{F}^{2}-\|\mathbf{Q}^{T}\mathbf{AQ}-\mathbf{N}\|_{F}^{2}}{2}.$

By Lemma 4.1, we have

(4.7)
$\displaystyle\|\mathbf{Q}^{T}\mathbf{A}\mathbf{Q}-\mathbf{N}\|_{F}^{2}\geq\sum_{i=1}^{n}(\lambda_{i}-\mu_{i})^{2}.$

Thus

(4.8)
$\displaystyle\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})\leq\dfrac{1}{2}\left(\sum_{i=1}^{n}\lambda_{i}^{2}+\mu_{i}^{2}-(\lambda_{i}-\mu_{i})^{2}\right)=\sum_{i=1}^{n}\mu_{i}\lambda_{i}.$

The equality holds if $\mathbf{Q}$ diagonalizes $\mathbf{A}$ with descending
order diagonal entries, i.e. (4.6) holds. ∎

Lemma 4.1 incorporates the meaning of introducing the weight $\mathbf{N}$: by
the rearrangement inequality, it distinguishs different eigenvectors by
assigning different weights to them and forces the column of $\mathbf{Q}$ to
align in a descending order of corresponding eigenvalues. This mechanism
facilitates the eigenspace decomposition.

Indeed, $-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ is a Lyapunov
function for (2.9):

######  Lemma 4.2.

Suppose that $\mathbf{Q}(t),\ t\geq 0$ is the solution to (2.9) for some
initial value $\mathbf{Q_{0}}\in O(n)$. Let
$\mathbf{N}=\mathrm{diag}\\{\mu_{1},\mu_{2},...,\mu_{n}\\}\in\mathbb{R}^{n\times
n}$ satisfy

(4.9) $\displaystyle\mu_{1}>\mu_{2}>...>\mu_{n}.$

Define

(4.10) $\displaystyle E:=\\{\mathrm{diag}\\{\varepsilon_{1},\
\varepsilon_{2},\ ...,\ \varepsilon_{n}\\}\mathbf{P}:\ \varepsilon_{i}=\pm 1,\
i=1,2,...,n,\ \mathbf{P}\ \mathrm{is}\ \mathrm{a}\ \mathrm{permutation}\
\mathrm{matrix.}\\}.$

Then

  1. (i)

(Lyapunov function) Consider the weighted Rayleigh quotient
$\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ defined in (4.2). Then

(4.11)
$\displaystyle\dfrac{\mathrm{d}(-\mathcal{E}(\mathbf{Q}(t);\mathbf{A},\mathbf{N}))}{\mathrm{d}t}=-2\sum_{j=1}^{n}\sum_{k=1}^{j-1}(\mu_{k}-\mu_{j})(\mathbf{q_{k}}(t)\cdot\mathbf{Aq_{j}}(t))^{2}\leq
0;$

  2. (ii)

(equilibrium points) the ’=’ in (4.11) holds if and only if $\mathbf{Q}\in E$;

  3. (iii)

(stable points) $\mathbf{F}(\mathbf{Q})=0$ if and only if $\mathbf{Q}\in E$;

  4. (iv)

(finite set) $|E|=2^{n}\cdot n!$.

See Section 6 for the proof of this lemma.

Therefore, we will choose the weighted Rayleigh quotient as the free energy of
the LLG equation.

###  4.2. Construction of the Riemmanian metric

Another important feature for the LLG equation or a gradient flow structure is
the Riemmanian metric, which defines the geodesic and the distance between two
points on the manifold. In this section, we will construct a Riemannian metric
on the Stiefel manifold (after exclusion of a finite number of equilibria) so
that (2.9) can be interpreted as a LLG equation of the weighted Rayleigh
quotient (see Definition (4.2)), which is indeed a Lyapunov function of (2.9).
From this Lyapunov function, we will utilize the method in [1] to construct
the Riemannian manifold. We will first introduce the general method in [1] and
then apply it to our context.

######  Definition 4.1.

Let $\mathbf{F}$ be a continuous tangent vector field on a manifold $M$. A
continuously differentiable function $\mathcal{E}:M\to\mathbb{R}$ is a strict
Lyapunov function for

(4.12) $\displaystyle\dot{\mathbf{Q}}=\mathbf{F(Q)},$

if

(4.13)
$\displaystyle\langle\mathcal{E}^{\prime}(\mathbf{Q}),\mathbf{F(Q)}\rangle<0\
\mathrm{whenever}\ \mathbf{Q}\in M\ \mathrm{and}\
\mathbf{F(Q)}\neq\mathbf{0}.$

According to Theorem 1 in [1], we can now construct the Riemannian metirc
$\tilde{g}$.

######  Proposition 4.2.

(the gradient flow on M) Let $M_{0}$ be a manifold, $\mathbf{F}$ a continuous
tangent vector field on $M_{0}$ and $\mathcal{E}:M_{0}\to\mathbb{R}$ be a
continuously differentiable, strict Lyapunov function (see Definition 4.1) for
(4.12). Then there exists a Riemannian metric $\tilde{g}$ on the open set

(4.14) $\displaystyle M:=\\{\mathbf{Q}\in M_{0}:\ \mathbf{F}(\mathbf{Q})\neq
0\\}$

such that

(4.15) $\displaystyle\nabla_{\tilde{g}}\mathcal{E}=-\mathbf{F}.$

In particular, (4.12) is a gradient system on the Riemannian manifold
$(M,\tilde{g})$.

###### Proof.

Define the kernel of a cotangent vector field $\mathcal{F}$ at each tangent
space $T_{\mathbf{Q}}M$ is defined as

(4.16) $\displaystyle\mathrm{Ker}(\mathcal{F}(\mathbf{Q})):=\\{\Omega\in
T_{\mathbf{Q}}M:\ \langle\mathcal{F}(\mathbf{Q}),\Omega\rangle_{F}=0\\}.$

Equivalently, each cotangent vector field $\mathcal{F}$ is viewed as a bounded
linear functional on each tangent space, and the kernel of it is the kernel of
the linear functional.

For every $\mathbf{Q}\in M$, we have $\mathbf{F(Q)}\neq 0$ and
$\langle\mathcal{E}^{\prime}(\mathbf{Q}),\ \mathbf{F(Q)}\rangle_{F}<0$. Thus

(4.17)
$\displaystyle\mathbf{F(Q)}\notin\mathrm{Ker}\mathcal{E}^{\prime}(\mathbf{Q}),\
\mathrm{Ker}\mathcal{E}^{\prime}(\mathbf{Q})\neq T_{\mathbf{Q}}M.$

Because $\mathcal{E}^{\prime}$ is a non-trivial cotangent vector field, so
$\mathrm{dim}(\mathrm{Ran}(\mathcal{E}^{\prime}))=1$, which is the dimension
of $\mathbb{R}$ (as a linear space on $\mathbb{R}$), thus

$\displaystyle\dim(\mathrm{Ker}(\mathcal{E}^{\prime}(\mathbf{Q})))$
$\displaystyle=\dim(T_{\mathbf{Q}}M)-\dim(\mathrm{Ran}(\mathcal{E}^{\prime}(\mathbf{Q})))$
$\displaystyle=\dim(T_{\mathbf{Q}}M)-1.$

Denote $\langle\mathbf{F}\rangle$ be the linear space generated by
$\mathbf{F}$, which is a subspace of $T_{\mathbf{Q}}M$, thus

$\displaystyle\dim(\mathrm{Ker}(\mathcal{E}^{\prime}(\mathbf{Q})))+\dim(\langle\mathbf{F(\mathbf{Q})}\rangle)=\dim(T_{\mathbf{Q}}M).$

Thus by (4.17), we know that the tangent bundle $TM$ is the direct sum of the
bundle $\mathrm{Ker}\mathcal{E}^{\prime}$ and the bundle gerenated by the
vector field $\mathbf{F(Q)}$:

(4.18) $\displaystyle
TM=\mathrm{Ker}(\mathcal{E}^{\prime})\oplus\langle\mathbf{F}\rangle.$

For every tangent vector field $\mathbf{X}$ on $M$, we define the
corresponding decomposition, which are well-defined and continuous:

(4.19) $\displaystyle\mathbf{X}=\mathbf{X_{0}}+\mathbf{X_{1}},\
\mathbf{X_{0}}:=\mathbf{X}-\dfrac{\langle\mathcal{E}^{\prime},\
\mathbf{X}\rangle_{F}}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\mathbf{F}\in\mathrm{Ker}\mathcal{E}^{\prime},\
\mathbf{X_{1}}:=\dfrac{\langle\mathcal{E}^{\prime},\
\mathbf{X}\rangle_{F}}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\mathbf{F}\in\langle\mathbf{F}\rangle.$

Suppose that $g$ is an arbitrary Riemannian metric on (4.23), for every
tangent vector fields $\mathbf{X,Y}$ on $M$, define

(4.20) $\displaystyle\langle\mathbf{X},\mathbf{Y}\rangle_{\tilde{g}}$
$\displaystyle:=\langle\mathbf{X_{0}},\mathbf{Y_{0}}\rangle_{g}-\dfrac{1}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\langle\mathcal{E}^{\prime},\
\mathbf{X}\rangle_{F}\langle\mathcal{E}^{\prime},\ \mathbf{Y}\rangle_{F}$
$\displaystyle=\langle\mathbf{X_{0}},\mathbf{Y_{0}}\rangle_{g}-\dfrac{1}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\langle\mathcal{E}^{\prime},\
\mathbf{X_{1}}\rangle_{F}\langle\mathcal{E}^{\prime},\
\mathbf{Y_{1}}\rangle_{F}.$

Because $g,\mathcal{E}^{\prime}$ and $\mathbf{F}$ are contiuous, so is
$\tilde{g}$. Moreover, for each $\mathbf{Q}\in M$ and any $\mathbf{X}\in
T_{\mathbf{Q}}M$, we have

(4.21) $\displaystyle\langle\mathbf{X},\mathbf{X}\rangle_{\tilde{g}}$
$\displaystyle=\langle\mathbf{X_{0}},\mathbf{X_{0}}\rangle_{g}-\dfrac{\langle\mathcal{E}^{\prime},\
\mathbf{X}\rangle_{F}^{2}}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\geq 0.$

Here we use the fact that
$\langle\mathcal{E}^{\prime},\mathbf{F}\rangle_{F}\leq 0$. The quality holds
if and only if $\mathbf{X}=0$. So $\tilde{g}$ is also positive definite. Thus
$\tilde{g}$ is a Riemannian metric on $M$. Notice that $\tilde{g}$ is
independent with the manifold structure of $M$, so $(M,\tilde{g})$ is a
Riemannian manifold.

Finally, we check the compatibility condition. For every tangent vector field
on $M$, we have

(4.22) $\displaystyle\langle\mathbf{F},\
\mathbf{X}\rangle_{\tilde{g}}=0-\dfrac{1}{\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}}\langle\mathcal{E}^{\prime},\
\mathbf{F}\rangle_{F}\langle\mathcal{E}^{\prime},\
\mathbf{X}\rangle_{F}=-\langle\mathcal{E}^{\prime},\ \mathbf{X}\rangle_{F}.$

Thus $\nabla_{\tilde{g}}\mathcal{E}=-\mathbf{F}$. ∎

######  Remark 4.1.

In fact, the Riemannian metric $\tilde{g}$ is not defined at points where
$\langle\mathcal{E}^{\prime}(\mathbf{Q}),\mathbf{F}(\mathbf{Q})\rangle_{F}=0$.
This discontinuous point can be removable or irremovable depending on
selection of the original metric $g$, see [1]. Meanwhile, the Riemannian
metric is not degenerate around equilibrium points due to
$\langle\mathcal{E}^{\prime},\ \mathbf{F}\rangle_{F}$ on the denominator.

Now we focus on the context of (2.9) and the Stiefel manifold. Remember that
we select $-\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ in (4.2) as the
free energy. Now consider

(4.23) $\displaystyle M:=O(n)\setminus E.$

Here $E$ is defined in (4.10).

By Lemma 4.2, we have: first, because $|E|<\infty$, $M$ is an open subset of
the Stiefel manifold. Second, we know that for all $\mathbf{Q}\in M$, we have
$\mathbf{F}(\mathbf{Q})\neq 0$ and
$\langle\mathcal{E}^{\prime}(\mathbf{Q};\mathbf{A},\mathbf{N}),\mathbf{F(Q)}\rangle<0$.
Thus $\mathcal{E}$ is a strict Lyapunov function on $M$. These properties
allow us to apply Proposition 4.2 on $M$ and
$\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ in (4.2):

######  Theorem 4.1.

Consider $\mathbf{F}$ in (2.8), $M$ in (4.23) and
$\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ in (4.2). Then there exists a
Riemannian metric $\tilde{g}$ on $M$ such that

(4.24) $\displaystyle\nabla_{\tilde{g}}\mathcal{E}=-\mathbf{F}.$

###### Proof.

Let $g$ be the Euclidean metric, i.e., for each $\mathbf{Q}\in M$ and
$\mathbf{X,Y}\in T_{\mathbf{Q}}M$, let

$\displaystyle\langle\mathbf{X},\mathbf{Y}\rangle_{g}=\langle\mathbf{X},\mathbf{Y}\rangle_{F}.$

Notice that $\mathcal{E}$ in (4.2) is a smooth strict Lyapunov function of
(2.9), and $\mathbf{F}$ defines a tangent vector field on $M$, so by
Proposition (4.2), $\tilde{g}$ defined in (4.20) is a Riemannian metric on $M$
and (4.24) holds. ∎

######  Remark 4.2.

The choice of $g$ in the proof above is arbitrary as long as it is a
Riemmanian metric on the Stiefel manifold $O(n)$.

###  4.3. The LLG equation

Now after selecting the free energy and assigning the Riemannian metric, we
can formulate the ODE of the SGA method, i.e. (2.9) as a LLG equation on the
Stiefel manifold with a Riemmanian metric.

Suppose $\mathbf{Q}\in O(n)$. Then the tangent space at $\mathbf{Q}$, i.e.
$T_{\mathbf{Q}}O(n)$ is

(4.25) $\displaystyle T_{\mathbf{Q}}O(n)=\\{\mathbf{Q}\bm{\Omega}:\
\bm{\Omega}+\ \bm{\Omega}^{T}=0\\}.$

Then any tangent vector field on $O(n)$ can be written in the form
$\mathbf{QS(Q)}$, where $\mathbf{S(Q)}\in\mathbb{R}^{n\times n}$ satisfies
$\mathbf{S(Q)}+\mathbf{S(Q)}^{T}=0$. Then by Proposition 3.1, we know that

(4.26)
$\displaystyle\left\\{\begin{split}\dfrac{\mathrm{d}\mathbf{Q}(t)}{\mathrm{d}t}&=\widehat{\mathbf{Q}\mathbf{S(Q)}}+\nabla_{\tilde{g}}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N}),\\\
\widehat{\mathbf{Q}\mathbf{S(Q)}}&=\|\nabla_{\tilde{g}}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})\|^{2}_{\tilde{g}}\cdot\mathbf{Q}\mathbf{S(Q)}-\langle\nabla_{\tilde{g}}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N}),\mathbf{Q}\mathbf{S(Q)}\rangle_{\tilde{g}}\cdot\nabla_{\tilde{g}}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N}),\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n)\end{split}\right.$

is an LLG equation on the Riemannian manifold $(O(n)\setminus E,\tilde{g})$.
Here the hat operator is defined in (3.7), $E$ is the stable points of
$\mathbf{F}(\mathbf{Q})$ defined in (4.10) and $\tilde{g}$ is the Riemannian
metric defined in (4.20) in Proposition 4.2, by taking $g$ as the Euclidean
metric. If we substitute the explicit expression of $\tilde{g}$ in (4.20),
(4.26) then reads as

(4.27)
$\displaystyle\left\\{\begin{split}\dfrac{\mathrm{d}\mathbf{Q}(t)}{\mathrm{d}t}&=\widehat{\mathbf{Q}\mathbf{S(Q)}}+\mathbf{F(Q)},\\\
\widehat{\mathbf{Q}\mathbf{S(Q)}}&=2\mathrm{tr}(\mathbf{NQ}^{T}\mathbf{AF(Q)})\mathbf{QS(Q)}-2\mathrm{tr}(\mathbf{NQ}^{T}\mathbf{AQS(Q)})\mathbf{F(Q)},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n).\end{split}\right.$

Again, $\mathbf{S(Q)}\in\mathbb{R}^{n\times n}$ is arbitrary smooth (w.r.t.
$\mathbf{Q}$) skew-symmetric matrices field. The convergence behavior of this
system will be thoroughly studied in the following sections.

Equation (2.9) is a special case of (4.26) by taking $\mathbf{S}=0$, which
reduces to the gradient flow system. Although the Riemmanian metric
$\tilde{g}$ is constructed on $O(n)\setminus E$ (see Definition (4.10)), the
LLG equation is well defined on the whole Stiefel manifold: if
$\mathbf{Q_{0}}\in E$, then the flow stays invariant, i.e.
$\mathbf{Q}(t)=\mathbf{Q_{0}}$.

In [7], the Oja-Brockett flow was also formulated as the gradient flow of the
weighted Rayleigh quotient on $(O(n),g_{e})$, where $g_{e}$ is the Euclidean
metric. The gradient of $\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$ on
$(O(n),g_{e})$ is given by

(4.28)
$\displaystyle\nabla_{g_{e}}\mathcal{E}(\mathbf{Q};\mathbf{A},\mathbf{N})$
$\displaystyle=\dfrac{1}{2}(\mathcal{E}^{\prime}(\mathbf{Q};\mathbf{A},\mathbf{N})-\mathbf{Q}(\mathcal{E}^{\prime}(\mathbf{Q};\mathbf{A},\mathbf{N}))^{T}\mathbf{Q})$
$\displaystyle=\mathbf{AQN}-\mathbf{QNQ}^{T}\mathbf{AQ}.$

See Section 6 for calculation of gradients in $(O(n),g_{e})$. If we replace
$\tilde{g}$ by the Euclidean metric in (4.26), then we have

(4.29)
$\displaystyle\left\\{\begin{split}\dfrac{\mathrm{d}\mathbf{Q}(t)}{\mathrm{d}t}&=\widehat{\mathbf{Q}\mathbf{S(Q)}}+\mathbf{AQN-
QN}\mathbf{Q}^{T}\mathbf{AQ},\\\
\widehat{\mathbf{Q}\mathbf{S(Q)}}&=2\mathrm{tr}(\mathbf{Q}^{T}\mathbf{AQN}^{2}-(\mathbf{Q}^{T}\mathbf{AQN})^{2})\mathbf{Q}\mathbf{S(Q)}\\\
&-\mathrm{tr}(\mathbf{NQ}^{T}\mathbf{AQS(Q)}-\mathbf{Q}^{T}\mathbf{AQNS(Q)})(\mathbf{AQN-
QN}\mathbf{Q}^{T}\mathbf{AQ}),\\\ \mathbf{Q}(0)&=\mathbf{Q_{0}}\in
O(n).\end{split}\right.$

If $\mathbf{S}=\mathbf{0}$, then we recover the Oja-Brockett flow [4].

##  5\. Convergence analysis of the LLG equation on the Stiefel manifold

###  5.1. Convergence to equilibria

Notice that the number of stable points of (2.9) is finite, we can prove the
following convergence theorem:

######  Theorem 5.1.

(convergence) For any initial value $\mathbf{Q_{0}}\in O(n)$, suppose that
$\mathbf{Q}(t),\ t\geq 0$ is the solution of (4.26) or (4.29). Then there
exists $\mathbf{Q}^{*}\in E$ such that

(5.1) $\displaystyle\lim\limits_{t\to\infty}\mathbf{Q}(t)=\mathbf{Q}^{*}.$

Here $E$ is defined in (4.23).

###### Proof.

If $\mathbf{Q_{0}}\in E$, then $\mathbf{Q}(t)=\mathbf{Q_{0}}$ for any $t\geq
0$, hence $\mathbf{Q}^{*}=\mathbf{Q_{0}}$. If $\mathbf{Q_{0}}\notin E$, by
Proposition 3.1 and Lemma (4.2), we have

$\displaystyle-\dfrac{\mathrm{d}\mathcal{E}(\mathbf{Q}(t);\mathbf{A},\mathbf{N})}{\mathrm{d}t}=-\|\nabla_{g}\mathcal{E}(\mathbf{Q}(t);\mathbf{A},\mathbf{N})\|_{g}^{2}\leq
0,$

where $g=\tilde{g}$ for (4.26) and $g=g_{e}$ for (4.29).

By Lasalle’s invariance principle [9], the $\omega$-limit set of (4.26) or
(4.29) is a subset of zeros of
$\nabla_{g}\mathcal{E}(\mathbf{Q}(t);\mathbf{A},\mathbf{N})$. For (4.26),
$\nabla_{g}\mathcal{E}(\mathbf{Q}(t);\mathbf{A},\mathbf{N})=\mathbf{F(Q)}$, so
by Lemma 4.2, zeros of the gradient is exactly the set $E$; for (4.29), by
(4.24), any point in the $\omega-$limit set should solve

$\displaystyle\mathbf{AQN}-\mathbf{QN}\mathbf{Q}^{T}\mathbf{AQ}=0,$

or equivalently
$\mathbf{Q^{T}}\mathbf{AQN}=\mathbf{N}\mathbf{Q}^{T}\mathbf{AQ}$. Because
$\mathbf{N}$ is a diagonal matrix with different diagonal entries, so
$\mathbf{Q}^{T}\mathbf{AQ}$ is also diagonal, thus $\mathbf{Q}\in E$.
Therefore, for both (4.26) and (4.29), the $\omega-$limit set should be a
subset of $E$.

Because $O(n)$ is compact, so the $\omega-$limit set should be connected. Thus
$\mathbf{Q}(t)$ converges to a connected component of $E$. However, by Lemma
4.2, $E$ consists of only finite points (in fact $n!2^{n}$ of them), so there
exists $\mathbf{Q}^{*}\in E$ such that (5.1) holds. ∎

###  5.2. Stable manifolds

From now on, we focus on (2.9) instead of the general case (4.26) and proceed
to discuss the stable manifolds of each equilibrium. Even though a complete
characterization of the stable manifolds for the Oja-Brockett flow [4]
($\mathbf{S=0}$ in (4.29)) is still an open problem [3], the stable manifolds
of (2.9) will be fully determined in this section. The critical technique
which facilitates our characterization of the stable manifolds is the
following solution formula.

####  5.2.1. The solution formula by the Cholesky decomposition

First of all, we rewrite (2.9) as

(5.2)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{Q}}&=\mathbf{AQ}-\mathbf{QT(Q)},\\\
\mathbf{Q}(0)&=\mathbf{Q_{0}}\in O(n),\end{split}\right.$

where $\mathbf{T}$ is defined as

(5.3) $\displaystyle\mathbf{T(Q)}$
$\displaystyle:=\mathbf{Q}^{T}\mathbf{AQ}+\sum_{k=1}^{n}\sum_{j=1}^{k-1}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}$
(5.4) $\displaystyle=\mathrm{diag}\\{\mathbf{q_{1}}\cdot\mathbf{Aq_{1}},\
\mathbf{q_{2}}\cdot\mathbf{Aq_{2}},\ ...,\
\mathbf{q_{n}}\cdot\mathbf{Aq_{n}}\\}+2\sum_{k=1}^{n}\sum_{j=1}^{k-1}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}.$

So $\mathbf{T}$ is an upper triangular matrix. This property determines that
the evolution of $\mathbf{q_{i}}$ only depends on $\mathbf{q_{j}},\ 1\leq j<i$
but is independent with $\mathbf{q_{k}},\ i<k\leq n$. This independence is
helpful to the description of the stable manifolds.

We first derive a representation formula of $\mathbf{Q}$ from (5.2).

######  Lemma 5.1.

Consider (2.9) (or the alternative form (5.2)). Then there exists an upper
triangular matrix $\mathbf{G}(t)$ such that

(5.5) $\displaystyle\mathbf{Q}(t)=e^{\mathbf{A}t}\mathbf{Q_{0}}\mathbf{G}(t).$

Moreover, diagonal entries of $\mathbf{G}(t)$ are positive and $\mathbf{G}(t)$
satisfies

(5.6)
$\displaystyle\mathbf{G}(t)\mathbf{G}^{T}(t)=\mathbf{Q_{0}}^{T}e^{-2\mathbf{A}t}\mathbf{Q_{0}}.$

See Section 6 for the proof of this lemma.

Notice that (6.27) should just be understood as an iteration formula coupled
with $\mathbf{Q}$ instead of an explicit closed formula of the solution.
However, $\mathbf{G}(t)$ is uniquely determined by (5.6), since (5.6) implies
that $(\mathbf{G}^{T})^{-1}(t)$ is the Cholesky decomposition of
$\mathbf{Q_{0}}^{T}e^{2\mathbf{A}t}\mathbf{Q_{0}}$, i.e.

(5.7)
$\displaystyle(\mathbf{G}^{T})^{-1}(t)\mathbf{G}^{-1}(t)=\mathbf{Q_{0}}^{T}e^{2\mathbf{A}t}\mathbf{Q_{0}}.$

####  5.2.2. Characterization of stable manifolds

We first introduce the definition of stable manifolds for readers’
convenience.

######  Definition 5.1.

Let $M$ be a smooth manifold and $a\in M$ is an equilibrium of a smooth vector
vector field $X:M\to TM$. The stable manifold of $a$ (or the inset of $a$) is

(5.8) $\displaystyle W^{s}(a):=\\{x\in M:\ L_{\omega}(x)=\\{a\\}\\}.$

Here $L_{\omega}(x)$ is the $\omega$ limit set of $x$.

We will exploit (5.5) in this section to determine the stable manifolds of
each equilibria of $\mathcal{E}$. For $\mathbf{Q}\in\mathbb{R}^{n\times n}$
which is invertible, we recursively define $\sigma(i),\ i=1,2,...,n$ which is
a permutation of $\\{1,2,...,n\\}$:

(5.9) $\displaystyle\sigma_{1}(\mathbf{Q}):=\inf\limits_{1\leq k\leq n}\\{k:\
q_{1,k}\neq 0\\},$ (5.10)
$\displaystyle\sigma_{m}(\mathbf{Q}):=\inf\limits_{1\leq k\leq n}\\{k:\
\mathrm{det}(\mathbf{Q}[1,2,3,...,m;\ \sigma(1),\ \sigma(2),\ ...,\ k])\neq
0\\}.$

For instance, if

(5.11) $\displaystyle\mathbf{M}=\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&1\end{pmatrix},$

then $\sigma_{1}(\mathbf{M})=2,\ \sigma_{2}(\mathbf{M})=1,\
\sigma_{3}(\mathbf{M})=3$.

A complete characterization of the stable manifolds of (2.9) will be derived
upon a dedicate analysis of properties of $\sigma$. We state the following
lemmas for $\sigma$ whose proofs are attached to the Appendix (see Section
6.2). In the following sections, by default, if the proof is not presented
after a lemma then it is summarized in Section 6.2.

######  Lemma 5.2.

Let $\mathbf{Q}\in\mathbb{R}^{n\times n}$ be invertible. Then for any $1\leq
m,j\leq n$, we have

(5.12) $\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right)+\mathrm{card}(\\{k:\ \sigma_{k}(\mathbf{Q})>j,\ 1\leq k\leq
m\\})=m.$

See Section 6 for proof of this lemma. From Lemma 5.2, we derive the following
lemma which will be directly used to determine the stable manifolds:

######  Lemma 5.3.

Let $\mathbf{Q}\in\mathbb{R}^{n\times n}$ be invertible. Define permutation
$\sigma_{m},\ m=1,2,...n$ as in (5.9), and

(5.13) $\displaystyle\mathbf{i_{m}}:=(i_{m,1},\ i_{m,2},\ ...,\ i_{m,m}),1\leq
m\leq n$

where $\mathbf{i_{m}}$ is a permutation of
$\\{\sigma_{k}(\mathbf{Q_{0}})\\}_{1\leq k\leq m}$ such that
$i_{m,1}<i_{m,2}<i_{m,3}<...<i_{m,m}$. Then

(5.14) $\displaystyle\mathbf{Q}[m;\
1,2,3,...,\sigma_{m}(\mathbf{Q})]=\left\\{\begin{split}&\mathbf{c_{m}}\mathbf{Q}[1,2,...,m-1;\
1,2,3,...,\sigma_{m}({\mathbf{Q}})]+(0,0,0,...,0,z_{m}),\ 2\leq m\leq n,\\\
&(0,0,...,0,z_{m}),\ m=1.\end{split}\right.$

where $\mathbf{c_{m}}=(c_{m,1},\ c_{m,2},\ ...,\ c_{m,m-1}),2\leq m\leq n$ is
the unique solution of

(5.15) $\displaystyle\mathbf{Q}[m;\
\mathbf{i_{m-1}}]=\mathbf{c_{m}}\mathbf{Q}[1,2,...,m-1;\ \mathbf{i_{m-1}}],$

and $z_{m}$ is defined as

(5.16) $\displaystyle
z_{m}=\left\\{\begin{split}&\dfrac{\det(\mathbf{Q}[1,2,3,...,m;\
\mathbf{i_{m}}])}{\det(\mathbf{Q}[1,2,3,...,m-1;\ \mathbf{i_{m-1}}])},\ 2\leq
m\leq n,\\\ &q_{1,\sigma_{1}(\mathbf{Q})},\ m=1.\end{split}\right.$

See Section 6 for proof of this lemma. By Lemma 5.3, we have the following
theorem on stable manifolds:

######  Theorem 5.2.

Suppose that $\mathbf{Q}(t),t\geq 0$ solves (2.9) with initial value
$\mathbf{Q_{0}}\in O(n)$. Define permutation $\sigma_{m},\ m=1,2,...n$ as in
(5.9), and

(5.17) $\displaystyle\mathbf{i_{m}}:=(i_{m,1},\ i_{m,2},\ ...,\ i_{m,m}),1\leq
m\leq n$

where $\mathbf{i_{m}}$ is a permutation of
$\\{\sigma_{k}(\mathbf{Q_{0}})\\}_{1\leq k\leq m}$ such that
$i_{1}<i_{2}<i_{3}<...<i_{m}$. Then

(5.18)
$\displaystyle\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}(t)\to\mathrm{sgn}(z_{m})\mathbf{e_{m}},\
m=1,2,...,n,$

as $t\to\infty$. Here $z_{m}$ is defined in $\eqref{def:z_m}$

###### Proof.

We prove by induction. For $m=1$, we know that $q_{1,j}(0)=0$ for all $1\leq
j<\sigma_{1}(\mathbf{Q_{0}})$ while $q_{1,\sigma_{1}(\mathbf{Q_{0}})}\neq 0$.
By Lemma 5.1, we have

(5.19) $\displaystyle
q_{1,\sigma_{1}(\mathbf{Q_{0}})}(t)=\dfrac{e^{\lambda_{1}t}q_{1,\sigma_{1}(\mathbf{Q_{0}})}(0)}{f_{\sigma_{1}(\mathbf{Q_{0}})}(t)}.$

Here $f$ is defined in (6.25). By Theorem 5.1, we know that
$q_{1,\sigma_{1}(\mathbf{Q_{0}})}(t)$ converges to one of $1,\ -1$ and $0$.
Because
$\mathbf{q_{\sigma_{1}(\mathbf{Q_{0}})}}\cdot\mathbf{Aq_{\sigma_{1}(\mathbf{Q_{0}})}}\leq\lambda_{1}$,
so

(5.20)
$\displaystyle\dfrac{e^{\lambda_{1}t}}{f_{\sigma_{1}(\mathbf{Q_{0}})}(t)}=e^{\int_{0}^{t}\lambda_{1}-\mathbf{q_{\sigma_{1}(\mathbf{Q_{0}})}}(s)\cdot\mathbf{Aq_{\sigma_{1}(\mathbf{Q_{0}})}}(s)\mathrm{d}s}\geq
e^{0}=1.$

Thus $q_{1,\sigma_{1}(\mathbf{Q_{0}})}(t)$ does not converge to 0, and it
converges to either $1$ or $-1$ which depends on the sign of the initial value
$q_{1,\sigma_{1}(\mathbf{Q_{0}})}$. This implies

(5.21)
$\displaystyle\mathbf{q_{\sigma_{1}(\mathbf{Q_{0}})}}\to\mathrm{sgn}(q_{1,\sigma_{1}(\mathbf{Q_{0}})})\mathbf{e_{1}}$

as $t\to\infty$. Thus the claim holds for $m=1$.

Now suppose that (5.18) holds for all $j=1,2,...,m-1$, we prove it for $j=m$.
By (5.14) in Lemma 5.3 and (6.24) in Lemma 5.1, we have

$\displaystyle q_{m,\sigma_{m}(\mathbf{Q_{0}})}(t)$
$\displaystyle=e^{\lambda_{m}t}\sum_{l=1}^{\sigma_{m}(\mathbf{Q_{0}})}q_{m,l}(0)g_{l,\sigma_{m}(\mathbf{Q_{0}})}(t)$
$\displaystyle=e^{\lambda_{m}t}\mathbf{Q_{0}}[m;\
1,2,3,...,\sigma_{m}(\mathbf{Q_{0}})]\mathbf{G}(t)[1,2,...,\sigma_{m}(\mathbf{Q_{0}});\
\sigma_{m}(\mathbf{Q_{0}})]$
$\displaystyle=e^{\lambda_{m}t}\mathbf{c_{m}}\mathbf{Q_{0}}[1,2,...,m-1;\
1,2,3,...,\sigma_{m}({\mathbf{Q_{0}}})]\mathbf{G}(t)[1,2,...,\sigma_{m}(\mathbf{Q_{0}});\
\sigma_{m}(\mathbf{Q_{0}})]$
$\displaystyle+\dfrac{e^{\lambda_{m}t}z_{m}}{f_{\sigma_{m}(\mathbf{Q_{0}})}(t)}$
(5.22)
$\displaystyle=\sum_{j=1}^{m-1}e^{(\lambda_{m}-\lambda_{j})t}c_{m,j}q_{j,\sigma_{m}(\mathbf{Q_{0}})}(t)+\dfrac{e^{\lambda_{m}t}z_{m}}{f_{\sigma_{m}(\mathbf{Q_{0}})}(t)}.$

Here $\mathbf{c_{m}}$ and $z_{m}$ are defined in Lemma 5.3. By the hypothesis
of induction, we know that $\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}$ does not
converge to $\pm\mathbf{e_{1}},\ \pm\mathbf{e_{1}},\ ,...\pm\mathbf{e_{m-1}}$,
so

(5.23) $\displaystyle q_{j,\sigma_{m}(\mathbf{Q_{0}})}(t)\to 0\ \mathrm{as}\
t\to\infty.$

Thus

(5.24)
$\displaystyle\sum_{j=1}^{m-1}e^{(\lambda_{m}-\lambda_{j})t}c_{m,j}q_{j,\sigma_{m}(\mathbf{Q_{0}})}(t)\to
0\ \mathrm{as}\ t\to\infty.$

Again, by Theorem 5.1, $q_{m,\sigma_{m}(\mathbf{Q_{0}})}$ converges to one of
$1,-1$ and 0. Now we prove that it converges to either $1$ or $-1$. Otherwise
$\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}$ can only converge to
$\pm\mathbf{e_{p}},p\geq m+1$, hence
$\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}\cdot\mathbf{A\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}}\to\lambda_{p}<\lambda_{m}$,
which results in

$\displaystyle\dfrac{e^{\lambda_{m}t}}{f_{\sigma_{m}(\mathbf{Q_{0}})}(t)}=e^{\int_{0}^{t}\lambda_{m}-\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}(s)\cdot\mathbf{A\mathbf{q_{\sigma_{m}(\mathbf{Q_{0}})}}}(s)\mathrm{d}t}\to\infty\
\mathrm{as}\ t\to\infty,$

a contradiction. Thus
$\dfrac{e^{\lambda_{m}t}z_{m}}{f_{\sigma_{m}(\mathbf{Q_{0}})}(t)}$ converges
to either $1$ or $-1$ and

(5.25) $\displaystyle
q_{m,\sigma_{m}(\mathbf{Q_{0}})}(t)=\dfrac{e^{\lambda_{m}t}z_{m}}{f_{\sigma_{m}(\mathbf{Q_{0}})}(t)}\to\mathrm{sgn}(z_{m}).$

Thus (5.18) holds for $m$. By induction, it holds for all $m=1,2,...,n$. ∎

As a direct corollary, we can see that $\sigma_{m}(\mathbf{Q}(t)),m=1,2,...,n$
does not change for all $t>0$.

######  Corollary 5.1.

Suppose that $\mathbf{Q}(t)$ is a solution of (2.9). Then for all
$m=1,2,...,n$, we have

$\displaystyle\sigma_{m}(\mathbf{Q}(t))=\sigma_{m}(\mathbf{Q}(0)),t>0.$

###### Proof.

For arbitrary $t_{0}>0$, consider $\mathbf{R}(t):=\mathbf{Q}(t+t_{0}),t>0$.
Then $\mathbf{R}(t),t>0$ is the unique solution to (2.9) with initial value
$\mathbf{Q}(t_{0})$. Then by Theorem 5.2, for all $m=1,2,...,n$, we have

$\displaystyle\mathbf{r}_{\sigma_{m}(\mathbf{R}(0))}(t)\to\mathrm{sgn}(z_{m}^{\prime})\mathbf{e_{m}}\
\mathrm{as}\ t\to\infty.$

Remember that

$\displaystyle\mathbf{q}_{\sigma_{m}(\mathbf{Q}(0))}(t)\to\mathrm{sgn}(z_{m})\mathbf{e_{m}}\
\mathrm{as}\ t\to\infty.$

Moreover, $\mathbf{R}(t)=\mathbf{Q}(t+t_{0})$, so

$\displaystyle\lim\limits_{t\to\infty}\mathbf{R}(t)=\lim\limits_{t\to\infty}\mathbf{Q}(t).$

Thus

$\displaystyle\lim\limits_{t\to\infty}\mathbf{r}_{\sigma_{m}(\mathbf{R}(0))}(t)=\lim\limits_{t\to\infty}\mathbf{q}_{\sigma_{m}(\mathbf{Q}(0))}(t).$

So

$\displaystyle\sigma_{m}(\mathbf{R}(0))=\sigma_{m}(\mathbf{Q}(0)).$

By orthogonality of $\mathbf{Q}$ and $\mathbf{R}(t)=\mathbf{Q}(t+t_{0})$, we
know

$\displaystyle\sigma_{m}(\mathbf{Q}(t_{0}))=\sigma_{m}(\mathbf{R}(0))=\sigma_{m}(\mathbf{Q}(0)).$

∎

####  5.2.3. An example

To illustrate Theorem 5.2, we discuss an example here. Consider

$\displaystyle\mathbf{Q_{1}}=\begin{pmatrix}0&\dfrac{\sqrt{2}}{2}&-\dfrac{\sqrt{3}}{3}&\dfrac{\sqrt{6}}{6}\\\
0&\dfrac{\sqrt{2}}{2}&\dfrac{\sqrt{3}}{3}&-\dfrac{\sqrt{6}}{6}\\\
-\dfrac{\sqrt{2}}{2}&0&\dfrac{\sqrt{6}}{6}&\dfrac{\sqrt{3}}{3}\\\
\dfrac{\sqrt{2}}{2}&0&\dfrac{\sqrt{6}}{6}&\dfrac{\sqrt{3}}{3}\end{pmatrix}.$

Suppose that $\mathbf{Q}(t),t\geq 0$ solves (2.9) with initial value
$\mathbf{Q_{1}}$. According (5.9), we have

$\displaystyle\sigma_{1}(\mathbf{Q_{1}})=2,\ \sigma_{2}(\mathbf{Q_{1}})=3,\
\sigma_{3}(\mathbf{Q_{1}})=1,\ \sigma_{4}(\mathbf{Q_{1}})=4.$

by Theorem 5.2, we have

$\displaystyle\mathbf{q_{1}}(t)$
$\displaystyle\to\mathrm{sgn}(z_{3})\mathbf{e_{3}},\
\mathbf{q_{2}}(t)\to\mathrm{sgn}(z_{1})\mathbf{e_{1}},$
$\displaystyle\mathbf{q_{3}}(t)$
$\displaystyle\to\mathrm{sgn}(z_{2})\mathbf{e_{2}},\
\mathbf{q_{4}}(t)\to\mathrm{sgn}(z_{4})\mathbf{e_{4}}.$

A direction computation according to (5.16) yields
$z_{1}=\sqrt{2}/2,z_{2}=2\sqrt{3}/3,z_{3}=-\sqrt{2}/2,z_{4}=2\sqrt{3}$. Thus
as $t\to\infty$,

$\displaystyle\mathbf{q_{1}}(t)\to-\mathbf{e_{3}},\
\mathbf{q_{2}}(t)\to\mathbf{e_{1}},\ \mathbf{q_{3}}(t)\to\mathbf{e_{2}},\
\mathbf{q_{4}}(t)\to\mathbf{e_{4}}.$

###  5.3. Convergence rate to the stable equilibria

Applying the technique of linearization, we can prove that an equilibrium
$\mathbf{Q}^{*}\in E$ (defined in (4.23)) is stable if and only if

(5.26)
$\displaystyle\mathbf{Q}^{*}=\mathrm{diag}\\{\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},...,\varepsilon_{n}\\},\
\varepsilon_{i}=\pm 1,\ i=1,2,...,n.$

See Section 6.3 for details. By Theorem 5.2, the solution of (2.9) with
initial value $\mathbf{Q_{0}}$ converges to a stable equilibrium if and only
if

(5.27) $\displaystyle\sigma_{k}(\mathbf{Q_{0}})=k,\ k=1,2,...,n,$

or equivalently, every leading principle submatrix is invertible, i.e.

(5.28) $\displaystyle\mathrm{det}(\mathbf{Q_{0}}[1,2,...,k;\ 1,2,...,k])\neq
0,\ k=1,2,...,n.$

Now we proceed to prove the exponential convergence to a stable equilibrium.
The main idea of this proof is similar to the one of power method: utilizing
$(\lambda_{i}/\lambda_{1})^{n}\to 0$ (here $e^{(\lambda_{i}-\lambda_{1})t}\to
0$), we derive convergence rate.

######  Theorem 5.3.

(exponential convergence) Suppose that $\mathbf{Q_{0}}$ is in the stable
manifold of a stable equilibrium, i.e., (5.28) holds. Let $\mathbf{Q}(t),t\geq
0$ be the solution of (2.9) with initial data $\mathbf{Q_{0}}$. Define

(5.29) $\displaystyle\nu_{k}$ $\displaystyle:=\min\\{\lambda_{1}-\lambda_{2},\
\lambda_{2}-\lambda_{3},\ ...,\ \lambda_{k}-\lambda_{k+1}\\}>0,\
k=1,2,...,n-1,$ $\displaystyle\nu_{n}$ $\displaystyle:=\nu_{n-1}.$

Then there exists a constant $C>0$ such that

(5.30) $\displaystyle|q_{i,j}^{2}(t)-\delta_{i,j}|\leq Ce^{-2\nu_{i\wedge
j}t},\ i,j=1,2,...,n.$

Here $\delta_{i,j},\ i,j=1,2,...,n$ are Kronecker functions.

###### Proof.

In the following proof, $C$ just represents for a general constant which needs
not to be invariant through the whole proof.

We will prove by induction. As we pointed out in introduction, alignment of
$\mathbf{q_{1}}$ is the fastest, so we consider it first.

We first prove (5.30) for $q_{k,1},q_{1,k}\ k=1,2,...,n$. By the computation
in Lemma 5.1, we have

$\displaystyle\mathbf{q_{1}}(t)=\dfrac{e^{\mathbf{A}t}\mathbf{q_{1,0}}}{\|e^{\mathbf{A}t}\mathbf{q_{1,0}}\|}.$

Remember that $\mathbf{Q_{0}}$ satisfies (5.28), so $q_{1,1}(0)\neq 0$ and

$\displaystyle\dfrac{q_{k,1}(t)}{q_{1,1}(t)}=\dfrac{q_{k,1}(0)}{q_{1,1}(0)}e^{(\lambda_{k}-\lambda_{1})t},\
k=1,2,...,n.$

So $|q_{k,1}(t)|\leq Ce^{(\lambda_{k}-\lambda_{1})t}|q_{1,1}(t)|$. Summing up
$k$ from $2$ to $n$ yields

$\displaystyle 1-|q_{1,1}(t)|^{2}$
$\displaystyle=\sum_{k=2}^{n}q_{k,1}^{2}(t)=C\sum_{k=2}^{n}e^{2(\lambda_{k}-\lambda_{1})t}|q_{1,1}(t)|^{2}\leq
Ce^{2(\lambda_{2}-\lambda_{1})t}|q_{1,1}(t)|^{2}.$

Thus

$\displaystyle
1-|q_{1,1}^{2}(t)|\leq\dfrac{Ce^{2(\lambda_{2}-\lambda_{1})t}}{1+Ce^{2(\lambda_{2}-\lambda_{1})t}}\leq
Ce^{2(\lambda_{2}-\lambda_{1})t}.$

Thus by orthogonality, we have

$\displaystyle|q_{1,k}(t)|^{2}\leq 1-|q_{1,1}(t)|^{2}\leq
Ce^{2(\lambda_{2}-\lambda_{1})t},\ k=2,...,n,$
$\displaystyle|q_{k,1}(t)|^{2}\leq 1-|q_{1,1}(t)|^{2}\leq
Ce^{2(\lambda_{2}-\lambda_{1})t},\ k=2,...,n.$

Thus (5.30) holds for $i,j$ such that $i\wedge j=1$.

After alignment of first $k-1$ components, we consider that of
$\mathbf{q_{k}}$. Suppose that (5.30) holds for $i,j$ such that $i\wedge
j=1,2,...,k-1,k\geq 2$, we prove that (5.30) also holds for $i,j$ such that
$i\wedge j=k$. By (6.24), we have for $m=1,2,...,n$

$\displaystyle
e^{-\lambda_{m}t}q_{m,k}(t)=\sum_{j=1}^{n}g_{j,k}(t)q_{m,j}(0)=\mathbf{Q_{0}}[m;\
1,2,...,k]\mathbf{G}(t)[1,2,...,k;\ m].$

Here $\mathbf{G}(t)$ is defined in Lemma 5.1 (see (5.5)). Remember that
$\mathbf{Q_{0}}[1,2,...,k;\ 1,2,...,k]$ is invertible, so there exist
constants $c_{m,j},j=1,2,...,k$ such that for all $k+1\leq m\leq n$,

$\displaystyle\mathbf{Q_{0}}[m;\
1,2,...,k]=\sum_{j=1}^{k}c_{m,j}\mathbf{Q_{0}}[j;\ 1,2,...,k].$

Thus for $k+1\leq m\leq n$,

$\displaystyle e^{-\lambda_{m}t}q_{m,k}(t)$
$\displaystyle=\sum_{j=1}^{k}c_{m,j}\mathbf{Q_{0}}[j;\
1,2,...,k]\mathbf{G}(t)[1,2,...,k;\ m]$
$\displaystyle=\sum_{j=1}^{k}c_{m,j}e^{-\lambda_{j}t}q_{j,k}(t),$

so by induction hypothesis (i.e., (5.30) holds for $i\wedge j=k-1$)

$\displaystyle|q_{m,k}(t)|\leq
Ce^{(\lambda_{m}-\lambda_{k})t}|q_{k,k}(t)|+Ce^{-\nu_{k-1}t}.$

Summing up from $m=k+1$ to $n$ yields

$\displaystyle
1-\sum_{j=1}^{k}q_{j,k}^{2}(t)-q_{k,k}^{2}(t)=\sum_{m=k+1}^{n}|q_{m,k}(t)|^{2}\leq
Ce^{2(\lambda_{k+1}-\lambda_{k})t}|q_{k,k}(t)|^{2}+Ce^{-2\nu_{k-1}t}.$

Thus

$\displaystyle
1-q_{k,k}^{2}(t)\leq\sum_{j=1}^{k}q_{j,k}^{2}(t)+Ce^{2(\lambda_{k+1}-\lambda_{k})t}q_{k,k}^{2}(t)+Ce^{-2\nu_{k-1}t}\leq
Ce^{2(\lambda_{k+1}-\lambda_{k})t}q^{2}_{k,k}(t)+Ce^{-2\nu_{k-1}t},$

hence

$\displaystyle
1-q_{k,k}^{2}(t)\leq\dfrac{Ce^{2\nu_{k-1}t}+Ce^{2(\lambda_{k+1}-\lambda_{k})t}}{1+e^{2(\lambda_{k+1}-\lambda_{k})t}}\leq
Ce^{-2\min\\{\nu_{k-1},\lambda_{k}-\lambda_{k+1}\\}t}=Ce^{-2\nu_{k}t}.$

So for any $j\geq k+1$, by orthogonality,

$\displaystyle q_{j,k}^{2}(t)\leq 1-q^{2}_{k,k}(t)\leq Ce^{-2\nu_{k}t},$
$\displaystyle q_{k,j}^{2}(t)\leq 1-q^{2}_{k,k}(t)\leq Ce^{-2\nu_{k}t}.$

Thus (5.30) also holds for $i,j$ such that $i\wedge j=k+1$. So by induction,
(5.30) holds for $i,j=1,2,...,n$. ∎

###  5.4. The case of $p<n$

Now we discuss the case $p<n$. In this case, the ODE system is reduced to

(5.31)
$\displaystyle\left\\{\begin{split}\mathbf{\dot{q_{1}}}&=\mathbf{Aq_{1}}-\mathbf{(q_{1}\cdot
Aq_{1})q_{1}},\\\ \mathbf{\dot{q_{j}}}&=\mathbf{Aq_{j}}-\mathbf{(q_{j}\cdot
Aq_{j})q_{j}}-2\sum_{i=1}^{j-1}\mathbf{(q_{i}\cdot Aq_{j})q_{i}},\ j=2,\ 3,\
...,\ p.\\\ \mathbf{q_{i}}(0)&=\mathbf{q_{i,0}},\ i=1,\ 2,\ ...,\
p.\end{split}\right.$

The initial value $\mathbf{Q_{0}}\in\mathbb{R}^{n\times p}$ satisfies
$\mathbf{Q_{0}}^{T}\mathbf{Q_{0}}=\mathbf{I_{p}}$. Denote the solution of
(5.31) with initial value $\mathbf{Q_{0}}$ as $\mathbf{Q}(t),t>0$.

To derive the convergence property of $\mathbf{Q}(t)$, we just need to
complete $\mathbf{Q_{0}}$, make it in $O(n)$. Now let
$\mathbf{R_{0}}\in\mathbb{R}^{n\times(n-p)}$ satisfy that

(5.32) $\displaystyle\mathbf{Q_{1}}=[\mathbf{Q_{0}},\mathbf{R_{0}}]\in O(n).$

Let $\mathbf{Q_{1}}(t)$ be the solution of (2.9) with initial value
$\mathbf{Q_{1}}$. Notice that (5.31) is exactly the same as the first $p$
components of (1.5) and (2.9), thus it is also semi-decoupling. Therefore, by
uniqueness the solution,

$\displaystyle\mathbf{Q}(t)=\mathbf{Q_{1}}[1,2,...,n;1,2,...,p](t),\ t>0,$

i.e., the first $p$ columns of $\mathbf{Q_{1}}$ and $\mathbf{Q}$ are the same.
So the same convergence properties in Theorem 5.2 and Theorem 5.3 hold for
$\mathbf{Q}$.

## Acknowledgement

Jian-Guo Liu was supported in part by the National Science Foundation (NSF)
under award DMS-2106988.

## References

  * [1] Tomáš Bárta, Ralph Chill, and Eva Fašangová.  Every ordinary differential equation with a strict lyapunov function is a gradient system.  Monatshefte für Mathematik, 166(1):57–72, 2012.
  * [2] Sergio Bittanti, Alan J Laub, and Jan C Willems.  The Riccati Equation.  Springer Science & Business Media, 2012.
  * [3] Vincent D Blondel, Alexandre Megretski, and Vincent DD Blondel.  Unsolved problems in mathematical systems and control theory.  Princeton University Press Princeton, NJ, 2004.
  * [4] Roger W Brockett.  Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems.  Linear Algebra and its applications, 146:79–91, 1991.
  * [5] Thomas L Gilbert.  A lagrangian formulation of the gyromagnetic equation of the magnetization field.  Phys. Rev., 100:1243, 1955.
  * [6] Thomas L Gilbert.  A phenomenological theory of damping in ferromagnetic materials.  IEEE transactions on magnetics, 40(6):3443–3449, 2004.
  * [7] Uwe Helmke and John B Moore.  Optimization and dynamical systems.  Springer Science & Business Media, 2012.
  * [8] L Landau and W Lifshitz.  Theory of the dispersion of magnetic permeability in ferromagnetic bodies. phys. z. sowjetunion 8, 153 (1935). a. hubert and r. schäfer, magnetics domains, 1998.
  * [9] Joseph LaSalle.  Some extensions of liapunov’s second method.  IRE Transactions on circuit theory, 7(4):520–527, 1960.
  * [10] Chris Junchi Li, Mengdi Wang, Han Liu, and Tong Zhang.  Diffusion approximations for online principal component estimation and global convergence.  Advances in Neural Information Processing Systems, 30, 2017.
  * [11] Erkki Oja.  Simplified neuron model as a principal component analyzer.  Journal of mathematical biology, 15(3):267–273, 1982.
  * [12] Erkki Oja.  Principal components, minor components, and linear neural networks.  Neural networks, 5(6):927–935, 1992.
  * [13] Erkki Oja and Juha Karhunen.  On stochastic approximation of the eigenvectors and eigenvalues of the expectation of a random matrix.  Journal of mathematical analysis and applications, 106(1):69–84, 1985.
  * [14] Terence D Sanger.  Optimal unsupervised learning in a single-layer linear feedforward neural network.  Neural networks, 2(6):459–473, 1989.
  * [15] Wei-Yong Yan, Uwe Helmke, and John B Moore.  Global analysis of oja’s flow for neural networks.  IEEE Transactions on Neural Networks, 5(5):674–683, 1994.
  * [16] Shintaro Yoshizawa, Uwe Helmke, and Konstantin Starkov.  Convergence analysis for principal component flows.  2001\.

##  6\. Appendix

###  6.1. Riemannian manifolds and the Stiefel manifold

We denote the tangent space at $\mathbf{m}$ on manifold $\mathcal{M}$ as
$T_{\mathbf{m}}\mathcal{M}$, the tangent vector field on $\mathcal{M}$ as
$\Gamma(T\mathcal{M})$. The tangent bundle (the disjoint union of the tangent
spaces) is denoted as $T\mathcal{M}$.

######  Definition 6.1.

(Riemannian manifolds) Suppose that $\mathcal{M}$ is a smooth manifold. A
Riemannian manifold $(\mathcal{M},g)$ is a smooth mainfold equipped with an
inner product $g_{\mathbf{m}}$ on $T_{\mathbf{m}}\mathcal{M}$ at each
$\mathbf{m}\in\mathcal{M}$. Moreover, for any tangent vector field
$\dot{\mathbf{x}}$ and $\dot{\mathbf{y}}$, the function

(6.1)
$\displaystyle\langle\dot{\mathbf{x}}(\mathbf{m}),\dot{\mathbf{y}}(\mathbf{m})\rangle_{g_{\mathbf{m}}}:\
\mathcal{M}\to\mathbb{R}$

is smooth.

Given a Riemannian metric $g$ on $M$, the gradient of a smooth function
$\mathcal{E}$ on $\mathcal{M}$ is defined as

######  Definition 6.2.

(the gradient on the Riemannian manifold) A tangent vector field
$\nabla_{g}\mathcal{E}$ on $M$ is called the gradient of $\mathcal{E}$ w.r.t.
the metric $g$ if for every tangent vector field $\dot{\mathbf{x}}$ on $M$,

(6.2) $\displaystyle\langle\mathcal{E}^{\prime},\
\dot{\mathbf{x}}\rangle_{F}=\langle\nabla_{g}\mathcal{E},\
\dot{\mathbf{x}}\rangle_{g}.$

Here $\mathcal{E}^{\prime}$ is the derivative of $\mathcal{E}$, which is a
cotangent vector field.

Now we consider the Stiefel manifold with the Euclidean metric $g_{e}$, under
the global coordinate $\mathbf{Q}\in\mathbb{R}^{n\times n}$. We first
introduce several important properties of $O(n)$. See [7] for more details.

######  Lemma 6.1.

(the Stiefel manifold) The Stiefel manifold $O(n)$ is a smooth, compact
manifold of dimension $n(n-1)/2$. The tangent space at $\mathbf{Q}$ is given
by

(6.3) $\displaystyle T_{\mathbf{Q}}O(n)=\\{\mathbf{Q}\bm{\Omega}\ |\
\bm{\Omega}\in\mathbb{R}^{n\times n},\
\bm{\Omega}+\bm{\Omega}^{T}=\mathbf{0}\\},$

while the normal space at $\mathbf{Q}$ is given by

(6.4) $\displaystyle T_{\mathbf{Q}}O(n)^{\perp}=\\{\mathbf{Q}\bm{\Omega}\ |\
\bm{\Omega}\in\mathbb{R}^{n\times n},\ \bm{\Omega}=\bm{\Omega}^{T}\\}.$

See [7] for proof of Lemma 6.1. By Lemma 6.1, we can prove that for any
$\mathbf{M}\in\mathbb{R}^{n\times n}$, the projection on the tangent spaces
and the normal spaces are respectively:

(6.5)
$\displaystyle\mathcal{P}_{T_{\mathbf{Q}}O(n)}\mathbf{M}:=\dfrac{1}{2}(\mathbf{M}-\mathbf{Q}\mathbf{M}^{T}\mathbf{Q}),\
\mathcal{P}_{T_{\mathbf{Q}}O(n)^{\perp}}\mathbf{M}:=\dfrac{1}{2}(\mathbf{M}+\mathbf{Q}\mathbf{M}^{T}\mathbf{Q}).$

######  Lemma 6.2.

(Gradient on $(O(n),g_{e})$) Suppose that
$\varphi(\mathbf{Q}):O(n)\to\mathbb{R}$ is a restriction of a smooth function
(still denoted as $\varphi:\mathbb{R}^{n\times n}\to\mathbb{R}$) on $O(n)$.
Then gradient of $\varphi$ w.r.t. $g_{e}$ at point $\mathbf{Q}$ is given by

(6.6)
$\displaystyle\nabla_{g_{e}}\varphi:=\mathcal{P}_{T_{\mathbf{Q}}O(n)}(\nabla\varphi)=\dfrac{1}{2}(\nabla\varphi-\mathbf{Q}(\nabla\varphi)^{T}\mathbf{Q}),$

here $\nabla\varphi\in\mathbb{R}^{n\times n}$ is the gradient of $\varphi$ in
$\mathbb{R}^{n\times n}$, i.e.
$(\nabla\varphi)_{ij}=\dfrac{\partial\varphi}{\partial q_{ij}}$.

###### Proof.

We just need to prove that for any $\mathbf{Q}\in O(n)$ and any tangent vector
at $\mathbf{Q}$, (6.2) holds. By Lemma 6.1, a tangent vector at $\mathbf{Q}$
can be represented as $\mathbf{Q}\bm{\Omega}$ where $\bm{\Omega}$ is skew-
symmetric. Therefore, for any $\bm{\Omega}$ that is skew-symmetric, we have

(6.7)
$\displaystyle\langle\nabla\varphi,\mathbf{Q}\bm{\Omega}\rangle_{F}=\langle\mathcal{P}_{T_{\mathbf{Q}}O(n)}(\nabla\varphi),\
\mathbf{Q}\bm{\Omega}\rangle_{F}+\langle\mathcal{P}_{T_{\mathbf{Q}}O(n)^{\perp}}(\nabla\varphi),\
\mathbf{Q}\bm{\Omega}\rangle_{F}.$

Notice that

$\displaystyle\langle\mathcal{P}_{T_{\mathbf{Q}}O(n)^{\perp}}(\nabla\varphi),\
\mathbf{Q}\bm{\Omega}\rangle_{F}$
$\displaystyle=\langle\mathbf{Q}(\mathbf{Q}^{T}\nabla\varphi+(\nabla\varphi)^{T}\mathbf{Q}),\
\mathbf{Q}\bm{\Omega}\rangle_{F}$
$\displaystyle=\langle\mathbf{Q}^{T}\nabla\varphi+(\nabla\varphi)^{T}\mathbf{Q},\
\bm{\Omega}\rangle_{F}$ $\displaystyle=0,$

since $\langle\mathbf{M},\ \mathbf{N}\rangle_{F}=0$ if $\mathbf{M}$ is
symmetric while $\mathbf{N}$ is skew-symmetric. Therefore,

$\displaystyle\langle\nabla\varphi,\mathbf{Q}\bm{\Omega}\rangle_{F}=\langle\mathcal{P}_{T_{\mathbf{Q}}O(n)}(\nabla\varphi),\
\mathbf{Q}\bm{\Omega}\rangle_{F}=\langle\mathcal{P}_{T_{\mathbf{Q}}O(n)}(\nabla\varphi),\
\mathbf{Q}\bm{\Omega}\rangle_{g_{e}}.$

So $\nabla_{g_{e}}\varphi=\mathcal{P}_{T_{\mathbf{Q}}O(n)}(\nabla\varphi)$. ∎

Under this global coordinate, the divergence on $(O(n),g_{e})$ can also be
explicitly computed. The $ij$ entry of $\nabla_{g_{e}}\varphi$ is given by

(6.8)
$\displaystyle(\nabla_{g_{e}}\varphi)_{ij}=\dfrac{1}{2}\left(\dfrac{\partial\varphi}{\partial
q_{ij}}-\sum_{k,l}q_{ik}q_{lj}\dfrac{\partial\varphi}{\partial
q_{lk}}\right).$

So given a tangent vector field $\mathbf{H}(\mathbf{Q})=(h_{ij}(\mathbf{Q}))$,
the divergence of it is defined by

(6.9)
$\displaystyle\nabla_{g_{e}}\cdot\mathbf{H}(\mathbf{Q}):=\sum_{i,j}(\nabla_{g_{e}}(h_{ij}(\mathbf{Q})))_{ij}.$

The Laplace-Beltrami operator is then defined as:

(6.10)
$\displaystyle\Delta_{g_{e}}\varphi:=\nabla_{g_{e}}\cdot\nabla_{g_{e}}\varphi.$

Explicit expression of the Laplace-Beltrami operator is given by the following
lemma:

######  Lemma 6.3.

(Laplace-Beltrami operator) The Laplace-Beltrami operator on $(O(n),g_{e})$ is
given by

(6.11)
$\displaystyle\Delta_{g_{e}}\varphi=\dfrac{1}{2}\left(\sum_{i,j}\dfrac{\partial^{2}\varphi}{\partial{q_{ij}^{2}}}-(n-1)\sum_{i,j}q_{ij}\dfrac{\partial\varphi}{\partial
q_{ij}}-\sum_{i,j,k,l}q_{ik}q_{lj}\dfrac{\partial^{2}\varphi}{\partial
q_{ij}\partial q_{lk}}\right).$

###### Proof.

By the expression of the divergence and gradient, we have

$\displaystyle\Delta_{g_{e}}\varphi=\sum_{i,j}\left(\dfrac{\partial
h_{ij}}{\partial q_{ij}}-\sum_{k,l}q_{ik}q_{lj}\dfrac{\partial
h_{ij}}{\partial q_{lk}}\right),\
h_{ij}=\sum_{i,j}\left(\dfrac{\partial\varphi}{\partial
q_{ij}}-\sum_{k,l}q_{ik}q_{lj}\dfrac{\partial\varphi}{\partial
q_{lk}}\right).$

Thus

$\displaystyle\dfrac{\partial h_{ij}}{\partial q_{ij}}$
$\displaystyle=\dfrac{1}{2}\left(\dfrac{\partial^{2}\varphi}{\partial
q_{ij}^{2}}-\sum_{k,l}\delta_{kl}q_{lj}\dfrac{\partial\varphi}{\partial
q_{lk}}-\sum_{k,l}\delta_{il}q_{ik}\dfrac{\partial\varphi}{\partial
q_{lk}}-\sum_{k,l}q_{ik}q_{lj}\dfrac{\partial^{2}\varphi}{\partial
q_{ij}q_{lk}}\right)$
$\displaystyle=\dfrac{1}{2}\left(\dfrac{\partial^{2}\varphi}{\partial
q_{ij}^{2}}-\sum_{l}q_{lj}\dfrac{\partial\varphi}{\partial
q_{lj}}-\sum_{k}q_{ik}\dfrac{\partial\varphi}{\partial
q_{ik}}-\sum_{k,l}q_{ik}q_{lj}\dfrac{\partial^{2}\varphi}{\partial
q_{ij}q_{lk}}\right),$

and

$\displaystyle\dfrac{\partial h_{ij}}{\partial
q_{lk}}=\dfrac{1}{2}\left(\dfrac{\partial^{2}\varphi}{\partial q_{ij}\partial
q_{lk}}-\sum_{k^{\prime},l^{\prime}}\delta_{il}\delta_{kk^{\prime}}q_{l^{\prime}j}\dfrac{\partial\varphi}{\partial
q_{l^{\prime}k^{\prime}}}-\sum_{k^{\prime},l^{\prime}}\delta_{ll^{\prime}}\delta_{kj}q_{ik^{\prime}}\dfrac{\partial\varphi}{\partial
q_{l^{\prime}k^{\prime}}}-\sum_{k^{\prime},l^{\prime}}q_{ik^{\prime}}q_{l^{\prime}j}\dfrac{\partial^{2}\varphi}{\partial
q_{l^{\prime}k^{\prime}}\partial q_{lk}}\right).$

Substituting the above formulas into the Laplacian operator, we can derive
(6.11). ∎

###  6.2. Proofs of lemmas and omitted calculations

####  6.2.1. Section 2

We will provide main steps of proof of (2.6). Details can be found in [15].
Remember that $\mathbf{Q}(t)$ is now not necessarily be orthogonal, but only
invertible.

First, by Lemma 2.1 in [15], an explicit solution formula of (2.5) is given by

(6.12)
$\displaystyle\mathbf{P}(t)=e^{\mathbf{A}t}\mathbf{P}(0)[\mathbf{I_{n}}+(e^{2\mathbf{A}t}-\mathbf{I_{n}})\mathbf{P}(0)]^{-1}e^{\mathbf{A}t}.$

Moreover, notice that
$\mathbf{Q}^{T}(0)(e^{2\mathbf{A}t}-\mathbf{I_{n}})\mathbf{Q}(0)$ is positive
definite for all $t>0$, so it only adimts non-negative eigenvalues. Remember
that $\mathbf{AB}$ and $\mathbf{BA}$ have same non-zero eigenvalues, so
eigenvalues of
$(e^{2\mathbf{A}t}-\mathbf{I_{n}})\mathbf{P}(0)=(e^{2\mathbf{A}t}-\mathbf{I_{n}})\mathbf{Q}(0)\mathbf{Q}^{T}(0)$
are non-negative. Thus
$\mathbf{I_{n}}+(e^{2\mathbf{A}t}-\mathbf{I_{n}})\mathbf{P}(0)$ is invertible,
and solution of (1.6) exists for all time. Moreover, (6.12) implies that
rank$(\mathbf{P}(t))=n$ for all $t>0$.

Second, we prove that under the assumption that $\mathbf{Q}(0)$ is invertible,
the singular values of $\mathbf{Q}(t)$ are all positive and bounded from below
uniformly for all $t>0$. Let orthogonal $\mathbf{U}$ satisfy

$\displaystyle\mathbf{Q}(0)\mathbf{Q}^{T}(0)=\mathbf{U}\mathbf{L}\mathbf{U}^{T},$

here $\mathbf{L}$ is diagonal. Let
$\mathbf{M}=\mathbf{U}e^{2\mathbf{A}t}\mathbf{U}^{T}$, then by (6.12), direct
calculation yields

$\displaystyle\mathbf{P}(t)=e^{\mathbf{A}t}\mathbf{U}^{T}(\mathbf{L}^{-1}+\mathbf{M}-\mathbf{I_{n}})^{-1}\mathbf{U}e^{\mathbf{A}t}.$

Because $\mathbf{Q}(0)$ is invertible and
$e^{\mathbf{A}t}\succ\mathbf{I_{n}}$, thus both $\mathbf{M-I_{n}}$ and
$\mathbf{L}^{-1}$ are positive definite, so are
$(\mathbf{L}^{-1}+\mathbf{M}-\mathbf{I_{n}})^{-1}$ and $\mathbf{P}(t)$. So
$\mathbf{P}(t)$ is positive definite. Following the steps in Theorem 2.2 of
[15], we can prove that singular values of $\mathbf{Q}(t)$ lower bounded
uniformly for $t>0$. Let $\alpha>0$ be the lower bound.

Finally, let $v(t)=\|\mathbf{I_{n}}-\mathbf{Q}^{T}(t)\mathbf{Q}(t)\|_{F}^{2}$,
direct calculation yields

(6.13)
$\displaystyle\dot{v}(t)=-4\mathrm{tr}((\mathbf{I_{n}}-\mathbf{Q}^{T}(t)\mathbf{Q}(t))\mathbf{Q}^{T}(t)\mathbf{A}\mathbf{Q}(t)(\mathbf{I_{n}}-\mathbf{Q}^{T}(t)\mathbf{Q}(t)))\leq-4\alpha^{2}\lambda_{n}v(t).$

By Gronwall’s inequality, we proved exponential convergence under the
Frobenius norm.

####  6.2.2. Section 4

###### Proof of Lemma 4.1.

Suppose that $\mathbf{P_{1}}$ and $\mathbf{P_{2}}$ are orthogonal matrices
which diagonalizes $\mathbf{M}$ and $\mathbf{N}$ respectively:

$\displaystyle\mathbf{M}=\mathbf{P_{1}}^{T}\bm{\Lambda_{1}}\mathbf{P_{1}},\
\mathbf{N}=\mathbf{P_{2}}^{T}\bm{\Lambda_{2}}\mathbf{P_{2}},$

here $\bm{\Lambda_{1}}$ and $\bm{\Lambda_{2}}$ are diagonal matrices and the
entries on the diagonal line in descending order:

$\displaystyle\bm{\Lambda_{1}}=\mathrm{diag}\\{\lambda_{1,1},\lambda_{1,2},...\lambda_{1,n}\\},\
\lambda_{1,1}\geq\lambda_{1,2}\geq...\geq\lambda_{1,n},$
$\displaystyle\bm{\Lambda_{2}}=\mathrm{diag}\\{\lambda_{2,1},\lambda_{2,2},...\lambda_{2,n}\\},\
\lambda_{2,1}\geq\lambda_{2,2}\geq...\geq\lambda_{2,n}.$

Notice that $\|\mathbf{PM}\|_{F}=\|\mathbf{MP}\|_{F}=\|\mathbf{M}\|$ holds for
any orthogonal matrix $\mathbf{P}$, therefore

$\displaystyle\|\mathbf{M}-\mathbf{N}\|_{F}^{2}$
$\displaystyle=\|\mathbf{P_{1}}^{T}\bm{\Lambda_{1}}\mathbf{P_{1}}-\mathbf{P_{2}}^{T}\bm{\Lambda_{2}}\mathbf{P_{2}}\|_{F}^{2}$
$\displaystyle=\|\mathbf{P_{1}}^{T}(\bm{\Lambda_{1}}\mathbf{P_{1}P_{2}}^{T}-\mathbf{P_{1}P_{2}}^{T}\bm{\Lambda_{2}})\mathbf{P_{2}}\|_{F}^{2}$
$\displaystyle=\|\bm{\Lambda_{1}}\mathbf{P_{1}P_{2}}^{T}-\mathbf{P_{1}P_{2}}^{T}\bm{\Lambda_{2}}\|_{F}^{2}.$

Let $\mathbf{Q}=\mathbf{P_{1}P_{2}}^{T}$, then $\mathbf{Q}$ is also an
orthogonal matrix and

(6.14) $\displaystyle\|\mathbf{M}-\mathbf{N}\|_{F}^{2}=\sum_{1\leq i,j\leq
n}|q_{i,j}|^{2}(\lambda_{1,i}-\lambda_{2,j})^{2}.$

Let $\mathbf{W}\in\mathbb{R}^{n\times n}$ be $w_{i,j}=|q_{i,j}|^{2}$ for
$i,j=1,2,...,n$. Then $\mathbf{W}$ is a doubly stochastic matrix since
$\mathbf{Q}$ is orthogonal. Here a matrix $\mathbf{W}$ is doubly stochastic if
and only if both $\mathbf{W}$ and $\mathbf{W}^{T}$ are transition matrices.

Denote the set of all doubly matrices in $\mathbb{R}^{n\times n}$ as
$\mathcal{D}_{n}$. Then $\mathcal{D}_{n}$ is a compact and convex set in
$\mathbb{R}^{n\times n}$ under the Frobenius norm. Define the following
functional for all doubly stochastic matrix $\mathbf{W}$:

(6.15) $\displaystyle\phi(\mathbf{W}):=\sum_{1\leq i,j\leq
n}w_{i,j}(\lambda_{1,i}-\lambda_{2,j})^{2},$

then $\phi$ is a linear (hence convex) functional w.r.t. to $\mathbf{W}$ and

$\displaystyle\|\mathbf{M}-\mathbf{N}\|_{F}^{2}\geq\min\limits_{\mathbf{W}\in\mathcal{D}_{n}}\phi(\mathbf{W}).$

By the celebrated theorem due to Birkhoff and von Neumann which states that
the convex hull of permutation matrices in $\mathbb{R}^{n\times n}$ is
$\mathcal{D}_{n}$, we know the that minimum of $\phi(\mathbf{W})$ is attained
when $\mathbf{W}$ is a permutation matrix. Thus

(6.16)
$\displaystyle\min\limits_{\mathbf{W}\in\mathcal{D}_{n}}\phi(\mathbf{W})=\min\limits_{\sigma\in
S_{n}}\sum_{i=1}^{n}(\lambda_{1,\sigma(i)}-\lambda_{2,i})^{2}.$

Here $S_{n}$ is the $n$th order permutation group. Remember that
$\lambda_{1,i}$ and $\lambda_{2,i}$ are in descending order, by the
rearrangement inequality,

(6.17)
$\displaystyle\sum_{i=1}^{n}\lambda_{1,\sigma(i)}\lambda_{2,i}\leq\sum_{i=1}^{n}\lambda_{1,i}\lambda_{2,i}.$

Thus

$\displaystyle\|\mathbf{M}-\mathbf{N}\|_{F}^{2}\geq\min\limits_{\mathbf{W}\in\mathcal{D}_{n}}\phi(\mathbf{W})=\sum_{i=1}^{n}(\lambda_{1,i}-\lambda_{2,i})^{2}.$

This proves the inequality. The equality holds if and only if
$\mathbf{W}=\mathbf{I_{n}}$, i.e. $\mathbf{M}$ and $\mathbf{N}$ can be
diagonalized simultaneously and the eigenvalues are paired in order. ∎

####  6.2.3. Section 5

###### Proof of Lemma 4.2.

We prove $(i)$. First, we have

(6.18)
$\displaystyle-\mathcal{E}^{\prime}(\mathbf{Q};\mathbf{A},\mathbf{N})=-2\mathbf{AQN}.$

Then, using $\mathrm{tr}(\mathbf{MN})=\mathrm{tr}(\mathbf{MN}),\
\mathbf{M,N}\in\mathbb{R}^{n\times n}$ and
$\mathrm{tr}(\mathbf{M})=\mathrm{tr}(\mathbf{M}^{T})$, we know that

$\displaystyle\mathrm{tr}(\mathbf{N\dot{Q}}^{T}\mathbf{AQ})=\mathrm{tr}(\mathbf{AQ}\mathbf{N\dot{Q}}^{T})=\mathrm{tr}(\mathbf{\dot{Q}N}\mathbf{Q}^{T}\mathbf{A})=\mathrm{tr}(\mathbf{NQ}^{T}\mathbf{A}\mathbf{\dot{Q}}),$

thus according to (2.8), we have

$\displaystyle\dfrac{\mathrm{d}\mathcal{E}(t)}{\mathrm{d}t}$
$\displaystyle=-\mathrm{tr}(\mathbf{N\dot{Q}}^{T}\mathbf{AQ})-\mathrm{tr}(\mathbf{NQ}^{T}\mathbf{A}\mathbf{\dot{Q}})$
$\displaystyle=-2\mathrm{tr}(\mathbf{AQN}\mathbf{\dot{Q}}^{T})$
$\displaystyle=-2\mathrm{tr}\left(\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{AQNE_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}\mathbf{Q}^{T}-\mathbf{AQNE_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}\mathbf{Q}^{T}\right)$
$\displaystyle=-2\mathrm{tr}\left(\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{j}Q}^{T}\mathbf{AQNE_{k}}\mathbf{Q}^{T}\mathbf{AQ}-\mathbf{E_{k}Q}^{T}\mathbf{AQNE_{j}}\mathbf{Q}^{T}\mathbf{AQ}\right)$
$\displaystyle=-2\mathrm{tr}\left(\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{j}^{2}Q}^{T}\mathbf{AQNE_{k}^{2}}\mathbf{Q}^{T}\mathbf{AQ}-\mathbf{E_{k}^{2}Q}^{T}\mathbf{AQNE_{j}^{2}}\mathbf{Q}^{T}\mathbf{AQ}\right)$
$\displaystyle=-2\mathrm{tr}\left[\sum_{j=1}^{n}\sum_{k=1}^{j-1}\left(\mathbf{E_{j}Q}^{T}\mathbf{AQE_{k}}\right)\mathbf{N}\left(\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}\right)-\left(\mathbf{E_{k}Q}^{T}\mathbf{AQE_{j}}\right)\mathbf{N}\left(\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}\right)\right]$
$\displaystyle=-2\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mu_{k}(\mathbf{q_{k}}(t)\cdot\mathbf{Aq_{j}}(t))^{2}-\mu_{j}(\mathbf{q_{k}}(t)\cdot\mathbf{Aq_{j}}(t))^{2}$
$\displaystyle=-2\sum_{j=1}^{n}\sum_{k=1}^{j-1}(\mu_{k}-\mu_{j})(\mathbf{q_{k}}(t)\cdot\mathbf{Aq_{j}}(t))^{2}\leq
0.$

The equality holds if and only if $\mathbf{q_{k}}\cdot\mathbf{Aq_{j}}=0$ for
all $j\neq k,j,k=1,2,...,n$. Equivalently, $\mathbf{q_{j}},\ j=1,2,...,n$ are
unit eigenvectors of $\mathbf{A}$, i.e., $\mathbf{Q}\in E$. This also proves
$(ii)$

Then we prove $(iii)$. Suppose that $\mathbf{Q}\in O(n)$ satisfies
$\mathbf{F}(\mathbf{Q})=0$, then

$\displaystyle\mathbf{F(Q)}=\mathbf{Q}\left(\sum_{j=1}^{n}\sum_{k=1}^{j-1}\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}-\mathbf{E_{k}}\mathbf{Q}^{T}\mathbf{AQE_{j}}\right)=\mathbf{0}.$

Thus for all $j\neq k,j,k=1,2,...,n$, we have
$\mathbf{E_{j}}\mathbf{Q}^{T}\mathbf{AQE_{k}}=0$, or equivalently
$\mathbf{q_{j}}\cdot\mathbf{Aq_{k}}=0$. Thus $\mathbf{F}(\mathbf{Q})=0$ if and
only if $\mathbf{Q}\in E$.

∎

###### Proof of Lemma 5.1.

From (5.2), we have

(6.19)
$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}(e^{-\mathbf{A}t}\mathbf{Q})+e^{-\mathbf{A}t}\mathbf{Q}\mathbf{T}(\mathbf{Q})=0,$

or column-wisely

(6.20)
$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}e^{-\mathbf{A}t}\mathbf{q_{k}}+2\sum_{j=1}^{k-1}(\mathbf{q_{j}}\cdot\mathbf{Aq_{k}})e^{-\mathbf{A}t}\mathbf{q_{j}}+(\mathbf{q_{k}}\cdot\mathbf{Aq_{k}})e^{-\mathbf{A}t}\mathbf{q_{k}}=0,\
k=1,2,...,n.$

First compute $\mathbf{q_{1}}$. Let $k=1$, we have

(6.21)
$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}e^{-\mathbf{A}t}\mathbf{q_{1}}+(\mathbf{q_{1}}\cdot\mathbf{Aq_{1}})e^{-\mathbf{A}t}\mathbf{q_{1}}=0.$

Denote
$g_{1,1}(t)=e^{-\int_{0}^{t}\mathbf{q_{1}(s)}\cdot\mathbf{Aq_{1}}(s)\mathrm{d}s}$,
then we have

(6.22) $\displaystyle
e^{-\mathbf{A}t}\mathbf{q_{1}}(t)=g_{1,1}(t)\mathbf{q_{1,0}}.$

Here $\mathbf{q_{1,0}}$ is the first column of the initial value
$\mathbf{Q_{0}}$ (see Section 2). By Lemma 1.1, we know that
$\|\mathbf{q_{1}}(t)\|=1$, so

(6.23) $\displaystyle
g_{1,1}(t)=\dfrac{1}{\|e^{\mathbf{A}t}\mathbf{q_{1,0}}\|}.$

This solves $\mathbf{q_{1}}$ and we know that the following claim holds for
$j=1$: there exist $g_{l,j}(t),1\leq l\leq j$ such that

(6.24) $\displaystyle
e^{-\mathbf{A}t}\mathbf{q_{j}}(t)=\sum_{l=1}^{j}\mathbf{q_{l,0}}g_{l,j}(t).$

We will prove that (6.24) holds for all $j=1,2,...,n$. Define

(6.25) $\displaystyle
f_{k}(t):=e^{\int_{0}^{t}\mathbf{q_{k}}(s)\cdot\mathbf{Aq_{k}}(s)\mathrm{d}s},\
k=1,2,...,n.$

Suppose that (6.24) holds for $j=k-1,k\geq 2$, then by (6.20) and (6.24), we
have

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}\left(f_{k}(t)e^{-\mathbf{A}t}\right)$
$\displaystyle=-2\sum_{j=1}^{k-1}f_{k}(t)(\mathbf{q_{j}}\cdot\mathbf{Aq_{k}})e^{-\mathbf{A}t}\mathbf{q_{j}}$
$\displaystyle=-2f_{k}(t)\sum_{l=1}^{k-1}\left[\sum_{j=l}^{k-1}g_{l,j}(t)(\mathbf{q_{j}}\cdot\mathbf{Aq_{k}})\right]\mathbf{q_{l,0}}.$

Integrating on both sides yields

(6.26) $\displaystyle
e^{-\mathbf{A}t}\mathbf{q_{k}}(t)=\dfrac{1}{f_{k}(t)}\mathbf{q_{k,0}}-\sum_{l=1}^{k-1}\left[\int_{0}^{t}\dfrac{2f_{k}(s)}{f_{k}(t)}\sum_{j=l}^{k-1}g_{j,l}(s)(\mathbf{q_{j}}(s)\cdot\mathbf{Aq_{k}}(s))\mathrm{d}s\right]\mathbf{q_{l,0}}.$

Thus (6.24) also holds for $j=k$ and we derive a iteration formula for
$g_{l,j},\ 1\leq l\leq j$:

(6.27) $\displaystyle g_{l,k}(t)=\left\\{\begin{split}&\dfrac{1}{f_{k}(t)},\ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ l=k,\\\
&-2\int_{0}^{t}\dfrac{f_{k}(s)}{f_{k}(t)}\sum_{j=l}^{k-1}g_{j,l}(s)(\mathbf{q_{j}}(s)\cdot\mathbf{Aq_{k}}(s))\mathrm{d}s,\
1\leq l\leq k-1.\end{split}\right.$

In fact, $g_{l,k}(t),1\leq l\leq k\leq n$ are exactly elements of
$\mathbf{G}(t)$. Thus (5.5) holds, and $g_{k,k}(t)=1/f_{k}(t)>0$, i.e.,
elements of $\mathbf{G}(t)$ on the diagonal line are all positive.

Finally, by Lemma 1.1, we have

(6.28)
$\displaystyle\mathbf{Q}(t)\mathbf{Q}^{T}(t)=e^{\mathbf{A}t}\mathbf{Q_{0}G}(t)\mathbf{G}^{T}(t)\mathbf{Q_{0}}^{T}e^{\mathbf{A}t}=\mathbf{I_{n}}.$

Thus

$\displaystyle\mathbf{G}(t)\mathbf{G}^{T}(t)=\mathbf{Q_{0}}^{T}e^{-2\mathbf{A}t}\mathbf{Q_{0}}.$

So (5.6) holds. ∎

###### Proof of Lemma 5.2.

We prove by induction. Apparently (5.12) holds for $m=1$ by definition.
Suppose that (5.12) holds for $m$, we prove that (5.12) also holds for $m+1$.
Define

(6.29) $\displaystyle I_{m,j}:=\\{\sigma_{k}(\mathbf{Q}):\
\sigma_{k}(\mathbf{Q})\leq j,\ 1\leq k\leq m\\},\ 1\leq m,j\leq n.$

Notice that

$\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m+1;\
1,2,...,j]\right)-\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right)=0\ \mathrm{or}\ 1,$

so we discuss the following two cases:

  1. (1)

If $j\geq\sigma_{m+1}(\mathbf{Q})$, we just need to prove that

$\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m+1;\
1,2,...,j]\right)-\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right)=1.$

By definition of $\sigma$, we know that column vectors of
$\mathbf{Q}[1,2,...,m;\ I_{m,j}]$ are linearly independent. Notice that by
induction hypothesis,

$\displaystyle\ \ \ \ \mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right)$ $\displaystyle=m-\mathrm{card}(\\{k:\
\sigma_{k}(\mathbf{Q})>j,\ 1\leq k\leq m\\})$
$\displaystyle=\mathrm{card}(I_{m,j}).$

Thus column vectors of $\mathbf{Q}[1,2,...,m;\ I_{m,j}]$ form a basis of the
column space of $\mathbf{Q}[1,2,...,m;\ 1,2,...,j]$. However, because
$j\geq\sigma_{m+1}(\mathbf{Q})$, thus $\mathbf{Q}[1,2,...,m+1;\
\sigma_{m+1}(\mathbf{Q})]$ is linearly independent with column vectors of
$\mathbf{Q}[1,2,...,m+1;\ I_{m,j}]$, thus

$\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m+1;\
1,2,...,j]\right)-\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right)=1.$

  2. (2)

If $j<\sigma_{m+1}(\mathbf{Q})$, then

$\displaystyle\mathrm{card}(\\{k:\ \sigma_{k}(\mathbf{Q})>j,\ 1\leq k\leq
m+1\\})-\mathrm{card}(\\{k:\ \sigma_{k}(\mathbf{Q})>j,\ 1\leq k\leq m\\})=1.$

Consider integer $j^{\prime}$ such that $1\leq j^{\prime}\leq j$. Because
$j^{\prime}<\sigma_{m+1}(\mathbf{Q})$, so $\mathbf{Q}[1,2,...,m+1;\
j^{\prime}]$ can be linearly represented by column vectors of
$\mathbf{Q}[1,2,...,m+1;\ \sigma_{1}(\mathbf{Q}),\ \sigma_{2}(\mathbf{Q}),\
...,\ \sigma_{m}(\mathbf{Q})]$. Suppose that

(6.30) $\displaystyle\mathbf{Q}[1,2,...,m+1;\
j^{\prime}]=\sum_{k=1}^{m}y_{k}\mathbf{Q}[1,2,...,m+1;\
\sigma_{k}(\mathbf{Q})].$

If focusing on the first $m$ rows, we have

$\displaystyle\mathbf{Q}[1,2,...,m;\
j^{\prime}]=\sum_{k=1}^{m}y_{k}\mathbf{Q}[1,2,...,m;\
\sigma_{k}(\mathbf{Q})].$

By definition of $\sigma$, we know that for any $1\leq j^{\prime}\leq j$, the
column vector $\mathbf{Q}[1,2,...,m;\ j^{\prime}]$ can be linearly represented
by column vectors of $\mathbf{Q}[1,2,...,m;\ I_{m,j}]$ (see the argument in
(i)). Uniqueness of linear representations determines that

(6.31) $\displaystyle y_{k}=0\ \mathrm{if}\ \sigma_{k}(\mathbf{Q})\notin
I_{m,j},\ k=1,2,...,m.$

Back to (6.30), we know that $\mathbf{Q}[1,2,...,m+1;\ j^{\prime}]$ can be
linearly represented by $\mathbf{Q}[1,2,...,m+1;\ I_{m,j}]$. Thus

$\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m+1;\
1,2,...,j]\right)=\mathrm{card}(I_{m,j})=\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m;\
1,2,...,j]\right).$

So

$\displaystyle\mathrm{rank}\left(\mathbf{Q}[1,2,3,...,m+1;\
1,2,...,j]\right)+\mathrm{card}(\\{k:\ \sigma_{k}(\mathbf{Q})>j,\ 1\leq k\leq
m+1\\})=m+1.$

∎

###### Proof of Lemma 5.3.

If $m=1$, by definition $q_{1,j}=0$ for all $1\leq j<\sigma_{1}(\mathbf{Q})$.
So (5.14) holds.

If $2\leq m\leq n$, by definition, $\mathbf{Q}[1,2,3,...,m-1;\
\mathbf{i_{m-1}}]\in\mathbb{R}^{(m-1)\times(m-1)}$ is invertible, so
$\mathbf{c_{m}}$ is uniquely determined. Moreover, because
$\mathbf{Q}[1,2,...,m;\ \sigma_{m}(\mathbf{Q})]$ is linearly independent with
$\mathbf{Q}[1,2,...,m;\ \sigma_{k}(\mathbf{Q})],1\leq k\leq m-1$, so there
exists $z_{m}$ such that

(6.32) $\displaystyle\mathbf{Q}[m;\
\mathbf{i_{m}}]=\mathbf{c_{m}}\mathbf{Q}[1,2,...,m-1;\
\mathbf{i_{m-1}}]+(0,0,...,0,z_{m},0,...,0).$

Taking determinant on both sides yields

(6.33) $\displaystyle z_{m}=\dfrac{\det(\mathbf{Q}[1,2,3,...,m;\
\mathbf{i_{m}}])}{\det(\mathbf{Q}[1,2,3,...,m-1;\ \mathbf{i_{m-1}}])}.$

The only thing left to prove is that

(6.34) $\displaystyle\mathbf{Q}[m;\
1,2,3,...,\sigma_{m}(\mathbf{Q})]-\mathbf{c_{m}}\mathbf{Q}[1,2,...,m-1;\
1,2,3,...,\sigma_{m}(\mathbf{Q})]=(0,0,...,0,z_{m}).$

If $\sigma_{m}(\mathbf{Q})=1$, then $\mathbf{c_{m}}=\mathbf{0}$ and
$z_{m}=q_{m,1}$, (6.34) holds. If $\sigma_{m}(\mathbf{Q})\geq 1$, denote
$K_{m}=\\{k:\sigma_{k}(\mathbf{Q})\leq\sigma_{m}(\mathbf{Q})-1,1\leq k\leq
m\\}$. By (6.32), all entries on $k$-th column where $k\in K_{m}$ are zero. So
we just need to consider those entries are not on these columns.

Otherwise, suppose that there exists $j\notin K_{m},1\leq
j<\sigma_{m}(\mathbf{Q})$ such that the entry on the $j$-th column is non-
zero. Denote $\mathbf{B}=\mathbf{C_{m}}\mathbf{Q}[1,2...,m;\
1,2,3,...,\sigma_{m}(\mathbf{Q})-1]$ where

(6.35) $\displaystyle\mathbf{C_{m}}=\begin{pmatrix}1&0&0&...&0\\\
0&1&0&...&0\\\ 0&0&1&...&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
-c_{m,1}&-c_{m,2}&-c_{m,3}&...&1\end{pmatrix}$

is in $\mathbb{R}^{m\times m}.$ Then the L.H.S. of (6.34) is exactly the
$m$-th row of $\mathbf{B}$. Meanwhile, $\mathbf{B}$ has same rank with
$\mathbf{Q}[1,2,...,m;\ 1,2,3,...,\sigma_{m}(\mathbf{Q})-1]$ since
$\mathbf{C_{m}}$ is invertible. By Lemma 5.2, we know that

$\displaystyle\mathrm{rank}(\mathbf{B})$
$\displaystyle=\mathrm{rank}(\mathbf{Q}[1,2,...,m;\
1,2,3,...,\sigma_{m}(\mathbf{Q})-1])$
$\displaystyle=m-\mathrm{card}(\\{k:\sigma_{k}(\mathbf{Q})>\sigma_{m}(\mathbf{Q})-1,1\leq
k\leq m\\})$
$\displaystyle=\mathrm{card}(\\{k:\sigma_{k}(\mathbf{Q})\leq\sigma_{m}(\mathbf{Q})-1,1\leq
k\leq m\\})$ $\displaystyle=\mathrm{card}(K_{m}).$

By definition of $\sigma$, we know that $\mathbf{B}[1,2,...,m;\
\sigma_{k}(\mathbf{Q})]$ are linearly independent if $k\in K$. Meanwhile, for
all $k\in K_{m}$, by (6.32), $\mathbf{B}[m;\
\sigma_{k}(\mathbf{Q})]=\mathbf{0}$. Remember that we assume $b_{m,j}$ is non-
zero for some $j<\sigma_{m}(\mathbf{Q})$, so $\mathbf{B}[1,2,3,...,m;j]$ is
linearly independent with $\mathbf{B}[1,2,...,m;\ \sigma_{k}(\mathbf{Q})],\
k\in K_{m}$, hence

(6.36) $\displaystyle\mathrm{card}(K_{m})=\mathrm{rank}(\mathbf{B})\geq
1+\mathrm{card}(K_{m})$

which is a contradiction. So (6.34) holds. ∎

######  Remark 6.1.

The matrix $\mathbf{C_{m}}$ can be interpreted as row transformations on
$\mathbf{Q}[1,2,...,m;\ 1,2,...,\sigma_{m}(\mathbf{Q})-1]$ to clear all
entries on the $m$-th row, the $k$-th column, $k\in K_{m}$. The aim of the
last part of the proof is to prove that not only these columns, but also all
columns are cleared.

###  6.3. Linearization near stable points

We will prove that the asymptotic stable points of (1.5) are $\mathbf{Q}^{*}$
defined in (5.26). Without loss of generality, we just consider the case where
$\epsilon_{i}=1,\ i=1,2...,n$ in (1.5).

Consider equilibrium
$\mathbf{Q_{1}}=[e_{\tau(1)},e_{\tau(2)},...,e_{\tau(n)}]$, here
$\tau(i),i=1,2,...,n$ is a permutation of $1,2,...,n$. We will prove that if
$\mathbf{Q_{1}}$ is table, then $\tau(i)=i$ holds for all $i=1,2,...,n$. We
first consider linearization of $\mathbf{q}_{1}$ which reads as

$\displaystyle\dot{\mathbf{q_{1}}}$
$\displaystyle=\mathbf{A}(\mathbf{q_{1}}-\mathbf{e}_{\tau(1)})-[(\mathbf{e}_{\tau(1)}^{T}\mathbf{A}\mathbf{e}_{\tau(1)})\mathbf{I_{n}}+2\lambda_{\tau(1)}\mathbf{e}_{\tau(1)}\mathbf{e}_{\tau(1)}^{T}](\mathbf{q_{1}}-\mathbf{e}_{\tau(1)})$
$\displaystyle=(\mathbf{A}-\lambda_{\tau(1)}\mathbf{I_{n}}-2\lambda_{\tau(1)}\mathbf{E}_{\tau(1)})(\mathbf{q_{1}}-\mathbf{e}_{\tau(1)}).$

Remember that $\mathbf{A}$ is assumed to be diagonal with entries on diagonal
line aligned in a descending order, so
$\mathbf{B_{1}}=\mathbf{A}-\lambda_{\tau(1)}\mathbf{I_{n}}-2\lambda_{\tau(1)}\mathbf{E}_{\tau(1)}$
is also diagonal. Let
$\mathbf{B_{1}}=\mathrm{diag}\\{b_{1,1},b_{1,2},...,b_{1,n}\\}$, then

(6.37) $\displaystyle
b_{j}=\left\\{\begin{split}&\lambda_{j}-\lambda_{\tau(1)},\ \ \ \ \
j\neq\tau(1),\\\ &-2\lambda_{\tau(1)},\ \ \ \ \ \
j=\tau(1).\end{split}\right.$

Because this equilibrium is stable, we have $\lambda_{j}<\lambda_{\tau(1)}$
for all $j\neq\tau(1)$. Thus $\lambda_{\tau(1)}$ is the largest one in
$\lambda_{i},\ i=1,2,...,n$. Thus $\tau(1)=1.$

Then, we prove $\tau(j)=j$ for $j\geq 2$ by induction. Suppose that
$\tau(i)=i$ holds for all $1\leq i\leq j-1$. Therefore, linearization of
$\mathbf{q_{j}}$ reads as

(6.38)
$\displaystyle\dot{\mathbf{q_{j}}}=\left(\mathbf{A}-\lambda_{\tau(j)}\mathbf{I_{n}}-2\lambda_{\tau(j)}\mathbf{E}_{\tau(j)}-2\sum_{i=1}^{j-1}\lambda_{i}\mathbf{E_{i}}\right)(\mathbf{q_{j}}-\mathbf{e}_{\tau(j)})-2\left(\sum_{i=1}^{j-1}\lambda_{\tau(j)}\mathbf{E}_{\tau(j),i}(\mathbf{q_{i}}-\mathbf{e_{i}})\right).$

Again,
$\mathbf{B_{j}}=\mathbf{A}-\lambda_{\tau(j)}\mathbf{I_{n}}-2\lambda_{\tau(j)}\mathbf{E}_{\tau(j)}-2\sum_{i=1}^{j-1}\lambda_{i}\mathbf{E_{i}}$
is a diagonal matrix. Let $\mathbf{B_{j}}=\\{b_{j,1},b_{j,2},...,b_{j,n}\\}$,
then

(6.39) $\displaystyle
b_{j,i}=\left\\{\begin{split}&-\lambda_{\tau(j)}-\lambda_{i},\ i<j,\\\
&\lambda_{i}-\lambda_{\tau(j)},\ \ \ \ \ i\geq j,i\neq\tau(j),\\\
&-2\lambda_{\tau(j)},\ \ \ \ \ \ i=\tau(j).\end{split}\right.$

Again, $b_{j,i},i=1,2,...,n$ are negative due to asymptotic stability. Thus
$\lambda_{i}<\lambda_{\tau(j)}$ holds for all $i\geq j,i\neq\tau(j)$. Thus
$\tau(j)=j$.

Therefore, by induction, we proved that $\tau(i)=1,i=1,2,...,n.$ So
$\mathbf{Q_{1}}$ is asymptotic stable if and only if it is in the set defined
in (5.26).

######  Remark 6.2.

In (6.38), we see that residues $\mathbf{q_{i}}-\mathbf{e_{i}}.i=1,2,...,j-1$
also contribute to the evolution of $\mathbf{q_{j}}$. This indicates that the
convergence rate of $\mathbf{q_{j}}$ depends on those of
$\mathbf{q_{i}},i=1,2,...,j-1$ and can not exceed them. This is clearly
observed in Theorem 5.3