Jean Jacod
email: jean.jacod@upmc.fr

Mathieu Rosenbaum
email: mathieu.rosenbaum@upmc.fr

Estimation of volatility functionals: the case of a √ n window

Keywords: semimartingale, high frequency data, volatility estimation, central limit theorem, efficient estimation MSC2010: 60F05, 60G44, 62F12

 this was done with k n → ∞ and k n √ ∆ n → 0, and a central limit theorem was given after suitable de-biasing. Here we do the same with the choice k n ≍ 1/ √ ∆ n . This results in a smaller bias, although more difficult to eliminate.

Introduction

Consider an Itô semimartingale X t , whose squared volatility c t (a d × d matrices-valued process if X is d-dimensional) is itself another Itô semimartingale. The process X is observed at discrete times i∆ n for i = 0, 1, • • • , the time lag ∆ n being small (high-frequency setting) and eventually going to 0. The aim is to estimate integrated functionals of the volatility, that is t 0 g(c s ) ds for arbitrary (smooth enough) functions g, on the basis of the observations at stage n and within the time interval [0, t].

In [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF], to which we refer for detailed motivations for this problem, we have exhibited estimators which are consistent, and asymptotically optimal, in the sense that they asymptotically achieve the best rate 1/ √ ∆ n , and also the minimal asymptotic variance in the cases where optimality is well-defined (namely, when X is continuous and has a Markov type structure, in the sense of [START_REF] Clément | An infinite dimensional convolution theorem with applications to the efficient estimation of the integrated volatility[END_REF]). These estimators have this rate and minimal asymptotic variance as soon as the jumps of X are summable, plus some mild technical conditions. The aim of this report is to complement [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF] with another estimator, of the same type, but using spot volatility estimators based on a different window size. In this introduction we explain the differences between the estimator in that paper and the one presented here.

For the sake of simplicity we consider the case when X is continuous and one-dimensional (the discontinuous and multi-dimensional case is considered later), that is of the form (X (i+j)∆n -X (i+j-1)∆n ) 2 (1. [START_REF] Alvarez | Estimation of the instantaneous volatility[END_REF] for an arbitrary sequence of integers such that k n → ∞ and k n ∆ n → 0: One knows that V n (g) t P -→ V (g) t (when g is continuous and of polynomial growth).

X t = X 0 +
The variables c n i are spot volatility estimators, and according to [START_REF] Jacod | Discretization of Processes[END_REF] we know that c n [t/∆n] estimate c t , with a rate depending on the "window size" k n . The optimal rate 1/∆ 1/4 n is achieved by taking k n ≍ 1/ √ ∆ n . When k n is smaller, the rate is √ k n and the estimation error is a purely "statistical error"; when k n is bigger, the rate is 1/ √ k n ∆ n and the estimation error is due to the variability of the volatility process c t itself (its volatility and its jumps). With the optimal choice k n ≍ 1/ √ ∆ n the estimation error is a mixture of the statistical error and the error due to the variability of c t .

In [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF] we have used a "small" window, that is k n ≪ 1/ √ ∆ n . Somewhat surprisingly, this allows for optimality in the estimation of t 0 g(c s ) ds (rate 1/ √ ∆ n and minimal asymptotic variance). However, the price to pay is the need of a de-biasing term to be subtracted from V n (g), without which the rate is smaller and no Central Limit Theorem is available.

Here, we consider the window size k n ≍ 1/ √ ∆ n . This leads to a convergence rate 1/ √ ∆ n for V n (g) itself, and the limit is again conditionally Gaussian with the "minimal" asymptotic variance, but with a bias that depends on the volatility of the volatility c t , and on its jumps. It is however possible to subtract from V n (g) a de-biasing term again, so that the limit becomes (conditionally) centered.

Section 2 is devoted to presenting assumptions and results, and all proofs are gathered in Section 3. The reader is referred to [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF] for motivation and various comments and a detailed discussion of optimality. However, in order to make this report readable, we basically give the full proofs, even though a number of partial results have already been proved in the above-mentioned paper, and with the exception of a few well designated lemmas.

The results

Setting and Assumptions

The underlying process X is d-dimensional, and observed at the times i∆ n for i = 0, 1, • • • , within a fixed interval of interest [0, t]. For any process we write ∆ n i Y = Y i∆n -Y (i-1)∆n for the increment over the ith observation interval. We assume that the sequence ∆ n goes to 0. The precise assumptions on X are as follows:

First, X is an Itô semimartingale on a filtered space (Ω, F, (F t ) t≥0 , P). It can be written in its Grigelionis form, as follows, using a d-dimensional Brownian motion W and a Poisson random measure µ on R + ×E, with E is an auxiliary Polish space and with the (non-random) intensity measure ν(dt, dz) = dt ⊗ λ(dz) for some σ-finite measure λ on E:

X t = X 0 + t 0 b s ds + t 0 σ s dW s + t 0 E δ(s, z) 1 { δ(s,z) ≤1} (µ -ν)(ds, dz) + t 0 E δ(s, z) 1 { δ(s,z) >1} µ(ds, dz). (2.1) This is a vector-type notation: the process b t is R d -valued optional, the process σ t is R d ⊗R d - valued optional, and δ = δ(ω, t, z) is a predictable R d -valued function on Ω × R + × E.
The spot volatility process c t = σ t σ * t ( * denotes transpose) takes its values in the set M + d of all nonnegative symmetric d × d matrices. We suppose that c t is again an Itô semimartingale, which can be written as

c t = c 0 + t 0 b s ds + t 0 σ s dW s + t 0 E δ(s, z) 1 { δ(s,z) ≤1} (µ -ν)(ds, dz) + t 0 E δ(s, z) 1 { δ(s,z) >1} µ(ds, dz). (2.2)
with the same W and µ as in (2.1). This is indeed not a restriction: if X and c are two Itô semimartingales, we have a representation as above for the pair (X, c) and, if the dimension of W exceeds the dimension of X one can always add fictitious component to X, arbitrarily set to 0, so that the dimensions of X and W agree.

In (2.2), b and σ are optional and δ is as δ; moreover b and δ are R d 2 -valued. Finally, we need the spot volatility of the volatility and "spot covariation" of the continuous martingale parts of X and c, which are

c ij,kl t = d m=1 σ ij,m t σ kl,m t , c ′i,jk t = d l=1 σ il t σ jk,l t .
The precise assumption on the coefficients are as follows, with r a real in [0, 1):

Assumption (A'-r): There are a sequence (J n ) of nonnegative bounded λ-integrable functions on E and a sequence (τ n ) of stopping times increasing to ∞, such that

t ≤ τ n (ω) =⇒ δ(ω, t, z) r ∧ 1 + δ(ω, t, z) 2 ∧ 1 ≤ J n (z) (2.3) Moreover the processes b ′ t = b t -δ(t, z) 1 { δ(t,z) ≤1} λ(dz) (which is well defined), c t and c ′
t are càdlàg or càglàd, and the maps t → δ(ω, t, z) are càglàd (recall that δ should be predictable), as well as the processes b t + δ(t, z) κ( δ(t, z) ) -1 { δ(t,z) ≤1} ) λ(dz) for one (hence for all) continuous function κ on R + with compact support and equal to 1 on a neighborhood of 0.

2

The bigger r is, the weakest Assumption (A-r) is, and when (A-0) holds the process X has finitely many jumps on each finite interval. The part of (A-r) concerning the jumps of X implies that s≤t ∆X s r < ∞ a.s. for all t < ∞, and it is in fact "almost" implied by this property. Since r < 1, this implies s≤t ∆X s < ∞ a.s. Remark 2.1 (A'-r) above is basically the same as Assumption (A-r) in [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF], albeit (slightly) stronger (hence its name): some degree of regularity in time seems to be needed for b, c, c ′ , δ in the present case.

A First Central Limit Theorem.

For defining the estimators of the spot volatility, we first choose a sequence k n of integers which satisfies, as n → ∞:

k n ∼ θ √ ∆ n , θ ∈ (0, ∞), (2.4)
and a sequence u n in (0, ∞]. The M + d -valued variables c n i are defined, componentwise, as

c n,lm i = 1 k n ∆ n kn-1 j=0 ∆ n i+j X l ∆ n i+j X m 1 { ∆ n i+j X ≤un} , (2.5)
and they implicitly depend on ∆ n , k n , u n .

One knows that c n n : the choice (2.4) is optimal, in the sense that it allows us to have the fastest possible rate by a balance between the involved "statistical error" which is of order 1/ √ k n , and the variation of c t over the interval [t, t + k n ∆ n ], which is of order √ k n ∆ n because c t is an Itô semimartingale (and even when it jumps), see [START_REF] Alvarez | Estimation of the instantaneous volatility[END_REF][START_REF] Jacod | Discretization of Processes[END_REF].

By Theorem 9.4.1 of [START_REF] Jacod | Discretization of Processes[END_REF], one also knows that under (A'-r) and if

u n ≍ ∆ ̟ n for some ̟ ∈ p-1 2p-r , 1 2 we have V (g) n t := ∆ n [t/∆n]-kn+1 i=1 g( c n i ) u.c.p. =⇒ V (g) t := t 0 g(c s ) ds (2.6)
(convergence in probability, uniform over each compact interval; by convention b i=a v i = 0 if b < a), as soon as the function g on M + d is continuous with |g(x)| ≤ K(1 + x p ) for some constants K, p. Actually, for this to hold we need much weaker assumptions on X, but we do not need this below. Note also that when X is continuous, the truncation in (2.5) is useless: one may use (2.5) with u n ≡ ∞, which reduces to (1.1) in the one-dimensional case. Now, we want to determine at which rate the convergence (2.6) takes place. This amounts to proving an associated Central Limit Theorem. For an appropriate choice of the truncation levels, such a CLT is available for V (g) n , with the rate 1/ √ ∆ n , but the limit exhibits a bias term. Below, g is a smooth function on M + d , and the two first partial derivatives are denoted as ∂ jk g and ∂ 2 jk,lm g, since any x ∈ M + d has d 2 components x jk . The family of all partial derivatives of order j is simply denoted as ∂ j g. Theorem 2.2 Assume (A'-r) for some r < 1. Let g be a C 3 function on M + d such that

∂ j g(x) ≤ K(1 + x p-j ), j = 0, 1, 2, 3 (2.7)
for some constants K > 0, p ≥ 3. Either suppose that X is continuous and u n /∆ ε n → ∞ for some ε < 1/2 (for example, u n ≡ ∞, so there is no truncation at all), or suppose that

u n ≍ ∆ ̟ n , 2p -1 2(2p -r) ≤ ̟ < 1 2 .
(2.8)

Then we have the finite-dimensional (in time) stable convergence in law

1 √ ∆ n (V (g) n t -V (g) t ) L f -s -→ A 1 t + A 2 t + A 3 t + A 4 t + Z t , (2.9)
where Z is a process defined on an extension ( Ω, F , ( F t ) t≥0 , P) of (Ω, F, (F t ) t≥0 , P), which conditionally on F is a continuous centered Gaussian martingale with variance

E (Z t ) 2 | F = d j,k,l,m=1 t 0 ∂ jk g(c s ) ∂ lm g(c s ) c jl s c km s + c jm s c kl s ds, (2.10)
and where, with the notation

G(x, y) = 1 0 g(x + wy) -wg(x + y) -(1 -w)g(x) dw, (2.11)
we have

A 1 t = -θ 2 g(c 0 ) + g(c t ) A 2 t = 1 2θ d j,k,l,m=1 t 0 ∂ 2 jk,lm g(c s ) c jl s c km s + c jm s c kl s ds A 3 t = -θ 12 d j,k,l,m=1 t 0 ∂ 2 jk,lm g(c s ) c jk,lm s ds A 4 t = θ s≤t G(c s-, ∆c s ).
(2.12)

Note that |G(x, y)| ≤ K(1 + x ) p y 2 , so the sum defining A 4 t is absolutely convergent, and vanishes when c t is continuous.

Remark 2.3 The bias has four parts:

1) The first one is due to a border effect: the formula giving V n (g) t contains [t/∆ n ]-k n + 1 summands only, whereas the natural (unfeasible) approximation ∆ n [t/∆n] i=1 g(c (i-1)∆n ) contains [t/∆ n ] summands. The sum of the lacking k n summands is of order of magnitude (k n -1)∆ n , which goes to 0 and thus does not impair consistency, but it creates an obvious bias after normalization by 1/ √ ∆ n . Removing this source of bias is straightforward: since g(c s ) is "under-represented" when s is close to 0 or to t, we add to V n (g) t the variable

(k n -1)∆ n 2 g( c n 1 ) + g( c n [t/∆n]-kn+1 . (2.13)
Of course, other weighted averages of g( c n i ) for i close to 0 or to [t/∆ n ]k n + 1 would be possible.

2) The second part A 2 is continuous in time and is present even for the toy model X t = √ c W t with c a constant and ∆ n = 1 n and T = 1. In this simple case it can be interpreted as follows: instead of taking the "optimal" g( c n ) for estimating g(c), with c n = n i=1 (∆ n i X) 2 , one takes 1 n n i=1 g( c n i ) with c n i a "local" estimator of c. This adds a statistical error which results in a bias. Note that, even in the general case, this bias would disappear, were we taking in (2.4) the (forbidden) value θ = ∞ (with still k n ∆ n → 0, at the expense of a slower rate of convergence.

3) The third and fourth parts A 3 and A 4 are respectively continuous and purely discontinuous, and due to the continuous part and to the jumps of the volatility process c t itself. These two biases disappear if we take θ = 0 in (2.4) (with still k n → ∞), again a forbidden value, and again at the expense of a slower rate of convergence. The only test function g for which the last three biases disappear is the identity g(x) = x. This is because, in this case, and up to the border terms, V (g) n t is nothing but the realized quadratic variation itself and the spot estimators c n i actually merge together and disappear as such.

Remark 2.4 Observe that (2.8) implies r < 1. This restriction is not a surprise, since one needs r ≤ 1 in order to estimate the integrated volatility by the (truncated) realized volatility, with a rate of convergence 1/ √ ∆ n . When r = 1 it is likely that the CLT still holds fore an appropriate choice of the sequence u n , and with another additional bias, see e.g. [START_REF] Vetter | Limit theorems for bipower variation of semimartingales[END_REF] for a slightly different context. Here we let this borderline case aside.

Estimation of the Bias.

Now we proceed to "remove" the bias, which means subtracting consistent estimators for the bias from V ′ (g) n t . As written before, we have

A n,1 t = - k n √ ∆ n 2 g( c n 1 ) + g( c n [t/∆n]-kn+1 P -→ A 1 t (2.14) (this comes from c n 1 P -→ c 0 and c n [t/∆n]-kn+1 P -→ c t-, plus c t-= c t a.s.). Next, observe that A 2 = 1 θ V (h) for the test function h defined on M + d by h(x) = 1 2 d j,k,l,m=1 ∂ 2 jk,lm g(x) x jl x km + x jm x kl . (2.15) Therefore A n,2 t = 1 k n √ ∆ n V (h) n t P -→ A 2 t . ( 2

.16)

The term A 3 t involves the volatility of the volatility, for which estimators have been provided in the one-dimensional case by M. Vetter in [START_REF] Vetter | Estimation of integrated volatility of volatility with applications to a goodness-of-fit testing[END_REF]; namely, if d = 1 and under suitable technical assumptions (slightly stronger than here), plus the continuity of X t and c t , he proves that

3 2k n [t/∆n]-2kn+1 i=1 ( c n i+kn -c n i ) 2 converges to t 0 c s + 6 θ 2 (c s ) 2 ds.
Of course, we need to modify this estimator here, in order to include the function ∂ 2 g in the limit and account for the possibilities of having d ≥ 2 and having jumps in X. We propose to take

A n,3 t = - √ ∆ n 8 [t/∆n]-2kn+1 i=1 d j,k,l,m=1 ∂ 2 jk,lm g( c n i ) ( c n,jk i+kn -c n,jk i ) ( c n,lm i+kn -c n,lm i
).

(2.17)

When X and c are continuous one may expect the convergence to

A 3 t -1 2 A 2 t (observe that √ ∆n 4 ∼ 3 2kn θ 12
), and one may expect the same when X jumps and c is still continuous, because in (2.5) the truncation basically eliminates the jumps of X. In contrast, when c jumps, the limit should rather be related to the "full" quadratic variation of c, and indeed we have: Theorem 2.5 Under the assumptions of Theorem 2.2, for all t ≥ 0 we have

A n,3 t P -→ - 1 2 A 2 t + A 3 t + A ′4 t , (2.18)
where

A ′4 t = θ s≤t G ′ (c s-, ∆c s ) (2.19)
and

G ′ (x, y) = - 1 8 j,k,l,m 1 0 ∂ 2 jk,lm g(x) + ∂ 2 jk,lm g(x + (1 -w)y) w 2 y jk y lm dw. (2.20)
At this stage, it remains to find consistent estimators for A 4 t -A ′4 t , which has the form

A 4 t -A ′4 t = θ s≤t G ′′ (c s-, ∆c s ), where G ′′ = G -G ′ .
More generally, we aim at estimating The solution to this problem is not as simple as it might appear at first glance. We first truncate from below, taking any sequence u ′ n of truncation levels satisfying

V(F ) t = s≤t F (c s-, ∆c s ), ( 2
u ′ n → 0, u ′ n ∆ ̟ ′ n → ∞ for some ̟ ′ ∈ 0, 1 8 (2.22)
Second, we resort on the following trick. Since c n i is "close" to the average of c t over the interval (i∆ n , (i + k n )∆ n ], we (somehow wrongly) pretend that, for all j:

∃s ∈ ((j -1)k n ∆ n , jk n ∆ n ] with ∆c s > u ′ n ⇔ c n jkn -c n (j-2)kn > u ′ n ∆c s ∼ c n jkn -c n (j-2)kn , c n (j-1)kn -c n (j-3)kn c n (j+1)kn -c n (j-1)kn < c n jkn -c n (j-2)kn . (2.23)
The condition (2.22) implies that for n large enough there is at most one jump of size bigger than u

′ n in each interval (i -1)∆ n , (i -1 + k n )∆ n ] within [0, t],
and no two consecutive intervals of this form contain such jumps. Despite this, the statement above is of course not true, the main reason being that c n i and c n i do not exactly agree. However it is "true enough" to allow for the next estimators to be consistent for V(F ) t :

V(F ) n t = [t/kn∆n]-3 j=3 F ( c n (j-3)kn+1 , δ n j c) 1 { δ n j-1 c ∨ δ n j+1 c ∨u ′ n < δ n j c } , where δ n j c = c n jkn+1 -c n (j-2)kn+1 . (2.24) Since this is a sum of approximately [t/k n ∆ n ] terms, the rate of convergence of V(F ) n t toward V(F ) t is law, probably 1/∆ 1/4 n
only. However, here we are looking for consistent estimators, and the rate is not of concern to us. Note that, again, the upper limit in the sum above is chosen in such a way that V(F ) n t is computable on the basis of the observations within the interval [0, t].

Theorem 2.6 Assume all hypotheses of Theorem 2.2, and let F be a continuous function on R + × R satisfying, with the same p ≥ 3 as in (2.8),

|F (x, y)| ≤ K(1 + x + y ) p-2 y 2 .
(2.25)

Then for all t ≥ 0 we have

V(F ) n t P -→ V(F ) t .
(2.26)

An Unbiased Central Limit Theorem.

At this stage, we can set, with the notation (2.15), (2.16), (2.17) and (2.24), and also (2.12) and (2.20) for G and G ′ :

V (g) n t = V (g) n t + k n ∆ n 2 g( c n 1 ) + g( c n [t/∆n]-kn+1 -∆ n 3 2 A n,2 t + A n,3 t -k n ∆ n V(G -G ′ ) n t . (2.27)
We then have the following, which is a straightforward consequence of the three previous theorems and of k n √ ∆ n → θ, plus (2.14) and (2.16) and the fact that the function G -G ′ satisfies (2.25) when g satisfies (2.7):

Theorem 2.7 Under the assumptions of Theorem 2.2, and with Z as in this theorem, for all t ≥ 0 we have the finite-dimensional stable convergence in law

1 √ ∆ n (V (g) n t -V (g) t ) L f -s -→ Z t .
(2.28)

Note that θ no longer explicitly appears in this statement, so one can replace (2.4) by the weaker statement

k n ≍ 1 √ ∆ n (2.29)
(this is easily seen by taking subsequences n l such that k n l ∆ n l converge to an arbitrary limit in (0, ∞)).

It is simple to make this CLT "feasible", that is usable in practice for determining a confidence interval for V (g) t at any time t > 0. Indeed, we can define the following function on M + d :

h(x) = d j,k,l,m=1
∂ jk g(x) ∂ lm g(x) x jl x km + x jm x kl .

(2.30)

We then have V (h) n u.c.p. =⇒ V (h), whereas V (h) t is the right hand side of (2.9). Then we readily deduce: Corollary 2.8 Under the assumptions of the previous theorem, for any t > 0 we have the following stable convergence in law, where Y is an N (0, 1) variable:

V (g) n t -V (g) t ∆ n V (h) n t L-s -→ Y, in restriction to the set {V (h) t > 0}, (2.31)
Finally, let us mention that the estimators V (g) n t enjoy excatly the same asymptotic efficiency properties as the estimators in [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF], and we refer to this paper for a discussion of this topic.

Example 2.9 (Quarticity) Suppose d = 1 and take g(x) = x 2 , so we want tho estimate the quarticity t 0 c 2 s ds. In this case we have

h(x) = 2x 2 , G(x, y) -G ′ (x, y) = - y 2 6 .
Then the "optimal" estimator for the quarticity is

∆ n 1- 3 k n [t/∆n]-kn+1 i=1 ( c n i ) 2 + ∆ n 4 [t/∆n]-2kn+1 i=1 ( c n i+kn -c n i ) 2 + (k n -1)∆ n 2 ( c n 1 ) 2 +( c n [t/∆n]-kn+1 ) 2 .
The asymptotic variance is 8 t 0 c 4 s ds, to be compared with the asymptotic variance of the more usual estimators 1 3∆n

[t/∆n] i=1 (∆ n i X) 4 , which is 32 3 t 0 c 4 s ds.
3 Proofs

Preliminaries.

According to the localization lemma 4.4.9 of [START_REF] Jacod | Discretization of Processes[END_REF] (for the assumption (K) in that lemma), it is enough to show all four Theorems 2.2, 2.5, 2.6, 2.7 under the following stronger assumption:

Assumption (SA'-r): We have (A'-r). Moreover we have, for a λ-integrable function J on E and a constant A:

b , b ′ , b , c , c , c ′ , J ≤ A, δ(ω, t, z) r ≤ J(z), δ(ω, t, z) 2 ≤ J(z), (3.1)
In the sequel we thus suppose that X satisfies (SA'-r), and also that (2.4) holds: these assumptions are typically not recalled. Below, all constants are denoted by K, and they vary from line to line. They may implicitly depend on the process X (usually through A in (3.1)). When they depend on an additional parameter p, we write K p .

We will usually replace the discontinuous process X by the continuous process

X ′ t = t 0 b ′ s ds + t 0 σ s dW s , (3.2)
connected with X by X t = X 0 + X ′ t + s≤t ∆X s . Note that b ′ is bounded, and without loss of generality we will use below its càdlàg version. Note also that, since the jumps of c are bounded, one can rewrite (2.2) as

c t = c 0 + t 0 b s ds + t 0 σ s dW s + t 0 E δ(s, z) (µ -ν)(ds, dz). (3.3) This amounts to replacing b in (2.2) by b t+ + E δ(t+, z)(κ( δ(t+, z) )-1 { δ(t+,z) ≤1} ) λ(dz),
where κ is a continuous function with compact support, equal to 1 on the set [0, A]. Note that the new process b is bounded càdlàg.

With any process Z we associate the variables

η(Z) t,s = E sup v∈(t,t+s] Z t+v -Z t 2 | F t , (3.4)
and we recall Lemma 4.2 of [5]:

Lemma 3.1 For all t > 0, all bounded càdlàg processes Z, and all sequences v n ≥ 0 of reals tending to 0, we have ∆ n E

[t/∆n] i=1 η(Z) (i-1)∆n,vn → 0, and for all 0 ≤ v ≤ s we have E(η(Z) n t+v,s | F t ) ≤ η(Z) t,s .

An Auxiliary Result on Itô Semimartingales

In this subsection we give some simple estimates for a d-dimensional semimartingale

Y t = t 0 b Y s ds + t 0 σ Y s dW s + t 0 E δ Y (s, z) (µ -ν)(ds, dz)
on some space (Ω, F, (F t ) t≥0 , P), which may be different from the one on which X is defined, as well as W and µ, but we still suppose that the intensity measure ν is the same. Note that Y 0 = 0 here. We assume that for some constant A and function J Y we have, with

c Y = σ Y σ Y, * : b Y ≤ A, c Y ≤ A 2 , δ Y (ω, t, z) 2 ≤ J Y (z) ≤ A 2 , E J Y (z) λ(dz) ≤ A 2 . (3.5)
The compensator of the quadratic variation of Y is of the form

t 0 c Y s ds, where c Y t = c Y t + E δ Y (t, z) δ Y (t, z) * λ(dz).
E(Y t | F 0 ) -tb Y 0 ≤ t η(b Y ) 0,t ≤ Kt E(Y j t Y m t | F 0 ) -tc Y,jm 0 ≤ Kt(t + √ t η(b Y ) 0,t + η(c Y ) 0,t ) ≤ Kt, (3.6
)

and if further E(c Y t -c Y 0 | F 0 ) ≤ A 2
t for all t, we also have

E(Y j t Y m t | F 0 ) -tc Y,jm 0 ≤ 2 t 3/2 (2A 2 √ t + Aη(b Y ) 0,t ) ≤ Kt 3/2 . (3.7) b) When Y is continuous, and if E( c Y t -c Y 0 2 | F 0 ) ≤ A 4 t for all t, we have E Y j t Y k t Y l t Y m t | F 0 -t 2 (c Y,jk 0 c Y,lm 0 + c Y,jl 0 c Y,km 0 + c Y,jm 0 c Y,kl 0 ) ≤ Kt 5/2 . ( 3
E (Y j t Y k t -tc Y,jk 0 )(c Y,lm t -c Y,lm 0 ) | F 0 ≤ Kt 3/2 ( √ t + η( c ′Y ) 0,t ). (3.9)
Proof. The first part of (3.6) follows by taking the F 0 -conditional expectation in the decomposition

Y t = M t + tb Y 0 + t 0 (b Y s -b Y 0 ) ds,
where M is a d-dimensional martingale with M 0 = 0. For the second part, we deduce from Itô's formula that Y j Y m is the sum of a martingale vanishing at 0 and of

b j 0 t 0 Y m s ds+b m 0 t 0 Y j s ds+ t 0 Y m s (b j s -b j 0 ) ds+ t 0 Y j s (b m s -b m 0 ) ds+c Y,jm 0 t+ t 0 (c Y,jm s -c Y,jm 0 ) ds. Since E( Y t | F 0 ) ≤ KA √ t,
as in (3.10), we deduce the second part of (3.6) and also (3.7) by taking again the conditional expectation and by using the Cauchy-Schwarz inequality and the first part.

(3.8) is a part of Lemma 5.1 of [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF]. For (3.9), we first observe that

Y j t Y k t -tc Y,jk 0 = B t +M t and c Y,lm t -c Y,lm 0 = B ′ t + M ′ t ,
|b s | ≤ 2 Y s b Y s + c Y s -c Y 0 , |b ′ s | ≤ K, |h s | ≤ K Y s 2 , |h ′ s | ≤ K, whereas h ′′ s = Y j s c ′Y,k,lm + Y k s c ′Y,j;lm . As seen before, E( Y t q | F 0 ) ≤ K q t q/2 for all q, and E( c Y t -c Y 0 2 | F 0 ) ≤ Kt. This yields E(B 2 t | F 0 ) ≤ Kt 3 and E(M 2 t | F 0 ) ≤ Kt 2 . Since |B ′ t | ≤ Kt and E(M ′2 t | F 0 ) ≤ Kt, we deduce that the F 0 -conditional expectations of B t B ′ t and B t M ′ t and M t B ′ t are smaller than Kt 2 . Finally E(M t M ′ t | F 0 ) = E( M, M ′ t | F 0 ), and M, M ′ t is the sum of c ′Y,k,lm 0 t 0 Y j s ds + t 0 Y j s ( c ′Y,k,lm s -c ′Y,k,lm 0
) ds and a similar term with k and j exchanged. Then using again

E( Y t 2 | F 0 ) ≤ Kt, plus E(Y t | F 0 )
≤ Kt and Cauchy-Schwarz inequality, we obtain that the above conditional expectation is smaller than K(t 2 + t 3/2 η( c ′Y ) t ). This completes the proof of (3.9).

Some Estimates.

1) We begin with well known estimates for X ′ and c, under (3.1) and for s, t ≥ 0 and q ≥ 0:

E sup w∈[0,s] X ′ t+w -X ′ t q | F t ≤ K q s q/2 , E(X ′ t+s -X ′ t | F s ) ≤ Ks E sup w∈[0,s] c t+w -c t q | F t ≤ K q s 1∧(q/2) , E(c t+s -c t | F s ) ≤ Ks. (3.10)
Next, it is much easier (although unfeasible in practice) to replace c n i in (2.6) by the estimators based on the process X ′ given by (3.2). Namely, we will replace c n i by the following:

c ′n i = 1 k n ∆ n kn-1 j=0 ∆ n i+j X ′ ∆ n i+j X ′ * .
The difference between c n i and c ′n i is estimated by the following inequality, valid when u n ≍ ∆ ̟ n and q ≥ 1, and where a n denotes a sequence of numbers (depending on u n ), going to 0 as n → ∞ (this is Equation 4.8 of [START_REF] Jacod | Quarticity and other functionals of volatility: efficient estimation[END_REF]):

E c n i -c ′n i q ≤ K q a n ∆ (2q-r)̟+1-q n . (3.11)
2) The jumps of c also potentially cause troubles. So we will eliminate the "big" jumps as follows. For any ρ > 0 we consider the subset E ρ = {z : J(z) > ρ}, which satisfies λ(E ρ ) < ∞, and we denote by G ρ the σ-field generated by the variables µ([0, t] × A), where t ≥ 0 and A runs through all Borel subsets of E ρ . The process

N ρ t = µ((0, t] × E ρ ) (3.12)
is a Poisson process and we let S ρ 1 , S ρ 2 , • • • be its successive jump times, and Ω n,t,ρ be the set on which S ρ j / ∈ {i∆ n : i ≥ 1} for all j ≥ 1 such that S ρ j < t, and S ρ j+1 > t ∧ S ρ j + (6k n + 1)∆ n for all j ≥ 0 (with the convention S ρ 0 = 0; taking 6k n here instead of the more natural k n will be needed in the proof of Theorem 2.6, and makes no difference here). All these objects are G ρ -measurable, and P(Ω n,t,ρ ) → 1 as n → ∞, for all t, ρ > 0.

We define the processes By Lemma 2.1.5 and Proposition 2.1.10 in [START_REF] Jacod | Discretization of Processes[END_REF] applied to each components of X ′ and c (2) (ρ), plus the property b(ρ) ≤ K/ρ, for all t ≥ 0, s ∈ [0, 1], ρ ∈ (0, 1], q ≥ 2, we have

b(ρ) t = b t - Eρ δ(t+, z) λ(dz), c(ρ) t = σ t σ * t + (Eρ) c δ(t+, z) δ(t+, z) * λ(dz) c(ρ) t = c t - t 0 Eρ δ(s, z) µ(ds, dz) = c (1) (ρ) t + c (2) (ρ) t , where c (1) (ρ) t = c 0 + t 0 b(ρ) s ds + t 0 σ s dW s c (2) (ρ) t = t 0 (Eρ) c δ(t-, z) (µ -ν)(ds, dz), (3.13) so c(ρ), which is R d 2 ⊗ R d 2 -valued,
E sup w∈[0,s] X ′ t+w -X ′ t q | F t ∨ G ρ ≤ K q s q/2 E(X ′ t+s -X ′ t | F s ∨ G ρ ) + E(c(ρ) t+s -c(ρ) t | F s ∨ G ρ ) ≤ Ks E sup w∈[0,s] c (2) (ρ) t+w -c (2) (ρ) t q | F t ∨ G ρ ≤ K q φ ρ (s + s q/2 ) E sup w∈[0,s] c(ρ) t+w -c(ρ) t q | F t ∨ G ρ ≤ K q φ ρ s + s q/2 + s q ρ q ≤ K q,ρ s. (3.14)
where φ ρ = (Eρ) c J(z) λ(dz) → 0 as ρ → 0. Note also that b(ρ) t ≤ K/ρ.

3) For convenience, we put

b n i = b (i-1)∆n , c n i = c (i-1)∆n b(ρ) n i = b(ρ) (i-1)∆n , c(ρ) n i = c(ρ) (i-1)∆n , c(ρ) n i = c(ρ) (i-1)∆n F n i = F (i-1)∆n , F n,ρ i = F n i ∨ G ρ .
(3.15)

All the above variables are F n,ρ i -measurable. Recalling (3.4), and writing η(Z, (H t )) t,s if we use the filtration (H t ) instead of (F t ), we also set

η(ρ) n i,j = max(η(Y, (G ρ F t )) (i-1)∆n,j∆n : Y = b ′ , b(ρ), c, c(ρ), c ′ , η(ρ) n i = η(ρ) n i,i+2kn .
Therefore, Lemma 3.1 yields for all t, ρ > 0 and j, k such that j + k ≤ 2k n :

∆ n E [t/∆n] i=1 η(ρ) n i → 0, E(η(ρ) n i+j,k | F n,ρ i ) ≤ η(ρ) n i . (3.16)
We still need some additional notation. First, define G ρ -measurable (random) set of integers:

L(n, ρ) = {i = 1, 2, • • • : N ρ (i+2kn)∆n -N ρ (i-1
)∆n = 0} (3.17) (taking above, 2k n instead of k n , is necessary for the proof of Theorem 2.5). Observe that

i ∈ L(n, ρ), 0 ≤ j ≤ 2k n + 1 ⇒ c n i+j -c n i = c(ρ) n i+j -c(ρ) n i . (3.18)
Second, we define the following R d ⊗ R d -valued variables

α n i = ∆ n i X ′ ∆ n i X ′ * -c n i ∆ n β n i = c ′n i -c n i = 1 kn∆n kn-1 j=0 α n i+j + (c n i+j -c n i )∆ n γ n i = c ′n i+kn -c ′n i = β n i+kn -β n i + c n i+kn -c n i . (3.19)
4) Now we proceed with estimates. (3.14) yields, for all q ≥ 0:

E( α n i q | F n,ρ i ) ≤ K q ∆ q n , E(α n i | F n,ρ i ) ≤ K∆ 3/2 n E kn-1 j=0 α n i+j q | F n,ρ i ≤ K q ∆ 3q/4 n , E( c ′n i q | F n,ρ i ) ≤ K q (3.20)
the third inequality following from the two first one, plus Burkholder-Gundy and Hölder inequalities, and the last inequality form the third one and the boundedness of c t . Moreover, since the set {i ∈ L(n, ρ)} is G ρ -measurable, the last part of (3.14), (3.18), and Hölder's inequality, readily yield

q ≥ 2, i ∈ L(n, ρ) ⇒ E β n i q | F n,ρ i ) ≤ K q ∆ n φ ρ + ∆ q/4 n + ∆ q/2 n ρ q . (3.21)
5) The previous estimates are not enough for us. We will apply the estimates of Lemma 3.2 with Y t = X ′ (i-1)∆n+t -X ′ (i-1)∆n for any given pair n, i, and with the filtration (F (i-1)∆n+t ∨ G ρ ) t≥0 . We observe that on the set A(ρ, n, i) = {∃j ≤ 2k n : ij ∈ L(n, ρ)}, which is G ρmeasurable, and because of (3.18), the process c Y coincide with c(ρ

) (i-1)∆n +t -c(ρ) (i-1)∆n if t ∈ [0, ∆ n ].
Then in restriction to this set, by (3.7) and (3.8) and by the definition of η(ρ) n i,1 , we have

E(∆ n i X ′j ∆ n i X ′m | F n,ρ i ) -c n,jm i ∆ n ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i,1 ) E ∆ n i X ′j ∆ n i X ′k ∆ n i X ′l ∆ n i X ′m | F n,ρ i -(c n,jk i c n,lm i + c n,jl i c n,km i + c n,jm i c n,kl i )∆ 2 n ≤ K ρ ∆ 5/2 n
(the constant above depends on ρ, through the bound K/ρ for the drift of c(ρ)). Then a simple calculation gives us

E(α n i | F n,ρ i ) ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i,1 ) E α n,jk i α n,lm i | F n,ρ i ) -(c n,jl i c n,km i + c n,jm i c n,kl i )∆ 2 n ≤ K ρ ∆ 5/2 n on A(ρ, n, i) (3.22)
Next, we apply Lemma 3.2 to the process Y t = c(ρ) (i-1)∆n+tc(ρ) (i-1)∆n for any given pair n, i, and with the filtration (F (i-1)∆n+t ∨ G ρ ) t≥0 . We then deduce from (3.6), plus again (3.18), that

i ∈ L(n, ρ), 0 ≤ t ≤ k n ∆ n ⇒ E((c jk (i-1)∆n+t -c jk (i-1)∆n )(c lm (i-1)∆n +t -c lm (i-1)∆n ) | F n,ρ i ) -tc(ρ) n,jklm i ≤ K ρ t η(ρ) n i,kn E(c (i-1)∆n+t -c (i-1)∆n | F n,ρ i ) -t b(ρ) n i ≤ K ρ t η(ρ) n i,kn ≤ K p t.
(3.23) Moreover, the Cauchy-Schwarz inequality and (3.20) on the one hand, and (3.9) applied with the process Y t = X ′ (i-1)∆n+t -X ′ (i-1)∆n on the other hand, give us

i ∈ L(n, ρ) ⇒ E α n,kl i ∆ n i b(ρ) ms | F n,ρ i ≤ K∆ n η(ρ) n i,1 E α n,kl i ∆ n i c ms | F n,ρ i ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i,1
).

(3.24)

6)

We now proceed to estimates on β n i :

Lemma 3.3 We have on the set where i belongs to L(n, ρ):

E(β n,jk i β n,lm i | F n,ρ i ) -1 kn (c n,jl i c n,km i + c n,jm i c n,kl i ) -kn∆n 3 c(ρ) n,jklm i ≤ K ρ √ ∆ n (∆ 1/4 n + η(ρ) n i ) E(β n,jk i (c n,lm i+kn -c n,lm i ) | F n,ρ i ) -kn∆n 2 c(ρ) n,jklm i ≤ K ρ √ ∆ n ( √ ∆ n + η(ρ) n i ).
Proof. We set

ζ n i,j = α n i+j + (c n i+j -c n i )∆ n and write β n,jk i β n,lm i as 1 k 2 n ∆ 2 n kn-1 u=0 ζ n,jk i,u ζ n,lm i,u + 1 k 2 n ∆ 2 n kn-2 u=0 kn-1 v=u+1 ζ n,jk i,u ζ n,lm i,v + 1 k 2 n ∆ 2 n kn-2 u=0 kn-1 v=u+1 ζ n,lm i,u ζ n,jk i,v . (3.25)
For the estimates below, we implicitly assume i ∈ L(n, ρ) and u, v

∈ {0, • • • , k n -1}.
First, we deduce from (3.22) and (3.23), plus (3.24) and successive conditioning, that

E(ζ n,jk i,u ζ n,lm i,u | F n,ρ i ) -(c n,jl i c n,km i + c n,jm i c n,kl i )∆ 2 n ≤ K∆ 5/2 n . (3.26)
Second, if u < v, the same type of arguments and the boundedness of b(ρ)) t and c t yield

|E(ζ n,jk i,v | F n,ρ i+u+1 ) -(c n,jk i+u+1 -c n,jk i )∆ n -b(ρ) n,jk i+u+1 ∆ 2 n (v -u -1)| ≤ K∆ 3/2 n (k n √ ∆ n + η(ρ) n i+v,1 ) |E(α n,lm i+u (c n,jk i+u+1 -c n,jk i+u ) | F n,ρ i+u )| ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i+u,1 ) |E(α n,lm i+u (c n,jk i+u -c n,jk i ) | F n,ρ i+u )| ≤ K∆ 3/2 n ( √ ∆ n + η n i+u,1 ) |E(α n,lm i+u ( b(ρ) n,jk i+u+1 -b(ρ) n,jk i+u ) | F n,ρ i+u )| ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i+u,1 ) |E(α n,lm i+u b(ρ) n,jk i+u | F n,ρ i+u )| ≤ K ρ ∆ 3/2 n ( √ ∆ n + η(ρ) n i+u,1 ) |E((c n,lm i+u -c n,lm i ) (c n,jk i+u+1 -c n,jk i ) | F n,ρ i ) -c(ρ) n,jklm i ∆ n u| ≤ K ρ ∆ n η(ρ) n i |E((c n,lm i+u -c n,lm i ) b(ρ) n,jk i+u+1 | F n,ρ i )| ≤ K ρ ∆ 1/4 n . Since kn-2 u=0 kn-1 v=u+1 u = k 3 n /6+ O(k 2 n )
, we easily deduce that the F n,ρ i -conditional expectation of the last term in (3.25) is 1 6 c(ρ) n,jklm i k n ∆ n , up to a remainder term which is

O( √ ∆ n (∆ 1/4 n + η(ρ) n i ))
, and the same is obviously true of the second term. The first claim of the lemma readily follows from this and (3.25) and (3.26).

The proof of the second claim is similar. Indeed, we have

β n,jk i (c n,lm i+kn -c n,lm i ) = 1 k n ∆ n kn-1 u=0 α n,jk i,u + (c n,jk i+u -c n,jk i )∆ n c n,lm i+kn -c n,lm i ) and E(c n,lm i+kn -c n,lm i | F n,ρ i+u+1 ) -c n,lm i+u+1 -c n,lm i -b(ρ) n,lm i+u+1 ∆ n (k n -u -1) ≤ K∆ n η(ρ) n i+u+1,kn-u .
Using the previous estimates, we conclude as for the first claim. 2

Finally, we deduce the following two estimates on the variables γ n i of (3.19), for any q ≥ 2:

i ∈ L(n, ρ) ⇒      E γ n,jk i γ n,lm i | F n,ρ i ) -2 kn (c n,jl i c n,km i + c n,jm i c n,kl i ) -2kn∆n 3 c(ρ) n,jklm i ≤ K ρ √ ∆ n ∆ 1/8 n + η(ρ) n i E( γ n i q | F n,ρ i ) ≤ K q √ ∆ n φ ρ + ∆ q/4 n + ∆ q/2 n ρ q . (3.27)
To see that the first claim holds, one expands the product γ n,jk i γ n,lm i and use successive conditioning, the Cauchy-Schwarz inequality and (3.14), (3.18) and(3.23), and Lemma 3.3; the contributing terms are

β n,jk i β n,lm i + β n,jk i+kn β n,lm i+kn + (c n,jk i+kn -c n,jk i )(c n,lm i+kn -c n,lm i ) -β n,jk i (c n,lm i+kn -c n,lm i ) -β n,lm i (c n,jk i+kn -c n,jk i ).
For the second claim we use (3.14), (3.18) and (3.21), and it holds for all q ≥ 2.

The Behavior of Some Functionals of c(ρ).

For ρ > 0 we set

U (ρ) n t = [t/kn∆n]-3 j=3 µ(ρ) n j 2 1 { µ(ρ) n j >u ′ n /4} , where µ(ρ) n j = 1 kn kn-1 w=0 (c(ρ) n jkn+w -c(ρ) n (j-2)kn+w ).
(3.28)

The aim of this subsection is to prove the following lemma:

Lemma 3.4 Under (SA'-r) and ( 2.22) we have

lim ρ→0 lim sup n→∞ E U (ρ) n t = 0.
Assumption (SA'-r) is of course not fully used. What is needed is the assumptions concerning the process c t only.

Proof. With the notation (3.13), and for l = 1, 2 we define µ (l) (ρ) n j and U (l) (ρ) n t as above, upon substituting c(ρ) and u ′ n /4 with c (l) (ρ) and u ′ n /8. Since U (ρ) n t ≤ 4U (1) (ρ) n t +4U (2) (ρ) n t , it suffices to prove the result for each U (l) (ρ) n t . First, µ (1)

(ρ) n j 2 1 { µ (1) (ρ) n i >u ′ n /8} is smaller than K µ (1) (ρ) n j 4 /u ′2 n , whereas (recalling b(ρ) ≤ K/ρ) classical estimates yield E µ (1) (ρ) n j 4 ≤ K∆ n (1 + ∆ n /ρ). Thus the expec- tation of U (1) (ρ) n t is less than K∆ 1/2-2̟ ′ n
(1 + ∆ n /ρ), yielding the result for U (1) (ρ) n t . Secondly, we have U (2)

(ρ) n t ≤ [t/kn∆n] j=3 µ (2) (ρ) n i 2
and the first part of (3.14) yields

E µ (2) (ρ) n i 2 ≤ Kφ ρ √ ∆ n . Since φ ρ → 0 as ρ → 0, the result for U (1) (ρ) n t follows. 2
3.5 A Basic Decomposition.

We start the proof of Theorem 2.2 by giving a decomposition of V (g) n -V (g), with quite a few terms. It is based on the key property c ′n i = c n i + β n i and on the definition (3.19) of α n i and

β n i . A simple calculation shows that 1 √ ∆n (V ′ (g) n t -V (g) t ) = 5 j=1
V n,j t , as soon as t > k n ∆ n , where (the sums on components below always extend from 1 to d):

V n,1 t = √ ∆ n [t/∆n]-kn+1 i=1 g( c n i ) -g( c ′n i ) V n,2 t = 1 √ ∆n [t/∆n]-kn+1 i=1 i∆n (i-1)∆n (g(c n i ) -g(c s )) ds V n,3 t = 1 kn √ ∆n [t/∆n]-kn+1 i=1 l,m ∂ lm g(c n i ) kn-1 u=0 α n,lm i+u V n,4 t = √ ∆n kn [t/∆n]-kn+1 i=1 l,m ∂ lm g(c n i ) kn-1 u=1 (c n,lm i+u -c n,lm i ) -1 √ ∆n t ∆n([t/∆n]-kn+1 g(c s ) ds V n,5 t = √ ∆ n [t/∆n]-kn+1 i=1 g(c n i + β n i ) -g(c n i ) - l,m ∂ lm g(c n i ) β n,lm i .
The leading term is V n,3 , the bias comes from the terms V n,4 and V n,5 , and the first two terms are negligible, in the sense that they satisfy

j = 1, 2 ⇒ V n,j t P -→ 0 for all t > 0. (3.29)
We end this subsection with the proof of (3.29).

The case j = 1: (2.7) implies

|g( c n i )-g( c ′n i )| ≤ K(1+ c n i + c ′n i ) p-1 c n i -c ′n i ≤ K(1+ c ′n i ) p-1 c n i -c ′n i +K c n i -c ′n i p .
Recalling the last part of (3.20), we deduce from (3.11), from the fact that 1

-r̟ -p(1 - 2̟) < (2-r)̟ 2q
for all q > 1 small enough, and from Hölder's inequality, that E(|g(

c n i ) - g( c ′n i )|) ≤ Ka n ∆ (2p-r)̟+1-p n
. Therefore

E sup s≤t |V n,1 s | ≤ Kta n ∆ (2p-r)̟+1/2-p
n and (3.29) for j = 1 follows.

The case j = 2: Since g is C 2 and c t is an Itô semimartingale with bounded characteristics, the convergence V n,2 u.c.p. =⇒ 0 is well known: see for example the proof of (5.3.24) in [START_REF] Jacod | Discretization of Processes[END_REF], in which one replaces ρ cs (f ) by g(c s ).

3.6

The Leading Term V n,3 .

Our aim here is to prove that V n,3 L-s =⇒ Z (3.30) (functional stable convergence in law), where Z is the process defined in Theorem 2.2.

A change of order of summation allows us to rewrite V n,3 as

V n,3 t = 1 √ ∆ n [t/∆n] i=1 l,m w n,lm i α n,lm i , where w n,lm i = 1 k n (i-1)∧(kn-1) j=(i-[t/∆n]+kn-1) + ∂ lm g(c n i-j ).
Observe that w n i and α n i are measurable with respect to F n i and F n i+1 , respectively, so by Theorem IX.7.28 of [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF] (with G = 0 and Z = 0 in the notation of that theorem) it suffices to prove the following four convergences in probability, for all t > 0 and all component indices:

1 √ ∆ n [t/∆n]-kn+1 i=1 w n,lm i E(α n,lm i | F n i ) P -→ 0 (3.31) 1 ∆ n [t/∆n]-kn+1 i=1 w n,jk i w n,lm i E(α n,jk i α n,lm i | F n i ) P -→ t 0 ∂ jk g(c s ) ∂ lm g(c s ) c jl s c km s + c jm s c kl s ds (3.32) 1 ∆ 2 n [t/∆n]-kn+1 i=1 w n i 4 E( α n i 4 | F n i ) P -→ 0 (3.33) 1 √ ∆ n [t/∆n]-kn+1 i=1 w n,lm i E(α n,lm i ∆ n i N | F n i ) P -→ 0, (3.34)
where N = W j for some j, or is an arbitrary bounded martingale, orthogonal to W .

For proving these properties, we pick a ρ bigger than the upper bound of the function J, so G ρ becomes the trivial σ-field and F n i = F n,ρ i and L(n, ρ) = N. In such a way, we can apply all estimates of the previous subsections with the conditioning σ-fields F n i . Therefore (3.20) and the property w n i ≤ K readily imply (3.31) and (3.33). In view of the form of α n i , a usual argument (see e.g. [START_REF] Jacod | Discretization of Processes[END_REF]) shows that in fact E(α n,lm i In view of the definition of w n i , for each t we have w n,jk i(n,t) → ∂ jk g(c t ) and c n,jk i(n,t) → c jk t almost surely if |i(n, t)∆ n -t| ≤ k n ∆ n (recall that c is almost surely continuous at t, for any fixed t), and the above convergence follows by the dominated convergence theorem, thus ending the proof of (3.30).

∆ n i N | F n i ) =

3.7

The Term V n, 4 .

In this subsection we prove that, for all t, V n,

(i -1)∆ n < s ≤ i∆ n ⇒ H(n, t) lm s = √ ∆ n k n kn-2 u=0 (k n -1 -u)∂ lm g(c n i-u )
(its expression on [0, k n ∆ n ) and on (tk n ∆ n , t] is more complicated, but not needed, apart from the fact that it is uniformly bounded). Now, since kn-2 u=0 (k n -1u) = k 2 n /2+ O(k n ) as n → ∞, we observe that H(n, t) lm s converges to θ 2 ∂ lm g(c s-) for all s ∈ (0, t). Since c is a.s. continuous at t, we deduce from the dominated convergence theorem for stochastic integrals that V ′n,4 t indeed converges in probability to the first term in the right side of (3.35).

3.8

The Term V n,5 .

The aim of this subsection is to prove the convergence

V n,5 t P -→ A 2 t -2A 3 t + θ s≤t 1 0 g(c s-+ w∆c s ) -g(c s-) -w l,m ∂ lm g(c s-) ∆c lm s dw (3.36) We have V n,5 t = [t/∆n]-kn+1 i=1
v n i , where

v n i = ∆ n g(c n i + β n i ) -g(c n i ) - l,m ∂ lm g(c n i ) β n,lm i .
We also set

α n i = 1 kn∆n kn-1 u=0 α n i+u , β n i = β n i -α n i = 1 kn kn-1 u=1 (c n i+u -c n i ) v ′n i = √ ∆ n g(c n i + β n i ) -g(c n i ) - l,m ∂ lm g(c n i ) β n,lm i , v ′′n i = v n i -v ′n i . (3.37)
We take ρ ∈ (0, 1], and will eventually let it go to 0. With the sets L(n, ρ) of (3.17), we associate

L(n, ρ, t) = {1, • • • , [t/∆ n ] -k n + 1} ∩ L(n, ρ) L(n, ρ, t) = {1, • • • , [t/∆ n ] -k n + 1}\L(n, ρ).
We split the sum giving V n,5 t into three terms:

U n,ρ t = i∈L(n,ρ,t) v n i , U ′n,ρ t = i∈L(n,ρ,t) v ′n i , U ′′n,ρ t = i∈L(n,ρ,t) v ′′n i . ( 3

.38)

A) The processes U n,ρ . A Taylor expansion and (2.7) give us

v n i = v(1) n i +v(2) n i +v(3) n i , where      v(1) n i = √ ∆n 2 j,k,l,m ∂ 2 jk,lm g(c n i ) E(β n,jk i β n,lm i | F n,ρ i ) v(2) n i = √ ∆n 2 j,k,l,m ∂ 2 jk,lm g(c n i ) β n,jk i β n,lm i -v(1) n i |v(3) n i | ≤ K √ ∆ n (1 + β n i ) p-3 β n i 3 . Therefore U n,ρ = 3 j=1 U (j) n,ρ , where U (j) n t = i∈L(n,ρ,t) v(j) n i . (3.39)
On the one hand, and letting

w(ρ) n i = j,k,l,m ∂ 2 jk,lm g(c n i ) 1 2k n √ ∆ n (c n,jl i c n,km i + c n,jm i c n,kl i ) + k n √ ∆ n 6 c(ρ) n,jklm i , the càdlàg property of c and c(ρ) and k n √ ∆ n → θ imply W (ρ) n t := ∆ n [t/∆n]-kn+1 i=1 w(ρ) n i P -→ U (1) ρ t := A 2 t + θ 6 j,k,l,m t 0 ∂ 2 jk,lm g(c s ) c(ρ) jklm s ds.
On the other hand, Lemma 3.

3 yields [v(1) n i -∆ n w(ρ) n i | ≤ K ρ ∆ n (∆ 1/4 n + η(ρ) n i ) when i ∈ L(n, ρ), whereas |w(ρ) n i | ≤ K always. Therefore E |U (1) n,ρ t -W (ρ) n t | ≤ K ρ ∆ n E [t/∆n] i=1 ( ∆ n + η(ρ) n i ) + K∆ n E #(L(n, ρ, t)) . Now, #(L(n, ρ, t)) is not bigger than (2k n + 1)N ρ t , implying that ∆ n E(#(L ′ (n, ρ, t))) ≤ K ρ √ ∆ n .
Taking advantage of (3.16), we deduce that the above expectation goes to 0 as n → ∞, and thus

U (1) n,ρ t P -→ U (1) ρ t . (3.40) Next, v(2) n i is F n,ρ i+kn -measurable, with vanishing F n,ρ i -conditional expectation, and each set {i ∈ L(n, ρ)} is F n,ρ 0 -measurable. It follows that E (U (2) n,ρ t ) 2 ≤ 2k n E i∈L(n,ρ,t) E |v(2) n i | 2 | F n,ρ i ≤ Kk n ∆ n E i∈L(n,ρ,t) E |β n i | 4 | F n,ρ i ≤ Ktφ ρ + K ρ t √ ∆ n ,
where we have applied (3.21) for the last inequality. Another application of the same estimate gives us

E |U (3) n t |) ≤ Ktφ ρ + K ρ t∆ 1/4 n .
These two results and the property φ ρ → 0 as ρ → 0 clearly imply

lim ρ→0 lim sup n→∞ E(|U (2) n,ρ t | + |U (3) n,ρ t |) = 0. (3.41)
B) The processes U ′n,ρ . We will use here the jump times S ρ 1 , S ρ 2 , • • • of the Poisson process N ρ , and will restrict our attention to the set Ω n,t,ρ defined before (3.13), whose probability goes to 1 as n → ∞. On this set, L(n, ρ, t) is the collection of all integers i which are between [S ρ q /∆ n ] -2k n + 2 and [S ρ q /∆ n ] + 1, for some q between 1 and N ρ t . Thus

U ′n,ρ t = N ρ t q=1 H(n, ρ, q), where H(n, ρ, q) = [S ρ q /∆n]+1 i=[S ρ q /∆n]-2kn+1 v ′n i . (3.42)
The behavior of each H(n, ρ, q) is a pathwise question. We fix q and set S = S ρ q and a n = [S/∆ n ], so S > a n ∆ n because S is not a multiple of ∆ n . For further reference we consider a case slightly more general than strictly needed here. We have c n i → c S-when a n -6k n + 1 ≤ i ≤ a n + 1 and c n i → c S when a n + 2 ≤ i ≤ a n + 6k n , uniformly in i (for each given outcome ω). Hence

β n i - (k n -a n + i -2) + ∧ (k n -1) k n ∆c S → 0 uniformly in i ∈ {a n -6k n + 2, • • • , a n + 5k n }.
(3.43) Thus, the following convergence holds, uniform in

i ∈ {a n -2k n + 1, • • • , a n + 1}: 1 √ ∆n v ′n i -g c S-+ kn-an+i-2 kn ∆c S -g(c S-) -l,m ∂ lm g(c S -) c jk S-+ kn-an+i-2 kn ∆c lm S → 0, which implies H(n, ρ, q) -∆ n kn-3 u=1 g c Sq-+ u k n ∆c Sq -g(c Sq -) - l,m ∂ lm g(c Sq -) u k n ∆c lm Sq → 0
and by Riemann integration this yields H(n, ρ, q) → θ ). On the one hand, V n,5 = U (1) n,ρ + U (2) n,ρ + U (3) n,ρ + U ′n,ρ + U ′′n,ρ ; on the other hand, the dominated convergence theorem (observe that c(ρ) t → σ 2 t for all t) yields that U (1) ρ E) Proof of Theorem 2.2. We are now ready to prove Theorem 2.2. Recall that 1 √ ∆n (V (g)n t -V (g)) = 5 j=1 V n,j . By virtue of (3.29), (3.30), (3.35), (3.36), it is enough to check that 3.9 Proof of Theorem 2.5.

A 1 t + A 3 t + A 4 t + A 5 t = θ
The proof of Theorem 2.5 follows the same line as in Subsection 3.8, and we begin with an auxiliary step.

Step 1) Replacing c n i by c ′n i . The summands in the definition (2.17 We prove here that A n,3 t -A ′n,3 t P -→ 0 (3.46) for all t, and this is done as in to the step j = 1 in Subsection 3.5. The function R is C 1 on R 2 + with ∂ j R(x, y) ≤ K(1 + x + y ) p-j for j = 0, 1, by (2.7). Thus Then, exactly as in the case afore-mentioned, we conclude (3.46), and it remains to prove that, for all t, we have

A ′n,3 t P -→ - 1 2 A 2 t + A 3 t + A ′4 t .
Step 2) From now on we use the same notation as in Subsection 3.8, although they denote different variables or processes. For any ρ ∈ (0, 1] we have A ′n,3 = U n,ρ + U ′n,ρ + U ′′n,ρ , as defined in (3.38), but with

v n i =

[

  s and c t = σ 2 t is the squared volatility. Natural estimators for t 0 g(c s ) ds are V n (g) t = ∆ n

  [t/∆n] P -→ c t for any t, and there is an associated Central Limit Theorem under (A-2), with rate 1/∆ 1/4

  .21) at least when the function F on M + d × M d , where M d is the set of all d × d matrices, is C 1 and |F (x, y)| ≤ K y 2 uniformly in x within any compact set, as is the function G ′′ above.

  Moreover, if the process c Y is itself an Itô semimartingale, the quadratic covariation of the continuous martingale parts of Y and c Y is also of the form t 0 c ′Y s ds for some process c ′Y , necessarily bounded if both Y and c Y satisfy (3.5) (and, if Y = X, we have c Y = c and c ′Y = c ′ ).

Lemma 3 . 2

 32 Below we assume(3.5), and the constant K only depends on A. a) We have for t ∈ [0, 1]:

  with M and M ′ martingales (M is continuous). The processes B, B ′ , M, M , M ′ , M ′ and M, M ′ are absolutely continuous, with densities b s , b ′ s , h s , h ′ s and h ′′ s satisfying, by (3.5) for Y and c Y :

  is the càdlàg version of the density of the predictable quadratic variation of c(ρ). Moreover G ρ = {∅, Ω} and ( b(ρ), c(ρ)) = ( b, c) when ρ exceeds the bound of the function J. Note also that b(ρ) and c(ρ) are càdlàg.

[ 0 ∂

 0 0 for all N as above, hence (3.34) holds.For (3.32), by(3.22) it suffices to prove that ∆ n jk g(c s ) ∂ lm g(c s ) c jl s c km s +c jm s c kl s ds.

1 0g 4 n

 14 (c Sq-+ w∆c Sq )g(c Sq -)w l,m ∂ lm g(c Sq -) ∆c lm S Sq -+w∆c Sq )-g(c Sq -)-w l,m ∂ lm g(c Sq -) ∆c lm Sq dw. (3.44) C) The processes U ′′n,ρ . Since |β n i | ≤ K we deduce from (2.7) that |v ′′n i | ≤ K √ ∆ n ( α n i + α n i p ). (3.20) yields E( α n i q | F n,ρ i ) ≤ K q ∆q/for all q > 0. ThereforeE |U ′′n,ρ t | ≤ K∆ 3/4 n E(#(L(, n, ρ, t))) ≤ K ρ ∆ 1/4 n ,by virtue of what precedes (3.40). We then deduce U

  s-+ w∆c s )g(c s-)w l,m ∂ lm g(c s-) ∆c lm s dw as ρ → 0 (for the latter convergence, note that |g(x + y)g(x)l,m ∂ lm g(x)y lm | ≤ K y 2 when x, y stay in a compact set). Then the property (3.36) follows from (3.40), (3.41), (3.44) and (3.45).

2 l,m t 0 ∂ 3 t + θ s≤t 1 0,m t 0 ∂ 1 0 w dw = 1 2 )

 031012 lm g(c s-) dc lm sθ g(c t ) -2A g(c s-+ w∆c s )g(c s-)w l,m ∂ lm g(c s-) ∆c lm s dw.To this aim, we observe that Itô's formula gives usg(c t ) = g(c 0 )+ llm g(c s-) dc lm ss-+∆c s )-g(c s-)l,m∂ lm g(c s-) ∆c lm s , so the desired equality is immediate (use also , and the proof of Theorem 2.2 is complete.

) of A n, 3 t

 3 are R( c n i , c n i+kn ), where R(x, y) = j,k,l,m ∂ 2jk,lm g(x)(y jkx jk )(y lmx lm ), and we set A

|R( c n

  i , c n i+kn ) -R( c ′n i , c ′n i+kn )| ≤ K(1 + c ′n i + c ′n i+kn ) ) p-1 ( c n ic ′n i + c n i+knc ′n i+kn ) +K c n ic ′n i p + K c n i+knc ′n i+kn p .

  .8) c) When c Y is a (possibly discontinuous) semimartingale satisfying the same conditions (3.5) as Y , and if Y itself is continuous, we have

  is a d × d-dimensional predictable process, bounded uniformly (in n, s, ω) and given on the set [k n ∆ n , tk n ∆ n ] by

	We first observe that c n i+u -c n i = summation, we can rewrite V ′n,4 t as	u-1 v=0 ∆ n i+v c. Then, upon changing the order of
	V n,4 t	=	[t/∆n]-1 i=1	l,m	w n,lm i	∆ n i c lm ,		w n,lm i	=	√ ∆ n k n	i-1)∧(kn-2) u=0∨(i+kn-1-[t/∆n] (k n -1 -u)∂ lm g(c n i-u ).
	In other words, recalling k n	√ ∆ n ≤ K and ∂g(c s ) ≤ K, we see that
								V n,4 t	=	l,m	0	t	H(n, t) lm s dc lm s ,
	where H(n, t)
						t	4	P -→	θ 2 l,m	0	t	∂ lm g(c s-) dc lm s -θ g(c t ).	(3.35)
	We call V ′n,4 t of V n,4 t . Since k n and V ′′n,4 t √ ∆ n → θ and c is a.s. continuous at t, it is obvious that V ′′n,4 , respectively, the first sum, and the last integral, in the definition t converges t to the first term almost surely to -θ g(c t ), and it remains to prove the convergence of V ′n,4 in the right side of (3.35).

s

  -

	√ i = -v ′n √ ∆n 8 R(c n i + β n i , c n i+kn + β n i+kn ) ∆n 8 R(c n i + β n i , c n i+kn + β n i+kn ), v ′′n i = v n i -v ′n i .
	Recalling γ n i in (3.19), the decomposition (3.39) holds with
	v(1) n i = -v(2) n i = -v(3) n i = v n √ √ i -v(1) n ∆n 8 j,l,k,m ∂ 2 jl,km g(c n i ) E γ n,jk i ∆n 8 j,l,k,m ∂ 2 jl,km g(c n i ) γ n,jk i γ n,lm γ n,lm i i -v(1) n | F n,ρ i i i -v(2) n i .
	Use c ′n i -c n i = β n i and (2.7) and a Taylor expansion to check that
	|v(3) n i | ≤ K ∆ n γ n i i | ≤ K We also have |v(2) n √ ∆ n γ n i 2 , hence (3.21) and (3.27) yield 2 β n i (1 + β n i ) p-3 .
	E(|v(3) n i | | G n +	∆ n ρ p ,
	and thus (3.41) holds here as well, by the same argument. Moreover, (3.27) again yields
	(3.40), with now
	U ρ t = -	j,k,l,m	0	t	∂ 2 jk,lm g(c s )	θ 12	c(ρ) jklm s	+	1 4θ	(c jl s c km s + c jm s c kl s ) ds.
	This goes to A 3 t -1 2 A 2 t as ρ → 0. Another application of (2.7) gives us
	|v ′′n i | ≤ K ∆ n 1 + γ n i	2	α n i + α n i+kn + α n i	p + α n i+kn	p .
	Then another application of (3.20), (3.21) and (3.27) yields E(|v ′′n i | | G ρ ) ≤ K∆ conclude (3.45) as previously. We are thus left to prove that	3/4 n and we
	ρ > 0 ⇒ U ′n,ρ t	P -→ U ′ρ

ρ ) + E(|v(2) n i | 2 | G ρ ) ≤ K∆ n φ ρ + ∆ 1/4 t , with, as ρ → 0, U ′ρ t P -→ A ′4 t .

(3.47)

Step 3) On the set Ω n,t,ρ we have (3.42) and we study H(n, ρ, q), in the same way as before, on the set Ω n,t,ρ . We fix q and set S = S q and a n = [S/∆ n ]. We then apply (3.43) and also c n i → c S-or c n i → c S , according to whether a n -2k n +1 ≤ i ≤ a n +1 or a n +2 ≤ i ≤ a n +k n , to obtain v ′n iv ′n i → 0, uniformly in i between a n -2k n + 1 and a n + 1, where

We then deduce, by Riemann integration, that

which is θG ′ (c Sq-, ∆c Sq ), hence the first part of (3.47), with

. The second part of (3.47) follows from the dominated convergence theorem, and the proof of Theorem 2.5 is complete.

Proof of Theorem 2.6.

The proof is once more somewhat similar to the proof of Subsection 3.8, although the way we replace c n i by c ′n i and further by α n i + β n i is different.

A) Preliminaries. The jth summand in (2.24) involves several estimators c n i , spanning the time interval ((j -3)k n ∆ n , (j + 2)k n ∆ n ]. It is thus convenient to replace the sets L(n, ρ), L(n, ρ, t) and L(n, ρ, t), for ρ, t > 0, by the following ones:

For any ρ ∈ (0, 1] we write V(F

We also set

(3.11) and the last part of (3.20) yield

Observe that δ n j c ′ is analogous to γ n i , with a doubled time lag, so it satisfies a version of (3.27) and, for q ≥ 2, we have

B) The processes V n,ρ . (2.25) yields

Thus a (tedious) computation shows that, with the notation

, where

First, (3.48) and Hólder's inequality give us E(a ′n j ) ≤ K q,v ∆ l(q,v) n

for any q > 1 and v > 0, where (recalling (2.8) and (2.22) for ̟ and ̟ ′ ) we have set l(q, v) = 1-r̟ q p(1 -2̟) ∨ v(1 -2̟ + ̟ ′ ) . Upon choosing v small enough and q close enough to 1, and in view of (2.8), we see that l(q, v) > 1 2 , thus implying

Next, we deduce from (3.49) that

and thus, since p ≥ 3, lim

The analysis of B n,ρ t is more complicated. We have δ n j c ′ = z n j + z ′n j , where

It easily follows that for all

On the one hand, (3.20) and Hölder's inequality yield E(a(1

On the other hand, observe that z ′n j = µ(ρ) n j , with the notation (3.28), and as soon as j ∈ L ′ (n, ρ), so Lemma 3.4 gives us

Finally, (3.14) shows that E( z ′n

∆ n for all q ≥ 2 and j ∈ L ′ (n, ρ), whereas c ′n (j-3)kn+1 is F n (j-2)kn+1 -measurable, so (3.14), (3.20) and successive conditioning yield E(a(3

Then, again as for (3.53), one obtains C) The processes V n,ρ . With the previous notation S ρ j and N ρ t , and on the set Ω n,ρ,t , we have converge in probability to c S-if j < 0 and to c S if j > 0, whereas for j = 0 we have

By virtue of the definition of v n j , and since u ′ n → 0 and also since w n is almost surely in (0, 1) and F is continuous and F (x, 0) = 0, one readily deduces that In view of (2.25), an application of the dominated convergence theorem gives V ρ t → V(F ) t . Then (2.26) follows from V(F ) n t = V n,ρ t + V n,ρ t and (3.56) and (3.58), and the proof of Theorem 2.6 is complete.


Nicolas Le Roux
email: nicolas@le-roux.name

Francis Bach
email: francis.bach@ens.fr

Local Component Analysis

Kernel density estimation, a.k.a. Parzen windows, is a popular density estimation method, which can be used for outlier detection or clustering. With multivariate data, its performance is heavily reliant on the metric used within the kernel. Most earlier work has focused on learning only the bandwidth of the kernel (i.e., a scalar multiplicative factor). In this paper, we propose to learn a full Euclidean metric through an expectation-minimisation (EM) procedure, which can be seen as an unsupervised counterpart to neighbourhood component analysis (NCA). In order to avoid overfitting with a fully nonparametric density estimator in high dimensions, we also consider a semi-parametric Gaussian-Parzen density model, where some of the variables are modelled through a jointly Gaussian density, while others are modelled through Parzen windows. For these two models, EM leads to simple closed-form updates based on matrix inversions and eigenvalue decompositions. We show empirically that our method leads to density estimators with higher test-likelihoods than natural competing methods, and that the metrics may be used within most unsupervised learning techniques that rely on local distances, such as spectral clustering or manifold learning methods. Finally, we present a stochastic approximation scheme which allows for the use of this method in a large-scale setting.

Introduction

Most unsupervised learning methods rely on a metric on the space of observations. The quality of the metric directly impacts the performance of such techniques and a significant amount of work has been dedicated to learning this metric from data when some supervised information is available [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF][START_REF] Goldberger | Neighbourhood components analysis[END_REF][START_REF] Bach | Learning spectral clustering[END_REF]. However, in a fully unsupervised scenario, most practitioners use the Ma-halanobis distance obtained from principal component analysis (PCA). This is an unsatisfactory solution as PCA is essentially a global linear dimension reduction method, while most unsupervised learning techniques, such as spectral clustering or manifold learning, are local.

In this paper, we cast the unsupervised metric learning as a density estimation problem with a Parzen windows estimator based on a Euclidean metric. Using the maximum likelihood framework, we derive in Section 3 an expectation-minimisation (EM) procedure that maximizes the leave-one-out log-likelihood, which may be considered as an unsupervised counterpart to neighbourhood component analysis (NCA) [START_REF] Goldberger | Neighbourhood components analysis[END_REF]. As opposed to PCA, which performs a whitening of the data based on global information, our new algorithm globally performs a whitening of the data using only local information, hence the denomination local component analysis (LCA).

Like all non-parametric density estimators, Parzen windows density estimation is known to overfit in high dimensions [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], and thus LCA should also overfit. In order to keep the modelling flexibility of our density estimator while avoiding overfitting, we propose a semiparametric Parzen-Gaussian model; following [START_REF] Blanchard | In search of non-Gaussian components of a high-dimensional distribution[END_REF], we linearly transform then split our variables in two parts, one which is modelled through a Parzen windows estimator (where we assume the interesting part of the data lies), and one which is modelled as a multivariate Gaussian (where we assume the noise lies). Again, in Section 4, an EM procedure for estimating the linear transform may be naturally derived and leads to simple closed-form updates based on matrix inversions and eigenvalue decompositions. This procedure contains no hyperparameters, all the parameters being learnt from data.

Since the EM formulation of LCA scales quadratically in the number of datapoints, making it impractical for large datasets, we introduce in Section 5 both a stochastic approximation and a subsampling technique allowing us to achieve a linear cost and thus to scale LCA to much larger datasets.

Finally, in Section 6, we show empirically that our method leads to density estimators with higher testlikelihoods than natural competing methods, and that the metrics may be used within unsupervised learning techniques that rely on such metrics, like spectral clustering.

Previous work

Many authors aimed at learning a Mahalanobis distance suited for local learning. While some techniques required the presence of labelled data [START_REF] Goldberger | Neighbourhood components analysis[END_REF][START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF][START_REF] Bach | Learning spectral clustering[END_REF], others proposed ways to learn the metric in a purely unsupervised way, e.g., [START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF] who used the distance to the k-th nearest neighbour as the local scaling around each datapoint. Most of the other attempts at unsupervised metric learning were developed in the context of kernel density estimation, a.k.a. Parzen windows. The Parzen windows estimator [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] is a nonparametric density estimation model which, given n datapoints {x 1 , . . . ,

x n } in R d , defines a mixture model of the form p(x) = 1 n n j=1 K(x, x j , θ)
where K is a kernel with compact support and parameters θ. We relax the compact support assumption and choose K to be the normal kernel, that is

p(x) = 1 n n j=1 N (x, x j , Σ) ∝ 1 n |Σ| n j=1 exp - 1 2 (x -x j ) ⊤ Σ -1 (x -x j ) ,
where Σ is the covariance matrix of each Gaussian.

As the performance of the Parzen windows estimator is more reliant on the covariance matrix than on the kernel, there has been a large body of work, originating from the statistics literature, attempting to learn this matrix. However, almost all attempts are focused on the asymptotic optimality of the estimators obtained with little consideration for the practicality in high dimensions. Thus, the vast majority of the work is limited to isotropic matrices, reducing the problem to finding a single scalar h [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF][START_REF] Duin | On the choice of smoothing parameters for Parzen estimators of probability density functions[END_REF][START_REF] Rudemo | Empirical choice of histograms and kernel density estimators[END_REF][START_REF] Chow | Consistent cross-validated density estimation[END_REF][START_REF] Bowman | An alternative method of crossvalidation for the smoothing of density estimates[END_REF][START_REF] Park | Comparison of datadriven bandwidth selectors[END_REF][START_REF] Sheather | A reliable databased bandwidth selection method for kernel density estimation[END_REF], the bandwidth, and the few extensions to the non-isotropic cases are numerically expensive [START_REF] Duong | Cross-validation bandwidth matrices for multivariate kernel density estimation[END_REF][START_REF] Jones | Maximum likelihood kernel density estimation: On the potential of convolution sieves[END_REF].

An exception is the approach proposed in [START_REF] Vincent | Manifold Parzen windows[END_REF], which is very similar to our method, as the authors learn the covariance matrix of the Parzen windows estimator using local neighbourhoods. However, their algorithm does not minimize a well-defined cost function, making it unsuitable for kernels other than the Gaussian one, and the locality used to compute the covariance matrix depends on parameters which must be handtuned or cross-validated. Also, the modelling of all the dimensions using the Parzen windows estimator makes the algorithm unsuitable when the data lie on a high-dimensional manifold. In an extension to [START_REF] Vincent | Manifold Parzen windows[END_REF], [START_REF] Bengio | Nonlocal manifold Parzen windows[END_REF] uses a neural network to compute the leading eigenvectors of the local covariance matrix at every point in the space, then uses these matrices to do density estimation and classification. Despite the algorithm's impressive performance, it does not correspond to a linear reparametrisation of the space and thus cannot be used as a preprocessing step.

Local Component Analysis

Seeing the density as a mixture of Gaussians, one can easily optimize the covariances using the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. However, maximizing the standard loglikelihood of the data would trivially lead to the degenerate solution where Σ goes to 0 to yield a sum of Dirac distributions. One solution to that problem is to penalize some norm of the precision matrix to prevent it from going to infinity. Another, more compelling, way is to optimize the leave-one-out log-likelihood, where the probability of each datapoint x i is computed under the distribution obtained when x i has been removed from the training set. This technique is not new and has already been explored both in the supervised [START_REF] Goldberger | Neighbourhood components analysis[END_REF][START_REF] Globerson | Metric learning by collapsing classes[END_REF] and in the unsupervised setting [START_REF] Duin | On the choice of smoothing parameters for Parzen estimators of probability density functions[END_REF]. However, in the latter case, the cross-validation was then done by hand, which explains why only one bandwidth parameter could be optimized1 . We will thus use the following criterion:

L(Σ) = - n i=1 log 1 n -1 j =i N (x i , x j , Σ) (1)
≤ cst - n i=1 j =i λ ij log N (x i , x j , Σ) + n i=1 j =i λ ij log λ ij , (2)
with the constraints ∀i , j =i λ ij = 1. This variational bound is obtained using Jensen's inequality.

The EM algorithm optimizes the right-hand side of Eq. ( 2) by alternating between the optimisations of λ and Σ in turn. The algorithm is guaranteed to converge, and does so to a stationary point of the true function over Σ defined in Eq. [START_REF] Bach | Diffrac: a discriminative and flexible framework for clustering[END_REF]. At each step, the optimal solutions are:

λ * ij = N (x i , x j , Σ) k =i N (x i , x k , Σ) if j = i (3) λ * ii = 0 (4) Σ * = ij λ ij (x i -x j )(x i -x j ) T n . (5)
The "responsibilities" λ * ij define the relative proximity of x j to x i (compared to the proximity of all the x k 's to x i ) and Σ * is the average of all the local covariance matrices.

This algorithm, which we coin LCA, for local component analysis, transforms the data to make it locally isotropic, as opposed to PCA which makes it globally isotropic. Fig. [START_REF] Bach | Diffrac: a discriminative and flexible framework for clustering[END_REF] shows a comparison of PCA and LCA on the word sequence "To be or not to be". Whereas PCA is highly sensitive to the length of the text, LCA is only affected by the local shapes, thus providing a much less distorted result. 2First, one may note at this point that Manifold Parzen Windows [START_REF] Vincent | Manifold Parzen windows[END_REF] is equivalent to LCA with only one step of EM. This makes Manifold Parzen Windows more sensitive to the choice of the original covariance matrix whose parameters must be carefully chosen. As we shall see later in the experiments, running EM to convergence is important to get good accuracy when using spectral clustering on the transformed data. Second, it is also worth noting that, similarly to Manifold Parzen Windows, LCA can straightforwardly be extended to the cases where each datapoint uses its own local covariance matrix (possibly with a smoothing term), or where the covariance Σ * is the sum of a low-rank matrix and some scalar multiplied by the identity matrix.

Not only may LCA be used to learn a linear transformation of the space, but it also defines a density model. However, there are two potential negative aspects associated with this method. First, in high dimensions, Parzen windows is prone to overfitting and must be regularized [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. Second, if there are some directions containing a small Gaussian noise, the local isotropy will blow them up, swamping the data with clutter. This is common to all the techniques which renormalise the data by the inverse of some variance. A solution to both of these issues is to consider a product of two densities: one is a low-dimensional Parzen windows estimator, which will model the interesting signal, and the other is a Gaussian, which will model the noise.

LCA with a multiplicative Gaussian component

We now assume that there are irrelevant dimensions in our data which can be modelled by a Gaussian. In other words, we consider an invertible linear transformation (B G , B L ) ⊤ x of the data, modelling B ⊤ G x as a multivariate Gaussian and B ⊤ L x through kernel density estimation, the two parts being independent, leading to p

(x) ∝ p(B ⊤ G x, B ⊤ L x) = p G (B ⊤ G x)p L (B ⊤ L x)
, where p G is a Gaussian and p L is the Parzen windows estimator, i.e.,

p(x i ) ∝ B G B ⊤ G + B L B ⊤ L 1 2 n -1 × exp - 1 2 (x i -µ) ⊤ B G B ⊤ G (x i -µ) ×   j =i exp - 1 2 (x i -x j ) ⊤ B L B ⊤ L (x i -x j )   ,
with (B G , B L ) a full-rank square matrix. Using EM, we can upper-bound the negative log-likelihood:

-2 i log p(x i ) ≤ tr(B ⊤ G C G B G ) + tr(B ⊤ L C L B L ) -log |B G B ⊤ G + B L B ⊤ L | , (6)
with

C G = 1 n i (x i -µ)(x i -µ) ⊤ , C L = 1 n ij λ ij (x i -x j )(x i -x j ) ⊤ .
The matrices B G and B L minimizing the right-hand side of Eq. ( 6) may be found using the following proposition (see proof in the appendix):

Proposition 1 Let B G ∈ R d×d1 and B L ∈ R d×d2 , with d = d 1 + d 2 and B = (B G , B L ) ∈ R d×d invertible. Consider two symmetric positive matrices M 1 and M 2 in R d×d . The problem min BG,BL tr B ⊤ G M 1 B G +tr B ⊤ L M 2 B L -log det(B G B ⊤ G +B L B ⊤ L ) (7)
has a finite solution only if M 1 and M 2 are invertible and, if these conditions are met, reaches its optimum at

B G = M -1/2 1 U + , B L = M -1/2 1 U -D -1/2 - , where U + are the eigenvectors of M -1/2 1 M 2 M -1/2 1
associated with eigenvalues greater than or equal to 1, U -are the eigenvectors of

M -1/2 1 M 2 M -1/2 1 associated Figure 1:
Results obtained when transforming "To be or not to be" using PCA (left) and LCA (right). Left: To make the data globally isotropic, PCA awkwardly compresses the letters horizontally. Right: Since LCA is insensitive to the spacing between letters and to the length of the text, there is no horizontal compression.

with eigenvalues smaller than 1 and D -is the diagonal matrix containing the eigenvalues of

M -1/2 1 M 2 M -1/2 1 smaller than 1.
The resulting procedure is described in Algorithm 1, where all dimensions are initially modelled by the Parzen windows estimator, which empirically yielded the best results.

Algorithm 1 LCA -Gauss Input: X (dataset), iterMax (maximum number of iterations), ν (regularisation) Output: B G (Gaussian part transformation), B L (Parzen windows transformation) C G ← cov(X) + νI d {Initialize C to the global covariance} I G ← C -1 2 G B G = 0, B L = chol(C -1 G
) {Assign all dimensions to the Parzen windows estimator} for iter = 1:iterMax do

M ij ← exp - (xi-xj) ⊤ B ⊤ L BL(xi-xj ) 2 , M ii ← 0 λ ij ← Mij k M ik C L ← ij λij (xi-xj)(xi-xj ) ⊤ n + νI d [V, D] ← eig(I G C L I G ) {Eigendecomposition of I G C L I G } t 1 = max z D(z, z) ≤ 1 {Cut-off between eigenval- ues smaller and larger than 1} t + = {t|t 1 ≤ t ≤ d} , t -= {t|1 ≤ t < t 1 } V + ← V (:, t + ), V -← V (:, t -), D -= D(t -, t -) B L = I G V -D -1/2 - , B G = I G V + end for
Relationship with ICA. Independent component analysis (ICA) can be seen as a density model where x = As and s has independent components (see, e.g., [START_REF] Hyvärinen | Independent component analysis[END_REF]). In the Parzen windows framework, this corresponds to modelling the density of s by a product of univariate kernel density estimators [START_REF] Boscolo | Independent component analysis based on nonparametric density estimation[END_REF]. This however causes two problems: first, while this assumption is appropriate in settings such as source separation, it is violated in most settings, and having a multivariate kernel density estimation is preferable. Second, most algorithms are dedicated to finding independent components which are non-Gaussian. In the presence of more than one Gaussian dimension, most ICA frameworks become unidentifiable, while our explicit modelling of such Gaussian components allows us to tackle this situation (a detailed analysis of the identifiability of our Parzen/Gaussian model is out of the scope of this paper).

Relationship with NGCA. NGCA [START_REF] Blanchard | In search of non-Gaussian components of a high-dimensional distribution[END_REF] makes an assumption similar to ours (they rather assume an additive Gaussian noise on top of a low-dimensional non-Gaussian signal) but uses a projection pursuit algorithm to iteratively find the directions of non-Gaussianity. Unlike in FastICA, the contrast functions used to find the interesting directions can be different for each direction. However, like all projection pursuit algorithms, the identification of interesting directions gets much harder in higher dimensions, as most of them will be almost Gaussian. Our use of a nonparametric density estimator with a log-likelihood cost allows us to globally optimize all directions simultaneously and does not rely on the model being correct. Finally, LCA estimates all its parameters from data as opposed to NGCA which requires the number of non-Gaussian directions to be set.

Escaping local optima. Though our model allows for the modification of the number of dimensions modelled by the Gaussian through the analysis of the spectrum of

C -1/2 G C L C -1/2 G
, it is sensitive to local optima. It is for instance rare that a dimension modelled by a Gaussian is switched to the Parzen windows estimator. Even though the algorithm will more easily switch from the Parzen windows estimator to the Gaussian model, it will typically stop too early, that is model many dimensions using the Parzen windows estimator rather than the better Gaussian. To solve these issues, we propose an alternate algorithm, LCA-Gauss-Red, which explores the space of dimensions modelled by a Gaussian more aggressively using a search algorithm, namely:

1. We run the algorithm LCA -Gauss for a few iterations (40 in our experiments);

2. We then "transfer" some columns from B L (the Parzen windows model) to B G (the Gaussian model), and rerun LCA -Gauss using these new matrices as initialisations;

3. We iterate step 2 using a dichotomic search of the optimal number of dimensions modelled by the Gaussian, until a local optimum is found;

4. Once we have a locally optimum number of dimensions modelled by the Gaussian model, we run LCA -Gauss to convergence.

Speeding up LCA

Computing the local covariance matrix of the points using Eq. ( 3), ( 4) and ( 5) has a complexity in O(dn 2 + d 2 n + d 3 ), with d the dimensionality of the data and n the number of training points. Since this is impractical for large datasets, we can resort to sampling to keep the cost linear in the number of datapoints. We may further use low-rank or diagonal approximation to achieve a complexity which grows quadratically with d instead of cubically.

Averaging a subset of the local covariance matrices

Instead of averaging the local covariances over all datapoints, we may only average them over a subset of datapoints. This estimator is unbiased and, if the local covariance matrices are not too dissimilar, which is the assumption underlying LCA, then its variance should remain small. This is equivalent to using a minibatch procedure: every time we have a new minibatch of size B, we compute its local covariance C L , which is then averaged with the previously computed C L using

C L ← γ B n C L + (1 -γ B n ) C L (8)
to yield the updated C L . The exponent B/n is so that γ, the discount factor, determines the weight of the old covariance matrix after an entire pass through the data, which makes it insensitive to the particular choice of batch size. As opposed to many such algorithms where the choice of γ is critical as it helps retaining the information of previous batches, the locality of the EM estimate makes it less so. However, if the number of datapoints used to estimate C L is not much larger than the dimension of the data, we need to set a higher γ to avoid degenerate covariance matrices. In simulations, we found that using a value of γ = .6 worked well. Similarly, the size of the minibatch influences only marginally the final result and we found a value of 100 to be large enough.

Computing the local covariance matrices using a subset of the datapoints

Rather than using only a subset of local covariance matrices, one may also wonder if using the entire dataset to compute these matrices is necessary. Also, as the number of datapoints grows, the chances of overfitting increase. Thus, one may choose to use only a subset of the datapoints to compute these matrices. This will increase the local covariances, yielding a biased estimate of the final result, but may also act as a regulariser. In practice, for very large datasets, one will want the largest neighbourhood size while keeping the computational cost tractable.

Denoting n i the number of locations at which we estimate the local covariance and n j the number of neighbours used to estimate this covariance, the cost per update is now O(d 2 [n i + n j ] + dn i n j + d 3 ). Since only n j should grow with n, this is linear in the total number of datapoints.

Though they may appear similar, these are not "landmark" techniques (see, e.g., [START_REF] Silva | Sparse multidimensional scaling using landmark points[END_REF]) as there is still one Gaussian component per datapoint, and the n i datapoints around which we compute the local covariances are randomly sampled at every iteration.

Experiments

LCA has three main properties: first, it transforms the data to make it locally isotropic, thus being well-suited for preprocessing the data before using a clustering algorithm like spectral clustering; second, it extracts relevant, non-Gaussian components in the data; third, it provides us with a good density model through the use of the Parzen windows estimator.

In the experiments, we will assess the performance of the following algorithms: LCA, the original algorithm; LCA-Gauss, using a multiplicative Gaussian component, as described in Section 4; LCA-Gauss-Red, the variant of LCA-Gauss using the more aggressive search to find a better number of dimensions to be modelled by the Gaussian component. The MATLAB code for LCA, LCA-Gauss and LCA-Gauss-Red is available at http://nicolas.le-roux.name/code.html.

Improving clustering methods

We first try to solve three clustering problems: one for which the clusters are convex and the direction of interest does not have a Gaussian marginal (Fig. 2), right). Following [START_REF] Bach | Diffrac: a discriminative and flexible framework for clustering[END_REF][START_REF] Bach | Learning spectral clustering[END_REF], the data is progressively corrupted by adding dimensions of white Gaussian noise, then whitened. We compare here the clustering accuracy, which is defined as 100 n min P tr(EP ) where E is the confusion matrix and P is the set of permutations over cluster labels, obtained with the following five techniques:

1. Spectral clustering (SC) [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] on the whitened data (using the code of [START_REF] Chen | Parallel spectral clustering in distributed systems[END_REF]);

2. SC on the projection on the first two components found by FastICA using the best contrast function and the correct number of components;

3. SC on the data transformed using the metric learnt with LCA;

4. SC on the data transformed using the metric learnt with the product of LCA and a Gaussian;

5. SC on the projection of the data found using NGCA [START_REF] Blanchard | In search of non-Gaussian components of a high-dimensional distribution[END_REF] with the correct number of components.

Our choice of spectral clustering stems from its higher clustering performance compared to K-means. Results are reported in Fig. [START_REF] Bengio | Nonlocal manifold Parzen windows[END_REF]. Because of the whitening, the Gaussian components in the first dataset are shrunk along the direction containing information. As a result, even with little noise added, the information gets swamped and spectral clustering fails completely. On the other hand, LCA and its variants are much more robust to the presence of irrelevant dimensions.

Though NGCA works very well on the first dataset, where there is only one relevant component, its performance drops quickly when there are two relevant components (note that, for all datasets, we provided the true number of relevant dimensions as input to NGCA). This is possibly due to the deflation procedure which is not adapted when no single component can be clearly identified in isolation. This is in contrast with LCA and its variants which circumvent this issue, thanks to their global optimisation procedure. Note also that LCA-Gauss allows us to perform unsupervised dimensionality reduction with the same performance as previously proposed supervised algorithms (e.g., [START_REF] Bach | Learning spectral clustering[END_REF]).

Figure [START_REF] Blanchard | In search of non-Gaussian components of a high-dimensional distribution[END_REF] shows the clustering accuracy on the three datasets for various numbers of EM iterations, one iteration corresponding to Manifold Parzen Windows [START_REF] Vincent | Manifold Parzen windows[END_REF] with a Gaussian kernel whose covariance matrix is the data covariance kernel. As one can see, running the EM algorithm to convergence yields a significant improvement in clustering accuracy. The performance of Manifold Parzen Windows could likely have been improved with a careful initialisation of the original kernel, but this would have been at the expense of the simplicity of the algorithm.

LCA as a density model

We now assess the quality of LCA as a density model. We build a density model of the USPS digits dataset, a 256-dimensional dataset of handwritten digits. We compared several algorithms:

• An isotropic Parzen windows estimator with the bandwidth estimated using LCA (replacing Σ * of Eq. ( 5) by λI so that the two matrices have the same trace);

• A Parzen windows estimator with diagonal metric (equal to the diagonal of Σ * in Eq. ( 5);

• A Parzen windows estimator with the full metric as obtained using LCA;

• A single Gaussian model;

• A product of a Gaussian and a Parzen windows estimator (as described in Section 4).

The models were trained on a set of 2000 datapoints and regularized by penalizing the trace of Σ -1 (in the case of the last model, both covariance matrices, local and global, were penalized). The regularisation parameter was optimized on a validation set of 1000 datapoints. For the last model, the regularisation parameter of the global covariance was set to the one yielding the best performance for the full Gaussian model on the validation set. Thus, we only had to optimize the regularisation parameter for the local covariance.

The final performance was then evaluated on a set of 3000 datapoints which had not been used for training nor validation. We ran the experiment 20 times, randomly selecting the training, validation and test set each time.

Fig. [START_REF] Borwein | Convex analysis and nonlinear optimization, theory and examples[END_REF] shows the mean and the standard error of the negative log-likelihood on the test set. As one can see, modelling all dimensions using the Parzen windows estimator leads to poor performance in high dimensions, despite the regulariser and the leave-one-out criterion.

On the other hand, LCA-Gauss and LCA-Gauss-Red clearly outperform all the other models, justifying our choice of modelling some dimensions using a Gaussian. Also, as opposed to the previous experiments, there is no performance gain induced by the use of LCA-Gauss-Red as opposed to LCA-Gauss, which we believe stems from the fact that the switch from one model to the other is easier to make when there are plenty of dimensions to choose from. The poor performance of LCA-Full is a clear indication of the problems suffered by Parzen windows in high dimensions.

Subsampling

We now evaluate the loss in performance incurred by the use the subsampling procedure described in Section 5, both on the train and test negative log- likelihoods. For that purpose, we used the USPS digit recognition dataset, which contains 8298 datapoints in dimension 256, which we randomly split into a training set of n = 6000 datapoints, using the rest as the test set. We tested the following hyperparameters:

• Discount factor γ = 0.3, 0.6, 0.9 ,

• Batch size B = 1000, 3000, 6000 ,

• Neighbourhood size N = 1000, 3000, 6000 .

Fig. [START_REF] Boscolo | Independent component analysis based on nonparametric density estimation[END_REF] show the log-likelihood differences induced by the use of smaller batch sizes and neighbourhood sizes. For each set of hyperparameters, 20 experiments were run using different training and test sets, and the means and standard errors are reported. The results for γ = 0.9 were very similar and are not included due to space constraints.

Three observations may be made. First, reducing the batchsize has little effect, except when γ is small. Second, reducing the neighbourhood size has a regularizing effect at first but drastically hurts the performance if reduced too much. Third, the value of γ, the discount factor, has little influence, but larger values proved to yield more consistent test performance, at the expense of slower convergence. The consistency of these results shows that it is safe to use subsampling (with values of γ = 0.6, B = 100 and N = 3000, for instance) especially if the training set is very large.

Conclusion

Despite its importance, the learning of local or global metrics is usually an overseen step in many practical algorithms. We have proposed an extension of the general bandwidth selection problem to the multidimensional case, with a generalisation to the case where several components are Gaussian. Additionally, we proposed an approximate scheme suited to large datasets which allows to find a local optimum in linear time.

We believe LCA can be an important preprocessing tool for algorithms relying on local distances, such as manifold learning methods or many semi-supervised algorithms. Another use would be to cast LCA within the mean-shift algorithm, which finds the modes of the Parzen windows estimator, in the context of image segmentation [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF]. In the future, we would like to extend this model to the case where the metric is allowed to vary with the position in space, to account for more complex geometries in the dataset.

  [START_REF] Bach | Learning spectral clustering[END_REF], left), one for which the clusters are not convex (Fig.(2), middle), and one for which the directions of interest have almost Gaussian marginals (Fig.(

Figure 2 :

 2 Figure 2: Noise-free data used to assess the robustness of K-means to noise. Left: mixture of two isotropic Gaussians of unit variance and means [-3, 0] ⊤ and [3, 0] ⊤ . Centre: two concentric circles with radii 1 and 2, with added isotropic Gaussian noise of standard deviation .1. Right: mixture of five Gaussians. The centre cluster contains four times as many datapoints as the other ones.

Figure 3 :

 3 Figure 3: Average clustering accuracy (100% = perfect clustering, chance is 50% for the first two datasets and 20% for the last one) on 100 runs for varying number of dimensions of noise added. The error bars represent one standard error. Left: mixture of isotropic Gaussians presented in Fig. (2) (left). Centre: two concentric circles presented in Fig. (2) (centre). Right: mixture of five Gaussians presented in Fig. (2) (right).

Figure 4 :

 4 Figure 4: Average clustering accuracy (100% = perfect clustering, chance is 50% for the first two datasets and 20% for the last one) on 100 runs for varying number of dimensions of noise added and varying number of EM iterations in the LCA algorithm (MPW = one iteration). The error bars represent one standard error. Left: mixture of isotropic Gaussians presented in Fig. (2) (left). Centre: two concentric circles presented in Fig. (2) (centre). Right: mixture of five Gaussians presented in Fig. (2) (right).

Figure 5 :

 5 Figure 5: Test negative log-likelihood on the USPS digits dataset, averaged over 20 runs.

Most of the literature on estimating the covariance matrix discards the log-likelihood cost because of its sensitivity to outliers and prefers AMISE (see, e.g.,[START_REF] Duong | Cross-validation bandwidth matrices for multivariate kernel density estimation[END_REF]). However, in all our experiments, the number of datapoints was large enough so that LCA did not suffer from the presence of outliers.

Since both methods are insensitive to any linear reparametrisation of the data, we do not include the original data in the figure.

Acknowledgements

Nicolas Le Roux and Francis Bach are supported in part by the European Research Council (SIERRA-ERC-239993). We would also like to thank Warith Harchaoui for its valuable input.

Appendix

We prove here Proposition 1.

Proof If M 1 is singular, then the minimum value is -∞, because we can have

The reasoning is similar for M 2 .

We thus assume that M 1 and M 2 are invertible.

We consider the eigendecomposition of

= U Diag(e)U ⊤ , which corresponds to the generalized eigendecomposition of the pair (M 1 , M 2 ).

By taking derivatives with respect to A 1 , we get

with

By left-multiplying both sides of Eq. ( 9) by A ⊤ 2 , we obtain

By left-multiplying by A ⊤ 1 , we get

Thus, we now need to minimize with respect to A 2 the following cost function

Let s be the vector of singular values of A 2 , ordered in decreasing order and let the e i be ordered in increasing order. We have:

with equality if and only if the eigenvectors of A 2 A ⊤ 2 are aligned with the ones of Diag(e) (the -e i being also in decreasing order) (Theorem 1.2.1, [START_REF] Borwein | Convex analysis and nonlinear optimization, theory and examples[END_REF]).

Thus, we have A 2 A ⊤ 2 = diag(s) 2 with only d 2 nonzero elements in s. Let J 2 be the index of non zeroelements. We thus need to minimize

with optimum s 2 j = e -1 j and value:

Thus, we need to take J 2 corresponding to the smallest eigenvalues e j . If we also optimize with respect to d 2 , then J 2 must only contain the elements smaller than 1.