# Scalable symmetric Tucker tensor decomposition

Ruhui Jin, Joe Kileel,11footnotemark: 1 Tamara G. Kolda, and Rachel
Ward11footnotemark: 1 Department of Mathematics, University of Texas, Austin,
TXMathSci.ai, Dublin, CA

###### Abstract

We study the best low-rank Tucker decomposition of symmetric tensors,
advocating a straightforward projected gradient descent (PGD) method for its
computation. The main application of interest is in decomposing higher-order
multivariate moments, which are symmetric tensors. We develop scalable
adaptations of the basic PGD method and higher-order eigenvalue decomposition
(HOEVD) to decompose sample moment tensors. With the help of implicit and
streaming techniques, we evade the overhead cost of building and storing the
moment tensor. Such reductions make computing the Tucker decomposition
realizable for large data instances in high dimensions. Numerical experiments
demonstrate the efficiency of the algorithms and the applicability of moment
tensor decompositions to real-world datasets. Last, we study the convergence
on the Grassmannian manifold, and prove that the update sequence derived by
the PGD solver achieves first and second-order criticality.

##  1 Introduction

Tensor decompositions [35] give powerful tools to represent and compress
higher-order arrays. In this work, we focus on the _symmetric rank- $r$ Tucker
decomposition_ [61]. Applied to a $d$-th order symmetric tensor of size
$n\times\ldots\times n$, it expresses the input using a symmetric core tensor
$\boldsymbol{\mathscr{C}}\in\mathbb{R}^{r\times\dots\times r}$ $(r\ll n)$ and
a matrix $\mathbf{Q}\in\mathbb{R}^{n\times r}$ with orthonormal columns. See
Figure 1 for an illustration.

$\boldsymbol{\mathscr{X}}$ $\approx$ $\boldsymbol{\mathscr{C}}$ $\mathbf{Q}$
$\mathbf{Q}$ $\mathbf{Q}$ Figure 1: Symmetric Tucker decomposition.

_Symmetric tensors_ abound in statistics as a means of representing higher-
order moments and cumulants of a random vector. They are commonly used to
analyze dependences in multivariate distributions [12, 29, 13], learn latent
variable models via the method-of-moments [4, 23] and solve important signal
processing problems [9, 10], among many others.

In most cases, the exact moment of the ground-truth distribution is not
available. The idea of moment estimation is to use the $d$-th order _sample
moment tensor_ $\boldsymbol{\mathscr{M}}^{d}$ instead. This is formed from
data observations
$\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{p}\\}\subseteq\mathbb{R}^{n}$ of a
random variable $\mathbf{x}$ as

$\boldsymbol{\mathscr{M}}^{d}=\frac{1}{p}\sum_{i=1}^{p}\mathbf{x}_{i}^{\otimes
d}\approx\mathbb{E}(\mathbf{x}^{\otimes d}).$ (1)

There is an urgent need for economical representations of
$\boldsymbol{\mathscr{M}}^{d}$, which otherwise takes $\mathcal{O}(n^{d})$
space to store. Using the symmetric Tucker decomposition, one approximates
$\boldsymbol{\mathscr{M}}^{d}$ solely by determining the factor matrix
$\mathbf{Q}$, which takes just $rn$ memory space. (The core
$\boldsymbol{\mathscr{C}}$ is a byproduct once $\mathbf{Q}$ is identified, see
Remark 4.) However, decomposing the moment still has upfront costs of
$\mathcal{O}(pn^{d})$ and $\mathcal{O}(n^{d})$ to initially construct and
store $\boldsymbol{\mathscr{M}}^{d}$. For large data volumes $p$ and high
dimensionalities $n$, we turn to _implicit_ (or _tensor-free_) and _streaming_
strategies to escape such a prohibitive overhead.

The symmetric Tucker decomposition of moment tensors has other motivations,
besides compression. In some scientific domains, one encounters data having
high volatility, e.g., in combustion simulations [2], outlier detection [24,
53] and financial asset allocation [47, 8, 6]. Moments and cumulants ($d\geq
3$) are sensitive to heavy tails and anomalies. Thus, as principal component
analysis (PCA) and the covariance matrix are well-suited to Gaussian data,
finding the key subspaces of variability in higher-order statistics is
necessary for non-normal data analysis. These subspaces are found by computing
the low-rank symmetric Tucker approximation.

In other situations, high-dimensional observations are assumed to lie in a
low-dimensional subspace. For example, this holds in the Fama-French model
[21] of finance and in the psychometrics [46] discipline. The resulting moment
tensors admit a symmetric Tucker format. Thus the decomposition enforces the
low-rank structure and improves moment estimation accuracy.

###  1.1 Related works

Due to the NP-hard nature of tensor decompositions [28], an approximate
symmetric Tucker decomposition is often computed by higher-order eigenvector
decomposition (HOEVD) [14]. It is the eigendecomposition of the flattened
tensor and can be viewed as the symmetric version of higher-order singular
value decomposition (HOSVD) [13] and sequentially truncated HOSVD [62].
Despite their popularity, HOSVD-based approaches offer no guarantee of
returning critical points for the problem. Recently, randomized projection
methods have been applied to accelerate the procedure, see [42, 59, 44, 11, 3]
for state-of-the-art implementations of HOSVD. But these methods may not be a
good fit for the symmetric setting. Also, they take the full tensor as input.

We instead turn to iterative search methods to try to achieve the best
possible solution, using HOEVD as an initialization. For Tucker decomposition,
Newton [20, 32], quasi-Newton [55] and trust-region [31] methods have been
considered. These second-order optimization approaches rely on sophisticated
updates on manifolds, based on Riemannian gradient and Hessian computations.
Higher-order orthogonal iterations (HOOI) [15], Jacobi rotations [30, 40] and
the power method [52] have been proposed for symmetric tensors. Among these,
the power method does not always converge, as shown in numerical examples of
[52, Section 6]. Regalia briefly mentions a shifted variant of the power
method in [52], which is akin to a gradient descent method, although strong
convergence analysis is omitted in the work.

Regarding the decomposition of sample moment tensors 1, streaming and implicit
techniques can be applied to address the computational burden in the data-rich
and high-dimensional regime. For covariance matrices which are centered
second-order moments, the streaming PCA scheme a.k.a. Oja’s method [49], and
its adaptations [45, 65, 27], have been thoroughly studied. Recently, a
similar strategy is applied to optimization with higher-order moment tensors
for the task of learning Gaussian mixture models, see the works [57, 50]. The
common ideas with the present work are that we pass over only a subset of
samples at a time and drive the optimization without forming gigantic moment
tensors.

###  1.2 Our contributions

In this paper, we apply the projected gradient descent (PGD) framework to the
symmetric Tucker tensor decomposition problem, see Algorithm 1. The approach
is led by a straightforward first-order gradient update governed by the Tucker
approximation cost and it is free of non-trivial higher-order update methods
on manifolds. We also establish convergence guarantees. Employing theory from
manifold optimization and dynamical systems, we prove that for any symmetric
tensor, Algorithm 1 converges monotonically to a first-order critical point
given arbitrary initialization. Almost surely, it converges to a second-order
critical point.

In Section 3, our attention focuses to decomposing sample moment tensors 1. By
exploiting the special input structure, we modify the PGD framework and
develop the much more efficient Algorithm 2. There are two main highlights.

Firstly, we adopt the implicit technique, so that there is no need to build
the moment tensor for the optimization procedure. We operate on the input
sample vectors directly, which are of much smaller size. Proposed in [57, 50],
the implicit strategy greatly reduces the memory cost and effort for gradient
computation. The second advantage is that Algorithm 2 is built on a streaming
model. By passing over small batches of data observations, this leaves no
storage requirement for huge numbers of samples. Similar ideas have been used
in CP moment tensor decomposition [57].

We also develop a scalable implementation of HOEVD for sample moment tensors
1, see Algorithm 3. This is akin to the streaming PCA idea. While Algorithm 3
may be of independent interest, we use the HOEVD approximation as the
initialization to Algorithm 2.

The savings from the streaming and implicit strategies allow Algorithms 2 and
3 to handle large-scale computations. The efficiency gains are summarized in
Table 1. Here $p$ is the total sample number and $b$ is the batch size at each
iteration. Usually, $p$ is hundreds to thousands times of $b$. These
advantages are corroborated by extensive numerical experiments in Section 4.
Moment tensor decomposition is also applied to real-world datasets for anomaly
detection and portfolio allocation.

Table 1: Computation and storage of Algorithms 2 and 3 | overhead | storage per iter.
---|---|---
streaming $\&$ implicit | none | $nr+nb$
full $\&$ implicit | none | $nr+np$
full $\&$ explicit | $n^{d}p$ | $nr+n^{d}$

##  2 Symmetric Tucker tensor decomposition

###  2.1 Notation and background

Scalars, vectors, matrices and tensors are represented as lowercase letters
$x$, lowercase bold letters $\mathbf{x}$, uppercase bold letters $\mathbf{X},$
and calligraphic bold letters $\boldsymbol{\mathscr{X}}$, respectively. The
$(i_{1},\dots,i_{d})$-th element of a $d$-way tensor is denoted by
$x_{i_{1},\dots,i_{d}}$. The $j$-th column of a matrix $\mathbf{X}$ is
$\mathbf{x}_{j}$. The spectral norm of a matrix is $\|\cdot\|_{2}$ . The
identity matrix is $\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$. The Moore-
Penrose pseudoinverse of a matrix is denoted using ${\dagger}$. The matrix
symmetrization operator is written $\operatorname{sym}:\mathbb{R}^{n\times
n}\to\mathbb{S}^{2}(\mathbb{R}^{n})$ and given by
$\operatorname{sym}(\mathbf{X})=\tfrac{1}{2}(\mathbf{X}+\mathbf{X}^{\top})$
for $\mathbf{X}\in\mathbb{R}^{n\times n}$. We write QR for the function that
selects the orthonormal factor $\mathbf{Q}\in\mathbb{R}^{n\times r}$ from the
QR decomposition of a full-rank $n\times r$ matrix.111For ease of analysis, we
ensure uniqueness of the QR factorization by using the one where all diagonal
entries of the $\mathbf{R}$ factor are positive. But none of our algorithms
require this choice. The notation $\mathbf{X}^{[d]}$ is the elementwise
exponentiation, where each coordinate of $\mathbf{X}$ is raised to the $d$-th
power. The set $\\{1,\ldots,n\\}$ is denoted by $[n]$ when $n$ is a positive
integer.

######  Definition 1 (symmetric tensors).

A tensor $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{n\times\dots\times n}$ is
symmetric if its entries are invariant under any permutation of indices, i.e.

$x_{i_{1},\dots,i_{d}}=x_{i_{\sigma(1)},\dots,i_{\sigma(d)}},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall(i_{1},\dots,i_{d})\in[n]^{d}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \sigma\in\Pi(d)$

where $\Pi(d)$ is the permutation group on $[d]$. Here $d$ and $n$ are called
the order and dimension of $\boldsymbol{\mathscr{X}}$ respectively. The space
of all real-valued symmetric tensors of order $d$ and dimension $n$ is denoted
by $\mathbb{S}^{d}(\mathbb{R}^{n})$.

######  Definition 2 (symmetric outer product).

The outer product of a vector $\mathbf{x}\in\mathbb{R}^{n}$ with itself $d$
times is denoted as $\mathbf{x}^{\otimes d}\in\mathbb{S}^{d}(\mathbb{R}^{n})$.
The $(i_{1},\dots,i_{d})$-th entry of this product is

$\left(\mathbf{x}^{\otimes d}\right)_{i_{1},\dots,i_{d}}=x_{i_{1}}\\!\cdots
x_{i_{d}},\leavevmode\nobreak\ \leavevmode\nobreak\
\forall(i_{1},\dots,i_{d})\in[n]^{d}.$

######  Definition 3 (symmetric Tucker product).

For a symmetric tensor
$\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ and a matrix
$\mathbf{Y}\in\mathbb{R}^{n\times m}$, their symmetric Tucker product is
defined as $\boldsymbol{\mathscr{X}}$ with $\mathbf{Y}$ multiplied along each
mode,
$\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Y},\dots,\mathbf{Y}\right)\in\mathbb{S}^{d}(\mathbb{R}^{m}).$
The $(j_{1},\dots,j_{d})$-th entry of this product is

$\left(\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Y},\dots,\mathbf{Y}\right)\right)_{i_{1},\dots,i_{d}}=\sum_{j_{1}=1}^{n}\cdots\sum_{j_{d}=1}^{n}x_{j_{1},\dots,j_{d}}\,y_{j_{1},i_{1}}\dots
y_{j_{d},i_{d}},\leavevmode\nobreak\ \leavevmode\nobreak\
\forall(i_{1},\dots,i_{d})\in[m]^{d}.$

######  Definition 4 (symmetric tensor matricization).

A symmetric tensor $\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})$
can be flattened into an as wide-as-possible matrix, denoted
$\operatorname{mat}(\boldsymbol{\mathscr{X}})\in\mathbb{R}^{n\times n^{d-1}}$.
The index $(i_{1},\dots,i_{d})\in[n]^{d}$ of $\boldsymbol{\mathscr{X}}$ is
mapped to the index
$\left(i_{1},1+\sum_{\ell=2}^{d}(i_{\ell}-2)\,n^{\ell-1}\right)\in[n]\times[n^{d-1}]$
of $\operatorname{mat}(\boldsymbol{\mathscr{X}})$.

######  Definition 5 (inner product).

For (not necessarily symmetric) tensors
$\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{Y}}\in\mathbb{R}^{n_{1}\times
n_{2}\dots\times n_{d}}$ of the same size, their inner product is

$\langle\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{Y}}\rangle=\sum_{i_{d}=1}^{n_{d}}\dots\sum_{i_{1}=1}^{n_{1}}x_{i_{1},\dots,i_{d}}y_{i_{1},\dots,i_{d}}\in\mathbb{R}.$

The norm of a tensor is
$\left\|\boldsymbol{\mathscr{X}}\right\|=\sqrt{\langle\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{X}}\rangle}.$

######  Definition 6.

For tensors $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{m_{1}\times
n_{2}\dots\times n_{d}}$ and
$\boldsymbol{\mathscr{Y}}\in\mathbb{R}^{m^{\prime}_{1}\times n_{2}\dots\times
n_{d}},$ their inner product in all modes except the first one [20, 16] is
defined as:

$\mathbf{Z}=\langle\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{Y}}\rangle_{-1}\in\mathbb{R}^{m_{1}\times
m^{\prime}_{1}},\leavevmode\nobreak\ \leavevmode\nobreak\
\text{where}\leavevmode\nobreak\
z_{i_{1},i^{\prime}_{1}}=\sum_{j_{d}=1}^{n_{d}}\dots\sum_{j_{2}=1}^{n_{2}}x_{i_{1},j_{2},\dots,j_{d}}y_{i^{\prime}_{1},j_{2},\dots,j_{d}},\leavevmode\nobreak\
\forall\leavevmode\nobreak\
i_{1}\in[m_{1}],i^{\prime}_{1}\in[m^{\prime}_{1}].$

In parts of the paper, we use the standard big $\mathcal{O}$ notation, see
e.g. [1, Appendix A.4].

###  2.2 Problem statement

Let $\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ be a symmetric
tensor. Fix a user-specified rank $r$ (typically $r\ll n$). Our goal is to
find the best approximation of $\boldsymbol{\mathscr{X}}$ represented by the
symmetric Tucker product of a core tensor
$\boldsymbol{\mathscr{C}}\in\mathbb{S}^{d}(\mathbb{R}^{r})$ and a rank-$r$
orthonormal matrix $\mathbf{Q}\in\mathbb{R}^{n\times r}$. The problem is
formulated as follows:

$\min_{\boldsymbol{\mathscr{C}}\in\mathbb{S}^{d}(\mathbb{R}^{r})\atop\mathbf{Q}\in\mathbb{R}^{n\times
r}}\left\|\boldsymbol{\mathscr{X}}-\boldsymbol{\mathscr{C}}\cdot\left(\mathbf{Q}^{\top},\dots,\mathbf{Q}^{\top}\right)\right\|^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \text{subject\leavevmode\nobreak\ to}\leavevmode\nobreak\
\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}_{r}.$ (2)

In fact, the minimization 2 is equivalent to solving for a rank-$r$ basis
$\mathbf{Q}$ that maximizes the following cost function:

$\max_{\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}_{r}}F(\mathbf{Q})\equiv\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\|^{2}.$
($\star$)

######  Remark 1.

The problem formulation transforms from 2 to $\star$ ‣ 2.2 when we eliminate
$\boldsymbol{\mathscr{C}}$ in 2. Solving the unconstrained linear least
squares in 2 for the core, we obtain

$\boldsymbol{\mathscr{C}}=\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right).$
(3)

Substituting 3 into 2 gives

$\begin{array}[]{ll}\left\|\boldsymbol{\mathscr{X}}-\boldsymbol{\mathscr{C}}\cdot\left(\mathbf{Q}^{\top},\dots,\mathbf{Q}^{\top}\right)\right\|^{2}&=\left\|\boldsymbol{\mathscr{X}}-\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q}\mathbf{Q}^{\top},\dots,\mathbf{Q}\mathbf{Q}^{\top}\right)\right\|^{2}\\\
&=\left\|\boldsymbol{\mathscr{X}}\right\|^{2}-\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q}\mathbf{Q}^{\top},\dots,\mathbf{Q}\mathbf{Q}^{\top}\right)\right\|^{2}\\\
&=\left\|\boldsymbol{\mathscr{X}}\right\|^{2}-\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\|^{2}.\end{array}$

As $\|\boldsymbol{\mathscr{X}}\|^{2}$ is a constant, minimizing the above
quantity is the same as maximizing
$\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\|^{2}$.

The feasible domain in $\star$ ‣ 2.2 for $\mathbf{Q}$ is the Stiefel manifold:

$\operatorname{St}(n,r):=\left\\{\mathbf{Q}\in\mathbb{R}^{n\times
r}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}_{r}\right\\},$ (4)

which is a compact Riemannian submanifold of $\mathbb{R}^{n\times r}$ with
dimension $\frac{r(2n-r-1)}{2}$. However the cost function $F$ is rotationally
invariant, in fact depending on the orthogonal projector
$\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}\in\mathbb{S}^{2}(\mathbb{R}^{n})$,
since
$\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\|=\left\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{P},\dots,\mathbf{P}\right)\right\|$.
Therefore, the optimization $\star$ ‣ 2.2 naturally runs on the Grassmannian
manifold, which consists of all orthogonal projectors $\mathbf{P}$:

$\operatorname{Gr}(n,r):=\left\\{\mathbf{P}\in\mathbb{S}^{2}(\mathbb{R}^{n})\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathbf{P}^{2}=\mathbf{P},\leavevmode\nobreak\
\text{rank}(\mathbf{P})=r\right\\}.$ (5)

Here, $\operatorname{Gr}(n,r)$ is a compact Riemannian submanifold of
$\mathbb{S}^{2}(\mathbb{R}^{n})$ with dimension $r(n-r)$.

######  Remark 2.

No simple characterization of the critical points of $\star$ ‣ 2.2 is known
when $d\geq 3$. This is in contrast to case of $d=2$, where the Eckart-Young
theorem [19] holds. Various possible extensions of Eckart-Young to tensors
have been considered and shown to fail, e.g. see [34].

###  2.3 PGD algorithm

This paper develops a PGD method to solve the constrained optimization problem
$\star$ ‣ 2.2. The PGD framework works as follows: given an initial point, we
move the current iterate $\mathbf{Q}_{t-1}$ along the negative of the
Euclidean gradient direction

$\nabla
F(\mathbf{Q}_{t-1})=2d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{Q}_{t-1},\dots,\mathbf{Q}_{t-1}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q}_{t-1},\dots,\mathbf{Q}_{t-1}\right)\right\rangle_{-1}\in\mathbb{R}^{n\times
r}$ (6)

by a step size of $\gamma_{t}>0$. Then we project the update back to the
Stiefel manifold 4 by selecting the $\mathbf{Q}$ factor from the QR
decomposition. The procedure is described in Algorithm 1. See Section 6 for a
convergence analysis of Algorithm 1.

Algorithm 1 PGD for symmetric Tucker decomposition

Input: symmetric tensor
$\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})\\\ $ number of
iteration T, step sizes $\\{\gamma_{t}\\}\subseteq\mathbb{R}$
initial guess $\mathbf{Q}_{0}\in\operatorname{St}(n,r)$ $\triangleright$
initialization
Output: rank-$r$ orthonormal basis $\mathbf{Q}\in\operatorname{St}(n,r)$

1:for $t=1,\dots,T$ do

2:
$\mathbf{Q}_{t}\leftarrow\mathbf{Q}_{t-1}+2d\gamma_{t}\,\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{Q}_{t-1},\dots,\mathbf{Q}_{t-1}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q}_{t-1},\dots,\mathbf{Q}_{t-1}\right)\right\rangle_{-1}$
$\triangleright$ gradient update

3: $\mathbf{Q}_{t}\leftarrow\text{qr}(\mathbf{Q}_{t},0)$ $\triangleright$
retraction

4:end for

5:return $\mathbf{Q}\leftarrow{\mathbf{Q}}_{T}$

###  2.4 HOEVD approximation

An approximate solution to the best symmetric Tucker decomposition $\star$ ‣
2.2 is the HOEVD [14] result:

$\displaystyle\mathbf{Q}_{\text{hoevd}}\in\operatorname{St}(n,r)\leavevmode\nobreak\
\text{given\leavevmode\nobreak\ by\leavevmode\nobreak\ the\leavevmode\nobreak\
}\text{leading\leavevmode\nobreak\ }r\leavevmode\nobreak\
\text{eigenvectors\leavevmode\nobreak\ of\leavevmode\nobreak\
}\operatorname{mat}(\boldsymbol{\mathscr{X}})\operatorname{mat}(\boldsymbol{\mathscr{X}})^{\top}.$
(7)

HOEVD is widely used for Tucker approximation, as it can be computed via a
direct matrix eigendecomposition.

Generally speaking, the HOEVD solution is not a critical point to the
symmetric tensor cost $\star$ ‣ 2.2 [14]. Its approximation error to a $d$-th
order tensor $\boldsymbol{\mathscr{X}}$ is estimated by [25, Theorems 6.9-10]
as

$\left\|\boldsymbol{\mathscr{X}}-\hat{\boldsymbol{\mathscr{X}}}_{\text{hoevd}}\right\|\leq\sqrt{d}\left\|\boldsymbol{\mathscr{X}}-\hat{\boldsymbol{\mathscr{X}}}_{\text{best}}\right\|.$

Here,
$\hat{\boldsymbol{\mathscr{X}}}_{\text{hoevd}}=\boldsymbol{\mathscr{X}}\cdot(\mathbf{Q}_{\text{hoevd}}\mathbf{Q}_{\text{hoevd}}^{\top},\dots,\mathbf{Q}_{\text{hoevd}}\mathbf{Q}_{\text{hoevd}}^{\top})$
as in Remark 4 and $\hat{\boldsymbol{\mathscr{X}}}_{\text{best}}$ comes from
the argmin of 2. We observe such suboptimality later in Section 4.4. However
according to [13, Section 3.4], HOEVD’s solution usually belongs to an
“attractive” region around the desired local optima in numerical experiments,
so it is a good candidate for initialization.

##  3 Scalable sample moment tensor decomposition

The main interest in this work is the decomposition of the sample moment
tensor $\boldsymbol{\mathscr{M}}^{d}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ 1. We
design scalable algorithms for the symmetric Tucker decomposition $\star$ ‣
2.2 as well as the HOEVD approximation 7 of $\boldsymbol{\mathscr{M}}^{d}$. An
average of symmetric outer products, the sample moment structure 1 is the key
to the scalability of the decomposition computation. Tensor operations with
outer products can be reduced to computations only involving their factor
vectors. Specifically we use the following identities: for vectors
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{n},\mathbf{z}\in\mathbb{R}^{m}$ and a
matrix $\mathbf{Y}\in\mathbb{R}^{n\times m},$ it holds

$\mathbf{x}^{\otimes
d}\cdot\left(\mathbf{Y},\dots,\mathbf{Y}\right)=\left(\mathbf{Y}^{\top}\mathbf{x}\right)^{\otimes
d}\in\mathbb{S}^{d}(\mathbb{R}^{m}),$ (8) $\langle\mathbf{x}^{\otimes
d},\mathbf{y}^{\otimes
d}\rangle=\langle\mathbf{x},\mathbf{y}\rangle^{d}\in\mathbb{R},$ (9)
$\langle\mathbf{z}\otimes\mathbf{x}^{\otimes(d-1)},\mathbf{y}^{\otimes
d}\rangle_{-1}=\operatorname{mat}\left(\mathbf{z}\otimes\mathbf{x}^{\otimes(d-1)}\right)\operatorname{mat}\left(\mathbf{y}^{\otimes
d}\right)^{\\!\top}=\mathbf{z}\,\langle\mathbf{x},\mathbf{y}\rangle^{d-1}\mathbf{y}^{\top}\in\mathbb{R}^{m\times
n}.$ (10)

###  3.1 Scalable PGD algorithm

Here we introduce a scalable version of Algorithm 1 when the input tensor is a
sample moment 1, see pseudocode in Algorithm 2. The symmetric Tucker cost
$\star$ ‣ 2.2 now depends solely on the data observations
$\mathbf{X}=[\mathbf{x}_{1},\dots,\mathbf{x}_{p}]\in\mathbb{R}^{n\times p}.$
Due to 8,

$F(\mathbf{Q}):=\left\|\left(\frac{1}{p}\sum_{i=1}^{p}\mathbf{x}_{i}^{\otimes
d}\right)\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\|^{2}=\left\|\frac{1}{p}\sum_{i=1}^{p}\left(\mathbf{Q}^{\top}\mathbf{x}_{i}\right)^{\otimes
d}\right\|^{2}.$ (11)

We actually do not need to create or store the sample moment. Such an idea is
referred to as the implicit technique in [57, 50].

Commonly, however, there is a huge amount of data in exceedingly high
dimensions. We hence turn to the streaming strategy to speed up computations.
Here data is revealed as a sequence of sample batches
$\\{\mathbf{X}_{1},\dots,\mathbf{X}_{t},\dots\\}\subseteq\mathbb{R}^{n\times
b}$ with batch size $b$ (where $b\ll p$). We read each batch as it arrives but
discard it afterwards. This streaming model is called the turnstile model [48,
41, 60].

Therefore, as the $t$-th $(t\geq 1)$ innovation
$\mathbf{X}_{t}=[\mathbf{x}_{t,1},\dots,\mathbf{x}_{t,b}]^{\top}\in\mathbb{R}^{n\times
b}$ arrives, the cost function 11 is recast stochastically:

$F_{t}(\mathbf{Q}):=\left\|\frac{1}{b}\sum_{i=1}^{b}\left(\mathbf{Q}^{\top}\mathbf{x}_{t,i}\right)^{\otimes
d}\right\|^{2}.$ (12)

In light of 10, we can compute the stochastic gradient in a tensor-free
manner:

$\begin{array}[]{ll}\nabla
F_{t}(\mathbf{Q})&\displaystyle=2d\,\langle\frac{1}{b}\sum_{i=1}^{b}\mathbf{x}_{t,i}\otimes\left(\mathbf{Q}^{\top}\mathbf{x}_{t,i}\right)^{\otimes(d-1)},\frac{1}{b}\sum_{j=1}^{b}\left(\mathbf{Q}^{\top}\mathbf{x}_{t,j}\right)^{\otimes
d}\rangle_{-1}\\\
&\displaystyle=\frac{2d}{b^{2}}\,\sum_{i=1}^{b}\sum_{j=1}^{b}\mathbf{x}_{t,i}\langle\mathbf{Q}^{\top}\mathbf{x}_{t,i},\mathbf{Q}^{\top}\mathbf{x}_{t,j}\rangle^{d-1}\mathbf{x}_{t,j}^{\top}\mathbf{Q}\\\
&\displaystyle=\frac{2d}{b^{2}}\,\mathbf{X}_{t}\left(\mathbf{X}_{t}^{\top}\mathbf{Q}\mathbf{Q}^{\top}\mathbf{X}_{t}\right)^{[d-1]}\mathbf{X}_{t}^{\top}\mathbf{Q}\,\,\in\,\,\mathbb{R}^{n\times
r}.\end{array}$ (13)

Here, $\left(\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ \right)^{[d-1]}$
denotes raising each element of the input matrix to the $(d-1)$-th power.
Computing the gradient implicitly via online samples reduces the cost to
$\mathcal{O}(rnb+rb^{2})$ flops per iteration, improving upon
$\mathcal{O}(rnp+rp^{2})$ using the full amount of samples or a cost that is
exponential in $d$ if we operate on an explicit tensor.

Moreover, a scalable HOEVD solver, s-hoevd in line 1 of Algorithm 2, is
utilized for good initialization; see Section 3.2 for details. We also use
AdaGrad [18] to automatically tune step sizes. Specifically, we implement the
column-wise extension of AdaGrad proposed in [27], since it is natural to
independently update the orthonormal basis vectors in $\mathbf{Q}$. See lines
4 and 5, where $./$ is the columnwise division.

All together, Algorithm 2 is a two-phased iterative algorithm. It consists of

  * •

Phase I: HOEVD approximation (line 1),

  * •

Phase II: symmetric Tucker decomposition (lines 2-7).

The alterations with respect to Algorithm 1 are highlighted in blue.

Algorithm 2 Scalable PGD for sample moment decomposition

Input: numbers of iterations $T_{1},T_{2}$, batch sizes $b_{1},b_{2}$
data streams
$\\{\mathbf{X}_{1},\dots,\mathbf{X}_{T_{1}}\\}\subseteq\mathbb{R}^{n\times
b_{1}}$ for initialization,
$\\{\mathbf{X}_{T_{1}+1},\dots,\mathbf{X}_{T_{2}}\\}\subseteq\mathbb{R}^{n\times
b_{2}}$ for PGD
hyperparameters:
$(\boldsymbol{\gamma}_{1})_{0}=(\boldsymbol{\gamma}_{2})_{0}=10^{-5}\text{ones}(1,r),$
$c_{1},c_{2}>0$
Output: rank-$r$ orthonormal basis $\mathbf{Q}\in\operatorname{St}(n,r)$

1:$\mathbf{Q}_{T_{1}}\leftarrow\textsc{s-hoevd}(T_{1},b_{1},\\{\mathbf{X}_{1},\dots,\mathbf{X}_{T_{1}}\\},(\boldsymbol{\gamma}_{1})_{0},c_{1})$
$\triangleright$ initialization

2:for $t=T_{1}+1,\dots,T_{2}$ do

3:
$\mathbf{G}_{t}\leftarrow\frac{2d}{b_{2}^{2}}\,\mathbf{X}_{t}\left(\mathbf{X}_{t}^{\top}\mathbf{Q}_{t-1}\mathbf{Q}_{t-1}^{\top}\mathbf{X}_{t}\right)^{[d-1]}\mathbf{X}_{t}^{\top}\mathbf{Q}_{t-1}$
$\triangleright$ gradient computation

4:
$(\boldsymbol{\gamma}_{2})_{t}\leftarrow\sqrt{(\boldsymbol{\gamma}_{2})_{t-1}^{[2]}+\text{colnorm}(\mathbf{G}_{t})^{[2]}}$
$\triangleright$ step size

5:
$\mathbf{Q}_{t}\leftarrow\mathbf{Q}_{t-1}+c_{2}\,\mathbf{G}_{t}\,./(\boldsymbol{\gamma}_{2})_{t}$

6: $\mathbf{Q}_{t}\leftarrow\text{qr}(\mathbf{Q}_{t},0)$

7:end for

8:return $\mathbf{Q}\leftarrow{\mathbf{Q}}_{T_{2}}$

###  3.2 Scalable HOEVD

We present the scalable version of HOEVD approximation shown in Algorithm 3.
This serves as the Phase I initialization in Algorithm 2, but it may be of
independent interest. As HOEVD is an eigendecomposition, we follow the
streaming PCA idea. The solver uses the PGD framework, but applied to a
different cost function than 11. Specifically, we solve HOEVD iteratively to
find $\mathbf{Q}\in\operatorname{St}(n,r)$ that maximizes the following
stochastic cost formed by a data batch $\mathbf{X}_{t}\in\mathbb{R}^{n\times
b_{1}}$:

$\left\|\operatorname{mat}\left(\frac{1}{b_{1}}\sum_{i=1}^{b_{1}}\mathbf{x}_{t,i}^{\otimes
d}\right)^{\top}\\!\mathbf{Q}\right\|^{2}\\!.$

Using 10, the gradient can be calculated implicitly:

$\begin{array}[]{ll}\displaystyle
2\,\operatorname{mat}\left(\frac{1}{b_{1}}\sum_{i=1}^{b_{1}}\mathbf{x}_{t,i}^{\otimes
d}\right)\operatorname{mat}\left(\frac{1}{b_{1}}\sum_{i=1}^{b_{1}}\mathbf{x}_{t,i}^{\otimes
d}\right)^{\top}\\!\mathbf{Q}&\displaystyle=\frac{2d}{b_{1}^{2}}\,\sum_{i=1}^{b_{1}}\sum_{j=1}^{b_{1}}\mathbf{x}_{t,i}\langle\mathbf{x}_{t,i},\mathbf{x}_{t,j}\rangle^{d-1}\mathbf{x}_{t,j}^{\top}\mathbf{Q}\\\
&\displaystyle=\frac{2}{b_{1}^{2}}\,\mathbf{X}_{t}\left(\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\right)^{[d-1]}\mathbf{X}_{t}^{\top}\mathbf{Q}\in\mathbb{R}^{n\times
r}.\end{array}$ (14)

It costs $\mathcal{O}(rnb_{1}+nb_{1}^{2}+rb_{1}^{2})$ flops to compute the
gradient. See the whole procedure in Algorithm 3.

Algorithm 3 Scalable HOEVD for sample moment decomposition

Input: number of iterations $T_{1}$, batch size $b_{1}$
data stream
$\\{\mathbf{X}_{1},\dots,\mathbf{X}_{T_{1}}\\}\subseteq\mathbb{R}^{n\times
b_{1}}$
hyperparameters: $(\boldsymbol{\gamma}_{1})_{0}=10^{-5}\text{ones}(1,r),$
$c_{1}>0$
Output: top-$r$ eigenvectors
$\mathbf{Q}_{\text{hoevd}}\in\operatorname{St}(n,r)$

1:function
s-hoevd($\\{\mathbf{X}_{1},\dots,\mathbf{X}_{T_{1}}\\},(\boldsymbol{\gamma}_{1})_{0},c_{1}$)

2:$\mathbf{Q}_{0}\leftarrow[\mathcal{N}(0,1)]^{n\times r}$

3:$\mathbf{Q}_{0}\leftarrow\text{qr}(\mathbf{Q}_{0},0)$ $\triangleright$
initialization

4: for $t=1,\dots,T_{1}$ do

5:
$\mathbf{G}_{t}\leftarrow\frac{2}{b_{1}^{2}}\,\mathbf{X}_{t}\left(\mathbf{X}_{t}^{\top}\mathbf{X}_{t}\right)^{[d-1]}\mathbf{X}_{t}^{\top}\mathbf{Q}_{t-1}$

6:
$(\boldsymbol{\gamma}_{1})_{t}\leftarrow\sqrt{(\boldsymbol{\gamma}_{1})_{t}^{[2]}+\text{colnorm}(\mathbf{G}_{t})^{[2]}}$

7:
$\mathbf{Q}_{t}\leftarrow\mathbf{Q}_{t-1}+c_{1}\,\mathbf{G}_{t}\,./(\boldsymbol{\gamma}_{1})_{t}$
$\triangleright$ gradient update

8: $\mathbf{Q}_{t}\leftarrow\text{qr}(\mathbf{Q}_{t},0)$ $\triangleright$
retraction

9: end for

10:return $\mathbf{Q}_{\text{hoevd}}\leftarrow\mathbf{Q}_{T_{1}}$

11:end function

######  Remark 3.

In Algorithms 2 and 3, the implicit method obviates the need to form the
tensor, which would require $\mathcal{O}(pn^{d})$ flops. Each iteration
temporarily stores an updated matrix $\mathbf{Q}\in\mathbb{R}^{n\times r}$ and
an online sample batch $\mathbf{X}_{t}\in\mathbb{R}^{n\times b}$, taking only
$nr+nb$ memory. This explains the complexity counts shown in Table 1.

######  Remark 4.

The core tensor $\boldsymbol{\mathscr{C}}\in\mathbb{S}^{d}(\mathbb{R}^{r})$ 3
can be evaluated as follows. After obtaining
$\mathbf{Q}\in\operatorname{St}(n,r)$ by Algorithm 2 or 3, we draw fresh data
$\mathbf{X}_{\text{c}}=[(\mathbf{x}_{\text{c}})_{1},\dots,(\mathbf{x}_{\text{c}})_{p_{\text{c}}}]^{\top}\in\mathbb{R}^{n\times
p_{\text{c}}}$ and then set:

$\boldsymbol{\mathscr{C}}=\frac{1}{p_{\text{c}}}\sum_{i=1}^{p_{\text{c}}}\left(\mathbf{Q}^{\top}(\mathbf{x}_{\text{c}})_{i}\right)^{\otimes
d}.$

This costs $\mathcal{O}(r^{d})$ in storage and
$\mathcal{O}(p_{c}rn+p_{c}r^{d})$ in flops.

##  4 Numerical tests

We utilize synthetic tests to study the proposed algorithms for sample moment
tensor decomposition. We aim to show:

  * •

The speedup of implicit compared to explicit optimization algorithms in Figure
2.

  * •

The speedup of streaming compared to full data usage in Figure 3.

  * •

The robustness of the HOEVD approximation and the symmetric Tucker
decomposition to noise and target rank choice in Figures 4 and 5.

###  4.1 Set-up

In these tests, we generate moment tensors with approximate low-rank symmetric
Tucker structure 3. One way to do this is to assume data instance
$\mathbf{x}\in\mathbb{R}^{n}$ follows a linear factor model:

$\mathbf{x}=\mathbf{B}\mathbf{f}+\boldsymbol{\boldsymbol{\varepsilon}}.$ (15)

Here we fix the factor loading matrix $\mathbf{B}\in\mathbb{R}^{n\times
r_{\text{true}}}$. The latent factor
$\mathbf{f}\in\mathbb{R}^{r_{\text{true}}}$ is drawn from a coordinate-
independent, zero-mean and unit-variant Normal-inverse Gaussian (NIG)
distribution [5, 6]. This distribution is skewed and heavy-tailed, which is
aligned with our interest mentioned in Section 1. The Gaussian noise
$\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{n})$ is referred to
${\boldsymbol{\boldsymbol{\varepsilon}}}\in\mathbb{R}^{n}$. Note that the
$d$-th moment of $\mathbf{x}$ roughly admits a low-rank Tucker decomposition
format:

$\begin{array}[]{lll}\boldsymbol{\mathscr{M}}^{d}\\!\\!&\approx\mathbb{E}\left(\mathbf{x}^{\otimes
d}\right)\\!\\!&\approx\,\,\,\displaystyle\mathbb{E}\left((\mathbf{B}\mathbf{f})^{\otimes
d}\right)\\\
&&\displaystyle=\,\,\,\mathbb{E}\left((\mathbf{R}^{*}\mathbf{f})^{\otimes
d}\right)\cdot\left(\mathbf{Q}^{*\top},\dots,\mathbf{Q}^{*\top}\right).\end{array}$

The matrices $\mathbf{Q}^{*}\in\mathbb{R}^{n\times
r_{\text{true}}},\mathbf{R}^{*}\in\mathbb{R}^{r_{\text{true}}\times
r_{\text{true}}}$ are from the QR factorization of $\mathbf{B}$.

We define the inverse signal-to-noise ratio to quantify the noise level of the
model:

$\texttt{SNR}^{-1}:=\frac{\|\mathbb{E}\left(\boldsymbol{\boldsymbol{\varepsilon}}\right)\|}{\|\mathbb{E}\left(\mathbf{B}\mathbf{f}\right)\|}=\frac{\sigma\sqrt{n}}{\|\mathbf{B}\|}.$

######  Remark 5.

We emphasize that our goal is simply to perform symmetric Tucker decomposition
of the sample moment tensor
$\boldsymbol{\mathscr{M}}^{d}\in\mathbb{S}^{d}(\mathbb{R}^{n}).$ This is
different from inference procedures tailored to the linear factor model, e.g.
independent component analysis (ICA) [13].

######  Remark 6.

Besides the model 15 and heavy-tailed data, we note that our proposed
algorithms are able to efficiently decompose $\boldsymbol{\mathscr{M}}^{d}$
where the sample vectors are from any probabilistic models.

All tests are implemented using MATLAB R$2022$a on a MacBook Pro with 16GB RAM
memory. Each reported plot is averaged over $10$ independent simulations. The
matrix $\mathbf{B}$ in 15 is randomly generated with i.i.d. standard normal
entries and then fixed for all runs. We generate new $\mathbf{f}$ and
$\boldsymbol{\boldsymbol{\varepsilon}}$ to form the data samples in each
stochastic simulation.

###  4.2  Implicit vs. explicit algorithms

We illustrate the efficiency of implicit over explicit implementations. We
compare the tensor-free Algorithm 2 to explicit algorithms for $\star$ ‣ 2.2,
namely the explicit Algorithm 1, HOOI [14], quasi-Newton [55] and trust-region
[31] methods. For the last three algorithms, we use the code “tucker-als.m” in
Tensor Toolbox [36], “symalgQNlc.m” in package [54] and “lmlra3-rtr.m” in
TensorLab [63], respectively. HOOI and trust-region apply to general
asymmetric tensors, because there do not exist variants just for symmetric
tensors.

The first experiment is _non-streaming_. The implicit Algorithm 2 runs on a
fixed set of observations $\mathbf{X}\in\mathbb{R}^{n\times p}$. The explicit
approaches require constructing the sample moment tensor
$\boldsymbol{\mathscr{M}}^{d}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ as overhead,
and then operating on it. The moment formation uses the “symktensor” package
in [36]. Due to limited resources for testing explicit methods, the sample
size $p$ and dimension $n$ are taken to be small here. Note that in this test
Algorithm 2 matches Algorithm 1 iteration-by-iteration since there is no
streaming. However, they compute the updates differently. All algorithms start
from the explicit HOEVD solution 7.

Figure 2 shows the running time of Phase II - the symmetric Tucker
decomposition. The vertical axes plot the relative gradient as a measure of
convergence precision, as used in [20, 55, 31]:

$\texttt{relative\leavevmode\nobreak\
gradient}:=\frac{\|\text{grad}\leavevmode\nobreak\
F(\mathbf{Q})\|}{F(\mathbf{Q})}.$

Here
$\operatorname{grad}F(\mathbf{Q})=(\mathbf{I}_{n}-\mathbf{Q}\mathbf{Q}^{\top})\nabla
F(\mathbf{Q})\in\mathbb{R}^{n\times r}$ denotes the Riemannian gradient on the
Stiefel manifold 4 where $\nabla F(\mathbf{Q})\in\mathbb{R}^{n\times r}$ is
the Euclidean gradient 6. We use $\operatorname{grad}F(\mathbf{Q})$, because
the possible update directions on the Stiefel manifold are constrained to the
tangent space at $\mathbf{Q}.$

The number of iterations are set to be: in (a),
$T_{\text{PGD}}=7,T_{\text{HOOI}}=5,T_{\text{QN}}=4,T_{\text{TR}}=8$; in (b),
$T_{\text{PGD}}=8,T_{\text{HOOI}}=6,T_{\text{QN}}=5$. We omit the trust-region
method in (b) since its code does not apply to order-$4$ tensors. The PGD step
size is $c_{2}=1000$ in both (a) and (b).

(a)  (b)

Figure 2: Implicit Algorithm 2 vs. explicit methods (including Algorithm 1).
The relative gradient is plotted against runtime. The size of the problem,
explicit overhead and initialization costs are given in the figure titles. The
markers indicate iterations.

It is evident from Figure 2 that the implicit Algorithm 2 is much faster than
the explicit algorithms. It enjoys an up to $50\times$ speedup compared to the
basic PGD Algorithm 1, before even taking into account the upfront cost in
Algorithm 1.

###  4.3 Streaming vs. full data usage

Here we stress the gains provided by the streaming technique. This is
demonstrated in a stochastic experiment for (a): scalable HOEVD Algorithm 3
and (b): Phase II in Algorithm 2. See Figure 3. The vertical axes represent
the approximation accuracy evaluated by the metric:

$\texttt{subspace\leavevmode\nobreak\
error}:=\|\mathbf{Q}\mathbf{Q}^{\top}-\mathbf{Q}^{*}\mathbf{Q}^{*\top}\|.$

The measure is rotationally invariant.

To compare streaming and full data methods, we first generate a pool of
observations $\mathbf{X}\in\mathbb{R}^{n\times p}$. Each streaming method
selects a subset of samples $\mathbf{X}(:,b(t-1)+1:bt)$ from the pool at step
$t$. We ensure the total sample number $p$ is just enough for the streaming
case with the largest batch size $b_{\max}$. The non-streaming counterparts
use all the samples $\mathbf{X}$ at once. The number of streaming iterations
in both (a) and (b) is $T=250$, and the total sample size is
$p=Tb_{\max}=25,000.$

In (a), the full HOEVD method applies direct eigendecomposition to
$\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{X}\right)\mathbf{X}^{\top}\\!/p^{2}$.
It implicitly computes
$\operatorname{mat}(\boldsymbol{\mathscr{M}}^{d})\operatorname{mat}(\boldsymbol{\mathscr{M}}^{d})^{\top}$,
which is already more efficient than explicit HOEVD. The step size for all
streaming cases is $c_{1}=0.05.$ In (b), all plots start at the same
initialization reached by Algorithm 3 when $T_{1}=20$ and $c_{1}=1.$ The
streaming and full cases use step size $c_{2}=1$ and $10$ respectively.

(a)  (b)

Figure 3: Streaming vs. full data usage: We plot the subspace error against
runtime. The final error values for all methods are marked on the plot.

It is seen in Figure 3 that the streaming Algorithms 3 and 2 have limited loss
of accuracy compared to the methods that use the full dataset. On the other
hand, the runtimes for all streaming cases are nearly negligible compared to
timing for the full algorithms.

###  4.4 Robustness

We test how the noise level in the observations and the choice of target rank
affect HOEVD and the symmetric Tucker decomposition. In the two experiments
below, we plot the evolution of both phases of Algorithm 2, see Figures 4 and
5. As the HOEVD phase (in dashed lines) plateaus, its last iterate gives the
initialization for Phase II (in solid lines). The number of iterations are
shown on the horizontal axes. The parameters in Algorithm 2 are set to
$b_{1}=b_{2}=50,c_{1}=c_{2}=1.$

Consider the algorithms’ sensitivity to noise in the data samples. To test
this, we generate a sequence of clean samples and then plant Gaussian noise
with different levels of $\texttt{SNR}^{-1}$. The algorithms’ results are
shown in Figure 4. In particular, note the improvement in Phase II over the
HOEVD phase. The discrepancy is enlarged in the high noise regime
($\texttt{SNR}^{-1}\geq 1$), when the HOEVD solver barely outputs a qualified
solution. In contrast, relatively good approximation is obtained using the
symmetric Tucker decomposition (Phase II) when $\texttt{SNR}^{-1}$ is no
bigger than $1$.

Next, we test various rank choices in Figure 5. All cases pass over the same
set of data observations. The algorithms’ recovery qualities are reported in
terms of the objective value $F$, and compared to the “ground-truth” value
$F(\mathbf{Q}^{*})$. As there is no deterministic moment
$\boldsymbol{\mathscr{M}}^{d}$ in the streaming case, we estimate $F$ using an
additional set of testing samples
$\mathbf{X}_{\text{test}}\in\mathbb{R}^{n\times p_{\text{test}}}$ drawn from
15. We calculate $F$ implicitly using 9:

$\begin{array}[]{ll}\texttt{objective\leavevmode\nobreak\
value:\leavevmode\nobreak\
}F(\mathbf{Q})&\displaystyle=\frac{1}{p^{2}_{\text{test}}}\left\|\sum_{i=1}^{p_{\text{test}}}\left(\mathbf{Q}^{\top}\mathbf{X}_{\text{test}}\right)_{i}^{\otimes
d}\right\|^{2}=\displaystyle\frac{1}{p^{2}_{\text{test}}}\left\|\left(\mathbf{X}_{\text{test}}^{\top}\mathbf{Q}\mathbf{Q}^{\top}\mathbf{X}_{\text{test}}\right)^{[d]}\right\|.\end{array}$

When $r\geq r_{\text{true}},$ the solutions obtained by Algorithm 2 fully
approximate the moment tensor in the sense that their objective values match
the ground truth. When $r<r_{\text{true}}$, we can only partially recover the
tensor. The less the rank is, the lower the objective value. Also note that as
$r$ decreases, HOEVD’s initialization moves further away from the ground truth
and the gap between the phases is more pronounced.

Figure 4: Two phases with different noise levels.

Figure 5: Two phases with different target ranks.

Figures 4 and 5 show that the robustness of the symmetric Tucker decomposition
is stronger than that of the HOEVD approximation. For “hard” problem instances
where noise is high or the rank is under-specified, the HOEVD solution is
clearly improved upon during Phase II.

##  5 Applications

We illustrate the applicability of moment tensor decompositions to the real
datasets from anomaly detection and portfolio allocation.

###  5.1 Anomaly detection in hyperspectral imaging

Outliers can be discovered using higher-order statistics, when the mean and
covariance information are not sufficient [17, 2]. In hyperspectral imaging
(HSI) [51], the detection goal is to find a low-rank feature subspace, denoted
by $\mathbf{Q}\in\operatorname{St}(n,r)$. We wish that after transforming HSI
pixels by $\mathbf{Q},$ unusual targets in the image are revealed [24].
Intuitively if the anomalies are captured by $\text{colspan}(\mathbf{Q})$,
they would stand out as bright spots in the transformed image.

In this experiment, we test three methods to extract the feature subspace from
HSI data: Tucker factorization of the skewness (order-3 standardized moment)
[43], HOEVD applied to the skewness, and eigendecompostion of the covariance.
Here, skewness and covariance are formed from HSI pixels. The results are
evaluated by their performance in revealing anomalies, as in [24].

Figure 6 shows the data structure of a typical HSI image. It is a $3$-way
hypercube $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{n\times p_{1}\times p_{2}}.$
The data consists of $n$ band images in $\mathbb{R}^{p_{1}\times p_{2}}$ and
$p=p_{1}p_{2}$ HSI pixels in $\mathbb{R}^{n}$.

band $1\in\mathbb{R}^{p_{1}\times p_{2}}$ $\cdots$ band $n$ pixel
$\mathbf{x}_{i}=[x_{i,1},\dots,x_{i,n}]^{\top}\in\mathbb{R}^{n}$
$\bullet$$\bullet$$\bullet$ Figure 6: HSI hypercube.

We use a real hyperspectral image from the Airport-Beach-Urban dataset.222The
dataset is available at http://xudongkang.weebly.com/data-sets.html. The scene
was taken at Gulfport, Florida by the Airborne Visible/Infrared Imaging
Spectrometer (AVIRIS) sensor. The top two images of Figure 7 are the false-
color composite image and the detection reference image (with white spots
depicting the anomalous aircraft). The number of bands is $n=191$ and the size
of each band is $p_{1}\times p_{2}=100\times 100.$ The $d$-th moment would
take $10000\times 191^{d}$ flops and $191^{d}$ space to form and store.
Instead we use scalable algorithms for the moment decomposition task.

First, the 3D HSI cube $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{n\times
p_{1}\times p_{2}}$ is reshaped into an $n\times p$ matrix. This is whitened
[22] to give $\mathbf{X}\in\mathbb{R}^{n\times p}$. We choose rank $r=4$, and
seek the leading rank-$r$ subspaces $\mathbf{Q}\in\mathbb{R}^{n\times r}$ from
the skewness tensor
$\boldsymbol{\mathscr{M}}^{3}\in\mathbb{S}^{3}(\mathbb{R}^{n})$ and the
covariance matrix $\mathbf{M}^{2}\in\mathbb{S}^{2}(\mathbb{R}^{n})$, both
formed using the columns in $\mathbf{X}$. The subspaces are found using
Algorithm 2, Algorithm 3 and PCA, respectively. The parameters for Algorithms
2 and 3 are $b_{1}=b_{2}=100$, $c_{1}=0.35,c_{2}=0.5$ and
$T_{1}=500,T_{2}=1500.$ For each subspace $\mathbf{Q}$ obtained, $\mathbf{X}$
is transformed into $\mathbf{Q}^{\top}\mathbf{X}\in\mathbb{R}^{r\times p}$.
Each row of $\mathbf{Q}^{\top}\mathbf{X}$ is reshaped to a $p_{1}\times p_{2}$
matrix. The resulting bands 1 - 4 are shown in Figure 7.

Since the reference map in Figure 7 is sparse, we further pursue a sparsest
detection image by optimizing over weights on the bands. This is achieved by a
constrained $\ell_{1}$-minimization:

$\mathbf{w}^{*}=\arg\min_{\mathbf{w}\in\mathbb{R}^{r}}\|\mathbf{X}^{\top}\mathbf{Q}\mathbf{w}\|_{1},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{subject\leavevmode\nobreak\ to}\leavevmode\nobreak\ \leavevmode\nobreak\
\|\mathbf{w}\|_{2}=1.$

This is solved by the built-in Matlab tool “fmincon”. Then we reshape
$\mathbf{X}^{\top}\mathbf{Q}\mathbf{w}^{*}$ into a band image. The resulting
sparse images are in the first column of Figure 7.

Inspecting Figure 7, we see that Algorithms 2 and 3 based on skewness
distinguish between abnormal and normal pixels well. The subspace provided by
Algorithm 2 enables especially clear detection, since the targets are the most
pronounced. It also outlines both aircraft individually in the sparse image.
The covariance information, however, fails to identify anomalies.

Figure 7: Qualitative detection comparison.

In Figure 8, we show a quantitative evaluation in anomaly detection, via the
receiver operating characteristic-area under curve (ROC-AUC). This quantifies
how accurate a method is in classifying normal and abnormal objects, where a
higher AUC value indicates better detection. To draw the ROC curve, we whiten
the matrix $\mathbf{X}^{\top}\mathbf{Q}$ and assign the anomaly score for each
pixel as its $\ell_{2}$-norm. The dashed curve indicates a uniformly random
classifier. Note that the ROC of skewness-based detection schemes lie strictly
above that of PCA. The AUC values corresponding to Algorithms 2 and 3 and PCA
are $0.9947,0.9900$ and $0.8893$, respectively.

Figure 8: Quantitative detection comparison.

###  5.2 Moment estimation for portfolio allocation

Moment estimation is sometimes improved by enforcing low-rank structure [6].
In portfolio allocation, we seek a distribution of investments to maximize
future returns. Judicious predictions can be obtained by including higher-
order moments in allocation optimization, as financial data is heavy-tailed
[33].

In this experiment, we test two different ways to estimate moments from
financial data: the empirical sample moment $\boldsymbol{\mathscr{M}}^{d}$ 1
and the low-rank symmetric Tucker approximation to
$\boldsymbol{\mathscr{M}}^{d}$. We evaluate the moments by deriving portfolios
based on them and comparing the portfolios’ performances, as in [6]. We
emphasize that, in theory, higher moments of financial data should roughly
admit a low-rank symmetric Tucker format, because of the Fama-French model
[21].

We use the real HFR hedge fund dataset 333The HFRX indices dataset can be
found at: https://www.hfr.com/family-indices/hfrx.. It consists of daily
return rates (in percentage) of $n=30$ HFRX indices from January 2nd, 2009 to
December 31st, 2019. The time series is split into in-sample data (2009 -
2016) $\mathbf{X}_{\text{in}}\in\mathbb{R}^{n\times
p_{\text{in}}}\leavevmode\nobreak\ (p_{\text{in}}=2016)$ and out-of-sample
data (2017 - 2019) $\mathbf{X}_{\text{out}}\in\mathbb{R}^{n\times
p_{\text{out}}}\leavevmode\nobreak\ (p_{\text{out}}=754)$. The in-sample data
is centered. We calculate the averaged sample moments 1:
$\mathbf{M}^{2}\in\mathbb{S}^{2}(\mathbb{R}^{n}),\boldsymbol{\mathscr{M}}^{3}\in\mathbb{S}^{3}(\mathbb{R}^{n})$
and $\boldsymbol{\mathscr{M}}^{4}\in\mathbb{S}^{4}(\mathbb{R}^{n})$. Then
their low-rank Tucker approximations are computed by Algorithm 2. The rank is
set to $r=15$, and 3 is used to solve for the cores once the subspaces
$\mathbf{Q}$ are determined. We run Phase II of Algorithm 2 thrice. Starting
from a random initialization, we obtain $\mathbf{Q}$ for $\mathbf{M}^{2}$.
This is used an initialization for decomposing $\boldsymbol{\mathscr{M}}^{3}$.
Finally we tackle $\boldsymbol{\mathscr{M}}^{4}$, starting from the
$\mathbf{Q}$ for $\boldsymbol{\mathscr{M}}^{3}$. The parameters are set to
$c_{2}=1,10^{-2},10^{-5}$; $T_{2}=1000,1000,200$; and $b_{2}=500,500,500$ for
the three runs respectively. This produces
$\mathbf{M}_{\operatorname{Tucker}}^{2}\in\mathbb{S}^{2}(\mathbb{R}^{n}),\boldsymbol{\mathscr{M}}_{\operatorname{Tucker}}^{3}\in\mathbb{S}^{3}(\mathbb{R}^{n})$
and
$\boldsymbol{\mathscr{M}}_{\operatorname{Tucker}}^{4}\in\mathbb{S}^{4}(\mathbb{R}^{n})$.

Given moment estimates, we select the allocation weights that maximize the
expected utility (EU) [6, Supplementary Material, Equation (17)]:

$\begin{array}[]{ll}\displaystyle{\mathbf{w}^{*}}:=\arg\max_{\begin{subarray}{c}\sum_{i}w_{i}=1\\\
w_{i}\geq
0\end{subarray}}\text{EU}(\mathbf{w})\,\equiv&\displaystyle-\frac{\mu}{2}\left\langle\mathbb{E}\left((\mathbf{x}-\overline{\mathbf{x}})^{\otimes
2}\right),\mathbf{w}^{\otimes 2}\right\rangle\\\
&\displaystyle+\frac{\mu\,(\mu+1)}{6}\left\langle\mathbb{E}\left((\mathbf{x}-\overline{\mathbf{x}})^{\otimes
3}\right),\mathbf{w}^{\otimes 3}\right\rangle\\\
&\displaystyle-\frac{\mu\,(\mu+1)\,(\mu+2)}{24}\left\langle\mathbb{E}\left((\mathbf{x}-\overline{\mathbf{x}})^{\otimes
4}\right),\mathbf{w}^{\otimes 4}\right\rangle,\end{array}$ (16)

where the risk aversion parameter $\mu$ is set to be 1. We estimate
$\mathbb{E}\left((\mathbf{x}-\overline{\mathbf{x}})^{\otimes
d}\right)\approx\boldsymbol{\mathscr{M}}^{d}$ in 16, and then apply “fmincon”
in Matlab to output the optimal weights. Meanwhile, estimating
$\mathbb{E}\left((\mathbf{x}-\overline{\mathbf{x}})^{\otimes
d}\right)\approx\boldsymbol{\mathscr{M}}_{\operatorname{Tucker}}^{d}$ results
in another weight vector.

The two portfolios are evaluated on the out-sample data. The acquired daily
returns $\\{y_{t}\\}\subseteq\mathbb{R}$ are collected in a vector
$\mathbf{y}=\mathbf{X}_{\text{out}}^{\top}\mathbf{w}^{*}\in\mathbb{R}^{p_{\text{out}}}.$
In Figure 9, we plot the cumulative daily return rates
$(\prod_{t^{\prime}=1}^{t}(1+\frac{y_{t^{\prime}}}{100})-1)\times 100$ at days
$t=1,\dots,p_{\text{out}}$.

Figure 9: Portfolio returns under different weights.

We observe that the low-rank moment estimates produce a better return rate
than the plain sample moments do. Despite mild fluctuations, the two
strategies always produce positive return rates. A naive strategy, the
equally-weighted portfolio, is the most unstable. It sometimes yields negative
returns.

Other statistical measures are shown in Table 2. The annual geometric mean of
returns over the years $2017-2019$ is
$\left(\left(\prod_{t=1}^{p_{\text{out}}}(1+\frac{y_{t}}{100})\right)^{1/3}-1\right)\times
100.$ The centered $d$-th moment ($d=2,3,4$) of daily returns $\\{y_{t}\\}$ is
$m^{d}=\sum_{t=1}^{p_{\text{out}}}\frac{(y_{t}-\overline{y})^{d}}{p_{\text{out}}}\in\mathbb{R}.$

Table 2: Statistical performance of portfolios. | Equally weighted | Sample | Low rank ($r=15$)
---|---|---|---
annual geometric mean | $9.4143\times 10^{-1}$ | $2.2434$ | $\bf 2.3268$
$m^{2}$ | $2.6914\times 10^{-2}$ | $8.9233\times 10^{-3}$ | $\bf 8.5933\times 10^{-3}$
$m^{3}$ | $-2.9724\times 10^{-3}$ | $-7.7193\times 10^{-4}$ | $\bf-6.0692\times 10^{-4}$
$m^{4}$ | $3.4507\times 10^{-3}$ | $5.4339\times 10^{-4}$ | $\bf 4.7906\times 10^{-4}$

Large geometric mean and low absolute values of moments are desirable
allocation outcomes. In each row, we highlight the preferred score in bold.
They are all achieved by the portfolio based on low-rank Tucker approximations
to the moments.

##  6 Convergence analysis

This section presents a convergence theory for the basic PGD method (Algorithm
1) applied to symmetric Tucker decomposition $\star$ ‣ 2.2. Here
$\boldsymbol{\mathscr{X}}$ need not arise as a moment tensor; even when it
does, we do not analyze the effects of streaming. A basic analysis is
warranted first, because $\star$ ‣ 2.2 is non-convex, with a complex landscape
if $d\geq 3$ (Remark 2). We also emphasize that the theory does _not_
immediately follow from generic properties of PGD. We use the specific
structure of the objective function and QR retraction.

###  6.1 Main statements

Recall the set-up in Section 2.2. Because Tucker factorization 2 is non-
unique, as is $\star$ ‣ 2.2 since

$\left\|\boldsymbol{\mathscr{X}}\cdot(\mathbf{Q},\ldots,\mathbf{Q})\right\|^{2}=\left\|\boldsymbol{\mathscr{X}}\cdot(\mathbf{Q}\mathbf{O},\ldots,\mathbf{Q}\mathbf{O})\right\|^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\ r\times
r\leavevmode\nobreak\ \text{orthogonal matrices}\leavevmode\nobreak\
\mathbf{O},$

the convergence of $\left\\{\mathbf{Q}_{t}\right\\}$ is not relevant. Instead
we focus on the convergence of the projectors,
$\mathbf{P}_{t}:=\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\in\operatorname{Gr}(n,r).$
Thus, while Algorithm 1 runs on the Stiefel manifold and generates the
sequence of _bases_ $\\{\mathbf{Q}_{t}\\}\subseteq\operatorname{St}(n,r)$, we
analyze the sequence of _subspaces_
$\\{\mathbf{P}_{t}\\}\subseteq\operatorname{Gr}(n,r)$ in the Grassmannian
manifold.

Given scalar step sizes $\\{\gamma_{t}\\}\subseteq\mathbb{R}_{\geq 0}$, the
update scheme starts from the PGD iteration on $\mathbf{Q}$:

$\mathbf{Q}_{t+1}=\Phi(\mathbf{Q}_{t}):=\framebox{Q}\text{R}\left(\mathbf{Q}_{t}+\gamma_{t}\nabla
F(\mathbf{Q}_{t})\right)\in\operatorname{St}(n,r).$ (17)

The corresponding update on $\mathbf{P}$ is defined as

$\mathbf{P}_{t+1}=\phi(\mathbf{P}_{t}):=\mathbf{Q}_{t+1}\mathbf{Q}_{t+1}^{\top}=\Phi(\mathbf{Q}_{t})\Phi(\mathbf{Q}_{t})^{\top}\in\operatorname{Gr}(n,r).$
(18)

${\mathbf{Q}_{t+1}=\Phi(\mathbf{Q}_{t})\atop\displaystyle\operatorname{St}(n,r)}$${\mathbf{P}_{t+1}=\mathbf{Q}_{t+1}\mathbf{Q}_{t+1}^{\top}=\phi(\mathbf{P}_{t})\atop\displaystyle\operatorname{Gr}(n,r)}$${\displaystyle\operatorname{St}(n,r)\atop\mathbf{Q}_{t}}$${\displaystyle\operatorname{Gr}(n,r)\atop\mathbf{P}_{t}=\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}}$$\Phi$$\phi$
Figure 10: Update diagram.

The sequence $\left\\{\mathbf{P}_{t}\right\\}$ is assessed in terms of the
cost function

$f(\mathbf{P}):=\|\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{P},\dots,\mathbf{P}\right)\\!\|^{2}.$
(19)

As noted in Section 2.2, this is equivalent to the cost in $\star$ ‣ 2.2, as
$f(\mathbf{P})=F(\mathbf{Q})$ when $\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$.
We prove that the iterates $\\{\mathbf{P}_{t}\\}$ converge to a critical point
of $f$ on the Grassmannian. Note that this is stronger than just proving
convergence of the objective values $\\{f(\mathbf{P}_{t})\\}$ in $\mathbb{R}$.

The first and second-order criticality conditions are stated as follows in the
manifold setting.

######  Definition 7.

A point $\mathbf{P}\in\operatorname{Gr}(n,r)$ is first-order critical for the
cost function $f$ 19 if

$\operatorname{grad}f(\mathbf{P})={\bf 0}.$ (20)

Here,
$\operatorname{grad}f(\mathbf{P})\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\subseteq\mathbb{S}^{2}(\mathbb{R}^{n})$
denotes the Riemannian gradient, which is defined as the Euclidean gradient
$\nabla f(\mathbf{P})\in\mathbb{S}^{2}(\mathbb{R}^{n})$ projected to
$\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$. The tangent space at
$\mathbf{P}$ is given by

$\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r):=\left\\{\Delta\mathbf{P}\in\mathbb{S}^{2}(\mathbb{R}^{n})\leavevmode\nobreak\
|\leavevmode\nobreak\
\mathbf{P}\Delta\mathbf{P}+(\Delta\mathbf{P})\mathbf{P}=\Delta\mathbf{P},\,\text{trace}(\Delta\mathbf{P})=0\right\\}.$
(21)

######  Definition 8.

A point $\mathbf{P}\in\operatorname{Gr}(n,r)$ is second-order critical for the
cost function $f$ 19 if

$\operatorname{grad}f(\mathbf{P})={\bf 0}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\left\langle\Delta\mathbf{P},\operatorname{Hess}f(\mathbf{P})[\Delta\mathbf{P}]\right\rangle\leq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\
\leavevmode\nobreak\
\Delta\mathbf{P}\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r).$ (22)

Here,
$\operatorname{Hess}f(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
denotes the Riemannian Hessian. It is calculated in [26, Eqn. (2.109)].

The first theoretical result is a global guarantee: Algorithm 1 always
converges to a first-order critical point.

######  Theorem 6.1.

Let $\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ be any
symmetric tensor, and consider the objective function $f$ defined by 19. There
exists a constant $\Gamma^{*}(\boldsymbol{\mathscr{X}},r)>0$ such that if the
step sizes $\left\\{\gamma_{t}\right\\}\subseteq\mathbb{R}_{\geq 0}$ satisfy
$\sup_{t}\gamma_{t}\leq\Gamma^{*}$ and $\inf_{t}\gamma_{t}>0$, the following
holds. For all initializations $\mathbf{Q}_{0}\in\operatorname{St}(n,r)$, the
sequence $\left\\{\mathbf{P}_{t}\right\\}\subseteq\operatorname{Gr}(n,r)$
generated via 17-18 by Algorithm 1 converges monotonically to a first-order
critical point 20 of $f$. The convergence is at no less than an algebraic
rate.

Next, in fact, Algorithm 1 converges to a second-order critical point, for
almost all initializations.

######  Theorem 6.2.

Let $\boldsymbol{\mathscr{X}}\in\mathbb{S}^{d}(\mathbb{R}^{n})$ be any
symmetric tensor, and consider the cost $f$ defined by 19. Use a constant step
size $\gamma$ for all iterates. There exist a constant
$\Gamma^{**}(\boldsymbol{\mathscr{X}},r)>0$ and a measure-zero subset
$\mathcal{I}\subseteq\operatorname{St}(n,r)$ such that if
$0<\gamma\leq\Gamma^{**}$, then for all initializations
$\mathbf{Q}_{0}\in\operatorname{St}(n,r)\setminus\mathcal{I}$, the sequence
$\left\\{\mathbf{P}_{t}\right\\}\subseteq\operatorname{Gr}(n,r)$ generated via
17-18 by Algorithm 1 converges to a second-order critical point 22 of $f$.

We build the proofs of Theorems 6.1 and 6.2 respectively in Sections 6.2 and
6.3. The main technical tools are a convergence theorem for real-analytic cost
functions and the center-stable manifold theorem. These are stated in Appendix
G.

###  6.2 Monotonic convergence to first-order critical points

To establish Theorem 6.1, we apply the guarantee stated in Theorem G.1,
originally from [56, Theorem 2.3]. It requires us to prove
$\\{\mathbf{P}_{t}\\}$ meets the conditions in Theorem G.1.

We first compute the Riemannian gradient of the cost $f$ 19 on
$\operatorname{Gr}(n,r)$. Although the conclusion of Theorem 6.1 is in terms
of $\mathbf{P}\in\operatorname{Gr}(n,r)$, our proof uses
$\mathbf{Q}\in\operatorname{St}(n,r)$ due to the explicit update rule in 17.
Thus, we state the relation between $\operatorname{grad}f(\mathbf{P})$ and
$\mathbf{Q},\nabla F(\mathbf{Q})$. As short-hand, define
$w:\operatorname{Gr}(n,r)\to\mathbb{S}^{2}(\mathbb{R}^{n})$ by

$w(\mathbf{P})=d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\in\mathbb{S}^{2}(\mathbb{R}^{n}).$
(23)

Note that $w(\mathbf{P})$ is a positive semidefinite (PSD) matrix.

######  Proposition 6.3.

For $\mathbf{Q}\in\operatorname{St}(n,r)$ and
$\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}\in\operatorname{Gr}(n,r),$ the
following holds

$\operatorname{grad}f(\mathbf{P})=\operatorname{sym}\Big{(}2\left(\mathbf{I}_{n}-\mathbf{P}\right)w(\mathbf{P})\mathbf{P}\Big{)}=\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}\mathbf{Q}^{\top}\right)\nabla
F(\mathbf{Q})\mathbf{Q}^{\top}\Big{)}.$ (24)

Proposition 6.3 is proven in Appendix A.

Next, we analyze the update rule on $\mathbf{P}\in\operatorname{Gr}(n,r)$. For
$t\geq 0$ and step size $\gamma_{t}>0$, define
$\alpha_{t}:\operatorname{Gr}(n,r)\to\mathbb{S}^{2}(\mathbb{R}^{n})$ by

$\alpha_{t}(\mathbf{P})=\mathbf{P}+2\gamma_{t}\mathbf{P}w(\mathbf{P})\mathbf{P}\in\mathbb{S}^{2}(\mathbb{R}^{n}).$
(25)

The expression $\alpha_{t}(\mathbf{P})$ helps to represent the update on
$\mathbf{P}$.

######  Lemma 6.4.

The difference between the current iterate $\mathbf{P}_{t}$ and the next
$\mathbf{P}_{t+1}$ 17-18 is approximated as follows:

$\begin{array}[]{ll}\mathbf{P}_{t+1}-\mathbf{P}_{t}&=4\,\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
+\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|^{2}\right)+\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right).\end{array}$
(26)

Lemma 6.4 is proven in Appendix B.

The following lemma gives properties of $\alpha_{t}(\mathbf{P})$ and shows
that when the step size is sufficiently small,
$\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P})^{\dagger}\right\|$
and $\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|$ are on the same
order.

######  Lemma 6.5.

Let $t\geq 0$ and $\mathbf{P}\in\operatorname{Gr}(n,r)$. The pseudoinverse
$\alpha_{t}(\mathbf{P})^{\dagger}\in\mathbb{S}^{2}(\mathbb{R}^{n})$ satisfies
$\alpha_{t}(\mathbf{P})^{\dagger}\alpha_{t}(\mathbf{P})=\mathbf{P}.$ Moreover,
there exists a constant $\Gamma_{1}(\boldsymbol{\mathscr{X}},r)>0$ such that
if $\gamma_{t}$ satisfies $0<\gamma_{t}\leq\Gamma_{1},$ then
$\|\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}\|\leq\frac{1}{2}$ and

$\frac{\sqrt{2}}{3}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|\leq\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\|\leq\frac{\sqrt{2r}}{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|.$
(27)

Please refer to Appendix C for the proof of Lemma 6.5.

###### Proof of Theorem 6.1.

With the above useful lemmas and proposition in place, in this proof, we show
that our update scheme for $\mathbf{P}\in\operatorname{Gr}(n,r)$ 17-18
satisfies the conditions C1 \- C4 in Theorem G.1. Fix a step $t\geq 0$. Assume
$0<\gamma_{t}\leq\Gamma_{1}$ where $\Gamma_{1}$ is as in Lemma 6.5.

  1. (C1)

We apply 27 that
$\mathcal{O}\left(\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|\right)$
has the same scale as
$\mathcal{O}\left(\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|\right)$.
By Lemma 6.4, it follows

$\mathbf{P}_{t+1}-\mathbf{P}_{t}=4\,\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}+\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right).$
(28)

It implies that there exists a constant $C_{1}(\boldsymbol{\mathscr{X}},r)>0$
such that

$\left\|\mathbf{P}_{t+1}-\mathbf{P}_{t}-4\,\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\right\|\leq
C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}.$

By 51 in Appendix B,

$\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}=(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\mathbf{P}_{t}\alpha_{t}(\mathbf{P}_{t})^{\dagger}=(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger},$

and so

$\begin{array}[]{l}\left\langle\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger},\alpha_{t}(\mathbf{P}_{t})^{\dagger}\operatorname{grad}f(\mathbf{P}_{t})\right\rangle\\\
=\text{trace}\left((\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right)=0.\end{array}$

It implies

$\left\|\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\right\|=\displaystyle\frac{\sqrt{2}}{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|.$

By the triangle inequality,

$\begin{array}[]{ll}\left\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\right\|&\displaystyle\geq
4\times\frac{\sqrt{2}}{2}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|-C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\\\
&\geq\displaystyle
2\sqrt{2}\times\frac{\sqrt{2}}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|-C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\\\
&=\displaystyle\frac{4}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|-C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}.\end{array}$
(29)

Here, the second inequality is due to the first half of 27. We further
constrain the step size as follows

$0<\gamma_{t}\leq\Gamma_{2}:=\min\left(\Gamma_{1},\frac{\frac{2}{3}}{C_{1}\sup_{\mathbf{P}\in\operatorname{Gr}(n,r)}\|\operatorname{grad}f(\mathbf{P})\|}\right)\\!,$
(30)

so that

$C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\leq\frac{2}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|.$
(31)

Note that
$\sup_{\mathbf{P}\in\operatorname{Gr}(n,r)}\|\operatorname{grad}f(\mathbf{P})\|$
is finite because $\operatorname{grad}f$ is continuous on
$\operatorname{Gr}(n,r)$ and the Grassmann is compact. (We also exclude the
trivial case when $\operatorname{grad}f(\mathbf{P})\equiv 0$.) Continuing from
29 and substituting 31, we have

$\left\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\right\|\geq(\frac{4}{3}-\frac{2}{3})\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|=\frac{2}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|.$

Letting $\kappa=\inf_{t\to\infty}\frac{2}{3}\gamma_{t}>0$ ensures that C1 76
holds.

In preparation for C2, we derive an upper bound for
$\left\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\right\|$. The idea is similar to 29:

$\begin{array}[]{ll}\left\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\right\|&\leq
2\sqrt{2}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|+C_{1}\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\\\
&\displaystyle\leq
2\sqrt{2}\times\frac{\sqrt{2r}}{2}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|+\frac{2}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|\\\
&\displaystyle\leq\frac{8\sqrt{r}}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|.\end{array}$
(32)

The second inequality is obtained from the latter half of 27 and 31.

  2. (C2)

Keep the assumption 30 on $\gamma_{t}$. Consider the first-order Taylor’s
expansion for $f$ on $\operatorname{Gr}(n,r)$ [7, Section 10.7]. Inserting 28,
the difference of the current and next $f$ is:

$\begin{array}[]{ll}f(\mathbf{P}_{t+1})-f(\mathbf{P}_{t})&=\langle\operatorname{grad}f(\mathbf{P}_{t}),\mathbf{P}_{t+1}-\mathbf{P}_{t}\rangle+\mathcal{O}\left(\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\|^{2}\right)\\\
&=\left\langle\operatorname{grad}f(\mathbf{P}_{t}),4\gamma_{t}\,\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\right\rangle\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
+\mathcal{O}\left(\gamma_{t}^{2}\max\left(\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{3},\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right)\right).\end{array}$
(33)

The first term in the result of 33, since
$\operatorname{grad}f(\mathbf{P}_{t})$ is symmetric, is

$\begin{array}[]{l}\left\langle\operatorname{grad}f(\mathbf{P}_{t}),4\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\right\rangle\\\
=4\gamma_{t}\left\langle\operatorname{grad}f(\mathbf{P}_{t}),\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\rangle\\\
=4\gamma_{t}\left\langle\operatorname{grad}f(\mathbf{P}_{t})\mathbf{P}_{t},\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}\right\rangle\\\
=4\gamma_{t}\left\langle\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\alpha_{t}(\mathbf{P}_{t}),\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}\right\rangle.\end{array}$

The second equation is due to
$\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}=\alpha_{t}(\mathbf{P}_{t})^{\dagger}.$
The definition 25 gives
$\alpha_{t}(\mathbf{P}_{t})=\mathbf{P}_{t}(\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P}_{t}))\mathbf{P}_{t}$.
Continuing from above,

$\begin{array}[]{l}\left\langle\operatorname{grad}f(\mathbf{P}_{t}),4\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}\right\rangle\\\
=4\gamma_{t}\left\langle\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}(\mathbf{I}_{n}+2w(\mathbf{P}_{t})),\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}\right\rangle\\\
\geq
4\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\mathbf{P}_{t}\right\|^{2}=4\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|^{2}\\\
\geq\displaystyle\frac{8}{9}\gamma_{t}\|\operatorname{grad}f(\mathbf{P}_{t})\|^{2}.\end{array}$
(34)

Here we used that
$\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P}_{t})\succeq\mathbf{I}_{n}$ for
$\gamma_{t}\geq 0$, since $w(\mathbf{P}_{t})\in\mathbb{S}^{2}(\mathbb{R}^{n})$
is PSD. The second inequality is derived from 27.

Considering the remainder term in 33, there is a constant
$C_{2}(\boldsymbol{\mathscr{X}},r)>0$ such that

$f(\mathbf{P}_{t+1})-f(\mathbf{P}_{t})\geq\frac{8}{9}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}-C_{2}\gamma_{t}^{2}\max\left(\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{3},\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right).$
(35)

If the step size $\gamma_{t}$ is small enough to satisfy

$\begin{array}[]{ll}0<\gamma_{t}\\!\\!\\!&\leq\Gamma^{*}:=\min\left(\Gamma_{2},\frac{4}{9C_{2}}\right),\end{array}$
(36)

where $\Gamma_{2}>0$ is defined in 30, then

$C_{2}\gamma_{t}^{2}\max\left(\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{3},\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right)\leq\frac{4}{9}\gamma_{t}\left\|(\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}.$

Inserting into 35,

$\begin{array}[]{ll}f(\mathbf{P}_{t+1})-f(\mathbf{P}_{t})&\displaystyle\geq\left(\frac{8}{9}-\frac{4}{9}\right)\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\\\
&=\displaystyle\frac{1}{6\sqrt{r}}\left(\frac{8\sqrt{r}}{3}\gamma_{t}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|\right)\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|\\\\[10.00002pt]
&\geq\displaystyle\frac{1}{6\sqrt{r}}\,\|\mathbf{P}_{t+1}-\mathbf{P}_{t}\|\|\operatorname{grad}f(\mathbf{P}_{t})\|.\end{array}$

The final inequality is by 32. Thus C2 holds if we set
$\sigma=\frac{1}{6\sqrt{r}}$ in 77.

  3. (C3)

“$\Rightarrow$”: From $\operatorname{grad}f(\mathbf{P})=\bf 0$ and 24, we know

$\begin{array}[]{ll}\operatorname{grad}f(\mathbf{P})\mathbf{Q}={\bf
0}=\operatorname{grad}F(\mathbf{Q})=(\mathbf{I}_{n}-\mathbf{Q}\mathbf{Q}^{\top})\nabla
F(\mathbf{Q}),\end{array}$

where $\mathbf{Q}\in\operatorname{St}(n,r)$ is a basis of
$\mathbf{P}\in\operatorname{Gr}(n,r)$. It implies that $\text{colspan}(\nabla
F(\mathbf{Q}))$ is in $\text{colspan}(\mathbf{Q}).$ Note that
$\mathbf{Q}+\gamma\nabla F(\mathbf{Q})=\left(\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\right)\mathbf{Q}$ has full column rank, like $\mathbf{Q}$,
because $\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\in\mathbb{S}^{2}(\mathbb{R}^{n})$ is positive definite. Hence
$\text{colspan}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)=\text{colspan}(\mathbf{Q})$. Then the next iterate
$\Phi(\mathbf{Q})$ obtained by QR decomposition still outputs a basis for
$\text{colspan}(\mathbf{Q}).$ So the orthogonal projectors remain the same,
i.e. $\phi(\mathbf{P})=\mathbf{P}$.

“$\Leftarrow$”: This directly follows from C1 if $\gamma_{t}$ is sufficiently
small to satisfy 30.

  4. (C4)

Note the Grassmannian is a real-analytic submanifold of
$\operatorname{Gr}(n,r)$, because $\operatorname{Gr}(n,r)$ satisfies property
$(3)$ in [37, Proposition 2.7.3]. Also, the cost function $f$ 19 is a
polynomial in $\mathbf{P}\in\operatorname{Gr}(n,r)$, so real-analytic as well.
By [56, Proposition 2.2], all points $\mathbf{P}\in\operatorname{Gr}(n,r)$
satisfy a Łojasiewicz inequality for $f$.

We have shown that the sequence $\\{\mathbf{P}_{t}\\}$ meets all of the
conditions C1 - C4 provided $\\{\gamma_{t}\\}$ satisfies $\inf\gamma_{t}>0$
and $\sup\gamma_{t}\leq\Gamma^{*}$, with $\Gamma^{*}$ given in 36. By Theorem
G.1, it follows Algorithm 1 achieves convergence to a first-order critical
point of $f$. The proof of Theorem 6.1 is complete. ∎

###  6.3 Convergence to second-order critical points

We study second-order criticality when the step size is constant. Theorem 6.2
relies on the center-stable manifold theorem [58], which we recall in Theorem
G.2. It describes the local structure of a smooth local diffeomorphism around
a fixed point.

The calculation of the Riemannian Hessian,
$\operatorname{Hess}f(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$,
and the differential of the update map,
$\text{D}\phi(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\text{T}_{\phi(\mathbf{P})}\operatorname{Gr}(n,r)$,
are essential to the second-order guarantee. As short-hand, define a linear
operator
$v(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\mathbb{R}^{n\times
n}$ by

$v(\mathbf{P})[\Delta\mathbf{P}]=d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\Delta\mathbf{P},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{P}\dots,\mathbf{P}\right)\right\rangle_{-1}\mathbf{P}\in\mathbb{R}^{n\times
n}.$ (37)

######  Proposition 6.6.

For $\mathbf{P}\in\operatorname{Gr}(n,r)$ and
$\Delta\mathbf{P}\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r),$ the
Riemannian Hessian of the cost $f$ 19 on $\operatorname{Gr}(n,r)$ is given by

$\operatorname{Hess}f(\mathbf{P})[\Delta\mathbf{P}]=2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\mathbf{P})v(\mathbf{P})[\Delta\mathbf{P}]+w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}-\Delta\mathbf{P}w(\mathbf{P})\mathbf{P}\Big{)}.$
(38)

The differential of $\phi$ 18 on $\operatorname{Gr}(n,r)$ is given by

$\text{D}\phi(\mathbf{P})[\Delta\mathbf{P}]=2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\phi(\mathbf{P}))\left(\Delta\mathbf{P}+2\gamma
w(\mathbf{P})\Delta\mathbf{P}+2\gamma
v(\mathbf{P})[\Delta\mathbf{P}]\right)(\mathbf{P}+2\gamma
w(\mathbf{P})\mathbf{P})^{{\dagger}}\Big{)}.$ (39)

Proposition 6.6 is proven in Appendix D.

We next relate the eigenvalues of $\operatorname{Hess}f(\mathbf{P})$ and
$\text{D}\phi(\mathbf{P})$ when $\mathbf{P}$ is a critical point of $f$. In
this case, $\mathbf{P}$ is a fixed point of $\phi$ by C3 in the proof of
Theorem 6.1.

######  Lemma 6.7.

Let $\mathbf{P}$ be a first-order critical point of $f$ 20. If
$\operatorname{Hess}f(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
has a positive eigenvalue, then
$\text{D}\phi(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\to\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
has an eigenvalue which is greater than $1$.

Lemma 6.7 implies that any first-order, but not second-order, critical point
of $f$ is an unstable fixed point of $\phi$. Its proof is in Appendix E.

The following lemma guarantees that the update map $\phi$ has a nice local
property given a sufficiently small step size.

######  Lemma 6.8.

There exists a constant
$0<\Gamma^{**}(\boldsymbol{\mathscr{X}},r)\leq\Gamma^{*}$ (cf. Theorem 6.1)
such that if the step size satisfies $0<\gamma\leq\Gamma^{**}$, then the
update map $\phi:\operatorname{Gr}(n,r)\to\operatorname{Gr}(n,r)$ 18 is a
local diffeomorphism.

The proof of Lemma 6.8 is in Appendix F.

###### Proof of Theorem 6.2.

Fix $0<\gamma\leq\Gamma^{**}$. From Theorem 6.1, we know that the sequence
$\\{\mathbf{P}_{t}\\}$ always converges to a first-order critical point of $f$
20, thus a fixed point of $\phi$ 18 by C3 in Section 6.2. In this proof, we
consider the case when the limit $\mathbf{P}\in\operatorname{Gr}(n,r)$ is a
“bad” critical point, meaning it is a first-order critical point but fails
second-order criticality 22. The idea is to apply Theorem G.2 to show that the
set of $\mathbf{P}_{0}$ that converge to bad critical points has zero measure
in the Grassmannian. We then pull this back to $\operatorname{St}(n,r)$.

Because $\phi$ is a local diffeomorphism by Lemma 6.8, we are eligible to
apply Theorem G.2. It gives an open neighborhood
$B_{\mathbf{P}}\subseteq\operatorname{Gr}(n,r)$ around each bad critical point
$\mathbf{P}$ satisfying the conditions in Theorem G.2. As
$\operatorname{Gr}(n,r)$ is second-countable, it has the Lindelöf property,
i.e. every open cover has a countable subcover. So, there exists a countable
set $\mathcal{C}\subseteq\operatorname{Gr}(n,r)$ consisting of some of the bad
critical points such that

$\bigcup_{\mathbf{P}\in\operatorname{Gr}(n,r)\text{\leavevmode\nobreak\
is\leavevmode\nobreak\ a\leavevmode\nobreak\ bad\leavevmode\nobreak\
critical\leavevmode\nobreak\
point}}B_{\mathbf{P}}\quad\quad=\quad\quad\bigcup_{\mathbf{P}\in\mathcal{C}}B_{\mathbf{P}}.$
(40)

Define the subset

$\mathcal{J}\,:=\,\bigcup_{\mathbf{P}\in\mathcal{C}}\bigcup_{t=0}^{\infty}\phi^{-t}(W_{\mathbf{P}})\,\,\,\subseteq\,\,\,\operatorname{Gr}(n,r).$
(41)

We claim that

$\left\\{\mathbf{P}_{0}\in\operatorname{Gr}(n,r):\phi^{t}(\mathbf{P}_{0})\text{\leavevmode\nobreak\
converges\leavevmode\nobreak\ to\leavevmode\nobreak\ a\leavevmode\nobreak\
bad\leavevmode\nobreak\ critical\leavevmode\nobreak\ point\leavevmode\nobreak\
as\leavevmode\nobreak\
}t\rightarrow\infty\right\\}\,\,\subseteq\,\,\mathcal{J}.$ (42)

To see this, assume
$\lim_{t\rightarrow\infty}\phi^{t}(\mathbf{P}_{0})=\mathbf{P}$ is bad. By 40,
$\mathbf{P}\in B_{\mathbf{P}^{\prime}}$ for some
$\mathbf{P}^{\prime}\in\mathcal{C}$. Since $B_{\mathbf{P}^{\prime}}$ is open,
there is $t_{0}\geq 0$ such that $\phi^{t}(\mathbf{P}_{0})\in
B_{\mathbf{P}^{\prime}}$ for all $t\geq t_{0}$. As $\mathbf{P}^{\prime}$ is a
fixed point of $\phi$, by 80 in Theorem G.2, we obtain
$\phi^{t_{0}}(\mathbf{P}_{0})\in W_{\mathbf{P}^{\prime}}$. In other words,
$\mathbf{P}_{0}\in\phi^{-t_{0}}(W_{\mathbf{P}^{\prime}})$. So,
$\mathbf{P}_{0}\in\mathcal{J}$ as claimed.

We claim that $\mathcal{J}$ is measure-zero in $\operatorname{Gr}(n,r)$. For
each $\mathbf{P}\in\mathcal{C},$ at least one eigenvalue of
$\operatorname{Hess}f(\mathbf{P})$ is greater than 0. By Lemma 6.7, we know at
least one eigenvalue of $\text{D}\phi(\mathbf{P})$ is greater than 1. By the
first result of Theorem G.2, the dimension of $W_{\mathbf{P}}$ is strictly
less than the dimension of $\operatorname{Gr}(n,r)$. So $W_{\mathbf{P}}$ is
measure-zero in $\operatorname{Gr}(n,r)$ by [39, Corollary 6.12]. Given that
the step size satisfies $0<\gamma\leq\Gamma^{**}$, the map $\phi$ is a local
diffeomorphism by Lemma 6.8, and so is $\phi^{t}$ for all $t\geq 0$. Thus
their pre-images of measure-zero sets have zero measure. As
$(\phi^{t})^{-1}\equiv\phi^{-t}$, the set $\phi^{-t}(W_{\mathbf{P}})$ has zero
measure for each $t\geq 0$. Then $\mathcal{J}$ is a countable union of
measure-zero sets, so has zero measure.

To conclude, define $\mathcal{I}\subseteq\operatorname{St}(n,r)$ as the pre-
image of $\mathcal{J}\subseteq\operatorname{Gr}(n,r)$ under the map
$\mathbf{Q}\mapsto\mathbf{Q}\mathbf{Q}^{\top}$. Since the map is a submersion
from the Stiefel manifold onto the Grassmannian, $\mathcal{I}$ has zero
measure in $\operatorname{St}(n,r)$ as $\mathcal{J}$ does in
$\operatorname{Gr}(n,r)$. For all initializations $\mathbf{Q}_{0}$ in
$\operatorname{St}(n,r)$ but outside the null set $\mathcal{I}$, the sequence
$\\{\phi^{t}(\mathbf{P}_{0})\\}=\\{\mathbf{P}_{t}\\}$ converges to a second-
order critical point of $f$ if $0<\gamma\leq\Gamma^{**}$. The proof of Theorem
6.2 is complete. ∎

##  7 Conclusion and future work

In this paper, we developed the PGD method (Algorithm 1) for symmetric Tucker
tensor decomposition. It has a simple update scheme compared to other
iterative solvers running on Riemannian manifolds. We further designed a
scalable PGD (Algorithm 2) and HOEVD solver (Algorithm 3) to decompose sample
moment tensor 1. The two algorithms are free of constructing and storing
moment tensors; the iterations are conducted by streaming only a small subset
of samples. Tremendous computational savings were seen in numerical
experiments without sacrifices on accuracy. We also showed the applicability
of the Tucker decomposition of moment tensors to detect anomalies and allocate
assets. Finally, utilizing manifold optimization theory, we derived
theoretical guarantees for the PGD sequence on the Grassmannian. We proved it
always converges to a first-order critical point, and almost surely to a
second-order critical point.

There are a number of directions that would be interesting to research in the
future.

  1. 1.

Implementation for cumulants. Cumulant tensors are central to non-Gaussian
data analysis. They are nonlinear combinations of moments. It would be useful
to derive implicit and streaming implementations of our algorithms for sample
cumulants.

  2. 2.

Extending the analysis. It is certainly worth investigating convergence
guarantees in the streaming setting of Algorithms 2 and 3. We also want to
study convergence properties when the algorithms are led by adaptive step
sizes.

  3. 3.

Estimating the core without fresh samples. Currently, the core tensor
$\boldsymbol{\mathscr{C}}$ is computed using independent data samples after
the subspace $\mathbf{Q}$ has been calculated. We would like to figure out a
good way of updating $\boldsymbol{\mathscr{C}}$ while updating $\mathbf{Q}$,
using the same data stream.

## Acknowledgements

The authors are grateful to Dries Cornilly and Xiurui Geng for their help with
the real data applications. We also thank Hemanth Kolla, Amit Singer and Mihai
Sirbu for useful discussions.

R.J. was supported by NSF DMS 1952735. J.K. was supported by start-up grants
from the College of Natural Sciences and Oden Institute for Computatational
Engineering and Sciences at UT Austin. R.W. was supported by AFOSR MURI
FA9550-19-1-0005, NSF DMS 1952735, NSF HDR-1934932 and NSF 2019844.

## References

  * [1] P.-A. Absil, Robert Mahony, and Rodolphe Sepulchre.  Optimization Algorithms on Matrix Manifolds, volume 78.  Dec. 2008.
  * [2] Konduri Aditya, Hemanth Kolla, W. Philip Kegelmeyer, Timothy M. Shead, Julia Ling, and Warren L. Davis.  Anomaly detection in scientific data using joint statistical moments.  Journal of Computational Physics, 387:522–538, 2019.
  * [3] Salman Ahmadi-Asl, Stanislav Abukhovich, Maame G. Asante-Mensah, Andrzej Cichocki, Anh Huy Phan, Tohishisa Tanaka, and Ivan Oseledets.  Randomized algorithms for computation of Tucker decomposition and higher order SVD (HOSVD).  IEEE Access, 9:28684–28706, 2021.
  * [4] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky.  Tensor decompositions for learning latent variable models.  Journal of Machine Learning Research, 15(80):2773–2832, 2014.
  * [5] Ole E. Barndorff-Nielsen.  Normal inverse Gaussian distributions and stochastic volatility modelling.  Scandinavian Journal of Statistics, 24(1):1–13, 1997.
  * [6] Kris Boudt, Dries Cornilly, and Tim Verdonck.  Nearest comoment estimation with unobserved factors.  Journal of Econometrics, 217:381–397, 2020.
  * [7] Nicolas Boumal.  _An Introduction to Optimization on Smooth Manifolds_.  To appear in Cambridge University Press, Jan. 2022.
  * [8] Walter Briec, Kristiaan Kerstens, and Octave Jokung.  Mean-variance-skewness portfolio performance gauging: A general shortage function and dual approach.  Management Science, 53(1):135–149, 2007.
  * [9] Jean-Francois Cardoso.  Super-symmetric decomposition of the fourth-order cumulant tensor. blind identification of more sources than sensors.  In Proceedings ICASSP 91: 1991 International Conference on Acoustics, Speech, and Signal Processing, pages 3109–3112, vol.5, 1991.
  * [10] Jean-Francois Cardoso.  Blind signal separation: Statistical principles.  Proceedings of the IEEE, 86(10):2009–2025, 1998.
  * [11] Maolin Che, Yimin Wei, and Hong Yan.  The computation of low multilinear rank approximations of tensors via power scheme and random projection.  SIAM Journal on Matrix Analysis and Applications, 41(2):605–636, 2020.
  * [12] Pierre Comon.  Independent component analysis, a new concept?  Signal Processing, 36(3):287–314, 1994.
  * [13] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  An introduction to independent component analysis.  Journal of Chemometrics, 14(3):123–149, 2000.
  * [14] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  A multilinear singular value decomposition.  SIAM Journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000.
  * [15] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  On the best rank-1 and rank-($R_{1},R_{2},...,R_{n}$) approximation of higher-order tensors.  SIAM Journal on Matrix Analysis and Applications, 21(4):1324–1342, 2000.
  * [16] Vin de Silva and Lek-Heng Lim.  Tensor rank and the ill-posedness of the best low-rank approximation problem.  SIAM Journal of Matrix Analysis and Applications, 30:1084–1127, 2008.
  * [17] Krzysztof Domino, Piotr Gawron, and Lukasz Pawela.  Efficient computation of higher-order cumulant tensors.  SIAM Journal on Scientific Computing, 40(3):A1590–A1610, 2018.
  * [18] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12(61):2121–2159, 2011.
  * [19] Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1:211–218, 1936.
  * [20] Lars Eldén and Berkant Savas.  A Newton-Grassmann method for computing the best multilinear rank-$(r_{1},r_{2},r_{3})$ approximation of a tensor.  SIAM Journal on Matrix Analysis and Applications, 31(2):248–271, 2009.
  * [21] Eugene F. Fama and Kenneth R. French.  Common risk factors in the returns on stocks and bonds.  Journal of Financial Economics, 33(1):3–56, 1993.
  * [22] Jerome H. Friedman.  Exploratory projection pursuit.  Journal of the American Statistical Association, 82(397):249–266, 1987.
  * [23] Rong Ge, Qingqing Huang, and Sham M. Kakade.  Learning mixtures of Gaussians in high dimensions.  In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing, pages 761–770, 2015.
  * [24] Xiurui Geng, Kang Sun, Luyan Ji, and Yongchao Zhao.  A high-order statistical tensor based algorithm for anomaly detection in hyperspectral imagery.  Scientific Reports, 4:6869, 11 2014.
  * [25] Wolfgang Hackbusch.  Numerical tensor calculus.  Acta Numerica, 23:651–742, 05 2014.
  * [26] Uwe Helmke, Knut Hüper, and Jochen Trumpf.  Newton’s method on Grassmann manifolds.  arXiv preprint arXiv:0709.2205, 2007.
  * [27] Amelia Henriksen and Rachel Ward.  AdaOja: Adaptive learning rates for streaming PCA.  arXiv preprint arXiv:1905.12115, 2019.
  * [28] Christopher J. Hillar and Lek-Heng Lim.  Most tensor problems are NP-hard.  Journal of the ACM (JACM), 60(6):1–39, 2013.
  * [29] Aapo Hyvärinen and Errki Oja.  Independent component analysis: algorithms and applications.  Neural Networks, 13(4):411–430, 2000.
  * [30] Mariya Ishteva, P.-A. Absil, and Paul Van Dooren.  Jacobi algorithm for the best low multilinear rank approximation of symmetric tensors.  SIAM Journal on Matrix Analysis and Applications, 34(2):651–672, 2013.
  * [31] Mariya Ishteva, P.-A. Absil, Sabine Van Huffel, and Lieven De Lathauwer.  Best low multilinear rank approximation of higher-order tensors, based on the Riemannian trust-region scheme.  SIAM Journal on Matrix Analysis and Applications, 32(1):115–135, 2011.
  * [32] Mariya Ishteva, Lieven De Lathauwer, P.-A Absil, and Sabine Huffel.  Differential-geometric Newton method for the best rank-($r_{1},r_{2},r_{3}$) approximation of tensors.  Numerical Algorithms, 51:179–194, 06 2009.
  * [33] Eric Jondeau and Michael Rockinger.  Optimal portfolio allocation under higher moments.  European Financial Management, 12:29–55, 01 2006.
  * [34] Tamara G. Kolda.  A counterexample to the possibility of an extension of the Eckart–Young low-rank approximation theorem for the orthogonal rank tensor decomposition.  SIAM Journal on Matrix Analysis and Applications, 24(3):762–767, 2003.
  * [35] Tamara G. Kolda and Brett W. Bader.  Tensor decompositions and applications.  SIAM Review, 51(3):455–500, 2009.
  * [36] Tamara G. Kolda, Brett W. Bader, and et al.  _Tensor Toolbox for MATLAB_ v. 3.2.1, 2021.
  * [37] Steven Krantz and Harold Parks.  A Primer of Real Analytic Functions.  Birkhauser, 2002.
  * [38] Jason D. Lee, Ioannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I. Jordan, and Benjamin Recht.  First-order methods almost always avoid strict saddle points.  Mathematical Programming, 176(1):311–337, 2019.
  * [39] John M. Lee.  Smooth Manifolds.  Springer, 2012.
  * [40] Jianze Li, Konstantin Usevich, and Pierre Comon.  Globally convergent Jacobi-type algorithms for simultaneous orthogonal symmetric tensor diagonalization.  SIAM Journal on Matrix Analysis and Applications, 39(1):1–22, 2018\.
  * [41] Yi Li, Huy L Nguyen, and David P Woodruff.  Turnstile streaming algorithms might as well be linear sketches.  In Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, pages 174–183, 2014.
  * [42] Osman Asif Malik and Stephen Becker.  Low-rank Tucker decomposition of large tensors using TensorSketch.  In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
  * [43] K. V. Mardia.  Measures of multivariate skewness and kurtosis with applications.  Biometrika, 57(3):519–530, 12 1970.
  * [44] Rachel Minster, Arvind K. Saibaba, and Misha E. Kilmer.  Randomized algorithms for low-rank tensor decompositions in the Tucker format.  SIAM Journal on Mathematics of Data Science, 2(1):189–215, 2020\.
  * [45] Ioannis Mitliagkas, Constantine Caramanis, and Prateek Jain.  Memory limited, streaming PCA.  In Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
  * [46] Ab Mooijaart.  Factor analysis for non-normal variables.  Psychometrika, 50:323–342, 1985.
  * [47] Jason Morton and Lek-Heng Lim.  Principal cumulant component analysis.  2009\.
  * [48] S. Muthukrishnan.  Data streams: Algorithms and applications.  Foundations and Trends in Theoretical Computer Science, 1(2):117–236, 08 2005.
  * [49] Erkki Oja.  Simplified neuron model as a principal component analyzer.  Journal of Mathematical Biology, 15:267–273, 1982.
  * [50] João M. Pereira, Joe Kileel, and Tamara G. Kolda.  Tensor moments of Gaussian mixture models: Theory and applications.  arXiv preprint arXiv:2202.06930, 2022.
  * [51] Antonio Plaza, Jon Atli Benediktsson, Joseph W. Boardman, Jason Brazile, Lorenzo Bruzzone, Gustavo Camps-Valls, Jocelyn Chanussot, Mathieu Fauvel, Paolo Gamba, Anthony Gualtieri, Mattia Marconcini, James C. Tilton, and Giovanna Trianni.  Recent advances in techniques for hyperspectral image processing.  Remote Sensing of Environment, 113:S110–S122, 2009.  Imaging Spectroscopy Special Issue.
  * [52] Phillip A. Regalia.  Monotonically convergent algorithms for symmetric tensor approximation.  Linear Algebra and its Applications, 438(2):875–890, 2013.
  * [53] Venkatraman Renganathan, Navid Hashemi, Justin Ruths, and Tyler H. Summers.  Higher-order moment-based anomaly detection.  IEEE Control Systems Letters, 6:211–216, 2022.
  * [54] Berkant Savas.  Best low rank tensor approximation, 2009.
  * [55] Berkant Savas and Lek-Heng Lim.  Quasi-Newton methods on Grassmannians and multilinear approximations of tensors.  SIAM Journal on Scientific Computing, 32(6):3352–3393, 2010.
  * [56] Reinhold Schneider and André Uschmajew.  Convergence results for projected line-search methods on varieties of low-rank matrices via Łojasiewicz inequality.  SIAM Journal on Optimization, 25(1):622–646, 2015.
  * [57] Samantha Sherman and Tamara G. Kolda.  Estimating higher-order moments using symmetric tensor decomposition.  SIAM Journal on Matrix Analysis and Applications, 41(3):1369–1387, 2020.
  * [58] Michael Shub.  Global Stability of Dynamical Systems.  Springer, 1987.
  * [59] Yiming Sun, Yang Guo, Charlene Luo, Joel Tropp, and Madeleine Udell.  Low-rank Tucker approximation of a tensor from streaming data.  SIAM Journal on Mathematics of Data Science, 2(4):1123–1150, 2020\.
  * [60] Joel Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher.  Streaming low-rank matrix approximation with an application to scientific simulation.  SIAM Journal on Scientific Computing, 41(4):A2430–A2463, 2019.
  * [61] Ledyard R. Tucker.  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31:279–311, 1966.
  * [62] Nick Vannieuwenhoven, Raf Vandebril, and Karl Meerbergen.  A new truncation strategy for the higher-order singular value decomposition.  SIAM Journal on Scientific Computing, 34:1027–1052, 04 2012.
  * [63] N. Vervliet, O. Debals, L. Sorber, M. Van Barel, and L. De Lathauwer.  _Tensorlab_ 3.0, 3 2016.
  * [64] Sebastian F. Walter, Lutz Lehmann, and Rene Lamour.  On evaluating higher-order derivatives of the QR decomposition of tall matrices with full column rank in forward and reverse mode algorithmic differentiation.  Optimization Methods and Software, 27(2):391–403, 2012.
  * [65] Puyudi Yang, Cho-Jui Hsieh, and Jane-Ling Wang.  History PCA: A new algorithm for streaming PCA.  arXiv preprint arXiv:1802.05447, 2018.
  * [66] Ralf Zimmermann.  Hermite interpolation and data processing errors on Riemannian matrix manifolds.  SIAM Journal on Scientific Computing, 42(5):A2593–A2619, 2020.

##  Appendix A Proof of Proposition 6.3

###### Proof.

We first compute the Riemannian gradient,
$\operatorname{grad}f(\mathbf{P})\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$.
It is the orthogonal projection of the Euclidean gradient $\nabla
f(\mathbf{P})\in\mathbb{R}^{n\times r}$ onto the tangent space at
$\mathbf{P}$. The projection
$\pi_{\mathbf{P}}:\mathbb{S}^{2}(\mathbb{R}^{n})\to\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
is computed as follows [26, Proposition 2.1]:

$\pi_{\mathbf{P}}(\mathbf{X})=\text{ad}^{2}_{\mathbf{P}}(\mathbf{X})\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\
\text{ad}_{\mathbf{P}}^{2}=\text{ad}_{\mathbf{P}}\circ\text{ad}_{\mathbf{P}}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
\text{ad}_{\mathbf{P}}(\mathbf{X})=\mathbf{P}\mathbf{X}-\mathbf{X}\mathbf{P}$

for $\mathbf{X}\in\mathbb{S}^{2}(\mathbb{R}^{n})$. Here, ad stands for the
adjoint operator.

The Euclidean gradient of $f$ 19 is

$\begin{array}[]{ll}\nabla
f(\mathbf{P})&=d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}+d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\\\
&=w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P}).\end{array}$ (43)

Hence the Riemannian gradient equals

$\begin{array}[]{ll}\operatorname{grad}f(\mathbf{P})&=\text{ad}^{2}_{\mathbf{P}}\left(\nabla
f(\mathbf{P})\right)\\\ &=\text{ad}_{\mathbf{P}}\left(\mathbf{P}\nabla
f(\mathbf{P})-\nabla f(\mathbf{P})\mathbf{P}\right)\\\
&=\mathbf{P}\left(\mathbf{P}\nabla f(\mathbf{P})-\nabla
f(\mathbf{P})\mathbf{P}\right)-\left(\mathbf{P}\nabla f(\mathbf{P})-\nabla
f(\mathbf{P})\mathbf{P}\right)\mathbf{P}\\\ &=\mathbf{P}\left(\nabla
f(\mathbf{P})\right)-2\mathbf{P}\left(\nabla
f(\mathbf{P})\right)\mathbf{P}+\left(\nabla
f(\mathbf{P})\right)\mathbf{P}\end{array}$

Substituting 43 into the above equation,

$\begin{array}[]{ll}\operatorname{grad}f(\mathbf{P})&=\mathbf{P}\left(w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})\right)-2\mathbf{P}\left(w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})\right)\mathbf{P}+\left(w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})\right)\mathbf{P}\\\
&=\mathbf{P}w(\mathbf{P})\left(\mathbf{I}_{n}-\mathbf{P}\right)+\left(\mathbf{I}_{n}-\mathbf{P}\right)w(\mathbf{P})\mathbf{P}=\operatorname{sym}\Big{(}2\left(\mathbf{I}_{n}-\mathbf{P}\right)w(\mathbf{P})\mathbf{P}\Big{)}.\end{array}$
(44)

Next, we relate the gradients of $f(\mathbf{P})$ and $F(\mathbf{Q})$. The
Euclidean gradient of $F$ is

$\begin{array}[]{ll}\nabla
F(\mathbf{Q})&=2d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{Q},\dots,\mathbf{Q}\right)\right\rangle_{-1}\\\
&=2d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{Q}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\mathbf{Q}\\\
&=2d\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\mathbf{Q}=2w(\mathbf{P})\mathbf{Q}.\end{array}$
(45)

Therefore,

$\left(\mathbf{I}_{n}-\mathbf{Q}\mathbf{Q}^{\top}\right)\nabla
F(\mathbf{Q})\mathbf{Q}^{\top}=2\left(\mathbf{I}_{n}-\mathbf{P}\right)w(\mathbf{P})\mathbf{Q}\mathbf{Q}^{\top}=2\left(\mathbf{I}_{n}-\mathbf{P}\right)w(\mathbf{P})\mathbf{P}.$

By 44, we have

$\operatorname{grad}f(\mathbf{P})=\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}\mathbf{Q}^{\top}\right)\nabla
F(\mathbf{Q})\mathbf{Q}^{\top}\Big{)}.$

The proof of Proposition 6.3 is complete. ∎

##  Appendix B Proof of Lemma 6.4

###### Proof.

We derive the update on $\mathbf{P}_{t}$ 18 based on the iterate
$\mathbf{Q}_{t}$ 17. For a step $t\geq 0$, decompose the update direction
$\gamma_{t}\nabla F(\mathbf{Q}_{t})$ into two parts respectively in
$\text{colspan}(\mathbf{Q}_{t})$ and $\text{colspan}(\mathbf{Q}_{t})^{\perp}$,

$\gamma_{t}\nabla
F(\mathbf{Q}_{t})=\gamma_{t}\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\nabla
F(\mathbf{Q}_{t})+\gamma_{t}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t}).$

We approximate
$\mathbf{Q}_{t+1}=\framebox{Q}\text{R}\left(\mathbf{Q}_{t}+\gamma_{t}\nabla
F(\mathbf{Q}_{t})\right)$ by the first-order Taylor polynomial of QR centered
at $\mathbf{Q}_{t}+\gamma_{t}\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\nabla
F(\mathbf{Q}_{t})$:

$\begin{array}[]{ll}\mathbf{Q}_{t+1}&=\framebox{Q}\text{R}\Big{(}\left(\mathbf{Q}_{t}+\gamma_{t}\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\nabla
F(\mathbf{Q}_{t})\right)+\gamma_{t}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\Big{)}\\\
&=\widetilde{\mathbf{Q}_{t}}+\text{D}\framebox{Q}\text{R}(\widetilde{\mathbf{Q}_{t}})[\gamma_{t}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})]+\mathcal{O}\left(\gamma_{t}^{2}\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\|^{2}\right),\end{array}$

where

$\widetilde{\mathbf{Q}_{t}}\mathbf{R}=\mathbf{Q}_{t}+\gamma_{t}\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\nabla
F(\mathbf{Q}_{t})\in\mathbb{R}^{n\times r}$ (46)

denotes the QR factorization. (Note that
$\framebox{Q}\text{R}:\operatorname{St}(n,r)\to\operatorname{St}(n,r)$ is
$C^{\infty}$ at full-rank inputs [64].)

By the formula for the derivative of the QR factorization given in [66,
Algorithm 1] (originally [64]),

$\begin{array}[]{ll}\mathbf{Q}_{t+1}&=\widetilde{\mathbf{Q}_{t}}+\gamma_{t}\left(\mathbf{I}_{n}-\widetilde{\mathbf{Q}_{t}}\widetilde{\mathbf{Q}_{t}}^{\top}\right)\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\\\ &\leavevmode\nobreak\ \leavevmode\nobreak\
+\widetilde{\mathbf{Q}_{t}}\,s\left(\gamma_{t}\widetilde{\mathbf{Q}_{t}}^{\top}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\right)+\mathcal{O}\left(\gamma_{t}^{2}\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\|^{2}\right)\\\
&=\widetilde{\mathbf{Q}_{t}}+\gamma_{t}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}+\mathcal{O}\left(\gamma_{t}^{2}\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\|^{2}\right).\end{array}$ (47)

To see the second equality, note
$\mathbf{I}_{n}-\widetilde{\mathbf{Q}_{t}}\widetilde{\mathbf{Q}_{t}}^{\top}\\!=\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}$,
because $\widetilde{\mathbf{Q}_{t}}$ and $\mathbf{Q}_{t}$ are both the
orthonormal bases for $\text{colspan}(\mathbf{Q}_{t}).$ Moreover, given the
linear operator $s:\mathbb{R}^{r\times r}\to\mathbb{R}^{r\times r}$,444The
operator $s$ is a composition of element-wise multiplication and skew-
symmetrization. See [66, Algorithm 1] for details. the term
$s\left(\gamma_{t}\widetilde{\mathbf{Q}_{t}}^{\top}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\right)=\bf 0$ as
$\widetilde{\mathbf{Q}_{t}}^{\top}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)=\bf
0.$

From 47,

$\begin{array}[]{ll}\mathbf{P}_{t+1}-\mathbf{P}_{t}=\mathbf{Q}_{t+1}\mathbf{Q}_{t+1}^{\top}-\widetilde{\mathbf{Q}_{t}}\widetilde{\mathbf{Q}_{t}}^{\top}\\!\\!\\!\\!\\!&=2\gamma_{t}\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}\Big{)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
+\mathcal{O}\left(\gamma_{t}^{2}\left\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\right\|^{2}\right)\\\ &\leavevmode\nobreak\
\leavevmode\nobreak\
+\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{Q}_{t}^{\top}\Big{)}\right\|^{2}\right).\end{array}$
(48)

Here the second line refers to the symmetric term

$\gamma_{t}^{2}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\mathbf{R}^{-\top}\nabla
F(\mathbf{Q}_{t})^{\top}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right),$

whose norm is indeed less than
$\gamma_{t}^{2}\left\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\right\|^{2}$.

We transform the right-hand side of 48 from an expression in terms of
$\mathbf{Q}$ to one with $\mathbf{P}$.

For the first component of 48, we separate
$\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}$ into two
parts:

$\texttt{Factor\leavevmode\nobreak\
1}:=\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t}){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathbf{Q}_{t}^{\top}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\texttt{Factor\leavevmode\nobreak\
2}:={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathbf{Q}_{t}}\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top},$

using $\mathbf{Q}_{t}^{\top}\mathbf{Q}_{t}=\mathbf{I}_{r}$. By 45,

$\texttt{Factor\leavevmode\nobreak\
1}=2(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}=2(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\mathbf{P}_{t}.$
(49)

Also,

$\texttt{Factor\leavevmode\nobreak\
2}=\mathbf{Q}_{t}\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}=(\mathbf{Q}_{t}^{\top})^{\dagger}\mathbf{R}^{-1}(\widetilde{\mathbf{Q}_{t}})^{\dagger}=\left(\widetilde{\mathbf{Q}_{t}}\mathbf{R}\mathbf{Q}_{t}^{\top}\right)^{\dagger}.$
(50)

by the fact that pseudoinversion satisfies
$\mathbf{X}^{\dagger}\mathbf{Y}^{\dagger}=(\mathbf{Y}\mathbf{X})^{\dagger}$ if
$\mathbf{X}$ and $\mathbf{Y}$ respectively have orthonormal rows and columns.
Inserting 46 gives

$\texttt{Factor\leavevmode\nobreak\
2}=\left(\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}+\gamma_{t}\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\nabla
F(\mathbf{Q}_{t})\mathbf{Q}_{t}^{\top}\right)^{\dagger}=\left(\mathbf{P}_{t}+2\gamma_{t}\mathbf{P}_{t}w(\mathbf{P}_{t})\mathbf{P}_{t}\right)^{\dagger}=\alpha_{t}(\mathbf{P}_{t})^{\dagger},$

where the last equality is by 45 and the definition of $\alpha_{t}$ 25.
Moreover,

$\texttt{Factor\leavevmode\nobreak\
1}^{\top}\times\texttt{Factor\leavevmode\nobreak\ 2}=\mathbf{Q}_{t}\nabla
F(\mathbf{Q}_{t})^{\top}\underbrace{\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\mathbf{Q}_{t}}_{\bf
0}\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}=\bf 0.$

Combining 49 and 50,
$\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}$ equals:

$\begin{array}[]{l}\texttt{Factor\leavevmode\nobreak\
1}\times\texttt{Factor\leavevmode\nobreak\
2}=2(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\mathbf{P}_{t}\alpha_{t}(\mathbf{P}_{t})^{\dagger}\\\
=(\texttt{Factor\leavevmode\nobreak\ 1}+\texttt{Factor\leavevmode\nobreak\
1}^{\top})\times\texttt{Factor\leavevmode\nobreak\
2}=2\operatorname{sym}\Big{(}2(\mathbf{I}_{n}-\mathbf{P}_{t})w(\mathbf{P}_{t})\mathbf{P}_{t}\Big{)}\alpha_{t}(\mathbf{P}_{t})^{\dagger}\\\
=2\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}.\end{array}$
(51)

The formula of $\operatorname{grad}f$ 44 is used in the last equality.
Therefore, the first component of 48 reads

$2\gamma_{t}\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}\Big{)}=4\gamma_{t}\operatorname{sym}\Big{(}\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\Big{)}.$
(52)

As for second component of 48, note

$\begin{array}[]{ll}\left\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\right\|^{2}&=\left\|\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{R}^{-1}\widetilde{\mathbf{Q}_{t}}^{\top}\right\|^{2}\\\
&=4\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|^{2}.\end{array}$

The orthogonality of $\widetilde{\mathbf{Q}_{t}}$, which preserves norms,
gives the first equality. The second is from 51.

The third component of 48 is dealt with directly by 24,

$\begin{array}[]{ll}\left\|\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\mathbf{Q}_{t}\mathbf{Q}_{t}^{\top}\right)\nabla
F(\mathbf{Q}_{t})\mathbf{Q}_{t}^{\top}\Big{)}\right\|^{2}&=\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}.\end{array}$

Thus, 48 consists of the main component 52 and the residuals
$\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\alpha_{t}(\mathbf{P}_{t})^{\dagger}\right\|^{2}\right)$,
$\mathcal{O}\left(\gamma_{t}^{2}\left\|\operatorname{grad}f(\mathbf{P}_{t})\right\|^{2}\right)$
as needed in 26. This completes the proof of Lemma 6.4. ∎

##  Appendix C Proof of Lemma 6.5

###### Proof.

Consider the step size $\gamma_{t}>0$ and projector
$\mathbf{P}\in\operatorname{Gr}(n,r).$ Definition 25 reads

$\alpha_{t}(\mathbf{P})=\mathbf{P}+2\gamma_{t}\mathbf{P}w(\mathbf{P})\mathbf{P}=\mathbf{P}\left(\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P})\right)\mathbf{P}.$
(53)

So
$\text{colspan}(\alpha_{t}(\mathbf{P}))\subseteq\text{colspan}(\mathbf{P})$.
Also for any vector $\mathbf{x}\in\text{colspan}(\mathbf{P})$,

$\mathbf{x}^{\top}\alpha_{t}(\mathbf{P})\mathbf{x}=\mathbf{x}^{\top}\mathbf{P}(\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P}))\mathbf{P}\mathbf{x}=\mathbf{x}^{\top}(\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P}))\mathbf{x}\geq\|\mathbf{x}\|^{2},$

as $\mathbf{P}\mathbf{x}=\mathbf{x}$ and
$\mathbf{I}_{n}+2\gamma_{t}w(\mathbf{P})\succeq\mathbf{I}_{n}.$ It follows
$\text{colspan}(\alpha_{t}(\mathbf{P}))=\text{colspan}(\mathbf{P})$, and
$\alpha_{t}(\mathbf{P})$ has only $r$ non-zero singular values, each greater
than $1$. By properties of pseudoinverses,

$\alpha_{t}(\mathbf{P})^{\dagger}\alpha_{t}(\mathbf{P})=\textbf{Proj}_{\text{range}(\alpha_{t}(\mathbf{P}))}=\mathbf{P},$
(54)

and $\alpha_{t}(\mathbf{P})^{\dagger}$ has only $r$ non-zero singular values,
each less than 1. Hence

$\|\alpha_{t}(\mathbf{P})^{\dagger}\|\leq\sqrt{r}.$ (55)

Constrain the step size by

$0<\gamma_{t}\leq\Gamma_{1}:=\frac{1}{4\sqrt{r}\sup_{\mathbf{P}\in\operatorname{Gr}(n,r)}\|\mathbf{P}w(\mathbf{P})\mathbf{P}\|},$
(56)

where
$\sup_{\mathbf{P}\in\operatorname{Gr}(n,r)}\|\mathbf{P}w(\mathbf{P})\mathbf{P}\|$
is finite because the function
$\mathbf{P}\mapsto\mathbf{P}w(\mathbf{P})\mathbf{P}$ is continuous and
$\operatorname{Gr}(n,r)$ is compact. (We exclude the trivial case
$\mathbf{P}w(\mathbf{P})\mathbf{P}\equiv 0$ when
$\boldsymbol{\mathscr{X}}=0$.) From 54,

$\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}=\alpha_{t}(\mathbf{P})^{\dagger}\alpha_{t}(\mathbf{P})-\alpha_{t}(\mathbf{P})^{\dagger}\mathbf{P}=-\alpha_{t}(\mathbf{P})^{\dagger}(\mathbf{P}-\alpha_{t}(\mathbf{P})).$
(57)

Then by 55 and 56,

$\begin{array}[]{ll}\|\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}\|&=\|-\alpha_{t}(\mathbf{P})^{\dagger}\left(\mathbf{P}-\alpha_{t}(\mathbf{P})\right)\|\\\
&\leq\|\alpha_{t}(\mathbf{P})^{\dagger}\|\|\mathbf{P}-\alpha_{t}(\mathbf{P})\|\\\
&\displaystyle\leq
2\gamma_{t}\sqrt{r}\|\mathbf{P}w(\mathbf{P})\mathbf{P}\|\leq\frac{1}{2}.\end{array}$
(58)

We now compare $\|\operatorname{grad}f(\mathbf{P})\|$ and
$\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|$.
By 24,

$\operatorname{grad}f(\mathbf{P})\mathbf{P}=(\mathbf{I}_{n}-\mathbf{P})w(\mathbf{P})\mathbf{P}.$

Inserting the above equation into 24,

$\operatorname{grad}f(\mathbf{P})=\operatorname{sym}\Big{(}2(\mathbf{I}_{n}-\mathbf{P})w(\mathbf{P})\mathbf{P}\Big{)}=\operatorname{sym}\Big{(}2\operatorname{grad}f(\mathbf{P})\mathbf{P}\Big{)}.$

Since

$\left\langle\operatorname{grad}f(\mathbf{P})\mathbf{P},\mathbf{P}\operatorname{grad}f(\mathbf{P})\right\rangle=\text{trace}\left((\mathbf{I}_{n}-\mathbf{P})w(\mathbf{P})\underbrace{\mathbf{P}(\mathbf{I}_{n}-\mathbf{P})}_{\bf
0}w(\mathbf{P})\mathbf{P}\right)=0,$

we have

$\left\|\operatorname{grad}f(\mathbf{P})\right\|=\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\mathbf{P}\right\|.$
(59)

Then due to triangle inequality and 57,

$\begin{array}[]{ll}\|\operatorname{grad}f(\mathbf{P})\|&\leq\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|+\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\left(\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}\right)\right\|\\\
&\displaystyle=\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|+\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\left(\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}\right)\right\|\\\
&\displaystyle\leq\sqrt{2}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|\left(1+\left\|\mathbf{P}-\alpha_{t}(\mathbf{P})^{\dagger}\right\|\right)=\frac{3\sqrt{2}}{2}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|.\end{array}$

Inserting 58 shows the last inequality. Thus,

$\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|\geq\frac{\sqrt{2}}{3}\left\|\operatorname{grad}f(\mathbf{P})\right\|.$

On the other hand, as
$\alpha_{t}(\mathbf{P})^{\dagger}=\mathbf{P}\alpha_{t}(\mathbf{P})^{\dagger},$
it holds

$\begin{array}[]{ll}\left\|\operatorname{grad}f(\mathbf{P})\alpha_{t}(\mathbf{P})^{\dagger}\right\|&\leq\left\|\operatorname{grad}f(\mathbf{P})\mathbf{P}\right\|\left\|\alpha_{t}(\mathbf{P})^{\dagger}\right\|\\\
&\leq\displaystyle\frac{\sqrt{2r}}{2}\left\|\operatorname{grad}f(\mathbf{P})\right\|.\end{array}$

The results 59 and 55 give the last derivation. Combining the above two
inequalities completes the proof of Lemma 6.5. ∎

##  Appendix D Proof of Proposition 6.6

###### Proof.

Let $\mathbf{P}\in\operatorname{Gr}(n,r)$ and
$\Delta\mathbf{P}\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$. We first
show the computation of the Riemannian Hessian. Applying [26, Equation
(2.109)],

$\begin{array}[]{ll}\operatorname{Hess}f(\mathbf{P})[\Delta\mathbf{P}]&=\text{ad}_{\mathbf{P}}^{2}\left(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}]\right)-\text{ad}_{\mathbf{P}}\text{ad}_{\nabla
f(\mathbf{P})}(\Delta\mathbf{P})\\\
&=\text{ad}_{\mathbf{P}}\Big{(}\mathbf{P}\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}]-(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}])\mathbf{P}-\nabla
f(\mathbf{P})\Delta\mathbf{P}+\Delta\mathbf{P}\nabla f(\mathbf{P})\Big{)}\\\
&=\mathbf{P}(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}])(\mathbf{I}_{n}-\mathbf{P})+(\mathbf{I}_{n}-\mathbf{P})(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}])\mathbf{P}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ -\mathbf{P}\nabla
f(\mathbf{P})\Delta\mathbf{P}+\mathbf{P}\Delta\mathbf{P}\nabla
f(\mathbf{P})+\nabla
f(\mathbf{P})(\Delta\mathbf{P})\mathbf{P}-\Delta\mathbf{P}\nabla
f(\mathbf{P})\mathbf{P}\\\
&=2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\mathbf{P})(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}])\mathbf{P}+\nabla
f(\mathbf{P})(\Delta\mathbf{P})\mathbf{P}-\Delta\mathbf{P}\nabla
f(\mathbf{P})\mathbf{P}\Big{)}.\end{array}$ (60)

Recall 43 that $\nabla
f(\mathbf{P})=w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})$. We
differentiate $\nabla f(\mathbf{P})$ to obtain the Euclidean Hessian of $f$:

$\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}]=\left(\text{D}w(\mathbf{P})[\Delta\mathbf{P}]\right)\mathbf{P}+w(\mathbf{P})\Delta\mathbf{P}+\Delta\mathbf{P}w(\mathbf{P})+\mathbf{P}\left(\text{D}w(\mathbf{P})[\Delta\mathbf{P}]\right),$
(61)

where the differential is

$\begin{array}[]{ll}\text{D}w(\mathbf{P})[\Delta\mathbf{P}]&=2d(d-1)\operatorname{sym}\Big{(}\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\Delta\mathbf{P},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{P},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\Big{)}\\\
&=d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},(\Delta\mathbf{P})\mathbf{P}+\mathbf{P}\Delta\mathbf{P},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}\\\
&=d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\Delta\mathbf{P},\mathbf{P},\dots,\mathbf{P}\right),\boldsymbol{\mathscr{X}}\cdot\left(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{P},\dots,\mathbf{P}\right)\right\rangle_{-1}.\end{array}$
(62)

The last equality follows from
$(\Delta\mathbf{P})\mathbf{P}+\mathbf{P}\Delta\mathbf{P}=\Delta\mathbf{P}$ in
21.

For the Riemannian Hessian 60, we first consider the term involving
$\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}]$:

$\begin{array}[]{l}(\mathbf{I}_{n}-\mathbf{P})(\nabla^{2}f(\mathbf{P})[\Delta\mathbf{P}])\mathbf{P}\\\
=(\mathbf{I}_{n}-\mathbf{P})\left(\text{D}w(\mathbf{P})[\Delta\mathbf{P}]\right)\mathbf{P}+(\mathbf{I}_{n}-\mathbf{P})w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}+(\mathbf{I}_{n}-\mathbf{P})\Delta\mathbf{P}w(\mathbf{P})\mathbf{P}\\\
=(\mathbf{I}_{n}-\mathbf{P})v(\mathbf{P})[\Delta\mathbf{P}]+(\mathbf{I}_{n}-\mathbf{P})w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}+(\mathbf{I}_{n}-\mathbf{P})\Delta\mathbf{P}w(\mathbf{P})\mathbf{P}.\end{array}$
(63)

The last equality is due to
$v(\mathbf{P})[\Delta\mathbf{P}]=\left(\text{D}w(\mathbf{P})[\Delta\mathbf{P}]\right)\mathbf{P}$
by 37 and 62.

We then compute the remaining terms in 60, which involve $\nabla
f(\mathbf{P})$:

$\begin{array}[]{l}\nabla
f(\mathbf{P})(\Delta\mathbf{P})\mathbf{P}-\Delta\mathbf{P}\nabla
f(\mathbf{P})\mathbf{P}\\\
=\left(w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})\right)\left(\Delta\mathbf{P}\right)\mathbf{P}-\Delta\mathbf{P}\left(w(\mathbf{P})\mathbf{P}+\mathbf{P}w(\mathbf{P})\right)\mathbf{P}\\\
=\mathbf{P}w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}-\Delta\mathbf{P}w(\mathbf{P})\mathbf{P}-\left(\Delta\mathbf{P}\right)\mathbf{P}w(\mathbf{P})\mathbf{P}.\end{array}$
(64)

In the last equality, we eliminate the term that includes
$\mathbf{P}(\Delta\mathbf{P})\mathbf{P}$, since this is $\bf 0$ due to 21.

We add 63 and 64 to obtain

$\begin{array}[]{rl}\operatorname{Hess}f(\mathbf{P})[\Delta\mathbf{P}]=&2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\mathbf{P})v(\mathbf{P})[\Delta\mathbf{P}]+(\mathbf{I}_{n}-\mathbf{P}+\mathbf{P})w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}\\\
&+(\mathbf{I}_{n}-\mathbf{P}-\mathbf{I}_{n})\left(\Delta\mathbf{P}\right)w(\mathbf{P})\mathbf{P}-\left(\Delta\mathbf{P}\right)\mathbf{P}w(\mathbf{P})\mathbf{P}\Big{)}\\\
=&2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\mathbf{P})v(\mathbf{P})[\Delta\mathbf{P}]+w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}-(\mathbf{P}\Delta\mathbf{P}+\left(\Delta\mathbf{P}\right)\mathbf{P})w(\mathbf{P})\mathbf{P}\Big{)}\\\
=&2\operatorname{sym}\Big{(}(\mathbf{I}_{n}-\mathbf{P})v(\mathbf{P})[\Delta\mathbf{P}]+w(\mathbf{P})\left(\Delta\mathbf{P}\right)\mathbf{P}-\Delta\mathbf{P}w(\mathbf{P})\mathbf{P}\Big{)}.\end{array}$

This is the formula for the Riemannian Hessian in 38.

Next, we compute the differential
$\text{D}\phi(\mathbf{P})[\Delta\mathbf{P}].$ This is the term that is linear
in $\Delta\mathbf{P}$ in the Taylor expansion of the difference:

$\phi(\mathbf{P}+\Delta\mathbf{P})-\phi(\mathbf{P})=\Phi(\mathbf{Q}+\Delta\mathbf{Q})\Phi(\mathbf{Q}+\Delta\mathbf{Q})^{\top}-\Phi(\mathbf{Q})\Phi(\mathbf{Q})^{\top}.$
(65)

So we start with $\Phi(\mathbf{Q}+\Delta\mathbf{Q})$ and $\Phi(\mathbf{Q})$.
From the first-order Taylor expansion of QR,

$\begin{array}[]{ll}\Phi(\mathbf{Q}+\Delta\mathbf{Q})&=\framebox{Q}\text{R}\left(\mathbf{Q}+\Delta\mathbf{Q}+\gamma\nabla
F(\mathbf{Q}+\Delta\mathbf{Q})\right)\\\
&=\framebox{Q}\text{R}\Big{(}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)+\Delta\mathbf{Q}+\gamma\left(\nabla
F(\mathbf{Q}+\Delta\mathbf{Q})-\nabla F(\mathbf{Q})\right)\Big{)}\\\
&=\framebox{Q}\text{R}\Big{(}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)+\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]+\mathcal{O}\left(\|\Delta\mathbf{Q}\|^{2}\right)\Big{)}\\\
&=\framebox{Q}\text{R}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)+\text{D}\framebox{Q}\text{R}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)\left[\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]+\mathcal{O}\left(\|\Delta\mathbf{Q}\|^{2}\right)\right].\end{array}$

Write $\Phi(\mathbf{Q})\mathbf{R}$ for the QR decomposition of
$\mathbf{Q}+\gamma\nabla F(\mathbf{Q})$. We again apply the differentiation
formula for QR in [66, Algorithm 1] and get

$\begin{array}[]{ll}\Phi(\mathbf{Q}+\Delta\mathbf{Q})&=\Phi(\mathbf{Q})+\left(\mathbf{I}_{n}-\Phi(\mathbf{Q})\Phi(\mathbf{Q})^{\top}\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right)\mathbf{R}^{-1}+\Phi(\mathbf{Q})\mathbf{S}+\mathcal{O}\left(\|\Delta\mathbf{Q}\|^{2}\right).\end{array}$
(66)

Here $\mathbf{S}\in\mathbb{R}^{r\times r}$ is a skew-symmetric matrix that is
linear in $\Delta\mathbf{Q}$. The Euclidean gradient $\nabla F(\mathbf{Q})$ is
given by 45, and the Euclidean Hessian of $F$ is calculated as

$\begin{array}[]{ll}\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]&=2d\left\langle\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q}),\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q})\right\rangle_{-1}\Delta\mathbf{Q}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
+2d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},(\Delta\mathbf{Q})\mathbf{Q}+\mathbf{Q}\Delta\mathbf{Q},\mathbf{Q},\dots,\mathbf{Q}),\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q})\right\rangle_{-1}\mathbf{Q}\\\
&=2w(\mathbf{P})\Delta\mathbf{Q}+2d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\Delta\mathbf{P},\mathbf{Q},\dots,\mathbf{Q}),\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q})\right\rangle_{-1}\mathbf{Q},\end{array}$
(67)

due to 23 and
$\Delta\mathbf{P}=\left(\Delta\mathbf{Q}\right)\mathbf{Q}^{\top}+\mathbf{Q}\left(\Delta\mathbf{Q}\right)^{\top}.$

We substitute 66 into 65. We only keep the terms that are linear in
$\Delta\mathbf{Q}$ and get

$\begin{array}[]{ll}\text{D}\phi(\mathbf{P})[\Delta\mathbf{P}]&=2\,\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\Phi(\mathbf{Q})\Phi(\mathbf{Q})^{\top}\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right)\mathbf{R}^{-1}\Phi(\mathbf{Q})^{\top}\Big{)}\\\
&=2\,\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right)\left(\mathbf{R}^{-1}\Phi(\mathbf{Q})^{\top}\right)\Big{)}\\\
&=2\,\operatorname{sym}\Big{(}\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right)\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)^{\dagger}\Big{)}.\end{array}$ (68)

For the first equality, we use
$\Phi(\mathbf{Q})\mathbf{S}\Phi(\mathbf{Q})^{\top}+\Phi(\mathbf{Q})\mathbf{S}^{\top}\Phi(\mathbf{Q})^{\top}=\bf
0$ and the fact $\mathbf{S}$ is linear in $\Delta\mathbf{Q}$

It remains to transform the differential expression 68 in
$\mathbf{Q},\Delta\mathbf{Q}$ to be in terms of $\mathbf{P},\Delta\mathbf{P}$.
To do this, we split 68 into two factors:

$\texttt{Factor\leavevmode\nobreak\
1}:=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathbf{Q}^{\top}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\texttt{Factor\leavevmode\nobreak\
2}:={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathbf{Q}}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)^{\dagger}.$

We first tackle Factor 2. As $\mathbf{Q}^{\top}$ and $\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})$ respectively have orthonormal rows and columns, we have

$\mathbf{Q}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)^{\dagger}=\left(\mathbf{Q}^{\top}\right)^{\dagger}\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)^{\dagger}=\left((\mathbf{Q}+\gamma\nabla
F(\mathbf{Q}))\mathbf{Q}^{\top}\right)^{\dagger}.$

Inserting 45 into the above equation,

$\begin{array}[]{ll}\texttt{Factor\leavevmode\nobreak\
2}&=\left(\left(\mathbf{Q}+2\gamma
w(\mathbf{P})\mathbf{Q}\right)\mathbf{Q}^{\top}\right)^{{\dagger}}=\left(\mathbf{P}+2\gamma
w(\mathbf{P})\mathbf{P}\right)^{{\dagger}}.\end{array}$ (69)

Next, we consider Factor 1. By 67, it equals

$\begin{array}[]{l}\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\Delta\mathbf{Q}+\gamma\nabla^{2}F(\mathbf{Q})[\Delta\mathbf{Q}]\right)\mathbf{Q}^{\top}\\\
=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\right)(\Delta\mathbf{Q})\mathbf{Q}^{\top}\\\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(2\gamma
d(d-1)\left\langle\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\Delta\mathbf{P},\mathbf{Q},\dots,\mathbf{Q}),\boldsymbol{\mathscr{X}}\cdot(\mathbf{I}_{n},\mathbf{I}_{n},\mathbf{Q},\dots,\mathbf{Q})\right\rangle_{-1}\mathbf{Q}\mathbf{Q}^{\top}\right)\\\
=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\left(\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\right)(\Delta\mathbf{P}-\mathbf{Q}(\Delta\mathbf{Q})^{\top})+2\gamma
v(\mathbf{P})[\Delta\mathbf{P}]\right).\end{array}$ (70)

The last equality is from the definition of $v(\mathbf{P})[\Delta\mathbf{P}]$
37.

However, note that

$\begin{array}[]{ll}\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\right)\mathbf{Q}(\Delta\mathbf{Q})^{\top}&=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\mathbf{Q}+2\gamma
w(\mathbf{P})\mathbf{Q}\right)(\Delta\mathbf{Q})^{\top}\\\
&=\left(\mathbf{I}_{n}-\Phi(\mathbf{Q})\Phi(\mathbf{Q})^{\top}\right)\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)(\Delta\mathbf{Q})^{\top}=\bf 0.\end{array}$

Here, we replace $2w(\mathbf{P})\mathbf{Q}$ by $\nabla F(\mathbf{Q})$ by 45.
Then from $\text{colspan}(\mathbf{Q}+\gamma\nabla F(\mathbf{Q}))$ equals
$\text{colspan}(\Phi(\mathbf{Q}))$ by 17, it holds
$\left(\mathbf{I}_{n}-\Phi(\mathbf{Q})\Phi(\mathbf{Q})^{\top}\right)\left(\mathbf{Q}+\gamma\nabla
F(\mathbf{Q})\right)=\bf 0$.

Returning to 70,

$\begin{array}[]{ll}\texttt{Factor\leavevmode\nobreak\
1}&=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\left(\mathbf{I}_{n}+2\gamma
w(\mathbf{P})\right)(\Delta\mathbf{P})+2\gamma
v(\mathbf{P})[\Delta\mathbf{P}]\right)\\\
&=\left(\mathbf{I}_{n}-\phi(\mathbf{P})\right)\left(\Delta\mathbf{P}+2\gamma
w(\mathbf{P})(\Delta\mathbf{P})+2\gamma
v(\mathbf{P})[\Delta\mathbf{P}]\right).\end{array}$ (71)

Putting 68, 69 and 71 together, we obtain the formula for the differential
$\text{D}\phi(\mathbf{P})[\Delta\mathbf{P}]$ 39. The proof of Proposition 6.6
is complete. ∎

##  Appendix E Proof of Lemma 6.7

###### Proof.

Let $\mathbf{P}\in\operatorname{Gr}(n,r)$ be a first-order critical point of
$f$ 19. The Riemannian Hessian
$\operatorname{Hess}f(\mathbf{P}):\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)\rightarrow\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
is diagonalizable because it is self-adjoint. The dimension of the tangent
space $\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$ is
$\text{dim}(\operatorname{Gr}(n,r))=r(n-r)$. Let
$\Delta\mathbf{P}_{1},\dots,\Delta\mathbf{P}_{r(n-r)}\in\text{T}_{\mathbf{P}}\operatorname{Gr}(n,r)$
be a basis of eigenmatrices of $\operatorname{Hess}f(\mathbf{P})$ with
corresponding eigenvalues $\lambda_{1},\dots,\lambda_{r(n-r)},$ i.e.

$\operatorname{Hess}f(\mathbf{P})[\Delta\mathbf{P}_{i}]=\lambda_{i}\Delta\mathbf{P}_{i}$

for $i\in[r(n-r)]$. Assume there is at least one strictly positive eigenvalue
of $\operatorname{Hess}f(\mathbf{P})$, w.l.o.g. $\lambda_{1}>0$.