The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
117,161 characters
CP-Factorization for High Dimensional Tensor Time Series and Double Projection Iterations
\if11
{
\spacingset{1.25}
\title{\bf \Large CP-Factorization for High Dimensional Tensor Time Series and Double Projection Iterations}
\author[1,2]{Jinyuan Chang}
\author[1]{Guanglin Huang}
\author[3]{Qiwei Yao}
\author[4]{Long Yu}
\affil[1]{\it \small Joint Laboratory of Data Science and Business
Intelligence, Institute of Statistical
Interdisciplinary Research, Southwestern University of Finance and Economics, Chengdu, China}
\affil[2]{\it \small State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China}
\affil[3]{\it \small Department of Statistics, The London School of Economics and Political Science, London, U.K.}
\affil[4]{\it \small School of Statistics and Data Science, Institute of Data Science and Statistics, Shanghai University of Finance and Economics, Shanghai, China}
\setcounter{Maxaffil}{0}
\date{\vspace{-5ex}}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{
\Large \bf CP-Factorization for High Dimensional Tensor Time Series and Double Projection Iterations
}
\end{center}
\medskip
} \fi
\spacingset{1.5}
\begin{abstract}
We adopt the canonical polyadic (CP) decomposition to model high-dimensional tensor time series. Our primary goal is to identify and estimate the factor loadings in the CP decomposition. We propose a one-pass estimation procedure through standard eigen-analysis for a matrix constructed based on the serial dependence structure of the data. The asymptotic properties of the proposed estimator are established under a general setting as long as the factor loading vectors are linearly independent, allowing the factors to be correlated and the factor loading vectors to be not nearly orthogonal. The procedure adapts to the sparsity of the factor loading vectors, accommodates weak factors, and demonstrates strong performance across a wide range of scenarios. To further reduce estimation errors, we also introduce an iterative algorithm based on a novel double projection approach. We theoretically justify the improved convergence rate of the iterative estimator, and derive the associated limiting distribution. A consistent estimator of the asymptotic variance is also provided, which plays a key role in the related inference problems. All results are validated through extensive simulations and two real data applications.
\end{abstract}
\noindent {\sl Keywords}: CP decomposition; dimension reduction; double projection iteration; statistical inference; tensor time series.
\spacingset{1.69}
\setlength{\abovedisplayskip}{0.2\baselineskip}
\setlength{\belowdisplayskip}{0.2\baselineskip}
\setlength{\abovedisplayshortskip}{0.2\baselineskip}
\setlength{\belowdisplayshortskip}{0.2\baselineskip}
\section{Introduction}\label{sec: introduction}
Due to recent advances in information technology and data science,
the demand for tensor analysis arises in a variety of fields, including but not limited to neuroimaging \citep{zhou2013tensor}, recommendation systems \citep{bi2018}, and dynamic transportation networks \citep{chen2021autoregressive,chen2022factor}. This motivates the rapid development of efficient inference tools and the associated theory for tensor analysis.
The size of tensor data is often large or extremely large. Low-rank approximation remains one of the most powerful tools for balancing between computational and statistical efficiencies.
In tensor analysis, there are two types of frequently used rank decompositions: the Tucker decomposition and the canonical polyadic (CP) decomposition, and both can be viewed as a natural extension of the singular value decomposition (SVD) for matrices \citep{kolda2009tensor}.
The Tucker decomposition is often achieved by the SVD on the unfolded matrices \citep{de2000multilinear}.
By contrast, computing the CP decomposition is
NP-hard, and the alternating least squares iteration remains as the workhorse method
\citep{wang2017tensor}.
In real applications, tensor data are often recorded in chronological order, and the dynamics of the data are often driven by a small number of factors.
Therefore, it is natural to extend the vector-valued factor models
\citep{bai2003inferential, lam2012factor, chang2015high} for tensor time series, leading to
two types of tensor factor models based on, respectively, the Tucker decomposition and the CP decomposition. Most existing studies focus on tensor Tucker-factor models, which can be traced back to \cite{wang2019factor}, where a two-way factor structure was introduced for matrix time series (i.e. a tensor with two modes).
See also \cite{elynn2020}, \cite{yu2022projected}, and \cite{chen2023statistical}.
Extensions to higher-order tensor Tucker-factor models have been studied in
\cite{chen2024rank}, \cite{han2024tensor},
\cite{chen2024semi}, \cite{barigozzi2023statistical}, and
\cite{HE2026105557}. Robust estimation methods for tensor Tucker-factor models are further considered in \cite{barigozzi2023robust} and \cite{barigozzi2025tail}.
Note that both the factors and the factor loadings in the Tucker decomposition are not uniquely defined: the decomposition is invariant under general invertible linear transformations. In empirical practice, rotations are often applied to the estimated loadings and factors to enhance interpretability.
In contrast, the factor loadings in tensor CP-factor models are uniquely defined up to the reflection and permutation indeterminacy, where the reflection indeterminacy is also referred to as the sign indeterminacy. This facilitates a straightforward and practically meaningful interpretation of the fitted models.
See, for example, the real data illustration in Section \ref{sec:application} below.
For tensor CP-factor models, \cite{han2024cp} propose an algorithm for the so-called High-Order Projection Estimator (HOPE), which consists of a principal component analysis (PCA) based initialization followed by recursive iterations. Remarkably, the estimation accuracy improves progressively with each iteration even though the initial estimate is not consistent. \cite{chen2026estimation} extend the HOPE using a contemporary covariance matrix and randomized projection, and also derive the limiting distribution of the associated estimator.
From the theoretical perspective, the HOPE is constructed under the following requirements:
(a) the factor loading vectors are nearly orthogonal, and
(b) the factors are almost uncorrelated.
Note that both factor
loadings and factors are uniquely defined in CP decomposition (up to the reflection and permutation indeterminacy), and there is no guarantee that those two requirements fulfill.
Free from requirement (b) and replacing ``nearly orthogonal" by ``linearly independent" in requirement (a),
\cite{chang2023modelling} propose a one-pass estimation procedure (i.e. without iterations) for matrix CP-factor models. Furthermore the ``linearly independent" requirement is freed in \cite{chang2024unified} which propose another one-pass estimation procedure applicable to matrix CP-factor models with rank-deficient factor loadings.
Note that the two estimation procedures are radically different: the method of \cite{chang2023modelling} is based on a generalized eigen-analysis while
the procedure of \cite{chang2024unified} is more involved and a key step is
to identify the factor loadings by a joint diagonalization of several symmetric matrices defined by the basis vectors of a linear system. However, neither \cite{chang2023modelling} nor \cite{chang2024unified} provide results on statistical inference, and it remains unclear how to generalize their procedures to higher-order tensor settings.
In this paper, we develop a unified framework for estimating CP-factor models for tensor time series, including the models for matrix time series as special cases. We do not impose requirements (a) and (b) stated above in our analysis. Under the mild assumption that the factor loading vectors are linearly independent, we propose two new estimation methods. Different from the method of \cite{chang2023modelling} which relies on the generalized eigen-analysis of certain matrices, our first method is established through the standard eigen-analysis of a well-designed matrix constructed based on the serial dependence structure of the
data, which can substantially attenuate the plug-in errors, and adapt to the higher (than two) mode tensor structures.
The second new estimation, termed as double projection method, is an iterative algorithm.
This algorithm substantially outperforms state-of-the-art iterative approaches \citep{han2024cp, chen2026estimation} in terms of both statistical and computational efficiency, particularly in the presence of correlated factors. The superiority of the iterative estimator is rigorously established through theoretical analysis and extensive numerical studies. Furthermore, we derive a tractable limiting representation for this new iterative estimator, along with its explicit asymptotic distribution. An estimator for the asymptotic variance is also provided, which makes the related statistical inference feasible. The \textsf{R}-function \texttt{CP\_TTS} for implementing our newly proposed methods is available publicly in the \texttt{HDTSA} package \citep{chang2024hdtsa}.
The rest of the paper is organized as follows. Section \ref{sec: model} introduces the tensor CP-factor models and the basic settings. Section \ref{sec: methodology} presents our methodology, including the one-pass and iterative estimators together with the inference procedure.
Sections \ref{sec:numerical} and \ref{sec:application} validate the performance of our proposed methods through simulation studies and a real data analysis, respectively. Section \ref{sec: theoretical} develops the theoretical guarantees of the proposed estimators. Section \ref{sec: discuss} provides some discussion.
All technical details, and additional simulation and empirical results, are provided in the supplementary material. The replication code for both the simulations and the real data analysis is
available at the GitHub repository: \url{https://github.com/JinyuanChang-Lab/CPTensorTimeSeries}.
\emph{Notation}. For any integer $p\ge 1$, let $[p]=\{1,\ldots,p\}$, and $\mathbf{I}_p$ be the $p\times p$ identity matrix. Denote by $I(\cdot)$ the indicator function. For a vector $\mathbf{a}=(a_1,\ldots,a_p)^\mathrm{\scriptscriptstyle \top }$, define $|\mathbf{a}|_2=(\sum_{i=1}^p a_i^2)^{1/2}$ and $|\mathbf{a}|_0 = \sum_{i=1}^p I(a_i \neq 0)$. For a complex vector $\mathbf{a}$, $\operatorname{Re}(\mathbf{a})$ denotes the vector of its real parts taken entry-wise. For a matrix $\mathbf{A}=(a_{i,j})_{p_1 \times p_2}$, denote by $\sigma_i(\mathbf{A})$, $\sigma_{\max}(\mathbf{A})$, and $\sigma_{\min}(\mathbf{A})$, respectively, its $i$-th largest, maximum, and minimum singular values. We write $\|\mathbf{A}\|_2=\sqrt{\sigma_1(\mathbf{A}^\mathrm{\scriptscriptstyle \top }\mathbf{A})}$, $\|\mathbf{A}\|_{\rm F}=\sqrt{\sum_{i=1}^{p_1}\sum_{j=1}^{p_2}a_{i,j}^2}$, and $ | \mathbf{A} |_{\max}=\max_{i\in[p_1],j\in[p_2]}|a_{i,j}|$. For an $m$-mode tensor $\mathcal{A}=(a_{i_1,\ldots,i_m})_{p_1\times\cdots \times p_m}$, we write $[\mathcal{A}]_{i_1,\ldots,i_m}=a_{i_1,\ldots,i_m}$. The vectorization $\textup{vec}(\mathcal{A})$ is an ${\textstyle{\prod}}_{j=1}^m p_j$-dimensional vector, with the $\{ 1 + \sum_{k = 1}^m(i_k - 1)\prod_{\ell=1}^{k-1}p_\ell\}$-th element being $a_{i_1,\ldots,i_m}$. The matricization $\textup{Mat}_q(\mathcal{A})$ is a $p_q \times \prod_{j \neq q}p_j$ matrix, with the $\{i_{q}, 1 + \sum_{k \neq q}^m(i_k - 1)\prod_{\ell \neq q}^{k-1}p_\ell\}$-th element being $a_{i_1,\ldots,i_m}$. For two sequences of positive numbers $\{a_n\}_{n\ge 1}$ and $\{b_n\}_{n\ge 1}$, we write $a_n\lesssim b_n$ or $b_n\gtrsim a_n$ if $\limsup_{n\rightarrow\infty}a_n/b_n< \infty$, $a_n\asymp b_n$ if and only if $b_n\gtrsim a_n$ and $a_n\gtrsim b_n$ hold simultaneously, and $a_n\ll b_n$ or $b_n\gg a_n$ if $\limsup_{n\rightarrow\infty}a_n/b_n=0$. For any $a, b\in\mathbb{R}$, let $a \vee b=\max(a,b)$ and $a \wedge b=\min(a,b)$. Denote by $\lfloor x \rfloor$ the largest integer less than or equal to $x$. Let $\circ$ and $\otimes$ denote the vector outer product and the Kronecker product, respectively.
\section{Model}\label{sec: model}
Let $\mathcal{Y}_t\in \mathbb{R}^{d_1\times \cdots \times d_m}$ be an observed $m$-mode tensor with
$m\geq2$. We consider the tensor CP-factor model \citep{han2024cp} as follows:
\begin{equation}\label{model cp}
\mathcal{Y}_t= \sum_{i=1}^rw_if_{t,i}\,\mathbf{a}_{i,1} \circ \mathbf{a}_{i,2} \circ \cdots \circ \mathbf{a}_{i,m}+\mathcal{E}_t\,,\quad t\geq1\,,
\end{equation}
where $1\leq r \leq \min_{j\in[m]}d_j$ is a fixed but unknown constant, $\mathcal{E}_t \in \mathbb{R}^{d_1\times \cdots \times d_m}$ is the idiosyncratic error tensor, $\mathbf{f}_t = (f_{t,1},\ldots,f_{t,r})^\mathrm{\scriptscriptstyle \top }$ is the $r$-dimensional factor vector, $w_i$ is the strength of the $i$-th factor, and $\mathbf{a}_{i,j}$ is a $d_j$-dimensional factor loading vector corresponding to the $i$-th factor and $j$-th mode.
Without loss of generality, we assume $|\mathbf{a}_{i,j}|_2=1$ for any $i\in[r]$ and $j\in[m]$. When $w_i=1$ for all $i\in[r]$, model \eqref{model cp} is an extension of the model considered in \cite{chang2023modelling} from matrix to tensor regimes. In practice, $w_i$ can be either a constant or grow with the dimensions. To understand this, let us consider a toy example that $[\mathcal{Y}_{t}]_{h_1,\ldots,h_m} = \beta f_{t}+[\mathcal{E}_t]_{h_1,\ldots,h_m}$ for any $h_j\in[d_j]$ with $j\in[m]$. This example can be formulated as model \eqref{model cp} with $r = 1$, $w_1=\beta(\prod_{j=1}^{m}d_j)^{1/2}$ and $\mathbf{a}_{1,j}=(d_j^{-1/2},\ldots,d_j^{-1/2})^\mathrm{\scriptscriptstyle \top }$ for each $j \in [m]$, where $w_1$ will diverge if at least one $d_j$ grows to infinity as $n\rightarrow\infty$.
The key problem of interest for the tensor CP-factor model \eqref{model cp} is to identify the factor loading vectors $\{\mathbf{a}_{i,j}\}_{i\in[r],j\in[m]}$ and also provide the suitable estimates for them. To do this, we first impose the following regularity assumption on the tensor CP-factor model \eqref{model cp}.
\begin{assumption}\label{error}
It holds that $\mathbb{E}(\mathcal{E}_t)={\bf 0}$ for any $t\in[n]$, $\mathbb{E}(\mathcal{E}_t\otimes\mathcal{E}_s)={\bf 0}$ for any $t\ne s$, and $\mathbb{E}(f_{t,i}\mathcal{E}_s)={\bf 0}$ for any $i\in[r]$ and $t,s \in [n]$.
\end{assumption}
\begin{remark}
\textup{(a)} Assumption \textup{\ref{error}} is significantly weaker than the assumptions imposed in \textup{\cite{han2024cp}}, which is satisfied automatically under the assumptions of \textup{\cite{han2024cp}}. More specifically, \textup{\cite{han2024cp}} require the error process $\{\mathcal{E}_t\}_{t \ge 1}$ to be independent Gaussian tensors conditional on the factor process $\{\mathbf{f}_{t}\}_{t \ge 1}$. Furthermore, \textup{\cite{han2024cp}} require the factor process $\{\mathbf{f}_{t}\}_{t \ge 1}$ to be stationary with zero mean and also to satisfy $\mathbb{E}(f_{t,i}^2) = 1$ and $\mathbb{E}(f_{t,i}f_{t-k,j})=0$ for all $i\neq j$ and $k\ge 1$, while the stationarity and zero mean are not necessarily required in our framework. \textup{(b)} For the relationship between the factor process $\{\mathbf{f}_t\}_{t \ge 1}$ and the error process $\{\mathcal{E}_t\}_{t \ge 1}$, the theoretical analysis of the one-pass estimator introduced in Section \textup{\ref{sec: initial}} only requires $\mathbb{E}(f_{t,i}\mathcal{E}_s)={\bf 0}$ for any $i\in[r]$ and $t \le s$, while the iterative estimator introduced in Section \textup{\ref{sec: Double projection iterations}} necessitates the stronger condition as stated in Assumption \textup{\ref{error}}. \textup{(c)} Write $\sigma_{t,i}^2 = \mathbb{E}(f^2_{t,i})$. Different from \textup{\cite{han2024cp}}, we do not require $\sigma_{t,i}^2$
equal to $1$, which allows $\sigma_{t,i}^2$ to vary with $t$ for each given $i$. If $\sigma_{t,i}^2 \equiv \sigma_{i}^2$ for all $t \in [n]$, we can assume $\mathbb{E}(f^2_{t,i}) = 1$ without loss of generality. From this perspective, our model setting is more general than that in \textup{\cite{han2024cp}}. For more general scenarios, $w_i$ and $f_{t,i}$ cannot be identified separately. Nevertheless, the loading vectors $\mathbf{a}_{i,1},\ldots,\mathbf{a}_{i,m}$ remain identifiable up to the reflection and permutation indeterminacy. \textup{(d)} The idiosyncratic errors are assumed to be serially uncorrelated, which enables a direct separation of the signal part and the noise part through the auto-covariances of the observed data. This is a common assumption in the literature on factor modeling; see, for example, \textup{\cite{lam2012factor}} and \textup{\cite{han2024cp}}. By contrast, the factors are allowed to be serially correlated; see Assumption \textup{\ref{mixing}} in Section \textup{\ref{sec: asmp}}.
\end{remark}
For each $j\in[m]$, write $d_{\scalebox{1}{\text{-}} j} = \prod_{j^\prime\ne j}^md_{j^\prime}$. Following the tensor matricization in \cite{kolda2009tensor}, we can reshape $\mathcal{Y}_t$ into a $d_j\times d_{\scalebox{1}{\text{-}} j}$ matrix as follows:
\begin{equation}\label{unfold X}
\mathbf{Y}_{t,j}=\text{Mat}_j(\mathcal{Y}_t)= \underbrace{ {\textstyle\sum\nolimits}_{i=1}^r w_if_{t,i}\mathbf{a}_{i,j}\mathbf{b}_{i,j}^\mathrm{\scriptscriptstyle \top }}_{\mathbf{C}_{t,j}} +\underbrace{\text{Mat}_j(\mathcal{E}_t)}_{\mathbf{E}_{t,j}}\,,~~~~j\in[m]\,,
\end{equation}
where $ \mathbf{b}_{i,j}=\mathbf{a}_{i,m}\otimes\cdots\otimes \mathbf{a}_{i,j+1}\otimes \mathbf{a}_{i,j-1}\otimes\cdots\otimes \mathbf{a}_{i,1} \in \mathbb{R}^{d_{\scalebox{1}{\text{-}} j}}$. Write $\mathbf{A}_j=(\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j})$, $\mathbf{B}_j=(\mathbf{b}_{1,j},\ldots,\mathbf{b}_{r,j})$, and $\mathbf{X}_t = \textup{diag}(w_1f_{t,1},\ldots,w_r f_{t,r})$. Then model \eqref{unfold X} can be rewritten as the matrix CP-factor model considered in \cite{chang2023modelling}:
\begin{equation}\label{eq:matrixform}
\mathbf{Y}_{t,j} = \mathbf{A}_j \mathbf{X}_{t} \mathbf{B}_j^{\mathrm{\scriptscriptstyle \top }} + \mathbf{E}_{t,j}\,, ~~~~j\in[m]\,.
\end{equation}
We assume $\textup{rank}(\mathbf{A}_j) = r$ for $j \in [m]$, which is also required in
\cite{han2024cp}.
Under this assumption, Lemma \ref{pro:rank-B} in the supplementary material shows that $\mathbf{B}_1,\ldots,\mathbf{B}_m$ also have full column rank.
We could estimate $\mathbf{a}_{i,j}$ by the method proposed in \cite{chang2023modelling}. However, for the tensors with more than two modes, the matrix $\mathbf{B}_j$ in \eqref{eq:matrixform} involves a Kronecker product structure. Direct application of the method in \cite{chang2023modelling} would discard this structural information, leading to substantial statistical inefficiency. Meanwhile, the two-stage estimation approach of \cite{chang2023modelling} introduces cross-step plug-in errors, which significantly complicates statistical inference for $\mathbf{a}_{i,j}$ in high-dimensional settings. In this paper, we propose a new one-pass estimation method that only requires the eigen-decomposition of a $d_j \times d_j$ matrix to obtain estimates of $\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j}$. This novel one-pass design eliminates cross-stage plug-in errors, and further motivates an iterative estimation procedure introduced in Section \ref{sec: Double projection iterations}. Moreover, simulation studies in Section \ref{sec:numerical} show that even for matrix-variate cases, our newly proposed methods perform significantly better than the procedure in \cite{chang2023modelling}.
\section{Methodology}\label{sec: methodology}
\subsection{One-pass estimation of the factor loading vectors}\label{sec: initial}
Let $\xi_{t}$ be a linear combination of the components of $\mathcal{Y}_t$. For any $k\ge 1$ and $t \ge k+1$, we define $\bm{\Xi}_{t,k,j}=\mathbb{E}[\{\mathbf{Y}_{t,j}-\mathbb{E}(\bar\mathbf{Y}_j)\}\{\xi_{t-k}-\mathbb{E}(\bar\xi)\}]$ with $\bar \mathbf{Y}_j=n^{-1}\sum_{t=1}^n\mathbf{Y}_{t,j}$ and $\bar\xi=n^{-1}\sum_{t=1}^n\xi_{t}$. Given observations $\{\mathcal{Y}_t\}_{t=1}^n$, for any $k\geq1$, we write
\begin{equation}\label{Xi and Sigma Y}
\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)=\frac{1}{n-k}\sum_{t=k+1}^n\bm{\Xi}_{t,k,j}\,,
\end{equation}
and let $\mathbf{G}_{k,\xi}=\textup{diag}(g_{k,1,\xi},\ldots,g_{k,r,\xi})$ be an $r\times r$ diagonal matrix with
\begin{equation}\label{gki}
g_{k,i,\xi}=\frac{1}{n-k}\sum_{t=k+1}^n w_i\mathbb{E}[\{f_{t,i}-\mathbb{E}(\bar f_i)\}\{\xi_{t-k}-\mathbb{E}(\bar\xi)\}]\,,
\end{equation}
where $\bar f_i=n^{-1}\sum_{t=1}^n f_{t,i}$. Then $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k) = \mathbf{A}_j \mathbf{G}_{k,\xi} \mathbf{B}_j^{\mathrm{\scriptscriptstyle \top }}$. For each given $j \in [m]$, by singular value decomposition of $\mathbf{B}_j$, there exist a $d_{\scalebox{1}{\text{-}} j} \times r$ column-orthogonal matrix $\mathbf{Q}_j$ and an $r\times r$ invertible matrix $\mathbf{V}_j$ such that $\mathbf{B}_j=\mathbf{Q}_j\mathbf{V}_j$. If $\textup{rank}(\mathbf{G}_{1,\xi}) = r = \textup{rank}(\mathbf{G}_{2,\xi})$, we define
\begin{align}\label{Kbj}
\mathbf{K}_{1,2,j}&=\bm{\Sigma}_{\mathbf{Y}_j,\xi}(1)\mathbf{Q}_j \{\mathbf{Q}_j^\mathrm{\scriptscriptstyle \top }\bm{\Sigma}_{\mathbf{Y}_j,\xi}(2)^\mathrm{\scriptscriptstyle \top }\bm{\Sigma}_{\mathbf{Y}_j,\xi}(2)\mathbf{Q}_j \}^{-1}\mathbf{Q}_j^\mathrm{\scriptscriptstyle \top }\bm{\Sigma}_{\mathbf{Y}_j,\xi}(2)^\mathrm{\scriptscriptstyle \top }\,.
\end{align}
Since $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k) = \mathbf{A}_j \mathbf{G}_{k,\xi} \mathbf{B}_j^{\mathrm{\scriptscriptstyle \top }}$, we have $\mathbf{K}_{1,2,j} = \mathbf{A}_j\mathbf{G}_{1,\xi}\mathbf{G}_{2,\xi}^{-1}(\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j)^{-1}\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }$, which can be used to identify $\mathbf{A}_j$. Write $\bar\lambda_i=g_{2,i,\xi}^{-1}g_{1,i,\xi}$ with $g_{k,i,\xi}$ defined in \eqref{gki}. For each given $j\in[m]$, since $(\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j)^{-1}\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j=\mathbf{I}_r$, then
$
\mathbf{K}_{1,2,j}\mathbf{a}_{i,j}=\bar\lambda_i\mathbf{a}_{i,j}$ for any $i\in[r]$,
which implies that, as long as $\bar\lambda_1,\ldots,\bar\lambda_r$ are distinct, $\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j}$ can be identified uniquely up to the reflection and permutation indeterminacy by solving the eigen-equation
\begin{equation}\label{generalized eigenequation}
\mathbf{K}_{1,2,j}\mathbf{a}=\lambda\mathbf{a}\,.
\end{equation}
In practice, $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ and $\mathbf{Q}_j$ in \eqref{Kbj} are unknown. Given observations $\{\mathcal{Y}_t\}_{t=1}^n$, in the spirit of \cite{bickel2008covariance}, we can estimate $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ by
\begin{equation}\label{hat Sigma kj}
\tilde\bm{\Sigma}_{k,j}=T_{\delta_1}\{\tilde\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)\}~~\textrm{with}~~\tilde\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)=\frac{1}{n-k}\sum_{t=k+1}^n(\mathbf{Y}_{t,j}-\bar\mathbf{Y}_j)(\xi_{t-k}-\bar\xi)\,,
\end{equation}
where $T_{\delta_1}(\cdot)$ is a thresholding operator, i.e., $[T_{\delta_1}(\mathbf{W})]_{i,j} = W_{i,j}\, I(|W_{i,j}| \ge \delta_1)$ for any matrix $\mathbf{W} = (W_{i,j})$ with the threshold level $\delta_1 \ge 0$. To estimate $\mathbf{Q}_j$, define
\begin{equation}\label{tilde Mj}
\mathbf{M}_{j}=\sum_{k=1}^K \bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)^\mathrm{\scriptscriptstyle \top }\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)
\end{equation}
for some predetermined integer $K\ge 1$. Under Assumption \ref{error} and $\textup{rank}(\mathbf{A}_{j}) = r = \textup{rank}(\mathbf{G}_{1,\xi})$, we have $\mathbf{M}_{j} = \mathbf{B}_j (\sum_{k=1}^K\mathbf{G}_{k,\xi}\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j\mathbf{G}_{k,\xi})\mathbf{B}_j^\mathrm{\scriptscriptstyle \top }$ with $\textup{rank}(\mathbf{M}_{j}) = r$, which implies the columns of $\mathbf{Q}_j$ are in the linear space spanned by the $r$ orthonormal eigenvectors of $\mathbf{M}_{j}$ corresponding to its $r$ largest eigenvalues. Let $\tilde{r}$ be a consistent estimate of $r$, which will be specified in Section \ref{sec:tuning}. Then we select $\tilde\mathbf{Q}_j$ as a $d_{\scalebox{1}{\text{-}} j} \times \tilde{r}$ matrix of which the columns are the $\tilde{r}$ orthonormal eigenvectors of $\tilde{\mathbf{M}}_{j} =\sum_{k=1}^K \tilde\bm{\Sigma}_{k,j}^\mathrm{\scriptscriptstyle \top }\tilde\bm{\Sigma}_{k,j}$ corresponding to its $\tilde{r}$ largest eigenvalues.
By plugging $\tilde\mathbf{Q}_j$ and $\tilde\bm{\Sigma}_{k,j}$ into \eqref{Kbj}, we can estimate $\mathbf{K}_{1,2,j}$ by
\begin{equation}\label{tilde K21j}
\tilde\mathbf{K}_{1,2,j}=\tilde\bm{\Sigma}_{1,j}\tilde\mathbf{Q}_j(\tilde\mathbf{Q}_j^\mathrm{\scriptscriptstyle \top }\tilde\bm{\Sigma}_{2,j}^\mathrm{\scriptscriptstyle \top }\tilde\bm{\Sigma}_{2,j}\tilde\mathbf{Q}_j)^{-1}\tilde\mathbf{Q}_j^\mathrm{\scriptscriptstyle \top }\tilde\bm{\Sigma}_{2,j}^\mathrm{\scriptscriptstyle \top }\,.
\end{equation}
Let $\tilde\lambda_{i,j}$ and $\tilde\mathbf{a}_{i,j}$ be the $i$-th largest eigenvalue (in terms of absolute value) and the associated eigenvector of $\tilde\mathbf{K}_{1,2,j}$, respectively. Note that $\mathbf{K}_{1,2,j}$ and $\tilde\mathbf{K}_{1,2,j}$ are in general nonsymmetric matrices. Although \eqref{generalized eigenequation} indicates that the eigenvectors of $\mathbf{K}_{1,2,j}$ are real vectors, those of the estimator $\tilde\mathbf{K}_{1,2,j}$ are not guaranteed to always be real vectors in practice. If $\tilde\mathbf{a}_{i,j}$ is a complex vector, we replace it by $\operatorname{Re}(\tilde\mathbf{a}_{i,j})/|\operatorname{Re}(\tilde\mathbf{a}_{i,j})|_2$, which has a negligible effect on the consistency of the estimator. The loading matrix $\mathbf{A}_j=(\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j})$ is then estimated by $(\tilde{\mathbf{a}}_{1,j},\ldots,\tilde{\mathbf{a}}_{\tilde{r},j})$. When $\tilde r=r$, for each given $j\in[m]$, Theorem \ref{thm: aij} in Section \ref{sec: iterative theorem} shows that the proposed one-pass estimator $\{\tilde\mathbf{a}_{i,j}\}_{i\in[\tilde r]}$ is consistent to $\{\mathbf{a}_{i,j}\}_{i\in[r]}$ up to the reflection and permutation indeterminacy.
\subsection{Double projection estimation for the factor loading vectors}\label{sec: Double projection iterations}
The one-pass estimation procedure entails thresholding the large $d_j\times d_{-j}$ matrix $\tilde\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ in \eqref{hat Sigma kj}, where the thresholding errors of all the elements in the matrix will accumulate. Moreover, its performance depends on the choice of the linear combination $\xi_t$. To address these issues, we introduce a novel double projection iteration method that achieves high accuracy without relying on the uncorrelated factor assumption or the near-orthogonality condition on factor loadings required in \cite{han2024cp}. More specifically, when a consistent initial estimator is available, by projecting the tensor-valued data into lower dimensions, it will suffice to perform thresholding on a $d_j$-dimensional vector rather than a large $d_j\times d_{-j}$ matrix. Meanwhile, the initial estimator can be leveraged to construct a specific linear combination of
$\mathcal{Y}_t$, denoted by $\tilde{\xi}_t$, to further reduce the estimation error.
Write $(\mathbf{b}_{1,j}^{\mathrm{\scriptscriptstyle + }},\ldots,\mathbf{b}_{r,j}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}=(\mathbf{B}_j^\mathrm{\scriptscriptstyle \top }\mathbf{B}_j)^{-1}\mathbf{B}_j^{\mathrm{\scriptscriptstyle \top }}$. When $\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}$ is given, it follows from \eqref{eq:matrixform} that, for each $j\in[m]$, the $d_j\times d_{-j}$ matrix $\mathbf{Y}_{t,j}$ can be projected into the $d_j$-dimensional vector
\begin{equation}\label{ytij}
\mathbf{y}_{t,i,j}=\mathbf{Y}_{t,j}\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}=w_if_{t,i}\mathbf{a}_{i,j}+\mathbf{e}_{t,i,j}\,,
\end{equation}
where $\mathbf{e}_{t,i,j}=\mathbf{E}_{t,j}\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}$ with $\mathbf{E}_{t,j}$ specified in \eqref{unfold X}. Then, it reduces to a standard vector-variate factor model with only one factor and a much lower dimension. Given $\xi_t$, a linear combination of $\mathcal{Y}_t$, for each $i \in [r]$ and $j \in [m]$, similarly to \eqref{Xi and Sigma Y}, we let
\begin{equation*}\label{Sigma small y}
\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)=\frac{1}{n-1}\sum_{t=2}^n\mathbb{E}[\{\mathbf{y}_{t,i,j}-\mathbb{E}(\bar\mathbf{y}_{i,j})\}\{\xi_{t-1}-\mathbb{E}(\bar \xi)\}]\,,
\end{equation*}
where $\bar\mathbf{y}_{i,j}=n^{-1}\sum_{t=1}^n\mathbf{y}_{t,i,j}$.
Under Assumption \ref{error} and $\min_{i \in [r]}|g_{1,i,\xi}| > 0$ for $g_{1,i,\xi}$ defined in \eqref{gki}, it holds that $ \bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)=g_{1,i,\xi}\mathbf{a}_{i,j}$ and $ \bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)/|\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)|_2 \in \{ \mathbf{a}_{i,j}, -\mathbf{a}_{i,j}\}$. In practice, $\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}$ is unknown. Based on the estimates $\{\tilde\mathbf{a}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}$, we can plug them into the definition of $\mathbf{B}_j$ to obtain $\tilde{\mathbf{B}}_j$, the estimate of $\mathbf{B}_j$. Set
$(\tilde\mathbf{b}_{1,j}^{\mathrm{\scriptscriptstyle + }},\ldots, \tilde\mathbf{b}_{\tilde r,j}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}
= (\tilde\mathbf{B}_j^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{B}_j)^{-1}\tilde\mathbf{B}_j^{\mathrm{\scriptscriptstyle \top }}$ and define
\[
\begin{split}
\tilde\mathbf{y}_{t,i,j}=\mathbf{Y}_{t,j}\tilde\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}= \underbrace{w_if_{t,i}\mathbf{a}_{i,j}(\mathbf{b}_{i,j}^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }})}_{\textup{``target" factor}} +\underbrace{\textstyle{\sum\nolimits_{\ell\ne i}}w_{\ell}f_{t,\ell}\mathbf{a}_{\ell, j}(\mathbf{b}_{\ell, j}^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }})}_{\textup{``noisy" factors}}+\underbrace{\mathbf{E}_{t,j}\tilde\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}}_{\textup{error term}}\,, ~ i\in[\tilde{r}],\, j\in[m]\,.
\end{split}
\]
Based on a similar projection, \cite{han2024cp} update their estimators using the eigenvector of $(n-k)^{-1}\sum_{t=k+1}^n(\tilde\mathbf{y}_{t-k,i,j}\tilde\mathbf{y}_{t,i,j}^\mathrm{\scriptscriptstyle \top }+\tilde\mathbf{y}_{t,i,j}\tilde\mathbf{y}_{t-k,i,j}^\mathrm{\scriptscriptstyle \top })$ associated with the largest eigenvalue for some $k\ge 1$. However, the ``target" and ``noisy" factors will interact with each other in their procedure. Therefore, they require the assumption of uncorrelated factors, i.e., $\mathbb{E}(f_{t,i}f_{t-k,j})=0$ for all $i\neq j$ and $k\ge 1$, to ensure the iteration works. When the factors are correlated, the iterative method in \cite{han2024cp} becomes inefficient, as shown in Figure \ref{fig:iter-step} in Section \ref{sec:numerical}. This motivates us to explore a new approach.
To reduce the effect of the ``noisy" factors for estimating $\mathbf{a}_{i,j}$, we need to involve a new linear combination of $\mathcal{Y}_t$, denoted by $\tilde{\xi}_{t,i}$, such that $\tilde{\xi}_{t-1,i}$ is correlated to the ``target" factor $f_{t,i}$ but almost uncorrelated to the ``noisy'' factors $f_{t,\ell}$ for $\ell \ne i$. To this end, we first estimate the factor series. Notice that $(\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }} \otimes \cdots \otimes \mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}\textup{vec}(\mathcal{Y}_t)=w_i f_{t,i}+(\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }} \otimes \cdots \otimes \mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}\textup{vec}(\mathcal{E}_t)$, where $(\mathbf{a}_{1,j}^{\mathrm{\scriptscriptstyle + }},\ldots,\mathbf{a}_{r,j}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}=(\mathbf{A}_j^\mathrm{\scriptscriptstyle \top } \mathbf{A}_j)^{-1}{\mathbf{A}}_j^{\mathrm{\scriptscriptstyle \top }}$. We therefore estimate $w_i f_{t,i}$ by $\check f_{t,i}=(\tilde\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }} \otimes \cdots \otimes \tilde\mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}\textup{vec}(\mathcal{Y}_t)$, where $(\tilde\mathbf{a}_{1,j}^{\mathrm{\scriptscriptstyle + }},\ldots,\tilde\mathbf{a}_{\tilde r,j}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }}=(\tilde\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{A}_j)^{-1}\tilde{\mathbf{A}}_j^{\mathrm{\scriptscriptstyle \top }}$ with $\tilde\mathbf{A}_j=(\tilde\mathbf{a}_{1,j},\ldots,\tilde\mathbf{a}_{\tilde r,j})$. Standardize the series and write
\begin{equation}\label{eq: standardized series}
\tilde f_{t,i}= (\check f_{t,i}- \bar{\check{f}}_{i})/\tilde{\sigma}_{\check f,i}\,,
\end{equation}
where $ \bar{\check{f}}_{i} = n^{-1}\sum_{t=1}^n \check f_{t,i}$ and $\tilde{\sigma}^2_{\check f,i} = (n-1)^{-1}\sum_{t=1}^n(\check f_{t,i}-\bar{\check{f}}_{i})^2$. If $\tilde r=1$, we let $\tilde \xi_{t,i}=\tilde f_{t,i}$ for $t\in[n]$ and $i\in[\tilde r]$. If $\tilde r\ge 2$,
let $\tilde\mathbf{f}_{i}=(\tilde f_{1,i},\ldots,\tilde f_{n-1,i})^\mathrm{\scriptscriptstyle \top }$ and $\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i}$ be a $(n-1)\times (\tilde r-1)$ matrix of which the columns are composed of $(\tilde f_{2,\ell},\ldots,\tilde f_{n,\ell})^\mathrm{\scriptscriptstyle \top }$ for $\ell \ne i$. We project $\tilde\mathbf{f}_{i}$ into the complementary space of $\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i}$ and obtain
\begin{equation}\label{eq: tilde xi}
(\tilde\xi_{1,i},\ldots,\tilde\xi_{n-1,i})^{\mathrm{\scriptscriptstyle \top }} = \{\mathbf{I}_{n-1}-\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i}(\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i}^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i})^{-1}\tilde\mathbf{F}_{\scalebox{1}{\text{-}} i}^{\mathrm{\scriptscriptstyle \top }}\}\tilde\mathbf{f}_{i}\,,
\end{equation}
where $\{\tilde\xi_{t,i}\}_{t=1}^{n-1}$ satisfies $ \sum_{t=2}^n \tilde\xi_{t-1,i} \tilde f_{t,\ell} = 0$ for $\ell\ne i$. Define
\begin{equation*}\label{eq:tildebSigmatildeybij}
\tilde\bm{\Sigma}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1) =\frac{1}{n-1}\sum_{t=2}^n(\tilde\mathbf{y}_{t,i,j}-\bar{\tilde\mathbf{y}}_{i,j})\tilde\xi_{t-1,i}
\end{equation*}
with $\bar{\tilde\mathbf{y}}_{i,j} = n^{-1}\sum_{t=1}^n\tilde\mathbf{y}_{t,i,j}$.
Considering that the loading vector $\mathbf{a}_{i,j}$ may be sparse, we can update the estimator $\tilde\mathbf{a}_{i,j}$ with $T_{\delta_{2,j}}\{\tilde\bm{\Sigma}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1)\}/|T_{\delta_{2,j}}\{\tilde\bm{\Sigma}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1)\}|_2$, where $\delta_{2,j}\geq0$ is the threshold level.
Such a double projection refinement can be naturally designed into an iterative procedure, as shown in Algorithm \ref{alg1}. When $\tilde r=r$, for each given $j\in[m]$, Theorem \ref{thm: iterative} in Section \ref{sec: iterative theorem} shows that the iterative estimator $\{\hat\mathbf{a}_{i,j}\}_{i\in[\tilde r]}$ obtained in Algorithm \ref{alg1} is consistent to $\{\mathbf{a}_{i,j}\}_{i\in[r]}$ up to the reflection and permutation indeterminacy.
\begin{algorithm}[H]
\footnotesize
\caption{Double projection iterations for estimating $\{\mathbf{A}_j\}_{j=1}^m$}\label{alg1}
\begin{algorithmic}[1]
\Require tensor observations $\{\mathcal{Y}_t\}_{t=1}^n$, number of factors $\tilde r$, initial estimates $\tilde\mathbf{A}^{(0)}_j = (\tilde\mathbf{a}^{(0)}_{1,j},\ldots,\tilde\mathbf{a}^{(0)}_{\tilde r,j} )$ for $j \in [m]$, maximal number of iterations $L$, upper error bound $\epsilon_0$, threshold levels $\{\delta_{2,j}\}_{j=1}^m$.
\Ensure Iterative estimates $\{\hat\mathbf{A}_j\}_{j=1}^m$
\State (\textbf{Initialization}) $\{ (\tilde\mathbf{a}_{1,j}^{(0)})^\mathrm{\scriptscriptstyle + },\ldots,(\tilde\mathbf{a}_{\tilde r,j}^{(0)})^\mathrm{\scriptscriptstyle + }\}^{\mathrm{\scriptscriptstyle \top }} \gets \{(\tilde\mathbf{A}^{(0)}_j)^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{A}^{(0)}_j\}^{-1}(\tilde\mathbf{A}^{(0)}_j)^{\mathrm{\scriptscriptstyle \top }}$ for $j\in[m]$, $v=1$, $\epsilon^{(0)}=2\epsilon_0$;
\While {$v \le L$ and $\epsilon^{(\textit{v}-1)}>\epsilon_0$}
\For {$j=1$ to $m$}
\For {$t=1$ to $n$, $i=1$ to $\tilde r$}
\STATE $
\check f^{(\textit{v},j)}_{t,i} \gets \{(\tilde\mathbf{a}_{i,m}^{(\textit{v}-1)})^{\mathrm{\scriptscriptstyle + }} \otimes \cdots \otimes (\tilde\mathbf{a}_{i,j}^{(\textit{v}-1)})^{\mathrm{\scriptscriptstyle + }} \otimes (\tilde\mathbf{a}_{i,j-1}^{(\textit{v})})^{\mathrm{\scriptscriptstyle + }} \otimes \cdots \otimes (\tilde\mathbf{a}_{i,1}^{(\textit{v})})^{\mathrm{\scriptscriptstyle + }}\}^{\mathrm{\scriptscriptstyle \top }}\textup{vec}(\mathcal{Y}_t)\,$;
\STATE replace $\check f_{t,i}$ in \eqref{eq: standardized series} with $\check f^{(\textit{v},j)}_{t,i}$ to obtain the standardized factors $\tilde f^{(\textit{v},j)}_{t,i}$;
\EndFor
\STATE $\tilde\mathbf{b}^{(\textit{v})}_{i,j} \gets \tilde\mathbf{a}^{(\textit{v}-1)}_{i,m} \otimes \cdots\otimes \tilde\mathbf{a}^{(\textit{v}-1)}_{i,j+1} \otimes \tilde\mathbf{a}^{(\textit{v})}_{i,j-1} \otimes \cdots \otimes \tilde\mathbf{a}^{(\textit{v})}_{i,1}$, $i \in [\tilde r]$;
\STATE $\tilde \mathbf{B}^{(\textit{v})}_j \gets (\tilde{\mathbf{b}}^{(\textit{v})}_{1,j} ,\ldots,\tilde{\mathbf{b}}^{(\textit{v})}_{\tilde r,j})$, $\{(\tilde\mathbf{b}^{(\textit{v})}_{1,j})^{\mathrm{\scriptscriptstyle + }} ,\ldots,(\tilde\mathbf{b}^{(\textit{v})}_{\tilde r,j})^{\mathrm{\scriptscriptstyle + }}\}^{\mathrm{\scriptscriptstyle \top }}\gets\{(\tilde\mathbf{B}^{(\textit{v})}_j)^\mathrm{\scriptscriptstyle \top } \tilde\mathbf{B}^{(\textit{v})}_j \}^{-1}(\tilde\mathbf{B}^{(\textit{v})}_j)^{\mathrm{\scriptscriptstyle \top }} $;
\For {$i=1$ to $\tilde r$}
\STATE if $\tilde r=1$, $\tilde \xi_{t,i}^{(\textit{v},j)}\gets \tilde f^{(\textit{v},j)}_{t,i},t\in[n],i\in[\tilde r]$; if $\tilde r\ge 2$, replace $\tilde f_{t,i}$ in \eqref{eq: tilde xi} with $\tilde f^{(\textit{v},j)}_{t,i}$ to obtain $\{\tilde\xi^{(\textit{v},j)}_{t,i}\}_{t=1}^{n-1}$;
\STATE $\tilde\mathbf{y}^{(\textit{v})}_{t,i,j} \gets\mathbf{Y}_{t,j}(\tilde\mathbf{b}_{i,j}^{(\textit{v})})^{\mathrm{\scriptscriptstyle + }}$ for $t\in[n]$;
\STATE $\tilde\bm{\Sigma}^{(\textit{v},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1) \gets(n-1)^{-1}\sum_{t=2}^n\{\tilde\mathbf{y}^{(\textit{v})}_{t,i,j}-n^{-1}\sum_{s=1}^n\tilde\mathbf{y}^{(\textit{v})}_{s,i,j}\}\tilde\xi^{(\textit{v},j)}_{t-1,i}$;
\STATE $ \tilde\mathbf{a}^{(\textit{v})}_{i,j}\gets T_{\delta_{2,j}}\{\tilde\bm{\Sigma}^{(\textit{v},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1)\} / |T_{\delta_{2,j}}\{\tilde\bm{\Sigma}^{(\textit{v},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1) \}|_2$;
\EndFor
\State $ \tilde\mathbf{A}^{(\textit{v})}_j \gets (\tilde\mathbf{a}^{(\textit{v})}_{1,j} ,\ldots,\tilde\mathbf{a}^{(\textit{v})}_{\tilde r,j} )$, $\{(\tilde\mathbf{a}^{(\textit{v})}_{1,j})^{\mathrm{\scriptscriptstyle + }} ,\ldots,(\tilde\mathbf{a}^{(\textit{v})}_{\tilde r,j})^{\mathrm{\scriptscriptstyle + }}\}^{\mathrm{\scriptscriptstyle \top }} \gets \{(\tilde\mathbf{A}^{(\textit{v})}_j)^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{A}^{(\textit{v})}_j\}^{-1}(\tilde\mathbf{A}^{(\textit{v})}_j)^{\mathrm{\scriptscriptstyle \top }}$;
\EndFor
\STATE $\epsilon^{(\textit{v})}\gets \max_{j\in[m]}\max_{\ell\in[\tilde r]}\min_{i\in[\tilde r]}\{1-|(\tilde\mathbf{a}^{(\textit{v})}_{i,j})^\mathrm{\scriptscriptstyle \top }\tilde\mathbf{a}^{(\textit{v} - 1)}_{\ell, j}|^2\}$, $v\gets v+1$;
\EndWhile
\STATE $\hat\mathbf{a}_{i,j}\gets \tilde\mathbf{a}^{(\textit{v})}_{i,j} $, $\hat\mathbf{A}_j\gets (\hat\mathbf{a}_{1,j},\ldots,\hat\mathbf{a}_{\tilde r,j})$ for $i\in[\tilde r]$ and $j\in[m]$.
\end{algorithmic}
\end{algorithm}
\subsection{Inference procedure for the factor loading vectors}\label{sec:inference procedure}
To introduce the main idea of our inference procedure based on $\{\hat\mathbf{a}_{i,j}\}_{i\in[\tilde{r}]}$, we assume $\hat{\mathbf{a}}_{i,j}$ is consistent to $\mathbf{a}_{i,j}$ for each $i\in[\tilde{r}]$ to simplify the notation. We consider a function $\mathbf{R}_{i,j}^*(\cdot): \mathbb{R}^{d_j}\rightarrow\mathbb{R}^{d_j}$ defined as
$$
\mathbf{R}_{i,j}^*(\mathbf{a})= \bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1) - \{\mathbf{a}^{\mathrm{\scriptscriptstyle \top }}\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)\}\mathbf{a}\,,\quad\mathbf{a} \in \mathbb{R}^{d_j}\,.
$$
Since $\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1) = g_{1,i,\xi}\mathbf{a}_{i,j}$ and $ | \mathbf{a}_{i,j} |_2 = 1$, we have $\mathbf{R}^*_{i,j}(\mathbf{a}_{i,j}) = {\bf 0}$.
For any deterministic vector $\mathbf{h}\in\mathbb{R}^{d_j}$, under some regularity conditions, it follows from the Taylor expansion that
\begin{align}\label{eq: debias taylor iterative}
\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg\{\frac{\partial \mathbf{R}^*_{i,j}(\mathbf{a}_{i,j})}{\partial \mathbf{a}^\mathrm{\scriptscriptstyle \top }}\bigg\}^{-1}\mathbf{R}_{i,j}^*(\hat\mathbf{a}_{i,j})&= \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg\{\frac{\partial \mathbf{R}^*_{i,j}(\mathbf{a}_{i,j})}{\partial \mathbf{a}^\mathrm{\scriptscriptstyle \top }}\bigg\}^{-1}\mathbf{R}_{i,j}^*(\mathbf{a}_{i,j})+\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j})\,\\
&\quad +O(|\mathbf{h}|_2|\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}|_2^2)\,. \nonumber
\end{align}
Since $\mathbf{R}^*_{i,j}(\mathbf{a}_{i,j}) = {\bf 0}$,
we have
\begin{equation}\label{eq: iterative debis expansion 0}
\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\hat\mathbf{a}_{i,j} -\mathbf{a}_{i,j}) = \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg\{\frac{\partial \mathbf{R}^*_{i,j}(\mathbf{a}_{i,j})}{\partial \mathbf{a}^\mathrm{\scriptscriptstyle \top }}\bigg\}^{-1} \mathbf{R}^*_{i,j}(\hat\mathbf{a}_{i,j})+O(|\mathbf{h}|_2|\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}|_2^2)\,.
\end{equation}
Therefore, the asymptotic representation of $\mathbf{h}^\mathrm{\scriptscriptstyle \top }\hat{\mathbf{a}}_{i,j}$ is primarily driven by the leading term on the right-hand side of \eqref{eq: iterative debis expansion 0}. However, directly deriving the asymptotic distribution of this term creates significant difficulties. The thresholding technique involved in defining $T_{\delta_{2,j}}\{\tilde\bm{\Sigma}^{(\textit{v},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_{i}}(1)\}$ in Algorithm \ref{alg1} introduces additional bias whose impact is difficult to characterize, making the derivation of an asymptotic distribution intractable. To guarantee a tractable asymptotic distribution, we construct a quantity $\hat{\bm{\vartheta}}_{i,j}$ to account for the bias induced by thresholding, and consider the asymptotic distribution of $\mathbf{h}^\mathrm{\scriptscriptstyle \top } (\hat\mathbf{a}_{i,j} -\mathbf{a}_{i,j}- \hat{\bm{\vartheta}}_{i,j})$, where
\begin{equation}\label{eq:iterative debias}
\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\hat\mathbf{a}_{i,j} -\mathbf{a}_{i,j}-\hat{\bm{\vartheta}}_{i,j})=\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg\{\frac{\partial \mathbf{R}^*_{i,j}(\mathbf{a}_{i,j})}{\partial \mathbf{a}^\mathrm{\scriptscriptstyle \top }}\bigg\}^{-1} \mathbf{R}^*_{i,j}(\hat\mathbf{a}_{i,j})-\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\hat{\bm{\vartheta}}_{i,j}+O(|\mathbf{h}|_2|\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}|_2^2)\,.
\end{equation}
Notice that
\begin{align}\label{eq: iterative debis expansion}
\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg\{\frac{\partial \mathbf{R}^*_{i,j}(\mathbf{a}_{i,j})}{\partial \mathbf{a}^\mathrm{\scriptscriptstyle \top }}\bigg\}^{-1} \mathbf{R}^*_{i,j}(\hat\mathbf{a}_{i,j})&= \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}\bigg[\frac{\{\hat\mathbf{a}_{i,j}^{\mathrm{\scriptscriptstyle \top }}\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)\}\hat\mathbf{a}_{i,j}-\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)}{\mathbf{a}_{i,j}^{\mathrm{\scriptscriptstyle \top }}\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)}\bigg]\, \\
&\quad +O(|\mathbf{h}|_2|\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}|_2^2)\,. \nonumber
\end{align}
To obtain a tractable asymptotic distribution in \eqref{eq:iterative debias},
we construct the bias-correction term $ \hat{\bm{\vartheta}}_{i,j}$ based on the leading term on the right-hand side of
\eqref{eq: iterative debis expansion}. Specifically, we replace the unknown quantities
$\mathbf{a}_{i,j}$ and $\bm{\Sigma}_{\mathbf{y}_{i,j},\xi}(1)$ with their plug-in estimators
$\hat\mathbf{a}_{i,j}$ and
$\tilde\bm{\Sigma}^{(\textit{v}_{\max},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_i}(1)$, respectively,
where $\textit{v}_{\max}$ denotes the stopping iteration of Algorithm~\ref{alg1}.
This leads to the following estimator
\[
\hat{\bm{\vartheta}}_{i,j}=
\frac{
\{\hat\mathbf{a}_{i,j}^{\mathrm{\scriptscriptstyle \top }} \tilde\bm{\Sigma}^{(\textit{v}_{\max},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_i}(1)\}
\hat\mathbf{a}_{i,j}
- \tilde\bm{\Sigma}^{(\textit{v}_{\max},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_i}(1)
}{
\hat\mathbf{a}_{i,j}^{\mathrm{\scriptscriptstyle \top }} \tilde\bm{\Sigma}^{(\textit{v}_{\max},j)}_{\tilde\mathbf{y}_{i,j},\tilde\xi_i}(1)
}\,.
\]
Theorem \ref{thm: debias iterative} in Section \ref{sec: iterative theorem} shows that
$\sqrt{n}\{w_{i}\bar{\tau}^{-1}_{i,j}(\mathbf{h})\} \mathbf{h}^\mathrm{\scriptscriptstyle \top } (\hat\mathbf{a}_{i,j} -\mathbf{a}_{i,j}- \hat{\bm{\vartheta}}_{i,j})$ is asymptotically standard normal for $\bar{\tau}_{i,j}(\mathbf{h})$ specified in \eqref{iterative not degenerate}.
Finally, we provide two estimators of the asymptotic variance $w_i^{-2}\bar{\tau}_{i,j}^2(\mathbf{h})$ so that statistical inference based on the iterative estimator can be implemented in practice. The estimation of $w_i^{-2} \bar{\tau}_{i,j}^{2}(\mathbf{h})$ is essentially a long-run variance estimation problem. Section~\ref{sec:variance iter} in the supplementary material provides an estimator $\hat{w}^{-2}_{i,j}\tilde{\tau}^{2}_{i,j}(\mathbf{h})$ for $w_i^{-2} \bar{\tau}_{i,j}^{2}(\mathbf{h})$ based on the kernel-type long-run variance estimator $\tilde{\tau}^{2}_{i,j}(\mathbf{h})$. The consistency of such kernel-type long-run variance estimator is well-known. See, for example, \cite{andrews1991heteroskedasticity} and \cite{chang2018confidence}. Therefore,
\begin{equation}\label{eq:normality-plugin-longrun}
\sqrt{n}\, \{\hat w_{i,j}\tilde\tau^{-1}_{i,j}(\mathbf{h}) \} \mathbf{h}^\mathrm{\scriptscriptstyle \top }(\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}-\hat{\bm{\vartheta}}_{i,j})
\overset{{\rm d}}{\rightarrow}
\mathcal{N}(0,1)\,.
\end{equation}
Furthermore, if the error process $\{\mathcal{E}_t\}_{t\ge 1}$ is independent of the factor process $\{\mathbf{f}_t\}_{t\ge 1}$, the asymptotic variance $w_i^{-2}\bar{\tau}_{i,j}^2(\mathbf{h})$ admits a simple form, which motivates a plug-in estimation method. Section~\ref{sec:variance iter} in the supplementary material further provides such plug-in estimator $\hat{w}^{-2}_{i,j}\hat{\tau}^{2}_{i,j}(\mathbf{h})$ for $w_i^{-2}\bar{\tau}_{i,j}^2(\mathbf{h})$.
Theorem \ref{thm: estimation of iteration variance} in the supplementary material establishes the consistency of this plug-in estimator, and hence
\begin{equation}\label{eq:normality-plugin}
\sqrt{n}\, \{\hat w_{i,j}\hat\tau^{-1}_{i,j}(\mathbf{h}) \} \mathbf{h}^\mathrm{\scriptscriptstyle \top }(\hat\mathbf{a}_{i,j}-\mathbf{a}_{i,j}-\hat{\bm{\vartheta}}_{i,j})
\overset{{\rm d}}{\rightarrow}
\mathcal{N}(0,1)\,.
\end{equation}
The simulation results in Table~\ref{table:var-iter-all} in the supplementary material further demonstrate the effectiveness of the proposed estimators for the asymptotic variance.
\subsection{Selection of tuning parameters}\label{sec:tuning}
There are some tuning parameters that need to be determined in our proposed methods. The key quantities include the number of factors $r$ specified in
\eqref{model cp}, the linear combination $\xi_t$ used to construct $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$
in \eqref{Xi and Sigma Y}, the lag parameter $K$ specified in \eqref{tilde Mj}, and two threshold levels:
$\delta_1$, used in the one-pass estimation as defined in
\eqref{hat Sigma kj}, and $\delta_{2,j}$, employed in the
iterative procedure described in Algorithm \ref{alg1}.
Write $d_{\min}=\min_{j\in[m]}d_j$ and $\mathbf{Y} = \{\text{vec}(\mathcal{Y}_1),\ldots,\text{vec}(\mathcal{Y}_n)\}^{\mathrm{\scriptscriptstyle \top }} \in \mathbb{R}^{n \times \prod_{j=1}^m d_j}$.
First, we determine the lag parameter $K$. As discussed in Remark 3 of \cite{chang2023modelling}, choosing a larger $K$ makes it more likely that the condition $\textup{rank}(\mathbf{M}_j)=r$ holds, since more lagged information is incorporated. On the other hand, as shown in Section 5.1 of \cite{chang2023modelling}, an excessively large $K$ may reduce the estimation accuracy of both the number of factors and the factor loading vectors. Balancing these two considerations, \cite{chang2023modelling} recommend choosing $K \le 10$ and show through simulations that the estimation performance is robust to the choice of $K$ within a moderate range. Our additional simulations, reported in Figures \ref{fig:Krobust-acc} and \ref{fig:Krobust-error} in the supplementary material, further support this recommendation. Specifically, the estimation accuracy improves as $K$ increases initially and then stabilizes, with almost no visible change once $K > 10$. Therefore, in practice, we recommend setting $K = 10$.
Second, we introduce how to determine $r$. When $\xi_t$ is specified, for given $\delta_1$ and $j\in[m]$, \cite{chang2023modelling} employ the eigenvalue-ratio (ER) method to estimate $r$ in the matrix CP-factor model ($m = 2$):
\begin{equation}\label{hat rj}
\tilde r^{(\textup{er})}_j(\delta_1) =\arg\min_{1 \le i \le \lfloor 0.5d_{\min} \rfloor}\frac{ \sigma_{i+1}(\tilde{\mathbf{M}}_{j})+c_n}{ \sigma_i(\tilde{\mathbf{M}}_{j})+c_n}\,,\quad j\in[m]\,,
\end{equation}
where $c_n \rightarrow 0^{\mathrm{\scriptscriptstyle + }}$ as $n\rightarrow\infty$, and $\tilde{\mathbf{M}}_{j}$ is the plug-in estimator of $\mathbf{M}_j$ specified above \eqref{tilde K21j}. Such defined ER method has also been used in \cite{chang2015high,chang2018principal,chang2025modeling} for solving other problems. In practice, we can set $c_n=n^{-1}\hat\sigma_0^2$ with $\hat\sigma_0^2=(n\prod_{j=1}^m d_j)^{-1}\|\mathbf{Y}\|_{\rm F}^2$.
Notice that $\textup{rank}(\mathbf{M}_j) = r$, and Theorem 1 of \cite{chang2023modelling} implies that $\mathbb{P}\{\tilde r^{(\textup{er})}_j(\delta_1) = r\} \to 1$ as $n \to \infty$ under certain regularity conditions for each $j \in [m]$. For the tensor CP-factor model \eqref{model cp} with more than two modes ($m > 2$), to aggregate the information from the estimators $\tilde r^{(\textup{er})}_j(\delta_1)$ across all modes, we may consider selecting $\tilde{r}$ as $\max_{j \in [m]}\tilde r^{(\textup{er})}_j(\delta_1)$.
However, when the factor loading vectors or the factor processes are highly correlated, the largest eigenvalue of $\tilde{\mathbf{M}}_{j}$ may be inflated relative to the remaining eigenvalues. As pointed out by \cite{brown1989number}, this may lead to the so-called ``one-factor bias'', under which the conventional ER method tends to favor a one-factor model even when the true number of factors is larger than one. This phenomenon is particularly pronounced when the sample size $n$ is small and is also consistent with the simulation results reported in Table \ref{table: rank} in Section \ref{sec:numerical}.
To avoid this issue, we suggest estimating $r$ by $\max_{j \in [m]}\tilde r_j^{(\log)}(\delta_1)$, where
\begin{equation}\label{hat rj-log}
\tilde r^{(\text{log})}_j(\delta_1)=\arg\min_{1 \le i \le \lfloor 0.5d_{\min} \rfloor}\frac{ \log\{1+\sigma_{i+1}(\tilde{\mathbf{M}}_{j})\}+c_n}{ \log\{1+\sigma_{i}(\tilde{\mathbf{M}}_{j})\}+c_n}\,,\quad j\in[m]\,,
\end{equation}
with the same setting as in \eqref{hat rj}. Table \ref{table: rank} in Section \ref{sec:numerical} shows that the logarithmic eigenvalue-ratio (log-ER) method \eqref{hat rj-log} exhibits better finite-sample performance than the ER method \eqref{hat rj}. Specifically, when there is a high degree of correlation among factor loading vectors, the ER method tends to underestimate the number of factors, whereas the log-ER method performs stably across all scenarios. Theorem~\ref{thm: factor number consistency} in the supplementary material establishes the consistency of the ER and log-ER estimators. Section~\ref{sec:misspecify r} in the supplementary material further examines the robustness of the proposed estimation procedures in Sections \ref{sec: initial} and \ref{sec: Double projection iterations} to misspecification of $r$.
Next, we consider how to select $\xi_t$. For the special case of the tensor CP-factor model \eqref{model cp} with $m = 2$, \cite{chang2023modelling} suggest selecting $\xi_t$ as the average of the principal components of $\{\text{vec}(\mathcal{Y}_t)\}_{t=1}^n$. Here, we propose a randomized projection approach to select $\xi_t$, which can be viewed as the extension of the method suggested by \cite{chang2023modelling}.
For a prescribed integer $p > 1$, perform PCA on $\mathbf{Y}$ and then obtain the first $p$ principal components, denoted by $\{\tilde{\eta}_{t,1},\ldots,\tilde{\eta}_{t,p}\}_{t=1}^n$. We then randomly generate a set of $p \times p$ orthonormal matrices $\bm{\Omega}^{(1)},\ldots,\bm{\Omega}^{(M)}$ and define $\xi^{(l)}_t = p^{-1} \sum_{i = 1}^{p}\eta^{(l)}_{t,i}$ with $(\eta^{(l)}_{t,1},\ldots,\eta^{(l)}_{t,p})^{\mathrm{\scriptscriptstyle \top }} = \bm{\Omega}^{(l)}(\tilde{\eta}_{t,1},\ldots,\tilde{\eta}_{t,p})^{\mathrm{\scriptscriptstyle \top }}$. Our goal is to choose the optimal candidate from $\xi^{(1)}_t,\ldots,\xi^{(M)}_{t}$ as the final $\xi_t$. For each $l \in [M]$ and a given $\breve r\geq r$, we obtain $\tilde{\mathbf{a}}_{i,j}(l)$ in the same manner as $\tilde{\mathbf{a}}_{i,j}$ defined in Section \ref{sec: initial} for $i \in [\breve{r}]$ and $j \in [m]$ but with replacing $(\xi_t,\tilde{r})$ by $(\xi^{(l)}_t,\breve{r})$.
Our guiding principle is to choose the index
$l$ for which the associated estimates $\{\tilde{\mathbf{a}}_{i,j}(l)\}_{i\in[\breve{r}],j\in[m]}$ are most similar to the other estimates $\{\tilde{\mathbf{a}}_{i,j}(\tilde{l})\}_{i\in[\breve{r}],j\in[m],\tilde{l}\neq l}$.
For any $l\in[M]$, we consider the measure
$$D(l) = \sum_{\tilde{l} \neq l}\sum_{i = 1}^{\breve{r}} I\bigg\{\max_{j\in[m]}\min_{\ell \in[\breve{r}]}[1-|\{\tilde{\mathbf{a}}_{i,j}(l)\}^\mathrm{\scriptscriptstyle \top } \tilde{\mathbf{a}}_{\ell,j}(\tilde{l}) |^2] < \varepsilon \bigg\}
\,,$$
where $\varepsilon>0$ is a prescribed distance threshold. The measure $D(l)$ quantifies the similarity between $\{\tilde{\mathbf{a}}_{i,j}(l)\}_{i\in[\breve{r}],\,j\in[m]}$ and $\{\tilde{\mathbf{a}}_{i,j}(\tilde{l})\}_{i\in[\breve{r}],\,j\in[m]}$ with $\tilde{l}\neq l$, where larger values of $D(l)$ indicate higher similarity. We then select $\xi_t$ as $\xi^{(l^*)}_t$ with $l^* = \arg\max_{l \in [M]} D(l)$. In practice, we set $p = 10$, $\breve r = 2 \tilde{r}^{*}$, $M = 50$ and $\varepsilon = 0.1$, where $\tilde{r}^{*}$ is an initial estimate of $r$ obtained via the log-ER method \eqref{hat rj-log} with $\xi_t$ selected using the approach proposed in \cite{chang2023modelling}.
Finally, we determine the threshold levels $\delta_{1}$ and $\{\delta_{2,j}\}_{j = 1}^m$.
Let $\tilde v^{(\textup{log})}_j(\delta_1)$ be the minimal ratio in \eqref{hat rj-log} corresponding to the $j$-th mode for a given $\delta_1$.
We can select $\delta_1$ as
\[
\delta_1 = \arg\min_{0<\delta < 0.1\hat\sigma_0 (n^{-1}\sum_{j = 1}^m\log d_j)^{1/2}} \frac{1}{m}\sum_{j=1}^m \tilde v^{(\textup{log})}_j(\delta)\,.
\]
Additionally, we suggest setting $\delta_{2,j} = \tilde C_* \hat\sigma_0 (n^{-1}\log d_j)^{1/2}$ with some prescribed constant $\tilde C_* \ge 0$. Extensive simulation studies demonstrate that the performance of Algorithm \ref{alg1} introduced in Section \ref{sec: Double projection iterations} with such selected $\delta_{2,j}$ is robust with respect to $\tilde C_* \in [0,1]$. We therefore recommend setting $\tilde C_* = 1$ in practice.
\section{Numerical studies}\label{sec:numerical}
We generate the observations $\{\mathcal{Y}_t\}_{t=1}^n$ via the tensor CP-factor model \eqref{model cp}.
For each $j \in [m]$, we generate $\mathbf{A}^*_j = (\mathbf{a}^*_{1,j},\ldots,\mathbf{a}^*_{r,j}) \in \mathbb{R}^{d_j \times r}$ with elements drawn independently from the uniform distribution $U(-1,1)$, subject to the restriction $\text{rank}(\mathbf{A}^*_j) = r$, and let $\breve{\mathbf{a}}_{1,j} = \mathbf{a}^*_{1,j}$ and $\breve{\mathbf{a}}_{i,j} = \mathbf{a}^*_{i,j} + \phi \mathbf{a}^*_{i-1,j}$ for $2 \le i \le r$. For each $i \in [r]$ and $j\in[m]$, we obtain $\bar{\mathbf{a}}_{i,j}$ based on $\breve{\mathbf{a}}_{i,j}$ by randomly setting its $\lfloor s d_j \rfloor$ components to be zero, and let $\mathbf{a}_{i,j} = \bar{\mathbf{a}}_{i,j}/|\bar{\mathbf{a}}_{i,j}|_2$. Here, the parameters $s$ and $\phi$, respectively, control the sparsity of $\mathbf{a}_{i,j}$ and the correlations among $(\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j})$. We generate $\{f^*_{t,i}\}_{t = 1}^n$ for $i\in [r]$ as $r$ independent AR(1) sequences, i.e.\ $f^*_{t,i} = \beta_i f^*_{t-1,i} + v_{t,i}$, where the innovations $v_{t,i}$ are independently drawn from the standard normal distribution $\mathcal{N}(0,1)$, and let $(f_{t,1},\ldots,f_{t,r})^{\mathrm{\scriptscriptstyle \top }} = \mathbf{J}^{1/2}(f^*_{t,1},\ldots,f^*_{t,r})^{\mathrm{\scriptscriptstyle \top }}$, where $\mathbf{J}$ is an $r \times r$ matrix with $[\mathbf{J}]_{i,j} = I(i = j) + \rho I(i \neq j)$. Here, the parameter $\rho$ governs the correlation among the factor processes. The elements of the error term sequence $\{\mathcal{E}_t\}_{t=1}^n$ are independently drawn from $\mathcal{N}(0,1)$. We set $m=2$ (matrix time series), $r=3$ (three factors), $w_i = 15$, $d_j = 20$, $\beta_i = 0.85 - 0.05i$, $n \in \{400,800\}$, $s \in \{0,0.3,0.6\}$, $\phi \in \{0.25,0.75\}$ and $\rho \in \{0,0.75\}$. We follow the methods described in Section \ref{sec:tuning} to select the tuning parameters involved in our proposed methods.
Table \ref{table: rank} compares the performance of two estimation methods (the ER estimator and the log-ER estimator) introduced in Section \ref{sec:tuning} and the unfolded eigenvalue-ratio (Unfolded-ER) estimator considered in \cite{chen2026estimation} for estimating $r$. We can find that log-ER outperforms ER in estimating $r$ across all scenarios, and that, except for the case $(\rho,\phi,s) = (0.75,0.75,0)$, the performance of log-ER is comparable to that of Unfolded-ER. When $(\rho,\phi,s) = (0.75,0.75,0)$, both ER and Unfolded-ER tend to underestimate $r$, whereas log-ER still maintains high accuracy in estimating $r$.
\begin{table}[htbp]
\caption{
Relative frequency estimates of $\mathbb{P}(\tilde{r} < r)$, $\mathbb{P}(\tilde{r} = r)$ and $\mathbb{P}(\tilde{r} > r)$ with $\tilde{r}$ determined by the ER estimator \eqref{hat rj}, the log-ER estimator \eqref{hat rj-log} and the Unfolded-ER estimator based on 2000 repetitions. All numbers reported below are multiplied by 100.
}
\label{table: rank}
\footnotesize
\centering
\begin{tabular}{c|c|c|c|ccc|ccc|ccc}
\hline\hline
\multirow{2}{*}{$\rho$} & \multirow{2}{*}{$\phi$} & \multirow{2}{*}{$s$} & \multirow{2}{*}{$n$} & \multicolumn{3}{c|}{\textbf{log-ER}} & \multicolumn{3}{c|}{\textbf{ER}} & \multicolumn{3}{c}{\textbf{Unfolded-ER}} \\ \cline{5-13}
& & & & $\mathbb{P}(\tilde{r} < r)$ & $\mathbb{P}(\tilde{r} = r)$ & $\mathbb{P}(\tilde{r} > r)$ & $\mathbb{P}(\tilde{r} < r)$ & $\mathbb{P}(\tilde{r} = r)$ & $\mathbb{P}(\tilde{r} > r)$ & $\mathbb{P}(\tilde{r} < r)$ & $\mathbb{P}(\tilde{r} = r)$ & $\mathbb{P}(\tilde{r} > r)$ \\ \hline
\multirow{12}{*}{0} & \multirow{6}{*}{0.25} & \multirow{2}{*}{0} & 400 & 0.25 & 99.75 & 0.00 & 6.30 & 93.70 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.65 & 99.35 & 0.00 & 7.55 & 92.45 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.3} & 400 & 0.40 & 99.60 & 0.00 & 4.80 & 95.20 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.40 & 99.60 & 0.00 & 6.85 & 93.15 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.6} & 400 & 0.10 & 99.90 & 0.00 & 2.80 & 97.20 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.55 & 99.45 & 0.00 & 6.05 & 93.95 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{2-13}
& \multirow{6}{*}{0.75} & \multirow{2}{*}{0} & 400 & 4.55 & 95.45 & 0.00 & 44.60 & 55.40 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 1.60 & 98.40 & 0.00 & 39.05 & 60.95 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.3} & 400 & 0.60 & 99.40 & 0.00 & 14.05 & 85.95 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.35 & 99.65 & 0.00 & 11.85 & 88.15 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.6} & 400 & 0.10 & 99.90 & 0.00 & 4.25 & 95.75 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.20 & 99.80 & 0.00 & 5.20 & 94.80 & 0.00 & 0.00 & 100.00 & 0.00 \\ \hline
\multirow{12}{*}{0.75} & \multirow{6}{*}{0.25} & \multirow{2}{*}{0} & 400 & 0.10 & 99.90 & 0.00 & 5.30 & 94.70 & 0.00 & 0.10 & 99.90 & 0.00 \\
& & & 800 & 0.00 & 100.00 & 0.00 & 1.25 & 98.75 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.3} & 400 & 0.00 & 100.00 & 0.00 & 2.05 & 97.95 & 0.00 & 0.00 & 100.00 & 0.00 \\
& & & 800 & 0.00 & 100.00 & 0.00 & 0.80 & 99.20 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.6} & 400 & 0.15 & 99.85 & 0.00 & 1.75 & 98.25 & 0.00 & 0.05 & 99.95 & 0.00 \\
& & & 800 & 0.00 & 100.00 & 0.00 & 0.95 & 99.05 & 0.00 & 0.00 & 100.00 & 0.00 \\ \cline{2-13}
& \multirow{6}{*}{0.75} & \multirow{2}{*}{0} & 400 & 13.80 & 86.20 & 0.00 & 71.20 & 28.80 & 0.00 & 66.25 & 33.75 & 0.00 \\
& & & 800 & 0.55 & 99.45 & 0.00 & 17.50 & 82.50 & 0.00 & 37.90 & 62.10 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.3} & 400 & 0.65 & 99.35 & 0.00 & 15.95 & 84.05 & 0.00 & 7.40 & 92.60 & 0.00 \\
& & & 800 & 0.00 & 100.00 & 0.00 & 2.55 & 97.45 & 0.00 & 0.50 & 99.50 & 0.00 \\ \cline{3-13}
& & \multirow{2}{*}{0.6} & 400 & 0.15 & 99.85 & 0.00 & 3.45 & 96.55 & 0.00 & 0.50 & 99.50 & 0.00 \\
& & & 800 & 0.00 & 100.00 & 0.00 & 1.10 & 98.90 & 0.00 & 0.05 & 99.95 & 0.00 \\ \hline\hline
\end{tabular}
\end{table}
We also compare the performance of our proposed one-pass initial estimate (Pro.init) introduced in Section \ref{sec: initial} and iterative estimate (Pro.iter) introduced in Section \ref{sec: Double projection iterations} with the composite PCA method (cPCA) and High-Order Projection Estimator (HOPE) proposed by \cite{han2024cp}, the methods of Randomized Projection PCA (RP-PCA) and Contemporary Covariance-based Iterative Simultaneous Orthogonalization (CC-ISO) proposed by \cite{chen2026estimation}, and the refined estimate for the matrix CP-factor model (RCP) proposed by \cite{chang2023modelling}. Notice that cPCA, RP-PCA and RCP are one-pass estimates, and HOPE and CC-ISO are iterative estimates. We set the tuning parameter $h = 1$ in cPCA and HOPE, as in the simulation studies of \cite{han2024cp}, and the tuning parameter $K = 10$ in RCP as suggested by \cite{chang2023modelling}. As shown in Section \ref{sec: addtional simulation results in main paper} in the supplementary material, our proposed methods are robust to the selection of $K$. For each method, the estimation error between the obtained estimates $\{\check{\mathbf{a}}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}$ and
the true factor loading vectors $\{\mathbf{a}_{i,j}\}_{i\in[r],j\in[m]}$ is measured by
\begin{equation}\label{eq:estimation error}
\psi^2(\{\check{\mathbf{a}}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}, \{\mathbf{a}_{i,j}\}_{i\in[r],j\in[m]}) = \max_{j\in[m]}\max_{\ell\in[r]}\min_{i\in[\tilde{r}]}(1-|\check{\mathbf{a}}_{i,j}^\mathrm{\scriptscriptstyle \top } \mathbf{a}_{\ell, j} |^2)\,,
\end{equation}
where $\tilde{r}$ is the associated estimate of $r$. For methods without a dedicated procedure for estimating $r$, we substitute the value obtained from the log-ER estimator when implementing their methods. As shown in Table \ref{table:rf-all}, when $\rho = 0$, Pro.iter performs comparably to CC-ISO and significantly outperforms the other methods. When $\rho = 0.75$, both CC-ISO and HOPE exhibit poor performance, whereas Pro.iter remains effective across all scenarios. Moreover, Pro.init outperforms all other one-pass estimators in all scenarios. Given $\tilde{r} = r$, we further evaluate the iterative efficiency of Pro.iter against CC-ISO and HOPE. For Pro.iter, we consider three choices of the initialization: Pro.init, cPCA, and RP-PCA. As shown in Figure \ref{fig:iter-step}, the estimation errors of Pro.iter converge to nearly zero in very few iterations across different scenarios, irrespective of the initial estimates used. However, HOPE and CC-ISO require more steps for iterative convergence. When $\rho=0.75$, the estimation errors of HOPE and CC-ISO cannot converge to zero even after a large number of iterations. This suggests that these two methods break down under such scenarios, whereas our iterative algorithm remains effective.
Recall that $\rho$ measures the degree of correlation among factors, with larger values corresponding to stronger factor correlations. The simulation results demonstrate that HOPE and CC-ISO perform poorly in scenarios with highly correlated factors. Notice that HOPE proposed by \cite{han2024cp} explicitly requires the uncorrelated factor assumption, i.e. $\mathbb{E}(f_{t,i}f_{t-k,j})=0$ for all $i\neq j$ and $k\ge 1$, while our proposed methods do not rely on this assumption.
\begin{table}
\scriptsize
\caption{
The averages and standard deviations (in parentheses) of the estimation errors \eqref{eq:estimation error} for different methods based on 2000 repetitions. Bold numbers indicate the smallest average estimation error among all competing methods. All numbers reported below are multiplied by 100.}
\centering
\setlength{6pt}{3pt}
\label{table:rf-all}
\begin{tabular}{c|c|c|c|ccc|cccc}
\hline\hline
\multirow{2}{*}{\textbf{$\rho$}} & \multirow{2}{*}{\textbf{$\phi$}} & \multirow{2}{*}{\textbf{$s$}} & \multirow{2}{*}{\textbf{$n$}} & \multicolumn{3}{c|}{\textbf{Iterative estimates}} & \multicolumn{4}{c}{\textbf{One-pass estimates}} \\ \cline{5-11}
& & & & \textbf{Pro.iter} & \textbf{HOPE} & \textbf{CC-ISO} & \textbf{Pro.init} & \textbf{cPCA} & \textbf{RP-PCA} & \textbf{RCP} \\ \hline
\multirow{12}{*}{0} & \multirow{6}{*}{0.25} & \multirow{2}{*}{0} & 400 & \textbf{0.26} (4.47) & 0.67 (7.42) & 0.75 (7.67) & 4.44 (8.52) & 17.01 (17.16) & 19.28 (17.99) & 31.77 (38.92) \\
& & & 800 & 0.63 (7.59) & 0.78 (8.34) & \textbf{0.39} (5.52) & 2.80 (8.82) & 14.40 (16.01) & 16.79 (16.27) & 27.06 (37.09) \\ \cline{3-11}
& & \multirow{2}{*}{0.3} & 400 & \textbf{0.40} (5.87) & 0.93 (9.05) & 0.47 (6.43) & 4.23 (8.86) & 14.94 (17.15) & 16.54 (17.88) & 29.74 (38.40) \\
& & & 800 & \textbf{0.36} (5.53) & 0.47 (6.23) & 0.48 (6.30) & 2.24 (6.81) & 11.72 (14.86) & 13.47 (15.81) & 26.22 (37.14) \\ \cline{3-11}
& & \multirow{2}{*}{0.6} & 400 & \textbf{0.12} (2.89) & 0.82 (8.46) & 0.69 (7.64) & 3.46 (6.20) & 13.61 (17.82) & 14.85 (17.92) & 28.83 (37.90) \\
& & & 800 & 0.50 (6.66) & 0.63 (7.46) & \textbf{0.48} (6.67) & 2.14 (7.65) & 9.99 (15.05) & 10.98 (15.13) & 25.11 (36.79) \\ \cline{2-11}
& \multirow{6}{*}{0.75} & \multirow{2}{*}{0} & 400 & 1.55 (7.31) & 1.99 (7.91) & \textbf{0.53} (3.88) & 12.15 (15.88) & 32.07 (10.40) & 33.59 (10.33) & 52.67 (31.25) \\
& & & 800 & \textbf{0.53} (4.39) & 1.35 (6.33) & 0.83 (5.01) & 5.09 (10.28) & 33.49 (10.27) & 34.49 (10.27) & 54.65 (31.55) \\ \cline{3-11}
& & \multirow{2}{*}{0.3} & 400 & \textbf{0.45} (5.50) & 0.74 (6.73) & 0.54 (5.48) & 6.59 (11.90) & 25.73 (15.45) & 28.57 (14.94) & 40.10 (38.10) \\
& & & 800 & \textbf{0.27} (4.26) & 0.57 (5.52) & 0.64 (5.76) & 2.92 (7.37) & 26.25 (14.52) & 30.18 (14.13) & 42.32 (39.31) \\ \cline{3-11}
& & \multirow{2}{*}{0.6} & 400 & \textbf{0.12} (2.62) & 0.72 (7.72) & 0.37 (5.10) & 4.17 (7.60) & 17.99 (17.77) & 19.52 (17.66) & 31.90 (38.83) \\
& & & 800 & \textbf{0.20} (4.13) & 0.40 (5.54) & 0.30 (4.73) & 2.10 (7.01) & 15.18 (16.13) & 17.47 (16.75) & 28.32 (37.64) \\ \hline
\multirow{12}{*}{0.75} & \multirow{6}{*}{0.25} & \multirow{2}{*}{0} & 400 & \textbf{0.37} (2.51) & 24.74 (37.81) & 27.41 (38.74) & 8.75 (13.39) & 48.35 (14.97) & 49.65 (14.48) & 22.31 (25.08) \\
& & & 800 & \textbf{0.12} (0.05) & 23.64 (37.38) & 24.29 (37.29) & 4.50 (8.96) & 48.54 (14.26) & 48.89 (13.65) & 21.08 (25.36) \\ \cline{3-11}
& & \multirow{2}{*}{0.3} & 400 & \textbf{0.22} (0.11) & 27.39 (39.77) & 29.57 (40.26) & 7.80 (12.18) & 49.57 (16.08) & 50.56 (15.15) & 21.85 (25.60) \\
& & & 800 & \textbf{0.09} (0.04) & 29.66 (40.68) & 29.79 (40.63) & 3.58 (7.33) & 50.31 (15.64) & 51.27 (15.06) & 19.13 (23.97) \\ \cline{3-11}
& & \multirow{2}{*}{0.6} & 400 & \textbf{0.32} (3.66) & 30.52 (41.50) & 31.97 (41.82) & 6.74 (10.93) & 51.31 (18.01) & 52.42 (16.77) & 20.94 (24.68) \\
& & & 800 & \textbf{0.08} (0.03) & 32.94 (42.23) & 32.56 (41.97) & 3.15 (6.75) & 51.05 (17.24) & 52.10 (16.47) & 20.31 (25.42) \\ \cline{2-11}
& \multirow{6}{*}{0.75} & \multirow{2}{*}{0} & 400 & \textbf{4.07} (9.69) & 6.57 (15.32) & 30.07 (21.34) & 21.52 (18.99) & 38.28 (8.65) & 42.07 (9.09) & 27.61 (20.44) \\
& & & 800 & \textbf{0.34} (1.39) & 4.02 (14.52) & 18.10 (21.57) & 10.20 (13.27) & 38.59 (9.06) & 39.92 (8.32) & 23.73 (20.67) \\ \cline{3-11}
& & \multirow{2}{*}{0.3} & 400 & \textbf{0.49} (2.52) & 13.68 (28.22) & 17.83 (30.08) & 12.23 (16.41) & 44.22 (10.84) & 45.84 (10.59) & 23.49 (23.76) \\
& & & 800 & \textbf{0.12} (0.05) & 12.97 (27.65) & 14.70 (28.94) & 5.03 (8.92) & 44.37 (10.81) & 45.29 (10.61) & 20.89 (23.60) \\ \cline{3-11}
& & \multirow{2}{*}{0.6} & 400 & \textbf{0.28} (2.41) & 22.73 (37.14) & 24.11 (37.76) & 8.39 (13.31) & 49.42 (14.94) & 50.52 (14.37) & 22.16 (25.61) \\
& & & 800 & \textbf{0.09} (0.04) & 25.04 (38.18) & 25.42 (38.16) & 3.27 (6.39) & 50.41 (14.25) & 50.68 (13.61) & 19.13 (23.68) \\ \hline\hline
\end{tabular}
\end{table}
\begin{figure}[htbp]
\centerline{\includegraphics[width= 12cm]{plot/iter-step-20260526.png}}
\caption{The lineplots for the averages of estimation errors \eqref{eq:estimation error} in the first 20 iterative steps based on 2000 repetitions. The sample size $n = 400$. The legend is defined as follows: (i) Pro.iter initialized with Pro.init (\color{black}{$-$\scalebox{0.75}{$\blacksquare$}$-$}), (ii) Pro.iter initialized with cPCA ($\color{green}{-\vartriangle-}$), (iii) Pro.iter initialized with RP-PCA ({\color{brown}$- \circ-$}), (iv) HOPE ($\color{red}{-\blacktriangle-}$), and (v) CC-ISO ($\color{blue}{-\bullet-}$).}
\label{fig:iter-step}
\end{figure}
We also evaluate the asymptotic normality \eqref{eq:normality-plugin-longrun} and \eqref{eq:normality-plugin} for the iterative estimator $\hat{\mathbf{a}}_{i,j}$ obtained in Algorithm \ref{alg1} for two choices of $\mathbf{h}$: (i) $\mathbf{h}_{1} = (1,0,\ldots,0)^\mathrm{\scriptscriptstyle \top }$ and (ii) $\mathbf{h}_{2} = (d_j^{-1/2},\ldots, d_j^{-1/2})^\mathrm{\scriptscriptstyle \top }$. It should be noted that there exists the reflection and permutation indeterminacy between the estimates and the true factor loadings. Here we set $(i,j) = (1,1)$ and impose $z_1 = \arg\min_{i \in [\tilde r]}\{1 - |\mathbf{a}_{1,1}^{\mathrm{\scriptscriptstyle \top }}\hat{\mathbf{a}}_{i,1}|^2\}$, thereby eliminating the reflection and permutation indeterminacy between $\hat{\mathbf{a}}_{z_1,1}$ and $\text{sign}(\mathbf{a}_{1,1}^{\mathrm{\scriptscriptstyle \top }}\hat{\mathbf{a}}_{z_1,1}) \cdot \mathbf{a}_{1,1}$. We exclude the replications with $\tilde{r} \neq r$ to avoid outliers.
Figures \ref{fig:normality-h1-hist-iter-est} and \ref{fig:normality-h2-hist-iter-est}
present the histograms of $\{\hat w_{z_1,1}\hat \tau_{z_1,1}^{-1}(\mathbf{h}_k)\}\sqrt{n}\,\mathbf{h}_{k}^{\mathrm{\scriptscriptstyle \top }}\{\hat{\mathbf{a}}_{z_1,1}-\text{sign}(\mathbf{a}_{1,1}^{\mathrm{\scriptscriptstyle \top }}\hat{\mathbf{a}}_{z_1,1})\cdot \mathbf{a}_{1,1}-\hat{\bm{\vartheta}}_{z_1,1}\}$ for $k \in \{1,2\}$ based on 2000 repetitions, which verify the asymptotic normality of our iterative estimator based on the asymptotic variance estimation $\hat w_{z_1,1}^{-2} \hat\tau_{z_1,1}^{2}(\mathbf{h})$. Figures
\ref{fig:normality-h1-hist-iter-est-longrun} and
\ref{fig:normality-h2-hist-iter-est-longrun} in the supplementary material also verify the asymptotic normality of our iterative estimator based on the asymptotic variance estimation
$\hat w_{z_1,1}^{-2}\tilde\tau_{z_1,1}^{2}(\mathbf{h})$.
\begin{figure}[htbp]
\centerline{\includegraphics[width= 12cm]{plot/hist-var-iter-h1-20260526.png}}
\caption{The histograms of $\{\hat{w}_{z_1,1} \hat \tau_{z_1,1}^{-1}(\mathbf{h}_1)\} \sqrt{n}\, \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}_{1}\{\hat{\mathbf{a}}_{z_1,1} - \text{sign}(\mathbf{a}_{1,1}^{\mathrm{\scriptscriptstyle \top }}\hat{\mathbf{a}}_{z_1,1} ) \cdot \mathbf{a}_{1,1} - \hat{\bm{\vartheta}}_{z_1,1} \}$ based on 2000 repetitions. The sample size $n = 400$. The red curve plots the density of $\mathcal{N}(0,1)$.}
\label{fig:normality-h1-hist-iter-est}
\end{figure}
\begin{figure}[htbp]
\centerline{\includegraphics[width= 12cm]{plot/hist-var-iter-h2-20260526.png}}
\caption{The histograms of $\{\hat{w}_{z_1,1} \hat \tau_{z_1,1}^{-1}(\mathbf{h}_2)\} \sqrt{n}\, \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}_{2}\{\hat{\mathbf{a}}_{z_1,1} - \text{sign}(\mathbf{a}_{1,1}^{\mathrm{\scriptscriptstyle \top }}\hat{\mathbf{a}}_{z_1,1} ) \cdot \mathbf{a}_{1,1} - \hat{\bm{\vartheta}}_{z_1,1} \}$ based on 2000 repetitions. The sample size $n = 400$. The red curve plots the density of $\mathcal{N}(0,1)$.}
\label{fig:normality-h2-hist-iter-est}
\end{figure}
We finally evaluate the computational speed and cost of the proposed iterative method initialized with our one-pass estimator. We fix $n=400$ and vary $(d_1,d_2) \in \{(20,20),(40,40),$ $(60,60),(80,80)\}$.
Across multiple scenarios, we benchmark the proposed Pro.iter (initialized with Pro.init) against HOPE and CC-ISO in terms of runtime and peak RAM, accounting for both initialization and iterative-phase costs.
Figure~\ref{fig:runtime} plots the averages of runtime with standard deviation bands (shaded regions) for the proposed method, HOPE, and CC-ISO across different scenarios, based on 100 replications. Results for peak RAM are similar and can be found in Figure~\ref{fig:peak RAM} in the supplementary material. When the target tensor is low dimensional, the runtime and peak RAM of Pro.iter are comparable to HOPE and CC-ISO. As dimensionality increases, Pro.iter attains markedly shorter runtime and lower peak RAM usage than HOPE and CC-ISO.
\begin{figure}[htbp]
\centerline{\includegraphics[width= 12cm]{plot/runtime-20251102.png}}
\caption{The lineplots for the averages and standard deviations (shaded region) of runtime based on 100 repetitions. The sample size $n = 400$. The legend is defined as follows: (i) Pro.iter initialized with Pro.init (\color{black}{$-$\scalebox{0.75}{$\blacksquare$}$-$}), (ii) HOPE ($\color{red}{-\blacktriangle-}$), and (iii) CC-ISO ($\color{blue}{-\bullet-}$).}
\label{fig:runtime}
\end{figure}
\section{Real data analysis: Air pollution data}\label{sec:application}
In this section, we analyze the spatio-temporal structure of air pollution in Beijing using a multi-dimensional representation of the monitoring data. This dataset contains six hourly air-pollution variables ($\textup{PM}_{2.5},\textup{PM}_{10},\textup{SO}_2,\textup{NO}_2,\textup{CO}$ and $\textup{O}_3$) from 12 nationally controlled air-quality monitoring stations in Beijing, which can be downloaded from \url{https://archive.ics.uci.edu/dataset/501/beijing+multi+site+air+quality+data}.
The observation period spans from March 1, 2013, to February 28, 2017.
The dataset contains some missing and extreme values, which were handled through interpolation. Since our focus is on the spatio-temporal structure of air-pollution variations, we apply differencing to the hourly pollutant observations for each monitoring site and pollutant type. All series are standardized to remove the impact of different measurement scales. The aforementioned procedures result in a tensor time series $\mathcal{Y}_t = (y_{t,\ell_1,\ell_2,\ell_3})_{12 \times 6 \times 24}$ for $t \in [1461]$ (i.e. $m=3,d_1 = 12, d_2 =6, d_3 = 24, n = 1461$), where $y_{t,\ell_1,\ell_2,\ell_3}$ records the concentration change of pollutant $\ell_2$ at station $\ell_1$ during the $\ell_3$-th hour of day $t$. Figure \ref{fig:app-air-timeseries} in the supplementary material shows the time series plots of $\{\mathcal{Y}_t\}_{t=1}^n$. This tensor representation enables the exploration of multi-way dependencies in Beijing’s air quality data, revealing how pollution intensity co-varies across space, time, and pollutant dimensions.
We use the tensor CP-factor model \eqref{model cp} to fit $\{\mathcal{Y}_t\}_{t=1}^n$, where $\mathbf{a}_{i,j}$ represents the factor loading vector of the $i$-th factor in the $j$-th mode. To estimate the factor loading vectors based on our proposed method, we set the tuning parameters following Section~\ref{sec:tuning}. Using the proposed log-ER method, we obtain $\tilde{r} = 2$, indicating the presence of two latent factors. Initialized with the one-pass estimator introduced in Section~\ref{sec: initial}, the proposed iterative estimator in Algorithm~\ref{alg1} converges successfully.
Table \ref{table:app-loading-a2} presents the estimations of the factor loadings $\mathbf{a}_{i,2} \in \mathbb{R}^6$ based on Pro.iter, which reveal two main patterns of pollutant variation. The first loading vector ($i = 1$) has a very high value for O$_3$ (0.953) but small values for other pollutants, indicating that this factor mainly reflects changes in ozone concentration, which vary differently from other pollutants.
The second loading vector ($i = 2$) has positive values for PM$_{2.5}$, PM$_{10}$, SO$_2$, NO$_2$, and CO, suggesting a common pollution pattern where several pollutants increase or decrease together. Therefore, we refer to the first estimated factor as the \textit{ozone-related factor},
which mainly captures variations driven by O$_3$,
and the second as the \textit{general pollution factor},
representing the joint fluctuation of multiple pollutants.
\begin{table}[htbp]
\centering
\setlength{6pt}{6pt}
\caption{Estimations of the loading vectors $\mathbf{a}_{i,2}\in\mathbb{R}^6$ for the pollution-variable mode based on Pro.iter. Standard errors reported in parentheses are calculated based on the asymptotic variance estimation $\hat{w}_{i,j}^{-2}\hat \tau_{i,j}^{2}(\mathbf{h})$. $^{*}$, $^{**}$, and $^{***}$ indicate significance at the levels 5\%, 1\%, and 1\textperthousand, respectively, based on two-sided $t$-tests. }
\label{table:app-loading-a2}
\begin{tabular}{c|cc}
\hline\hline
Pollutant & $i=1$ & $i=2$ \\
\hline
PM$_{2.5}$ & 0.008 (0.015) & 0.659$^{***}$ (0.035) \\
PM$_{10}$ & $-$0.021 (0.013) & 0.430$^{***}$ (0.025) \\
SO$_2$ & 0.049$^{**}$ (0.017) & 0.304$^{***}$ (0.036) \\
NO$_2$ & $-$0.236$^{***}$ (0.016) & 0.289$^{***}$ (0.053) \\
CO & 0.182$^{***}$ (0.012) & 0.452$^{***}$ (0.031) \\
O$_3$ & 0.953$^{***}$ (0.002) & 0.009 (0.083) \\
\hline\hline
\end{tabular}
\end{table}
\begin{figure}[htbp]
\centering
\subfigure[ozone-related factor $(i = 1)$]{\includegraphics[width=0.45\textwidth]{plot/Rplot-1.png}}
\subfigure[general pollution factor $(i = 2)$]{\includegraphics[width=0.45\textwidth]{plot/Rplot-2.png}}
\caption{Estimations of the loading vectors $\mathbf{a}_{i,1} \in \mathbb{R}^{12}$ for the monitoring-station mode based on Pro.iter.}
\label{fig:beijing-station}
\end{figure}
Figure \ref{fig:beijing-station} illustrates the estimations of the factor loadings $\mathbf{a}_{i,1} \in \mathbb{R}^{12}$ based on Pro.iter for the monitoring-station mode. Figure \ref{fig:beijing-station}(a) shows the factor loadings of the \textit{ozone-related factor} to be nearly uniform across stations, indicating the dominance of regional photochemistry and synoptic meteorology rather than local emissions. Figure \ref{fig:beijing-station}(b) reveals pronounced spatial heterogeneity for the \textit{general pollution factor}: Dingling and Changping (northern mountains) have small loadings due to sparse population, limited sources, and effective ventilation, whereas Shunyi (eastern plain, downwind of the urban core and influenced by airport, traffic, and industry) exhibits the largest loading, consistent with higher emissions and advective transport. Moreover, as reported in Table \ref{table:app-loading-a1-all} in the supplementary material, all estimated loadings based on Pro.iter for the \textit{ozone-related factor} are statistically significant, while for the \textit{general pollution factor}, all estimated loadings are statistically significant except those corresponding to Dingling and Changping. These results provide additional inferential support for the corresponding spatial interpretation.
\begin{figure}[htbp]
\centering
\subfigure[ozone-related factor $(i = 1)$]{\includegraphics[width=0.45\textwidth]{plot/24h-a1-inf.png}}
\subfigure[general pollution factor $(i = 2)$]{\includegraphics[width=0.45\textwidth]{plot/24h-a2-inf.png}}
\caption{Estimations of the loading vectors $\mathbf{a}_{i,3} \in \mathbb{R}^{24}$ for the diurnal mode based on Pro.iter. The gray shaded region represents the pointwise 95\% confidence interval for the estimated loadings. Standard errors are calculated based on the asymptotic variance estimation $\hat{w}_{i,j}^{-2}\hat \tau_{i,j}^{2}(\mathbf{h})$.}
\label{fig:beijing-hourly}
\end{figure}
Figure \ref{fig:beijing-hourly} illustrates the estimations of the factor loadings $\mathbf{a}_{i,3}\in\mathbb{R}^{24}$ based on Pro.iter for the diurnal mode. The loading vector of the \textit{ozone-related factor} rises after sunrise, peaks around 13:00–14:00, and turns negative at night, tracking the canonical photochemical cycle of daytime production and nocturnal loss via deposition \citep{li2015diurnal}. By contrast, the \textit{general pollution factor} is distinctly bimodal—minimal before dawn, a first peak near 10:00–11:00, and a higher evening peak (around 19:00–20:00). This bimodal profile accords with established diurnal emission and mixing cycles: a morning peak from traffic and industrial start-up, and a higher evening peak from rush-hour emissions combined with boundary-layer stabilization that suppresses dispersion. Analogous morning–evening bimodality for urban aerosols (e.g., $\textup{PM}_{2.5}$ and $\textup{PM}_{10}$ in Beijing) is well documented \citep{Liu2015}. Overall, the two factors capture complementary diurnal dynamics—one driven by photochemical reactions (ozone-related) and the other by human emission activities (general pollution).
\begin{figure}[htbp]
\centering
\subfigure[ozone-related factor $(i = 1)$]{\includegraphics[width=0.45\textwidth]{plot/time-series-f1.png}}
\subfigure[general pollution factor $(i = 2)$]{\includegraphics[width=0.45\textwidth]{plot/time-series-f2.png}}
\caption{The time series plots of the two estimated latent factors based on Pro.iter. The dark blue solid line represents the 30-day one-sided simple moving average.}
\label{fig:beijing-seasonal}
\end{figure}
Figure \ref{fig:beijing-seasonal} presents the time series of the two estimated latent factors based on Pro.iter; see Section \ref{sec:factor and cp estimation} in the supplementary material for details on the factor estimation procedure. Figure \ref{fig:beijing-seasonal}(a) shows the \textit{ozone-related factor}, which exhibits a distinct seasonal cycle with higher values in summer and lower values in winter. This pattern aligns with ozone's photochemical formation mechanism, which is strongly dependent on solar radiation and temperature. Figure \ref{fig:beijing-seasonal}(b) depicts the \textit{general pollution factor}, displaying the opposite seasonal pattern—higher values in cold seasons and lower values in warm seasons. This arises because cold months see increased coal/fuel combustion for heating (releasing more fine particles and gases) and stable atmospheric conditions; in contrast, summer features stronger air movement, higher wind speeds, and frequent rainfall, which stir, disperse, and scavenge pollutants to improve air quality.
Overall, by representing the multi-site, multi-pollutant, and hourly observations as a tensor time series, the proposed CP-factor estimation method successfully identifies two dominant latent components that capture distinct physical mechanisms underlying Beijing’s air pollution. The \textit{ozone-related factor} reflects photochemical processes, whereas the \textit{general pollution factor} represents anthropogenic emission activities. Together, these two factors provide a concise and interpretable characterization of the complementary seasonal dynamics of photochemically and emission-driven pollution in Beijing.
Section~\ref{sec: addtional empirical results in main paper} in the supplementary material reports the results of the real data analysis based on the other three methods (Pro.init, HOPE, and CC-ISO). The findings suggest that the estimates of Pro.iter are more interpretable and more consistent with well-established pollution mechanisms than those produced by Pro.init, HOPE, and CC-ISO. To assess the reliability of the above empirical results, we also show in Section~\ref{sec: addtional empirical results in main paper} in the supplementary material that these results remain essentially unchanged when using the winsorized data, indicating the robustness of our conclusions.
Section \ref{sec: app-famafrench} in the supplementary material provides another real data analysis for financial data. It is well known that financial data exhibit strong cross-sectional dependence, and often involve highly correlated latent factors. Table \ref{table:app-forecast} in the supplementary material reports the average forecasting errors for financial returns based on different methods. It can be observed that the tensor CP-factor methods without uncorrelated factor assumption consistently outperform the methods that rely on this assumption, which provides further evidence for the applicability of our proposed methods in practice.
\section{Theoretical analysis}\label{sec: theoretical}
\subsection{Assumptions}\label{sec: asmp}
We first present some technical assumptions for our theoretical analysis. Assumptions \ref{tail}--\ref{gap new} are imposed to guarantee the consistency of the one-pass estimator introduced in Section \ref{sec: initial}.
Given a general consistent initial estimator, Assumptions \ref{tail}--\ref{sparsity} and \ref{cross} are required to establish the theoretical guarantees for the associated iterative estimator introduced in Section \ref{sec: Double projection iterations}.
\begin{assumption}\label{tail}
There exist universal constants $C_1>1$, $C_2>0$, and $c_1\in(0,2]$ such that $\max_{i\in[r]}\max_{t\in[n]}\mathbb{P}(|f_{t,i}|>x)\le C_1\exp(-C_2x^{c_1})$, $\max_{j\in[m]}\max_{p_j\in[d_j]}\max_{t\in[n]}\mathbb{P}(|[\mathcal{E}_t]_{p_1,\ldots,p_m}|>x)\le C_1\exp(-C_2x^{c_1})$, and $\max_{t\in[n]}\mathbb{P}(|\xi_{t}|>x)\le C_1\exp(-C_2x^{c_1})$
for any $x>0$. There also exists a universal constant $C_3>0$ such that $w_1\ge\cdots \ge w_r\ge C_3$.
\end{assumption}
\begin{assumption}\label{mixing}
Define the $\alpha$-mixing coefficients associated with the factors and error processes $\{f_{t,1},\ldots,f_{t,r},\mathcal{E}_t\}_{t\geq1}$ as
$$\alpha(k)=\sup_t\sup_{A\in \mathcal{F}_{-\infty}^t,\,B\in\mathcal{F}_{t+k}^\infty}|\mathbb{P}(A\cap B)-\mathbb{P}(A)\mathbb{P}(B)|\,, ~~~k \ge 1\,,$$
where $\mathcal{F}_{t_1}^{t_2}$ is the $\sigma$-field generated by $\{f_{t,1},\ldots,f_{t,r},\mathcal{E}_t\}_{t=t_1}^{t_2}$. There exist some universal constants $C_4>1$, $C_5> 0$ and $c_2\in(0,1]$ such that $\alpha(k)\le C_4\exp(-C_5 k^{c_2})$ for any $k\ge 1$.
\end{assumption}
\begin{assumption}\label{sparsity}
There exists a universal constant $C_6>1$ such that $C_6^{-1}\le \sigma_r(\mathbf{A}_j)\le \sigma_1(\mathbf{A}_j)\le C_6$ for any $j\in[m]$. It holds that
$ \max_{i\in[r]}|\mathbf{a}_{i,j}|_0 \le s_{j}$ for any $j \in [m]$.
\end{assumption}
As pointed out in \cite{chang2023modelling,chang2021central}, Assumptions \ref{tail} and \ref{mixing} are standard in the literature on high-dimensional data analysis, which are satisfied for a wide range of time series models. These assumptions ensure exponential-type upper bounds for the tail probabilities of the statistics concerned. Assumption \ref{tail} focuses on the cases where all the factors, idiosyncratic errors, and the linear combination $\xi_t$ have exponentially decaying tails. Assumption \ref{mixing} is a standard $\alpha$-mixing condition allowing the data to be serially dependent but not necessarily stationary. As shown in \eqref{Xi and Sigma Y}, our procedure is based on lag-$k$ auto-covariances averaged over the sampling periods, rather than on a fixed stationary auto-covariance structure. Therefore, the factor process need not be stationary where $\mathbb{E}(f_{t,i}^2)$ is allowed to vary with $t$.
Assumption \ref{sparsity} can accommodate both sparse and dense loadings within a unified formulation. For example, if $s_j \ll d_j$, Assumption \ref{sparsity} corresponds to the sparse case; if $s_j=d_j$, it covers the dense case. Sparse loadings arise naturally in several important settings in the recent factor model literature \citep{uematsu2022estimation}.
Properly handling the sparsity via thresholding can improve estimation efficiency. We only require that $\sigma_1(\mathbf{A}_j), \ldots, \sigma_r(\mathbf{A}_j)$ are uniformly bounded away from $0$ and $\infty$ for $j \in [m]$, which is weaker than the requirement $\max_{j \in [m]}\|\mathbf{A}_j^{\mathrm{\scriptscriptstyle \top }} \mathbf{A}_j - \mathbf{I}_r\|_2 < 1$ imposed in \cite{han2024cp} and \cite{chen2026estimation}.
\begin{assumption}\label{eigenvalue}
Define $\ubar\sigma_{\xi}^2=\min_{i\in[r],j\in[m]}\sigma_{i}(\mathbf{M}_{j})$ and $\bar\sigma_{\xi}^2=\max_{i\in[r],j\in[m]}\sigma_{i}(\mathbf{M}_{j})$. Assume that $\ubar\sigma_{\xi}\ge C_7$ for some universal constant $C_7>0$ and $ n^{-1/2}\bar\sigma_{\xi}w_1\ll\ubar\sigma_{\xi}^2$.
\end{assumption}
We can regard $\xi_{t}$ as a projection of $\text{vec}(\mathcal{Y}_t)$ to a lower dimension. Certainly, we expect that the projection can retain the signal of the factor process, which is guaranteed by Assumption \ref{eigenvalue}. Specifically, by Theorem 7 of \cite{horn2026positivity}, we have
$$\ubar\sigma_{\xi}^2 \ge \bigg(\min_{i \in [r]} \sum_{k = 1}^K |g_{k,i,\xi}|^2 \bigg) \bigg\{\min_{j \in [m]} \frac{\sigma_{r}^2(\mathbf{A}_j)\sigma_{r}^2(\mathbf{B}_j)}{\sigma_{1}^{2}(\mathbf{A}_j)}\bigg\}\,.$$
Therefore, for each $i \in [r]$, if there exists some $k \in [K]$ such that $|g_{k,i,\xi}|$ is bounded away from $0$, then it follows from Assumption \ref{sparsity} that $\ubar\sigma_{\xi}$ is bounded away from $0$.
According to the definition of $\mathbf{M}_{j}$, we have $\ubar\sigma_{\xi} \lesssim \max_{k\in[K],i\in[r]}|g_{k,i,\xi}|\lesssim \bar\sigma_{\xi} \lesssim w_1$. If $w_1,\ldots,w_r$ are fixed constants, Assumption \ref{eigenvalue} holds when both $\bar\sigma_{\xi}$ and $\ubar\sigma_{\xi}$ are uniformly bounded away from $0$ and $\infty$.
For the toy example mentioned in Section \ref{sec: model} with $\beta\ne 0$, if we select $\xi_t=(\prod_{j=1}^{m}d_j)^{-1}\sum_{h_1 = 1 }^{d_1}\cdots \sum_{h_m = 1 }^{d_m} [\mathcal{Y}_t]_{h_1,\ldots,h_m}$, then $\bar\sigma_{\xi}\asymp \ubar\sigma_{\xi}\asymp w_1\asymp (\prod_{j=1}^{m}d_j)^{1/2}$ provided that $|(n-k)^{-1}\sum_{t = k+1}^n\mathbb{E}[\{f_t - \mathbb{E}(\bar{f})\}\{f_{t-k} - \mathbb{E}(\bar{f})\}]|$ with $\bar{f} = n^{-1}\sum_{t=1}^n f_t$ is uniformly bounded away from $0$ and $\infty$, which implies that Assumption \ref{eigenvalue} holds automatically. For more general scenarios, as long as $\xi_{t}$ is properly selected, it is expected that $\ubar\sigma_{\xi}\asymp w_r$ and $\bar\sigma_{\xi}\asymp w_1$. Then, Assumption \ref{eigenvalue} will hold when $w_1^2\ll w_r^2 \sqrt{n}$, which is a requirement on the relative strength of the factors. A similar condition also appears in Theorem 2 of \cite{han2024cp}.
\begin{assumption}\label{gap new}
All the eigenvalues $\bar\lambda_1,\ldots,\bar\lambda_r$ in \eqref{generalized eigenequation} are uniformly bounded away from $0$ and $\infty$. Moreover, $\min_{i \neq \ell}|\bar\lambda_i - \bar\lambda_{\ell}| \ge C_8$ and $\min_{i\in[r]}|g_{2,i,\xi}|\ge C_8\ubar\sigma_{\xi}$ for some universal constant $C_8>0$, where $g_{2,i,\xi}$ is defined in \eqref{gki}.
\end{assumption}
To identify the $r$ eigenvectors of $\mathbf{K}_{1,2,j}$ defined in \eqref{Kbj} corresponding to its $r$ nonzero eigenvalues, it is crucial that these eigenvalues are distinct. This is a typical assumption in the literature on eigen-analysis; see also Condition 5 in \cite{chang2023modelling}.
\begin{assumption}\label{cross} For any deterministic vector $\bm{\beta}\in\mathbb{R}^{\prod_{j=1}^{m}d_j}$, it holds that
\[
\max_{t\in[n]}\mathbb{P}\{|\bm{\beta}^\mathrm{\scriptscriptstyle \top }{\rm vec}(\mathcal{E}_t)|>x|\bm{\beta}|_2\}\le C_1\exp(-C_2x^{c_1})
\]
for any $x > 0$, where $C_1$, $C_2$, and $c_1$ are the same constants as those in Assumption \textup{\ref{tail}}.
\end{assumption}
Assumption \ref{cross} provides a tail bound for any linear combination of the idiosyncratic error tensor. It allows for cross-sectional dependence among the errors. This assumption can hold under very general scenarios, e.g., when $\textup{vec}(\mathcal{E}_t)$ follows a multivariate Gaussian distribution with a covariance matrix bounded in spectral norm; see also Assumption 1 in \cite{han2024cp}.
\subsection{Theoretical guarantees of the proposed methods}\label{sec: iterative theorem}
Let $D_n=\prod_{j=1}^{m}d_j$ and $S_{n}=\prod_{j=1}^{m}s_{j}$, and further define
\begin{equation*}
\Pi_{n} = \frac{\bar\sigma_{\xi} }{\ubar\sigma_{\xi}^2}\bigg(\frac{S_n \log D_n}{n}\bigg)^{1/2} \,.
\end{equation*}
Set the threshold level $\delta_1=C_* (n^{-1} \log D_n)^{1/2}$ in \eqref{hat Sigma kj} for some constant $C_*>0$.
Theorem \ref{thm: aij} shows the consistency (up to the reflection and permutation indeterminacy) of the one-pass estimator $\{\tilde\mathbf{a}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}$
introduced in Section \ref{sec: initial}.
\begin{theorem}\label{thm: aij}
Under Assumptions \textup{\ref{error}}--\textup{\ref{gap new}}, if $\Pi_{n} \ll 1$ and $\tilde r=r$, there exists a permutation of $[r]$, denoted by $\{z_1,\ldots,z_r\}$, such that
$$|\tilde\mathbf{a}_{z_i,j}- \tilde\kappa_{i,j} \mathbf{a}_{i,j}|_2 = \ubar\sigma_{\xi}\,\bar\sigma_{\xi}^{-1}\, O_{\rm p}(\Pi_{n})$$ for any $i\in[r]$ and $j\in[m]$ with some constants $ \tilde\kappa_{i,j} \in\{-1,1\}$, provided that $\log D_n\ll n^c$ for some constant $c\in(0,1)$ depending only on $c_1$ and $c_2$ specified in Assumptions {\rm\ref{tail}} and {\rm\ref{mixing}}.
\end{theorem}
In Algorithm \ref{alg1}, the estimated factor series $\{\check f_{t,i}^{(\textit{v},j)}\}_{t=1}^n$ can be regarded as linear combinations of $\{\text{vec}(\mathcal{Y}_t)\}_{t=1}^n$ with some plug-in estimators of the coefficients $(\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }}\otimes \cdots\otimes\mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})$, where $(\mathbf{a}_{1,j}^{\mathrm{\scriptscriptstyle + }},\ldots,\mathbf{a}_{ r,j}^{\mathrm{\scriptscriptstyle + }})^{\mathrm{\scriptscriptstyle \top }} = (\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j)^{-1}{\mathbf{A}}_j^{\mathrm{\scriptscriptstyle \top }}$. To study the statistical error of the iterative estimator, we define the oracle linear combinations as $\xi_{t,i}=(\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }}\otimes \cdots\otimes\mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})^\mathrm{\scriptscriptstyle \top } \text{vec}(\mathcal{Y}_t)=w_if_{t,i}+(\mathbf{a}_{i,m}^{\mathrm{\scriptscriptstyle + }}\otimes \cdots\otimes\mathbf{a}_{i,1}^{\mathrm{\scriptscriptstyle + }})^\mathrm{\scriptscriptstyle \top } \text{vec}(\mathcal{E}_t)$, and write $$\xi_{t,i}^{\textup{s}}=\bigg[\frac{1}{n}\sum_{s=1}^n\mathbb{E}\{(\xi_{s,i}-\bar\xi_i)^2\}\bigg]^{-1/2}(\xi_{t,i}-\bar\xi_i)$$ with $\bar \xi_i=n^{-1}\sum_{t=1}^n \xi_{t,i}$. If $r=1$, let $\xi_{t,i}^{\textup{sp}}=\xi_{t,i}^{\textup{s}}$ for $t\in[n]$ and $i\in[r]$. If $r\ge 2$, let $\bm{\xi}_i^{\textup{s}} = (\xi_{1,i}^{\textup{s}},\ldots,\xi_{n-1,i}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}$, and $\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}}$ be a $(n-1)\times (r-1)$ matrix of which the columns are composed of $(\xi_{2,\ell}^{\textup{s}},\ldots,\xi_{n,\ell}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}$ for $\ell \ne i$. Then, following the double projection step, we define
\begin{equation}\label{xi it sp}
(\xi_{1,i}^{\textup{sp}},\ldots,\xi_{n-1,i}^{\textup{sp}})^{\mathrm{\scriptscriptstyle \top }}=\bm{\xi}_i^{\textup{s}}-\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}}[\mathbb{E}\{(\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}}\}]^{-1}\mathbb{E}\{(\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}\bm{\xi}_i^{\textup{s}}\}\,.
\end{equation}
Let $\bar\bm{\varphi}_i=(\bar\varphi_{i,1},\ldots,\bar\varphi_{i,r})^{\mathrm{\scriptscriptstyle \top }}$ be the $r$-dimensional vector with the $i$-th entry equal to 1, while the remaining $r-1$ entries form the vector $-[\mathbb{E}\{(\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}}\}]^{-1}
\mathbb{E}\{(\mathbf{F}_{\xi,\scalebox{1}{\text{-}} i}^{\textup{s}})^{\mathrm{\scriptscriptstyle \top }}\bm{\xi}_i^{\textup{s}}\}$. Then, $\xi_{t,i}^{\textup{sp}}=\xi_{t,i}^{\textup{s}}+\sum_{\ell\ne i}\bar\varphi_{i,\ell}\xi_{t+1,\ell}^{\textup{s}}$ for $t\in[n-1]$.
Set the threshold levels $\delta_{2,j}=\tilde C_*(n^{-1}\log d_j)^{1/2}$ for $j\in[m]$ in Algorithm \ref{alg1} with some sufficiently large constant $\tilde C_*>0$. Write
\[
\Phi_{n,j}=\frac{1 }{w_r}\sqrt{\frac{s_j\log d_j}{n}}~~\textrm{and}~~\bm{\Upsilon}_k = (\Upsilon_{k,i,\ell})_{r \times r}\,,
\]
where $\Upsilon_{k,i,\ell}= (n-k)^{-1} \sum_{t=k+1}^n\mathbb{E}\{(f_{t,i}-\bar f_i)(f_{t-k,\ell}-\bar f_{\ell})\}$ for $i,\ell\in[r]$ and $k\in \{0,1\}$. Let
\begin{equation}\label{auto cross f xi}
\begin{split}
&\gamma_{\max}=\max_{i\ne \ell}|\Upsilon_{1,i,\ell}|\,, ~~ \sigma_{f_i,\xi_i} = \mathbb{E}\bigg\{\frac{1}{n-1}\sum_{t=2}^n(f_{t,i}-\bar f_i)\xi_{t-1,i}^{\textup{sp}}\bigg\}\,,\\
& ~~~~~~\textup{and}~~ L_n= \bigg(\frac{\sum_{j=1}^m d_j\log d_j}{n}\bigg)^{1/2}+\frac{(\sum_{j=1}^m d_j)^{1/\tilde c}}{n} \,,
\end{split}
\end{equation}
where $\tilde c^{-1}=1+2c_1^{-1}+c_2^{-1}$ for $c_1$ and $c_2$ in Assumptions \ref{tail} and \ref{mixing}. Theorem \ref{thm: iterative} gives the convergence rate of the iterative estimator obtained by Algorithm \ref{alg1}.
\begin{theorem}\label{thm: iterative}
Let Assumptions \textup{\ref{error}}--\textup{\ref{sparsity}} and \textup{\ref{cross}} hold. Assume that $\tilde r=r$ and the initial estimates in Algorithm \textup{\ref{alg1}} satisfy $\max_{i\in[r],j\in[m]}w_r^{-1}w_1|\tilde\mathbf{a}_{z_i,j}^{(0)}- \tilde\kappa_{i,j} \mathbf{a}_{i,j}|_2 = o_{\rm p}(1)$ for some permutation $\{z_1,\ldots,z_r\}$ of $[r]$ and some constants $\tilde\kappa_{i,j} \in \{-1,1\}$. If $D_n\rightarrow \infty$ as $n\rightarrow \infty$, $C_{9}^{-1}\le\sigma_r(\bm{\Upsilon}_0)\le \sigma_1(\bm{\Upsilon}_0)\le C_{9}$ for some universal constant $C_{9}>1$,
\begin{equation}\label{strong factor condition}
\frac{w_1}{w_r^2} \bigg( \frac{\gamma_{\max} }{w_r}+\frac{1}{w_r\sqrt{n}}+ L_n\bigg) \ll 1 ~~\text{and}~~|\sigma_{f_i,\xi_i}| \ge C_{9}^{-1}
\end{equation}
for all $i\in[r]$, then we have
\[
\max_{i\in[r],j\in[m]}| \hat\mathbf{a}_{z_i,j}- \kappa_{i,j} \mathbf{a}_{i,j}|_2 = O_{\rm p}\bigg(\max_{j\in[m]}\Phi_{n,j}+\frac{\gamma_{\max}}{w_r^2}\bigg)\,
\]
for some constants $\kappa_{i,j} \in \{-1,1\}$, provided that
the number of iterations satisfies
$\textit{v}_{\max}\gtrsim-\log(\max_{j\in[m]}\Phi_{n,j}+\gamma_{\max} w_r^{-2})$ and $\max_{j \in [m]}\log d_j \ll n^c$ for some constant $c\in(0,1)$ depending only on $c_1$ and $c_2$ specified in Assumptions {\rm\ref{tail}} and {\rm\ref{mixing}}.
\end{theorem}
The requirement on the convergence rate of $\{\tilde\mathbf{a}_{i,j}^{(0)}\}_{i\in[\tilde{r}],j\in[m]}$ can be easily satisfied if taking the one-pass estimator introduced in Section \ref{sec: initial} as the initial estimator of Algorithm \ref{alg1}. As discussed below Assumption \ref{eigenvalue}, when the initial linear combination is properly selected, it is expected that $\underaccent{\bar}{\sigma}_{\xi}\asymp w_r$ and $\bar{\sigma}_{\xi}\asymp w_1$. Theorem \ref{thm: aij} implies that the convergence rate of the one-pass estimator satisfies this requirement automatically. The first part of condition \eqref{strong factor condition} is mainly to control the plug-in error of $\{\tilde{\mathbf{a}}^{(\textit{v})}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}$ in the iterations.
The requirement $w_1w_r^{-2}n^{-1}(\sum_{j=1}^m d_j)^{1/\tilde c} \ll 1$ originates from the serial dependence of the error process $\{\mathcal{E}_t\}_{t\ge 1}$, and is unnecessary if $\{\mathcal{E}_t\}_{t\ge 1}$ are serially independent sub-Gaussian tensors, and are also independent of the factor process $\{\mathbf{f}_t\}_{t \ge 1}$.
The second part of condition \eqref{strong factor condition} is similar to Assumption \ref{eigenvalue}, which requires the lag-one cross-correlation between $\{\xi_{t,i}^{\textup{sp}}\}_{t=1}^n$ and $\{f_{t,i}\}_{t=1}^n$ to be non-vanishing. The convergence rate of the iterative estimator includes two parts. The first part depends on $\Phi_{n,j}$, which is a typical rate under sparsity. The second part depends on the lag-one cross-correlations of the factors $\gamma_{\max}$ and the factor strength $w_r$, which is mainly from the estimation error of the factors when we decorrelate them in the double projection step. Under Assumption \ref{eigenvalue}, we can show that $\underaccent{\bar}{\sigma}_{\xi}\lesssim w_r$. When $\mathbb{E}(f_{t,i})=0$ for all $t,i$ and $\mathbb{E}(f_{t,i}f_{t-k,j})=0$ for all $i\ne j, k\ge 1$ as assumed in \cite{han2024cp}, we have $\gamma_{\max}=O(n^{-1})$ and the convergence rate in Theorem \ref{thm: iterative} can be simplified as $O_{\rm p}(\max_{j\in[m]}\Phi_{n,j})$, which implies the iterative estimator in this scenario is more accurate in comparison to the one-pass estimator introduced in Section \ref{sec: initial}.
Let $\bar\mathbf{e}_{i,j}=n^{-1}\sum_{t=1}^n\mathbf{e}_{t,i,j}$ with $\mathbf{e}_{t,i,j}$ defined in \eqref{ytij}. With $\xi_{t,i}^{\textup{s}}$ defined above \eqref{xi it sp}, write
\begin{equation}\label{sigma ei xi}
\begin{split}
\tilde\bm{\Sigma}_{\mathbf{e}_{\ell,j},\xi_i}(1)&=\frac{1}{n-1}\sum_{t=2}^n(\mathbf{e}_{t,\ell,j}-\bar \mathbf{e}_{\ell,j})\xi_{t-1,i}^{\textup{s}}\,,\\
\tilde\bm{\Sigma}_{\mathbf{e}_{\ell,j},\xi_i}(0)&=\frac{1}{n-1}\sum_{t=2}^n[(\mathbf{e}_{t,\ell,j}-\bar \mathbf{e}_{\ell,j})\xi_{t,i}^{\textup{s}}-\mathbb{E}\{(\mathbf{e}_{t,\ell,j}-\bar \mathbf{e}_{\ell,j})\xi_{t,i}^{\textup{s}}\}]
\end{split}
\end{equation}
for $j \in [m]$ and any $i,\ell \in [r]$. Theorem \ref{thm: debias iterative} provides a limiting representation for the iterative estimator $\hat\mathbf{a}_{i,j}$.
\begin{theorem}\label{thm: debias iterative}
Let $\mathbf{h} \in\mathbb{R}^{d_j}$ be any non-random vector satisfying $|\mathbf{h}|_2 = 1$, and the conditions in Theorem \textup{\ref{thm: iterative}} hold. For any $i\in[r]$ and $j\in[m]$, it holds that
\begin{align*}
\mathbf{h}^{\mathrm{\scriptscriptstyle \top }} ( \hat\mathbf{a}_{z_i,j} & - \kappa_{i,j} \mathbf{a}_{i,j}- \hat{\bm{\vartheta}}_{z_i,j})\\
=\,& \frac{ \kappa_{i,j}}{w_{i}\sigma_{f_{i},\xi_{i}}}\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\mathbf{I}_{d_j}-\mathbf{a}_{{i},j}\mathbf{a}_{{i},j}^{\mathrm{\scriptscriptstyle \top }})\bigg\{\tilde\bm{\Sigma}_{\mathbf{e}_{i,j},\xi_{i}}(1)+\sum_{\ell \neq i} \bar\varphi_{i,\ell} \tilde\bm{\Sigma}_{\mathbf{e}_{i,j},\xi_\ell}(0)\bigg\}\\
&+ O_{\rm p}\bigg(\frac{\gamma_{\max}}{w_{i}w_r}+\frac{1}{\sqrt{n}w_{i}w_r}+\frac{w_1}{w_{i}w_r} L_n \max_{j\in[m]}\Phi_{n,j}\bigg)
+o_{\rm p}\bigg(\frac{1}{w_i\sqrt{n}}\bigg)\,,
\end{align*}
where $z_i$ and $\kappa_{i,j}$ are specified in Theorem \textup{\ref{thm: iterative}}, and $\sigma_{f_{i},\xi_{i}}$ and $\bar\varphi_{i,\ell}$ are defined, respectively, in \eqref{auto cross f xi} and below \eqref{xi it sp}. Furthermore, if
\begin{equation}\label{iterative not degenerate}
\lim_{n\rightarrow\infty}{\rm Var}\bigg[\frac{\sqrt{n}}{\sigma_{f_{i},\xi_{i}}}\mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\mathbf{I}_{d_j}-\mathbf{a}_{i,j}\mathbf{a}_{i,j}^{\mathrm{\scriptscriptstyle \top }})\bigg\{\tilde\bm{\Sigma}_{\mathbf{e}_{i,j},\xi_{i}}(1)+\sum_{\ell\ne i}\bar\varphi_{i,\ell}\tilde\bm{\Sigma}_{\mathbf{e}_{{i},j},\xi_\ell}(0)\bigg\}\bigg] = \bar{\tau}^2_{i,j}(\mathbf{h})
\end{equation}
for some deterministic positive number $\bar{\tau}^2_{i,j}(\mathbf{h})$ and
\begin{align}
&\frac{\gamma_{\max}}{w_r}+\frac{1}{\sqrt{n}w_r}+ \frac{w_1}{w_r} L_n \max_{j\in[m]} \Phi_{n,j} \ll n^{-1/2}\,,\label{negligible error 2}
\end{align}
then
\[
\sqrt{n}\{w_{i}\bar{\tau}^{-1}_{i,j}(\mathbf{h})\} \mathbf{h}^{\mathrm{\scriptscriptstyle \top }}(\hat\mathbf{a}_{z_i,j}-\kappa_{i,j}\mathbf{a}_{i,j}- \hat{\bm{\vartheta}}_{z_i,j})\overset{{\rm d}}{\rightarrow} \mathcal{N}(0,1)\,.
\]
\end{theorem}
Condition \eqref{iterative not degenerate} is to ensure that the asymptotic variance is not degenerate.
Condition \eqref{negligible error 2} is to control the estimation error of the factors and the plug-in error of $\{\tilde{\mathbf{a}}^{(\textit{v})}_{i,j}\}_{i\in[\tilde{r}],j\in[m]}$ in the iterations. If all the factors are strong factors such that $ w_1\asymp w_r\asymp \sqrt{D_n}$, condition \eqref{negligible error 2} holds provided that $\gamma^2_{\max}\ll D_n/n$, $D_n\gg 1$, $(\max_{j\in[m]} s_j\log d_j)(\sum_{j=1}^md_j\log d_j)\ll nD_n$, and $(\max_{j\in[m]} s_j\log d_j)(\sum_{j=1}^md_j)^{2/\tilde c}\ll n^2D_n$.
\section{Discussion}\label{sec: discuss}
In this paper, we develop new estimation methods for tensor CP-factor models that explicitly exploit the tensor structure and allow for correlated factors and loadings, thereby providing useful tools for analyzing high-dimensional tensor-valued data. Several assumptions adopted in this paper can be further relaxed. Assumption \ref{error}, which requires the error process $\{\mathcal{E}_t\}_{t \ge 1}$ in \eqref{model cp} to be serially uncorrelated, is a key condition for the validity of our auto-covariance-based procedures. Once serial correlation is present in the error process, extending the proposed methods is challenging. In particular, the key identity $\mathbf{K}_{1,2,j} = \mathbf{A}_j\mathbf{G}_{1,\xi}\mathbf{G}_{2,\xi}^{-1}(\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }\mathbf{A}_j)^{-1}\mathbf{A}_j^\mathrm{\scriptscriptstyle \top }$ with $\mathbf{A}_j=(\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j})$ used to identify the factor loading vectors $\mathbf{a}_{1,j},\ldots,\mathbf{a}_{r,j}$ does not hold. How to identify and estimate $\{\mathbf{a}_{i,j}\}_{j\in[m],i \in [r]}$ in the setting with serially correlated error process deserves further investigation. Section \ref{sec: relax error serial dependence} in the supplementary material provides some further discussion for this. We also discuss in Section \ref{sec: relax tail condition} in the supplementary material that our theoretical results can be extended from the exponential-decay assumptions in Assumptions \ref{tail} and \ref{mixing} to polynomial-decay conditions.
Assumption \ref{mixing} requires weak serial dependence among the observed tensor process $\{\mathcal{Y}_t\}_{t \ge 1}$ which does not cover the cases with unit-root tensor process. It would be interesting to extend the proposed methods to handle unit-root tensor process $\{\mathcal{Y}_t\}_{t \ge 1}$. We plan to investigate it in our future research.
\bibliographystyle{jasa}
\begingroup
\setlength{\bibsep}{2pt}
\linespread{0.9}\selectfont
\bibliography{Ref-abbreviation}
\endgroup
\clearpage
\newpage