This paper studies the multi-task high-dimensional linear regression models where the noise among different tasks is correlated, in the moderately high dimensional regime where sample size $n$ and dimension $p$ are of the same order.
Our goal is to estimate the covariance matrix of the noise random vectors, or equivalently the correlation of the noise variables on any pair of two tasks. Treating the regression coefficients as a nuisance parameter, we leverage the multi-task elastic-net and multi-task lasso estimators to estimate the nuisance. By precisely understanding the bias of the squared residual matrix and by correcting this bias, we develop a novel estimator of the noise covariance that converges in Frobenius norm at the rate $n^{-1/2}$ when the covariates are Gaussian. This novel estimator is efficiently computable.

Under suitable conditions, the proposed estimator of the noise covariance attains the same rate of convergence as the “oracle” estimator that knows in advance the regression coefficients of the multi-task model. The Frobenius error bounds obtained in this paper also illustrate the advantage of this new estimator compared to a method-of-moments estimator that does not attempt to estimate the nuisance.

As a byproduct of our techniques, we obtain an estimate of the generalization error of the multi-task elastic-net and multi-task lasso estimators. Extensive simulation studies are carried out to illustrate the numerical performance of the proposed method.

§ INTRODUCTION

§.§ Model and estimation target

Consider a multi-task linear model with $T$ tasks and $n$ observations $(\bx_i, Y_{i1}, Y_{i2},\dots, Y_{iT})$, $\forall i=1,...,n$, where $\bx_i\in \R^p$ is a random feature vector and $Y_{i1}, \ldots, Y_{iT}$ are responses
in the model

\begin{equation}\label{eq: model}
\begin{aligned}
Y_{it} &= \bx_i^\top \bbeta^{(t)} + E_{it} &&\text{for each } t = 1, ..., T; i=1,...,n
&&\text{(scalar form)},
\\
\by\smash{{}^{(t)}}                             &= \bX \bbeta\smash{{}^{(t)}} + \bep\smash{{}^{(t)}}
&&\text{for each } t = 1, ..., T
&&\text{(vector form)},
\\
\bY &= \bX\bB^* + \bE
&&\text{(matrix form)},
\end{aligned}
\end{equation}

where $\bX\in\R^{n\times p}$ is the design matrix with rows $(\bx_i^{\top})_{i=1,...,n}$,
$\by^{(t)} = (Y_{1t},..., Y_{nt})^\top$ is the response vector for task $t$,
$\bep^{(t)} = (E_{1t},...,E_{nt})^\top$ is the noise vector for task $t$,
$\bbeta^{(t)} \in \R^{p}$ is an unknown fixed coefficient vector for task $t$.
In matrix form,
$\bY\in \R^{n\times T}$ is the response matrix with columns $\by^{{(1)}},...,\by^{(T)}$,
$\bE\in \R^{n\times T}$ has columns $\bep^{{(1)}},...,\bep^{(T)}$,
and $\bB^*\in\R^{p\times T}$ is an unknown coefficient matrix with
columns $\bbeta^{{(1)}},...,\bbeta^{(T)}$. The three forms in (<ref>) are equivalent.

While the $n$ vectors $(\bx_i^\top, y_i^{(1)}, \ldots, y_i^{(T)})_{i=1,...,n}$ of dimension $p+T$ are ,
we assume that for each observation $i=1,...,n$, the noise random variables $E_{i1},...,E_{iT}$ are centered and correlated.

The focus of the present paper is on estimation of the
noise covariance matrix $\bS\in\R^{T\times T}$, which has entries
$\bS_{tt'} = \E[\varepsilon_1^{(t)}\varepsilon_1^{(t')}]$ for any pair $t,t'=1,\ldots,T$, or equivalently

\bS = \E[\tfrac1n\bE^\top\bE].

The noise covariance plays a crucial role in multi-task linear models because it characterizes the noise level and correlation between different tasks:
if tasks $t=1,...,T$ represent time this captures temporal correlation;
if tasks $t=1,...,T$ represent different activation areas in the brain (, [Bertrand et al., 2019]) this captures spatial correlation.

Since $\bS$ is the estimation target,
we view $\bB^*$ as an unknown nuisance parameter.
If $\bB^*=\mathbf0$,  then $\bY = \bE$, hence $\bE$ is directly observed and a natural estimator is the sample covariance
There are other possible choices for the sample covariance;
ours coincides with the maximum likelihood estimator of the centered Gaussian model
where the $n$ samples are from $\mathcal N_T(\bf 0, \bS)$.

In the presence of a nuisance parameter $\bB^*\neq \bf 0$,
the above sample covariance is not computable since we only observe $(\bX, \bY)$ and do not have access to $\bE$.
Thus we will refer to
$\frac1n\bE^\top\bE \in\R^{T\times T}$
as the oracle estimator for $\bS$, and its error
$\frac1n\bE^\top\bE - \bS$
will serve as a benchmark.

The nuisance parameter $\bB^*$ is not of interest by itself, but if an estimator $\hbB$ is available that provides good estimation
of $\bB^*$, we would hope to leverage $\hbB$ to estimate the nuisance
and improve estimation of $\bS$.
For instance given an estimate $\hbB$ such that
$\fnorm{\bX(\hbB-\bB^*)}^2/n\to0$, one may use the estimator

\begin{equation}
\textstyle
\label{naive}
\hbS_{(\text{naive})} = \frac1n (\bY - \bX\hbB)^\top(\bY - \bX\hbB)
\end{equation}

to consistently estimate $\bS$ in Frobenius norm.
We refer to this estimator as the naive estimator since it is obtained by simply replacing the noise $\bE$ in the oracle estimator $\frac1n \bE^\top\bE$ with the residual matrix $\bY - \bX\hbB$.
in the regime $p/n\to\gamma$ of interest in the present paper,
the convergence $\fnorm{\bX(\hbB-\bB^*)}^2/n\to0$ is not true even for $T=1$
and common high-dimensional estimators such as Ridge regression [Dobriban and Wager, 2018] or the Lasso [Bayati and Montanari, 2012, Miolane and Montanari, 2018].
Simulations in <Ref> will show that
(<ref>) presents a major bias for estimation of $\bS$.
One goal of this paper is to develop estimator $\hbS$ of $\bS$
by exploiting a commonly used estimator $\hbB$ of the nuisance,

so that in the regime $p/n\to\gamma$ the error $\hbS-\bS$ is comparable to the benchmark $\frac1n\bE^\top\bE - \bS$.

§.§ Related literature

If $T=1$, the above model (<ref>) reduces to the standard linear model
with $\bX\in\R^{n\times p}$ and response vector $\by^{(1)}\in \R^n$.
We will refer to the $T=1$ case as the single-task linear model and drop the superscript $^{(1)}$ for brevity, ,

$y_i = \bx_i^\top \bbeta^* + \ep_i$,

where $\ep_i$ are with mean $0$, and unknown variance $\sigma^2$. The coefficient vector $\bbeta^*$ is typically assumed to be $s$-sparse, , $\bbeta^*$ has at most $s$ nonzero entries. In this single-task linear model, estimation of noise covariance $\bS$ reduces to estimation of the noise variance  $\sigma^2=\E[\ep_i^2]$, which has been studied in the literature. Fan et al., 2012 proposed a consistent estimator for $\sigma^2$ based on a refitted cross validation method, which assumes the support of $\bbeta^*$ is correctly recovered; [Belloni et al., 2011] and  [Sun and Zhang, 2012] introduced square-root Lasso (scaled Lasso) to jointly estimate the coefficient $\bbeta^*$ and noise variance $\sigma^2$ by

\begin{equation}\label{eq: scaled-lasso}
\textstyle
(\hbbeta, \hsigma) =\argmin_{\bbeta\in \R^p, \sigma>0} \frac{\norm{\by - \bX\bbeta}^2}{2n\sigma} + \frac{\sigma}{2} + \lambda_0\norm{\bbeta}_1.
\end{equation}

This estimator $\hsigma$ is consistent only when the prediction error $\norm{\bX(\hbbeta - \bbeta^*)}^2 /n$ goes to 0, which requires $s\log(p)/n\to 0$.
Estimation of $\sigma^2$ without assumption on $\bX$ was proposed in [Yu and Bien, 2019] by utilizing natural parameterization of the penalized likelihood of the linear model. Their estimator can be expressed as the minimizer of the Lasso problem:
$\hsigma^2_{\lambda} = \min_{\bbeta\in \R^p} \frac{1}{n} \norm{\by - \bX\bbeta}^2 + 2\lambda\norm{\bbeta}_1.$ Consistency of these estimators [Sun and Zhang, 2012, Belloni et al., 2011, Belloni et al., 2014, Yu and Bien, 2019] requires $s\log(p)/n \to 0$ and does not hold in the high-dimensional proportional regime $p/n\to\gamma\in (0, \infty)$.
For this proportional regime $p/n \to \gamma \in (0, \infty)$, [Dicker, 2014] introduced a method-of-moments estimator $\hsigma^2$ of $\sigma^2$,

\begin{equation}\label{eq: est-dicker14}
\hsigma^2 = \frac{n+p+1}{n(n+1)} \norm{\by}^2 - \frac{1}{n(n+1)}\norm{\bSigma^{-\frac12}\bX^\top\by}^2,
\end{equation}

which is unbiased, consistent, and asymptotically normal in high-dimensional linear models with Gaussian predictors and errors. Moreover, [Janson et al., 2017] developed an EigenPrism procedure for the same task as well as confidence intervals for $\sigma^2$.
The estimation procedures in these two papers don't attempt to estimate the nuisance parameter $\bbeta^*$, and require no sparsity on $\bbeta^*$ and isometry structure on $\bSigma$, but assume $\|\bSigma^{\frac 12}\bbeta^*\|^2$ is bounded.
Maximum Likelihood Estimators (MLEs) were studied in
[Dicker and Erdogdu, 2016] for joint estimation of noise level and signal strength in high-dimensional linear models with fixed effects; they showed that a classical MLE for random-effects models may also be used effectively in fixed-effects models.

In the proportional regime,
[Bayati et al., 2013, Miolane and Montanari, 2018] used the Lasso to estimate the nuisance $\bbeta^*$ and produce estimator for $\sigma^2$. Their approach requires an uncorrelated Gaussian design assumption with $\bSigma = \bI_p$.
Bellec, 2020 provided consistent estimators of a similar nature for $\sigma^2$ using more general M-estimators with convex penalty without requiring $\bSigma = \bI_p$. In the special case of the squared loss, this estimator has the form [Bayati et al., 2013, Miolane and Montanari, 2018, Bellec, 2020]

\begin{equation}\label{eq: est-bellec20}
\hsigma^2 = (n -\df)^{-2}\big\{ \norm{\by-\bX\hbbeta}^2 (n+p -2\df) - \norm{\bSigma^{-\frac12}(\by-\bX\hbbeta)}\big\},
\end{equation}

where $\df = \trace[(\partial/\partial \by) \bX\hbbeta]$ denotes the degrees of freedom. This estimator coincides with the method-of-moments estimator in [Dicker, 2014] when $\hbbeta = \bf0$.

For multi-task high-dimensional linear model (<ref>) with $T\ge 2$, the estimation of $\bB^*$ is studied in [Lounici et al., 2011], [Obozinski et al., 2011], [Simon et al., 2013]. These works suggest to use a joint convex optimization problem over the tasks to estimate $\bB^*$. A popular choice is the multi-task elastic-net, which
solves the convex optimization problem

\begin{equation}\label{eq: hbB}
\hbB=\argmin_{\bB\in\R^{p\times T}}
\Big(
\frac{1}{2n}\fnorm*{\bY - \bX\bB }^2 + \lambda \norm{\bB}_{2,1}
+ \frac{\tau}{2} \fnorm{\bB}^2
\Big),
\end{equation}

where $\|\bB\|_{2,1} = \sum_{j=1}^p \|{\bB^{\top} \be_j}\|_2$, and $\fnorm{\cdot}$ denotes the Frobenius norm of a matrix.
This optimization problem can be efficiently solved by existing statistical packages, for instance, scikit-learn [Pedregosa et al., 2011], and glmnet  [Friedman et al., 2010].
Note that  (<ref>) is also referred to as multi-task (group) Lasso and multi-task Ridge if $\tau = 0$ and $\lambda =0$, respectively.
van de Geer and Stucky, 2016 extended square-root Lasso [Belloni et al., 2011] and scaled Lasso [Sun and Zhang, 2012] to multi-task setting by solving the following problem

\begin{align}\label{eq: multi-ScaledLasso}
(\hbB, \hbS) =
\argmin_{\bB,\bS \succ 0} \Big\{\frac{1}{n} \trace\big((\bY - \bX\bB)\bS^{-\frac 12}(\bY - \bX\bB)^\top\big) + \trace(\bS^{\frac 12}) + 2\lambda_0\norm{\bB}_1 \Big\},
\end{align}

where $\norm{\bB}_1 = \sum_{j,t} |B_{jt}|$.
Note that the covariance estimator in (<ref>) is constrained to be positive definite.
Molstad, 2019 studied the same problem and proposed to estimate $\bS$ by (<ref>) with $\hbB$ in (<ref>), which is consistent under Frobenius norm loss when $\fnorm{\bX(\hbB-\bB^*)}^2/n\to0$.
In a recent paper, Bellec and Romon, 2021 studied the multi-task Lasso problem and proposed confidence intervals for single entries of $\bB^*$ and confidence ellipsoids for single rows of $\bB^*$ under the assumption that $\bS$ is proportional to the identity, which may be restrictive in practice. This literature generalizes degrees of freedom adjustments from single-task to multi-task models, which we will illustrate in <Ref>.

Noise covariance estimation in the high dimensional multi-task linear model is a difficult problem. If the estimand $\bS$ is known to be diagonal,
estimating $\bS$ reduces to the estimation of noise variance for each task,
in which the existing methods for single-task high-dimensional linear models can be applied.
Nonetheless, for general positive semi-definite matrix $\bS$, the noise among different tasks may be correlated,
hence the existing methods are not readily applicable, and a more careful analysis is called for to incorporate the correlation between different tasks. Fourdrinier et al., 2021 considered estimating $\bS$ for the multi-task model (<ref>) where rows of $\bE$ have elliptically symmetric distribution and in the classical regime $p\le n$. However, their estimator has no statistical guarantee under Frobenius norm loss.

Recently, for the proportional regime $p/n \to \gamma \in (0, \infty)$, [Celentano and Montanari, 2021] generalized the estimator $\hsigma^2$ in [Bayati et al., 2013] to the multi-task setting with $T=2$. Their work covers correlated Gaussian designs, where a Lasso or Ridge regression is used to estimate $\bbeta^{(1)}$ for the first task, and another Lasso or Ridge regression is used to estimate $\bbeta^{(2)}$ for the second task.
In other words, they estimate the coefficient vector for each task separately instead of using a multi-task estimator like (<ref>).
It is not trivial to adapt their estimator from the setting $T=2$ to larger $T$, and allow $T$ to increase with $n$.
This present paper takes a different route and aims to fill this gap by proposing a novel noise covariance estimator with theoretical guarantees.
Of course, our method applies directly to the 2-task linear model considered in [Celentano and Montanari, 2021].

§.§ Main Contributions

The present paper introduces a novel estimator $\hbS$ in (<ref>) of the noise covariance $\bS$, which provides consistent estimation of $\bS$ in Frobenius norm,
in the regime where $p$ and $n$ are of the same order.
The estimator $\hbS$ is based on the multi-task elastic-net
estimator $\hbB$ in (<ref>) of the nuisance, and can be seen
as a de-biased version of the naive estimator (<ref>).
The naive estimator (<ref>) suffers from a strong bias in the regime
where $p$ and $n$ are of the same order,
and the estimator $\hbS$
is constructed by precisely understanding this bias and correcting it.

After introducing this novel estimator $\hbS$ in <Ref> below,
we prove several rates of convergence for the Frobenius error
$\fnorm{\hbS-\bS}$, which is comparable, in terms of rate of convergence,
to the benchmark $\fnorm{\frac 1 n \bE^\top\bE - \bS}$ under suitable assumptions.

As a by-product of the techniques developed for the construction of $\hbS$,
we obtain estimates of the generalization error

of $\hbB$, which are of independent interest and can be used for parameter tuning.

§.§ Notation

Basic notation and definitions that will be used in the rest of the paper are given here. Let $[n] = \{1, 2,\ldots, n \}$ for all $n\in\N$.
The vectors $\be_i\in \R^n,\be_j\in\R^p, \be_t\in\R^T$ denote the canonical basis vector of the corresponding index.
We consider restrictions of vectors (, of matrices)
by zeroing the corresponding entries (, columns).
More precisely, for $\bv\in\R^p$ and index set $B\subset [p]$, $\bv_B\in\R^p$ is the vector
with $(\bv_B)_j = 0$ if $j\notin B$ and $(\bv_B)_j = v_j$ if $j\in B$.
If $\bX\in\R^{n\times p}$ and $B\subset[p]$, $\bX_B\in\R^{n\times p}$
such that $(\bX_B)\be_j = {\mathbf 0}$ if $j\notin B$ and
$(\bX_B)\be_j= \bX\be_j$ if $j\in B$.

For a real vector $\ba \in \R^p$,
denotes its Euclidean norm.
For any matrix $\bA$, $\bA^\dagger$ is its Moore–Penrose inverse;
denote its Frobenius, operator and nuclear norm, respectively.
Let $\norm{\bA}_0$ be the number of non-zero rows of $\bA$.
Let $\bA\otimes \bB$ be the Kronecker product of $\bA$ and $\bB$, and $\langle \bA, \bB\rangle = \trace(\bA^\top\bB)$
is the Frobenius inner product for matrices of identical size.
For $\bA$ symmetric, $\phi_{\min}(\bA)$ and $\phi_{\max}(\bA)$ denote its smallest and largest eigenvalues, respectively.
Let $\bI_n$ denote the identity matrix of size $n$ for all $n\in\N$.
For a random sequence $\xi_n$, we write $\xi_n = O_P(a_n)$ if $\xi_n/a_n$ is stochastically bounded.
$C$ denotes an absolute constant and
$C(\tau, \gamma)$ stands for a generic positive constant depending on $\tau,\gamma$;
their expression may vary from place to place.

§.§ Organization

The rest of the paper is organized as follows.
<Ref> introduces our proposed estimator for noise covariance.
<Ref> presents our main theoretical results on proposed estimator and some relevant estimators.
<Ref> demonstrates through numerical experiments that our estimator outperforms several existing methods in the literature, which corroborates our theoretical findings in <Ref>.
<Ref> provides discussion and points out some future research directions.
Proofs of all the results stated in the main body are given in the supplementary, which starts with an outline for ease of navigation.

§ ESTIMATING NOISE COVARIANCE, WITH POSSIBLY DIVERGING NUMBER OF TASKS  T

Before we can define our noise covariance estimator, we need to introduce the following building blocks.
Let $\hat{\mathscr{S}} = \{k\in [p]: {\hbB{}^{\top} \be_k} \neq 0\}$

denote the set of nonzero rows of $\hbB$ in (<ref>), and  let $|\hat{\mathscr{S}}|$ denote the cardinality of $\hat{\mathscr{S}}$.
For each $k\in \hat{\mathscr{S}}$, define $\bH^{(k)}=\lambda\|\hbB{}^\top \be_k\|^{-1}(\bI_T - \hbB{}^\top\be_k \be_k^\top\hbB ~  \|\hbB{}^\top\be_k\|^{-2} )$, which is the Hessian of the map $\bu \mapsto \lambda\norm*{\bu}$ at $\bu = \hbB{}^\top \be_k$ when $\bu\ne \mathbf0$.

Define $\bM,\bM_1\in\R^{pT \times pT}$ by

\begin{equation}
\textstyle
\bM_1 = \bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})
\qquad
\bM = \bM_1 + n\sum_{k\in\hat{\mathscr{S}}} (\bH^{(k)} \otimes \be_k\be_k^\top)
\end{equation}

where $\bP_{\hat{\mathscr{S}}} = \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top\in\R^{p\times p}$.

Define the residual matrix $\bF$, the error matrix $\bH$, and $\bN$ by

\begin{equation}
\bF=\bY-\bX\hbB,
\qquad
\bH = \bSigma^{1/2}(\hbB - \bB),
\qquad
\bN = (\bI_T \otimes \bX)\bM^\dagger (\bI_T \otimes \bX^\top)
\in \R^{Tn \times Tn}.
\label{eq:def F H N}
\end{equation}

To construct our estimator we also make use of the so-called interaction matrix $\hbA\in \R^{T\times T}$.

The interaction matrix $\hbA\in \R^{T\times T}$ of the estimator $\hbB$ in  (<ref>) is defined by

\begin{align}\label{eq: hbA-matrix}
\hbA
= \sum_{i=1}^n (\bI_T \otimes \be_i^\top \bX) \bM^{\dagger} (\bI_T \otimes \bX^\top\be_i)
= \sum_{i=1}^n (\bI_T \otimes \be_i^\top)\bN(\bI_T \otimes \be_i).
\end{align}

The matrix $\hbA$ was introduced in [Bellec and Romon, 2021], where it is used alongside the multi-task Lasso estimator ($\tau=0$ in (<ref>)).
It generalizes the degrees of freedom from Stein, 1981 to the multi-task case.
Intuitively, it captures the correlation between the residuals on different tasks <cit.>.
Our definition of the noise covariance estimator involves $\hbA$,
although our statistical purposes differ greatly from the confidence intervals developed in [Bellec and Romon, 2021].

We are now ready to introduce our estimator $\hbS$ of the noise covariance $\bS$.

With $\bF=\bY-\bX\hbB$ and $\hbA$ as above, define

\begin{equation}\label{eq: hbS}
\hbS = (n\bI_T - \hbA)^{-1} \Bigl[\bF^\top \big( (p+n)\bI_n - \bX \bSigma^{-1}\bX^\top\big) \bF - \hbA\bF^\top\bF
- \bF^\top\bF\hbA
\Bigr](n\bI_T - \hbA)^{-1}.
\end{equation}

Efficient solvers (,  in <cit.>) are available to compute $\hbB$.
Computation of $\bF$ is then straightforward, and computing the matrix $\hbA$
only requires inverting a matrix of size $|\hat{\mathscr{S}}|$ <cit.>.
The estimator $\hbS$ generalizes the scalar estimator
(<ref>) to the multi-task setting
in the sense that for $T=1$, $\hbS$ is exactly equal to (<ref>).

Note that unlike in (<ref>), here $\bF^\top\bF$, $\hbA$ and $(n\bI_T - \hbA)$ are matrices of size $T\times T$: the order of matrix multiplication in $\hbS$ matters and should not be switched.
This non-commutativity is not present for $T=1$ in (<ref>) where matrices in $\R^{T\times T}$ are reduced to scalars.
Another special case of $\hbS$ can be seen in [Celentano and Montanari, 2021] for $T=2$ where the matrix $\hbA\in\R^{2\times 2}$ is diagonal and the two columns of $\hbB\in\R^{p\times 2}$ are two Lasso or Ridge estimators computed independently of each other, one for each task. Except in these two special cases — (<ref>) for $T=1$, [Celentano and Montanari, 2021] for $T=2$ and two Lasso/Ridge — we are not aware of previously proposed estimators of the same form as $\hbS$.

§ THEORETICAL ANALYSIS

§.§ Oracle and method-of-moments estimator

Before moving on to the theoretical analysis of $\hbS$, we state our randomness assumptions for $\bE, \bX$ and we study two preliminary estimators:
the oracle $\frac 1 n \bE^\top\bE$ and another estimator obtained by the method of moments.

[Gaussian noise]
$\bE\in \R^{n\times T}$ is a Gaussian noise matrix with $\mathcal N_T(\bf0,\bS)$ rows, where $\bS\in \R^{T\times T}$ is an unknown positive semi-definite matrix.

An oracle with access to the noise matrix $\bE$ may compute the oracle estimator $\hbS_{\rm{(oracle)}} \defas \frac 1n \bE^\top \bE$, with convergence rate given by the following theorem, which will serve as a benchmark.

Under <Ref>,

\begin{equation}
\E\big[	\fnorm{\hbS_{\rm{(oracle)}}- \bS}^2 \big] = \tfrac1n [(\trace(\bS))^2 + \trace(\bS^2)].
\label{eq: bound oracle}
\end{equation}

Consequently, $n^{-1} (\trace(\bS))^2 \le
\E\big[ \fnorm{\hbS_{\rm{(oracle)}} - \bS}^2 \big]  \le 2 n^{-1} (\trace(\bS))^2$.

The next assumption concerns the design matrix $\bX$ with rows $\bx_1^\top,\ldots,\bx_n^\top$.
[Gaussian design]
$\bX\in \R^{n\times p}$ is a Gaussian design matrix with $\mathcal N_p(\mathbf 0,\bSigma)$ rows, where  $\bSigma$ is a known positive definite matrix. The matrices $\bE$ and $\bX$ are independent.

Under the preceding assumptions, we obtain the following method-of-moments estimator, which extends the estimator for noise variance in [Dicker, 2014] to the multi-task setting. Its error will also serve as a benchmark.

Under <Ref>, the method-of-moments estimator defined as

\begin{equation}
\label{hbS_mm}
\hbS_{\rm{(mm)}} = \frac{(n+1+p)}{n(n+1)} \bY^\top \bY - \frac{1}{n(n+1)}  \bY^\top\bX \bSigma^{-1}\bX^\top\bY
\end{equation}

is unbiased for $\bS$, , $\E [\hbS_{\rm{(mm)}} ] = \bS.$ Furthermore, the Frobenius error is bounded from below as

\begin{align}
\E [\fnorm{\hbS_{\rm{(mm)}} - \bS}^2 ] \ge \frac{p-2}{(n+1)^2} \big[\trace(\bS) + \fnorm{\bSigma^{\frac12}\bB^*}^2\big]^2.
\label{eq:lower-boud-mom}
\end{align}

By (<ref>), a larger norm $\fnorm{\bSigma^{1/2}\bB^*}$ induces a larger variance for $\hbS_{\rm{(mm)}}$.
Our goal with an estimate $\hbS$, when a good estimator $\hbB$ of the nuisance is available, is to improve upon the right-hand side of (<ref>)
when the estimation error $\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}$ is smaller than

A high-probability upper bound of the form
$\fnorm{\hbS_{\rm{(mm)}} - \bS}^2
\le C \frac{n+p}{n^2}[\trace(\bS) + \fnorm{\bSigma^{\frac12}\bB^*}^2]^2
that matches the lower bound (<ref>) when $p>n$,
is a consequence of our main result below.
Indeed, when $\hbB=\bf 0$ then $\hbA=\bf0$ and our estimator $\hbS$ from <Ref>
coincides with $\hbS_{\rm{(mm)}}$ up to the minor modification
of replacing $n+1$ by $n$ in (<ref>).
This replacement is immaterial
compared to the right-hand side in (<ref>).
Furthermore, such $\hbS$
corresponds to one of $\tau$ or $\lambda$ being $+\infty$ in (<ref>)
and the aforementioned upper bound

follows by taking $\tau=+\infty$ in the proof of <Ref> below.
The empirical results in <Ref> confirm that $\hbS$ has smaller variance compared to $\hbS_{\rm{(mm)}}$ in simulations.

§.§ Theoretical results for proposed estimator

We have established lower bounds for the oracle estimator and the
method-of-moments estimator that will serve as benchmarks.
We turn to the analysis of the estimator $\hbS$ from <Ref>
under the following additional assumptions.

[High-dimensional regime]
$n,p$ satisfy $p/n \le \gamma$ for a constant $\gamma\in (0, \infty)$.

For asymptotic statements such as those involving the stochastically
bounded notation $O_p(\cdot)$ or the convergence in probability in (<ref>) below, we implicitly consider a sequence
of multi-task problems indexed by $n$ where $p,T,\bB^*,\hbB,\bS$ all implicitly
depend on $n$. The Assumptions, such as $p/n\le\gamma$ above, are required to hold at
all points of the sequence. In particular, $p/n\to\gamma'$ is allowed for any limit $\gamma'\le \gamma$ under <Ref>, although our results do not require a specific value for the limit.

Assume either one of the following:

* $\tau>0$ in the penalty of estimator (<ref>), and let $\tau' = \tau/\opnorm{\bSigma}$.

* $\tau=0$ and for $c>0$, $P(U_1) \ge 1-\frac 1T$ and $P(U_1)\to1$ as $n\to\infty$, where $U_1 = \{\norm{\hbB}_0 \le n(1-c)/2 \}$ is the event that $\hbB$ has at most $n(1-c)/2$ nonzero rows. Finally, $T \le e^{\sqrt{n}}$.

<Ref>(i) requires that the Ridge penalty in (<ref>) be enforced,
so that the objective function is strongly convex.
<Ref>(ii), on the other hand, does not require strong convexity
but that the number of nonzero rows of $\hbB$ is small enough with high-probability,
which is a reasonable assumption when the tuning parameter $\lambda$ in (<ref>) is large enough and $\bB^*$ is sparse enough. While we do not prove
in the present paper that $\P(U_1)\to1$ under assumptions on the tuning parameter $\lambda$ and the sparsity of $\bB^*$, results of a similar nature have been obtained previously in several group-Lasso settings

Suppose that <Ref> hold for all $n, p$ as $n\to\infty$, then almost surely

\begin{equation}
\fnorm{(\bI_T - \hbA/n) (\hbS -\bS)(\bI_T - \hbA/n)} \le \Theta_1 n^{-\frac 12} \big( \fnorm{\bF}^2/n + \fnorm{\bH}^2 +\trace(\bS)\big)
\label{eq:thm33}
\end{equation}

for some non-negative random variable $\Theta_1$ of constant order, in the sense that $\E [\Theta_1^2] \le C(\tau')(T \wedge (1 + \frac pn))(1 + \frac pn))\le C(\gamma, \tau')$
under <Ref>(i),
and ${\E [I(\Omega)\Theta_1^2]} \le C(\gamma,c)$ under <Ref>(ii), where $I(\Omega)$ is the indicator function of an event $\Omega$
with $\P(\Omega)\to 1$.

Above, $\Theta_1\ge 0$ is said to be of constant order because $\Theta_1=O_P(1)$
follows from $\E[\Theta_1^2]\le C(\gamma,\tau')$ or from $\E[I(\Omega)\Theta_1^2]\le C(\gamma,c)$ if the stochastically bounded notation $O_P(1)$ is allowed to
hide constants depending on $(\gamma,\tau')$ or $(\gamma,c)$ only.
In the left-hand side of (<ref>), multiplication by $\bI_T - \hbA/n$ on both sides of the error $\hbS -\bS$ can be further removed, as

\begin{equation}
\fnorm{\hbS -\bS}
\le
\fnorm{(\bI_T - \hbA/n) (\hbS -\bS)(\bI_T - \hbA/n)}
\opnorm{(\bI_T - \hbA/n)^{-1}}^2
\label{eq: ineq op norm}
\end{equation}

and the fact that $\opnorm{(\bI_T - \hbA/n)^{-1}}$ is bounded from above with high probability by a constant depending on $\gamma, \tau', c$ only.
Upper bounds on $\opnorm{(\bI_T - \hbA/n)^{-1}}$ are formally stated in the supplementary material.

§.§ Understanding the right-hand side of (14), and the multi-task generalization error

Before coming back to upper bounds on the error $\fnorm{\hbS-\bS}$,
let us study the quantities appearing in the right-hand side of (<ref>). By (<ref>), $\fnorm{\bF}^2/n$
is the mean squared norm of the residuals and is observable, while the squared error
$\fnorm{\bH}^2=\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$ and $\trace[\bS]$ are unknown.
By analogy with single task models, we define the generalization error as the matrix $\bH^\top\bH + \bS$ of size $T\times T$, whose $(t,t')$-th entry is
$\E[(Y^{new}_t - \bx_{new}^T\hbB\be_t)(Y^{new}_{t'} - \bx_{new}^T\hbB\be_{t'})|(\bX,\bY)]$
where $(Y^{new}_t,Y^{new}_{t'},\bx_{new})$ is independent of $(\bX,\bY)$
and has the same distribution
as $(Y_{it},Y_{it'},\bx_i)$ for some $i=1,...,n$.
Estimating the generalization error is useful for parameter tuning:

\begin{equation}
\trace[\bH^T\bH + \bS] = \fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2 + \trace[\bS],
\label{eq:trace-gen}
\end{equation}

minimizing an estimator of $\trace[\bH^T\bH + \bS]$ is a useful proxy
to minimize the Frobenius error $\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$
of $\hbB$.
The following theorem gives an estimate for the generalization error matrix
as well as a consistent estimator for its trace (<ref>).

Let <Ref> be fulfilled. Then

\begin{align*}
\fnorm{ \bF^\top\bF/n - (\bI_T - \hbA/n) (\bH^\top\bH + \bS) (\bI_T - \hbA/n)} \le \Theta_2 n^{-\frac 12} \big(\fnorm*{\bF}^2/n + \fnorm{\bH}^2 + \trace(\bS)\big),
\end{align*}

for some non-negative random variable $\Theta_2$ of constant order, in the sense that $\E [\Theta_2] \le C(\gamma,\tau')$ under <Ref>(i),
and with ${\E [I(\Omega)\Theta_2]} \le C(\gamma,c)$ under <Ref>(ii),  where $I(\Omega)$ is the indicator function of an event $\Omega$
with $\P(\Omega)\to 1$.

Furthermore, if $T = o(n)$ as $n, p\to \infty$ while $\tau', \gamma, c$ stay constant,

\begin{equation}
\frac{\trace(\bS) + \fnorm{\bH}^2}{\fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n} \overset{p}{\to} 1.
\label{eq: convergence proba generalization error}
\end{equation}

In the above theorem, $\bS$ and $\bH$ are unknown, while $\hbA$ and $\bF$ can be computed from the observed data $(\bX, \bY)$.
Thus (<ref>) shows that $\fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n$ is a consistent estimate for
the unobserved quantity $\trace(\bS) + \fnorm{\bH}^2$.

§.§ Back to bounds on estimation error

We are now ready to present our main result on the error bounds for $\hbS$.
It is a consequence of (<ref>), (<ref>) and (<ref>).

Let <Ref> be fulfilled and $T = o(n)$. Then

\begin{align}
\fnorm{\hbS - \bS} &\le O_P(n^{-\frac 12}) (\fnorm*{\bF}^2/n),\\
\fnorm{\hbS - \bS} &\le O_P(n^{-\frac 12}) [\trace(\bS) + \fnorm{\bH}^2].
\label{eq:upper bound trS+H2}
\end{align}

Here the $O_P(n^{-\frac12})$ notation involves constants depending on $\gamma,\tau',c$.

It is instructive at this point to compare (<ref>)
with the lower bound (<ref>) on the Frobenius error
of the method-of-moments estimator. When $p\ge n$ then $\E[\fnorm{\hbS_{\rm{(mm)}}-\bS}^2]\ge \frac{c}{n}[\trace[\bS] + \fnorm{\bSigma^{1/2}\bB^*}^2]^2$;
this is the situation where the Statistician does not attempt to estimate
$\bB^*$, and pays a price of $[\trace[\bS] + \fnorm{\bSigma^{1/2}\bB^*}^2]^2/n$.
On the other hand, by definition of $\bH$ in (<ref>), the right-hand side of (<ref>),
when squared, is of order
$n^{-1}[\trace[\bS] + \fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2]^2$.

Here the error bound only depends on $\bB^*$ through
the estimation error for the nuisance $\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$.
This explains that when $\hbB$ is a good estimator of $\bB^*$
and $\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$ is smaller compared to $\fnorm{\bSigma^{1/2}\bB^*}^2$, the estimator $\hbS$ that leverages $\hbB$ will outperform
the method-of-moments estimator $\hbS_{\rm{(mm)}}$ which does not attempt to estimate the nuisance.

Finally, the next results show that under additional assumptions,
the estimator $\hbS$ enjoys Frobenius error bounds similar to the oracle estimator $\frac1n\bE^\top\bE$.

$\text{SNR }\le \mathfrak{snr}$ for some positive constant $\mathfrak{snr}$ independent of $n, p, T$, where  $\text{SNR}=  \fnorm{\bSigma^{\frac12}\bB^*}^2/\trace(\bS)$ denotes the signal-to-noise ratio of the multi-task linear model (<ref>).

Suppose that Assumptions <ref>, <ref>, <ref>,  <ref>(i), <ref> and $T=o(n)$ hold, then

\begin{equation}
\fnorm{\hbS - \bS} \le O_P(n^{-\frac 12}) \trace(\bS),
\end{equation}

where $O_P(\cdot)$ hides constants depending on $\gamma,\tau',\mathfrak{snr}$.

\begin{align*}
&\fnorm{\hbS - \bS}^2 \le O_P(T/n) \fnorm{\bS}^2 = o_P(1)\fnorm{\bS}^2,\\
&\big|\norm{\hbS}_* - \trace(\bS)\big| \le O_P(\sqrt{T/n}) \trace(\bS) = o_P(1) \trace(\bS).
\end{align*}

Suppose that Assumptions <ref>, <ref>, <ref>, <ref>(ii)
and $T = o(n)$ hold.
If $\norm{\bB^*}_0\le (1-c)n/2$ and the tuning parameter $\lambda$ is of the form $\lambda=\mu\sqrt{\trace(\bS)/n}$ for some positive constant $\mu$,

\begin{equation}
\fnorm{\hbS - \bS} \le O_P(n^{-\frac12}) (1 + \mu^2) \trace(\bS),
\label{eq: cor37}
\end{equation}

where $O_P(\cdot)$ hides constants depending on $c,\gamma, \phi_{\min}(\bSigma)$.

Comparing <Ref> with <Ref>, we conclude
that $\fnorm{\hbS-\bS}^2$ is of the same order as the Frobenius error of the oracle estimator
in (<ref>) up to constants depending on the signal-to-noise ratio, $\gamma$, and $\tau'$ under <Ref>(i), and up to constants depending on $\mu$, $c, \gamma, \phi_{\min}(\bSigma)$
under <Ref>(ii).

The error bounds in (<ref>)-(<ref>) are measured in Frobenius norm, similarly to existing works on
noise covariance estimation [Molstad, 2019].
Outside the context of linear regression models, much work has been devoted to covariance estimation
in the operator norm.
By the loose bound $\opnorm*{\bM}\leq \fnorm*{\bM}$, our upper bounds carry over to the operator norm.
The same cannot be said for lower bounds, since for instance
$\E\big[\opnorm{\hbS_{\rm{(oracle)}}- \bS}^2 \big] \asymp n^{-1} \opnorm{\bS} \trace(\bS)$ (see, , <cit.>).

[Boxplot for estimating a full rank $\bS$ ]

[Boxplot for estimating a low rank $\bS$ ]
Boxplots for Frobenius norm loss over 100 repetitions.

§ NUMERICAL EXPERIMENTS

Regarding parameters for our simulations,
we set $T=20$, $p=1.5n$ and $n$ equals successively $1000,1500,2000$.
We consider two settings for the noise covariance matrix:
$\bS$ is full-rank and $\bS$ is low-rank.
The complete construction of $\bS$, $\bB^*$ and $\bX$,
as well as implementation details are given in the supplementary material.

We compare our proposed estimator $\hbS$ with relevant estimators including
(1) the naive estimate $\hbS_{\text{(naive)}} = n^{-1}\bF^\top\bF$,
(2) the method-of-moments estimate $\hbS_{\rm{(mm)}}$ defined in Proposition <ref>, and
(3) the oracle estimate $\hbS_{\rm{(oracle)}} = n^{-1}\bE^\top\bE$.
The performance of each estimator is measured in Frobenius norm: for instance,  $\fnorm{\hbS-\bS}$ is the loss for proposed estimator $\hbS$.
<Ref> displays the boxplots of the Frobenius loss from the different methods over 100 repetitions.
<Ref> shows that, besides the oracle estimator, our proposed estimator has the best performance with significantly smaller loss  compared to the naive and method-of-moments estimators.

Since the estimation target is a $T\times T$ matrix, we also want to compare different estimators in terms of the bias and standard deviation for each entry of $\bS$.
<Ref> presents the heatmaps of bias and standard deviation from different estimators for full rank $\bS$ with $n=1000$.
The remaining heatmaps for different $n$ and for estimation of low rank $\bS$ are available in the supplementary material.

As expected, the oracle estimator has best performance in <Ref>
and smallest bias and variance in <Ref>.

The naive estimator has large bias as we see in <Ref>,
though it has small standard deviation.
The method-of-moments estimator is unbiased but its variance is relatively large,
which means its performance is not stable, as was reflected in <Ref>.
Our proposed estimator improves on both the naive and method-of-moments estimators
because it has much smaller bias than the former,
while having smaller standard deviation than the latter.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation for each entry ]
Heatmaps for estimation of full rank $\bS$ with $n=1000$ over 100 repetitions.

§ LIMITATIONS AND FUTURE WORK

One limitation of the proposed estimator $\hbS$ is
that its construction necessitates the knowledge of $\bSigma$.
Let us first mention that the estimator
$n^{-1}\fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}$ of $\trace(\bS) + \fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$ in <Ref> does not require knowing $\bSigma$. Thus, this estimator can further be used as a proxy of the error $\fnorm{\bSigma^{1/2}(\hbB-\bB^*)}^2$,
say for parameter tuning, without the knowledge of $\bSigma$.
The problem of estimating $\bS$ with known
$\bSigma$ was studied in [Celentano and Montanari, 2021] for $T=2$: in this
inaccurate covariate model
and for $p/n\le \gamma$, our results yield the convergence rate $n^{-1/2}$
for $\bS$ which improves upon the rate $n^{-c_0}$ for a non-explicit constant $c_0>0$ in <cit.>.

In order to use $\hbS$ when $\bSigma$ is unknown, one may plug-in
an estimator $\hbSigma$ in <Ref>, resulting in an extra term
of order $\opnorm{\hbSigma{}^{-1} - \bSigma^{-1}}\fnorm{\bF}$ for the Frobenius error. See <cit.> for related discussions in the $T=1$ (single-task) case. While, under the proportional regime $p/n\to \gamma$,
no estimator is consistent for all covariance matrices $\bSigma$ in operator norm, consistent estimators do exist under additional structural assumptions [Bickel and Levina, 2008, El Karoui, 2008, Cai et al., 2010].
If available, additional unlabeled samples $(\bx_i)_{i\ge n+1}$ can also
be used to construct norm-consistent estimator of $\bSigma$.

Future directions include extending estimator $\hbS$ to utilize

other estimators of the nuisance $\bB^*$
than the multi-task elastic-net (<ref>); for instance
(<ref>) or the estimators studied in [van de Geer and Stucky, 2016, Molstad, 2019, Bertrand et al., 2019]. In the simpler case where
columns of $\bB^*$ are estimated independently on each task, e.g.,
if the $T$ columns of $\hbB$ are Lasso estimators
$(\hbbeta^{(t)})_{t\in[T]}$ each computed from $\by{}^{(t)}$,
then minor modifications of our proof yield that the estimator (<ref>)
with $\hbA=\diag(\|\hbbeta^{(1)}\|_0,...,\|\hbbeta^{(T)}\|_0)$ enjoys similar
Frobenius norm bounds of order $n^{-1/2}$.

[Bayati and Montanari, 2012]
Mohsen Bayati and Andrea Montanari.
The lasso risk for gaussian matrices.
IEEE Transactions on Information Theory, 580
(4):0 1997–2017, 2012.

[Bayati et al., 2013]
Mohsen Bayati, Murat A Erdogdu, and Andrea Montanari.
Estimating lasso risk and noise level.
In NIPS, volume 26, pages 944–952, 2013.

[Bellec and Kuchibhotla, 2019]
Pierre Bellec and Arun Kuchibhotla.
First order expansion of convex regularized estimators.
In Advances in Neural Information Processing Systems, pages
3462–3473, 2019.

[Bellec, 2020]
Pierre C Bellec.
Out-of-sample error estimate for robust m-estimators with convex
arXiv preprint arXiv:2008.11840, 2020.

[Bellec and Romon, 2021]
Pierre C Bellec and Gabriel Romon.
Chi-square and normal inference in high-dimensional multi-task
arXiv preprint arXiv:2107.07828, 2021.

[Bellec and Tsybakov, 2017]
Pierre C Bellec and Alexandre B Tsybakov.
Bounds on the prediction error of penalized least squares estimators
with convex penalty.
In Vladimir Panov, editor, Modern Problems of Stochastic
Analysis and Statistics, Selected Contributions In Honor of Valentin
Konakov. Springer, 2017.
URL <https://arxiv.org/pdf/1609.06675.pdf>.

[Bellec and Zhang, 2019]
Pierre C Bellec and Cun-Hui Zhang.
De-biasing convex regularized estimators and interval estimation in
linear models.
arXiv preprint arXiv:1912.11943, 2019.

[Bellec and Zhang, 2021]
Pierre C Bellec and Cun-Hui Zhang.
Second-order stein: Sure for sure and other applications in
high-dimensional inference.
The Annals of Statistics, 490 (4):0
1864–1903, 2021.

[Belloni et al., 2011]
Alexandre Belloni, Victor Chernozhukov, and Lie Wang.
Square-root lasso: pivotal recovery of sparse signals via conic
Biometrika, 980 (4):0 791–806, 2011.

[Belloni et al., 2014]
Alexandre Belloni, Victor Chernozhukov, and Lie Wang.
Pivotal estimation via square-root lasso in nonparametric regression.
Ann. Statist., 420 (2):0 757–788, 04 2014.
URL <https://doi.org/10.1214/14-AOS1204>.

[Bertrand et al., 2019]
Quentin Bertrand, Mathurin Massias, Alexandre Gramfort, and Joseph Salmon.
Handling correlated and repeated measurements with the smoothed
multivariate square-root lasso.
In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural
Information Processing Systems, volume 32. Curran Associates, Inc., 2019.

[Bickel and Levina, 2008]
Peter J Bickel and Elizaveta Levina.
Covariance regularization by thresholding.
The Annals of Statistics, 360 (6):0
2577–2604, 2008.

[Boucheron et al., 2013]
Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.
Concentration inequalities: A nonasymptotic theory of
Oxford University Press, 2013.

[Cai et al., 2010]
T Tony Cai, Cun-Hui Zhang, and Harrison H Zhou.
Optimal rates of convergence for covariance matrix estimation.
The Annals of Statistics, 380 (4):0
2118–2144, 2010.

[Celentano and Montanari, 2021]
Michael Celentano and Andrea Montanari.
Cad: Debiasing the lasso with inaccurate covariate model.
arXiv preprint arXiv:2107.14172, 2021.

[Davidson and Szarek, 2001]
Kenneth R Davidson and Stanislaw J Szarek.
Local operator theory, random matrices and banach spaces.
Handbook of the geometry of Banach spaces, 10
(317-366):0 131, 2001.

[Dicker, 2014]
Lee H Dicker.
Variance estimation in high-dimensional linear models.
Biometrika, 1010 (2):0 269–284, 2014.

[Dicker and Erdogdu, 2016]
Lee H Dicker and Murat A Erdogdu.
Maximum likelihood for variance estimation in high-dimensional linear
In Artificial Intelligence and Statistics, pages 159–167.
PMLR, 2016.

[Dobriban and Wager, 2018]
Edgar Dobriban and Stefan Wager.
High-dimensional asymptotics of prediction: Ridge regression and
The Annals of Statistics, 460 (1):0 247–279,

[El Karoui, 2008]
Noureddine El Karoui.
Operator norm consistent estimation of large-dimensional sparse
covariance matrices.
The Annals of Statistics, 360 (6):0
2717–2756, 2008.

[Fan et al., 2012]
Jianqing Fan, Shaojun Guo, and Ning Hao.
Variance estimation using refitted cross-validation in ultrahigh
dimensional regression.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 740 (1):0 37–65, 2012.

[Fourdrinier et al., 2021]
Dominique Fourdrinier, Anis M Haddouche, and Fatiha Mezoued.
Covariance matrix estimation under data–based loss.
Statistics & Probability Letters, page 109160, 2021.

[Friedman et al., 2010]
Jerome Friedman, Trevor Hastie, and Rob Tibshirani.
Regularization paths for generalized linear models via coordinate
Journal of statistical software, 330 (1):0 1,

[Janson et al., 2017]
Lucas Janson, Rina Foygel Barber, and Emmanuel Candes.
Eigenprism: inference for high dimensional signal-to-noise ratios.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 790 (4):0 1037–1065, 2017.

[Koltchinskii and Lounici, 2017]
Vladimir Koltchinskii and Karim Lounici.
Concentration inequalities and moment bounds for sample covariance
Bernoulli, 230 (1):0 110–133, 2017.

[Laurent and Massart, 2000]
B. Laurent and P. Massart.
Adaptive estimation of a quadratic functional by model selection.
Ann. Statist., 280 (5):0 1302–1338, 10 2000.
URL <http://dx.doi.org/10.1214/aos/1015957395>.

[Liu and Zhang, 2009]
Han Liu and Jian Zhang.
Estimation consistency of the group lasso and its applications.
In Artificial Intelligence and Statistics, pages 376–383.
PMLR, 2009.

[Lounici et al., 2011]
Karim Lounici, Massimiliano Pontil, Sara Van De Geer, and Alexandre B Tsybakov.
Oracle inequalities and optimal inference under group sparsity.
The annals of statistics, 390 (4):0
2164–2204, 2011.

[Miolane and Montanari, 2018]
Léo Miolane and Andrea Montanari.
The distribution of the lasso: Uniform control over sparse balls and
adaptive parameter tuning.
arXiv preprint arXiv:1811.01212, 2018.

[Molstad, 2019]
Aaron J Molstad.
New insights for the multivariate square-root lasso.
arXiv preprint arXiv:1909.05041, 2019.

[Obozinski et al., 2011]
Guillaume Obozinski, Martin J Wainwright, and Michael I Jordan.
Support union recovery in high-dimensional multivariate regression.
The Annals of Statistics, 390 (1):0 1–47,

[Pedregosa et al., 2011]
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos,
D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.
Scikit-learn: Machine learning in Python.
Journal of Machine Learning Research, 12:0 2825–2830,

[Simon et al., 2013]
Noah Simon, Jerome Friedman, and Trevor Hastie.
A blockwise descent algorithm for group-penalized multiresponse and
multinomial regression.
arXiv preprint arXiv:1311.6529, 2013.

[Stein, 1981]
Charles M Stein.
Estimation of the mean of a multivariate normal distribution.
The annals of Statistics, pages 1135–1151, 1981.

[Sun and Zhang, 2012]
Tingni Sun and Cun-Hui Zhang.
Scaled sparse linear regression.
Biometrika, 990 (4):0 879–898, 2012.

[van de Geer and Stucky, 2016]
Sara van de Geer and Benjamin Stucky.
$\chi$ 2-confidence sets in high-dimensional regression.
In Statistical analysis for high-dimensional data, pages
279–306. Springer, 2016.

[Yu and Bien, 2019]
Guo Yu and Jacob Bien.
Estimating the error variance in a high-dimensional linear model.
Biometrika, 1060 (3):0 533–546, 2019.

§ SUPPLEMENT

This supplement is organized as follows:

* In <Ref> we provide details regarding the setting of our simulations, as well as additional experiment results.

* In <Ref> we establish an upper bound for the out-of-sample error, which we could not put in the full paper due to page length limit.

* <Ref> provides the upper bound for $\opnorm{(\bI_T - \hbA/n)^{-1}}$ mentioned after <Ref> in the full paper, and
Appendices <ref>, <ref>, <ref> contain preliminary theoretical statements, which will be useful for proving our main results in the full paper.

* <Ref> contains proofs of main results in <Ref> of the full paper and <Ref>.

* <Ref> contains proofs of preliminary results in Appendices <ref> to <ref>.

Here we introduce basic notations that will be used throughout this supplement.
We use indexes $i$ and $l$ only to loop or sum over $[n] =  \{1, 2, \ldots, n\}$, use $j$ and $k$ only to loop or sum over $[p] =  \{1, 2, \ldots, p\}$,
use $t$ and $t'$ only to loop or sum over $[T] =  \{1, 2, \ldots, T\}$, so that $\be_i$ (and $\be_l$) refer to the $i$-th (and $l$-th) canonical basis vector in $\R^n$, $\be_j$ (and $\be_k$) refer to the $j$-th (and $k$-th) canonical basis vector in $\R^p$, $\be_t$ (and $\be_{t'}$) refer to the $t$-th (and $t'$-th) canonical basis vector in $\R^T$.
For any two real numbers $a$ and $b$, let $a\vee b = \max(a,b)$, and $a\wedge b = \min(a,b)$.
Positive constants that depend on $\gamma, \tau'$ only are denoted by $C(\gamma, \tau')$, and positive constants that depend on $\gamma, c$ only are denoted by $C(\gamma, c)$. The values of these constants
may vary from place to place.

§ EXPERIMENT DETAILS AND ADDITIONAL SIMULATION RESULTS

§.§ Experimental details

This section provides more experimental detail for <Ref> of the full paper.

We consider two types of noise covariance matrix:
(i) $\bS$ is full-rank with $(t,t')$-th entry $\bS_{t,t'} = \frac{\cos(t-t')}{1 + \sqrt{|t-t'|}}$;
(ii) $\bS$ is low-rank with $\bS = \bu\bu^\top$, where $\bu\in \R^{T\times 10}$ has entries from $\calN(0, 1/T)$.

To build the coefficient matrix $\bB^*$, we first set its sparsity pattern, , we define the support $\mathscr{S}$ of cardinality $|\mathscr{S}| = 0.1 p$,
then we generate an intermediate matrix $\bB\in \R^{p\times T}$.
The $j$-th row of $\bB$ is sampled from $\calN_T(\mathbold0, p^{-1}\bI_T)$ if $j\in \mathscr{S}$,
otherwise we set it to be the zero vector.
Finally we let $\bB^* =\bB [\trace(\bS)/ \trace(\bB^\top\bSigma \bB)]^{\frac 12}$,
which forces a signal-to-noise ratio of exactly $1$.

The design matrix $\bX$ is constructed by independently sampling its rows from $\calN_p(\mathbold 0, \bSigma)$ with $\bSigma_{jk} = |j-k|^{0.5}$.

The Python library Scikit-learn [Pedregosa et al., 2011] is used to calculate $\hbB$ in (<ref>).
More precisely we invoke  to
obtain $\hbB$ by 5-fold cross-validation with parameters , .
To compute the interaction matrix $\hbA$
we used the efficient implementation described in <cit.>.

The full code needed to reproduce our experiments is part of the supplementary material.
A detailed Readme file is located in the corresponding folder.

The simulations results reported in the full paper and this supplementary material were run on a cluster of 50 CPU-cores (each is an Intel Xeon E5-2680 v4 @2.40GHz) equipped with a total of 150 GB of RAM.
It takes approximately six hours to obtain all of our simulation results.

§.§ Numerical results of Frobenius norm loss

While <Ref> in the full paper provides boxplots of Frobenius norm loss for 100 repetitions,
we report in following <Ref> the mean and standard deviation of the Frobenius norm loss over 100 repetitions.

Comparison of different methods for estimation of $\bS$

2c|$n=1000$     2c|$n=1500$      2c$n=2000$

$\bS$    method     mean     sd     mean     sd     mean     sd

4*full rank
naive  			   2.593    0.090       2.572   0.076    2.562    0.070

mm     			    2.030    0.616       1.554    0.421      1.413    0.405

proposed	    1.207    0.119       0.984    0.084      0.858    0.072

oracle 			    0.652    0.061      0.534    0.052      0.469    0.045

4*low rank
naive  			   2.942    0.119       2.912   0.102    2.908    0.094

mm     			    2.027    0.686       1.561    0.435      1.423    0.414

proposed	    1.216    0.172       0.989    0.125      0.854    0.118

oracle 			    0.654    0.096      0.531    0.081      0.464    0.065

The numerical results in  <Ref> are consistent with the boxplots in <Ref>.
It is clear from <Ref> that our proposed estimator has significantly smaller loss than the  naive estimator and method-of-moments estimator.
These comparisons again show the superior performance of our proposed estimator.

§.§ Additional heatmaps for estimating full rank S

When estimating the full rank $\bS$ with $(t,t')$-th entry $\bS_{t,t'} = \frac{\cos(t-t')}{1 + \sqrt{|t-t'|}}$, the heatmaps for different estimators from $n=1500$ and $n=2000$ are presented in <Ref>, respectively.
The comparison patterns in <Ref> are similar to the case $n=1000$ in <Ref> of the full paper; our proposed estimator outperforms the naive estimator and method-of-moments estimator.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation in each entry ]
Heatmaps for estimation of full rank $\bS$ with $n=1500$ over 100 repetitions.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation in each entry ]
Heatmaps for estimation of full rank $\bS$ with $n=2000$ over 100 repetitions.

§.§ Heatmaps for estimating low rank S

When estimating the low rank with $\bS = \bu\bu^\top$, and $\bu\in \R^{T\times 10}$ with entries are from $\calN(0, 1/T)$.
We present the heatmaps for different estimators with $n=1000, 1500, 2000$ in <Ref> below. All of these figures convince us that besides the oracle estimator, the proposed estimator has the best performance.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation in each entry ]
Heatmaps for estimation of low rank $\bS$ with $n=1000$ over 100 repetitions.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation in each entry ]
Heatmaps for estimation of low rank $\bS$ with $n=1500$ over 100 repetitions.

[Heatmap of bias for each entry ]

[Heatmap of standard deviation in each entry ]
Heatmaps for estimation of low rank $\bS$ with $n=2000$ over 100 repetitions.

§ OUT-OF-SAMPLE ERROR ESTIMATION

In this appendix, we present a by-product of our techniques for estimating the noise covariance.
Suppose we wish to evaluate the performance of a regression method on a new data, we define the out-of-sample error for the multi-task linear model (<ref>) as

\E \big[(\hbB - \bB^*)^\top\bx_{\text{new}} \bx_{\text{new}}^\top (\hbB - \bB^*)|(\bX, \bY)\big] = \bH^\top \bH,

where $\bx_{\text{new}}$ is independent of the data $(\bX; \bY)$ with the same distribution as any row of $\bX$.

The following theorem on estimation of out-of-sample error is an by-product of our technique for constructing $\hbS$.

Under the same conditions of <Ref>, with $\bZ = \bX \bSigma^{-\frac12}$, we have

\begin{align*}
(\bI_T -\hbA/n)\bH^\top\bH (\bI_T -\hbA/n) -  \frac{1}{n^2} \big( \bF^\top \bZ \bZ^\top \bF + \hbA \bF^\top \bF + \bF^\top \bF \hbA -
p\bF^\top \bF\big)}\\
\le &  \Theta_3 n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)
\end{align*}

for some non-negative random variable $\Theta_3$ of constant order, in the sense that $\E [\Theta_3] \le C(\gamma,\tau')$ under <Ref>(i),
and with ${\E [I(\Omega)\Theta_3]} \le C(\gamma,c)$ under <Ref>(ii),  where $I(\Omega)$ is the indicator function of an event $\Omega$
with $\P(\Omega)\to 1$.

<Ref> generalizes the result in [Bellec, 2020] to multi-task setting.
While the out-of-sample error $\bH^\top\bH$ is unknown, the quantities $\bZ$, $\bF$, $\hbA$ are observable. Since typically the quantity $(\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)$ is of a constant order, <Ref> suggests the following estimate of $\bH^\top\bH$:

\[
\frac{1}{n^2} (\bI_T - \hbA/n)^{-1} \big( \bF^\top \bZ\bZ^\top \bF + \hbA \bF^\top \bF + \bF^\top \bF \hbA -
p\bF^\top \bF\big) (\bI_T - \hbA/n)^{-1},
\]

which can further be used for parameter tuning in multi-task linear model.

We present the proof of <Ref> in <Ref>

§ USEFUL OPERATOR NORM BOUNDS

Let us first introduce two events besides the event $U_1 = \big\{ \norm{\hbB}_0 \le n(1-c)/2 \big\}$ in <Ref>(ii), we define events $U_2$ and $U_3$ as below,

\begin{align*}
%U_1 &= \big\{ \norm{\hbB}_0 \le n(1-c)/2 \big\},\\
U_2 &= \big\{\inf_{\bb\in \R^p: \| \bb\|_0 \le (1-c)n}  \|\bX \bb\|^2/(n \|\bSigma^{\frac 12} \bb\|^2)  > \eta\big\},\\
U_3 &= \big\{ \opnorm{\bX\bSigma^{-\frac 12}} < 2 \sqrt n  + \sqrt p\big\}.
%~ \opnorm{\bE\bS^{-\frac 12}} < 2 \sqrt n  + \sqrt T \big\}.
%U_4 &= \big\{\bE : \opnorm{\bE\bS^{-\frac 12}} < 2 \sqrt n  + \sqrt T \big\}
\end{align*}

Under <Ref>, <cit.> guarantees $\P(U_2) \ge 1 - C(\gamma, c)e^{-C(\gamma, c)n}$ for some constant $\eta$ that only depends on constants $\gamma, c$.
Under <Ref>, <cit.> guarantees $\P(U_3) > 1 - e^{-n/2}$ and there exists a random variable $z\sim \calN(0,1)$ s.t. $\opnorm{\bX\bSigma^{-\frac 12}} \le \sqrt{n} + \sqrt{p} + z$ almost surely. Therefore, under <Ref>, we have

\begin{equation}\label{eq: opnorm-Z}
\E[ \opnorm{n^{-\frac 12}\bX\bSigma^{-\frac 12}}^2] \le (1 + \sqrt{p/n})^2 + n^{-1} \le C(\gamma).
\end{equation}

Furthermore, under <Ref>(ii),  $\P(U_1\cap U_2\cap U_3)\to1$ by a union bound, and for large enough $n$,

\begin{align}
\P\big\{(U_1\cap U_2\cap U_3)^c\big\} \notag
<& 1/T + C(\gamma, c)e^{-n/C(\gamma, c)} \notag\\
=& \frac 1T (1 + TC(\gamma, c)e^{-n/C(\gamma, c)} )\notag\\
<& \frac 1T  (1 + C(\gamma, c)e^{\sqrt{n} -n/C(\gamma, c)})\notag \\
<&  \frac 1T C(\gamma, c). \label{eq: P-Omega}
\end{align}

Now we provide the operator norm bounds for $\bI_T - \hbA/n$ and $(\bI_T - \hbA/n)^{-1}$.

Suppose that <Ref> holds. If $\tau >0$ in (<ref>) with $\tau' = \tau /\opnorm{\bSigma}$, then

* $\opnorm{\bI_T - \hbA/n} \le 1$.

* In the event $U_3$, we have $\opnorm{(\bI_T - \hbA/n)^{-1}} \le 1 + (\tau')^{-1} (2 + \sqrt{p/n})^2$. Furthermore, $\E [\opnorm{(\bI_T - \hbA/n)^{-1}} ]\le 1 + (\tau')^{-1} [(1 + \sqrt{p/n})^2 +n^{-1}].$

Suppose that <Ref> holds.
If $\tau=0$ in (<ref>), then

* In the event $U_1$, we have $\opnorm{\bI_T - \hbA/n} \le 1$.

* In the event $U_1$, $\opnorm{(\bI_T - \hbA/n)^{-1}} \le C(c)$. Hence, $\E [I(U_1)\opnorm{(\bI_T - \hbA/n)^{-1}} ]\le C(c).$

With $\bN = (\bI_T \otimes \bX) \bM^\dagger (\bI_T \otimes \bX^\top)$, we have $\opnorm{\bN}\le 1$.

§ LIPSCHITZ AND DIFFERENTIAL PROPERTIES FOR A GIVEN, FIXED NOISE MATRIX E

We need to study Lipschitz and differential properties of certain mappings when the noise matrix $\bE$ is fixed.
$g:\R^{p\times T}\to \R$ defined by
$g(\bB)=\tau \fnorm{\bB}^2/2 + \lambda \norm{\bB}_{2,1}$ be the penalty in
For a fixed value of $\bE$,
define the mappings

\begin{align}
\bZ
\bH(\bZ) = \argmin_{\bar\bH\in\R^{p\times T}} \frac{1}{2n}\fnorm{\bE - \bZ\bar\bH}^2 + g(\bSigma^{-1/2}\bar\bH)
\R^{n\times p}
\to
\R^{p\times T}
\\
\bZ
\bF(\bZ) = \bE - \bZ\bH(\bZ)
\R^{n\times p}
\to
\R^{n\times T}
\\
\bZ
D(\bZ) = (\fnorm{\bH(\bZ)}^2 + \fnorm{\bF(\bZ)}^2/n)^{1/2}
\R^{n\times p}
\to
\R
\end{align}

Next, define the random variable $\bZ = \bX\bSigma^{-\frac12}\in\R^{n\times p}$, and let us use the convention that if arguments of $\bH,\bF$ or $D$ are omitted then these mappings are implicitly taken at the realized value of the random variable $\bZ = \bX\bSigma^{-\frac12}\in\R^{n\times p}$ where $\bX$ is the observed design matrix. With this convention and by definition of the above mappings, we then have
$\bH = \bH(\bZ) = \bSigma^{1/2}(\hbB - \bB^*)$ as well as
$\bF = \bF(\bZ) = \bY - \bX\hbB$
and $D = [\fnorm{\bH}^2 + \fnorm{\bF}^2/n]^{1/2}$
so that the notation is consistent with the rest of the paper
(in particular, with (<ref>)).

Finally, denote the $(i,j)$-th entry of $\bZ$ by $z_{ij}$
throughout this appendix, and the corresponding partial derivatives of the above mappings
by $\frac{\partial}{\partial z_{ij}}$.

§.§ Lipschitz properties

For multi-task elastic-net (, $\tau>0$ in (<ref>)), the mapping $\bZ \mapsto D^{-1}\bF/\sqrt{n}$ is $n^{-\frac12}L$-Lipschitz with $L = 8 \max(1, (2\tau')^{-1})$, where $\tau' = \tau/\opnorm{\bSigma}$.

For multi-task group Lasso (, $\tau=0$ in (<ref>)).

we have

(1) In $U_1\cap U_2$, the map $\bZ \mapsto D^{-1}\bF/\sqrt{n}$ is $n^{-\frac12}L$-Lipschitz with $L = 8 \max(1, (2\eta)^{-1})$.

(2) In $U_1\cap U_2 \cap U_3$,  the map $\bZ \mapsto D^{-1}\bZ^\top\bF/n$ is $n^{-1/2} (1 + (2 +\sqrt{p/n})L)$-Lipschitz, where $L = 8 \max(1, (2\eta)^{-1})$ as in (1).

Suppose that <Ref> holds, then

(1) Under <Ref>(i) that $\tau>0$ and $\tau'=\tau/\opnorm{\bSigma}$, we have

\begin{align*}
\sum_{ij}\Big(\frac{\partial D}{\partial z_{ij}}\Big)^2 \le n^{-1}D^2 [4 \max(1, (2\tau')^{-1})]^2.
\end{align*}

This implies that
$ nD^{-2}\sum_{ij}(\frac{\partial D}{\partial z_{ij}})^2 \le C(\tau').

(2) Under <Ref>(ii) that $\tau=0$ and $\P(U_1)\to 1$, in the event $U_1\cap U_2$, we have

\begin{align*}
\sum_{ij}\Big(\frac{\partial D}{\partial z_{ij}}\Big)^2 \le n^{-1}D^2 [4 \max(1, (2\eta)^{-1})]^2.
\end{align*}

This implies that
$ nD^{-2}\sum_{ij}(\frac{\partial D}{\partial z_{ij}})^2 \le C(\eta) = C(\gamma, c)$ since $\eta$ is a constant that only depends on $\gamma, c$.

§.§ Derivative formulae

Note that with a fixed noise $\bE$, <Ref> guarantee that the map $\bZ \mapsto \bF$ is Lipschitz, hence the derivative exists almost everywhere by Rademacher’s theorem.
We present the formula for derivative of this map in <Ref>.

Recall $\bF = \bY - \bX\hbB$ with $\hbB$ defined in (<ref>). Under <Ref>(i) $\tau>0$, or under <Ref>(ii) $\tau =0$ and in the event $U_1\cap U_2$, for each $i,l\in [n], j\in [p], t\in [T]$, the following derivative exists almost everywhere and has the expression

\begin{equation*}
\frac{\partial F_{lt}}{\partial z_{ij}} = D_{ij}^{lt} + \Delta_{ij}^{lt},
\end{equation*}

= -(\be_j^\top\bH \otimes \be_i^\top) (\bI_{nT} - \bN) (\be_t\otimes \be_l),$ and
= -(\be_t^\top \otimes \be_l^\top)(\bI_T\otimes \bX)
\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})
\bigl(\bF^\top \otimes \bI_{p}\bigr)(\be_i \otimes\be_j).
Furthermore, a straightforward calculation yields

$$\sum_{i=1}^n D_{ij}^{it} = -\be_j^\top \bH (n\bI_T - \hbA)\be_t.$$

Suppose that <Ref> holds.

(1) Under <Ref>(i) that $\tau>0$ and $\tau' = \tau/\opnorm{\bSigma}$, we have

\begin{align*}
\frac 1n \sum_{ij}\norm*{\frac{\partial (\bF/D)}{\partial z_{ij}}}^2_{\rm F} \le \underbrace{4 \max(1, (\tau')^{-1} (T\wedge \frac{p}{n})) + 2  n^{-1}  [4 \max(1, (2\tau')^{-1})]^2}_{f(\tau', T, n, p)}.
\end{align*}

(2) Under <Ref>(ii) that $\tau=0$ and $\P(U_1)\to 1$, in the event $U_1\cap U_2$, we have

\begin{align*}
\frac 1n \sum_{ij}\norm*{\frac{\partial (\bF/D)}{\partial z_{ij}}}^2_{\rm F} \le \underbrace{4 \max(1, (\eta)^{-1} (T\wedge \frac{p}{n})) + 2  n^{-1}  [4 \max(1, (2\eta)^{-1})]^2}_{f(\eta, T, n, p)}.
\end{align*}

Furthermore, the right-hand side in (1) can be bounded from above by $C(\tau') (T\wedge \frac pn)$, and the  right-hand side in (2) can be bounded from above by $C(\gamma, c)$ in the regime $p/n\le \gamma$.

§ LIPSCHITZ AND DIFFERENTIAL PROPERTIES FOR A GIVEN, FIXED DESIGN MATRIX

We also need to study Lipschitz and derivative properties of functions
of the noise $\bE$ when the design $\bX$ is fixed. Formally, for a given and fixed design matrix $\bX$, define the function
$\bE\mapsto \bF(\bE)$ by
the value $\bY-\bX\hbB$ of the residual matrix
when the observed data $(\bX,\bY)$ is $(\bX,\bX\bB^* + \bE)$
and with $\hbB$ the estimator (<ref>).
Note that this map is 1-Lipschitz by <cit.>. Rademacher’s theorem thus guarantees this map is differentiable almost everywhere. We denote its partial derivative by $\frac{\partial}{\partial E_{it}}$
for each entry $(E_{it})_{i\in[n],t\in[T]}$ of the noise matrix $\bE$.
We present its derivative formula in <Ref> below.

For each $i,l\in [n], t,t'\in [T]$, the following derivative exists almost everywhere and has the expression

\begin{align*}
\frac{\partial F_{lt}}{\partial E_{it'}} = \be_l^\top\be_i\be_t^\top\be_{t'}
- \be_l^\top (\be_t^\top \otimes \bX)\bM^\dagger (\be_{t'} \otimes \bX^\top) \be_i.
\end{align*}

As a consequence, we further have

\begin{equation*}
\sum_{i=1}^n \frac{\partial F_{it}}{\partial E_{it'}} = \be_t^\top (n\bI_T -\hbA)\be_{t'},\quad \sum_{i=1}^n \frac{\partial \be_i^\top\bZ \bH\be_{t}}{\partial E_{it'}} = \be_t^\top \hbA \be_{t'}.
\end{equation*}

§ PROBABILISTIC TOOLS

We first list some useful variants of Stein's formulae and Gaussian-Poincaré inequalities. Let $f'$ denote the derivative of a differentiable univariate function. For a differentiable vector-valued function $\bff(\bz): \R^n \to \R^n$, denote its Jacobian (derivative) and divergence respectively by $\nabla \bff(\bz)$ and $\div \bff(\bz)$, , $[\nabla \bff(\bz)]_{i,l} = \frac{\partial f_i(\bz)}{\partial z_l}$ for $i,l\in [n]$, and $\div \bff(\bz) =  \trace(\nabla \bff(\bz))$.

The following identities hold provided the involved derivatives exist a.e. and the expectations are finite.

* $z\sim \calN(0, 1)$, $f: \R \to \R$, then

$$\E [(z f(z) - f'(z))^2] = \E [f(z)^2] + \E[(f'(z))^2].$$

* $\bz \sim \calN_n(\mathbold 0, \bI_n)$, $f: \R^n \to \R^n$, then

\E [(\bz^\top \bff(\bz) - \div\bff(\bz))^2] = \E \big[\norm{\bff(\bz)}^2 + \trace[( \nabla\bff(\bz))^2]\big] \le  \E \big[\norm{\bff(\bz)}^2 + \fnorm{ \nabla\bff(\bz)}^2\big],

where the inequality uses Cauchy-Schwarz inequality.

* More generally, for $\bz \sim \calN_n(\mathbold 0, \bSigma)$, $\bff: \R^n \to \R^n$, then

\begin{align*}
\E [(\bz^\top \bff(\bz) - \trace(\bSigma \nabla\bff(\bz))^2]
&= \E \big[\norm{\bSigma^{\frac12}\bff(\bz)}^2 + \trace[(\bSigma \nabla\bff(\bz))^2]\big]\\
&\le \E \big[\norm{\bSigma^{\frac12}\bff(\bz)}^2 + \fnorm{(\bSigma \nabla\bff(\bz)}^2]\big],
\end{align*}

where the inequality uses Cauchy-Schwarz inequality.

The following inequalities hold provided the right-hand side derivatives exist a.e. and the expectations are finite.

* $z \sim \calN(0, 1)$, $f: \R \to \R$, then

$$\text{Var} [f(z)] \le \E [(f'(z))^2].$$

* $\bz \sim \calN_n(\mathbold 0, \bI_n)$, $f: \R^n \to \R$, then

$$\text{Var} [f(\bz)] \le \E [\norm{\nabla f(\bz)}^2].$$

* $\bz \sim \calN_n(\mathbold 0, \bI_n)$, $\bff: \R^n \to \R^m$, then

$$\E[\norm{\bff(\bz) - \E[\bff(\bz)]}^2] \le \E [\fnorm{\nabla \bff(\bz)}^2].$$

* More generally, for $\bz \sim \calN_n(\mathbold 0, \bSigma)$, $\bff: \R^n \to \R^m$, then

$$\E[\norm{\bff(\bz) - \E[\bff(\bz)]}^2] \le \E [\fnorm{\bSigma^{\frac12}\nabla \bff(\bz)}^2].$$

Now we present a few important lemmas, whose proofs are based on <Ref> and <Ref>.

Assume that Assumption <ref> holds. For fixed $\bX$, we have

\begin{equation*}
\E\Bigl[
\fnorm{\bE^\top \bF/\tD - \bS (n\bI_T - \hbA )/\tD}^2\Bigr]
\le 4\trace(\bS),
\end{equation*}

where $\tD = \big(\fnorm*{\bF}^2 + n \trace(\bS)\big)^{\frac 12}$.

Let $\bU, \bV:\R^{n\times p} \to \R^{n\times T}$ be two locally Lipschitz functions of $\bZ$ with $\calN(0,1)$ entries, then

\begin{align*}
&\E\Big[\norm[\Big]{\bU^\top \bZ \bV -
\sum_{j=1}^p\sum_{i=1}^n
\frac{\partial}{\partial z_{ij} }\Bigl(\bU^\top \be_i \be_j^\top \bV \Bigr)
}_{\rm F}^2\Big]\\
\le~& \E \fnorm*{\bU}^2 \fnorm*{\bV}^2+ \E \sum_{ij}\Big[
2\fnorm*{\bV}^2\fnorm*{ \frac{\partial \bU}{\partial z_{ij}} }^2
+ 2\fnorm*{\bU}^2\fnorm*{ \frac{\partial \bV}{\partial z_{ij}} }^2\Big].
\end{align*}

Assume the same setting as <Ref>.
If on some open set $\Omega\subset \R^{n\times p}$ with $\P(\Omega^c)\le C/T$ for some constant $C$, we have (i) $\bU$ is $n^{-1/2} L_1$ -Lipschitz and $\fnorm{\bU}\le 1$, (ii) $\bV$ is $n^{-1/2}L_2$ -Lipschitz and $\fnorm{\bV}\le K$. Then

\begin{align*}
&\E\Big[I(\Omega) \Big\|\bU^\top \bZ \bV -
\sum_{j=1}^p\sum_{i=1}^n
\frac{\partial}{\partial z_{ij} }\Bigl(\bU^\top \be_i \be_j^\top \bV \Bigr)
\Big\|_{\rm F}^2\Big]\\
\le~& K^2 + 2C( K^2 L_1^2  + L_2^2 ) + 2\E \Big[I(\Omega)
\sum_{ij}\Big( K^2\fnorm*{\frac{\partial \bU}{\partial z_{ij}} }^2 + \fnorm*{\frac{\partial \bV}{\partial z_{ij}} }^2\Big) \Big].
\end{align*}

Let $\bU,\bV: \R^{n\times p} \to \R^{n\times T}$
be two locally Lipschitz
functions of $\bZ$ with $\calN(0,1)$ entries.
Assume also that $\fnorm{\bU} \vee \fnorm{\bV} \le 1$ almost surely.

\begin{align*}
\fnorm[\Big]{
p \bU^\top \bV - \sum_{j=1}^p
\Bigl(\sum_{i=1}^n \partial_{ij} \bU^\top \be_i - \bU^\top \bZ \be_j\Bigr)
\Bigl(\sum_{i=1}^n \partial_{ij} \be_i^T \bV  - \be_j^T \bZ^T \bV\Bigr)
\Bigr]
\\
\le~&
2 \|\bU\|_\partial
\|\bV\|_\partial
\sqrt p
\bigl(
\sqrt 2 + (3+\sqrt{2})(\|\bU\|_{\partial} + \|\bV\|_{\partial})
\bigr),
\end{align*}

where $ \partial_{ij} \bU \defas \partial \bU /\partial z_{ij}$,
and $\|\bU\|_\partial \defas \E[\sum_{i=1}^n\sum_{j=1}^p
\fnorm{ \partial_{ij} \bU}^2]^{\frac 12}$.

Suppose that <Ref> holds. Let $\bQ_1 = \frac{\frac{1}{n}\big(\bF^\top\bF + \bH^\top\bZ^\top\bF - \bS(n\bI_T - \hbA) \big)}{\fnorm{\bS^{\frac 12}}(\fnorm*{\bF}^2/n + \trace(\bS))^{\frac 12} n^{-\frac 12} }$,
then  $\E [\fnorm*{\bQ_1}^2] \le 4$.

Suppose that <Ref> hold.
\[\bQ_2 = \frac{\frac{1}{n^2} \big( \bF^\top \bZ\bZ^\top \bF - \bF^\top\bF (p\bI_T - \hbA)+ (n\bI_T -\hbA) \bH^\top\bZ^\top\bF \big)}{(\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12}},
\]

then $\E [\fnorm*{\bQ_2}^2] \le C(\tau') (T\wedge (1+\frac{p}{n})) (1 + \frac pn)$ under <Ref>(i), and
$\E [I(\Omega)\fnorm*{\bQ_2}^2] \le C(\gamma, c)$ under <Ref>(ii) for some set $\Omega$ with $\P(\Omega)\to 1$.

Suppose that <Ref> hold.
Let $\Xi= (n\bI_T -\hbA)\bH^\top\bZ^\top\bF$, and

\[
\bQ_3 = \frac{\frac{1}{n^2} \big( p\bF^\top \bF -\bF^\top \bZ\bZ^\top \bF - (n\bI_T -\hbA)\bH^\top\bH (n\bI_T -\hbA) -  \Xi - \Xi^\top\big)}{(\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12}},
\]

then $\E [\fnorm*{\bQ_3}] \le C(\gamma,\tau')$ under <Ref>(i), and $\E [I(\Omega)\fnorm*{\bQ_3}] \le C(\gamma, c)$ under <Ref>(ii) for some set $\Omega$ with $\P(\Omega)\to 1$.

§ PROOFS OF MAIN RESULTS

In this appendix, we provide proofs of the theoretical results in <Ref> of the pull paper and  <Ref> of this supplement.

§.§ Proof of Propositionprop: oracle

With $\bS=\sum_{t=1}^T \sigma_t^2 \bu_t \bu_t^T$
the spectral decomposition of $\bS$,
we have $\fnorm{\bE^\top\bE-n\bS}^2
[\bu_{t'}^T((\bE^\top\bE-n\bS) \bu_t]^2$.
We now compute the expectation of one term
indexed by $(t,t')$.
The random variable $\bu_{t'}^T(\bE^\top\bE-n\bS) \bu_t$
is the sum of $n$ mean zero random variables
with the same distribution as $z_{t'} z_{t}-\bu_{t'}^T\bS\bu_t$
where $(z_t,z_{t'})\sim \calN_2(\mathbf{0},\diag(\sigma_t^2,\sigma_{t'}^2))$. Thus

\E[(\bu_{t'}^T(\bE^\top\bE-n\bS) \bu_t)^2]
n\text{Var}[z_{t'} z_{t}-\bu_{t'}^T\bS\bu_t]
n\sigma_t^2\sigma_{t'}^2 I_{t\ne t'}

due to $\text{Var}[\chi^{2}_1]=2$
if $t=t'$ and independence if $t\ne t'$.
Summing over all $(t,t')\in[T] \times [T]$
$2n\sum_{t=1}^T \sigma_t^4
+n \sum_{t\ne t'}\sigma_t^2\sigma_{t'}^2
=n\sum_{t=1}^T \sigma_t^4
= n \fnorm{\bS}^2
+ n [\trace(\bS)]^2
as desired.

The inequality simply follows from
$\fnorm{\bS}^2\le [\trace(\bS)]^2$ since $\bS$ is positive semi-definite.

§.§ Proof of Propositionprop: mom

Without of loss of generality, we assume $\bSigma = \bI_p$. For general positive definite $\bSigma$, the proof follows by replacing $(\bX, \bB^*)$ with $(\bX \bSigma^{-\frac12}, \bSigma^{\frac12}\bB^*)$.

We first derive the method-of-moments estimator $\hbS_{\rm{(mm)}}$. Under <Ref> with $\bSigma = \bI_p$, $\bX$ has rows from $\calN_p(\mathbf 0, \bI_p)$, $\bE$ has rows from $\calN_T(\bf0, \bS)$, and $\bX$ and $\bE$ are independent. Then, the expectations of $\bY^\top \bY$ and  $\bY^\top\bX \bX^\top\bY$ are given by

\begin{align}\label{eq:yy}
\E (\bY^\top \bY) &= \E \big[(\bX\bB^* + \bE)^\top (\bX\bB^* + \bE)\big] = n (\bB^{*\top}\bB^* + \bS),
\end{align}

\begin{align}\label{eq:yxxy}
\E (\bY^\top\bX\bX^\top\bY ) &= \E \big[(\bX\bB^* + \bE)^\top \bX \bX^\top (\bX\bB^* + \bE)\big]\notag\\
&= \E (\bB^{*\top}\bX^\top\bX \bX^\top\bX\bB^*) + \E (\bE^\top\bX \bX^\top\bE)\notag\\
&= \bB^{*\top}\E(\bX^\top\bX \bX^\top\bX)\bB^* + \E (\bE^\top\bX \bX^\top\bE)\notag\\
&= n(n+p+1)  \bB^{*\top}\bB^* + np\bS,
\end{align}

where the last line uses

\begin{align*}
&\E (\bX^\top\bX \bX^\top\bX) \\
=~& \E \Big[\sum_{i=1}^n (\bx_i\bx_i^\top) \sum_{l=1}^n(\bx_l\bx_l^\top)\Big]\\
=~& \sum_{i\ne l}\E(\bx_i\bx_i^\top \bx_l\bx_l^\top) + \sum_{i= l} \E (\bx_i\bx_i^\top \bx_l\bx_l^\top)\\
=~& n(n-1) \bI_p^2 + n\E (\bx_1\bx_1^\top \bx_1\bx_1^\top)\\
=~& n(n-1) \bI_p + n [2\bI_p^2 + \trace(\bI_p) \bI_p]\\
=~& n(n+p+1) \bI_p,
\end{align*}

\begin{align*}
\E (\bE^\top\bX \bX^\top\bE)
= \E\big[ \E (\bE^\top\bX \bX^\top\bE |\bE)\big]
= \E\big[ \bE^\top \E(\bX \bX^\top)\bE\big]
= np \bS.
%&= \E\big[ \bE^\top \trace(\bI_p) \bI_n \bE \big]\\
\end{align*}

Solving for $\bS$ from the system of equations (<ref>) and (<ref>), we obtain the method-of-moments estimator

\begin{align*}
\hbS_{\rm{(mm)}} = \frac{(n+p+1)}{n(n+1)} \bY^\top \bY - \frac{1}{n(n+1)}  \bY^\top\bX \bX^\top\bY,
\end{align*}

and $\E [\hbS_{\rm{(mm)}} ]= \bS$.

Now we derive the variance lower bound for $\hbS_{\rm{(mm)}}$. Since $\E[\hbS_{\rm{(mm)}} ] = \bS$, $\E \big[\fnorm{\hbS_{\rm{(mm)}} - \bS}^2\big] = \sum_{t, t'} \text{Var}\big\{[\hbS_{\rm{(mm)}} ]_{t,t'}\big\}.$ By definition of $\hbS_{\rm{(mm)}} $,

\begin{align*}
[\hbS_{\rm{(mm)}} ]_{t,t'} = \frac{n+p+1}{n(n+1)} [\by^{(t)}]^\top \by^{(t')}  - \frac{1}{n(n+1)}[\by^{(t)}]^\top \bX\bSigma^{-1}\bX^\top\by^{(t')}.
\end{align*}

Since $\by^{(t)} = \bX \bbeta^{(t)} + \bep^{(t)},\quad \by^{(t')} = \bX \bbeta^{(t')} + \bep^{(t')},$ for $t\ne t'$,

without loss of generality, we assume $\bbeta^{(t)} =  a_0\be_1$ and $\bbeta^{(t')} =  a_1\be_1 + a_2\be_2$ for some constants $a_0, a_1, a_2$.
If necessary, we could let $\bu_1 = \bbeta^{(t)}/\norm{\bbeta^{(t)}}$, and $\bu_2 = \tbu_2/\norm{\tbu_2}$ where $\tbu_2 =\bbeta^{(t')} - \bP_{\bu_1}\bbeta^{(t')}$, and completing the basis to obtain an orthonormal basis $\{\bu_1, \bu_2, \ldots, \bu_p\}$ for $\R^p$. Let $\bU =[\bu_1, \bu_2, \ldots, \bu_p]$, then $\bU$ is an orthogonal matrix, hence $\bX\bU$ and $\bX$ have the same distribution, only the first coordinate of $\bU^\top\bbeta^{(t)}$ is nonzero, and only the first two coordinates of $\bU^\top\bbeta^{(t')}$ are be nonzero. That is, we could perform change of variables by replacing $(\bX, \bbeta^{(t)}, \bbeta^{(t')})$ with $(\bX\bU, \bU^\top\bbeta^{(t)}, \bU^\top\bbeta^{(t')})$.

Therefore,  $\by^{(t)}$ and $\by^{(t')}$ are independent of $\{\bX\be_j: 3\le j\le p\}$.
Let $\calF = \sigma(\by^{(t)}, \by^{(t')}, \bX\be_1, \bX\be_2)$ be the $\sigma-$field generated by $(\by^{(t)}, \by^{(t')}, \bX\be_1, \bX\be_2)$, then

\begin{align*}
\text{Var}\big\{[\hbS_{\rm{(mm)}} ]_{t,t'}\big\}
&\ge \E\big[ \text{Var}\big\{[\hbS_{\rm{(mm)}} ]_{t,t'}|\calF\big\}\big]= \frac{1}{n^2(n+1)^2} \E\big[ \text{Var}\big\{ [\by^{(t)}]^\top \bX\bX^\top\by^{(t')}|\calF\big\}\big].
\end{align*}

Note that in the above display,

\begin{align*}
&[\by^{(t)}]^\top \bX\bX^\top\by^{(t')}
= \sum_{j=1}^2[\by^{(t)}]^\top \bX\be_j\be_j^\top\bX^\top\by^{(t')}
+ \sum_{j=3}^p[\by^{(t)}]^\top \bX\be_j\be_j^\top\bX^\top\by^{(t')},
\end{align*}

where the first term is measurable with respect to $\calF$, and the second term is a quadratic form

\begin{align*}
&\sum_{j=3}^p[\by^{(t)}]^\top \bX\be_j\be_j^\top\bX^\top\by^{(t')} = \sum_{j=3}^p \be_j^\top \bX^\top\by^{(t')} [\by^{(t)}]^\top \bX \be_j =  \bxi^\top \bLambda \bxi,
\end{align*}

here $\bxi = [\be_3^\top\bX^\top, \ldots, \be_p^\top\bX^\top]^\top\sim \calN(\mathbold{0}, \bI_{n(p-2)})$, and $\bLambda = \bI_{p-2} \otimes \by^{(t')} [\by^{(t)}]^\top$.
Thus, for $t\ne t'$,

\begin{align*}
\text{Var}\big\{[\hbS_{\rm{(mm)}} ]_{t,t'}\big\}
\ge~& \frac{1}{n^2(n+1)^2} \E \Big\{\text{Var}\big\{ \bxi^\top \bLambda \bxi|\calF\big\}\Big\}\\
=~& \frac{1}{n^2(n+1)^2} \E \Big\{ \fnorm{\bLambda}^2 + \trace(\bLambda^2)\Big\}\\
\ge~& \frac{1}{n^2(n+1)^2} \E [ \fnorm{\bLambda}^2 ]\\
=~& \frac{p-2}{n^2(n+1)^2} \E[\norm{\by^{(t)}}^2\norm{\by^{(t')}}^2].
\end{align*}

For $t=t'$, using a similar argument we obtain

\begin{align*}
\text{Var}\big\{[\hbS_{\rm{(mm)}} ]_{t,t'}\big\}
\ge~  \frac{p-1}{n^2(n+1)^2} \E[\norm{\by^{(t)}}^2\norm{\by^{(t')}}^2].
\end{align*}

Summing over all $(t,t')\in [T]\times [T]$ yields

\begin{align*}
\E \big[\fnorm{\hbS_{\rm{(mm)}} - \bS}^2\big]
&\ge \frac{p-2}{n^2(n+1)^2}  \sum_{t,t'}\E[\norm{\by^{(t)}}^2\norm{\by^{(t')}}^2]\\
&= \frac{p-2}{n^2(n+1)^2} \E [\fnorm{\bY}^4]\\
&\ge \frac{p-2}{n^2(n+1)^2} (\E [\fnorm{\bY}^2])^2\\
&= \frac{p-2}{(n+1)^2} [\trace(\bS) + \norm{\bB^*}^2]^2.
\end{align*}

§.§ Proof of Theoremthm:covariance

Recall definition of $\hbS$ in <Ref>, and let $\bQ_1$, $\bQ_2$ be defined as in <Ref>. With $\bZ = \bX\bSigma^{-1/2}$, we obtain

\begin{align*}
&n^2\Big[\bQ_2 (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12} - n^{-1} (n\bI_T - \hbA) \bQ_1 \big(\fnorm{\bS^{\frac 12}}\big(\fnorm*{\bF}^2/n + \trace(\bS)\big)^{\frac 12} n^{-\frac 12}\big) \Big]\\
=~& \big( \bF^\top \bZ\bZ^\top \bF - \bF^\top\bF(p\bI_T - \hbA) + (n\bI_T -\hbA)\bH^\top\bZ^\top\bF\big) \\
& - \big[(n\bI_T - \hbA) (\bF^\top\bF + \bH^\top\bZ^\top\bF- \bS (n\bI_T - \hbA))\big]\\
=~& \big( \bF^\top \bZ\bZ^\top \bF - \bF^\top\bF (p\bI_T - \hbA)\big) - (n\bI_T - \hbA) (\bF^\top\bF  - \bS (n\bI_T - \hbA))\\
=~& \bF^\top \bZ\bZ^\top \bF + \bF^\top\bF\hbA + \hbA\bF^\top\bF -(n+p)\bF^\top\bF + (n\bI_T - \hbA) \bS (n\bI_T - \hbA)
\\
=~& (n\bI_T - \hbA) \bS (n\bI_T - \hbA) + \bF^\top\bF\hbA + \hbA\bF^\top\bF - \bF^\top((n+p) \bI_T - \bZ\bZ^\top )\bF\\
=~& (n\bI_T - \hbA) \bS (n\bI_T - \hbA) - (n\bI_T - \hbA) \hbS (n\bI_T - \hbA)\\
=~& (n\bI_T - \hbA) (\bS - \hbS) (n\bI_T - \hbA).
\end{align*}

Therefore, by triangle inequality and $\opnorm*{\bI_T - \hbA/n}\le 1$ in <Ref>,

\begin{align*}
&\fnorm[\big]{(\bI_T - \hbA/n) (\bS - \hbS) (\bI_T - \hbA/n)}\\
\le~& \fnorm*{\bQ_2}n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n) +\opnorm*{\bI_T - \hbA/n}\fnorm*{\bQ_1}n^{-\frac 12} \fnorm{\bS^{\frac 12}}\big(\fnorm*{\bF}^2/n + \trace(\bS)\big)^{\frac 12} \\
\le~& \fnorm*{\bQ_2}n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n) +\fnorm*{\bQ_1}n^{-\frac 12} \fnorm{\bS^{\frac 12}}\big(\fnorm*{\bF}^2/n + \trace(\bS)\big)^{\frac 12} \\
\le~& \fnorm*{\bQ_2}n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n) +\fnorm*{\bQ_1}n^{-\frac 12} \frac12 \big[\trace(\bS)+ \big(\fnorm*{\bF}^2/n + \trace(\bS)\big)\big] \\
\le~ & (\fnorm*{\bQ_2} + \fnorm*{\bQ_1}) n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n + \trace(\bS)).
\end{align*}

\fnorm[\big]{(\bI_T - \hbA/n) (\bS -\hbS)(\bI_T - \hbA/n)} \le \Theta_1 n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n + \trace(\bS)),

where $\Theta_1 = \fnorm*{\bQ_1} + \fnorm*{\bQ_2}$.
Note that we have $\E[ \fnorm*{\bQ_1}^2] \le 4$ from <Ref>.
By <Ref>, we have

(1) under <Ref>(i), $\E [\fnorm*{\bQ_2}^2] \le C(\tau')(T \wedge (1 + \frac pn))(1 + \frac pn)$.

\begin{align*}
\E [\Theta_1^2] \le 2\E[\fnorm*{\bQ_1}^2 +\fnorm*{\bQ_2}^2]
&\le 2 [4 + C(\tau')(T \wedge (1 + \frac pn))(1 + \frac pn)]\\
&\le C(\tau')(T \wedge (1 + \frac pn))(1 + \frac pn).
\end{align*}

(2) under <Ref>(ii), $\E [I(\Omega)\fnorm*{\bQ_2}^2] \le C(\gamma, c)$ with $\P(\Omega)
\to 1$. Thus,
$\E [I(\Omega)\Theta_1^2] \le 2\E[\fnorm*{\bQ_1}^2 +I(\Omega)\fnorm*{\bQ_2}^2]\le %C(\eta)\frac pn \le
C(\gamma, c)

§.§ Proof of Theoremthm:generalization error

From the definitions of $\bQ_1, \bQ_2, \bQ_3$ in <Ref>, we have

\begin{align*}
&(\bQ_2^\top + \bQ_3) n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n) + (\bI_T - \hbA/n) \bQ_1 n^{-\frac 12}\fnorm{\bS^{\frac 12}}\big(\fnorm*{\bF}^2/n + \trace(\bS)\big)^{\frac 12}\\
=~& \frac{1}{n} \bF^\top \bF - (\bI_T - \hbA/n) (\bH^\top\bH + \bS) (\bI_T - \hbA/n)\\
=~& (\bI_T - \hbA/n)\big[n^{-1}(\bI_T - \hbA/n)^{-1} \bF^\top \bF (\bI_T - \hbA/n)^{-1} - (\bH^\top\bH + \bS)\big](\bI_T - \hbA/n)\\
=~& (\bI_T - \hbA/n) (\hbR - \bR) (\bI_T - \hbA/n),
\end{align*}

where $\hbR \defas n^{-1}(\bI_T - \hbA/n)^{-1}\bF^\top \bF (\bI_T - \hbA/n)^{-1}$, and $\bR \defas \bH^\top\bH + \bS$.

Therefore, by triangle inequality and $\opnorm*{\bI_T - \hbA/n}\le 1$ from <Ref>,

\begin{align*}
&\fnorm[\big]{(\bI_T - \hbA/n) (\hbR - \bR) (\bI_T - \hbA/n)} \\
\le~& (\fnorm*{\bQ_2} + \fnorm*{\bQ_3}) n^{-\frac 12} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n) + \fnorm*{\bQ_1}n^{-\frac 12} \fnorm{\bS^{\frac 12}}\big(\fnorm*{\bF}^2/n + \trace(\bS)\big)^{\frac 12} \\
\le~& (\fnorm*{\bQ_2} + \fnorm*{\bQ_3} + \fnorm*{\bQ_1}) n^{-\frac 12} (\fnorm*{\bF}^2/n +\fnorm{\bH}^2 + \trace(\bS))\\
=~& \Theta_2 n^{-\frac 12} (\fnorm*{\bF}^2/n +\fnorm{\bH}^2 + \trace(\bS)),
\end{align*}

where $\Theta_2 = \fnorm*{\bQ_1} + \fnorm*{\bQ_2} + \fnorm*{\bQ_3}$.
By <Ref>, we obtain
$\E [\Theta_2] \le C(\gamma, \tau')$ under <Ref>(i)
$\E [I(\Omega)\Theta_2] \le C(\gamma, c)$ with $P(\Omega)\to1$ under <Ref>(ii).

Furthermore, since $\Theta_2 = O_P(1)$, and $\opnorm*{(\bI_T - \hbA/n)^{-1}} = O_P(1)$ from <Ref>,

\begin{align*}
\fnorm{\hbR - \bR} &\le \opnorm*{(\bI_T - \hbA/n)^{-1}}^2 \Theta_2 n^{-\frac 12} (\fnorm*{\bF}^2/n +\fnorm{\bH}^2 + \trace(\bS))\\
&= O_P(n^{-\frac 12})  (\fnorm*{\bF}^2/n +\fnorm{\bH}^2 + \trace(\bS)).
\end{align*}

Since $\frac{1}{n} \bF^\top \bF = (\bI_T - \hbA/n) \hbR (\bI_T - \hbA/n)$, taking trace of both sides gives $\frac 1n \fnorm*{\bF}^2 \le \norm{\hbR}_*$ thanks to $\opnorm{(\bI_T - \hbA/n)} \le1$. Note that $\norm{\bR}_* = \fnorm{\bH}^2 + \trace(\bS)$ by definition of $\bR$, we obtain

\begin{align}\label{eq: R-R}
\fnorm{\hbR - \bR} &\le O_P(n^{-\frac 12})  (\norm{\hbR}_* + \norm{\bR}_*).
\end{align}

Since $\hbR$ and $\bR$ are both $T\times T$ positive semi-definite matrices, whose ranks are at most $T$,

\begin{align*}
&\big|\|\hbR\|_* - \|\bR\|_*\big| \le  \|\hbR - \bR\|_*
\le \sqrt{2T}\fnorm{\hbR - \bR}\\
\le~& O_P( (T/n)^{\frac 12}) (\norm{\hbR}_* + \norm{\bR}_*) = o_P(1) (\norm{\hbR}_* + \norm{\bR}_*),
\end{align*}

thanks to $T = o(n)$. That is,

\begin{align*}
\frac{\big|\|\hbR\|_* - \|\bR\|_*\big|}{\|\hbR\|_*+ \|\bR\|_*} \le O_P( (T/n)^{\frac 12}),
\end{align*}

which implies
$\frac{\|\bR\|_*}{\|\hbR\|_*} -1 = O_P( (T/n)^{\frac 12})$, ,

\[
\frac{\trace(\bS)+ \fnorm{\bH}^2}{\fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n} -1 = O_P( (T/n)^{\frac 12}) = o_P(1).
\]

§.§ Proof of Theoremthm:main

This proof is based on results of <Ref>.
We begin with the result of  <Ref>,

\begin{equation*}
\frac{\trace(\bS)+ \fnorm{\bH}^2}{\fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n} \overset{p}\to 1.
\end{equation*}

In other words,

\[
\trace(\bS)+ \fnorm{\bH}^2 = (1 + o_P(1)) \fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n.
\]

Thus, the upper bound in <Ref> can be bounded from above as follows

\begin{align*}
&\fnorm{(\bI_T - \hbA/n) (\hbS -\bS)(\bI_T - \hbA/n)} \\
\le~&\Theta_1 n^{-\frac 12} ( \fnorm*{\bF}^2/n + \fnorm{\bH}^2 + \fnorm{\bS^{\frac 12}}^2) \\
\le~& \Theta_1 n^{-\frac 12} (\fnorm*{\bF}^2/n + (1 + o_P(1)) \fnorm{(\bI_T - \hbA/n)^{-1}\bF^\top}^2/n)\\
\le~& \Theta_1 n^{-\frac 12} \big(1 + (1 + o_P(1)) \opnorm{(\bI_T -\hbA)^{-1}}^2\big)\fnorm*{\bF}^2/n\\
=~& O_P(n^{-\frac12}) \fnorm*{\bF}^2/n,
\end{align*}

Using $\opnorm{(\bI_T -\hbA)^{-1}}= O_P(1)$ again, it follows

\begin{align}\label{eq: pf1}
\fnorm{\hbS -\bS} \le O_P(n^{-\frac12})  \fnorm{\bF}^2/n.
\end{align}

A similar argument leads to

\begin{align}\label{eq: pf2}
\fnorm{\hbS -\bS} \le O_P(n^{-\frac12}) (\trace(\bS) + \fnorm{\bH}^2).
\end{align}

§.§ Proof of Corollarycor36

Under <Ref>(i)  and <ref>, we proceed to bound $\fnorm{\bF}^2/n$ in terms of $\trace(\bS)$.
Let $L(\bB)  = \frac{1}{2n}\fnorm*{\bY - \bX\bB }^2 + \lambda \norm{\bB}_{2,1} + \frac{\tau}{2} \fnorm{\bB}^2$ be the objective function in (<ref>), then $L(\hbB) \le L(\bf0)$ by definition of $\hbB$. Thus,

\begin{align*}
&\frac{1}{2n}\fnorm{\bF }^2 \le \frac{1}{2n}\fnorm{\bF }^2 + \lambda \norm{\hbB}_{2,1}
+ \frac{\tau}{2} \fnorm{\hbB}^2 \le\frac{1}{2n}\fnorm{\bY }^2.
\end{align*}

Now we bound $\frac1n\fnorm{\bY}^2$ by Hanson-Wright inequality.
Since $\bY = \bX\bB^* + \bE$, the rows of $\bY$ are $\calN_T(\bf0, \bSigma_{\by})$ with $\bSigma_{\by} = (\bB^*)^\top \bSigma \bB^*+ \bS$, then $\vec(\bY^\top) \sim \calN(\bf0, \bI_n \otimes \bSigma_{\by})$, and $\bxi \defas [\bI_n \otimes \bSigma_{\by}]^{-\frac12} \vec(\bY^\top)\sim \calN(\mathbold{0}, \bI_{nT})$.
$\fnorm{\bY}^2 = [\vec(\bY^\top)]^\top \vec(\bY^\top) = \bxi^\top (\bI_n \otimes \bSigma_{\by}) \bxi$,
we apply the following variant of Hanson-Wright inequality.

For $\bxi\sim \calN(\mathbold{0}, \bI_N)$, then

\begin{align*}
\P(\bxi^\top \bA\bxi - \trace(\bA) \le 2 \sqrt{x}\fnorm{\bA} + 2x\opnorm{\bA}) \ge1 - \exp(-x).
\end{align*}

In our case, take $\bA = (\bI_n \otimes \bSigma_{\by})$, then $\trace(\bA) = n\trace(\bSigma_{\by})$, $\fnorm{\bA}  =  \sqrt{n}\fnorm{\bSigma_{\by}}\le \sqrt{n} \trace(\bSigma_{\by})$, $\opnorm{\bA} = \opnorm{\bSigma_{\by}}\le \trace(\bSigma_{\by})$,
thus with probability at least $1 - \exp(-x)$,

\begin{align*}
\fnorm{\bY}^2 - n\trace(\bSigma_{\by})
\le 2 \sqrt{nx}\trace(\bSigma_{\by}) + 2x \trace(\bSigma_{\by}).
\end{align*}

Take $x=n$, then with probability at least $1 - \exp(-n)$,

\begin{align*}
\fnorm{\bF}^2/n \le  \fnorm{\bY}^2/n \le  5\trace(\bSigma_{\by}).
\end{align*}

Thus, $\fnorm{\bF}^2/n = O_P(1)\trace(\bSigma_{\by}).$ Together with (<ref>), we obtain

\begin{equation*}
\fnorm{\hbS -\bS} \le O_P(n^{-\frac12}) \trace(\bSigma_{\by}).
\end{equation*}

Note that by <Ref>, $\trace(\bSigma_{\by}) = \fnorm{\bSigma^{\frac12} \bB^*}^2 + \trace(\bS) \le (1 + \mathfrak{snr})\trace(\bS).$
Therefore, we obtain

\begin{equation*}
\fnorm{\hbS -\bS} \le O_P(n^{-\frac12}) \trace(\bS).
\end{equation*}

Furthermore, since $\trace(\bS)\le \sqrt{T}\fnorm{\bS}$ and $T = o(n)$, we have

\fnorm{\hbS -\bS} \le  O_P(\sqrt{ T/n}) \fnorm{\bS} = o_P(1) \fnorm{\bS}.

Finally, since $\norm{\bS}_* = \trace(\bS)$, by triangular inequality

\begin{equation*}
\big|\norm{\hbS}_* - \trace(\bS)\big|
\le \norm{\hbS - \bS}_*
\le \sqrt{T}\fnorm{\hbS - \bS}
\le O_P(\sqrt{T/n}) \trace(\bS)
= o_P(1) \trace(\bS).
\end{equation*}

§.§ Proof of Corollarycor37

For $\tau=0$, by the optimality of $\hbB$ in (<ref>),

\frac{1}{2n}\fnorm{\bF}^2 + \lambda\|\hbB\|_{2,1}
\le
\frac{1}{2n}\fnorm{\bE}^2  + \lambda \|\bB^*\|_{2,1}.

Note that $\bF = \bE - \bX(\hbB - \bB^*)= \bE - \bZ\bH$, expanding the squares and rearranging terms yields

\begin{equation}\label{eq:boundH1}
\fnorm{\bZ\bH}^2 \le 2\langle \bE, \bZ\bH\rangle + 2n\lambda (\|\bB^*\|_{2,1} - \|\hbB\|_{2,1})  \le 2\langle \bE, \bZ\bH\rangle + 2n\lambda \|\hbB - \bB^*\|_{2,1}.
\end{equation}

From assumptions in this corollary, $\hbB - \bB^*$ has at most $(1-c)n$ rows. Thus, in the event $U_2$, we have

n\eta \fnorm{\bH}^2 = n\eta \fnorm{\bSigma^{1/2}(\hbB - \bB^*)}^2\le \fnorm{\bX(\hbB - \bB^*)}^2 = \fnorm{\bZ\bH}^2.

We bound the right-hand side two terms in (<ref>) by Cauchy-Schwarz inequality,

\quad \|\hbB - \bB^*\|_{2,1}\le \sqrt{(1-c)n} \fnorm{\hbB - \bB^*}\le  \frac{\sqrt{(1-c)n}}{\sqrt{\phi_{\min}(\bSigma)}} \fnorm{\bH} \le \frac{\sqrt{1-c}}{\sqrt{\eta \phi_{\min}(\bSigma)}} \fnorm{\bZ\bH},

and $	\langle \bE, \bZ\bH\rangle \le \fnorm{\bE}  \fnorm{\bZ\bH} \le \fnorm{\bS^{\frac 12}} \opnorm{\bE\bS^{-\frac 12}} \fnorm{\bZ\bH}.$

Therefore, by canceling a factor $\fnorm{\bZ\bH}$ from both sides of (<ref>), we have

\begin{align*}
\sqrt{n\eta} \fnorm{\bH} \le \fnorm{\bZ\bH} \le 2 \fnorm{\bS^{\frac 12}} \opnorm{\bE\bS^{-\frac 12}}  + \frac{2\sqrt{(1-c)}n\lambda}{\sqrt{\eta \phi_{\min}(\bSigma)}}.
\end{align*}

Using $(a + b)^2 \le 2a^2 + 2b^2$,

\begin{align*}
\fnorm{\bH}^2 \le \frac{4}{n\eta} \trace(\bS) \opnorm{\bE\bS^{-\frac 12}}^2 + \frac{4(1-c)n\lambda^2}{\eta^2 \phi_{\min}(\bSigma)}.% \le \frac{36}{\eta} \trace(\bS) + \frac{4(1-c)n\lambda^2}{\eta^2 \phi_{\min}(\bSigma)}.
\end{align*}

Hence, using $\lambda$ is of the form $\mu\sqrt{\trace(\bS)/n}$, we have

\begin{align}
&\trace(\bS) + \fnorm{\bH}^2 \\
\le~& (1 +4 \eta^{-1} n^{-1}\opnorm{\bE\bS^{-\frac 12}}^2) \trace(\bS) + \frac{4(1-c)\mu^2}{\eta^2 \phi_{\min}(\bSigma)}\trace(\bS)\\
\le~& O_P(1) (1 + \mu^2)\trace(\bS),
\end{align}

where we used that $n^{-1}\opnorm{\bE\bS^{-\frac 12}} = O_P(1)$ by <cit.> and $T = o(n)$.
Now, by <Ref>,

\begin{align*}
\fnorm{\hbS - \bS} \le O_P(n^{-\frac12}) [\trace(\bS) + \fnorm{\bH}^2]\le O_P(n^{-\frac12}) (1 + \mu^2)\trace(\bS),
\end{align*}

where the $O_P(\cdot) $ hides constants depending on $\gamma, c, \phi_{\min}(\bSigma)$ since $\eta$ is a constant that only depends on $\gamma, c$.

§.§ Proof of Theoremthm:out-of-sample

From the definitions of $\bQ_2, \bQ_3$ in <Ref>, we have

\begin{align*}
&\bQ_2 + \bQ_2^\top + \bQ_3\\
=~& \frac{n^{-2}\big( \bF^\top \bZ\bZ^\top \bF +
\hbA \bF^\top \bF + \bF^\top \bF \hbA -
p\bF^\top \bF -  (n\bI_T -\hbA)\bH^\top\bH (n\bI_T -\hbA)}{(\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12}}.
\end{align*}

\begin{align*}
(\bI_T -\hbA/n)\bH^\top\bH (\bI_T -\hbA/n) -  n^{-2} \big( \bF^\top \bZ\bZ^\top \bF + \hbA \bF^\top \bF + \bF^\top \bF \hbA -
p\bF^\top \bF\big)}\\
=~& \fnorm{\bQ_2 +\bQ_2^\top + \bQ_3} (\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12}\\
\le ~&\Theta_3(\fnorm{\bH}^2 + \fnorm*{\bF}^2/n)n^{-\frac 12},
\end{align*}

where $\Theta_3= 2\fnorm*{\bQ_2} + \fnorm*{\bQ_3}$. The conclusion thus follows by <Ref>.

§ PROOFS OF PRELIMINARY RESULTS

§.§ Proofs of results in Appendixsec:opnorm-bound

(i) For any $\bu\in \R^T$, by defintion (<ref>),

\begin{align*}
\bu^\top\hbA\bu &= \trace\big[ (\bu^{\top} \otimes \bX_{\hat{\mathscr{S}}})
\bM^{\dagger}
(\bu \otimes \bX^\top_{\hat{\mathscr{S}}})\big]\\
&\le \trace\big[ (\bu^{\top} \otimes \bX_{\hat{\mathscr{S}}})
[\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})^{\dagger}]
(\bu \otimes \bX_{\hat{\mathscr{S}}}^\top)\big]\\
&= \trace\big[
(\bu^{\top}\bI_T \bu) \otimes [\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})^{\dagger}\bX_{\hat{\mathscr{S}}}^\top]\big]\\
&= \norm*{\bu}^2\trace[\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})^{\dagger}\bX_{\hat{\mathscr{S}}}^\top]\\
&= \norm*{\bu}^2\trace[\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})^{\dagger}].
\end{align*}

Let $r = \rank(\bX_{\hat{\mathscr{S}}})\le \min(n, |\hat{\mathscr{S}}|)$ be the rank of $\bX_{\hat{\mathscr{S}}}$, and $\hphi_1\ge\cdots\ge \hphi_{r}>0$ be the nonzero eigenvalues of $\frac 1n \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}$. We have

\begin{align*}
\opnorm{\hbA/n}&\le \frac 1n \trace[\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})^{\dagger}]\\
&= \frac 1n \trace[\frac 1n \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}(\frac 1n \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau\bP_{\hat{\mathscr{S}}})^{\dagger}]\\
&\le \frac rn \opnorm[\Big]{\frac 1n \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}(\frac 1n \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau\bP_{\hat{\mathscr{S}}})^{\dagger}}\\
&\le \frac{\hphi_1}{\hphi_1 + \tau}\le 1.
\end{align*}

Thus, $\opnorm{\bI - \hbA/n}\le1$ as $\hbA$ is positive semi-definite.

(ii) Note that
\opnorm{(\bI_T - \hbA/n)^{-1}} = (1- \opnorm*{\hbA/n})^{-1} \le 1 +  \frac{\hphi_1}{\tau},
\hphi_1
= \opnorm{\frac{1}{n} \bX_{\hat{\mathscr{S}}}^\top \bX_{\hat{\mathscr{S}}}}
\le  \opnorm{\frac{1}{n} \bX^\top \bX}\le \frac{1}{n}\opnorm{ \bX^\top\bSigma^{-\frac 12}}^2\opnorm{\bSigma}.

(1) in the event $\{\opnorm{\bX\bSigma^{-\frac12}} <  2\sqrt{n}+\sqrt{p}\}$,  we have

\opnorm{(\bI_T - \hbA/n)^{-1}} \le 1 + \tau^{-1} (2 + \sqrt{p/n})^2\opnorm{\bSigma} = 1 + (\tau')^{-1} (2 + \sqrt{p/n})^2.

$\E[\hphi_1] \le  \E[n^{-1}\opnorm{ \bX^\top\bSigma^{-\frac 12}}^2\opnorm{\bSigma}] \le [(1 + \sqrt{p/n})^2 +n^{-1}]\opnorm{\bSigma}$  by (<ref>).

\E \opnorm{(\bI_T - \hbA/n)^{-1}}\le 1 + \tau^{-1} \E[\hphi_1]  \le 1 + (\tau')^{-1} [(1 + \sqrt{p/n})^2 +n^{-1}].

(i) For $\tau=0$, using the same arguement as proof of <Ref>, we obtain

\begin{align*}
\bu^\top\hbA\bu
\le
\norm{\bu}^2\trace[\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})^{\dagger}]
\le \norm{\bu}^2 |\hat{\mathscr{S}}|.
\end{align*}

Thus, in the event $U_1$, we have $\opnorm{\hbA}/n \le |\hat{\mathscr{S}}|/n\le (1-c)/2<1$, hence

\opnorm{\bI_T - \hbA/n} \le 1.

(ii) In the event $U_1$, we have
$\opnorm{(\bI_T - \hbA/n)^{-1}} = (1- \opnorm*{\hbA/n})^{-1}\le (1- (1-c)/2)^{-1}$. Furthermore, $\E [ I(U_1) \opnorm{(\bI_T - \hbA/n)^{-1}}]  \le (1- (1-c)/2)^{-1}.$

Since $\bM^\dagger\preceq \bM_1^\dagger = \bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger$,

\begin{align*}
\opnorm{\bN} &= \opnorm{(\bI_T \otimes \bX)\bM^\dagger (\bI_T \otimes \bX^\top)}\\
&\le \opnorm{(\bI_T \otimes \bX)(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger) (\bI_T \otimes \bX^\top)}\\
&=\opnorm{\bX_{\hat{\mathscr{S}}}(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger \bX_{\hat{\mathscr{S}}}^\top}\\
&\le 1,
\end{align*}

where the first inequality uses $\opnorm{\bA\bB\bA^\top} \le \opnorm{\bA\bC\bA^\top}$ for $0\preceq \bB \preceq \bC$.

§.§ Proofs of results in Appendixsec:Lipschitz-fix-E

Fixing $\bE$, if $\bX,\bar\bX$ are two design matrices, and $\hbB, \bar\bB$ are the two corresponding multi-task elastic net estimates.
Let $\bZ = \bX \bSigma^{-\frac 12}$,
$\bar\bZ = \bar\bX \bSigma^{-\frac 12}$,
$\bar\bH = \bSigma^{\frac12}(\bar\bB - \bB^*)$,
$\bar\bF = \bY - \bar\bX\bar\bB$,
and $\bar D = [\fnorm{\bar\bH}^2 + \fnorm{\bar\bF}^2/n]^{\frac12}$.
Without loss of generality, we assume $\bar D \le D$.
Recall the multi-task elastic net estimate $\hbB = \argmin_{\bB\in\R^{p\times T}}
\big( \frac{1}{2n}\fnorm{\bY - \bX\bB }^2 + g(\bB) \big)$, where  $g(\bB) = \lambda \norm{\bB}_{2,1}
+ \frac{\tau}{2} \fnorm{\bB}^2$.
Define $\varphi:\bB\mapsto \frac 1{2n}
\fnorm{\bE+\bX(\bB^*-\bB)}^2 + g(\bB)$,
$\psi:\bB\mapsto \frac 1{2n} \fnorm{\bX(\hbB-\bB)}^2$
and $\zeta:\bB\mapsto \varphi(\bB) - \psi(\bB)$.
When expanding the squares, it is clear that $\zeta$ is the sum of a linear function and a $\tau$-strong convex penalty, thus $\zeta$ is $\tau$-strongly convex of $\bB$. Additivity of subdifferentials yields $\partial \varphi (\hbB) = \partial \zeta(\hbB) + \partial \psi(\hbB) = \partial \zeta(\hbB)$. By optimality of $\hbB$ we have ${\mathbf 0}_{p\times T}\in \partial \varphi(\hbB)$, thus ${\mathbf 0}_{p\times T}\in \partial \zeta(\hbB)$. By strong convexity of $\zeta$,
\zeta(\bar\bB) - \zeta(\hbB)
\ge \langle \partial \zeta(\hbB), \bar\bB - \hbB \rangle + \frac{\tau}{2} \fnorm{\bar\bB - \hbB}^2
=  \frac{\tau}{2} \fnorm{\bar\bB - \hbB}^2,
which can further be rewritten as

\fnorm{\bX (\hbB - \bar\bB)}^2 + n\tau \fnorm{\hbB -\bar\bB}^2
\le \fnorm{\bE - \bX (\bar\bB - \bB^*)}^2 -  \fnorm{\bE - \bX (\hbB - \bB^*)}^2 + 2n(g(\bar\bB) - g(\hbB)),

$$\fnorm{\bZ (\bH - \bar\bH)}^2 + n\tau \fnorm{\bSigma^{-\frac12}(\bH - \bar\bH)}^2
\le \fnorm{\bE - \bZ \bar\bH}^2 -  \fnorm{\bE - \bZ \bH}^2 + 2n(g(\bar\bB) - g(\hbB)).$$

Summing the above inequality with its counterpart obtained by replacing $(\bX, \hbB, \bH)$ with $(\bar\bX, \bar\bB, \bar\bH)$, we have

\begin{align*}
&(LHS) \\
\defas~&\fnorm{\bZ (\bH - \bar\bH)}^2 + \fnorm{\bar\bZ (\bH - \bar\bH)}^2  + 2n\tau' \fnorm{\bH - \bar\bH}^2\\
\le~& \fnorm{\bZ (\bH - \bar\bH)}^2 + \fnorm{\bar\bZ (\bH - \bar\bH)}^2  + 2n\tau \fnorm{\bSigma^{-\frac12}(\bH - \bar\bH)}^2\\
\le~& \fnorm{\bE - \bZ \bar\bH}^2 -  \fnorm{\bE - \bZ \bH}^2  + \fnorm{\bE - \bar\bZ \bH}^2 -  \fnorm{\bE - \bar\bZ \bar\bH}^2\\
=~& \langle \bZ(\bH - \bar\bH),\bF+\bar\bF +(\bar\bZ- \bZ)\bar\bH\rangle
+ \langle -\bar\bZ(\bH - \bar\bH),\bF+\bar\bF +(\bZ-\bar\bZ)\bH\rangle\\
=~& \langle (\bZ-\bar\bZ)(\bH - \bar\bH),\bF+\bar\bF\rangle + \langle \bZ(\bH - \bar\bH),(\bar\bZ- \bZ)\bar\bH\rangle + \langle \bar\bZ(\bar\bH - \bH), (\bZ-\bar\bZ)\bH\rangle\\
\le~& \opnorm{\bZ-\bar\bZ}\fnorm{\bH - \bar\bH} (\fnorm{\bF} + \fnorm{\bar\bF})
+ \opnorm{\bZ-\bar\bZ}\fnorm{\bZ (\bH - \bar\bH)}   \fnorm{\bar\bH}\\
&+ \opnorm{\bZ-\bar\bZ}\fnorm{\bar\bZ(\bH - \bar\bH)} \fnorm{\bH}\\
\le~&  \opnorm{\bZ-\bar\bZ}
\Big[
\sqrt{\frac{(LHS)}{2n\tau'}} (\fnorm{\bF} + \fnorm{\bar\bF}) + \sqrt{(LHS)} (\fnorm{\bar\bH} + \fnorm{\bH})
\Big]\\
\le~& \opnorm{\bZ-\bar\bZ}  \sqrt{(LHS)} (D + \bar D)\max(1, (2\tau')^{-\frac12})
\end{align*}

where $\tau' = \tau \phi_{\min}(\bSigma^{-1}) = \tau/\opnorm{\bSigma}$.
That is,

\begin{align*}
\sqrt{(LHS)}
\le
\opnorm{\bZ-\bar\bZ} 2D \max(1, (2\tau')^{-\frac12}).
\end{align*}

\begin{align*}
&n^{-\frac12}\fnorm{\bF - \bar\bF} = n^{-\frac12}\fnorm{\bZ\bH - \bar\bZ\bar\bH}\\
\le~& n^{-\frac12}[\fnorm{\bZ(\bH-\bar\bH)} + \fnorm{(\bZ- \bar\bZ)\bar\bH}]\\
\le~&n^{-\frac12}[\fnorm{\bZ(\bH-\bar\bH)} + \opnorm{\bZ- \bar\bZ}\fnorm{\bH}]\\
\le~&n^{-\frac12}[\sqrt{(LHS)}+ \opnorm{\bZ- \bar\bZ}D]\\
\le~&n^{-\frac12}  \opnorm{\bZ- \bar\bZ} D[2 \max(1, (2\tau')^{-\frac12}) + 1].
\end{align*}

So far we obtained

\begin{align*}
\fnorm{\bH - \bar\bH}
\sqrt{\frac{(LHS)}{2n\tau'}}\le
n^{-\frac12} \opnorm{\bZ-\bar\bZ}D (2\tau')^{-\frac12} 2\max(1, (2\tau')^{-\frac12}),\\
n^{-\frac12}\fnorm{\bF - \bar\bF}
&\le n^{-\frac12}  \opnorm{\bZ- \bar\bZ} D[2\max(1, (2\tau')^{-\frac12}) + 1].
\end{align*}

Let $\bQ = [\bH^\top, \bF^\top/\sqrt{n}]^\top$ and $\bar \bQ = [\bar\bH^\top, \bar\bF^\top/\sqrt{n}]^\top$, then $D = \fnorm{\bQ}$, $\bar D = \fnorm{\bar \bQ}$. By triangular inequality,

\begin{align*}
|D - \bar D|\le \fnorm{\bQ-\bar \bQ} \le~& \fnorm{\bH-\bar\bH} + \fnorm{\bF-\bar\bF}/\sqrt{n} \\
\le~ & n^{-\frac12}\opnorm{\bZ- \bar\bZ} D [4 \max(1, (2\tau')^{-1})],
\end{align*}

where the last inequality uses the elementary inequality $\max(a,b)(a+b)\le 2 [\max(a,b)]^2$ for $a,b>0$ with $a = 1, b = (2\tau')^{-\frac12}$.
Let $\frac{\partial D}{\partial \bZ}\defas \frac{\partial D}{\partial \vec(\bZ)} \in \R^{1\times np}$, then $\norm*{\frac{\partial D}{\partial \bZ}} \le n^{-\frac12}D L_1$ with $L_1 = [4 \max(1, (2\tau')^{-1})]$. Hence,

\begin{align*}
\sum_{ij} \Big(\frac{\partial D}{\partial z_{ij}}\Big)^2
=  \norm*{\frac{\partial D}{\partial \bZ}}^2
\le n^{-1}D^2 L_1^2.
\end{align*}

Furthermore, by triangle inequality

\begin{align*}
\fnorm[\Big]{\frac{\bQ}{D} - \frac{\bar \bQ}{\bar D}}
& \le \frac1D \fnorm{\bQ-\bar \bQ} + \Big|\frac1D - \frac{1}{\bar D} \Big| \fnorm{\bar \bQ}\\
& = \frac1D \fnorm{\bQ-\bar \bQ} + \frac{|D-\bar D|}{D\bar D} \fnorm{\bar \bQ}\\
&\le \frac1D \fnorm{\bQ-\bar \bQ} + \frac{1}{D} \fnorm{\bQ - \bar \bQ}\\
&\le n^{-\frac12}\opnorm{\bZ- \bar\bZ} L,
\end{align*}

where $L = 8 \max(1, (2\tau')^{-1})$.
Therefore, when $\tau >0$, we obtain the two mappings $\bZ \mapsto D^{-1}\bF/\sqrt{n}$, and $\bZ \mapsto D^{-1}\bH$ are both $n^{-\frac12} L$-lipschitz with $L = 8 \max(1, (2\tau')^{-1})$, where $\tau' = \tau/\opnorm{\bSigma}$.

The proof of <Ref> uses a similar argument as proof of <Ref>, we present it here for completeness.

For multi-task group Lasso ($\tau=0$), we restrict our analysis in the event $U_1\cap U_2$, where $U_1 = \big\{ \norm{\hbB}_0 \le n(1-c)/2 \big\}$, $U_2 = \big\{\inf_{\bb\in \R^p: \| \bb\|_0 \le (1-c)n}  \|\bX \bb\|^2/(n \|\bSigma^{\frac 12} \bb\|^2)  > \eta\big\}.$

Since the only randomness of the problem comes from $\bX$ and $\bE$, there exists a measurable set $\calU$ such that $U_1\cap U_2 =\{ (\bX, \bE)\in \calU\}$.
For some noise matrix $\bE$, consider  $\bX,\bar\bX$ two design matrices such that $(\bX, \bE)\in \calU$ and $(\bar\bX, \bE)\in \calU$.
We slightly abuse the notation and let $\hbB, \bar\bB$ denote the two corresponding multi-task group-Lasso estimates.
Thus, the row sparsity of $ \hbB-\bar\bB$ is at most $n(1-c)$.
$\bar\bH = \bar\bB - \bB^*$, $\bar\bF = \bY - \bar\bX\bar\bB$, and $\bar D = [\fnorm{\bar\bH}^2 + \fnorm{\bar\bF}^2/n]^{\frac12}$.
Without loss of generality, we assume $\bar D \le D$.
Since when $\tau=0$, the multi-task group Lasso estimate is $\hbB = \argmin_{\bB\in\R^{p\times T}}
\big( \frac{1}{2n}\fnorm{\bY - \bX\bB }^2 + g(\bB) \big)$, where  $g(\bB) = \lambda \norm{\bB}_{2,1}$.
Define $\varphi:\bB\mapsto \frac 1{2n}
\fnorm{\bE+\bX(\bB^*-\bB)}^2 + g(\bB)$,
$\psi:\bB\mapsto \frac 1{2n} \fnorm{\bX(\hbB-\bB)}^2$
and $\zeta:\bB\mapsto \varphi(\bB) - \psi(\bB)$.
Under $\tau=0$, by the same arguments in proof of <ref> with the same functions $\varphi(\cdot), \psi(\cdot), \zeta(\cdot)$, we obtain

\fnorm{\bX (\hbB - \bar\bB)}^2
\le \fnorm{\bE - \bX (\bar\bB - \bB^*)}^2 -  \fnorm{\bE - \bX (\hbB - \bB^*)}^2 + 2n(g(\bar\bB) - g(\hbB)).

Summing the above inequality with its counterpart obtained by replacing $(\bX, \hbB, \bH)$ with $(\bar\bX, \bar\bB, \bar\bH)$, we have

\begin{align*}
&\fnorm{\bX (\hbB - \bar\bB)}^2 + \fnorm{\bar\bX (\hbB - \bar\bB)}^2\\
\le~& \fnorm{\bE - \bZ \bar\bH}^2 - \fnorm{\bE - \bZ \bH}^2
+ \fnorm{\bE - \bar\bZ \bH}^2 - \fnorm{\bE - \bar\bZ \bar\bH}^2.
\end{align*}

Note that in event $U_1\cap U_2$, we have

\begin{align*}
\eta n\fnorm{\bSigma^{\frac12}(\hbB - \bar\bB)}^2
\le \fnorm{\bX (\hbB - \bar\bB)}^2, \quad \eta n\fnorm{\bSigma^{\frac12}(\hbB - \bar\bB)}^2 \le \fnorm{\bar\bX (\hbB - \bar\bB)}^2.
\end{align*}

2\eta n\fnorm{(\hbH - \bar\bH)}^2 \le \fnorm{\bZ(\bH - \bar\bH)}^2 + \fnorm{\bar\bZ(\bH - \bar\bH)}^2
$, and

\begin{align*}
\defas~& \max (2\eta n\fnorm{\bH - \bar\bH}^2,\fnorm{\bZ(\bH - \bar\bH)}^2 + \fnorm{\bar\bZ(\bH - \bar\bH)}^2)\\
%	&2\eta n\fnorm{\bSigma^{\frac 12}(\bH - \bar\bH)}^2\\
=~&\fnorm{\bZ(\bH - \bar\bH)}^2 + \fnorm{\bar\bZ(\bH - \bar\bH)}^2\\
\le~& \fnorm{\bE - \bZ \bar\bH}^2 - \fnorm{\bE - \bZ \bH}^2
+ \fnorm{\bE - \bar\bZ \bH}^2 - \fnorm{\bE - \bar\bZ \bar\bH}^2.
\end{align*}

Now, in $U_1\cap U_2$, the Lipschitz property of the map $\bZ \mapsto D^{-1}\bF/\sqrt{n}$ follows from the same arguments in proof of <Ref>, with $\tau'$ in <ref> replaced by $\eta$ in this proof.

Furthermore, in the event $U_1\cap U_2\cap U_3$, the Lipschitz property of $\bZ \mapsto D^{-1}\bZ^\top\bF/n$ follows by triangle inequality. To see this, let $\bU = D^{-1}\bF/\sqrt{n}$, and $\bV = D^{-1}\bZ^\top\bF/n = n^{-1/2} \bZ^\top \bU$, thus by  triangle inequality

\begin{align*}
\opnorm{\bV - \bar\bV} &= n^{-1/2} \opnorm{\bZ^\top \bU - \bar\bZ^\top \bar\bU}\\
&= n^{-1/2}[ \opnorm{(\bZ - \bar\bZ)^\top \bU} + \opnorm{\bar\bZ^\top (\bU - \bar\bU)}]\\
&\le n^{-1/2}[ \opnorm{\bZ - \bar\bZ} + \opnorm{\bar\bZ}\opnorm{\bU - \bar\bU}]\\
&\le n^{-1/2}( 1 + n^{-1/2}\opnorm{\bar\bZ}L)
\opnorm{\bZ - \bar\bZ}\\
&\le n^{-1/2} (1 + (2 +\sqrt{p/n})L).
\end{align*}

where the last line uses $\opnorm{\bar\bZ}\le 2\sqrt{n} +\sqrt{p} $ in the event $U_3$.

<Ref> (1) is a direct consequence of the intermediate result $|D - \bar D| \le  n^{-\frac12}\opnorm{\bZ- \bar\bZ} D [4 \max(1, (2\tau')^{-1})]$ in proof of <Ref>, while <Ref> (2) is a direct consequence of the intermediate result $|D - \bar D| \le  n^{-\frac12}\opnorm{\bZ- \bar\bZ} D [4 \max(1, (2\eta)^{-1})]$ in proof of <Ref>.

Before proving the derivative formula, we restate $\hbB$ (defined in (<ref>) of the full paper) below,

\begin{equation}\label{eq: hbB-1}
\hbB=\argmin_{\bB\in\R^{p\times T}}
\Big(
\frac{1}{2n}\fnorm*{\bY - \bX\bB }^2 + \lambda \norm{\bB}_{2,1}
+ \frac{\tau}{2} \fnorm{\bB}^2
\Big),
\end{equation}

where $\|\bB\|_{2,1} = \sum_{j=1}^p \|{\bB^{\top} \be_j}\|_2$.

For the reader's convenience, we recall some useful notations. $\bP_{\hat{\mathscr{S}}} = \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top$. For each $k\in \hat{\mathscr{S}}$, $\bH^{(k)}=\lambda\|\hbB{}^\top \be_k\|_2^{-1}\left(\bI_T - \hbB{}^\top\be_k \be_k^\top\hbB ~  \|\hbB{}^\top\be_k\|_2^{-2} \right)$. $\tbH = \sum_{k\in\hat{\mathscr{S}}} (\bH^{(k)} \otimes \be_k\be_k^\top).$ $\bM_1 = \bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})$, $\bM = \bM_1 + n\tbH\in \R^{pT\times pT}$,  and $\bN = (\bI_T \otimes \bX)\bM^\dagger (\bI_T \otimes \bX^\top)$.

We first derive $\frac{\partial F_{lt}}{\partial x_{ij}}$.
Since $\bF = \bY - \bX\hbB = \bE - \bX(\hbB -\bB^*)$, by product rule,

\begin{align*}
\frac{\partial F_{lt}}{\partial x_{ij} } = \be_l^\top \frac{\partial \bE - \bX(\hbB -\bB^*) }{\partial x_{ij} } \be_t
= - \be_l^\top (\dot\bX (\hbB -\bB^*) + \bX\dot\bB) \be_t,
\end{align*}

$\dot\bX \defas \frac{\partial \bX}{\partial x_{ij}} = \be_i\be_j^\top$,
and $\dot\bB \defas \frac{\partial \hbB}{\partial x_{ij}}$.

Now we derive $\vec(\dot\bB)$ from KKT conditions for $\hbB$ defined in (<ref>):

1) For $k\in \hat{\mathscr{S}}$, , $\hbB{}^\top \be_k \ne\mathbf{0}$,

$$\be_k^\top\bX^\top\big[\bE-\bX(\hbB-\bB^*)\big] -n\tau\be_k^\top \hbB= \frac{n \lambda}{\|\hbB{}^\top\be_k\|_2} \be_k^\top \hbB
\quad \in\R^{1\times T}.$$

2) For $k\notin \hat{\mathscr{S}}$, , $\hbB{}^\top \be_k = \mathbf 0$,

$$\norm*{\be_k^\top\bX^\top\big[\bE-\bX(\hbB-\bB^*)\big] -n\tau\be_k^\top \hbB}< n\lambda.$$

Here the strict inequality is guaranteed by Proposition 2.3  of [Bellec, 2020].

Keeping $\bE$ fixed, differentiation of the above display for $k\in\hat{\mathscr{S}}$ w.r.t. $x_{ij}$ yields

$$\be_k^\top\Big[\dot\bX{}^\top\bF - \bX^\top[ \dot\bX(\hbB-\bB^*) +\bX\dot\bB]-n\tau\dot \bB\Big]= n\be_k^\top \dot\bB \bH^{(k)},$$

with $\bH^{(k)}
\lambda
\|\hbB{}^\top \be_k\|_2^{-1}\left(\bI_T - \hbB{}^\top\be_k \be_k^\top\hbB ~  \|\hbB{}^\top\be_k\|_2^{-2} \right)\in\R^{T\times T}$.
and using $\dot\bX = \be_i\be_j^\top$,

\[
\be_k^\top\Big[\be_j\be_i^\top\bF - \bX^\top \be_i \be_j^\top(\hbB-\bB^*)\Big]= \be_k^\top [(\bX^\top\bX + n\tau\bI_{p}) \dot\bB + n\dot\bB \bH^{(k)}].
\]

Recall $\bP_{\hat{\mathscr{S}}} = \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top$. Multiplying by $\be_k$ to the left and summing over $k\in \hat{\mathscr{S}}$, we obtain

\[
\bP_{\hat{\mathscr{S}}}\Big[\be_j\be_i^\top\bF - \bX^\top \be_i \be_j^\top(\hbB-\bB^*)\Big]= \bP_{\hat{\mathscr{S}}} (\bX^\top\bX + n\tau\bI_{p}) \dot\bB + n \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top \dot\bB \bH^{(k)}.
\]

Since $\hat{\mathscr{S}}$ is locally constant in a small neighborhood of $\bX$, $\hbB_{\hat{\mathscr{S}}{}^c}=0$, $\supp(\dot\bB)\subseteq \hat{\mathscr{S}}$. Hence, $\bP_{\hat{\mathscr{S}}}\dot\bB = \dot\bB$, and $\bX\dot\bB = \bX_{\hat{\mathscr{S}}}\dot\bB$. The above display can be rewritten as

\[
\bP_{\hat{\mathscr{S}}} \be_j\be_i^\top\bF - \bX_{\hat{\mathscr{S}}}^\top \be_i \be_j^\top(\hbB-\bB^*)
(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}}) \dot\bB + n \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top \dot\bB \bH^{(k)}.
\]

Vectorizing the above display using property $\vec(\bA\bB\bC) = (\bC^\top \otimes \bA)\vec(\bA)$ yields

\begin{align*}
&(\bF^\top \otimes \bP_{\hat{\mathscr{S}}} \be_j) \vec(\be_i^\top) -
((\hbB-\bB^*)^\top\be_j\otimes \bX_{\hat{\mathscr{S}}}^\top)\vec(\be_i)
\\=~&
[\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}) +  n\sum_{k\in\hat{\mathscr{S}}} (\bH^{(k)} \otimes \be_k\be_k^\top)] \vec(\dot\bB)\\
=~& (\bM_1 + n \tbH)\vec(\dot\bB)\\
=~& \bM \vec(\dot\bB),
\end{align*}

where $\bM_1 = \bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})$, and $\tbH = \sum_{k\in\hat{\mathscr{S}}} (\bH^{(k)} \otimes \be_k\be_k^\top)$.

Under <Ref>(i) that $\tau>0$, it's obviously that $\rank(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})) = T |\hat{\mathscr{S}}|$.
Under <Ref>(ii) that $\tau=0$ with $\P(U_1)\to 1$. In the event $U_1\cap U_2$, we know $\rank(\bX_{\hat{\mathscr{S}}}) = |\hat{\mathscr{S}}|$ from <cit.>, hence $\rank(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})) = T |\hat{\mathscr{S}}|$.
In either of the above two scenarios, we thus have $\dim(\ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}))) = T (p - |\hat{\mathscr{S}}|)$ by rank-nullity theorem. Since $[\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})] (\be_t\otimes\be_k) = \bf0$ for $t\in[T], k\in \hat{\mathscr{S}}^c$.
Let $V = \{(\be_t\otimes \be_k): t\in[T], k\in\hat{\mathscr{S}}^c\}$ be a vector space, then the elements of $V$ are linear independent, and $\dim(V) = T (p -|\hat{\mathscr{S}}|)$. Thus, $V$ forms a basis for $\ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})$. Since for any $(\be_t\otimes \be_k) \in V$,  we also have $\tbH (\be_t\otimes \be_k) = \bf0$, $\ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})) \subseteq \ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}) + n\tbH)$. On the other hand, if any vector $\bv$ s.t. $[\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}) + n\tbH]\bv = \bf0$, since these matrices are all positive semi-definite, we have $\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})\bv = \bf0$, which implies that $\ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})+ n\tbH) \subseteq \ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}))$. Therefore,

\begin{align*}
\ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})+ n\tbH) &= \ker(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}}))\\
&= \mathrm{span}\{(\be_t\otimes \be_k): t\in[T], k\in\hat{\mathscr{S}}^c\},
\end{align*}

\begin{align*}
\range(\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}+ n\tau\bP_{\hat{\mathscr{S}}})+ n\tbH)
&= \mathrm{span}\{(\be_t\otimes \be_k): t\in[T], k\in\hat{\mathscr{S}}\}.
\end{align*}

Since $\dot\bB = \bP_{\hat{\mathscr{S}}}\dot\bB$, $\vec(\dot\bB) = (\bI_T \otimes \bP_{\hat{\mathscr{S}}})\vec(\dot\bB)$, then $\vec(\dot\bB) \in \mathrm{col}(\bI_T\otimes \bP_{\hat{\mathscr{S}}}) = \range (\bM)$. Since $\bM$ is symmetric, $\bM^{\dagger} \bM$ is the orthogonal projection on the range of $\bM$.

\begin{align}\label{eq: vecB}
\vec(\dot\bB) = \bM^{\dagger} \bM \vec(\dot\bB) = \bM^{\dagger} [(\bF^\top\otimes \be_j) - ((\hbB-\bB^*)^\top\be_j\otimes\bX^\top)] \be_i.
\end{align}

Since $\supp(\dot\bB)\subseteq \hat{\mathscr{S}}$, $\bX\dot\bB = \bX_{\hat{\mathscr{S}}}\dot\bB$, we have

\begin{align*}
\frac{\partial F_{lt}}{\partial x_{ij} }
&= - \be_l^\top (\dot\bX(\hbB-\bB^*) + \bX\dot\bB) \be_t\\
&= - (\be_l^\top \be_i\be_j^\top(\hbB-\bB^*)\be_t +
\be_l^\top\bX_{\hat{\mathscr{S}}}\dot\bB\be_t)\\
&= -(\be_l^\top \be_i\be_j^\top(\hbB-\bB^*)\be_t + (\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\vec(\dot\bB))\\
&= - \be_l^\top \be_i\be_j^\top(\hbB-\bB^*)\be_t - (\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\bM^{\dagger} [(\bF^\top\otimes \be_j) - ((\hbB-\bB^*)^\top\be_j\otimes\bX^\top)] \be_i\\
&= - (e_j^\top(\hbB-\bB^*) \otimes \be_i^\top) (\be_t\otimes \be_l) +
(\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\bM^{\dagger}((\hbB-\bB^*)^\top\be_j\otimes\bX^\top\be_i) \\
&\quad -
(\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\bM^{\dagger} (\bF^\top\otimes \be_j)\be_i\\
&= - (e_j^\top(\hbB-\bB^*) \otimes \be_i^\top) (\be_t\otimes \be_l) +
(e_j^\top(\hbB-\bB^*) \otimes \be_i^\top)\bN (\be_t\otimes \be_l)\\
&\quad -
(\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\bM^{\dagger} (\bF^\top\otimes \bI_p) (\be_i\otimes \be_j)\\
&= - (e_j^\top(\hbB-\bB^*) \otimes \be_i^\top) (\bI_{nT} -\bN)(\be_t\otimes \be_l)
(\be_t^\top \otimes \be_l^\top\bX)\bM^{\dagger} (\bF^\top\otimes \bI_p) (\be_i\otimes \be_j)
\end{align*}

Now we calculate $\frac{\partial F_{lt}}{\partial z_{ij}}$.
Since $\bX = \bZ \bSigma^{\frac 12}$,
$x_{ik} = \sum_{j=1}^p z_{ij} (\bSigma^{\frac 12})_{jk}$, $ \frac{\partial x_{ik}}{\partial z_{ij}} = (\bSigma^{\frac 12})_{jk}$,

\begin{align*}
\frac{\partial F_{lt}}{\partial z_{ij}}
= \sum_{k=1}^p \frac{\partial F_{lt}}{\partial x_{ik}}  \frac{\partial x_{ik}}{\partial z_{ij}}
= \sum_{k=1}^p \frac{\partial F_{lt}}{\partial x_{ik}} (\bSigma^{\frac 12})_{jk}
= D_{ij}^{lt} + \Delta_{ij}^{lt},
\end{align*}

\begin{align*}
D_{ij}^{lt} %&= \sum_{k=1}^p D_{ik}^{lt}  (\bSigma^{\frac 12})_{jk}\\
&= -\sum_{k=1}^p (\be_k^\top(\hbB-\bB^*) \otimes \be_i^\top) (\bI_{nT} - \bN) (\be_t\otimes \be_l) (\bSigma^{\frac 12})_{jk}\\
&= -(\be_j^\top \bSigma^{\frac 12} (\hbB-\bB^*) \otimes \be_i^\top) (\bI_{nT} - \bN) (\be_t\otimes \be_l)\\
&= -(\be_j^\top\bH \otimes \be_i^\top) (\bI_{nT} - \bN) (\be_t\otimes \be_l),
\end{align*}

\begin{align*}
\Delta_{ij}^{lt}
&=-\sum_{k=1}^p  (\be_t^\top \otimes \be_l^\top)(\bI_T\otimes \bX )
\bM^\dagger\bigl(\bF^\top \otimes \bI_{p}\bigr)(\be_i \otimes\be_k) (\bSigma^{\frac 12})_{jk}\\
&=- (\be_t^\top \otimes \be_l^\top)(\bI_T\otimes \bX)
\bM^\dagger\bigl(\bF^\top \otimes \bI_{p}\bigr)(\be_i \otimes \bSigma^{\frac 12}\be_j)\\
&=- (\be_t^\top \otimes \be_l^\top)(\bI_T\otimes \bX)
\bM^\dagger (\bI_T\otimes \bSigma^{\frac 12}) \bigl(\bF^\top \otimes \bI_{p}\bigr)(\be_i \otimes \be_j)\\
%	&= -(\be_t^\top \otimes \be_l^\top)(\bI_T\otimes \bZ )
%	\tbM^\dagger\bigl(\bF^\top \otimes \bI_{p}\bigr)(\be_i \otimes\be_j).
\end{align*}

It follows that

\begin{align*}
\sum_{i=1}^n D_{ij}^{it}
&=-\sum_{i=1}^n  (\be_j^\top\bH \otimes \be_i^\top) (\bI_{nT} - \bN) (\be_t\otimes \be_i)\\
&=-  \be_j^\top\bH \big[ \sum_{i=1}^n  (\bI_T \otimes \be_i^\top) (\bI_{nT} - \bN) (\bI_T \otimes \be_i)\big] \be_t\\
&= -\be_j^\top \bH (n\bI_T - \hbA)\be_t,
\end{align*}

where the last line follows from definition of $\hbA$ in (<ref>).

(1) For $\tau>0$, by formula of $\frac{\partial F_{lt}}{\partial z_{ij} }$ in <Ref>, we have

\begin{align*}
&\sum_{ij}\norm*{\frac{\partial \bF}{\partial z_{ij}}}^2_{\rm F} = \sum_{ij}\sum_{lt} \Big(\frac{\partial F_{lt}}{\partial z_{ij} } \Big)^2=\sum_{ij}\sum_{lt} \Big( D_{ij}^{lt} + \Delta_{ij}^{lt} \Big)^2 \\
\le~ & 2\sum_{ij,lt} ( D_{ij}^{lt})^2 + 2\sum_{ij,lt} (\Delta_{ij}^{lt})^2 \\
=~ & 2 \fnorm*{(\bH\otimes \bI_n)(\bI_{nT} - \bN)}^2 + 2 \fnorm{(\bI_T\otimes \bX)\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2\\
\le~ & 2n\fnorm{\bH}^2 + 2 \fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2.
\end{align*}

Since $0\preceq \bM^\dagger \preceq  \bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger$,

\begin{align*}
&\fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2\\
\le~& \opnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})}^2
\fnorm{(\bF^\top \otimes \bI_{p})}^2\\
\le~& p\opnorm{\bSigma} \fnorm{\bF}^2
\opnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger}^2 \\
\le~ &p\opnorm{\bSigma}\fnorm{\bF}^2\opnorm{(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger\bX_{\hat{\mathscr{S}}}^\top}^2\\
\le~& \frac{p}{n\tau}\opnorm{\bSigma} \fnorm*{\bF}^2\\
=~& \frac{p}{n\tau'} \fnorm*{\bF}^2,
\end{align*}

where the last inequality uses $\opnorm{(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger\bX_{\hat{\mathscr{S}}}^\top}\le (n\tau)^{-1}$.

On the other hand, we also have

\begin{align*}
&\fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2\\
\le~& \fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger }^2
\opnorm{(\bI_T\otimes \bSigma^{\frac12})(\bF^\top \otimes \bI_{p})}^2\\
\le~& \fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}}) (\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger )}^2
\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T\fnorm{\bX_{\hat{\mathscr{S}}} (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger}^2
\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T\trace\big[
(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger
\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}
(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger\big]
\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T\trace\big[
(\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + \tau n \bP_{\hat{\mathscr{S}}})^\dagger\big]
\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T (\tau)^{-1}\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T\trace\big[
(\tau n \bP_{\hat{\mathscr{S}}})^\dagger\big]
\fnorm{\bF}^2 \opnorm{\bSigma}\\
\le~& T (\tau)^{-1}\fnorm{\bF}^2 \opnorm{\bSigma}\\
=~& \frac{T}{\tau'} \fnorm*{\bF}^2,
\end{align*}

\begin{align*}
\frac{1}{n}\sum_{ij}\norm*{\frac{\partial \bF}{\partial z_{ij}}}^2_{\rm F}
&\le 2\fnorm{\bH}^2 + 2
(\tau')^{-1} (T\wedge \frac{p}{n})\fnorm*{\bF}^2/n \\
&\le 2\max(1, (\tau')^{-1} (T\wedge \frac{p}{n})) (\fnorm*{\bF}^2/n + \fnorm{\bH}^2)\\
&= 2\max(1, (\tau')^{-1} (T\wedge \frac{p}{n})) D^2.
\end{align*}

Now by product rule and triangle inequality

\begin{align*}
&\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF/D}{\partial z_{ij}}}^2 \\
\le~&
2 D^{-2}\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF}{\partial z_{ij}}}^2
2 \frac{1}{n}\sum_{ij}\fnorm*{\bF\frac{\partial D^{-1}}{\partial z_{ij}}}^2\\
=~& 2 D^{-2}\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF}{\partial z_{ij}}}^2
2 D^{-4}\frac{1}{n}\fnorm{\bF}^2 \sum_{ij}\Big(\frac{\partial D}{\partial z_{ij}}\Big)^2\\
\le~& 2 D^{-2}\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF}{\partial z_{ij}}}^2
2 D^{-4}\frac{1}{n}\fnorm{\bF}^2 n^{-1} D^2 [4 \max(1, (2\tau')^{-1})]^2\\
\le~& 2 D^{-2}\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF}{\partial z_{ij}}}^2
2  n^{-1}  [4 \max(1, (2\tau')^{-1})]^2\\
\le~& 4 \max(1, (\tau')^{-1} (T\wedge \frac{p}{n})) + 2  n^{-1}  [4 \max(1, (2\tau')^{-1})]^2\\
:=~& f(\tau', T, n, p),
\end{align*}

where the second inequality is by <Ref>.

(2) For $\tau=0$, by <Ref>, in the event $U_1\cap U_2$, we obtain the same upper bounds as in the first case (1) with $\tau'$ replaced by $\eta$. To see this,

\begin{align*}
&\fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2\\
\le~& \opnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})}^2
\fnorm{(\bF^\top \otimes \bI_{p})}^2\\
=~&  \opnorm{(\bI_T\otimes \bSigma^{\frac12})\bM^\dagger(\bI_T\otimes \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})}
\le~&  \opnorm{(\bI_T\otimes \bSigma^{\frac12})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})}
\le~& p\ \fnorm{\bF}^2 \frac{1}{n\eta}\\
=~& \frac{p}{n\eta} \fnorm*{\bF}^2,
\end{align*}

where the third inequality is by <Ref>.
Also, we have

\begin{align*}
&\fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12}) (\bF^\top \otimes \bI_{p})}^2\\
\le~& \fnorm{(\bI_T\otimes \bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma^{\frac12})}^2 \opnorm{(\bF^\top \otimes \bI_{p})}^2\\
\le~& \trace\big[(\bI_T\otimes \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})\bM^\dagger (\bI_T\otimes \bSigma_{\hat{\mathscr{S}},\hat{\mathscr{S}}})\bM^\dagger\big]
\fnorm{\bF}^2 \\
\le~& \trace\big[ (\bI_T\otimes \bSigma_{\hat{\mathscr{S}},\hat{\mathscr{S}}})\bM^\dagger\big]
\fnorm{\bF}^2 \\
\le~& \trace\big[ (\bI_T\otimes \bSigma_{\hat{\mathscr{S}},\hat{\mathscr{S}}}) (\bI_T\otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})^\dagger)\big]
\fnorm{\bF}^2 \\
=~& T \trace\big[  \bSigma_{\hat{\mathscr{S}},\hat{\mathscr{S}}} (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})^\dagger\big]
\fnorm{\bF}^2 \\
\le~& T \trace\big[ (n\eta)^{-1} \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}})^\dagger \big]
\fnorm{\bF}^2 \\
\le~& \frac{T}{\eta} \fnorm*{\bF}^2,
\end{align*}

where the penultimate inequality uses $\bSigma_{\hat{\mathscr{S}},\hat{\mathscr{S}}} \preceq (n\eta)^{-1} \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}$ in the event $U_1\cap U_2$.
Therefore, on $U_1\cap U_2$, we have

\begin{align*}
\frac{1}{n}\sum_{ij}\fnorm*{\frac{\partial \bF/D}{\partial z_{ij}}}^2
&\le 4 \max(1, (\eta)^{-1} (T\wedge \frac{p}{n})) + 2  n^{-1}  [4 \max(1, (2\tau')^{-1})]^2 \\
&:= f(\eta, T, n, p),
\end{align*}

where the function $f$ is the same as in case (1). The only difference is that $\tau'$ in the upper bound for case (1) is replaced by $\eta$ in case (2).

§.§ Proofs of results in Appendixsec:Lipschitz-fix-X

The following proof of <Ref> relies on a similar argument as proof of <Ref>, we present the proof here for completeness.

Recall the KKT condtions for $\hbB$ defined in (<ref>):

1) For $k\in \hat{\mathscr{S}}$, , $\hbB{}^\top \be_k \ne\mathbf{0}$,

\be_k^\top\bX^\top\big[\bE-\bX(\hbB-\bB^*)\big] -n\tau\be_k^\top \hbB= \frac{n \lambda}{\|\hbB{}^\top\be_k\|_2} \be_k^\top \hbB
\quad \in\R^{1\times T}.

2) For $k\notin \hat{\mathscr{S}}$, , $\hbB{}^\top \be_k = \mathbf 0$,

\norm*{\be_k^\top\bX^\top\big[\bE-\bX(\hbB-\bB^*)\big] -n\tau\be_k^\top \hbB}< n\lambda.

Here the strict inequality is guaranteed by Proposition 2.3  of [Bellec, 2020].

Let $\ddot{\bB} = \frac{\partial \hbB}{\partial E_{it'}}$, $\dot\bE = \frac{\partial \bE}{\partial E_{it'}}$. Differentiation of the above display for $k\in\hat{\mathscr{S}}$ w.r.t. $E_{it'}$ yields

$$\be_k^\top\bX^\top(\dot\bE-\bX \ddot\bB) - n\tau \be_k^\top \ddot\bB
= n\be_k^\top \ddot\bB \bH^{(k)}$$

with $\bH^{(k)}
\lambda
\|\hbB{}^\top \be_k\|_2^{-1}\left(\bI_T - \hbB{}^\top\be_k \be_k^\top\hbB ~  \|\hbB{}^\top\be_k\|_2^{-2} \right)\in\R^{T\times T}$. Rearranging
and using $\dot\bE = \be_i \be_{t'}^\top$,

$$\be_k^\top \bX^\top \be_i \be_{t'}^\top = \be_k^\top[n\ddot \bB \bH^{(k)} + (\bX^\top\bX+n\tau\bI_{p\times p})\ddot\bB].$$

Recall $\bP_{\hat{\mathscr{S}}} = \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top\in\R^{p\times p}$. Multiplying by $\be_k$ to the left and summing over $k\in \hat{\mathscr{S}}$, we obtain

$$\bP_{\hat{\mathscr{S}}} \bX^\top \be_i \be_{t'}^\top =  n \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top \ddot \bB \bH^{(k)} + \bP_{\hat{\mathscr{S}}} (\bX^\top\bX+n\tau\bI_{p\times p})\ddot\bB,$$

which reduces to the following by $\supp(\ddot \bB)\subseteq \hat{\mathscr{S}}$ and $\bX\ddot\bB = \bX_{\hat{\mathscr{S}}}\ddot\bB$,

\begin{align*}
\bX_{\hat{\mathscr{S}}}^\top \be_i \be_{t'}^\top &= n \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top\ddot \bB \bH^{(k)} + \bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}}\ddot\bB \bI_T + n\tau\bP_{\hat{\mathscr{S}}} \ddot\bB \bI_T\\
&= n \sum_{k\in\hat{\mathscr{S}}} \be_k\be_k^\top\ddot \bB \bH^{(k)} + (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})
\ddot\bB \bI_T.
\end{align*}

Vectorizing the above yields

\begin{align*}
(\be_{t'} \otimes \bX_{\hat{\mathscr{S}}}^\top) \vec(\be_i) &= [n\sum_{k\in\hat{\mathscr{S}}} (\bH^{(k)}\otimes \be_k\be_k^\top) +\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}})] \vec(\ddot\bB) \\
&= (n \tbH +\bI_T \otimes (\bX_{\hat{\mathscr{S}}}^\top\bX_{\hat{\mathscr{S}}} + n\tau\bP_{\hat{\mathscr{S}}}) )\vec(\ddot\bB)\\
&= \bM\vec(\ddot\bB).
\end{align*}

A similar argument as in Proof of <Ref> leads to

\begin{align*}
\vec(\ddot \bB) = \bM^\dagger\bM \vec(\ddot \bB)
=\bM^{\dagger} (\be_{t'} \otimes \bX_{\hat{\mathscr{S}}}^\top) \be_i.
\end{align*}

Therefore, by $\bX\ddot{\bB}=\bX_{\hat{\mathscr{S}}}\ddot{\bB}$,

\begin{align*}
\frac{\partial F_{lt}}{\partial E_{it'} }
&= \be_l^\top \frac{\partial \bE - \bX(\hbB - \bB^*)}{\partial E_{it'}}\be_t\\
&= \be_l^\top \big( \be_i\be_{t'}^\top - \bX \ddot{\bB}\big) \be_t\\
&= \be_l^\top \be_i\be_{t'}^\top\be_t- \be_l^\top\bX \ddot{\bB}\be_t\\
&= \be_l^\top \be_i\be_{t'}^\top\be_t- (\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}})\vec(\ddot\bB)\\
&= \be_l^\top \be_i\be_{t'}^\top\be_t- (\be_t^\top \otimes \be_l^\top\bX_{\hat{\mathscr{S}}}) \bM^{\dagger} (\be_{t'} \otimes \bX_{\hat{\mathscr{S}}}^\top) \be_i\\
&= \be_l^\top \be_i\be_{t'}^\top\be_t- \be_l^\top(\be_t^\top \otimes \bX_{\hat{\mathscr{S}}}) \bM^{\dagger} (\be_{t'} \otimes \bX_{\hat{\mathscr{S}}}^\top) \be_i\\
&= \be_l^\top \be_i\be_{t'}^\top\be_t- \be_l^\top(\be_t^\top \otimes \bX) \bM^{\dagger} (\be_{t'} \otimes \bX^\top) \be_i,
\end{align*}

where the last equality is due to $\bM^\dagger = (\bI_T \otimes \bP_{\hat{\mathscr{S}}}) \bM^\dagger  (\bI_T \otimes \bP_{\hat{\mathscr{S}}})$.

Now the calculation of $\sum_{i=1}^n\frac{\partial F_{it}}{\partial E_{it'}}$ is straightforward,

\begin{align*}
\sum_{i=1}^n\frac{\partial F_{it}}{\partial E_{it'}} &= \sum_{i=1}^n \big[ \be_i^\top\be_i\be_t^\top\be_{t'}
- \be_i^\top (\be_t^\top \otimes \bX)\bM^\dagger (\be_{t'} \otimes \bX^\top) \be_i\big]\\
&= n\be_t^\top\be_{t'} - \trace[(\be_t^\top \otimes \bX)\bM^\dagger (\be_{t'} \otimes \bX^\top)]\\
&= n\be_t^\top\be_{t'} - \be_t^\top \hbA \be_{t'}\\
&= \be_t^\top (n\bI_T -\hbA)\be_{t'},
\end{align*}

where the third equality is due to the formula of $\hbA$ in (<ref>).

Noting that $\bF = \bE - \bZ\bH$, it follows that
$\sum_{i=1}^n\frac{\partial \be_i^\top\bZ\bH \be_t}{\partial E_{it'}} = \be_t^\top \hbA\be_{t'}$.

§.§ Proofs of results in Appendixsec:proba-tools

Let $\bz = \vec(\bE)$, then $\bz\sim \calN(\mathbf{0}, \bK)$ with $\bK = \bS \otimes \bI_n$ by Assumption <ref>.
For each $t_0, t_0' \in [T]$, let
$\bG^{(t_0, t_0')} = \bF  \be_{t_0'}\be_{t_0}^\top$,  and $\bff(\bz)^{(t_0, t_0')}  = \vec(\bG) \tD^{-1} $.
For convenience, we will drop the superscript ${(t_0, t_0')}$ from $\bG^{(t_0, t_0')}$  and $\bff(\bz)^{(t_0, t_0')}$	in this proof.
By $\trace(\bA^\top\bB) = \vec(\bA)^\top \vec(\bB)$, we obtain

\begin{align}\label{eq: a1}
\be_{t_0}^\top \bE^\top \bF \tD^{-1} \be_{t_0'} = \trace(\bE^\top \bF  \be_{t_0'}\be_{t_0}^\top) \tD^{-1}= \trace(\bE^\top \bG\tD^{-1}) = \bz^\top\bff(\bz).
\end{align}

By product rule, we have

\begin{align}\label{eq: nabla-f}
\nabla \bff(\bz) = \frac{\partial \vec(\bG) }{\partial \vec(\bE) } \tD^{-1} + \underbrace{\vec(\bG) \frac{\partial \tD^{-1}}{\partial \vec(\bE) }}_{\Rem},
\end{align}

where $\Rem = \bu\bv^\top$ with $\bu = \vec(\bG)\in \R^{nT\times 1}$, $\bv^\top = \frac{\partial \tD^{-1}}{\partial \vec(\bE) }\in \R^{1\times nT}$.
It follows that

\begin{align}\label{eq: stein-tr1}
\trace(\bK \nabla \bff(\bz)) = \trace\Big(\bK \frac{\partial \vec(\bG) }{\partial \vec(\bE)}\Big) \tD^{-1} + \trace(\bK \Rem).
\end{align}

Since $\bK = \bS \otimes \bI_n$ and $\bG = \bF  \be_{t_0'}\be_{t_0}^\top$,
$\bK_{it, lt'}= S_{tt'}I(i=l)$, and $G_{it} = F_{it_0'} I(t=t_0)$. It follows

\begin{equation}\label{eq: stein-tr2}
\begin{aligned}
\trace\Big(\bK \frac{\partial \vec(\bG) }{\partial \vec(\bE)}\Big)
= \sum_{i,t}\sum_{l,t'} \bK_{it, lt'} \frac{\partial G_{it}}{\partial E_{lt'} }
%&= \sum_{i}\sum_{t,t'} S_{tt'} \frac{\partial G_{it} }{\partial E_{it'}}\\
= \sum_{t'} S_{t_0t'} \sum_{i}\frac{\partial F_{it_0'}}{\partial E_{it'}}
%&= \sum_{t'} S_{t_0t'}  \be_{t_0'}^\top(n\bI_T - \hbA)\be_{t'}\\
= \be_{t_0}^\top \bS (n\bI_T - \hbA) \be_{t_0'},
\end{aligned}
\end{equation}

where the last equality used <Ref> and that $\hbA$ is symmetric.

Now we rewrite the quantity we want to bound as

\begin{align}
\fnorm{\bE^\top \bF/\tD - \bS (n\bI_T - \hbA )/\tD}^2\Bigr] \nonumber\\
=~& \sum_{t_0,t_0'} \E \Big[\Big( \be_{t_0}^\top \bE^\top \bF \tD^{-1} \be_{t_0'} - \be_{t_0}^\top \bS (n\bI_T - \hbA) \be_{t_0'}\tD^{-1} \Big)^2\Big]\nonumber\\
=~& \sum_{t_0,t_0'} \E \Big[ \big(\bz^\top \bff(\bz) - \trace(\bK \nabla \bff(\bz)) + \trace(\bK\Rem)\big)^2\Big] \nonumber\\
\le~& 2\sum_{t_0,t_0'} \Big\{ \E \Big[ \big(\bz^\top \bff(\bz) - \trace(\bK \nabla \bff(\bz)) \big)^2\Big] +  \E\Big[ \big(\trace(\bK\Rem)\big)^2 \Big] \Big\}\label{eq: LHS} ,
\end{align}

where the second equality follows from (<ref>), (<ref>) and (<ref>), and the last inequality uses elementary inequality $(a+b)^2 \le 2(a^2 + b^2)$.
We next bound the two terms in (<ref>).

First term in (<ref>).
By second-order Stein formula in <Ref>,

\begin{equation}\label{eq: stein}
\sum_{t_0,t_0'} \E \big(\bz^\top \bff(\bz) - \trace(\bK \nabla \bff(\bz))\big)^2
= \sum_{t_0,t_0'}
\E \Big[\fnorm{\bK^{\frac12}\bff(\bz)}^2  +  \trace\big[\big(\bK \nabla \bff(\bz)\big)^2\big] \Big].
\end{equation}

Now we bound the two terms in the right-hand side of (<ref>).
For the first term, recall  $\bff(\bz) = \vec(\bG) \tD^{-1} $, and $\bG = \bF  \be_{t_0'}\be_{t_0}^\top$, we obtain

\begin{align*}
\fnorm{\bK^{\frac12}\bff(\bz)}^2 =  \tD^{-2} \fnorm{(\bS^{\frac12}\otimes \bI_n)\vec(\bG)}^2  =  \tD^{-2} \fnorm{\bG\bS^{\frac12}}^2  = \tD^{-2} \fnorm{\bS^{\frac12}\be_{t_0}}^2  \fnorm{\bF\be_{t_0'}}^2.
\end{align*}

Summing over all $(t_0, t_0')\in [T] \times [T]$, we obtain

\begin{equation}\label{eq: stein-RHS1}
\sum_{t_0, t_0'} \fnorm{\bK^{\frac12}\bff(\bz)}^2 =  \tD^{-2} \fnorm{\bF}^2 \trace(\bS).
\end{equation}

For the second term in RHS of (<ref>), recall $ \nabla\bff(\bz) = \frac{\partial \vec(\bG) }{\partial \vec(\bE) } \tD^{-1} + \Rem$,

\begin{align}
&\trace\big[\big(\bK \nabla \bff(\bz)\big)^2\big] \nonumber\\
%&= \trace \Big\{\Big[\bK\Big(  \frac{\partial \vec(\bG) }{\partial \vec(\bE) } \tD^{-1} + \Rem\Big)\Big]^2\Big\}\nonumber\\
%=~& \trace \Big[\Big( \bK\frac{\partial \vec(\bG) }{\partial \vec(\bE) } \tD^{-1}\Big)^2\Big]  + \trace[(\bK\Rem)^2] + 2\trace\Big[ \bK\frac{\partial \vec(\bG) }{\partial \vec(\bE) } \tD^{-1} \bK \Rem\Big]\nonumber\\
=~& \tD^{-2}\trace \Big[\Big( \bK\frac{\partial \vec(\bG) }{\partial \vec(\bE)} \Big)^2\Big]  + \trace[(\bK\Rem)^2] + 2\tD^{-1}\trace\Big[ \bK\frac{\partial \vec(\bG) }{\partial \vec(\bE) } \bK \Rem\Big].\label{eq: *2s}
\end{align}

By property of vectorization operation,
$\vec(\bG) = \vec(\bF\be_{t_0'}\be_{t_0}^\top) = (\be_{t_0}\be_{t_0'}^\top \otimes \bI_n) \vec(\bF)$,

$$\frac{\partial \vec(\bG) }{\partial \vec(\bE)}
= (\be_{t_0}\be_{t_0'}^\top \otimes \bI_n)  \frac{\partial \vec(\bF) }{\partial \vec(\bE)}, $$

where $\opnorm{ \frac{\partial \vec(\bF) }{\partial \vec(\bE)} } \le 1$ since the map $\vec(\bE) \mapsto \vec(\bF)$ is 1-Lipschitz by <cit.>.

Now  we bound the three terms in (<ref>). For the first term, by Cauchy-Schwarz inequality,

\begin{align*}