Serge Cohen
email: serge.cohen@math.univ-toulouse.fr

Alexander Lindner
email: a.lindner@tu-bs.de

A central limit theorem for the sample autocorrelations of a Lévy driven continuous time moving average process

Keywords: Bartlett's formula, continuous time moving average process, estimation of the Hurst index, fractional Lévy process, Lévy process, limit theorem, sample autocorrelation, sample autocovariance, sample mean

In this article we consider Lévy driven continuous time moving average processes observed on a lattice, which are stationary time series. We show asymptotic normality of the sample mean, the sample autocovariances and the sample autocorrelations. A comparison with the classical setting of discrete moving average time series shows that in the last case a correction term should be added to the classical Bartlett formula that yields the asymptotic variance. An application to the asymptotic normality of the estimator of the Hurst exponent of fractional Lévy processes is also deduced from these results.

Introduction

Statistical models are often written in a continuous time setting for theoretical reasons (e.g. diffusions). But if one wants to estimate the parameters of these models, one usually assumes only the observation of a discrete sample. At this point a very general question, the answer of which depends on the model chosen, is to know if the estimation should not have been performed with an underlying discrete model in the beginning. In this article we will consider this for moving average processes and we refer to the classical moving average time series models as a discrete counterpart of this continuous model.

To be more specific, let L = (L t ) t∈R be a two sided one-dimensional Lévy process, i.e. a stochastic process with independent and stationary increments, càdlàg sample paths and which satisfies L 0 = 0. Assume further that L has finite variance and expectation zero, and let f : R → R be in L 2 (R). Let µ ∈ R. Then the process (X t ) t∈R , given by

X t = µ + R f (t -s) dL s , t ∈ R, (1.1)
can be defined in the L 2 sense and is called a continuous time moving average process with mean µ and kernel function f , driven by L. See also [START_REF] Cohen | Springer Verlag collection Lévy Matters[END_REF] for more information on such processes, in particular fractional Lévy processes. The process (X t ) t∈R is then strictly stationary. Equation (1.1) is the natural continuous time analogue of discrete time moving average processes

X t = µ + i∈Z ψ t-i Z i , t ∈ Z, (1.2)
where (Z t ) t∈Z is an independent and identically distributed (i.i.d.) noise sequence with finite variance and expectation zero, and (ψ i ) i∈Z is a square summable sequence of real coefficients. The asymptotic behaviour of the sample mean and sample autocorrelation function of X t in (1.2) has been studied for various cases of noise sequences (Z i ) i∈Z , such as regularly varying noise (cf. Davis and Mikosch [START_REF] Davis | The sample autocorrelations of heavy-tailed processes with applications to ARCH[END_REF]), martingale difference sequences (cf. Hannan [START_REF] Hannan | The asymptotic distribution of serial covariances[END_REF]), or i.i.d. sequences with finite fourth moment or finite variance but more restrictive conditions on the decay of the sequence (ψ i ) i∈Z (cf. Section 7 of Brockwell and Davis [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]).

Another approach to obtain limit theorems for sample autocovariances is to prove strong mixing properties of the time series under consideration, and provided it has finite (4 + δ)-moment, use the corresponding central limit theorems (such as in Ibragimov and Linnik [START_REF] Ibragimov | Independent and Stationary Sequences of Random Variables[END_REF], Theorem 18.5.3). If even stronger strong mixing conditions hold, then existence of a fourth moment may be enough. Observe however that processes with long memory are often not strongly mixing, and in this paper we are aiming also at applications with respect to the fractional Lévy noise, which is not strongly mixing.

In this paper we shall study the asymptotic behaviour as n → ∞ of the sample mean

X n;∆ := n -1 n i=1 X i∆ , (1.3)
of the process (X t ) t∈R defined in (1.1) when sampled at (∆n) n∈N , where ∆ > 0 is fixed, and of its sample autocovariance and sample autocorrelation function

γ n;∆ (∆h) := n -1 n-h i=1 (X i∆ -X n;∆ )(X (i+h)∆ -X n;∆ ), h ∈ {0, . . . , n -1}, (1.4)
ρ n;∆ (∆h) := γ n;∆ (∆h)/ γ n;∆ (0), h ∈ {0, . . . , n -1}.

(1.5)

We write N = {0, 1, 2, . . .}. Under appropriate conditions on f and L, in particular assuming L to have finite fourth moment for the sample autocorrelation functions, it will be shown that X n;∆ and ( ρ n;∆ (∆), . . . , ρ n;∆ (h∆)) are asymptotically normal for each h ∈ N as n → ∞. This is similar to the case of discrete time moving average processes of the form (1.2) with i.i.d. noise, but unlike for those, the asymptotic variance of the sample autocorrelations of model (1.1) will turn out to be given by Bartlett's formula plus an extra term which depends explicitly on the fourth moment of L, and in general this extra term does not vanish. This also shows that the "naive" approach of trying to write the sampled process (X n∆ ) n∈Z as a discrete time moving average process as in (1.2) with i.i.d. noise does not work in general, since for such processes the asymptotic variance would be given by Bartlett's formula only. If µ = 0, then further natural estimators of the autocovariance and autocorrelation are given by

γ * n;∆ (∆h) := n -1 n i=1 X i∆ X (i+h)∆ , h ∈ {0, . . . , n -1}, (1.6)
ρ * n;∆ (∆h) := γ * n;∆ (∆h)/γ * n;∆ (0), h ∈ {0, . . . , n -1}, (1.7)
and the conditions we have to impose to get asymptotic normality of γ * n;∆ and ρ * n;∆ are less restrictive than those for γ n;∆ and ρ n;∆ . We will be particularly interested in the case when f decays like a polynomial, which is e.g. the case for fractional Lévy noises. For a given Lévy process with expectation zero and finite variance, and a parameter d ∈ (0, 1/2), the (moving average) fractional Lévy process (M t;d ) t∈R with Hurst parameter H := d + 1/2 is given by Marquardt [11]). A process also called fractional Lévy process was introduced before by Benassi et al. [START_REF] Benassi | On roughness indices for fractional fields[END_REF], where (x) + = max(x, 0) is replaced by an absolute value in (1.8),

M 1 t;d := 1 Γ(d + 1) ∞ -∞ (t -s) d + -(-s) d + dL s , t ∈ R (1.8) (cf.
M 2 t;d := ∞ -∞ |t -s| d -|s| d dL s , t ∈ R. (1.9)
Although both processes have different distributions, they enjoy similar properties.

For instance the sample paths of both versions are Hölder continuous, have the same pointwise Hölder exponent, and they are both locally self-similar (see [START_REF] Benassi | On roughness indices for fractional fields[END_REF] for the definition of this local property of their distributions). The corresponding fractional Lévy noises based on increments of length ∆ > 0 are given by

X i t = M i t;d -M i t-∆;d , t ∈ R i = 1, 2.
Hence the fractional Lévy noise is a Lévy driven moving average process with kernel function

f 1 d,∆ (s) = 1 Γ(d + 1) s d + -(s -∆) d + , s ∈ R, (1.10) or f 2 d,∆ (s) = |s| d -|s -∆| d , s ∈ R. (1.11)
While the kernel functions f i d,∆ , i ∈ {1, 2}, do not satisfy the assumptions we will impose for the theorems regarding the sample mean X n;∆ and the sample autocorrelation function ρ n;∆ , for d ∈ (0, 1/4) they do satisfy the assumptions we impose for the asymptotic behaviour of ρ * n;∆ , so that an asymptotically normal estimator of the autocorrelation and hence of the Hurst index can be obtained if d ∈ (0, 1/4). For general d ∈ (0, 1/2), one may take the differenced fractional Lévy noises M i t;d -2M i t-∆;d + M i t-2∆;d , t ∈ R, and our theorems give asymptotically normal estimators for the autocorrelation function of these processes. Please note that asymptotically normal estimators of the Hurst exponent for M 2 t;d are already described in [START_REF] Benassi | On roughness indices for fractional fields[END_REF] but they use fill-in observations of the sample paths

X 2 (k/2 n ) for k = 1, . . . , 2 n -1. If L is a Brownian motion, then X 1 d = CX 2
, where d = means equality in distribution for processes and C is a constant, is the fractional Brownian motion and it is self-similar. Except in this case, fractional Lévy processes are not self-similar and therefore observations on a grid k/2 n do not yield the same information as the time series X i (t), t ∈ Z.

The paper is organised as follows: in the next section we will derive asymptotic normality of the sample mean. Then, in Section 3 we will derive central limit theorems for the sample autocovariance γ n;∆ and the sample autocorrelation ρ n;∆ , as well as for the related estimators γ * n;∆ and ρ * n;∆ of (1.6) and (1.7). As a byproduct of the asymptotic normality, these quantities are consistent estimators of the autocovariance and autocorrelation. Section 4 presents an application of our results to the estimation of the parameters of fractional Lévy noises, where the underlying Hurst parameter is estimated. We also recall there that fractional Lévy noises are mixing in the ergodic-theoretic sense, and we prove that they fail to be strongly mixing.

Throughout the paper, unless indicated otherwise, L will be a Lévy process with mean zero and finite variance σ 2 = EL 2 1 , and X = (X t ) t∈R denotes the process defined in (1.1) with kernel f ∈ L 2 (R), f : R → R. Its autocovariance at lag h ∈ R will be denoted by

γ(h) = γ f (h) = Cov (X 0 , X h ) = σ 2 R f (-s)f (h -s) ds, (1.12)
where the last equation follows from the Itô isometry.

Let us set some notations used in the sequel. If v is vector or A a matrix the transposed is denoted by v ′ , respectively by A ′ . Convergence in distribution is denoted by d →. The function 1 A for a set A is one for x ∈ A, and vanishing elsewhere. The autocorrelation of X at lag h will be denoted by ρ(h) = ρ f (h) = γ(h)/γ(0).

Asymptotic normality of the sample mean

The sample mean X n;∆ of the moving average process X of (1.1) behaves like the sample mean of a discrete time moving average process with i.i.d. noise, in the sense that it is asymptotically normal with variance σ 2 ∞ k=-∞ γ(k∆), provided the latter is absolutely summable.

Theorem 2.1. Let L have zero mean and variance σ 2 , let µ ∈ R and ∆ > 0. Suppose that   F ∆ : [0, ∆] → [0, ∞], u → F ∆ (u) = ∞ j=-∞ |f (u + j∆)|   ∈ L 2 ([0, ∆]). (2.1)
Then ∞ j=-∞ |γ(∆j)| < ∞, ∞ j=-∞ γ(∆j) = σ 2 ∆ 0   ∞ j=-∞ f (u + j∆)   2 du, (2.2)
and the sample mean of X ∆ , . . . , X n∆ is asymptotically normal as n → ∞, more precisely

√ n X n;∆ d → N   µ, σ 2 ∆ 0   ∞ j=-∞ f (u + ∆j)   2 du   as n → ∞.
Remark 2.2. Throughout the paper the assumption that L has zero mean can be dropped very often. For instance, if

f ∈ L 1 (R) ∩ L 2 (R)
, then the assumption of zero mean of L presents no restriction, for in that case L ′ t := L t -tE(L 1 ), t ∈ R, defines another Lévy process with mean zero and the same variance, and it holds

X t = µ + E(L 1 ) R f (s) ds + R f (t -s) dL ′ s , t ∈ Z, (2.3
)

which has mean µ + E(L 1 ) R f (s) ds.
Proof. For simplicity in notation, assume that ∆ = 1, and write F = F 1 . Continue F periodically on R by setting

F (u) = ∞ j=-∞ |f (u + j)|, u ∈ R. Since |γ f (h)| ≤ σ 2 ∞ -∞ |f (-s)| |f (h -s)| ds
by (1.12), we have

1 σ 2 ∞ h=-∞ |γ f (h)| ≤ ∞ -∞ |f (-s)| ∞ h=-∞ |f (h -s)| ds = ∞ -∞ |f (s)|F (s) ds = ∞ j=-∞ 1 0 |f (s + j)|F (s) ds = 1 0 F (s) F (s) ds < ∞. (2.4)
The same calculation without the modulus gives (2.2). The proof for asymptotic normality is now much in the same spirit as for discrete time moving average processes, by reducing the problem to m-dependent sequences first and then applying an appropriate variant of Slutsky's theorem. By subtracting the mean we may assume without loss of generality that µ = 0. For m ∈ N, let f m := f 1 (-m,m) , and denote

X (m) t := R f m (s) dL s = t+m t-m f (t -s) dL s , t ∈ Z.

Observe that (X (m) t

) t∈Z is a (2m-1)-dependent sequence, i.e. (X (m) j

) j≤t and (X (m) j ) j≥t+2m are independent for each t ∈ Z. From the central limit theorem for strictly stationary (2m -1)-dependent sequences (cf. Theorem 6.4.2 in Brockwell and Davis [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]) we then obtain that

√ n X (m) n;1 = n -1/2 n t=1 X (m) t d → Y (m) , n → ∞, (2.5)
where Y (m) is a random variable such that

Y (m) d = N (0, v m ) with v m = 2m j=-2m γ fm (j). Since lim m→∞ γ fm (j) = γ f (j) for each j ∈ Z by (1.12), since |γ f 1 (-m,m) (j)|≤σ 2 ∞ -∞ |f (-s)| |f (j -s)| ds, and ∞ j=-∞ ∞ -∞ |f (-s)| |f (j -s)| < ∞ by (2.4), it follows from Lebesgue's domi- nated convergence theorem that lim m→∞ v m = ∞ j=-∞ γ f (j). Hence by (2.2), Y (m) d → Y, m → ∞, where Y d = N   0, σ 2 1 0   ∞ j=-∞ f (u + j)   2 du   . (2.6) A similar argument gives lim m→∞ ∞ j=-∞ γ f -fm (j) = 0, so that lim m→∞ lim n→∞ Var n 1/2 (X n;1 -X (m) n;1 ) = lim m→∞ lim n→∞ n Var n -1 n t=1 ∞ -∞ (f (t -s) -f m (t -s)) dL s , = lim m→∞ ∞ j=-∞ γ f -fm (j) = 0,
where we used Theorem 7.1.1 in Brockwell and Davis [START_REF] Brockwell | Time Series: Theory and Methods[END_REF] for the second equality. An application of Chebychef's inequality then shows that lim m→∞ lim sup

n→∞ P (n 1/2 |X n;1 -X (m) n;1 | > ε) = 0
for every ε > 0. Together with (2.5) and (2.6) this implies the claim by a variant of Slutsky's theorem (cf. [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], Proposition 6.3.9).

Remark 2.3. Let us start with an easy remark on a necessary condition on the kernel f to apply the previous theorem. Obviously

F ∆ ∈ L 1 ([0, ∆]) is equivalent to f ∈ L 1 (R). Hence F ∆ ∈ L 2 ([0, ∆]) ⇒ f ∈ L 1 (R).
Remark 2.4. Unlike for the discrete time moving average process of (1.2), where absolute summability of the autocovariance function is guaranteed by absolute summability of the coefficient sequence, for the continuous time series model (1.1) it is not enough to assume that the kernel satisfies f ∈ L 1 (R) ∩ L 2 (R). An example is given by taking ∆ = 1 and

f (u) :=          0, u ≤ 0, 1, u ∈ [0, 1), 1•3•••(2j-1) 2 j j! (u -j) j , u ∈ [j, j + 1), j ∈ N.
For then the function F 1 is given by

F 1 (u) = j∈Z f (u + j) = (1 -u) -1/2 , u ∈ [0, 1), so that F 1 ∈ L 1 ([0, 1]) \ L 2 ([0, 1]). But F 1 ∈ L 1 ([0, 1]) is equivalent to f ∈ L 1 (R), and since |f (u)| ≤ 1 for all u ∈ R, this implies also f ∈ L 2 (R).
Observe further that for non-negative f , condition (2.1) is indeed necessary and sufficient for absolute summability of the autocovariance function.

Asymptotic normality of the sample autocovariance

As usual, we consider the stationary process

X t = ∞ -∞ f (t -s) dL s , t ∈ R. (3.1)
We recall that

γ * n;∆ (h∆) = n -1 n t=1 X t∆ X (t+h)∆ , h ∈ N
and first we establish an asymptotic result for Cov (γ * n;∆ (p∆), γ * n;∆ (q∆)).

Proposition 3.1. Let L be a (non-zero) Lévy process, with expectation zero, and finite fourth moment, and denote σ 2 := EL 2 1 and η := σ -4 EL 4 1 . Let ∆ > 0, and suppose further that f ∈ L 2 (R) ∩ L 4 (R) and that

[0, ∆] → R, u → ∞ k=-∞ f (u + k∆) 2 ∈ L 2 ([0, ∆]). (3.2)
For q ∈ Z denote

g q;∆ : [0, ∆] → R, u → ∞ k=-∞ f (u + k∆)f (u + (k + q)∆),
which belongs to L 2 ([0, ∆]), by the previous assumption. If further

∞ h=-∞ |γ(h∆)| 2 < ∞, (3.3)
then we have for each p, q ∈ N

lim n→∞ nCov (γ * n;∆ (p∆), γ * n;∆ (q∆)) = (η -3)σ 4 ∆ 0 g p;∆ (u)g q;∆ (u) du+ ∞ k=-∞ γ(k∆)γ((k -p + q)∆) + γ((k + q)∆)γ((k -p)∆) . (3.4)
Proof. For simplicity in notation we assume that ∆ = 1. The general case can be proved analogously or reduced to the case ∆ = 1 by a simple time change. We shall first show that for t, p, h, q

∈ Z E(X t X t+p X t+h+q X t+h+p+q ) = (η -3)σ 4 ∞ -∞ f (u)f (u + p)f (u + h + p)f (u + h + p + q) du +γ(p)γ(q) + γ(h + p)γ(h + q) + γ(h + p + q)γ(h). (3.5)
To show this, assume first that f is of the form

f (s) = f m,ǫ (s) = m/ǫ i=-m/ǫ ψ i 1 (iǫ,(i+1)ǫ] (s), (3.6)
where m ∈ N, ǫ > 0 such that 1/ǫ ∈ N, and

ψ i ∈ R, i = -m/ǫ, . . . , m/ǫ. Denote X t;m,ǫ := ∞ -∞ f m,ǫ (t -s) dL s = m/ǫ i=-m/ǫ ψ i (L t-iǫ -L t-(i+1)ǫ ), t ∈ R.

Denote further

Z i := L iǫ -L (i-1)ǫ , i ∈ Z.
Then (Z i ) i∈Z is i.i.d. and we have

X tǫ;m,ǫ = m/ǫ i=-m/ǫ ψ i Z t-i , t ∈ Z.
At this point we will need to compute the fourth moment of integrals of the Lévy process. Let us state an elementary result that yields a formula for this moment.

Lemma 3.2. Let φ ∈ L 2 (R) ∩ L 4 (R)
, then, with the assumptions and notations on L used in Proposition 3.1,

E( R φ(s)dL s ) 4 = (η -3)σ 4 R φ 4 (s)ds + 3σ 4 ( R φ 2 (s)ds) 2 . ( 3.7)
Proof. If ν is the Lévy measure of L and A its Gaussian variance, then by the Lévy Khintchine formula we get

ξ(u) = E exp iu R φ(s)dL s = exp - 1 2 Au 2 R φ 2 (s)ds + R×R [e iuφ(s)x -1 -iuφ(s)x]ν(dx)ds .
Then E( R φ(s)dL s ) 4 is obtained as the fourth derivative of ξ at u = 0. If we recall that (η -3)σ 4 = R x 4 ν(dx), and σ 2 = A+ R x 2 ν(dx), we get (3.7), after elementary but tedious computations.

To continue with the proof of Proposition 3.1, we now apply (3.7) to the special case where f (s) = 1 (0,ǫ] (s) and we get

EZ 2 i = EL 2 ǫ = σ 2 ǫ, EZ 4 i = EL 4 ǫ = ησ 4 ǫ -3σ 4 ǫ + 3σ 4 ǫ 2 . (3.8)
As shown in the proof of Proposition 7.3.1 in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], we then have

E(X t;m,ǫ X t+p;m,ǫ X t+h+p;m,ǫ X t+h+p+q;m,ǫ ) = EZ 4 i -3(EZ 2 i ) 2 m/ǫ i=-m/ǫ ψ i ψ i+p/ǫ ψ i+h/ǫ+p/ǫ ψ i+h/ǫ+p/ǫ+q/ǫ +γ m,ǫ (p)γ m,ǫ (q) + γ m,ǫ (h + p)γ m,ǫ (h + q) + γ m,ǫ (h + p + q)γ m,ǫ (h),
where γ m,ǫ (u) = E(X 0;m,ǫ X u;m,ǫ ), u ∈ R. By (3.8),

EZ 4 i -3(EZ 2 i ) 2 = (η -3)σ 4 ǫ, and
ǫ m/ǫ i=-m/ǫ ψ i ψ i+p/ǫ ψ i+h/ǫ+p/ǫ ψ i+h/ǫ+p/ǫ+q/ǫ = ∞ -∞ f (u)f (u + p)f (u + h + p)f (u + h + p + q) du, so that (3.5) follows for f of the form f = f m,ǫ . Now let f ∈ L 2 (R) ∩ L 4 (R)
f m k ,ǫ k converge to f both in L 2 (R)
and L 4 (R) and observing that γ m k ,ǫ k (u) → γ(u) for each u ∈ R. From (3.5) we conclude that, with p, q ∈ N,

Cov (γ * n;1 (p), γ * n;1 (q)) = n -1 |k|<n (1 -n -1 |k|)T k , (3.9)
where

T k = γ(k)γ(k -p + q) + γ(k + q)γ(k -p) +(η -3)σ 4 ∞ -∞ f (u)f (u + p)f (u + k)f (u + q + k) du. Now by (3.3), ∞ k=-∞ |T k | < ∞ if ∞ k=-∞ ∞ -∞ f (u)f (u + p)f (u + k)f (u + q + k) du (3.10) is finite. Denote G r (u) := ∞ k=-∞ |f (u + k)f (u + k + r)|, u ∈ R, r ∈ N.
Then G r is periodic, and by assumption, G r restricted to [0, 1] is square integrable.

Hence we can estimate (3.10) by

∞ k=-∞ ∞ -∞ |f (u)f (u + p)| |f (u + k)f (u + q + k)| du = ∞ h=-∞ h+1 h |f (u)f (u + p)| G q (u) du = ∞ h=-∞ 1 0 |f (u + h)f (u + p + h)| G q (u) du = 1 0 G p (u)G q (u) du < ∞.
The same calculation without the modulus and an application of the dominated convergence theorem to (3.9) then shows (3.4).

Remark 3.3. A sufficient condition for (3.3) is that ∞ h=-∞ |γ(h∆)| < ∞, which is implied by the function F ∆ in Theorem 2.1 belonging to L 2 ([0, ∆]). Another suf- ficient condition is that Φ : u → k∈Z |F(f )(u + 2πk/∆)| 2 , u ∈ [0, 2π/∆] is in L ∞ ([0, 2π/∆]), where F(f ) is the Fourier transform of f ∈ L 2 (R) in the form z → ∞ -∞ e izt f (t) dt (for L 1 -functions). For if Φ ∞ ≤ B, then (f (• + h∆)) h∈Z is a Bessel sequence in L 2 (R) with bound B/∆, i.e. ∞ h=-∞ ∞ -∞ ϕ(u)f (u + h∆) du 2 ≤ B∆ -1 ∞ -∞ ϕ(u) 2 du ∀ ϕ ∈ L 2 (R),
see e.g. Theorem 7.2.3 in Christensen [START_REF] Christensen | An Introduction to Frames and Riesz Bases[END_REF]. Taking ϕ = f then gives the square summability of the autocovariance functions by (1.12).

Please remark that ∞ h=-∞ γ(h∆) 2 < ∞ cannot be deduced from the condition that u → ∞ k=-∞ f (u + k∆) 2 is in L 2 ([0, ∆]). One can take ∆ = 1 and f (s) = i≥1

1 (i,i+1] (s) i H for 1 2 < H ≤ 3 4 , to get the latter condition but not ∞ h=-∞ γ(h) 2 < ∞.
Remark 3.4. By (1.12), the condition (3.3) can be written as

∞ k=-∞ ∞ -∞ f (s)f (s + k∆) ds

Then we have for each

h ∈ N √ n(γ * n;∆ (0) -γ(0), . . . , γ * n;∆ (h∆) -γ(h)) ′ d → N (0, V ), n → ∞, (3.12)
where V = (v pq ) p,q=0,...,h ∈ R h+1,h+1 is the covariance matrix defined by

v pq = (η -3)σ 4 ∆ 0 g p;∆ (u)g q;∆ (u) du+ ∞ k=-∞ γ(k∆)γ((k -p + q)∆) + γ((k + q)∆)γ((k -p)∆) . (3.13)
(b) In addition to the assumptions of (a), assume that the function

u → ∞ j=-∞ |f (u + j∆)| is in L 2 ([0, ∆]
). Denote by

γ n;∆ (j∆) = n -1 n-j t=1 (X t∆ -X n;∆ )(X (t+j)∆ -X n;∆ ), j = 0, 1, . . . , n -1,
the sample autocovariance, as defined in (1.4). Then we have for each

h ∈ N √ n( γ n;∆ (0) -γ(0), . . . , γ n;∆ (h∆) -γ(h)) ′ d → N (0, V ), n → ∞,
where V = (V pq ) p,q=0,...,h is defined by (3.13).

(c) For j ∈ N let ρ * n;∆ (j∆) = γ * n;∆ (j∆)/γ * n;∆ (0) and ρ n (j∆) = γ n;∆ (j∆)/ γ n;∆ (0), the latter being the sample autocorrelation at lag j∆. Suppose that f is not almost everywhere equal to zero. Then, under the assumptions of (a), we have for each

h ∈ N, that √ n(ρ * n;∆ (∆) -ρ(∆), . . . , ρ * n;∆ (h∆) -ρ(h∆)) ′ d → N (0, W ), n → ∞, (3.14)
where W = W ∆ = (w ij;∆ ) i,j=1,...,h is given by

w ij;∆ = w ij;∆ + (η -3)σ 4 γ(0) 2 ∆ 0 g i;∆ (u) -ρ(i∆)g 0;∆ (u) (g j;∆ (u) -ρ(j∆)g 0;∆ (u) du, and
w ij;∆ = ∞ k=-∞ ρ((k + i)∆)ρ((k + j)∆) + ρ((k -i)∆)ρ((k + j)∆) + 2ρ(i∆)ρ(j∆)ρ(k∆) 2 -2ρ(i∆)ρ(k∆)ρ((k + j)∆) -2ρ(j∆)ρ(k∆)ρ((k + i)∆) = ∞ k=1 ρ((k + i)∆) + ρ((k -i)∆) -2ρ(i∆)ρ(k∆) × ρ((k + j)∆) + ρ((k -j)∆) -2ρ(j∆)ρ(k∆)
is given by Bartlett's formula. If additionally the function

u → ∞ j=-∞ |f (u + j∆)| is in L 2 ([0, ∆]), then it also holds that √ n( ρ n;∆ (∆) -ρ(∆), . . . , ρ n;∆ (h∆) -ρ(h∆)) ′ d → N (0, W ), n → ∞. (3.15)
Proof. For simplicity in notation we assume again ∆ = 1 in this proof.

(a) Using Proposition 3.1 it follows as in the proof of Proposition 7.3.2 in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], that the claim is true if f has additionally compact support. For general f and m ∈ N let f m := f 1 (-m,m) . Hence we have that

n 1/2 (γ * n;(m) (0) -γ m (0), . . . , γ * n;(m) (h) -γ m (h)) ′ d → Y m , m → ∞,
where γ m is the autocovariance function of the process

X t;m = ∞ -∞ f m (t -s) dL s , γ * n;(m) (p) = n -1 n t=1 X t;
g p;(m) (u)g q;(m) (u) du+ ∞ k=-∞ γ m (k)γ m (k-p+q)+γ m (k+q)γ m (k-p) .
Here,

g p;(m) (u) = ∞ k=-∞ f m (u + k)f m (u + k + p), u ∈ [0, 1]. Next, we want to show that lim m→∞ V m = V . Observe first that g p;(m) (u) = ∞ k=-∞ f m (u+k)f m (u+k+p) → ∞ k=-∞ f (u+k)f (u+k+p) = g p,1 (u) =: g p (u) almost surely in the variable u as m → ∞ by Lebesgue's dominated convergence theorem, since u → ∞ k=-∞ |f (u + k)f (u + k + p)| is in L 2 ([0, 1]) by (3.
2) and hence is almost surely finite. Further we have

|g p;(m) (u)| ≤ ∞ k=-∞ |f (u + k)f (u + k + p)|
uniformly in u and m, so that again by the dominated convergence theorem we have that g

p;(m) → g p in L 2 ([0, 1]) as m → ∞. Next, observe that |γ m (k)| ≤ ∞ -∞ |f (s)f (s + k)| ds ∀ m ∈ N ∀ k ∈ Z.
Since lim m→∞ γ m (k) = γ(k) for every k ∈ Z, it follows from the dominated convergence theorem and (3.11) that (γ m (k)) k∈Z converges in l 2 (Z) to (γ(k)) k∈Z . This together with the convergence of g p;(m) gives the desired lim m→∞ V m = V , so that

Y m d → Y, m → ∞, where Y d = N (0, V ). Finally, that lim m→∞ lim sup n→∞ P (n 1/2 |γ * n;(m) (p)-γ m (p)-γ * (p)+γ(p)| > ε) = 0 ∀ ε > 0, p ∈ {0, . . . , h}
follows as in Equation (7.3.9) in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]. An application of a variant of Slutsky's theorem (cf. [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], Proposition 6.3.9) then gives the claim.

(b) This follows as in the proof of Proposition 7.3.4 in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]. One only has to observe that by Theorem 2.1, √ n X n;1 converges in distribution to a normal random variable as n → ∞. In particular, X n;1 must converge to 0 in probability as n → ∞.

(c) The limit theorem follows as in the proof of Theorem 7.2.1 in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF], and for w ij we have the representation

w ij;∆ = (v ij -ρ(i)v 0j -ρ(j)v i0 + ρ(i)ρ(j)v 00 )/γ(0) 2 = w ij;∆ + (η -3)σ 4 γ(0) 2 × 1 0 g i (u)g j (u) -ρ(i)g 0 (u)g j (u) -ρ(j)g i (u)g 0 (u) + ρ(i)ρ(j)g 0 (u) 2 du,
giving the claim.

Remark 3.6. It is easy to check that

w ij;∆ = w ij;∆ if f is of the form f = ∞ i=-∞ ψ i 1 (i∆,(i+1)∆]
, in accordance with Bartlett's formula, since then (X t∆ ) t∈Z has a discrete time moving average representation with i.i.d. coefficients. Remark 3.7. Another case when w ij;∆ = w ij;∆ is when η = 3, happens if and only if L is Brownian motion. However, in general we do not have w ij;∆ = w ij;∆ . An example is given by f = 1 (0,1/2] + 1 [START_REF] Ash | Topics in Stochastic Processes[END_REF][START_REF] Benassi | On roughness indices for fractional fields[END_REF] and ∆ = 1, in which case g 1;1 = 1 (0,1/2] and g 0;

1 = 2 • 1 (0,1/2] + 1 (1/2,1]
, and it is easy to see that g 1;1 -ρ(1)g 0;1 is not almost everywhere zero, so that w 11;1 = w 11;1 if η = 3. The latter example corresponds to a moving average process, which is varying at the scale 1 2 , but sampled at integer times. Observe however that w 11;1/2 = w 11;1/2 by Remark 3.6. A more detailed study of such phenomena in discrete time can be found in Niebuhr and Kreiss [START_REF] Niebuhr | Asymptotics for autocovariances and integrated periodograms for linear processes observed at lower frequencies[END_REF].

Remark 3.8. Recently, sophisticated and powerful results on the normal approximation of Poisson functionals using Malliavan calculus have been obtained. E.g., Peccati and Taqqu [15, Theorems 2, 3 and 5] prove a central limit theorem for double Poisson integrals and apply this to a specific quadratic functional of a Lévy driven Ornstein-Uhlenbeck process, and Peccati et al. [START_REF] Peccati | Stein's method and normal approximation of Poisson functionals[END_REF]Section 4] obtain bounds for such limit theorems, to name just of few of some recent publications on this subject. It may be possible to apply the results of [START_REF] Peccati | Stein's method and normal approximation of Poisson functionals[END_REF][START_REF] Peccati | Central limit theorems for double Poisson integrals[END_REF] to obtain another proof of Theorem 3.5 under certain conditions such as finite 6th moment, but we have not investigated this issue further. Note that our proof uses only basic knowledge of stochastic integrals and methods from time series analysis.

An application to fractional Lévy noise

We will now apply the previous results to fractional Lévy processes. Recall from (1.8) and (1.9) that these were denoted by

M 1 t;d := 1 Γ(d + 1) ∞ -∞ (t -s) d + -(-s) d + dL s , t ∈ R, and
M 2 t;d = ∞ -∞ |t -s| d -|s| d dL s , t ∈ R,
respectively, and the corresponding fractional Lévy noises based on increments of length ∆ > 0 by

X i t = M i t;d -M i t-∆;d , t ∈ R, i = 1, 2.
Hence the fractional Lévy noises are Lévy driven moving average processes with kernel functions

f 1 d,∆ (s) = 1 Γ(d + 1) s d + -(s -∆) d + , s ∈ R, and
f 2 d,∆ (s) = d -|s -∆| d , s ∈ R, respectively. Neither f 1 d,∆ nor f 2 d,∆ are in L 1 (R)
, so Theorem 2.1 cannot be applied because of Remark 2.3. Please note that for the same reason the assumptions for (b) of Theorem 3.5 and for (3.15) are not fulfilled.

For simplicity in notation we assume ∆ = 1, and drop the subindex ∆. Although the fractional noises X i have different distributions for i = 1 and i = 2, they are both stationary with the autocovariance

E(X i t+h X i t ) = γ X i (h) = C i (d)σ 2 2 |h + 1| 2d+1 -2|h| 2d+1 + |h -1| 2d+1 , (4.1)
where C i (d) is a normalising multiplicative constant depending on d. Both processes X i are infinitely divisible and of moving average type, hence we know from [START_REF] Cambanis | Chaotic behavior of infinitely divisible processes[END_REF][START_REF] Fuchs | Mixing conditions for multivariate infinitely divisible processes with an application to mixed moving averages and the sup ou stochastic volatility model[END_REF] that (X i t ) t∈Z is mixing in the ergodic-theoretic sense. For fixed h ∈ Z, define the function

F : R Z → R, (x n ) n∈Z → x 0 x h .
If T denotes the forward shift operator, then

F (T k (X i t ) t∈Z ) = X i k X i k+h ,
and from Birkhoff's ergodic theorem (e.g. Ash and Gardner [START_REF] Ash | Topics in Stochastic Processes[END_REF], Theorems 3.3.6 and 3.3.10) we know that

1 n n k=1 X i k X i k+h → E F ((X i t ) t∈Z ) = EX i 0 X i h , n → ∞, for i = 1, 2,
(h) = C i (d)σ 2 2 |h + 1| 2d+1 -2|h| 2d+1 + |h -1| 2d+1 a.s. Since γ(0) = C i (d)σ 2 and γ(1) = C i (d)σ 2 (2 2d -1), ρ * n (1) = γ * n (1) γ *
n (0) is a strongly consistent estimator for 2 2d -1. Hence,

d := 1 2 log(ρ * n (1) + 1) log 2 (4.2)
is a strongly consistent estimator for d.

The question of the asymptotic normality of these estimators arises naturally. There are many classical techniques to show the asymptotic normality of an ergodic stationary sequence by assuming some stronger mixing assumption. As far as we know, they do not work in our setting. To illustrate this point, we would like to show that fractional Lévy noises are not strongly mixing. Let us first recall the definition.

Definition 4.1. Let (X n ) n∈Z be a stationary sequence, and let

α X (n) = sup{|P (A ∩ B) -P (A)P (B)|, A ∈ σ(X k , k ≤ 0), B ∈ σ(X k , k ≥ n)}.
The sequence (X n ) n∈Z is strongly mixing if lim n→∞ α X (n) = 0.

In our case we know the weak mixing property lim n→∞ |P (A∩B)-P (A)P (B)| = 0 for A ∈ σ(X 0 ), B ∈ σ(X n ), because of [START_REF] Cambanis | Chaotic behavior of infinitely divisible processes[END_REF][START_REF] Fuchs | Mixing conditions for multivariate infinitely divisible processes with an application to mixed moving averages and the sup ou stochastic volatility model[END_REF]. There are classical central limit theorems for strongly mixing sequences, see e.g. [START_REF] Merlevéde | Recent advances in invariance principles for stationary sequences[END_REF] for an overview. The following result, which is stated as Proposition 34 in [START_REF] Merlevéde | Recent advances in invariance principles for stationary sequences[END_REF], will be particularly useful for us. Theorem 4.2. Suppose that (X t ) t∈Z is a mean zero, strongly mixing sequence and that there exists some δ > 0 and a constant K > 0 such that

E|X 0 | 2+δ < ∞, (4.3) lim m→∞ Var ( m i=1 X i ) = ∞, (4.4)
E| m i=1 X i | 2+δ ≤ K Var ( m i=1 X i ) 1+δ/2 ∀ m ∈ N. (4.5) Write S (n) (t) := ⌊nt⌋ i=1 X i , n ∈ N, t ∈ [0, 1],
where ⌊x⌋ is the integer part of the real number x. Then

Var ( n i=1 X i ) -1/2 S (n) d → B weakly in D[0, 1], (4.6)
where B is a standard Brownian motion.

Lemma 4.3. Let L be a two sided non-zero Lévy process with expectation zero and finite fourth moment, and let d ∈ (0, 1/2). Let

X i t := M i t;d -M i t-1;d , t ∈ Z,
be the corresponding fractional Lévy noises, for i = 1, 2. Then (X i t ) t∈Z satisfies (4.3) -(4.5) with δ = 2, for i = 1, 2.

Proof. The proof is only written for the fractional noise X 1 denoted by X but it is similar for X 2 . Equation (4.3) holds since L has finite fourth moment and since the kernel function of fractional noise is in L 2 (R) ∩ L 4 (R), and (4.4) follows from the fact that m i=1

X i = M m;d , (4.7)
and

Var (M m;d ) = Cm 2d+1 ∀ m ∈ N (4.8)
for some constant C. To see (4.5) for δ = 2, we use (4.7) and

f m (s) := 1 Γ(d + 1) [(m -s) d + -(-s) d + ], s ∈ R.
Then by Lemma 3.2,

E|M d (m)| 4 = (η -3)σ 4 R f 4 m (s)ds + 3σ 4 R f 2 m (s)ds 2 . Observe that R f 2 m (s) ds 2 = Cm 4d+2 and that R f 4 m (s) ds ≤ C ′ m 4d+1
for positive constants C, C ′ , which gives the claim.

A consequence of the previous theorem and lemma is the following negative result.

Corollary 4.4. Assume the assumptions of Lemma 4.3. Then the fractional Lévy noises X 1 and X 2 are not strongly mixing.

Proof. If fractional Lévy noises were strongly mixing, then (4.6) would follow. Please remark that, since fractional noises are increments, S (n) (t) = M i d (⌊nt⌋) -M i d (0), and (Var ( n i=1 X i )) = Cn 2d+1 by (4.8). Owing to the asymptotic self-similarity of fractional Lévy processes (Proposition 3.1 in [START_REF] Benassi | On roughness indices for fractional fields[END_REF]), we know that

M 2 d (⌊nt⌋) -M 2 d (0) n d+1/2 d → B d+1/2 (t),
where the limit is a fractional Brownian motion with Hurst exponent d + 1/2. A similar asymptotic self-similarity holds for M 1 d . Hence (4.6) is violated and the corollary is proved by contradiction.

Nevertheless one can apply Theorem 3.5 to get the asymptotic normality of the estimator (4.2). Proof. The proof is only written for the fractional noise X 1 denoted by X but it is similar for X 2 . We shall apply (3.14) to get convergence of √ n(ρ * n (1) -ρ(0)) to a Gaussian random variable. First, observe that

f 1 d,1 ∈ L 2 (R) ∩ L 4 (R), (4.9)
and ∞ k=-∞ γ(k) 2 < ∞. The latter inequality is classical for the fractional Gaussian noise, when d < 1/4, and holds for fractional Lévy noises since they have the same autocorrelation as fractional Gaussian noise. This also implies (3.11)

Z i t = X i t -X i t-1 , t ∈ Z,
for which it holds ∞ k=-∞ γ Z i (k) 2 < ∞. We conclude from Birkhoff's ergodic theorem that γ * n,Z (h) := Proof. The proof is only written for the fractional noise X 1 denoted by X but it is similar for X 2 . Let us remark that

1 n n k=2 Z k Z k+h → E(Z 0 Z h ), n → ∞,
Z t = ∞ -∞ f 1 d,1 (t -s) dL s , where f 1 d,1 (s) = f 1 d,1 (s) -f 1 d,1 (s -1),
To apply Theorem 3.5, we have to check that

f 1 d,1 (s) ∈ L 2 (R) ∩ L 4 (R),
which is obvious from (4.9). Moreover we already know that ∞ k=-∞ γ Z i (k) 2 < ∞. This time, however, the kernel function f 1 d,1 is not nonnegative, but it is easy to see that | f

2 <

 2 ∞. The assumption (3.11) used in Theorem 3.5 below is slightly stronger than (3.3), but equivalent to (3.3) if f ≥ 0. The following theorem gives asymptotic normality of the sample autocovariance and sample autocorrelation and the related estimators γ * n;∆ and ρ * n;∆ . Theorem 3.5. (a) Suppose the assumptions of Proposition 3.1 are satisfied and suppose further that∞

  m X t+p;m the corresponding autocovariance estimate, and Y m d = N (0, V m ) with V m = (v pq;m ) p,q=0,...,h and v pq;m = (η-3)σ 4 1 0

Proposition 4 . 5 .

 45 Assume the assumptions of Lemma 4.3, and let d be defined by (4.2). If d ∈ (0, 1/4) then √ n( d -d) converges in distribution to a Gaussian random variable as n → ∞.

  by Remark 3.4 since f 1 d,1 ≥ 0. Let us check that g 0 := g 0;1 ∈ L 2 (0, 1). SinceΓ 2 (d + 1)g 0 (u) = ∞ k=-∞ ((u + k) d + -(u + k -1) d + ) 2 ,it follows for all u ∈ (0, 1) thatΓ 2 (d + 1)g 0 (u) ≤ g 0 ∈ L ∞ ([0, 1]).Hence the assumptions of Theorem 3.5 (a) are fulfilled, and the result follows from(3.14).If d ≥ 1/4 then ∞ k=-∞ γ(k) 2 = ∞,since fractional Lévy noises have the same autocorrelation as the fractional Gaussian noise. Hence we consider

7 8-2 2d+2 .. 10 ) 4 . 6 .

 72d+21046 is a strongly consistent estimator ofE(Z 0 Z h ) = C i (d)σ 2 2 -|h + 2| 2d+1 + 4|h + 1| 2d+1 -6|h| 2d+1 + 4|h -1| 2d+1 -|h -2| 2d+1 . Therefore ρ * n,Z (1) = γ * n,Z (1) γ * n,Z(0) is a strongly consistent estimator for φ(d) = -3 2d+1 +4•2 2d+1 -It turns out that φ is increasing on (0, 1/2). Therefore we can define the estimatord := φ -1 (ρ * n,Z (1)). (4PropositionAssume the assumptions of Lemma 4.3, and let d be defined by (4.10). If d ∈ (0, 1/2) then √ n( d-d) converges in distribution to a Gaussian random variable as n → ∞.

1 d, 1

 1 (t)| ≤ C min(1, |t| d-2 ) and hence that ∞ -∞ | f 1 d,1 (t) f 1 d,1 (t + k)| dt ≤ C ′ min(1, |k| d-1 ), ∀ k ∈ Z,for some constants C, C ′ , giving(3.11). Finally,Γ 2 (d + 1)g 0 (u) = ∞ k=-∞ ((u + k) d + + (u + k -2) d + -2(u + k -1) d + ) 2 ,and estimating the summands separately for k < 0, k = 0, 1 and k ≥ 2 we obtain for u ∈ [0, 1]Γ 2 (d + 1)g 0 (u) ≤ 1 + (2 d + 2) 2 + ∞ k=2 (k d -(k -1) d ) 2 < ∞, so that g 0 ∈ L ∞ ([0, 1]) ⊂ L 2 ([0,1]). The claim now follows from Theorem 3.5, using(3.14).

  and X t , t ∈ R, defined by(3.1). Then there is a sequence of functions (f m k ,ǫ k ) k∈N of the form (3.6) such that f m k ,ǫ k converges to f both in L 2 (R) and in L 4 (R) as k → ∞. Then for each fixed t ∈ R, we have that X t;m k ,ǫ k → X t in L 2 (P ) (P the underlying probability measure) as k → ∞, where we used the Itô isometry. Further, by Lemma 3.2, and convergence of f m k ,ǫ k both in L 2 (R) and in L 4 (R), we get convergence of X t;m k ,ǫ k to X t in L 4 (P ). This then shows(3.5), by letting

Acknowledgement

Major parts of this research where carried out while AL was visiting the Institut de Mathématiques de Toulouse as guest professor in 2008 and 2009. He takes pleasure in thanking them for the kind hospitality and their generous financial support.


Keywords: noyaux à valeurs opérateurs, apprentissage de noyaux multiples, analyse des données fonctionnelles, espace de Hilbert à noyau reproduisant

Apprentissage de Noyaux à Valeurs Opérateurs Multiples Résumé : Dans cet article, nous proposons une méthode d'apprentissage de noyaux multiples à valeurs opérateurs dans le cas d'une régression ridge à réponse fonctionnelle. Notre méthode est basée sur la résolution d'un système d'équations linéaires d'opérateurs en utilisant une procédure de type Iterative Coordinate Descent. Nous validons expérimentalement notre approche sur un problème de prédiction de mouvement de doigt dans un contexte d'Interface Cerveau-Machine.

Introduction

During the past decades, a large number of algorithms have been proposed to deal with learning problems in the case of single-valued functions (e.g., binaryoutput function for classification or real output for regression). Recently, there has been considerable interest in estimating vector-valued functions [START_REF] Micchelli | On learning vector-valued functions[END_REF][START_REF] Caponnetto | Universal multitask kernels[END_REF][START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF]. Much of this interest has arisen from the need to learn tasks where the target is a complex entity, not a scalar variable. Typical learning situations include multi-task learning [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF], functional regression [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], and structured output prediction [START_REF] Brouard | Semi-supervised penalized output kernel regression for link prediction[END_REF].

In this paper, we are interested in the problem of functional regression with functional responses in the context of brain-computer interface (BCI) design. More precisely, we are interested in finger movement prediction from electrocorticographic signals [START_REF] Schalk | Decoding two-dimensional movement trajectories using electrocorticographic signals in humans[END_REF]. Indeed, from a set of signals measuring brain surface electrical activity on d channels during a given period of time, we want to predict, for any instant of that period whether a finger is moving or not and the amplitude of the finger flexion. Formally, the problem consists in learning a functional dependency between a set of d signals and a sequence of labels (a step function indicating whether a finger is moving or not) and between the same set of signals and vector of real values (the amplitude function). While, it is clear that this problem can be formalized as functional regression problem, from our point of view, such problem can benefit from the multiple operator-valued kernel learning framework. Indeed, for these problems, one of the difficulties arises from the unknown latency between the signal related to the finger movement and the actual movement [START_REF] Pistohl | Prediction of arm movement trajectories from ecog-recordings in humans[END_REF]. Hence, instead of fixing in advance some value for this latency in the regression model, our framework allows to learn it from the data by means of several operator-valued kernels.

If we wish to address functional regression problem in the principled framework of reproducing kernel Hilbert spaces (RKHS), we have to consider RKHSs whose elements are operators that map a function to another function space, possibly source and target function spaces being different. Working in such RKHSs, we are able to draw on the important core of work that has been performed on scalar-valued and vector-valued RKHSs [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF][START_REF] Micchelli | On learning vector-valued functions[END_REF]. Such a functional RKHS framework and associated operator-valued kernels have been introduced recently [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF][START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF]. A basic question with reproducing kernels is how to build these kernels and what is the optimal kernel choice for a given application. In order to overcome the need for choosing a kernel before the learning process, several works have tried to address the problem of learning the scalar-valued kernel jointly with the decision function [START_REF] Lanckriet | Learning the kernel matrix with semi-definite programming[END_REF][START_REF] Sonnenburg | Large scale multiple kernel learning[END_REF]. Since these seminal works, many efforts have been carried out in order to theoretically analyze the kernel learning framework [START_REF] Cortes | Generalization bounds for learning kernels[END_REF][START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF] or in order to provide efficient algorithms [23,[START_REF] Aflalo | Variable sparsity kernel learning[END_REF][START_REF] Kloft | ℓ p -norm multiple kernel learning[END_REF]. While many works have been devoted to multiple scalar-valued kernel learning, this problem of kernel learning have been barely investigated for operator-valued kernels. One motivation of this work is to bridge the gap between multiple kernel learning (MKL) and operator-valued kernels by proposing a framework and an algorithm for learning a finite linear combination of operator-valued kernels. While each step of the scalar-valued MKL framework can be extended without major difficulties to operator-valued kernels, technical challenges arise at all stages because we deal with infinite dimensional spaces. It should be pointed out that in a recent work [START_REF] Dinuzzo | Learning output kernels with block coordinate descent[END_REF], the problem of learning the output kernel was formulated as an optimization problem over the cone of positive semidefinite matrices, and proposed a block-coordinate descent method to solve it. However, they did not focus on learning the input kernel. In contrast, our multiple operator-valued kernel learning formulation can be seen as a way of learning simultaneously input and output kernels, although we consider a linear combination of kernels which are fixed in advance.

In this paper, we make the following contributions:

• we introduce a novel approach to infinite-dimensional multiple operatorvalued kernel learning (MovKL) suitable for learning the functional dependencies and interactions between continuous data,

• we extend the original formulation of ridge regression in dual variables to the functional data analysis domain, showing how to perform nonlinear functional regression with functional responses by constructing a linear regression operator in an operator-valued kernel feature space (Section 2),

• we derive a dual form of the MovKL problem with functional ridge regression, and show that a solution of the related optimization problem exists (Section 2),

• we propose a block-coordinate descent algorithm to solve the MovKL optimization problem which involves solving a challenging linear system with a sum of block operator matrices, and show its convergence in the case of compact operator-valued kernels (Section 3),

• we provide an empirical evaluation of MovKL performance which demonstrates its effectiveness on a BCI dataset (Section 4).

Problem Setting

Before describing the multiple operator-valued kernel learning algorithm that we will study and experiment with in this paper, we first review notions and properties of reproducing kernel Hilbert spaces with operator-valued kernels, show their connection to learning from multiple response data (multiple outputs; see [START_REF] Micchelli | On learning vector-valued functions[END_REF] for discrete data and [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF] for continuous data), and describe the optimization problem for learning kernels with functional response ridge regression.

Notations and Preliminaries

We start by some standard notations and definitions used all along the paper. Given a Hilbert space H, •, • H and • H refer to its inner product and norm, respectively. We denote by G x and G y the separable real Hilbert spaces of input and output functional data, respectively. In functional data analysis domain, continuous data are generally assumed to belong to the space of square integrable functions L 2 . In this work, we consider that G x and G y are the Hilbert space L 2 (Ω) which consists of all equivalence classes of square integrable functions on a finite set Ω. Ω being potentially different for G x and G y . We denote by F (G x , G y ) the vector space of functions f : G x -→ G y , and by L(G y ) the set of bounded linear operators from G y to G y . We consider the problem of estimating a function f such that f (x i ) = y i when observed functional data (x i , y i ) i=1,...,n ∈ (G x , G y ). Since G x and G y are INRIA spaces of functions, the problem can be thought of as an operator estimation problem, where the desired operator maps a Hilbert space of factors to a Hilbert space of targets. We can define the regularized operator estimate of f ∈ F as:

f λ arg min f ∈F 1 n n i=1 y i -f (x i ) 2 Gy + λ f 2 F . (1)
In We now introduce (function) G y -valued reproducing kernel Hilbert spaces and show the correspondence between such spaces and positive definite (operator) L(G y )-valued kernels. These extend the traditional properties of scalarvalued kernels.

Definition 1 (function-valued RKHS)

A Hilbert space F of functions from G x to G y is called a reproducing kernel Hilbert space if there is a positive definite

L(G y )-valued kernel K F (w, z) on G x × G x such that: i. the function z -→ K F (w, z)g belongs to F , ∀z ∈ G x , w ∈ G x , g ∈ G y , ii. ∀f ∈ F , w ∈ G x , g ∈ G y , f, K F (w, •)g F = f (w), g Gy (reproducing property). Definition 2 (operator-valued kernel) An L(G y )-valued kernel K F (w, z) on G x is a function K F (•, •) : G x × G x -→ L(G y ); furthermore: i. K F is Hermitian if K F (w, z) = K F (z, w) *
, where * denotes the adjoint operator,

ii. K F is positive definite on G x if it is Hermitian and for every natural number r and all {(w i , u i ) i=1,...,r } ∈ G x × G y , i,j K F (w i , w j )u j , u i Gy ≥ 0.

Theorem 1 (bijection between function valued RKHS and operator-valued kernel)

An L(G y )-valued kernel K F (w, z) on G x is the reproducing kernel of some Hilbert space F , if and only if it is positive definite.

The proof of Theorem 1 can be found in [START_REF] Micchelli | On learning vector-valued functions[END_REF]. For further reading on operatorvalued kernels and their associated RKHSs, see, e.g., [START_REF] Caponnetto | Universal multitask kernels[END_REF][START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and mercer theorem[END_REF][START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF].

Functional Response Ridge Regression in Dual Variables

We can write the ridge regression with functional responses optimization problem (1) as follows:

min f ∈F 1 2 f 2 F + 1 2nλ n i=1 ξ i 2 Gy with ξ i = y i -f (x i ).
(2)

RR n°7900

Now, we introduce the Lagrange multipliers α i , i = 1, . . . , n which are functional variables since the output space is the space of functions G y . For the optimization problem (2), the Lagrangian multipliers exist and the Lagrangian function is well defined. The method of Lagrange multipliers on Banach spaces, which is a generalization of the classical (finite-dimensional) Lagrange multipliers method suitable to solve certain infinite-dimensional constrained optimization problems, is applied here. For more details, see [START_REF] Kurcyusz | On the existence and nonexistence of lagrange multipliers in Banach spaces[END_REF]. Let α = (α i ) i=1,...,n ∈ G n y the vector of functions containing the Lagrange multipliers, the Lagrangian function is defined as

L(f, α, ξ) = 1 2 f 2 F + 1 2nλ ξ 2 G n y + α, y -f (x) -ξ G n y , (3)
where

α = (α 1 , . . . , α n ) ∈ G n y , y = (y 1 , . . . , y n ) ∈ G n y , ξ = (ξ 1 , . . . , ξ n ) ∈ G n y , f (x) = (f (x 1 ), . . . , f (x n )) ∈ G n y , and ∀a, b ∈ G n y , a, b G n y = n i=1 a i , b i Gy .
Differentiating (3) with respect to f ∈ F and setting to zero, we obtain

f (.) = n i=1 K(x i , .)α i , (4)
where

K : G x × G x -→ L(G y ) is the operator-valued kernel of F .
Substituting this into (3) and minimizing with respect to ξ, we obtain the dual of the functional response ridge regression problem

max α - nλ 2 α 2 G n y - 1 2 Kα, α G n y + α, y G n y , (5)
where K = [K(x i , x j )] n i,j=1 is the block operator kernel matrix.

MovKL in Dual Variables

Let us now consider that the function

f (•) is sum of M functions {f k (•)} M k=1
where each f k belongs to a G y -valued RKHS with kernel K k (•, •). Similarly to scalar-valued multiple kernel learning, we adopt the convention that x 0 = 0 if x = 0 and ∞ otherwise, and we can cast the problem of learning these functions

f k as min d∈D min f k ∈F k M k=1 f k 2 F k 2d k + 1 2λ n i=1 ξ i 2 Gy with ξ i = y i -M k=1 f k (x i ), (6)
where

d = [d 1 , • • • , d M ], D = {d : ∀k, d k ≥ 0 and k d r k ≤ 1} and 1 ≤ r ≤ ∞.
Note that this problem can equivalently be rewritten as an unconstrained optimization problem. Before deriving the dual of this problem, it can be shown by means of the generalized Weierstrass theorem [START_REF] Kurdila | Convex Functional Analysis[END_REF] that this problem admits a solution (a detailed proof is provided in the supplementary material). Now, following the lines of [23], a dualization of this problem leads to the following equivalent one

min d∈D max α∈G n y - nλ 2 α 2 G n y - 1 2 Kα, α G n y + α, y G n y , (7)

INRIA

where K = M k=1 d k K k and K k is the block operator kernel matrix associated to the operator-valued kernel K k . The KKT conditions also state that at optimality we have

f k (•) = n i=1 d k K k (x i , •)α i .

Solving the MovKL Problem

After having presented the framework, we now devise an algorithm for solving this MovKL problem.

Block-coordinate descent algorithm

Since the optimization problem ( 6) has the same structure as a multiple scalarvalued kernel learning problem, we can build our MovKL algorithm upon the MKL literature. Hence, we propose to borrow from [START_REF] Kloft | ℓ p -norm multiple kernel learning[END_REF], and consider a blockcoordinate descent method. The convergence of a block coordinate descent algorithm which is related closely to the Gauss-Seidel method was studied in works of [START_REF] Tseng | Convergence of block coordinate descent method for nondifferentiable minimization[END_REF] and others. The difference here is that we have operators and block operator matrices rather than matrices and block matrices, but this doesn't increase the complexity if the inverse of the operators are computable (typically analytically or by spectral decomposition). Our algorithm iteratively solves the problem with respect to α with d being fixed and then with respect to d with α being fixed (see Algorithm 1). After having initialized {d k } to non-zero values, this boils down to the following steps : 1. with {d k } fixed, the resulting optimization problem with respect to α has a simple closed-form solution:

(K + λI)α = 2y, (8)
where

K = M k=1 d k K k .
While the form of solution is rather simple, solving this linear system is still challenging in the operator setting and we propose below an algorithm for its resolution.

2. with {f k } fixed, according to problem (6), we can rewrite the problem as

min d∈D M k=1 f k 2 F k d k (9)
which has a closed-form solution and for which optimality occurs at [START_REF] Micchelli | Learning the kernel function via regularization[END_REF]:

d k = f k 2 r+1 ( k f k 2r r+1 ) 1/r . ( 10
)
This algorithm is similar to that of [START_REF] Cortes | L 2 regularization for learning kernels[END_REF] and [START_REF] Kloft | ℓ p -norm multiple kernel learning[END_REF] both being based on alternating optimization. The difference here is that we have to solve a linear system involving a block-operator kernel matrix which is a combination of basic kernel matrices associated to M operator-valued kernels. This makes the system very challenging, and we present an algorithm for solving it in the next paragraph.

Algorithm 1 ℓr-norm MovKL Input K k for k = 1, . . . , M d 1 k ←- 1 M for k = 1, . . . , M α ←-0 for t = 1, 2, . . . do α ′ ←-α K ←-k d t k K k α ←-solution of (K + λI)α = 2y if α -α ′ < ǫ then break end if d t+1 k ←- f k 2 r+1 ( k f k 2r r+1 ) 1/r for k = 1, . . . , M end for
A detailed proof of convergence of the MovKL algorithm, in the case of compact operator-valued kernels, is given in the supplementary material. It proceeds by showing that the sequence of functions {f (n) k } generated by the above alternating optimization produces a non-increasing sequence of objective values of Equation [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and mercer theorem[END_REF]. Then, using continuity and boundedness arguments, we can also show that the sequence {f (n) k } is bounded and thus converges to a minimizer of Equation ( 6). The proof is actually an extension of results obtained by [START_REF] Argyriou | Convex multi-task feature learning[END_REF] and [START_REF] Rakotomamonjy | l(p)-l(q) penalty for sparse linear and sparse multiple kernel multitask learning[END_REF] for scalar-valued kernels. However, the extension is not straighforward since infinite-dimensional Hilbert spaces with operator-valued reproducing kernels raise some technical issues that we have leveraged in the case of compact operators.

Solving a linear system with multiple block operatorvalued kernels

One common way to construct operator-valued kernels is to build scalar-valued ones which are carried over to the vector-valued (resp. function-valued) setting by a positive definite matrix (resp. operator). In this setting an operator-valued kernel has the following form:

K(w, z) = G(w, z)T,
where G is a scalar-valued kernel and T is an operator in L(G y ). In multi-task learning, T is a finite dimensional matrix that is expected to share information between tasks [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF][START_REF] Caponnetto | Universal multitask kernels[END_REF]. More recently and for supervised functional output learning problems, T is chosen to be a multiplication or an integral operator [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF][START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF]. This choice is motivated by the fact that functional linear models for functional responses [START_REF] Ramsay | Functional Data Analysis[END_REF] are based on these operators and then such kernels provide an interesting alternative to extend these models to nonlinear contexts. In addition, some works on functional regression and structured-output learning consider operator-valued kernels constructed from the identity operator as in [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF] INRIA Algorithm 2 Gauss-Seidel Method choose an initial vector of functions α (0)

repeat for i = 1, 2, . . . , n α (t) i ←-sol. of (13): [K(xi, xi) + λI]α (t) i = si end for until convergence
and [START_REF] Brouard | Semi-supervised penalized output kernel regression for link prediction[END_REF]. In this work we adopt a functional data analysis point of view and then we are interested in a finite combination of operator-valued kernels constructed from the identity, multiplication and integral operators. A problem encountered when working with operator-valued kernels in infinite-dimensional spaces is that of solving the system of linear operator equations [START_REF] Cortes | L 2 regularization for learning kernels[END_REF]. In the following we show how to solve this problem for two cases of operator-valued kernel combinations.

Case 1: multiple scalar-valued kernels and one operator. This is the simpler case where the combination of operator-valued kernels has the following form

K(w, z) = M k=1 d k G k (w, z)T, (11)
where G k is a scalar-valued kernel, T is an operator in L(G y ), and d k are the combination coefficients. In this setting, the block operator kernel matrix K can be expressed as a Kronecker product between the multiple scalar-valued kernel matrix

G = M k=1 d k G k , where G k = [G k (x i , x j )] n i,j=1
, and the operator T . Thus we can compute an analytic solution of the system of equations ( 8) by inverting K + λI using the eigendecompositions of G and T as in [START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF].

Case 2: multiple scalar-valued kernels and multiple operators. This is the general case where multiple operator-valued kernels are combined as follows

K(w, z) = M k=1 d k G k (w, z)T k , (12)
where G k is a scalar-valued kernel, T k is an operator in L(G y ), and d k are the combination coefficients. Inverting the associated block operator kernel matrix K is not feasible in this case, that is why we propose a Gauss-Seidel iterative procedure (see Algorithm 2) to solve the system of linear operator equations [START_REF] Cortes | L 2 regularization for learning kernels[END_REF]. Starting from an initial vector of functions α (0) , the idea is to iteratively compute, until a convergence condition is satisfied, the functions α i according to the following expression

[K(x i , x i ) + λI]α (t) i = 2y i - i-1 j=1 K(x i , x j )α (t) j - n j=i+1 K(x i , x j )α (t-1) j , ( 13
)
where t is the iteration index. This problem is still challenging because the kernel K(•, •) still involves a positive combination of operator-valued kernels. Our algorithm is based on the idea that instead of inverting the finite combination of operator-valued kernels [K(x i , x i )+λI], we can consider the following variational formulation of this system min

α (t) i 1 2 M+1 k=1 K k (x i , x i )α (t) i , α (t) i Gy -s i , α (t) i Gy ,
where

s i = 2y i - i-1 j=1 K(x i , x j )α (t) j - n j=i+1 K(x i , x j )α (t-1) j , K k = d k G k T k , ∀k ∈ {1, .
. . , M }, and K M+1 = λI. Now, by means of a variable-splitting approach, we are able to decouple the role of the different kernels. Indeed, the above problem is equivalent to the following one :

min α (t) i 1 2 K(x i , x i )α (t) i , α (t) i G M y -s i , α (t) i G M y with α (t) i,1 = α (t) i,k for k = 2, . . . , M + 1, where K(x i , x i ) is the (M + 1) × (M + 1) diagonal matrix [K k (x i , x i )] M+1 k=1 . α (t) i is the vector (α (t) i,1 , . . . , α (t)
i,M+1 ) and the (M + 1)-dimensional vector s i = (s i , 0, . . . , 0). We now have to deal with a quadratic optimization problem with equality constraints. Writing down the Lagrangian of this optimization problem and then deriving its first-order optimality conditions leads us to the following set of linear equations

   K 1 (x i , x i )α i,1 -s i + M k=1 γ k = 0 K k (x i , x i )α i,k -γ k = 0 α i,1 -α i,k = 0 ( 14
)
where k = 2, . . . , M + 1 and {γ k } are the Lagrange multipliers related to the M equality constraints. Finally, in this set of equations, the operator-valued kernels have been decoupled and thus, if their inversion can be easily computed (which is the case in our experiments), one can solve the problem [START_REF] Kloft | ℓ p -norm multiple kernel learning[END_REF] with respect to {α i,k } and γ k by means of another Gauss-Seidel algorithm after simple reorganization of the linear system.

Experiments

In order to highlight the benefit of our multiple operator-valued kernel learning approach, we have considered a series of experiments on a real dataset, involving functional output prediction in a brain-computer interface framework. The problem we addressed is a sub-problem related to finger movement decoding from Electrocorticographic (ECoG) signals. We focus on the problem of estimating if a finger is moving or not and also on the direct estimation of the finger movement amplitude from the ECoG signals. The development of the full BCI application is beyond the scope of this paper and our objective here is to prove that this problem of predicting finger movement can benefit from multiple kernel learning. To this aim, the fourth dataset from the BCI Competition IV [START_REF] Miller | Prediction of finger flexion: 4th brain-computer interface data competition[END_REF] was used. The subjects were 3 epileptic patients who had platinium electrode grids placed on the surface of their brains. The number of electrodes varies between 48 to 64 depending on the subject, and their position on the cortex was unknown. ECoG signals of the subject were recorded at a 1KHz sampling using BCI2000 [START_REF] Schalk | BCI2000: a general-purpose brain-computer interface system[END_REF].

A band-pass filter from 0.15 to 200Hz was applied to the ECoG signals. The finger flexion of the subject was recorded at 25Hz and up-sampled to 1KHz by means of a data glove which measures the finger movement amplitude. Due to the acquisition process, a delay appears between the finger movement and the measured ECoG signal [START_REF] Miller | Prediction of finger flexion: 4th brain-computer interface data competition[END_REF]. One of our hopes is that this time-lag can be properly learnt by means of multiple operator-valued kernels. Features from the ECoG signals are built by computing some band-specific amplitude modulation features, which is defined as the sum of the square of the band-specific filtered ECoG signal samples during a fixed time window. For our finger movement prediction task, we have kept 5 channels that have been manually selected and split ECoG signals in portions of 200 samples. For each of these time segments, we have the label of whether at each time sample, the finger is moving or not as well as the real movement amplitudes. The dataset is composed of 487 couples of input-output signals, the output signals being either the binary movement labels or the real amplitude movement. An example of input-output signals are depicted in Figure 1. In a nutshell, the problem boils down to be a functional regression task with functional responses.

To evaluate the performance of the multiple operator-valued kernel learning approach, we use both: (1) the percentage of labels correctly recognized (LCR) defined by (W r /T n ) × 100%, where W r is the number of well-recognized labels and T n the total number of labels to be recognized; (2) the residual sum of squares error (RSSE) as evaluation criterion for curve prediction

RSSE = i {y i (t) -y i (t)} 2 dt, ( 15
)
where y i (t) is the prediction of the function y i (t) corresponding to real finger movement or the finger movement state.

For the multiple operator-valued kernels having the form [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], we have used a Gaussian kernel with 5 different bandwidths and a polynomial kernel of degree 1 to 3 combined with three operators T : identity T y(t) = y(t), multiplication operator associated with the function e -t 2 defined by T y(t) = e -t 2 y(t), and the integral Hilbert-Schmidt operator with the kernel e -|t-s| proposed in [START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF], T y(t) = e -|t-s| y(s)ds. The inverses of these operators can be computed ana- lytically. While the inverses of the identity and the multiplication operators are easily and directly computable from the analytic expressions of the operators, the inverse of the integral operator is computed from its spectral decomposition as in [START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF]. The number of eigenfunctions as well as the regularization parameter λ are fixed using "one-curve-leave-out cross-validation" [START_REF] Rice | Estimating the mean and covariance structure nonparametrically when the data are curves[END_REF] with the aim of minimizing the residual sum of squares error. Empirical results on the BCI dataset are summarized in Tables 1 and2 . The dataset was randomly partitioned into 65% training and 35% test sets. We compare our approach in the case of ℓ 1 and ℓ 2 -norm constraint on the combination coefficients with: (1) the baseline scalar-valued kernel ridge regression algorithm by considering each output independently of the others, (2) functional response ridge regression using an integral operator-valued kernel [START_REF] Kadri | Functional regularized least squares classification with operator-valued kernels[END_REF], (3) kernel ridge regression with an evenly-weighted sum of operator-valued kernels, which we denote by ℓ ∞ -norm MovKL.

As in the scalar case, using multiple operator-valued kernels leads to better results. By directly combining kernels constructed from identity, multiplication and integral operators we could reduce the residual sum of squares error and enhance the label classification accuracy. Best results are obtained using the MovKL algorithm with ℓ 2 -norm constraint on the combination coefficients. RSSE and LCR of the baseline kernel ridge regression are significantly outperformed by the operator-valued kernel based functional response regression. These results confirm that by taking into account the relationship between outputs we can improve performance. This is due to the fact that an operatorvalued kernel induces a similarity measure between two pairs of input/output.

Conclusion

In this paper we have presented a new method for learning simultaneously an operator and a finite linear combination of operator-valued kernels. We have extended the MKL framework to deal with functional response kernel ridge regression and we have proposed a block coordinate descent algorithm to solve the resulting optimization problem. The method is applied on a BCI dataset to predict finger movement in a functional regression setting. Experimental

INRIA

Appendix A Existence of Minimizers

We discuss in this section the existence of minimizers of problems ( 1) and [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and mercer theorem[END_REF].

Because in both problems, we deal with infinite dimensional spaces in the optimization problem, we have to consider appropriate tools for doing so. Existence of f λ in the problem given in Equation ( 1) is guaranteed, for λ > 0 by the generalized Weierstrass Theorem and one of its corollary that we both remind below [START_REF] Kurdila | Convex Functional Analysis[END_REF].

Theorem 2 Let X be a reflexive Banach space and C ⊆ X a weakly closed and bounded set. Suppose h : C → R is a proper lower semi-continuous function. Then h is bounded from below and has a minimizer on C.

Corollary 3 Let H be a Hilbert space and h : H → R is a strongly lower semicontinuous, convex and coercive function. Then h is bounded from below and attains a minimizer.

We can straighforwadly apply this corollary to Problem (1) by defining

h(f ) = n i=1 y i -f (x i ) 2
Gy + λ f 2 F with f ∈ F (which is an Hilbert space). It is easy to note that h is continuous and convex. Besides, h is coercive for λ > 0 since f 2 F is coercive and the sum involves only positive terms. Hence f λ exists.

Regarding the MKL problem given in (6), we show existence of a solution in d and {f k } by defining the function, for fixed {d k }

h 1 (f 1 , • • • , f k ; {d k } k ) = n i=1 y i - k f k (x i ) 2 Gy + λ k f k 2 F d k RR n°7900
and by rewriting problem [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and mercer theorem[END_REF] as

min d∈D J(d) with J(d) = min {f k } h 1 (f 1 , •, , f k ; {d k } k )
Using similar arguments as above, it can be shown that h 1 is proper, strictly convex and coercive for fixed non-negative {d k } and λ > 0 (remind the convention that x 0 = 0 if x = 0 and ∞ otherwise). Hence, minimizers of h 1 w.r.t. {f 1 , • • • , f k } exists and are unique. Since the function J(d) which is equal to

h 1 (f ⋆ 1 , • • • , f ⋆ k ; {d k } k
) is continuous over the compact subset of R M defined by the constraints on d, it also attains its minimum. This conclude the proof that a solution of problem (6) exists.

B Dual Formulation of Functional Ridge Regression

Essential computational details regarding the dual formulation of functional ridge regression presented in Section 2 are discussed here.

The functional response ridge regression optimization problem has the following form:

min f ∈F 1 2 f 2 F + 1 2nλ n i=1 ξ i 2 Gy with ξ i = y i -f (x i ). (16)
where (x i , y i ) i=1,...,n ∈ (G x , G y ). G x and G y are the Hilbert space L 2 (Ω) which consists of all equivalence classes of square integrable functions on a finite set Ω, and F is a RKHS whose elements are continuous linear operators on G x with values in G y . K is the L(G y )-valued reproducing kernel of F .

Since G x and G y are functional spaces, to derive a "dual version" of problem [START_REF] Kurdila | Convex Functional Analysis[END_REF] we use the method of Lagrange multipliers on Banach spaces which is suitable to solve certain infinite-dimensional constrained optimization problems. The method is a generalization of the classical method of Lagrange multipliers. The existence of Lagrangian multipliers for the problem (16) which involves an equality constraint is guaranteed by Theorem 4.11 in [START_REF] Kurcyusz | On the existence and nonexistence of lagrange multipliers in Banach spaces[END_REF]. As consequence, the Lagrangian function associated to ( 16) is well defined and Fréchet-differentiable. Let α = (α i ) i=1,...,n ∈ G n y the vector of functions containing the Lagrange multipliers, the Lagrangian function is given by

L(f, α, ξ) = 1 2 f 2 F + 1 2nλ ξ 2 G n y + α, y -f (x) -ξ G n y , (17)
where

α = (α 1 , . . . , α n ) ∈ G n y , y = (y 1 , . . . , y n ) ∈ G n y , ξ = (ξ 1 , . . . , ξ n ) ∈ G n y , f (x) = (f (x 1 ), . . . , f (x n )) ∈ G n y , and ∀a, b ∈ G n y , a, b G n y = n i=1 a i , b i Gy .
Now, we compute L ′ (f ) the derivative of L(f, α, ξ) with respect to f using the Gâteaux derivative (generalization of the directional derivative) which can be defined for the direction h ∈ F by:

D h L(f ) = lim τ -→0 L(f + τ h) -L(f ) τ Using the fact that D h L(f ) = L ′ (f ), h F , we obtain i. G(f ) = f 2 F lim τ -→0 f + τ h 2 F -f 2 F τ = 2 f, h =⇒ G ′ (f ) = 2f ii. H(f ) = α, y -f (x) -ξ G n y lim τ -→0 α, y -f (x) -τ h(x) -ξ G n y -α, y -f (x) -ξ G n y τ = -α, h(x) G n y = -i α i , h(x i ) Gy = -i K F (x i , •)α i , h F (using the reproducing property) =⇒ H ′ (f ) = -i K F (x i , •)α i (i), (
ii), and L ′ (f ) = 0, we obtain the (representer theorem) solution:

f (•) = n i=1 K F (x i , •)α i (18)
Substituting this into [START_REF] Lanckriet | Learning the kernel matrix with semi-definite programming[END_REF], the problem ( 16) becomes

min ξ max α - 1 2 Kα, α G n y + 1 2nλ ξ 2 G n y + α, y -ξ G n y (19)
where K = [K(x i , x j )] n i,j=1 is the block operator kernel matrix. Differentiating [START_REF] Micchelli | Learning the kernel function via regularization[END_REF] with respect to ξ using the same procedure as described above, we obtain ξ = nλα and then the dual of the functional response ridge regression problem is given by max

α - nλ 2 α 2 G n y - 1 2 Kα, α G n y + α, y G n y (20)

C Convergence of Algorithm 1

In this section, we present a proof of convergence of Algorithm 1. The proof is an extension of results obtained by [START_REF] Argyriou | Convex multi-task feature learning[END_REF] and [START_REF] Rakotomamonjy | l(p)-l(q) penalty for sparse linear and sparse multiple kernel multitask learning[END_REF] to infinite dimensional Hilbert spaces with operator-valued reproducing kernels. Let R(f, d) be the objective function of the MovKL problem defined by ( 6

S(f (n+1) ) ≤ g(f (n) ) ≤ S(f (n) )
where n is the number of iteration. So, the sequence {S(f (n) ), n ∈ N} is nonincreasing and then it is bounded since L is bounded from below. Thus, as n -→ ∞, S(f (n) ) converges to a number which we denote by S * . {S(f (n) )} is convergent and S is a coercive function, then the sequence { f (n) , n ∈ N} is bounded. Consequently, the sequence {f (n) , n ∈ N} is bounded.

Next we show the following subsequence convergence property which underlies the convergence of Algorithm 1.

Proposition 4 If F is a RKHS associated to a compact-operator-valued kernel, the sequence {f (n) ∈ F , n ∈ N}, since it is bounded, has a convergent subsequence.

Proof. The analogue of the Bolzano-Weierstrass theorem2 in Hilbert spaces states that there exists a weakly convergent subsequence {f (n l ) , l ∈ N} of the bounded sequence {f (n) } which (weakly) converges to f ∈ F . By definition of weakly convergence, we have ∀g ∈ F (F is a RKHS with the operator-valued kernel K): lim At this stage, we have shown that any convergent subsequence of {f (n) , n ∈ N} converges to the minimizer of R. Since the sequence {f (n) , n ∈ N} is bounded, it follows that the whole sequence converges to minimizer of R.

n l →∞ f (n l ) (•), g F = f (•), g F (21)

INRIA

Figure 1 :

 1 Figure 1: Example of a couple of input-output signals in our BCI task. (left) Amplitude modulation features extracted from ECoG signals over 5 pre-defined channels. (middle) Signal of labels denoting whether the finger is moving or not. (right) Real amplitude movement of the finger.

f k (x i ) 2

 2 Gy . Substituting Equation[START_REF] Dinuzzo | Learning output kernels with block coordinate descent[END_REF] in R we obtain the objective function:S(f ) := R(f, d(f )) = L +The function S is strictly convex since L is convex and the function definedby f -→ M k=1 f k 2r r+1 r+1 ris strictly convex (this follows directly from strict convexity of the function x -→ x p when x ≥ 0 and p > 1). Thus, S(f ) admits a unique minimizer. Now let us define the function g by:g(f ) = min u {R(u, d(f ))}.The function g is continuous. This comes from the fact that the function:G(d) = min u {R(u, d)} is continuous. Indeed, Gis the minimal of value of a functional response kernel ridge regression problem in a function-valued RKHS associated to an operatorvalued kernel K. So, G(d) = R(d, u * ) with u * = K(d) + λI -1 y (see Equation (8)). u * is continuous, and hence G(d) is also continuous. By definition we have S(f ) = R(f, d(f )), and since d(f ) minimizes R(f, •), we obtain that:

  Let g = K(x, •)β. Using the reproducing property, we obtain limn l →∞ f (n l ) (•), g F = lim n l →∞ f (n l ) (x), β Gy and f (•), g F = f (x), β Gy ⇒ lim n l →∞ f (n l ) (x), β Gy = f (x), β Gy using (21)Thus the subsequence {f (n l ) (•)} is (weakly) pointwise convergent.Now we show that:lim n l →∞ f (n l ) (•) F = f (•) F ( * ) of g and S, g(f * ) = S(f * ). This implies that f * is a minimizer of R(•, d(f * )), because R(f * , d(f * )) = S(f * ). Moreover, d(f * ) is the minimizer of R(•, f * )subject to the constraints on d. Thus, since the objective function R is smooth, the pair (f * , d(f * )) is the minimizer of R.

  this work, we are looking for a solution to this minimization problem in a function-valued reproducing kernel Hilbert space F . More precisely, we mainly focus on the RKHS F whose elements are continuous linear operators on G x with values in G y . The continuity property is obtained by considering a special class of reproducing kernels called Mercer kernels [7, Proposition 2.2]. Note that in this case, F is separable since G x and G y are separable [6, Corollary 5.2].

Table 1 :

 1 Results for the movement state prediction. Residual Sum of Squares Error (RSSE) and the percentage number of Labels Correctly Recognized (LCR) of : (1) baseline KRR with the Gaussian kernel, (2) functional response KRR with the integral operator-valued kernel, (3) MovKL with ℓ ∞ , ℓ 1 and ℓ 2 -norm constraint.

	Algorithm	RSSE LCR(%)
	KRR -scalar-valued -	68.32	72.91
	KRR -functional response -	49.40	80.20
	MovKL -ℓ ∞ norm -	45.44	81.34
	MovKL -ℓ 1 norm -	48.12	80.66
	MovKL -ℓ 2 norm -	39.36	84.72

Table 2 :

 2 Residual Sum of Squares Error (RSSE) results for finger movement prediction.

	Algorithm	RSSE
	KRR -scalar-valued -	88.21
	KRR -functional response -	79.86
	MovKL -ℓ ∞ norm -	76.52
	MovKL -ℓ 1 norm -	78.24
	MovKL -ℓ 2 norm -	75.15
	results show that our algorithm achieves good performance outperforming ex-
	isting methods. It would be interesting for future work to thoroughly compare
	the proposed MKL method for operator estimation with previous related meth-
	ods for multi-class and multi-label MKL in the contexts of structured-output
	learning and collaborative filtering.

RR n°7900

This theorem considers only equality constraint, but it is a particular case of Theorem 3.1 in[START_REF] Kurcyusz | On the existence and nonexistence of lagrange multipliers in Banach spaces[END_REF] which deals with more general context.INRIA

The Bolzano-Weierstrass theorem states that each bounded sequence in R n has a convergent subsequence. For infinite-dimensional spaces, strong convergence of the subsequence is not reached and only weak convergence is obtained. Proposition 4 shows that strong convergence can be reached for the sequence {f(n) , n ∈ N} solution of our MovKL optimization problem in Hilbert spaces with reproducing compact operator-valued kernels.INRIA

Since f (n l ) ∈ F is solution of the minimization of the optimization problem [START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces of integrable functions and mercer theorem[END_REF] with the kernel combination coefficients d k fixed, it can be written as

converges weakly to α ∈ G n y . Indeed, ∀β ∈ G y we have:

, α

(n l ) j

Gy (using the reproducing property)

= lim

Using ( * ) and weak convergence, we obtain the strong convergence of the subsequence {f (n l ) }.

F (using weak convergence) = 0 (using( * ))

⇒ f (n l ) converges strongly to f By Proposition 4, there exists a convergent subsequence {f (n l ) , l ∈ N} of the bounded sequence {f (n) , n ∈ N}, whose limit we denote by f * . Since S(f (n+1) ) ≤ g(f (n) ) ≤ S(f (n) ) , g(f (n) ) converges to S * . Thus, by the continuity 3 Compact operator maps weakly convergent sequences into strongly convergent sequences.

RR n°7900


Alexandre Baksic

Pierre Nataf

Cristiano Ciuti

Superradiant phase transitions with three-level systems

We determine the phase diagram of N identical three-level systems interacting with a single photonic mode in the thermodynamical limit (N → ∞) by accounting for the so-called diamagnetic term and the inequalities imposed by the Thomas-Reich-Kuhn (TRK) oscillator strength sum rule.

The key role of transitions between excited levels and the occurrence of first-order phase transitions is discussed. We show that, in contrast to two-level systems, in the three-level case the TRK inequalities do not always prevent a superradiant phase transition in presence of a diamagnetic term.

The collective superradiant coupling of a large number of atoms (or artificial atoms) has attracted a considerable interest since the pioneering work of Dicke [1] and is now the focus of many recent studies in cavity [2][3][4][5][6] and circuit [7][8][9] quantum electrodynamics. The well-known Dicke model describes the coupling between a collection of two-level systems and a single photon mode. Remarkably, for increasing light-matter coupling such a model predicts a superradiant phase transition [10][11][12], with a doubly degenerate ground state above a critical vacuum Rabi coupling. The so-called superradiant phase is characterized by a spontaneous polarization of the atoms and a spontaneous coherence of the cavity field in the ground state. In the case of time-independent Hamiltonians, photons in the ground state are 'virtual' (i.e., bound in the cavity) and they cannot be radiated out of the cavity [13] unless some non-adiabatic modulation of the Hamiltonian is applied [14] (in analogy with the dynamical Casimir effect [15]). In this sense, the term 'superradiant' currently used for the Hepp-Lieb [10] ground state is unfortunate, because it was originally introduced by Dicke [1] for collective excited state radiative decay and not for ground state properties. Even in spite of no extracavity emission from the ground state, occurrence of a superradiant critical point can be in principle monitored by measuring the dispersion of the collective excitations via standard optical techniques (e.g., through transmission spectroscopy). In the case of non-equilibrium phase transitions for pumped open systems, there is not a true ground state [2][3][4][5][6], but instead a stationary state, which is accompanied by the emission of real photons.

However, for the case of electric dipole transitions, the Dicke model does not include the so-called diamagnetic term, which is proportional to the squared electromagnetic vector potential present in the minimal coupling Hamiltonian describing light-matter interaction in the non-relativistic regime. It is known that in the case of two-level real atoms, such a diamagnetic term is crucial, because it forbids the phase transition as a result of the TRK oscillator strength sum rule [8,16]. We point out that such no-go theorem cannot be necessarily applied to time-dependent Hamiltonians with applied dressing fields [3], to the case of magnetic dipole coupling [17] and to more complex effective systems simulating the Dicke Hamiltonian with different degrees of freedom [3,4]. Moreover, the no-go theorem [8,16] is formulated for twolevel systems, while the study of the multilevel case has been initiated only recently [18][19][20] A generalized model including three-level atoms in the lambda configuration coupled to two photon modes has been theoretically investigated in a recent paper [20] showing a very rich phase diagram with superradiant transitions of both second and first order (in the case of two-level atoms the superradiant transition is of second order). In such an interesting work [20], the diamagnetic term has not been included and the particular two-color lambda configuration has been considered, so the general case of an arbitrary three-level system needs to be explored. In a recent letter [18], it has been claimed that for multilevel atoms coupled to a single photon mode the no-go theorem can be generalized using again the TRK oscillator strength sum rule: in the proof reported in such a work [18], through a perturbative argument, it is assumed that transitions between excited states can be always neglected in the thermodynamical limit: this is rather surprising [19], since in the work by Hayn et al.[20], transitions between excited states play instead a crucial role in the thermodynamical limit and are responsible for the first-order transition boundaries. It is therefore a fundamental problem to explore superradiant phase transitions with arbitrary three-level systems.

In this letter, we investigate the existence of superradiant phase transitions of a system consisting of N three-level systems coupled to a single photon mode including the diamagnetic term in the Hamiltonian and TRK inequalities. We show the rich phase diagram in the thermodynamical limit (N → ∞), obtained via a multilevel Holstein-Primakoff approach. We have studied particular configurations (ladder, lambda, V-type) and also the general three-level coupling configuration. We find that for three-level systems TRK inequalities do not always prevent superradiant phase transitions and that excited state transitions play a crucial role in the thermodynamical limit.

As sketched in Fig. 1 the atomic index. Each atom is assumed to be independent from the others and to interact identically with a single bosonic, photonic mode. The energies of the three levels are ω 0 < ω 1 < ω 2 .

By introducing the collective operators Σij =

H/ = ω cav â † â + 2 j=0 ω j Σjj + D(â + â † ) 2 + 2 i,j=0 (i =j) Ω ij 1 √ N ( Σij + Σji )(â + â † ). (1)
The cavity mode is described by the frequency ω cav and by its creation (annihilation) bosonic operator a † (a). The coupling between the cavity mode and the atomic i → j transition is quantified by the collective vacuum Rabi frequency Ω ij . Ω ij / √ N is the light-matter coupling per atom. The term proportional to D and quadratic in the photon operators is the so-called diamagnetic term. The area below the red-dashed line corresponds to the region satisfying the TRK inequalities. Bottom panel: the black line represents the spatial profile of an illustrative artificial one-dimensional potential with squared wells (as those realizable with semiconductor heterostructures) in units of the energy Ec = 2 2mL 2 where the spatial coordinate is expressed in units of the length L. The horizontal lines depict the energies of the three bound states with their corresponding wavefunctions (other states are in the continuum). For this potential shape, one obtains the oscillator strengths (f01 = 0.3995, f02 = 0.4069, f12 = 0.735) and anharmonicity ratio ω10/ω21 = 0.1709 corresponding to the point depicted by the blue cross in the middle panel.

It originates from the (p-q Â) 2 2m form of the non-relativistic electron-light interaction in the so-called p •  gauge, where p is the electron momentum operator, q the electron charge and  is the electromagnetic vector potential operator. As it can be deduced by the treatment in Ref. [8], the value of the vacuum Rabi frequencies Ω ij are linked to the diamagnetic term amplitude via the oscillator strengths as follows:

Ω 2 ij ω ji = f ij D (2)
where ω ji = ω j -ω i is the transition frequency and

f ij = 2m (ω j -ω i )|d ij | 2 (3)
is the transition oscillator strength with d ij the corresponding atomic electric dipole matrix element (along the direction of the photon mode polarization). Indeed, as shown in Ref. [8] , D = q 2 2m n el N A 2 0 where q is the electron charge, n el the number of electrons per atom and  = A 0 (a+a † ) is the electromagnetic vector potential in the position where the atoms are assumed to be coupled (the spatial variation of the cavity field is neglected in the region occupied by the atoms, i.e., the electric dipole approximation has been considered). For a given collection of atoms and for a given cavity, D is a constant depending on the density of atoms in the cavity volume via the factor N A 2 0 . A related similar model can be obtained by considering a two-dimensional electron gas in a semiconductor quantum well heterostructure [21][22][23][24] (replacing the bare electron mass with the semiconductor conduction band effective mass).

The TRK oscillator strength reads j f ij = 1. Note that by definition the oscillator strengths f ij > 0 if ω j > ω i , while it is negative if ω j < ω i . Hence, for a transition from the ground level |0 to the first excited level |1 , we have always 0 ≤ f 01 ≤ 1 as f 0j ≥ 0. For two-level system models, the presence of the diamagnetic term and the constraint 0 ≤ f 01 ≤ 1 is sufficient to prevent the superradiant phase transition (no-go theorem) [8,16]. Here, we will analyze what happens with three-level systems. By applying the method detailed in Ref. [25] , we can express the generalized collective transition operators by using a multilevel Holstein-Primakoff transformation in terms of bosonic annihilation (creation) operators ĉk ( ĉ † k ). When we use this procedure for the three-level case, we have to choose a state of reference that we will call r. The multilevel boson mapping reads:

Σkk = ĉ † k ĉk (k = r) , Σrr = N - k =r ĉ † k ĉk Σir = ĉ † i N - k =r ĉ † k ĉk , Σij = ĉ † i ĉj (i, j = r) . (4)
The collective transition operators Σij defined above do not commute and are such that Σij , Σkl = δ jk Σilδ il Σkj . The Hamiltonian can be rewritten as

H/ =ω cav â † â + k =r (ω k -ω r )ĉ † k ĉk + ω r N + k =r Ω kr √ N   ĉ † k N - l ĉ † l ĉl + N - l ĉ † l ĉl ĉk   + i>j (i,j =r) Ω ij √ N (ĉ † i ĉj + ĉ † j ĉi ) (â + â † ) + D(â + â † ) 2 (5)
In the following, we choose r = 1 as reference state (choosing another reference state leads to the same results). In order to determine the phase diagram in the thermodynamical limit, we can use the mean-field approach as in Refs. [12,20], by the replacement a = a † = α , c j = c † j = β j with j ∈ {0, 2}. It can be shown that due to the form of Eq. ( 1) the solutions for α and β j are necessarily real numbers. Therefore, we obtain the mean-field expression of the ground state energy ( 5)

E G / =(ω cav + 4D)α 2 + ω 01 β 2 0 + ω 21 β 2 2 + ω 1 N + 4α (Ω 10 β 0 + Ω 21 β 2 ) N -β 2 0 -β 2 2 + Ω 20 β 0 β 2 , (6)
whose global minimization will give the values of α (photonic coherence), β 0 and β 2 (collective atomic coherences) in the ground state. A photonic order parameter α = 0 implies that the system in the Normal (N) phase. If α = 0, a SuperRadiant (SR) phase occurs [12]. In order to minimize the ground state energy in the thermodynamical limit (N → ∞) it is convenient to introduce the rescaled quantities

x = α √ N , y = β0 √ N and z = β2 √ N .
The ground state energy can therefore be rewritten as

E G / N =(ω cav + 4D)x 2 + ω 01 y 2 + ω 21 z 2 + ω 1 + 4x (Ω 10 y + Ω 21 z) 1 -y 2 -z 2 + Ω 20 yz . (7)
As ∂ 2 EG ∂x 2 = 2(ω cav + 4D) > 0, the minimization with respect to x is therefore trivial, giving x as a function of y and z. In order to minimize the ground state energy, we therefore need to minimize with respect to y and z the following energy function (obtained from Eq. ( 7)):

E G =ω 01 y 2 + ω 21 z 2 + ω 1 - 4 ω cav + 4D (Ω 10 y + Ω 21 z) 1 -y 2 -z 2 + Ω 20 yz 2 . (8)
Note that β 2 0 +β 2 2 ≤ N (which is equivalent to y 2 +z 2 ≤ 1) must be fulfilled within the Holstein-Primakoff framework. Hence, the global minimum of E G has to be looked for in a disk with unity radius centered in the origin of the 0yz plane, which can be done with a straightforward numerical calculation.

We start by considering the ladder configuration (f 02 = 0) for the three-level systems. In Fig. 2, we have plotted results for the photonic coherence x = α/ √ N = a / √ N as a function of the oscillator strengths f 01 and f 12 of the ladder coupling for D = 3ω cav , ω 21 = ω cav and ω 10 /ω 21 = 0.1. We point out that these parameters correspond to a strong anharmonicity with the photon mode in resonance with the excited transition. The value of the collective vacuum Rabi frequencies is given by Eq. ( 2). The normal phase occurs when α = 0 (black region). The results in Fig. 2 show that in the considered system a superradiant phase (α = 0) do occur. The phase transition boundary can be of both of first and second order. The first-order phase transition occurs when α has a discontinuous jump from 0 to a finite value (upper part of the frontier). The first-order transition boundary is due to the excited transition 1 → 2 as in the case recently studied by Hayn et al. [20]: at the transition there is a macroscopic occupation of the intermediate state |1 . The second-order phase transition occurs when α is not discontinuous, but its gradient is. The important point to consider here is the fact that due to the TRK sum rule the oscillator strengths are subject to constraints. For the ladder configuration, we have 0 ≤ f 01 ≤ 1 and 0 ≤ f 10 + f 12 ≤ 1, where f 10 = -f 01 ≤ 0. In Fig. 2 the TRK boundaries imposed by such inequalities have been indicated: the area compatible with the TRK sum rule is below the red-dashed line. It is therefore apparent that there is an overlap between the superradiant part of the diagram and the region compatible with the TRK sum rule. Note that the relevant superradiant phase region compatible with the TRK inequalities contains a firstorder transition boundary.

In Fig. 3, we consider instead the V-type configuration, where f 12 = 0, i.e., there is no coupling between the excited states. In such a configuration, the photonic order parameter x = α/ √ N is plotted as a function of f 01 and f 02 (parameters in the caption). The TRK sum rule imposes the inequality 0 ≤ f 01 + f 02 ≤ 1, again indicated by the red-dashed line. In such a V-type configuration, there is a superradiant phase boundary, which is always of second order. In fact, the transition between excited states is by definition inactive in such a configuration. Importantly, we notice that here there is no overlap between the superradiant part of the diagram and the area compatible with the TRK sum rule. Indeed, the TRK sum rule always prevents the superradiant phase transition in the V-type configuration as in the two-level case [8].

In Fig. 4, we show results for the lambda configuration (f 01 = 0): the photonic coherence is shown versus f 12 and f 02 (parameters in the caption). Here, the TRK inequalities 0 ≤ f 12 ≤ 1 and 0 ≤ f 02 ≤ 1, whose boundaries are once again delimited by the red-dashed line. As in the ladder case, in the lambda configuration there is a superradiant phase in the region compatible with the TRK inequalities (are below the red-dashed line) with the transition boundaries being of the first order.

So far, we have shown the simplest three-level configurations (ladder, V-type and lambda). In Fig. 5, we show results for a generic three level system where all the three oscillator strengths are finite (detailed parameters in the caption). The superradiant part satisfying the TRK inequalities is depicted by the red filled volume in the (f 01 , f 02 , f 12 ) space (top panel). The middle panel shows the photonic order parameter on a planar section of such oscillator strength three-dimensional space. We point out that in the generic three-level case the interesting superradiant part compatible with the TRK inequalities is widened due to the additional freedom associated to the third oscillator strengths . In the bottom panel of Fig. 5, we give an illustrative example of a spatial artificial potential providing a three-level system with oscillator strengths (f 01 , f 02 , f 12 ) and anharmonic spectrum corresponding to the cross in the middle panel, for which a superradiant transition is possible while satisfying the TRK constraints. As a perspective, it will be interesting in the future to address more complex multilevel structures and to investigate also the role of direct Coulomb interactions [26].

In conclusion, we have determined the phase diagram for a model system consisting of N three-level systems coupled to a single photonic boson mode, by including the diamagnetic contribution in the light-matter coupling. We have demonstrated that in the considered model system superradiant phase transitions can occur while preserving the TRK inequalities, in stark contrast to the case of two-level systems. We have found that the transition between excited levels have a key role in such superradiant phase transition in contrast to what assumed in Ref. [18]. Our results show that the physics of superradiant phase transitions with multilevel systems might be achievable in a broad range of physical systems, especially those where it is possible to engineer the spectra and oscillator strengths. We would like to thank T. Brandes, I. Carusotto, C. Emary and J. Keeling for stimulating discussions. C. C. is member of Institut Universitaire de France (IUF). We acknowledge support from the ANR project QPOL.

12 FIG. 1 :

 121 FIG.1:A sketch of the considered system consisting of N identical three-level atoms identically coupled to a single cavity mode. Each three-level system is represented by the transition frequencies ω10 = ω1 -ω0 > 0, ω21 = ω2 -ω1 > 0, ω20 = ω21 + ω10 > 0 and by the oscillator strengths f01 ≥ 0, f12 ≥ 0, f02 ≥ 0 (see definition(3) in the text).

FIG. 3 :

 3 FIG. 3: V-type configuration (f12 = 0). Photonic order parameter as a function of the oscillator strengths f01 and f02. The red-dashed lines indicate the boundaries imposed by the TRK sum rule, namely 0 ≤ f01 + f02 ≤ 1. Parameters: D = ωcav, ω10 = ωcav, ω21 = 0.1ωcav. In this configuration, the superradiant phase is always incompatible with the TRK boundaries.

FIG. 5 :

 5 FIG.5: Top panel: the filled red volume in the (f01, f02, f12) space represents the superradiant part of the diagram compatible with the TRK inequalities 0 ≤ f01 + f02 ≤ 1 and 0 ≤ f10 + f12 ≤ 1. Parameters: D = 5ωcav, ω10 = 0.17ωcav, ω21 = ωcav. Middle panel: the photonic order parameter x as a function of f01 and f02 for a fixed value f12 = 0.735. The area below the red-dashed line corresponds to the region satisfying the TRK inequalities. Bottom panel: the black line represents the spatial profile of an illustrative artificial one-dimensional potential with squared wells (as those realizable with semiconductor heterostructures) in units of the energy Ec =