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

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 (2005)’s estimator selection framework is more amenable to loss functions that operate over random vectors. Accordingly, we propose the observation-level Frobenius loss:

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 $\bm{\psi}\in\mathbf{\Psi}$, and $\eta_{0}$ is a $J\times J$ matrix acting as a scaling factor, that is, 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

Ledoit and Wolf (2004), Bickel and Levina (2008a), 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; however, when the scaling factor is constant, it only impacts the interpretation of the loss. Constant scaling factors have no impact on our asymptotic analysis. Since it need not be estimated, it is not a nuisance parameter in the conventional sense.

When the scaling factor of Equation (4) is constant, the risk minimizers are identical for the cross-validated observation-level Frobenius risk and the common cross-validated Frobenius risk (used by, for example, Bickel and Levina, 2008a; Rothman et al., 2009; Fan et al., 2013; Fang et al., 2016).

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-diagonal elements’ errors by their corresponding diagonal entries, especially useful when the random variables are of different scales. Such a scaling factor might offer a more equitable evaluation across all entries of the parameter:

Here, $\eta_{0}=\text{diag}(\bm{\psi}_{0})$ is a genuine nuisance parameter which can be estimated via the diagonal entries of the sample covariance matrix.

Finally, the covariance matrix $\bm{\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 $\bm{\psi}$ are interchangeable.

The proof is provided in the Online Supplement. Our proposed loss therefore satisfies the condition of Equation (1). 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}=\bm{1}$, and suppress dependence of the loss function on the scaling factor wherever possible throughout the remainder of the text.

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 (2)). This theorem extends the general framework of Dudoit and van der Laan (2005) 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 (3). All proofs are provided in Section LABEL:supp:proofs of the Online Supplement.

The proof relies on special properties of the random variable $Z_{k}~{}\coloneqq~{}L(X;\hat{\Psi}_{k}(P_{n}))~{}-~{}L(X;\bm{\psi}_{0})$ and on an application of Bernstein’s inequality (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 and the set of possible estimates 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 is technical in nature — it makes the proofs tractable. While it may appear stringent, and, for instance, is not satisfied by Gaussian distributions, we believe it to be generally applicable. We stress that parametric data-generating distributions, like those exhibiting Gaussianity, rarely reflect reality, that is, they are merely mathematical conveniences¹¹1Anecdotally, one cannot help but find themself reminded that “Everyone is sure of this [that errors are normally distributed] $\ldots$ since the experimentalists believe that it is a mathematical theorem, and the mathematicians that it is an experimentally determined fact.” (Poincaré, 1912, p. 171). Most random variables, or transformations thereof, arising in scientific practice 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, constraining them to $[0,2^{16})$. Similarly, the measurement of immunologic markers, of substantial interest in vaccine efficacy trials of HIV-1, COVID-19, and other infectious diseases, are bounded by the limits of detection and/or quantitation imposed by assay biotechnology.

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 (7) and Equation (8) 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. This intuition is corroborated by our extensive simulation study in the following section. Consistent estimation in terms of the Frobenius risk is therefore not possible in high-dimensions without additional assumptions about the data-generating process.

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 (2008a) and accompanying hard-thresholding estimators (see Section LABEL:supp:thresh-est of the Online Supplement):

Proposition 3 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 (2008a)’s Theorem 2 and its subsequent discussion. Plainly, the true covariance matrix cannot be too sparse if the collection of considered estimators contains the hard-thresholding estimator with a correctly specified thresholding hyperparameter value.

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. 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. In situations where the true covariance matrix is known to possess this level sparsity, practitioners might instead appeal to Equation (39) of Bickel and Levina (2008a) to support their use of a cross-validated estimator selection procedure. This result, coupled with that of Proposition 1, likely explains the aforementioned simulation findings of the banded and tapered covariance matrices.

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 (3), however, we would like our cross-validated procedure to be asymptotically equivalent to the oracle over the entire data set. The conditions to obtain this desired result are provided in Corollary 1, a minor adaptation of previous work by Dudoit and van der Laan (2005). This extension accounts for increasing $J$, thereby permitting its use in high-dimensional asymptotics.

The proof 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 (9) is that there exists a $\gamma>0$ such that

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 (9) 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.

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 (7) and (8) are achieved, respectively.

Finally, we have assumed thus far that $\mathbb{E}_{P_{0}}[X]$ is known. This is generally not the case in practice. In place of a random vector centered at zero, we might instead consider the set of $n$ demeaned random vectors $\tilde{\mathbf{X}}_{n\times J}$ where $\tilde{X}_{i}=X_{i}-\bar{X}$ and $\bar{X}^{(j)}=1/n\sum X_{i}^{(j)}$. It follows from the details given in Remark LABEL:supp:remark:unknown-mean of the Online Supplement that the asymptotic results of Theorem 1 and Corollary 1 apply to $\tilde{\mathbf{X}}_{n\times J}$.

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 the Online Supplement’s Equation (LABEL:supp:factor-model). Let $\bm{\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 $\bm{\psi}=\bm{\beta}\bm{\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 data sets consisting of $n\in\{50,100,200,500\}$ i.i.d. multivariate Gaussian, mean-zero observations. The uniform boundedness condition of Theorem 1’s Assumption 1 is therefore not satisfied; we do this purposefully to further stress that this assumption is not limiting in many practical settings. 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 the Online Supplement’s Table LABEL:supp:table:sim-hyperparams, which includes details on these estimators’ possible hyperparameters. Seventy-four estimators make up the library of candidates. 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.

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}$. The Monte Carlo expectations of the risk differences stratified by $n$, $J/n$, and the covariance model were computed from the cross-validated conditional risk differences. The ratios of the expected risk differences are presented in Figure 1A. 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.

A stronger result, corresponding to Equation (8) of Theorem 1, is presented in the Online Supplement’s Figure LABEL:supp:sim:conv-in-prob. 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 1B 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 and spectral norms 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 (Figures LABEL:supp:mean-frobenius and LABEL:supp:mean-spectral of the Online Supplement). Our data-adaptive procedure is found to perform at least as well as the best alternative estimation strategy across all simulation scenarios under both norms.

We thank Brian Collica and Jamarcus Liu for significant contributions to the development of the cvCovEst R package. PB is funded by the Fonds de recherche du Québec - Nature et technologies, and the Natural Sciences and Engineering Research Council of Canada. NSH gratefully acknowledges the support of the National Science Foundation (award no. DMS 2102840).