# Blessing of dimension in Bayesian inference on covariance matrices

Shounak Chattopadhyay<EMAIL_ADDRESS>Department of
Biostatistics, University of California, Los Angeles  Anru R. Zhang
Department of Biostatistics and Bioinformatics, Duke University  Department of
Computer Science, Duke University  Department of Statistical Science, Duke
University  David Dunson  Department of Statistical Science, Duke University
Department of Mathematics, Duke University

###### Abstract

Bayesian factor analysis is routinely used for dimensionality reduction in
modeling of high-dimensional covariance matrices. Factor analytic
decompositions express the covariance as a sum of a low rank and diagonal
matrix. In practice, Gibbs sampling algorithms are typically used for
posterior computation, alternating between updating the latent factors,
loadings, and residual variances. In this article, we exploit a blessing of
dimensionality to develop a provably accurate pseudo-posterior for the
covariance matrix that bypasses the need for Gibbs or other variants of Markov
chain Monte Carlo sampling. Our proposed Factor Analysis with BLEssing of
dimensionality (FABLE) approach relies on a first-stage singular value
decomposition (SVD) to estimate the latent factors, and then defines a jointly
conjugate prior for the loadings and residual variances. The accuracy of the
resulting pseudo-posterior for the covariance improves with increasing
dimensionality. We show that FABLE has excellent performance in high-
dimensional covariance matrix estimation, including producing well calibrated
credible intervals, both theoretically and through simulation experiments. We
also demonstrate the strength of our approach in terms of accurate inference
and computational efficiency by applying it to a gene expression data set.

Keywords: Bayes, Factor analysis, High-dimensional, Large p, Posterior
approximation, Scalable, Singular value decomposition

##  1 Introduction

Inference on covariance matrices in high-dimensional data is a key focus in
many application areas, motivating a rich literature on associated statistical
methods. One thread of this literature avoids modeling of the data and instead
focuses on high-dimensional covariance matrix estimators under various
assumptions on the inherent low-dimensional structure in the data, including
(but not limited to) banded covariance (Bickel and Levina,, 2008), low rank
structure (Shikhaliev et al.,, 2019), low rank with sparsity (Richard et al.,,
2012), sparse covariance (Bien and Tibshirani,, 2011), and sparse inverse
precision matrix estimation (Zhang and Zou,, 2014). Our interest is instead in
model-based Bayesian approaches, which have advantages in terms of naturally
accommodating complexities in the data and quantifying uncertainty, while
having disadvantages in terms of computational efficiency, particularly in
implementing Markov chain Monte Carlo (MCMC) algorithms to sample from the
posterior distribution.

Our particular focus is on factor models, which express the covariance matrix
in a low rank plus diagonal form, and represent a popular probabilistic model-
based alternative to principal components analysis (PCA). Factor analysis
introduces a lower-dimensional latent factor $\eta_{i}\in\mathbb{R}^{k}$
corresponding to each observation $y_{i}\in\mathbb{R}^{p}$ with $k\ll p$, such
that $y_{i}=\Lambda\eta_{i}+\epsilon_{i}$ for a factor loadings matrix
$\Lambda\in\mathbb{R}^{p\times k}$ and independent idiosyncratic errors
$\epsilon_{i}\sim N_{p}(0,\Sigma)$. Following typical practice, we let the
latent factors have independent Gaussian priors $\eta_{i}\sim
N_{k}(0,\mathbb{I}_{k})$. Marginalizing out the latent factors leads to
$y_{i}\sim N_{p}(0,\Psi)$, with the covariance matrix
$\Psi=\Lambda\Lambda^{\top}+\Sigma$. This provides a convenient decomposition
of the covariance, which dramatically reduces the number of free parameters
required to model the $p\times p$ covariance matrix from $p(p+1)/2$ to
$(p-(k-1)/2)k+p$.

There is a vibrant recent literature improving upon and expanding the scope of
Bayesian factor analysis (Schiavon et al.,, 2022; De Vito et al.,, 2021;
Frühwirth-Schnatter,, 2023; Roy et al.,, 2021; Ma and Liu,, 2022; Bolfarine et
al.,, 2022; Xie et al.,, 2022). Even with increasingly rich classes of priors
and data types, the canonical approach for posterior computation remains Gibbs
samplers that iterate between updating latent factors, factor loadings,
residual variances, hyperparameters controlling the hierarchical prior, and
other model parameters. This approach is simple to implement in broad model
classes, but commonly faces problems with slow mixing, particularly as data
dimensionality and complexity increase. Issues with slow mixing and
bottlenecks with large datasets have motivated alternative posterior sampling
algorithms that apply Hamiltonian Monte Carlo (HMC) after marginalizing out
the latent factors, exploiting sufficient statistics and parallel computation
to facilitate implementation for large sample sizes (Chandra et al.,, 2023;
Man and Culpepper,, 2022). Nonetheless, all current MCMC sampling algorithms
for large covariance matrices face substantial computational hurdles as the
number of dimensions increases.

The lack of sufficiently scalable MCMC algorithms has motivated a rich
literature on developing more computationally efficient algorithms for
Bayesian inference in factor models. Variational Bayes approaches (Hansen et
al.,, 2023; Wang and Stephens,, 2021) attempt to approximate the posterior but
often the accuracy of the approximation is poor, with a tendency to
dramatically underestimate posterior uncertainty. A concise review comparing
variational approximations and posterior sampling approaches for factor
analysis is provided in Foo and Shim, (2021). Maximum a posteriori (MAP)
estimation under sparsity priors (Srivastava et al.,, 2017) and expectation-
maximization (EM) algorithms (Ročková and George,, 2016; Avalos-Pacheco et
al.,, 2022; Zhao et al.,, 2016) provide computationally efficient point
estimates of the covariance matrix. Sabnis et al., (2016) develop a
distributed computing framework aggregating local estimates of the covariance
across different cores. Srivastava et al., (2017); Ročková and George, (2016),
and Sabnis et al., (2016) provide theoretical results regarding estimation
accuracy and posterior concentration, without providing any guarantees on
uncertainty quantification. A key open problem in the literature remains how
to obtain an efficient posterior approximation for a high-dimensional
covariance that has guarantees in terms of valid frequentist coverage of
credible intervals.

Fast algorithms for Bayesian factor analysis capable of scaling efficiently to
high-dimensional data typically sacrifice the ability to provide an accurate
characterization of uncertainty. The focus of this article is on proposing a
simple approach for overcoming this limitation, providing computationally
efficient inference on a high-dimensional covariance matrix. With the
increasing dimension, we have more and more variables loading on the same
small number of latent factors. As a consequence, we obtain a blessing of
dimensionality phenomenon allowing us to first pre-estimate the latent factors
and then leverage on this point estimate to obtain a pseudo-posterior
distribution of the factor loadings $\Lambda$ and error variances $\Sigma$
under conjugate priors on these parameters. This in turn induces a pseudo-
posterior distribution for the covariance matrix
$\Psi=\Lambda\Lambda^{\top}+\Sigma$. Our proposed Factor Analysis with
BLEssing of dimensionality (FABLE) approach completely bypasses MCMC. To
ensure that the pseudo-posterior credible intervals for the entries of the
covariance are asymptotically normal with guarantees on valid frequentist
coverage, we develop a data-adaptive coverage correction (CC); we refer to the
adjusted FABLE approach as CC-FABLE.

Fan et al., (2023) also pre-estimate latent factors, but in the context of
regression with high-dimensional correlated covariates. In their approach, a
first-stage principal components analysis is carried out for the high-
dimensional covariates. The results are then used to obtain a modified
covariate vector consisting of principal component scores and de-correlated
versions of the original covariates. Using this new covariate vector in the
second stage high-dimensional regression avoids multicollinearity. This focus,
and the corresponding methodology, are substantially different from ours. A
blessing of dimensionality phenomenon is also explored in Li et al., (2018),
where estimation of the covariance between a subset of the variables of
interest is improved by incorporating auxiliary information from other
variables.

In Section 2, we develop our FABLE methodology with the CC-FABLE modification
and provide an approach to choose key hyperparameter values. In Section 3, we
provide theoretical results on pseudo-posterior contraction rates and
uncertainty quantification, leveraging on a blessing of dimensionality
phenomenon. In Section 4, we validate our approach with comparison to other
state-of-the-art approaches from the viewpoint of estimation error,
frequentist coverage, and computational runtime in different numerical
experiments. In Section 5, we present a relevant application of our approach
to a large gene expression data set and highlight distinct advantages over
competitors. In Section 6, we provide a discussion on future research
directions.

##  2 Proposed Methodology

###  2.1 Initial Approach

The observed data consist of
$\mathbf{Y}=[y_{1},\ldots,y_{n}]^{\top}\in\mathbb{R}^{n\times p}$, where
$y_{i}\in\mathbb{R}^{p}$ for $i=1,\ldots,n$. We consider a factor model for
$y_{i}$ given by

$y_{i}=\Lambda\eta_{i}+\epsilon_{i},\qquad\epsilon_{i}\overset{iid}{\sim}N_{p}(0,\Sigma),$
(1)

where we have omitted the intercept term, assuming the data have been centered
prior to analysis. Here, $\Lambda\in\mathbb{R}^{p\times k}$ is an unknown
matrix of factor loadings,
$\eta_{i}\overset{iid}{\sim}N_{k}(0,\mathbb{I}_{k})$ is the $i$th latent
factor, and $\epsilon_{i}$ is a zero-mean idiosyncratic error having diagonal
covariance $\Sigma=\mbox{diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$ for
$i=1,\ldots,n$. It is typical to assume that the dimension of the latent
factors $k$ is much smaller than the dimension of the data $p$. Integrating
out the latent factors leads to

$y_{i}\mid\Lambda,\Sigma\overset{iid}{\sim}N_{p}(0,\Lambda\Lambda^{\top}+\Sigma).$
(2)

The covariance matrix $\Psi=\Lambda\Lambda^{\top}+\Sigma$ of the $y_{i}$s in
(2) is decomposed as the sum of a low-rank and a diagonal matrix. In this
paper, our goal is to estimate $L=\Lambda\Lambda^{\top},\Sigma,$ and
$\Psi=\Lambda\Lambda^{\top}+\Sigma$. We first illustrate our methodology
assuming $k$ is known. We discuss an approach to estimate $k$ in Section 2.3.

Let $\mathbf{M}=[\eta_{1},\ldots,\eta_{n}]^{\top}\in\mathbb{R}^{n\times k}$
and $\Lambda=[\lambda_{1},\ldots,\lambda_{p}]^{\top}\in\mathbb{R}^{p\times
k}$, with $\lambda_{j}^{\top}$ denoting the $j$th row of $\Lambda$. We also
denote the $j$th column of $\mathbf{Y}$ as $y^{(j)}$, so that
$\mathbf{Y}=[y^{(1)},\ldots,y^{(p)}].$ The latent factor model (1) can be
alternatively expressed as

$y^{(j)}=\mathbf{M}\lambda_{j}+\epsilon^{(j)},$ (3)

where $\epsilon^{(j)}$ is the $j$th column of the matrix
$E=[\epsilon_{1},\ldots,\epsilon_{n}]^{\top}=[\epsilon^{(1)},\ldots,\epsilon^{(p)}],$
with
$\epsilon^{(j)}\overset{ind}{\sim}N_{n}(0_{n},\sigma_{j}^{2}\mathbb{I}_{n})$
for $j=1,\ldots,p$. Writing in matrix form, we obtain

$\mathbf{Y}=\mathbf{M}\Lambda^{\top}+E.$ (4)

Our strategy to obtain computationally efficient pseudo-posterior samples for
$\Psi$ is outlined as follows. We first estimate $\mathbf{M}$ by a plug-in
estimator $\widehat{\mathbf{M}}$ based on the singular value decomposition
(SVD) of $\mathbf{Y}$ and then consider a surrogate model to (3) where
$\mathbf{M}$ is replaced by $\widehat{\mathbf{M}}$. Conditional on the
observed data, $\widehat{\mathbf{M}}$ is fixed. Thus, the surrogate model is
equivalent to $p$ parallel regression problems. We then let the regression
coefficient and the error variance for the $j$th regression have conjugate
normal-inverse gamma prior distributions for $j=1,\ldots,p$ to obtain their
pseudo-posterior samples across the $p$ different regressions in parallel.
Finally, we combine the pseudo-posterior samples of the regression
coefficients and the error variances to obtain a pseudo-posterior sample of
$\Psi$.

To motivate our choice for $\widehat{\mathbf{M}}$, we start with the SVD of
$\mathbf{Y}$, given by

$\mathbf{Y}=UDV^{\top}+U_{\perp}D_{\perp}V_{\perp}^{\top},$ (5)

where $U\in\mathbb{R}^{n\times k},D\in\mathbb{R}^{k\times
k},V\in\mathbb{R}^{p\times
k},U_{\perp}\in\mathbb{R}^{n\times(r-k)},D_{\perp}\in\mathbb{R}^{(r-k)\times(r-k)},$
and $V_{\perp}\in\mathbb{R}^{p\times(r-k)}$, with $r=p\wedge n.$ The columns
of $U,U_{\perp},V,V_{\perp}$ consist of orthonormal vectors, with
$U^{\top}U_{\perp}=V^{\top}V_{\perp}=\mathbf{O}_{k\times(r-k)}$, where
$\mathbf{O}_{p_{1}\times p_{2}}$ denotes the $p_{1}\times p_{2}$ matrix with
all entries equal to $0$. Let $\mathbf{A}=\mathbf{Y}V/\sqrt{p}=UD/{\sqrt{p}}$
and $\widehat{C}\in\mathbb{R}^{k\times k}$ be an invertible matrix satisfying

$\widehat{C}\widehat{C}^{\top}=\dfrac{1}{n}\mathbf{A}^{\top}\mathbf{A}=\dfrac{D^{2}}{np}.$
(6)

An immediate choice of $\widehat{C}$ satisfying (6) is given by
$\widehat{C}=D/\sqrt{np}$. Lastly, we let

$\widehat{\mathbf{M}}=\mathbf{A}(\widehat{C}^{\top})^{-1}.$ (7)

For the choice of $\widehat{C}=D/\sqrt{np},$ it is straightforward from (7)
that $\widehat{\mathbf{M}}=\sqrt{n}U.$ We will denote this particular choice
of $\widehat{\mathbf{M}}=\sqrt{n}U$ as the canonical choice of
$\widehat{\mathbf{M}}$, observing that other choices are also possible for
different choices of $\widehat{C}$ satisfying (6). We provide a heuristic
leading to (7) at the end of this Section.

Given a particular choice of $\widehat{\mathbf{M}}$, we consider the following
surrogate model:

$y^{(j)}=\widehat{\mathbf{M}}\widetilde{\lambda}_{j}+\widetilde{\epsilon}^{(j)},\quad\widetilde{\epsilon}^{(j)}\overset{ind}{\sim}N_{n}(0,\widetilde{\sigma}_{j}^{2}\mathbb{I}_{n})$
(8)

for $j=1,\ldots,p$. The model (8) could be interpreted as a version of (3)
with $\widehat{\mathbf{M}}$ substituted for the original matrix of latent
factors $\mathbf{M}$ and new parameters
$\widetilde{\lambda}_{j}\in\mathbb{R}^{k}$, $\widetilde{\sigma}_{j}^{2}>0$ for
$j=1,\ldots,p$. Next, we endow
$(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})_{j=1}^{p}$ with
independent normal-inverse gamma (NIG) priors
$(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})\overset{iid}{\sim}\mbox{NIG}(0_{k},\tau^{2}\mathbb{I}_{k},\gamma_{0}/2,\gamma_{0}\delta_{0}^{2}/2)$.
That is, we let

$\displaystyle\widetilde{\lambda}_{j}\mid\widetilde{\sigma}_{j}^{2}$
$\displaystyle\sim
N_{k}(0_{k},\widetilde{\sigma}_{j}^{2}\,\tau^{2}\,\mathbb{I}_{k}),\quad\widetilde{\sigma}_{j}^{2}\sim\mbox{IG}\left(\dfrac{\gamma_{0}}{2},\dfrac{\gamma_{0}\delta_{0}^{2}}{2}\right).$
(9)

The global shrinkage parameter $\tau^{2}$ allows us to a priori shrink the
factor loadings towards zero, regularizing
$\widetilde{\Lambda}=[\widetilde{\lambda}_{1},\ldots,\widetilde{\lambda}_{p}]^{\top}$
and favoring sparsity. We discuss a data-driven strategy to estimate
$\tau^{2}$ in Section 2.3.

The surrogate model (8) and the conjugate prior specification (9) motivate the
pseudo-posterior densities $\widetilde{\Pi}_{j}$ of
$(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})$ for each
$j=1,\ldots,p$, given by

$\displaystyle\begin{split}\widetilde{\Pi}_{j}(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})&=\mbox{NIG}\left(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2}\mid\mu_{j},\mathbf{K},\gamma_{n}/2,\gamma_{n}\delta_{j}^{2}/2\right)\\\
&=N_{k}\left(\widetilde{\lambda}_{j}\mid\mu_{j},\widetilde{\sigma}_{j}^{2}\mathbf{K}\right)\mbox{IG}\left(\widetilde{\sigma}_{j}^{2}\,\middle|\,\dfrac{\gamma_{n}}{2},\dfrac{\gamma_{n}\delta_{j}^{2}}{2}\right),\end{split}$
(10)

with the updated hyperparameters given by

$\displaystyle\begin{split}\mu_{j}&=\left(\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}+\dfrac{\mathbb{I}_{k}}{\tau^{2}}\right)^{-1}\widehat{\mathbf{M}}^{\top}y^{(j)},\\\
\mathbf{K}&=\left(\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}+\dfrac{\mathbb{I}_{k}}{\tau^{2}}\right)^{-1},\\\
\gamma_{n}&=\gamma_{0}+n,\\\
\gamma_{n}\delta_{j}^{2}&=\gamma_{0}\delta_{0}^{2}+\left(y^{(j)\top}y^{(j)}-\mu_{j}^{\top}\mathbf{K}^{-1}\mu_{j}\right).\end{split}$
(11)

The pseudo-posteriors are simply obtained by treating $\widehat{\mathbf{M}}$
as a fixed design matrix for the $j$th surrogate regression and then applying
Bayes’ rule. To obtain independent pseudo-posterior samples of $\Psi$, we
first draw independent samples
$(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})\overset{ind}{\sim}\widetilde{\Pi}_{j}$
for $j=1,\ldots,p$. We then let
$\widetilde{\Lambda}=[\widetilde{\lambda}_{1},\ldots,\widetilde{\lambda}_{p}]^{\top},\widetilde{L}=\widetilde{\Lambda}\widetilde{\Lambda}^{\top},\widetilde{\Sigma}=\mbox{diag}(\widetilde{\sigma}_{1}^{2},\ldots,\widetilde{\sigma}_{p}^{2}),$
and $\widetilde{\Psi}=\widetilde{L}+\widetilde{\Sigma}$ denote the pseudo-
posterior samples of $\Lambda,L=\Lambda\Lambda^{\top}$, $\Sigma$, and $\Psi$,
respectively. Our proposed Factor Analysis with BLEssing of dimensionality
(FABLE) approach obtains independent pseudo-posterior samples of the
covariance matrix $\Psi$ in an embarrassingly parallel fashion, bypassing the
need to carry out MCMC.

A natural concern is whether different choices of $\widehat{\mathbf{M}}$
satisfying (7) would affect the distribution of $\widetilde{\Psi}$. The
following result, with proof in Section A of the Supplementary Material,
ensures that this is not the case.

######  Proposition 2.1.

(i) For any $\widehat{\mathbf{M}}$ satisfying (7), we have
$\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}=n\mathbb{I}_{k}$ and
$\widehat{\mathbf{M}}\widehat{\mathbf{M}}^{\top}=nUU^{\top}.$ (ii) The pseudo-
posterior distribution of $\widetilde{L},\widetilde{\Sigma},$ and
$\widetilde{\Psi}$ obtained from the FABLE approach only depends on
$\widehat{\mathbf{M}}$ through
$\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}$ and
$\widehat{\mathbf{M}}\widehat{\mathbf{M}}^{\top}$.

In our experience, the pseudo-posterior samples $\widetilde{\Psi}$ perform
well when estimating the true covariance matrix, both in terms of simulation
studies and theoretical posterior contraction rates. However, we observed that
the entry-wise pseudo-posterior credible intervals obtained from the samples
$\widetilde{\Psi}$ tend to underestimate the uncertainty associated with
estimating the entry-wise elements of the true covariance matrix, from the
perspective of frequentist coverage. In Section 2.2, we build on the initial
FABLE approach developed in this subsection and propose a coverage correction
(CC) of FABLE called CC-FABLE. We provide provably accurate guarantees on
posterior contraction rates and frequentist coverage of the CC-FABLE procedure
in Section 3.

Lastly, we describe the heuristic leading to (7). From (4), the matrix
$\mathbf{A}=\mathbf{Y}V/\sqrt{p}$ satisfies

$\mathbf{A}=\mathbf{M}\dfrac{\Lambda^{\top}V}{\sqrt{p}}+\dfrac{EV}{\sqrt{p}}.$

Let $C=V^{\top}\Lambda/\sqrt{p}\in\mathbb{R}^{k\times k}$. Based on the
consistency of spectral estimates, we expect $E$ to be approximately
independent of $V$ as both $n,p$ grow. As a result, we expect
$EV/\sqrt{p}\approx 0$ for increasing $p$. This leads us to

$\mathbf{A}\approx\mathbf{M}C^{\top}$ (12)

or equivalently, $a_{i}\approx C\eta_{i}$ for $i=1,\ldots,n$, where $a_{i}$ is
the $i$th row of $\mathbf{A}$ for $i=1,\ldots,n$. Since $\eta_{i}\sim
N_{k}(0,\mathbb{I}_{k}),$ the marginal density of $a_{i}$ is approximately
given by $a_{i}\overset{iid}{\sim}N_{k}(0_{k},CC^{\top})$, for $i=1,\ldots,n$.
This motivates the following estimator $\widehat{C}\widehat{C}^{\top}$ of
$CC^{\top}$:

$\widehat{C}\widehat{C}^{\top}=\dfrac{1}{n}\mathbf{A}^{\top}\mathbf{A}=\dfrac{1}{np}D^{2},$

leading us to (6). Given any $\widehat{C}$ satisfying (6), we use (12) to
propose $\widehat{\mathbf{M}}=\mathbf{A}(\widehat{C}^{\top})^{-1}$ as an
estimate of $\mathbf{M}.$ In particular, the canonical choice
$\widehat{\mathbf{M}}=\sqrt{n}U$ is the same as the estimator obtained by
carrying out principal components analysis (PCA) on the matrix
$\mathbf{Y}\mathbf{Y}^{\top}$ (Fan et al.,, 2023). However, (7) permits other
choices corresponding to different choices of $\widehat{C}$ satisfying (6),
providing a general framework for obtaining estimates of the latent factors.

###  2.2 Coverage Correction

Although the pseudo-posterior of $\widetilde{\Psi}$ obtained from (10)
performed well in terms of point estimation across simulation studies, we
observed under-coverage of the entry-wise credible intervals, particularly
when the signal-to-noise ratios $R_{j}=\|\lambda_{j}\|_{2}^{2}/\sigma_{j}^{2}$
for $j=1,\ldots,p$ are large. To improve coverage, we propose a simple
coverage correction (CC) of the FABLE approach, which we denote by CC-FABLE.
We provide a theoretical justification of this approach in Section 3. We found
that the correction leads to improved performance in simulation studies,
particularly with increasing $R_{j}$s.

To improve coverage, we include an additional factor $\rho^{2}$ in the
conditional pseudo-posterior variance when generating
$\widetilde{\lambda}_{j}$ for $j=1,\ldots,p$ following (10). That is, we
modify $\widetilde{\Pi}_{j}$ in (10) to

${\Pi}_{j}^{C}\left(\lambda_{C,j},\sigma_{C,j}^{2}\right)=N_{k}\left(\lambda_{C,j}\middle|\mu_{j},\rho^{2}\sigma_{C,j}^{2}\mathbf{K}\right)\mbox{IG}\left(\sigma_{C,j}^{2}\middle|\dfrac{\gamma_{n}}{2},\dfrac{\gamma_{n}\delta_{j}^{2}}{2}\right),$
(13)

for $j=1,\ldots,p$, where
$\mathbf{K}=\\{\mathbf{M}^{\top}\mathbf{M}+(\mathbb{I}_{k}/\tau^{2})\\}^{-1}=\mathbb{I}_{k}/(n+\tau^{-2})$
and $\gamma_{n},\mu_{j},\delta_{j}^{2}$ are as in (11). We then proceed as
before to draw coverage-corrected pseudo-posterior samples
$(\lambda_{C,j},\sigma_{C,j}^{2})\overset{ind}{\sim}\Pi_{j}^{C}$ for
$j=1,\ldots,p$. Letting
$\Lambda_{C}=[\lambda_{C,1},\ldots,\lambda_{C,p}]^{\top}$ and
$\Sigma_{C}=\mbox{diag}(\sigma_{C,1}^{2},\ldots,\sigma_{C,p}^{2})$, we
construct a coverage-corrected (CC) pseudo-posterior sample of the covariance
matrix as $\Psi_{C}=\Lambda_{C}\Lambda_{C}^{\top}+\Sigma_{C}.$ Having computed
the hyperparameters $\gamma_{n}$ and $\delta_{j}^{2}$ for $j=1,\ldots,p$ as in
(11), the pseudo-posterior mean $\widehat{\Psi}_{C}$ of $\Psi_{C}$ is
available analytically, with its $uv$th entry for $1\leq u,v\leq p$ given by

$\displaystyle\begin{split}\widehat{\Psi}_{C,uv}&=\left\\{\begin{array}[]{ll}\mu_{u}^{\top}\mu_{v},&\text{if
$u\neq v,$}\\\
\|\mu_{u}\|_{2}^{2}+\left(1+\dfrac{k\rho^{2}}{n+\tau^{-2}}\right)\dfrac{\gamma_{n}\delta_{u}^{2}}{\gamma_{n}-2},&\text{if
$u=v,$}\end{array}\right.\end{split}$ (14)

provided $\gamma_{n}>2$. For finite $k,\rho^{2},$ and $\tau^{2}$, the pseudo-
posterior mean $\widehat{\Psi}_{C}$ can approximated as
$\widehat{\Psi}_{C}=G_{0}G_{0}^{\top}+\Delta,$ where
$G_{0}=[\mu_{1},\ldots,\mu_{p}]^{\top}$ and
$\Delta=\mbox{diag}(\delta_{1}^{2},\ldots,\delta_{p}^{2})$, as $n$ increases.
From here on, we take $\widehat{\Psi}_{C}=G_{0}G_{0}^{\top}+\Delta$ to be the
pseudo-posterior mean of $\Psi$ under CC-FABLE, unless mentioned otherwise.

We choose $\rho$ to ensure that a summary statistic of the entrywise
asymptotic coverages is $1-\alpha$, such as their average or minimum. Let $B$
denote the $p\times p$ symmetric matrix with entries $b_{uv}$, where

$\displaystyle\begin{split}b_{uv}&=\left\\{\begin{array}[]{ll}\left(1+\dfrac{\|\mu_{u}\|_{2}^{2}\|\mu_{v}\|_{2}^{2}+(\mu_{u}^{\top}\mu_{v})^{2}}{\mathcal{V}_{u}^{2}\|\mu_{v}\|^{2}+\mathcal{V}_{v}^{2}\|\mu_{u}\|_{2}^{2}}\right)^{1/2},&\text{if
$u\neq v,$}\\\
\left(1+\dfrac{\|\mu_{u}\|_{2}^{2}}{2\mathcal{V}_{u}^{2}}\right)^{1/2},&\text{if
$u=v,$}\end{array}\right.\end{split}$ (15)

with $\mathcal{V}_{j}^{2}=\|(\mathbb{I}_{n}-UU^{\top})y^{(j)}\|_{2}^{2}/n$ for
$j=1,\ldots,p$. For example, we can ensure that all the entrywise asymptotic
coverages are at least $1-\alpha$ by letting $\rho=\sup_{u,v}b_{uv}$. Another
possibility is to let $\rho$ be the root of a nonlinear equation, which
ensures the average of the asymptotic entrywise coverages is $1-\alpha$; we
discuss the associated nonlinear equation along with the rationale behind
possible choices of $\rho$ at the end of Section 3.3. Based on our experiences
with finite sample simulations, we recommend choosing $\rho$ to the mean of
$b_{uv}$ for $1\leq u\leq v\leq p$ as a computationally efficient strategy.
This performed well with both median and average entrywise coverages close to
$1-\alpha$. In what follows, we only consider the CC-FABLE version of our
approach for theoretical justification, simulation results, and the real data
application.

###  2.3 Hyperparameter Choice

Before obtaining pseudo-posterior samples of the covariance matrix, one needs
to carefully choose two key hyperparameters. These are the number of factors
$k$ and the global variance of the factor loadings $\tau^{2}$. To estimate
$k$, we use the approach described in Chen and Li, (2022). The authors
construct an information criterion $\mbox{JIC}(k)$ by adding the penalty
$\nu(k,n,p)=k(n\vee p)\log(n\wedge p)$ to the model deviance. The proposed
$\mbox{JIC}(k)$ criterion is similar in spirit to the Akaike information
criterion (AIC) or the Bayesian information criterion (BIC) (Akaike,, 1974;
Schwarz,, 1978). Given an upper bound $\mathcal{K}_{0}$ for the possible
values of $k$, we estimate $\widehat{k}=\min_{1\leq
k\leq\mathcal{K}_{0}}\mbox{JIC}(k).$ To choose an upper bound
$\mathcal{K}_{0}$ for plausible values of $k$, we consider the singular values
$s_{1}\geq s_{2}\geq\ldots\geq s_{n\wedge p}$ of $\mathbf{Y}$ and let
$\mathcal{K}_{0}$ be the smallest integer $\mathcal{K}$ such that
$\left(\sum_{j=1}^{\mathcal{K}}s_{j}\right)/\left(\sum_{j=1}^{n\wedge
p}s_{j}\right)\geq S_{0}$, where we typically let $S_{0}$ be $0.5$ or $0.75$.
In our experience across simulation studies and applications to gene
expression datasets, we found this choice of $\widehat{k}$ to perform well.

We next proceed to estimate $\tau^{2}$, which can be interpreted as the global
shrinkage parameter of the loadings. For this, we employ an empirical Bayes
approach. Since
$\widetilde{\lambda}_{j}\mid\widetilde{\sigma}_{j}^{2},\tau^{2}\sim
N_{k}(0,\tau^{2}\widetilde{\sigma}_{j}^{2}\mathbb{I}_{k})$ for $j=1,\ldots,p$
a priori, we obtain an estimate of $\tau^{2}$ by conditioning on
$(\widetilde{\lambda}_{j},\widetilde{\sigma}_{j}^{2})$ and then maximizing the
conditional likelihood, leading to
$\widehat{\tau}_{0}^{2}=\sum_{j=1}^{p}\left({\|\widetilde{\lambda}_{j}\|_{2}^{2}}/{\widetilde{\sigma}_{j}^{2}}\right)/(kp).$
Following from the pseudo-posterior updates (10) and (11), we simply estimate
$k$ by $\widehat{k}$ as before, $\|\widetilde{\lambda}_{j}\|_{2}^{2}$ by
$\mathcal{L}_{j}^{2}=\|U^{\top}y^{(j)}\|_{2}^{2}/n$, and
$\widetilde{\sigma}_{j}^{2}$ by
$\mathcal{V}_{j}^{2}=\|(\mathbb{I}_{n}-UU^{\top})y^{(j)}\|_{2}^{2}/n$ for
$j=1,\ldots,p$. We remark here that $U$ is dependent on the estimate
$\widehat{k}$. This leads us to the plug-in estimate of $\tau^{2}$ given by

$\widehat{\tau}^{2}=\dfrac{1}{p\widehat{k}}\sum_{j=1}^{p}\dfrac{\mathcal{L}_{j}^{2}}{\mathcal{V}_{j}^{2}}.$
(16)

The consistency of the point estimators $\mathcal{L}_{j}^{2}$ and
$\mathcal{V}_{j}^{2}$ for $1\leq j\leq p$ is proven in Section E, Lemma E.4 of
the Supplementary Material. Lastly, we observed that inference using the CC-
FABLE procedure is not sensitive to $(\gamma_{0},\delta_{0}^{2})$, and thus we
use the default values $\gamma_{0}=\delta_{0}^{2}=1$ in practice. Once the
number of factors $k$ and the common variance $\tau^{2}$ have been estimated
using $\widehat{k}$ and $\widehat{\tau}^{2}$, respectively, we then proceed to
obtain pseudo-posterior samples of $\Psi$ using CC-FABLE as described before.

###  2.4 Final Algorithm

In this subsection, we provide an algorithm for implementing CC-FABLE on a
given data set to obtain pseudo-posterior samples of the high-dimensional
covariance matrix. We first tune the relevant hyperparameters $k$ and
$\tau^{2}$ as in Section 2.3 and proceed to obtain coverage-corrected samples
as in Section 2.2. The relevant steps are described in Algorithm 1.

Input: The data matrix $\mathbf{Y}\in\mathbb{R}^{n\times p}$, number of Monte
Carlo (MC) samples $N_{0}$, inverse-gamma hyperparameters
$(\gamma_{0},\sigma_{0}^{2}),$ and the upper bound on the cumulative singular
value proportion ${S}_{0}$. Let $r=n\wedge p$.

Step 1: Compute the SVD of $\mathbf{Y}$ as $\mathbf{Y}=U^{*}D^{*}V^{*\top}$
with $U\in\mathbb{R}^{n\times r},D\in\mathbb{R}^{r\times
r},V\in\mathbb{R}^{p\times r}$, and $D=\mbox{diag}(s_{1},\ldots,s_{r})$ such
that $s_{1}\geq s_{2}\geq\ldots\geq s_{r}\geq 0.$

Step 2: Let
$\widehat{k}={\mbox{argmin}}\,\\{\mbox{JIC}(k^{*})\,|\,k^{*}=1,\ldots,\mathcal{K}_{0}\\},$
where $\mathcal{K}_{0}=\mbox{argmin}_{1\leq\mathcal{K}\leq
r}\left(\sum_{j=1}^{\mathcal{K}}s_{j}\right)/\left(\sum_{j=1}^{r}s_{j}\right).$

Step 3: Let $U$ consist of the first $\widehat{k}$ columns of $U^{*}$. For
$1\leq j\leq p$, let $y^{(j)}$ denote the $j$th column of $\mathbf{Y}$, and
obtain $\mathcal{L}_{j}^{2}=\|U^{\top}y^{(j)}\|_{2}^{2}/n$ and
$\mathcal{V}_{j}^{2}=\|(\mathbb{I}_{n}-UU^{\top})y^{(j)}\|_{2}^{2}/n.$

Step 4: Estimate $\tau^{2}$ by
$\widehat{\tau}^{2}=\dfrac{1}{p\widehat{k}}\sum_{j=1}^{p}\dfrac{\mathcal{L}_{j}^{2}}{\mathcal{V}_{j}^{2}}.$

Step 5: For $1\leq j\leq p$, let
$\mu_{j}=\sqrt{n}U^{\top}y^{(j)}/\left(n+\widehat{\tau}^{-2}\right).$ For
$1\leq u,v\leq p$, let
$\displaystyle\begin{split}b_{uv}&=\left\\{\begin{array}[]{ll}\left(1+\dfrac{\|\mu_{u}\|_{2}^{2}\|\mu_{v}\|_{2}^{2}+(\mu_{u}^{\top}\mu_{v})^{2}}{\mathcal{V}_{u}^{2}\|\mu_{v}\|^{2}+\mathcal{V}_{v}^{2}\|\mu_{u}\|_{2}^{2}}\right)^{1/2},&\text{if
$u\neq v,$}\\\
\left(1+\dfrac{\|\mu_{u}\|_{2}^{2}}{2\mathcal{V}_{u}^{2}}\right)^{1/2},&\text{if
$u=v,$}\end{array}\right.\end{split}$ and let
$\rho=\overline{b}=\left(\sum_{1\leq u\leq v\leq
p}b_{uv}\right)/\\{p(p+1)/2\\}$ be the coverage-correction factor.

Step 6: Let $\gamma_{n}=\gamma_{0}+n$. For $j=1,\ldots,p$, evaluate
$\gamma_{n}\delta_{j}^{2}=\gamma_{0}\delta_{0}^{2}+y^{(j)\top}\left(\mathbb{I}_{n}-\dfrac{nUU^{\top}}{n+\widehat{\tau}^{-2}}\right)y^{(j)}.$

Step 7: For each $t=1,\ldots,N_{0}$, independently sample
$({\lambda}_{C,j}^{(t)},{\sigma}_{C,j}^{(t)2})$ for $j=1,\ldots,p$ from
${\sigma}_{C,j}^{(t)2}\sim\mbox{IG}\left(\dfrac{\gamma_{n}}{2},\dfrac{\gamma_{n}\delta_{j}^{2}}{2}\right),\qquad{\lambda}_{C,j}^{(t)}\mid{\sigma}_{C,j}^{(t)2}\sim
N_{k}\left(\mu_{j},\dfrac{\rho^{2}{\sigma}_{C,j}^{(t)2}}{n+\widehat{\tau}^{-2}}\mathbb{I}_{k}\right).$

Form
${\Lambda}_{C}^{(t)}=[{\lambda}_{C,1}^{(t)},\ldots,{\lambda}_{C,p}^{(t)}]^{\top}$
and
${\Sigma}_{C}^{(t)}=\mbox{diag}({\sigma}_{C,1}^{(t)2},\ldots,{\sigma}_{C,p}^{(t)2})$.

Step 8: For each $t=1,\ldots,N_{0}$, let
$L_{C}^{(t)}=\Lambda_{C}^{(t)}\Lambda_{C}^{(t)\top}$, and compute the $t$th
coverage-corrected sample as  $\Psi_{C}^{(t)}=L_{C}^{(t)}+\Sigma_{C}^{(t)}.$

Output: The $N_{0}$ MC samples of the covariance matrix
$\Psi_{C}^{(1)},\ldots,\Psi_{C}^{(N_{0})}$.

Algorithm 1 Steps to obtain $N_{0}$ pseudo-posterior samples of the covariance
matrix using CC-FABLE.

##  3 Theoretical Support

###  3.1 Setup and Assumptions

In this section, we provide theoretical guarantees of the proposed CC-FABLE
procedure in estimating the high-dimensional covariance matrix. Most of the
existing literature (Pati et al.,, 2014; Bhattacharya and Dunson,, 2011) on
obtaining provable guarantees for Bayesian factor models utilizes the
machinery of Bayes’ rule in deriving posterior contraction rates. There is a
clear need for further theoretical work on uncertainty quantification in such
models. Xie and Wu, (2023) consider a signal-plus-noise model and derive
theoretical properties of a quasi-posterior based approach, with guarantees on
asymptotic frequentist coverage of the quasi-posterior intervals. However,
their model only considers a square signal matrix, with the underlying focus
different from high-dimensional covariance estimation. In contrast, the
proposed method carries out inference on the underlying covariance matrix,
without directly relying on Bayes’ rule to obtain the pseudo-posterior
samples, which makes theoretical analysis of the proposed approach
challenging. We overcome this challenge by leveraging on a blessing of
dimensionality phenomenon which requires both $n$ and $p=p_{n}$ to grow,
providing results on both the pseudo-posterior contraction and uncertainty
quantification of our approach. The blessing of dimensionality phenomenon
allows accurate estimation of the latent factor subspace up to rotational
ambiguity, which serves as an important component behind the theoretical
results. The proofs of all the theorems can be found in Sections B-D of the
Supplementary Material.

We assume the following data-generating model:

$y_{i}=\Lambda_{0}\eta_{0i}+\epsilon_{i},$ (17)

where $\epsilon_{i}\overset{iid}{\sim}N_{p}(0,\Sigma_{0})$ for $i=1,\ldots,n$
with $\Sigma_{0}=\mbox{diag}(\sigma_{01}^{2},\ldots,\sigma_{0p}^{2}),$ and
$\eta_{0i}\overset{iid}{\sim}N_{k}(0,\mathbb{I}_{k})$ for $i=1,\ldots,n$.
Here, $\Lambda_{0}=[\lambda_{01},\ldots,\lambda_{0p}]^{\top}$ is the true
factor loadings matrix while $\eta_{0i}$ are the true latent factors;
integrating them out provides us
$y_{i}\overset{iid}{\sim}N_{p}(0,\Lambda_{0}\Lambda_{0}^{\top}+\Sigma_{0})$
for $i=1,\ldots,n$ as the marginal distribution of the data. Let
$\mathbf{Y}=[y_{1},\ldots,y_{n}]^{\top}\in\mathbb{R}^{n\times p}$ and
$M_{0}=[\eta_{01},\ldots,\eta_{0n}]^{\top}\in\mathbb{R}^{n\times k}$ be the
data matrix and the true matrix of latent factors, respectively, so that the
true data generating model may be written as

$\mathbf{Y}=M_{0}\Lambda_{0}^{\top}+E,$ (18)

where $E=[\epsilon_{1},\ldots,\epsilon_{n}]^{\top}.$ Our primary goal is the
inference for the covariance matrix
$\Psi_{0}=\Lambda_{0}\Lambda_{0}^{\top}+\Sigma_{0}.$

For a matrix $A\in\mathbb{R}^{n_{1}\times n_{2}}$, suppose the singular values
of $A$ are given by $s_{1}(A)\geq\ldots\geq s_{n_{1}\wedge n_{2}}(A)$; let
$\|A\|=\underset{\|x\|_{2}=1}{\sup}\|Ax\|_{2}=s_{1}(A)$ denote the operator
norm of $A$ and $\|A\|_{\infty}=\max_{ij}|A_{ij}|$ denote the max-norm of $A$,
with $A_{ij}$ denoting the $ij$th entry of $A$. For two non-negative sequences
$a_{m},b_{m}$, we say $a_{m}\asymp b_{m}$ if $a_{m}=\mathcal{O}(b_{m})$ and
$b_{m}=\mathcal{O}(a_{m})$ as $m\to\infty$. We assume the following conditions
on the true data-generating model:

######  Assumption 3.1.

We have $p_{n}\to\infty$ and $(\log p_{n})/n=o(1)$ as $n\to\infty$.

######  Assumption 3.2.

$s_{k}(\Lambda_{0})\asymp\|\Lambda_{0}\|\asymp\sqrt{p_{n}}$ as $n\to\infty$,
$\|\Lambda_{0}\|_{\infty}<\infty$, and $\underset{1\leq j\leq
p_{n}}{\min}\|\lambda_{0j}\|_{2}\geq c_{1}$ for some constant $c_{1}>0$.

######  Assumption 3.3.

The true error variances satisfy $\underset{1\leq j\leq
p_{n}}{\max}\sigma_{0j}^{2}=\mathcal{O}(1)$ and $\underset{1\leq j\leq
p_{n}}{\min}\sigma_{0j}^{2}\geq c_{2}$ for some constant $c_{2}>0$.

######  Assumption 3.4.

The hyperparameters $k,\tau^{2},\gamma_{0},\delta_{0}^{2}$, and $\rho$ are
fixed constants.

Such assumptions are standard in the existing literature on asymptotic
theoretical properties of latent factor models. Assumption 3.1 allows the
number of dimensions $p_{n}$ to scale as any polynomial function of $n$.
Assumption 3.2 ensures that the true loadings matrix $\Lambda_{0}$ is well-
conditioned, with the low-rank portion $\Lambda_{0}\Lambda_{0}^{\top}$
identifiable from noise in the asymptotic regime, and also ensures that none
of the columns of $\mathbf{Y}$ consist purely of noise. Assumptions 3.3 and
3.4 assume the scalar error variances and model hyperparameters are finite.
For our theoretical requirements, we will assume that $k$ and $\tau^{2}$ are
known and fixed. To ease exposition, we will often suppress the dependence of
$p=p_{n}$ on $n$ and express it when needed.

Let $M_{0}\Lambda_{0}^{\top}=U_{0}D_{0}V_{0}^{\top}$ be the singular value
decomposition of the signal, with $U_{0}\in\mathbb{R}^{n\times
k},V_{0}\in\mathbb{R}^{p\times k}$ having orthonormal columns and
$D_{0}\in\mathbb{R}^{k\times k}$ a diagonal matrix of positive singular
values. Suppose the singular value decomposition of $\mathbf{Y}$ is

$\mathbf{Y}=UDV^{\top}+U_{\perp}D_{\perp}V_{\perp}^{\top},$

where $U\in\mathbb{R}^{n\times k},V\in\mathbb{R}^{p\times k}$ have orthonormal
columns, and $D\in\mathbb{R}^{k\times k}$ contain the $k$ largest singular
values of $\mathbf{Y}$. We first provide a result showcasing the blessing of
dimensionality when estimating $U_{0}U_{0}^{\top}$ by $UU^{\top}$, which forms
a key part of the results that follow. Let us denote the induced CC-FABLE
pseudo-posterior measure, the true data-generating measure, and the
expectation under the true data-generating measure by ${\Pi}_{C}$, $P_{0}$,
and $E_{0}$, respectively.

######  Proposition 3.5.

Suppose Assumptions 3.1 – 3.4 hold. Then, there exists a constant $G_{1}>0$
such that

$\lim_{n\to\infty}P_{0}\left\\{\|UU^{\top}-U_{0}U_{0}^{\top}\|>G_{1}\left(\dfrac{1}{n}+\dfrac{1}{p_{n}}\right)\right\\}=0.$

The proof of Proposition 3.5 is available in Section E.1 of the Supplementary
Material. First, in Theorem 3.6, we provide psuedo-posterior contraction rates
when estimating $\Lambda_{0}\Lambda_{0}^{\top},\Sigma_{0},$ and $\Psi_{0}$
using $L_{C},{\Sigma}_{C}$, and ${\Psi}_{C}=L_{C}+\Sigma_{C}$, respectively,
obtained from the CC-FABLE methodology in Section 2.2. Later, using Theorems
3.7 and 3.8, we justify the uncertainty quantification of the entrywise
elements of $\Psi_{0}$ via pseudo-posterior credible intervals obtained from
CC-FABLE.

###  3.2 Pseudo-posterior Contraction Rates

######  Theorem 3.6.

Suppose Assumptions 3.1 – 3.4 hold. Then, as $n\to\infty$,

  1. (a)

There exists a constant $D_{1}>0$ such that

$E_{0}\left[{\Pi}_{C}\left\\{\dfrac{\|L_{C}-\Lambda_{0}\Lambda_{0}^{\top}\|}{\|\Lambda_{0}\Lambda_{0}^{\top}\|}>D_{1}\left(\dfrac{1}{\sqrt{n}}+\dfrac{1}{\sqrt{p_{n}}}\right)\right\\}\right]\to
0.$

  2. (b)

There exists a constant $D_{2}>0$ such that

$E_{0}\left[{\Pi}_{C}\left\\{\|\Sigma_{C}-\Sigma_{0}\|>D_{2}\left(\left(\dfrac{\log
p_{n}}{n}\right)^{1/3}+\dfrac{1}{p_{n}}\right)\right\\}\right]\to 0.$

  3. (c)

There exists a constant $D>0$ such that

$E_{0}\left[{\Pi}_{C}\left\\{\dfrac{\|\Psi_{C}-\Psi_{0}\|}{\|\Psi_{0}\|}>D\left(\dfrac{1}{\sqrt{n}}+\dfrac{1}{\sqrt{p_{n}}}\right)\right\\}\right]\to
0.$

Theorem 3.6 shows pseudo-posterior concentration of the relevant quantities of
the CC-FABLE procedure. The relative errors when estimating each of these true
parameters of interest converge to $0$ at the rate $n^{-1/2}+p_{n}^{-1/2}$,
which showcases the blessing of dimensionality acting on the FABLE procedure.
One obtains pseudo-posterior concentration if and only if both the number of
samples $n$ and the number of dimensions $p_{n}$ increase. The proof of
Theorem 3.6 is available in Section B of the Supplementary Material.

###  3.3 Uncertainty Quantification

We now consider uncertainty quantification of entrywise elements of the true
covariance matrix $\Psi_{0}$. For $1\leq u,v\leq p,$ let ${\Psi}_{C,uv}$ and
$\Psi_{0,uv}$ denote the $uv$th element of ${\Psi}_{C}$ and $\Psi_{0}$,
respectively, where $\Psi_{C}$ denotes the coverage-corrected pseudo-posterior
sample obtained from the CC-FABLE discussed in Section 2.2. For $1\leq u,v\leq
p$, we have
$\Psi_{0,uv}=\lambda_{0u}^{\top}\lambda_{0v}+\sigma_{0u}^{2}\mathbbm{1}(u=v)$,
where $\mathbbm{1}$ denotes the indicator function. In what follows, we
discuss results for a fixed value of the coverage-correction factor $\rho>0$
and then provide an approach to appropriately choose $\rho$.

We now let $T_{uv}=\mu_{u}^{\top}\mu_{v}+\delta_{u}^{2}\mathbbm{1}(u=v)$ be an
estimator of $\Psi_{0,uv}$. Then, we have the following result approximating
the pseudo-posterior distribution of ${\Psi}_{C,uv}$ by a suitable Gaussian
distribution centered at $T_{uv}$ as both $n,p_{n}$ increase with $\rho>0$
fixed.

######  Theorem 3.7.

Suppose Assumptions 3.1 – 3.4 hold and let $\sqrt{n}/p_{n}=o(1)$. For $1\leq
u,v\leq p$, let

$\displaystyle\begin{split}l_{0,uv}^{2}(\rho)&=\left\\{\begin{array}[]{ll}\rho^{2}(\sigma_{0v}^{2}\|\lambda_{0u}\|_{2}^{2}+\sigma_{0u}^{2}\|\lambda_{0v}\|_{2}^{2}),&\text{for
$u\neq v$},\\\
2\sigma_{0u}^{4}+4\rho^{2}\,\sigma_{0u}^{2}\|\lambda_{0u}\|_{2}^{2},&\text{for
$u=v$.}\end{array}\right.\end{split}$

Then, as $n\to\infty,$

$\sup_{x\in\mathbb{R}}\,\left|{\Pi}_{C}\left\\{\dfrac{\sqrt{n}({\Psi}_{C,uv}-T_{uv})}{\l_{0,uv}(\rho)}\leq
x\right\\}-\Phi(x)\right|\overset{P_{0}}{\to}0\quad\text{for all $1\leq
u,v\leq p.$}$

Theorem 3.7 allows us to approximate the asymptotic pseudo-posterior
distribution of each element of the covariance matrix after suitable centering
and scaling using a Gaussian distribution with mean $0$ and variance
$l_{0,uv}^{2}(\rho)$. The proof of Theorem 3.7 is in Section C of the
Supplementary Material. We next state a result regarding the asymptotic law of
the quantity $\sqrt{n}\,(T_{uv}-\Psi_{0,uv})/\mathcal{S}_{0,uv}$, demonstrate
how this result can be used to obtain asymptotic entrywise frequentist
coverage for $\Psi_{0,uv}$ using CC-FABLE, and outline different choices of
$\rho$. For general random variables $X_{n}$ and $X$, we denote $X_{n}$
converging in distribution to $X$ by $X_{n}\implies X.$

######  Theorem 3.8.

Suppose the Assumptions of Theorem 3.7 hold. For $1\leq u,v\leq p$, let

$\displaystyle\begin{split}\mathcal{S}_{0,uv}^{2}&=\left\\{\begin{array}[]{ll}\sigma_{0v}^{2}\|\lambda_{0u}\|_{2}^{2}+\sigma_{0u}^{2}\|\lambda_{0v}\|_{2}^{2}+\|\lambda_{0u}\|_{2}^{2}\|\lambda_{0v}\|_{2}^{2}+(\lambda_{0u}^{\top}\lambda_{0v})^{2},&\text{for
$u\neq v$},\\\ 2(\|\lambda_{0u}\|_{2}^{2}+\sigma_{0u}^{2})^{2},&\text{for
$u=v$.}\end{array}\right.\end{split}$

Then, as $n\to\infty$, one has
$\sqrt{n}(T_{uv}-\Psi_{0,uv})/\mathcal{S}_{0,uv}\implies N(0,1)$ for all
$1\leq u,v\leq p.$

The proof of Theorem 3.8 is available in Section D of the Supplementary
Material. Under Theorem 3.7, the $100(1-\alpha)\%$ asymptotic pseudo-posterior
credible interval of ${\Psi}_{C,uv}$ is given by

$\mathcal{C}_{uv}(\rho)=\left[T_{uv}\pm
z_{1-(\alpha/2)}\dfrac{l_{0,uv}(\rho)}{\sqrt{n}}\right],$

where $z_{1-(\alpha/2)}=\Phi^{-1}\left\\{1-(\alpha/2)\right\\}.$ Thus, the
probability of $\mathcal{C}_{uv}(\rho)$ covering $\Psi_{0,uv}$ under repeated
sampling is given by

$\displaystyle P_{0}\left\\{\Psi_{0,uv}\in\mathcal{C}_{uv}(\rho)\right\\}$
$\displaystyle=P_{0}\left\\{\dfrac{\sqrt{n}\,\left|T_{uv}-\Psi_{0,uv}\right|}{\mathcal{S}_{0,uv}}\leq
z_{1-(\alpha/2)}\dfrac{l_{0,uv}(\rho)}{\mathcal{S}_{0,uv}}\right\\}$
$\displaystyle\to
q_{uv}(\rho):=2\Phi\left\\{z_{1-(\alpha/2)}\dfrac{l_{0,uv}(\rho)}{\mathcal{S}_{0,uv}}\right\\}-1,$

as $n\to\infty,$ using Theorem 3.8.

The form of $q_{uv}(\rho)$ suggests a general data-adaptive framework to
choose $\rho$. For example, we can ensure that all the entrywise asymptotic
coverages are at least $1-\alpha$ by letting $\rho=\sup_{u,v}b_{uv}$ for
$1\leq u\leq v\leq p$, with $b_{uv}$ defined in Section 2.2. To see why, first
let $\hat{q}_{uv}(\rho)$ be the consistent estimator of $q_{uv}(\rho)$
obtained by substituting $\sigma_{0j}^{2}$ and $\lambda_{0j}$ with
$\mathcal{V}_{j}^{2}$ and $\mu_{j}$, respectively, for $j=1,\ldots,p$; these
quantities are defined in Section 2. Next, observe that $\hat{q}_{uv}(\rho)$
is a non-decreasing function of $\rho>0$ for all $1\leq u,v\leq p$. It is
clear that $\hat{q}_{uv}(b_{uv})=1-\alpha$ for all $1\leq u\leq v\leq p.$
Thus,
$\hat{q}_{uv}\left(\sup_{u^{\prime},v^{\prime}}b_{u^{\prime}v^{\prime}}\right)\geq\hat{q}_{uv}(b_{uv})=1-\alpha.$
Alternatively, we can ensure the average of the entrywise asymptotic coverages
is $1-\alpha$ by numerically solving

$\dfrac{2}{p(p+1)}\displaystyle\sum_{1\leq u\leq v\leq
p}\hat{q}_{uv}(\rho)=1-\alpha$ (19)

for $\rho$. We observed in simulation experiments that simply letting

$\rho=\overline{b}=\left(\sum_{1\leq u\leq v\leq
p}b_{uv}\right)/\left\\{\dfrac{p(p+1)}{2}\right\\}$

is a close approximation to the root of (19). In practice, we thus recommend
choosing $\rho=\overline{b}$ for computational convenience, as this avoids
solving the nonlinear equation. This worked well in the simulations we carried
out, both from the perspective of median and average asymptotic entrywise
coverages. Since Lemma E.5 in Section E of the Supplementary Material ensures
that the matrix $B=(b_{uv})_{1\leq u,v\leq p}$ satisfies
$\|B\|_{\infty}=O_{P_{0}}(1)$, all the aforementioned estimates of $\rho$ are
bounded.

##  4 Simulation Results

###  4.1 Setup

In this Section, we compare the proposed CC-FABLE approach with competitors
when judging their performance from the viewpoint of estimation error,
uncertainty quantification (UQ), and computational efficiency. With
$\Lambda_{0}=(\Lambda_{0,uv})_{1\leq u,v\leq p}$ the true matrix of loadings
and $\Sigma_{0}=\mbox{diag}(\sigma_{01}^{2},\ldots,\sigma_{0p}^{2})$ the true
matrix of idiosyncratic error variances, the data generating covariance matrix
is given by $\Psi_{0}=\Lambda_{0}\Lambda_{0}^{\top}+\Sigma_{0}$. For both the
estimation error and UQ experiments, we proceed as follows. We vary the number
of samples $n\in\\{500,1000\\}$, the number of dimensions
$p\in\\{1000,5000\\}$, and let the true number of factors be $k=10.$ For each
possible $(n,p)\in\\{500,1000\\}\times\\{1000,5000\\}$, we generate
$(\Lambda_{0},\Sigma_{0})$ as

$\displaystyle\Lambda_{0,jl}$
$\displaystyle\overset{ind}{\sim}0.5\,\widetilde{\delta}_{0}+0.5\,N(0,0.5^{2})$
for all $1\leq j\leq p,1\leq l\leq k$, and $\displaystyle\sigma_{0j}^{2}$
$\displaystyle\overset{ind}{\sim}\mathcal{U}(0.5,5),$ for all $1\leq j\leq p,$

where $\widetilde{\delta}_{0}$ denotes a point mass at $0$. For each
combination $(n,p)$, we use the generated
$\Psi_{0}=\Lambda_{0}\Lambda_{0}^{\top}+\Sigma_{0}$ and replicate the data
generating process $R=100$ times. Given an estimator $\widehat{\Psi}_{0}$ of
$\Psi_{0}$, we judge the estimator using the relative spectral error, defined
as

$\mathcal{L}(\Psi_{0},\widehat{\Psi}_{0})=\dfrac{\|\Psi_{0}-\widehat{\Psi}_{0}\|}{\|\Psi_{0}\|},$

where $\|A\|=\sup_{\|x\|_{2}=1}\|Ax\|_{2}$. When using a Bayesian method, we
consider $\widehat{\Psi}_{0}$ to be the posterior mean. For a given $(n,p)$
and a particular replicate, we obtain
$\mathcal{L}(\Psi_{0},\widehat{\Psi}_{0})$ given an estimator
$\widehat{\Psi}_{0}$ and then average this quantity over the replicates. For
the UQ experiment, we simply consider the $95\%$ entrywise credible/pseudo-
credible intervals obtained and compute the frequentist coverage over the
$R=100$ replicates. To ease the computational burden, we consider the coverage
of a randomly chosen $100\times 100$ submatrix of $\Psi_{0}$, corresponding to
the covariance of $100$ randomly chosen variables held fixed across
replicates, for each $(n,p).$ Lastly, we present a runtime analysis of the
different methods considered in the simulation studies. First, we compute the
theoretical computational complexity of the proposed CC-FABLE method. Next, we
carry out numerical experiments to compare the computational runtimes for the
different approaches.

As competitors to the proposed CC-FABLE approach, we consider the
multiplicative gamma shrinkage prior formulation (Bhattacharya and Dunson,,
2011) and the approach in Ročková and George, (2016) and denote them by MGSP
and ROTATE, respectively, across all the simulations. The approach in
Bhattacharya and Dunson, (2011) implements posterior sampling with Markov
chain Monte Carlo (MCMC), yielding posterior samples of the covariance matrix,
from which we extract the posterior mean and entrywise $95\%$ credible
intervals. In contrast, the approach of Ročková and George, (2016) only
provides a point estimate of the covariance matrix, using expectation-
maximization (EM) algorithms to approximate the posterior mean corresponding
to spike and slab prior distributions on the factor loadings. As a result, we
do not consider ROTATE for our UQ experiments. For CC-FABLE, we consider
$1000$ Monte Carlo iterates; for estimation error experiments, we use the
pseudo-posterior mean, which is analytically available. For MGSP, we consider
$3000$ MCMC iterates and discard the first $1000$ as burn-in. As described
earlier, CC-FABLE does not require MCMC for inference, instead providing
direct Monte Carlo samples. To implement MGSP, we used the infinitefactor
package (Poworoznek et al.,, 2021), while code to implement ROTATE was
obtained from http://veronikarock.com/FACTOR_ANALYSIS.zip. All the simulations
were carried out in the R programming language (R Core Team,, 2021).

###  4.2 Estimation Performance

We show the results from the simulation experiments in Table 1. For all four
different combinations of $(n,p)$, the performance of CC-FABLE is either the
best or very close to being the best. Both CC-FABLE and MGSP perform uniformly
better than the ROTATE approach across the cases. For the most challenging
case with $n=500$ and $p=5000$, CC-FABLE performs noticeably better than MGSP.
We believe the well-documented issues in mixing when implementing MCMC when
sampling from the posterior for MGSP to be a contributing factor. Furthermore,
MGSP seems to be more affected than either CC-FABLE or ROTATE when the number
of dimensions increases.

Table 1: Comparison of estimation error between CC-FABLE, MGSP, and ROTATE.  Cases | CC-FABLE | MGSP | ROTATE |
---|---|---|---|---
$n=500,p=1000$ | 0.31 | 0.34 | 0.43 |
$n=1000,p=1000$ | 0.22 | 0.21 | 0.28 |
$n=500,p=5000$ | 0.32 | 0.40 | 0.43 |
$n=1000,p=5000$ | 0.24 | 0.25 | 0.29 |

###  4.3 Frequentist Coverage

As described earlier, for each $(n,p)$, we first consider $100$ randomly
chosen variables and hold them fixed across $R=100$ replicates of the data.
Next, we obtain the average coverage and interval width of $95\%$
posterior/pseudo-posterior credible intervals corresponding to the relevant
$100\times 100$ submatrix of $\Psi_{0}$. For all the simulation instances, we
construct the entries $b_{uv}$ as described in Section 2.2, and then proceed
to set the coverage-correction factor $\rho$ to be the average of
$\\{b_{uv}:1\leq u\leq v\leq p\\}$ for coverage-correction. The results,
provided in Table 2, suggest that CC-FABLE provides superior entrywise
coverage when compared to MGSP, with comparable interval width. In all the
cases, the average coverage of entrywise intervals obtained from MGSP falls
short of the nominal value $0.95$, while CC-FABLE meets the nominal coverage.
Furthermore, the average entrywise coverage obtained from MGSP decreases when
the number of dimensions increases for a fixed sample size, while the results
of CC-FABLE are not affected.

Table 2: Comparison of average frequentist coverages and average interval widths across entrywise pseudo-credible intervals obtained from CC-FABLE and MGSP. Cases | Coverage | Coverage | Width | Width
---|---|---|---|---
$(n,p)$ | (CC-FABLE) | (MGSP) | (CC-FABLE) | (MGSP)
$n=500,p=1000$ | 0.95 | 0.84 | 0.48 | 0.48
$n=1000,p=1000$ | 0.95 | 0.72 | 0.35 | 0.31
$n=500,p=5000$ | 0.95 | 0.81 | 0.50 | 0.46
$n=1000,p=5000$ | 0.95 | 0.69 | 0.35 | 0.31

###  4.4 Computational Efficiency

For the purposes of computing order complexity, we assume $k$ is fixed and
$n,p$ vary. We first consider the order complexity when generating $N_{0}$
pseudo-posterior samples of the covariance matrix. The first stage SVD of the
$n\times p$ matrix $\mathbf{Y}$ has $\mathcal{O}(np\min(n,p))$ complexity
(Vasudevan and Ramakrishna,, 2017), which needs to be only computed once.
Computing the $p$ hyperparameters $\mu_{j}$ and $\delta_{j}^{2}$ for
$j=1,\ldots,p$ has $\mathcal{O}(np)$ complexity, while generating each Monte
Carlo sample of $({\Lambda}_{C},{\Sigma}_{C})$ for $j=1,\ldots,p$ has
$\mathcal{O}(p)$ complexity since $({\lambda}_{C,j},{\sigma}_{C,j}^{2})$ are
independently generated for $j=1,\ldots,p$. Lastly, obtaining
$\Psi_{C}=\Lambda_{C}\Lambda_{C}^{\top}+\Sigma_{C}$ has $\mathcal{O}(p^{2})$
complexity for each sample, leading to $\mathcal{O}(p^{2}N_{0})$ complexity to
obtain $N_{0}$ samples. Thus, obtaining $N_{0}$ pseudo-posterior samples of
the covariance matrix after coverage-correction has
$\mathcal{O}\\{p^{2}N_{0}+pN_{0}+np\min(n,p)\\}=\mathcal{O}\\{p^{2}N_{0}+np\min(n,p)\\}$
complexity. Next, we consider the order complexity to obtain the pseudo-
posterior mean of the covariance. As before, we would require computing the
SVD with $\mathcal{O}\\{np\min(n,p)\\}$ complexity, while computing the
hyperparameters $(\mu_{j},\delta_{j}^{2})$ for $j=1,\ldots,p$ has
$\mathcal{O}(np)$ complexity. Obtaining the matrix $G=G_{0}G_{0}^{\top}$ from
$G_{0}$ has $\mathcal{O}(p^{2})$ complexity. Thus, computing the pseudo-
posterior mean of $\Psi$ has asymptotic complexity
$\mathcal{O}\\{p^{2}+np\min(n,p)\\}$.

The runtime for both sampling and pseudo-posterior mean computation may be
substantially improved by considering computationally efficient
implementations of the SVD step such as the randomized SVD (Halko et al.,,
2011) when either $n$ or $p$ is very large. Additionally, since the CC-FABLE
algorithm employs embarrassingly parallel sampling, the computational runtime
would be further improved by parallelization of the computations across cores.

In Figure 1, we compare the three methods CC-FABLE, MGSP, and ROTATE in terms
of obtaining (pseudo) posterior samples and point estimates of the covariance
matrix. We consider the same data-generating scenario as in Section 4.1. We
fix $n=500$ and vary $p$ from $500$ to $5000$ in increments of $500$,
obtaining runtime measurements across $R=20$ independent replicates for each
choice of $p$. As before, we consider $3000$ MCMC iterations for MGSP and
discard the first $1000$ iterates as burn-in. The experiments were carried out
on an M1 MacBook Pro with 32 GB of RAM (random-access memory). We compare CC-
FABLE with MGSP when obtaining samples of the factor loadings and the error
variances, and with ROTATE when obtaining a point estimate of the covariance
matrix. In both cases, CC-FABLE is faster compared to its competitors. When
obtaining samples, CC-FABLE is about $150-200$ times faster than MGSP, with
this ratio remaining fairly constant as $p$ increases. When obtaining a point
estimator of the covariance, CC-FABLE is at least twice as fast as ROTATE.

|
---|---
Figure 1: Plots comparing the runtimes in seconds (in $\log_{10}$ scale) of
CC-FABLE, MGSP, and ROTATE approaches. The first plot compares the runtimes of
CC-FABLE and MGSP when obtaining (pseudo) posterior samples of the covariance
matrix. The second plot compares runtimes of CC-FABLE and ROTATE when
obtaining a point estimate of the covariance matrix.

##  5 Application

In this Section, we apply our method to identify associations between
expressed genes in a high-throughput sequencing immune cell dataset.
Identifying gene expression associations in varying immune cell populations
and among multiple lineages is critical to understanding the inflammatory
response to infectious and chronic diseases (Ota et al.,, 2021; Saini et al.,,
2022). This enables subsequent detection of corresponding cellular and
molecular pathways underlying pathological conditions. For this application,
we consider the GSE109125 dataset, which contains RNAseq data from 127 highly
purified immune cell lineages including adaptive and innate lymphocytes,
myeloid cells, mast cells, and neutrophils. The lineages along with the
associated differentiation cascades were profiled using the ImmGen ULI
pipeline (Yoshida et al.,, 2019).

We first log-transformed raw data counts $R_{ij}$ for the $i$th cell and $j$th
gene as $y^{*}_{ij}=\log_{2}(R_{ij}+1)$. Next, we proceed to center the data
for analysis before fitting the latent factor model. Finally, we use the
genefilter package in R (Gentleman et al.,, 2023) to filter relevant genes. In
this case, we filtered genes according to their variance and considered the
top $10\%$ of genes with the highest variances. There were no cells with
missing entries. After carrying out appropriate pre-processing of the dataset,
we obtain expression data $\mathbf{Y}\in\mathbb{R}^{n\times p}$ for $p=5300$
genes measured on $n=205$ cells.

We compare results from applying the three methods CC-FABLE, MGSP, and ROTATE
on the data. The fitted log-likelihoods using the whole data $\mathbf{Y}$ are
$-490006.4$, $-503574$, and $-728959.9$ for CC-FABLE, MGSP, and ROTATE,
respectively. This suggests that CC-FABLE fits the data the best among the
three methods. The estimated rank when using the $\mbox{JIC}(k)$ criterion as
highlighted in Section 2.3 is $\widehat{k}=30$. As a posterior predictive
check, we investigated the average coverage of entrywise predictive intervals
across all the entries of $\mathbf{Y}$. When using CC-FABLE, the average
predictive coverage is $0.97$. To further investigate model fit, we also
looked into the proportion of variance explained for each gene using CC-FABLE.
The average proportion of variance explained across variables when using CC-
FABLE is $0.83$, indicating that the factor model is a good fit.

We now focus on estimating a $100\times 100$ submatrix of the $5300\times
5300$ covariance matrix, corresponding to the $100$ genes with the highest
variability. Suppose the indices of these genes are
$\mathcal{I}_{0}\subset\\{1,2,\ldots,p\\}$. We assess estimation performance
using out-of-sample fit on a randomly chosen $50\times 100$ submatrix
$\mathbf{Y}_{T}$ of $\mathbf{Y}$, corresponding to expression data on the
$100$ genes in $\mathcal{I}_{0}$ for $50$ cells. To estimate the covariance of
genes in $\mathcal{I}_{0}$, it is typical to only use expression data for
those genes. Alternatively, we infer the covariance of these genes jointly
with additional genes with indices in the set $\mathcal{A}(p_{S})$ of length
$p_{S}>0$ such that $\mathcal{A}(p_{S})\cap\mathcal{I}_{0}=\phi$. We vary
$p_{S}\in\\{100j:j=1,\ldots,10\\}\cup\\{2000,4000\\}$ and extract the relevant
$100\times 100$ covariance $\widehat{\Psi}(p_{S})$ of the genes with indices
in $\mathcal{I}_{0}$. We then use $\widehat{\Psi}(p_{S})$ to evaluate the out-
of-sample log-likelihood (OOSLL) on the held-out test set $\mathbf{Y}_{T}$. We
carry out this analysis for $R=10$ independent replicates of the train-test
split and compare the methods CC-FABLE, MGSP, and ROTATE on the basis of
average OOSLL. The results are illustrated in Figure 2.

In Figure 2, it is clear that both CC-FABLE and MGSP exhibit superior
performance to ROTATE. Both MGSP and ROTATE exhibit an essentially
monotonically decreasing trend as extra genes are added, so including
additional genes does not improve inferences on the genes of interest in using
these methods. In sharp contrast, CC-FABLE has rapidly improved performance
initially as extra genes are added, exhibiting a clear blessing of
dimensionality. After a few hundred genes, the gain levels off, and eventually
there is a modest decline, perhaps due to the need to add additional factors
not related to the genes of interest when large numbers of additional genes
are added.

The proposed method thus shows promising results in this high-dimensional gene
expression data application and compares favourably to the two state-of-the-
art gold standard approaches to Bayesian factor analysis. The results obtained
from CC-FABLE come at a fraction of the computational budget compared to MGSP
due to its embarrassingly parallel sampling scheme requiring no MCMC. When
fitted on the full dataset on an M1 MacBook Pro with $32$ GB of RAM, CC-FABLE
only took $1.1$ seconds to compute the pseudo-posterior mean while running
MGSP with $3000$ MCMC and $1000$ burn-in iterates took $\sim 27$ minutes. This
translates to a speed-up of close to $1600$ times when using CC-FABLE.

Figure 2: Figure showing log-likelihood of the test set (y-axis) averaged over
$R=10$ replicates of the train-test split, when using additional genes
(x-axis) to estimate the covariance of the relevant set of $100$ variables,
using CC-FABLE, MGSP, and ROTATE. The training set has $155$ cells and the
test set has $50$ cells for each of the $R=10$ train-test split replicates.

##  6 Discussion

In this paper, we developed a computationally scalable approach to fit high-
dimensional Bayesian latent factor models with Gaussian data. Leveraging on a
blessing of dimensions, the proposed FABLE approach completely bypasses the
need to carry out MCMC to provide pseudo-posterior samples of the covariance
matrix for inference. By avoiding MCMC altogether, our approach averts the
pitfalls of MCMC in higher dimensions, such as slow mixing and unsatisfactory
convergence. Due to its embarrassingly parallel nature, the FABLE approach
leads to immense computational benefits over current approaches dependent on
MCMC. We developed a coverage-correction modification to FABLE called CC-FABLE
with provable guarantees on asymptotic frequentist coverage of its pseudo-
credible intervals. The CC-FABLE pseudo-posterior also enjoys desirable
theoretical properties such as consistency and compares favorably to state-of-
the-art approaches to Bayesian factor analysis from the viewpoint of
estimation error, uncertainty quantification, and computational efficiency in
our numerical experiments and gene expression data application.

This work naturally leads to several important directions for future research.
The first direction is to model non-Gaussian data, such as binary or count
data. Such data are routinely collected in ecological studies. For example,
$y_{ij}$ for $i=1,\ldots,n$ and $j=1,\ldots,p$ could indicate the
presence/absence or abundance of species $j$ in sample $i$. For
presence/absence data, let
$\mbox{P}(y_{ij}=1\mid\alpha_{0},\eta_{i},\lambda_{j})=\mbox{F}\left(\alpha_{0}+\eta_{i}^{\top}\lambda_{j}\right),$
with $\eta_{i}\overset{iid}{\sim}N_{k}(0,\mathbb{I}_{k})$. Here, $F$ is a link
function, such as the probit or the logistic link. The observations $y_{ij}$
are mutually independent conditional on the latent factors $\eta_{i}$, but
integrating them out induces dependence between $y_{i1},\ldots,y_{ip}.$
Bayesian inference in such models typically proceeds via a data augmentation
Gibbs sampler by introducing auxiliary continuous variables $z_{ij}$ and
alternating between sampling $\eta_{i},\lambda_{j},$ and $z_{ij}$ for
$i=1,\ldots,n$ and $j=1,\ldots,p$ given the data $\mathbf{Y}=(y_{ij})$ using
MCMC (Albert and Chib,, 1993; Polson et al.,, 2013). However, in our
motivating ecological settings involving insects and fungi, $p$ is on the
order of $10,000-100,000$, rendering current MCMC-based Bayesian factor
analytic approaches computationally infeasible. To extend FABLE to this
setting, we can pre-estimate the latent factors via an SVD or an appropriate
generalized linear model (GLM) variant, and then fit the resulting Bayesian
GLMs in parallel for each species, potentially relying on Laplace
approximations.

An additional important direction is to modify the FABLE inferential framework
to allow more intricate hierarchical modeling. For example in Bayesian factor
analysis, it is often of interest to (1) include covariate effects, (2) use
more elaborate priors on the loadings to favor sparsity and adaptive selection
of the number of factors, (3) consider more flexible latent factor
distributions, and/or (4) allow nonlinear measurement structure. In
exploratory analyses, we have observed good performance for ad hoc FABLE
modifications to include covariates and for simulated data having non-Gaussian
latent factor distributions. Problem (2) can potentially be addressed in the
second stage via MCMC algorithms implemented in parallel for the different
sparse regressions after inferring the latent factors using a sparse SVD in a
first stage. Problem (4) may necessitate non-linear dimensionality reduction
in the first stage. Obtaining concrete implementations with theory support for
problems (1)-(4) is of future interest.

The current work also hints at further investigation into the theoretical
properties of Bayesian factor models. Although there are existing results
focusing on posterior consistency with carefully chosen shrinkage priors on
the factor loadings, these are primarily in the case where the idiosyncratic
model errors are Gaussian. It would thus be interesting to explore how the
blessing of dimensions influences both posterior consistency and uncertainty
quantification in different settings, such as when the data are non-Gaussian.
Furthermore, there remains a clear need to extend existing theoretical results
on Bayesian factor models to results on uncertainty quantification of the
covariance matrix.

## Acknowledgments

This work was partially supported by Merck, the United States Office of Naval
Research (ONR) under grant N00014-21-1-2510, the European Research Council
(ERC) under grant 856506, and the National Institutes of Health (NIH) under
grant R01ES035625. Code used to implement CC-FABLE used in the simulation
experiments and the gene expression data application is available to download
at https://github.com/shounakch/FABLE.

## Data Availability

The gene expression data set used in Section 5 of this paper is publicly
available at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109125.

## References

  * Akaike, (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on Automatic Control, 19(6):716–723.
  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  Journal of the American Statistical Association, 88(422):669–679.
  * Avalos-Pacheco et al., (2022) Avalos-Pacheco, A., Rossell, D., and Savage, R. S. (2022).  Heterogeneous large datasets integration using Bayesian factor regression.  Bayesian Analysis, 17(1):33–66.
  * Bandeira and Van Handel, (2016) Bandeira, A. S. and Van Handel, R. (2016).  Sharp nonasymptotic bounds on the norm of random matrices with independent entries.  Annals of Probability, 44:2479–2506.
  * Bhattacharya and Dunson, (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse bayesian infinite factor models.  Biometrika, 98(2):291–306.
  * Bickel and Levina, (2008) Bickel, P. J. and Levina, E. (2008).  Regularized estimation of large covariance matrices.  The Annals of Statistics, 36(1):199–227.
  * Bien and Tibshirani, (2011) Bien, J. and Tibshirani, R. J. (2011).  Sparse estimation of a covariance matrix.  Biometrika, 98(4):807–820.
  * Bolfarine et al., (2022) Bolfarine, H., Carvalho, C. M., Lopes, H. F., and Murray, J. S. (2022).  Decoupling shrinkage and selection in Gaussian linear factor analysis.  Bayesian Analysis, 1(1):1–23.
  * Chandra et al., (2023) Chandra, N. K., Dunson, D. B., and Xu, J. (2023).  Inferring covariance structure from multiple data sources via subspace factor analysis.  arXiv preprint arXiv:2305.04113.
  * Chen and Li, (2022) Chen, Y. and Li, X. (2022).  Determining the number of factors in high-dimensional generalized latent factor models.  Biometrika, 109(3):769–782.
  * De Vito et al., (2021) De Vito, R., Bellio, R., Trippa, L., and Parmigiani, G. (2021).  Bayesian multistudy factor analysis for high-throughput biological data.  The Annals of Applied Statistics, 15(4):1723–1741.
  * Eckart and Young, (1936) Eckart, C. and Young, G. (1936).  The approximation of one matrix by another of lower rank.  Psychometrika, 1(3):211–218.
  * Fan et al., (2023) Fan, J., Lou, Z., and Yu, M. (2023).  Are latent factor regression and sparse regression adequate?  Journal of the American Statistical Association, pages 1–13.
  * Foo and Shim, (2021) Foo, Y. S. and Shim, H. (2021).  A comparison of Bayesian inference techniques for sparse factor analysis.  arXiv preprint arXiv:2112.11719.
  * Frühwirth-Schnatter, (2023) Frühwirth-Schnatter, S. (2023).  Generalized cumulative shrinkage process priors with applications to sparse Bayesian factor analysis.  Philosophical Transactions of the Royal Society A, 381(2247):20220148.
  * Gentleman et al., (2023) Gentleman, R., Carey, V., Huber, W., and Hahne, F. (2023).  Genefilter: methods for filtering genes from high-throughput experiments.  R package version, 1.84.0.
  * Halko et al., (2011) Halko, N., Martinsson, P.-G., and Tropp, J. A. (2011).  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review, 53(2):217–288.
  * Hansen et al., (2023) Hansen, B., Avalos-Pacheco, A., Russo, M., and De Vito, R. (2023).  Fast variational inference for Bayesian factor analysis in single and multi-study settings.  arXiv preprint arXiv:2305.13188.
  * Laurent and Massart, (2000) Laurent, B. and Massart, P. (2000).  Adaptive estimation of a quadratic functional by model selection.  The Annals of Statistics, pages 1302–1338.
  * Li et al., (2018) Li, Q., Cheng, G., Fan, J., and Wang, Y. (2018).  Embracing the blessing of dimensionality in factor models.  Journal of the American Statistical Association, 113(521):380–389.
  * Luo et al., (2021) Luo, Y., Han, R., and Zhang, A. R. (2021).  A Schatten-q low-rank matrix perturbation analysis via perturbation projection error bound.  Linear Algebra and its Applications, 630:225–240.
  * Ma and Liu, (2022) Ma, Y. and Liu, J. S. (2022).  On posterior consistency of Bayesian factor models in high dimensions.  Bayesian Analysis, 17(3):901–929.
  * Man and Culpepper, (2022) Man, A. X. and Culpepper, S. A. (2022).  A mode-jumping algorithm for Bayesian factor analysis.  Journal of the American Statistical Association, 117(537):277–290.
  * Ota et al., (2021) Ota, M., Nagafuchi, Y., Hatano, H., Ishigaki, K., Terao, C., Takeshima, Y., Yanaoka, H., Kobayashi, S., Okubo, M., Shirai, H., et al. (2021).  Dynamic landscape of immune cell-specific gene regulation in immune-mediated diseases.  Cell, 184(11):3006–3021.
  * 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(3):1102–1130.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using pólya–gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Poworoznek et al., (2021) Poworoznek, E., Ferrari, F., and Dunson, D. (2021).  Efficiently resolving rotational ambiguity in Bayesian matrix sampling with matching.  arXiv preprint arXiv:2107.13783.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Richard et al., (2012) Richard, E., Savalle, P.-A., and Vayatis, N. (2012).  Estimation of simultaneously sparse and low rank matrices.  arXiv preprint arXiv:1206.6474.
  * Ročková and George, (2016) Ročková, V. and George, E. I. (2016).  Fast Bayesian factor analysis via automatic rotations to sparsity.  Journal of the American Statistical Association, 111(516):1608–1622.
  * Roy et al., (2021) Roy, A., Lavine, I., Herring, A. H., and Dunson, D. B. (2021).  Perturbed factor analysis: Accounting for group differences in exposure profiles.  Annals of Applied Statistics, 15(3):1386–1404.
  * 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.
  * Saini et al., (2022) Saini, A., Ghoneim, H. E., Lio, C.-W. J., Collins, P. L., and Oltz, E. M. (2022).  Gene regulatory circuits in innate and adaptive immune cells.  Annual Review of Immunology, 40:387–411.
  * Schiavon et al., (2022) Schiavon, L., Canale, A., and Dunson, D. B. (2022).  Generalized infinite factorization models.  Biometrika, 109(3):817–835.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, pages 461–464.
  * Shikhaliev et al., (2019) Shikhaliev, A. P., Potter, L. C., and Chi, Y. (2019).  Low-rank structured covariance matrix estimation.  IEEE Signal Processing Letters, 26(5):700–704.
  * Srivastava et al., (2017) Srivastava, S., Engelhardt, B. E., and Dunson, D. B. (2017).  Expandable factor analysis.  Biometrika, 104(3):649–663.
  * Vasudevan and Ramakrishna, (2017) Vasudevan, V. and Ramakrishna, M. (2017).  A hierarchical singular value decomposition algorithm for low rank matrices.  arXiv preprint arXiv:1710.02812.
  * Vershynin, (2010) Vershynin, R. (2010).  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027.
  * Vershynin, (2011) Vershynin, R. (2011).  Spectral norm of products of random and deterministic matrices.  Probability Theory and Related Fields, 150(3-4):471–509.
  * Wang and Stephens, (2021) Wang, W. and Stephens, M. (2021).  Empirical bayes matrix factorization.  The Journal of Machine Learning Research, 22(1):5332–5371.
  * Xie et al., (2022) Xie, F., Cape, J., Priebe, C. E., and Xu, Y. (2022).  Bayesian sparse spiked covariance model with a continuous matrix shrinkage prior.  Bayesian Analysis, 17(4):1193–1217.
  * Xie and Wu, (2023) Xie, F. and Wu, D. (2023).  An eigenvector-assisted estimation framework for signal-plus-noise matrix models.  Biometrika.
  * Yoshida et al., (2019) Yoshida, H., Lareau, C. A., Ramirez, R. N., Rose, S. A., Maier, B., Wroblewska, A., Desland, F., Chudnovskiy, A., Mortha, A., Dominguez, C., et al. (2019).  The cis-regulatory atlas of the mouse immune system.  Cell, 176(4):897–912.
  * Zhang and Zhou, (2020) Zhang, A. R. and Zhou, Y. (2020).  On the non-asymptotic and sharp lower tail bounds of random variables.  Stat, 9(1):e314.
  * Zhang and Zou, (2014) Zhang, T. and Zou, H. (2014).  Sparse precision matrix estimation via lasso penalized D-trace loss.  Biometrika, 101(1):103–120.
  * Zhao et al., (2016) Zhao, S., Gao, C., Mukherjee, S., and Engelhardt, B. E. (2016).  Bayesian group factor analysis with structured sparsity.  The Journal of Machine Learning Research.

## Supplementary Material

Details on the content of this supplement are as follows.

  1. 1.

Section A contains the proof of Proposition 2.1 in the main paper.

  2. 2.

Section B contains the proof of Theorem 3.6 in the main paper.

  3. 3.

Section C contains the proof of Theorem 3.7 in the main paper.

  4. 4.

Section D contains the proof of Theorem 3.8 in the main paper.

  5. 5.

Section E contains the proof of Proposition 3.5 and relevant lemmas used to
prove Theorem 3.6 in the main paper.

  6. 6.

Section F contains the proof of relevant lemmas used to prove Theorems 3.7 and
3.8 in the main paper.

##  Appendix A Proof of Proposition 2.1

###### Proof.

  1. 1.

We first have
$\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}=\widehat{C}^{-1}\mathbf{A}^{\top}\mathbf{A}(\widehat{C}^{\top})^{-1}=\widehat{C}^{-1}(n\widehat{C}\widehat{C}^{\top})(\widehat{C}^{\top})^{-1}=n\mathbb{I}_{k}$.
Next, we have
$\widehat{\mathbf{M}}\widehat{\mathbf{M}}^{\top}=\mathbf{A}(\widehat{C}\widehat{C}^{\top})^{-1}\mathbf{A}^{\top}=(UD/\sqrt{p})(npD^{-2})(DU^{\top}/\sqrt{p})=nUU^{\top}$.

  2. 2.

For any $\widehat{\mathbf{M}}$ satisfying the condition,
$\mathbf{K}=\\{\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}+(\mathbb{I}_{k}/\tau^{2})\\}^{-1}=(n+\tau^{-2})^{-1}\mathbb{I}_{k}.$
It is clear that the pseudo-posterior of $\widetilde{\sigma}_{j}^{2}$ only
depends on $\widehat{\mathbf{M}}$ through
$\widehat{\mathbf{M}}^{\top}\widehat{\mathbf{M}}$ and
$\widehat{\mathbf{M}}\widehat{\mathbf{M}}^{\top}$ for all $1\leq j\leq p$. For
$1\leq u\leq p,$ a pseudo-posterior draw of $\tilde{\lambda}_{u}$ may be
represented as

$\widetilde{\lambda}_{u}=\mu_{u}+\widetilde{\sigma}_{u}\mathbf{K}^{1/2}Z_{u},$

where $Z_{u}\overset{ind}{\sim}N_{k}(0_{k},\mathbb{I}_{k})$ for $1\leq u\leq
p$ and $Z_{u}$ is independent of $\widetilde{\sigma}_{u}^{2}$. This implies
that for $1\leq u,v\leq p$,

$\widetilde{\lambda}_{u}^{\top}\widetilde{\lambda}_{v}=\mu_{u}^{\top}\mu_{v}+\widetilde{\sigma}_{v}\mu_{u}^{\top}\mathbf{K}^{1/2}Z_{v}+\widetilde{\sigma}_{u}\mu_{v}^{\top}\mathbf{K}^{1/2}Z_{u}+\widetilde{\sigma}_{u}\widetilde{\sigma}_{v}Z_{u}^{\top}\mathbf{K}Z_{v}.$

We now consider each expression one-by-one. The first term is
$\mu_{u}^{\top}\mu_{v}=y^{(u)\top}\widehat{\mathbf{M}}\mathbf{K}^{2}\widehat{\mathbf{M}}^{\top}y^{(v)}=(n+\tau^{-2})^{-2}y^{(u)\top}\widehat{\mathbf{M}}\widehat{\mathbf{M}}^{\top}y^{(v)}.$
When $u\neq v$, the second and third terms have Gaussian distributions with
means $0$ and variances
$(n+\tau^{-2})^{-1}\widetilde{\sigma}_{u}^{2}\|\mu_{v}\|_{2}^{2}$ and
$(n+\tau^{-2})^{-1}\widetilde{\sigma}_{v}^{2}\|\mu_{u}\|_{2}^{2}$,
respectively. When $u=v$, both terms are equal and the sum has a Gaussian
distribution with mean $0$ and variance
$4\widetilde{\sigma}_{u}^{2}\|\mu_{u}\|_{2}^{2}$. The final term is equal to
$(n+\tau^{-2})^{-1}\widetilde{\sigma}_{u}\widetilde{\sigma}_{v}Z_{u}^{\top}Z_{v}$.
Thus, all the terms depend on $\widehat{\mathbf{M}}$ only as hypothesized,
implying the result for $\widetilde{L}$. The result for $\widetilde{\Psi}$ is
obtained by simply observing that
$\widetilde{\Psi}=\widetilde{L}+\widetilde{\Sigma}.$

This proves the desired result. ∎

##  Appendix B Proof of Theorem 3.6

###  B.1 Proof of part (a)

###### Proof.

We first start with the pseudo-posterior contraction of $L_{C}$ to $L_{0}.$
Let $G_{0}=\sqrt{n}\mathbf{Y}^{\top}U/(n+\tau^{-2})$ and
$G=G_{0}G_{0}^{\top}=n\mathbf{Y}^{\top}UU^{\top}\mathbf{Y}/(n+\tau^{-2})^{2}.$
We have

$\begin{split}L_{C}-L_{0}&={\Lambda}_{C}{\Lambda}_{C}^{\top}-\Lambda_{0}\Lambda_{0}^{\top}\\\
&=\left[\frac{\mathbf{Y}^{\top}\widehat{M}\widetilde{E}^{\top}+\widetilde{E}\widehat{M}^{\top}\mathbf{Y}}{n+\tau^{-2}}+\widetilde{E}\widetilde{E}^{\top}\right]+\left(G-\Lambda_{0}\Lambda_{0}^{\top}\right),\end{split}$

by expressing $\Lambda_{C}=G_{0}+\widetilde{E}$, where
$\widetilde{E}=[\tilde{e}_{1},\ldots,\tilde{e}_{p}]^{\top}$ with
$\tilde{e}_{j}\overset{ind}{\sim}N_{k}\left(0,\frac{\rho^{2}\sigma_{C,j}^{2}}{n+\tau^{-2}}\mathbb{I}_{k}\right)$
for $j=1,\ldots,p$. We decompose

$\begin{split}&G-\Lambda_{0}\Lambda_{0}^{\top}=\frac{n}{(n+\tau^{-2})^{2}}\left(\mathbf{Y}^{\top}\mathbf{Y}-\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}\right)-\Lambda_{0}\Lambda_{0}^{\top}\\\
&=\frac{-n}{(n+\tau^{-2})^{2}}\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}+\frac{n}{(n+\tau^{-2})^{2}}\left(M_{0}\Lambda_{0}^{\top}+E\right)^{\top}\left(M_{0}\Lambda_{0}^{\top}+E\right)-\Lambda_{0}\Lambda_{0}^{\top}\\\
&=\frac{-n}{(n+\tau^{-2})^{2}}\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}+\frac{n}{(n+\tau^{-2})^{2}}\Lambda_{0}(M_{0}^{\top}M_{0}-nI)\Lambda_{0}^{\top}+\left(\frac{n^{2}}{(n+\tau^{-2})^{2}}-1\right)\Lambda_{0}\Lambda_{0}^{\top}+\\\
&\frac{n}{(n+\tau^{-2})^{2}}(\Lambda_{0}M_{0}^{\top}E+E^{\top}M_{0}\Lambda_{0}^{\top}+E^{\top}E).\end{split}$

To develop an upper bound of $\|L_{C}-L_{0}\|$, we aim to develop an upper
bound for the spectral norm of each term in the inequality above. We enumerate
them as follows.

  * (i)

First, we develop a probabilistic upper bound for ${\sigma}_{C,j}^{2}.$ Recall
that
$\|y^{(j)}\|_{2}^{2}=\sum_{i=1}^{n}y_{ij}^{2}\sim(\sigma_{0j}^{2}+\|\lambda_{0j}\|_{2}^{2})\,\chi_{n}^{2}.$
We apply the tail inequality bound of the $\chi^{2}$ distribution following
Lemma 1 in Laurent and Massart, (2000) to obtain $\|y^{(j)}\|_{2}\lesssim n$
with probability at least $1-o(1)$. Thus, with probability at least $1-o(1)$
we have

$\frac{\gamma_{0}\delta_{0}^{2}+\|y^{(j)}\|_{2}^{2}-(n+\tau^{-2})\|\mu_{j}\|_{2}^{2}}{2}\lesssim
n.$

As

${\sigma}_{C,j}^{2}\sim\mbox{IG}\left(\frac{\gamma_{0}+n}{2},\frac{\gamma_{0}\delta_{0}^{2}+\|y^{(j)}\|_{2}^{2}-(n+\tau^{-2})\|\mu_{j}\|_{2}^{2}}{2}\right),$

the Gamma distribution tail bound in Theorem 5 of Zhang and Zhou, (2020)
implies there exists a constant $C^{\prime}>0$ such that

$\begin{split}\Pi_{C}\left(\max_{1\leq j\leq p}{\sigma}_{C,j}^{2}\leq
C^{\prime}\right)=1-o_{P_{0}}(1).\end{split}$ (20)

  * (ii)

Since $\|\Lambda_{0}\|\asymp\sqrt{p_{n}}$, we have

$\left\|\left\\{\frac{n^{2}}{(n+\tau^{-2})^{2}}-1\right\\}\Lambda_{0}\Lambda_{0}^{\top}\right\|=\left|\frac{n^{2}-(n+\tau^{-2})^{2}}{(n+\tau^{-2})^{2}}\right|\|\Lambda_{0}\|^{2}\lesssim\dfrac{p_{n}}{n}.$
(21)

  * (iii)

Recall $M_{0}$ is an $n$-by-$k$ matrix with i.i.d. Gaussian entries. By Lemma
E.1, we have with probability at least $1-o(1)$ the following:

$\displaystyle\left\|\frac{n}{(n+\tau^{-2})^{2}}\Lambda_{0}(M_{0}^{\top}M_{0}-nI)\Lambda_{0}^{\top}\right\|$
$\displaystyle\lesssim\frac{1}{n}\|\Lambda_{0}\|^{2}\|M_{0}^{\top}M_{0}-nI\|$
$\displaystyle\lesssim\frac{1}{n}\|\Lambda_{0}\|^{2}\|M_{0}^{\top}M_{0}-nI\|$
$\displaystyle\lesssim\frac{p_{n}(k+\sqrt{nk})}{n}\lesssim\dfrac{p_{n}}{\sqrt{n}}.$

  * (iv)

Next, we look at

$\begin{split}&\left\|\frac{n}{(n+\tau^{-2})^{2}}(\Lambda_{0}^{\top}M_{0}E+E^{\top}M_{0}\Lambda_{0}^{\top}+E^{\top}E)\right\|\lesssim\frac{1}{n}\left(2\left\|\Lambda_{0}^{\top}M_{0}E\right\|+\left\|E^{\top}E\right\|\right).\end{split}$

By Lemma E.2, $\|E^{\top}E\|=\|E\|^{2}\lesssim n+p_{n}+\sqrt{np_{n}}$ with
probability at least $1-o(1)$ and
$\|\Lambda_{0}M_{0}^{\top}E\|\leq\|\Lambda_{0}\|\,\|M_{0}\|\,\|E\|\lesssim\sqrt{np_{n}}(\sqrt{n}+\sqrt{p_{n}})$
with probability at least $1-o(1)$. Thus, with probability at least $1-o(1)$,

$\left\|\frac{n}{(n+\tau^{-2})^{2}}(\Lambda_{0}^{\top}M_{0}E+E^{\top}M_{0}\Lambda_{0}^{\top}+E^{\top}E)\right\|\lesssim\frac{n+p_{n}+p_{n}\sqrt{n}+n\sqrt{p_{n}}+\sqrt{np_{n}}}{n}\lesssim\dfrac{p_{n}\sqrt{n}+n\sqrt{p_{n}}}{n}.$

  * (v)

Since $U$ corresponds to the top $k$ left singular vectors of $\mathbf{Y}$, we
have
$\|\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}\|=\|\mathbf{Y}^{\top}U_{\perp}\|^{2}=s_{k+1}^{2}(\mathbf{Y})$,
where $s_{k+1}^{2}(\mathbf{Y})$ is the $(k+1)$-th singular value of
$\mathbf{Y}$. By Eckart–Young Theorem (Lemma E.3), we have

$\|\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}\|=s_{k+1}^{2}(\mathbf{Y})=\min_{B:\text{rank}(B)\leq
k}\|\mathbf{Y}-B\|^{2}\leq\|\mathbf{Y}-M_{0}\Lambda_{0}^{\top}\|^{2}=\|E\|^{2}\lesssim
n+p_{n}+\sqrt{np_{n}}.$

Then,

$\left\|\frac{-n}{(n+\tau^{-2})^{2}}\mathbf{Y}^{\top}U_{\perp}U_{\perp}^{\top}\mathbf{Y}\right\|\lesssim\frac{n+p_{n}+\sqrt{np_{n}}}{n}.$

  * (vi)

Finally, we consider

$\begin{split}\left\|\frac{\mathbf{Y}^{\top}\widehat{M}\widetilde{E}^{\top}+\widetilde{E}\widehat{M}^{\top}\mathbf{Y}}{n+\tau^{-2}}+\widetilde{E}\widetilde{E}^{\top}\right\|\lesssim\dfrac{2}{n}\left\|\mathbf{Y}^{\top}\widehat{M}\widetilde{E}^{\top}\right\|+\|\widetilde{E}\widetilde{E}^{\top}\|.\end{split}$

By Lemma E.2, the following hold with probability at least $1-o_{P_{0}}(1)$:

$\|\widetilde{E}\widetilde{E}^{\top}\|=\|\widetilde{E}\|^{2}\lesssim\frac{p_{n}}{n}\rho^{2}\left(\max_{1\leq
j\leq p}{\sigma}_{C,j}^{2}\right)\lesssim\frac{p_{n}}{n}.$
$\displaystyle\frac{1}{n}\left\|\mathbf{Y}^{\top}\widehat{M}\widetilde{E}^{\top}\right\|$
$\displaystyle\lesssim\frac{1}{n}\|\mathbf{Y}^{\top}\sqrt{n}U\|\left\|\widetilde{E}\right\|\lesssim\frac{\rho}{\sqrt{n}}\|\mathbf{Y}\|\sqrt{\frac{p_{n}}{n}\left(\max_{1\leq
j\leq p}{\sigma}^{2}_{C,j}\right)}$
$\displaystyle\lesssim\left(\frac{\sqrt{np_{n}}+\sqrt{n}+\sqrt{p_{n}}}{\sqrt{n}}\right)\sqrt{\frac{p_{n}}{n}}\lesssim\frac{p_{n}}{\sqrt{n}}.$

Thus, with probability at least $1-o_{P_{0}}(1)$, we have

$\left\|\frac{\mathbf{Y}^{\top}\widehat{M}\widetilde{E}^{\top}+\widetilde{E}\widehat{M}^{\top}\mathbf{Y}}{(n+\tau^{-2})}+\widetilde{E}\widetilde{E}^{\top}\right\|\lesssim\dfrac{p_{n}}{\sqrt{n}}.$

By combining the previous steps (i)-(vi) and
$\|\Lambda_{0}\|\asymp\sqrt{p_{n}}$, we have proved the desired result. ∎

###  B.2 Proof of part (b)

###### Proof.

We now show the contraction result for the pseudo-posterior of $\Sigma_{C}$.
We start out by observing that $||{\Sigma}_{C}-\Sigma_{0}||=\max_{1\leq j\leq
p}|D_{j}|$, where $D_{j}={\sigma}_{C,j}^{2}-\sigma_{0j}^{2}$ for
$j=1,\ldots,p.$ Let $\kappa_{j}={\sigma}_{C,j}^{-2}$, so that
$\kappa_{j}\overset{ind}{\sim}G(\gamma_{n}/2,\gamma_{n}\delta_{j}^{2}/2)$
under $\Pi_{C}$. Let

$U_{j}=\dfrac{\gamma_{n}\delta_{j}^{2}}{2}\kappa_{j}-\dfrac{\gamma_{n}}{2}.$

We can now express $D_{j}$ as

$D_{j}=\left(1+\dfrac{2U_{j}}{\gamma_{n}}\right)^{-1}\,\left\\{(\delta_{j}^{2}-\sigma_{0j}^{2})-\dfrac{2U_{j}}{\gamma_{n}}\sigma_{0j}^{2}\right\\}.$

Using Lemma E.7, we have $\max_{1\leq j\leq
p}|U_{j}|/\gamma_{n}\lesssim\\{(\log p_{n})/n\\}^{1/3}$ with probability at
least $1-o_{P_{0}}(1).$ Thus, $\min_{1\leq j\leq
p}|1+(2U_{j}/\gamma_{n})|\gtrsim 1/2$ with probability at least
$1-o_{P_{0}}(1).$ Therefore with probability at least $1-o_{P_{0}}(1),$

$\max_{1\leq j\leq p}|D_{j}|\lesssim 2\max_{1\leq j\leq
p}|\delta_{j}^{2}-\sigma_{0j}^{2}|+4\sigma_{0j}^{2}\max_{1\leq j\leq
p}\dfrac{|U_{j}|}{\gamma_{n}}.$

Let $Q_{j}=\delta_{j}^{2}-\sigma_{0j}^{2}$. Using Lemma E.6, we can represent

$Q_{j}=\dfrac{\sigma_{0j}^{2}}{n}\left\\{\dfrac{Z_{j}}{\sigma_{0j}^{2}}-(n-k)\right\\}-\dfrac{k\sigma_{0j}^{2}}{n}+F_{j},$

where $Z_{j}/\sigma_{0j}^{2}\sim\chi_{n-k}^{2}\equiv G\\{(n-k)/2,1/2\\}$ and
$\max_{1\leq j\leq p}|F_{j}|\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}}$ with
probability at least $1-o(1).$ Using Lemma E.7, we obtain

$\max_{1\leq j\leq p}\left|\dfrac{Z_{j}}{\sigma_{0j}^{2}}-(n-k)\right|\lesssim
n\left(\dfrac{\log p_{n}}{n}\right)^{1/3}$

with probability at least $1-o(1)$. Thus, with probability at least $1-o(1)$,
we get

$\max_{1\leq j\leq p}|Q_{j}|\lesssim\left(\dfrac{\log
p_{n}}{n}\right)^{1/3}+\dfrac{1}{p_{n}}.$

From Lemma E.7, we have $\max_{1\leq j\leq
p}|U_{j}|/\gamma_{n}\lesssim\left({\log p_{n}}/{n}\right)^{1/3}$ with
probability at least $1-o_{P_{0}}(1)$. Thus, we obtain

$\max_{1\leq j\leq p}|D_{j}|\lesssim\left(\dfrac{\log
p_{n}}{n}\right)^{1/3}+\dfrac{1}{p_{n}}$

with probability at least $1-o_{P_{0}}(1),$ for all $j=1,\ldots,p$ or
equivalently,

$\|{\Sigma}_{C}-\Sigma_{0}\|\lesssim\left(\dfrac{\log
p_{n}}{n}\right)^{1/3}+\dfrac{1}{p_{n}}$

with probability at least $1-o_{P_{0}}(1)$. This proves the result. ∎

###  B.3 Proof of part (c)

###### Proof.

Under the conditions of Theorem 1, we have $\|\Psi_{0}\|\asymp\|L_{0}\|\asymp
p_{n}.$ By the triangle inequality,
$\|\Psi_{C}-\Psi_{0}\|\leq\|L_{C}-L_{0}\|+\|\Sigma_{C}-\Sigma_{0}\|$. Thus, we
have under $\Pi_{C}$ that

$\dfrac{\|{\Psi}_{C}-\Psi_{0}\|}{\|\Psi_{0}\|}\lesssim\dfrac{1}{\sqrt{p_{n}}}+\dfrac{1}{\sqrt{n}}+o\left(\dfrac{1}{p_{n}}\right)\lesssim\dfrac{1}{\sqrt{p_{n}}}+\dfrac{1}{\sqrt{n}}$

with probability at least $1-o_{P_{0}}(1)$, proving the result. ∎

##  Appendix C Proof of Theorem 3.7

###### Proof.

First let $u\neq v$, so that $\Psi_{C,uv}=L_{C,uv}$, where $L_{C,uv}$ is the
$uv$th element of $L_{C}=\Lambda_{C}\Lambda_{C}^{\top}$. According to the
pseudo-posterior generation mechanism after coverage correction,

$L_{C,uv}=\mu_{u}^{\top}\mu_{v}+\dfrac{\rho\left({\sigma}_{C,v}\mu_{u}^{\top}\tilde{e}_{C,v}+{\sigma}_{C,u}\mu_{v}^{\top}\tilde{e}_{C,u}\right)}{\sqrt{n+\tau^{-2}}}+\dfrac{\rho^{2}{\sigma}_{C,u}{\sigma}_{C,v}\tilde{e}_{C,u}^{\top}\tilde{e}_{C,v}}{n+\tau^{-2}},$

for
$\widetilde{e}_{C,u},\widetilde{e}_{C,v}\overset{ind}{\sim}N_{k}(0,\mathbb{I}_{k})$.
We can express

$\displaystyle\rho\left(\sigma_{C,v}\mu_{u}^{\top}\tilde{e}_{C,v}+\sigma_{C,u}\mu_{v}^{\top}\tilde{e}_{C,u}\right)$
$\displaystyle=l_{0,uv}(\rho)R_{uv}+d_{uv}(\rho)R_{uv},$

for $R_{uv}\sim N(0,1)$ such that $R_{uv}$ is independent of
$\mathbf{Y},\sigma_{C,u}^{2},$ and $\sigma_{C,v}^{2}$, with
$l_{0,uv}^{2}(\rho)=\rho^{2}\left(\sigma_{0v}^{2}\|\lambda_{0u}\|_{2}^{2}+\sigma_{0u}^{2}\|\lambda_{0v}\|_{2}^{2}\right)$
and
$d_{uv}(\rho)=\rho(\sigma_{C,v}^{2}||\mu_{u}||_{2}^{2}+\sigma_{C,u}^{2}||\mu_{v}||_{2}^{2})^{1/2}-l_{0,uv}(\rho)$.
We can now express

$\displaystyle\sqrt{n}(L_{C,uv}-\mu_{u}^{\top}\mu_{v})=\dfrac{\sqrt{n}}{\sqrt{n+\tau^{-2}}}l_{0,uv}(\rho)R_{uv}+d_{uv}^{*}(\rho),$

where

$d_{uv}^{*}(\rho)=\sqrt{\dfrac{n}{n+\tau^{-2}}}d_{uv}(\rho)R_{uv}+\sqrt{n}\left(\dfrac{\rho^{2}\sigma_{C,v}\sigma_{C,u}\tilde{e}_{C,u}^{\top}\tilde{e}_{C,v}}{n+\tau^{-2}}\right).$

Due to pseudo-posterior concentration of $\sigma_{C,j}^{2}$ to
$\sigma_{0j}^{2}$ and convergence in probability of $||\mu_{j}||_{2}^{2}$ to
$||\lambda_{0j}||_{2}^{2}$, we have pseudo-posterior concentration of
$d_{uv}(\rho)$ around $0$ as $n\to\infty$ for any finite $\rho>0$.
Furthermore, with probability at least $1-o_{P_{0}}(1)$,

$\sqrt{n}\left(\dfrac{{\sigma}_{C,u}{\sigma}_{C,v}\tilde{e}_{C,u}^{\top}\tilde{e}_{C,v}}{n+\tau^{-2}}\right)\lesssim\dfrac{1}{\sqrt{n}}.$

Thus, the pseudo-posterior of $d_{uv}^{*}(\rho)$ concentrates around $0$
asymptotically. As $n\to\infty$, Lemma F.2 leads to

$\underset{x}{\sup}\,\left|{\Pi}_{C}\left\\{\dfrac{\sqrt{n}(L_{C,uv}-\mu_{u}^{\top}\mu_{v})}{l_{0,uv}(\rho)}\leq
x\right\\}-\Phi(x)\right|\overset{P_{0}}{\to}0.$

For $u=v$, we first observe that
$\sqrt{n}(\Psi_{C,uu}-T_{uu})=\sqrt{n}(L_{C,uu}-\|\mu_{u}\|_{2}^{2})+\sqrt{n}(\sigma_{C,u}^{2}-\delta_{u}^{2}).$
For a given $1\leq j\leq p$, using the decomposition in part (b) of the proof
of Theorem 3.6, we can express

$\sqrt{n}(\sigma_{C,j}^{2}-\delta_{j}^{2})=-\dfrac{2U_{j}}{\gamma_{n}}\sigma_{0j}^{2}+\mathcal{G}_{j},$

for all $1\leq j\leq p$, where
$|\mathcal{G}_{j}|\lesssim(1/\sqrt{n})+(\sqrt{n}/p_{n})$ with probability at
least $1-o_{P_{0}}(1).$ Letting
$\mathcal{N}_{j}=-\dfrac{2U_{j}}{\gamma_{n}}\sigma_{0j}^{2},$ we have
$\mathcal{N}_{j}/\sqrt{2\sigma_{0u}^{4}}\implies N(0,1)$, using a normal
approximation to the Gamma distribution. Thus, we can write

$\sqrt{n}\left(\Psi_{C,uu}-T_{uu}\right)=\left(2\rho\,\sigma_{0u}\|\lambda_{0u}\|_{2}\right)R_{uu}+\mathcal{N}_{u}+\mathcal{G}_{u}.$

By construction, $R_{uu}$ is independent of $\mathcal{N}_{u}$ for $1\leq u\leq
p$. Combining all these observations and using Lemmas F.1 and F.2, we have the
desired result for $u=v$, with
$l_{0,uu}(\rho)=2\sigma_{0u}^{4}+4\rho^{2}\sigma_{0u}^{2}\|\lambda_{0u}\|_{2}^{2}$
for $1\leq u\leq p$. ∎

##  Appendix D Proof of Theorem 3.8

###### Proof.

We first assume $u\neq v$. Let

$S_{uv}=\dfrac{n}{(n+\tau^{-2})^{2}}y^{(u)\top}UU^{\top}y^{(v)}-\lambda_{0u}^{\top}\lambda_{0v}$

and $f_{n}=n/(n+\tau^{-2})^{2}.$ Since
$y^{(u)}=M_{0}\lambda_{0u}+\epsilon^{(u)}$, we can decompose $S_{uv}$ as

$\displaystyle S_{uv}=$
$\displaystyle\lambda_{0u}^{\top}(f_{n}M_{0}^{\top}UU^{\top}M_{0}-\mathbb{I}_{k})\lambda_{0v}$
$\displaystyle+f_{n}(\lambda_{0u}^{\top}M_{0}^{\top}UU^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}UU^{\top}\epsilon^{(u)}+\epsilon^{(u)\top}UU^{\top}\epsilon^{(v)}).$

Observe that $f_{n}=(1/n)+\Delta_{n}$, where $|\Delta_{n}|\asymp 1/n^{2}$. We
first break up $S_{uv}$ into two parts $S_{uv}=D_{uv}+R_{uv},$ where

$\displaystyle D_{uv}$
$\displaystyle=\lambda_{0u}^{\top}\left(\dfrac{1}{n}M_{0}^{\top}U_{0}U_{0}^{\top}M_{0}-\mathbb{I}_{k}\right)\lambda_{0v}+\dfrac{1}{n}(\lambda_{0u}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(u)}),$
$\displaystyle R_{uv}=$
$\displaystyle\Delta_{n}\lambda_{0u}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}M_{0}\lambda_{0v}+f_{n}\lambda_{0u}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})M_{0}\lambda_{0v}$
$\displaystyle+\Delta_{n}(\lambda_{0u}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(u)})$
$\displaystyle+f_{n}(\lambda_{0u}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})\epsilon^{(u)})$
$\displaystyle+f_{n}\epsilon^{(u)\top}UU^{\top}\epsilon^{(v)}.$

We first deal with $R_{uv}$ and show that $\sqrt{n}R_{uv}=o_{P_{0}}(1)$. We
use throughout that $\|M_{0}\|\asymp\sqrt{n}$ and
$\|\epsilon^{(u)}\|_{2}\lesssim\sqrt{n}$ with probability at least $1-o(1),$
using Lemma E.8. First consider bounding
$\sqrt{n}f_{n}\epsilon^{(u)\top}UU^{\top}\epsilon^{(v)}$ as follows:

$\displaystyle\left|\sqrt{n}f_{n}\epsilon^{(u)\top}UU^{\top}\epsilon^{(v)}\right|$
$\displaystyle\leq\sqrt{n}f_{n}\left|\epsilon^{(u)\top}U_{0}U_{0}^{\top}\epsilon^{(v)}\right|+\sqrt{n}f_{n}\left|\epsilon^{(u)\top}(UU^{\top}-U_{0}U_{0}^{\top})\,\epsilon^{(v)}\right|.$

The first term is
$\sqrt{n}f_{n}\left|\epsilon^{(u)}U_{0}U_{0}^{\top}\epsilon^{(v)}\right|\,=O_{P_{0}}(\sqrt{n}f_{n})=O_{P_{0}}(1/\sqrt{n})$,
since $U_{0}^{\top}\epsilon^{(u)}\sim N(0,\sigma_{0u}^{2}\mathbb{I}_{k})$ and
$k$ is finite. The second term is handled by observing that

$\displaystyle\sqrt{n}f_{n}\left|\epsilon^{(u)\top}(UU^{\top}-U_{0}U_{0}^{\top})\epsilon^{(v)}\right|$
$\displaystyle\leq\sqrt{n}f_{n}\|UU^{\top}-U_{0}U_{0}^{\top}\|\|\epsilon^{(u)}\|_{2}\|\epsilon^{(v)}\|_{2}$
$\displaystyle\lesssim
n\left(\dfrac{1}{n}+\dfrac{1}{p_{n}}\right)\sqrt{n}f_{n}$
$\displaystyle=o_{P_{0}}(1).$

Thus,
$\sqrt{n}\left|f_{n}\epsilon^{(u)\top}UU^{\top}\epsilon^{(v)}\right|=o_{P_{0}}(1)$.
We now consider the remaining terms.

  * (i)

We first observe that with probability at least $1-o(1)$,

$\displaystyle\sqrt{n}\Delta_{n}\left|\lambda_{0u}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}M_{0}\lambda_{0v}\right|$
$\displaystyle=\sqrt{n}\Delta_{n}\left|\lambda_{0u}M_{0}^{\top}M_{0}\lambda_{0v}\right|$
$\displaystyle\leq\sqrt{n}\Delta_{n}\|M_{0}^{\top}M_{0}\|||\lambda_{0u}||_{2}||\lambda_{0v}||_{2}\lesssim\dfrac{1}{\sqrt{n}}.$

  * (ii)

Using Proposition 3.5, we next consider

$\displaystyle\sqrt{n}f_{n}\left|\lambda_{0u}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})M_{0}\lambda_{0v}\right|$
$\displaystyle\lesssim\sqrt{n}f_{n}\|M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})M_{0}\|$
$\displaystyle\leq\sqrt{n}f_{n}\|M_{0}\|^{2}\|UU^{\top}-U_{0}U_{0}^{\top}\|$
$\displaystyle=o_{P_{0}}(1).$

  * (iii)

We next observe that with probability at least $1-o(1),$

$\displaystyle\sqrt{n}\Delta_{n}\left|\lambda_{0u}^{\top}U_{0}U_{0}^{\top}M_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}U_{0}U_{0}^{\top}M_{0}^{\top}\epsilon^{(u)}\right|$
$\displaystyle=\sqrt{n}\Delta_{n}\left|\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}\epsilon^{(u)}\right|$
$\displaystyle\lesssim\dfrac{2}{n^{2}}(\sqrt{n})^{3}\lesssim\dfrac{1}{\sqrt{n}}.$

  * (iv)

Using Proposition 3.5, we have lastly

$\displaystyle\sqrt{n}f_{n}\left|\lambda_{0u}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})\epsilon^{(u)}\right|$
$\displaystyle\lesssim$
$\displaystyle\,\dfrac{2}{\sqrt{n}}\|M_{0}^{\top}\|\|UU^{\top}-U_{0}U_{0}^{\top}\|\,\|\epsilon^{(v)}\|_{2}$
$\displaystyle=$ $\displaystyle\,o_{P_{0}}(1).$

Putting everything together implies $\sqrt{n}|R_{uv}|=o_{P_{0}}(1)$ for all
$1\leq u,v\leq p$. We now look at $\sqrt{n}D_{uv}$. By the central limit
theorem, we first observe that

$\displaystyle\sqrt{n}\lambda_{0u}^{\top}\left(\dfrac{1}{n}M_{0}^{\top}U_{0}U_{0}^{\top}M_{0}-\mathbb{I}_{k}\right)\lambda_{0v}$
$\displaystyle=\sqrt{n}\lambda_{0u}^{\top}\left(\dfrac{1}{n}M_{0}^{\top}M_{0}-\mathbb{I}_{k}\right)\lambda_{0v}$
$\displaystyle=\sqrt{n}(\overline{V}_{n}-\lambda_{0u}^{\top}\lambda_{0v})$
$\displaystyle{\implies}N(0,\xi_{0,uv}^{2}),$

as $n\to\infty$, where $\overline{V}_{n}=(1/n)\sum_{i=1}^{n}V_{i},$ with
$V_{i}:=(\lambda_{0u}^{\top}\eta_{0i})(\lambda_{0v}^{\top}\eta_{0i})$ for
$i=1,\ldots,n$, and

$\displaystyle\begin{split}\mathcal{\xi}_{0,uv}^{2}=\mbox{var}\left[(\lambda_{0u}^{\top}\eta_{01})(\lambda_{0v}^{\top}\eta_{01})\right]&=\left\\{\begin{array}[]{ll}(\lambda_{0u}^{\top}\lambda_{0v})^{2}+||\lambda_{0u}||_{2}^{2}||\lambda_{0v}||_{2}^{2},&\text{for
$u\neq v$},\\\ 2\|\lambda_{0u}\|_{2}^{4},&\text{for
$u=v$.}\end{array}\right.\end{split}$

For $u=v$, we use the fact that $V_{i}/\|\lambda_{0u}\|_{2}^{2}$ are
independent for $i=1,\ldots,n$ and marginally $\chi_{1}^{2}$ and use the
central limit theorem accordingly. We next observe that for $u\neq v,$

$\sqrt{n}\left\\{\dfrac{1}{n}\left(\lambda_{0u}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}U_{0}U_{0}^{\top}\epsilon^{(u)}\right)\right\\}=\sqrt{n}\left\\{\dfrac{1}{n}\left(\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}\epsilon^{(u)}\right)\right\\}.$

Let

$l_{uv}^{2}(M_{0})=\sigma_{0v}^{2}\lambda_{0u}^{\top}\dfrac{M_{0}^{\top}M_{0}}{n}\lambda_{0u}+\sigma_{0u}^{2}\lambda_{0v}^{\top}\dfrac{M_{0}^{\top}M_{0}}{n}\lambda_{0v}.$

Since $\epsilon^{(u)}$ and $\epsilon^{(v)}$ are independent for $u\neq v$, we
have

$\left.\dfrac{\sqrt{n}}{n}\left(\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}\epsilon^{(u)}\right)\,\right|M_{0}\sim
N(0,l_{uv}^{2}(M_{0}))\overset{d}{=}l_{uv}(M_{0})Z_{uv},$

where $Z_{uv}\sim N(0,1)$ and $Z_{uv}$ is independent of $M_{0}.$ Let

$l_{0,uv}^{2}=\sigma_{0v}^{2}||\lambda_{0u}||_{2}^{2}+\sigma_{0u}^{2}||\lambda_{0v}||_{2}^{2}.$

Then,

$\displaystyle\sqrt{n}D_{uv}$
$\displaystyle=\sqrt{n}\lambda_{0u}^{\top}\left(\dfrac{1}{n}M_{0}^{\top}U_{0}U_{0}^{\top}M_{0}-\mathbb{I}_{k}\right)\lambda_{0v}+l_{uv}(M_{0})Z_{uv}$
$\displaystyle=\sqrt{n}\lambda_{0u}^{\top}\left(\dfrac{1}{n}M_{0}^{\top}M_{0}-\mathbb{I}_{k}\right)\lambda_{0v}+l_{0,uv}Z_{uv}+\\{l_{uv}(M_{0})-l_{0,uv}\\}Z_{uv}.$

Since $\|M_{0}\|\asymp\sqrt{n}$, we have
$\left|l_{uv}(M_{0})+l_{0,uv}\right|=O_{P_{0}}(1)$ and
$\left|l_{uv}(M_{0})+l_{0,uv}\right|\geq l_{0,uv}\geq 2c_{1}c_{2}>0$, from
Assumptions 3.2 and 3.3. Thus, the third term can be handled by observing

$\displaystyle\\{l_{uv}(M_{0})-l_{0,uv}\\}Z_{uv}$
$\displaystyle=\dfrac{l_{uv}^{2}(M_{0})-l_{0,uv}^{2}}{l_{uv}(M_{0})+l_{0,uv}}Z_{uv}$
$\displaystyle=O_{P_{0}}\left(\dfrac{1}{\sqrt{n}}\right),$

as
$\left\|\dfrac{M_{0}^{\top}M_{0}}{n}-\mathbb{I}_{k}\right\|=O_{P_{0}}\left(\dfrac{1}{\sqrt{n}}\right)$
from Lemma E.1. Since $Z_{uv}$ is independent of $M_{0}$, Lemma F.1
immediately implies that

$\sqrt{n}D_{uv}\overset{d}{=}N(0,\xi_{0,uv}^{2}+l_{0,uv}^{2})+o_{P_{0}}(1).$

Let $\mathcal{S}_{0,uv}^{2}=l_{0,uv}^{2}+\xi_{0,uv}^{2}$. Putting all the
previous results back together and invoking Slutsky’s theorem, we obtain

$\sqrt{n}S_{uv}\overset{}{\implies}N(0,\mathcal{S}_{0,uv}^{2}).$

When $u=v$, we first observe that $\delta_{u}^{2}$ is approximately
independent of $\|\mu_{u}\|_{2}^{2}$ by replacing $UU^{\top}$ by
$U_{0}U_{0}^{\top}$, as
$U_{0}U_{0}^{\top}(\mathbb{I}_{n}-U_{0}U_{0}^{\top})=\mathbf{O}_{n\times n}$,
where $\mathbf{O}_{n\times n}$ denotes the $n\times n$ matrix with all entries
equal to zero. Finally, we use Proposition 3.5 as earlier to argue the
remainder terms are small. Using the central limit theorem for the
$\chi_{n-k}^{2}$ distribution, we obtain the desired result. ∎

##  Appendix E Proofs of Proposition 3.5 and Related Lemmas for Theorem 3.6

###  E.1 Proof of Proposition 3.5

###### Proof.

Suppose $X_{0}=U_{0}D_{0}V_{0}^{\top}=M_{0}\Lambda_{0}^{\top}$ is the true
signal. We start with
$\|X_{0}^{\top}UU^{\top}X_{0}-X_{0}^{\top}U_{0}U_{0}^{\top}X_{0}\|$ and
observe that

$\displaystyle\|X_{0}^{\top}UU^{\top}X_{0}-X_{0}^{\top}U_{0}U_{0}^{\top}X_{0}\|$
$\displaystyle=\|\Lambda_{0}M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})M_{0}\Lambda_{0}^{\top}\|$
$\displaystyle\geq
s_{\min}^{2}(\Lambda_{0})\|M_{0}^{\top}(UU^{\top}-U_{0}U_{0}^{\top})M_{0}\|$
$\displaystyle\geq
s_{\min}^{2}(\Lambda_{0})s_{\min}^{2}(M_{0})\|UU^{\top}-U_{0}U_{0}^{\top}\|,$

since $\|AB\|\geq\|A\|s_{\min}(B)$ for any matrices $A,B$, with $s_{min}(B)$
denoting the smallest singular value of $B$. Both $\Lambda_{0}$ and $M_{0}$
have exactly $k$ positive singular values since the rank is $k$. Thus, we have

$\|X_{0}^{\top}UU^{\top}X_{0}-X_{0}^{\top}U_{0}U_{0}^{\top}X_{0}\|\geq
s_{\min}^{2}(\Lambda_{0})s_{\min}^{2}(M_{0})\|UU^{\top}-U_{0}U_{0}^{\top}\|.$

Since the Schatten $\infty$-norm is identical to the spectral norm, Theorem 2
in Luo et al., (2021) implies

$\|X_{0}^{\top}UU^{\top}X_{0}-X_{0}^{\top}U_{0}U_{0}^{\top}X_{0}\|=\|X_{0}^{\top}UU^{\top}X_{0}-X_{0}^{\top}X_{0}\|=\|(\mathbb{I}_{n}-UU^{\top})X_{0}\|^{2}\leq
4\|E\|^{2}.$

Thus with probability at least $1-o(1)$,

$\|UU^{\top}-U_{0}U_{0}^{\top}\|\,s_{\min}^{2}(\Lambda_{0})s_{\min}^{2}(M_{0})\leq{4\|E\|^{2}}.$

Using Lemma E.8, we have $s_{\min}(M_{0})\asymp\|M_{0}\|\ \asymp\sqrt{n}$ with
probability at least $1-o(1)$, since $k=o(n)$. Under Assumption 3.2,
$s_{\min}(\Lambda_{0})\asymp\|\Lambda_{0}\|\asymp\sqrt{p_{n}}$. Suppose
$\sigma_{\text{sum}}^{2}=\sum_{j=1}^{p}\sigma_{0j}^{2}$ and
$\sigma_{\max}^{2}=\max_{1\leq j\leq p}\sigma_{0j}^{2}=\mathcal{O}(1)$ by
Assumption 3.3. Lemma E.2 implies
$\|E\|\lesssim(\sigma_{\text{sum}}+\sqrt{n}\sigma_{\max})$ with probability at
least $1-o(1)$. Since $\sigma_{\text{sum}}\leq\sqrt{p_{n}}\sigma_{\max}$, one
has $\|E\|\lesssim\sqrt{n}+\sqrt{p_{n}}$ with probability at least $1-o(1)$,
as $n\to\infty$. Thus, with probability at least $1-o(1)$,

$\displaystyle\|UU^{\top}-U_{0}U_{0}^{\top}\|$
$\displaystyle\lesssim\dfrac{(\sqrt{n}+\sqrt{p_{n}})^{2}}{np_{n}}\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}},$

which proves the desired result. ∎

###  E.2 Other Lemmas for Theorem 1

######  Lemma E.1.

Let $E_{0}\in\mathbb{R}^{n_{1}\times n_{2}}$ have i.i.d. standard Gaussian
entries. Then for every $t>0$, one has

$\mathbb{P}\left(\|E_{0}\|\leq\sqrt{n_{1}}+\sqrt{n_{2}}+t\right)\geq
1-2\exp(-t^{2}/2).$
$\mathbb{E}\|E_{0}^{\top}E_{0}-n_{1}\mathbb{I}_{n_{2}}\|\leq\left(n_{1}^{1/2}+n_{2}^{1/2}\right)^{2}-n_{1}.$

###### Proof.

Refer to Corollary 5.35 in Vershynin, (2010). ∎

######  Lemma E.2.

Suppose $E_{0}\in\mathbb{R}^{n_{1}\times n_{2}}$ has independent entries such
that $E_{0,ij}=g_{ij}b_{ij}$ for $g_{ij}\overset{iid}{\sim}N(0,1)$ and
$\\{b_{ij}:1\leq i\leq n_{1},1\leq j\leq n_{2}\\}$ are fixed scalars. Let
$\sigma_{1}=\underset{i}{\max}\sqrt{\sum_{j}b_{ij}^{2}}$,
$\sigma_{2}=\underset{j}{\max}\sqrt{\sum_{i}b_{ij}^{2}},$ and
$\sigma_{*}=\underset{i,j}{\max}|b_{ij}|.$ Then for every $\epsilon\in(0,1/2]$
there exists a $c_{\epsilon}^{{}^{\prime}}>0$ such that for all $t\geq 0$,

$\mathbb{P}\left\\{\|E_{0}\|\geq(1+\epsilon)(\sigma_{1}+\sigma_{2})+t\right\\}\leq(n_{1}\wedge
n_{2})\exp\\{-t^{2}/(c_{\epsilon}^{{}^{\prime}}\sigma_{*}^{2})\\}.$

###### Proof.

Refer to Corollary 3.11 in Bandeira and Van Handel, (2016). ∎

######  Lemma E.3.

For any matrix $Z\in\mathbb{R}^{n_{1}\times n_{2}}$, let its $l$th largest
singular value be $s_{l}(Z)$ for $l=1,\ldots,n\wedge p$. Then for all $1\leq
l\leq n\wedge p,$

$s_{l+1}(Z)=\min_{\text{rank}(M)\leq l}\|Z-M\|.$

###### Proof.

Refer to Eckart and Young, (1936). ∎

######  Lemma E.4.

Suppose Assumptions 3.1-3.4 hold. For all $1\leq u,v\leq p,$ we have as
$n\to\infty$,

$\dfrac{1}{n}y^{(u)\top}UU^{\top}y^{(v)}\overset{P_{0}}{\to}\lambda_{0u}^{\top}\lambda_{0v}$

and

$\dfrac{1}{n}\|(\mathbb{I}_{n}-UU^{\top})y^{(u)}\|_{2}^{2}\overset{P_{0}}{\to}\sigma_{0u}^{2}.$

###### Proof.

We first note the identity

$\displaystyle\dfrac{1}{n}y^{(u)\top}UU^{\top}y^{(v)}-\lambda_{0u}^{\top}\lambda_{0v}$
$\displaystyle=$
$\displaystyle\,\lambda_{0u}^{\top}\left(\dfrac{M_{0}^{\top}M_{0}}{n}-\mathbb{I}_{k}\right)\lambda_{0v}+\dfrac{1}{n}(\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}+\lambda_{0v}^{\top}M_{0}^{\top}\epsilon^{(u)})\,$
$\displaystyle+\dfrac{1}{n}\epsilon^{(u)\top}U_{0}U_{0}^{\top}\epsilon^{(v)}+\dfrac{1}{n}y^{(u)\top}(UU^{\top}-U_{0}U_{0}^{\top})y^{(v)}.$

he first term is bounded by
$\|\lambda_{0u}\|_{2}\|\lambda_{0v}\|_{2}\|(M_{0}^{\top}M_{0}/n)-\mathbb{I}_{k}\|\lesssim
1/\sqrt{n}$ with probability at least $1-o(1)$, using Lemma E.1. The second
term is handled by first observing that
$\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}\mid M_{0}\sim
N\left(0,\sigma_{0v}^{2}\lambda_{0u}^{\top}M_{0}^{\top}M_{0}\lambda_{0u}\right)$
and then Lemma E.1 to argue that
$|\lambda_{0u}^{\top}M_{0}^{\top}\epsilon^{(v)}|\lesssim 1/\sqrt{n}$ with
probability at least $1-o(1).$ The third term is seen to be $O_{P_{0}}(1/n)$
and the fourth term is seen to be $O_{P_{0}}(n^{-1}+p_{n}^{-1})$ using
Proposition 3.5 along with $\|y^{(u)}\|_{2}\lesssim\sqrt{n}$, with probability
at least $1-o(1)$. Therefore, one has the desired result for all $1\leq
u,v\leq p$.

To derive the second result, we first observe that
$\|y^{(u)}\|_{2}^{2}\sim(\|\lambda_{0u}\|_{2}^{2}+\sigma_{0u}^{2})\,\chi_{n}^{2}$.
An application of the above result along with the weak law of large numbers
provides that this quantity converges in probability to
$\|\lambda_{0u}\|_{2}+\sigma_{0u}^{2}-\|\lambda_{0u}\|_{2}^{2}=\sigma_{0u}^{2}$
for any $1\leq u\leq p.$ ∎

######  Lemma E.5.

Let $B=(b_{uv})_{1\leq u,v\leq p}$ denote the matrix of coefficients as in
Section 2.2. Then, $\|B\|_{\infty}=O_{P_{0}}(1)$ under Assumptions 3.1-3.4.

###### Proof.

First, consider $u,v$ such that $u\neq v$. Following from Assumption 3.2,
there exists $W_{1}>0$ such that $\|\lambda_{0u}\|_{2}>c_{1}$ for all $1\leq
u\leq p$. Let $c_{2}=\min_{u}\sigma_{0u}^{2}>0$ from Assumption 3.3. Following
from Lemma E.4, we have for sufficiently large $n$ with probability at least
$1-o(1)$: (1) $\|\mu_{u}\|_{2}^{2}>c_{1}/2$ and $\mathcal{V}_{u}^{2}>c_{2}/2,$
(2)
$\left|\lambda_{0u}^{\top}\lambda_{0v}\right|\leq\|\lambda_{0u}\|_{2}\|\lambda_{0v}\|_{2}\leq
k\|\Lambda_{0}\|_{\infty}^{2}$. This implies that with probability at least
$1-o(1)$,

$b_{uv}^{2}\leq
1+\dfrac{4k^{2}\|\Lambda_{0}\|_{\infty}^{4}}{c_{1}c_{2}}<\infty.$

The result is proved analogously for $b_{uu}$ with the bound independent of
$u$. This proves the result. ∎

######  Lemma E.6.

Suppose Assumptions 3.1-3.4 hold. For each $j=1,\ldots,p$, we have

$\delta_{j}^{2}=\dfrac{Z_{j}}{n}+F_{j},$

where $Z_{j}/\sigma_{0j}^{2}\sim\chi_{n-k}^{2}$ and $\max_{1\leq j\leq
p}|F_{j}|\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}}$ with probability at least
$1-o(1)$.

###### Proof.

Suppose the SVD of $X_{0}=M_{0}\Lambda_{0}^{\top}=U_{0}D_{0}V_{0}^{\top}$. Let
$P_{U_{0}}=U_{0}(U_{0}^{\top}U_{0})^{-1}U_{0}^{\top}=U_{0}U_{0}^{\top}$ denote
the projection matrix onto the column space of $U_{0}$. We can express
$\delta_{j}^{2}$ as $\delta_{j}^{2}=\dfrac{Z_{j}}{n}+F_{j},$ where we have

$\displaystyle Z_{j}$
$\displaystyle=y^{(j)\top}\left(\mathbb{I}_{n}-P_{U_{0}}\right)y^{(j)},$
$\displaystyle F_{j}$
$\displaystyle=\dfrac{\gamma_{0}\delta_{0}^{2}}{\gamma_{n}}+\dfrac{1}{n}y^{(j)\top}\left(P_{U_{0}}-\dfrac{\widehat{M}\widehat{M}^{\top}}{n+\tau^{-2}}\right)y^{(j)}-\dfrac{\gamma_{0}}{n\gamma_{n}}y^{(j)\top}\left(\mathbb{I}_{n}-\dfrac{\widehat{M}\widehat{M}^{\top}}{n+\tau^{-2}}\right)y^{(j)}.$

Letting $\mathbf{Y}=M_{0}\Lambda_{0}^{\top}+E$ where
$E=[\epsilon^{(1)}\ldots\epsilon^{(p)}]$, we have
$(\mathbb{I}_{n}-P_{U_{0}})\mathbf{Y}=(\mathbb{I}_{n}-P_{U_{0}})E$, implying
$(\mathbb{I}_{n}-P_{U_{0}})y^{(j)}=(\mathbb{I}_{n}-P_{U_{0}})\epsilon^{(j)}\sim
N_{n}(0,\sigma_{0j}^{2}(\mathbb{I}_{n}-P_{U_{0}}))$. Since
$\mathbb{I}_{n}-P_{U_{0}}$ is idempotent (i.e.,
$(\mathbb{I}_{n}-P_{U_{0}})^{2}=\mathbb{I}_{n}-P_{U_{0}}$), we have

$\dfrac{Z_{j}}{\sigma_{0j}^{2}}=\dfrac{\|\left(\mathbb{I}_{n}-P_{U_{0}}\right)\,y^{(j)}\|_{2}^{2}}{\sigma_{0j}^{2}}\sim\chi_{\mbox{tr}(\mathbb{I}_{n}-P_{U_{0}})}^{2}\equiv\chi_{n-k}^{2}.$

We now obtain the stated probabilistic upper bound on $|F_{j}|.$ We observe

$\max_{1\leq j\leq
p}|F_{j}|\leq\dfrac{\gamma_{0}\delta_{0}^{2}}{\gamma_{n}}+\dfrac{1}{n}\left\|P_{U_{0}}-\dfrac{nUU^{\top}}{n+\tau^{-2}}\right\|\max_{1\leq
j\leq
p}||y^{(j)}||_{2}^{2}+\dfrac{\gamma_{0}}{n\gamma_{n}}\left\|\mathbb{I}_{n}-\dfrac{nUU^{\top}}{n+\tau^{-2}}\right\|\max_{1\leq
j\leq p}||y^{(j)}||_{2}^{2},$

as $\widehat{M}\widehat{M}^{\top}=nUU^{\top}$. Next, we have $\max_{1\leq
j\leq p}||y^{(j)}||_{2}^{2}\lesssim n$ with probability at least $1-o(1)$,
following Lemma 1 in Laurent and Massart, (2000). We start with

$\displaystyle G_{1}$
$\displaystyle:=\left\|P_{U_{0}}-\dfrac{nUU^{\top}}{n+\tau^{-2}}\right\|$
$\displaystyle\leq\left\|U_{0}U_{0}^{\top}-UU^{\top}\right\|+\left\|UU^{\top}-\dfrac{nUU^{\top}}{n+\tau^{-2}}\right\|$
$\displaystyle\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}}+\dfrac{1}{n}\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}},$

with probability at least $1-o(1)$. The second inequality is obtained using
Proposition 3.5. Next, we work with

$\displaystyle G_{2}:=$
$\displaystyle\,\left\|\mathbb{I}_{n}-\dfrac{nUU^{\top}}{n+\tau^{-2}}\right\|$
$\displaystyle\leq$
$\displaystyle\,\|\mathbb{I}_{n}-UU^{\top}\|+\dfrac{\tau^{-2}}{n}\|UU^{\top}\|$
$\displaystyle\lesssim$ $\displaystyle\,1.$

Combining the bounds for $G_{1}$ and $G_{2}$, we have with probability at
least $1-o(1),$ $\max_{1\leq j\leq
p}|F_{j}|\lesssim\dfrac{1}{n}+\dfrac{1}{p_{n}}.$ ∎

######  Lemma E.7.

Let $V_{n}\sim G(\gamma_{n}/2,1)$ such that $\gamma_{n}\asymp n$ and let
$U_{n}=V_{n}-(\gamma_{n}/2)$. For any $a_{n}$ satisfying $a_{n}\to 0$ and
$\sqrt{n}a_{n}\to\infty$,

$P(|U_{n}|\geq\gamma_{n}a_{n})\lesssim\exp(-na_{n}^{2}).$

###### Proof.

Immediate from Theorem 5 in Zhang and Zhou, (2020). ∎

######  Lemma E.8.

Let $M_{0}\in\mathbb{R}^{n\times k}$ be a matrix of iid $N(0,1)$ entries with
$n>k$ and let $s_{min}(M_{0})$ and $s_{max}(M_{0})$ be the smallest and
largest singular values of $M_{0}$, respectively. Then,

$\displaystyle\sqrt{n}-\sqrt{k-1}\lesssim s_{min}(M_{0})\leq
s_{max}(M_{0})\equiv\|M_{0}\|$ $\displaystyle\lesssim\sqrt{n}+\sqrt{k}$

with probability at least $1-o(1).$ In particular, if $k=o(n)$, we have
$\|M_{0}\|\asymp\sqrt{n}$ with probability at least $1-o(1).$

###### Proof.

Refer to Sections 1.1 and 1.3 in Vershynin, (2011). ∎

##  Appendix F Proofs of Related Lemmas for Theorems 3.7 and 3.8

######  Lemma F.1.

Suppose $X_{n}=Y_{n}+Z_{n},$ where
$Y_{n}\overset{}{\implies}N(0,\sigma_{1}^{2})$ and $Z_{n}\sim
N(0,\sigma_{2}^{2}),$ with $Y_{n}$ independent of $Z_{n}.$ Then,
$X_{n}\overset{}{\implies}N(0,\sigma_{1}^{2}+\sigma_{2}^{2}).$

###### Proof.

It is immediate from considering the characteristic function of $X_{n}$ and
taking limits. ∎

######  Lemma F.2.

Suppose in the setup of Theorem 3.7, we have $T_{n}=Z_{n}+Y_{n}$ for random
variables $T_{n},Z_{n},Y_{n}$ such that the pseudo-posterior of $Y_{n}$
concentrates around $0$ as $n,p\to\infty$ and the pseudo-posterior of $Z_{n}$
is the $N(0,1)$ density. Then as $n,p\to\infty,$

$\sup_{x}\,\left|{\Pi}_{C}(T_{n}\leq x)-\Phi(x)\right|\overset{P_{0}}{\to}0.$

###### Proof.

Fix $\epsilon>0$. The concentration property of $Y_{n}$ implies that
${\Pi}_{C}(|Y_{n}|>\epsilon)\overset{P_{0}}{\to}0$ as $n,p\to\infty.$ Fix
$x\in\mathbb{R}$. Using the triangle inequality, one obtains

$\displaystyle\left|{\Pi}_{C}(T_{n}\leq x)-\Phi(x)\right|$
$\displaystyle\leq\left|{\Pi}_{C}(T_{n}\leq
x,|Y_{n}|\leq\epsilon)-\Phi(x)\right|+{\Pi}_{C}(T_{n}\leq
x,|Y_{n}|>\epsilon),$

with the second term bounded by
${\Pi}_{C}(|Y_{n}|>\epsilon)\overset{P_{0}}{\to}0.$ The first term may be
decomposed as

$\displaystyle\left|{\Pi}_{C}(T_{n}\leq x,|Y_{n}|\leq\epsilon)-\Phi(x)\right|$
$\displaystyle=S_{1}+S_{2},$

where $S_{1}=|{\Pi}_{C}(T_{n}\leq
x,|Y_{n}|\leq\epsilon)-\Phi(x)|\mathbbm{1}(Y_{n}>0)$ and
$S_{2}=|{\Pi}_{C}(T_{n}\leq
x,|Y_{n}|\leq\epsilon)-\Phi(x)|\mathbbm{1}(Y_{n}\leq 0).$ We first consider

$\displaystyle S_{1}=$ $\displaystyle\left|{\Pi}_{C}(T_{n}\leq
x,|Y_{n}|\leq\epsilon)-\Phi(x)\right|\mathbbm{1}(Y_{n}>0)$ $\displaystyle\leq$
$\displaystyle\left|{\Pi}_{C}(Z_{n}\leq
x-Y_{n},Y_{n}\leq\epsilon)-{\Pi}_{C}(Z_{n}\leq
x,Y_{n}\leq\epsilon)\right|\mathbbm{1}(Y_{n}>0)+$
$\displaystyle{\Pi}_{C}(Z_{n}\leq x,Y_{n}>\epsilon)\mathbbm{1}(Y_{n}>0)$
$\displaystyle\leq$ $\displaystyle\,{\Pi}_{C}(x-Y_{n}\leq Z_{n}\leq
x,Y_{n}\leq\epsilon)\mathbbm{1}(Y_{n}>0)+{\Pi}_{C}(|Y_{n}|>\epsilon),$

where the inequality in the second line above is obtained by observing
$\Phi(x)=\Pi_{C}(Z_{n}\leq x)=\Pi_{C}(Z_{n}\leq
x,Y_{n}\leq\epsilon)+\Pi_{C}(Z_{n}\leq x,Y_{n}>\epsilon)$ and then using the
triangle inequality. The first term can be bounded by
$\Phi(x)-\Phi(x-\epsilon)\leq L\epsilon$ for a universal constant $L>0$, as
$\Phi$ is Lipschitz continuous. The second term goes to $0$ in
$P_{0}$-probability, and thus is less than $L\epsilon$ with probability at
least $1-o(1)$ for sufficiently large $n$. Thus, for sufficiently large $n$,
$S_{1}\leq 2L\epsilon$ with probability at least $1-o(1)$. Similarly,
$S_{2}\leq 2L\epsilon$ with probability at least $1-o(1)$ for sufficiently
large $n$. Combining, we have

$\sup_{x}\,\left|{\Pi}_{C}(T_{n}\leq x)-\Phi(x)\right|\leq 4L\epsilon$

with probability at least $1-o(1),$ for any fixed $\epsilon>0$. This shows the
result. We remark here that the result may be easily extended to the case
where $Z_{n}\implies N(0,1)$ as $n\to\infty$ analogously. ∎