# On the Multiway Principal Component Analysis∗

Jialin Ouyang and Ming Yuan
Department of Statistics
Columbia University

(())

###### Abstract

Multiway data are becoming more and more common. While there are many
approaches to extending principal component analysis (PCA) from usual data
matrices to multiway arrays, their conceptual differences from the usual PCA,
and the methodological implications of such differences remain largely
unknown. This work aims to specifically address these questions. In
particular, we clarify the subtle difference between PCA and singular value
decomposition (SVD) for multiway data, and show that multiway principal
components (PCs) can be estimated reliably in absence of the eigengaps
required by the usual PCA, and in general much more efficiently than the usual
PCs. Furthermore, the sample multiway PCs are asymptotically independent and
hence allow for separate and more accurate inferences about the population
PCs. The practical merits of multiway PCA are further demonstrated through
numerical, both simulated and real data, examples.

11footnotetext:  This research was supported in part by NSF Grants DMS-2015285
and DMS-2052955.

##  1 Introduction

More and more often in practice, we need to deal with data of rich and complex
structures that are more appropriately organized as multiway arrays rather
than the usual data matrices. Examples of such multiway data are ubiquitous in
many fields such as chemometrics, economics, psychometrics, and signal
processing among others (see, e.g., Kroonenberg, 2008; Anandkumar et al.,
2014; Zhang and Xia, 2018; Chen et al., 2020a, b; Han et al., 2020; Xia et
al., 2020; Bi et al., 2021; Chen et al., 2021; Han et al., 2022). In this
paper, we investigate the methodological implications and statistical
properties of principal component analysis (PCA) for this type of data and
pinpoint the benefits and challenges of doing so.

PCA is among the most popular statistical methods for multivariate data
analysis when data are organized as matrices. See, e.g., Anderson (1984);
Jolliffe (2002). With each column vector of a data matrix as an observation,
PCA seeks orthogonal linear transformations of these vectors into a new
coordinate system so that the variance of each coordinate is maximized
successively. It allows us to represent most of the variation in the data by a
small number of coordinates and therefore can guide us in reducing the
dimensionality. As such, PCA often serves as a critical first step to capture
the essential features in a dataset for many downstream analyses and is widely
used in many scientific and engineering fields. Moving beyond matrices, for
multiway data, each observation itself forms a matrix or more generally a
multiway array. For example, when repeated measurements are made across
different combinations of location and time, each observation can be more
naturally organized as a matrix with each row corresponding to a certain
location and each column a time point. To apply PCA to this type of data, it
is tempting to neglect the multiway nature of the observations and treat each
observation as a vector nonetheless, a practice often referred to as
_stringing_. However, as observed in numerous practical applications,
appropriately accounting for the additional structure when applying PCA can
greatly enhance interpretability and improve efficiency. See, e.g.,
Kroonenberg (2008).

There is a long and illustrious history of developing suitable methods for
such a purpose and it can be traced back at least to the pioneering work of
Tucker, Harshman, and Carroll in the 1960s. Since then, numerous approaches
have also been developed. Examples include Kroonenberg and De Leeuw (1980); De
Lathauwer et al. (2000); Vasilescu and Terzopoulos (2002); Yang et al. (2004);
Kong et al. (2005); Zhang and Zhou (2005); Lu et al. (2006, 2008); Li et al.
(2010); Liu et al. (2017); Taguchi (2018) among many others. See, e.g., Lu et
al. (2011); Cichocki et al. (2015) for recent surveys of existing techniques.
Most of these developments are outside the mainstream statistics literature
and often with a strong algorithmic flavor and exploratory data analysis
focus. These approaches are intuitive and often yield more interpretable
insights than naively applying PCA after stringing. However, their statistical
underpinnings are largely unknown. The main goal of this article is to fill in
this void. Indeed, as we shall demonstrate, a careful and rigorous statistical
treatment allows for a better understanding of the operating characteristics
of multiway PCA, leads to improved methodology, and reveals new opportunities
and challenges in analyzing multiway data.

More specifically, we focus on a simple and natural approach to multiway PCA:
when seeking linear transformations that maximize the variance, we impose the
additional constraint that they conform to the multiway structure of the data.
Doing so not only allows for enhanced interpretability but also inherits many
nice properties of the usual PCA. Just as the usual principal components (PCs)
are the eigenvectors of the covariance matrix, the multiway PCs can be
identified with certain eigenvectors of the covariance operator. To better
understand the impact of multiway structure on our ability to recover and make
inferences about the multiway PCs, we also investigate the properties of
multiway PCA under a spiked covariance model.

Statistical properties of the usual PCA are well understood in the classical
setting where the sample size is large whereas the number of variables is
small (see, e.g., Anderson, 1984). More and more often in today’s
applications, however, the dimensionality can also be large. There are
abundant theoretical results concerning the usual PCA in such a high-
dimensional setting as well, especially in the context of the spiked
covariance model. For example, Johnstone and Lu (2009) first demonstrated the
critical role of dimensionality in PCA by showing that, with fixed signal
strength, the sample PCA is consistent if and only if the number of variables
is of a smaller order than the sample size. In another influential paper, Paul
(2007) established the asymptotic distribution of sample PCs. Other related
treatments include Baik and Silverstein (2006); Nadler (2008); Johnstone and
Lu (2009); Jung and Marron (2009); Bai and Silverstein (2010); Lee et al.
(2010); Benaych-Georges and Nadakuditi (2011); Bai and Yao (2012); Shen et al.
(2013); Koltchinskii and Lounici (2014); Koltchinskii et al. (2017); Wang and
Fan (2017); Koltchinskii et al. (2020) among numerous others. In a sense, our
results naturally extend these earlier works to multiway PCA. However, the
need to work with higher-order covariance operators rather than covariance
matrices creates new and fundamental challenges and requires us to develop a
different proof strategy and several new technical tools. More importantly,
our analysis also reveals fundamental differences in behavior between the
usual PCA and multiway PCA and inspires new methodological development for the
latter.

Firstly, we establish the rates of convergence for the sample multiway PCs
under mild regularity conditions. These rates explain why it is essential that
we account for the inherent data structure when applying PCA to multiway data,
and why naively applying PCA after stringing could be problematic.
Intuitively, multiway PCA uses fewer parameters than the usual PCA and
therefore is easier for estimation. This is described precisely by our result
in that the estimation error of multiway PCs is determined by the dimension of
each mode of the data array rather than the total number of entries, and
therefore multiway PCs can be estimated accurately even if the latter far
exceeds the sample size. But a more important observation is that how well a
multiway PC can be estimated is determined by the corresponding eigenvalue of
the covariance operator, and not the gap between its eigenvalues like the
usual PCA. This somewhat surprising finding has far-reaching implications. In
particular, it means that for multiway data the PCs can be estimated well even
if their corresponding eigenvalues are not simple.

Moreover, to facilitate making statistical inferences about multiway PCs, we
derive asymptotic distributions of the sample multiway PCs. Our results again
reveal unexpected but important distinctions between multiway PCA and usual
PCA. For example, the estimated multiway PCs are asymptotically independent of
each other, and their asymptotic distribution is determined by their
corresponding eigenvalues instead of eigengaps. Furthermore, we show that bias
correction is important for the sample multiway PCs. Similar to the usual PCA,
sample multiway PC can exhibit significant bias when the dimension (of each
mode) is high. But there is also another source of bias that may arise due to
the inherent ambiguity in ordering the PCs in absence of eigengaps.
Nonetheless, we show that both types of bias can be eliminated, enabling us to
make inferences about and construct confidence intervals for the multiway PCs.

The rest of the paper is organized as follows. In Section 2, we introduce the
notion of multiway PCA both at a population level and how it works on a finite
sample. Section 3 investigates the rates of convergence for the sample
multiway PCs. Turning our attention to the asymptotic distribution of multiway
PCA in section 4, we show how to make valid inferences about the multiway PCs.
The merits of the multiway PCA and our proposed approaches are further
demonstrated through numerical experiments, both simulated and real, in
Section 5. We conclude with a summary in Section 6. Due to the space limit,
all proofs are relegated to supplementary material.

##  2 Multiway PCA

Multiway PCA can be viewed through the lens of usual PCA with the additional
multiway structure imposed on the PCs. Let
$\mathscr{X}\in\mathbb{R}^{d_{1}\times\cdots\times d_{p}}$ be an order-$p$
random array. To simplify, we shall assume in what follows that $\mathscr{X}$
is centered, i.e., $\mathbb{E}\mathscr{X}=0$, unless otherwise indicated. The
idea behind PCA is to look for a linear transformation of $\mathscr{X}$ that
maximizes the variance:

$\max_{\mathscr{W}\in\mathbb{R}^{d_{1}\times\cdots\times
d_{p}}:\|\mathscr{W}\|_{\rm F}=1}{\rm
var}(\langle\mathscr{X},\mathscr{W}\rangle).$ (1)

Here $\|\mathscr{W}\|_{\rm F}=\langle\mathscr{W},\mathscr{W}\rangle^{1/2}$ and

$\langle\mathscr{X},\mathscr{W}\rangle=\sum_{j_{1}=1}^{d_{1}}\cdots\sum_{j_{p}=1}^{d_{p}}x_{j_{1},\ldots,j_{p}}w_{j_{1},\ldots,j_{p}}.$

Denote by $\mathscr{U}_{1}$ the solution to (1). The basic premise of multiway
PCA is that $\mathscr{U}_{1}$ conforms to the multiway structure underlying
$\mathscr{X}$ in that it is a rank-one tensor and can be expressed as

$\mathscr{U}_{1}=\mathbf{u}_{1}^{(1)}\otimes\mathbf{u}_{1}^{(2)}\otimes\cdots\otimes\mathbf{u}_{1}^{(p)},$
(2)

where $\mathbf{u}_{1}^{(q)}$ is a unit length vector $\mathbb{R}^{d_{q}}$ and
$\otimes$ stands for the outer product, i.e., the $(i_{1},\ldots,i_{p})$ entry
of $\mathscr{U}_{1}$ is given by

$\left[\mathscr{U}_{1}\right]_{i_{1},\ldots,i_{p}}=u_{1i_{1}}^{(1)}u_{1i_{2}}^{(2)}\cdots
u_{1i_{p}}^{(p)}.$

In other words, $\mathscr{U}_{1}$ is also the solution to

$\max_{\mathscr{W}\in\Theta}{\rm var}(\langle\mathscr{X},\mathscr{W}\rangle),$
(3)

where $\Theta$ is the collection of all unit length rank-one tensors of
conformable dimensions, i.e.,

$\Theta=\\{\mathscr{W}=\mathbf{w}^{(1)}\otimes\mathbf{w}^{(2)}\otimes\cdots\otimes\mathbf{w}^{(p)}:\mathbf{w}^{(q)}\in\mathbb{R}^{d_{q}},\|\mathbf{w}^{(q)}\|=1,\forall
q=1,\ldots,p\\}.$

Even if the solution to (1) is not strictly rank-one as described by (2),
imposing such a constraint when seeking variance-maximizing transformation can
nonetheless be desirable because of the enhanced interpretability: the
additional rank-one constraint allows us to separate the effect along each
mode, and help address questions such as “who does what to whom and when”
which are often central to multiway data analysis. See, e.g., Kroonenberg
(2008) for further discussion and numerous motivating examples.

Subsequent PCs can be defined successively:

$\max_{\begin{subarray}{c}\mathscr{W}\in\mathbb{R}^{d_{1}\times\cdots\times
d_{p}}:\|\mathscr{W}\|_{\rm F}=1\\\
\mathscr{W}\perp\mathscr{U}_{l},l=1,\ldots,k-1\end{subarray}}{\rm
var}(\langle\mathscr{X},\mathscr{W}\rangle).$ (4)

As before, we shall consider the case when the solution has rank one. A key
requirement in defining PCs is that the $k$th PC is orthogonal to all other
PCs, i.e., $\mathscr{W}\perp\mathscr{U}_{l}$. In vector case, i.e., $p=1$,
this simply means that $\langle\mathscr{W},\mathscr{U}_{l}\rangle=0$. In
multiway case, however, there are many different notions of orthogonality.
See, e.g., Kolda (2001) for a detailed discussion on this subject. Each notion
has its own subtleties and caveats that may have different statistical
implications. In this work we shall focus on the notion of _complete
orthogonality_ : two rank-one tensors
$\mathscr{W}_{1}=\mathbf{w}_{1}^{(1)}\otimes\mathbf{w}_{1}^{(2)}\otimes\cdots\otimes\mathbf{w}_{1}^{(p)}$
and
$\mathscr{W}_{2}=\mathbf{w}_{2}^{(1)}\otimes\mathbf{w}_{2}^{(2)}\otimes\cdots\otimes\mathbf{w}_{2}^{(p)}$
are complete orthogonal if and only if
$\langle\mathbf{w}_{1}^{(q)},\mathbf{w}_{2}^{(q)}\rangle=(\mathbf{w}_{1}^{(q)})^{\top}\mathbf{w}_{2}^{(q)}=0$
for all $q=1,\ldots,p$. More specifically, the $k$th multiway PC, denoted by
$\mathscr{U}_{k}$, solves

$\max_{\mathscr{W}\in\Theta:\mathscr{W}\perp_{c}\mathscr{U}_{l},\forall
l<k}{\rm var}(\langle\mathscr{X},\mathscr{W}\rangle),$ (5)

where $\perp_{c}$ stands for complete orthogonality.

As in the case of the usual PCA, multiway PCs can also be equivalently defined
using the covariance matrix of ${\rm vec}(\mathscr{X})$. In fact, it is more
convenient to think of a covariance operator when it comes to multiway data.
More specifically, we shall view

$\Sigma:={\rm cov}(\mathscr{X})=\mathbb{E}(\mathscr{X}\otimes\mathscr{X})$

as a $d_{1}\times d_{2}\times\cdots\times d_{p}\times d_{1}\times\cdots\times
d_{p}$ array. Then for any $\mathscr{W}\in\Theta$,

${\rm
var}(\langle\mathscr{X},\mathscr{W}\rangle)=\langle\Sigma,\mathscr{W}\otimes\mathscr{W}\rangle=\langle\Sigma,\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(p)}\otimes\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(p)}\rangle.$

Write

$\lambda_{k}={\rm var}(\langle\mathscr{X},\mathscr{U}_{k}\rangle).$

Because of the symmetry of $\Sigma$,

$\mathscr{U}_{1}\otimes\mathscr{U}_{1}=\operatorname{argmax}_{\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(2p)}:\|\mathbf{w}^{(q)}\|=1}\langle\Sigma,\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(2p)}\rangle$

so that $\lambda_{1}\mathscr{U}_{1}\otimes\mathscr{U}_{1}$ is also the best
rank-one approximation to $\Sigma$ (see, e.g., Friedland, 2013). Similarly,

$\mathscr{U}_{k}\otimes\mathscr{U}_{k}=\operatorname{argmax}_{\begin{subarray}{c}\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(2p)}:\|\mathbf{w}^{(q)}\|=1\\\
\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(2p)}\perp_{c}\mathscr{U}_{l}\otimes\mathscr{U}_{l},\quad
l<k\end{subarray}}\langle\Sigma,\mathbf{w}^{(1)}\otimes\cdots\otimes\mathbf{w}^{(2p)}\rangle$

In vector case, e.g. $p=1$, $\\{(\lambda_{k},\mathscr{U}_{k}):k\geq 1\\}$ are
the eigenpairs of the covariance matrix $\Sigma$ and

$\Sigma_{r}:=\sum_{k=1}^{r}\lambda_{k}\mathscr{U}_{k}\otimes\mathscr{U}_{k}$

is the best rank-$r$ approximation to $\Sigma$, i.e.,

$\Sigma_{r}=\operatorname{argmin}_{A\in\mathbb{R}^{d_{1}\times d_{1}}:{\rm
rank}(A)\leq r}\|A-\Sigma\|.$

When $p>1$, this characterization becomes tenuous because the notion of best
low-rank approximation becomes precarious. For matrices, best low-rank
approximations can be identified with singular value decomposition thanks to
the Eckart-Young theorem. Low-rank approximation to tensors is much more
subtle and the best low-rank approximation may not exist in general. See,
e.g., Hackbusch (2012). Nonetheless, by construction, $\Sigma_{r}$ is the so-
called best rank-$r$ _greedy orthogonal approximation_ to $\Sigma$. See, e.g.,
Kolda (2001). In particular, when the multiway structure does manifest itself
in a way such that the usual PCs are rank-one tensors, for example, the
solution to (1) and (4) has rank one, then $\Sigma_{r}$ is the best low-rank
approximation to $\Sigma$.

Sample multiway PCs can also be defined in a similar fashion. Specifically,
given a sample $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ of independent copies
of $\mathscr{X}$, $\mathscr{U}_{k}$s can be estimated by maximizing the sample
variances:

$\widehat{\mathscr{U}}_{k}:=\operatorname{argmax}_{\mathscr{W}\in\Theta:\mathscr{W}\perp_{c}\mathscr{U}_{l},\forall
l<k}\frac{1}{n}\sum_{i=1}^{n}\langle\mathscr{X}_{i},\mathscr{W}\rangle^{2}.$
(6)

Let

$\widehat{\Sigma}=\frac{1}{n}\sum_{i=1}^{n}\mathscr{X}_{i}\otimes\mathscr{X}_{i}$

be the sample covariance operator. Then $\widehat{\mathscr{U}}_{1}$ can be
defined via the best rank-one approximation to $\widehat{\Sigma}$

$\widehat{\mathscr{U}}_{1}=\operatorname{argmax}_{\mathscr{W}\in\Theta}\langle\widehat{\Sigma},\mathscr{W}\otimes\mathscr{W}\rangle.$

And other PCs can also be equivalently defined as

$\widehat{\mathscr{U}}_{k}=\operatorname{argmax}_{\mathscr{W}\in\Theta:\mathscr{W}\perp_{c}\widehat{\mathscr{U}}_{l},\forall
l<k}\langle\widehat{\Sigma},\mathscr{W}\otimes\mathscr{W}\rangle.$

Note that $\widehat{\mathscr{U}}_{k}$ can also be identified with the best
rank-one approximation to a deflated covariance operator:

$\widehat{\mathscr{U}}_{k}=\operatorname{argmax}_{\mathscr{W}\in\Theta}\langle\check{\Sigma}_{k},\mathscr{W}\otimes\mathscr{W}\rangle,$

where

$\check{\Sigma}_{k}=\widehat{\Sigma}\times_{1}\widehat{{\cal
P}}_{k}^{(1)}\times_{2}\cdots\times_{p}\widehat{{\cal
P}}_{k}^{(p)}\times_{p+1}\widehat{{\cal
P}}_{k}^{(1)}\times_{p+2}\cdots\times_{2p}\widehat{{\cal P}}_{k}^{(p)}$

and $\widehat{{\cal P}}_{k}^{(1)}$ is the projection matrix of the linear
subspace spanned by $\\{\widehat{\mathbf{u}}_{l}^{(q)}:1\leq l<k\\}$.
Hereafter $\times_{q}$ represents the mode $q$ product between a tensor
$\mathscr{T}\in\mathbb{R}^{d_{1}\times d_{2}\times\dots\times d_{k}}$ and a
matrix $A\in\mathbb{R}^{m\times d_{q}}$ so that
$\mathscr{T}\times_{q}A\in\mathbb{R}^{d_{1}\times\dots d_{q-1}\times m\times
d_{q+1}\dots\times d_{k}}$ with elements

$[\mathscr{T}\times_{q}A]_{i_{1}\dots i_{q-1}ji_{q+1}\dots
i_{k}}=\sum_{i_{q}=1}^{d_{q}}\mathscr{T}_{i_{1}\dots i_{q}\dots
i_{k}}A_{ji_{q}}.$

Computing the best rank-one approximation to a tensor is a classical problem
in numerical linear algebra, and casting the sample multiway PCA as such
allows us to take advantage of the many existing algorithms for doing so. In
this work, we focus on the statistical properties of multiway PCA. Readers
interested in further discussions about the computational aspect are referred
to, e.g., Zhang and Golub (2001); Hackbusch (2012); Janzamin et al. (2019) and
references therein.

Similar to the usual PCs, multiway PCs can be used to construct low-rank
approximations of the original data. However, there are also fundamental,
albeit sometimes subtle, differences between the two types of PCA. The usual
sample PCs coincide with the leading singular vectors of the data matrix after
appropriate centering and therefore can be computed via singular value
decomposition (SVD). In contrast, multiway PCA is, while closely related to,
not equivalent to the best low-rank approximations of the original data array
in general. More specifically, consider stacking the observations into a
higher-order tensor $\mathbf{X}\in\mathbb{R}^{n\times d_{1}\times\cdots\times
d_{p}}$ whose $i$th frontal slice is $\mathscr{X}_{i}$. In the case when
$\mathscr{X}$ is a vector, i.e., $p=1$, $\mathbf{X}$ is a matrix and the
sample PC, $\widehat{\mathscr{U}}_{k}$ as defined above, is its $k$th right
singular vector. It is therefore tempting to do the same and estimate
$\mathscr{U}_{k}$s by seeking the best orthogonal low-rank approximation to
$\mathbf{X}$ directly:

$\min_{\begin{subarray}{c}\mathscr{W}_{1},\ldots,\mathscr{W}_{r}\in\Theta,\mathbf{a}_{1},\ldots,\mathbf{a}_{r}\in\mathbb{R}^{n}\\\
\mathscr{W}_{l}\perp_{c}\mathscr{W}_{k},\mathbf{a}_{l}\perp\mathbf{a}_{k},\forall
l\neq
k\end{subarray}}\left\|\mathbf{X}-\left(\mathbf{a}_{1}\otimes\mathscr{W}_{1}+\cdots+\mathbf{a}_{r}\otimes\mathscr{W}_{r}\right)\right\|_{\rm
F}$ (7)

See, e.g., Harshman and Lundy (1984). This problem, often known as the tensor
SVD problem, has attracted a lot of attention in recent years. See, e.g.,
Richard and Montanari (2014); Hopkins et al. (2015); Liu et al. (2017); Zhang
and Xia (2018); Auddy and Yuan (2020). However, the sample multiway PCs are
generally not the solution to (7). First of all, the difference between the
best orthogonal rank-$r$ and rank-$(r-1)$ approximations to $\mathbf{X}$ is
generally not a rank-one tensor and therefore cannot be associated with a
multiway PC. See, e.g., Hackbusch (2012). To overcome this challenge, one may
consider solving (7) in a greedy fashion, i.e, optimizing (7) over
$\mathscr{W}_{k}$ and $\mathbf{a}_{k}$ only while fixing the other ones. In
general, however, this still results in a different set of PCs because of the
extra orthogonality constraint on $\mathbf{a}_{k}$s imposed by (7). As we
shall see, this subtle distinction between multiway PCA and low-rank
approximations to a data tensor not only means that a treatment different from
that for the tensor SVD is needed for multiway PCA but also leads to different
statistical behavior between the two.

##  3 Rates of Convergence

A natural question one first asks is how well $\mathscr{U}_{k}$ and its
components $\mathbf{u}_{k}^{(q)}$s can be estimated by their sample
counterparts. We shall now turn our attention to this question and study the
rate of convergence for the sample multiway PCs. On the one hand, we provide
further justification for the superiority of multiway PCA to the usual PCA
with stringing, in addition to enhanced interpretability. On the other hand,
our investigation also leads to new insights into the operating
characteristics of sample multiway PCA and its intriguing distinction from the
usual PCA. To fix ideas, we shall consider the so-called spiked covariance
model as a working model for our theoretical development.

Suppose that a random array $\mathscr{X}\in\mathbb{R}^{d_{1}\times\cdots\times
d_{p}}$ follows a linear factor model:

$\mathscr{X}=\sum_{k=1}^{r}\sigma_{k}\theta_{k}\mathscr{U}_{k}+\sigma_{0}\mathscr{E},$
(8)

where $(\theta_{1},\ldots,\theta_{r})^{\top}\sim N(0,I_{r})$ are the random
factor loadings, $\mathscr{U}_{k}$s ($\in\Theta$) are unit length rank-one
principal components such that $\mathscr{U}_{k}\perp_{c}\mathscr{U}_{l}$ for
any $k\neq l$, and $\mathscr{E}$ is a noise tensor with independent $N(0,1)$
entries. It is worth pointing out that our results and arguments can be
extended beyond normality and applied to general subgaussian distributions. We
opt for the normality assumption for ease of presentation. Without loss of
generality, we shall also assume that eigenvalues of the signal are nontrivial
and sorted in non-increasing order, i.e.,
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{r}>0$. Note that we do not
require $\sigma_{k}$s to be distinct. It is not hard to see that the
covariance operator of the aforementioned $\mathscr{X}$ is given by

$\Sigma=\sum_{k=1}^{r}\sigma_{k}^{2}\mathscr{U}_{k}\otimes\mathscr{U}_{k}+\sigma_{0}^{2}\mathscr{I},$

where $\mathscr{I}$ is the identity tensor, i.e., $\mathscr{I}_{j_{1}\ldots
j_{p}j_{1}^{\prime}\ldots j_{p}^{\prime}}=1$ if $j_{q}=j_{q}^{\prime}$ for all
$q=1,\ldots,p$ and $0$ otherwise. The spiked covariance model such as (8) is
widely used as a working model to study PCA in the case of vector
observations, i.e., $p=1$. See, e.g., Johnstone (2001) and Paul (2007).

In this section, we shall establish the rates of convergence of the sample
multiway PCs. To this end, denote by $\angle(\mathbf{w}_{1},\mathbf{w}_{2})$
the angle between two vectors $\mathbf{w}_{1}$ and $\mathbf{w}_{2}$ taking
value in $[0,\pi/2]$, and similarly for two arrays $\mathscr{W}_{1}$ and
$\mathscr{W}_{2}$, $\angle(\mathscr{W}_{1},\mathscr{W}_{2})$ denotes the angle
between their vectorizations ${\rm vec}(\mathscr{W}_{1})$ and ${\rm
vec}(\mathscr{W}_{2})$.

It is instructive to begin with the classical setting where the dimensionality
$d_{1},\ldots,d_{p}$ as well as all other parameters, e.g. $\sigma_{0}$,
$\sigma_{k}$s and $r$, are held fixed as the sample size $n$ diverges. Our
first result shows that the sample PC $\widehat{\mathscr{U}}_{k}$ and its
components $\widehat{\mathbf{u}}_{k}^{(q)}$s are root-$n$ consistent in this
regime.

######  Theorem 3.1.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Assume that all parameters are fixed as the sample size
$n$ increases. Let
$\widehat{\mathscr{U}}_{k}=\widehat{\mathbf{u}}_{k}^{(1)}\otimes\cdots\otimes\widehat{\mathbf{u}}_{k}^{(p)}$
be the sample multiway PC as defined by (6). Then there exists a permutation
$\pi$ over $[r]:=\\{1,\ldots,r\\}$ such that

$\max_{1\leq q\leq
p}\sin\angle(\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{u}_{\pi(k)}^{(q)})=O_{p}(n^{-1/2}),$
(9)

for all $k\in[r]$, and hence

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})=O_{p}(n^{-1/2}),\qquad
k=1,\ldots,r$

as $n\to\infty$.

The most notable difference between the above result and those for the usual
PCA (e.g., Anderson, 1984) is that fact that the root-$n$ consistency of the
sample multiway PCs does not require that the eigenvalues
($(\sigma_{k}^{2}+\sigma_{0}^{2})$s or equivalently $\sigma_{k}$s) of the
covariance matrix be simple, i.e., $\sigma_{k}\neq\sigma_{k+1}$. Note that,
without the multiway structural constraint, the usual PCs are only uniquely
defined and hence can possibly be estimated if their corresponding eigenvalues
are simple. As Theorem 3.1 indicates, such a restriction is not necessary for
multiway PCA. For multiway PCA, each sample PC is root-$n$ consistent
regardless of the other eigenvalues. It is also worth noting that, since we do
not require the $\sigma_{k}$s to be distinct, there is no guarantee that
$\widehat{\mathscr{U}}_{k}$ estimates $\mathscr{U}_{k}$. This is not a
deficiency of multiway PCA, but rather a necessity due to the possible
indeterminacy of the $k$th largest eigenvalue. In fact, if
$\sigma_{k+1}<\sigma_{k}<\sigma_{k-1}$, then we can choose $\pi(k)=k$ in
Theorem 3.1. In general, Theorem 3.1 shows that each of the sample PCs is
necessarily a root-$n$ consistent estimate of one of the multiway PCs.

To further understand the operating characteristics and merits of multiway
PCA, we now consider the more general case and further highlight the role of
dimensionality and signal-to-noise ratio. For brevity, in what follows, we
shall assume that $\mathscr{X}$ is “nearly cubic” in that there exist
constants $0<c_{1},c_{2}<\infty$ such that $c_{1}d\leq d_{1},\ldots,d_{p}\leq
c_{2}d$ for some natural number $d$ which may diverge with $n$. General cases
can be treated similarly but incur considerably more cumbersome notation and
tedious derivation.

######  Theorem 3.2.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$.
Let
$\widehat{\mathscr{U}}_{k}=\widehat{\mathbf{u}}_{k}^{(1)}\otimes\cdots\otimes\widehat{\mathbf{u}}_{k}^{(p)}$
be the sample multiway PC as defined by (6). Suppose that

$r\log r\leq c_{0}\min\\{n,d\\}\quad{\rm
and}\quad\left(\frac{\sigma_{0}}{\sigma_{r}}+\frac{\sigma_{0}^{2}}{\sigma_{r}^{2}}\right)\cdot\max\left\\{\sqrt{\frac{d}{n}},{\frac{d}{n}}\right\\}\leq\frac{c_{0}}{\sqrt{r}},$
(10)

for a sufficiently small constant $c_{0}>0$. Then there exist a constant $C>0$
and a permutation $\pi$ over $[r]$ such that

$\max_{1\leq q\leq
p}\sin\angle(\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{u}_{\pi(k)}^{(q)})\leq
C\left(\frac{\sigma_{0}}{\sigma_{\pi(k)}}+\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}\right)\cdot\max\left\\{\sqrt{\frac{d}{n}},{\frac{d}{n}}\right\\},$
(11)

for all $k\in[r]$, and hence

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})\leq
C\left(\frac{\sigma_{0}}{\sigma_{\pi(k)}}+\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}\right)\cdot\max\left\\{\sqrt{\frac{d}{n}},{\frac{d}{n}}\right\\},\qquad
k=1,\ldots,r,$

with probability tending to one as $n$ diverges.

Theorem 3.2 can be viewed as a generalization of Theorem 3.1. Its proof is
rather involved and we shall brief discuss some of the challenges and the main
ideas for resolving them. The proof proceeds by induction over $k$. Special
attention is needed to deal with the case when an eigenvalue is not simple or
the eigengap is small. This creates difficulty in identifying which multiway
PC a sample multiway PC estimates, or equivalently the permutation $\pi$. To
this end, we shall define

$\pi(1)=\operatorname{argmax}_{1\leq l\leq
r}\left\\{\sigma_{l}^{2}\left|\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{1}^{(q)}\rangle\right|\right\\},$

and for $k>1$,

$\pi(k):=\operatorname{argmax}_{l\notin\pi([k-1])}\left\\{\sigma_{l}^{2}\left|\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right|\right\\}.$

To remove the influence of eigengaps altogether, we need to carefully quantify
the impact of estimation error of
$\widehat{\mathscr{U}}_{1},\ldots,\widehat{\mathscr{U}}_{k-1}$ on the $k$th
sample multiway PC. To this end, we shall derive bounds for both

$\max_{1\leq q\leq
p}\sin\angle(\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}),$

and

$\max_{1\leq q\leq
p}\max_{l\notin\pi([k])}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle,$

and leverage the fact that the latter can be much smaller than the former.

When $d=O(n)$, the convergence rate given in Theorem 3.2 is

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})\leq
C\left(\frac{\sigma_{0}}{\sigma_{\pi(k)}}+\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}\right)\cdot\sqrt{\frac{d}{n}};$

and when $d\gg n$, we have

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})\leq
C\left(\frac{\sigma_{0}}{\sigma_{\pi(k)}}+\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}\right)\cdot\frac{d}{n}.$

In particular, $\widehat{\mathscr{U}}_{k}$ is consistent, e.g.,

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})=o_{p}(1),$

whenever $\sigma_{\pi(k)}/\sigma_{0}\gg\max\\{d/n,(d/n)^{1/4}\\}$.

Of particular interest here is the role of dimensionality. The rates of
convergence given by Theorem 3.2 depend on the dimensionality through $d$
rather than the ambient dimension $D:=d_{1}d_{2}\cdots d_{p}$. This is because
multiway PCA restricts PCs to be rank-one tensors and therefore has fewer
parameters. Such dimensionality reduction is especially important for multiway
data. Consider, for example, the case when $\sigma_{0}$ and $\sigma_{k}$s are
fixed, then by virtue of the results from Johnstone and Lu (2009), direct
application of the usual PCA after stringing necessarily leads to an
inconsistent estimate of $\mathscr{U}_{k}$ whenever $D\gg n$. Yet, our result
indicates that multiway PCA is consistent whenever $d\ll n$.

To draw further comparisons with the usual PCA, we now focus on the case when
$d\ll n$ and $r$, $\sigma_{0},\sigma_{1},\ldots,\sigma_{r}$ are fixed. As
shown by Birnbaum et al. (2013), in this regime, the usual PCA (i.e., $p=1$)
satisfies

$\sin\angle(\widehat{\mathscr{U}}_{k},\mathscr{U}_{\pi(k)})\asymp\left(\frac{\sigma_{0}}{\sigma_{\pi(k)}}+\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}\right)\cdot\sqrt{\frac{d}{n}}+{1\over\sqrt{n}}\left(\sum_{k^{\prime}\neq
k}{(\sigma_{0}+\sigma_{\pi(k)})(\sigma_{0}+\sigma_{\pi(k^{\prime})})\over\sigma_{\pi(k)}^{2}-\sigma_{\pi(k^{\prime})}^{2}}\right)$

with probability tending to one. Comparing the above rate with that from
Theorem 3.2, it is clear that the difference between the two lie at the second
term on the right hand side. Its presence for the usual PCA dictates that
there should be no ties among $\sigma_{k}$s. Even if the $\sigma_{k}$s are all
distinct, how well we can estimate a PC crucially depends on the gap between
its corresponding eigenvalue and the other eigenvalues when $p=1$. In
contrast, the bounds given by Theorem 3.2 are determined by $\sigma_{\pi(k)}$
alone and not the eigengap
$\min\\{\sigma_{\pi(k)-1}^{2}-\sigma_{\pi(k)}^{2},\sigma_{\pi(k)}^{2}-\sigma_{\pi(k)+1}^{2}\\}$
as in the usual PCA case.

It is also instructive to compare the convergence rate for the multiway PCA
from Theorem 3.2 with those for tensor SVD. Recall that

$\mathbf{X}=\sum_{k=1}^{r}\sigma_{k}\Theta_{k}\otimes\mathscr{U}_{k}+\sigma_{0}\mathbf{E},$

where $\Theta_{k}=(\theta_{1k},\ldots,\theta_{nk})^{\top}$ is a vector
containing the $n$ realizations of $\theta_{k}$ and $\mathbf{E}$ is a $n\times
d_{1}\times\cdots\times d_{p}$ tensor whose $i$th frontal slice is
$\mathscr{E}_{i}$. In contrast, $\Theta_{k}$s are deterministic in a tensor
SVD model. If $\Theta_{k}$s are orthogonal to each other, then
$\mathscr{U}_{k}$ can be estimated at the rate of
$\sigma_{0}\|\mathbf{E}\|/(\sigma_{k}\|\Theta_{k}\|)$ which is of the order
$(\sigma_{0}/\sigma_{k})\max\\{\sqrt{d/n},1\\}$. This is a direct consequence
of the perturbations bounds from Auddy and Yuan (2020) and a similar bound was
also derived by Richard and Montanari (2014) in the rank-one case, i.e.,
$r=1$. In our case, however, $\Theta_{k}$ and $\Theta_{l}$ are random and in
general not orthogonal to each other. As a result, the rates we obtained are
different in their dependence on the signal-to-noise ratio
$\sigma_{k}/\sigma_{0}$. Similar phenomenon has also been observed for the
usual PCA (see, e.g., Birnbaum et al., 2013).

##  4 Asymptotic Normality and Bias Correction

We now turn to the distributional properties of multiway PCA. This requires us
to further delineate the role of bias in the sample PCs. It is known that the
usual PCA is biased when the dimension ($D$) is large when compared with the
sample size. See, e.g., Koltchinskii and Lounici (2014); Koltchinskii et al.
(2020) and the references therein. The same phenomenon is observed for the
sample multiway PCs and a non-negligible bias arise when the dimension of each
mode ($d$) is large when compared with the sample size. In addition, there is
a more subtle source of bias for the sample multiway PCs due to the ambiguity
in ordering the multiway PCs in the absence of eigengaps. As noted before, the
lack of an eigengap means that the $k$th PC may not necessarily be estimated
by the $k$th sample multiway PC. As a more concrete example, consider the case
when $r=2$ and $\lambda_{1}=\lambda_{2}$. Then $\mathscr{U}_{1}$ can be
estimated by either $\widehat{\mathscr{U}}_{1}$ or
$\widehat{\mathscr{U}}_{2}$, and as Theorem 3.2 shows, the rate of convergence
remains the same in both cases. But the asymptotic distribution may differ
between the two scenarios: $\widehat{\mathscr{U}}_{2}$ is required to be
orthogonal to $\widehat{\mathscr{U}}_{1}$ and estimating $\mathscr{U}_{1}$ by
$\widehat{\mathscr{U}}_{2}$ may incur extra bias.

In this section, we shall introduce ways to correct for both types of bias and
establish the asymptotic normality of the bias-corrected sample PCs. As is
customary in the literature, we shall assume that $r$ and
$\sigma_{1},\ldots,\sigma_{r}$ are fixed for brevity. In light of the results
from the previous section, the sample PCs are consistent if $d\ll n$ in this
setting. We shall therefore focus on this regime in the current section.

###  4.1 When $d=o(\sqrt{n})$

When $d$ is not too large, the bias is solely due to the possibility of
repeated eigenvalues and thus ambiguity of the ordering of PCs. Indeed if
$\sigma_{k}$s are distinct, then there is no need for bias correction when
$d=o(\sqrt{n})$ and all of our results in this subsection will hold for the
sample multiway PCs. But in practice, we may not know or want to assume that
the eigenvalues are simple. Fortunately, we can remove any possible bias
fairly easily by a simple one-step update of the sample PCs. More
specifically, we shall consider estimating $\mathbf{u}_{\pi(k)}^{(q)}$ by
$\tilde{\mathbf{u}}_{k}^{(q)}$, the leading eigenvector of

$\widehat{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)}).$

The additional step frees up the orthogonality constraints imposed on the
$k$th sample multiway PC and therefore allows us to suppress any adverse
influence of $\widehat{\mathscr{U}}_{1},\ldots,\widehat{\mathscr{U}}_{k-1}$.

We now consider the asymptotic distribution of the bias-corrected sample PCs.
We again start with the classical regime when all parameters are fixed as $n$
increases.

######  Theorem 4.1.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Assume that all parameters are fixed as the sample size
$n$ increases. Let
$\tilde{\mathbf{u}}_{1}^{(q)},\ldots,\tilde{\mathbf{u}}_{r}^{(q)}$ be defined
as above. Then there exists a permutation $\pi:[r]\to[r]$ such that

$\displaystyle\sqrt{n}\left[{\rm vec}(\tilde{\mathbf{U}}^{(q)})-{\rm
vec}(\mathbf{U}^{(q)}_{\pi})\right]$
$\displaystyle\overset{d}{\to}N\left(0,{\rm
diag}\left(\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(1)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(1)}^{4}}\right){\cal
P}_{\mathbf{u}_{\pi(1)}^{(q)}}^{\perp},\ldots,\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(r)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(r)}^{4}}\right){\cal
P}_{\mathbf{u}_{\pi(r)}^{(q)}}^{\perp}\right)\right),$

as $n\to\infty$, where
$\tilde{\mathbf{U}}^{(q)}=[\tilde{\mathbf{u}}_{1}^{(q)},\ldots,\tilde{\mathbf{u}}_{r}^{(q)}]$,
$\mathbf{U}_{\pi}^{(q)}=[\mathbf{u}_{\pi(1)}^{(q)},\ldots,\mathbf{u}_{\pi(r)}^{(q)}]$
and ${\cal
P}_{\mathbf{u}_{k}^{(q)}}^{\perp}=I_{d_{q}}-\mathbf{u}_{k}^{(q)}\otimes\mathbf{u}_{k}^{(q)}$.

Theorem 4.1 indicates that

$n\cdot{\rm
var}\left(\tilde{\mathbf{u}}_{k}^{(q)}\right)\to{\sigma_{0}^{2}(\sigma_{\pi(k)}^{2}+\sigma_{0}^{2})\over\sigma_{\pi(k)}^{4}}\left(I-\mathbf{u}_{\pi(k)}^{(q)}\otimes\mathbf{u}_{\pi(k)}^{(q)}\right),$

and

$n\cdot{\rm
cov}\left(\tilde{\mathbf{u}}_{k}^{(q)},\tilde{\mathbf{u}}_{l}^{(q)}\right)\to
0.$

Namely, all estimates of the multiway PCs are asymptotically normal and
independent of each other. Note also that the asymptotic distribution of
$\tilde{\mathbf{u}}_{k}^{(q)}$ does not depend on other eigenvalues or PCs. In
other words, it can be estimated to the same precision as if all other
components $\mathscr{U}_{l}$, $l\neq k$ are known! This is to be contrasted
with the usual PCA where the asymptotic distribution of $\mathbf{u}_{k}^{(q)}$
depends on all other eigenvectors and eigenvalues.

More specifically, it is well known that in vector case, i.e., when $p=1$,
under the additional assumption that $\sigma_{1}^{2},\ldots,\sigma_{r}^{2}$
are distinct, the sample PCs satisfy

$\displaystyle n\cdot{\rm var}\left(\widehat{\mathbf{u}}_{k}^{(1)}\right)$
$\displaystyle\to\sum_{1\leq l\leq r,l\neq
k}{(\sigma_{k}^{2}+\sigma_{0}^{2})(\sigma_{l}^{2}+\sigma_{0}^{2})\over(\sigma_{k}^{2}-\sigma_{l}^{2})^{2}}\mathbf{u}_{l}^{(1)}\otimes\mathbf{u}_{l}^{(1)}+{(\sigma_{k}^{2}+\sigma_{0}^{2})\sigma_{0}^{2}\over\sigma_{k}^{4}}\left(I-\sum_{1\leq
l\leq r}\mathbf{u}_{l}^{(1)}\otimes\mathbf{u}_{l}^{(1)}\right)$

and for any $1\leq l\leq r$ and $l\neq k$,

$n\cdot{\rm
cov}\left(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{l}^{(1)}\right)\to-{(\sigma_{k}^{2}+\sigma_{0}^{2})(\sigma_{l}^{2}+\sigma_{0}^{2})\over(\sigma_{k}^{2}-\sigma_{l}^{2})^{2}}\cdot\mathbf{u}_{k}^{(1)}\otimes\mathbf{u}_{l}^{(1)}.$

See, e.g., Anderson (1984). It is clear that the sample PCs are always
correlated with each other. Moreover, note that

${\sigma_{0}^{2}(\sigma_{k}^{2}+\sigma_{0}^{2})\over\sigma_{k}^{4}}\leq{(\sigma_{k}^{2}+\sigma_{0}^{2})(\sigma_{l}^{2}+\sigma_{0}^{2})\over(\sigma_{k}^{2}-\sigma_{l}^{2})^{2}}.$

and the strict inequality holds for any $k\neq l\leq r$. This suggests that
the estimated multiway PCs have smaller variations than the usual PCs with the
same set of eigenvalues.

We now turn our attention to the more general case when the dimensionality and
other parameters are allowed to diverge with $n$. Because the PCs now may have
different dimensions for different sample sizes, it is more natural to
consider their linear forms, e.g.
$\langle\mathbf{u}^{(q)}_{k},\mathbf{v}\rangle$, for some fixed vector
$\mathbf{v}\in\mathbb{R}^{d_{q}}$. If the dimensions are fixed, Theorem 4.1
immediately suggests that
$\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle$ estimates
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$, and

$\sqrt{n}\left(\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\to_{d}N\left(0,\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)\|{\cal
P}_{\mathbf{u}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\|^{2}\right)$

The following result shows that this continues to hold as long as
$d=o(\sqrt{n})$.

######  Theorem 4.2.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Assume that (10) holds, $d=o(\sqrt{n})$. Then there exists
a permutation $\pi:[r]\to[r]$ such that

$\sqrt{n}\left(\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\to_{d}N\left(0,\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)\|{\cal
P}_{\mathbf{u}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\|^{2}\right)$

as $n\to\infty$.

Theorem 4.2 shows that the same asymptotic behavior of
$\tilde{\mathbf{u}}_{j}^{(q)}$ as in the fixed dimension case can be expected
whenever $d=o(\sqrt{n})$.

###  4.2 When $d=o(n)$

For higher dimension, the simple bias-correction described above is no longer
sufficient and a close inspection reveals that $\tilde{\mathbf{u}}_{j}^{(q)}$
still incurs a non-negligible bias when $d\gg\sqrt{n}$. Thankfully, both types
of bias can be corrected with a sample-splitting approach similar in spirit to
the scheme developed by Koltchinskii and Lounici (2014) for the usual PCA.

Without loss of generality, assume that $n$ is an even number and we randomly
split the $n$ observations into two halves:
$(\mathscr{X}_{1},\ldots,\mathscr{X}_{n/2})$ and
$(\mathscr{X}_{n/2+1},\ldots,\mathscr{X}_{n})$. Denote by
$\widehat{\Sigma}^{[1]}$ and $\widehat{\Sigma}^{[2]}$ the sample covariance
operator based on the two halves of data respectively. Similarly, we shall
write
$\widehat{\mathscr{U}}_{k}^{[h]}=\widehat{\mathbf{u}}_{k}^{(1),[h]}\otimes\cdots\otimes\widehat{\mathbf{u}}_{k}^{(p),[h]}$
the $k$th sample PC based on the $h$ ($=1$ or $2$) halves of the data.
However, as noted before, $\widehat{\mathscr{U}}_{k}^{[1]}$ and
$\widehat{\mathscr{U}}_{k}^{[2]}$ may not estimate the same PC. To this end,
we shall reorder $\widehat{\mathscr{U}}_{k}^{[1]}$s and
$\widehat{\mathscr{U}}_{k}^{[2]}$s with $\widehat{\mathscr{U}}_{k}$s (i.e.,
the estimators derived from the entire dataset) as reference points.
Specifically, without loss of generality, we assume that

$\widehat{\mathscr{U}}_{k}^{[1]}=\operatorname{argmin}_{k\leq l\leq
r}\sin\angle\left(\widehat{\mathscr{U}}_{l}^{[1]},\widehat{\mathscr{U}}_{k}\right).$

The same procedure is applied to relabel $\widehat{\mathscr{U}}_{k}^{[2]}$s.
Note also that the sign of a PC is irrelevant in that $\mathscr{U}_{k}$ and
$-\mathscr{U}_{k}$ represent the same transformation. We shall therefore also
assume hereafter, without loss of generality, that
$\langle\widehat{\mathbf{u}}_{k}^{(q),[1]},\widehat{\mathbf{u}}_{k}^{(q),[2]}\rangle\geq
0$.

Recall that

$\displaystyle\sigma_{k}^{2}\mathbf{u}_{k}^{(q)}{\mathbf{u}_{k}^{(q)}}^{\top}+\sigma_{0}^{2}I_{d_{q}}=\Sigma(\mathbf{u}_{k}^{(1)},\ldots,\mathbf{u}_{k}^{(q-1)},\cdot,\mathbf{u}_{k}^{(q+1)},\ldots,\mathbf{u}_{k}^{(p)},$
$\displaystyle\mathbf{u}_{k}^{(1)},\ldots,\mathbf{u}_{k}^{(q-1)},\cdot,\mathbf{u}_{k}^{(q+1)},\ldots,\mathbf{u}_{k}^{(p)}).$

We shall then update the sample PC using the above identity with $\Sigma$ and
$\mathbf{u}_{k}^{(q)}$s estimated from separate halves. Denote by
$\check{\mathbf{u}}_{k}^{(q),[1]}$ the leading eigenvector of

$\displaystyle\widehat{\Sigma}^{[1]}$
$\displaystyle(\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]},$
$\displaystyle\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]}),$

and similarly $\check{\mathbf{u}}_{k}^{(q),[2]}$ the leading eigenvector of

$\displaystyle\widehat{\Sigma}^{[2]}$
$\displaystyle(\widehat{\mathbf{u}}_{k}^{(1),[1]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[1]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[1]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[1]},$
$\displaystyle\widehat{\mathbf{u}}_{k}^{(1),[1]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[1]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[1]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[1]}).$

To avoid losing efficiency due to sample splitting, we consider a new estimate
$\check{\mathscr{U}}_{k}=\check{\mathbf{u}}_{k}^{(1)}\otimes\cdots\otimes\check{\mathbf{u}}_{k}^{(p)}$
where

$\check{\mathbf{u}}_{k}^{(q)}={\check{\mathbf{u}}_{k}^{(q),[1]}+\check{\mathbf{u}}_{k}^{(q),[2]}\over\left\|\check{\mathbf{u}}_{k}^{(q),[1]}+\check{\mathbf{u}}_{k}^{(q),[2]}\right\|}.$

The following theorem shows that we can construct an unbiased estimate of
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$ by appropriately
rescaling $\langle\check{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle$, as long as
$d=o(n^{2/3})$.

######  Theorem 4.3.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Let
$\check{\mathscr{U}}_{k}=\check{\mathbf{u}}_{k}^{(1)}\otimes\cdots\otimes\check{\mathbf{u}}_{k}^{(p)}$
be the estimated PC as defined above. Assume $r$ and
$\sigma_{1},\dots,\sigma_{r}$ are fixed, and $d=o(n^{2/3})$. Then there exists
a permutation $\pi:[r]\to[r]$ such that

$\sqrt{n}\left((1+b_{k}^{(q)})\langle\check{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\to_{d}N\left(0,\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)\|{\cal
P}_{\mathbf{u}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\|^{2}\right)$

as $n\to\infty$ where

$b_{k}^{(q)}=\sqrt{1+\frac{d_{q}}{n}\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)}-1.$
(12)

It is worth pointing out that when $d=o(n^{1/2})$, the bias correction factor
described by (12) obeys $b_{k}^{(q)}=o(n^{-1/2})$ and therefore can be
neglected. This agrees with our earlier observation and of course also
suggests that sample-splitting is unnecessary if $d\ll n^{1/2}$. When $d\gg
n^{1/2}$, bias correction becomes essential. In particular, Theorem 4.3
suggests that, as long as $d=o(n^{2/3})$, an explicit bias correction factor
can be applied. For higher dimensions, it is unclear if a similar explicit
expression exists for the debiasing factor. Nonetheless, we can derive a
suitable bias correction factor for all $d\ll n$ via additional sample
splitting.

More specifically, we first randomly split the observations into two halves.
The first half of the data is then further split into two equal-sized groups
to compute the sample covariance operators $\widehat{\Sigma}^{[1][1]}$ and
$\widehat{\Sigma}^{[1][2]}$, then we compute
$\widehat{\mathbf{u}}_{k}^{(q),[1][1]}$ and
$\widehat{\mathbf{u}}_{k}^{(q),[1][2]}$ as the leading eigenvectors of

$\displaystyle\widehat{\Sigma}^{[1][1]}$
$\displaystyle(\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]},$
$\displaystyle\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]}),$
$\displaystyle\widehat{\Sigma}^{[1][2]}$
$\displaystyle(\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]},$
$\displaystyle\widehat{\mathbf{u}}_{k}^{(1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(q-1),[2]},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1),[2]},\dots,\widehat{\mathbf{u}}_{k}^{(p),[2]}),$

Similarly, we used the second half of the data to compute
$\widehat{\mathbf{u}}_{k}^{(q),[2][1]}$s, and
$\widehat{\mathbf{u}}_{k}^{(q),[2][2]}$s. As before, we shall sort these
estimates in compatible order and sign. Let

$\widehat{b}_{k}^{(q)}=\frac{\left\|\check{\mathbf{u}}_{k}^{(q),[1]}+\check{\mathbf{u}}_{k}^{(q),[2]}\right\|}{\sqrt{\left\langle\widehat{\mathbf{u}}_{k}^{(q),[1][1]},\widehat{\mathbf{u}}_{k}^{(q),[1][2]}\right\rangle}+\sqrt{\left\langle\widehat{\mathbf{u}}_{k}^{(q),[2][1]},\widehat{\mathbf{u}}_{k}^{(q),[2][2]}\right\rangle}}-1.$
(13)

######  Theorem 4.4.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Let
$\check{\mathscr{U}}_{k}=\check{\mathbf{u}}_{k}^{(1)}\otimes\cdots\otimes\check{\mathbf{u}}_{k}^{(p)}$
be the estimated PC as defined above. Assume $r$ and
$\sigma_{1},\dots,\sigma_{r}$ are fixed, and $d=o(n)$. Then there exists a
permutation $\pi:[r]\to[r]$ such that

$\sqrt{n}\left((1+\widehat{b}_{k}^{(q)})\langle\check{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\to_{d}N\left(0,\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)\|{\cal
P}_{\mathbf{u}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\|^{2}\right),$

as $n\to\infty$ where $\widehat{b}_{k}^{(q)}$ is given by (13). Moreover,

$\displaystyle\widehat{b}_{k}^{(q)}=\sqrt{1+\frac{d_{q}}{n}\left(\frac{\sigma_{0}^{2}}{\sigma_{\pi(k)}^{2}}+\frac{\sigma_{0}^{4}}{\sigma_{\pi(k)}^{4}}\right)}-1+O_{p}\left(\frac{d^{3/2}}{n^{3/2}}\right)+o_{p}\left(\frac{1}{\sqrt{n}}\right).$

In light of Theorem 4.4, the double sample splitting approach can be employed
to derive confidence intervals for linear forms of the multiway PCs as long as
$d=o(n)$. This robustness, however, comes at the expense of increased
computational cost and could incur a loss of efficiency in finite samples. In
practice, one may still prefer the explicit bias correction as described by
Theorem 4.3 if $d$ is not very large, or the one-step update if $d$ is small.

###  4.3 Inference about multiway PCs

The asymptotic normality we showed earlier in the section forms the basis for
making inferences about linear forms
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$. In particular, one of
the most interesting and also simplest examples of linear forms of PCs is
their coordinates, i.e., $\mathbf{v}$ is a column vector of the identity
matrix. To derive confidence intervals of or testing hypotheses about
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$, however, we need to also
estimate its variance. Specifically, its asymptotic distribution depends only
on $\sigma_{0}$, $\sigma_{\pi(k)}$, and $\mathbf{u}_{k}^{(q)}$, all of which
can be consistently estimated by their sample counterpart. Let

$\widehat{\sigma}_{0}^{2}={1\over\prod_{q=1}^{p}(d_{q}-r)}\sum_{1\leq
i_{q}\leq d_{q},1\leq q\leq p}[\check{\Sigma}_{r+1}]_{i_{1}\cdots
i_{p}i_{1}\cdots i_{p}}$

and

$\widehat{\sigma}_{\pi(k)}^{2}=\langle\widehat{\Sigma},\widehat{\mathscr{U}}_{k}\otimes\widehat{\mathscr{U}}_{k}\rangle-\widehat{\sigma}_{0}^{2}.$

The following theorem suggest that the asymptotic normality remains valid if
we replace the variance of linear forms
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$ with these estimates:

######  Theorem 4.5.

Let $\mathscr{X}_{1},\ldots,\mathscr{X}_{n}$ be independent observations
following the spiked covariance model (8) with $p>1$ such that
$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\cdots\otimes\mathbf{u}_{k}^{(p)}$
and $\sigma_{k}>0$. Assume $r$ and $\sigma_{1},\dots,\sigma_{r}$ are fixed.
There exists a permutation $\pi:[r]\to[r]$ such that

  * (a)

If $d=o(\sqrt{n})$, then

${\sqrt{n}\left(\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\over\sqrt{\frac{\widehat{\sigma}_{0}^{2}}{\widehat{\sigma}_{\pi(k)}^{2}}+\frac{\widehat{\sigma}_{0}^{4}}{\widehat{\sigma}_{\pi(k)}^{4}}}\
\big{\|}{\cal
P}_{\tilde{\mathbf{u}}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\big{\|}}\to_{d}N\left(0,1\right),$

  * (b)

If $d=o(n^{2/3})$, then

${\sqrt{n}\left((1+b_{k}^{(q)})\langle\check{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\over\sqrt{\frac{\widehat{\sigma}_{0}^{2}}{\widehat{\sigma}_{\pi(k)}^{2}}+\frac{\widehat{\sigma}_{0}^{4}}{\widehat{\sigma}_{\pi(k)}^{4}}}\
\big{\|}{\cal
P}_{\check{\mathbf{u}}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\big{\|}}\to_{d}N\left(0,1\right),$

where $b_{k}^{(q)}$ is given by (12).

  * (c)

If $d=o(n)$, then

${\sqrt{n}\left((1+\widehat{b}_{k}^{(q)})\langle\check{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\over\sqrt{\frac{\widehat{\sigma}_{0}^{2}}{\widehat{\sigma}_{\pi(k)}^{2}}+\frac{\widehat{\sigma}_{0}^{4}}{\widehat{\sigma}_{\pi(k)}^{4}}}\
\big{\|}{\cal
P}_{\check{\mathbf{u}}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\big{\|}}\to_{d}N\left(0,1\right),$

where $\widehat{b}_{k}^{(q)}$ is given by (13).

Theorem 4.5 is an immediate consequence of Slutsky’s Theorem and Theorems
4.2-4.4. It allows us to make inference or construct confidence intervals for
$\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$. Consider, for example,
testing hypothesis that

$H_{0}:\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle=0\qquad vs\qquad
H_{a}:\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\neq 0,$

when $d=o(\sqrt{n})$. We can proceed to reject $H_{0}$ if and only if

$\left|\sqrt{n}\left(\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle-\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle\right)\right|\geq
z_{\alpha/2}\sqrt{\frac{\widehat{\sigma}_{0}^{2}}{\widehat{\sigma}_{\pi(k)}^{2}}+\frac{\widehat{\sigma}_{0}^{4}}{\widehat{\sigma}_{\pi(k)}^{4}}}\big{\|}{\cal
P}_{\tilde{\mathbf{u}}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\big{\|},$

where $z_{\alpha/2}$ is the upper $\alpha/2$ quantile of the standard normal
distribution. Theorem 4.5 guarantees this is a level-$\alpha$ test
asymptotically. Similarly, we can also construct $(1-\alpha)$ confidence
interval for $\langle\mathbf{u}^{(q)}_{\pi(k)},\mathbf{v}\rangle$:

$\left(\langle\tilde{\mathbf{u}}^{(q)}_{k},\mathbf{v}\rangle\pm\frac{z_{\alpha/2}}{\sqrt{n}}\sqrt{\frac{\widehat{\sigma}_{0}^{2}}{\widehat{\sigma}_{\pi(k)}^{2}}+\frac{\widehat{\sigma}_{0}^{4}}{\widehat{\sigma}_{\pi(k)}^{4}}}\big{\|}{\cal
P}_{\tilde{\mathbf{u}}_{\pi(k)}^{(q)}}^{\perp}\mathbf{v}\big{\|}\right).$

In particular, by taking
$\mathbf{v}\in\\{\mathbf{e}_{1},\ldots,\mathbf{e}_{d_{q}}\\}$, we can use the
above formula to derive confidence intervals for the coordinates of
$\mathbf{u}_{\pi(k)}^{(q)}$. Situations with larger $d$ can also be treated
accordingly.

##  5 Numerical Experiments

To complement our theoretical analyses and further demonstrate the merits of
multiway PCA, we conducted several sets of numerical experiments.

###  5.1 Simulation Studies

We first present a set of simulation studies to illustrate the finite-sample
behavior of the sample PCs. These experiments are specifically designed to
assess the role of bias correction, and robustness to deviation from the
normal distribution. Throughout this subsection, unless otherwise noted,
samples were generated according to the spike covariance model (8) with $p=2$,
e.g., each $\mathscr{X}_{i}$ is a matrix. Since the two modes are
exchangeable, we only focus on the first mode $q=1$ for brevity. We also fixed
the number of spikes at $r=2$. In each case, we shall set the singular values
$\sigma_{1}=\sigma_{2}$. In other words, for each of our examples, the usual
PCA (with stringing) will not be able to identify the PCs because of the
multiplicity. As mentioned before, without loss of generality and for the sake
of brevity, we reordered $\check{\mathbf{u}}_{1}^{(1)}$ and
$\check{\mathbf{u}}_{2}^{(1)}$ such that
$\sin\angle(\check{\mathbf{u}}_{1}^{(1)},\mathbf{u}_{1})\leq\sin\angle(\check{\mathbf{u}}_{2}^{(1)},\mathbf{u}_{1})$.
In addition, we replaced $\check{\mathbf{u}}_{k}^{(1)}$ with
$-\check{\mathbf{u}}_{k}^{(1)}$ whenever
$\langle\check{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{k}^{(1)}\rangle<0$. For low-
dimensional setup, $\tilde{\mathbf{u}}_{1}^{(1)}$ and
$\tilde{\mathbf{u}}_{2}^{(1)}$ are treated similarly.

In the first set of experiments, we considered a low-dimensional setup with
$d_{1}=d_{2}=10$, $n=200$, $\sigma_{1}=\sigma_{2}=2$, and the true PCs were
given by

$\displaystyle\mathbf{u}_{1}^{(1)}=(\sqrt{3}/2,1/2,0,\dots,0)^{\top},\quad\mathbf{u}_{1}^{(2)}=(1,0,\dots,0)^{\top},$
$\displaystyle\mathbf{u}_{2}^{(1)}=(-1/2,\sqrt{3}/2,0,\dots,0)^{\top},\quad\mathbf{u}_{2}^{(2)}=(0,1,0,\dots,0)^{\top}.$
(14)

Figure 1(a) reports the histograms of the first two (nonzero) entries of
$\tilde{\mathbf{u}}_{1}^{(1)}$ based on 300 simulation runs. The histograms
are overlaid with the asymptotic distributions derived in Theorem 4.1. The
agreement between the two confirms the accuracy of the asymptotic distribution
when the dimensionality is low.

(a) $d_{1}=d_{2}=10$

(b) $d_{1}=d_{2}=50$

Figure 1: Multiway PCA for data generated from normal distribution.

To demonstrate the need and effectiveness of bias correction, we increased the
dimension to $d_{1}=d_{2}=50$. Correspondingly we set $n=400$ and
$\sigma_{1}=\sigma_{2}=3$. We repeated the experiment another 300 times and as
before, Figure 1(b) reports the histograms of the first two entries of
$\check{\mathbf{u}}_{1}^{(1)}$ along with the asymptotic distribution derived
in Theorem 4.3, plotted in red lines. The dashed black line overlaid with the
histogram of the first entries corresponds to the asymptotic distribution
without bias correction as given by Theorem 4.1. It is clear that in this
setting, debiasing is necessary and the bias correction of Theorem 4.3 indeed
leads to a more precise approximation of the finite sample distribution.

Our next set of simulations aims to explore the robustness of our approach to
deviation from normality. To this end, $\\{\theta_{k},k\in[r]\\}$ and the
entries of $\mathscr{E}$ were simulated independently from ${\rm
Poisson}(1)-1$ (so that they still have mean $0$ and variance $1$). Again we
set $n=400$ and $\sigma_{1}=\sigma_{2}=3$. Figures 2(a) and 2(b) summarize
results based on 300 runs, for dimensions $d_{1}=d_{2}=10$ and
$d_{1}=d_{2}=50$, respectively. We overlay them with the theoretical
asymptotic distributions given by Theorems 4.1 and 4.3. The results are
qualitatively similar to those from the previous setting.

(a) $d_{1}=d_{2}=10$

(b) $d_{1}=d_{2}=50$

Figure 2: Multiway PCA for data generated from Poisson distribution.

###  5.2 World Bank Data

We now consider a real-world data example – the open source global development
data from the World Bank***https://data.worldbank.org/. The world Bank offers
access to annual country-level data of a number of development indicators. In
particular, we shall focus on the following nine most common and important
economic and demographic indicators:

  * `GDP`: gross domestic product (GDP) based on purchasing power parity;

  * `Import`: import volume index (year 2000=100);

  * `Export`: export volume index (year 2000=100);

  * `CO2`: total CO2 emissions, in kilo-ton;

  * `CPI`: Consumer price index (year 2010 = 100);

  * `Life Span`: Life expectancy at birth;

  * `Urban Population`: Urban population, percentage of total population;

  * `Tourism`: number of international inbound tourists;

  * `Birth Rate`: Birth rate, crude (per 1,000 people).

Yearly data for these indicators have been recorded and we focus on data from
Year 2000 through 2018, as considerable data are missing outside this range.
We also discarded countries that have more than 5% of missing data in our
analysis, resulting in a total of 160 countries under consideration.

These indicators are all positive but of vastly different magnitudes. To this
end, a log transformation was first applied. Each log-transformed indicator
was then standardized so that the log-transformed indicator has a mean $0$ and
a mean absolute deviation $1$ for all countries. The use of mean absolute
deviation, instead of variance, for standardization allows more robust
analysis in the presence of outlying observations. Denote by
$\mathbf{X}_{k,t,i}$ the resulting indicator $i$ for country $k$ at time $t$.
There remain a handful of missing values and for convenience, they are
replaced with $0$ in our analysis. The data tensor $\mathbf{X}$ of dimensions
$160\times 19\times 9$. Each frontal slice

$\mathscr{X}_{k}=\mathbf{X}_{k,\cdot,\cdot}$

corresponds to a country and is a $19\times 9$ matrix. Note that its ambient
dimension is $19\times 9=171$ and greater than the number of countries so it
is problematic to apply the usual PCA with stringing. Accounting for the
multiway structure, we can consider the multiway PCs of the form

$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\mathbf{u}_{k}^{(2)}\in\mathbb{R}^{19\times
9}.$

These PC carry a clear meaning: each $\mathscr{U}_{k}$ represents a shared
_development pattern_ , where $\mathbf{u}_{k}^{(1)}$ is the corresponding
shared _temporal trend_ , and $\mathbf{u}_{k}^{(2)}$ is the corresponding
_comovement pattern_.

Figure 3 plots the estimated leading PC along both modes, namely
$\check{\mathbf{u}}_{1}^{(1)}$ and $\check{\mathbf{u}}_{k}^{(2)}$, together
with the $95\%$ confidence intervals for each of their coordinates. It is by
far the most significant component, explaining 57.6% of the total variation.
It is also evident from the temporal component that the first PC describes a
roughly constant growth trend. The only year with a decrease is 2008 when the
Global Financial Crisis took place. Correspondingly, except for the entry
corresponding to birth rate, all other entries of
$\check{\mathbf{u}}_{k}^{(2)}$ are positive. This suggests a general economic
development during this period, with the birth rate in decline.

Figure 3: The first PC, general economic development: $\mathbf{u}_{1}^{(1)}$
and $\mathbf{u}_{1}^{(2)}$ plotted with 95% confidence intervals.

Similarly, Figure 4 shows the second multiway PC in the two modes along with
their $95\%$ confidence bands. This PC captures a change of developmental
direction at the year of 2008. In particular, CPI, life span, urban
population, and tourism steadily decreased prior to 2008 but reversed course
after the financial crisis. In contrast, GDP, import, export, CO2 emission and
birth rate followed an opposite pattern. There are many plausible explanations
for this pattern. It is possible that the quantitative easing policies applied
by most major economies since 2008 led to growth in the domestic market, thus
enhancing the life-quality indicators. It is also possible that the growing
inequality after 2008, also caused by quantitative easing among other factors
(see, e.g., Montecino and Epstein, 2015), has in turn caused the increase in
life quality among the upper and the upper middle class. The tourism indicator
is the number of international inbound tourists, which most likely is driven
by the upper middle class and beyond. The continuous increase in life
expectancy in the USA is also reported to be driven primarily by the well-off
(see, e.g., Chetty et al., 2016).

Figure 4: The second PC, life quality: $\mathbf{u}_{2}^{(1)}$ and
$\mathbf{u}_{2}^{(2)}$ plotted with 95% confidence intervals.

Finally, Figure 5 shows the third multiway PC. We begin to see much wider
confidence intervals as the signal becomes weaker. In fact, only the period
around 2008 are significantly different from zero, and likewise, the entries
corresponding to life span, urbanization, and tourism are statistically
insignificant. This indicates that these patterns likely focus on the impact
of the 2008 financial crisis: it caused an immediate economic downturn but
recovered not long after.

Figure 5: The third PC, international trade: $\mathbf{u}_{3}^{(1)}$ and
$\mathbf{u}_{3}^{(2)}$ plotted with 95% confidence intervals.

###  5.3 NYC Bike Rental Data

Another data example we considered is the Citibike trip
data†††https://ride.citibikenyc.com/system-data. In particular, all the
Citibike trips from January 1, 2018 to December 31, 2019, on weekdays (522
days in total) that started in Manhattan and lasted for at least 60 seconds
were used in our analysis. During this period, there are 35 zip codes in
Manhattan with at least one Citibike station. There are a total of 29,515,527
trips and we form a data tensor $\mathbf{Y}$ of dimension $522\times 24\times
35$ where $Y_{kij}$ denotes the number of trips starting during the $i$th hour
of the $k$th day from the $j$th zip code.

The number of counts at different zip codes are of drastically different
magnitudes, and the total counts during the two years also display a clear
seasonal trend. To facilitate our analysis, we standardized the counts from
each zip code $j$ at each day $k$ so that they have mean $0$ and mean absolute
deviation $1$. As in the previous example, a direct application of the usual
PCA can be misleading as the ambient dimension of the daily observation is
$24\times 35=840$ and greater than the number of days. Nonetheless, it is
helpful to consider multiway PCs of the form

$\mathscr{U}_{k}=\mathbf{u}_{k}^{(1)}\otimes\mathbf{u}_{k}^{(2)}\in\mathbb{R}^{24\times
35},$

where $\mathbf{u}_{k}^{(1)}$ captures the time-of-the-day effect of bike
rental, and $\mathbf{u}_{k}^{(2)}$ the location pattern.

Figure 6 plots the first multiway PC. The spatial pattern clearly indicates
that this represents an overall pattern across Manhattan with all 35 entries
of $\mathbf{u}_{k}^{(2)}$ being estimated as positive. The temporal pattern
indicates that bike rental strongly coincides with the rush hours with two
peaks during the morning and afternoon rush hours. The blank area downtown is
zip code 10006, the big blank rectangular is Central Park, the small blank
underneath is zip code 10020, and the blank area to the north of Central Park
has zip codes 10030 and 10031. At the time of the recorded period, no Citibike
station existed in these areas.

Figure 6: The first PC: overall pattern.

The second PC, as shown in Figure 7, reveals differences in rental patterns
across neighborhoods. While the first PC suggests increased rental activities
both in the morning and afternoon rush hours, the second PC captures the
difference between morning and evening rental patterns as indicated by the
positive peak during the evening rush hours and the negative peak during the
morning rush hours. As such, a neighborhood with positive loadings may see
more evening rentals than morning rentals. These are the downtown Financial
District, Lower Manhattan, and Midtown, largely corresponding to the business
area of Manhattan. On the other hand, zip codes corresponding to negative
loadings represent mostly residential areas of Manhattan, including the East
Village, Upper West Side, and Upper East Side.

Figure 7: The second PC: rush hour differences.

Figure 8 depicts the third PC. The temporal pattern has a narrow and tall peak
during the afternoon rush hours suggesting that this PC captures the subtle
spatial difference during this time of the day. In particular, the zip codes
with large positive values (purple color) are the area around Wall Street (the
small purple block in Lower Manhattan), the area around Grand Central
Terminal, and an area in Upper East Side. The negative zip codes in this
pattern include the areas around SoHo, Greenwich Village, and Harlem.

Figure 8: The third PC: afternoon rush hour details.

##  6 Summary

In this paper, we study PCA under the settings that each observation is a
matrix or more generally a multiway array. We investigate how to extract
multiway PCs and study their statistical properties. In addition to the
obvious advantages of increased efficiency and enhanced interpretability, our
analysis provides a number of new insights into the operating characteristics
of multiway PCA and their methodological implications.

First, we show that multiway PCs can be estimated without the eigengap
requirement. Specifically, under a spike covariance model, we establish rates
of convergence for the sample multiway PCs. In particular, they are consistent
whenever the signal-to-noise ratio
$\sigma_{k}/\sigma_{0}\gg\max\\{d/n,(d/n)^{1/4}\\}$ where $d$ is the dimension
of one mode. Perhaps more interestingly, we prove that the sample multiway PCs
are asymptotically independent of each other, at least when the dimension
$d=o(\sqrt{n})$. In higher dimensions, the sample PCs can be biased and the
bias can be corrected via sample-splitting to lead to asymptotically normal
estimates of the multiway PCs, which enables us to construct confidence
intervals or conduct hypothesis testing for linear forms of the PCs.

Our theoretical developments are complemented by numerical experiments, both
simulated and real. In particular, meaningful findings can be inferred when
applying our methods to two real-world datasets, further demonstrating the
merits of our methodology.

## References

  * Anandkumar et al. (2014) 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:2773–2832, 2014\.
  * Anderson (1984) T. W. Anderson.  _An Introduction to Multivariate Statistical Analysis_.  Wiley, New York, NY, second edition, 1984.
  * Auddy and Yuan (2020) Arnab Auddy and Ming Yuan.  Perturbation bounds for (nearly) orthogonally decomposable tensors.  _arXiv preprint arXiv:2007.09024_ , 2020.
  * Bai and Silverstein (2010) Zhidong Bai and Jack W Silverstein.  _Spectral analysis of large dimensional random matrices_ , volume 20.  Springer, 2010.
  * Bai and Yao (2012) Zhidong Bai and Jianfeng Yao.  On sample eigenvalues in a generalized spiked population model.  _Journal of Multivariate Analysis_ , 106:167–177, 2012\.
  * Baik and Silverstein (2006) Jinho Baik and Jack W Silverstein.  Eigenvalues of large sample covariance matrices of spiked population models.  _Journal of multivariate analysis_ , 97(6):1382–1408, 2006.
  * Benaych-Georges and Nadakuditi (2011) Florent Benaych-Georges and Raj Rao Nadakuditi.  The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices.  _Advances in Mathematics_ , 227(1):494–521, 2011\.
  * Bi et al. (2021) Xuan Bi, Xiwei Tang, Yubai Yuan, Yanqing Zhang, and Annie Qu.  Tensors in statistics.  _Annual review of statistics and its application_ , 8:345–368, 2021.
  * Birnbaum et al. (2013) Aharon Birnbaum, Iain M Johnstone, Boaz Nadler, and Debashis Paul.  Minimax bounds for sparse pca with noisy high-dimensional data.  _Annals of statistics_ , 41(3):1055, 2013.
  * Chen et al. (2020a) Elynn Y Chen, Jianqing Fan, and Ellen Li.  Statistical inference for high-dimensional matrix-variate factor model.  _arXiv preprint arXiv:2001.01890_ , 2020a.
  * Chen et al. (2020b) Elynn Y Chen, Dong Xia, Chencheng Cai, and Jianqing Fan.  Semiparametric tensor factor analysis by iteratively projected svd.  _arXiv preprint arXiv:2007.02404_ , 2020b.
  * Chen et al. (2021) Rong Chen, Dan Yang, and Cun-Hui Zhang.  Factor models for high-dimensional tensor time series.  _Journal of the American Statistical Association_ , pages 1–23, 2021\.
  * Chetty et al. (2016) Raj Chetty, Michael Stepner, Sarah Abraham, Shelby Lin, Benjamin Scuderi, Nicholas Turner, Augustin Bergeron, and David Cutler.  The association between income and life expectancy in the united states, 2001-2014.  _Jama_ , 315(16):1750–1766, 2016.
  * Cichocki et al. (2015) Andrzej Cichocki, Danilo Mandic, Lieven De Lathauwer, Guoxu Zhou, Qibin Zhao, Cesar Caiafa, and Huy Anh Phan.  Tensor decompositions for signal processing applications: From two-way to multiway component analysis.  _IEEE signal processing magazine_ , 32(2):145–163, 2015.
  * De Lathauwer et al. (2000) 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.
  * Friedland (2013) Shmuel Friedland.  Best rank one approximation of real symmetric tensors can be chosen symmetric.  _Frontiers of Mathematics in China_ , 8(1):19–40, 2013.
  * Hackbusch (2012) Wolfgang Hackbusch.  _Tensor spaces and numerical tensor calculus_ , volume 42.  Springer, 2012.
  * Han et al. (2022) Rungang Han, Rebecca Willett, and Anru R Zhang.  An optimal statistical and computational framework for generalized tensor estimation.  _The Annals of Statistics_ , 50(1):1–29, 2022\.
  * Han et al. (2020) Yuefeng Han, Rong Chen, Dan Yang, and Cun-Hui Zhang.  Tensor factor model estimation by iterative projection.  _arXiv preprint arXiv:2006.02611_ , 2020.
  * Harshman and Lundy (1984) Richard A Harshman and Margaret E Lundy.  The parafac model for three-way factor analysis and multidimensional scaling.  _Research methods for multimode data analysis_ , 46:122–215, 1984.
  * Hopkins et al. (2015) Samuel B Hopkins, Jonathan Shi, and David Steurer.  Tensor principal component analysis via sum-of-square proofs.  In _Conference on Learning Theory_ , pages 956–1006, 2015.
  * Janzamin et al. (2019) Majid Janzamin, Rong Ge, Jean Kossaifi, Anima Anandkumar, et al.  Spectral learning on matrices and tensors.  _Foundations and Trends® in Machine Learning_ , 12(5-6):393–536, 2019.
  * Johnstone (2001) Iain M Johnstone.  On the distribution of the largest eigenvalue in principal components analysis.  _Annals of statistics_ , pages 295–327, 2001.
  * Johnstone and Lu (2009) Iain M Johnstone and Arthur Yu Lu.  On consistency and sparsity for principal components analysis in high dimensions.  _Journal of the American Statistical Association_ , 104(486):682–693, 2009.
  * Jolliffe (2002) I. Jolliffe.  _Principal Component Analysis_.  Springer, 2002.
  * Jung and Marron (2009) Sungkyu Jung and J Stephen Marron.  Pca consistency in high dimension, low sample size context.  _The Annals of Statistics_ , 37(6B):4104–4130, 2009.
  * Kolda (2001) Tamara G Kolda.  Orthogonal tensor decompositions.  _SIAM Journal on Matrix Analysis and Applications_ , 23(1):243–255, 2001.
  * Koltchinskii and Lounici (2014) Vladimir Koltchinskii and Karim Lounici.  Asymptotics and concentration bounds for spectral projectors of sample covariance.  _arXiv preprint arXiv:1408.4643_ , 2014.
  * Koltchinskii et al. (2017) Vladimir Koltchinskii, Karim Lounici, et al.  Concentration inequalities and moment bounds for sample covariance operators.  _Bernoulli_ , 23(1):110–133, 2017.
  * Koltchinskii et al. (2020) Vladimir Koltchinskii, Matthias Löffler, and Richard Nickl.  Efficient estimation of linear functionals of principal components.  _The Annals of Statistics_ , 48(1):464–490, 2020\.
  * Kong et al. (2005) Hui Kong, Lei Wang, Eam Khwang Teoh, Xuchun Li, Jian-Gang Wang, and Ronda Venkateswarlu.  Generalized 2d principal component analysis for face image representation and recognition.  _Neural Networks_ , 18(5-6):585–594, 2005.
  * Kroonenberg (2008) Pieter M Kroonenberg.  _Applied multiway data analysis_.  John Wiley & Sons, 2008.
  * Kroonenberg and De Leeuw (1980) Pieter M Kroonenberg and Jan De Leeuw.  Principal component analysis of three-mode data by means of alternating least squares algorithms.  _Psychometrika_ , 45(1):69–97, 1980.
  * Lee et al. (2010) Seunggeun Lee, Fei Zou, and Fred A Wright.  Convergence and prediction of principal component scores in high-dimensional settings.  _Annals of statistics_ , 38(6):3605, 2010.
  * Li et al. (2010) Xuelong Li, Yanwei Pang, and Yuan Yuan.  L1-norm-based 2dpca.  _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , 40(4):1170–1175, 2010.
  * Liu et al. (2017) Tianqi Liu, Ming Yuan, and Hongyu Zhao.  Characterizing spatiotemporal transcriptome of human brain via low rank tensor decomposition.  _arXiv preprint arXiv:1702.07449_ , 2017.
  * Lu et al. (2006) Haiping Lu, Konstantinos N Plataniotis, and Anastasios N Venetsanopoulos.  Multilinear principal component analysis of tensor objects for recognition.  In _18th International Conference on Pattern Recognition (ICPR’06)_ , volume 2, pages 776–779. IEEE, 2006.
  * Lu et al. (2008) Haiping Lu, Konstantinos N Plataniotis, and Anastasios N Venetsanopoulos.  Mpca: Multilinear principal component analysis of tensor objects.  _IEEE transactions on Neural Networks_ , 19(1):18–39, 2008.
  * Lu et al. (2011) Haiping Lu, Konstantinos N Plataniotis, and Anastasios N Venetsanopoulos.  A survey of multilinear subspace learning for tensor data.  _Pattern Recognition_ , 44(7):1540–1551, 2011\.
  * Montecino and Epstein (2015) Juan Montecino and Gerald Epstein.  Did quantitative easing increase income inequality?  _Institute for New Economic Thinking working paper series_ , (28), 2015.
  * Nadler (2008) Boaz Nadler.  Finite sample approximation results for principal component analysis: A matrix perturbation approach.  _The Annals of Statistics_ , 36(6):2791–2817, 2008.
  * Paul (2007) Debashis Paul.  Asymptotics of sample eigenstructure for a large dimensional spiked covariance model.  _Statistica Sinica_ , pages 1617–1642, 2007.
  * Richard and Montanari (2014) Emile Richard and Andrea Montanari.  A statistical model for tensor pca.  In _Advances in Neural Information Processing Systems_ , pages 2897–2905, 2014.
  * Shen et al. (2013) Dan Shen, Haipeng Shen, Hongtu Zhu, and JS Marron.  Surprising asymptotic conical structure in critical sample eigen-directions.  _arXiv preprint arXiv:1303.6171_ , 2013.
  * Taguchi (2018) Y-H Taguchi.  Tensor decomposition-based and principal-component-analysis-based unsupervised feature extraction applied to the gene expression and methylation profiles in the brains of social insects with multiple castes.  _BMC bioinformatics_ , 19(4):99, 2018.
  * Vasilescu and Terzopoulos (2002) M Alex O Vasilescu and Demetri Terzopoulos.  Multilinear analysis of image ensembles: Tensorfaces.  In _European conference on computer vision_ , pages 447–460. Springer, 2002.
  * Vershynin (2010) Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Wang and Fan (2017) Weichen Wang and Jianqing Fan.  Asymptotics of empirical eigenstructure for high dimensional spiked covariance.  _Annals of statistics_ , 45(3):1342, 2017.
  * Xia et al. (2020) Dong Xia, Anru R Zhang, and Yuchen Zhou.  Inference for low-rank tensors–no need to debias.  _arXiv preprint arXiv:2012.14844_ , 2020.
  * Yang et al. (2004) Jian Yang, David Zhang, Alejandro F Frangi, and Jing-yu Yang.  Two-dimensional pca: a new approach to appearance-based face representation and recognition.  _IEEE transactions on pattern analysis and machine intelligence_ , 26(1):131–137, 2004.
  * Zhang and Xia (2018) Anru Zhang and Dong Xia.  Tensor svd: Statistical and computational limits.  _IEEE Transactions on Information Theory_ , 64(11):7311–7338, 2018.
  * Zhang and Zhou (2005) Daoqiang Zhang and Zhi-Hua Zhou.  (2d) 2pca: Two-directional two-dimensional pca for efficient face representation and recognition.  _Neurocomputing_ , 69(1-3):224–231, 2005.
  * Zhang and Golub (2001) Tong Zhang and Gene H Golub.  Rank-one approximation to high order tensors.  _SIAM Journal on Matrix Analysis and Applications_ , 23(2):534–550, 2001.

##  Appendix A Notations and Preliminary Bounds

Write $a\vee b:=\max\\{a,b\\}$ and $a\wedge b:=\min\\{a,b\\}$. For a positive
integer $n$, let $[n]:=\\{1,2,\dots,n\\}$. For a vector
$\mathbf{x}\in\mathbb{R}^{d}$, denote $\|\mathbf{x}\|$ to be its
$\ell_{2}$-norm, $\|\mathbf{x}\|_{1}$ to be its $\ell_{1}$-norm, and
$\|\mathbf{x}\|_{\infty}=\max_{i}|x_{i}|$ to be its $\ell_{\infty}$-norm. For
two sequences of real numbers $\\{a_{n}\\}$ and $\\{b_{n}\\}$, write
$a_{n}=O(b_{n})$ if $\exists C,\exists M$, such that $\forall n>M$,
$|a_{n}|\leq C|b_{n}|$. Write $a_{n}=o(b_{n})$ if
$\lim_{n\to\infty}a_{n}/b_{n}=0$. For two sequences of real-valued random
variables $X_{n}$ and $Y_{n}$, write $X_{n}=O_{p}(Y_{n})$ if
$X_{n}=R_{n}Y_{n}$ and $R_{n}$ is uniformly tight. Write $X_{n}=o_{p}(Y_{n})$
if $X_{n}=R_{n}Y_{n}$ and $R_{n}\overset{p}{\to}0$. For linear subspace $U$ of
$\mathbb{R}^{d}$, denote $P_{U}$ and $P_{U}^{\perp}$ to be the orthogonal
projection onto $U$ and its orthogonal complement $U^{\perp}$, respectively.
For a non-zero vector $u\in\mathbb{R}^{d}$, denote $P_{u}:=P_{{\rm
span}\\{u\\}}$ and $P_{u}^{\perp}:=P_{{\rm span}\\{u\\}}^{\perp}$.

For an order-$k$ tensor $\mathscr{T}\in\mathbb{R}^{d_{1}\times
d_{2}\times\dots\times d_{k}}$, define its _tensor operator norm_ as:

$\displaystyle\|\mathscr{T}\|:=\sup_{\mathbf{u}_{j}\in\mathbb{R}^{d_{j}},\|\mathbf{u}_{j}\|=1}\mathscr{T}(\mathbf{u}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{k}).$
(15)

Specifically, when $k=2$ so that $\mathscr{T}\in\mathbb{R}^{d_{1}\times
d_{2}}$ is a matrix, $\|\mathscr{T}\|$ is the matrix spectral norm of
$\mathscr{T}$. For tensor $\mathscr{T}\in\mathbb{R}^{d_{1}\times
d_{2}\times\dots\times d_{k}}$, write

$\|\mathscr{T}\|_{\max}=\max_{i_{1},\dots,i_{k}}\left|\mathscr{T}_{i_{1},\dots,i_{k}}\right|$

to be its $\ell_{\infty}$-norm.

With a slight abuse of notation, the mode $q$ product of
$\mathscr{T}\in\mathbb{R}^{d_{1}\times d_{2}\times\dots\times d_{k}}$ with a
vector $\mathbf{a}\in\mathbb{R}^{d_{q}}$, denoted by
$\mathscr{T}\times_{q}\mathbf{a}$, is defined as an order-$(k-1)$ tensor of
size $d_{1}\times\dots d_{q-1}\times d_{q+1}\dots\times d_{k}$, with elements

$[\mathscr{T}\times_{q}\mathbf{a}]_{i_{1}\dots i_{q-1}i_{q+1}\dots
i_{k}}=\sum_{i_{q}=1}^{d_{q}}\mathscr{T}_{i_{1}\dots i_{q}\dots
i_{k}}\mathbf{a}_{i_{q}}.$

Write

$\widehat{\Sigma}_{\theta}={1\over
n}\sum_{i=1}^{n}\theta_{i}\otimes\theta_{i},\qquad\widehat{\Sigma}_{\mathscr{E}}={1\over
n}\sum_{i=1}^{n}\mathscr{E}_{i}\otimes\mathscr{E}_{i},$

and

$\widehat{\Sigma}_{\theta,\mathscr{E}}={1\over
n}\sum_{i=1}^{n}\theta_{i}\otimes\mathscr{E}_{i},$

the sample covariance matrices of $\theta$, $\mathscr{E}$ and between them
respectively. Correspondingly denote by $\Sigma_{\theta}$,
$\Sigma_{\mathscr{E}}$ and $\Sigma_{\theta,\mathscr{E}}$ their population
counterpart. It is clear $\Sigma_{\theta,\mathscr{E}}=0$. Recall also that

$\widehat{\Sigma}={1\over
n}\sum_{i=1}^{n}\mathscr{X}_{i}\otimes\mathscr{X}_{i}.$

and

$\Sigma=\sum_{l=1}^{r}\sigma_{l}\mathbf{u}_{l}^{(1)}\otimes\cdots\otimes\mathbf{u}_{l}^{(p)}\otimes\mathbf{u}_{l}^{(1)}\otimes\cdots\otimes\mathbf{u}_{l}^{(p)}+\sigma_{0}^{2}\mathscr{I}$

are the sample and population covariance matrices of $\mathscr{X}$.

The proof relies on the following technical lemmas.

######  Lemma 1.

There exists a numerical constant $C>0$ such that for any $t\geq 1$,

$\|\widehat{\Sigma}-\Sigma\|\leq
C(\sigma_{1}^{2}+\sigma_{0}^{2})\max\left\\{\sqrt{d\over n},{d\over
n},\sqrt{t\over n},{t\over n}\right\\},$
$\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\leq
C\sigma_{0}\max\left\\{\sqrt{d\over n},{d\over n},\sqrt{t\over n},{t\over
n}\right\\},$ $\|\widehat{\Sigma}_{\mathscr{E}}-\Sigma_{\mathscr{E}}\|\leq
C\sigma_{0}^{2}\max\left\\{\sqrt{d\over n},{d\over n},\sqrt{t\over n},{t\over
n}\right\\},$

and

$\|\widehat{\Sigma}_{\theta}-\Sigma_{\theta}\|\leq C\max\left\\{\sqrt{r\over
n},{r\over n},\sqrt{t\over n},{t\over n}\right\\},$

with probability at least $1-e^{-t}$.

Note that we shall use $C$ to denote a constant that may take different values
at each appearance. We shall also make use the following bounds:

######  Lemma 2.

There exists a numerical constant $C>0$ such that for any $t>0$,

$\|\widehat{\Sigma}_{\theta}-\Sigma_{\theta}\|_{\max}\leq
C\max\left\\{\sqrt{\log r\over n},{\log r\over n},\sqrt{t\over n},{t\over
n}\right\\},$ $\max_{1\leq l_{1},l_{2}\leq
r}\left|\widehat{\Sigma}_{\mathscr{E},\theta}(\mathbf{u}_{l_{1}}^{(1)},\mathbf{u}_{l_{2}}^{(2)},\ldots,\mathbf{u}_{l_{2}}^{(p)},{\bf
e}_{l_{2}})\right|\leq C\sigma_{0}\max\left\\{\sqrt{\log r\over n},{\log
r\over n},\sqrt{t\over n},{t\over n}\right\\}$

where ${\bf e}_{l_{2}}$ is the $l_{2}$th canonical basis of
${\mathbb{R}}^{r}$, and

$\max_{1\leq l\leq
r}\left|(\widehat{\Sigma}_{\mathscr{E}}-\Sigma_{\mathscr{E}})(\mathbf{u}_{l}^{(1)},\ldots,\mathbf{u}_{l}^{(p)},\mathbf{u}_{l}^{(1)},\ldots,\mathbf{u}_{l}^{(p)})\right|\leq
C\sigma_{0}^{2}\max\left\\{\sqrt{\log r\over n},{\log r\over n},\sqrt{t\over
n},{t\over n}\right\\},$

with probability at least $1-e^{-t}$.

Both Lemmas are well known and follow immediately from an application of union
bounds and $\chi^{2}$ tail bounds. See, e.g., Vershynin (2010).

##  Appendix B Proof of Theorems 3.1 and 3.2

Theorem 3.1 follows immediately from Theorem 3.2 and it suffices to prove the
latter. For brevity, we shall focus on the case when $d\leq n$ and $r$
diverges with $n$. Denote by ${\cal E}$ the event that

$\|\widehat{\Sigma}_{\theta}-\Sigma_{\theta}\|\leq C\sqrt{r\over n},\quad{\rm
and}\quad\sigma_{0}^{-2}\|\widehat{\Sigma}_{\mathscr{E}}-\Sigma_{\mathscr{E}}\|,\sigma_{0}^{-1}\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\leq
C\sqrt{d\over n}$

and

$\|\widehat{\Sigma}_{\theta}-\Sigma_{\theta}\|_{\max},\
\sigma_{0}^{-1}\max_{1\leq l_{1},l_{2}\leq
r}\left|\widehat{\Sigma}_{\mathscr{E},\theta}(\mathbf{u}_{l_{1}}^{(1)},\mathbf{u}_{l_{2}}^{(2)},\ldots,\mathbf{u}_{l_{2}}^{(p)},{\bf
e}_{l_{2}})\right|\leq C\sqrt{\log r\over n}$

By Lemmas 1 and 2, ${\cal E}$ holds with probability tending to one. It
suffices to proceed conditional on the event ${\cal E}$.

As noted, the $k$th sample PCs may not correspond to the $k$th population PCs
because we do not assume the existence of eigengap and $\sigma_{k}$s may not
even be distinct. Nonetheless, we can match the sample PCs with population PCs
as follows. Define

$\pi(1)=\operatorname{argmax}_{1\leq l\leq
r}\left\\{\sigma_{l}^{2}\left|\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{1}^{(q)}\rangle\right|\right\\},$

and for $k>1$,

$\pi(k):=\operatorname{argmax}_{l\notin\pi([k-1])}\left\\{\sigma_{l}^{2}\left|\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right|\right\\}.$

The goal is to show that with high probability,

$\eta_{k}:=\max_{1\leq q\leq
p}\sin\angle(\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)})\leq
C\left({\sigma_{0}\over\sigma_{\pi(k)}}+{\sigma_{0}^{2}\over\sigma_{\pi(k)}^{2}}\right)\max\left\\{\sqrt{d\over
n},{d\over n}\right\\}=:\delta_{k},$ (16)

for $k=1,\ldots,r$. Our proof proceeds by induction over $k$. To facilitate
the induction, we shall also prove that

$\displaystyle\tilde{\eta}_{k}$ $\displaystyle:=$ $\displaystyle\max_{1\leq
q\leq
p}\max_{l\notin\pi([k])}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle$
(17) $\displaystyle\leq$ $\displaystyle
C\left({\sigma_{0}\over\sigma_{\pi(k)}}+{\sigma_{0}^{2}\over\sigma_{\pi(k)}^{2}}\right)^{2}\max\left\\{{d\over
n},{d^{2}\over
n^{2}}\right\\}+C\left({\sigma_{0}\over\sigma_{\pi(k)}}+{\sigma_{0}^{2}\over\sigma_{\pi(k)}^{2}}\right)\sqrt{\log
r\over n}$ $\displaystyle=:$ $\displaystyle\tilde{\delta}_{k}.$

In addition to (16) and (17), we shall also prove that

$\sigma_{\pi(k)}^{2}\geq\max_{l\notin\pi([k])}\sigma_{l}^{2}(1-C\delta_{k}^{2}).$
(18)

This immediately implies that

$\sum_{l=1}^{k}\tilde{\delta}_{k}^{2}\leq C\delta_{k}^{2}\quad{\rm
and}\quad\max_{1\leq l\leq k}\\{\sigma_{\pi(l)}\delta_{l}\\}\leq
C\sigma_{\pi(k)}\delta_{k},$

by the taking the constant $c_{0}$ in (10) small enough. We shall make use of
these bounds repeatedly.

As noted, we shall proceed by induction over $k$. In particular, we shall
denote by $\delta_{0}=\tilde{\delta}_{0}=0$ so that the the base case holds
trivially when $k=0$. Now assume the induction hypotheses (16) and (17) holds
for $1,\ldots,k-1$. We want to show that they continue to hold for $k$. The
general architect of the argument is similar to that for the base case, but
additional challenges arise with the need to control the impact of estimation
error of $\widehat{\mathbf{u}}_{l}^{(q)}$s for $1\leq l<k$.

Denote by $\widehat{{\cal P}}^{(q)}$ the projection matrix onto the linear
space spanned by
$\\{\widehat{\mathbf{u}}_{1}^{(q)},\dots,\widehat{\mathbf{u}}_{k-1}^{(q)}\\}$,
for $q\in[p]$. Note that in the case when $k=1$, $\widehat{{\cal
P}}^{(q)}={\bf 0}_{d\times d}$. Then

$\displaystyle(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})$
$\displaystyle=$
$\displaystyle\operatorname{argmax}_{\begin{subarray}{c}\|\mathbf{w}^{(q)}\|\leq
1,\langle\mathbf{w}^{(q)},\widehat{\mathbf{u}}_{l}^{(q)}\rangle=0,\\\ \forall
l\leq
k-1,q\in[p]\end{subarray}}\widehat{\Sigma}(\mathbf{w}^{(1)},\dots,\mathbf{w}^{(p)},\mathbf{w}^{(1)},\dots,\mathbf{w}^{(p)})$
$\displaystyle=$
$\displaystyle\operatorname{argmax}_{\begin{subarray}{c}\|\mathbf{w}^{(q)}\|=1,\forall
q\in[p]\end{subarray}}\widehat{\Sigma}(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)})$

where $\widehat{{\cal P}}^{(q)}_{\perp}=I-\widehat{{\cal P}}^{(q)}$. Observe
that

$\widehat{{\cal
P}}_{\perp}^{(q)}\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\cdot,\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\propto\widehat{\mathbf{u}}_{k}^{(q)},$

where $\tilde{\Sigma}=\widehat{\Sigma}-\sigma_{0}^{2}\mathscr{I}$. This
implies that

$\langle\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{w}\rangle={\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\over\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})}.$

In particular, for $l\neq k$,

$\sin\angle(\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{u}_{l}^{(q)})={\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\over\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})},$

and

$\sin\angle(\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{u}_{k}^{(q)})={\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(q-1)},\widehat{{\cal
P}}_{\perp}^{(q)}(\widehat{\mathbf{u}}_{k}^{(q)}-\langle\widehat{\mathbf{u}}_{k}^{(q)},\mathbf{u}_{\pi(k)}^{(q)}\rangle\mathbf{u}_{\pi(k)}^{(q)}),\widehat{\mathbf{u}}_{k}^{(q+1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\over\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})},$

We shall derive lower bounds for the nominators and an upper bound for the
denominator. It suffices to consider the case $q=1$. Other indices can be
treated in an identical fashion.

###  B.1 Lower Bound for
$\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})$

Denote by $\tilde{\Sigma}=\Sigma-\sigma_{0}^{2}\mathscr{I}$. Observe that

$\displaystyle\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})$
(19) $\displaystyle\geq$
$\displaystyle\max_{l\notin\pi([k-1])}\tilde{\Sigma}(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)})$ $\displaystyle\geq$
$\displaystyle\max_{l\notin\pi([k-1])}(\Sigma-\sigma_{0}^{2}\mathscr{I})\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)})$
$\displaystyle-\sup_{\|\mathbf{w}^{(q)}\|\leq 1,1\leq q\leq
p}\left|(\widehat{\Sigma}-\Sigma)(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})\right|.$

Next we bound the two terms on the rightmost hand side.

Starting with the first term, note that for any $l\notin\pi([k-1])$,

$\displaystyle(\Sigma-\sigma_{0}^{2}\mathscr{I})(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)})$ $\displaystyle=$
$\displaystyle\sum_{1\leq l^{\prime}\leq
r}\left[\sigma_{l^{\prime}}^{2}\prod_{q=1}^{p}\langle\widehat{{\cal
P}}^{(q)}_{\perp}\mathbf{u}_{l}^{(q)},\mathbf{u}_{l^{\prime}}^{(q)}\rangle^{2}\right]$
$\displaystyle\geq$
$\displaystyle\sigma_{l}^{2}\prod_{q=1}^{p}\langle\widehat{{\cal
P}}^{(q)}_{\perp}\mathbf{u}_{l}^{(q)},\mathbf{u}_{l}^{(q)}\rangle^{2}$
$\displaystyle=$ $\displaystyle\sigma_{l}^{2}\prod_{q=1}^{p}\|\widehat{{\cal
P}}^{(q)}_{\perp}\mathbf{u}_{l}^{(q)}\|^{4}.$

By the induction hypothesis (17),

$\|\widehat{{\cal P}}^{(q)}\mathbf{u}_{l}^{(q)}\|^{2}=\sum_{1\leq
l^{\prime}<k}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{l^{\prime}}^{(q)}\rangle^{2}\leq\sum_{1\leq
l^{\prime}<k}\tilde{\delta}^{2}_{l^{\prime}}\leq C\delta_{k}^{2}.$

By taking the constant $c_{0}$ of (10) small enough, we can ensure that

$\max_{l\notin\pi([k-1])}(\Sigma-\sigma_{0}^{2}\mathscr{I})(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{u}_{l}^{(1)},\ldots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{u}_{l}^{(p)})\geq(1-C\delta_{k}^{4})\tau^{2}.$ (20)

where

$\tau^{2}=\max_{l\notin\pi([k-1])}\sigma_{l}^{2}.$

Next we derive a bound for

$\sup_{\|\mathbf{w}^{(q)}\|\leq 1,1\leq q\leq
p}(\Sigma-\widehat{\Sigma})(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)}).$

Note that

$\displaystyle(\Sigma-\widehat{\Sigma})(\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)},\widehat{{\cal
P}}^{(1)}_{\perp}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}^{(p)}_{\perp}\mathbf{w}^{(p)})$ (21) $\displaystyle=$
$\displaystyle\sum_{l_{1}=1}^{r}\sum_{l_{2}=1}^{r}\sigma_{l_{1}}\sigma_{l_{2}}\left(\widehat{\Sigma}_{\theta,l_{1}l_{2}}-\Sigma_{\theta,l_{1}l_{2}}\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)$
$\displaystyle+\frac{2}{n}\sum_{i=1}^{n}\sum_{l=1}^{r}\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})$
$\displaystyle+\left(\frac{1}{n}\sum_{i=1}^{n}[\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})]^{2}-\sigma_{0}^{2}\prod_{q=1}^{p}\|\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\|^{2}\right).$

We bound the three terms on the right hand side separately.

The first term can be bounded by

$\displaystyle\biggl{|}\sum_{l_{1}=1}^{r}\sum_{l_{2}=1}^{r}\sigma_{l_{1}}\sigma_{l_{2}}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il_{1}}\theta_{il_{2}}\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)$
$\displaystyle\qquad-\sum_{l=1}^{r}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}\biggr{|}$
$\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta}-I_{r}\|\left[\sum_{l=1}^{r}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}\right]$ $\displaystyle=$
$\displaystyle\|\widehat{\Sigma}_{\theta}-I_{r}\|\left[\sum_{l\in\pi([k-1])}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}+\sum_{l\notin\pi([k-1])}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}\right],$

Recall that for any $l\in\pi([k-1])$,

$|\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle|=|\langle\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{u}_{l}^{(q)},\mathbf{w}^{(q)}\rangle|\leq\|\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{u}_{l}^{(q)}\|\leq\delta_{l}.$

Therefore,

$\displaystyle\sum_{l\in\pi([k-1])}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}$ $\displaystyle\leq$
$\displaystyle\max_{l\in\pi([k-1])}\sigma_{l}^{2}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}$ $\displaystyle\leq$
$\displaystyle\max_{1\leq
l<k}\\{\sigma_{\pi(l)}^{2}\delta_{l}^{2p-2}\\}\leq\max_{1\leq
l<k}\\{\sigma_{\pi(l)}^{2}\delta_{l}^{2}\\}\leq\tau^{2},$

by taking $c_{0}$ of (10) small enough. On the other hand,

$\sum_{l\notin\pi([k-1])}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}\leq\max_{l\notin\pi([k-1])}\sigma_{l}^{2}=\tau^{2}.$

This implies that

$\left|\sum_{l_{1}=1}^{r}\sum_{l_{2}=1}^{r}\sigma_{l_{1}}\sigma_{l_{2}}\left(\widehat{\Sigma}_{\theta,l_{1}l_{2}}-\Sigma_{\theta,l_{1}l_{2}}\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\right|\leq
C\tau^{2}\sqrt{r\over n}.$ (22)

Similarly, the second term can be bounded by

$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}\sum_{l=1}^{r}\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})\right|$ (23) $\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\left[\sum_{l=1}^{r}\sigma_{l}^{2}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\rangle\right)^{2}\right]^{1/2}$
$\displaystyle\leq$ $\displaystyle C\tau\sigma_{0}\sqrt{d\over n}.$

Finally, the third term can be bounded by

$\left|\frac{1}{n}\sum_{i=1}^{n}[\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})]^{2}-\sigma_{0}^{2}\prod_{q=1}^{p}\|\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{w}^{(q)}\|^{2}\right|\leq\|\widehat{\Sigma}_{\mathscr{E}}-\Sigma_{\mathscr{E}}\|\leq
C\sigma_{0}^{2}\sqrt{d\over n}.$ (24)

Combing (21)-(24), we get

$\displaystyle\sup_{\|\mathbf{w}^{(q)}\|\leq 1,1\leq q\leq
p}\left|(\widehat{\Sigma}-\Sigma)(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}^{(1)},\dots,\widehat{{\cal
P}}_{\perp}^{(p)}\mathbf{w}^{(p)})\right|$ $\displaystyle\leq$ $\displaystyle
C\tau^{2}\sqrt{r\over n}+C(\sigma_{0}\tau+\sigma_{0}^{2})\sqrt{d\over n}.$

Together with (19), this implies

$\tilde{\Sigma}(\widehat{\mathbf{u}}_{k}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\geq\tau^{2}\left(1-C\delta_{k}^{4}-C\sqrt{r\over
n}\right)-C(\sigma_{0}\tau+\sigma_{0}^{2})\sqrt{d\over
n}-C\sigma_{0}^{2}\delta_{k}^{2},$ (25)

by taking $c_{0}$ of (10) small enough.

###  B.2 Upper Bounds for $\tilde{\Sigma}(\widehat{{\cal
P}}_{\perp}^{(1)}(\widehat{\mathbf{u}}_{k}^{(1)}-\langle\widehat{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{\pi(k)}^{(1)}\rangle\mathbf{u}_{\pi(k)}^{(1)}),\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})$

Observe that for any $\mathbf{w}$ orthogonal to $\mathbf{u}_{\pi(k)}^{(1)}$,
we have

$\displaystyle\tilde{\Sigma}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})$
(26) $\displaystyle=$
$\displaystyle\sigma_{\pi(k)}^{2}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}^{2}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle$
$\displaystyle+\sum_{l\neq\pi(k)}\sigma_{\pi(k)}\sigma_{l}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}\theta_{il}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle$
$\displaystyle+\sum_{l\neq\pi(k)}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle$
$\displaystyle+\sum_{l_{1},l_{2}\neq\pi(k)}\sigma_{l_{1}}\sigma_{l_{2}}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il_{1}}\theta_{il_{2}}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l_{1}}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l=1}^{r}\sigma_{l}\theta_{il}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l=1}^{r}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})-\sigma_{0}^{2}\langle\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(p)}\rangle.$

Again each term on the right hand side needs to be bounded carefully.

The first term on the right hand side of (26) can be bounded by

$\displaystyle\left|\sigma_{\pi(k)}^{2}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}^{2}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right|$
$\displaystyle\leq\tau^{2}\|\widehat{\Sigma}_{\theta}\|_{\max}|\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}\mathbf{w}\rangle|.$

In particular, when
$\mathbf{w}=\widehat{\mathbf{u}}_{k}^{(1)}-\langle\widehat{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{\pi(k)}^{(1)}\rangle\mathbf{u}_{\pi(k)}^{(1)}$,
we have

$\displaystyle|\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle|$ $\displaystyle=$
$\displaystyle\left|\langle\widehat{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{\pi(k)}^{(1)}\rangle(1-\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{\pi(k)}^{(1)}\rangle)\right|$ $\displaystyle\leq$
$\displaystyle 1-\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{\pi(k)}^{(1)}\rangle$ $\displaystyle=$
$\displaystyle\|\widehat{{\cal
P}}^{(1)}\mathbf{u}_{\pi(k)}^{(1)}\|^{2}\leq\sum_{1\leq
l<k}\tilde{\delta}_{l}^{2}\leq C\delta_{k}^{2},$

by taking $c_{0}$ small enough. Thus,

$\left|\sigma_{\pi(k)}^{2}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}^{2}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right|\leq C\tau^{2}\delta_{k}^{2}.$ (27)

The second term on the right hand side of (26) can be bounded as follows:

$\displaystyle\left|\sum_{l\neq\pi(k)}\sigma_{\pi(k)}\sigma_{l}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}\theta_{il}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right|$ $\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\left(\sum_{l\neq\pi(k)}\langle\mathbf{u}_{l}^{(1)},\widehat{\mathbf{u}}_{k}^{(1)}\rangle^{2}\right)^{1/2}\max_{l\neq\pi(k)}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}$
$\displaystyle=$
$\displaystyle\tau\delta_{k}\|\widehat{\Sigma}_{\theta}-I\|\max_{l\neq\pi(k)}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}$
$\displaystyle\leq$
$\displaystyle\tau\delta_{k}\|\widehat{\Sigma}_{\theta}-I\|\max\left\\{\max_{l\in\pi([k-1])}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\},\max_{l\notin\pi([k])}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}\right\\}.$

Note that

$\max_{l\in\pi([k-1])}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}\leq\max_{l\in\pi([k-1])}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\widehat{{\cal
P}}_{\perp}^{(q)}\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}\leq\max_{1\leq
l<k}\\{\sigma_{\pi(l)}\delta_{l}^{p-1}\\}\leq 2\tau\delta_{k}^{p-1},$

and

$\max_{l\notin\pi([k])}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle|\right\\}\leq\tau\delta_{k}^{p-1}.$

We get

$\left|\sum_{l\neq\pi(k)}\sigma_{\pi(k)}\sigma_{l}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}\theta_{il}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right|\leq
C\tau^{2}\delta_{k}^{p}\sqrt{r\over n}.$ (28)

And similarly, when
$\mathbf{w}=\widehat{\mathbf{u}}_{k}^{(1)}-\langle\widehat{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{\pi(k)}^{(1)}\rangle\mathbf{u}_{\pi(k)}^{(1)}$,
we bound the third term on the right hand side of (26) by

$\displaystyle\left|\sum_{l\neq\pi(k)}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right|$ (29) $\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\left(\sum_{l\neq\pi(k)}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle^{2}\right)^{1/2}\max_{l\neq\pi(k)}\left\\{\sigma_{l}\prod_{q=2}^{p}|\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right\\}$
$\displaystyle\leq$ $\displaystyle C\tau^{2}\delta_{k}^{p}\sqrt{r\over n},$

where in the last inequality we used the fact that

$\sum_{l\neq\pi(k)}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle^{2}\leq\|\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\|^{2}\leq\|\mathbf{w}\|^{2}\leq\delta_{k}^{2};$

the fourth term by

$\displaystyle\left|\sum_{l_{1},l_{2}\neq\pi(k)}\sigma_{l_{1}}\sigma_{l_{2}}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il_{1}}\theta_{il_{2}}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l_{1}}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right|$ (30) $\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta}\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle C\tau^{2}\delta_{k}^{2(p-1)};$

the fifth term by

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l\neq\pi(k)}\sigma_{l}\theta_{il}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\right]$ (31) $\displaystyle=$
$\displaystyle\sum_{l\neq\pi(k)}\left\\{\sigma_{l}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle\left[{1\over
n}\sum_{i=1}^{n}\theta_{il}\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\right]\right\\}$
$\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\cdot\left[\sum_{l\neq\pi(k)}\left(\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)\right]^{1/2}$
$\displaystyle\leq$ $\displaystyle C\sigma_{0}\tau\delta_{k}^{p-1}\sqrt{d\over
n};$

and the sixth term by

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l=1}^{r}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
(32) $\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{\pi(k)}\theta_{i\pi(k)}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l\neq\pi(k)}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
$\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|+\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle C\sigma_{0}\tau\sqrt{d\over n}.$

Finally the last term can be bounded by

$\left|\frac{1}{n}\sum_{i=1}^{n}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})-\sigma_{0}^{2}\langle\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{w},\widehat{\mathbf{u}}_{k}^{(1)}\rangle\right|\leq\|\widehat{\Sigma}_{\mathscr{E}}-\Sigma_{\mathscr{E}}\|\leq
C\sigma_{0}^{2}\sqrt{d\over n}.$

Together with (27)-(32), we get

$\displaystyle\tilde{\Sigma}(\widehat{{\cal
P}}_{\perp}^{(1)}(\widehat{\mathbf{u}}_{k}^{(1)}-\langle\widehat{\mathbf{u}}_{k}^{(1)},\mathbf{u}_{\pi(k)}^{(1)}\rangle\mathbf{u}_{\pi(k)}^{(1)}),\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})\leq
C\tau^{2}\delta_{k}^{2}+C(\sigma_{0}^{2}+\sigma_{0}\tau)\sqrt{d\over n}.$

###  B.3 Upper Bounds for $\max_{l\notin\pi([k])}\tilde{\Sigma}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)})$

To derive the helper bound (17), we also need an upper bound for

$\max_{l\notin\pi([k])}\tilde{\Sigma}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)},\widehat{\mathbf{u}}_{k}^{(1)},\ldots,\widehat{\mathbf{u}}_{k}^{(p)}).$

We shall follow a similar step by bound each term on the right hand side of
(26), but now with $\mathbf{w}=\mathbf{u}_{m}^{(1)}$ ($m\notin\pi([k])$).

Specifically, the first term can be bounded by

$\displaystyle\left|\sigma_{\pi(k)}^{2}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}^{2}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$
$\displaystyle\leq\tau^{2}\|\widehat{\Sigma}_{\theta}\|_{\max}|\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle|.$

Note that

$\langle\mathbf{u}_{\pi(k)}^{(1)},\mathbf{u}_{m}^{(1)}\rangle=\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}^{(1)}\mathbf{u}_{m}^{(1)}\rangle+\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle=0.$

We get

$|\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle|=|\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}^{(1)}\mathbf{u}_{m}^{(1)}\rangle|\leq\|\widehat{{\cal
P}}\mathbf{u}_{\pi(k)}^{(1)}\|\|\widehat{{\cal
P}}\mathbf{u}_{m}^{(1)}\|\leq\sum_{1\leq l<k}\tilde{\delta}_{l}^{2}\leq
C\delta_{k}^{2},$

by Cauchy-Schwartz inequality. This implies that

$\left|\sigma_{\pi(k)}^{2}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}^{2}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|\leq
C\tau^{2}\delta_{k}^{2}.$

by taking $c_{0}$ small enough.

The second term can also be bounded by

$\displaystyle\left|\sum_{l\neq\pi(k)}\sigma_{\pi(k)}\sigma_{l}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{i\pi(k)}\theta_{il}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$ $\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\||\langle\mathbf{u}_{\pi(k)}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle\tau^{2}\delta_{k}^{p+1}\sqrt{r\over n}.$

We bound the third term by

$\displaystyle\left|\sum_{l\neq\pi(k)}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$ $\displaystyle\leq$
$\displaystyle\left|\sigma_{m}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{im}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{m}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{m}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$
$\displaystyle+\left|\sum_{l\notin\\{\pi(k),m\\}}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$

The first term on the right hand side can be further bounded by

$\tau^{2}\|\widehat{\Sigma}_{\theta}-\Sigma_{\theta}\|_{\max}\prod_{q=2}^{p}\left|\langle\mathbf{u}_{m}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right|\leq
C\tau^{2}\delta_{k}^{p-1}\sqrt{\log r\over n}.$

Now consider the second term:

$\displaystyle\left|\sum_{l\notin\\{\pi(k),m\\}}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$ $\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\left(\sum_{l\notin\\{\pi(k),m\\}}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\tau\|\widehat{\Sigma}_{\theta}-I\|\Biggl{(}\sum_{l\in\pi([k-1])}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}$
$\displaystyle\hskip
100.0pt+\sum_{l\notin\pi([k])\cup\\{m\\}}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\Biggr{)}^{1/2}.$

The first term in the bracket on the rightmost hand side can be bounded by

$\displaystyle\sum_{l\in\pi([k-1])}\sigma_{l}^{2}\|\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{l}^{(1)}\|^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}$
$\displaystyle\leq$
$\displaystyle\left(\sum_{l\in\pi([k-1])}\langle\mathbf{u}_{l}^{(2)},\widehat{\mathbf{u}}_{k}^{(2)}\rangle^{2}\right)\left(\max_{l\in\pi([k-1])}\sigma_{l}^{2}\|\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{l}^{(1)}\|^{2}\prod_{q=3}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)$
$\displaystyle\leq$
$\displaystyle\delta_{k}^{2}\left(\max_{l\in\pi([k-1])}\sigma_{l}^{2}\delta_{l}^{2}\delta_{k}^{2(p-2)}\right)\leq
C\tau^{2}\delta_{k}^{2p};$

the second term by

$\|\widehat{{\cal
P}}^{(1)}\mathbf{u}_{m}^{(1)}\|^{2}\left(\max_{l\notin\pi([k])\cup\\{m\\}}\sigma_{l}^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)\leq
C\tau^{2}\delta_{k}^{2p},$

so that

$\displaystyle\left|\sum_{l\neq\pi(k)}\sigma_{l}\sigma_{\pi(k)}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il}\theta_{i\pi(k)}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$ $\displaystyle\leq$
$\displaystyle C\left(\tau^{2}\delta_{k}^{p-1}\sqrt{\log r\over
n}+\tau^{2}\delta_{k}^{p}\sqrt{r\over n}\right).$

Similar to before, the fourth term on the right hand side of (26) can be
bounded by

$\displaystyle\left|\sum_{l_{1},l_{2}\neq\pi(k)}\sigma_{l_{1}}\sigma_{l_{2}}\left(\frac{1}{n}\sum_{i=1}^{n}\theta_{il_{1}}\theta_{il_{2}}\right)\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l_{1}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l_{2}}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l_{1}}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right|$ $\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta}\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle C\tau^{2}\delta_{k}^{2(p-1)};$

and the fifth term by

$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l\neq\pi(k)}\sigma_{l}\theta_{il}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\right]\right|$ $\displaystyle=$
$\displaystyle\left|\sum_{l\neq\pi(k)}\left\\{\sigma_{l}\left(\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle\left[{1\over
n}\sum_{i=1}^{n}\theta_{il}\mathscr{E}_{i}(\widehat{\mathbf{u}}_{k}^{(1)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\right]\right\\}\right|$
$\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\cdot\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\langle\mathbf{u}_{l}^{(1)},\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)}\rangle^{2}\prod_{q=2}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle C\sigma_{0}\tau\delta_{k}^{p-1}\sqrt{d\over
n}.$

We now turn to the sixth term on the right hand side of (26). Write

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l=1}^{r}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l\neq\pi(k)}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]$
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\left[\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{\pi(k)}\theta_{i\pi(k)}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right].$

The first term again can be bounded by

$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{l\neq\pi(k)}\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{l}\theta_{il}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]\right|$
$\displaystyle\leq$
$\displaystyle\|\widehat{\Sigma}_{\theta,\mathscr{E}}\|\left(\sum_{l\neq\pi(k)}\sigma_{l}^{2}\prod_{q=1}^{p}\langle\mathbf{u}_{l}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle^{2}\right)^{1/2}\leq
C\sigma_{0}\tau\delta_{k}^{p-1}\sqrt{d\over n}.$

For the second term, note that

$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}\left[\mathscr{E}_{i}(\widehat{{\cal
P}}_{\perp}^{(1)}\mathbf{u}_{m}^{(1)},\widehat{\mathbf{u}}_{k}^{(2)},\dots,\widehat{\mathbf{u}}_{k}^{(p)})\sigma_{\pi(k)}\theta_{i\pi(k)}\left(\prod_{q=1}^{p}\langle\mathbf{u}_{\pi(k)}^{(q)},\widehat{\mathbf{u}}_{k}^{(q)}\rangle\right)\right]\right|$