2.1 Diagonal estimation

Method 1: Diagonal Estimator (DE)

Under the “large p small n” setting, one naive approach is to estimate Σ by the diagonal sample covariance matrix, i.e., diag(Sn). This estimator was first considered in [37] to propose a diagonal linear discriminant analysis. It was further considered in [38] where the authors demonstrated that a diagonal covariance matrix estimation can be sometimes reasonable when p is much larger than n. Let diag(Σ)=diag(σ12,…,σp2) where σj2 are the covariate-specific variances for j=1,…,p, and diag(Sn)=diag(s12,…,sp2) where sj2 are the sample variances of σj2, respectively. By letting Σˆ=diag(Sn), we define the first estimator of θ as

We refer to θˆ(1) as the diagonal estimator (DE). To be specific, DE is proposed to estimate log|diag(Σ)| rather than log|Σ|.

Method 2: Improved Diagonal Estimator (IDE)

It is noteworthy that DE may not perform well as an estimate of log|diag(Σ)| when the sample size is small, mainly due to the unreliable estimates of the sample variances. Various approaches have been proposed to improving the variance estimation in the literature. See, for example, [39, 40, 4142], and [43].

To improve DE, we consider the optimal shrinkage estimator in [42],

where spool2=∏i=1p(sj2)1/p, hp(1)=(ν/2)Γ(ν/2)/Γ(ν/2+1/p)p with ν=n−1, Γ(⋅) is the Gamma function, and α∈[0,1] is the shrinkage parameter. Replacing sj2 in DE by σˆj2, we have

where C=log{hpαp(1)h1(1−α)p(1)} is a constant.

The estimation structure in eq. (2) shows that the DE estimator, θˆ(1), can be further improved. Specifically, if we have C0 such that E(θˆ(1)+C0)=log|diag(Σ)|, then C0 is defined as the optimal C value so that the estimator θˆ(1)+C0 minimizes the mean squared error in the family of estimators {θˆ(1)+C: C∈(−∞,∞)}.

The proof of Theorem 1 is given in the Appendix. By eq. (2) and Theorem 1, we define the second estimator of θ as

We refer to θˆ(2) as the improved diagonal estimator (IDE).

2.2 Shrinkage estimation

Recall that the sample covariance matrix Sn is singular when the dimension is larger than the sample size. To overcome the singularity problem, other than the diagonal methods in Section 2.1, one may also estimate the covariance matrix by the following convex combination:

where T is the target matrix, and δ∈[0,1] is the shrinkage parameter. Both the target matrix and the shrinkage parameter play an important role in the shrinkage estimation. For instance, if we let T=diag(Sn) and δ=1, then S∗ reduces to the DE estimator.

The appropriate choice of the target matrix has been extensively studied in the literature. See, for example, [6, 33, 44, 45], and [7] and the references therein. Note that T is often chosen to be positive definite and well-conditioned, and consequently, the final estimate S∗ is also guaranteed positive definite and well-conditioned for any dimensionality. As suggested in [33] and [7], we consider a popular target matrix for nonhomogeneous variances: the “diagonal, unequal variance” matrix, i.e., the diagonal sample covariance matrix diag(Sn).

We also note that, given the target matrix, the estimation of the shrinkage parameter δ is also crucial to the final estimation. The available estimation methods for the shrinkage parameter are mainly: (1) the unbiased estimation, and (2) the consistent estimation. The unbiased estimation is replacing unknown terms in the optimal value by their unbiased estimators [33]. Whereas, the consistent estimation is replacing the unknown terms in the optimal shrinkage parameter with (n,p)-consistent estimators [7]. Taken together, we present below the four methods for estimating the covariance matrix and consequently for estimating θ, respectively.

Method 3: Unbiased Shrinkage Estimator with T=I (USIE)

Letting the target matrix be T=I, [33] proposed an unbiased estimator for the shrinkage parameter, denoted by δˆ1∗. This leads to S∗=δˆ1∗I+(1−δˆ1∗)Sn. We then define the third estimator of θ as

Method 4: Consistent Shrinkage Estimator with T=I (CSIE)

Letting the target matrix be T=I, [7] proposed a consistent estimator for the shrinkage parameter, denoted by δˆ2∗. This leads to S∗=δˆ2∗I+(1−δˆ2∗)Sn. We then define the fourth estimator of θ as

Method 5: Unbiased Shrinkage Estimator with T=diag(Sn) (USDE)

Letting T=diag(Sn), [33] also proposed an unbiased estimator for the shrinkage parameter, denoted by δˆ3∗. This leads to S∗=δˆ3∗diag(Sn)+(1−δˆ3∗)Sn. We then define the fifth estimator of θ as

Method 6: Consistent Shrinkage Estimator with T=diag(Sn) (CSDE)

Letting T=diag(Sn), [7] also proposed a consistent estimator for the shrinkage parameter, denoted by δˆ4∗. This leads to S∗=δˆ4∗diag(Sn)+(1−δˆ2∗)Sn. We then define the sixth estimator of θ as

2.3 Sparse estimation

When p is much larger than n, the shrinkage methods in Section 2.2 may not achieve a significant improvement over Sn. In such settings, to have a good estimate of Σ, one may have to impose some structural assumptions such as sparsity in the parameters. Recently, [15] reviewed some methods on estimating structured high-dimensional covariance and precision matrix. A typical sparsity is to assume that most of the off-diagonal elements in the covariance matrix are zero. To estimate the covariance matrix under a sparsity condition, various thresholding-based methods have been proposed in the literature that aim to locate some “large” off-diagonal elements. See, for example, [8, 9, 46, 47, 48, 49, 5051], and [52]. Particularly, the adaptive thresholding estimator proposed by [49] achieves the optimal rate of convergence over a large class of sparse covariance matrix under a wide spectral norms. Besides, it can be shown that the adaptive thresholding estimator also attains the optimal convergence rate under Bregman divergence losses over a large parameter class [15, 50]. Therefore, we also consider the sparsity methods as a representative and use them to estimate θ, i.e., the log-determinant of the covariance matrix.

Method 7: Adaptive Thresholding Estimator (ATE)

Bickel and Levina [8] proposed a universal thresholding method where all entries in the sample covariance matrix are thresholded by a common value γ. They required that the variances σj2 are uniformly bounded by a constant K, and consequently, the variances of the entries of the sample covariance matrix are also uniformly bounded. However, it was shown that a universal thresholding method is suboptimal over a certain class of sparse covariance matrices.

To improve the method above, [49] proposed an adaptive thresholding estimator for the covariance matrix:

where γij is the corresponding threshold of σ˜ij∗, and sγij(⋅) is a generalized thresholding operator [47], which is specified as the soft thresholding throughout simulations. With the proper γij, the estimator Σˆ∗ adaptively achieves the optimal rate of convergence over a large class of sparse covariance matrix under the spectral norm. Now by Σˆ∗, the seventh estimator of θ is

We refer to θˆ(7) as the adaptive thresholding estimator (ATE).

2.4 Factor model estimation

The sparsity condition on the covariance matrix assumes that most of covariates are uncorrelated to each other. Note that, however, this assumption may not be realistic in practice. Recently, under the assumption of conditional sparsity, [54] introduced a principle orthogonal complement thresholding method using the factor model. In this section, we briefly review their method and then apply it to estimate the log-determinant of the covariance matrix.

Method 8: Principal Orthogonal Complement Thresholding Estimator (POET)

Fan et al. [54] considered the approximate factor model:

where yg=(y1g,…,ypg)T is the observed response, B=(b1,…,bp)T is the loading matrix, fg is a Q×1 vector of common factors, and ug=(u1g,…,upg)T is the error vector. In this model, we can only observe yg. Let

where Σu is the covariance matrix of ug. To estimate Σ, [54] applied the spectral decomposition on the sample covariance matrix:

where λˆ1≥λˆ1≥…≥λˆp are eigenvalues of Sn, ξˆj, j=1,…,p are the corresponding eigenvectors, and RˆQ=∑j=Q+1pλˆjξˆjξˆjT is the principal orthogonal complement. For this decomposition, the first Q principal components were kept and the thresholding was applied on RˆQ. Here, the generalized thresholding operator can be used. In addition, [54] also introduced a method to obtain an estimation of Q, denoted by Qˆ. Their final estimator of Σ is

where RˆQˆT is the thresholding result of RˆQ. Now by eq. (8), we define the last estimator of θ as

We refer to θˆ(8) as the principal orthogonal complement thresholding estimator (POET).

3.1 Normal data

In this setup, we consider a block diagonal structure for the covariance matrix. This structure is widely adopted in the literature, e.g., [55] and [56]. Specifically, we let

where D=diag(σ12,…,σp2) with σj2 being i.i.d. from the distribution χ52/5, and R follows a block diagonal structure:

In our simulations, we consider Σρ=(σij(ρ))q×q with σij(ρ)=ρ|i−j| for 1≤i,j≤q. In addition, we set ρ=0, 0.3, 0.6 or 0.9, to represent different levels of dependence, and (p,q)=(50,5) or (300,10), respectively.

Figures 1 and 2 display the log(MSE) of the eight methods for different levels of dependence, dimension and sample size. From these figures, we have the following findings. When the covariates are uncorrelated, IDE gives the best performance under a high dimension (e.g., p=300). However, if the dimension is not large (e.g., p=50), and the covariates are uncorrelated or weakly correlated, shrinking the covariance matrix toward an identity matrix leads to a better performance under a small sample size. This is because when the sample size is small, the variances of the entries of the sample covariance matrix are large. Hence, CSIE and USIE stabilize both diagonal and off-diagonal entries and, at the same time, an identity target possesses an explicit structure which in turn requires little data to fit. Consequently, the resulting estimators have a good bias–variance tradeoff. In addition, when the correlation and dimension are both large, imposing additional structure assumptions is necessary. Under this situation, ATE and POET turn out to be the best two methods among the eight methods unless the sample size is relatively small. When the sample size is small, the pattern of ATE is very similar to that of DE. When the sample size and dimension are both large, ATE outperforms all other methods except for POET.

Figure 3 displays the performance of the eight methods for different levels of dependence with p=300. The pattern is consistent with Figure 2. In particular, as the correlation and sample size are large, the performance of POET is satisfactory. From Figures 1 and 2, however, we note that the log(MSE) of POET tends to be oscillating as the sample size increases. This may due to that POET depends on the estimated number of factors K. In [54], the authors used a consistent estimator for K and showed that POET is robust to over-estimated number of factors under the spectral norm. Our simulations in Table 1, however, show that the robustness for estimating the covariance matrix may not hold any more when the purpose is to estimate the determinant. In particular for small sample sizes, either over-estimated or under-estimated K leads to a large bias for the determinant estimator.

Figure 4

Log MSEs for data from mixture normal distribution with p=50, and ρ ranging from 0 to 0.9. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.

Figure 5

Log MSEs for data from mixture normal distribution with p=300, and ρ ranging from 0 to 0.9. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.

3.2 Mixture normal data

In this setup, we consider a mixture model where the random vectors are generated from

where f1(X) and f2(X) are the density functions of Np(μ3,Σ3) and Np(μ4,Σ4), respectively. For the covariance matrices, we consider a sparse block diagonal structure as follows:

where D=diag(σ12,…,σp2) with σj2 being i.i.d. from the distribution (1/5)χ52, and R(ρ) being the same as in Setup II. For simplicity, we set α1=α2=1/2 and μ1=μ2=0. Under this setting, the covariance matrix of X is simplified as (Σ3+Σ4)/2, which results in a highly sparse matrix where the odd off-diagonal parts in diagonal blocks are zeros. We set (p,q)=(50,5) or (300,10), and ρ=0, 0.3, 0.6 or 0.9.

Figure 6

Log MSEs for data from heavy-tailed distribution with p=50, and ρ ranging from 0 to 0.9. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.

Figure 7

Log MSEs for data from heavy-tailed distribution with p=300, and ρ ranging from 0 to 0.9. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.

Figures 4 and 5 display the log(MSE) of the eight methods under different levels of dependence and sample size. When the sample size is large and the covariates are uncorrelated, IDE gives the best performance. When the sample size is small and the dimension is not large (e.g., n=5,p=50), shrinking the covariance matrix toward an identity matrix (e.g., USIE and CSIE) outperforms the other methods except that the correlation is very large (e.g., ρ=0.9). However, as the sample size and dimension are both large, the shrinkage methods will become suboptimal. Instead,if the correlation is also large (e.g., ρ=0.6), ATE and POET outperform the other methods in most settings. As aforementioned, the performance of POET is not stable and may not be satisfactory when the sample size is not large.

3.3 Heavy-tailed data

In this setup, we consider to simulate heavy-tailed data from a log-normal distribution, lnN(μ,σ2), where the mean and variance are eμ+σ2/2 and (eσ2−1)e2μ+σ2, respectively. First of all, we generate n independent random vectors Zi=(zi1,…,zip)T, where all the components of Zi are sampled independently from lnN(0,1). Let Xi=Σ1/2Zi∗ with Zi∗=(zi1−e1/2,…,zip−e1/2)T/{e(e−1)}1/2, and Σ is a p×p positive definite matrix. Consequently, the mean vector and covariance matrix of Xi are 0p×1 and Σp×p, respectively. For the covariance matrix, we consider the block diagonal structure as described in Section 3.1. We set (p,q)=(50,5) or (300,10), and ρ=0, 0.3, 0.6 or 0.9.

Figures 6 and 7 display the log(MSE) of the eight methods under different levels of dependence and sample size. When the dimension and correlation are both small, USIE and CSIE outperform the other methods. The reason is similar as the discussion in Section 3.1, the heavy-tailed data may lead to unstable estimates of the entries of Σ, hence shrinking towards a simple identity target, which requires little data to fit, stabilizes the sample covariance matrix. In addition, as shown in Figure 7, when the dimension is large and the correlation is not small, ATE and POET are the only two methods that have a better performance than the other methods except that the sample size is small. Finally, we also note that IDE cannot provide a satisfactory performance even if the covariates are uncorrelated. As demonstrated in Theorem 1, IDE estimator is derived under the normal distribution and may not be robust to heavy-tailed data.

3.4 Degenerate normal data

To further investigate the performance of the eight methods, we consider a non-positive definite covariance matrix in which the positive definite assumption of the covariance matrix is violated. Note that this new setting is often overlooked in the literature. To construct a non-positive definite covariance matrix, we define the affine transformation C as

We then apply the affine transformation to the covariance matrix in Setup II and form

Figure 8

Log MSEs for data from degenerate normal distribution with p=50. The sample size ranges from 5 to 50. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.

It is obvious that |Σ5|=0 since |C|=0. We set (p,q)=(50,5), and ρ=0, 0.3, 0.6 or 0.9. Note that he log-determinant of Σ5 is negative infinity. Hence, for this degenerate setting, the MSE is defined on the determinant rather than on the log-determinant. Specifically, it is

Figure 8 shows the log(MSE) of all eight methods for different levels of dependence and sample size. We can see that the simulation results are different from those in the previous three setups. POET gives the best performance among the eight methods. In addition, we note that, under the non-positive definite setting, POET performs extremely well when the sample size is very small. For this phenomenon, we explore the possible reasons in the next paragraph.

To estimate Σ, [54] applied the spectral decomposition on the sample covariance matrix:

If the sample size is much smaller than the dimension p, most eigenvalues of Sn are zeros. This leads to RˆQ, the principal orthogonal complement of the largest Q eigenvalues, is nearly a zero matrix. And consequently, the final estimator of POET, ΣˆQˆ=∑j=1QˆλˆjξˆjξˆjT+RˆQˆT, tends to be highly degenerate for small sample sizes rather than for large sample sizes.

Finally, it is noteworthy that when the correlation is strong, the log(MSE) of POET is also fluctuant as the sample size increases. This again verifies that both the correlation and sample size have a large impact on the performance of POET.

3.5 Real data

In this setup, we consider to generate a realistic covariance matrix from the Myeloma data [57], which is a real microarray data set including a total of 54, 675 genes, with 351 samples in the first group and 208 samples in the second group. To generate the covariance matrix, we first select 100 genes randomly from the first group and then compute the sample covariance matrix using the selected genes, denoted by Σr. Next, to evaluate the performance of the estimators under different levels of dependence, we follow [58] and define the true covariance matrix as

where ρ controls the level of dependence. We set ρ=0, 1/3, 2/3 or 1. Note that ρ=0 corresponds to a diagonal covariance matrix, and ρ=1 treats the generated sample covariance matrix as the true covariance matrix.

Figure 9 shows the log(MSE) of the eight methods for different levels of dependence and sample size. The comparison results are summarized as follows. When the sample size and correlation are both small, the methods that shrinking the covariance matrix toward the identity matrix (e.g., USIE and CSIE) lead to a good performance. When the covariates are uncorrelated and the sample size is large, IDE has the best performance. In addition, when the sample size is large and the correlation is moderate (e.g., n=80 and ρ=2/3), shrinking the sample covariance matrix toward a diagonal target matrix (e.g., USDE and CSDE) has a good performance. When the correlation and sample size are both large, ATE outperforms or is at least comparable to USDE and CSDE. Finally, POET is not stable and very sensitive to both the correlation and the sample size. When the correlation and sample size is not large, POET may fail to provide a satisfactory performance owing to the largely increased bias compared with the other methods.

Figure 9

Log MSEs for real data with p=100. The sample size ranges from 10 to 80. In all figures, “1” to “8” represent the eight methods: DE [38], IDE, USIE [33], CSIE [7], USDE [33], CSDE [7], ATE [49], and POET [54], respectively.