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

Haddouche, Anis M.

Lu, Wei

Comptes Rendus Mathématique

English

arXiv

Anis  M. Haddouche

Wei Lu

Comptes RendusMathématiqueAnis M. Haddouche and Wei LuA unified approach for covariance matrix estimation under Stein lossVolume 360 (2022), p. 1093-1098https://doi.org/10.5802/crmath.356This article is licensed under theCreative Commons Attribution 4.0 International License.http://creativecommons.org/licenses/by/4.0/Les Comptes Rendus. Mathématique sont membres duCentre Mersenne pour l’édition scientifique ouvertewww.centre-mersenne.orge-ISSN : 1778-3569Comptes RendusMathématique2022, Vol. 360, p. 1093-1098https://doi.org/10.5802/crmath.356Statistics / StatistiquesA unified approach for covariance matrixestimation under Stein lossAnis M. Haddouche ∗, a andWei Luaa INSA Rouen, Normandie Univ, UNIROUEN, UNIHAVRE, LITIS and LMI, avenue del’Université, BP 8, Saint-Étienne-du-Rouvray, 76801, FranceE-mails: mohamed.haddouche@insa-rouen.fr (A. M. Haddouche),wei.lu@insa-rouen.fr (W. Lu)Abstract. In this paper, we address the problem of estimating a covariance matrix of a multivariate Gaussiandistribution, from a decision theoretic point of view, relative to a Stein type loss function. We investigate thecase where the covariance matrix is invertible and the case when it is non–invertible in a unified approach.Résumé. Dans cet article, nous abordons le problème de l’estimation d’une matrice de covariance d’unedistribution gaussienne multivariée, du point de vue de la théorie de la décision, par rapport à une fonction decoût de type Stein. Nous étudions dans une approche unifiée le cas où la matrice de covariance est inversibleet le cas où elle n’est pas inversible.Mathematical subject classification (2010). 62H12, 62F10, 62C99.Manuscript received 13 July 2021, revised 4 January 2022, accepted 13 March 2022.1. IntroductionLet X be an observed p ×n matrix of the formX = B Z , (1)where B is a p × r matrix of unknown parameters, with r ≤ p, and Z is a r × n randommatrix. Assume that r is known and that the columns of Z are identically and independentlydistributed as the r -dimensional multivariate normal distribution Nr (0r , I r ). Then the columnsof X are identically and independently distributed from the p-dimensional multivariate normalNp (0p ,Σ), where Σ= B B T is the unknown p ×p covariance matrix withrank(Σ) = r ≤ p .It follows that the p × p sample covariance matrix S = X X T has a singular Wishart distribution(see Srivastava [9]) such thatrank(S) = min(n,r ) = q ≤ p ,∗Corresponding author.ISSN (electronic) : 1778-3569 https://comptes-rendus.academie-sciences.fr/mathematique/1094 Anis M. Haddouche and Wei Luwith probability one. Note that, as the matrices S and Σ are non–invertible, we denote in thefollowing by S+ and Σ+ the Moore-Penrose inverses of S and Σ respectively.We consider the problem of estimating the covariance matrix Σ under the Stein type lossfunctionL(Σ̂,Σ)= tr(Σ+Σ̂)− ln∣∣Λ(Σ+Σ̂)∣∣−q , (2)where Σ̂ estimates Σ and Λ(Σ+Σ̂) is the diagonal matrix of the q positives eigenvalues of Σ+Σ̂.Note that this loss function is based on the Kullback–Leibler divergence between two multivariatenormal distributions with the same means and with covariances Σ̂ and Σ. It is an extension of theusual Stein loss function (see Stein [10]). Thus, the corresponding risk function is given byR(Σ̂,Σ)= E [L (Σ̂,Σ)] ,where E(·) 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 [11]) to the context of the model (1)(see Tsukuma [13] for more details).As related by Ledoit and Wolf [8], the difficulty of covariance estimation is commonly charac-terized by the ratio p/n. In fact, the usual estimators of the formΣ̂a = a S with a > 0, (3)perform poorly when n, p → ∞ with p/n → c > 0. Indeed, the larger (smaller) eigenvalues ofΣ are overestimated (underestimated) by those estimators. Hence, in this asymptotic context,alternative estimators dominating (with lower risk functions) these usual estimators are needed.A possible approach is to regularize the eigenvalues of Σ̂a , which gives rise to the class oforthogonally invariant estimators (see Takemura [12]) as those in (5) below.Note that, as we consider the model (1), we deal in a unified approach, with the followingcases.(i) n < r = p: Σ is invertible of rank p and S is non–invertible of rank n;(ii) r = p ≤ n: Σ and S are invertible;(iii) r < p ≤ n: Σ and S are non–invertible of rank r ;(iv) r ≤ n < p: Σ and S are non–invertible of rank r ;(v) n < r < p: Σ and 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 [11], Dey and Srivastava [2] and Haff [5] for the case (i), Konno [6] and Haddouche etal. [4] for the cases (i) and (ii). See also Chételat and Wells [1] for the cases (iii) and (iv). Recently,Tsukuma and Kubokawa [13] extended the Stein [11] estimator to the five possible cases, citedabove, in a unified approach. Similarly, we extend here the class of Haff [5] estimators to thecontext of the model (1).The rest of this paper is organized as follows. In Section 2, we derive the domination result ofthe Haff estimators over the usual estimators. We study numerically the behavior of the proposedestimators in Section 3.2. Main resultImproving the class of the natural estimators in (3) relies on improving the optimal estimatoramong this class, that is, the one which minimizes the loss function (2).Proposition 1 (Tsukuma [13]). Under the Stein type loss function (2), the optimal estimatoramong the class (3) is given byΣ̂ao = ao S, where ao =1mand m = max(n,r ). (4)Anis M. Haddouche and Wei Lu 1095As mentioned in Section 1, we consider the class of orthogonally invariant estimators asalternative estimators. Let S = H L H> be the eigenvalue decomposition of S where H is a p × qsemi–orthogonal matrix of eigenvectors and L = diag(l1, . . . , lq ), with l1 >, . . . ,> lq , is the diagonalmatrix of the q positive corresponding eigenvalues (see Kubokawa and Srivastava [7] for moredetails). The class of orthogonally invariant estimators is of the formΣ̂Ψ = ao(S +H LΨ(L) H>) (5)with Ψ(L) = diag(ψ1(L), . . . ,ψq (L)), where ψi (L) (i = 1, . . . , q) is a differentiable function of L.More precisely, we consider an extension of the empirical Bayes Haff type estimators in [5], to thecontext of the model (1), defined asΣ̂α = ao(S +HLΨ(L)H>) where Ψ(L) = b L−αtr(L−α), with α≥ 1 and b > 0, (6)where ao is given in (4). Note that, for α = 1, this is the estimator considered by Konno [6],who deals with the cases (i) and (ii) under a quadratic loss. An extension to elliptical settingwas considered recently by Haddouche et al. [4]. Note also that Tsukuma [13] considered alsoempirical Bayes estimators but not the one in (6). We give in the following theorem our mainresult.Theorem 2. Let Ψ(L) = diag(ψ1(L), . . . ,ψq (L)) where, for any i = 1, . . . , q, ψi (L) = blαi /tr(L−α)with α ≥ 1 and b > 0, be a differentiable function of L. The Haff type estimators in (6) dominatesthe optimal estimator in (4), under the loss function (2), as soon as0 < b ≤ bo =2(q −1)m −q +1 .Proof. The proof consists to showing that the risk difference between the Haff type estimatorsin (6) and the optimal estimator in (4), namely,∆(α,ao ) = R(Σ̂α,Σ)−R (Σ̂ao ,Σ) , (7)is non–positive. Note that Σ̂α can be written asΣ̂α = ao HLΦ(L)H> with Φ(L) = I q +Ψ(L).The risk of these estimators under the Stein type loss function (2) is given byR(Σ̂α,Σ)= E (tr(Σ+Σ̂α))−E (ln∣∣Λ(Σ+Σ̂α)∣∣)−q . (8)First, dealing with E(tr(Σ̂αΣ+)), we apply in Tsukuma [13, Lemma A.2] in order to get rid of theunknown parameter Σ+. It follows that,E(tr(Σ+Σ̂α))= aoE(q∑i=1{(m −q +1)ϕi +2li ∂ϕi∂li+2q∑j > iliϕi − l jϕ jli − l j})(9)where, for any i = 1, . . . , q ,φi = 1+bl−αitr(L−α),∂ϕi∂li= bα 1− tr(L−α) lαitr2 (L−α) l 1+2αiandq∑j > iliϕi − l jϕ jli − l j=q∑j > i(1+ btr(L−α){l 1−αi − l 1−αjli − l j}).Using the fact, for any j > i , l j > li , it can be shown thatq∑j > iliϕi − l jϕ jli − li≤ (q − i ) . (10)1096 Anis M. Haddouche and Wei LuTherefore, using (10), we obtainE(tr(Σ+Σ̂α))≤ ao E(q∑i=1{(m −q +1)(1+b li −αtr(L−α))+2bα 1− tr(L−α) lαitr2 (L−α) l 2αi+2(q − i )})= ao m q +ao b E(q∑i=1{(m −q +1) li −αtr(L−α)+2α1− tr(L−α) lαitr2 (L−α) l 2αi})= ao(m q +b (m −q +1))+2αE ( tr(L−2α)tr2 (L−α)−1).Now, from the submultiplicativity of the trace for semi–definite positive matrices, we havetr(L−2α) ≤ tr2(L−α). Then, an upper bound for (9) is given byE(tr(Σ+Σ̂α))≤ ao (m q +b (m −q +1)) . (11)Secondly, dealing with E(ln |Λ(Σ̂αΣ+)|) in (8), it can be shown thatΛ(Σ+Σ̂α)= aoΛ(Σ+HLΦ(L)H>)= aoΛ(L1/2H>Σ+HL1/2Φ(L)).Note that L1/2H>Σ+HL1/2 andΦ(L) are full rank q ×q matrices. It follows that∣∣Λ(Σ̂αΣ+)∣∣= aqo ∣∣L1/2H>Σ+HL1/2Φ(L)∣∣ .ThereforeE(ln∣∣Λ(Σ̂αΣ+)∣∣)= q ln(ao)+E (ln ∣∣L1/2H TΣ+HL1/2∣∣)+E (ln |Φ(L)|) . (12)Using the fact that ln(1+x) ≥ 2 x/(2+x), for any positive constant x, we haveln |Φ(L)| = ln∣∣∣∣I q + b L−αtr(L−α)∣∣∣∣=q∑i=1ln(1+ b l−αitr(L−α))≥q∑i=12b l−αi /tr(L−α)2+b l−αi /tr(L−α).Thenln |Φ(L)| ≥ 2b2+b , (13)since, for any i = 1, . . . , q , l−αi ≤ tr(L−α). Consequently, thanks to (13), a lower bound for (12) isgiven byE(ln∣∣Λ(Σ̂αΣ+)∣∣)≥ q ln(ao)+E (ln ∣∣L1/2H TΣ+HL1/2∣∣)+ 2b2+b . (14)Now, relying on the proof of [13, Proposition 2.1] in Tsukuma, it can be shown thatR(Σ̂ao ,Σ)=−q ln(ao)−E (ln∣∣L1/2H TΣ+HL1/2∣∣) . (15)Finally, combining (11), (14) and (15), an upper bound for the risk difference (7) is given by∆(α,ao ) ≤ ao(m −q +1) b − 2b2+b= bb +2( (m −q −1) (b +2)m−2),since ao = 1/m, which is non-positive as soon as0 < b ≤ bo =2(q −1)m −q +1 . Anis M. Haddouche and Wei Lu 10973. Numerical studyWe study here numerically the performance of the Haff type estimators in (6). As shown numeri-cally by Fourdrinier et al. [3] that the best constant b is bo . So that, we consider in this section thefollowing Haff type estimatorsΣ̂α = ao(S +HLΨ(L)H>) where Ψ(L) = bo L−αtr(L−α),with α≥ 1 and bo =2(q −1)m −q +1 . (16)As for the population covariance matrixΣ, we consider the following structures: (i) the identitymatrix I p and (ii) an autoregressive structure with coefficient 0.9. We set their p − r smallesteigenvalues to zero in order to construct matrices of rank r ≤ p.To assess the performance of the proposed estimators, we compute the Percentage ReductionIn Average Loss (PRIAL), for some values of p, n, r and α, defined asPRIAL(Σ̂α)= R (Σ̂ao ,Σ)−R (Σ̂α,Σ)R(Σ̂ao ,Σ) ×100,where Σ̂ao and Σ̂α are respectively defined in (4) and (3).Table 1. Effect of α on PRIAL’s for the structures(i) and (ii) of Σ.Σ (p,n) r Σ̂1 Σ̂2 Σ̂3 Σ̂4 Σ̂5(i)(30,50)10 6.85 12.45 15.53 16.95 17.5320 9.20 13.91 14.88 14.68 12.4730 11.81 14.33 13.41 12.43 11.71(50,30)20 18.31 19.65 17.75 16.44 15.6340 17.12 16.33 14.07 12.78 12.0250 11.80 14.23 13.29 12.30 11.59(150,30)20 18.17 19.69 17.87 16.56 15.7140 17.08 16.31 14.07 12.76 11.9960 8.88 12.27 12.33 11.75 11.19150 2.83 5.09 6.50 7.25 7.61(ii)(30,50)10 6.06 8.70 9.62 9.89 9.9220 8.81 11.66 12.12 11.93 11.6130 11.46 13.15 12.35 11.48 10.82(50,30)20 17.18 17.45 15.85 14.76 14.0740 16.34 15.18 13.19 12.00 11.2850 11.33 12.82 12.02 11.19 10.57(150,30)20 17.30 18.00 16.38 15.19 14.4240 16.27 15.04 13.04 11.84 11.1360 8.48 10.43 10.29 9.81 9.35150 2.73 4.03 4.61 4.86 4.95Table 1 shows that the proposed estimators improve over Σ̂ao for any possible ordering ofp,n and r . Compared to other cases, the Haff type estimators Σ̂α (for α = 1, . . . ,5) have betterperformances in the case where p > n > r , with PRIAL’s higher than 14.07% for both structures (i)and (ii) ofΣ. We report that the optimal value ofα, which maximizes the PRIAL’s, depends on p,nand r .1098 Anis M. Haddouche and Wei LuWe then compare the Prial’s of the Haff type estimators in and the well known James–Steintype estimator Σ̂JS considered by Tsukuma [13], given byΣ̂JS = HLD JSq H>,where D JSq = diag(dJS1 , . . . , dJSq ) with, for i = 1, . . . , q , d JSi = (n + r −2i +1)−1.Table 2. PRIAL’s of the Haff type estimators and the Stein type estimator for the struc-tures (ii) of Σ.(p,n) r Σ̂1 Σ̂2 Σ̂3 Σ̂4 Σ̂JS(30,40)7 9.06 12.34 13.36 13.68 9.8714 9.60 13.09 13.78 13.74 20.45(40,30)7 11.87 15.57 16.41 16.55 16.0414 13.91 16.41 16.13 15.55 26.35(20,100)5 4.04 5.61 6.14 6.29 4.6815 3.67 5.45 6.38 7.12 11.80Table 2 shows that in some cases the PRIALs of the Haff-type estimators are higher thanthose of the Stein-type estimator, while in other cases the Stein-type estimator has a betterperformance. Consequently, their risk functions intersect and thus, in our opinion, a theoreticaldominance results between these two estimators cannot be derived.References[1] D. Chételat, M. T. Wells, “Improved second order estimation in the singular multivariate normal model”, J. Multivari-ate Anal. 147 (2016), p. 1-19.[2] D. K. Dey, C. Srinivasan, “Estimation of a covariance matrix under Stein’s loss”, Ann. Stat. 13 (1985), p. 1581-1591.[3] D. Fourdrinier, A. M. Haddouche, F. Mezoued, “Covariance matrix estimation under data-based loss”, Stat. Probab.Lett. 177 (2021), article no. 109160.[4] A. M. Haddouche, D. Fourdrinier, F. Mezoued, “Scale matrix estimation of an elliptically symmetric distribution inhigh and low dimensions”, J. Multivariate Anal. 181 (2021), article no. 104680.[5] L. R. Haff, “Empirical Bayes estimation of the multivariate normal covariance matrix”, Ann. Stat. 8 (1980), no. 3,p. 586-597.[6] Y. Konno, “Shrinkage estimators for large covariance matrices in multivariate real and complex normal distributionsunder an invariant quadratic loss”, J. Multivariate Anal. 100 (2009), no. 10, p. 2237-2253.[7] T. Kubokawa, M. S. Srivastava, “Estimation of the precision matrix of a singular Wishart distribution and its applica-tion in high-dimensional data”, J. Multivariate Anal. 99 (2008), no. 9, p. 1906-1928.[8] O. Ledoit, M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices”, J. Multivariate Anal.88 (2004), no. 2, p. 365-411.[9] M. S. Srivastava, “Singular Wishart and multivariate Beta distributions”, Ann. Stat. 31 (2003), no. 5, p. 1537-1560.[10] C. M. Stein, “Estimation of a covariance matrix”, in Rietz Lecture, 39th Annual Meeting of the IMS, Atlanta, GA, 1975,1975.[11] ——— , “Lectures on the theory of estimation of many parameters”, J. Sov. Math. 34 (1986), no. 1, p. 1373-1403.[12] A. Takemura, “An orthogonally invariant minimax estimator of the covariance matrix of a multivariate normalpopulation”, Tsukuba J. Math. 8 (1984), p. 367-376.[13] H. Tsukuma, “Estimation of a high-dimensional covariance matrix with the Stein Loss”, J. Multivariate Anal. 148(2016), p. 1-17.

A unified approach for covariance matrix estimation under Stein loss

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

Haddouche, Anis  M.

Directory of Open Access Journals

Numérisation de Documents Anciens Mathématiques

http://www.numdam.org/item/10.5802/crmath.356.pdf

A unified approach for covariance matrix estimation under Stein loss

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mathématique

Directory of Open Access Journals

Numérisation de Documents Anciens Mathématiques