# Meta Sparse Principal Component Analysis

Imon Banerjee111Departments of Statistics and Jean Honorio222Departments of
Computer Science
Purdue University
West-Lafayette, IN 47906
<EMAIL_ADDRESS>

###### Abstract

We study the meta-learning for support (i.e. the set of non-zero entries)
recovery in high-dimensional Principal Component Analysis. We reduce the
sufficient sample complexity in a novel task with the information that is
learned from auxiliary tasks. We assume each task to be a different random
Principal Component (PC) matrix with a possibly different support and that the
support union of the PC matrices is small. We then pool the data from all the
tasks to execute an improper estimation of a single PC matrix by maximising
the $l_{1}$-regularised predictive covariance to establish that with high
probability the true support union can be recovered provided a sufficient
number of tasks $m$ and a sufficient number of samples
$O\left(\frac{\log(p)}{m}\right)$ for each task, for $p$-dimensional vectors.
Then, for a novel task, we prove that the maximisation of the
$l_{1}$-regularised predictive covariance with the additional constraint that
the support is a subset of the estimated support union could reduce the
sufficient sample complexity of successful support recovery to $O(\log|J|)$,
where $J$ is the support union recovered from the auxiliary tasks. Typically,
$|J|$ would be much less than $p$ for sparse matrices. Finally, we demonstrate
the validity of our experiments through numerical simulations.

##  1 Introduction

Principal component analysis (PCA) is an important problem in high-dimensional
statistical learning with several applications in information theory [3],
image classification [21], brain computer interface [8], speech recognition
[29] among others. PCA seeks a linear transformation that maximises the
predictive covariance of the data, while depending on a small number of
variables. The learner faces several challenges in PCA, such as high dimension
or severe heterogeneity of the data. To be precise, the dimension of the data,
$p$, could be much higher than the sample size $n$, making estimation
infeasible. Or there could be limited samples from the distribution of
interest but a large amount of samples from multiple multivariate
distributions with different principal component (PC) matrices.

#### Practical Applications.

Even amongst the varied applications of PCA, there exists an important niche
concerned with estimating multiple PC matrices simultaneously. Here we present
some immediate examples. The recent paper by [18] produces a data driven
approach to map the temporal lobes of brain where they calculated $41$
separate covariance matrices each of which can be analysed using PCA. [27]
calculates the PC matrices of solar insolation over the Pacific region for
$264$ months from July 1983 to June 2005, simultaneously calculating $264$ PC
matrices. These matrices are often sparse, having only one or two significant
eigenvalues. Therefore, it is obvious that there are far-reaching consequences
of being able to simultaneously estimate multiple sparse PC matrices.

#### Motivation.

Theoretical results for consistency, rates of convergence, minimax risk bounds
for estimating eigenvectors and principal subspaces have been considered in
various prior works [30, 6, 5, 10, 9, 16]. It is well documented that the
classical PCA has major practical and theoretical drawbacks when it is applied
to high-dimensional data. The estimates PCA produces in high-dimensional
situations are often inconsistent [17, 14, 6, 7]. The loadings of PC matrices
are typically non-zero. This often makes it difficult to interpret the PC
matrix and identify important variables. For this reason, estimating the PC
matrix to recover its support set, which is the set of non-zero entries, is a
common strategy of structure learning. An estimate of the PC matrix will be
called sign consistent if it has the same support and sign of entries with
respect to the true PC matrix.

For the first challenge, we assume the PC matrices are sparse and consider a
general class of distributions (sub-Gaussian). To address the challenge of
heterogeneity, we consider a meta-learning learning problem where the learner
treats each different distribution as a task with a related PC matrix and
solves every task simultaneously. We suppose there are $m$ “auxiliary" tasks
and $n$ samples with dimension $p$ per task. When there is only one task
($m=1$), [7] proved that, under general assumptions, $n\in O(\log p)$ is
sufficient for the sign consistency of the PC matrix under $l_{1}$
regularisation. Similarly, when $m=1$, using a mixed convex norm, [20] proved
the consistency of PCA. However, neither of them performs support recovery
while simultaneously accounting for the heterogeneity of data.

In this paper, we solve the heterogeneity challenge with meta-learning. A
single task denotes the recovery of the $k$-dimensional principal subspace of
a given covariance matrix. We recover the support of the PC matrix in a novel
task with the information learned from other auxiliary tasks. We also use
improper estimation in our meta-learning method to obtain novel theoretical
guarantees for support recovery. Instead of estimating every PC matrix in the
auxiliary tasks, which might require more samples than feasible, we pool all
the samples from the auxiliary tasks together to estimate a single common PC
matrix and recover the “support union” (termed $J$). This allows us to provide
statistical guarantees for the sign consistency when $m\in O(\log p)$.
Moreover, given the recovered support union, we also provide statistical
guarantees for estimating any extra unimportant variables in a novel task when
the number of samples in the novel task is $O(\log|J|)$. This implies
substantial reductions in the number of samples required for the support
recovery of the novel task.

To the best of our knowledge, this is the first work concerning the sign-
consistency of sparse PC matrices under a meta-learning setup. We were able to
introduce randomness in the PC matrices of different tasks. Our theoretical
results hold for a wide class of distributions of the PC matrices under some
conditions, which implies broad application scenarios of our method. The use
of improper estimation in our method is innovative for the problem of support
recovery of high-dimensional PC matrices. Thereby, our work fills the gap
between theory and methodology of meta-learning in high-dimensional PCA.

Table 1: Notations used in the paper Notation | Description | Notation | Description
---|---|---|---
$[a]$ |  The first $a$ integers, i.e. $\left\\{1,2,\dots,a\right\\}$. | $|J|$ |  Cardinality of a set $J$.
$A_{i,j}$ |  The $i,j^{th}$ cell of matrix $A$. | $U_{i}$ |  The $i^{th}$ coordinate of a vector $u$.
$supp(A)$ |  Support of $diag(A)$, i.e. $\left\\{i:A_{i,i}\neq 0\right\\}$. | $A_{*,i}$ |  The $i^{th}$ column of matrix $A$.
$A_{i,*}$ |  The $i^{th}$ row of $A$. | $sign(x)$ |  The sign of $x,x\in\mathbb{R}$.
$\Sigma$ |  The true covariance $A$. | $\Sigma^{(i)}$ |  Covariance of auxiliary task $i$.
$\lambda_{k}(A)$ |  The $k^{th}$ largest eigenvalue of $A$ | $\|u\|_{q}$ |  $l_{q}$-norm of the vector $u\in^{p}$
$\|A\|_{q_{1},q_{2}}$ |  $\|\|A_{1,\infty}\|_{q_{1}},\dots,\|A_{p,*}\|_{q_{1}}\|_{q_{2}}$. | $\|A\|_{F}$ | The Frobenius norm of a matrix A
$\langle A,B\rangle$ |  $\mathrm{Trace}(A^{T}B)$. | $diag(A)$ |  Diagonal of matrix $A$
$\mathbb{D}(A)$ |  The diagonal matrix with elements $A_{i,i}$ | $\mathbb{I}_{p}$ |  The $p$-dimensional identity matrix
$A_{I,J}$ |  Submatrix of $A$ with rows $I$ and columns $J$ |  |

#### Contributions.

This paper has the following four contributions. In section 3, we propose a
meta-learning approach by introducing multiple auxiliary learning tasks for
support recovery of high-dimensional PC matrices with improper estimation. We
define a generative model for population covariance matrices that is amenable
for analysis. In section 4, we prove that $O\left(\frac{\log p}{m}\right)$ is
a sufficient sample complexity for support union recovery for $p$-dimensional
multivariate sub-Gaussian random vectors and $m$ auxiliary tasks with support
union $J$. In section 5, we prove that given a recovered support union from
the auxiliary tasks, we can reduce the dimension of the novel task from
$p\times p$ to $|J|\times|J|$. We then use this fact to establish a sample
complexity of $O(\log(|J|))$ for support recovery of the novel task. This
provides the theoretical basis for introducing more tasks for meta-learning in
support recovery of PC matrices. Lastly, we conduct synthetic experiments to
validate our theory. We calculate the support union recovery rates of our
meta-learning approach for different sizes of samples and tasks. For a fixed
task size $m$, our approach achieves high support union recovery rates when
the sample size per task has the order $O(\frac{\log p}{m})$. For a fixed
sample size per task, our method performs the best when the task size $m$ is
large.

##  2 Preliminaries

In this section, we introduce the definitions of the concepts important to our
meta-learning problem. A summary of notations used in the paper is illustrated
in Table 1.

#### Multivariate sub-Gaussian Distribution.

This paper is concerned with meta-learning of multiple principal component
matrices, each of which are generated randomly. To make inferences about those
matrices in sufficient generality, we must first formally define the class of
sub-Gaussian matrices with random covariance matrices.

######  Definition 2.1.

For $j\in[n^{(i)}]$ and $i\in[m]$, we say $X_{j}^{(i)}\in^{p}$ comes from a
family of random $p$-dimensional multivariate sub-Gaussian with randomised
covariance matrices $\left\\{\mathcal{G}^{(i)}\right\\}$ distribution of size
$m$ and parameter $\sigma$ if

  1. 1.

$\mathbb{E}[X_{j}^{(i)}|\mathcal{G}^{(i)}]=0$,
$\mathrm{Cov}(X_{j}^{(i)}|\mathcal{G}^{(i)})=\mathcal{G}^{(i)}$.

  2. 2.

Given $\mathcal{G}^{(i)}$, $\frac{X_{l,j}^{(i)}}{\mathcal{G}_{l,l}^{(i)}}$ is
sub-Gaussian with parameter $\sigma$, $\forall i\in[m]$ and $j\in[n^{(i)}]$
and $l\in[p]$.

  3. 3.

Given $\mathcal{G}^{(i)}$, $X_{j}^{(i)}$ are independent random variables.

######  Remark 2.2.

We would like to point out that such definitions are standard in the
literature and have been used extensively in prior work (see, for example,
[26]) for modelling sub-Gaussian random variables with randomised covariance
matrices under the meta-learning setup. We have deferred the definition of
Sub-Gaussian and Sub-exponential distributions to Appendix A.

#### Problem Formulation.

In this paper, we have focused on estimating the support of the PC matrix of a
multivariate sub-Gaussian distribution. We first estimate a superset of the
support of the PC matrices. Then we solve a novel task utilising the recovered
support union.

To be precise, let for $i\in[m]$, $\Sigma$ be the true underlying covariance
matrix and $\Sigma^{(i)}$ be perturbations of $\Sigma$, with each
$\Sigma^{(i)}$ forming the covariance matrix of an auxiliary task $i$. Then,
we are interested in using those auxiliary tasks to recover support union $J$
of the PC matrix $\Pi$ corresponding to $\Sigma$. It is noteworthy that by
support union we refer to $supp(\Pi)$ and not a union of supports of each
individual auxiliary task. This is because we are interested in a novel task
where the support of its PC matrix is assumed to be a subset of the support of
$\Pi$. This allows us to achieve a substantial reduction in sample complexity
from $O(\log(p))$ to $O(\log|J|)$ while estimating the support in the novel
task.

##  3 Our Meta-Learning Model

In this section, we define a generative model for the population covariance
matrices. Let $\Sigma$ be a positive semi-definite, real symmetric matrix. For
some orthonormal matrix U, and positive diagonal matrix $\Lambda$, the
spectral decomposition of $\Sigma$ gives us,

$\displaystyle\Sigma=U\Lambda
U^{T}:=\sum_{l=1}^{p}\lambda_{l}(\Sigma)U_{*,l}U_{*,l}^{T},$

where $U_{*,l}$ are orthonormal vectors and $\lambda_{i}$ are the diagonal
elements of $\Lambda$. If
$\lambda_{l}(\Sigma)>\lambda_{l+1}(\Sigma)\forall\enspace l$, then,

$\displaystyle~{}\Pi=\sum_{l=1}^{k}U_{*,l}U_{*,l}^{T}$ (1)

gives us the unique $k$-dimensional principal subspace of $\Sigma$. With
random matrices $R^{(i)}$ and random matrices $D^{(i)}$, we generate a
sequence of random matrices $\Sigma^{(1)},\dots,\Sigma^{(m)}\in^{p}$. For each
$i\in[m]$, the covariance matrix $\Sigma^{(i)}$ corresponding to the $i$-th
auxiliary task is generated as,

$\displaystyle\Sigma^{(i)}$
$\displaystyle=R^{(i)}U(\Lambda+D^{(i)})U^{T}\left(R^{(i)}\right)^{T}~{}.$ (2)

Referring to Definition 2.1, we assume that for $j\in[n^{(i)}]$, our data,
$X_{j}^{(i)}$ is a random p-dimensional multivariate sub-Gaussian random
variable with randomised covariance matrix $\Sigma^{(i)}$. Let
$S^{(i)}:=\frac{1}{n^{(i)}}\sum_{j=1}^{n^{(i)}}X_{j}^{(i)}(X_{j}^{(i)})^{T}$
be the sample covariance matrix corresponding to the auxiliary task
$\Sigma^{(i)}$. Recall from [13] that an improper estimator of a parameter is
defined to be any out-of-model predictor for that parameter. A “pooled"
estimator $S$ forms the improper estimation of $\Sigma$ as defined below:

$\displaystyle~{}S:=\frac{1}{m}\sum_{i=1}^{m}S^{(i)}.$ (3)

###  3.1 The Predictive Covariance Loss Function

In this subsection, we describe the objective function that we maximise to
obtain the improper estimator. Recall from Theorem 1 in [15] that PCA can be
interpreted as a covariance maximisation technique maximising the predictive
covariance $\mathrm{Trace}\left(\mathrm{Cov}(X,HX|H)\right)$. Observe that,

$\displaystyle\mathrm{Trace}\left(\mathrm{Cov}\left(X,HX|H\right)\right)=\langle\Sigma,H\rangle,$

for $\mathrm{Cov}(X)=\Sigma$ and a rank $k$ orthogonal matrix $H$.
Accordingly, the $k$-dimensional sparse principal component seeks the rank-$k$
orthogonal projection matrix $H$ such that for a penalty $\rho>0$ the
penalised predictive covariance

$\displaystyle~{}\langle\Sigma,H\rangle-\rho\|H\|_{1,1}$ (4)

is maximised. Crucial to getting a consistent estimate of the principal
component matrix is its identifiability. Unless the PC matrix is unique, there
is no hope to correctly estimate the PC matrix. For example, if
$\Sigma=\mathbb{I}_{p}$, then every $k$-dimensional subspace of p is a valid
principal component. Slightly abusing notation, by $\lambda_{k}$ we denote
$\lambda_{k}(\Sigma)$ and by $\lambda_{k+1}$ we denote
$\lambda_{k+1}(\Sigma)$. We make our first assumption about the uniqueness of
the PC matrix, denoted in terms of the spectral gap of the covariance matrix
$\Sigma$.

######  Assumption 3.1 (Uniqueness).

Let $\lambda_{k}$ and $\lambda_{k+1}$ be the $k$-th and $k+1$-th eigenvalues
of $\Sigma$. Then, $\lambda_{k}-\lambda_{k+1}>0$.

######  Remark 3.2.

This condition ensures that the principal subspace is identifiable.
Subsequently, we need to define the notion of sparsity of the PC matrix.
Recall from Equation 1 that $\Pi$ is a positive semi-definite matrix.
Consequently, the $i$-th row/column of $\Pi$ is $0$ is non-zero if and only if
$\Pi_{i,i}\neq 0$. Hence, it is sufficient to recover the set of non-zero
elements in the diagonal of $\Pi$. Hence, we can make our sparsity assumption.

######  Assumption 3.3 (Sparsity).

Let $\Pi$ be the $k$-dimensional PC matrix corresponding to $\Sigma$. Then for
some support $J$, we have $supp(\Pi)\subseteq J.$

Throughout Section 4, we will assume that $\Sigma$ satisfies both 3.1 and 3.3.

###  3.2 The Improper Objective Function

As illustrated in Section 2, we recover the support union by estimating the
true covariance matrix. To be specific, we pool together all the samples from
$m$ tasks and maximise the $l_{1}$-regularised predictive covariance between
the true and the estimated covariance matrices, i.e., we solve the following
optimisation problem with regularisation constant $\rho>0$

$\displaystyle\hat{\Pi}:=\arg\sup_{H}\sum_{i=1}^{m}T_{i}\left\\{\langle
S^{(i)},H\rangle-\rho\|H\|_{1,1}\right\\}\qquad\text{ subject to
}H\in\mathcal{F}^{k},\text{ where, }$ (5)
$\displaystyle\mathcal{F}^{k}:=\\{M\in^{p\times p}:0\preceq
M\preceq\mathbb{I}_{p}\text{ and }\mathrm{Trace}(M)=k\\},$

and $T_{i}$’s are weights proportional to the number of samples $n^{(i)}$ for
each $i$. For clarity of exposition, we assume that $n^{(i)}=n$ for all values
of $i$ for the reminder of this paper. This allows us to assume
$T^{(i)}=\frac{1}{m}$ without losing generality. Consequently, we can rewrite
the objective function in (5) as

$\displaystyle\hat{\Pi}=\arg\sup_{H}\langle S,H\rangle-\rho\|H\|_{1,1}\text{
subject to }H\in\mathcal{F}^{k},$ (6)

where $S$ is as defined in Equation 3. It is noteworthy that Equation 5 is an
improper estimation technique since it estimates a single PC matrix with data
from different distributions. This enables us to efficiently recover the true
support with an optimal number of sample size per task (see Section 4).

For the next step, suppose that we have successfully recovered the true
support union $J$ in the first step. Then for the novel task, ($m+1$-th task),
we can assume that the support of its PC matrix is also a subset of the
recovered support union. Then, for a noisy estimate $S^{(m+1)}$ of
$\Sigma^{(m+1)}$ we propose the following $l_{1}$-regularised predictive
covariance loss function

$\displaystyle\langle S^{(m+1)},H\rangle-\rho\|H\|_{1,1}\text{ subject to
}H\in\mathcal{F}^{k}\text{ and }supp(H)\subseteq J.~{}$ (7)

We note that Equation 7 is also a form of improper estimation due to presence
of extra constraint. This improper estimation technique reduces the dimension
of the optimisation problem from $p\times p$ to $|J|\times|J|$.

##  4 Support Union Recovery for the Auxiliary Tasks

In this section, we formally show that a sample complexity of
$O\left(\frac{\log p}{m}\right)$ is sufficient for the recovery of the support
union via optimising Equation 6. One key task in achieving the said sample
complexity is to carefully upper bound the tail probabilities of
$\left\|S-\Sigma\right\|_{\infty,\infty}$. To that end, we prove the following
tail probability for a family of multivariate sub-Gaussian with randomised
covariance matrices. Let
$S^{(i)}=\frac{1}{n}\sum_{t=1}^{n}X_{t}^{(i)}(X_{t}^{(i)})^{\text{T}}$.

######  Lemma 4.1.

Let $\Sigma^{(i)}$ be a sequence of random covariance matrices such that
$\left\|\Sigma^{(i)}\right\|_{\infty,\infty}\leq\eta$. Then, for
$\\{X_{j}^{(i)}\\}_{1\leq j\leq n,1\leq i\leq m}$ following a family of random
$p$-dimensional multivariate sub-Gaussian distributions with parameter
$\sigma$ and random covariance matrices $\Sigma^{(i)}$,

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}\geq\epsilon\right)\leq\frac{p(p+1)}{2}e^{-\frac{nm\epsilon^{2}}{512\sigma^{4}\eta}},$

whenever $0<\epsilon<32\eta\sigma^{2}$.

The proof of this lemma is presented in Section B.4. Recall from Equation 2
that we defined each auxiliary task $\Sigma^{(i)}$ in terms of random matrices
$R^{(i)}$ and random matrices $D^{(i)}$. Key to the process of improper
estimation is the control of the maximal loss between the auxiliary tasks
$\Sigma^{(i)}$ and the true underlying task $\Sigma$. In other words, we need
to find theoretical guarantees for
$\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}$. We achieve it by
controlling the variance and support of the perturbations $R^{(i)}$ and
$D^{(i)}$. To that end, we make the following assumption.

######  Assumption 4.2.

Assume that for each $1\leq i\leq m$, $R^{(i)}$ is a random matrix, and
$D^{(i)}$ is a random matrix. Then,

  1. 1.

Each $R^{(i)}$ is independent and identically distributed. Moreover, for a
constant $\mathbb{C}_{R}\geq 0$,

$\displaystyle\mathbb{E}\|R^{(i)}-\mathbb{I}_{p}\|_{\infty,1}^{2}\leq\frac{\mathbb{C}_{R}}{p}\text{
and,}$ (8)
$\displaystyle\|R^{(i)}-\mathbb{I}_{p}\|_{\infty,1}\leq\mathbb{C}_{R}\text{
almost surely.}$ (9)

  2. 2.

Each $D^{(i)}$ is independent and identically distributed with
$\mathbb{E}[D^{(i)}]=0$. Moreover, for some constant $L>0$,

$\displaystyle\lambda_{1}(D^{(i)})<L\text{ almost surely.}~{}$ (10)

######  Remark 4.3.

White noise additive models have been used extensively in prior literature to
study precision and covariance matrices [25, 26]. It is noteworthy that under
our assumptions, we can recover additive models from Equation 2 when
$R^{(i)}=\mathbb{I}_{p}$ for every $i$. However, $R^{(i)}$ suffers from curse
of dimensionality. Equation 8 provides a necessary shrinkage property to
$R^{(i)}$. It is noteworthy that such a property is necessary to ensure
consistency of $\frac{1}{m}\sum_{i=1}^{m}\Sigma^{(i)}$ as an estimator of
$\Sigma$. Our final technical assumption is given below.

######  Assumption 4.4.

If $\Sigma$ is a covariance satisfying 3.3. Then,

$\displaystyle\frac{8|J|}{\lambda_{k}-\lambda_{k+1}}\lVert\Sigma_{J^{c},J}\rVert_{2,\infty}<1.$

This condition previously appears in [7] and recovers the conditions by [1]
when $\Sigma_{J^{c},J}=0$. It is also similar to the irrepresentability
condition in the seminal papers [28, 11, 12] for $l_{1}$ penalised sparse
regression and lasso type recovery. As far as we know, this is the most
general condition for the support recovery of PC matrices and encompasses the
conditions when $\Sigma$ is block diagonal. For ease of exposition, we define
the following notation:

$\displaystyle\lambda_{diff}:=\lambda_{k}-\lambda_{k+1},$
$\displaystyle\lambda^{\dagger}:=2\left(\mathbb{C}_{r}+1\right)^{2}\left(\lambda_{1}(\Sigma)+L\right),$
$\displaystyle\rho_{1}:=4\lambda^{\dagger}\sqrt{\frac{\log(p)}{m}}+\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma),$
$\displaystyle\rho_{2}:=16\sqrt{2\sigma^{4}\lambda^{\dagger}\frac{\log\left(p+1\right)}{nm}}.$

The following theorem establishes the $O(\frac{\log p}{m})$ sample complexity
for the support union recovery from multiple auxiliary tasks.

######  Theorem 4.5 (Support Union Recovery).

Let 4.2 and 4.4 hold. Then there exists a large enough constant
$\mathbb{C}_{\pi}>0$ and a constant $\alpha\in(0,1)$ such that, whenever

$\displaystyle\alpha\max\left\\{\rho_{1},\rho_{2}\right\\}<\rho<\min\left\\{\frac{\left(\lambda_{diff}\right)\min_{j\in
J}{\Pi_{j,j}}}{16|J|},\frac{\left(\lambda_{diff}\right)^{2}}{4|J|\left(\lambda_{diff}+8\lambda_{1}(\Sigma)\right)}\right\\}$

then, with probability at least $1-\frac{2}{p^{2}}-\frac{1}{2(p+1)^{2}}$ the
solution $\hat{\Pi}$ to the objective function Equation 6 satisfies,

$\displaystyle 1.\enspace supp(\hat{\Pi})=J,$ $\displaystyle
2.\left\|\hat{\Pi}-\Pi\right\|_{\infty,\infty}\leq\mathbb{C}_{\pi}\rho.$

In other words, Theorem 4.5 gives us our required sample complexity for sign
consistency and error bound.

###  4.1 Sketch of Proof for Theorem 4.5

The complete proof of this theorem may be found at Section B.1 and we only
make the outline here. The proof of Theorem 4.5 is via primal-dual witness
method. It starts by writing down the objective function Equation 6.

$\displaystyle\max_{H\in\mathcal{F}^{k}}\langle S,H\rangle-\rho\|H\|_{1,1}.$

Define the subspace

$\mathbb{B}_{p}:=\left\\{Z\in^{p\times
p}:diag(Z)=0,Z=Z^{T},\left\|Z\right\|_{\infty,\infty}\leq 1\right\\}.$

Using strong-duality of $\|\cdot\|_{1,1}$, we can rewrite the objective
function as,

$\displaystyle\max_{\begin{subarray}{c}H\in\mathcal{F}^{k}\\\
Z\in\mathbb{B}\end{subarray}}\langle S,H\rangle-\rho\langle
H,Z+\mathbb{I}_{p}\rangle\equiv\
\max_{\begin{subarray}{c}H\in\mathcal{F}^{k}\\\
Z\in\mathbb{B}\end{subarray}}\langle S,H\rangle-\rho\langle H,Z\rangle-\rho
k\equiv\ $ $\displaystyle\max_{\begin{subarray}{c}H\in\mathcal{F}^{k}\\\
Z\in\mathbb{B}\end{subarray}}\langle H,S-\rho Z\rangle,$ (11)

where the last step follows from the fact that $k$ is a constant. To maximise
Equation 11, we need to find solutions $(\hat{H},\hat{Z})$ such that Karush-
Kuhn-Tucker (KKT) [2] conditions

$\displaystyle\hat{Z}_{i,j}=sign(\hat{H}_{i,j})$ $\displaystyle\forall\enspace
i\neq j,\hat{H}_{i,j}\neq 0,~{}$ (12)
$\displaystyle\hat{Z}_{i,j}\in[-1,1]\forall$ $\displaystyle\enspace i\neq
j,\hat{H}_{i,j}=0,~{}$ (13) $\displaystyle\hat{H}=\arg\max$
$\displaystyle\langle S-\rho\hat{Z},H\rangle,~{}$ (14)

are satisfied. Proceeding according to the primal-dual witness method, we
construct the constrained optimisation problem where the $supp(H)\subseteq J$.

$\displaystyle\max_{H\in\mathcal{F}^{k},supp(H)\subseteq J}\langle
S,H\rangle-\rho\|H\|_{1,1}.~{}$ (15)

Let $(\tilde{H},\tilde{Z})$ be the corresponding primal and dual variables.
The proof selects an appropriate $|J|\times|J|$ orthonormal matrix $Q$
according to Lemma A.6 with the correct rotation such that the Frobenius norm
between the true eigenvectors $U_{J,*}$ and the estimated eigenvectors
$\hat{U}_{J,*}$ is small. Then we propose the following solution
$(\hat{H},\hat{Z})$ to the dual objective Equation 11. It is noteworthy that
this construction is similar to that in Theorem 1 in [7].

$\displaystyle\hat{H}=\tilde{H},$
$\displaystyle\hat{Z}_{J,J}=\tilde{Z}_{J,J},$ (16)
$\displaystyle\hat{Z}_{ij}=\frac{1}{\rho}\bigl{\\{}S_{ij}-\langle
Q_{i,*},\Sigma_{J,j}\rangle\bigr{\\}},(i,j)\in J\times J^{c},$ (17)
$\displaystyle\hat{Z}_{ij}=\frac{1}{\rho}(S-\Sigma)_{ij},(i,j)\in J^{c}\times
J^{c},i\neq j.$ (18)

The rest of the proof follows in $4$ parts.

  1. 1.

Finding an appropriate $Q$ matrix and an appropriate tail probability bound
under Assumption 4.2 by using a combination of Matrix Chernoff bound (see
Appendix A) and Lemma 4.1.

  2. 2.

Proving that the proposed solution satisfies the KKT conditions (12), (13),
and (14) with high probability for our choices of $\rho$ and is unique.

  3. 3.

Find an upper bound for estimation error. We achieve this through the
following upper bound,
$\left\|\hat{H}-\Pi\right\|_{\infty,\infty}<2\frac{\rho|J|}{\lambda_{diff}}$,
which is a mild variation of a similar upper bound from Lemma 2 in [24]. Thus,
for sufficiently small enough $\rho$, we prove that
$sign(\hat{\Pi})=sign(\Pi)$, and our estimate is sign consistent.

  4. 4.

Showing that the proposed solution is sign consistent using a similar
technique to part 3.

##  5 Support Recovery for the Novel Task

In this section, we show that $O\left(\log|J|\right)$ sample complexity is
sufficient for the support recovery for a novel task where $J$ is the
recovered support union in Section 4. We also show that for the novel task,
the estimation problem can be reduced from estimating a $p\times p$
dimensional matrix to a $|J|\times|J|$ dimensional matrix. Let $n^{(m+1)}$ be
the number of samples in the novel task. We begin by establishing the
reduction in dimension. Let $\Sigma^{(m+1)}$ be the variance matrix
corresponding to the novel task with spectral decomposition

$\displaystyle U^{(m+1)}\Lambda^{(m+1)}\left(U^{(m+1)}\right)^{T},$

and the $k$-dimensional principal component matrix $\Pi^{(m+1)}$. Throughout
the rest of this section, we assume that $\Sigma^{(m+1)}$ satisfies 3.1 and
3.3. We write $\Sigma^{(m+1)}$ and its estimator $S^{(m+1)}$ as the following
block matrices

$\displaystyle~{}\Sigma^{(m+1)}=\begin{bmatrix}\Sigma_{J,J}^{(m+1)}&\Sigma_{J,J^{c}}^{(m+1)}\\\
\Sigma_{J^{c},J}^{(m+1)}&\Sigma_{J^{c},J^{c}}^{(m+1)}\end{bmatrix},$
$\displaystyle\quad
S^{(m+1)}=\begin{bmatrix}S_{J,J}^{(m+1)}&S_{J,J^{c}}^{(m+1)}\\\
S_{J^{c},J}^{(m+1)}&S_{J^{c},J^{c}}^{(m+1)}\end{bmatrix}.$ (19)

We obtain the $k$-dimensional principal component of the novel task
$\Pi^{(m+1)}$ by optimising the objective function,

$\displaystyle\max_{H\in\mathcal{F}^{k}}\langle
S^{(m+1)},H\rangle-\rho\|H\|_{1,1}.~{}~{}$ (20)

However, we have already recovered the support union $J$. Hence, we can add
another constraint on $H$: $supp(H)=J$. Moreover, since the row and column $i$
of the principal component $\Pi^{(m+1)}$ is $0$ if and on
$\Pi_{i,i}^{(m+1)}=0$, we can rewrite this additional constraint as,

$\displaystyle H=\begin{bmatrix}H_{J,J}&0\\\ 0&0\end{bmatrix}.$

Since $H\in\mathcal{F}^{k}$, the maximisation problem only makes sense if
$|J|>k$. This also automatically implies that $H\in\mathcal{F}_{J}^{k}$, where
$\mathcal{F}_{J}^{k}:=\left\\{M\in^{|J|\times|J|}:0\preceq
M\preceq\mathbb{I}_{|J|}\text{ and }\mathrm{Trace}(M)=k\right\\}$. Therefore,
the objective corresponding to the novel task maximises $\langle
S^{(m+1)},H\rangle-\rho\|H\|_{1,1}$ subject to $H\in\mathcal{F}^{k}$ and
$supp(H)\subseteq J$. This is equivalent to maximising

$\displaystyle{\hat{\Pi}^{(m+1)}}:=\arg\sup_{H_{J,J}}\ \langle
S_{J,J}^{(m+1)},H_{J,J}\rangle-\rho\|H_{J,J}\|_{1,1},\text{ for
}H_{J,J}\in\mathcal{F}_{J}^{k}~{}$ (21)

where $supp(\hat{\Pi})=J$ from Theorem 4.5 is the recovered support union, and
$H_{J,J},S_{J,J}\in^{|J|\times|J|}$. The following proposition on the
eigenvalues and eigenvectors of $\Sigma_{J,J}^{(m+1)}$ shows that the
estimation problem can be reduced from estimating a $p\times p$ dimensional
matrix to a $|J|\times|J|$ dimensional matrix.

######  Proposition 5.1.

For $1\leq l\leq k$, $U_{l,J}^{(m+1)}$ are the first $k$ eigenvectors of
$\Sigma_{J,J}^{(m+1)}$ with corresponding eigenvalues $\lambda_{l}^{(m+1)}$.

The proof of this proposition is given in Section B.2. Following from the fact
that $\Sigma^{(m)}$ satisfies 3.1 and 3.3, Proposition 5.1 establishes as an
immediate consequence that

$\displaystyle\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}>0,\text{ and,
}supp(\Pi^{(m+1)})\subseteq J^{(m+1)}.$

This establishes that the $k$-dimensional PC matrix of $\Sigma_{J,J}^{(m+1)}$
is sparse and unique, as well as the solution to a lower dimensional objective
function.

Overloading notation, define
$\lambda_{diff}^{(m+1)}:=\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}$. We make
the following technical assumption on the elements of $\Sigma_{J,J}^{(m+1)}$.

######  Assumption 5.2.

Let $\Sigma^{(m+1)}$ be the covariance matrix corresponding to the novel task
and $J$ be the recovered support union. Then,

$\displaystyle\frac{8|J^{(m+1)}|}{\lambda_{diff}^{(m+1)}}\lVert\Sigma^{(m+1)}_{(J^{(m+1)})^{c},J^{(m+1)}}\rVert_{2,\infty}<1.$

We define

$\displaystyle\mathbb{C}_{n1}:=\frac{1}{32\lambda_{1}\left(\Sigma^{(m+1)}\sigma^{2}\right)},\quad\mathbb{C}_{n2}:=\frac{4|J^{(m+1)}|\left(1+\frac{8\lambda_{1}^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\right)}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}},$
$\displaystyle\mathbb{C}_{n3}:=\frac{\left(\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}\right)\min_{j\in\left[|J^{(m+1)}|\right]}\Pi_{j,j}^{(m+1)}}{4|J^{(m+1)}|}.$

The following theorem shows the $O(\log|J|)$ sample complexity for support
recovery of the novel task.

######  Theorem 5.3 (Novel Task Support Recovery).

Let $\Sigma^{(m+1)}$ be the covariance matrix corresponding to the novel task
satisfying 5.2 and $J$ be the recovered support union. Then, there exists a
constant $\alpha^{(m+1)}\in(0,1)$ such that if

$\displaystyle
n^{\frac{1}{3}}>\max\left\\{\mathbb{C}_{n1},\mathbb{C}_{n2},\mathbb{C}_{n3}\right\\}\text{
and,
}\frac{1}{n^{1/3}\alpha^{(m+1)}}<\rho<\frac{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}{4|J^{(m+1)}|\left(1+\frac{8\lambda_{1}^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\right)},$

then, with probability,
$1-\frac{|J|(|J|+1)}{2}e^{-\frac{n^{1/3}}{8\mathbb{C}\sigma^{2}\left(2\lambda_{1}(\Sigma^{(m+1)})\right)}}$,
the maximiser $\hat{\Pi}^{(m+1)}$ to the objective function Equation 7
satisfies, $supp(\hat{\Pi}^{(m+1)})=J^{(m+1)}$.

The proof of this theorem is deferred to Section B.3

##  6 Experiments

In this section, we empirically validate our theoretical claims in Sections 4
and 5.

#### Gaussian Distribution Setting.

For all the experiments in this subsection, we let $p=50$, $k=5$ and $|J|=5$
and perform 100 repetitions of each setting. We first consider the setting of
different sample sizes. We choose $n\in\left\\{3,5,7,9\right\\}$ and use
$\rho=\sqrt{\frac{\log(p+1)}{mn}}$ for all pairs of $(m,n)$. The number of
tasks $m$ is rescaled to $T$ defined by $\frac{mn}{\log(p+1)}$. For each
$j\in[n]$, $X^{(i)}_{j}$ is generated from a Gaussian distribution with mean a
$p$-dimensional $0$ vector and covariance matrix $\Sigma^{(i)}$. Figure 1
depicts the outcome of our experiments. For different choices of $n$, the
graphs overlap each other perfectly (both $\mathbb{P}(\hat{J}=J)$ and
$\|{\hat{\Pi}-\Pi}\|_{\infty,\infty}$). Then, we took the recovered support
union and derived the probability of support recovery for a novel task.

We perform the experiments when there are $2,3$, or $4$ extra zeros and all of
those have greater than $95\%$ probability of accurately identifying the extra
zeros. Moreover, in Appendix C we provide more details and show extra results
for the uniform and the mixture of sub-Gaussian settings.

#### Real World Experiments.

In Appendix D, we show that our method obtains anywhere between $85-95\%$
accuracy in a real world brain-imaging dataset.

(a) Probability of Support Union Recovery (b) Maximal Error Bound (c) Support
Recovery of Novel Task

Figure 1: Left, Center: Simulations for Theorem 4.5 on the probability of
support union recovery and maximal error under various settings of $n$.
$\rho=\sqrt{\frac{\log(p+1)}{mn}}$. The x-axis is set by
$T:=\frac{mn}{\log(p+1)}$ varying from $\left\\{1,\dots,15\right\\}$ for
support union recovery. Right: Probability of support recovery for novel task.
$\rho=\sqrt{\frac{\log(|J|+1)}{n}}$. The x-axis is set by
$T:=\frac{n}{\log(|J|+1)}$ varying from $\left\\{10,\dots,25\right\\}$

##  7 Future Studies

There still exists some open problems regarding Meta Sparse PCA that were not
addressed in this work. Dependencies in the data, such as when the data is
from a time series makes an interesting such case we do not address.
Similarly, it is also unclear how estimation should proceed when there is
missing data. Finally, we would like to point out that there are newer
techniques such as Kernel PCA which can describe non-linear boundaries and it
is unclear how Meta sparse PCA matrices should proceed under such
circumstances.

## References

  * [1] Arash A Amini and Martin J Wainwright.  High-dimensional analysis of semidefinite relaxations for sparse principal components.  In 2008 IEEE international symposium on information theory, pages 2454–2458. IEEE, 2008.
  * [2] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [3] Yash Deshpande and Andrea Montanari.  Information-theoretically optimal sparse pca.  In 2014 IEEE International Symposium on Information Theory, pages 2197–2201, 2014.
  * [4] Jean Honorio and Tommi Jaakkola.  Tight bounds for the expected risk of linear classifiers and pac-bayes finite-sample guarantees.  In Artificial Intelligence and Statistics, pages 384–392. PMLR, 2014.
  * [5] Rodolphe Jenatton, Guillaume Obozinski, and Francis Bach.  Structured sparse principal component analysis.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 366–373. JMLR Workshop and Conference Proceedings, 2010.
  * [6] Iain M Johnstone and Arthur Yu Lu.  On consistency and sparsity for principal components analysis in high dimensions.  Journal of the American Statistical Association, 104(486):682–693, 2009.
  * [7] Jing Lei and Vincent Q Vu.  Sparsistency and agnostic inference in sparse pca.  The Annals of Statistics, 43(1):299–322, 2015.
  * [8] Chin-Teng Lin, Yu-Chieh Chen, Teng-Yi Huang, Tien-Ting Chiu, Li-Wei Ko, Sheng-Fu Liang, Hung-Yi Hsieh, Shang-Hwa Hsu, and Jeng-Ren Duann.  Development of wireless brain computer interface with embedded multitask scheduling and its application on real-time driver’s drowsiness detection and warning.  IEEE Transactions on Biomedical Engineering, 55(5):1582–1591, 2008\.
  * [9] Karim Lounici.  Sparse principal component analysis with missing observations.  In High dimensional probability VI, pages 327–356. Springer, 2013\.
  * [10] Zongming Ma.  Sparse principal component analysis and iterative thresholding.  The Annals of Statistics, 41(2):772–801, 2013.
  * [11] Nicolai Meinshausen and Peter Bühlmann.  High-dimensional graphs and variable selection with the lasso.  The annals of statistics, 34(3):1436–1462, 2006.
  * [12] Nicolai Meinshausen and Bin Yu.  Lasso-type recovery of sparse representations for high-dimensional data.  The annals of statistics, 37(1):246–270, 2009.
  * [13] Jaouad Mourtada and Stéphane Gaïffas.  An improper estimator with optimal excess risk in misspecified density estimation and logistic regression.  J. Mach. Learn. Res., 23:31–1, 2022.
  * [14] Boaz Nadler.  Finite sample approximation results for principal component analysis: A matrix perturbation approach.  The Annals of Statistics, 36(6):2791–2817, 2008.
  * [15] Michael L Overton and Robert S Womersley.  On the sum of the largest eigenvalues of a symmetric matrix.  SIAM Journal on Matrix Analysis and Applications, 13(1):41–45, 1992\.
  * [16] Seyoung Park and Hongyu Zhao.  Sparse principal component analysis with missing observations.  The Annals of Applied Statistics, 13(2):1016–1042, 2019.
  * [17] Debashis Paul.  Asymptotics of sample eigenstructure for a large dimensional spiked covariance model.  Statistica Sinica, pages 1617–1642, 2007.
  * [18] Andrew S Persichetti, Joseph M Denning, Stephen J Gotts, and Alex Martin.  A data-driven functional mapping of the anterior temporal lobes.  Journal of Neuroscience, 41(28):6038–6049, 2021.
  * [19] David Pollard.  PROCESSES: THEORY AND APPLICATIONS.  1984\.
  * [20] Xin Qi, Ruiyan Luo, and Hongyu Zhao.  Sparse principal component analysis by choice of norm.  Journal of multivariate analysis, 114:127–160, 2013.
  * [21] Yaoqi Sun, Liang Li, Liang Zheng, Ji Hu, Wenchao Li, Yatong Jiang, and Chenggang Yan.  Image classification base on pca of multi-view deep representation.  Journal of Visual Communication and Image Representation, 62:253–258, 2019.
  * [22] Joel A Tropp.  An introduction to matrix concentration inequalities.  arXiv preprint arXiv:1501.01571, 2015.
  * [23] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027, 2010.
  * [24] Vincent Q Vu, Juhee Cho, Jing Lei, and Karl Rohe.  Fantope projection and selection: A near-optimal convex relaxation of sparse pca.  Advances in neural information processing systems, 26, 2013.
  * [25] Zhanyu Wang and Jean Honorio.  The sample complexity of meta sparse regression.  In International Conference on Artificial Intelligence and Statistics, pages 2323–2331. PMLR, 2021.
  * [26] Qian Zhang, Yilin Zheng, and Jean Honorio.  Meta learning for support recovery in high-dimensional precision matrix estimation.  In International Conference on Machine Learning, pages 12642–12652. PMLR, 2021.
  * [27] Taiping Zhang, William S Chandler, James M Hoell, David Westberg, Charles H Whitlock, and Paul W Stackhouse.  A global perspective on renewable energy resources: Nasa’s prediction of worldwide energy resources (power) project.  In Proceedings of ISES World Congress 2007 (Vol. I–Vol. V), pages 2636–2640. Springer, 2008.
  * [28] Peng Zhao and Bin Yu.  On model selection consistency of lasso.  The Journal of Machine Learning Research, 7:2541–2563, 2006.
  * [29] Hao Zheng, Zhanlei Yang, Liwei Qiao, Jianping Li, and Wenju Liu.  Attribute knowledge integration for speech recognition based on multi-task learning neural networks.  In Sixteenth Annual Conference of the International Speech Communication Association, 2015.
  * [30] Hui Zou, Trevor Hastie, and Robert Tibshirani.  Sparse principal component analysis.  Journal of computational and graphical statistics, 15(2):265–286, 2006.

##  Appendix A Technical Desiderata

In this section of the appendix, we state some technical definitions and
lemmas. We begin with the definition of a real valued sub-Gaussian random
variable [26].

######  Definition A.1.

A random variable $X$ is called sub-Gaussian with parameter $\sigma\geq 0$ if,
$\forall\rho\in\mathbb{R}$

$~{}\mathbb{E}[e^{\rho X}]\leq e^{\frac{\rho^{2}\sigma^{2}}{2}}.$ (22)

As noted in Example 5.8 [23], all bounded or Gaussian random variables are
also sub-Gaussian. We can extend the definition of univariate sub-Gaussian
distributions to define a general class of multivariate sub-Gaussian
distribution on p. Let $\Sigma\in^{p\times p}$ be a positive semi-definite,
real symmetric matrix.

######  Definition A.2.

We say that a random vector $X\in^{p}$ follows a multivariate sub-Gaussian
distribution with covariance matrix $\Sigma$ and parameter $\sigma$ if

  1. 1.

$\mathbb{E}[X]=0$, and $\mathrm{Cov}(X)=\Sigma$.

  2. 2.

For $1\leq i\leq p$, $\frac{X_{i}}{\sqrt{\Sigma_{i,i}}}$ is a sub-Gaussian
random variable with parameter $\sigma$.

Our last definition of this section is the sub-exponential distribution [19].

######  Definition A.3.

A random variable $X$ with mean $0$ is said to follow a sub-exponential
distribution with parameters $(\eta,b)$ if,

$\displaystyle\mathbb{E}[e^{\lambda X}]\leq
e^{\frac{\lambda^{2}\eta^{2}}{2}}\text{ for }|\lambda|<\frac{1}{b}.$

######  Remark A.4.

It is well established that (see for instance Equation (6) in [4]) the square
of a sub-Gaussian random variable with parameter $\sigma$ is sub-exponential
with parameters $(4\sqrt{2}\sigma^{2},4\sigma^{2})$.

Next, we recall the matrix Chernoff bound (Theorem 5.1.1 in [22]).

######  Lemma A.5.

Let $A_{i}$ be random matrices such that $\mathbb{E}[A_{i}]=0$. If
$\left\|A_{i}\right\|_{\infty,\infty}\leq\gamma$ almost surely, then for every
$\epsilon>0$,

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{m}\sum
A_{i}\right\|_{\infty,\infty}\geq\epsilon\right)\leq
2pe^{-\frac{3\epsilon^{2}m}{8\gamma^{2}}}.$

We then recall the following lemma from [7].

######  Lemma A.6.

Under the assumptions in Theorem 4.5, let $\tilde{H}$ be the solution to the
further constrained problem (15). Then $\tilde{H}$ is rank $k$ and unique.
Furthermore, there exist $|J|\times k$ orthonormal matrices $U_{J}$,
$\hat{U}_{J}$ such that:

  1. 1.

$\begin{bmatrix}U_{J}\\\ 0\end{bmatrix}$ and $\begin{bmatrix}\hat{U}_{J}\\\
0\end{bmatrix}$ span the $k$-dimensional principal subspaces of $\Sigma$ and
$S-\rho\tilde{Z}$, respectively.

  2. 2.

There exists a $|J|\times|J|$ orthonormal matrix $Q$ such that

$\displaystyle\hat{U}_{J}$ $\displaystyle=$ $\displaystyle QU_{J},$
$\displaystyle\lVert Q-\mathbb{I}_{|J|}\rVert_{F}$ $\displaystyle\leq$
$\displaystyle\frac{8\rho|J|}{\lambda_{diff}}.$

We also recall the following lemma about the principal subspace projection of
a matrix. Theorem 1 in [15].

######  Lemma A.7.

Let $A$ be a symmetric matrix with eigenvalues
$\gamma_{1}\geq\cdots\geq\gamma_{p}$ and orthonormal eigenvectors
$v_{1},\ldots,v_{p}$.

$\max_{H\in\mathcal{F}^{k}}\langle A,H\rangle=\gamma_{1}+\cdots+\gamma_{k}$
and the maximum is achieved by the projector of a $k$-dimensional principal
subspace of $A$. Moreover, the maximiser is unique if and only if
$\gamma_{k}>\gamma_{k+1}$.

We also produce the following lemma along with the proof on the upper bound of
the infinite norm of the product of two matrices.

######  Lemma A.8.

Let $A$ and $B$ be two matrices with conformable dimensions. Then,

$\displaystyle\left\|AB\right\|_{\infty,\infty}\leq\left\|A\right\|_{\infty,\infty}\|B\|_{1,\infty}$

###### Proof.

$\displaystyle\left\|AB\right\|_{\infty,\infty}=\max_{i,j\in[p]}|\sum_{k\in[p]}A_{i,k}B_{k,j}|\leq\max_{i,j,k_{0}\in[p]}|A_{i,k_{0}}||\sum_{k\in[p]}B_{k,j}|\leq\left\|A\right\|_{\infty,\infty}\|B\|_{1,\infty}.$

∎

##  Appendix B Proofs of Main Theorems

###  B.1 Proof of Theorem 4.5

###### Proof.

Continue from the setup in Section 4.1. We now have to complete the proof by
verifying the remaining 4 steps.

STEP 1. We begin by finding a probabilistic upper bound on
$\|S-\Sigma\|_{\infty,\infty}$. As seen in Equation 3,
$S=\frac{1}{m}\sum_{i=1}^{m}S^{(i)}$, where
$S^{(i)}=\frac{1}{n}\sum_{j=1}^{n}X_{j}^{(i)}(X_{j}^{(i)})^{T}$ are estimated
variance matrix from auxiliary tasks.

$\displaystyle\left\|S-\Sigma\right\|_{\infty,\infty}$
$\displaystyle\leq\left\|\frac{1}{m}\sum_{i=1}^{m}\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}+\left\|\frac{1}{m}\sum_{i=1}^{m}\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}$
(23) $\displaystyle\text{Term 1}+\text{Term 2}.$

TERM 2. We begin by analysing the second term in Equation 23. We note that for
any fixed $i$,

$\displaystyle\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}$
$\displaystyle=\left\|R^{(i)}U(\Lambda+D^{(i)})(R^{(i)}U)^{T}-U\Lambda
U^{T}\right\|_{\infty,\infty}$
$\displaystyle\leq\left\|R^{(i)}U(\Lambda+D^{(i)})(R^{(i)}U)^{T}\right\|_{\infty,\infty}+\left\|U\Lambda
U^{T}\right\|_{\infty,\infty}$
$\displaystyle=\left\|(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})U(\Lambda+D^{(i)})U^{T}(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})^{T}\right\|_{\infty,\infty}+\left\|U\Lambda
U^{T}\right\|_{\infty,\infty}$
$\displaystyle=\left\|(R^{(i)}-\mathbb{I}_{p})U(\Lambda+D^{(i)})U^{T}(R^{(i)}-\mathbb{I}_{p})^{T}\right\|_{\infty,\infty}+\left\|U(\Lambda+D^{(i)})U^{T}(R^{(i)}-\mathbb{I}_{p})^{T}\right\|_{\infty,\infty}$
$\displaystyle\qquad+\left\|(R^{(i)}-\mathbb{I}_{p})U(\Lambda+D^{(i)})U^{T}\right\|_{\infty,\infty}+\left\|U\Lambda
U^{T}\right\|_{\infty,\infty}.$

Applying Lemma A.8 we get,

$\displaystyle\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}$
$\displaystyle\leq\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}^{2}\left\|U(\Lambda+D^{(i)})U^{T}\right\|_{\infty,\infty}+2\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}\left\|U(\Lambda+D^{(i)})U^{T}\right\|_{\infty,\infty}$
$\displaystyle\qquad+\left\|U\Lambda U^{T}\right\|_{\infty,\infty},$

which can be upper bounded by,

$\displaystyle\left(\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}+1\right)^{2}\left(\left\|U\Lambda
U^{T}\right\|_{\infty,\infty}+\left\|UD^{(i)}U^{T}\right\|_{\infty,\infty}\right).~{}$
(24)

Since $U$ is an orthogonal matrix,
$\left\|UD^{(i)}D^{T}\right\|_{\infty,\infty}\leq\lambda_{1}(D^{(i)})$. Then,
from 4.2, we get,

$\displaystyle\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}\leq\left(\mathbb{C}_{r}+1\right)^{2}\left(\lambda_{1}(\Sigma)+L\right).$

It is easy to see, by an application of triangle inequality and Jensen’s
inequality that,

$\displaystyle\left\|\Sigma^{(i)}-\Sigma-\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}\leq
2\left(\mathbb{C}_{r}+1\right)^{2}\left(\lambda_{1}(\Sigma)+L\right).$

To upper bound
$\left\|\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}$, we follow a
similar decomposition.

$\displaystyle\left\|\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}$
$\displaystyle=\left\|\mathbb{E}[R^{(i)}U(\Lambda+D^{(i)})(R^{(i)}U)^{T}]-U\Lambda
U^{T}\right\|_{\infty,\infty}$
$\displaystyle=\left\|\mathbb{E}[(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})U(\Lambda+D^{(i)})U^{T}(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})^{T}]-U\Lambda
U^{T}\right\|_{\infty,\infty}$ Since $D^{(i)}$ is independent of $R^{(i)}$, we
can decompose the expectation of products into the product of expectations.
Consequently, the terms $D^{(i)}$ vanishes, and we are left with the following
equality,
$\displaystyle\left\|\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}$
$\displaystyle=\left\|\mathbb{E}[(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})U\Lambda
U^{T}(R^{(i)}-\mathbb{I}_{p}+\mathbb{I}_{p})^{T}]-U\Lambda
U^{T}\right\|_{\infty,\infty}$
$\displaystyle=\left\|\mathbb{E}\left[(R^{(i)}-\mathbb{I}_{p})U\Lambda
U^{T}(R^{(i)}-\mathbb{I}_{p})^{T}+U\Lambda
U^{T}(R^{(i)}-\mathbb{I}_{p})^{T}+(R^{(i)}-\mathbb{I}_{p})U\Lambda
U^{T}\right]\right\|_{\infty,\infty}.$

Repeatedly applying Lemma A.8 and Jensen’s inequality, we get

$\displaystyle\left\|\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}$
$\displaystyle\leq\mathbb{E}[\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}]^{2}\left\|U\Lambda
U^{T}\right\|_{\infty,\infty}+2\mathbb{E}[\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}]\left\|U\Lambda
U^{T}\right\|_{\infty,\infty},$ which, upon applying Chebyshev’s inequality on
$\mathbb{E}[\|R^{(i)}-\mathbb{I}_{p}\|_{1,\infty}]$ yields,
$\displaystyle\left\|\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}$
$\displaystyle\leq\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma).$
(25)

Now, putting $A_{i}=\Sigma^{(i)}-\Sigma-\mathbb{E}[\Sigma^{(i)}-\Sigma]$
applying Lemma A.5, we get that for any $\epsilon>0$

$\displaystyle\mathbb{P}\left(\frac{1}{m}\sum_{i=1}^{m}\left\|\Sigma^{(i)}-\Sigma-\mathbb{E}[\Sigma^{(i)}-\Sigma]\right\|_{\infty,\infty}>\epsilon\right)\leq
2pe^{-\frac{3\epsilon^{2}m}{8\left(\lambda^{\dagger}\right)^{2}}}.$

Consequently, from Equation 25 and an application of triangle inequality we
get,

$\displaystyle\mathbb{P}\left(\frac{1}{m}\sum_{i=1}^{m}\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}>\epsilon\right)\leq
2pe^{-\frac{3m\left(\epsilon-\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)-2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma)\right)^{2}}{8\left(\lambda^{\dagger}\right)^{2}}}$

Let

$\displaystyle\epsilon=4\lambda^{\dagger}\sqrt{\frac{\log(p)}{m}}+\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma).$

Then, with probability at most $\frac{2}{p^{2}}$,

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\left\|\Sigma^{(i)}-\Sigma\right\|_{\infty,\infty}\geq
4\lambda^{\dagger}\sqrt{\frac{\log(p)}{m}}+\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma).$

TERM 1. In this part of the proof, we provide an upper bound for the first
term from Equation 23. We begin by finding an upper bound to
$\left\|\Sigma^{(i)}\right\|_{\infty,\infty}$. Similarly, as to the derivation
of Equation 24 we can also write,

$\displaystyle\left\|\Sigma^{(i)}\right\|_{\infty,\infty}\leq\lambda^{\dagger}.$

By using Lemma 4.1 we get,

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}\geq\epsilon\right)\leq\frac{p(p+1)}{2}e^{-\frac{nm\epsilon^{2}}{512\sigma^{4}\lambda^{\dagger}}}.$

putting
$\epsilon=16\sqrt{2\sigma^{4}\lambda^{\dagger}\frac{\log\left(p+1\right)}{nm}}$,
we get with probability $1-\frac{1}{2(p+1)^{2}}$,

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}\leq
16\sqrt{2\sigma^{4}\lambda^{\dagger}\frac{\log\left(p+1\right)}{nm}}\right).$

Thus,

$\displaystyle\mathbb{P}\left(\left\|S-\Sigma\right\|_{\infty,\infty}>\max\left\\{4\lambda^{\dagger}\sqrt{\frac{\log(p)}{m}}+\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma),16\sqrt{2\sigma^{4}\lambda^{\dagger}\frac{\log\left(p+1\right)}{nm}}\right\\}\right)$
$\displaystyle\qquad\leq\frac{2}{p^{2}}+\frac{1}{2(p+1)^{2}}$

Since,

$\frac{8|J|}{\lambda_{k}-\lambda_{k+1}}\lVert\Sigma_{J^{c},J}\rVert_{2,\infty}<1,$

there exists a constant $\alpha\in(0,1)$ such that if,

$\rho>\alpha\max\left\\{4\lambda^{\dagger}\sqrt{\frac{\log(p)}{m}}+\frac{\mathbb{C}_{R}}{p}\lambda_{1}(\Sigma)+2\sqrt{\frac{\mathbb{C}_{R}}{p}}\lambda_{1}(\Sigma),16\sqrt{2\sigma^{4}\lambda^{\dagger}\frac{\log\left(p+1\right)}{nm}}\right\\},$

then with probability at least $1-\frac{2}{p^{2}}-\frac{1}{2(p+1)^{2}}$,

$\displaystyle\rho^{-1}\lVert
S-\Sigma\rVert_{\infty,\infty}+\frac{8|J|}{\lambda_{diff}}\lVert\Sigma_{J^{c},J}\rVert_{2,\infty}\leq
1.$ (26)

It only remains to choose an appropriate matrix $Q$. By the application of
Lemma A.6 there exists a $Q$ matrix such that

$\lVert Q-\mathbb{I}_{|J|}\rVert_{F}\leq\frac{8\rho|J|}{\lambda_{diff}},$

and

$U_{J}=Q\hat{U}_{J},$

where $U_{J}$ and $\hat{U}_{J}$ are defined as in Lemma A.6. We select this
$Q$ for the remainder of the proof.

STEP 2. In this section of the proof, we show that our proposed estimator is
feasible, unique, and optimal, i.e., satisfies the KKT conditions (12,13,14).

FEASIBLITY. The feasibility of $\hat{H}$ is obvious from construction. To
check the feasibility of $\hat{Z}$, we need to verify that for a given choice
of $\rho$,

$\hat{Z}_{ij}\in[-1,1]\forall(i,j)\in(J\times J)^{c}.$

Therefore, it is sufficient to verify

$\displaystyle\frac{1}{\rho}\bigl{[}|S_{ij}-\Sigma_{ij}|+\bigl{|}\Sigma_{ij}-\langle
Q_{i*},\Sigma_{J,j}\rangle\bigr{|}\bigr{]}\leq 1.$

By an application of Cauchy-Schwarz inequality, we get

$\displaystyle\frac{1}{\rho}\bigl{[}|S_{ij}-\Sigma_{ij}|+\bigl{|}\Sigma_{ij}-\langle
Q_{i*},\Sigma_{J,j}\rangle\bigr{|}\bigr{]}$
$\displaystyle\leq\frac{1}{\rho}\bigl{[}|\left\|S-\Sigma\right\|_{\infty,\infty}+\|Q-\mathbb{I}_{|J|}\|_{F}\|\Sigma\|_{2,\infty}\bigr{]}$
$\displaystyle\leq\frac{1}{\rho}\left[\left\|S-\Sigma\right\|_{\infty,\infty}+\frac{8\rho|J|}{\lambda_{diff}}\|\Sigma\|_{2,\infty}\right].$

Thus, with probability at least $1-\frac{2}{p^{2}}-\frac{1}{2(p+1)^{2}}$,

$\displaystyle\rho^{-1}\lVert
S-\Sigma\rVert_{\infty,\infty}+\frac{8|J|}{\lambda_{diff}}\lVert\Sigma_{J^{c},J}\rVert_{2,\infty}\leq
1.$ (27)

OPTIMALITY. $\hat{H}$ only has non-zero elements in the sub-matrix
$\hat{H}_{J,J}$. By construction
$(\hat{H}_{i,j},\hat{Z}_{i,j})=(\tilde{H}_{i,j},\tilde{Z}_{i,j})\forall(i,j)\in
J\times J$. From the optimality of $(\tilde{H},\tilde{Z})$, Equation 12 is
satisfied.

The fact that Equation 13 is satisfied follows from the optimality of
$(\tilde{H},\tilde{Z})$ when $(i,j)\in J\times J$ and feasibility of
$\tilde{Z}$ otherwise.

Note that $\hat{U}_{J}$ spans the $k$-dimensional principal subspace of
$S-\rho\tilde{Z}$. We now show that $\hat{U}$ spans one (not necessarily
principal) $k$-dimensional subspace of $S-\rho\hat{Z}$. In other words, it is
enough to show that, for some diagonal matrix $D_{0}$,

$(S-\rho\hat{Z})\begin{bmatrix}\hat{U}_{J}\\\
0\end{bmatrix}=D\begin{bmatrix}\hat{U}_{J}\\\ 0\end{bmatrix}.~{}$

Analysing the term on the left-hand side of the previous equation, we get,

$\displaystyle(S-\rho\widehat{Z})\begin{bmatrix}\hat{U}_{J}\\\ 0\end{bmatrix}$
$\displaystyle=\begin{bmatrix}S_{J,J}-\rho\tilde{Z}_{J,J}&Q\Sigma_{J,J^{c}}\\\
\Sigma_{J^{c},J}Q^{T}&\Sigma_{J^{c},J^{c}}+\mathbb{D}(S_{J^{c},J^{c}}-\Sigma_{J^{c},J^{c}})\end{bmatrix}\begin{bmatrix}\hat{U}_{J}\\\
0\end{bmatrix}$
$\displaystyle=\begin{bmatrix}(S_{J,J}-\rho\tilde{Z}_{J,J})\hat{U}_{J}\\\
\Sigma_{J^{c}J}Q^{T}\hat{U}_{J}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}(S_{J,J}-\rho\tilde{Z}_{J,J})\hat{U}_{J}\\\
\Sigma_{J^{c}J}Q^{T}QU_{J}\end{bmatrix}$

Recall from Lemma A.6 that there exists a diagonal matrix $D_{1}$ such that
$(S_{J,J}-\rho\tilde{Z}_{J,J})\hat{U}_{J}=D_{1}\hat{U}_{J}$. Also, the fact
that $\begin{bmatrix}U_{J}\\\ 0\end{bmatrix}$ spans $k$-dimensional principal
subspace of $\Sigma$ implies that $\Sigma_{J^{c}J}U_{J}=0$. Therefore,

$\displaystyle\begin{bmatrix}(S_{J,J}-\rho\tilde{Z}_{J,J})\hat{U}_{J}\\\
\Sigma_{J^{c}J}Q^{T}QU_{J}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}D_{1}\hat{U}_{J}\\\
\Sigma_{J^{c}J}U_{J}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}D_{1}\hat{U}_{J}\\\ 0\end{bmatrix}$
$\displaystyle=\begin{bmatrix}D_{1}\hat{U}_{J}\\\ 0\end{bmatrix}$
$\displaystyle=\begin{bmatrix}D_{1}&0\\\
0&0\end{bmatrix}\begin{bmatrix}\hat{U}_{J}\\\ 0\end{bmatrix}$

Setting

$D_{0}=\begin{bmatrix}D_{1}&0\\\ 0&0\end{bmatrix}$

we successfully show that Section B.1 holds. We will further prove that
$diag(D_{1})$ are the first $k$ eigenvalues of $(S-\rho\widehat{Z})$, and
$\lambda_{k}(S-\rho\widehat{Z})-\lambda_{k+1}(S-\rho\widehat{Z})>0$. Then
uniqueness will follow by Lemma A.6.

$\displaystyle(S-\rho\hat{Z})$
$\displaystyle=\begin{bmatrix}S_{JJ}-\rho\tilde{Z}_{JJ}-Q\Sigma_{JJ}Q^{T}&0\cr
0&\Sigma_{J^{c},J^{c}}+\mathbb{D}(S_{J^{c},J^{c}}-\Sigma_{J^{c},J^{c}})\end{bmatrix}+\begin{bmatrix}Q\Sigma_{JJ}Q^{T}&Q\Sigma_{JJ^{c}}\vspace*{3pt}\cr\Sigma_{J^{c}J}Q^{T}&\Sigma_{J^{c}J^{c}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}S_{JJ}-\rho\tilde{Z}_{JJ}-Q\Sigma_{JJ}Q^{T}&0\cr
0&\Sigma_{J^{c},J^{c}}+\mathbb{D}(S_{J^{c},J^{c}}-\Sigma_{J^{c},J^{c}})\end{bmatrix}+\begin{bmatrix}Q&0\cr
0&I\end{bmatrix}\times\Sigma\times\begin{bmatrix}Q^{T}&0\cr 0&I\end{bmatrix}$

However, since we just proved that $(S-\rho\widehat{Z})$ satisfies the
conditions of Proposition 5.1, the problem reduces to finding an upper bound
to $\lambda_{l}(S_{J,J}-\rho\tilde{Z}_{J,J}-Q\Sigma_{J,J}Q^{T})$ for each
$l\in[|J|]$. Note that,

$\displaystyle\lambda_{l}(S_{J,J}-\rho\tilde{Z}_{J,J}-Q\Sigma_{J,J}Q^{T})$
$\displaystyle\leq\lambda_{l}(S_{J,J}-\rho\tilde{Z}_{J,J}-\Sigma_{J,J})+\lambda_{l}(\Sigma_{J,J}-Q\Sigma_{J,J}Q^{T})$
$\displaystyle\leq\lVert
S_{J,J}-\rho\tilde{Z}_{J,J}-\Sigma_{J,J}\rVert_{F}+\lVert\Sigma_{J,J}-Q\Sigma_{J,J}Q^{T}\rVert_{F}$
$\displaystyle\leq\lVert
S_{J,J}-\rho\tilde{Z}_{J,J}-\Sigma_{J,J}\rVert_{F}+2\lVert\Sigma_{J,J}\rVert_{F}\times\lVert
Q-I\rVert_{F}$ $\displaystyle\leq 2\rho s+2\lambda_{1}(\Sigma)\times
8\rho|J|/(\lambda_{k}-\lambda_{k+1}).$

Where the last inequality follows from the combination of part 1 and Lemma
A.6. The rest of the proof since for a large enough $\rho$ we have,

$\displaystyle
4\rho|J|+\frac{32\lambda_{1}(\Sigma)\rho|J|}{(\lambda_{k}-\lambda_{k+1})}\leq\lambda_{k}-\lambda_{k+1}.$

This completes the proof of step 2.

STEP 3. In this section of the proof, we prove that
$\left\|\hat{H}-\Pi\right\|_{\infty,\infty}<2\frac{\rho|J|}{\lambda_{diff}}$.
By step 2, we have,

$\displaystyle\left\|\hat{H}-\Pi\right\|_{\infty,\infty}$
$\displaystyle=\left\|\begin{bmatrix}\hat{U}_{J}\hat{U}_{J}^{T}&0\\\
0&0\end{bmatrix}-\begin{bmatrix}U_{J}U_{J}^{T}&0\\\
0&0\end{bmatrix}\right\|_{\infty,\infty}$
$\displaystyle=\left\|\begin{bmatrix}QU_{J}U_{J}^{T}Q^{T}&0\\\
0&0\end{bmatrix}-\begin{bmatrix}U_{J}U_{J}^{T}&0\\\
0&0\end{bmatrix}\right\|_{\infty,\infty}$
$\displaystyle=\left\|QU_{J}U_{J}^{T}Q^{T}-U_{J}U_{J}^{T}\right\|_{\infty,\infty}$
$\displaystyle\leq\lambda_{1}(U_{J}U_{J}^{T})\|\mathbb{I}-Q\|_{F}$
$\displaystyle\leq\frac{8|J|\rho}{\lambda_{diff}}.$

Setting $\mathbb{C}_{\pi}=\frac{8|J|}{\lambda_{diff}}$, completes the proof of
the upper bound for $\left\|\hat{H}-\Pi\right\|_{\infty,\infty}$.

STEP 4. Set

$\rho<\min_{i\in J}\frac{\Pi_{i,i}\lambda_{diff}}{16|J|}.$

The rest of the proof now follows from the step 3. ∎

###  B.2 Proof of Proposition 5.1

###### Proof.

Let $\Pi^{(m+1)}$ be the $k$-th dimensional principal component corresponding
to $\Sigma^{(m+1)}$. Hence,

$\displaystyle\Pi^{(m+1)}=\sum_{i=1}^{k}U_{*,i}^{(m+1)}\left(U_{*,i}^{(m+1)}\right)^{T}.$

We also know from the recovered support union that $\Pi_{i,i}^{(m+1)}=0$
whenever $i\notin J$. However,

$\displaystyle\Pi_{i,i}^{(m+1)}=\sum_{j=1}^{k}\left(U_{j,i}^{(m+1)}\right)^{2}.$

Hence, $\forall\enspace j\in J^{c}$, $u_{i,j}^{(m+1)}=0$. Therefore,

$U_{i}^{(m+1)}=\begin{pmatrix}U_{J,i}^{(m+1)}\\\ 0\end{pmatrix}.$

Moreover, we also know that

$\displaystyle\Sigma^{(m+1)}U_{i}^{(m+1)}=\lambda_{i}U_{i}^{(m+1)}.$

In other words,

$\displaystyle\begin{bmatrix}\Sigma_{J,J}^{(m+1)}&\Sigma_{J,J^{c}}^{(m+1)}\\\
\Sigma_{J^{c},J}^{(m+1)}&\Sigma_{J^{c},J^{c}}^{(m+1)}\end{bmatrix}\begin{pmatrix}U_{J,i}^{(m+1)}\\\
0\end{pmatrix}=\lambda_{i}\begin{pmatrix}U_{J,i}^{(m+1)}\\\ 0\end{pmatrix}.$

Which in turn implies,

$\displaystyle\Sigma_{J,J}^{(m+1)}U_{J,i}^{(m+1)}=\lambda_{i}U_{J,i}^{(m+1)}.$

This completes the proof. ∎

###  B.3 Proof of Theorem 5.3

###### Proof.

As given, $\Sigma^{(m+1)}$ is the novel covariance matrix with the following
block structure,

$\displaystyle\Sigma^{(m+1)}=\begin{bmatrix}\Sigma_{J,J}^{(m+1)}&\Sigma_{J,J^{c}}^{(m+1)}\\\
\Sigma_{J^{c},J}^{(m+1)}&\Sigma_{J^{c},J^{c}}^{(m+1)}\end{bmatrix}.$

For $1\leq i\leq k$, let $U_{i}^{(m+1)}$ be the eigenvectors corresponding to
$\Sigma^{(m+1)}$, and $\lambda_{i}^{(m+1)}$ be the corresponding eigenvalues.
Let $J$ be the recovered support union. $U_{J,i}^{(m+1)}$ denotes those
coordinates of $U_{i}^{(m+1)}$ that belong to set $J$. Similarly to the proof
of Theorem 4.5, if the penalty parameter $\rho$ satisfies,

$\displaystyle\rho^{-1}\lVert
S_{J,J}^{(m+1)}-\Sigma_{J,J}^{(m+1)}\rVert_{\infty,\infty}+\frac{8|J|^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\lVert\Sigma_{(J^{(m+1})^{c}J^{(m+1)}}\rVert_{2,\infty}\leq
1\quad\text{and }$ $\displaystyle
4\rho|J^{(m+1)}|\biggl{(}1+\frac{8\lambda_{1}^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\biggr{)}<\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)},$

then the solution $\widehat{H}^{(m+1)}$ of Equation 20 is unique and
satisfies${supp}(\widehat{H}^{(m+1)})\subseteq J^{(m+1)}$. It is noteworthy
that if for each $j\in[n^{(m+1)}]$, $X^{(m+1)}_{j}$ is a $p$-dimensional
multivariate sub-Gaussian distribution with covariance matrix $\Sigma$ and
parameter $\sigma$, then for any set of sub-indices $J$, $X_{j,J}^{(m+1)}$ is
a $|J|$-dimensional multivariate sub-Gaussian distribution with covariance
matrix $\Sigma_{J,J}$ and parameter $\sigma$. Following from Equation 19, we
can also see that by construction

$S_{J,J}^{(m+1)}=\frac{1}{n^{(m+1)}}\sum_{l=1}^{n^{(m+1)}}X_{l,J}^{(m+1)}\left(X_{l,J}^{(m+1)}\right)^{T}.$

By using Lemma 4.1, we know that as long as $\epsilon\leq
32\lambda_{1}\left(\Sigma^{(m+1)}\right)\sigma^{2}$,

$\displaystyle\mathbb{P}\left(\left\|\left(S^{(m+1)}-\Sigma^{(m+1)}\right)\right\|_{\infty,\infty}\geq\epsilon\right)\leq\frac{|J|(|J|+1)}{2}e^{-\frac{n\epsilon^{2}}{512\sigma^{4}\lambda_{1}\left(\Sigma^{(m+1)}\right)}}.$

Define
$\alpha^{(m+1)}:=1-\frac{8|J|^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\lVert\Sigma_{(J^{(m+1})^{c}J^{(m+1)}}\rVert_{2,\infty}$.
By 5.2, $\alpha^{(m+1)}\in(0,1)$. Therefore, as long as,

$\displaystyle
n^{\frac{1}{3}}>\max\left\\{\frac{1}{32\lambda_{1}\left(\Sigma^{(m+1)}\sigma^{2}\right)},\left(\frac{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}{4|J^{(m+1)}|\left(1+\frac{8\lambda_{1}^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\right)}\right)^{-1}\right\\}\text{
and,}$
$\displaystyle\frac{1}{n^{1/3}\alpha^{(m+1)}}<\rho<\frac{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}{4|J^{(m+1)}|\left(1+\frac{8\lambda_{1}^{(m+1)}}{\lambda_{k}^{(m+1)}-\lambda_{k+1}^{(m+1)}}\right)},$

with probability
$1-\frac{|J|(|J|+1)}{2}e^{-\frac{n^{\frac{1}{3}}}{512\sigma^{4}\lambda_{1}\left(\Sigma^{(m+1)}\right)}}$,
we have $supp(\Pi^{(m+1)})=J^{(m+1)}$. ∎

###  B.4 Proof of Lemma 4.1

###### Proof.

We note that $S^{(i)}-\Sigma^{(i)}$ are i.i.d with
$\mathbb{E}\left[S^{(i)}-\Sigma^{(i)}\right]=0$. We begin by finding an
entrywise bound on the elements of
$\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)$. Since
$S^{(i)}=\frac{1}{p}\sum X_{p}^{(i)}\left(X_{p}^{(i)}\right)^{T}$,

$\displaystyle\mathbb{P}\left(\left|\frac{1}{mn}\sum_{i,p}X_{k,p}^{(i)}X_{l,p}^{(i)}-\sigma_{k,l}^{(i)}\right|\geq\epsilon\right).$

Let $\sigma_{k}:=\max_{1\leq i\leq K}{\sigma}_{k,k}^{(i)}$,
$\tilde{X}_{k,p}^{(i)}:=\frac{X_{k,p}^{(i)}}{\sqrt{\sigma_{k}}}$,
$\tilde{\rho}_{ij}^{(k)}:=\frac{{\sigma}_{k,l}^{(i)}}{\sqrt{\sigma_{k}\sigma_{l}}}$.
We have

$\mathbb{P}\left(\left|\frac{1}{mn}\sum_{i,p}X_{k,p}^{(i)}X_{l,p}^{(i)}-\sigma_{k,l}^{(i)}\right|\geq\epsilon\right)=\mathbb{P}\left[4\left|\sum_{i,p}\left(\tilde{X}_{k,p}^{(i)}\tilde{X}_{l,p}^{(i)}-\tilde{\rho}_{k,l}^{(i)}\right)\right|>\frac{4nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right]$

Define $U^{(i)}_{p,k,l}:=\tilde{X}_{p,k}^{(i)}+\tilde{X}_{p,l}^{(i)}$,
$V^{(i)}_{p,k,l}:=\tilde{X}_{p,k}^{(i)}-\tilde{X}_{p,l}^{(i)}$. Then for any
$r\in\mathbb{R}$,

$4\sum_{i,p}\left(\tilde{X}_{k,p}^{(i)}\tilde{X}_{l,p}^{(i)}-\tilde{\rho}_{k,l}^{(i)}\right)=\sum_{k,t}\left\\{\left(U^{(i)}_{p,k,l}\right)^{2}-2\left(r+\tilde{\rho}^{(i)}_{k,l}\right)\right\\}-\sum_{k,t}\left\\{\left(V^{(i)}_{p,k,l}\right)^{2}-2\left(r-\tilde{\rho}^{(i)}_{k,l}\right)\right\\}$
(29)

Thus,

$\displaystyle\mathbb{P}\left(\left|\frac{1}{mn}\sum_{i,p}X_{k,p}^{(i)}X_{l,p}^{(i)}-\sigma_{k,l}^{(i)}\right|\geq\epsilon\right)\leq$
$\displaystyle\mathbb{P}\left[\Big{|}\sum_{k,t}\left\\{\left(U^{(k)}_{t,ij}\right)^{2}-2\left(r+\tilde{\rho}^{(k)}_{ij}\right)\right\\}\Big{|}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right]$
(30)
$\displaystyle+\mathbb{P}\left[\Big{|}\sum_{k,t}\left\\{\left(V^{(k)}_{t,ij}\right)^{2}-2\left(r-\tilde{\rho}^{(k)}_{ij}\right)\right\\}\Big{|}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right].$

Consider the first term

$\displaystyle\mathbb{P}\left[\Big{|}\sum_{i,p}\left\\{\left(U^{(i)}_{p,k,l}\right)^{2}-2\left(r+\tilde{\rho}^{(k)}_{ij}\right)\right\\}\Big{|}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right],$

and let
$r=\frac{\sigma_{k,k}^{(i)}}{\sigma_{k}}+\frac{\sigma_{l,l}^{(i)}}{\sigma_{l}}$.
Then,

$\mathbb{E}\left[\left(U^{(i)}_{p,k,l}\right)^{2}-2\left(r+\tilde{\rho}^{(k)}_{ij}\right)\right]=0.$

Let
$Z_{p,k,l}^{(i)}:=\left(U^{(i)}_{p,k,l}\right)^{2}-2\left(r+\tilde{\rho}^{(k)}_{ij}\right)$.
Consequently,

$\displaystyle\mathbb{P}\left(\Big{|}\sum_{i,p}\left\\{\left(U^{(i)}_{p,k,l}\right)^{2}-2\left(r+\tilde{\rho}^{(k)}_{ij}\right)\right\\}\Big{|}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$
$\displaystyle=\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)~{}$
(31)
$\displaystyle\qquad+\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}<-\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$

To see that $U_{p,k,l}^{(i)}$ is sub-Gaussian, we note that,

$\displaystyle\mathbb{E}(e^{\lambda U_{p,k,l}^{(i)}})$
$\displaystyle=\mathbb{E}\left(e^{\lambda\left(X_{p,k}^{(i)}+X_{p,l}^{(i)}\right)}\right)$
$\displaystyle\leq\left[\mathbb{E}e^{2\lambda\tilde{X}_{p,k}^{(i)}}\right]^{0.5}\left[\mathbb{E}e^{2\lambda\tilde{X}_{p,l}^{(i)}}\right]^{0.5}$
$\displaystyle\leq
e^{\lambda^{2}\left(\frac{\sigma_{k,k}^{(i)}}{\sigma_{k}}\sigma\right)^{2}}e^{\lambda^{2}\left(\frac{\sigma_{l,l}^{(i)}}{\sigma_{l}}\sigma\right)^{2}}.$
Since $\frac{\sigma_{k,k}^{(i)}}{\sigma_{k}}$ and
$\frac{\sigma_{l,l}^{(i)}}{\sigma_{l}}$ are upper bounded by 1, we
consequently have, $\displaystyle\mathbb{E}(e^{\lambda U_{p,k,l}^{(i)}})$
$\displaystyle\leq e^{2\lambda^{2}\sigma^{2}},$

which shows that $U_{p,k,l}^{(i)}$ is sub-Gaussian with parameter $2\sigma$.
We also know, from [23, Section 5.2.4], that if $U_{p,k,l}^{(i)}$ is a sub-
Gaussian random variable with then $\left(U_{p,k,l}^{(i)}\right)^{2}$ is a
sub-exponential random variable. As noted in Section 2 if $X$ is a sub-
Gaussian with parameter $\omega$, then, $X^{2}$ is a sub-exponential
distribution with parameters $(4\sqrt{2}\omega^{2},4\omega^{2})$.
Consequently, if $|t|<\frac{1}{4\omega^{2}}$,

$\displaystyle\mathbb{E}[e^{tX^{2}}]\leq e^{16\omega^{4}t^{2}}.$

Returning to Equation 31, we observe that each $Z_{p,k,l}^{(i)}$ is i.i.d.
subexponential with parameter $(16\sqrt{2}\sigma^{2},16\sigma^{2})$. Hence for
any $\frac{1}{16\sigma^{2}}>\theta>0$,

$\displaystyle\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$
$\displaystyle=\mathbb{P}\left(\theta\sum_{i,p}Z_{p,k,l}^{(i)}>\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$
$\displaystyle=\mathbb{P}\left(e^{\theta\sum_{i,p}Z_{p,k,l}^{(i)}}>e^{\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}\right).$

Thus, by Markov’s inequality, it follows that,

$\displaystyle\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$
$\displaystyle\leq\mathbb{E}\left[e^{\theta\sum_{i,p}Z_{p,k,l}^{(i)}}\right]e^{-\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}$
$\displaystyle=\prod_{i,p}\mathbb{E}\left[e^{\theta
Z_{p,k,l}^{(i)}}\right]e^{-\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}$
$\displaystyle\leq\prod_{i,p}e^{512\theta^{2}\sigma^{4}}e^{-\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}$
$\displaystyle=e^{512nm\theta^{2}\sigma^{4}}e^{-\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}$
$\displaystyle=e^{512mn\theta^{2}\sigma^{4}-\theta\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}}$

Minimising over $\theta>0$, we get the derivative of the exponent with respect
to $\theta$ to be equal to $0$. In other words,

$\displaystyle
512mn\theta\sigma^{4}-\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}=0,$

which in turn implies that whenever

$\frac{\epsilon}{512\sigma^{4}\sqrt{\sigma_{k}\sigma_{l}}}<\frac{1}{16\sigma^{2}},$

the minimum is achieved when
$\theta=\frac{\epsilon}{512\sigma^{4}\sqrt{\sigma_{k}\sigma_{l}}}$. In other
words, if

$\epsilon<32\sigma^{2}\sqrt{\sigma_{k}\sigma_{l}},$

the minimum is,
$e^{-\frac{nm\epsilon^{2}}{512\sigma^{4}\sqrt{\sigma_{k}\sigma_{l}}}}$.
Consequently, our upper bound can be further refined as,

$\displaystyle\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)\leq
e^{-\frac{nm\epsilon^{2}}{512\sqrt{\sigma_{k}\sigma_{l}}\sigma^{4}}}.$

Since for each $k\in[m]$,
$\sigma_{k}<\max_{i\in[m]}\left\|\Sigma^{(i)}\right\|_{\infty,\infty}<\eta$,
this in turn implies that the upper bound is independent of the value of
$\Sigma^{(i)}$’s. In particular, this upper bound is distribution free. Hence,
for $\epsilon<32\eta\sigma^{2}$,

$\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}>\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)\leq
e^{-\frac{nm\epsilon^{2}}{512\eta\sigma^{4}}}.$

Similarly, the upper bound to
$\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}<-\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)$,
can be achieved; and,

$\displaystyle\mathbb{P}\left(\sum_{i,p}Z_{p,k,l}^{(i)}<-\frac{2nm\epsilon}{\sqrt{\sigma_{k}\sigma_{l}}}\right)\leq
e^{-\frac{nm\epsilon^{2}}{512\eta\sigma^{4}}}.$

Hence, we have found the entrywise error bound to
$\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)$. To find an upper bound to
$\left\|\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}$,
we use a union bound.

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{m}\sum\left(S^{(i)}-\Sigma^{(i)}\right)\right\|_{\infty,\infty}>\epsilon\right)$
$\displaystyle\leq\sum_{k,l}\mathbb{P}\left(\left|\frac{1}{mn}\sum_{i,p}X_{k,p}^{(i)}X_{l,p}^{(i)}-\sigma_{k,l}^{(i)}\right|\geq\epsilon\right).$
$\displaystyle\leq\frac{p(p+1)}{2}e^{-\frac{nm\epsilon^{2}}{512\eta\sigma^{4}}}.$

This completes the proof. ∎

##  Appendix C Additional Experiments

To validate our results, we first generate random matrices $\Sigma^{(i)}$ for
each $i\in[m]$ and then generate $X^{(i)}_{j}$ for each $j\in[n]$. We now
describe the generative procedure for $\Sigma^{(i)}$. We first produce a
random orthogonal matrix $U$ where $U_{l,*}$ is supported on $[J]$ for each
$l\in[k]$. We then produce a random diagonal matrix with each entry produced
uniformly from $[0,1]$ and add $500$ to the first $k$ many diagonal elements
to produce $\Lambda$. $\Sigma=U\Lambda(U)^{T}.$ For each $i\in[m]$, we
generate randomised uniform white noise matrices $\mathcal{A}^{(i)}$
$\mathcal{A}^{(i)}_{l,j}\overset{iid}{\sim}Uniform(0,1)$. $R^{(i)}$ is
produced by adding scaled $\mathcal{A}^{(i)}$’s to $\mathbb{I}_{p}$; i.e.
$R^{(i)}=\mathcal{A}^{(i)}/p^{2}+\mathbb{I}_{p}.$ We generate $D^{(i)}$’s, as
white noise matrices where for each $l,j\in[p];l>j,D^{(i)}_{l,j}$ is generated
uniformly from $(0,1)$. This ensures 4.2 is satisfied. $\Sigma^{i}$ is then
generated as,
$\Sigma^{(i)}=R^{(i)}U(\Lambda+D^{(i)})U^{T}\left(R^{(i)}\right)^{T}.$ We
compute the empirical probability of successful support recovery
$\mathbb{P}(\hat{J}=J)$ as the number of times we obtain the exact support
recovery among 100 repetitions, divided by 100. We compute the standard
deviation as $({\mathbb{P}(\hat{J}=J)(1-\mathbb{P}(\hat{J}=J))})/{100}$.
Performance is measured by the empirical probability of correct classification
$\mathbb{P}(\hat{J}=J)$ and maximal estimation error
$\|\hat{\Pi}-\Pi\|_{\infty,\infty}$. To create error bars we use a empirical
95% confidence interval around the estimate.

#### Uniform Distribution.

For all the experiments in this sub-section, we let $p=80$, $k=6$ and $|J|=6$
and perform 100 repetitions of each setting. We now consider the setting of
different dimensions. We choose $n\in\left\\{3,5,7,9\right\\}$ and use
$\rho=\sqrt{\frac{\log(p+1)}{mn}}$ for all pairs of $(m,n)$. $T$ is defined as
above. For each $j\in[n]$, $X^{(i)}_{j}\sim
0.5\mathcal{N}(0,\Sigma^{(i)})+0.5\Delta_{i}$, where $\Delta_{i}$ is a
$p\times 1$ vector of independent $Uniform(0,1)$ random variables. Figure 2
depicts the outcome of our experiments. For different choices of $n$, the
graphs overlap each other perfectly (both $\mathbb{P}(\hat{J}=J)$ and
$\left\|\hat{\Pi}-\Pi\right\|_{\infty,\infty}$). Then, we took the recovered
support and derived the probability of support recovery for a novel task.

|  |
---|---|---
(a) Support Recovery | (b) Maximal Error | (c) Novel task

Figure 2: Left, Center: Simulations for Theorem 4.5 on the probability of
support union recovery and maximal error under various settings of $n$. Let
$\rho=\sqrt{\frac{\log(p+1)}{mn}}$. The x-axis is set by
$T:=\frac{mn}{\log(p+1)}$ varying from $\left\\{1,\dots,15\right\\}$ for
support recovery. Right: Probability of support recovery for novel task. We
let $\rho=\sqrt{\frac{\log(|J|+1)}{n}}$. The x-axis is set by
$T:=\frac{n}{\log(|J|+1)}$ varying from $\left\\{10,\dots,25\right\\}$

We perform the experiments when there are $2,3$, or $4$ extra zeros and in all
of those experiments we have found greater than $95\%$ probability of
accurately identifying the extra zeros.

#### Sub-Gaussian Distribution.

For all the experiments in this sub-section, we let $n=5$, $k=6$ and $|J|=6$
and perform 100 repetitions of each setting. We now consider the setting of
different dimensions. We choose $p\in\left\\{40,45,50,55\right\\}$. $T$ and
$\rho$ are defined as above. For each $j\in[n]$, $X^{(i)}_{j}\sim
0.5\mathcal{N}(0,\Sigma^{(i)})+0.5\Delta_{i}$, where $\Delta_{i}$ is a
$p\times 1$ vector of independent $Exponential(0,1)$ random variables. Figure
3 depicts the outcome of our experiments. For different choices of $p$, the
graphs overlap each other perfectly (both $\mathbb{P}(\hat{J}=J)$ and
$\left\|\hat{\Pi}-\Pi\right\|_{\infty,\infty}$). Then, we took the recovered
support union and derived the probability of support recovery for a novel
task. We perform the experiments when there are $2,3$, or $4$ extra zeros and
in all of those experiments we have found greater than $90\%$ probability of
accurately identifying the extra zeros.

(a) Support Union Recovery (b) Maximal Error (c) Novel Task

Figure 3: Left, Center: Simulations for Theorem 4.5 on the probability of
support union recovery and maximal error under various settings of $p$. We let
$\rho=\sqrt{\frac{\log(p+1)}{mn}}$. The x-axis is set by
$T:=\frac{mn}{\log(p+1)}$ varying from $\left\\{1,\dots,15\right\\}$ for
support union recovery. Right: Probability of support recovery for novel task.
We take $\rho=\sqrt{\frac{\log(|J|+1)}{n}}$. The x-axis is set by
$T:=\frac{n}{\log(|J|+1)}$ varying from $\left\\{10,\dots,25\right\\}$

##  Appendix D Real World Datasets

#### Brain Imaging Data.

The 1000 functional connectomes dataset is from nitrc.org and contains
resting-state functional magnetic resonance signals from 1128 subjects, 41
sites around the world, and 157 brain regions. We use a set of independent
samples to construct the true covariance matrix, and thus the true support
union. We then constructed the pooled estimator using $m$ auxiliary tasks,
where $m\in\left\\{1,2,3,5,8,10,12,15,18,20,25,30,35,41\right\\}$. Then we
compared the recovered support with the true support to calculate the mean
percentage and standard deviation of error in the recovered support. The error
bars were obtained by creating a 95% confidence interval around the mean. As
we clearly see the probability of incurring an error decreases as we increase
the number of tasks.

Figure 4: Probability of error of support recovery for brain imaging dataset.
We took an increasing number of auxiliary tasks
$m\in\left\\{1,2,3,5,8,\dots,30,35,41\right\\}$. The probability of error
decreases.

#### Cancer Genetics Data.

The cancer genome atlas program dataset is available at gdc.cancer.gov and
contains the genetic sequencing of 1576 patients with 11 different categories
of cancer. We use a set of independent samples to construct the true
covariance matrix and the true support union. We then constructed the pooled
estimator using $m$ auxiliary tasks where $m\in\\{1,\dots 10\\}$. Then we
compared the recovered support with the true support to calculate the mean
percentage and the standard deviation of the error in recovered support. The
error bars were obtained by empirically calculating the $95\%$ confidence
interval around the mean. As we can see the probability of incurring an error
decreases as we increase the number of tasks.

Figure 5: Probability of error of support recovery for cancer genetics
dataset. We took an increasing number of auxiliary tasks
$m\in\left\\{1,2,\dots,10\right\\}$. The probability of error can be seen to
be decreasing.

##  Appendix E Comparisons

###  E.1 Multiple Independent Sparse PCA.

In this section, we compare the meta-learning method against a multiple
independent learning method when applied to support union recovery under a
Gaussian distribution setting. The method of generating the dataset and meta-
learning is as described in Section 6. As before, we let $p=50$, $|J|=5$, and
fix $n=3$. The values we choose for $m$ are chosen from the set
$\mathcal{M}\in\\{2,3,4,6,7,8,10,11,12,14\\}$. We perform 16 replications and
report the mean probability and the 95% confidence interval as the error bars
around the mean. The results are plotted in Figure 7.

#### Multiple Independent Learning.

For the multiple independent learning setup there are $m$ many objective
functions being maximised simultaneously. Let $S^{(i)}$ be the covariance
matrix associated with the task $i$. Then the $i$-th objective is to maximise,

$\displaystyle\langle S^{(i)},H\rangle-\rho\|H\|_{1,1}\text{ subject to
}H\in\mathcal{F}^{k}.$ (32)

Let $\hat{H}^{(i)}$ be the maximiser of this objective function. Let
$\hat{J}^{(i)}=supp(H^{(i)})$ and take
$\hat{J}=\bigcup_{i\in\mathcal{M}}\hat{J}^{(i)}$ to be the natural estimator
for the support union. We apply $\sqrt{\frac{\log(p+1)}{n}}$ as the penalty
for each of the $m$ tasks.

#### Meta-Learning.

For the meta-learning setup, we proceed according to Section 6 and recover the
support union from the pooled covariance estimator. We apply
$\sqrt{\frac{\log(p+1)}{mn}}$ as the penalty when $m$ is the total number of
tasks.

Multiple independent learning of sparse PCA can be seen to be completely
useless whenever the number of tasks get large and fails to recover the
correct support union after a while. On the other hand, meta-learning improves
steadily with the increasing number of tasks and the results coincide with our
experiments in Section 6.

Figure 6: Comparison showing the difference in the probability of recovering
the correct support union of meta-learning versus multiple independent sparse
PCA. y-axis : Probability of support union recovery. x-axis : Number of
auxiliary tasks.

###  E.2 Comparison With Multi-Task PCA.

In this section, we compare the meta-learning method against a multi-task
learning method when applied to support union recovery under a Gaussian
distribution setting. The method of generating the dataset and meta-learning
is as described in Section 6. As before, we let $p=50$, $|J|=5$, and fix
$n=3$. The values we choose for $m$ are chosen from the set
$\mathcal{M}\in\\{2,3,4,6,7,8,10,11,12,14\\}$. We perform 16 replications and
report the mean probability and the 95% confidence interval as the error bars
around the mean. The results are plotted in Figure 7.

Figure 7: Comparison showing the difference in the probability of recovering
the correct support union of meta-learning versus multi-task learning. y-axis
: Probability of support union recovery. x-axis : Number of auxiliary tasks.

#### Multi-Task Learning.

For the multi-task learning method, we maximise the following objective
function:

$\displaystyle\sum_{i=1}^{m}\langle
S^{(i)},H^{(i)}\rangle-\rho\|H^{(1)},\dots,H^{(m)}\|_{1}\text{ subject to
$H^{(i)}\in\mathcal{F}^{k}\ \forall\ i\in{1,\dots,m}$ },$

where $S^{(i)}$ is the covariance matrix corresponding to task $i$ and
$\rho\|H^{(1)},\dots,H^{(m)}\|_{1}=\sum_{j,k}\sup_{i\in\\{1,\dots,m\\}}|H_{j,k}^{(i)}|$.
We apply $\sqrt{\frac{\log(p+1)}{mn}}$ as the penalty for the tasks. We
recover the support union by taking the union of the recovered support.

#### Meta-Learning.

For the meta-learning setup, proceed according to Section 6 and recover the
support union from the pooled covariance estimator. We apply
$\sqrt{\frac{\log(p+1)}{mn}}$ as the penalty when $m$ is the total number of
tasks.

As we can see, meta-learning outperforms multi-task learning with 100%
accuracy for large values of $m$.