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.
118,154 characters · 10 sections · 66 citation commands
CP Factor Model for Dynamic Tensors
{\bf Keywords}: Tensor Factor Model, CP Decomposition, Tensor Time Series, Dimension Reduction, Orthogonal Projection.
In recent years, information technology has made tensors or high-order arrays observations routinely available in applications. For example, such data arises naturally from genomics alter2005, omberg2007, neuroimaging analysis zhou2013tensor, sun2017store, recommender systems bi2018, computer vision liu2012, community detection anandkumar2014community, longitudinal data analysis hoff2015multilinear, among others. Most of the developed tensor-based methods were designed for independent and identically distributed (i.i.d.) tensor data or tensor data with i.i.d. noise.
On the other hand, in many applications, the tensors are observed over time, and hence form a tensor-valued time series. For example, the monthly import export volumes of multi-categories of products (e.g. Chemical, Food, Machinery and Electronic, and Footwear and Headwear) among countries naturally form a dynamic sequence of 3-way tensor-variates, each of which representing a weighted directional transportation network. Another example is functional MRI, which typically consists hundreds of thousands of voxels observed over time. A sequence of 2-D or 3-D images can also be modeled as matrix or tensor time series to preserve temporal structure. Development of statistical methods for analyzing such large scale tensor-valued time series is still in its infancy.
In many settings, although the observed tensors are of high order and high dimension, there is often hidden low-rank structures in the tensors that can be exploited to facilitate the data analysis. Such a low-rank condition provides convenient decomposable structures and has been widely used in tensor data analysis. Two common choices of low-rank tensor structures are CANDECOMP/PARAFAC (CP) structure and multilinear/Tucker structure, and each of them has their respective benefits; see the survey in kolda2009tensor.
In dynamic data, the low-rank structures are often realized through factor models, one of the most effective and popular dimension reduction tools. Over the past few decades, there has been a large body of literature in the statistics and econometrics communities on factor models for vector time series. An incomplete list of the publications includes chamberlain1983, bai2002, stock2002, bai2003, fan2011, fan2013, forni2000, forni2004generalized, forni2005generalized, fan2016, pena1987identifying, pan2008, lam2011, lam2012. Recently, the factor model approach has been developed for analyzing high dimensional dynamic tensor time series wang2019, chen2019constrained, chen2023statistical, chen2024semi, chen2022factor, han2020iterative, han2022rank, chang2023modelling. These existing works utilize the Tucker low-rank structure in formulating the factor models. Such Tucker type tensor factor model is also closely related to separable factor analysis in fosdick2014separable under the array Normal distribution of hoff2011separable.
In this paper, we investigate a tensor factor model with a CP type low-rank structure, called TFM-cp. Specifically, let ${\cal X}_t$ be an order $K$ tensor of dimensions $d_1\times d_2\times\ldots\times d_K$. We assume
where $\otimes$ denotes tensor product, $w_i>0$ represents the signal strength, ${\mbox{\boldmath $ a$}}_{ik}$, $i=1,\ldots, r$, are unit vectors of dimension ${d_k}$, with $\|{\mbox{\boldmath $ a$}}_{ik}\|_2=1$, ${\cal E}_t$ is a noise tensor of the same dimension as ${\cal X}_t$, and \{$f_{it}$, $i=1,\ldots,r$\} is a set of uncorrelated univariate latent factor processes. That is, the signal part of the observed tensor at time $t$ is a linear combination of $r$ rank-one tensors, $w_i{\mbox{\boldmath $ a$}}_{i1}\otimes{\mbox{\boldmath $ a$}}_{i2}\otimes\cdots\otimes{\mbox{\boldmath $ a$}}_{iK}$. These rank-one tensors are fixed and do not change over time. Here, $\{{\mbox{\boldmath $ a$}}_{ik}, 1\le i\le r, 1\le k \le K\}$ are called loading vectors and the loading vectors for each mode, $\{{\mbox{\boldmath $ a$}}_{ik}, 1\le i\le r\}$, are not necessarily orthogonal. The dynamics of the tensor time series are driven by the $r$ univariate latent processes $f_{it}$. By stacking the fibers of the tensor ${\cal X}_t$ into a vector, the TFM-cp can be written as a vector factor model, with $r$ factors and a $d\times r$ (where $d=d_1\ldots d_K$) loading matrix of a special structure induced by the TFM-cp. More detailed discussion of the model is given in Section (ref).
A standard approach for dynamic factor model estimation is through the analysis of the covariance or autocovariance of the observed process. The autocovariance of a TFM-cp process in (ref) is also a tensor with a low-rank CP structure. Hence, potentially the estimation of (ref) can be done with a tensor CP decomposition procedure. However, tensor CP decomposition is well known to be a notoriously challenging problem as it is in general NP hard to compute and the CP rank is not lower semi-continuous haastad1990tensor, kolda2009tensor, hillar2013most. There are a number of works on tensor CP decomposition, which is often called tensor principal component analysis (PCA) in the literature, including alternating least squares comon2009tensor, robust tensor power methods with orthogonal components anandkumar2014tensor, tensor unfolding approaches montanari2014statistical, wang2017tensor, rank-one alternating least squares anandkumar2014guaranteed, sun2017provable, and simultaneous matrix diagonalization kuleshov2015tensor. See also zhou2013tensor,wangmy2017tensor,hao2020sparse, wang2020learning, auddy2023perturbation, han2023guaranteed, among others. Although these methods can be used directly to obtain the low-rank CP components of the autocovariance tensors, they have been designed for general tensors and do not utilize the special structure embedded in the TFM-cp.
In this paper, we develop a new estimation procedure, named as {\bf H}igh-{\bf O}rder {\bf P}rojection {\bf E}stimators (HOPE), for TFM-cp in (ref). The procedure includes a warm-start initialization using a newly developed composite principal component analysis (cPCA), and an iterative simultaneous orthogonalization scheme to refine the estimator. The procedure is designed to take the advantage of the special structure of TFM-cp whose autocovariance tensor has a specific CP structure with components close to being orthogonal and of a high-order coherence in a multiplicative form. The proposed cPCA takes advantage of this feature so the initialization is better than using random projection initialization often used in generic CP decomposition algorithms. The refinement step makes use of the multiplicative coherence again and is better than the alternating least squares, the iterative projection algorithm han2020iterative, and other forms of the high order orthogonal iteration (HOOI) de2000, liu2014, zhang2018tensor. Our theoretical analysis provides details of these improvements.
In the theoretical analysis, we establish statistical upper bounds on the estimation errors of the factor loading vectors for the proposed algorithms. The cPCA yields useful and good initial estimators with less restrictive conditions, and the iterative algorithm provides faster statistical error rates under weaker conditions than the generic CP decomposition algorithms. For cPCA, the number of factors $r$ can increase with the dimensions of the tensor time series and is allowed to be larger than $\max_k d_k$. We also derive the statistical guarantees of the iterative algorithm under the settings where the tensor is (sufficiently) undercomplete ($r \ll \min_k d_k$). It is worth noting that the iterative refinement algorithm has much sharper upper bounds for the statistical error than the cPCA initial estimators.
The TFM-cp in (ref) can also be written as a tensor factor model with a Tucker form (TFM-tucker) of a special structure. See (ref) for the definition of TFM-tucker and Remark (ref) below for the comparison between TFM-cp and TFM-tucker from the perspectives of modeling assumptions and interpretations. Potentially, the iterative estimation procedures designed for TFM-tucker can also be used here han2020iterative, ignoring the special TFM-cp structure. However, HOPE {has lower computational complexity per iteration}, requires less restrictive conditions and exhibits faster convergence rate, by fully utilizing the structure of TFM-cp. See Remark (ref) for further discussion. They also share the nice properties that the increase in either the dimensions $d_1,\ldots, d_k$, or the sample size can improve the estimation of the factor loading vectors or spaces.
The rest of the paper is organized as follows. After a brief introduction of the basic notations and preliminaries of tensor analysis in Section (ref), we introduce a tensor factor model with CP low-rank structure in Section (ref). The estimation procedures of the factors and the loading vectors are presented in Section (ref). Section (ref) investigates the theoretical properties of the proposed methods. Section (ref) develops some alternative algorithms to tensor factor models, which extend existing popular CP methods to the auto-covariance tensors with cPCA as initialization, and provides some simulation studies to demonstrate the numerical performance of all the estimation procedures. Section (ref) illustrates the model and its interpretations in real data applications. Section (ref) provides a short concluding remark. All technical details and more simulation results are relegated to the supplementary materials.
The following basic notations and preliminaries will be used throughout the paper. Define $\|x\|_q = (x_1^q+...+x_p^q)^{1/q}$, $q\ge 1$, for any vector $x=(x_1,...,x_p)^\top$. The matrix spectral norm is denoted as $$\|A\|_{\rm S}= \max_{\|x\|_2=1,\|y\|_2= 1} \|x^\top A y\|_2.$$ 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$) holds 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 exists a constant $C$ such that $a_n\le Cb_n$ (resp. $a_n\ge Cb_n$).
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} .$$ The $k$-mode product of ${\cal A}\in\mathbb{R}^{r_1\times r_2\times \cdots \times r_K}$ with a matrix $U\in\mathbb{R}^{m_k\times r_k}$ is an order $K$-tensor of size $r_1\times \cdots \times r_{k-1} \times m_k\times r_{k+1} \times \cdots \times r_K$ and will be denoted as ${\cal A}\times_k U$, such that $$ ({\cal A}\times_k U)_{i_1,...,i_{k-1},j,i_{k+1},...,i_K}=\sum_{i_k=1}^{r_k} {\cal A}_{i_1,i_2,...,i_K} U_{j,i_k}. $$
Given ${\cal A}\in \mathbb{R}^{m_1\times\cdots\times m_K}$ and $m=\prod_{j=1}^K m_j$, let ${\rm{vec}}({\cal A})\in \mathbb{R}^m$ be vectorization of the matrix/tensor ${\cal A}$, $\hbox{\rm mat}_k({\cal A})\in \mathbb{R}^{m_k\times(m/m_k)}$ the mode-$k$ matrix unfolding of ${\cal A}$, and $\hbox{\rm mat}_k(\hbox{\rm vec}({\cal A}))=\hbox{\rm mat}_k({\cal A})$.
Again, we specifically consider the following tensor factor model with CP low-rank structure (TFM-cp) for observations ${\cal X}_t\in\mathbb{R}^{d_1\times \cdots \times d_K}$, $1\le t\le T$,
where $f_{it}$ is the unobserved latent factor process and ${\mbox{\boldmath $ a$}}_{ik}$ are the fixed unknown factor loading vectors. We assume without loss of generality, $\mathbb{E} f_{it}^2=1$, $\|{\mbox{\boldmath $ a$}}_{ik}\|_2=1$, for all $1\le i\le r$ and $1\le k\le K$. Then, all the signal strengths are contained in $w_i$. A key assumption of TFM-cp is that the factor process $f_{it}$ is assumed to be uncorrelated across different factor processes, e.g., $\mathbb{E} f_{it-h}f_{jt}=0$ for $i\neq j$ and $h\ge 1$. In addition, we assume that the noise tensor ${\cal E}_t$ are uncorrelated (white) across time, but with an arbitrary contemporary covariance structure, following lam2012, chen2022factor. In this paper, we consider the case that the order of the tensor $K$ is fixed but the dimensions $d_1,...,d_K\to\infty$ and rank $r$ can be fixed or diverging.
In this section, we focus on the estimation of the factors and loading vectors of model (ref). The proposed procedure includes two steps: an initialization step using a new composite PCA (cPCA) procedure, presented in Algorithm (ref), and an iterative refinement step using a new iterative simultaneous orthogonalization (ISO) procedure, presented in Algorithm (ref). We call this two-step procedure HOPE ({\bf H}igh-{\bf O}rder {\bf P}rojection {\bf E}stimators) as it repeatedly perform high order projections on high order moments of the tensor observations. It utilizes the special structure of the model and leads to higher statistical and computational efficiency, which will be demonstrated later.
For ${\cal X}_t$ following (ref), the lagged cross-product operator, denoted by ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h$, is the $(2K)$-tensor satisfying
for a given $h\ge 1$, where $\lambda_{i,h}=w_i^2\mathbb{E} f_{i,t-h} f_{i,t}$. Note that the tensor ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h$ is expressed in a CP-decomposition form with each ${\mbox{\boldmath $ a$}}_{ik}$ used twice. Let $\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_h$ be the sample version of ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h$,
When ${\cal X}_t$ is weakly stationary and ${\cal E}_t$ is white noise, a natural approach to estimating the loading vectors is via minimizing the empirical squared loss
where the Hilbert Schmidt norm for a tensor ${\cal A}$ is defined as $ \|{\cal A}\|_{{\rm HS}}=\| \hbox{\rm vec}({\cal A})\|_2.$ In other words, ${\mbox{\boldmath $ a$}}_{i1}\otimes{\mbox{\boldmath $ a$}}_{i2}\otimes\cdots\otimes {\mbox{\boldmath $ a$}}_{iK}$ can be estimated by the leading {\it principal component} of the sample auto-covariance tensor $\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_h$. However, due to the non-convexity of (ref) or its variants, a straightforward implementation of many local search algorithms, such as gradient descent and alternating minimization, may easily get trapped into local optima and result in sub-optimal statistical performance. As shown by auffinger2013random, there could be an exponential number of local optima and the great majority of these local optima are far from the best low rank approximation. However, if we start from an appropriate initialization not too far from the global optimum, then a local optimum reached may be as good an estimator as the global optimum. A critical task in estimating the factor loading vectors is thus to obtain good initialization.
We develop a warm initialization procedure, the composite PCA (cPCA) procedure. Note that, if we unfold ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h$ into a $d\times d$ matrix ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h^*$, where $d=d_1\ldots d_K$, then (ref) implies that
a sum of $r$ rank-one matrices, each of the form ${\mbox{\boldmath $ a$}}_i{\mbox{\boldmath $ a$}}_i^\top$, where ${\mbox{\boldmath $ a$}}_i = \hbox{\rm vec}(\otimes_{k=1}^K {\mbox{\boldmath $ a$}}_{ik})$. This is very close to the principal component decomposition of ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h^*$, except that ${\mbox{\boldmath $ a$}}_i$'s are not necessarily orthogonal in this case. However, the following intuition provides a solid justification of using PCA to obtain an estimate of ${\mbox{\boldmath $ a$}}_i$. We call this estimator the cPCA estimator.
The accuracy of using the principal components of ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h^*$ as the estimate of ${\mbox{\boldmath $ a$}}_i$ heavily depends on the coherence of the components, defined as $\vartheta = \max_{1\le i < j\le r}| {\mbox{\boldmath $ a$}}_{i}^\top {\mbox{\boldmath $ a$}}_{j}|$, the maximum pairwise correlation among the ${\mbox{\boldmath $ a$}}_i$'s. When the components are orthogonal ($\vartheta=0$), there is no error in using PCA. The main idea of cPCA is to take advantage of the special structure of TFM-cp, which leads to a multiplicative high-order coherence of the CP components. In the following, we provide an analysis of $\vartheta$ under TFM-cp.
Let ${\mbox{\boldmath $ A$}}_k = ({\mbox{\boldmath $ a$}}_{1k},\ldots,{\mbox{\boldmath $ a$}}_{rk})\in \mathbb{R}^{d_k\times r}$ be the matrix with ${\mbox{\boldmath $ a$}}_{ik}$ as its columns, and ${\mbox{\boldmath $ A$}}_k^\top {\mbox{\boldmath $ A$}}_k = (\sigma_{ij,k})_{r\times r}$. As $\sigma_{ii,k}=\|{\mbox{\boldmath $ a$}}_{ik}\|_2^2=1$, the correlation among columns of ${\mbox{\boldmath $ A$}}_k$ can be measured by
Similarly we use
to measure the correlation of the matrix ${\mbox{\boldmath $ A$}} = ({\mbox{\boldmath $ a$}}_1,\ldots,{\mbox{\boldmath $ a$}}_r)\in \mathbb{R}^{d\times r}$ with ${\mbox{\boldmath $ a$}}_i = \hbox{\rm vec}(\otimes_{k=1}^K {\mbox{\boldmath $ a$}}_{ik})$ and $d=\prod_{k=1}^K d_k$. It can be seen that the coherence $\vartheta$ has the bound $\vartheta\le \prod_{k=1}^K\vartheta_k\le \vartheta_{\max}^K$, due to ${\mbox{\boldmath $ a$}}_{i}^\top {\mbox{\boldmath $ a$}}_{j} = \prod_{k=1}^K {\mbox{\boldmath $ a$}}_{ik}^\top {\mbox{\boldmath $ a$}}_{jk} =\prod_{k=1}^K \sigma_{ij,k}$. The spectrum norm $\delta$ is also bounded by the multiplicative of correlation measures in (ref). More specifically, we have the following proposition.
When (most of) the quantities in (ref) are small, the products in (ref) would be very small so that the ${\mbox{\boldmath $ a$}}_i$'s are nearly orthogonal. For example, if $({\mbox{\boldmath $ a$}}_{11},{\mbox{\boldmath $ a$}}_{21})$ and $({\mbox{\boldmath $ a$}}_{12},{\mbox{\boldmath $ a$}}_{22})$ both have i.i.d. bi-variate random rows with correlation coefficients $\rho_1$ and $\rho_2$, and independent, then the population correlation coefficient of $\hbox{\rm vec}({\mbox{\boldmath $ a$}}_{11}\otimes{\mbox{\boldmath $ a$}}_{12})$ and $\hbox{\rm vec}({\mbox{\boldmath $ a$}}_{21}\otimes{\mbox{\boldmath $ a$}}_{22})$ is $\rho_1\rho_2$, though the variation of the sample correlation coefficient depends on the length of the ${\mbox{\boldmath $ a$}}_{ik}$'s.
The pseudo-code of cPCA is provided in Algorithm (ref). Though ${\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^*$ is symmetric, its sample version $\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^*$ in general is not. We use $(\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^*+\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^{*\top})/2$ to ensure symmetry and reduce the noise. The cPCA produces definitive initialization vectors up to the sign change.
After obtaining a warm start via cPCA (Algorithm (ref)), we engage an iterative simultaneous orthogonalization (ISO) algorithm (Algorithm (ref)) to refine the solution of ${\mbox{\boldmath $ a$}}_{ik}$ and obtain estimations of the factor process $f_{it}$ and the signal strength $w_i$. Algorithm (ref) can be viewed as an extension of HOOI de2000,zhang2018tensor and the iterative projection algorithm in han2020iterative to undercomplete ($r < d_{\min}$) and non-orthogonal CP decompositions. It is motivated by the following observation. Define ${\mbox{\boldmath $ A$}}_k=({\mbox{\boldmath $ a$}}_{1k},\ldots,{\mbox{\boldmath $ a$}}_{rk})$ and ${\mbox{\boldmath $ B$}}_k = {\mbox{\boldmath $ A$}}_k({\mbox{\boldmath $ A$}}_k^{\top} {\mbox{\boldmath $ A$}}_k)^{-1} = ({\mbox{\boldmath $ b$}}_{1k},...,{\mbox{\boldmath $ b$}}_{rk}) \in\mathbb{R}^{d_k\times r}$. Let
Since ${\mbox{\boldmath $ a$}}_{jk}^\top{\mbox{\boldmath $ b$}}_{ik}=I_{\{i=j\}}$, model (ref) implies that
Here ${\cal Z}_{t,ik}$ is a vector, and (ref) is in a factor model form with a univariate factor. The estimation of ${\mbox{\boldmath $ a$}}_{ik}$ can be done easily and much more accurately than dealing with the much larger ${\cal X}_t$. The operation in (ref) achieves two objectives. First, by multiplying a vector on every mode except the $k$-th mode to ${\cal X}_t$, it reduces the tensor to a vector. It also serves as an averaging operation to reduce the noise variation. Second, as ${\mbox{\boldmath $ b$}}_{ik}$ is orthogonal to all ${\mbox{\boldmath $ a$}}_{jk}$ except ${\mbox{\boldmath $ a$}}_{ik}$, it is an orthogonal projection operation that eliminates all $\otimes_{k=1}^K{\mbox{\boldmath $ a$}}_{jk}$ terms in (ref) except the $i$-th term, resulting in (ref). If the matrix ${\mbox{\boldmath $ A$}}_k^{\top} {\mbox{\boldmath $ A$}}_k$ is not ill-conditioned, i.e. $\{{\mbox{\boldmath $ a$}}_{ik}, 1\le i\le r\}$ are not highly correlated, then ${\mbox{\boldmath $ B$}}_k$ and all individual ${\mbox{\boldmath $ b$}}_{jk}$ are well defined and this procedure shall work well. Under proper conditions on the combined noise tensor ${\cal E}_{t,ik}^*$, estimation of the loading vectors ${\mbox{\boldmath $ a$}}_{ik}$ based on ${\cal Z}_{t,ik}$ can be made significantly more accurate, as the statistical error rate now depends on $d_k$ rather than $d_1d_2\ldots d_k$. Intuitively, ${\mbox{\boldmath $ b$}}_{ik}$ can also be viewed as a form of (normalized) residuals of ${\mbox{\boldmath $ a$}}_{ik}$ projected onto the space spanned by $\{{\mbox{\boldmath $ a$}}_{jk}, j\neq i, 1\le j\le r \}$.
In practice we do not know ${\mbox{\boldmath $ b$}}_{il}$, for $1\le i\le r$, $1\le l\le K$ and $l\neq k$. Similar to back-fitting algorithms, we iteratively estimate the loading vector ${\mbox{\boldmath $ a$}}_{ik}$ at iteration number $m$ based on
using the estimate $\widehat {\mbox{\boldmath $ b$}}_{il}^{(m-1)},~ k<l\le K$, obtained in the previous iteration and the estimate $\widehat {\mbox{\boldmath $ b$}}_{il}^{(m)},~ 1\le l< k$, obtained in the current iteration. As we shall show in the next section, such an iterative procedure leads to a much improved statistical rate in the high dimensional tensor factor model scenarios, as if all ${\mbox{\boldmath $ b$}}_{il}$, $1\le i\le r$, $1\le l\le K$, $l\neq k$, are known and we indeed observe ${\cal Z}_{t,ik}$ that follows model (ref). Note that the projection error is \[ {\cal Z}_{t,ik}^{(m)}-{\cal Z}_{t,ik} =\sum_{j=1}^r w_{j}f_{j,t}\xi_{ij}^{(m)}{\mbox{\boldmath $ a$}}_{jk} +{\cal E}_{t,ik}^{*(m)} - {\cal E}_{t,ik}^* \] where
and ${\cal E}_{t,ik}^{*(m)}$ is that in (ref) with ${\mbox{\boldmath $ b$}}_{ik}$ replaced with $\widehat{\mbox{\boldmath $ b$}}_{ik}^{(m-1)}$ or $\widehat{\mbox{\boldmath $ b$}}_{ik}^{(m)}$. The multiplicative measure of projection error $|\xi_{ij}^{(m)}|$ decays rapidly since, for $j\neq i$, ${\mbox{\boldmath $ a$}}_{j\ell}^\top \widehat{\mbox{\boldmath $ b$}}_{i\ell}^{(m)}$ goes to zero quickly as the iteration $m$ increases, and $\xi_{ij}^{(m)}$ is a product of $K-1$ such terms. In fact, the higher the tensor order $K$ is, the faster the error goes to zero.
In this section, we shall investigate the statistical properties of the proposed algorithms described in the last section. Our theories provide theoretical guarantees for consistency and present statistical error rates in the estimation of the factor loading vectors ${\mbox{\boldmath $ a$}}_{ik}$, $1\le i\le r, 1\le k\le K$, under proper regularity conditions. As the loading vector ${\mbox{\boldmath $ a$}}_{ik}$ is identifiable only up to the sign change, we use
to measure the distance between $\widehat{\mbox{\boldmath $ a$}}_{ik}$ and ${\mbox{\boldmath $ a$}}_{ik}$.
Recall ${\mbox{\boldmath $ \mathnormal\Sigma$}}_{h} = \mathbb{E}\widehat{\mbox{\boldmath $ \mathnormal\Sigma$}}_{h} =\sum_{i=1}^r \lambda_{i,h} ({\mbox{\boldmath $ a$}}_{i1}\otimes{\mbox{\boldmath $ a$}}_{i2}\otimes\cdots\otimes {\mbox{\boldmath $ a$}}_{iK})^{\otimes2}, $ as in (ref) and $\lambda_{i,h}= w_i^2\mathbb{E} f_{i,t-h}f_{i,t}$. We will also continue to use the notations ${\mbox{\boldmath $ A$}}_k=({\mbox{\boldmath $ a$}}_{1k},...,{\mbox{\boldmath $ a$}}_{rk})\in\mathbb{R}^{d_k\times r}$, and ${\mbox{\boldmath $ B$}}_k = {\mbox{\boldmath $ A$}}_k({\mbox{\boldmath $ A$}}_k^{\top} {\mbox{\boldmath $ A$}}_k)^{-1} = ({\mbox{\boldmath $ b$}}_{1k},...,{\mbox{\boldmath $ b$}}_{rk}) \in\mathbb{R}^{d_k\times r}$. Let $d=\prod_{k=1}^K d_k$, $d_{\min}=\min\{d_1,...,d_K\}$, $d_{\max}=\max\{d_1,...,d_K\}$ and $d_{-k}=\prod_{j\ne k}d_j$.
To present theoretical properties of the proposed procedures, we impose the following assumptions.
Assumption (ref) is similar to those on the noise imposed in lam2011, lam2012, han2020iterative. It accommodates general patterns of dependence among individual time series fibers, but also allows a presentation of the main results with manageable analytical complexity. The normality assumption, which ensures fast statistical error rates in our analysis, is imposed for technical convenience. In theory, it can be supplanted by a sub-Gaussian condition, or replaced with more general distributions with heavier tails. However, adopting such weaker conditions would significantly complicate the formulae, statistical outcomes, and requisite conditions within our time series framework. This added complexity would likely detract from the paper's readability, without providing additional statistical insights. Our primary goal is to maintain the readers' focus on the core content of the paper, and thus, we have chosen to assume additive Gaussian errors. This choice simplifies the exposition without compromising the fundamental tenets and the findings of our study.
Assumption (ref) is standard. It allows a very general class of time series models, including causal ARMA processes with continuously distributed innovations; see also tong1990non, bradley2005, tsay2005analysis, fan2008nonlinear, rosenblatt2012markov, tsay2018nonlinear, among others. The restriction $\gamma_1\le 1$ is introduced only for presentation convenience. Assumption (ref) requires that the tail probability of $f_{it}$ decays exponentially fast. In particular, when $\gamma_2=2$, $f_{it}$ is sub-Gaussian.
Assumption (ref) is sufficient to guarantee that all the factor loading vectors ${\mbox{\boldmath $ a$}}_{ik}$ can be uniquely identified up to the sign change. The parameters $\lambda_i$ can be viewed as an analogue of eigenvalues in the order-$2K$ tensor ${\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}$. Similar to the eigen decomposition of a matrix, if some $\lambda_i$ are equal, estimation of the loading vectors ${\mbox{\boldmath $ a$}}_{ik}$ may suffer from label shift across $i$. When $h$ is fixed and $\mathbb{E} f_{i,t} f_{i,t-h}\asymp 1$, $\lambda_i\asymp w_i^2$. The signal strength of each factor is measured by $\lambda_i$.
Let us first study the behavior of the cPCA estimators in Algorithm (ref). Theorem (ref) presents the performance bounds, which depends on the coherence (the degree of non-orthogonality) of the factor loading vectors.
Let
with $\lambda_0=\infty$, $\lambda_{r+1}=0$, be the minimum gap between the signal strengths of the factors.
The first term in the upper bound (ref) is induced by the non-orthogonality of the loading vectors ${\mbox{\boldmath $ a$}}_{ik}$, which can be viewed as bias. The second term in (ref) comes from a concentration bound for the random noise, and thus can be interpreted as stochastic error. By Proposition (ref), it implies that a larger $K$ (e.g. higher order tensors) leads to smaller bias and higher statistical accuracy of cPCA. If $\delta\gtrsim R^{(0)}/\lambda_1$, then the error bound (ref) is dominated by the bias related to $\delta$, otherwise it is dominated by the stochastic error. Equation (ref) shows that $R^{(0)}$ in the stochastic error comes from the fluctuation of the factor process $f_{it}$ (the first two terms) and the noise ${\cal E}_t$ in (ref) (the other two terms). When $\sum_{i=1}^r\lambda_i\asymp r\lambda_1 \asymp r w_1^2$ and $R^{(0)}/\lambda_1 +T/d \lesssim 1$, the terms related to the noise becomes $\sqrt{r/T}/\sqrt{{\rm SNR}}$, where the signal-to-noise ratio (SNR) is
Roughly speaking though not completely correct, the term $\lambda_i-\lambda_{i+1}$ can be viewed as the gap of $i$-th and $(i+1)$-th largest eigenvalues of ${\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^*$ with ${\mbox{\boldmath $ \mathnormal\Sigma$}}_{h}^*$ given in (ref). In particular, if $\lambda_1\asymp ...\asymp \lambda_r \asymp w_1^2$, then $\lambda_*\asymp w_1^2/r$. In this case, the bound (ref) can be simplified to
Then, by (ref) and Proposition (ref), the consistency of the cPCA estimators only requires the incoherence parameter to be at most $\vartheta_{\max}\lesssim r^{-2/K}$.
Next, let us consider the statistical performance of the iterative algorithm (Algorithm (ref)) after cPCA initialization, i.e. HOPE estimators. As discussed earlier, the operation in (ref) achieves dimension reduction by projecting ${\cal X}_t$ into a vector and retains only one of the $r$ factor terms, hence eliminates the interaction effects between different factors. As we update the estimation of each individual loading vector ${\mbox{\boldmath $ a$}}_{ik}$ separately in the algorithm, ideally this would remove the bias part in (ref) which is due to the non-orthogonality of the loading vectors, and replace the eigengap $\lambda_*$ in (ref) by $\lambda_i$, as (ref) only involves one eigenvector. It also leads to the elimination of the first two terms of $R^{(0)}$. As mentioned in Section (ref), when updating $\widehat {\mbox{\boldmath $ a$}}_{ik}^{(m)}$, we take advantages of the multiplicative nature of the project error $\xi_{ij}^{(m)}$ in (ref), and the rapid growth of such benefits as the iteration number $m$ grows. Thus we expect that the rate of HOPE estimators would become
where
Note that $R_{k,i}^{(\text{\footnotesize ideal})}$ replaces all $d=d_1\ldots d_K$ in the noise component of the stochastic error in (ref) by $d_k$ due to dimension reduction. The following theorem provides conditions under which this ideal rate is indeed achieved.
Let the statistical error bound of the initialization used in Algorithm (ref) be $\psi_0$. For cPCA,
where $\lambda_*$ is the eigengap defined in (ref) and $R^{(0)}$ is defined in (ref).
The detailed proof of the theorem is in Appendix (ref). The key idea of the analysis of HOPE is to show that the iterative estimator has an error contraction effect in each iteration. Theorem (ref) implies that HOPE will achieve a faster statistical error rate than the typical $O_{\mathbb{P}}(T^{-1/2})$ whenever $\lambda_r\gg \sigma^2 \max_k d_k$. As $\hbox{\rm vec}({\cal X}_t)$ has $d$ elements, the strong factors setting in the literature lam2011,chen2022factor,han2020iterative typically assumes ${\rm SNR}\asymp 1$. In our case it is similar to assuming the signal strength $\mathbb{E}\big\|\sum_{i=1}^r w_i f_{it}\otimes_{k=1}^K {\mbox{\boldmath $ a$}}_{ik}\big\|_{\rm HS}^2 \asymp \sigma^2 d$. When $r$ is fixed and $\lambda_1 \asymp \cdots \asymp \lambda_r$, the statistical error rate will be reduced to $O_{\mathbb{P}}(T^{-1/2}d_{-k}^{-1/2})$, where $d_{-k}=\prod_{j\ne k}d_j$.
Theorem (ref) specifies the convergence rate for the estimated factors $f_{it}$. When $\lambda_r \gg \sigma^2 d_{\max} +T$, $w_i^{-1} | \widehat w_i^{{\rm\tiny iso}} \widehat f_{it}^{{\rm\tiny iso}} - w_i f_{it}|$ is much smaller than the parametric rate $T^{-1/2}$. If all the factors are strong lam2011 such that $\lambda_1\asymp\lambda_r \asymp \sigma^2 d$, (ref) implies that $w_i^{-1} | \widehat w_i^{{\rm\tiny iso}} \widehat f_{it}^{{\rm\tiny iso}} - w_i f_{it} |=O_{\mathbb{P}}(d^{-1/2}+d_{\max}^{1/2}d^{-1/2}T^{-1/2})$. Then, as long as $d_k\to \infty$ and $K\ge 2$, the estimated factors are consistent, even under a fixed $T$. In comparison, the convergence rate of the estimated factors in Theorem 1 of bai2003 for vector factor models is $O_{\mathbb{P}}(d^{-1/2}+T^{-1})$. Moreover, (ref) shows that the error rates for the sample auto-cross-moment of the estimated factors to the true sample auto-cross-moment is also $o_{\mathbb{P}}(T^{-1/2})$ when $\lambda_r \gg \sigma^2 d_{\max}$. This implies that it is a valid option to use the estimated factor processes as the true factor processes to model the dynamics of the factors. When the estimation of these time series models only {required} auto-correlation and partial auto-correction functions, the results are expected to be the same as using the true factor process, without loss of efficiency. The statistical rates in Theorem (ref) lay a foundation for further modeling of the estimated factor processes with vast repository of linear and nonlinear options.
Here we present two alternative estimation algorithms for TFM-cp, by extending the popular rank one alternating least square (ALS) algorithm of anandkumar2014guaranteed and orthogonalized alternating least square (OALS) of sharan2017orthogonalized designed for CP decomposition of noisy tensors, because $\widehat {\mbox{\boldmath $ \mathnormal\Sigma$}}_h$ in (ref) is indeed in a CP form, but with repeated components. In addition, we use cPCA estimates for initialization, instead of randomized initialization used for general CP decomposition. We will denote the algorithms as cALS (Algorithm (ref)) and cOALS (Algorithm (ref)), respectively. The simulation study below shows that, although cALS and cOALS perform better than the straightforward implementation of ALS and OALS with randomized initialization, they do not perform as well as the proposed HOPE algorithm. Hence we do not investigate their theoretical properties in this paper.
In this section, we compare the empirical performance of different procedures of estimating the loading vectors of TFM-cp, under various simulation setups. We consider the cPCA initialization (Algorithm (ref)) alone, the iterative procedure HOPE, and the intermediate output from the iterative procedure when the number of iteration is 1 after initialization. The one step procedure will be denoted as 1HOPE. We also check the performance of the alternative algorithms ALS, OALS, cALS, and cOALS as described above. The estimation error shown is given by $\max_{i,k}\|\widehat {\mbox{\boldmath $ a$}}_{ik}\widehat {\mbox{\boldmath $ a$}}_{ik}^\top - {\mbox{\boldmath $ a$}}_{ik} {\mbox{\boldmath $ a$}}_{ik}^\top \|_{\rm S}$.
We demonstrate the performance of all procedures under TFM-cp with $K=2$ (matrix time series) with
For $K=2$ with model (ref), we consider the following three experimental configurations:
Results from an additional simulation settings under $K=2$ and $K=3$ cases are given in Appendix (ref). We repeat all the experiments 100 times. For simplicity, we set $h=1$.
The loading vectors are generated as follows. First, the elements of matrices $\widetilde {\mbox{\boldmath $ A$}}_{k}=(\widetilde {\mbox{\boldmath $ a$}}_{1k},..., \widetilde {\mbox{\boldmath $ a$}}_{rk})\in \mathbb{R}^{d_k\times r}$, $1\le k\le K$, are generated from i.i.d. $N(0,1)$ and then orthonormalized through QR decomposition. Then if $\delta=0$, set ${\mbox{\boldmath $ A$}}_{k}=\widetilde {\mbox{\boldmath $ A$}}_{k}$, otherwise, set ${\mbox{\boldmath $ a$}}_{1k}=\widetilde {\mbox{\boldmath $ a$}}_{1k}$ and ${\mbox{\boldmath $ a$}}_{ik}= (\widetilde {\mbox{\boldmath $ a$}}_{1k}+\theta \widetilde {\mbox{\boldmath $ a$}}_{ik})/\|\widetilde {\mbox{\boldmath $ a$}}_{1k}+\theta \widetilde {\mbox{\boldmath $ a$}}_{ik}\|_2$ for all $i\ge 2$ and $1\le k\le K$, with $\vartheta=\delta/(r-1)$ and $\theta=(\vartheta^{-2/K}-1)^{1/2}$. The commonly used incoherence measure anandkumar2014guaranteed,hao2020sparse under this construction is $\vartheta_{\max}=(1+\theta^2)^{-1/2}=\vartheta^{1/K}$.
The noise ${\cal E}_t$ in the model is white ${\cal E}_t\perp {\cal E}_{t+h},h>0$, and generated according to ${\cal E}_t=\Psi_1^{1/2} Z_t\Psi_2^{1/2}$ where all of the elements in the $d_1\times d_2$ matrix $Z_t$ are i.i.d. $N(0,1)$. Furthermore, $\Psi_1, ~\Psi_2$ are the covariance matrices along each mode with the diagonal elements being $1$ and all the off-diagonal elements being $\psi_1,~\psi_2$. Throughout this section, we set the off-diagonal entries of the covariance matrices of the noise as $\psi_1=\psi_2$.
Under Configurations I and II with $r=2$, the factor processes $f_{1t}$ and $f_{2t}$ are generated as two independent AR(1) processes, following $f_{1t}=0.8f_{1t-1}+e_{1t}$, $f_{2t}=0.6f_{2t-1}+e_{2t}$. Under Configuration III and Configurations IV and V in Appendix (ref), with $r=3$, $f_{1t},f_{2t},f_{3t}$ are generated as independent AR(1) processes, with $f_{1t}=0.8f_{1t-1}+e_{1t}, f_{2t}=0.7f_{2t-1}+e_{2t}, f_{3t}=0.6f_{3t-1}+e_{3t}$. Here, all of the innovations follow i.i.d. $N(0,1)$. The factors are not normalized.
Figure (ref) shows the boxplots of the estimation errors for cPCA and HOPE under configuration I, for different $\delta$. It can be seen that the performance of cPCA deteriorates as $\delta$ increases, while that of HOPE remains almost unchanged. The median of the cPCA estimation errors increases almost linearly with $\delta$, with a $R^2$ of $0.977$. This linear effect of $\delta$ on the performance bounds of cPCA is confirmed by the theoretical results in (ref).
The experiment of Configuration II is conducted to verify the theoretical bounds on different sample sizes $T$ and signal strengths $w$. Figures (ref) and (ref) show the logarithm of the estimation errors under different ($w,T$) combinations. It can be seen from Figure (ref) that the estimation error of cPCA decreases to a lower bound as $w$ and $T$ increases. The lower bound is associated with the bias term in (ref) that cannot be reduced by a larger $w$ and $T$. This is the baseline error due to the non-orthogonality. In contrast, the phenomenon of HOPE is very different. Figure (ref) shows that the performance improves monotonically as $w$ or $T$ increases. Again, this is consistent with the theoretical bounds in (ref).
Figure (ref) shows the boxplots of the logarithm of the estimation errors for 7 different methods with choices of $\delta$ under configuration III. ALS and OALS are implemented with $L=200$ random initiations. It can be seen that HOPE outperforms all the other methods. Again, the choice of $\delta$ does not affect the performance of HOPE significantly. One-step method (1HOPE) is better than the cPCA alone, and the iterative method HOPE is in turn better than the one-step method. When the coherence $\delta$ decreases, all methods perform better, but the advantage of HOPE over one-step method and the advantage of one-step method over the cPCA initialization become smaller. For the extremely small $\delta=0.01$, all loading vectors are almost orthogonal to each other. In this case, all the iterative procedures, including the one-step HOPE, perform similarly. In addition, ALS and cALS are always the worst under the cases $\delta \ge 0.1$. The hybrid methods cALS and cOALS improve the original randomized initialized ALS and OALS significantly, showing the advantages of the cPCA initialization. It is worth noting that cOALS has comparable performance with 1HOPE and HOPE when $\delta$ is small.
In this section, we demonstrate the use of TFM-cp model using the taxi traffic data set used in chen2022factor. The data set was collected by the Taxi & Limousine Commission of New York City, and published at {\rm https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page.} Within Manhattan Island, it contains 69 predefined pick-up and drop-off zones and 24 hourly period for each day from January 1, 2009 to December 31, 2017. The total number of rides moving among the zones within each hour is recorded, yielding a ${\cal X}_t \in \mathbb{R}^{69\times 69\times 24}$ tensor for each day, using the hour of day as the third dimension.
One natural way to model the hourly traffic data is to consider ${\mbox{\boldmath $ X$}}_{t}\in\mathbb{R}^{69\times 69}$ for $t=1,\dots,24\times N$, where $N$ is the total number of days. It is apparent that such hourly data have two types of seasonality: weekly seasonality and daily seasonality. We handle the weekly seasonality by separating the time series into two parts: business-day series and non-business-day series. The length of the business-day series is 2,262 days, and that of the non-business-day is 1,025 days. The daily seasonality is a more interesting and important issue. It is clear that taxi usage heavily depends on time of the day (morning and evening rush hours, lunch hours etc). Seasonal time series models have been extensively studied for univariate time series Box&Jenkins76, tsay2005analysis,Shumway&Stoffer06,reinsel2003elements, but these parametric models are difficult to be extended to deal with high dimensional matrix time series. Segmenting the 24-hour day into distinct intervals, such as morning rush hours and business hours, loses the detailed hourly information and requires pre-determined segmentation scheme zhang2019spectral,zhu2022learning. Here we adopt the nonparametric approach by stacking the hourly observations into a (multi-dimensional) daily observation. This is equivalent to turning hourly observations into daily 24-dimensional vector observations in the univariate time series case, a commonly used approach. By jointly modeling the 24-hourly observations within the day together, the detailed daily pattern can be captured more accurately and more flexibly without a parametric model assumption. This way, the daily pattern is built into the model simultaneously, and the interaction among the geographic and temporal patterns will be revealed.
After some exploratory analysis, we decide to use the TFM-cp with $r=4$ factors for both business-day series and non-business-day series, and estimate the model with $h=1$. For the non-business-day series, TFM-cp explains 63.0% of the variability in the data. In comparison, treating the tensor time series as a $114,264$ {($=69\times 69\times 24$)} dimensional vector time series, the traditional vector factor model with 4 factors explains about 90.0% variability, but uses $4\times 114,264$ parameters for the loading matrix. chen2022factor used TFM-tucker with $4\times 4\times 4$ core factor tensor process. Using iTIPUP estimator of han2020iterative, TFM-tucker explains 80.1% variability. Similarly, for the business-day series, the explained fractions of variability by the TFM-cp {with 4 factors}, vector factor model with 4 factors, and TFM-tucker with $4\times 4\times 4$ core factor tensor are 68.1%, 90.9%, 84.0%, respectively. Comparison of model complexity between TFM-cp and TFM-tucker is substantially more complex compared to that of traditional tensor decomposition, owing to the stochastic nature of the latent factor process. Since both models are estimated through the sample autocovariance tensor, we may count, under each model, the number of parameters required in the population version of the autocovariance tensor. Both models {require} $(69+69+24)\times 4$ number of parameters in terms of loading matrices or vectors, though the {degree} of freedom for TFM-tucker is slightly smaller due to orthonormal requirements of the loading matrices. However, for the lag-1 autocovariance of the factor processes, TFM-cp only requires 4 parameters (for the four uncorrelated factor process) while TFM-tucker requires $(4\times4\times 4)^2$ parameters (minus certain savings from rotation ambiguity). More importantly, TFM-cp harbors a much smaller dimensional factor process (${\mbox{\boldmath $ f$}}_t\in \mathbb{R}^4$) than TFM-tucker (${\cal F}_t\in\mathbb{R}^{4\times 4\times 4}$), thereby simplifying subsequent modeling of the latent factor process.
The literature on Markov process models, for example zhang2019spectral,zhu2022learning, regards each trip as a transition from the pickup location to the drop-off location, rendering the data as a collection of fragmented sample paths representative of a city-wide Markov process. When one aggregates the individual trips within a time period to form a traffic volume matrix, the model, with an assumed fixed (reduced rank) Markov transition matrix, essentially induces an order-1 autoregressive model on the volume matrix time series. Therefore, the distinction between the Markov process models and the CP factor models parallels that between autoregressive models and factor models.
Figures (ref) and (ref) show the heatmap of the estimated loading vectors $({\mbox{\boldmath $ a$}}_{11},\ldots, {\mbox{\boldmath $ a$}}_{41})$ (related to pick-up locations) and $({\mbox{\boldmath $ a$}}_{12},\ldots, {\mbox{\boldmath $ a$}}_{42})$ (related to drop-off locations) of the 69 zones in Manhattan, respectively, for the business-day series. Table (ref) shows the corresponding loading vectors $({\mbox{\boldmath $ a$}}_{13},\ldots, {\mbox{\boldmath $ a$}}_{43})$ on the time of day dimension. For a more meaningful interpretation, we have re-scaled the loading vectors ${\mbox{\boldmath $ a$}}_{ik}$ and the factors $w_i f_{it}$ such that $\|{\mbox{\boldmath $ a$}}_{ik}\|_1=1$, {for} $1\le i\le 4, 1\le k\le 3$. Figure (ref) shows the estimated four factors ($w_i f_{it}$) for business day series in 1,000. (Please note the significant difference in scale). For a more detailed examination, we show the four factor series in the third year (year 2011) in Figure (ref) in Appendix (ref).
It is seen that the estimated loading vectors and the factors are predominantly positive, although there are a few small negative values which we will ignore. When the loading vectors are scaled to sum to 1 (hence percentages), the model has the following interesting interpretation. First, the expected daily total volume ($\sum_{i,j,k}X_{t,ijk}$) is the sum of the four factors $w_1f_{1t}+\ldots+w_4f_{4t}$. Hence the daily traffic volumes essentially consist of taxi rides following four different patterns, each corresponding to the rank-1 tensor ${\mbox{\boldmath $ a$}}_{i1}\otimes{\mbox{\boldmath $ a$}}_{i2}\otimes{\mbox{\boldmath $ a$}}_{i3}$, $i=1,\ldots,4$. One may also imagine that there are four types of taxi users in the city, each following one specific traffic pattern (of course an individual may take multiple trips in a day and follow different patterns for each trip).
It is interesting to study the component of the rank-1 tensor ${\mbox{\boldmath $ a$}}_{i1}\otimes{\mbox{\boldmath $ a$}}_{i2}\otimes{\mbox{\boldmath $ a$}}_{i3}$. Specifically, ${\mbox{\boldmath $ a$}}_{i3}$ shows how the total volume of traffic pattern $i$ ($w_if_{it}$) is distributed to different hours of the day. For example, Table (ref) shows that 16% of pattern 2 volume $w_2f_{2t}$ is allotted to between 7am and 8am, while only 3% of pattern 4 volume $w_4f_{4t}$ is allotted to that hour. The rank-1 matrix ${\mbox{\boldmath $ a$}}_{i1}{\mbox{\boldmath $ a$}}_{i2}^\top$ shows the spatial pattern of traffic pattern $i$, where ${\mbox{\boldmath $ a$}}_{i1}$ shows the percentage of pattern $i$ volume ($w_if_{it}$) at each hour being picked up in each location, and ${\mbox{\boldmath $ a$}}_{i2}$ shows the percentage of pattern $i$ volume $w_if_{it}$ being dropped-off in each location, and $a_{i1k} a_{i2\ell}$ is the percentage of pattern $i$ volume $w_if_{it}$ from location $k$ to location $\ell$. This spatial pattern does not change through the day, but the volume in each hour is controlled by ${\mbox{\boldmath $ a$}}_{i3}w_if_{it}$.
From Table (ref), it is seen that Factor 1 (or traffic pattern 1) roughly corresponds to the evening hours of 6pm to 12am, by the loading vector ${\mbox{\boldmath $ a$}}_{13}$, with main activities in the SoHu and lower east side as both the pick-up and drop-off locations. From the estimated factor series plot in Figure (ref), it seems that this traffic pattern (pattern 1) has the largest overall volume, but with a very strong yearly seasonal pattern and a large daily variation. Intuitively people use less taxi service when the weather is nice, hence the volume is relatively small in summer and early fall, even though there are more evening activities in the summer. The large daily variation is due to a weekly effect. Figure (ref) shows the 3-month business-day period from January 1 to March 31 in year 2011, in which the vertical line marks the end of working week (Friday or the day before holiday). It is clearly seen that, for this mainly evening-activity traffic pattern, the volume in the end of working week is almost twice as large as that in the beginning of the working week.
Again from Table (ref), it is seen that Factor 2 (or traffic pattern 2) roughly corresponds to the morning rush hours of 6am to 12am, by the loading vector ${\mbox{\boldmath $ a$}}_{23}$, with main activities in the midtown area as the pick-up locations, and Times square and 5th Avenue as the drop-off locations. About $23.1\%$ (defined as $\sum_t w_2f_{2t}/\sum_t\sum_i w_if_{it}$) of the total traffic follows this pattern. From Figure (ref), it is seen that Factor 2 is quite stable throughout the year, which is again intuitively understandable as the traffic pattern is mainly used by the steady population of people commuting to work in morning rush hours. There is a large number of (small value) outliers, most of them corresponding to the business days before or after major holidays. It can be seen more clearly from Figure (ref) in Appendix (ref).
For Factors 3 and 4, the areas that load heavily on the factors for pick-up are quite similar to that for drop-off, i.e., upper east side (with affluent neighborhoods and museums) on Factor 3, and upper west side (with affluent neighborhoods and performing arts) on Factor 4. The conventional business hours are heavily and almost exclusively loaded on these factors. From Figure (ref), it seems that both patterns have a yearly seasonal effect, small in the summer and early fall, which can be seen more clearly in Figure (ref) in Appendix (ref). Their volumes are relatively small than that of Factors 1 and 2.
We note that TFM-cp representation is unique which facilitates a more “unique” interpretation. On the other hand, TFM-tucker is subject to arbitrary rotation. Using TFM-tucker to analyze the same data set, chen2022factor used varimax rotation to obtain one specific representation of their estimated model and provided interesting interpretations. Their results are quite different from that of TFM-cp. First, since TFM-tucker representation requires orthonormal loading matrices, the discovered patterns in the loading matrices are forced to be different. For example, the daily patterns revealed in chen2022factor have quite distinct periods, while {Table} (ref) shows more intertwined (non-orthogonal) patterns. Second, TFM-tucker requires $4\times 4\times 4$ factor processes. The column loading vectors in each loading matrices work on all these factors, instead of on only one factor as in TFM-cp. The interpretation of these loading vectors are more convoluted. For example, in TFM-tucker, the volume from all four heavily loaded pick-up areas identified by ${\mbox{\boldmath $ a$}}_{i1}, i=1,\ldots, 4$ can be traveling to all four heavily loaded drop-off areas ${\mbox{\boldmath $ a$}}_{i2}, i=1,\ldots, 4$. But in TFM-cp, the rank-1 matrix ${\mbox{\boldmath $ a$}}_{i1}{\mbox{\boldmath $ a$}}_{i2}^\top$ shows the exact proportion of pattern $i$ traffic from each of the pick-up area identified by ${\mbox{\boldmath $ a$}}_{i1}$ to the drop-off area ${\mbox{\boldmath $ a$}}_{i2}$. In particular, it is seen that ${\mbox{\boldmath $ a$}}_{i1}$ is very similar to ${\mbox{\boldmath $ a$}}_{i2}$ for $i=1,3,4$. This observation suggests that our TFM-cp model may offer better intuitive understanding as it aligns with the expectation that most taxi traffic activities are likely confined within specific areas. This comparison further underscores the distinct analytical insights offered by the TFM-cp model in capturing the spatial-temporal dynamics of urban taxi traffic.
For the non-business day series, the estimated loading vectors $({\mbox{\boldmath $ a$}}_{11},\ldots, {\mbox{\boldmath $ a$}}_{41})$ (related to pick-up locations), $({\mbox{\boldmath $ a$}}_{12},\ldots, {\mbox{\boldmath $ a$}}_{42})$ (related to drop-off locations), $({\mbox{\boldmath $ a$}}_{13},\ldots, {\mbox{\boldmath $ a$}}_{43})$ (on the hour of day dimension), the estimated factors $(w_1f_{1t},\ldots,w_4f_{4t})$ are showed in Figures (ref), (ref), Table (ref), and Figure (ref), respectively. Understandably the morning rush hour pattern in the business day series (Factor 2) disappears here but the night-time pattern (Factor 1) now lasts deep into the early hours, comparing Tables (ref) and (ref). From Figure (ref), it is seen that there exist two different yearly seasonal patterns. Factors 3 and 4 are similar to that of business day series, with small volumes in the summer and fall, again confirming that the use of taxi service is relatively low when the weather is good for walking in the city. On the other hand, the volume of Factor 2 is typically small in the winter time. The volumes of night-life pattern in Factor 1 remain to be volatile. It has many small-value outliers, mostly on the day before a business day (Sundays or the end of holiday.) These can be seen more clearly in the more detailed Figure (ref), which shows the estimated factors of all the non-business days in Year 2011 (year 3), with vertical lines indicating the day before a business day (dashed lines for Sundays and solid lines for Mondays of long weekend when Tuesday is the start of business week.) This is again intuitively understandable, because people tend not to stay out too late if they need to work the next day.
The pick-up and drop-off locations that heavily load on Factors 1, 3, 4 are similar to that for Factors 1, 3, 4 in the business day series. The daytime hours load on Factors 3 and 4, and the night life hours from 12am to 4am load on Factor 1. As for the second factor, it loads heavily on midtown area for pick-up, on the lower west side near Chelsea (with many restaurants and bars) for drop-off, on the afternoon/evening hours between 1pm to 8pm as the dominating periods.
We remark that this example is just for illustration and showcasing the interpretation of the proposed tensor factor model. Again we note that for the TFM-tucker model, one needs to identify a proper representation of the loading space in order to interpret the model. In chen2022factor, varimax rotation was used to find the most sparse loading matrix representation to model interpretation. For TFM-cp, the model is unique hence interpretation can be made directly. Interpretation is impossible for the vector factor model in such a high dimensional case.
In this paper, we propose a tensor factor model with a low rank CP structure and develop its corresponding estimation procedures. The estimation procedure takes advantage of the special structure of the model, resulting in faster convergence rate and more accurate estimations comparing to the standard procedures designed for the more general TFM-tucker, and the more general tensor CP decomposition. Numerical study illustrates the finite sample properties of the proposed estimators. The results show that HOPE uniformly outperforms the other methods, when the observations follow the specified TFM-cp.
The HOPE in this paper is based on CP decomposition of the second moment tensor ${\mbox{\boldmath $ \mathnormal\Sigma$}}_h=\sum_{i=1}^r \lambda_i (\otimes_{k=1}^K{\mbox{\boldmath $ a$}}_{ik})^{\otimes 2}$, an order $2K$ tensor. The intuition that higher order tensors tend to have smaller coherence among the CP components leads to the consideration of using higher order cross-moments to have more orthogonal CP components. For example, let the $m$-th cross moment tensor with lags $0=h_1<\cdots < h_m$ be
When the factor processes $f_{it}$, $i=1,\ldots, r$ are independent across different $i$ in TFM-cp, a naive 4-th cross moment tensor to estimate ${\mbox{\boldmath $ a$}}_{ik}$ is
with $\{{\cal X}_t^*\}$ being an independent coupled process of $\{{\cal X}_t\}$ and when $h_j=(j-1)h$,
This naive 4-th cross moment tensor has more orthogonal CP bases. In light tailed case, simulation shows that it is much worse than the second moment tensor, due to the reduced signal strength $\lambda_{i,h_1h_2h_3h_4}^{(4)}$. However, for heavy tailed and skewed data, this procedure would be helpful. It would be an interesting and challenging problem to develop an efficient higher cross moment tensor to improve the statistical and computational performance. We leave this for future research.
Our primary consideration was directed towards the CP factor model in a time series setting, as the need of effectively analyzing tensor time series has arisen in many applications, and CP factor model is an efficient approach for such analysis. Without a specified (parametric) model for the latent factor processes (an topic currently under investigation), we focus on the auto-covariance and auto-cross-moment tensor for effective estimation of the proposed model. This is the main contribution of the paper. However, the proposed cPCA and ISO can be used directly in or be extended to many other problems involving CP decomposition of certain type of tensors. For example, in many problems where higher order moments can be introduced to reduce incoherence, the cPCA and ISO algorithms may offer better initialization and outperform conventional tensor power iteration methods. One specific example is the kurtosis tensor in independent component analysis in auddy2023large. Another possible extension is highlighted in Remark (ref) where we pointed out how cPCA can be modified to deal with situations when a few leading singular values of the auto-cross-moment tensor are the same.