Pascal Bégout
email: begout@ann.jussieu.fr

Convergence to Scattering States in the Nonlinear Schr ödinger Equation

In this paper, we consider global solutions of the following nonlinear Schrödinger equation iut

|x| 2 ; dx). We show that, under suitable conditions, if the solution u satisfies e -it∆ u(t) -u± → 0 in X as t → ±∞ then u(t)e it∆ u± → 0 in X as t → ±∞. We also study the converse. Finally, we estimate | u(t) Xe it∆ u± X | under some less restrictive assumptions.

Introduction and notations

We consider the following Cauchy problem,

     i ∂u ∂t + ∆u + λ|u| α u = 0, (t, x) ∈ (-T * , T * ) × R N , u(0) = ϕ, in R N , (1.1)
where λ ∈ R, 0 α < 4 N -2 (0 α < ∞ if N = 1) and ϕ a given initial data.

It is well-known that if λ < 0, α > 4 N and ϕ ∈ H 1 (R N ), then there exists u ± ∈ H 1 (R N ) such that lim t→±∞ T (-t)u(t)-u ± H 1 = 0 (Ginibre and Velo [START_REF] Ginibre | Scattering theory in the energy space for a class of nonlinear Schrödinger equations[END_REF], Nakanishi [START_REF] Nakanishi | Energy scattering for nonlinear Klein-Gordon and Schrödinger equations in spatial dimensions 1 and 2[END_REF][START_REF] Nakanishi | Remarks on the energy scattering for nonlinear Klein-Gordon and Schrödinger equations[END_REF]). Since (e it∆ ) t∈R is an isometry on H 1 (R N ), we also have lim

t→±∞ u(t) -T (t)u ± H 1 = 0. Furthermore, if α > -(N -2)+ √ N 2 +12N +4 2N
and if ϕ ∈ X ≡ H 1 (R N ) ∩ L 2 (|x| 2 ; dx), then there exist u ± ∈ X such that lim t→±∞ T (-t)u(t) -u ± X = 0 (Tsutsumi [START_REF] Tsutsumi | Scattering problem for nonlinear Schrödinger equations[END_REF]). The same result holds without assumption on the λ's sign if the initial data is small 2000 Mathematics Subject Classification: 35Q55 (35B40, 35B45, 35P25) 1 enough in X and if α > 4 N + 2 (Cazenave and Weissler [START_REF] Cazenave | Rapidly decaying solutions of the nonlinear Schrödinger equation[END_REF]). Note that to have these limits, we have to make a necessary assumption on α (Barab [START_REF] Barab | Nonexistence of asymptotically free solutions for a nonlinear Schrödinger equation[END_REF], Strauss [START_REF] Strauss | Nonlinear scattering theory[END_REF][START_REF] Strauss | Nonlinear scattering theory at low energy[END_REF], Tsutsumi and Yajima [START_REF] Tsutsumi | The asymptotic behavior of nonlinear Schrödinger equations[END_REF]), 2

N < α < 4 N -2 (2 < α < ∞ if N = 1).
The purpose of this paper is to study the asymptotic behavior of u(t) -T (t)u ± X under the assumption lim t→±∞ T (-t)u(t) -u ± X = 0, and the converse. In the linear case (i.e. : λ = 0) or if the initial data is 0, the answer is trivial since T (-t)u(t) -u ± ≡ u(t) -T (t)u ± ≡ 0, for all t ∈ R. Since (e it∆ ) t∈R is an isometry on H 1 (R N ), the equivalence on H 1 (R N ) is trivial. But (e it∆ ) t∈R is not an isometry on X and so it is natural to wonder whether or not we have lim t→±∞ u(t) -T (t)u ± X = 0 when lim t→±∞ T (-t)u(t) -u ± X = 0 and conversely.

This paper is organized as follows. In Section 2, we give the main results. In Section 3, we establish some a priori estimates. In Section 4, we prove Theorems 2.1, 2.4, 2.5 and Proposition 2.8. In Section 5, we prove Theorem 2.10.

Before closing this section, we give some notations which will be used throughout this paper and we recall some properties of the solutions of the nonlinear Schrödinger equation.

z is the conjugate of the complex number z; Re z and Im z are respectively the real and imaginary part of the complex number z; ∆ = N j=1 ∂ 2 ∂x 2 j ; for 1 p ∞, p ′ is the conjugate of the real number p defined by 1 p + 1 p ′ = 1 and

L p = L p (R N ) = L p (R N ; C) with norm . L p ; H 1 = H 1 (R N ) = H 1 (R N ; C) with norm . H 1 ; for all (f, g) ∈ L 2 × L 2 , (f, g) = Re R N f (x)g(x)dx; X = ψ ∈ H 1 (R N ; C); ψ X < ∞ with norm ψ 2 X = ψ 2 H 1 (R N ) + R N |x| 2 |ψ(x)| 2 dx; (T (t)
) t∈R is the group of isometries (e it∆ ) t∈R generated by i∆ on L 2 (R N ; C); C are auxiliary positive constants and C(a 1 , a 2 , . . . , a n ) indicates that the constant C depends only on parameters a 1 , a 2 , . . . , a n and that the dependence is continuous.

It is well-known that for every ϕ ∈ X, (1.1) has a unique solution u ∈ C((-T * , T * ); X) which satisfies the conservation of charge and energy, that is for all t ∈ (-T * , T * ), u(t [START_REF] Cazenave | An introduction to nonlinear Schrödinger equations, volume 26 of Textos de Métodos Matemáticos[END_REF], Ginibre and Velo [START_REF] Ginibre | On a class of nonlinear Schrödinger equations. III. Special theories in dimensions 1, 2 and 3[END_REF][START_REF] Ginibre | On a class of nonlinear Schrödinger equations. I. The Cauchy problem, general case[END_REF][START_REF] Ginibre | On a class of nonlinear Schrödinger equations. II. Scattering theory, general case[END_REF][START_REF] Ginibre | The global Cauchy problem for the nonlinear Schrödinger equation revisited[END_REF], Kato [START_REF] Kato | On nonlinear Schrödinger equations[END_REF]). Definition 1.1. We say that (q, r) is an admissible pair if the following holds.

) L 2 = ϕ L 2 and E(u(t)) = E(ϕ) def = 1 2 ∇ϕ 2 L 2 -λ α+2 ϕ α+2 L α+2 . Moreover, if λ 0, if α < 4 N or if ϕ H 1 is small enough then T * = T * = ∞ and u L ∞ (R;H 1 ) < ∞ (see for example Cazenave
(i) 2 r 2N N -2 (2 r < ∞ if N = 2, 2 r ∞ if N = 1), (ii) 2 q = N 1 2 -1 r .
Note that in this case 2 q ∞ and q = 4r N (r -2) .

Definition 1.2. We say that a solution u ∈ C((-T * , T * ); X) of (1.1) has a scattering state u + at

+∞ (respectively u -at -∞) if T * = ∞ and if u + ∈ X is such that lim t→∞ T (-t)u(t) -u + X = 0 (respectively if T * = ∞ and if u -∈ X is such that lim t→-∞ T (-t)u(t) -u -X = 0).
We recall the Strichartz' estimates. Let I ⊆ R, be an interval, let t 0 ∈ I, let (q, r) and (γ, ρ)

be two admissible pairs, let ϕ ∈ L 2 (R N ) and let f ∈ L γ ′ (I; L ρ ′ (R N )). Then the following integral equation defined for all t ∈ I, u(t) = T (t)ϕ + i t t0
T (t -s)f (s)ds, satisfies the following inequality

u L q (I,L r ) C 0 ϕ L 2 +C 1 f L γ ′ (I;L ρ ′ ) , where C 0 = C 0 (N, r) and C 1 = C 1 (N, r, ρ).
For more details, see Keel and Tao [START_REF] Keel | Endpoint Strichartz estimates[END_REF].

The main results

Theorem 2.1.

Let λ = 0, 2 N < α < 4 N -2 (2 < α < ∞ if N = 1)
, ϕ ∈ X and let u be the solution of (1.1) such that u(0) = ϕ. We assume that u has a scattering state u ± at ±∞ (see Definition 1.2).

Then the following holds.

1. (a) If N 2 and if α > 4 N then lim t→±∞ u(t) -T (t)u ± X = 0. (b) If 3 N 5 and if α > 8 N + 2 then lim t→±∞ u(t) -T (t)u ± X = 0.
2. If N = 1 and α = 4 or if 3 N 5 and α = 8 N + 2 then we have,

sup t 0 u(t) -T (t)u + X < ∞ and sup t 0 u(t) -T (t)u -X < ∞.
Remark 2.2. Remark that in Theorem 2.1, no hypothesis on the λ ′ s sign is made.

Remark 2.3. N ∈ {3, 4, 5} =⇒ 4 N < 8 N + 2 < 6 N < 4 N -2 . N 6 =⇒ 4 N -2 8 N + 2
.

Despite the fact we do not know if lim

t→±∞ u(t) -T (t)u ± X = 0 when α 8 N + 2 (α 4/N if N 2)
or when N 6, we can give an estimate of the difference of the norms, as shows the following theorem, without any restriction on the dimension space N and on α (except α > 2 N ). Since under the scattering state assumption we always have lim

t→±∞ u(t) -T (t)u ± H 1 = 0, it is sufficient to estimate | xu(t) L 2 -xT (t)u ± L 2 | as t -→ ±∞. N < α < 4 N -2 (2 < α < ∞ if N = 1)
, ϕ ∈ X and let u be the associated solution of (1.1). Assume that u has a scattering state u ± at ±∞ (see Definition 1.2). We define for

all t ∈ R, A ± (t) = xu(t) L 2 -xT (t)u ± L 2 and h(t) = xu(t) 2 L 2 .
Then we have the following result.

sup t 0 | u(t) X -T (t)u + X | < ∞ and sup t 0 | u(t) X -T (t)u -X | < ∞,
with the following estimates.

1. If α < 4 N then - C ∇u ± L 2 lim inf t→±∞ A ± (t) lim sup t→±∞ A ± (t) ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 . 2. If α > 4 N then ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 lim inf t→±∞ A ± (t) lim sup t→±∞ A ± (t) C ∇u ± L 2 . 3. If α = 4 N then lim t→±∞ A ± (t) = ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 . Furthermore, h ′ (0) = 4Im R N ϕ(x)x.∇ϕ(x)dx and C = C(sup t∈R T (-t)u(t) X , N, α, λ). Theorem 2.5. Let λ > 0, 2 N < α < 4 N -2 (2 < α < ∞ if N = 1)
, ϕ ∈ X and let u be the associated solution of (1.1). Assume that u has a scattering state u ± at ±∞ (see Definition 1.2). We define for

all t ∈ R, A ± (t) = xu(t) L 2 -xT (t)u ± L 2 and h(t) = xu(t) 2 L 2 .
Then we have the following result.

sup t 0 | u(t) X -T (t)u + X | < ∞ and sup t 0 | u(t) X -T (t)u -X | < ∞,
with the following estimates.

1. If α < 4 N then ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 lim inf t→±∞ A ± (t) lim sup t→±∞ A ± (t) C ∇u ± L 2 . 2. If α > 4 N then - C ∇u ± L 2 lim inf t→±∞ A ± (t) lim sup t→±∞ A ± (t) ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 . 3. If α = 4 N then lim t→±∞ A ± (t) = ± h ′ (0) + 4(xu ± , i∇u ± ) 4 ∇u ± L 2 . Furthermore, h ′ (0) = 4Im R N ϕ(x)x.∇ϕ(x)dx and C = C(sup t∈R T (-t)u(t) X , N, α, λ).
Remark 2.6. By Theorem 2.1 and 2 of Theorems 2.4 and 2.5

, if α > 4 N when N 2 or if α > 8 N + 2 when N ∈ {3, 4, 5}, we have Im R N u -(x)x.∇u -(x)dx Im R N ϕ(x)x.∇ϕ(x)dx Im R N u + (x)x.∇u + (x)dx, if λ < 0, Im R N u + (x)x.∇u + (x)dx Im R N ϕ(x)x.∇ϕ(x)dx Im R N u -(x)x.∇u -(x)dx, if λ > 0. Remark 2.7. If ϕ ∈ H 1 (R N ; C) satisfies ϕ ≡ aψ with ψ ∈ H 1 (R N ; R) and a ∈ C, then we have h ′ (0) ≡ d dt xu(t) 2 L 2 |t=0 ≡ 4Im R N ϕ(x)x.∇ϕ(x)dx = 0.
The following proposition offers others estimates.

Proposition 2.8. Let λ = 0, 2 N < α < 4 N -2 (2 < α < ∞ if N = 1)
, ϕ ∈ X and let u be the solution of (1.1) such that u(0) = ϕ. Assume that u has a scattering state u ± at ±∞. Then the following estimates hold.

1. If λ < 0, lim sup t→±∞ ( xu(t) L 2 -xT (t)u ± L 2 ) xu ± L 2 ∇u ± L 2 ± (xu ± , i∇u ± ) ∇u ± L 2 . 2. If λ > 0, lim inf t→±∞ ( xu(t) L 2 -xT (t)u ± L 2 ) - xu ± L 2 ∇u ± L 2 ∓ (xu ± , i∇u + ) ∇u ± L 2 .
Remark 2.9. By 3 of Theorems 2.4 and 2.5 and by Proposition 2.8, if α = 4 N then we have,

-xu -L 2 ∇u -L 2 1 4 d dt xu(t) 2 L 2 |t=0 xu + L 2 ∇u + L 2 , if λ < 0, -xu + L 2 ∇u + L 2 1 4 d dt xu(t) 2 L 2 |t=0 xu -L 2 ∇u -L 2 , if λ > 0.
Now we give the result concerning the converse.

Theorem 2.10.

Let λ = 0, 2 N < α < 4 N -2 (2 < α < ∞ if N = 1)
, ϕ ∈ X and u be the associated solution of (1.1). Assume that u is global in time and there exists u + ∈ X and u -∈ X such that

lim t→±∞ u(t) -T (t)u ± X = 0. Let α 0 = -(N -2) + √ N 2 + 12N + 4 2N
. Then, we have the following result.

1. If λ < 0 and if α α 0 (α > α 0 if N = 2) then lim t→±∞ T (-t)u(t) -u ± X = 0. 2. If λ > 0 and if α > 4 N then lim t→±∞ T (-t)u(t) -u ± X = 0. 3. If α > 4 N + 2 and if ϕ X is small enough then lim t→±∞ T (-t)u(t) -u ± X = 0.
Remark 2.11. Note that in the case 2, no hypothesis on the ϕ X ' s size is made.

Remark 2.12. Assume there exists u

± , v ± ∈ X such that lim t→±∞ T (-t)u(t) -u ± X = 0 and lim t→±∞ u(t) -T (t)v ± X = 0. Then we have, u + = v + and u -= v -. Indeed, since X ֒→ L 2 (R N ) and
T (t) is an isometry on L 2 (R N ), we have lim t→±∞ T (-t)u(t) -u ± L 2 = lim t→±∞ T (-t)u(t) -v ± L 2 = 0.
Hence the result.

Remark 2.13.

α 0 ∈ 4 N +2 , 4 N (α 0 ∈ 2 N , 4 N if N = 1).

A priori estimates

Throughout this section, we make the following assumptions.

   λ = 0, 2 N < α < 4 N -2 (2 < α < ∞) if N = 1), ϕ ∈ X, u ∈ C(R; X)
is the associated solution of (1.1) and has a scattering state u ± ∈ X at ± ∞.

(3.1)

We define the following real. T (-t)u(t) X < ∞). Let (q, r) be an admissible pair (see Definition 1.1). Then the following holds.

γ * =        α -2 2 , if N = 1, α(N + 2) -4 4 , if N 2. ( 3
1. For all t = 0, u(t

) L r C|t| -2 q , where C = C(sup t∈R T (-t)u(t) X , N, r). 2. If furthermore α > 4 N + 2 then u ∈ L q (R; W 1,r (R N )).
Proof. We follow the method of Cazenave [START_REF] Cazenave | An introduction to nonlinear Schrödinger equations, volume 26 of Textos de Métodos Matemáticos[END_REF] (see Theorem 7.2.1 and Corollary 7.2.4). We set w(t, x) = e -i |x| 2 4t u(t, x) and f (u) = λ|u| α u. We already know that for every admissible pair (q, r), u ∈ L q loc (R; W 1,r (R N )) (see for example Cazenave [START_REF] Cazenave | An introduction to nonlinear Schrödinger equations, volume 26 of Textos de Métodos Matemáticos[END_REF]; Theorem 5.3.1 and Remark 5.3). We only prove the case t > 0, the case t < 0 following by applying the result for t > 0 to u(-t) solution of (1.1) with initial value ϕ. We proceed in 2 steps.

Step

1. u(t) L r C(sup t∈R T (-t)u(t) X , N, r)|t| -2
q , for every admissible pair (q, r) and for all t = 0.

We have xT

(-t)u(t) L 2 = (x+2it∇)u(t) L 2 C. Furthermore, (x+2it∇)u(t, x) = 2ite -i |x| 2 4t ∇w(t, x).
Using the Gagliardo-Nirenberg's inequality, we obtain

u(t) L r ≡ w(t) L r C ∇w(t) N ( 1 2 -1 r ) L 2 w(t) 1-N ( 1 2 -1 r ) L 2 C (x + 2it∇)u(t) L 2 |t| -1 N ( 1 2 -1 r ) C|t| -N ( 1 2 -1 r ) .
Hence the result.

Step 2. u ∈ L q (R; W 1,r (R N )) for every admissible pair (q, r).

By the Strichartz' estimates and by Hölder's inequality (applying twice), we have

f (u) L q ′ ((0,∞);W 1,r ′ ) C u L qα q-2 ((0,1);L rα r-2 ) + u L qα q-2 ((1,∞);L rα r-2 ) α u L q ((0,∞);W 1,r ) , (3.3) u L q ((S,∞);W 1,r ) C + C u α L qα q-2 ((S,∞);L rα r-2 ) u L q ((S,∞);W 1,r ) , (3.4)
for all S > 0 and for every admissible pair (q, r).

Case N 3. We set r = 4N 2N -α(N -2) . Since α ∈ 0, 4 N -2 , we have r ∈ 2, 2N N -2 .
So we can take q such that (q, r) is an admissible pair. For this choice of r, we have rα r-2 = 2N N -2 and q q-2 = 4 4-α(N -2) . By (3.4) and the first step we have for all S > 0,

u L q ((S,∞);W 1,r ) C + C ∞ S t -4α 4-α(N -2) dt q-2 q u L q ((S,∞);W 1,r ) . And 4α 4 -α(N -2) > 1 ⇐⇒ α > 4 N + 2
. Thus, there exists S 0 > 0 large enough such that

C ∞ S0 t -4α 4-α(N -2) dt q q-2 1 2 ,
and then,

u L q ((S0,∞);W 1,r ) 2C.
For this choice of (q, r), we deduce from (3.3) that f (u) L q ′ ((0,∞);W 1,r ′ ) < ∞. Hence the result for every admissible pair by the Strichartz' estimates.

Case N=2. Since α > 1 is fixed, we take r > 2 sufficiently close to 2 to have α >

2(r -1) r . So, in particular, rα r -2 > 2. Moreover, q q -2 = r 2
where q is such that (q, r) is an admissible pair. So by Hölder's inequality (twice), (3.4) and the first step, we have for all S > 0,

u L q ((S,∞);W 1,r ) C + C ∞ S t -rα-2(r-2) 2 dt 2 r u L q ((S,∞);W 1,r ) . And rα -2(r -2) 2 > 1 ⇐⇒ α > 2(r -1) r .
And we conclude exactly as the case N 3.

Case N=1. We take (3.3) with the admissible pair (∞, 2) and apply the first step. So,

f (u) L 1 ((0,∞);H 1 ) C u L α ((0,1);L ∞ ) + u L α ((1,∞);L ∞ ) α u L ∞ ((0,∞);H 1 ) C + C u α L α ((1,∞),L ∞ ) C + C ∞ 1 t -α 2 dt < ∞.
Hence the result for every admissible pair by the Strichartz' estimates.

Remark 3.2. We set v(t, x) = (x + 2it∇)u(t, x). In the same way as the above proof, we can show under the assumptions of Proposition 3.

1 that if α > 4 N + 2
, then for every admissible pair (q, r), . Then, the following estimates hold.

v ∈ L q (R; L r (R N )) (
1. If N = 1 then for all t = 0 we have T (-t)u(t) -u ± H 1 C|t| -γ * .

2. If N = 2 then for all t = 0 and for any γ < γ * , we have T (-t)u(t) -u ± H 1 C|t| -γ .

3. If N 3 then for all t = 0 we have

T (-t)u(t) -u ± H 1 C|t| -γ * .
Proof. Denote f (u) = λ|u| α u. We only prove the case t > 0, the case t < 0 following by applying the result for t > 0 to v(t) = u(-t) solution of (1.1) with v(0) = ϕ. In this case, v + = u -and the result follows. By applying the Strichartz' estimates and Hölder's inequality (twice), we have

u(t) -T (t)u + H 1 C f (u) L q ′ (t,∞);W 1,r ′ ) C u α L qα q-2 ((t,∞);L rα r-2 ) u L q (R;W 1,r ) ,
for every admissible pair (q, r) and for all t > 0. Thus, by Proposition 3.1, 2, we have

u(t) -T (t)u + H 1 C u α L qα q-2 ((t,∞);L rα r-2 )
, for every admissible pair (q, r) and for all t > 0. Now, we conclude by the same way than for the step 2 of the proof of Proposition 3.1, using 1 of this proposition, considering separately the three cases N = 1, N = 2, N 3, and using the same admissible pairs. This achieves the proof.

u(t) -T (t)u + X < ∞ and sup t 0 u(t) -T (t)u -X < ∞ when N = 2 and α = 4 N .
We define the following function h by

∀t ∈ (-T * , T * ), h(t) = xu(t) 2 L 2 . (4.1)
Lemma 4.2. Let u satisfying (3.1) and let h be defined by (4.1). Then h ∈ C 2 (R), and we have

∀t ∈ R,      h ′ (t) = 4Im R N u(t, x)x.∇u(t, x)dx, h ′′ (t) = 2N α ∇u + 2 L 2 -2(N α -4) ∇u(t) 2 L 2 . Furthermore, if λ < 0, ∀t ∈ R, ∇u(t) L 2 ∇u + L 2 , (4.2)
and if λ > 0,

∀t ∈ R, ∇u(t) L 2 ∇u + L 2 . (4.3)
Proof. See Ginibre and Velo [START_REF] Ginibre | On a class of nonlinear Schrödinger equations. II. Scattering theory, general case[END_REF] or Cazenave [START_REF] Cazenave | An introduction to nonlinear Schrödinger equations, volume 26 of Textos de Métodos Matemáticos[END_REF], Proposition 6.4.2 to have h ∈ C 2 (R), the expression of h ′ and ∀t ∈ R, h ′′ (t) = 4N αE(ϕ) -2(N α -4) ∇u(t) 2 L 2 . Furthermore, using the conservation of energy and lim t→±∞ u(t) L α+2 = 0 (Proposition 3.1), we obtain ∇u + 2 L 2 = 2E(ϕ). Which gives, with the above identity, the desired expression of h ′′ . Finally, with the equality ∇u + 2 L 2 = 2E(ϕ) and the conservation of energy, we obtain (4.2) if λ < 0 and (4.3) if λ > 0. Now, we establish 2 lemmas which will be used to prove Theorem 2.4. Lemma 4.3. Let u satisfying (3.1) and let h be defined by (4.1). Assume that λ < 0. Then the following holds.

(i) If α < 4 N then lim sup t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) h ′ (0) + 4(xu + , i∇u + ) 4 ∇u + L 2 . (ii) If α > 4 N then lim inf t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) h ′ (0) + 4(xu + , i∇u + ) 4 ∇u + L 2 . (iii) If α = 4 N then lim t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) = h ′ (0) + 4(xu + , i∇u + ) 4 ∇u + L 2 .
Proof. We proceed in 4 steps.

Step

1. (a) If α < 4 N then ∀t 0, h ′ (t) h ′ (0) + 2N α ∇u + 2 L 2 t. (b) If α > 4 N then ∀t 0, h ′ (t) h ′ (0) + 8 ∇u + 2 L 2 t. (c) If α = 4 N then ∀t ∈ R, h ′ (t) = h ′ (0) + 8 ∇u + 2 L 2 t.
We integrate between 0 and t 0 the function h ′′ of Lemma 4.2.

α 4 N =⇒ -2(N α -4) 0 =⇒ (a). α 4 N =⇒ -2(N α -4) 0 with (4.2) =⇒ (b). α = 4 N =⇒ -2(N α -4) = 0 =⇒(c). Step 2. (a) If α < 4 N then ∀t 0, h ′ (t) h ′ (0) + 8 ∇u + 2 L 2 t. (b) If α > 4 N then ∀t 0, h ′ (t) h ′ (0) + 2N α ∇u + 2 L 2 t.
We integrate between 0 and t 0 the function h ′′ of Lemma 4.2.

α < 4 N =⇒ -2(N α -4) > 0 and (4.2) =⇒(a). α > 4 N =⇒(b).
Step 3. (a) If α < 4 N then ∀t 0,

xϕ 2 L 2 + h ′ (0)t + N α ∇u + 2 L 2 t 2 h(t) xϕ 2 L 2 + h ′ (0)t + 4 ∇u + 2 L 2 t 2 . (b) If α > 4 N then ∀t 0, xϕ 2 L 2 + h ′ (0)t + 4 ∇u + 2 L 2 t 2 h(t) xϕ 2 L 2 + h ′ (0)t + N α ∇u + 2 L 2 t 2 . (c) If α = 4 N then ∀t ∈ R, h(t) = xϕ 2 L 2 + h ′ (0)t + 4 ∇u + 2 L 2 t 2 .
It is sufficient to integrate between 0 and t 0 the formulas of steps 1 and 2 to obtain the step 3.

Step 4. Conclusion.

We set :

g(t) = xϕ 2 L 2 + h ′ (0)t + 4 ∇u + 2 L 2 t 2
, t > 0 large enough. Then for t > 0 large enough, we have the following asymptotic development:

g(t) = 2 ∇u + L 2 t + h ′ (0) 4 ∇u + L 2 + xϕ 2 L 2 4 ∇u + L 2 t + o 1 t . (4.4)
In the same way, for t > 0 large enough, we have :

xT (t)u + L 2 = 2 ∇u + L 2 t - (xu + , i∇u + ) ∇u + L 2 + xu + 2 L 2 4 ∇u + L 2 t + o 1 t . (4.5)
And, applying the step 3 (a), (4.4) and (4.5) and taking lim sup t→∞ , we get (i). Indeed,

xu(t) L 2 -xT (t)u + L 2 2 ∇u + L 2 t + h ′ (0) 4 ∇u + L 2 + xϕ 2 L 2 4 ∇u + L 2 t -2 ∇u + L 2 t + (xu + , i∇u + ) ∇u + L 2 - xu + 2 L 2 4 ∇u + L 2 t + o 1 t = h ′ (0) 4 ∇u + L 2 + xϕ 2 L 2 4 ∇u + L 2 t + (xu + , i∇u + ) ∇u + L 2 - xu + 2 L 2 4 ∇u + L 2 t + o 1 t .
Hence (i) by taking lim sup t→∞ in the above expression.

By applying the step 3 (b), (4.4) and (4.5) and taking lim inf t→∞ , we get (ii) by the same way.

By applying the step 3 (c), (4.4) and (4.5) and letting t -→ ∞, we get (iii) by the same way.

Lemma 4.4. Let u satisfying (3.1). Assume that λ < 0. Then the following holds.

(i) If α < 4 N then lim inf t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) - C ∇u ± L 2 . (ii) If α > 4 N then lim sup t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) C ∇u ± L 2 . Furthermore, C = C(sup t∈R T (-t)u(t) X , N, α, λ).
Proof. We proceed in 2 steps. Let h be defined by (4.1).

Step 1. ∀t 1, -

t 1 s 1 ∇u(σ) 2 L 2 dσds C + Ct -N α-4 2 + Ct + ∇u + 2 t - 1 2 ∇u + 2 L 2 t 2 , where C = C(sup t∈R T (-t)u(t) X , N, α, λ). By Proposition 3.1, we have u(σ) α+2 L α+2 C(sup t∈R T (-t)u(t) X , N, α)σ -N α 2
, for all σ > 0. With the conservation of energy and the formula ∇u + 2 L 2 = 2E(ϕ), we deduce that for all σ > 0,

∇u + 2 L 2 -∇u(σ) 2 L 2 C(sup t∈R T (-t)u(t) X , N, α, λ)σ -N α 2 (since λ < 0). Integrating this expression over [1, t] × [1, s],
we obtain the desired result.

Step 2. Conclusion.

(i) Lemma 4.2 implies that for every t 0,

h(t) = xϕ 2 L 2 + h ′ (0)t + N α ∇u + 2 L 2 t 2 -2(N α -4) t 0 s 0 ∇u(σ) 2 L 2 dσds. (4.6)
By (4.6) and step 1, we obtain,

∀t 1, h(t) C + h ′ (0)t + N α ∇u + 2 L 2 t 2 + Ct 4-N α 2 -Ct + 2(N α -4) ∇u + 2 L 2 t -(N α -4) ∇u + 2 t 2 .
And so for all t > 1, h(t) C + Ct

4-N α 2 + (h ′ (0) -2(4 -N α) ∇u + 2 L 2 -C)t + 4 ∇u + 2 L 2 t 2
. By an asymptotic development on this last inequality, we obtain

xu(t) L 2 2 ∇u + L 2 t + Ct -1 + Ct -N α-2 2 + h ′ (0) -2(4 -N α) ∇u + 2 L 2 -C 4 ∇u + L 2 + o 1 t ,
for all t > 0 large enough. From this last expression and (4.5), we obtain (i) (see the end of the proof of (i) of Lemma 4.3).

(ii) From (4.6) and step 1, we obtain, for all t > 1,

h(t) C + f ′ (0)t + N α ∇u + 2 L 2 t 2 + Ct 4-N α 2 + Ct + 2(N α -4) ∇u + 2 L 2 t -(N α -4) ∇u + 2 L 2 t 2 .
And so, for all t > 1, h(t) C + Ct

4-N α 2 + (h ′ (0) + 2(N α -4) ∇u + 2 L 2 + C)t + 4 ∇u + 2 L 2 t 2
. By an asymptotic development on this last inequality, we obtain Proof. We proceed as for the proof of Lemma 4.3, using (4.3) instead of (4.2).

xu(t) L 2 2 ∇u + L 2 t + Ct -1 + Ct -N α-2 2 + h ′ (0) -2(4 -N α) ∇u + 2 L 2 + C 4 ∇u + L 2 + o 1 t
Lemma 4.6. Let u satisfying (3.1). Assume that λ > 0. Then the following holds.

(i) If α < 4 N then lim sup t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) C ∇u ± L 2 . (ii) If α > 4 N then lim inf t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) - C ∇u ± L 2 .
Furthermore, C = C(sup t∈R T (-t)u(t) X , N, α, λ).

Proof. We proceed in 2 steps.

Step 1. ∀t 1, -

t 1 s 1 ∇u(σ) 2 L 2 dσds -C -Ct -N α-4 2 -Ct + ∇u + 2 L 2 t - 1 2 ∇u + 2 L 2 t 2 , where C = C(sup t∈R T (-t)u(t) X , N, α, λ). By Proposition 3.1, we have u(σ) α+2 L α+2 C(sup t∈R T (-t)u(t) X , N, α)σ -N α 2
, for all σ > 0. With the conservation of energy and the equality ∇u + 2 L 2 = 2E(ϕ), we deduce that for all σ > 0,

∇u + 2 L 2 -∇u(σ) 2 L 2 -C(sup t∈R T (-t)u(t) X , N, α, λ)σ -N α 2 (since λ > 0). Integrating this expression over [1, t] × [1, s],
we get the desired result.

Step 2. Conclusion.

First, note that we have for all t > 1,

t 0 s 0 ∇u(σ) 2 L 2 dσds = t 0 1 0 ∇u(σ) 2 L 2 dσds + 1 0 s 1 ∇u(σ) 2 L 2 dσds + t 1 s 1 ∇u(σ) 2 L 2 dσds t 0 1 0 ∇u(σ) 2 L 2 dσds + t 1 s 1 ∇u(σ) 2 L 2 dσds.
And so, for all t ∈ R and for every λ = 0.

t 0 s 0 ∇u(σ) 2 L 2 dσds C( ϕ H 1 , N, α, λ)t + t 1 s 1 ∇u(σ)
We have also the following estimate.

xT (t)u + L 2 = 2t ∇u + L 2 - (xu + , i∇u + ) ∇u + L 2 + xu + 2 L 2 4 ∇u + L 2 t + o 1 t , (4.8)
for all t > 0 large enough and for every λ = 0.

We first establish 1 in the case t > 0 that we note in this proof 1 + . 1 in the case t < 0 is obviously noted 1 -. By (4.7), (4.8) and (4.2), we have

xu(t) L 2 -xT (t)u + L 2 2t ∇u(t) L 2 + xT (-t)u(t) L 2 -2t ∇u + L 2 + (xu + , i∇u + ) ∇u + L 2 - xu + 2 L 2 4 ∇u + L 2 t + o 1 t xT (-t)u(t) L 2 + (xu + , i∇u + ) ∇u + L 2 - xu + 2 L 2 4 ∇u + L 2 t + o 1 t
for all t > 0 large enough.

Hence 1 + by taking lim sup t→∞ in the above expression. 1 -follows by applying 1 + to v solution of (1.1) such that v(0) = ϕ. Indeed, by uniqueness v(t) = u(-t) and so v + = u -. Then, using 1 + and the identity T (t)ψ = T (-t)ψ which holds for all t ∈ R and ψ ∈ L 2 , we obtain 1 -.

Using (4.3) instead of (4.2), we obtain 2 by the same way.

5 Proof of Theorem 2.10

Proof of Theorem 2.10. We proceed in 2 steps.

Step 1. We have 1 and 3.

It is well-known that there exist v ± ∈ X such that T (-t)u(t) X ----→ t→±∞ v ± (Cazenave and Weissler [START_REF] Cazenave | Rapidly decaying solutions of the nonlinear Schrödinger equation[END_REF]).

The result follows from Remark 2.12.

Step 2. We have 2.

We set v(t, x) = (x + 2it∇)u(t, x), w(t, x) = e -i |x| 2 4t u(t, x) and f (u) = λ|u| α u. Since T (t) is an isometry on H 1 (R N ), we only have to show that xT (-t)u(t) v ∈ L q (R; L r (R N )), (5.2) for every admissible pair (q, r) (see for example Cazenave [START_REF] Cazenave | An introduction to nonlinear Schrödinger equations, volume 26 of Textos de Métodos Matemáticos[END_REF], Theorem 7.5.3 for (5.1); following the proof of this theorem with v instead of u, yields (5.2)).

From (5.2) we have in particular, v ∈ L ∞ (R; L 2 (R N )) and so by Proposition 3.1, 1,

u(t) L r C|t| -2 q , (5.3)
for every admissible pair (q, r) and for all t = 0.

We have the following integral equation. For all t ∈ R, u(t) = T (t)u ± -i We also have (x + 2it∇)u(t, x) = 2ite i |x| 2 4t ∇w(t, x). Moreover, e -i |x| 2 4t f (u(t, x)) = f (w(t, x)). Thus, (x + 2it∇)f (u(t, x)) = 2ite i , for every admissible pair (q, r) and for all t > 0.

We use this inequality with the admissible pair (q, r) such that r = α + 2 and we apply (5.3). Then, xT (-t)u(t) -xu + L 2 C ∞ t s -2α q-2 ds q-2 q

Ct -2α-(q-2) q t→∞ ---→ 0, since α > 4 N =⇒ 2α > q -2. Hence the result. The case t < 0 follows with the same method.

Proposition 3 . 3 .

 33 follow the step 2 of the proof of Corollary 7.2.4 of Cazenave [2]; consider separately the three cases N = 1, N = 2, N 3 and replace the admissible pairs therein by those of the proof of Proposition 3.1). Let γ * be defined by (3.2). Assume u satisfies (3.1) and α > 4 N + 2

L 2 (

 2 R N ) -----→ t→±∞ xu ± . We have X ֒→ L α+2 (R N ), thus u(t) -T (t)u ± L α+2 (R N ) ------→ t→±∞ 0 and so lim t→±∞ u(t) L α+2 = lim t→±∞ T (t)u ± L α+2 = 0. Therefore, since α > 4 N , u ∈ L q (R; W 1,r (R N )),(5.1)

T

  (t -s)f (u(s))ds, from which we deduce,∀t ∈ R, T (t)(xT (-t)u(t) -xu ± ) = -i ±∞ t T (t -s)(x + 2is∇)f (u(s))ds.(5.4)

2 C

 2 |x| 2 4t ∇f (w(t, x)), and so|(x + 2it∇)f (u(t, x))| = 2|t||∇f (w(t, x))|. Finally, (x+2it∇)f (u) L r ′ = 2|t| ∇f (w) L r ′ C|t| |w| α ∇w L r ′ .From this inequality, by using the Hölder's inequality twice, we deduce that (note that |w| = |u|)(x + 2it∇)f (u) L q ′ (I,L r ′ ) 2it∇)u L q (I,L r ) ,(5.5)for every admissible pair (q, r) and for any interval I ⊆ R. From (5.4), from the Strichartz' estimates, from (5.5) and from (5.2), we havexT (-t)u(t) -xu + L 2 = T (t)(xT (-t)u(t) -xu + ) L (x + 2is∇)f (u) L q ′ ((t,∞);L r ′ )

  Remark 4.1. Since we do not have the estimate 2 of Proposition 3.3 for γ = γ

	t < 0. So to conclude, it is sufficient to prove that lim t→∞	xu(t) -xT (t)u + L 2 = 0 to obtain 1 and
	sup t 0	xu(t) -xT (t)u + L 2 < ∞ to obtain 2. We have
		xu(t) -xT (t)u + = xu(t) -T (t)xu + + 2itT (t)∇u +
		= xu(t) + 2it∇u(t) -T (t)xu + + 2itT (t)∇u + -2it∇u(t)
	4 Proof of Theorems 2.1, 2.4, 2.5 and Proposition 2.8
	Throughout this section, we assume that u satisfies (3.1).
	Proof of Theorem 2.1. Since lim

t→±∞

T (-t)u(t) -u ± X = 0 and T (t) is an isometry on H 1 , we have lim t→±∞ u(t) -T (t)u ± H 1 = 0. Thus, it is sufficient to prove that lim

t→±∞ xu(t) -xT (t)u ± L 2 = 0

to obtain 1 and that sup

t 0 xu(t) -xT (t)u + L 2 < ∞ and sup t 0 xu(t) -xT (t)u -L 2 < ∞ to obtain 2.

Suppose that the result is proved for t > 0. Then we apply it to v(t) = u(-t) solution of

(1.1)

with initial data v(0) = ϕ. Then v + = u -is the scattering state at +∞ of u(-t). And using the identity T (t)ψ = T (-t)ψ which holds for all t ∈ R and for every ψ ∈ L 2 , we obtain the result for

= T (t) [(xT (-t)u(t) -xu + ) + 2it(∇u + -T (-t)∇u(t))] ,

for all t > 0. With Proposition 3.3, we obtain

xu(t) -xT (t)u + L 2 xT (-t)u(t) -xu + L 2 + 2t T (-t)∇u(t) -∇u + L 2 xT (-t)u(t) -xu + L 2 + Ct -(γ-

1)

, for all t > 0 and for all γ ∈ (0, γ * ] if N = 2 and for all γ ∈ (0, γ * ) if N = 2, where γ * is defined by (

3.2)

. And by assumption, lim t→∞ xT (-t)u(t) -xu + L 2 = 0 and γ * -1 > 0 if and only if α > 8 N +2 if N 3 and α > 4 N if N 2. Thus, in the above expression, it is sufficient to choose γ = γ * if N = 2, and γ ∈ (0, γ * ) close enough to γ * if N = 2. Hence the result. * , we do not know whether or not sup t 0

  Proof of Theorem 2.4. As for the proof of Theorem 2.1, it is sufficient to show that | xu(t) L 2 -xT (t)u + L 2 | remains bounded as t -→ ∞. Lemmas 4.3 and 4.4 achieve the proof and give the desired estimates.

	Now, we establish 2 lemmas which will be used to prove Theorem 2.5. The proof is very similar
	to the Theorem 2.4 one.
	Lemma 4.5. Let u satisfying (3.1) and let h be defined by (4.1). Assume that λ > 0. Then the
	following holds.
	(i) If α < (ii) If α > (iii) If α =	4 N 4 N 4 N	then lim inf t→∞ then lim sup ( xu(t) L 2 -xT (t)u + L 2 ) t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) then lim t→∞ ( xu(t) L 2 -xT (t)u + L 2 ) =	h ′ (0) + 4(xu + , i∇u + ) 4 ∇u + L 2 h ′ (0) + 4(xu + , i∇u + ) . 4 ∇u + L 2 h ′ (0) + 4(xu + , i∇u + ) 4 ∇u + L 2 .	.
					,
	for all t > 0 large enough. With this last expression and (4.5), we obtain (ii) (see the end of the proof
	of (i) of Lemma 4.3).
	Now, we are able to prove Theorem 2.4.

  2 L 2 dσds, for all t > 1. And we proceed as for the proof of Lemma 4.4.

	Proof of Proposition 2.8. By Cauchy-Schwarz' inequality, we have
	| xu(t) L 2 -2t ∇u(t) L 2 |	xT (-t)u(t) L 2 ,	(4.7)

Now, we are able to prove Theorem 2.5.

Proof of Theorem 2.5. As for the proof of Theorem 2.1, it is sufficient to show that

| xu(t) L 2 -xT (t)u + L 2 | remains bounded as t -→ ∞.

Lemmas 4.5 and 4.6 achieve the proof and give the desired estimates.

Acknowledgments

I would like to express my great gratitude to Professor Thierry Cazenave, my thesis adviser, for having suggested me this work and for encouragement and helpful advices.


Fabienne Comte

Charles-André Cuenod

Marianna Pensky

Yves Rozenholc

Laplace deconvolution and its application to Dynamic Contrast Enhanced imaging

Keywords: AMS 2010 subject classifications. 62G05, 62G20, 62P10 Laplace deconvolution, complexity penalty, Dynamic Contrast Enhanced imaging

In the present paper we consider the problem of Laplace deconvolution with noisy discrete observations. The study is motivated by Dynamic Contrast Enhanced imaging using a bolus of contrast agent, a procedure which allows considerable improvement in evaluating the quality of a vascular network and its permeability and is widely used in medical assessment of brain flows or cancerous tumors. Although the study is motivated by medical imaging application, we obtain a solution of a general problem of Laplace deconvolution based on noisy data which appears in many different contexts. We propose a new method for Laplace deconvolution which is based on expansions of the convolution kernel, the unknown function and the observed signal over Laguerre functions basis. The expansion results in a small system of linear equations with the matrix of the system being triangular and Toeplitz. The number m of the terms in the expansion of the estimator is controlled via complexity penalty. The advantage of this methodology is that it leads to very fast computations, does not require exact knowledge of the kernel and produces no boundary effects due to extension at zero and cut-off at T . The technique leads to an estimator with the risk within a logarithmic factor of m of the oracle risk under no assumptions on the model and within a constant factor of the oracle risk under mild assumptions. The methodology is illustrated by a finite sample simulation study which includes an example of the kernel obtained in the real life DCE experiments. Simulations confirm that the proposed technique is fast, efficient, accurate, usable from a practical point of view and competitive.

Introduction

Cancers and vascular diseases present major public health concerns. Considerable improvement in assessing the quality of a vascular network and its permeability have been achieved through Dynamic Contrast Enhanced (DCE) imaging using a bolus of contrast agent at high frequency such This imaging procedure has great potential for cancer detection and characterization, as well as for monitoring in vivo the effects of treatments. It is also used, for example, after a stroke for prognostic purposes or for occular blood flow evaluation.

As an example, below we consider a DCE-CT experiment that follows the diffusion of a bolus of a contrast agent injected into a vein. At the microscopic level, for a given voxel of interest having unit volume, the number of arriving particles at time t is given by β AIF(t), where the Arterial Input Function (AIF) measures concentration within a unit volume voxel inside the aorta and β is a proportion of the AIF which enters the tissue voxel. Denote the number of particles in the voxel at time t by Y (t) and the random lapse of time during which a particle sojourns in the voxel by S.

Assuming sojourn times for different particles to be i.i.d. with c.d.f. F , one obtains the following equation for the average number of contrast agent particles at the moment t EY (t) = where the expectation is taken under the unknown distribution of the sojourn times. In reality, one does not know EY (t) and has discrete noisy observations

Y (t i ) = EY (t i ) + σ i .
Medical doctors are interested in a reproducible quantification of the blood flow inside the tissue which is characterized by f (t) = β(1 -F (t)) since this quantity is independent of the concentration of particles of contrast agent within a unit volume voxel inside the aorta described by AIF(t).

The sequential imaging acquisition is illustrated by Figure 1. The contrast agent arrives with the oxygenated blood through the aorta (red arrow) where its concentration, AIF, within unit volume voxel is first measured when it passes through the CT cross section (red box). Subsequently, the contrast agent enters the arterial system, and it is assumed that its concentration does not change during this phase. The exchange within the tissue of both oxygen and contrast agent occurs after the arterial phase and the concentration of contrast agent during this exchange is measured in all tissue voxels (grey voxel in the zoom) inside the CT cross section. Later the contrast agent returns to the venous system with the de-oxygenated blood (blue arrow).

To complete description of this experiment, one has to take into account that there is a delay δ between the measurement of the contrast agent concentration inside the aorta (first cross of the CT section) and its arrival inside the tissue. This leads to the following complete model:

Y (t i ) = t i -δ 0 βAIF(t i -τ )(1 -F (τ ))dτ + σ i , i = 1, ..., n.
(

The value of delay δ can be measured with a small error using the decay between the jumps after the injection of the contrast agent inside the aorta and the tissue. Unfortunately, evaluation of the proportion β is a much harder task which is realized with a larger error. In the spirit of complete model (1.1) for DCE-CT experiments, one can consider a more general model of Laplace convolution equation based on noisy observations which presents a necessary theoretical step before obtaining medical answers provided by model (1.1).

Indeed, for a known value of δ, equation (1.1) reduces to a noisy version of a Laplace convolution equation

y(t i ) = t i 0 g(t i -τ )f (τ )dτ + σ i , i = 1, . . . , n, (1.2)
where function g is considered to be known, f is a function of interest, measurements y(t i ) are taken at points 0 ≤ t 1 ≤ ... ≤ t n ≤ T < ∞ and i are i.i.d. N (0, 1). The corresponding noiseless version of this equation can be written as

q(t) = t 0 g(t -τ )f (τ )dτ, t ≥ 0. (1.3)
Formally, by setting g(t) = f (t) ≡ 0 for t < 0, equation (1.3) can be viewed as a particular case of the Fredholm convolution equation The issue of having measurements only on the part t ≤ T of half line (0, ∞) does not affect the Laplace deconvolution since it exhibits causality property: the values of q(t) for 0 ≤ t ≤ T depend on values of f (t) for 0 ≤ t ≤ T only and vice versa.

q(t) = b a g(t -τ )f (τ )dτ, (1.4
y(t i ) = b a g(t i -τ )f (τ )dτ + σ i , i = 1, ..., n, (1.5
The mathematical theory of (noiseless) convolution type Volterra equations is well developed (see, e.g., [START_REF] Gripenberg | Volterra Integral and Functional Equations[END_REF]) and the exact solution of (1.3) can be obtained through Laplace transform. However, direct application of Laplace transform for discrete measurements faces serious conceptual and numerical problems. The inverse Laplace transform is usually found by application of tables of inverse Laplace transforms, partial fraction decomposition or series expansion (see, e.g., [START_REF] Polyanin | Handbook of Integral Equations[END_REF], neither of which is applicable in the case of the discrete noisy version of Laplace deconvolution. Only few applied mathematicians and researchers in natural sciences took an effort to solve the problem using discrete measurements in the left hand side of (1.4).

Since the problem arises in medical imaging, few scientists put an effort to solve equation (1.1) using singular value decomposition (SVD) with the subsequent application of Tikhonov regularization (see, e.g., [START_REF] Axel | Cerebral blood flow determination by rapid-sequence computed tomography: theoretical analysis[END_REF], [START_REF] Ostergaard | High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis[END_REF] and an extensive review in Fieselmann et al.

( 3) satisfies a particular linear differential equation and, hence, can be recovered using q(t) and its derivative q (t). For this particular kind of kernel, the authors derived convergence rates for the quadratic risk of the proposed estimators, as n increases, under

the assumption that the s-th derivative of f is continuous on (0, ∞). However, in Dey, Martin and Ruymgaart (1998) it is assumed that data are available on the whole positive half-line (i.e. T = ∞)

and that s is known (i.e., the estimator is not adaptive).

Recently, Abramovich et al. (2012) studied the problem of Laplace deconvolution based on discrete noisy data. The idea of the method is to reduce the problem to estimation of the unknown regression function, its derivatives and, possibly, some linear functionals of these derivatives. The estimation is carried out using kernel method with the Lepskii technique for the choice of the bandwidth (although it is mentioned in the paper that other methodologies for the choice of bandwidth can also be applied). The method has an advantage of reducing a new statistical problem to a well studied one. However, the shortcoming of the technique is that it requires meticulous boundary correction and is strongly dependent on the knowledge of the kernel g. Indeed, small change in the kernel may produce significant changes in the expression for the estimator.

In the present paper we suggest a method which is designed to overcome limitations of the previously developed techniques. The new methodology is based on expansions of the kernel, unknown function f and the right-hand side in equation (1.2) over the Laguerre functions basis.

The expansion results in a small system of linear equations with the matrix of the system being triangular and Toeplitz. The number of the terms in the expansion of the estimator is controlled via complexity penalty. The advantage of this methodology is that it leads to very fast computations and produces no boundary effects due to extension at zero and cut-off at T . The technique does not require exact knowledge of the kernel since it is represented by its Laguerre coefficients only and leads to an estimator with the risk within a logarithmic factor of m of the oracle risk under no assumptions on the model and within a constant factor of the oracle risk under mild assumptions.

Another merit of the new methodology includes the fact that, since the unknown functions are represented by a small number of Laguerre coefficients, it is easy to cluster or classify them for various groups of patients. Simulation study shows that the method is very accurate and stable and easily outperforms SVD and kernel-based technique of Abramovich, Pensky, and Rozenholc (2012).

The rest of the paper is organized as follows. In Section 2 we derive the system of equations resulting from expansion of the functions over the Laguerre basis, study the effect of discrete, possible irregularly spaced data and introduce selection of model size via penalization. Corollary 1 indeed confirms that the risk of the penalized estimator lies within a logarithmic factor of m of the minimal risk. In Section 3 we obtain asymptotic upper bounds for the risk and prove the risk lies within a constant factor of an oracle risk. The proof of this fact rests on nontrivial facts of the theory of Toeplitz matrices. Section 4 provides a finite sample simulation studies. Finally, Section 5 discusses results obtained in the paper. Section 6 contains proofs of the results in the earlier sections.

2 Laplace deconvolution via expansion over Laguerre functions basis

Relations between coefficients of the Laguerre expansion

One of the possible solution of the problem (1.2) is to use Galerkin method with the basis represented by a system of Laguerre functions. Laguerre functions are defined as

φ k (t) = √ 2ae -at L k (2at), k = 0, 1, . . . , (2.1)
where L k (t) are Laguerre polynomials (see, e.g., [START_REF] Gradshtein | Tables of integrals, series, and products[END_REF])

L k (t) = k j=0 (-1) j k j t j j! , t ≥ 0.
It is known that functions φ k (•), k = 0, 1, . . ., form an orthonormal basis of the L 2 (0, ∞) space and, therefore, functions f (•), g(•), q(•) and y(•) can be expanded over this basis with coefficients f (k) , g (k) , q (k) and y (k) , k = 0, . . . , ∞, respectively. By plugging these expansions into formula (1.3), we obtain the following equation

∞ k=0 q (k) φ k (t) = ∞ k=0 ∞ j=0 f (k) g (j) t 0 φ k (x)φ j (t -x)dx. (2.
2)

It turns out that coefficients of interest f (k) , k = 0, 1, . . . , can be represented as a solution of an infinite triangular system of linear equations. Indeed, it is easy to check that (see, e.g., 7.411.4

in Gradshtein and Ryzhik (1980))

t 0 φ k (x)φ j (t -x)dx = 2ae -at t 0 L k (2at)L j (2a(t -x))dx = (2a) -1/2 [φ k+j (t) -φ k+j+1 (t)].
Hence, equation (2.2) can be re-written as

∞ k=0 q (k) φ k (t) = ∞ k=0 φ k (t)[(2a) -1/2 f (k) g (0) + k-1 l=0 (2a) -1/2 (g (k-1) -g (k-l-1) )f (l) ].
Equating coefficients for each basis function, we obtain an infinite triangular system of linear equations. In order to use this system for estimating f , we define

f m (x) = m-1 k=0 f (k) φ k (x), (2.3)
approximation of f based on the first m Laguerre functions. The following Lemma states how the coefficients in (2.3) can be recovered.

Lemma 1. Let f m , g m and q m be m-dimensional vectors with elements f (k) , g (k) and q (k) , k = 0, 1, . . . , m -1, respectively. Then, for any m, one has q m = G m f m where G m is the lower triangular Toeplitz matrix with elements

G (i,j) =        (2a) -1/2 g (0) , if i = j, ( 2a
) -1/2 (g (i-j) -g (i-j-1) ), if j < i, 0, if j > i.
(2.4) Hence, f (x) can be estimated by

fm (x) = m-1 k=0 f (k) φ k (x) (2.5)
where f m = G -1 m qm and qm is an unbiased estimator of the unknown vector of coefficients q m .

Recovering Laguerre coefficients from discrete noisy data

Unfortunately, unlike some other linear ill-posed problems, data does not come in the form of unbiased estimators of the unknown coefficients q (k) , k = 0, 1, . . . , m -1. Below, we examine how the length of the observation interval T and spacing of observations t i in equation (1.2) affect the system of equations in Lemma 1.

Let P : [0, T ] → [0, T ] be a function generating observations in (1.2) such that P is a continuously differentiable strictly increasing function

P (0) = 0, P (T ) = T, P (t i ) = iT /n, i = 1, . . . , n. (2.6)
Under conditions (2.6), P is a one-to-one function and, therefore, has an inverse P -1 .

Choose M large enough that the bias in representation (2.3) of f by f M is very small and form

an (n × M ) matrix Φ M with elements Φ (i,k) = φ k (t i ), i = 1, . . . , n, k = 0, . . . , M -1. Let z M be
the M -dimensional vector with elements z (j) = y, φ j , j = 0, . . . , M -1. Then, it follows that

y(t i ) = M -1 l=0 z (l) φ l (t i ) = (Φ M z M ) (i) , i = 1, . . . , n.
If y and h are n-dimensional vectors with components y(t i ) and q(t i ), i = 1, . . . , n, respectively, then the vectors q M and z M of the true and the estimated Laguerre coefficients of q(x) can be represented, respectively, as

q M = (Φ T M Φ M ) -1 Φ T M h, z M = (Φ T M Φ M ) -1 Φ T M y. (2.7)
Let us examine matrix Φ T M Φ M . Note that, for any k and l,

(Φ T M Φ M ) (k,l) = n i=1 φ k (P -1 (iT /n))φ l (P -1 (iT /n)) ≈ nT -1 T 0 φ k (P -1 (x))φ l (P -1 (x))dx = nT -1 T 0 φ k (t)φ l (t)p(t)dt
where p(t) = P (t). It follows from the above that matrix Φ T M Φ M should be normalized by a factor n -1 T . Indeed, if points t i are equispaced on the interval (0, T ], then, for n and T large enough,

(Φ T M Φ M ) ≈ nT -1 I M where I M is the M -dimensional identity matrix.
Hence, in what follows, we are going to operate with matrix A M = T n -1 (Φ T M Φ M ) and its inverse

Ω M = (A M ) -1 = nT -1 (Φ T M Φ M ) -1 . (2.8)
Let be the vector with components (t i ), i = 1, . . . , n, and

ξ M = n/T (Φ T M Φ M ) -1 Φ T M .
Then, the vector f M of the true Laguerre coefficients of the unknown function f satisfies the following equation

z M = G M f M + σ T /n ξ M , ξ M ∼ N (0, Ω M ).
(2.9)

If points t i are equispaced on the interval (0, T ] and both n and T are large, then, in (2.9), Ω M ≈ I M .

Model selection and oracle risk

Equation (2.9) implies that one can estimate unknown vector f M by f M = G -1 M z M . However, since the value of M is large, the variance of this estimator,

E f M -G -1 M q M 2 = σ 2 T n Tr(G -1 M Ω M G -T M ),
where A -T = (A T ) -1 = (A -1 ) T for any matrix A, will be too large while the bias in the representation (2.3) of f by f M will be very small. Hence, in order to balance the bias and the variance components of the error, one needs to choose the best possible number m of Laguerre functions in the representation (2.3) of f , i.e., choose the model size.

In order to achieve a required balance between the bias and the variance components of the error, consider a collection of integer indices M n = {1, . . . , M } where M < n may depend on n and, for m ∈ M n , the associated subspaces S m ⊆ R M defined by t ∈ S m if t = (t (0) , t (1) , . . . , t (m-1) , 0, 0, . . . , 0) T .

Let us denote by z m , q m , f m , ξ m and f m the M -dimensional vectors where the first m elements coincide with the elements of m-dimensional vectors z m , q m , f m , ξ m and f m respectively, and the last (M -m) elements are identical zeros. For each m ∈ M n , evaluate

f m = (G m ) -1 z m = G -1 M z m
and denote

Q m = G -1 m Ω m G -T m .
(2.10)

For the estimator fm of f given by (2.5) with the vector of coefficients f m , the bias-variance decomposition of the mean squared error is of the form

E fm -f 2 = f m -f 2 + σ 2 T n -1 Tr(Q m ), (2.11)
where the bias term

f m -f 2 = ∞ j=m (f (j)
) 2 is decreasing and the variance term σ 2 T n -1 Tr(Q m ) is growing with m. The smallest possible risk, the so-called oracle risk, is obtained by minimizing the right-hand side of expression (2.11) with respect to m:

R oracle = min m E fm -f 2 = min m f m -f 2 + σ 2 T n -1 n Tr(Q m ) .
(2.12)

Hence, the objective is to choose a value of m which delivers an estimator of the unknown function f (x) with the risk as close to the oracle risk (2.12) as possible. Since the bias in the right-hand side of expression (2.11) is unknown, in order to attain this goal, one can use a penalized version of estimator (2.5) as it is described in the next section.

Selection of model size via penalization

For any vector t ∈ R M , we define contrast as

γ 2 n (t) = t 2 -2 t, G -1 M z M (2.13)
and note that for t ∈ S m one has, thanks to the nul coordinates of t and the lower triangular form

of G M and G m , t, G -1 M z M = t, G -1 M z m = t, f m .
Let A 2 = Tr(A T A) and A = λ max (A T A) be, respectively, the Frobenius and the spectral norm of a matrix A, where λ max (U ) is the largest eigenvalue of U . Denote

v 2 m = Q m 2 2 = Tr(Q m ), ρ 2 m = Q m 2 = λ max (Q m ) (2.14)
where Q m is a lower triangular matrix such that (

Q m ) T Q m = Q m .
Assume that ρ 2 m grows at most polynomially in m, i.e. there exist positive constants α and C ρ such that

ρ 2 m ≤ C ρ m α .
(2.15)

Choose any constant B > 0 and introduce a penalty

pen(m) = 4σ 2 T n -1 (1 + B)v 2 m + (1 + B -1 )(2α + 2)ρ 2 m log m . (2.16)
For each m = 1, . . . , M , construct estimator fm (x) of f (x) of the form (2.5) with coefficients

f m = (G m ) -1 z m , the augmented version f m satisfies f m = arg min t∈Sm γ 2 n (t). Now, choose m = m where m = arg min m ∈ M n : γ 2 n ( f m ) + pen(m) .
(2.17)

The following statement holds.

Theorem 1. Let condition (2.15) hold for some positive constants α and C ρ . Then, for any B > 0,

one has R( f m ) := E( f m -f 2 ) ≤ min m∈Mn 3 f m -f 2 + 4pen(m) + 16C 2 ρ σ 2 (1 + B -1 ) T mn . (2.18)
The proof of this and later statements are given in Section 6.

Note that the upper bound in Theorem 1 is non-asymptotic and holds for any values of T and n and any distribution of points t i , i = 1, . . . , n.

In order to evaluate relative precision of the estimator f m constructed above, we shall compare its risk with the oracle risk (2.12). Since ρ 2 m ≤ v 2 m for any value of m, it follows from Theorem 1 that, for any value of m, the risk of the estimator f m lies within a logarithmic factor of the oracle risk, i.e., the estimator is optimal within a logarithmic factor of m. In particular, the following corollary holds.

Corollary 1. Under conditions of Theorem 1, R( f m ) ≤ 16[(1 + B) + (1 + B -1 )(2α + 2) log m 0 ]R oracle + 16C 2 ρ σ 2 (1 + B -1 ) T m -1 0 n -1 , (2.19)
where m 0 = m 0 (n, T ) is the value of m delivering the minimum in the right-hand side of (2.18).

3 Asymptotic upper bounds for the risk and optimality of the estimator

Assumptions

Corollary 1 is valid for any function g and any distribution of sampling points, hence, it is true in the "worst case scenario". In majority of practical situations, however, v 2 m increases much faster with m than ρ 2 m and the risk of the estimator f m can exceed the oracle risk only by a finite factor independent of m 0 and n. In particular, in what follows, we shall show that, under certain conditions, for n large enough and T = T n , the ratio between R( f m ) and R oracle is bounded by a constant independent of n.

For this purpose, assume that function g(x), its Laplace transform G(s) and matrix Ω m defined in (2.8) satisfy the following conditions (A1) There exists an integer r ≥ 1 such that

d j g(t) dt j t=0 = 0, if j = 0, ..., r -2, B r = 0, if j = r -1. (3.1) (A2) g ∈ L 1 [0, ∞) is r times differentiable with g (r) ∈ L 1 [0, ∞).
(A3) Laplace transform G(s) of g has no zeros with nonnegative real parts except for zeros of the form s = ∞ + ib.

(A4) There exists n 0 such that for n > n 0 , eigenvalues of matrix Ω m are uniformly bounded, i.e.

0 < λ 1 ≤ λ min (Ω m ) ≤ λ max (Ω m ) ≤ λ 2 < ∞ (3.2)
for some absolute constants λ 1 and λ 2 .

Introduction to theory of banded Toeplitz matrices

The proof of asymptotic optimality of the estimator f m relies heavily on the theory of banded Toeplitz matrices developed in Böttcher andGrudsky (2000, 2005). In this subsection, we review some of the facts about Toeplitz matrices which we shall use later.

Consider a sequence of numbers

{b k } ∞ k=-∞ such that ∞ k=-∞ |b k | < ∞
(z) = ∞ k=-∞ b k z k , z ∈ C. (3.3) Since, B(θ) = b(e iθ ) = ∞ k=-∞ b k e ikθ ,
-(z) = ∞ k=0 b - -k z -k , b + (z) = ∞ k=0 b + k z k (see
(z) = z -J b K J 0 j=1 (z -µ j ) K 0 k=1 (z -ν k ) with |µ j | < 1, |ν k | > 1. ( 3
+ e iθ ) (1 -e iθ ) = ∞ k=0 b k e iθk , (3.7)
where G(s) is the Laplace transform of the kernel g(x).

For any function w(z) with an argument on a unit circle C denote

w circ = max |z|=1 w(z).
The following lemma shows that indeed ρ 2 m log m = o(v 2 m ) as m → ∞.

w(z) = (1 -z) -r b(z), w -1 (z) = (1 -z) r b -1 (z), z = 1. (3.8)
Then, under assumptions (A1)-(A4), w(z) and w -1 (z) have no zero on the complex unit circle and, for m large enough, one has

C r 2λ 1 ( w circ ) -1 m 2r+1 ≤ v 2 m ≤ 2C r λ 2 w -1 circ m 2r+1 , (3.9)
m ρ 2 m ≤ C(r, w) v 2 m , (3.10)
where ρ 2 m and v 2 m are defined in (2.14), λ 1 and λ 2 are given by (3.2) and C(r, w) is an absolute constant which depends only on w and r:

C(r, w) = 2 4r+1 [(r -1)!] 2 w circ w -1 circ 2 λ 2 /λ 1 .

Asymptotic optimality of the estimators

Note that Lemma 3 implies that, in (2.16), ρ 2 m log m = o(v 2 m ) as m → ∞, so that the second term in (2.16) is of smaller asymptotic order than the first term. Consequently, as n → ∞, T /n → 0, the right-hand side of (2.18) is of the same asymptotic order as the oracle risk (2.12), so that, combination of Theorem 1 and Lemma 3 leads to the following statement.

Theorem 2. Let condition (2.15) hold for some positive constants α and C ρ . Then, under assumptions (A1)-(A4), for an estimator f m of f with penalty given by equation (2.16) with B > 0,

as n → ∞, R( f m ) R oracle ≤ 16(1 + B)(1 + o(1)), (3.11)
provided T /n → 0 as n → ∞.

Proof Let m 0 = arg min m [ f m -f 2 + σ 2 T n -1 v 2 m ].
Then, due to bounds (3.9) on v 2 m , one has m 0 → ∞ and Remark 1. The theory above is valid for T being finite as well as for T = T n → ∞ as long as T n /n → 0 as n → ∞. Indeed, the natural consequence of T being finite is that the bias term f -f m 2 might be relatively large due to mis-representation of f for t > T . However, since both the risk of the estimator R( f m ) and the oracle risk R oracle are equally affected, Theorem 2 remains valid whether T = T n grows with n or not.

Remark 2. The right hand side of formula (3.11) is strictly increasing in B, so, the smaller B is, the closer the risk to the optimal oracle risk as n → ∞. Note, however, that choosing asymptotically small value for B (e.g, B = 1/n) can make the second term in the penalty (2.16) dominant, so that (3.11) will become invalid.

Simulation study

In order to evaluate finite sample performance of the methodology presented above, we carried out a simulation study. We chose three versions of the kernel g, normalized to have their maximum equal to 1:

• g 1 (t) which coincides with the fit of an arterial input function (AIF) for real data obtained in the REMISCAN (2012) study. The real-life observations of an AIF corresponding to kernel g 1 coming from one patient in the REMISCAN study [START_REF]STIC 2006; Research Ethics Board (REB) approved-cohort funding by INCa ([END_REF] and fitted estimator of g 1 using an expansion over the system of the Laguerre functions with M = 18 are presented in Figure 2. One can see clearly two behavioral patterns : initial high frequency behavior caused by injection of the contrast agent as a bolus and subsequent slow decrease with regular fluctuations due to the recirculation of the contrast agent inside the blood system. one patient in the REMISCAN study [START_REF]STIC 2006; Research Ethics Board (REB) approved-cohort funding by INCa ([END_REF] and fitted estimator of g using an expansion over the system of the Laguerre functions with M = 17.

• g 2 (t) = t 2 e -0.1t which aims to reproduce a long injection of contrast agent;

• g 3 (t) = t 7 (100 + t) -1 exp -0.9t 3/4 which describes an injection with a recirculation of the contrast agent inside the blood network.

Simulations were carried out for five different test functions f :

• f 1 (x) = exp(-0.1x),
• f 2 (x) = exp(-0.6x),

• f 3 (x) = 0.5 exp(-0.1x) + 0.5 exp(-0.6x),

• f 4 (x) = 1 -IG(2; 0.5) where IG(2; 0.5) is the cdf of the gamma distribution with the shape parameter 2 and the scale parameter 0.5,

• f 5 (x) = (x + 1) -1/3 .
The value of a in formula (2.1) was chosen so that to provide the best possible fit for the kernel g when the number of terms in the expansion of g is maximum, i.e. m = M .

The functions f and g are shown in Figure 3. We illustrate performance of our methodology using kernel g 1 and test functions f 1 , ..., f 4 .

Figure 4 shows the observations and the true convolution for a medium signal-to-noise ratio 8. The associated estimators are presented in Figure 5. Here SN R is defined as

SN R = Var(f ) (σ 2 Var(g))
where, for any function ϕ, we define Var(ϕ) as

Var(ϕ) = T 0 ϕ 2 (x)dx - T 0 ϕ(x)dx 2 .
The idea of defining of SN R in this manner is to remove the effect of convolution with g. This corresponds to SNR of Abramovich and Silverman (1998).

For simulations with g 1 (t) we chose β = 1 in (1.1). We should mention that the value of β is usually unknown in real-life situations. However, since in equation (1.1), f (t) = β(1 -F (t)) where F (t) is a cdf, one knows that f (0) = β and, therefore, can estimate β as β = f (0). We executed simulations with T = 100, M = 11, two values of sample sizes, n = 100 and n = 200, and three signal-to-noise ratios (SNR), namely, SN R = 5, 8 and 15. The value 5 corresponds to real-life conditions, smaller values 8 and 15 correspond to noise level attained after the first denoising step as described in [START_REF] Rozenholc | Preserving time structures while denoising a dynamical image, Mathematical Methods for Signal and Image Analysis and Representation[END_REF].

For a given trajectory, the empirical risk was evaluated as

r( f ) = n -1 n i=1 f (t i ) -f (t i ) 2 (4.
f 3 *g 1 f 1 *g 1 f 4 *g 1 f 2 *g 1
Figure 4: Observations and true convolutions of kernel g 1 with (unknown) functions f 1 , . . . , f 4

Table 1: The values of empirical risk increased by a factor of 100: 100 × R( f ). Empirical risks are computed for 400 samples for g 1 , g 2 and g 3 and five functions of interest f 1 , . . . , f 5 . Other fine dashed lines represent the estimates for m = 0, . . . , 5. The selected value of m is given by m.

100 × R( f ) n = 100 n = 200 SN R f 1 f 2 f 3 f 4 f 5 f 1 f 2 f 3 f 4 f 5 5
Results of simulations are presented in Table 1. Table 1 verifies that indeed the methodology proposed in the paper works exceptionally well for functions f i , i = 1, . . . , 4, and is still quite precise for test function f 5 for which Fourier transform does not even exist. The table demonstrates the effect of choosing parameter a: for function f 3 and n = 100, the average empirical risk does not decline when SNR grows. This is due to the bias problem arising from the fact that f 3 is the sum of two exponentials and we fit only one value of a.

We also carried out a limited comparison of the method suggested above with the technique presented in Abramovich, Pensky and Rozenholc (2012). The comparison is performed using just one simple example where f (x) = 0.2 exp(-0.5x) + 0.8 exp(-2x) and g(t) = t 2 (t + 1)e -t (see Figure 6). In this example, the value of r in (3.1) is r = 3 and we used n = 200, σ = 0.025 and T = 15.

It is easy to see from Figure 7 We tried various regularization methods including thresholding and Tikhonov regularization with rectangular or trapezoid rules for approximation of the convolution integral and played with the constant of regularization in order to find manually the best possible tuning in each case. In Figure 8 we display one of the best reconstructions which we managed to achieve with the SVD approach.

One can clearly see how this technique fail to adequately recover unknown function f : first, it introduces a shift, second, it produces estimators which fails to be a decreasing functions (recall

that the function of interest in DCE imaging experiments is f (t) = β(1 -F (t)) where F (t)

Discussion

In the present paper, we study a noisy version of a Laplace convolution equation. Equations of this type frequently occur in various kinds of DCE imaging experiments. We propose an estimation technique for the solutions of such equation based on expansion of the unknown solution, the kernel and the measured right-hand side over a system of the Laguerre functions. The number of the terms in the expansion of the estimator is controlled via complexity penalty. The technique leads to an estimator with the risk within a logarithmic factor of m of the oracle risk under no assumptions on the model and within a constant factor of the oracle risk under mild assumptions.

The major advantage of the methodology presented above is that it is usable from a practical point of view. Indeed, the expansion results in a small system of linear equations with the matrix of the system being triangular and Toeplitz. The exact knowledge of the kernel is not required:

the AIF curve can be fitted using data from DCE-CT experiments as it is shown in Figure 2. This distinguishes the present technique with the method of Abramovich, Pensky and Rozenholc (2012) (referenced later as APR) which strongly depends on the knowledge of the kernel in general and the value of r in (3.1), in particular. After that, the method can be applied to any voxel of interest, either at the voxel level or using ROI (region of interest) manually drawn by a doctor or obtained using any clustering technique.

The method is computationally very easy and fast (requires solution of a small triangular system of linear equations) and produces no boundary effects due to extension at zero and cut-off at T . Moreover, application of the technique to discrete data does not require re-fitting the model for each model size separately. On the contrary, the vector of the Laguerre coefficients of the observed function is fitted only once, for the largest model size, and then is truncated for models of smaller sizes. The complexity of representation of g adjusts to the complexity of representation of f and the noise level. Moreover, if g can be represented by a finite expansion over Laguerre functions with k terms, the matrix of the system is k-diagonal.

The method performs very well in simulations. It is much more precise than the APR technique as Figure 7 confirms. In fact, the absence of exhaustive comparisons between the two methods is due to the fact that it is very tricky to produce estimators by the APR method, especially, in the case of g 1 which represents real life AIF. Similarly, as our study and Figure 8 show, the method is much more accurate than the SVD-based techniques.

There are few more advantages which are associated with the use of Laguerre functions basis.

Since one important goal of future analysis of DCE-CT data is classification of the tissues and The complexity of representation of g is controlled by the choice of parameter a. Parameter a is a non-asymptotic constant which does not affect the convergence rates. In practice, one can choose a in order to minimize g -ĝM where ĝM is a fitted version of g using the first M Laguerre functions. Then, the same value of a can be used in representation of the solution f . Our choice of a provides a reasonable trade-off between the bias and the variance for majority of kernels considered above, including a real life AIF kernel coming from the REMISCAN (2007) study. However, our limited experimentation with choices of a shows that there is room for improvement: undeniably, fine tuning parameter a can improve estimation precision, especially, in the case when kernel g has a strong exponential decay. However, this issue is a matter of future investigation.

γ n (t) -γ n (s) = t -f 2 -s -f 2 -2 t -s, G -1 M η m * , (6.1)
where f = f M is the vector of the true M first coefficients of function f . Note that, due to orthonormality of the Laguerre system, for any m,

fm -f 2 = f m -f 2 + ∞ j=M f (j) 2 and f m -f 2 = f m -f 2 + ∞ j=M f (j) 2 . (6.2)
Now, the definition of m yields that for any m ∈ M n one has

γ n ( f m) + pen( m) ≤ γ n ( f m ) + pen(m),
which with (6.1), implies

f m -f 2 ≤ f m -f 2 + pen(m) + ∆ m, m.
Here

∆ m, m = 2 ˆ f m -f m , G -1 M η m -pen( m)
, where m = m ∨ m. Therefore, using (6.2), we obtain that, for any m ∈ M n ,

f m -f 2 ≤ f m -f 2 + pen(m) + ∆ m, m (6.3)
Note that, due to 2xy ≤ (x 2 /4) + 4y 2 for all x > 0, y > 0,

∆ m, m ≤ 2 ˆ f m -f m sup t∈S m t =1 t, G -1 M η m -pen( m) ≤ 1 4 ˆ f m -f m 2 + 4 sup t∈S m t =1 t, G -1 M η m 2 -pen( m) Now, denote τ (m, m ) = σ 2 T n (1 + B)v 2 m * + 2(1 + B -1 )(α + 1) log(m * )ρ 2 m * , (6.4)
where m * = m ∨ m . Since, for any m,

ˆ f m -f m 2 ≤ 2 f m -f 2 + 2 f m -f 2 , then ∆ m, m ≤ 1 2 f m -f 2 + 1 2 f m -f 2 + 4 sup t∈S m∨ m t, G -1 M η m 2 -τ (m, m) + (6.5) + 4τ (m, m) -pen( m).
Using the fact that 4τ (m, m) ≤ pen(m) + pen( m), combining (6.3), (6.4) and (6.5), derive

f m -f 2 ≤ f m -f 2 + pen(m) + 1 2 f m -f 2 + 1 2 f m -f 2 + 4 sup t∈S m∨ m t, G -1 M η m 2 -τ (m, m) + .
Finally, subtracting f m -f 2 /2 from both sides of the last equation and multiplying both sides by 2, obtain

f m -f 2 ≤ 3 f m -f 2 + 4pen(m) (6.6) + 8 sup t∈S m∨ m t, G -1 M η m 2 -τ (m, m) + .
Hence, validity of Theorem 1 rests on the following lemma which will be proved later.

E sup t∈S m∨ m , t =1 t, G -1 M η m 2 -τ (m, m) + ≤ 2C 2 ρ σ 2 T mn 1 + 1 B .
Proof of Lemma 4 is given in Section 6.3.

Proofs of Lemmas 2 and 3

Proof of Lemma 2. To prove this statement, we shall follow the theory of Wiener-Hopf integral equations described in [START_REF] Gohberg | Convolution equations and projection methods for their solution[END_REF]. Denote Fourier transform of a function p(x)

by p(ω) = ∞ -∞ e iωx p(x)dx and observe that φk (ω) = (-1) k √ 2a (a + iω) k (a -iω) k+1 .
Therefore, elements of the infinite Toeplitz matrix G in (2.4) are generated by the sequence b j , j ≥ 0, where i(e iθ -1) .

b j = (2a) -1/2 (g (j) -g (j-1) ) = 1 2π ∞ -∞ ĝ(ω)[ φj (ω) -φj-1 (ω)]dω = a π ∞ -∞ ĝ(ω) iω -a iω + a j dω a 2 + ω 2 , j = 0, 1, . . . . ( 6
Now, let us show that b j = 0 for j < 0. Indeed, if j = -k, k > 0, then

b j = a π ∞ -∞ ĝ(ω) iω + a iω -a k dω a 2 + ω 2 = a π ∞ -∞ ĝ(ω) i(-ω) -a i(-ω) + a k dω a 2 + ω 2 = 1 2π ∞ -∞ ĝ(ω)[ φj (-ω) -φj-1 (-ω)]dω = ∞ -∞ g(x) [φ k (-x) -φ k-1 (-x)] dx = 0 since g(x) = 0 if x < 0 and φ k (-x) = 0 if x > 0. Hence, function B(θ) = b(e iθ
) has only coefficients b j , j ≥ 0, in its Fourier series. Now, to complete the proof, one just needs to note that G(s) = ĝ(is)

for any s such that Laplace transform G(s) of g exists.

Proof of Lemma 3. Let us first find upper and lower bounds on Then, b(z) can be written as b(z) = (1 -z) r w(z) where w(z) is defined by formula (3.8) and all zeros of w(z) lie outside the complex unit circle. Therefore, w(z) can be written as

G -1 m 2 2 = Tr(G -T m G -1 m ) and G -1 m 2 = λ max (G -T m G -1 m ).
w(z) = C w N j=1 (z -ζ j ), 0 ≤ N ≤ ∞, |ζ j | > 1, (6.8)
where C w is an absolute constant. Since b(z) does not contain any negative powers of z in its representation, J 0 = 0 and J = 0 in (3.5) and, consequently, wind(w) = 0. Also, by (3.4) and (3.8),

one has T -1 (b) = T (b -1 ) where b -1 (z) = w -1 (z)(1 -z) -r . Now, recall that G -1 m 2 2 = T m (b -1 ) 2 2 and for G -1 m 2 = T m (b -1 ) 2 .

Using relation between

Frobenius and spectral norms A 1 A 2 2 ≤ A 1 2 A 2 for any matrices A 1 and A 2 (see, e.g., [START_REF] Böttcher | Toeplitz Matrices, Asymptotic Linear Algebra, and Functional Analysis[END_REF], page 116), obtain

T m (b -1 ) 2 ≤ T m ((1 -z) -r ) 2 T m (w -1 ) , T m (b -1 ) ≤ T m ((1 -z) -r ) T m (w -1
) ,(6.9) Then, it follows from (3.6), (6.9) and (6.10) that, for m large enough, 0.5 ( w circ ) -2 ν 2 f (m) ≤ T m (b -1 ) 2 2 ≤ 2 w -1 2 circ ν 2 f (m), (6.12) 0.5 ( w circ ) -2 ν 2 s (m) ≤ T m (b -1 ) 2 ≤ 2 w -1 2 circ ν 2 s (m). (6.13) In order to finish the proof, we need to evaluate ν 2 f (m) and ν 2 s (m) and also to derive a relation between v 2 m , ρ 2 m , T m (b -1 ) 2 2 and T m (b -1 ) 2 . The first task is accomplished by the following lemma.

T m ((1 -z) -r ) 2 ≤ T m (b -1 ) 2 T m (w) , T m ((1 -z) -r ) ≤ T m (b -1 ) T m (w) . ( 6
Lemma 5. Let ν f (m) and ν s (m) be defined in (6.11). Then, Proof of Lemma 5 is given in Section 6.3. Now, to complete the proof, recall that matrix Ω m given by (2.8) is symmetric positive definite, so that there exist an orthogonal matrix U m and a diagonal matrix D m , with eigenvalues of Ω so that

ρ 2 m ≤ λ 2 T m (b -1 ) 2 , λ 1 T m (b -1 ) 2 2 ≤ v 2 m ≤ λ 2 T m (b -1 ) 2 2 .
(6.16)

Combination of (6.12) -(6.16) and Lemma 5 complete the proof.

Proofs of supplementary Lemmas

Proof of Lemma 4.

The proof of

sup t∈Sm+S m , t =1 t, G -1 M η m * 2 = Q m * ζ m * 2 .
Thus, it follows from (6.17 [(r -1)!] 2 , which proves validity of (6.14). To show that (6.15) holds, observe that, by formula 0.151.1 of [START_REF] Gradshtein | Tables of integrals, series, and products[END_REF],

m-1 j=0 r + j -1 r -1 = r + m -1 r , m r r! < r + m -1 r < m r .

  Dynamic Contrast Enhanced Computer Tomography (DCE-CT), Dynamic Contrast Enhanced Magnetic Resonance Imaging (DCE-MRI) and Dynamic Contrast Enhanced Ultra Sound (DCE-US). Such techniques are widely used in medical assessment of brain flows or cancerous tumors (see, e.g., Cao et al., 2010; Goh et al., 2005; Goh and Padhani, 2007; Cuenod et al., 2006; Cuenod et al., 2011; Miles, 2003; Padhani and Harvey, 2005 and Bisdas et al., 2007).

t 0 β 0 β 0 β

 000 AIF(t -τ ) dτ arrived before time t -t AIF(t -τ ) P (S ≤ τ ) dτ left before time t = t AIF(t -τ )(1 -F (τ ))dτ,

  )where a = -∞ and b = ∞ for Fourier convolution on a real line and -∞ < a < b < ∞ for circular convolution. Discrete stochastic version of equation(1.4)

  )). In fact, SVD has been widely used in the context of DCE imaging since mid-nineties. The technique, however, is very computationally unstable, especially, in the presence of recirculation of contrast agent. For this reason, SVD has been mostly used in the simplified framework of brain imaging due to the presence of white barrier which prevents circulation of contrast agent outside blood vessels.[START_REF] Ameloot | Extension of the performance of Laplace deconvolution in the analysis of fluorescence decay curves[END_REF] applied Laplace deconvolution for the analysis of fluorescence curves and used a parametric presentation of the solution f as a sum of exponential functions with parameters evaluated by minimizing discrepancy with the right-hand side. In a somewhat similar manner, Maleknejad et al. (2007) proposed to expand the unknown solution over a wavelet basis and find the coefficients via the least squares algorithm. Lien et al. (2008), following Weeks (1966), studied numerical inversion of the Laplace transform using Laguerre functions. Finally, Lamm (1996) and Cinzori and Lamm (2000) used discretization of the equation (1.3) and applied various versions of the Tikhonov regularization technique. However, in all of the above papers, the noise in the measurements was either ignored or treated as deterministic. The presence of random noise in (1.2) makes the problem even more challenging. For the reasons listed above, estimation of f from discrete noisy observations y in (1.2) requires extensive investigation. Unlike Fourier deconvolution that has been intensively studied in statistical literature (see references above), Laplace deconvolution received very little attention within statistical framework. To the best of our knowledge, the only paper which tackles the problem is Dey, Martin and Ruymgaart (1998) which considers a noisy version of Laplace deconvolution with a very specific kernel of the form g(t) = be -at . The authors use the fact that, in this case, the solution of the equation (1.

.

  An infinite Toeplitz matrix T = T (b) is the matrix with elements T i,j = b i-j , i, j = 0, 1, . . .. Let C = {z ∈ C : |z| = 1} be the complex unit circle. With each Toeplitz matrix T (b) we can associate its symbol b

  numbers b k are Fourier coefficients of function B(θ) = b(e iθ ). There is a very strong link between properties of a Toeplitz matrix T (b) and function b(z). In particular, if b(z) = 0 for z ∈ C and wind(b) = J b , then b(z) allows Wiener-Hopf factorization b(z) = b -(z) b + (z) z J b where b + and b -have the following forms b

  Theorem 1.8 of Böttcher and Grudsky (2005)). If T (b) is a lower triangular Toeplitz matrix, then b(z) ≡ b + (z) with b + k = b k . In this case, the product of two Toeplitz matrices can be obtained by simply multiplying their symbols and the inverse of a Toeplitz matrix can be obtained by taking the reciprocal of function b + (z): T (b + d + ) = T (b + )T (d + ), T -1 (b + ) = T (1/b + ). (3.4) Let T m (b) = T m (b + ) ∈ R m×m be a banded lower triangular Toeplitz matrix corresponding to the Laurent polynomial b(z) = m-1 k=0 b k z k . In practice, one usually use only finite, banded, Toeplitz matrices with elements T i,j , i, j = 0, 1, . . . , m -1. In this case, only a finite number of coefficients b k do not vanish and function b(z) in (3.3) reduces to a Laurent polynomial b(z) = K k=-J b k z k , z ∈ C, where J and K are nonnegative integers, b -J = 0 and b K = 0. If b(z) = 0 for z ∈ C, then b(z) can be represented in a form b

. 5 )between ρ 2 m and v 2 mLemma 2 .

 522 In this case, the winding number of b(z) is wind(b) = J 0 -J.LetT m (b) = T m (b + ) ∈ R m×m be a banded lower triangular Toeplitz matrix corresponding to the Laurent polynomial b(z) = m-1 k=0 b k z k . If b has no zeros on the complex unit circle C and wind(b) = 0, then, due to Theorem 3.7 of Böttcher and Grudsky (2005), T (b) is invertible and lim m→∞ sup T -1 m (b) < ∞. Moreover, by Corollary 3In order to apply the theory surveyed above, we first need to examine function b(z) associated with the infinite lower triangular Toeplitz matrix G defined by (2.4) and the Laurent polynomial associated with its banded version G m . It turns out that b(z) can be expressed via the Laplace transform G(s) of the kernel g(t). In particular, the following statement holds. Consider a sequence {b k } ∞ k=0 with elements b 0 = g (0) and b k = g (k) -g (k-1) , k = 1, 2, . . . where g (k) are Laguerre coefficients of the kernel g in (1.2). Then, b k , k ≥ 0, are Fourier coefficients of the function b(e iθ ) = G a(1

Lemma 3 .

 3 Let b(z) be given by (3.7), i.e., b(z) = G(a(1 + z)/(1 -z)), z = 1. Denote

  T /n → 0. Hence, it follows from Lemma 3 that ρ 2 m log m = o(v 2 m ) as m → ∞ which, in combination with Theorem 1, completes the proof.

Figure 2 :

 2 Figure 2: Observations of an arterial input function (AIF) corresponding to kernel g coming from

Figure 3 :

 3 Figure 3: Test functions : (left) the kernel functions g -(right) the estimated functions f

Figure 5 :

 5 Estimators (thick plain line) and (unknown) functions (thick dotted line) f 1 , . . . , f 4 .

  that the Laguerre functions based estimator outperforms the kernel estimator of Abramovich, Pensky and Rozenholc (2012) and also it does not exhibit boundary effects. Finally, we compared our method to Singular Value Decomposition (SVD) techniques as described in the context of DCE imaging in Ostergaard et al. (1996) and Fieselmann et al. (2011).

Figure 6 :

 6 Figure 6: Functions g (left) and f (right) for comparison between deconvolution methods using adaptive kernels and penalized Laguerre functions.

Figure 7 :

 7 Figure 7: Comparison between deconvolution using kernel method and using penalized Laguerre functions : (left) q = f * g and observations; (right) estimates f : Penalized Laguerre deconvolution (thick plain line) -adaptive kernel estimation (thick dashed line) -true function f (dotted line).

  is a cdf and we use β = 1 in our simulations). One reason for these shortcoming is that SVD estimates are smooth and degenerate at 0. As it is noted in the papers on DCE imaging (see, e.g., Fieselmann et al. (2011)), for convolution kernels corresponding to recirculation of the contrast agent (which is a common real-life scenario), SVD fails completely and needs some extra tuning in order to obtain quite poor results similar to those presented in Figure8.

observationsFigure 8 :

 8 Figure 8: Comparison between deconvolution using SVD method with Tikhonov regularization and penalized Laguerre functions. SNR=8, n = 200, f = f 1 , (upper-left) the kernel g(t) = t 3 exp(-t/3); (bottom-left) observations; (right) estimates f : Penalized Laguerre deconvolution (thick plain line) -SVD for various regularization constants (fine line) and true function f (dotted line).

  clustering of curves f (t) = β(1-F (t)) which characterize their blood flow properties, representation of the curves via Laguerre basis allows to replace the problem of classification of curves by classification of relatively low-dimensional vectors. In addition, due to the absence of boundary effects, the method allows to estimate classical medical parameters of interest β which describes the perfusion of blood flow, and also I f = f (s) ds which characterizes the vascular mean transit time. These parameters can be estimated by β = 1/ f (0) and Îf = f (s) ds, respectively.

1

 1 Let m, m ∈ M n , t ∈ S m and s ∈ S m . Denote m * = max(m, m ) = m ∨ m , η m = σ T /n ξ m and observe that

Lemma 4 .

 4 Let condition (2.15) hold for some positive constants α and C ρ . Then, for any m and any B > 0, one has

. 7 )

 7 Note that |(iω -a)/(iω + a)| = 1, so that we can use the following substitution in the integral (6.7):iω -a iω + a = e -iθ =⇒ ω = a(e iθ + 1) i(e iθ -1) = a sin θ cos θ -1 , 0 ≤ θ ≤ 2π. iθ + 1) i(+e iθ -1)e -iθj dθ, so that b j , j ∈ Z, are Fourier coefficients of the function B(θ) = b(e iθ ) = ĝ a(e iθ + 1)

-a y + a 2 = 1

 21 For this purpose, examine the function b(z) = ĝ a(z + 1) i(z -1) = G a(z + 1) 1 -z , |z| = 1. Denote y = a(z + 1)/(1 -z), so that z = (y -a)/(y + a) and G(y) = b((y -a)/(y + a)). Let us show that, under Assumptions (A1)-(A4), b(z) has a zero of order r at z = 1 and all other zeros of b(z) lie outside the unit circle. For this purpose, assume that y = α + iβ is a zero of G, i.e. G(α + iβ) = 0. Simple calculus yields y -4αa (α + a) 2 + β 2 , so that |z| = |(y -a)/(y + a)| ≤ 1 iff α ≥ 0. But, by Assumption (A3), G(y) has no zeros with nonnegative real parts, so that α < 0 and |z| = |(y -a)/(y + a)| > 1. Therefore, all zeros of b(z), which correspond to finite zeros of G, lie outside the complex unit circle C. Assumptions (A1), (A2) and properties of Laplace transform imply that G(s) = s -r (B r +G r (s)) where G r (s) is the Laplace transform of g (r) (t). Hence, lim Re s→∞ s j G(s) = 0, if j = 0, ..., r -1, B r = 0, if j = r, so that y = ∞ + iβ is zero of order r of G(y). Since lim Re y→∞ (y -a)/(y + a) = 1, b(z) has zero of order r at z = 1.

. 10 )

 10 Note that (see[START_REF] Böttcher | Spectral Properties of Banded Toeplitz Matrices[END_REF], page 13)lim m→∞ T m (w -1 ) = w -1 circ , lim m→∞ T m (w) = w circ ,Also, due to representation (6.8), both w and w -1 are bounded, and, therefore, 0< w -1 circ < ∞ and 0 < w circ < ∞. Denote ν f (m) = T m ((1 -z) -r ) 2 , ν s (m) = T m ((1 -z) -r ) . (6.11)

2 -

 2 (4r-1) [(r -1)!] -2 m 2r+1 ≤ ν 2 f (m) ≤ 0.5m 2r+1 , (6.14) (r!) -2 m 2r ≤ ν 2 s (m) ≤ m 2r . (6.15)

  Lemma 4 has two steps. The first one is the application of a χ 2 -type deviation inequality stated in[START_REF] Laurent | Adaptive estimation of a quadratic functional by model selection[END_REF], and improved by Gendre (see Lemma 3.10 of Gendre (2009)). The second step consists of integrating this deviation inequality.The χ 2 -inequality is formulated as follows. Let A be a p × p matrix A ∈ M p (R) and ζ be a standard Gaussian vector. Denote v 2 A = T r(A T A) and ρ 2 (A) = λ max (A T A). Then, for any x > 0,P Aζ 2 ≥ v 2 A + 2 v 2 A ρ 2 (A)x + ρ 2 (A)x ≤ e -x . (6.17)Now, recall that for t ∈ S m + S m = S m * where m * = m ∨ m , one hast, G -1 M η m * L = σ 2 T n t m * , Q m * ζ m *where t m * is the m * -dimensional vector formed by the first m * coordinates of t and ζ m * is a standard m * -dimensional Gaussian vector. Moreover,

) that P n T σ 2 QP Q m * ζ m * 2 - σ 2 T n ( 1 +≤ ( 1 + B - 1 ) ρ 2 m * σ 2 T n +∞ 2 ( 2 ≤ C ρ ( 1 + B - 1 )- 2 ≤ 2 = 2 Proof of Lemma 5 .r 1 2 + (m - 2 ) 2 r + 1 2 2 +-1 2 = m- 1 j=0r + j -1 j 2 (- 1 j=0r + j -1 j z j = m- 1 j=0r + j -1 r - 1 . 1 j=0r + j -1 j 2 (m -j) 2 =- 1 + 1 )j

 2211122211222522221211112211 m * ζ m * 2 ≥ v 2 m * + 2 ρ 2 m * v 2 m * x + ρ 2 m * x ≤ e -x .(6.18)For any B > 0, one has 2ρ 2 m * v 2 m * x ≤ Bv 2 m * + B -1 ρ 2 m * x so that P n T σ 2 Q m * ζ m * 2 ≥ (1 + B)v 2 m * + (1 + B -1 )ρ 2 m * x ≤ e -x .Therefore, using definition (6.4) of τ (m, m ), obtainE sup t∈Sm+S m , t =1 t, G -1 M η m * 2 -τ (m, m ) + = E Q m * ζ m * 2 -τ (m, m ) B)v 2 m * + 2(1 + B -1 )(α + 1) log(m * )ρ 2 m * ≥ ξ dξ.Changing variables2σ 2 T n (α + 1)(1 + B -1 ) log(m * )ρ 2 m * + ξ = σ 2 T n (1 + B -1 )ρ 2 m * xand application of (6.18) yieldE sup t∈Sm+S m , t =1 t, G -1 M η m * 2 -τ (m, m ) + α+1) log(m * ) e -x dx = (1 + B -1 ) ρ 2 m * σ 2 T n (m * ) -2(α+1) ≤ C ρ (1 + B -1 ) σ 2 T n (m * ) -2 .Recall that m = m ∨ m and obtainE sup t∈S m∨ m t, G -1 M η m 2 -τ (m, m) t∈Sm+S m , t =1 t, G -1 M η m * 2 -τ (m, m ) + ≤ C ρ (1 + B -1 ) σ 2 T n m ∈Mn (m ∨ m ) -C ρ (1 + B -1 ) 2C ρ (1 + B -1 ) σ 2 T nm ,which concludes the proof. Note that, by formula 1.110 of[START_REF] Gradshtein | Tables of integrals, series, and products[END_REF],(1 -z) -r = ∞ j=0 r + j -1 j z j ,so that, by definition of Frobenius norm,T m ((1 -z) -r ) 2 2 = m 2 + (m -1) 2 . . . + r + m -2 m m -j) 2 , T m ((1 -z) -r ) 2 = max |z|=1 mIf r = 1, then m-. . . (r -1 + j) (r -1)! ≤ (j + 1) r-1 ,so that, for m ≥ 4,ν 2 f (m) = T m ((1 -z) -r ) 2 2 ≤ 0.5 m 2r+1 , 2r-2 [(r -1)!] 2 (m -j) 2 ≥ m 2r+1 2 -(4r-1)

  [START_REF] Abramovich | Laplace deconvolution with noisy observations[END_REF]

	1 2 3 4	+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +	1.2 0.4 0.8	+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
	0	+ + + +						0.0	+ + +
		0	20	40	60	80	100		0	20	40	60	80	100
	0.0 0.5 1.0 1.5 2.0	+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +	0.0 0.5 1.0 1.5 2.0 2.5	+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
		0	20	40	60	80	100		0	20	40	60	80	100

Acknowledgments

Marianna Pensky was partially supported by National Science Foundation (NSF), grant DMS-1106564. The authors want to express sincere gratitude to Sergei Grudski for his invaluable help in the proof of Lemma 3 and very helpful discussions.