A unified approach for covariance matrix estimation under Stein loss

Anis M. Haddouche¹¹footnotemark: 1 Wei Lu ²²footnotemark: 2

Abstract

In this paper, we address the problem of estimating a covariance matrix of a multivariate Gaussian distribution, relative to a Stein loss function, from a decision theoretic point of view. We investigate the case where the covariance matrix is invertible and the case when it is non–invertible in a unified approach.

keywords:

Orthogonally invariant estimators , singular covariance matrix , statistical decision theory , high–dimensional statistics.

MSC:

[2010]
62H12 , 62F10 , 62C99.

^†^†journal: Statistics and Probability Letters

\affiliation

[label1]organization=INSA Rouen, Normandie Univ, UNIROUEN, UNIHAVRE, INSA Rouen, LITIS and LMI, addressline=avenue de l’Université, BP 8, city=Saint-Étienne-du-Rouvray, postcode=76801, country=France., ead= Mohamed.haddouche@insa-rouen.fr \affiliation[label2]organization=INSA Rouen, Normandie Univ, UNIROUEN, UNIHAVRE, INSA Rouen, LMI, addressline=avenue de l’Université, BP 8, city=Saint-Étienne-du-Rouvray, postcode=76801, country=France., ead= wei.lu@insa-rouen.fr

1 Introduction

Let $\boldsymbol{X}$ be an observed $p\times n$ matrix of the form

\displaystyle\boldsymbol{X}=\boldsymbol{B}\boldsymbol{Z},

(1)

where $\boldsymbol{B}$ is a $p\times r$ matrix of unknown parameters, with $r\leq p$ , and $\boldsymbol{Z}$ is a $r\times n$ random matrix. Assume that $r$ is known and that the columns of $\boldsymbol{Z}$ are identically and independently distributed as the $r$ -dimensional multivariate normal distribution $\mathcal{N}_{r}({0}_{r},\boldsymbol{I}_{r})$ . Then the columns of $\boldsymbol{X}$ are identically and independently distributed from the $p$ -dimensional multivariate normal $\mathcal{N}_{p}({0}_{p},\boldsymbol{\Sigma})$ , where $\boldsymbol{\Sigma}=\boldsymbol{B}\boldsymbol{B}^{T}$ is the unknown $p\times p$ covariance matrix with

{\rm rank}(\boldsymbol{\Sigma})=r\leq p\,.

It follows that the $p\times p$ sample covariance matrix $\boldsymbol{S}=\boldsymbol{X}\boldsymbol{X}^{T}$ has a singular Wishart distribution (see Srivastava (2003)) such that

{\rm rank}(\boldsymbol{S})=\min(n,r)=q\leq p\,,

with probability one. We denote in the following by $\boldsymbol{S}^{+}$ and $\boldsymbol{\Sigma}^{+}$ the Moore-Penrose inverses of $\boldsymbol{S}$ and $\boldsymbol{\Sigma}$ respectively.

We consider the problem of estimating the covariance matrix $\boldsymbol{\Sigma}$ under the Stein type loss function

\displaystyle L(\boldsymbol{\hat{\Sigma}},\boldsymbol{\Sigma})={\rm tr}({\boldsymbol{\Sigma^{+}{\hat{\Sigma}}}})-\ln\!|\Lambda(\boldsymbol{\Sigma^{+}{\hat{\Sigma}}})|-q\,,

(2)

where $\boldsymbol{\hat{\Sigma}}$ estimates $\boldsymbol{\Sigma}$ and $\Lambda(\boldsymbol{\Sigma^{+}{\hat{\Sigma}}})$ is the diagonal matrix of the $q$ positives eigenvalues of $\boldsymbol{\Sigma^{+}{\hat{\Sigma}}}$ . The corresponding risk function is denoted by

\displaystyle R(\boldsymbol{\hat{\Sigma}},\boldsymbol{\Sigma})={E}[L(\boldsymbol{\hat{\Sigma}},\boldsymbol{\Sigma})]\,,

where ${E}(\cdot)$ denotes the expectation with respect to the model (1). Note that the loss function (2) is an adaptation of the original Stein loss function (see Stein (1986)) to the context of the model (1) (see Tsukuma (2016) for more details).

The difficulty of covariance estimation is commonly characterized by the ratio $p/n$ . The usual estimators of the form

\displaystyle\boldsymbol{\hat{\Sigma}}_{a}={a}\,\boldsymbol{S},\quad\text{with}\quad a>0,

(3)

perform poorly when $n,p\to\infty$ with $p/n\to c>0$ (see Ledoit and Wolf (2004)). Hence, in this situation, alternative estimators are needed. Indeed, as pointed out by James and Stein (1961), the larger (smaller) eigenvalues of $\boldsymbol{\Sigma}$ are overestimated (underestimated) by those estimators. Therefore, a possible approach to derive an improved estimators is to regularize the eigenvalues of $\boldsymbol{\hat{\Sigma}}_{a}$ . This fact suggest to consider the class of orthogonally invariant estimators (see Takemura (1984)) in (5) below.

Considering the model (1), we deal, in a unified approach, with the following cases.

1.

$n<r=p$ : $\boldsymbol{\Sigma}$ is invertible of ${\rm rank}$ $p$ and $\boldsymbol{S}$ is non–invertible of ${\rm rank}$ $n$ ;
2.

$r=p\leq n$ : $\boldsymbol{\Sigma}$ and $\boldsymbol{S}$ are invertible;
3.

$r<p\leq n$ : $\boldsymbol{\Sigma}$ and $\boldsymbol{S}$ are non–invertible of rank $r$ ;
4.

$r\leq n<p$ : $\boldsymbol{\Sigma}$ and $\boldsymbol{S}$ are non–invertible of rank $r$ ;
5.

$n<r<p$ : $\boldsymbol{\Sigma}$ and $\boldsymbol{S}$ are non–invertible of ranks $r$ and $n$ respectively.

The class of orthogonally invariant estimators was considered by various authors. For instance, see Stein (1986), Dey and Srinivasan (1985) and Haff (1980) for the case $(i)$ , Konno (2009) and Haddouche et al. (2021) for the cases $(i)$ and $(ii)$ . See also Chételat and Wells (2016) for the cases $(iii)$ and $(iv)$ . Recently Tsukuma (2016) extend the Stein (1986) estimator to the five possible cases above in a unified approach. Similarly, we extend here the class of Haff (1980) estimators to the context of the model (1).

The rest of this paper is organized as follows. In Section 2, we derive the improvement result of the proposed estimators over the usual estimators. We study numerically the behavior of the proposed estimators in Section 3.

2 Main result

Improving the class of the natural estimators in (3) relies on improving the optimal estimator among this class, that is, the one which minimizes the loss function (2).

Proposition 1 (Tsukuma (2016))

Under the Stein loss function (2), the optimal estimator among the class (3) is given by

\displaystyle\boldsymbol{\hat{\Sigma}}_{a_{o}}={a_{o}}\boldsymbol{S},\quad\text{where}\quad a_{o}=\frac{1}{m}\quad\text{and}\quad m=\max{(n,r)}.

(4)

As mentioned in Section 1, we consider the class of orthogonnally invariant estimators. Let $\boldsymbol{S}=\boldsymbol{H}\,\boldsymbol{L}\,\boldsymbol{H}^{\top}$ be the eigenvalue decomposition of $\boldsymbol{S}$ where $\boldsymbol{H}$ is a $p\times q$ semi–orthogonal matrix of eigenvectors and $\boldsymbol{L}={\rm diag}(l_{1},\dots,l_{q})$ , with $l_{1}>,\dots,>l_{q}$ , is the diagonal matrix of the $q$ positive corresponding eigenvalues (see Kubokawa and Srivastava (2008) for more details). The class of orthogonally invariant estimators is of the form

\displaystyle\hat{\boldsymbol{\Sigma}}_{\Psi}

\displaystyle=a_{o}\,\big{(}\boldsymbol{S}+\boldsymbol{H}\,\boldsymbol{L}\,\boldsymbol{\Psi}(\boldsymbol{L})\,\boldsymbol{H}^{\top})\,

(5)

with $\boldsymbol{\Psi}(\boldsymbol{L})={\rm diag}(\psi_{1}(\boldsymbol{L}),\dots,\psi_{q}(\boldsymbol{L}))$ , where $\psi_{i}(\boldsymbol{L})$ ( $i=1,\dots,q$ ) is a differentiable function of $\boldsymbol{L}$ .

More precisely, we consider an extension of the class of Haff (1980) estimators, to the context of the model (1), defined as

\displaystyle\hat{\boldsymbol{\Sigma}}_{\alpha}=a_{o}\big{(}\boldsymbol{S}+\boldsymbol{H}\boldsymbol{L}\boldsymbol{\Psi}(\boldsymbol{L})\boldsymbol{H}^{\top}\big{)}\,\text{with }\,\,\alpha\geq 1\,,\,b>0\,\,\text{and}\,\,\,\boldsymbol{\Psi}(\boldsymbol{L})=b\frac{\boldsymbol{L}^{-\alpha}}{{\rm tr}(\boldsymbol{L}^{-\alpha})},

(6)

where $a_{o}$ is given in (4). We give in the following proposition our main result.

Proposition 2

The Haff type estimators in (6) improves over the optimal estimator in (4), under the loss function (2), as soon as

\displaystyle 0<b\leq b_{o}=\frac{2\,(q-1)}{m-q+1}\,.

Proof 1

We aim to show that the risk difference between the Haff type estimators in (6) and the optimal estimator in (4), namely,

\displaystyle\Delta_{(\alpha,a_{o})}=R(\boldsymbol{\hat{\Sigma}}_{\alpha},\boldsymbol{\Sigma})-R(\boldsymbol{\hat{\Sigma}}_{a_{o}},\boldsymbol{\Sigma}),

(7)

is non–positive. Note that $\boldsymbol{\hat{\Sigma}}_{\alpha}$ can be written as

\displaystyle\hat{\boldsymbol{\Sigma}}_{\alpha}=a_{o}\,\boldsymbol{H}\boldsymbol{L}\boldsymbol{\Phi}(\boldsymbol{L})\boldsymbol{H}^{\top}\quad\text{with}\quad\boldsymbol{\Phi}(\boldsymbol{L})=\boldsymbol{I}_{q}+\boldsymbol{\Psi}(\boldsymbol{L}).

The risk of these estimators under the Stein loss function (2) is given by

\displaystyle R(\hat{\boldsymbol{\Sigma}}_{\alpha},\boldsymbol{\Sigma})={E}\big{(}{\rm tr}(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha})\big{)}-{E}\big{(}\!\ln\!|\Lambda(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha})|\big{)}-q\,.

(8)

First, dealing with ${E}\big{(}{\rm tr}(\hat{\boldsymbol{\Sigma}}_{\alpha}\boldsymbol{\Sigma}^{+})\big{)}$ , we apply Lemma A.2 in Tsukuma (2016) in order to get rid of the unknown parameter $\boldsymbol{\Sigma}^{+}$ . It follows that,

\displaystyle{E}\big{(}{\rm tr}(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha})\big{)}

\displaystyle=a_{o}{E}\left(\sum_{i=1}^{q}\!\left\{\!(m-q+1)\varphi_{i}+2l_{i}\frac{\partial\varphi_{i}}{\partial l_{i}}+2\sum_{j>i}^{q}\frac{l_{i}\varphi_{i}-l_{j}\varphi_{j}}{l_{i}-l_{j}}\!\right\}\!\right),

(9)

where, for $i=1,\dots,q$ ,

	$\displaystyle\phi_{i}=1+b\,\frac{l_{i}^{-\alpha}}{{\rm tr}(L^{-\alpha})},\quad\frac{\partial\varphi_{i}}{\partial l_{i}}=b\,\alpha\,\frac{1-{\rm tr}(\boldsymbol{L}^{-\alpha})\,l_{i}^{\alpha}}{{\rm tr}^{2}(\boldsymbol{L}^{-\alpha})\,l_{i}^{1+2\alpha}}$
and
	$\displaystyle\sum_{j>i}^{q}\frac{l_{i}\varphi_{i}-l_{j}\varphi_{j}}{l_{i}-l_{j}}=\sum_{j>i}^{q}\left(1+\frac{b}{{\rm tr}(\boldsymbol{L}^{-\alpha})}\left\{\frac{l_{i}^{1-\alpha}-l_{j}^{1-\alpha}}{l_{i}-l_{j}}\right\}\right)\,.$

Using the fact, for $j>i$ , $l_{j}>l_{i}$ , it can be shown that

\displaystyle\sum_{j>i}^{q}\frac{l_{i}\varphi_{i}-l_{j}\varphi_{j}}{l_{i}-l_{i}}\leq(q-i)\,.

(10)

Therefore, using (10), we obtain

	$\displaystyle{E}({\rm tr}(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha}))$	$\displaystyle\leq a_{o}\,{E}\Biggr{(}\sum_{i=1}^{q}\Biggr{\{}(m-q+1)\left(1+b\frac{{l_{i}}^{-\alpha}}{{\rm tr}(\boldsymbol{L}^{-\alpha})}\right)$
		$\displaystyle\hskip 76.82234pt+2\,b\,\alpha\,\frac{1-{\rm tr}(\boldsymbol{L}^{-\alpha})\,l_{i}^{\alpha}}{{\rm tr}^{2}(\boldsymbol{L}^{-\alpha})\,l_{i}^{2\alpha}}+2(q-i)\Biggl{\}}\Biggl{)}$
		$\displaystyle=a_{o}\,m\,q+a_{o}\,b\,{E}\Biggr{(}\sum_{i=1}^{q}\Biggr{\{}(m-q+1)\frac{{l_{i}}^{-\alpha}}{{\rm tr}(\boldsymbol{L}^{-\alpha})}$
		$\displaystyle\hskip 142.26378pt+2\,\alpha\frac{1-{\rm tr}(\boldsymbol{L}^{-\alpha})\,l_{i}^{\alpha}}{{\rm tr}^{2}(\boldsymbol{L}^{-\alpha})\,l_{i}^{2\alpha}}\Biggl{\}}\Biggl{)}$
		$\displaystyle=a_{o}\left(m\,q+b\,(m-q+1)\right)+2\,\alpha\,{E}\left(\frac{{\rm tr}(\boldsymbol{L}^{-2\alpha})}{{\rm tr}^{2}(\boldsymbol{L}^{-\alpha})}-1\right)\,.$

From the submultiplicativity of the trace for semi–definite positive matrices, we have ${\rm tr}(\boldsymbol{L}^{-2\alpha})\leq{\rm tr}^{2}(\boldsymbol{L}^{-\alpha})$ . Then, an upper bound for (9) is given by

\displaystyle{E}({\rm tr}(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha}))

\displaystyle\leq a_{o}\left(m\,q+b\,(m-q+1)\right)\,.

(11)

Secondly, dealing with ${E}\big{(}\!\ln|\Lambda(\hat{\boldsymbol{\Sigma}}_{\alpha}\boldsymbol{\Sigma}^{+})|\big{)}$ in (8), it can be shown that

\displaystyle\Lambda(\boldsymbol{\Sigma}^{+}\hat{\boldsymbol{\Sigma}}_{\alpha})=a_{o}\,\Lambda\bigr{(}\boldsymbol{\Sigma}^{+}\boldsymbol{H}\boldsymbol{L}\boldsymbol{\Phi}(\boldsymbol{L})\boldsymbol{H^{\top}}\bigl{)}=a_{o}\,\Lambda\bigr{(}\boldsymbol{L}^{1/2}\boldsymbol{H^{\top}}\boldsymbol{\Sigma^{+}}\boldsymbol{H}\boldsymbol{L}^{1/2}\boldsymbol{\Phi}(\boldsymbol{L})\bigl{)}.

Note that $\boldsymbol{L}^{1/2}\boldsymbol{H^{\top}}\boldsymbol{\Sigma^{+}}\boldsymbol{H}\boldsymbol{L}^{1/2}$ and $\boldsymbol{\Phi}(\boldsymbol{L})$ are full rank $q\times q$ matrices. It follows that

|\Lambda(\hat{\boldsymbol{\Sigma}}_{\alpha}\boldsymbol{\Sigma}^{+})|=a^{q}_{o}\,|\boldsymbol{L}^{1/2}\boldsymbol{H^{\top}}\boldsymbol{\Sigma^{+}}\boldsymbol{H}\boldsymbol{L}^{1/2}\boldsymbol{\Phi}(\boldsymbol{L})|.

Therefore

\displaystyle{E}(\ln|\Lambda(\hat{\boldsymbol{\Sigma}}_{\alpha}\boldsymbol{\Sigma}^{+})|)

\displaystyle=q\,\ln(a_{o})+{E}(\ln|\boldsymbol{L}^{1/2}\boldsymbol{H}^{T}\boldsymbol{\Sigma}^{+}\boldsymbol{H}\boldsymbol{L}^{1/2}|)+{E}(\ln|\boldsymbol{\Phi}(\boldsymbol{L})|).

(12)

Using the fact that $\ln(1+x)\geq 2\,x/(2+x)$ , for any positive constant $x$ , then

	$\displaystyle\ln\|\boldsymbol{\Phi}(\boldsymbol{L})\|=\ln\|\boldsymbol{I}_{q}+\frac{b\,\boldsymbol{L}^{-\alpha}}{{\rm tr}(\boldsymbol{L}^{-\alpha})}\|$	$\displaystyle=\sum_{i=1}^{q}\ln\left(1+\frac{b\,l^{-\alpha}_{i}}{{\rm tr}(\boldsymbol{L}^{-\alpha})}\right)$
		$\displaystyle\geq\sum_{i=1}^{q}\frac{2\,b\,l^{-\alpha}_{i}/{\rm tr}(\boldsymbol{L}^{-\alpha})}{2+b\,l^{-\alpha}_{i}/{\rm tr}(\boldsymbol{L}^{-\alpha})}.$

Thus

\displaystyle\ln|\boldsymbol{\Phi}(\boldsymbol{L})|\geq\frac{2b}{2+b}\,,

(13)

since, for $i=1,\dots,q$ , $l^{-\alpha}_{i}\leq{\rm tr}(\boldsymbol{L}^{-\alpha})$ . Consequently, thanks to (13), a lower bound for (12) is given by

\displaystyle{E}(\ln|\Lambda(\hat{\boldsymbol{\Sigma}}_{\alpha}\boldsymbol{\Sigma}^{+})|)\geq q\,\ln(a_{o})+{E}(\ln|\boldsymbol{L}^{1/2}\boldsymbol{H}^{T}\boldsymbol{\Sigma}^{+}\boldsymbol{H}\boldsymbol{L}^{1/2}|)+\frac{2\,b}{2+b}.

(14)

Now, relying on the proof of Proposition 2.1 in Tsukuma (2016), it can be shown that

\displaystyle R(\hat{\boldsymbol{\Sigma}}_{a_{o}},\boldsymbol{\Sigma})=-q\ln(a_{o})-{E}(\ln\!|\boldsymbol{L}^{1/2}\boldsymbol{H}^{T}\boldsymbol{\Sigma}^{+}\boldsymbol{H}\boldsymbol{L}^{1/2}|).

(15)

Finally, combining (11), (14) and (15), an upper bound for the risk difference in (7) is given by

\displaystyle\Delta_{(\alpha,a_{o})}

\displaystyle\leq a_{o}\,(m-q+1)\,b-\frac{2\,b}{2+b}=b\,\Bigr{(}a_{o}\,(m+q-1)\,(b+2)-2\Bigl{)},

since $a_{o}=1/m$ , which is non-positive as soon as

\displaystyle 0<b\leq b_{o}=\frac{2\,(q-1)}{m-q+1},

\qed

3 Numerical study

We study here numerically the performance of the proposed estimators of the form

\displaystyle\hat{\boldsymbol{\Sigma}}_{\alpha}=a_{o}\big{(}\boldsymbol{S}+\boldsymbol{H}\boldsymbol{L}\boldsymbol{\Psi}(\boldsymbol{L})\boldsymbol{H}^{\top}\big{)}\quad\text{with}\quad\alpha>0,\quad\boldsymbol{\Psi}(\boldsymbol{L})=b_{o}\frac{\boldsymbol{L}^{-\alpha}}{{\rm tr}(\boldsymbol{L}^{-\alpha})},

(16)

where $b_{o}$ is given in Proposition 2.

We consider the following structures of $\boldsymbol{\Sigma}$ : $(i)$ the identity matrix $\boldsymbol{I}_{p}$ and $(ii)$ an autoregressive structure with coefficient $0.9$ . We set their $p-r$ smallest eigenvalues to zero in order to construct matrices of rank $r\leq p$ .

To assess the performance of the proposed estimators, we compute the Percentage Reduction In Average Loss (PRIAL), for some values of $p$ , $n$ , $r$ and $\alpha$ , defined as

\displaystyle\textrm{PRIAL}(\hat{\boldsymbol{\Sigma}}_{\alpha})=\frac{R(\hat{\boldsymbol{\Sigma}}_{a_{o}},\boldsymbol{\Sigma})-{R(\hat{\boldsymbol{\Sigma}}_{\alpha},\boldsymbol{\Sigma})}}{R(\hat{\boldsymbol{\Sigma}}_{a_{o}},\boldsymbol{\Sigma})}\times 100,

where $\hat{\boldsymbol{\Sigma}}_{a_{o}}$ and $\hat{\boldsymbol{\Sigma}}_{\alpha}$ are respectively defined in (4) and (16).

$\boldsymbol{\Sigma}$	$(p,n)$	$r$	$\hat{\boldsymbol{\Sigma}}_{1}$	$\hat{\boldsymbol{\Sigma}}_{2}$	$\hat{\boldsymbol{\Sigma}}_{3}$	$\hat{\boldsymbol{\Sigma}}_{4}$	$\hat{\boldsymbol{\Sigma}}_{5}$
$(i)$	$(30,50)$	10	6.85	12.45	15.53	16.95	17.53
		20	9.20	13.91	14.88	14.68	12.47
		30	11.81	14.33	13.41	12.43	11.71
	$(50,30)$	20	18.31	19.65	17.75	16.44	15.63
		40	17.12	16.33	14.07	12.78	12.02
		50	11.80	14.23	13.29	12.30	11.59
	$(150,30)$	20	18.17	19.69	17.87	16.56	15.71
		40	17.08	16.31	14.07	12.76	11.99
		60	8.88	12.27	12.33	11.75	11.19
		150	2.83	5.09	6.50	7.25	7.61
$(ii)$	$(30,50)$	10	6.06	8.70	9.62	9.89	9.92
		20	8.81	11.66	12.12	11.93	11.61
		30	11.46	13.15	12.35	11.48	10.82
	$(50,30)$	20	17.18	17.45	15.85	14.76	14.07
		40	16.34	15.18	13.19	12.00	11.28
		50	11.33	12.82	12.02	11.19	10.57
	$(150,30)$	20	17.30	18.00	16.38	15.19	14.42
		40	16.27	15.04	13.04	11.84	11.13
		60	8.48	10.43	10.29	9.81	9.35
		150	2.73	4.03	4.61	4.86	4.95

Table 1: Effect of

\alpha=1,\dots,5

on PRIAL’s for the structures

(i)

and

(ii)

\boldsymbol{\Sigma}

Table 1 shows that the proposed estimators improve over $\hat{\boldsymbol{\Sigma}}_{a_{o}}$ for any possible ordering of $p,n$ and $r$ . Compared to other cases, the Haff type estimators $\hat{\boldsymbol{\Sigma}}_{\alpha}$ (for $\alpha=1,\dots,5$ ) have better performances in the case where $p>n>r$ , with PRIAL’s higher than $14.07\%$ for both structures $(i)$ and $(ii)$ of $\boldsymbol{\Sigma}$ . We report that the optimal value of $\alpha$ , which maximizes the PRIAL’s, depends on $p,n$ and $r$ .

Acknowledgement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

References

Srivastava (2003) M. S. Srivastava, Singular Wishart and multivariate Beta distributions, Ann. Statis. 31 (2003) 1537–1560.
Stein (1986) C. Stein, Lectures on the theory of estimation of many parameters, J. Sov. Math. 34 (1986) 1373–1403.
Tsukuma (2016) H. Tsukuma, Estimation of a high-dimensional covariance matrix with the Stein loss, J. Multivar. Anal. 148 (2016) 1–17.
Ledoit and Wolf (2004) O. Ledoit, M. Wolf, A well-conditioned estimator for large-dimensional covariance matrices, J. Multivar. Anal. 88 (2004) 365–411.
James and Stein (1961) W. James, C. Stein, Estimation with quadratic loss, in: Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics, Berkeley, California, 1961, pp. 361–379.
Takemura (1984) A. Takemura, An orthogonally invariant minimax estimator of the covariance matrix of a multivariate normal population, Tsukuba J. Math. 8 (1984) 367–376.
Dey and Srinivasan (1985) D. K. Dey, C. Srinivasan, Estimation of a covariance matrix under Stein’s loss, Ann. of Statis. 13 (1985) 1581–1591.
Haff (1980) L. Haff, Empirical Bayes estimation of the multivariate normal covariance matrix, Ann. Statis. 8 (1980) 586–597.
Konno (2009) Y. Konno, Shrinkage estimators for large covariance matrices in multivariate real and complex normal distributions under an invariant quadratic loss, J. Multivar. Anal. 100 (2009) 2237–2253.
Haddouche et al. (2021) A. M. Haddouche, D. Fourdrinier, F. Mezoued, Scale matrix estimation of an elliptically symmetric distribution in high and low dimensions, J. Multivar. Anal. 181 (2021) 104680.
Chételat and Wells (2016) D. Chételat, M. T. Wells, Improved second order estimation in the singular multivariate normal model, J. Multivar. Anal. 147 (2016) 1–19.
Kubokawa and Srivastava (2008) T. Kubokawa, M. Srivastava, Estimation of the precision matrix of a singular Wishart distribution and its application in high-dimensional data, J. Multivar. Anal. 99 (2008) 1906–1928.