This paper studies a general framework for high-order tensor SVD. We propose a new computationally efficient algorithm, tensor-train orthogonal iteration (TTOI), that aims to estimate the low tensor-train rank structure from the noisy high-order tensor observation. The proposed TTOI consists of initialization via TT-SVD [62] and new iterative backward/forward updates. We develop the general upper bound on estimation error for TTOI with the support of several new representation lemmas on tensor matricizations. By developing a matching information-theoretic lower bound, we also prove that TTOI achieves the minimax optimality under the spiked tensor model. The merits of the proposed TTOI are illustrated through applications to estimation and dimension reduction of high-order Markov processes, numerical studies, and a real data example on New York City taxi travel records. The software of the proposed algorithm is available online (link: <https://github.com/Lili-Zheng-stat/TTOI>).

Address for Correspondence: Anru R. Zhang, Departments of Biostatistics & Bioinformatics, Computer Science, Mathematics, and Statistical Science, Duke University, Durham, NC 27710, USA (email: ).

The research of Yuchen Zhou and Anru R. Zhang were supported in part by NSF under Grants CAREER-1944904, DMS-1811868, and NIH under Grant R01 GM131399; the research of Yazhen Wang was supported in part by NSF under Grants DMS-1707605 and DMS-1913149. This work was done while Yuchen Zhou, Anru R. Zhang, and Lili Zheng were at the University of Wisconsin-Madison.

Yuchen Zhou is with the Department of Statistics and Data Science, The Wharton School, University of Pennsylvania, Philadelphia, PA 19104, USA (email: ).

Lili Zheng is with Department of Electrical and Computer Engineering, Rice University, Houston, TX 77005, USA (email: ).

Yazhen Wang is with the Department of Statistics, University of Wisconsin-Madison, Madison, WI 53706, USA (email: ).

§ INTRODUCTION

Tensors, or high-order arrays, have attracted increasing attention in modern machine learning, computational mathematics, statistics, and data science. Some specific examples include recommender systems [10, 54], neuroimaging analysis [85, 94], latent variable learning [4], multidimensional convolution [63], signal processing [20], neural network [52, 93], computational imaging [42, 89], contingency table [9, 29]. In addition to low-order tensors (e.g., tensor with a relatively small value of order number), the high-order tensors also commonly arise in applications in statistics and machine learning. For example, in convolutional neural networks, parameters in fully connected layers can be represented as high-order tensors [16, 56]. In an order-$d$ Markov process, where the future states depend on jointly the current and $(d-1)$ previous states, the transition probabilities form an order-$(d+1)$ tensor. For an order-$d$ Markov decision process, the transition probabilities can be represented by an order-$(2d+1)$ tensor, with additional $d$ directions representing past $d$ actions. High-order tensors are also used to represent the joint probability in Markov random fields [57].

Compared to the low-order tensors, high-order tensors encompass much more parameters and sophisticated structure, while leading to inhibitive cost in storage, processing, and analysis: an order-$d$ dimension-$p$ tensor contains $p^d$ parameters. To address this issue, some low-dimensional parametrization is usually considered to capture the most informative subspaces in the tensor.
In particular, the tensor-train (TT) decomposition [30, 59, 60, 62, 58] introduced a classic low-dimensional parameterization to model the subspaces and latent cores in high-order tensor structures. TT decomposition has been used in a wide range of applications in physics and quantum computation [12, 30, 58, 68, 71], signal processing [20], and supervised learning [76] among many others. For example, the TT decomposition framework is utilized in quantum information science for modeling complex quantum states and handling the quantum mean value problem [12, 30, 58, 68]. The TT-decomposition of a tensor $\bfc{X} \in \bbR^{p_1 \times \cdots \times p_d}$ is defined as below:

\begin{equation}\label{TT}
\begin{split}
\bfc{X}_{i_1, \cdots, i_d} = & G_{1, [i_1, :]}\bfc{G}_{2, [:, i_2, :]}\cdots \bfc{G}_{d-1, [:, i_{d-1}, :]}G_{d, [i_d, :]}^\top\\
= & \sum_{\alpha_1 = 1}^{r_1}\cdots \sum_{\alpha_{d-1} = 1}^{r_{d-1}}G_{1, [i_1, \alpha_1]}\bfc{G}_{2, [\alpha_1, i_2, \alpha_2]}\cdots \bfc{G}_{d-1, [\alpha_{d-2}, i_{d-1}, \alpha_{d-1}]}G_{d, [i_d, \alpha_{d-1}]}.
\end{split}
\end{equation}

Here, the smallest values of $r_1, \dots, r_{d-1}$ that enable the decomposition (<ref>) are called the TT-rank of $\bfc{X}$. [62] shows that the TT-rank $r_k = \rank([\bfc{X}]_k)$, i.e., the rank of the $k$th sequential unfolding of $\bfc{X}$ (see formal definition of sequential unfolding in Section <ref>). $G_1 \in \bbR^{p_1 \times r_1}$, $\bfc{G}_k \in \bbR^{r_{k-1} \times p_k \times r_k},  G_d \in \bbR^{p_d \times r_{d-1}}$ are the TT-cores that multiply sequentially like a “train": $\bfc{X}_{i_1, \cdots, i_d}$ equals the product of $i_1$th vector in $G_1$, $i_2$th matrix in $\bfc{G}_{2}$, …, $i_{d-1}$th matrix in $\bfc{G}_{d-1}$, and $i_d$th vector in $G_d$.  For convenience of presentation, we simplify (<ref>) to

$$\bfc{X} = \llbracket G_1, \bfc{G}_2, \ldots, \bfc{G}_{d-1}, G_d\rrbracket$$

and denote $r_0 = r_d = 1$ throughout the paper. In particular, the TT rank and TT decomposition reduce to the regular matrix rank and decomposition when $d=2$. If all dimensions $p$ and ranks $r$ are the same, the TT-parametrization involves $O(2pr + (d-2)pr^2)$ values, which can be significantly smaller than the ones for Tucker-decomposition $O(r^d + dpr)$ and the regular parameterization $O(p^d)$.

In most of the existing literature, the TT-decomposition was considered under the deterministic settings, and the central goal was often to approximate the nonrandom high-order tensors by low-dimensional structures [11, 61, 62]. However, in modern applications in data science such as Markov processes, Markov decision processes, and Markov random fields, the (transition) probability tensor computed based on data is often a random realization of the underlying true tensor. In these cases, the estimation of the underlying low-dimensional parameters hidden in the noisy observations can be more important: an accurate estimation of the transition tensor renders reliable prediction for future states in high-order Markov chains and better decision-making in high-order Markov decision processes; an accurate estimation of probability tensor sheds light on the underlying relationship among different variables in a random system [57]. To achieve such a goal, it is crucial to develop dimension reduction methods that can incorporate TT-decomposition into probabilistic models. Since singular value decomposition (SVD) is one of the most important dimension reduction methods involving probabilistic models for matrices, and there is no counterpart of it for high-order tensors, we aim to fill this void by developing a statistical framework and a computationally feasible method for high-order tensor SVD in this paper.

§.§ Problem Formulation

This paper focuses on the following high-order tensor SVD model. Suppose we observe an order-$d$ tensor $\bfc{Y}$ that contains a hidden tensor-train (TT) low-rank structure:

\begin{equation}\label{eq:obs_model}
\bfc{Y} = \bfc{X} + \bfc{Z}, \quad \bfc{Y}, \bfc{X}, \bfc{Z} \in \bbR^{\otimes_{k=1}^d p_k}.
\end{equation}

Here, $\bfc{X}$ is TT-decomposable as (<ref>) and $\bfc{Z}$ is a noise tensor. Our goal is to estimate $\bfc{X}$ and the TT cores of $\bfc{X}$ based on $\bfc{Y}$. To this end, a straightforward idea is to minimize the approximation error as follows,

\begin{equation}\label{eq:minimization}
\widehat{\bfc{X}} = \argmin_{\bfc{A} \text{ is decomposable as } \eqref{TT}}\left\|\bfc{Y} - \bfc{A}\right\|_{\F}^2.
\end{equation}

However, the approximation error minimization (<ref>) is highly non-convex and finding the global optimal solution, even if the rank $r_1 = \cdots = r_{d-1}= 1$, is NP-hard in general [34]. Instead, a variety of computationally feasible methods have been proposed to approximate the best tensor-train low-rank decomposition in the literature. TT-SVD, a sequential singular value thresholding scheme, was introduced by [62] to be discussed in detail later. [62] also proposed TT-rounding via sequential QR decompositions, which reduces the TT-rank while ensuring approximation accuracy. [24] introduced the alternating minimal energy algorithm to reconstruct a TT-low-rank tensor approximately based on only a small proportion of revealed entries of the target tensor. <cit.> proposed a sketching-based algorithm for fast low TT rank approximation of arbitrary tensors. [11] studied the tensor-train decomposition for functional tensors. [41] proposed the FastTT algorithm for fast sparse tensor decomposition based on parallel vector rounding and TT-rounding. [48] studied dynamical approximation with TT format for time-dependent tensors. [32] proposed the alternating least squares for tensor completion in the TT format. [6] studied the completion of low TT rank tensor and the applications to color image and video recovery. [75] studied the Riemannian optimization methods for TT decomposition and completion.
Also see [55] for a TT decomposition library in TensorFlow. To our best knowledge, the estimation performance of most procedures here remains unclear. Departing from these existing work, in this paper, we make a first attempt to minimize the estimation error of $\bfc{X}$ in addition to achieving the minimal approximation error under possibly random settings.

§.§ Our Contributions

Under Model (<ref>), we make the following contributions to high-order tensor SVD in this paper.

First, we propose a new algorithm, Tensor-Train Orthogonal Iteration (TTOI), that provides a computationally efficient estimation of the low-rank TT structure from the noisy observation. The proposed algorithm includes two major steps. First, we obtain initial estimates $\widehat{G}_1^{(0)}, \widehat{\bfc{G}}_2^{(0)},\dots, \widehat{\bfc{G}}_{d-1}^{(0)}, \widehat{G}_d$ by performing forward sequential SVD based on matricizations and projections. This step was known as TT-SVD in the literature [62]. Next, we utilize the initialization and perform the newly developed backward updates and forward updates alternatively and iteratively. The TTOI procedure will be discussed in detail in Section <ref>.

To see why the TTOI iterations yield better estimation than the classic TT-SVD method, recall that TT-SVD first performs singular value thresholding on $[\bfc{Y}]_1$, i.e., the unfolding of $\bfc{Y}$, without any additional updates (see detailed procedure of TT-SVD and formal definition of $[\bfc{Y}]_1$ in Section <ref>), which can be inaccurate since $[\bfc{Y}]_1$, a $p_1$-by-$\prod_{k=2}^dp_k$ matrix, has a great number of columns. In contrast, TTOI iteration utilizes the intermediate outcome of the previous iteration to substantially reduce the dimension of $[\bfc{Y}]_1$ while performing singular value thresholding. In Figure <ref>, we provide a simple simulation example to show that even one TTOI iteration can significantly improve the estimation of the left singular subspace of $G_1$ (left panel) and the overall tensor $\bfc{X}$ (right panel). Therefore, a one-step TTOI, i.e., the initialization with one TTOI iteration, can be used in practice when the computational cost is a concern.

Average estimation error (dots) and standard deviation (bars) of $\|\sin\Theta(\widehat{U}_1, U_1)\|$ and $\|\widehat{\bfc{X}} - \bfc{X}\|_\F$ by  TT-SVD and one-step TTOI.  Both algorithms are performed based on the observation $\bfc{Y}$ generated from (<ref>), where $\bfc{Z}\overset{\text{i.i.d.}}{\sim}N(0,\sigma^2)$, $\bfc{X}$ is a randomly generated order-5 tensor based on (<ref>) with $p=20, r=1$, $G_1,\bfc{G}_2,\dots,\bfc{G}_{d-1}, G_d \overset{\text{i.i.d.}}{\sim} N(0,1)$.

We develop theoretical guarantees for TTOI. In particular, we introduce a series of representation lemmas for tensor matricizations with TT format. Based on them, we develop a deterministic upper bound of estimation error for both forward and backward updates in TTOI iterations. Under the benchmark setting of spiked tensor model, we develop matching upper/lower bounds and prove that the proposed TTOI algorithm achieves the minimax optimal rate of estimation error. To the best of our knowledge, this is the first statistical optimality results for high-order tensors with TT format. We also prove for any high-order tensor, TTOI iteration has monotone decreasing approximation error with respect to the iteration index.

Moreover, to break the curse of dimensionality in high-order Markov processes, we study the state aggregatable high-order Markov processes and establish a key connection to TT decomposable tensors. We propose a TTOI estimator for the transition probability tensor in high-order state-aggregatable Markov processes and establish the theoretical guarantee. We conduct simulation experiments to demonstrate the performance of TTOI and validate our theoretical findings. We also apply our method to analyze a New York taxi dataset. By modeling taxi trips as trajectories realized from a citywide Markov chain, we found that the Manhattan traffic zone exhibits high-order Markovian dependence and the proposed TTOI reveals latent traffic patterns and meaningful partition of Manhattan traffic zones.  Finally, we discuss several applications that our proposed algorithm is applicable to, including transition probability tensor estimation in high-order Markov decision processes and joint probability tensor estimation in Markov random fields.

§.§ Related Literature

In addition to the aforementioned literature on TT decomposition, our work is also related to a substantial body of work on matrix/tensor decomposition and SVD, spiked tensor model, etc. These literature are from a range of communities including  applied mathematics, information theory, machine learning, scientific computing, signal processing, and statistics. Here we try to review existing literature in these communities without claiming this literature survey is exhaustive.

First, the matrix singular value thresholding was commonly used and extensively studied in various problems in data science, including matrix denoising [15, 17, 25], matrix completion [13, 19, 37, 90], principal component analysis (PCA) [53], Markov chain state aggregation [87]. Such the task was also widely considered for tensors of order-3 or higher. In particular, to perform SVD and decomposition for tensors with Tucker low-rank structures, [22, 23] introduced the higher-order SVD (HOSVD) and higher-order orthogonal iteration (HOOI). [88] established the statistical and computational limits of tensor SVD, compared the theoretical properties of HOSVD and HOOI, and proved that HOOI achieves both statistical and computational optimality. [79] introduced the sequentially truncated higher-order singular value decomposition (ST-HOSVD). [86] introduced a thresholding & projection based algorithm for sparse tensor SVD. A non-exhaustive list of methods for SVD and decomposition for tensors with CP low-rank structures include alternating least squares [38, 72], eigendecomposition-based approach [40], enhanced line search [67], power iteration with SVD-based initialization [4], simultaneous diagonalization and higher-order SVD [21].

In addition, the spiked tensor model and tensor principal component analysis (tensor PCA) are widely discussed in the literature. [3, 5, 35, 51, 64, 69] considered the statistical and computational limits of rank-$1$ spiked tensor model. [39] studied the statistical and computational phase transitions and theoretical properties of the approximate message passing algorithm (AMP) under a Bayesian spiked tensor model. [1, 2] developed the regularization-based methods for tensor PCA. [44, 46, 47, 95] studied the robust tensor PCA to handle the possible outliers from the tensor observation.

Different from Tucker and CP decompositions, which have been a pinpoint in the enormous existing literature on tensors, we focus on the TT-structure associated with high-order tensors for the following reasons: (1) Tucker and CP decompositions do not involve the sequential structure of different modes, i.e., the Tucker and CP decompositions still hold if the $d$ modes are arbitrarily permuted. While in applications such as high-order Markov process, high-order Markov decision process, and fully connected layers of deep neural networks, the order of different modes can be crucial; (2) the number of entries involved in the low-Tucker-rank parameterization grows exponentially with respect to the order $d$ ($r^d$); (3) methods that explore CP low-rank structure can be numerically unstable for high-order tensors in computation as pointed out by [61]. In comparison, the TT-structure incorporates the order of different modes sequentially and involves much fewer parameters for high-order tensors, which renders it more suitable in many scenarios.

In Section <ref>, we will further discuss the application of TTOI on high-order Markov processes and state aggregation. This problem is related to a body of literature on dimension reduction and state aggregation for Markov processes that we will discuss in Section <ref>.

§.§ Organization

The rest of the article is organized as follows. In Section <ref>, after a brief  introduction of the notation and preliminaries, we introduce the procedure of the tensor-train orthogonal iteration. The theoretical results, including three representation lemmas, a general estimation error bound, and the minimax optimal upper and lower bounds under the spiked tensor model, are provided in Sections <ref> and <ref>. The application to high-order Markov chains is discussed in Section <ref>. The simulation and real data analysis are provided in Sections <ref> and <ref>, respectively. Discussions and further applications to Markov random fields and high-order Markov decision processes are briefly discussed in Section <ref>. All technical proofs are provided in Section <ref>.

§ PROCEDURE OF TENSOR-TRAIN ORTHOGONAL ITERATION

§.§ Notation and Preliminaries

We first introduce the notation and preliminaries to be used throughout the paper. We use the lowercase letters, e.g., $x, y, z,$ to denote scalars or vectors. We use $C, c, C_0, c_0, \dots$ to denote generic constants, whose actual values may change from line to line. A random variable $z$ is $\sigma$-sub-Gaussian if $\bbE e^{t(z - \bbE z)} \leq e^{\sigma^2 t^2/2}$ for any  $t \in \bbR.$ We say $a \lesssim b$ or $a = O(b)$ if $a \leq Cb$ for some uniform constant $C > 0$. We write $a = \widetilde{O}(b)$ if $a = O(b\log^{C'}(b))$ for constant $C'>0$. The capital letters, e.g., $X, Y, Z$, are used to denote matrices. Specifically, $\mathbb{O}_{p, r} := \{U \in \bbR^{p \times r}: U^\top U = I_r\}$ is the set of all $p$-by-$r$ matrices with orthogonal columns. For $U \in \mathbb{O}_{p, r}$, let $U_{\perp} \in \mathbb{O}_{p, p-r}$ be the orthonormal complement of $U$, and let $P_{U} = UU^\top$ denote the projection matrix onto the column space of $U$. For any matrix $A \in \bbR^{p_1 \times p_2}$, let $A = \sum_{i=1}^{p_1\wedge p_2}s_i u_iv_i^\top$ be the singular value decomposition, where $s_{1}(A) \geq \cdots \geq s_{p_1 \wedge p_2}(A) \geq 0$ are the singular values of $A$ in non-increasing order. Define $s_{\min}(A) = s_{p_1 \wedge p_2}(A)$, $\text{SVD}_{r}^L(A) = [u_1 \ldots u_r] \in \mathbb{O}_{p_1, r}$, and $\text{SVD}_{r}^R(A) = [v_1 \ldots v_r] \in \mathbb{O}_{p_2, r}$ be the smallest non-trivial singular value, leading $r$ left singular vectors, and leading $r$ right singular vectors of $A$, respectively. We also write $\text{SVD}^L(A) = \text{SVD}_{p_1 \wedge p_2}^L(A)$ and $\text{SVD}^R(A)=\text{SVD}_{p_1 \wedge p_2}^L(A)$ as the collection of all left and right singular vectors of $A$, respectively. Define the Frobenius and spectral norms of $A$ as $\|A\|_{\F} = \sqrt{\sum_{i=1}^{p_1}\sum_{j=1}^{p_2}A_{ij}^2} = \sqrt{\sum_{i=1}^{p_1 \wedge p_2}s_i^2(A)}$ and $\|A\| = s_{1}(A) = \max_{x \in \bbR^{p_2}}\|Ax\|_2/\|x\|_2$.
For any two matrices $U \in \bbR^{m_1 \times n_1}$ and $V \in \bbR^{m_2 \times n_2}$, let

\begin{equation*}
U \otimes V =
\begin{bmatrix}
U_{11}\cdot V & \dots & U_{1n_1}\cdot V\\
\vdots &  & \vdots\\
U_{m_1 1}\cdot V & \dots & U_{m_1 n_1}\cdot V
\end{bmatrix} \in\mathbb{R}^{(m_1m_2)\times(n_1n_2)}
\end{equation*}

be their Kronecker product. To quantify the distance among subspaces, we define the principle angles between $U, \widehat{U} \in \mathbb{O}_{p ,r}$ as an $r$-by-$r$ diagonal matrix: $\Theta(U, \widehat{U}) = \diag(\text{arccos}(s_1), \dots, \text{arccos}(s_r))$, where $s_1 \geq \cdots \geq s_r \geq 0$ are the singular values of $U^\top \widehat{U}$. Define the sin$\Theta$ norm as

\begin{equation*}
\|\sin\Theta(U, \widehat{U})\| = \|\diag\left(\sin(\text{arccos}(s_1)), \dots, \sin(\text{arccos}(s_r))\right)\| = \sqrt{1 - s_r^2}.
\end{equation*}

The boldface calligraphic letters, e.g., $\bfc{X}, \bfc{Y}, \bfc{Z}$, are used to denote tensors. For an order-$d$ tensor $\bfc{X} \in \mathbb{R}^{\otimes_{i=1}^d p_i}$ and $1 \leq k \leq d-1$, we define $[\bfc{X}]_k \in \bbR^{(p_1\times \cdots \times p_k) \times (p_{k+1}\cdots p_d)}$ as the sequential unfolding of $\bfc{X}$ with rows enumerating all indices in Modes $1, \dots, k$ and columns enumerating all indices in Modes $(k + 1), \cdots, d$, respectively. That is, for any $1 \leq k \leq d$ and $1\leq i_k\leq p_k$,

$$([\bfc{X}]_k)_{(i_k-1)p_1\cdots p_{k-1} + (i_{k-1}-1)p_1\cdots p_{k-2} + \cdots + i_1, (i_d - 1)p_{k+1}\cdots p_{d-1} + (i_{d-1}-1)p_{k+1}\cdots p_{d-2} + \cdots + i_{k+1}} = \bfc{X}_{i_1 \dots i_d}.$$

Following the convention of  function in MATLAB, we define the reshape of any matrix $X$ of dimension $p_1\cdots p_k\times p_{k+1}\cdots p_d$ as an inverse operation of tensor matricization: $\bfc{X}=\mathrm{Reshape}(X,p_1,p_2,\dots,p_d)$ if $X=[\bfc{X}]_k$. For any two matrices $A\in \mathbb{R}^{q_1\times q_2q_3}$ and $\widetilde{A}\in \mathbb{R}^{q_1q_2\times q_3}$, we denote $\widetilde{A} = \mathrm{Reshape}(A, q_1q_2, q_3)$ and $A = \mathrm{Reshape}(\widetilde{A}, q_1, q_2q_3)$ if and only if
$$\widetilde{A}_{(i_2-1)p_1 + i_1, i_3} = A_{i_1, (i_3-1)p_2 + i_2}, \quad \forall 1 \leq i_j \leq q_j, \quad j =1,2,3.$$

We also define the tensor Frobenius norm of $\bfc{X}$ as $\|\bfc{X}\|_{\F}^2 = \sum_{i_1 = 1}^{p_1}\cdots\sum_{i_d = 1}^{p_d} \bfc{X}_{i_1, \dots, i_d}^2.$ For any matrix $A \in \bbR^{p_1 \times p_2}$ and any tensor $\bfc{B} \in \bbR^{p_1 \times \cdots \times p_d}$, let vec$(A)$ and vec$(\bfc{B})$ be the vectorization of $A$ and $\bfc{B}$, respectively. Formally,

$$\left(\text{vec}(\bfc{B})\right)_{(i_d-1)p_1\cdots p_{d-1} + (i_{d-1}-1)p_1\cdots p_{d-2} + \cdots + i_1} = \bfc{B}_{i_1, \dots, i_d}, \quad 1\leq i_k\leq p_k, \quad k=1,\ldots, d.$$

§.§ Procedure of Tensor-Train Orthogonal Iteration

We are now in position to introduce the procedure of Tensor-Train Orthogonal Iteration (TTOI). The pseudocode of the overall procedure is given in Algorithm <ref>. TTOI includes three main parts: we first run initialization, then perform backward update and forward update alternatively and iteratively.
Tensor-Train Orthogonal Iteration (TTOI)
Input: $\bfc{Y}, \{p_k\}_{k=1}^d, \{r_k\}_{k=1}^{d-1}$, increment tolerance $\varepsilon > 0$, maximum number of iterations $t_{\max}$
[1]
Obtain Initialization $\widetilde{R}_1^{(0)}, \dots, \widetilde{R}_{d-1}^{(0)}, \widehat{\bfc{X}}^{(0)}$ by Algorithm <ref>
$t = 1, \ldots, t_{\max}$
$t$ is odd
Apply Algorithm <ref> with input $\widetilde{R}_1^{(t-1)}, \dots, \widetilde{R}_{d-1}^{(t-1)}$ to obtain $\widehat{V}^{(t)}_1, \dots, \widehat{V}^{(t)}_d, \widehat{\bfc{X}}^{(t)}$
Apply Algorithm <ref> with input $\widehat{V}^{(t-1)}_1, \dots, \widehat{V}^{(t-1)}_d$ to obtain $\widetilde{R}_1^{(t)}, \dots, \widetilde{R}_{d-1}^{(t)}, \widehat{\bfc{X}}^{(t)}$
If  $\|\widehat{\bfc{X}}^{(t)}\|_{\F}^2 - \|\widehat{\bfc{X}}^{(t-1)}\|_{\F}^2 \leq \varepsilon$ then    break from the for loop

Output: $\widehat{\bfc{X}}=\widehat{\bfc{X}}^{(t)}$

* Part 1: Initialization. First, we obtain an initial estimate of TT-cores $G_1, \bfc{G}_2, \dots, \bfc{G}_{d-1}, G_d$. This step is the tensor-train-singular value decomposition (TT-SVD) originally introduced by [62].

(i) Let $R_1^{(0)}$ be the unfolding of $\bfc{Y}$ along Mode 1. We compute the top-$r_1$ SVD of $R_1^{(0)}$. Let $\widehat{U}_1^{(0)} \in \mathbb{O}_{p_1, r_1}$ be the first $r_1$ left singular vectors of $R_1^{(0)}$ and calculate $\widetilde{R}_1^{(0)} = (\widehat{U}_1^{(0)})^\top R_1^{(0)} \in \bbR^{r_1 \times (p_2\dots p_d)}$. Then, $\widehat{U}_1^{(0)}$ is an initial estimate of the subspace that $G_1$ lies in and $\widetilde{R}_1^{(0)}$ can be seen as the projection residual.

(ii) Next, we realign the entries of $\widetilde{R}_1^{(0)} \in \bbR^{r_1 \times (p_2\dots p_d)}$ to $R_2^{(0)} \in \bbR^{(r_1p_2) \times (p_3\dots p_d)}$, where the rows and columns of $R_2^{(0)}$ correspond to indices of Modes-1, 2 and Modes-$3,\ldots, d$, respectively. Then, we evaluate the top-$r_2$ SVD of $R_2^{(0)}$. Let $\widehat{U}_2^{(0)}$ be the first $r_2$ left singular vectors of $R_2^{(0)}$ and evaluate $\widetilde{R}_2^{(0)} = (\widehat{U}_2^{(0)})^\top R_2^{(0)} \in \bbR^{r_2 \times p_3\dots p_d}$. Again, $\widehat{U}_2^{(0)}$ is an estimate of the singular subspace that $\mathcal{G}_2$ lies on and $\widetilde{R}_2^{(0)}$ is the projection residual for the next calculation.

(iii) We apply Step (ii) on $\widetilde{R}_2^{(0)}$ to obtain $\widehat{U}_3^{(0)}\in \mathbb{O}_{r_2p_3, r_3}$ and $\widetilde{R}_3^{(0)}\in\mathbb{R}^{r_3\times (p_4\cdots p_d)}$; $\ldots$; apply Step (ii) on $\widetilde{R}_{d-2}^{(0)}$ to obtain $\widehat{U}_{d-1}^{(0)}\in \mathbb{O}_{r_{d-2}p_{d-1}, r_{d-1}}$ and $\widetilde{R}_{d-1}^{(0)}\in \mathbb{R}^{r_{d-1}\times p_d}$. Then we reshape matrix $\widehat{U}_k^{(0)} \in \bbR^{(p_kr_{k-1}) \times r_k}$ to tensor $\widehat{\bfc{U}}_k^{(0)} \in \bbR^{r_{k-1} \times p_k \times r_k}$ for $k = 2, \dots, d - 1$. Now,
$\left(\widehat{U}_1^{(0)}, \widehat{\bfc{U}}_2^{(0)}, \dots, \widehat{\bfc{U}}_{d-1}^{(0)}, \widetilde{R}_{d-1}^{(0)\top}\right)$
yield the initial estimates of TT-cores of $\bfc{X}$ and we expect that

$$\bfc{X}\approx \bfc{X}^{(0)} = \llbracket\widehat{U}_{1}^{(0)}, \widehat{\bfc{U}}_{2}^{(0)}, \cdots,  \widehat{\bfc{U}}_{d-1}^{(0)}, \widetilde{R}_{d-1}^{(0)}\rrbracket.$$

The initialization step is summarized to Algorithm <ref> and illustrated in Figure <ref>. In summary, we perform SVD on some “residual" $R_k^{(0)}$ sequentially for $k=1,\ldots, d-1$. As will be shown in Lemma <ref>, $R_k^{(0)}$ satisfies

$$R_k^{(0)} = (I_{p_k} \otimes \widehat{U}^{(0)\top}_{k - 1})\cdots (I_{p_2\cdots p_k} \otimes \widehat{U}^{(0)\top}_{1})[\bfc{Y}]_k,$$

where $[\bfc{Y}]_k \in \bbR^{(p_1\cdots p_k) \times (p_{k+1}\cdots p_d)}$ is the $k$th sequential unfolding of $\bfc{Y}$ (see definition in Section <ref>). This quantity plays a key role in the backward update next.
Initialization (TT-SVD [62])
Input: $\bfc{Y}$, $\{r_k\}_{k=1}^{d-1}$, $\{p_k\}_{k=1}^d$
[1]
Calculate $R_1^{(0)} = [\bfc{Y}]_1$
$\widehat{U}_k^{(0)} = \text{SVD}_{r_k}^L(R_k^{(0)})$
If $k=1$   then   $U_{\text{prod}}^{(0)} = \widehat{U}_k^{(0)}$   else   $U_{\text{prod}}^{(0)} = (I_{p_k} \otimes U_{\text{prod}}^{(0)})\widehat{U}_k^{(0)}$
$\widetilde{R}_k^{(0)} = \widehat{U}_k^{(0)\top} R_k^{(0)}$
If $k < d - 1$   then  $R_{k+1}^{(0)} = \mathrm{reshape}(\widetilde{R}_k^{(0)}, r_kp_{k+1}, p_{k+2}\cdots p_d)$
$[\widehat{X}^{(0)}]_{d-1} = U_{\text{prod}}^{(0)}\widetilde{R}_{d-1}^{(0)}$
Reshape $[\widehat{X}^{(0)}]_{d-1} \in \bbR^{(p_1\cdots p_{d-1}) \times p_d}$ to $\widehat{\bfc{X}}^{(0)} \in \bbR^{p_1 \times \cdots \times p_d}$

Output: $\widetilde{R}_1^{(0)}, \dots, \widetilde{R}_{d-1}^{(0)}, \widehat{\bfc{X}}^{(0)}$

A Pictorial Illustration of Initialization (Algorithm <ref>, $d=3$)

TT-Backward Update
Input: $\bfc{Y}, \{r_k\}_{k=1}^{d-1}, \{p_k\}_{k=1}^d, \widetilde{R}_1^{(t-1)}, \dots, \widetilde{R}_{d-1}^{(t-1)}$ for odd iteration number $t$
[1]
$\widehat{V}_{d-k+1}^{(t)} = \text{SVD}_{r_{d-k}}^{R}\left(\widetilde{R}_{d-k}^{(t-1)}\right)$,   $V_{\text{prod}}^{(t)} = \widehat{V}_{d-k+1}^{(t)}$
$\widehat{V}_{d-k+1}^{(t)} = \text{SVD}_{r_{d-k}}^{R}\left( \widetilde{R}_{d-k}^{(t-1)}(V_{\text{prod}}^{(t)} \otimes I_{p_{d-k+1}})\right)$,  $V_{\text{prod}}^{(t)} = (V_{\text{prod}}^{(t)} \otimes I_{p_{d-k+1}})\widehat{V}_{d-k+1}^{(t)}$
$\widehat{V}_{1}^{(t)} = [\bfc{Y}]_1V_{\text{prod}}^{(t)}$,   $[\widehat{X}^{(t)}]_1 = \widehat{V}_{1}^{(t)}V_{\text{prod}}^{{(t)}\top}$,   reshape $[\widehat{X}^{(t)}]_1 \in \bbR^{p_1 \times (p_2\cdots p_d)}$ to $\widehat{\bfc{X}}^{(t)} \in \bbR^{p_1 \times \cdots \times p_d}$

Output: $\widehat{V}_1^{(t)}, \dots, \widehat{V}_d^{(t)}, \widehat{\bfc{X}}^{(t)}$

TT-Forward Update
Input: $\bfc{Y}, \{r_k\}_{k=1}^{d-1}, \{p_k\}_{k=1}^d, \widehat{V}_1^{(t-1)}, \dots, \widehat{V}_d^{(t-1)}$ for even iteration number $t$
[1]
$R_1^{(t)} = [\bfc{Y}]_1$
$\widehat{U}_{1}^{(t)} = \text{SVD}_{r_1}^{L}\left(\widehat{V}_{1}^{(t)}\right)$,   $U_{\text{prod}}^{(t)} = \widehat{U}_{1}^{(t)}$
$\widehat{U}_{k}^{(t)} = \text{SVD}_{r_{k}}^{R}\left(R_k^{(t)}(\widehat{V}_d^{(t-1)} \otimes I_{p_{k+1}\dots p_{d-1}})\cdots (\widehat{V}_{k+2}^{(t-1)} \otimes I_{p_{k+1}})\widehat{V}_{k+1}^{(t-1)}\right)$
$U_{\text{prod}}^{(t)} = (I_{p_k} \otimes U_{\text{prod}}^{(t)})\widehat{U}_{k}^{(t)}$
$\widetilde{R}_k^{(t)} = \widehat{U}_k^{(t)\top} R_k^{(t)}$
If $k < d - 1$   then   $R_{k+1}^{(t)} = \text{reshape}\left(\widetilde{R}_k^{(t)}, r_kp_{k+1}, p_{k+2}\cdots p_d\right)$
$[\widehat{X}^{(t)}]_{d-1} = U_{\text{prod}}^{(t)}\widetilde{R}_{d-1}^{(0)}$,  reshape $[\widehat{X}^{(t)}]_{d-1} \in \bbR^{(p_1\cdots p_{d-1}) \times p_d}$ to $\widehat{\bfc{X}}^{(t)} \in \bbR^{p_1 \times \cdots \times p_d}$

Output: $\widetilde{R}_1^{(t)}, \dots, \widetilde{R}_{d-1}^{(t)}, \widehat{\bfc{X}}^{(t)}$

A pictorial illustration of TT-Backward update (Algorithm <ref>,  $d=3$)

The initialization step mainly focuses on the left singular spaces of $[\bfc{X}]_k$ while ignoring the information included in the right singular spaces. Due to this fact, we develop the following new backward update that utilizes both the left and right singular space estimates from the previous step to refine our estimates. Similarly, we can also perform a forward update to further improve the outcome of backward update, and then iteratively alternate between backward and forward updates. The detailed descriptions of these two updates are presented as follows, and a further explanation is given in Remark <ref>.

* Part 2: Backward update. For iterations $t=1,3,5, \ldots$, we perform backward update, i.e., to sequentially obtain $\widehat{V}^{(t)}_d, \ldots, \widehat{V}^{(t)}_2$ based on the intermediate results from the $(t-1)$st iteration (0th iteration is the initialization). The pseudocode of backward update is provided in Algorithm <ref>. The calculation in Algorithm <ref> is equivalent to

$$\widehat{V}_d^{(t)} = \SVD^R\left(\widetilde{R}_{d-1}^{(t-1)}\right),$$

$$\widehat{V}_k^{(t)} = \SVD^R\left(\widetilde{R}_{k-1}^{(t-1)}(\widehat{V}_d^{(t)} \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k + 1}^{(t)} \otimes I_{p_k})\right), \quad k=d-1,\ldots, 2,$$

$$\widehat{V}_1^{(t)} = [\bfc{Y}]_1(\widehat{V}_d^{(t)} \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3}^{(t)} \otimes I_{p_2})\widehat{V}_{2}^{(t)} \in \bbR^{p_1 \times r_1}.$$

$$\widetilde{R}_k^{(t-1)} = (\widehat{U}_k^{(t-1)})^\top(I_{p_k} \otimes \widehat{U}^{(t-1)\top}_{k - 1})\cdots (I_{p_2\cdots p_k} \otimes \widehat{U}^{(t-1)\top}_{1})[\bfc{Y}]_k$$

are the projection residual term in the intermediate outcome of the $(t-1)$st iteration. Then, we reshape $\widehat{V}_k^{(t)\top} \in \bbR^{r_{k-1} \times (p_kr_k)}$ to  $\widehat{\bfc{V}}_k^{(t)} \in \bbR^{r_{k-1} \times p_k \times r_k}$. The backward updated estimate is

$$\widehat{\bfc{X}}^{(t)} =  \llbracket\widehat{V}_1^{(t)}, \widehat{\bfc{V}}_2^{(t)}, \dots, \widehat{\bfc{V}}_{d-1}^{(t)}, \widehat{V}_{d}^{(t)}\rrbracket.$$

[Interpretation of backward update]
The backward updates utilize and extract the right singular vectors of the intermediate products of the $(t-1)$st iteration,

$$\widetilde{R}_k^{(t-1)} = (\widehat{U}_k^{(t-1)})^\top(I_{p_k} \otimes \widehat{U}^{(t-1)\top}_{k - 1})\cdots (I_{p_2\cdots p_k} \otimes \widehat{U}^{(t-1)\top}_{1})[\bfc{Y}]_k,$$

as opposed to the entire data $[\bfc{Y}]_k$. Such a dimension reduction scheme is the key to the backward update: it can simultaneously reduce the dimension of the matrix of interest, $[\bfc{Y}]_k$, and the noise therein, while preserving the signal strength.
Different from the initialization in Step 1, the backward update utilizes the information from both the forward and backward singular subspaces of the tensor-train structure of $\bfc{X}$. See Section <ref> for more illustration.

* Part 3: Forward Update. For iteration $t=2,4,6,\ldots$, we perform forward update, i.e., to sequentially obtain $\widehat{U}_1^{(t)}, \ldots, \widehat{U}_d^{(t)}$ based on the intermediate results from the $(t-1)$st iteration. Essentially, the forward update can be seen as a reversion of the backward update by flipping all modes of tensor $\bfc{Y}$. The pseudocode of this procedure is collected in Algorithm <ref>. Recall $[\bfc{Y}]_1(\widehat{V}_d^{(t-1)} \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3}^{(t-1)} \otimes I_{p_2})\widehat{V}_{2}^{(t-1)}$ is the intermediate product from the $(t-1)$st update. We sequentially compute

$$\widehat{U}^{(t)}_1 = \SVD^L \left([\bfc{Y}]_1(\widehat{V}_d^{(t - 1)} \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3}^{(t-1)} \otimes I_{p_2})\widehat{V}_{2}^{(t-1)}\right);$$

$$\widehat{U}_k^{(t)} = \SVD^L\left((I_{p_k} \otimes \widehat{U}^{(t)\top}_{k - 1})\cdots (I_{p_2\cdots p_k} \otimes \widehat{U}^{(t)\top}_{1})[\bfc{Y}]_k(\widehat{V}_d^{(t-1)} \otimes I_{p_{k+1}\dots p_{d-1}})\cdots (\widehat{V}_{k + 2}^{(t-1)} \otimes I_{p_{k+1}})\widehat{V}_{k + 1}^{(t - 1)}\right)$$

for $k=2,\ldots, d-1$, and

$$\widehat{U}_d^{(t)} = [(\widehat{U}_{d-1}^{(t)})^\top(I_{p_{d-1}} \otimes (\widehat{U}_{d-2}^{(t)})^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes (\widehat{U}_{1}^{(t)})^\top)[\bfc{Y}]_{d-1}]^\top \in \bbR^{p_d \times r_{d-1}}.$$

Reshape $\widehat{U}_k^{(t)} \in \bbR^{(p_kr_{k-1}) \times r_k}$ to $\widehat{\bfc{U}}_k^{(t)} \in \bbR^{r_{k-1} \times p_k \times r_k}$ for $k = 2, \dots, d - 1$.
Then, compute

$$\widehat{\bfc{X}}^{(t)} = \llbracket\widehat{U}_1^{(t)}, \widehat{\bfc{U}}_2^{(t)}, \dots, \widehat{\bfc{U}}_{d-1}^{(t)}, \widehat{U}_{d}^{(t)}\rrbracket.$$

We will explain the algebraic schemes in the TTOI procedure through several representation lemmas in Section <ref>. We will also show in Theorem <ref> that the objective function $\|\bfc{Y} - \widehat{\bfc{X}}^{(t)}\|_{\F}^2$ is monotone decreasing with respect to the iteration index $t$. In the large-scale scenarios that performing iterations is beyond the capacity of computing, we can reduce the number of iterations, and even to $t_{\max}=1$, i.e., the one-step iteration, which have often yielded sufficiently accurate estimation as we will illustrate in both theory and simulation studies. Such the phenomenon has been recently discovered for HOOI in the Tucker low-rank tensor decomposition [50].

[Computational and storage costs of TTOI]
We consider the computational and storage costs of TTOI on the $p$-dimensional, rank-$r$, order-$d$, and dense tensor. Since computing the first $r$ singular vectors of an $m \times n$ matrix via block power method requires $\widetilde{O}(mnr)$ operations, initialization costs $\widetilde{O}(p^dr)$ operations, each iteration of TTOI, including  forward and backward updates, costs $O(p^dr)$. Therefore, the total number of operations of TTOI  with $T$ iterations is $\widetilde{O}(p^dr) + O(Tp^dr)$, which is not significantly more than the number of elements of the target tensor. Moreover, TTOI requires $O(p^d)$ storage cost, which is not significantly more than the storage cost of the original tensor.

§ THEORETICAL ANALYSIS

This section is devoted to the theoretical analysis of the proposed procedure. For convenience, we introduce the following two abbreviations for matrix sequential products: for $M_i \in \bbR^{(p_ir_{i-1}) \times r_i}, 1 \leq i \leq d-1$ and $B_j \in \bbR^{(r_jp_j) \times r_{j-1}}, 2 \leq j \leq d$, we denote

$$M_{\text{prod}, k}^{(L)} = (I_{p_2\cdots p_k} \otimes M_1)\cdots(I_{p_k} \otimes M_{k-1})M_k \in \bbR^{(p_1\cdots p_k) \times r_k}, \quad \forall 1 \leq k \leq d-1,$$

$$B_{\text{prod}, k}^{(R)} = (B_d \otimes I_{p_k\cdots p_{d-1}})\cdots(B_{k+1} \otimes I_{p_k})B_k \in \bbR^{(p_k\cdots p_d) \times r_{k-1}}, \quad \forall 2 \leq k \leq d.$$

Equivalently, $M_{\text{prod}, k}^{(L)}$ and $B_{\text{prod}, k}^{(R)}$ can be defined sequentially as

\begin{equation*}
M_{\text{prod}, 1}^{(L)} = M_1, \quad M_{\text{prod}, k+1}^{(L)} = (I_{p_{k+1}} \otimes M_{\text{prod}, k}^{(L)})M_{k+1}, \quad 1 \leq k \leq d-2,
\end{equation*}

\begin{equation*}
B_{\text{prod}, d}^{(R)} = B_d, \quad B_{\text{prod}, k}^{(R)} = (B_{\text{prod}, k+1}^{(R)} \otimes I_{p_k})B_k, \quad 2 \leq k \leq d-1.
\end{equation*}

§.§ Representation Lemmas for high-order tensors

Since the computation of high-order tensors with tensor-train structures involves extensive tensor algebra, we introduce the following three lemmas on the matrix representation of high-order tensors. These lemmas play a fundamental role in the later theoretical analysis.

Suppose $\bfc{X} = \llbracket G_1,\bfc{G}_2, \ldots, \bfc{G}_{d-1}, G_d \rrbracket$. Then the sequential matricization of $\bfc{X}$ can be written as

\begin{equation}
\begin{split}
[\bfc{X}]_k = &  (I_{p_2\cdots p_k} \otimes G_1)(I_{p_3\cdots p_k} \otimes \left[\bfc{G}_2\right]_2)\cdots (I_{p_k} \otimes \left[\bfc{G}_{k-1}\right]_2)\left[\bfc{G}_k\right]_2\left[\bfc{G}_{k+1}\right]_1\\
& \cdot \left(\left[\bfc{G}_{k+2}\right]_1 \otimes I_{p_{k+1}}\right)\cdots\left([\bfc{G}_{d-1}]_1 \otimes I_{p_{k+1}\cdots p_{d-2}}\right)\left(G_d^\top \otimes I_{p_{k+1}\cdots p_{d-1}}\right).
\end{split}
\end{equation}

For any tensor $\bfc{T} \in \RR^{\otimes_{k=1}^d p_k}$ and $1 \leq i < j \leq d-1$, we have

\begin{equation*}
[\bfc{T}]_j = (I_{p_{i+1}\cdots p_j} \otimes [\bfc{T}]_i)A^{(p_{i+1}\cdots p_j, p_{j+1}\cdots p_d)},\quad
[\bfc{T}]_i = A^{(p_{i+1}\cdots p_j,p_1\cdots p_i)\top}([\bfc{T}]_j \otimes I_{p_{i+1}\cdots p_{j}}).
\end{equation*}

Here, we define $e_{k}^{(ij)}$ as the $k$th canonical basis of $\bbR^{ij}$ and

\begin{equation}\label{eq:realignment}
A^{(i, j)} = \begin{bmatrix}
e_{1}^{(ij)} & e_{i+1}^{(ij)} & \cdots & e_{i\left(j - 1\right) + 1}^{(ij)}\\
e_{2}^{(ij)} & e_{i+2}^{(ij)} & \cdots & e_{i(j-1) + 2}^{(ij)}\\
\vdots & \vdots & \ddots & \vdots\\
e_{i}^{(ij)} & e_{2i}^{(ij)} & \cdots & e_{ij}^{(ij)}
\end{bmatrix} \in \bbR^{(i^2j) \times j}.
\end{equation}

Lemmas <ref> and <ref> can be proved by checking each entry of the corresponding matricizations. In addition, the following lemma provides a representation of sequential reshaping tensor, in particular for $R_k^{(t)}$ and $\widetilde{R}_k^{(t)}$, the key intermediate outcomes in TTOI procedure.

Suppose $\bfc{T}\in \mathbb{R}^{\otimes_{k=1}^d p_k}, M_i \in \bbR^{(r_{i-1}p_i) \times r_i}$ for $1 \leq i \leq d-1$, $B_i \in \bbR^{(p_ir_i) \times r_{i-1}}$ for $2 \leq i \leq d$, where $r_0 =  r_d = 1$. Consider the following sequential multiplication:

Forward sequential multiplication: Let $S_1=[\bfc{T}]_1$. For $k=1,\ldots, d-1$, calculate

\begin{equation*}
\begin{split}
& \widetilde{S}_k = M_k^\top S_k \in \bbR^{r_k \times (p_{k+1}\cdots p_d)},\\
& S_{k+1} = \mathrm{Reshape} (\widetilde{S}_k, r_kp_{k+1}, p_{k+2}\cdots p_d) \quad \text{if } k < d-1.
\end{split}
\end{equation*}

Then for any $1 \leq k \leq d-1$,

\begin{equation}\label{eq:realignment_forward}
\begin{split}
&S_k = (I_{p_k} \otimes M_{\text{prod}, k-1}^{(L)\top})[\bfc{T}]_k, \quad
\widetilde{S}_k = M_{\text{prod}, k}^{(L)\top}[\bfc{T}]_k.
\end{split}
\end{equation}

Here, $I_{p_k} \otimes M_{\text{prod}, k-1}^{(L)\top} = I_{p_1}$ if $k = 1$.

Backward sequential multiplication: Let $W_{d-1} = [\bfc{T}]_{d-1}$. For $k=d-1,\ldots,1$, calculate

\begin{equation*}
\begin{split}
&\widetilde{W}_k = W_kB_{k+1} \in \bbR^{(p_1\cdots p_k) \times r_k}, \\
&W_{k-1} = \mathrm{Reshape} (\widetilde{W}_k, p_1\cdots p_{k-1}, p_{k}r_{k})  \quad \text{if } k > 1.
\end{split}
\end{equation*}

Then for any $1 \leq k \leq d-1$,

\begin{equation*}
\begin{split}
&W_k = [\bfc{T}]_k(B_{\text{prod}, k+2}^{(R)} \otimes I_{p_{k+1}}), \quad
\widetilde{W}_k = [\bfc{T}]_kB_{\text{prod}, k+1}^{(R)}.
\end{split}
\end{equation*}

Here, $B_{\text{prod}, k+2}^{(R)} \otimes I_{p_{k+1}} = I_{p_d}$ if $k = d-1$.

In particular, $R_k^{(0)}, \widetilde{R}_k^{(0)}$ in Algorithm <ref> and $R_k^{(t)}, \widetilde{R}_k^{(t)} (t \in \{2,4,6,\ldots\})$ in Algorithm <ref> satisfy

\begin{equation}\label{eq1}
\begin{split}
&R_k^{(t)} = \left(I_{p_k} \otimes (\widehat{U}^{(t)})_{\text{prod}, k-1}^{(L)\top}\right)[\bfc{Y}]_k, \quad \widetilde{R}_k^{(t)} = (\widehat{U}^{(t)})_{\text{prod}, k}^{(L)\top}[\bfc{Y}]_k, \quad \forall 1 \leq k \leq d-1.
\end{split}
\end{equation}

The proof of Lemma <ref> is provided in Section <ref>.

§.§ Deterministic Upper Bounds for Estimation Error of TTOI

Now we are in position to analyze the performance of TTOI. The following Theorem <ref> introduces an upper bound on estimation error of $\widehat{\bfc{X}}^{(2t+1)}$ (backward update) and $\widehat{\bfc{X}}^{(2t+2)}$ (forward update).

Suppose we observe $\bfc{Y} = \bfc{X} + \bfc{Z}$, where $\bfc{X}$ admits a TT decomposition as (<ref>).

(A deterministic estimation error bound for backward updates) Let $\widetilde{U}_1^{(2t)} = U_1 \in \bbR^{p_1 \times r_1}$ be the left singular space of $[\bfc{X}]_1$. For $2 \leq k \leq d-1$, define $\widetilde{U}_k^{(2t)} \in \bbR^{p_kr_{k-1} \times r_k}$ as the left singular subspace of $\left(I_{p_k} \otimes (\widehat{U}^{(2t)})_{\text{prod}, k-1}^{(L)\top}\right)[\bfc{X}]_k$. If for some constant $c_0 \in (0, 1)$,

\begin{equation}\label{ineq:initial_condition}
\begin{split}
& \left\|\sin\Theta\left(\widehat{U}_k^{(2t)}, \widetilde{U}_k^{(2t)}\right)\right\| \leq c_0, \quad \forall 1 \leq k \leq d-1,
\end{split}
\end{equation}

then there exists a constant $C_d > 0$ that only depends on $d$ such that the outcome of Algorithm <ref> satisfies

\begin{equation}\label{ineq:upper_deterministic}
\begin{split}
\left\|\widehat{\bfc{X}}^{(2t+1)} - \bfc{X}\right\|_{\F}^2 \leq C_d\left(\sum_{k=1}^{d-1}A_k^{(2t+1)} + B^{(2t+1)}\right),	\end{split}
\end{equation}

$$A_k^{(2t+1)} = \left\|(\widehat{U}^{(2t)})_{\text{prod}, k}^{(L)\top}[\bfc{Z}]_{k}\left((\widehat{V}^{(2t+1)})_{\text{prod}, k+2}^{(R)} \otimes I_{p_{k+1}}\right)\right\|_{\F}^2,$$

$$B^{(2t+1)}=\left\|[\bfc{Z}]_1(\widehat{V}^{(2t+1)})_{\text{prod}, 2}^{(R)}\right\|_{\F}^2.$$

Here, $(\widehat{V}^{(2t+1)})_{\text{prod}, k+2}^{(R)} \otimes I_{p_{k+1}} = I_{p_d}$ if $k = d-1$.

(A deterministic estimation error bound for forward updates) For $2 \leq k \leq d - 1$, let $\widetilde{V}_{k}^{(2t+1)} \in \bbR^{(p_kr_k) \times r_{k-1}}$ be the right singular space of $[\bfc{X}]_{k-1}\left((\widehat{V}^{(2t+1)})_{\text{prod}, k+1}^{(R)} \otimes I_{p_k}\right)$ and let $\widetilde{V}_{d}^{(2t+1)} = V_d \in \bbR^{p_d \times r_{d-1}}$ be the right singular space of $[\bfc{X}]_{d-1}$. If for some constant $c_0 \in (0, 1)$,

\begin{equation*}
\begin{split}
& \left\|\sin\Theta\left(\widehat{V}_k^{(2t+1)}, \widetilde{V}_k^{(2t+1)}\right)\right\| \leq c_0, \quad \forall 2 \leq k \leq d,
\end{split}
\end{equation*}

then there exists a constant $C_d > 0$ that only depends on $d$ such that the outcome of Algorithm <ref> satisfies

\begin{equation}\label{ineq:upper_deterministic_X^2t+2}
\left\|\widehat{\bfc{X}}^{(2t+2)} - \bfc{X}\right\|_{\F}^2 \leq C_d\left(\sum_{k=1}^{d-1}A_k^{(2t+2)} + B^{(2t+2)}\right),
\end{equation}

$$A_k^{(2t+2)} = \left\|\left(I_{p_{k}} \otimes (\widehat{U}^{(2t+2)})_{\text{prod}, k-1}^{(L)\top}\right)[\bfc{Z}]_{k}(\widehat{V}^{(2t+1)})_{\text{prod}, k+1}^{(R)}\right\|_{\F}^2,$$

$$B^{(2t+2)}=\left\|(\widehat{U}^{(2t+2)})_{\text{prod}, d-1}^{(L)\top}[\bfc{Z}]_{d-1}\right\|_{\F}^2.$$

Here, $I_{p_{k}} \otimes (\widehat{U}^{(2t+2)})_{\text{prod}, k-1}^{(L)\top} = I_{p_{1}}$ if $k = 1$.

The proof of Theorem <ref> is provided in Section <ref>. Theorem <ref> shows the estimation error $\|\widehat{\bfc{X}}^{(t+1)} - \bfc{X}\|_{\F}^2$ can be bounded by the projected noise $\bfc{Z}$, i.e., $A^{(t+1)}_k$ and $B^{(t+1)}$, if the estimates in initialization ($t=0$) or the previous iteration $(t\geq 1)$,  $\{\widehat{U}_k^{(t)}\}_{k=1}^{d-1}$ or $\{\widehat{V}_k^{(t)}\}_{k=2}^d$, are within constant distance to the true underlying subspaces.
The developed upper bound can be significantly smaller than $C\left\|\bfc{Z}\right\|_{\rm F}^2$, the classic upper bound induced from the approximation error (e.g., Theorem 2.2 in [62]), especially in the high-dimensional setting ($p\gg r$).

[Interpretation of error bounds in Theorem <ref>]
Here, we provide some explanation for $A_{k}^{(2t+1)}$ and $B^{(2t+1)}$ in the error bound (<ref>). By algebraic calculation, the TT-core estimation via backward update can be written as

\begin{align*}
\begin{split}
\widehat{V}_{k+1}^{(2t+1)} =&  SVD^R\Big\{(\widehat{U}^{(2t)})_{\text{prod}, k}^{(L)\top}([\bfc{X}]_k + [\bfc{Z}]_k)\Big((\widehat{V}^{(2t+1)})_{\text{prod}, k+2}^{(R)} \otimes I_{p_{k+1}}\Big)\Big\}, \quad \forall 1 \leq k \leq d-1
\end{split}
\end{align*}

\begin{equation*}
\widehat{V}_1^{(2t+1)} = ([\bfc{X}]_1 + [\bfc{Z}]_1)(\widehat{V}^{(2t+1)})_{\text{prod}, 2}^{(R)}.
\end{equation*}

From the definition of $A_k^{(2t+1)}$, we have see  $A_k^{(2t+1)}$ quantifies the error of the singular subspace estimate  $\widehat{V}_{k+1}^{(2t+1)}$ and $B^{(2t+1)}$ quantifies the error of the projected residual $\widehat{V}_1^{(2t+1)}$. By symmetry, similar interpretation also applies to $A_k^{(2t+2)}$ and $B^{(2t+2)}$ for the error bound of forward update (<ref>).

[Proof Sketch of Theorem <ref>]
While the complete proof of Theorem <ref> is provided in Section <ref>, we provide a brief proof sketch here.

Without loss of generality, we focus on (<ref>) for $t = 0$ while other cases follows similarly. For convenience, we simply let $\widehat{U}_i, \widehat{V}_i$ denote $\widehat{U}_i^{(0)}, \widehat{V}_i^{(1)}$, respectively. First, by Lemma <ref>, we can transform $[\widehat{X}^{(1)}]_1$, the outcome of backward update, to

$$[\widehat{X}^{(1)}]_1 = [\bfc{Y}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}}.$$

Then we can further bound the estimation error of $\widehat{\bfc{X}}^{(1)}$ as

\begin{equation*}
\begin{split}
\|\widehat{\bfc{X}}^{(1)}-\bfc{X}\|_\F^2 \leq & C \left\|[\bfc{Z}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}\right\|_{\F}^2\\
& + C_d\sum_{k=2}^{d}\left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}})\right\|_{\F}^2.
\end{split}
\end{equation*}

Next, based on Lemma <ref> and (<ref>), we can prove

\begin{equation*}
\begin{split}
&\left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}})\right\|_{\F}\\
=&\left\|[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& C_d\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}.
\end{split}
\end{equation*}

Finally, we apply the perturbation projection error bound (Lemma <ref>) to prove that

\begin{equation*}
\begin{split}
&C_d\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& C_d\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{Z}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\right\|_{\F}.
\end{split}
\end{equation*}

Theorem (<ref>) is proved by combing all inequalities above.

Next, we establish a decomposition formula for the approximation error, i.e., the objective function in (<ref>) $\left\|\bfc{Y} - \bfc{X}^{(t)}\right\|_{\F}^2$, and show that the approximation error is monotone decreasing through TTOI iterations.

We implement TTOI on $\bfc{Y}$. Let $\widehat{\bfc{X}}^{(t)}$ be the outcome after the $t$th iteration. For any $k\geq 1$, we have

\begin{equation}\label{ineq:convergence}
\text{(Approximation error decay)}\quad 	\|\bfc{Y}\|_\F^2 - \|\widehat{\bfc{X}}^{(t+1)}\|_{\F}^2 \leq \|\bfc{Y}\|_\F^2 - \|\widehat{\bfc{X}}^{(t)}\|_{\F}^2,
\end{equation}

\begin{equation}
\text{(Approximation error decomposition)} \quad \|\bfc{Y} - \widehat{\bfc{X}}^{(t+1)}\|_{\F}^2 = \|\bfc{Y}\|_\F^2 - \|\widehat{\bfc{X}}^{(t+1)}\|_{\F}^2.
\end{equation}

See Section <ref> for the proof of Theorem <ref>.

§ TTOI FOR TENSOR-TRAIN SPIKED TENSOR MODEL

In this section, we further focus on a probabilistic setting, spiked tensor model, where the noise tensor $\bfc{Z}$ has independent, mean zero, and $\sigma$-sub-Gaussian entries (see definition in Section <ref>). The spiked tensor model has been widely studied as a benchmark setting for tensor PCA/SVD and dimension reduction in recent literature in machine learning, information theory, statistics, and data science [39, 69, 64, 84, 88]. The central goal therein is to discover the underlying low-rank tensor $\bfc{X}$. Most of the existing works focused on tensors with Tucker or CP decomposition.

Under the spiked tensor model, we can verify that the initialization step of TTOI gives sufficiently good initial estimations with high probability that matches the required condition in Theorem <ref>.

Suppose $\bfc{X}$ is TT-decomposable as (<ref>) and $\bfc{Z}$ have independent zero mean and $\sigma$-sub-Gaussian random variables. Denote  $p = \min\{p_1, \cdots,p_d\}$. If there exists a constant $C_{gap}$ such that $\lambda_k = s_{r_k}([\bfc{X}]_k) \geq C_{gap}\left((\sum_{i=1}^{d}p_ir_{i-1}r_i)^{1/2} + (p_{k+1}\cdots p_d)^{1/2}\right)\sigma$ for $1 \leq k \leq d-1$, then there exist some constants $C, c > 0$ and $C_d > 0$ that only depends on $d$, with probability at least $1 - C\exp(-cp)$,

\begin{equation}\label{ineq36}
\max_{k=1,\ldots, d-1} \left\|\sin\Theta\left(\widehat{U}_k^{(0)}, \widetilde{U}_k^{(0)}\right)\right\| \leq \frac{1}{2},
\end{equation}

\begin{equation}\label{ineq37}
\max_{\substack{k=1,\ldots, d-1\\t=2, 4, 6, \ldots}} \left\|\sin\Theta\left(\widehat{U}_k^{(t)}, \widetilde{U}_k^{(t)}\right)\right\|, \max_{\substack{k=2,\ldots, d\\t=1,3, 5, \ldots}} \left\|\sin\Theta\left(\widehat{V}_k^{(t)}, \widetilde{V}_k^{(t)}\right)\right\|
\leq \frac{1}{2},
\end{equation}
and for all $t \geq 1$,

\begin{equation}\label{ineq35}
\begin{split}
\max\{A_k^{(t)}, B^{(t)}\} \leq C_d\sigma^2 \sum_{i=1}^d p_i r_ir_{r-1}.
\end{split}
\end{equation}

Here, $\widetilde{U}_k^{(t)}$, $\widetilde{V}_k^{(t)}$, $A_k^{(t)}$ and $B^{(t)}$ are defined in Theorem <ref>.

The proof of Theorem <ref> is provided in Section <ref>. Based on Theorems <ref> and <ref>, we can further prove:

Suppose $\bfc{X}$ can be decomposed as (<ref>), $\bfc{Z}_{i_1, \dots, i_d}$ are independent zero mean and $\sigma$-sub-Gaussian random variables, $p = \min\{p_1, \cdots,p_d\}$. Suppose there exists a constant $C_{gap}$ such that $\lambda_k = s_{r_k}([\bfc{X}]_k) \geq C_{gap}\left((\sum_{i=1}^{d}p_ir_{i-1}r_i)^{1/2} + (p_{k+1}\cdots p_d)^{1/2}\right)\sigma$ for $1 \leq k \leq d-1$. Then with probability at least $1 - Ce^{-cp}$, for all $t \geq 1$,

\begin{equation}\label{ineq:upper bound-spiked-model}
\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\|_{\F}^2 \leq  C_d\sigma^2\sum_{i=1}^{d}p_ir_ir_{i-1}.
\end{equation}

The proof of Corollary <ref> is provided in Section <ref>.
[Interpretation of Corollary <ref>]
Note that the TT-cores $G_1, \bfc{G}_i, G_d$ respectively have $p_1r_1,  p_ir_ir_{i-1}, p_dr_{d-1}$ free parameters, the upper bound (<ref>) can be seen as the noise level $\sigma^2$ times the degrees of freedom of the low TT rank tensors.

Next, we develop a minimax lower bound for the low TT rank structure estimation. Consider the following general class of tensors with dimension $\bp = (p_1, \dots, p_d)$ and TT rank $\br = (r_1, \dots, r_{d-1})$,

\begin{equation}\label{lower:class}
\mathcal{F}_{\bp, \br}(\bm\lambda) = \left\{\bfc{X} \in \bbR^{p_1 \times \cdots \times p_d}, \bfc{X} \text{ can be decomposed as } \eqref{TT}, s_{r_k}\left([\bfc{X}]_k\right) \geq \lambda_k, 1 \leq k \leq d-1\right\},
\end{equation}

and a class of distributions of $\sigma$-sub-Gaussian noise tensors

\begin{equation}\label{eq:distribution-class-D}
\mathcal{D} = \{D: \text{ if } \bfc{Z} \sim D, \text{ then } \bfc{Z}_{i_1, \dots, i_d} \text{ are indep. zero mean and } \sigma \text{ sub-Gaussian random variables}\}.
\end{equation}

Here, the constraints on the least singular value of $[\bfc{X}]_k$ and the $\sigma$-sub-Gaussian assumption correspond to the conditions required for upper bound in Theorem <ref>.

Consider the order-$d$ TT spiked tensor model (<ref>) and distribution class $\mathcal{D}$ in (<ref>).
Assume $p = \min\{p_1, \dots, p_d\} \geq C_0$ for some large constant $C_0$, $r_1\leq p_1/2, r_i\leq p_ir_{i-1}/2, r_{i-1}\leq p_ir_{i}/2$ for $2 \leq i \leq d-1$, $r_{d-1} \leq p_d$, and $\lambda_i > 0$. Also assume $r_1r_2 \leq p_1$ if $d = 3$. Then there exists a constant $c_d > 0$ that only depends on $d$ such that

\begin{equation}\label{ineq:lower_bound_Gaussian}
\inf_{\widehat{\bfc{X}}}\sup_{\bfc{X} \in \mathcal{F}_{\bp, \br}(\bm\lambda), D \in \mathcal{D}}\bbE_{\bfc{Z} \sim D}\left\|\widehat{\bfc{X}} - \bfc{X}\right\|_{\F}^2 \geq c_d\sigma^2\sum_{i=1}^{d}p_ir_ir_{i-1}.
\end{equation}

See Section <ref> for the proof of Theorem <ref>.

§ TTOI FOR DIMENSION REDUCTION AND STATE AGGREGATION IN HIGH-ORDER MARKOV CHAIN

Since the introduction at the beginning of the 20th century, the Markov process has been ubiquitous in a variety of disciplines. In the literature, the first order Markov process, i.e., the future observation at $(t+1)$ is conditionally independent of those at times $1, \ldots, (t-1)$ given the immediate past observation at time $t$, has been commonly used and extensively studied. Moreover, the high-order Markov process often appear in many scenarios, where the future observation is affected by a longer history. For example, in the taxi travel trajectory, the future stop of a taxi not only depends on the current location but also the past path that reveals the direction this taxi is heading to [7]. The high-order Markov processes have also been applied to inter-personal relationship [66], financial econometrics [78], traffic flow [92], among many other applications.

We specifically consider an ergodic, time-invariant, and $(d-1)$st order Markov process on a finite state space $\{1, \ldots, p\}$. That is, the future state $X_{t+d}$ depends on the current state $X_{t+d-1}$ and the previous $(d - 2)$ states $(X_{t+d-2}, \ldots, X_{t+1})$ jointly:

\begin{equation}\label{eq:transition-tensor}
\bbP\left(X_{t+d}| X_1, \ldots, X_{t+d-1}\right) = \bbP\left(X_{t+d}| X_{t + 1}, \ldots, X_{t+d-1}\right) = \bfc{P}_{[X_{t+1}, \ldots, X_{t+d}]}.
\end{equation}

Our goal is to achieve a reliable estimation of the transition tensor $\bfc{P}$ and to predict the future state $X_{t+d}$ based on an observable trajectory. Since the total number of free parameters in a $(d-1)$st order Markov transition tensor $\bfc{P}$ is $O(p^d)$ without further assumptions, it may be prohibitively difficult to infer $\bfc{P}$ in both statistics and computation even if $p$ and $d$ are only of moderate scale. Instead, a sufficient dimension reduction for high-order Markov processes is in demand.

To enable the statistical inference and dimension reduction for high-order Markov processes, a powerful tool, mixed transition distribution model (MTD), was introduced [66]. The MTD model assumes that the distribution of future state is a linear combination of the distributions associated with the $(d-1)$ immediate past states.
The readers are also referred to [8] for a survey on mixed transition distribution model. The linear assumption, however, does not take into account the potential interactions of past states that commonly appear in practice. For example in the New York taxi trip data, the interaction among past locations of a taxi indicates its potential future direction.

On the other hand, there is a recent surge of development in dimension reduction and state aggregation for first order Markov chains. For example, [31] considered the Markov chain aggregation and the application to biology; [26] considered the rank-reduced Markov model and mode clustering; [87] considered Markov rank, aggregagability, and lumpability of Markov processes and proposed the dimension reduction and state aggregation methods through spectral decomposition with theoretical guarantees; [70] proposed clustering block model and proposed efficient algorithm to solve it; [96] introduced a convex and non-convex methods to estimate the rank-reduced low-rank Markov transition matrix.

Inspired by these work, we propose and study the state aggregation model for the discrete-time high-order Markov processes as follows.
[$(d-1)$st order state aggregatable Markov process]
Suppose there exist maps $G_1: [p]\to \mathbb{R}^{r_1}$, $G_k: [p]\times \mathbb{R}^{r_{k-1}} \to \mathbb{R}^{r_k}$, $G_d: [p]\times\mathbb{R}^{r_{d-1}}\to \mathbb{R}$ such that $G_2,\ldots, G_d$ are linear: $G_k(X, \lambda_1 u+\lambda_2 v) = \lambda_1 G_k(X, u) + \lambda_2 G_k(X, v)$ for any vectors $u, v$, scalars $\lambda_1, \lambda_2\in \mathbb{R}$. We say a Markov process $\{X_1,X_2,\ldots\}$ is $(d-1)$st order state aggregatable if for all $t\geq 0$, the transition can be sequentially generated as follows,

\begin{equation*}
\begin{split}
&\widetilde{P}_1(X_{t+1}) =  G_1(X_{t+1}) \in \bbR^{r_1},\\
&\widetilde{P}_k(X_{t+1}, \dots, X_{t+k}) = G_k(X_{t+k}, \widetilde{P}_{k-1}(X_{t+1}, \dots, X_{t+k-1})) \in \bbR^{r_k}, \quad k= 2,\ldots, d-1,\\
& \bbP\left(X_{t+d} | X_1, \dots, X_{t+d-1}\right) = \bbP\left(X_{t+d} | X_{t+1}, \dots, X_{t+d-1}\right) = G_d(X_{t+d}, \widetilde{P}_{d-1}(X_{t+1}, \dots, X_{t+d-1})).
\end{split}
\end{equation*}

In a $(d-1)$st order state aggregatable Markov process, the future state $X_{t+d}$ relies on a sequential aggregation of the previous $d-1$ states $X_{t+1},\ldots, X_{t+d-1}$ as follows: we first project $X_{t+1}$ to a $r_1$-dimensional vector $\widetilde{P}_1(X_{t+1})$ via $G_1$, then project $\widetilde{P}_1(X_{t+1})$ jointly with $X_{t+2}$ to a $r_2$-dimensional vector $\widetilde{P}_1(X_{t+1}, X_{t+2})$ via $G_2$. We repeat such the projection sequentially for $X_{t+3}, \ldots, X_{t+d}$ and yield the transition probability $\mathbb{P}\left(X_{t+d}|X_{t+1},\ldots, X_{t+d-1}\right)$. Also, see Figure <ref> for a pictorial illustration.

Based on the definition of the state aggregatable Markov chain, we can prove the corresponding probability transition tensor $\bfc{P}$ will have low TT rank.

The transition tensor $\bfc{P}$ of the rank reduced high-order Markov model in Definition <ref> has TT-rank no more than $(r_1,\dots,r_{d-1})$. In other words, $\bfc{P}$ satisfies $\rank([\bfc{P}]_k)\leq r_k$.

The proof of Proposition <ref> is provided in Section <ref>.

A pictorial illustration of a $(d-1)$st order state aggregatable Markov chain

Next, we focus on a synchronous or generative setting, which can be seen as a high-order generalization of the classic observation model for the analysis of Markov (decision/reward) processes (see [36] for an introduction), for the high-order Markov process. To be specific, for each sample index $k = 1, \dots, n$ and previous states $(i_1, \dots, i_{d-1}) \in [p]^{d-1}$, suppose we observe the next state $X(i_1, \dots, i_{d-1}; k)$ drawn from the Markov transition tensor $\bfc{P}$. It is natural to estimate $\bfc{P}$ via the empirical transition tensor:

\begin{equation*}
\widehat{\bfc{P}}^{\rm{emp}}_{i_1,\dots,i_d}=\sum_{k=1}^{n}1_{\{X(i_1, \dots, i_{d-1}; k) = i_d\}}\Big/n, \quad i_1, \dots, i_d \in \{1, \dots, p\}^{d}.
\end{equation*}

Then, $\widehat{\bfc{P}}^{\rm{emp}}$ is an unbiased estimator of $\bfc{P}$. However, if the entries of $\bfc{P}$ are approximately balanced, the mean squared error of $\widehat{\bfc{P}}^{\rm{emp}}$ satisfies

\begin{equation}\label{ineq:P^emp}
\begin{split}
& \bbE\left\|\widehat{\bfc{P}}^{\rm{emp}} - \bfc{P}\right\|_{\rm F}^2 = \sum_{i_1, \dots, i_d}\Var\left(\widehat{\bfc{P}}_{i_1, \dots, i_d}^{\rm{emp}}\right) \\
= & \sum_{i_1, \dots, i_{d-1}}\sum_{i_d}\frac{\bbP\left(i_d|i_1, \dots, i_{d-1}\right)\left(1 - \bbP\left(i_d|i_1, \dots, i_{d-1}\right)\right)}{n} \asymp \frac{p^{d-1}}{n},
\end{split}
\end{equation}

To obtain a more accurate estimator, we propose to first perform TTOI on $\widehat{\bfc{P}}^{\rm emp}$ to obtain $\widehat{\bfc{P}}^{(1)}$, then project each row of $[\widehat{\bfc{P}}^{(1)}]_{d-1}$, or equivalently, each mode-$d$ fiber of $\widehat{\bfc{P}}^{(1)}$, onto the simplex $S^{p-1} = \{x \in \bbR^{p}: \sum_{i=1}^{p}x_i = 1, x_i \geq 0 \text{ for all } 1 \leq i \leq p\}$ via probability simplex projection (see an implementation in [28]) and obtain $\widehat{\bfc{P}}$.

We establish an upper bound on estimation error for the TTOI estimator $\widehat{\bfc{P}}$.

Consider the synchronous or generative model for a $(d-1)$st order state aggregatable Markov process described above. Suppose the initialization condition (<ref>) in Theorem <ref> holds. Then with probability at least $1 - Ce^{-cp}$, the output of one-step TTOI followed by the probability simplex projection satisfies

\begin{equation*}
\left\|\widehat{\bfc{P}} - \bfc{P}\right\|_{\rm F}^2 \leq C\left(\max_{1 \leq i \leq d-1}r_i\right)\sum_{i=1}^dp_ir_ir_{i-1}\Big/n.
\end{equation*}

The proof of Proposition <ref> is provided in Section <ref>.
Compared to the estimation error rate of $\widehat{\bfc{P}}^{\rm emp}$ in (<ref>), Proposition <ref> shows TTOI achieves significantly reduced estimation error by exploiting the low TT rank structure of the high-order Markov process.

If the observations form one transition trajectory $\{X_0,\dots,X_{N}\}$, we can work on the following empirical transition tensor

\begin{equation}\label{eq:MC_empirical_est}
\widehat{\bfc{P}}^{\rm{emp}}_{i_1,\dots,i_d}=\begin{cases}
\frac{\sum_{t=0}^{N-d+1} 1_{\{X_t=i_1,\dots,X_{t+d-1}=i_d\}}}{\sum_{t=0}^{N-d+1} 1_{\{X_t=i_1,\dots,X_{t+d-2}=i_{d-1}\}}},&\sum_{t=1}^{N-d+1} 1_{\{X_t=i_1,\dots,X_{t+d-2}=i_{d-1}\}}>0;\\
1/p,&\sum_{t=1}^{N-d+1} 1_{\{X_t=i_1,\dots,X_{t+d-2}=i_{d-1}\}}=0.
\end{cases}
\end{equation}

Then $\widehat{\bfc{P}}^{\rm{emp}}$ can be a nearly unbiased and strongly consistent estimator for $\bfc{P}$. When the Markov process is $(d-1)$st order state aggregatable, we can apply TTOI to obtain a better estimate. As will be explored by numerical studies in Section <ref>, the TTOI estimator achieves favorable performance on the estimation of $\bfc{P}$.

§ NUMERICAL STUDIES

In this section, we investigate the numerical performance of TTOI.

§.§ Simulation

In each simulation setting, we present the numerical results in both average estimation error (denoted by dots) and standard deviation (denoted by bars) based on 100 repetitions. We assume the true TT-ranks are known in the first three settings. Afterwards, we introduce a BIC-type data-driven scheme for TT-rank selection and present its numerical performance. All experiments are conducted by a quad-core 2.3 GHz Intel Core i5  processor.

We first consider the tensor-train spiked tensor model (<ref>) discussed in Section <ref>. Specifically, we randomly generate  $G_1,\bfc{G}_2,\dots,\bfc{G}_{d-1}, G_d$ with i.i.d. standard normal entries, and generate $\bfc{Z}$ with i.i.d. $\mathcal{N}(0,\sigma^2)$ or $\text{Unif}(-b,b)$ entries. Let $p_1=\dots=p_d=p$, $r_1=\dots=r_{d-1}=r$, and consider four settings: (1) $p = 100, d = 3, r = 1$; (2) $p = 50, d = 4, r = 1$; (3) $p = 20, d = 5, r = 1$; (4) $p = 20, d = 5, r = 2$. For varying values of $\sigma \in[1, 19]$ and $b \in [3, 30]$, we evaluate the estimation error $\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}$ of the TT-SVD and TTOI estimators with 1 or 2 iterations, i.e., $t_{\max}=0,1,2$. From the results summarized in Figure <ref> (normal noise) and Figure <ref> (uniform noise), we can see TTOI, even with one iteration, performs significantly better than TT-SVD, and the advantage becomes more significant as the noise level $\sigma, b$ grows. This suggests that the proposed TTOI is effective for high-order tensor SVD compared to the classic TT-SVD, especially when the observations are corrupted by substantial noise. Table <ref> summarizes the runtime of TT-SVD and TTOI, which suggests that the additional computational cost incurred by the backward and forward updates in TTOI is negligible compared to the runtime of the original TT-SVD.

Estimation error of TT-SVD and TTOI for high-order spiked tensor model. Here, $\bfc{Z} \overset{\text{i.i.d.}}{\sim}N(0, \sigma^2)$.

Estimation error of TT-SVD and TTOI for high-order spiked tensor model. Here, $\bfc{Z} \overset{\text{i.i.d.}}{\sim}{\rm Unif}(-b,b)$.

To understand the influence of TT-rank to the performance of the TT-SVD and TTOI estimators, we conduct numerical experiments under the spiked tensor model (<ref>) with $r_1=\cdots=r_{d-1}=r$ for various values of $r$. In particular, $G_1, \bfc{G}_2, \dots, \bfc{G}_{d-1}, G_d$ are still generated with i.i.d. standard normal entries, and $\bfc{Z}$ has i.i.d. $\mathcal{N}(0,\sigma^2)$ entries. Letting $p_1=\dots=p_d=p$, we consider two settings: (1) $p = 100, d = 3, \sigma=20$; (2) $p = 500, d = 3,\sigma=100$. For $r=1,\dots,10$, we evaluate the average estimation error $\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}$ of TT-SVD, TTOI with 1 iteration, and TTOI with 2 iterations (i.e., $t_{\max}=0,1,2$), and present the results in Figure <ref>. Figure <ref> suggests that the estimation errors increase as the rank increases, while TTOI with 1 or 2 iterations both performs better than TT-SVD. The improvement of TTOI over TT-SVD is more significant under larger $p$ or smaller $r$. An intuitive explanation for this phenomenon is as follows: the key idea of TTOI is to utilize the previous updates to reduce the dimension of the sequential unfolding $[\bfc{Y}]_k$ before performing singular value thresholding; such the dimension reduction is more significant for large $p$ or small $r$.

Estimation error of TT-SVD and TTOI for high-order spiked tensor model with varying TT-ranks

Setting   TT-SVD   TTOI ($t_{\max}=1$)    TTOI ($t_{\max}=2$)

$p = 100$, $d = 3$, $r=1$   0.3316 (0.0708)   0.3338 (0.0710)   0.3402 (0.0744)

$p = 50$, $d = 4$, $r=1$   1.1652 (0.1728)   1.1690 (0.1720)   1.2012 (0.1709)

$p = 20$, $d = 5$, $r=1$   0.7254 (0.0926)   0.7300 (0.0919)   0.7514 (0.0949)

$p = 20$, $d = 5$, $r=2$   0.6718 (0.1002)   0.6764 (0.1005)   0.7080 (0.1032)

Runtime (in seconds) of TT-SVD, TTOI with $1$ iteration, and TTOI with $2$ iterations under the high-order spiked tensor model with $\bfc{Z}\overset{\text{i.i.d.}}{\sim}N(0, 400)$. The mean runtime of $50$ independent replicates are presented and the standard deviations are listed in parentheses.

Next, we demonstrate the performance of TTOI on transition tensor estimation for the high-order state-aggregatable Markov chains studied in Section <ref>. We consider the $(d-1)$st order Markov chain on $p$ states. To generate the transition tensor $\bfc{P}$, we first draw $\widetilde{G}_1\in \bbR^{p\times r}, \widetilde{\bfc{G}}_2\in \bbR^{r\times p\times r},\dots, \widetilde{G}_d\in \bbR^{r\times p}$ with i.i.d. standard normal entries, then normalize the rows of $\widetilde{G}_1, \widetilde{\bfc{G}}_2,\dots, \widetilde{G}_d$ in absolute values as

$$G_{1, [i, j]} = \frac{|\widetilde{G}_{1, [i, j]}|}{\sum_{j'} |\widetilde{G}_{1, [i, j']}|}, \quad \bfc{G}_{k, [i_1, i_2, j]} = \frac{|\widetilde{\bfc{G}}_{k, [i_1, i_2, j]}|}{\sum_{j'}|\widetilde{\bfc{G}}_{k, [i_1, i_2, j']}|}, \quad G_{d, [i, j]} = \frac{|\widetilde{G}_{d, [i, j]}|}{\sum_{j'} |\widetilde{G}_{d, [i, j']}|}.$$

By this means, $\bfc{P} = \llbracket G_1, \bfc{G}_2,\ldots, \bfc{G}_{d-1}, G_d\rrbracket$ satisfies $\bfc{P}_{i_1,\ldots, i_d}\geq 0$, $\sum_{i_d=1}^p \bfc{P}_{i_1,\ldots, i_d} = 1$ for any $(i_1,\ldots, i_{d-1})$, so $\bfc{P}$ forms a Markov transition tensor. To generate the trajectory $\{X_1,\dots,X_{N}\}$, we generate the initial $d-1$ states $X_1,\dots,X_{d-1}$ i.i.d. uniformly from $[p]$, then generate $X_d,\dots, X_{N}$ sequentially according to (<ref>). To estimate $\bfc{P}$, we construct the empirical probability tensor $\widehat{\bfc{P}}^{\mathrm{emp}}$ by (<ref>), then apply TT-SVD and TTOI with input $\widehat{\bfc{P}}^{\mathrm{emp}}$ as detailed in Section <ref> to obtain  $\widehat{\bfc{P}}$. We consider two numerical settings: (1) $p=100, d=3, r=1$; (2) $p=50, d=4, r=1$. We evaluate the estimation error $\|\widehat{\bfc{P}}^{(i)}-\bfc{P}\|_{\F}$ for each setting and summarize the results to Figure <ref>. Again, TTOI exhibits clear advantage over the existing methods in all simulation settings.

Estimation error of the transition tensor versus length of the observable trajectory in high order state-aggregatable Markov chain estimation.

Selection of TT-ranks. The proposed TTOI algorithm requires specifying TT-ranks $r_1, \ldots, r_{d-1}$ as inputs and the appropriate choices of $r_1, \ldots, r_{d-1}$ are crucial in practice. We propose a data-driven scheme to select the TT-ranks: we choose $r_1,\dots,r_{d-1}\geq 1$ such that the following Bayesian information criterion (BIC) under the spiked tensor model is minimized:

\begin{equation}\label{eq:BIC}
\begin{split}
\mathrm{BIC}(r_1,\ldots, r_{d-1}):= & \prod_{k=1}^dp_k\log\|\bfc{Y}-\widehat{\bfc{X}}(r_1,\dots,r_{d-1})\|_F^2\\
& +\left(p_1r_1 + \sum_{k=2}^{d-1} p_kr_{k-1}r_k + p_dr_{d-1}\right)\left(\sum_{k=1}^d \log p_k\right).
\end{split}
\end{equation}

Here, $\widehat{\bfc{X}}(r_1,\dots,r_{d-1})$ is the output of TTOI (Algorithm <ref>) with the input TT-ranks $r_1,\dots,r_{d-1}$. This BIC-type criterion was also adopted in prior works on tensor clustering [33].

Then we conduct numerical experiments under the same setting as the bottom two plots in Figure <ref> on the spiked tensor model with Gaussian noise. Figure <ref> summarizes the estimation errors of TT-SVD and TTOI with $1$ and $2$ iterations, respectively, with the ranks selected based on the proposed BIC criterion (<ref>). Comparing Figure <ref> to the bottom two plots in Figure <ref>, we can see the proposed criterion can select the true ranks accurately and the performance of both TT-SVD and TTOI with tuned ranks is very similar to the one by inputting the true ranks.

Average estimation error of TT-SVD and TTOI for high-order spiked tensor model with BIC-tuned ranks.

§.§ Real Data Experiments

We apply the proposed method to investigate the Manhattan taxi data[2013 Trip Data, available at  <https://chriswhong.com/open-data/foil_nyc_taxi/>]. This dataset contains the New York City taxi trip records from 14,144 drivers in 2013. We treat each travel record as a transition among different locations at New York City, then the overall dataset can be organized as a collection of fragmented sample trajectories of a Markov chain on New York City traffic.
Some recent analysis on such data can be seen at, e.g., [7, 45, 87].

Due to the high-dimensional spatiotemporal nature of the dataset, a sufficient dimension reduction or state aggregation is often a crucial first step to study a metropolitan-wide traffic pattern. To this end, we apply the high-order Markov model as described in Section <ref>. Specifically, we discretize the Manhattan region into a grid of $p=119$ states that forms a state space. Then, we collect all travel records in Manhattan of each driver from the dataset, sort them by time, and form into Markovian transition trajectories. In particular, each travel record is treated as a transition from the pickup to the drop-off location. If the drop-off location $i$ of the previous trip is different from the pickup location $j$ of the next trip by the same driver, we also form a transition from states $i$ to $j$.
Based on the trajectories, we can construct a high-order Markov chain with an order $d$ empirical transition probability tensor $\widehat{\bfc{P}}^{\mathrm{emp}}\in\bbR^{\otimes_{k=1}^d p}$ as described in Section <ref>. Assuming the true probability tensor is state aggregatable (Definition <ref>), we apply one-step TTOI  proposed in Section <ref> and obtain $\widehat{\bfc{P}}$. It is noteworthy if $d=2$, the described procedure of $\widehat{\bfc{P}}$ is equivalent to the classic matrix spectral decomposition in the literature. Figure <ref> plots the singular values of the sequential unfolding matrices of $\widehat{\bfc{P}}^{\rm{emp}}$ for $d=3$, which clearly demonstrates the low-TT-rankness of the probability transition tensor $\bfc{P}$. In the following experiments, we focus on the order-2 Markov model and analyze all consecutive two transitions: $i\to j \to k$, corresponding to the $d=3$ case.

Singular values of sequential unfolding matrices $[\widehat{\bfc{P}}^{\rm{emp}}]_1$ (left panel) and $[\widehat{\bfc{P}}^{\rm{emp}}]_2$ (right panel)

Inspired by the classic methods of matrix spectral decomposition, we aggregate all location states in Manhattan into a few clusters via both $\widehat{\bfc{P}}$ and $\widehat{\bfc{P}}^{\rm emp}$. Specifically, we calculate $\widehat{G}_d^{\top}$, i.e., the last TT-core of $\widehat{\bfc{P}}$, and $[\widehat{\bfc{P}}^{\mathrm{emp}}]_{d-1}$, i.e., the matricization of $\widehat{\bfc{P}}^{\mathrm{emp}}$ whose columns correspond to the last mode. Then we perform $k$-means on all columns of $\widehat{G}_d^{\top}$ and $[\widehat{\bfc{P}}^{\mathrm{emp}}]_{d-1}$, record the cluster index, associate the index to each location state, and plot the results in Figure <ref> (Panels (a)(b) are for TTOI and Panels (c)(d) are for empirical estimate). From Figure <ref> (a)(b), we can clearly identify four regions: (i) lower Manhattan (orange), (ii) midtown (dark blue),  (iii) upper west side (green), and (iv) upper east side (brown or  black). In contrast, direct clustering on $\bfc{P}^{\rm emp}$ yields less interpretable results as the majority points go to one cluster. It is also worth noting even the location information is not provided to this experiment, the resulting clusters in Figures <ref> (a)(b) show good spatial proximity between locations. This illustrates the effectiveness of TTOI in dimension reduction and state-aggregation for high-order Markov processes.

[$\widehat{G}_d$, $r=6,k=6$]

[$\widehat{G}_d$, $r=7,k=7$]

[$[\widehat{\bfc{P}}^{\rm{emp}}]_{d-1}$, $k=6$]

[$[\widehat{\bfc{P}}^{\rm{emp}}]_{d-1}$, $k=7$]

State aggregation based on TTOI and empirical estimate

Next, we illustrate the high-order nature of the city-wide taxi trip through the following experiment. For each initial state $i \in [p]$, we apply $k$-means to cluster the column span of $\widehat{\bfc{P}}_{[i, :, :]}$, where $\widehat{\bfc{P}}$ is the outcome of TTOI. We present the results in Figure <ref>, where the red triangles denote the given first state $i$ and $r=k=7$. If the city-wide taxi trips do not have significant high-order effects, $\widehat{\bfc{P}}$ should be reducible to a first order Markov process and $\widehat{\bfc{P}}_{[i,:,:]}$ should have similar values for different $i$. However, as we can see from Figure <ref> that the clustering results highly depends on the first state $i$, the high-order effects exist in the city-wide taxi trip Markov process. In addition, the states in different directions of $i$ are often clustered to different regions, which shows that the taxi drivers may tend to move to the same direction in consecutive trips, which yields the high-order effects in the driving trajectories.

Based on second order Markov model, state aggregation results are different with different initial state (the red triangle denotes the initial state $i$ in each subfigure)

§ DISCUSSIONS AND ADDITIONAL APPLICATIONS

In this paper, we propose a general framework for high-order SVD. We introduce a novel procedure, tensor-train orthogonal iteration (TTOI), that efficiently estimates the low tensor train rank structure from the high-order tensor observation. TTOI has significant advantages over the classic ones in the literature. We establish a general deterministic error bound for TTOI with the support of several new representation lemmas for tensor matricizations. Under the commonly studied spiked tensor model, we establish an upper bound for TTOI and a matching information-theoretic lower bound. We also illustrate the merits of TTOI through simulation studies and a real data example in New York City taxi trips.

In addition to the high-order Markov processes, the proposed TTOI can also be applied to the Markov random field (MRF) estimation. We give a brief description of MRF below. Consider an undirected graph $G = (V, E)$, where $V = \{1, \dots, d\}$ is a set of vertices and $E \subseteq V \times V$ is a collection of edges. Each vertex $i\in V$ is associated with a random variable $X_i$, taking values in $\{s_1,\dots,s_p\}$. In an MRF model, the distribution of $\{X_1,\dots, X_d\}$ can be factorized as

\begin{equation*}
\bbP(X_1,\dots,X_d)=\frac{1}{Z}\prod_{C\in \mathcal{C}}\psi_C(X_C),
\end{equation*}

where $\mathcal{C}$ is a collection of subgraphs of $G$ and $X_C=(X_v, v\in C)$ denotes the random vector corresponding to vertices in $C$. The joint probability function $\bbP(\cdot)$ can be written as a tensor $\bfc{P}\in\bbR^{\otimes_{k=1}^d p}$, where $\bfc{P}_{i_1,\dots,i_d}=\bbP(X_1=s_{i_1},\dots,X_d=s_{i_d})$. The MRFs have a wide range of applications, including image analysis [43, 91], genomic study [83], and natural language processing [18]. The readers are referred to, e.g., [82] for an introduction to MRFs.

A central problem of MRF is how to estimate the population density $\bfc{P}$ based on a limited number of samples $\{X^{(i)}_1,\dots,X^{(i)}_d\}_{i=1}^n$. It is straightforward to estimate $\bfc{P}$ via the empirical probability tensor $\widehat{\bfc{P}}^{\mathrm{emp}}$:

$$\widehat{\bfc{P}}^{\mathrm{emp}}_{i_1,\dots,i_d}=\sum_{i=1}^n\prod_{k =1}^d1(X^{(i)}_k=s_{i_k})\Big/n.$$

We can show that $\widehat{\bfc{P}}^{\mathrm{emp}}$ is unbiased for $\bfc{P}$. Recently, [57] pointed out that $\bfc{P}$ is often approximately low tensor-train rank in practice. To further exploit such the structure, we can conduct TTOI on $\widehat{\bfc{P}}^{\mathrm{emp}}$. Under regularity conditions, it can be shown that the entries of $\bfc{Z}$ are bounded and weakly independent, then Corollary <ref> suggests the following estimation error rate of the TTOI estimator: $\|\widehat{\bfc{P}}-\bfc{P}\|_F^2 \leq C\sum_{i=1}^{d}r_ir_{i-1}/(np^{2d-1})$, which can be significantly smaller than the estimation error of original empirical estimator $\widehat{\bfc{P}}^{\rm emp}$.

Moreover, the proposed framework can be also applied to high-order Markov decision process (high-order MDP). MDP has been commonly used as a baseline in control theory and reinforcement learning [27, 65, 73, 77]. Despite the wide applications of MDPs, most of the existing work focus on the first-order Markov processes. However, the high-order effects often appear, i.e., the transition probability at the current time depends not only on current, but also the past $(d-1)$ states and actions. See Figure <ref> for an example. Since the number of free parameters in such MDPs can be huge, a sufficient dimension reduction for the state and action space can be a crucial first step.
Similarly to the example of high-order Markov process in Section <ref>, the TTOI can be applied to achieve better dimension reduction and state aggregation for the high-order Markov decision processes.

Illustration of a high-order state aggregatable Markov decision process

§ PROOFS

We collect all technical proofs of this paper in this section.

§.§ Proof of Theorem <ref>

For convenience, let $\widehat{U}_i$, $\widehat{V}_i$, $R_i$ and $\widetilde{R}_i$ denote $\widehat{U}_i^{(0)}$, $\widehat{V}_i^{(1)}$, $R_i^{(0)}$ and $\widetilde{R}_i^{(0)}$, respectively. By Lemma <ref> and

\begin{equation*}
\begin{split}
&I_{p_2\cdots p_d} - P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}}\\ =& P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2\perp}} + P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_4 \otimes I_{p_2p_3})(\widehat{V}_{3\perp} \otimes I_{p_2})} + \cdots + P_{\widehat{V}_{d\perp} \otimes I_{p_2\dots p_{d-1}}},
\end{split}
\end{equation*}

we have

\begin{equation}\label{ineq14}
\begin{split}
&\left\|\widehat{\bfc{X}}^{(1)} - \bfc{X}\right\|_{\F}^2\\ =& \left\|\left[[\bfc{Y}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}\right]\widehat{V}_{2}^{\top}(\widehat{V}_{3}^{\top} \otimes I_{p_2})\cdots(\widehat{V}_d^{\top} \otimes I_{p_2\dots p_{d-1}}) - [\bfc{X}]_1\right\|_{\F}^2\\
=& \left\|[\bfc{Z}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}} + [\bfc{X}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}} - [\bfc{X}]_1\right\|_{\F}^2\\
\leq& C\bigg(\left\|[\bfc{Z}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}}\right\|_{\F}^2 + \left\|[\bfc{X}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2\perp}}\right\|_{\F}^2\\
&+ \left\|[\bfc{X}]_1P_{(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_4 \otimes I_{p_2p_3})(\widehat{V}_{3\perp} \otimes I_{p_2})}\right\|_{\F}^2 + \cdots + \left\|[\bfc{X}]_1P_{\widehat{V}_{d\perp} \otimes I_{p_2\dots p_{d-1}}}\right\|_{\F}^2\bigg)\\
\leq&  C\bigg(\left\|[\bfc{Z}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2}\right\|_{\F}^2 + \left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2\perp}\right\|_{\F}^2\\
&+ \left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_4 \otimes I_{p_2p_3}) (\widehat{V}_{3\perp} \otimes I_{p_2})\right\|_{\F}^2 +\cdots+ \left\|[\bfc{X}]_1(\widehat{V}_{d\perp} \otimes I_{p_2\dots p_{d-1}})\right\|_{\F}^2\bigg).
\end{split}
\end{equation}

To prove (<ref>), we only need to show that for all $2 \leq k \leq d$,

\begin{equation}\label{ineq13}
\begin{split}
&\left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}})\right\|_{\F}\\ \leq& C\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{Z}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\right\|_{\F},
\end{split}
\end{equation}

\begin{equation*}
[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}}) = [\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3} \otimes I_{p_2})\widehat{V}_{2\perp}
\end{equation*}

if $k = 2$ and

\begin{equation*}
[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}}) = [\bfc{X}]_1(\widehat{V}_{d\perp} \otimes I_{p_2\dots p_{d-1}})
\end{equation*}

if $k = d$.

By Lemma <ref>, we have

\begin{equation}\label{eq3}
\begin{split}
&\left\|[\bfc{X}]_1(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}})\right\|_{\F}\\
=& \left\|[A^{(p_2\cdots p_{k-1}, p_1)}]^\top([\bfc{X}]_{k-1} \otimes I_{p_2\cdots p_{k-1}})(\widehat{V}_d \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_2\cdots p_k})(\widehat{V}_{k\perp} \otimes I_{p_2\cdots p_{k-1}})\right\|_{\F}\\
=& \left\|[A^{(p_2\cdots p_{k-1}, p_1)}]^\top\left(\left([\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right)\otimes I_{p_2\cdots p_{k-1}}\right) \right\|_{\F}\\
=& \left\|[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}.
\end{split}
\end{equation}

The third equation holds since the realignment doesn't change the Frobenious norm.

Moreover, recall that $U_1 \in \bbR^{p_1 \times r_1}$ is the left singular space of $[\bfc{X}]_1$, and $\widetilde{U}_j \in \bbR^{p_jr_{j-1} \times r_j}$ is the left singular space of $(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j$ for $2 \leq j \leq d - 1$, by Lemma <ref>, for any $2 \leq k \leq d-1$,

\begin{equation}\label{eq5}
\begin{split}
[\bfc{X}]_k =& (I_{p_2\cdots p_k} \otimes [\bfc{X}]_1)A^{(p_2\cdots p_k, p_{k+1}\cdots p_d)}\\ =& (I_{p_2\cdots p_k} \otimes P_{U_1}[\bfc{X}]_1)A^{(p_2\cdots p_k, p_{k+1}\cdots p_d)}\\ =& (I_{p_2\cdots p_k} \otimes P_{U_1})(I_{p_2\cdots p_k} \otimes [\bfc{X}]_1)A^{(p_2\cdots p_k, p_{k+1}\cdots p_d)}\\ =& (I_{p_2\cdots p_k} \otimes P_{U_1})[\bfc{X}]_k,
\end{split}
\end{equation}

and for any $2 \leq j < k$,

\begin{equation}\label{eq6}
\begin{split}
&(I_{p_{j}\cdots p_k} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}\cdots p_{k}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{k}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_k\\
=&(I_{p_{j}\cdots p_k} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}\cdots p_{k}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{k}} \otimes \widehat{U}_{1}^{\top})(I_{p_{j+1}\cdots p_k} \otimes [\bfc{X}]_j)A^{(p_{j+1}\cdots p_k, p_{k+1}\cdots p_d)}\\
=&\left(I_{p_{j+1}\cdots p_k} \otimes \left[(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right]\right)A^{(p_{j+1}\cdots p_k, p_{k+1}\cdots p_d)}\\
=& \left(I_{p_{j+1}\cdots p_k} \otimes \left[P_{\widetilde{U}_j}(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right]\right)A^{(p_{j+1}\cdots p_k, p_{k+1}\cdots p_d)}\\
=& (I_{p_{j+1}\cdots p_k} \otimes P_{\widetilde{U}_j})(I_{p_{j}\cdots p_k} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}\cdots p_{k}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{k}} \otimes \widehat{U}_{1}^{\top})(I_{p_{j+1}\cdots p_k} \otimes [\bfc{X}]_j)A^{(p_{j+1}\cdots p_k, p_{k+1}\cdots p_d)}\\
=&(I_{p_{j+1}\cdots p_k} \otimes P_{\widetilde{U}_j})(I_{p_{j}\cdots p_k} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}\cdots p_{k}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{k}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_k,
\end{split}
\end{equation}

where $A^{(i,j)}$ is defined in (<ref>) for any $i, j > 0$.

Therefore, by (<ref>),

\begin{equation}\label{ineq26}
\begin{split}
&\left\|[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
=& \left\|(I_{p_2\cdots p_{k-1}} \otimes P_{U_1})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
=& \left\|(I_{p_2\cdots p_{k-1}} \otimes U_1^\top)[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& \left\|(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})(I_{p_2\cdots p_{k-1}} \otimes U_1)(I_{p_2\cdots p_{k-1}} \otimes U_1^\top)[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\&\cdot s_{\min}^{-1}\left((I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})(I_{p_2\cdots p_{k-1}} \otimes U_1)\right)\\
=& \left\|(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\cdot s_{\min}^{-1}(\widehat{U}_1^{\top}U_1).
\end{split}
\end{equation}

The inequality holds since $\|B\|_{\F} \leq \|AB\|_{\F}\cdot  s_{\min}^{-1}(A)$ for any invertible matrix $A \in \bbR^{m_1 \times m_1}$ and $B \in \bbR^{m_1 \times m_2}$; in the last step, we used $(I_{p_2\cdots p_{k-1}} \otimes U_1)(I_{p_2\cdots p_{k-1}} \otimes U_1^\top)[\bfc{X}]_{k-1} = (I_{p_2\cdots p_{k-1}} \otimes P_{U_1})[\bfc{X}]_{k-1} = [\bfc{X}]_{k-1}$.
Similarly to (<ref>), by (<ref>), for $1 \leq j \leq k-2$,

\begin{equation}\label{ineq27}
\begin{split}
&\left\|(I_{p_{j+1}\cdots p_{k-1}} \otimes \widehat{U}_{j}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
=& \left\|(I_{p_{j+2}\cdots p_{k-1}} \otimes P_{\widetilde{U}_{j+1}})(I_{p_{j+1}\cdots p_{k-1}} \otimes \widehat{U}_{j}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
=& \left\|(I_{p_{j+2}\cdots p_{k-1}} \otimes \widetilde{U}_{j+1}^\top)(I_{p_{j+1}\cdots p_{k-1}} \otimes \widehat{U}_{j}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& \left\|(I_{p_{j+2}\cdots p_{k-1}} \otimes \widehat{U}_{j+1}^\top)(I_{p_{j+1}\cdots p_{k-1}} \otimes \widehat{U}_{j}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
& \cdot  s_{\min}^{-1}(\widehat{U}_{j+1}^\top\widetilde{U}_{j+1}).
\end{split}
\end{equation}

By (<ref>) and (<ref>),

\begin{equation}\label{ineq12}
\begin{split}
&\left\|[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& \left\|(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F} s_{\min}^{-1}(\widehat{U}_1^{\top}U_1)\\
\leq& \left\|(I_{p_3\cdots p_{k-1}} \otimes \widehat{U}_2^{\top})(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\&\cdot  s_{\min}^{-1}(U_1^\top \widehat{U}_1) s_{\min}^{-1}(\widetilde{U}_2^{\top}\widehat{U}_2)\\ \leq&\dots\\\leq&  \left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}
\\&\cdot  s_{\min}^{-1}(U_1^\top \widehat{U}_1) s_{\min}^{-1}(\widetilde{U}_2^{\top}\widehat{U}_2)\cdots  s_{\min}^{-1}(\widetilde{U}_{k-1}^{\top}\widehat{U}_{k-1})\\
\leq& \left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\&\cdot \left(\frac{1}{\sqrt{1 - c_0^2}}\right)^{k-1}\\
\leq& C\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}.
\end{split}
\end{equation}

By the definition of $\widehat{V}_{k} \in \bbR^{(p_kr_k) \times r_{k-1}}$ and Lemma <ref>, we know that $\widehat{V}_{k}$ is the right singular space of

\begin{equation*}
\begin{split}
&\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{Y}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\\
=& \widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\\
&+ \widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{Z}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k}),
\end{split}
\end{equation*}

Lemma <ref> shows that

\begin{equation}\label{ineq11}
\begin{split}
&\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{X}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\widehat{V}_{k\perp}\right\|_{\F}\\
\leq& 2\left\|\widehat{U}_{k-1}^{\top}(I_{p_{k-1}} \otimes \widehat{U}_{k-2}^{\top})\cdots(I_{p_2\cdots p_{k-1}} \otimes \widehat{U}_1^{\top})[\bfc{Z}]_{k-1}(\widehat{V}_d \otimes I_{p_k\dots p_{d-1}})\cdots (\widehat{V}_{k+1} \otimes I_{p_k})\right\|_{\F}.
\end{split}
\end{equation}

Combine (<ref>), (<ref>) and (<ref>) together, we know that (<ref>) holds for all $2 \leq k \leq d$, which has finished the proof of Theorem <ref>.

§.§ Proof of Theorem <ref>

For $i \geq 1$, by the definition of $\bfc{X}^{(2i)}$ and Lemma <ref>, we have

\begin{equation*}
\begin{split}
\left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i)}\right\|_{\F}^2 =& \left\|\left(I_{p_1\cdots p_{d-1}} - P_{(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)})\cdots (I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)})\widehat{U}_{d-1}^{(2i)}}\right)[\bfc{Y}]_{d-1}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_{d-1}\right\|_{\F}^2 - \left\|P_{(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)})\cdots (I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)})\widehat{U}_{d-1}^{(2i)}}[\bfc{Y}]_{d-1}\right\|_{\F}^2\\
=& \left\|\bfc{Y}\right\|_{\F}^2 - \left\|\widehat{\bfc{X}}^{(2i)}\right\|_{\F}^2.
\end{split}
\end{equation*}

Similarly, we have

\begin{equation*}
\left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i - 1)}\right\|_{\F}^2 = \left\|\bfc{Y}\right\|_{\F}^2 - \left\|\widehat{\bfc{X}}^{(2i - 1)}\right\|_{\F}^2.
\end{equation*}

In addition, we have

\begin{equation*}
\begin{split}
\left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i)}\right\|_{\F}^2
=& \left\|[\bfc{Y}]_{d-1}\right\|_{\F}^2 - \left\|P_{(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)})\cdots (I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)})\widehat{U}_{d-1}^{(2i)}}[\bfc{Y}]_{d-1}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_{d-1}\right\|_{\F}^2 - \left\|\widehat{U}_{d-1}^{(2i)\top}(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-1}^{(2i)\top}(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_d^{(2i-1)}\right\|_{\F}^2\\ &- \left\|\widehat{U}_{d-1}^{(2i)\top}(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_{d\perp}^{(2i-1)}\right\|_{\F}^2\\
\leq& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-1}^{(2i)\top}(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_d^{(2i-1)}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})(I_{p_{d-2}p_{d-1}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_d^{(2i-1)}\right\|_{\F}^2.
\end{split}
\end{equation*}

The last equation holds since $\widehat{U}_{d-1}^{(2i)}$ is the left singular space of $(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})(I_{p_{d-2}p_{d-1}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_d^{(2i-1)}$.

For any $B \in \bbR^{n \times r}$ and $1 \leq l \leq r$, we can check that the $l$-th columns of $A^{(m,n)}B$ and $(I_m \otimes B \otimes I_m)A^{(m,r)}$ are equal:

$$(A^{(m,n)}B)_{[:, l]} = \sum_{j=1}^{n}B_{j,l}\sum_{k=1}^{m}e_{(k-1)mn+(j-1)m+k}^{(m^2n)} = ((I_m \otimes B \otimes I_m)A^{(m,r)})_{[:, l]}$$

where $e_{(k-1)mn+(j-1)m+k}^{(m^2n)}$ is the $((k-1)mn+(j-1)m+k)$-th canonical basis of $\bbR^{m^2n}$ and $A^{(i,j)}$ is defined in (<ref>). Therefore,

\begin{equation*}
A^{(m,n)}B = (I_m \otimes B \otimes I_m)A^{(m,r)}.
\end{equation*}

By the last equation and Lemma <ref>, we have

\begin{equation*}
\begin{split}
&(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})(I_{p_{d-2}p_{d-1}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-1}\widehat{V}_d^{(2i-1)}\\
=& (I_{p_{d-1}} \otimes \widehat{U}_{d-2}^{(2i)\top})(I_{p_{d-2}p_{d-1}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-1}} \otimes \widehat{U}_1^{(2i)\top})(I_{p_{d-1}} \otimes [\bfc{Y}]_{d-2})A^{(p_{d-1}, p_d)}\widehat{V}_d^{(2i-1)}\\
=& \left(I_{p_{d-1}} \otimes \left(\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}\right)\right)\left(I_{p_{d-1}} \otimes (\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\right)A^{(p_{d-1}, r_{d-1})}\\
=& \left(I_{p_{d-1}} \otimes \left(\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\right)\right)A^{(p_{d-1}, r_{d-1})}\\
=& \mathrm{Reshape}\big(\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}}), r_{d-2}p_{d-1}, r_{d-1}\big).
\end{split}
\end{equation*}

Since the realignment does not change the Frobenius norm, we have

\begin{equation}\label{ineq32}
\left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i)}\right\|_{\F}^2 \leq \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\right\|_{\F}^2.
\end{equation}

By similar proof of (<ref>), we have

\begin{equation*}
\begin{split}
\left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i)}\right\|_{\F}^2 \leq& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\widehat{V}_{d-1}^{(2i-1)}\right\|_{\F}^2\\
& - \left\|\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\widehat{V}_{d-1\perp}^{(2i-1)}\right\|_{\F}^2\\
\leq& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|\widehat{U}_{d-2}^{(2i)\top}(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\widehat{V}_{d-1}^{(2i-1)}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|(I_{p_{d-2}} \otimes \widehat{U}_{d-3}^{(2i)\top})\cdots(I_{p_2\cdots p_{d-2}} \otimes \widehat{U}_1^{(2i)\top})[\bfc{Y}]_{d-2}(\widehat{V}_d^{(2i-1)} \otimes I_{p_{d-1}})\widehat{V}_{d-1}^{(2i-1)}\right\|_{\F}^2\\
\leq& \cdots \\
\leq& \left\|[\bfc{Y}]_1\right\|_{\F}^2 - \left\|[\bfc{Y}]_1(\widehat{V}_d^{(2i-1)} \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3}^{(2i-1)} \otimes I_{p_2})\widehat{V}_{2}^{(2i-1)}\right\|_{\F}^2\\
=& \left\|[\bfc{Y}]_1\left(I_{p_2\cdots p_d} - P_{(\widehat{V}_d^{(2i-1)} \otimes I_{p_2\dots p_{d-1}})\cdots (\widehat{V}_{3}^{(2i-1)} \otimes I_{p_2})\widehat{V}_{2}^{(2i-1)}}\right)\right\|_{\F}^2\\
= & \left\|\bfc{Y} - \widehat{\bfc{X}}^{(2i-1)}\right\|_{\F}^2.
\end{split}
\end{equation*}

Similarly, we can prove (<ref>) holds for $k = 2i, i \geq 0$.

§.§ Proof of Theorem <ref>

Without loss of generality, we assume $\sigma^2 = 1$. We still let $\widehat{U}_i$, $\widehat{V}_i$, $R_i$ and $\widetilde{R}_i$ denote $\widehat{U}_i^{(0)}$, $\widehat{V}_i^{(1)}$, $R_i^{(0)}$ and $\widetilde{R}_i^{(0)}$, respectively.

Lemma <ref> Part 4 immediately shows that (<ref>) holds with probability at least $1 - Ce^{-cp}$. Next, we show that with probability at least $1 - Ce^{-cp}$,

\begin{equation}\label{ineq:perturbation_initial}
\begin{split}
\left\|\sin\Theta(\widehat{U}_k, \widetilde{U}_k)\right\| \leq C\frac{\sqrt{\sum_{i=1}^{k-1}p_ir_{i-1}r_i} + \sqrt{p_kr_{k-1}} + \sqrt{p_{k+1}\cdots p_d}}{\lambda_k} \leq \frac{1}{2}, \quad \forall 1 \leq k \leq d-1.
\end{split}
\end{equation}

Recall that

\begin{equation*}
\widehat{U}_1 = \text{SVD}_{r_1}^L([\bfc{Y}]_1), \quad [\bfc{Y}]_1 = [\bfc{X}]_1 + [\bfc{Z}]_1,
\end{equation*}

where $[\bfc{X}]_1 \in \bbR^{p_1 \times p_{-1}}$ satisfying rank$([\bfc{X}]_1) = r_1$, $[\bfc{Z}]_1 \in \bbR^{p_1 \times p_{-1}}$, by Lemmas <ref> and <ref>, with probability $1 - Ce^{-cp}$, we have

\begin{equation*}
\|\widehat{U}_{1\perp}^\top [\bfc{X}]_1\| \leq 2\|[\bfc{Z}]_1\| \leq C(p_1^{1/2} + (p_2\cdots p_d)^{1/2}).
\end{equation*}

Therefore, with probability at least $1 - Ce^{-cp}$,

\begin{equation*}
\left\|\sin\Theta(\widehat{U}_1, U_1)\right\| \leq \frac{\left\|\widehat{U}_{1\perp}^\top U_1U_1^\top [\bfc{X}]_1\right\|}{s_{r_1}(U_1^\top [\bfc{X}]_1)} = \frac{\left\|\widehat{U}_{1\perp}^\top [\bfc{X}]_1\right\|}{s_{r_1}([\bfc{X}]_1)} \leq C\frac{\sqrt{p_1} + \sqrt{p_2\cdots p_d  }}{\lambda_1}.
\end{equation*}

For $2 \leq i \leq  j \leq d - 1$, by the definition of $\widetilde{U}_i$ and Lemma <ref>, we have

\begin{equation}\label{eq7}
\begin{split}
=& (I_{p_{2}\cdots p_j} \otimes [\bfc{X}]_1)A^{(p_{2}\cdots p_j, p_{j+1}\cdots p_d)}
= (I_{p_{2}\cdots p_j} \otimes (P_{U_1}[\bfc{X}]_1))A^{(p_{2}\cdots p_j, p_{j+1}\cdots p_d)}\\
=& (I_{p_{2}\cdots p_j} \otimes P_{U_1})(I_{p_{2}\cdots p_j} \otimes [\bfc{X}]_1)A^{(p_{2}\cdots p_j, p_{j+1}\cdots p_d)}
= (I_{p_{2}\cdots p_j} \otimes U_1)(I_{p_{2}\cdots p_j} \otimes U_1^\top)[\bfc{X}]_j
\end{split}
\end{equation}

\begin{equation}\label{eq2}
\begin{split}
&\left(I_{p_i\cdots p_j} \otimes \widehat{U}_{i-1}^{\top}\right)\cdots \left(I_{p_2\cdots p_{j}} \otimes \widehat{U}_{1}^{\top}\right)[\bfc{X}]_j\\
=& \left(I_{p_{i+1}\cdots p_j} \otimes (I_{p_i} \otimes \widehat{U}_{i-1}^{\top}) \right)\cdots \left(I_{p_{i+1}\cdots p_j} \otimes (I_{p_2\cdots p_i} \otimes \widehat{U}_{1}^{\top})\right)(I_{p_{i+1}\cdots p_j} \otimes [\bfc{X}]_i)A^{(p_{i+1}\cdots p_j, p_{j+1}\cdots p_d)}\\
=& \left(I_{p_{i+1}\cdots p_j} \otimes \left((I_{p_i} \otimes \widehat{U}_{i-1}^{\top})\cdots (I_{p_2\cdots p_i} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_i\right)\right)A^{(p_{i+1}\cdots p_j, p_{j+1}\cdots p_d)}\\
=& \left(I_{p_{i+1}\cdots p_j} \otimes \left(P_{\widetilde{U}_i}(I_{p_i} \otimes \widehat{U}_{i-1}^{\top})\cdots (I_{p_2\cdots p_i} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_i\right)\right)A^{(p_{i+1}\cdots p_j, p_{j+1}\cdots p_d)}\\
=& \left(I_{p_{i+1}\cdots p_j} \otimes P_{\widetilde{U}_i}\right)\left(I_{p_{i+1}\cdots p_j} \otimes \left((I_{p_i} \otimes \widehat{U}_{i-1}^{\top})\cdots (I_{p_2\cdots p_i} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_i\right)\right)A^{(p_{i+1}\cdots p_j, p_{j+1}\cdots p_d)}\\
=& \left(I_{p_{i+1}\cdots p_j} \otimes \widetilde{U}_i\right)\left(I_{p_{i+1}\cdots p_j} \otimes \widetilde{U}_i^\top\right)\left(I_{p_i\cdots p_j} \otimes \widehat{U}_{i-1}^{\top}\right)\cdots \left(I_{p_2\cdots p_{j}} \otimes \widehat{U}_{1}^{\top}\right)[\bfc{X}]_j,
\end{split}
\end{equation}

where $I_{p_{i+1}\cdots p_j} = 1$ if $i = j$.
$$L_k = \left\|\sin \Theta \left(\widetilde{U}_k, \widehat{U}_k\right)\right\|, \quad 2 \leq k \leq d - 1.$$

For $k = 2$, by (<ref>) and Lemma <ref>, with probability at least $1 - Ce^{-cp}$,

\begin{equation*}
\begin{split}
s_{r_2}\left((I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\right) \geq& s_{\min}\left((I_{p_2} \otimes \widehat{U}_{1}^{\top})(I_{p_2} \otimes U_{1})\right)s_{r_2}([\bfc{X}]_2)\\
=& s_{\min}(\widehat{U}_{1}^{\top}U_{1})\lambda_2\\
=& \sqrt{1 - \|\sin\Theta (\widehat{U}_1, U_1)\|^2}\lambda_2\\
\geq& \sqrt{\frac{3}{4}}\lambda_2.
\end{split}
\end{equation*}

Since $\widehat{U}_{2} = \SVD_{r_2}^L((I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_2)$, and $(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_2 = (I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2 + (I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_2$, by Lemma <ref> and Lemma <ref>, we know that with probability at least $1 - Ce^{-cpr}$,

\begin{equation*}
\begin{split}
\|\widehat{U}_{2\perp}^{\top}(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\| \leq& 2\|(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_2\| \leq C(\sqrt{p_2r_1} + (p_3\cdots p_d)^{1/2} + \sqrt{p_1r_1}).
\end{split}
\end{equation*}

Combine the two previous inequalities together and recall that $\widetilde{U}_2$ is the left singular space of $(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2$, we have

\begin{equation*}
\begin{split}
\left\|\sin \Theta \left(\widehat{U}_2, \widetilde{U}_2\right)\right\| \leq& \frac{\|\widehat{U}_{2\perp}^{\top}\widetilde{U}_2\widetilde{U}_2^{\top}(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\|}{s_{r_2}\left(\widetilde{U}_2^{\top}(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\right)}\\ =& \frac{\|\widehat{U}_{2\perp}^{\top}(I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\|}{s_{r_2}\left((I_{p_2} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_2\right)}\\ \leq& C\frac{\sqrt{p_1r_1} + \sqrt{p_2r_1} + (p_3\cdots p_d)^{1/2}}{\lambda_2}
\end{split}
\end{equation*}

with probability at least $1 - Ce^{-cp}$.
Assume that (<ref>) holds for $k \leq j - 1$ with probability $1 - Ce^{-cp}$. For $k = j$,
by Lemma <ref> and (<ref>), with probability at $1 - Ce^{-cp}$, we have

\begin{equation}\label{ineq4}
\begin{split}
&s_{r_j}\left((I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)\\
\geq& s_{\min}\left((I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_j} \otimes \widetilde{U}_{j-1})\right)s_{r_j}\left((I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)\\
=& s_{\min}\left(\widehat{U}_{j-1}^{\top} \widetilde{U}_{j-1}\right)s_{r_j}\left((I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)\\
\geq& s_{\min}\left(\widehat{U}_{j-1}^{\top} \widetilde{U}_{j-1}\right)s_{\min}\left((I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})(I_{p_{j-1}p_{j}} \otimes \widetilde{U}_{j-2})\right)\\&\cdot s_{r_j}\left((I_{p_{j-2}p_{j-1}p_{j}} \otimes \widehat{U}_{j-3}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)\\
\geq& \cdots\\
\geq& s_{\min}\left(\widehat{U}_{j-1}^{\top} \widetilde{U}_{j-1}\right)\cdots s_{\min}\left(\widehat{U}_{1}^{\top} \widetilde{U}_{1}\right)s_{r_j}([\bfc{X}]_j)\\
=& \sqrt{1 - L_{j-1}^{2}}\cdots \sqrt{1 - L_1^{2}}\lambda_j\\
\geq& (\sqrt{3/4})^{j-1}\lambda_j
\geq c\lambda_j.
\end{split}
\end{equation}

In the last inequality, we used the fact that $d$ is a fixed number and $(\sqrt{3/4})^{j-1} \geq (\sqrt{3/4})^{d-1} \geq c$.

By the definition of $\widehat{U}_{j}$ and Lemma <ref>, we have
$$\widehat{U}_{j} = \SVD_{r_j}^L\left((I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_j\right).$$
Note that

\begin{equation*}
\begin{split}
&(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_j\\= &(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\\ &+ (I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_j,
\end{split}
\end{equation*}

by Lemma <ref>, with probability at least $1 - e^{-cpr^2}$,

\begin{equation*}
\begin{split}
&\left\|\widehat{U}_{j\perp}^{\top}(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right\|\\
\leq& 2\left\|(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_j\right\|\\
\leq& C\left(\left({\sum_{i=1}^{j-1}p_ir_{i-1}r_i}\right)^{1/2} + (p_{j}r_{j-1})^{1/2} +  (p_{j+1}\cdots p_d)^{1/2}\right).
\end{split}
\end{equation*}

Therefore, with probability at least $1 - Ce^{-cp}$,

\begin{equation*}
\begin{split}
\left\|\sin\Theta \left(\widehat{U}_j, \widetilde{U}_j\right)\right\|
\leq& \frac{\left\|\widehat{U}_{j\perp}^{\top}\widetilde{U}_j\widetilde{U}_j^{\top}(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right\|}{s_{r_j}\left(\widetilde{U}_j^{\top}(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)}\\
= & \frac{\left\|\widehat{U}_{j\perp}^{\top}(I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right\|}{s_{r_j}\left((I_{p_{j}} \otimes \widehat{U}_{j-1}^{\top})(I_{p_{j-1}p_{j}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}p_{j}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_j\right)}\\
\leq& C\frac{\left({\sum_{i=1}^{j-1}p_ir_{i-1}r_i}\right)^{1/2}  + (p_{j}r_{j-1})^{1/2} + (p_{j+1}\cdots p_d)^{1/2}}{\lambda_j}.
\end{split}
\end{equation*}

Therefore, (<ref>) holds with probability $1 - Ce^{-cp}$.

Finally, we consider (<ref>). Let $\mathcal{E}_0 = \{\eqref{ineq36} \text{ and }\eqref{ineq35} \text{ hold}\}$. Without loss of generality, we only show that under $\mathcal{E}_0$,

\begin{equation}\label{ineq9}
\left\|\sin\Theta\left(\widehat{V}_k, \widetilde{V}_{k}\right)\right\| \leq C\frac{\sqrt{\sum_{i=1}^{d}p_ir_{i-1}r_i}}{\lambda_{k-1}} \leq \frac{1}{2}, \quad \forall 2 \leq k \leq d.
\end{equation}

In fact, (<ref>) can be proved by induction. Let $V_d \in \bbR^{p_d \times r_{d-1}}$ be the right singular space of $[\bfc{X}]_{d-1}$. Then there exists an orthogonal matrix $\widetilde{Q}_{d-1} \in \mathbb{O}_{r_{d-1}}$ such that
$$V_d\widetilde{Q}_{d-1} = \text{SVD}^R\big(\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{X}]_{d-1}\big).$$
Similarly to (<ref>), under $\mathcal{E}_0$,

\begin{equation*}
s_{r_{d-1}}\left(\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{X}]_{d-1}\right) \geq \left(\sqrt{3/4}\right)^{d-1}\lambda_{d-1} \geq c\lambda_{d-1}.
\end{equation*}

Therefore, by Lemma <ref>, under $\mathcal{E}_0$,

\begin{equation*}
\begin{split}
\left\|\sin \Theta\left(\widehat{V}_d, V_d\right)\right\| =& \left\|\sin \Theta\left(\widehat{V}_d, V_d\widetilde{Q}_{d-1}\right)\right\|\\
\leq& \frac{\left\|\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{X}]_{d-1}\widehat{V}_{d\perp}^{\top}\right\|}{s_{r_{d-1}}\left(\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{X}]_{d-1}\right)}\\
\leq& \frac{2\left\|\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{Z}]_{d-1}\right\|}{s_{r_{d-1}}\left(\widehat{U}_{d-1}^\top(I_{p_{d-1}} \otimes \widehat{U}_{d-2}^\top)\cdots(I_{p_{d-1}\dots p_2} \otimes \widehat{U}_{1}^\top)[\bfc{X}]_{d-1}\right)}\\
\leq& C\frac{\sqrt{\sum_{i=1}^{d}p_ir_{i-1}r_i}}{\lambda_{d-1}}.
\end{split}
\end{equation*}

Suppose (<ref>) holds for $j+1 \leq k \leq d$. For $k = j$, since $\widetilde{V}_{j}$ is the right singular space of $[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})$, there exists $\widetilde{Q}_{j-1} \in \mathbb{O}_{r_{j-1}}$ such that
$$\widetilde{V}_{j}\widetilde{Q}_{j-1} = SVD^R\big(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\big).$$
By Lemma <ref>, (<ref>), (<ref>) and (<ref>), under $\mathcal{E}_0$,

\begin{equation*}
\begin{split}
&s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\right)\\
\geq& s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 2} \otimes I_{p_jp_{j+1}})(\widetilde{V}_{j + 1} \otimes I_{p_j})\right)\\&\cdot s_{\min}\left((\widetilde{V}_{j + 1}^\top \otimes I_{p_j})(\widehat{V}_{j + 1} \otimes I_{p_j})\right)\\
=& s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 2} \otimes I_{p_jp_{j+1}})\right)\\&\cdot s_{\min}(\widetilde{V}_{j + 1}^\top \widehat{V}_{j + 1})\\
\geq& \cdots\\
\geq& s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}\right)s_{\min}(\widetilde{V}_{d}^\top \widehat{V}_{d})\cdots s_{\min}(\widetilde{V}_{j + 1}^\top \widehat{V}_{j + 1})\\
\geq& s_{\min}(\widehat{U}_{j-1}^{\top}\widetilde{U}_{j-1})s_{r_{j-1}}\left((I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}\right)s_{\min}(\widetilde{V}_{d}^\top \widehat{V}_{d})\cdots s_{\min}(\widetilde{V}_{j + 1}^\top \widehat{V}_{j + 1})\\
\geq& \left(\sqrt{\frac{3}{4}}\right)^{j-1}\lambda_{j-1}\cdot\left(\sqrt{\frac{3}{4}}\right)^{d-j} \geq c\lambda_{j-1}.
\end{split}
\end{equation*}

Note that $\widehat{V}_j \in \mathbb{O}_{p_jr_j, r_{j-1}}$ is the right singular space of $\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})$ and

\begin{equation*}
\begin{split}
&\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{Y}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\\=&\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\\
&+\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j}),
\end{split}
\end{equation*}

By Lemma <ref>, under $\mathcal{E}_0$,

\begin{equation*}
\begin{split}
&\left\|\sin \Theta\left(\widehat{V}_j, \widetilde{V}_j\right)\right\| = \left\|\sin \Theta\left(\widehat{V}_j, \widetilde{V}_{j}\widetilde{Q}_{j-1}\right)\right\|\\
\leq& \frac{\left\|\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\widehat{V}_{j\perp}\right\|}{s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\right)}\\
\leq& \frac{2\left\|\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{Z}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\right\|}{s_{r_{j-1}}\left(\widehat{U}_{j-1}^{\top}(I_{p_{j-1}} \otimes \widehat{U}_{j-2}^{\top})\cdots(I_{p_2\cdots p_{j-1}} \otimes \widehat{U}_{1}^{\top})[\bfc{X}]_{j-1}(\widehat{V}_d \otimes I_{p_j\dots p_{d-1}})\cdots (\widehat{V}_{j + 1} \otimes I_{p_j})\right)}\\
\leq& C\frac{\left(\sum_{i=1}^d p_ir_ir_{i-1}\right)^{1/2}}{\lambda_{j-1}}.
\end{split}
\end{equation*}

Therefore, under $\mathcal{E}_0$, (<ref>) holds.

Thus, we have finished the proof of Theorem <ref>.

§.§ Proof of Corollary <ref>

Let $Q = \{\eqref{ineq35}, \eqref{ineq:perturbation_initial}\text{ hold}\}$, then $\bbP(Q^c) \leq C\exp(-cp)$ and

$$\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\|_{\F}^2 \leq C\sum_{i=1}^{d}p_ir_ir_{i-1} \quad \text{under } Q.$$

Under $Q^c$, due to the property of projection matrices, we know that

\begin{equation*}
\left\|\widehat{\bfc{X}}^{(t)}\right\|_{\F} \leq \|\bfc{Y}\|_{\F} \leq \|\bfc{X}\|_{\F} + \|\bfc{Z}\|_{\F}.
\end{equation*}

\begin{equation*}
\begin{split}
&\bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^4 \leq C\left(\bbE\left\|\widehat{\bfc{X}}^{(t)}\right\|_{\F}^4 + \|\bfc{X}\|_{\F}^4\right) \leq C\|\bfc{X}\|_{\F}^4 + C\bbE\|\bfc{Z}\|_{\F}^4\\
\leq& C\exp(4c_0p) + C\bbE\left(\chi_{p_1\cdots p_d}^2\right)^2 \leq C\exp(4c_0p) + C(p_1\cdots p_d)^2\\
\leq& C\exp(4c_0p) + C\exp\left(2c_0p\right) \leq C\exp(4c_0p).
\end{split}
\end{equation*}

Therefore, we have the following upper bound for the Frobenius norm risk of $\widehat{\bfc{X}}$:

\begin{equation*}
\begin{split}
&\bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^2 = \bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^21_Q + \bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^21_{Q^c}\\ \leq& C\sum_{i=1}^{d}p_ir_ir_{i-1} + \sqrt{\bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^4\cdot \bbP(Q^c)}\\
\leq& C\sum_{i=1}^{d}p_ir_ir_{i-1} + C\exp\left((4c_0 - c)p/2\right).
\end{split}
\end{equation*}

By selecting $c_0 < c/4$, we have

\begin{equation*}
\bbE\left\|\widehat{\bfc{X}}^{(t)} - \bfc{X}\right\|_{\F}^2 \leq C\sum_{i=1}^{d}p_ir_ir_{i-1}.
\end{equation*}

Therefore, we have finished the proof of Corollary <ref>.

§.§ Proof of Theorem <ref>

Since the i.i.d. Gaussian distribution, $\bfc{Z}\sim N(0, \sigma^2)$, is a special case of $\mathcal{D}$ and

$$ \inf_{\widehat{\bfc{X}}}\sup_{\bfc{X} \in \mathcal{F}_{\bp, \br}(\bm\lambda), D \in \mathcal{D}}\bbE_{\bfc{Z} \sim D}\left\|\widehat{\bfc{X}} - \bfc{X}\right\|_{\F}^2 \geq \inf_{\widehat{\bfc{X}}}\sup_{\bfc{X} \in \mathcal{F}_{\bp, \br}(\bm\lambda), \bfc{Z}\overset{\text{i.i.d.}}{\sim}N(0, \sigma^2)}\bbE_{\bfc{Z} \sim D}\left\|\widehat{\bfc{X}} - \bfc{X}\right\|_{\F}^2,$$

we only need to focus on the setting that $\bfc{Z}\sim N(0, \sigma^2)$ while developing the lower bound result.

Without loss of generality, assume $\sigma^2 = 1$. Since $d$ is a fixed number, we only need to show that for any $1 \leq i \leq d$,

\begin{equation}\label{ineq:lower}
\inf_{\widehat{\bfc{X}}}\sup_{\bfc{X} \in \mathcal{F}_{\bp, \br}(\bm\lambda)}\bbE\left\|\widehat{\bfc{X}} - \bfc{X}\right\|_{\F}^2 \geq cp_ir_ir_{i-1}.
\end{equation}

Suppose $\bfc{X}$ can be written as (<ref>), $U_j \in \bbR^{(p_jr_{j-1}) \times r_j}$ and $V_j \in \bbR^{(p_jr_{j}) \times r_{j-1}}$ are reshaped from $\bfc{G}_j \in \bbR^{r_{j-1} \times p_j \times r_j}$, $G_1 = U_1$, $G_d = V_d$. For any $1 \leq i \leq d-1$, by Lemma <ref>, we have

\begin{equation}\label{TT_form}
[\bfc{X}]_i = \left(I_{p_2\cdots p_i} \otimes U_1\right)\cdots (I_{p_i} \otimes U_{i-1})U_iV_{i+1}^\top\left(V_{i+2}^\top \otimes I_{p_{i+1}}\right)\cdots \left(V_d^\top \otimes I_{p_{i+1}\cdots p_{d-1}}\right).
\end{equation}

For all $j \neq i, 1 \leq j \leq d-1$, let $U_j \stackrel{\text{i.i.d.}}{\sim} N(0,1), V_d \stackrel{\text{i.i.d.}}{\sim} N(0,1)$ and $U_1, \dots, U_{i-1}, U_{i+1}, \dots, U_{d-1}, V_d$ are all independent. By Lemma <ref>, for any $1 \leq j \leq d-1$, we have

\begin{equation*}
s_{r_j}\left(\left(I_{p_2\cdots p_j} \otimes U_1\right)\cdots (I_{p_j} \otimes U_{j-1})U_j\right) \geq  s_{\min}\left(I_{p_2\cdots p_j} \otimes U_1\right)\cdots  s_{\min}(U_j) =  s_{r_1}(U_1)\cdots  s_{r_j}(U_j).
\end{equation*}

\begin{equation*}
s_{r_j}\left(V_{j+1}^\top\left(V_{j+2}^\top \otimes I_{p_{j+1}}\right)\cdots \left(V_d^\top \otimes I_{p_{j+1}\cdots p_{d-1}}\right)\right) \geq  s_{r_j}(V_{j+1})\cdots  s_{r_{d-1}}(V_d).
\end{equation*}

Moreover, Lemma <ref> Part 1 tells us

\begin{equation}\label{ineq21}
\begin{split}
& s_{r_j}\left(\left(I_{p_2\cdots p_j} \otimes U_1\right)\cdots (I_{p_j} \otimes U_{j-1})U_jV_{j+1}^\top\left(V_{j+2}^\top \otimes I_{p_{j+1}}\right)\cdots \left(V_d^\top \otimes I_{p_{j+1}\cdots p_{d-1}}\right)\right)\\
\geq&  s_{r_j}\left(\left(I_{p_2\cdots p_j} \otimes U_1\right)\cdots (I_{p_j} \otimes U_{j-1})U_j\right) s_{r_j}\left(V_{j+1}^\top\left(V_{j+2}^\top \otimes I_{p_{j+1}}\right)\cdots \left(V_d^\top \otimes I_{p_{j+1}\cdots p_{d-1}}\right)\right)\\
\geq&  s_{r_1}(U_1)\cdots  s_{r_j}(U_j) s_{r_j}(V_{j+1})\cdots  s_{r_{d-1}}(V_d).
\end{split}
\end{equation}

Recall that $V_j$ is reshaped from $U_j$ for all $1 \leq j \leq d - 1$, by [80][Corollary 5.35], we know that with probability at least $1 - Ce^{-cp}$, for all $1 \leq j \leq d-1, j \neq i$,

\begin{equation}\label{ineq16}
\begin{split}
&\frac{\sqrt{p_jr_{j-1}}}{4} \leq \sqrt{p_jr_{j-1}} - \sqrt{r_j} - \frac{\sqrt{p_jr_{j-1}}}{25} \leq  s_{r_j}(U_j) \leq s_1(U_j) \leq \sqrt{p_jr_{j-1}} + \sqrt{r_j} + \frac{\sqrt{p_jr_{j-1}}}{25} \leq 2\sqrt{p_jr_{j-1}},\\
&\frac{\sqrt{p_{j}r_{j}}}{4} \leq  s_{r_{j-1}}(V_j) \leq s_1(V_j) \leq 2\sqrt{p_{j}r_{j}}, \quad \text{ and } \quad \frac{\sqrt{p_d}}{4} \leq  s_{r_{d-1}}(V_d) \leq  s_{r_1}(V_d) \leq 2\sqrt{p_d}.
\end{split}
\end{equation}

For a fixed $U_0 \in \mathbb{O}_{p_ir_{i-1}, r_i}$, define the following ball with radius $\varepsilon > 0$,

\begin{equation*}
B(U_0, \varepsilon) = \left\{U' \in \mathbb{O}_{p_ir_{i-1}, r_i}: \|\sin\Theta(U', U_0)\|_{\F} \leq \varepsilon\right\}.
\end{equation*}

By Lemma 1 in [14], for $0 < \alpha < 1$ and $0 < \varepsilon \leq 1$, there exist $\widetilde{U}_i^{(1)'}, \dots, \widetilde{U}_i^{(m)'} \subseteq B(U_0, \varepsilon)$ such that

\begin{equation*}
m \geq \left(\frac{c_0}{\alpha}\right)^{r_i(p_ir_{i-1} - r_i)}, \quad \min_{1 \leq j \neq k \leq m}\left\|\sin\Theta\left(\widetilde{U}_i^{(j)'}, \widetilde{U}_i^{(k)'}\right)\right\|_{\F} \geq \alpha\varepsilon.
\end{equation*}

By Lemma 1 in [15], one can find a rotation matrix $O_k \in \mathbb{O}_{r_i}$  such that

\begin{equation*}
\|U_0 - \widetilde{U}_i^{(k)'}O_k\|_{\F} \leq \sqrt{2}\left\|\sin\Theta\left(U_0, \widetilde{U}_i^{(k)'}\right)\right\|_{\F} \leq \sqrt{2}\varepsilon.
\end{equation*}

Let $\widetilde{U}_i^{(k)} = \widetilde{U}_i^{(k)'}O_k$, we have

\begin{equation*}
\left\|\widetilde{U}_i^{(k)} - U_0\right\|_{\F} \leq \sqrt{2}\varepsilon, \quad \left\|\sin\Theta\left(\widetilde{U}_i^{(j)}, \widetilde{U}_i^{(k)}\right)\right\|_{\F} \geq \alpha\varepsilon, \quad 1 \leq j < k \leq m.
\end{equation*}

Let $U_i^{(k)} = S + \widetilde{U}_i^{(k)}$, where $S \stackrel{\text{i.i.d.}}{\sim} N(0, \tau^2)$.
Set $\tau \ge 8/\sqrt{p_i}$, [80][Corollary 5.35] shows that with probability at least $1 - Ce^{-cp}$,

\begin{equation}\label{ineq19}
\begin{split}
&\frac{\tau\sqrt{p_ir_{i-1}}}{8} \leq \tau\left(\sqrt{p_ir_{i-1}} - \sqrt{r_{i}} - \frac{\sqrt{p_ir_{i-1}}}{25}\right) - 1 \leq  s_{r_i}\left(S\right) -  s_{1}\left(\widetilde{U}_i^{(k)}\right) \leq  s_{r_i}\left(U_i^{(k)}\right)\\ \leq&  s_{1}\left(U_i^{(k)}\right) \leq  s_{1}\left(S\right) +  s_{1}\left(\widetilde{U}_i^{(k)}\right) \leq \tau\left(\sqrt{p_ir_{i-1}} + \sqrt{r_{i}} + \frac{\sqrt{p_ir_{i-1}}}{25}\right) + 1 \leq 2\tau\sqrt{p_ir_{i-1}}.
\end{split}
\end{equation}

If $2 \leq i \leq d - 1$, since $V_i^{(k)}$ is reshaped from $U_i^{(k)}$, we know that $V_i^{(k)} = T + \widetilde{V}_i^{(k)}$, where $T \stackrel{\text{i.i.d.}}{\sim} N(0, \tau^2)$, and $\widetilde{V}_i^{(k)}$ is realigned from $\widetilde{U}_i^{(k)}$. Notice that

\begin{equation*}
s_{1}(\widetilde{V}_i^{(k)}) = \|\widetilde{V}_i^{(k)}\| \leq \|\widetilde{V}_i^{(k)}\|_{\F} = \|\widetilde{U}_i^{(k)}\|_{\F} = r_{i},
\end{equation*}

Since $\tau \ge 8/\sqrt{p_i}$, by [80][Corollary 5.35], with probability at least $1 - Ce^{-cp_ir_i}$,

\begin{equation}\label{ineq20}
\begin{split}