Andreas Neuenkirch

Samy Tindel

A LEAST SQUARE-TYPE PROCEDURE FOR PARAMETER ESTIMATION IN STOCHASTIC DIFFERENTIAL EQUATIONS WITH ADDITIVE FRACTIONAL NOISE

Keywords: November 8, 2011. 2010 Mathematics Subject Classification. 60G15, 62M09, 62F12 fractional Brownian motion, parameter estimation, least square procedure, ergodicity

We study a least square-type estimator for an unknown parameter in the drift coefficient of a stochastic differential equation with additive fractional noise of Hurst parameter H > 1/2. The estimator is based on discrete time observations of the stochastic differential equation, and using tools from ergodic theory and stochastic analysis we derive its strong consistency.

Introduction

In this article, we will consider the following R d -valued stochastic differential equation (SDE)

Y t = y 0 + t 0 b(Y s ; ϑ 0 ) ds + m j=1 σ j B j t , t ∈ [0, T ]. (1)
Here y 0 ∈ R d is a given initial condition, B = (B 1 , . . . , B m ) is an m-dimensional fractional Brownian motion (fBm) with Hurst parameter H ∈ (0, 1), the unknown parameter ϑ 0 lies in a certain set Θ which will be specified later on, {b(•; ϑ), ϑ ∈ Θ} is a family of drift coefficients with b(•; ϑ) : R d → R, and σ 1 , . . . , σ m ∈ R d are assumed to be known diffusion coefficients.

Let us recall that B is a centred Gaussian process defined on a complete probability space (Ω, F , P). Its law is thus characterized by its covariance function, which is defined by

E B i t B j s = 1 2 t 2H + s 2H -|t -s| 2H 1 (i=j) , s, t ∈ R.
The variance of the increments of B is then given by

E |B i t -B i s | 2 = |t -s| 2H , s, t ∈ R, i = 1, . . . , m,
and this implies that almost surely the fBm paths are γ-Hölder continuous for any γ < H. Furthermore, for H = 1/2, fBm coincides with the usual Brownian motion, converting the family {B H , H ∈ (0, 1)} into the most natural generalization of this classical process. In the current paper we assume that the Hurst coefficient satisfies H > 1/2 and we focus on the estimation of the parameter ϑ 0 ∈ Θ. Note that the Hurst parameter and the diffusion coefficients can be estimated via the quadratic variation of Y , see e.g. [START_REF] Bégyn | Quadratic variations along irregular subdivisions for Gaussian processes[END_REF][START_REF] Coeurjolly | Estimating the Parameters of a Fractional Brownian Motion by Discrete Variations of its Sample Paths[END_REF][START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF].

Estimators for the unknown parameter in equation [START_REF] Bégyn | Quadratic variations along irregular subdivisions for Gaussian processes[END_REF] based on continuous observation of Y have been studied e.g. in [START_REF] Belfadli | Parameter Estimation for Fractional Ornstein-Uhlenbeck Processes: Non-ergodic Case[END_REF][START_REF] Hu | Parameter estimation for fractional Ornstein-Uhlenbeck processes[END_REF][START_REF] Kleptsyna | Statistical analysis of the fractional Ornstein-Uhlenbeck type process[END_REF][START_REF] Breton | Filtering and parameter estimation in a simple linear system driven by a fractional Brownian motion[END_REF][START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF][START_REF] Rao | Statistical inference for fractional diffusion processes[END_REF][START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF]. Estimators based on discrete time data, which are important for practical applications, are then obtained via discretization. However, to the best of our knowledge no genuine estimators based on discrete time data have been analyzed yet.

We propose here a least square estimator for ϑ 0 based on discrete observations of the process Y at times {t k ; 0 ≤ k ≤ n}. For simplicity, we shall take equally spaced observation times with t k+1 -t k = κ n -α := α n with given α ∈ (0, 1), κ > 0. We call our method least square-type procedure, insofar as we consider a quadratic statistics of the form

Q n (ϑ) = 1 nα 2 n n-1 k=0 δY t k t k+1 -b(Y t k ; ϑ)α n 2 -σ 2 α 2H n , (2)
where

δY u 1 u 2 := Y u 2 -Y u 1 for any 0 ≤ u 1 ≤ u 2 ≤ T and σ 2 = m j=1 |σ j | 2 .
Let us now describe the assumptions under which we shall work, starting from a standard hypothesis on the parameter set Θ: Hypothesis 1.1. The set Θ is compactly embedded in R q for a given q ≥ 1.

In order to describe the assumptions on our coefficients b, we will use the following notation for partial derivatives:

Notation 1.2. Let f : R d × Θ → R be a C p 1 ,p 2 function for p 1 , p 2 ≥ 1. Then for any tuple (i 1 , . . . i p 1 ) ∈ {1, . . . , d} p 1 , we set ∂ i 1 ...ip 1 f for ∂ p 1 f ∂x i 1 ...∂x ip 1
. Moreover, we will write ∂ x f resp.

∂ ϑ f for the Jacobi-matrices (∂ x 1 f, . . . , ∂ x d f ) and (∂ ϑ 1 f, . . . , ∂ ϑq f ) .
With this notation in mind, our drift coefficients and their derivatives will satisfy a polynomial growth condition, plus an inward condition which is traditional for estimation procedures in the Brownian diffusion case (see e.g [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF][START_REF] Kasonga | The consistency of a non-linear least squares estimator from diffusion processes[END_REF]): (ii) For every x ∈ R d and ϑ ∈ Θ the following growth bounds are satisfied:

Hypothesis 1.3. We have b ∈ C 1,1 (R d × Θ; R d )
|b(x; ϑ)| ≤ c 2 1 + |x| N , |∂ x b(x; ϑ)| ≤ c 2 1 + |x| N , |∂ ϑ b(x; ϑ)| ≤ c 2 1 + |x| N .
As a consequence of the above assumptions on the drift coefficient and the initial condition, for given ϑ 0 ∈ Θ the solution of equation ( 1) converges for t → ∞ to a stationary and ergodic stochastic process (Y t , t ≥ 0), see the next section.

Finally, we also assume that our-drift coefficient is of gradient-type, i.e.:

Hypothesis 1.4. There exists a function

U ∈ C 2,1 (R d × Θ; R) such that ∂ x U(x; ϑ) = b(x; ϑ), x ∈ R d , ϑ ∈ Θ.
With those assumptions in mind, we obtain the following convergence result:

Theorem 1.5. Assume that the Hypotheses 1.1, 1.3 and 1.4 are satisfied for equation [START_REF] Bégyn | Quadratic variations along irregular subdivisions for Gaussian processes[END_REF] and that we moreover have H > 1/2. Let Q n (ϑ) be defined by [START_REF] Belfadli | Parameter Estimation for Fractional Ornstein-Uhlenbeck Processes: Non-ergodic Case[END_REF]. Then we have

sup ϑ∈Θ Q n (ϑ) -E |b(Y 0 ; ϑ 0 )| 2 -E|b(Y 0 ; ϑ)| 2 → 0 (3)
in the P-almost sure sense.

This convergence is in contrast to the case H = 1/2, i.e. to the case of SDEs with additive Brownian noise. There it holds

sup ϑ∈Θ Q n (ϑ) -E |b(Y 0 ; ϑ 0 ) -b(Y 0 ; ϑ)| 2 → 0 (4)
in the P-almost sure sense, see Remark 3.6, and usually the consistent least squares estimator

argmin ϑ∈Θ n-1 k=0 δY t k t k+1 -b(Y t k ; ϑ)α n 2 ,
is considered, see e.g. [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF][START_REF] Kasonga | The consistency of a non-linear least squares estimator from diffusion processes[END_REF]. The difference in the limits (3) and ( 4) is due to the higher smoothness and long-range dependence of fractional Brownian motion for H > 1/2. Our estimator can thus be seen as a "zero squares" estimator instead of a classical least square estimator. In order to show its convergence, we shall work under the following natural identifiability assumption:

Hypothesis 1.6. For any ϑ 0 ∈ Θ, we have

E |b(Y 0 ; ϑ 0 )| 2 = E|b(Y 0 ; ϑ)| 2 iff ϑ = ϑ 0 .
With this additional Hypothesis, the main result of the current article is the consistency of the zero squares estimator based on the statistics Q n : Theorem 1.7. Assume that the Hypotheses 1.1, 1.3, 1.4 and 1.6 are satisfied for equation [START_REF] Bégyn | Quadratic variations along irregular subdivisions for Gaussian processes[END_REF] and let H > 1/2. Let Q n (ϑ) be defined by [START_REF] Belfadli | Parameter Estimation for Fractional Ornstein-Uhlenbeck Processes: Non-ergodic Case[END_REF], and let ϑ n = argmin ϑ∈Θ |Q n (ϑ)|. Then for any ϑ 0 ∈ Θ, we have lim n→∞ ϑ n = ϑ 0 in the P-almost sure sense.

Note that minimizing |Q n (ϑ)| is of course equivalent to finding the zero of Q n (ϑ). Let us shortly compare Theorem 1.7 with the existing literature on estimation procedures for fBm driven equations: (i) Most of the previous results, see e.g. [START_REF] Belfadli | Parameter Estimation for Fractional Ornstein-Uhlenbeck Processes: Non-ergodic Case[END_REF][START_REF] Hu | Parameter estimation for fractional Ornstein-Uhlenbeck processes[END_REF][START_REF] Kleptsyna | Statistical analysis of the fractional Ornstein-Uhlenbeck type process[END_REF][START_REF] Rao | Statistical inference for fractional diffusion processes[END_REF], deal with the one-dimensional fractional Ornstein-Uhlenbeck process in a continuous observation setting. In particular, for this process simple continuous time least-square estimators are obtained in [START_REF] Belfadli | Parameter Estimation for Fractional Ornstein-Uhlenbeck Processes: Non-ergodic Case[END_REF][START_REF] Hu | Parameter estimation for fractional Ornstein-Uhlenbeck processes[END_REF], for which also covergence rates and asymptotic error distributions are derived. Compared to these results our estimation procedure covers a broad class of ergodic multi-dimensional equations and relies on discrete data only. (ii) A general estimation procedure based on moment matching is established in [START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF]. However, the main assumption in [START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF] is that many independent observations of sample paths over a short time interval are available, which is not the case in many practical situations where rather one sample path is discretely observed for a long time period. Let us also mention the article [START_REF] Chronopoulou | On inference for fractional differential equations[END_REF], in which a general discrete data maximum likelihood type procedure has been designed for parameter estimation in both the drift and diffusion coefficients, however without proof of consistency.

(iii) Our current work probably compares best with the maximum likelihood estimator analyzed in [START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF]. The latter pioneering reference focused on one-dimensional SDEs of the form dY t = ϑ 0 h(Y t ) dt + dB t with h : R → R satisfying suitable regularity assumptions. Strong consistency is obtained for the continuous time estimator and also for a discretized version of the estimator. However, the discretized estimator involves rather complicated operators related to the kernel functions arising in the Wiener-integral representation of fBm, which are avoided in our approach. Moreover, in contrast to [START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF] the consistency proof for our estimator does not rely on Malliavin calculus methods. So, in view of the existing results in the literature, Theorem 1.7 can be seen a step towards simple and implementable parameter estimation procedures for SDEs driven by fBm.

Finally, let us comment on the assumptions we have imposed on the drift coefficient and on the Hurst parameter: (a) The hypothesis of Theorem 1.5 are standard for the case H = 1/2, except Hypothesis 1.4 which restricts us to gradient-type drift coefficients. We require this condition to show an ergodic-type result for weighted sums of the increments of fBm, see Lemma 3.2. However, this Hypothesis 1.4 is also implicitly present in the additional condition of [START_REF] Kasonga | The consistency of a non-linear least squares estimator from diffusion processes[END_REF]Theorem 1]. (b) It can easily be shown that whenever ϑ is a one-dimensional coefficient (namely for q = 1), Hypothesis 1.6 is satisfied if the drift coefficient is of the form b(x; ϑ) = ϑh(x) for some h : R d → R d and if the stationary solution is non-degenerate, i.e. we have E|Y 0 | 2 = 0. The latter conditions hold in particular in the case of the ergodic fractional Ornstein-Uhlenbeck process. It would be nice to obtain criteria for richer classes of examples, but this would rely on differentiability and non-degeneracy properties of the map ϑ → E|b(Y 0 ; ϑ)| 2 (see [START_REF] Hairer | A simple framework to justify linear response theory[END_REF] in the Markovian case). We wish to investigate this question in future works. (c) Even if the noise enters additively in our equation, we still need the assumption H > 1/2 in order to prove Theorem 1.7. Indeed, this hypothesis ensures the convergence of some deterministic and stochastic Riemann sums in the computations below (see Remark 3.3 for further details). Whether an adaptation of the proposed zero squares estimator is also convergent in the case H < 1/2 remains an open problem.

Let us finish this introduction with the simplest example of an equation which satisfies the above assumptions: namely the one-dimensional fractional Ornstein-Uhlenbeck process given by

dY t = ϑ 0 Y t dt + dB t , Y 0 = y 0 ∈ R with ϑ 0 < 0.
The solution of this SDE reads as

Y t = y 0 exp(ϑ 0 t) + exp(ϑ 0 t) t 0 exp(-ϑ 0 s) dB s .
For t → ∞ this process converges to the stationary fractional Ornstein-Uhlenbeck process

exp(ϑ 0 t) t -∞
exp(-ϑ 0 s) dB s , t ≥ 0, see e.g. [START_REF] Garrido-Atienza | Discretization of stationary solutions of stochastic systems driven by fractional Brownian motion[END_REF]. Here straightforward computations yield the explicit estimator

ϑ n = n-1 k=0 Y t k δY t k t k+1 n-1 k=0 Y 2 t k α n - n-1 k=0 Y t k δY t k t k+1 n-1 k=0 Y 2 t k α n 2 - n-1 k=0 |δY t k t k+1 | 2 -α 2H n n-1 k=0 Y 2 t k α 2 n .
Notice that even in the case of the fractional Ornstein-Uhlenbeck process, we could neither proof the consistency nor show the inconsistency of our estimator for H < 1/2.

The remainder of this paper is structured as follows: In Section 2 we give some auxiliary results on stochastic calculus for fractional Brownian motion. Section 3 is then devoted to the proof of our main theorems.

Auxiliary Results

Ergodic Properties of the SDE.

To deduce the ergodic properties of SDE (1) we will work without loss of generality on the canonical probability space (Ω, F , P), i.e. Ω = C 0 (R, R m ) equipped with the compact open topology, F is the corresponding Borel-σ-algebra and P is the distribution of the fractional Brownian motion B, which is consequently given here by the canonical process B t (ω) = ω(t), t ∈ R. Together with the shift operators θ t : Ω → Ω defined by

θ t ω(•) = ω(• + t) -ω(t), t ∈ R, ω ∈ Ω,
the canonical probability space defines an ergodic metric dynamical system, see e.g. [START_REF] Garrido-Atienza | Ergodicity of the infinite-dimensional fractional Brownian motion[END_REF]. In particular, the measure P is invariant to the shift operators θ t , i.e. the shifted process (B s (θ t •)) s∈R is still an m-dimensional fractional Brownian motion and for any integrable random variable F : Ω → R we have

lim T →∞ 1 T T 0 F (θ t (ω)) dt = E F
for P-almost all ω ∈ Ω. Owing to the results in Section 4 of [START_REF] Garrido-Atienza | Discretization of stationary solutions of stochastic systems driven by fractional Brownian motion[END_REF] we have the following:

Theorem 2.1. Let Hypothesis 1.3 hold. Then for any ϑ 0 ∈ Θ the following holds: (i) Equation (1) admits a unique solution Y in C λ (R + ; R d ) for all λ < H.

(ii) There exists a random variable

Y : Ω → R d such that lim t→∞ |Y t (ω) -Y (θ t ω)| = 0 for P-almost all ω ∈ Ω. Moreover, we have E|Y | p < ∞ for all p ≥ 1.
Note that the law of Y must coincide with the attracting invariant measure for (1) given in [START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF], see also [START_REF] Hairer | Ergodicity theory of SDEs with extrinsic memory[END_REF][START_REF] Hairer | Ergodicity of hypoelliptic SDEs driven by fractional Brownian motion[END_REF]. Moreover, proceeding as in [START_REF] Garrido-Atienza | Discretization of stationary solutions of stochastic systems driven by fractional Brownian motion[END_REF] we have: Proposition 2.2. Assume Hypothesis 1.3 holds true. Then for any ϑ 0 ∈ Θ and p ≥ 1 there exist constants c p , k p > 0 such that

E |Y t | p ≤ c p , E |Y t -Y s | p ≤ k p |t -s| pH , for all s, t ≥ 0.
The integrability of Y now implies the ergodicity of equation ( 1):

Proposition 2.3. Assume Hypothesis 1.3 holds true. Then for any ϑ 0 ∈ Θ and any

f ∈ C 1 (R d ; R) such that |f (x)| + |∂ x f (x)| ≤ c 1 + |x| N , x ∈ R d ,
for some c > 0, N ∈ N, we have lim

T →∞ 1 T T 0 f (Y t ) dt = Ef (Y ) P-a.s.
Proof. Since the shift operator is ergodic and f has polynomial growth, we have lim

T →∞ 1 T T 0 f (Y (θ t )) dt = Ef (Y ) P-a.s.
Moreover, since lim

t→∞ |Y t (ω) -Y (θ t ω)| = 0
by Theorem 2.1 and f is polynomially Lipschitz, the assertion easily follows.

2.2. Generalized Riemann-Stieltjes Integrals. We set

f ∞;[a,b] = sup t∈[a,b] |f (t)|, |f | λ;[a,b] = sup s,t∈[a,b] |f (t) -f (s)| |t -s| λ
where f : R → R n and λ ∈ (0, 1).

Now, let f ∈ C λ ([a, b]; R) and g ∈ C µ ([a, b]; R) with λ + µ > 1.
Then it is well known that the Riemann-Stieltjes integral b a f (x) dg(x) exists, see e.g. [START_REF] Young | An inequality of Hölder type connected with Stieltjes integration[END_REF]. Also, the classical chain rule for the change of variables remains valid, see e.g. [START_REF] Zähle | Stochastic differential equations with fractal noise[END_REF]: Let f ∈ C λ ([a, b]; R) with λ > 1/2 and F ∈ C 1 (R; R). Then we have

F (f (y)) -F (f (a)) = y a F ′ (f (x)) df (x), y ∈ [a, b]. (5)
Moreover, one has a density type formula: let f, h

∈ C λ ([a, b]; R) and g ∈ C µ ([a, b]; R) with λ + µ > 1. Then for ϕ : [a, b] → R, ϕ(y) = y a f (x) dg(x), y ∈ [a, b],
we have

b a h(x) dϕ(x) = b a h(x)f (x) dg(x). (6)
For later use, we also note the following estimate, which can be found e.g. in [START_REF] Young | An inequality of Hölder type connected with Stieltjes integration[END_REF].

2.3.

The Garcia-Rademich-Rumsey Lemma. We will use the following variant of the Garcia-Rademich-Rumsey Lemma [8]:

Lemma 2.5. Let q > 1, α ∈ (1/q, 1) and f : [0, ∞) → R be a continuous function. Then there exists a constant c α,q > 0, depending only on α, q, such that

|f | q [s,t];α-1/q ≤ c α,q t s t s |f (u) -f (v)| q
|u -v| 1+qα du dv.

2.4.

A Lemma on Pathwise Convergence Rates. The following Lemma (see e.g. [START_REF] Kloeden | The pathwise convergence of approximation schemes for stochastic differential equations[END_REF]), which is a direct consequence of the Borel-Cantelli Lemma, allows us to turn convergence rates in the p-th mean into pathwise convergence rates.

Lemma 2.6. Let α > 0, p 0 ∈ N and c p ∈ [0, ∞) for p ≥ p 0 . In addition, let Z n , n ∈ N, be a sequence of random variables such that

(E|Z n | p ) 1/p ≤ c p • n -α
for all p ≥ p 0 and all n ∈ N. Then for all ε > 0 there exists a random variable η ε such that

|Z n | ≤ η ε • n -α+ε a.s.
for all n ∈ N. Moreover, E|η ε | p < ∞ for all p ≥ 1.

Quadratic Variations of Fractional Brownian Motion.

The following result for the behavior of the quadratic variations of a one-dimensional fractional Brownian motion β with Hurst parameter H is well known, see e.g [START_REF] Tudor | Variations and estimators for self-similarity parameters via Malliavin calculus[END_REF]. Indeed, for H < 3/4 we have

E 1 n n-1 k=0 |δ kk+1 β| 2 -1 2 ≤ c H • 1 n , (7)
while for

H = 3 4 , n > 1, it holds E 1 n n-1 k=0 |δ kk+1 β| 2 -1 2 ≤ c 3/4 • log(n) n . (8)
Finally, if

H ∈ ( 3 4 , 1) then E 1 n n-1 k=0 |δ kk+1 β| 2 -1 2 ≤ c H • 1 n 4-4H . ( 9
)
Here, c H > 0 denotes a constant depending only on H.

3.

Proof of Theorems 1.5 and 1.7

We will denote constants, whose particular value is not important (and which do not depend on ϑ or n) by c, regardless of their value.

Recall that

Q n (ϑ) = 1 nα 2 n n-1 k=0 δY t k t k+1 -b(Y t k ; ϑ)α n 2 -σ 2 α 2H n .
For t ≥ 0, setting

F t = m j=1 σ j B (j)
t , and r k =

t k+1 t k (b(Y u ; ϑ 0 ) -b(Y t k ; ϑ 0 )) du
and moreover using the notation

δ ϑ 0 ϑ b(x) = b(x; ϑ) -b(x; ϑ 0 ), and δF t k t k+1 = F t k+1 -F t k ,
it is readily checked that

Q n (ϑ) = 1 nα 2 n n-1 k=0 |δ ϑ 0 ϑ b(Y t k )| 2 α 2 n - 2 nα 2 n n-1 k=0 δ ϑ 0 ϑ b(Y t k ), δF t k t k+1 α n + 1 nα 2 n n-1 k=0 |δF t k t k+1 | 2 -σ 2 α 2H n + 1 nα 2 n n-1 k=0 |r k | 2 (10)
- 2 nα 2 n n-1 k=0 δ ϑ 0 ϑ b(Y t k ), r k α n + 2 nα 2 n n-1 k=0 δF t k t k+1 , r k .
Note that our assumptions on the drift coefficient imply that

sup ϑ∈Θ |b(x; ϑ) -b(y; ϑ)| ≤ c 1 + |x| N + |y| N • |x -y| for all x, y ∈ R d and |b(x; ϑ 1 ) -b(x; ϑ 2 )| ≤ c 1 + |x| N • |ϑ 1 -ϑ 2 |
for all x ∈ R d and ϑ 1 , ϑ 2 ∈ Θ. So, straightforward estimations using Proposition 2.2 give

E|r k | p ≤ c • α p(1+H) n .
Hence for all p ≥ 1 it holds

E n-1 k=0 |r k | 2 p ≤ c • n p α 2p(1+H) n ,
and Lemma 2.6 implies

lim n→∞ 1 nα 2 n n-1 k=0 |r k | 2 = 0 P-a.s. ( 11
)
Using Proposition 2.2 and Lemma 2.6 again, it follows similarly

lim n→∞ sup ϑ∈Θ 2 nα 2 n n-1 k=0 δ ϑ 0 ϑ b(Y t k ), r k α n = 0 P-a.s. ( 12
)
and, since H > 1/2, we also have

lim n→∞ 2 nα 2 n n-1 k=0 δF t k t k+1 , r k = 0 P-a.s. ( 13
)
Plugging relations ( 11)-( 13) into [START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF], we have obtained that

Q n (ϑ) = Q (1) n (ϑ) -2Q (2) n (ϑ) + Q (3) n + R n (ϑ), (14)
where lim n→∞ sup ϑ∈Θ |R n (ϑ)| = 0 in the P-almost sure sense and

Q (1) n (ϑ) = 1 n n-1 k=0 |δ ϑ 0 ϑ b(Y t k )| 2 , Q (2) n (ϑ) = 1 nα n n-1 k=0 δ ϑ 0 ϑ b(Y t k ), δF t k t k+1
and

Q (3) n = 1 nα 2 n n-1 k=0 |δF t k t k+1 | 2 -σ 2 α 2H n .
The treatment of the terms Q

n (ϑ), Q (1)
n (ϑ) and Q

n will be carried out in the following series of Lemmata. We first show a discrete version of Proposition 2.3:

Lemma 3.1. Let f ∈ C 1,1 (R d × Θ; R d ) be a function such that |f (x; ϑ)| ≤ c 1 + |x| N , |∂ x f (x; ϑ)| ≤ c 1 + |x| N , |∂ ϑ f (x; ϑ)| ≤ c 1 + |x| N for some c > 0, N ∈ N, independent of ϑ ∈ Θ. Then we have sup ϑ∈Θ 1 n n-1 k=0 |f (Y t k ; ϑ)| 2 -E|f (Y ; ϑ)| 2 → 0 P-a.s.
In particular, we have

sup ϑ∈Θ Q (1) n (ϑ) -E|δ ϑ 0 ϑ b(Y )| 2 → 0 P-a.s.
Proof. Let T n = nα n and set

V n (ϑ) = 1 T n Tn 0 |f (Y s ; ϑ)| 2 ds.
The ergodicity of Y yields that there exists a set A 1 ∈ F with full measure such that

lim n→∞ V n (ϑ)(ω) = E|f (Y ; ϑ)| 2
for all ϑ ∈ Θ ∩ Q q and all ω ∈ A 1 . The assumptions on f give

|V n (ϑ 1 ) -V n (ϑ 2 )| ≤ c • 1 + 1 T n Tn 0 |Y s | 2N ds • |ϑ 1 -ϑ 2 |, (15)
so V n is Lipschitz continuous in ϑ and thus

sup ϑ∈Θ V n (ϑ) -E|f (Y ; ϑ)| 2 = sup ϑ∈Θ∩Q q V n (ϑ) -E|f (Y ; ϑ)| 2 .
However, from [START_REF] Istas | Quadratic variations and estimation of the local Hölder index of a Gaussian process[END_REF] and the ergodicity of Y , it also follows that there exists a set A 2 ∈ F with P(A 2 ) = 1 in which the family of random functions V n : Θ → R, n ∈ N, is equicontinuous, and hence the Arzela-Ascoli Theorem yields the desired uniform convergence, i.e.

lim n→∞ sup ϑ∈Θ∩Q q V n (ϑ) -E|f (Y ; ϑ)| 2 = 0 P-a.s. (16)
Setting

G n (t; ϑ) = |f (Y t ; ϑ)| 2 -|f (Y t k ; ϑ)| 2 , t ∈ [t k , t k+1 ), k = 0, 1, . . . , it remains to show that 1 T n Tn 0 sup ϑ∈Θ |G n (t; ϑ)| dt → 0 P-a.s.
To this aim, the assumptions on f imply that

sup ϑ∈Θ |G n (t; ϑ)| ≤ c • (1 + Y 2N t + Y 2N t k ) • |Y t -Y t k |.
Using Proposition 2.2 and Hölder's inequality we obtain

sup t≥0 E sup ϑ∈Θ |G n (t; ϑ)| p ≤ c • α pH n ( 17
)
for all p ≥ 1. Now, Jensen's inequality gives for n → ∞.

E 1 T n Tn 0 sup ϑ∈Θ |G n (t; ϑ)| dt p ≤ 1 T n Tn 0 E sup ϑ∈Θ |G n (t; ϑ)| p dt, so (17) yields
We have a similar ergodic result for weighted sums of the increments of the process F .

Lemma 3.2. Let f ∈ C 1,1 (R d × Θ; R d ) be a function such that |f (x; ϑ)| ≤ c 1 + |x| N , |∂ x f (x; ϑ)| ≤ c 1 + |x| N , |∂ ϑ f (x; ϑ)| ≤ c 1 + |x| N
for some c > 0, N ∈ N, independent of ϑ ∈ Θ. Assume moreover that there exists a function

U ∈ C 2,1 (R d × Θ; R) such that ∂ x U(x; ϑ) = f (x; ϑ), x ∈ R d , ϑ ∈ Θ,
i.e. f is of gradient type. Then, for H > 1/2, we have

sup ϑ∈Θ 1 nα n n-1 k=0 f (Y t k ; ϑ), δF t k t k+1 + E b(Y ; ϑ 0 ), f (Y ; ϑ) → 0 P-a.s.
In particular,

sup ϑ∈Θ Q (2) n (ϑ) + E b(Y ; ϑ 0 ), δ ϑ 0 ϑ b(Y ) → 0 P-a.s.
Proof. Let T n = nα n . First note that the chain of variable and density formula for Riemann-Stieltjes integrals, see ( 5) and( 6) in Subsection 2.2, gives that

1 T n (U(Y Tn ; ϑ) -U(y 0 ; ϑ)) = 1 T n Tn 0 f (Y u ; ϑ), b(Y u ; ϑ 0 ) du + 1 T n Tn 0 f (Y u ; ϑ), dF u .
Now the properties of f , Proposition 2.2 and Lemma 2.6 imply that

sup ϑ∈Θ 1 T n |U(Y Tn ; ϑ) -U(y 0 ; ϑ)| → 0 P-a.s.
Moreover, we have

sup ϑ∈Θ 1 T n Tn 0 f (Y u ; ϑ), b(Y u ; ϑ 0 ) du -E f (Y ; ϑ), b(Y ; ϑ 0 ) → 0 P-a.s.,
which can be derived completely analogously to [START_REF] Kasonga | The consistency of a non-linear least squares estimator from diffusion processes[END_REF]. It follows

sup ϑ∈Θ 1 T n Tn 0 f (Y u ; ϑ), dF u + E f (Y ; ϑ), b(Y ; ϑ 0 ) → 0 P-a.s.
So, it remains to show that

sup ϑ∈Θ 1 T n Tn 0 G n (t; ϑ), dF t → 0 P-a.s. ( 18
)
where

G n (t; ϑ) = f (Y t ; ϑ) -f (Y t k ; ϑ), t ∈ [t k , t k+1 ), k = 0,
E|B (j) | p λ;[t k ;t k+1 ] 1/p ≤ c • α H-λ n . Since moreover sup t∈[t k ,t k+1 ] (1 + |Y t | N ) ≤ c • Y N t k + α λN n • |Y | N λ;[t k ,t k+1 ] and sup t≥0 E|Y t | p < ∞ for all p ≥ 1, it follows that E sup ϑ∈Θ 1 T n Tn 0 G n (t; ϑ), dF t p 1/p ≤ c • α 2H-1 n (19)
Now Lemma 2.6 implies [START_REF] Breton | Filtering and parameter estimation in a simple linear system driven by a fractional Brownian motion[END_REF], since H > 1/2.

Remark 3.3. As mentioned in the introduction, the condition H > 1/2 is used in our proofs. Specifically, it is invoked in the convergence of the weighted stochastic integral [START_REF] Kloeden | The pathwise convergence of approximation schemes for stochastic differential equations[END_REF] and also in order to derive [START_REF] Hairer | Ergodicity of hypoelliptic SDEs driven by fractional Brownian motion[END_REF].

The following Lemma deals with the remaining term, i.e. the quadratic variations of the process F : Lemma 3.4. We have

lim n→∞ Q (3) n = lim n→∞ 1 nα 2 n n-1 k=0 |δF t k t k+1 | 2 -σ 2 α 2H n = 0 P-a.s. with σ 2 = m j=1 |σ j | 2 . Proof. We have |δF t k t k+1 | 2 -σ 2 α 2H n = m j=1 |σ j | 2 |δB (j) t k t k+1 | 2 -α 2H n =I (1) k + m i,j=1, i =j σ i , σ j δB (i) t k t k+1 δB (j) t k t k+1 =I (2) k
.

Owing to the scaling property of fBm it follows that

E n-1 k=0 I (1) k p = α 2Hp n E n-1 k=0 m j=1 |σ j | 2 |δB (j) kk+1 | 2 -1 p .
Since all moments of random variables in a finite Gaussian chaos are equivalent, it follows from ( 7)-( 9) that

E n-1 k=0 |δB (j) kk+1 | 2 -1 p ≤ c • | log(n)| p • n p(2H-1)
and consequently 1) .

E n-1 k=0 I (1) k p ≤ c • α 2Hp n • | log(n)| p • n p(2H-
Since α n = κ • n -α with α ∈ (0, 1) Lemma 2.6 now implies that

lim n→∞ 1 nα 2 n n-1 k=0 I (1) k
= 0 P-a.s.

So it remains to consider the off-diagonal terms, i.e. I

k . Here we can exploit the following trick: Let β and β be two independent fractional Brownian motions with the same Hurst index. From ( 7)-( 9) we clearly have that 1) .

V n = n-1 k=0 |δ t k t k+1 β| 2 -|δ t k t k+1 β| 2 satisfies E|V n | p ≤ c • α 2Hp n • | log(n)| p • n p(2H -
However, setting B (i) = (β + β)/ √ 2 and B (j) = (β -β)/ √ 2, then B (i) and B (j) are two independent fractional Brownian motions and (1) With probability one, L n (ϑ) → L(ϑ) uniformly in ϑ ∈ Θ as n → ∞.

V n L = 2
(2) The limit L is non-random and L(ϑ 0 ) ≤ L(ϑ) for all ϑ ∈ Θ.

  and there exist constants c 1 , c 2 > 0 and N ∈ N such that: (i) For every x, y ∈ R d and ϑ ∈ Θ we have b(x; ϑ) -b(y; ϑ), x -y ≤ -c 1 |x -y| 2

Proposition 2 . 4 .

 24 Let f, g as above. There exists a constant c λ,µ (independent of a, b) such that b a (f (s) -f (a))dg(s) ≤ c λ,µ |f | λ;[a,b] |g| µ;[a,b] |b -a| λ+µ holds for all a, b ∈ [0, ∞).

  t; ϑ)| dt → 0 P-a.s.

n- 1 k=0δ 5 .

 15 t k t k+1 B (i) δ t k t k+1 B(j) . Now we can easily conclude that lim Let us go back to expression[START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF] and[START_REF] Hu | Parameter estimation for fractional Ornstein-Uhlenbeck processes[END_REF]. Applying Lemma 3.1 and 3.2 we obtain thatlim n→∞ sup ϑ∈Θ Q (1) n (ϑ) -2Q (2) n (ϑ) -E|b(Y ; ϑ)| 2 -E|b(Y ; ϑ 0 )| 2 = 0almost surely. Furthermore, recall that Lemma 3.4 asserts that lim n→∞ Q(3) n = 0 almost surely. The proof is now finished.Proof of Theorem 1.7. Let us first recall the following result ([START_REF] Frydman | A proof of the consistency of maximum likelihood estimators of non-linear regression models with autocorellated errors[END_REF][START_REF] Kasonga | The consistency of a non-linear least squares estimator from diffusion processes[END_REF]): Proposition 3.5. Assume that the family of random variables L n (ϑ), n ∈ N, ϑ ∈ Θ, satisfies:

( 3 )nα 2 n n- 1 k=0 |δF t k t k+1 | 2 -σ 2 α n + 1 nα 2 n

 32122 It holds L(ϑ) = L(ϑ 0 ) if and only if ϑ = ϑ 0 . Then, we have ϑ n → ϑ 0 P-a.s. for n → ∞, whereL n ( ϑ n ) = min ϑ∈Θ L n (ϑ).The strong consistency of the zero squares estimator follows now from Theorem 1.5 and an application of Proposition 3.5 to |Q n (ϑ)|.Remark 3.6. In the case H = 1/2 we have under similar assumptions thatsup ϑ∈Θ Q n (ϑ) -Q n (ϑ 0 ) -E |b(Y 0 ; ϑ 0 ) -b(Y 0 ; ϑ)| 2 → 0in the P-almost sure sense, see e.g.[START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF]. SinceQ n (ϑ 0 ) = 1 t k t k+1 , r k ,an application of Lemma 3.1 and 3.4 (which are also valid for H = 1/2) yield thatlim n→∞ Q n (ϑ 0 ) = lim t k+1 , r k P-a.s.However, using the Itô-isometry and the Burkholder-Davis-Gundy inequality we haver k , δ t k t k+1 F p ≤ c • 1 n p α 2p n • T 1+p/2 n α 3p/2 n ≤ c • n -p/2+1-αand Lemma 2.6 thus gives lim n→∞ Q n (ϑ 0 ) = 0 P-a.s.Hence we end up withsup ϑ∈Θ Q n (ϑ) -E |b(Y 0 ; ϑ 0 ) -b(Y 0 ; ϑ)| 2 → 0,so the limit of the statistics Q n is different for H = 1/2, where one obtains the standard least square estimator.

  1, . . . . Applying Proposition 2.4 and using the polynomial Lipschitz continuity of f yields, for all λ < H, |Y t | N )|Y | λ;[t k ;t k+1 ] |B (j) | λ;[t k ;t k+1 ] .

	Tn	m	n-1
	0 (1 + From the Garcia-Rademich-Rumsey inequality, see Lemma 2.5, and Proposition 2.2 we G n (t; ϑ), dF t ≤ c • α 2λ n • j=1 k=0 sup t∈[t k ,t k+1 ]
	have that and also	E|Y | p λ;[t k ;t k+1 ]	1/p	n ≤ c • α H-λ


Hebrew in the Universities 1

Pablo Kirtchuk

Hebrew has been studied in the most important Universities ever since the Renaissance, though the process leading to its institutionalization had begun already in the Middle Age. In those universities as well as in more recent ones it is still studied as a language of culture, prestige and research. In addition, many other offer courses of Contemporary Hebrew without necessarily disposing of dedicated departments. By the middle of the XIII century, Roger Bacon of Oxford (c. 1240-90) advocated the study of Hebrew for other than missionary purposes, and though he head an imperfect command of the language he wrote a grammar of it. Cambridge sees its first important Hebraist with Henry de Costesy, appointed reader in 1326, who had probably learned under Nicholas de Lyra in Paris and was interested in the literal meaning of Scripture. One of his successors was Richard Brinkley (1480Brinkley ( -1518)). Immanuel Tremellius (1510-80), a Swiss Protestant, taught Hebrew in Cambridge whence, upon Queen Mary's accession, he fled to Heidelberg, where he published some scholarly books including a Latin translation to Targum Jonathan of the twelve minor prophets; an Aramaic and Syriac Grammar, and his magnum opus, a Latin translation of the Hebrew Bible. To acquire official status, however, Hebrew studies in England had to wait until the first third of the XVI century, when they were promoted by Richard Fox, who founded Corpus Christi College in Oxford, and by John Fisher, founder of St John's in Cambridge. He had learnt the language from Robert Wakefield, the first English hebraist worthy of the name, who had taught in some of the most famous universities: Paris, Louvain and Tübingen. At the dissolution of the Monasteries he obtained a Hebrew dictionary compiled by Laurence Holbeach (d. c. 1420) on the basis of work begun by a former prior, Gregor of Huntingdon (d. c. 1290). In 1524, Fisher delivered an oration in praise of Semitic studies under the title Oratio de laudibus et utilitate trium linguarum, Arabicae, Chaldaicae et Hebraicae, and six years later he produced the Syntagma de hebraeorum codicum incorruptione. In addition to promoting their study he proved their kinship on grammatical criteria, and showed some Latin and Greek words to be borrowings from Hebrew. In Israel itself, where it is the main language and one of the three official ones along with Arabic and English, there are departments of Hebrew language at the 5 universities (Jerusalem, Tel Aviv, Haifa, Beer Sheva and Bar Ilan); even the Haifa Technion has a small section devoted to Hebrew, the Weizmann institute being the only Scientific Institution without any trace of Hebrew presence (even most of its publications and courses are in English). The thirty-odd colleges also teach Hebrew language to Hebrew teachers in the making. The departments of Hebrew Language in Israeli Universities count philologists and linguists who do teaching and research on Hebrew in several contexts: Semitic linguistics, so-called Jewish languages, General linguistics, Typology, Synchrony, Diachrony as well as specialists in specific fields such as Phonetics, Semantics, Morphology, Syntax, Pragmatics, Discourse, Language contact and Language acquisition. Tel Aviv University leads the constitution of a Corpus of Contemporary Hebrew, Ben Gurion University in Beer Sheva leads one on Biblical Hebrew lexicography, and the Hebrew Language Academy, which along with its official status as a regulating and normalizing institution concerning Israeli Hebrew carries out 5 investigation as well, is embarked in the most ambitious endeavor of the Historical Dictionary, an on-line ever increasing publication of huge dimensions that will eventually be the most complete source of lexico-semantic, grammatical and textual information on Hebrew words and roots through the ages.

  It is probably there that Calvin studied Hebrew during his stay in Paris (1530-3) under the influence of his cousin Robert Olivetan, the first reformer to publish a translation of the Hebrew Bible into French. In 1751, a Chair of Hebrew was created in the Sorbonne by the Duc d'Orléans, 'for the purpose of reviving Oriental learning in the University of Paris and explaining the Hebrew Scriptures'. In this Parisian atmosphere of intellectual freedom and cultural opportunities was educated the Scottish Hebraist Alexander Geddes, a great and open-minded scholar. In 1533 Calvin fled to Basle and continued his training there with Münster and another renowned Hebraist, Simon Grynaeus. During his exile from Geneva (1538-41), Calvin spent most of his time in Strasbourg where Hebrew was taught by Wolfgang Capito and Martin Bucer. A third hebraist was Gregor Caselius. The city of Strasbourg adopted a University program designed by I. Sturm (1507-89), a graduate of the trilingual school in Louvain. Greek and Latin were obligatory and Hebrew was an option eventually chosen by many students. Hebrew scholarship in the Christian world, studied there. Alonso de Zamora's articles on grammar and lexicography appended to the polyglot Bible whose edition he supervised proved that he was an accomplished Hebraist, whose intellectual honesty led him to prefer the Bible translation of Pagninus (1470-1536) to Jerôme's Vulgate. The University of Louvain, Flanders, inaugurated its 'Collegium Trilinguie' in 1520, under the guidance of Jerôme Busleiden. The first Hebrew teacher was Mathaeus Adrianus, one of the first Christian scholars to draw up a systematic course of instruction in Hebrew, and Louvain was the example of what could be actually achieved in order to give Hebraic studies the high place they deserve. In the Netherlands, Hebrew was included in the syllabus of theology students at the University of Leyde from its foundation in 1575, when the Chair was conferred to Johann Drusius, whose edition of Martinius's Grammar was the first Hebrew book to be published in the country. Three other outstanding Dutch Hebraists of the time were Thomas Erpenius (1584-1624), Sixtinus Amama (1593-1629) andConstantijn L'Empereur (1591-1648). Incidentally, the Compendium grammatices linguae hebraeae by Baruch Spinoza (1632-77) appeared in the Netherlands at that period. On the Mayflower there were two Hebraists of competence, William Bradford and William Brewster. The Puritans of New England were avid students of the Hebrew Bible, and the most famous of their scholars, Cotton Mather, used many Hebrew words and phrases in his prolific writings. This enthusiasm led to the inclusion of the language in the curricula of the ten American Colleges founded before 1776. At Harvard, where the first two presidents, Dunster and Chauncey were scholars of Hebrew, all students had to study the language ever since its foundation in 1636. A converted Jew, Judah Monis, appointed lecturer of Hebrew in 1722, published in Boston in 1735 a Grammar of the Hebrew Tongue, considered the first of its kind meant for instruction at an American university. Hebrew ceased to be an obligatory subject at Harvard in 1787. Yale, founded in 1701, offered Hebrew in its early years. It was to become a center of interest in Hebrew under the presidency of Ezra Stiles, who assumed the Chair in 1777, but it was no longer a required subject after 1789. A group of Yale graduates founded Darmouth College in 1769, where Hebrew was given a special prominence until 1828; it was re-established in the high position it deserves in the 1980s. Seventy years after it appeared, it was deemed relevant enough as to be published in an augmented English version made by T. Muraoka(Leyde University, 1991-6). By the end of the 19 th century the SPH model was well established in Harvard, Columbia, Penn State, Johns Hopkins, Chicago, and Berkeley; (3) American Wissenschaft Hebrew, (AWH) which follows the tradition of the Wissenschaft des Judentums, and whose most outstanding representatives were S. Baron and H. Wolfsohn; (4) Tarbut Ivrit Hebrew (TIH), named after the Tarbut movement, which regarded Hebrew chiefly as a Modern language; TIH promoters were chiefly Jews driven to America by Hitler who were attached to Jewish culture but not particularly to Judaism as such; (5) Area Studies Hebrew (ASH), which follows the great amount of research made by H. Rosén and his students on all stages of Hebrew as an object of modern linguistic investigation; (6) Jewish Studies Hebrew (JSH) which offers a holistic view of the language with a rather theoretical focus; since many students are admitted into Jewish Studies without being required to actually learn Hebrew in any of its forms, and finally (7) Israeli Hebrew (IH) which only teaches the language as it is spoken today, with no attention whatsoever to earlier layers -while of course, Israeli Hebrew, unlike most other languages, did not develop from older layers but is a projection thereof and to this extent, paradoxically as it may seem, in order to master thoroughly Israeli Hebrew one cannot but have at least a glimpse at earlier stages. Hebrew is offered by most of the important Universities across Europe, including Eastern Europe and Russia, where it is taught and investigated in world famous institutions such as St Petersburg's Oriental Institute, the Academy of Sciences in Moscow and the University of Novosibirsk among others. France has departments of Hebrew at the Universities of Paris 3, Paris 4, Paris 8, l'Institut National des Langues et Civilisations Orientales, Strasbourg 2, Lyon 3, Aix-Marseille and Lille 3. In Australia it is taught and investigated mainly at Monash University in Melbourne, but also at La Trobe (Victoria) among others. In North America, it is taught in 606 colleges and universities. In Latin America, it is taught mainly at the main universities (Argentina: Buenos Aires and Córdoba; Brazil: Sao Paulo, Campinas and Rio de Janeiro, Chile: Santiago). In the Arab world, Hebrew is taught in several universities of Morocco, Egypt, Syria, Lebanon and Irak, among others, mainly in the context of Ancient History. An Israeli organization is dedicaed to assist European universities in appointing Readers (or Lectors) of Contemporary Hebrew. In France, for example, such positions exist at the above mentioned universities as well as at Lyon 2, Montpellier, Bordeaux and Toulouse. Such positions exist outside of Europe too, e.g. in Japan They are devoted to the teaching of Contemporary Hebrew and are usually filled by Israeli Ph.D. students.

	1547), Agathias Guidacerius (1477-1540) and Paul Paradis (d. 1549). One of their
	immediate successors, Quanquarbre, published in 1554 De grammatica Hebraeorum
	syllabus for divinity students at both universities. Permanent Chairs of Hebrew were created at Oxford and Cambridge as a result of Henry VIII's split with the Catholic Church, which brought the universities increasingly under royal control. In 1536 an act of parliament was approved leading opus. At North American universities the study of Hebrew is nowadays organized according to
	eventually to the creation of those Chairs, and Hebrew was taught on the pattern of the seven different models (Band 1993): (1) Divinity School Hebrew (DSH), which is
	trilingual Colleges at Alcalá and Louvain. Thomas More himself was aware of the Protestant-theologically oriented; (2) Semitic Philology Hebrew (SPH), whose roots
	importance of Hebrew and he encouraged all students of both the Hebrew Bible and the harken back to 19 th century Germany where the commanding figure was W. Gesenius
	Gospel to go ad purissimas fontes, i.e. to the Hebrew and Greek originals. (1786-1842) author of both a Hebrew and Aramaic Lexicon (1810-12) of the Hebrew
	The Collège de France, founded by king François I in 1530 for the promotion of the Bible as well as a Grammar (1813), an English version of the Lexicon by Brown, Driver
	Humanities, became the leading institution in Europe as far as the teaching of Hebrew & Briggs appeared in 1907 and one of the Grammar, by A.E. Cowley (1910) ;. H. Ewald
	was concerned. It boasted no less than three teachers of Hebrew: François Vatable (d. (1803-1875) had published one in 1838 and E. König (1846-1936) puiblished his in 1881.
	1

The statuses of Oxford and Cambridge gave an important place to the study of Hebrew in the curriculum, beside that of Greek and Latin, and in 1535 Melanchton's works replaced those of Duns Scotus as part of the In 1559, Calvin provided his hometown Geneva with a school and an Academy, whose principal was Theodor Beza, a promoter of Hebrew whose Icones, pen-portraits of famous scholars, which he dedicated in 1580 to James VI of Scotland, included six Christian Hebraists, the brightest of whom was famous Johannes Reuchlin. The Chair of Hebrew at Calvin's Academy was filled by A. Chevalier, who gave eight hours of instruction per week, five of which were devoted to Hebrew grammar. The Academy's library possessed number of valuable Hebrew books, including three copies of Daniel Bomberghi's rabbinic Bible edited in Venice. Eventually, this Academy became the leading institution in Europe for the study of Hebrew during the whole of the XVII century and it drew students from all around the continent. A renowned Huguenot Hebraist in this context was the poet Agrippa d'Aubigné (1552-1630), who could read Biblical and post-Biblical Hebrew 'without points' (= vowel signs). Yet it is in Alcalá de Henares, Spain, with the foundation in 1498 by Francisco Ximénez de Cisneros of the Colegio de San Ildefonso, that the first attempt was made towards the patterns of trilingual instruction, which eventually spread around the continent. In 1508 the first professor of Hebrew, Alonso de Zamora, was appointed. He was a recognized authority in all things Hebrew and made the University of Alcalá an acknowledged centre in the field. Benito Arias Montano (1527-98), an outstanding figure of However, the study of Hebrew as it exists presently at universities in Europe and beyond reflects Ernest Renan's (France, 1823-1892) scientific position, which maintained that the language can and indeed should be studied as such, independently from any theological consideration. In 1855, he published his Histoire générale et système comparé des langues sémitiques. Thereafter he was appointed Professor at the Collège de France and promptly removed since in his opening lesson he dare call Jesus an admirable man. Both facts are significant. on one hand the man who, despite his Christian education, rejects the apparently indissoluble link between Hebrew and the dogma; on the other hand, the Establishment, which finds this impossible to accept. Even in our day the Chair at the College de France is common to Hebrew and Aramaic, as is its vis-à-vis at another prestigious Parisian High Education institution, L'Ecole Pratique des Hautes Etudes. Eventually, however, Renan triumphed: he was rehabilitated and appointed director of the same

University, and his compatriot Paul Joüon (1871-1940)

published in 1923 his Grammaire de l'Hébreu Biblique, a masterpiece of linguistic thought, not in the least influenced by religious considerations despite the author's being a clergyman.

This article is based to a fair extent on Kirtchuk

(see Bibliography).


Keywords: D.1.3 [Concurrent Programming]: Parallel Programing, I.6.6 [Simulations and Modeling]: Simulation Output Analysis, E.5 [Files]: Optimization Design, Experimentation, Performance Exascale Computing, Multicore Architectures, I/O, Dedicated Cores

SRC: Damaris -Using Dedicated I/O Cores for Scalable Post-petascale HPC Simulations

Matthieu Dorier ENS Cachan, Brittany -IRISA Rennes, France matthieu.dorier@eleves.bretagne.ens-cachan.fr

Advisor: Gabriel Antoniu, INRIA Rennes, France, gabriel.antoniu@inria.fr

As we enter the post-petascale era, scientific applications running on large-scale platforms generate increasingly larger amounts of data for checkpointing or offline visualization, which puts current storage systems under heavy pressure. Unfortunately, I/O scalability rapidly fades behind the increasing computation power available, and thereby reduced the overall application performance scalability. We consider the common case of large-scale simulations who alternate between computation phases and I/O phases. Two main approaches have been used to handle these I/O phases: 1) each process writes an individual file, leading to a very large number of files from which it is hard to retrieve scientific insights; 2) processes synchronize and use collective I/O to write to the same shared file. In both cases, because of mandatory communications betweens processes during the computation phase, all processes enter the I/O phase at the same time, which leads to huge access contention and extreme performance variability.

Previous research efforts have focused on improving each layer of the I/O stack separately: at the highest level scientific data formats like HDF5 [START_REF] Ncsa | Hierarchical Data Format HDF5[END_REF] allow to keep a high degree of semantics within files, while leveraging MPI-IO optimizations. Parallel file systems like GPFS [START_REF] Schmuck | GPFS: A Shared-Disk File System for Large Computing Clusters[END_REF] or PVFS [START_REF] Carns | PVFS: A parallel file system for Linux clusters[END_REF] are also subject to optimization efforts, as they usually represent the main bottleneck of this I/O stack.

As a step forward, we introduce Damaris (Dedicated Adaptable Middleware for Application Resources Inline Steering), an approach targeting large-scale multicore SMP supercomputers. The main idea is to dedicate one or a few cores on each node to I/O and data processing to provide an efficient, scalable-by-design, in-compute-node data processing service. Damaris takes into account user-provided information related to the application, the file system and the intended use of the datasets to better schedule data transfers and processing. It may also respond to visualization tools to allow in-situ visualization without impacting the simulation.

We tested our implementation of Damaris as an I/O backend for the CM1 atmospheric model1 , one of the application intended to run on next generation supercomputer BlueWa-Copyright is held by the author/owner(s). ICS'11, May 31-June 4, 2011, Tuscon, Arizona, USA. ACM 978-1-4503-0102-2/11/05.

ters [START_REF]The Blue Waters Project[END_REF] at NCSA. CM1 is a typical MPI application, originally writing one file per process at each checkpoint using HDF5. Deployed on 1024 cores on BluePrint, the Blue-Water's interim system at NCSA with GPFS as underlying filesystem, this approach induces up to 10 seconds overhead in checkpointing phases every 2 minutes, with a high variability in the time spent by each process to write its data (from 1 to 10 seconds). Using one dedicated I/O core in each 16-cores SMP node, we completely remove this overhead. Moreover, the time spared by the I/O core enables a better compression level, thus reducing both the number of files produced (by a factor of 16) and the total data size. Experiments conducted on the French Grid5000 [3] testbed with PVFS as underlying filesystem and a 24 cores/node cluster emphasized the benefit of our approach, which allows communication and computation to overlap, in a context involving high network contention at multiple levels.

This work was done in the framework of a collaboration between the KerData INRIA -ENS Cachan/Brittany team (Rennes, France) and the NCSA (Urbana-Champaign, USA) within the Joint INRIA-UIUC Laboratory for Petascale Computing.


Hélène Langet
email: helene.langet@ge.com.

Cyril Riddell

Yves Trousset

Elisabeth Lahalle

Arthur Tenenhaus

Gilles Fleury

Nikos Paragios

Sparsity constraints and dedicated acquisition protocols for improved Digital Subtraction Rotational Angiography

Keywords: Digital Subtraction Angiography, Rotational angiography, Compressed sensing, 1-minimization, Iterative reconstruction

Digital Subtraction Rotational Angiography (DSRA) allows reconstruction of three-dimensional vascular structures from two spins: the contrast is acquired after injecting vessels with a contrast medium, whereas the mask is acquired in the absence of injection. The vessels are then detected by subtraction of the mask from the contrast. Standard DSRA protocol samples the same set of equiangular-spaced positions for both spins. Due to technical limitations of C-arm systems, streak artifacts degrade the quality of all three reconstructed volumes.

Recent developments of compressed sensing have demonstrated that it is possible to recover a signal that is sparse in some basis under limited sampling conditions. In this paper, we propose to improve the reconstruction quality of non-sparse volumes when there exists a sparse combination of these volumes. To this purpose, we develop an extension of iterative filtered backprojection that jointly reconstructs the mask and contrast volumes via 1-minimization of sparse priors. A dedicated protocol based upon interleaving both spins is shown to further benefit from the sparsity assumptions, while using the same total number of measurements.

Our approach is evaluated in parallel geometry on simulated phantom data.

I. INTRODUCTION

In interventional radiology, one of the main purposes is the visualization of vascular structures. Since vessels inherit contrast comparable to the one of encompassing tissues, they are imaged by injecting a radio-opaque contrast medium into the blood.

In rotational angiography, a tomographic acquisition (spin) of two-dimensional (2D) X-Ray projection views is used to reconstruct a three-dimensional (3D) model of the injected vessels in their environment. However, it might be difficult to separate the vascular structures from surrounding bones or dense devices such as coils. This can be circumvented by performing two acquisitions in a single protocol, similarly to 2D Digital Subtracted Angiography (DSA) [START_REF] Brody | Digital subtraction angiography[END_REF]. The first acquisition, called mask, is performed without injection, while the second acquisition, called contrast, is performed after vessel opacification. All structures but the vessels are removed by digital subtraction of the mask volume from the contrast volume. One can thus visualize either the vessels alone (subtracted volume), or the vessels and their context (contrast volume) or the context alone (mask volume of bones, tissues and devices).

DSRA is more challenging than DSA, because the injection must be such that all vessels located inside the field of view are fully opacified from the beginning to the end of the contrast spin. The contrast medium is however rapidly flushed in the blood flow. To minimize contrast use, the rotation of the C-arm has to be as fast as possible, while its acquisition frame rate is limited, which significantly restricts the angular sampling. Subsampling has little incidence on visualization of highly contrasted structures, but it generates streak artifacts that hide weakly contrasted structures such as soft tissues. To reconstruct a satisfying subtracted volume, the mask and contrast spins are traditionnally acquired with identical parameters, which allows the straightforward removal of redundant background structures and their associated streaks.

Standard reconstruction of the mask and contrast volumes is obtained independently. Each volume suffers from the same undersampling artifacts and the same noise level. They cannot be combined a posteriori to either reduce noise or artifacts because vessel streaks would then propagate into the mask volume.

In the following, we propose to use recent compressed sensing results to jointly reconstruct the mask and contrast volumes. Based on the assumption that vessels are sparse, we show how improved image quality can be obtained for the non-sparse background structures as well. Through a temporal interpretation of the subtraction problem, we define a coupling between the mask and contrast volumes that allows for their joint reconstruction via imposing sparse priors. We extend the iterative filtered backprojection so that vessel sparsity can be applied and promote redundance of the non-opacified structures captured in both spins. We show that adopting interleaved spins rather than identical spins allows for increasing the angular sampling of the mask and contrast volumes. Two sparsity constraints and two acquisition protocols are combined and evaluated in a numerical simulation phantom study. Potential applications of this work are finally discussed.

II. METHOD Let

f = f C f M be the vector containing the contrast volume (f C ∈ R J ) and the mask volume (f M ∈ R J ),
where J is the number of voxels. Let p = p C p M be the vector containing the injected projections (p C ∈ R I ) and the mask projections (p M ∈ R I ), where I is the total number of measurements. Let R = diag{R C , R M } with R C , R M ∈ R I×J be the block-diagonal matrix describing the trajectory of the contrast and mask spins respectively. The reconstruction problem consists in recovering f given that p = Rf .

A. Least-squares approach

The reconstruction used in practice consists in a discretization of the analytical inverse of the system, i.e. backprojection of the filtered projections: f = R T Dp, where D refers to the ramp filter. This method is known as filtered backprojection (FBP).

There exists an iterative version of FBP (iFBP), which comes round to use the euclidian norm as measure of the distance between p and the filtered projections of f :

f * = argmin f Q(f ) Q(f ) = 1 2 D 1/2 (Rf -p) 2 2 (1)
iFBP scheme is a standard gradient descent with step τ > 0:

f (0) = 0 f (n+1) = f (n) -τ ∇Q(f (n) )
FBP and iFBP provide valid solutions if the system is well-determined. When the system is underdetermined (i.e. I << J), typically when the number of projections is low, these methods set to zero the Fourier components that are unobservable given the measures, which is not justified from a physical point of view, and results in a volume that is deteriorated by streak artifacts.

In a general way, if f is the solution of the fully sampled problem, it can be expressed as:

f = f * + f ⊥ where f * is the iFBP solution f ⊥ ∈ Ker(R)
To reduce the magnitude of f ⊥ , which is the error image made of complementary streak artifacts, it is necessary to constrain f by introducing additional regularization terms to the quadratic functional.

B. Compressed sensing approach

In [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF] [3] notably, it is demonstrated that it is possible to recover a signal from a small number of measurements if it is sparse in some basis, that is, if most of the signal energy is concentrated in a few coefficients only. In practice, minimization of the 1 -norm promotes sparsity [START_REF] Donoho | For most large underdetermined systems of linear equations, the minimal 1 solution is also the sparsest solution[END_REF] and is applied for generating sparse approximations of undersampled signals.

We define a functional J as the sum of a quadratic data fidelity term Q and a sparsity penalty ϕ, i.e J(f ) = Q(f ) + ϕ(f ). The function Q is convex and differentiable, whereas ϕ is convex and non-differentiable due to the singularity of the 1 -norm. This excludes the use of conventional smooth optimization techniques.

Combettes and Pesquet [START_REF] Combettes | Split convex minimization algorithm for signal recovery[END_REF] have shown that it is possible to reach convergence with a simple implementation that consists in an explicit gradient step for minimizing Q (in our case iFBP) and an implicit step applying the constraint ϕ:

f (n+1/2) = f (n) -τ ∇Q(f (n) ) f (n+1) = prox τ,ϕ f (n+1/2) (2)
where the penalization is applied via its proximity operator defined as:

∀x ∈ R J , prox γ,ϕ : x → argmin y∈R J ϕ(y) + 1 2γ x -y 2 2
where γ > 0 and is set equal to the gradient step τ in [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF]. Several accelerated implementations of ( 2) have been proposed. In practice, we used the Fast Iterative Shrinkage Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

C. Temporal interpretation

DSRA can be viewed as a temporal acquisition with two time points t O and t M , such that f C = f (t C ) and f M = f (t M ). It is possible to define an invertible temporal transform H t such that at least one of the components of H t f is sparse.

Here H t is the operator associated with the one-dimensional (1D) Haar wavelet transform:

H t = 1 1 1 -1 Define w by w = w 1 w 2 = H t f = f C + f M f C -f M .
The vector w 2 corresponds to the subtracted volume which contains vascular structures only, and as such, is intrinsically sparse. The use of H t thus allows for isolating the sparse structures induced by contrast injection. 1) Soft-thresholded (ST) reconstruction: 1 penalization through w 2 1 has already been successfully applied to the dynamic reconstruction of naturally sparse volumes such as the coronary arteries [START_REF] Hansis | Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography[END_REF]. In our case of DSRA, we expect that the sparsity constraint on w 2 will promote in f the redundance of the non-opacified structures captured in both spins. We thus minimize the following functional:

J ST (f ) = Q(f ) + λ w 2 1 ( 3
)
where λ is the regularization parameter that allows for controlling the strength of the 1 -penalization. In this case, the proximal operator corresponds to soft-thresholding: S λτ (w 2 ) of threshold λτ where τ is the step defined in [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF]. By minimizing J ST , the maximal sampling that can be achieved for w 1 cannot exceed the addition of the mask and the contrast sampling.

2) Total Variation (TV) reconstruction: The addition of the mask and the contrast sampling may still result in an undersampled problem. Since w 1 is not naturally sparse, an efficient method consists in promoting sparsity of T V (w 1 ) = ∇w 1 1 [START_REF] Candès | Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[END_REF]. The corresponding solution is a piece-wise constant approximation of w 1 . Furthermore, the subtracted volume often contains not only vessels, but also perfused structures such as parenchyma that can be approximated by piecewise constant distributions. Thus, we propose to minimize T V (w 2 ) as well. The corresponding functional is:

J T V (f ) = Q(f ) + λ 1 ∇w 1 1 + λ 2 ∇w 2 1 (4)
To minimize (4), the digital TV filter [START_REF] Chan | The Digital TV Filter and Nonlinear Denoising[END_REF] is used as proximal operator of ∇ • 1 . This filter is applied to w 1 with strength λ 1 τ and to w 2 with strength λ 2 τ . In order to minimize bias in the estimation of w 1 , we set λ 1 << λ 2 .

3) Data acquisition: We consider N equiangular-spaced views of the contrast and mask volumes with two protocols:

• Homogeneous protocol: this is the standard DSRA acquisition where the trajectories of each spin samples the same set of angular positions:

A C = A M = {α n = n N -1 •α range | n = 0 • • • N -1}
• Heterogeneous protocol: alternatively each spin samples two sets of interleaved angular positions defined by:

A M = {α n = n N -1 • α range | n = 0 • • • N -1} A C = {α n = n+δ N -1 • α range | n = 0 • • • N -1}
where δ =]0, 1[ is the angular shift.

Standard DSRA uses redundance of the non-opacified structures in the mask and contrast views to obtain the vessels without background streaks, which is only possible using the homogeneous protocol.

Introduction of sparsity assumptions based on the minimization of J ST or J T V will also remove background streaks in the subtracted volume, whatever the acquisition protocol. However, when interleaving the contrast and mask spins, the angular sampling of the non-opacified structures is increased. In particular, with δ = 0.5, background structures are sampled with 2N equiangular-spaced views, whereas the vessels are sampled with N equiangular-spaced views.

III. RESULTS

Minimization of J ST and J T V was evaluated for each acquisition protocol in parallel geometry on a 2D numerical simulation phantom.

A. Simulations

In order to simulate a realistic anatomy, we used a true abdominal CT cross-section (Fig. 1(a) and Fig. 1(b)) as mask image. We obtained the contrast image by adding to the mask four intense disks (Fig. 1(c)) that represent strongly opacified arteries. Intensities are given in Hounsfield Unit (HU). The value of simulated injected vessels is about 3000 HU, while values of soft tissues are those of the original CT slice, around 0 HU. Acquisitions were simulated with 150 views over 180 • . The angular shift δ = 0.5 was used for the heterogeneous protocol.

FISTA for iFBP was run for 100 iterations with step τ = 0.1. We set λ = 10 in J ST and λ 2 = 10 in J T V . λ 1 in J T V was set to 1, i.e. one tenth of the penalization strength applied on w 2 .

B. Evaluation

Images are visualized with a windowing of window center (WC) 60 HU and window width (WW) 75 HU.

FBP reconstruction from 150 views, which is the current level of image quality in clinical routine, is displayed in Fig. 2. Fig. 2(a) shows details of the vessel streaks that actually degrade the entire field of view, while Fig. 2(b) shows details of the background streaks.

Both J ST and J T V perfectly restored the vessels, since they verify the sparsity constraints, by setting to zero all nonvessel pixels. Consequently mask and contrast reconstructions shared the same background. Artifact level in the background depended on both acquisition protocol and sparsity constraint. This level was quantified by computing over the J b pixels of the background structures the root mean square deviation (RMSD) between the reconstruction and the reference

d = 1 J b J b j=1 f j -fj 2 .
Fig. 3(a) and Fig. 3(b) show that the minimization of J ST from homogeneous spins did not improve the quality of the background over FBP. RMSD values are equivalent in these two cases: d = 18 HU. On the opposite, Fig. 3(c) and Fig. 3(d) show that the heterogeneous protocol allowed for reducing streak artifacts and benefited from a sampling of twice as many views. The RMSD d was reduced to 11 HU.

Fig. 4 shows that the use of TV resulted in less streaks in the homogeneous protocol (Fig. 4(a)) as well as in the heterogeneous protocol (Fig. 4(c)). The RMSD was d = 13 HU in the homogeneous case and was reduced to d = 11 HU in the heterogeneous case. 4(d)) also highlight that the piecewise-constant approximation implied by TV affected the overall texture of the background.

IV. DISCUSSION Soft-thresholded reconstruction did not improve the image quality when spins are homogeneous, even though the simulated vessels verified the sparsity hypothesis. In this latter case, the use of an heterogeneous protocol was required to raise image quality without introducing a bias. However, natural sparsity does not generalize to contrast distributions that may be found in perfused tissues. Reconstruction penalized by TV gave decreased artifact levels, whatever the acquisition protocol, at the expense of a change in the overall texture of the background that may not be clinically acceptable. Still, these results suggest dose reduction strategies. On the one hand, it should be possible to reduce the noise (simultaneously to the reduction of the streaks) and raise image quality, at constant dose level. On the other hand, one could reduce the dose by half without deteriorating image quality.

Confirmation of the applicability of the proposed algorithms in cone-beam geometry and on clinical data is on-going.

Compressed sensing penalties were selected in this work according to the acquisition pattern at hand. An alternative sparsity penalty has been proposed in PICCS [START_REF] Chen | Prior image constrained compressed sensing (PICCS)[END_REF], which relies on a prior image. If the mask image is taken as the prior image, the redundant background structures will be set to the mask and will not be improved further. The sum of the mask and contrast volumes cannot be taken as prior, since it is degraded by vessel streaks. The PICCS approach rather suggests an unbalanced acquisition protocol where many views would be dedicated to sampling the mask, while the contrast volume would be estimated from very few contrast projections and a prior image equal to the mask.

V. CONCLUSION Extension of iterative filtered backprojection with 1minimization has been presented for X-ray Digital Subtraction Rotational Angiography. Evaluation on a simulated phantom in parallel geometry showed that sparsity constraints on vessels allowed to generate streak-free vessel images. More interestingly, the non-sparse non-injected structures were estimated from both mask and contrast acquisition, thus benefiting from an increased angular sampling that reduced the streak level. Overall image quality depended on the selected penalties and acquisition protocols. A simple angular shift between the mask and the contrast acquisition resulted in artifact reduction, whatever the sparsity constraint. This suggests strategies for either reducing the dose or improving the image quality of current clinical X-ray DSRA exams.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Numerical simulation phantom. (a) Mask image; (b) Details of (a); (c) Simulated vessels.

  Fig.3(a) and Fig.3(b)show that the minimization of J ST from homogeneous spins did not improve the quality of the background over FBP. RMSD values are equivalent in these two cases: d = 18 HU. On the opposite, Fig.3(c) and Fig.3(d)show that the heterogeneous protocol allowed for reducing streak artifacts and benefited from a sampling of twice as many views. The RMSD d was reduced to 11 HU.Fig.4shows that the use of TV resulted in less streaks in the homogeneous protocol (Fig.4(a)) as well as in the heterogeneous protocol (Fig.4(c)). The RMSD was d = 13 HU in the homogeneous case and was reduced to d = 11 HU in the heterogeneous case. Figure details (Fig. 4(b) and Fig.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Reconstruction by minimization of J ST from 150 views. WC = 60 HU and WW = 75 HU. (a) Reconstruction of the contrast from homogeneous spins; (b) Details of (a); (c) Reconstruction of the contrast from heterogeneous spins; (d) Details of (c).


Virgile Garcia
email: virgile.garcia@insa-lyon.fr

Nikolai Lebedev
email: lebedev@cpe.fr

Jean-Marie Gorce
email: jean-marie.gorce@insa-lyon.fr

Capacity-Fairness Trade-off Using Coordinated Multi-Cell Processing

Keywords: Multicell processing, Capacity distribution, macro-diversity, interferences coordination, fairness, small cells, Rayleigh fading

Multi-cell processing, also called Coordinated Multiple Point (CoMP), is a promising distributed technique that uses neighbour cells' antennas. It is expected to be the part of next generation cellular standards such as LTE-A. Small cell networks in dense urban environments are limited by interferences and CoMP can strongly take advantage of this fact to improve cell-edge users' throughput. The present study introduces a distributed criterion for mobiles to select their optimal set of Base Stations (BS) to perform CoMP, and evaluates the impact of this association on the fairness and the total cell throughput. For that, we use a known theoretical expression for the capacity outage probability of CoMP under Rayleigh fading and evaluate the goodputs of antennas associations. The proposed criterion is used in combination with α-fair resource allocation to perform a joint double-objective optimization of fairness and efficiency.

I. INTRODUCTION

Today, it is commonly anticipated that cellular networks will further evolve to denser and smaller cells, with a strong spatial frequency reuse, possibly with the whole spectrum resource made available within each cell (reuse 1) [START_REF] Necker | Towards frequency reuse 1 cellular FDM/TDM systems[END_REF]. However, this increasing number of base stations (BS) with shorter coverage range critically increases the inter-cell interference level, especially at cell edges. Next generation (4G) systems intend to use small cells in urban environment to create local coverage to serve only few users at a time.

This dense and uncoordinated cells deployment makes the resource management a hard task, and previously used for 2G/3G systems pattern-based resource allocation schemes turn out to be suboptimal, or too restrictive. More promising approaches for OFDMA systems rely on dynamic fractional frequency reuse [START_REF] Necker | Towards frequency reuse 1 cellular FDM/TDM systems[END_REF], [START_REF] Stolyar | Self-organizing dynamic fractional frequency reuse in ofdma systems[END_REF]. However, these approaches still suffer from different open problems including back-hauling overload, inter cells synchronization and signal coordination. A lot of work is devoted to reuse-1 schemes and to the development of Self-Optimized and Self-Organized Networks (SON) [START_REF]E-UTRAN. Self-configuring and selfoptimizing network (SON) use cases and solutions[END_REF] with minimal human involvement.

As a complementary solution, distributed MIMO in multicell environment appeared recently as a breakthrough technology to improve the energy-capacity trade-off of cellular networks [START_REF] Jing | Multicell downlink capacity with coordinated processing[END_REF]. CoMP was introduced, not for providing the huge increase of the total throughput [START_REF] Toskala | LTE advanced in 3GPP Release 10[END_REF], but rather to improve the SINR of cell edge users through spatial diversity. Edge-users will be able to decode a joint transmission from two or more BS in the same frequency band. That way, a better coverage is ensured with less users experiencing outage, and thus, a more intense reuse of the spectrum in dense areas is allowed. Several cooperation techniques of CoMP have been studied in [START_REF] Bjornson | Distributed multicell and multiantenna precoding: Characterization and performance evaluation[END_REF], [START_REF] Hoydis | On the optimal number of cooperative base stations in network MIMO systems[END_REF]. The maximal sum-rate capacity has been derived in [START_REF] Jing | Multicell downlink capacity with coordinated processing[END_REF] providing an upper-bound and revealing the potential of distributed MIMO.

The major issue related to these approaches is the need for coordination between cells: extensive data exchanges of channel state information (CSI) fed back by the mobiles, algorithms related parameters as well as synchronization data, making it impractical due to backhaul limitations. The evolved MIMO-like coordination, such as distributed beamforming or dirty-paper-coding requires a perfect time synchronisation in order to distributively generate in-phase signals or efficient precoding patterns, which require the knowledge of channel gains [START_REF] Bjornson | Distributed multicell and multiantenna precoding: Characterization and performance evaluation[END_REF]. In [START_REF] Simeone | Downlink multicell processing with limited-backhaul capacity[END_REF], the authors evaluate the impact of a limited feedback. In [START_REF] Papadogiannis | A dynamic clustering approach in wireless networks with Multi-Cell cooperative processing[END_REF], clustering methods are used to pre-define which BSs are working together.

We propose a distributed method with the reduced need of feedback by using less constrained techniques. As described below, the proposed method solely requires average users' SINR computed from powers received from the BS, and fed back by the mobile. The underlying technique is referred to as open-loop MIMO. In an optimal approach, the subband powers are weighted to control the interference level but it turns out that the interference levels cannot be predicted without a strict coordination between BS. This is why we study here a simplified resource allocation mechanism relying on constant transmission powers. Note that Wifi standard follows this principle. At a first glance, this assumption may appear strongly suboptimal, but the main advantage relies on interference stabilisation. On the other hand, adapting a station power to reduce interference to a given user has a limited effect due to the presence of numerous other interfering stations.

If it is easy to find work on multicell processing and clustering [START_REF] Papadogiannis | A dynamic clustering approach in wireless networks with Multi-Cell cooperative processing[END_REF], [START_REF] Zhang | Networked MIMO with clustered linear precoding[END_REF], it is harder to find a practical criterion to decide on cooperation. In most cases, a predefined separation (based on distance or received signal strength) between inner and outer-users is considered, using Reuse1 or CoMP (e.g. [START_REF] Jing | Multicell downlink capacity with coordinated processing[END_REF]). Also, in most of works related to CoMP, the objective is formulated as a sum-rate capacity under power constraints, whereas the fairness is neglected.

The novelties of this paper rely on: 1) A double-objective analysis of CoMP-to our knowledge, the problem of capacity-fairness trade-off has not been studied yet in the context of CoMP algorithms. 2) The proposal of a new distributed criterion for BS association, trading off fairness and throughput. This criterion may be used at mobile (or at its primary BS) to determine whether it requires other BS to cooperate. 3) A model that accounts for uncertainty due to fading, through our recent result on capacity outage distribution in CoMP with Rayleigh fading [START_REF] Garcia | Outage Probability for Multi-Cell Processing under Rayleigh Fading[END_REF].

Section II details assumptions and system modelling. The proposed BS association procedure is presented. The distribution of the capacity under Rayleigh fading [START_REF] Garcia | Outage Probability for Multi-Cell Processing under Rayleigh Fading[END_REF] is also derived to evaluate the resulting users' goodputs. In Section III, the resource allocation algorithm is provided, together with the fairness measure. In Section IV we describe a new cooperation criterion tunable according to different fairness constraints. We further evaluate the benefit of this approach for the efficiencyfairness trade-off, as illustrated in simulation results.

II. SYSTEM MODEL

The present study focuses only on downlink and is based on a dense small cell scenario, which typically corresponds to urban conditions with a high density of base-stations to ensure coverage continuity and user rate requirements. The whole system is considered as a snapshot and the time dimension is therefore omitted. The users are supposed to be randomly distributed on a 2D plan. We aim at optimizing the association of users to base station(s) and the resource allocation to users with different fairness objectives.

A. Assumptions

The proposed approach relies on the following rules and assumptions:

1) Each mobile is firstly associated with one BS referred to as the primary BS. 2) Similar to Wimax and LTE-A standards, a frame-based OFDMA modulation, with a set of time-frequency resource blocks (RB), is considered. 3) Each BS uses the whole RB set at maximal power for its users. Thus, the fairness level is managed by the policy defining the blocks allocated to each user. 4) Each RB experiences an independent Rayleigh fading. 5) Due to the previous assumptions, the interference level at a mobile user is almost stable within the block (averaged over the fading). 6) Each mobile feeds back a vector of average powers received from BSs in its vicinity. 7) We assume that any serving BS can request a neighbour BS to help serving one of its users; typically, this will be the strongest interfering BS for this user. 8) Since the exact and instanteneous CSI is not known, and no clock synchronisation between the BS assumed, the beamforming is not possible, and the Alamouti like precoding is used to maximize the mobile received SINR. 9) We assume single antenna receiver.

B. Multi-cell coordination principle

The BSs distributed coordination for the CoMP algorithm holds in four steps:

1) Power sensing: Each mobile measures the average power received from each neighbour BS during the beacon period. It returns a powers vector

P [P 1 , • • • , P N ]
to its primary BS. 2) CoMP capacity computation: Based on P , the BS computes the SINR of different CoMP associations, and selects the best policy. This procedure is detailed in next subsection. 3) Multiple BS association: According to the selected CoMP mode, a BS can request help to another BSs to simultaneously serve a given mobile. It is anticipated that even if only a small fraction of users adopts a CoMP approach, it can result in a capacity gain for all users. 4) RB blocks allocation: a negotiation phase is needed to coordinate the RB selection of CoMP users. BSs need to know the transmission rate and the set of cooperating stations used by each user. In this paper, for the sake of consistency, we don't address this algorithmic issue as we rather focus on expected performance. We thus assume that BSs always accept any cooperation request, and automatically adapt its resource allocation so that it respects neighbour's will and total available resources. The exact implementation issue and algorithmic performance will be addressed in a future work.

C. SINR and downlink CoMP

The combination of signals from two or more BSs allows users to improve their SINR level: the total interference power is decreased significantly by using the strongest interferer(s) as useful signal(s). M users are considered in each cell. We note N i the set of BSs that serve user i a coordinated manner, to create a distributed MISO link. At a given mobile i, on a given channel l, the received complex signal is:

y i,l = n∈Ni h i,l,n P i,n x i,l,n + n / ∈Ni h i,l,n P i,n x i,l,n +z (1)
where n is BS index, z is the thermal noise, modeled as normal random variable with a σ z standard deviation h i,l,n ∼ CN (0, 1) are the fast fading random variables with Rayleigh distributed envelope; x i,l,n are the corresponding transmitted symbols; P i,n is the received power, averaged over the fast fading, including pathloss and shadowing effects from BS n to mobile i, P i,n is supposed common to all considered frequencies. The shadowing is considered slowly varying with respect to the period of measurement reports from mobiles, with constant gain over the RB. We then suppose that P i,n is known at BS. Assuming a perfect CSI at the receiver side and coherent detection, the received SINR is:

γ i,l (N i ) = n∈Ni |h i,l,n | 2 P i,n n / ∈Ni |h i,l,n | 2 P i,n + σ 2 z , (2)
This SINR formulation results from the sum of random phase signals. In addition to the contribution of a (strong) interferer, CoMP offers a macro-diversity due to distant locations of the BSs. One of the objectives of this publication is to determine the set N i that optimizes the system performances (spectral efficiency and fairness).

D. Capacity outage probability

To simplify notations, we drop the i and l indexes when no confusion is possible. According to Theorem 1 in [START_REF] Garcia | Outage Probability for Multi-Cell Processing under Rayleigh Fading[END_REF], the probability for a channel of a given user with a set N of cooperating stations to be in capacity outage (i.e. that fails to reach a target capacity R th = log 2 (1 + γ th )) is:

P out N (R th ) = P log 2 (1 + γ(N )) ≤ R th = (3) 1 - n∈N e -γ th σ 2 z Pn j∈N ,j =n P n P n -P j k / ∈N P n P k γ th + P n , (4)

E. Parallel channels distribution

LTE uses OFDMA, with sub-carriers grouped into resources blocks (RB) which form a time-frequency allocation slot. RB are assumed to experience a flat-fading across all sub-carriers they contain. Uncorrelated fading can be assumed between the RBs, which therefore can be considered as parallel channels, when allocated to the same user. Capacity outage probability of a given user with L parallel channels [START_REF] Tse | Fundamentals of wireless communication[END_REF] is:

P out N ,L (R th ) = P 1 L L l=1 log 2 (1 + γ l ) < R th (5)
= • • • R L-1 F N L • R th - L-1 l=1 x l L-1 l=1 f N (x l )dx l (6)
Where f N (x l ) are the probability density functions (PDF) of the capacity, obtained by derivation of (4).

Considering a system without perfect CSIT, the best way to allocate L resource blocks to a user is to spread the RBs all over the available spectrum (to obtain frequency diversity) with an equal power.

F. Goodput

We use the goodput as capacity metric C (bps/Hz). The goodput is defined as the user's spectral efficiency resulting from successful transmissions. For a given user i, once the set N i is chosen and L i RBs allocated, his goodput G i can be maximized with respect to the transmission rate R.

G i = max R R 1 -P out Ni,Li (R) (7)
Independently of N i , G i grows with L i , since parallel channels provide a frequency diversity and mitigate fading. But the allocation of the number of channels for each user depends on the required fairness among users and their respective SINR. This is introduced in Section III.

An optimization of the set N i is proposed in Section IV taking into account the cost in resource of multiple BS transmission and depending on fairness criteria.

III. CAPACITY-FAIRNESS TRADE-OFF

To quantify the trade-off between fairness and efficiency, it appears necessary to introduce a double-objective characterisation of this system. As pointed out, the goodput of the user i is the function of the amount of resources L i allocated to it, and of channels' SINRs γ l (via the capacity outage probability). Under the obvious constraint of limited total quantity of resources, this allocation need to be done fairly among users.

We choose to represent the Jain's index of fairness [START_REF] Jain | A quantitative measure of fairness and discrimination for resource allocation in shared computer systems[END_REF] and the cell total throughput (bps/Hz) as the two measures to balance on a per cell basis. A clear distinction has to be made between the α-fairness allocation, which is the criterion for the resources sharing and the Jain's index, which is a measure of how fair this allocation is.

A. Jain's index fairness

The Jain's index J ∈ [0, 1], is defined as:

J = M i=1 L i G i 2 M M i=1 (L i G i ) 2 . ( 8
)
This index is advantageous by its fixed bounds and intuitive numerical interpretation (higher index means more fair): J = 1/M for the least fair allocation, meaning only the best user is served by the station. J = 1 for the fairest case (i.e. all the users have the same capacity). Fairness can be computed either locally for each cell, or in a global way for the system.

B. α-fairness

This study uses the well known α-fairness [START_REF] Mo | Fair end-to-end window-based congestion control[END_REF] policy defining the resources partitioning between the users. α = 0 offers the maximum throughput, α = 1 is the proportional fairness, α = 2 is the harmonic fairness, α → ∞ is the maximum fairness (or "uniform capacity").

To be "α-fair", an allocation L = [L 1 , ..., L M ] of resources among M users has to satisfy:

arg max L i 1 1-α (L i C i ) 1-α for α ≥ 0, α = 1 arg max L i log(L i C i ) for α = 1 (9)
In the following, we use the α-fair allocation of resources among users, since it allows a practical and easily parametrisable allocation method (through the adaptation of α), straightforward for the desired tradeoff.

IV. COMP AND FAIRNESS

A. Cooperation criterion

Mobiles are aware of their neighbour cells mean powers, so a practical algorithm based on these simple measurements should enable the system to select the cooperating BSs directly by the mobile or by the attached BS in a fully distributed manner. An appropriate criterion has to be defined to make this selection.

To fairly compare the performance of mobiles attached to one BS with those attached to several BSs, we need to take into account the multiplication of the resource cost induced the latter: we need to look at users' per-BS spectral efficiency S i (bps/Hz/BS):

S i = G i card(N i ) (10)
To determine the set of cooperating stations N i , every user then just need to maximize his G i , with respect to N i . This intuitive utility function was used for example in [START_REF] Zhang | Networked MIMO with clustered linear precoding[END_REF], [START_REF] Garcia | Multi-cell processing for uniform capacity improvement in full spectral reuse system[END_REF]. Figure 1 is an illustration of the trade-off that has to be made, using simulation environment described in subsection IV-B. Each point of the curves is taken for a particular value of α and its corresponding allocation L, from which the total throughput is deduced. The arrows represent how the Jain-Throughput couples evolve using CoMP, for given values of α with criterion (10) applied at every mobile. This criterion, referred to as simple CoMP, can be seen advantageous when more fairness is required, compared to no CoMP case.

However, this somewhat arbitrary criterion used for uniform capacity (i.e., when the objective is to have the same throughput for every user) may not be optimal for other fairness policies: the performance in both efficiency and fairness depends on the BS association. One can see for instance that the gain in both fairness and efficiency is not the same for all fairness policies and using this simple coordination criterion is not worthy when looking at low fairness (α < 1). This phenomenon is explained by the fact that low fairness tends to serve outer users with more resources than inner users. This decreases the cell efficiency for an increase of fairness.

We therefore propose in this paper another, more general criterion for users' association to BSs, in order to balance the impact of CoMP on the goodput and the number of cooperating stations to use. The proposed criterion has been motivated so that a simple parameter β provides a large class of functions, including concave (0 ≤ β < 1) and convex (β > 1) ones :

N i = arg max N (G i ) β card(N ) , (11)
Therefore, when β is high, more users use cooperative antennas: a slight increase of the spectral efficiency using more BSs will easily compensate the cost in resources induced. On the contrary, with β → 0, CoMP is never used, since the criterion judges any improvement of capacity as not significant enough. Each different value of β has its own impact on the network and associates differently BSs to users, and thus modifies the resource allocation result, regardless of the value of α. This general formulation provides the criterion which can be adjusted jointly with the resources allocation to optimize the trade-off between capacity and fairness.

B. Simulation results

Simulation parameters can be found in Table I. A twodimensional area with 120 BSs has been generated, but only 60 cells are taken into account to avoid border effects-the interference being the most important limiting effect, it is important not to limit the study to a very simplistic model (for example considering only first-tier interferers). The cells' 1: Jain's index versus total throughput. locations are randomly generated around the regular hexagonal grid. The users are uniformly distributed within each cell.

We are looking at the evolution of the couple Efficiency-Fairness compared to Reuse1 mode. The objective is obviously to improve either the Jain's index or the total throughput, or even both when possible.

In the global system allocation (BS and resources), the parameters α and β are not interdependent, as they are used to tune two separate parts of the optimization procedure. Nevertheless, as shown previously, they both have an impact on the performance and so it is interesting to investigate their joint effect. The objective is then to obtain the best possible couples (α,β) for the system.

Consider the Figure 2. The purple circles pattern shows the impact of β for only one value of α = 2. The three curves represent no-CoMP (lower continuous blue line), simple CoMP (10) (black dashed line) and CoMP case with β criterion (red dotted line) obtained as the envelope corresponding to the best couples (α, β), of the cloud of patterns representing all possible fairness-efficiency couples, obtained by extensive simulations. (In other words, the optimal couples are such combinations (α, β) that no other couple can improve both the fairness and total throughput).

We can see that the variation in both fairness and cell efficiency is hardly tractable, since it depends on the specific SINRs of users in each cell. However, it can be seen that when β = 0, no user will ever choose multiple antennas, so that we will have the same result as for Reuse1. For growing values of β, CoMP first improves both fairness and efficiency before reaching a limit and starts decreasing the efficiency, while still improving the fairness. When β is too high, the performance of the system is even worse than without any BS cooperation, and we would rather use no coordination with a different α value instead. Carrying out simulations for other values of α, revealed the importance of the selected cooperation criterion value β, compared to the simple coordination criterion. When the fairness is important, the envelope tends to be similar to [START_REF] Zhang | Networked MIMO with clustered linear precoding[END_REF], but when lower fairness is required, an important gain in both metric appears, and keeps the CoMP more efficient and fair than without coordination. Some other interesting questions arise. For example, can we describe among all the optimal couples, the interaction between the values of α and β. Specifically, one could ask whether for every α-fair allocation, the corresponding optimal point is delivered by the same unique value of β ? There is no simple parametrization rule and both α and β may have several optimal points, again, depending on users' SINRs. We noted that the optimal value of β tends to 1 when extreme fairness is required, but decreases with less fairness, tending to 0 for maximum throughput. This can be explained by the fact that low fairness criterion tends to give the most of resources to users with higher SINRs. Reallocating a part of good users' bandwidth to worse users strongly decreases their capacity, while only slightly increasing the low-rate users' throughput. The total throughput of the cell is then obviously reduced at a price of only slight increase of fairness. V. CONCLUSION This paper introduces a generalized distributed criterion for mobiles to select their optimal set of Base Stations (BS) to perform the CoMP. We compared its performance with a simple CoMP scheme and no cooperation cases, and shown its benefit for most fairness-efficiency scenarios. We showed that our approach can improve both the fairness and the total throughput, especially for low SINR users which otherwise can consume too many resources.

Fig. 2 :

 2 Fig. 2: Jain's index vs total throughput. Blue line for no cooperation, black dashed line is the simple CoMP, and the red dotted curve is the envelope of the best (α, β) couples. Purple points correspond to the α = 2, for 0 ≤ β ≤ 5 case.

TABLE I :

 I Simulation parameters.

ACKNOWLEDGEMENT

work has been carried out in the frame of the joint lab between INRIA and Alcatel-Lucent Bell Labs on "Self Organizing Networks".