Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
117,161 characters · 13 sections · 80 citation commands
CP-Factorization for High Dimensional Tensor Time Series and Double Projection Iterations
\if11 { \spacingset{1.25}
\affil[1]{\it Joint Laboratory of Data Science and Business Intelligence, Institute of Statistical Interdisciplinary Research, Southwestern University of Finance and Economics, Chengdu, China} \affil[2]{\it State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China} \affil[3]{\it Department of Statistics, The London School of Economics and Political Science, London, U.K.} \affil[4]{\it School of Statistics and Data Science, Institute of Data Science and Statistics, Shanghai University of Finance and Economics, Shanghai, China}
\setcounter{Maxaffil}{0}
} \fi \if01 {
} \fi
\spacingset{1.5}
{\sl Keywords}: CP decomposition; dimension reduction; double projection iteration; statistical inference; tensor time series.
\spacingset{1.69} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip}
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 zhou2013tensor, recommendation systems bi2018, and dynamic transportation networks 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 kolda2009tensor. The Tucker decomposition is often achieved by the SVD on the unfolded matrices de2000multilinear. By contrast, computing the CP decomposition is NP-hard, and the alternating least squares iteration remains as the workhorse method 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 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 wang2019factor, where a two-way factor structure was introduced for matrix time series (i.e. a tensor with two modes). See also elynn2020, yu2022projected, and chen2023statistical. Extensions to higher-order tensor Tucker-factor models have been studied in chen2024rank, han2024tensor, chen2024semi, barigozzi2023statistical, and HE2026105557. Robust estimation methods for tensor Tucker-factor models are further considered in barigozzi2023robust and 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) below. For tensor CP-factor models, 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. 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), chang2023modelling propose a one-pass estimation procedure (i.e. without iterations) for matrix CP-factor models. Furthermore the “linearly independent" requirement is freed in 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 chang2023modelling is based on a generalized eigen-analysis while the procedure of 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 chang2023modelling nor 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 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 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 R-function CP_TTS for implementing our newly proposed methods is available publicly in the HDTSA package chang2024hdtsa.
The rest of the paper is organized as follows. Section (ref) introduces the tensor CP-factor models and the basic settings. Section (ref) presents our methodology, including the one-pass and iterative estimators together with the inference procedure. Sections (ref) and (ref) validate the performance of our proposed methods through simulation studies and a real data analysis, respectively. Section (ref) develops the theoretical guarantees of the proposed estimators. Section (ref) 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}.
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.
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 han2024cp as follows:
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 (ref) is an extension of the model considered in 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 (ref) 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 (ref) 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 (ref).
For each $j\in[m]$, write $d_{\scalebox{1}{\text{-}} j} = \prod_{j^\prime\ne j}^md_{j^\prime}$. Following the tensor matricization in kolda2009tensor, we can reshape $\mathcal{Y}_t$ into a $d_j\times d_{\scalebox{1}{\text{-}} j}$ matrix as follows:
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 (ref) can be rewritten as the matrix CP-factor model considered in chang2023modelling:
We assume $\textup{rank}(\mathbf{A}_j) = r$ for $j \in [m]$, which is also required in han2024cp. Under this assumption, Lemma (ref) 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 chang2023modelling. However, for the tensors with more than two modes, the matrix $\mathbf{B}_j$ in (ref) involves a Kronecker product structure. Direct application of the method in chang2023modelling would discard this structural information, leading to substantial statistical inefficiency. Meanwhile, the two-stage estimation approach of 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). Moreover, simulation studies in Section (ref) show that even for matrix-variate cases, our newly proposed methods perform significantly better than the procedure in chang2023modelling.
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
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
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
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 (ref). 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
In practice, $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ and $\mathbf{Q}_j$ in (ref) are unknown. Given observations $\{\mathcal{Y}_t\}_{t=1}^n$, in the spirit of bickel2008covariance, we can estimate $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ by
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
for some predetermined integer $K\ge 1$. Under Assumption (ref) 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). 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 (ref), we can estimate $\mathbf{K}_{1,2,j}$ by
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 (ref) 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) in Section (ref) 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.
The one-pass estimation procedure entails thresholding the large $d_j\times d_{-j}$ matrix $\tilde\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ in (ref), 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 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 (ref) 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
where $\mathbf{e}_{t,i,j}=\mathbf{E}_{t,j}\mathbf{b}_{i,j}^{\mathrm{\scriptscriptstyle + }}$ with $\mathbf{E}_{t,j}$ specified in (ref). 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 (ref), we let
where $\bar\mathbf{y}_{i,j}=n^{-1}\sum_{t=1}^n\mathbf{y}_{t,i,j}$. Under Assumption (ref) and $\min_{i \in [r]}|g_{1,i,\xi}| > 0$ for $g_{1,i,\xi}$ defined in (ref), 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 \[
\] Based on a similar projection, 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 han2024cp becomes inefficient, as shown in Figure (ref) in Section (ref). 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
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
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
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). When $\tilde r=r$, for each given $j\in[m]$, Theorem (ref) in Section (ref) shows that the iterative estimator $\{\hat\mathbf{a}_{i,j}\}_{i\in[\tilde r]}$ obtained in Algorithm (ref) is consistent to $\{\mathbf{a}_{i,j}\}_{i\in[r]}$ up to the reflection and permutation indeterminacy.
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
Since $\mathbf{R}^*_{i,j}(\mathbf{a}_{i,j}) = {\bf 0}$, we have
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 (ref). 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) 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
Notice that
To obtain a tractable asymptotic distribution in (ref), we construct the bias-correction term $ \hat{\bm{\vartheta}}_{i,j}$ based on the leading term on the right-hand side of (ref). 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). 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) in Section (ref) 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 (ref).
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) 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, andrews1991heteroskedasticity and chang2018confidence. Therefore,
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) 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) in the supplementary material establishes the consistency of this plug-in estimator, and hence
The simulation results in Table (ref) in the supplementary material further demonstrate the effectiveness of the proposed estimators for the asymptotic variance.
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 (ref), the linear combination $\xi_t$ used to construct $\bm{\Sigma}_{\mathbf{Y}_j,\xi}(k)$ in (ref), the lag parameter $K$ specified in (ref), and two threshold levels: $\delta_1$, used in the one-pass estimation as defined in (ref), and $\delta_{2,j}$, employed in the iterative procedure described in Algorithm (ref). 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 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 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, 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) and (ref) 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]$, chang2023modelling employ the eigenvalue-ratio (ER) method to estimate $r$ in the matrix CP-factor model ($m = 2$):
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 (ref). Such defined ER method has also been used in 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 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 (ref) 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 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) in Section (ref). To avoid this issue, we suggest estimating $r$ by $\max_{j \in [m]}\tilde r_j^{(\log)}(\delta_1)$, where
with the same setting as in (ref). Table (ref) in Section (ref) shows that the logarithmic eigenvalue-ratio (log-ER) method (ref) exhibits better finite-sample performance than the ER method (ref). 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) in the supplementary material establishes the consistency of the ER and log-ER estimators. Section (ref) in the supplementary material further examines the robustness of the proposed estimation procedures in Sections (ref) and (ref) to misspecification of $r$.
Next, we consider how to select $\xi_t$. For the special case of the tensor CP-factor model (ref) with $m = 2$, 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 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) 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 (ref) with $\xi_t$ selected using the approach proposed in 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 (ref) 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) introduced in Section (ref) 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.
We generate the observations $\{\mathcal{Y}_t\}_{t=1}^n$ via the tensor CP-factor model (ref). 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) to select the tuning parameters involved in our proposed methods.
Table (ref) compares the performance of two estimation methods (the ER estimator and the log-ER estimator) introduced in Section (ref) and the unfolded eigenvalue-ratio (Unfolded-ER) estimator considered in 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$.
We also compare the performance of our proposed one-pass initial estimate (Pro.init) introduced in Section (ref) and iterative estimate (Pro.iter) introduced in Section (ref) with the composite PCA method (cPCA) and High-Order Projection Estimator (HOPE) proposed by han2024cp, the methods of Randomized Projection PCA (RP-PCA) and Contemporary Covariance-based Iterative Simultaneous Orthogonalization (CC-ISO) proposed by chen2026estimation, and the refined estimate for the matrix CP-factor model (RCP) proposed by 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 han2024cp, and the tuning parameter $K = 10$ in RCP as suggested by chang2023modelling. As shown in Section (ref) 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
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), 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), 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 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.
We also evaluate the asymptotic normality (ref) and (ref) for the iterative estimator $\hat{\mathbf{a}}_{i,j}$ obtained in Algorithm (ref) 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) and (ref) 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) and (ref) 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})$.
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) 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) 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.
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) 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 (ref) 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). 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), the proposed iterative estimator in Algorithm (ref) converges successfully.
Table (ref) 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 ozone-related factor, which mainly captures variations driven by O$_3$, and the second as the general pollution factor, representing the joint fluctuation of multiple pollutants.
Figure (ref) 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)(a) shows the factor loadings of the ozone-related factor to be nearly uniform across stations, indicating the dominance of regional photochemistry and synoptic meteorology rather than local emissions. Figure (ref)(b) reveals pronounced spatial heterogeneity for the 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) in the supplementary material, all estimated loadings based on Pro.iter for the ozone-related factor are statistically significant, while for the 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.
Figure (ref) 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 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 li2015diurnal. By contrast, the 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 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).
Figure (ref) presents the time series of the two estimated latent factors based on Pro.iter; see Section (ref) in the supplementary material for details on the factor estimation procedure. Figure (ref)(a) shows the 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)(b) depicts the 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 ozone-related factor reflects photochemical processes, whereas the 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) 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) in the supplementary material that these results remain essentially unchanged when using the winsorized data, indicating the robustness of our conclusions.
Section (ref) 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) 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.
We first present some technical assumptions for our theoretical analysis. Assumptions (ref)--(ref) are imposed to guarantee the consistency of the one-pass estimator introduced in Section (ref). Given a general consistent initial estimator, Assumptions (ref)--(ref) and (ref) are required to establish the theoretical guarantees for the associated iterative estimator introduced in Section (ref).
As pointed out in chang2023modelling,chang2021central, Assumptions (ref) and (ref) 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) focuses on the cases where all the factors, idiosyncratic errors, and the linear combination $\xi_t$ have exponentially decaying tails. Assumption (ref) is a standard $\alpha$-mixing condition allowing the data to be serially dependent but not necessarily stationary. As shown in (ref), 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) can accommodate both sparse and dense loadings within a unified formulation. For example, if $s_j \ll d_j$, Assumption (ref) 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 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 han2024cp and chen2026estimation.
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). Specifically, by Theorem 7 of 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) 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) 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) 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) 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) 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 han2024cp.
To identify the $r$ eigenvectors of $\mathbf{K}_{1,2,j}$ defined in (ref) 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 chang2023modelling.
Assumption (ref) 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 han2024cp.
Let $D_n=\prod_{j=1}^{m}d_j$ and $S_{n}=\prod_{j=1}^{m}s_{j}$, and further define
Set the threshold level $\delta_1=C_* (n^{-1} \log D_n)^{1/2}$ in (ref) for some constant $C_*>0$. Theorem (ref) 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).
In Algorithm (ref), 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
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) 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
where $\tilde c^{-1}=1+2c_1^{-1}+c_2^{-1}$ for $c_1$ and $c_2$ in Assumptions (ref) and (ref). Theorem (ref) gives the convergence rate of the iterative estimator obtained by Algorithm (ref).
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) as the initial estimator of Algorithm (ref). As discussed below Assumption (ref), 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) implies that the convergence rate of the one-pass estimator satisfies this requirement automatically. The first part of condition (ref) 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 (ref) is similar to Assumption (ref), 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), 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 han2024cp, we have $\gamma_{\max}=O(n^{-1})$ and the convergence rate in Theorem (ref) 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).
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 (ref). With $\xi_{t,i}^{\textup{s}}$ defined above (ref), write
for $j \in [m]$ and any $i,\ell \in [r]$. Theorem (ref) provides a limiting representation for the iterative estimator $\hat\mathbf{a}_{i,j}$.
Condition (ref) is to ensure that the asymptotic variance is not degenerate. Condition (ref) 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 (ref) 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$.
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), which requires the error process $\{\mathcal{E}_t\}_{t \ge 1}$ in (ref) 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) in the supplementary material provides some further discussion for this. We also discuss in Section (ref) in the supplementary material that our theoretical results can be extended from the exponential-decay assumptions in Assumptions (ref) and (ref) to polynomial-decay conditions. Assumption (ref) 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.
\begingroup {2pt} \linespread{0.9}\selectfont \endgroup