lemmasection

# High-Dimensional Low-Rank Tensor Autoregressive Time Series Modeling

Di Wang†, Yao Zheng‡ and Guodong Li†
†University of Hong Kong and ‡University of Connecticut

###### Abstract

Modern technological advances have enabled an unprecedented amount of
structured data with complex temporal dependence, urging the need for new
methods to efficiently model and forecast high-dimensional tensor-valued time
series. This paper provides the first practical tool to accomplish this task
via autoregression (AR). By considering a low-rank Tucker decomposition for
the transition tensor, the proposed tensor autoregression can flexibly capture
the underlying low-dimensional tensor dynamics, providing both substantial
dimension reduction and meaningful dynamic factor interpretation. For this
model, we introduce both low-dimensional rank-constrained estimator and high-
dimensional regularized estimators, and derive their asymptotic and non-
asymptotic properties. In particular, by leveraging the special balanced
structure of the AR transition tensor, a novel convex regularization approach,
based on the sum of nuclear norms of square matricizations, is proposed to
efficiently encourage low-rankness of the coefficient tensor. A truncation
method is further introduced to consistently select the Tucker ranks.
Simulation experiments and real data analysis demonstrate the advantages of
the proposed approach over various competing ones.

Keywords: dimension reduction; high-dimensional time series; nuclear norm;
tensor decomposition; tensor-valued data

##  1 Introduction

Modern technological development has made available enormous quantities of
data, many of which are structured and collected over time. Tensor-valued time
series data, namely observations on a set of variables structured in a tensor
form collected over time, have become increasingly common in a wide variety of
areas. Examples include multiple macroeconomic indices time series for
multiple countries (Chen et al., 2020c, ), dynamic inter-regional transport
flow data (Chen and Chen,, 2019), and sequential image and video processing
(Walden and Serroukh,, 2002), among many others.

To model temporal dependencies of tensor-valued time series data, naturally
one might resort to vectorization of the tensor-valued observations so that
conventional vector-valued time series models can be directly applied. For
instance, let
$\mbox{\boldmath$\mathscr{Y}$}_{t}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$ be the tensor-valued observation at time $t$, for $t=1,\dots,T$, where
$T$ is the sample size. Then an obvious first step is to conduct the following
vector autoregression (VAR) for its vectorization:

$\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})=\boldsymbol{A}\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})+\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t}),$
(1)

where $\boldsymbol{A}\in\mathbb{R}^{p\times p}$ is the unknown transition
matrix, with $p=\prod_{i=1}^{d}p_{i}$, and
$\mbox{\boldmath$\mathscr{E}$}_{t}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$ is the tensor-valued random error at time $t$. Although model (1) is
potent in that it takes into account the linear association between every
variable in $\mbox{\boldmath$\mathscr{Y}$}_{t}$ and that in
$\mbox{\boldmath$\mathscr{Y}$}_{t-1}$, it clearly suffers from two major
drawbacks:

  * •

The vectorization will destroy the important multidimensional structural
information inherently embedded in the tensor-valued observations, resulting
in a lack of interpretability;

  * •

The number of unknown parameters, $p^{2}=(\prod_{i=1}^{d}p_{i})^{2}$, can be
formidable even for small $d$ and $p_{i}$; e.g., even when $d=3$ and
$p_{1}=p_{2}=p_{3}=3$, the number of unknown parameters will be as large as
$(3\times 3\times 3)^{2}=729$, while the sample size $T$ often has a similar
magnitude in practice; see, e.g., the real data examples in Section 7.

In this work, motivated by the idea of Tucker decomposition (Tucker,, 1966),
we propose the Low-Rank Tensor Autoregressive (LRTAR) model through folding
the $p\times p$ transition matrix $\boldsymbol{A}$ in (1) into the $2d$-th-
order transition tensor
$\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}\times p_{1}\times\cdots\times p_{d}}$ which is assumed to have low
multilinear ranks $(r_{1},\dots,r_{2d})$. That is, $r_{i}$ can be much smaller
than $p_{i}$, where $p_{d+i}=p_{i}$ for $i=1,\dots,d$. As a result, the
aforementioned drawbacks of the vectorization approach can be overcome by the
proposed model:

  * •

As we will show in Section 3, the proposed model can flexibly retain the
distinct structural information across all modes of observed tensor process,
encapsulating an interpretable low-dimensional tensor dynamic relationship
between $\mbox{\boldmath$\mathscr{Y}$}_{t}$ and
$\mbox{\boldmath$\mathscr{Y}$}_{t-1}$.

  * •

As the transition tensor $\mathscr{A}$ is assumed to have low multilinear
ranks, the parameter space is simultaneously constrained in $2d$ directions,
reducing its dimension from $(\prod_{i=1}^{d}p_{i})^{2}$ drastically to
$\prod_{i=1}^{2d}r_{i}+\sum_{i=1}^{d}r_{i}(p_{i}-r_{i})+\sum_{i=1}^{d}r_{d+i}(p_{i}-r_{d+i})$.

Recently there has been a rapidly emerging interest in high-dimensional
matrix- and tensor-valued time series analysis, particularly through factor
models, e.g., the matrix factor model proposed in Wang et al., (2019), the
constrained matrix factor model in Chen et al., 2020a , and the tensor factor
model in Chen et al., 2020c . Factor models are powerful tools for dimension
reduction with great interpretability. However, unlike autoregression, factor
models do not seek to explicitly model the temporal dependency and thus by
themselves cannot be directly used for forecasting. On the other hand, despite
the extensive literature on high-dimensional VAR models in recent years (e.g.
Negahban and Wainwright,, 2011; Basu and Michailidis,, 2015; Han et al.,
2015a, ; Wang et al.,, 2020; Zheng and Cheng,, 2020), counterparts able to
meet the particular needs of matrix- and tensor-valued time series modeling
tasks have been rarely explored. The most relevant existing work in this
direction so far might be the matrix autoregressive (MAR) model in Chen et
al., 2020b , where the focus is on low-dimensional modeling; see also Hoff,
(2015). As we will discuss in Section 3, the proposed model includes the MAR
model as a special case, but enjoys greater flexibility and a more drastic
reduction of the dimensionality; see also Figure 1 for an illustration.

The estimation of the proposed model is studied under both low- and high-
dimensional scaling. When the sample size $T$ is sufficiently large and the
transition tensor $\mathscr{A}$ is exactly low-rank, we consider the rank-
constrained estimation method and prove the asymptotic normality for the
proposed low-Tucker-rank (LTR) estimator. Under the high-dimensional setup, we
consider a more general and natural setting where $\mathscr{A}$ can be well
approximated by a low-Tucker-rank tensor, and develop regularized estimation
methods based on nuclear-norm-type penalties. We first study the Sum of
Nuclear (SN) norm regularizer, defined as the sum of nuclear norms of all one-
mode matricizations of $\mathscr{A}$, and derive the non-asymptotic estimation
error bound for the corresponding estimator. The SN norm regularizer has been
widely used in the literature for various low-rank tensor problems (Gandy et
al.,, 2011; Tomioka et al.,, 2011; Liu et al.,, 2013; Raskutti et al.,, 2019).
Its major strength lies in the fact that the summation of nuclear norms allows
enforcing the low-rankness simultaneously across all modes of the tensor. In
contrast, if only a single nuclear norm is used, the low-rankness of only one
mode will be accounted for, obviously leading to a much less efficient
estimator.

However, although penalizing multiple one-mode matricizations of $\mathscr{A}$
simultaneously is far more efficient than penalizing only one of them, the SN
regularized estimator suffers from serious efficiency loss due to the fat-and-
short shape of the one-mode matricizations. Note that the low-rankness in fact
can also be encouraged by penalizing nuclear norms of multi-mode
matricizations. For instance, the conventional Matrix Nuclear (MN) norm
regularized estimator (Negahban and Wainwright,, 2011) simply penalizes the
nuclear norm of the transition matrix $\boldsymbol{A}$ in representation (1),
which under the proposed LRTAR model actually is a square-shaped multi-mode
matricization, namely square matricization, of the transition tensor
$\mathscr{A}$. As we will show in Theorem 3 and simulations, even though the
MN regularizer incorporates only one single square matricization, it still
clearly beats the SN regularizer, since the former avoids the efficiency
bottleneck caused by the imbalance of the one-mode matricization.

Indeed, due to the autoregressive form of the proposed model, the transition
tensor $\mathscr{A}$ has a special balanced structure in the sense that its
first $d$ modes have exactly the same dimensions as its other $d$ modes. As a
result, actually many different square matricizations of $\mathscr{A}$ can be
formed by appropriately pairing up its modes; see Section 5.2 for details.
This naturally motivates us to propose a new regularizer that combines the
strengths of both SN and MN norms. Specifically, for the proposed tensor
autoregression, we introduce a novel Sum of Square-matrix Nuclear (SSN) norm
regularizer, defined as the sum of nuclear norms of all the $p\times p$ square
matricizations of the transition tensor $\mathscr{A}$. The SSN regularizer is
expected to be superior to the MN, since it simultaneously encourages the low-
rankness across all possible square matricizations of $\mathscr{A}$ rather
than only one of them. Moreover, thanks to the use of square matricizations,
the SSN regularized estimator is provably more efficient than the SN
regularized one; see Theorem 4 and the last simulation experiment in Section
6. Note that the adoption of a more balanced (square) matricization to improve
estimation performance was proposed in Mu et al., (2014) for low-rank tensor
completion problems, where only a single square matricization was penalized,
similarly to the MN regularizer.

This work is also related to the literature of matrix-variate regression and
tensor regression for independent data. The matrix-variate regression model in
Ding and Cook, (2018) has the same basic bilinear form as that of the MAR
model mentioned earlier, while an envelope method was introduced to further
reduce the dimension. Raskutti et al., (2019) proposed a multi-response tensor
regression model, where they mainly studied the third-order coefficient tensor
and the SN regularization which is known to be statistically sub-optimal for
higher-order tensor estimation. In contrast, we study the model for general
higher-order tensor-valued time series. Moreover, our SSN regularized
estimator has a much faster statistical convergence rate than the SN
estimator. Recently, Chen et al., (2019) and Han et al., (2020) studied non-
convex projected gradient descent methods for tensor regression. While their
non-convex approaches require exact low-rankness with known Tucker ranks, our
methods can handle both exact and approximate low-rankness and select the
unknown ranks consistently in the exactly low-rank case. In addition, existing
literature on tensor regression has only considered independent data and
Gaussian time series data, whereas we allow sub-Gaussianity of the time
series. This is a non-trivial relaxation, since unlike the Gaussian case, sub-
Gaussian time series cannot be linearly transformed into independent samples.

We summarize the most important contributions of this paper as follows:

  1. (i)

This paper provides the first practical tool to model and forecast general
structured, high-dimensional data with complex temporal dependence via tensor
autoregression. By flexibly and efficiently capturing the underlying low-
dimensional tensor dynamics, the proposed model delivers significant dimension
reduction, meaningful structural interpretation and favorable forecast
performance.

  2. (ii)

By exploiting the special balanced structure of the transition tensor
$\mathscr{A}$, a novel SSN regularization approach is introduced to
simultaneously and efficiently encourage low-rankness across all square
matricizations of $\mathscr{A}$, outperforming both the SN and MN methods
under both exact and approximate low-rankness. For exactly low-rank
$\mathscr{A}$, a truncated estimator is further introduced for consistent rank
selection.

  3. (iii)

On the technical front, by establishing a novel martingale-based concentration
bound, this paper relaxes the conventional Gaussian assumption in the
literature on high-dimensional time series to sub-Gaussianity. This technique
is generally applicable to the non-asymptotic estimation theory for high-
dimensional time series models with a VAR representation and hence is of
independent interest.

The rest of the paper is organized as follows. Section 2 introduces basic
notation and tensor algebra. Section 3 presents the proposed LRTAR model.
Section 4 studies the low-dimensional least squares estimator and its
asymptotic properties. The high-dimensional regularized estimation is covered
in Section 5, where we develop the non-asymptotic theory for three regularized
estimators and rank selection consistency for the truncated estimator.
Sections 6 and 7 present simulation studies and real data analysis,
respectively. Section 8 concludes with a brief discussion. All technical
proofs are relegated to the Appendix.

##  2 Preliminaries: Notation and Tensor Algebra

Tensors, also known as multidimensional arrays, are natural higher-order
extensions of matrices. The order of a tensor is known as the dimension, way
or mode, so a multidimensional array
$\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$
is called a $d$-th-order tensor. We refer readers to Kolda and Bader, (2009)
for a detailed review of basic tensor algebra.

Throughout this paper, we denote vectors by boldface small letters, e.g.
$\boldsymbol{x}$, $\boldsymbol{y}$, matrices by boldface capital letters, e.g.
$\boldsymbol{X}$, $\boldsymbol{Y}$, and tensors by boldface Euler capital
letters, e.g. $\mathscr{X}$, $\mathscr{Y}$. For any two real-valued sequences
$x_{k}$ and $y_{k}$, we write $x_{k}\gtrsim y_{k}$ if there exists a constant
$c>0$ such that $x_{k}\geq cy_{k}$ for all $k$, and write $x_{k}\gg y_{k}$ if
$\lim_{k\rightarrow\infty}y_{k}/x_{k}=0$. In addition, write $x_{k}\asymp
y_{k}$ if $x_{k}\gtrsim y_{k}$ and $y_{k}\gtrsim x_{k}$. We use $C$ to denote
a generic positive constant, which is independent of the dimensions and the
sample size.

For a generic matrix $\boldsymbol{X}$, we let $\boldsymbol{X}^{\top}$,
$\|\boldsymbol{X}\|_{\text{F}}$, $\|\boldsymbol{X}\|_{\text{op}}$,
$\|\boldsymbol{X}\|_{*}$, $\text{vec}(\boldsymbol{X})$, and
$\sigma_{j}(\boldsymbol{X})$ denote its transpose, Frobenius norm, operator
norm, nuclear norm, vectorization, and $j$-th largest singular value,
respectively. For any matrix $\boldsymbol{X}\in\mathbb{R}^{m\times n}$, recall
that the nuclear norm and its dual norm, the operator norm, are defined as

$\|\boldsymbol{X}\|_{*}=\sum_{j=1}^{\min(m,n)}\sigma_{j}(\boldsymbol{X})~{}~{}\text{and}~{}~{}\|\boldsymbol{X}\|_{\text{op}}=\sigma_{1}(\boldsymbol{X}).$
(2)

For any square matrix $\boldsymbol{X}$, we let
$\lambda_{\min}(\boldsymbol{X})$ and $\lambda_{\max}(\boldsymbol{X})$ denote
its minimum and maximum eigenvalues. For any real symmetric matrices
$\boldsymbol{X}$ and $\boldsymbol{Y}$, we write
$\boldsymbol{X}\leq\boldsymbol{Y}$ if $\boldsymbol{Y}-\boldsymbol{X}$ is a
positive semidefinite matrix.

Matricization, also known as unfolding, is the process of reordering the
elements of a third- or higher-order tensor into a matrix. The most commonly
used matricization is the one-mode matricization defined as follows. For any
$d$-th-order tensor
$\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$,
its mode-$s$ matricization
$\mbox{\boldmath$\mathscr{X}$}_{(s)}\in\mathbb{R}^{p_{s}\times p_{-s}}$, with
$p_{-s}=\prod_{i=1,i\neq s}^{d}p_{i}$, is the matrix obtained by setting the
$s$-th tensor mode as its rows and collapsing all the others into its columns,
for $s=1,\dots,d$. Specifically, the $(i_{1},\dots,i_{d})$-th element of
$\mathscr{X}$ is mapped to the $(i_{s},j)$-th element of
$\mbox{\boldmath$\mathscr{X}$}_{(s)}$, where

$j=1+\sum_{\begin{subarray}{c}k=1\\\ k\neq
s\end{subarray}}^{d}(i_{k}-1)J_{k}~{}~{}\text{with}~{}~{}J_{k}=\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq s\end{subarray}}^{k-1}p_{\ell}.$ (3)

The above one-mode matricization can be extended to the multi-mode
matricization by combining multiple modes to rows and combining the rest to
columns of a matrix. For any index subset $S\subset\\{1,\dots,d\\}$, the
multi-mode matricization $\mbox{\boldmath$\mathscr{X}$}_{[S]}$ is the
$\prod_{i\in S}p_{i}$-by-$\prod_{i\notin S}p_{i}$ matrix whose $(i,j)$-th
element is mapped from the $(i_{1},\dots,i_{d})$-th element of $\mathscr{X}$,
where

$i=1+\sum_{k\in S}(i_{k}-1)I_{k}~{}\text{and}~{}j=1+\sum_{k\notin
S}(i_{k}-1)J_{k},~{}\text{with}~{}I_{k}=\prod_{\begin{subarray}{c}\ell\in S\\\
\ell<k\end{subarray}}p_{\ell}~{}\text{and}~{}J_{k}=\prod_{\begin{subarray}{c}\ell\notin
S\\\ \ell<k\end{subarray}}p_{\ell}.$ (4)

Note that the modes in the multi-mode matricization are collapsed following
their original order $1,\dots,d$. Moreover, it holds
$\mbox{\boldmath$\mathscr{X}$}_{[S]}=\mbox{\boldmath$\mathscr{X}$}_{[S^{\complement}]}^{\top}$,
where $S^{\complement}=\\{1,\dots,d\\}\setminus S$ is the complement of $S$.
In addition, the one-mode matricization $\mbox{\boldmath$\mathscr{X}$}_{(s)}$
defined above is simply $\mbox{\boldmath$\mathscr{X}$}_{[\\{s\\}]}$.

We next review the concepts of tensor-matrix multiplication, tensor
generalized inner product and norm. For any $d$-th-order tensor
$\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$
and matrix $\boldsymbol{Y}\in\mathbb{R}^{q_{k}\times p_{k}}$ with $1\leq k\leq
d$, the mode-$k$ multiplication
$\mbox{\boldmath$\mathscr{X}$}\times_{k}\boldsymbol{Y}$ produces a $d$-th-
order tensor in $\mathbb{R}^{p_{1}\times\cdots\times p_{k-1}\times q_{k}\times
p_{k+1}\times\cdots\times p_{d}}$ defined by

$\left(\mbox{\boldmath$\mathscr{X}$}\times_{k}\boldsymbol{Y}\right)_{i_{1}\cdots
i_{k-1}ji_{k+1}\dots
i_{d}}=\sum_{i_{k}=1}^{p_{k}}\mbox{\boldmath$\mathscr{X}$}_{i_{1}\cdots
i_{d}}\boldsymbol{Y}_{ji_{k}}.$ (5)

For any two tensors $\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times
p_{2}\times\cdots\times p_{n}}$ and
$\mbox{\boldmath$\mathscr{Y}$}\in\mathbb{R}^{p_{1}\times p_{2}\times\cdots
p_{m}}$ with $n\geq m$, their generalized inner product
$\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{Y}$}\rangle$ is
the $(n-m)$-th-order tensor in $\mathbb{R}^{p_{m+1}\times\dots\times p_{n}}$
defined by

$\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{Y}$}\rangle_{i_{m+1}\dots
i_{n}}=\sum_{i_{1}=1}^{p_{1}}\sum_{i_{2}=1}^{p_{2}}\dots\sum_{i_{m}=1}^{p_{m}}\mbox{\boldmath$\mathscr{X}$}_{i_{1}i_{2}\dots
i_{m}i_{m+1}\dots i_{n}}\mbox{\boldmath$\mathscr{Y}$}_{i_{1}i_{2}\dots
i_{m}},$ (6)

where $1\leq i_{m+1}\leq p_{m+1},\dots,1\leq i_{n}\leq p_{n}$. In particular,
when $n=m$, it reduces to the conventional real-valued inner product. In
addition, the Frobenius norm of any tensor $\mathscr{X}$ is defined as
$\|\mbox{\boldmath$\mathscr{X}$}\|_{\text{F}}=\sqrt{\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{X}$}\rangle}$.

Some basic properties of the tensor generalized inner product are as follows.
Let $\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times
p_{2}\times\cdots\times p_{n}}$,
$\mbox{\boldmath$\mathscr{Y}$}\in\mathbb{R}^{p_{1}\times p_{2}\times\cdots
p_{m}}$ and
$\mbox{\boldmath$\mathscr{Z}$}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{k-1}\times q_{k}\times p_{k+1}\cdots\times p_{m}}$ be tensors with $n\geq
m\geq k\geq 1$. If $\boldsymbol{Y}\in\mathbb{R}^{q_{k}\times p_{k}}$, then

$\langle\mbox{\boldmath$\mathscr{X}$}\times_{k}\boldsymbol{Y},\mbox{\boldmath$\mathscr{Z}$}\rangle=\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{Z}$}\times_{k}\boldsymbol{Y}^{\top}\rangle.$
(7)

If $\boldsymbol{Z}\in\mathbb{R}^{q_{m+j}\times p_{m+j}}$ with $1\leq j\leq
n-m$, then

$\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{Y}$}\rangle\times_{j}\boldsymbol{Z}=\langle\mbox{\boldmath$\mathscr{X}$}\times_{m+j}\boldsymbol{Z},\mbox{\boldmath$\mathscr{Y}$}\rangle.$
(8)

Moreover,

$\text{vec}(\langle\mbox{\boldmath$\mathscr{X}$},\mbox{\boldmath$\mathscr{Y}$}\rangle)=\mbox{\boldmath$\mathscr{X}$}_{[S]}\text{vec}(\mbox{\boldmath$\mathscr{Y}$}),$
(9)

where $S=\\{m+1,\dots,n\\}$, and when $m=n$,
$\mbox{\boldmath$\mathscr{X}$}_{[\emptyset]}=\text{vec}(\mbox{\boldmath$\mathscr{X}$})^{\top}$.

Finally, we summarize some concepts and useful results of the Tucker
decomposition (Tucker,, 1966; De Lathauwer et al.,, 2000). For any tensor
$\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$,
its multilinear ranks $(r_{1},\dots,r_{d})$ are defined as the matrix ranks of
its one-mode matricizations, namely
$r_{i}=\text{rank}(\mbox{\boldmath$\mathscr{X}$}_{(i)})$, for $i=1,\dots,d$.
Note that $r_{i}$’s are analogous to the row and column ranks of a matrix, but
are not necessarily equal for third- and higher-order tensors. Suppose that
$\mathscr{X}$ has multilinear ranks $(r_{1},\dots,r_{d})$. Then $\mathscr{X}$
has the following Tucker decomposition:

$\mbox{\boldmath$\mathscr{X}$}=\mbox{\boldmath$\mathscr{Y}$}\times_{1}\boldsymbol{Y}_{1}\times_{2}\boldsymbol{Y}_{2}\cdots\times_{d}\boldsymbol{Y}_{d}=\mbox{\boldmath$\mathscr{Y}$}\times_{i=1}^{d}\boldsymbol{Y}_{i},$
(10)

where $\boldsymbol{Y}_{i}\in\mathbb{R}^{p_{i}\times r_{i}}$ for $i=1,\dots,d$
are the factor matrices and
$\mbox{\boldmath$\mathscr{Y}$}\in\mathbb{R}^{r_{1}\times\cdots\times r_{d}}$
is the core tensor. Hence, the multilinear ranks are also called Tucker ranks.

If $\mathscr{X}$ has the Tucker decomposition in (10), then we have the
following results for its one- and multi-mode matricizations:

$\mbox{\boldmath$\mathscr{X}$}_{(s)}=\boldsymbol{Y}_{s}\mbox{\boldmath$\mathscr{Y}$}_{(s)}(\boldsymbol{Y}_{d}\otimes\cdots\otimes\boldsymbol{Y}_{s+1}\otimes\boldsymbol{Y}_{s-1}\cdots\otimes\boldsymbol{Y}_{1})^{\top}=\boldsymbol{Y}_{s}\mbox{\boldmath$\mathscr{Y}$}_{(s)}(\otimes_{i\neq
s}\boldsymbol{Y}_{i})^{\top},\quad s=1,\dots,d,$ (11)

and

$\mbox{\boldmath$\mathscr{X}$}_{[S]}=(\otimes_{i\in
S}\boldsymbol{Y}_{i})\mbox{\boldmath$\mathscr{Y}$}_{[S]}(\otimes_{i\notin
S}\boldsymbol{Y}_{i})^{\top},\quad S\subset\\{1,\dots,d\\},$ (12)

where $\otimes_{i\neq s},\otimes_{i\in S}$ and $\otimes_{i\notin S}$ are
matrix Kronecker products operating in the reverse order within the
corresponding index sets.

Note that for any nonsingular matrices
$\boldsymbol{O}_{i}\in\mathbb{R}^{r_{i}\times r_{1}}$ for $i=1,\dots,d$, it
holds

$\mbox{\boldmath$\mathscr{Y}$}\times_{i=1}^{d}\boldsymbol{Y}_{i}=(\mbox{\boldmath$\mathscr{Y}$}\times_{i=1}^{d}\boldsymbol{O}_{i})\times_{i=1}^{d}(\boldsymbol{Y}_{i}\boldsymbol{O}_{i}^{-1}).$

This indicates that the Tucker decomposition in (10) is not unique unless
appropriate identifiability constraints are imposed. In this paper, to fix
ideas, we will focus on a special Tucker decomposition called the higher-order
singular value decomposition (HOSVD). In the HOSVD, the factor matrix
$\boldsymbol{Y}_{i}$ in (10) is defined as the tall orthonormal matrix
consisting of the top $r_{i}$ left singular vectors of
$\mbox{\boldmath$\mathscr{X}$}_{(i)}$, for $i=1,\dots,d$. Consequently, the
core tensor
$\mbox{\boldmath$\mathscr{Y}$}=\mbox{\boldmath$\mathscr{X}$}\times_{i=1}^{d}\boldsymbol{Y}_{i}^{\top}$
has the all-orthogonal property that
$\mbox{\boldmath$\mathscr{Y}$}_{(i)}\mbox{\boldmath$\mathscr{Y}$}_{(i)}^{\top}$
is a diagonal matrix for $i=1,\dots,d$.

##  3 Low-Rank Tensor Autoregression

For the tensor-valued time series
$\\{\mbox{\boldmath$\mathscr{Y}$}_{t}\\}_{t=1}^{T}$, we propose the following
Low-Rank Tensor Autoregressive (LRTAR) model:

$\mbox{\boldmath$\mathscr{Y}$}_{t}=\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle+\mbox{\boldmath$\mathscr{E}$}_{t},$
(13)

where $\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}\times p_{1}\times\cdots\times p_{d}}$ is the $2d$-th-order transition
tensor which is assumed to have multilinear ranks $(r_{1},\dots,r_{2d})$, with
$r_{i}=\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{(i)})$,
$\langle\cdot,\cdot\rangle$ is the generalized tensor inner product defined in
(6), and
$\mbox{\boldmath$\mathscr{E}$}_{t}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$ is the mean-zero random error at time $t$ with possible dependencies
among its contemporaneous elements. Denote $S_{1}=\\{1,2,\dots,d\\}$ and
$S_{2}=\\{d+1,d+2,\dots,2d\\}$. Note that by (9), model (13) can be written
into the VAR form in (1) with transition matrix
$\boldsymbol{A}=\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$.

Then, we have the higher-order singular value decomposition (HOSVD) of
$\mathscr{A}$,

$\mbox{\boldmath$\mathscr{A}$}=\mbox{\boldmath$\mathscr{G}$}\times_{i=1}^{2d}\boldsymbol{U}_{i},$
(14)

where $\mbox{\boldmath$\mathscr{G}$}\in\mathbb{R}^{r_{1}\times\cdots\times
r_{2d}}$ is the core tensor, and each
$\boldsymbol{U}_{i}\in\mathbb{R}^{p_{i}\times r_{i}}$ is the orthonormal
factor matrix defined as the top $r_{i}$ left singular vectors of
$\mbox{\boldmath$\mathscr{A}$}_{(i)}$, for $1\leq i\leq 2d$. Thus, by (12),
the VAR representation of model (13) can be written as

$\underbrace{\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})}_{\Large{\boldsymbol{y}_{t}}}=\underbrace{(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top}}_{\Large{\normalsize{\mbox{\boldmath$\mathscr{A}$}}_{[S_{2}]}}}\underbrace{\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})}_{\Large{\boldsymbol{y}_{t-1}}}+\underbrace{\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})}_{\Large{\boldsymbol{e}_{t}}},$
(15)

where $\boldsymbol{y}_{t}=\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})$ and
$\boldsymbol{e}_{t}=\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})$.

In contrast to the conventional VAR model in (1) which has $p^{2}$ unknown
parameters, where $p=\prod_{i=1}^{d}p_{i}$, the dimension of the parameter
space for model (13) is reduced substantially to

$\prod_{i=1}^{2d}r_{i}+\sum_{i=1}^{d}r_{i}(p_{i}-r_{i})+\sum_{i=1}^{d}r_{d+i}(p_{i}-r_{d+i}),$
(16)

which is computed by subtracting the number of constraints due to the
orthonormality of $\boldsymbol{U}_{i}$’s and the all-orthogonal property of
$\mathscr{G}$ from the total number of parameters.

By the VAR representation in (15), we immediately have the necessary and
sufficient condition for the existence of a unique strictly stationary
solution to model (13) as follows.

######  Assumption 1.

The spectral radius of $\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$ is strictly
less than one.

The proposed model implies an interesting low-dimensional tensor dynamical
structure. To be specific, by (7), (8) and the orthonormality of
$\boldsymbol{U}_{i}$, it can be shown that (13) together with (14) implies

$\underbrace{\mbox{\boldmath$\mathscr{Y}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{d+i}^{\top}}_{r_{d+1}\times
r_{d+2}\times\cdots\times
r_{2d}}=\big{\langle}\mbox{\boldmath$\mathscr{G}$},~{}\underbrace{\mbox{\boldmath$\mathscr{Y}$}_{t-1}\times_{i=1}^{d}\boldsymbol{U}_{i}^{\top}}_{r_{1}\times
r_{2}\times\cdots\times
r_{d}}\big{\rangle}+\mbox{\boldmath$\mathscr{E}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{d+i}^{\top}.$
(17)

Note that in (17), $\mbox{\boldmath$\mathscr{Y}$}_{t}$ and
$\mbox{\boldmath$\mathscr{E}$}_{t}$ are both projected onto a low-dimensional
space via the $\boldsymbol{U}_{d+i}$’s, while
$\mbox{\boldmath$\mathscr{Y}$}_{t-1}$ is projected onto another low-
dimensional space via the $\boldsymbol{U}_{i}$’s, with $1\leq i\leq d$. Hence,
(17) is essentially a low-dimensional tensor regression defined on the
projections of $\mbox{\boldmath$\mathscr{Y}$}_{t}$ and
$\mbox{\boldmath$\mathscr{Y}$}_{t-1}$. Element-wisely, the low-dimensional
tensor
$\mbox{\boldmath$\mathscr{Y}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{d+i}^{\top}$
can be interpreted as $\prod_{i=1}^{d}r_{d+i}$ multilinear response factors,
$\mbox{\boldmath$\mathscr{Y}$}_{t-1}\times_{i=1}^{d}\boldsymbol{U}_{i}^{\top}$
as $\prod_{i=1}^{d}r_{i}$ multilinear predictor factors, and
$\mbox{\boldmath$\mathscr{E}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{d+i}^{\top}$
as multilinear error factors. For this reason, we call
$\boldsymbol{U}_{d+1},\dots,\boldsymbol{U}_{2d}$ the response factor matrices,
and $\boldsymbol{U}_{1},\dots,\boldsymbol{U}_{d}$ the predictor factor
matrices.

We discuss some special cases of the proposed model below.

######  Example 1.

For simplicity, we first consider the case with $d=2$, so
$\mbox{\boldmath$\mathscr{Y}$}_{t}\equiv\boldsymbol{Y}_{t},\mbox{\boldmath$\mathscr{E}$}_{t}\equiv\boldsymbol{E}_{t}\in\mathbb{R}^{p_{1}\times
p_{2}}$ are matrices. Then the VAR representation in (15) becomes

$\textup{vec}(\boldsymbol{Y}_{t})=(\boldsymbol{U}_{4}\otimes\boldsymbol{U}_{3})\mbox{\boldmath$\mathscr{G}$}_{[\\{3,4\\}]}(\boldsymbol{U}_{2}^{\top}\otimes\boldsymbol{U}_{1}^{\top})\textup{vec}(\boldsymbol{Y}_{t-1})+\textup{vec}(\boldsymbol{E}_{t}),$
(18)

and the low-dimensional representation in (17) becomes

$\boldsymbol{U}_{3}^{\top}\boldsymbol{Y}_{t}\boldsymbol{U}_{4}=\left\langle\mbox{\boldmath$\mathscr{G}$},\boldsymbol{U}_{1}^{\top}\boldsymbol{Y}_{t-1}\boldsymbol{U}_{2}\right\rangle+\boldsymbol{U}_{3}^{\top}\boldsymbol{E}_{t}\boldsymbol{U}_{4},$

where $\mbox{\boldmath$\mathscr{G}$}\in\mathbb{R}^{r_{1}\times\cdots\times
r_{4}}$. It is interesting to compare this model with the matrix
autoregressive (MAR) model in Chen et al., 2020b  and Hoff, (2015), which is
defined by

$\boldsymbol{Y}_{t}=\boldsymbol{B}_{1}\boldsymbol{Y}_{t-1}\boldsymbol{B}_{2}^{\top}+\boldsymbol{E}_{t},$
(19)

where $\boldsymbol{B}_{1}\in\mathbb{R}^{p_{1}\times p_{1}}$ and
$\boldsymbol{B}_{2}\in\mathbb{R}^{p_{2}\times p_{2}}$, whose vector form is

$\textup{vec}(\boldsymbol{Y}_{t})=(\boldsymbol{B}_{2}\otimes\boldsymbol{B}_{1})\textup{vec}(\boldsymbol{Y}_{t-1})+\textup{vec}(\boldsymbol{E}_{t}).$
(20)

It can be easily seen that if $r_{1}=r_{3}=p_{1}$, $r_{2}=r_{4}=p_{2}$,
$\boldsymbol{U}_{3}=\boldsymbol{I}_{p_{1}}$,
$\boldsymbol{U}_{4}=\boldsymbol{I}_{p_{2}}$, and
$\mbox{\boldmath$\mathscr{G}$}_{[\\{3,4\\}]}=(\boldsymbol{B}_{2}\otimes\boldsymbol{B}_{1})(\boldsymbol{U}_{2}\otimes\boldsymbol{U}_{1})$,
then (18) becomes exactly (20). Thus, the MAR model in (19) can be viewed as a
special case of the proposed model without reducing dimensions $p_{i}$’s to
$r_{i}$’s and without transforming $\boldsymbol{Y}_{t}$; see Figure 1 for an
illustration. The above comparison also applies to the general case with
$d\geq 3$. The tensor version of the MAR model is considered in Hoff, (2015)
and is defined as

$\mbox{\boldmath$\mathscr{Y}$}_{t}=\mbox{\boldmath$\mathscr{Y}$}_{t-1}\times_{i=1}^{d}\boldsymbol{B}_{i}+\mbox{\boldmath$\mathscr{E}$}_{t},$
(21)

where $\boldsymbol{B}_{i}\in\mathbb{R}^{p_{i}\times p_{i}}$ for $i=1,\dots,d$.
We call (21) the multilinear tensor autoregressive (MTAR) model. Note that its
vector form is

$\textup{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})=(\boldsymbol{B}_{d}\otimes\cdots\otimes\boldsymbol{B}_{1})\textup{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})+\textup{vec}(\mbox{\boldmath$\mathscr{E}$}_{t}).$
(22)

Similarly, (22) is a special case of (15) with $r_{i}=r_{d+i}=p_{i}$,
$\boldsymbol{U}_{d+i}=\boldsymbol{I}_{p_{i}}$, for $i=1,\dots,d$, and
$\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}=(\otimes_{i\in
S_{1}}\boldsymbol{B}_{i})(\otimes_{i\in S_{1}}\boldsymbol{U}_{i})$. Obviously,
the number of unknown parameters in the MTAR model, $\sum_{i=1}^{d}p_{i}^{2}$,
is much larger than that of the proposed model as shown in (16). Also note
that Chen et al., 2020b  focuses on the low-dimensional estimation and its
asymptotic theory, while Hoff, (2015) considers a Bayesian estimation method.

Figure 1: Illustration of the MAR model and the proposed LRTAR model in the
case of $d=2$.

######  Example 2.

In the special case where $\boldsymbol{U}_{d+i}=\boldsymbol{U}_{i}$ and
$r_{d+i}=r_{i}$ for $i=1,\dots,d$, the proposed model may be understood from
the perspective of dynamic factor modeling (Stock and Watson,, 2011; Bai and
Wang,, 2016) for tensor-valued time series. Specifically, consider the
following model:

$\mbox{\boldmath$\mathscr{Y}$}_{t}=\mbox{\boldmath$\mathscr{F}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{i},\quad\mbox{\boldmath$\mathscr{F}$}_{t}=\langle\mbox{\boldmath$\mathscr{G}$},\mbox{\boldmath$\mathscr{F}$}_{t-1}\rangle+\mbox{\boldmath$\mathscr{H}$}_{t},$
(23)

where $\mbox{\boldmath$\mathscr{Y}$}_{t}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$ is the observed tensor-valued time series,
$\mbox{\boldmath$\mathscr{F}$}_{t}\in\mathbb{R}^{r_{1}\times\cdots\times
r_{d}}$ represents $\prod_{i=1}^{d}r_{i}$ factors, and
$\boldsymbol{U}_{i}\in\mathbb{R}^{p_{i}\times r_{i}}$ are orthonormal matrices
for $i=1,\dots,d$. Here $\mbox{\boldmath$\mathscr{F}$}_{t}$ follows the tensor
autoregression (TAR) with transition tensor
$\mbox{\boldmath$\mathscr{G}$}\in\mathbb{R}^{r_{1}\times\cdots\times
r_{d}\times r_{1}\times\cdots\times r_{d}}$ and random error
$\mbox{\boldmath$\mathscr{H}$}_{t}$. Note that (23) can be rewritten as

$\mbox{\boldmath$\mathscr{Y}$}_{t}=\left\langle\mbox{\boldmath$\mathscr{G}$}\times_{i=1}^{d}\boldsymbol{U}_{i}\times_{i=d+1}^{2d}\boldsymbol{U}_{i},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\right\rangle+\mbox{\boldmath$\mathscr{H}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{i}.$

Thus, model (23) is a special case of the proposed model with
$\boldsymbol{U}_{d+i}=\boldsymbol{U}_{i}$ and $r_{d+i}=r_{i}$ for
$i=1,\dots,d$, and
$\mbox{\boldmath$\mathscr{E}$}_{t}=\mbox{\boldmath$\mathscr{H}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{i}$.
Chen et al., 2020c  introduces the tensor factor model in the form of
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\mbox{\boldmath$\mathscr{F}$}_{t}\times_{i=1}^{d}\boldsymbol{U}_{i}+\mbox{\boldmath$\mathscr{E}$}_{t}$
without an explicit modeling of the latent factors
$\mbox{\boldmath$\mathscr{F}$}_{t}$. Hence, model (23) may be regarded as a
special tensor factor model with autoregressive dynamic factors, but without
any random error in the model equation of $\mbox{\boldmath$\mathscr{Y}$}_{t}$.

##  4 Low-Dimensional Least Squares Estimation

We consider the parameter estimation for the low-dimensional case where the
sample size $T$ is sufficiently large such that the dimension of the parameter
space can be assumed fixed. Throughout this section, we assume that the data
are generated from the proposed model in Section 3, where the true transition
tensor $\mathscr{A}$ is exactly low-Tucker rank with multilinear ranks
$(r_{1},\dots,r_{2d})$. This assumption will be relaxed under the high-
dimensional setup in the next section.

Suppose that the true ranks $(r_{1},\dots,r_{2d})$ of the exactly low-rank
tensor $\mathscr{A}$ are known; we relegate the rank selection to Section 5.3.
Then the parameters can be estimated by

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}=\mbox{\boldmath$\mathscr{\widehat{G}}$}\times_{i=1}^{2d}\boldsymbol{\widehat{U}}_{i}=\underset{\scalebox{0.7}{\mbox{\boldmath$\mathscr{G}$}},\boldsymbol{U}_{i}}{\arg\min}\sum_{t=1}^{T}\left\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{G}$}\times_{i=1}^{2d}\boldsymbol{U}_{i},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\right\|_{\textup{F}}^{2}.$
(24)

We call $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$ the low-Tucker-
rank (LTR) estimator. Note that the minimization in (24) is unconstrained, so
the Tucker decomposition of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$ is not unique. Indeed,
there are more than one solution of $\mathscr{\widehat{G}}$ and
$\boldsymbol{\widehat{U}}_{i}$’s corresponding to the same
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$. Due to the lack of
identifiability of the Tucker decompositions, standard asymptotics of the
maximum likelihood estimation cannot apply directly. Nevertheless, we can
still derive the asymptotic distribution of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$ using the asymptotic
theory for overparameterized models in Shapiro, (1986).

Recall that $S_{1}=\\{1,2,\dots,d\\}$ $S_{2}=\\{d+1,d+2,\dots,2d\\}$,
$\boldsymbol{y}_{t}=\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})$, and
$\boldsymbol{e}_{t}=\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})$. Let
$\boldsymbol{\theta}=(\text{vec}(\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]})^{\top},\text{vec}(\boldsymbol{U}_{1})^{\top},\cdots,\text{vec}(\boldsymbol{U}_{2d})^{\top})^{\top}$
be the parameter vector, and let
$\boldsymbol{h}=\boldsymbol{h}(\boldsymbol{\theta})=\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})=\text{vec}((\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top})$ be the vectorization of the transition
matrix. Denote
$\boldsymbol{\Sigma}_{\boldsymbol{y}}=\text{var}(\boldsymbol{y}_{t})$,
$\boldsymbol{\Sigma}_{\boldsymbol{e}}=\text{var}(\boldsymbol{e}_{t})$, and
$\boldsymbol{J}=\boldsymbol{\Sigma}_{\boldsymbol{e}}^{-1}\otimes\boldsymbol{\Sigma}_{\boldsymbol{y}}$.
In addition, for $i=1,\dots,2d$, let $\boldsymbol{P}_{[S_{2}]}^{(i)}$ be the
$p^{2}\times p^{2}$ permutation matrix such that
$\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})=\boldsymbol{P}_{[S_{2}]}^{(i)}\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{(i)})$.
Then, it can be shown that the Jacobian matrix
$\boldsymbol{H}:={\partial\boldsymbol{h}(\boldsymbol{\theta})}/{\partial\boldsymbol{\theta}}$
is given by

$\begin{split}\boldsymbol{H}=\big{(}&(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})\otimes(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i}),\boldsymbol{P}^{(1)}_{[S_{2}]}\left\\{\left[(\otimes_{i=1,i\neq
1}^{2d}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{(1)}^{\top}\right]\otimes\boldsymbol{I}_{p_{1}}\right\\},\\\
&\boldsymbol{P}^{(2)}_{[S_{2}]}\left\\{\left[(\otimes_{i=1,i\neq
2}^{2d}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{(2)}^{\top}\right]\otimes\boldsymbol{I}_{p_{2}}\right\\},\dots,\boldsymbol{P}^{(2d)}_{[S_{2}]}\left\\{\left[(\otimes_{i=1,i\neq
2d}^{2d}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{(2d)}^{\top}\right]\otimes\boldsymbol{I}_{p_{2d}}\right\\}\big{)},\end{split}$
(25)

where $p_{d+i}=p_{i}$ for $i=1,\dots,d$.

######  Theorem 1.

Suppose that the time series $\\{\mbox{\boldmath$\mathscr{Y}$}_{t}\\}$ is
generated by model (13) with
$\mathbb{E}\|\boldsymbol{e}_{t}\|_{2}^{4}<\infty$, and Assumption 1 holds.
Then,

$\sqrt{T}\textup{vec}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{LTR}}-\mbox{\boldmath$\mathscr{A}$})_{[S_{2}]})\rightarrow
N(\boldsymbol{0},\boldsymbol{\boldsymbol{\Sigma}}_{\textup{LTR}}),$ (26)

in distribution as $T\rightarrow\infty$, where
$\boldsymbol{\boldsymbol{\Sigma}}_{\textup{LTR}}=\boldsymbol{H}(\boldsymbol{H}^{\top}\boldsymbol{JH})^{\dagger}\boldsymbol{H}^{\top}$,
and $\dagger$ is the Moore-Penrose inverse.

Since the asymptotic theory for overparameterized models in Shapiro, (1986)
allows for unidentifiability of the components $\mathscr{G}$ and
$\boldsymbol{U}_{i}$’s in the decomposition of $\mathscr{A}$, Theorem 1 does
not require that the Tucker decomposition of $\mathscr{A}$ is unique; see the
proof of Theorem 1 in Appendix A for more details.

Next we compare the result in Theorem 1 to those of two other consistent
estimators for the proposed model in the low-dimensional setup. Note that the
rank of $\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$ in (15) is at most
$s_{1}:=\min(\prod_{i=1}^{d}r_{i},\prod_{i=d+1}^{2d}r_{i})$. Thus,
$\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$ can be estimated by both the
reduced-rank regression (RRR) and ordinary least squares (OLS) methods,

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{RRR}}=\underset{\text{rank}(\scalebox{0.7}{\mbox{\boldmath$\mathscr{B}$}}_{[S_{2}]})\leq
s_{1}}{\operatorname*{arg\,min}}~{}\frac{1}{T}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{B}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}$
(27)

and

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{OLS}}=\underset{\scalebox{0.7}{\mbox{\boldmath$\mathscr{B}$}}}{\operatorname*{arg\,min}}~{}\frac{1}{T}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{B}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}.$

Naturally, under model (13),
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$ is asymptotically more
efficient than $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{RRR}}$ and
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{OLS}}$:

######  Corollary 1.

If the conditions of Theorem 1 hold, then
$\sqrt{T}\textup{vec}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{OLS}}-\mbox{\boldmath$\mathscr{A}$})_{[S_{2}]})\to
N(0,\boldsymbol{\Sigma}_{\textup{OLS}})$ and
$\sqrt{T}\textup{vec}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{RRR}}-\mbox{\boldmath$\mathscr{A}$})_{[S_{2}]})\to
N(0,\boldsymbol{\Sigma}_{\textup{RRR}})$ in distribution as $T\to\infty$.
Moreover, it holds that
$\boldsymbol{\Sigma}_{\textup{LTR}}\leq\boldsymbol{\Sigma}_{\textup{RRR}}\leq\boldsymbol{\Sigma}_{\textup{OLS}}$.

To solve the minimization problem in (24), we propose an alternating least
squares (ALS) method. Specifically, by the vector representation of the
proposed model in (15), the loss function in (24) can be rewritten as

$L(\mbox{\boldmath$\mathscr{G}$},\boldsymbol{U}_{1},\dots,\boldsymbol{U}_{2d})=\sum_{t=1}^{T}\Big{\|}\boldsymbol{y}_{t}-(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top}\boldsymbol{y}_{t-1}\Big{\|}_{2}^{2}=\sum_{t=1}^{T}\Big{\|}\boldsymbol{y}_{t}-\boldsymbol{m}(\boldsymbol{\theta})\Big{\|}_{2}^{2},$

where $\boldsymbol{m}(\boldsymbol{\theta})=(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top}\boldsymbol{y}_{t-1}$ is the conditional mean
of $\boldsymbol{y}_{t}=\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t})$ and is a
function of the parameter vector
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{0}^{\top},\boldsymbol{\theta}_{1}^{\top},\dots,\boldsymbol{\theta}_{2d}^{\top})^{\top}$,
with
$\boldsymbol{\theta}_{0}=\text{vec}(\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]})$
and $\boldsymbol{\theta}_{i}=\text{vec}(\boldsymbol{U}_{i})$ for
$i=1,\dots,2d$. Note that $\boldsymbol{m}(\boldsymbol{\theta})$ is linear in
each $\boldsymbol{\theta}_{i}$ while keeping the other components
$\boldsymbol{\theta}_{j}$ with $0\leq j\neq i\leq 2d$ fixed. This is indeed
because

$\displaystyle\boldsymbol{m}(\boldsymbol{\theta})$
$\displaystyle=(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}\boldsymbol{P}_{k,1}\left\\{\left[(\otimes_{i\in
S_{1},i\neq
k}\boldsymbol{U}_{i})^{\top}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})_{(k)}^{\top}\right]\otimes\boldsymbol{I}_{r_{k}}\right\\}\boldsymbol{P}_{k,3}\text{vec}(\boldsymbol{U}_{k})$
(28) $\displaystyle=\boldsymbol{P}_{k,2}\left\\{\left[(\otimes_{i\in
S_{2},i\neq
d+k}\boldsymbol{U}_{i})(\mbox{\boldmath$\mathscr{M}$}_{t-1})_{(k)}^{\top}\right]\otimes\boldsymbol{I}_{p_{k}}\right\\}\text{vec}(\boldsymbol{U}_{d+k})$
(29) $\displaystyle=\left\\{\left[\boldsymbol{y}_{t-1}^{\top}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})\right]\otimes(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\right\\}\text{vec}(\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}),$
(30)

for $k=1,\dots,d$, where
$\boldsymbol{P}_{k,1}\in\mathbb{R}^{\prod_{i=1}^{d}r_{i}\times\prod_{i=1}^{d}r_{i}}$,
$\boldsymbol{P}_{k,2}\in\mathbb{R}^{p\times p}$ and
$\boldsymbol{P}_{k,3}\in\mathbb{R}^{r_{k}p_{k}\times r_{k}p_{k}}$ are
permutation matrices defined such that
$\boldsymbol{P}_{k,1}\text{vec}(\mbox{\boldmath$\mathscr{T}$}_{(k)})=\text{vec}(\mbox{\boldmath$\mathscr{T}$})$
for any $\mbox{\boldmath$\mathscr{T}$}\in\mathbb{R}^{r_{1}\times\cdots\times
r_{d}}$,
$\boldsymbol{P}_{k,2}\text{vec}(\mbox{\boldmath$\mathscr{T}$}_{(k)})=\text{vec}(\mbox{\boldmath$\mathscr{T}$})$
for any $\mbox{\boldmath$\mathscr{T}$}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$, and
$\boldsymbol{P}_{k,3}\text{vec}(\boldsymbol{U}_{k})=\text{vec}(\boldsymbol{U}_{k}^{\top})$,
and
$\mbox{\boldmath$\mathscr{M}$}_{t-1}\in\mathbb{R}^{r_{d+1}\times\cdots\times
r_{2d}}$ is defined such that
$\text{vec}(\mbox{\boldmath$\mathscr{M}$}_{t})=\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top}\boldsymbol{y}_{t-1}$. Therefore, we can
update each component parameter vector $\boldsymbol{\theta}_{i}$, and hence
the corresponding $\mathscr{G}$ and $\boldsymbol{U}_{i}$’s, iteratively by the
least squares method. The resulting ALS algorithm is shown in Algorithm 1.

Algorithm 1 ALS algorithm for LTR estimator

Initialize:
$\mbox{\boldmath$\mathscr{A}$}^{(0)}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{RRR}}$
or $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{OLS}}$
HOSVD:
$\mbox{\boldmath$\mathscr{A}$}^{(0)}\approx\mbox{\boldmath$\mathscr{G}$}^{(0)}\times_{i=1}^{2d}\boldsymbol{U}_{i}^{(0)}$,
with multilinear ranks $(r_{1},\dots,r_{2d})$.
repeat $s=0,1,2,\dots$
for $k=1,\dots,d$
$\boldsymbol{U}_{k}^{(s+1)}\leftarrow\underset{\boldsymbol{U}_{k}}{\arg\min}\sum_{t=1}^{T}\Big{\|}\boldsymbol{y}_{t}-(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i}^{(s)})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}\boldsymbol{P}_{k,1}\left\\{\left[(\otimes_{i\in
S_{1},i\neq
k}\boldsymbol{U}_{i}^{(s)})^{\top}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})_{(k)}^{\top}\right]\otimes\boldsymbol{I}_{r_{k}}\right\\}$
$\text{vec}(\boldsymbol{U}_{k})\Big{\|}_{2}^{2}$
end for
for $k=1,\dots,d$
$\boldsymbol{U}_{d+k}^{(s+1)}\leftarrow\underset{\boldsymbol{U}_{d+k}}{\arg\min}\sum_{t=1}^{T}\Big{\|}\boldsymbol{y}_{t}-\boldsymbol{P}_{k,2}\big{\\{}\left[(\otimes_{i\in
S_{2},i\neq
d+k}\boldsymbol{U}_{i}^{(s)})(\mbox{\boldmath$\mathscr{M}$}^{(s+1)}_{t-1})_{(k)}^{\top}\right]\otimes\boldsymbol{I}_{p_{k}}\big{\\}}\text{vec}(\boldsymbol{U}_{d+k})\Big{\|}_{2}^{2}$
end for
$\mbox{\boldmath$\mathscr{G}$}^{(s+1)}\leftarrow\underset{\scalebox{0.75}{\mbox{\boldmath$\mathscr{G}$}}}{\arg\min}\sum_{t=1}^{T}\Big{\|}\boldsymbol{y}_{t}-\left\\{\left[\boldsymbol{y}_{t-1}^{\top}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i}^{(s+1)})\right]\otimes(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i}^{(s+1)})\right\\}\text{vec}(\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]})\Big{\|}_{2}^{2}$
$\mbox{\boldmath$\mathscr{A}$}^{(s+1)}\leftarrow\mbox{\boldmath$\mathscr{G}$}^{(s+1)}\times_{i=1}^{2d}\boldsymbol{U}_{i}^{(s+1)}$
until convergence

As mentioned, the minimization in (24) is unconstrained. Accordingly, to
obtain $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$, no
orthogonality constraint of $\mathscr{G}$ and $\boldsymbol{U}_{i}$’s is needed
in Algorithm 1. Instead, we compute the final unique estimates of
$\mathscr{\widehat{G}}$ and $\boldsymbol{\widehat{U}}_{i}$’s by the HOSVD
operation based on the unconstrained solution
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}}$ obtained from Algorithm
1. Specifically, we compute each $\boldsymbol{\widehat{U}}_{i}$ uniquely as
the top $r_{i}$ left singular vectors of
$(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}})_{(i)}$ such that the
first element in each column of $\boldsymbol{\widehat{U}}_{i}$ is positive,
and set
$\mbox{\boldmath$\mathscr{\widehat{G}}$}=[\\![\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{LTR}};\boldsymbol{\widehat{U}}_{1}^{\top},\cdots,\boldsymbol{\widehat{U}}_{2d}^{\top}]\\!]$.
Similar alternating algorithms without imposing identification constraints can
be found in the literature of tensor decomposition; see, e.g. Zhou et al.,
(2013) and Li et al., (2018).

##  5 High-Dimensional Regularized Estimation

###  5.1 Regularization via One-Mode Matricization

The methods in the previous section rely on the assumptions that the dimension
is fixed and that the true transition tensor is exactly low-rank with known
Tucker ranks. We next relax both assumptions. Specifically, under the high-
dimensional setup, we consider regularized estimation of model (13) via
different nuclear-norm-type penalties and develop the corresponding non-
asymptotic theory under only an approximately low-Tucker-rank assumption on
the underlying true transition tensor.

In model (13), the exactly low-rank transition tensor
$\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{2d}}$
is subject to the constraints
$r_{i}=\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{(i)})$, for $i=1,\dots,2d$.
A commonly used convex relaxation of such multilinear rank constraints is the
regularization via the sum of nuclear (SN) norms of all the one-mode
matricizations,

$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SN}}=\sum_{i=1}^{2d}\|\mbox{\boldmath$\mathscr{A}$}_{(i)}\|_{*}.$
(31)

The SN norm has been widely used in the literature (Gandy et al.,, 2011;
Tomioka et al.,, 2011; Liu et al.,, 2013; Raskutti et al.,, 2019) to
simultaneously encourage the low-rankness for all modes of a tensor. This
leads us to the SN norm regularized estimator

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SN}}=\underset{\scalebox{0.7}{\mbox{\boldmath$\mathscr{A}$}}}{\operatorname*{arg\,min}}\left\\{\frac{1}{T}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}+\lambda_{\textup{SN}}\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SN}}\right\\},$
(32)

where $\lambda_{\text{SN}}$ is the tuning parameter. Note that if instead of
$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SN}}$, only one single nuclear norm,
say $\|\mbox{\boldmath$\mathscr{A}$}_{(1)}\|_{*}$, is penalized, then the
resulting estimator will only enforce the low-rankness for the first mode of
$\mathscr{A}$, while failing to do so for all the other $2d-1$ modes, and
hence will be less efficient than the above SN estimator.

To derive the estimation error bound for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SN}}$, we make the following
assumption on the random error
$\boldsymbol{e}_{t}=\textup{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})$.

######  Assumption 2.

Let
$\boldsymbol{e}_{t}=\boldsymbol{\Sigma}_{\boldsymbol{e}}^{1/2}\boldsymbol{\xi}_{t}$,
where $\\{\boldsymbol{\xi}_{t}\\}$ is a sequence of $i.i.d.$ random vectors,
with $\mathbb{E}(\boldsymbol{\xi}_{t})=\boldsymbol{0}$ and
$\textup{var}(\boldsymbol{\xi}_{t})=\boldsymbol{I}_{p}$, and
$\boldsymbol{\Sigma_{e}}=\textup{var}(\boldsymbol{e}_{t})$ is a positive
definite matrix. In addition, the entries $(\boldsymbol{\xi}_{it})_{1\leq
i\leq p}$ of $\boldsymbol{\xi}_{t}$ are mutually independent and
$\kappa^{2}$-sub-Gaussian, i.e., $\mathbb{E}(e^{\mu\xi_{it}})\leq
e^{\kappa^{2}\mu^{2}/2}$, for any $\mu\in\mathbb{R}$ and $i=1,\dots,p$.

The sub-Gaussianity condition in Assumption 2 is milder than the commonly used
normality assumption in the literature on high-dimensional stationary vector
autoregressive models (Basu and Michailidis,, 2015; Raskutti et al.,, 2019).
This relaxation is made possible through establishing a novel martingale-based
concentration bound in the proof of the deviation bound; see Lemma B.5 in
Appendix B.4. The covariance matrix $\boldsymbol{\Sigma_{e}}$ captures the
contemporaneous dependency in $\mbox{\boldmath$\mathscr{E}$}_{t}$, and the
constant $\kappa$ controls the tail heaviness of the marginal distributions.

For any $z\in\mathbb{C}$, let
$\mathcal{A}(z)=\boldsymbol{I}_{p}-\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}z$
be a matrix polynomial, where $\mathbb{C}$ is the set of complex numbers. Let
$\mu_{\min}(\mathcal{A})=\underset{|z|=1}{\min}\,\lambda_{\min}(\mathcal{A}^{*}(z)\mathcal{A}(z))$
and
$\mu_{\max}(\mathcal{A})=\underset{|z|=1}{\max}\,\lambda_{\max}(\mathcal{A}^{*}(z)\mathcal{A}(z))$,
where $\mathcal{A}^{*}(z)$ is the conjugate transpose of $\mathcal{A}(z)$; see
Basu and Michailidis, (2015) for more discussions on the connection between
the spectral density of the VAR process and the two quantities. In addition,
let

$\alpha_{\textup{RSC}}=\frac{\lambda_{\min}(\boldsymbol{\Sigma_{e}})}{\mu_{\max}(\mathcal{A})},\quad
M_{1}=\frac{\lambda_{\max}(\boldsymbol{\Sigma_{e}})}{\mu^{1/2}_{\min}(\mathcal{A})},\quad\text{and}\quad
M_{2}=\frac{\lambda_{\min}(\boldsymbol{\Sigma}_{\boldsymbol{e}})\mu_{\max}(\mathcal{A})}{\lambda_{\max}(\boldsymbol{\Sigma_{e}})\mu_{\min}(\mathcal{A})}.$

The exactly low-rankness of the true transition tensor $\mathscr{A}$ assumed
in Section 4 could be too stringent in real-world applications. In what
follows, we relax it to the following approximately low-rank assumption: We
assume that all one-mode matricizations of the underlying true transition
tensor $\mathscr{A}$ belong to the set of approximately low-rank matrices,
namely
$\mbox{\boldmath$\mathscr{A}$}_{(i)}\in\mathbb{B}_{q}(r_{q}^{(i)};p_{i},p_{-i}p)$
for some $q\in[0,1)$, where

$\mathbb{B}_{q}(r;d_{1},d_{2}):=\left\\{\boldsymbol{M}\in\mathbb{R}^{d_{1}\times
d_{2}}:\sum_{i=1}^{\min(d_{1},d_{2})}\sigma_{i}(\boldsymbol{M})^{q}\leq
r\right\\},$ (33)

and $p_{-i}=p/p_{i}=\prod_{j=1,j\neq i}^{d}p_{j}$. For the convenience of
notation, here we require that $0^{0}=0$. Note that when $q=0$,
$\mathbb{B}_{0}(r;d_{1},d_{2})$ is the set of $d_{1}$-by-$d_{2}$ rank-$r$
matrices. For $q>0$, the restriction on
$\sum_{i=1}^{\min(d_{1},d_{2})}\sigma_{i}(\boldsymbol{M})^{q}\leq r$ requires
that the singular values decay fast, and it is more general and natural than
the exactly low-rank assumption.

######  Theorem 2.

Under Assumptions 1 and 2, if $T\gtrsim\max_{1\leq i\leq
d}p_{-i}p+\max(\kappa^{2},\kappa^{4})M_{2}^{-2}p$,
$\lambda_{\textup{SN}}\gtrsim\kappa^{2}M_{1}d^{-2}\sum_{i=1}^{d}\sqrt{p_{-i}p/T}$,
and
$\mbox{\boldmath$\mathscr{A}$}_{(i)}\in\mathbb{B}_{q}(r_{q}^{(i)};p_{i},p_{-i}p)$
for some $q\in[0,1)$ and all $i=1,\dots,2d$, then with probability at least
$1-2\sum_{i=1}^{d}\exp(-Cp_{-i}p)-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$,

$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\textup{F}}\lesssim\sqrt{r_{q}}\left(\frac{2d\cdot\lambda_{\textup{SN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2}$
(34)

where $r_{q}=(2d)^{-1}\sum_{i=1}^{2d}r_{q}^{(i)}$.

By Theorem 2, when
$\lambda_{\textup{SN}}\asymp\kappa^{2}M_{1}d^{-2}\sum_{i=1}^{d}\sqrt{p_{-i}p/T}$,
the estimation error bound scales as $\sqrt{r_{q}}(\max_{1\leq i\leq
d}p_{-i}p/T)^{1/2-q/4}$; note that the factor $d$ in the error bounds is
canceled by the $d^{-2}$ in the rate of $\lambda_{\textup{SN}}$. When $q=0$,
namely $\mathscr{A}$ is exactly low-rank with Tucker ranks
$(r_{0}^{(1)},\dots,r_{0}^{(2d)})$, the error bound reduces to
$\sqrt{r_{0}\max_{1\leq i\leq d}p_{-i}p/T}$ and it is comparable to that in
Tomioka et al., (2011) for $i.i.d.$ tensor regression.

However, recent research (e.g., Mu et al.,, 2014; Raskutti et al.,, 2019) has
shown that the SN norm regularization approach is suboptimal, mainly because
$\mbox{\boldmath$\mathscr{A}$}_{(i)}$ is an unbalanced fat-and-short
matricization of a higher-order tensor. Technically, in the proof of Theorem
2, an essential intermediate step is to establish the deviation bound, where
we need to upper bound the operator norm of a sub-Gaussian random matrix with
the same dimensions as $\mbox{\boldmath$\mathscr{A}$}_{(i)}$; see Lemma B.5 in
Appendix B.3 for details. Undesirably, the order of this operator norm will be
dominated by the larger of the row and column dimensions of the matrix
$\mbox{\boldmath$\mathscr{A}$}_{(i)}\in\mathbb{R}^{p_{i}\times p_{-i}p}$, and
hence by the column dimension $p_{-i}p$, which eventually appears in the error
bound. This indicates that the imbalance of the matricization leads to the
efficiency bottleneck of the SN estimator.

On the other hand, similarly to (27), since the the reduced-rank VAR model can
be regarded as an overparameterization of the proposed LRTAR model,
alternatively one may focus on the approximately low-rank structure of the
transition matrix $\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$ in the VAR
representation in (15), and adopt the matrix nuclear (MN) estimator (Negahban
and Wainwright,, 2011) to estimate $\mathscr{A}$,

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{MN}}=\underset{\scalebox{0.7}{\mbox{\boldmath$\mathscr{A}$}}}{\operatorname*{arg\,min}}\left\\{\frac{1}{T}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}+\lambda_{\textup{MN}}\|\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}\|_{*}\right\\}.$
(35)

Note that the multi-mode matricization
$\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}=\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}^{\top}$
is a $p\times p$ square matrix. Thus, the loss of efficiency due to the
unbalanced matricization can be avoided, which is confirmed by the following
theorem.

######  Theorem 3.

Under Assumptions 1 and 2, if
$T\gtrsim[1+\max(\kappa^{2},\kappa^{4})M_{2}^{-2}]p$,
$\lambda_{\textup{MN}}\gtrsim\kappa^{2}M_{1}\sqrt{p/T}$, and
$\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}\in\mathbb{B}_{q}(s_{q}^{(1)};p,p)$
for some $q=[0,1)$, then with probability at least
$1-\exp(-Cp)-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$,

$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\textup{F}}\lesssim\sqrt{s_{q}^{(1)}}\left(\frac{\lambda_{\textup{MN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2}.$

Theorem 3 shows that, with
$\lambda_{\textup{MN}}\asymp\kappa^{2}M_{1}\sqrt{p/T}$, the estimation error
bound for $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{MN}}$ scales as
$\sqrt{s_{q}^{(1)}}(p/T)^{1/2-q/4}$. This result is comparable to that in
Negahban and Wainwright, (2011) for reduced-rank VAR models, yet we relax both
the constraint $\|\boldsymbol{A}\|_{\text{op}}<1$ on the transition matrix
$\boldsymbol{A}$ and the normality assumption on the random error in their
paper. This estimation error bound is clearly smaller than that in Theorem 2,
as $(\max_{1\leq i\leq d}p_{-i}p/T)^{1/2-q/4}$ in general can be much larger
than $(p/T)^{1/2-q/4}$ when $s_{q}^{(1)}\asymp r_{q}$. Therefore, adopting
square matricization can indeed improve the estimation performance.

The idea of using square matricization to improve efficiency was adopted by Mu
et al., (2014) in low-rank tensor completion problems. Their proposed method,
called the square deal, is to first unfold a general higher-order tensor into
a matrix with similar numbers of rows and columns, and then use the MN norm as
the regularizer. However, for our estimation problem, despite the advantage of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ over
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}}$, Theorem 3 reveals
another drawback of $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$.
That is, the error bounds for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ depend on the
$\ell_{q}$ radius $s_{q}^{(1)}$ of the singular values of
$\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$, suggesting that
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ may perform badly when
$s_{q}^{(1)}$ is relatively large. In other words, unless we have prior
knowledge that the $\ell_{q}$ “norm” of singular values of particular
matricization $\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$ is truly small,
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ may not be desirable
in practice.

On the other hand, although the SN regularizer in (31) suffers from
inefficiency due to the imbalance of one-mode matricizations, it has the
attractive feature of simultaneously enforcing the low-rankness across all
modes of the tensor $\mathscr{A}$, and thus is more efficient than its
counterpart which considers only one single one-mode matricization, say,
$\|\mbox{\boldmath$\mathscr{A}$}_{(1)}\|_{*}$. Similarly, if we can enforce
the approximate low-rankness across all possible square matricizations of
$\mathscr{A}$, the estimation performance may be further improved upon
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$. This motivates us to
propose a new regularization approach in the next subsection.

###  5.2 Regularization via Square Matricization

In this subsection, we propose a novel convex regularizer which improves upon
both the SN and MN regularizers in (31) and (35), respectively, by
simultaneously encouraging the low-rankness across all possible square
matricizations of the transition tensor $\mathscr{A}$.

For any $2d$-th-order tensor
$\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}\times p_{1}\times\cdots\times p_{d}}$, its multi-mode matricization
$\mbox{\boldmath$\mathscr{A}$}_{[I]}$ will be a $p\times p$ square matrix,
with $p=\prod_{i=1}^{d}p_{i}$, if the index set is chosen as

$I=\\{\ell_{1},\dots,\ell_{d}\\},$

where each index $\ell_{i}$ is set to either $i$ or $d+i$, for $i=1,\dots,d$.
For instance, $\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$ is the square
matricization formed by setting $\ell_{i}=i$ for all $i=1,\dots,d$. Moreover,
if $\mathscr{A}$ has multilinear ranks $(r_{1},\dots,r_{2d})$, then the rank
of the matricization $\mbox{\boldmath$\mathscr{A}$}_{[I]}$ is at most
$\min(\prod_{i=1,i\in I}^{2d}r_{i},\prod_{i=1,i\notin I}^{2d}r_{i})$.
Therefore, if we penalize the sum of nuclear norms of all such squares
matricizations, which we call the sum of square-matrix nuclear (SSN) norms for
simplicity, then the resulting estimator would enjoy the efficiency gain from
both the use of square matricizations and simultaneous incorporation of many
rank constraints.

Obviously, there are $2^{d}$ possible choices of the index set $I$ that
corresponds to a square matricization $\mbox{\boldmath$\mathscr{A}$}_{[I]}$.
However, since
$\mbox{\boldmath$\mathscr{A}$}_{[I]}=\mbox{\boldmath$\mathscr{A}$}_{[I^{\complement}]}^{\top}$,
when defining the SSN norm, we only need to include one of $I$ and its
complement $I^{\complement}$. A simple way to do so is to choose only sets
containing the index one. That is, fix $\ell_{1}=1$ and choose $\ell_{i}=i$ or
$d+i$ for $i=2,\dots,d$. This results in totally $2^{d-1}$ chosen index sets,
denoted by $I_{1},I_{2},\dots,I_{2^{d-1}}$. Note that
$I_{1}=S_{1}=\\{1,\dots,d\\}$. For example, when $d=3$, we have four chosen
index sets, $I_{1}=\\{1,2,3\\},I_{2}=\\{1,5,3\\},I_{3}=\\{1,2,6\\}$ and
$I_{4}=\\{1,5,6\\}$.

Based on the above choice of the $2^{d-1}$ index sets, we introduce the
following SSN norm,

$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}=\sum_{k=1}^{2^{d-1}}\left\|\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}\right\|_{*}.$
(36)

The corresponding estimator is defined as

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}=\underset{\scalebox{0.7}{\mbox{\boldmath$\mathscr{A}$}}}{\operatorname*{arg\,min}}\left\\{\frac{1}{T}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}+\lambda_{\textup{SSN}}\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}\right\\},$
(37)

where $\lambda_{\textup{SSN}}$ is the tuning parameter.

If the rank of one-mode matricizations
$\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{(i)})=r_{i}$, each square
matricization $\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}$ is also low-rank with
$\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]})\leq\min(\prod_{i=1,i\in
I_{k}}^{2d}r_{i},\prod_{i=1,i\notin I_{k}}^{2d}r_{i})$. Similarly, if all
$\mbox{\boldmath$\mathscr{A}$}_{(i)}$s are approximately low-rank, the square
matricizations are approximately low-rank as well. In contrast to the SN norm
in (31) which directly matches the multilinear ranks
$\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{(i)})$ for $i=1,\dots,d$, the SSN
norm encourages the multilinear low-rank structure of $\mathscr{A}$ by
simultaneously enforcing the low-rankness of all the square matricizations
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}$. The following theorem gives the
theoretical results for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}$.

######  Theorem 4.

Under Assumptions 1 and 2, if
$T\gtrsim[1+\max(\kappa^{2},\kappa^{4})M_{2}^{-2}]p$,
$\lambda_{\textup{SSN}}\gtrsim\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, and
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}\in\mathbb{B}(s_{q}^{(k)};p,p)$ for
some $q\in[0,1)$ and all $k=1,\dots,2^{d-1}$, with probability at least
$1-\exp[-C(p-d)]-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$,

$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\textup{F}}\lesssim\sqrt{s_{q}}\left(\frac{2^{d-1}\lambda_{\textup{SSN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2}$
(38)

where $s_{q}=2^{1-d}\sum_{k=1}^{2^{d-1}}s_{q}^{(k)}$.

| SN | MN | SSN
---|---|---|---
Sample size | $T\gtrsim(\max_{1\leq i\leq d}p_{-i}+M_{2}^{-2})p$ | $T\gtrsim(1+M_{2}^{-2})p$ | $T\gtrsim(1+M_{2}^{-2})p$
Estimation error | $\sqrt{r_{q}}(\max_{1\leq i\leq d}p_{-i}p/T)^{1/2-q/4}$ | $\sqrt{s_{q}^{(1)}}(p/T)^{1/2-q/4}$ | $\sqrt{s_{q}}(p/T)^{1/2-q/4}$
Table 1: Summary of the sample size conditions and error upper bounds in
Theorems 2–4, where $p_{-i}=\prod_{j=1,j\neq i}^{d}p_{j}$,
$r_{q}=(2d)^{-1}\sum_{i=1}^{2d}r_{q}^{(i)}$, and
$s_{q}=2^{1-d}\sum_{k=1}^{2^{d-1}}s_{q}^{(k)}$.

By Theorem 4, when
$\lambda_{\textup{SSN}}\asymp\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, the estimation
error bound scales as $\sqrt{s_{q}}(p/T)^{1/2-q/4}$ and reduces to
$\sqrt{s_{0}p/T}$ in the exactly low-rank setting for $q=0$. For a clearer
comparison among the three estimators
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}},\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$
and $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$, we summarize the
main results of Theorems 2–4 in Table 1. First, both
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ and
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ have much smaller
error bounds and less stringent sample size requirements than
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}}$, due to the diverging
dimension $p_{-i}$ in the results of the latter. This reaffirms the advantage
of the square matricizations.

Secondly, comparing $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$
to $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$, since the factor
$s_{q}$ in the error bounds of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ is the average of all
$s_{q}^{(k)}$ for $k=1,\dots,2^{d-1}$,
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ can protect us from
the bad scenarios where the $\ell_{q}$ “norm” of the singular values of
$\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$ is relatively large. If all the
$s_{q}^{(k)}$’s are of the same order, then the error upper bounds for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ and
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ in Table 1 will be
similar. However, our simulation results in Section 6 show that
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ clearly outperforms
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ under various
settings, even when $s_{q}^{(1)}=\cdots=s_{q}^{(2d)}$. Indeed, the error
bounds for $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ in Theorem
4 is likely to be loose, which is believed to be caused by taking the upper
bounds on the dual norm of the SSN norm in the proof of Lemma B.3; see
Appendix B.3 for details. By contrast, the error bounds for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ are minimax-optimal
(Han et al., 2015b, ). Therefore, although our theoretical results are not
sharp enough to distinguish clearly between the error rates of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ and
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$, we conjecture that
the actual rate of the former is generally smaller than that of the latter.
Methodologically, this is also easy to understand because, unlike
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$,
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ simultaneously
enforces the low-rankness across all square matricizations of $\mathscr{A}$
rather than just on $\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$.

######  Remark 1.

While our SSN regularization is proposed in the time series context, the idea
of imposing joint penalties on all (close to) square matricizations of the
coefficient tensor can potentially be extended to general higher-order tensor
estimation problems. Moreover, it can be refined to accommodate particular
structures of the data. For example, if some of the $d$ modes of the tensor-
value time series $\mbox{\boldmath$\mathscr{Y}$}_{t}$, namely
$p_{1},\dots,p_{d}$, are equal, then even a greater number of possible square
matricizations of $\mathscr{A}$ can be formed, resulting in improved
estimation efficiency.

###  5.3 Truncated Regularized Estimation

While the proposed regularized estimation methods do not require exact low-
rankness of the true transition tensor $\mathscr{A}$, sometimes imposing exact
low-rankness may be more desirable if one wants to interpret the underlying
low-dimensional tensor dynamics. As discussed in Section 3, the Tucker ranks
can be interpreted as the numbers of dynamic factors in each mode. In this
section, we propose a truncation method to consistently estimate the true
multilinear ranks $(r_{1},\dots,r_{2d})$ under the exact low-rank assumption.

Let $\gamma$ be a threshold parameter to be chosen properly. Given the
estimator $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}$, for each
$i=1,\dots,2d$, we calculate the singular value decomposition (SVD) of the
mode-$i$ matricization
$(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}})_{(i)}$ with the
singular values arranged in descending order. Next we truncate the SVD by
retaining only singular values greater than $\gamma$, and take their
corresponding left singular vectors to define the matrix
$\boldsymbol{\widetilde{U}}_{i}$. Then, the truncated core tensor is defined
as

$\mbox{\boldmath$\mathscr{\widetilde{G}}$}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}\times_{i=1}^{2d}\boldsymbol{\widetilde{U}}_{i}^{\top},$

based on which we propose the truncated sum of square-matrix nuclear (TSSN)
estimator

$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{TSSN}}=\mbox{\boldmath$\mathscr{\widetilde{G}}$}\times_{i=1}^{2d}\boldsymbol{\widetilde{U}}_{i}.$

To derive the theoretical results on rank selection, we make the following
assumption on the exact Tucker ranks and the magnitude of the singular values.

######  Assumption 3.

For all $i=1,\dots,2d$, $\sigma_{r}(\mbox{\boldmath$\mathscr{A}$}_{(i)})=0$
for all $r>r_{i}$, and there exists a constant $C>1$ such that $\min_{1\leq
i\leq 2d}\sigma_{r_{i}}\left(\mbox{\boldmath$\mathscr{A}$}_{(i)}\right)\geq
C\gamma$. As $T\rightarrow\infty$, the threshold parameter satisfies
$\gamma\gg(\kappa^{2}M_{1}/\alpha_{\textup{RSC}})\sqrt{s_{0}p/T}$, where
$s_{0}=2^{1-d}\sum_{k=1}^{2^{d-1}}\textup{rank}(\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]})$.

Assumption 3 requires that $\mathscr{A}$ has exact Tucker ranks
$(r_{1},\dots,r_{2d})$ which do not diverge too fast. The smallest positive
singular value for each $\mbox{\boldmath$\mathscr{A}$}_{(i)}$ is assumed to be
bounded away from the threshold $\gamma$ when the sample size is sufficiently
large. Since Assumption 3 involves unknown quantities, it cannot be used
directly for determining $\gamma$ in practice. Instead, we recommend using the
data-driven threshold parameter $\gamma$ described at the end of the next
subsection.

The rank selection consistency of the truncation method and the asymptotic
estimation error rate of
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{TSSN}}$ are given by the
following theorem.

######  Theorem 5.

Under the conditions of Theorem 4 and Assumption 3, if the tuning parameter
$\lambda_{\textup{SSN}}\asymp\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, then

$\mathbb{P}\left\\{\textup{rank}\left((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{TSSN}})_{(i)}\right)=\textup{rank}(\mbox{\boldmath$\mathscr{A}$}_{(i)}),~{}\text{for}~{}i=1,\dots,2d\right\\}\to
1,$

as $T\to\infty$, and for any fixed $d$,

$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{TSSN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\textup{F}}=O_{p}\left(\sqrt{s_{0}p/T}\right),$

where $s_{0}$ is defined as in Assumption 3.

###  5.4 ADMM Algorithm

This subsection presents the algorithm for the proposed (T)SSN regularized
estimator. The algorithm for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}}$ can be developed
analogously, while $\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}$ can
be obtained easily as in Negahban and Wainwright, (2011).

The objective function for the estimator
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}$ in (37) can be
rewritten as

$\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$})+\lambda_{\text{SSN}}\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}=\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$})+\lambda_{\text{SSN}}\sum_{k=1}^{2^{d-1}}\|\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}\|_{*},$
(39)

where
$\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$})=T^{-1}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}$
is the quadratic loss function. In (39), the regularizer
$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}$ involves $2^{d-1}$ nuclear
norms $\|\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}\|_{*}$, which are challenging
to handle at the same time. A similar difficulty also occurs in low-rank
tensor completion, for which Gandy et al., (2011) applied the alternating
direction method of multipliers (ADMM) algorithm (Boyd et al.,, 2011) to
efficiently separate the different nuclear norms. Borrowing the idea of Gandy
et al., (2011), we develop an ADMM algorithm for the miminization of (39).

Algorithm 2 ADMM algorithm for (T)SSN estimator

Initialize: $\mbox{\boldmath$\mathscr{C}$}_{k}^{(0)}$,
$\mbox{\boldmath$\mathscr{W}$}_{k}^{(0)}=\mbox{\boldmath$\mathscr{A}$}^{(0)}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{MN}}$,
for $k=1,\dots,2^{d-1}$, threshold parameter $\gamma$
for $j\in\\{0,1,\dots,J-1\\}$ do
$\mbox{\boldmath$\mathscr{A}$}^{(j+1)}\leftarrow\arg\min\Big{\\{}\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$})+\sum_{k=1}^{2^{d-1}}\rho\|\mbox{\boldmath$\mathscr{A}$}-\mbox{\boldmath$\mathscr{W}$}_{k}^{(j)}+\mbox{\boldmath$\mathscr{C}$}_{k}^{(j)}\|_{\text{F}}^{2}\Big{\\}}$
for $k\in\\{1,2,\dots,2^{d-1}\\}$ do
$\mbox{\boldmath$\mathscr{W}$}_{k}^{(j+1)}\leftarrow\arg\min\Big{\\{}\rho\|\mbox{\boldmath$\mathscr{A}$}^{(j+1)}-\mbox{\boldmath$\mathscr{W}$}_{k}+\mbox{\boldmath$\mathscr{C}$}_{k}^{(j)}\|_{\text{F}}^{2}+\lambda_{\text{SSN}}\|(\mbox{\boldmath$\mathscr{W}$}_{k})_{[I_{k}]}\|_{*}\Big{\\}}$
$\mbox{\boldmath$\mathscr{C}$}^{(j+1)}_{k}\leftarrow\mbox{\boldmath$\mathscr{C}$}_{k}^{(j)}+\mbox{\boldmath$\mathscr{A}$}^{(j+1)}-\mbox{\boldmath$\mathscr{W}$}^{(j+1)}_{k}$
end for
end for
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}\leftarrow\mbox{\boldmath$\mathscr{A}$}^{(J)}$
for $i\in\\{1,2,\dots,2d\\}$ do
$\boldsymbol{\widetilde{U}}_{i}\leftarrow\text{Truncated\\_SVD}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}})_{(i)},\gamma)$
end for
$\mbox{\boldmath$\mathscr{\widetilde{G}}$}\leftarrow\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}\times_{i=1}^{2d}\boldsymbol{\widetilde{U}}_{i}^{\top}$
$\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{TSSN}}\leftarrow\mbox{\boldmath$\mathscr{\widetilde{G}}$}\times_{i=1}^{2d}\boldsymbol{\widetilde{U}}_{i}$

To separate the $2^{d-1}$ nuclear norms in
$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}$, for each
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}$, we introduce a different dummy
variable $\mbox{\boldmath$\mathscr{W}$}_{k}$ as a surrogate for $\mathscr{A}$,
where $k=1,\dots,2^{d-1}$. Then the augmented Lagrangian is

$\mathcal{L}(\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{W}$},\mbox{\boldmath$\mathscr{C}$})=\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$})+\sum_{k=1}^{2^{d-1}}\Big{[}\lambda_{\text{SSN}}\|(\mbox{\boldmath$\mathscr{W}$}_{k})_{[I_{k}]}\|_{*}+2\rho\langle\mbox{\boldmath$\mathscr{C}$}_{k},\mbox{\boldmath$\mathscr{A}$}-\mbox{\boldmath$\mathscr{W}$}_{k}\rangle+\rho\|\mbox{\boldmath$\mathscr{A}$}-\mbox{\boldmath$\mathscr{W}$}_{k}\|_{\text{F}}^{2}\Big{]},$

where $\mbox{\boldmath$\mathscr{C}$}_{k}$ are the Lagrangian multipliers, for
$k=1,\dots,2^{d-1}$, and $\rho$ is the regularization parameter. Then we can
iteratively update
$\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{W}$}_{k}$ and
$\mbox{\boldmath$\mathscr{C}$}_{k}$ by the ADMM, as shown in Algorithm 2.

In Algorithm 2, the $\mathscr{A}$-update step is an $\ell_{2}$-regularized
least squares problem. Similarly to Gandy et al., (2011), the
$\mbox{\boldmath$\mathscr{W}$}_{k}$-update step can be solved by applying the
explicit soft-thresholding operator to the singular values of
$(\mbox{\boldmath$\mathscr{A}$}+\mbox{\boldmath$\mathscr{C}$}_{k})_{[I_{k}]}$.
Both subproblems have close-form solutions. Thus, the miminization of (39) can
be solved efficiently.

For the tuning parameter selection, since the cross-validation method is
unsuitable for time series or intrinsically ordered data, we apply the
Bayesian information criterion (BIC) to select the optimal
$\lambda_{\text{SSN}}$ from a sequence of tuning parameter values, where the
number of degrees of freedom is defined as
$2^{-(d-1)}\sum_{k=1}^{2^{d-1}}s_{k}(2p-s_{k})$. For the threshold parameter
$\gamma$ of the truncated estimator, we recommend
$\gamma=2^{d-1}\lambda_{\text{SSN}}/4$ to practitioners, where
$\lambda_{\text{SSN}}$ is the optimal tuning parameter selected by the BIC.

##  6 Simulation Studies

We conduct three simulation experiments to examine the finite-sample
performance of the proposed low- and high-dimensional estimation methods for
the LRTAR model in previous sections. Throughout, we generate the data from
model (13) with
$\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})\overset{i.i.d.}{\sim}N(\boldsymbol{0},\boldsymbol{I}_{p})$.
The entries of the core tensor $\mathscr{G}$ are generated independently from
$N(0,1)$ and rescaled such that
$\|\mbox{\boldmath$\mathscr{G}$}\|_{\text{F}}=5$. The factor matrices
$\boldsymbol{U}_{i}$’s are generated by extracting the leading singular
vectors from Gaussian random matrices while ensuring the stationarity
condition in Assumption 1.

The first experiment focuses on the proposed low-dimensional estimation method
considered in Section 4. We consider four cases of data generating processes.
In both cases (1a) and (1b), we consider $d=2$ and multilinear ranks
$(r_{1},r_{2},r_{3},r_{4})=(1,1,1,1)$, $(2,2,2,2)$, or $(2,3,2,3)$. In both
cases (2a) and (2b), we consider $d=3$ and multilinear ranks
$(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})=(1,1,1,1,1,1)$, $(2,2,2,1,1,1)$, or
$(2,2,2,2,2,2)$. Both pairs of cases differ in the setting of the dimensions
$p_{i}$’s: (1a) $p_{1}=p_{2}=5$; (1b) $p_{1}=p_{2}=10$; (2a)
$p_{1}=p_{2}=p_{3}=5$; and (2b) $p_{1}=p_{2}=p_{3}=7$. We repeat each data
generating process 300 times, and conduct the estimation using true
multilinear ranks. Figure 2 displays the average estimation error against
$T\in[1000,2000]$. The LTR estimator clearly outperforms both the RRR and OLS
estimators throughout all cases, confirming the theoretical efficiency
comparison in Corollary 1.

The second experiment aims to verify the estimation error bound for the
proposed SSN estimator in high dimensions. By Theorem 4, when the true
transition tensor $\mathscr{A}$ is exactly low-rank, we have
$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}^{2}=O_{p}(s_{0}p/T)$,
where $s_{0}=2^{1-d}\sum_{k=1}^{2^{d-1}}s_{0}^{(k)}$ with
$s_{0}^{(k)}=\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]})$. Note that
this rate is dependent on the overall dimension, $p=\prod_{i=1}^{d}p_{i}$, but
independent of the individual dimensions $p_{i}$. To examine the relationship
between the estimation error and the overal dimension $p$, sample size $T$,
multilinear ranks $r_{i}$ and individual dimensions $p_{i}$, we generate the
data under the eight settings listed in Table 2, and plot
$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}^{2}$
against the varying parameter in each case in Figure 3, including
$p,1/T,s_{0}$ and different settings of $p_{i}$’s under a fixed overall
dimension $p$. The first three columns of Figure 3 show that
$\|\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}-\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}^{2}$
roughly scales linearly in $p,T$ and $s_{0}$, while the last column suggests
that the estimation error is invariant under different settings of individual
$p_{i}$ as long as $p$ and the other parameters are fixed. This lends support
to the theoretical error bound for the SSN estimator.

In the third experiment, we compare the performance of all the four high-
dimensional estimators discussed in Section 5, namely the SN, MN, SSN and TSSN
estimators. We consider four cases of data generating processes. In cases (3a)
and (3b), we consider $d=2$ and multilinear ranks
$(r_{1},r_{2},r_{3},r_{4})=(1,1,1,1)$, $(2,2,1,1)$, or $(2,2,2,2)$. In cases
(4a) and (4b), we consider $d=3$ and multilinear ranks
$(r_{1},r_{2},r_{3},r_{4},r_{5},r_{6})=(1,1,1,1,1,1)$, $(2,2,2,1,1,1)$, or
$(2,2,2,2,2,2)$. Similarly to the first experiment, both pairs of cases differ
in the setting for $p_{i}$’s: (3a) $p_{1}=p_{2}=5$; (3b) $p_{1}=p_{2}=10$;
(4a) $p_{1}=p_{2}=p_{3}=5$; and (4b) $p_{1}=p_{2}=p_{3}=7$. For each setting,
we repeat 300 times, and conduct the estimation using SN, MN, SSN and TSSN.
The tuning parameter and the truncation threshold parameter are selected by
the method described in Section 5.4. In Figure 4, the average estimation error
is plotted against $T\in[400,1000]$ for cases (3a) and (3b), and
$T\in[600,1200]$ for cases (4a) and (4b). First, it can be seen that the SN
estimator is much inferior to the other three estimators, which is due to its
use of the unbalanced one-mode matricizations. Secondly, the SSN and TSSN
estimators outperform or are at least as good as the MN estimator in all
cases, and their advantage is remarkably clear even when
$r_{1}=\cdots=r_{2d}$. In addition, the TSSN estimator generally performs
better than the SSN, probably because the former yields a more parsimonious
model which further improves the estimation efficiency.

##  7 Real Data Analysis

###  7.1 Matrix-Valued Time Series

We first consider the modeling of a matrix-valued time series. The data is the
monthly market-adjusted return series of Fama–French $10\times 10$ portfolios
from January 1979 to December 2019, obtained from French, (2020). The
portfolios are constructed as the intersections of 10 portfolios formed by the
book-to-market (B/M) ratio and 10 portfolios formed by the size (market value
of equity). Hence, $d=2$, $p_{1}=p_{2}=10$ and $T=492$ months. We remove the
market effect for each portfolio by subtracting its average return from the
original return series.

Let $\boldsymbol{Y}_{t}\in\mathbb{R}^{10\times 10}$ be the observed matrix-
valued time series, with its rows and columns corresponding to different B/M
ratios (sorted from lowest to highest) and sizes (sorted from smallest to
largest), respectively, and denote
$\boldsymbol{y}_{t}=\text{vec}(\boldsymbol{Y}_{t})$. For comparison, we
consider the following five candidate models:

  * •

Vector autoregression (VAR):
$\boldsymbol{y}_{t}=\boldsymbol{A}\boldsymbol{y}_{t-1}+\boldsymbol{e}_{t}$,
where $\boldsymbol{A}\in\mathbb{R}^{100\times 100}$. The model is estimated by
the least squares method.

  * •

Vector factor model (VFM):
$\boldsymbol{y}_{t}=\boldsymbol{\Lambda}\boldsymbol{f}_{t}+\boldsymbol{e}_{t}$,
where $\boldsymbol{f}_{t}$ is the low-dimensional vector-valued latent factor,
and $\boldsymbol{\Lambda}$ is the loading matrix. The model is estimated by
the method in Lam et al., (2012), and for prediction, the estimated factors
$\boldsymbol{\widehat{f}}_{t}$ are then fitted by a VAR(1) model.

  * •

Matrix autoregression (MAR):
$\boldsymbol{Y}_{t}=\boldsymbol{B}_{1}\boldsymbol{Y}_{t-1}\boldsymbol{B}_{2}^{\top}+\boldsymbol{E}_{t}$,
where $\boldsymbol{B}_{1},\boldsymbol{B}_{2}\in\mathbb{R}^{10\times 10}$ are
coefficient matrices; see Example 1. The model is estimated by the iterated
least squares method in Chen et al., 2020b .

  * •

Matrix factor model (MFM):
$\boldsymbol{Y}_{t}=\boldsymbol{R}\boldsymbol{F}_{t}\boldsymbol{C}^{\top}+\boldsymbol{E}_{t}$,
where $\boldsymbol{F}_{t}$ is the low-dimensional matrix-valued latent factor,
and $\boldsymbol{R}$ and $\boldsymbol{C}$ are the loading matrices. The model
is estimated by the method in Wang et al., (2019), and for prediction, the
estimated factors $\boldsymbol{\widehat{F}}_{t}$ are then fitted by a VAR(1)
model.

  * •

The proposed LRTAR model:
$\boldsymbol{Y}_{t}=\langle\mbox{\boldmath$\mathscr{A}$},\boldsymbol{Y}_{t-1}\rangle+\boldsymbol{E}_{t}$,
with
$\mbox{\boldmath$\mathscr{A}$}=\mbox{\boldmath$\mathscr{G}$}\times_{i=1}^{4}\boldsymbol{U}_{i}$.
The model is estimated using the SSN and TSSN regularized estimation methods
in Section 5.

We first focus on the results for proposed LRTAR model. By the proposed
truncation method, we obtain the estimated multilinear ranks
$(\widehat{r}_{1},\widehat{r}_{2},\widehat{r}_{3},\widehat{r}_{4})=(8,8,2,2)$.
As shown in Example 1, the model can be written equivalently as
$\boldsymbol{U}_{3}^{\top}\boldsymbol{Y}_{t}\boldsymbol{U}_{4}=\left\langle\mbox{\boldmath$\mathscr{G}$},\boldsymbol{U}_{1}^{\top}\boldsymbol{Y}_{t-1}\boldsymbol{U}_{2}\right\rangle+\boldsymbol{U}_{3}^{\top}\boldsymbol{E}_{t}\boldsymbol{U}_{4}$.
Thus, $\boldsymbol{U}_{3}$ and $\boldsymbol{U}_{1}$ can be viewed as the
factor loadings across different B/M ratios, while $\boldsymbol{U}_{4}$ and
$\boldsymbol{U}_{2}$ represent those across different sizes. By the factor
interpretation below (17), this result indicates that the information in
$\boldsymbol{Y}_{t}$ can be effectively summarized into the $2\times 2$
response factors
$\boldsymbol{U}_{3}^{\top}\boldsymbol{Y}_{t}\boldsymbol{U}_{4}$, while the
low-rank structure associated with the predictor $\boldsymbol{Y}_{t-1}$ is not
strong. Moreover, the estimated multilinear ranks suggest that the MAR model
might be inefficient for the data, since the MAR imposes that
$r_{1}=r_{2}=r_{3}=r_{4}=10$ and
$\boldsymbol{U}_{3}=\boldsymbol{U}_{4}=\boldsymbol{I}_{10}$; see Example 1. In
addition, by an argument similar to the discussion in Example 2, the MFM can
be regarded as a special case of the proposed model with
$(\boldsymbol{U}_{1},\boldsymbol{U}_{2})=(\boldsymbol{U}_{3},\boldsymbol{U}_{4})$
and $(r_{1},r_{2})=(r_{3},r_{4})$. Thus, the fact that $\widehat{r}_{3}$ and
$\widehat{r}_{4}$ are much smaller than $\widehat{r}_{1}$ and
$\widehat{r}_{2}$ suggests that the MFM may be too restrictive for the data.

The TSSN estimates of the factor matrices are shown in Figure 5. The patterns
of the estimated response factor matrices
$\boldsymbol{\widetilde{U}}_{3},\boldsymbol{\widetilde{U}}_{4}\in\mathbb{R}^{10\times
2}$ are particularly interesting. First, for both
$\boldsymbol{\widehat{U}}_{3}$ and $\boldsymbol{\widehat{U}}_{4}$, the uniform
pattern of the first column indicate that portfolios across different B/M
ratios (or sizes) contribute approximately equally to the first B/M (or size)
response factor. Thus, this factor represents the component of the market
performance which is invariant to the size and B/M ratio of the portfolios.
The significance of this factor may be partially because we fit the model to
market-adjusted returns, where the average return is subtracted for all
portfolios. Meanwhile, for both response factor matrices, the second column
has a smoothly increasing pattern, suggesting that part of the return
variation in the market has a monotonic relationship with the size and B/M
ratio. Moreover, the above interpretations are consistent with the famous
Fama-French three-factor model, where the return of a portfolio is expected to
be affected by the market premium, outperformance of small versus big
companies, and outperformance of high B/M versus small B/M companies.

The performance of the five models are compared through both average in-sample
and out-of-sample forecasting errors. The average in-sample forecasting error
is calculated based on the fitted models for the entire data, while the
average out-of-sample forecasting error is calculated based on the rolling
forecast procedure as follows. From January 2016 ($t=445$) to December 2019
($t=492$), we fit the models using all the available data until time $t-1$ and
obtain the one-step-ahead forecast $\boldsymbol{\widehat{Y}}_{t}$. Then, we
obtain the average of the rolling forecasting errors for this period.

The average forecasting errors in $\ell_{2}$ and $\ell_{0}$ norms are
presented in Table 3. Firstly, the VAR model has the smallest in-sample
forecasting error among all models, which is as expected because the VAR model
is highly overparametrized. The bad in-sample performance of the MFM agrees
with our discussion about its restrictiveness for the data due to the mismatch
of the multilinear ranks. It is worth noting that the proposed LRTAR model has
competitive in-sample forecasting performance among all models.

The out-of-sample forecasting results provides a fuller picture of the
efficiency of different methods. It can be seen that the VAR and VFM models
perform worst among all, as they both completely ignore the matrix structure
of the data. Notably the proposed LRTAR model, fitted by either the SSN or the
TSSN methods, have the smallest out-of-sample forecasting errors. This
suggests that the proposed LRTAR model can indeed efficiently capture the
dynamic structural information in the data.

###  7.2 Three-Way Tensor-Valued Time Series

Fama and French, (2015) extended the classical Fama-French three-factor model
to a five-factor model which further incorporates the effect of operating
profitability (OP) and investment (Inv). This motivates tensor-valued stock
returns data formed according to the size (small and big), B/M ratio (four
groups sorted from lowest to highest), OP (four groups sorted from lowest to
highest), and Inv (four groups sorted from lowest to highest). We consider two
datasets retrieved from French, (2020). The first dataset consists of monthly
market-adjusted returns series of $4\times 4\times 2$ portfolios from July
1963 to December 2019, formed by the OP, B/M ratio and size. The second
dataset consists of those of $4\times 4\times 2$ portfolios formed by the Inv,
B/M ratio and size. Hence, both are tensor-valued time series with $d=3$,
$p_{1}=p_{2}=4$, $p_{3}=2$ and $T=678$ months.

Similar to the analysis in the previous subsection, five candidate models are
considered:

  * •

Vector autoregression (VAR):
$\boldsymbol{y}_{t}=\boldsymbol{A}\boldsymbol{y}_{t-1}+\boldsymbol{e}_{t}$,
where $\boldsymbol{A}\in\mathbb{R}^{32\times 32}$.

  * •

Vector factor model (VFM):
$\boldsymbol{y}_{t}=\boldsymbol{\Lambda}\boldsymbol{f}_{t}+\boldsymbol{e}_{t}$,
where $\boldsymbol{f}_{t}$ is the low-dimensional vector-valued latent factor,
and $\boldsymbol{\Lambda}$ is the loading matrix.

  * •

Multilinear tensor autoregression (MTAR):
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\mbox{\boldmath$\mathscr{Y}$}_{t-1}\times_{i=1}^{3}\boldsymbol{B}_{i}+\mbox{\boldmath$\mathscr{E}$}_{t}$,
where $\boldsymbol{B}_{1},\boldsymbol{B}_{2}\in\mathbb{R}^{4\times 4}$ and
$\boldsymbol{B}_{3}\in\mathbb{R}^{2\times 2}$ are coefficient matrices; see
Example 1.

  * •

Tensor factor model (TFM):
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\mbox{\boldmath$\mathscr{F}$}_{t}\times_{i=1}^{3}\boldsymbol{U}_{i}+\mbox{\boldmath$\mathscr{E}$}_{t}$,
where $\mbox{\boldmath$\mathscr{F}$}_{t}$ is the low-dimensional tensor-valued
latent factor, and $\boldsymbol{U}_{i}$’s are the loading matrices.

  * •

The proposed LRTAR model:
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle+\mbox{\boldmath$\mathscr{E}$}_{t}$,
with
$\mbox{\boldmath$\mathscr{A}$}=\mbox{\boldmath$\mathscr{G}$}\times_{i=1}^{6}\boldsymbol{U}_{i}$.

The TFM is estimated by the method in Chen et al., 2020c , and for prediction,
the estimated factors $\mbox{\boldmath$\mathscr{\widehat{F}}$}_{t}$ are fitted
by a VAR(1) model. The other four models are fitted by the same methods as in
the previous subsection. The multilinear ranks selected by the truncation
method are
$(\widehat{r}_{1},\widehat{r}_{2},\widehat{r}_{3},\widehat{r}_{4},\widehat{r}_{5},\widehat{r}_{6})=(3,3,2,1,1,1)$
and $(2,2,2,1,1,1)$ for the OP-BM-Size portfolio return series and the Inv-BM-
Size portfolio return series, respectively. Note that similar to the BM-Size
$10\times 10$ series in the previous subsection, the low-rank structure for
the response in these two tensor-valued datasets is more evident than that for
the predictor. Figure 6 shows the TSSN estimates of the factor matrices. We
find that the estimated response factors have a uniform pattern similar to
that in Figure 5 for the matrix-valued data. The fact that only one response
factor is extracted in each direction suggests that there might not be
substantial effect of OP, Inv, B/M ratio or size on the returns for these
datasets.

Finally, using the same methods as in the previous subsection, we calculate
the average in-sample and out-of-sample forecasting errors for both datasets
fitted with the five models. As shown in Table 4, the comparison results for
the two datasets are quite similar. It can be clearly observed that the VAR
model always has the smallest in-sample forecasting error, yet the largest
out-of-sample forecasting error. On the contrary, the proposed LRTAR model,
fitted by either the SSN or the TSSN methods, has the smallest out-of-sample
forecasting error. Moreover, the in-sample forecasting performance of the
LRTAR model is competitive even compared to the VAR model. Similarly to the
MFM in the previous subsection, the TFM model has poor in-sample performance,
possibly due to the discrepancy between $(r_{1},r_{2},r_{3})$ and
$(r_{4},r_{5},r_{6})$, as reflected by the estimated multilinear ranks. In
sum, the results support the efficiency and flexibility of the proposed model
and estimation methods for tensor-valued time series data.

##  8 Conclusion and Discussion

Efficient modeling and forecasting of general structured (tensor-valued),
high-dimensional time series data is an important research topic which however
has been rarely explored in the literature so far. This paper makes the first
thorough attempt to address this problem by introducing the low-rank tensor
autoregressive model. By assuming the exactly or approximately low-Tucker-rank
structure of the transition tensor, the model exploits the low-dimensional
tensor dynamic structure of the high-dimensional time series data, and
summarizes the complex temporal dependencies into interpretable dynamic
factors.

Asymptotic and non-asymptotic properties are derived for the proposed low- and
high-dimensional estimators, respectively. For the latter, we relax the
conventional Gaussian assumption in the high-dimensional time series
literature to sub-Gaussianity via a new martingale-based concentration
technique. Moreover, based on the special structure of the transition tensor,
a novel convex regularizer, the SSN, is proposed, gaining efficiencies from
both the square matricization and simultaneous penalization across modes. A
truncation method, the TSSN, is further introduced to consistently select the
multilinear ranks and enhance model interpretability.

We discuss several directions for future research. First, the proposed
estimators cannot adapt to the contemporaneous correlation among elements of
the random error $\mbox{\boldmath$\mathscr{E}$}_{t}$, which may lead to
efficiency loss. This issue can be addressed by considering the generalized
least squares loss function,
$\mathcal{L}_{T}(\mbox{\boldmath$\mathscr{A}$};\boldsymbol{\Sigma}_{\boldsymbol{e}})=T^{-1}\sum_{t=1}^{T}\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle)^{\top}\boldsymbol{\Sigma}^{-1}_{\boldsymbol{e}}\text{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t}-\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle)$,
where
$\boldsymbol{\Sigma}_{\boldsymbol{e}}=\text{var}(\text{vec}(\mbox{\boldmath$\mathscr{E}$}_{t}))$.
Then $\mathscr{A}$ may be estimated jointly with
$\boldsymbol{\Sigma}_{\boldsymbol{e}}$ or by a two-step approach based on a
consistent estimator $\boldsymbol{\widehat{\Sigma}}_{\boldsymbol{e}}$ (Basu
and Michailidis,, 2015; Davis et al.,, 2016).

Secondly, the proposed methods can be generalized to the LRTAR model of finite
lag order $L$, defined as
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\langle\mbox{\boldmath$\mathscr{A}$}_{1},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle+\cdots+\langle\mbox{\boldmath$\mathscr{A}$}_{L},\mbox{\boldmath$\mathscr{Y}$}_{t-L}\rangle+\mbox{\boldmath$\mathscr{E}$}_{t}$,
where
$\mbox{\boldmath$\mathscr{A}$}_{1},\dots,\mbox{\boldmath$\mathscr{A}$}_{L}$
are $2d$-th-order multilinear low-rank coefficient tensors. Then, one may
consider the SSN regularized estimation by minimizing
$T^{-1}\sum_{t=1}^{T}\|\mbox{\boldmath$\mathscr{Y}$}_{t}-\sum_{j=1}^{L}\langle\mbox{\boldmath$\mathscr{A}$}_{j},\mbox{\boldmath$\mathscr{Y}$}_{t-j}\rangle\|_{\text{F}}^{2}+\sum_{j=1}^{L}\lambda_{j}\|\mbox{\boldmath$\mathscr{A}$}_{j}\|_{\text{SSN}}$.
Alternatively, the $\mbox{\boldmath$\mathscr{A}$}_{j}$’s can be combined into
a $(2d+1)$-th-order tensor
$\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots p_{d}\times
p_{1}\times\cdots\times p_{d}\times L}$ whose mode-$(2d+1)$ matricization is
$\mbox{\boldmath$\mathscr{A}$}_{(2d+1)}=(\mbox{\boldmath$\mathscr{A}$}_{1},\dots,\mbox{\boldmath$\mathscr{A}$}_{L})$;
see Wang et al., (2020) for a similar idea. Even though $\mathscr{A}$ may not
have exactly square matricizations for $L>1$, the proposed SSN can still be
adapted by employing approximately square matricizations. For instance,
consider the $pL\times p$ multi-mode matricizations
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}\cup\\{2d+1\\}]}$, where the index sets
$I_{k}$ for $k=1,\dots,2^{d-1}$ are defined as in this paper. Then, a
generalized SSN norm can be constructed as
$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{GSSN}}=\sum_{k=1}^{2^{d-1}}\|\mbox{\boldmath$\mathscr{A}$}_{[I_{k}\cup\\{2d+1\\}]}\|_{*}$,
which will reduce to $\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{SSN}}$ when
$L=1$.

Thirdly, the LRTAR model can be readily extended to incorporate exogenous
tensor-valued predictors, giving rise to the LRTAR-X model,
$\mbox{\boldmath$\mathscr{Y}$}_{t}=\langle\mbox{\boldmath$\mathscr{A}$},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle+\langle\mbox{\boldmath$\mathscr{B}$},\mbox{\boldmath$\mathscr{X}$}_{t}\rangle+\mbox{\boldmath$\mathscr{E}$}_{t}$,
where $\mbox{\boldmath$\mathscr{X}$}_{t}$ is an $m$-th-order tensor of
exogenous variables, and $\mathscr{B}$ is a $(d+m)$-th-order coefficient
tensor. When the dimension of $\mbox{\boldmath$\mathscr{X}$}_{t}$ is high, a
low-dimensional structure, such as sparsity, group sparsity or low-rankness,
can be imposed on $\mathscr{B}$ to improve the estimation efficiency.

Moreover, an open question for matrix- or tensor-valued time series models is
how to incorporate additional structures or constraints; for instance,
transport or trade flow data have unspecified diagonal entries (Chen and
Chen,, 2019), and realized covariance matrix data are subject to positive
definite constraints. Beyond the time series context, it is also worth
investigating the generalization of the SSN regularization method in higher-
order tensor estimation and completion applications, such as neuroimaging
analysis (Li et al.,, 2018), recommender sytem (Bi et al.,, 2018) and natural
language processing (Frandsen and Ge,, 2019).

##  Appendix A: Proofs for Low-Dimensional Estimation

Below we give the proofs of Theorem 1 and Corollary 1 in Section 4, which
generally follow from Proposition 4.1 in Shapiro, (1986) for overparameterized
models.

###### Proof of Theorem 1.

The proposed model in (15) can be written in the matrix form

$\underbrace{\begin{bmatrix}\boldsymbol{y}_{1}^{\top}\\\
\boldsymbol{y}_{2}^{\top}\\\ \vdots\\\
\boldsymbol{y}_{T}^{\top}\end{bmatrix}}_{\boldsymbol{Y}}=\underbrace{\begin{bmatrix}\boldsymbol{y}_{0}^{\top}\\\
\boldsymbol{y}_{1}^{\top}\\\ \vdots\\\
\boldsymbol{y}_{T-1}^{\top}\end{bmatrix}}_{\boldsymbol{X}}\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}^{\top}+\underbrace{\begin{bmatrix}\boldsymbol{e}_{1}^{\top}\\\
\boldsymbol{e}_{2}^{\top}\\\ \vdots\\\
\boldsymbol{e}_{T}^{\top}\end{bmatrix}}_{\boldsymbol{E}}.$ (A.1)

Let
$\boldsymbol{\widehat{h}}_{\text{OLS}}=\text{vec}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{OLS}})_{[S_{2}]})$.
Under Assumption 1, by the classical asymptotic theory for stationary VAR
models, as $T\to\infty$, we have

$\sqrt{T}(\boldsymbol{\widehat{h}}_{\text{OLS}}-\boldsymbol{h})\overset{d}{\to}N(\boldsymbol{0},\boldsymbol{\Sigma}_{\textup{OLS}}),$
(A.2)

where
$\boldsymbol{\Sigma}_{\textup{OLS}}=\boldsymbol{\Sigma}_{\boldsymbol{e}}\otimes\boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1}$,
with $\boldsymbol{\Sigma}_{\boldsymbol{e}}=\text{var}(\boldsymbol{e}_{t})$ and
$\boldsymbol{\Sigma}_{\boldsymbol{y}}=\text{var}(\boldsymbol{y}_{t})$.

Following Shapiro, (1986), consider the discrepancy function

$F(\boldsymbol{\widehat{h}}_{\text{OLS}},\boldsymbol{h})=\|\text{vec}(\boldsymbol{Y})-(\boldsymbol{I}_{p}\otimes\boldsymbol{X})\boldsymbol{h}\|_{2}^{2}-\|\text{vec}(\boldsymbol{Y})-(\boldsymbol{I}_{p}\otimes\boldsymbol{X})\boldsymbol{\widehat{h}}_{\text{OLS}}\|_{2}^{2}.$
(A.3)

Note that $F(\boldsymbol{\widehat{h}}_{\text{OLS}},\boldsymbol{h})$ is
nonnegative and twice continuously differentiable, and equals zero if and only
if $\boldsymbol{\widehat{h}}_{\text{OLS}}=\boldsymbol{h}$.

The Jacobian matrix
$\boldsymbol{H}={\partial\boldsymbol{h}(\boldsymbol{\theta})}/{\partial\boldsymbol{\theta}}$
in (25) can be verified by noting that

$\boldsymbol{h}=\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})=\text{vec}((\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}(\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})^{\top})=\left((\otimes_{i\in
S_{1}}\boldsymbol{U}_{i})\otimes(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\right)\text{vec}(\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]})$

and that for any $1\leq i\leq 2d$,

$\begin{split}\boldsymbol{h}=\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})=\boldsymbol{P}^{(i)}_{[S_{2}]}\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{(i)})&=\boldsymbol{P}_{[S_{2}]}^{(i)}\text{vec}\left(\boldsymbol{U}_{i}\mbox{\boldmath$\mathscr{G}$}_{(i)}(\otimes_{j=1,j\neq
i}^{2d}\boldsymbol{U}_{j}^{\top})\right)\\\
&=\boldsymbol{P}^{(i)}_{[S_{2}]}\left\\{\left[(\otimes_{j=1,j\neq
i}^{2d}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{(i)}^{\top}\right]\otimes\boldsymbol{I}_{p_{i}}\right\\}\text{vec}(\boldsymbol{U}_{i}).\end{split}$
(A.4)

Then, by Proposition 4.1 in Shapiro, (1986), we obtain that the minimizer of
$F(\boldsymbol{\widehat{h}}_{\text{OLS}},\cdot)$, namely the LTR estimator,
has the asymptotic normality,

$\sqrt{T}(\boldsymbol{\widehat{h}}_{\text{LTR}}-\boldsymbol{h})\overset{d}{\to}N(\boldsymbol{0},\boldsymbol{\Sigma}_{\text{LTR}})$
(A.5)

and
$\boldsymbol{\Sigma}_{\text{LTR}}=\boldsymbol{P}\boldsymbol{\Sigma}_{\textup{OLS}}\boldsymbol{P}^{\top}$,
where
$\boldsymbol{P}=\boldsymbol{H}(\boldsymbol{H}^{\top}\boldsymbol{JH})^{\dagger}\boldsymbol{H}^{\top}\boldsymbol{J}$
is the projection matrix,
$\boldsymbol{J}=\boldsymbol{\Sigma}_{\boldsymbol{e}}^{-1}\otimes\boldsymbol{\Sigma}_{\boldsymbol{y}}$
is the Fisher information matrix of $\boldsymbol{h}$, and $\dagger$ denotes
the Moore-Penrose inverse. Since
$\boldsymbol{\Sigma}_{\textup{OLS}}=\boldsymbol{J}^{-1}$, we can obtain that
$\boldsymbol{\Sigma}_{\text{LTR}}=\boldsymbol{H}(\boldsymbol{H}^{\top}\boldsymbol{JH})^{\dagger}\boldsymbol{H}^{\top}$.
∎

###### Proof of Corollary 1.

As discussed in the proof of Theorem 1,
$\boldsymbol{\Sigma}_{\text{OLS}}=\boldsymbol{J}^{-1}$, and observe that

$\boldsymbol{\Sigma}_{\text{LTR}}=\boldsymbol{H}(\boldsymbol{H}^{\top}\boldsymbol{JH})^{\dagger}\boldsymbol{H}^{\top}=\boldsymbol{J}^{-1/2}\boldsymbol{Q}_{\boldsymbol{J}^{1/2}\boldsymbol{H}}\boldsymbol{J}^{-1/2},$
(A.6)

where $\boldsymbol{Q}_{\boldsymbol{J}^{1/2}\boldsymbol{H}}$ is the projection
matrix onto the orthogonal compliment of
$\text{span}(\boldsymbol{J}^{1/2}\boldsymbol{H})$.

On the other hand, under the proposed model, the transition matrix can be
decomposed as
$\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}=\boldsymbol{V}_{1}\boldsymbol{V}_{2}^{\top}$,
where $\boldsymbol{V}_{1}=\otimes_{i\in S_{1}}\boldsymbol{U}_{i}$,
$\boldsymbol{V}_{2}=(\otimes_{i\in
S_{2}}\boldsymbol{U}_{i})\mbox{\boldmath$\mathscr{G}$}_{[S_{2}]}$, and
$\text{rank}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})\leq
s_{1}=\min(\prod_{i=1}^{d}r_{i},\prod_{i=d+1}^{2d}r_{i})$. Thus, we can write
$\boldsymbol{h}=\text{vec}(\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]})=\boldsymbol{h}(\boldsymbol{\phi})$,
where
$\boldsymbol{\phi}=(\text{vec}(\boldsymbol{V}_{1})^{\top},\text{vec}(\boldsymbol{V}_{2})^{\top})^{\top}$
is the parameter vector for the RRR estimator. Then, similarly to the proof of
Theorem 1, we denote the Jacobian matrix by
$\boldsymbol{R}=\partial\boldsymbol{h}/\partial\boldsymbol{\phi}$. By similar
arguments, we have
$\sqrt{T}(\boldsymbol{\widehat{h}}_{\text{RRR}}-\boldsymbol{h})\overset{d}{\to}N(\boldsymbol{0},\boldsymbol{\Sigma}_{\text{RRR}})$,
where
$\boldsymbol{\widehat{h}}_{\text{RRR}}=\text{vec}((\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{RRR}})_{[S_{2}]})$,
and

$\boldsymbol{\Sigma}_{\text{RRR}}=\boldsymbol{R}(\boldsymbol{R}^{\top}\boldsymbol{JR})^{\dagger}\boldsymbol{R}^{\top}=\boldsymbol{J}^{-1/2}\boldsymbol{Q}_{\boldsymbol{J}^{1/2}\boldsymbol{R}}\boldsymbol{J}^{-1/2},$
(A.7)

where $\boldsymbol{Q}_{\boldsymbol{J}^{1/2}\boldsymbol{R}}$ is the projection
matrix onto the orthogonal compliment of
$\text{span}(\boldsymbol{J}^{1/2}\boldsymbol{R})$. Hence, it is clear that
$\boldsymbol{\Sigma}_{\text{RRR}}\leq\boldsymbol{J}^{-1}=\boldsymbol{\Sigma}_{\text{OLS}}$.

Moreover, since the Tucker decomposition can be viewed as a further
decomposition of the low-rank decomposition $\boldsymbol{VU}^{\top}$ for
$\mbox{\boldmath$\mathscr{A}$}_{[S_{2}]}$, we have
$\boldsymbol{H}=\partial\boldsymbol{h}/\partial\boldsymbol{\theta}=\boldsymbol{R}\cdot\partial\boldsymbol{\phi}/\partial\boldsymbol{\theta}$.
By (A.6) and (A.7), we have
$\boldsymbol{\Sigma}_{\text{LTR}}\leq\boldsymbol{\Sigma}_{\text{RRR}}$, since
$\text{span}(\boldsymbol{J}^{1/2}\boldsymbol{H})\subset\text{span}(\boldsymbol{J}^{1/2}\boldsymbol{R})$.
∎

##  Appendix B: Proofs for High-Dimensional Estimation

In this appendix, we provide the proofs of Theorems 2–5 in Section 5. We start
with a preliminary analysis in Appendix B.1 which lays out the common
technical framework for proving the estimation and prediction error bounds of
the SN, MN and SSN regularized estimators, and four lemmas, Lemmas B.1–B.4,
are introduced herein. Then in Appendix B.2 we give the proofs of Theorems
2–5. The proofs of Lemmas B.1–B.4 are provided in Appendix B.3, and three
auxiliary lemmas are collected in Appendix B.4

###  B.1 Preliminary Analysis

The technical framework for proving the error bounds in Theorem 2–4 consists
of two main steps, a deterministic analysis and a stochastic analysis, given
in Sections B.1.1 and B.1.2, respectively. The goal of the first one is to
derive the error bounds given the deterministic realization of the time
series, assuming that the parameters satisfy certain regularity conditions.
The goal of the second one is to verify that under stochasticity these
regularity conditions are satisfied with high probability.

####  B.1.1 Deterministic Analysis

Throughout the appendix, we adopt the following notations. We use $C$ to
denote a generic positive constant, which is independent of the dimensions and
the sample size. For any matrix $\boldsymbol{M}$ and a compatible subspace
$\mathcal{S}$, we denote by $\boldsymbol{M}_{\mathcal{S}}$ the projection of
$\boldsymbol{M}$ onto $\mathcal{S}$. In addition, let
$\text{col}(\boldsymbol{M})$ be the column space of $\boldsymbol{M}$, and let
$\mathcal{S}^{\perp}$ be the complement of the subspace $\mathcal{S}$. For a
generic tensor
$\mbox{\boldmath$\mathscr{W}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{2d}}$,
the dual norms of its SSN norm and SN norm, denoted by
$\|\mbox{\boldmath$\mathscr{W}$}\|_{\text{SSN}^{*}}$ and
$\|\mbox{\boldmath$\mathscr{W}$}\|_{\text{SN}^{*}}$, respectively, are defined
as

$\begin{split}\|\mbox{\boldmath$\mathscr{W}$}\|_{\text{SSN}^{*}}=\sup_{\scalebox{0.7}{\mbox{\boldmath$\mathscr{T}$}}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{2d}},\|\scalebox{0.7}{\mbox{\boldmath$\mathscr{T}$}}\|_{\text{SSN}}\leq
1}\langle\mbox{\boldmath$\mathscr{W}$},\mbox{\boldmath$\mathscr{T}$}\rangle,~{}~{}\text{and}~{}~{}\|\mbox{\boldmath$\mathscr{W}$}\|_{\text{SN}^{*}}=\sup_{\scalebox{0.7}{\mbox{\boldmath$\mathscr{T}$}}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{2d}},\|\scalebox{0.7}{\mbox{\boldmath$\mathscr{T}$}}\|_{\text{SN}}\leq
1}\langle\mbox{\boldmath$\mathscr{W}$},\mbox{\boldmath$\mathscr{T}$}\rangle.\end{split}$
(B.1)

Moreover, for any two tensors
$\mbox{\boldmath$\mathscr{X}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{m}}$
and $\mbox{\boldmath$\mathscr{Y}$}\in\mathbb{R}^{p_{m+1}\times\cdots\times
p_{m_{n}}}$, their tensor outer product is defined as
$(\mbox{\boldmath$\mathscr{X}$}\circ\mbox{\boldmath$\mathscr{Y}$})\in\mathbb{R}^{p_{1}\times\cdots\times
p_{m}\times p_{m+1}\times\cdots\times p_{m+n}}$ where

$(\mbox{\boldmath$\mathscr{X}$}\circ\mbox{\boldmath$\mathscr{Y}$})_{i_{1}\dots
i_{m}i_{m+1}\dots i_{m+n}}=\mbox{\boldmath$\mathscr{X}$}_{i_{1}\dots
i_{m}}\mbox{\boldmath$\mathscr{Y}$}_{i_{m+1}\dots i_{m+n}},$ (B.2)

for any $1\leq i_{1}\leq p_{1}$, $\dots$, $1\leq i_{m+n}\leq p_{m+n}$.

For the theory of regularized $M$-estimators, restricted error sets and
restricted strong convexity are essential definitions. To define the former,
we need to first introduce the following restricted model subspaces.

For $i=1,\dots,2d$, denote by $\widetilde{\mathcal{U}}_{i}$ and
$\widetilde{\mathcal{V}}_{i}$ the spaces spanned by the first $r_{i}$ left and
right singular vectors in the SVD of $\mbox{\boldmath$\mathscr{A}$}_{(i)}$,
respectively. Define the collections of subspaces

$\mathcal{N}=(\mathcal{N}_{1},\dots,\mathcal{N}_{2d})\quad\text{and}\quad\overline{\mathcal{N}}^{\perp}=(\overline{\mathcal{N}}_{1}^{\perp},\dots,\overline{\mathcal{N}}_{2d}^{\perp}),$

where

$\begin{split}\mathcal{N}_{i}&=\\{\boldsymbol{M}\in\mathbb{R}^{p_{i}\times
p_{-i}p}|\text{col}(\boldsymbol{M})\subset\widetilde{\mathcal{U}}_{i},\text{col}(\boldsymbol{M}^{\top})\subset\widetilde{\mathcal{V}}_{i}\\},\\\
\overline{\mathcal{N}}_{i}^{\perp}&=\\{\boldsymbol{M}\in\mathbb{R}^{p_{i}\times
p_{-i}p}|\text{col}(\boldsymbol{M})\perp\widetilde{\mathcal{U}}_{i},\text{col}(\boldsymbol{M}^{\top})\perp\widetilde{\mathcal{V}}_{i}\\},\end{split}$
(B.3)

for $i=1,\dots,2d$. Note that
$\mathcal{N}_{i}\subset\overline{\mathcal{N}}_{i}$.

Furthermore, for $k=1,\dots,2^{d-1}$, denote by $\mathcal{U}_{k}$ and
$\mathcal{V}_{k}$ the spaces spanned by the first $s_{k}^{*}$ left and right
singular vectors in the SVD of the square matricization
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}$, respectively, where
$s_{k}^{*}=\textup{rank}(\mbox{\boldmath$\mathscr{A}$})_{[I_{k}]}$. Similarly,
define the collections of subspaces

$\mathcal{M}:=(\mathcal{M}_{1},\dots,\mathcal{M}_{2^{d-1}})\quad\text{and}\quad\overline{\mathcal{M}}^{\perp}=(\overline{\mathcal{M}}_{1}^{\perp},\dots,\overline{\mathcal{M}}_{2^{d-1}}^{\perp}),$

where

$\begin{split}&\mathcal{M}_{k}=\\{\boldsymbol{M}\in\mathbb{R}^{p\times
p}|\text{col}(\boldsymbol{M})\subset\mathcal{U}_{k},~{}\text{col}(\boldsymbol{M}^{\top})\subset\mathcal{V}_{k}\\},\\\
&\overline{\mathcal{M}}_{k}^{\perp}=\\{\boldsymbol{M}\in\mathbb{R}^{p\times
p}|\text{col}(\boldsymbol{M})\perp\mathcal{U}_{k},~{}\text{col}(\boldsymbol{M}^{\top})\perp\mathcal{V}_{k}\\},\end{split}$
(B.4)

for $k=1,\dots,2^{d-1}$. In particular, as described in Section 5.2,
$I_{1}=S_{1}=\\{1,\dots,d\\}$. Thus, $\mathcal{M}_{1}$ and
$\overline{\mathcal{M}}_{1}^{\perp}$ are the subspaces associated with the
square matricization $\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}$.

Then, for simplicity, for any
$\mbox{\boldmath$\mathscr{W}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{2d}}$,
we denote

$\begin{split}&\mbox{\boldmath$\mathscr{W}$}_{\mathcal{N}}^{(i)}=(\mbox{\boldmath$\mathscr{W}$}_{(i)})_{\mathcal{N}_{i}},\quad\mbox{\boldmath$\mathscr{W}$}_{\mathcal{N}^{\perp}}^{(i)}=(\mbox{\boldmath$\mathscr{W}$}_{(i)})_{\mathcal{N}_{i}^{\perp}},\quad\mbox{\boldmath$\mathscr{W}$}_{\overline{\mathcal{N}}}^{(i)}=(\mbox{\boldmath$\mathscr{W}$}_{(i)})_{\overline{\mathcal{N}}_{i}},\quad\mbox{\boldmath$\mathscr{W}$}_{\overline{\mathcal{N}}^{\perp}}^{(i)}=(\mbox{\boldmath$\mathscr{W}$}_{(i)})_{\overline{\mathcal{N}}_{i}^{\perp}}\\\
&\mbox{\boldmath$\mathscr{W}$}_{\mathcal{M}}^{(k)}=(\mbox{\boldmath$\mathscr{W}$}_{[I_{k}]})_{\mathcal{M}_{k}},\quad\mbox{\boldmath$\mathscr{W}$}_{\mathcal{M}^{\perp}}^{(k)}=(\mbox{\boldmath$\mathscr{W}$}_{[I_{k}]})_{\mathcal{M}_{k}^{\perp}},\quad\mbox{\boldmath$\mathscr{W}$}_{\overline{\mathcal{M}}}^{(k)}=(\mbox{\boldmath$\mathscr{W}$}_{[I_{k}]})_{\overline{\mathcal{M}}_{k}},\quad\mbox{\boldmath$\mathscr{W}$}_{\overline{\mathcal{M}}^{\perp}}^{(k)}=(\mbox{\boldmath$\mathscr{W}$}_{[I_{k}]})_{\overline{\mathcal{M}}_{k}^{\perp}},\end{split}$
(B.5)

where $i=1,\dots,2d$ and $k=1,\dots,2^{d-1}$. Based on the subspaces defined
in (B.3) and (B.4), we can define the restricted error sets corresponding to
the three regularized estimators as follows.

######  Definition 1.

The restricted error set corresponding to $\overline{\mathcal{M}}$ is defined
as

$\mathbb{C}_{\textup{SSN}}(\overline{\mathcal{M}}):=\left\\{\boldsymbol{\Delta}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{2d}}:\sum_{k=1}^{2^{d-1}}\|\boldsymbol{\Delta}_{\overline{\mathcal{M}}^{\perp}}^{(k)}\|_{*}\leq
3\sum_{k=1}^{2^{d-1}}\|\boldsymbol{\Delta}^{(k)}_{\overline{\mathcal{M}}}\|_{*}+4\sum_{k=1}^{2^{d-1}}\|\mbox{\boldmath$\mathscr{A}$}^{(k)}_{\mathcal{M}^{\perp}}\|_{*}\right\\}.$
(B.6)

The restricted error set corresponding to $\overline{\mathcal{N}}$ is defined
as

$\mathbb{C}_{\textup{SN}}(\overline{\mathcal{N}}):=\left\\{\boldsymbol{\Delta}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{2d}}:\sum_{i=1}^{2d}\|\boldsymbol{\Delta}_{\overline{\mathcal{N}}^{\perp}}^{(i)}\|_{*}\leq
3\sum_{i=1}^{2d}\|\boldsymbol{\Delta}^{(i)}_{\overline{\mathcal{N}}}\|_{*}+4\sum_{i=1}^{2d}\|\mbox{\boldmath$\mathscr{A}$}^{(i)}_{\mathcal{N}^{\perp}}\|_{*}\right\\}.$
(B.7)

The restricted error set corresponding to $\overline{\mathcal{M}}_{1}$ is
defined as

$\mathbb{C}_{\textup{MN}}(\overline{\mathcal{M}}_{1}):=\left\\{\boldsymbol{\Delta}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{2d}}:\|\boldsymbol{\Delta}_{\overline{\mathcal{M}}^{\perp}}^{(1)}\|_{*}\leq
3\|\boldsymbol{\Delta}^{(1)}_{\overline{\mathcal{M}}}\|_{*}+4\|\mbox{\boldmath$\mathscr{A}$}^{(1)}_{\mathcal{M}^{\perp}}\|_{*}\right\\}.$
(B.8)

The first lemma shows that if the tuning parameter is well chosen for each
regularized estimator, the estimation error belongs to the corresponding
restricted error set.

######  Lemma B.1.

For the SSN estimator, if the regularization parameter
$\lambda_{\textup{SSN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\textup{SSN}^{*}}$,
the error
$\boldsymbol{\Delta}_{\textup{SSN}}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SSN}}-\mbox{\boldmath$\mathscr{A}$}$
belongs to the set $\mathbb{C}_{\textup{SSN}}(\overline{\mathcal{M}})$.

For the SN estimator, if the regularization parameter
$\lambda_{\textup{SN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\textup{SN}^{*}}$,
the error
$\boldsymbol{\Delta}_{\textup{SN}}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{SN}}-\mbox{\boldmath$\mathscr{A}$}$
belongs to the set $\mathbb{C}_{\textup{SN}}(\overline{\mathcal{N}})$.

For the MN estimator, if the regularization parameter
$\lambda_{\textup{MN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\textup{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})\textup{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})^{\top}\|_{*}$,
the error
$\boldsymbol{\Delta}_{\textup{MN}}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\textup{MN}}-\mbox{\boldmath$\mathscr{A}$}$
belongs to the set $\mathbb{C}_{\textup{MN}}(\overline{\mathcal{M}}_{(1)})$.

Following Negahban and Wainwright, (2012) and Negahban et al., (2012), a
restricted strong convexity (RSC) condition for the square loss function can
be defined as follows.

######  Definition 2.

The loss function satisfies the RSC condition with curvature
$\alpha_{\textup{RSC}}>0$ and restricted error set $\mathbb{C}$, if

$\frac{1}{T}\sum_{t=1}^{T}\|\langle\boldsymbol{\Delta},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}\geq\alpha_{\textup{RSC}}\|\boldsymbol{\Delta}\|_{\textup{F}}^{2},~{}~{}~{}\forall\boldsymbol{\Delta}\in\mathbb{C}.$
(B.9)

Based on the restricted error sets and RSC conditions, the estimation errors
have the following deterministic upper bounds.

######  Lemma B.2.

Suppose that $\lambda_{\textup{SSN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\textup{SSN}^{*}}$,
the RSC condition holds with the parameter $\alpha_{\textup{RSC}}$ and
restricted error set $\mathbb{C}_{\textup{SSN}}(\overline{\mathcal{M}})$, and
$\mbox{\boldmath$\mathscr{A}$}_{[I_{k}]}\in\mathbb{B}_{q}(s_{q}^{(k)};p,p)$
for some $q\in[0,1)$ and all $k=1,\dots,2^{d-1}$,

$\|\boldsymbol{\Delta}_{\textup{SSN}}\|_{\textup{F}}\lesssim\sqrt{s_{q}}\left(\frac{2^{d-1}\lambda_{\textup{SSN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2},$
(B.10)

where $s_{q}=2^{1-d}\sum_{k=1}^{2^{d-1}}s_{q}^{(k)}$.

Suppose that $\lambda_{\textup{SN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\textup{SN}^{*}}$,
the RSC condition holds with the parameter $\alpha_{\textup{RSC}}$ and
restricted error set $\mathbb{C}_{\textup{SN}}(\overline{\mathcal{N}})$, and
$\mbox{\boldmath$\mathscr{A}$}_{(i)}\in\mathbb{B}_{q}(r_{q}^{(i)};p_{i},p_{-i}p)$
for some $q\in[0,1)$ and all $i=1,\dots,2d$,

$\|\boldsymbol{\Delta}_{\textup{SN}}\|_{\textup{F}}\lesssim\sqrt{r_{q}}\left(\frac{2d\cdot\lambda_{\textup{SN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2},$
(B.11)

where $r_{q}=(2d)^{-1}\sum_{i=1}^{2d}r_{q}^{(i)}$.

Suppose that $\lambda_{\textup{MN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\textup{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})\textup{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})\|_{*}$,
the RSC condition holds with the parameter $\alpha_{\textup{RSC}}$ and
restricted error set $\mathbb{C}_{\textup{MN}}(\overline{\mathcal{M}}_{1})$,
and
$\mbox{\boldmath$\mathscr{A}$}_{[S_{1}]}\in\mathbb{B}_{q}(s_{q}^{(1)};p,p)$
for some $q\in[0,1)$,

$\|\boldsymbol{\Delta}_{\textup{MN}}\|_{\textup{F}}\lesssim\sqrt{s_{q}^{(1)}}\left(\frac{\lambda_{\textup{MN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2}.$
(B.12)

Note that Lemma B.2 is deterministic and the radius $s_{q}$, $r_{q}$, and
$s_{q}^{(1)}$ can also diverge to infinity.

####  B.1.2 Stochastic Analysis

We continue with the stochastic analysis to show that the deviation bound and
the RSC condition hold simultaneously with high probability.

######  Lemma B.3 (Deviation bound).

Suppose that Assumptions 1 and 2 hold. If $T\gtrsim p$ and
$\lambda_{\textup{SSN}}\gtrsim\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, with
probability at least $1-\exp[-C(p-d)]$,

$\left\|\frac{1}{T}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\right\|_{\textup{SSN}^{*}}\leq\frac{\lambda_{\textup{SSN}}}{4}$
(B.13)

where
$M_{1}=\lambda_{\max}(\boldsymbol{\Sigma_{e}})/\mu^{1/2}_{\min}(\mathcal{A})$.

If $T\gtrsim\max_{1\leq i\leq d}p_{-i}p$ and
$\lambda_{\textup{SN}}\gtrsim\kappa^{2}M_{1}d^{-2}\sum_{i=1}^{d}\sqrt{p_{-i}p/T}$,
with probability at least $1-2\sum_{i=1}^{d}\exp(-Cp_{-i}p)$,

$\left\|\frac{1}{T}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\right\|_{\textup{SN}^{*}}\leq\frac{\lambda_{\textup{SN}}}{4}.$
(B.14)

Moreover, if $T\gtrsim p$ and
$\lambda_{\textup{MN}}\gtrsim\kappa^{2}M_{1}\sqrt{p/T}$, with probability at
least $1-\exp(-Cp)$,

$\left\|\frac{1}{T}\sum_{t=1}^{T}\textup{vec}(\mbox{\boldmath$\mathscr{Y}$}_{t-1})\textup{vec}(\mbox{\boldmath$\mathscr{E}$}_{t})^{\top}\right\|_{*}\leq\frac{\lambda_{\textup{MN}}}{4}.$
(B.15)

Next, we prove the restricted strong convexity for regularized estimators.
According to Lemma B.3, we need the sample size $T\gtrsim p$ for all three
estimators. In this case, we can establish the strong convexity condition that
is stronger than the RSC condition.

######  Lemma B.4 (Strong convexity).

Under Assumptions 1 and 2, for
$T\gtrsim\max(\kappa^{2},\kappa^{4})M_{2}^{-2}p$, with probability at least
$1-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$,

$\frac{1}{T}\sum_{t=1}^{T}\|\langle\boldsymbol{\Delta},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}\geq\alpha_{\textup{RSC}}\|\boldsymbol{\Delta}\|_{\textup{F}}^{2},$
(B.16)

where
$M_{2}=[\lambda_{\min}(\boldsymbol{\Sigma}_{\boldsymbol{e}})\mu_{\max}(\mathcal{A})]/[\lambda_{\max}(\boldsymbol{\Sigma}_{\boldsymbol{e}})\mu_{\min}(\mathcal{A})]$
and
$\alpha_{\textup{RSC}}=\lambda_{\min}(\boldsymbol{\Sigma}_{\boldsymbol{e}})/(2\mu_{\max}(\mathcal{A}))$.

###  B.2 Proofs of Theorems 2–5

###### Proof of Theorems 2 and 3.

Theorems 2 and 3 can be proved based on Lemmas B.2–B.4 following the same line
of the proof of Theorem 4 given below. Therefore, we omit the details here. ∎

###### Proof of Theorem 4.

The proof of Theorem 4 has been split into Lemmas B.2–B.4. By Lemma B.2, for
deterministic realization with sample size $T$ of a tensor autoregressive
process, if we choose $\lambda_{\text{SSN}}\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\text{SSN}^{*}}$
and RSC condition holds for the square loss with the parameter
$\alpha_{\text{RSC}}$, the following error upper bound can be established

$\|\boldsymbol{\Delta}\|_{\textup{F}}\lesssim\sqrt{s_{q}}\left(\frac{2^{d-1}\lambda_{\textup{SSN}}}{\alpha_{\textup{RSC}}}\right)^{1-q/2}.$
(B.17)

Denote the events $E_{1}(\beta)=\\{\beta\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\text{SSN}^{*}}\\}$
and
$E_{2}(\alpha)=\\{\lambda_{\min}(\boldsymbol{X}\boldsymbol{X}^{\top}/T)\geq\alpha\\}$.
If we take $\lambda_{\text{SSN}}\gtrsim\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, it
suffices to show that $E_{1}(C\kappa^{2}M_{1}2^{1-d}\sqrt{p/T})$ and
$E_{2}(\alpha_{\text{RSC}}/2)$ occur simultaneously with high probability.

By Lemma B.3, when $T\gtrsim p$,

$\left\|\frac{1}{T}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\right\|_{\text{SSN}^{*}}\lesssim\kappa^{2}M_{1}2^{1-d}\sqrt{\frac{p}{T}}$
(B.18)

with probability at least $1-\exp[-C(p-d)]$.

By Lemma B.4, when $T\gtrsim\max(\kappa^{2},\kappa^{4})M_{2}^{-2}p$, for any
$\boldsymbol{\Delta}\in\mathbb{R}^{p_{1}\times\cdots\times p_{2d}}$,

$\frac{1}{T}\sum_{t=1}^{T}\|\langle\boldsymbol{\Delta},\mbox{\boldmath$\mathscr{Y}$}_{t-1}\rangle\|_{\text{F}}^{2}\geq\frac{\lambda_{\min}(\boldsymbol{\Sigma_{e}})}{2\mu_{\max}(\mathcal{A})}\|\boldsymbol{\Delta}\|_{\text{F}}^{2}$
(B.19)

with probability at least $1-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$.

Hence, when $T\gtrsim[1+\max(\kappa^{2},\kappa^{4})M_{2}^{-2}]p$ and
$\lambda\gtrsim\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, with probability at least
$1-\exp[-C(p-d)]-\exp[-C\min(\kappa^{-2},\kappa^{-4})M_{2}^{2}p]$, the
condition $\lambda\geq
4\|T^{-1}\sum_{t=1}^{T}\mbox{\boldmath$\mathscr{Y}$}_{t-1}\circ\mbox{\boldmath$\mathscr{E}$}_{t}\|_{\text{SSN}^{*}}$
and the RSC condition with the parameter
$\alpha_{\text{RSC}}=\lambda_{\min}(\boldsymbol{\Sigma_{e}})/\mu_{\max}(\mathcal{A})$
hold. ∎

###### Proof of Theorem 5.

Theorem 4 gives the Frobenius estimation error bound. For simplicity, we write
$\mbox{\boldmath$\mathscr{\widehat{A}}$}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{SSN}}$
and
$\mbox{\boldmath$\mathscr{\widetilde{A}}$}=\mbox{\boldmath$\mathscr{\widehat{A}}$}_{\text{TSSN}}$
in this proof. By definition, for any tensor
$\mbox{\boldmath$\mathscr{A}$}\in\mathbb{R}^{p_{1}\times\cdots\times p_{2d}}$,

$\|\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}^{2}=\|\mbox{\boldmath$\mathscr{A}$}_{(i)}\|_{\text{F}}^{2}=\sum_{j=1}^{p_{i}}\sigma_{j}^{2}(\mbox{\boldmath$\mathscr{A}$}_{(i)}),~{}~{}i=1,2,\dots,2d.$
(B.20)

In other words, the Frobenius norm of the error tensor is equivalent to the
$\ell_{2}$ norm of singular values of the one-mode matricization. By Mirsky’s
singular value inequality (Mirsky,, 1960),

$\sum_{j=1}^{p_{i}}[\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})-\sigma_{j}(\mbox{\boldmath$\mathscr{A}$}_{(i)})]^{2}\leq\sum_{j=1}^{p_{i}}\sigma_{j}^{2}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)}-\mbox{\boldmath$\mathscr{A}$}_{(i)})=\|\mbox{\boldmath$\mathscr{\widehat{A}}$}-\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}^{2},~{}~{}i=1,2,\dots,2d.$
(B.21)

Obviously, the $\ell_{\infty}$ error bound is smaller than the $\ell_{2}$
error bound, so it follows the same upper bound. By Theorem 4, when
$\lambda_{\textup{SSN}}\asymp\kappa^{2}M_{1}2^{1-d}\sqrt{p/T}$, with
probability approaching one,

$\displaystyle\max_{1\leq i\leq 2d}\max_{1\leq j\leq
p_{i}}|\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})-\sigma_{j}(\mbox{\boldmath$\mathscr{A}$}_{(i)})|\leq$
$\displaystyle\max_{1\leq i\leq
2d}\left\\{\sum_{j=1}^{p_{i}}[\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})-\sigma_{j}(\mbox{\boldmath$\mathscr{A}$}_{(i)})]^{2}\right\\}^{1/2}$
$\displaystyle\leq$
$\displaystyle\|\mbox{\boldmath$\mathscr{\widehat{A}}$}-\mbox{\boldmath$\mathscr{A}$}\|_{\text{F}}\lesssim\frac{\kappa^{2}M_{1}}{\alpha_{\text{RSC}}}\sqrt{\frac{s_{0}p}{T}}.$
(B.22)

Therefore, by Assumption 3, as $T\rightarrow\infty$,

$\gamma\gg\max_{1\leq i\leq 2d}\max_{1\leq j\leq
p_{i}}|\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})-\sigma_{j}(\mbox{\boldmath$\mathscr{A}$}_{(i)})|.$
(B.23)

Then, for any $j>r_{i}$, since
$\sigma_{j}(\mbox{\boldmath$\mathscr{A}$}_{(i)})=0$, we have
$\gamma\gg\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})$. Thus,
for all $i=1,\dots,2d$,
$\sigma_{j}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})$ will be truncated
for all $j>r_{i}$. Meanwhile, by Assumption 3 and (B.22), we have
$\sigma_{r_{i}}(\mbox{\boldmath$\mathscr{\widehat{A}}$}_{(i)})>\gamma$ for $T$
sufficiently large, for all $i=1,\dots,2d$. Therefore, the rank selection
consistency of the truncated estimator $\mathscr{\widetilde{A}}$ can be
established.

Denote the event