suppReferences

# Cross-Validated Loss-Based Covariance Matrix Estimator Selection in High
Dimensions

Philippe Boileau
Graduate Group in Biostatistics,
University of California, Berkeley
<EMAIL_ADDRESS>
&Nima S. Hejazi
Graduate Group in Biostatistics, and
Center for Computational Biology,
University of California, Berkeley
<EMAIL_ADDRESS>
&Mark J. van der Laan
Division of Biostatistics,
Department of Statistics, and
Center for Computational Biology,
University of California, Berkeley
<EMAIL_ADDRESS>
&Sandrine Dudoit
Department of Statistics,
Division of Biostatistics, and
Center for Computational Biology,
University of California, Berkeley
<EMAIL_ADDRESS>

###### Abstract

The covariance matrix plays a fundamental role in many modern exploratory and
inferential statistical procedures, including dimensionality reduction,
hypothesis testing, and regression. In low-dimensional regimes, where the
number of observations far exceeds the number of variables, the optimality of
the sample covariance matrix as an estimator of this parameter is well-
established. High-dimensional regimes do not admit such a convenience,
however. As such, a variety of estimators have been derived to overcome the
shortcomings of the sample covariance matrix in these settings. Yet, the
question of selecting an optimal estimator from among the plethora available
remains largely unaddressed. Using the framework of cross-validated loss-based
estimation, we develop the theoretical underpinnings of just such an estimator
selection procedure. In particular, we propose a general class of loss
functions for covariance matrix estimation and establish finite-sample risk
bounds and conditions for the asymptotic optimality of the cross-validated
estimator selector with respect to these loss functions. We evaluate our
proposed approach via a comprehensive set of simulation experiments and
demonstrate its practical benefits by application in the exploratory analysis
of two single-cell transcriptome sequencing datasets. A free and open-source
software implementation of the proposed methodology, the cvCovEst R package,
is briefly introduced.

##  1 Introduction

The covariance matrix underlies numerous exploratory and inferential
statistical procedures central to analyses regularly performed in fields as
diverse as computational biology and financial economics. For instance, in
statistical genomics, the covariance matrix serves as a key ingredient in many
popular dimensionality reduction, clustering, and classification methods which
are regularly relied upon in quality control assessments, exploratory data
analysis, and, recently, the discovery and characterization of novel types of
cells. In financial economics, covariance matrices lie at the heart of
techniques used to gauge the excess returns of financial assets over time,
including in such applications as portfolio optimization and returns
forecasting. Still other important applications in which the covariance matrix
is critical include climate modeling and weather forecasting, imaging data
processing and analysis, probabilistic graphical modeling, and text corpora
compression and retrieval. Even more fundamentally, this statistical parameter
plays a key role in assessing the strengths of linear relationships within
multivariate data, in generating simultaneous confidence bands and regions,
and in the construction and evaluation of hypothesis tests. Accurate
estimation of the covariance matrix is therefore an essential step in the
successful application of many of the most prominent tools of modern
exploratory and inferential statistics.

When the number of observations in a dataset far exceeds the number of
features, the estimator of choice for the covariance matrix is the sample
covariance matrix: it is an efficient estimator under minimal regularity
assumptions on the data-generating distribution. In high-dimensional regimes,
however, this simple estimator has undesirable properties. When the number of
features outstrips the number of observations, the sample covariance matrix is
singular. Even when the number of observations slightly exceeds the number of
features, the sample covariance matrix is numerically unstable on account of
an overly large condition number (Golub and Van Loan, 1996). Its eigenvalues
are also generally over-dispersed when compared to those of the population
covariance matrix (Johnstone, 2001; Ledoit and Wolf, 2004): the leading
eigenvalues are positively biased, while the trailing eigenvalues are
negatively biased.

High-dimensional data have become increasingly widespread in myriad scientific
domains, with many notable examples arising from challenges posed by cutting-
edge biological sequencing technologies. Accordingly, researchers have turned
to developing novel covariance matrix estimators to remediate the sample
covariance matrix’s deficiencies. Under certain sparsity assumptions, Bickel
and Levina (2008a, b); Rothman et al. (2009); Lam and Fan (2009); Cai et al.
(2010), and Cai and Liu (2011), among others, have demonstrated that the true
covariance matrix can be estimated consistently under certain losses through
the application of element-wise thresholding or tapering functions to the
sample covariance matrix. Another thread of the literature, which includes
notable contributions by Stock and Watson (2002); Bai (2003); Fan et al.
(2008, 2013, 2016, 2019), and Onatski (2012), has championed methods employing
factor models in covariance matrix estimation. Other popular proposals include
the families of estimators formulated by Schäfer and Strimmer (2005), and
Ledoit and Wolf (2004, 2012, 2015, 2018), which use an approach inspired by
the empirical Bayes framework (Robbins, 1964; Efron, 2012). For convenience,
we briefly review and detail several of these estimator families in Section
S1.

Despite the flexibility provided by an apparent wealth of candidate
estimators, this variety poses many practical issues. Namely, identifying the
most appropriate estimator from among a collection of candidates is itself a
significant challenge. A partial answer to this problem has come in the form
of data-adaptive approaches designed to select the optimal estimator within a
particular class (e.g., Bickel and Levina, 2008a, b; Cai and Liu, 2011; Fang
et al., 2016; Fan et al., 2013). Such approaches, however, are inherently
limited by their focus on relatively narrow families of covariance matrix
estimators. The successful application of such estimator selection frameworks
requires, as a preliminary step, that the practitioner make a successful
choice among estimator families, thus injecting a degree of subjectivity in
their deployment. On the other hand, the broader question of selecting an
optimal estimator from among a diverse library of candidates has remained
unaddressed. In answer to this, we offer a general loss-based framework
building upon the seminal work of Dudoit and van der Laan (2003); van der Laan
and Dudoit (2003); van der Vaart et al. (2006) for asymptotically optimal
covariance matrix estimator selection based upon cross-validation.

In the cross-validated loss-based estimation framework, the parameter of
interest is defined as the risk minimizer, with respect to the data-generating
distribution, based on a loss function chosen to reflect the problem at hand.
Candidate estimators may be generated in a variety of manners, including as
empirical risk minimizers with respect to an empirical distribution over
parameter subspaces corresponding to models for the data-generating
distribution. One would ideally select as optimal estimator that which
minimizes the “true” risk with respect to the data-generating distribution.
However, as this distribution is typically unknown, one turns to cross-
validation for risk estimation. Dudoit and van der Laan (2003); van der Laan
and Dudoit (2003); van der Vaart et al. (2006) have shown that, under general
conditions on the data-generating distribution and loss function, the cross-
validated estimator selector is asymptotically optimal, in the sense that it
performs asymptotically as well in terms of risk as an optimal oracle selector
based on the true, unknown data-generating distribution.

Focusing specifically on the covariance matrix as the parameter of interest,
we address the choice of loss function and candidate estimators, and derive
new, high-dimensional asymptotic optimality results for multivariate cross-
validated estimator selection procedures. Requiring generally nonrestrictive
assumptions about the structure of the true covariance matrix, the proposed
framework accommodates arbitrary families of covariance matrix estimators. The
method therefore enables the objective selection of an optimal estimator while
fully taking advantage of the plethora of available estimators. As such, it
generalizes existing, but more limited, data-adaptive estimator selection
frameworks where the library of candidate estimators is narrowed based on
available subject matter knowledge, or, as is more commonly the case, for
convenience’s sake.

The remainder of the present manuscript is organized as follows. In Section 2,
we formulate the problem within the general framework of cross-validated loss-
based estimation. In Section 3, we propose a general class of loss functions
for covariance matrix estimation and establish finite-sample risk bounds and
conditions for the asymptotic optimality of the cross-validated estimator
selector with respect to these loss functions. Section 4 provides a thorough
interrogation of the proposed cross-validated estimator selector in a series
of simulation experiments. The practical benefits of using the cross-validated
estimator selection strategy are then demonstrated by its application in a
dimensionality reduction pipeline for novel cell-type identification with two
distinct single-cell RNA-seq datasets in Section 5. We conclude the manuscript
with a brief discussion in Section 6.

##  2 Problem Formulation and Background

Consider a dataset $\mathbf{X}_{n\times
J}=\\{X_{1},\ldots,X_{n}:X_{i}\in\mathbb{R}^{J}\\}$, comprising $n$
independent and identically distributed (i.i.d.) random vectors, where
$n\approx J$ or $n<J$. Let $X_{i}\sim P_{0}\in\mathcal{M}$, where $P_{0}$
denotes the true data-generating distribution and $\mathcal{M}$ the
statistical model, i.e., a collection of possible data-generating
distributions $P$ for $X_{i}$. We assume a nonparametric statistical model
$\mathcal{M}$ for $P_{0}$, minimizing assumptions on the form of $P_{0}$. We
denote by $P_{n}$ the empirical distribution of the $n$ random vectors forming
$\mathbf{X}_{n\times J}$. Letting $\mathbb{E}[X_{i}]=0$ without loss of
generality and defining $\boldsymbol{\psi}_{0}\coloneqq\text{Var}[X_{i}]$, we
take as our goal the estimation of the covariance matrix
$\boldsymbol{\psi}_{0}$. This is accomplished by identifying the “optimal”
estimator of $\boldsymbol{\psi}_{0}$ from among a collection of candidates,
where, as detailed below, optimality is defined in terms of risk.

For any distribution $P\in\mathcal{M}$, define its covariance matrix as
$\boldsymbol{\psi}=\Psi(P)$, in terms of a mapping $\Psi$ from the model
$\mathcal{M}$ to the space of $J\times J$ symmetric, positive semi-definite
matrices. Given this mapping, define the parameter space as
$\boldsymbol{\Psi}\coloneqq\\{\boldsymbol{\psi}=\Psi(P):P\in\mathcal{M}\\}$,
the set of all covariance matrices corresponding to the data-generating
distributions $P$ belonging to the model $\mathcal{M}$. Furthermore, candidate
estimators of the covariance matrix are defined as
$\hat{\boldsymbol{\psi}}_{k}\coloneqq\hat{\Psi}_{k}(P_{n})$, in terms of
mappings $\hat{\Psi}_{k}$ from the empirical distribution $P_{n}$ to the
parameter space $\boldsymbol{\Psi}$, $k=1,\ldots,K$. While this notation
suggests that the number of candidate estimators $K$ is fixed (and we treat it
as such throughout this work), we note that this framework may be extended
such that $K$ grows as a function of $n$ and $J$.

In order to assess the optimality of estimators in the set
$\boldsymbol{\Psi}$, we introduce a generic loss function
$L(X;\boldsymbol{\psi},\eta)$ characterizing a cost applicable to any
$\boldsymbol{\psi}\in\boldsymbol{\Psi}$ and $X\sim P\in\mathcal{M}$, and where
$\eta$ is a (possibly empty) nuisance parameter. Specific examples of loss
functions for the covariance estimation setting are proposed in Section 3.1.
Define $H$ as the mapping from the model $\mathcal{M}$ to the nuisance
parameter space $\mathbf{H}\coloneqq\\{\eta=H(P):P\in\mathcal{M}\\}$ and let
$\hat{\eta}\coloneqq\hat{H}(P_{n})$ denote a generic nuisance parameter
estimator, where $\hat{H}$ is a mapping from $P_{n}$ to $\mathbf{H}$. Given
any $\eta\in\mathbf{H}$, the risk under $P\in\mathcal{M}$ for any
$\boldsymbol{\psi}\in\boldsymbol{\Psi}$ is defined as the expected value of
$L(X;\boldsymbol{\psi},\eta)$ with respect to $P$:

$\begin{split}\Theta(\boldsymbol{\psi},\eta,P)&\coloneqq\int
L(x;\boldsymbol{\psi},\eta)dP(x)\\\
&=\mathbb{E}_{P}\left[L(X;\boldsymbol{\psi},\eta)\right].\end{split}$ (1)

We then place the additional constraint on the loss function that as long as a
risk minimizer exists under the true data-generating distribution $P_{0}$, the
minimizer is given by the parameter of interest

$\boldsymbol{\psi}_{0}\coloneqq\operatorname*{\arg\\!\min}_{\boldsymbol{\psi}\in\mathbf{\Psi}}\Theta(\boldsymbol{\psi},\eta_{0},P_{0}),$
(2)

where $\eta_{0}\coloneqq H(P_{0})$. The risk minimizer need not be unique. The
optimal risk under $P_{0}$ is then given by

$\theta_{0}\coloneqq\min_{\boldsymbol{\psi}\in\mathbf{\Psi}}\;\Theta(\boldsymbol{\psi},\eta_{0},P_{0}),$

which is to say that a risk minimizer $\boldsymbol{\psi}_{0}$ attains risk
$\theta_{0}$.

For any given estimator $\hat{\boldsymbol{\psi}}_{k}$ of
$\boldsymbol{\psi}_{0}$, its conditional risk given $P_{n}$ with respect to
the true data-generating distribution $P_{0}$ is defined as

$\begin{split}\tilde{\theta}_{n}(k,\eta_{0})&\coloneqq\mathbb{E}_{P_{0}}[L(X;\hat{\Psi}_{k}(P_{n}),\eta_{0})\mid
P_{n}]\\\ &=\int L(x;\hat{\Psi}_{k}(P_{n}),\eta_{0})dP_{0}(x)\\\
&=\Theta(\hat{\boldsymbol{\psi}}_{k},\eta_{0},P_{0}).\end{split}$ (3)

Defining the risk difference of the $k^{\text{th}}$ estimator as
$\tilde{\theta}_{n}(k,\eta_{0})-\theta_{0}$, the index of the estimator that
achieves the minimal risk difference is denoted by

$\tilde{k}_{n}\coloneqq\operatorname*{\arg\\!\min}_{k\in\\{1,\ldots,K\\}}\;\tilde{\theta}_{n}(k,\eta_{0})-\theta_{0}.$

The subscript $n$ emphasizes that the risk and optimal estimator index are
conditional on the empirical distribution $P_{n}$. They are therefore random
variables.

Given the high-dimensional nature of the data, it is generally most convenient
to study the performance of estimators of $\boldsymbol{\psi}_{0}$ using
Kolmogorov asymptotics — that is, in the setting in which both
$n\rightarrow\infty$ and $J\rightarrow\infty$ such that $J/n\rightarrow
m<\infty$. Historically, estimators have been derived within this high-
dimensional asymptotic regime to improve upon the finite sample results of
estimators brought about by traditional asymptotic arguments. After all, the
sample covariance matrix retains its asymptotic optimality properties when $J$
is fixed, even though it is known to perform poorly in high-dimensional
settings.

Naturally, it would be desirable for an estimator selection procedure to
select the estimator indexed by $\tilde{k}_{n}$; however, this quantity
depends on the true, unknown data-generating distribution $P_{0}$. As a
substitute for the candidates’ true conditional risks, we employ instead the
cross-validated estimators of these same conditional risks.

Cross-validation (CV) consists of randomly, and possibly repeatedly,
partitioning a dataset into a training set and a validation set. The
observations in the training set are fed to the candidate estimators and the
observations in the validation set are used to evaluate the performance of
these estimators (Breiman and Spector, 1992; Friedman et al., 2001). A range
of CV schemes have been proposed and investigated, both theoretically and
computationally; Dudoit and van der Laan (2005) provide a thorough review of
popular CV schemes and their properties. Among the variety, $V$-fold and Monte
Carlo CV stand out as approaches that have gained traction on account of their
relative computational feasibility and good performance. Any CV scheme can be
expressed in terms of a binary random vector $B_{n}$, which assigns
observations into either the training or validation set. Observation $i$ is
said to lie in the training set when $B_{n}(i)=0$ and in the validation set
otherwise. The training set therefore contains
$\sum_{i}(1-B_{n}(i))=n(1-p_{n})$ observations and the validation set
$\sum_{i}B_{n}(i)=np_{n}$ observations, where $p_{n}$ is the fixed validation
set proportion corresponding to the chosen CV procedure. The empirical
distributions of the training and validation sets are denoted by
$P_{n,B_{n}}^{0}$ and $P_{n,B_{n}}^{1}$, respectively, for any given
realization of $B_{n}$.

Using this general definition of CV, we may now define the cross-validated
estimator of a candidate $\hat{\Psi}_{k}$’s risk as

$\begin{split}\hat{\theta}_{p_{n},n}(k,\hat{H}(P_{n,B_{n}}^{0}))&\coloneqq\mathbb{E}_{B_{n}}\left[\Theta(\hat{\Psi}_{k}(P_{n,B_{n}}^{0}),\hat{H}(P_{n,B_{n}}^{0}),P_{n,B_{n}}^{1})\right]\\\
&=\mathbb{E}_{B_{n}}\left[\frac{1}{np_{n}}\sum_{i=1}^{n}\mathbb{I}(B_{n}(i)=1)L(X_{i};\hat{\Psi}_{k}(P^{0}_{n,B_{n}}),\hat{H}(P_{n,B_{n}}^{0}))\right],\end{split}$
(4)

for some nuisance parameter estimator mapping $\hat{H}$. The corresponding
cross-validated selector is then given by

$\hat{k}_{p_{n},n}\coloneqq\operatorname*{\arg\\!\min}_{k\in\\{1,\ldots,K\\}}\hat{\theta}_{p_{n},n}(k,\hat{H}(P_{n,B_{n}}^{0})).$

As a benchmark, the unknown cross-validated conditional risk of the $k$th
estimator is

$\tilde{\theta}_{p_{n},n}(k,\eta_{0})\coloneqq\mathbb{E}_{B_{n}}\left[\Theta(\hat{\Psi}_{k}(P_{n,B_{n}}^{0}),\eta_{0},P_{0})\right].$

The cross-validated oracle selector is then given by

$\tilde{k}_{p_{n},n}\coloneqq\operatorname*{\arg\\!\min}_{k\in\\{1,\ldots,K\\}}\tilde{\theta}_{p_{n},n}(k,\eta_{0}).$

As before, the $p_{n}$ and $n$ subscripts highlight the dependence of these
objects on the CV procedure and the empirical distribution $P_{n}$,
respectively (thus making them random variables).

Ideally, the cross-validated estimator selection procedure should identify a
$\hat{k}_{p_{n},n}$ that is asymptotically (in $n$, $J$, and possibly $K$)
equivalent in terms of risk to the oracle $\tilde{k}_{p_{n},n}$, under a set
of nonrestrictive assumptions based on the choice of loss function, target
parameter space, and, if applicable, nuisance parameter space, in the sense
that

$\frac{\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n},\eta_{0})-\theta_{0}}{\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n},\eta_{0})-\theta_{0}}\overset{P}{\longrightarrow}1\text{
as }n,J\rightarrow\infty.$ (5)

That is, the estimator selected via CV is equivalent in terms of risk to the
CV scheme’s oracle estimator chosen from among all candidates.

When Equation (5) holds, a further step may be taken by relating the
performance of the cross-validated selector to that of the full-dataset oracle
selector, $\tilde{k}_{n}$:

$\frac{\tilde{\theta}_{n}(\hat{k}_{p_{n},n},\eta_{0})-\theta_{0}}{\tilde{\theta}_{n}(\tilde{k}_{n},\eta_{0})-\theta_{0}}\overset{P}{\longrightarrow}1\text{
as }n,J\rightarrow\infty.$ (6)

When the cross-validated selection procedure’s full-dataset conditional risk
difference converges in probability to that of the full-dataset oracle’s, we
say that its chosen estimator is asymptotically optimal. In other words, the
data-adaptively selected estimator performs asymptotically as well, with
respect to the chosen loss, as the candidate that would be picked from the
collection of estimators if the true data-generating distribution were known.

Next, we focus on identifying a loss function with which to perform this
procedure in the covariance estimation setting, and on determining the
conditions under which such a procedure generates optimal results.

##  3 Loss Functions and Estimator Selection

###  3.1 Proposed Loss Function

The choice of loss function should reflect the goals of the estimation task.
While loss functions based on the sample covariance matrix and either the
Frobenius or the spectral norms are often employed in the covariance matrix
estimation literature, Dudoit and van der Laan (2003)’s estimator selection
framework is more amenable to loss functions that operate over random vectors.
Accordingly, we propose the observation-level Frobenius loss:

$\begin{split}L(X;\boldsymbol{\psi},\eta_{0})&\coloneqq\lVert
XX^{\top}-\boldsymbol{\psi}\rVert^{2}_{F,\eta_{0}}\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}\eta_{0}^{(jl)}(X^{(j)}X^{(l)}-\psi^{(jl)})^{2},\end{split}$
(7)

where $X^{(j)}$ is the $j$th element of a random vector $X\sim
P\in\mathcal{M}$, $\psi^{(jl)}$ is the entry in the $j$th row and $l$th column
of an arbitrary covariance matrix $\boldsymbol{\psi}\in\mathbf{\Psi}$, and
$\eta_{0}$ is a $J\times J$ matrix acting as a scaling factor (i.e., a
potential nuisance parameter). For an estimator $\hat{\eta}$ of $\eta_{0}$,
the cross-validated risk estimator of the $k$th candidate estimator
$\hat{\Psi}_{k}$ under the observation-level Frobenius loss is

$\begin{split}\hat{\theta}_{p_{n},n}(k,\hat{H}(P_{n,B_{n}}^{0}))&=\mathbb{E}_{B_{n}}\left[\Theta(\hat{\Psi}_{k}(P_{n,B_{n}}^{0}),\hat{H}(P_{n,B_{n}}^{0}),P_{n,B_{n}}^{1})\right]\\\
&=\mathbb{E}_{B_{n}}\left[\frac{1}{np_{n}}\sum_{i=1}^{n}\mathbb{I}(B_{n}(i)=1)\lVert
X_{i}X^{\top}_{i}-\hat{\Psi}_{k}(P_{n,B_{n}}^{0})\rVert^{2}_{F,\hat{H}(P_{n,B_{n}}^{0})}\right].\end{split}$
(8)

Ledoit and Wolf (2004), Bickel and Levina (2008b), and Rothman et al. (2009),
among others, have employed analogous (scaled) Frobenius losses to prove
various optimality results, defining $\eta_{0}^{(jl)}=1/J$, $\forall j,l$.
This particular choice of scaling factor is such that whatever the value of
$J$, $\lVert\mathbf{I}_{J\times J}\rVert_{F,\eta_{0}}=1$. With such a scaling
factor, the loss function may be viewed as a relative loss whose yardstick is
the $J\times J$ identity matrix. A similarly reasonable option for when the
true covariance matrix is assumed to be dense is $\eta_{0}^{(jl)}=1/J^{2}$.
This weighting scheme effectively computes the average squared error across
every entry of the covariance matrix. Note, however, that when the scaling
factor is constant, it only impacts the interpretation of the loss. Constant
scaling factors have no impact on our asymptotic analysis. Thus, it is not a
nuisance parameter in the conventional sense since it need not be estimated.

Further, when the scaling factor is constant, the risk minimizers are
identical for the cross-validated observation-level Frobenius risk and the
common cross-validated Frobenius risk.

######  Proposition 1.

Define the cross-validated Frobenius risk for an estimator $\hat{\Psi}_{k}$ as

$\hat{R}_{n}(\hat{\Psi}_{k},\eta_{0})\coloneqq\mathbb{E}_{B_{n}}\left[\lVert\mathbf{S}_{n}(P_{n,B_{n}}^{1})-\hat{\Psi}_{k}(P_{n,B_{n}}^{0})\rVert^{2}_{F,\eta_{0}}\right],$
(9)

where $\mathbf{S}_{n}(P_{n,B_{n}}^{1})$ is the sample covariance matrix
computed over the validation set $P_{n,B_{n}}^{1}$, and $\eta_{0}$ is some
constant scaling matrix. Then,

$\begin{split}\hat{k}_{p_{n},n}&=\operatorname*{\arg\\!\min}_{k\in\\{1,\ldots,K\\}}\hat{\theta}_{p_{n},n}(k,\eta_{0})\\\
&=\operatorname*{\arg\\!\min}_{k\in\\{1,\ldots,K\\}}\hat{R}_{n}(\hat{\Psi}_{k},\eta_{0}).\end{split}$

Note that the traditional Frobenius loss corresponds to the sum of the squared
eigenvalues of the difference between the sample covariance matrix and the
estimate. Proposition 1 therefore implies the existence of a similar
relationship for our observation-level Frobenius loss. It may therefore serve
as a surrogate for a loss based on the spectral norm.

We are not restricted to a constant scaling factor matrix. One might consider
weighting the covariance matrix’s off-diagnoal elements’ errors by their
corresponding diagonal entries, especially useful when the random variables
are of different magnitudes. Such a scaling factor might offer a more
equitable evaluation across all entries of the parameter:

$\begin{split}L_{\text{weighted}}(X;\boldsymbol{\psi},\eta_{0})&\coloneqq\sum_{j=1}^{J}\sum_{l=1}^{J}\eta_{0}^{(jl)}(X^{(j)}X^{(l)}-\psi^{(jl)})^{2}\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}\frac{1}{\sqrt{\psi_{0}^{(jj)}\psi_{0}^{(ll)}}}(X^{(j)}X^{(l)}-\psi^{(jl)})^{2}.\end{split}$
(10)

Here, $\eta_{0}$ is a genuine nuisance parameter which can be estimated via
the diagonal entries of the sample covariance matrix.

Finally, the covariance matrix $\psi_{0}$ is the risk minimizer of the
observation-level Frobenius loss if the integral with respect to $X$ and the
partial differential operators with respect to $\boldsymbol{\psi}$ are
interchangeable. Our proposed loss therefore satisfies the condition of
Equation (2). The main results of the paper, however, relate only to the
constant scaling factor case. In a minor abuse of notation, we set
$\eta_{0}=\boldsymbol{1}$, and suppress dependence of the loss function on the
scaling factor wherever possible throughout the remainder of the text.

###  3.2 Optimality of the Cross-validated Estimator Selector

Having defined a suitable loss function, we turn to a discussion of the
theoretical properties of the cross-validated estimator selection procedure.
Specifically, we present, in Theorem 1, sufficient conditions under which the
method is asymptotically equivalent in terms of risk to the commensurate CV
oracle selector (as per Equation (5)). This theorem extends the general
framework of Dudoit and van der Laan (2003) for use in high-dimensional
multivariate estimator selection. Adapting their existing theory to this
setting requires a judicious choice of loss function, new assumptions, and
updated proofs reflecting the use of high-dimensional asymptotics. Corollary 1
then builds on Theorem 1 and details conditions under which the procedure
produces asymptotically optimal selections in the sense of Equation (6).

######  Theorem 1.

Let $X_{1},\ldots,X_{n}$ be a random sample of $n$ i.i.d. random vectors of
dimension $J$, such that $X_{i}\sim P_{0}\in\mathcal{M},i=1,\ldots,n$. Assume,
without loss of generality, that $\mathbb{E}[X_{i}]=0$, and define
$\boldsymbol{\psi}_{0}\coloneqq\text{Var}[X_{i}]$. Additionally, let
$\mathbf{\Psi}\coloneqq\\{\psi=\Psi(P):P\in\mathcal{M}\\}$ be the parameter
space of all possible covariance matrices for the model $\mathcal{M}$. Denote
the set of $K$ candidate estimator mappings as
$\\{\hat{\Psi}_{k}(\cdot):k=1,\ldots,K\\}$, and assume that the empirical
distribution $P_{n}\in\mathcal{M}$. Next, for
$\boldsymbol{\psi}\in\mathbf{\Psi}$, define the observation-level Frobenius
loss function as $L(X;\boldsymbol{\psi})\coloneqq\lVert
X^{\top}X-\boldsymbol{\psi}\rVert^{2}_{F,1}$. Finally, designate $p_{n}$ as
the proportion of observations in any given cross-validated validation set.

Consider the following assumptions:

######  Assumption 1.

For any $P\in\mathcal{M}$ and $X\sim P$, $\max_{j=1,\ldots,J}\>(\lvert
X^{(j)}\rvert)<\sqrt{M_{1}}<\infty$ almost surely (a.s.).

######  Assumption 2.

For each $\boldsymbol{\psi}\in\mathbf{\Psi}$,
$\lvert{\psi}^{(jl)}\rvert<M_{2}<\infty$, $\forall\;j,l=1,\ldots,J$.

Finite-Sample Result. Let $\overline{M}(J)\coloneqq 4(M_{1}+M_{2})^{2}J^{2}$
and $c(\delta,\overline{M}(J))\coloneqq
2(1+\delta)^{2}\overline{M}(J)(1/3+1/\delta)$. Then, for any $\delta>0$,

$0\leq\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}]\leq(1+2\delta)\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]+2c(\delta,\overline{M}(J))\frac{1+\log(K)}{np_{n}}.$
(11)

High-Dimensional Asymptotic Result. The finite-sample result in Equation (11)
has the following asymptotic implications: If
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}])\rightarrow
0$ and $J/n\rightarrow m<\infty$ as $n,J\rightarrow\infty$, then

$\frac{\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}]}{\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]}\rightarrow
1.$ (12)

Further, if
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}))\overset{P}{\rightarrow}0$
as $n,J\rightarrow\infty$, then

$\frac{\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}}{\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}}\overset{P}{\rightarrow}1.$
(13)

The proof is provided in Section S2. It relies on special properties of the
random variable
$Z_{k}~{}\coloneqq~{}L(X;\hat{\Psi}_{k}(P_{n}))~{}-~{}L(X;\boldsymbol{\psi}_{0})$
and on an application of Bernstein’s inequality (Bernstein, 1946; Bennett,
1962). Together, they are used to show that
$2c(\delta,\overline{M}(J))(1+\log(K))/(np_{n})$ is a finite-sample bound for
comparing the performance of the cross-validated selector,
$\hat{k}_{p_{n},n}$, against that of the oracle selector over the training
sets, $\tilde{k}_{p_{n},n}$. Once this bound is established, the high-
dimensional asymptotic results follow immediately.

Only a few sufficient conditions are required to provide finite-sample bounds
on the expected risk difference of the estimator selected via our CV
procedure. First, that each element of the random vector $X$ be bounded, and,
second, that the entries of all covariance matrices in the parameter space be
bounded. Together, these assumptions allow for the definition of
$\overline{M}(J)$, the object permitting the extension of the loss-based
estimation framework to the high-dimensional covariance matrix estimation
problem.

The first assumption may appear stringent. It is not satisfied, for instance,
by Gaussian distributions. However, we stress that parametric data-generating
distributions rarely reflect reality — they are merely mathematically
convenient. We posit that most random variables arising in scientific
practice, or transformations thereof, are bounded by limitations of the
physical, electronic, or biological measurement process; thus, our method
remains widely applicable. For example, in microarray and next-generation
sequencing experiments, the raw data are images on a 16-bit scale and are
therefore constrained to belong to $[0,2^{16})$. Similarly, the measurement of
immune correlates — of substantial interest in vaccine efficacy trials of
HIV-1, COVID-19, and other diseases — are bounded by the limits of detection
and quantitation imposed by the assays used in their measurement.

Verifying that the additional assumptions required by Theorem 1’s asymptotic
results hold proves to be more challenging. We write $f(y)=O(g(y))$ if $\lvert
f\rvert$ is bounded above by $g$, $f(y)=o(g(y))$ if $f$ is dominated by $g$,
$f(y)=\Omega(g(y))$ if $f$ is bounded below by $g$, and $f(y)=\omega(g(y))$ if
$f$ dominates $g$, all in asymptotics with respect to $n$ and $J$. Further, a
subscript “P” might be added to these bounds, denoting a convergence in
probability. Now, note that
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n})=O(J)$ for fixed $p_{n}$ and as
$J/n\rightarrow m>0$. Then the conditions associated with Equation (12) and
Equation (13) hold so long as
$\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]=\omega(J)$
and $\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}=\omega_{P}(J)$,
respectively.

These requirements do not seem particularly restrictive given that the
complexity of the problem generally increases as a function of the number of
features. There are many more entries in the covariance matrix requiring
estimation than there are observations. Consistent estimation in terms of the
Frobenius risk is therefore not possible in high-dimensions without additional
assumptions about the data-generating process. This intuition is corroborated
by our extensive simulation study in the following section.

Some additional insight might be gained by identifying conditions under which
these assumptions are unmet for popular structural beliefs about the true
covariance matrix. In particular, we consider the sparse covariance matrices
defined in Bickel and Levina (2008b) and accompanying hard-thresholding
estimators (see Section S1.1):

######  Proposition 2.

In addition to Assumptions 1 and 2 of Theorem 1, assume that the hard-
thresholding estimator is in the library of candidates, that there is a
correspondingly large grid of possible thresholding hyperparameters, and that
$\psi_{0}$ is a member of the following set of matrices:

$\left\\{\psi:\psi^{(jj)}<M_{2},\sum_{l=1}^{J}I(\psi^{(jl)}\neq 0)<s(J)\text{
for all },j=1,\ldots,J\right\\},$ (14)

where $s(J)$ is the row sparsity. Then, by Theorems 2 of Bickel and Levina
(2008b),
$\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]=o(J)$
if $s(J)/J=o(1/\log J)$ asymptotically in $n$ and $J$.

Proposition 2 states that the conditions for achieving the asymptotic results
of Theorem 1 are not met if the proportion of non-zero elements in the
covariance matrix’s row with the most non-zero elements converges to zero
faster than $1/\log J$ and the library of candidates possesses a hard-
thresholding estimator whose thresholding hyperparameter is reasonable in the
sense of Bickel and Levina (2008b)’s Theorem 2 and its subsequent discussion.
Plainly, the true covariance matrix cannot be too sparse if the collection of
considered estimators contains a diversity of hard-thresholding estimators.

This implies that banded covariance matrices whose number of bands are fixed
for $J$ do not meet the criteria for our theory to apply, assuming that one of
the candidate estimators correctly specifies the number of bands.
Nevertheless, we observe empirically in Section 4 that our cross-validated
procedure selects an optimal estimator when the true covariance matrix is
banded or tapered more quickly in terms of $n$ and $J$ than any other type of
true covariance matrix.

These results are likely explained by the relatively low complexity of the
estimation problem in this setting. High-dimensional asymptotic arguments are
perhaps unnecessary when the proportion of entries needing to be estimated in
the true covariance matrix quickly converges to zero. In any case, these
limitations of our theory reflect stringent – and typically unverifiable –
structural assumptions about the estimand. We reiterate that the conditions of
Theorem 1 are generally satisfied.

Now, Theorem 1’s high-dimensional asymptotic results relate the performance of
the cross-validated selector to that of the oracle selector for the CV scheme.
As indicated by the expression in Equation (6), however, we would like our
cross-validated procedure to be asymptotically equivalent to the oracle over
the entire dataset. The conditions to obtain this desired result are provided
in Corollary 1, a minor adaptation of previous work by Dudoit and van der Laan
(2003). This extension accounts for increasing $J$, thereby permitting its use
in high-dimensional asymptotics.

######  Corollary 1.

Building upon Assumptions 1 and 2 of Theorem 1, we introduce the additional
assumptions that, as $n,J\rightarrow\infty$ and $J/n\rightarrow m<\infty$,
$p_{n}\rightarrow 0$,
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}))\overset{P}{\rightarrow}0$,
and

$\frac{\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}}{\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0}}\overset{P}{\rightarrow}1.$
(15)

Under these assumptions, it follows that

$\frac{\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}}{\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0}}\overset{P}{\longrightarrow}1.$
(16)

The proof is provided in Section S2 and is a direct application of the
asymptotic results of Theorem 1.

As before, the assumption that
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}))\overset{P}{\rightarrow}0$
remains difficult to verify, but essentially requires the estimation error of
the oracle to increase quickly as the number of features grows. That is,
$np_{n}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0})=\omega_{P}(J)$.
We posit that this condition is generally satisfied, similarly to the
asymptotic results of Theorem 1.

Now, a sufficient condition for Equation (15) is that there exists a
$\gamma>0$ such that

$\left(n^{\gamma}(\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0}),\;(n(1-p_{n}))^{\gamma}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0})\right)\overset{d}{\rightarrow}(Z,Z),$
(17)

for a single random variable $Z$ with $\mathbb{P}(Z>a)=1$ for some $a>0$. For
single-split validation, where $\mathbb{P}(B_{n}=b)=1$ for some
$b\in\\{0,1\\}^{n}$, it suffices to assume that there exists a $\gamma>0$ such
that
$n^{\gamma}(\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0})\overset{d}{\rightarrow}Z$
for a random variable $Z$ with $\mathbb{P}(Z>a)=1$ for some $a>0$.

Equation (15) essentially requires that the (appropriately scaled)
distributions of the cross-validated and full-dataset conditional risk
differences of their respective oracle selections converge in distribution as
$p_{n}\rightarrow 0$. Again, this condition is unrestrictive. As
$p_{n}\rightarrow 0$, the composition of each training set becomes
increasingly similar to that of the full-dataset. The resulting estimates
produced by each candidate estimator over the training sets and the full-
dataset will correspondingly converge. Naturally, so too will the cross-
validated and full-dataset conditional risk difference distributions of their
respective selections.

Finally, while the number of candidates in the estimator library $K$ has been
assumed to be fixed in this discussion of the proposed method’s asymptotic
results, it may be allowed to grow as a function of $n$ and $J$. Of course,
this will negatively impact the convergence rates of
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}])$
and
$c(\delta,\overline{M}(J))(1+\log(K))/(np_{n}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}))$.
The sufficient conditions outlined in the asymptotic results of Theorem 1 are
achieved so long as the library of candidates does not grow too aggressively.
That is, we can make the additional assumptions that
$K=o(\text{exp}\\{\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]/J\\})$
and
$K=o_{P}(\text{exp}\\{(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0})/J\\})$
such that the results of Equations (12) and (13) are achieved, respectively.

Having provided a formal description of the procedure’s asymptotic behavior,
what remains is to provide empirical evidence that these results may be
approximately achieved under a variety of covariance matrix structures and
estimators, and for values of $n$ and $J$ encountered in practice. This is
explored via a simulation study in Section 4.

###  3.3 Software Implementation

The proposed cross-validated covariance matrix estimator selection procedure
has been implemented in cvCovEst (Boileau et al., 2021), a free and open-
source software package for the R language and environment for statistical
computing (R Core Team, 2021). Internally, the cvCovEst package relies upon
the generalized CV framework of the origami R package (Coyle and Hejazi,
2018). The observation-level Frobenius loss with a constant scaling factor (as
in Equation (7)) is used to compute the cross-validated risks, though a few
other choices of loss functions are provided as well. In particular, users may
opt for the matrix-based Frobenius loss of Equation (9), which has the benefit
of carrying a reduced computational burden, to select the asymptotically
optimal estimator.

Results for the simulation experiments and real data examples, presented in
Sections 4 and 5, respectively, were produced using version 0.1.3 of the
cvCovEst R package. The code and data used to produce these analyses are
publicly available in a GitHub repository at
https://github.com/PhilBoileau/pub_cvCovEst. The cvCovEst package is available
via the Comprehensive R Archive Network (CRAN) at
https://CRAN.R-project.org/package=cvCovEst.

##  4 Simulation Study

###  4.1 Simulation Study Design

We conducted a series of simulation experiments using prominent covariance
models to verify the theoretical results of our cross-validated estimator
selection procedure. These models are described below.

Model 1:

A dense covariance matrix, where

$\psi^{(jl)}=\begin{cases}1,&j=l\\\ 0.5,&\text{otherwise}\\\ \end{cases}.$

Model 2:

An AR(1) model, where $\psi^{(jl)}=0.7^{\lvert j-l\rvert}$. This covariance
matrix, corresponding to a common timeseries model, is approximately sparse
for large $J$, since the off-diagonal elements quickly shrink to zero.

Model 3:

An MA(1) model, where $\psi^{(jl)}=0.7^{\lvert
j-l\rvert}\cdot\mathbb{I}(\lvert j-l\rvert\leq 1)$. This covariance matrix,
corresponding to another common timeseries model, is truly sparse. Only the
diagonal, subdiagonal, and superdiagonal contain non-zero elements.

Model 4:

An MA(2) model, where

$\psi^{(jl)}=\begin{cases}1,&j=l\\\ 0.6,&\lvert j-l\rvert=1\\\ 0.3,&\lvert
j-l\rvert=2\\\ 0,&\text{otherwise}\end{cases}.$

This timeseries model is similar to Model 3, but slightly less sparse.

Model 5:

A random covariance matrix model. First, a $J\times J$ random matrix whose
elements are i.i.d. $\text{Uniform}(0,1)$ is generated. Next, entries below
$1/4$ are set to $1$, entries between $1/4$ and $1/2$ are set to $-1$, and the
remaining entries are set to $0$. The square of this matrix is then computed
and added to the identity matrix $\mathbf{I}_{J\times J}$. Finally, the
corresponding correlation matrix is computed and used as the model’s
covariance matrix.

Model 6:

A Toeplitz covariance matrix, where

$\psi^{(jl)}=\begin{cases}1,&j=l\\\ 0.6\lvert
j-l\rvert^{-1.3},&\text{otherwise}\end{cases}.$

Like the AR(1) model, this covariance matrix is approximately sparse for large
$J$. However, the off-diagonal entries decay less quickly as their distance
from the diagonal increases.

Model 7:

A Toeplitz covariance matrix with alternating signs, where

$\psi^{(jl)}=\begin{cases}1,&j=l\\\ (-1)^{\lvert j-l\rvert}0.6\lvert
j-l\rvert^{-1.3},&\text{otherwise}\end{cases}.$

This model is almost identical to Model 6, though the signs of the covariance
matrix’s entries are alternating.

Model 8:

A covariance matrix inspired by the latent variable model described in
Equation (S3). Let $\boldsymbol{\beta}_{J\times
3}~{}=~{}(\beta_{1},\ldots,\beta_{J})^{\top}$, where $\beta_{j}$ are randomly
generated using a $N(0,\mathbf{I}_{3\times 3})$ distribution for
$j=1,\ldots,J$. Then
$\boldsymbol{\psi}=\boldsymbol{\beta}\boldsymbol{\beta}^{\top}+\mathbf{I}_{J\times
J}$ is the covariance matrix of a model with three latent factors.

Each covariance model was used to generate datasets consisting of
$n\in\\{50,100,200,500\\}$ i.i.d. multivariate Gaussian, mean-zero
observations. For each model and sample size, five data dimension ratios were
considered: $J/n\in\\{0.3,0.5,1,2,5\\}$. Together, the eight covariance
models, four sample sizes, and five dimensionality ratios result in $160$
distinct simulation settings. For each such setting, the performance of the
cross-validated selector with respect to the various oracle selectors and
several well-established estimators is evaluated based on aggregation across
$200$ Monte Carlo repetitions.

We applied our estimator selection procedure, which we refer to as cvCovEst,
using a 5-fold CV scheme. The library of candidate estimators is provided in
Table 1, which includes details on these estimators’ possible hyperparameters.
Seventy-four estimators make up the library of candidates. The considered
classes of covariance matrix estimators are briefly introduced in Section S1
of the supplement.

We also compared the estimates generated by our method against those of the
individual candidate procedures using the simulated datasets. Our CV scheme
was used to select hyperparameters of these competing approaches where
necessary. As stated previously, our CV approach is generally equivalent to
these estimators’ hyperparameter selection procedures. The hyperparameters
considered are provided in Table S1. Where appropriate, the competing methods’
hyperparameters are more varied than those used by cvCovEst, reflecting more
aggressive estimation procedures that one might employ when only using a
single family of estimators.

Table 1: Families of candidate estimators used by cvCovEst in the simulation study ($74$ distinct estimators in total). Estimator | Hyperparameters
---|---
Sample covariance matrix | Not applicable
Hard thresholding (Bickel and Levina, 2008b) | Thresholds $=\\{0.1,0.2,\ldots,1.0\\}$
SCAD thresholding (Fan and Li, 2001; Rothman et al., 2009) | Thresholds $=\\{0.1,0.2,\ldots,1.0\\}$
Adaptive LASSO (Rothman et al., 2009) | Thresholds $=\\{0.1,0.2,\ldots,0.5\\}$; exponential weights $=\\{0.1,0.2,\ldots,0.5\\}$
Banding (Bickel and Levina, 2008a) | Bands $=\\{1,2,\ldots,5\\}$
Tapering (Cai et al., 2010) | Bands $=\\{2,4,\ldots,10\\}$
Linear shrinkage (Ledoit and Wolf, 2004) | Not applicable
Dense linear shrinkage (Schäfer and Strimmer, 2005) | Not applicable
Nonlinear shrinkage (Ledoit and Wolf, 2018) | Not applicable
POET (Fan et al., 2013) using hard thresholding | Latent factors $=\\{1,2,\ldots,5\\}$; thresholds $=\\{0.1,0.2,0.3\\}$

###  4.2 Simulation Study Results

Figure 1: Comparison of the cross-validated selection and cross-validated
oracle selection’s mean cross-validated conditional risk differences. In all
cases, the mean risk difference ratios rapidly converge to one, even when the
true data-generating process possesses a dense covariance matrix. Note the
differing y-axis scales for the different models.

To examine empirically the optimality results of Theorem 1, we computed
analytically, for each replication, the cross-validated conditional risk
differences of the cross-validated selection, $\hat{k}_{p_{n},n}$, and the
cross-validated oracle selection, $\tilde{k}_{p_{n},n}$, using as candidates
the estimators in Table 1. We note that no penalty is attributed to estimators
generating rank-deficient estimates, like the sample covariance matrix when
$J>n$, though this limitation is generally of practical importance. When the
situation dictates that the resulting estimate must be positive-definite, the
library of candidates should be assembled accordingly.

From the cross-validated conditional risk differences, the Monte Carlo
expectations of the risk differences stratified by $n$, $J/n$, and the
covariance model were computed. The ratios of the expected risk differences
are presented in Figure 1. These results make clear that, for virtually all
models considered, the estimator chosen by our proposed cross-validated
selection procedure has a risk difference asymptotically identical on average
to that of the cross-validated oracle.

Figure 2:  Comparison of the cross-validated selection and oracle selection’s
full-dataset conditional risk differences. Asymptotic optimality is nearly
achieved for small to moderate sample sizes and numbers of features in all
considered models. Note the differing y-axis scales for the different models.
Figure 3:  Comparison of competing, bespoke covariance matrix estimation
procedures to our cross-validated selection approach in terms of the Monte
Carlo mean Frobenius norm under a variety of data-generating processes. Note
that the scales of the y-axis are tailored to the covariance matrix model.

A stronger result, corresponding to Equation (13) of Theorem 1, is presented
in Figure S1. For all but Models 1 and 8, we find that our algorithm’s
selection is virtually equivalent to the cross-validated oracle selection for
$n\geq 200$ and $J/n\geq 0.5$. Even for Model 8, in which the covariance
matrices are more difficult to estimate due to their dense structures, we find
that our selector identifies the optimal estimator with probability tending to
$1$ for $n\geq 200$ and $J/n=5$.

More impressive still are the results presented in Figure 2 that characterize
the full-dataset conditional risk difference ratios. For all covariance matrix
models considered, with the exception of Model 1, our procedure’s selections
attain near asymptotic optimality for moderate values of $n$ and $J/n$. This
suggests that our loss-based estimator selection approach’s theoretical
guarantee, as outlined in Corollary 1, is achievable in many practical
settings.

In addition to verifying our method’s asymptotic behavior, we compared its
estimates against those of competing approaches. We computed the Frobenius
norm of each procedure’s estimate against the corresponding true covariance
matrix. The mean norms over all simulations were then computed for each
covariance matrix estimation procedure, again stratified by $n$, $J/n$, and
the covariance matrix model (Figure 3). We repeat this benchmarking experiment
under the spectral norm even though our theoretical results do not relate to
it. We nonetheless expect similar outcomes given the relationship between
these two norms for reasons previously described. The accuracy of our
estimation procedure with respect to the leading eigenvalue of the covariance
matrix is presented in Figure S5.

Altogether, the results of these simulation experiments speak plainly. Our
proposed cross-validated loss-based covariance matrix estimator selection
procedure should be employed when estimating covariance matrices in high-
dimensional settings. An estimator selected by our procedure tends to perform
no worse than the single best estimator contained in the library of
candidates, while avoiding the need for making subjective (i.e., possibly
unwarranted) choices among estimator families. By letting the data inform the
choice of covariance estimator, significant gains are to be expected relative
to even informed choices among popular or state-of-the-art estimators.

##  5 Real Data Examples

Single-cell transcriptome sequencing (scRNA-seq) has taken the biological
sciences by storm, allowing researchers to study the gene expression profiles
of individual cells. The fine-grained transcriptomic data that it provides
have been used to identify rare cell populations and to elucidate the
developmental relationships between diverse cellular states. These
developments will likely lead to new and exciting clinical applications in the
near future (Haque et al., 2017).

Given that a typical scRNA-seq dataset possesses tens of thousands of features
(i.e., genes), most workflows prescribe a dimensionality reduction step
following data preprocessing. In addition to reducing the amount of
computational resources needed to analyze the data, reducing dimensionality
often has the effect of mitigating noise corrupting the biological signal of
interest. The lower-dimensional embedding may then subsequently be used for
downstream analysis, such as novel cell-type discovery via clustering.

One of the most popular methods used to reduce the dimensionality of scRNA-seq
data is uniform manifold approximation and projection (UMAP) (McInnes et al.,
2018). This method captures the most salient non-linear relationships among a
high-dimensional dataset’s features and projects them onto a reduced-
dimensional space. Instead of applying UMAP directly, the scRNA-seq dataset’s
leading principal components (PCs) are often used as an initialization. This
initial dimensionality reduction by PCA is believed to play a helpful role in
denoising.

However, PCA typically relies on the sample covariance matrix, and so when the
dataset is high-dimensional, the resulting principal components are known to
be poor estimates of those of the population (Johnstone and Lu, 2009). We
hence posit that our cross-validated estimator selection procedure could form
a basis for an improved PCA. That is, we hope that the eigenvectors resulting
from the eigendecomposition of our estimated covariance matrix could be used
to generate a set of estimates closer to the true PCs in terms of risk. These
PCs could then be used as initializations for UMAP, giving rise to low-
dimensional embeddings that better reflect a dataset’s biological signal.
Indeed, additional simulation results provided in Section S5 suggest that
cvCovEst produces estimates of the leading eigenvalue at least as well as
those produced by the sample covariance matrix, in terms of the spectral norm.

To verify this conjecture, we applied our procedure to two scRNA-seq datasets
for which the cell types are known a priori. These data were obtained from the
scRNAseq Bioconductor R package (Risso and Cole, 2020), and prepared for
analysis using a workflow outlined in Amezquita et al. (2020). A 5-fold CV
scheme was used; the library of candidate estimators is provided in Table S2.
We expect that cells of the same class will form tight, distinct clusters
within the low-dimensional representations. The resulting embeddings, which we
refer to as the cvCovEst-based embeddings, were then compared to those
produced by UMAP using traditional PCA for initialization, which we refer to
as the PCA-based embeddings. For each embedding, the 20 leading PCs were fed
to UMAP. The first dataset is a collection of 285 mouse visual cortex cells
(Tasic et al., 2016), and the second dataset consists of 2,816 mouse brain
cells (Zeisel et al., 2015). The ratios of features to observations, 3.51 and
0.36 for the Tasic et al. (2016) and Zeisel et al. (2015) datasets,
respectively, fall closely in line with those considered in the simulation
study presented in Section 4. The 1,000 most variable genes of each dataset
were used to compute the PCs of both embeddings.

Figure 4:  Comparisons of scRNA-seq datasets’ UMAP embeddings based on vanilla
PCA or PCA with the cross-validated selection’s covariance matrix estimate,
for the 1,000 most variable genes. The datasets consist of (A) 285 cells
collected from the visual cortex of mice and (B) 2,816 mouse brain cells. The
cvCovEst-based PCA embedding improves upon that based on traditional PCA with
respect to the average silhouette width in both datasets.

The resulting UMAP plots are presented in Figure 4. Though the two embeddings
generated for each dataset are qualitatively similar, the low-dimensional
representations relying on our loss-based approach are refined. In Figure 4A,
there are a number of cells that are erroneously clustered with the wrong cell
types in the PCA-based embedding, but not in the cvCovEst-based embedding.
This is further captured by the 41% increase in average silhouette width of
our method relative to the traditional approach. The improvements are less
pronounced in Figure 4B: some cells are more densely clustered in the
cvCovEst-based embedding and there are fewer mis-clustered cells, but the
average silhouette width is only marginally better.

That the relative improvement of the cvCovEst-based approach over the
traditional procedure is more striking in the Tasic et al. (2016) dataset
should not be surprising. The estimation problem is truly high-dimensional,
whereas this is not the case for the Zeisel et al. (2015) dataset. As the
ratio of the number of features to observations decreases, there is simply
less to gain from the use of estimators other than the sample covariance
matrix. Ideally, data-adaptive selection procedures should be cognizant of
this. Indeed, cvCovEst, when applied to the Zeisel et al. (2015) dataset,
selects an estimator whose cross-validated empirical risk is only slightly
smaller than that of the sample covariance matrix, and whose leading PCs are
virtually identical. Alternatively, it may be that the underlying, population
covariance matrix is a member of the class of parameters described in Bunea
and Xiao (2015). The diagnostics of the datasets’ corresponding covariance
matrix estimation procedures are presented in Figures S2 and S3 to provide
further intuition.

These results show that our loss-based covariance matrix estimation procedure
can improve the resolution of scRNA-seq data analyses while promoting
objectivity. We expect that similar outcomes are possible when including our
method in other procedures used to analyze high-dimensional data generated by
similar high-throughput biological assays.

##  6 Discussion

This work extends Dudoit and van der Laan (2003)’s framework for
asymptotically optimal, data-adaptive estimator selection to the problem of
covariance matrix estimation in high-dimensional settings. We provide
sufficient conditions under which our cross-validated procedure is
asymptotically optimal in terms of risk, and show that it generalizes the
cross-validated hyperparameter selection procedures employed by existing
estimation approaches. We further demonstrate its applicability to a wide
range of data-generating processes and estimators in both comprehensive
simulation experiments and two real-world data analysis examples.

The simulation study provides evidence that near-optimal results are achieved
in datasets with relatively modest numbers of observations and many features
across models indexed by diverse covariance matrix structures. These results
also establish that our cross-validated procedure performs as well as the best
bespoke estimation procedure in a variety of settings. Our scRNA-seq data
examples further illustrate the utility of our approach in fields where high-
dimensional data are collected routinely.

Practitioners need no longer rely upon intuition alone when deciding which
candidate estimator is best from among a library of diverse estimators. We
expect that a variety of analytical tasks relying upon the accurate estimation
of the covariance matrix beyond the exploratory analyses considered here, like
clustering, hypothesis testing, and latent variable estimation, stand to
benefit from the application of this framework.

## Acknowledgements

We thank Brian Collica and Jamarcus Liu for significant contributions to the
development of the cvCovEst R package.

## References

  * Amezquita et al. [2020] R. A. Amezquita, A. T. L. Lun, E. Becht, V. J. Carey, L. N. Carpp, L. Geistlinger, F. Marini, K. Rue-Albrecht, D. Risso, C. Soneson, L. Waldron, H. Pagès, M. L. Smith, W. Huber, M. Morgan, R. Gottardo, and S. C. Hicks.  Orchestrating single-cell analysis with bioconductor.  _Nature Methods_ , 17(2):137–145, 2020.  doi: 10.1038/s41592-019-0654-x.  URL https://doi.org/10.1038/s41592-019-0654-x.
  * Bai [2003] J. Bai.  Inferential theory for factor models of large dimensions.  _Econometrica_ , 71(1):135–171, 2003.  doi: 10.1111/1468-0262.00392.
  * Bennett [1962] G. Bennett.  Probability inequalities for the sum of independent random variables.  _Journal of the American Statistical Association_ , 57(297):33–45, 1962.
  * Bernstein [1946] S. N. Bernstein.  The theory of probabilities.  1946\.
  * Bickel and Levina [2008a] P. J. Bickel and E. Levina.  Regularized estimation of large covariance matrices.  _Annals of Statistics_ , 36(1):199–227, 2008a.  URL http://www.jstor.org/stable/25464621.
  * Bickel and Levina [2008b] P. J. Bickel and E. Levina.  Covariance regularization by thresholding.  _Annals of Statistics_ , 36(6):2577–2604, 2008b.  URL http://www.jstor.org/stable/25464728.
  * Boileau et al. [2021] P. Boileau, N. S. Hejazi, B. Collica, M. J. van der Laan, and S. Dudoit.  ‘cvcovest‘: Cross-validated covariance matrix estimator selection and evaluation in ‘r‘.  _Journal of Open Source Software_ , 6(63):3273, 2021.  doi: 10.21105/joss.03273.  URL https://doi.org/10.21105/joss.03273.
  * Breiman and Spector [1992] L. Breiman and P. Spector.  Submodel selection and evaluation in regression. the X-random case.  _International Statistical Review_ , pages 291–319, 1992.
  * Bunea and Xiao [2015] F. Bunea and L. Xiao.  On the sample covariance matrix estimator of reduced effective rank population matrices, with applications to fPCA.  _Bernoulli_ , 21(2):1200 – 1230, 2015.  doi: 10.3150/14-BEJ602.  URL https://doi.org/10.3150/14-BEJ602.
  * Cai and Liu [2011] T. Cai and W. Liu.  Adaptive thresholding for sparse covariance matrix estimation.  _Journal of the American Statistical Association_ , 106(494):672–684, 2011.  doi: 10.1198/jasa.2011.tm10560.  URL https://doi.org/10.1198/jasa.2011.tm10560.
  * Cai et al. [2010] T. T. Cai, C.-H. Zhang, and H. H. Zhou.  Optimal rates of convergence for covariance matrix estimation.  _Annals of Statistics_ , 38(4):2118–2144, 08 2010\.  doi: 10.1214/09-AOS752.
  * Coyle and Hejazi [2018] J. Coyle and N. Hejazi.  origami: A generalized framework for cross-validation in R.  _Journal of Open Source Software_ , 3(21):512, 2018.  doi: 10.21105/joss.00512.  URL https://doi.org/10.21105/joss.00512.
  * Dudoit and van der Laan [2003] S. Dudoit and M. J. van der Laan.  Asymptotics of cross-validated risk estimation in estimator selection and performance assessment.  Working Paper 126, University of California, Berkeley, Berkeley, 2003\.  URL https://biostats.bepress.com/ucbbiostat/paper126.
  * Dudoit and van der Laan [2005] S. Dudoit and M. J. van der Laan.  Asymptotics of cross-validated risk estimation in estimator selection and performance assessment.  _Statistical Methodology_ , 2(2):131–154, 2005\.
  * Efron [2012] B. Efron.  _Large-scale inference: empirical Bayes methods for estimation, testing, and prediction_.  Cambridge University Press, 2012.
  * Fan and Li [2001] J. Fan and R. Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ , 96(456):1348–1360, 2001.  doi: 10.1198/016214501753382273.
  * Fan et al. [2008] J. Fan, Y. Fan, and J. Lv.  High dimensional covariance matrix estimation using a factor model.  _Journal of Econometrics_ , 147(1):186 – 197, 2008.  doi: https://doi.org/10.1016/j.jeconom.2008.09.017.  URL http://www.sciencedirect.com/science/article/pii/S0304407608001346.  Econometric modelling in finance and risk management: An overview.
  * Fan et al. [2013] J. Fan, Y. Liao, and M. Mincheva.  Large covariance estimation by thresholding principal orthogonal complements.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 75(4):603–680, 8 2013.  doi: 10.1111/rssb.12016.
  * Fan et al. [2016] J. Fan, Y. Liao, and W. Wang.  Projected principal component analysis in factor models.  _Annals of Statistics_ , 44(1):219–254, 02 2016\.  doi: 10.1214/15-AOS1364.
  * Fan et al. [2019] J. Fan, W. Wang, and Y. Zhong.  Robust covariance estimation for approximate factor models.  _Journal of Econometrics_ , 208(1):5 – 22, 2019\.  doi: https://doi.org/10.1016/j.jeconom.2018.09.003.  Special Issue on Financial Engineering and Risk Management.
  * Fang et al. [2016] Y. Fang, B. Wang, and Y. Feng.  Tuning-parameter selection in regularized estimations of large covariance matrices.  _Journal of Statistical Computation and Simulation_ , 86(3):494–509, 2016.  doi: 10.1080/00949655.2015.1017823.  URL https://doi.org/10.1080/00949655.2015.1017823.
  * Friedman et al. [2001] J. Friedman, T. Hastie, and R. Tibshirani.  _The elements of statistical learning_ , volume 1.  Springer series in statistics New York, 2001.
  * Golub and Van Loan [1996] G. H. Golub and C. F. Van Loan.  _Matrix computations_.  Johns Hopkins Universtiy Press, 3rd edition, 1996.
  * Haque et al. [2017] A. Haque, J. Engel, S. A. Teichmann, and T. Lönnberg.  A practical guide to single-cell rna-sequencing for biomedical research and clinical applications.  _Genome Medicine_ , 9(1):75, 2017.  doi: 10.1186/s13073-017-0467-4.
  * Johnstone [2001] I. M. Johnstone.  On the distribution of the largest eigenvalue in principal components analysis.  _Annals of Statistics_ , pages 295–327, 2001.
  * Johnstone and Lu [2009] I. M. Johnstone and A. Y. Lu.  On consistency and sparsity for principal components analysis in high dimensions.  _Journal of the American Statistical Association_ , 104(486):682–693, 2009.  doi: 10.1198/jasa.2009.0121.  PMID: 20617121.
  * Lam and Fan [2009] C. Lam and J. Fan.  Sparsistency and rates of convergence in large covariance matrix estimation.  _Annals of Statistics_ , 37(6B):4254–4278, 12 2009.  doi: 10.1214/09-AOS720.
  * Ledoit and Wolf [2004] O. Ledoit and M. Wolf.  A well-conditioned estimator for large-dimensional covariance matrices.  _Journal of Multivariate Analysis_ , 88(2):365 – 411, 2004.  doi: https://doi.org/10.1016/S0047-259X(03)00096-4.
  * Ledoit and Wolf [2012] O. Ledoit and M. Wolf.  Nonlinear shrinkage estimation of large-dimensional covariance matrices.  _Annals of Statistics_ , 40(2):1024–1060, 04 2012\.  doi: 10.1214/12-AOS989.
  * Ledoit and Wolf [2015] O. Ledoit and M. Wolf.  Spectrum estimation: A unified framework for covariance matrix estimation and PCA in large dimensions.  _Journal of Multivariate Analysis_ , 139:360 – 384, 2015\.  ISSN 0047-259X.  doi: https://doi.org/10.1016/j.jmva.2015.04.006.
  * Ledoit and Wolf [2018] O. Ledoit and M. Wolf.  Analytical nonlinear shrinkage of large-dimensional covariance matrices.  ECON - Working Papers 264, Department of Economics - University of Zurich, 2018.  URL https://EconPapers.repec.org/RePEc:zur:econwp:264.
  * McInnes et al. [2018] L. McInnes, J. Healy, and J. Melville.  UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction.  2 2018.  URL http://arxiv.org/abs/1802.03426.
  * Onatski [2012] A. Onatski.  Asymptotics of the principal components estimator of large factor models with weakly influential factors.  _Journal of Econometrics_ , 168(2):244 – 258, 2012.  doi: https://doi.org/10.1016/j.jeconom.2012.01.034.
  * R Core Team [2021] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Risso and Cole [2020] D. Risso and M. Cole.  _scRNAseq: Collection of Public Single-Cell RNA-Seq Datasets_ , 2020\.  R package version 2.3.17.
  * Robbins [1964] H. Robbins.  The empirical bayes approach to statistical decision problems.  _Annals of Mathematical Statistics_ , 35(1):1–20, 1964.
  * Rothman et al. [2009] A. J. Rothman, E. Levina, and J. Zhu.  Generalized thresholding of large covariance matrices.  _Journal of the American Statistical Association_ , 104(485):177–186, 2009.  doi: 10.1198/jasa.2009.0101.  URL https://doi.org/10.1198/jasa.2009.0101.
  * Schäfer and Strimmer [2005] J. Schäfer and K. Strimmer.  A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics.  _Statistical Applications in Genetics and Molecular Biology_ , 4(1), 2005.  doi: https://doi.org/10.2202/1544-6115.1175.
  * Stock and Watson [2002] J. H. Stock and M. W. Watson.  Forecasting using principal components from a large number of predictors.  _Journal of the American Statistical Association_ , 97(460):1167–1179, 2002.  URL http://www.jstor.org/stable/3085839.
  * Tasic et al. [2016] B. Tasic, V. Menon, T. N. Nguyen, T. K. Kim, T. Jarsky, Z. Yao, B. Levi, L. T. Gray, S. A. Sorensen, T. Dolbeare, D. Bertagnolli, J. Goldy, N. Shapovalova, S. Parry, C. Lee, K. Smith, A. Bernard, L. Madisen, S. M. Sunkin, M. Hawrylycz, C. Koch, and H. Zeng.  Adult mouse cortical cell taxonomy revealed by single cell transcriptomics.  _Nature Neuroscience_ , 19(2):335–346, 2016.  doi: 10.1038/nn.4216.  URL https://doi.org/10.1038/nn.4216.
  * van der Laan and Dudoit [2003] M. J. van der Laan and S. Dudoit.  Unified cross-validation methodology for selection among estimators and a general cross-validated adaptive epsilon-net estimator: Finite sample oracle inequalities and examples.  Working Paper 130, University of California, Berkeley, Berkeley, 2003\.  URL https://biostats.bepress.com/ucbbiostat/paper130/.
  * van der Vaart et al. [2006] A. W. van der Vaart, S. Dudoit, and M. J. van der Laan.  Oracle inequalities for multi-fold cross validation.  _Statistics and Decisions_ , 24:351–371, 2006.
  * Zeisel et al. [2015] A. Zeisel, A. B. Muñoz-Manchado, S. Codeluppi, P. Lönnerberg, G. La Manno, A. Juréus, S. Marques, H. Munguba, L. He, C. Betsholtz, C. Rolny, G. Castelo-Branco, J. Hjerling-Leffler, and S. Linnarsson.  Cell types in the mouse cortex and hippocampus revealed by single-cell rna-seq.  _Science_ , 347(6226):1138–1142, 2015.  ISSN 0036-8075.  doi: 10.1126/science.aaa1934.

SUPPLEMENTARY MATERIAL

##  S1 Candidate Covariance Matrix Estimators

Using the proposed cross-validated selection procedure effectively requires a
large and diverse set of candidate covariance matrix estimators. In this
spirit, we provide in the sequel a brief overview of select covariance matrix
estimators that have proven to be effective in a variety of settings. Note,
however, that the proposed selection framework need not be limited to those
described here. Thorough reviews of estimators have been provided by
pourahmadi:2013, fan:2016, and ledoit:2020a, to name only a few.

###  S1.1 Thresholding Estimators

An often natural, simplifying assumption about the true covariance matrix’s
structure is that it is sparse, that is, a non-negligible portion of its off-
diagonal elements have a value near zero or equal to zero. This assumption is
not altogether unreasonable: Given a system of numerous variables, it seems
unlikely in many settings that a majority of these variables would depend on
one another.

The class of generalized thresholding estimators bickelb:2008,
rothman:2009,cai:2011 is one collection of covariance matrix estimators based
upon this structural assumption. Given the sample covariance matrix
$\mathbf{S}_{n}$, the generalized thresholding estimator is defined as

$\hat{\Psi}_{\text{gt}}(P_{n};t)\coloneqq\left\\{t\left(S_{n}^{(jl)}\right):\;j,l=1,\ldots,J\right\\},$
(S1)

where $t(\cdot)$ is a thresholding function that often requires one or more
hyperparameters dictating the amount of regularization. The hard, soft,
smoothly clipped absolute deviation (SCAD) fan:2001, and adaptive LASSO
rothman:2009 functions are among the collection of suitable thresholding
operators. As an example, the hard-thresholding function is defined as
$t_{u}(S_{n}^{(jl)})\equiv S_{n}^{(jl)}I(S_{n}^{(jl)}>u)$ for some threshold
$u>0$. cai:2011 also demonstrated that element-specific thresholding functions
might be useful when the features’ variances are highly variable. The
hyperparameters for any specific thresholding function are often selected via
CV. Regardless of the choice of $t(\cdot)$, these estimators preserve the
symmetry of the sample covariance matrix and are invariant under permutations
of the features’ order.

bickelb:2008 and rothman:2009 have shown that these estimators are consistent
under the spectral norm (which is defined as a square matrix’s largest
absolute eigenvalue), assuming that $\log J/n\rightarrow 0$, that the
observations’ marginal distributions satisfy a tail condition, and that the
true covariance matrix is a member of the class of matrices satisfying a
particular notion of “approximate sparsity.” cai:2011 have derived similar
results for their entry-specific thresholding estimator over an even broader
parameter space of sparse covariance matrices.

###  S1.2 Banding and Tapering Estimators

A family of estimators related to thresholding estimators are banding and
tapering estimators bickela:2008,cai:2010. Like thresholding estimators, these
estimators rely on the assumption that the true covariance matrix is sparse.
However, the structural assumption of these estimators on
$\boldsymbol{\psi}_{0}$ is much more rigid than that of thresholding
estimators. Specifically, such estimators assume that the true covariance
matrix is a band matrix, that is, a sparse matrix whose non-zero entries are
concentrated about the diagonal. These estimators therefore require a natural
ordering of the observations’ features, operating under the assumption that
“distant” variables are uncorrelated. Such structure is often encountered in
longitudinal and time-series data.

Given the sample covariance matrix $\mathbf{S}_{n}$, the banding estimator of
bickela:2008 is defined as

$\hat{\Psi}_{\text{band}}(P_{n};b)\coloneqq\left\\{S_{n}^{(jl)}\mathbb{I}(\lvert
j-l\rvert\leq b):\;j,l=1,\ldots,J\right\\},$

where $b$ is a hyperparameter that determines the number of bands to retain
from the sample covariance matrix and is chosen via a CV procedure. For the
class of “bandable” covariance matrices, i.e., the set of well-conditioned
matrices whose elements not in the central bands of the matrix decay rapidly,
this banding estimator has been shown to be uniformly consistent in the
spectral norm so long as $\log(J)/n\rightarrow 0$.

The tapering estimator of cai:2010 is the smooth generalization of the banding
estimator, gradually shrinking the off-diagonal bands of the sample covariance
matrix towards zero. It is defined as

$\hat{\Psi}_{\text{tap}}(P_{n};b)\coloneqq\mathbf{W}_{b}\circ\mathbf{S}_{n},\\\
$

for some weight matrix $\mathbf{W}_{b}$. Here, “$\circ$” denotes the Hadamard
(element-wise) matrix product. Clearly, letting
$W^{(jl)}_{b}=\mathbb{I}(\lvert j-l\rvert\leq b)$ for some positive integer
$b$ results in the banding estimator. A popular weighting scheme derived by
cai:2010 is

$W^{(jl)}_{b}\coloneqq\begin{cases}1,&\text{when $\lvert
j-l\rvert\leq\frac{b}{2}$}\\\ 2-\frac{\lvert j-l\rvert}{b},&\text{when
$\frac{b}{2}<\lvert j-l\rvert\leq b$}\\\ 0,&\text{otherwise}\end{cases},$

which we use in our simulation study presented in Section 4. cai:2010 also
derived the optimal rates of convergence for this estimator under the
Frobenius and spectral norms, considering a class of bandable covariance
matrices that is more general than that considered by bickela:2008: The
smallest eigenvalue of the covariance matrices in this class can take on a
value of zero. However, this estimator does not improve upon the bounds set by
the banding estimator.

###  S1.3 Shrinkage Estimators

We next consider the linear and non-linear shrinkage estimators inspired by
Stein’s work on empirical Bayesian methods. These estimators are rotation-
equivariant, shrinking the eigenvalues of the sample covariance matrix towards
a set of target values, whilst leaving its eigenvectors untouched. In doing
so, the resultant estimators are better-conditioned than the sample covariance
matrix in a manner guaranteeing that the resultant covariance matrix estimator
be non-singular. Further, these estimators do not rely on sparsity assumptions
about the true covariance matrix, setting them apart from those previously
discussed.

One of the first shrinkage estimators, the linear shrinkage estimator of the
sample covariance matrix, was proposed by ledoit:2004. This estimator is
defined as the convex combination of the sample covariance matrix and the
identity matrix. Hence, it represents a compromise between $\mathbf{S}_{n}$,
an unbiased but highly variable estimator of $\boldsymbol{\psi}_{0}$ in high
dimensions, and $\mathbf{I}_{J\times J}$, a woefully biased but fixed
estimator. ledoit:2004 found that, under mild conditions, the asymptotically
optimal shrinkage intensity with respect to the scaled (by $J$) Frobenius norm
can be estimated consistently in high dimensions. This estimator is defined as

$\hat{\Psi}_{\text{ls}}(P_{n})\coloneqq\frac{b^{2}_{n}}{d^{2}_{n}}m_{n}\mathbf{I}+\frac{a_{n}^{2}}{d_{n}^{2}}\mathbf{S}_{n},$
(S2)

for $m_{n}=\text{tr}(\mathbf{S}_{n})/J$,
$d^{2}_{n}=\lVert\mathbf{S}_{n}-m_{n}\mathbf{I}\rVert^{2}_{F,1/J}$,
$\bar{b}^{2}_{n}=n^{-2}\sum_{i}\lVert
X_{i}X_{i}^{\top}-\mathbf{S}_{n}\rVert^{2}_{F,1/J}$,
$b^{2}_{n}=\min(\bar{b}^{2}_{n},d^{2}_{n})$, and
$a^{2}_{n}=d^{2}_{n}-b^{2}_{n}$.

Bespoke shrinkage targets may be used in place of the identity. For example,
one might consider a dense matrix target whose diagonal elements are the
average of the sample covariance matrix’s diagonal elements and whose off-
diagonal elements are equal to the average of all the sample covariance
matrix’s off-diagonal elements. For the sake of brevity, discussion of such
estimators is omitted, but examples are provided by, among others, ledoit:2003
and shafer:2005, particularly for use in financial economics and statistical
genomics applications, respectively.

When assumptions about the true covariance matrix’s structure are unfounded,
it can become impossible to select an appropriate linear shrinkage target.
Instead, one might consider generalizing these estimators to shrink the
eigenvalues of the sample covariance matrix in a non-linear fashion. That is,
an estimator that shrinks the sample covariance matrix’s eigenvalues not by a
common shrinkage factor (as with linear shrinkage estimators) but with
shrinkage factors tailored to each sample eigenvalue. As with the
aforementioned linear shrinkage estimators, such non-linear shrinkage
estimators produce positive-definite estimates so long as the shrunken sample
eigenvalues are positive and rotation-equivariant. These estimators belong to
a class initially introduced by stein:1986 and have since seen a resurgence in
the work of ledoit:2012,ledoit:2015. More recently, ledoit:2020b derived an
analytical non-linear shrinkage estimator that is asymptotically optimal in
high dimensions and more computationally efficient than their previously
formulated estimators.

###  S1.4 Estimators Based on Factor Models

Covariance matrix estimators based on factor models form another broad family
of estimators that do not assume sparsity of the true covariance matrix. Often
encountered in econometrics and finance, these estimators utilize the
operating assumption that the dataset’s observations are functions of a few
common, often latent, factors. The factor model can be described as follows:

$X=\mu+\boldsymbol{\beta}F+\epsilon,$ (S3)

where $X_{J\times 1}$ represents a random observation, $\mu_{J\times 1}$ a
mean vector, $\boldsymbol{\beta}_{J\times L}$ a matrix of factor coefficients,
$F_{L\times 1}$ a random vector of $L$ common factors, and $\epsilon_{J\times
1}$ a random error vector. Assuming that $F$ and $\epsilon$ are uncorrelated,
the covariance matrix of $X$ is given by

$\boldsymbol{\psi}=\boldsymbol{\beta}\text{Cov}(F)\boldsymbol{\beta}^{\top}+\boldsymbol{\psi}_{\epsilon},$
(S4)

where $\psi_{\epsilon}$ is the covariance matrix of the random error.

For a review on estimating the covariance matrix in the presence of observable
factors, see fan:2016. We now briefly discuss the estimation of
$\boldsymbol{\psi}$ when the factors are unobservable. Notice that when
$\boldsymbol{\psi}_{\epsilon}$ is not assumed to be diagonal, the
decomposition of $\boldsymbol{\psi}$ in Equation (S4) is not identifiable for
fixed $n$ and $J$, since the random vector $X$ is the only observed component
of the model. By letting $J\rightarrow\infty$, and assuming that the
eigenvalues of $\boldsymbol{\psi}_{\epsilon}$ are uniformly bounded or grow at
a slow rate relative to $J$ and that the eigenvalues of
$(1/J)\;\boldsymbol{\beta}^{\top}\boldsymbol{\beta}$ are uniformly bounded
away from zero and infinity, it can be shown that
$\boldsymbol{\beta}\text{Cov}(F)\boldsymbol{\beta}^{\top}$ is asymptotically
identifiable cai:2010. It follows from these assumptions that the signal in
the factors increases as the number of features increases, while the noise
contributed by the error term remains constant. The eigenvalues associated
with $\boldsymbol{\beta}\text{Cov}(F)\boldsymbol{\beta}^{\top}$ therefore
become easy to differentiate from those of $\boldsymbol{\psi}_{\epsilon}$.

Now, even with $\boldsymbol{\beta}\text{Cov}(F)\boldsymbol{\beta}^{\top}$
being asymptotically identifiable, $\boldsymbol{\beta}$ and $F$ cannot be
distinguished. As a solution, the following constraint is imposed on $F$:
$\text{Cov}(F)=\mathbf{I}_{L\times L}$. It then follows that

$\boldsymbol{\psi}=\boldsymbol{\beta}\boldsymbol{\beta}^{\top}+\boldsymbol{\psi}_{\epsilon}.$
(S5)

Under the additional assumption that the columns of $\boldsymbol{\beta}$ be
orthogonal, fan:2013 found that the leading $L$ eigenvalues of
$\boldsymbol{\psi}$ are spiked, meaning that they are bounded below by some
constant johnstone2001distribution, and grow at rate $O(J)$ as the dimension
of $\boldsymbol{\psi}$ increases. The remaining $J-L$ eigenvalues are then
either bounded above or grow slowly. This implies that the latent factors and
their loadings can be approximated via the eigenvalues and eigenvectors of
$\boldsymbol{\psi}$.

fan:2013 therefore proposed the Principal Orthogonal compleEment Thresholding
(POET) estimator of $\boldsymbol{\psi}$, which was motivated by the spectral
decomposition of the sample covariance matrix

$\begin{split}\mathbf{S}_{n}&=\sum_{j=1}^{J}\lambda_{j}V_{\cdot,j}V_{\cdot,j}^{\top}\\\
&\approx\sum_{j=1}^{L}\beta_{\cdot,j}\beta_{\cdot,j}^{\top}+\sum_{j=L+1}^{J}\lambda_{j}V_{\cdot,j}V_{\cdot,j}^{\top},\end{split}$

where $\lambda_{j}$ and $V_{\cdot,j}$ are the $j^{\text{th}}$ eigenvalues and
eigenvectors of $\mathbf{S}_{n}$, respectively, and $\beta_{\cdot,j}$ is the
$j^{\text{th}}$ column of $\boldsymbol{\beta}$. For ease of notation, we
denote the second term by $\mathbf{S}_{\epsilon}$ and refer to this matrix as
the principal orthogonal complement. The estimator for $L$ latent variables is
then defined as

$\hat{\Psi}_{\text{POET}}(P_{n};L,s)\coloneqq\sum_{j=1}^{L}\lambda_{j}V_{\cdot,j}V_{\cdot,j}^{\top}+\mathbf{T}_{\epsilon,s},$
(S6)

where $\mathbf{T}_{\epsilon,s}$ is the generalized thresholding matrix of
$\mathbf{S}_{\epsilon}$

$T_{\epsilon,s}^{(jl)}\coloneqq\begin{cases}S_{\epsilon}^{(jj)},&\text{when
$j=l$}\\\ s\left(S_{\epsilon}^{(jl)}\right),&\text{otherwise}\end{cases},$

for some generalized thresholding function $s$.

Although this estimator is computationally efficient, the assumptions encoding
the factor based model under which it is derived are such that the latent
features’ eigenvalues grow in $J$. This results in a poor convergence rate
under the spectral norm fan:2016.

##  S2 Proofs

###### Proof.

Proposition 1.

Assume without loss of generality that $\mathbb{E}_{P_{0}}[X]=0$ and
$\eta_{0}=1$. Then,

$\begin{split}\hat{\theta}_{p_{n},n}(k,1)&=\mathbb{E}_{B_{n}}\left[\frac{1}{np_{n}}\sum_{i=1}^{n}\mathbb{I}(B_{n}(i)=1)\lVert
X_{i}X^{\top}_{i}-\hat{\Psi}_{k}(P_{n,B_{n}}^{0})\rVert^{2}_{F,1}\right]\\\
&=\mathbb{E}_{B_{n}}\left[\frac{1}{np_{n}}\sum_{i=1}^{n}\mathbb{I}(B_{n}(i)=1)\sum_{j=1}^{J}\sum_{l=1}^{J}(X_{i}^{(j)}X^{(l)}_{i}-\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)})^{2}\right]\\\
&=\mathbb{E}_{B_{n}}\left[\frac{1}{np_{n}}\sum_{j=1}^{J}\sum_{l=1}^{J}\left(\sum_{\\{i:B_{n}(i)=1\\}}\left((X_{i}^{(j)}X^{(l)}_{i})^{2}-2X_{i}^{(j)}X^{(l)}_{i}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}\right)+(\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)})^{2}\right)\right]\\\
&=\mathbb{E}_{B_{n}}\Bigg{[}\sum_{j=1}^{J}\sum_{l=1}^{J}\bigg{(}(\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)})^{2}-2S(P_{n,B_{n}}^{1})^{(jl)}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}\\\
&\qquad\qquad\qquad\qquad\qquad+\frac{1}{np_{n}}\sum_{\\{i:B_{n}(i)=1\\}}(X_{i}^{(j)}X^{(l)}_{i})^{2}\bigg{)}\Bigg{]}\\\
&=\mathbb{E}_{B_{n}}\left[\sum_{j=1}^{J}\sum_{l=1}^{J}\bigg{(}\big{(}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}\big{)}^{2}-2S(P_{n,B_{n}}^{1})^{(jl)}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}\bigg{)}\right]+C_{1},\end{split}$
(S7)

where $C_{1}$ is constant with respect to $\hat{\Psi}_{k}(P_{n,B_{n}}^{0})$.
Now, notice that

$\begin{split}\hat{R}_{n}(\hat{\Psi}_{k},1)&=\mathbb{E}_{B_{n}}\left[\lVert\hat{\Psi}_{k}(P_{n,B_{n}}^{0})-\mathbf{S}_{n}(P_{n,B_{n}}^{1})\rVert^{2}_{F,1}\right]\\\
&=\mathbb{E}_{B_{n}}\left[\sum_{j=1}^{J}\sum_{l=1}^{J}\left((\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)})^{2}-2S(P_{n,B_{n}}^{1})^{(jl)}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}+(S(P_{n,B_{n}}^{1})^{(jl)})^{2}\right)\right]\\\
&=\mathbb{E}_{B_{n}}\left[\sum_{j=1}^{J}\sum_{l=1}^{J}\left((\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)})^{2}-2S(P_{n,B_{n}}^{1})^{(jl)}\hat{\Psi}_{k}(P_{n,B_{n}}^{0})^{(jl)}\right)\right]+C_{2},\end{split}$
(S8)

where $C_{2}$ is constant with respect to $\hat{\Psi}_{k}(P_{n,B_{n}}^{0})$.

Thus,

$\begin{split}\hat{k}_{p_{n},n}&=\text{arg
min}_{k\in\\{1,\ldots,K\\}}\hat{\theta}_{p_{n},n}(k,1)\\\ &=\text{arg
min}_{k\in\\{1,\ldots,K\\}}\hat{R}_{n}(\hat{\Psi}_{k},1).\end{split}$

∎

###### Proof.

Proposition 2. Let the hard-thresholding estimator minimizing the expected
cross-validated conditional risk under $P_{0}$ be indexed by $k_{0}$. Then

$\begin{split}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]&\leq\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(k_{0})-\theta_{0}]\\\
&=\mathbb{E}_{P_{0}}\left[\mathbb{E}_{B_{n}}\left[\mathbb{E}_{P_{0}}\left[\lVert\hat{\Psi}_{k_{0}}(P_{n,B_{n}}^{0})-XX^{\top}\rVert^{2}_{F,1}-\lVert\boldsymbol{\psi}_{0}-XX^{\top}\rVert^{2}_{F,1}\mid
P_{n,B_{n}}^{0}\right]\right]\right]\\\
&=\mathbb{E}_{P_{0}}\left[\mathbb{E}_{B_{n}}\left[\mathbb{E}_{P_{0}}\left[\lVert\hat{\Psi}_{k_{0}}(P_{n,B_{n}}^{0})-\boldsymbol{\psi}_{0}\rVert^{2}_{F,1}\mid
P_{n,B_{n}}^{0}\right]\right]\right]\\\ &=O(s(J)\log J).\end{split}$

The first inequality follows from the definition of the cross-validated
conditional risk under $P_{0}$; the last equality follows from Theorem 2 (and
its subsequent discussion) of Bickel and Levina (2008b), since $X$ is sub-
Gaussian by the boundedness condition of Assumption 1 and $J/n\rightarrow m$
as $n,J\rightarrow\infty$. Then,

$\frac{c(\delta,\overline{M}(J))(1+\log(K)}{np_{n}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}]}=\Omega\left(\frac{J}{s(J)\log
J}\right).$

∎

######  Lemma A1.

Let $Z_{k}\coloneqq
L(X;\hat{\Psi}_{k}(P^{0}_{n,B_{n}}))-L(X;\boldsymbol{\psi}_{0})$. Then,
assuming $\mathbb{E}_{P_{0}}[Z_{k}\mid P^{0}_{n,B_{n}},B_{n}]>0$,

$\text{Var}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]\leq\overline{M}(J)\mathbb{E}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right],$

where $\overline{M}(J)\coloneqq 4(M_{1}+M_{2})^{2}J^{2}$.

###### Proof.

$\begin{split}\mathbb{E}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]&=\mathbb{E}_{P_{0}}\left[L(X;\hat{\Psi}_{k}(P^{0}_{n,B_{n}}))-L(X;\boldsymbol{\psi}_{0})\Big{|}P^{0}_{n,B_{n}},B_{n}\right]\\\
&=\mathbb{E}_{P_{0}}\left[\sum_{j=1}^{J}\sum_{l=1}^{J}(X^{(j)}X^{(l)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}-(X^{(j)}X^{(l)}-\psi_{0}^{(jl)})^{2}\bigg{|}P_{n,B_{n}}^{0},B_{n}\right]\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}\mathbb{E}_{P_{0}}\bigg{[}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))(2X^{(j)}X^{(l)}-\psi_{0}^{(jl)}\\\
&\qquad\qquad\qquad\qquad\qquad-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))\bigg{|}P_{n,B_{n}}^{0},B_{n}\bigg{]}\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))\\\
&\qquad\qquad\qquad\qquad\qquad\mathbb{E}_{P_{0}}\bigg{[}2X^{(j)}X^{(l)}-\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0})\bigg{|}P_{n,B_{n}}^{0},B_{n}\bigg{]}\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}\\\
&=\sum_{j=1}^{J}\sum_{l=1}^{J}\mathbb{E}_{P_{0}}\left[(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}\Big{|}P_{n,B_{n}}^{0},B_{n}\right],\\\
\end{split}$

where the second to last equality follows by noting that
$\mathbb{E}_{P_{0}}[X^{(j)}X^{(l)}|P_{n,B_{n}}^{0},B_{n}]=\mathbb{E}_{P_{0}}[X^{(j)}X^{(l)}]$
and that, by definition, $\mathbb{E}_{P_{0}}[X^{(j)}X^{(l)}]=\psi_{0}^{(jl)}$.
Then,

$\begin{split}\text{Var}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]&\leq\mathbb{E}_{P_{0}}\left[Z_{k}^{2}\mid
P^{0}_{n,B_{n}},B_{n}\right]\\\
&=\mathbb{E}_{P_{0}}\Bigg{[}\bigg{(}\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))(2X^{(j)}X^{(l)}-\psi_{0}^{(jl)}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))\bigg{)}^{2}\Bigg{|}P^{0}_{n,B_{n}},B_{n}\Bigg{]}\\\
&=J^{4}\mathbb{E}_{P_{0}}\Bigg{[}\bigg{(}\frac{1}{J^{2}}\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))(2X^{(j)}X^{(l)}-\psi_{0}^{(jl)}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))\bigg{)}^{2}\Bigg{|}P^{0}_{n,B_{n}},B_{n}\Bigg{]}\\\
&\leq
J^{4}\mathbb{E}_{P_{0}}\Bigg{[}\frac{1}{J^{2}}\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}(2X^{(j)}X^{(l)}-\psi_{0}^{(jl)}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}\Bigg{|}P^{0}_{n,B_{n}},B_{n}\Bigg{]}\\\
&\leq
J^{2}4(M_{1}+M_{2})^{2}\mathbb{E}_{P_{0}}\left[\sum_{j=1}^{J}\sum_{l=1}^{J}(\psi_{0}^{(jl)}-\hat{\Psi}^{(jl)}_{k}(P_{n,B_{n}}^{0}))^{2}\Bigg{|}P^{0}_{n,B_{n}},B_{n}\right]\\\
&=4(M_{1}+M_{2})^{2}J^{2}\mathbb{E}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]\\\
&=\overline{M}(J)\mathbb{E}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right].\\\ \end{split}$

Here, the second inequality holds from the application of Jensen’s Inequality
to the square of the double sum, which is effectively an expectation of a
discrete uniform random variable when scaled by $J^{2}$. The final inequality
results from Assumptions 1 and 2, and concludes the proof. ∎

The following lemma is a result taken directly from laan-dudoit:2003. It is
restated here for convenience.

######  Lemma A2.

Let $Y_{1},Y_{2},\ldots$ be a sequence of random variables. If
$\mathbb{E}[\lvert Y_{n}\rvert]=O(g(n))$ for some positive function
$g(\cdot)$, then $Y_{n}=O_{P}(g(n))$.

###### Proof.

We must show that, for each $\epsilon>0$, there exists an $N$ and $B>0$ such
that $\mathbb{P}(\lvert Y_{n}\rvert/g(n)>B)<\epsilon$ for all $n>N$. By
assumption, there exists an $N$ and a $C>0$ such that $\mathbb{E}[\lvert
Y_{n}\rvert]/g(n)<C$ for all $n>N$. By defining $C/B=\epsilon$ and making use
of Markov’s Inequality, we find that

$\mathbb{P}\left(\frac{\lvert
Y_{n}\rvert}{g(n)}>B\right)\leq\frac{\mathbb{E}\left[\lvert
Y_{n}\rvert\right]}{g(n)B}\leq\frac{C}{B}=\epsilon.$

∎

Having derived Lemma A1, the remainder of the proof for Theorem 1 closely
follows the proof of the first theorem in Dudoit and van der Laan (2005).

###### Proof.

Theorem 1, Finite-Sample Result.

$\begin{split}0&\leq\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}\\\
&=\mathbb{E}_{B_{n}}\bigg{[}\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP_{0}(x)\\\
&\qquad\qquad-(1+\delta)\int
L((x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\\\
&\qquad\qquad+(1+\delta)\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\bigg{]}\\\
&\leq\mathbb{E}_{B_{n}}\bigg{[}\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP_{0}(x)\\\
&\qquad\qquad-(1+\delta)\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\\\
&\qquad\qquad+(1+\delta)\int(L(x;\hat{\Psi}_{\tilde{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\bigg{]}\\\
&=\mathbb{E}_{B_{n}}\bigg{[}\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\psi_{0}))dP_{0}(x)\\\
&\qquad\qquad-(1+\delta)\int(L(x;\hat{\Psi}_{\hat{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\\\
&\qquad\qquad+(1+\delta)\int(L(x;\hat{\Psi}_{\tilde{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\\\
&\qquad\qquad-(1+2\delta)\int(L(x;\hat{\Psi}_{\tilde{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP_{0}(x)\\\
&\qquad\qquad+(1+2\delta)\int(L(x;\hat{\Psi}_{\tilde{k}_{p_{n},n}}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP_{0}(x)\bigg{]}.\\\
\end{split}$ (S9)

The first inequality is by assumption, and the second is by definition of
$\hat{k}_{p_{n},n}$ such that
$\hat{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})\leq\hat{\theta}_{p_{n},n}(k)$
$\forall k$. For simplicity in the remainder of the proof, we replace
$\hat{k}_{p_{n},n}$ and $\tilde{k}_{p_{n},n}$ with $\hat{k}$ and $\tilde{k}$,
respectively, in a slight abuse of notation.

Now, let the first two terms of the last expression in Equation (S9) be
denoted by $R_{\hat{k},n}$, and the third and fourth terms by
$T_{\tilde{k},n}$. The last term is the cross-validated oracle’s risk
difference: $(1+2\delta)(\tilde{\theta}_{p_{n},n}(\tilde{k})-\theta_{0})$.
Thus,

$0\leq\tilde{\theta}_{p_{n},n}(\hat{k})-\theta_{0}\leq(1+2\delta)(\tilde{\theta}_{p_{n},n}(\tilde{k})-\theta_{0})+R_{\hat{k},n}+T_{\tilde{k},n}.$
(S10)

We next show that
$\mathbb{E}_{P_{0}}[R_{\hat{k},n}+T_{\tilde{k},n}]~{}\leq~{}2c(\delta,\overline{M}(J))(1+\text{log}(K))/(np_{n})$,
where $c(\delta,\overline{M}(J))=2(1+\delta)^{2}\overline{M}(J)(1/\delta+1/3)$
for some $\delta>0$. For convenience, let

$\begin{split}\hat{H}_{k}&\coloneqq\int(L(x;\hat{\Psi}_{k}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP^{1}_{n,B_{n}}(x)\\\
\tilde{H}_{k}&\coloneqq\int(L(x;\hat{\Psi}_{k}(P_{n,B_{n}}^{0}))-L(x;\boldsymbol{\psi}_{0}))dP_{0}(x)\\\
R_{k,n}(B_{n})&\coloneqq(1+\delta)(\tilde{H}_{k}-\hat{H}_{k})-\delta\tilde{H}_{k}\\\
T_{k,n}(B_{n})&\coloneqq(1+\delta)(\hat{H}_{k}-\tilde{H}_{k})-\delta\tilde{H}_{k},\end{split}$

so that $R_{k,n}=\mathbb{E}_{B_{n}}[R_{k,n}(B_{n})]$ and
$T_{k,n}=\mathbb{E}_{B_{n}}[T_{k,n}(B_{n})]$.

Given $B_{n}$ and $P_{n,B_{n}}^{0}$, let $Z_{k,i},\;1\leq i\leq np_{n}$,
denote the $np_{n}$ i.i.d. copies of $Z_{k}$ corresponding with the validation
set, i.e., with $\\{X_{i}:B_{n}(i)=1\\}$ (as defined in Lemma A1). Then,
$\hat{H}_{k}=\sum_{i}Z_{k,i}/np_{n}$ and
$\tilde{H}_{k}=\mathbb{E}_{P_{0}}[Z_{k,i}\mid P^{0}_{n,B_{n}},B_{n}]$. Hence,
$\tilde{H}_{k}-\hat{H}_{k}$ is an empirical mean of $np_{n}$ i.i.d. centered
random variables. Further, by Assumptions 1 and 2, $\lvert
Z_{k,i}\rvert<2(M_{1}+M_{2})^{2}J^{2}$ a.s.. Next, we apply Bernstein’s
Inequality to the centered empirical mean $\tilde{H}_{k}-\hat{H}_{k}$, using
the property of the $Z_{k,i}$’s derived in Lemma A1, to obtain a bound for the
tail probabilities of $R_{k,n}(B_{n})$ and $T_{k,n}(B_{n})$:

$\sigma^{2}_{k}\coloneqq\text{Var}_{P_{0}}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]\leq\overline{M}(J)\mathbb{E}\left[Z_{k}\mid
P^{0}_{n,B_{n}},B_{n}\right]=\overline{M}(J)\tilde{H}_{k}.$

Then, for $s>0$, Bernstein’s Inequality yields

$\begin{split}\mathbb{P}_{P_{0}}\left(R_{k,n}(B_{n})>s\mid
P_{n,B_{n}}^{0},B_{n}\right)&=\mathbb{P}_{P_{0}}\left(\tilde{H}_{k}-\hat{H}_{k}>\frac{s+\delta\tilde{H}_{k}}{1+\delta}\Bigg{|}P_{n,B_{n}}^{0},B_{n}\right)\\\
&\leq\mathbb{P}_{P_{0}}\left(\tilde{H}_{k}-\hat{H}_{k}>\frac{s+\delta\sigma_{k}^{2}/\overline{M}(J)}{1+\delta}\Bigg{|}P_{n,B_{n}}^{0},B_{n}\right)\\\
&\leq\text{exp}\left\\{-\frac{np_{n}}{2(1+\delta)^{2}}\frac{\left(s+\delta\sigma_{k}^{2}/\overline{M}(J)\right)^{2}}{\sigma_{k}^{2}+\frac{\overline{M}(J)}{3(1+\delta)}\left(s+\delta\sigma_{k}^{2}/\overline{M}(J)\right)}\right\\}.\end{split}$

Note that

$\frac{\left(s+\delta\sigma_{k}^{2}/\overline{M}(J)\right)^{2}}{\sigma_{k}^{2}+\frac{\overline{M}(J)}{3(1+\delta)}\left(s+\delta\sigma_{k}^{2}/\overline{M}(J)\right)}=\frac{s+\delta\sigma_{k}^{2}/\overline{M}(J)}{\frac{\sigma_{k}^{2}}{s+\delta\sigma_{k}^{2}/\overline{M}(J)}+\frac{\overline{M}(J)}{3(1+\delta)}}\geq\frac{s+\delta\sigma_{k}^{2}/\overline{M}(J)}{\overline{M}(J)\left(\frac{1}{\delta}+\frac{1}{3}\right)}\geq\frac{s}{\overline{M}(J)\left(\frac{1}{\delta}+\frac{1}{3}\right)}.$

And so, for $s>0$,

$\mathbb{P}_{P_{0}}\left(R_{k,n}(B_{n})>s\mid
P_{n,B_{n}}^{0},B_{n}\right)\leq\text{exp}\left\\{-\frac{np_{n}}{c\left(\delta,\overline{M}(J)\right)}s\right\\}\leq
K\;\text{exp}\left\\{-\frac{np_{n}}{c\left(\delta,\overline{M}(J)\right)}s\right\\},$

where
$c(\delta,\overline{M}(J))=2(1+\delta)^{2}\overline{M}(J)(\frac{1}{\delta}+\frac{1}{3})$.
The same bound applies for the marginal probabilities of
$\mathbb{P}_{P_{0}}(R_{k,n}(B_{n})>s)$ since they hold for arbitrary choices
of $B_{n}$ and $P_{n,B_{n}}^{0}$. The second inequality follows from $K$ being
larger than or equal to $1$ by definition.

Finally, for any $u>0$, we have by the properties of expectations and the
previously derived result that

$\begin{split}\mathbb{E}_{P_{0}}[R_{\hat{k},n}]&=\int_{0}^{\infty}\mathbb{P}_{P_{0}}(R_{\hat{k},n}>s)ds-\int_{-\infty}^{0}\mathbb{P}_{P_{0}}(R_{\hat{k},n}\leq
s)ds\\\ &\leq\int_{0}^{\infty}\mathbb{P}_{P_{0}}(R_{\hat{k},n}>s)ds\\\ &\leq
u+\int_{u}^{\infty}\mathbb{P}_{P_{0}}(R_{\hat{k},n}>s)ds\\\ &\leq
u+\int_{u}^{\infty}K\text{exp}\left\\{-\frac{np_{n}}{c\left(\delta,\overline{M}(J)\right)}s\right\\}ds.\end{split}$
(S11)

Since the expression on the right-hand side of the inequality above achieves
its minimum value of $c(\delta,\overline{M}(J))(1+\text{log}(K)/(np_{n})$ at
$u_{n}=c(\delta,\overline{M}(J))\text{log}(K)/(np_{n})$, then

$\mathbb{E}_{P_{0}}[R_{\hat{k},n}]\leq
c\left(\delta,\overline{M}(J)\right)\frac{1+\text{log}(K)}{np_{n}}.$

The same bound applies to $\mathbb{E}_{P_{0}}[T_{\tilde{k},n}]$. Therefore,
taking the expected values of the inequality in Equation (S10), we produce the
desired finite-sample result in Equation (11). ∎

###### Proof.

Theorem 1, High-Dimensional Asymptotic Result.

The expected risk differences ratio’s convergence follows directly from
Equation (11) for some $\delta>0$, so long as
$c(\delta,\overline{M}(J))(1+\text{log}(K))/(np_{n}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}])\rightarrow
0$ as $n,J\rightarrow\infty$. Given the assumption in Kolmogorov asymptotics
that $J/n\rightarrow m<\infty$ as $J,n\rightarrow\infty$, an equivalent
condition is that
$m(M_{1}+M_{2})^{2}J(1+\text{log}(K))/(p_{n}\mathbb{E}_{P_{0}}[\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}])\rightarrow
0$ as $n,J\rightarrow\infty$. Convergence in probability then follows from
Lemma A2.

∎

Though there are minor adaptations to the assumptions to reflect the use of
high-dimensional asymptotics, the proof of Corollary 1 follows that of
Corollary 1 in laan-dudoit:2003.

###### Proof.

Corollary 1.

The asymptotic statement of Equation (16) is an immediate result of Theorem
1’s Equation (13).

$\frac{\tilde{\theta}_{p_{n},n}(\hat{k}_{p_{n},n})-\theta_{0}}{\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}}\frac{\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0}}{\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0}}\overset{P}{\rightarrow}1.$

Letting
$Z_{1,n}\coloneqq(n(1-p_{n}))^{\gamma}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0})$
and $Z_{2,n}\coloneqq
n^{\gamma}(\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0})$, and assuming that
the sufficient condition in Equation (17) holds, we find that
$Z_{1,n}/Z_{2,n}\overset{d}{\rightarrow}1$ by the Continuous Mapping Theorem.
Then, notice that

$\frac{Z_{1,n}}{Z_{2,n}}=\frac{(1-p_{n})^{\gamma}(\tilde{\theta}_{p_{n},n}(\tilde{k}_{p_{n},n})-\theta_{0})}{\tilde{\theta}_{n}(\tilde{k}_{n})-\theta_{0}},$

which yields the desired result of Equation (15) when $p_{n}\rightarrow 0$. In
the case of single-split validation,
$Z_{2,n}\overset{d}{=}Z_{1,\frac{n}{1-p_{n}}}$, and so
$Z_{2,n}\overset{d}{\rightarrow}Z$ implies that
$(Z_{1,n},Z_{2,n})\overset{d}{\rightarrow}(Z,Z)$.

∎

##  S3 Supplementary Figures

Figure S1:  Comparison of the cross-validated ($\hat{k}_{p_{n},n}$) and cross-
validated oracle ($\tilde{k}_{p_{n},n}$) selections’ cross-validated
conditional risk differences. The proposed cross-validated selection procedure
achieves asymptotic equivalence in most settings for relatively small sample
sizes and numbers of features. Figure S2: Tasic dataset: Diagnostic plots and
tables generated using the cvCovEst R package. (A) The top-left plot presents
the cross-validated Frobenius risk of the estimator selected by our method,
which in this case is the POET estimator with $k=5$ and lambda = 0.3, and
alongside the other hyperparameter combinations considered. Here, $k$
represents the number of potential latent factors, and lambda the thresholding
value used. The top-right panel contains a line plot of the selected
estimator’s eigenvalues. The bottom-left plot displays the absolute values of
the estimated correlation matrix output by the cvCovEst selection, and the
bottom-right table lists the best performing estimators from all classes of
estimators considered. The table presented in (A) makes it clear that our
method’s selection substantially improves upon the sample covariance matrix in
terms of the cross-validated Frobenius risk. The comparison in (B) suggests
that the sample covariance matrix over-estimates the 20 leading eigenvalues,
resulting in a poorer UMAP representation. Figure S3: Zeisel dataset:
Diagnostic plots and tables generated using the cvCovEst R package. The
components in (A) can be interpreted in the same manner as the previous
figure, with the exception of the top-left panel. In this table, a five number
summary of the cross-validated Frobenius risk is given for each class of
estimators considered across all possible combinations of hyperparameters, if
any. It is clear from the tables in (A) that the cvCovEst selection is
essentially equivalent in terms of cross-validated risk to the sample
covariance matrix. The plots in (B) further highlight that the 20 leading
eigenvalue of the cvCovEst estimate and the sample covariance matrix are
indistinguishable.

##  S4 Supplementary Tables

Table S1: Families of candidate estimators compared to the cross-validated loss-based estimator selection procedure. Note that the library of candidate estimators used by the proposed method is provided in Table 1. Estimator | Hyperparameters
---|---
Sample covariance matrix | NA
Hard thresholding | Thresholds $=\\{0.05,0.10,\ldots,1.00\\}$
SCAD thresholding | Thresholds $=\\{0.05,0.10,\ldots,1.00\\}$
Adaptive LASSO | Thresholds $=\\{0.1,0.2,\ldots,0.5\\}$; exponential weights $=\\{0.1,0.2,\ldots,0.5\\}$
Banding | Bands $=\\{1,2,\ldots,10\\}$
Tapering | Bands $=\\{2,4,\ldots,10\\}$
Linear shrinkage | NA
Dense linear shrinkage | NA
Nonlinear shrinkage | NA
POET using hard thresholding | Latent factors $=\\{1,2,\ldots,10\\}$; thresholds $=\\{0.1,0.2,\ldots,1.0\\}$
Table S2: Families of candidate estimators used in single-cell transcriptomic data analyses. Estimator | Hyperparameters
---|---
Sample covariance matrix | NA
Hard thresholding | Thresholds $=\\{0.05,0.10,\ldots,0.30\\}$
SCAD thresholding | Thresholds $=\\{0.05,0.10,\ldots,0.50\\}$
Adaptive LASSO | Thresholds $=\\{0.1,0.2,\ldots,0.5\\}$; exponential weights $=\\{0.1,0.2,\ldots,0.5\\}$
Linear shrinkage | NA
Dense linear shrinkage | NA
POET using hard thresholding | Latent factors$=\\{5,6,\ldots,10\\}$; thresholds $=\\{0.05,0.10,\ldots,0.3\\}$

##  S5 Simulation Results: Spectral Norm

We performed an additional simulation analysis analogous to that presented in
Figure 3, but using the spectral norm instead of the Frobenius norm to assess
the accuracy of our estimation procedure with respect to the leading
eigenvalue of the covariance matrix. Recall that the spectral norm of a square
matrix is defined as it’s largest absolute eigenvalue. We applied the
estimation procedures listed in Table S1 to all simulated datasets along with
our data-adaptive selection procedure (still under the Frobenius loss) using
the candidates given in Table 1. When necessary, our cross-validation
framework was used to select hyperparameters within classes of estimators
using the Frobenius loss. We then computed the mean spectral norm against the
true covariance matrix stratified by $n$, $J/n$, and the covariance matrix
model.

Figure S4: Comparison of competing, bespoke covariance matrix estimation
procedures to our cross-validated selection approach in terms of the Monte
Carlo mean spectral norm under a variety of data-generating processes. Note
that the scales of the y-axis are tailored to the covariance matrix model.

The results, presented in Figure S4, demonstrate that our estimator selection
procedure performs at least as well as the best alternative estimation
strategy. This suggests that procedures dedicated to or relying upon the
accurate estimation of leading eigenvalues and eigenvectors, like principal
component analysis and latent variable estimation, might benefit from the
integration of our cross-validated covariance matrix estimation framework.

abbrvnat references