latexReferences

Inferring Covariance Structure from
Multiple Data Sources via
Subspace Factor Analysis

Noirrit Kiran Chandra†, David B. Dunson‡, Jason Xu‡,∗

†Department of Mathematical Sciences,

The University of Texas at Dallas, Richardson, TX

Email<EMAIL_ADDRESS>‡Department of Statistical Science, Duke
University, Durham, NC ∗Department of Biostatistics and Bioinformatics, Duke
University, Durham, NC

Abstract

Factor analysis provides a canonical framework for imposing lower-dimensional
structure such as sparse covariance in high-dimensional data. High-dimensional
data on the same set of variables are often collected under different
conditions, for instance in reproducing studies across research groups. In
such cases, it is natural to seek to learn the shared versus condition-
specific structure. Existing hierarchical extensions of factor analysis have
been proposed, but face practical issues including identifiability problems.
To address these shortcomings, we propose a class of SUbspace Factor Analysis
(SUFA) models, which characterize variation across groups at the level of a
lower-dimensional subspace. We prove that the proposed class of SUFA models
lead to identifiability of the shared versus group-specific components of the
covariance, and study their posterior contraction properties. Taking a
Bayesian approach, these contributions are developed alongside efficient
posterior computation algorithms. Our sampler fully integrates out latent
variables, is easily parallelizable and has complexity that does not depend on
sample size. We illustrate the methods through application to integration of
multiple gene expression datasets relevant to immunology.

_Keywords_ : Data-augmented Markov chain Monte Carlo, Data Integration,
Gradient-Based Sampling, Latent Variable Models, Multi-Study Factor Analysis.

##  Introduction

With increasing calls for reproducibility and transparency in science, it has
become standard to make datasets widely available to the scientific community.
This motivates interest in aggregating different but related datasets
together. A prominent example arises in gene network analyses where it is of
interest to merge datasets from studies that consider a common set of genes.
Combining datasets increases sample size in studying covariance structure
among the genes. To conduct principled inference, it is important to account
for differences across studies, and not simply pool the data. This article
focuses on developing statistical methods for inferring common versus study-
specific covariance structures, with a particular motivation to multi-study
applications in the high-dimensional setting. These methods immediately apply
to data from a single study featuring hierarchical, multi-group structure as
well.

Consider the analysis of gene expression levels in immune cells. Integrating
data from multiple studies serves to (1) increase statistical precision in
making inferences on covariance structure between genes; (2) yield results
that are more robust to study-to-study variability and hence more
generalizable; and (3) obtain insight into shared versus study-specific
contributors to the covariance. We focus on datasets from the Immunological
Genome Project (Heng et al., 2008), comprising microarray assays as well as
bulk RNA sequencing data as shown in Figure S.7 in the supplementary
materials. Substantial similarity between the datasets can be observed albeit
with significant amount of heterogeneity. Such heterogeneity is typical and
arises due to differences between the subject populations, and variation in
data collection technologies inducing platform-specific effects in the
respective datasets.

Given our interest in the covariance structure of high-dimensional data, it is
natural to consider Bayesian sparse factor analysis (FA) models, which achieve
state-of-the-art performance in the single-study case (Fan et al., 2008;
Knowles and Ghahramani, 2011). Even when different studies collect the same
variables for closely-related populations, applying factor analysis separately
to each dataset can lead to very different inferred factor structures. There
are recent approaches extending factor analysis to the multi-study
setting—notably, De Vito et al. (2019, 2021) proposed a conceptually simple
multi-study FA model. This appealing model includes shared and study-specific
components via an additive expansion but is known to face identifiability
issues, discussed further below. Alternatively, Roy et al. (2021) proposed a
multiplicative ‘perturbed’ FA model that focuses on inferring the shared
structure while making use of subject-specific perturbations to resolve
identifiability issues in post-processing steps.

Motivated by surmounting the identifiability issues arising in multi-study
factor analysis particularly in the high-dimensional setting, we propose a
novel shared SUbspace Factor Analysis (SUFA) model. In particular, SUFA
assumes that there is a common lower-dimensional subspace shared across the
multiple studies under consideration. The data from different studies provide
partial information about the shared subspace, providing a mechanism for
borrowing of information. Exploiting the subspace structure we propose almost
sure solutions to identify the study-specific variations from the shared
covariance structure across all studies. Related work by Franks and Hoff
(2019) focuses on identifying the best shared subspace across groups. In
contrast, we focus on learning the shared versus study-specific contributors
to covariance structure. Our factor analysis approach allows one to infer
latent factors jointly with the structure of their loadings, providing
valuable interpretations in applications including gene expression studies
(Iacob et al., 2016).

We focus on identifiability up to the usual rotational ambiguity encountered
in factor modeling. Then, following standard practice in the Bayesian factor
modeling literature, we rotationally align the factors across MCMC iterations
in a post-processing stage. By using a fast algorithm recently introduced by
Poworoznek et al. (2021), we maintain sparsity in the post-processed samples
of the loadings by leveraging a key sparsity-inducing property of the varimax
rotation (Kaiser, 1958; Rohe and Zeng, 2020). Sparse loadings have
considerable interpretability advantages; for example, a subset of genes that
loads onto a particular factor is commonly interpreted as belonging to the
same pathway.

To provide theoretical support in terms of covariance estimation, we derive
favorable posterior contraction rates for subspace factor models in a high-
dimensional setting. Our analysis shows that the shared study-specific
covariance structures can be separately recovered with added precision for
each additional study, even when the marginal distributions of data differ
substantially across studies. Although the classical factor analysis
literature relies heavily on choosing the number of latent factors correctly,
we show that our covariance matrix estimation is robust to this choice.

Our model is developed alongside computational contributions to enable
efficient posterior inference. Latent variable-based Gibbs sampling
(Bhattacharya and Dunson, 2011; Sabnis et al., 2016) or expectation-
maximization (EM) algorithms (Ročková and George, 2016) are the most popular
approaches for Bayesian factor models. However, the conditional updating can
become a computational bottleneck for massive sample sizes, further
exacerbated in multi-study settings. Dai et al. (2020) proposed a fast matrix-
free approach for exploratory factor analysis but their approach does not
apply in Bayesian contexts. We develop a scalable (large $n$) Hamiltonian
Monte Carlo (HMC, Neal, 2011) sampler that jointly updates parameters while
marginalizing out latent factors. This leads to significant improvements in
mixing relative to Gibbs sampling. The proposed sampler depends only on the
study-specific sample covariance matrices, which can be computed and cached
prior to running MCMC. The computational complexity is then essentially
invariant of the sample size, and we design a distributed computing framework
to parallelize each MCMC iteration.

The remainder of the article is organized as follows: Section 2 describes the
SUFA model, and 2.1 discusses and resolves the potential identifiability
issues common to multi-study factor models. Section 2.2 discusses our choice
of prior and other important specifications. Section 3 establishes posterior
contraction rates in high-dimensional settings. These contributions are made
practical in Section 4, where we propose a scalable HMC algorithm for
posterior sampling. Section 5.1 compares SUFA with existing approaches via a
suite of simulation studies, and the methods are applied in an integrative
gene network analysis case study in Section 5.2. We close by discussing our
contributions and future directions in Section 6.

##  Bayesian SUbspace Factor Analysis (SUFA)

In this article, we consider the setting where data consist of $S$ studies
each comprising $d$-variate observations
${\mathbf{Y}}_{s,i}=(Y_{s,i,1},\dots,Y_{s,i,d})^{\rm T}$, $i=1,\dots,n_{s}$,
$s=1,\dots,S$ measured on the same set of features with
${\mathbf{Y}}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(\mbox{\boldmath$\mu$}_{s},\bm{\Sigma}_{s})$.
Without loss of generality, we assume $\mbox{\boldmath$\mu$}_{s}={\mathbf{0}}$
following standard practice to center the data. In our motivating application,
we seek to jointly learn from two microarray datasets and a bulkRNASeq
dataset: we have $S=3$ studies each analyzing $d=474$ genes. In studying the
correlation structure between genes, we expect a common fundamental
association between features across the studies, along with study-specific
variations. Hence, we let
$\bm{\Sigma}_{s}=\bm{\Sigma}+\mbox{\boldmath$\Gamma$}_{s}$ where $\bm{\Sigma}$
is a positive definite matrix quantifying the shared structure, and
$\mbox{\boldmath$\Gamma$}_{s}$ accounts for the respective study-specific
dependencies.

Bayesian factor analysis often yields state-of-the-art performance in single-
study high-dimensional covariance estimation (Bhattacharya and Dunson, 2011;
Pati et al., 2014; Ročková and George, 2016). Adopting a usual factor-analytic
covariance factorization as the sum of a low rank and diagonal matrix to the
shared $\bm{\Sigma}$, we let $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Delta}$, where $\bm{\Lambda}$ is a $d\times q$ factor loading matrix,
with $q\ll d$ the number of latent factors, and
$\bm{\Delta}=\mathrm{diag}(\delta_{1}^{2},\dots,\delta_{d}^{2})$ is a diagonal
matrix of idiosyncratic variances. Evidence in the literature suggests that
expression-level dependent error variances tend to produce better results
(Kepler et al., 2002), so we do not assume $\bm{\Delta}=\sigma^{2}{\bf
I}_{d}$.

To model the study-specific deviations, we let
$\mbox{\boldmath$\Gamma$}_{s}=\bm{\Lambda}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}\bm{\Lambda}^{\rm T}$, where ${\mathbf{A}}_{s}$ is a $q\times q_{s}$ matrix.
This leads to the following expression for the covariance specific to study
$s$: $\bm{\Sigma}_{s}=\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Lambda}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}\bm{\Lambda}^{\rm
T}+\bm{\Delta}.$ This can equivalently be written

${\mathbf{Y}}_{s,i}=\bm{\Lambda}\boldsymbol{\eta}_{s,i}+\bm{\Lambda}{\mathbf{A}}_{s}\mbox{\boldmath$\zeta$}_{s,i}+\mbox{\boldmath$\epsilon$}_{s,i},~{}\boldsymbol{\eta}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{q}),~{}\mbox{\boldmath$\zeta$}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q_{s}}({\mathbf{0}},{\mathbf{I}}_{q_{s}}),~{}\mbox{\boldmath$\epsilon$}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\bm{\Delta}),$
(1)

where $\boldsymbol{\eta}_{s,i}$ is a $q$ dimensional latent factor in the
shared subspace, $\mbox{\boldmath$\zeta$}_{s,i}$ are study-specific latent
factors of dimension $q_{s}$, and $\mbox{\boldmath$\epsilon$}_{s,i}$ are mean
zero Gaussian error terms. Because $\boldsymbol{\eta}_{s,i}$ are supported on
the same subspace, our hierarchical model allows borrowing of information
across studies in estimating $\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Delta}$. This is formalized in Section 3. Model (1) implies the
following marginal distributions

${\mathbf{Y}}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Lambda}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}\bm{\Lambda}^{\rm
T}+\bm{\Delta})\text{ for }s=1,\dots,S.$ (2)

Typically we take $q_{s}<q$, as integrative analyses of multiple studies are
meaningful when the signals are mostly shared across studies. Imposing that
$q_{s}<q$ allows $\bm{\Lambda}$ to be the dominant term explaining the
variation across the studies. The following sections make this intuition
rigorous, and show how it ensures identifiability.

##### Related works:

Under the same data organization, De Vito et al. (2019, 2021) define the
multi-study factor analysis (MSFA) model:

${\mathbf{Y}}_{s,i}=\bm{\Lambda}\boldsymbol{\eta}_{s,i}+\bm{\Phi}_{s}\mbox{\boldmath$\zeta$}_{s,i}+\mbox{\boldmath$\epsilon$}_{s,i},~{}~{}\boldsymbol{\eta}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{q}),~{}~{}\mbox{\boldmath$\zeta$}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q_{s}}({\mathbf{0}},{\mathbf{I}}_{q_{s}}),~{}~{}\mbox{\boldmath$\epsilon$}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\bm{\Delta}_{s}),$
(3)

where $\bm{\Lambda}$ is the $d\times q$ shared factor loading matrix and
$\bm{\Phi}_{s}$’s are the $d\times q_{s}$ study-specific loading matrices,
with the remaining notation matching ours above. The formulation is clearly
related to ours—in fact, it is strictly more general in that (1) can be viewed
as the case where the study-specific loadings take the form
$\bm{\Phi}_{s}=\bm{\Lambda}{\mathbf{A}}_{s}$. Indeed, our work is motivated by
observing that (3) is too flexible; in particular, allowing $\bm{\Phi}_{s}$ to
be arbitrary produces critical identifiability issues down the road. To see
this, consider its implied marginal distribution
${\mathbf{Y}}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Phi}_{s}\bm{\Phi}_{s}^{\rm T}+\bm{\Delta}_{s})$. Setting
$\bm{\Lambda}={\mathbf{0}}$ here, for instance, can always fit the data
equally well as the full model, as the effect of
$\bm{\Lambda}\bm{\Lambda}^{\rm T}$ can be absorbed entirely or partially into
the study-specific covariance $\bm{\Phi}_{s}\bm{\Phi}_{s}^{\rm
T}+\bm{\Delta}_{s}$ rendering a critical identifiability issue in interpreting
$\bm{\Lambda}\bm{\Lambda}^{\rm T}$ as the shared covariance part. In contrast,
(2) is more restrictive in not allowing $\bm{\Delta}$ to vary per study, and
with $\bm{\Lambda}$ parametrizing both the shared and study-specific terms.
Our contributions in Section 2.1 make precise how our proposed, more
parsimonious approach avoids these pitfalls. The combinatorial FA model
developed by Grabski et al. (2023) borrows information across related studies
by identifying latent factors that any combination of studies or groups can
share albeit without notions of shared and study-specific components.

Perturbed factor analysis (PFA, Roy et al., 2021) takes an altogether
different approach with the same aim of learning shared covariance structure:

${\mathbf{Q}}_{s}{\mathbf{Y}}_{s,i}=\bm{\Lambda}\boldsymbol{\eta}_{s,i}+\mbox{\boldmath$\epsilon$}_{s,i},\qquad\boldsymbol{\eta}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}({\mathbf{0}},{\mathbf{I}}_{q}),\qquad\mbox{\boldmath$\epsilon$}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}({\mathbf{0}},\bm{\Delta}),$
(4)

where $\bm{\Lambda}$ is the $d\times q$ common factor loading matrix,
${\mathbf{Q}}_{s}$ is a $d\times d$ perturbation matrix, and the remaining
notation as above. Though it overcomes some of the pitfalls of the MSFA
models, the study-specific effects cannot be recovered under PFA, and the
introduction of ${\mathbf{Q}}_{s}$ makes it difficult to scale beyond a few
hundred dimensions. Avalos-Pacheco et al. (2022) proposed the Bayesian factor
regression (BFR) model with additional hierarchies to handle different types
of batch effects but do not consider separate covariance structures across
studies.

###  Model Identifiablity Guarantees

Factor analytic models such as SUFA are prone to two key identifiability
issues:

  1. (I)

Information Switching: If there exists $q\times q$ non-null symmetric matrices
${\mathbf{C}}_{1}$ and ${\mathbf{C}}_{2}$ such that
${\mathbf{I}}_{q}-{\mathbf{C}}_{1}\succ{\mathbf{0}}$,
$\bm{\Delta}-\bm{\Lambda}{\mathbf{C}}_{2}\bm{\Lambda}^{\rm T}$ is a diagonal
matrix with positive entries and ${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}+{\mathbf{C}}_{1}+{\mathbf{C}}_{2}=\widetilde{{\mathbf{A}}}_{s}\widetilde{{\mathbf{A}}}_{s}^{\rm
T}$ for $q\times q_{s}$ matrices $\widetilde{{\mathbf{A}}}_{s}$’s for all
$s=1,\dots,S$, then substituting
$\widetilde{\bm{\Lambda}}=\bm{\Lambda}({\mathbf{I}}_{q}-{\mathbf{C}}_{1})^{\frac{1}{2}}$,
$\widetilde{\bm{\Delta}}=\bm{\Delta}-\bm{\Lambda}{\mathbf{C}}_{2}\bm{\Lambda}^{\rm
T}$ and $\widetilde{{\mathbf{A}}}_{s}$ yield the same marginal distribution in
(2), leading to an identifiability issue.

  2. (II)

Rotational Ambiguity:  Let $\widetilde{\bm{\Lambda}}=\bm{\Lambda}{\mathbf{H}}$
and $\widetilde{{\mathbf{A}}}_{s}={\mathbf{H}}^{\rm
T}{\mathbf{A}}_{s}{\mathbf{H}}_{s}$ where ${\mathbf{H}}$ and
${\mathbf{H}}_{s}$’s are orthogonal matrices of order $q\times q$ and
$q_{s}\times q_{s}$ respectively. Then substituting $\bm{\Lambda}$ and
${\mathbf{A}}_{s}$’s by $\widetilde{\bm{\Lambda}}$ and
$\widetilde{{\mathbf{A}}}_{s}$’s respectively in (2) yields identical marginal
distributions.

This information switching allows $\bm{\Lambda}{\mathbf{A}}_{s}$’s to have
shared column(s) across all studies posing a critical issue in interpreting
$\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\rm T}+\bm{\Delta}$ as the shared
covariance term, since $\bm{\Sigma}$ will not be identifiable. To establish
some intuition, we consider an example where $d=5$, $S=2$, $q=3$ with
$q_{1}=q_{2}=2$,

$\bm{\Lambda}=\begin{bmatrix}7&5&6\\\ 6&6&7\\\ 6&9&4\\\ 5&5&6\\\
4&6&6\end{bmatrix},\quad{\mathbf{A}}_{1}=\begin{bmatrix}3&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\\\
0&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}\\\
0&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{bmatrix}\quad\text{
and
}{\mathbf{A}}_{2}=\begin{bmatrix}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}&0\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&4\end{bmatrix}.$
(5)

Equation (1) implies that for each of the two “studies”, ${\mathbf{Y}}_{s,i}$
is given by

${\mathbf{Y}}_{1,i}=\left[\begin{array}[]{rrr|rrr}7&5&6&21&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}\\\
6&6&7&18&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\\\
6&9&4&18&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}18}\\\
5&5&6&15&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}\\\
\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}4&6&6\end{matrix}}}_{\text{$\bm{\Lambda}$}}}$}4&6&6&\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}12&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\end{matrix}}}_{\text{$\bm{\Lambda}{\mathbf{A}}_{1}$}}}$}12&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\\\
\end{array}\right]\begin{bmatrix}\boldsymbol{\eta}_{1,i}\\\
\mbox{\boldmath$\zeta$}_{1,i}\end{bmatrix}+\mbox{\boldmath$\epsilon$}_{1,i}=\left[\begin{array}[]{rrrr|rr}7&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&21\\\
6&6&7&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&18\\\
6&9&4&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}18}&18\\\
5&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&15\\\
\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}4&6&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\end{matrix}}}_{\text{$\widetilde{\bm{\Lambda}}$}}}$}4&6&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}12\end{matrix}}}_{\text{$\widetilde{\bm{\Lambda}}\widetilde{{\mathbf{A}}}_{1}$}}}$}12\\\
\end{array}\right]\begin{bmatrix}\widetilde{\boldsymbol{\eta}}_{1,i}\\\
\widetilde{\mbox{\boldmath$\zeta$}}_{1,i}\end{bmatrix}+\mbox{\boldmath$\epsilon$}_{1,i},$
(6)
${\mathbf{Y}}_{2,i}=\left[\begin{array}[]{rrr|rrr}7&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&24\\\
6&6&7&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&28\\\
6&9&4&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}18}&16\\\
5&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&24\\\
\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}4&6&6\end{matrix}}}_{\text{$\bm{\Lambda}$}}}$}4&6&6&\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&24\end{matrix}}}_{\text{$\bm{\Lambda}{\mathbf{A}}_{2}$}}}$}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&24\\\
\end{array}\right]\begin{bmatrix}\boldsymbol{\eta}_{2,i}\\\
\mbox{\boldmath$\zeta$}_{2,i}\end{bmatrix}+\mbox{\boldmath$\epsilon$}_{2,i}=\left[\begin{array}[]{rrrr|rr}7&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&24\\\
6&6&7&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&28\\\
6&9&4&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}18}&16\\\
5&5&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&24\\\
\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}4&6&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\end{matrix}}}_{\text{$\widetilde{\bm{\Lambda}}$}}}$}4&6&6&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&\makebox[0.0pt][l]{$\smash{\underbrace{\phantom{\begin{matrix}24\end{matrix}}}_{\text{$\widetilde{\bm{\Lambda}}\widetilde{{\mathbf{A}}}_{2}$}}}$}24\\\
\end{array}\right]\begin{bmatrix}\widetilde{\boldsymbol{\eta}}_{2,i}\\\
\widetilde{\mbox{\boldmath$\zeta$}}_{2,i}\end{bmatrix}+\mbox{\boldmath$\epsilon$}_{2,i}.$
(7)

That the above equations can each be written in two equivalent decompositions
illustrates the information switching problem: we cannot distinguish whether
the column
$\bm{\Lambda}\times{\scriptsize\begin{bmatrix}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\end{bmatrix}^{\rm
T}}={\scriptsize\begin{bmatrix}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}18}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}12}\end{bmatrix}^{\rm
T}}$ is a part of the shared effect or the study-specific effect. Although the
${\mathbf{A}}_{s}$ matrices should take into account the study-specific
variations only, the individual effect of one study in this example can be at
least “partially explained” by the others. From the perspective of the
statistical model, ${{\mathbf{A}}_{s}}$’s are no longer study-specific, but
absorb a part of the shared variation which should be captured entirely by
$\bm{\Sigma}$. In Section 2.1.1 we identify the necessary and sufficient
condition causing this issue and propose an almost sure solution.

Rotational ambiguity has been well documented in the prior literature,
affecting both MSFA and PFA. In addition, MSFA can suffer from information
switching—briefly, for a $d\times d$ order symmetric matrix
$\widetilde{{\mathbf{C}}}$ such that $\bm{\Lambda}\bm{\Lambda}^{\rm
T}-\widetilde{{\mathbf{C}}}\succ{\mathbf{0}}$,
$\bm{\Phi}_{s}\bm{\Phi}_{s}^{\rm
T}+\widetilde{{\mathbf{C}}}=\widetilde{\bm{\Phi}}_{s}\widetilde{\bm{\Phi}}_{s}^{\rm
T}$ for some $d\times q_{s}$ order matrix $\widetilde{\bm{\Phi}}_{s}$. To
resolve these, De Vito et al. (2019) restrict $\bm{\Lambda}$ and
$\bm{\Phi}_{s}$’s to be lower-triangular matrices while enforcing the
augmented matrix
$\begin{bmatrix}\bm{\Lambda}&\bm{\Phi}_{1}&\cdots&\bm{\Phi}_{S}\end{bmatrix}$
to be of full column-rank while obtaining the MLE. Lower-triangular structures
incur the tradeoff of introducing an order dependence in the variables
(Frühwirth-Schnatter and Lopes, 2018); when no natural ordering exists,
results can be very sensitive to permutations of the variables (Carvalho et
al., 2008). These structural constraints also can affect mixing and yield
inconsistent estimates (Millsap, 2001; Erosheva and Curtis, 2017). In high-
dimensional applications often sparsity is assumed on the loading matrices
(Daniele et al., 2019); with exact zero entries in $\bm{\Lambda}$ and
$\bm{\Phi}_{s}$’s, ensuring full column-rank of
$\begin{bmatrix}\bm{\Lambda}&\bm{\Phi}_{1}&\cdots&\bm{\Phi}_{S}\end{bmatrix}$
is non-trivial. The Bayesian MSFA by De Vito et al. (2021) assumes
multiplicative gamma process priors on the loading matrices without any
restrictions on the column dimensions and thus can be susceptible to both (I)
and (II). On the other hand, Roy et al. (2021) set
${\mathbf{Q}}_{1}={\mathbf{I}}_{d}$ to impose identifiability in PFA, which
can again be sensitive to the choice of the “first” study. In the following
section, we provide almost sure solutions to the identifiability issues under
much milder conditions, avoiding any such structural assumptions.

####  Resolving Information Switching

We now derive the necessary and sufficient conditions for information
switching.

######  Lemma 1.

Assume that the data admits the marginal distribution in (2) for all the
studies. Then information switching occurs _if and only if_ there exists
$\bm{\Lambda}$, $\bm{\Delta}$ and ${\mathbf{A}}_{s}$’s such that
${\mathbf{Y}}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\bm{\Lambda}\bm{\Lambda}^{\rm
T}+\bm{\Lambda}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}\bm{\Lambda}^{\rm
T}+\bm{\Delta})$ and
$\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ is non-null or
$\mathrm{rank}({\mathbf{A}}_{s})<q_{s}$ for all $s$, where
$\mathbb{C}\left({\mathbf{A}}\right)$ denotes the column space of a matrix
${\mathbf{A}}$.

Intuitively, a non-null
$\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ implies that the
${\mathbf{A}}_{s}$ matrices are no longer study-specific. As we saw in example
(5), ${\mathbf{A}}_{s}$ can absorb part of the shared variation which is meant
to be captured by $\bm{\Sigma}$. Rank deficiencies in the columns leave space
for the ${\mathbf{A}}_{s}$’s to partially absorb the shared information. The
Lemma provides requisite insights to avoid this with a very simple condition
to ensure $\bm{\Sigma}$ fully explains the shared covariance.

######  Theorem 1.

If $\sum_{s=1}^{S}q_{s}\leq q$ then information switching has zero support
under any non-degenerate continuous prior on ${\mathbf{A}}_{s}$.

Intuitively, aggregating multiple studies is a fruitful pursuit only when the
datasets share enough similarity or structure; it is natural to expect that
our subspace factor model is well-posed when there are fewer study-specific
latent factors than shared factors. Theorem 1 makes this intuition precise,
showing that the potential problem of information switching is resolved almost
surely, i.e., generic identifiability of $\bm{\Sigma}$ (Allman et al., 2009)
is ensured under any continuous prior on ${\mathbf{A}}_{s}$ as long as
$\sum_{s=1}^{S}q_{s}\leq q$. The assumption is therefore a natural one, and
$\bm{\Sigma}$ is completely identifiable as a result, allowing for us to study
and interpret the interaction between variables. To illustrate concretely in
the context of the above example, consider setting the shared
$\widetilde{\bm{\Lambda}}$ to be a $5\times 4$ matrix and $q_{1}=q_{2}=1$,
i.e., right-hand side of (6)-(7). Now, $\widetilde{\bm{\Lambda}}$ accounts for
the variability in the marginal distribution (2), and in turn the parts
explained by ${\mathbf{A}}_{1}$ and ${\mathbf{A}}_{2}$ no longer have
intersection. Notably, the strategy of restricting the augmented matrix
$\begin{bmatrix}\bm{\Lambda}&\bm{\Phi}_{1}&\cdots&\bm{\Phi}_{S}\end{bmatrix}$
to be of full column-rank by De Vito et al. (2019) also requires a similar
condition $dq-q(q-1)/2+\sum_{s=1}^{S}\\{dq_{s}-q_{s}(q_{s}-1)/2\\}\leq
Sd(d-1)/2$.

####  Resolving Rotational Ambiguity

Identifiability with respect to orthogonal rotation (II) is not essential for
estimating the shared covariance $\bm{\Sigma}$. However, rotational ambiguity
does create obstacles to inferring interpretable lower-dimensional factors
(Russell, 2002; Iacob et al., 2016). Notably $\bm{\Lambda}\bm{\Lambda}^{\rm
T}$ is identifiable from the diagonal $\bm{\Delta}$ if $q<(2d-\sqrt{8d+1})/2$
(Bekker and ten Berge, 1997) and with this condition it is thus common to
address (II) via post-processing the $\bm{\Lambda}$ samples (Roy et al., 2021;
De Vito et al., 2021; Papastamoulis and Ntzoufras, 2022, and others). Such
approaches avoid order dependence arising from structural assumptions.

Interpretability of the factor loadings is greatly enhanced by sparsity—for
example, the sparse subset of genes with nonzero loadings onto a common factor
can be associated with a common biological pathway. Even when one places a
shrinkage prior on the loadings to favor (near) sparsity, this desired
structure is potentially destoyed after applying post-processing algorithms to
the MCMC samples of the loadings. To avoid this, we adopt the approach of
Poworoznek et al. (2021). The method first tackles the generic rotational
invariance across the MCMC samples using an orthogonalization step leveraging
the varimax transformation (Kaiser, 1958) to resolve ambiguity up to switching
of the column labels and signs. These ambiguities are both resolved in the
next step by matching each MCMC sample to a reference matrix called a pivot,
aligning samples via a greedy maximization scheme.

The post-processed MCMC samples can thus be considered matched and aligned
with respect to a common orthogonal transformation for inference downstream.
Critically, the varimax rotation in the orthogonalization step implicitly
induces sparsity (Rohe and Zeng, 2020). In more detail, the varimax criterion
is given by

$\textstyle{{\mathbf{H}}_{\mathrm{VARIMAX}}={\arg\max}_{{\mathbf{H}}}\left[\frac{1}{d}\sum_{h=1}^{q}\sum_{j=1}^{d}(\bm{\Lambda}{\mathbf{H}})_{j,h}^{4}-\sum_{h=1}^{q}\left\\{\frac{1}{d}\sum_{j=1}^{d}(\bm{\Lambda}{\mathbf{H}})_{j,h}^{2}\right\\}^{2}\right]},$

and describes the optimal rotation maximizing the sum of the variances of
squared loadings. Intuitively, this is achieved if (i) any given variable has
a high loading on a single factor but near-zero loadings on the others, and
(ii) any given factor is constituted by a few variables with very high
loadings but near-zero support from the others. As a result, the summary
matrix obtained from the posterior MCMC samples of $\bm{\Lambda}$ is sparse
upon applying varimax rotation, while the marginal distribution of the data is
unaffected. Section S.4.2 of the supplementary materials contains thorough
simulations showing see that this strategy accurately recovers $\bm{\Lambda}$
consistently across several realistic scenarios. As
$\bm{\Lambda}^{(s)}=\bm{\Lambda}{\mathbf{A}}_{s}$ can be defined to denote the
study-specific loading matrices, the above prescription applies similarly to
obtaining study-specific loadings from samples of $\bm{\Lambda}^{(s)}$.

###  Latent Dimensions and Prior Specification

##### Latent Dimensions:

In most practical applications, we do not know the column dimensions of
$\bm{\Lambda}$ and ${\mathbf{A}}_{s}$, denoted by $q$ and $q_{s}$,
respectively. Although one could choose priors on $q$ and $q_{s}$ and
implement reversible-jump algorithms (Green, 1995), this will often lead to
inefficient computation, particularly in large $d$ cases. Instead, it has
become common practice in the single-study literature to use overfitted factor
models by first fixing $q$ at some upper bound and then leveraging appropriate
priors to shrink the extra columns (Legramanti et al., 2020; Schiavon et al.,
2022, and related others). This strategy substantially simplifies MCMC
implementations. To choose upper bounds on the number of factors in practice,
we employ augmented implicitly restarted Lanczos bidiagonalization (Baglama
and Reichel, 2005) to obtain approximate singular values and eigenvectors of
the pooled dataset, choosing the smallest $\widehat{q}$ that explains at least
$95\%$ of the variability in the data. To ensure identifiablity of
$\bm{\Lambda}\bm{\Lambda}^{\rm T}$ from $\bm{\Delta}$, we restrict
$\widehat{q}<(2d-\sqrt{8d+1})/2$ (Bekker and ten Berge, 1997). After doing so,
the simple choice $\widehat{q_{s}}=\widehat{q}/S$ satisfies the conditions in
Theorem 1. We later show in Section 3 that recovering the marginal covariance
is asymptotically robust with respect to the choice of $q$ and $q_{s}$ under
appropriate priors on $\bm{\Lambda}$.

##### Prior Specifications:

In high-dimensional applications, it is important to reduce the number of
parameters in the factor loadings matrices for statistical efficiency in
estimating the shared and study-specific contributions to the large covariance
matrix. Among a wide variety of appropriate shrinkage priors, we use the
Dirichlet-Laplace (DL, Bhattacharya et al., 2015) for its near-minmax optimal
contraction rates in the single-study context (Pati et al., 2014) and
computational simplicity. We let
$\mathrm{vec}(\bm{\Lambda})\sim\mathrm{DL}(a)$ where $a$ is a suitably chosen
hyperparameter; details appear in Section S.1 of the supplementary materials.
For the study-specific terms ${\mathbf{A}}_{s}$, we let its entries
$a_{s,i,j}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(0,b_{{\mathbf{A}}})$,
a choice that avoids information switching almost surely. We avoid shrinkage
priors on ${\mathbf{A}}_{s}$’s as in practice it can potentially lead to
concentration close to parameter values that violate the identifiability
conditions of Lemma 1. Regarding the idiosyncratic variances, we consider
$\log\delta_{j}^{2}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(\mu_{\delta},\sigma^{2}_{\delta})$
as log-normal distributions have modes bounded away from zero and tend to
produce more numerically stable estimates compared to commonly used inverse-
gamma and half-Cauchy priors (Gelman, 2006).

##  Posterior Contraction Rates

We now analyze the posterior contraction rates for recovering the shared and
study-specific covariance matrices. Before detailing the assumptions on the
true data-generating mechanism and the postulated SUFA model, we fix
notational conventions. We let $\Pi_{n}(\cdot)$ denote the prior and
$\Pi_{n}(\cdot\mid\mathcal{D}_{n})$ the corresponding posterior given data
$\mathcal{D}_{n}=\\{{\mathbf{Y}}_{1,1},\dots,{\mathbf{Y}}_{1,n_{1}},\dots,{\mathbf{Y}}_{S,1},\dots,{\mathbf{Y}}_{S,n_{S}}\\}$.
Throughout, $\norm{{\mathbf{A}}}_{F}$ and $\norm{{\mathbf{A}}}_{2}$ denote the
Frobenius and spectral norms of a matrix ${\mathbf{A}}$, respectively. For
real sequences $\\{a_{n}\\}$, $\\{b_{n}\\}$, $a_{n}=o(b_{n})$ implies that
$\lim\absolutevalue{{a_{n}}/{b_{n}}}=0$ and $a_{n}\succ b_{n}$ implies that
$\liminf{a_{n}}/{b_{n}}>0$. A $d$-dimensional vector $\theta$ is said to be
$s$-sparse if it has only $s$ nonzero elements, and we denote the set of all
$s$-sparse vectors in $\mathbb{R}^{d}$ by $\ell_{0}\left[s,d\right]$.

##### Data-generating mechanism:

We assume that the data are generated according to
${\mathbf{Y}}_{s,1:n_{s}}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Sigma}_{0sn})$
for each study $s=1,\dots,S$. Here,
$\bm{\Sigma}_{0sn}=\bm{\Lambda}_{0n}\bm{\Lambda}_{0n}^{\rm
T}+\bm{\Lambda}_{0n}{\mathbf{A}}_{0sn}{\mathbf{A}}_{0sn}^{\rm
T}\bm{\Lambda}_{0n}^{\rm T}+\mbox{\boldmath$\Delta$}_{0n}$ with
$\bm{\Lambda}_{0n}$ a $d_{n}\times q_{0n}$ sparse matrix, ${\mathbf{A}}_{0sn}$
a $q_{0n}\times q_{0sn}$ matrix with real entries, and
$\mbox{\boldmath$\Delta$}_{0n}:=\mathrm{diag}(\delta^{2}_{01},\dots,\delta^{2}_{0d_{n}})$.
Compared to the setup in prior work by Pati et al. (2014), we consider
multiple studies and allow heterogeneous idiosyncratic errors, which will
provide more realism in modeling the data as well as require a more nuanced
analysis. Denoting $n=\sum_{s=1}^{S}n_{s}$ the combined sample size across
studies, we will allow the model parameters to increase in dimension with $n$
as suggested by the subscripts in our notations. Under this setup, we now
state the required sparsity conditions on the true parameters in order to
recover the shared as well as the study specific covariance structures in
high-dimensional settings.

  1. (C1)

For $s=1,\dots,S$, $\liminf_{n\to\infty}\frac{n_{s}}{n}>0$.

  2. (C2)

Let $\\{q_{0n}\\}$, $\\{q_{0sn}\\}$ for $s=1,\dots,S$, $\\{s_{n}\\}$ and
$\\{d_{n}\\}$ be increasing sequences of positive integers such that
$\sum_{s=1}^{S}q_{0sn}\leq q_{0n}<d_{n}$.

  3. (C3)

$\bm{\Lambda}_{0n}$ is a $d_{n}\times q_{0n}$ full rank matrix such that each
column of $\bm{\Lambda}_{0n}$ belongs to $\ell_{0}\left[s_{n},d_{n}\right]$,
$\norm{\bm{\Lambda}_{0n}}_{2}^{2}=o(\frac{c_{n}}{q_{0n}})$ and $\max_{1\leq
s\leq S}\norm{{\mathbf{A}}_{0sn}}_{2}^{2}=o(q_{0n})$ where $\\{c_{n}\\}$ is an
increasing sequence of positive real numbers.

  4. (C4)

$\max_{j}\delta^{2}_{0d_{n}}=o(c_{n})$ and
$\min_{j}\delta^{2}_{0d_{n}}>{\delta}_{\min}^{2}$ where ${\delta}_{\min}^{2}$
is a positive constant.

  5. (C5)

$\frac{d_{n}\max\\{(\log c_{n})^{2},\log
d_{n}\\}}{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{0n})}=o(1)$.

In less technical terms, (C1) ensures that none of the studies has a
negligible proportion of the data. (C2) and (C3) specify the requisite
sparsity conditions. Together, (C3) and (C4) imply that the marginal variances
grow in $o(c_{n})$ while ensuring that they are also bounded away from 0. The
technical condition (C5) specifies the relative rate of increment of the
parameters.

##### Specifics of the posited SUFA model:

In practice, the latent dimensions are unknown, and it is desirable to
establish theory for the strategy of overspecifying a model and then shrinking
redundant parameters via shrinkage priors. This section identifies conditions
under which we may derive guarantees in this overparametrized setting. To this
end, let
${\mathbf{Y}}_{s,1:n_{s}}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Lambda}_{n}\bm{\Lambda}_{n}^{\rm
T}+\bm{\Lambda}_{n}{\mathbf{A}}_{sn}{\mathbf{A}}_{sn}^{\rm
T}\bm{\Lambda}_{n}^{\rm T}+\mbox{\boldmath$\Delta$}_{n})$ where
$\bm{\Lambda}_{n}$ and ${\mathbf{A}}_{sn}$ are $d_{n}\times q_{n}$ and
$q_{n}\times q_{sn}$ matrices, respectively, and
$\mbox{\boldmath$\Delta$}_{n}=\mathrm{diag}(\delta^{2}_{1},\dots,\delta^{2}_{d_{n}})$.
The parameters in the posited (possibly misspecified) model do not have $0$
subscripts, to distinguish from their ground truth counterparts.

We let $\\{q_{n}\\}$, $\\{q_{sn}\\}$ be increasing sequences of positive
integers such that $q_{0n}\leq q_{n}$, $q_{0sn}\leq q_{sn}$ so that the chosen
number of shared and study-specific latent factors upper bound the respective
true (also unknown) numbers resulting in a misspecified and over-parametrized
factor model in the postulation. We additionally assume
$\sum_{s=1}^{S}q_{sn}\leq q_{n}$ to ensure the identifiability condition from
Theorem 1. As discussed in Section 2.2 we let
$\mathrm{vec}(\bm{\Lambda}_{n})\sim\mathrm{DL}(a_{n})$ with
$a_{n}={1}/{d_{n}q_{n}}$ to impose sparsity on $\bm{\Lambda}_{n}$. Below, we
provide sufficient conditions on the overspecified dimensions to ensure
consistency of the class of postulated models even when they are misspecified.

  1. (D1)

$q_{n}^{2}=o(c_{n}^{2}s_{n}q_{0n})$ and
$\lim_{n\to\infty}{\frac{c_{n}^{12}}{n}}{\left\\{s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{3}}=0$.

  2. (D2)

$\log(d_{n}q_{n})\min\left\\{\frac{c_{n}}{s_{n}q_{0n}^{3}}(d_{n}q_{n})^{c_{n}^{2}},s_{n}(d_{n}q_{n})^{\frac{c_{n}^{2}}{q_{n}^{2}}s_{n}q_{0n}},\frac{s_{n}q_{0n}}{(\log
c_{n})^{2}}(d_{n}q_{n})^{\frac{c_{n}^{2}}{d_{n}}s_{n}q_{0n}}\right\\}\succ n$.

The first condition (D1) makes explicit how much overparametrization can be
tolerated, while the second condition specifies the signal-to-noise ratio with
respect to the true column dimension $q_{0n}$, upper bound $c_{n}$ of the
marginal variances and number of non-zero elements $s_{n}$ in each column of
$\bm{\Lambda}$. The second condition in (D1) and (C5) jointly imply that
$d_{n}\log d_{n}$ grows in $o(n^{\alpha})$ for some $0<\alpha<1$. (D2) can be
interpreted as a lower bound on the signal strength in order to recover the
covariance structures from the data.

##### Contraction results:

Let
$\mbox{\boldmath$\Theta$}_{0n}=\\{\bm{\Lambda}_{0n},\mbox{\boldmath$\Delta$}_{0n},{\mathbf{A}}_{01n},\dots,{\mathbf{A}}_{0Sn}\\}$
be the true data-generating parameters and $\mathbb{P}_{0n}$ denote the joint
distribution of $\mathcal{D}_{n}$ under $\mbox{\boldmath$\Theta$}_{0n}$. We
define $\bm{\Sigma}_{0n}=\bm{\Lambda}_{0n}\bm{\Lambda}_{0n}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n}$ as the true shared covariance structure and
$\bm{\Sigma}_{n}=\bm{\Lambda}_{n}\bm{\Lambda}_{n}^{\rm
T}+\mbox{\boldmath$\Delta$}_{n}$ to be the shared covariance in the postulated
model. The following result formally characterizes the posterior contraction
properties around $\bm{\Sigma}_{0n}$.

######  Theorem 2 (Contraction rate).

For any $\mbox{\boldmath$\Theta$}_{0n}$ satisfying (C1)-(C5) and priors
satisfying (D1)-(D2), we have
$\mathbb{E}_{\mathbb{P}_{0n}}\lim_{n\to\infty}\Pi_{n}\left(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n}\right)=0$,
where
$\varepsilon_{n}=c_{n}^{6}\sqrt{\frac{\\{s_{n}q_{0n}\log(d_{n}q_{n})\\}^{3}}{n_{1}+\dots+n_{S}}}$
and $M>0$ is a large enough constant.

Theorem 2 has several important implications for modeling and computation.
First, the combined sample size appears in the denominator of the contraction
rate $\varepsilon_{n}$, which is consistent with the intuition that learning
the shared covariance structure improves inference by borrowing strength
across multiple studies or data views. Because the true values of $q_{0n}$ and
$q_{0sn}$ are unknown, we made the weak assumption that the specified column
dimensions are larger than the ground truth dimensions. Although the best
contraction rate is unsurprisingly achieved when correctly setting
$q_{n}=q_{0n}$, Theorem 2 validates the practice of beginning with
overparametrized models, showing that it is still possible to recover
$\bm{\Sigma}_{0n}$. This result provides an important theoretical foundation
to the common heuristic strategy of setting the column dimension of
$\bm{\Lambda}$ to some upper bound in FA model implementations.

Note that $\bm{\Lambda}_{0n}{\mathbf{A}}_{0sn}{\mathbf{A}}_{0sn}^{\rm
T}\bm{\Lambda}_{0n}^{\rm T}$s are the study-specific covariance structures.
The following result shows that these study-specific terms can also be
recovered for all $s$.

######  Corollary 1 (Recovering individual structures).

Under the same conditions as Theorem 2,
$\lim_{n\to\infty}\mathbb{E}_{\mathbb{P}_{0n}}\Pi_{n}\left(\norm{\bm{\Lambda}_{n}{\mathbf{A}}_{sn}{\mathbf{A}}_{sn}^{\rm
T}\bm{\Lambda}_{n}^{\rm
T}-\bm{\Lambda}_{0n}{\mathbf{A}}_{0sn}{\mathbf{A}}_{0sn}^{\rm
T}\bm{\Lambda}_{0n}^{\rm
T}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n}\right)=0$.

##  Gradient-based Posterior Computation

Gibbs sampling is the most popular approach for posterior inference in
Bayesian factor models. While such an approach is straightforward for our
model (1), it is well known that alternately updating latent factors and
covariance structure parameters can lead to slow mixing. In addition,
instantiating the latent factors $\boldsymbol{\eta}_{s,i}$ and
$\mbox{\boldmath$\zeta$}_{s,i}$ quickly increases computational cost for large
samples. It has been observed that marginalizing out auxiliary parameters can
often dramatically improve MCMC performance (Robert and Roberts, 2021). In
view of this, we develop a Hamiltonian Monte Carlo (HMC, Neal, 2011) within
Gibbs sampler that makes use of the marginal posterior, after integrating out
all latent factors. The HMC sampler makes smarter proposals utilizing the
gradient of the $\log$-target density. We will see that our proposed algorithm
confers substantial gains in real computation times compared to standard Gibbs
approaches.

Denote the log-likelihood of the marginal SUFA model (2) by ${\cal
L}:=K-\frac{1}{2}\sum_{s=1}^{S}\\{n_{s}\log\absolutevalue{\bm{\Sigma}_{s}}+\tr(\bm{\Sigma}_{s}^{-1}{\mathbf{W}}_{s})\\}$
with marginal covariance matrix
$\bm{\Sigma}_{s}=\bm{\Lambda}({\mathbf{I}}_{q}+{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T})\bm{\Lambda}^{\rm T}+\bm{\Delta}$, sample sum of squares matrix
${\mathbf{W}}_{s}=\sum_{i=1}^{n_{s}}{\mathbf{Y}}_{s,i}{\mathbf{Y}}_{s,i}^{\rm
T}$ for study $s$ and a constant $K$. The DL prior on $\bm{\Lambda}$ admits a
Gaussian distribution conditionally on the hyperparameters $\tau$, $\psi$ and
$\phi$ (see (S.2) in the supplementary materials for details). Thus from the
prior specifications outlined in Section 2.2, the posterior density of
$\mbox{\boldmath$\Theta$}:=(\bm{\Lambda},\bm{\Delta},{\mathbf{A}}_{1},\dots,{\mathbf{A}}_{S})$
given $\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$},\tau$ and the observed
data is given by

$\textstyle{\Pi(\mbox{\boldmath$\Theta$}\mid-)=\exp({\cal
L})\times\Pi(\bm{\Lambda}\mid\tau,\mbox{\boldmath$\phi$},\mbox{\boldmath$\psi$})\times\Pi(\bm{\Delta})\times\prod_{s=1}^{S}\Pi({\mathbf{A}}_{s}).}$
(8)

To implement HMC, gradients of $\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$ with
respect to $\Theta$ are required (see Section S.3 in the supplementary
materials for a general discussion on HMC algorithms). In Section A.2 of the
Appendix we show that the gradients can be analytically expressed and
executable via distributed parallel operations. Thus we avoid expensive
automatic differentiation algorithms (Carpenter et al., 2015) used by
probabilistic samplers like Stan. Equipped with these gradients, we propose an
HMC-within-Gibbs sampler where $\Theta$ is updated given
$(\tau,\mbox{\boldmath$\phi$},\mbox{\boldmath$\psi$})$ using an HMC step, and
then the hyperparameters
$(\tau,\mbox{\boldmath$\phi$},\mbox{\boldmath$\psi$})$ are updated
conditionally on $\bm{\Lambda}$ using a Gibbs step following Bhattacharya et
al. (2015).

1 Initialize $\Theta$, $\psi$, $\phi$ and $\tau$.

2 for  _$i=1,\dots,N$_ do

3HMC Update $\Theta$:

Draw ${\mathbf{p}}\sim\mathrm{N}({\mathbf{0}},{\mathbf{M}})$ and initiate
$\\{\mbox{\boldmath$\Theta$}(0),{\mathbf{p}}(0)\\}=(\mbox{\boldmath$\Theta$}_{i-1},{\mathbf{p}})$.
Repeat the following $L$ leapfrog steps to numerically solve Hamiltonian
equations:

4 for  _$\ell=0,\dots,L-1$_ do

5 $\displaystyle\textstyle{\mbox{\boldmath$\Theta$}\\{t+(\ell+1)\delta
t\\}=\mbox{\boldmath$\Theta$}(t+\ell\delta t)+\delta
t{\mathbf{M}}^{-1}\left[{\mathbf{p}}(t+\ell\delta t)-\frac{\delta t}{2}\nabla
V\left\\{\mbox{\boldmath$\Theta$}(t+\ell\delta t)\right\\}\right]},$
$\displaystyle\textstyle{{\mathbf{p}}\\{t+(\ell+1)\delta
t\\}={\mathbf{p}}(t+\ell\delta t)-\frac{\delta t}{2}\left[\nabla
V\left\\{\mbox{\boldmath$\Theta$}(t+\ell\delta t)\right\\}+\nabla
V\left\\{\mbox{\boldmath$\Theta$}(t+(\ell+1)\delta t)\right\\}\right]}.$

6 end for

Letting
$(\mbox{\boldmath$\Theta$}^{\prime},{\mathbf{p}}^{\prime})=(\mbox{\boldmath$\Theta$}(t+L\delta
t),{\mathbf{p}}(t+L\delta t))$, set
$\mbox{\boldmath$\Theta$}_{i}=\mbox{\boldmath$\Theta$}^{\prime}$ with
probability
$\min\left\\{1,e^{H\left(\mbox{\boldmath$\Theta$}^{\prime},{\mathbf{p}}^{\prime}\right)-H\left(\mbox{\boldmath$\Theta$}_{i-1},{\mathbf{p}}\right)}\right\\},$
where
$H(\mbox{\boldmath$\Theta$},{\mathbf{p}})=-\log\Pi(\mbox{\boldmath$\Theta$}\mid-)+{\mathbf{p}}^{\rm
T}{\mathbf{M}}^{-1}{\mathbf{p}}/2$.

* 7Gibbs Update $\psi$, $\phi$ and $\tau$:

  1. (i)

For $j=1,\dots,d$ and $h=1,\dots q$ sample $\widetilde{\psi}_{j,h}$
independently from inverse-Gaussian distributions
$\mathrm{iG}\left(\tau_{i-1}{\phi_{i-1,j,h}}/{\absolutevalue{\lambda_{i,j,h}}},1\right)$
and set $\psi_{i,j,h}=1/\widetilde{\psi}_{j,h}$.

  2. (ii)

Draw $T_{j,h}$ independently from from a generalized inverse Gaussian
distributions $\mathrm{giG}(a-1,1,2\absolutevalue{\lambda_{i,j,h}})$ and set
$\phi_{i,j,h}=T_{j,h}/T$ with $T=\sum_{j,h}T_{j,h}$.

  3. (iii)

Sample
$\tau_{i}\sim\mbox{giG}\left\\{dq(1-a),1,2\sum_{j,h}\frac{\absolutevalue{\lambda_{i,j,h}}}{\phi_{i,j,h}}\right\\}$.

8 end for

Algorithm 1 HMC-within-Gibbs Sampler for SUFA

The vanilla HMC algorithm requires the following tuning parameters: a positive
real number $\delta t$, a positive integer $L$ and a positive definite matrix
${\mathbf{M}}$ of the same order as $\Theta$. Letting $\nabla
V(\mbox{\boldmath$\Theta$})=-\partialderivative{\mbox{\boldmath$\Theta$}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$
be the gradient of the negative $\log$-posterior, we outline the HMC-within-
Gibbs sampler in Algorithm 1 to obtain $N$ MCMC samples from the joint
posterior distribution of
$(\mbox{\boldmath$\Theta$},\tau,\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$})$.
Implementing Algorithm 1 requires computing
$\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$ as well as
$\partialderivative{\mbox{\boldmath$\Theta$}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$
at each MCMC step. In Section A.2 of the Appendix we show that these
quantities can be computed very efficiently by distributing the calculations
over parallel processes.

Before validating theory and methods in empirical studies, we highlight two
advantages of our sampler. First, it can be adapted to any prior with
(conditionally) differentiable density, including the spike-and-slab (Ishwaran
and Rao, 2005), spike-and-slab lasso (Ročková and George, 2016), horseshoe
(Carvalho et al., 2009), multiplicative gamma (Bhattacharya and Dunson, 2011),
generalized double Pareto (Armagan et al., 2013), cumulative shrinkage
(Legramanti et al., 2020), and others. Almost all of these admit a
conditionally Gaussian prior on $\bm{\Lambda}$, so that adapting Algorithm 1
requires only modifying the Gibbs step. Second, a crucial advantage of
Algorithm 1 is its scalability with respect to the sample size. In particular,
the log-likelihood ${\cal L}$ and its gradient depend on the data only through
the study-specific sample sum of squares ${\mathbf{W}}_{s}$, which just needs
to be computed and cached once. Choices of the tuning parameters are discussed
in Section S.3.1 of the supplementary materials.

##  Empirical Study

###  Simulated Experiments

To assess the proposed methodology, we generate synthetic data according to
the true model (3) under three scenarios induced by different $\bm{\Lambda}$
matrices, which we refer to as FM1, FM2 and FM3. The ground truth shared
covariance structures are displayed in the first column of Figure 2; complete
details of the simulation setup and further figures appear in Section S.4.1 of
the supplementary materials. For each scenario, we also examine performance
under two misspecified cases: (i) slight misspecification, where
$\bm{\Phi}_{s}=\bm{\Lambda}{\mathbf{A}}_{s}+{\mathbf{E}}_{s}$ and
${\mathbf{E}}_{s}$ is a matrix of randomly generated errors, and (ii) complete
misspecification, where $\bm{\Phi}_{s}$’s are in the null space of
$\mathbb{C}\left(\bm{\Lambda}\right)$ for all $s=1,\dots,S=5$. We vary
$d=50,200$ and $450$ and set
$n_{1},\dots,n_{S}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\max\\{\mathrm{Poisson}(\frac{d}{S}),\frac{d}{S}\\}$
for each $d$. We generate $q_{s}\sim\mathrm{Poisson}(\frac{q}{S})$ where
$q=10$ for the $d=50$ and $q=20$ otherwise.

##### Peer methods:

We compare the performance of our proposed SUFA model with B-MSFA (De Vito et
al., 2021), BFR (Avalos-Pacheco et al., 2022) and PFA (Roy et al., 2021).
However, PFA did not scale beyond $d>50$ and therefore we could only apply it
for $d=50$. Implementations of all these methods require the user to input an
upper bound on the number of latent factors; we choose the bound following the
strategy in Section 2.2. We repeat BFR and SUFA over 100 replicates. PFA and
B-MSFA codes did not scale across independent parallel replications and we
repeat them over 30 replications.

To assess model fit, we consider the widely applicable Bayesian information
criterion (WBIC, Watanabe, 2013), with lower values implying better fit.
Boxplots of the WBIC values across independent replicates are shown in the top
panel of Figure 1. Notably, posterior samples are required to obtain WBIC
values. The BFR implementation provides the MAP estimate only and we can not
report WBIC values for BFR. For $d=50$, PFA performs slightly better than SUFA
only in the completely misspecifed case. As dimension grows, performance
improvements under SUFA become more evident. Identifiability issues with
B-MSFA may degrade parameter estimation performance but are not expected to
adversely impact goodness-of-fit measures such as WBIC. As SUFA was motivated
by improving identifiability, the better model fit is surprising, and may be
explained by (i) improved learning of the model parameters in high-dimensional
low sample size cases by SUFA due to parsimony and (ii) more efficient
posterior sampling using the proposed HMC-within-Gibbs sampler (itself in part
due to improved parameter identifiability).

We compute the Frobenius norms between the simulation truth and estimated
values of the shared covariance matrix $\bm{\Sigma}$. Boxplots of the norms
across independent replicates are shown in the bottom panel of Figure 1. In
almost every case, SUFA yields the lowest estimation errors.

Figure 1: Comparing our proposed SUFA with B-MSFA, BFR and PFA across
different simulation scenarios (FM1, FM2 and FM3) and model misspecifications:
The top and bottom panels show the boxplots of the WBIC values and the
boxplots of Frobenius norms between the true and estimated shared covariance
$\bm{\Sigma}$, respectively.

Figure 2: The true and estimated shared correlation matrices: The rows
correspond to the three different scenarios considered. In each row, the true
correlation matrix is shown in the leftmost plot; the middle and right plots
are the point estimates in the slight and complete misspecified cases,
respectively.

In Figure 2, we plot the shared correlation matrices
${\mathbf{R}}=\mathrm{diag}(\bm{\Sigma})^{-\frac{1}{2}}\bm{\Sigma}\mathrm{diag}(\bm{\Sigma})^{-\frac{1}{2}}$
recovered by SUFA. In each setting, as a representative of the 100 independent
replicates, we choose the one with the median Frobenius norm from the truth.
We use the posterior mean across the MCMC samples as the point estimate. Due
to page limits, we only show results for $d=200$ in Figure 2. For easier
visualization of the correlation structure, in Figure 2 we use gray color for
correlations having absolute value less then $0.10$. The heatplots under the
“True” panel of Figure 2 show the simulation truths of the three different
correlation structures; the “Slight” and “Complete” misspecification panels
show the recovered ${\mathbf{R}}$ by the SUFA model for the two model
misspecification types under consideration. The figure clearly indicates that
SUFA accurately recovers the shared correlation structures under misspecified
scenarios as well. Additional favorable results on inferring the loadings
$\bm{\Lambda}$ are included in Section S.4.2 of the supplementary materials.
Finally, we compare the total runtimes for $7,500$ MCMC iterations in Figure
3. BFR is fastest since its implementation provides the MAP estimate only
without any uncertainty quantification. PFA is slowest due to requiring
updates of $S$ dense $d\times d$ matrices. However, owing to the distributed
computing implementation of our HMC-within-Gibbs sampler, we observe dramatic
improvements in computational efficiency using SUFA. Complete details of the
parallelized implementation and a complexity analysis showing that each HMC
step is of order $\mathcal{O}(Lqd^{2})$ appear in Section A.2 of the Appendix.

Figure 3: Comparing execution times (in minutes) of B-MSFA, BFR, PFA and SUFA
across increasing dimensions.

###  Application to Gene Expression Data

We now turn to a case study on gene associations among immune cells, which
serve specialized roles in innate and adaptive immune responses that function
to eliminate antigens (Gonzalez et al., 2018). Understanding gene associations
in immune cells is of particular current interest as an essential step in
developing cancer therapeutics in immunotheraphy (Tan et al., 2020). Here, we
integrate data from three studies analyzing gene expressions. The first is the
GSE109125 bulkRNAseq dataset, collected from 103 highly purified immunocyte
populations representing all lineages and several differentiation cascades and
profiled using the ImmGen pipeline (Yoshida et al., 2019). The second study is
a microarray dataset GSE15907 (Painter et al., 2011; Desch et al., 2011),
measured on multiple ex-vivo immune lineages, primarily from adult B6 male
mice. Finally, we include the GSE37448 (Elpek et al., 2014) microarray
dataset, also part of the Immgen project. After standard pre-processing
(detailed in Section S.5.2 of the supplementary material), we work with 474
common genes measured on 156, 146 and 628 cells from the respective datasets.

(a) Shared correlation matrix.

(b) Shared correlation network between genes.

(a) Shared and study-specific loading matrices. Figure 5: Results on ImmGen
data: Panel (a) shows the transpose of the shared loading matrix
$\bm{\Lambda}$ and correlation matrix derived from $\bm{\Sigma}$ from left to
right; panel (b) is the circos plot of the dependency structure between genes,
blue (red) connection implies a positive (negative), their opacities being
proportional to the corresponding association strengths. Associations are only
plotted for absolute correlation $\geq 0.25$; panel (c) plots the top six
columns with highest sum-of-squares values from the shared $\bm{\Lambda}$, and
repeat for the top three columns of each study-specific loading matrix.

We apply SUFA and B-MSFA to the integrated datasets, resulting in WBIC values
of $424278$ and $445916$, respectively. This suggests that SUFA provides a
better fit, and we hence focus on interpreting the SUFA results. First, we
visualize the estimated shared correlation structure between genes in Figure
4(a). We use the uncertainty quantified by the posterior samples to infer
whether correlations are zero—following common practice in the literature on
Bayesian inference on sparsity patterns in covariance matrices (Ksheera Sagar
et al., 2021), we encode the off-diagonals of the correlation matrix as zero
if the respective 95% posterior credible interval contains zero. This leads to
a distinctive sparsity pattern.

Next, we focus on identifying important hub genes that have absolute
correlation $\geq 0.25$ with at least 10 other genes. The resulting dependency
network is summarized as a circos plot (Gu et al., 2014) in Figure 4(b). Our
findings are quite consistent with results from the prior literature: we
observe strong positive correlation between the genes Ly6a, Ly6c1 and Ly6c2
identified within the murine Ly6 complex on chromosome 15 (Lee et al., 2013);
similar behavior is observed within the membrane-spanning 4A class of genes,
namely Ms4a4c, Ms4a4b and Ms4a6b (Liang et al., 2001). Validation of our
findings via a thorough literature review is discussed in Section S.5.3 of the
supplementary materials. The joint analysis of multiple datasets using our
SUFA framework thus reveals interesting findings adding to the relevant
immunology literature, discovering potential relationships that warrant
further scientific study.

Even though gene expression datasets are high-dimensional, often the
variability in the data is governed by relatively few latent factors often
interpreted as meta-genes or pathways (Brunet et al., 2004). It is thus
possible to potentially reduce the dimension of expression data from hundreds
of genes to a handful of meta-genes to identify distinct molecular patterns.
Accordingly the shared $\bm{\Lambda}$ can be interpreted as the common
activation profiles of the meta-genes (Pournara and Wernisch, 2007) across
studies. Differences in gene-networks across studies are captured by the
study-specific latent factors ($\mbox{\boldmath$\zeta$}_{s,i}$’s in equation
(1)); while they are also loaded via the $\bm{\Lambda}$ matrix, the respective
${\mathbf{A}}_{s}$ matrices allow varying profiles.

We visualize the six columns having largest sum-of-squares from the shared
loading matrix, and repeat for the top three columns of each study-specific
loading matrix. These appear in Figure 5(a). We observe large loadings for the
structurally related genes Il18rap and Il18r1 within the Il18 receptor complex
(Parnet et al., 1996) in the first columns consistently across the shared and
study-specific loading matrices. The membrane-spanning 4A class of genes
(Ms4a4c, Ms4a4b and Ms4a6b) exhibit similar behavior (see column 3 of the
shared loading in Figure 5(a)). These results indicate potential existence of
meta-genes controlling the co-expression profiles of genes in the same
pathway. Additional discussions on interpretation of factor loadings are
provided in Section S.5.4 of the supplementary materials. Coherent with
empirical studies (Wang et al., 2014), the strongest signals appear in the
bulkRNAseq data despite having much smaller sample size than the microarrays.
The loadings of Microarray 2 show much stronger signals than those in
Microarray 1, possibly due to the larger sample size of the former. These
intuitive findings, corroborated by results from the prior literature,
illustrate the interpretability of the inferred study-specific loadings under
SUFA.

##  Discussion

This article proposes a new factor analytic approach for covariance learning
by integrating multiple related studies in high-dimensional settings. We
provide both practical solutions and theoretical guarantees for learning the
shared and study-specific covariance structures, ensuring generic
identifiability. Moreover, we quantify the utility of data integration under
the proposed model by way of improved learning rates with each additional
study in our analysis of posterior contraction properties. The theoretical
investigations on posterior concentration provide support for the study-
specific inferences as well. Finally, these contributions are developed
alongside a scalable HMC approach to posterior inference capable of handing
high-dimensional datasets with massive sample size. In realistic simulation
studies, we show that our proposed method outperforms several existing
approaches in multiple aspects. Our approach yields new insights when applied
to a gene network estimation problem, integrating immune cell gene expression
datasets.

Leveraging low-rank matrix decompositions, factor models have widespread scope
in contemporary data analysis regimes beyond our focus on estimating
covariance structure. These extensions include measurement error problems
(Sarkar et al., 2021), model-based clustering (Chandra et al., 2023),
conditional graph estimation (Chandra et al., 2021), among many others. Our
proposed SUFA approach can be readily adopted in these domains. Integrating
data from multiple sources is also of interest in many applications involving
non-Gaussian data. Building on copula-based factor models (Murray et al.,
2013), our approach can be hierarchically adopted to handle such data. When
additional information is available regarding groups of features, it can be of
interest to use group sparsity priors, potentially building on the approach of
Schiavon et al. (2022) for more refined inference on covariance structures.
These extensions represent several examples of promising future directions
building upon the ideas in this article.

### Supplementary Materials

Complete prior specifications, proofs and details of the theoretical results
and technical lemmas, extended simulation results and further details on the
gene expression datasets.

### Acknowledgments

The authors are grateful for partial funding support from NIH grants
R01-ES027498 and R01-ES028804, NSF grants DMS-2230074 and PIPP-2200047, and
ERC Horizon 2020 grant agreement No. 856506.

## Appendix

###  Proofs of Main Results

###### Proof of Lemma 1.

We first show that information switching or the existence of matrices
${\mathbf{C}}_{1}$ and ${\mathbf{C}}_{2}$ as mentioned in (I) implies
$\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ is non-null. We let
${\mathbf{C}}={\mathbf{C}}_{1}+{\mathbf{C}}_{2}$. Since ${\mathbf{C}}$ is
symmetric, ${\mathbf{C}}={\mathbf{H}}{\mathbf{D}}{\mathbf{H}}^{\rm T}$ where
${\mathbf{H}}=\begin{bmatrix}{\mathbf{h}}_{1}&\cdots&{\mathbf{h}}_{q}\end{bmatrix}$
is an orthogonal matrix and ${\mathbf{D}}=\mathrm{diag}(d_{1},\dots,d_{q})$ is
a diagonal matrix with real entries. Clearly
${\mathbf{C}}=\sum_{i=1}^{q}d_{i}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{\rm T}$.
Letting
${\mathbf{C}}^{+}=\sum_{i=1}^{q}d_{i}\mathbbm{1}_{d_{i}>0}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{\rm
T}$ and
${\mathbf{C}}^{-}=-\sum_{i=1}^{q}d_{i}\mathbbm{1}_{d_{i}<0}{\mathbf{h}}_{i}{\mathbf{h}}_{i}^{\rm
T}$, we have ${\mathbf{C}}={\mathbf{C}}^{+}-{\mathbf{C}}^{-}$. We establish
the result in a case-by-case manner:

Case 1: ${{\mathbf{C}}}^{+}={\mathbf{0}}$

Let ${\mathbf{x}}$ be a non-zero vector such that ${\mathbf{x}}^{\rm
T}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}{\mathbf{x}}=0$. Since by assumption
${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}+{\mathbf{C}}$ is positive semi-
definite (psd), ${\mathbf{x}}^{\rm T}({\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}+{\mathbf{C}}){\mathbf{x}}\geq 0$. Now ${\mathbf{x}}^{\rm
T}({\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}+{\mathbf{C}}){\mathbf{x}}=-{\mathbf{x}}^{\rm
T}{\mathbf{C}}^{-}{\mathbf{x}}\geq 0$ implies that ${\mathbf{x}}^{\rm
T}{\mathbf{C}}^{-}{\mathbf{x}}=0$ as ${\mathbf{C}}^{-}$ is a psd matrix by
construction. Now for any psd matrix ${\mathbf{B}}$ and non-zero vector
${\mathbf{x}}$, ${\mathbf{x}}^{\rm T}{\mathbf{B}}{\mathbf{x}}=0$ implies that
${\mathbf{x}}\in\mathbb{N}\left({\mathbf{B}}\right)$ where
$\mathbb{N}\left({\mathbf{B}}\right)$ denotes the null space of
$\mathbb{C}\left({\mathbf{B}}\right)$. Hence
$\mathbb{N}\left({\mathbf{A}}_{s}\right)\subset\mathbb{N}\left({{\mathbf{C}}^{-}}\right)$
for all $s=1,\dots,S$ implying that
$\mathbb{C}\left({\mathbf{C}}^{-}\right)\subseteq\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$.

Case 2: ${{\mathbf{C}}}^{-}={\mathbf{0}}$

Due to information switching ${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}+{\mathbf{C}}=\widetilde{{\mathbf{A}}}_{s}\widetilde{{\mathbf{A}}}_{s}^{\rm
T}$ for $q\times q_{s}$ matrix $\widetilde{{\mathbf{A}}}_{s}$ for all $s$. Now
${\mathbf{C}}={\mathbf{C}}^{+}\Rightarrow\mathbb{C}\left({\mathbf{A}}_{s}\right)\subset\mathbb{C}\left(\widetilde{{\mathbf{A}}}_{s}\right)$.
Note that $\mathrm{rank}(\widetilde{{\mathbf{A}}}_{s})\leq q_{s}$. Hence

$\mathrm{rank}({\mathbf{A}}_{s})=q_{s}\Rightarrow\mathbb{C}\left({\mathbf{A}}_{s}\right)=\mathbb{C}\left(\widetilde{{\mathbf{A}}}_{s}\right)\Rightarrow\mathbb{C}\left({\mathbf{C}}^{+}\right)\subset\mathbb{C}\left({\mathbf{A}}_{s}\right)\Rightarrow\cap_{s=1}^{s}\mathbb{C}\left({\mathbf{A}}_{s}\right)\supset\mathbb{C}\left({\mathbf{C}}^{+}\right)\neq\phi.$

Case 3: ${{\mathbf{C}}}^{+},{{\mathbf{C}}}^{-}\neq{\mathbf{0}}$

If possible let $\cap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)=\phi$.
This assumption implies that
$\cup_{s=1}^{S}\mathbb{N}\left({\mathbf{A}}_{s}\right)=\mathbb{R}^{q}$. Note
that ${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}+{\mathbf{C}}=\widetilde{{\mathbf{A}}}_{s}\widetilde{{\mathbf{A}}}_{s}^{\rm
T}+{\mathbf{C}}^{-}$ where
$\widetilde{{\mathbf{A}}}_{s}=\begin{bmatrix}{\mathbf{A}}_{s}&{\mathbf{C}}^{+\frac{1}{2}}\end{bmatrix}$.
Hence
$\mathbb{C}\left({\mathbf{A}}_{s}\right)\subset\mathbb{C}\left(\widetilde{{\mathbf{A}}}_{s}\right)$.
Using the argument from case 1 we have
$\mathbb{C}\left({\mathbf{C}}^{-}\right)\subset\mathbb{C}\left(\widetilde{{\mathbf{A}}}_{s}\right)$.
Since ${\mathbf{C}}^{-}$ is non-null, we can find a non-zero vector
${\mathbf{x}}\in\mathbb{C}\left({\mathbf{C}}^{-}\right)$. Note that by
construction $\mathbb{C}\left({\mathbf{C}}^{-}\right)$ and
$\mathbb{C}\left({\mathbf{C}}^{+}\right)$ are orthogonal and hence
${\mathbf{x}}^{\rm T}{\mathbf{C}}^{+}{\mathbf{x}}=0$. Also
${\mathbf{x}}\in\mathbb{N}\left({\mathbf{A}}_{\widetilde{s}}\right)$ for some
$\widetilde{s}$ as
$\cup_{s=1}^{S}\mathbb{N}\left({\mathbf{A}}_{s}\right)=\mathbb{R}^{q}$
implying that ${\mathbf{x}}^{\rm
T}{\mathbf{A}}_{\widetilde{s}}{\mathbf{A}}_{\widetilde{s}}^{\rm
T}{\mathbf{x}}=0$. Now
${\mathbf{A}}_{\widetilde{s}}{\mathbf{A}}_{\widetilde{s}}^{\rm
T}+{\mathbf{C}}$ is a psd matrix implies that ${\mathbf{x}}^{\rm
T}\left({\mathbf{A}}_{\widetilde{s}}{\mathbf{A}}_{\widetilde{s}}^{\rm
T}+{\mathbf{C}}\right){\mathbf{x}}\geq 0$. By definition of ${\mathbf{x}}$

${\mathbf{x}}^{\rm
T}\left({\mathbf{A}}_{\widetilde{s}}{\mathbf{A}}_{\widetilde{s}}^{\rm
T}+{\mathbf{C}}\right){\mathbf{x}}={\mathbf{x}}^{\rm
T}{\mathbf{A}}_{\widetilde{s}}{\mathbf{A}}_{\widetilde{s}}^{\rm
T}{\mathbf{x}}+{\mathbf{x}}^{\rm
T}{\mathbf{C}}^{+}{\mathbf{x}}-{\mathbf{x}}^{\rm
T}{\mathbf{C}}^{-}{\mathbf{x}}=-{\mathbf{x}}^{\rm
T}{\mathbf{C}}^{-}{\mathbf{x}}<0.$

Since ${\mathbf{x}}\in\mathbb{C}\left({\mathbf{C}}^{-}\right)$,
${\mathbf{x}}^{\rm T}{\mathbf{C}}^{-}{\mathbf{x}}>0$ leading to the
contradiction in the above equation. Hence
$\cap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ can not be null.

Next we show that a non-null
$\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ guarantees the
existence of a (not necessarily unique) psd matrix ${\mathbf{C}}$ such that
${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}-{\mathbf{C}}\succeq{\mathbf{0}}$ for
all $s$ resulting in information switching. Let
${\mathbf{H}}=\begin{bmatrix}{\mathbf{h}}_{1}&\cdots&{\mathbf{h}}_{r}\end{bmatrix}$
be a $q\times r$ matrix with orthonormal columns such that the columns form a
basis of $\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$. Let
${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}={\mathbf{H}}_{s}{\mathbf{D}}_{s}{\mathbf{H}}_{s}^{\rm T}$ be the spectral
decomposition of ${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}$. Clearly
${\mathbf{D}}_{s}$ are diagonal matrices with non-negative entries. Let
$\varsigma$ be the minimum among all positive entries of ${\mathbf{D}}_{s}$
across all $s=1,\dots,S$. Define
${\mathbf{C}}=\varsigma{\mathbf{H}}{\mathbf{I}}_{r}{\mathbf{H}}^{\rm T}$. By
construction ${\mathbf{x}}^{\rm T}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm
T}{\mathbf{x}}=0$ implies that ${\mathbf{x}}^{\rm
T}{\mathbf{C}}{\mathbf{x}}=0$. Also for
${\mathbf{x}}\in\mathbb{C}\left({\mathbf{C}}\right)$, we have
$0<{\mathbf{x}}^{\rm
T}{\mathbf{C}}{\mathbf{x}}\leq\varsigma\norm{{\mathbf{x}}}^{2}\leq{\mathbf{x}}^{\rm
T}{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}{\mathbf{x}}$. Therefore,
${\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}-{\mathbf{C}}\succeq{\mathbf{0}}$ for
all $s=1,\dots,S$, establishing the result. ∎

###### Proof of Theorem 1.

when information switching occurs, Lemma 1 guarantees that
$\bigcap_{s=1}^{S}\mathbb{C}\left({\mathbf{A}}_{s}\right)$ is non-null, so
there exist vectors ${\mathbf{x}}_{1},\dots,{\mathbf{x}}_{S}$ such that
${\mathbf{A}}_{1}{\mathbf{x}}_{1}=\dots={\mathbf{A}}_{S}{\mathbf{x}}_{S}\neq{\mathbf{0}}$.
Now, for any non-degenerate continuous prior on
${\mathbf{A}}_{s},s=1,\dots,S$, the augmented matrix
$\begin{bmatrix}{\mathbf{A}}_{1}&\cdots&{\mathbf{A}}_{S}\end{bmatrix}$ is full
rank if and only if $\sum_{s=1}^{S}q_{s}\leq q$. Hence, the necessary and
sufficient condition in Lemma 1 has zero prior support. ∎

###### Proof of Theorem 2.

For two sequences $a_{n},b_{n}\geq 0$ we let $a_{n}\asymp b_{n}$ imply
$0<\liminf\absolutevalue{{a_{n}}/{b_{n}}}\leq\limsup\absolutevalue{{a_{n}}/{b_{n}}}<\infty$.
We define
$\mbox{\boldmath$\Theta$}_{n}=\\{\bm{\Lambda}_{n},\mbox{\boldmath$\Delta$}_{n},{\mathbf{A}}_{1n},\dots,{\mathbf{A}}_{Sn}\\}$
as the set of all parameters, and let $\mathbb{P}_{n}$ and $\mathbb{P}_{0n}$
denote the joint distributions of $\mathcal{D}_{n}$ under
$\mbox{\boldmath$\Theta$}_{n}$ and the true value
$\mbox{\boldmath$\Theta$}_{0n}$, respectively. to be the Kullback-Leibler (KL)
divergence between $\mathbb{P}_{0n}$ and $\mathbb{P}_{n}$.

Let $\mathcal{P}_{n}=\mathcal{P}_{1,n}\cup\mathcal{P}_{2,n}$ be a partition of
the parameter space $\mathcal{P}_{n}$ where
$\mathcal{P}_{1,n}=\\{\mbox{\boldmath$\Theta$}_{n}:\norm{{\mathbf{A}}_{sn}}_{2}^{2}\leq
Cn\tau_{n}^{2}\text{ for all }s=1,\dots,S\\}$ for some large enough constant
$C>0$ with $n\tau_{n}^{2}=c_{n}^{2}q_{0n}s_{n}\log(d_{n}q_{n})$ and
$\mathcal{P}_{2,n}=\mathcal{P}_{1,n}^{\complement}$. Note that
$\Pi_{n}\left(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n}\right)\leq\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})+\Pi_{n}\left(\mathcal{P}_{2,n}\mid\mathcal{D}_{n}\right)$.

To prove the theorem, we show that, as $n\to\infty$, the posterior
distribution on $\mathcal{P}_{1,n}$ concentrates around $\bm{\Sigma}_{0n}$
while the remaining mass assigned to $\mathcal{P}_{2,n}$ diminishes to zero.
We define
$B_{n,0}(\mbox{\boldmath$\Theta$}_{0n},\epsilon)=\left\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{n}:\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq
n\epsilon^{2}\right\\}$ as the set of probability measures characterized by
$\mbox{\boldmath$\Theta$}_{n}$ within $\epsilon$-radius of $\mathbb{P}_{0n}$
in KL divergence. We state the following lemma.

######  Lemma 2.

Let $e_{n}$ be a semimetric on $\mathcal{P}_{n}$ such that
$e_{n}(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n})=\frac{1}{n\tau_{n}^{2}c_{n}^{3}}(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}-c_{n}^{2}\epsilon_{n})$
where $\epsilon_{n}\asymp c_{n}n\tau_{n}^{3}$. Then
$\Pi_{n}\left\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:e_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)>M\tau_{n}\mid\mathcal{D}_{n}\right\\}\to
0$ in $\mathbb{P}_{0n}$-probability for a sufficiently large $M$ as long as
the following conditions hold for sufficiently large $j\in\mathbb{N}$:

  1. (I)

For some $C>0$,
$\Pi_{n}\left\\{B_{n,0}(\mbox{\boldmath$\Theta$}_{0n},\tau_{n})\right\\}\geq
e^{-Cn\tau_{n}^{2}}$.

  2. (II)

Define the set
$G_{j,n}=\left\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}<e_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)\leq
2j\tau_{n}\right\\}$. There exist test functions $\varphi_{n}$ such that for
some $K>0$,
$\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{0n}}\varphi_{n}=0$ and
$\sup_{\mbox{\boldmath$\Theta$}_{n}\in
G_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{n})\leq\exp(-Knj^{2}\tau_{n}^{2})$.

The remaining technical arguments are included in full in the Supplemental
Materials. We prove Lemma 2 in Section S.2.1, and verify condition (I) in
Theorem 3 and show the existence of a sequence of test functions satisfying
condition (II) in Theorem 4 there. Note that
$e_{n}\left(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n}\right)>M\tau_{n}\Leftrightarrow\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}$.
This implies that
$\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\to
0$ in $\mathbb{P}_{0n}$-probability. Subsequently applying the dominated
convergence theorem (DCT) establishes that
$\lim_{n\to\infty}\mathbb{E}_{\mathbb{P}_{0n}}\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})=0$.
It remains to show that the remaining mass assigned to $\mathcal{P}_{2,n}$
goes to 0, which is established in Theorem 5. ∎

###### Proof of Corollary 1.

Let ${\mathbf{F}}_{sn}=\bm{\Lambda}_{n}{\mathbf{A}}_{sn}{\mathbf{A}}_{sn}^{\rm
T}\bm{\Lambda}_{n}^{\rm T}$ and
${\mathbf{F}}_{0sn}=\bm{\Lambda}_{0n}{\mathbf{A}}_{0sn}{\mathbf{A}}_{0sn}^{\rm
T}\bm{\Lambda}_{0n}^{\rm T}$. Using the notations defined in the proof of
Theorem 2, we have
$\Pi_{n}(\norm{{\mathbf{F}}_{sn}-{\mathbf{F}}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\leq\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{{\mathbf{F}}_{sn}-{\mathbf{F}}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})+\Pi_{n}\left(\mathcal{P}_{2,n}\mid\mathcal{D}_{n}\right)$.
As Theorem 5 in the supplementary materials shows that
$\Pi_{n}\left(\mathcal{P}_{2,n}\mid\mathcal{D}_{n}\right)\to 0$ in
$\mathbb{P}_{0n}$-probability we focus on the first expression in the previous
display in this proof.

Note that $\bm{\Sigma}_{0sn}=\bm{\Sigma}_{0n}+{\mathbf{F}}_{0sn}$ and
$\bm{\Sigma}_{sn}=\bm{\Sigma}_{n}+{\mathbf{F}}_{sn}$, and therefore
$\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{{\mathbf{F}}_{sn}-{\mathbf{F}}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\leq\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}>\frac{M}{2}\varepsilon_{n}\mid\mathcal{D}_{n})+\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}>\frac{M}{2}\varepsilon_{n}\mid\mathcal{D}_{n})$.
In Theorem 2 we established that the first term in the previous display
diminishes to 0 in $\mathbb{P}_{0n}$-probability; in Theorem 6 of the
supplementary materials we show that
$\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}>\frac{M}{2}\varepsilon_{n}\mid\mathcal{D}_{n})\to
0$ in $\mathbb{P}_{0n}$-probability. Hence
$\Pi_{n}(\norm{{\mathbf{F}}_{sn}-{\mathbf{F}}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\to
0$ in $\mathbb{P}_{0n}$-probability. Subsequently applying DCT we conclude the
proof. ∎

###  Sampler Details and Distributed Computation

In this section, instead of $\bm{\Delta}$ we use its one-to-one transformation
$\widetilde{\bm{\delta}}=(\widetilde{\delta}_{1},\dots,\widetilde{\delta}_{d})^{\rm
T}=\log\mathrm{diag}(\bm{\Delta})$. It allows updating the
$\widetilde{\bm{\delta}}$ vector in the unrestricted space $\mathbb{R}^{d}$
resulting in simplified numerical operations. Accordingly, we redefine
$\mbox{\boldmath$\Theta$}:=(\bm{\Lambda},\widetilde{\bm{\delta}},{\mathbf{A}}_{1},\dots,{\mathbf{A}}_{S})$.

##### Gradients of the $\log$-posterior:

Following equation (8), we have

$\textstyle{\log\Pi(\mbox{\boldmath$\Theta$}\mid-)={\cal
L}+\log\Pi(\bm{\Lambda}\mid\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$},\tau)+\log\Pi(\widetilde{\bm{\delta}})+\sum_{s=1}^{S}\log\Pi({\mathbf{A}}_{s})},$

where
$\log\Pi(\bm{\Lambda}\mid\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$},\tau)=K_{\bm{\Lambda}}-\frac{1}{2\tau^{2}}\sum_{j=1}^{d}\sum_{h=1}^{q}\frac{\lambda_{j,h}^{2}}{\psi_{j,h}\phi_{j,h}^{2}}$,
$\log\Pi(\widetilde{\bm{\delta}})=K_{\widetilde{\bm{\delta}}}-\frac{1}{2\sigma^{2}_{\delta}}\sum_{j=1}^{d}(\widetilde{\delta}_{j}-\mu_{\delta})^{2}$
and
$\log\Pi({\mathbf{A}}_{s})=K_{{\mathbf{A}}_{s}}-\frac{1}{2b_{{\mathbf{A}}}}\sum_{j=1}^{q}\sum_{h=1}^{q_{s}}a_{s,j,h}^{2}$.
The constants $K_{\bm{\Lambda}}$, $K_{\widetilde{\bm{\delta}}}$ and
$K_{{\mathbf{A}}_{s}}$’s are never required as they cancel out in the
Metropolis-Hastings ratio. Analytical expressions for the partials of
$\log\Pi(\bm{\Lambda}\mid-)$ are given below:

$\displaystyle\textstyle{\partialderivative{}{\bm{\Lambda}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)=-\sum_{s=1}^{S}{\mathbf{G}}_{s}\bm{\Lambda}{\mathbf{C}}_{s}-\left(\frac{\lambda_{j,h}}{\psi_{j,h}\phi_{j,h}^{2}\tau^{2}}\right)_{d\times
q}},$ (9)
$\displaystyle\textstyle{\partialderivative{}{\widetilde{\bm{\delta}}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)=-\frac{1}{2}\sum_{s=1}^{S}\mathrm{diag}({\mathbf{G}}_{s})\odot\exp(\widetilde{\bm{\delta}})-\left(\frac{\widetilde{\delta}_{1}-\mu_{\delta}}{\sigma^{2}_{\delta}},\dots,\frac{\widetilde{\delta}_{d}-\mu_{\delta}}{\sigma^{2}_{\delta}}\right)^{\rm
T}},$ (10)
$\displaystyle\textstyle{\partialderivative{}{{\mathbf{A}}_{s}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)=-\bm{\Lambda}^{\rm
T}{\mathbf{G}}_{s}\bm{\Lambda}{\mathbf{A}}_{s}-\frac{1}{b_{{\mathbf{A}}}}{\mathbf{A}}_{s}},$
(11)

where
${\mathbf{G}}_{s}=n_{s}\bm{\Sigma}_{s}^{-1}-\bm{\Sigma}_{s}^{-1}{\mathbf{W}}_{s}\bm{\Sigma}_{s}^{-1}$,
${\mathbf{C}}_{s}={\mathbf{I}}_{q}+{\mathbf{A}}_{s}{\mathbf{A}}_{s}^{\rm T}$
and $\odot$ denotes the Hadamard product between two vectors/matrices.

##### Distributed gradient and likelihood computation:

In each HMC step of Algorithm 1, calculating
$\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$ and its gradients
$\partialderivative{\mbox{\boldmath$\Theta$}}\log\Pi(\mbox{\boldmath$\Theta$}\mid-)$
are required for $L$ many leapfrog operations. Note that the first terms in
the RHS of (9)-(11) are gradients of ${\cal L}$ with respect to parameters,
while the second terms are gradients of the prior densities. Although
computing the gradients of ${\cal L}$ are seemingly intensive, these can be
obtained very efficiently by dividing the calculations into parallel sub-
operations across different studies. These sub-operations can be distributed
over parallel processes in multi-thread computing systems. We elaborate the
procedure in Algorithm 2, and illustrate it in a schematic diagram in Figure
6. The figure sketches out a flowchart to avoid repetitive calculations by
storing variables and subsequently using them in following steps. The next
lemma establishes its complexity.

1 parallel for  _$s=1,\dots,S$_ do

2 Compute & store$~{}{\mathbf{C}}_{s}$, ${\mathbf{C}}_{s}^{\frac{1}{2}}$,
$\widetilde{\bm{\Lambda}}_{s}=\bm{\Lambda}{\mathbf{C}}_{s}^{\frac{1}{2}}$
where ${\mathbf{C}}_{s}^{\frac{1}{2}}=\texttt{Cholesky}({\mathbf{C}}_{s})$ and
$\bm{\Sigma}_{s}^{-1}=(\widetilde{\bm{\Lambda}}_{s}\widetilde{\bm{\Lambda}}_{s}^{\rm
T}+\bm{\Delta})^{-1}=\bm{\Delta}^{-1}$
$-\bm{\Delta}^{-1}\widetilde{\bm{\Lambda}}_{s}({\mathbf{I}}_{q}+\widetilde{\bm{\Lambda}}_{s}^{\rm
T}\bm{\Delta}^{-1}\widetilde{\bm{\Lambda}}_{s})^{-1}\widetilde{\bm{\Lambda}}_{s}^{\rm
T}\bm{\Delta}^{-1}$.

3 Compute & store$~{}\bm{\Sigma}_{s}^{^{-1}}{\mathbf{W}}_{s}$ and subsequently
${\mathbf{G}}_{s}$.

4 Compute & store$~{}{\mathbf{G}}_{s}\bm{\Lambda}$ and
${\mathbf{G}}_{s}\bm{\Lambda}{\mathbf{C}}_{s}$.

5 Compute $\bm{\Lambda}^{\rm T}{\mathbf{G}}_{s}\bm{\Lambda}{\mathbf{A}}_{s}$
to obtain $\partialderivative{{\cal L}}{{\mathbf{A}}_{s}}$.

6 end for

7Sum the vectors $\mathrm{diag}({\mathbf{G}}_{s})$ from Step 2 across
$s=1,\dots,S$ to obtain $\partialderivative{{\cal
L}}{\widetilde{\bm{\delta}}}$.

Sum the matrices ${\mathbf{G}}_{s}\bm{\Lambda}{\mathbf{C}}_{s}$ from Step 2
across $s=1,\dots,S$ to obtain $\partialderivative{{\cal L}}{\bm{\Lambda}}$.

Algorithm 2 Distributed Gradient Computation of the $\log$-likelihood ${\cal
L}$

######  Lemma 3.

The runtime complexity of each HMC step is $\mathcal{O}(Lqd^{2})$ with $L$
being the number of leapfrog steps.

Figure 6: Schematic diagram of the distributed gradient computation.

## References

  * Allman et al. (2009) Allman, E. S., Matias, C., and Rhodes, J. A. (2009).  Identifiability of parameters in latent structure models with many observed variables.  The Annals of Statistics, 37, 3099–3132.
  * Armagan et al. (2013) Armagan, A., Dunson, D. B., and Lee, J. (2013).  Generalized double Pareto shrinkage.  Statistica Sinica, 23, 119–143.
  * Avalos-Pacheco et al. (2022) Avalos-Pacheco, A., Rossell, D., and Savage, R. S. (2022).  Heterogeneous large datasets integration using Bayesian factor regression.  Bayesian Analysis, 17, 33–66.
  * Baglama and Reichel (2005) Baglama, J. and Reichel, L. (2005).  Augmented implicitly restarted Lanczos bidiagonalization methods.  SIAM Journal on Scientific Computing, 27, 19–42.
  * Bekker and ten Berge (1997) Bekker, P. A. and ten Berge, J. M. (1997).  Generic global indentification in factor analysis.  Linear Algebra and its Applications, 264, 255–263.
  * Bhattacharya and Dunson (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse Bayesian infinite factor models.  Biometrika, 98, 291–306.
  * Bhattacharya et al. (2015) Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B. (2015).  Dirichlet-Laplace priors for optimal shrinkage.  Journal of the American Statistical Association, 110, 1479–1490.
  * Brunet et al. (2004) Brunet, J.-P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2004).  Metagenes and molecular pattern discovery using matrix factorization.  Proceedings of the National Academy of Sciences, 101, 4164–4169.
  * Carpenter et al. (2015) Carpenter, B., Hoffman, M. D., Brubaker, M. A., Lee, D. D., Li, P., and Betancourt, M. (2015).  The Stan math library: Reverse-mode automatic differentiation in C++.  arXiv 1509.07164.
  * Carvalho et al. (2008) Carvalho, C. M., Chang, J., Lucas, J. E., Nevins, J. R., Wang, Q., and West, M. (2008).  High-dimensional sparse factor modeling: Applications in gene expression genomics.  Journal of the American Statistical Association, 103, 1438–1456.
  * Carvalho et al. (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In Artificial Intelligence and Statistics, pages 73–80. PMLR.
  * Chandra et al. (2021) Chandra, N. K., Müller, P., and Sarkar, A. (2021).  Bayesian precision factor analysis for high-dimensional sparse Gaussian graphical models.  arXiv:2107.11316.
  * Chandra et al. (2023) Chandra, N. K., Canale, A., and Dunson, D. B. (2023).  Escaping the curse of dimensionality in Bayesian model-based clustering.  Journal of Machine Learning Research.  To appear.
  * Dai et al. (2020) Dai, F., Dutta, S., and Maitra, R. (2020).  A matrix-free likelihood method for exploratory factor analysis of high-dimensional Gaussian data.  Journal of Computational and Graphical Statistics, 29, 675–680.  PMID: 33041614.
  * Daniele et al. (2019) Daniele, M., Pohlmeier, W., and Zagidullina, A. (2019).  Sparse approximate factor estimation for high-dimensional covariance matrices.  arXiv:1906.05545.
  * De Vito et al. (2019) De Vito, R., Bellio, R., Trippa, L., and Parmigiani, G. (2019).  Multi-study factor analysis.  Biometrics, 75, 337–346.
  * De Vito et al. (2021) De Vito, R., Bellio, R., Trippa, L., and Parmigiani, G. (2021).  Bayesian multistudy factor analysis for high-throughput biological data.  The Annals of Applied Statistics, 15, 1723–1741.
  * Desch et al. (2011) Desch, A. N., Randolph, G. J., Murphy, K., et al. (2011).  CD103+ pulmonary dendritic cells preferentially acquire and present apoptotic cell–associated antigen.  Journal of Experimental Medicine, 208, 1789–1797.
  * Elpek et al. (2014) Elpek, K. G., Cremasco, V., Shen, H., et al. (2014).  The tumor microenvironment shapes lineage, transcriptional, and functional diversity of infiltrating myeloid cells.  Cancer Immunology Research, 2, 655–667.
  * Erosheva and Curtis (2017) Erosheva, E. A. and Curtis, S. M. (2017).  Dealing with reflection invariance in Bayesian factor analysis.  Psychometrika, 82, 295–307.
  * Fan et al. (2008) Fan, J., Fan, Y., and Lv, J. (2008).  High dimensional covariance matrix estimation using a factor model.  Journal of Econometrics, 147, 186–197.
  * Franks and Hoff (2019) Franks, A. M. and Hoff, P. (2019).  Shared subspace models for multi-group covariance estimation.  Journal of Machine Learning Research, 20, 1–37.
  * Frühwirth-Schnatter and Lopes (2018) Frühwirth-Schnatter, S. and Lopes, H. F. (2018).  Sparse Bayesian factor analysis when the number of factors is unknown.  arXiv preprint arXiv:1804.04231.
  * Gelman (2006) Gelman, A. (2006).  Prior distributions for variance parameters in hierarchical models.  Bayesian Analysis, 1, 515–534.
  * Gonzalez et al. (2018) Gonzalez, H., Hagerling, C., and Werb, Z. (2018).  Roles of the immune system in cancer: From tumor initiation to metastatic progression.  Genes & Development, 32, 1267–1284.
  * Grabski et al. (2023) Grabski, I. N., Vito, R. D., Trippa, L., and Parmigiani, G. (2023).  Bayesian combinatorial multistudy factor analysis.  The Annals of Applied Statistics, 17, 2212–2235.
  * Green (1995) Green, P. J. (1995).  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika, 82, 711–732.
  * Gu et al. (2014) Gu, Z., Gu, L., Eils, R., Schlesner, M., and Brors, B. (2014).  circlize implements and enhances circular visualization in R.  Bioinformatics, 30, 2811–2812.
  * Heng et al. (2008) Heng, T. S., Painter, M. W., et al. (2008).  The immunological genome project: networks of gene expression in immune cells.  Nature Immunology, 9, 1091–1094.
  * Iacob et al. (2016) Iacob, E., Light, A. R., Donaldson, G. W., et al. (2016).  Gene expression factor analysis to differentiate pathways linked to fibromyalgia, chronic fatigue syndrome, and depression in a diverse patient sample.  Arthritis Care & Research, 68, 132–140.
  * Ishwaran and Rao (2005) Ishwaran, H. and Rao, J. S. (2005).  Spike and slab variable selection: Frequentist and Bayesian strategies.  The Annals of Statistics, 33, 730–773.
  * Kaiser (1958) Kaiser, H. F. (1958).  The varimax criterion for analytic rotation in factor analysis.  Psychometrika, 23, 187–200.
  * Kepler et al. (2002) Kepler, T. B., Crosby, L., and Morgan, K. T. (2002).  Normalization and analysis of DNA microarray data by self-consistency and local regression.  Genome Biology, 3, 1–12.
  * Knowles and Ghahramani (2011) Knowles, D. and Ghahramani, Z. (2011).  Nonparametric Bayesian sparse factor models with application to gene expression modeling.  The Annals of Applied Statistics, 5, 1534–1552.
  * Ksheera Sagar et al. (2021) Ksheera Sagar, K. N., Banerjee, S., Datta, J., and Bhadra, A. (2021).  Precision matrix estimation under the horseshoe-like prior-penalty dual.  arXiv:2104.10750.
  * Lee et al. (2013) Lee, P. Y., Wang, J.-X., Parisini, E., Dascher, C. C., and Nigrovic, P. A. (2013).  Ly6 family proteins in neutrophil biology.  Journal of Leukocyte Biology, 94, 585–594.
  * Legramanti et al. (2020) Legramanti, S., Durante, D., and Dunson, D. B. (2020).  Bayesian cumulative shrinkage for infinite factorizations.  Biometrika, 107, 745–752.
  * Liang et al. (2001) Liang, Y., Buckley, T. R., et al. (2001).  Structural organization of the human MS4A gene cluster on chromosome 11q12.  Immunogenetics, 53, 357–368.
  * Millsap (2001) Millsap, R. E. (2001).  When trivial constraints are not trivial: The choice of uniqueness constraints in confirmatory factor analysis.  Structural Equation Modeling: A Multidisciplinary Journal, 8, 1–17.
  * Murray et al. (2013) Murray, J. S., Dunson, D. B., Carin, L., and Lucas, J. E. (2013).  Bayesian Gaussian copula factor models for mixed data.  Journal of the American Statistical Association, 108, 656–665.
  * Neal (2011) Neal, R. (2011).  MCMC using Hamiltonian dynamics. In Chapter 5 of the Handbook of Markov Chain Monte Carlo Edited by Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng.
  * Painter et al. (2011) Painter, M. W., Davis, S., Hardy, R. R., et al. (2011).  Transcriptomes of the B and T lineages compared by multiplatform microarray profiling.  The Journal of Immunology, 186, 3047–3057.
  * Papastamoulis and Ntzoufras (2022) Papastamoulis, P. and Ntzoufras, I. (2022).  On the identifiability of bayesian factor analytic models.  Statistics and Computing, 32, 23.
  * Parnet et al. (1996) Parnet, P., Garka, K. E., et al. (1996).  IL-1Rrp is a novel receptor-like molecule similar to the type I Interleukin-1 receptor and its homologues T1/ST2 and IL-1R AcP.  Journal of Biological Chemistry, 271, 3967–3970.
  * Pati et al. (2014) Pati, D., Bhattacharya, A., Pillai, N. S., and Dunson, D. (2014).  Posterior contraction in sparse Bayesian factor models for massive covariance matrices.  The Annals of Statistics, 42, 1102–1130.
  * Pournara and Wernisch (2007) Pournara, I. and Wernisch, L. (2007).  Factor analysis for gene regulatory networks and transcription factor activity profiles.  BMC Bioinformatics, 8, 1–20.
  * Poworoznek et al. (2021) Poworoznek, E., Ferrari, F., and Dunson, D. (2021).  Efficiently resolving rotational ambiguity in Bayesian matrix sampling with matching.  arXiv:2107.13783.
  * Robert and Roberts (2021) Robert, C. P. and Roberts, G. (2021).  Rao–Blackwellisation in the Markov chain Monte Carlo era.  International Statistical Review, 89, 237–249.
  * Ročková and George (2016) Ročková, V. and George, E. I. (2016).  Fast Bayesian factor analysis via automatic rotations to sparsity.  Journal of the American Statistical Association, 111, 1608–1622.
  * Rohe and Zeng (2020) Rohe, K. and Zeng, M. (2020).  Vintage factor analysis with varimax performs statistical inference.  Journal of the Royal Statistical Society: Series B (Statistical Methodology).  To appear.
  * Roy et al. (2021) Roy, A., Lavine, I., Herring, A. H., and Dunson, D. B. (2021).  Perturbed factor analysis: Accounting for group differences in exposure profiles.  The Annals of Applied Statistics, 15, 1386–1404.
  * Russell (2002) Russell, D. W. (2002).  In search of underlying dimensions: The use (and abuse) of factor analysis in personality and social psychology bulletin.  Personality and Social Psychology Bulletin, 28(12), 1629–1646.
  * Sabnis et al. (2016) Sabnis, G., Pati, D., Engelhardt, B., and Pillai, N. (2016).  A divide and conquer strategy for high dimensional Bayesian factor models.  arXiv:1612.02875.
  * Sarkar et al. (2021) Sarkar, A., Pati, D., Mallick, B. K., and Carroll, R. J. (2021).  Bayesian copula density deconvolution for zero-inflated data in nutritional epidemiology.  Journal of the American Statistical Association, 116, 1075–1087.
  * Schiavon et al. (2022) Schiavon, L., Canale, A., and Dunson, D. B. (2022).  Generalized infinite factorization models.  Biometrika, 109, 817–835.
  * Tan et al. (2020) Tan, S., Li, D., and Zhu, X. (2020).  Cancer immunotherapy: Pros, cons and beyond.  Biomedicine & Pharmacotherapy, 124, 109821.
  * Wang et al. (2014) Wang, C., Gong, B., et al. (2014).  The concordance between RNA-seq and microarray data depends on chemical treatment and transcript abundance.  Nature Biotechnology, 32, 926–932.
  * Watanabe (2013) Watanabe, S. (2013).  A widely applicable Bayesian information criterion.  Journal of Machine Learning Research, 14, 867–897.
  * Yoshida et al. (2019) Yoshida, H., Lareau, C. A., et al. (2019).  The cis-regulatory atlas of the mouse immune system.  Cell, 176, 897–912.

Supplementary Materials for

Inferring Covariance Structure from
Multiple Data Sources via
Subspace Factor Analysis

Noirrit Kiran Chandra†, David B. Dunson‡, Jason Xu‡,∗

†Department of Mathematical Sciences,

The University of Texas at Dallas, Richardson, TX

Email<EMAIL_ADDRESS>‡Department of Statistical Science, Duke
University, Durham, NC ∗Department of Biostatistics and Bioinformatics, Duke
University, Durham, NC

Supplementary materials provide further details on Hamiltonian Monte Carlo,
complete prior specifications, proofs and details of the theoretical results
and technical lemmas, extended simulation results, and details on the gene
expression datasets.

##  Appendix S.1 Details on Prior Specifications

##### Specifics of the Dirichlet-Laplace (DL) prior:

On a $d$-dimensional vector $\theta$, the DL prior with parameter $a$, denoted
by $\mathrm{DL}(a)$, can be specified hierarchically as

$\textstyle{\theta_{j}\mid\mbox{\boldmath$\phi$},\tau\stackrel{{\scriptstyle\mathrm{\scriptsize{ind}}}}{{\sim}}\mathrm{DE}(\phi_{j}\tau),~{}\mbox{\boldmath$\phi$}\sim\hbox{Dir}(a,\ldots,a),~{}\tau\sim\hbox{Ga}\left(da,\frac{1}{2}\right)},$
(S.1)

where $\theta_{j}$ is the $j^{th}$ element of $\theta$, $\tau\in\mathbb{R}$,
$\mbox{\boldmath$\phi$}\in\mathbb{R}^{d}$, $\mathrm{DE}(b)$ is a Laplace
distribution with mean 0 and variance $2b^{2}$,
$\hbox{Dir}(a_{1},\dots,a_{d})$ is the $d$-dimensional Dirichlet distribution,
and $\hbox{Ga}(a,b)$ is the gamma distribution with mean $a/b$ and variance
$a/b^{2}$. For straightforward Gibbs sampling, often the following equivalent
representation of model (S.1) is used

$\textstyle{\theta_{j}\mid\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$},\tau\stackrel{{\scriptstyle\mathrm{\scriptsize{ind}}}}{{\sim}}\mathrm{N}(0,\psi_{j}\phi_{j}^{2}\tau^{2}),~{}\psi_{j}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\hbox{Exp}\left(\frac{1}{2}\right),~{}\mbox{\boldmath$\phi$}\sim\hbox{Dir}(a,\ldots,a),~{}\tau\sim\hbox{Ga}\left(da,\frac{1}{2}\right)},$
(S.2)

where $\mbox{\boldmath$\psi$}\in\mathbb{R}^{d}$ and $\hbox{Exp}(a)$ is an
exponential distribution with mean $1/a$.

##### Choice of hyperparameters:

Following the suggestions by dir_laplace we set $a=\frac{1}{2}$ as the DL
hyperparameter. Recall that we have assumed
$a_{s,j,h}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(0,b_{{\mathbf{A}}})$
and
$\log\delta_{j}^{2}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(\mu_{\delta},\sigma_{\delta}^{2})$
where ${\mathbf{A}}_{s}=((a_{s,j,h}))$ and
$\bm{\Delta}=\mathrm{diag}(\delta_{1}^{2},\dots,\delta_{d}^{2})$. To specify
weakly informative priors, we set $b_{{\mathbf{A}}}=1$ and choose
$\mu_{\delta},\sigma_{\delta}^{2}$ such that $\mathbb{E}(\delta_{j}^{2})=1$
and $\mathrm{var}(\delta_{j}^{2})=7$ a priori for all $j=1,\dots,d$.

##  Appendix S.2 Proofs of Theoretical Results

##### Notations:

For two sequences $a_{n},b_{n}\geq 0$, $a_{n}\precsim b_{n}$ implies that
$a_{n}\leq Cb_{n}$ for some constant $C>0$; $a_{n}\asymp b_{n}$ implies that
$0<\liminf\absolutevalue{{a_{n}}/{b_{n}}}\leq\limsup\absolutevalue{{a_{n}}/{b_{n}}}<\infty$.
$\absolutevalue{{\mathbf{A}}}$ denotes the determinant of the square matrix
${\mathbf{A}}$. For a set $S$, $\absolutevalue{S}$ denotes its cardinality.
Let $\norm{{\mathbf{x}}}$ be the Euclidean norm of a vector ${\mathbf{x}}$. We
define
$\mbox{\boldmath$\Theta$}_{n}=\\{\bm{\Lambda}_{n},\mbox{\boldmath$\Delta$}_{n},{\mathbf{A}}_{1n},\dots,{\mathbf{A}}_{Sn}\\}$
as the set of all parameters and
$\mbox{\boldmath$\Theta$}_{0n}=\\{\bm{\Lambda}_{0n},\mbox{\boldmath$\Delta$}_{0n},{\mathbf{A}}_{01n},\dots,{\mathbf{A}}_{0S}\\}$
as the true data-generating values, and let $\mathbb{P}_{n}$ and
$\mathbb{P}_{0n}$ denote the joint distributions of
$\mathcal{D}_{n}=\\{{\mathbf{Y}}_{1,1},\dots,{\mathbf{Y}}_{1,n_{1}},\dots,{\mathbf{Y}}_{S,1},\dots,{\mathbf{Y}}_{S,n_{S}}\\}$
under $\mbox{\boldmath$\Theta$}_{n}$ and the true value
$\mbox{\boldmath$\Theta$}_{0n}$, respectively. We define
$\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)$ to be the
Kullback-Leibler (KL) divergence between $\mathbb{P}_{0n}$ and
$\mathbb{P}_{n}$.

For brevity of notations, we reuse $C$, $\widetilde{C}$, $C^{\prime\prime}$,
etc. in the proofs to denote constants whose values may not be the same
throughout the same proof. Nevertheless, we are careful to make sure that
these quantities are indeed constants.

We define the following quantities extensively used in the proofs

$\displaystyle\widetilde{\bm{\Lambda}}_{0n}$
$\displaystyle=\begin{bmatrix}\bm{\Lambda}_{0n}&{\mathbf{0}}_{d_{n}\times{(q_{n}-q_{0n})}}\end{bmatrix}:\bm{\Lambda}_{0n}\text{
padded with ${\mathbf{0}}$ columns to have the same order as
$\bm{\Lambda}_{n}$},$ $\displaystyle\widetilde{{\mathbf{A}}}_{0sn}$
$\displaystyle=\begin{bmatrix}{\mathbf{A}}_{0sn}&{\mathbf{0}}_{q_{0sn}\times{(q_{sn}-q_{0sn})}}\end{bmatrix}:{\mathbf{A}}_{0sn}\text{
padded with ${\mathbf{0}}$ columns to have the same order as
${\mathbf{A}}_{sn}$},$ $\displaystyle\widetilde{\bm{\Lambda}}_{0sn}$
$\displaystyle=\widetilde{\bm{\Lambda}}_{0n}\widetilde{{\mathbf{A}}}_{0sn}:\text{product
of the padded matrices to have the same order as
$\bm{\Lambda}_{n}{\mathbf{A}}_{sn}$}.$

######  Theorem 3 (KL support).

$\Pi_{n}\left\\{B_{n,0}(\mathbb{P}_{0n},\tau_{n})\right\\}\geq
e^{-Cn\tau_{n}^{2}}$ for $n\tau_{n}^{2}\asymp
c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})$.

###### Proof.

Since the priors on $\bm{\Lambda}_{n}$, ${\mathbf{A}}_{sn}$ and
$\mbox{\boldmath$\Delta$}_{n}$ are independent, from Lemma 9, we have

$\displaystyle\Pi_{n}\left\\{B_{n,0}(\mathbb{P}_{0n},\tau_{n})\right\\}\geq\Pi_{n}\left\\{\sum_{s=1}^{S}\frac{n_{s}\norm{\bm{\Sigma}_{0sn}-\bm{\Sigma}_{sn}}_{F}^{2}}{n{\delta}_{\min}^{2}s_{\min}\left(\bm{\Sigma}_{sn}\right)}\leq\tau_{n}^{2}\right\\}$
$\displaystyle\geq$
$\displaystyle\Pi_{n}\left(\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}<\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}\right)\times\Pi_{n}\left(\max_{1\leq
s\leq
S}\norm{{\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn}}_{F}<\frac{C\tau_{n}}{c_{n}\sqrt{q_{0n}}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\Pi_{n}\left\\{\norm{\mbox{\boldmath$\Delta$}_{0n}-\mbox{\boldmath$\Delta$}_{n}}_{F}\leq
C\tau_{n},s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)\geq\nu\right\\}.$
(S.3)

We handle the ${\mathbf{A}}_{s}$, $\bm{\Lambda}$ and $\bm{\Delta}$ parts in
(S.3) separately. We then conclude the proof by showing that each part
individually exceeds $e^{-Cn\tau_{n}^{2}}$ for some constant $C>0$.

The ${\mathbf{A}}_{sn}$ term in (S.3):

Note that the priors on ${\mathbf{A}}_{sn}$ are independent. Also from (C3)
$\norm{\widetilde{{\mathbf{A}}}_{0sn}}_{F}<\sqrt{q_{0sn}}\norm{{\mathbf{A}}_{0sn}}_{2}=o(\sqrt{q_{0sn}q_{0n}})$.
Further using Lemma 7 and the conditions $\sum_{s=1}^{S}q_{0sn}\leq q_{0n}$
and $\sum_{s=1}^{S}q_{sn}\leq q_{n}$, we obtain

$\displaystyle\Pi_{n}\left(\max_{1\leq s\leq
S}\norm{{\mathbf{A}}_{sn}-\widetilde{{\mathbf{A}}}_{0sn}}_{F}\leq\frac{C\tau_{n}}{c_{n}\sqrt{q_{0n}}}\right)$
$\displaystyle\geq
e^{-C^{\prime}\sum_{s=1}^{S}\max\\{q_{n}q_{sn}\log\frac{c_{n}\sqrt{q_{0n}}}{\tau_{n}},q_{n}q_{sn}\log(q_{n}q_{sn}),q_{0n}q_{0sn}\\}}$
$\displaystyle\geq
e^{-C^{\prime\prime}\max(q_{n}^{2}\log\frac{c_{n}\sqrt{q_{0n}}}{\tau_{n}},q_{n}^{2}\log
q_{n},q_{0n}^{2})}.$ (S.4)

Using (D1) we get $\max(q_{n}^{2}\log q_{n},q_{0n}^{2})=o(n\tau_{n}^{2})$, and
(D2) implies that
$q_{n}^{2}\log\frac{c_{n}\sqrt{q_{0n}}}{\tau_{n}}=o\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\\}$.
Thus the quantity in (S.4) exceeds $e^{-Cn\tau_{n}^{2}}$ for some $C>0$.

The $\bm{\Lambda}_{n}$ term in (S.3):

Using Pati et al. (2014, Lemma 7.1) we obtain

$\Pi_{n}\left(\norm{\bm{\Lambda}_{n}-\widetilde{\bm{\Lambda}}_{0n}}_{F}<\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}\right)\geq
e^{-C\max\left\\{\norm{\bm{\Lambda}_{0n}}_{F}^{2},s_{n}q_{0n}\log\frac{s_{n}q_{0n}}{\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}},\log(d_{n}q_{n})\right\\}}.$
(S.5)

Note that from (C3) $\norm{\bm{\Lambda}_{0n}}_{F}^{2}\leq
q_{0n}\norm{\bm{\Lambda}_{0n}}_{2}^{2}=o(c_{n})$ and hence
$\max\left\\{\norm{\bm{\Lambda}_{0n}}_{F}^{2},\log(d_{n}q_{n})\right\\}=o(n\tau_{n}^{2})$.
Additionally using (D2), we obtain that
$s_{n}q_{0n}\log\frac{s_{n}q_{0n}}{\frac{C\tau_{n}}{q_{0n}\sqrt{c_{n}}}}=\frac{1}{2}s_{n}q_{0n}\log\frac{ns_{n}q_{0n}^{3}}{c_{n}\log(d_{n}q_{n})}=o\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\\}$.
Hence the RHS of (S.5) exceeds $e^{-Cn\tau_{n}^{2}}$ for some $C>0$.

The $\mbox{\boldmath$\Delta$}_{n}$ term in (S.3):

Note that for a differentiable function $f(\cdot)$ in a compact interval
$[a,b]$, using the mean value theorem we have
$\absolutevalue{f(w)-f(y)}=f^{\prime}(x)(b-a)$ where $a<w<x<y<b$ implying that
$\absolutevalue{f(w)-f(y)}\leq(b-a)\times\sup_{x\in(a,b)}\absolutevalue{f^{\prime}(x)}$.
Thus we have $\nu\leq\min_{j}\delta_{jn}^{2}\leq\max_{j}\delta_{jn}^{2}\leq
c_{n}+M\Rightarrow\absolutevalue{\delta_{jn}^{2}-\delta_{0jn}^{2}}\leq(c_{n}+M)\log(c_{n}+M)\absolutevalue{\log\delta_{jn}^{2}-\log\delta_{0jn}^{2}}\leq
Cc_{n}\log c_{n}\absolutevalue{\log\delta_{jn}^{2}-\log\delta_{0jn}^{2}}$
where $C,M>0$ are positive large enough constants. Hence

$\displaystyle\Pi_{n}\left\\{\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq
C\tau_{n},s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)\geq\nu\right\\}$
$\displaystyle\geq$
$\displaystyle\Pi_{n}\left(\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq
C\tau_{n},\min_{1\leq j\leq d_{n}}\delta_{jn}^{2}\geq\nu,\max_{1\leq j\leq
d_{n}}\delta_{jn}^{2}\leq c_{n}+M\right)$ $\displaystyle\geq$
$\displaystyle\Pi_{n}\left(\norm{\widetilde{\boldsymbol{\delta}}_{n}-\widetilde{\boldsymbol{\delta}}_{0n}}\leq\frac{C^{\prime}\tau_{n}}{c_{n}\log
c_{n}},\min_{1\leq j\leq d_{n}}\delta_{jn}^{2}\geq\nu,\max_{1\leq j\leq
d_{n}}\delta_{jn}^{2}\leq c_{n}+M\right),$ (S.6)

where
$\widetilde{\boldsymbol{\delta}}_{n}=(\log\delta_{1n}^{2},\dots,\log\delta_{d_{n}n}^{2})^{\rm
T}$ and
$\widetilde{\boldsymbol{\delta}}_{0n}=(\log\delta_{01n}^{2},\dots,\log\delta_{0d_{n}n}^{2})^{\rm
T}$. Recall that
$\log\delta_{jn}^{2}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}(\mu_{\delta},\sigma_{\delta}^{2})$.
Adding a constant to the elements of $\widetilde{\boldsymbol{\delta}}_{0n}$
does not violate any assumptions and therefore without loss of generality the
parameter $\mu_{\delta}$ can be assumed to be zero.

Note that
$\norm{\widetilde{\boldsymbol{\delta}}_{n}-\widetilde{\boldsymbol{\delta}}_{0n}}\leq\frac{C^{\prime}\tau_{n}}{c_{n}\log
c_{n}}\Rightarrow\max_{j}\absolutevalue{\log\delta_{jn}^{2}-\log\delta_{0jn}^{2}}\leq\frac{C^{\prime}\tau_{n}}{c_{n}\log
c_{n}}$. Since $\tau_{n}=o(1)$ the last display in conjunction with (C4)
further implies that
$\max_{j}\delta_{jn}^{2}<\max_{j}\delta_{0jn}^{2}+C^{\prime\prime}\tau_{n}<c_{n}+M$
and
$\min_{j}\delta_{jn}^{2}>\min_{j}\delta_{0jn}^{2}-C^{\prime\prime}\tau_{n}>\nu$.
Hence the set
$\left\\{\widetilde{\boldsymbol{\delta}}_{n}:\norm{\widetilde{\boldsymbol{\delta}}_{n}-\widetilde{\boldsymbol{\delta}}_{0n}}\leq\frac{C^{\prime}\tau_{n}}{c_{n}\log
c_{n}}\right\\}$ is a subset of $\left\\{\min_{1\leq j\leq
d_{n}}\delta_{jn}^{2}\geq\nu,\max_{1\leq j\leq d_{n}}\delta_{jn}^{2}\leq
c_{n}+M\right\\}$. Therefore using Lemma 7

$\displaystyle\Pi_{n}\left\\{\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq
C\tau_{n},s_{\min}\left(\mbox{\boldmath$\Delta$}_{n}\right)\geq\nu\right\\}\geq\Pi_{n}\left(\norm{\widetilde{\boldsymbol{\delta}}_{n}-\widetilde{\boldsymbol{\delta}}_{0n}}\leq\frac{C^{\prime}\tau_{n}}{c_{n}\log
c_{n}}\right)$ $\displaystyle\geq$
$\displaystyle\exp{-C\max\left(d_{n}\log\frac{\sqrt{\sigma_{\delta}}c_{n}\log
c_{n}}{\tau_{n}},d_{n}\log
d_{n},\frac{1}{\sigma_{\delta}^{2}}\norm{\widetilde{\boldsymbol{\delta}}_{0n}}^{2}\right)}.$
(S.7)

Using (C3) we have
$\norm{\widetilde{\boldsymbol{\delta}}_{0n}}^{2}=o\\{d_{n}(\log
c_{n})^{2}\\}$. Hence $\max\left(d_{n}\log
d_{n},\frac{1}{\sigma_{\delta}^{2}}\norm{\widetilde{\boldsymbol{\delta}}_{0n}}^{2}\right)=o\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\\}$
using (C5). Additionally from (D2) $d_{n}\frac{n\sigma_{\delta}(\log
c_{n})^{2}}{s_{n}q_{0n}\log(d_{n}q_{n})}=o\\{c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})\\}$.
Thus the RHS of (S.7) exceeds $e^{-Cn\tau_{n}^{2}}$ for some $C>0$.

As we have shown that each of the product terms in (S.3) individually exceeds
$e^{-Cn\tau_{n}^{2}}$, we conclude the proof. ∎

######  Theorem 4 (Test function).

Recalling the definitions from the proof of Theorem 2, we let the set
$G_{j,n}=\\{\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}\leq
e_{n}(\mbox{\boldmath$\Theta$}_{n},\mbox{\boldmath$\Theta$}_{0n})\leq
2j\tau_{n}\\}$ denote an annulus of inner and outer radii
$n\tau_{n}^{2}c_{n}^{3}(c_{n}^{2}\epsilon_{n}+j\tau_{n})$ and
$n\tau_{n}^{2}c_{n}^{3}(c_{n}^{2}\epsilon_{n}+2j\tau_{n})$, respectively, in
spectral norm around $\bm{\Sigma}_{0n}$ for positive integer
$j=o\left[\frac{\sqrt{n}}{c_{n}^{5}}\left\\{s_{n}q_{0n}\log(d_{n}q_{n})\right\\}^{-\frac{3}{2}}\right]$.
Based on observed data $\mathcal{D}_{n}$ consider the following hypothesis
testing problem
$H_{0}:\mbox{\boldmath$\Theta$}_{n}=\mbox{\boldmath$\Theta$}_{0n}$ versus
$H_{1}:\mbox{\boldmath$\Theta$}_{n}\in G_{j,n}$.

Define
$\bm{\Lambda}_{0sn}=\bm{\Lambda}_{0n}\begin{bmatrix}{\mathbf{I}}_{q_{0n}}&{\mathbf{A}}_{0sn}\end{bmatrix}$,
$\widetilde{q}_{0,sn}=q_{0n}+q_{0sn}$,
$\bm{\Delta}_{\eta,0sn}=\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}$,
$\bm{\Sigma}_{\eta,0sn}=\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}\bm{\Delta}_{\eta,0sn}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}$
and
$\mbox{\boldmath$\Upsilon$}_{s,in}=({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn})^{-1}\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}{\mathbf{Y}}_{s,i}$. Further define
$\mbox{\boldmath$\Xi$}_{sn}=\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\left(\sum_{i=1}^{n_{s}}\mbox{\boldmath$\Upsilon$}_{s,in}\mbox{\boldmath$\Upsilon$}_{s,in}^{\rm
T}\right)\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}$ and the indicator functions
$\varphi_{sn}=\mathbbm{1}\left(\norm{\frac{1}{n_{s}}\mbox{\boldmath$\Xi$}_{sn}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}>\epsilon_{n}\right)$
for $s=1,\dots,S$. We define the test function:
$\varphi_{n}=1-\prod_{s=1}^{S}(1-\varphi_{sn})$. Then for some absolute
constant $K>0$

$\mathbb{E}_{H_{0}}(\varphi_{n})\to
0;\qquad\sup_{\mbox{\boldmath$\Theta$}_{n}\in
G_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{n})\leq\exp(-Knj^{2}\tau_{n}^{2}).$

###### Proof.

Type-I error:

Note that under $H_{0}$

$\displaystyle\mathrm{var}_{H_{0}}(\mbox{\boldmath$\Upsilon$}_{s,in})$
$\displaystyle=$
$\displaystyle\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}$
(S.8) $\displaystyle=$
$\displaystyle\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}(\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}+\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn})\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}$
$\displaystyle=$
$\displaystyle\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}(\bm{\Delta}_{\eta,0sn}^{2}+\bm{\Delta}_{\eta,0sn})\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}$
$\displaystyle=$
$\displaystyle\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}\bm{\Delta}_{\eta,0sn}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}=\bm{\Sigma}_{\eta,0sn},$

implying that
$\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\mbox{\boldmath$\Upsilon$}_{s,in}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{\widetilde{q}_{0,sn}})\Rightarrow\frac{1}{n_{s}}\mathbb{E}_{H_{0}}(\mbox{\boldmath$\Xi$}_{sn})={\mathbf{I}}_{\widetilde{q}_{0,sn}}$
for all $s=1,\dots,S$. From [Corollary 5.50]vershynin_book,
$\mathbb{E}_{H_{0}}(\varphi_{sn})\leq
2\exp(-\widetilde{C}\widetilde{q}_{0,sn}t_{n}^{2})$ for a universal constant
$\widetilde{C}>0$ and any sequence $t_{n}$ satisfying
$\widetilde{q}_{0,sn}t_{n}^{2}\leq n_{s}\epsilon_{n}^{2}$. From (C1), (C2) and
(D1) $\widetilde{q}_{0,sn}\to\infty$ and
$n_{s}\epsilon_{n}^{2}/\widetilde{q}_{0,sn}>c_{n}^{2}q_{0n}s_{n}\log(d_{n}q_{n})\to\infty$
as $n\to\infty$. Hence $t_{n}$ can be constructed such that $t_{n}\succ 1$ and
hence $\lim_{n\to\infty}\mathbb{E}_{H_{0}}(\varphi_{sn})=0$. Since the
$\varphi_{sn}$s are independent across $s$,
$\lim_{n\to\infty}\mathbb{E}_{H_{0}}(\varphi_{n})=1-\prod_{s=1}^{S}\lim_{n\to\infty}\\{1-\mathbb{E}_{H_{0}}(\varphi_{sn})\\}=0$.

Type-II error:

For data generating parameters $\mbox{\boldmath$\Theta$}_{n}\in G_{j,n}$, we
define
$\bm{\Sigma}_{\eta,sn}=\mathrm{cov}(\mbox{\boldmath$\Upsilon$}_{s,in})$. Then

$\bm{\Sigma}_{\eta,sn}=\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{sn}\bm{\Lambda}_{sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}.$

Hence for
${\mathbf{z}}_{s,i}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{\widetilde{q}_{0,sn}}({\mathbf{0}},{\mathbf{I}}_{\widetilde{q}_{0,sn}})$
for $i=1,\dots,n_{s}$,

$\displaystyle
1-\varphi_{sn}=\mathbbm{1}\left\\{\norm{\frac{1}{n_{s}}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\left(\sum_{i=1}^{n_{s}}\mbox{\boldmath$\Upsilon$}_{s,in}\mbox{\boldmath$\Upsilon$}_{s,in}^{\rm
T}\right)\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\leq\epsilon_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\frac{1}{n_{s}}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\left(\sum_{i=1}^{n_{s}}\mbox{\boldmath$\Upsilon$}_{s,in}\mbox{\boldmath$\Upsilon$}_{s,in}^{\rm
T}-\bm{\Sigma}_{\eta,sn}\right)\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}}_{2}\geq\norm{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}-\epsilon_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left(\norm{\bm{\Sigma}_{\eta,0sn}^{-1}\bm{\Sigma}_{\eta,sn}}_{2}\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\norm{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}-\epsilon_{n}\right).$
(S.9)

Note that,

$\displaystyle{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}={\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}(\bm{\Sigma}_{\eta,sn}-\bm{\Sigma}_{\eta,0sn})\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}}\text{
and}$
$\displaystyle\bm{\Sigma}_{\eta,sn}-\bm{\Sigma}_{\eta,0sn}=\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}{\bm{\Lambda}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-1}}(\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}){\mbox{\boldmath$\Delta$}_{0n}^{-1}}\bm{\Lambda}_{0sn}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-1}$
$\displaystyle\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}=\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}\left\\{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\right\\}^{-\frac{1}{2}}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}.$

We obtain the last identity in the above display from (S.8). Therefore,

$\displaystyle\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}=\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}\left\\{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\right\\}^{-\frac{1}{2}}$
$\displaystyle\times\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}{\bm{\Lambda}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-1}}(\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}){\mbox{\boldmath$\Delta$}_{0n}^{-1}}{\bm{\Lambda}_{0sn}}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}$
$\displaystyle\times\left\\{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\right\\}^{-\frac{1}{2}}\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}.$

Lemma 6 implies that

$\norm{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}\times
s_{\min}^{2}\left[\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}\times\right.\\\
\left.\left\\{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}^{\rm
T}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}\right\\}^{-\frac{1}{2}}\times\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}{\bm{\Lambda}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-1}}\right].$ (S.10)

Let us define
$\mbox{\boldmath\large$\chi$}_{0sn}=\mbox{\boldmath$\Delta$}_{0n}^{-\frac{1}{2}}\bm{\Lambda}_{0sn}$.
Then $\bm{\Delta}_{\eta,0sn}=\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}$ and the term inside the $s_{\min}^{2}$
in (S.10) simplifies to

$\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}\right)^{\frac{1}{2}}\times\left\\{\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}(\mbox{\boldmath\large$\chi$}_{0sn}\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}+{\mathbf{I}}_{d_{n}})\mbox{\boldmath\large$\chi$}_{0sn}\right\\}^{-\frac{1}{2}}\\\
\times\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}\right)^{-\frac{1}{2}}{\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-\frac{1}{2}}}.$ (S.11)

Letting ${\chi}_{0s,rn},r=1,\dots,\widetilde{q}_{0,sn}$ denote the eigenvalues
of $\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}$ and using Lemma 6, we have

$s_{\min}^{2}\left[\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}\right)^{\frac{1}{2}}\times\left\\{\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}(\mbox{\boldmath\large$\chi$}_{0sn}\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}+{\mathbf{I}}_{d_{n}})\mbox{\boldmath\large$\chi$}_{0sn}\right\\}^{-\frac{1}{2}}\right.\\\
\left.\times\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}\mbox{\boldmath\large$\chi$}_{0sn}\right)^{-\frac{1}{2}}{\mbox{\boldmath\large$\chi$}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-\frac{1}{2}}}\right]\geq
s_{\min}\left(\mbox{\boldmath$\Delta$}_{0n}^{-1}\right)\times\\\
\left[\min_{r=1,\dots,\widetilde{q}_{0,sn}}\left\\{(1+{\chi}_{0s,rn})^{\frac{1}{2}}\left({{\chi}_{0s,rn}}{(1+{\chi}_{0s,rn})}\right)^{-\frac{1}{2}}(1+{\chi}_{0s,rn})^{-\frac{1}{2}}{\chi}_{0s,rn}^{\frac{1}{2}}\right\\}\right]^{2}\\\
\geq\frac{1}{\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}}\times\min_{r=1,\dots,\widetilde{q}_{0,sn}}\frac{1}{1+{\chi}_{0s,rn}}\geq\frac{1}{\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}}\times\frac{1}{1+\norm{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}}_{2}}.$ (S.12)

From (C3)-(C4) $\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}=o(c_{n})$,
$\norm{\bm{\Lambda}_{0sn}}_{2}^{2}=o(c_{n})$. Now combining (S.10)-(S.12), we
get

$\displaystyle\norm{\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}\bm{\Sigma}_{\eta,sn}\bm{\Sigma}_{\eta,0sn}^{-\frac{1}{2}}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}$
$\displaystyle\geq\frac{\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}}{\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}(1+\norm{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\bm{\Lambda}_{0sn}}_{2})}$
$\displaystyle\geq\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}}{\norm{\mbox{\boldmath$\Delta$}_{0n}}_{2}\left(1+\frac{\norm{\bm{\Lambda}_{0sn}}_{2}^{2}}{{\delta}_{\min}^{2}}\right)}\geq\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}}{c_{n}^{2}}.$
(S.13)

Similarly, we have

$\displaystyle\norm{\bm{\Sigma}_{\eta,0sn}^{-1}\bm{\Sigma}_{\eta,sn}}_{2}\leq\norm{\bm{\Sigma}_{sn}}_{2}\times\left\|\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{\frac{1}{2}}\times\right.$
$\displaystyle\quad\left.\left\\{\bm{\Lambda}_{0sn}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-\frac{1}{2}}(\bm{\Lambda}_{0sn}\bm{\Lambda}_{0sn}+\mbox{\boldmath$\Delta$}_{0n})\mbox{\boldmath$\Delta$}_{0n}^{-\frac{1}{2}}\bm{\Lambda}_{0sn}\right\\}^{-\frac{1}{2}}\times\left({\mathbf{I}}_{\widetilde{q}_{0,sn}}+\bm{\Delta}_{\eta,0sn}\right)^{-\frac{1}{2}}{\bm{\Lambda}_{0sn}^{\rm
T}}{\mbox{\boldmath$\Delta$}_{0n}^{-1}}\right\|_{2}^{2}$ $\displaystyle\leq$
$\displaystyle\norm{\bm{\Sigma}_{sn}}_{2}\times\left[\max_{1\leq
r\leq\widetilde{q}_{0,sn}}\left\\{(1+{\chi}_{0s,rn})^{\frac{1}{2}}\left({{\chi}_{0s,rn}}{(1+{\chi}_{0s,rn})}\right)^{-\frac{1}{2}}(1+{\chi}_{0s,rn})^{-\frac{1}{2}}{\chi}_{0s,rn}^{\frac{1}{2}}\right\\}\right]^{2}$
(S.14) $\displaystyle\leq$
$\displaystyle\norm{\bm{\Sigma}_{n}}_{2}+\norm{\bm{\Lambda}_{n}{\mathbf{A}}_{sn}{\mathbf{A}}_{sn}^{\rm
T}\bm{\Lambda}_{n}^{\rm
T}}_{2}\leq\norm{\bm{\Sigma}_{n}}_{2}\left(1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\right).$
(S.15)

Note that
$\norm{\bm{\Sigma}_{n}}_{2}\leq\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}+c_{n}$.
Combining (S.9), (S.13) and (S.15) we get

$\displaystyle 1-\varphi_{sn}$ $\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\bm{\Sigma}_{n}}_{2}\left(1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\right)\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}}{c_{n}^{2}}-\epsilon_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\frac{\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}-\epsilon_{n}c_{n}^{2}}{c_{n}^{2}\left(\norm{\bm{\Sigma}_{n}-\bm{\Sigma}_{0n}}_{2}+c_{n}\right)\left(1+\norm{{\mathbf{A}}_{sn}}_{2}^{2}\right)}\right\\},$
(S.16)

where
${\mathbf{z}}_{s,1:n_{s}}\stackrel{{\scriptstyle\mathrm{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{\widetilde{q}_{0,sn}}({\mathbf{0}},{\mathbf{I}}_{\widetilde{q}_{0,sn}})$.
Using Lemma 10, we see that the RHS of (S.16) is bounded below by
$\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{1}{\sqrt{n_{s}}}$
for some absolute constant $\widetilde{C}>0$ if
$\mbox{\boldmath$\Theta$}_{n}\in G_{j,n}$. Further applying [Equation
5.26]vershynin_book, for all $\mbox{\boldmath$\Theta$}_{n}\in G_{j,n}$ we get

$\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq\Pr\left\\{\norm{\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}{\mathbf{z}}_{s,i}{\mathbf{z}}_{s,i}^{\rm
T}-{\mathbf{I}}_{\widetilde{q}_{0,sn}}}_{2}\geq\left(\widetilde{C}\sqrt{\widetilde{q}_{0,sn}}+j\sqrt{n}\tau_{n}\right)\frac{1}{\sqrt{n_{s}}}\right\\}\leq
e^{-Kj^{2}n_{s}\tau_{n}^{2}}.$

Finally, $\sup_{\mbox{\boldmath$\Theta$}_{n}\in
G_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{n})=\sup_{\mbox{\boldmath$\Theta$}_{n}\in
G_{j,n}}\prod_{s=1}^{S}\mathbb{E}_{\mbox{\boldmath$\Theta$}_{n}}(1-\varphi_{sn})\leq
e^{-K\left(\sum_{s=1}^{S}n_{s}\right)j^{2}\tau_{n}^{2}}$.

Hence the proof. ∎

######  Theorem 5 (Remaining probability mass).

$\lim_{n\to\infty}\mathbb{E}_{\mathbb{P}_{0n}}\Pi_{n}\left(\mathcal{P}_{2,n}\mid\mathcal{D}_{n}\right)=0$.

###### Proof.

For densities ${\mathbf{p}}_{0n}$ and ${\mathbf{p}}_{n}$ corresponding to
$\mathbb{P}_{0n}$ and $\mathbb{P}_{n}$, respectively, we define the average KL
variation
$\mathbb{V}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)=\int\left\\{\log\frac{{\mathbf{p}}_{0n}}{{\mathbf{p}}_{n}}-\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\right\\}^{2}\mathrm{d}{\mathbf{p}}_{0n}$
and the set
$B_{n,2}(\mathbb{P}_{0n},\tau)=\\{\mathbb{P}_{n}:\mathbb{KL}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq
n\tau^{2},~{}\mathbb{V}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq
n\tau^{2}\\}$.

From [Theorem 3.1]banerjee2015bayesian we find that the KL variation between
$\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Sigma}_{0sn})$ and
$\mathrm{N}_{d_{n}}({\mathbf{0}},\bm{\Sigma}_{sn})$ is
$\frac{1}{2}\sum_{j=1}^{d_{n}}(1-\psi_{j})^{2}$ where $\psi_{j}$’s are
eigenvalues of
$\bm{\Sigma}_{sn}^{-\frac{1}{2}}\bm{\Sigma}_{0sn}\bm{\Sigma}_{sn}^{-\frac{1}{2}}$.
Hence
$\mathbb{V}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)=\sum_{s=1}^{S}n_{s}\norm{\bm{\Sigma}_{sn}^{-1}(\bm{\Sigma}_{0sn}-\bm{\Sigma}_{sn})^{2}\bm{\Sigma}_{sn}^{-1}}_{F}^{2}\leq\sum_{s=1}^{S}\frac{n_{s}\norm{\bm{\Sigma}_{0sn}-\bm{\Sigma}_{sn}}_{F}^{2}}{2s_{\min}^{2}(\bm{\Sigma}_{sn})}$.
Therefore, the same three conditions in Lemma 9 imply that
$2\mathbb{V}\left(\mathbb{P}_{0n}\parallel\mathbb{P}_{n}\right)\leq
n\tau_{n}^{2}$ and accordingly
$\Pi_{n}\left\\{B_{n,2}(\mathbb{P}_{0n},\tau_{n})\right\\}>e^{-\widetilde{C}n\tau_{n}^{2}}$
for some absolute constant $\widetilde{C}>0$.

Note that the elements of ${\mathbf{A}}_{sn}$ are iid
$\mathrm{N}(0,b_{{\mathbf{A}}})$ across all $s$. Thus using [Theorem
5.39]vershynin_book we get that $\Pi_{n}(\norm{{\mathbf{A}}_{sn}}_{2}\geq
C\sqrt{n\tau_{n}^{2}})\leq e^{-C^{\prime}n\tau_{n}^{2}}$ with $C^{\prime}$
depending only on the constants $b_{{\mathbf{A}}}$ and $C$. Therefore
$\Pi_{n}(\mathcal{P}_{2,n})\leq\sum_{s=1}^{S}\Pi_{n}(\norm{{\mathbf{A}}_{sn}}_{2}\geq
C\sqrt{n\tau_{n}^{2}})\leq e^{-C^{\prime\prime}n\tau_{n}^{2}}$ where
$C^{\prime\prime}$ is a constant that depends only on the constants
$b_{{\mathbf{A}}}$, $S$ and $C$. Hence,
$\frac{\Pi_{n}(\mathcal{P}_{2,n})}{\Pi_{n}\left\\{B_{n,2}(\mathbb{P}_{0n},\tau_{n})\right\\}}\leq
e^{-2(C^{\prime\prime}-\widetilde{C})n\tau_{n}^{2}}=o(e^{-2n\tau_{n}^{2}})$.
The last display holds if the constant $C$ in the definition of
$\mathcal{P}_{2,n}$ is large enough, however, it is independent of $n$. Using
[Theorem 8.20]ghosal_book and subsequent application of DCT we conclude the
proof. ∎

######  Theorem 6.

Under the same conditions of Theorem 2,
$\Pi_{n}(\mbox{\boldmath$\Theta$}_{n}\in\mathcal{P}_{1,n}:\norm{\bm{\Sigma}_{sn}-\bm{\Sigma}_{0sn}}_{2}>M\varepsilon_{n}\mid\mathcal{D}_{n})\to
0$ in $\mathbb{P}_{0n}$-probability.

###### Proof.

Recall $\epsilon_{n}\asymp c_{n}n\tau_{n}^{3}$ and
$n\tau_{n}^{2}=c_{n}^{2}s_{n}q_{0n}\log(d_{n}q_{n})$ defined in the proof of
Theorem 2. Using these we now state a variation of Lemma 2. stated earlier.

######  Lemma 4.

Let ${\widetilde{e}}^{(s)}_{n}$ be a semimetric on $\mathcal{P}_{n}$ such that