Benoit Oriol
email: benoit.oriol@dauphine.eu

Alexandre Miot

Benoît Oriol

Ledoit-Wolf linear shrinkage with unknown mean

Keywords: covariance matrix estimation, linear shrinkage, Ledoit-Wolf estimator, unknown mean, general asymptotics

come

Introduction and related work

The covariance matrix plays a major role in numerous machine learning algorithms and statistics. Just to cite a few, the PCA [START_REF] Pearson | on lines and planes of closest fit to systems of points in space[END_REF] in machine learning, Markowitz portfolio management [START_REF] Markowitz | Portfolio selection[END_REF] in finance, or generalized method of moments estimators [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF] in statistics. However, those algorithms are designed to use the true covariance matrix, which is often unaccessible. Even if the sample covariance matrix seems to be a simple and appealing choice, it severely fails in many applications: for instance, the use of the sample covariance matrix for Markowitz portfolio management does not beat a naive uniform distribution among the assets [START_REF] Demiguel | Optimal versus naive diversification: How inefficient is the 1/n portfolio strategy[END_REF].

In the context of Kolmogorov asymptotics, where the ratio of the dimension p n and the number of samples n tends to a finite positive constant pn n → c > 0, this estimator fails to converge quadratically. Moreover, its eigenvalue spectrum is biased: high eigenvalues tends at 2 Notations, definitions and hypotheses Let us introduce the following notations. Notation 1. In the following we consider a sequence of observation matrices (X n ) n∈N * with X n ∈ R pn×n of n iid observations on a system of p n dimensions. Decomposing the covariance matrix, we denote Σ n = Γ n Λ n Γ t n , where Λ n is a diagonal matrix and Γ n a rotation matrix. The diagonal elements of Λ n are the eigenvalues λ n 1 , ..., λ n pn , and the columns of Γ n are the eigenvectors γ n 1 , ..., γ n pn . Y n = Γ t n X n is a p n × n matrix of n iid observations of p n uncorrelated random variables (y n 1 , ..., y n n ).

Notation 2. Let A n and B n two p n × p n matrices. We consider the Frobenius norm: A n = tr(A n A t n )/p n , and the associated inner product: A n , B n n = tr(A n B t n )/p n . Dividing by the dimension is not standard, it is done to fix the norm of the identity as 1 regardless of the dimension.

Notation 3. Let (E n ) n a sequence of euclidean spaces with associated norm • n . The quadratic convergence of a random variable Z n ∈ E n , i.e. E[ Z n

2 n ] → 0, is denoted as Z n -→ q.m 0.

We describe now several assumptions, the same used in the linear shrinkage of Ledoit and Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], that will be used in the following. Assumption 1. There exists a constant K 1 independent of n such that p n /n ≤ K 1 .

Assumption 2. There exists a constant K 2 independent of n such that

ỹn l1 ]) 2 |Q n | = 0,
where Q n denotes the set of all the quadruples that are made of four distinct integers between 1 and p n , and for all i ∈ 1, p n , ỹn

i1 = y n i1 -E[y n i1 ].
We need some definitions to properly define the problem and the asymptotics.

Definition 1 (Empirical covariance). For an observation matrix X n of size p n × n, we define the empirical covariance as:

S n = Xn Xt n /(n -1), with ( Xn ) ik = (X n ) ik -1 n n k =1 (X n ) ik . Definition 2 (Scalars (µ n , α 2 n , β 2 n , δ 2 n )).
We define four scalars:

µ n = Σ n , I pn n , α 2 n = Σ n -µ n I pn 2 n , β 2 n = E[ S n -Σ n 2 n ], δ 2 n = E[ S n -µ n I pn 2 n ].
(Lemma 2.1 in [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]) proves that α 2 n + β 2 n = δ 2 n .

The oracle linear shrinkage estimator is given by the following minimization problem. The following corollary is the central point of the linear shrinkage methods.

Corollary 1 (Corollary of theorem 2.1 from [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]). Consider the optimization problem:

minimize ρ1,ρ2 s.t. Σ * n =ρ1Ip n +ρ2Sn E[ Σ * n -Σ n 2 n ],
where the coefficients ρ 1 and ρ 2 are not random. Its solution Σ * n verifies:

Σ * n = β 2 n δ 2 n µ n I pn + α 2 n δ 2 n S n , E[ Σ * n -Σ n 2 n ] = α 2 n β 2 n δ 2 n .
Remark 1. Corollary 1 remains true for any unbiased estimator Σn instead of S n .

(µ n , α 2 n , β 2 n , δ 2 n ) depends on the true covariance Σ n , and thus can't be used directly in the estimation of Σ * n . The central issue of this work is to find estimators (m n , a 2 n , b 2 n , d 2 n ) of (µ n , α 2 n , β 2 n , δ 2 n ) in order to compute an estimation S * n of Σ * n . As the mean is unknown, those estimators differ from Ledoit and Wolf work [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], particularly when p n is higher than n.

Theoretical results

All the following results extend the work of Ledoit and Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] in the case where the empirical mean is used as estimator of the mean. All proofs are shown in appendix A.

Remark 2. In the following, as all the estimators are invariant by change of mean, resulting from the definition of X = Xk X •,k , we can assume E[X] = 0 for the simplicity of notations.

We present a sequence of lemmata, that naturally define estimators with suitable asymptotic properties for the scalars (µ n , α 2 n , β 2 n , δ 2 n ).

Lemma 1. Under assumptions 1 and 2, µ n , α 2 n , β 2 n , δ 2 n remain bounded as n → ∞.

Theorem 1. Under assumptions 1 and 2, define

θ 2 n = V 1 pn pn i=1 (y n i1 ) 2 ]
. θ 2 n is bounded as n → ∞, and we have:

lim n→∞ E[ S n -Σ n 2 n ] - p n n (µ 2 n + θ 2 n ) = 0.
In particular, taking p n = p constant, we see that S n -→ q.m Σ n . But, when p n is of the same order of magnitude than n, the sample covariance generally fails to converge as the error is at least of the same order of magnitude as µ 2 n = Σ n , I pn

2 n .
Lemma 2 (Estimator of µ n ). Define m n = S n , I pn n . Then, under assumptions 1 and 2, E[m n ] = µ n for all n, and m n -µ n converges to zero in quartic mean (fourth moment) as n goes to infinity.

Corollary 2. Under assumptions 1 and 2, m 2 n -E[m 2 n ] converges to zero in quadratic mean as n goes to infinity.

Lemma 3 (Estimator of δ 2 n ). Define d 2 n = S n -m n I pn 2
n . Then, under assumptions 1, 2 and 3,

d 2 n -δ 2 n -→ q.m 0. It follows that d 2 n -E[d 2 n ] -→ q.m 0.
The following lemmata aim at defining an unbiased estimator of β 2 n that quadratically converges. We work around b2 n , inspired from the estimator of β 2 n in the case where the mean is known, and write residual terms in the expectation as a combination of m n and d 2 n .

Lemma 4. Define:

b2 n = 1 n 2 n k=1 n n -1 xn •k (x n •k ) t -S n 2 n ,
where xn

•k = x n •k -1 n n k =1 x n •k and (x n •1 , ..., x n •k , ..., x n •n ) are the independent samples forming X. Then, under assumption 1, E[ b2 n ] = c 0 β 2 n + c 1 δ 2 n + c 2 µ 2 n , with γ n = n(n-1) n 2 -3n+3 , λ n = n 2 (n-2) (n-1)(n 2 -3n+3) , c 0 = 1 γn -1 n -λn γnn 2 , c 1 = λn γnn 2 , c 2 = (p + 1)c 1 .
Lemma 5. Under assumption 1, we have:

E[ b2 n ] = c 0 β 2 n + c 1 E[d 2 n ] + c 2 E[m 2 n ] + (c 1 -c 2 )V[m n ]. Lemma 6.
Under assumption 1, we have:

V[ b2 n ] -→ n→∞ 0.
We need to compute V[m n ] which is unknown in the development of E[ b2 n ] in Lemma 5.

Lemma 7.

Under assumption 1, we have:

V[m n ] = q 0 β 2 n + q 1 δ 2 n -q 2 µ 2 n , with q 0 = n-2 p(n-1) , q 1 = 1 p(n-1) , q 2 = p-1 p(n-1) .
Lemma 8. Under assumption 1, we have:

V[m n ] = 1 1-q1-q2 q 0 β 2 n + q 1 E[d 2 n ] -q 2 E[m 2 n ] . Lemma 9. Define: b 2 n = 1 c f 0 b2 n -c f 1 d 2 n -c f 2 m 2 n , with c f 0 = c 0 + (c 1 -c 2 ) q0 1-q1-q2 , c f 1 = c 1 + (c 1 -c 2 ) q1 1-q1-q2 , c f 2 = c 2 -(c 1 -c 2 ) q2 1-q1-q2 . Then, under assumptions 1, 2 and 3, b 2 n is an unbiased estimator of β 2 n , i.e. E[b 2 n ] = β 2 n , and b 2 n -β 2 n -→ q.m 0.
For notation consistency with the estimators in Ledoit-Wolf linear shrinkage [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], we keep the notation b 2 n even if its value can be negative. Lemma 10

(Estimator of β 2 n ). Define: b 2 n,u = min((b 2 n ) + , d 2 n ) and a 2 n,u = d 2 n -b 2 n,u . Under assumptions 1, 2 and 3, b 2 n,u -β 2 n -→ q.m 0 and a 2 n,u -α 2 n -→ q.m 0.
We can now define our linear shrinkage estimator S * n and prove its asymptotic properties. Definition 3 (Final estimator LW u). Let's define our estimator:

S * n = b 2 n,u d 2 n,u m n,u I pn + a 2 n,u d 2 n,u S n , with m n,u = m n , d 2 n,u = d 2 n .
Theorem 2. Under assumptions 1, 2 and 3,

E[ S * n -Σ * n 2 ] → 0. As a consequence, S * n has the same asymptotic expected loss as Σ * n , i.e. E[ S * n -Σ n 2 n ] -E[ Σ * n -Σ n 2
n ] → 0. The following lemma gives an asymptotic estimation of the optimal error

E Σ * n -Σ n 2 n = α 2 n β 2 n δ 2 n . Lemma 11. E a 2 n,u b 2 n,u d 2 n - α 2 n β 2 n δ 2 n → 0.
The last results easily make possible to extend the Theorems 3.3 and 3.4 of ( [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] in our situation where the mean is unknown. Previously, we showed that our estimator's loss converge to the optimal one in the class of linear combinations of S n and I pn with non random coefficients, the optimal estimator of this class being Σ * n . In the following, we show that our estimator is still asymptotically optimal with respect to a bigger class, where the coefficients can be random. Formally, we are looking for the following optimal loss (this time, there is no expectation in the minimization). Let Σ * * n be the linear combination of S n and I pn solving:

minimize ρ1,ρ2 s.t. Σ * * n =ρ1Ip n +ρ2Sn Σ * * n -Σ n 2 n .
By construction, Σ * * n has a lower loss than S * n , but we show that the difference converges to 0. Theorem 3. S * n converges to Σ * * n in quadratic mean, i.e. S * n -Σ * * n -→ q.m 0. As a consequence, S * n has the same asymptotic expected loss as Σ * * n , more precisely we have:

E S * n -Σ n 2 n -Σ * * n -Σ n 2 n → 0.
Theorem 4. For any sequence of linear combinations Σn of I n and S n , the estimator S * n verifies:

lim N →∞ inf n≥N E Σn -Σ n 2 n -E S * n -Σ n 2 n ≥ 0.
In addition, every Σn that performs as well as S * n is identical to S * n in the limit:

lim N →∞ E Σn -Σ n 2 n -E S * n -Σ n 2 n = 0 ⇐⇒ Σn -S * n n -→ q.m 0.
We introduce three other estimators to compare with, which are implemented, recommended, or natural to define. We prove that their asymptotic behavior is similar, and, through different experiments, show the differences in performance.

Definition 4 (Ledoit-Wolf recommended estimators). The estimators recommended by Ledoit and Wolf [START_REF] Ledoit | The power of (non-)linear shrinking: A review and guide to covariance matrix estimation[END_REF], indexed by the letter "r", are:

m n,r = m n , d 2 n,r = d 2 n , b 2 n,r = min   1 (n -1) 2 n k=1 xn •k (x n •k ) t -S n 2 n + , d 2 n,r   , a 2 n,r = d 2 n,r -b 2 n,r .
Theorem 5 (Ledoit-Wolf recommended estimators). Under Assumptions 1, 2 and 3, m n,r -µ n converges to 0 in quartic mean, and that

d 2 n,r -δ 2 n , b 2 n,r -β 2
n and a 2 n,r -α 2 n converge in quadratic mean to 0 as n goes to infinity.

Moreover, the conclusions of Theorem 2 remain true with the estimated matrix

S * n,r = b 2 n,r d 2 n,r m n,r I pn + a 2 n,r d 2 n,r S, i.e. E[ S * n,r -Σ * n 2 ] → 0 and E S * n,r -Σ n 2 n -Σ * * n -Σ n 2 n → 0. From the proof, b2 n,r = b2 n + 1 n(n-1) 2 S n 2
n , it is then natural to define the following estimator.

Definition 5 ("Natural" estimators). The estimators that naturally emerge, indexed by the letter "m", are:

m n,m = m n , d 2 n,m = d 2 n , b 2 n,m = min b2 n + , d 2 n,m , a 2 n,m = d 2 n,m -b 2 n,m .
Theorem 6 ("Natural" estimators). Under Assumptions

m n,s = n -1 n m n , d 2 n,s = d 2 n , b 2 n,s = min b2 n + , d 2 n,s , a 2 n,s = n -1 n d 2 n,s -b 2 n,s , S * n,s = n -1 n S *

Experimental results

The experimental estimations are compared to the theoretical value of Σ * in the Ledoit-Wolf setting, the implementation in ScikitLearn 1.2.2, the implementation recommended by Ledoit and Wolf [START_REF] Ledoit | The power of (non-)linear shrinking: A review and guide to covariance matrix estimation[END_REF], and to the other algorithms implemented in ScikitLearn 1.2.2, for multivariate Gaussian and Student-t distributions.

We first derive the exact values of Σ * for those two distributions.

Oracle estimators

4.1.1 Gaussian distribution Lemma 12. Let (X n ) •,k ∼ N (0, Σ n ), k ∈ 1, n , n iid samples.
Then, the analytical oracle estimators are:

µ n = Σ n , I n , α 2 n = Σ n 2 n -µ 2 n , β 2 n = p+1 n-1 µ 2 n + 1 n-1 α 2 n , δ 2 n = α 2 n + β 2 n . 4.1.2 t-distribution oracle estimator Lemma 13. Let (X n ) •,k ∼ t ν (0, Σn ), k ∈ 1, n , n iid samples with scale matrix Σn = ν-2 ν Σ n and covariance V[X n ] = Σ n .
The density of the multivariate t-distribution is:

f n (x) = Γ[(ν + p)/2] Γ(ν + p)ν p/2 π p/2 | Σn | 1/2 1 + 1 ν x T Σ-1 n x -(ν+p)/2
.

Then, the analytical oracle estimators are:

µ n = Σ n , I , α 2 n = Σ n 2 n -µ 2 n , β 2 n = 1 n ν ν-4 + 1 n-1 α 2 n + (p + 1)µ 2 n -2p n(ν-4) µ 2 n , δ 2 n = α 2 n +β 2 n .

Experimental setup

We considered 2 settings:

• a Monte-Carlo computation of the loss on a 2d-grid of the parameters (p n , n) ∈ 5, 100 2 , with a step size of 2, to visualize the effect of changing the ratio p n /n, and see the domains where our algorithm is most suited;

• a Monte-Carlo computation of the loss with a fixed ratio p n /n = c ∈ {1, 2, 4}, to compare the rate of convergence of each algorithm.

In both cases, the Monte-Carlo is computed with n M C = 10000 iterations. Three different distributions are explored: the multivariate Gaussian, and the t-distribution with ν = 10 and ν = 8.5. Note that we have to ensure ν > 8 to respect Assumption 2. Two different way of choosing Σ n are explored: fixing Σ n = I pn -particular case where the oracle Ledoit-Wolf loss is null -, and drawing at each iteration a covariance matrix Σ n from a Wishart distribution with p n degrees of freedom, and normalizing it by Σ n Σ T n n -to respect the assumption 2. Note that when drawn from a Wishart, Σ n Σ T n n > 0 almost surely.

Assumptions check

For the second study at fixed p n /n in order to compare the rate of convergence, we check that we are under the three assumptions that guarantee the theoretical results on convergence proved in section 3.

Assumption 1 As we fixed the ratio c = p n /n, Assumption 1 is trivially respected.

Assumption 2 -Gaussian distribution Let (X n ) •,k ∼ N (0, Σ n ), k ∈ 1, n , n iid samples. As previously, we denote Σ n = Γ n Λ n Γ t n
, where Λ n is a diagonal matrix and Γ n a rotation matrix, and Y n = Γ t n X n is a p n × n matrix of n iid observations of p n uncorrelated random variables. Using the fact that for z ∼ N (0, λ), we have E[z 8 ] = 105λ 4 , we deduce:

1 p n pn i=1 E[y 8 i1 ] = 1 p n pn i=1 105λ 4 ii = 105 Λ n Λ T n 2 n = 105 Σ n Σ T n 2 n . (1)
In the case where we fix Σ n = I pn , we obviously have Σ n Σ T n 2 n = 1, so Assumption 2 is respected.

In the case where we draw Σ n from a Wishart distribution with p n degrees of freedom, and normalize it by Σ n Σ T n n , we have by construction Σ n Σ T n 2 n = 1, so Assumption 2 is respected here too.

Assumption 2 -t-distribution Let (X n ) •,k ∼ t ν (0, Σn ), k ∈ 1, n , n iid samples with ν > 8, scale matrix Σn = ν-2 ν Σ n and covariance V[X n ] = Σ n . As previously, we denote Σ n = Γ n Λ n Γ t n
, where Λ n is a diagonal matrix and Γ n a rotation matrix, and Y n = Γ t n X n is a p n × n matrix of n iid observations of p n uncorrelated random variables. From a characterization of multivariate t-distributions, for each k ∈ 1, n , there exist 2 independent random variables U k and Z •,k such that:

U k ∼ χ 2 ν , Z •,k ∼ N 0, ν -2 ν Λ n , (Y n ) •,k = Σ -1/2 (X n ) •,k = ν U k Z •,k .
Moreover, we notice that:

E 1 U 4 1 = R (1/2) ν/2-1 Γ(ν/2) x (ν-8)/2-1 e -x/2 dx = (1/2) ν/2 (1/2) (ν-8)/2 Γ((ν -8)/2) Γ(ν/2) R (1/2) (ν-8)/2-1 Γ((ν -8)/2) x (ν-8)/2-1 e -x/2 pdf of χ 2 ν-8 dx E 1 U 4 1 = 1 (ν -8)(ν -6)(ν -4)(ν -2) . (2)
This 2 previous points lead to:

1 p n pn i=1 E[y 8 i1 ] = 1 p n pn i=1 E[z 8 i1 ]E ν 4 U 4 1 = 1 p n pn i=1 105 ν -2 ν λ ii 4 ν 4 (ν -8)(ν -6)(ν -4)(ν -2) = 105 (ν -8)(ν -6)(ν -4) (ν -2) 3 ΛΛ T 2 n 1 p n pn i=1 E[y 8 i1 ] = 105 (ν -8)(ν -6)(ν -4) (ν -2) 3 ΣΣ T 2 n . (3)
Similarly as the Gaussian case, when we fix Σ n = I pn , we obviously have Σ n Σ T n 2 n = 1, so Assumption 2 is respected, and when we draw Σ n from a Wishart distribution with p n degrees of freedom, and normalize it by Σ n Σ T n n , we have by construction Σ n Σ T n 2 n = 1, so Assumption 2 is respected here too.

Assumption 3 Let Y a p n -dimensional random variables drawn from a centered multivariate Gaussian distribution, then for all (i, j) ∈ 1, p n 2 , we have the following property:

Cov(Y i , Y j ) = 0 =⇒ Y i , Y j independent.
Moreover, when drawing X n from n iid multivariate Gaussian, we have that Y n , using the previous notations, is made of n iid samples of an uncorrelated p n -dimensional centered multivariate Gaussian distribution. So, for all (i, j) ∈ 1, p n 2 , i = j, we have that y i1 , y j1 independent. Finally, for all (i, j, k, l) ∈ 1, p n 4 where i, j, k, l are all different, we have:

Cov[y n i1 y n j1 , y n k1 y n l1 ] = 0.

In the case where X n are n p n -dimensional iid samples drawn from a centered multivariate t-distribution, we have that Y n , using the previous notations, is made of n iid samples of an uncorrelated p n -dimensional centered multivariate t-distribution. Then we use the decomposition Y = ν U Z where Z is drawn from a multivariate Gaussian distribution independent from U , drawn from a χ 2 ν distribution. As for all i = j, Cov(y n i1 , y n j1 ) = 0, then we trivially have Cov(z n i1 , z n j1 ) = 0. So z i1 and z j1 are independent, which immediately leads to the fact that for all (i, j, k, l) ∈ 1, p n 4 where i, j, k, l are all different, we have:

Cov[y n i1 y n j1 , y n k1 y n l1 ] = 0.
This proves that Assumption 3 is respected in all the experimental cases we studied.

Results

In the following, we will use abbreviations to refer the different expected losses of each algorithms.

Concerning the variants of Ledoit-Wolf shrinkage estimators with unknown mean, we denote:

• LW u for the estimator we propose in this paper,

• LW r for the implementation recommended by Ledoit and Wolf in 2020 [START_REF] Ledoit | The power of (non-)linear shrinking: A review and guide to covariance matrix estimation[END_REF],

• LW s for the implementation of ScikitLearn 1.2.2,

• LW m for the natural estimator,

• LW ex for the oracle estimator Σ * ,

• LW op for the optimal estimator Σ * * , Concerning the other baseline algorithms implemented in ScikitLearn, we have:

• EC for the Empirical Covariance estimator,

• SC for the Shrunk Covariance estimator,

• OAS for the Oracle Approximated Shrinkage estimator,

We didn't run the Elliptic Envelope, GLasso and MinCovDet estimators present in ScikitLearn, due to time complexity: in our setup, the computing time of those ones exceeds by a factor at least 10 the computing time of the shrinkage estimators listed before. Consequently, for reason of feasibility, we chose not to compare to them, considering that the latter algorithms are part of a different class of estimators.

Constant covariance Σ n = I pn

Study on a grid over (p,n)

As they often show similar behaviors, we only show a subset of the experimental results for brevity. The three estimators LW s, LW r and LW m have a very similar behavior compared to LW u in this scenario, that's why we will only show the comparison with LW m, having the best performance among the three. The results are shown in figure 1. The black contour on the surface plots is the iso-line at level 0, where the expected losses are equal. In this scenario, LW u is constantly better than the other estimators, and the important difference is in the part p > n, where the mean estimation affects a lot the overall covariance estimation.

Convergence study

We now fix c = p n /n and study the convergence of the different algorithms we cited in the experimental setup. We only show the results with c = 1 as the other cases only widen the differences but do not change the order. The cases t 10 and t 8.5 -distributions are very similar, that's why we show only the t 8.5 one. The key difference between the Gaussian case and the t-distribution, is that the OAS doesn't converge in the latter, while not being so efficient in the Gaussian case which is tailored for it. The results are shown in figure 2.

Covariance drawn from normalized Wishart

Study on a grid over (p,n)

The three estimators LW r and LW m have a very similar behavior compared to LW u in this scenario, that's why we will only show the comparison with LW m, having the best performance among the two. Moreover, the results between the t 10 and t 8.5 -distributions are very similar, so we will only show the t 8.5 case. The results are shown in figures 3 and 4. The black contour on the surface plots is the iso-line at level 0, here it is where the losses are equal. In the case p > n, LW u is far better than the other estimators, where the mean estimation affects a lot the overall covariance estimation. In a finite subset of the part n > p, LW s is slightly better. LW m presents no significant advantage compared to the two others. Here, note that LW ex and LW op have both null loss, we didn't plot them.

Convergence study

We fix c = p n /n and study the convergence of the different algorithms we cited in the experimental setup. The cases t 10 and t 8.5 -distributions being very similar, we only show t 8.5 . The key difference between the Gaussian case and the t-distribution, is that the OAS doesn't converge in the latter, while not being so efficient in the Gaussian case which is tailored for it. The results are shown in figures 5 and 6.

Conclusion

In this work, we extended the linear shrinkage approach of Ledoit and Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] for covariance matrix estimation to the case where the mean of the distribution is unknown. Theoretically, we showed that in this case we have similar asymptotic properties as in the situation when the mean is known. Four different estimators emerged, three around those implemented in Scik-itLearn or recommended by Ledoit and Wolf, and one naturally emerging from the theoretical proofs. Experimentally, the latter showed improved performances in a large spectrum of situations compared to ScikitLearn 1.2.2 baselines and to the three other estimators presented in the theoretical part. The gain in performance is particularly high when the dimension is bigger than the number of samples, while the differences are comparably low when dimension is smaller than the number of samples. Ledoit and Wolf developed several non-linear shrinkage estimators where the mean is known [START_REF] Ledoit | Nonlinear shrinkage estimation of large-dimensional covariance matrices[END_REF][START_REF] Ledoit | Optimal estimation of a large-dimensional covariance matrix under Stein's loss[END_REF][START_REF] Ledoit | Analytical nonlinear shrinkage of large-dimensional covariance matrices[END_REF]. Work needs to be conducted to investigate if similar approach can be used to extend their non-linear frameworks.

Appendix A: Proofs of the technical results

For brevity, we omit the subscript n; but it is understood that everything depends on n. Coefficients of Λ are denoted λ ij , and if not stated otherwise, sum indices i, j, l, m are in 1, p and k 1 , k 1 , k 2 , k 2 , ... are in 1, n . Moreover, we denote ȳ = 1 n k y •,k . We recall that, from Remark 2, as all the estimators are invariant by change of mean from the definition of X = Xk X •,k , we assume E[X] = 0 for the simplicity of notations.

Technical lemma

Frequently used identities are proven in preamble of the other proofs here.

Identity 1

Let (i, j) ∈ 1, p 2 , then:

k (y ik -ȳi )(y jk -ȳj ) = k y ik - 1 n k y ik y jk - 1 n k y jk = 1 n 2 k,k ,k (y ik -y ik )(y jk -y jk ) = 1 n 2 k,k ,k (y ik y jk -y ik y jk -y ik y jk + y ik y jk ) = 1 n 2 k,k ,k (y ik y jk -y ik y jk -y ik y jk + y ik y jk ) (re-indexing) = 1 n 2 k,k ,k y ik (y jk -y jk ) k (y ik -ȳi )(y jk -ȳj ) = 1 n k,k =k y ik (y jk -y jk ).
(4)

6.1.2 Identity 2 1 p i,j λ 2 ij = Λ 2 p = Σ 2 p = Σ -µI 2 p + µI 2 p = α 2 + µ 2 .
(5)

6.1.3 Identity 3

1 p i,j λ ii λ jj = 1 p i λ ii 2 = 1 p tr(Λ) 2 = 1 p tr(Σ) 2 = p Σ, I 2 p = pµ 2 . ( 6
)

Proof of Lemma 1

We have:

Σ 2 = Λ 2 = 1 p p i=1 E[y 2 i1 ] 2 ≤ 1 p p i=1 E[y 4 i1 ] ≤ 1 p p i=1 E[y 8 i1 ] ≤ K 2 . (7)
As µ = Σ, I ≤ Σ , µ remains bounded as n goes to infinity. Also, α 2 = Σ -µI 2 = Σ 2 -µ 2 , so remains bounded as n goes to infinity too. For β 2 , we will deeply decompose the expectation. This is not absolutely necessary to prove the boundedness, but the decomposition will be of utter importance in the following proofs. So, we have:

β 2 = E[ S -Σ 2 ] = 1 p i,j E   1 n -1 k (y ik -ȳi )(y jk -ȳj ) -λ ij 2   = 1 p i,j E      1 n(n -1) k,k =k y ik (y jk -y jk ) -λ ij   2    (Identity 1) β 2 = 1 p i,j 1 n 2 (n -1) 2 k1 k 1 =k1 k2 k 2 =k2 E y ik1 (y jk1 -y jk 1 ) -λ ij y ik2 (y jk2 -y jk 2 ) -λ ij . (8) We denote, for k 1 , k 1 = k 1 , k 2 , k 2 = k 2 : E ij (k 1 , k 1 , k 2 , k 2 ) = E y ik1 (y jk1 -y jk 1 ) -λ ij y ik2 (y jk2 -y jk 2 ) -λ ij . (9)
• If |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 0: E ij (k 1 , k 1 , k 2 , k 2 ) = 0. • If |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 1: -If k 1 = k 2 : then k 1 = k 2 , and k 1 = k 2 . So, E ij (k 1 , k 1 , k 1 , k 2 ) = E y ik1 (y jk1 -y jk 1 ) -λ ij y ik1 (y jk1 -y jk 2 ) -λ ij = E (y ik1 y jk1 -λ ij ) 2 = V[y ik1 y jk1 ] E ij (k 1 , k 1 , k 1 , k 2 ) = V[y i1 y j1 ]. (10)
Moreover, the number of terms in the initial sum on

k 1 , k 1 = k 1 , k 2 , k 2 = k 2 satisfying
the conditions of this case on the indices is:

(k 1 , k 1 , k 2 , k 2 ) ∈ 1, n 4 |k 1 = k 1 , k 2 = k 2 , |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 1, k 1 = k 2 = n(n -1)(n -2). ( 11
)
-If k 1 = k 2 (or similarly k 1 = k 2 ): then k 1 = k 2 and k 1 = k 2 . So, E ij (k 1 , k 1 , k 2 , k 1 ) = 0. -If k 1 = k 2 : E ij (k 1 , k 1 , k 2 , k 1 ) = 0. • If |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 2: -If k 1 = k 2 and k 1 = k 2 : then k 1 = k 1 . So, E ij (k 1 , k 1 , k 1 , k 1 ) = E y ik1 (y jk1 -y jk 1 ) -λ ij 2 = V[y ik1 y jk1 ] + λ ii λ jj E ij (k 1 , k 1 , k 1 , k 1 ) = V[y i1 y j1 ] + λ ii λ jj . (12)
Moreover, the number of terms in the initial sum on

k 1 , k 1 = k 1 , k 2 , k 2 = k 2 satisfying
the conditions of this case on the indices is:

(k 1 , k 1 , k 2 , k 2 ) ∈ 1, n 4 |k 1 = k 1 , k 2 = k 2 , |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 2, k 1 = k 2 , k 1 = k 2 = n(n -1). ( 13
)
-If k 1 = k 2 and k 1 = k 2 : then k 1 = k 1 . So, E ij (k 1 , k 1 , k 1 , k 1 ) = E y ik1 (y jk1 -y jk 1 ) -λ ij y ik 1 (y jk 1 -y jk1 ) -λ ij E ij (k 1 , k 1 , k 1 , k 1 ) = λ 2 ij . (14)
Moreover, the number of terms in the initial sum on

k 1 , k 1 = k 1 , k 2 , k 2 = k 2 satisfying
the conditions of this case on the indices is:

(k 1 , k 1 , k 2 , k 2 ) ∈ 1, n 4 |k 1 = k 1 , k 2 = k 2 , |{k 1 , k 1 } ∩ {k 2 , k 2 }| = 2, k 1 = k 2 , k 1 = k 2 = n(n -1). ( 15
)
Using the latter decomposition on (k 1 , k 1 , k 2 , k 2 ), we deduce:

β 2 = 1 p i,j 1 n 2 (n -1) 2 n(n -1)(n -2)V[y i1 y j1 ] + n(n -1)(V[y i1 y j1 ] + λ ii λ jj ) + n(n -1)λ 2 ij β 2 = 1 pn i,j V[y i1 y j1 ] + p + 1 n(n -1) µ 2 + 1 n(n -1) α 2 (Identities 2 & 3). ( 16
)
• µ 2 is bounded when n goes to infinity, so p+1 n(n-1) µ 2 ≤ K1+1 n-1 µ 2 remains bounded too.

• α 2 is bounded when n goes to infinity, so 1 n(n-1) α 2 remains bounded too.

• Finally, following [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [10] proof of Lemma 3.1,

1 pn i,j V[y i1 y j1 ] ≤ 1 pn i,j E[y 2 i1 y 2 j1 ] ≤ 1 pn i,j E[y 4 i1 ]E[y 4 j1 ] (Cauchy-Schwarz) ≤ 1 pn i,j 4 E[y 8 i1 ]E[y 8 j1 ] (Cauchy-Schwarz) = p n i 1 p 4 E[y 8 i1 ] 2 ≤ p n i 1 p E[y 8 i1 ] (Jensen for 4 √ • concave) 1 pn i,j V[y i1 y j1 ] ≤ K 1 K 2 . ( 17
)
So, β 2 remains bounded as n goes to infinity. Finally, δ 2 = α 2 + β 2 is also bounded as n goes to infinity, which conclude the proof of the lemma.

Proof of Theorem 1

We have:

µ 2 + θ 2 = E 1 p i y 2 i1 2 + V 1 p i y 2 i1 = E   1 p i y 2 i1 2   = 1 p 2 ij E y 2 i1 y 2 j1 = 1 p 2 ij V [y i1 y j1 ] + λ 2 ij = 1 p (α 2 + µ 2 ) + 1 p 2 ij V [y i1 y j1 ] ( Identity 2)
µ 2 + θ 2 ≤ 1 p (α 2 + µ 2 ) + K 2 (Proof of Lemma 1). (18)
As α 2 and µ 2 remains bounded as n goes to infinity, so is θ 2 . And, from the proof of Lemma 1:

β 2 = 1 pn i,j V[y i1 y j1 ] + p + 1 n(n -1) µ 2 + 1 n(n -1) α 2 . ( 19
)
So,

β 2 = p n µ 2 + θ 2 - 1 p (α 2 + µ 2 ) + p + 1 n(n -1) µ 2 + 1 n(n -1) α 2 = p n µ 2 + θ 2 - 1 n(n -1) α 2 + p -n + 2 n(n -1) µ 2 . ( 20
)
As α 2 , µ 2 and p n remains bounded as n goes to infinity, immediately we have:

β 2 - p n µ 2 + θ 2 -→ n→∞ 0. ( 21
)

Proof of Lemma 2

By linearity of the inner product, we trivially have:

E[m] = µ.
For the quartic mean convergence, we write:

E[(m -µ) 4 ] = E   1 p i 1 n -1 k (y ik -ȳi ) 2 - n -1 n λ ii 4   = E   1 n -1 k 1 p i (y ik -ȳi ) 2 - n -1 n λ ii 4   . (22)
Using Identity 1, we obtain:

E[(m -µ) 4 ] = 1 (n -1) 4 k1,k2,k3,k4 k 1 ,k 2 ,k 3 ,k 4 E 1 p i 1 n y 2 ik1 -λ ii - 1 n y ik1 y ik 1 -δ k1=k 1 λ ii × 1 p i 1 n y 2 ik2 -λ ii - 1 n y ik2 y ik 2 -δ k2=k 2 λ ii × 1 p i 1 n y 2 ik3 -λ ii - 1 n y ik3 y ik 3 -δ k3=k 3 λ ii × 1 p i 1 n y 2 ik4 -λ ii - 1 n y ik4 y ik 4 -δ k4=k 4 λ ii . (23) E[(m -µ) 4 ] = 1 (n -1) 4 k1,k2,k3,k4 k 1 ,k 2 ,k 3 ,k 4 E 4 s=1 1 pn i (y 2 iks -λ ii ) -4E 3 s=1 1 pn i (y 2 iks -λ ii ) 1 pn i (y ik4 y ik 4 -δ k4=k 4 λ ii ) + 6E 2 s=1 1 pn i (y 2 iks -λ ii ) 4 t=3 1 pn i (y ikt y ik t -δ kt=k t λ ii ) -4E 1 pn i (y 2 ik1 -λ ii ) 4 t=2 1 pn i (y ikt y ik t -δ kt=k t λ ii ) + E 4 t=1 1 pn i (y ikt y ik t -δ kt=k t λ ii ) . ( 24
) If k 1 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 ,
E y ik y ik -δ kt=k t λ ii 4 = E (y ik y ik -E[y ik y ik ]) 4 = E[(y ik y ik ) 4 ] -4E[(y ik y ik ) 3 ]E[y ik y ik ] + 6E[(y ik y ik ) 2 ]E[y ik y ik ] 2 -3E[y ik y ik ] 4 . ( 25
) If k = k then E[(y ik y ik ) 3 ]E[y ik y ik ] ≥ 0 and if k = k , E[(y ik y ik ) 3 ]E[y ik y ik ] = 0 ≥ 0, so: E y ik y ik -δ kt=k t λ ii 4 ≤ E[(y ik y ik ) 4 ] + 6E[(y ik y ik ) 2 ]E[y ik y ik ] 2 ≤ E[(y ik y ik ) 4 ] + 6E[(y ik y ik ) 4 ] (Cauchy-Schwarz) E y ik y ik -δ kt=k t λ ii 4 ≤ 7E[(y ik y ik ) 4 ]. (26)
Back to the bound of our expectation, let N ∈ 1, 4 , and we have:

E N s=1 1 pn i (y 2 iks -λ ii ) 4 t=N +1 1 pn i (y ikt y ik t -δ kt=k t λ ii ) ≤ 1 n 4 N s=1 4 E   1 p i (y 2 iks -λ ii ) 4   4 t=N +1 4 E   1 p i (y ikt y ik t -δ kt=k t λ ii ) 4   (2 Cauchy-Schwarz) ≤ 1 n 4 N s=1 4 1 p i E[(y 2 iks -λ ii ) 4 ] 4 t=N +1 4 1 p i E[(y ikt y ik t -δ kt=k t ) 4 ] (Jensen) ≤ 7 n 4 N s=1 4 1 p i E[y 8 iks ] 4 t=N +1 4 1 p i E[y 4 ikt y 4 ik t ] (Previous remark) ≤ 7 n 4 N s=1 4 1 p i E[y 8 iks ] 4 t=N +1 4 1 p i E[y 8 ikt ] (Jensen if k t = k t ) ≤ 7 n 4 4 s=1 4 1 p i E[y 8 iks ] ≤ 7 n 4 1 p i E[y 8 ik1 ] ≤ 7K 2 n 4 . (27)
So, in conclusion of this proof,

E[(m -µ) 4 ] ≤ 7(1 + 4 + 6 + 4 + 1) 28n 7 + O(n 6 ) n 4 (n -1) 4 K 2 -→ n→∞ 0.
(28)

Proof of Corollary 2

We have:

V[m 2 ] = E[m 4 ] -E[m 2 ] 2 = E[(m -µ + µ) 4 ] -E[(m -µ + µ) 2 ] 2 = E[(m -µ) 4 ] + 4E[(m -µ) 3 ]µ + 6E[(m -µ) 2 ]µ 2 + µ 4 -E[(m -µ) 2 ] 2 -2E[(m -µ) 2 ]µ 2 -µ 4 V[m 2 ] = E[(m -µ) 4 ] + 4E[(m -µ) 3 ]µ + 4E[(m -µ) 2 ]µ 2 -E[(m -µ) 2 ] 2 . ( 29
)
And, from Lemma 2, E[(m -µ) 4 ]→0, and so goes for the smaller moments by Jensen inequality, and µ is bounded from Lemma 1. So, V[m 2 ]→0.

6.6 Proof of Lemma 3

Preliminary combinatorial result

Let K ∈ N * , and

K indices (k 1 , ..., k K ) ∈ 1, n K .
Let's associate a graph with K vertices V = {1, ..., K} to this set of indices. The set of edges E is built as following: there is an edge between the node a ∈ V and b ∈ V, a = b (we don't allow self-loops), if the corresponding indices are equal, i.e if k a = k b . We finally define our graph G = (V, E).

Proposition 1. Let G = (V, E) a graph with K vertices generated from some indices (k

(0) 1 , ..., k (0)
K ) ∈ 1, n K with the procedure described previously. Suppose G has C ∈ 1, K connected components. Then, there are

C-1 i=0 (n -i) set of indices (k 1 , ..., k K ) ∈ 1, n K which have the associated graph G.
For each node v ∈ V, v belongs to a unique connected component that we denote c(v) ∈ 1, C . Then, the function:

x ∈ x ∈ 1, n C , x 1 , ..., x C all different → (x c(1) , ...., x c(K) ) is a bijection between x ∈ 1, n C , x 1 , ..., x C all different and
(k 1 , ..., k K ) ∈ 1, n K which have the associated graph G . Immediately, we deduce that its cardinal is equal to

C-1 i=0 (n -i).

Proof of Lemma 3

From the proof of Lemma 3.3 in (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], we have:

d 2 -δ 2 = -(m -µ) 2 + 2µ(µ -m) + S 2 -E[ S 2 ] . (30)
The first two terms converge to 0 in quadratic mean thanks to Lemma 2. Let's show that the last term converges to 0 in quadratic mean too, i.e V[ S 2 ] -→ n→∞ 0 .

21

Decomposing S 2 , we have:

S 2 = 1 p 1 n -1 k X •k -X X •k -X T 2 = 1 p 1 n -1 k Y •k -Ȳ Y •k -Ȳ T 2 = 1 p(n -1) 2 i,j k (y ik -ȳi )(y jk -ȳj ) 2 = 1 p(n -1) 2 i,j   1 n k,k y ik (y jk -y jk )   2 (Identity 1) = 1 pn 2 (n -1) 2 i,j   k1,k 1 =k1 y ik1 (y jk1 -y jk 1 )     k2,k 2 =k2 y ik2 (y jk2 -y jk 2 )   = p n 2 (n -1) 2 k1,k 1 =k1 k2=k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk2 -y jk1 y jk 2 -y jk 1 y jk2 + y jk 1 y jk 2 + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk2 -y jk1 y jk 2 -y jk 1 y jk2 + y jk 1 y jk = p n 2 (n -1) 2 k1,k 1 =k1 k2=k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 1 -y jk1 y jk 2 -y 2 jk 1 + y jk 1 y jk 2 + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk2 -y jk1 y jk 2 -y jk 1 y jk2 + y jk 1 y jk = p n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1   1 p i y ik1 y ik 1 2 - 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 + y 2 jk 1 -y jk 1 y jk   + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2   1 p i y ik1 y ik2 2 - 1 p 2 i,j y ik1 y ik2 y jk1 y jk 2 + y jk 1 y jk2 -y jk y jk 2   . (31)
S 2 = p n 2 (n -1) k1,k 1 =k1 1 p i y ik1 y ik 1 2 (t1) - p n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 (t2) - p n 2 (n -1) k1,k 1 =k1 1 p 2 i,j y ik1 y ik 1 y 2 jk 1 (t3) + p n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1 1 p 2 i,j y ik1 y ik 1 y jk 1 y jk 2 (t4) + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p i y ik1 y ik2 2 (t5) - p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk 2 (t6) - p n 2 (n -1) k1,k 1 =k1,k2 =k 1 1 p 2 i,j y ik1 y ik2 y jk 1 y jk2 (t7) + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk 1 y jk 2 (t8). ( 32
)
We notice that, re-injecting the missing terms into the sum:

(t2) = - p n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 (t2) = - p n 2 (n -1) 2 k1,k 1 ,k 2 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 - 1 n y jk 1 - 1 n y 2 ik1 y jk1 y jk 2 - 1 n y jk1 . (33)
Similarly, we have:

(t4) = p n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1 1 p 2 i,j y ik1 y ik 1 y jk 1 y jk 2 (t4) = p n 2 (n -1) 2 k1,k 1 ,k 2 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 - 1 n y jk1 - 1 n y 2 ik1 y jk1 y jk 2 - 1 n y jk1 . (34) So, (t2) + (t4) = p n 2 (n -1) 2 k1,k 1 =k1 1 p 2 i,j y ik1 y ik 1 y jk1 (y jk 1 -y jk1 ) (t2) + (t4) = 1 n -1 ((1) + (3)) . (35)
And,

(t5) = p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p i y ik1 y ik2 2 = p(n -2) n 2 (n -1) k1,k2 =k1 1 p i y ik1 y ik2 2 + p n 2 k1 1 p i y 2 ik1 2 (t5) = (n -2) × (t1) + p n 2 k1 1 p i y 2 ik1 2 . (36)
Moreover,

(t6) = - p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk 2 = - p(n -2) n 2 (n -1) 2 k1,k2 =k1 k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk1 y jk 2 - p n 2 (n -1) k1 k 2 =k1 1 p 2 i,j y 2 ik1 y jk1 y jk 2 (t6) = (n -2) × (t2) + (t3). (37) And, (t7)
= - p n 2 (n -1) k1,k 1 =k1,k2 =k 1 1 p 2 i,j y ik1 y ik2 y jk 1 y jk2 = - p n 2 (n -1) k2,k 1 =k2,k1 =k 1 1 p 2 i,j
y ik1 y ik2 y jk 1 y jk2 (re-indexing) (t7) = (n -1) × (t2).

(38) Finally, we obtain,

S 2 = p n 2 k1 1 p i y 2 ik1 2 (t1 ) + p(n 2 -2n + 2) n 2 (n -1) 2 k1,k 1 =k1 1 p i y ik1 y ik 1 2 (t2 ) - p(2n -3) n 2 (n -1) 2 k1,k 1 =k1,k 2 =k 1 1 p 2 i,j y ik1 y ik 1 y jk1 y jk 2 (t3 ) - p(2n -1) n 2 (n -1) 2 k1,k 1 =k1 1 p 2 i,j y ik1 y ik 1 y 2 jk 1 (t4 ) + p n 2 (n -1) 2 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 1 p 2 i,j y ik1 y ik2 y jk 1 y jk 2 (t5 ). ( 39
)
It is sufficient to show that the variance of each of the 5 term converges to 0 as n goes to infinity in order to prove that V[ S 2 ] converges to 0 as n goes to infinity.

• V[(t1 )] → 0 immediately from the proof of Lemma 3.3 in (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF].

• V[(t2 )] → 0 immediately from the proof of Lemma 3.3 in (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF].

• Let's prove that V[(t3 )] → 0. V[(t3 )] = p 2 (2n -3) 2 n 4 (n -1) 4 k1,k 1 =k1 k 2 =k 1 k3,k 3 =k3 k 4 =k 3 i,j,l,m 1 p 4
Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ).

(40)

Let (k 1 , k 1 , k 2 , k 3 , k 3 , k 4 ) ∈ 1, n 6
respecting the conditions given in the sums. Suppose there exists (i, j, l, m) ∈ 1, p 4 such that Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = 0. Let consider the graph G = (V, E) built from (k 1 , k 1 , k 2 , k 3 , k 3 , k 4 ) following the procedure described in the preliminary combinatorial result.

As for some (i, j, l, m) ∈ 1, p 4 we have Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = 0, by independence, the nodes 1', 2', 3', and 4' can't be isolated. As a consequence, G has at least 2 edges.

-When the graph G has only 2 edges, we have either of the following conditions, that we denote the (*) conditions:

* (k 1 = k 3 ) ∧ (k 2 = k 4 ) ∧ (k 1 = k 3 ) ∧ (k 1 = k 2 ) ∧ (k 3 = k 4 ), * (k 1 = k 4 ) ∧ (k 3 = k 2 ) ∧ (k 1 = k 3 ) ∧ (k 1 = k 2 ) ∧ (k 3 = k 4 ).
Note that the case where (k Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = λ ij λ lm λ im λ jl .

1 = k 2 ) ∧ (k 3 = k 4 )
(42)

Using the fact that i = j =⇒ λ ij = 0, in both cases we have:

i,j,l,m
Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) =

i λ 4 ii = p i 1 p E[y 2 i1 ] 4 ≤ p i 1 p E[y 8 i1 ] ≤ p i 1 p E[y 8 i1 ] i,j,l,m
Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) ≤ pK 2 .

(

) 43
Under the conditions (*), G has exactly 4 connected components. So, from the preliminary combinatorial result, there are 2n(n -1)(n -2)(n -3) different 6-uples of indices k respecting the (*) conditions. We finally have:

p 2 (2n -3) 2 n 4 (n -1) 4 (k1,k 1 ,k 2 ,k3,k 3 ,k 4 ) under (*) conditions i,j,l,m 1 p 4 Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) ≤ p 2 (2n -3) 2 n 4 (n -1) 4 × 2n(n -1)(n -2)(n -3) 1 p 4 pK 2 ≤ 2 (2n -3) 2 pn 3 (n -1) K 2 ≤ 8 n(n -1) K 2 -→ n→∞ 0. (44)
-Otherwise, G has 3 edges or more: we denote it as the (**) condition. As it has only 6 vertices, there are at most 3 connected component. So, from the preliminary combinatorial result, there are n(n -1)(n -2) + n(n -1) + n = n((n -1)

1 p 4 V[y ik1 y ik 1 y jk1 y jk 2 ]V[y lk3 y lk 3 y mk3 y mk 4 ] ≤ i,j,l,m 1 p 4 E[y 2 ik1 y 2 ik 1 y 2 jk1 y 2 jk 2 ]E[y 2 lk3 y 2 lk 3 y 2 mk3 y 2 mk 4 ] ≤ i,j,l,m 1 p 4 8 E[y 8 ik1 ]E[y 8 ik 1 ]E[y 8 jk1 ]E[y 8 jk 2 ]E[y 8 lk3 ]E[y 8 lk 3 ]E[y 8 mk3 ]E[y 8 mk 4 ] ≤ i,j,l,m 1 p 4 4 E[y 8 i1 ]E[y 8 j1 ]E[y 8 l1 ]E[y 8 m1 ] ≤ i 1 p 4 E[y 8 i1 ] 4 ≤ i 1 p E[y 8 i1 ] ≤ K 2 .
(45)

So, using both of the previous inequalities,

p 2 (2n -3) 2 n 4 (n -1) 4 (k1,k 1 ,k 2 ,k3,k 3 ,k 4 ) under (**) condition i,j,l,m 1 p 4 Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) ≤ p 2 (2n -3) 2 n 4 (n -1) 4 × N n((n -1) 2 + n)K 2 ≤ (2n -3) 2 n 2 (n -1) 4 × N n((n -1) 2 + n)K 2 1 K 2 -→ n→∞ 0. ( 46
)
So, from both previous cases, we immediately have:

V[(t3 )] -→ n→∞ 0. ( 47
) • Let's prove that V[(t4 )] → 0. V[(t4 )] = p 2 (2n -1) 2 n 4 (n -1) 4 k1 k 1 =k1 k3 k 3 =k3 i,j,l,m 1 p 4 Cov(y ik1 y ik 1 y 2 jk 1 , y lk3 y lk 3 y 2 mk 3 ). ( 48
)
Similarly as the previous case, let (k 1 , k 1 , k 3 , k 3 ) ∈ 1, n 4 respecting the conditions given in the sums. Suppose there exists (i, j, l, m) ∈ 1, p 4 such that Cov(y ik1 y ik 1 y 2 jk 1

, y lk3 y lk 3 y 2 mk 3 ) = 0. Let consider the associated graph G = (V, E). By independence, the nodes 1 and 3 can't be isolated. As a consequence, G has at most 3 connected components. Denoting N the number of graphs with at most 3 connected components and 4 vertices, from the preliminary combinatorial result, there are at most N n((n -1) and(i, j, l, m) ∈ 1, p 4 such that Cov(y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 ) = 0. Let's consider the associated graph G, built following the procedure described in the preliminary combinatorial result. For each node a ∈ V, the connected component of G containing a contains at least 2 nodes (a and at least an other one). Otherwise, a is isolated, which means that k a is different from all the other indices, and so by independence, Cov(y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 ) = 0, which is in contradiction which our hypothesis. As each connected component contains at least 2 nodes, there are at most 4 connected components in G. So, from the preliminary combinatorial result, there are at most n(n-1)(n-2)(n-3)+n(n-

2 + 1) different indices (k 1 , k 1 , k 3 , k 3 ) ∈ 1, n 4 such that i,j,l,m Cov(y ik1 y ik 1 y 2 jk 1 , y lk3 y lk 3 y 2 mk 3 ) = 0. Moreover, i,j,l,m 1 p 4 |Cov(y ik1 y ik 1 y 2 jk 1 , y lk3 y lk 3 y 2 mk 3 )| ≤ i,j,l,m 1 p 4 V[y 2 jk 1 y ik1 y ik 1 ]V[y 2 mk 3 y lk3 y lk 3 ] ≤ i,j,l,m 1 p 4 E[y 4 jk 1 y 2 ik1 y 2 ik 1 ]E[y 4 mk 3 y 2 lk3 y 2 lk 3 ] ≤ i,j,l,m 1 p 4 4 E[y 8 jk 1 ] E[y 8 ik1 ]E[y 8 ik 1 ]E[y 8 mk 3 ] E[y 8 lk3 ]E[y 8 lk 3 ] ≤ i,j,l,m 1 p 4 4 E[y 8 j1 ] E[y 8 i1 ]E[y 8 i1 ]E[y 8 m1 ] E[y 8 l1 ]E[y 8 l1 ] ≤ i,j,l,m 1 p 4 4 E[y 8 i1 ]E[y 8 j1 ]E[y 8 l1 ]E[y 8 m1 ] ≤ i 1 p 4 E[y 8 i1 ] 4 ≤ i 1 p E[y 8 i1 ] i,j,l,m 1 p 4 |Cov(y ik1 y ik 1 y 2 jk 1 , y lk3 y lk 3 y 2 mk 3 )| ≤ K 2 . (49) So, V[(t4 )] ≤ p 2 (2n -1) 2 n 4 (n -1) 4 N n((n -1) 2 + 1)K 2 (2n -1) 2 n(n -1) 4 N ((n -1) 2 + 1)K 2 1 K 2 V[(t4 )] -→ n→∞ 0. ( 50
) • Let's prove that V[(t5 )] → 0. V[(t5 )] = p 2 n 4 (n -1) 4 k1,k 1 =k1 k2 =k 1 ,k 2 =k2 k3,k 3 =k3 k4 =k 3 ,k 4 =k4 i,j,l,m 1 p 4 Cov(y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 ). (51) Let (k 1 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 , k 4 ) ∈ 1, n 8
1)(n -2) + n(n -1) + n = n (n -1)(n -2) 2 + n different (k 1 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 , k 4 ) ∈
1, n 8 which have the same associated graph G.

Moreover, there is a finite number N ∈ N (independent of n) of graphs with 8 nodes and at least 4 connected components.

Finally, combining the previous 2 counting, we deduce that there are at most N n (n -1)(n -2) 2 + n terms such that i,j,l,m 1 p 4 Cov(y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 ) = 0. As previously, we have also, i,j,l,m

1 p 4 |Cov(y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 )| ≤ K 2 . (52)
Finally,

V[(t5 )] ≤ p 2 n 4 (n -1) 4 N n (n -1)(n -2) 2 + n K 2 ≤ N ((n -1)(n -2) 2 + n) n(n -1) 3 K 2 1 K 2 V[(t5 )] -→ n→∞ 0. (53)
We showed that each of the 5 terms of S 2 have a variance that converges to 0 as n goes to infinity. So,

V[ S 2 ] -→ n→∞ 0. ( 54
)
Which concludes the proof of the first part of the lemma:

d 2 -δ 2 q.m -→ n→∞ 0. (55)
Finally, by property of the expectation, we have that

E[(d 2 -E[d 2 ]) 2 ] ≤ E[(d 2 -δ 2 ) 2 ]
, so it follows that:

d 2 -E[d 2 ] q.m -→ n→∞ 0. (56)
6.7 Proof of Lemma 4

E[ b2 ] = E 1 n 2 k n n -1 x•k (x •k ) t -S 2 = E 1 n 2 k n n -1 x•k (x •k ) t -Σ + Σ -S 2 = E 1 n 2 k n n -1 x•k (x •k ) t -Σ 2 - 1 n E S -Σ 2 = E 1 n 2 k n n -1 ỹ•k (ỹ •k ) t -Λ 2 - 1 n β 2 = 1 pn 2 i,j k E n n -1 (y ik -ȳi ) (y jk -ȳj ) -λ ij 2 - 1 n β 2 = n p(n -1) 2 i,j E (y i1 -ȳi ) (y j1 -ȳj ) - n -1 n λ ij 2 - 1 n β 2 = n p(n -1) 2 i,j E (y i1 y j1 -λ ij ) -ȳi y j1 - 1 n λ ij -ȳj y i1 - 1 n λ ij + ȳi ȳj - 1 n λ ij 2 - 1 n β 2 E[ b2 ] = n p(n -1) 2 i,j V [y i1 y j1 ] + V [ȳ i y j1 ] + V [ȳ j y i1 ] + V [ȳ i ȳj ] -2E (y i1 y j1 -λ ij ) ȳi y j1 - 1 n λ ij -2E (y i1 y j1 -λ ij ) ȳj y i1 - 1 n λ ij + 2E (y i1 y j1 -λ ij ) ȳi ȳj - 1 n λ ij + 2E ȳi y j1 - 1 n λ ij ȳj y i1 - 1 n λ ij -2E ȳi y j1 - 1 n λ ij ȳi ȳj - 1 n λ ij -2E ȳj y i1 - 1 n λ ij ȳi ȳj - 1 n λ ij - 1 n β 2 . (57)
We simplify each of the terms into manageable quantities:

• V [ȳ i y j1 ] = 1 n 2 V[y i1 y j1 ] + n-1 n 2 λ ii λ jj , • V [ȳ j y i1 ] = 1 n 2 V[y i1 y j1 ] + n-1 n 2 λ ii λ jj , • V [ȳ i ȳj ] = 1 n 3 V[y i1 y j1 ] + (n -1)(λ ii λ jj + λ 2 ij ) , • E (y i1 y j1 -λ ij ) ȳi y j1 -1 n λ ij = 1 n V[y i1 y j1 ], • E (y i1 y j1 -λ ij ) ȳj y i1 -1 n λ ij = 1 n V[y i1 y j1 ], • E (y i1 y j1 -λ ij ) ȳi ȳj -1 n λ ij = 1 n 2 V[y i1 y j1 ], • E ȳi y j1 -1 n λ ij ȳj y i1 -1 n λ ij = 1 n 2 V[y i1 y j1 ] + n-1 n 2 λ 2 ij , • E ȳi y j1 -1 n λ ij ȳi ȳj -1 n λ ij = 1 n 3 V[y i1 y j1 ] + (n -1)(λ ii λ jj + λ 2 ij ) , • E ȳj y i1 -1 n λ ij ȳi ȳj -1 n λ ij = 1 n 3 V[y i1 y j1 ] + (n -1)(λ ii λ jj + λ 2 ij ) .
Adding up all the terms, we obtain:

E b2 = 1 p(n -1) i,j n 2 -3n + 3 n 2 V[y i1 y j1 ] + 2n -3 n 2 λ ii λ jj + 2n -3 n 2 λ 2 ij - 1 n β 2 = 1 p(n -1)   n 2 -3n + 3 n 2 i,j V[y i1 y j1 ] + p 2 2n -3 n 2 Σ, I 2 + p 2n -3 n 2 Σ 2   - 1 n β 2 E b2 = 1 p(n -1)   n 2 -3n + 3 n 2 i,j V[y i1 y j1 ] + p 2 2n -3 n 2 µ 2 + p 2n -3 n 2 (α 2 + µ 2 )   - 1 n β 2 . (58)
From the proof of Lemma 1, we have:

1 pn i,j V[y i1 y j1 ] = β 2 - p + 1 n(n -1) µ 2 - 1 n(n -1) α 2 . (59)
Denoting:

γ n = n(n-1) n 2 -3n+3 , λ n = n 2 (n-2) (n-1)(n 2 -3n+3) , c 0 = 1 γn -1 n -λn γnn 2 , c 1 = λn γnn 2 , c 2 = (p + 1
)c 1 , we notice that:

c 0 + c 1 = n 2 -3n + 3 n(n -1) - 1 n , c 1 = - n 2 -3n + 3 n 2 (n -1) 2 + 2n -3 n 2 (n -1) , c 2 = - (n 2 -3n + 3)(p + 1) n 2 (n -1) 2 + (2n -3)(p + 1) n 2 (n -1) . (60)
We obtain from the previous lines:

E[ b2 ] = (c 0 + c 1 )β 2 + c 1 α 2 + c 2 µ 2 .
So, we can conclude using

δ 2 = α 2 + β 2 : E[ b2 ] = c 0 β 2 + c 1 δ 2 + c 2 µ 2 .

Proof of Lemma 5

We compute the expectations of m 2 and d 2 .

E[m 2 ] = µ 2 + V[m]. (61)
And,

E[d 2 ] = E[ S -mI 2 ] = E[ S -µI 2 + µI -mI 2 + 2 S -µI, µI -mI ] = δ 2 + V[m] -2V[m] E[d 2 ] = δ 2 -V[m]. (62)
Moreover, from Lemma 4, we have:

E[ b2 ] = c 0 β 2 n + c 1 δ 2 n + c 2 µ 2 n .
So, combining the last 3 equations, we obtain:

E[ b2 ] = c 0 β 2 + c 1 E[d 2 ] + c 2 E[m 2 ] + (c 1 -c 2 )V[m]. 6.9 Proof of Lemma 6 b2 -E[ b2 ] = 1 n 2 k n n -1 x•k xT •k -Σ 2 -E[ b2 ] + 1 n 2 k n n -1 x•k xT •k -S 2 - 1 n 2 k n n -1 x•k xT •k -Σ 2 . (63)
Firstly, we want to show that the variance of the first term converges to 0 a n goes to infinity. Following the decomposition developed in Lemma 4, we have:

V 1 n 2 k n n -1 x•k xT •k -Σ 2 = V 1 (n -1) 2 k ỹ•k ỹT •k - n -1 n Λ 2 = V   1 p(n -1) 2 k i,j (y ik -ȳi )(y jk -ȳj ) - n -1 n λ ij 2   = V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) -(ȳ i y jk - 1 n λ ij ) -(y ik ȳj - 1 n λ ij ) + (ȳ i ȳj - 1 n λ ij ) 2   = V 1 p(n -1) 2 k i,j [y ik y jk -λ ij ] 2 + ȳi y jk - 1 n λ ij 2 + ȳj y ik - 1 n λ ij 2 + ȳi ȳj - 1 n λ ij 2 -2 (y ik y jk -λ ij ) ȳi y jk - 1 n λ ij -2 (y ik y jk -λ ij ) ȳj y ik - 1 n λ ij + 2 (y ik y jk -λ ij ) ȳi ȳj - 1 n λ ij + 2 ȳi y jk - 1 n λ ij ȳj y ik - 1 n λ ij -2 ȳi y jk - 1 n λ ij ȳi ȳj - 1 n λ ij -2 ȳj y ik - 1 n λ ij ȳi ȳj - 1 n λ ij . (64)
We then prove that the variance of each of the 10 separated sums converges to 0. For that, let's show a useful inequality.

Preliminary inequality

Let (k 1 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 , k 4 ) ∈ 1,
1 p 4 Cov y ik1 y jk 1 -λ ij 1 k1=k 1 y ik2 y jk 2 -λ ij 1 k2=k 2 , y lk3 y mk 3 -λ lm 1 k3=k 3 y lk4 y mk 4 -λ lm 1 k4=k 4 ≤ 9K 2 .
(66)

Variance of the first term

Now, we can prove that the variance of each of the 10 separated sums converges to 0.

• We have:

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) 2   = p 2 (n -1) 4 k1,k2 i,j,l,m 1 p 4 Cov (y ik1 y jk1 -λ ij )
2 , (y lk2 y mk2 -λ lm ) 2 .

(67)

If k 1 = k 2 , the covariance is trivially null. So,

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) 2   = p 2 (n -1) 4 k1 i,j,l,m 1 p 4 Cov (y ik1 y jk1 -λ ij ) 2 , (y lk1 y mk1 -λ lm ) 2 ≤ K 2 1 n 3 (n -1) 4 i,j,l,m 1 p 4 Cov (y i1 y j1 -λ ij )
2 , (y l1 y m1 -λ lm ) 2 .

(68)

And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov (y i1 y j1 -λ ij ) 2 , (y l1 y m1 -λ lm ) 2 ≤ 9K 2 . (69)
Finally, we can conclude,

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) 2   ≤ n 3 (n -1) 4 × 9K 2 1 K 2 V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) 2   -→ n→∞ 0.
(70)

• We have:

V   1 p(n -1) 2 k i,j ȳi y jk - 1 n λ ij 2   = p 2 (n -1) 4 n 4 k1,k 1 ,k 1 k2,k 2 ,k 2 i,j,l,m 1 p 4 Cov y ik 1 y jk1 -λ ij 1 k 1 =k1 y ik 1 y jk1 -λ ij 1 k 1 =k1 , y lk 2 y mk2 -λ lm 1 k 2 =k2 y lk 2 y mk2 -λ lm 1 k 2 =k2 . ( 71
) If k 1 , k 1 , k 1 , k 2 , k 2 , k 2
are all different, then, the covariance is null. So there are at most n 6 -n(n -1)(n -2)(n -3)(n -4)(n -5) non-zero terms in the sum over k-indices.

And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov y ik 1 y jk1 -λ ij 1 k 1 =k1 y ik 1 y jk1 -λ ij 1 k 1 =k1 , y lk 2 y mk2 -λ lm 1 k 2 =k2 y lk 2 y mk2 -λ lm 1 k 2 =k2 ≤ 9K 2 . (72) So, V   1 p(n -1) 2 k i,j ȳi y jk - 1 n λ ij 2   ≤ p 2 (n 6 -n(n -1)(n -2)(n -3)(n -4)(n -5)) (n -1) 4 n 4 × 9K 2 ≤ n 2 (n 6 -n(n -1)(n -2)(n -3)(n -4)(n -5)) (n -1) 4 n 4 × 9K 2 1 K 2 V   1 p(n -1) 2 k i,j ȳi y jk - 1 n λ ij 2   -→ n→∞ 0. ( 73
)
• Immediately from the previous point,

V   1 p(n -1) 2 k i,j ȳj y ik - 1 n λ ij 2   -→ n→∞ 0. (74)
• We have:

V   1 p(n -1) 2 k i,j ȳi ȳj - 1 n λ ij 2   = p 2 (n -1) 4 n 6 k1,k 1 ,k2,k 2 k3,k 3 ,k4,k 4 i,j,l,m 1 p 4 Cov y ik1 y jk 1 -λ ij 1 k1=k 1 y ik2 y jk 2 -λ ij 1 k2=k 2 , y lk3 y mk 3 -λ lm 1 k3=k 3 y lk4 y mk 4 -λ lm 1 k4=k 4 . ( 75
) If k 1 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 , k 4
are all different, then, the covariance is null. So there are at most n 8 -n(n -1)(n -2)(n -3)(n -4)(n -5)(n -6)(n -7) non-zero terms in the sum over k-indices.

And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov y ik1 y jk 1 -λ ij 1 k1=k 1 y ik2 y jk 2 -λ ij 1 k2=k 2 , y lk3 y mk 3 -λ lm 1 k3=k 3 y lk4 y mk 4 -λ lm 1 k4=k 4 ≤ 9K 2 . (76) So, V   1 p(n -1) 2 k i,j ȳi ȳj - 1 n λ ij 2   ≤ p 2 (n 8 -n(n -1)(n -2)(n -3)(n -4)(n -5)(n -6)(n -7)) (n -1) 4 n 6 × 9K 2 ≤ n 2 (n 8 -n(n -1)(n -2)(n -3)(n -4)(n -5)(n -6)(n -7)) (n -1) 4 n 6 × 9K 2 1 K 2 -→ n→∞ 0.
(77)

• We have:

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) ȳi y jk - 1 n λ ij   = p 2 (n -1) 4 n 2 k1,k 1 ,k2,k 2 i,j,l,m 1 p 4 Cov (y ik1 y jk1 -λ ij ) y ik 1 y jk1 -λ ij 1 k1=k 1 , (y lk2 y mk2 -λ lm ) y lk 2 y mk2 -λ lm 1 k2=k 2 . (78)
If k 1 / ∈ {k 1 , k 2 , k 2 } , then, the covariance is null. So there are at most n 4 -n(n -1) 3 non-zero terms in the sum over k-indices. And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov (y ik1 y jk1 -λ ij ) y ik 1 y jk1 -λ ij 1 k1=k 1 , (y lk2 y mk2 -λ lm ) y lk 2 y mk2 -λ lm 1 k2=k 2 ≤ 9K 2 . (79) So, V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) ȳi y jk - 1 n λ ij   ≤ p 2 (n 4 -n(n -1) 3 ) (n -1) 4 n 2 × 9K 2 ≤ n 2 (n 4 -n(n -1) 3 ) (n -1) 4 n 2 × 9K 2 1 K 2 -→ n→∞ 0. (80)
• Immediately from the previous point,

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) ȳj y i1 - 1 n λ ij   -→ n→∞ 0. (81)
• We have:

V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) ȳi ȳj - 1 n λ ij   = p 2 (n -1) 4 n 4 k1,k 1 ,k 1 k2,k 2 ,k 2 i,j,l,m 1 p 4 Cov (y ik1 y jk1 -λ ij ) y ik 1 y jk 1 -λ ij 1 k 1 =k 1 , (y lk2 y mk2 -λ lm ) y lk 2 y mk 2 -λ lm 1 k 2 =k 2 . ( 82
) If k 1 / ∈ {k 1 , k 1 , k 2 , k 2 , k 2
} , then, the covariance is null. So there are at most n -n(n -1) 5 non-zero terms in the sum over k-indices. And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov (y ik1 y jk1 -λ ij ) y ik 1 y jk 1 -λ ij 1 k 1 =k 1 , (y lk2 y mk2 -λ lm ) y lk 2 y mk 2 -λ lm 1 k 2 =k 2 ≤ 9K 2 . (83) So, V   1 p(n -1) 2 k i,j (y ik y jk -λ ij ) ȳi ȳj - 1 n λ ij   ≤ p 2 (n 6 -n(n -1) 5 (n -1) 4 n 4 × 9K 2 ≤ n 2 (n 6 -n(n -1) 5 ) (n -1) 4 n 4 × 9K 2 1 K 2 -→ n→∞ 0. (84)
• We have:

V   1 p(n -1) 2 k i,j ȳi y jk - 1 n λ ij y ik ȳj - 1 n λ ij   = p 2 (n -1) 4 n 4 k1,k 1 ,k 1 k2,k 2 ,k 2 i,j,l,m 1 p 4 Cov y ik 1 y jk1 -λ ij 1 k 1 =k1 y ik1 y jk 1 -λ ij k1=k 1 , y lk 2 y mk2 -λ lm 1 k2=k 2 y lk2 y mk 2 -λ lm 1 k2=k 2 . (85) If k 1 / ∈ {k 1 , k 1 , k 2 , k 2 , k 2
} , then, the covariance is null. So there are at most n -n(n -1) 5 non-zero terms in the sum over k-indices.

And we have from the preliminary inequality, i,j,l,m

1 p 4 Cov y ik 1 y jk1 -λ ij 1 k 1 =k1 y ik1 y jk 1 -λ ij 1 k1=k 1 , y lk 2 y mk2 -λ lm 1 k2=k 2 y lk2 y mk 2 -λ lm 1 k2=k 2 ≤ 9K 2 . (86) So, V   1 p(n -1) 2 k i,j (ȳ i y jk -λ ij ) y ik ȳj - 1 n λ ij   ≤ p 2 (n 6 -n(n -1) 5 ) (n -1) 4 n 4 × 9K 2 ≤ n 2 (n 6 -n(n -1) 5 ) (n -1) 4 n 4 × 9K 2 1 K 2 -→ n→∞ 0. (87)
• Immediately from the fourth point,

V   1 p(n -1) 2 k i,j ȳi y jk - 1 n λ ij ȳi ȳj - 1 n λ ij   -→ n→∞ 0. (88)
• Immediately from the previous point,

V   1 p(n -1) 2 k i,j ȳj y ik - 1 n λ ij ȳi ȳj - 1 n λ ij   -→ n→∞ 0. (89)
So, the first term has its variance converging to 0:

V 1 n 2 k n n -1 x•k xT •k -Σ 2 -→ n→∞ 0. (90)

Variance of the second term

The variance of the second term converges to 0 as n goes to infinity. Indeed, from the proof of Lemma 3.4 in [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [10]:

1 n 2 k n n -1 x•k xT •k -S 2 - 1 n 2 k n n -1 x•k xT •k -Σ 2 = 1 n S -Σ 2 . ( 91
)
And we have:

E S -Σ 4 ≤ E ( S -µI + µI -Σ ) 4 . (92)
E[ S-µI 4 ] and µI-Σ are bounded, from Lemma 2 and Lemma 1 respectively, so

E S -Σ 4 is bounded. Consequently, V 1 n S -Σ 2 -→ n→∞ 0. To conclude, V b2 -→ n→∞ 0.
6.10 Proof of Lemma 7

V[m] = E   1 p i 1 n -1 k (y ik -ȳi ) 2 - 1 p i λ ii 2   = E   1 p i 1 n -1 k (y ik -ȳi ) 2 - n -1 n λ ii 2   = 1 (n -1) 2 E   k1 k2 1 p i (y ik1 -ȳi ) 2 - n -1 n λ ii   1 p j (y jk2 -ȳj ) 2 - n -1 n λ jj     = 1 p 2 (n -1) 2 k1,k2 i,j E y 2 ik1 -2y ik1 ȳi + ȳ2 i - n -1 n λ ii y 2 jk2 -2y jk2 ȳj + ȳ2 j - n -1 n λ jj . (93)
Fully developing the terms, we obtain,

V[m] = 1 p 2 (n -1) 2 k1 k2 i,j E y 2 ik1 -λ ii y 2 jk2 -λ jj (t1) -2E y 2 ik1 -λ ii y jk2 ȳj - 1 n λ jj (t2) + E y 2 ik1 -λ ii ȳ2 j - 1 n λ jj (t3) -2E y ik1 ȳi - 1 n λ ii y 2 jk2 -λ jj (t4) + 4E y ik1 ȳi - 1 n λ ii y jk2 ȳj - 1 n λ jj (t5) -2E y ik1 ȳi - 1 n λ ii ȳ2 j - 1 n λ jj (t6) + E ȳ2 i - 1 n λ ii y 2 jk2 -λ jj (t7) -2E ȳ2 i - 1 n λ ii y jk2 ȳj - 1 n λ jj (t8) + E ȳ2 i - 1 n λ ii ȳ2 j - 1 n λ jj (t9). (94)
µ is bounded by Lemma 1 and V[m] → 0 by Lemma 2, so µ 2 V[m] → 0. Considering the other term, we have:

V[ S, Σ ] = V 1 n -1 Ỹ Ỹ T , Λ = V 1 p(n -1) i λ ii k (y ik -ȳi ) 2 = V   1 pn(n -1) i λ ii k,k =k y ik (y ik -y ik )   (Identity 1) = 1 n 2 (n -1) 2 k1,k 1 =k1 k2,k 2 =k2 Cov 1 p i λ ii y ik1 (y ik1 -y ik 1 ), 1 p i λ ii y ik2 (y ik2 -y ik 2 ) = 1 n 2 (n -1) 2 k1,k 1 =k1 k2,k 2 =k2 Cov 1 p i λ ii y 2 ik1 , 1 p i λ ii y 2 ik2 + 1 n 2 (n -1) 2 k1,k 1 =k1 k2,k 2 =k2 Cov 1 p i λ ii y ik1 y ik 1 , 1 p i λ ii y ik2 y ik 2 = 1 n 2 k1 Cov 1 p i λ ii y 2 ik1 , 1 p i λ ii y 2 ik1 + 1 n 2 (n -1) 2 k1,k 1 =k1 Cov 1 p i λ ii y ik1 y ik 1 , 1 p i λ ii y ik1 y ik 1 + 1 n 2 (n -1) 2 k1,k 1 =k1 Cov 1 p i λ ii y ik1 y ik 1 , 1 p i λ ii y ik 1 y ik1 = 1 n V 1 p i λ ii y 2 i1 + 2 n(n -1) V 1 p i λ ii y i1 y i2 ≤ 1 n E   1 p i λ ii y 2 i1 2   + 2 n(n -1) E   1 p i λ ii y i1 y i2 2   ≤ 1 n E 1 p i λ 2 ii 1 p i y 4 i1 + 2 n(n -1) E 1 p i λ 2 ii 1 p i y 2 i1 y 2 i2 = 1 n 1 p i E y 2 i1 2 1 p i E y 4 i1 + 2 n(n -1) 1 p i E y 2 i1 2 2 ≤ n + 1 n(n -1) 1 p i E y 4 i1 2 ≤ n + 1 n(n -1) 1 p i E y 8 i1 V[ S, Σ ] ≤ n + 1 n(n -1) K 2 . (116) So, V[ S, Σ ] → 0, and so E[(α 2 -α 2 ) 2 ] → 0. As E[(a 2 u -α 2 ) 2 ] → 0 by Lemma 9, it comes that E[(a 2 u -α 2 ) 2 ] → 0.
Therefore, the assumptions of Lemma A.1 of ( [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [10] are verified by u 2 = (a 2 u -α 2 ) 2 , τ 1 = 2 and τ 2 = 0. It proves that:

E (a 2 u -α 2 ) 2 d 2 → 0. ( 117
)
Backing up, we have shown that E[ S * -Σ * * 2 ] → 0. We complete the proof of the theorem with the following inequality:

E | S * -Σ 2 -Σ * * -Σ 2 ] = E [| S * -Σ * * , S * + Σ * * -2Σ |] ≤ E [ S * -Σ * * 2 ] E [ S * + Σ * * -2Σ 2 ]. ( 118
)
The first term converges to 0 as we showed above, and the second term is bounded because E[ S * -Σ 2 ] is bounded. So, the product converges to 0, which completes the proof.

Proof of Theorem 4

The proof from (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] can be applied as it is, because it uses only the results of Theorem 3 which are the same as Theorem 3.3 in [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF].

Proof of Theorem 5

Respectively from Lemma 2 and 3, m r -µ converges to 0 in quartic mean and d 2 r -δ 2 converges to 0 in quadratic mean. Let's define: (120) b2 -β 2 converges to 0 in quadratic mean by Lemma 6. Moreover, S 2 -E[ S 2 ] converges to 0 in quadratic mean by Lemma 3, and E[ S 2 ] = β 2 + α 2 + µ 2 , so E[ S 2 ] is bounded by Lemma 1. So, 1 n(n-1) 2 S 2 converges to 0 in quadratic mean, which finally proves that b2 r -β 2 converges to 0 in quadratic mean. Following the idea of proof of Lemma 10, we have:

-max(| b2 r -β 2 |, |d 2 r -δ 2 |) ≤ b 2 r -β 2 ≤ max(| b2 r -β 2 |, |d 2 r -δ 2 |). (121)
So,

E[(b 2 r -β 2 ) 2 ] ≤ E[max(| b2 r -β 2 |, |d 2 r -δ 2 |) 2 ] ≤ E[( b2 r -β 2 ) 2 ] + E[(d 2 r -δ 2 ) 2 ]. (122)
Which, as previously, leads to the 2 following results: E[(b 2 r -β 2 ) 2 ] → 0, E[(a 2 r -α 2 ) 2 ] → 0. Finally, we check that we have the set of hypotheses required by the proof of Theorem 2 to work:

• α 2 , β 2 , δ 2 are non-negative, bounded, and α 2 + β 2 = δ 2 ,

• m r -µ converges to 0 in quartic mean,

• d 2 r is nonnegative, and d 2 r -δ 2 -→ q.m 0,

• 0 ≤ a 2 r ≤ d 2 s and a 2 r -α 2 -→ q.m 0,

• a 2 r + b 2 r = d 2 r , with b 2 r ≥ 0.

Then, we can apply the result of the theorem 3.2 from (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], so identically we have that E[ S * r -Σ * 2 ] → 0 and S * r has the same asymptotic expected loss as E S * r -Σ 2 n -Σ * * -Σ 2 → 0.

Proof of Theorem 6

Respectively from Lemma 2 and 3, m m -µ converges to 0 in quartic mean and d 2 m -δ 2 converges to 0 in quadratic mean. b2 -β 2 converges to 0 in quadratic mean by Lemma 6, and following the idea of proof of Lemma 10, we have: Then, we can apply the result of the theorem 3.2 from (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], so identically we have that E[ S * m -Σ * 2 ] → 0 and S * m has the same asymptotic expected loss as Σ * , i.e. E S * m -Σ 2 -Σ * -Σ 2 → 0.

Proof of Theorem 7

As d 2 s = d 2 and b 2 s = b 2 m , we trivially have the quadratic convergence to δ 2 and β 2 respectively with Lemma 3 and Theorem 6. m s = n-1 n m, and m -µ converges in quartic mean to 0 with µ bounded, so m s -µ converges in quartic mean to 0. Similarly, a 2 s = n-1 n a 2 m and a 2 m -α 2 converges in quadratic mean with α 2 bounded, so a 2 s -α 2 converges in quadratic mean. Finally, as S * s = n-1 n S * m , we have:

E[ S * s -Σ * 2 ] = (n -1) 2 n 2 E S * m -Σ * - 1 n -1 Σ * 2 = (n -1) 2 n 2 E S * m -Σ * 2 - 2 n -1 S * m -Σ * , Σ * + 1 n 2 Σ * 2 E[ S * s -Σ * 2 ] ≤ (n -1) 2 n 2 E S * m -Σ * 2 + 1 n 2 Σ * 2 + 2(n -1) n 2 E[ S * m -Σ * 2 ]E[ Σ * 2 ]. (125)
Using Theorem 6 and that E[ Σ * 2 ] bounded, we have that E[ S * s -Σ * 2 ] → 0 and S * s has the same asymptotic expected loss as Σ * , i.e. E S * s -Σ 2 -Σ * -Σ 2 → 0.

Proof of Lemma 12

Let X •,k ∼ N (0, Σ), k ∈ 1, n , n iid samples. Trivially, µ = 1 p tr(Σ), α 2 = Σ -µI 2 . For β 2 , we will use the equation from the proof of Lemma 1.

β 2 = 1 pn i,j V[y i1 y j1 ] + p + 1 n(n -1) µ 2 + 1 n(n -1) α 2 . ( 126
)
As X is Gaussian, we have for all (i, j) ∈ 1, p

β 2 = p + 1 n -1 µ 2 + 1 n -1 α 2 . (128)
And, of course, δ 2 = α 2 + β 2 .

Proof of Lemma 13

Let ν > 8, Σ a covariance matrix, X •,k ∼ t ν (0, Σ), k ∈ 1, n , n iid samples, with scale matrix Σ = ν-2 ν Σ. With this setup, we have as expected: V[X] = Σ. Obviously, µ = 1 p tr(Σ), α 2 = Σ -µI 2 . For β 2 , we will use the equation from the proof of Lemma 1, as in the Gaussian case.

Figure 1 :

 1 Figure 1: With p the dimension, n the number of samples, the axis z shows the difference of the log 10 of the expected losses between LW m and LW u. The losses are relative to the theoretical bound LW op. The samples are drawn with a Gaussian (left)/t 8.5 (right) distribution, Σ = I p .

Figure 2 :

 2 Figure 2: oss comparison to LW op, Gaussian (left)/t 8.5 (right) distribution, Σ = I p , c = 1.Here, note that LW ex and LW op have both null loss, we didn't plot them.

Figure 3 :

 3 Figure 3: With p the dimension, n the number of samples, the axis z shows the difference of the log 10 of the expected losses between LW m (left)/LW s (right) and LW u. The losses are relative to the theoretical bound LW op. Samples are drawn from a Gaussian distribution, random Σ. The black contour is the iso-line at level 0, where the expected losses are equal.

Figure 4 :

 4 Figure 4: With p the dimension, n the number of samples, the axis z shows the difference of the log 10 of the expected losses between LW m (left)/LW s (right) and LW u. The losses are relative to the theoretical bound LW op. Samples are drawn from a t 8.5 -distribution, random Σ. The black contour is the iso-line at level 0, where the expected losses are equal.

Figure 5 : 1 Figure 6 :

 516 Figure 5: Loss comparison to LW op, Gaussian (left)/t 8.5 (right) distribution, random Σ, c = 1

  (x •k ) t -S 2 .(119)Let's show that b2 r -β 2 converges to 0 in quadratic mean. We use the following decomposition: -1) 2 S 2 .

- 2 m -β 2 ) 2 ]• 0 ≤ a 2 m ≤ d 2 s and a 2 m -α 2

 222222 max(| b2 -β 2 |, |d 2 m -δ 2 |) ≤ b 2 m -β 2 ≤ max(| b2 -β 2 |, |d 2 m -δ 2 |). ≤ E[max(| b2 -β 2 |, |d 2 m -δ 2 |) 2 ] ≤ E[( b2 -β 2 ) 2 ] + E[(d 2 m -δ 2 ) 2 ].(124)Which, as previously, leads to the 2 following results:E[(b 2 m -β 2 ) 2 ] → 0, E[(a 2 m -α 2 ) 2 ] → 0.Finally, we check that we have the set of hypotheses required by the proof of Theorem 2 to work:• α 2 , β 2 , δ 2 are non-negative, bounded, andα 2 + β 2 = δ 2 ,• m m -µ converges to 0 in quartic mean,• d 2 m is nonnegative, and d 2 m -δ 2 -→

  k 4 are all different, then the expectation in the sum equals 0. So there are at most 28n 7 + O(n 6 ) non-zero terms in the sum.

Now, let's find a bound of those expectations. Let's first note that, for all (k, k ) ∈ 1, n 2 and i ∈ 1, p :

  is impossible due to the constraints on the indices in the sum. In the first case, we have: Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = E y ik1 y ik 1 y jk1 y jk 2 y lk3 y lk 3 y mk3 y mk 4= E y ik1 y jk1 y lk3 y mk3 y ik 1 y lk 3 y jk 2 y mk 4 = E y ik1 y jk1 y lk3 y mk3 y ik 1 y lk 1 y jk 2 y mk 2 = E y ik1 y jk1 ]E[ y lk3 y mk3 ]E[ y ik 1 y lk 1 ]E[ y jk 2 y mk 2Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = λ ij λ lm λ il λ jm . ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 ) = E y ik1 y ik 1 y jk1 y jk 2 y lk3 y lk 3 y mk3 y mk 4 = E y ik1 y jk1 y lk3 y mk3 y ik 1 y mk 4 y jk 2 y lk 3 = E y ik1 y jk1 y lk3 y mk3 y ik 1 y mk 1 y jk 2 y lk 2 = E y ik1 y jk1 ]E[ y lk3 y mk3 ]E[ y ik 1 y mk 1 ]E[ y jk 2 y lk 2

	(41)
	And in the second case, we have:
	Cov(y

  |Cov(y ik1 y ik 1 y jk1 y jk 2 , y lk3 y lk 3 y mk3 y mk 4 )|

	And we have:
	1
	i,j,l,m p 4 ≤
	i,j,l,m

2

+ n) different 6-uples of indices k that have G as associated graph. Moreover, there are a finite number N ∈ N (independent of n) of graphs with 6 vertices that have at most 3 connected components. So, there are at most N n((n -1) 2 + n) different 6-uples of indices k such that the associated graph G has 3 edges or more.

  n 8 . Then, using multiple Cauchy-Schwarz and Jensen inequalities, we have:Cov y ik1 y jk 1 -λ ij 1 k1=k 1 y ik2 y jk 2 -λ ij 1 k2=k 2 , y lk3 y mk 3 -λ lm 1 k3=k 3 y lk4 y mk 4 -λ lm 1 k4=k 4 = Cov y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 -λ ij 1 k1=k 1 Cov y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 -λ ij 1 k2=k 2 Cov y ik1 y jk 1 , y lk3 y mk 3 y lk4 y mk 4 -λ lm 1 k3=k 3 Cov y ik1 y jk 1 y ik2 y jk 2 , y lk4 y mk 4 -λ lm 1 k4=k 4 Cov y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 + λ ij λ lm 1 k1=k 1 1 k3=k 3 Cov y ik2 y jk 2 , y lk4 y mk 4 + λ ij λ lm 1 k1=k 1 1 k4=k 4 Cov y ik2 y jk 2 , y lk3 y mk 3 + λ ij λ lm 1 k2=k 2 1 k3=k 3 Cov y ik1 y jk 1 , y lk4 y mk 4 + λ ij λ lm 1 k2=k 2 1 k4=k 4 Cov y ik1 y jk 1 , y lk3 y mk 3 = Cov y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 -E[y ik1 y jk 1 ]Cov y ik2 y jk 2 , y lk3 y mk 3 y lk4 y mk 4 -E[y ik2 y jk 2 ]Cov y ik1 y jk 1 , y lk3 y mk 3 y lk4 y mk 4 -E[y lk3 y mk 3 ]Cov y ik1 y jk 1 y ik2 y jk 2 , y lk4 y mk 4 -E[y lk4 y mk 4 ]Cov y ik1 y jk 1 y ik2 y jk 2 , y lk3 y mk 3 + E[y ik1 y jk 1 ]E[y lk3 y mk 3 ]Cov y ik2 y jk 2 ,y lk4 y mk 4 + E[y ik1 y jk 1 ]E[y lk4 y mk 4 ]Cov y ik2 y jk 2 , y lk3 y mk 3 + E[y ik2 y jk 2 ]E[y lk3 y mk 3 ]Cov y ik1 y jk 1 , y lk4 y mk 4 + E[y ik2 y jk 2 ]E[y lk4 y mk 4 ]Cov y ik1 y jk 1 , y lk3 y mk 3 , k 1 , k 2 , k 2 , k 3 , k 3 , k 4 , k 4 ) ∈ 1, n 8 , we have,

	≤ 9 4 E[y 8 i1 ]E[y 8 j1 ]E[y 8 l1 ]E[y 8 m1 ].
	(65)
	So, for all (k 1 i,j,l,m

  2 ,V[y i1 y j1 ] = λ ii λ jj + λ 2 ij . So,

	1 pn	i,j	V[y i1 y j1 ] =	1 n	Σ 2 +	1 p	tr(Σ) 2	=	1 n	(α 2 + (p + 1)µ 2 ).	(127)
	Which finally leads to,
	β 2 =	p + 1 n(n -1)	µ 2 +	1 n(n -1)	α 2 +	1 n	(α 2 + (p + 1)µ 2 )

Acknowledgments

We want to thank Gabriel Turinici for its advices all along the work, and Olivier Ledoit and Michael Wolf for the impressive work they have done in the field of covariance estimation for twenty years.

Funding

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

We compute the expectations separately for simplicity.

ik1 -λ ii y jk2 ȳj -

Injecting those expectations in the main equation, we obtain:

From the proof Lemma 1, we have:

So, combining those two equations, we obtain:

which concludes the proof.

Proof of Lemma 8

From the proof of Lemma 5, we have:

And from Lemma 7, we have:

Which immediately finishes the proof:

6.12 Proof of Lemma 9

From Lemma 5:

And from Lemma 8:

So,

Then, we deduce:

Concerning the quadratic mean convergence, we use the fact that the variances of m 2 , d 2 and b2 converge to 0 as n goes to infinity, from Corollary 2, and Lemma 3 and 6 respectively. So,

(106)

Proof of Lemma 10

For the upper bound:

For the lower bound:

So, from Lemma 3 and 9,

Which leads trivially to:

Proof of Theorem 2

We will use the proof of Theorem 3.2 in (Ledoit and Wolf, 2004) [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] to prove ours. We check that we have the set of hypotheses required by the proof to work:

• α 2 , β 2 , δ 2 are non-negative, bounded, and

• m -µ converges to 0 in quartic mean,

• d 2 u is nonnegative, and

, with b 2 u ≥ 0. Then, we can apply the result of the theorem 3.2 from [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], so

Proof of Lemma 11

The proof from [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] can be applied here, as the hypotheses of their Lemma A.1 [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] are verified for

, from the same arguments they used.

Proof of Theorem 3

Solving the convex minimization problem, we obtain easily that Σ * * = µI + α2 d 2 (S -mI), with α 2 = S, Σ -mµ, is a minimizer. We have then:

(m -µ) 2 converges to 0 in quadratic mean by Lemma 2. For the second term, we will use the Lemma A.1 from Ledoit and Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] with u 2 = (a 2 u -α 2 ) 2 , τ 1 = 2 and τ 2 = 0. In the following, we check the assumptions of the Lemma A.1, i.e.

It comes that:

To easily compute the variance term, we will use a characterization of multivariate t-distributions.

In fact, as for all k ∈ 1, n , X •,k ∼ t ν (0, Σ), there exists two independant random variables U k and Z •,k such that:

We deduce then,

We have:

So,

And,

We can conclude,

And, of course, δ 2 = α 2 + β 2 .


Nicolas Ragot
email: nicolas.ragot@esigelec.fr

Benchmark of Visual SLAM Algorithms: ORB-SLAM2 vs RTAB-Map* *

Keywords: Visual SLAM, ORB SLAM, RTAB SLAM, Localization, Mapping, Visual Odometry

. Different motion scenarios are used to test and benchmark the SLAM algorithms in various configurations: straight-line, straight-line and back, circular path with loop closure, etc.

I. INTRODUCTION

Estimating the camera trajectory while simultaneously reconstructing the environment is a key and well-know problem in robotics and computer vision. This can be achieved by the use of different measuring means (LIDAR, RADAR) and algorithms (particle filter, extended Kalman filter, Graph-SLAM). When cameras are used, SLAM techniques are named as visual SLAM (vSLAM). vSLAM algorithms have been widely proposed in the field of computer vision, robotics. Taketomi & al. [6] proposes an exhaustive survey of real-time This work is carried out as part of the INTERREG VA FMA ADAPT project "Assistive Devices for empowering disAbled People through robotic Technologies" http://adapt-project.com/index.php. The Interreg FCE Programme is a European Territorial Cooperation programme that aims to fund high quality cooperation projects in the Channel border region between France and England. The Programme is funded by the European Regional Development Fund (ERDF) vSLAM algorithms developed from 2010 to 2016. In [START_REF] Ouerghi | Embedded Architecture for Vision-based Mobile Localization[END_REF], a very exhaustive state of the art about the vSLAMis presented. Obtaining a complete vSLAM system essentially consists in adding the cartography or the geometric consistency of the visual odometry map. The pioneering approach of visual odometry was presented by Nister et al. [START_REF] Ouerghi | Embedded Architecture for Vision-based Mobile Localization[END_REF] in 2003 for monocular and stereo cases. They presented estimation within RANSAC for the first time as well as the estimation from 2D-3D matches instead of triangulating for each stereo pair which allowed to have more robustness [START_REF] Scaramuzza | Visual Odometry [Tutorial[END_REF]. Indeed, the majority of the methods follow the indirect scheme where they carry out features detection and movement estimation within RANSAC. However, direct methods have also been implemented such as SVO [START_REF] Forster | SVO: fast semi-direct monocular visual odometry[END_REF] and DSO [START_REF] Engel | Direct sparse odometry[END_REF].

Location-based on existing maps could also be done without the on-line mapping process or building a map beforehand. In recent years, specific card formats have been proposed to meet the need for prior knowledge for the autonomous car such as OpenStreetMap (OSM). Guo et al. proposed in [START_REF] Guo | Automatic lane-level map generation for advanced driver assistance systems using lowcost sensors[END_REF] a low-cost approach based on OSM maps, INS, GPS and spelling images generated from a camera. A second optimization using visual odometry is performed [START_REF] Guo | A Low-Cost Solution for Automatic Lane-Level Map Generation Using Conventional In-Car Sensors[END_REF].

Location-based on recognition of places: usually done in two phases: an off-line learning phase where a geotagged database is built and the location is done online by finding the image that most closely resembles the current one. Site recognition is also an essential part of a visual SLAM system to establish loop closures on one hand and to relocate on the other hand in the event of loss of tracking. The visual points of interest are further subdivided into two categories: those that select parts of the image and those that use the entire image without resorting to a selection step. Examples of the first category are the SIFT [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF] and SURF [START_REF] Bay | Surf: Speeded up robust features[END_REF] descriptors, and the second category global descriptors such as HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] and GIST [START_REF] Oliva | Building the gist of a scene: the role of global image features in recognition[END_REF]. bag-of-words separates the point of interest space into a finite number of visual words. FabMap [START_REF] Cummins | FAB-MAP: Probabilistic Localization and Mapping in the Space of Appearance[END_REF] [START_REF] Cummins | Appearance-only SLAM at large scale with FAB-MAP 2.0[END_REF] which is a probabilistic approach that has shown its effectiveness on a trajectory of 1000 km in length. Other work has been based on global descriptors such as HOG in [START_REF] Mcmanus | Scene Signatures: Localised and Pointless Features for Localisation[END_REF] and GIST in [START_REF] Siagian | Biologically Inspired Mobile Robot Vision Localization[END_REF]. Global descriptors are more dependent on the point of view than local descriptors, but they do not work well with changes in brightness. The most representative approach in this category is SeqSLAM [START_REF] Milford | SeqSLAM: visual route-based navigation for sunny summer days and stormy winter nights[END_REF] which has demonstrated a much better performance than FabMAP under conditions of extreme atmospheric changes and brightness. SeqSLAM is based on the use of subsampled whole images by searching in a database that is built in the sequence that most closely resembles the current sequence acquired. An extension of SeqSLAM for 2D interior environments has been proposed in [START_REF] Mount | 2D visual place recognition for domestic service robots at night[END_REF].

We therefore present in this paper a benchmark comparative study between two visual SLAM approaches: ORB SLAM and RTAB SLAM. The ORB SLAM has been implemented taking into account different kinds of sensors such as monocular, stereo and RGB-D camera. The RTAB SLAM, meanwhile, it has been experimented taking into account both monocular and RGB-D images. The wheelchair has been modified to fit the developed platform. A new motor controller has been developed, the wheelchair is connected within a wireless communication system, and the instrumentation is carried out within a Realsense camera including a monocular, stereo and a RGB-D sensors. A ground truth data based VICON system is used to validate the results obtained.

This paper is organized as follows: Section 1 introduces the motivation of the paper. Section 2 presents the state of the art about SLAM. Section 3 presents the wheelchair architecture and how to control the wheelchair. Section 4 illustrates the comparison between ORB SLAM and RTAB SLAM within different kind of camera like monocular, stereo and RGB-D camera. In the same section, we present how we have validated the approach developed by comparison between SLAM localization and ground truth localization carried out by VICON system. The results obtained will be presented in Section 5. Finally, in Section 6, we will conclude this paper.

II. RELATED WORKS

The literature review suggests to classify the methods in two main categories: one is related to the indirect methods whereas the other deals with the direct ones. The indirect methods are based on the detection of points of interest in the image (features) and the establishment of 2D-3D correspondences between the points of interest and the map by matching/tracking. The pose of the camera is then obtained by solving the Perspective-n-Point (PnP) problem. It is noted that the majority of vSLAM methods assume that the camera is previously calibrated, which implies that the intrinsic parameters of the camera are previously known [START_REF] Ouerghi | Mobile Robot Localization State of the Art[END_REF].

In robotic mapping and navigation, Simultaneous Localization and Mapping (SLAM) is the computational problem of constructing or updating a map of an unknown environment while simultaneously keeping track of an agent's location within it. While this initially appears to be a chicken-and-egg problem there are several algorithms known for solving it, at least approximately, in tractable time for certain environments. Popular approximate solution methods include the particle filter, extended Kalman filter, and GraphSLAM.

Indirect methods-based feature detection belong to other subcategories of vSLAM namely filter-based SLAM and optimization-based SLAM. A pioneering method is MonoSLAM presented in 2007, allowing to simultaneously estimate the pose of the camera and the 3D structure of an unknown environment using an Extended Kalman filter (EKF) [START_REF] Davison | MonoSLAM: Real-Time Single Camera SLAM[END_REF]. The 3D camera movement (6 DOF) and the 3D positions of the features are represented by a state vector of EKF. PTAM offers an improvement by separating tracking and mapping in different threads of the CPU [START_REF] Klein | Parallel Tracking and Mapping on a Camera Phone[END_REF]. PTAM also introduced the use of keyframes for mapping as well as a feature classifier to search for the keyframe closest to the input image. ORB-SLAM can be considered as an extension of PTAM including features such as separate threads and keyframes, the detection of closing loops and the optimization of the poses graph [START_REF] Mur-Artal | ORB-SLAM: A Versatile and Accurate Monocular SLAM System[END_REF]. ORB-SLAM [START_REF] Mur-Artal | ORB-SLAM: A Versatile and Accurate Monocular SLAM System[END_REF], introduced in 2015, is the most complete indirect monocular vSLAM method and has been extended for stereo systems and RGB-D systems in ORB-SLAM2 in 2017 [START_REF] Mur-Artal | Orb-slam2: An open-source slam system for monocular, stereo, and rgb-d cameras[END_REF].

Unlike feature-based indirect methods, use directly pixel intensity of input images. A completely direct method is DTAM which, like PTAM, has two distinct parts: pose tracking and mapping [START_REF] Newcombe | DTAM: Dense tracking and mapping in real-time[END_REF]. However, PTAM tracks a set of 3D points while DTAM maintains dense depth maps for a selection of keyframes. The most representative method in the category of direct methods is LSD-SLAM (Large-Scale Direct SLAM) which is a monocular SLAM algorithm [START_REF] Engel | LSD-SLAM: Large-Scale Direct Monocular SLAM[END_REF]. The map is constructed based on keyframes composed of camera image, a map of the inverse depths (inverse depth map) and the variance of the inverse of depths.

III. ROBOTICS WHEELCHAIR ARCHITECTURE

The robotics electrical powered wheelchair from the ESIG-ELEC's lab is an Invacare, model Bora from which all the proprietary electronics has been removed and replaced by:

• Embedded PC running a Linux Ubuntu 16.04 LTS distribution • Motor driver from Roboteq • Xbox joystick in place of the original one which can be connected to the wheelchair through a USB or Bluetooth connection.

• WIFI router for getting a wireless local network on the wheelchair • Embedded HMI with a touch-screen panel The wheelchair software is fully developed under the ROS robotics middleware. Figure 1 shows the wheelchair description. An embedded computer is used as a heart of processor in the wheelchair. It possesses 8 GB of RAM and 250 GB of SSD. It runs Ubuntu 16.04 as operating system. A tactile screen is connected to the embedded computer using the HDMI cable to present the user interface to the user. A motor driver from Roboteq is used to drive motors. The wheelchair also has its own Wi-Fi router. An Xbox One joystick is used to control the wheelchair. A RealSense D435 depth camera from Intel is used for embedded vision. Figure 2 shows the hardware architecture of the wheelchair.

IV. VSLAM COMPARISON BETWEEN ORB-SLAM AND RTABMAP SLAM A. ORB-SLAM

ORB-SLAM2 is a real-time SLAM library for monocular, stereo and RGB-D cameras that computes the camera trajec-tory and a sparse 3D reconstruction. It is able to detect loops and re-localize the camera in real time. The system works in real-time on standard CPUs in a wide variety of environments from small hand-held indoors sequences, to drones ying in industrial environments and cars driving around a city. The back-end based on bundle adjustment with monocular and stereo observations allows for accurate trajectory estimation with metric scale. The system includes a lightweight localization mode that leverages visual odometry tracks for unmapped regions and matches to map points that allow for zero-drift localization. The main functionalities of ORB SLAM are: feature tracking, mapping, loop closure and localization.

1) Monocular ORB-SLAM 2: Monocular SLAM requires just a single camera. The ORB SLAM converts the image to grayscale for its application. The first step is to detect features and initialize the map and its position. Once it gets initialized, it starts creating map. Monocular SLAM requires a procedure to create an initial map because depth cannot be recovered from a single image. One way to solve the problem is to initially track a known structure. In the context of ltering approaches, points can be initialized with high uncertainty in depth using an inverse depth parametrization, which hopefully will later converge to their real positions. The figure below shows the flow of process in ORB-SLAM monocular.

2) Stereo ORB-SLAM 2: Stereo ORB-SLAM requires a stereo camera. The feature detection in stereo SLAM is better than monocular SLAM. The initialization is fast and has many advantages over monocular SLAM. For stereo cameras, we extract ORB in both images and for every left ORB we search for a match in the right image. This can be done very efficiently assuming stereo rectified images, so that epipolar lines are horizontal. We then generate the stereo keypoint with the coordinates of the left ORB and the horizontal coordinate of the right match, which is subpixel refined by patch correlation.

3) RGB-Depth ORB SLAM 2: RGB-D SLAM requires an RGB image with its depth image. For RGB-D cameras, we extract ORB features on the RGB image, for each feature with coordinates (uL; vL) we transform its depth value d into a virtual right coordinate:

uR = uL - f x * b d (1)
where f x is the horizontal focal length and b is the baseline between the structured light projector and the infrared camera.

The uncertainty of the depth sensor is represented by the uncertainty of the virtual right coordinate. In this way, features from stereo and RGB-D input are handled equally by the rest of the system. 4) MAP Save and Load: A feature to save the map generated by the ORB SLAM is developed. Ability to load the map for localization purpose was also developed. All keyframes, features and its descriptors are saved to a binary file. After every mapping session the algorithm asks the user whether to save the map or not. At the beginning of SLAM, the algorithm asks the user whether to load the saved map to begin localization or mapping within the environment which is already built. RTAB-Map supports both visual and lidar SLAM, providing in one package a tool allowing users to implement and compare a variety of 3D and 2D solutions for a wide range of applications with dierent robots and sensors. It uses depth image with RGB images to construct maps. The graph is created here, where each node contains RBG and depth images with corresponding odometry pose. The links are transformation between each node. When the graph is updated, RTAB-Map compares the new image with all previous ones in the graph to find a loop closure. When a loop closure is found, graph optimization is done to correct the poses in the graph. For each node in the graph, we generate a point cloud from the RGB and depth images. This point cloud is transformed using the pose in the node. The 3D map is then created [START_REF] Scaramuzza | Visual Odometry [Tutorial[END_REF].

V. RESULTS ANALYSIS

The tests were carried out in LNA lab of IRSEEM for indoor scenarios and in car parking of IRSEEM for outdoor scenarios. LNA lab is equipped with Vicon camera system for tracking the wheelchair position and will be used as ground truth to evaluate the distance measurement of SLAM.

A. ORB SLAM

ORB SLAM works perfectly except with monocular camera as it loses points during mapping and also initializes very slowly. It can estimate its odometry and track features efficiently. The ORB SLAM creates a sparse 3D map as shown in the Figure 3.

The ORB SLAM creates a sparse 3D map. ORB SLAM is successful in creating a map and generating its trajectory in the environment. It can also detect closed loops successfully. The Figure 4 shows the trajectory created by ORB SLAM in indoor environments. Figure 5 shows the trajectory in outdoor environments.

The distance estimation is also accurate to the ground truth measurement by the Vicon. The only problem with ORB SLAM with monocular camera was that it was slow at initialization and also it loses points during the mapping. Due to this the ORB SLAM loses its position and will not be able to continue with mapping. The Figure 6 shows an example where a trajectory is generated by ORB SLAM with monocular in comparison with ORB SLAM with stereo and RGBD.

The Table I below shows the distance measurement for various paths. For the indoor scenarios we have used Vicon as the ground truth and for outdoor scenarios we have used measuring tape as the ground truth (approximate distance). The ORB SLAM with monocular camera had trouble initializing and also could not complete the trajectory. This is a limitation with monocular cameras and only alternative is to use either stereo camera or RGB-D camera. From Table I we can observe the results of ORB SLAM with stereo and RGBD. In indoor scenarios, ORB SLAM was accurate with both stereo and RGBD cameras but in outdoor scenarios, stereo camera has shown better results compared to the RGBD. A comparison of trajectories in two indoor scenarios by ORB SLAM with ground truth trajectory by Vicon was done. The trajectories generated by ORB SLAM was closely accurate to that of the ground truth trajectory by Vicon. The Figure 7 shows the trajectory comparison. The Table II shows the RATBMAP SLAM measurement of distance in comparison with Vicon measurement. From above results you can see that the mapping and trajectory estimation of RTABMAP is very good but the distance measurement is not. The distance estimated by the RTABMAP is not better than that of ORB SLAM.

C. ORB SLAM vs RTABMAP

The ORB SLAM was compared with RTABAMP in same scenarios as RTABMAP was tested. The Table III shows the distance measurement comparison between two SLAM algorithms. From the above table we can say that the distance measurements by ORB SLAM are more accurate than that of RTABMAP. From the above table we can say that the distance measurements by ORB SLAM are more accurate than that of RTABMAP.

VI. CONCLUSION

In this paper, we have presented a wheelchair localization and mapping-based comparison of two vSLAM algorithms, and within both of indoor and outdoor environment. Different approaches in VSLAM were developed around ORB SLAM and RTAB SLAM. In order to improve the quality of localization and mapping obtained by the basic SLAM, we have updated the ORB SLAM to our system by taking into account some parameters. We have tested ORB SLAM using Monocular camera, Stereo camera and RGBD camera.

we also added a feature where ORB SLAM prompts the user for saving the map and loading the map for pure localization. The ORB SLAM is able to save map and later load the map for localization purpose. The ORB SLAM has also been improved to estimate an accurate trajectory and also provide accurate odometry estimations. However, ORB SLAM with Monocular camera could not provide good trajectory or distance estimations due to the limitations of Monocular camera. The RTABMAP is also tested and used to compare with ORB SLAM. The trajectory estimation by RTABMAP is accurate but the odometry estimation is not as accurate as ORB SLAM. Few tests show that RTABMAP is better at trajectory estimation when compared to ORB SLAM but when it comes to distance measurement ORB SLAM is better. We have validated the whole developments under wheelchair platform within comparison between SLAM data, RTAB data, and ground truth data-based VICON obtained in indoor environment. Other scenario was carried out in outdoor environment for the validation of the platform.

Fig. 1 .

 1 Fig. 1. ADAPT Real Wheelchair Platform.

Fig. 2 .

 2 Fig. 2. Hardware Architecture of the Wheelchair.

  B. RTABMAP SLAM RTAB-MAP stands for Real Time Appearance Based Mapping. It is distributed as an open source library since 2013. RTAB-Map started as an appearance-based loop closure detection approach with memory management (shown in below figure) to deal with large-scale and long-term online operation. It then grew to implement Simultaneous Localization and Mapping (SLAM) on various robots and mobile platforms.

Fig. 3 .Fig. 5 .

 35 Fig. 3. An example of a map created by ORB SLAM

Fig. 6 .

 6 Fig. 6. ORB SLAM trajectory comparison between Monocular, Stereo and RGBD cameras.

Fig. 7 .

 7 Fig. 7. ORB SLAM and Vicon trajectory comparison.

Fig. 8 .

 8 Fig. 8. Map and trajectory estimation by RTABMAP SLAM.

TABLE I ORB

 I SLAM AND VICON DISTANCE (IN METERS) COMPARISON.

	Scenario	Ground Truth	ORB Stereo ORB RGBD
	Indoor 1	7.82	7.69	7.16
	Indoor 2	6.98	7.38	6.81
	Indoor 3	15.05	15.49	14.32
	Indoor 4	21.75	21.17	21.31
	Indoor 5	29.69	21.47	22.13
	Indoor 6	37.69	38.84	38.79
	Outdoor 1	28	28.85	19.27
	Outdoor 2	28	27.29	17.83
	Outdoor 3	56	56.63	37.25
	Outdoor 4	92	94.55	65.87
	Outdoor 5	92	95.8	66.97
	Outdoor 6	80	84.38	66.24
	Outdoor 7	80	81.9	63.36
	Outdoor 8	160	172.6	129.41

TABLE II RTABMAP

 II DISTANCE COMPARISON WITH VICON.

	Scenario	Ground Truth	RTABMAP RGBD
	Indoor 1	6.48	11.83
	Indoor 2	24.5	37.61
	Indoor 3	23.26	31.31

TABLE III ORB

 III SLAM AND RTABMAP COMPARISON.

	Scenario	Ground Truth	RTAB	ORBSLAM	ORBSLAM
			RGBD	Stereo	RGBD
	Indoor 1	6.48	11.83	5.33	8.94
	Indoor 2	24.5	37.61	30.64	27.78
	Indoor 3	23.26	31.31	20.82	21.43

ACKNOWLEDGMENT

This research is supported by ADAPT Project (cofinanced by the European Regional Development Fund within the framework of the INTERREG VA France (Channel) England programme). Many thanks to the engineers of Autonomous Navigation Laboratory of IRSEEM for their support with VICON system. Many thanks to Technological Resources Center (TRC) of ESIGELEC for his help in testing phase.