Random Matrix Improved Covariance Estimation for a Large Class of Metrics

TL;DR

This paper introduces an improved covariance estimation method using random matrix theory that outperforms traditional estimates and is computationally efficient, with demonstrated benefits in machine learning classification tasks.

Contribution

It proposes a novel covariance estimation technique based on random matrix theory applicable to a broad class of metrics, improving accuracy and computational simplicity.

Findings

Outperforms sample covariance matrix estimates in various metrics.

Competes with state-of-the-art methods while being simpler to compute.

Enhances classification performance in discriminant analysis.

Abstract

Relying on recent advances in statistical estimation of covariance distances based on random matrix theory, this article proposes an improved covariance and precision matrix estimation for a wide family of metrics. The method is shown to largely outperform the sample covariance matrix estimate and to compete with state-of-the-art methods, while at the same time being computationally simpler. Applications to linear and quadratic discriminant analyses also demonstrate significant gains, therefore suggesting practical interest to statistical machine learning.

Tables

Table 1. Table 1: Distances d 𝑑 d and divergences δ 𝛿 \delta , and their corresponding f ​ ( z ) 𝑓 𝑧 f(z) functions.

Divergences	$f\left(z\right)$
$d_{\text{R}}^2$	${\log }^2\left(z\right)$
$d_{\text{B}}^2$	$-\frac{1}{4}\log \left(z\right)+\frac{1}{2}\log \left(1+z\right)-\frac{1}{2}\log \left(2\right)$
${\delta }_{\text{KL}}$	$\frac{1}{2}z-\frac{1}{2}\log \left(z\right)-\frac{1}{2}$
${\delta }_{\alpha \text{R}}$	$\frac{-1}{2\left(\alpha -1\right)}\log \left(\alpha +\left(1-\alpha \right)z\right)+\frac{1}{2}\log \left(z\right)$

Table 2. Table 2: Values of G ​ ( z ) 𝐺 𝑧 G(z) and F ​ ( z ) 𝐹 𝑧 F(z) for “atomic” f ​ ( z ) 𝑓 𝑧 f(z) functions used in most distances and divergences under study; here s > 0 𝑠 0 s>0 and z ∈ ℂ 𝑧 ℂ z\in{\mathbb{C}} .

$f\left(z\right)$	$G\left(z\right)$
${\log }^2\left(z\right)$	$z\left({\log }^2\left(z\right)-2\log \left(z\right)+2\right)$
$\log \left(z\right)$	$-z\log \left(z\right)+z$
$\log \left(1+sz\right)$	$s\log \left(s+z\right)+z\log \left(\frac{s+z}{z}\right)$
$z$	$\log \left(z\right)$

Equations

$$C\equiv {\mathrm{argmin}}_{M\succ 0}\delta (M,C)$$

$$\delta (M,C)$$

$$m_{\theta }(z)=\int \frac{\theta (dt)}{t-z}.$$

$${\mu }_p\equiv \frac{1}{p}\sum _{i=1}^p{\delta }_{{\lambda }_i(M^{-1}\hat{C})}\text{ and }{\nu }_p\equiv \frac{1}{p}\sum _{i=1}^p{\delta }_{{\lambda }_i(M^{-1}C)}.$$

$$\delta (M,C)=\frac{1}{2\pi ı}{\oint }_{\Gamma }f(z)m_{{\nu }_p}(z;M)dz$$

$$\delta (M,C)-\hat{\delta }(M,X)\to 0$$

$$\hat{\delta }(M,X)\equiv \frac{1}{2\pi ıc}{\oint }_{\hat{\Gamma }}G\left(-m_{{\tilde{\mu }}_p}(z;M)\right)dz$$

$${\mathrm{argmin}}_{M\succ 0}\{h_X(M)\},\text{ where }{h}_X(M)\equiv (\hat{\delta }(M,X){)}^2.$$

$$\mathrm{D}{h}_X(M)[\xi ]=⟨\nabla h_X(M),\xi {⟩}_M^{S_n^{++}}$$

$$⟨\eta ,\xi {⟩}_M^{S_n^{++}}=\mathrm{tr}\left(M^{-1}\eta M^{-1}\xi \right).$$

$$\mathrm{D}{h}_X(M)[\xi ]$$

$$=\frac{-\hat{\delta }(M,X)}{\pi ic}{\oint }_{\Gamma }g(-m_{\widetilde{{\mu }_p}}\left(z,M\right))\mathrm{D}{m}_{\widetilde{{\mu }_p}}\left(z,M\right)[\xi ]dz.$$

$$\mathrm{D}{m}_{\widetilde{{\mu }_p}}\left(z,M\right)[\xi ]$$

$$=\frac{c}{p}\mathrm{D}\mathrm{tr}\left({\left[M^{-1}\hat{C}-z{I}_p\right]}^{-1}\right)[\xi ]$$

$$=\frac{c}{p}\mathrm{tr}\left(M^{-1}\hat{C}{\left[M^{-1}\hat{C}-z{I}_p\right]}^{-2}{M}^{-1}\xi \right)$$

$$={⟨\frac{c}{p}\mathrm{sym}\left(\hat{C}{\left[M^{-1}\hat{C}-z{I}_p\right]}^{-2}\right),\xi ⟩}_M^{S_n^{++}}$$

$$-ı\pi p\frac{\nabla h_X(M)}{\hat{\delta }(M,X)}$$

$$={\oint }_{\hat{\Gamma }}g\left(-m_{{\tilde{\mu }}_p}(z;M)\right)\mathrm{sym}\left(\hat{C}(M^{-1}\hat{C}-z{I}_p{)}^{-2}\right)dz$$

$${\mathbb{R}}_{+}$$

$$t$$

$$M_{k+1}$$

$$\delta (M^{(p)},C^{(p)})-\hat{\delta }(M^{(p)},X^{(p)})\to 0$$

$$[{\omega }_{k+1}{]}_i=[{\omega }_k{]}_i\exp \left(-t\frac{[{\delta }_k{]}_i}{[{\omega }_k{]}_i}\right)$$

$$M_k$$

$$M_{k+1}$$

$$+\frac{t^2}{2}\nabla h_X(M_k)M_k^{-1}\nabla h_X(M_k)+O(t^3).$$

$$\delta (M,C^{-1})-{\hat{\delta }}^{\mathrm{inv}}(M,X)\to 0$$

$${\hat{\delta }}^{\mathrm{inv}}(M,X)\equiv \frac{1}{2\pi ıc}{\oint }_{\hat{\Gamma }}F\left(-m_{{\tilde{\mu }}_p^{\mathrm{inv}}}(z;M)\right)dz$$

$$ı\pi p\frac{\nabla h_X^{\mathrm{inv}}(M)}{{\hat{\delta }}^{\mathrm{inv}}(M,X)}$$

$$={\oint }_{\hat{\Gamma }}f\left(-m_{{\tilde{\mu }}_p^{\mathrm{inv}}}(z;M)\right)\mathrm{sym}\left(M\hat{C}(M\hat{C}-z{I}_p{)}^{-2}M\right)dz.$$

$$\nabla h_X(M)\equiv 2\hat{\delta }(M,X)\cdot \mathrm{sym}\left(\hat{C}\cdot V{\Lambda }_{\nabla }{V}^{-1}\right)$$

$$\Lambda -\frac{1}{n}\sqrt{\lambda }{\sqrt{\lambda }}^{\mathsf{T}}$$

$$m_{{\tilde{\mu }}_p}(t)=-s.$$

$$[{\Lambda }_{\nabla }{]}_{kk}$$

Full text

covariance matrix, precision matrix, LDA, QDA

1 Introduction

Covariance and precision matrix estimation is a fundamental and simply posed, yet still largely considered, key problems of statistical data analysis, with countless applications in statistical inference. In machine learning, it is notably at the core of elementary methods as linear (LDA) and quadratic discriminant analysis (QDA) (McLachlan, 2004).

Estimation of the covariance matrix $C\in{\mathbb{R}}^{p\times p}$ based on $n$ independent (say zero mean) samples $x_{1},\ldots,x_{n}\in{\mathbb{R}}^{p}$ is conventionally performed using the sample covariance matrix (SCM) $\hat{C}\equiv\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\sf T}$ (and its inverse using $\hat{C}^{-1}$). The estimate is however only consistent for $n\gg p$ and only invertible for $n\geq p$. Treating the important practical cases where $n\sim p$ and even $n\ll p$ has recently spurred a series of parallel lines of research. These directions rely either on structural constraints, such as “toeplitzification” procedures for Toeplitz covariance models (particularly convenient for time series) (Bickel et al., 2008; Wu & Pourahmadi, 2009; Vinogradova et al., 2015), on sparse constraints with LASSO and graphical LASSO-based approaches (Friedman et al., 2008) or, more interestingly for the present article, on exploiting the statistical independence in the entries of the vectors $x_{i}$.

(Ledoit & Wolf, 2004) proposes to linearly “shrink” $\hat{C}$ as $\hat{C}(\rho)\equiv\rho I_{p}+\sqrt{1-\rho^{2}}\hat{C}$ for $\rho>0$ chosen to minimize the expected Frobenius distance ${\rm E}[\|C-\hat{C}(\rho)\|_{F}]$ in the asymptotic $p,n\to\infty$ limit with $p/n\to c>0$. Basic results from random matrix theory (RMT) are used here to estimate $\rho$ consistently. This procedure is simple and quite flexible and has been generalized in various directions (e.g., in (Couillet & McKay, 2014) with a robust statistics approach). However, the method only applies a naive homothetic map to each $\lambda_{i}(\hat{C})$ of $\hat{C}$ in order to better estimate $\lambda_{i}(C)$. A strong hope to recover a better approximation of the $\lambda_{i}(C)$’s then arose from (Silverstein & Bai, 1995; Silverstein & Choi, 1995) that provide a random matrix result relating directly the limiting eigenvalue distributions of $C$ and $\hat{C}$. Unfortunately, while estimating the $\lambda_{i}(\hat{C})$’s from the $\lambda_{i}(C)$’s is somewhat immediate, estimating the $\lambda_{i}(C)$’s backward from the $\lambda_{i}(\hat{C})$’s is a difficult task. (El Karoui et al., 2008) first proposed an optimization algorithm to numerically solve this problem, however with little success as the method is quite unstable and has rarely been efficiently reproduced. (Mestre, 2008) later offered a powerful idea, based on contour integral, to consistently estimate linear functionals $\frac{1}{n}\sum_{i=1}^{n}f(\lambda_{i}(C))$ from the $\lambda_{i}(\hat{C})$’s. But $f$ is constrained to be very smooth (complex analytic) which prevents the estimation of the individual $\lambda_{i}(C)$’s. Recently, Ledoit and Wolf took over the work of El Karoui, which they engineered to obtain a more efficient numerical method, named QuEST (Ledoit & Wolf, 2015). Rather than inverting the Bai–Silverstein equations, the authors also proposed, with the same approach, to estimate the $\lambda_{i}(C)$’s by minimizing a Frobenius norm distance (Ledoit & Wolf, 2015) (named QuEST1 in the present article) or a Stein loss (Ledoit et al., 2018) (QuEST2 here).

These methods, although more stable than El Karoui’s initial approach, however suffer several shortcomings: (i) they are still algorithmically involved as they rely on a series of fine-tuned optimization schemes, and (ii) they are only adaptable to few error metrics (Frobenius, Stein).

Inspired by Mestre’s approach and the recent work (Couillet et al., 2018), this article proposes a different procedure consisting in (i) writing $C$ as the solution to $\operatorname{argmin}_{M\succ 0}\delta(M,C)$ for a wide range of metrics $\delta$ (Fisher, Batthacharyya, Stein’s loss, Wasserstein, etc.), (ii) based on (Couillet et al., 2018), using the fact that $\delta(M,C)-\hat{\delta}(M,X)\to 0$ for some consistent estimator $\hat{\delta}$, valid for all deterministic $M$ and samples $X=[x_{1},\ldots,x_{n}]\in{\mathbb{R}}^{p\times n}$ having zero mean and covariance $C$, and (iii) proceeding to a gradient descent on $\hat{\delta}$ rather than on the unknown $\delta$ itself. With appropriate adaptations, the estimation of $C^{-1}$ is similarly proposed by solving instead $\operatorname{argmin}_{M\succ 0}\delta(M,C^{-1})$.

While only theoretically valid for matrices $M$ independent of $X$, and thus only in the first steps of the gradient descent, the proposed method has several advantages: (i) it is easy to implement and technically simpler than QuEST, (ii) it is adaptable to a large family of distances and divergences, and, most importantly, (iii) simulations suggest that it systematically outperforms the SCM and is competitive with, if not better than, QuEST.

The remainder of the article is organized as follows. Section 2 introduces preliminary notions and concepts on which are hinged our proposed algorithms, thereafter described in Section 3. Section 4 provides experimental validations and applications, including an improved version of LDA/QDA based on the proposed enhanced estimates.

Reproducibility. Matlab codes for the proposed estimation algorithms are available as supplementary materials and are based on Manopt, a Matlab toolbox for optimization on manifolds (Boumal et al., 2014).

2 Preliminaries

Let $n$ random vectors $x_{1},\ldots,x_{n}\in{\mathbb{R}}^{p}$ be of the form ${x_{i}=C^{\frac{1}{2}}z_{i}}$ for some positive definite matrix ${C\in{\mathbb{R}}^{p\times p}}$ and $z_{1},\ldots,z_{n}\in{\mathbb{R}}^{p}$ independent random vectors of independent entries with ${\rm E}[[z_{i}]_{j}]=0$ and ${\rm E}[|[z_{i}]_{j}|^{2}]=1$. We further assume the following large dimensional regime for $n$ and $p$.

Assumption 1** (Growth Rate).**

As $n\to\infty$, ${p/n\to c\in(0,1)}$ and $\limsup_{p}\max\{\|C^{-1}\|,\|C\|\}<\infty$ for $\|\cdot\|$ the matrix operator norm.

Our objective is to estimate $C$ and $C^{-1}$ based on $x_{1},\ldots,x_{n}$ under the above large $p,n$ regime. For simplicity of exposition and readability, we mostly focus on the estimation of $C$ and more briefly discuss that of $C^{-1}$.

Our approach relies on the following elementary idea:

[TABLE]

where, for some function $f$,

[TABLE]

is a divergence (possibly a squared distance $\delta=d^{2}$) between the positive definite matrices $M$ and $C$, depending only on the eigenvalues of $M^{-1}C$. Among divergences satisfying this condition, we find the natural Riemannian distance $d_{\textup{R}}^{2}$ (Bhatia, 2009), which corresponds to the Fisher metric for the multivariate normal distribution (Skovgaard, 1984); the Battacharyya distance $d_{\textup{B}}^{2}$ (Sra, 2013), which is close to the natural Riemannian distance while numerically less expensive; the Kullback-Leibler divergence $\delta_{\textup{KL}}$, linked to the likelihood and studied for example in (Moakher, 2012); the Rényi divergence $\delta_{\alpha\textup{R}}$ for Gaussian $x_{i}$’s (Van Erven & Harremos, 2014); etc.111The Frobenius distance does not fall into this setting but has already largely been investigated and optimized. Table 1 reports the explicit values of $f$ for these divergences.

Since $\delta(M,C)$ is not accessible as $C$ is unknown, our approach exploits an estimator for $\delta(M,C)$ which is consistent in the large $n,p$ regime of Assumption 1.

To this end, our technical arguments are fundamentally based on random matrix theory, and notably rely on the so-called Stieljes transform of eigenvalue distributions. For an arbitrary real-supported probability measure $\theta$, the Stieltjes transform $m_{\theta}:{\mathbb{C}}\setminus{\rm supp}(\theta)\to{\mathbb{C}}$ is defined as

[TABLE]

The key interest of the Stieltjes transform in this article is that it allows one to relate the distributions of the eigenvalues of $C$ and $\hat{C}$ as $p,n\to\infty$ (Silverstein & Bai, 1995). More specifically here, for arbitrary deterministic matrices $M$, Stieltjes transform relations connect the empirical spectral (i.e., eigenvalue) distribution $\nu_{p}$ of $M^{-1}C$ to the empirical spectral distribution $\mu_{p}$ of $M^{-1}\hat{C}$ (Couillet et al., 2018), defined as

[TABLE]

The connecting argument goes as follows: first, from Cauchy’s integral formula (stating that $f(t)=\frac{1}{2\pi\imath}\oint_{\Gamma}f(z)/(t-z)dz$ for $\Gamma$ a complex contour enclosing $t$), the metric $\delta(M,C)$ in (1) relates to the Stieltjes transform $m_{\nu_{p}}(z;M)$ through

[TABLE]

for $\Gamma\subset{\mathbb{C}}$ a (positively oriented) contour surrounding the eigenvalues of $M^{-1}C$. The notation $m_{\nu_{p}}(z;M)$ reminds the dependence of $m_{\nu_{p}}(z)$ in the matrix $M$.

Then, by exploiting the relation between the Stieltjes transforms $\mu_{\nu_{p}}$ and $m_{\mu_{p}}$, it is shown in (Couillet et al., 2018) that, under Assumption 1, for all deterministic $M$ of bounded operator norm,

[TABLE]

almost surely, where $X=[x_{1},\ldots,x_{n}]$ and

[TABLE]

with $G$ such that $G^{\prime}(z)\equiv g(z)=f(1/z)$, $\hat{\Gamma}$ a contour surrounding the support of the almost sure limiting eigenvalue distribution of $M^{-1}\hat{C}$ and $\tilde{\mu}_{p}=\frac{p}{n}\mu_{p}+(1-\frac{p}{n})\delta_{0}$ (and thus $m_{\tilde{\mu}_{p}}(z)=cm_{\mu_{p}}(z)+(1-\frac{p}{n})/z$). Note that, by the linearity of $G$ in (4), it is sufficient in practice to evaluate $\hat{\delta}(M,X)$ for elementary functions (such as $f(z)=z$, $f(z)=\log(z)$, etc.) in order to cover most distances and metrics of interest (see again Table 1). Table 2 reports the values of $G$ for such atomic functions $f$.

Our main idea is to estimate $C$ by minimizing the approximation $\hat{\delta}(M,X)$ of $\delta(M,C)$ over $M$. However, it is important to note that, as discussed in (Couillet et al., 2018), the random quantity $\hat{\delta}(M,X)$ may be negative with non-zero probability. As such, minimizing $\hat{\delta}(M,X)$ over $M$ may lead to negative solutions. Our proposed estimation method therefore consists in approximating $C$ by the solution to the optimization problem

[TABLE]

3 Methodology and Main Results

3.1 Estimation Method

We solve (5) via a gradient descent algorithm in the Riemannian manifold $S_{n}^{++}$ of positive definite $n\times n$ matrices.

To evaluate the gradient $\nabla h_{X}(M)$ of $h_{X}$ at $M\in S_{n}^{++}$, recall that on $S_{n}^{++}$ the differential ${\rm D}h_{X}(M)[\xi]$ of the functional $h_{X}:S_{n}^{++}\to{\mathbb{R}}^{+}$, at position $M\in S_{n}^{++}$ and in the direction of $\xi\in S_{n}$ (the Riemannian manifold of symmetric $n\times n$ matrices), is given by (Absil et al., 2009)

[TABLE]

where $\langle\cdot,\cdot\rangle_{.}^{S_{n}^{++}}$ is the Riemannian metric defined through

[TABLE]

Differentiating $\hat{\delta}^{2}(M,X)$ at $M$ in the direction $\xi$ yields:

[TABLE]

By using the fact that

[TABLE]

where ${\rm sym}(A)=\frac{1}{2}(A+A^{\sf T})$ is the symmetric part of $A\in{\mathbb{R}}^{p\times p}$, we retrieve the gradient of $h_{X}(M)$ as

[TABLE]

(recall that the right-hand side still depends on $X$ implicitly through $\tilde{\mu}_{p}$ and $\hat{C}$).

Once $\nabla h_{X}$ estimated, every gradient descent step in $S_{n}^{++}$ corresponds to a small displacement on the geodesic starting at $M$ and towards $-\nabla h_{X}(M)$, defined as the curve

[TABLE]

where, for $A=U\Lambda U^{\sf T}\in S_{n}^{++}$ in its spectral decomposition, $\exp(A)\equiv U\exp(\Lambda)U^{\sf T}$ (with $\exp$ understood here applied entry-wise on the diagonal elements of $\Lambda$).

That is, letting $M_{0},M_{1},\ldots$ and $t_{0},t_{1},\ldots$ be the successive iterates and step sizes of the gradient descent, we have, for some given initialization $M_{0}\in S_{n}^{++}$,

[TABLE]

Our proposed method is summarized as Algorithm 1.

We conclude this section by an important remark on the fundamental limitations of the proposed algorithm.

Remark 1** (Approximation of $\delta(M_{k},C)$ by $\hat{\delta}(M_{k},X)$).**

It is fundamental to understand the result from (Couillet et al., 2018) at the heart of the proposed method. There, it is precisely shown that, for every deterministic sequence of matrices $\{M^{(p)},~{}p=1,2,\ldots\}$ and $\{C^{(p)},~{}p=1,2,\ldots\}$, with $M^{(p)},C^{(p)}\in{\mathbb{R}}^{p\times p}$ and $\max(\|C^{(p)}\|,\|M^{(p)}\|)<K$ for some constant $K$ independent of $p$, we have that, for $X^{(p)}=[x^{(p)}_{1},\ldots,x^{(p)}_{n}]$ with $x_{i}^{(p)}=C^{(p)\frac{1}{2}}z_{i}^{(p)}$ and $z_{i}^{(p)}$ i.i.d. vectors of i.i.d. zero mean and unit variance entries,

[TABLE]

almost surely as $n,p\to\infty$ and $p/n\to c\in(0,1)$. This result seems to suggest that $\hat{\delta}(M_{k},X)$ in our algorithm is a good approximation for the sought for $\delta(M_{k},C)$. This however only holds true so long that $M_{k}$ is independent of $X$ which clearly does not stand when proceeding to successive gradient descent steps in the direction of $\nabla h_{X}(M)$ which depends explicitly on $X$. As such, while initializations with, say, $M_{0}=I_{p}$, allow for a close approximation of $\delta(M_{k},C)$ in the very first steps of the descent, for larger values of $k$, the descent is likely to drive the optimization in less accurate directions.

Remark 1 is in fact not surprising. Indeed, finding the minimum of $\delta(M,C)$ over $M\succ 0$ would result in finding $C$, which cannot be achieved for unconstrained matrices $C$ and for non vanishing values of $p/n$. Figure 1 provides a typical evolution of the distance $\delta(M_{k},C)$ versus its approximation $\hat{\delta}(M_{k},X)$ at the successive steps $k=1,2,\ldots$ of Algorithm 1, initialized at $M_{0}=I_{p}$. As expected, the difference $|\hat{\delta}(M_{k},X)-\delta(M_{k},C)|$, initially small (at $k=1$, $\hat{\delta}(I_{p},X)\simeq\delta(I_{p},C)$), increases with $k$, until the gradient vanishes and the divergence $\delta(M_{k},C)$ converges.

3.2 Practical Implementation

In order to best capture the essence of Algorithm 1, as well as its various directions of simplification and practical fast implementation, a set of important remarks are in order.

First note that the generic computation of $M_{k+1}$ in Equation (7) may be numerically expensive, unless $M_{k}$ and $\nabla h_{X}(M)$ share the same eigenvectors. In this case, letting $M_{k}=U\Omega_{k}U^{\sf T}$ and $\nabla_{X}h(M_{k})=U\Delta_{k}U^{\sf T}$, we have the recursion

[TABLE]

where $\delta_{k}=\operatorname{diag}(\Delta_{k})$ and $\omega_{k}=\operatorname{diag}(\Omega_{k})$.

In particular, if $M_{0}=\alpha I_{n}+\sqrt{1-\alpha^{2}}\hat{C}$ is a linear shrinkage for some $\alpha\in[0,1]$, we immediately find that, for all $k\geq 0$,

[TABLE]

where $\hat{U}\in{\mathbb{R}}^{p\times p}$ are the eigenvectors of $\hat{C}$ (i.e., in its spectral decomposition, $\hat{C}=\hat{U}\hat{\Lambda}\hat{U}^{\sf T}$) and $\Omega_{k}$ is recursively defined through (8).

This shows that, initialized as such, the ultimate limiting estimate $M_{\infty}$ (i.e., the limit of $M_{k}$) of $C$ shares the same eigenvectors as $\hat{C}$, and thus reduces to a “non-linear shrinkage” procedure, similar to that of (Ledoit & Wolf, 2015). Extensive simulations in fact suggest that, if initialized randomly (say with $M_{0}$ a random Wishart matrix), after a few iterates, the eigenvectors of $M_{k}$ do converge to those of $\hat{C}$ (see further discussions in Section 4). As such, for computational ease, we suggest to initialize the algorithm with $M_{0}=I_{p}$ or with $M_{0}$ a linear shrinkage of $\hat{C}$.

A further direction of simplification of Algorithm 1 relates to the fact that, for generic values of $M_{0}$ (notably having eigenvectors different from those of $\hat{C}$), Equation (7) is computationally expensive to evaluate. A second-order simplification for small $t$ is often used in practice (Jeuris et al., 2012), as follows

[TABLE]

Simulations with this approximation suggest almost no difference in either the number of steps until convergence or accuracy of the solution.

3.3 Estimation of $C^{-1}$

In our framework, estimating $C^{-1}$ rather than $C$ can be performed by minimizing $\delta(M,C^{-1})$ instead of $\delta(M,C)$. In this case, under Assumption 1, (3) now becomes

[TABLE]

almost surely, for every deterministic $M$ of bounded operator norm and $X=[x_{1},\ldots,x_{n}]$, where

[TABLE]

for $F$ such that $F^{\prime}(z)\equiv f(z)$, $\hat{\Gamma}$ a contour surrounding the support of the almost sure limiting eigenvalue distribution of $M\hat{C}$ and $\tilde{\mu}_{p}^{\rm inv}=\frac{p}{n}\mu_{p}^{\rm inv}+(1-\frac{p}{n})\delta_{0}$, where $\mu_{p}^{\rm inv}\equiv\frac{1}{p}\sum_{i=1}^{p}\delta_{\lambda_{i}(M\hat{C})}$. The cost function to minimize under this setting is now given by $h^{\rm inv}(M)\equiv(\hat{\delta}^{\rm inv}(M,X))^{2}$ with gradient $\nabla h_{X}^{\rm inv}(M)$ satisfying

[TABLE]

With these amendments, Algorithm 1 can be adapted to the estimation of $C^{-1}$. Table 2 provides the values of $F$ for the atomic functions $f$ of interest.

3.4 Application to Explicit Metrics

Algorithm 1 is very versatile as it merely consists in a gradient descent method for various metrics $f$ through adaptable definitions of the function $h_{X}(M)=\hat{\delta}(M,X)^{2}$ and its resulting gradient. Yet, because of the integral form assumed by the gradient (Equation (6)), a possibly computationally involved complex integration needs to be numerically performed at each gradient descent step.

In this section, we specify closed-form expressions for the gradient for the atomic $f$ functions of Table 2 (which is enough to cover the list of divergences in Table 1).

3.4.1 Estimation of $C$

Let us denote

[TABLE]

where $V$ are the eigenvectors of $M^{-1}\hat{C}$ and we will determine $\Lambda_{\nabla}$ for each function $f$. Again, we recall from the discussion in Section 3.2, that $V=\hat{U}$ the eigenvectors of $\hat{C}$ if $M$ shares the same eigenvectors as $\hat{C}$ (which thus avoids evaluating the eigenvectors $V$ of $M_{k}^{-1}\hat{C}$ at each step $k$ of the algorithm).

For readability in the following, let us denote $\lambda_{i}\equiv\lambda_{i}(M^{-1}\hat{C})$, $i\in\{1,\ldots,p\}$, the eigenvalues of the matrix $M^{-1}\hat{C}$ and $\xi_{1},\ldots,\xi_{p}$ the eigenvalues of

[TABLE]

with $\Lambda={\rm diag}(\lambda_{1},\ldots,\lambda_{p})$ and $\lambda=(\lambda_{1},\ldots,\lambda_{p})^{\sf T}$. Finally, for $s>0$, let $\kappa_{s}\in(-1/(s(1-p/n)),0)$ be the unique negative number $t$ solution of the equation (see (Couillet et al., 2018) for details)

[TABLE]

With these notations at hand, following the derivations in (Couillet et al., 2018) (detailed in supplementary material), we have the following determinations for $\Lambda_{\nabla}$.

Proposition 1** (Case $f(t)=t$).**

For $f(t)=t$,

[TABLE]

with $m_{\tilde{\mu_{p}}}^{\prime}$ the derivative of $m_{\tilde{\mu_{p}}}$.

Proposition 2** (Case $f(t)=\log(t)$).**

For $f(t)=\log(t)$,

[TABLE]

Proposition 3** (Case $f(t)=\log(1+st)$).**

For $s>0$ and $f(t)=\log(1+st)$,

[TABLE]

Proposition 4** (Case $f(t)=\log^{2}(t)$).**

For $f(t)=\log^{2}(t)$,

[TABLE]

These results are mostly achieved by residue calculus for entire analytic functions $f$ or by exploiting more advanced complex integration methods (in particular branch-cut methods) for more challenging functions (involving logarithms in particular).

Combining these formulas provides an analytical expression for the gradient of all aforementioned divergences and square distances, for the estimation of $C$.

3.4.2 Estimation of $C^{-1}$

Similarly, for the problem of estimating $C^{-1}$, recalling Remark 3.3, we may denote

[TABLE]

with $V_{\rm inv}$ the eigenvectors of $M\hat{C}$.

We redefine in this section $\lambda_{i}\equiv\lambda_{i}(M\hat{C})$, and $\xi_{1},\ldots,\xi_{p}$ the eigenvalues of $\Lambda-\frac{1}{n}\sqrt{\lambda}\sqrt{\lambda}^{\sf T}$ with $\Lambda={\rm diag}(\lambda_{1},\ldots,\lambda_{p})$ and $\lambda=(\lambda_{1},\ldots,\lambda_{p})^{\sf T}$. Again, for $s>0$, let $\kappa_{s}<0$ be the only negative real number $t$ solution of

[TABLE]

With the same approach as in the previous section, we here obtain the following values for $\Lambda^{\rm inv}_{\nabla}$.

Proposition 5** (Case $f(t)=t$).**

For $f(t)=t$,

[TABLE]

Proposition 6** (Case $f(t)=\log(t)$).**

For $f(t)=\log(t)$,

[TABLE]

Proposition 7** (Case $f(t)=\log(1+st)$).**

For $s>0$ and $f(t)=\log(1+st)$,

[TABLE]

Proposition 8** (Case $f(t)=\log^{2}(t)$).**

For $f(t)=\log^{2}(t)$,

[TABLE]

4 Experimental Results

This section introduces experimental results on the direct application of our proposed method to the estimation of $C$ and $C^{-1}$ as well as on its use as a plug-in estimator in more advanced procedures, here in the scope of linear and quadratic discriminant analyses (LDA/QDA).

4.1 Validation on synthetic data

In this first section, we provide a series of simulations on the estimation of $C$ and $C^{-1}$ based on the Fisher distance and for several examples of genuine matrices $C$. Similar results and conclusions were obtained for the other metrics discussed above (the KL divergence and the square Battacharrya distance especially) which are thus not explicitly reported here. The interested reader can refer to the code provided by the authors for self experimentation as supplementary material.

Our preference for the Fisher distance for fair comparisons lies in the fact that it is the “natural” Riemannian distance to compare covariance matrices in $S_{n}^{++}$, therefore in entire agreement with the proposed estimation strategy through gradient descents in $S_{n}^{++}$. Besides, for this specific case, Theorem 4 in (Smith, 2005) establishes an exact and very straightforward formula for the Cramer-Rao bound (CRB) on unbiased estimators of $C$. Although the compared estimators of $C$ are likely all biased and that $M_{0}$ initializations may by chance bring additional information disrupting a formally fair CRB comparison, the CRB at least provides an indicator of relevance of the estimators.

The examples of covariance matrix $C$ under consideration in the following are:

(i) [Wishart] a random ($p$-dimensional) standard Wishart matrix with $2p$ degrees of freedom,

(ii) [Toeplitz $a$] the Toeplitz matrix defined by $C_{ij}=a^{|i-j|}$,

(iii) [Discrete] a matrix $C$ with uniform eigenvector distribution and eigenvalues equal to $.1$, $1$, $3$, $4$ each with multiplicity $p/4$.

Figure 2 (for the estimation of $C$) and Figure 3 (for $C^{-1}$) report comparative performances on the aforementioned $C$ matrices for the SCM, QuEST1, QuEST2, and our proposed estimator, the latter three being initialized at $M_{0}$ the shrinkage estimation from (Ledoit & Wolf, 2004) (consistent with the choice made for QuEST1, QuEST2 in (Ledoit & Wolf, 2015; Ledoit et al., 2018)). In the figures, “SCM th” refers to the asymptotic analytical approximation of $\delta(C,\hat{C})$ as defined in Remark 2. It is observed that, while the SCM never reaches the unbiased CRB, in many cases the QuESTx estimators and our proposed method overtake the CRB, sometimes significantly so. The Wishart matrix case seems more challenging from this perspective. In terms of performances, both our proposed method and QuESTx perform competitively and systematically better than the sample covariance matrix.

Remark 2** (Consistent estimator for $\delta(C,\hat{C})$).**

With the same technical tools from (Couillet et al., 2018), it is straightforward to estimate the distance $\delta(C,\hat{C})$. Indeed, $\delta(C,\hat{C})=\frac{1}{2\pi\imath}\oint_{\Gamma}f(z)m_{\gamma}(z)$ for $m_{\gamma}(z)$ the Stieljes transform of the eigenvalue distribution of $C^{-1}\hat{C}$; the limiting distribution of the latter is the popular Marcenko-Pastur law (Marc̆enko & Pastur, 1967), the expression of which is well known. The estimate is denoted “SCM th” in Figures 2–3. The observed perfect match between limiting theory and practice confirms the consistency of the random matrix approach even for not too large $p,n$.

4.2 Application to LDA/QDA

As pointed out in the introduction, the estimation of the covariance and inverse covariance matrices of random vectors are at the core of a wide range of applications in statistics, machine learning and signal processing. As a basic illustrative example, we focus here on linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA). Both exploit estimates covariance matrices of the data or their inverse in order to perform the classification.

Suppose $x^{(1)}_{1},\ldots,x^{(1)}_{n_{1}}\sim N(\mu_{1},C_{1})$ and $x^{(2)}_{1},\ldots,x^{(2)}_{n_{2}}\sim N(\mu_{2},C_{2})$ are two sets of random independent $p$-dimensional training vectors forming two classes of a Gaussian mixture. The objective of LDA and QDA is to estimate the probability for an arbitrary random vector $x$ to belong to either class by replacing the genuine means $\mu_{a}$ and covariances $C_{a}$ by sample estimates, with in the case of LDA the underlying (possibly erroneous) assumption that $C_{1}=C_{2}$. Defining $C$ as $C\equiv\frac{n_{1}}{n_{1}+n_{2}}C_{1}+\frac{n_{2}}{n_{1}+n_{2}}C_{2}$, the classification rules for LDA and QDA for data point $x$ depend on the signs of the respective quantities:

[TABLE]

where $\hat{\mu}_{a}\equiv\frac{1}{n_{a}}\sum_{i=1}^{n_{a}}x_{i}^{(a)}$ is the sample estimate of $\mu_{a}$ and $\check{C}_{a}^{-1}$ are some estimate of $C_{a}^{-1}$, while $\check{C}\equiv\frac{n_{1}}{n_{1}+n_{2}}\check{C}_{1}+\frac{n_{2}}{n_{1}+n_{2}}\check{C}_{2}$ with $\check{C}_{a}$ the estimation of $C_{a}$. As such, in the following simulations, LDA will exclusively exploit estimations of $C_{1}$ and $C_{2}$ (before inverting their estimated average), while QDA will focus on estimating directly the inverses $C_{1}^{-1}$ and $C_{2}^{-1}$.

The first three displays in Figures 4 and 5 compare the accuracies of the LDA/QDA algorithms for $C_{1}$ and $C_{2}$ chosen among Wishart and Toeplitz matrices, and for $\mu_{2}=\mu_{1}+\frac{80}{p}$ for the LDA and $\mu_{2}=\mu_{1}+\frac{1}{p}$ for the QDA settings (in order to avoid trivial classification). The bottom right displays are applications to real EEG data extracted from the dataset in (Andrzejak et al., 2001). The dataset contains five subsets (denoted A-E). Sets A and B were collected from healthy volunteers while C, D, E were collected from epileptic patients. The graph presents all combinations of binary classes between healthy volunteers and epilectic subjects (e.g., A/E for subsets A and E). There we observe that, for most considered settings, our proposed algorithm almost systematically outperforms competing methods, with QuEST1 and QuEST2 exhibiting a much less stable behavior and particularly weak performances in all synthetic scenarios.

5 Discussion and Concluding Remarks

Based on elementary yet powerful contour integration techniques and random matrix theory, we have proposed in this work a systematic framework for the estimation of covariance and precision matrices. Unlike alternative state-of-the-art techniques that attempt to invert the fundamental Bai–Silverstein equations (Silverstein & Bai, 1995), our proposed method relies on a basic gradient descent approach in $S_{n}^{++}$ that, in addition to performing competitively (if not better), is computationally simpler.

While restricted to metrics depending on the eigenvalues of products of covariance matrices, our approach may be flexibly adapted to further matrix divergences solely depending on eigenvalue relations. The same framework can notably be applied to the Wasserstein distance between zero-mean Gaussian laws. A reservation nonetheless remains on the need for ${n>p}$ in settings involving inverse matrices that, when not met, does not allow to relate the sought-for eigenvalues to the (undefined) inverse sample covariance matrices. (Couillet et al., 2018) shows that this problem can be partially avoided for some divergences (not for the Fisher distance though). A systematic treatment of the ${n<p}$ case is however lacking.

Our approach also suffers from profound theoretical limitations that need be properly addressed: the fact that $\hat{\delta}(M,X)$ only estimates the sought-for $\delta(M,C)$ for $M$ independent of $X$ poses a formal problem when implemented in the gradient descent approach. This needs be tackled: (i) either by estimating the introduced bias so to estimate the loss incurred or, better, (ii) by accounting for the dependence to provide a further estimator $\hat{\hat{\delta}}(M(X),X)$ of $\delta(M(X),C)$ for all $X$-dependent matrices $M(X)$ following a specific form. Notably, given that, along the gradient descent initialized at $M_{0}=\hat{U}D_{0}\hat{U}^{\sf T}$, with $\hat{U}$ the eigenvectors of $\hat{C}$ and $D_{0}$ some diagonal matrix, all subsequent $M_{k}$ matrices share the same profile (i.e., $M_{k}=\hat{U}D_{k}\hat{U}^{\sf T}$ for some diagonal $D_{k}$), a first improvement would consist in estimating consistently $\delta(\hat{U}D\hat{U}^{\sf T},C)$ for deterministic diagonal matrices $D$.

The stability of the eigenvectors when initialized at $M_{0}=\hat{U}D_{0}\hat{U}^{\sf T}$, with $\hat{U}$ the eigenvectors of $\hat{C}$, turns our algorithm into a “non-linear shrinkage” method, as called by (Ledoit & Wolf, 2015). Parallel simulations also suggest that arbitrary initializations $M_{0}$ which do not follow this structure tend to still lead to solutions converging to the eigenspace of $\hat{C}$. However, this might be a consequence of the inconsistency of $\hat{\delta}(M,X)$ for $X$-dependent matrices $M$: it is expected that more consistent estimators might avoid this problem of “eigenvector of $\hat{C}$” attraction, thereby likely leading to improved estimations of the eigenvectors of $C$.

We conclude by emphasizing that modern large dimensional statistics have lately realized that substituting large covariance matrices by their sample estimators (or even by improved covariance estimators) is in general a weak approach, and that one should rather focus on estimating some ultimate functional (e.g., the result of a statistical test) involving the covariance (see, e.g., (Mestre & Lagunas, 2008) in array processing or (Yang et al., 2015) in statistical finance). It is to be noted that our proposed approach is consistent with these considerations as various functionals of $C$ can be obtained from Equation (2), from which similar derivations can be performed.

Acknowledgement

This work is supported by the ANR Project RMT4GRAPH (ANR-14-CE28-0006) and by the IDEX GSTATS Chair at University Grenoble Alpes.