# Dynamic Matrix Recovery

Ziyuan Chen
School of Mathematical Sciences, Center for Statistical Science,
Peking University, Beijing, China
Ying Yang
Academy of Mathematics and Systems Science,
Chinese Academy of Sciences, Beijing, China
Fang Yao
School of Mathematical Sciences, Center for Statistical Science,
Peking University, Beijing, China
Ziyuan Chen and Ying Yang are the co-first authors, and contribute equally to
this work. Fang Yao is the corresponding author, E-mail: fyaomath.pku.edu.cn.
Fang Yao’s research is partially supported by the National Key R&D Program of
China (No. 2022YFA1003801, 2020YFE0204200), the National Natural Science
Foundation of China (No. 12292981, 11931001), the LMAM and the Fundamental
Research Funds for the Central Universities and the LMEQF. Ying Yang’s
research is partially supported by China Postdoctoral Science Foundation (No.
2022TQ0360, 2022M723334), the National Natural Science Foundation of China
(No. 71988101) and the Guozhi Xu Posdoctoral Research Foundation.

###### Abstract

Matrix recovery from sparse observations is an extensively studied topic
emerging in various applications, such as recommendation system and signal
processing, which includes the matrix completion and compressed sensing models
as special cases. In this work, we propose a general framework for dynamic
matrix recovery of low-rank matrices that evolve smoothly over time. We start
from the setting that the observations are independent across time, then
extend to the setting that both the design matrix and noise possess certain
temporal correlation via modified concentration inequalities. By pooling
neighboring observations, we obtain sharp estimation error bounds of both
settings, showing the influence of the underlying smoothness, the dependence
and effective samples. We propose a dynamic fast iterative shrinkage-
thresholding algorithm that is computationally efficient, and characterize the
interplay between algorithmic and statistical convergence. Simulated and real
data examples are provided to support such findings.

Keywords: compressed sensing, local smoothing, low rank, matrix completion

##  1 Introduction

Matrix recovery from sparse observations has been extensively studied in a
wide range of problems, such as signal processing (Li et al., 2019),
recommendation system (Koren et al., 2009) and the quantum state tomography
(Gross et al., 2010). The classical model is to link the response $Y$ with the
random matrix $X\in\mathbb{R}^{m_{1}\times m_{2}}$ via

$Y=\operatorname{Tr}(X^{\top}M^{0})+\xi$ (1)

to recover the low-rank matrix $M^{0}\in\mathbb{R}^{m_{1}\times m_{2}}$ from
the independent observations $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$, where
$\operatorname{Tr}(\cdot)$ is the trace of a matrix, $\xi$ is a mean zero
random error and $n$ denotes the sample size, which is typically much smaller
than $m_{1}m_{2}$ (i.e., $n\ll m_{1}m_{2}$). By imposing different assumptions
on $X$, the model includes a variety of matrix recovery models, in which
matrix completion and compressed sensing are two important cases. Most
existing works assume that the target matrix $M^{0}$ is static (Klopp, 2011;
Koltchinskii & Xia, 2015; Fan et al., 2021; Xia & Yuan, 2021). However, there
are various scenarios where $M^{0}$ is dynamic and observed at multiple time
points, such as in time-varying user interests in recommendation systems (Xu
et al., 2020), temporally-changed graph models (Qiu et al., 2017) and dynamic
quantum dot matrices (Csurgay et al., 2000). To formulate a general framework
accounting the dynamic changes in numerous applications, we propose
investigating the recovery of the dynamic matrix as follows,

$Y_{t}=\operatorname{Tr}(X_{t}^{\top}M_{t}^{0})+\xi_{t},\ t=1,\ldots,T,$ (2)

where $t$ is the timestamp, $\\{Y_{t},X_{t},\xi_{t}\\}$ is the tuple of the
response, the design matrix and the noise at time $t$, respectively, and
$M_{t}^{0}$ is a dynamic matrix which is assumed smooth across time $t$ with a
low-rank structure.

A prevalent line of research for solving the classical static model (1) is to
impose strict low-rank constraints based on matrix factorization. For
instance, Keshavan et al. (2009) proposed singular value decomposition to make
the recovery, Zhao et al. (2015) used QR decomposition instead, and Zheng &
Lafferty (2016) adopted Burer-Monteiro factorization for rectangle matrix
completion. Another stream is the convex approaches, thanks to Candès & Recht
(2009) that first relaxed the non-convex rank restriction to nuclear norm for
efficient optimization. A variety of convex algorithms and theoretical results
are further established (Recht et al., 2010; Candes & Plan, 2011; Agarwal et
al., 2012). Though $M_{t}^{0}$ in (2) can be recovered by adopting the static
matrix recovery at each time $t$, these single-stage estimates are less
desirable, as they only use data at time $t$ to reconstruct $M_{t}^{0}$ and
fail to take advantage of the neighboring observations and smoothness of the
underlying matrix. Related work also includes tensor regression to rebuild
$M_{t}^{0}$, which views the dynamic model (2) from a static perspective by
letting $N\in\mathbb{R}^{T\times m_{1}\times m_{2}}$ be a three order tensor
whose $(t,i,j)$th element is the $(i,j)$th element of $M_{t}^{0}$. Gandy et
al. (2011) and Liu et al. (2012) proposed algorithms for tensor completion and
Zhang et al. (2020) considered the low-rank tensor trace regression via the
importance sketching algorithm. We mention that their low-rank assumption of
$N$ is different from model (2) that $M_{t}^{0}$ is dynamically and smoothly
evolving along time $t$. Consequently, one needs to process the whole data to
recover $M_{t}$ in tensor regression setting, which is time- and space-
consuming.

A more efficient approach is to adopt the dynamic modeling/algorithm that
avoids computing the whole data for update and has been adopted in some
relevant problems. Lois & Vaswani (2015) designed a memory-saving algorithm
for batch robust principle component analysis; Gao & Wang (2018) studied
dynamic robust principal component analysis and proposed a fast iterative
algorithm; Hu & Yao (2022) proposed a unified framework to directly estimate
dynamic principal subspace in high dimension. The dynamic algorithm for matrix
recovery has been explored in Xu & Davenport (2016), which assumed that
$\xi_{t}$ are independently and identically distributed (i.i.d.) zero-mean
Gaussian error and imposed a parametric assumption on structure of the
underlying matrix, $M_{t}^{0}=M_{t-1}^{0}+\varepsilon_{t}$, where
$\varepsilon_{t}$ are i.i.d. zero-mean Gaussian noises. Another line of work,
namely the dynamic tensors, also finds extensive applications in various
domains, and several works have been noted in this direction. For instance,
Zhang et al. (2021) proposed a time-varying coefficient model for temporal
tensor factorization, employing a polynomial spline approximation, and Zhou et
al. (2013) introduced a family of rank-R generalized linear tensor regression
models. Most of the current research on dynamic tensors relies on tensor
decomposition techniques and assumes that only certain decomposition factors
or coefficients are changing (Zhang et al., 2021; Bi et al., 2021; Koren,
2009; Wang et al., 2016).

Different from the existing work, we consider the general scenario where the
dynamics are not limited to decomposition components or a parametric structure
priori. There exist two major challenges. First, one needs to maintain the
low-rank structure in the whole procedure of dynamic algorithm. Second, the
combination of neighboring information often inflates the estimation error
when the observations are correlated along time. To resolve these, we utilize
local smoothing which has the advantages that the manifold structure of low-
dimensional matrices can be approximately retained and temporal dependence has
less influence when $T$ increases. Based on these considerations, we propose a
dynamic matrix recovery algorithm without imposing structural assumption
besides smoothness and allowing the measurement error $\xi_{t}$ and/or design
matrix $X_{t}$ to be possibly dependent across time. The idea is to pool
together the observations that lie in the local window to estimate
$M_{t}^{0}$. Noting that $M_{t}^{0}$ is low-rank, we set the target function
composed of a locally weighted loss function and the nuclear norm penalty of
$M_{t}^{0}$ and propose a Dynamic Fast Iterative Shrinkage-Thresholding
Algorithm (DFISTA) for optimization.

The main contributions of this work are summarized as follows. First, we
establish a general framework for dynamic matrix recovery with theoretical
guarantees. We allow the underlying matrix varying smoothly under the low-rank
constraint and allow the observations to be dependent across time, which has
more flexibility and applicability. While existing works only study the
estimation error bounds for specific problems such as matrix completion and
compressed sensing (Koltchinskii et al., 2011a; Negahban & Wainwright, 2011;
Xu & Davenport, 2016; Fan et al., 2021), we derive the error bounds under this
general framework, which can be readily adopted to such matrix recovery
problems. Second, the proposed method attains a faster convergence rate than
classical static methods in the following sense. Take the matrix completion
for instance. The ratio of error bounds of the proposed estimator compared to
the classical one (i.e., single-stage estimates at each $t$) is
$O_{p}\big{(}\\{m\log m/n_{t}\\}^{-1/10}T^{-2/5}\big{)}$, where $n_{t}$ is the
number of observations at time $t$ and $m=m_{1}\vee m_{2}$ (”$\vee$” indicates
the greater of two quantities). When the observations are sparse in the sense
that $n_{t}\ll m\log m$, the ratio becomes $o_{p}\big{(}T^{-2/5}\big{)}$.
Specifically, the dynamic method uses only $n^{\prime}=\big{(}m\log
m/n\big{)}^{-1/4}n/T=o(n/T)$ observations to attain the same convergence rate
as the static method using $n$ observations. Last but not least, a
computationally efficient algorithm is devised for dynamic matrix recovery.
This algorithm iterates only a subset of the data to update the estimate at
each step. By using the estimate of the last step as the initial of the
current optimization procedure, the algorithm converges faster than the
single-stage recovery. Specifically, the ratio of computational complexities
of our algorithm compared to the single-stage recovery is $O(T^{-3/5}\vee
n^{-1/2})$.

The rest of the paper is organized as follows. In Section 2, we present the
proposed method and algorithm for dynamic matrix recovery, and establish
theoretical guarantees in Section 3 for both independent and dependent
observations with algorithmic analysis. Finally, we conduct simulation and
real data experiments in Section 4 and 5 to show the favorable performance of
the proposed dynamic method/algorithm compared with other existing methods.
Additional theoretical results and technical proofs are deferred to an online
Supplementary Material due for space economy.

##  2 Methodology and Algorithm

We consider the general framework in the form of dynamic trace regression
whose coefficient matrix is changing smoothly along time. Specifically, at
time $t=1,\ldots,T$, we observe $(Y_{ti},X_{ti})$ following model (2),

$Y_{ti}=\operatorname{Tr}(X_{ti}^{\top}M_{t}^{0})+\xi_{ti},\quad
i=1,\ldots,n_{t},$

where $M_{t}^{0}$ is the time-varying coefficient matrix, $X_{ti}$ and
$\xi_{ti}$ are design matrices and zero-mean noises, respectively,
$\\{X_{ti},\xi_{ti}\\}_{i}$ (given $t$) are mutually independent across
subject and $\\{X_{ti},\xi_{ti}\\}_{t}$ (given $i$) may be correlated across
time. Parallel to the static model (1), the dynamic trace regression includes
several problems such as dynamic matrix completion and dynamic compressed
sensing that can be applied to dynamic recommendation systems and video signal
processing as illustrated in numerical experiments of Section 5.

###### Example 1

(Dynamic matrix completion). For each $t$, suppose that the design matrices
$X_{ti}$ are i.i.d. uniformly distributed on the set

$\mathcal{E}=\\{e_{j}(m_{1})e_{k}^{\top}(m_{2}):1\leq j\leq m_{1},1\leq k\leq
m_{2}\\},$ (3)

where $e_{j}(m)$ are the canonical basis vectors in $\mathbb{R}^{m}$ and then
$\mathcal{E}$ forms an orthonormal basis in the space $\mathbb{R}^{m_{1}\times
m_{2}}$. Then estimating $M_{t}^{0}$ is equivalent to the problem of matrix
completion under uniform sampling at random (USR).

###### Example 2

(Dynamic compressed sensing). For each $t$, suppose that the design matrices
$X_{ti}$ are i.i.d. replicates of a random matrix $X_{t}$ such that $\langle
M,X_{t}\rangle$ is a sub-gaussian random variable for any
$M\in\mathbb{R}^{m_{1}\times m_{2}}$. In this work, we focus on the concrete
example that each design matrix $X_{ti}$ is a random matrix whose elements are
independent mean-zero sub-gaussian random variables with variance
$\sigma_{X_{t}}^{2}$. Moreover, if the design $X_{t}$ is homogeneous across
$t$, we denote the variance as $\sigma_{X}^{2}$.

We first briefly review the classical static method in the context of dynamic
modeling, i.e., adopting the static method at each time $t$, which is referred
to as the single-stage estimates in the sequel. Define the matrix inner
product by $\langle M,X\rangle=\operatorname{Tr}(M^{\top}X)$. Koltchinskii et
al. (2011a) proposed to incorporate the knowledge of the distribution of
$X_{ti}$ and estimate $M_{t}^{0}$ by

$\widehat{M}_{t}^{\lambda}=\mathop{\rm
argmin}_{M\in\mathbb{M}}\left[\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\mathbb{E}(\langle
M,X_{ti}\rangle)^{2}-\left\langle\frac{2}{n_{t}}\sum_{i=1}^{n_{t}}Y_{ti}X_{ti},M\right\rangle+\lambda\|M\|_{1}\right],$
(4)

where $\mathbb{M}\subset\mathbb{R}^{m_{1}\times m_{2}}$ is a convex set of
matrices, $\lambda>0$ is a regularization parameter and $\|M\|_{1}$ is the
nuclear norm of $M$. They proved that in the setting of matrix completion in
Example 1, if $\lambda$ is appropriately set, with probability at least
$1-3/(m_{1}+m_{2})$,

$(m_{1}m_{2})^{-1/2}\|\widehat{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}\leq
C\left\\{\frac{\log(m_{1}+m_{2})\max(m_{1},m_{2})r_{t}}{n_{t}}\right\\}^{1/2},$

where $r_{t}=\mathop{\rm rank}(M_{t}^{0})$ and $\|M\|_{2}$ is the Frobenius
norm of matrix $M$. Though they have proved that the rate is optimal up to
logarithmic factors comparing to the minimax lower bounds in classical trace
regression, this approach fails to incorporate temporal information and hence
is not efficient in multi-stage case.

A direct modification is to utilize the two-step procedure that first obtains
the single-stage estimates $\widehat{M}^{\lambda}_{t}$ and then applies local
smoothing to $\\{\widehat{M}_{t}^{\lambda}\\}_{t}$ to achieve the final
estimate as follows,

$\widehat{M}_{t}^{\text{t-s}}=\sum_{j=1}^{T}\omega_{h}(j-t)\widehat{M}_{j}^{\lambda},$
(5)

where the local smoothing weights
$\omega_{h}(j-t)=K_{Th}\left(j-t\right)/\sum_{k=1}^{T}K_{Th}\left(k-t\right)$
with the kernel function $K(x)$ and bandwidth $h$. The two-step estimation
approach computes a weighted average of low-rank matrices that have been
estimated at distinct time points. However, this method encounters a
significant challenge when dealing with sparse observations. In such
scenarios, the estimated matrix for each time point could potentially disrupt
the true low-rank structure, displaying significant deviation from the actual
value. Unfortunately, this disparity cannot be corrected through smoothing
techniques.

Motivated by the observation that the manifold structure of low-rank matrices
can be approximately preserved locally without matrix factorization, we design
the dynamic target function with nuclear norm penalty, i.e.,

$\widetilde{M}_{t}^{\lambda}=\mathop{\rm
argmin}_{M\in\mathbb{M}}\sum_{j=1}^{T}\omega_{h}(j-t)\left[\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}(\left\langle
M,X_{ji}\right\rangle)^{2}-2\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},M\right\rangle\right]+\lambda\|M\|_{1},$
(6)

which combines the observations across time and employs local smoothing for
incorporating the temporal information. It is worth noting that in many
applications including our Examples 1 and 2, the distribution of $X_{t}$ is
known, and yet this information has not been exploited in empirical risk. Thus
we may utilize the distribution knowledge by substituting the empirical form
with true values $\mathbb{E}\langle M,X_{ji}\rangle$. The temporal smoothness
of $M_{t}^{0}$ guarantees that information in adjacent time points can be
utilized to reduce the estimation error. We mention that the classical
estimate $\widehat{M}_{t}^{\lambda}$ in (4) can be regarded as a special case
of the proposed local estimate $\widetilde{M}_{t}^{\lambda}$ in (6) as the
bandwidth $h\rightarrow 0$. By choosing $h$ appropriately, our estimate pools
more observations to reduce the estimation variance while controlling the bias
caused by the varying of $M_{t}^{0}$, as discussed in detail in Section 3.

For implementation, we adopt the Fast Iterative Shrinkage-Thresholding
Algorithm (FISTA) which is an accelerated proximal method proposed by Beck &
Teboulle (2009). It improves traditional proximal methods and terminates in
$\mathcal{O}(\varepsilon^{-1/2})$ iterations with a $\varepsilon$-optimal
solution. To handle with the trace norm penalty optimization problem (6), we
use the matrix form FISTA established by Ji & Ye (2009) and Toh & Yun (2010).
Define the objective function in (6) as $F_{t}(M)=f_{t}(M)+g(M)$, where

$\displaystyle
f_{t}(M)=\sum_{j=1}^{T}\omega_{h}(j-t)\left[\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}(\left\langle
M,X_{ji}\right\rangle)^{2}-2\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},M\right\rangle\right].$

We mention that the gradient function of $f_{t}$ is equipped with Lipschitz
constant $L_{f}$ presented in Section 3.4, and that $g(M)=\lambda\|M\|_{1}$ is
a continuous convex function. Then the matrix form FISTA can be applied to
solve the dynamic trace regression problem

$\widetilde{M}^{\lambda}_{t}=\arg\min_{M\in\mathbb{M}}f_{t}(M)+g(M).$

Note that $\widetilde{M}^{\lambda}_{t}$ need to be calculated for a sequence
$t=1,2,\cdots,T$. The algorithm is accelerated at time $t$ by utilizing the
estimate from the timestamp $(t-1)$. Specifically, suppose that at time $t$ we
apply the matrix form FISTA for $k_{t}$ iteration times to numerically
approximate the estimate
$\widetilde{M}_{t}^{\lambda}=\arg\min_{M\in\mathbb{M}}f_{t}(M)+g(M)$ denoted
by $M_{t}^{(k_{t})}$. Then input $M_{t}^{(k_{t})}$ as the initial matrix in
the procedure of estimating
$\widetilde{M}_{t+1}^{\lambda}=\arg\min_{M\in\mathbb{M}}f_{t+1}(M)+g(M)$ at
time $(t+1)$. Such initialization method accelerates the algorithm, as
demonstrated in Section 3.4. We summarize the implementations of the proposed
DFISTA in Algorithm 2. In the algorithm, “svd” means calculating the single
value decomposition components $U,D,V$ for a matrix
$G\in\mathbb{R}^{m_{1}\times m_{2}}$ and “$(\cdot)_{+}$” means the threshold
operator such as $(D-\gamma)_{+}=\operatorname{diag}((\lambda_{1}-\gamma)\vee
0,\cdots,(\lambda_{m_{1}}-\gamma)\vee 0)$ where
$D=\operatorname{diag}(\lambda_{1},\cdots,\lambda_{m_{1}})$. The learning rate
$s_{k}$ is updated according to Nesterov (1983).

Algorithm 1 Dynamic Low-Rank Matrix Recovery

0: $\sum_{t=1}^{T}n_{t}$ matrices $X_{ti}\in\mathbb{R}^{m_{1}\times m_{2}}$
and $\sum_{t=1}^{T}n_{t}$ scalars $Y_{ti}$, $t=1,2,\dots,T$,
$i=1,2,\cdots,n_{j}$, max iteration steps $K$ and torrent tor.

0: $T$ matrices $M_{t}\in\mathbb{R}^{m_{1}\times m_{2}},\ t=1,2,\dots,T$.

1: Choose a initial matrix $M_{0}\in\mathbb{R}^{m_{1}\times m_{2}}$.

2: while $t\leq T$ do

3: $M_{t}^{(-1)}\leftarrow M_{t-1},\ M_{t}^{(0)}\leftarrow M_{t-1}$;

4: $s_{0}\leftarrow 1,\ s_{-1}\leftarrow 1$;

5: for $k=0,1,\ldots,K$ do

6: $N_{t}^{(k)}\leftarrow
M_{t}^{(k)}+\frac{s_{k-1}-1}{s_{k}}\left(M_{t}^{(k)}-M_{t}^{(k-1)}\right)$;

7: $G_{t}^{(k)}\leftarrow N_{t}^{(k)}-L_{f}^{-1}\sum_{j=(t-Th/2)\vee
0}^{(t+Th/2)\wedge
T}\frac{w_{h}(j-t)}{n_{j}}\sum_{i=1}^{n_{j}}X_{ji}\left(\langle
X_{ji},N_{t}^{(k)}\rangle-Y_{j}\right)$;

8: $U,D,V\leftarrow$ svd($G_{t}^{(k)}$);

9: $M_{t}^{(k+1)}\leftarrow U(D-2\lambda/L_{f})_{+}V^{T}$;

10: $s_{k+1}\leftarrow\frac{1+\sqrt{1+4s_{k}^{2}}}{2}$;

11: if $|F_{t}(M_{t}^{(k+1)})-F_{t}(M_{t}^{(k)})|\leq{\tt tor}$ then

12: $\operatorname{break}$

13: end if

14: end for

15: $M_{t}=M_{t}^{(k+1)}$;

16: end while

17: return $M_{t},t=1,2,\cdots,T$

Unlike most factorization based algorithms that are challenging to analyze for
their convergence properties, we are able to analyze the error bounds of the
algorithm’s results with respect to the true value. This analysis helps bridge
the gap between the statistical estimator and the estimator obtained from the
algorithm, as discussed in detail in Section 3.4.

##  3 Theoretical Guarantees

In this section, we study both the statistical properties of the estimation
error $\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}$ and the algorithmic
convergence of the proposed method. We first present the error bound of the
estimate $\widetilde{M}_{t}^{\lambda}$ for general dynamic trace regressions.
Then we derive the explicit expressions of the upper bound of
$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}$ when $\\{X_{t}\\}_{t}$ and
$\\{\xi_{t}\\}_{t}$ are independent across $t$ and when $\\{X_{t}\\}_{t}$ and
$\\{\xi_{t}\\}_{t}$ are $\phi$-mixing processes, respectively. Noting that
$X_{t}$ and/or $\xi_{t}$ may be unbounded, the assumptions for classical
concentration inequalities are violated (Merlevède et al., 2009, 2011; Hang &
Steinwart, 2017). We thus modify the Talagrand’s inequalities and adopt
truncation technique appropriately to derive the convergence rate for the
dependent setting. In both scenarios, we shall explain how the number of time
points $T$, the smoothness of $M_{t}^{0}$ and the dependence of $X_{t}$ and
$\xi_{t}$ across time influence the estimation error explicitly.

In the sequel, we denote $\|M\|_{0}$, $\|M\|_{1}$, $\|M\|_{2}$,
$\|M\|_{\infty}$ and $\langle M,N\rangle$ as the rank, nuclear norm, Frobenius
norm, maximum singular value of $M$ and the trace of $M^{\top}N$,
respectively.

###  3.1 General error bounds

Let $\mathop{\rm vec}(X)\in\mathbb{R}^{m_{1}m_{2}}$ be the vectorization of
$X\in\mathbb{R}^{m_{1}\times m_{2}}$ and define the second moment matrix of
$X_{t}$ by $\Sigma_{t}\in\mathbb{S}^{m_{1}m_{2}\times m_{1}m_{2}}$ such that

$\Sigma_{t}=\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\mathbb{E}[\mathop{\rm
vec}({X}_{ti})\mathop{\rm vec}({X}_{ti})^{\top}].$

It is necesary to guarantee the identifiability of matrix recovery that the
positive semi-definite matrix $\Sigma_{t}$ is invertible. To simplify, we
consider the homogeneous scenario such that at each time point $t$, the design
$\\{X_{t}\\}_{t}$ have the same second moment matrix as stated below.

###### Assumption 1

There exists a positive definite matrix
$\Sigma\in\mathbb{S}_{+}^{m_{1}m_{2}\times m_{1}m_{2}}$ with smallest
eigenvalue $\mu>0$ such that $\Sigma_{j}=\Sigma,\ j=1,2,\dots,T$.

We remark that in the heterogeneous situation, the convergence rate can be
derived similarly, which is at the same order as the homogeneous setting with
a more tedious presentation. An important indication of Assumption 1 is that,
for all $M\in\mathbb{R}^{m_{1}\times m_{2}}$ and $j=1,2,\cdots,T$,

$\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}(\langle
M,X_{ji}\rangle)^{2}=\mathop{\rm vec}({M})^{\top}\Sigma\mathop{\rm
vec}({M})\geq\mu\|M\|_{2}^{2},$

which is a standard assumption in matrix recovery. The interpretation is that
the smallest eigenvalue $\mu>0$ determines the information contained in the
design matrix $X_{t}$. Intuitively, the larger $\mu$ is, the more information
is revealed by data, thus the more precise the estimator
$\widetilde{M}_{t}^{\lambda}$ is. In most cases, $\mu$ can be calculated
directly. For instances, in matrix completion of Example 1, the smallest
eigenvalue of $\Sigma$ is $\mu=1/m_{1}m_{2}$, and in compressed sensing of
Example 2, $\mu$ is the smallest variance of elements in $X_{ji}$, which is
equal to $\sigma_{X}^{2}$ under Assumption 1.

For brevity, we further denote
$\Delta_{j}=1/n_{j}\sum_{i=1}^{n_{j}}\left(Y_{ji}X_{ji}-\mathbb{E}(Y_{ji}X_{ji})\right)$
and

$\delta_{h}M_{t}=\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{h}(j-t)M_{j}^{0}\right\|_{2},\quad\mathcal{W}_{h}\Delta_{t}=\left\|\sum_{j=1}^{T}\omega_{h}(j-t)\Delta_{j}\right\|_{\infty}.$

The following theorem gives the general upper bound of estimation error
$\|\widetilde{M}_{t}^{\lambda}-M^{0}_{t}\|_{2}$.

###### Theorem 1

Under Assumption 1, if $\lambda\geq 2\mathcal{W}_{h}\Delta_{t}$, then

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\min\left\\{\left((\delta_{h}M_{t})^{2}+\frac{2\lambda}{\mu}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2},\frac{1+\sqrt{2}}{2}\frac{\lambda}{\mu}\sqrt{r_{t}}\right\\}.$
(7)

When selecting $\lambda=2\mathcal{W}_{h}\Delta_{t}$, we have

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\min\left\\{\left((\delta_{h}M_{t})^{2}+2\mu^{-1}\mathcal{W}_{h}\Delta_{t}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2},{(1+\sqrt{2})}{\mu^{-1}}{\mathcal{W}_{h}\Delta_{t}}\sqrt{r_{t}}\right\\}.$
(8)

The term $\delta_{h}M_{t}$ corresponds to the bias caused by kernel smoothing
which is affected by kernel $K(x)$ and the bandwidth $h$, and
$\mathcal{W}_{h}\Delta_{t}$ corresponds to the variance term caused by the
effective samples in the local window and the measurement error $\xi$. When
$h=0$, the result in Theorem 1 coincides with the classical result of
Corollary 1 in Koltchinskii et al. (2011a) such that when
$\lambda=2\|\Delta_{t}\|_{\infty}$,

$\left\|\widehat{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\min\left\\{\left(2\mu^{-1}\left\|\Delta_{t}\right\|_{\infty}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2},{(1+\sqrt{2})}{\mu^{-1}}{\|\Delta_{t}\|_{\infty}}\sqrt{r_{t}}\right\\}.$
(9)

Compared to (9), the bound (8) replaces the error term $\Delta_{t}$ by the
smoothed version $\mathcal{W}_{h}\Delta_{t}$ at the cost of the bias
$\delta_{h}M_{t}$.

To obtain an explicit expression of the error bound that may be used for
bandwidth selection, we further investigate the properties of
$\delta_{h}M_{t}$ and $\mathcal{W}_{h}\Delta_{t}$. Assumption 2 and 3 are
standard conditions for temporal smoothness and the kernel function.

###### Assumption 2

There exists a matrix function $M(t)$ supported on $[0,1]$, satisfying that

  * 1.

for any $j=1,2,\dots,T$, $M_{j}^{0}=M(j/T)$ and $\|M_{j}^{0}\|_{1}\leq C_{*},\
\max_{j,s,l}|[M_{j}^{0}]_{s,l}|\leq C_{M}$, for some constants $C_{*}$ and
$C_{M}$;

  * 2.

the first and second derivatives $\nabla M(t),\nabla^{2}M(t)$ are continuous
and bounded by some constants $D_{1}$ and $D_{2}$, respectively.

###### Assumption 3

The kernel function $K(x)$ is a symmetric probability density function on [-1,
1] and satisfies that

$\displaystyle\alpha(K)=\int x^{2}K(x)\ dx<\infty,\quad R(K)=\int K(x)^{2}\
dx<\infty.$ (10)

Define the $\psi-$norm $\|Z\|_{\psi(\eta)}$ for random variable $Z$ with
$0<\eta<\infty$ as

$\|Z\|_{\psi(\eta)}=\inf\left\\{t>0,\mathbb{E}\left[\exp\left(\frac{\|Z\|_{\infty}^{\eta}}{t^{\eta}}\right)\right]\leq
2\right\\}$

and
$\|Z\|_{\psi(\infty)}=\inf\left\\{t>0,\mathbb{E}\mathbf{1}_{\left\\{\|Z\|_{\infty}>t\right\\}}=0\right\\}.$
Then we impose the following assumption on the distributions of $\xi_{ji}$ and
$X_{ji}$.

###### Assumption 4

The designs and noises at the same time point $\\{X_{ji},\xi_{ji}\\},\
i=1,2,\dots,n_{j}$ are mutually independent for each fixed $j$ and there
exists constants $K_{1},K_{2}>0,\ \alpha\geq 1,\ \beta\geq 2$ that for any
$j=1,2\dots,T,\ i=1,2,\dots,n_{j}$

  * 1.

$\xi_{ji}$ is sub-exponential distribution with
$\|\xi_{ji}\|_{\psi(\alpha)}\leq K_{1}$;

  * 2.

$X_{ji}$ is sub-gaussian distribution with $\|X_{ji}\|_{\psi(\beta)}\leq
K_{2}$,

  * 3.

There exists $\gamma\geq 1$ such that $1/\alpha+1/\beta=1/\gamma$.

The tail distributions of $\xi_{ji}X_{ji}$ and $\langle M,X_{ji}\rangle
X_{ji}$ are both sub-exponential when $X_{ji},\xi_{ji}$ satisfy Assumption 4.
Denote $\mathcal{D}_{\alpha,\beta,K_{1},K_{2}}$ as the distribution set of
$\\{X_{ji},\xi_{ji}\\}_{ji}$ satisfying Assumption 4. We remark that
$\mathcal{D}_{\alpha,\beta,K_{1},K_{2}}$ is indeed general enough for
applications. In the matrix completion problem of Example 1, $X_{ji}$ are
bounded with $\|X_{ji}\|_{\psi(\infty)}=1$. Then the distributions of
$X_{ji},\xi_{ji}$ lie in $\mathcal{D}_{1,\infty,K_{1},1}$ as long as the
noises $\xi_{ji}$ follow sub-exponential distributions with $\alpha\geq 1$. In
the compressed sensing problem of Example 2, from Theorem 4.4.5 in Vershynin
(2018), we can set $X_{ji}$ to be sub-guassian distributed with $\beta=2$ and
$K_{2}\asymp(m_{1}\vee m_{2})^{1/2}\sigma_{X}$. Then the distributions of
$X_{ji},\xi_{ji}$ lie in $\mathcal{D}_{2,2,K_{1},\sigma_{X}}$ if noises
$\xi_{ji}$ follow a sub-guassian distributions with $\alpha=2$.

Next we consider the temporal dependence between
$\\{X_{ji},\xi_{ji}\\},i=1,2,\dots,n_{j}$ and
$\\{X_{si},\xi_{si}\\},i=1,2,\dots,n_{s}$ that decays when the distance
between time points $j$ and $s$ increases. To characterize the structure of
such dependence, we impose the strongly dependent assumption that $X_{ji}$ and
$\xi_{ji}$ are $\phi$-mixing processes across $j$. Our definition of
$\phi$-mixing is the same as in Doukhan (2012), and the detail definition is
given in Supplementary Material S.4. In the following subsections, we first
derive the estimation error bound in the ideal case that $\xi_{ji}$ and
$X_{ji}$ are independent across time $j$, and then extend to the dependent
case. These two settings are described in the following assumptions,
respectively.

###### Assumption 5

(Independent case) The designs and noises $\\{X_{ji},\xi_{ji}\\}$ are mutually
independent across $j$ for each $i$, $j=1,2,\dots,T,\ i=1,2,\dots,N_{j}$.

###### Assumption 6

($\phi$-mixing processes) There exist two independent sequences of
$\sigma$-fields $\mathcal{X}_{j}$ and $\mathcal{Y}_{j},j\in\mathbb{N}^{*}$
such that

$\sigma((X_{j1},X_{j2},\dots,X_{jn_{j}})^{\top})\subseteq\mathcal{X}_{j},\quad\sigma((\xi_{j1},\xi_{j2},\dots,\xi_{jn_{j}})^{\top})\subseteq\mathcal{Y}_{j}.$

The $\sigma$-fields $\mathcal{X}_{j}$ and $\mathcal{Y}_{j}$,
$j\in\mathbb{N}^{*}$ are both $\phi$-mixing with coefficients
$\phi_{\mathcal{X}}(k)$ and $\phi_{\mathcal{Y}}(k)$, $k\in\mathbb{N}$
satisfying that

$\Phi_{\mathcal{X}}\overset{\triangle}{=}\sum_{k=0}^{\infty}\sqrt{\phi_{\mathcal{X}}(k)}<\infty,\quad\Phi_{\mathcal{Y}}\overset{\triangle}{=}\sum_{k=0}^{\infty}\sqrt{\phi_{\mathcal{Y}}(k)}<\infty.$

###  3.2 Independent case

The theorem below states the upper bound of
$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}$ under the independent case
when the distributions of $X_{ji},\xi_{ji}$ belong to
$\mathcal{D}_{\alpha,\beta,K_{1},K_{2}}$. For conciseness, we assume that the
numbers of observations $n_{j},j=1,2,\cdots,T$ have the same order denoted by
$n$, i.e., $n_{j}\asymp n$ for $j=1,2,\cdots,T$. The conclusions for $n_{j}$
of different orders can be derived similarly. In the following, the notations
$a_{n}\ll b_{n}$, $a_{n}\gg b_{n}$ and $a_{n}\asymp b_{n}$ mean that the ratio
$a_{n}/b_{n}$ approaches to 0, infinity and a constant, respectively, and
“$\lceil\cdot\rceil$” denotes the ceiling function.

###### Theorem 2

Under Assumption 1–5, let $n\lceil
Th\rceil\gg(K_{*}/\sigma_{*})^{2}\log(m_{1}+m_{2})$, $h\rightarrow 0$ and
$(m_{1}m_{2})^{1/2}\lceil T^{2}h^{2}\rceil T^{-1}\rightarrow\infty$ as
$n,m_{1},m_{2},T\rightarrow\infty$, when

$\lambda=2C_{1}\sigma_{*}\sqrt{\frac{\log(m_{1}+m_{2})}{n\lceil Th\rceil}},$

where $C_{1}$ and $\sigma_{*}$ are defined in Lemma S.4 of Supplementary
Material, then with probability at least $1-3/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+\left(1+\sqrt{2}\right)C_{1}\sigma_{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}n\lceil
Th\rceil}\right)^{1/2}+o(h^{2}),$ (11)

where $\alpha(K)$ is defined in (10) and $D_{2}$ in Assumption 2.

We remark that $n\lceil Th\rceil$ is the number of samples used in a local
window and the condition $n\lceil
Th\rceil\gg(K_{*}/\sigma_{*})^{2}\log(m_{1}+m_{2})$ is the requirement for the
effective sample size to guarantee the upper bound (11) tending to zero. The
error bound (11) is a trade-off between the bias and variance which is
controlled by the bandwidth $h$. When
$n\mu^{2}m_{1}m_{2}/(T^{4}\sigma_{*}^{2}r_{t}\log(m_{1}+m_{2}))=o(1)$, we can
choose the optimal bandwidth

$h=C_{h}\left(\frac{\sigma_{*}^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nT}\right)^{1/5},$
(12)

where $C_{h}=\left[(2+2\sqrt{2})C_{1}/(\alpha(K)D_{2})\right]^{2/5}$ is a
constant. When
$n\mu^{2}m_{1}m_{2}/(T^{4}\sigma_{*}^{2}r_{t}\log(m_{1}+m_{2}))=\mathcal{O}(1)$,
we choose $h=0$ which is degenerate to the classical situation with the
convergence rate $n^{-1/2}$.

We focus on a large enough value for $T$ in the subsequent discussion, which
allows us to omit the ceiling operator $\lceil\cdot\rceil$. For situations
where $T$ is small, we can simply choose $h=0$, which yields the same results
as static methods. With an appropriately selected $h$, we have the following
corollary.

###### Corollary 1

Under assumptions of Theorem 2, when

$nT\gg\max\\{K_{*}^{5/2}\sigma_{*}^{-3}(\mu^{2}m_{1}m_{2})^{1/4}\log(m_{1}+m_{2}),\
\sigma_{*}^{2}(\mu^{2}m_{1}m_{2})^{-1}\log(m_{1}+m_{2})\\},$

with probability at least $1-3/(m_{1}+m_{2})$,

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{2}\left(\frac{\sigma_{*}^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nT}\right)^{2/5},$
(13)

where
$C_{2}=1/2\left[(2+2\sqrt{2})C_{1}\right]^{4/5}\left[\alpha(K)D_{2}\right]^{1/5}$.

The bound requirement for $nT$ is composed of two terms: the first is to make
$nTh\gg(K_{*}/\sigma_{*})^{2}$ $\log(m_{1}+m_{2})$ hold when $h$ takes the
form (12), and the second is to guarantee the optimal bandwidth (12) tends to
0 as $n,T\rightarrow\infty$. In (13), the upper bound is proportional to
$(nT)^{-2/5}$ where $nT$ is the total sample size. This indicates that our
dynamic method makes efficient use of information from adjacent time points to
improve the estimation quality. The number of observations in a single time
point $n$ and the number of time points $T$ complement each other to reduce
the upper bound (13). When the samples are extremely sparse at each time
point, it is infleasible for the classical static method to reconstruct the
eigenvalues and eigenvectors of $M_{t}^{0}$. However, the proposed method can
still control the estimation error
$\|\widetilde{M}^{\lambda}_{t}-M^{0}_{t}\|_{2}$ to a desirable level, as long
as there are sufficient time points. We also mention that the two-step
smoothing method performs poorly, because the first-step estimate
$\widehat{M}_{t}^{\lambda}$ and the reconstruction of eigen-subspace at each
time point are far from satisfaction, which cannot be fixed by the second
smoothing step regardless how dense the time points are. We emphasize that
Theorem 2 and Corollary 1 can be applied directly to any dynamic trace
regression problems as long as $X_{ji}$ and $\xi_{ji}$ lie in
$\mathcal{D}_{\alpha,\beta,K_{1},K_{2}}$. Taking dynamic matrix completion for
example, we have the following corollary. The corresponding result for dynamic
compressed sensing can be derived similarly, which is presented in Corollary
S.3 of the Supplement.

###### Corollary 2

Under Assumption 1–5, when $X_{ji}$ are i.i.d. uniformly distributed on
$\mathcal{E}$, $\xi_{ji}$ are independently follow sub-exponential mean-zero
distributions, $nTh\gg(m_{1}\wedge m_{2})\log^{1+2/\alpha}(m_{1}+m_{2})$, and
$(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as $h\rightarrow
0,n,m_{1},m_{2},T\rightarrow\infty$, then with probability at least
$1-3/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}$
$\displaystyle\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+(1+\sqrt{2})C_{1}(C_{M}\vee\sigma_{\xi})\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nTh}\right)^{1/2}+o(h^{2}).$

When $nT\gg(m_{1}\vee m_{2})\log^{1+5/(2\alpha)}(m_{1}+m_{2})$, let

$h=C_{h}\left(\frac{(C_{M}\vee\sigma_{\xi})^{2}r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nT}\right)^{1/5},$ (14)

then

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{2}(C_{M}\vee\sigma_{\xi})^{4/5}\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nT}\right)^{2/5},$ (15)

where $C_{1},C_{h},C_{2}$ are the same constants as above.

Recall that the classical static result for the error bound of matrix
completion in Koltchinskii et al. (2011a) is

$(m_{1}m_{2})^{-1/2}\left\|\widehat{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq(1+\sqrt{2})C_{1}(C_{M}\vee\sigma_{\xi})\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{n}\right)^{1/2}.$

The dynamic method provides a sharper upper bound than the classic one. In
other words, when the sample size in each single time point is sparse, our
method can borrow information from adjacent time points to improve the
estimation, as long as the number of time points $T$ is sufficiently large.

###  3.3 Dependent case

In parallel to the independent case, we derive the bounds of
$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}$ when $\xi_{ji},X_{ji}$ are
$\phi$-mixing processes, as stated in Theorem 3 below.

###### Theorem 3

With Assumption 1–4 and 6, when $Th\gg\log(m_{1}+m_{2})$,
$nTh\gg\left(K_{*}/\sigma_{*}\right)^{2}\log^{3}(m_{1}+m_{2})$, $h\rightarrow
0$ and $(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as
$n,m_{1},m_{2},T\rightarrow\infty$, with probability at least
$1-3/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+(1+\sqrt{2})\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\sigma_{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nTh}\right)^{1/2}+o(h^{2})$
(16)

where $\mathcal{C}_{1}>0$ is a constant defined in Lemma S.5 of Supplementary
Material.

We use $\Phi_{\mathcal{X}}$ and $\Phi_{\mathcal{Y}}$ to measure the temporal
dependence of design matrices $X_{ji}$ and noises $\xi_{ji}$, respectively,
where larger $\Phi_{\mathcal{X}}$ and $\Phi_{\mathcal{Y}}$ mean stronger
dependence. When $\Phi_{\mathcal{X}}=\Phi_{\mathcal{Y}}=1$, it indicates that
$X_{ji},\xi_{ji}$ are independent across $j$ and the bound (16) is the same as
(11) in Theorem 2 in the independent case. When
$nT\gg(\mu^{2}m_{1}m_{2})^{-1}\sigma_{*}^{2}\log(m_{1}+m_{2})$, if we choose

$h=\mathcal{C}_{h}\left(\frac{\sigma_{*}^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nT}\right)^{1/5},$
(17)

where
$\mathcal{C}_{h}=\left[(2+\sqrt{2})\mathcal{C}_{1}/(\alpha(K)D_{2})\right]^{2/5}$,
then we have the following corollary.

###### Corollary 3

Under assumptions of Theorem 3, when

$n^{-1}T^{4}\gg\sigma_{*}^{-2}\mu^{2}m_{1}m_{2}\log^{4}(m_{1}+m_{2})$

and

$nT\gg\max\\{K_{*}^{5/2}\sigma_{*}^{-3}(\mu^{2}m_{1}m_{2})^{1/4}\log(m_{1}+m_{2}),\
\sigma_{*}^{2}(\mu^{2}m_{1}m_{2})^{-1}\log(m_{1}+m_{2})\\},$

with probability at least $1-3/(m_{1}+m_{2})$,

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\mathcal{C}_{2}\left(\frac{\sigma_{*}^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nT}\right)^{2/5},$
(18)

where
$\mathcal{C}_{2}=1/2\left[(2+\sqrt{2})\mathcal{C}_{1}\right]^{4/5}\left[\alpha(K)D_{2}\right]^{1/5}$.

When $X_{ji}$ and/or $\xi_{ji}$ have strong long-term dependence across $j$,
$\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}}$ becomes large and inflates the
bound (18), which means that the convergence rate of
$\widetilde{M}_{t}^{\lambda}$ is usually slower than that in independent case.
Note that the bound for the single-stage estimates at each $t$ using static
method can also be derived from Theorem 2 by using a degenerated kernel
$\delta_{0}(x)$, given by

$(m_{1}m_{2})^{-1/2}\left\|\widehat{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq(1+\sqrt{2})C_{1}\sigma_{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}n}\right)^{1/2}.$

Then when

$\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}}\ll
T^{1/2}\left(\frac{\sigma^{*2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}n}\right)^{1/8},$
(19)

our dynamic method yields a sharper bound than the static one though strong
dependence across time exists. With the expression of $h$ in (17), the
condition (19) is equal to
$Th\gg(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}$. Note that $Th$
corresponds to the effective time points used in estimation at each $t$. This
implies that as the number of effective time points increases, the negative
effect caused by temporal correlation would be offset to some extent. For
space economy, the applications of Theorem 3 and Corollary 3 to dynamic matrix
completion and dynamic compressed sensing are presented as Corollary S.2 and
Corollary S.4 in Supplementary Material.

###  3.4 Empirical error and complexity analysis

Recall that the sequence generated by the iteration of Algorithm 2 at time $t$
are denoted by $\\{M^{(k)}_{t}\\}_{k},k\in\mathbb{N}^{*}$, where $M^{(k)}_{t}$
is the output matrix after $k$-th iteration and $M_{t}^{(0)}$ is the initial
input matrix. In practice, we use the output $M^{(k)}_{t}$ to approximate
$M_{t}^{0}$ for a certain $k$. Thus, we derive the upper bound of the
empirical error $\big{\|}M_{t}^{(k)}-M_{t}^{0}\big{\|}_{2}$ in the following
theorem, which takes both the algorithmic error between $M^{(k)}_{t}$ and
$\widetilde{M}_{t}^{\lambda}$ and the estimation error between
$\widetilde{M}_{t}^{\lambda}$ and $M_{t}^{0}$ into account.

###### Theorem 4

Under assumptions of Theorem 1, for any $k>1$ and $\lambda\geq
2\|\mathcal{W}_{h}\Delta_{t}\|_{\infty}$, we have

$\displaystyle\left\|M_{t}^{(k)}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\min\left\\{\left((\delta_{h}M_{t})^{2}+\frac{2\lambda}{\mu}\|M_{t}^{0}\|_{1}+\frac{2L_{f}\gamma_{t}^{2}}{\mu(k+1)^{2}}\right)^{1/2},\right.$
$\displaystyle\left.\frac{1+\sqrt{2}}{2}\frac{\lambda}{\mu}\sqrt{r_{t}}+\left(\frac{2L_{f}\gamma_{t}^{2}}{\mu(k+1)^{2}}\right)^{1/2}\right\\}$

where $\gamma_{t}=\big{\|}M_{t}^{(0)}-\widetilde{M}^{\lambda}_{t}\big{\|}_{2}$
and
$L_{f}=2\big{\|}\sum_{j=1}^{T}\omega_{h}(j-t)\cdot\frac{1}{n_{j}}\sum_{i=1}^{n_{t}}\left(\|X_{ji}\|_{2}X_{ji}\right)\big{\|}_{2}$
is the Lipschitz constant of gradient function $\nabla f_{t}(M)$ that is
finite under Assumption 4.

The empirical error bound in Theorem 4 has an additional term
$2L_{f}\gamma_{t}^{2}/\mu(k+1)^{2}$ compared to the estimation error bound (7)
in Theorem 1, which is the so-called algorithmic error. The optimal choice is
to stop iterating when the algorithmic error attains the same order as the
estimation error. Let $k_{t}$ be the minimal number of iteration steps needed
to attain the desired order. From the expression of the algorithmic error,
$k_{t}$ depends on the distance between the initial value $M_{t}^{(0)}$ and
the target $\widetilde{M}^{\lambda}_{t}$, denoted by $\gamma_{t}$. Since less
iterations are required for a smaller initial distance, the naive strategy to
initialize $M_{t}^{(0)}$ randomly at each time $t$ would lead a large
$\gamma_{t}$. Noting that $\widetilde{M}_{t}^{\lambda}$ is close to
$\widetilde{M}_{t-1}^{\lambda}$, we thus make use of the neighboring
information in the initialization to reduce computation cost. Hence we propose
to use the output matrix at time $(t-1)$ as the initial value for estimating
$\widetilde{M}_{t}^{\lambda}$ at time $t$, i.e.
$M_{t}^{(0)}=M_{t-1}^{(k_{t-1})}$. Take the dependent case in Section 3.3 for
instance. Recall that the optimal estimation error bound is derived in
Corollary 3 with appropriate choice of the bandwidth $h$. Then desired
empirical error bound for Theorem 4 can be derived by

$(m_{1}m_{2})^{-1/2}\left\|M_{t}^{(k_{t})}-M_{t}^{0}\right\|_{2}\leq
2\mathcal{C}_{2}\left(\frac{\sigma_{*}^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nT}\right)^{2/5}.$
(20)

To quantify the improvement of our initialization strategy, in Corollary 4, we
compute the total number of iteration steps to attain the bound (20) at each
$t$, and compare with that of the random initialization. The result for
independent case can be derived similarly, which is presented in Corollary S.5
of Supplementary Material for space economy.

###### Corollary 4

Under the assumptions of Theorem 3, when the total number of iterations of
random initialization $K_{0}$ satisfies

$K_{0}\geq\mathcal{C}_{2}^{-1}\left(\frac{\sigma_{*}^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}}\right)^{-2/5}\left(\frac{2L_{f}}{\mu
m_{1}m_{2}}\right)^{1/2}\gamma_{1}n^{2/5}T^{7/5},$

or total number of iterations of the proposed initial strategy $K_{1}$
satisfies

$\displaystyle K_{1}\geq$
$\displaystyle\mathcal{C}_{2}^{-1}\left(\frac{\sigma_{*}^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}}\right)^{-2/5}\left(\frac{2L_{f}}{\mu
m_{1}m_{2}}\right)^{1/2}(\gamma_{1}+D_{1})(nT)^{2/5}$
$\displaystyle+3\sqrt{2}\left(\frac{L_{f}}{\mu}\right)^{1/2}T.$

the empirical error is lower than the optimal upper bound (20) for each
$t=1,\ldots,T$.

From Corollary 4, once the total number of iterations of random initialization
attains the level $\mathcal{O}(n^{2/5}T^{7/5})$ and that of Algorithm 2
attains $\mathcal{O}(n^{2/5}T^{2/5}+T)$, the error is bounded by (20). Then
the ratio of computational cost of our proposed algorithm compared to the
random initialization method is $\mathcal{O}(T^{-1}\vee n^{-2/5}T^{-2/5})$. We
finally mention that the classical static method need $\mathcal{O}(n^{1/2})$
iterations to attain the optimal rate $\mathcal{O}(n^{-1/2})$ for each $t$,
and hence the single-stage method (at all $T$ times) need
$\mathcal{O}(n^{1/2}T)$ iterations in total. Therefore the ratio of
computational cost of Algorithm 2 compared to the single-stage method is
$\mathcal{O}(T^{-3/5}\vee n^{-1/2})$, which leads to the conclusion that our
algorithm is computationally more efficient than the random initial method and
the classical static methods.

##  4 Simulation Studies

We take the dynamic matrix completion as the instance to verify the usefulness
of our method. To demonstrate the advantages, we compare the proposed
estimates with three benchmark methods: (i) the classical static estimator
$\widehat{M}_{t}^{\lambda}$ defined as in (4), abbreviated as Static; (ii) the
two-step smoothing estimator $\widehat{M}_{t}^{\text{t-s}}$ given in (5),
abbreviated as TwoStep; (iii) the low-rank tensor completion estimator,
abbreviated as Tensor,

$\widehat{N}^{\lambda}=\mathop{\rm argmin}_{N\in\mathbb{R}^{T\times
m_{1}\times
m_{2}}}\|\mathcal{P}_{\Omega}(N)-Y\|_{2}^{2}+\lambda\sum_{i=1}^{3}\|N_{(i)}\|_{1},$

where $\mathcal{P}_{\Omega}$ is the projection to the observation subspace
$\Omega$, Y is the observed data tensor with element

$\displaystyle[Y]_{j,s,l}=\left\\{\begin{aligned}
&\frac{1}{|\mathcal{A}_{j,s,l}|}\sum_{i\in\mathcal{A}_{j,s,l}}Y_{ji},\ \text{
if }\mathcal{A}_{j,s,l}\neq\emptyset\\\ &0,\ \text{
otherwise}\end{aligned}\right.$

for $\mathcal{A}_{j,s,l}=\left\\{i\mid
X_{ji}=e_{s}(m_{1})e^{\top}_{l}(m_{2}),1\leq i\leq n_{j}\right\\}$, and
$N_{(i)}$ is the mode-$i$ unfolding of tensor $N$ which is a matrix that
arranges the mode-$i$ fibers to be columns. We use matrix form FISTA algorithm
in Toh & Yun (2010) to approximate $\widehat{M}_{t}^{\lambda}$, and
$\widehat{M}^{\text{t-s}}_{t}$ has a closed form solution. In Tensor, we use
AMD-TR(E) algorithm in Gandy et al. (2011) to numerically approximate
$\widehat{N}^{\lambda}$.

We set $m_{1}=500,m_{2}=300,r_{t}=10$ and the matrix function $M(t)$ in
Assumption 2 is constructed by $M(t)=U(t)D(t)V^{T}(t)$ with

$\displaystyle U(t)$
$\displaystyle=\cos\left(\frac{t\pi}{2}\right)U_{0}+\sin\left(\frac{t\pi}{2}\right)U_{1},$
$\displaystyle D(t)$
$\displaystyle=10\left(\operatorname{diag}\left\\{10^{2},9^{2},\dots,1\right\\}+t\operatorname{diag}\left\\{10,9,\dots,1\right\\}\right),$
$\displaystyle V(t)$
$\displaystyle=\cos\left(\frac{t\pi}{2}\right)V_{0}+\sin\left(\frac{t\pi}{2}\right)V_{1},$

where $U_{0},U_{1}\in\mathbb{R}^{m_{1}\times r_{t}}$,
$V_{0},V_{1}\in\mathbb{R}^{m_{2}\times r_{t}}$ are composed of standard
orthogonal basis. Here we consider the situation that each $n_{t}$ is the
same, i.e., $n_{t}=n,t=1,2,\cdots,T$. Denote the rate of observation samples
as $\rho=n/(m_{1}m_{2})$ and the density of time points as $\tau=1/T$. We use
a 5-fold cross validation to select $\lambda$ in the numerical experiments.
Please refer to Section S.9.4 of the Supplement for more details.

For the independent case, we sample $X_{ji},\ i=1,2,\dots,n$ from
$\mathcal{E}$ defined in (3) independently with uniform probabilities, and let
$\xi_{ji},j=1,2,\dots,T,\ i=1,2,\dots,n$ i.i.d follow the normal distribution
$N(0,\sigma_{\xi}^{2})$, where $T=100$ and $\sigma_{\xi}=1$. We use the
Epanechnikov kernel and set the bandwidth $h$ by the optimal selection form
(14) in which we approximate $C_{M}\vee\sigma_{\xi}$ by mean of the top 10%
largest $Y_{ji}$ and $D_{2}$ by
$\sum_{j=2}^{T}|1/n_{j}\sum_{i=1}^{n_{j}}Y_{ji}-1/n_{j-1}\sum_{i=1}^{n_{j-1}}Y_{j-1i}|$.
We set $\rho=0.2$ for our proposed dynamic low-rank (DLR) and Tensor methods.
In Static and TwoStep, $\widehat{M}_{t}^{\lambda}$ and
$\widehat{M}_{t}^{\text{t-s}}$ fail to recover the matrix structure under the
same setting $\rho=0.2$, so we set $\rho=0.8$ for feasibility. Figure 1
illustrates the mean square errors

$\mbox{MSE}_{t}=(m_{1}m_{2})^{-1}\left\|\widehat{M}_{t}-M^{0}_{t}\right\|_{2}^{2},$
(21)

for $t=1,2,\cdots,T$. For the proposed DLR and TwoStep methods, the MSEt
floats up at the boundaries due to slight edge effect. The average MSE, i.e.,
$T^{-1}\sum_{t=1}^{T}\mbox{MSE}_{t}$ for the DLR estimate is $0.085$ and for
Tensor is $0.44$ when $\rho=0.2$. This implies that the tensor recovery is not
effective at least when the temporal structure is smooth. For Static and
TwoStep, when $\rho=0.8$, the average MSE is $0.16$ and $0.13$ respectively.
That is, these two methods with four times data still perform worse than the
proposed DLR method. This verifies the theoretical finding that when the
sample size $n$ at a single time point is small, the classical trace
regression fails to estimate $M^{0}_{t}$, while our dynamic method can borrow
adjacent information to improve the estimation quality.

Figure 1: The comparisons of MSEt for different methods. In DLR and Tensor,
$\rho=0.2$, and in Static and TwoStep, $\rho=0.8$.

Now we study the empirical influence of $\rho$ and $\tau$ on the estimation
results. One can see from (15) that given the underlying matrix, the bound of
MSE is in proportion to $(\rho/\tau)^{-4/5}$. We set eight different values
for $(\rho,\tau)$ which can be divided into three groups, i.e.
$\rho/\tau=5,10,20$. The $MSE_{t}$ across $t$ under these settings are shown
in the left panel of Figure 2, which indicates that given the matrix
dimension, the empirical MSE depends on $\rho/\tau$. We further plot the
logarithm of average MSE across time versus the logarithm of $\rho/\tau$ in
the right panel of Figure 2, which reveals the linear relationship of the
slope -4/5 between them. This further validates the conclusion (15).

Figure 2: The left panel corresponds to MSEt of DLR estimates under different
settings of $\rho$ and $\tau$. The curves are clustered into three classes
corresponding to $\rho/\tau=5,10,20$. The right panel shows the logarithm of
average MSE versus the logarithm of $\rho/\tau$ showing a linear trend with
the slope -4/5.

We also investigate the influence of the dependent noise on the proposed
estimator. The noises $\xi_{ji}$ are generated by

$\displaystyle E_{0}$ $\displaystyle=U_{0},\quad
U_{0}\sim\operatorname{WN}(0,\sigma_{\xi}^{2}\mathbf{1}_{m_{1}}\mathbf{1}^{T}_{m_{2}}),$
$\displaystyle E_{t}$ $\displaystyle=\beta
E_{t-1}+\sqrt{1-\beta^{2}}U_{t},\quad
U_{t}\sim\operatorname{WN}(0,\sigma_{\xi}^{2}\mathbf{1}_{m_{1}}\mathbf{1}^{T}_{m_{2}}),$
(22) $\displaystyle\xi_{ji}$ $\displaystyle=\langle E_{j},X_{ji}\rangle,\quad
j=1,2,\dots,T,\ i=1,2,\dots,n_{t},$

where
$\operatorname{WN}(0,\sigma_{\xi}^{2}\mathbf{1}_{m_{1}}\mathbf{1}^{T}_{m_{2}})$
generates the $m_{1}\times m_{2}$ random matrix whose elements are i.i.d.
normal with mean zero and variance $\sigma_{\xi}^{2}$. The larger $\beta$ is,
the stronger dependence exists among $\\{\xi_{ji}\\}_{j}$. The other
parameters are set are the same as those in the independent case. Then
$\Phi_{\mathcal{X}}=1$ and $\Phi_{\mathcal{Y}}$ can be calculated from (4)
which is larger than 1 when $\beta>0$. We apply the proposed method to recover
$M_{t}^{0}$ for $1\leq t\leq T$ under the combinations of $\sigma_{\xi}=1,2$
and $\beta=0,0.3,0.6,0.9$. As illustrated in the left panel of Figure 3, the
MSE increases as $\sigma_{\xi}^{2}$ and/or $\beta$ enlarge(s). Note that
according to (S.26) in Supplement, the MSE is proportional to
$\Phi_{\mathcal{Y}}^{8/5}$ given $M_{t},n,T$. We plot the average MSE versus
$\Phi_{\mathcal{Y}}^{8/5}$ under different $\sigma_{\xi}^{2}$ in the right
panel of Figure 3 which verifies the conclusion.

Figure 3: The left panel illustrates the MSEt across time with different
$\sigma_{\xi}$ and $\beta$ and the right panel shows that the average MSE
versus $\Phi_{\mathcal{Y}}^{8/5}$ under different settings when
$\Phi_{\mathcal{X}}=1$.

Figure 4: The left panel illustrates the MSEt across time with different
$\sigma_{\xi}$ and $\alpha$ and the right panel shows that the average MSE
versus $\Phi_{\mathcal{X}}^{8/5}$ under different settings when
$\Phi_{\mathcal{Y}}=1$.

Finally we study the influence of the dependence of design matrices on the
estimation. We sample $X_{1i},i=1,\dots,n$ from $\mathcal{E}$ with uniform
probabilities independently and generate $X_{ji},j\geq 2,i=1,2,\dots,n$ by
preserving $\alpha n$ elements randomly in $X_{j-1,i},i=1,2,\dots,n$ and
choosing $(1-\alpha)n$ new elements i.i.d. from $\mathcal{E}$. Then the
dependence of design matrices across time becomes stronger as $\alpha$
increases. We set the remaining parameters the same as those in the
independent case and have $\Phi_{\mathcal{Y}}=1$, $\Phi_{\mathcal{X}}\geq 1$,
where the equality holds if and only if $\alpha=0$. We plot the
$\operatorname{MSE}_{t}$ across time under different $\sigma_{\xi}=1,2$ and
$\alpha=0,0.3,0.6,0.9$ and verify the linear relationship between the average
MSE and $\Phi_{\mathcal{X}}^{8/5}$ in Figure 4. We remark that the slopes of
linear relationships in the right panels of Figure 3 and 4 are different. This
is because the bound (S.26) in Supplement is obtained by the technique that
$C_{M}\Phi_{\mathcal{X}}+\sigma_{\xi}\Phi_{\mathcal{Y}}\leq(C_{M}\vee\sigma_{\xi})(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})$.
Hence the curves of the average MSE versus $\Phi_{\mathcal{Y}}^{8/5}$ become
sharper as $\sigma_{\xi}$ increases, as shown in the right panel of Figure 3,
while the slopes of the average MSE versus $\Phi_{\mathcal{X}}^{8/5}$ under
different $\sigma_{\xi}$ in the right panel of Figure 4 are similar.

##  5 Real Data Examples

In this section, we apply the proposed method to two real data examples for
the dynamic matrix completion and compressed sensing problems, respectively.

The first example is the Netflix Prize Dataset (Netflix, 2006), in which users
ratings for movies received by Netflix were collected from October 1998 to
December 2005 and can be downloaded at
https://www.kaggle.com/datasets/netflix-inc/netflix-prize-data. The ratings
are integer-valued from 1 to 5. We preprocess the dataset to remove those
movies which are watched less than 25000 times.

For user selection, we consider two different filters. Filter 1 is to choose
users whose rating times are more than 700 because users who watch a larger
number of movies tend to offer more reliable ratings and make better
comparisons across different movies due to their broader exposure. There
remain $3036$ users and $1034$ movies with total 2337997 ratings. We split the
observations into $T$ time intervals in chronological order, then randomly
choose $4/5$ of them in each interval as training data and set the rest as
test data. The observations within each time interval is collapsed as observed
on a time point. At each time point $t$, there exists an underlying matrix
$M_{t}\in\mathbb{R}^{m_{1}\times m_{2}}$ to be recovered, where the rows
represent users and the columns represent movies. The observed sparse matrix
$Y_{t}\in\mathbb{R}^{m_{1}\times m_{2}}$ has elements at position $(j,k)$,
representing the rating given by user $j$ to movie $k$ at time point $t$, or
zero if no rating is available. Thus we can convert it to the trace regression
model. The average sample size at each time point is about $17803$ and thus
the rate of observation samples in each time point is $\rho\approx 0.0057$. We
use again the Epanechnikov kernel $K(x)$ and select the bandwidth $h$ in the
same way as in the simulations. Filter 2 is to select users whose rating times
are more than 30 because those users can represent the target groups of users
who often using Netflix to watch movies. Then we conduct a random selection of
3000 individuals from the pool of eligible users for ease of computation.

We use 5-fold cross validation to choose tuning parameter $\lambda$ and
evaluate the performance of $\widetilde{M}^{\lambda}_{t}$ by calculating

$\mbox{MSE}^{*}_{t}=\frac{1}{n_{t}^{\text{test}}}\sum_{i=1}^{n_{t}^{\text{test}}}\left(\langle
X_{t,i}^{\text{test}},\widetilde{M}^{\lambda}_{t}\rangle-
Y_{t,i}^{\text{test}}\right)^{2},$ (23)

using the test data. To investigate the influence of different values of $T$,
we conduct an experiment in which $T$ is varied from 1 to 1000. Detailed
information about this experiment is presented in Supplement S.9.1. The
results indicate that for sufficiently large $T$, the variation in $T$ has a
negligible impact while when $T$ is too small, increasing its value leads to a
reduction in estimation error. This observation also corresponds to the
dynamic nature inherent in the underlying low-rank matrix, signifying its
dynamic essence rather than a static state. Then we set $T=100$ as usual and
compare the proposed method with those three benchmarks mentioned in Section
4. Figure 5 shows that our dynamic method is more accurate than those
benchmarks across the time domain. And the average MSE∗, i.e.
$T^{-1}\sum_{t=1}^{T}\mbox{MSE}^{*}_{t}$, of our dynamic method and benchmarks
Static, TwoStep and Tensor methods are compared in Table 1

Figure 5: The MSE${}^{*}_{t}$ (23) at different time points for the proposed method and three benchmarks for the Netflix dataset. The left plot is the result for Filter 1 and the right plot is the result for Filter 2. Table 1: The average MSE∗ for Netflix dataset. Setting | DLR | Static | TwoStep | Tensor
---|---|---|---|---
Filter 1 MSE∗ | 0.781 | 1.431 | 1.397 | 1.082
Filter 2 MSE∗ | 0.796 | 1.601 | 1.565 | 1.360

Another example is the compression and recovery of videos. We use the lion
video from Davis 2017 dataset (Pont-Tuset et al., 2017) and treat each frame
as a matrix at each time point.The dataset is publicly available at
https://davischallenge.org/davis2017/code.html#unsupervised. For compression,
we first separate the matrices $M^{0}_{t}$ into the sparse part $S_{t}$ and
the low-rank part $L_{t}$ using robust principal component analysis (Candès et
al., 2011). The design matrices $X_{ti}$ are the convolution kernel matrices
with random centers $(a_{ti},b_{ti})$, i.e.,
$[X_{ti}]_{a_{ti},b_{ti}}=4,[X_{ti}]_{a_{ti}\pm
1,b_{ti}}=[X_{ti}]_{a_{ti},b_{ti}\pm 1}=2,[X_{ti}]_{a_{ti}\pm 1,b_{ti}\pm
1}=1$. It is straightforward to verify that the distributions of
$\\{X_{ti}\\}_{t,i}$ satisfy Assumption 4. Then we compress the low-rank part
to obtain corresponding output $Y_{ti}=\langle X_{ti},L_{t}\rangle$. For
recovery, we use the dynamic trace regression to reconstruct the low-rank
part, denoted as $\widetilde{L}_{t}^{\lambda}$. Then we recover each frame by
$\widetilde{M}_{t}^{\lambda}=S_{t}+\widetilde{L}^{\lambda}_{t}$, while the
sparse part $S_{t}$ is retained. We emphasize that we only need store the
centers $(a_{ti},b_{ti})$, the output $Y_{ti}$ and the nonzero part of $S_{t}$
to recover the matrices $M_{t}^{0}$, which is considerably space-saving
compared to storing the whole video. In the lion video, each frame has
$480\times 854$ pixels. We generate $Y_{ti}$ with the rate of observation
samples $\rho=0.146$ and the file volume is reduced by 70% after compression.
Here we use Static and TwoStep mentioned in Section 4 to reconstruct
$\widehat{L}_{t}^{\lambda},\widehat{L}_{t}^{\text{t-s}}$ and recover
$\widehat{M}^{\lambda}_{t}=S_{t}+\widehat{L}_{t}^{\lambda}$,
$\widehat{M}^{\text{t-s}}=S_{t}+\widehat{L}_{t}^{\text{t-s}}$ as benchmarks.
Figure 6 shows the original frames $M_{t}^{0}$ and the estimates
$\widetilde{M}^{\lambda}_{t},\widehat{M}_{t}^{\lambda},\widehat{M}^{\text{t-s}}$
at different time points. The background and the smooth motions of all three
lions are recovered well by the proposed DLR method, which can be seen
visually better than two benchmark methods. Recall the definition (23), the
average MSE${}^{*}=T^{-1}\sum_{t=1}^{T}\mbox{MSE}^{*}_{t}$ is $2.9\times
10^{-3}$ for our dynamic method while is $7.6\times 10^{-2}$ and $4.7\times
10^{-2}$ for the two benchmarks Static and TwoStep repsectively.

Figure 6: The top five pictures are the $5$th, $25$th, $45$th, $65$th and
$85$th original frames in the lions video, the second rows are the
corresponding estimates of the proposed DLR method and the third and fourth
rows are the benchmarks Static and TwoStep, respectively.

## References

  * (1)
  * Agarwal et al. (2012) Agarwal, A., Negahban, S. & Wainwright, M. J. (2012), ‘Noisy matrix decomposition via convex relaxation: Optimal rates in high dimensions’, The Annals of Statistics 40(2), 1171–1197.
  * Beck & Teboulle (2009) Beck, A. & Teboulle, M. (2009), ‘A fast iterative shrinkage-thresholding algorithm for linear inverse problems’, SIAM journal on imaging sciences 2(1), 183–202.
  * Bi et al. (2021) Bi, X., Tang, X., Yuan, Y., Zhang, Y. & Qu, A. (2021), ‘Tensors in statistics’, Annual review of statistics and its application 8, 345–368.
  * Candès et al. (2011) Candès, E. J., Li, X., Ma, Y. & Wright, J. (2011), ‘Robust principal component analysis?’, Journal of the ACM (JACM) 58(3), 1–37.
  * Candes & Plan (2011) Candes, E. J. & Plan, Y. (2011), ‘Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements’, IEEE Transactions on Information Theory 57(4), 2342–2359.
  * Candès & Recht (2009) Candès, E. J. & Recht, B. (2009), ‘Exact matrix completion via convex optimization’, Foundations of Computational mathematics 9(6), 717–772.
  * Csurgay et al. (2000) Csurgay, Á. I., Porod, W. & Lent, C. S. (2000), ‘Signal processing with near-neighbor-coupled time-varying quantum-dot arrays’, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 47(8), 1212–1223.
  * Doukhan (2012) Doukhan, P. (2012), Mixing: properties and examples, Vol. 85, Springer Science & Business Media.
  * Fan et al. (2021) Fan, J., Wang, W. & Zhu, Z. (2021), ‘A shrinkage principle for heavy-tailed data: High-dimensional robust low-rank matrix recovery’, Annals of statistics 49(3), 1239.
  * Gandy et al. (2011) Gandy, S., Recht, B. & Yamada, I. (2011), ‘Tensor completion and low-n-rank tensor recovery via convex optimization’, Inverse problems 27(2), 025010.
  * Gao & Wang (2018) Gao, P. & Wang, M. (2018), Dynamic matrix recovery from partially observed and erroneous measurements, in ‘2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)’, IEEE, pp. 4089–4093.
  * Gross et al. (2010) Gross, D., Liu, Y.-K., Flammia, S. T., Becker, S. & Eisert, J. (2010), ‘Quantum state tomography via compressed sensing’, Physical review letters 105(15), 150401.
  * Hang & Steinwart (2017) Hang, H. & Steinwart, I. (2017), ‘A bernstein-type inequality for some mixing processes and dynamical systems with an application to learning’, The Annals of Statistics 45(2), 708–743.
  * Hu & Yao (2022) Hu, X. & Yao, F. (2022), ‘Dynamic principal component analysis in high dimensions’, Journal of the American Statistical Association pp. 1–12.
  * Ji & Ye (2009) Ji, S. & Ye, J. (2009), An accelerated gradient method for trace norm minimization, in ‘Proceedings of the 26th annual international conference on machine learning’, pp. 457–464.
  * Keshavan et al. (2009) Keshavan, R., Montanari, A. & Oh, S. (2009), ‘Matrix completion from noisy entries’, Advances in neural information processing systems 22.
  * Klopp (2011) Klopp, O. (2011), ‘Rank penalized estimators for high-dimensional matrices’, Electronic Journal of Statistics 5, 1161–1183.
  * Koltchinskii (2011b) Koltchinskii, V. (2011b), ‘Von neumann entropy penalization and low-rank matrix estimation’, The Annals of Statistics 39(6), 2936–2973.
  * Koltchinskii et al. (2011a) Koltchinskii, V., Lounici, K. & Tsybakov, A. B. (2011a), ‘Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion’, The Annals of Statistics 39(5).
  * Koltchinskii & Xia (2015) Koltchinskii, V. & Xia, D. (2015), ‘Optimal estimation of low rank density matrices.’, J. Mach. Learn. Res. 16(53), 1757–1792.
  * Koren (2009) Koren, Y. (2009), Collaborative filtering with temporal dynamics, in ‘Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining’, pp. 447–456.
  * Koren et al. (2009) Koren, Y., Bell, R. & Volinsky, C. (2009), ‘Matrix factorization techniques for recommender systems’, Computer 42(8), 30–37.
  * Li et al. (2019) Li, X. P., Huang, L., So, H. C. & Zhao, B. (2019), ‘A survey on matrix completion: Perspective of signal processing’, arXiv preprint arXiv:1901.10885 .
  * Liu et al. (2012) Liu, J., Musialski, P., Wonka, P. & Ye, J. (2012), ‘Tensor completion for estimating missing values in visual data’, IEEE transactions on pattern analysis and machine intelligence 35(1), 208–220.
  * Lois & Vaswani (2015) Lois, B. & Vaswani, N. (2015), Online matrix completion and online robust pca, in ‘2015 IEEE International Symposium on Information Theory (ISIT)’, IEEE, pp. 1826–1830.
  * Merlevède et al. (2009) Merlevède, F., Peligrad, M. & Rio, E. (2009), Bernstein inequality and moderate deviations under strong mixing conditions, in ‘High dimensional probability V: the Luminy volume’, Institute of Mathematical Statistics, pp. 273–292.
  * Merlevède et al. (2011) Merlevède, F., Peligrad, M. & Rio, E. (2011), ‘A bernstein type inequality and moderate deviations for weakly dependent sequences’, Probability Theory and Related Fields 151(3), 435–474.
  * Negahban & Wainwright (2011) Negahban, S. & Wainwright, M. J. (2011), ‘Estimation of (near) low-rank matrices with noise and high-dimensional scaling’, The Annals of Statistics 39(2), 1069–1097.
  * Nesterov (1983) Nesterov, Y. E. (1983), A method of solving a convex programming problem with convergence rate $o(\frac{1}{k^{2}})$, in ‘Doklady Akademii Nauk’, Russian Academy of Sciences, pp. 543–547.
  * Netflix (2006) Netflix (2006), ‘Netflix prize data’, http://netflixprize.com/index.html.
  * Pont-Tuset et al. (2017) Pont-Tuset, J., Perazzi, F., Caelles, S., Arbeláez, P., Sorkine-Hornung, A. & Van Gool, L. (2017), ‘The 2017 davis challenge on video object segmentation’, arXiv:1704.00675 .
  * Qiu et al. (2017) Qiu, K., Mao, X., Shen, X., Wang, X., Li, T. & Gu, Y. (2017), ‘Time-varying graph signal reconstruction’, IEEE Journal of Selected Topics in Signal Processing 11(6), 870–883.
  * Recht et al. (2010) Recht, B., Fazel, M. & Parrilo, P. A. (2010), ‘Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization’, SIAM review 52(3), 471–501.
  * Samson (2000) Samson, P.-M. (2000), ‘Concentration of measure inequalities for markov chains and phi-mixing processes’, The Annals of Probability 28(1), 416–461.
  * Toh & Yun (2010) Toh, K.-C. & Yun, S. (2010), ‘An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems’, Pacific Journal of optimization 6(615-640), 15.
  * Vershynin (2018) Vershynin, R. (2018), High-dimensional probability: An introduction with applications in data science, Vol. 47, Cambridge university press.
  * Wang et al. (2016) Wang, X., Donaldson, R., Nell, C., Gorniak, P., Ester, M. & Bu, J. (2016), Recommending groups to users using user-group engagement and time-dependent matrix factorization, in ‘Proceedings of the AAAI Conference on Artificial Intelligence’.
  * Watson (1992) Watson, G. A. (1992), ‘Characterization of the subdifferential of some matrix norms’, Linear algebra and its applications 170(0), 33–45.
  * Xia & Yuan (2021) Xia, D. & Yuan, M. (2021), ‘Statistical inferences of linear forms for noisy matrix completion’, Journal of the Royal Statistical Society Series B 83(1), 58–77.
  * Xu & Davenport (2016) Xu, L. & Davenport, M. (2016), ‘Dynamic matrix recovery from incomplete observations under an exact low-rank constraint’, Advances in Neural Information Processing Systems 29.
  * Xu et al. (2020) Xu, X., Dong, F., Li, Y., He, S. & Li, X. (2020), Contextual-bandit based personalized recommendation with time-varying user interests, in ‘Proceedings of the AAAI Conference on Artificial Intelligence’, Vol. 34, pp. 6518–6525.
  * Zhang et al. (2020) Zhang, A. R., Luo, Y., Raskutti, G. & Yuan, M. (2020), ‘Islet: Fast and optimal low-rank tensor regression via importance sketching’, SIAM journal on mathematics of data science 2(2), 444–479.
  * Zhang et al. (2021) Zhang, Y., Bi, X., Tang, N. & Qu, A. (2021), ‘Dynamic tensor recommender systems’, The Journal of Machine Learning Research 22(1), 3032–3066.
  * Zhao et al. (2015) Zhao, T., Wang, Z. & Liu, H. (2015), ‘A nonconvex optimization framework for low rank matrix estimation’, Advances in Neural Information Processing Systems 28.
  * Zheng & Lafferty (2016) Zheng, Q. & Lafferty, J. (2016), ‘Convergence analysis for rectangular matrix completion using burer-monteiro factorization and gradient descent’, arXiv preprint arXiv:1605.07051 .
  * Zhou et al. (2013) Zhou, H., Li, L. & Zhu, H. (2013), ‘Tensor regression with applications in neuroimaging data analysis’, Journal of the American Statistical Association 108(502), 540–552.

Supplementary Material for
“Dynamic Matrix Recovery”

##  S.1 Proof of Theorem 1

We need the following lemma for proving Theorem 1.

###### Lemma S.1

Let $\Sigma\in\mathbb{S}_{+}^{m\times m}$ be a positive definite matrix with
the smallest eigenvalue $\mu>0$. Then the following inequality

$\alpha^{\top}\Sigma(\alpha+\beta)\geq\mu\left(\|\alpha\|_{2}^{2}-\|\alpha\|_{2}\|\beta\|_{2}\right)$
(S.1)

holds for any $\alpha,\beta\in\mathbb{R}^{m}$ with
$\|\alpha\|_{2}\geq\|\beta\|_{2}$.

Proof of Lemma S.1: Note that for a fixed $\alpha$,

$\alpha^{\top}\Sigma(\alpha+\beta)=\alpha^{\top}\Sigma\alpha+\langle\Sigma\alpha,\beta\rangle\geq\alpha^{\top}\Sigma\alpha-\|\beta\|_{2}\|\Sigma\alpha\|_{2}.$

Let $\lambda_{1},\ldots,\lambda_{m}$ be the eigenvalues of $\Sigma$ satisfying
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{m}=\mu>0$. Using the
representation $\alpha=\sum_{s=1}^{m}a_{s}\psi_{s}$ where $\psi_{s}$ are
orthogonal eigenvectors corresponding to $\lambda_{s}$, we have

$\alpha^{\top}\Sigma\alpha-\|\beta\|_{2}\|\Sigma\alpha\|_{2}=\sum_{s=1}^{m}a_{s}^{2}\lambda_{s}-\|\beta\|_{2}\left(\sum_{s=1}^{m}a_{s}^{2}\lambda^{2}_{s}\right)^{1/2}.$

Consider the minimization problem

$\displaystyle\min$
$\displaystyle\quad\sum_{s=1}^{m}a_{s}^{2}\lambda_{s}-\|\beta\|_{2}\left(\sum_{s=1}^{m}a_{s}^{2}\lambda^{2}_{s}\right)^{1/2}$
$\displaystyle\operatorname{s.t.}$
$\displaystyle\quad\sum_{s=1}^{m}a_{s}^{2}=\|a\|_{2}^{2}.$

This is a convex problem with

$\mathcal{L}(a_{1},\dots,a_{m},\eta)=\sum_{s=1}^{m}a_{s}^{2}\lambda_{s}-\|\beta\|_{2}\left(\sum_{s=1}^{m}a_{s}^{2}\lambda^{2}_{s}\right)^{1/2}-\eta\left(\sum_{s=1}^{m}a_{s}^{2}-\|a\|_{2}^{2}\right)$

and

$\frac{\partial\mathcal{L}}{\partial
a_{s}}=2a_{s}\left(\lambda_{s}-\frac{1}{2}\|\beta\|_{2}\lambda^{2}_{s}\left(\sum_{s=1}^{m}a_{s}^{2}\lambda^{2}_{s}\right)^{-1/2}-\eta\right).$

With Karush-Kuhn-Tucker condition, we have minimum values
$a_{s},s=1,2,\dots,m$ satisfying that only one of $a_{s}$ is equal to
$\|a\|_{2}$ and others are all zero, which means that

$\alpha^{\top}\Sigma(\alpha+\beta)\geq\min_{s=1,2,\dots,m}\left\\{\lambda_{s}\left(\|\alpha\|_{2}^{2}-\|\beta\|_{2}\|\alpha\|_{2}\right)\right\\}\geq\mu(\|\alpha\|_{2}^{2}-\|\alpha\|_{2}\|\beta\|_{2}).$

$\square$

Proof of Theorem 1: For brevity, we denote

$\omega_{j}=\omega_{h}(j-t),\quad\mathcal{C}_{j}(M,N)=\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}\langle
M,X_{ji}\rangle\langle N,X_{ji}\rangle$

for all $M,N\in\mathbb{R}_{m_{1}\times m_{2}}$. With the definition of
$\widetilde{M}_{t}^{\lambda}$ ,

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda},\widetilde{M}_{t}^{\lambda})-2\sum_{j=1}^{T}\omega_{j}\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},\widetilde{M}_{t}^{\lambda}\right\rangle+\lambda\|\widetilde{M}_{t}^{\lambda}\|_{1}$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M,M)-2\sum_{j=1}^{T}\omega_{j}\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},M\right\rangle+\lambda\|M\|_{1}$

holds for any $M\in\mathbb{R}^{m_{1}\times m_{2}}$. With the identity

$\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},M\right\rangle=\langle\Delta_{j},M\rangle+\mathcal{C}_{j}(M_{j}^{0},M),$

we have

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda},\widetilde{M}_{t}^{\lambda})-2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M_{j}^{0},\widetilde{M}_{t}^{\lambda})$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M,M)-2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M_{j}^{0},M)+2\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle+\lambda(\|M\|_{1}-\|\widetilde{M}_{t}^{\lambda}\|_{1}).$

Using $\langle\cdot,*\rangle\leq\|\cdot\|_{\infty}\|*\|_{1}$, one can obtain

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})+2\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\|\widetilde{M}_{t}^{\lambda}-M\|_{1}+\lambda(\|M\|_{1}-\|\widetilde{M}_{t}^{\lambda}\|_{1}).$

Under the assumption that
$2\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\leq\lambda$, we have

$\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})\leq\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})+2\lambda\|M\|_{1}.$

Set $M=M_{t}^{0}$, then

$\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})\leq\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M_{t}^{0}-M_{j}^{0},M_{t}^{0}-M_{j}^{0})+2\lambda\|M_{t}^{0}\|_{1}.$

With the identity

$\displaystyle\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})$
$\displaystyle=$
$\displaystyle\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+\mathcal{C}_{j}(M_{t}^{0}-M_{j}^{0},M_{t}^{0}-M_{j}^{0})+2\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},M_{t}^{0}-M_{j}^{0}),$

we have

$\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},M_{t}^{0}-M_{j}^{0})\leq
2\lambda\|M_{t}^{0}\|_{1}.$ (S.2)

If
$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}<2\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}$,
we can maintain the bound

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}<2\delta_{h}M_{t}.$
(S.3)

If $\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\geq
2\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}$, with
Assumption 1, we know that

$\mathcal{C}_{j}(M,N)=\mathop{\rm vec}({M})^{\top}\Sigma\mathop{\rm
vec}({N}).$ (S.4)

So (S.2) can be expressed as

$\displaystyle\left(\mathop{\rm vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)^{\top}\Sigma\left[\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)+\left(2\mathop{\rm
vec}({M}_{t}^{0})-2\sum_{j=1}^{T}\omega_{j}\mathop{\rm
vec}({M}_{j}^{0})\right)\right]$ $\displaystyle\leq$ $\displaystyle
2\lambda\|M_{t}^{0}\|_{1}.$

From Lemma S.1,

$\displaystyle\left(\mathop{\rm vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)^{\top}\Sigma\left[\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)+\left(2\mathop{\rm
vec}({M}_{t}^{0})-2\sum_{j=1}^{T}\omega_{j}\mathop{\rm
vec}({M}_{j}^{0})\right)\right]$ (S.5) $\displaystyle\geq$
$\displaystyle\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}.$

Thus we have

$\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}-2\lambda\|M_{t}^{0}\|_{1}\leq
0.$

By solving this inequality, we have

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{\mu\delta_{h}M_{t}+\sqrt{\mu^{2}(\delta_{h}M_{t})^{2}+8\mu\lambda\|M_{t}^{0}\|_{1}}}{2\mu}\leq\delta_{h}M_{t}+\sqrt{\frac{2\lambda\|M_{t}^{0}\|_{1}}{\mu}}.$
(S.6)

Combining (S.3) and (S.6),

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\max\left\\{\delta_{h}M_{t},\sqrt{\frac{2\lambda\|M_{t}^{0}\|_{1}}{\mu}}\right\\}\leq\delta_{h}M_{t}+\left((\delta_{h}M_{t})^{2}+2\mu^{-1}\lambda\|M_{t}^{0}\|_{1}\right)^{1/2}.$
(S.7)

Because $\widetilde{M}_{t}^{\lambda}$ is the minimizer of (6), there exists a
sub-gradient matrix $B$ satisfying that for all $M\in\mathbb{R}^{m_{1}\times
m_{2}}$,

$\left\langle B,\widetilde{M}_{t}^{\lambda}-M\right\rangle\leq 0,$

which equals to that there exists
$\widehat{V}_{t}\in\partial\|\widetilde{M}_{t}^{\lambda}\|_{1}$ satisfying
that

$2\sum_{j=1}^{T}\omega_{j}\left(\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda},\widetilde{M}_{t}^{\lambda}-M)-\left\langle\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji},\widetilde{M}_{t}^{\lambda}-M\right\rangle\right)+\lambda\langle\widehat{V}_{t},\widetilde{M}_{t}^{\lambda}-M\rangle\leq
0.$

It is equivalent to

$\displaystyle
2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M)+\lambda\langle\widehat{V}_{t}-V,\widetilde{M}_{t}^{\lambda}-M\rangle$
$\displaystyle\leq$ $\displaystyle-\lambda\langle
V,\widetilde{M}_{t}^{\lambda}-M\rangle+2\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle,$

where $V\in\partial\|M\|_{1}$.
Denote that $M=\sum_{j=1}^{r}\sigma_{j}u_{j}v_{j}^{\top}$ with support
$S_{1},S_{2}$. $V$ have representation from Watson (1992):

$V=\sum_{j=1}^{r}u_{j}v_{j}^{\top}+P_{S_{1}^{\perp}}WP_{S_{2}^{\perp}},\
\|W\|_{\infty}\leq 1.$

We choose a $W$ subject to

$\langle
P_{S_{1}^{\perp}}WP_{S_{2}^{\perp}},\widetilde{M}_{t}^{\lambda}-M\rangle=\langle
W,P_{S_{1}^{\perp}}\widehat{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\rangle=\|P_{S_{1}^{\perp}}\widehat{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\|_{1}.$

With the monotonicity of sub-gradients, we have
$\langle\widehat{V}_{t}-V,\widetilde{M}_{t}^{\lambda}-M\rangle\geq 0$. So

$\displaystyle
2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M)+\lambda\|P_{S_{1}^{\perp}}\widetilde{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\|_{1}$
$\displaystyle\leq$
$\displaystyle-\lambda\langle\sum_{j=1}^{r}u_{j}v_{j}^{\top},\widetilde{M}_{t}^{\lambda}-M\rangle+2\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle.$

Using the identity

$\displaystyle
2\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\left(\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})+\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M,\widetilde{M}_{t}^{\lambda}-M)-\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})\right)$

and the inequality

$-\langle\sum_{j=1}^{r}u_{j}v_{j}^{\top},\widetilde{M}_{t}^{\lambda}-M\rangle=-\langle\sum_{j=1}^{r}u_{j}v_{j}^{\top},P_{S_{1}}(\widetilde{M}_{t}^{\lambda}-M)P_{S_{2}}\rangle\leq\|P_{S_{1}}(\widetilde{M}_{t}^{\lambda}-M)P_{S_{2}}\|_{1},$

we have

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})+\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M,\widetilde{M}_{t}^{\lambda}-M)+\lambda\|P_{S_{1}^{\perp}}\widetilde{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\|_{1}$
(S.8) $\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})+\lambda\|P_{S_{1}}(\widetilde{M}_{t}^{\lambda}-M)P_{S_{2}}\|_{1}+2\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle.$

With identity

$\displaystyle\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle$
$\displaystyle=\langle
P_{M}(\Delta_{j}),\widehat{M}^{\lambda}_{t}-M\rangle+\langle
P_{S_{1}^{\perp}}\Delta_{j}P_{S_{2}^{\perp}},\widehat{M}^{\lambda}_{t}-M\rangle$
$\displaystyle=\langle
P_{M}(\Delta_{j}),\widehat{M}^{\lambda}_{t}-M\rangle+\langle
P_{S_{1}^{\perp}}\Delta_{j}P_{S_{2}^{\perp}},\widehat{M}^{\lambda}_{t}\rangle,$

we have that

$\displaystyle\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle$
$\displaystyle=$ $\displaystyle\left\langle
P_{M}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right),\widehat{M}^{\lambda}_{t}-M\right\rangle+\left\langle
P_{S_{1}^{\perp}}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)P_{S_{2}^{\perp}},\widehat{M}^{\lambda}_{t}\right\rangle$
$\displaystyle\leq$
$\displaystyle\left\|P_{M}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)\right\|_{2}\left\|\widetilde{M}_{t}^{\lambda}-M\right\|_{2}+\left\|P_{S_{1}^{\perp}}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)P_{S_{2}^{\perp}}\right\|_{\infty}\left\|P_{S_{1}^{\perp}}\widetilde{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\right\|_{1}.$

With

$\displaystyle\left\|P_{M}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)\right\|_{2}$
$\displaystyle\leq\sqrt{\mathop{\rm
rank}\left[P_{A}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)\right]}\left\|P_{M}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)\right\|_{\infty}$
$\displaystyle\leq\sqrt{2\mathop{\rm
rank}(M)}\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\ $

and

$\left\|P_{S_{1}^{\perp}}\left(\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right)P_{S_{2}^{\perp}}\right\|_{\infty}\leq\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty},$

we have

$\sum_{j=1}^{T}\omega_{j}\langle\Delta_{j},\widetilde{M}_{t}^{\lambda}-M\rangle\leq\sqrt{2\mathop{\rm
rank}(M)}\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\|\widetilde{M}_{t}^{\lambda}-M\|_{2}+\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\|P_{S_{1}^{\perp}}\widetilde{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\|_{1}.$
(S.9)

Also we have

$\displaystyle\|P_{S_{1}}(\widetilde{M}_{t}^{\lambda}-M)P_{S_{2}}\|_{1}$
$\displaystyle\leq\sqrt{\mathop{\rm
rank}(M)}\|P_{S_{1}}(\widetilde{M}_{t}^{\lambda}-M)P_{S_{2}}\|_{2}$ (S.10)
$\displaystyle\leq\sqrt{\mathop{\rm
rank}(M)}\|\widetilde{M}_{t}^{\lambda}-M\|_{2}.$

With (LABEL:eq_proof_error_bound_1), (S.9) and (S.10),

$\displaystyle\sum_{j=1}^{T}\omega_{j}\left[\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})+\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M,\widetilde{M}_{t}^{\lambda}-M)\right]$
$\displaystyle+$
$\displaystyle\left(\lambda-2\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\right)\|P_{S_{1}^{\perp}}\widetilde{M}_{t}^{\lambda}P_{S_{2}^{\perp}}\|_{1}$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})+\sqrt{\mathop{\rm
rank}(M)}\left(\lambda+2\sqrt{2}\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}\right)\|\widetilde{M}_{t}^{\lambda}-M\|_{2}.$

With $\lambda\geq 2\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\|_{\infty}$,

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})+\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M,\widetilde{M}_{t}^{\lambda}-M)$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M-M_{j}^{0},M-M_{j}^{0})+\lambda(1+\sqrt{2})\sqrt{\mathop{\rm
rank}(M)}\|\widetilde{M}_{t}^{\lambda}-M\|_{2}.$

Set $M=M_{t}^{0}$, so

$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})+\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{T}\omega_{j}\mathcal{C}_{j}(M_{t}^{0}-M_{j}^{0},M_{t}^{0}-M_{j}^{0})+\lambda(1+\sqrt{2})\sqrt{\mathop{\rm
rank}(M)}\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}.$

With the identity

$\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{j}^{0},\widetilde{M}_{t}^{\lambda}-M_{j}^{0})=\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+2\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M^{0}_{t},M_{t}^{0}-M_{j}^{0})+\mathcal{C}_{j}(M_{t}^{0}-M_{j}^{0},M_{t}^{0}-M_{j}^{0}),$

we have

$\displaystyle\sum_{j=1}^{T}\omega_{j}\left[\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M^{0}_{t},M_{t}^{0}-M_{j}^{0})\right]$
$\displaystyle\leq$ $\displaystyle\frac{1+\sqrt{2}}{2}\lambda\sqrt{\mathop{\rm
rank}(M_{t}^{0})}\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}.$

If
$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}<\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}$,
we maintain the bound

$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}<\delta_{h}M_{t}.$ (S.11)

If
$\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\|_{2}\geq\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}$,
using (S.4) and Lemma S.1, we know that

$\displaystyle\sum_{j=1}^{T}\omega_{j}\left[\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M^{0}_{t},M_{t}^{0}-M_{j}^{0})\right]$
(S.12) $\displaystyle=$ $\displaystyle\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)^{\top}\Sigma\left(\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)+\left(M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right)\right)$
$\displaystyle\geq$
$\displaystyle\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}.$

So we have

$\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\mu\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}\leq\frac{1+\sqrt{2}}{2}\lambda\sqrt{\mathop{\rm
rank}(M_{t}^{0})}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}.$
(S.13)

Combining (S.11) and (S.13), we have

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\frac{1+\sqrt{2}}{2}\frac{\lambda}{\mu}\sqrt{r_{t}},$
(S.14)

where $r_{t}=\mathop{\rm rank}(M_{t}^{0})$. With (S.7) and (S.14), the proof
is completed.

$\square$

##  S.2 Proofs of Theorem 2 and Corollary 1, 2

###### Proposition 1

Under Assumption 2 and 3, if $h\rightarrow 0$ and
$(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as
$m_{1},m_{2},T\rightarrow\infty$, then

$\delta_{h}M_{t}\leq\frac{1}{2}\alpha(K)(m_{1}m_{2})^{1/2}D_{2}h^{2}+o\left((m_{1}m_{2})^{1/2}h^{2}\right).$

Proof of Proposition 1: Note that $T^{-1}\ll(m_{1}m_{2})^{-1/2}h^{2}$. The
conclusion follows from the classical nonparametric results (of fixed design)
that

$\displaystyle\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}$
$\displaystyle=\left\|\sum_{j=1}^{T}\frac{1}{Th}K\left(\frac{j-t}{Th}\right)\left\\{M\left(\frac{j}{T}\right)-M\left(\frac{t}{T}\right)\right\\}\right\|_{2}$
$\displaystyle=\left\|\int\frac{1}{h}K\left(\frac{u-t/T}{h}\right)\left\\{M(u)-M\left(\frac{t}{T}\right)\right\\}du\right\|_{2}+O\left(T^{-1}\right)$
$\displaystyle=\left\|\int
K(z)\left\\{M\left(zh+\frac{t}{T}\right)-M\left(\frac{t}{T}\right)\right\\}dz\right\|_{2}+O\left(T^{-1}\right)$
$\displaystyle=\frac{1}{2}\alpha(K)\left\|\nabla^{2}M(\frac{t}{T})\right\|_{2}h^{2}+o\left((m_{1}m_{2})^{1/2}h^{2}\right)$
$\displaystyle\leq\frac{1}{2}\alpha(K)(m_{1}m_{2})^{1/2}D_{2}h^{2}+o\left((m_{1}m_{2})^{1/2}h^{2}\right).$

$\square$

Now we prove two lemmas which are needed to finish the proofs Lemma S.4 and
Theorem 2. For brevity, define that
$\|\cdot\|\overset{\triangle}{=}\|\cdot\|_{\infty}$. We extend proposition 2
in Koltchinskii (2011b) to obtain the following Lemma S.2.

###### Lemma S.2

Let $Z_{1},Z_{2},\dots,Z_{n}\in(\mathcal{B}(\mathbb{R}^{m_{1}\times
m_{2}}),\mathbb{R}^{m_{1}\times m_{2}},\mu^{\otimes m_{1}m_{2}})$ be
independent, mean-zero random variables satisfying that

$\|Z_{i}\|_{L^{2}}\leq\sigma,\quad\|Z_{i}\|_{\psi(\alpha)}\leq
K,\quad\|\cdot\|_{\psi(\alpha)}:=\inf\left\\{t>0,\mathbb{E}\left[\exp\left(\frac{\|\cdot\|^{\alpha}}{t^{\alpha}}\right)\right]\leq
2\right\\}.$

Let $\alpha>1$ and $a=(a_{1},a_{2},\dots,a_{n})^{\top}\in\mathbb{R}^{n}$.
There exists a constant $C>0$ such that, for all $t>0$, with probability at
least $1-e^{-t}$

$\left\|\sum_{i=1}^{n}a_{i}Z_{i}\right\|\leq
C\max\left(\sqrt{\|a\|_{2}^{2}\sigma^{2}(t+\log(2m))},\|a\|K\left(\log\frac{K}{\sigma}\right)^{1/\alpha}(t+\log(2m))\right).$
(S.15)

Proof of Lemma S.2: We follow the proof of proposition 2 in Koltchinskii
(2011b) here.
Let

$S_{n}:=\sum_{i=1}^{n}a_{i}\widetilde{Z}_{i},\quad\widetilde{Z}=\left(\begin{matrix}0&Z\\\
Z^{\top}&0\\\ \end{matrix}\right)\in\mathbb{R}^{m\times m}.$

We know that $\|S_{n}\|<s$ if and only if $-sI_{m}\prec S_{n}\prec sI_{m}$.
Therefore,

$P(\|S_{n}\|\geq s)=P(S_{n}\not\preceq sI_{m})+P(S_{n}\not\succeq-sI_{m})$

and

$P(S_{n}\not\preceq sI_{m})=P\left(e^{\lambda S_{n}}\not\preceq e^{\lambda
sI_{m}}\right)\leq P\left(\operatorname{Tr}(e^{\lambda S_{n}})\geq e^{\lambda
s}\right)\leq e^{-\lambda s}\mathbb{E}\operatorname{Tr}(e^{\lambda S_{n}}).$

To bound $\mathbb{E}\operatorname{Tr}(e^{\lambda S_{n}})$, we use independence
and Golden-Thompson inequality

$\displaystyle\mathbb{E}\operatorname{Tr}(e^{\lambda S_{n}})$
$\displaystyle=\mathbb{E}\operatorname{Tr}\left(e^{\lambda S_{n-1}+\lambda
a_{n}Z_{n}}\right)$
$\displaystyle\leq\mathbb{E}\operatorname{Tr}\left(e^{\lambda
S_{n_{1}}}e^{\lambda a_{n}Z_{n}}\right)$
$\displaystyle=\operatorname{Tr}\left(\mathbb{E}e^{\lambda
S_{n-1}}\mathbb{E}e^{\lambda a_{n}Z_{n}}\right)$
$\displaystyle\leq\mathbb{E}\operatorname{Tr}\left(e^{\lambda
S_{n-1}}\right)\left\|\mathbb{E}e^{\lambda a_{n}Z_{n}}\right\|.$

By induction, we have

$\mathbb{E}\operatorname{Tr}(e^{\lambda S_{n}})\leq
m\prod_{i=1}^{n}\left\|\mathbb{E}e^{\lambda a_{i}Z_{i}}\right\|.$

It remains to bound the norm $\|\mathbb{E}e^{\lambda a_{i}Z_{i}}\|$. Using
Taylor expansion and $\mathbb{E}(Z_{i})=0$, we have

$\displaystyle\mathbb{E}e^{\lambda a_{i}Z_{i}}$
$\displaystyle=I_{m}+\mathbb{E}\lambda^{2}a_{i}^{2}Z_{i}^{2}\left[\frac{1}{2!}+\frac{\lambda
a_{i}Z_{i}}{3!}+\frac{\lambda^{2}a_{i}^{2}Z_{i}^{2}}{4!}+\ldots\right]$
$\displaystyle\leq
I_{m}+\lambda^{2}a_{i}^{2}\mathbb{E}Z_{i}^{2}\left[\frac{1}{2!}+\frac{\lambda
a_{i}\|Z_{i}\|}{3!}+\frac{\lambda^{2}a_{i}^{2}\|Z_{i}\|^{2}}{4!}+\ldots\right]$
$\displaystyle=I_{m}+\lambda^{2}a_{i}^{2}\mathbb{E}Z_{i}^{2}\left[\frac{e^{\lambda
a_{i}\|Z_{i}\|}-1-\lambda
a_{i}\|Z_{i}\|}{\lambda^{2}a_{i}^{2}\|Z_{i}\|^{2}}\right].$

Therefore, for all $\tau>0$ we have

$\displaystyle\left\|\mathbb{E}e^{\lambda a_{i}Z_{i}}\right\|$
$\displaystyle\leq
1+\lambda^{2}a_{i}^{2}\left\|\mathbb{E}Z_{i}^{2}\left[\frac{e^{\lambda
a_{i}\|Z_{i}\|}-1-\lambda
a_{i}\|Z_{i}\|}{\lambda^{2}a_{i}^{2}\|Z_{i}\|^{2}}\right]\right\|$
$\displaystyle\leq
1+\lambda^{2}\left\|\mathbb{E}a_{i}^{2}Z_{i}^{2}\right\|\left[\frac{e^{\lambda\tau}-1-\lambda\tau}{\lambda^{2}\tau^{2}}\right]+\lambda^{2}\mathbb{E}\|a_{i}Z_{i}\|^{2}\left[\frac{e^{\lambda
a_{i}\|Z_{i}\|}-1-\lambda
a_{i}\|Z_{i}\|}{\lambda^{2}a_{i}^{2}\|Z_{i}\|^{2}}\right]I(\|a_{i}Z_{i}\|\geq\tau).$

Let

$\tau=C_{2}\|a_{i}\|_{\infty}K(\log\frac{K}{\sigma})^{1/\alpha},\quad\lambda\tau\leq
1,$

we have

$\left\|\mathbb{E}e^{\lambda a_{i}Z_{i}}\right\|\leq
1+C_{1}\lambda^{2}a_{i}^{2}\sigma^{2}\leq\exp\\{C_{1}\lambda^{2}a_{i}^{2}\sigma^{2}\\}.$

So when

$P(\|S_{n}\|\geq s)\leq 2m\exp\\{-\lambda
s+C_{1}\lambda^{2}\|a\|_{2}^{2}\sigma^{2}\\}$

and $\lambda$ is chosen to be

$\min\left\\{\frac{1}{C_{2}\|a\|_{\infty}K\left(\log
K/\sigma\right)^{1/\alpha}},\frac{s}{2C_{1}\|a\|_{2}^{2}\sigma^{2}}\right\\},$

which can directly deduce (S.15).

$\square$

###### Lemma S.3

If X,Y are random variables with $\|X\|_{\psi(\alpha)}\leq K_{1}$,
$\|Y\|_{\psi(\beta)}\leq K_{2}$, and there exists $\gamma\geq 1$ that
$\frac{1}{\alpha}+\frac{1}{\beta}=\frac{1}{\gamma}\leq 1$, then we have

$\|XY\|_{\psi(\gamma)}\leq K_{1}K_{2}.$ (S.16)

Proof of Lemma S.3: Using Young’s inequality we know that for all
$\theta\in(0,1)$

$\frac{\left(\|X\|^{\gamma}\right)^{\alpha/\gamma}}{t^{\alpha\theta}\alpha/\gamma}+\frac{\left(\|Y\|^{\gamma}\right)^{\beta/\gamma}}{t^{\beta(1-\theta)}\beta/\gamma}\geq\frac{\left(\|XY\|\right)^{\gamma}}{t^{\gamma}}.$

Thus we have

$\displaystyle\mathbb{E}\left[\exp\left(\frac{\|XY\|^{\gamma}}{t^{\gamma}}\right)\right]$
$\displaystyle\leq\mathbb{E}\left[\exp\left(\frac{\left(\|X\|^{\gamma}\right)^{\alpha/\gamma}}{t^{\alpha\theta}\alpha/\gamma}+\frac{\left(\|Y\|^{\gamma}\right)^{\beta/\gamma}}{t^{\beta(1-\theta)}\beta/\gamma}\right)\right]$
$\displaystyle\leq\frac{\gamma}{\alpha}\mathbb{E}\left[\exp\left(\frac{(\|X\|^{\alpha}}{t^{\alpha\theta}}\right)\right]+\frac{\gamma}{\beta}\mathbb{E}\left[\exp\left(\frac{(\|Y\|^{\beta}}{t^{\beta(1-\theta)}}\right)\right]$
$\displaystyle\leq 2,$

where we choose $t,\theta$ by $t^{\theta}=K_{1},t^{1-\theta}=K_{2}$.
So we have

$\mathbb{E}\left[\exp\left(\frac{\|XY\|^{\gamma}}{(K_{1}K_{2})^{\gamma}}\right)\right]\leq
2.$

$\square$

We construct the following theorem to bound
$\|\sum_{j=1}^{T}\omega_{h}(j-t)\Delta_{j}\|$ and Lemma S.4 is an immediate
result of the following theorem. For brevity, we denote
$\omega_{j}=\omega_{h}(j-t)$.

###### Theorem S.1

Under Assumption 4-5, for all $t>0$, with probability at least $1-3e^{-t}$,

$\displaystyle\left\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right\|$
$\displaystyle\leq
C_{1}\max\left(\sigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{K\left(\log\frac{K}{\sigma}\right)^{1/\gamma}(t+\log(2m))}{W}\right)$
$\displaystyle+C_{2}\mu_{X}\max\left(\sigma_{\xi}\sqrt{\frac{t}{N}},\frac{K_{1}(\log\frac{K_{1}}{\sigma_{\xi}})^{1/\alpha}t}{W}\right)$
$\displaystyle+C_{3}\max\left(\varsigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{\mathcal{K}\left(\log\frac{\mathcal{K}}{\varsigma}\right)^{2/\beta}(t+\log(2m))}{W}\right),$

where

$N=\frac{1}{\sum_{j=1}^{T}\omega_{j}^{2}/n_{j}},\quad
W=\frac{1}{\max_{j}\omega_{j}/n_{j}}$

and constants with respect to $\omega_{j},n_{j}$ such that

$\displaystyle
K=\max\|\xi_{ji}(X_{ji}-\mathbb{E}X_{ji})\|_{\psi(\gamma)}\lesssim
K_{1}K_{2},$
$\displaystyle\sigma=\max\|\xi_{ji}(X_{ji}-\mathbb{E}X_{ji})\|_{L^{2}}\lesssim
K_{1}K_{2}2^{1/\gamma},$
$\displaystyle\mu_{X}=\max\|\mathbb{E}X_{ji}\|\lesssim K_{2},$
$\displaystyle\sigma_{\xi}=\max\|\xi_{ji}\|_{L^{2}}\lesssim
K_{1}2^{1/\alpha},$ $\displaystyle\mathcal{K}=\max\|\langle
M_{j}^{0},X_{ji}\rangle X_{ji}-\mathbb{E}(\langle M^{0}_{j},X_{ji}\rangle
X_{ji})\|_{\psi(\beta/2)}\lesssim C_{*}K_{2}^{2},$
$\displaystyle\varsigma=\max\|\langle M_{j}^{0},X_{ji}\rangle
X_{ji}-\mathbb{E}(\langle M^{0}_{j},X_{ji}\rangle X_{ji})\|_{L^{2}}\lesssim
C_{*}K_{2}^{2}2^{2/\beta}.$

Proof of Theorem S.1: Because that

$\displaystyle\left\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right\|$
$\displaystyle=\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\left(Y_{ji}X_{ji}-EY_{ji}X_{ji}\right)\right\|$
(S.17)
$\displaystyle\leq\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\xi_{ji}X_{ij}\right\|$
$\displaystyle+\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\left\\{\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}-E\left(\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}\right)\right\\}\right\|$
$\displaystyle\leq\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\xi_{ji}(X_{ij}-EX_{ij})\right\|+\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\xi_{ji}\mathbb{E}X_{ij}\right\|$
$\displaystyle+\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{1}{n_{j}}\omega_{j}\left\\{\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}-E\left(\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}\right)\right\\}\right\|$
$\displaystyle:=\mathcal{D}_{1}+\mathcal{D}_{2}+\mathcal{D}_{3}.$

  * a.

For $\mathcal{D}_{1}$.

We denote that $Z_{ji}:=\xi_{ji}(X_{ji}-\mathbb{E}X_{ji})$. Using Lemma S.3,
$Z_{ji}$ are independent, mean-zero random variables satisfying that

$\|Z_{ji}\|_{\psi(\gamma)}\leq
K,\quad\|Z_{ji}\|_{L^{2}}\leq\sigma(2)^{1/\gamma}.$

Set $a_{ji}={\omega_{j}}/{n_{j}}$. Using Lemma S.2, we have that for all
$t>0$, with probability at least $1-e^{-t}$

$\mathcal{D}_{1}=\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}a_{ji}Z_{ji}\right\|\leq
C\max\left(\sigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{K\left(\log\frac{K}{\sigma}\right)^{1/\gamma}(t+\log(2m))}{W}\right).$
(S.18)

  * b.

For $\mathcal{D}_{2}$.

Directly using Lemma S.2 with
$\|\mathbb{E}X_{ji}\|\leq\mu_{X},\|\xi_{ji}\|_{L^{2}}\leq\sigma_{\xi}$ and
$a_{ji}=\frac{\omega_{j}}{n_{j}}$ , we know that for all $t\geq 0$, with
probability at least $1-e^{-t}$

$\mathcal{D}_{2}=\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}a_{ji}\mu_{X}\xi_{ji}\right\|\leq
C\mu\max\left(\sigma_{\xi}\sqrt{\frac{t}{N}},\frac{K_{1}(\log\frac{K_{1}}{\sigma_{\xi}})^{1/\alpha}t}{W}\right).$
(S.19)

  * c.

For $\mathcal{D}_{3}$. Denote that $Z_{ji}=\langle M_{j}^{0},X_{ji}\rangle
X_{ji}-\mathbb{E}(\langle M^{0}_{j},X_{ji}\rangle X_{ji})$ and
$a_{ji}=\frac{\omega_{j}}{n_{j}}$ with

$\|Z_{ji}\|_{\psi(\beta/2)}\leq\mathcal{K},\quad\|Z_{ji}\|_{L^{2}}\leq\varsigma.$

Using Lemma S.2, we know that for all $t\geq 0$, with probability at least
$1-e^{-t}$

$\mathcal{D}_{3}=\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}a_{ji}Z_{ji}\right\|\leq
C\max\left(\varsigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{\mathcal{K}\left(\log\frac{\mathcal{K}}{\varsigma}\right)^{2/\beta}(t+\log(2m))}{W}\right).$
(S.20)

From (S.18), (S.19), (S.20), we know that for all $t>0$, with probability at
least $1-3e^{-t}$,

$\displaystyle\left\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right\|$
$\displaystyle\leq
C_{1}\max\left(\sigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{K\left(\log\frac{K}{\sigma}\right)^{1/\gamma}(t+\log(2m))}{W}\right)$
$\displaystyle+C_{2}\mu_{X}\max\left(\sigma_{\xi}\sqrt{\frac{t}{N}},\frac{K_{1}(\log\frac{K_{1}}{\sigma_{\xi}})^{1/\alpha}t}{W}\right)$
$\displaystyle+C_{3}\max\left(\varsigma\sqrt{\frac{(t+\log(2m))}{N}},\frac{\mathcal{K}\left(\log\frac{\mathcal{K}}{\varsigma}\right)^{2/\beta}(t+\log(2m))}{W}\right).$

$\square$

###### Lemma S.4

Under Assumption 4 and 5, if $h\rightarrow 0$ and
$(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as
$m_{1},m_{2},T\rightarrow\infty$,
$nTh\gg(K_{*}/\sigma_{*})^{2}\log(m_{1}+m_{2})$, then with probability at
least $1-3/(m_{1}+m_{2})$,

$\mathcal{W}_{h}\Delta_{t}\leq
C_{1}\sigma_{*}\sqrt{\frac{\log(m_{1}+m_{2})}{nTh}},$ (S.21)

where $C_{1}>0$ is a constant independent to
$\sigma_{*},K_{*},T,h,m_{1},m_{2}$ such that

$\displaystyle\sigma_{*}$
$\displaystyle=\max\left\\{\sigma,\mu_{X}\sigma_{\xi},\varsigma\right\\},$
$\displaystyle K_{*}$
$\displaystyle=\max\left\\{\mathcal{K}\left(\log\frac{\mathcal{K}}{\varsigma}\right)^{2/\beta},\mu_{X}K_{1}\left(\log\frac{K_{1}}{\sigma_{\xi}}\right)^{1/\alpha},K\left(\log\frac{K}{\sigma}\right)^{1/\gamma}\right\\},$

in which $\mu_{X}=\max_{i,j}\|\mathbb{E}X_{ji}\|_{\infty},\
\sigma_{\xi}=\max_{i,j}\|\xi_{ji}\|_{L^{2}}$ and

$\displaystyle
K=\max_{i,j}\left\|\xi_{ji}(X_{ji}-\mathbb{E}X_{ji})\right\|_{\psi(\gamma)},\quad\mathcal{K}=\max_{i,j}\|\langle
M_{j}^{0},X_{ji}\rangle X_{ji}-\mathbb{E}(\langle M^{0}_{j},X_{ji}\rangle
X_{ji})\|_{\psi(\beta/2)},$
$\displaystyle\sigma=\max_{i,j}\|\xi_{ji}(X_{ji}-\mathbb{E}X_{ji})\|_{L^{2}},\quad\varsigma=\max_{i,j}\|\langle
M_{j}^{0},X_{ji}\rangle X_{ji}-\mathbb{E}(\langle M^{0}_{j},X_{ji}\rangle
X_{ji})\|_{L^{2}},$

with
$\|X\|_{L^{2}}=\|\mathbb{E}XX^{T}\|_{\infty}^{1/2}\vee\|\mathbb{E}X^{T}X\|_{\infty}^{1/2}$.

Proof of Lemma S.4:

Define that

$\mathcal{N}_{h}(t)=\\{j:t-Th<j<t+Th\\},\quad\bar{n}_{t}=\frac{1}{|\mathcal{N}_{h}(t)|}\sum_{j\in\mathcal{N}_{h}(t)}n_{j},\quad\check{n}_{t}=\min_{j\in\mathcal{N}_{h}(t)}n_{j}.$

Recall the definition of $\omega_{j}$ that
$\omega_{j}=\omega_{h}(j-t)=\frac{K\left(\frac{j-t}{Th}\right)}{\sum_{k=1}^{T}K\left(\frac{k-t}{Th}\right)},$
we have that

$\displaystyle N$
$\displaystyle=\frac{1}{\sum_{j=1}^{T}\omega_{h}(j-t)^{2}/n_{j}}=\frac{\left\\{\sum_{k=1}^{T}\frac{1}{Th}K\left(\frac{k-t}{Th}\right)\right\\}^{2}}{\sum_{j=1}^{T}\frac{1}{T^{2}h^{2}}K\left(\frac{j-t}{Th}\right)^{2}/n_{j}}$
$\displaystyle\asymp\frac{1}{\sum_{j=1}^{T}\frac{1}{n_{j}T^{2}h^{2}}K\left(\frac{j-t}{Th}\right)^{2}}\asymp\left\\{\frac{1}{\bar{n}_{t}Th}\sum_{j=1}^{T}\frac{1}{Th}K\left(\frac{j-t}{Th}\right)^{2}\right\\}^{-1}$
$\displaystyle=\bar{n}_{t}ThR(K)^{-1},$

and

$\displaystyle W$
$\displaystyle=\frac{1}{\max_{j}\omega_{h}(j-t)/n_{j}}=\left\\{\max_{j}\frac{\frac{1}{Th}K\left(\frac{j-t}{Th}\right)}{n_{j}\frac{1}{Th}\sum_{k=1}^{T}K\left(\frac{k-t}{Th}\right)}\right\\}^{-1}$
$\displaystyle\asymp\left\\{\max_{j}{\frac{1}{n_{j}Th}K\left(\frac{j-t}{Th}\right)}\right\\}^{-1}\asymp\check{n}_{t}ThK(0)^{-1}.$

With the assumption $\bar{n}_{t},\check{n}_{t}\asymp n$, for
$t\leq\log(m_{1}+m_{2})$, when
$nTh\gg\\{K_{*}/\sigma_{*}\\}^{2}\log(m_{1}+m_{2})$, the first terms in
(S.18), (S.19) and (S.20) dominate those bounds which gives

$\mathcal{W}_{h}\Delta_{t}\leq{C_{1}}\sigma_{*}\sqrt{\frac{s+\log(m_{1}+m_{2})}{nTh}}\ll
1.$

By setting $t=\log(m_{1}+m_{2})$, we can deduce Lemma S.4 directly.

$\square$

Proof of Theorem 2: With Lemma S.4 and Theorem 1, we have

$\left(\delta_{h}M_{t}^{2}+2\mu^{-1}\mathcal{W}_{h}\Delta_{t}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2}\gg{(1+\sqrt{2})}{\mu^{-1}}{\mathcal{W}_{h}\Delta_{t}}\sqrt{r_{t}},$

which leads to the following conclusion. From Lemma S.4, we know that for all
$t>0$, with probability at least $1-3e^{-t}$,

$\left\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right\|_{\infty}\leq{C_{1}}\max\left\\{\sigma^{*}\sqrt{\frac{t+\log(m_{1}+m_{2})}{\bar{n}_{t}Th}},K^{*}\frac{t+\log(m_{1}+m_{2})}{\check{n}_{t}Th}\right\\}.$

Choose that

$t=\log(m_{1}+m_{2}),\quad\lambda=2\sqrt{2}C_{1}\sigma^{*}\sqrt{\frac{\log(m_{1}+m_{2})}{\bar{n}_{t}Th}}.$

If
$(\check{n}_{t}^{2}/\bar{n}_{t})Th\gg\\{K^{*}/\sigma^{*}\\}^{2}\log(m_{1}+m_{2})$,
than

$\sigma^{*}\sqrt{\frac{t+\log(m_{1}+m_{2})}{\bar{n}_{t}Th}}\gg
K^{*}\frac{t+\log(m_{1}+m_{2})}{\check{n}_{t}Th}.$

From Proposition 1

$(m_{1}m_{2})^{-1/2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}\leq\frac{1}{2}\alpha(K)D_{2}h^{2}.$

With Theorem 1, so

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{1}{2}C_{2}\alpha(K)D_{2}h^{2}+\left(2+\sqrt{2}\right)C_{1}\sigma^{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}\bar{n}_{t}Th}\right)^{1/2}$

holds with probability at least $1-3/(m_{1}+m_{2})$.

$\square$

Proof of Corollary 1: Plug in the upper bound of $\mathcal{W}_{h}\Delta_{t}$
in Lemma S.4, Corollary 1 can be obtained immediately. $\square$

Proof of Corollary 2:

From the distribution of $X_{ji}$ and Assumption 4, we have that

$\gamma=\alpha,\quad\beta=\infty,\quad\mu=(m_{1}m_{2})^{-1},$

and

$\displaystyle\sigma\leq\sigma_{\xi}\sqrt{\frac{1}{m_{1}\wedge
m_{2}}},\quad\mu_{X}\sigma_{\xi}=\sigma_{\xi}\sqrt{\frac{1}{m_{1}\wedge
m_{2}}},\quad\varsigma\leq M\sqrt{\frac{1}{m_{1}\wedge m_{2}}},$
$\displaystyle K\asymp\sigma_{\xi},\quad
K_{1}\asymp\sigma_{\xi},\quad\mathcal{K}\asymp M.$

Then it can be obtained that

$\displaystyle\sigma^{*}=(C_{M}\vee\sigma_{\xi})\sqrt{\frac{1}{m_{1}\wedge
m_{2}}},$ $\displaystyle K^{*}/\sigma^{*}\asymp(m_{1}\wedge
m_{2})^{1/2}\log^{1/\alpha}(m_{1}\wedge m_{2}).$

Then $(K^{*}/\sigma^{*})^{2}\log(m_{1}+m_{2})\asymp(m_{1}\wedge
m_{2})\log^{2/\alpha}(m_{1}\wedge m_{2})\log(m_{1}+m_{2})$. Note that
$\log(m_{1}+m_{2})>\log(m_{1}\wedge m_{2})$, when
$(\check{n}_{t}^{2}/\bar{n}_{t})Th\gg(m_{1}\wedge
m_{2})\log^{1+2/\alpha}(m_{1}+m_{2})$, the first term in (S.21) dominates the
bound. Apply Theorem 2 and Corollary 1 and the proof is finished.

$\square$

##  S.3 Heterogeneous Scenario

Consider the heterogeneous assumption first.

###### Assumption S.1

(Heterogeneous Assumption) The second moment matrix $\Sigma_{t}$ of $X_{t}$ is
positive definite with smallest eigenvalue $\mu_{t}>0$ for $t=1,2,\cdots,T$.

Here we provide the corresponding results of Theorem 1, Theorem 2 and
Corollary 1 under the heterogeneous scenario. The key difference between
homogeneous and heterogeneous scenarios is that the covariance matrix used for
time points $t$ is $\Sigma$ in homogeneous case and
$\sum_{j=1}^{T}\omega_{h}(j-t)\Sigma_{j}$ in heterogeneous case. Define
$\nu_{t}=\sum_{j=1}^{T}\omega_{h}(j-t)\mu_{t}$, and the difference is using
$\nu_{t}$ to replace $\mu$ in each theorems and corollaries for the upper
bound of $\widetilde{M}_{t}^{\lambda}$.

###### Theorem S.2

(Heterogeneous Theorem 1) Under Assumption S.1, if $\lambda\geq
2\mathcal{W}_{h}\Delta_{t}$, then

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\min\left\\{\left((\delta_{h}M_{t})^{2}+\frac{2\lambda}{\nu_{t}}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2},\frac{1+\sqrt{2}}{2}\frac{\lambda}{\nu_{t}}\sqrt{r_{t}}\right\\}.$
(S.22)

When selecting $\lambda=2\mathcal{W}_{h}\Delta_{t}$, we have

$\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\delta_{h}M_{t}+\min\left\\{\left((\delta_{h}M_{t})^{2}+2\nu_{t}^{-1}\mathcal{W}_{h}\Delta_{t}\left\|M_{t}^{0}\right\|_{1}\right)^{1/2},{(1+\sqrt{2})}{\nu_{t}^{-1}}{\mathcal{W}_{h}\Delta_{t}}\sqrt{r_{t}}\right\\}.$
(S.23)

Proof In the Proof of Theorem 1, the whole proof is not changed but changing
(S.5) as

$\displaystyle\left(\mathop{\rm vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)^{\top}\Sigma\left[\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)+\left(2\mathop{\rm
vec}({M}_{t}^{0})-2\sum_{j=1}^{T}\omega_{j}\mathop{\rm
vec}({M}_{j}^{0})\right)\right]$ $\displaystyle\geq$
$\displaystyle\nu_{t}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\nu_{t}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}.$

and changing (LABEL:eq:homo2) as

$\displaystyle\sum_{j=1}^{T}\omega_{j}\left[\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M_{t}^{0},\widetilde{M}_{t}^{\lambda}-M_{t}^{0})+\mathcal{C}_{j}(\widetilde{M}_{t}^{\lambda}-M^{0}_{t},M_{t}^{0}-M_{j}^{0})\right]$
$\displaystyle=$ $\displaystyle\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)^{\top}\Sigma\left(\left(\mathop{\rm
vec}({\widehat{M}}_{t}^{\lambda})-\mathop{\rm
vec}({M}_{t}^{0})\right)+\left(M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right)\right)$
$\displaystyle\geq$
$\displaystyle\nu_{t}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}^{2}-\nu_{t}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\left\|M_{t}^{0}-\sum_{j=1}^{T}\omega_{j}M_{j}^{0}\right\|_{2}.$

$\square$

###### Theorem S.3

(Heterogeneous Theorem 2) Under Assumption 1-5, let
$nTh\gg(K_{*}/\sigma_{*})^{2}\log(m_{1}+m_{2})$, $h\rightarrow 0$ and
$(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as
$n,m_{1},m_{2},T\rightarrow\infty$, when

$\lambda=2C_{1}\sigma_{*}\sqrt{\frac{\log(m_{1}+m_{2})}{nTh}},$

where $C_{1}$ and $\sigma_{*}$ are defined in Lemma S.4 of Supplementary
Material, then with probability at least $1-3/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+\left(1+\sqrt{2}\right)C_{1}\sigma_{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\nu_{t}^{2}m_{1}m_{2}nTh}\right)^{1/2}+o(h^{2}),$
(S.24)

where $\alpha(K)$ is defined in (10) and $D_{2}$ in Assumption 2.

###### Corollary S.1

(Heterogeneous Corollary 1) Under assumptions of Theorem S.3, when

$nT\gg\max\\{K_{*}^{5/2}\sigma_{*}^{-3}(\nu_{t}^{2}m_{1}m_{2})^{1/4}\log(m_{1}+m_{2}),\
\sigma_{*}^{2}(\nu_{t}^{2}m_{1}m_{2})^{-1}\log(m_{1}+m_{2})\\},$

with probability at least $1-3/(m_{1}+m_{2})$,

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{2}\left(\frac{\sigma_{*}^{2}r_{t}\log(m_{1}+m_{2})}{\nu_{t}^{2}m_{1}m_{2}nT}\right)^{2/5},$
(S.25)

where
$C_{2}=1/2\left[(2+2\sqrt{2})C_{1}\right]^{4/5}\left[\alpha(K)D_{2}\right]^{1/5}$.

##  S.4 Definition of $\phi-$mixing

The dependence of two $\sigma$-field $\mathcal{A}$ and $\mathcal{B}$ is
measured by

$\phi(\mathcal{A},\mathcal{B})=\sup\left\\{\left|P(A)-\frac{P(A\cap
B)}{P(B)}\right|;A\in\mathcal{A},B\in\mathcal{B},P(B)\not=0\right\\}.$

And for a sequence of $\sigma$-field $\mathcal{A}_{j},j\in\mathbb{N}^{*}$, the
$\phi$-coefficients $\phi_{\mathcal{A}}(k),k\in\mathbb{N}$ is defined as

$\phi_{\mathcal{A}}(k)=\sup_{|j-s|\leq
k}\phi(\mathcal{A}_{j},\mathcal{A}_{s}).$

If $\lim_{k\rightarrow\infty}\phi_{\mathcal{A}}(k)=0$, the sequence of
$\sigma$-field $\mathcal{A}_{j},j\in\mathbb{N}^{*}$ is said to be
$\phi$-mixing.

##  S.5 Proofs of Theorem 3 and Corollary 3, S.2

###### Lemma S.5

Under Assumption 4 and 6, if $h\rightarrow 0$ as $n,T\rightarrow\infty$,
$Th\gg\log(m_{1}+m_{2})$ and
$nTh\gg\left(K_{*}/\sigma_{*}\right)^{2}\log^{3}(m_{1}+m_{2})$, with
probability at least $1-3/(m_{1}+m_{2})$,

$\mathcal{W}_{h}\Delta_{t}\leq\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\sigma_{*}\sqrt{\frac{\log(m_{1}+m_{2})}{nTh}},$

where $\mathcal{C}_{1}>0$ is a constant independent to
$\sigma_{*},K_{*},T,h,m_{1},m_{2}$.

Proof of Lemma S.5 and Theorem 3: Denote
$\mathbf{X}_{j}=(X_{j1},X_{j2},\dots,X_{jn_{j}})^{\top},\mathbf{\xi}_{j}=(\xi_{j1},\xi_{j2},\dots,\xi_{jn_{j}})^{\top}$.
We know that
$\sigma(\mathbf{X}_{j}-\mathbb{E}\mathbf{X}_{j})\subseteq\mathcal{X}_{j}$.
Define $\mathcal{Z}_{j}$ such that
$\mathcal{Z}_{j}=\sigma((\mathbf{X}^{\top}_{j},\mathbf{\xi}^{\top}_{j})^{\top})$,
then

$\displaystyle\left|P(A_{x}\otimes A_{y})-\frac{P(A_{x}\otimes A_{y}\cap
B_{x}\otimes B_{y})}{P(B_{x}\otimes B_{y})}\right|$ $\displaystyle=$
$\displaystyle\left|P(A_{x})\left(P(A_{y})-\frac{P(A_{y}\cap
B_{y})}{P(B_{y})}\right)+\frac{P(A_{y}\cap
B_{y})}{P(B_{y})}\left(P(A_{x})-\frac{P(A_{x}\cap
B_{x})}{P(B_{x})}\right)\right|$ $\displaystyle\leq$
$\displaystyle\phi_{\mathcal{X}}(|j-l|)+\phi_{\mathcal{Y}}(|j-l|)-\phi_{\mathcal{X}}(|j-l|)\phi_{\mathcal{Y}}(|j-l|),$

where $A_{x}\otimes A_{y}\in\mathcal{Z}_{j},\ B_{x}\otimes
B_{y}\in\mathcal{Z}_{l}$. So we know that
$\phi_{\mathcal{Z}}\leq\phi_{\mathcal{X}}+\phi_{\mathcal{Y}}-\phi_{\mathcal{X}}\phi_{\mathcal{Y}}$,
which means that

$\Phi_{\mathcal{Z}}=\sum_{k=1}^{\infty}\sqrt{\phi_{\mathcal{Z}}(k)}\leq\Phi_{\mathcal{X}}+\Phi_{\mathcal{Y}}<\infty.$

Let
$Z_{ji}=\frac{\omega_{j}}{n_{j}}\xi_{ji}\left(X_{ji}-\mathbb{E}X_{ji}\right)$
and $Z_{j}=\sum_{i=1}^{n_{j}}Z_{ji}$. Using Lemma S.2, we know that

$P\left(\left\|Z_{j}\right\|\geq s\right)\leq
2m\exp\left\\{-C\min\left\\{\frac{n_{j}s^{2}}{w^{2}_{j}\sigma^{2}},\frac{n_{j}s}{\omega_{j}K(\log\frac{K}{\sigma})^{1/\gamma}}\right\\}\right\\}.$

Denote the bounded function

$g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})=\frac{\omega_{j}}{n_{j}}(\mathbf{X}_{j}-\mathbb{E}\mathbf{X}_{j})^{\top}\mathbf{\xi}_{j}\mathbf{1}_{\left\|\frac{\omega_{j}}{n_{j}}(\mathbf{X}_{j}-\mathbb{E}\mathbf{X}_{j})^{\top}\mathbf{\xi}_{j}\right\|\leq
s},$

we know that

$\displaystyle P\left(\left\|\sum_{j=1}^{T}Z_{j}\right\|\geq t\right)$
$\displaystyle\leq\sum_{j=1}^{T}P(\|Z_{j}\|\geq
s)+P\left(\left\|\sum_{j=1}^{T}Z_{j}\mathbf{1}_{\|Z_{j}\|\leq s}\right\|\geq
t\right)$ $\displaystyle=\sum_{j=1}^{T}P(\|Z_{j}\|\geq
s)+P\left(\left\|\sum_{j=1}^{\top}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\geq
t\right).$

With the tail distribution assumptions for $X_{j},\xi_{j}$, when
$s\leq\frac{\omega_{j}\sigma^{2}}{K(\log\frac{K}{\sigma})^{1/\gamma}}$, we
have the inequalities

$\displaystyle\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|$
$\displaystyle=\left\|\sum_{j=1}^{T}Z_{j}\mathbf{1}_{\|Z_{j}\|\geq s}\right\|$
$\displaystyle\leq\sum_{j=1}^{T}\left|\int_{s}^{\infty}P(\|Z_{j}\|\geq
s)dx+sP(\|Z_{j}\|\geq s)\right|$
$\displaystyle\leq\left(s+\frac{c}{snT^{2}h^{2}}\right)\sum_{j=1}^{T}P(\|Z_{j}\|\geq
s).$

And when
$s\geq\frac{\omega_{j}\sigma^{2}}{K(\log\frac{K}{\sigma})^{1/\gamma}}$, we
have that

$\displaystyle\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\leq\left(s+\frac{cK(\log\frac{K}{\sigma})^{1/\gamma}}{nTh}\right)\sum_{j=1}^{T}P(\|Z_{j}\|\geq
s).$

Meanwhile we have the bound

$\left\|\mathbb{E}\sum_{j=1}^{T}g^{2}((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\leq\sum_{j=1}^{T}\frac{\omega_{j}^{2}}{n_{j}}\sigma_{\xi}^{2}\sigma_{X}^{2}=\frac{\sigma^{2}}{N}.$

Similar to Theorem 3 in Samson (2000), with the matrix value function
$g(\cdot)$ and norm $\|\cdot\|_{\infty}$

$\displaystyle
P\left(\left\|\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\geq
t\right)$ $\displaystyle\leq$ $\displaystyle
P\left(\left\|\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})-\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\geq
t-\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|\right)$
$\displaystyle\leq$ $\displaystyle
2m\exp\left\\{-\frac{1}{8\Phi_{Z}^{2}}\min\left\\{\frac{t}{2s},\frac{Nt^{2}}{16\sigma^{2}}\right\\}\right\\},$

when
$t\geq\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|$.
Thus we have

$\displaystyle P\left(\left\|\sum_{j=1}^{T}Z_{j}\right\|\geq t\right)$
$\displaystyle\leq$ $\displaystyle\sum_{j=(t-h/2)\vee 0}^{(t+h/2)\wedge
T}2m\exp\left\\{-C\min\left\\{\frac{n_{j}s^{2}}{w^{2}_{j}\sigma^{2}},\frac{n_{j}s}{\omega_{j}K(\log\frac{K}{\sigma})^{1/\gamma}}\right\\}\right\\}$
$\displaystyle+$ $\displaystyle
2m\exp\left\\{-\frac{1}{8\Phi_{\mathcal{Z}}^{2}}\min\left\\{\frac{t}{2s},\frac{Nt^{2}}{16\sigma^{2}}\right\\}\right\\}$
$\displaystyle\leq$ $\displaystyle\sum_{j=(t-h/2)\vee 0}^{(t+h/2)\wedge
T}2m\exp\left\\{-\frac{Cn_{j}s^{2}}{\omega_{j}^{2}\sigma^{2}}\right\\}+2m\exp\left\\{-\frac{1}{8\Phi_{\mathcal{Z}}^{2}}\min\left\\{\frac{t}{2s},\frac{Nt^{2}}{16\sigma^{2}}\right\\}\right\\}$
$\displaystyle\lesssim$ $\displaystyle
2m\exp\left\\{-\frac{CnT^{2}h^{2}s^{2}}{\sigma^{2}}+\log
Th\right\\}+2m\exp\left\\{-\frac{1}{8\Phi_{\mathcal{Z}}^{2}}\min\left\\{\frac{t}{2s},\frac{Nt^{2}}{16\sigma^{2}}\right\\}\right\\}$

when
$t\geq\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|$
and $s\leq\frac{\omega_{j}\sigma^{2}}{K(\log\frac{K}{\sigma})^{1/\gamma}}$.
And similarly

$\displaystyle P\left(\left\|\sum_{j=1}^{T}Z_{j}\right\|\geq t\right)$
$\displaystyle\lesssim
2m\exp\left\\{-C\frac{nThs}{K(\log\frac{K}{\sigma})^{1/\gamma}}+\log
Th\right\\}+2m\exp\left\\{-\frac{1}{8\Phi_{\mathcal{Z}}^{2}}\min\left\\{\frac{t}{2s},\frac{Nt^{2}}{16\sigma^{2}}\right\\}\right\\}$

when
$t\geq\left\|\mathbb{E}\sum_{j=1}^{T}g((\mathbf{X}_{j}^{\top},\mathbf{\xi}_{j}^{\top})^{\top})\right\|$
and $s\geq\frac{\omega_{j}\sigma^{2}}{K(\log\frac{K}{\sigma})^{1/\gamma}}$.
When

$n\gtrsim\left(K(\log K/\sigma)^{1/\gamma}/\sigma\right)^{2}\log(Thm),$

choose

$s\asymp n^{-1/2}(Th)^{-1}\sigma\log^{1/2}(Thm),\quad t\asymp
n^{-1/2}(Th)^{-1/2}\sigma\Phi_{\mathcal{Z}}\log^{1/2}m.$

If

$Th\gg\Phi_{\mathcal{Z}}^{2}\log(Thm)\log m,$

we have that $t\gg s\Phi_{\mathcal{Z}}^{2}\log m.$
Meanwhile, when

$n\lesssim\left(K(\log K/\sigma)^{1/\gamma}/\sigma\right)^{2}\log(Thm),$

choose

$s\asymp n^{-1}(Th)^{-1}K(\log K/\sigma)^{1/\gamma}\log(Thm),\quad t\asymp
n^{-1/2}(Th)^{-1/2}\sigma\Phi_{\mathcal{Z}}\log^{1/2}m.$

If

$nTh\gg(K(\log
K/\sigma)^{1/\gamma}/\sigma)^{2}\Phi_{\mathcal{Z}}^{2}\log^{2}(Thm)\log m,$

we have that $t\gg s\Phi_{\mathcal{Z}}^{2}\log m$.
Thus when

$Th\gg\Phi_{\mathcal{Z}}^{2}\log(Thm)\log m$

and

$nTh\gg(K(\log
K/\sigma)^{1/\gamma}/\sigma)^{2}\Phi_{\mathcal{Z}}^{2}\log^{2}(Thm)\log m,$

the inequality

$P\left(\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{\omega_{j}}{n_{j}}\xi_{ji}(X_{ji}-\mathbb{E}(X_{ji}))\right\|\geq
t\right)\leq\frac{1}{m}.$

holds where $t\asymp n^{-1/2}(Th)^{-1/2}\sigma\Phi_{\mathcal{Z}}\log^{1/2}m.$
Similarly, we choose $t\asymp
n^{-1/2}(Th)^{-1/2}\sigma_{\xi}\mu_{X}\Phi_{\mathcal{Y}}\log^{1/2}m$ and

$P\left(\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{\omega_{j}}{n_{j}}\xi_{ji}\mathbb{E}(X_{ji})\right\|\geq
t\right)\leq\frac{1}{m}$

holds when $Th\gg\Phi_{\mathcal{Y}}^{2}\log(Thm)\log m$ and

$nTh\gg(K_{1}(\log
K_{1}/\sigma_{\xi})^{1/\alpha}/\sigma_{\xi})^{2}\Phi_{\mathcal{Y}}^{2}\log^{2}(Thm)\log
m.$

And choose $t\asymp
n^{-1/2}(Th)^{-1/2}\varsigma\Phi_{\mathcal{X}}\log^{1/2}m$, we have

$P\left(\left\|\sum_{j=1}^{T}\sum_{i=1}^{n_{j}}\frac{\omega_{j}}{n_{j}}\left\\{\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}-E\left(\left\langle
M_{0}(j),X_{ji}\right\rangle X_{ji}\right)\right\\}\right\|\geq
t\right)\leq\frac{1}{m}$

holds when $Th\gg\Phi_{\mathcal{X}}^{2}\log(Thm)\log m$ and
$nTh\gg(\mathcal{K}(\log\mathcal{K}/\varsigma)^{2/\beta}/\varsigma)^{2}\Phi_{\mathcal{X}}^{2}\log^{2}(Thm)\log
m$.

Summarizing the above results with
$\Phi_{\mathcal{Z}}\leq\Phi_{\mathcal{X}}+\Phi_{\mathcal{Y}}$, we know that
when $Th\gg(\Phi^{2}_{\mathcal{X}}\vee\Phi_{\mathcal{Y}}^{2})\log(Thm)$ and
$nTh\gg\left(K^{*}/\sigma^{*}\right)^{2}(\Phi_{\mathcal{X}}^{2}\vee\Phi_{\mathcal{Y}}^{2})\log^{2}(Thm)\log
m$, there exists a constant $\mathcal{C}_{1}>0$, with probability at least
$1-3/m$,

$\left\|\sum_{j=1}^{T}\omega_{j}\Delta_{j}\right\|\leq\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\sigma^{*}\sqrt{\frac{\log(m_{1}+m_{2})}{nTh}}.$

Choose

$\lambda=2\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\sigma^{*}\sqrt{\frac{\log(m_{1}+m_{2})}{nTh}}$

and using Theorem 1 and proof of Theorem 2, we have that

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+(1+\sqrt{2})\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\sigma^{*}\left(\frac{r_{t}\log(m_{1}+m_{2})}{\mu^{2}m_{1}m_{2}nTh}\right)^{1/2}.$

Note that $Th\gg\log(Th)$ when $Th\rightarrow\infty$, we can represent the
condition $Th\gg(\Phi^{2}_{\mathcal{X}}\vee\Phi_{\mathcal{Y}}^{2})\log(Thm)$
and
$nTh\gg\left(K^{*}/\sigma^{*}\right)^{2}(\Phi_{\mathcal{X}}^{2}\vee\Phi_{\mathcal{Y}}^{2})\log^{2}(Thm)\log
m$ as $Th\gg(\Phi^{2}_{\mathcal{X}}\vee\Phi_{\mathcal{Y}}^{2})\log m$ and
$nTh\gg\left(K^{*}/\sigma^{*}\right)^{2}(\Phi_{\mathcal{X}}^{2}\vee\Phi_{\mathcal{Y}}^{2})\log^{3}m$
respectively.

$\square$

###### Corollary S.2

Under assumptions of Theorem 3, when $X_{ji}$ are i.i.d. uniformly distributed
on $\mathcal{E}$, $\xi_{ji}$ are independently follow sub-exponential mean-
zero distributions, $Th\gg\log(m_{1}+m_{2})$ and $nTh\gg(m_{1}\wedge
m_{2})\log^{3+2/\alpha}(m_{1}+m_{2})$, then with probability at least
$1-3/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2}$
$\displaystyle)^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}$
$\displaystyle\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+(1+\sqrt{2})\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})(C_{M}\vee\sigma_{\xi})\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nTh}\right)^{1/2}.$

When $nT\gg(m_{1}\vee m_{2})\log^{(1+5/2\alpha)}(m_{1}+m_{2})$ and
$n^{-1}T^{4}\gg(m_{1}\vee m_{2})^{-1}\log^{4}(m_{1}+m_{2})$, we select

$h=\mathcal{C}_{h}\left(\frac{(C_{M}\vee\sigma_{\xi})^{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nT}\right)^{1/5},$

then

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\mathcal{C}_{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{4/5}(C_{M}\vee\sigma_{\xi})^{4/5}\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nT}\right)^{2/5},$ (S.26)

where $\mathcal{C}_{1},\mathcal{C}_{h},\mathcal{C}_{2}$ are the same constants
as above.

The proofs of Corollary 3 and S.2 are similar to those of Corollary 1 and 2,
which we do not repeat.

##  S.6 Application to Compressed Sensing

###### Corollary S.3

(Independent case) Under Assumption 1-5, when $X_{ji}$ are random matrices
with independent mean-zero sub-gaussian elements with variance
$\sigma_{X}^{2}$, $\xi_{ji}$ are independently follow sub-gaussian
distributions, $nTh\gg m_{1}\vee m_{2}$, $h\rightarrow 0$ and
$(m_{1}m_{2})^{1/2}Th^{2}\rightarrow\infty$ as
$n,m_{1},m_{2},T\rightarrow\infty$, with probability at least
$1-2/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}$
$\displaystyle\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+\left(1+\sqrt{2}\right)C_{1}\left[\eta\vee\\{M(m_{1}m_{2})^{1/2}\\}\right]\left(\frac{r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nTh}\right)^{1/2}+o(h^{2}),$

When $nT\gg(m_{1}\vee m_{2})\log(m_{1}+m_{2})$, let

$h=C_{h}\left(\frac{\left\\{\eta^{2}\vee(M^{2}m_{1}m_{2})\right\\}r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nT}\right)^{1/5},$

then

$(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{2}\left(\frac{\left\\{\eta^{2}\vee(M^{2}m_{1}m_{2})\right\\}r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nT}\right)^{2/5},$ (S.27)

where $\eta=\sigma_{\xi}/\sigma_{X}$ and $C_{1},C_{h},C_{2}$ are the same
constants as above.

Theorem 6 in Koltchinskii (2011b) presented the error bound for static matrix
compressed sensing which is

$(m_{1}m_{2})^{-1/2}\left\|\widehat{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{4}\eta\left(\frac{r_{t}}{(m_{1}\wedge m_{2})n}\right)^{1/2}.$

Similarly to the matrix completion, when
$T\gg\eta^{-1}\\{\eta^{2}\vee(M^{2}m_{1}m_{2})\\}\\{n(m_{1}\wedge
m_{2})\\}^{1/4}\log(m_{1}+m_{2})$, our dynamic method gives a sharper bound
than the static method.

From Corollary 2 and Corollary S.3, when the variances of observation errors
$\xi_{ji}$ are small enough comparing to the variances of $X_{ji}$, i.e.,
$\sigma_{\xi}\leq M$ in matrix completion setting and $\eta\leq
M(m_{1}m_{2})^{1/2}$ in compressed sensing setting, the upper bound (15) and
(S.27) have the same order such that

$(m_{1}m_{2})^{1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq
C_{2}M^{4/5}\left(\frac{r_{t}(m_{1}\vee
m_{2})\log(m_{1}+m_{2})}{nT}\right)^{2/5}.$

Proof of Corollary S.3: Because $\mathbb{E}X_{ji}=0$, we know that $\mu_{X}=0$
and (b) in (S.17) becomes zero. From the distribution of $X_{ji}$ and
Assumption 4, we know that $\beta=2$ and $K_{2}\asymp(m_{1}\vee
m_{2})^{1/2}\sigma_{X}$ by Theorem 4.4.5 in Vershynin (2018). Also, it is easy
to check that $K\asymp\sigma_{\xi}\sigma_{X}(m_{1}\vee m_{2})^{1/2},\
\sigma=\sigma_{\xi}\sigma_{X}(m_{1}\vee m_{2})^{1/2}$. With that

$\displaystyle\left\|\mathbb{E}(\langle M_{j}^{0},X_{ji}\rangle
X_{ji}-\mathbb{E}(\langle M_{j}^{0},X_{ji}\rangle X_{ji})(\langle
M_{j}^{0},X_{ji}\rangle X_{ji}-\mathbb{E}(\langle M_{j}^{0},X_{ji}\rangle
X_{ji})^{T}\right\|_{\infty}$ $\displaystyle\leq$ $\displaystyle
M^{2}m_{1}m_{2}(m_{1}\vee m_{2})\sigma^{4}_{X}$

and with Lemma S.3

$\displaystyle\left\|\langle M_{j}^{0},X_{ji}\rangle X_{ji}-\mathbb{E}\langle
M_{j}^{0},X_{ji}\rangle X_{ji}\right\|_{\psi(1)}$
$\displaystyle\asymp\left\|\langle M_{j}^{0},X_{ji}\rangle
X_{ji}\right\|_{\psi(1)}$ $\displaystyle\leq\left\|\langle
M_{j}^{0},X_{ji}\rangle\right\|_{\psi(2)}\left\|X_{ji}\right\|_{\psi(2)}$
$\displaystyle\leq M(m_{1}m_{2})^{1/2}\sigma_{X}^{2},$

we have that $\varsigma\leq M(m_{1}m_{2})^{1/2}(m_{1}\vee
m_{2})^{1/2}\sigma^{2}_{X},\mathcal{K}\lesssim
M(m_{1}m_{2})^{1/2}\sigma_{X}^{2}$.

Then it can be obtained that

$\displaystyle\sigma^{*}=\left[\sigma_{\xi}\vee
M\sigma_{X}(m_{1}m_{2})^{1/2}\right](m_{1}\vee m_{2})^{1/2}\sigma_{X}$
$\displaystyle K^{*}/\sigma^{*}\asymp(m_{1}\vee m_{2})^{-1/2}\log(1/(m_{1}\vee
m_{2})).$

Apply Theorem 2 and Corollary 1, the proof is completed. $\square$

###### Corollary S.4

(Dependent case.) Under assumptions of Theorem 3, when $X_{ji}$ are random
matrices with independent mean-zero sub-gaussian elements with variance
$\sigma_{X}^{2}$ and $\xi_{ji}$ are independently follow sub-gaussian
distributions, $Th\gg\log(m_{1}+m_{2})$ and $nTh\gg(m_{1}\vee m_{2})$, with
probability at least $1-2/(m_{1}+m_{2})$,

$\displaystyle(m_{1}m_{2})^{-1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}$
$\displaystyle\leq\frac{1}{2}\alpha(K)D_{2}h^{2}+\left(1+\sqrt{2}\right)\mathcal{C}_{1}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})\left[\eta\vee
M(m_{1}m_{2})^{1/2}\right]\left(\frac{r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nTh}\right)^{1/2}.$

When $nT\gg(m_{1}\vee m_{2})\log(m_{1}+m_{2})$ and $n^{-1}T^{4}\gg(m_{1}\vee
m_{2})^{-1}\log^{4}(m_{1}+m_{2})$, we select

$h=\mathcal{C}_{h}\left(\frac{\left(\eta^{2}\vee
M^{2}m_{1}m_{2}\right)(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{2}r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nT}\right)^{1/5},$

then

$(m_{1}m_{2})^{1/2}\left\|\widetilde{M}_{t}^{\lambda}-M_{t}^{0}\right\|_{2}\leq\mathcal{C}_{2}(\Phi_{\mathcal{X}}\vee\Phi_{\mathcal{Y}})^{4/5}\left(\frac{\left(\eta^{2}\vee
M^{2}m_{1}m_{2}\right)r_{t}\log(m_{1}+m_{2})}{(m_{1}\wedge
m_{2})nT}\right)^{2/5},$

where $\eta=\sigma_{\xi}/\sigma_{X}$ and
$\mathcal{C}_{1},\mathcal{C}_{h},\mathcal{C}_{2}$ are the same constants as
above.

The proof of Corollary S.4 is similar to that of Corollary S.3, which we do
not repeat.

##  S.7 Proof of Theorem 4

Proof of Theorem 4: First we know that

$\displaystyle\nabla
f_{t}(M)=2\sum_{j=1}^{T}\omega_{j}\left[\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}\langle
M,X_{ji}\rangle X_{ji}-\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}Y_{ji}X_{ji}\right].$

This gives

$\displaystyle\left\|\nabla f_{t}(M)-\nabla f_{t}(N)\right\|_{2}$
$\displaystyle=\left\|2\sum_{j=1}^{T}\frac{\omega_{j}}{n_{j}}\sum_{i=1}^{n_{j}}\mathbb{E}\langle
M-N,X_{ji}\rangle X_{ji}\right\|_{2}$ $\displaystyle\leq
2\left\|\sum_{j=1}^{T}\frac{\omega_{j}}{n_{j}}\sum_{i=1}^{n_{t}}\left(\|X_{ji}\|_{2}X_{ji}\right)\right\|_{2}\|M-N\|_{2}.$

Define that
$2L_{f}=4\left\|\sum_{j=1}^{T}\omega_{h}(j-t)/n_{j}\sum_{i=1}^{n_{j}}\left(\|X_{ji}\|_{2}X_{ji}\right)\right\|_{2}$,
Toh & Yun (2010) proved that for any $k>1$,

$F_{t}(M_{t}^{(k)})-F_{t}(\widetilde{M}_{t}^{\lambda})\leq\frac{2L_{f}\|M_{t}^{(0)}-\widetilde{M}_{t}^{\lambda}\|_{2}^{2}}{(k+1)^{2}}.$
(S.28)

Note that (S.28) focuses on the convergence rate of object function