latexReferences

Bayesian Scalable Precision Factor
Analysis for Massive Sparse
Gaussian Graphical Models

Noirrit Kiran Chandraa<EMAIL_ADDRESS>
Peter Müllerb,c<EMAIL_ADDRESS>
Abhra Sarkarb<EMAIL_ADDRESS>
aDepartment of Mathematical Sciences,

The University of Texas at Dallas,

800 W. Campbell Rd, Richardson, TX 75080-3021, USA

bDepartment of Statistics and Data Sciences,

The University of Texas at Austin,

2317 Speedway D9800, Austin, TX 78712-1823, USA

cDepartment of Mathematics,

The University of Texas at Austin,

2515 Speedway, PMA 8.100, Austin, TX 78712-1823, USA

Abstract

We propose a novel approach to estimating the precision matrix of multivariate
Gaussian data that relies on decomposing them into a low-rank and a diagonal
component. Such decompositions are very popular for modeling large covariance
matrices as they admit a latent factor based representation that allows easy
inference. The same is however not true for precision matrices due to the lack
of computationally convenient representations which restricts inference to
low-to-moderate dimensional problems. We address this remarkable gap in the
literature by building on a latent variable representation for such
decomposition for precision matrices. The construction leads to an efficient
Gibbs sampler that scales very well to high-dimensional problems far beyond
the limits of the current state-of-the-art. The ability to efficiently explore
the full posterior space also allows the model uncertainty to be easily
assessed. The decomposition crucially additionally allows us to adapt sparsity
inducing priors to shrink the insignificant entries of the precision matrix
toward zero, making the approach adaptable to high-dimensional small-sample-
size sparse settings. Exact zeros in the matrix encoding the underlying
conditional independence graph are then determined via a novel posterior false
discovery rate control procedure. A near minimax optimal posterior
concentration rate for estimating precision matrices is attained by our method
under mild regularity assumptions. We evaluate the method’s empirical
performance through synthetic experiments and illustrate its practical utility
in data sets from two different application domains.

Key Words: Factor models, False discovery rate control, Gaussian graphical
models, Markov chain Monte Carlo, Posterior concentration, Precision matrix
estimation, Scalable computation, Shrinkage priors

Short/Running Title: Precision Factor Analysis

Corresponding Author: Abhra Sarkar<EMAIL_ADDRESS>
##  1 Introduction

For multivariate Gaussian distributed data
${\mathbf{y}}=(y_{1},\dots,y_{d})^{\rm
T}\sim\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Sigma$})$, all conditional
dependence information is contained in the inverse covariance matrix
$\mbox{\boldmath$\Sigma$}^{-1}=\mbox{\boldmath$\Omega$}=((\omega_{j,j^{\prime}}))$,
or the precision. Two ‘nodes’ $y_{j}$ and $y_{j^{\prime}}$ are conditionally
independent given the rest if and only if $\omega_{j,j^{\prime}}=0$. The
underlying conditional (in)dependence graph is then obtained by connecting the
pairs of nodes $\\{(j,j^{\prime}):\omega_{j,j^{\prime}}\neq 0\\}$ by
undirected ‘edges’. Estimating the precision matrix, including especially its
sparsity patterns, for such data is therefore an important statistical problem
(Lauritzen, 1996; Koller and Friedman, 2009).

In this article, we model $\Omega$ via a low-rank and diagonal (LRD)
decomposition. Bhattacharya et al. (2016) introduced a representation for
efficient sampling from Gaussian distributions with known LRD structured
precision matrices. In this article, we adapt the representation in a
different way to obtain a novel factor analytic framework for unknown
precision matrices modeled via LRD decompositions with applications to
graphical models. While such decomposition is already a widely used standard
tool for modeling high-dimensional $\Sigma$, to our knowledge, the
construction has not been utilized before for modeling $\Omega$. Our research
addresses this remarkable gap in the literature. Although any positive
definite matrix can always be factorized this way, the main challenge is to
introduce such constructions for $\Omega$ in a way that allows efficient and
scalable posterior inference. Going significant steps further, we also adapt
this approach to sparse high-dimensional settings. Noting that any sparse
$\Omega$ can always be represented with a sparse factorization, we impose
sparsity in $\Omega$ by using sparsity-inducing priors for the factorization.

Existing Methods for Covariance and Precision Matrix Estimation: The existing
literature on sparse covariance and precision matrix estimation is vast. In
sparse covariance matrix estimation problems, the entire matrix $\Sigma$ is
often directly penalized (Levina et al., 2008; Bien and Tibshirani, 2011). In
an alternative approach, $\Sigma$ is assumed to admit a LRD structure
$\mbox{\boldmath$\Sigma$}=\widetilde{\mbox{\boldmath$\Lambda$}}\widetilde{\mbox{\boldmath$\Lambda$}}^{\rm
T}+\widetilde{\mbox{\boldmath$\Delta$}}$ where
$\widetilde{\mbox{\boldmath$\Lambda$}}$ is a $d\times q$ order matrix and
$\widetilde{\mbox{\boldmath$\Delta$}}$ is a diagonal matrix with all positive
entries. In theory, all positive definite matrices admit such a representation
for some $0\leq q\leq d$, and in practice, $q\ll d$ often suffice to produce
good approximations. Also importantly, this model admits the latent variable
representation
${\mathbf{y}}=\widetilde{\mbox{\boldmath$\Lambda$}}\widetilde{{\mathbf{u}}}+\widetilde{{\mathbf{v}}}$
where
$\widetilde{{\mathbf{u}}}\sim\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{q})$
and
$\widetilde{{\mathbf{v}}}\sim\mathrm{N}_{d}({\mathbf{0}},\widetilde{\mbox{\boldmath$\Delta$}})$.
The formulation allows massive scalability in computation, making LRD based
methods popular in the high-dimensional covariance matrix estimation
literature (Fan et al., 2011, 2018; Daniele et al., 2019), especially in
Bayesian settings (Bhattacharya and Dunson, 2011; Pati et al., 2014; Zhu et
al., 2014; Kastner, 2019, and others). A sparse $\Sigma$, however, does not
usually produce a sparse $\Omega$ and the strategy of inverting the estimated
$\Sigma$ to obtain an estimate of $\Omega$ tends to exhibit poor empirical
performance (Pourahmadi, 2013).

As in covariance matrix estimation problems, penalized likelihood based
methods that directly penalize the number and/or absolute values of non-zero
entries in $\Omega$ have also been developed in the frequentist setting (Yuan
and Lin, 2007; Banerjee et al., 2008; Rothman et al., 2008; d’Aspremont et
al., 2008; Friedman et al., 2008; Witten et al., 2011; Mazumder and Hastie,
2012; Zhang and Zou, 2014). Alternatively, Meinshausen and Bühlmann (2006);
Peng et al. (2009) developed neighborhood selection methods that learn the
edges by regressing each variable on the rest with penalties on large
regression coefficients.

When the main focus is on estimating the underlying dependence graph, Bayesian
approaches instead rely on defining a hierarchical prior on $\Omega$ preceded
by a prior on the graph. Choices for the latter include uniform priors over
graph sizes (Armstrong et al., 2009), priors centered around some informed
location (Mitra et al., 2013), priors with edge inclusions following a
binomial distribution (Dobra et al., 2004; Carvalho and Scott, 2009), a hyper-
Markov distribution on decomposable graphs (Dawid and Lauritzen, 1993), its
generalizations to non-decomposable settings (Roverato, 2002; Khare et al.,
2018), etc. However, such hierarchical construction with a separate model
layer for the underlying graph structure makes posterior exploration quite
challenging. Markov chain Monte Carlo (MCMC) algorithms have been designed
specifically for such models (Dellaportas et al., 2003; Atay-Kayis and Massam,
2005; Carvalho et al., 2007; Dobra et al., 2011; Green and Thomas, 2013;
Lenkoski, 2013; Mohammadi and Wit, 2015, and others) but these strategies
still rely on expensive local exploration moves, often involving trans-
dimensional proposals in the graph space and/or approximations of intractable
normalizing constants, hence remaining computationally infeasible beyond only
a few tens of dimensions (Jones et al., 2005).

Bayesian methods that directly penalize $\Omega$, thereby avoiding to have to
specify a separate prior for the underlying graph, have started to get some
attention but the literature remains sparse. Yoshida and West (2010) proposed
a factor model with complex constraints enforcing identical sparsity patterns
in the covariance and precision matrices which may be restrictive in practice
while also having scalability issues. Banerjee and Ghosal (2015) studied spike
and slab type priors (Ishwaran and Rao, 2005) to shrink irrelevant off-
diagonals to zero. For point mass mixture priors, MCMC based model space
exploration can generally be daunting and may lead to slow mixing and
convergence even in simple mean and linear regression problems. Continuous
shrinkage priors (Polson and Scott, 2010) that allow fast and efficient
posterior exploration have also been adapted for precision matrices. Wang
(2012); Khondker et al. (2013) and Li et al. (2019a) designed block Gibbs
samplers that allow updating entire columns of $\Omega$ at once. Mohammadi et
al. (2021) proposed an approximated sampler that can scale up to a few
hundreds but the problem remains infeasible for modern applications with many
thousands of nodes. More recent developments along these lines (Gan et al.,
2019; Li et al., 2019b; Deshpande et al., 2019) have focused on fast
deterministic Expectation-Maximization (EM) algorithms instead, which scale
well to problems with a few hundred dimensional nodes but only estimate the
posterior mode (MAP) and not the full posterior. Ksheera Sagar et al. (2021)
studied both MCMC and MAP estimation for element-wise horseshoe like priors on
the precision matrix and studied their convergence properties. The approach,
however, suffers from similar scalability issues. Also, in many of these
approaches (Friedman et al., 2008; Peng et al., 2009), the estimated $\Omega$
is not guaranteed to be positive definite, requiring post-hoc analysis to fix
the estimate.

Owing to the lack of computationally tractable hierarchical structures, in
high-dimensions, posterior explorations in existing Bayesian precision matrix
and graph estimation methods thus still remain prohibitively expensive if not
entirely impossible. With the few exceptions, existing Bayesian approaches
also do not come with rigorous theoretical guarantees.

Our Proposed Precision Factor Model: In contrast to covariance matrices, there
are currently no flexible LRD decomposition methods that admit easily
tractable latent variable representations for precision matrices, posing
significant methodological and computational challenges, especially for MCMC
based Bayesian inference.

In this article, we derive a factor analytic framework building on a
representation for Gaussian distributions with an LRD decomposed precision
matrix
$\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$}$ for some $d\times q$ order matrix $\Lambda$ and
some diagonal matrix $\Delta$ (Bhattacharya et al., 2016). The representation
is immediately useful in facilitating efficient posterior simulation in
Bayesian inference of $\Omega$ and easily scales to problems with dimensions
$d$ far beyond the limits of the current state-of-the-art.

Since all positive definite matrices admit LRD representations ($q=0$ and
$q=d$ being the two extremes), the proposed approach imposes no restrictive
assumptions on the precision matrix or the underlying graph. Additionally, we
discuss a constructive approach to find an LRD representation of arbitrary
sparse $\Omega$. In simulation examples, we show that $q\ll d$ often suffices
to produce good approximations of $\Omega$ even when they do not exactly admit
an LRD for such $q$. Conversely, since the representation always produces a
positive definite matrix, unlike many existing procedures such as Meinshausen
and Bühlmann (2006); Friedman et al. (2008); Gan et al. (2019), we obtain
positive definite estimates simply by design.

As the off-diagonals elements of $\Omega$ are contributed entirely by
$\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm T}$, a sparse
$\Lambda$ is expected to produce a sparse $\Omega$ (see Figure 3 and Section
2.2.2). A suitably chosen penalty on $\Lambda$ therefore allows to flexibly
adapt its non-zero elements such that insignificant off-diagonals of $\Omega$
are shrunk towards zero. The strategy has been successfully employed in high-
dimensional sparse covariance matrix estimation literature in both Bayesian
(Pati et al., 2014) and frequentist paradigms (Daniele et al., 2019).

In this article, we adapt the Dirichlet-Laplace shrinkage priors (Bhattacharya
et al., 2015) on $\Lambda$ for their theoretical and computational
tractability. As an artifact of Bayesian methods with continuous shrinkage
priors, exact zeroes do not appear in the posterior samples. Thus, we address
the problem of non-zero off-diagonal/edge selection from the posterior MCMC
samples through a novel multiple hypothesis testing approach that allows the
posterior false edge discovery rate (FDR) (Chandra and Bhattacharya, 2019) to
be controlled at any desired level. This is another salient feature of our
proposed method that properly accommodates posterior uncertainty in reporting
a point estimate for the graph.

With our model and prior specifications, we establish near minimax-optimal
contraction rates of the posterior around the true $\Omega$ under mild
assumptions when the number of nodes increases exponentially with the sample
size. We evaluate the proposed method’s finite sample efficacy through
simulation experiments where it either outperforms or is competitive with
previously existing methods in moderately high-dimensional problems while also
scaling to dimensions far beyond the reach of many of those methods. We
illustrate our method’s practical utility in real data sets from two different
application domains.

Our Key Contributions: Overall, our main contributions to the literature
include (a) proposing a likelihood based approach built on a low-rank and
diagonal decomposition of the precision matrix, (b) build on a latent factor
representation for such decomposition that provides new interpretations for
such models, (c) designing a Gibbs sampling algorithm that exploits this
latent factor representation and scales very well to high-dimensional
problems, (d) adapting shrinkage priors for the low-rank component that
further makes the method applicable to high-dimensional small-sample-size
sparse settings, (e) developing a novel FDR control procedure for graph
selection from the posterior samples of the precision matrix, and (f)
establishing rigorous asymptotic properties of the posterior of the proposed
approach.

Outline of the Article: The rest of this article is organized as follows.
Section 2 details our matrix decomposition based model. Section 2.1 discusses
our novel factor analytic representation; Section 2.2 discusses how a low-rank
approach can be used to learn arbitrary sparse precision matrices; Section 2.3
discusses the priors; Section 2.4 describes the posterior sampling algorithm;
Section 2.5 presents our graph selection procedure via FDR control; Section
2.6 discusses the posterior’s asymptotic properties. Section 3 summarizes the
results of simulation experiments. Section 4 presents the results for two real
data sets. Section 5 contains concluding remarks.

##  2 Gaussian Precision Factor Models

We consider a random sample of $n$ observations assumed to be independently
and identically distributed (iid) $d$-dimensional random vectors following a
multivariate Normal distribution with mean zero and precision matrix $\Omega$
as

$\displaystyle{\mathbf{y}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Omega$}^{-1}),~{}~{}~{}i=1,\dots,n,$
(1)

where ${\mathbf{y}}_{i}=(y_{i,1},\dots,y_{i,d})^{\rm T}$. The primary goal is
to estimate $\Omega$, especially identifying its sparsity pattern that
characterizes conditional independence relationships between different
components of ${\mathbf{y}}$.

To this end, we consider an LRD decomposition of $\Omega$ as

$\displaystyle\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$},$ (2)

where $\mbox{\boldmath$\Lambda$}^{d\times q}=((\lambda_{r,c}))$ and
$\mbox{\boldmath$\Delta$}=\mathrm{diag}(\delta_{1}^{2},\dots,\delta_{d}^{2})$.
Factorization (2) is completely flexible in the sense that a matrix is
positive definite if and only if it admits such a representation for some
$0\leq q\leq d$. Any precision matrix $\Omega$ can therefore be written as (2)
for a sufficiently large value of $q$. For most practical cases, values of
$q\ll d$ suffice to approximate $\Omega$ well, resulting in a huge reduction
in dimensions and allowing massive scalability in computation. In model (2),
$\Lambda$ is not strictly identifiable since, for example,
$\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}=\mbox{\boldmath$\Lambda$}{\mathbf{Q}}{\mathbf{Q}}^{\rm
T}\mbox{\boldmath$\Lambda$}^{\rm T}$ for any orthogonal matrix ${\mathbf{Q}}$.
Inference on $\Omega$ being the primary interest, individual identifiability
or interpretation of the model parameters in (2) is, however, not required.

###  2.1 Factor Analytic Representation

One main advantage of modeling covariance matrices via LRD decompositions is
the existence of a latent factor representation that greatly facilitates
computation. Tractable latent factor representations are, however, not known
to exist for precision matrices, presenting major barriers against efficient
inference. Recently Bhattacharya et al. (2016) introduced a representation
that greatly facilitates fast sampling from multivariate Gaussian
distributions with known precision matrices with an LRD structure. We show
that, by exploiting this representation in a different way, a latent factor
model can in fact be obtained for LRD decomposed precision matrices as well
that allows efficient posterior inference of its unknown components in much
the same way as classical factor models facilitate computation for unknown LRD
decomposed covariance matrices. The representation, formalized in Proposition
1 below and referred to in this article as the ‘precision factor model’, also
provides novel insights into the very construction of Gaussian graphical
models.

######  Proposition 1.

The model
${\mathbf{y}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Omega$}^{-1}),i=1,\dots,n$,
with
$\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$}$, where $\Lambda$ is ${d\times q}$ with $q\leq d$
and
$\mbox{\boldmath$\Delta$}=\mathrm{diag}(\delta_{1}^{2},\dots,\delta_{d}^{2})$,
admits the equivalent representation

$\displaystyle{\mathbf{y}}_{i}=-\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}({\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$})^{-1}{\mathbf{u}}_{i}+{\mathbf{v}}_{i},$
(3) $\displaystyle\text{where}~{}\begin{bmatrix}{\mathbf{u}}_{i}\\\
{\mathbf{v}}_{i}\end{bmatrix}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q+d}\left({\mathbf{0}},\begin{bmatrix}{\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}&\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\\\
\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}&\mbox{\boldmath$\Delta$}^{-1}\end{bmatrix}\right).$
(4)

The proposition follows straightforwardly using Sherman-Woodbury identity for
the covariance matrix $\mbox{\boldmath$\Sigma$}=\mbox{\boldmath$\Omega$}^{-1}$
given by

$\displaystyle\mbox{\boldmath$\Omega$}^{-1}=(\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$})^{-1}=\mbox{\boldmath$\Delta$}^{-1}-\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}({\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$})^{-1}\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}.$

Note that $\hbox{cov}({\mathbf{y}}_{i},{\mathbf{u}}_{i})={\mathbf{0}}$ in (3).
A reverse-engineered construction that is particularly useful for posterior
simulation in Bayesian settings utilizes this fact and generates
${\mathbf{u}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{P}})$
with ${\mathbf{P}}=({\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$})$ independent of
${\mathbf{y}}_{1:n}$ first and then sets
${\mathbf{v}}_{i}=\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}+{\mathbf{y}}_{i}$.
Clearly then,
${\mathbf{v}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Delta$}^{-1})$.
Also,

$\displaystyle\begin{bmatrix}{\mathbf{u}}_{i}\\\
{\mathbf{y}}_{i}\end{bmatrix}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q+d}\left({\mathbf{0}},\begin{bmatrix}{\mathbf{P}}&{\mathbf{0}}\\\
{\mathbf{0}}&\mbox{\boldmath$\Omega$}^{-1}\end{bmatrix}\right)\text{ then
implies }\begin{bmatrix}{\mathbf{u}}_{i}\\\
{\mathbf{v}}_{i}\end{bmatrix}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q+d}\left({\mathbf{0}},\begin{bmatrix}{\mathbf{P}}&\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\\\
\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}&\mbox{\boldmath$\Delta$}^{-1}\end{bmatrix}\right),$

as in (4) in Proposition 1. Importantly, we can also write

$\displaystyle{\mathbf{u}}_{i}=\mbox{\boldmath$\Lambda$}^{\rm
T}{\mathbf{v}}_{i}+\mbox{\boldmath$\varepsilon$}_{i},\text{ where
}\mbox{\boldmath$\varepsilon$}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{q}).$
(5)

Although mathematically simple, Proposition 1 has far-reaching implications.
An efficient and highly scalable Gibbs sampler follows immediately from the
construction by first generating the latent vectors
${\mathbf{u}}_{i},{\mathbf{v}}_{i}$ given
$\mbox{\boldmath$\Lambda$},\mbox{\boldmath$\Delta$}$ and ${\mathbf{y}}_{i}$ as
described above, and then updating the rows of $\Lambda$ given
${\mathbf{u}}_{1:n}$ and ${\mathbf{v}}_{1:n}$ using (5).

Figure 1: Graphical view of the precision factor model: Plots of observed
${\mathbf{y}}_{i}$, latent
${\mathbf{w}}_{i}=\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}$
and latent ${\mathbf{v}}_{i}={\mathbf{w}}_{i}+{\mathbf{y}}_{i}$ for a model
with dimension $d=3$ and rank $q=2$ from different viewing angles. The
$d$-dimensional vectors ${\mathbf{w}}$ are supported on a lower
$q$-dimensional plane. Also, ${\mathbf{u}}$ and ${\mathbf{y}}$ are independent
of each other. Since ${\mathbf{w}}$ are linear transformations of
${\mathbf{u}}$, they are also independent of ${\mathbf{y}}$. In the plots,
${\mathbf{w}}$ and the average values of ${\mathbf{y}}$ can be seen to be
living on planes that are orthogonal to each other. The vectors ${\mathbf{v}}$
are more scattered than ${\mathbf{y}}$ and, in obtaining the variance-
covariance of ${\mathbf{y}}={\mathbf{v}}-{\mathbf{w}}$, the larger variances
of ${\mathbf{v}}$ are perfectly compensated by the negative covariances
between ${\mathbf{v}}$ and ${\mathbf{w}}$.

Continuing the thread at the beginning of this subsection, parallels can be
drawn with latent factor models for covariance matrices where a similar
decomposition
$\mbox{\boldmath$\Sigma$}=\widetilde{\mbox{\boldmath$\Lambda$}}\widetilde{\mbox{\boldmath$\Lambda$}}^{\rm
T}+\widetilde{\mbox{\boldmath$\Delta$}}$ arises from the model
${\mathbf{y}}_{i}=\widetilde{\mbox{\boldmath$\Lambda$}}\widetilde{\mathbf{u}}_{i}+\widetilde{\mathbf{v}}_{i}$
with independent latent components $\widetilde{\mathbf{u}}_{i}$ and
$\widetilde{\mathbf{v}}_{i}$, where the latent factors
$\widetilde{\mathbf{u}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{I}}_{q})$,
with $q$ typically $\ll d$, and the errors
$\widetilde{\mathbf{v}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\widetilde{}\mbox{\boldmath$\Delta$})$.
The covariance between the components of ${\mathbf{y}}_{i}$ and a part of the
variance of ${\mathbf{y}}_{i}$ are thus explained by (a) the variance-
covariance of the latent factors $\widetilde{\mathbf{u}}_{i}$ (b) while the
remaining unexplained variance is attributed to the errors
$\widetilde{\mathbf{v}}_{i}$.

In contrast, for the precision factor model depicted in Figure 1, we have from
model (3) that
${\mathbf{y}}_{i}=-\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}+{\mathbf{v}}_{i}$
with dependent latent components ${\mathbf{u}}_{i}$ and ${\mathbf{v}}_{i}$,
where the latent factors
${\mathbf{u}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{P}})$,
with $q$ expected again to be $\ll d$, and the ‘errors’
${\mathbf{v}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Delta$}^{-1})$.
The variance-covariance of ${\mathbf{y}}_{i}$ is thus explained by (a) the
variance-covariance of the latent factors ${\mathbf{u}}_{i}$, (b) the variance
of the ‘errors’ ${\mathbf{v}}_{i}$, and (c) the covariance between
${\mathbf{u}}_{i}$ and ${\mathbf{v}}_{i}$.

If independence between the variance contributing components is desired, an
alternative view
${\mathbf{v}}_{i}=\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}+{\mathbf{y}}_{i}$
of model (3) is to see the latent ${\mathbf{v}}_{i}$’s be composed of
independent components ${\mathbf{u}}_{i}$ and ${\mathbf{y}}_{i}$, where the
latent factors
${\mathbf{u}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{P}})$,
with $q\ll d$ as before, and the ‘errors’
${\mathbf{y}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Sigma$})$.
The ${\mathbf{v}}_{i}$’s can thus be represented in an orthogonal
decomposition with components
${\mathbf{w}}_{i}=\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}$
and ${\mathbf{y}}_{i}$. The larger variances $\mbox{\boldmath$\Delta$}^{-1}$
of ${\mathbf{v}}_{i}$ are now being explained by (a) the variance-covariance
of the latent factors ${\mathbf{u}}_{i}$ and (b) the variance-covariance of
the ‘errors’ ${\mathbf{y}}_{i}$, the off-diagonal covariance terms of these
components perfectly cancelling each other to produce the diagonal matrix
$\mbox{\boldmath$\Delta$}^{-1}$.

In yet another view based on equation (5), the vectors ${\mathbf{v}}_{i}$ can
instead be interpreted as independent heterogeneous latent factors and
${\mathbf{u}}_{i}$ the associated response vectors subject to white noises
$\mbox{\boldmath$\varepsilon$}_{i}$. Contrary to the classical factor model,
in this view, we have the factor dimension $d\gg$ the response dimension $q$.
Clearly, $\mbox{\boldmath$\Omega$}^{-1}$ is now the conditional variance-
covariance of ${\mathbf{v}}_{i}$ given ${\mathbf{u}}_{i}$, that is, given the
response variables ${\mathbf{u}}_{i}$, we can study the behavior of the latent
factors ${\mathbf{v}}_{i}$ by examining $\Omega$.

Although the precision factor model described here has immediate practical
implications for Bayesian inference of Gaussian precision matrices considered
in this article, the representation itself is not specific to the chosen
inferential paradigm but is a mathematical identity that may be of broad
general interest in understanding the basic construction of such models.

###  2.2 Low-rank Modeling of Sparse Precision Matrices

In this section we discuss how a sparse precision matrix can indeed admit an
LRD decomposition. Then we discuss our strategy of inducing sparsity in high-
dimensional precision matrix estimation problems.

####  2.2.1 LRD Decomposition of Sparse Precision Matrices

We first discuss how a sparse $\Omega$ can admit an LRD decomposition. To get
some insights, consider the stylized example from Figure 2 where $d=5$ and
$\Omega$ has $\widetilde{s}=3$ non-zero off-diagonals. We recall the factor
analytic representation of our model from Section 2.1 where we write
${\mathbf{u}}=\mbox{\boldmath$\Lambda$}^{\rm
T}{\mathbf{v}}+\mbox{\boldmath$\varepsilon$}$ in equation (5) and interpret
$\mbox{\boldmath$\Omega$}=((\omega_{j,h}))$ as the precision matrix of
${\mathbf{v}}$ conditionally on ${\mathbf{u}}$. Under this representation,
$\omega_{j,h}\neq 0$ if and only if $\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{h}\neq 0$, that is, if $v_{j}$ and $v_{h}$ are
connected to the same $u_{\ell}$, $1\leq\ell\leq q$.

Figure 2: Constructing an LRD decomposition of a sparse
$\mbox{\boldmath$\Omega$}^{5\times 5}$: For each non-zero $\omega_{j,h}$, we
connect $v_{j}$ and $v_{h}$ to $u_{\ell}$. Notably, $y_{1}$ is marginally
independent to every other variable and therefore we do not connect $v_{1}$ to
any $u_{\ell}$.

Likewise, we can construct a sparse $\Lambda$ as follows

$\underbrace{\begin{bmatrix}u_{1}\\\ u_{2}\\\
u_{3}\end{bmatrix}}_{{\mathbf{u}}}=\underbrace{\begin{bmatrix}0&\lambda_{2,1}&0&0&\lambda_{5,1}\\\
0&0&\lambda_{3,2}&\lambda_{4,2}&0\\\
0&0&0&\lambda_{4,3}&\lambda_{5,3}\end{bmatrix}}_{\mbox{\boldmath$\Lambda$}^{\rm
T}}\underbrace{\begin{bmatrix}v_{1}\\\ v_{2}\\\ v_{3}\\\ v_{4}\\\ v_{5}\\\
\end{bmatrix}}_{{\mathbf{v}}}+\underbrace{\begin{bmatrix}\varepsilon_{1}\\\
\varepsilon_{2}\\\
\varepsilon_{3}\end{bmatrix}}_{\mbox{\boldmath$\varepsilon$}},$

subject to $\lambda_{2,1}\lambda_{5,1}=\omega_{2,5}$,
$\lambda_{3,2}\lambda_{4,2}=\omega_{3,4}$,
$\lambda_{4,3}\lambda_{5,3}=\omega_{4,5}$ and
$\delta_{j}^{2}+\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{j}=\omega_{j,j}$ for all $j=1,\dots,d$. Although
this construction is not unique, we see that a sparse LRD decomposition indeed
exists for our concerned $\Omega$.

We can generalize this argument for arbitrary sparse $\Omega$ with
$\widetilde{s}$ number of non-zero off-diagonals. We start with a $d\times q$
order $\Lambda$ with all entries equal to zero. Let ${\cal E}$ be the set of
edges in the conditional dependence graph. Note that each edge corresponds to
a partial correlation between $Y_{j}$ and $Y_{h}$, $1\leq j,h\leq d$. For an
edge $\omega_{j,h}\in{\cal E}$, we set $\lambda_{j,\ell_{1}}$ and
$\lambda_{h,\ell_{1}}$ to be non-zero for some $\ell_{1}$ so that
$\mbox{\boldmath$\lambda$}_{j}^{\rm T}\mbox{\boldmath$\lambda$}_{h}\neq 0$.
For the next edge $\omega_{j^{\prime},h^{\prime}}\in{\cal E}$, suppose
$\omega_{j,j^{\prime}}=\omega_{j,h^{\prime}}=\omega_{h,h^{\prime}}=0$. Then we
set $\lambda_{j,\ell_{2}}$ and $\lambda_{h,\ell_{2}}$ to be non-zero for some
$\ell_{2}\neq\ell_{1}$ so that $\mbox{\boldmath$\lambda$}_{j^{\prime}}^{\rm
T}\mbox{\boldmath$\lambda$}_{h^{\prime}}\neq 0$ but
$\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{j^{\prime}}=\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{h^{\prime}}=\mbox{\boldmath$\lambda$}_{h}^{\rm
T}\mbox{\boldmath$\lambda$}_{h^{\prime}}=0$, and so on. Note that in this
construction (which need not be unique), for each edge
$\omega_{j,j^{\prime}}\neq 0$, we need to connect $v_{j}$ and $v_{j^{\prime}}$
to a $u_{\ell}$. Therefore, we do not need $q$ to be larger than
$\widetilde{s}$. The existence of such a $\Lambda$ and $\Delta$ is subject to
the solution of the following system of quadratic equations on
$\\{\mathrm{vec}(\mbox{\boldmath$\Lambda$}),\mbox{\boldmath$\delta$}_{1:d}\\}$
with $dq+d$ unknowns and $d+\widetilde{s}$ equations

$\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{h}=\omega_{j,h}\text{ for all
}\omega_{j,h}\in{\cal
E},\quad\delta_{j}^{2}+\mbox{\boldmath$\lambda$}_{j}^{\rm
T}\mbox{\boldmath$\lambda$}_{j}=\omega_{j,j}\text{ for all }j=1,\dots,d.$

In practical high-dimensional applications, even though $d$ is large, the
conditional dependence graph is often sparse, that is, $\widetilde{s}\ll d$.
Barvinok and Rudelson (2022) noted that when the number of unknown variables
are considerably larger than the number of equations, a (not necessarily
unique) solution to quadratic system of equations exists under very generic
conditions.

####  2.2.2 Inducing Sparsity in $\Omega$ via Sparsity in $\Lambda$

$\times$

$=$

Figure 3: Sparse $\Lambda$ producing a sparse
$\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm T}$. Here, the gray
cells represent exact zeros, and purple to yellow represent smaller (negative)
to larger (positive) values.

The off-diagonals elements of
$\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$}$ are contributed entirely by
$\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm T}$. This allows to
achieve sparsity in $\Omega$ by inducing sparsity in $\Lambda$. Following the
discussion in Section 2.2.1, a carefully designed data-adaptive penalty on
$\Lambda$ (e.g., a shrinkage prior on $\Lambda$ that sufficiently increases
the probability of obtaining a sparse $\Omega$) would induce sparsity in
$\Lambda$ in such a (not necessarily unique) way that the sparsity patterns in
$\Omega$ are also accurately recovered. In this section, we show that by
inducing sparsity in $\mbox{\boldmath$\Lambda$}^{d\times
q}=((\lambda_{j,h}))$, sparsity is also induced in
$\mbox{\boldmath$\Omega$}^{d\times d}=((\omega_{j,j^{\prime}}))$. We begin
with the spike-and-slab prior (Ishwaran and Rao, 2005) that allows exact zeros
via a spike at zero and large nonzero elements via a continuous slab
$g(\cdot)$ as

$\lambda_{j,h}\mid\pi,\theta_{j,h}\sim\pi\mathbbm{1}_{\\{0\\}}(\cdot)+(1-\pi)g(\cdot),\quad\pi\sim\hbox{Beta}(a_{\pi},b_{\pi}),$

where $\pi$ is the prior probability of observing a zero for $\lambda_{j,h}$.
For such priors $\Pi(\omega_{j,j^{\prime}}=0\mid\pi)$ can be analytically
reduced to a simple form that helps gain insights into how inducing sparsity
in $\Lambda$ can induce sparsity in $\Omega$ with high probability.
Specifically, since
$\omega_{j,j^{\prime}}=\mbox{\boldmath$\lambda$}_{j}\mbox{\boldmath$\lambda$}_{j^{\prime}}^{\rm
T}$, where $\mbox{\boldmath$\lambda$}_{j}$ is the $j^{th}$ row of $\Lambda$,
$\omega_{j,j^{\prime}}=0$ can only happen when the $r^{th}$ entries of
$\mbox{\boldmath$\lambda$}_{j}$ and $\mbox{\boldmath$\lambda$}_{j^{\prime}}$
are both not from the slab distribution $g()$ for all $r=1,\dots,q$.
Therefore, we have

$\textstyle{\Pi(\omega_{j,j^{\prime}}=0\lvert\pi)=\sum_{r=0}^{q}\left\\{{q\choose
r}\pi^{r}(1-\pi)^{q-r}\right\\}\times\pi^{q-r}=\pi^{q}(2-\pi)^{q}}.$

The behavior of $\Pi(\omega_{j,j^{\prime}}=0\mid\pi)$ for varying values of
$\pi$ and $q$ are shown in Figure 4. It can be seen that by controlling $\pi$,
it is possible to induce any desired level of sparsity in $\Omega$. The
hierarchical Beta prior on $\pi$ allows for data-adaptive learning and shrink
the elements of $\Omega$ accordingly (Scott and Berger, 2010).

Figure 4: $\pi$ versus $\Pi(\omega_{j,j^{\prime}}=0\mid\pi)$ for different
values of $q$. Clearly, by controlling $\pi$, it is possible to induce desired
level of sparsity in $\Omega$. A hierarchical prior on $\pi$ adaptively learns
from the data and shrinks the elements of $\Omega$ accordingly.

Implementation of spike-and-slab type mixture priors is computationally
challenging. Thus, continuous global-local shrinkage priors have gained
popularity in the Bayesian sparse estimation literature as they often greatly
simplify posterior computation (Polson and Scott, 2010) while retaining almost
similar statistical properties of the classical spike-and-slab prior. In
particular, we use the Dirichlet-Laplace (DL, Bhattacharya et al., 2015) prior
that has equivalent asymptotic properties with respect to the spike-and-slab
in a similar context (Pati et al., 2014, Theorem 5.1) while being amenable to
scalable posterior computation. The results for the spike-and-slab discussed
here provides the intuitions on how continuous shrinkage priors like the DL
induce sparsity in $\Omega$ by inducing shrinkage on $\Lambda$.

####  2.2.3 Advantages of the LRD Decomposition

For high-dimensional sparse covariance matrix estimation, this LRD
decomposition strategy is extensively used in both Bayesian (Pati et al.,
2014) and frequentist paradigms (Daniele et al., 2019) where the associated
latent factor formulation allows scalable computation in big data problems
(Bhattacharya and Dunson, 2011; Fan et al., 2011; Sabnis et al., 2016; Fan et
al., 2018; Kastner, 2019, and others). Penalizing $\Lambda$ is thus a sensible
approach to induce sparsity in $\Omega$. It also comes with many practical
advantages described below.

Cholesky factorization based covariance and precision matrix estimation
methods also use a similar strategy and penalize the lower triangular matrix
${\mathbf{L}}$ to induce sparsity in
$\mbox{\boldmath$\Omega$}={\mathbf{L}}{\mathbf{L}}^{\rm T}$ (Dallakyan and
Pourahmadi, 2020). However, for undirected graphs, the estimates can vastly
differ depending on the ordering of the variables (Kang and Deng, 2020). In
contrast, our proposed LRD decomposition based approach with a penalty on
$\Lambda$ is invariant to the ordering of the variables.

In most existing approaches that directly penalize $\Omega$, frequentist or
Bayesian, ensuring the positive definiteness of $\Omega$ is also non-trivial,
particularly in high dimensional settings. In frequentist regimes, for
example, the graphical lasso (Friedman et al., 2008) does not guarantee
positive definiteness of $\Omega$ (Mazumder and Hastie, 2012); similar
behavior has also been observed for neighborhood selection methods
(Meinshausen and Bühlmann, 2006; Peng et al., 2009) and post-hoc treatments
are required for positive definiteness. In Bayesian MAP approaches such as Gan
et al. (2019), special care is needed in designing the optimization algorithm
to ensure positive definiteness. For many existing Bayesian methods relying on
MCMC (Wang, 2012; Khondker et al., 2013; Li et al., 2019a), constrained update
of $\Omega$ in each step is difficult and expensive. Applications of these
methods thus remain fairly limited to small to moderate dimensional problems.
In contrast, the LRD formulation is always guaranteed to produce a positive
definite estimate by construction.

As also discussed in the Introduction, posterior exploration is extremely
challenging in Bayesian approaches that penalize the underlying the dense
graphs first and then assign a prior on $\Omega$ conditional on the graph,
often requiring restrictive assumptions on the graph (Dawid and Lauritzen,
1993) while still remaining computationally infeasible beyond only a few tens
of dimensions (Jones et al., 2005). Our proposed approach on the other hand,
works directly in the $\Omega$ space, albeit via the LRD decomposition,
avoiding having to define separate complex priors on the graph space thereby
also avoiding restrictive assumptions on the graph while also achieving
scalability far beyond such approaches via its latent factor representation.

###  2.3 Prior Specification

To accommodate high-dimensional data, with $d\gg n$, it is crucial to reduce
the effective number of parameters in the $d\times q$ loadings matrix
$\Lambda$. A wide variety of sparsity inducing shrinkage priors for $\Lambda$
can be considered for this purpose. In this article, we employ a two-parameter
generalization of the original Dirichlet-Laplace (DL) prior from Bhattacharya
et al. (2015) that allows more flexible tail behavior. On a $d$-dimensional
vector $\theta$, our DL prior with parameters $a$ and $b$, denoted by
$\mbox{DL}(a,b)$, can be specified in the following hierarchical manner

$\displaystyle\theta_{j}\lvert\mbox{\boldmath$\psi$},\mbox{\boldmath$\phi$},\tau\stackrel{{\scriptstyle\mbox{\scriptsize{ind}}}}{{\sim}}\hbox{N}(0,\psi_{j}\phi_{j}^{2}\tau^{2}),~{}~{}\psi_{j}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\hbox{Exp}(1/2),~{}~{}\mbox{\boldmath$\phi$}\sim\hbox{Dir}(a,\ldots,a),~{}~{}\tau\sim\hbox{Ga}(da,b),$

where $\theta_{j}$ is the $j^{th}$ element of $\theta$, $\phi$ is a vector of
same length as $\theta$, $\hbox{Exp}(a)$ is an exponential distribution with
mean $1/a$, $\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}$. The original DL prior is a special case with $b=1/2$.
We let $\mbox{vec}(\mbox{\boldmath$\Lambda$})\sim\mbox{DL}(a,b)$.

The column dimension $q$ of $\Lambda$ will almost always be unknown. Assigning
a prior on $q$ and implementing a reversible jump MCMC (Green, 1995) type
algorithm can be inefficient and expensive. In this paper, we adopt an
empirical Bayes type approach to set $q$ to a large value determined from the
data and let the prior shrink the extra columns to zeros, substantially
simplifying the computation. The strategy is discussed in greater details
later in this section. We show in Section 2.6 that this approach is sufficient
for the recovery of true $\Omega$ under very mild conditions.

When $d\gg n$, $d$ distinct $\delta_{j}^{2}$’s also result in over-
parametrization. To reduce the number of parameters, we assume the
$\delta_{j}^{2}$’s to comprise a small number of unique values. To achieve
this in a data adaptive way, we use a Dirichlet process (DP) prior (Ferguson,
1973) on the $\delta_{j}^{2}$’s as

$\displaystyle\delta_{j}^{2}\lvert
G\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}G,~{}~{}~{}G\lvert\alpha\sim\mathrm{DP}(\alpha,G_{0})\text{
with
}G_{0}=\hbox{Ga}(a_{\delta},b_{\delta}),~{}~{}~{}\alpha\sim\hbox{Ga}(a_{\alpha},b_{\alpha}),$
(6)

where $\alpha$ is the concentration parameter and $G_{0}$ is the base measure.
Samples drawn from a DP are almost surely discrete, inducing a clustering of
the $\delta_{j}^{2}$’s and thus reducing the dimension. Integrating out $G$,
model (6) leads to a recursive Polya urn scheme illustrative of the clustering
mechanism while also being convenient for posterior computation (Escobar and
West, 1995; Neal, 2000). Specifically, we have

$\displaystyle\textstyle\delta_{j+1}^{2}\lvert\alpha,\mbox{\boldmath$\delta$}_{1:j}^{2}\propto\sum_{\ell=1}^{k_{j}}d_{j,\ell}\mathbbm{1}(\delta_{j+1}^{2}=\delta_{\ell}^{*2})+\alpha
G_{0}(\delta_{j+1}^{2}),$

where $\\{\delta_{1}^{*2},\dots,\delta_{k_{j}}^{*2}\\}$ denote the unique
values among $\mbox{\boldmath$\delta$}_{1:j}^{2}$ and
$d_{j,r}=\sum_{\ell=1}^{j}\mathbbm{1}(c_{\ell}=r)$ denote their
multiplicities, $c_{1},\dots,c_{d}$ being the latent variables such that
$c_{j}=r$ if $\delta_{j}^{2}$ belongs to the $r^{th}$ cluster.

Choice of Hyperparameters: To specify the value of the latent dimension $q$,
we adopt a principal component analysis (PCA) based empirical Bayes type
approach (Bai and Ng, 2008). First, we perform a sparse PCA (Baglama and
Reichel, 2005) on the data matrix ${\mathbf{y}}_{1:n}$ and compute the
reciprocals of the singular values. We set $q$ to be the number of inverse
singular values in decreasing order that adds up to $95\%$ of the total sum.
We set $a=0.5$ and $b=2.0$ in all our simulation studies and real data
applications. For the prior on the residual variances and the DP concentration
parameter, we set $a_{\delta}=b_{\delta}=a_{\alpha}=b_{\alpha}=0.1$.

###  2.4 Posterior Computation

The latent factor construction discussed in Section 2.1 leads to a novel,
elegant Gibbs sampler that is operationally simple and free of any tuning
parameter. There is also no need for additional constraints to ensure positive
definiteness which used to be a major setback for MCMC based methods for
precision matrix estimation. These features allow the sampler to be applied to
dimensions far beyond the reach of the current state-of-the-art. Scalable
posterior computation in LRD decomposed covariance matrix models is pretty
standard in the literature but remained a major barrier for precision matrix
models. The Gibbs sampler described here addresses this significant gap in the
literature.

Starting with some initial values of
$\mbox{\boldmath$\Lambda$},\mbox{\boldmath$\Delta$}$ and other parameters, our
sampler iterates between the following steps. Other parameters and
hyperparameters being implicitly understood in the conditioning, Step 1
defines a transition for
${\mathbf{u}},{\mathbf{v}}|{\mathbf{y}},\mbox{\boldmath$\Lambda$},\mbox{\boldmath$\Delta$}$;
Step 2 for $\mbox{\boldmath$\Lambda$}|{\mathbf{u}},{\mathbf{v}}$ (which is
identical to
$\mbox{\boldmath$\Lambda$}|{\mathbf{y}},{\mathbf{u}},{\mathbf{v}},\mbox{\boldmath$\Delta$}$);
Step 3 for $\mbox{\boldmath$\Delta$}|{\mathbf{v}}$ (which is identical to
$\mbox{\boldmath$\Delta$}|{\mathbf{y}},{\mathbf{u}},{\mathbf{v}},\mbox{\boldmath$\Lambda$}$);
and Steps 4 and 5 are standard updates for the parameters and hyper-parameters
of the DL and DP priors, respectively.

Step 1

Generate
${\mathbf{u}}_{1},\dots,{\mathbf{u}}_{n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{P}})$
with ${\mathbf{P}}=({\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$})$ independently from
${\mathbf{y}}_{1:n}$ and let
${\mathbf{v}}_{i}={\mathbf{y}}_{i}+\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{i}$.

Step 2

We have
${\mathbf{u}}_{i}=\sum_{r=1}^{d}\mbox{\boldmath$\lambda$}_{r}v_{r,i}+\mbox{\boldmath$\varepsilon$}_{i}$,
where $\mbox{\boldmath$\lambda$}_{r}=(\lambda_{r,1},\dots,\lambda_{r,q})$ is
the $r^{th}$ row of $\Lambda$ and
${\mathbf{v}}_{i}=(v_{1,i},\dots,v_{d,i})^{\rm T}$. Define
${\mathbf{u}}_{i}^{(j)}={\mathbf{u}}_{i}-\sum_{r\neq
j}\mbox{\boldmath$\lambda$}_{r}v_{r,i}$. Then
${\mathbf{u}}_{i}^{(j)}=\mbox{\boldmath$\lambda$}_{j}v_{j,i}+\mbox{\boldmath$\varepsilon$}_{i}$.
Conditioned on ${\mathbf{u}}_{i}^{(j)}$, ${\mathbf{v}}_{i}$ and the associated
hyper-parameters, $\mbox{\boldmath$\lambda$}_{j}$’s can be updated
sequentially for $j=1,\dots,d$ from the distribution

$\displaystyle\mbox{\boldmath$\lambda$}_{j}\sim\mathrm{N}_{q}\\{({\mathbf{D}}_{j}^{-1}+\norm{{\mathbf{v}}^{(j)}}^{2}{\mathbf{I}}_{q})^{-1}{\mathbf{w}}_{j},({\mathbf{D}}_{j}^{-1}+\norm{{\mathbf{v}}^{(j)}}^{2}{\mathbf{I}}_{q})^{-1}\\},$

where
${\mathbf{D}}_{j}=\tau^{2}\mathrm{diag}\left(\psi_{j,1}\phi^{2}_{j,1},\dots,\psi_{j,q}\phi^{2}_{j,q}\right)$,
${\mathbf{v}}^{(j)}=(v_{j,1},\dots,v_{j,n})^{\rm T}$ and
${\mathbf{w}}_{j}=\sum_{i=1}^{n}v_{j,i}{\mathbf{u}}_{i}^{(j)}$.

Step 3

Sample the $\delta_{j}^{2}$’s through the following steps.

  1. (i)

Let $d_{r,-j}=\sum_{\ell\neq j}\mathbbm{1}(c_{\ell}=r)$ and
${\mathbf{v}}^{(-j)}$ to be the collection of all ${\mathbf{v}}^{(\ell)}$’s,
$\ell=1,\dots,d$, excluding ${\mathbf{v}}^{(j)}$. For $j=1,\dots,d$, sample
the cluster indicators sequentially from the distribution

$\displaystyle
p(c_{j}=r)\propto\begin{cases}d_{r,-j}\int\mathrm{N}({\mathbf{v}}^{(j)};0,\delta_{r}^{*-2})\mathrm{d}G_{0}\left(\delta_{r}^{*2}\lvert{\mathbf{v}}^{(-j)}\right)\text{
for }r\in\\{c_{\ell}\\}_{\ell\neq j};\\\
\alpha\int\mathrm{N}({\mathbf{v}}^{(j)};0,\delta_{r}^{*-2})\mathrm{d}G_{0}(\delta_{r}^{*2})\text{
for }r\neq c_{\ell}\text{ for all }\ell\neq j.\end{cases}$

The above integrals are analytically available and involves the density of a
multivariate central Student’s $t$-distribution for
$G_{0}=\hbox{Ga}(a_{\delta},b_{\delta})$.

  2. (ii)

Let the unique values in ${\mathbf{c}}_{1:d}$ be {$1,\dots,k\\}$. For
$r=1,\dots,k$, set $d_{r}=\sum_{j}\mathbbm{1}(c_{j}=r)$ and
${\mathbf{V}}_{r}=\sum_{j:c_{j}=r}\norm{{\mathbf{v}}^{(j)}}^{2}$, and
independently sample
$\delta_{r}^{*2}\sim\hbox{Ga}\left(a_{\delta}+{nd_{r}}/{2},b_{\delta}+{{\mathbf{V}}_{r}}/{2}\right)$.

  3. (iii)

Set $\delta_{j}^{2}=\delta_{c_{j}}^{*2}$.

Step 4

Sample the hyper-parameters in the priors on $\Lambda$ through the following
steps.

  1. (i)

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

  2. (ii)

Sample the full conditional posterior distribution of $\tau$ from a
generalized inverse Gaussian
$\mbox{giG}\left\\{dq(1-a),2b,2\sum_{j,h}{\absolutevalue{\lambda_{j,h}}}/{\phi_{j,h}}\right\\}$
distribution.

  3. (iii)

Draw $T_{j,h}$ independently with
$T_{j,h}\sim\mbox{giG}(a-1,1,2\absolutevalue{\lambda_{j,h}})$ and set
$\phi_{j,h}=T_{j,h}/T$ with $T=\sum_{j,h}T_{j,h}$.

Step 5

Following West (1992), first generate $\varphi\sim\hbox{Beta}(\alpha+1,d)$,
evaluate
$\pi/(1-\pi)=(a_{\alpha}+k-1)/\left\\{d(b_{\alpha}-\log\varphi)\right\\}$ and
then generate

$\displaystyle\alpha\lvert\varphi,k\sim\begin{cases}\hbox{Ga}(\alpha+k,b_{\alpha}-\log\varphi)\text{
with probability }\pi,\\\ \hbox{Ga}(\alpha+k-1,b_{\alpha}-\log\varphi)\text{
with probability }1-\pi.\end{cases}$

######  Remark 1.

Note that the main strategies underlying the algorithm above are not specific
to the DL prior considered here. We are free to choose any other shrinkage
prior that admits a conditionally Gaussian hierarchical representation for the
entries of $\Lambda$ and modify Steps 2 and 4 accordingly. This is still a
very large class of priors (Polson and Scott, 2010), including, e.g.,
horseshoe (Carvalho et al., 2009), multiplicative gamma (Bhattacharya and
Dunson, 2011), etc. For non-conjugate priors, Steps 3 and 4 can also be
modified with appropriate Metropolis-Hastings schemes. The other steps remain
the same, making it a very broadly adaptable algorithm. An MCMC scheme for
generic priors is outlined in Section S.1 of the supplementary materials.

######  Remark 2.

In the above sampler, the conditional posterior covariance matrix of
$\mbox{\boldmath$\lambda$}_{j}$ is diagonal which allows us to update it with
linear complexity. Thus, in each MCMC iteration, we only need a single small
dimensional $q\times q$ order matrix factorization operation in Step 1 to
simulate ${\mathbf{u}}_{1:n}$. While high-dimensional matrix factorization
operations are usally numerically very expensive, we are able to completely
avoid that, facilitating substantial scalability.

######  Remark 3.

All but Steps 1 and 3 can be divided into parallel operations in a
straightforward manner. Additionally, in recent versions of many popular
statistical software, including R, matrix operations are inherently
parallelized and hence Step 1 is also highly scalable in any decent computing
system.

We implemented the Gibbs sampler in C++ and ported to R using the Rcpp package
(Eddelbuettel and Francois, 2011). In each case considered in this article,
synthetic or real, we ran $5,500$ iterations which takes approximately 37
minutes on a system with an i9-10900K CPU and 64GB memory for a $d=1,000$
dimensional problem with sample size $n=1,000$. The initial $1,250$ samples
were discarded as burn-in and the remaining samples were thinned by an
interval of $5$. In all our experiments, convergence was swift and mixing was
excellent. A comparison of the runtimes of our method and a few other existing
methods is presented in Figure 5(c) below.

###  2.5 Graph Selection

As discussed in Section 2.2.2, the off-diagonals of $\Omega$ are penalized by
inducing shrinkage on $\Lambda$. The theoretical results in Section 2.6 and
numerical experiments in Section 3 show that, with our carefully constructed
data adaptive shrinkage priors on $\Lambda$, the inferred sparsity patterns in
$\Lambda$ are such that a sparse $\Omega$ is also accurately recovered (see
Figures S.1 and 4(a) for estimates of sparse precision matrices obtained by
our method). Exact zero estimates are, however, not obtained even for the
insignificant off-diagonal elements of $\Omega$ for finite samples. This is an
artifact of continuous shrinkage priors, since the probabilities of exact
zeroes are almost surely null for finite samples although the posterior
probabilities of arbitrary sets around zeroes are very high.

We address the issue of non-zero edge selection through a novel multiple
hypothesis testing based approach. For $i=1,\dots,d,~{}j=i+1,\dots,d$ and some
$\epsilon>0$, we consider testing

$\displaystyle H_{0,i,j}:\absolutevalue{\rho_{i,j}}\leq\epsilon~{}~{}~{}\text{
versus }~{}~{}~{}H_{1,i,j}:\absolutevalue{\rho_{i,j}}>\epsilon,$ (7)

where $\rho_{i,j}$ is the $(i,j)^{th}$ element of
$\mathrm{diag}(\mbox{\boldmath$\Omega$})^{-\nicefrac{{1}}{{2}}}~{}\mbox{\boldmath$\Omega$}~{}\mathrm{diag}(\mbox{\boldmath$\Omega$})^{-\nicefrac{{1}}{{2}}}$,
the partial correlation matrix derived from $\Omega$. Here we follow Berger
(1985, Chapter 4, pp. 148) in replacing the point nulls
$H_{0,i,j}:\rho_{i,j}=0$ by reasonable interval nulls
$H_{0,i,j}:\absolutevalue{\rho_{i,j}}\leq\epsilon$. If $H_{0,i,j}$ is rejected
in favor of $H_{1,i,j}$, we conclude that there is an edge between nodes $i$
and $j$.

We utilize posterior uncertainty to resolve these testing problems.
Specifically, we define
$d_{i,j}=\mathbbm{1}\left\\{\Pi(H_{1,i,j}|{\mathbf{y}}_{1:n})>\beta\right\\}$
as the decision rule which controls the posterior FDR defined as

$\displaystyle\hbox{FDR}_{{\mathbf{y}}}=\frac{\sum_{i,j}d_{i,j}\Pi(H_{0,i,j}|{\mathbf{y}}_{1:n})}{\max(\sum_{i,j}d_{i,j},1)},$

at the level $1-\beta$. Importantly, the decision rule also incurs the lowest
false non-discovery rate (Müller et al., 2004). For a fixed $\beta$, the
$\hbox{FDR}_{{\mathbf{y}}}$ depends on the choice of $\epsilon$. To obtain the
optimal $\epsilon$, we compute the $\hbox{FDR}_{{\mathbf{y}}}$’s on a grid of
$\epsilon$ values in $(0,1)$ and then set
$\epsilon=\inf_{\epsilon^{\prime}}\mathrm{FDR}_{{\mathbf{y}}}(\epsilon^{\prime})\leq
1-\beta$. This way, we have the FDR, that is, a quantified statistical
uncertainty associated with the estimated graph. In all simulation experiments
and real data applications in this paper, we control
$\hbox{FDR}_{{\mathbf{y}}}$ at the $0.10$ level of significance.

######  Remark 4.

Although this FDR control procedure is widely applicable, efficient posterior
exploration is crucial for this step and hence may not be achievable or
scalable to high dimensional problems when adapted to previously existing
Bayesian approaches.

###  2.6 Posterior Concentration

Preliminaries and Notation: We let $\norm{{\mathbf{A}}}_{F}$ and
$\norm{{\mathbf{A}}}_{2}$ denote the Frobenius and spectral norms of a matrix
${\mathbf{A}}$, respectively; $s^{2}_{\min}{({\mathbf{A}})}$ be the smallest
singular value of ${\mathbf{A}}^{\rm T}{\mathbf{A}}$; and
$\lambda_{1}({\mathbf{A}}),\dots,\lambda_{d}({\mathbf{A}})$ be the eigenvalues
of ${\mathbf{A}}$ in decreasing order when ${\mathbf{A}}$ is a $d$-dimensional
diagonalizable matrix. Also, $a_{n}=o(b_{n})$ and $a_{n}=O(b_{n})$ imply that
$\lim\absolutevalue{{a_{n}}/{b_{n}}}=0$ and
$\limsup\absolutevalue{{a_{n}}/{b_{n}}}<\infty$, respectively. Throughout,
$C,C^{\prime}$ and $\widetilde{C}$ are used to denote positive constants whose
values might change from one line to the next but are independent from
everything else. A $d$-dimensional vector $\theta$ is said to be $s$-sparse if
only $s$ among the $d$ elements of $\theta$ are non-zero. We denote the set of
all $s$-sparse vectors in $\mathbb{R}^{d}$ by $\ell_{0}\left[s,d\right]$.

We allow the model parameters to increase in dimensions with sample size $n$,
indicated by associating them with the suffix $n$. We let $\Pi_{n}(\cdot)$
denote the prior and $\Pi_{n}(\cdot\lvert{\mathbf{y}}_{1:n})$ the
corresponding posterior given data ${\mathbf{y}}_{1:n}$, respectively.

Assumptions on the Data Generating Process: We assume that
${\mathbf{y}}_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Omega$}_{0n}^{-1}),i=1,\dots,n$.
In Section 2.2.1 we discussed how an LRD decomposition of a sparse
$\mbox{\boldmath$\Omega$}_{0n}$ can be constructed. Hence, we assume that
$\mbox{\boldmath$\Omega$}_{0n}$ admits the factor representation
$\mbox{\boldmath$\Omega$}_{0n}=\mbox{\boldmath$\Lambda$}_{0n}\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}+\delta_{0n}^{2}{\mathbf{I}}_{d_{n}}$, where
$\mbox{\boldmath$\Lambda$}_{0n}$ is a $d_{n}\times q_{0n}$ order sparse matrix
and $\delta_{0n}^{2}>0$ is a scalar. To simplify the theoretical analysis, we
deviate here slightly from the proposed model and assume that the
$\delta^{2}$’s all come from a single cluster with the common value
$\delta_{0n}^{2}$. We show that the precision factor model with an appropriate
DL prior can recover $\mbox{\boldmath$\Omega$}_{0n}$ in ultra high-dimensional
settings. As discussed in Section 2.3, we fix $q_{n}$ to a liberal large value
and let the prior shrink the extra columns to zeros. We show that with
appropriate sparsity conditions, the posterior of the precision factor model
with the DL prior then concentrates around $\mbox{\boldmath$\Omega$}_{0n}$,
even when the data dimension $d_{n}$ increases in exponential order with $n$.
The requisite conditions are stated below.

  1. (C1)

Let $\\{q_{0n}\\}_{n=1}^{\infty}$ and $\\{s_{n}\\}_{n=1}^{\infty}$ be
increasing sequences of positive integers such that $s_{n}^{2}=O(\log d_{n})$,
$q_{0n}=O\\{s_{n}\log(d_{n}q_{0n})\\}$ and
$s_{n}q_{0n}\log(d_{n}q_{0n})=o(n)$.

  2. (C2)

$\mbox{\boldmath$\Lambda$}_{0n}$ is a $d_{n}\times q_{0n}$ order full rank
matrix such that each column of $\mbox{\boldmath$\Lambda$}_{0n}$ belongs to
$\ell_{0}\left[s_{n},d_{n}\right]$, $\liminf
s_{\min}\left(\mbox{\boldmath$\Lambda$}_{0n}\right)>0$ and
$\norm{\mbox{\boldmath$\Lambda$}_{0n}}_{2}=O(1)$.

  3. (C3)

The scalar $\delta_{0n}^{2}$ lies in some compact set
$[\delta_{\min}^{2},\delta_{\max}^{2}]$.

Specifics of the Postulated Precision Factor Model: Let
$\\{q_{n}\\}_{n=1}^{\infty}$ be an increasing sequence of positive integers
such that $q_{0n}\leq q_{n}=o(n)$ and $s_{n}q_{0n}\log(d_{n}q_{n})=o(n)$. We
assume that
${\mathbf{y}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Omega$}_{n}^{-1})$,
where
$\mbox{\boldmath$\Omega$}_{n}=\mbox{\boldmath$\Lambda$}_{n}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}+\delta_{n}^{2}{\mathbf{I}}_{d_{n}}$ and $\mbox{\boldmath$\Lambda$}_{n}$ is
a $d_{n}\times q_{n}$ matrix. We consider a $\mathrm{DL}(a_{n},b_{n})$ prior
on $\mathrm{vec}(\mbox{\boldmath$\Lambda$}_{n})$ with $a_{n}={1}/{d_{n}q_{n}}$
and $b_{n}=\log^{3/2}(d_{n}q_{n})$. We assume a gamma prior on
$\delta_{n}^{2}\sim\hbox{Ga}(a_{\delta},b_{\delta})$ independent of
$\mbox{\boldmath$\Lambda$}_{n}$ and truncated to the compact set
$[\delta_{\min}^{2},\delta_{\max}^{2}]$.

Condition (C1) specifies the requisite sparsity conditions. It also imposes a
condition on $d_{n}$. Specifically, it can be seen that $d_{n}$ can be of
exponential order of $n$, allowing the recovery of massive precision matrices
based on relatively small sample sizes. Note that, subject to appropriate
orthogonal transformation, marginally
$\hbox{var}(y_{j})=1/\left\\{\lambda_{j}\left(\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Lambda$}_{0n}\right)+\delta_{0n}^{2}\right\\}$. Conditions
(C2) and (C3) ensure that these variances lie in a compact set. The prior
specification provides a set of sufficient conditions on the class of proposed
models. The true number of latent factors $q_{0n}$ is also assumed smaller
than the number of latent factors $q_{n}$ in the postulated models.

We let $\mathcal{P}_{0n}$ denote the class of precision matrices satisfying
(C1)-(C3) and $\mathcal{P}_{n}$ denote the class of positive definite matrices
parametrized by $(q_{n},\delta_{n}^{2},\mbox{\boldmath$\Lambda$}_{n})$ in our
precision factor model. Notably, $\mathcal{P}_{0n}$ can also be parametrized
by $(q_{0n},\delta_{0n}^{2},\mbox{\boldmath$\Lambda$}_{0n})$. However, since
$q_{0n}\leq q_{n}$, for a fixed $q_{n}$, $\mathcal{P}_{n}$ may be
overparametrized for $\mathcal{P}_{0n}$. In Theorem 1, we show that with
increasing sample size, the posterior distribution supported on
$\mathcal{P}_{n}$ still concentrates around the true
$\mbox{\boldmath$\Omega$}_{0n}$ from $\mathcal{P}_{0n}$. Using general results
from Ghosal and van der Vaart (2017), we establish a convergence rate in terms
of the operator norm. A rigorous proof is detailed in Section S.2 in the
supplementary materials.

######  Theorem 1.

For any $\mbox{\boldmath$\Omega$}_{0n}\in\mathcal{P}_{0n}$,
$\epsilon_{n}=(s_{n}q_{0n})^{4}\log(d_{n}q_{n})\sqrt{q_{0n}s_{n}\log(d_{n}q_{n})/n}$
and any sequence $M_{n}\to\infty$,

$\displaystyle\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{0n}}\Pi_{n}\left(\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}>M_{n}\epsilon_{n}\lvert{\mathbf{y}}_{1:n}\right)=0.$

In Section 2.2.1 we discussed a constructive way to find a sparse LRD
decomposition of arbitrary sparse $\mbox{\boldmath$\Omega$}_{0n}$ where the
column dimension $q_{0n}$ of $\mbox{\boldmath$\Lambda$}_{0n}$ need not exceed
the number of non-zero off-diagonals of $\mbox{\boldmath$\Omega$}_{0n}$ or
edges in the conditional dependence graph. Theorem 1 implies that the
precision matrix is recoverable if (C1) holds among others. Note that (C1)
implies $q_{0n}=o(n)$. Hence, $\mbox{\boldmath$\Omega$}_{0n}$ can be learned
from the data using the precision factor model as long as the number of edges
is bounded by the sample size. Such assumptions are required to recover
massive dimensional precision matrices from relatively smaller amount of data
(Meinshausen and Bühlmann, 2006; Ksheera Sagar et al., 2021).

######  Remark 5.

We derive a minimax lower bound of $\sqrt{s_{n}\log d_{n}/n}$ for the
precision matrix estimation problem in Theorem 4 in the supplementary
materials when $q_{0n}\leq q_{n}=O(1)$ which is a special case of (C1). The
rate in Theorem 1 is ${s_{n}^{4}\log d_{n}}$ times the lower bound. If we
further let $d_{n}=O(n^{r})$ for some fixed positive integer $r$, the rate
attains the minimax lower bound up to a $(\log n)^{3}$ term.

######  Remark 6.

An interesting implication of our theoretical results is the robustness with
respect to the choice of the latent dimension $q$ under very mild conditions.
For the class of our postulated precision factor models, we err on the side of
overestimating $q$ and let the DL prior shrink the extra factors to zero. The
theoretical results imply that this strategy is sufficient to recover
$\mbox{\boldmath$\Omega$}_{0n}$ efficiently. The precise recovery rate,
however, naturally depends on the choice of $q_{n}$, a key parameter that
distinguishes the class of postulated models from the class of true models
(Shalizi, 2009). The closer $q_{n}$ is to $q_{0n}$, the smaller the space to
search for the truth, and the better the rate.

##  3 Simulation Studies

In this section, we discuss the results of some synthetic numerical
experiments. We evaluate the performance of estimating the precision matrix
$\Omega$ itself as well as the underlying graph. We simulate data from three
different cases - (i) a Gaussian autoregressive (AR) process of order 2, (ii)
a multivariate Gaussian distribution with banded precision matrix, and (iii) a
multivariate Gaussian distribution with randomly generated arbitrarily
structured sparse (RSM) precision matrix. In all these cases, we take the mean
vector to be zero. For plotting purposes in the RSM case, we assume two nodes
to have an edge between them if the absolute value of the associated partial
correlation exceeds $0.1$ and refer to this as the ‘true’ graph. We perform
experiments for dimensions $d=$ 50, 100, 200 and 1,000. For $d\leq 200$ we
take $n=100$ and for $d=1,000$ we take $n=1,000$. We plot the true precision
matrices and the associated true graphs for $d=50$ in Figure S.1 in the
supplementary materials and in Figure 6(a) here in the main paper. For $d\leq
200$ and $d=1,000$ we consider 50 and 20 independent replications for each
scenario, respectively.

We apply our precision factor (PF) model to recover the precision matrices,
using the posterior mean as our Bayesian point estimate. We compare it with
the methods Bayesian graphical model under shrinkage (Bagus, Gan et al.,
2019), the graphical Lasso (Glasso, Friedman et al., 2008) and the
neighborhood selection by Meinshausen and Bühlmann (M&B, 2006). We do not
consider any previously existing Bayesian posterior sampling based methods
here as they do not scale well beyond only a few tens of dimensions. We assess
the comparative efficacy of these methods using qualitative graphical
summaries as well as quantitative performance measures. For the proposed PF
method, the graphs are estimated following the FDR control procedure outlined
in Section 2.5. For Bagus, we follow the prescription in Section 4.3 of Gan et
al. (2019). For $d=1,000$ in the banded case, the codes for Bagus did not work
and we could not consider it as a competitor in that particular setup. For the
other methods, the non-zero entries in the estimated precision matrix are
considered as edges.

We compute the Frobenius norm between the true and estimated precision
correlation matrices. To assess the accuracy of the derived graphs, we compute
specificity $=\frac{TN}{TN+FP}$ and sensitivity $=\frac{TP}{TP+FN}$, where
$TP$ (true positives), $FP$ (false positives), $TN$ (true negatives) and $FN$
(false negatives) are based on the detection of the edges in the estimated
graphs and comparing them with the corresponding true graphs. In Figure 5, the
average values of these criteria across the replications are plotted for the
competing methods. We also compare the execution times of the different
methods.

From Figure 5, we see that, for banded precision matrices, the Frobenius norms
obtained by our proposed PF method are very similar to those produced by the
competitors. In the case of the arbitrarily structured precision matrix,
however, the PF method performs much better, especially in higher dimensions.
Also, the PF method yields much higher sensitivity in almost all the scenarios
at the expense of slightly smaller specificity in some cases. This is expected
since we are allowing a small margin of error by controlling the FDR at
$0.10$. Notably, the PF method is much more powerful in detecting the true
edges. The other methods seem to be quite conservative in that regard, as seen
in Figure 5(b). This explains the high Frobenius norm produced by the PF
method in the AR(2) case with $d=1,000$. The true precision matrix being
extremely sparse with only two off-diagonals being non-zero, the conservative
competitors are performing better in recovering the precision matrix in this
particular setup. We see in Figure 5(c) that the PF method is slower than M&B
and Glasso but is much faster than Bagus. Note however that, unlike the
competitors, we are exploring the full posterior, not just providing a point
estimate. On a related important note, since we are able to sample from the
full posterior, unlike the other methods, we are able to provide a natural way
of quantifying posterior uncertainty. Specifically, posterior credible
intervals for each element of the precision matrix can be obtained from the
MCMC samples. For $d=50$ and sample size $n=100$, we plot the lower $2.5\%$
and upper $97.5\%$ quantiles of the entry-wise partial correlations along with
the simulation truths in Figure S.1 in the supplementary materials. The plots
indicate that the simulation truths are well-within the confidence bounds.

Circos plots (Gu et al., 2014) of the true and estimated graphs are shown in
Figure 7. For the AR(2) case, none of the methods are doing well in recovering
the graph. For the other cases, the proposed PF method outperforms the
competitors. Figures S.3 and S.5 in the supplementary materials provide
alternate graphical representations useful in visually discerning the superior
performance of the PF method.

(a) Frobenius norm

(b) Sensitivity and specificity

(c) Execution time

Figure 5: Results of simulation experiments: Panel (a) shows the Frobenius
norms between the true and estimated partial correlation matrices; panel (b)
shows the sensitivity and specificity; and panel (c) shows the execution times
in seconds.

(a) True graphs.

(b) Estimated graphs by the PF method.

(c) Estimated graphs by Bagus.

(d) Estimated graphs by M&B.

(a) Estimated graphs by Glasso. Figure 7: Results of simulation experiments:
Graph recovery: Panel (a) shows the true graphs for AR(2), banded and RSM
structures from left to right; panels (b), (c), (d) and (e) show the
corresponding estimated graphs for our proposed PF and Bagus, M&B and Glasso
methods, respectively. Positive (negative) associations are represented by
blue (red) links, their opacities being proportional to the corresponding
association strengths. The link widths are inversely proportional to the
number of edges associated with the corresponding nodes. Figure S.3 in the
supplementary materials provides a useful alternative graphical
representation.

##  4 Application

We applied our proposed method to two real data sets from two different
application domains, namely genomics and finance. To meet space constraints,
the finance application is presented separately in Section S.4 in the
supplementary materials.

Immune cells serve specialized roles in innate and adaptive bodily responses
to eliminate antigens. To understand the cell biology of carcinogenic
processes, study of immune cells and the genes therein are thus of immense
importance.

(a) Graph of all genes.

(b) Heatmap of the connected genes.

(c) Graph of the connected genes.

Figure 8: Results for ImmGen microarray data: Panel (a) shows the graph for
all genes; panel (b) shows the heatmap of the partial correlation matrix for
the genes having at least 5 edges; panel (c) zooms into the graph of these
genes. Positive (negative) associations are represented by blue (red) links,
their opacities being proportional to the corresponding association strengths.
The link widths are inversely proportional to the number of edges associated
with the corresponding nodes.

We obtain the data from the Immunological Genome Project (ImmGen) data browser
(Heng et al., 2008). ImmGen is a collaborative scientific research project
that is currently building a gene-expression database for all characterized
immune cells in mice. In particular, we use the GSE15907 microarray dataset
(Painter et al., 2011; Desch et al., 2011) comprising of multiple immune cell
lineages which were isolated ex-vivo, primarily from young adult B6 male mice
and double-sorted to $>99\%$ purity. The cell population includes all adaptive
and innate lymphocytes (B, abT, gdT, Innate-Like Lymphocytes), myeloid cells
(dendritic cells, macrophages, monocytes), mast cells and neutrophils. The
already normalized dataset has more than $21,000$ gene expressions from
$n=653$ immune cells. We made a $\log_{2}$ transformation of the data and
filtered the top $2.5\%$ genes with highest variances using the genefilter R
package (Gentleman et al., 2020), resulting in a $d=544$ dimensional problem.
Since different cell-types exhibit very different gene expression profiles, we
centered the gene-expressions separately within each cell type. The estimated
graph corresponding to all genes is shown in Figure 8(a). We can see clear
evidence towards conditional independence between most of the genes,
indicating that only a small subset of genes are functionally responsible for
the variability.

For more insights, we zoom into the connected genes and plot the corresponding
partial correlation matrix and graph, limited to the genes possessing at least
5 edges, in Figures 8(b) and 8(c) respectively. We notice overall positive
partial correlations between the histone class of genes such as Hist2h3b,
Hist1h3h, Hist2h3c1, Hist1h3c, etc (Wolffe, 2001). The positive correlations
between these protein coding genes indicate that these genes are generally
expressed together in the mechanism. Bcl2a1a and Bcl2a1d, two functional
isoforms of the B cell leukemia 2 family member A1 exerting important pro-
survival functions, show strong positive association. We also observe strong
positive association between Ly6c1 and Ly6c2 genes. Lee et al. (2013) noted
these genes to be located adjacent to each other in the mid‐section of the Ly6
complex. They share $>95\%$ similarity in their genomic and protein sequences
and these two genes has been considered synonymous to each other and hence
exhibit strong positive association. Although positive correlations are often
observed between the membrane-spanning 4A class of genes such as Ms4a4c,
Ms4a4b, Ms4a6b, etc. (Liang et al., 2001), we find the genes to be
conditionally almost independent in our analysis. It indicates that their
expressions might be regulated by other genes.

From Figure 8(b), it thus seems that there exist two blocks where within each
block the genes are positively associated whereas between the blocks the genes
are negatively associated. This might again be an indication that either of
these blocks of genes are generally expressed together in immune cells.

##  5 Discussion

In this article, we proposed a novel flexible statistical model for Gaussian
precision matrices that relies on decomposing them into a low-rank and a
diagonal component. The decomposition is theoretically and practically highly
flexible and is thus often used in covariance matrix models, arising naturally
in their computationally tractable latent factor based representations. The
approach has, however, not been popular in precision matrix and related graph
estimation problems as it poses daunting computational challenges when applied
to such settings. We addressed this issue in this article by exploring a
previously under-utilized latent variable construction, leading to a highly
scalable Gibbs sampler that allows efficient posterior inference. The
decomposition based strategy also allowed us to use sparsity inducing priors
to shrink insignificant off-diagonal entries toward zero while also making the
approach adaptable to high-dimensional sparse settings. We specifically
adapted the Dirichlet-Laplace prior for sparse precision matrix estimation. We
developed a novel posterior FDR control based method to perform graph
selection that properly accommodates posterior uncertainty. We also
established theoretical convergence guarantees for the proposed model in high-
dimensional sparse settings. In synthetic experiments, the proposed method
vastly outperformed its competitors in receiving the true underlying graphs.
We illustrated the method’s practical utility through real data examples from
genomics and finance.

Aside from providing fundamentally new probabilistic perspectives on Gaussian
precision matrix and related graph estimation problems, our work also bridges
the gap between frequentist penalized likelihood based strategies and Bayesian
shrinkage prior based ideas for such models. The simple and highly scalable
computational algorithms resulting from our latent factor representation
should also free up Bayesians from having to put restrictive assumptions on
the graph structures to achieve computational tractability, opening up new
opportunities to adapt the basic models to more complex and more realistic
data structures and study designs. A few such methodological extensions we are
pursuing as topics of separate ongoing research include dynamic Gaussian
graphical models (Huang and Chen, 2017), covariate dependent Gaussian
graphical models, nonparanormal (Liu et al., 2009) and Gaussian copula
graphical models (Pitt et al., 2006), etc.

## Supplementary Materials

Supplementary materials discuss a general strategy for posterior computation
under broad classes of generic priors, proofs of the theoretical results, some
additional figures, and an application to a NASDAQ-100 stock price dataset.

## References

  * Armstrong et al. (2009) Armstrong, H., Carter, C. K., Wong, K. F. K., and Kohn, R. (2009).  Bayesian covariance matrix estimation using a mixture of decomposable graphical models.  Statistics and Computing, 19, 303–316.
  * Atay-Kayis and Massam (2005) Atay-Kayis, A. and Massam, H. (2005).  A Monte Carlo method for computing the marginal likelihood in nondecomposable Gaussian graphical models.  Biometrika, 92, 317–335.
  * Baglama and Reichel (2005) Baglama, J. and Reichel, L. (2005).  Augmented implicitly restarted Lanczos bidiagonalization methods.  SIAM Journal on Scientific Computing, 27, 19–42.
  * Bai and Ng (2008) Bai, J. and Ng, S. (2008).  Large dimensional factor analysis.  Foundations and Trends in Econometrics, 3, 89–163.
  * Banerjee et al. (2008) Banerjee, O., El Ghaoui, L., and d’Aspremont, A. (2008).  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  JMLR, 9, 485–516.
  * Banerjee and Ghosal (2015) Banerjee, S. and Ghosal, S. (2015).  Bayesian structure learning in graphical models.  Journal of Multivariate Analysis, 136, 147–162.
  * Barvinok and Rudelson (2022) Barvinok, A. and Rudelson, M. (2022).  When a system of real quadratic equations has a solution.  Advances in Mathematics, 403, 108391.
  * Berger (1985) Berger, J. O. (1985).  Statistical decision theory and Bayesian analysis.  Springer series in statistics. Springer-Verlag, New York, 2nd edition.
  * 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.
  * Bhattacharya et al. (2016) Bhattacharya, A., Chakraborty, A., and Mallick, B. K. (2016).  Fast sampling with Gaussian scale mixture priors in high-dimensional regression.  Biometrika, 103, 985–991.
  * Bien and Tibshirani (2011) Bien, J. and Tibshirani, R. J. (2011).  Sparse estimation of a covariance matrix.  Biometrika, 98, 807–820.
  * Carvalho and Scott (2009) Carvalho, C. M. and Scott, J. G. (2009).  Objective Bayesian model selection in Gaussian graphical models.  Biometrika, 96, 497–512.
  * Carvalho et al. (2007) Carvalho, C. M., Massam, H., and West, M. (2007).  Simulation of hyper-inverse Wishart distributions in graphical models.  Biometrika, 94, 647–659.
  * 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 and Bhattacharya (2019) Chandra, N. K. and Bhattacharya, S. (2019).  Non-marginal decisions: A novel Bayesian multiple testing procedure.  Electronic Journal of Statistics, 13, 489–535.
  * Dallakyan and Pourahmadi (2020) Dallakyan, A. and Pourahmadi, M. (2020).  Fused-lasso regularized Cholesky factors of large nonstationary covariance matrices of longitudinal data.  arXiv:2007.11168.
  * Daniele et al. (2019) Daniele, M., Pohlmeier, W., and Zagidullina, A. (2019).  Sparse approximate factor estimation for high-dimensional covariance matrices.  arXiv:1906.05545.
  * d’Aspremont et al. (2008) d’Aspremont, A., Banerjee, O., and El Ghaoui, L. (2008).  First-order methods for sparse covariance selection.  SIAM Journal on Matrix Analysis and Applications, 30, 56–66.
  * Dawid and Lauritzen (1993) Dawid, A. P. and Lauritzen, S. L. (1993).  Hyper Markov laws in the statistical analysis of decomposable graphical models.  The Annals of Statistics, 21, 1272–1317.
  * Dellaportas et al. (2003) Dellaportas, P., Giudici, P., and Roberts, G. (2003).  Bayesian inference for nondecomposable graphical Gaussian models.  Sankhyā: The Indian Journal of Statistics, pages 43–55.
  * Desch et al. (2011) Desch, A. N., Randolph, G. J., et al. (2011).  CD103+ pulmonary dendritic cells preferentially acquire and present apoptotic cell–associated antigen.  Journal of Experimental Medicine, 208, 1789–1797.
  * Deshpande et al. (2019) Deshpande, S. K., Ročková, V., and George, E. I. (2019).  Simultaneous variable and covariance selection with the multivariate spike-and-slab lasso.  Journal of Computational and Graphical Statistics, 28, 921–931.
  * Dobra et al. (2004) Dobra, A., Hans, C., Jones, B., Nevins, J. R., Yao, G., and West, M. (2004).  Sparse graphical models for exploring gene expression data.  Journal of Multivariate Analysis, 90, 196–212.
  * Dobra et al. (2011) Dobra, A., Lenkoski, A., and Rodriguez, A. (2011).  Bayesian inference for general Gaussian graphical models with application to multivariate lattice data.  Journal of the American Statistical Association, 106, 1418–1433.
  * Eddelbuettel and Francois (2011) Eddelbuettel, D. and Francois, R. (2011).  Rcpp: Seamless R and C++ integration.  Journal of Statistical Software, 40, 1–18.
  * Escobar and West (1995) Escobar, M. D. and West, M. (1995).  Bayesian density estimation and inference using mixtures.  Journal of the American Statistical Association, 90, 577–588.
  * Fan et al. (2011) Fan, J., Liao, Y., and Mincheva, M. (2011).  High-dimensional covariance matrix estimation in approximate factor models.  The Annals of Statistics, 39, 3320–3356.
  * Fan et al. (2018) Fan, J., Liu, H., and Wang, W. (2018).  Large covariance estimation through elliptical factor models.  The Annals of Statistics, 46, 1383–1414.
  * Ferguson (1973) Ferguson, T. S. (1973).  A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1, 209–230.
  * Friedman et al. (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9, 432–441.
  * Gan et al. (2019) Gan, L., Narisetty, N. N., and Liang, F. (2019).  Bayesian regularization for graphical models with unequal shrinkage.  Journal of the American Statistical Association, 114, 1218–1231.
  * Gentleman et al. (2020) Gentleman, R., Carey, V., Huber, W., and Hahne, F. (2020).  genefilter: methods for filtering genes from high-throughput experiments.  R package version 1.70.0.
  * Ghosal and van der Vaart (2017) Ghosal, S. and van der Vaart, A. (2017).  Fundamentals of nonparametric Bayesian inference.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Green (1995) Green, P. J. (1995).  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika, 82, 711–732.
  * Green and Thomas (2013) Green, P. J. and Thomas, A. (2013).  Sampling decomposable graphs using a Markov chain on junction trees.  Biometrika, 100, 91–110.
  * 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., Elpek, K., Lukacs-Kornek, V., Mauermann, N., Turley, S. J., Koller, D., Kim, F. S., Wagers, A. J., Asinovski, N., et al. (2008).  The immunological genome project: networks of gene expression in immune cells.  Nature Immunology, 9, 1091–1094.
  * Huang and Chen (2017) Huang, F. and Chen, S. (2017).  Learning dynamic conditional Gaussian graphical models.  IEEE Transactions on Knowledge and Data Engineering, 30, 703–716.
  * 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.
  * Jones et al. (2005) Jones, B., Carvalho, C., Dobra, A., Hans, C., Carter, C., and West, M. (2005).  Experiments in stochastic computation for high-dimensional graphical models.  Statistical Science, 20, 388–400.
  * Kang and Deng (2020) Kang, X. and Deng, X. (2020).  An improved modified Cholesky decomposition approach for precision matrix estimation.  Journal of Statistical Computation and Simulation, 90, 443–464.
  * Kastner (2019) Kastner, G. (2019).  Sparse Bayesian time-varying covariance estimation in many dimensions.  Journal of Econometrics, 210, 98–115.
  * Khare et al. (2018) Khare, K., Rajaratnam, B., and Saha, A. (2018).  Bayesian inference for Gaussian graphical models beyond decomposable graphs.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, 80, 727–747.
  * Khondker et al. (2013) Khondker, Z. S., Zhu, H., Chu, H., Lin, W., and Ibrahim, J. G. (2013).  The Bayesian covariance lasso.  Statistics and its Interface, 6, 243.
  * Koller and Friedman (2009) Koller, D. and Friedman, N. (2009).  Probabilistic graphical models: Principles and techniques.  MIT press.
  * 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.
  * Lauritzen (1996) Lauritzen, S. L. (1996).  Graphical models.  Clarendon Press.
  * Lee et al. (2013) Lee, P. Y., Wang, J.-X., et al. (2013).  Ly6 family proteins in neutrophil biology.  Journal of Leukocyte Biology, 94, 585–594.
  * Lenkoski (2013) Lenkoski, A. (2013).  A direct sampler for G-Wishart variates.  Stat, 2, 119–128.
  * Levina et al. (2008) Levina, E., Rothman, A., and Zhu, J. (2008).  Sparse estimation of large covariance matrices via a nested Lasso penalty.  The Annals of Applied Statistics, 2, 245–263.
  * Li et al. (2019a) Li, Y., Craig, B. A., and Bhadra, A. (2019a).  The graphical horseshoe estimator for inverse covariance matrices.  Journal of Computational and Graphical Statistics, 28, 747–757.
  * Li et al. (2019b) Li, Z., Mccormick, T., and Clark, S. (2019b).  Bayesian joint spike-and-slab graphical lasso.  In International Conference on Machine Learning, pages 3877–3885. PMLR.
  * 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.
  * Liu et al. (2009) Liu, H., Lafferty, J., and Wasserman, L. (2009).  The nonparanormal: Semiparametric estimation of high dimensional undirected graphs.  Journal of Machine Learning Research, 10, 2295–2328.
  * Mazumder and Hastie (2012) Mazumder, R. and Hastie, T. (2012).  The graphical lasso: New insights and alternatives.  Electronic Journal of Statistics, 6, 2125–2149.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the Lasso.  The Annals of Statistics, 34, 1436–1462.
  * Mitra et al. (2013) Mitra, R., Müller, P., Liang, S., Yue, L., and Ji, Y. (2013).  A Bayesian graphical model for chip-seq data on histone modifications.  Journal of the American Statistical Association, 108, 69–80.
  * Mohammadi and Wit (2015) Mohammadi, A. and Wit, E. C. (2015).  Bayesian structure learning in sparse Gaussian graphical models.  Bayesian Analysis, 10, 109–138.
  * Mohammadi et al. (2021) Mohammadi, R., Massam, H., and Letac, G. (2021).  Accelerating Bayesian structure learning in sparse Gaussian graphical models.  Journal of the American Statistical Association.  To appear.
  * Müller et al. (2004) Müller, P., Parmigiani, G., Robert, C., and Rousseau, J. (2004).  Optimal sample size for multiple testing: The case of gene expression microarrays.  Journal of the American Statistical Association, 99, 990–1001.
  * Neal (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 9, 249–265.
  * Painter et al. (2011) Painter, M. W., Davis, S., Hardy, R. R., Mathis, D., Benoist, C., Consortium, I. G. P., et al. (2011).  Transcriptomes of the B and T lineages compared by multiplatform microarray profiling.  The Journal of Immunology, 186, 3047–3057.
  * 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.
  * Peng et al. (2009) Peng, J., Wang, P., Zhou, N., and Zhu, J. (2009).  Partial correlation estimation by joint sparse regression models.  Journal of the American Statistical Association, 104, 735–746.
  * Pitt et al. (2006) Pitt, M., Chan, D., and Kohn, R. (2006).  Efficient Bayesian inference for Gaussian copula regression models.  Biometrika, 93, 537–554.
  * Polson and Scott (2010) Polson, N. G. and Scott, J. G. (2010).  Shrink globally, act locally: Sparse Bayesian regularization and prediction.  Bayesian Statistics, 9, 1–24.
  * Pourahmadi (2013) Pourahmadi, M. (2013).  High-dimensional covariance estimation.  John Wiley & Sons.
  * Rothman et al. (2008) Rothman, A. J., Bickel, P. J., Levina, E., Zhu, J., et al. (2008).  Sparse permutation invariant covariance estimation.  Electronic Journal of Statistics, 2, 494–515.
  * Roverato (2002) Roverato, A. (2002).  Hyper inverse Wishart distribution for non-decomposable graphs and its application to Bayesian inference for Gaussian graphical models.  Scandinavian Journal of Statistics, 29, 391–411.
  * 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 preprint arXiv:1612.02875.
  * Scott and Berger (2010) Scott, J. G. and Berger, J. O. (2010).  Bayes and empirical-Bayes multiplicity adjustment in the variable-selection problem.  The Annals of Statistics, 38, 2587–2619.
  * Shalizi (2009) Shalizi, C. R. (2009).  Dynamics of Bayesian updating with dependent data and misspecified models.  Electronic Journal of Statistics, 3, 1039–1074.
  * Wang (2012) Wang, H. (2012).  Bayesian graphical lasso models and efficient posterior computation.  Bayesian Analysis, 7, 867–886.
  * West (1992) West, M. (1992).  Hyperparameter estimation in Dirichlet process mixture models.  Duke University ISDS Discussion Paper# 92-A03.
  * Witten et al. (2011) Witten, D. M., Friedman, J. H., and Simon, N. (2011).  New insights and faster computations for the graphical lasso.  Journal of Computational and Graphical Statistics, 20, 892–900.
  * Wolffe (2001) Wolffe, A. (2001).  Histone genes.  In Encyclopedia of Genetics, pages 948–952. Academic Press, New York.
  * Yoshida and West (2010) Yoshida, R. and West, M. (2010).  Bayesian learning in sparse graphical factor models via variational mean-field annealing.  Journal of Machine Learning Research, 11, 1771–1798.
  * Yuan and Lin (2007) Yuan, M. and Lin, Y. (2007).  Model selection and estimation in the Gaussian graphical model.  Biometrika, 94, 19–35.
  * Zhang and Zou (2014) Zhang, T. and Zou, H. (2014).  Sparse precision matrix estimation via lasso penalized D-trace loss.  Biometrika, 101, 103–120.
  * Zhu et al. (2014) Zhu, H., Khondker, Z., Lu, Z., and Ibrahim, J. G. (2014).  Bayesian generalized low rank regression models for neuroimaging phenotypes and genetic markers.  Journal of the American Statistical Association, 109, 977–990.

Supplementary Materials for
Bayesian Scalable Precision Factor
Analysis for Massive Sparse
Gaussian Graphical Models

Noirrit Kiran Chandraa<EMAIL_ADDRESS>
Peter Müllerb,c<EMAIL_ADDRESS>
Abhra Sarkarb<EMAIL_ADDRESS>
aDepartment of Mathematical Sciences,

The University of Texas at Dallas,

800 W. Campbell Rd, Richardson, TX 75080-3021, USA

bDepartment of Statistics and Data Sciences,

The University of Texas at Austin,

2317 Speedway D9800, Austin, TX 78712-1823, USA

cDepartment of Mathematics,

The University of Texas at Austin,

2515 Speedway, PMA 8.100, Austin, TX 78712-1823, USA

Supplementary materials discuss a general strategy for posterior computation
under a broad class of generic priors, proofs of the theoretical results, some
additional figures, and an application to a NASDAQ-100 stock price dataset.

##  S.1 General Strategy for Posterior Computation

In this section, we sketch out a general strategy for posterior inference for
the precision factor model not limited to the prior specifications in Section
2.3.

Sampling model:

The data generating mechanism, we recall, is

$\displaystyle{\mathbf{y}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d}({\mathbf{0}},\mbox{\boldmath$\Omega$}^{-1})~{}~{}~{}\text{
with
}~{}~{}~{}\mbox{\boldmath$\Omega$}=\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\Lambda$}^{\rm
T}+\mbox{\boldmath$\Delta$},$

where $\Lambda$ is a $d\times q$ dimensional matrix with $q\leq d$ and
$\mbox{\boldmath$\Delta$}=\mathrm{diag}(\delta_{1}^{2},\dots,\delta_{d}^{2})$.

Generic priors:

We consider the following broad classes of hierarchical priors on $\Lambda$
and $\Delta$ with $\xi$ and $\zeta$ being the respective associated
hyperparameters

$\displaystyle\mathrm{vec}(\mbox{\boldmath$\Lambda$})\lvert\mbox{\boldmath$\xi$}\sim\mathrm{N}_{dq}\\{{\mathbf{0}},g(\mbox{\boldmath$\xi$})\\},\quad\mbox{\boldmath$\xi$}\sim\Pi_{\mbox{\boldmath$\xi$}},\quad\quad\mbox{\boldmath$\Delta$}\lvert\mbox{\boldmath$\zeta$}\sim\Pi_{\mbox{\boldmath$\Delta$}\lvert\mbox{\boldmath$\zeta$}}(\mbox{\boldmath$\Delta$}\lvert\mbox{\boldmath$\zeta$}),\quad\mbox{\boldmath$\zeta$}\sim\Pi_{\mbox{\boldmath$\zeta$}}.$

Varying choices of
$g(\mbox{\boldmath$\xi$}),\Pi_{\mbox{\boldmath$\xi$}},\Pi_{\mbox{\boldmath$\Delta$}\lvert\mbox{\boldmath$\zeta$}},\Pi_{\mbox{\boldmath$\zeta$}}$
then produce a broad class of priors for the model parameters.

The Gibbs sampler then iterates through the following steps.

Step 1. Sample the latent variables:

Generate
${\mathbf{u}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q}({\mathbf{0}},{\mathbf{P}})$
with ${\mathbf{P}}=({\mathbf{I}}_{q}+\mbox{\boldmath$\Lambda$}^{\rm
T}\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$})$ independently from
${\mathbf{y}}_{1:n}$ and let
${\mathbf{v}}_{1:n}={\mathbf{y}}_{1:n}+\mbox{\boldmath$\Delta$}^{-1}\mbox{\boldmath$\Lambda$}{\mathbf{P}}^{-1}{\mathbf{u}}_{1:n}$.

Step 2. Sample the rows of $\Lambda$:

We have
${\mathbf{u}}_{i}=\sum_{r=1}^{d}\mbox{\boldmath$\lambda$}_{r}v_{r,i}+\mbox{\boldmath$\varepsilon$}_{i}$,
where $\mbox{\boldmath$\lambda$}_{r}=(\lambda_{r,1},\dots,\lambda_{r,q})$ is
the $r^{th}$ row of $\Lambda$ and
${\mathbf{v}}_{i}=(v_{1,i},\dots,v_{d,i})^{\rm T}$. Define
${\mathbf{u}}_{i}^{(j)}={\mathbf{u}}_{i}-\sum_{r\neq
j}\mbox{\boldmath$\lambda$}_{r}v_{r,i}$. Then
${\mathbf{u}}_{i}^{(j)}=\mbox{\boldmath$\lambda$}_{j}v_{j,i}+\mbox{\boldmath$\varepsilon$}_{i}$.
Conditioned on ${\mathbf{u}}_{i}^{(j)}$, ${\mathbf{v}}_{i}$ and the associated
hyper-parameters, $\mbox{\boldmath$\lambda$}_{j}$’s can be updated
sequentially $j=1,\dots,d$ from the following posterior distribution

$\mbox{\boldmath$\lambda$}_{j}\sim\mathrm{N}_{q}\\{({\mathbf{D}}_{j}^{-1}+\norm{{\mathbf{v}}^{(j)}}^{2}{\mathbf{I}}_{q})^{-1}{\mathbf{w}}_{j},({\mathbf{D}}_{j}^{-1}+\norm{{\mathbf{v}}^{(j)}}^{2}{\mathbf{I}}_{q})^{-1}\\},$

where ${\mathbf{D}}_{j}$ is the prior covariance matrix of
$\mbox{\boldmath$\lambda$}_{j}$,
${\mathbf{v}}^{(j)}=(v_{j,1},\dots,v_{j,n})^{\rm T}$ and
${\mathbf{w}}_{j}=\sum_{i=1}^{n}v_{j,i}{\mathbf{u}}_{i}^{(j)}$.

Step 3. Sample $\Delta$:

Sample $\Delta$ from
$\Pi_{\mbox{\boldmath$\Delta$}\lvert{\mathbf{v}}^{(1:d)},\mbox{\boldmath$\zeta$}}$.

Step 4. Sample the hyperparameters of $\Lambda$:

Sample $\xi$ from
$\mbox{\boldmath$\xi$}\sim\Pi_{\mbox{\boldmath$\xi$}\lvert\mbox{\boldmath$\Lambda$}}$.

Step 5. Sample the hyperparameters of $\Delta$:

Sample $\zeta$ from
$\mbox{\boldmath$\zeta$}\sim\Pi_{\mbox{\boldmath$\zeta$}\lvert\mbox{\boldmath$\Delta$}}$.

For non-conjugate priors, appropriate Metropolis-within-Gibbs type MCMC
algorithms can be employed in Steps 2, 3, 4 or 5.

##  S.2 Proofs of Theoretical Results

Notations: In what follows, for an operation ‘$*$’, $\overline{a*b}$ is
sometimes used for $(a*b)$. Also, 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}}}$ is the Euclidean norm of a vector ${\mathbf{x}}$. We
denote the Kullback-Leibler (KL) divergence between two mean zero Gaussian
distributions with precision matrices $\Omega$ and
$\mbox{\boldmath$\Omega$}^{\prime}$ by
$\mathbb{KL}\left(\mbox{\boldmath$\Omega$}\parallel\mbox{\boldmath$\Omega$}^{\prime}\right)$.
We borrow the following result from pati2014 (Lemma 1.1. from the supplement).
For brevity of notations, we reuse $C$, $\widetilde{C}$, $C^{\prime\prime}$,
etc. in the proofs to denote constants and their values may not be the same
throughout the same proof. Nevertheless, we were careful to make sure that
these quantities are indeed constants.

######  Lemma 1.

For any two matrices ${\mathbf{A}}$ and ${\mathbf{B}}$,

  1. (i)

$s_{\min}\left({\mathbf{A}}\right)\norm{{\mathbf{B}}}_{F}\leq\norm{{\mathbf{A}}{\mathbf{B}}}_{F}\leq\norm{{\mathbf{A}}}_{2}\norm{{\mathbf{B}}}_{F}$.

  2. (ii)

$s_{\min}\left({\mathbf{A}}\right)\norm{{\mathbf{B}}}_{2}\leq\norm{{\mathbf{A}}{\mathbf{B}}}_{2}\leq\norm{{\mathbf{A}}}_{2}\norm{{\mathbf{B}}}_{2}$.

  3. (iii)

$s_{\min}\left({\mathbf{A}}\right)s_{\min}\left({\mathbf{B}}\right)\leq
s_{\min}\left({\mathbf{A}}{\mathbf{B}}\right)\leq\norm{{\mathbf{A}}}_{2}s_{\min}\left({\mathbf{B}}\right)$.

###### Theorem 1.

Let $\mathcal{P}_{n}=\mathcal{P}_{1,n}\cup\mathcal{P}_{2,n}$ be an arbitrary
partition. Note that,

$\displaystyle\Pi_{n}\left(\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}>M_{n}\epsilon_{n}\lvert{\mathbf{y}}_{1:n}\right)\leq\Pi_{n}\left(\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}>M_{n}\epsilon_{n}\lvert{\mathbf{y}}_{1:n}\right)+\Pi_{n}\left(\mathcal{P}_{2,n}\lvert{\mathbf{y}}_{1:n}\right).$

To prove the theorem, we show that, as $n\to\infty$, the posterior
distribution on $\mathcal{P}_{1,n}$ concentrates around
$\mbox{\boldmath$\Omega$}_{0n}$ while the remaining mass assigned to
$\mathcal{P}_{2,n}$ diminishes to zero. We adopt a variation of Theorem 8.22
by ghosal_book on posterior contraction rates. Define
$B_{n,0}(\mbox{\boldmath$\Omega$}_{0n},\epsilon)=\left\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{n}:\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)\leq
n\epsilon^{2}\right\\}$. We now formally state the result tailored towards our
application.

######  Theorem 2.

Let $\mathcal{P}_{n}$ be a class of distributions prametrized by
$\mbox{\boldmath$\Omega$}_{n}$, $\mathbb{P}_{\mbox{\boldmath$\Omega$}_{n}}$
the corresponding distribution and $\mbox{\boldmath$\Omega$}_{0n}$ be the true
data-generating value of the parameter. Let $e_{n}$ be a metric on
$\mathcal{P}_{n}$ and $\mathcal{P}_{1,n}\subset\mathcal{P}_{n}$. For constants
$\tau_{n}$ with $n\tau_{n}^{2}\geq 1$ and every sufficiently large
$j\in\mathbb{N}$, assume the following conditions hold.

  1. (i)

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

  2. (ii)

Define the set
$G_{j,n}=\left\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}<e_{n}\left(\mbox{\boldmath$\Omega$}_{n},\mbox{\boldmath$\Omega$}_{0n}\right)\leq
2j\tau_{n}\right\\}$. There exists tests $\phi_{n}$ such that, for some $K>0$,

$\displaystyle\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{0n}}\phi_{n}=0;\qquad\sup_{\mbox{\boldmath$\Omega$}_{n}\in
G_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{n}}(1-\phi_{n})\leq\exp(-Knj^{2}\tau_{n}^{2}).$

Then,
$\Pi_{n}\left\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:e_{n}\left(\mbox{\boldmath$\Omega$}_{n},\mbox{\boldmath$\Omega$}_{0n}\right)>M_{n}\tau_{n}\lvert{\mathbf{y}}_{1:n}\right\\}\to
0$ in $\mathbb{P}_{\mbox{\boldmath$\Omega$}_{0n}}$-probability for any
$M_{n}\to\infty$.

We define
$\mathcal{P}_{1,n}=\\{\mbox{\boldmath$\Omega$}_{n}:\norm{\mbox{\boldmath$\Omega$}_{n}}_{2}<C^{\prime\prime}(q_{0n}s_{n})^{4}\log(d_{n}q_{n})\\}$,
$\mathcal{P}_{2,n}=\mathcal{P}_{n}\backslash\mathcal{P}_{1,n}$,
$e_{n}\left(\mbox{\boldmath$\Omega$}_{n},\mbox{\boldmath$\Omega$}_{0n}\right)=\frac{\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}}{C^{\prime\prime}(q_{0n}s_{n})^{4}\log(d_{n}q_{n})}$
and $n\tau_{n}^{2}=q_{0n}s_{n}\log(d_{n}q_{n})$ for some large enough constant
$C^{\prime\prime}>0$. We verify (i) in Lemma 3 and show the existence of a
sequence of test functions satisfying (ii) in Lemma 4. Hence,

$\displaystyle\Pi_{n}\left\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:e_{n}\left(\mbox{\boldmath$\Omega$}_{n},\mbox{\boldmath$\Omega$}_{0n}\right)>M_{n}\tau_{n}\lvert{\mathbf{y}}_{1:n}\right\\}\to
0~{}\text{in}~{}\mathbb{P}_{\mbox{\boldmath$\Omega$}_{0n}}\text{-probability
for every}~{}M_{n}\to\infty.$

Subsequently applying the dominated convergence theorem (DCT), we get
$\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{0n}}\Pi_{n}(\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}>M_{n}\tau_{n}\lvert{\mathbf{y}}_{1:n})=0$.
To conclude the proof, we show that the remaining mass assigned to
$\mathcal{P}_{2,n}$ goes to 0 in the following theorem.

######  Theorem 3.

$\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{0n}}\Pi_{n}\left(\mathcal{P}_{2,n}\lvert{\mathbf{y}}_{1:n}\right)=0$.

∎

###### Theorem 2.

The theorem closely resembles Theorem 8.22 from ghosal_book. In the discussion
following equation (8.22) in the book, the authors noted that for the iid
case, simpler theorems are obtained by using an absolute lower bound on the
prior mass and by replacing the local entropy by the global entropy. In
particular, (i) implies Theorem 8.19(i) in the book. Similarly, (ii) is the
same condition in Theorem 8.20 and hence the proof. ∎

######  Lemma 2.

(Kullback-Leibler upper bound) Let $\mbox{\boldmath$\Omega$}_{n}$ and
$\mbox{\boldmath$\Omega$}_{0n}$ be $d_{n}\times d_{n}$ order positive definite
matrices. Then,

$\displaystyle
2\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)$
$\displaystyle=-\log\absolutevalue{\mbox{\boldmath$\Omega$}_{0n}^{-1}\mbox{\boldmath$\Omega$}_{n}}+\mathrm{trace}(\mbox{\boldmath$\Omega$}_{0n}^{-1}\mbox{\boldmath$\Omega$}_{n}-{\mathbf{I}}_{d_{n}})\leq\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}^{2}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}/2{\delta}_{\min}^{6}.$

###### Proof.

Let
${\mathbf{H}}=\mbox{\boldmath$\Omega$}_{0n}^{-\nicefrac{{1}}{{2}}}\mbox{\boldmath$\Omega$}_{n}\mbox{\boldmath$\Omega$}_{0n}^{-\nicefrac{{1}}{{2}}}$.
Letting $\psi_{1},\dots,\psi_{d_{n}}$ be the eigenvalues of ${\mathbf{H}}$, we
note that

$\displaystyle
2\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)=\sum_{j=1}^{d_{n}}\left\\{(\psi_{j}-1)-\log\psi_{j}\right\\}.$
(S.1)

Consider the function $h_{\beta}(x)=\log x-(x-1)+\beta(x-1)^{2}/2$ on
$(0,\infty)$ for $\beta>1$. Note that $h_{\beta}(x)\geq 0$ for all $x\geq
1/\beta$. From (C3) and using Lemma 1, $\psi_{j}\geq
s_{\min}\left({\mathbf{H}}\right)\geq
s_{\min}\left(\mbox{\boldmath$\Omega$}_{n}\right)s_{\min}\left(\mbox{\boldmath$\Omega$}_{0n}^{-1}\right)={\delta}_{\min}^{2}/\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}$.
Noting that ${\delta}_{\min}^{2}/\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}<1$,
we set $\beta=\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}/{\delta}_{\min}^{2}$.
Therefore,

$\displaystyle\sum_{j=1}^{d_{n}}\left\\{\log\psi_{j}-(\psi_{j}-1)\right\\}\geq-\frac{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}{2{\delta}_{\min}^{2}}\sum_{j=1}^{d_{n}}(\psi_{j}-1)^{2}=-\frac{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}{2{\delta}_{\min}^{2}}\norm{{\mathbf{H}}-{\mathbf{I}}_{d_{n}}}_{F}^{2}.$
(S.2)

Again, using Lemma 1 and (C3), we have

$\displaystyle\norm{{\mathbf{H}}-{\mathbf{I}}_{d_{n}}}_{F}\leq\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}\norm{\mbox{\boldmath$\Omega$}_{0n}^{-1}}_{2}\leq\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}/{\delta}_{\min}^{2}.$
(S.3)

Combing (S.1)-(S.3), we conclude the lemma. ∎

######  Lemma 3.

Define the set
$B_{n,0}(\mbox{\boldmath$\Omega$}_{0n},\tau)=\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{n}:\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)\leq
n\tau^{2}\\}$. Then,
$\Pi_{n}\left\\{B_{n,0}(\mbox{\boldmath$\Omega$}_{0n},\tau_{n})\right\\}\geq
e^{-Cn\tau_{n}^{2}}$ for $n\tau_{n}^{2}=s_{n}q_{0n}\log(d_{n}q_{n})$ and some
absolute constant $C>0$.

###### Proof.

From Lemma 2, we note that

$\displaystyle\Pi_{n}\left\\{B_{n,0}(\mbox{\boldmath$\Omega$}_{0n},\tau_{n})\right\\}\geq\Pi_{n}\left(\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}^{2}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}/4{\delta}_{\min}^{6}\leq
n\tau_{n}^{2}\right)$ $\displaystyle\geq$
$\displaystyle\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}-\mbox{\boldmath$\Lambda$}_{0n}\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}}_{F}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}^{\nicefrac{{1}}{{2}}}}\right)\Pi_{n}\left(\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}^{\nicefrac{{1}}{{2}}}}\right).$
(S.4)

We handle the $\Lambda$ and $\Delta$ parts in (S.4) separately and conclude
the proof by showing that each of the terms individually exceeds
$e^{-Cn\tau_{n}^{2}}$ for some constant $C>0$.

The $\Lambda$ part in (S.4):  We define
$\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}=[\mbox{\boldmath$\Lambda$}_{0n}\quad{\mathbf{0}}^{d_{n}\times\overline{q_{n}-q_{0n}}}]$,
after augmenting $q_{n}-q_{0n}$ null columns to
$\mbox{\boldmath$\Lambda$}_{0n}$. Also
$\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}^{\rm
T}=\mbox{\boldmath$\Lambda$}_{0n}\mbox{\boldmath$\Lambda$}_{0n}^{\rm T}$.
Invoking (C2) and Lemma 1, we have that
$\norm{\mbox{\boldmath$\Lambda$}_{n}-\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}}_{F}<\epsilon$
implies that
$\norm{\mbox{\boldmath$\Lambda$}_{n}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}-\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}^{\rm
T}}_{F}\leq C\epsilon^{2}$ for some $C>0$. Note that, from (C2) and (C3),
$\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}=O(1)$. Let ${\mathbf{A}}^{S}$ denote
the vector of elements of ${\mathbf{A}}$ corresponding to an index set $S$,
and $S_{0}$ denote the set of nonzero elements in
$\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}$ . Then, for some absolute
constant $\widetilde{C}$, we have

$\displaystyle\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}-\mbox{\boldmath$\Lambda$}_{0n}\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}}_{F}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}^{\nicefrac{{1}}{{2}}}}\right)\geq\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}-\widetilde{\mbox{\boldmath$\Lambda$}}_{0n}}_{F}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}\right)$
$\displaystyle\geq\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}}-\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}^{C}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)$
(S.5)

The first term in (S.5) ensures that the Bayesian model assigns requisite
prior mass around the signals whereas the second term ensures that appropriate
shrinkage is envisaged through the prior. We show that each of these products
exceeds $e^{-Cn\tau_{n}^{2}}$.

The shrinkage part in (S.5):

From [equation (10)]bhattacharya2015dirichlet, we get that
$\mbox{\boldmath$\Lambda$}_{n}=\frac{2}{b_{n}}((\lambda_{n,j,h}^{*}))$ where
$\lambda_{n,j,h}^{*}$’s are marginally iid random variables with a priori with
pdf
$\Pi(\lambda_{n,j,h})=\absolutevalue{\lambda_{n,j,h}}^{(a_{n}-1)/2}K_{1-a_{n}}(\sqrt{2\absolutevalue{\lambda_{n,j,h}}})/\\{2^{(1+a_{n})/2}\Gamma(a_{n})\\}$
where
$K_{\nu}(x)=\frac{\Gamma(\nu+1/2)(2x)^{\nu}}{\sqrt{\pi}}\int_{0}^{\infty}\frac{\cos
t}{(t^{2}+x^{2})^{\nu+1/2}}\mathrm{d}t$ is the modified Bessel function of the
second kind. Using Lemma 3.2 of the same paper and letting
$H_{n}=b_{n}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}$, we have

$\displaystyle\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}^{C}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)\geq\left\\{\Pr(\absolutevalue{\lambda_{n,j,h}^{*}}\leq\frac{CH_{n}}{\absolutevalue{S_{0}^{C}}})\right\\}^{\absolutevalue{S_{0}^{C}}}\geq\left\\{1-\frac{C}{\Gamma(a_{n})}\log\frac{\absolutevalue{S_{0}^{C}}}{H_{n}}\right\\}^{\absolutevalue{S_{0}^{C}}}.$

Now, $\absolutevalue{S_{0}^{C}}=d_{n}q_{n}-s_{n}q_{0n}$ and
$\Gamma(a_{n})=\Gamma(1+a_{n})/a_{n}$. Since $a_{n}=1/(d_{n}q_{n})$ and
$H_{n}$ is an increasing sequence,
$\frac{1}{\Gamma(a_{n})}\log\frac{\absolutevalue{S_{0}^{C}}}{H_{n}}=\frac{1}{d_{n}q_{n}\Gamma(1+a_{n})}\log\frac{d_{n}q_{n}-s_{n}q_{0n}}{H_{n}}\leq
1$ and hence

$\displaystyle\begin{split}\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}^{C}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)&\geq\left\\{1-\frac{C}{\Gamma(a_{n})}\log\frac{\absolutevalue{S_{0}^{C}}}{H_{n}}\right\\}^{\absolutevalue{S_{0}^{C}}}\geq\left\\{1-\frac{C\log(d_{n}q_{n})}{d_{n}q_{n}\Gamma(1+a_{n})}\right\\}^{d_{n}q_{n}}\\\
&\geq
e^{-C\log(d_{n}q_{n})}\left\\{1-\frac{C^{2}\log^{2}(d_{n}q_{n})}{d_{n}q_{n}}\right\\}\geq
e^{-C^{\prime}q_{0n}s_{n}\log(d_{n}q_{n})}.\end{split}$

The second last inequality in the previous equation follows since, for any
$n\geq 1$ and $\absolutevalue{x}\leq n$, we have
$\left(1+\frac{x}{n}\right)^{n}\geq e^{x}\left(1-\frac{x^{2}}{n}\right)$.

The signal part in (S.5):

Let us define $v_{q}(r)$ to be the $q$-dimensional Euclidean ball of radius
$r$ centered at zero and $\absolutevalue{v_{q}(r)}$ denotes its volume. For
the sake of brevity, denote $v_{q}=\absolutevalue{v_{q}(1)}$, so that
$\absolutevalue{v_{q}(r)}=r^{q}v_{q}$. Letting
$\mbox{\boldmath$\Lambda$}_{n}^{*}=((\lambda^{*}_{n,j,h}))^{d_{n}\times
q_{n}}$ ($\lambda^{*}_{n,j,h}$’s defined earlier) and
$t_{n}=b_{n}\left(\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2+\norm{\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}\right)$,
we have

$\displaystyle\begin{split}\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}}-\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)&=\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{*S_{0}}-b_{n}\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}\leq\widetilde{C}b_{n}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)\\\
&\geq\absolutevalue{v_{\absolutevalue{S_{0}}}(t_{n})}\inf_{v_{\absolutevalue{S_{0}}}(t_{n})}\Pi_{n}(\mbox{\boldmath$\Lambda$}_{n}^{*S_{0}}).\end{split}$
(S.6)

Note that, $\absolutevalue{S_{0}}=s_{n}q_{0n}$ and
$\norm{\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}=\norm{\mbox{\boldmath$\Lambda$}_{0n}}_{F}$
and using [Lemma 5.3]castillo2012needles, $v_{q}\asymp(2\pi
e)^{q/2}q^{-\overline{q+1}/2}$. Note also that for $x>0$,
$\frac{\log(1+x)}{x}\leq\frac{1}{\sqrt{1+x}}$ which implies that $x^{x}\leq
e^{x^{3/2}}$. Hence,

$\displaystyle\absolutevalue{v_{\absolutevalue{S_{0}}}(t_{n})}=t_{n}^{s_{n}q_{0n}}v_{s_{n}q_{0n}}\geq\exp\\{s_{n}q_{0n}\log
t_{n}-C^{\prime}(q_{0n}s_{n})^{3/2}/2\\}\geq\exp\\{-Cs_{n}q_{0n}\log(d_{n}q_{n})\\}.$

The last inequality follows from (C1) and the prior specifications. From
[Lemma 3.1]bhattacharya2015dirichlet we have

$\displaystyle\inf_{v_{\absolutevalue{S_{0}}}(t_{n})}\Pi_{n}(\mbox{\boldmath$\Lambda$}_{n}^{*S_{0}})=\exp[-C\left\\{s_{n}q_{0n}\log(1/a_{n})+(s_{n}q_{0n})^{3/4}\sqrt{t_{n}}\right\\}]$
$\displaystyle=\exp\left[-C\left\\{s_{n}q_{0n}\log(d_{n}q_{n})+(s_{n}q_{0n})^{3/4}b_{n}^{\nicefrac{{1}}{{2}}}\left(\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2+\norm{\mbox{\boldmath$\Lambda$}_{0n}}_{F}\right)^{\nicefrac{{1}}{{2}}}\right\\}\right]$
$\displaystyle\geq\exp(-Cs_{n}q_{0n}\log\overline{d_{n}q_{n}}),$

since, $b_{n}=\log^{3/2}(d_{n}q_{n})$,
$n\tau_{n}^{2}=s_{n}q_{0n}\log\overline{d_{n}q_{n}}$ and
$\norm{\mbox{\boldmath$\Lambda$}_{0n}}_{F}^{2}=O(q_{0n})$ from (C2). Hence
from (S.6),
$\Pi_{n}\left(\norm{\mbox{\boldmath$\Lambda$}_{n}^{S_{0}}-\mbox{\boldmath$\Lambda$}_{0n}^{S_{0}}}\leq\widetilde{C}n^{1/4}\tau_{n}^{\nicefrac{{1}}{{2}}}/2\right)\geq\exp(-Cs_{n}q_{0n}\log\overline{d_{n}q_{n}})$.

The $\Delta$ part in (S.4): Note that,
$\Pi_{n}\left(\norm{\mbox{\boldmath$\Delta$}_{n}-\mbox{\boldmath$\Delta$}_{0n}}_{F}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}^{\nicefrac{{1}}{{2}}}}\right)\geq\Pi_{n}\left(\absolutevalue{\delta_{0n}^{2}-\delta_{n}^{2}}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\sqrt{d_{n}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}}\right)$.
Now, $\Pi_{n}(\absolutevalue{\delta_{0n}^{2}-\delta_{n}^{2}}\leq x)\geq
e^{-\delta_{0n}^{2}}(1-e^{-2x})$. Using $1-e^{-2x}\geq x$ for
$x\in(0,\nicefrac{{1}}{{2}})$ and since
$\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\sqrt{d_{n}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}}\to
0$ as $n\to\infty$, we have for some $C^{\prime}>0$

$\displaystyle\Pi_{n}\left(\absolutevalue{\delta_{0n}^{2}-\delta_{n}^{2}}\leq\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\sqrt{d_{n}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}}\right)\geq\exp(-\delta_{0}^{2}+\log\frac{\sqrt{n}\tau_{n}{\delta}_{\min}^{3}}{\sqrt{d_{n}\norm{\mbox{\boldmath$\Omega$}_{0n}}_{2}}})\geq\exp(-C^{\prime}s_{n}q_{0n}\log\overline{d_{n}q_{n}}).$

∎

######  Corollary 1.

For arbitrary $\widetilde{\tau}_{n}>n^{-\nicefrac{{1}}{{2}}}$, we have
$\Pi_{n}\left\\{B_{n,0}(\mbox{\boldmath$\Omega$}_{0n},\widetilde{\tau}_{n})\right\\}\geq
e^{-Cq_{0n}s_{n}\log(d_{n}q_{n})}$ for some constant $C>0$.

Following the proof of Lemma 3, it can be seen that the above corollary holds.

######  Lemma 4.

Let $\mbox{\boldmath$\Omega$}_{0n}\in\mathcal{P}_{0n}$ and the set
$G_{j,n}=\left\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}:j\tau_{n}<\frac{1}{\zeta_{n}}\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}\leq
2j\tau_{n}\right\\}$ denote an annulus of inner radius $j\tau_{n}\zeta_{n}$
and outer radius $(j+1)\tau_{n}\zeta_{n}$, where
$\zeta_{n}=C^{\prime\prime}(q_{0n}s_{n})^{4}\log(d_{n}q_{n})$, in operator
norm around $\mbox{\boldmath$\Omega$}_{0n}$ for some integer $j>1$. Based on
iid samples ${\mathbf{y}}_{1:n}$ from
$\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Omega$}_{0n}^{-1})$,
consider the following hypothesis testing problem

$\displaystyle
H_{0}:\mbox{\boldmath$\Omega$}_{n}=\mbox{\boldmath$\Omega$}_{0n}\text{ versus
}H_{1}:\mbox{\boldmath$\Omega$}_{n}\in G_{j,n}.$

We simulate
${\mathbf{u}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q_{0n}}({\mathbf{0}},{\mathbf{P}}_{0n})$,
with
${\mathbf{P}}_{0n}=({\mathbf{I}}_{q_{0n}}+\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\mbox{\boldmath$\Lambda$}_{0n})$
independently from ${\mathbf{y}}_{1:n}$ and define
${\mathbf{v}}_{i}={\mathbf{y}}_{i}+\mbox{\boldmath$\Delta$}_{0n}^{-1}\mbox{\boldmath$\Lambda$}_{0n}{\mathbf{P}}_{0n}^{-1}{\mathbf{u}}_{i}$.
Letting ${\mathbf{V}}^{\rm T}=({\mathbf{v}}_{1},\dots,{\mathbf{v}}_{n})$, we
define the following test function
$\phi_{n}=\mathbbm{1}\left\\{\norm{\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\left(\frac{1}{n}{\mathbf{V}}^{\rm
T}{\mathbf{V}}-\mbox{\boldmath$\Delta$}_{0n}^{-1}\right)\mbox{\boldmath$\Lambda$}_{0n}}_{2}>\tau_{n}\right\\}$.
Then

$\displaystyle\lim_{n\to\infty}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{0n}}\phi_{n}=0;\qquad\sup_{\mbox{\boldmath$\Omega$}_{n}\in
G_{j,n}}\mathbb{E}_{\mbox{\boldmath$\Omega$}_{n}}(1-\phi_{n})\leq\exp(-Knj^{2}\tau_{n}^{2}),$

for some absolute constant $K>0$.

###### Proof.

Type-I error:

Under $H_{0}$, we have
${\mathbf{v}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Delta$}_{0n}^{-1})$.
Letting
$\mbox{\boldmath$\Xi$}_{0n}=\hbox{cov}\left({\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}{\mathbf{v}}_{i}}\right)=\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\mbox{\boldmath$\Lambda$}_{0n}$, we have,
$\phi_{n}\leq\mathbbm{1}\left(\norm{\mbox{\boldmath$\Xi$}_{0n}}_{2}\norm{\frac{1}{n}\sum_{i=1}^{n}{\mathbf{z}}_{i}{\mathbf{z}}_{i}^{\rm
T}-{\mathbf{I}}_{q_{0n}}}_{2}>\tau_{n}\right)$, where
$\norm{\mbox{\boldmath$\Xi$}_{0n}}_{2}=O(1)$ from (C2). From [Corollary
5.50]vershynin_book, $\mathbb{E}_{H_{0}}\phi_{n}\leq
2\exp(-\widetilde{C}q_{0n}t_{n}^{2})$ for a universal constant
$\widetilde{C}>0$ and any sequence $t_{n}$ satisfying $q_{0n}t_{n}^{2}\leq
n\tau_{n}^{2}$. Since from (C1), $q_{0n}\to\infty$ and
$n\tau_{n}^{2}/q_{0n}\to\infty$ as $n\to\infty$, $t_{n}$ can be constructed
such that $t_{n}\succsim 1$ and hence
$\lim_{n\to\infty}\mathbb{E}_{H_{0}}\phi_{n}=0$.

Type-II error:

If
${\mathbf{y}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Omega$}_{n})$
with
$\mbox{\boldmath$\Omega$}_{n}=\mbox{\boldmath$\Lambda$}_{n}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}+\mbox{\boldmath$\Delta$}_{n}$, then
${\mathbf{v}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Delta$}_{n}^{-1}+\mbox{\boldmath$\nabla$}_{n})$
where
$\mbox{\boldmath$\nabla$}_{n}=\mbox{\boldmath$\Delta$}_{0n}^{-1}\mbox{\boldmath$\Lambda$}_{0n}({\mathbf{I}}_{q_{0n}}+\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}\mbox{\boldmath$\Lambda$}_{0n})^{-1}\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Delta$}_{0n}^{-1}-\mbox{\boldmath$\Delta$}_{n}^{-1}\mbox{\boldmath$\Lambda$}_{n}({\mathbf{I}}_{q_{n}}+\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}\mbox{\boldmath$\Delta$}_{n}^{-1}\mbox{\boldmath$\Lambda$}_{n})^{-1}\mbox{\boldmath$\Lambda$}_{n}^{\rm
T}\mbox{\boldmath$\Delta$}_{n}^{-1}$. Notably,
$\mbox{\boldmath$\Delta$}_{n}^{-1}-\mbox{\boldmath$\Delta$}_{0n}^{-1}+\mbox{\boldmath$\nabla$}_{n}=\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}$.
Now,

$\displaystyle\hskip-25.0pt1-\phi_{n}=\mathbbm{1}\left\\{\norm{\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\left(\frac{1}{n}{\mathbf{V}}^{\rm
T}{\mathbf{V}}-\mbox{\boldmath$\Delta$}_{0n}^{-1}\right)\mbox{\boldmath$\Lambda$}_{0n}}_{2}\leq\tau_{n}\right\\}$
(S.7) $\displaystyle=$
$\displaystyle\mathbbm{1}\left\\{\norm{\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\left\\{\left(\frac{1}{n}{\mathbf{V}}^{\rm
T}{\mathbf{V}}-\mbox{\boldmath$\Delta$}_{n}^{-1}-\mbox{\boldmath$\nabla$}_{n}\right)+(\mbox{\boldmath$\Delta$}_{n}^{-1}-\mbox{\boldmath$\Delta$}_{0n}^{-1}+\mbox{\boldmath$\nabla$}_{n})\right\\}\mbox{\boldmath$\Lambda$}_{0n}}_{2}\leq\tau_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\left(\frac{1}{n}{\mathbf{V}}^{\rm
T}{\mathbf{V}}-\mbox{\boldmath$\Delta$}_{n}^{-1}-\mbox{\boldmath$\nabla$}_{n}\right)\mbox{\boldmath$\Lambda$}_{0n}}_{2}>\norm{\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}(\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n})\mbox{\boldmath$\Lambda$}_{0n}}_{2}-\tau_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\left\\{\frac{1}{n}\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}{\mathbf{V}}^{\rm
T}{\mathbf{V}}\mbox{\boldmath$\Lambda$}_{0n}-\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}(\mbox{\boldmath$\Delta$}_{n}^{-1}+\mbox{\boldmath$\nabla$}_{n})\mbox{\boldmath$\Lambda$}_{0n}\right\\}}_{2}>{s_{\min}\left(\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Lambda$}_{0n}\right)}\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}-\tau_{n}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbbm{1}\left\\{\norm{\frac{1}{n}\sum_{i=1}^{n}{\mathbf{z}}_{i}{\mathbf{z}}_{i}^{\rm
T}-{\mathbf{I}}_{q_{n}}}_{2}>\frac{s_{\min}\left(\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Lambda$}_{0n}\right)}{\norm{\mbox{\boldmath$\Delta$}_{n}^{-1}+\mbox{\boldmath$\nabla$}_{n}}_{2}\norm{\mbox{\boldmath$\Lambda$}_{0n}}_{2}^{2}}\left(\norm{\mbox{\boldmath$\Omega$}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}-\frac{\tau_{n}}{s_{\min}\left(\mbox{\boldmath$\Lambda$}_{0n}^{\rm
T}\mbox{\boldmath$\Lambda$}_{0n}\right)}\right)\right\\},~{}~{}~{}~{}~{}$

where
${\mathbf{z}}_{1:n}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathrm{N}_{q_{n}}({\mathbf{0}},{\mathbf{I}}_{q_{n}})$.
By construction of $\mathcal{P}_{1,n}$,
$\norm{\mbox{\boldmath$\Delta$}_{n}^{-1}+\mbox{\boldmath$\nabla$}_{n}}_{2}\leq\zeta_{n}$
for $\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{1,n}$. Hence, for
$\mbox{\boldmath$\Omega$}_{n}\in G_{j,n}$, the RHS of (S.7) is bounded below
by $j\tau_{n}$ for sufficiently large $j$. Using [Eqn 5.26]vershynin_book
again, we have for all $\mbox{\boldmath$\Omega$}_{n}\in G_{j,n}$,
$\mathbb{E}_{\mbox{\boldmath$\Omega$}_{n}}(1-\phi_{n})\leq
e^{-Knj^{2}\tau_{n}^{2}}$.

Hence the proof. ∎

###### Theorem 3.

For densities $p$ and $q$, define $\mathbb{V}\left(p\parallel
q\right)=\int\left\\{\log\frac{p}{q}-\mathbb{KL}\left(p\parallel
q\right)\right\\}^{2}\mathrm{d}p$. When $p$ and $q$ are mean zero multivariate
Gaussian densities with precision matrices $\Omega$ and
$\mbox{\boldmath$\Omega$}^{\prime}$, we simply denote
$\mathbb{V}\left(\mbox{\boldmath$\Omega$}\parallel\mbox{\boldmath$\Omega$}^{\prime}\right)$.
We define the set
$B_{n,2}(\mbox{\boldmath$\Omega$}_{0n},\epsilon)=\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{n}:\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)\leq
n\epsilon^{2},~{}\mathbb{V}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)\leq
n\epsilon^{2}\\}$.

From [Theorem 3.1]banerjee2015bayesian and the proof of Lemma 2, we have
$\mathbb{V}\left(\mbox{\boldmath$\Omega$}_{0n}\parallel\mbox{\boldmath$\Omega$}_{n}\right)=\frac{1}{2{\delta}_{\min}^{4}}\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}^{2}$
and therefore
$B_{n,2}(\mbox{\boldmath$\Omega$}_{0n},\epsilon)=\\{\mbox{\boldmath$\Omega$}_{n}\in\mathcal{P}_{n}:\norm{\mbox{\boldmath$\Omega$}_{0n}-\mbox{\boldmath$\Omega$}_{n}}_{F}^{2}\leq
C\epsilon\\}$ for some constant $C>0$.

Let $\\{t_{n}\\}_{n=1}^{\infty}$ be an increasing sequence of positive
numbers. Then, $\Pi_{n}(\norm{\mbox{\boldmath$\Lambda$}_{n}}_{2}\geq
t_{n})\leq\Pi_{n}(\norm{\mbox{\boldmath$\Lambda$}_{n}}_{F}\geq
t_{n})\leq\Pi_{n}(\norm{\mathrm{vec}(\mbox{\boldmath$\Lambda$}_{n}^{*})}_{\ell_{1}}\geq
b_{n}t_{n})\leq 2e^{-C\sqrt{b_{n}t_{n}}}$ for some constant $C>0$ where
$\norm{\cdot}_{\ell_{1}}$ is the $\ell_{1}$ norm. The last inequality follows
from [Lemma 7.4]pati2014. Hence, for
$\zeta_{n}=C^{\prime\prime}(q_{0n}s_{n})^{4}\log(d_{n}q_{n})$,

$\displaystyle\Pi_{n}(\mathcal{P}_{2,n})\leq\Pi_{n}\left\\{\delta_{\max}^{2}+\norm{\mbox{\boldmath$\Lambda$}_{n}}_{2}^{2}\geq\zeta_{n}\right\\}\geq\Pi_{n}\left\\{\norm{\mbox{\boldmath$\Lambda$}_{n}^{*}}_{2}^{2}\geq
K^{\prime}C^{\prime\prime}(q_{0n}s_{n}\log\overline{d_{n}q_{n}})^{4}\right\\}\leq
e^{-CC^{\prime\prime}q_{0n}s_{n}\log(d_{n}q_{n})}.$

From Corollary 1, for arbitrary
$\widetilde{\tau}_{n}>n^{-\nicefrac{{1}}{{2}}}$,
$\Pi_{n}\left\\{B_{n,2}(\mbox{\boldmath$\Omega$}_{0n},\widetilde{\tau}_{n})\right\\}\geq
e^{-Cq_{0n}s_{n}\log(d_{n}q_{n})}$. Hence,
$\frac{\Pi_{n}(\mathcal{P}_{2,n})}{\Pi_{n}\left\\{B_{n,2}(\mbox{\boldmath$\Omega$}_{0n},\widetilde{\tau}_{n})\right\\}}\leq
e^{-2Cq_{0n}s_{n}\log(d_{n}q_{n})}=o(e^{-2n\widetilde{\tau}_{n}^{2}})$ for
$\widetilde{\tau}_{n}<\tau_{n}$. The last display holds for suitable large
enough choice of $C^{\prime\prime}$, which, however, can be chosen independent
of $n$. Using [Theorem 8.20]ghosal_book and subsequent application of DCT we
conclude the proof. ∎

######  Theorem 4 (minimax rate).

If $\widehat{\mbox{\boldmath$\Omega$}}_{n}$ is a sequence of estimators of
$\mbox{\boldmath$\Omega$}_{0n}\in\mathcal{P}_{0n}$ with $q_{0n}=O(1)$, then
for some absolute constant $C>0$

$\displaystyle\inf_{\widehat{\mbox{\boldmath$\Omega$}}_{n}}\sup_{\mbox{\boldmath$\Omega$}_{0n}\in\mathcal{P}_{0n}}\norm{\widehat{\mbox{\boldmath$\Omega$}}_{n}-\mbox{\boldmath$\Omega$}_{0n}}_{2}\geq
C\sqrt{\frac{s_{n}\log d_{n}}{n}}.$

###### Proof.

We use Fano’s lemma to derive a lower bound for the minimax risk. Let ${\cal
F}=\\{\mbox{\boldmath$\Omega$}_{(1)},\dots,\mbox{\boldmath$\Omega$}_{(m_{n})}\\}$,
$m_{n}\geq 2$ be a finite subset of $\mathcal{P}_{0n}$ and let
$\widehat{\mbox{\boldmath$\Omega$}}_{n}$ be an estimate of
$\mbox{\boldmath$\Omega$}_{0n}$ based on iid samples ${\mathbf{y}}_{1:n}$.
Suppose for all $j\neq j^{\prime}$, we have
$e(\mbox{\boldmath$\Omega$}_{(j)},\mbox{\boldmath$\Omega$}_{(j^{\prime})})\geq
e_{m_{n}}$ for some metric $e$ and
$\mathbb{KL}\left(\mbox{\boldmath$\Omega$}_{(j)}\parallel\mbox{\boldmath$\Omega$}_{(j^{\prime})}\right)\leq
K_{m_{n}}$. Letting $\mathbb{E}_{j}$ denote the expectation under
$\mathrm{N}_{d_{n}}({\mathbf{0}},\mbox{\boldmath$\Omega$}_{(j)}^{-1})$, Fano’s
lemma yu1997 implies

$\displaystyle\max_{1\leq j\leq
m_{n}}\mathbb{E}_{j}e(\mbox{\boldmath$\Omega$}_{(j)},\mbox{\boldmath$\Omega$}_{(j^{\prime})})\leq\frac{e_{m_{n}}}{2}\left(1-\frac{K_{m_{n}}+\log
2}{\log m_{n}}\right).$ (S.8)

We now construct the finite class ${\cal F}$. Let $r_{n}=d_{n}-1$. Define
${\cal M}=\\{{\mathbf{x}}\in\mathbb{R}^{r_{n}}:x_{j}\in\\{0,1\\}\text{ for all
}j,~{}\sum_{j}x_{j}=s_{n}\\}$ to be the collection of all binary vectors of
length $r_{n}$ with exactly $s_{n}$ ones and let
$\left\langle{{\mathbf{x}}},{{\mathbf{y}}}\right\rangle_{H}$ be the Hamming
distance between two binary vectors ${\mathbf{x}}$ and ${\mathbf{y}}$. Let
${\mathbf{g}}_{j}=({\mathbf{x}}_{j},0)$ denote the $d_{n}$-dimensional vector
obtained by appending zero at the end of ${\mathbf{x}}_{j}$ with
${\mathbf{x}}_{j}\in{\cal M}$. With this definitions, set
$\mbox{\boldmath$\Omega$}_{(j)}=\beta_{n}{\mathbf{I}}_{d_{n}}+\gamma_{n}{\mathbf{g}}_{j}{\mathbf{g}}_{j}^{\rm
T}+\kappa_{n}\mbox{\boldmath$\varepsilon$}_{d_{n}}\mbox{\boldmath$\varepsilon$}_{d_{n}}^{\rm
T}$ where $\mbox{\boldmath$\varepsilon$}_{d_{n}}$ is the vector with 1 in the
$d_{n}^{th}$ coordinate and zero elsewhere, and
$\gamma_{n}\leq\beta_{n}\leq\kappa_{n}$ are positive sequences to be chosen
below.

Observe that if
$\left\langle{{\mathbf{g}}_{j}},{{\mathbf{g}}_{j^{\prime}}}\right\rangle_{H}=s_{n}-p_{n}$
for $j\neq j$, then ${\mathbf{g}}_{j}^{\rm T}{\mathbf{g}}_{j^{\prime}}=p_{n}$.
Note also that
$\mbox{\boldmath$\Omega$}_{(j)}-\mbox{\boldmath$\Omega$}_{(j^{\prime})}=\gamma_{n}({\mathbf{g}}_{j}{\mathbf{g}}_{j}^{\rm
T}-{\mathbf{g}}_{j^{\prime}}{\mathbf{g}}_{j^{\prime}}^{\rm T})$. The nonzero
eigenvalues of the matrix ${\mathbf{B}}={\mathbf{g}}_{j}{\mathbf{g}}_{j}^{\rm
T}-{\mathbf{g}}_{j^{\prime}}{\mathbf{g}}_{j^{\prime}}^{\rm T}$ are
$(\sqrt{s_{n}^{2}-p_{n}^{2}},-\sqrt{s_{n}^{2}-p_{n}^{2}})$, since
$\mathrm{rank}({\mathbf{B}})=2$, $\mathrm{trace}({\mathbf{B}})=0$ and
$\mathrm{trace}({\mathbf{B}}^{2})=2(s_{n}^{2}-p_{n}^{2})$. This implies that