domain: math stat info

text:
Random Matrix-Improved Estimation

Random Matrix-Improved Estimation of the Wasserstein Distance between two Centered Gaussian Distributions

Malik Tiomoko 1 , Romain Couillet 1,2, * 1 CentraleSupélec, Université ParisSaclay, 2 GIPSA-lab, Université Grenoble-Alpes Abstract-This article proposes a method to consistently estimate functionals 1 ), when the size p and number na of the (zero mean) samples x (a) i are similar. As a corollary, a consistent estimate of the Wasserstein distance (related to the case f (t) = √ t) between centered Gaussian distributions is derived.

The new estimate is shown to largely outperform the classical sample covariance-based "plug-in" estimator. Based on this finding, a practical application to covariance estimation is then devised which demonstrates potentially significant performance gains with respect to state-of-the-art alternatives.

I. INTRODUCTION

Many machine learning and signal processing applications require an adequate framework to compare statistical objects, starting with probability distributions. The Wasserstein distance, initially inspired by Monge [START_REF] Gaspard Monge | Mémoire sur la théorie des déblais et des remblais[END_REF] and later by Kantorovich [START_REF] Kantorovich | On the translocation of masses[END_REF] in a transport theory analogy, provides a natural notion of dissimilarity for probability measures and finds a wide spectrum of applications in image analysis [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], shape matching [START_REF] Su | Optimal mass transport for shape matching and comparison[END_REF], computer vision [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF], etc.

However, computing the Wasserstein distance is expensive as it requires to minimize a cost function taking the form of an integral over the space of probability measures. Despite recent advances [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], where regularized approximations that reduce this numerical cost are proposed, the latter is still involved in general. Special cases exist for which the Wasserstein distance assumes a closed form, particularly when the underlying distributions are zero-mean Gaussian with covariance matrices C 1 and C 2 . The closed-form formula however involves the eigenvalues of C 1 C 2 and thus depends on the unknown population covariance matrices C 1 and C 2 . Assuming the observation of n 1 , n 2 p samples with covariances C 1 , C 2 , respectively, C 1 C 2 is conventionally approximated by its empirical version Ĉ1 Ĉ2 . As we will show, this induces a dramatic estimation bias in practical applications where p is rather large or, equivalently, n 1 , n 2 rather small, a standard assumption in big data applications.

Based on recent advances in random matrix theory, this article proposes a new consistent estimate for the Wasserstein distance between two centered Gaussian distributions when the dimension p of the samples is of the same order of magnitude as their numbers n 1 , n 2 . This work enters the scope of Mestre's seminal ideas [START_REF] Mestre | On the asymptotic behavior of the sample estimates of eigenvalues and eigenvectors of covariance matrices[END_REF] on the estimation of *Couillet's work is supported by the ANR Project RMT4GRAPH (ANR-14-CE28-0006) and the IDEX GSTATS Chair at University Grenoble Alpes. functionals 1 p p i=1 f (λ i (C)) of the eigenvalue distribution of population covariance matrices C, which can be related to the (limiting) eigenvalue distribution of the sample estimates Ĉ via a complex integration trick. We recently extended this work to the estimation of functionals of the eigenvalue distribution of F-matrices in [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF], i.e., matrices of the form C -1 1 C 2 , and applied to the estimation of the natural geodesic Fisher distance, Battacharrya distance, and Rényi/Kullbach-Leibler divergences between Gaussian distributions.

Our main contribution is the extension of [START_REF] Mestre | On the asymptotic behavior of the sample estimates of eigenvalues and eigenvectors of covariance matrices[END_REF], [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] to functionals f of the eigenvalues of products C 1 C 2 of population covariance matrices. The Wasserstein distance falls within this scope for f (t) = √ t. Unlike [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF], where the functionals of interest (f (t) = t, log(t), log 2 (t)) are amenable to explicit evaluations of the complex integrals, the present f (t) = √ t scenario is more technically involved and gives rise to real non-explicit, yet numerically computable, integrals.

In the remainder of the article, Section II introduces the main model and assumptions, Section III provides our key technical result and its corollary to the Wasserstein distance estimation, and a practical application to covariance matrix estimation is finally proposed in Section IV.

Reproducibility. Matlab codes for the various estimators introduced and studied in this article are available at https://github.com/maliktiomoko/RMTWasserstein II. MODEL AND MAIN OBJECTIVE

For a ∈ {1, 2}, let X a = [x (a) 1 , . . . , x (a)
na ] be n a independent and identically distributed random vectors with x

(a) i = C 1 2 a x(a)
i , where x(a) i ∈ R p has zero mean, unit variance and finite fourth order moment entries. This holds in particular for x (a) i ∼ N (0, C a ). In order to control the growth rates of n 1 , n 2 , p, we make the following assumption:

Assumption 1 (Growth Rates). As n a → ∞, p/n a → c a ∈ (0, 1) and lim sup p max{ C -1 a , C a } < ∞ for • the operator norm.
We define the sample covariance estimate Ĉa of C a as

Ĉa ≡ 1 n a X a X T a = 1 n a na i=1 x (a) i x (a)T i . The Wasserstein distance D W (C 1 , C 2 )
between two zeromean Gaussian distributions with covariances C 1 and C 2 , respectively, assumes the form [9, Remark 2.31]:

D W (C 1 , C 2 ) = tr(C 1 ) + tr(C 2 ) -2tr (C 1 2 1 C 2 C 1 2 1 ) 1 2
.

(1)

arXiv:1903.03447v1 [stat.ML] 8 Mar 2019

It is easily shown that, under Assumption 1,

1 p tr Ĉa - 1 p trC a → 0 almost surely. But estimating tr(C 1 2 1 C 2 C 1 2
1 )

12 is more involved: this is the focus of the article. Up to a normalization by p, this term can be written under the functional form:

1 p tr(C 1 2 1 C 2 C 1 2
1 )

1 2 = 1 p n i=1 λ i (C 1 C 2 ) ≡ D(C 1 , C 2 ; √ •) (2) with λ i (X) the i-th smallest eigenvalue of X.
Our objective is to estimate the more generic form

D(C 1 , C 2 ; f ) ≡ 1 p n i=1 f (λ i (C 1 C 2 )) (3)
for f : R → R a real function admitting a complexanalytic extension. To this end, we shall relate the eigenval-

ues λ i (C 1 C 2 ) to λ i ( Ĉ1 Ĉ2 ) through the Stieltjes transform (m θ (z) ≡ dθ(λ)
λ-z for measure θ and z ∈ C) of their associated normalized counting measures

µ p = 1 p p i=1 δ λi( Ĉ1 Ĉ2) , ν p = 1 p p i=1 δ λi(C1C2) .
In particular,

m µp (z) = 1 p p i=1 1 λi-z for λ i = λ i ( Ĉ1 Ĉ2
). With these notations, we are in position to introduce our main results.

III. MAIN RESULTS

The following theorem provides a consistent estimate for the metric D(C 1 , C 2 ; f ) defined in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF].

Theorem 1. Let Γ ⊂ {z ∈ C, real[z] > 0} be a contour surrounding ∪ ∞ p=1 supp(µ p ). Then, under Assumption 1, D(C 1 , C 2 ; f ) -D(X 1 , X 2 ; f ) a.s. -→ 0 where D(X 1 , X 2 ; f ) = n 2 2πip Γ f ϕ p (z) ψ p (z) ϕ p (z) ϕ p (z) - ψ p (z) ψ p (z) ψ p (z)dz and, recalling m µp (z) = 1 p p i=1 1 λi-z for λ i = λ i ( Ĉ1 Ĉ2 ), ϕ p (z) = z 1-p n 1 -p n 1 zmµ p (z) , ψ p (z) = 1 -p n2 -p n2 zm µp (z). The result of Theorem 1 is very similar to [8, Theorem 1] established for functionals of the eigenvalues of C -1 1 C 2 .
The main difference lies in the expression of the function ϕ p (z).

Proof. The proof of Theorem 1 is based on the same approach as for [10, Theorem 1]. One first creates a link between the Stieltjes transform m νp and D(C 1 , C 2 ; f ) using Cauchy's integral formula:

1 p p i=1 f (λ i (C 1 C 2 )) = f (t)dν p (t) = 1 2πi Γν f (z) z -t dz dν p (t) = -1 2πi Γν f (z)m νp (z)dz (4)
with Γ ν a contour surrounding the support supp(ν p ) of ν p . To relate the unknown m νp to the observable m µp , we proceed as follows. By first conditioning on Ĉ1 ,

Ĉ 1 2 1 Ĉ2 Ĉ 1 2 1
is seen as a sample covariance matrix for the samples

Ĉ 1 2 1 C 1 2 2
x(2) i , for which [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF] allows one to relate m µp to the Stieltjes transform of the eigenvalue distribution

ζ p of C 1 2 2 Ĉ1 C 1 2
2 . The latter is yet another sample covariance matrix for the samples C x [START_REF] Gaspard Monge | Mémoire sur la théorie des déblais et des remblais[END_REF] i ; exploiting [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF] again creates the connection from m ζp to m νp . This entails the two equations:

m µp (z) = ϕ p (z)m ζp (ϕ p (z)) + o p (1) (5)
m νp z Ψ p (z) = m ζp (z)Ψ p (z) + o p (1). ( 6
)
where Ψ p (z) ≡ 1 -p n2 -p n2 zm ζp (z). Successively plugging ( 5)-( 6) into ( 4) by means of two successive appropriate changes of variables, we obtain Theorem 1.

Theorem 1 takes the form of a complex integral which, for generic choices of f , needs be numerically evaluated. In the specific case of present interest where f (z) = √ z, this complex integral can be evaluated as follows.

Theorem 2. Let λ 1 ≤ . . . ≤ λ p , with λ i ≡ λ i ( Ĉ1 Ĉ2 ), and define {ξ i } p i=1 and {η i } p i=1 the (increasing) eigenvalues of Λ -

1 n1 √ λ √ λ T and Λ -1 n2 √ λ √ λ T
, respectively, where λ = (λ 1 , . . . , λ p )

T , Λ = diag(λ) and √ . is understood entry wise. Then, under Assumption 1,

D(C 1 , C 2 ; √ •) -D(X 1 , X 2 ; √ •) a.s.
-→ 0

where, if

n 1 = n 2 , D(X 1 , X 2 ; √ •) = 2 √ n 1 n 2 1 p p j=1 λ j + 2n 2 πp p j=1 ηj ξj - ϕ p (x) ψ p (x) ψ p (x)dx
with ϕ p , ψ p defined in Theorem 1 and, if

n 1 = n 2 , D(X 1 , X 2 ; √ •) = 2n 1 p p j=1 λ j -ξ j .
While still assuming an integral form (when n 1 = n 2 ), this formulation no longer requires the arbitrary choice of a contour Γ and significantly reduces the computational time to estimate D(C 1 , C 2 , √ •). For n 1 = n 2 , a case of utmost practical interest, the expression is completely explicit and computationally only requires to evaluate the eigenvalues ξ j of Λ-

1 n1 √ λ √ λ T .
The latter being a (negative definite) rank- Proof. The ξ i and η i , as defined in the theorem statement, are the respective zeros of the rational functions 1 -p n1p n1 zm μp (z) and 1-p n2 -p n2 zm μp (z) (see [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF]Appendix B]). Thus, ϕ p and ψ p can be expressed under the rational form:

ϕ p (z) = z p i=1 z -λ i p i=1 z -η i , ψ p (z) = p i=1 z -ξ i p i=1 z -λ i .
Evaluating the estimate from Theorem 1 for f (z) = √ z then requires to evaluate a complex integral involving rational functions and square roots of rational functions. Since the complex square root is multivalued, a careful control of "branch-cuts" is required. To perform this calculus, we deform the integration contour Γ of Theorem 1 into Γ as per Figure 1. In the case n 1 = n 2 , the closed null-integral contour Γ (blue in Figure 1) is the sum of the sought-for integral over Γ and of four extra components: 1) Integrals over -radius circles around ξ i : those are null in the limit → 0, as confirmed by a change of variable z = ξ i + e ıθ which allows one to bound the integrand; 2) Integrals over the real axis (in the → 0 limit):

A 2 = n 2 πp p j=1 ηj - ξj + -(ϕ p ψ p )(x) 2 ψ p (z) ψ p (z) - ϕ p (z) ϕ p (z) + ψ p (z) ψ p (z) dx = 2n 2 πp p j=1 ηj ξj -ϕ p ψ p (x) ψ p (z) ψ p (z) dx - n 2 πp p j=1 ηj - ξj + -ϕ p ψ p (x) ϕ p ψ p (x) d dx (ϕ p (x)ψ p (x)) dx = 2n 2 πp p j=1 ηj ξj - ϕ p (x) ψ p (x) ψ p (x)dx -2 n 2 πp p j=1 1 d dx 1 (ϕpψp(x)) (η j ) + o( )
3) Integrals over the -radius circles around η j , with → 0

A 3 = 2 n 2 πp p j=1 1 d dx 1 (ϕpψp(x)) (η j ) + o( )
which thus compensates the last ( -diverging) term in A 2 . 4) Residues in the λ j poles

A 4 = 2 n 2 p lim z→λj p j=1 (ϕ p ψ p )(z) = 2 n 2 p n 1 n 2 p j=1 λ j .
Putting these terms together entails the result of the theorem for the case where n 1 = n 2 . For n 1 = n 2 , it suffices to take the limit of the expression as ξ j → η j . This yields:

D(X 1 , X 2 ; √ •) = 2n 1 p p j=1 λ j + 2n 1 p p j=1 1 π lim t→ξj t ξj - ϕ p (x) ψ p (x) ψ p (x)dx = 2n 1 p p j=1 λ j - 2n 1 p p j=1 1 2πı lim →0 Γ ξ j -ϕ p ψ p (x) ψ p (x) ψ p (x) dx
where Γ ξj is an -radius circular contour around ξ j . The second equality is obtained by deforming the real integral in the complex plane (see [START_REF] Saff | Fundamentals of complex analysis with applications to engineering and science[END_REF] for complex analysis details).

The result unfolds by letting z = ξ i + e ıθ .

Consequently, we obtain the following n, p-consistent estimate for the Wasserstein distance D W (C 1 , C 2 ) of (1).

Corollary 1 (Consistent Estimate of D W (C 1 , C 2 )). Under Assumption 1, 1 p D W (C 1 , C 2 ) - 1 p tr( Ĉ1 + Ĉ2 ) -2 D(X 1 , X 2 ; √ •) a.s.
-→ 0

for D(X 1 , X 2 ; √ •) given by Theorem 2.

Remark 1

(Estimation of C 1 -C 2 2
F
). The Frobenius distance between two covariance matrices also falls under the scope of the present article for the function f (z) = z. Indeed,

D F (C 1 , C 2 ) = C 1 -C 2 2 F = tr C 2 1 + C 2 2 -2tr (C 1 C 2 )
. Then under Assumption 1 and along with the fact that 1 p trC 2 1 can be estimated consistently from

1 p tr Ĉ2 1 -1 n1p (tr Ĉ1 ) 2 , 1 p D F (C 1 , C 2 )- 1 p tr( Ĉ2 1 + Ĉ2 2 ) - p n 1 1 p tr Ĉ1 2 - p n 2 1 p tr Ĉ2 2 -2 D(X 1 , X 2 ; •) a.s.
-→ 0.

In this case, D(X 1 , X 2 ; •) assumes the simple expression

D(X 1 , X 2 ; •) = 1 p p j=1 λ j = 1 p tr Ĉ1 Ĉ2 which follows from 1 p tr Ĉ1 Ĉ2 -1 p trC 1 C 2 a.s.
-→ 0 (by elementary probability arguments) or equivalently from a residue calculus based on Theorem 1 for f (z) = z.

IV. SIMULATIONS AND APPLICATIONS

In this section, we first corroborate our theoretical findings by comparing the classical plug-in estimator to our proposed estimator on synthetic Gaussian data. We then provide an application of our results to improved covariance matrix estimation based on few samples.

A. Confirmation of our results on synthetic data

We here compare the classical plug-in estimate of the Wasserstein distance (that is (1) with C a replaced by Ĉa , a = 1, 2) with our proposed estimate in Corollary 1. Table I lists the results obtained for Toeplitz matrices C 1 , C 2 estimated based on various values of p, n 1 , n 2 . While our proposed estimator is designed under a large p, n 1 , n 2 assumption (as per Assumption 1), it achieves competitive performances even for small values of p, corroborating here our findings in [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] for other classes of covariance matrix distances.

p D W (C 1 , C 2 )
WITH [C 1 ] ij = .2 |i-j| , [C 2 ] ij = .4 |i-j| , x

B. Application to covariance matrix estimation

As a concrete application, Theorem 1 may be used to improve the actual estimation of covariance matrices under a small number n ∼ p of sample data, as similarly performed in [START_REF] Tiomoko | Random matrix improved covariance estimation for a large class of metrics[END_REF] for other covariance matrix distances.

The idea is as follows: we first particularize Theorem 1 and Theorem 2 to the case where one of the covariance matrices, say C 1 , is known by taking c 1 = 0 (i.e., n 1 → ∞ for all fixed p). This gives access to estimates for D W (M, C 2 ; √ •) for all deterministic positive definite matrix M . We then minimize this estimated distance over M in order to estimate C 2 by means of a gradient descent approach.

For C 1 known, we redefine µ p = 1 p p i=1 δ λi(C1 Ĉ2) and obtain, as a corollary of Theorem 1:

Theorem 3. Let Γ ⊂ {z ∈ C, real[z] > 0} a contour surrounding ∪ ∞ p=1 supp(µ p ). Then, D(C 1 , C 2 ; f ) - 1 2πic 2 Γ F -m μp (z) dz a.s. -→ 0 with m μp (z) = p n2 m µp (z) + p-n2 n2z and F (z) = f ( 1 z ).
Proof. For C 1 known (c 1 → 0), ϕ p (z) = z, and the estimator of Theorem 1 yields:

D(X 1 , X 2 ; f ) = 1 2πi Γ f z ψ p (z) ψ p (z) ψ p (z) - 1 z ψ p (z)dz c 2 .
Using the relation m μp (z) = -ψp(z) z , we then get

D(X 1 , X 2 ; f ) = - 1 2πic 2 Γ f - 1 m μp (z) m μp (z)zdz
and the result is immediate after an integration by parts.

For f (z) = √ z, one has F (z) = 2 √ z and we obtain, with a similar proof as for Theorem 2,

D(C 1 , C 2 ; √ •) -D(C 1 , X 2 ; √ •) a.s. -→ 0, D(C 1 , X 2 ; √ •) = 2 πc 2 p j=1 λj ξj √ m μp(x) dx.
Our objective is now to exploit the fact that

C 2 = argmin M 0 D W (M, C 2 ) (8)
where the minimization is over the open cone of positive definite matrices. Using the approximation D(M, C 2 ; √ •) D(M, X 2 ; √ •), we are then tempted to minimize

1 p tr(M + Ĉ2 ) -2 D(M, X 2 ; √ •) in place of D W (M, C 2 ).
The former quantity however has a non zero probability to be negative, and we thus instead propose to estimate C 2 as:

Č2 = argmin M h(M ) h(M ) = 1 p tr(M + Ĉ2 ) -2 D(M, X 2 ; √ •) 2 .
To compute the gradient ∇h(M ) of h at position M , one needs to evaluate the differential Dh(M )[ξ], at M and in the direction ξ, in the Riemmanian manifold of p × p symmetric positive definite matrices (see [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], [START_REF] Tiomoko | Random matrix improved covariance estimation for a large class of metrics[END_REF] to further technical details). We then use the relation Dh(M )[ξ] = ∇h(M ), ξ M where •, • . is the Riemmanian metric defined as η, ξ M = tr M -1 ηM -1 ξ . We obtain the relation

πıp ∇h(M ) 2 h(M ) = 1 p M 2 + p j=1 λj ξj 1 m μp (x) sym M Ĉ2 (M Ĉ2 -xI p ) -2 M dx
where sym(A) = 1 2 (A + A T ) is the symmetric part of A ∈ R p×p . We can write the latter as:

∇h(M ) = 2 h(M ) sym V Λ ∇ V -1 + 1 p M 2
where V is the orthogonal matrix of the eigenvectors of M Ĉ2 and Λ ∇ is the diagonal matrix with

[Λ ∇ ] kk = 1 πp j =k λj ξj 1 m μp (x) 1 (λ k -x) 2 dx + 1 πp j =k λ k ξ k 1 m μp (x) 1 (λ j -x) 2 dx. Algorithm 1 Proposed estimation algorithm. Require Positive definite initialization M = M 0 . Repeat M ← M 1 2 exp -tM -1 2 ∇h(M )M -1 2 M 1 2
with t either fixed or optimized by backtracking line search. Until Convergence.

Return M . This finally entails the gradient descent Algorithm 1.

Figure 2 depicts the results of the algorithm. There is displayed the Wasserstein distance D W (C, •) between a matrix C having four distinct eigenvalues of equal multiplicity (precisely, ν p = 1 4 (δ .1 + δ 3 + δ 4 + δ 5 )) and various estimators of C: the sample covariance matrix (SCM), the state-of-theart "non-linear shrinkage" estimators QuEST1 [START_REF] Ledoit | Spectrum estimation: A unified framework for covariance matrix estimation and pca in large dimensions[END_REF] (based on a Frobenius distance minimization) and QuEST2 [START_REF] Ledoit | Optimal estimation of a largedimensional covariance matrix under stein's loss[END_REF] (based on a Stein loss minimization), and the result of the gradient descent approach proposed in this section. For fair comparison, the iterative QuEST1, QuEST2 and our proposed method are all initialized at M 0 the linear shrinkage estimator from [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]. Note that our proposed choice of C is particularly suited to mimick an "optimal transport" problem of displacing the eigenvalues of M 0 to the discrete four positions of the eigenvalues of C.

In addition to the computational simplicity of our gradientdescent approach with respect to the QuEST estimators (see the numerical method details in [START_REF] Ledoit | Numerical implementation of the quest function[END_REF]), the figure demonstrates significant gains brought by our proposed approach for large values of p/n, where the SCM particularly fails.

V. CONCLUDING REMARKS Interestingly, while the Fisher distance or Kullbach-Liebler divergence, which depend on logarithms of inverse of covariance matrices, are understandably difficult to estimate in the n 1 , n 2 < p regime (see [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] for advanced discussions on this matter), the Wasserstein distance should not be confronted with this limitation. Yet, the invertibility of C 1 , C 2 and the request for c 1 , c 2 ∈ (0, 1) (i.e., p < n 1 , n 2 ) from Assumption 1 are fundamental to our proofs. Precisely, the variable changes exploited in the proof of Theorem 1 to reach a contour Γ ν correctly surrounding supp(ν p ) from a contour Γ surrounding supp(µ p ) are not satisfying if c 1 > 1 or c 2 > 1. These surprising difficulties need clarification.

Another point of interest lies in the comparative advantage of exploiting a particular covariance matrix distance in specific scenarios. For instance, it may seem that ill-conditioned matrices should be more tolerated by Wasserstein distance estimators than by Fisher distance estimators. Yet, this aspect is not obvious in our proofs and also deserves more insights.

  λi(C1C2)) of the eigenvalues of the product of two covariance matrices C1, C2 ∈ R p×p based on the empirical estimates λi( Ĉ1 Ĉ2) ( Ĉa = 1

Fig. 1 .

 1 Fig. 1. Deformation of the initial contour Γ (in black) into the new contour Γ (in blue). The branch cuts are represented in green (i.e., real z's for which the argument of ϕ(z)ψ(z) is negative).

  0, Ca); n 1 = 1024 AND n 2 = 2048 FOR DIFFERENT p. AVERAGED OVER 100 TRIALS.

TABLE I ESTIMATORS

 I OF THE WASSERSTEIN DISTANCE BETWEEN C 1 AND C 2

perturbation of Λ, by Weyl's interlacing lemma[START_REF] Joel N Franklin | Matrix theory[END_REF], the ξ j 's are interlaced with the λ j 's asξ 1 ≤ λ 1 ≤ ξ

≤ . . . ≤ ξ p ≤ λ p .As the λ j 's are of order O(1) with respect to p, |λ j -ξ j | ≤ |λ j -λ j-1 | = O(p -1), therefore explaining why the expression of D(X 1 , X 2 ; √ •) is of order O(1).

ACKNOWLEDGEMENT

We thank Pedro Rodrigues for helpful discussions and references on optimal transport.

) and (green) our proposed estimator, (blue) the sample covariance matrix, (red) and (light blue) the QuEST estimators proposed in [START_REF] Ledoit | Optimal estimation of a largedimensional covariance matrix under stein's loss[END_REF], [START_REF] Ledoit | Spectrum estimation: A unified framework for covariance matrix estimation and pca in large dimensions[END_REF]; for p = 100 and varying number of samples n averaged over 10 realizations.


domain: info math stat

text:
Malik Tiomoko

Romain Couillet

ESTIMATION OF COVARIANCE MATRIX DISTANCES IN THE HIGH DIMENSION LOW SAMPLE SIZE REGIME

Keywords: Random Matrix Theory, Statistical Inference, Covariance Matrix

A broad family of distances between two covariance matrices C1, C2 ∈ R p×p , among which the Frobenhius, Fisher, Battacharrya distances as well as the Kullback-Leibler, Rényi and Wasserstein divergence for centered Gaussian distributions can be expressed as functionals

Consistent estimates of such distances based on few (n1, n2) samples xi ∈ R p having covariance C1, C2 have been recently proposed using random matrix tools in the regime where n1, n2 ∼ p. These estimates however demand that n1, n2 > p for most functions f . The article proposes to alleviate this limitation using a polynomial approximation approach. The proposed method is supported by simulations in practical applications.

INTRODUCTION

Evaluating the distance between covariance matrices is at the core of many machine learning and signal processing applications. They are notably used for covariance features-based classification (for instance in brain signal or hyperspectral image classification), as well as for dimensionality reduction and representation of high dimensional points. Denote C1, C2 ∈ R p×p two large dimensional covariance matrices for which we would like to compute the distance D(C1, C2) based on few (p-dimensional) sample vectors. We assume that D can be written as a linear functional 1 p p i=1 f (li) of the eigenvalue distribution of either C - 1 1 C2 (li = λi(C -1 1 C2)) (as with the Fisher, Battacharrya, Kullback Leibler, Rényi divergences) or C1C2 (li = λi(C1C2)) (for the Wasserstein distance).

The natural estimation of D(C1, C2) based on na samples x a i , a = 1, 2, is traditionally performed via D( Ĉ1, Ĉ2), where Ĉa =

1 na na i=1 x (a) i x (a)T i
. This is justified by the law of large number since D( Ĉ1, Ĉ2) a.s.

-→ D(C1, C2) as n1, n2 p. However, the estimate induces dramatic erroneous when n1, n2 ∼ p and even diverges for na < p.

To deal with the n1, n2 ∼ p scenario, [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] and [START_REF] Tiomoko | Random matrix improved covariance estimation for a large class of metrics[END_REF] proposed a random matrix improved consistent estimate for D(C1, C2) for all aforementioned metrics. The main idea behind these works consists in three steps: i) write the sought-for distance D(C1, C2) as a complex integral involving the Stieltjes transform of the limiting eigenvalue distribution of C -1 1 C2 (or C1C2), ii) exploit a functional identity relating the Stieltjes transforms of the limiting population and sample eigenvalue distributions (e.g., of C -1 1 C2 and Ĉ-1 1 Ĉ2) using the results from [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF] and iii) find (if possible) a closed form solution for the resulting complex integral.

However, although consistent for n1, n2 ∼ p, these improved estimators still demand that n1, n2 > p for all functions f (z) having a singularity at z = 0 (e.g., 1/z, log(z), log 2 (z), √ z). For functions f having no singularity at zero, this constraint may be relaxed (with some extra care though) as discussed in [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF], but almost no usual covariance matrix distance falls into this scenario (except for the Frobenius distance, which can at any rate be treated easily without resorting to these elaborate methods).

Based on a polynomial approximation of the functions of interest, this article proposes to retrieve consistent estimates for the challenging n1, n2 < p scenario. Besides, a closed-form and numerically convenient expression of the proposed estimator is derived for any (local) polynomial approximation of the functions.

The article is organized as follows. In section 2, we will explain briefly the general idea of the estimation and the problem induced by the challenging case n1 < p and n2 < p. A solution based on polynomial approximation is then proposed in section 3 to cover the case n2 < p. As an application, a dimensionality reduction framework shows the consistency of the proposed estimator in section 4.

Reproducibility. A set of Matlab codes for the various estimators introduced and studied in this article are available at https://github.com/maliktiomoko/RMTCovDistHDLSS.

PRELIMINARIES

Models and Assumption

For a ∈ {1, 2}, let Xa = [x (a) 1 , . . . , x (a)
na ] be na independent and identically distributed random vectors with x

(a) i = C 1 2 a x(a)
i , where x(a) i ∈ R p has zero mean, unit variance and finite fourth order moment entries. This holds in particular for x (a) i ∼ N (0, Ca). In order to control the growth rates of n1, n2, p, we make the following assumption:

Assumption 1 (Growth Rates). As na → ∞, p/na → ca ∈ (0, ∞) and lim sup p max{ C -1 a , Ca } < ∞ for • the operator norm.

We define the sample covariance estimate Ĉa of Ca as

Ĉa ≡ 1 na XaX T a = 1 na na i=1 x (a) i x (a)T i .
Our objective is to estimate, under the difficult regime where p > na, the distance D(C1, C2) between the covariance matrices C1 and C2 of the form:

D(C1, C2) = 1 p p i=1 f (li)
where li = l - i are the eigenvalues of C -1 1 C2 or li = l + i are the eigenvalues of C1C2. As discussed at length in [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF], this form comprises, among others, the Fisher, Frobenius, Wasserstein and Battacharrya distances, along with the Kullbach-Liebler and Rényi divergences.

The proposed estimate relies on random matrix theory and particularly on the Stieltjes transform m θ (z) of a probability distribution θ defined as the complex function:

m θ : C \ supp(θ) → C, z → (λ -z) -1 dθ(λ).
The Stieltjes transform is here used to create a link between the population covariance eigenvalue distribution νp to the sample eigenvalue distribution µp defined by:

νp = 1 p p i=1 δ l i , µp = 1 p p i=1 δ li
where li = li are the eigenvalues of Ĉ-1 1 Ĉ2 or li = l+ i are the eigenvalues of Ĉ1 Ĉ2.

Previous Results

Theorem 1 from [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] and [START_REF] Tiomoko | Random matriximproved estimation of the wasserstein distance between two centered gaussian distributions[END_REF] provide an estimate of these distances in the "easy" regime where lim p/na < 1.

Theorem 1. Let f : C → C be analytic on a contour Γ ⊂ {z ∈ C, R[z] > 0} surrounding µp. Then, f dνp - 1 2πi Γ f ϕp(z) ψp(z) ψ p (z) ψp(z) - ϕ p (z) ϕp(z) ψp(z)dz c2 a.s.
-→ 0

where

ϕp(z) = z(1 + c1zmµ p (z)), µp = µ - p z 1-c 1 -c 1 zmµ p (z) , µp = µ + p ψp(z) = 1 -c2 -c2zmµ p (z).
To understand the generalization of Theorem 1 proposed in the present article, we need to recall the main steps of its proof (here for µp = µ - p ).

Sketch of Proof.

Using the Cauchy integral formula, we have

D(C1, C2) = 1 p p i=1 f (li) = f (t)νp(dt) (1)
= 1 2πı Γν f (z) z -t νp(dt) = -1 2πı Γν f (z)mν p (z)dz.
Thus, estimating D(C1, C2) is equivalent to relating mν p to mµ p . Since X1 and X2 are independent, we can condition first on X2. By [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF], the limiting eigenvalue distribution of C2 Ĉ-1 1 , denoted ζ, can be written as a function C2C -1 1 , and similarly for the limiting eigenvalue distributions of Ĉ2 Ĉ-1 1 and C2 Ĉ-1 1 . This entails the two equations:

zmµ p (z) = ϕp(z)m ζp (ϕp(z)) + op(1) (2) mν p (z/Ψp(z)) = m ζp (z)Ψp(z) + op(1). (3)
Through the changes of variable z → ϕp(z) and ω → Ψp(ω) applied in Γν f (z)mν p (z)dz, the result follows. (Bottom) n2 < p. For n2 < p, the left real crossing of Γν is necessarily negative (even if the mass at {0} of Sµ were not included in Γ). In case of singularities or branch-cuts (shown in red for the log(z) and √ z functions), the contours are invalid.

On the need for n2 > p and n1 > p in Theorem 1

Since distances involving the eigenvalues of C -1 1 C2 are estimated from the empirical matrix Ĉ-1 1 Ĉ2, the constraint n1 > p is inevitable to ensure the existence of Ĉ-1

1 . The requirement for n2 > p is less immediate. The two variable changes discussed in the proof of Theorem 1 are only licit if they realize a mapping from a contour Γ enclosing the limiting support Sµ of µp and a valid contour Γν enclosing the limiting support Sν of νp while enclosing no additional singular points of the function f (z) (otherwise the Cauchy formula used in (1) is incorrect). But for n2 < p, it can be proved (see details in [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF]) that the pre-image of Γµ by the variable changes wraps around Sν and around zero (with leftmost real crossing depending on the ratio n2/p). This is a problem for all functions f (z) singular at z = 0. In particular, 1/z, log(z), log 2 (z) and √ z are examples of such invalid functions which, for some, additionally have a branch-cut terminating at zero (that no valid contour may cross). This discussion is most conveniently illustrated in Figure 1.

Unfortunately, there seems to be no simple workaround to this situation. We propose in this article to (partially) solve the problem by introducing entire functions (thus analytical over C) as substitutes for the locally non-analytic functions log(z), log 2 (z) and √ z intervening in the distance D(C1, C2) formulation.

MAIN RESULTS

Proposed estimation

Our approach consists in approximating (arbitrarily closely) the analytic functions f under study that present singularities around zero by entire functions, and particularly degree-N polynomials fN (z) defined by fN (z) = N n=0 anz n . Our central argument relies in the fact that, since Ca and C -1 a are bounded (as per Assumption 1), the limiting support Sν of νp is a compact set strictly away from zero. As such, one needs not approximate f on the whole R + half-line (which would still pose problems in the vicinity of zero) but only on a subset [a, b] ⊂ (0, ∞) over which polynomials are universal approximators.

This gives rise to the following immediate extension of Theorem 1.

Theorem 2. Let Γ be a contour surrounding the limiting support Sµ of µp. Then, for all > 0, there exists N and fN (z) = N n=0 anz n a polynomial of order N , such that, under the notations of Theorem 1, lim sup

D(X1, X2; fN ) = n2 2πip Γ fN ϕp(z) ψp(z) ϕ p (z) ϕp(z) - ψ p (z) ψp(z) ψp(z)dz -a0 1-c2 c2 .
The result immediately follows from our discussion above and Theorem 1. To avoid the numerical integration and proper choice of a contour in the statement of Theorem 2, we next provide closedform expressions for the estimate D(X1, X2; fN ). For simplicity of exposition, we consider here only the case n2 < n1.

Theorem 3 (Case C1C2). Let 0 < λ1 ≤ . . . ≤ λn 2 be the n2 non zero eigenvalues of Ĉ1 Ĉ2, and define {ξi} n 2 -1 i=1 and {ηi} n 2 i=1 the (increasing non zero) eigenvalues of Λ -1

n 1 √ λ √ λ T and Λ -1 n 2 √ λ √ λ T
, respectively, where λ = (λ1, . . . , λn 2 ) T and Λ = diag(λ). Then,

D(X1, X2; fN ) = a0 + an c2 n 2 -1 j=1 ξj 1 + c1 c2 -c1 + N n=2 an c2 n 2 -1 j=1 1 (n -2)! ∂ n-2 ∂z n-2 A ξ j n (z)C ξ n (z) | z=ζ j + N n=1 an c2 n 2 j=1 1 (n -1)! ∂ n-1 ∂z n-1 A η j n (z)C η n (z) |z=η j where A ξ j n (z) = ϕ(z) ψ(z) (z -ξj) n ψ(z) z -ξj A η j n (z) = ϕ(z) ψ(z) (z -ηj) n ψ(z) C ξ n (z) = -1 n -1 ϕ (z) ϕ(z) , C η n (z) = -1 n ψ (z) ψ(z) .
Theorem 4 (Case C -1 1 C2). Let 0 < λ1 ≤ . . . ≤ λn 2 the n2 non zero eigenvalues of Ĉ-1 1 Ĉ2, and define {ξi} n 2 -1 i=1 and {ηi} n 2 i=1 the (increasing non zero) eigenvalues of Λ -

1 p-n 1 √ λ √ λ T and Λ -1 n 2 √ λ √ λ T
, respectively, where λ = (λ1, . . . , λn 2 ) T and Λ = diag(λ). Then,

D(X1, X2; fN ) = a0 + an c2 n 2 -1 j=1 ξj 1 + c1 c2 -c1 + N n=2 an c2 n 2 -1 j=1 1 (n -2)! ∂ n-2 ∂z n-2 A ξ j n (z)Cn(z) | z=ξ j + n 2 j=1 an p (1 -n1 + n2 -p) - c1 c2 λj n where A ξ j n (z) = ϕ(z) ψ(z) (z -ξj) n ψ(z) z-ξ j Cn(z) = -1 n-1 ϕ (z) ϕ(z) -1 z
Proof. The proofs of Theorems 3 and 4 are similar. For conciseness we only sketch the proof of the Theorem 4. The ξi and ηi are the respective zeros of the rational functions 1 -p n 2 -p n 2 zm(z) and 1 + p n 1 zm(z). Thus, ϕp and ψp can be expressed as the following rational functions:

ϕp(z) = (1 -c1) z n 2 i=1 z -ηi n 2 i=1 z -λi , ψp(z) = z n 2 -1 i=1 z -ξi n 2 i=1 z -λi
and their derivatives as

ϕ p (x) ϕp(x) - ψ p (x) ψp(x) = n 2 i=1 1 x -ηi - n 2 -1 i=1 1 x -ξi .
The complex integral can thus be written as N n=1 an 2πıc 2

In(z)dz with the integrand In(z) given by the rational function:

In(z) = (1 -c1) n n 2 i=1 (z -ηi) n n 2 -1 i=1 (z -ξi) n-1 n 2 i=1 (z -λi) n × n 2 i=1 1 z -ηi - n 2 -1 i=1 1 z -ξi
for which we need to evaluate the residue at poles:

• the (n -1)-th and n-th order pole at ξj with residue (for n ≥ 2) by:

R1 = n 2 -1 j=1 lim z→ξ j 1 (n -2)! ∂ n-2 ∂x n-2 In(z)(z -ξj) n-1 + n 2 -1 j=1 lim z→ξ j 1 (n -1)! ∂ n-1 ∂x n-1 [In(z)(z -ξj) n ]
• the 1st order pole in λj with residue:

R3 = n 2 j=1 lim z→λ j In(z)(z -λj)
• the 1st order pole in ηj with residue (for n = 0):

R4 = n 2 j=1 lim z→η j ψ(z).
Putting all residues together and exploiting further relations involving the functions ϕp and ψp and their derivatives (see [START_REF] Couillet | Random matrix-improved estimation of covariance matrix distances[END_REF] for details) entails the result.

Although Theorems 3 and 4 are difficult to interpret, they are numerically easy to implement: the terms involving the k-th derivative can indeed be evaluated iteratively by remarking the following identities for any rational functional function

A(z) = i z-a i i z-b i and C(z) = i 1 z-c i -1 z-d i (ai, bi, ci, di ∈ R): ∂ n ∂z n (A(z)C(z)) = n k=1 k n ∂ k ∂z k A(z) ∂ n-k ∂z n-k C(z) ∂ n ∂z n A(z) = ∂ n-1 ∂z n-1 A(z)B(z) = n-1 k=1 k n -1 ∂ k ∂x k A(z) ∂ n-1-k ∂x n-1-k B(z) where B(z) = i 1 z-a i -1
z-b i , from which we obtain a recursive expression for the derivatives of A(z).

SIMULATIONS AND APPLICATIONS

Applications to synthetic Gaussian data

We confirm the consistency of the proposed estimates through simulations on synthetic Gaussian data, here for the Wasserstein distance (that is, for f (z) = √ z and µp = µ + p ) and for the p = 2n2 = 1.5n1. Table 1 compares the proposed estimator to the traditional sample covariance matrix plugin estimate DSCM = D( Ĉ1, Ĉ2, f ). It is clear in this setting that, while DSCM to be largely erroneous, the proposed estimate performs well down to rather small values of p. Somewhat surprisingly, polynomials of high orders, which threaten the stability of the method (due to the not-so-low probability of occurrence of large eigenvalues for Ĉ1 Ĉ2), display a stable behavior even for small values of p. This suggests the possibility, in practice, to run the proposed estimator for rather large values of N .

∼ N (0, Ca); c1 = 1.5 and c2 = 2. Averaged over 10 trials. The function fN is the polynomial that best fits (in the least-squares sense)

√

D for D sampled logarithmically on the support of the eigenvalues of Ĉ1 Ĉ2. In boldface, distance estimates within 1% of relative error.

Application to dimensionality reduction

More interestingly in a statistical learning context, we now demonstrate the performance of the proposed estimator in a dimensionality reduction scenario. Consider m data points X1, . . . , Xm, each Xi ∈ R p×n i being obtained from ni independent p-dimensional centered Gaussian samples. For The objective is to depict the data in a two-dimensional projection space in order to identify the two classes based on their covariance matrices. To this end, we use kernel PCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF] for a Wasserstein distance metric, which consists in a (2D) principal component analysis of the kernel matrix K ∈ R m×m with Kij = exp(-1 2 DW (C (i) , C (j) )). Being unknown, the distances are approximated using either the classical estimate DW ( Ĉ(i) , Ĉ(j) ) or the proposed estimator of the Wasserstein distance.

i ≤ m/2, E[[Xi]•j[Xi] •j ] = C (i) = C1 =
Figure 2 shows that the proposed estimator preserves the local neighborhood structure of the data where the classical estimate dramatically fails. Fig. 2. First versus second eigenvectors of K for distinct uniformly random ni ∈ {p/2, . . . , p}, p = 512 and polynomial fN chosen similarly as in Table 1 with N = 6.

CONCLUDING REMARKS

The polynomial approximation approach proposed in this article breaks the stringent limitation that demands p < n for a host of estimators, here for covariance distance estimation. Yet, it presents some technical difficulties when dealing with functions with high variations near zero (such as log 2 (x)). Combining polynomials P (x) and other family of better behaving functions (e.g., exp(-P (x))) may leverage the problem.

Moreover, there is still a need to cover the case p > n in estimators involving covariance matrix inverses for which Ĉ1 -1 Ĉ2 is not even defined. Random projections and regularization methods can be devised to tackle this scenario, however possibly to the detriment of the estimator consistency.

Fig. 1 .

 1 Fig. 1. Illustration of the contours maps Γ → Γν (from right to left) by the variable changes leading up to Theorem 1. (Top) n2 > p.(Bottom) n2 < p. For n2 < p, the left real crossing of Γν is necessarily negative (even if the mass at {0} of Sµ were not included in Γ). In case of singularities or branch-cuts (shown in red for the log(z) and √ z functions), the contours are invalid.

  , C2; f ) -D(X1, X2; fN )| < almost surely, where

  Ip and for i > m/2, C (i) = C2 a random (pdimensional) standard Wishart matrix with 100p degrees of freedom (to avoid trivial classification).

Class 1

 1 (proposed) Class 2 (proposed) Class 1 (classical) Class 2 (classical)

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


Marion Poupard
email: marion.poupard@univ-tln.fr

Paul Best

Jan Schlüter

Helena Symonds

Paul Spong

Thierry Lengagne

Thierry Soriano

Hervé É Glotin
email: glotin@univ-tln.fr

Jan Schl Üter

Large-scale unsupervised clustering of Orca vocalizations: a model for describing Orca communication systems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The Orca (Orcinus orca) is a top-predator of the marine food chain [START_REF] Jefferson | A review of killer whale interactions with other marine mammals: Predation to co-existence[END_REF]. The Northern Resident Orcas community is composed of several "pods" composed of matrilines [START_REF] Bigg | Social organization and genealogy of resident killer whales (orcinus orca) in the coastal waters of british columbia and washington state[END_REF]. This cetacean can produce 3 different types of signals: clicks, whistles and pulsed calls [START_REF] Ford | Acoustic behaviour of resident killer whales (orcinus orca) off vancouver island, british columbia[END_REF]. This study focuses only on vocalizations (pulsed calls). Some biological studies describe the communication of Orcas [START_REF] Ford | A catalogue of underwater calls produced by killer whales (Orcinus orca) in British Columbia[END_REF][START_REF] Tyson | Nonlinear phenomena in the vocalizations of north atlantic right whales and killer whales[END_REF][START_REF] Weiß | Intra-and intergroup vocal behavior in resident killer whales, orcinus orca[END_REF][START_REF] Filatova | Call diversity in the north pacific killer whale populations: implications for dialect evolution and population history[END_REF], based on manual methods. Related work by [START_REF] Deecke | Quantifying complex patterns of bioacoustic variation: Use of a neural network to compare killer whale dialects[END_REF] compared dialects of Orcas using artificial neural networks and showed that acoustic similarities are significantly correlated with the group association patterns. In order to analyze the animal's communication in different spacial and temporal contexts, automated analysis for captured sound is crucial. For that purpose, the field of bioacoustics offers numerous approaches using neural networks and deep learning [START_REF] Glotin | Proc. nips4b : Neural information processing scaled for bioacoustics, from neurons to big data[END_REF]. The latter methods were explored to automatically detect orca calls emitted throughout 3 years of continuous recordings from 2015 to 2017 (Poupard et al., 2019a). In this study we build on these detections, and compute each vocalisation's pitch over time. This pitch analysis serves to differentiate vocalisations. In particular, we extract pitch features and cluster the vocalizations, partly recovering different call types as annotated by human experts.

MATERIAL

For 20 years, the NGO Orcalab developed and has maintained a unique multi-hydrophone recording station around Hanson Island (Northern Vancouver Island, Canada) to study Orcas. This acoustic station is composed of 5 hydrophones and extends over 50 km 2 of ocean (Fig. 1). In 2015, we have set up a continuous recording of all the hydrophones of this station (Fig. 1). The aim is to allow observation and modelling of bioacoustic activities for various species, at large spacial and temporal scales, including all details of their ecoacoustic niche, under various geophysical and anthropophonic conditions, more particularly in order to build new knowledge about Orcas.

DATA ANALYSIS

Vocalization Detector

We first designed an automatic acoustic event extractor (presented in [START_REF] Poupard | Passive acoustics to monitor flagship species near boat traffic in the UNESCO world heritage natural reserve of scandola[END_REF]). Its output helped us build a dataset composed of 872 Orca vocalization samples and 6801 noise samples (boats, rain, void. . . ), which we split randomly with 20% for the test set, 60% for the training set and 20% for the validation set.1 With that in hand, we trained a CNN (originally designed for a bird detection task [START_REF] Grill | Two convolutional neural networks for bird detection in audio signals[END_REF] to distinguish Orca vocalizations (not clicks) from boats and background noise (Poupard et al., 2019a). The model was trained with weakly annoted data (one label per file), originally using global max pooling to aggregate local predictions for comparison against the global label. After training, the global pooling was removed to obtain local probabilities for pitch and vocalization analysis. Max pooling lead to spiky local predictions (high precision, but low recall), which were unsuitable for our purposes. We found that training the model with global mean pooling instead gave much higher recall, covering the full length of each vocalization without sacrificing precision. The resulting Area Under the receiver operating characteristic Curve (AUC) of this detector is 89% (Poupard et al., 2019a).

Running this model on the summers of 2015, 2016 and 2017 results in 421,879 detected vocalizations across all five hydrophones.

Pitch Analysis

In order to describe call characteristics, the pitch (fundamental frequency) is often used [START_REF] Berthommier | A measure of speech and pitch reliability from voicing[END_REF]. The pitch is a property that describes the fundamental frequency of the speech wave [START_REF] Houtsma | Pitch and timbre: Definition, meaning and use[END_REF][START_REF] Babacan | A comparative study of pitch extraction algorithms on a large variety of singing sounds[END_REF]. Like humans, Orcas produce vocalizations that have several harmonic frequencies, combining into a multi-layered wave [START_REF] Ford | Acoustic behaviour of resident killer whales (orcinus orca) off vancouver island, british columbia[END_REF]. [START_REF] Foote | Variation in call pitch among killer whale ecotypes[END_REF] used pitch to differentiate different ecotypes of killer whales and [START_REF] Shapiro | A versatile pitch tracking algorithm: From human speech to killer whale vocalizations[END_REF] developed their own pitch tracker algorithm (PDA) based on human voice.

In this study, the Parselmouth Python library [START_REF] Jadoul | Introducing parselmouth: A python interface to praat[END_REF] was chosen as pitch estimator. It relies on the autocorrelation (AC) [START_REF] Boersma | Accurate short-term analysis of the fundamental frequency and the harmonics-to-noise ratio of a sampled sound[END_REF][START_REF] Berthommier | A measure of speech and pitch reliability from voicing[END_REF]. It is illustrated on a recording of Orca calls in Fig. 2. Computing all the pitches for one day of vocalizations on the 5 hydrophones took half an hour in average on GPU. The AC only outputs a pitch point if it has a certain confidence in it (using a threshold on the strength of the unvoiced candidate relative to the maximum possible AC). Thus, with some detected vocalizations, fewer points were output. This property allowed us to filter out false positives and too low Signal to Noise Ratio vocalization detections. Keeping only vocalizations with more than 200 points filtered out 284,791 noisy vocalizations and false detections.

We thus extracted the pitch of 137,088 vocalizations.

Unsupervised Clustering

Unsupervised clustering is often the solution to solve classification tasks for unannotated data. Our intuition was that the Orcas' call types [START_REF] Ford | A catalogue of underwater calls produced by killer whales (Orcinus orca) in British Columbia[END_REF][START_REF] Root-Gutteridge | Improving individual identification in captive eastern grey wolves (canis lupus lycaon) using the time course of howl amplitudes[END_REF]) could be automatically clustered by similarity in their pitch shape. A first step was thus to define the input to our unsupervised clustering algorithm. Thus features of the previously extracted pitch were selected to best describe the shape of the vocalizations with a minimum dimensionality.

The following features were chosen [START_REF] Ford | Call traditions and dialects of killer whales (Orcinus orca) in British Columbia[END_REF]: argminFreq, argmaxFreq, minVel, maxVel, meanVel, startVel, endVel, minAccel, maxAccel, argminAccel, argmaxAccel, deltaFreq. Here argmin/max stand for the position in time of the maximum/minimum relative to the total duration. Min/max stand for minimum/maximum values. Mean stands for the average value. Start/end stand for the mean of the first/last 5% of the call. Delta stands for the minimum value substracted from the maximum value. Freq stands for frequency values (the estimated pitch), Vel stands for velocity (the derivative of the pitch), and Accel stands for acceleration (the derivative of the velocity).

Having extracted those features, they were used as an input for the HDBSCAN algorithm [START_REF] Mcinnes | hdbscan: Hierarchical density based clustering[END_REF], which is a hierarchical implementation of the Density Based Spatial Clustering of Applications with Noise (DBSCAN) [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. Several minimum cluster sizes and minimum sample sizes were explored, to optimize the number of output clusters, and the strictness of the clustering. Eventually 30 and 3 were chosen for the latter parameters respectively.

RESULT

The clustering algorithm hardly works when applied to all the collected vocalisations together (coming from the 5 different hydrophones with different depth and sensitivity), whereas it works decently when applied to each hydrophone separately. Here we present the results for the H1 hydrophone (see map in Fig. 1), which represents 6796 vocalizations. Further work will focus on generalizing the clustering method to any hydrophone after some normalization.

The HDBSCAN found 13 clusters (0 to 12; Fig. 4). The '-1' cluster is the algorithm's output of classifying what it considers as noise. To measure the clustering's relevance, 2 trained persons annotated 50 samples (picked randomly) from each cluster, according to the Oca call types as defined in the literature [START_REF] Ford | A catalogue of underwater calls produced by killer whales (Orcinus orca) in British Columbia[END_REF]. We selected a subset of 6 call types (N1,N2,N4,N7,N9,N47): the ones most commonly found in our dataset (Fig. 3).

The distributions of call types among clusters (Fig. 4) show that our model was able in some clusters to isolate some type (N4 in clusters 0, 1, 4, and 5), to group calls with roughly similar upward types (N2, 3/6 Figure 3. Selected subset of call types as defined in the literature [START_REF] Ford | A catalogue of underwater calls produced by killer whales (Orcinus orca) in British Columbia[END_REF]. N7, N9, N47 in clusters 2 and 3), and to classify boat noise (clusters 9 and 12). Those results demonstrate a promising approach to classifying orca vocalizations, in approximately 20 days of computation for 3 years of pentaphonic continuous recording.

CONCLUSION AND DISCUSSION

Our primary results can be linked with the presence of particular pods in the area. In fact, British Columbia (BC) is composed of different "acoustic clans". An acoustic clan is a group of Orcas that share particular types of calls known as discrete calls [START_REF] Ford | Killer whale (orcinus orca) dialects as an indicator of stocks in british columbia[END_REF]. In the Northern and coastal BC, there are 4 main acoustic clans: the J, R, G and A Clan. For now, we will focus on the A Clan (Fig. 5), composed of several pods, themselves composed of groups of lineages called matrilines.

As shown in Fig. 5, there are 7 different pods in the A clan, having different call types [START_REF] Ford | Call traditions and dialects of killer whales (Orcinus orca) in British Columbia[END_REF]. For example the A4 pods can produce N2, N4, N7 and N9 call types. Some types of calls are shared among multiple pods within the clan. For example, the N7 call extends to all pods, however each pod produces an unique N7 call. By recognizing pods and recording the different calls, we can establish a link between the pods and our clusters. In fact clusters 0, 1, 4, 5, and -1 have a high proportion of the N4 call 4/6 type (Fig. 4), we can thus expect that the A1, A4, and A5 pod vocalizations are present in these 5 clusters. The N47 call type is produced only by the A1 pods and this type is very present in 2 clusters (3 and 6), so we can state the hypothesis that these 2 clusters correspond to the A1 pod.

With such reasoning, these clusters represent a first approach to acoustically classify pods in BC, and in the future, matrilines [START_REF] Weiß | Vocal behavior of resident killer whale matrilines with newborn calves: The role of family signatures[END_REF] and individual vocal signatures [START_REF] Weiß | Intra-and intergroup vocal behavior in resident killer whales, orcinus orca[END_REF]. [START_REF] Ford | Call traditions and dialects of killer whales (Orcinus orca) in British Columbia[END_REF]. Future work will improve the model at each of the 3 main steps: the learned vocalization detection, the pitch estimation that could be trained specifically to detect Orcas' pitch [START_REF] Kim | Crepe: A convolutional representation for pitch estimation[END_REF], and the unsupervised clustering of calls. An obvious improvement would consist in annotating more data for training. Parameter optimization is another possible enhancement, especially for the pitch estimation and the unsupervised clustering. For this purpose, relevant objective functions and accurate metrics need to be found. One could consider a global objective cost function to maximise the normalised mutual information of the bivariate distribution (Type, Cluster).

Once such an improved system is at hands, having a fully autonomous reliable Orca type call detector and classifier will open doors to many studies on Orca's communication and phonotactic regularities and divergence like in [START_REF] Malige | Use of recurrence plots to automatically extract songs in humpback whale recordings[END_REF]. It would also allow behavioural studies (ethoacoustics), within various environments, including increasing anthropophony or whale watching pressure like in [START_REF] Poupard | Ethoacoustic by bayesian non parametric and stochastic neighbor embedding to forecast anthropicpressure on dolphins[END_REF].

Figure 1 .

 1 Figure 1. Left: Map of the area and the listening range of the 5 hydrophones H1 to H5. Detection zones indicate which hydrophones can record Orca calls in a given area, w.r.t. thirty years of observations by Orcalab. Right: Representation of the data acquisition, from recording until storage at Toulon.

Figure 2 .

 2 Figure 2. Example of a pitch extraction (pitch floor=300, pitch ceiling=2500, voicing threshold=0.2).

Figure 4 .

 4 Figure 4. Distribution of call types among clusters found by the HDBSCAN algorithm. The numbers next to the cluster name show the amount of cluterized samples. The % are from the annotated subset.

Figure 5 .

 5 Figure 5. Selection of the 6 Call types produced by pods of the A clan inspired from[START_REF] Ford | Call traditions and dialects of killer whales (Orcinus orca) in British Columbia[END_REF].

A random split may sample train and test segments from nearby locations, giving an overly optimistic test error. We did not have enough annotated data for a chronological split avoiding this.

ACKNOWLEDGMENTS

We thank first the Orcalab direction and collaborators for their incredible inspired work. We thank Biosong SA for the PhD funding of M. Poupard. This research is partly funded by FUI 22 Abyssound, ANR-18-CE40-0014 SMILES, ANR-17-MRS5-0023 NanoSpike and MARITTIMO EUR GIAS projects on advanced studies on cetaceans. We thank MI CNRS MASTODONS SABIOD.org and EADM MADICS CNRS scaled bioacoustic research groups, IUF for Glotin's chair 2011-16 during which he installed the remote recording at Orcalab to UTLN DYNI, and SEAMED PACA project and CNRS platform support for J. Schlüter's Post doc grant.