# Guaranteed Functional Tensor Singular Value Decomposition3

Rungang Han1, Pixu Shi2, and Anru R. Zhang1,2

1Department of Statistics, University of Wisconsin-Madison
2Department of Biostatistics & Bioinformatics, Duke University

(())

###### Abstract

This paper introduces the functional tensor singular value decomposition
(FTSVD), a novel dimension reduction framework for tensors with one functional
mode and several tabular modes. The problem is motivated by high-order
longitudinal data analysis. Our model assumes the observed data to be a random
realization of an approximate CP low-rank functional tensor measured on a
discrete time grid. Incorporating tensor algebra and the theory of Reproducing
Kernel Hilbert Space (RKHS), we propose a novel RKHS-based constrained power
iteration with spectral initialization. Our method can successfully estimate
both singular vectors and functions of the low-rank structure in the observed
data. With mild assumptions, we establish the non-asymptotic contractive error
bounds for the proposed algorithm. The superiority of the proposed framework
is demonstrated via extensive experiments on both simulated and real data.

33footnotetext: The research of R. Han and A. R. Zhang was supported in part
by NSF Grant CAREER-1944904 and NIH Grant R01 GM131399.

##  1 Introduction

In recent decades, the analysis of tensor data has become an active research
topic. Datasets in the form of high-order tensors or multiway arrays arise
from various scientific applications, such as neuroimaging (Li et al.,, 2018;
Zhang et al.,, 2019), applied physics (Zhang et al.,, 2020), network modeling
(Bi et al.,, 2018; Jing et al.,, 2020), microscopy imaging (Zhang et al.,,
2020), and longitudinal microbiome study (Martino et al.,, 2021). Such high-
order data pose significant challenges in both theoretical analysis and
computation implementation due to their complicated structures and a large
number of entries involved, making it inappropriate to extend many existing
methods for matrices to the analysis of tensor data.

In real applications, the different modes of tensor datasets (or different
directions that the tensor arrays align in) can come in various formats. Two
prominent formats are (1) tabular modes, such as subject ID, genomics ID,
treatments, where shuffling the indices does not essentially change the data
structure; (2) functional modes, such as time, location, spectrum (in the
hyperspectral imaging), where the order of indices exhibits structures such as
continuity. Here we provide two scenarios to illustrate tensor data with both
types of modes.

  * •

_Multivariate functional data analysis._ Data with different functional
features appear in various applications and can be formatted into a tensor
with two tabular modes representing units and variables, and a functional mode
representing time point. One example is the longitudinal microbiome studies,
where microbiome samples are taken from multiple subjects (units) at multiple
time points to study the abundance of bacteria (variables) over time.
Depending on what taxonomic level is being studied, there can be hundreds and
thousands of bacterial taxa in the feature mode and many of these taxa have
strong correlation in their abundance.

  * •

_Dynamic networks._ In network analysis, one often observes multiple snapshots
of dynamic networks. The adjacency matrix in each snapshot can be stacked
together into an adjacency tensor, where two tabular modes correspond to the
vertices of the network, and a functional mode corresponds to time. The
dimension of the tabular modes, i.e., number of nodes, can be as large as
thousands, but the hidden community structure of the network can often be
modeled by low-rankness of adjacency tensor (Lei et al.,, 2020; Jing et al.,,
2020).

Just as in the examples above, the tabular modes may be very high-dimensional.
The functional mode is often continuous, which is essentially infinite-
dimensional. Therefore, dimension reduction tools are crucial for the
visualization and analysis of high-dimensional tensor data with both tabular
modes and functional modes. Various lines of methods can be applied, whereas
one prominent framework is the multivariate/multilevel functional principal
component analysis (MFPCA) in the functional data analysis (FDA) literature. A
more detailed discussions on the literature of multivariate functional
analysis is given in Section 1.2. The classic MFPCA methods often focus on a
moderately large number of variables and dimension reduction on the functional
mode. The resulting number of estimated eigenfunctions is often in the same
order as the number of variables, making the eigenfunctions hard to interpret
when the variable mode is high-dimensional. Another important class of works
utilize the high-order structure of the tensor data via different types of
low-rank tensor decomposition models, such as (sparse) CP low-rankness
(Anandkumar et al., 2014b, ; Sun et al.,, 2017), (sparse) Tucker low-rankness
(Zhang and Xia,, 2018; Zhang and Han,, 2019), tensor-train (Oseledets,, 2011;
Zhou et al., 2020b, ), or decomposable tensor covariance (Dawid,, 1981;
Tsiligkaridis and Hero,, 2013; Zhou,, 2014). While the low-rank models achieve
efficient dimension reduction and are interpretable when the modes are
tabular, they do not directly characterize or utilize the information included
in the functional modes such as time and location. Therefore, new methods and
theory are needed to tackle both the high-dimensional and functional aspects
of functional tensor data.

In this paper, we introduce a new framework for the dimension reduction of
functional tensor data which we refer to as _Functional Tensor Singular Value
Decomposition_ (FTSVD). Suppose ${\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times
p_{2}\times\mathcal{T}}$ is the underlying function of interest, where
$\mathcal{T}\subset{\mathbb{R}}$ is some compact set on the real axis. For any
fixed $(i,j)$ pair, ${\mathbfcal{Y}}_{ijt}$ is a function of
$t\in\mathcal{T}$. We aim to decompose ${\mathbfcal{Y}}$ into the following
canonical polyadic (CP) format:

${\mathbfcal{Y}}={\mathbfcal{X}}+{\mathbfcal{Z}},\quad\text{where}\quad{\mathbfcal{X}}=\sum_{l=1}^{r}\lambda_{l}a_{l}\otimes
b_{l}\otimes\xi_{l}\quad\text{or
equivalently}\quad{\mathbfcal{X}}_{ijt}=\sum_{l=1}^{r}\lambda_{l}(a_{l})_{i}(b_{l})_{j}\xi_{l}(t).$
(1)

Here, $a_{l}\in\mathbb{R}^{p_{1}}$, $b_{l}\in\mathbb{R}^{p_{2}}$ are singular
vectors of tabular modes; $\xi_{l}:\mathcal{T}\to\mathbb{R}$ is the singular
function, which corresponds to the eigenfunction in functional PCA;
${\mathbfcal{Z}}$ is the remainder term of the rank-$r$ CP decomposition of
${\mathbfcal{Y}}$. We assume the observed data $\widetilde{{\mathbfcal{Y}}}$
to be ${\mathbfcal{Y}}$ measured over a discrete grid
$\\{s_{k}\\}_{k=1}^{n}\subset\mathcal{T}$ with observational noise, i.e.,

$\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times
n},\qquad\widetilde{\mathbfcal{Y}}_{ijk}={\mathbfcal{Y}}_{ijs_{k}}+\varepsilon_{ijk}.$
(2)

Without loss of generality, we assume that $\mathcal{T}=[0,1]$ throughout the
paper.

The proposed model setting is flexible and adaptive. For example, in the
scenario of multivariate functional data analysis, different from many
existing literature on (multivariate) functional PCA (Di et al.,, 2009; Happ
and Greven,, 2018; Wang et al., 2020c, ), which assumed the samples are i.i.d.
distributed and try to estimate the covariance function, our model allows
significant heterogeneity among samples (which is characterized by the
singular vector of the sample mode). In the scenario of dynamic network, to
model the temporal connectivity dynamics, one can assume the latent factors
for all $p_{1}$ vertices
$Z=[\sqrt{\lambda_{1}}a_{1},\ldots,\sqrt{\lambda_{r}}a_{r}]$ are time-
invariant, and the similarity (connectivity probability) matrix at time
$s_{k}$ can be modeled as a weighted inner product:
$Z\Lambda_{s_{k}}Z^{\top}$, where the weight matrix $\Lambda_{s_{k}}={\rm
diag}\left(\xi_{1}(s_{k}),\ldots,\xi_{r}(s_{k})\right)$. The model can be
reduced to the regular symmetric CP decomposition with $a_{i}=b_{i}$. The
readers are referred to surveys on multilayer network in Goldenberg et al.,
(2009); Kivelä et al., (2014). Additional applications that can be formatted
under this framework include (hyper)imaging analysis (Zhang et al.,, 2020),
online recommendation system (Bi et al.,, 2018), etc.

Departing from the existing literature, the proposed framework FTSVD can
provide important new insights to the analysis of data tensors with both
tabular and functional modes. Firstly, through the lens of functional tensor
SVD, we can describe the trend of all variables ${\mathbfcal{X}}_{ij\cdot}$
over time using only a small set of nonparametric singular functions
$\xi_{l}(t)$. Such low-dimensional representation of the functional mode is
more interpretable than high-dimensional ones and may cast light on the
driving factors behind the trend of variables. Secondly, the singular vectors
$a_{l}$ and $b_{l}$ provide a low-dimensional depiction of the information in
the tabular modes, which can be used to improve the efficiency and
interpretability of subsequent analysis such as clustering, classification,
and regression.

###  1.1 Contributions

The main contributions of this paper are three-fold. First, to the best of our
knowledge, we are the first to study multivariate/multilevel functional data
analysis through a lens of a low-rank tensor decomposition model. We are also
the first to adopt the reproducing kernel Hilbert space (RKHS) framework for
the singular function in the low-rank decomposition. Compared to recently
emerged state-of-the-art tabular tensor models, the proposed model can better
characterize the longitudinal behaviors of the data when there exist temporal
modes. Compared to most classic FDA literature, our model does not assume the
samples to be i.i.d. Instead, our method may provide new insight on how to
analyze heteroskedastic functional effects for mixed-population samples.

Secondly, we propose a new power iteration algorithm for the estimation of our
model. In each iteration step, we use projection to update the estimators for
the tabular components and constrained/regularized empirical risk minimization
(ERM) to update the estimators for the functional component, respectively. The
novel inclusion of multiple-times ERM in our iterative algorithm imposes new
challenges on both statistical theory and numeric optimization.

Thirdly, we establish the finite-sample statistical error bounds for the
estimation of both singular vectors and functions (17). Our error bound
includes contribution from the remainder term of CP approximation (Eq. (1))
and observational noises (Eq. (2)) relative to the signal level in
${\mathbfcal{X}}$. It also reflects the effect of time grid density and
incoherence condition of the CP decomposition. To the best of our knowledge,
this is the first contractive error bound for iterative algorithms in
functional data analysis.

###  1.2 Related Work

This work is related to a series of papers on _low-rank tensor decomposition_
from statistics and machine learning, which aims at recovering the low-
dimensional linear representations of the noisy tabular tensors (De Lathauwer
et al., 2000a, ; De Lathauwer et al., 2000b, ; Richard and Montanari,, 2014;
Anandkumar et al., 2014a, ). Commonly-used low-rank models include CP
(Anandkumar et al., 2014b, ), Tucker (Zhang and Xia,, 2018), and tensor train
(Zhou et al., 2020b, ), where various power iteration methods were proposed on
the singular vector (or principal component) estimation with statistical
consistency guarantees. Motivated by real data applications in practice,
related statistical models were further studied under non-Gaussian or missing-
data scenarios. For example, Hong et al., (2020); Han et al., (2021)
considered the likelihood-based generalized low-rank tensor decomposition
under general exponential family; Cai et al., (2019); Tong et al., (2021)
studied the low-rank tensor estimation under incomplete measurements.

In addition, various dimension-reduced structures are imposed on the singular
vectors or subspace of the low-rank tensor decomposition model to better
capture other intrinsic properties of the data, such as sparsity (Sun et al.,,
2017; Zhang and Han,, 2019), blocking (Han et al.,, 2020), non-negativity (Xu
and Yin,, 2013) and many others. The most relevant paper to our work is the
spatial/temporal structure, which is incorporated to characterize
time/location-varying patterns for one or more of the tensor modes. There are
various formulations for modeling such structures based on different real data
motivations. For example, Sun and Li, (2019); Zhou et al., 2020a  considered a
fusion structure and assumed the time-varying trend to be piece-wise constant;
Han et al., (2019) considered the multidimensional isotonic regression, where
there are monotonous trends in each tensor mode. In most of such literature,
the temporal structures are still assumed on tabular tensors and make the
corresponding models restrictive. In contrast, the target of our framework is
the underlying non-parametric functional structures, which reflects the trend
of observed data over time.

Another related topic is the _functional data analysis_ (FDA), which is a
popular branch of statistical research. The readers are referred to surveys
and books in Ramsay and Silverman, (2006); Ferraty and Vieu, (2006); Brezis,
(2010); Hsing and Eubank, (2015); Wang et al., (2016). Starting from samples
with univariate functions, traditional FDA considers various statistical
learning tasks, such as mean estimation (Cai and Yuan,, 2011), regression (Yao
et al., 2005b, ; Cai and Yuan,, 2012; Raskutti et al.,, 2012; Sun et al.,,
2018), covariance function estimation (Rice and Silverman,, 1991; Li and
Hsing,, 2010), functional mixed effect model (Guo,, 2002; Morris and Carroll,,
2006), principal component analysis (PCA) (James et al.,, 2000; Yao et al.,
2005a, ), functional clustering (Abraham et al.,, 2003; Chiou and Li,, 2007;
Chiou,, 2012), and time series analysis (Aue et al.,, 2015; van Delft et al.,,
2017; Wang et al., 2020a, ).

Problems with multiple functional variables breed numerous works on
multivariate functional models, which is pioneered by Dubin and Müller,
(2005); Crainiceanu et al., (2009); Greven et al., (2010) among many others.
See Verbeke et al., (2014) for a survey on multivariate functional data
analysis. Under the multivariate regime, Di et al., (2009, 2014); Hasenstab et
al., (2017) considered the multilevel/multidimensional functional data PCA.
Wang et al., 2020c  studied the low-rank covariance estimation for
multidimensional functional data. Fan et al., (2015) considered the functional
additive regression method for high-dimensional functional regression. Hu and
Yao, (2021) considered the dimension reduction via dynamic principal subspace
for multivariate functional data. Allen, (2013) studied the general estimation
framework for arbitrary $p$-dimensional functions under regular or irregular
grid designs. Chen et al., (2017) studied the functional data analysis for
scenarios where the observations are functional at each location. The works on
multivariate functional principal component analysis (Greven et al.,, 2010;
Zipunnikov et al.,, 2011; Allen,, 2013; Chiou et al.,, 2014; Scheffler et
al.,, 2020; Berrendero et al.,, 2011; Happ and Greven,, 2018) also address
dimension reduction of multivariate functional data. Among these literature,
Happ and Greven, (2018); Wang et al., 2020c  are the most relevant to our
work. In particular, Wang et al., 2020c  focused on i.i.d. $p$-dimensional
function objects, and they impose tensor low-rankness on the corresponding
$(2p)$-dimensional covariance functions; Happ and Greven, (2018) considered
the multivariate functional PCA for i.i.d. samples from some multivariate
functional population and the data can be formulated as a typical functional
tensor stacked by “sample-feature-function.” Both papers assumed i.i.d.
samples and focused on the statistical inference in the covariance function.
Compared to these methods, our model does not assume independence of the
signal tensor along any mode. Such flexible modeling is able to encompass a
range of problems on heterogeneous or heteroskedastic data. The resulting
estimated singular vectors can also be used to perform many important tasks,
such as clustering, classification, and regression.

Our framework is related to a line of recent works on matrix/tensor-valued
times series arising from the applications in financial econometrics. See,
e.g., Wang et al., (2019); Chen et al., (2020, 2021) and the references
therein. Departing from their work that focuses on the dynamic and evolution
pattern of the functions over time, our work can reveal the structures under
each mode, including clusters in the tabular modes and variations in the
functional mode that are associated with other factors.

###  1.3 Organization

The rest of this paper is organized as follows. After an introduction of the
notation and preliminaries on tensor algebra and reproducing kernel Hilbert
space in Section 2.1, we present the statistical model of FTSVD in Section
2.2. We present the RKHS-based constrained power iteration method for
parameter estimation in Section 3. The theoretical guarantee is provided in
Section 4. Numeric experiments are collected in Section 5, with simulation
study and real data analysis presented in Sections 5.2 and 5.3, respectively.
Further discussions are given in Section 6. The proofs of main results are
presented in Section 7 and all the other proofs are collected in the
supplementary materials.

##  2 Low-rank Functional Tensor Decomposition Model

###  2.1 Notation and Preliminaries

We use the lowercase letters, e.g., $a,b,u,v$, to denote scalars or vectors.
For any $a,b\in{\mathbb{R}}$, let $a\wedge b$ and $a\vee b$ be the minimum and
maximum of $a$ and $b$, respectively. For a vector $u\in{\mathbb{R}}^{n}$,
$\|u\|_{2}$ denotes its Euclidean norm. The unit sphere in ${\mathbb{R}}^{n}$
is denoted as ${\mathbb{S}}^{n-1}$. Matrices are denoted as uppercase letters
such as $\boldsymbol{A}$ and $\boldsymbol{B}$. In addition, we let
${\mathbb{O}}_{p,r}$ be the collection of all $p$-by-$r$ matrices with
orthonormal columns:
${\mathbb{O}}_{p,r}:=\\{\boldsymbol{U}\in{\mathbb{R}}^{p\times
r}:\boldsymbol{U}^{\top}\boldsymbol{U}=\boldsymbol{I}\\}$, where
$\boldsymbol{I}$ is the identity matrix. We also denote
${\mathbb{O}}_{r}:={\mathbb{O}}_{r,r}$ as the set of all $r$-dimensional
orthogonal matrices. Let
$\lambda_{1}(\boldsymbol{A})\geq\lambda_{2}(\boldsymbol{A})\geq\cdots\geq 0$
be the singular values of $\boldsymbol{A}$ in descending order and let ${\rm
SVD}_{r}(\boldsymbol{A})$ be the matrix comprised of the top $r$ left singular
vectors of $\boldsymbol{A}$. Let
$\|\boldsymbol{A}\|=\lambda_{1}(\boldsymbol{A})$ and $\|\boldsymbol{A}\|_{\rm
F}=\sqrt{\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\boldsymbol{A}_{ij}^{2}}=\sqrt{\sum_{i=1}^{p_{1}\wedge
p_{2}}\lambda_{i}^{2}(\boldsymbol{A})}$ be the matrix spectral and Frobenius
norms, respectively. For any matrix
$\boldsymbol{A}=[a_{1},\ldots,a_{J}]\in\mathbb{R}^{I\times J}$ and
$\boldsymbol{B}\in\mathbb{R}^{K\times L}$, the _Kronecker product_ is defined
as the $(IK)$-by-$(JL)$ matrix
$\boldsymbol{A}\otimes\boldsymbol{B}=[a_{1}\otimes\boldsymbol{B}\cdots
a_{J}\otimes\boldsymbol{B}]$.

Without loss of generality, we assume the function domain of interest is
$\mathcal{T}=[0,1]$ and denote $\mathcal{L}^{2}([0,1])$ as the functional
space of all square-integrable functions on $[0,1]$, i.e.,

$\mathcal{L}^{2}\left([0,1]\right)=\left\\{f:[0,1]\rightarrow{\mathbb{R}},\|f\|_{\mathcal{L}^{2}}^{2}<\infty\right\\},\quad\text{where}\quad\|f\|_{\mathcal{L}^{2}}:=\left(\int_{0}^{1}f^{2}(t)dt\right)^{1/2}.$

Denote the inner product of any two functions $f,g\in\mathcal{L}^{2}([0,1])$
as $\langle f,g\rangle_{\mathcal{L}^{2}}=\int_{0}^{1}f(t)g(t)dt$. For a
sequence $\\{s_{k}\\}_{k=1}^{n}$ in $[0,1]$, we denote
$f_{n}:=(f(s_{1}),\ldots,f(s_{n}))\in{\mathbb{R}}^{n}$ and
$\|f\|_{n}:=\sqrt{\sum_{k=1}^{n}f^{2}(s_{k})/n}$. We use
$C,C_{0},C_{1},\ldots$ and $c,c_{0},c_{1},\ldots$ to represent generic large
and small positive constants, respectively. The actual values of these generic
symbols may differ from line to line. We denote $a\lesssim b$ if $a\leq Cb$
for a constant $C>0$ that does not depend on other model parameters; we say
$a\asymp b$ if $a\lesssim b$ and $b\lesssim a$ both hold.

Throughout this paper, the tensors are denoted by uppercase calligraphy
letters, such as ${\mathbfcal{X}},{\mathbfcal{Y}},{\mathbfcal{Z}}$. For a
tabular tensor ${\mathbfcal{A}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times
p_{3}}$, let ${\mathbfcal{A}}_{ijk}$ be the $(i,j,k)$th entry. For any
$u\in{\mathbb{R}}^{p_{1}},v\in{\mathbb{R}}^{p_{2}},w\in{\mathbb{R}}^{p_{3}}$,
the mode-1 tensor-vector product is defined as:

$\begin{split}{\mathbfcal{A}}\times_{1}u\in\mathbb{R}^{p_{2}\times
p_{3}},\qquad({\mathbfcal{A}}\times_{1}u)_{jk}:=\sum_{i=1}^{p_{1}}u_{i}{\mathbfcal{A}}_{ijk}.\end{split}$

The mode-2 and mode-3 tensor-vector products, $\mathcal{A}\times_{2}v$ and
$\mathcal{A}\times_{3}w$ for $v\in\mathbb{R}^{p_{2}}$ and
$w\in\mathbb{R}^{p_{3}}$, can be defined in the parallel way. The
multiplication along different modes is cumulative and commutative. For
example, ${\mathbfcal{A}}\times_{1}u\times_{2}v\in{\mathbb{R}}^{p_{3}}$ with
$\left({\mathbfcal{A}}\times_{1}u\times_{2}v\right)_{k}=\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}u_{i}v_{j}{\mathbfcal{A}}_{ijk}$;
and
${\mathbfcal{A}}\times_{1}u\times_{2}v\times_{3}w=\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{j=1}^{p_{3}}u_{i}v_{j}w_{k}{\mathbfcal{A}}_{ijk}\in{\mathbb{R}}$.
We also introduce the matricization operator that transforms tensors to
matrices. Particularly, the mode-$1$ matricization of
${\mathbfcal{A}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times p_{3}}$ is defined as

$\begin{split}\mathcal{M}_{1}({\mathbfcal{A}})\in{\mathbb{R}}^{p_{1}\times(p_{2}p_{3})},\quad\text{where}\quad[\mathcal{M}_{1}({\mathbfcal{X}})]_{i,j+p_{2}(k-1)}={\mathbfcal{A}}_{ijk}.\end{split}$

In other words, each row of $\mathcal{M}_{k}({\mathbfcal{A}})$ is the
vectorization of a mode-$k$ slice.

Next, we extend the notions of tabular tensors to tensors with a hybrid of
tabular and functional modes. We still denote the tensors with hybrid modes by
uppercase calligraphy letters. Suppose
${\mathbfcal{A}}=\\{{\mathbfcal{A}}_{ijt},i\in[p_{1}],j\in[p_{2}],t\in[0,1]\\}\in\mathbb{R}^{p_{1}\times
p_{2}\times[0,1]}$ is an order-3 tensor with two tabular modes and a
functional mode. The tensor-vector products on the tabular modes are defined
as convention:

${\mathbfcal{A}}\times_{1}u\in\mathbb{R}^{p_{2}\times[0,1]},\qquad({\mathbfcal{A}}\times_{1}u)_{jt}:=\sum_{i=1}^{n}u_{i}{\mathbfcal{A}}_{ijt}.$

Suppose $f\in\mathcal{L}^{2}([0,1])$, we define the tensor-functional product
as

${\mathbfcal{A}}\times_{3}f\in\mathbb{R}^{p_{1}\times
p_{2}},\qquad({\mathbfcal{A}}\times_{3}f)_{ij}:=\int_{0}^{1}{\mathbfcal{A}}_{ijt}f(t)dt.$

The multiple-mode tensor-vector/function multiplication can be similarly
defined by combining the operators of single-mode multiplication. For example,
${\mathbfcal{A}}\times_{1}u\times_{2}v$ yields a function
$\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}u_{i}v_{j}{\mathbfcal{A}}_{ijt}$ and
${\mathbfcal{A}}\times_{1}a\times_{2}b\times_{3}f$ yields a scalar
$\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\int_{0}^{1}u_{i}v_{j}{\mathbfcal{A}}_{ijt}f(t)dt$.

Next, we provide the preliminaries for the reproducing kernel Hilbert space
(RKHS). Consider a Hilbert space $\mathcal{H}\subset\mathcal{L}^{2}([0,1])$
with the associated inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$.
We assume there exists a continuous symmetric positive-semidefinite kernel
function $\mathbb{K}:[0,1]\times[0,1]\rightarrow\mathbb{R}_{+}$ that satisfies
the following RKHS conditions: (1) for any $s\in[0,1]$,
$\mathbb{K}(\cdot,s)\in\mathcal{H}$; (2) for each $g\in\mathcal{H}$,
$g(t)=\langle g,\mathbb{K}(\cdot,t)\rangle_{\mathcal{H}},\leavevmode\nobreak\
\forall t\in[0,1]$. By Mercer’s Theorem (Mercer,, 1909), there exists an
orthonormal basis $\\{\phi_{k}\\}_{k=1}^{\infty}$ of $\mathcal{L}^{2}([0,1])$,
such that $\mathbb{K}(\cdot,\cdot)$ admits the following eigendecomposition:

$\mathbb{K}(s,t)=\sum_{k=1}^{\infty}\mu_{k}\phi_{k}(s)\phi_{k}(t),\leavevmode\nobreak\
\leavevmode\nobreak\ \forall s,t\in[0,1].$

Here, $\mu_{1}\geq\mu_{2}\geq\cdots\geq 0$ are the non-negative eigenvalues of
$\mathbb{K}$. We assume $\mathbb{K}$ has a finite trace norm:
$\sum_{k=1}^{\infty}\mu_{k}<\infty$. Then any $f\in\mathcal{H}$ can be
decomposed on $\\{\phi_{k}\\}_{k=1}^{\infty}$ as

$f(t)=\sum_{k=1}^{\infty}a_{k}\phi_{k}(t),\qquad
a_{k}:=\int_{0}^{1}f(s)\phi_{k}(s)ds,\qquad t\in[0,1].$

The RKHS norm of $f$ can be then represented as
$\|f\|_{\mathcal{H}}=\sqrt{\langle
f,f\rangle_{\mathcal{H}}}=\sqrt{\sum_{k=1}^{\infty}a_{k}^{2}/\mu_{k}}$.
Accordingly, for any two functions $f=\sum_{k=1}^{\infty}a_{k}\phi_{k}$ and
$g=\sum_{k=1}^{\infty}b_{k}\phi_{k}$ in $\mathcal{H}$, the inner product with
respect to $\mathcal{H}$ can be calculate as $\langle
f,g\rangle_{\mathcal{H}}=\sum_{k=1}^{\infty}a_{k}b_{k}/\mu_{k}$.

We assume the following regularity conditions on the kernel function
${\mathbb{K}}$ throughout the paper.

######  Assumption 1 (Kernel Regularity Conditions).

(1) ${\mathbb{K}}$ satisfies
$\max\\{\mu_{1},\sup_{s\in[0,1]}{\mathbb{K}}(s,s)\\}\leq 1$;
(2) there exists an absolute constant $C_{\mathcal{H}}$, which only depends on
$\mathcal{H}$, such that

$\|fg\|_{\mathcal{H}}\leq
C_{\mathcal{H}}\|f\|_{\mathcal{H}}\|g\|_{\mathcal{H}},\qquad\forall
f,g\in\mathcal{H}.$ (3)

Assumption 1 (a) can be ensured by rescaling ${\mathbb{K}}$, i.e., multiplying
${\mathbb{K}}$ by some positive constant. This condition also implies

$\begin{split}\|f\|_{\mathcal{L}^{2}}\leq&\|f\|_{\infty}:=\sup_{t\in[0,1]}|f(t)|=\sup_{s\in[0,1]}\left|\langle
f(\cdot),{\mathbb{K}}(\cdot,s)\rangle_{\mathcal{H}}\right|\leq\|f\|_{\mathcal{H}}\sup_{s\in[0,1]}\|{\mathbb{K}}(\cdot,s)\|_{\mathcal{H}}\\\
=&\|f\|_{\mathcal{H}}\sup_{s\in[0,1]}\langle{\mathbb{K}}(\cdot,s),{\mathbb{K}}(\cdot,s)\rangle_{\mathcal{H}}^{1/2}=\|f\|_{\mathcal{H}}\sup_{s\in[0,1]}\sqrt{{\mathbb{K}}(s,s)}\leq\|f\|_{\mathcal{H}}.\end{split}$

Assumption 1 (b) is introduced for convenience of the later theoretical
analysis. This condition appears in the literature on high-dimensional
functional linear regression/auto-regression (Wang et al., 2020a, ; Wang et
al., 2020b, ) and also holds for some prevalent functional spaces of interest.
For example, Wang et al., 2020a ; Wang et al., 2020b  stated that the Sobolev
space

$W^{\alpha,2}=\left\\{f:f^{(r)},\text{the $r$-th derivative of $f$, is
absolutely continuous, $r=0,\ldots,\alpha$;
}f^{(\alpha)}\in\mathcal{L}^{2}([0,1])\right\\}$ (4)

satisfies (3) with some constant $C_{\mathcal{H}}$ that only depends on
$\alpha$ and showed in particular that $W^{1,2}$ satisfies (3) with
$C_{\mathcal{H}}=\sqrt{5}$.

Define the _effective dimension_ of $\mathcal{H}$ as
$p_{\mathcal{H}}:=\sum_{k=1}^{\infty}\mu_{k}^{2}/\mu_{1}^{2}.$ Intuitively
speaking, $p_{\mathcal{H}}$ captures the intrinsic dimension of the Hilbert
space $\mathcal{H}$ in the sense that the Gaussian processes over
$\mathcal{H}$ has the similar incoherence as the Gaussian random vectors in a
$p_{\mathcal{H}}$-dimensional Euclidean space (see the forthcoming Proposition
1 for details). As a concrete example, suppose $\mu_{k}=1$ for $k\leq d$ and
$\mu_{k}=0$ for $k>d$. Then, $p_{\mathcal{H}}=d$ and $\mathcal{H}$ is
equivalent to a $d$-dimensional Euclidean space.

###  2.2 Model

Let ${\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times[0,1]}$ be the
functional tensor of interest. Assume ${\mathbfcal{Y}}$ is approximately CP
rank-$r$:

${\mathbfcal{Y}}={\mathbfcal{X}}+{\mathbfcal{Z}},\qquad{\mathbfcal{X}}=\sum_{l=1}^{r}\lambda_{l}a_{l}\otimes
b_{l}\otimes\xi_{l}\in\mathbb{R}^{p_{1}\times p_{2}\times[0,1]}.$ (5)

Here, each tuple $(\lambda_{l},a_{l},b_{l},\xi_{l})$ corresponds to a singular
component of ${\mathbfcal{X}}$: $\lambda_{l}$ is the singular value,
$a_{l},b_{l}$ are the singular vectors, and $\xi_{l}$ is the singular
function. We also assume
$\|a_{l}\|_{2}=\|b_{l}\|_{2}=\|\xi_{l}\|_{\mathcal{L}_{2}}=1$ and
$\lambda_{l}>0$ for model identifiability. We denote
$\lambda_{max}=\max_{l\in[r]}\lambda_{l}$ and
$\lambda_{min}=\min_{l\in[r]}\lambda_{l}$. The functional tensor
${\mathbfcal{Z}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times[0,1]}$ is included to
model the unexplained remainder term. We also assume $\xi_{l}$ satisfies the
regularity condition, such that $\|\xi_{l}\|_{\mathcal{H}}\leq C_{\xi}$ for
some given RKHS $\mathcal{H}$ and some universal constant $C_{\xi}$.

Motivated by the previously discussed applications, we assume the functional
tensor is measured on a discrete grid of time points
$\\{s_{k}\\}_{k=1}^{n}\overset{i.i.d.}{\sim}\text{Unif}(0,1)$ and the
following tabular tensor
$\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times n}$ is
observed:

$\widetilde{\mathbfcal{Y}}_{ijk}:={\mathbfcal{Y}}_{ijs_{k}}+\varepsilon_{ijk},\quad
i\in[p_{1}],j\in[p_{2}],k\in[n].$ (6)

Here, $\varepsilon_{ijk}$’s are independent Gaussian random noises with mean
zero and variance $\tau^{2}$. Our goal is to perform functional tensor
singular value decomposition on ${\mathbfcal{Y}}$, that is to estimate
$\\{(a_{l},b_{l},\xi_{l})\\}_{l=1}^{r}$ based on $\widetilde{\mathbfcal{Y}}$
and $\\{s_{k}\\}_{k=1}^{n}$ from the model (5) and (6).

##  3 Methods for Functional Tensor SVD

Next, we discuss the estimation methods for FTSVD.

Power iterations. As a starting point, we approach the problem from an
optimization perspective:

$\begin{split}\min_{\lambda_{l},a_{l},b_{l},\xi_{l}}&\quad\sum_{i,j,k}\left(\widetilde{\mathbfcal{Y}}_{ijk}-\sum_{l=1}^{r}\lambda_{l}\cdot(a_{l})_{i}\cdot(b_{l})_{j}\cdot\xi_{l}(s_{k})\right)^{2}\\\
\text{subject to
}&\quad\|a_{l}\|_{2}=\|b_{l}\|_{2}=\|\xi_{l}\|_{\mathcal{L}^{2}}=1,\|\xi_{l}\|_{\mathcal{H}}\leq
C_{\xi}.\end{split}$ (7)

However, the computation of (7) is highly nontrivial and the exact solution is
computationally intractable in general since (a) the problem is highly non-
convex or even NP-hard due to the tensor product and the multilinear structure
(Hillar and Lim,, 2013); (b) (7) is essentially an infinite-dimensional
optimization problem due to the functional argument $\xi_{l}$. Therefore, we
introduce the following RKHS-based constrained power iteration method to
overcome these two difficulties. As a popular and powerful method for singular
value decomposition, power iteration has been successfully applied to tabular
tensor decomposition in the past two decades (Zhang and Golub,, 2001;
Anandkumar et al., 2014a, ; Anandkumar et al., 2014b, ; Sun et al.,, 2017).

We first discuss the power iteration for the one-component case (i.e., $r=1$),
while the multiple-component scenario is postponed later. Given estimates at
Step $t$: ($a^{(t)}$, $b^{(t)}$, $\xi^{(t)}$), power iteration updates
$a^{(t)}$, $b^{(t)}$, and $\xi^{(t)}$ alternatively. Specifically, recall
$\xi_{n}^{(t)}=(\xi^{(t)}(s_{1}),\ldots,\xi^{(t)}(s_{n}))$ is the
discretization of the function $\xi^{(t)}$. We update the singular vectors
$a^{(t)},b^{(t)}$ by projection-normalization:

$\begin{split}\tilde{a}^{(t+1)}&=\widetilde{\mathbfcal{Y}}\times_{2}b^{(t)}\times_{3}\xi_{n}^{(t)},\qquad
a^{(t+1)}=\tilde{a}^{(t+1)}/\|\tilde{a}^{(t+1)}\|_{2},\\\
\tilde{b}^{(t+1)}&=\widetilde{\mathbfcal{Y}}\times_{1}a^{(t)}\times_{3}\xi_{n}^{(t)},\qquad
b^{(t+1)}=\tilde{b}^{(t+1)}/\|\tilde{b}^{(t+1)}\|_{2}.\end{split}$ (8)

Since only a discrete subset of observations from the functional tensor are
accessible, we consider the following optimization to update the singular
function estimator:

$\begin{split}&\tilde{\xi}^{(t+1)}=\mathop{\rm
arg\min}_{\xi\in\mathcal{H}:\|\xi\|_{\mathcal{H}}\leq(C_{\xi}\lambda_{\max})}\leavevmode\nobreak\
\sum_{i,j,k}\left(\widetilde{\mathbfcal{Y}}_{ijk}-(a^{(t)})_{i}\cdot(b^{(t)})_{j}\cdot\xi(s_{k})\right)^{2},\quad\xi^{(t+1)}=\tilde{\xi}^{(t+1)}/\|\tilde{\xi}^{(t+1)}\|_{\mathcal{L}^{2}}.\end{split}$
(9)

Note that the first part of (9) is essentially a weighted mean functional
estimation problem, or a special case of functional linear regression (Yao et
al., 2005a, ; Cai and Yuan,, 2012). Since $\tilde{\xi}^{(t+1)}$ is essentially
a minimizer of a regularized empirical risk functional defined over the RKHS
$\mathcal{H}$, it admits a finite-dimensional closed-form solution by the
classic Representer Theorem (Kimeldorf and Wahba,, 1971) to be discussed in
Section 5.

Let $(\hat{a},\hat{b},\hat{\xi})$ be the final estimates after sufficient
iterations. The singular value is estimated via
$\hat{\lambda}=\widetilde{\mathbfcal{Y}}\times_{1}\hat{a}\times_{2}\hat{b}\times_{3}\hat{\xi}_{n}$.
The procedure of one-component power iteration is summarized to Algorithm 1.

Algorithm 1 Regularized Power Iteration

0: Tensor $\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times
n}$, initialization $(a^{(0)},b^{(0)})$, iteration number $T$

Calculate $\xi_{0}$ via:  $\tilde{\xi}^{(0)}=\mathop{\rm
arg\min}_{\|\xi\|_{\mathcal{H}}\leq(C_{\xi}\lambda_{\max})}\sum_{i,j,k}\left(\widetilde{\mathbfcal{Y}}_{ijk}-(a^{(0)})_{i}\cdot(b^{(0)})_{j}\cdot\xi(s_{k})\right)^{2},\qquad\xi^{(0)}=\tilde{\xi}^{(0)}/\|\tilde{\xi}^{(0)}\|_{\mathcal{L}^{2}}.$

for all $t=0,T-1$ do

Calculate $a^{(t+1)},b^{(t+1)}$ via (8) and calculate $\xi^{(t+1)}$ via (9).

end for

Calculate
$\lambda^{(T)}=\widetilde{\mathbfcal{Y}}\times_{1}a^{(T)}\times_{2}b^{(T)}\times_{3}\xi^{(T)}_{n}$.

return
$(\hat{\lambda},\hat{a},\hat{b},\hat{\xi})=(\hat{\lambda},a^{(T)},b^{(T)},\xi^{(T)})$

Initialization and Overall Algorithm. Next, we discuss initialization scheme
for $(a^{(0)},b^{(0)})$, which is required for the implementation of Algorithm
1. If ${\mathbfcal{X}}$ is a rank-1 tensor, the matricizations of the signal
tensor ${\mathbfcal{X}}$ admit the following rank-1 decompositions:

$\mathcal{M}_{1}({\mathbfcal{X}})=\lambda_{1}a_{1}\left(b_{1}\otimes(\xi_{1})_{n}\right)^{\top},\qquad\mathcal{M}_{2}({\mathbfcal{X}})=\lambda_{1}b_{1}\left((\xi_{1})_{n}\otimes
a_{1}\right)^{\top},$

i.e., $a_{1}$ and $b_{1}$ are the the left singular vectors of
$\mathcal{M}_{1}({\mathbfcal{X}})$ and $\mathcal{M}_{2}({\mathbfcal{X}})$
respectively. Given $\widetilde{\mathbfcal{Y}}$ being a noisy substitute of
${\mathbfcal{X}}$, it is natural to initialize $a_{1}$ and $b_{1}$ by the
first singular vectors of $\mathcal{M}_{1}(\widetilde{\mathbfcal{Y}})$ and
$\mathcal{M}_{2}(\widetilde{\mathbfcal{Y}})$, respectively:

$a^{(0)}={\rm SVD}_{1}(\mathcal{M}_{1}(\widetilde{\mathbfcal{Y}})),\quad
b^{(0)}={\rm SVD}_{1}(\mathcal{M}_{2}(\widetilde{\mathbfcal{Y}})).$ (10)

Our estimation procedure for the rank one scenario is then completed by
combining Algorithm 1 with Eq. (10).

When the rank $r>1$, there are multiple singular components to be estimated
and the power iteration will no longer be the exact alternating minimization
scheme. Nevertheless, when all the singular components
$\\{(a_{l},b_{l},\xi_{l})\\}_{l=1}^{r}$ satisfy the incoherent condition and
$(a^{(0)},b^{(0)})$ satisfies some initialization conditions (formal
definition and discussions are postponed to Section 4.1), Algorithm 1 still
yields estimators with guaranteed local convergence. Therefore, we propose to
perform spectral initialization scheme sequentially: we first let $a^{(0)}$
and $b^{(0)}$ be initialized via (10) and let
$(\hat{\lambda}_{1},\hat{a}_{1},\hat{b}_{1},\hat{\xi}_{1})$ be the first
estimated singular component by power iteration (Algorithm 1). Then we
subtract
$\hat{\lambda}_{1}\otimes\hat{a}_{1}\otimes\hat{b}_{1}\otimes(\hat{\xi}_{1})_{n}$
from the original observation $\widetilde{\mathbfcal{Y}}$. Note that this is
essentially equivalent to projecting $\widetilde{\mathbfcal{Y}}$ onto the low-
dimensional tensor space spanned by
$\hat{a}_{1}\otimes\hat{b}_{1}\otimes(\hat{\xi}_{1})_{n}$ and replacing the
data with the projection residue. We perform this procedure for $r$ times to
obtain $r$ estimated singular components. The pseudocode of this procedure is
summarized to Algorithm 2.

Algorithm 2 Functional Tensor SVD with Sequential Spectral Initialization

0: Tensor $\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times p_{2}\times
n}$, rank $r$

for $l=1,\ldots,r$ do

Calculate $a^{(0)}={\rm
SVD}_{1}\left(\mathcal{M}_{1}(\widetilde{\mathbfcal{Y}})\right)$,
$b^{(0)}={\rm
SVD}_{1}\left(\mathcal{M}_{2}(\widetilde{\mathbfcal{Y}})\right).$

Calculate $(\hat{\lambda}_{l},\hat{a}_{l},\hat{b}_{l},\hat{\xi}_{l})$ by
applying Algorithm 1 on $(a^{(0)},b^{(0)})$.

Update
$\widetilde{\mathbfcal{Y}}=\widetilde{\mathbfcal{Y}}-\hat{\lambda}_{l}\hat{a}_{l}\otimes\hat{b}_{l}\otimes(\hat{\xi}_{l})_{n}$.

end for

return
$\\{(\hat{\lambda}_{l},\hat{a}_{l},\hat{b}_{l},\hat{\xi}_{l})\\}_{l=1}^{r}$

Rank selection. The proposed algorithm requires the input of rank $r$ in
practice. We introduce the following BIC-type criteria to select the value of
$r$: $\hat{r}=\mathop{\rm arg\min}_{r}{\rm BIC}(r)$, where

$\text{BIC}(r):=2\log\left\\{\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}\left(\widetilde{\mathbfcal{Y}}_{ijk}-\sum_{l=1}^{r}\hat{\lambda}_{l}\hat{a}_{l}\otimes\hat{b}_{l}\otimes(\hat{\xi}_{l})_{n}\right)^{2}\right\\}+\frac{\log(p_{1}p_{2}n)}{p_{1}p_{2}n}\cdot(p_{1}+p_{2}+p_{\mathcal{H}})r.$
(11)

Here, $p_{\mathcal{H}}:=\sum_{k=1}^{\infty}\mu_{k}^{2}/\mu_{1}^{2}$ is the
effective dimension of the RKHS and $(p_{1}+p_{2}+p_{\mathcal{H}})r$ is the
effective dimensions of the rank-$r$ tensor model.

##  4 Theory

In this section, we present the theoretical results for the proposed
estimation algorithms. We particularly aim to study the statistical error
bound of $\hat{a},\hat{b},\hat{\xi}$. Since the singular vectors and functions
are identifiable up to sign flips, we focus on the following _sine values_ of
the pairs of vectors/functions:

$\begin{split}{\rm
dist}(u,v)&=\sqrt{1-\left(\frac{u^{\top}v}{\|u\|_{2}\|v\|_{2}}\right)^{2}},\qquad\forall
u,v\in{\mathbb{R}}^{p};\\\ {\rm
dist}(f,g)&=\sqrt{1-\left(\frac{\int_{0}^{1}f(t)g(t)dt}{\|f\|_{\mathcal{L}^{2}}\|g\|_{\mathcal{L}^{2}}}\right)^{2}},\qquad\forall
f,g\in\mathcal{L}^{2}([0,1]).\end{split}$ (12)

We introduce the following quantity associated with the RKHS $\mathcal{H}$:

$\zeta_{n}:=\inf\left\\{\zeta\geq\sqrt{\frac{\log
n}{n}}:Q_{n}(\delta)\leq\zeta\delta+\zeta^{2},\quad\forall\delta\in(0,1]\right\\},$
(13)

where

$Q_{n}(\delta):=\left(\sum_{k=1}^{\infty}\min\\{\delta^{2},\mu_{k}\\}/n\right)^{1/2},$
(14)

and $\mu_{1}\geq\mu_{2}\geq\ldots\geq 0$ are the eigenvalues of the
reproducing kernel $\mathbb{K}$ of $\mathcal{H}$. Essentially, $\zeta_{n}$
quantifies the information loss from only observing measurements on a discrete
grid, rather than the whole function, of $f$ in the RKHS $\mathcal{H}$. This
quantity is commonly used in the literature on functional data analysis (for
example, see (Koltchinskii and Yuan,, 2010, Eqn (8)); also see (Raskutti et
al.,, 2012, Eqn (13)) and (Wang et al., 2020b, , Eqn (16)) for almost
equivalent variants). As the grid density $n$ increases,
$\\{f(s_{k})\\}_{k=1}^{n}$ reveals more information on $f$ and the value of
$\zeta_{n}$ decreases. For some specific RKHS $\mathcal{H}$, the explicit rate
of $\zeta_{n}$ can be obtained. As an example, When $\mathcal{H}$ is the
Hilbert-Sobolev space $W^{\alpha,2}([0,1])$ defined in (4), we have
$\zeta_{n}\lesssim n^{-\alpha/(2\alpha+1)}$ (see Proposition 4). More
technical results on $Q_{n}(\delta)$ and $\zeta_{n}$ will be discussed in
Appendix B.

We also introduce the following quantity to measure the scale of the
unexplained remainder term ${\mathbfcal{Z}}$:

$\begin{split}&\mathcal{E}:=\sup_{s\in[0,1]}\left\|{\mathbfcal{Z}}_{\cdot\cdot
s}\right\|.\end{split}$ (15)

Here, $\|\cdot\|$ is the largest singular value of the matrix. $\mathcal{E}$
can be interpreted as an extension of the _tensor spectral norm_ from the
tabular tensors, which appears as the prominent term in the estimation error
bound in various statistical tensor models (Anandkumar et al., 2014b, ; Sun et
al.,, 2017; Han et al.,, 2021). In particular, the forthcoming Proposition 2
explicitly characterizes $\mathcal{E}$ for Gaussian processes
${\mathbfcal{Z}}_{ij\cdot}$ via the dimensions of the tabular modes
$p_{1},p_{2}$.

###  4.1 Assumptions

In this section, we present the technical assumptions that are used to
establish the theoretical guarantees for the proposed algorithm. We first
introduce the following incoherence conditions on the true singular components
$(a_{l},b_{l},\xi_{l})$ of the signal tensor ${\mathbfcal{X}}$.

######  Assumption 2 (Incoherence).

Consider model (5). Assume

$\mu:=\max\left\\{\max_{i\neq j}\left\\{|\langle
a_{i},a_{j}\rangle|\right\\},\leavevmode\nobreak\ \max_{i\neq
j}\left\\{|\langle b_{i},b_{j}\rangle|\right\\},\leavevmode\nobreak\
\max_{i\neq
j}\left\\{|\langle\xi_{i},\xi_{j}\rangle_{\mathcal{L}^{2}}|\right\\}\right\\}\leq
c/(\kappa^{2}r^{2}),$

where $\kappa:=\lambda_{max}/\lambda_{min}$ is the condition number.

Assumption 2 suggests the singular components of ${\mathbfcal{X}}$ are
pairwise incoherent. Such the condition is widely assumed in the literature on
tabular CP low-rank tensor decomposition (Anandkumar et al., 2014b, ; Sun et
al.,, 2017; Cai et al.,, 2019; Hao et al.,, 2020). By (Anandkumar et al.,
2014b, , Lemma 2), one can show that $\max_{i\neq j}\\{|\langle
a_{i},a_{j}\rangle|\\}\asymp 1/\sqrt{p_{1}}$ and $\max_{i\neq j}\\{|\langle
b_{i},b_{j}\rangle|\\}\asymp 1/\sqrt{p_{2}}$ with high probability if
$\\{a_{l}\\}_{l=1}^{r}$ and $\\{b_{l}\\}_{l=1}^{r}$ are i.i.d. uniformly
sampled from $\mathbb{S}^{p_{1}-1}$ and $\mathbb{S}^{p_{2}-1}$, respectively.
We can also establish the incoherence for singular functions under certain
sampling distributions descried as follows.

######  Proposition 1.

Suppose $r\geq 2$. Let $f_{i}\sim
GP(0,{\mathbb{G}}(\cdot,\cdot)),i=1,\ldots,r$ be i.i.d. mean-zero Gaussian
processes with covariance function
${\mathbb{G}}(s,t)=\sum_{k=1}^{\infty}\mu_{k}^{2}\phi_{k}(s)\phi_{k}(t)$.
Assume $\log r<c\left(\sum_{k=1}^{\infty}\mu_{k}^{4}/\mu_{1}^{4}\right)$ and
normalize $\xi_{i}=f_{i}/\|f_{i}\|_{\mathcal{L}_{2}}$. Recall the effective
dimension $p_{\mathcal{H}}:=\sum_{k=1}^{\infty}\mu_{k}^{2}/\mu_{1}^{2}$. Then,
with probability at least $1-Cr^{-8}$,

$\begin{split}\xi_{i}\in\mathcal{H},\quad|\langle\xi_{i},\xi_{j}\rangle|\leq
C\log r/\sqrt{p_{\mathcal{H}}},\quad\forall 1\leq i\neq j\leq r.\end{split}$

Our next assumption is on the grid density and signal-to-noise ratio (SNR).

######  Assumption 3 (Grid density and SNR).

The grid density $n$ and the least singular value $\lambda_{min}$ satisfy

$\frac{n}{\log^{2}n}\geq C\kappa(p_{1}+p_{2}),\qquad\zeta_{n}\leq
c/\kappa^{2}r,\qquad\lambda_{min}\geq C\kappa
r\left\\{\mathcal{E}+\tau\left(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n}\right)\right\\}$
(16)

for some sufficiently small constant $c>0$ and large constant $C>0$.

Recall that $\zeta_{n}$ characterizes the difficulty of estimating a function
in $\mathcal{H}$ under discretization and Gaussian observational errors. In
most interesting scenarios (which will be discussed later), these quantities
have polynomial decay with respect to $n$, and the first two conditions in
(16) is mild that allows $p,r$ to grow with $n$. The third condition in (16)
is similar to the widely-assumed SNR condition in the literature on tabular
tensor decomposition (Anandkumar et al., 2014b, ; Sun et al.,, 2017; Zhang and
Xia,, 2018; Han et al.,, 2021), where the ratios of noise (e.g.,
$\mathcal{E}_{2},\tau$) and the least singular value ($\lambda_{min}$) are
required to be sufficiently large. This guarantees that each singular
component has strong-enough signal and can be estimated consistently.

######  Assumption 4 (Initialization).

There exists some $l\in[r]$ and constant $c_{0}>0$ such that the
initialization error satisfies

$\max\left\\{{\rm dist}(a^{(0)},a_{l}),{\rm dist}(b^{(0)},b_{l})\right\\}\leq
c_{0}/(\kappa r).$

Assumption 4 requires that the initial tabular mode estimations are reasonably
close to the ground truth. Such the condition is widely used in the
theoretical analysis for power iteration (Anandkumar et al., 2014b, ; Sun et
al.,, 2017; Cai et al.,, 2019; Han et al.,, 2021). In particular, when
$\mathcal{X}$ is a well-conditioned tensor with a constant rank (i.e.,
$r,\kappa=O(1)$), this assumption means the initialization error is only
smaller than a constant. Later, we will show the initial estimator yielded by
the spectral initialization (Eqn. (10)) satisfies this assumption.

###  4.2 Local Convergence

The following theorem gives an estimation error bound for the proposed RKHS-
based constrained power iteration (Algorithm 1).

######  Theorem 1.

Suppose Assumptions 2, 3, and 4 hold. Let $a^{(t)},b^{(t)},\xi^{(t)}$ be
estimated singular vectors and function at step $t$ of Algorithm 1. Recall
$\mathcal{E}$ is defined in (15), $\zeta_{n}$ is defined in (13), $\mu$ is
defined in Assumption (2). Then with probability at least $1-2n^{-9}-C\log
n\cdot e^{-c(p_{1}\wedge p_{2})}$, we have

$\begin{split}&\max\\{{\rm dist}(a^{(t)},a_{l}),{\rm dist}(b^{(t)},b_{l}),{\rm
dist}(\xi^{(t)},\xi_{l})\\}\\\ &\qquad\leq
2^{-t}+C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau\left(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n}\right)}{\lambda_{min}}+\zeta_{n}+(r-1)\mu\right\\},\quad\forall
t\geq 0.\end{split}$ (17)

Specifically if $r=1$, the estimation error of $a^{(t)},b^{(t)}$ enjoys the
following better rate:

$\max\\{{\rm dist}(a^{(t)},a_{l}),{\rm dist}(b^{(t)},b_{l})\\}\leq
C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n})}{\lambda_{min}}\right\\},\quad\forall
t\geq 1.$

######  Remark 1 (Interpretation of Estimation Error Bound).

Theorem 1 shows if the number of iterations $T>\Omega(\log(1/\zeta_{n}))$, the
estimation error of $\xi_{s}$ can be upper bounded by the sum of the following
four terms: $\mathcal{E}/\lambda_{min}$,
$\tau\left(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n}\right)/\lambda_{min}$, $\zeta_{n}$
and $(r-1)\mu.$ Here, $\mathcal{E}/\lambda_{min}$ is induced by the
unexplained remainder term ${\mathbfcal{Z}}$, which can be interpreted as the
“spectral-norm-over-least-singular-value” term that widely appears in the
literature on tabular tensor decomposition (Anandkumar et al., 2014b, ; Han et
al.,, 2021). The second term
$\tau\left(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n}\right)/\lambda_{min}$ accounts for
the observational noise $\varepsilon_{ijk}$. On the other hand, the classic
analysis of tabular tensor decomposition (Zhang and Xia,, 2018, Theorem 1),
(Han et al.,, 2020, Proposition 1) only yields a much weaker rate of the
tabular estimation error,
$\tau/\lambda_{min}\cdot\left((p_{1}p_{2}/n)^{1/4}+\sqrt{(p_{1}+p_{2})/n}\right)$,
which can be much larger than the bound in Theorem 1 in the high-dimensional
region where $p_{1}p_{2}$ is comparable to $n$. This phenomenon also appears
in our simulation studies in Section 5.2. The third term corresponds to the
discretization error on estimating the function $f$ from a discrete grid,
which is unavoidable even in the noiseless orthogonal-decomposable case (i.e.,
$\mathcal{E}_{2}=\mathcal{E}_{\mathcal{H}}=\tau=\mu=0$). The last term
$(r-1)\mu$ corresponds to the coherence of the signal tensor
${\mathbfcal{X}}$: when different singular components are not exactly
orthogonal, the true singular vectors/functions $(a_{l},b_{l},\xi_{l})$ is not
a stationery point even in the noiseless cases. Nevertheless, this term often
has negligible rate compared to the other three terms in many scenarios (for
example, in the scenario of Proposition 1, $\mu\lesssim
1/\sqrt{\min\\{p_{1},p_{2},p_{\mathcal{H}}\\}}$; see the discussions therein).

It is also noteworthy that when the tensor ${\mathbfcal{X}}$ is rank-$1$, the
estimation error for singular functions does not contain the discretization
term $\zeta_{n}$ while this is not true for $r>1$. The reason is that when
there is only one singular component, the “inaccurate” estimation $\hat{\xi}$
will not introduce additional bias from other singular components when
projection is performed in power iteration. On the other hand, the estimation
error of $\hat{\xi}$ will be accumulated to the tabular modes when there are
more than one singular components in the signal tensor ${\mathbfcal{X}}$.

Note that the result of Theorem 1 is conditioning on fixed values of the
functional remainder term $\mathcal{E}$. When ${\mathbfcal{Z}}$ is random, we
can develop finer upper bounds on $\mathcal{E}$ as well as the estimation
error. For example, if ${\mathbfcal{Z}}$ is a collection of i.i.d. Gaussian
processes, we have the following result.

######  Assumption 5.

Suppose ${\mathbfcal{Z}}_{ij\cdot}$ are i.i.d. mean-zero Gaussian processes
with almost sure continuous path and

${\mathbb{E}}\sup_{s}|{\mathbfcal{Z}}_{ijs}|\leq m,\qquad\sup_{s}{\rm
Var}({\mathbfcal{Z}}_{ijs})\leq\sigma^{2}.$ (18)

######  Proposition 2.

Suppose Assumption 5 holds. Then with probability at least
$1-e^{-c(p_{1}+p_{2})}$,

$\begin{split}\mathcal{E}\leq m+\sigma\sqrt{p_{1}+p_{2}}.\end{split}$

To illustrate our theoretical results, we specifically consider two RKHSs of
major interests next. First, if the RKHS kernel ${\mathbb{K}}$ has $d$ non-
zero eigenvalues, it follows that $\zeta_{n}\lesssim\sqrt{d/n}$ (Proposition
4) and we can obtain the following estimation error upper bound.

######  Corollary 1 (Local Convergence: Finite-dimensional RKHS).

Suppose the reproducing kernel ${\mathbb{K}}$ of $\mathcal{H}$ has finitely
many, say $d$, non-zero eigenvalues. Suppose $\log n\leq d\leq n$, Assumptions
2, 3, 4, and 5 hold. Then, with probability at least $1-Cn^{-9}-C\log n\cdot
e^{-c(p_{1}\wedge p_{2})}$,

$\begin{split}&\max\\{{\rm dist}(\hat{a},a_{l}),{\rm dist}(\hat{b},b_{l}),{\rm
dist}(\hat{\xi},\xi_{l})\\}\\\
&\qquad\lesssim\sqrt{\frac{d}{n}}+\frac{m+\sigma\sqrt{p_{1}+p_{2}}}{\lambda_{min}}+\sqrt{\frac{p_{1}+p_{2}+d}{n}}\cdot\frac{\tau}{\lambda_{min}}+(r-1)\mu.\end{split}$

Another interesting example is the Sobolev-Hilbert space $W^{\alpha,2}$
defined in (4). If $\mathcal{H}=W^{\alpha,2}$ with constant $\alpha>1/2$, then
$\mathcal{H}$ is an RKHS with eigenvalues $\mu_{k}\asymp k^{-2\alpha}$ (see,
e.g., Micchelli and Wahba, (1981)) and $\zeta_{n}\lesssim
n^{-\alpha/(2\alpha+1)}$ (Proposition 4). The following estimation error bound
holds accordingly:

######  Corollary 2 (Local Convergence: Sobolev Space).

Suppose $\mathcal{H}=W^{\alpha,2}$ with $\alpha>1/2$ and Assumptions 2, 3, 4
and 5 hold. Then, with probability at least $1-Cn^{-9}-C\log n\cdot
e^{-c(p_{1}\wedge p_{2})}$,

$\begin{split}&\max\\{{\rm dist}(\hat{a},a_{l}),{\rm dist}(\hat{b},b_{l}),{\rm
dist}(\hat{\xi},\xi_{l})\\}\\\ &\qquad\lesssim
n^{-\frac{\alpha}{2\alpha+1}}+\frac{m+\sigma\sqrt{p_{1}+p_{2}}}{\lambda_{min}}+\left(n^{-\frac{\alpha}{2\alpha+1}}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\frac{\tau}{\lambda_{min}}+(r-1)\mu.\end{split}$

######  Remark 2 (Consistency Conditions for FTSVD in Sobolev Space and Its
Phase Transition).

When $m\ll\sigma\sqrt{p_{1}+p_{2}}$, ${\mathbfcal{X}}$ is orthogonally
decomposable (i.e., $\mu=0$) and let $n\rightarrow\infty$, Corollary 2 implies
that $\hat{a},\hat{b},\hat{\xi}$ are consistent estimators if

$\lambda_{min}/\sigma\gg\sqrt{p_{1}\vee p_{2}},\qquad\lambda_{min}/\tau\gg
n^{-\frac{\alpha}{2\alpha+1}}\vee\sqrt{(p_{1}\vee p_{2})/n}.$

Here, $\lambda_{min}/\sigma\gg\sqrt{p_{1}\vee p_{2}}$ is required for
consistent estimation of the singular vectors $a_{l},b_{l}$ (also see Cai et
al., (2019)); while $\lambda_{min}/\tau\gg
n^{-\alpha/(2\alpha+1)}\vee\sqrt{(p_{1}\vee p_{2})/n}$ is required for
consistent estimation/approximation of the singular functions $\xi_{l}$ and it
yields an interesting phase transition: when $p_{1}\vee p_{2}\ll
n^{1/(2\alpha+1)}$, the bound is dominated by the non-parametric rate
$n^{-\alpha/(2\alpha+1)}$; when $p_{1}\vee p_{2}\gg n^{1/(2\alpha+1)}$, the
bound is dominated by the parametric rate $\sqrt{(p_{1}\vee p_{2})/n}$. We
note that the phase transition between parametric and non-parametric rates
commonly appear in the study of high-dimensional functional regression with
sparsity, e.g., Wang et al., 2020b . To the best of our knowledge, we are the
first to establish such the phenomenon in the low-rank-based functional data
analyses.

###  4.3 Initialization

Next, we provide the theoretical guarantee for the initialization procedure
with rank $r=1$ described in (10). The theoretical analysis for initialization
of the rank $r>1$ case is technical challenging, which we leave as future
research. To simplify the notation, we omit the subscripts and denote
$\lambda=\lambda_{1}=\lambda_{min}$, $a=a_{1}$, etc.

######  Theorem 2.

Suppose $r=1$, $\zeta_{n}\leq c_{0}$ and Assumption 5 is satisfied. Recall
$\sigma$ and $\tau$ are the noise level of ${\mathbfcal{Z}}$ (defined in (18))
and $\varepsilon_{ijk}$ (defined in (6)), respectively. Suppose $n\geq C_{0}$
for some sufficiently large constant $C_{0}$. Then, for any $\delta\in(0,1)$,
there exits some constants $C_{\delta}$ and $C_{\delta}^{\prime}$ which depend
on $\delta$ such that as long as

$\begin{split}\frac{\lambda}{\sigma}&\geq C_{\delta}\left(\sqrt{p_{1}\vee
p_{2}}\right),\\\ \frac{\lambda}{\tau}&\geq
C_{\delta}^{\prime}\left(\sqrt{\frac{p_{1}\vee
p_{2}}{n}}+\left(\frac{p_{1}p_{2}}{n}\right)^{1/4}\right),\end{split}$ (19)

we have with probability at least $1-n^{-8}-e^{-\sqrt{p_{1}\wedge p_{2}}}$
that

$\begin{split}\max\left\\{{\rm dist}(a^{(0)},a),{\rm
dist}(b^{(0},b)\right\\}\leq\delta.\end{split}$

Moreover, if Assumptions 2, 3 and (19) hold, we have the same conclusion as
Theorem 1.

Theorem 2 implies that the initialization condition (Assumption 4) for local
convergence (Theorem 1) holds as long as the signal-to-noise ratio (SNR)
condition (19) holds.

##  5 Numeric Experiments

In this section, we investigate the performance of the proposed Functional
Tensor SVD estimators, which we refer as FTSVD for short later, on both
synthetic and real datasets. We use the rescaled Bernoulli polynomial as the
reproducing kernel throughout the section:

${\mathbb{K}}(x,y)=1+k_{1}(x)k_{1}(y)+k_{2}(x)k_{2}(y)-k_{4}(|x-y|),$

where $k_{1}(x)=x-.5$, $k_{2}(x)=(k_{1}^{2}(x)-1/12)/2$, and
$k_{4}(x)=(k_{1}^{4}(x)-k_{1}^{2}(x)/2+7/240)/24$ for any $x\in[0,1]$. Note
that ${\mathbb{K}}$ is the reproducing kernel for the Hilbert space $W^{2,2}$
(Gu,, 2013, Chapter 2.3.3).

###  5.1 Computation of (9)

Before proceeding further, we first discuss the computation strategy for (9),
which is a crucial step in updating the singular functions. By equivalence of
constrained and regularized optimization (Friedman et al.,, 2001), we
reformulate (9) to

$\tilde{\xi}^{(t+1)}=\mathop{\rm
arg\min}_{\xi\in\mathcal{L}^{2}([0,1])}\sum_{i,j,k}\left(\widetilde{\mathbfcal{Y}}_{ijk}-(a^{(t)})_{i}\cdot(b^{(t)})_{j}\cdot\xi(s_{k})\right)^{2}+C_{\lambda}\left\|\xi\right\|_{\mathcal{H}}^{2}$
(20)

for some regularization parameter $C_{\lambda}>0$. By Representer Theorem
(Kimeldorf and Wahba,, 1971; Wahba,, 1990), all minimizers
$\tilde{\xi}^{(t+1)}$ of (20) can be written as a finite linear combination of
kernel products, i.e., $\tilde{\xi}^{(t+1)}$ admits the following
representation:

$\tilde{\xi}^{(t+1)}(\cdot)=\sum_{k=1}^{n}\beta_{k}{\mathbb{K}}(\cdot,s_{k})=k(s)^{\top}\beta,$
(21)

where
$k(s)=\left({\mathbb{K}}(s,s_{1}),\ldots,{\mathbb{K}}(s,s_{n})\right)^{\top}\in{\mathbb{R}}^{n}$
and some
$\beta=\left(\beta_{1},\ldots,\beta_{n}\right)^{\top}\in{\mathbb{R}}^{n}$. We
further denote
$\boldsymbol{K}=\left[k(s_{1}),\cdots,k(s_{n})\right]\in{\mathbb{R}}^{n\times
n}$ as the discrete kernel matrix. Then the squared loss term in (20) can be
written as

$\sum_{i,j,k}\left(\widetilde{\mathbfcal{Y}}_{ijk}-(a^{(t)})_{i}\cdot(b^{(t)})_{j}\cdot\xi(s_{k})\right)^{2}=\left\|\mathcal{M}_{3}(\widetilde{\mathbfcal{Y}})-\boldsymbol{K}\beta\left(a^{(t)}\otimes
b^{(t)}\right)^{\top}\right\|_{\rm F}^{2}.$

Meanwhile,

$\left\|\xi\right\|_{\mathcal{H}}^{2}=\sum_{k_{1},k_{2}\in[n]}\beta_{k_{1}}\beta_{k_{2}}\left\langle{\mathbb{K}}(\cdot,s_{k_{1}}),{\mathbb{K}}(\cdot,s_{k_{2}})\right\rangle_{\mathcal{H}}=\sum_{k_{1},k_{2}\in[n]}\beta_{k_{1}}\beta_{k_{2}}{\mathbb{K}}(s_{k_{1}},s_{k_{2}})=\beta^{\top}\boldsymbol{K}\beta.$

Then, (20) has the following closed-form solution:

$\begin{split}\beta^{(t+1)}&=\mathop{\rm
arg\min}_{\beta}\left\|\mathcal{M}_{3}(\widetilde{\mathbfcal{Y}})-\boldsymbol{K}\beta\left(a^{(t)}\otimes
b^{(t)}\right)^{\top}\right\|_{\rm
F}^{2}+C_{\lambda}\beta^{\top}\boldsymbol{K}\beta\\\
&\overset{(a)}{=}\left(\boldsymbol{K}+C_{\lambda}\boldsymbol{I}\right)^{-1}\mathcal{M}_{3}(\widetilde{\mathbfcal{Y}})\left(a^{(t)}\otimes
b^{(t)}\right).\end{split}$ (22)

Here (a) is due to $\|a^{(t)}\otimes b^{(t)}\|_{2}=1$. Therefore, (21) and
(22) together yield a closed-form solution to the update of the singular
function estimator (9) in power iteration.

Note that $C_{\lambda}$ is a tuning parameter that controls the smoothness
level of the fitted function
$\tilde{\xi}^{(t+1)}(\cdot)=k(\cdot)^{\top}\beta^{(t+1)}$. In all the upcoming
experiments, we use a fixed value $C_{\lambda}=10^{-5}$. Unlike tuning
parameters for sparsity or low-rankness regularization, a suitable value of
$C_{\lambda}$ is hard to determine in theory as it is connected to the unknown
signal-to-noise ratio and intrinsic smoothness of the function. We thus leave
it as future research.

###  5.2 Simulation Studies

We simulate the $p_{1}\times p_{2}\times[0,1]$-dimensional rank-$r$ functional
tensor datasets as follows. For each $l\in[r]$, we sample $a_{l},b_{l}$
uniformly from the unit spheres $\mathbb{S}^{p_{1}-1},\mathbb{S}^{p_{2}-1}$
and generate $\xi_{l}$ from orthonormal basis functions
$\\{u_{i}(s)\\}_{i=1}^{10}\subset\mathcal{L}^{2}([0,1])$. Following Yuan and
Cai, (2010); Wang et al., 2020b , we set $u_{1}(s)=1$ and
$u_{i}(s)=\sqrt{2}\cos\left((i-1)\pi s\right)$ for $i=2,\ldots,10$. We
generate $x_{li}\sim\text{Unif}[-1/i,1/i]$ independently,
$\xi^{\prime}_{l}(\cdot)=\sum_{i=1}^{10}x_{li}u_{i}(\cdot)$, singular
functions $\xi_{l}=\xi^{\prime}_{l}/\|\xi^{\prime}_{l}\|_{\mathcal{L}^{2}}$,
and the signal tensor ${\mathbfcal{X}}=\sum_{l=1}^{r}\lambda_{l}a_{l}\otimes
b_{l}\otimes\xi_{l}\in{\mathbb{R}}^{p_{1}\times p_{2}\times[0,1]}$, where
$\lambda_{l}:=\lambda_{min}\cdot(r-l+1)$ and $\lambda_{min}>0$ is some pre-
specified least singular value. The remainder functions
${\mathbfcal{Z}}_{ij\cdot}$ are drawn independently in the same way as
$\xi_{l}$. We calculate
${\mathbfcal{Y}}={\mathbfcal{X}}+\sigma\cdot{\mathbfcal{Z}}$. Then, $\sigma$
essentially controls the amplitude of remainder functions. Finally, we
generate a discrete grid $\\{s_{k}\\}_{k=1}^{n}\sim\text{Unif}(0,1)$ and
obtain the observation $\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{p_{1}\times
p_{2}\times n}$ such that
$\widetilde{\mathbfcal{Y}}_{ijk}={\mathbfcal{Y}}_{ijs_{k}}+\varepsilon_{ijk}$,
where $\varepsilon_{ijk}\overset{i.i.d.}{\sim}N(0,\tau^{2})$. For each
simulation setting, the estimation errors are calculated by (12) and the
averages of 100-time repeated experiments are reported.

We first study the performance of FTSVD on rank-1 models. We set
$\lambda_{min}=20$, $\tau=1$, and $\sigma=0$. Note that this setting does not
include any functional remainder term ${\mathbfcal{Z}}$ in the observation and
can be taken as a particular tabular tensor CP-decomposition model. Therefore,
we compare our method with the classic method for rank-1 tabular tensor
decomposition (which we refer to as CP later). Since CP only yields the
estimation on the discrete grid for functional mode, we obtain the whole
function estimator by interpolation with minimal $\mathcal{H}$-norm. The
comparison between CP and the proposed FTSVD for dimension
$p=p_{1}=p_{2}\in\\{20,30,50\\}$ and various grid density $n$ is presented in
Figure 1. We can see the proposed algorithm has more accurate estimations than
the classic CP-decomposition method in all scenarios on both functional and
tabular modes. In particular, on the functional mode, our method is
significantly better than CP as we utilize the functional smoothness for
estimation at each step of power iteration, which further helps the estimation
of the other two singular vectors.

Figure 1: FTSVD versus CP under different dimension $p=p_{1}=p_{2}$ and grid
density $n$ for rank-$1$ functional tensor SVD. Upper and lower panels plot
the estimation errors for singular functions and vectors, respectively.

Next, we compare the performance of FTSVD with the existing FDA algorithms
when the functional remainder term ${\mathbfcal{Z}}$ exists, i.e., $\sigma>0$.
To this end, we unfold the original observational tensor
$\widetilde{\mathbfcal{Y}}$ to a matrix
$\widetilde{\boldsymbol{Y}}=\mathcal{M}_{3}(\widetilde{\mathbfcal{Y}})\in{\mathbb{R}}^{n\times(p_{1}p_{2})}$.
Then each of the $(p_{1}p_{2})$ columns of $\widetilde{\boldsymbol{Y}}$ can be
seen as a functional sample evaluated at a $n$-grid, which matches the
formulation of the functional principal component analysis (FPCA) in Yao et
al., 2005a . Here we only compare the estimation accuracy for the singular
function as the FDA-based method does not directly provide the singular vector
estimates. We take $\tau=0$ this time for simplicity, $p=20$ and $50$, and
vary the amplitude of the functional remainder $\sigma$. The grid density is
fixed to be $n=50$ and the least singular value $\lambda_{min}=20$. The result
is presented in Figure 2. We note that FPCA has uniformly higher estimation
errors than the other two tensor-based methods, as FPCA is designed for i.i.d.
samples but may not be suitable for the heterogeneous data considered in this
paper. On the other hand, one can see that CP and RKHS has almost the same
performance. This is because in this particular setting, the statistical rate
of FTSVD is dominated by $\mathcal{E}_{2}/\lambda_{min}$ as suggested by
Theorem 1, which is nearly the same as the one for CP (Anandkumar et al.,
2014b, ). However, it should be noted that in the general settings where both
functional remainder term and observational noises exist, CP may not be as
accurate as our method according to the results of the first simulation
setting.

Figure 2: Estimation error of singular functions under different dimension $p$
and amplitude of remainder $\sigma$ for rank-$1$ models.

We then explore the high-dimensional settings with imbalanced tabular
dimensions ($p_{1}\neq p_{2}$) and/or multiple components ($r>1$). We fix
$\lambda_{min}=80$, $\sigma=\tau=1$, $n=30$, and consider the following four
specific scenarios: i@: $p_{1}=20,p_{2}=500,r=1$; ii@:
$p_{1}=20,p_{2}=500,r=2$; iii@: $p_{1}=100,p_{2}=500,r=1$; iv@:
$p_{1}=100,p_{2}=100,r=3$. FPCA and CP, the two baseline methods, are
implemented as follows in the cases of $r>1$. For FPCA, we unfold the tensor
data to matrix similarly as the previous experiment and apply FPCA to obtain
the first $r$ eigenfunctions. For CP, we apply the classic power iteration
with random initialization proposed by Anandkumar et al., 2014b  (the number
of initializers is set to 20), and the same interpolation procedure is applied
to obtain a functional observation as we did in rank-$1$ situations. The
estimation errors of all singular vectors/functions and the corresponding
execution time are reported in Table 1. As we can see, FTSVD achieves better
performances on all the four scenarios than CP and FPCA, particularly in
singular function estimation. For the tabular modes, our method has nearly the
same performance as CP. This is because in the high-dimensional setting with
large $p$, our method shares the same dominant rate as CP, i.e.,
$O(\sqrt{p})$. In addition, our method has nearly the same execution time as
the classic CP.

Scenarios | Method | ${\rm dist}(a,\hat{a})$ | ${\rm dist}(b,\hat{b})$ | ${\rm dist}(\xi,\hat{\xi})$ | execution time (s)
---|---|---|---|---|---
i@ | FTSVD |  0.121 (0.025) |  0.482 (0.034) |  0.140 (0.070) | 0.67
| CP | 0.123 (0.025) | 0.485 (0.040) | 0.524 (0.224) | 0.67
| FPCA | NA | NA | 0.395 (0.122) | 1.96
ii@ | FTSVD |  0.262 (0.194) |  0.475 (0.118) |  0.432 (0.200) | 0.55
| CP | 0.377 (0.188) | 0.551 (0.118) | 0.840 (0.207) | 3.39
| FPCA | NA | NA | 0.578 (0.124) | 2.34
iii@ | FTSVD |  0.307 (0.158) |  0.513 (0.117) |  0.186 (0.189) | 3.34
| CP | 0.310 (0.159) | 0.516 (0.118) | 0.528 (0.245) | 3.28
| FPCA | NA | NA | 0.628 (0.238) | 12.82
iv@ | FTSVD |  0.225 (0.109) |  0.223 (0.107) |  0.397 (0.134) | 1.00
| CP | 0.520 (0.139) | 0.521 (0.138) | 0.709 (0.221) | 3.27
| FPCA | NA | NA | 0.691 (0.075) | 2.69
Table 1: Estimation errors and execution time of $\hat{a},\hat{b},\hat{\xi}$
in high-dimensional settings with imbalanced tabular dimensions and/or
multiple PCs. The standard errors are shown in parentheses and the minimal
error in each scenario is shown in bold. For multiple PCs scenarios, the
average estimation errors of all singular vectors/functions are presented. NAs
appear since FPCA does not yield estimates of $a$ or $b$.

###  5.3 Real Data Analysis

We illustrate the merit of our method through two real world datasets: a
longitudinal microbiome study and world-wide crop production data. It turns
out that the our algorithm yields several interesting results on these two
datasets.

####  5.3.1 Longitudinal Microbiome Study

We first study the change of fecal microbial composition for new-born infants
by applying our methods to the Early Childhood Antibiotics and the Microbiome
(ECAM) dataset published by Bokulich et al., (2016) (Qiita ID 10249). We
consider the 42 infants with multiple fecal microbiome measurements from birth
over the first 2 years of life. The infants have fecal microbiome sampled
monthly in the first year and bi-monthly in the second year. Among the 42
infants, 24 are vaginally delivered and 18 are Cesarean delivered. A natural
question is whether the delivery method affects the composition and
development of microbiome in the infants’ gut environment.

We focus on the 50 bacterial genera with non-zero read counts for more than
10% of all the samples. The data can be organized as an order-3 count tensor
$\bar{\mathbfcal{Y}}\in\mathbb{N}^{42\times 50\times 19}$, where the three
modes represent different subjects (i.e., infants), bacterial genus and
sampling time respectively. To account for the variation in sequencing depth,
we transform the count data to the log-composition after $.5$ is added to
every count:

$\widetilde{\mathbfcal{Y}}\in{\mathbb{R}}^{42\times 50\times
19},\qquad\widetilde{\mathbfcal{Y}}_{ijk}=\log\left(\frac{\bar{\mathbfcal{Y}}_{ijk}+.5}{\sum_{j^{\prime}=1}^{50}\left(\bar{\mathbfcal{Y}}_{ij^{\prime}k}+.5\right)}\right).$

We apply the proposed method with sequential spectral initialization
(Algorithm 2) to estimate the leading singular components for each mode. For
the purpose of illustration, we focus on the first three components $r=3$. The
singular vectors on the subject mode (i.e., $\\{\hat{a}_{s}\\}_{s=1}^{3}$) are
visualized using bi-plots in Figure 3, where a well separation of infants with
different delivery methods can be seen in the first three components,
particularly the second and the third component. We next present the estimated
singular functions on the time mode (i.e., $\\{\hat{\xi}_{s}\\}_{s=1}^{3}$) in
the left panel of Figure 4. The first singular function on the time mode is
nearly constant, suggesting that the first singular component
$\hat{a}_{1}\otimes\hat{b}_{1}\otimes\hat{\xi}_{1}$ may capture the time-
invariant scaling information of each genus. The second singular function on
the time mode is nearly monotone, which represents an increasing or decreasing
trend of some bacterial abundance that plateaued in the second year. The third
singular function on the time mode increases in the first three months and
decrease around 12 months, suggesting that it may characterize bacteria
abundance that is different between the first and second year of life. The two
panels on the right of Figure 4 demonstrate the mean and error bands of the
observed bacteria trajectories for all subjects grouped by delivery method,
where the trajectory of each subject is obtained by
$\sum_{j=1}^{50}\hat{b}_{lj}\widetilde{\mathbfcal{Y}}_{ijk}$, i.e., the
weighted average of the observed trajectories of all bacteria genera using the
feature singular vector of Component 2 or 3 as weights. The trajectories of
infants from different delivery methods are well separated in Component 3 and
in the latter part of Component 2. Note that all the singular
vectors/functions are learned without the information of delivery methods. Our
result indicates that the 50 bacteria genera can be reduced to two aggregated
bacteria genera using the singular vector on the bacterial mode and still
achieve high predictive ability of the delivery methods.

Figure 3: Biplot of three singular vectors on subject mode from ECAM data.
Note that each point represents an infant with color indicating the delivery
method.

Figure 4: Three singular functions from ECAM data and the aggregated observed
trajectory with respect to the second and third singular vector on bacteria
mode. Error bands of trajectories are obtained using mean $\pm 1.64\times$
(standard error of the mean).

###  5.4 Crop Production Data Analysis

We then use the world-wide crops production data (available on
http://www.fao.org/faostat/en/#data/QC) as the second illustration. The
dataset contains the annual production of 173 products between 1961 and 2019
from different countries and areas in the world. For simplicity, we only
consider 19 continent-based areas including Eastern Asian, Northern America,
West Europe, Central Africa, etc. We also select the most widely-planted 26
crops for our analysis. Therefore, the data we consider can be organized as a
$19$-by-$26$-by-$59$ tensor $\widetilde{\mathbfcal{Y}}$, with
$\widetilde{\mathbfcal{Y}}_{ijk}$ representing the production of the $j$th
crop item on the $i$th area in the $k$th year.

We apply the proposed algorithm with rank $r=2$ to
$\widetilde{\mathbfcal{Y}}$. The estimated singular vectors and functions for
each component are presented in Figure 5. The first component captures the
overall production of areas and items through the years. The magnitudes of the
first singular vectors on area and item modes (i.e., x-coordinate of the bi-
plot) coincide with the overall production of the areas and items,
respectively, with Northern America and Eastern Asia being the top two crop
production areas, and maze, rice paddy, wheat being the three major food
crops. The increasing first singular function on time mode coupled with
positive singular vectors on area and item modes reflects the increasing trend
in overall production for most areas and crops. The second component
characterizes the variation in production across different area-item pairs.
Such variation is quantified by the second singular function on time mode
coupled with the second item singular vectors on item and area modes (i.e.,
y-coordinate of the bi-plots). For example, Northern America has a large
negative value in Component 2, with the same sign as maize but different from
rice paddy. This implies maze takes a larger share in Northern America
compared to its share in other areas, while rice paddy takes a smaller share
in Northern America. An opposite conclusion can also be made on the three
Asian areas, whose large positive values in the second component have the same
sign as rice paddy, indicating the larger share of rice paddy in Asia compared
to other areas. The small difference between the two singular functions on
time mode indicate that the variation across area-item pairs are constant
across time.

Figure 5: Estimated singular vectors and functions from crop production data.
The left and central panels are biplots for the two estimated singular vectors
on area and item modes respectively. Points with significant scorings are
labelled. The right panel plots the estimated singular function on the time
mode.

##  6 Discussion

This paper proposed the functional tensor singular value decomposition for
high-dimensional high-order longitudinal tensor data. Compared with the
classical power iteration schemes for tensors with only tabular modes, the
proposed method can handle a functional mode based on the techniques from
reproducing kernel Hilbert space (RKHS) and achieve accurate statistical
estimation for singular vectors and functions in both theory and numerical
studies. We apply the algorithms to two real-life multivariate longitudinal
studies and obtain interesting interpretable findings that provide new
scientific insights.

Although the presentation of this paper focuses on order-3 tensor with two
tabular modes and one functional mode, the algorithms and all the theoretical
results can be generalized to arbitrary order-$d$ tensors with $(d-1)$ tabular
modes and one functional modes. It is an interesting future direction to
further study functional tensor SVD for tensors with multiple functional
modes, which may widely appear in spatial-temporal data analysis and imaging
processing. In addition, this paper focuses on the functional tensor SVD when
all the observations are from a grid of time points. It is interesting to
further study the same topic under irregular sampling schemes where the
sampling time points $t_{ijk}$ may be different across different units and
variables.

It is also worth mentioning that there are several recent results that
directly solve the tabular tensor decomposition by (accelerated) gradient
method with convergence guarantee (Cai et al.,, 2019; Han et al.,, 2021; Tong
et al.,, 2021). Such methods usually have an advantage over power iteration as
they directly aim at the optimization of likelihood, and can successfully
remove the effect of incoherence from the final estimation error bound.
However, it is unclear that whether and how these gradient-based algorithms
and analysis techniques can be applied to the functional setting, which might
be another interesting research direction in the future.

##  7 Proof

We provide the proof of Theorem 1 in this section and defer the proofs of
other technical results to the supplementary materials. In order to prove
Theorem 1, we introduce the following events:

  1. ($A_{1}$)

$\left|\int_{0}^{1}f(s)ds-\frac{1}{n}\sum_{k=1}^{n}f(s_{k})\right|\leq
C_{1}\left(\zeta_{n}\|f\|_{\mathcal{L}^{2}}+\zeta_{n}^{2}\|f\|_{\mathcal{H}}\right),\qquad\forall
f\in\mathcal{H}.$ (23)

  2. ($A_{2}$)

$\begin{split}\left|\left\|f\right\|_{\mathcal{L}^{2}}-\|f\|_{n}^{2}\right|\leq
C_{2}\|f\|_{\mathcal{H}}\left(\zeta_{n}\|f\|_{\mathcal{L}^{2}}+\zeta_{n}^{2}\|f\|_{\mathcal{H}}\right),\qquad\forall
f\in\mathcal{H}.\end{split}$ (24)

  3. ($A_{3}$)

$\begin{split}&\sup_{a\in{\mathbb{S}}^{p_{1}-1},b\in{\mathbb{S}}^{p_{2}-1}}\frac{1}{n}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}a_{i}b_{j}f(s_{k})\varepsilon_{ijk}\\\
&\leq C_{3}\tau\left\\{\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log
n}{n}\right)\|f\|_{\mathcal{H}}+\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\|f\|_{\mathcal{L}^{2}}\right\\},\qquad\forall
f\in\mathcal{H}.\end{split}$ (25)

Here recall that
$\|f\|_{n}:=\left(n^{-1}\sum_{k=1}^{n}f^{2}(s_{k})\right)^{1/2}$ and
$C_{1},C_{2},C_{3}$ are some universal constants. By Lemmas 4, 5, and 6, we
have

${\mathbb{P}}\left(A_{1}\right)\geq
1-Cn^{-9},\qquad{\mathbb{P}}\left(A_{2}\right)\geq
1-Cn^{-9},\qquad{\mathbb{P}}\left(A_{3}\right)\geq 1-Cn^{-9}-C\log
n\cdot\exp(-c(p_{1}+p_{2})).$

The following analyses will be conducted under the event $A:=A_{1}\cap
A_{2}\cap A_{3}$, which holds with probability at least $1-C\log
n\cdot\exp(-c(p_{1}+p_{2}))$.

We next introduce the following two Lemmas to establish the functional and
tabular error contractions for the proposed RKHS-constraint power iteration
(Algorithm 1).

######  Lemma 1.

Suppose the conditions in Theorem 1 and $(A_{1})-(A_{3})$ hold. For any
$t=0,\ldots,T-1$, if ${\rm dist}(a^{(t)},a_{l}),{\rm
dist}(b^{(t)},b_{l})\leq\delta^{(t)}$ for some $\delta^{(t)}\leq c_{0}/(\kappa
r)$, then $\|\xi^{(t+1)}\|_{\mathcal{H}}\leq 2\kappa C_{\xi}$ and

$\begin{split}\max\\{{\rm dist}(\xi^{(t+1)},\xi_{l}),{\rm
dist}(\xi^{(t+1)}_{n},(\xi_{l})_{n})\\}&\leq\frac{1}{2}\delta^{(t)}+C\kappa(\mu(r-1)+\zeta_{n})\\\
&\quad+C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}.\end{split}$

###### Proof.

See Section A.1. ∎

######  Lemma 2.

Suppose the conditions in Theorem 1 and $(A_{1})-(A_{3})$ hold. For any
$t=0,\ldots,T-1$, if ${\rm dist}(b^{(t)},b_{l}),{\rm
dist}((\xi^{(t)})_{n},(\xi_{l})_{n})\leq\delta^{(t)}$ for some
$\delta^{(t)}\leq c_{0}/(\kappa r)$ and $\|\xi^{(t)}\|_{\mathcal{H}}\leq
2\kappa C_{\xi}$. Then,

$\begin{split}{\rm dist}(a^{(t+1)},a_{l})\leq
C\sqrt{r-1}\kappa(\mu+\delta^{(t)})(\mu+\delta^{(t)}+\zeta_{n})+C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}.\end{split}$

###### Proof.

See Section A.1. ∎

Note that we can obtain almost the same error contraction for $b^{(t+1)}$ as
an analogue of Lemma 2. Now we are ready to prove Theorem 1.

###### Proof of Theorem 1.

Without loss of generality, we assume $l=1$. We first prove by induction that
for any $t=0,\ldots,T$,

$\max\left\\{{\rm dist}(a^{(t)},a_{1}),{\rm
dist}(b^{(t)},b_{1})\right\\}\leq\frac{c_{0}}{\kappa r}.$ (26) ${\rm
dist}((\xi^{(t)})_{n},(\xi_{1})_{n})\leq\frac{c_{0}}{\kappa r}.$ (27)
$\|\xi^{(t)}\|_{\mathcal{H}}\leq 2\kappa C_{\xi}.$ (28)

We claim that (26) implies (27) and (28) for any $t=0,\ldots,T$. By
Assumptions 2 and 3, we have

$\mu\leq\frac{c}{\kappa^{2}r^{2}},\qquad\zeta_{n}\leq\frac{c}{\kappa^{2}r},\quad\text{and}\quad\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\leq\frac{c_{0}}{4\kappa
r}.$

Assuming (26) and Applying Lemma 1, one immediately obtains

${\rm dist}(\xi_{n}^{(t)},(\xi_{1})_{n})\leq\frac{2c_{0}}{\kappa
r}+C\kappa\cdot\frac{3c}{\kappa^{2}r}+\frac{c_{0}}{4\kappa
r}\leq\frac{c_{0}}{\kappa r}$

and $\|\xi^{(t)}\|\leq 2\kappa C_{\xi}$ for any $t=0,\ldots,T$. So (27) and
(28) are proved. Now assuming (26), (27) and (28) hold for some $t=t_{0}<T$,
we show that (26) also holds for $t=t_{0}+1$. By Lemma 2, we have

$\begin{split}{\rm dist}(a^{(t_{0}+1)},a_{1})\leq
C\sqrt{r-1}\kappa\left(\frac{c}{\kappa^{2}r^{2}}+\frac{c_{0}}{\kappa
r}\right)\left(\frac{c}{\kappa^{2}r^{2}}+\frac{c_{0}}{\kappa
r}+\frac{c}{\kappa^{2}r}\right)+\frac{c_{0}}{4\kappa r}\leq\frac{c_{0}}{\kappa
r}.\end{split}$

Similarly one can also show that ${\rm
dist}(b^{(t_{0}+1)},b_{1})\leq\frac{c_{0}}{\kappa r}$ and thus (27) holds for
$t=t_{0}+1$. Finally, by Assumption 4, (26) holds at $t=0$. Therefore, we have
proved that (26) and (27) hold for all $t=0,\ldots,T$ by induction.

In addition, note that by the upper bounds of $\mu,\zeta_{n}$ and
$\delta^{(t)}$, we also have

$\sqrt{r-1}\kappa(\mu+\delta^{(t)})(\mu+\delta^{(t)}+\zeta_{n})\leq\frac{1}{2C}\delta^{(t)}+\kappa(\mu(r-1)+\zeta_{n}).$

Therefore, under the same condition of Lemma 2, we further have

$\begin{split}{\rm
dist}(a^{(t+1)},a_{1})&\leq\frac{1}{2}\delta^{(t)}+C\kappa(\mu(r-1)+\zeta_{n})\\\
&\quad+C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}.\end{split}$
(29)

Now we define

$e^{(t)}:=\max\left\\{{\rm dist}(a^{(t)},a_{1}),{\rm dist}(b^{(t)},b_{1}),{\rm
dist}(\xi^{(t)},\xi_{1})\right\\}.$

Combining (29) with Lemma 1, we obtain that

$e^{(t+1)}\leq\frac{1}{2}e^{(t)}+C\left\\{\kappa\left(\mu(r-1)+\zeta_{n}\right)+\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\},\qquad\forall
t\geq 0.$

Then it follows by induction that

$e^{(t)}\leq
2^{-t}e^{(0)}+2C\left\\{\kappa(\mu(r-1)+\zeta_{n})+\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\},\qquad\forall
t\geq 0.$

In particular, when $r=1$, Lemma 2 directly implies that

${\rm dist}(a^{(t)},a_{1})\leq
C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}$

for any $t\geq 1$. Now the proof of this theorem is completed. ∎

## References

  * Abraham et al., (2003) Abraham, C., Cornillon, P.-A., Matzner-Løber, E., and Molinari, N. (2003).  Unsupervised curve clustering using b-splines.  Scandinavian journal of statistics, 30(3):581–595.
  * Adamczak, (2008) Adamczak, R. (2008).  A tail inequality for suprema of unbounded empirical processes with applications to markov chains.  Electronic Journal of Probability, 13:1000–1034.
  * Adler and Taylor, (2009) Adler, R. J. and Taylor, J. E. (2009).  Random fields and geometry.  Springer Science & Business Media.
  * Allen, (2013) Allen, G. I. (2013).  Multi-way functional principal components analysis.  In 2013 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 220–223. IEEE.
  * (5) Anandkumar, A., Ge, R., Hsu, D., Kakade, S. M., and Telgarsky, M. (2014a).  Tensor decompositions for learning latent variable models.  The Journal of Machine Learning Research, 15(1):2773–2832.
  * (6) Anandkumar, A., Ge, R., and Janzamin, M. (2014b).  Guaranteed non-orthogonal tensor decomposition via alternating rank-$1$ updates.  arXiv preprint arXiv:1402.5180.
  * Aue et al., (2015) Aue, A., Norinho, D. D., and Hörmann, S. (2015).  On the prediction of stationary functional time series.  Journal of the American Statistical Association, 110(509):378–392.
  * Bartlett et al., (2005) Bartlett, P. L., Bousquet, O., and Mendelson, S. (2005).  Local rademacher complexities.  The Annals of Statistics, 33(4):1497–1537.
  * Bartlett and Mendelson, (2002) Bartlett, P. L. and Mendelson, S. (2002).  Rademacher and gaussian complexities: Risk bounds and structural results.  Journal of Machine Learning Research, 3(Nov):463–482.
  * Berrendero et al., (2011) Berrendero, J. R., Justel, A., and Svarc, M. (2011).  Principal components for multivariate functional data.  Computational Statistics & Data Analysis, 55(9):2619–2634.
  * Bi et al., (2018) Bi, X., Qu, A., and Shen, X. (2018).  Multilayer tensor factorization with applications to recommender systems.  The Annals of Statistics, 46(6B):3308–3333.
  * Bokulich et al., (2016) Bokulich, N. A., Chung, J., Battaglia, T., Henderson, N., Jay, M., Li, H., Lieber, A. D., Wu, F., Perez-Perez, G. I., and Chen, Y. (2016).  Antibiotics, birth mode, and diet shape microbiome maturation during early life.  Science translational medicine, 8(343):343ra82–343ra82.
  * Brezis, (2010) Brezis, H. (2010).  Functional analysis, Sobolev spaces and partial differential equations.  Springer Science & Business Media.
  * Cai et al., (2019) Cai, C., Li, G., Poor, H. V., and Chen, Y. (2019).  Nonconvex low-rank symmetric tensor completion from noisy data.  arXiv preprint arXiv:1911.04436.
  * Cai et al., (2020) Cai, T. T., Han, R., and Zhang, A. R. (2020).  On the non-asymptotic concentration of heteroskedastic wishart-type matrix.  arXiv preprint arXiv:2008.12434.
  * Cai and Yuan, (2011) Cai, T. T. and Yuan, M. (2011).  Optimal estimation of the mean function based on discretely sampled functional data: Phase transition.  The annals of statistics, 39(5):2330–2355.
  * Cai and Yuan, (2012) Cai, T. T. and Yuan, M. (2012).  Minimax and adaptive prediction for functional linear regression.  Journal of the American Statistical Association, 107(499):1201–1216.
  * Chen et al., (2020) Chen, E. Y., Tsay, R. S., and Chen, R. (2020).  Constrained factor models for high-dimensional matrix-variate time series.  Journal of the American Statistical Association, 115(530):775–793.
  * Chen et al., (2017) Chen, K., Delicado Useros, P. F., and Müller, H.-G. (2017).  Modelling function-valued stochastic processes, with applications to fertility dynamics.  Journal of the Royal Statistical Society. Series B, Statistical Methodology, 79(1):177–196.
  * Chen et al., (2021) Chen, R., Xiao, H., and Yang, D. (2021).  Autoregressive models for matrix-valued time series.  Journal of Econometrics, 222(1):539–560.
  * Chiou, (2012) Chiou, J.-M. (2012).  Dynamical functional prediction and classification, with application to traffic flow prediction.  The Annals of Applied Statistics, 6(4):1588–1614.
  * Chiou et al., (2014) Chiou, J.-M., Chen, Y.-T., and Yang, Y.-F. (2014).  Multivariate functional principal component analysis: A normalization approach.  Statistica Sinica, pages 1571–1596.
  * Chiou and Li, (2007) Chiou, J.-M. and Li, P.-L. (2007).  Functional clustering and identifying substructures of longitudinal data.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(4):679–699.
  * Crainiceanu et al., (2009) Crainiceanu, C. M., Staicu, A.-M., and Di, C.-Z. (2009).  Generalized multilevel functional regression.  Journal of the American Statistical Association, 104(488):1550–1561.
  * Dawid, (1981) Dawid, A. P. (1981).  Some matrix-variate distribution theory: notational considerations and a bayesian application.  Biometrika, 68(1):265–274.
  * (26) De Lathauwer, L., De Moor, B., and Vandewalle, J. (2000a).  A multilinear singular value decomposition.  SIAM journal on Matrix Analysis and Applications, 21(4):1253–1278.
  * (27) De Lathauwer, L., De Moor, B., and Vandewalle, J. (2000b).  On the best rank-1 and rank-(r 1, r 2,…, rn) approximation of higher-order tensors.  SIAM journal on Matrix Analysis and Applications, 21(4):1324–1342.
  * Di et al., (2014) Di, C., Crainiceanu, C. M., and Jank, W. S. (2014).  Multilevel sparse functional principal component analysis.  Stat, 3(1):126–143.
  * Di et al., (2009) Di, C.-Z., Crainiceanu, C. M., Caffo, B. S., and Punjabi, N. M. (2009).  Multilevel functional principal component analysis.  The annals of applied statistics, 3(1):458.
  * Dubin and Müller, (2005) Dubin, J. A. and Müller, H.-G. (2005).  Dynamical correlation for multivariate longitudinal data.  Journal of the American Statistical Association, 100(471):872–881.
  * Fan et al., (2015) Fan, Y., James, G. M., and Radchenko, P. (2015).  Functional additive regression.  The Annals of Statistics, 43(5):2296–2325.
  * Ferraty and Vieu, (2006) Ferraty, F. and Vieu, P. (2006).  Nonparametric functional data analysis: theory and practice.  Springer Science & Business Media.
  * Friedman et al., (2001) Friedman, J., Hastie, T., and Tibshirani, R. (2001).  The elements of statistical learning, volume 1.  Springer series in statistics New York.
  * Goldenberg et al., (2009) Goldenberg, A., Zheng, A. X., Fienberg, S. E., and Airoldi, E. M. (2009).  A survey of statistical network models.  arXiv preprint arXiv:0912.5410.
  * Greven et al., (2010) Greven, S., Crainiceanu, C., Caffo, B., and Reich, D. (2010).  Longitudinal functional principal component analysis.  Electronic Journal of Statistics, 4(none):1022 – 1054.
  * Gu, (2013) Gu, C. (2013).  Smoothing spline ANOVA models, volume 297.  Springer Science & Business Media.
  * Guo, (2002) Guo, W. (2002).  Functional mixed effects models.  Biometrics, 58(1):121–128.
  * Han et al., (2019) Han, Q., Wang, T., Chatterjee, S., and Samworth, R. J. (2019).  Isotonic regression in general dimensions.  Annals of Statistics, 47(5):2440–2471.
  * Han et al., (2020) Han, R., Luo, Y., Wang, M., and Zhang, A. R. (2020).  Exact clustering in tensor block model: Statistical optimality and computational limit.  arXiv preprint arXiv:2012.09996.
  * Han et al., (2021) Han, R., Willett, R., and Zhang, A. R. (2021).  An optimal statistical and computational framework for generalized tensor estimation.  The Annals of Statistics, to appear.
  * Hao et al., (2020) Hao, B., Zhang, A., and Cheng, G. (2020).  Sparse and low-rank tensor estimation via cubic sketchings.  IEEE Transactions on Information Theory.
  * Happ and Greven, (2018) Happ, C. and Greven, S. (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association, 113(522):649–659.
  * Hasenstab et al., (2017) Hasenstab, K., Scheffler, A., Telesca, D., Sugar, C. A., Jeste, S., DiStefano, C., and Şentürk, D. (2017).  A multi-dimensional functional principal components analysis of eeg data.  Biometrics, 73(3):999–1009.
  * Hillar and Lim, (2013) Hillar, C. J. and Lim, L.-H. (2013).  Most tensor problems are np-hard.  Journal of the ACM (JACM), 60(6):1–39.
  * Hong et al., (2020) Hong, D., Kolda, T. G., and Duersch, J. A. (2020).  Generalized canonical polyadic tensor decomposition.  SIAM Review, 62(1):133–163.
  * Hsing and Eubank, (2015) Hsing, T. and Eubank, R. (2015).  Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997.  John Wiley & Sons.
  * Hu and Yao, (2021) Hu, X. and Yao, F. (2021).  Dynamic principal subspaces with sparsity in high dimensions.  arXiv preprint arXiv:2104.03087.
  * James et al., (2000) James, G. M., Hastie, T. J., and Sugar, C. A. (2000).  Principal component models for sparse functional data.  Biometrika, 87(3):587–602.
  * Jing et al., (2020) Jing, B.-Y., Li, T., Lyu, Z., and Xia, D. (2020).  Community detection on mixture multi-layer networks via regularized tensor decomposition.  arXiv preprint arXiv:2002.04457.
  * Kimeldorf and Wahba, (1971) Kimeldorf, G. and Wahba, G. (1971).  Some results on tchebycheffian spline functions.  Journal of mathematical analysis and applications, 33(1):82–95.
  * Kivelä et al., (2014) Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J. P., Moreno, Y., and Porter, M. A. (2014).  Multilayer networks.  Journal of complex networks, 2(3):203–271.
  * Koltchinskii and Yuan, (2010) Koltchinskii, V. and Yuan, M. (2010).  Sparsity in multiple kernel learning.  The Annals of Statistics, 38(6):3660–3695.
  * Lei et al., (2020) Lei, J., Chen, K., and Lynch, B. (2020).  Consistent community detection in multi-layer network data.  Biometrika, 107(1):61–73.
  * Li et al., (2018) Li, X., Xu, D., Zhou, H., and Li, L. (2018).  Tucker tensor regression and neuroimaging analysis.  Statistics in Biosciences, 10(3):520–545.
  * Li and Hsing, (2010) Li, Y. and Hsing, T. (2010).  Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data.  The Annals of Statistics, 38(6):3321–3351.
  * Martino et al., (2021) Martino, C., Shenhav, L., Marotz, C. A., Armstrong, G., McDonald, D., Vázquez-Baeza, Y., Morton, J. T., Jiang, L., Dominguez-Bello, M. G., and Swafford, A. D. (2021).  Context-aware dimensionality reduction deconvolutes gut microbial community dynamics.  Nature biotechnology, 39(2):165–168.
  * Mendelson, (2002) Mendelson, S. (2002).  Geometric parameters of kernel machines.  In International Conference on Computational Learning Theory, pages 29–43. Springer.
  * Mercer, (1909) Mercer, J. (1909).  Xvi. functions of positive and negative type, and their connection the theory of integral equations.  Philosophical transactions of the royal society of London. Series A, containing papers of a mathematical or physical character, 209(441-458):415–446.
  * Micchelli and Wahba, (1981) Micchelli, C. A. and Wahba, G. (1981).  Design problems for optimal surface interpolation.  In Ziegler, Z., editor, Approximation Theory and Applications, pages 329–347. Academic Press, New York.
  * Morris and Carroll, (2006) Morris, J. S. and Carroll, R. J. (2006).  Wavelet-based functional mixed models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(2):179–199.
  * Oseledets, (2011) Oseledets, I. V. (2011).  Tensor-train decomposition.  SIAM Journal on Scientific Computing, 33(5):2295–2317.
  * Pisier, (1983) Pisier, G. (1983).  Some applications of the metric entropy condition to harmonic analysis.  In Banach Spaces, Harmonic Analysis, and Probability Theory, pages 123–154. Springer.
  * Ramsay and Silverman, (2006) Ramsay, J. and Silverman, B. (2006).  Functional Data Analysis.  Springer Science & Business Media.
  * Raskutti et al., (2012) Raskutti, G., J Wainwright, M., and Yu, B. (2012).  Minimax-optimal rates for sparse additive models over kernel classes via convex programming.  Journal of Machine Learning Research, 13(2).
  * Rice and Silverman, (1991) Rice, J. A. and Silverman, B. W. (1991).  Estimating the mean and covariance structure nonparametrically when the data are curves.  Journal of the Royal Statistical Society: Series B (Methodological), 53(1):233–243.
  * Richard and Montanari, (2014) Richard, E. and Montanari, A. (2014).  A statistical model for tensor PCA.  In Advances in Neural Information Processing Systems, pages 2897–2905.
  * Rudelson and Vershynin, (2013) Rudelson, M. and Vershynin, R. (2013).  Hanson-wright inequality and sub-gaussian concentration.  Electronic Communications in Probability, 18:1–9.
  * Scheffler et al., (2020) Scheffler, A., Telesca, D., Li, Q., Sugar, C. A., Distefano, C., Jeste, S., and Şentürk, D. (2020).  Hybrid principal components analysis for region-referenced longitudinal functional eeg data.  Biostatistics, 21(1):139–157.
  * Sun and Li, (2019) Sun, W. W. and Li, L. (2019).  Dynamic tensor clustering.  Journal of the American Statistical Association, 114(528):1894–1907.
  * Sun et al., (2017) Sun, W. W., Lu, J., Liu, H., and Cheng, G. (2017).  Provable sparse tensor decomposition.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3):899–916.
  * Sun et al., (2018) Sun, X., Du, P., Wang, X., and Ma, P. (2018).  Optimal penalized function-on-function regression under a reproducing kernel hilbert space framework.  Journal of the American Statistical Association, 113(524):1601–1611.
  * Tong et al., (2021) Tong, T., Ma, C., Prater-Bennette, A., Tripp, E., and Chi, Y. (2021).  Scaling and scalability: Provable nonconvex low-rank tensor estimation from incomplete measurements.  arXiv preprint arXiv:2104.14526.
  * Tsiligkaridis and Hero, (2013) Tsiligkaridis, T. and Hero, A. O. (2013).  Covariance estimation in high dimensions via kronecker product expansions.  IEEE Transactions on Signal Processing, 61(21):5347–5360.
  * van Delft et al., (2017) van Delft, A., Characiejus, V., and Dette, H. (2017).  A nonparametric test for stationarity in functional time series.  arXiv preprint arXiv:1708.05248.
  * Verbeke et al., (2014) Verbeke, G., Fieuws, S., Molenberghs, G., and Davidian, M. (2014).  The analysis of multivariate longitudinal data: a review.  Statistical methods in medical research, 23(1):42–59.
  * Vershynin, (2018) Vershynin, R. (2018).  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press.
  * Wahba, (1990) Wahba, G. (1990).  Spline models for observational data.  SIAM.
  * Wang et al., (2019) Wang, D., Liu, X., and Chen, R. (2019).  Factor models for matrix-valued high-dimensional time series.  Journal of econometrics, 208(1):231–248.
  * (79) Wang, D., Zhao, Z., Willett, R., and Yau, C. Y. (2020a).  Functional autoregressive processes in reproducing kernel hilbert spaces.  arXiv preprint arXiv:2011.13993.
  * (80) Wang, D., Zhao, Z., Yu, Y., and Willett, R. (2020b).  Functional linear regression with mixed predictors.  arXiv preprint arXiv:2012.00460.
  * (81) Wang, J., Wong, R. K., and Zhang, X. (2020c).  Low-rank covariance function estimation for multidimensional functional data.  Journal of the American Statistical Association, pages 1–14.
  * Wang et al., (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G. (2016).  Functional data analysis.  Annual Review of Statistics and Its Application, 3:257–295.
  * Wedin, (1972) Wedin, P.-Å. (1972).  Perturbation bounds in connection with singular value decomposition.  BIT Numerical Mathematics, 12(1):99–111.
  * Xu and Yin, (2013) Xu, Y. and Yin, W. (2013).  A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion.  SIAM Journal on imaging sciences, 6(3):1758–1789.
  * (85) Yao, F., Müller, H.-G., and Wang, J.-L. (2005a).  Functional data analysis for sparse longitudinal data.  Journal of the American statistical association, 100(470):577–590.
  * (86) Yao, F., Müller, H.-G., and Wang, J.-L. (2005b).  Functional linear regression analysis for longitudinal data.  The Annals of Statistics, pages 2873–2903.
  * Yuan and Cai, (2010) Yuan, M. and Cai, T. T. (2010).  A reproducing kernel hilbert space approach to functional linear regression.  The Annals of Statistics, 38(6):3412–3444.
  * Zhang and Han, (2019) Zhang, A. and Han, R. (2019).  Optimal sparse singular value decomposition for high-dimensional high-order data.  Journal of the American Statistical Association, pages 1708–1725.
  * Zhang and Xia, (2018) Zhang, A. and Xia, D. (2018).  Tensor SVD: Statistical and computational limits.  IEEE Transactions on Information Theory, 64(11):7311–7338.
  * Zhang and Zhou, (2020) Zhang, A. R. and Zhou, Y. (2020).  On the non-asymptotic and sharp lower tail bounds of random variables.  Stat, 9(1):e314.
  * Zhang et al., (2020) Zhang, C., Han, R., Zhang, A. R., and Voyles, P. M. (2020).  Denoising atomic resolution 4d scanning transmission electron microscopy data with tensor singular value decomposition.  Ultramicroscopy, 219:113123.
  * Zhang and Golub, (2001) Zhang, T. and Golub, G. H. (2001).  Rank-one approximation to high order tensors.  SIAM Journal on Matrix Analysis and Applications, 23(2):534–550.
  * Zhang et al., (2019) Zhang, Z., Allen, G. I., Zhu, H., and Dunson, D. (2019).  Tensor network factorizations: Relationships between brain structural connectomes and traits.  Neuroimage, 197:330–343.
  * (94) Zhou, J., Sun, W. W., Zhang, J., and Li, L. (2020a).  Partially observed dynamic tensor response regression.  arXiv preprint arXiv:2002.09735.
  * Zhou, (2014) Zhou, S. (2014).  Gemini: Graph estimation with matrix variate normal instances.  The Annals of Statistics, 42(2):532–562.
  * (96) Zhou, Y., Zhang, A. R., Zheng, L., and Wang, Y. (2020b).  Optimal high-order tensor svd via tensor-train orthogonal iteration.  arXiv preprint arXiv:2010.02482.
  * Zipunnikov et al., (2011) Zipunnikov, V., Caffo, B., Yousem, D. M., Davatzikos, C., Schwartz, B. S., and Crainiceanu, C. (2011).  Multilevel functional principal component analysis for high-dimensional data.  Journal of Computational and Graphical Statistics, 20(4):852–873.

Supplement to “Guaranteed Functional Tensor Singular

Value Decomposition”

Rungang Han,1 Pixu Shi,2 and Anru R. Zhang1,2

1Department of Statistics, University of Wisconsin-Madison

2Department of Biostatistics & Bioinformatics, Duke University

##  Appendix A Proofs of Main Results

###  A.1 Proofs of Lemmas in Theorem 1

###### Proof of Lemma 1.

Without loss of generality, we assume $l=1$. We focus on a particular step $t$
and assume $\langle a^{(t)},a_{1}\rangle,\langle b^{(t)},b_{1}\rangle\geq 0$
as the signs of $a^{(t)}$ and $b^{(t)}$ have no essential effect on the
algorithm and analysis.

Denote $\Delta(\cdot)=\tilde{\xi}^{(t+1)}(\cdot)-\lambda_{1}\xi_{1}(\cdot)$.
We first present several important inequalities that will be used throughout
the analysis. Firstly, since $\|\tilde{\xi}^{(t+1)}\|\leq\lambda_{max}C_{\xi}$
and $\|\xi_{1}\|_{\mathcal{H}}\leq C_{\xi}$, we have

$\left\|\Delta\right\|_{\mathcal{H}}\leq\|\tilde{\xi}^{(t+1)}\|_{\mathcal{H}}+\lambda_{1}\|\xi_{1}\|_{\mathcal{H}}\leq
2\lambda_{max}C_{\xi}.$ (30)

In addition, Condition $(A_{2})$ implies that for any $m\in[r]$,

$\begin{split}\|\Delta\|_{n}&\overset{\eqref{ineq:A-2}}{\leq}\|\Delta\|_{\mathcal{L}^{2}}+\sqrt{C_{2}\zeta_{n}\|\Delta\|_{\mathcal{H}}\cdot\|\Delta\|_{\mathcal{L}^{2}}}+\sqrt{C_{2}}\zeta_{n}\|\Delta\|_{\mathcal{H}}\\\
&\overset{(a)}{\leq}2\|\Delta\|_{\mathcal{L}^{2}}+C\zeta_{n}\|\Delta\|_{\mathcal{H}};\\\
\|\xi_{m}\|_{n}&\overset{\eqref{ineq:A-2}}{\leq}\|\xi_{m}\|_{\mathcal{L}^{2}}+\sqrt{C_{2}\zeta_{n}\|\xi_{m}\|_{\mathcal{H}}\cdot\|\xi_{1}\|_{\mathcal{L}^{2}}}+\sqrt{C_{2}}\zeta_{n}\|\xi_{m}\|_{\mathcal{H}}\\\
&\overset{(b)}{\leq}1+\sqrt{C_{2}C_{\xi}\zeta_{n}}+\sqrt{C_{2}}C_{\xi}\zeta_{n}\overset{(c)}{\leq}3/2.\end{split}$
(31)

Here $(a)$ comes from the Arithmetic-Geometric mean inequality; $(b)$ is due
to the assumption that $\|\xi_{m}\|_{\mathcal{L}^{2}}=1$ and
$\|\xi_{m}\|_{\mathcal{H}}\leq C_{\xi}$; and $(c)$ comes from the assumption
that $\zeta_{n}<c$ for sufficiently small constant.

Now we are ready for the proof. Since $\tilde{\xi}^{(t)}$ is the optimal
solution of (9) and $\lambda_{1}\xi_{1}$ is in the feasible set, we have

$\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}\left(\widetilde{\mathbfcal{Y}}_{ijk}-a_{i}^{(t)}b_{j}^{(t)}\tilde{\xi}^{(t)}(s_{k})\right)^{2}\leq\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}\left(\widetilde{\mathbfcal{Y}}_{ijk}-\lambda_{1}a_{i}^{(t)}b_{j}^{(t)}\xi_{1}(s_{k})\right)^{2}.$

Since

$\widetilde{\mathbfcal{Y}}_{ijk}=\sum_{m=1}^{r}\lambda_{m}(a_{m})_{i}(b_{m})_{j}\xi_{m}(s_{k})+{\mathbfcal{Z}}_{ijs_{k}}+\varepsilon_{ijk},$

the above inequality is equivalent to:

$\underbrace{\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}(a_{i}^{(t)})^{2}(b_{j}^{(t)})^{2}\Delta^{2}(s_{k})}_{(i)}\leq\underbrace{2\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}a_{i}^{(t)}b_{j}^{(t)}\Delta(s_{k})\tilde{\varepsilon}_{ij}(s_{k})}_{(ii)},$
(32)

where

$\tilde{\varepsilon}_{ij}(s_{k}):={\mathbfcal{Z}}_{ijs_{k}}+\varepsilon_{ijk}+\lambda_{1}\left((a_{1})_{i}(b_{1})_{j}-a_{i}^{(t)}b_{j}^{(t)}\right)\xi_{1}(s_{k})+\sum_{m=2}^{r}\lambda_{m}(a_{m})_{i}(b_{m})_{j}\xi_{m}(s_{k}).$
(33)

We start by providing a lower bound for $(i)$ in (32). By Condition $(A_{2})$,

$\begin{split}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}-\|\Delta\|_{n}^{2}&\overset{\eqref{ineq:A-2}}{\leq}C_{2}\left(\zeta_{n}\left\|\Delta\right\|_{\mathcal{H}}\cdot\left\|\Delta\right\|_{\mathcal{L}^{2}}+\zeta_{n}^{2}\|\Delta\|_{\mathcal{H}}^{2}\right)\\\
&\overset{(a)}{\leq}C_{2}\left(\frac{1}{4C_{2}}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}+(C_{2}+1)\zeta_{n}^{2}\|\Delta\|_{\mathcal{H}}^{2}\right)\\\
&\overset{\eqref{ineq:Delta-H-norm-
bound}}{\leq}\frac{1}{4}\|\Delta\|_{\mathcal{L}^{2}}^{2}+C\zeta_{n}^{2}\lambda_{max}^{2}.\end{split}$

Here $(a)$ comes from the Arithmetic-Geometric mean inequality. Therefore, we
have

$(i)=\left\|\Delta\right\|_{n}^{2}\geq\frac{3}{4}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}-C\zeta_{n}^{2}\lambda_{\max}^{2}.$
(34)

Now we give an upper bound for $(ii)$. By (33), we can further decompose
$(ii)=(ii_{1})+(ii_{2})+(ii_{3})+(ii_{4})$, where

$\begin{split}(ii_{1})&=2\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}a_{i}^{(t)}b_{j}^{(t)}\Delta(s_{k}){\mathbfcal{Z}}_{ijs_{k}},\\\
(ii_{2})&=2\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}a_{i}^{(t)}b_{j}^{(t)}\Delta(s_{k})\varepsilon_{ijk},\\\
(ii_{3})&=2\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}\lambda_{1}a_{i}^{(t)}b_{j}^{(t)}\left((a_{1})_{i}(b_{1})_{j}-a_{i}^{(t)}b_{j}^{(t)}\right)\Delta(s_{k})\xi_{1}(s_{k}),\\\
(ii_{4})&=2\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\frac{1}{n}\sum_{k=1}^{n}\sum_{m=2}^{r}\lambda_{l}a_{i}^{(t)}b_{j}^{(t)}(a_{m})_{i}(b_{m})_{j}\Delta(s_{k})\xi_{m}(s_{k}).\end{split}$

We bound these four terms separately.

  * •

$(ii_{1})$: First, recall the definitions of $\mathcal{E}$ and one immediately
has

$\left\|\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a_{i}^{(t)}b_{j}^{(t)}{\mathbfcal{Z}}_{ij\cdot}\right\|_{\infty}\leq\mathcal{E}.$

Then it follows that

$\begin{split}\frac{1}{2}\cdot(ii_{1})&\leq\left\|\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a_{i}^{(t)}b_{j}^{(t)}{\mathbfcal{Z}}_{ij\cdot}\right\|_{\infty}\cdot\left(\frac{1}{n}\sum_{k=1}^{n}\left|\Delta(s_{k})\right|\right)\\\
&\overset{(a)}{\leq}\mathcal{E}\cdot\|\Delta\|_{n}\\\
&\overset{\eqref{ineq:A-2}}{\leq}\mathcal{E}\cdot\sqrt{\|\Delta\|_{\mathcal{L}^{2}}^{2}+C_{2}\zeta_{n}\|\Delta\|_{\mathcal{L}^{2}}\cdot\|\Delta\|_{\mathcal{H}}+C_{2}\zeta_{n}^{2}\|\Delta\|_{\mathcal{H}}^{2}}\\\
&\overset{(b)}{\leq}\mathcal{E}\cdot\sqrt{2\|\Delta\|_{\mathcal{L}^{2}}^{2}+(C_{2}^{2}/4+C_{2})\zeta_{n}^{2}\|\Delta\|_{\mathcal{H}}^{2}}\\\
&\overset{(c)}{\leq}\frac{1}{8}\|\Delta\|_{\mathcal{L}}^{2}+(C_{2}^{2}/32+C_{2}/16)\zeta_{n}^{2}\|\Delta\|_{\mathcal{H}}^{2}+4\mathcal{E}^{2}\\\
&\overset{\eqref{ineq:Delta-H-norm-
bound}}{\leq}\frac{1}{8}\|\Delta\|_{\mathcal{L}^{2}}^{2}+C\left(\lambda_{max}^{2}\zeta_{n}^{2}+\mathcal{E}^{2}\right).\end{split}$
(35)

Here (a) comes from the inequality between $l_{1}$ and $l_{2}$ norm, and (b),
(c) come from the Arithmetic-Geometric mean inequality.

  * •

$(ii_{2})$: This term can be bounded by investigating the concentration of the
functional tensor spectral norm under the RKHS constraint given by Condition
$(A_{3})$:

$\begin{split}\frac{1}{2}\cdot(ii_{2})&\overset{\eqref{ineq:A-3}}{\leq}C\tau\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log
n}{n}\right)\|\Delta\|_{\mathcal{H}}+C\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\|\Delta\|_{\mathcal{L}^{2}}\\\
&\overset{\eqref{ineq:Delta-H-norm-
bound}}{\leq}C\lambda_{max}\tau\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log
n}{n}\right)+C\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\|\Delta\|_{\mathcal{L}^{2}}\\\
&\leq\frac{1}{8}\|\Delta\|_{\mathcal{L}^{2}}^{2}+C\tau^{2}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)^{2}+\lambda_{max}^{2}\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log^{2}n}{n}\right).\end{split}$
(36)

Here the last inequality is again an application of Arithmetic-Geometric mean
inequality.

  * •

$(ii_{3})$: Note that

$\frac{1}{2}\cdot(ii_{3})\leq\lambda_{max}\left|\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a_{i}^{(t)}b_{j}^{(t)}\left((a_{1})_{i}(b_{1})_{j}-a_{i}^{(t)}b_{j}^{(t)}\right)\right|\cdot\left|\frac{1}{n}\sum_{k=1}^{n}\Delta(s_{k})\xi_{1}(s_{k})\right|.$
(37)

On one hand, since $\langle a_{1},a^{(t)}\rangle,\langle
b_{1},b^{(t)}\rangle\geq 0$, we can bound

$\begin{split}&\left|\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a_{i}^{(t)}b_{j}^{(t)}\left((a_{1})_{i}(b_{1})_{j}-a_{i}^{(t)}b_{j}^{(t)}\right)\right|\\\
&=\left|\left(\sum_{i=1}^{n}(a_{1})_{i}a_{i}^{(t)}\right)\left(\sum_{j=1}^{n}(b_{1})_{j}b_{j}^{(t)}\right)-\left(\sum_{i=1}^{n}(a_{i}^{(t)})^{2}\right)\left(\sum_{j=1}^{n}(b_{j}^{(t)})^{2}\right)\right|\\\
&=1-\langle a_{1},a^{(t)}\rangle\langle b_{1},b^{(t)}\rangle\leq
1-\left(1-(\delta^{(t)})^{2}\right)=(\delta^{(t)})^{2}.\end{split}$ (38)

On the other hand,

$\begin{split}\left|\frac{1}{n}\sum_{k=1}^{n}\Delta(s_{k})\xi_{1}(s_{k})\right|&\leq\|\Delta\|_{n}\cdot\|\xi_{1}\|_{n}\\\
&\overset{\eqref{ineq:delta-xi-n-
bound}}{\leq}\frac{3}{2}\left(2\|\Delta\|_{\mathcal{L}^{2}}+C\zeta_{n}\|\Delta\|_{\mathcal{H}}\right)\\\
&\overset{\eqref{ineq:Delta-H-norm-
bound}}{\leq}3\|\Delta\|_{\mathcal{L}^{2}}+C\lambda_{max}\zeta_{n}.\end{split}$
(39)

Combining (37), (38) and (39), we obtain

$\begin{split}\frac{1}{2}\cdot(ii_{3})&\leq(\delta^{(t)})^{2}\left(3\lambda_{max}\|\Delta\|_{\mathcal{L}^{2}}+C\lambda_{max}^{2}\zeta_{n}^{2}\right)\\\
&\leq\frac{1}{8}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}+C\left((\delta^{(t)})^{4}\lambda_{max}^{2}+(\delta^{(t)})^{2}\lambda_{max}^{2}\zeta_{n}^{2}\right).\\\
&\leq\frac{1}{8}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}+\frac{c_{1}}{\kappa^{2}}(\delta^{(t)})^{2}\lambda_{max}^{2}.\end{split}$
(40)

Here the last inequality holds for some sufficiently small constant $c_{1}>0$
by the assumptions on $\delta^{(t)}$ and $\zeta_{n}$:

$\delta^{(t)}\leq c/r\kappa,\qquad\zeta_{n}\leq c/\sqrt{r}\kappa.$ (41)

  * •

$(ii_{4})$: For any $m=2,\ldots,r$,

$\begin{split}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a_{i}^{(t)}b_{j}^{(t)}(a_{m})_{i}(b_{m})_{j}&=\langle
a^{(t)},a_{m}\rangle\cdot\langle b^{(t)},b_{m}\rangle\\\ &\leq\left(|\langle
a,a_{m}\rangle|+\langle a^{(t)}-a,a_{m}\rangle\right)\cdot\left(|\langle
b,b_{m}\rangle|+\langle b^{(t)}-b,b_{m}\rangle\right)\\\ &\leq\left(|\langle
a,a_{m}\rangle|+\|a^{(t)}-a\|_{2}\right)\cdot\left(|\langle
b,b_{m}\rangle|+\|b^{(t)}-b\|_{2}\right)\\\ &\leq
C(\mu+\delta^{(t)})^{2}.\end{split}$

Here the last inequality comes from the incoherence assumption and the fact
that

$\left\|u-v\right\|_{2}\leq 2\sqrt{1-\langle u,v\rangle^{2}}=2{\rm dist}(u,v)$

for any two unit vectors $u,v$ with $\langle u,v\rangle\geq 0$. Then it
follows by the similar argument as (40) that

$\begin{split}\frac{1}{2}\cdot(ii_{4})&\leq
C(\mu+\delta^{(t)})^{2}(r-1)\left(\lambda_{max}\|\Delta\|_{\mathcal{L}^{2}}+\lambda_{max}^{2}\zeta_{n}\right)\\\
&\leq\frac{1}{8}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}+C\left((\mu+\delta^{(t)})^{4}(r-1)^{2}\lambda_{\max}^{2}+(\mu+\delta^{(t)})^{2}(r-1)\lambda_{max}^{2}\zeta_{n}\right)\\\
&\overset{\eqref{ineq:delta-zeta-
upper}}{\leq}\frac{1}{8}\left\|\Delta\right\|_{\mathcal{L}^{2}}^{2}+\frac{c_{1}}{\kappa^{2}}(\delta^{(t)})^{2}\lambda_{max}^{2}+C\mu^{2}(r-1)^{2}\lambda_{max}^{2}.\end{split}$
(42)

Combining (34) with (35), (36), (40) and (42), we obtain

$\begin{split}\|\Delta\|_{\mathcal{L}^{2}}^{2}&\leq
C\left(\lambda_{max}^{2}\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log^{2}n}{n}\right)+\mathcal{E}^{2}+\tau^{2}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)^{2}\right)\\\
&\quad+\frac{2c_{1}}{\kappa^{2}}(\delta^{(t)})^{2}\lambda_{max}^{2}+C\mu^{2}(r-1)^{2}\lambda_{max}^{2}.\end{split}$

Taking the square root, and further applying (31), we obtain that

$\begin{split}\|\Delta\|_{\mathcal{L}^{2}}&\leq
C\left(\lambda_{max}\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log^{2}n}{n}\right)+\mathcal{E}+\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right)\\\
&\quad+\frac{c}{\kappa}\delta^{(t)}\lambda_{max}+C\mu(r-1)\lambda_{max},\end{split}$
(43) $\begin{split}\|\Delta\|_{n}&\leq
2\|\Delta\|_{\mathcal{L}^{2}}+C\zeta_{n}\|\Delta\|_{\mathcal{H}}\\\ &\leq
C\left(\lambda_{max}\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log^{2}n}{n}\right)+\mathcal{E}+\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right)\\\
&\quad+\frac{c}{\kappa}\delta^{(t)}\lambda_{max}+C\mu(r-1)\lambda_{max}.\end{split}$
(44)

Combining (41) and the other assumptions on SNR and incoherence, i.e.,

$\frac{n}{\log^{2}n}\geq\kappa(p_{1}+p_{2}),\qquad\mu\leq\frac{c}{\kappa^{2}r^{2}},\qquad\lambda_{min}\geq
C\kappa
r\left\\{\mathcal{E}+\tau\left(\zeta_{n}+\sqrt{(p_{1}+p_{2})/n}\right)\right\\},$

we further conclude that

$\max\left\\{\left\|\Delta\right\|_{\mathcal{L}^{2}},\left\|\Delta\right\|_{n}\right\\}\leq\frac{1}{2}\lambda_{min}.$

Then it follows that

$\begin{split}\left\|\tilde{\xi}^{(t+1)}\right\|_{\mathcal{L}^{2}}&\geq\|\lambda_{1}\xi_{1}\|_{\mathcal{L}^{2}}-\|\Delta\|_{\mathcal{L}^{2}}\geq\lambda_{min}-\frac{1}{2}\lambda_{min}=\frac{1}{2}\lambda_{min}.\\\
\left\|\tilde{\xi}^{(t+1)}\right\|_{n}&\geq\|\lambda_{1}\xi_{1}\|_{n}-\|\Delta\|_{n}\\\
&\overset{\eqref{ineq:A-2}}{\geq}\lambda_{min}\left(\|\xi_{1}\|_{\mathcal{L}^{2}}-\sqrt{C_{2}\zeta_{n}\|\xi_{1}\|_{\mathcal{L}^{2}}\|\xi_{1}\|_{\mathcal{H}}}-\zeta_{n}\|\xi_{1}\|_{\mathcal{H}}\right)-\frac{1}{2}\lambda_{min}\\\
&\geq\frac{3}{4}\lambda_{min}-\frac{1}{2}\lambda_{min}=\frac{1}{4}\lambda_{min}.\end{split}$

Since
$\xi^{(t+1)}=\tilde{\xi}^{(t+1)}/\|\tilde{\xi}^{(t+1)}\|_{\mathcal{L}^{2}}$
and $\xi_{n}^{(t+1)}=\tilde{\xi}_{n}^{(t+1)}/\|\tilde{\xi}_{n}^{(t+1)}\|$, we
obtain

$\begin{split}{\rm
dist}(\xi^{(t+1)},\xi_{1})\leq\frac{\|\tilde{\xi}^{(t)}-\lambda_{1}\xi_{1}\|}{\|\tilde{\xi}^{(t+1)}\|_{\mathcal{L}^{2}}}\leq\frac{\|\Delta\|_{\mathcal{L}^{2}}}{\lambda_{min}},\\\
{\rm
dist}(\xi^{(t+1)}_{n},(\xi_{1})_{n})\leq\frac{\|\tilde{\xi}_{n}^{(t)}-\lambda_{1}(\xi_{1})_{n}\|}{\|\tilde{\xi}^{(t+1)}\|}\leq\frac{\|\Delta_{n}\|}{\lambda_{min}}.\end{split}$
(45)

Now the one-step error contraction is proved by combining (45) with (43) and
(44).

In addition, we also get that

$\left\|\xi^{(t+1)}\right\|_{\mathcal{H}}=\frac{\left\|\tilde{\xi}^{(t+1)}\right\|_{\mathcal{H}}}{\left\|\tilde{\xi}^{(t+1)}\right\|_{\mathcal{L}^{2}}}\leq\frac{C_{\xi}\lambda_{max}}{\frac{1}{2}\lambda_{min}}=2\kappa
C_{\xi}.$

∎

###### Proof of Lemma 2.

We follow the convention in the proof of Lemma 1 by assuming $l=1$ and
$\langle
b_{1},b^{(t)}\rangle,\langle\xi_{1},\xi^{(t)}\rangle_{\mathcal{L}^{2}}\geq 0$.
By definition, $a^{(t+1)}=\tilde{a}^{(t+1)}/\|\tilde{a}^{(t+1)}\|_{2}$, where

$\begin{split}(\tilde{a}^{(t+1)})_{i}&=\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k})\widetilde{\mathbfcal{Y}}_{ijk}\\\
&=\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k})\left(\lambda_{1}(a_{1})_{i}(b_{1})_{j}\xi_{1}(s_{k})+\sum_{m=2}^{r}\lambda_{m}(a_{m})_{i}(b_{m})_{j}\xi_{m}(s_{k})+{\mathbfcal{Z}}_{ijs_{k}}+\varepsilon_{ijk}\right).\end{split}$

Denote

$\bar{a}^{(t+1)}=\lambda_{1}\left(\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}(b_{1})_{j}\xi^{(t)}(s_{k})\xi_{1}(s_{k})\right)a_{1}=\left(\lambda_{1}\langle
b_{1},b^{(t)}\rangle\langle(\xi_{1})_{n},(\xi^{(t)})_{n}\rangle\right)a_{1}.$

Then it follows that

$\begin{split}\left\|\tilde{a}^{(t+1)}-\bar{a}^{(t+1)}\right\|_{2}&\leq\left\|\sum_{m=2}^{r}\lambda_{m}\langle
b_{m},b^{(t)}\rangle\left\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\right\rangle_{\mathcal{L}^{2}}a_{m}\right\|_{2}+\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k}){\mathbfcal{Z}}_{\cdot
js_{k}}\right\|_{2}\\\
&\qquad+\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k})\varepsilon_{\cdot
jk}\right\|_{2}.\end{split}$ (46)

We bound the three terms in the right-hand-side of (46) separately.

  * •

$\left\|\sum_{m=2}^{r}\lambda_{m}\langle
b_{m},b^{(t)}\rangle\left\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\right\rangle_{\mathcal{L}^{2}}a_{m}\right\|_{2}$.
By condition (23), for any $m,m^{\prime}=1,\ldots,r$, $m\neq m^{\prime}$, we
have

$\begin{split}\left|\frac{1}{n}\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\rangle-\langle\xi_{m},\xi^{(t)}\rangle_{\mathcal{L}^{2}}\right|&\leq
C\left(\zeta_{n}\left\|\xi_{m}\xi^{(t)}\right\|_{\mathcal{L}^{2}}+\zeta_{n}^{2}\left\|\xi_{m}\xi^{(t)}\right\|_{\mathcal{H}}\right)\leq
C\zeta_{n},\\\
\left|\frac{1}{n}\langle(\xi_{m})_{n},(\xi_{m^{\prime}})_{n}\rangle-\langle\xi_{m},\xi_{m^{\prime}}\rangle_{\mathcal{L}^{2}}\right|&\leq
C\left(\zeta_{n}\left\|\xi_{m}\xi_{m^{\prime}}\right\|_{\mathcal{L}^{2}}+\zeta_{n}^{2}\left\|\xi_{m}\xi_{m^{\prime}}\right\|_{\mathcal{H}}\right)\leq
C\zeta_{n}.\end{split}$

Here we use the assumptions that
$\|\xi_{m}\|_{\mathcal{L}^{2}}=\|\xi^{(t)}\|_{\mathcal{L}^{2}}=1$,
$\|\xi_{m}\|_{\mathcal{H}}\leq C_{\xi},\|\xi^{(t)}\|_{\mathcal{H}}\leq 2\kappa
C_{\xi}^{\prime}$ and the condition $\zeta_{n}\leq c/\kappa$. Then it follows
that

$\begin{split}\max_{m\neq
m^{\prime}}\frac{1}{n}\left|\langle(\xi_{m})_{n},(\xi_{m^{\prime}})_{n}\rangle\right|&\leq\max_{l\neq
l^{\prime}}\left|\langle\xi_{l},\xi_{l^{\prime}}\rangle_{\mathcal{L}^{2}}\right|+C\zeta_{n}\leq\mu+C\zeta_{n},\\\
\frac{1}{n}\langle(\xi_{1})_{n},(\xi^{(t)})_{n}\rangle&\geq\langle\xi_{1},\xi^{(t)}\rangle_{\mathcal{L}^{2}}-C\zeta_{n}\geq\sqrt{1-(\delta^{(t)})^{2}}-C\zeta_{n}\geq
1/2,\end{split}$ (47)

where the last inequality comes from the assumptions on $\delta^{(t)}$ and
$\zeta_{n}$. Therefore,

$\begin{split}\frac{1}{n}\left|\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\rangle\right|&\leq\langle\xi_{m},\xi^{(t)}\rangle_{\mathcal{L}^{2}}+C\zeta_{n}\leq\left|\langle\xi_{m},\xi_{1}\rangle_{\mathcal{L}^{2}}\right|+\left\|\xi_{1}-\xi^{(t)}\right\|_{\mathcal{L}^{2}}+C\zeta_{n}\\\
&\overset{\eqref{ineq:xi-N-
incoherence}}{\leq}\mu+2\delta+C\zeta_{n};\end{split}$

similarly, one can prove

$\langle b_{m},b^{(t)}\rangle\leq\mu+2\delta.$

Note that here we do not have the term $\zeta_{n}$ since $b^{(t)}$ corresponds
to a tabular mode and has error term $\zeta_{n}$ from discretization.

Combining the above results, we have

$\begin{split}&\left\|\sum_{m=2}^{r}\lambda_{m}\langle
b_{m},b^{(t)}\rangle\cdot\left\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\right\rangle
a_{m}\right\|_{2}^{2}\\\ \leq&\sum_{m=2}^{r}\lambda_{m}^{2}\langle
b_{m},b^{(t)}\rangle^{2}\left\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\right\rangle^{2}\\\
&\qquad+\sum_{\begin{subarray}{c}s,s^{\prime}\in[r]/\\{1\\}\\\ m\neq
m^{\prime}\end{subarray}}\lambda_{m}\lambda_{m^{\prime}}\langle
b_{m},b^{(t)}\rangle\langle
b_{m^{\prime}},b^{(t)}\rangle\left\langle(\xi_{m})_{n},(\xi^{(t)})_{n}\right\rangle\left\langle(\xi_{m^{\prime}})_{n},(\xi^{(t)})_{n}\right\rangle\langle
a_{m},a_{m^{\prime}}\rangle\\\ &\leq
Cn^{2}(r-1)\lambda_{max}^{2}(\mu+\delta)^{2}(\mu+\delta+\zeta_{n})^{2}\\\
&\qquad+Cn^{2}(r-1)^{2}\lambda_{max}^{2}\mu^{2}(\mu+\delta)^{2}(\mu+\delta+\zeta_{n})^{2}\\\
&\leq
Cn^{2}(r-1)\lambda_{max}^{2}(\mu+\delta)^{2}(\mu+\delta+\zeta_{n})^{2}.\end{split}$
(48)

Here we use the assumption $\mu\leq r^{-1/2}$ to obtain the final inequality.

  * •

$\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k}){\mathbfcal{Z}}_{\cdot
js_{k}}\right\|_{2}$. First, let $a^{*}\in{\mathbb{S}}^{p_{1}-1}$ such that

$\frac{1}{n}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}a^{*}_{i}(b^{(t)})_{j}\xi^{(t)}(s_{k}){\mathbfcal{Z}}_{ijs_{k}}=\left\|\frac{1}{n}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b_{1}^{(t)})_{j}\xi^{(t)}_{1}(s_{k}){\mathbfcal{Z}}_{\cdot
js_{k}}\right\|_{2}.$

Then following the same argument as (35), we have

$\begin{split}&\frac{1}{n}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}a^{*}_{i}(b^{(t)})_{j}\xi^{(t)}(s_{k}){\mathbfcal{Z}}_{ijs_{k}}\\\
\leq&\left\|\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}a^{*}_{i}(b^{(t)})_{j}{\mathbfcal{Z}}_{ijs_{k}}\right\|_{\infty}\cdot\left(\frac{1}{n}\sum_{k=1}^{n}\left|\xi^{(t)}(s_{k})\right|\right)\\\
\leq&\mathcal{E}\cdot\left\|\xi^{(t)}\right\|_{n}\\\
\overset{\eqref{ineq:A-2}}{\leq}&\mathcal{E}\cdot\sqrt{\|\xi^{(t)}\|_{\mathcal{L}^{2}}^{2}+C_{2}\zeta_{n}\|\xi^{(t)}\|_{\mathcal{L}^{2}}\cdot\|\xi^{(t)}\|_{\mathcal{H}}+C_{2}\zeta_{n}^{2}\|\xi^{(t)}\|_{\mathcal{H}}^{2}}\\\
\leq&\mathcal{E}\cdot\sqrt{1+C_{2}\zeta_{n}\|\xi^{(t)}\|_{\mathcal{H}}+C_{2}\zeta_{n}^{2}\|\xi^{(t)}\|_{\mathcal{H}}}\leq
2\mathcal{E}.\end{split}$

Note that the last inequality comes from the assumption that
$\|\xi_{1}^{(t)}\|_{\mathcal{H}}\leq 2\kappa C_{\xi}$ and $\zeta_{n}\leq
c/\kappa$ for a sufficiently small constant $c$.

Therefore,

$\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k}){\mathbfcal{Z}}_{\cdot
js_{k}}\right\|_{2}\leq 2n\mathcal{E}.$ (49)

  * •

$\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k})\varepsilon_{\cdot
jk}\right\|_{2}$. By Condition $(A_{3})$, we have

$\begin{split}\left\|\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}(b^{(t)})_{j}\xi^{(t)}(s_{k})\varepsilon_{\cdot
jk}\right\|_{2}&\leq\sup_{a\in{\mathbb{S}}^{p_{1}-1},b\in{\mathbb{S}}^{p_{2}-1}}\frac{1}{n}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{k=1}^{n}a_{i}b_{j}\xi^{(t)}(s_{k})\varepsilon_{ijk}\\\
&\overset{\eqref{ineq:A-3}}{\leq}Cn\tau\left\\{\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\|\xi^{(t)}\|_{\mathcal{L}^{2}}+\left(\zeta_{n}^{2}+\frac{(p_{1}+p_{2})\log
n}{n}\right)\|\xi^{(t)}\|_{\mathcal{H}}\right\\}\\\ &\leq
Cn\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right).\end{split}$ (50)

Combining (46), (48), (49) and (50), we obtain

$\left\|\tilde{a}^{(t+1)}-\bar{a}^{(t+1)}\right\|_{2}\leq
Cn\left\\{\sqrt{r-1}\lambda_{max}(\mu+\delta)(\mu+\delta+\zeta_{n})+\mathcal{E}+\tau\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}.$
(51)

Also note that

$\left\|\bar{a}_{1}^{(t+1)}\right\|_{2}=\lambda_{1}\left|\langle
b_{1},b^{(t)}_{1}\rangle\right|\cdot\left|\langle(\xi_{1})_{n},(\xi_{1}^{(t)})_{n}\rangle\right|\overset{\eqref{ineq:xi-
N-incoherence}}{\geq}\frac{n\lambda_{min}}{4}.$

Therefore, by Lemma 7, we obtain

${\rm dist}(a_{1},a_{1}^{(t+1)})\leq
C\sqrt{r-1}\kappa(\mu+\delta)(\mu+\delta+\zeta_{n})+C\left\\{\frac{\mathcal{E}}{\lambda_{min}}+\frac{\tau}{\lambda_{min}}\left(\zeta_{n}+\sqrt{\frac{p_{1}+p_{2}}{n}}\right)\right\\}.$

∎

###  A.2 Additional Proofs

###### Proof of Theorem 2.

We only need to prove the bound of ${\rm dist}(a,a^{(0)})$ as the proof for
the other tabular mode essentially follows. Denote
$\widetilde{\boldsymbol{Y}}:=\mathcal{M}_{1}(\widetilde{\mathbfcal{Y}})$ and
$\boldsymbol{Z}:=\mathcal{M}_{1}({\mathbfcal{Z}})$. Then we can write

$\widetilde{\boldsymbol{Y}}=\lambda
a(b\otimes\xi_{n})^{\top}+\boldsymbol{Z}+\boldsymbol{Z}_{\varepsilon},$

where $\boldsymbol{Z}_{\varepsilon}\in{\mathbb{R}}^{p_{1}\times p_{2}n}$ and
has i.i.d. $N(0,\tau^{2})$ entries.

Note that $a^{(0)}$ is the leading eigenvector of
$\widetilde{\boldsymbol{Y}}\widetilde{\boldsymbol{Y}}^{\top}$, which can be
decomposed as:

$\begin{split}\widetilde{\boldsymbol{Y}}\widetilde{\boldsymbol{Y}}^{\top}&=\lambda^{2}\cdot\|\xi_{n}\|_{2}^{2}\cdot
aa^{\top}+\boldsymbol{Z}\boldsymbol{Z}^{\top}+\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}+\boldsymbol{Z}\boldsymbol{Z}_{\varepsilon}^{\top}+\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}^{\top}\\\
&\qquad+\lambda
a(b\otimes\xi_{n})^{\top}(\boldsymbol{Z}+\boldsymbol{Z}_{\varepsilon})^{\top}+\lambda(\boldsymbol{Z}+\boldsymbol{Z}_{\varepsilon})(b\otimes\xi_{n})a^{\top}.\end{split}$

Then,

$\begin{split}&\left\|\left(\widetilde{\boldsymbol{Y}}\widetilde{\boldsymbol{Y}}^{\top}-{\mathbb{E}}\left[\boldsymbol{Z}\boldsymbol{Z}^{\top}\bigg{|}\\{s_{k}\\}_{k=1}^{n}\right]-{\mathbb{E}}\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}\right)-\lambda^{2}\|\xi_{n}\|^{2}\cdot
aa^{\top}\right\|\\\
\leq&\left\|\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}-{\mathbb{E}}\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|+2\lambda\left\|a(b\otimes\xi_{n})\boldsymbol{Z}_{\varepsilon}^{\top}\right\|\\\
&+\left\|\boldsymbol{Z}\boldsymbol{Z}^{\top}-{\mathbb{E}}\left[\boldsymbol{Z}\boldsymbol{Z}^{\top}\bigg{|}\\{s_{k}\\}_{k=1}^{n}\right]\right\|+2\lambda\left\|a(b\otimes\xi_{n})\boldsymbol{Z}^{\top}\right\|+2\left\|\boldsymbol{Z}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|.\end{split}$
(52)

Before we proceed, we first provide the probabilistic upper and lower bounds
for $\|\xi\|_{n}^{2}=n^{-1}\|\xi_{n}\|_{2}^{2}$. Recall that
$\|\xi\|_{\mathcal{L}^{2}}=1$ and $\|\xi\|_{\mathcal{H}}\leq C_{\xi}$. By
Lemma 5, with probability at least $1-n^{-9}$,

$\begin{split}\|\xi\|_{n}^{2}&\leq
1+C\zeta_{n}\leq\frac{3}{2},\quad\|\xi\|_{n}^{2}\geq
1-C\zeta_{n}\geq\frac{1}{2}.\\\ \end{split}$ (53)

Now we bound the five terms in (52) separately. First, by the tail bound of
Wishart-type random matrix (Lemma 9), we have

${\mathbb{P}}\left(\left\|\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}-{\mathbb{E}}\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|\geq
C\tau^{2}\left((\sqrt{p_{1}}+\sqrt{p_{2}n}+x)^{2}-p_{2}n\right)\right)\leq
e^{-x}.$

Let $x=\sqrt{p_{1}}$. Then we have with probability at least
$1-e^{-\sqrt{p_{1}}}$,

$\left\|\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}-{\mathbb{E}}\boldsymbol{Z}_{\varepsilon}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|\leq
C\tau^{2}\left(p_{1}+\sqrt{p_{1}p_{2}n}\right).$ (54)

Note that conditioning on $\\{s_{k}\\}_{k=1}^{n}$,
$\left\|a(b\otimes\xi_{n})^{\top}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|=\|(b\otimes\xi_{n})^{\top}\boldsymbol{Z}_{\varepsilon}^{\top}\|_{2}$
is the $l_{2}$ norm for a $p_{1}$-dimensional random vector with i.i.d.
$N(0,\tau^{2}\|\xi_{n}\|^{2})$ entries. Therefore, by Gaussian concentration,
we have with probability at least $1-e^{-p_{1}}$

$\left\|a(b\otimes\xi_{n})^{\top}\boldsymbol{Z}_{\varepsilon}^{\top}\right\|\leq
C\tau\|\xi_{n}\|_{2}\sqrt{p_{1}}\overset{\eqref{ineq:xi-n-
bound}}{\leq}C\tau\sqrt{np_{1}}.$ (55)

The analysis for
$\left\|\boldsymbol{Z}\boldsymbol{Z}^{\top}-{\mathbb{E}}\left[\boldsymbol{Z}\boldsymbol{Z}^{\top}\bigg{|}\\{s_{k}\\}_{k=1}^{n}\right]\right\|$
and $\|a(b\otimes\xi_{n})\boldsymbol{Z}^{\top}\|$ are more involved as the
entries of $\boldsymbol{Z}$ are dependent. We conduct our analysis
conditioning on fixed values of $\\{s_{k}\\}_{k=1}^{n}$. Note that for any
fixed pair $(i,j)\in[p_{1}]\times[p_{2}]$,

$\left({\mathbfcal{Z}}_{ijs_{1}},{\mathbfcal{Z}}_{ijs_{2}},\ldots,{\mathbfcal{Z}}_{ijs_{n}}\right)^{\top}\overset{i.i.d.}{\sim}N(0,\boldsymbol{\Sigma}),$

where
$(\boldsymbol{\Sigma})_{ll}=\text{Var}({\mathbfcal{Z}}_{ijs_{l}})\leq\sigma^{2}$
and it follows that

$\left\|\boldsymbol{\Sigma}\right\|\leq\left\|\boldsymbol{\Sigma}\right\|_{*}\leq
n\sigma^{2},$ (56)

where $\|\cdot\|$ and $\|\cdot\|_{*}$ are the matrix spectral and nuclear
norms respectively. Let
$\boldsymbol{\Sigma}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{\top}$
be the eigenvalue decomposition of $\boldsymbol{\Sigma}$ with
$\boldsymbol{U}\in{\mathbb{O}}_{n}$ and $\boldsymbol{\Lambda}$ being a
diagonal matrix. Let $\boldsymbol{W}$ be a $p_{1}$-by-$(np_{2})$ random matrix
with i.i.d. $N(0,1)$ entries. We can then equivalently represent

$\boldsymbol{Z}=\boldsymbol{W}\cdot\begin{bmatrix}\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}&&&\\\
&\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}&&\\\ &&\vdots&\\\
&&&\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}\end{bmatrix},$ (57)

and it follows that

$\boldsymbol{Z}\boldsymbol{Z}^{\top}=\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}$
(58)

for $\tilde{\boldsymbol{W}}\in{\mathbb{R}}^{p_{1}\times(np_{2})}$ with
independent mean-zero Gaussian random variables and
$\text{Var}(\boldsymbol{W}_{ij})=\boldsymbol{\Lambda}_{k+1,k+1}$ iff $(j\text{
mod }n)=k$. Now applying the Concentration of Heteroskedastic Wishart-type
Matrix (Lemma 9), we have

${\mathbb{P}}\left(\left\|\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}-{\mathbb{E}}\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}\right\|\geq
C\left((\sigma_{R}+\sigma_{C}+\sigma_{*}\log
p_{1}+\sigma_{*}x)^{2}-\sigma_{R}^{2}\right)\right)\leq e^{-x},$

where
$\sigma_{R}^{2}=p_{2}\sum_{k=1}^{n}\boldsymbol{\Lambda}_{kk}=p_{2}\|\boldsymbol{\Sigma}\|_{*}$,
$\sigma_{C}^{2}=p_{1}\|\boldsymbol{\Sigma}\|$ and
$\sigma_{*}=\|\boldsymbol{\Sigma}\|$. By specifying $x=\sqrt{p_{1}}$ and use
the bounds (56), we obtain that with probability at least
$1-n^{-9}-e^{-\sqrt{p_{1}}}$,

$\begin{split}\left\|\boldsymbol{Z}\boldsymbol{Z}^{\top}-{\mathbb{E}}\left[\boldsymbol{Z}\boldsymbol{Z}^{\top}\bigg{|}\\{s_{k}\\}_{k=1}^{n}\right]\right\|&=\left\|\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}-{\mathbb{E}}\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}\right\|\\\
&\leq
p_{1}\|\boldsymbol{\Sigma}\|+\sqrt{p_{1}p_{2}\|\boldsymbol{\Sigma}\|\cdot\|\boldsymbol{\Sigma}\|_{*}}\\\
&\leq n\sigma^{2}\left(p_{1}+\sqrt{p_{1}p_{2}}\right).\end{split}$ (59)

To bound $\left\|a(b\otimes\xi_{n})\boldsymbol{Z}^{\top}\right\|$, we also
utilize the representation (57) and it follows by the similar argument as (55)
that

$\begin{split}\left\|a(b\otimes\xi_{n})^{\top}\boldsymbol{Z}^{\top}\right\|&\leq
C\left\|(b\otimes\xi_{n})^{\top}{\rm
diag}(\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top},\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top},\ldots,\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top})\right\|_{2}\sqrt{p_{1}}\\\
&\leq
C\|\xi_{n}\|_{2}\cdot\|\Sigma\|^{1/2}\cdot\sqrt{p_{1}}\overset{\eqref{ineq:xi-
n-bound},\eqref{ineq:discrete-cov-spectral-
norm}}{\leq}n\sqrt{p_{1}}.\end{split}$ (60)

Finally, note that $\boldsymbol{Z}^{\top}\boldsymbol{Z}_{\varepsilon}$ has the