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.
104,910 characters · 18 sections · 85 citation commands
Supplementary Material to “Tensor Factor Model Estimation by Iterative Projection”
Motivated by a diverse range of modern scientific applications, analysis of tensors, or multi-dimensional arrays, has emerged as one of the most important and active research areas in statistics, computer science, and machine learning. Large tensors are encountered in genomics alter2005, omberg2007, neuroimaging analysis zhou2013, sun2017, recommender systems bi2018, computer vision liu2012, community detection anandkumar2014, among others. High-order tensors often bring about high dimensionality and impose significant computational challenges. For example, functional MRI produces a time series of 3-dimensional brain images, typically consisting of hundreds of thousands of voxels observed over time. Previous work has developed various tensor-based methods for independent and identically distributed (i.i.d.)\,tensor data or tensor data with i.i.d.\,noise. However, the statistical framework for general tensor time series data is much less studied in the literature.
Factor analysis is one of the most useful tools for understanding common dependence among multi-dimensional outputs. Over the past decades, vector factor models have been extensively studied in the statistics and economics communities. For instance, chamberlain1983, bai2002, stock2002 and bai2003 developed the static factor model using principal component analysis (PCA). They assumed that the common factors must have impact on most of the time series, and weak serial dependence is allowed for the idiosyncratic noise process. fan2011, fan2013, fan2018 established large covariance matrix estimation based on the static factor model. The static factor model has been further extended to the dynamic factor model in forni2000. In the dynamic factor model, the latent factors are assumed to follow a time series process, which is commonly taken to be a vector autoregressive process. fan2016 studied semi-parametric factor models through projected principal component analysis. pena1987identifying, pan2008, lam2011 and lam2012 adopted another type of factor model. They assumed that the latent factors {capture} all dynamics of the observed process, and thus the idiosyncratic noise process has no serial dependence. We will adopt this approach. We note that the factor process may have complex dynamic behavior, resulting in complex dynamics of the observed tensor, even with white additive noise process. Of course, when all the dynamics of the observed tensor process are `forced' to be included in the signal process induced by the factor process, situations may arise in which some factors are `weak' (or have impact on a small portion of the observed series in the tensor). {This} leads us to consider the `signal strength' in our investigation.
Although there have been significant efforts in developing methodologies and theories for vector factor models, there is a paucity of literature on matrix- or tensor-valued time series. wang2019 proposed a matrix factor model for matrix-valued time series, which explores the matrix structure. chen2020constrained established a general framework for incorporating domain and prior knowledge in the matrix factor model through linear constraints. chen2022modeling applied the matrix factor model to the dynamic transport network. chen2023statistical developed an inferential theory of the matrix factor model under a different setting from that in wang2019. chang2023modelling,han2023tensor,han2021cp studied factor models with CP type low rank structures.
Recently, chen2022factor introduced a factor approach for analyzing high dimensional dynamic tensor time series in the form
where ${\cal X}_1,...,{\cal X}_T\in\mathbb{R}^{d_1\times\cdots\times d_K}$ are the observed tensor time series, ${\cal M}_t$ and ${\cal E}_t$ are the corresponding signal and noise components of ${\cal X}_t$, respectively. The goal is to estimate the unknown signal tensor ${\cal M}_t$ from the tensor time series data. Following lam2012, it is assumed that the signal tensor accommodates all dynamics, making the idiosyncratic noise ${\cal E}_t$ uncorrelated (white) across time. It is further assumed that ${\cal M}_t$ lives in a lower dimensional space and has certain multilinear decomposition. Specifically, we assume that ${\cal M}_t$ satisfies a Tucker-type decomposition and model (ref) can be written as
where $A_k$ is the deterministic loading matrix of size $d_k\times r_k$ and $r_k\ll d_k$, and the core tensor ${\cal F}_t$ itself is a latent tensor factor process of dimension $r_1\times\ldots \times r_K$. Here the $k$-mode product of ${\cal X}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}$ with a matrix $U\in\mathbb{R}^{d_k'\times d_k}$, denoted as ${\cal X}\times_k U$, is an order $K$-tensor of size $d_1\times \cdots \times d_{k-1} \times d_k'\times d_{k+1}\times \cdots \times d_K$ such that $$ ({\cal X}\times_k U)_{i_1,...,i_{k-1},j,i_{k+1},...,i_K}=\sum_{i_k=1}^{d_k} {\cal X}_{i_1,i_2,...,i_K} U_{j,i_k}.$$
The core tensor ${\cal F}_t$ is usually much smaller than ${\cal X}_t$ in dimension. This structure provides an effective dimension reduction, as all the comovements of individual time series in ${\cal X}_t$ are driven by ${\cal F}_t$. Without loss of generality, assume that $A_k$ is of rank $r_k \ll d_k$. It should be noted that vector and matrix factor models can be viewed as special cases of our model since a vector time series is a tensor time series composed of a single fiber ($K=1$), and a matrix times series is one composed of a single slice ($K=2$).
chen2022factor proposed two estimation procedures, namely TOPUP and TIPUP, for estimating the column space spanned by the loading matrix $A_k$, for $k=1,\ldots,K$. The two procedures are based on different auto-cross-product operations of the observed tensors ${\cal X}_t$ to accumulate information, but they both utilize the assumption that the noise ${\cal E}_t$ and ${\cal E}_{t-h}, ~h>0$ are uncorrelated. The convergence rates of their estimators critically depend on $d=d_1d_2\ldots d_K$, a potentially very large number as $d_k,~ k=1,\ldots,K$, are large. Often a large $T$, the length of the time series, is required for accurate estimation of the loading spaces.
In this paper we propose extensions of the TOPUP and TIPUP procedures, motivated by the following observation. Suppose that the loading matrices $A_k$ are orthonormal with $A_k^\top A_k=I$, and we are given $A_2,\ldots, A_K$. Let \[ {\cal Z}_t={\cal X}_t\times_2 A_2^\top\times_3 \ldots \times_K A_K^\top; \mbox{\ \ and \ \ } {\cal E}_t^*={\cal E}_t\times_2 A_2^\top\times_3 \ldots \times_K A_K^\top; \] Then ((ref)) leads to
where ${\cal Z}_t$ is a $d_1\times r_2\times \ldots \times r_{K}$ tensor. Since $r_k\ll d_k$, ${\cal Z}_t$ is a much smaller tensor than ${\cal X}_t$. Under proper conditions on the combined noise tensor ${\cal E}_t^*$, the estimation of the loading space of $A_1$ based on ${\cal Z}_t$ can be made significantly more accurate, as the convergence rate now depends on $d_1r_2\ldots r_{K}$ rather than $d_1d_2\ldots d_{K}$.
Of course, in practice we do not know $A_2,\ldots, A_K$. Similar to backfitting algorithms, we propose an iterative algorithm. With a proper initial value, we iteratively estimate the loading space of $A_k$ at iteration $j$ based on \[ {\cal Z}_{t,k}^{(j)}={\cal X}_t\times_1\widehat{A}_1^{(j)\top}\times_2 \ldots \times_{k-1} \widehat{A}_{k-1}^{(j)\top}\times_{k+1} \widehat{A}_{k+1}^{(j-1)\top} \times_{k+2} \ldots \times_K \widehat{A}_{K}^{(j-1)\top}, \] using the estimate $\widehat{A}_{k'}^{(j-1)},~ k<k'\le K$ obtained in the previous iteration and the estimate $\widehat{A}_{k'}^{(j)},~ 1\le k'< k$, obtained in the current iteration. Our theoretical investigation shows that the iterative procedures for estimating $A_1$ can achieve the convergence rate as if all $A_2,\ldots,A_K$ are known and we indeed observe ${\cal Z}_t$ that follows model ((ref)). We call the procedure iTOPUP and iTIPUP, based on the matrix unfolding mechanism used, corresponding to TOPUP and TIPUP procedures. To be more specific, our algorithms have two steps: (i) We first use the estimated column space of factor loading matrices of TOPUP (resp. TIPUP) to construct the initial estimate of factor loading spaces; (ii) We then iteratively perform matrix unfolding of the auto-cross-moments of much smaller tensors ${\cal Z}_{t,k}^{(j)}$ to obtain the final estimator.
We note that the iterative procedure is related to higher order orthogonal iteration (HOOI) that has been widely studied in the literature; see, e.g., de2000, sheehan2007, liu2014, zhang2018tensor, among others. However, most of the existing works are not designed for tensor time series. They do not consider the special role of the time mode nor the covariance structure in the time direction. Typically HOOI treats the signal part as fixed or deterministic. In this paper we treat the signal as dynamic in the sense that the core tensor ${\cal F}_t$ in ((ref)) is dynamic and the relationship between ${\cal F}_t$ and the lagged ${\cal F}_{t-h}$ is of interest. Our setting requires special treatment although each iteration of our iterative procedures also consists of power up and orthogonal projection operations. While HOOI applies the SVD directly to the matrix unfolding of the iteratively projected data, in our approach the SVD is applied to the matrix unfolding of the outer- and inner-auto-cross-product of the iteratively projected data, respectively in iTOPUP and iTIPUP. Although the iTOPUP algorithm proposed here can be reformulated as a twist of HOOI on the auto-cross-moment tensor, the iTIPUP algorithm is different and cannot be recast equivalently as HOOI. More importantly, the theoretical analysis and theoretical properties of the estimators are fundamentally different from those of HOOI, due to the dynamic structure of tensor time series and the need to use the auto-cross-product operation between the SVD and data projection in each iteration. Different concentration inequalities are derived to study the performance bounds.
In this paper, we establish upper bounds on the estimation errors for both the iTOPUP and the iTIPUP, which are much sharper than the respective theoretical guarantees for TOPUP and TIPUP, demonstrating the benefits of using iterative projection. It is also shown that the number of iterations needed for convergence is of order no greater than $\log(d)$. We mainly focus on the cases where the tensor dimensions are large and of similar order. We also cover the cases where the ranks of the tensor factor process increase with the dimensions of the tensor time series.
chen2022factor showed that the TIPUP has a faster convergence rate in estimation error than the TOPUP, under a mild condition on the level of signal cancellation. In contrast, the theoretically guaranteed rate of convergence for the iTOPUP in this paper is of the same order or even faster than that for the iTIPUP under certain regularity conditions. Our results also suggest an interesting phenomenon. Using the iterative procedures, we find that the increase in either dimension or sample size can improve the estimation of the factor loading space of the tensor factor model with the tensor order $K\ge 2$. We believe that such a super convergence rate is new in the literature. Specifically, under proper regularity conditions, the convergence rate of the iterative procedures for estimating the space of $A_k$ is $O_{\mathbb{P}}(T^{-1/2}d_{-k}^{-1/2})$, where $d_{-k}=\prod_{j\ne k}d_j$, while the existing rate for non-iterative procedures is $O_{\mathbb{P}}(T^{-1/2})$ for the vector factor model lam2011 and the matrix/tensor factor models wang2019,chen2022factor. While the increase in the dimensions $d_k$ ($k=1,\ldots,K$) does not improve the performance of the non-iterative estimators, it significantly improves that of the proposed iterative estimators.
In addition, we establish the computational lower bound for the estimation of the loading spaces of tensor factor models under the hardness assumption of certain instances of hypergraphic planted clique detection problem. It shows that the sample size requirement (or signal to noise ratio condition) needed for using the TIPUP estimate as the initial values for the iterative procedures is unavoidable for any computationally manageable estimation procedure to achieve consistency, although the iterative procedures have faster convergence rates. Furthermore, we provide a statistical lower bound that matches the convergence rates of our iterative procedures under certain conditions, revealing a different effect of the ranks $r_k$ ($k=1,...,K$) compared to tensor Tucker decomposition zhang2018tensor.
{\it Related work.} We close this section by highlighting several recent papers on related topics. First, we draw attention to the work of foster1996time, fan2016 and chen2024semi. chen2024semi adopts an estimation precedure composed of a spectral initialization followed by an iterative refinement step, so that our methods are related to theirs. However, due to the differences in problem setting and model assumptions, their estimation procedures, performance bounds and analytic techniques are all significantly different from ours. foster1996time, fan2016 use the projection to the space spanned by the sieve bases without iteration. rogers2013multilinear assumes the tensor factor model in (ref), with an additional specific AR structure on the dynamic of the factor process. The additional model structure in their paper led to an EM type of estimation approach, quite different from the approach we develop here. wang2024high concerns low rank tensor AR model and uses a nuclear norm penalty to enforce the low rank structure and optimization algorithms for estimation, again quite different from our approach.
The paper is organized as follows. Section (ref) introduces basic notation and preliminaries of tensor analysis. We present the tensor factor model and the iTOPUP and iTIPUP procedures in Sections (ref) and (ref). Theoretical properties of the iTOPUP and iTIPUP are investigated in Section (ref). Section (ref) provides a brief summary. Numerical comparison of our iterative procedures and other methods, and all technical details are relegated to the Supplementary Material.
Throughout this paper, for a vector $x=(x_1,...,x_p)^\top$, define $\|x\|_q = (x_1^q+...+x_p^q)^{1/q}$, $q\ge 1$. For a matrix $A = (a_{ij})\in \mathbb{R}^{m\times n}$, write the SVD as $A=U\Sigma V^\top$, where $\Sigma=\text{diag}(\sigma_1(A), \sigma_2(A), ..., \sigma_{\min\{m,n\}}(A))$, with singular values $\sigma_{\max}(A) = $ $\sigma_1(A)\ge\sigma_2(A)\ge \cdots\ge \sigma_{\min\{m,n\}}(A)\ge 0$ in descending order. The matrix spectral norm is denoted as $\|A\|_{\rm S} =\sigma_1(A).$ Write $\sigma_{\min}(A)$ the smallest nontrivial singular value of $A$. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$, write $a_n=O(b_n)$ (resp. $a_n\asymp b_n$) if there exists a constant $C$ such that $|a_n|\leq C |b_n|$ (resp. $1/C \leq a_n/b_n\leq C$) for all sufficiently large $n$, and write $a_n=o(b_n)$ if $\lim_{n\to\infty} a_n/b_n =0$. Write $a_n\lesssim b_n$ (resp. $a_n\gtrsim b_n$) if there exist a constant $C$ such that $a_n\le Cb_n$ (resp. $a_n\ge Cb_n$). Write $a\wedge b=\min\{a,b\}$ and $a\vee b=\max\{a,b\}$. We use $C, C_1,c,c_1,...$ to denote generic constants, whose actual values may vary from line to line.
For any two $m\times r$ matrices with orthonormal columns, say, $U$ and $\widehat U$, suppose the singular values of $U^\top \widehat U$ are $\sigma_1\ge \sigma_2 \ge \cdots \ge \sigma_r\ge 0$. A natural measure of distance between the column spaces of $U$ and $\widehat U$ is then
which equals to the sine of the largest principle angle between the column spaces of $U$ and $\widehat U$. For any two matrices $A\in\mathbb{R}^{m_1\times r_1},B\in \mathbb{R}^{m_2\times r_2}$, denote the Kronecker product $\odot$ as $A\odot B\in \mathbb{R}^{m_1 m_2 \times r_1 r_2}$. For any two tensors ${\cal A}\in\mathbb{R}^{m_1\times m_2\times \cdots \times m_K}, {\cal B}\in \mathbb{R}^{r_1\times r_2\times \cdots \times r_N}$, denote the tensor product $\otimes$ as ${\cal A}\otimes {\cal B}\in \mathbb{R}^{m_1\times \cdots \times m_K \times r_1\times \cdots \times r_N}$, such that $$({\cal A}\otimes{\cal B})_{i_1,...,i_K,j_1,...,j_N}=({\cal A})_{i_1,...,i_K}({\cal B})_{j_1,...,j_N} .$$ Let ${\rm{vec}}(\cdot)$ be the vectorization of matrices and tensors. The mode-$k$ unfolding (or matricization) is defined as ${\rm{mat}}_k({\cal A})$, which maps a tensor ${\cal A}$ to a matrix ${\rm{mat}}_k({\cal A})\in\mathbb{R}^{m_k\times m_{-k}}$ where $m_{-k}=\prod_{j\neq k}^K m_j$. For example, if ${\cal A}\in\mathbb{R}^{m_1\times m_2\times m_3}$, then $$({\rm{mat}}_1({\cal A}))_{i,(j+m_2(k-1))}= ({\rm{mat}}_2({\cal A}))_{j,(k+m_3(i-1))}= ({\rm{mat}}_3({\cal A}))_{k,(i+m_1(j-1))} ={\cal A}_{ijk}. $$ For tensor ${\cal A}\in\mathbb{R}^{m_1\times m_2\times \cdots \times m_K}$, the Hilbert Schmidt norm is defined as $$ \|{\cal A}\|_{{\rm HS}}=\sqrt{\sum_{i_1=1}^{m_1}\cdots\sum_{i_K=1}^{m_K}({\cal A})_{i_1,...,i_K}^2 }. $$ For a matrix, the Hilbert Schmidt norm is just the Frobenius norm. Define the tensor operator norm for an order-4 tensor ${\cal A}\in\mathbb{R}^{m_1\times m_2\times m_3\times m_4}$, $$ \| {\cal A}\|_{\rm{op}} =\max\left\{ \sum_{i_1,i_2,i_3,i_4} u_{i_1,i_2} \cdot u_{i_3,i_4}\cdot ({\cal A})_{i_1,i_2,i_3,i_4}:\|U_1\|_{\rm HS}=\|U_2\|_{\rm HS}=1 \right\},$$ where {$U_1=(u_{i_1,i_2})\in\mathbb{R}^{m_1\times m_2}$ and $U_2=(u_{i_3,i_4})\in\mathbb{R}^{m_3\times m_4}$.}
Again, we consider {as in (ref)} \[ {\cal X}_t={\cal F}_t\times_1 A_1\times_2\ldots\times_K A_K+{\cal E}_t. \] Without loss of generality, assume that $A_k$ is of rank $r_k$. $A_k$ is not necessarily orthonormal, which is different from the classical Tucker decomposition tucker1966. Model (ref) is unchanged if we replace $(A_1,...,A_K, {\cal F}_t)$ by $(A_1H_1,...,A_KH_K, {\cal F}_t\times_{k=1}^K H_k^{-1})$ for any invertible $r_k\times r_k$ matrix $H_k$. Although $(A_1,...,A_K, {\cal F}_t)$ are not uniquely determined, the factor loading space, that is, the linear space spanned by the columns of $A_k$, is uniquely defined. Denote the orthogonal projection to the column space of $A_k$ as
{where $U_k$ is the left singular matrix in the SVD $A_k=U_k\Lambda_k V_k^\top$.} We use $P_k$ to represent the factor loading space of $A_k$. Thus, our objective is to estimate $P_k$.
The canonical representation of the tensor times series (ref) is written as $$ {\cal X}_t={\cal F}_t^{(\text{\footnotesize cano})} \times_{k=1}^K U_k+{\cal E}_t, $$ where the diagonal and right singular matrices of $A_k$ are absorbed into the canonical core tensor ${\cal F}_t^{(\text{\footnotesize cano})} = {\cal F}_t\times_{k=1}^K (\Lambda_k V_k^\top)$. In this canonical form, the loading matrices $U_k$ are identifiable up to a rotation in general and up to a permutation and sign changes of the columns of $U_k$ when the singular values are all distinct in the population version of the TOPUP or TIPUP methods, as we describe in Section (ref) below. In what follows, we may identify the tensor time series in its canonical form, i.e. $A_k=U_k$, without explicit declaration.
We do not impose any specific structure for the dynamics of the core tensor factor process ${\cal F}_t\in\mathbb{R}^{r_1\times \cdots\times r_K}$ beyond the independence between the core process and the noise process, and we do not require any additional structure on the correlation among different time series fibers of the noise process ${\cal E}_t$. Because of this generality, our estimator is based on the tensor version of the lagged sample cross product $\widehat \Sigma_h$, $h=1,...,h_0$, where
is an order-$2K$ tensor. The population version of this tensor autocovariance is
{Because} ${\cal M}_t={\cal M}_t\times_{k=1}^K P_k$ for all $t$,
with the notation $A_k=A_{k-K}$ and $P_k=P_{k-K}$ for all $k>K$.
In this paper, we consider iterative estimation procedures to achieve sharper convergence rates than the TOPUP and TIPUP procedures proposed in chen2022factor. We start with a quick description of their procedures as they serve as the starting point of our proposed iTOPUP and iTIPUP procedures. Note that the procedure in chen2022modeling and wang2019 is the non-iterative TOPUP.
{\bf (i) Time series Outer-Product Unfolding Procedure (TOPUP)}:
Let $\widehat \Sigma_h$ be the sample autocovariance of the data ${\cal X}_{1:T}=({\cal X}_1,\ldots,{\cal X}_T)$ as in (ref). Define
as a $d_k\times (dd_{-k}h_0)$ matrix, where $d=\prod_{k=1}^K d_k$, $d_{-k}=d/d_k$ and $h_0$ is a predetermined positive integer. Here we note that ${{\rm{TOPUP}}_k}$ is a function mapping a tensor time series to a matrix. In ${{\rm{TOPUP}}_k}$, the information from different time lags is accumulated, which is useful especially when the sample size $T$ is small. A relatively small $h_0$ is typically used, since the autocorrelation is often at its strongest with small time lags. See Remark (ref).
The TOPUP method performs SVD of (ref) to obtain the truncated left singular matrices
where LSVD$_m$ stands for the left singular matrix composed of the first $m$ left singular vectors corresponding to the largest $m$ singular values. Here $\text{$\widehat U_k$-TOPUP}$ is treated as an operator that maps a noisy tensor time series to a matrix of $m$ columns as an estimate of the mode-$k$ singular space of the low-rank signal tensor time series.
By (ref) and (ref), the expectation of (ref) satisfies
so that the TOPUP is expected to be consistent in estimating the column space of $A_k$.
{\bf (ii) Time series Inner-Product Unfolding Procedure (TIPUP)}:
Similar to (ref), define a $d_k\times (d_k h_0)$ matrix as
which replaces the tensor product by the inner product {through (ref) in} ((ref)). The TIPUP method performs SVD:
for $k=1,...,K$. Again, $\widehat U_k$-TIPUP is treated as an operator. We note that the TOPUP method in (ref) utilizes the entire auto-cross product tensor by applying the SVD to its mode $k$ unfolding, whereas the TIPUP only utilizes a matrix-valued linear mapping of the auto-cross product tensor by first taking the model-$k'$ trace operation for all $k'\neq k$. The trace operation cancels the noise but also possibly some signal.
{\bf (iii) iTOPUP and iTIPUP}: Next we describe a generic iterative procedure under the motivation described in Section 1. Its pseudo-code is provided in Algorithm (ref). It incorporates two estimators/operators $\widehat U_k$-INIT and $\widehat U_k$-ITER that map a tensor time series to an estimate of the loading matrix $U_k$. Respectively they stand for the procedures used for initialization and iteration. The $\widehat U_k$-TOPUP and $\widehat U_k$-TIPUP operators in (ref) and (ref) are examples of such operators.
When we use the $\widehat U_k$-TOPUP operator (ref) for both $\widehat U_k$-INIT and $\widehat U_k$-ITER in Algorithm (ref), it will be called iTOPUP procedure. Similarly, iTIPUP uses $\widehat U_k$-TIPUP operator (ref) for both $\widehat U_k$-INIT and $\widehat U_k$-ITER. Besides these two versions, we may also use \text{$\widehat U_k$-TIPUP} for \text{$\widehat U_k$-INIT} and \text{$\widehat U_k$-TOPUP} for \text{$\widehat U_k$-ITER}, named as TIPUP-iTOPUP. Similarly, TOPUP-iTIPUP uses \text{$\widehat U_k$-TOPUP} as \text{$\widehat U_k$-INIT} and \text{$\widehat U_k$-TIPUP} as \text{$\widehat U_k$-ITER}. These variants are sometimes useful, because TOPUP and TIPUP have different theoretical properties as the initializer or for iteration, as we will discuss in Section (ref). Other estimators of the loading spaces based on the tensor time series can also be used in place of \text{$\widehat U_k$-INIT} and \text{$\widehat U_k$-ITER}, such as the conventional high order SVD for tensor decomposition, which we refer to as Unfolding Procedure (UP), that simply performs SVD of the matricization along the appropriate mode of the $K+1$ order tensor $({\cal X}_1,\ldots, {\cal X}_T)$ with time dimension as the additional $(K+1)$-th mode.
In this section we present some theoretical properties of the iterative procedures. We first present the additional notation needed for the discussion, and then the error bounds for the iterative estimators under a minimum condition on the error process ${\cal E}_t$ in the model. These error bounds are quite general and cover many different models. To help decipher the general results, we present two concrete signal process models (or general sets of assumptions) with simpler and more explicit convergence rates.
Let $\overline\mathbb{E}[\cdot]=\mathbb{E}[\cdot|\{ {\cal F}_1,...,{\cal F}_T\}]$. Define $d=\prod_{k=1}^K d_k$, $d_{-k}=d/d_k$, $r=\prod_{k=1}^K r_k$ and $r_{-k}=r/r_k$. Define order-4 tensors
with $U_k$ from the SVD $A_k=U_k \Lambda_k V_k^\top$. We view $\Phi^{(\text{\footnotesize cano})}_{k,h}$ as the canonical version of the auto-covariance of the factor process. The noiseless version of the matrix TOPUP$_k$ in (ref) is
with $\Theta_{k,1:h_0} = (\Theta_{k,h}, h=1,\ldots,h_0)$. The canonical factor version of (ref) is $\text{mat}_1(\Phi^{{(\text{\footnotesize cano})}}_{k,1:h_0})\in\mathbb{R}^{r_k\times (rr_{-k}h_0)}$ with $\Phi^{{(\text{\footnotesize cano})}}_{k,1:h_0}=(\Phi^{{(\text{\footnotesize cano})}}_{k,h}, h=1,\ldots,h_0) \in \mathbb{R}^{r_k\times r_{-k}\times r_k \times r_{-k}\times h_0}$. Similarly define
The noiseless version of (ref) is
and its canonical factor version is $\Phi_{k,1:h_0}^{*{(\text{\footnotesize cano})}}=(\Phi_{k,h}^{*{(\text{\footnotesize cano})}}, h=1,\ldots,h_0) \in \mathbb{R}^{r_k\times (r_kh_0)}$. Let $\tau_{k,m}$ be the $m$-th singular value of the noiseless version of the TOPUP$_k$ matrix, $$\tau_{k,m}=\sigma_{m}\left(\,\overline\mathbb{E} \big[{{\rm{TOPUP}}_k}\big] \right) = \sigma_{m}\big({\text{mat}}_1(\Theta_{k,1:h_0})\big) = \sigma_{m}\big({\text{mat}}_1(\Phi^{(\text{\footnotesize cano})}_{k,1:h_0})\big). $$ The signal strength for iTOPUP can be characterized as
Similarly, let $$\tau_{k,m}^*=\sigma_{m}(\overline\mathbb{E} ({{\rm{TIPUP}}_k})) = \sigma_{m}\big( \Theta_{k,1:h_0}^*\big) = \sigma_{m}\big( \Phi_{k,1:h_0}^{*(\text{\footnotesize cano})}\big). $$ The signal strength for iTIPUP can be characterized as
We note that by (ref) and the Cauchy-Schwarz inequality,
Our general error bounds {for the proposed iTOPUP and iTIPUP} are established under the following assumption for the error process.
Assumption (ref) is used in chen2022factor for the theoretical investigation of the non-iterative TIPUP and TOPUP, and is similar to those on the noise imposed in lam2011, lam2012. The normality assumption, which ensures fast convergence rates in our analysis, is imposed for technical convenience. It accommodates general patterns of dependence among individual time series fibers, but also allows a presentation of the main results with manageable analytical complexity. In fact, direct extension is visible in our analysis under the sub-Gaussian and even more general tail probability conditions. Under Assumption (ref) the magnitude of the noise can be measured by the dimension $d_k$ before the projection and by the {rank} $r_k$ after the projection. The main theorems (Theorems (ref), (ref) and (ref)) in this section are based on this assumption on the noise alone, and cover all thereafter discussed settings of the signal ${\cal M}_t$.
Let us first study the behavior of iTOPUP procedure. By chen2022factor, the risk $\overline\mathbb{E}\big[\big\|\widehat U_k^{(0)}\widehat U_k^{(0)\top} - U_kU_k^\top\big\|_{\rm S}\big]$ of the TOPUP estimator for $U_k$, the initialization of iTOPUP, is no larger than a constant times
where $d_{-k} = \prod_{j\neq k}d_j$ and $r_{-k}=\prod_{j\neq k}r_j$. A variation of the wedin1972 perturbation theory, stated in Lemma (ref), provides a sharper bound for the TOPUP estimator as follows.
The aim of iTOPUP is to achieve dimension reduction by projecting the data in other modes of the tensor time series from $\mathbb{R}^{d_j}$ to $\mathbb{R}^{r_j}$, $j\neq k$. Ideally (e.g. when the true projection matrices $U_j$ are used), this would reduce the rate given in (ref) and (ref) to
by replacing all $d_j$ in $R_k^{(0)}$ with $r_j$, $j\neq k$, where
However, because the iteration uses the estimated $U_j, j\neq k,$ of total dimension $d^*_{-k} = \sum_{j\neq k} d_jr_j$, our analysis also involves the following additional error term,
The following theorem provides conditions under which the ideal rate is indeed achieved.
Now, let us consider the statistical performance of iTIPUP procedure. Again, by chen2022factor the TIPUP risk in the estimation of $P_k$ is bounded by
with $d_{-k}=\prod_{j\neq k}d_j$, and the aim of iTIPUP is to achieve the ideal rate
through dimension reduction, where $r_{-k}=\prod_{j\neq k}r_j$. As in the case of iTOPUP, our error bound for iTIPUP involves the additional error term
The following theorem, which allows the ranks $r_k$ to grow to infinity as well as $d_k$ when $T\to\infty$, provides sufficient conditions to guarantee the ideal convergence rate for iTIPUP.
We briefly discuss the conditions and conclusions of Theorem (ref) as the details are parallel to the remarks below Theorme (ref). By (ref), (ref) and the Cauchy-Schwarz inequality, $(1-h_0/T){\lambda}_k^{*2} \le \|\Theta_{k,0}^*\|_{\rm S}$, so that the first condition in (ref) guarantees a sufficiently small $R^{*(0)}$, which implies a sufficiently small error in the initialization of iTIPUP by (ref). The second condition in (ref) again has two terms respectively reflecting the ideal rate after dimension reduction by the true $U_{-k}=\odot_{j\neq k}U_j$ in the estimation of $U_k$ and the extra cost of estimating $U_{-k}$. The upper bound (ref) for the error of the $m$-step estimator is also comprised of two terms representing the cumulative iteration error and contracted initialization error. In Corollary (ref) below with fixed $r_k$, the smallest $\rho=C_{1,K}^{(\text{\footnotesize iter})}({R}^{*(\text{\footnotesize ideal})} + {R}^{*(\text{\footnotesize add})})$ is taken in (ref) to achieve (ref) in one iteration when ${R}_k^{*(\text{\footnotesize ideal})}$ dominates $R_k^{*(\text{\footnotesize add})}$. Moreover, Theorem (ref) allows diverging ranks $r_k$ and convergence rate $T^{-1/2}d_{-k}^{-1/2}$poly$(r_1,\ldots,r_K)$ under proper conditions as discussed in Remark (ref).
As discussed in Section (ref), we can mix the TOPUP and TIPUP operations for the initiation and iterative operations in Algorithm (ref). For example, the proof of Theorems (ref) yields the following error bound for the mixed TIPUP-iTOPUP algorithm.
We omit the statement of an analogous error bound for the \hbox{\rm TOPUP-iTIPUP} algorithm.
In this section we provide the convergence rate when the dimensions of the factors ${\cal F}_t$, or equivalently the ranks of the signal process ${\cal M}_t$, $r_1,\ldots,r_K$, are fixed, and the auto-cross-outer-product of the factor process is ergodic. Formally, we impose the following additional assumption.
Under Assumption (ref), the factor process has a fixed expected auto-cross-moment tensor with fixed dimensions. The assumption that the condition numbers of $A_k^\top A_k$ ($k=1,...,K$) are bounded corresponds to the pervasive condition (e.g., stock2002, bai2003). It ensures that all the singular values of $A_k$ are of the same order. Such conditions are commonly imposed in factor analysis.
As our methods are based on auto-cross-moment at nonzero lags, we do not need to assume any specific model for the latent process ${\cal F}_t$, except some rank conditions in Assumption (ref)(i) and (ii). Since the columns of $\Phi_{k,1:h_0}^{*(\text{\footnotesize cano})}$ are linear combinations of those of $\text{mat}_1(\Phi_{k,1:h_0}^{(\text{\footnotesize cano})})$ and $\mathbb{E}[\text{mat}_1(\Phi_{k,1:h_0})]$ and $\mathbb{E}[\text{mat}_1(\Phi^{(\text{\footnotesize cano})}_{k,1:h_0})]$ have the same rank, Assumption (ref)(ii) implies Assumption (ref)(i).
In order to provide a more concrete understanding of Assumption (ref)(i) and (ii), consider the case of $k=1$ and $K=2$. We write the factor process ${\cal F}_t = (f_{i,j,t})_{r_1\times r_2}$, and the stationary auto-cross-moments $\phi_{i_1,j_1,i_2,j_2,h} = \mathbb{E} (f_{i_1,j_1,t-h}f_{i_2,j_2,t})$. Hence $\mathbb{E}[\text{mat}_1(\Phi_{k,1:h_0})]$ is a $r_k\times(r_{-k}r_kr_{-k}h_0)$ matrix, with columns being $\phi_{\cdot, j_1,i_2,j_2,h}$. Since $\mathbb{E}[\text{mat}_1(\Phi_{k,1:h_0})]\mathbb{E}[\text{mat}_1(\Phi_{k,1:h_0})]^{\top}$ is a sum of many semi-positive definite $r_k\times r_k$ matrices, if any one of these matrices is full rank, then $\mathbb{E}[\text{mat}_1(\Phi_{k,1:h_0})]$ is of rank $r_k$. Hence Assumption (ref)(i) is relatively easy to fulfill. On the other hand, Assumption (ref)(ii) is quite different. First, the condition is imposed on the canonical form of the model as the inner product in TIPUP related procedures behaves differently. Let ${\cal F}_t^{(\text{\footnotesize cano})} = U_1^\top {\cal M}_t U_2 = (f_{i,j,t}^{(\text{\footnotesize cano})})_{r_1\times r_2}$, and $\phi^{(\text{\footnotesize cano})}_{i_1,j_1,i_2,j_2,h} = \mathbb{E} (f_{i_1,j_1,t-h}^{(\text{\footnotesize cano})}f_{i_2,j_2,t}^{(\text{\footnotesize cano})}$). Then $\|\Phi^{*{(\text{\footnotesize cano})}}_{1,1:h_0}\|_{\rm HS}^2 =\sum_{h=1}^{h_0}\sum_{i_1,i_2}\big(\sum_{j=1}^{r_2} \phi_{i_1,j,i_2,j,h}^{{(\text{\footnotesize cano})}}\big)^2$. As $\phi_{i_1,j,i_2,j,h}^{{(\text{\footnotesize cano})}}$ may be positive or negative for different $i_1,i_2,j,h$, the summation $\sum_{j=1}^{r_2} \phi_{i_1,j,i_2,j,h}^{{(\text{\footnotesize cano})}}$ is subject to potential signal cancellation for $h>0$. Assumption (ref)(ii) ensures that there is no complete signal cancellation that makes the rank of $\mathbb{E}[\Phi_{k,1:h_0}^{*{(\text{\footnotesize cano})}}]$ less than $r_k$. While the signal cancellation rarely causes the rank deficiency, the resulting loss of efficiency may still have an impact on the finite sample performance as our simulation results demonstrate. Of course complete signal cancellation is less likely with larger $h_0$.
The following corollary is a simplified version of Theorem (ref) under Assumption (ref)(i).
Corollary (ref) asserts that, in order to recover the factor loading space for $A_k$, the signal to noise ratio needs to satisfy $\lambda/\sigma \ge C_{0,K,r} \max_{k\le K}(d^{1/2}T^{-1/2}+d^{1/2}d_{k}^{-1/4}T^{-1/4})$ as in (ref), and the ideal rate (ref) can be achieved in one iteration. Under Assumptions (ref), (ref)(i), the error bound in Proposition (ref) yields the convergence rate
In comparison, the ideal rate is much sharper than the convergence rate of the non-iterative TOPUP in Proposition (ref) when $\lambda^2/\sigma^2\ll \min_{k\le K}\{d^{4/3}/(T^{1/3}d_{k}),d^2/(T^{1/2}d_{k}^{3/2})\}$.
The following corollary is a simplified version of Theorem (ref) under Assumption (ref)(ii), which excludes severe signal cancellation in iTIPUP.
Compared with the results in Corollary (ref) for iTOPUP, the achieved ideal rate (ref) is the same. However, the signal-to-noise ratio requirement (ref) is weaker but Assumption (ref)(ii) is stronger in Corollary (ref) for iTIPUP. Again, the ideal rate is much sharper than the convergence rate of the non-iterative TIPUP in chen2022factor.
The main theorems in Subsection 3.2 allow for the case where the dimensions of the core factor, $r_1,...,r_K$, diverge as the dimensions of the observed tensor $d_1,...,d_K$ grow to infinity. The following assumption provides a concrete set of conditions that can be used to provide some insights of the properties of iTOPUP and iTIPUP in such scenarios.
Assumption (ref) is similar to the signal strength condition of lam2012, and the pervasive condition on the factor loadings (e.g., stock2002 and bai2003). It is more general than Assumption (ref) in the sense that it allows $r_1,...,r_K$ to diverge and the latent process ${\cal F}_t$ does not have to be weakly stationary.
We take $\delta_0,\delta_1$ as measures of the strength of the signal process ${\cal M}_t$. They roughly indicate how much information is contained in the signals compared with the amount of noise, with respect to the dimensions and ranks, $d,r$ and $r_k$. In this sense, they reflect the signal to noise ratio. When $\delta_0=\delta_1=0$, the factors are called strong factors; otherwise, the factors are called weak factors.
We describe below the convergence rate of iTOPUP in terms of $d_k$, $r_k$ and $T$ under Assumption (ref)(i) when the dimensions of the core factor $r_1,...,r_K$ are allowed to diverge.
Note that the second part of Corollary (ref) says that when the condition is right, iTOPUP algorithm only needs a small number of iterations to converge, as $O(\log r)$ is typically very small. The noise level $\sigma$ does not appear directly in the rate since it is incorporated in the signal to noise ratio in the tensor form in Assumption (ref). In Corollary (ref), we show that as long as the sample size $T$ satisfies (ref), the iTOPUP achieves consistent estimation under proper regularity conditions. To digest the condition, we notice that (ref) becomes $T\ge C_0 \max_k (r_k r_{-k}^2) $ when the growth rate of $r_k$ is much slower than $d_k$ and the factors are strong with $\delta_0=\delta_1=0$.
The advantage of using index $\delta_0,\delta_1$ is to link the convergence rates of the estimated factor loading space explicitly to the strength of factors. It is clear that the stronger the factors are, the faster the convergence rate is. Equivalently, the stronger the factors are, the smaller the sample size is required.
When the ranks $r_k$ ($k=1,...,K$) also {diverge} and there is no severe signal cancellation in iTIPUP, we have the following convergence rate for iTIPUP under Assumption (ref)(ii).
When the average auto-correlation is of unit order and the signal cancellation for TIPUP has no impact on the order of the signal ($\delta_0=\delta_1$ and $\delta_2=0$ respectively), Corollary (ref) requires the sampling rate $T\gtrsim h_0+(d_kr_k+d^{\delta_0}r_k^2+d_{-k}^*r_k(1+r/d^{1-\delta_0}))/d^{1-\delta_0}$ and provides the convergence rate $(r_kd_k)^{1/2}(1+r/d^{1-\delta_0})^{1/2}/(Td^{1-\delta_0})^{1/2}$. For examples, $T\ge 4h_0+C_1$ gives the rate $(r_kd_k)^{1/2}/(Td^{1-\delta_0})^{1/2}$ when $\delta_0\le (K-2)/(2K)$ and $r_k^2\lesssim d_k\asymp d^{1/K}\,\forall k$, and the sample size requirement can be written as $T\gtrsim h_0+d^{\delta_0}r_k^2/d^{1-\delta_0}$ when $r_k^2\asymp r^{2/K}\lesssim d_k\asymp d^{1/K}\,\forall k$ regardless of $\delta_0\in [1/K,1]$. Thus, the side condition involving $R^{*(\text{\footnotesize add})}$ in the second part of (ref) is absorbed into the other components of (ref).
Corollary (ref) and Corollary (ref) offer comparison of the iTOPUP and iTIPUP when the ranks diverge from two perspectives: sample size requirements and convergence rates. The lower bounds on $T$ in (ref) in Corollary (ref) and (ref) in Corollary (ref) provide the sample complexity of the iTOPUP and iTIPUP respectively. In the case that the growth rate of $r_k$ is much slower than $d_k$ and the factors are strong with $\delta_0=\delta_1=0$, the required sample size of the iTIPUP reduces to $T\ge 4h_0+C_0\max_{j,k} \left( r_kr_{-k}^{2\delta_2} r_{-j}^{2\delta_2}/d_{-k}+ r_kr_{-k}^{2\delta_2} r_j/d_{-j}\right)$, where $r_{-k}=r/r_k$ and $d_{-k}=d/d_k$. By comparing with the comment after Corollary (ref), where the sample size requirement for the iTOPUP is $T\ge C_0\max_k(r_kr_{-k}^2)$ when $\delta_0=\delta_1=0$, it can be seen that the sample complexity for the iTIPUP is smaller, if $\delta_2$ is a small constant. From the perspective of convergence rate, let us compare (ref) in Corollary (ref) and (ref) in Corollary (ref). When ranks diverge, iTIPUP is slower than iTOPUP if $\delta_2>3/2$, or $\{0\le\delta\le 3/2, d_k\gtrsim rr_{-k}^{2-2\delta_2}, d_{-k}\gtrsim r r_{-k}^{3-2\delta_2}\}$, and faster if $d_k\lesssim r_kr_{-k}^{2-2\delta_2}$, no matter how strong the factor is or what values $\delta_0,\delta_1$ take. As expected, the convergence rate is slower in the presence of weak factors. See the simulation for more empirical evidence.
Similar to Corollaries (ref) and (ref), we have the following rate for TIPUP-iTOPUP.
Compared with Corollary (ref), Corollary (ref) provides the same error bound for smaller $T$ (possibly with bounded $T\gtrsim h_0$) when $r_{-k}^{2\delta_2} \lesssim r_{-k}d_{-k}$. The side condition involving $R^{(\text{\footnotesize add})}$ in the second part of (ref), corresponding to the last component of (ref) involving $d_{-k}^*$, is absorbed into the other components of (ref) when $r_k^{1/2}\le d^{\delta_1-\delta_0}\big(r_{-k}^{2\delta_2-1} + r_{-k}^2\big)\, \forall k\le K$.
Theorems (ref) and (ref) show that the convergence rates of the non-iterative estimators TOPUP and TIPUP can be improved by their iterative counterparts. Particularly, when the dimensions $r_k$ for the factor process are fixed and the respective signal strength conditions are fulfilled, the proposed iTOPUP and iTIPUP just need one-iteration to achieve the much sharper ideal rate $R^{(\text{\footnotesize ideal})}$ in (ref) and $R^{*(\text{\footnotesize ideal})}$ in (ref), compared with the rate (ref) of TOPUP and (ref) of TIPUP derived in chen2022factor, respectively. The improvement is achieved through replacing the much larger $d_{-k}$ by $r_{-k}$, via orthogonal projection. When the factors are strong with $\delta_0=\delta_1=0$ and the factor dimensions are fixed, the non-iterative TOPUP-based estimators of lam2011 for the vector factor model, wang2019 for the matrix factor and chen2022factor for tensor factor models all have the same $O_{\mathbb{P}}(T^{-1/2})$ convergence rate for estimating the loading space. In comparison, the convergence rate $O_{\mathbb{P}}(T^{-1/2}d_{-k}^{-1/2})$ of both iterative estimators, iTOPUP and iTIPUP (when there is no severe signal cancellation, with bounded $\delta_2$), is much sharper. Intuitively, when the signal is strong, the orthogonal projection operation helps to consolidate signals while potentially averaging out the noises, when the projection reduces {the dimension of} the mode-$k$ unfolded matrix from {$d_k\times d_{-k}$} for the tensor ${\cal X}_t$ {to $d_k\times r_{-k}$} for the projected tensor ${\cal Z}_t$, resulting in the improvement by a factor of $d_{-k}^{-1/2}$ in the convergence rate.
When $r_k$ are allowed to diverge, the iTOTUP and iTIPUP algorithms converge after at most $O(\log(d))$ iterations to achieve the ideal rate according to Theorems (ref) and (ref). The number of iterations needed can be as few as $O(\log(r))$ when the condition is right.
The inner product operation in (ref) for TIPUP-related procedures enjoys significant amount of noise cancellation comparing to the outer product operation in (ref) for TOPUP-related procedures. Compared with iTOPUP, the benefit of noise cancellation of the iTIPUP procedure is still visible through the reduction of $r_{-k}$ in (ref) to $\sqrt{r_{-k}}$ in (ref) in the ideal rates. However, this post-iteration benefit is much less pronounced compared with the reduction of $d_{-k}$ in (ref) for TOPUP to $\sqrt{d_{-k}}$ in (ref) for TIPUP in the non-iterative rates. Meanwhile, the potential {for} signal cancellation in the TIPUP related schemes persists as ${\lambda}^*_k$ and ${\lambda}_k$ are unchanged between the initial and ideal rates. We note that the signal strength can be viewed as ${\lambda}_k$ and ${\lambda}^*_k$ in Theorems (ref) and (ref) respectively for TOPUP/iTOPUP and TIPUP/iTIPUP, and that severe signal cancellation can be expressed as ${\lambda}_k^*\ll{\lambda}_k$. When $r_{-k}$ are allowed to diverge to infinity, the impact of signal cancellation is expressed in terms of $\delta_2$ in Assumption (ref): The iTOPUP has a faster rate than the iTIPUP when $\delta_2>3/2$, or $\{0\le\delta\le 3/2, d_k\gtrsim rr_{-k}^{2-2\delta_2}, d_{-k}\gtrsim r r_{-k}^{3-2\delta_2}\}$, and slower rate when $d_k\lesssim r_kr_{-k}^{2-2\delta_2}$, in view of Corollary (ref) and (ref). In Corollaries (ref) and (ref), iTOPUP and iTIPUP have the same convergence rate because Corollary (ref) assumes that signal cancellation does not change convergence rate.
Our results seem to suggest that the mixed TIPUP-iTOPUP procedure would strike a good balance between the benefit of noise cancellation (e.g. smaller $T$ for consistency) and the potential danger of signal cancellation (e.g. ${\lambda}^*_k\ll{\lambda}_k$) for the following four reasons: (1) The benefit of noise cancellation is much larger in the initialization, in term of $d_{-k}$, in view of the rates $R_k^{(0)}$ in (ref) and $R^{*(0)}$ in (ref). (2) The first part of condition (ref) for TIPUP-iTOPUP is weaker than the first part of condition (ref) for TIPUP-iTIPUP. (3) The signal strength ${\lambda}_k$ of the stronger TOPUP form is retained in the rate $R^{(\text{\footnotesize ideal})}$ after iTOPUP iteration. (4) As we will prove in Section (ref), the sample size requirement for the TIPUP initialization is optimal in the sense that it matches a computational lower bound under suitable conditions. Our simulation results support this recommendation, especially for relatively small $r_{-k}$. Of course if the sample size qualitatively justifies the condition $C_{1}^{(\text{\tiny TOPUP})}{R}^{(0)}\le (1-\rho)/4$ in (ref) and/or if a possible signal cancellation is a significant concern, the TOPUP initiation should be used.
The signal to noise ratio (SNR) condition, or equivalently the sample size requirement, is mainly used to ensure that the initial estimator has sufficiently small estimation error. Thus, the performance of iterative procedures is measured by both the SNR requirement and the error rate achieved. Consider fixed $h_0$ in the fixed rank case with $K=3$ and $d_{\max}\asymp d^{1/K}$. In the fixed signal model where ${\cal M}_t ={\cal M}$ is fixed and deterministic in (ref), applying HOOI to the average of ${\cal X}_t$ would require SNR $\lambda(T^{1/2}/\sigma)\ge C_0 d^{1/4}$ to achieve the loss of the order $(\sigma/T^{1/2}) d_k^{1/2}/{\lambda}$ according to zhang2018tensor, where $\sigma/T^{1/2}$ is viewed as the noise level for HOOI as it is the standard deviation of each element of the average tensor. In terms of the auto-crossproducts, taking the average over ${\cal X}_t$ roughly amounts to taking the average of all $T(T-1)/2$ lagged products between ${\cal X}_{t-h}$ and ${\cal X}_t$, $1\le t-h<t\le T$. However, in the tensor factor model (ref) where the signal part is random and serial correlated, the average is taken only over $T-h$ lagged products for each $h$. Thus, while the rate of the average of the signal-by-noise crossproducts in the factor model is heuristically expected to match that of HOOI at noise level $\sigma/T^{1/2}$, the rate of the average of the noise-by-noise crossproducts in the factor model is expected to only match that of HOOI with noise level $\sigma/T^{1/4}$. In Corollary (ref), the contribution of the noise-by-noise crossproducts dominates the initial estimation error as the SNR requirement $\lambda(T^{1/4}/\sigma)\ge C_0 d^{1/4}$ in (ref) matches that of HOOI with noise level $\sigma/T^{1/4}$; at the same time the contribution of the signal-by-noise crossproducts dominates the estimation error after iteration as the rate $(\sigma/T^{1/2}) d_k^{1/2}/{\lambda}$ in (ref) matches that of HOOI with noise level $\sigma/T^{1/2}$. Thus, if there is no severe signal cancellation, the signal to noise ratio requirement and convergence rate for iTIPUP and TIPIP-iTOPUP in the factor model are both comparable with those of HOOI in the simpler fixed signal setting, but the rate match is achieved in very different and subtle ways. We prove that this insight is intrinsic as the rates in (ref) and (ref) are both optimal according to the computational and statistical lower bounds in the following subsection.
In this subsection, we focus on the typical factor model setting that the condition numbers of $A_k^\top A_k$ are bounded. We shall prove that under the computational hardness assumption, the signal to noise ratio condition (ref) imposed on iTIPUP (also TIPUP-iTOPUP) in Corollary (ref) is unavoidable for computationally feasible estimators to be consistent. To be specific, we show that, if the signal to noise ratio condition is violated, then any computationally efficient and consistent estimator of the loading spaces leads to a computationally efficient and statistically consistent test for the Hypergraphic Planted Clique Detection problem in a regime where it is believed to be computationally intractable. In addition, we establish a statistical lower bound on the minimax risk of the estimators.
{\it Hypergraphic Planted Clique.} An $m$-hypergraph $G=(V(G),E(G))$ is a natural extension of regular graph, where $V(G)=[N]$ and each hyper-edge is represented by an unordered group of $m$ different vertices $i_j\in V(G)$ ($j=1,...,m$), denoted as $e=(i_1,...,i_m)\in E(G)$. Given a $m$-hypergraph its adjacency tensor ${\cal A}\in\{0,1\}^{N\times N\times\cdots\times N}$ is defined as
We denote by ${\cal G}_m(N, 1/2)$ the Erd\H{o}s–R\'enyi $m$-hypergraph on $N$ vertices where each hyper-edge $e$ is drawn independently with probability $1/2$, by ${\cal C}={\cal C}(N,\kappa)$ a random clique of size $\kappa$ where the $\kappa$ members are uniformly sampled from $[N]$ and $E({\cal C})$ is composed of all $e=(i_1,\ldots,i_m)$ with $i_j\in {\cal C}$, and by ${\cal G}_m(N, 1/2, \kappa)$ the random graph generated by first sampling independently ${\cal G}_m(N, 1/2)$ and ${\cal C}={\cal C}(N,\kappa)$ and then adding all the edges in $E({\cal C})$ to the set of edges in ${\cal G}_m(N,1/2)$. The Hypergraphic Planted Clique (HPC) detection problem of parameter $(N, \kappa, m)$ refers to testing the following hypotheses:
If $m=2$, the above HPC detection becomes the traditional planted clique (PC) detection problem. When $\kappa \ge c \sqrt{N}$, many computationally efficient algorithms have been developed for PC detection; see, alon1998finding, feige2000finding, feige2010finding, ames2011nuclear, dekel2014finding, deshpande2015finding, feldman2017statistical, among others. However, it has been widely conjectured that when $\kappa=o(\sqrt{N})$, the PC detection problem cannot be solved in randomized polynomial time, which is referred to as the hardness conjecture. Computational lower bounds in several statistical problems have been established by assuming the hardness conjecture of PC detection, including sparse PCA berthet2013complexity,berthet2013optimal, wang2016statistical, sparse CCA gao2017sparse, submatrix detection ma2015computational, cai2017computational, community detection hajek2015computational, etc.
Recently, motivated by tensor data analysis, hardness conjecture for HPC detection problem has been proposed; see, for example, zhang2018tensor, brennan2020reducibility, luo2022tensor, luo2020open, pananjady2022isotonic. Similar to the PC detection, they hypothesized that when $\kappa = O(N^{1/2-\delta})$ with $\delta>0$, the HPC detection problem (ref) cannot be solved by any randomized polynomial-time algorithm. Formally, the conjectured hardness of the HPC detection problem can be stated as follows.
Evidence supporting this hypothesis has been provided in zhang2018tensor,luo2022tensor. This version of the hypothesis is similar to the one in berthet2013complexity, ma2015computational, gao2017sparse for the PC detection problem.
For simplicity, we especially consider the factor model (ref) with each individual series of ${\cal F}_t$ being mean 0 and independent,
where $U_k\in\mathbb{R}^{d_k\times r_k}$, $U_k^\top U_k = I$ for $1\le k\le K$, and $0<c_1\le \sigma_{\min}(\mathbb{E} \mathrm{vec}1({\cal F}_t) \mathrm{vec}1^\top({\cal F}_t))\le \sigma_{\max}(\mathbb{E} \mathrm{vec}1({\cal F}_t) \mathrm{vec}1^\top({\cal F}_t))\le c_2<\infty$. The probability space we consider in this section is
The computational lower bound over ${\mathscr P}(T,d_1,...,d_K,\lambda)$ is then presented as below. In the case of $K=1$, the problem reduces to PCA of the spike covariance matrix. For general $K$, the auto-covariance tensor is of order $2K$.
Comparing (ref) with (ref), we see that the signal to noise ratio condition (ref) cannot be improved upon by a factor of $d^\vartheta$ with polynomial time complexity for any $\vartheta>0$. The condition $d_k \ge T$ is a technical requirement to use the theoretical tools in ma2015computational and brennan2020reducibility for the reduction from HPC.
Next, we establish the statistical lower bound for the tensor factor model problem. Again, we consider the probability space (ref).
In Lemma (ref) below, we provide an improvement of the matrix perturbation bound of wedin1972. The lemma, proved in Appendix (ref) in the supplementary material and used to prove Proposition (ref), is of independent interest due to wide applications of the wedin1972 bound.
The sharper perturbation bound in the middle of (ref) improves the commonly used version of the wedin1972 bound on the right-hand side, compared with Theorem 1 of cai2018 and Lemma 1 of chen2022rejoinder. As cai2018 pointed out, such variations of the wedin1972 bound provide sharper convergence rate when $\hbox{\rm error}_2\le \hbox{\rm error}_1$ in (ref), typically in the case of $d_1\ll d_2$, as in Proposition (ref).
In this paper we propose new estimation procedures for tensor factor model via iterative projection, and focus on two procedures: iTOPUP and iTIPUP. Theoretical analysis shows the asymptotic properties of the estimators. Simulation study presented in the supplementary material illustrates the finite sample properties of the estimators. While theoretical results are obtained under very general conditions, concrete specific cases are considered. In particular, under the typical factor model setting where the condition numbers of $A_k^\top A_k$ are bounded and the ranks $r_k$ are fixed, the proposed iterative procedures, iTOPUP method and iTIPUP method (with no severe signal cancellation) lead to a convergence rate $O_{\mathbb{P}}((Td_{-k})^{-1/2})$ under strong factors settings due to information pooling of the orthogonal projection of the other $d_{-k}$ dimensions. This rate is much sharper than the existing rate $O_{\mathbb{P}}(T^{-1/2})$ in the recent literature for non-iterative estimators for vector, matrix and tensor factor models. It implies that the accuracy can be improved by increasing the dimensions, and consistent estimation of the loading spaces can be achieved even with a fixed finite sample size $T$. This is in sharp contrast to the folklore based on the existing literature that only the sample size $T$ helps the estimation of the loading matrices in factor models. The proposed iterative estimation methods not only preserve the tensor structure, but also result in sharper convergence rate in the estimation of factor loading space.
The iterative procedure requires two operators, one for initialization and one for iteration. Under certain conditions of the signal to noise ratio (or the sample size requirement), we only need the initial estimator to have sufficiently small estimation errors but not the consistency of the initial estimator. Often, one iteration is sufficient. In more complicated general cases, at most $O(\log(d))$ iterations are needed to achieve the ideal rate of convergence. Based on the theoretical results and empirical evidence, we suggest to use iTOPUP for iteration when the ranks $r_k$ are small. In terms of initiation, the computational lower bound shows that the signal to noise ratio condition derived from TIPUP initialization is unavoidable for any computationally feasible estimation procedure to achieve consistency, while that from TOPUP initialization is not optimal. Based on this result, we suggest the use of TIPUP initialization. Of course, this should be done with precaution against potential signal cancellation, for example by using a slightly large $h_0$ as our empirical results show. By examination of the patterns of estimated singular values under different lag values $h_0$, using iTOPUP and iTIPUP, it is possible to detect signal cancellation, which has significant impact on iTIPUP estimators.
The proposed iterative procedure is similar to HOOI algorithms in spirit, but the detailed operations and the theoretical challenges are significantly different.
We would like to thank the Editor, the Associate Editor and the anonymous referees for their detailed reviews, which helped to improve the paper substantially.
Yuefeng Han's research is supported in part by National Science Foundation grant IIS-1741390. Rong Chen's research is supported in part by National Science Foundation grants DMS-1737857, IIS-1741390, CCF-1934924, DMS-2027855 and DMS-2319260. Dan Yang's research is supported in part by NSF grant IIS-1741390, Hong Kong grant GRF 17301620, Hong Kong grant CRF C7162-20GF and Shenzhen grant SZRI2023-TBRF-03. Cun-Hui Zhang's research is supported in part by NSF grants DMS-1721495, IIS-1741390, CCF-1934924, DMS-2052949 and DMS-2210850.