L Marino

S Menozzi

E Priola

Poisson process and sharp constants in L p and Schauder estimates for a class of degenerate Kolmogorov operators

Keywords: Degenerate Ornstein-Uhlenbeck operators, multidimensional parabolic equations, L p and Sobolev-space estimates, Schauder estimates, Poisson process. MSC: Primary: 35K10, 35K15, 60J76

We consider a possibly degenerate Kolmogorov-Ornstein-Uhlenbeck operator of the form L = Tr(BD 2 ) + Az, D , where A, B are N × N matrices, z ∈ R N , N ≥ 1, which satisfy the Kalman condition which is equivalent to the hypoellipticity condition. We prove the following stability result: the Schauder and Sobolev estimates associated with the corresponding parabolic Cauchy problem remain valid, with the same constant, for the parabolic Cauchy problem associated with a second order perturbation of L, namely for L + Tr(S(t)D 2 ) where S(t) is a non-negative definite N × N matrix depending continuously on t ∈ [0, T ]. Our approach relies on the perturbative technique based on the Poisson process introduced in [15].

Introduction

Let us first consider the following parabolic Cauchy problem:    ∂ t u(t, x, y) = ∆ x u(t, x, y) + x • ∇ y u(t, x, y) + f (t, x, y), u(0, x, y) = 0, (1.1) where (t, x, y) is in (0, T )×R 2d , for an integer d ≥ 1. The underlying differential operator

L K = ∆ x + x • ∇ y = d i=1 ∂ 2 x i x i + d i=1 x i ∂ y i
is the so-called Kolmogorov operator whose fundamental solution was derived in the seminal paper [START_REF] Kolmogorov | Zufällige Bewegungen (zur Theorie der Brownschen Bewegung)[END_REF]. This particular operator was also mentioned by Hörmander as the starting point for his theory of hypoelliptic operators [START_REF] Hörmander | Hypoelliptic second order differential equations[END_REF]. Let us write z = (x, y) ∈ R 2d and by ∂ z j and ∂ 2 z i z j we denote respectively the first and the second partial derivatives with i, j = 1, . . . , 2d.

We are interested in studying the influence of a second order perturbation on equation (1.1). Precisely, for a time-dependent matrix {S(t) : t ∈ [0, T ]} in R 2d ⊗ R 2d such that t → S(t) is continuous and S(t) is symmetric and non-negative definite for any fixed t, we consider the perturbed Cauchy problem:

       ∂ t u S (t, z) = L K u S (t, z) + 2d
i,j=1 S ij (t) ∂ 2 z i z j u S (t, z) + f (t, z) =: L K,S u S (t, z) + f (t, z), u S (0, z) = 0, z ∈ R 2d .

(1.2)

In particular, we will show that Sobolev (and Schauder) estimates which hold for solutions u of the Cauchy Problem (1.1) are also true, with the same constants, for solutions u S to (1.2). Clearly, the operator L K,S can be seen as a perturbation of L K involving second order partial derivatives with continuous time-dependent coefficients.

For now, let us explain our main results in a special form for equation (1.1) in the case of L p -estimates (or Sobolev estimates). For a statement of our results in the whole generality, we instead refer to Section 2. For a fixed final time T > 0 and a source f in C ∞ 0 ((0, T ) × R 2d ), it is known from the work of Bramanti et al. [START_REF] Bramanti | L p -estimates for some ultraparabolic operators with discontinuous coefficients[END_REF], Theorem 3.1 (see also Section 2.3 below), that equation (1.1) admits a unique classical bounded solution u which satisfies for p in (1, +∞) the following estimates: ∆ x u L p ((0,T )×R 2d ) ≤ C p f L p ((0,T )×R 2d ) = C p ∂ t u -L K u L p ((0,T )×R 2d ) . (1.3) Note that in this case C p = C p (d) > 0. We will actually manage to prove that the unique classical bounded solution u S to (1.2) satisfies the estimate ∆ x u S L p ((0,T )×R 2d ) ≤ C p f L p ((0,T )×R 2d ) = C p ∂ t u S -L K,S u S L p ((0,T )×R 2d ) , (1.4) with the same previous constant C p as in (1.3). This result seems to be new even in dimension 2d = 2 and even if we only consider S(t) = S, t ∈ [0, T ], where S is a 2 × 2 symmetric non-negative definite matrix. For a uniformly elliptic second order perturbation S(t) = S, t ∈ [0, T ], where S is positive definite, we could also have appealed to [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] to derive estimates like in (1.4). For related estimates in the uniformly elliptic case, see also Section 4 in Metafune et al. [START_REF] Metafune | The domain of the Ornstein-Uhlenbeck operator on an L p -space with invariant measure Ann[END_REF]. However, note that from [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] and [START_REF] Metafune | The domain of the Ornstein-Uhlenbeck operator on an L p -space with invariant measure Ann[END_REF] we could only deduce that the constant C p depends on the ellipticity constant of the perturbation (this is the first eigenvalue λ S of S) and on the maximum eigenvalue of S (on this respect, see also [START_REF] Krylov | Parabolic equations in L p -spaces with mixed norms[END_REF] and [START_REF] Priola | L p -parabolic regularity and non-degenerate Ornstein-Uhlenbeck type operators[END_REF]).

The remarkable point in (1.4) is that the L p -estimates are stable under second order perturbations, which can be possibly degenerate. Namely, the fact that S(t) might be degenerate for some t in (0, T ), or even in some non-empty sub-intervals of (0, T ), does not affect the estimates in (1.4).

To prove (1.4), we combine the results of [START_REF] Bramanti | L p -estimates for some ultraparabolic operators with discontinuous coefficients[END_REF] with a probabilistic perturbative approach based on the Poisson process inspired by [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF]. There, it was established in particular that the L p -estimates for non-degenerate parabolic heat equations with space homogeneous coefficients are valid with constants that are independent of the dimension.

Remark 1.1. Importantly, the approach of [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] turns out to be sufficiently robust to handle the estimates in the degenerate directions as well. We recall that the associated maximal L p -regularity was studied e.g. in [START_REF] Bouchut | Hypoelliptic regularity in kinetic equations[END_REF], [START_REF] Huang | L p estimates for degenerate non-local Kolmogorov operators[END_REF] or [START_REF] Chen | Propagation of regularity in L p -spaces for Kolmogorov-type hypoelliptic operators[END_REF]. Let p in (1, +∞), there exists Cp > 0 such that for f in C ∞ 0 ((0, T ) × R 2d ) the unique classical bounded solution u of (1.1) verifies

(∆ y ) 1 3 u L p ((0,T )×R 2d ) ≤ Cp f L p ((0,T )×R 2d ) = Cp ∂ t u -L K u L p ((0,T )×R 2d ) ,
(1.5) where (∆ y ) 1 3 denotes the fractional Laplacian with respect to the degenerate variables y in R d . It turns out that this estimate is also stable for the previously described second order perturbation. Namely, for u S solving (1.2), (∆ y )

1 3 u S L p ((0,T )×R 2d ) ≤ Cp f L p ((0,T )×R 2d ) = Cp ∂ t u S -L K,S u S L p ((0,T )×R 2d ) ,
(1.6) where again Cp is the same as in (1.5).

Remark 1.2. The same type of stability results will also hold for the corresponding global Schauder estimates, first established in the framework of anisotropic Hölder spaces for the solution of (1.1) by Lunardi [START_REF] Lunardi | Schauder estimates for a class of degenerate elliptic and parabolic operators with unbounded coefficients in R n[END_REF] (see also [START_REF] Marino | Schauder estimates for degenerate stable Kolmogorov equations[END_REF], [START_REF] Marino | Schauder estimates for degenerate Lévy Ornstein-Uhlenbeck operators[END_REF] and the references therein). We refer to estimate (4.17).

We point out that our results in Section 3 could also possibly be obtained by using the general theorems of Section 4 in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF]. This section in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] introduces a more general probabilistic approach and provides unexpected regularity results. However checking in our case all the assumptions given in that section is quite involved. On the other hand, we provide self-contained proofs inspired by Sections 2 and 3 of [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF].

It remains a challenging open problem to have a purely analytic proof of our regularity results.

L p -estimates for degenerate Ornstein-Uhlenbeck operators. Let us now describe the more general framework we are going to consider here.

Let R N = R d 0 × R d 1 where d 0 , d 1 are two non-negative integers such that d 0 + d 1 = N and d 0 ≥ 1. let us introduce the non-negative, symmetric matrix B in R N ⊗ R N given by B = B 0 0 0 0 , where B 0 is a symmetric, positive definite matrix in R d 0 ⊗ R d 0 such that ν d 0 i=1 ξ 2 i ≤ d 0 i,j=1 (B 0 ) ij ξ i ξ j ≤ 1 ν d 0 i=1 ξ 2 i ,
for all ξ ∈ R d 0 , for some ν > 0. We will use, as underlying proxy operators, the family of degenerate Ornstein-Uhlenbeck generators of the form

L ou f (z) = Tr(BD 2 f (z)) + Az, Df (z) , z = (x, y) ∈ R d 0 +d 1 = R N , (1.7)
for a matrix A in R N ⊗ R N , where •, • denotes the usual inner product in R N . Moreover, we assume the Kalman condition:

[K] There exists a non-negative integer k, such that

Rank[B, AB, . . . , A k B] = N, ( 1.8)
where [B, AB, ..., A k B] is the R N ⊗ R N (k+1) matrix whose blocks are B, AB, . . . , A k B. From the non-degeneracy of B 0 , the above condition amounts to say that the vectors

{e 1 , . . . , e d 0 , Ae 1 , . . . , Ae d 0 , . . . , A k e 1 , . . . , A k e d 0 } generate R N , ( 1.9)
where {e i } i∈{1,••• ,d 0 } are the first d 0 vectors of the canonical basis for R N . Assumption [K] (which also often appears in control theory; see e.g. [START_REF] Zabczyk | Mathematical control theory: an introduction[END_REF]) is equivalent to the Hörmander condition on the commutators (c.f. [START_REF] Hörmander | Hypoelliptic second order differential equations[END_REF]) ensuring the hypoellipticity of the operator ∂ t -L ou . In particular, it implies the existence and the smoothness of a distributional solution for the following equation:

   ∂ t u(t, z) = L ou u(t, z) + f (t, z), on (0, T ) × R N ; u(0, z) = 0, on R N , (1.10)
where

f is a function in C ∞ 0 ((0, T ) × R N ).
Similarly to [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF], we will prove below the existence and uniqueness of bounded regular solutions to (1.10) assuming that the source f belongs to the space B b 0, T ; C ∞ 0 (R N ) , which contains C ∞ 0 ((0, T ) × R N ), and that can be roughly described as the family of functions which are bounded measurable in time and compactly supported in space uniformly in time (see Section 1.2 for a precise definition). Equation (1.10) will be understood in an integral form (cf. formula (1.25)).

By Theorem 3 in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] and exploiting some explicit properties of the underlying heat kernel (see Section 2.3 below), it can be derived that for any fixed p in (1, +∞), there exists

C p = C p (ν, A, d 0 , d 1 , T ) such that D 2 x u L p ((0,T )×R N ) ≤ C p ∂ t u -L ou u L p ((0,T )×R N ) = C p f L p ((0,T )×R N ) , (1.11) where for any z in R N , t ∈ [0, T ], D 2
x u(t, z) stands for the Hessian matrix in R d 0 ⊗ R d 0 with respect to the variable x. We set

B I = I d 0 ,d 0 0 d 0 ,d 1 0 d 1 ,d 0 0 d 1 ,d 1
and note, in particular, that (1.11) can be rewritten in the following, equivalent way:

B I D 2 u B I L p ((0,T )×R N ) = D 2 x u L p ((0,T )×R N ) ≤ C p ∂ t u -L ou u L p ((0,T )×R N ) = C p f L p ((0,T )×R N ) ,
(1.12)

where D 2 u = D 2 z u represents instead the full Hessian matrix in R N ⊗ R N with respect to z.

Fixed a continuous mapping t → S(t) such that S(t) is a symmetric and non-negative definite matrix in R N ⊗ R N , t ∈ [0, T ], we consider again the following perturbation of L ou :

L ou,S t f (z) := Tr(BD 2 f (z)) + Tr(S(t)D 2 f (z)) + Az, Df (z) = L ou f (z) + Tr(S(t)D 2 f (z)), (1.13)
where

z = (x, y) is in R d 0 +d 1 = R N .
For the solution u S of the related Cauchy problem

   ∂ t u S (t, z) = L ou,S t u S (t, z) + f (t, z), on (0, T ) × R N ; u S (0, z) = 0, on R N , (1.14)
we will prove the following main theorem:

Theorem 1.1. Let us consider (1.14) with f ∈ B b 0, T ; C ∞ 0 (R N ) .
Then, there exists a unique solution u S of Cauchy Problem (1.14) which verifies, with the same constant C p as in (1.12),

D 2 x u S L p ((0,T )×R N ) = B I D 2 u S B I L p ((0,T )×R N ) (1.15) ≤ C p ∂ t u S -L ou,S t u S L p ((0,T )×R N ) = C p f L p ((0,T )×R N ) .
We point out that for time-homogeneous non-negative definite matrices S, the corresponding elliptic L p -estimates as in formula (5) of [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] (replacing A in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] with L ou,S := Tr(BD 2 •)+Tr(SD 2 •)+ Az, D• ) with a constant independent of S, could also be derived from (1.15) using an argument given in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF].

For more information on the OU operator L ou we also refer to the recent work by Fornaro et al. [START_REF] Fornaro | L p -spectrum of degenerate hypoelliptic Ornstein-Uhlenbeck operators[END_REF] about full description of the spectrum of degenerate OU operators in L p -spaces.

Independently from the constant preservation, we also emphasize that the L p -estimates in (1.15) for the perturbed operator seem, to the best of our knowledge, to be new and have some interest by their own.

Let us eventually mention that our stability results could turn out to be useful to investigate the well-posedness of some related stochastic differential equations through the corresponding martingale problem.

We could actually derive more general estimates, possibly involving the degenerate directions as well, dependingly on the structure of A. Some results in that direction are gathered in Section 4. Anyhow, to illustrate our approach we now briefly present the various steps to derive (1.15).

Strategy of the proof for estimate (1.15).

Fixed a classical bounded solution u to Cauchy Problem (1.10), let us introduce v(t, z) := u(t, e -tA z). This well-known transformation (cf. [START_REF] Da Prato | On the Ornstein-Uhlenbeck operator in spaces of continuous functions[END_REF]) precisely allows to get rid of the drift term in the PDE satisfied by v. Indeed, we have that u(t, z) = v(t, e tA z) and since u solves (1.10), it holds for any (t, z) in (0, T ) × R N , that:

f (t, z) = ∂ t u(t, z)-L ou u(t, z) = v t (
   ∂ t v(t, z) = Tr e tA Be tA * D 2 v(t, z) + f (t, z) on (0, T ) × R N ; v(0, z) = 0 on R N . (1.17)
In terms of the function v, the known estimates in (1.12) rewrites as:

B I e tA * D 2 v(t, e tA •) e tA B I L p ((0,T )×R N ) ≤ C p f (t, e tA •) L p ((0,T )×R N ) , (1.18)
where we used the notation B I e -tA * D 2 v(t, e tA •) e tA B I L p ((0,T )×R N ) to stress the dependence on t instead of the more precise formulation

B I e • A * D 2 v(•, e • A •) e • A B I L p ((0,T )×R N ) .
By changing variable in the integrals, control (1.18) is equivalent to

B I e tA * D 2 v(t, •) e tA B I L p ((0,T )×R N ,m) ≤ C p f L p ((0,T )×R N ,m) , (1.19)
where L p ((0, T )×R N , m) denotes the L p -norms w.r.t. the measure m(dt, dx) = det(e -At )dtdx.

Considering now the following more general Cauchy problem on

[0, T ] × R N    ∂ t w(t, z) = Tr e tA Be tA * D 2 w(t, z) + Tr e tA S(t)e tA * D 2 w(t, z) + f (t, z); w(0, z) = 0, ( 1
.20) we can establish the well-posedness of the Cauchy problem (1.20), exploting, for instance, probabilistic arguments, using the underlying Gaussian process. Now the crucial step consists in adapting some arguments from [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] based on the use of the Poisson process to derive that the same L p -estimates in (1.19) still hold for w, independently from the non-negative definite, symmetric matrices S(t). Precisely,

B I e tA * D 2 w(t, •) e tA B I L p ((0,T )×R N ,m) ≤ C p f (t, •) L p ((0,T )×R N ,m) , (1.21)
with the same constant C p appearing in (1.19).

The last step then consists in coming back to the Ornstein-Uhlenbeck operators framework. Namely, we introduce ũ(t, z) := w(t, e tA z) which solves, by definition, the following equation:

   ∂ t ũ(t, z) = L ou,S t ũ(t, z) + f (t, z), (t, z) ∈ (0, T ) × R N , ũ(0, z) = 0, z ∈ R N . Thus ũ = u S . Noticing that D 2 w(t, •) = D 2 [ũ(t, e -tA •)] = e -tA * D 2 ũ(t, e -tA •)e -tA
we thus get from (1.21) that the following estimates hold:

B I D 2 ũ B I L p ((0,T )×R N ) ≤ C p f L p ((0,T )×R N ) . (1.22)
Through the previous steps we have then constructed a solution ũ of Cauchy Problem (1.14) which indeed satisfies the estimates in (1.15) with the same C p , associated with the unperturbed or proxy operator. The maximum principle will eventually provide uniqueness for the solution ũ.

Remark 1.3. i)

We point out that we could also consider more general time-dependent Ornstein-Uhlenbeck operators like:

M = Tr(B(t)D 2 •) + Az, D• .
Arguing as before starting from L p -estimates (or Schauder estimates) for M we can derive the same L p -estimates (or Schauder estimates) for a perturbation of M like (1.13).

ii) We could extend the L p -estimates (or the Schauder estimates) related to L ou to more general operators like

L ou,S t f (z) + b(t), Df (z)
where b : R + → R N is continuous. We can even add to L t a possibly degenerate non-local perturbation (cf. Section 7 of [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF]). The L p -estimates (or Schauder estimates) are still preserved with the same constant. For the sake of simplicity in the sequel we will only consider b(t) = 0 and we will not deal with non-local perturbations of L ou,S t .

Organization of the paper. The article is organized as follows. At the end of the current section, we first fix some useful notations. In Section 2 we will then focus on driftless second order Cauchy problems associated with a nonnegative definite, possibly degenerate, diffusion matrix. We will also consider its relation to the Ornstein-Uhlenbeck dynamics. We will establish through the probabilistic perturbation approach of [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] that if some L p -estimate holds for a particular diffusion matrix so does it, with the same associated constant as explained before, for a non-negative perturbation of the diffusion matrix (see Section 3). Finally, by the arguments of Section 1.1 we will obtain (1.22). Stability results in anisotropic Sobolev space and Schauder estimates are given in Section 4.

Definition of solution and useful notations

Let us consider the following Cauchy problem:

   ∂ t v(t, z) = tr (Q(t)D 2 v(t, z)) + b(t, z), Dv(t, z) + f (t, z), on (0, T ) × R N ; v(0, z) = 0, on R N ; (1.23) where Q : [0, T ] → R N ⊗ R N is a continuous symmetric non-negative definite matrix and b : [0, T ] × R N → R N is a continuous function such that |b(t, z)| ≤ K T (1 + |z|), (t, z) ∈ [0, T ] × R N , for some constant K T > 0. The function f belongs to B b 0, T ; C ∞ 0 (R N ) , the space of all Borel bounded functions φ : [0, T ] × R N → R such that φ(t, •) is smooth and com- pactly supported for any t in [0, T ]; for any n in N the C n (R N )-norms of φ(t, •)
are bounded in time and the supports of the functions φ(t, •) are contained in the same ball. Moreover, we require that, for any z ∈ R N , the mapping:

t → φ(t, z) (1.24)
is a piece-wise continuous function on [0, T ], i.e. it is continuous except for a finite number of points.

Remark 1.4. Note that to perform the technique used in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] and based on the Poisson process we need to consider equations like (1.23) with a source f which is possibly discontinuous in time (cf. the proof in Section 2 of [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] and Section 3.2 below).

We interpret Cauchy Problem (1.23) in an integral form:

v(t, z) = t 0 f (s, z) + Tr(Q(s)D 2 v(s, z)) + b(s, z), Dv(s, z) ds. (1.25)
In particular, we say that a continuous and bounded function v

: [0, T ]×R N → R is a solution to equation (1.23) if v(t, •) belongs to C 2 (R N ), for any t ∈ [0, T ],
and (1.25) holds as well, for any (t, z) ∈ [0, T ] × R N . We finally note that, for any

z ∈ R N , the function t → v(t, z) is a C 1 -piece- wise function on [0, T ].
By Theorem 4.1 in [START_REF] Krylov | Elliptic and parabolic second-order PDEs with growing coefficients[END_REF] we deduce in a quite standard way that if a solution v exists then it is unique and the following maximum principle holds:

sup (t,z)∈[0,T ]×R N |v(t, z)| ≤ T sup (t,z)∈[0,T ]×R N |f (t, z)|. (1.26)
About the proof of (1.26) we only make some remarks. By considering v and -v we see that it is enough to prove that v(t, z)

≤ T f ∞ , for all (t, z) ∈ [0, T ] × R N . Moreover, setting ṽ = v -t f ∞ , we note that ṽ verifies (1.25) with f replaced by f -f ∞ ≤ 0.
Finally, by considering the equation verified by e -t ṽ, we can apply Theorem 4.1 in [START_REF] Krylov | Elliptic and parabolic second-order PDEs with growing coefficients[END_REF] to obtain the result.

Estimates for driftess second order operators and related perturbation

Throughout this section, we consider the following Cauchy problem:

   ∂ t v(t, z) = Tr (Q(t)D 2 v(t, z)) + f (t, z) on (0, T ) × R N ; v(0, z) = 0 on R N , ( 2.1)
which can be seen as a special case of (1.23) when b = 0. Moreover, we assume that Q is not identically zero.

Well-posedness

Proposition 2.1 (Well-posedness in integral form for the driftless Cauchy problem). Let f be in B b 0, T ; C ∞ 0 (R N ) . Then, there exists a unique solution v to Cauchy problem (2.1) in an integral sense, i.e., it solves for

(t, z) ∈ [0, T ] × R N : v(t, z) = t 0 f (s, z) + Tr(Q(s)D 2 v(s, z) ds.
(2.2)

We will denote in short v = P DE(Q, f ).

Proof. By the maximum principle (cf. equation (1.26)) uniqueness holds for Cauchy Problem (2.1). We can then focus on proving the existence of a solution. Let us introduce now

v(t, z) := t 0 E[f (s, z + I s,t )] ds
with the following notation:

I s,u := √ 2 u s Q(r) 1/2 dW r ,
where W is an Ndimensional Brownian motion on some probability space (Ω, F, (F t ) t≥0 , P) and

Q(r) 1/2 stands for a square root of Q(r), i.e. Q(r) = Q(r) 1/2 (Q(r) 1/2 ) * .
Applying the Itô formula in space to f (s, z + I s,u ) u∈[s,t] , we get that

Ef (s, z + I s,t ) = f (s, z) + E t s Tr(Q(u)D 2 f (s, z + I s,u ))du .
Hence,

v(t, z) = t 0 f (s, z) + E[ t s Tr(Q(u)D 2 f (s, z + I s,u ))du] ds,
from which it readily follows that

∂ t v(t, z) = f (t, z) + t 0 E[Tr(Q(t)D 2 f (s, z + I s,t ))]ds = f (t, z) + Tr Q(t)D 2 t 0 E[f (s, z + I s,t )]ds = f (t, z) + Tr Q(t)D 2 v(t, z) .
for almost every t ∈ [0, T ] and any z ∈ R N .

Relation to the Ornstein-Uhlenbeck dynamics

   ∂ t u(t, z) = L ou u(t, z) + f (t, z), on (0, T ) × R N ; u(0, z) = 0, on R N ; (2.3) with L ou as in (1.7), f (t, z) = f (t, e tA z). Precisely, for all (t, z) ∈ [0, T ] × R N , u(t, z) = t 0 f (s, z) + L ou v(s, z) ds. (2.4)
We have that u is a solution to (2.3). Let us also point out that the well-posedness of (2.3) could also have been obtained directly from Gaussian type calculations, similar to those in the proof of Proposition 2.1, introducing u ou (t, z) := t 0 E[ f (s, e (t-s)A z + I ou s,t )]ds where

I ou s,u := √ 2 u s e (u-v)A BdW v .

About the L p -estimate (1.11) for the OU operator

The aim of this section is to fully justify the estimates in (1.11). This is a consequence of the previous probabilistic representation and of Theorem 3 in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF]. For u solving (1.10) it holds that for all (t, z)

∈ [0, T ] × R N , u(t, z) = t 0 E[f (s, e A(t-s) z + I ou s,t )]ds = t 0 R N f (s, z )p ou (t -s, z, z )dz ds,
(2.5) where for v > 0, p ou (v, z, •) stands for the density at time v of the stochastic process

X ou u := e Au z + √ 2 u 0 e A(u-w) BdW w = z + u 0 AX ou w dw + √ 2BW u , u ≥ 0.
We recall from [START_REF] Lanconelli | On a class of hypoelliptic evolution operators[END_REF] that assumption [K] is equivalent to the fact that there exists k ∈ N and positive integers (

d i ) i∈{1,••• ,k} s.t. k i=1 d i = d 1
and for all i ∈ {1, • • • , k}, setting d 0 = d 0 and -1 m=0 = 0, the matrixes

A i := (A j, ) (j, )∈{ i-1 m=0 dm+1,••• , i m=0 dm}×{ i-1 m=1 dm+1,••• , i m=1 dm} , have rank d i . The matrix A writes: A =         
0 d k ,d 0 . . . 0 d k ,d k-1 A k *           . (2.6)
Following the proof of Lemma 5.5 in [START_REF] Delarue | Density estimates for a random noise propagating through a chain of differential equations[END_REF], where the case

d 0 = d, d i = d, k = n -1 is addressed, it can be derived that there exists C ≥ 1 s.t. for all (v, z, z ) ∈ (0, T ] × (R N ) 2 , |D 2 x p ou (v, z, z )| ≤ C v k-1 i=0 d i (i+ 1 2 )+1 exp -C -1 v|T -1 v (e Av z -z )| 2 , (2.7)
where d 0 = d 0 and

T v := diag(vI d 0 ×d 0 , v 2 I d 1 ×d 1 , . . . , v k+1 I d k ×d k ), v ≥ 0,
reflects the various scales of the system. For a given function f ∈ B b 0, T ; C ∞ 0 (R N ) , it is then clear from (2.5) and (2.7) that for all (t, z) ∈ (0, T ] × R N :

D 2 x u(t, z) = p.v. t 0 R N f (s, z )D 2 x p ou (t -s, z, z )dz ds. (2.8)
It indeed suffices to observe that:

|p.v. t 0 R N f (s, z )D 2 x p ou (t -s, z, z )dz ds| = |p.v. t 0 R N [f (s, z ) -f (s, e A(t-s) z)]D 2 x p ou (t -s, z, z )dz ds| ≤ (2.7) sup s∈[0,T ] Df (s, •) ∞ × t 0 R N C (t -s) k-1 i=0 d i (i+ 1 2 )+ 1 2 exp -C -1 (t -s)|T -1 t-s (e A(t-s) z -z )| 2 dz ds ≤ C sup s∈[0,T ] Df (s, •) ∞ T 1 2 .
The estimates in (1.11) now follows from the proof of Theorem 3 in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF], starting from (2.8) instead of ( 16) therein. The strategy is clear. It is necessary to introduce a cut-off function which separates the points (s, z ) which do not induce any singularity in (2.8) for the derivatives of the density, namely such that t -s ≥ c 0 or |e A(t-s) z -z | ≥ c 0 , for some fixed constant c 0 > 0, from those who are close to the singularity. For the non-singular part of the integral the expected L p -control readily follows from (2.7) and the Young inequality (see also Proposition 5 in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF]), whereas the derivation of the bound for the singular part requires some involved harmonic analysis, see Section 4 on the same reference. We can also refer to Theorem 11 and its proof in [START_REF] Priola | On weak uniqueness for some degenerate SDEs by global L p estimates[END_REF] for similar issues linked with the corresponding L p -estimates for degenerate Ornstein-Uhlenbeck operators in an elliptic setting.

The main result for equation (2.1): perturbation of second order driftless PDE

Let us fix p in (1, +∞) and assume that there exists R(t) ∈ R N ⊗ R N depending continuously on t ≥ 0 and a constant C p > 0, such that for any f in

B b 0, T ; C ∞ 0 (R N ) , the unique solution v = P DE(Q, f ) to equation (2.1) satisfies R(t) * D 2 v R(t) L p ((0,T )×R N ,m) ≤ C p f L p ((0,T )×R N ,m) , ( 2.9)
for some absolutely continuous measure m w.r.t. the Lebesgue measure on [0, T ] × R N such that m(dt, dx) = g(t)dtdx for some Borel bounded function g (note that in (1. [START_REF] Marino | Schauder estimates for degenerate Lévy Ornstein-Uhlenbeck operators[END_REF] we have R(t) = e tA B I , m(dt, dx) = g(t)dtdx = det(e -At )dtdx).

We would like to exhibit that a control like (2.9) also holds for the solution w to the following Cauchy Problem:

   ∂ t w(t, z) = tr (Q(t)D 2 w(t, z)) + tr (Q (t)D 2 w(t, z)) + f (t, z), on (0, T ) × R N ; w(0, z) = 0, on R N , ( 2
⊗ R N and f ∈ B b 0, T ; C ∞ 0 (R N )
. Assume that estimate (2.9) holds as explained above. Then the solution w to (2.10) verifies

R(t) * D 2 w R(t) L p ((0,T )×R N ,m) ≤ C p f L p ((0,T )×R N ,m) ,
(2.11) p ∈ (1, ∞) with the same constant C p as in (2.9).

From Theorem 2.2 using the argument of Section 1.1 we can easily derive Theorem 1.1.

A perturbation argument for proving Theorem 2.2

We aim here at applying the probabilistic perturbative approach considered in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF]. The key idea in that work was, for a well-posed PDE which enjoys some quantitative given estimates, to introduce a small random perturbation in the source f through a suitable Poisson type process and to investigate the properties of the associated PDE involving an unknown function v. After considering a small random perturbation of v, we arrive at the useful integral formula (3.8). Taking the expectation the contributions associated with the jumps yield, for an appropriate intensity of the underlying Poisson process, a finite difference operator. For the PDE satisfied by the expectation, involving the finite difference operator, the initial estimates are preserved. Repeating the previous argument we can obtain a PDE involving the composition of two finite difference operators. Compactness arguments then allow to derive that, the initial estimates still hold at the limit with the composition of two finite difference operators replaced by the corresponding differential operator of order two. Iterating this procedure we can obtain the result.

Below, we start recalling basic properties of Poisson type processes and corresponding stochastic integrals, which are needed for our approach.

Poisson stochastic integrals

We briefly recall here the very definition of the stochastic integral driven by a Poisson process. We start reminding the construction of such processes. Given a probability space (Ω, F, P) to be fixed from this point further, we start considering a sequence of independent real-valued random variables {τ n } n∈N on Ω whose distribution is exponential of parameter λ > 0:

P(τ n > r) = e -rλ , r ≥ 0.
We can then define the partial sums sequence {σ n } n∈N as follows:

σ 0 = 0; σ n = n i=1 τ i , n = 1, 2, . . .
For any fixed t ≥ 0, π t now denotes the number of consecutive sums of τ i which lie on [0, t], i.e.

π t = ∞ n=0 1 σn≤t , ( 3.1)
where 1 σn≤t represents the indicator function of the event {σ n ≤ t}. The process {π t } t≥0 we have just constructed is usually known in the literature as a Poisson process with intensity λ (see, for instance, [START_REF] Protter | Stochastic Integration and Differential Equations[END_REF]). Now, let c : [0, T ] → R N be a continuous function. We can define the Poisson stochastic integral as

b t := t 0 c(s)dπ s = σ k ≤t, k≥1 c(σ k ) = 0<s≤t c(s)(π s -π s-), t ∈ [0, T ], (3.2)
b 0 = 0 (as usual π s-(ω) denotes the left limit at s, for any ω, P-a.s.). We now recall the following formula for the expectation of the stochastic integral:

E[ t 0 c(s)dπ s ] = λ t 0 c(s)ds ∈ R N .
(3.3) (cf. Lemma 2.1 in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] for a direct proof; see also Theorem 16 in [START_REF] Protter | Stochastic Integration and Differential Equations[END_REF] and Theorem 5.3 in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] for a more general formula involving stochastic integrals of predictable processes against the Poisson process). We also recall the following more general result.

Proof of Theorem 2.2

According to the notations in Proposition 2.1, let v = P DE(Q, f ) and w = P DE(Q + Q , f ) be the unique solutions of equations (2.1) and (2.10), respectively. The proof of Theorem 2.2 will be obtained adapting the method developed in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] (see in particular Section 3, therein). Let e 1 be the first unit vector in R N . We define

X t = t 0 Q (r) e 1 dπ r
where Q (t) is the unique N × N symmetric non-negative definite square root of Q (t) and {π t } t≥0 is a Poisson Process of intensity λ (cf. (3.2)). The parameter λ will be chosen appropriately later on.

Recall that the solution v to (2.1) is given by

v(t, z) = t 0 ds R N [f (s, z + z )µ s,t (dz ) (3.5)
where µ s,t is the Gaussian law of the stochastic integral I s,t := √ 2 t s Q(r) dW r (see the proof of Proposition 2.1).

Let us fix > 0. We notice that the shifted source f (t, z) := f (t, z -X t ) (which also depends on ω; we have omitted to write such dependence on ω) is again in B b 0, T ; C ∞ 0 (R N ) . This is the reason why we have considered such a function space for the source. It precisely allows to take into account the time discontinuities coming from the jumps of the Poisson process.

For any fixed ω in Ω, Proposition 2.1 readily gives that there exists a unique solution v = PDE(Q, f (t, z -X t )), depending on and ω as parameters, such that sup

(t,z)∈[0,T ]×R N |v (t, z)| ≤ T sup (t,z)∈[0,T ]×R N |f (t, z)|. (3.6)
Moreover, thanks to the invariance for translations of the L p -norms, it follows from (2.9) that

R(t) * D 2 v R(t) L p ((0,T )×R N ,m) ≤ C p f L p ((0,T )×R N ,m) = C p f L p ((0,T )×R N ,m) .
(3.7) By Equation (3.5), we know that v is given by

v (t, z) = t 0 ds R N [f (s, z -X s + z )µ s,t (dz ).
For each z ∈ R N , the stochastic process (v (t, z)) t∈[0,T ] has continuous paths (P-a.s.) and it is F t -adapted where F t is the completed σ-algebra generated by the random variables π s , 0 ≤ s ≤ t.

For fixed z ∈ R N , > 0, let us introduce the process (v (t, z + X t )) t∈[0,T ] which is given by

v (t, z + X t ) = t 0 ds R N [f (s, z + X t -X s + z )µ s,t (dz ).
It is not difficut to check that it is F t -adapted and it has càdlàg paths.

Applying (2.2) on each interval [σ n , σ n+1 ∧ t), n ∈ {0, • • • , π t ] on which X s is constant, one then derives that:

v (t, z + X t ) = t 0 tr(Q(s)D 2 v (s, z + X s )) + f (s, z) ds + t 0 g (s, z) dπ s , (3.8) where g (s, z) = v (s, z + Q (s) e 1 + X s-) -v (s, z + X s-
) is precisely the contribution associated with the jump times. It is clear that g (s, z) = 0 if and only if π s has a jump at time s. We then have by Lemma 3.1:

E t 0 g (s, z) dπ s = λ t 0 v (s, z + Q (s) e 1 ) -v (s, z) ds, (3.9) where v (s, z) = E[v (s, z + X s )]. Let us denote l(t) := Q (t) e 1 .
Taking the expectation on both sides of equation (3.8), we find out that v is an integral solution of the following PDE:

∂ t v (t, z) = tr(Q(t)D 2 v (t, z)) + λ (v (t, z + l(t)) -v (t, z)) + f (t, z), (3.10)
with zero initial condition. Remark that uniqueness of bounded continuous solutions to (3.10) follows by the maximum principle, arguing as in the proof of Lemma 2.2 in [START_REF] Krylov | Poisson stochastic process and basic Schauder and Sobolev estimates in the theory of parabolic equations[END_REF] (first one considers the case λT ≤ 1/4 and then one iterates the procedure by steps of size 1/(4λ)).

Moreover by (3.7) we obtain (using also the Jensen inequality and the Fubini theorem)

R(t) * D 2 v R(t) p L p ((0,T )×R N ,m) = (0,T )×R N |R(t) * D 2 v (t, z) R(t)| p m(dt, dz) = (0,T )×R N |E[R(t) * D 2 v (t, z + X t )R(t)] | p dzg(t)dt ≤ [0,T ]×R N E[|R(t) * D 2 v (t, z + X t )R(t) | p ]dzg(t)dt = E [0,T ]×R N |R(t) * D 2 v (t, z + X t )R(t) | p dzg(t)dt = E [0,T ]×R N |R(t) * D 2 v (t, z)R(t) | p dzg(t)dt ≤ C p p f p L p ((0,T )×R N ,m) ,
using (3.7) for the last inequality (L p -estimate for the PDE with random source). Choosing λ = -2 we have from (3.10)

∂ t v (t, z) = tr(Q(t)D 2 v (t, z)) + -2 (v(t, z + l(t)) -v (t, z)) + f (t, z), ( 3

.11) with zero initial condition and moreover

R(t) * D 2 v R(t) p L p ((0,T )×R N ,m) ≤ C p p f p L p ((0,T )×R N ,m) . (3.12)
Now the idea is to apply again the same reasoning above to the equation (3.11) with respect to v , using f (t, z + X t ) again with λ = -2 . We obtain first a solution p to (3.11) corresponding to f (t, z + X t ) and then derive that

w (t, z) = E[p (t, z -X t )]
is the unique bounded continuous (integral) solution w of the following problem:

       ∂ t w (t, z) = tr(Q(t)D 2 w (t, z)) + -2 (w (t, z + l(t)) -2w (t, z) + w (t, z -l(t))) + f (t, z), w (0, z) = 0. (3.13)
The previous estimates still hold with w instead of v , i.e., sup

(t,z)∈[0,T ]×R N |w (t, z)| ≤ T sup (t,z)∈[0,T ]×R N |f (t, z)|; (3.14) R(t) * D 2 w R(t) L p ((0,T )×R N ,m) ≤ C p f L p ((0,T )×R N ,m) . ( 3.15)
We would like now to let goes to zero, possibly passing to a subsequence n → 0, and prove that the associated limit w solves To do so we will proceed by compactness. Namely, we are going to prove that the family of solutions w solving (3.13), indexed by the parameter , is equi-Lipschitz on any compact subset of [0, T ] × R N and the same holds for any derivative in space of w . Indeed, one can apply the finite difference operators with respect to z at any order in (3.13). We recall that for a smooth function φ : R N → R, the first finite difference δ h,i φ, i ∈ {1, • • • , N } of step h > 0 in the direction e i (i th basis vector) is given by

   ∂ t w(t, z) = tr(Q(t)D 2 w(t, z)) + D 2 w(t, z) Q (t)e 1 , Q (t)e 1 + f (t, z), w(0, z) = 0 (3.
δ h,i φ(z) = φ(z + he i ) -φ(z) h , z ∈ R N .
For a given multi-index γ ∈ N N , the γ-th order finite difference operator δ h,γ , is then defined, for any h > 0, through composition. Namely,

δ h,γ φ(z) = δ γ 1 h,1 δ γ 2 h,2 . . . δ γ N h,N φ(z),
where δ γ i h,i denotes the γ i -th times composition of δ h,i with itself. Since any spatial derivative of f belongs to B b 0, T ; C ∞ 0 (R N ) , using (3.14) we deduce first that any finite difference of any order of w is bounded. Consequently, w is infinitely differentiable in space with bounded derivatives on [0, T ] × R N . Equation (3.13), to be understood in its integral form similarly to (2.2), then gives that those derivatives are themselves Lipschitz continuous in time (uniformly in the space variable). This precisely gives the equi-Lipschitz on any compact subset of [0, T ] × R N of the family w and of any spatial derivative of w . We can now apply the Arzelà-Ascoli theorem to w showing the existence of a sub-sequence {w n } n∈N which converges uniformly on any compact set to a function w : [0, T ] × R N → R. Similarly, any derivative in space of w n tends to the respective derivatives of w, uniformly on the compact sets. Passing to the limit as n → ∞ along the sequence ( n ) n in equation (3.13) (written in the integral form), we can then conclude that w solves (3.16).

Moreover, estimates (3.14) and (3.15) holds with w replaced by w. Iterating the previous argument in N steps we finally prove that the unique solution w to

   ∂ t w(t, z) = tr(Q(t)D 2 w(t, z)) + N k=1 D 2 w(t, z) Q (t)e k , Q (t)e k + f (t, z), w(0, z) = 0
(3.17) verifies estimates (3.14) and (3.15) with w replaced by w. The proof is complete.

Additional stability results in anisotropic Sobolev space and Schauder estimates

In this section we extend the previous approach to derive the stability with respect to a second order perturbation of the OU operator in (1.7) under the Kalman condition [K]. Here we consider also L p -estimates involving the degenerate components of the OU operator and some associated Schauder estimates.

Anisotropic Sobolev spaces and maximal L p regularity

With the notations of Section 2.

3 we write z ∈ R N as z = (x, y) = (x, y 1 , • • • , y k ) with x ∈ R d 0 , y i ∈ R d i , i ∈ {1, • • • , k}, k i=1 d i = d 1 .
Given β in (0, 1) and i in 1, k , we want to introduce the β-fractional Laplacian ∆ β y i along the component y i . To do so, we follow [START_REF] Huang | L p estimates for degenerate non-local Kolmogorov operators[END_REF] by considering the orthogonal projection p i : R N → R d i such that p i (z) = p i (x, y) = y i and denoting its adjoint by E i : R d i → R N . We can now define the β-fractional Laplacian ∆ β y i as:

∆ β y i φ(z) := p.v. R d i [φ(z + E i w) -φ(z)] dw |w| d i +2β , z ∈ R N ,
for any sufficiently regular function φ : R N → R. Let p in (1, +∞), we recall that we have denoted by L p ((0, T ) × R N ) the standard L p -space with respect to the Lebesgue measure. We can now define the appropriate anisotropic Sobolev space to state our results. For notational simplicity, let us denote

α i := 1 1 + 2i . (4.1) Set now α := (α 1 , • • • , α k ) ∈ R k . The homogeneous space Ẇ 2,p α ([0, T ] × R N ) is composed by all the functions ϕ : [0, T ] × R N → R in L p ([0, T ] × R N ) such that (t, z) ∈ [0, T ] × R N → ∆ x ϕ(t, z) ∈ L p ([0, T ] × R N )
, where ∆ x ϕ is intended in distributional sense, and for any i in 1, k , ∆ α i y i ϕ(t, z) is well defined for almost every (t, z) and

∆ α i y i ϕ(t, z) := ∆ α i y i ϕ(t, •)(z) belongs to L p ([0, T ] × R N ).
It is endowed with the natural semi-norm ϕ Ẇ 2,p α where

ϕ p Ẇ 2,p α = ∆ x ϕ p L p + k i=1 ∆ α i y i ϕ p L p . (4.2)
The thresholds in (4.1) might seem awkward at first sight. They actually correspond to the indexes needed to get stability of the harmonic functions associated with the principal part of (1.7), that is considering A 0 consisting in the subdiagonal part of A only (i.e., considering (2.6) when the diagonal elements and the strictly upper diagonal elements are equal to zero) along an associated dilation operator. Namely, setting

L ou 0 f (z) = Tr(BD 2 f (z)) + A 0 z, Df (z) , z = (x, y) ∈ R d 0 +d 1 = R N , ( 4.3)
so that A 0 , B satisfy [K], if (∂ t -L ou 0 )u(t, z) = 0 then for all λ > 0 (∂ t -L ou 0 )u δ λ (t, z) = 0 where the dilation operator

δ λ (t, z) = (λ 1/2 t, λx, λ 1/3 y 1 , • • • , λ 1/(1+2k) y k ).
precisely exhibits the exponents in (4.1) for the degenerate components. In [START_REF] Huang | L p estimates for degenerate non-local Kolmogorov operators[END_REF], see also [START_REF] Chen | Propagation of regularity in L p -spaces for Kolmogorov-type hypoelliptic operators[END_REF] and [START_REF] Menozzi | Martingale problems for some degenerate Kolmogorov equations[END_REF] where time inhomogeneous coefficients are considered as well, it has been proven that if A, B satisfy [K] and the diagonal and the strictly upper diagonal elements of A in (2.6) are equal to zero (i.e., A = A 0 ) then the following Sobolev estimates hold:

u Ẇ 2,p α ≤ C p f L p , ( 4.4)
with

C p = C p (ν, A, d 0 , d 1 )
, where again u is the unique bounded solution to the corresponding Cauchy problem (1.10). In particular we get also the maximal smoothing effects w.r.t. the degenerate directions. Note that the solution u to (1.1) verifies (4.4). The specific structure assumed on A is actually due to the fact that for such matrices there is an underlying homogeneous space structure which makes easier to establish maximal regularity estimates (see e.g. [START_REF] Coifman | Analyse Harmonique non-commutative sur certains espaces homogènes[END_REF] in this general setting).

If A, B satisfy [K] with a general A as in (2.6), having non zero strictly upper diagonal entries (non zero entries in the diagonal should not create difficulties) we believe that the approach in [START_REF] Bramanti | Global L p estimates for degenerate Ornstein-Uhlenbeck operators[END_REF] could extend to show that (4.4) still holds in this general setting. However such estimates have not been, up to our best knowledge, proven yet.

L p -estimates for the degenerate directions of special OU operators. Setting, as in Section 1.1, u(t, z) = v(t, e tA z) and since u solves (1.10) we have that v in turn solves (1.17). From the previous computations, setting

B I = I d 0 ,d 0 0 d 0 ,d 1 0 d 1 ,d 0 0 d 1 ,d 1
and considering A as in [START_REF] Huang | L p estimates for degenerate non-local Kolmogorov operators[END_REF], with the diagonal and the strictly upper diagonal elements of A equal to zero in (2.6), we derive

D 2 x u L p ((0,T )×R N ) = B I e tA * D 2 v(t, e tA •) e tA B I L p ((0,T )×R N ) ≤C p f (t, e tA •) L p ((0,T )×R N ) .
On the other hand, for all i ∈ {1, • • • , k} and with α i as in (4.1), we get from the definition (4.2) that the estimate (4.4) rewrites in term of v as:

∆ α i y i u p L p ((0,T )×R N ) = T 0 dt R N dz p.v. R d i [u(t, z + E i w) -u(t, z)] dw |w| d i +2α i p = T 0 dt R N dz p.v. R d i [v(t, e tA (z + E i w)) -v(t, e tA z)] dw |w| d i +2α i p = T 0 dt R N dz p.v. R d i [v(t, z + e tA E i w)) -v(t, z)] dw |w| d i +2α i p =: ∆ α i ,i,A v p L p ((0,T )×R N ) ,
v p Ẇ 2,p,A α := k i=0 ∆ α i ,i,A v p L p ((0,T )×R N ) ≤ Cp f p L p ((0,T )×R N ) (4.5)
with Cp = C p p . We now want to prove that for w solving (1.20), namely

       ∂ t w(t, z) = Tr e tA Be tA * D 2 w(t, z) + Tr e tA S(t)e tA * D 2 w(t, z) + f (t, z), (t, z) ∈ (0, T ) × R N , w(0, z) = 0, z ∈ R N , it also holds that w p Ẇ 2,p,A α := k i=0 ∆ α i ,i,A w p L p ((0,T )×R N ) ≤ Cp f p L p ((0,T )×R N ) , ( 4.6)
with the same constants Cp as in (4.5). This can be done through the previous perturbative approach of Section 3.2 employed to prove Theorem 2.2, which actually gives the expected control for the second order derivatives contribution of the semi-norm • Ẇ 2,p,A α . For the other contributions and with the notations of Section 3.2, with Q (s) = e sA S(s)e sA * and with m which is the Lebesgue measure on [0, T ] × R N (indeed in the present case g(t) = det(e -At ) = 1, for all t) we would get

v p Ẇ 2,p,A α = k i=0 ∆ α i ,i,A v p L p ((0,T )×R N ) = k i=0 (0,T )×R N |∆ α i ,i,A v (t, z)| p dzdt = k i=0 (0,T )×R N |E[∆ α i ,i,A v (t, z + X t )] | p dzdt ≤ k i=0 [0,T ]×R N E[|∆ α i ,i,A v (t, z + X t ) | p ]dzdt = k i=0 E [0,T ]×R N |∆ α i ,i,A v (t, z + X t ) | p dzdt = k i=0 E [0,T ]×R N |∆ α i ,i,A v (t, z) | p dzdt ≤ Cp f p L p ((0,T )×R N ) ,
using for the last inequality that v also satisfies (4.5) (similarly to what had been established in (3.7)). The same previous procedure and the final compactness argument then yields (4.6). Setting eventually ũ(t, z) := w(t, e tA z), which is the unique integral solution (smooth in space) of

   ∂ t u S (t, z) = L ou,S t u S (t, z) + f (t, z), (t, z) ∈ (0, T ) × R N , u S (0, z) = 0, z ∈ R N ,
where L ou,S t introduced in (1.13) is the Ornstein-Uhlenbeck operator perturbed at second order, we derive that

u S Ẇ 2,p α ≤ C p f L p , ( 4.7)
with C p as in (4.4). We have thus extended the results of Theorem 1.1 with the anisotropic Sobolev semi-norm in (4.2). The estimate (4.4) is stable for a continuous, non-negative second order perturbation of the underlying degenerate Ornstein-Uhlenbeck operator.

Anisotropic Schauder estimates

Following Krylov [START_REF] Krylov | Lectures on elliptic and parabolic equations in Hölder spaces[END_REF], for some fixed in N 0 := N ∪ {0} and β in (0, 1], we introduce for a function φ : R N → R the Zygmund-Hölder semi-norm as

[φ] C +β :=      sup |ϑ|= sup x =y |D ϑ φ(x)-D ϑ φ(y)| |x-y| β , if β = 1; sup |ϑ|= sup x =y D ϑ φ(x)+D ϑ φ(y)-2D ϑ φ( x+y 2 ) |x-y| , if β = 1
(we are using usual multi-indices ϑ for the partial derivatives). Consequently, the Zygmund-Hölder space C +β b (R N ) is the family of bounded functions φ : R N → R such that φ and its derivatives up to order are continuous and the norm

φ C +β b := i=0 sup |ϑ|=i D ϑ φ ∞ + [φ] C +β is finite.
We can now define the anisotropic Zygmund-Hölder spaces associated with the current setting and which again reflect the various scales already introduced in (4.1). Let γ ∈ (0, 3), the space C γ b,d (R N ) is the family of functions φ : R N → R such that for any i in 0, k and any z 0 in R N , the real function

w ∈ R d i → φ(z 0 + E i (w)) belongs to C γ/(1+2i) b R d i ,
with a norm bounded by a constant independent from z 0 . In the above expression, we recall that the (E i ) i∈{1,••• ,k} have been defined in the previous paragraph, d 0 = d 0 and E 0 is the embedding matrix from with the same constant C β as in (4.11). We proceed through the previous perturbative approach of Section 3.2. With the notations employed therein, we deduce that there exists a unique solution v = PDE(Q, f (t, z -X t )), depending also on and ω as parameters such that sup By the translation invariance of the Hölder-norms, using also that X t = e tA e -tA X t , it is not difficult to prove that, for any ω, P-a.s., ].

R d 0 into R N . It is endowed with the norm φ C γ b,d := sup z 0 ∈R N φ z 0 + E 0 (•) C γ b (R d 0 ) + k i=1 sup z 0 ∈R N φ z 0 + E i (•) C γ/(1+2i) b (R d i ) . ( 4
We would get, similarly, using (4.15) for the last inequality. Now, continuing as in Section 3.2, using also a compactness argument, one would derive that (4.12) indeed holds. Going backwards, setting ũ(t, z) := w(t, e tA z), we find that ũ is the unique (integral) solution u S to (1.14); we finally derive that where C β is the same constant as in (4.9). Estimate (4.17) provides the extension of Theorem 1.1 for the anisotropic Schauder estimates.

Remark 4.1. Let us mention that for the perturbative method to work, roughly speaking, few properties were actually needed on the underlying norm. Namely, we used the translation invariance and some kind of commutation between the norm (or a function of the norm in the L p -case) and expectation. Hence, this approach could possibly be applied to a much wider class of estimates in other function spaces (like e.g. Besov spaces). This will concern further research.

  [START_REF] Lanconelli | On a class of hypoelliptic evolution operators[END_REF]) and estimates (3.14) and (3.15) hold with w replaced by w.

  using that T r(A) = 0. Hence, setting ∆ α 0 ,0,A v p L p ((0,T )×R N ) := Tr B I e tA * D 2 v(t, e tA •) e tA B I p L p ((0,T )×R N ) = Tr B I e tA * D 2 v e tA B I p L p ((0,T )×R N ) ,

. 8 ) 10 )

 810 We denote by C γ b,d this function space because the regularity exponents reflect again the multi-scale features of the system; the norm could equivalently be defined through the corresponding spatial parabolic distance d defined as follows. For all z = (x, y), z = (x , y ) ∈ R N :d(z, z ) := |x -x | + k i=1 |y i -y i | 1 1+2i ,where the exponents are again those who appeared in (4.1).Let f ∈ B b 0, T ; C ∞ 0 (R N ) . Under [K],by the results of Lunardi[START_REF] Lunardi | Schauder estimates for a class of degenerate elliptic and parabolic operators with unbounded coefficients in R n[END_REF] it follows that the unique bounded solution of the Cauchy Problem (1.10) (written in integral form) verifies the following anisotropic Schauder estimatesu L ∞ ((0,T ),C 2+β b,d ) ≤ C β f L ∞ ((0,T ),C β b,d ) ,(4.9)for some constant C β independent from f , i.e.,We again set as in the previous paragraph u(t, z) = v(t, e tA z)u L ∞ ((0,T ),C 2+β b,d ) = v(t, e tA •) L ∞ ((0,T ),C 2+β b,d ) =: v L ∞ ((0,T ),C 2+β b,d,A ) ≤ C β f L ∞ ((0,T ),C β b,d ) = C β f (t, e tA •) L ∞ ((0,T ),C β b,d ) =: C β f L ∞ ((0,T ),C β b,d,A ) , (4.11) denoting f (t, z) := f (t, e -tA z). We again want to prove as in Section 1.1 that for w solving (1.20), w L ∞ ((0,T ),C 2+β b,d,A ) ≤ C β f L ∞ ((0,T ),C β b,d,A ) (4.12)

(

  t,z)∈[0,T ]×R N |v (t, z)| ≤ T sup (t,z)∈[0,T ]×R N | f (t, z)|. (4.13)

fL

  ∞ ((0,T ),C β b,d,A ) = f (•, • -X • ) L ∞ ((0,T ),C β b,d,A ) .

(4. 14 )

 14 Thus it also holds from (4.11)v L ∞ ((0,T ),C 2+β b,d,A ) ≤ C β f L ∞ ((0,T ),C β b,d,A ) .

(4. 15 )

 15 Recalling now that v (s, z) = E[v (s, z + X s )] is an integral solution of ∂ t v (t, z) = tr(Q(t)D

[D 2

 2 x v (t, e At (z 0 + E 0 (•)))] C β ≤E[[D 2 x v (t, e At (z 0 + E 0 (•)))] C β ],and for all k ∈ {1, 2},D k x v (t, e At (z 0 + E 0 (•))) ∞ ≤ E[ D k x v (t, e At (z 0 + E 0 (•))) ∞ ].Summing those contributions, we thus derive from (4.8), (4.11) that:v L ∞ ((0,T ),C 2+β b,d,A ) ≤ sup 0≤t≤T E[ v (t, •) C 2+β b,d,A ) ] ≤ C β f L ∞ ((0,T ),C β b,d,A) , (4.16)

u

  S L ∞ ((0,T ),C 2+β b,d ) ≤ C β f L ∞ ((0,T ),C β b,d ) ,(4.17

  )

t, e tA z) + Dv(t, e tA z), Ae tA z -Tr e tA Be tA * D 2 v(t, e tA z) -Dv(t, e tA z), Ae tA z

  = v t (t, e tA z)-Tr e tA Be tA * D 2 v(t, e tA z) .

	(1.16)
	Denoting f (t, z) := f (t, e -tA z), it now follows that v satisfies the PDE:

  .10) Namely we have to prove the following result. Let us consider equations (2.1) and (2.10) where Q(t), Q (t) are two continuous in time, non-negative definite matrices in R N

	Theorem 2.2.

  At (z 0 + E i (w )) + e At e -At X t )|] ≤E[[v (t, e At (z 0 + E i (•)))] At (z 0 + E i (•)))] At (z 0 + E i (•)))]

			C	2+β 1+2i	] |w -w |	2+β 1+2i .
	Hence,
	[v (t, e C	2+β 1+2i	≤E[[v (t, e C	2+β 1+2i

2 z v (t, z)) + λ (v (t, z + l(t)) -v (t, z)) + f (t, z),

with zero initial condition, we write that for

i ∈ {1, • • • , k}, w, w ∈ R d i , (t, z 0 ) ∈ [0, T ] × R N , |v (t, e At (z 0 + E i (w)) -v (t, e At (z 0 + E i (w ))| ≤E[|v (t, e At (z 0 + E i (w)) + e At e -At X t ) -v (t, e

Acknowledgment

For the first author, the work was supported by a public grant (2018 -0024H) as part of the FMJH project. For the second author, the article was prepared within the framework of the HSE University Basic Research Program. The third author has been partially supported by GNAMPA of the Istituto Nazionale di Alta Matematica (INdAM).


Abstract-Genome-Wide Association Studies (GWAS) identify the genomic variations that are statistically associated with a particular phenotype (e.g., a disease). GWAS results, i.e., statistics, benefit research and personalized medicine. The confidence in GWAS increases with the number of genomes analyzed, which encourages federated computations where biocenters periodically include newly sequenced genomes. However, for legal and economical reasons, this collaboration can only happen if a release of GWAS results never jeopardizes the genomic privacy of data donors, if biocenters retain ownership and cannot learn each others' data. Therefore, two challenges need to be simultaneously addressed to enable federated GWAS: (i) the computation of GWAS statistics must rely on secure and privacy-preserving methods; and (ii) GWAS results that are publicly released should not allow any form of privacy attack. In this talk, we will introduce the problem we consider in more details and present DyPS [START_REF] Pascoal | Dyps: Dynamic, private and secure gwas[END_REF], the framework we have designed and recently presented at the Privacy Enhancing Technologies Symposium (PETS). We refer the reader to the full paper 1 for the details we cannot cover in this short version.

I. INTRODUCTION

GWAS allow the discovery of correlations between genetic variations, such as single-nucleotide polymorphisms (SNPs), and a phenotypic trait. For example, GWAS are widely used for the early detection of disease susceptibility. The GWAS results that we consider in this work include: (i) aggregate statistics, such as minor allele frequencies (MAF), single and pairwise allele frequencies; and (ii) test statistics, such as linkage disequilibrium (LD), χ 2 and p-values.

To obtain meaningful statistical findings, GWAS require large datasets, which motivates the development of federated environments where several biocenters can securely and collaboratively compute GWAS statistics over their genomic data. Nevertheless, enforcing secure distributed GWAS computation alone is not sufficient since the publication of GWAS results might enable sophisticated membership [START_REF] Homer | Resolving individuals contributing trace amounts of dna to highly complex mixtures using high-density snp genotyping microarrays[END_REF]- [START_REF] Humbert | De-anonymizing genomic databases using phenotypic traits[END_REF] or recovery attacks [START_REF] Zhou | To release or not to release: Evaluating information leaks in aggregate human-genome data[END_REF]- [START_REF] Ayday | Inference attacks against kin genomic privacy[END_REF]. Recent works have focused on either secure and privacy-preserving GWAS computation, or on privacy-preserving GWAS data releases, but not on both. However, we argue that those two problems need to be simultaneously addressed, since both can produce information leakage and lead to privacy attacks.

In addition, we observe that attacks can be facilitated by collusions among dataholders. Indeed, a subset of biocenters might be able to collude to infer other biocenters' genomic data. Furthermore, there is now a need to enable dynamic GWAS where new genomes can be added or removed (in order to comply with privacy regulations, such as GDPR) over time. If not crafted with proper care, GWAS updates 1 https://petsymposium.org/2021/files/papers/issue2/popets-2021-0025.pdf We introduce DyPS [START_REF] Pascoal | Dyps: Dynamic, private and secure gwas[END_REF], a novel framework that provides both secure GWAS computation and privacy-preserving public releasing of GWAS results. DyPS leverages a Trusted Execution Environment to secure dynamic GWAS computations. Moreover, DyPS uses a scaling mechanism to speed up the releases of GWAS results according to the evolving number of genomes used in the study, even if individuals retract their participation consent. Lastly, DyPS also tolerates up to allbut-one biocenters colluding without privacy leaks.

II. PRIVACY-PRESERVING GWAS COMPUTATION

DyPS leverages a TEE-enabled centralized server to allow the safe sharing of genomic data from several biocenters, and therefore allows individuals to stay in control of their own genomes, which are never revealed to other biocenters. DyPS's computations are performed inside an enclave, which is an processing environment that cannot be tampered with. Even privileged code (e.g., the hosting operating system) cannot access its data. DyPS uses an Intel SGX enclave that can be attested by biocenters to certify that trusted code is being run by a genuine SGX enclave. Once the remote attestation step succeeds, biocenters can securely outsource their genome data through authenticated secure channels, after encrypting and signing it. Moreover, it is also possible to store data outside the enclave (which only has 96 MB of available memory) by using a sealing mechanism. Such data is encrypted by the enclave and can only be accessed and verified by this enclave.

We illustrate DyPS's federated architecture in Figure 1. We assume the presence of B honest-but-curious biocenters, among which f = B -1 might collude. DyPS's algorithms only run inside the enclave, which enforces isolation, confidentiality, and integrity of both the data and operations. The enclave periodically collects the requests from the various biocenters and decides which requests are to be executed to dynamically update and release the GWAS results.

III. PROTECTING DYNAMIC GWAS RELEASES

DyPS mainly builds on Zhou et al.'s work [START_REF] Zhou | To release or not to release: Evaluating information leaks in aggregate human-genome data[END_REF] and on SecureGenome [START_REF] Sankararaman | Genomic privacy and limits of individual detection in a pool[END_REF], which define conditions for identifying when the release of static GWAS results is safe. In particular, DyPS extends such approaches in order to support dynamic updates of GWAS results.

To protect GWAS aggregate statistics against recovery attacks and test statistics against recovery and membership attacks, DyPS makes sure that enough genomes are being used so that any probabilistic polynomial-time (p.p.t.) adversary observing released results cannot successfully recover the complete solution space, and determine the right genotype sequences within the human-genome data set (i.e., individuals' real genome sequence) that participated in a given release.

In addition, to protect aggregate statistics against membership attacks, DyPS uses the Likelihood Ratio Test (LR-Test) to certify that no genomes participating in a release can be identified. For that purpose, DyPS leverages Se-cureGenome [START_REF] Sankararaman | Genomic privacy and limits of individual detection in a pool[END_REF], which is an LR-test approach that defines over which SNP positions can allele frequencies be safely released within a cohort of genomes.

In order to allow safe dynamic releases of GWAS, DyPS adopts the following pipeline:

(1) FIFO request pool. DyPS collects and treats genome addition and removal requests from biocenters in FIFO order.

(2) Requests selection. DyPS runs an algorithm that selects a safe batch of genome requests so that any possible subset of honest biocenters collectively have enough operations. This way, an external adversary or the remaining possibly colluding biocenters cannot isolate and attack a small number of requests.

(3) GWAS processing. DyPS computes and releases test statistics over the selected genomes.

(4) Membership tests. DyPS identifies over which subset of SNPs can aggregate statistics also be released or updated.

IV. PERFORMANCE EVALUATION

We used Graphene to run DyPS inside an SGX enclave. Biocenters adopted AES256 encryption and ECDSA signatures to outsource genome data. Genome addition or removal requests are generated following a Poisson distribution, with λ = 8 for additions and λ = 6 for removals. For larger-scale experiments, we proportionally increased λ. We evaluated DyPS under different GWAS scenarios. For instance, we considered a number of biocenters that ranges from 4 to 7, with up to f = B -1 of them colluding. We used simulated genomes (up to 6 million) and real genomes (up to 27,895 genomes).

First, we show that DyPS request selection mechanism runs with O(1) complexity while a Brute Force approach would run in O(2 N ). In addition, DyPS requires ≈ 75KB to encode 300,000 SNPs, and the overall request size is 258 bits, i.e., ≈ 48 Bytes after encryption. DyPS also presents a reasonable memory consumption (2 MB on average, which is within SGX's limits). DyPS also remains practical under a large-scale setting (5 biocenters, 5,000 SNPs, and approximately 28,000 real genomes), with a reasonable running time, even though it would increase with the number of releases due to the membership evaluation.

DyPS maintains genomic privacy, but releases GWAS results less frequently than a naive approach, which in turn would leave 4.98% of its releases vulnerable to attacks and 8% of the genomes exposed to privacy threats. In addition, the existence of colluding players only impacts the speed with which requests are selected and processed. As expected, the larger the number of colluding players the larger are the addition and removal request delays. Last but not least, we show that DyPS can regularly release GWAS results thanks to our dynamic scaling methods.

V. CONCLUSION

This paper shows the challenges and highlights the need to reconcile secure GWAS computation with privacypreserving GWAS releases. On one hand, an adversary might launch privacy attacks during the computation of GWAS. On the other hand, GWAS releases must preserve privacy, which requires computations over the distributed genomic dataset. Unfortunately, current works in the literature are not able to combine both aspects of a fully privacy-preserving GWAS (processing and releasing) in a homogenized form. Moreover, given the reduced cost of sequencing DNA nowadays, there is now a need to update GWAS results in a dynamic manner, while enabling donors to withdraw consent at any time. In this regard, we presented DyPS, a TEE-based framework that enables privacy-preserving processing and releasing of federated GWAS. DyPS is the first work capable of reconciling both challenges in a secure, scalable, and reliable fashion, and thus of protecting genomic privacy in terms of both biocenters' and donors' genome data. Future works include considering the existence of byzantine biocenters, and an oblivious version of DyPS to cope with side-channel attacks against SGX enclaves.

  DyPS: Dynamic, Private and Secure GWAS Talk proposal -GenoPri 2021 Túlio Pascoal * , Jérémie Decouchant † , Antoine Boutet ‡ , Paulo Esteves-Verissimo § * University of Luxembourg, † Delft University of Technology, ‡ Univ. Lyon, INSA Lyon, Inria, CITI, § KAUST -RC3

FederationFig. 1 :

 1 Fig. 1: DyPS' architecture to enable secure and privacypreserving federated GWAS.


Behzad Mirmahboub

Jérôme Moré

David Youssefi

Alain Giros

François Merciol

Sébastien Lefèvre

Fast Pattern Spectra using Tree Representation of the Image for Patch Retrieval

Keywords: content based image retrieval, patch retrieval, tree representation, pattern spectra, large satellite images

We extend the notion of content based image retrieval to patch retrieval where the goal is to find the similar patches to a query patch in a large image. Naive searching for similar patches by sequentially computing and comparing descriptors of sliding windows takes a lot of time in a large image. We propose a novel method to compute descriptors for all sliding windows independent from number of patches. We rely on tree representation of the image and exploit the histogram nature of pattern spectra to compute all the required descriptors in parallel. Computation time of the proposed method depends only on the number of tree nodes and is free from query selection. Experimental results show the effectiveness of the proposed method to reduce the computation time and its potential for object detection in large images.

Introduction

Content based image retrieval (CBIR) is the problem of finding images in a database that are similar to a query image [START_REF] Liu | A survey of content-based image retrieval with high-level semantics[END_REF]. This is generally done by two main components. The first step is feature extraction that computes discriminating descriptors for the images. Therefore, each image is represented by a descriptor. The second step is to find the similarities between the images by calculating the distances between their descriptors. Then, the dataset images are ranked based on their distances to the query image and the most similar images are retrieved.

Current datasets for image retrieval consist of the cropped images with specific content for each image and the task of CBIR is to find the similar images. However, in some applications such as remote sensing we have a patch in a large satellite image and we need to find similar patches in the same image or other images (Figure 1). CBIR can be trivially adapted to this problem by dividing the large image into small patches and extracting a descriptor for each image patch. Then, distances between patches can be computed to find the most similar patches to the query patch. This approach is time consuming and is not applicable in real time retrieval system for very large images.

In order to address this problem, we propose in this paper a patch retrieval system to quickly compute descriptors for all the patches in an image. As far as we know, it is the first work that tries to find the similar patches in a large image. We rely on Pattern Spectra (PS) based on tree representation of the image [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF]. It is a histogram-like morphological descriptor for images. We introduce fast-PS descriptor as an approximation of PS that is computed for all image patches in parallel. Therefore, the computation time will be independent from the number of patches which is important for processing large images. Fig. 1: Patch retrieval: User selects a query patch from an image and the goal is to find similar patches in the same or other images. This paper is organized as follows. Section 2 reviews previous works on image retrieval, tree representation of the image and pattern spectra. We explain our proposed method in Section 4 and present experimental results in Section 5. Section 6 concludes the paper.

Related works

Image Retrieval

CBIR is a challenging task to find the images in a database that are similar to a query image [START_REF] Liu | A survey of content-based image retrieval with high-level semantics[END_REF]. The general framework consists of two main steps as it is shown in Figure 2. First, a discriminant feature vector is extracted from each image. Second, these descriptors are used to compute the distance of each dataset image to the query image. The dataset images are sorted according to their distances to the query image and we expect to see the similar images at the beginning of the ranked list.

Various types of feature vectors are proposed in literature based on color or texture of the image such as Gabor Filter [START_REF] Li | Gabor-filtering-based nearest regularized subspace for hyperspectral image classification[END_REF], LBP (Local Binary Pattern) [START_REF] Li | Local binary patterns and extreme learning machine for hyperspectral imagery classification[END_REF], SIFT (Scale-Invariant Feature Transform) [START_REF] Paul | Remote sensing optical image registration using modified uniform robust sift[END_REF], MSER (Maximally Stable Extremal Regions) [START_REF] Yu | Mser based shadow detection in high resolution remote sensing image[END_REF] and recently deep CNN features [START_REF] Chen | Deep feature extraction and classification of hyperspectral images based on convolutional neural networks[END_REF]. An ideal feature vector must be rotation/translation/scale-invariant and robust to changes in illumination. Descriptor computation is the bottleneck in image retrieval. Therefore, Fig. 2: General framework of content based image retrieval. The images in the database are sorted based on their distances to a query image. (The images are selected from Merced dataset [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF] only for illustration. ) practical image retrieval systems pre-compute image descriptors and store them along the original images.

Distance between two feature vectors can be simply computed using their Euclidean distance. However, in many cases the feature vectors are high dimensional and different dimensions have different importance. An efficient solution is to learn a weighting matrix to emphasize the important dimensions [START_REF] Koestinger | Large scale metric learning from equivalence constraints[END_REF]. On the other hand, some distance metrics are specially designed to compare histograms. Since pattern spectra is a histogram-like feature vector, we use the well-known χ 2 distance to compare them. The χ 2 distance of two histograms x and y tries to reduce the effect of large bins:

X 2 px, yq " ÿ i px i ´yi q 2 px i `yi q (1)
Previous works on image retrieval mainly try to compare two whole images [START_REF] Cheng | Remote sensing image scene classification meets deep learning: Challenges, methods, benchmarks, and opportunities[END_REF]. To the best of our knowledge, it is the first time that the concept of CBIR is extended to find similar patches in a large image as depicted in Figure 1. We assume that the user selects a query patch with arbitrary size. Therefore, the query patch size is not known and the descriptors cannot be computed in advance. A closely related field of study exists in literature for patch matching [START_REF] Barnes | Patchmatch: A randomized correspondence algorithm for structural image editing[END_REF]. However, the solutions to this problem cannot meet our requirements since patch matching is generally applied between two similar scenes such as stereo images. Therefore, those methods exploit similarity between images to restrict the possible locations of matching patches and expedite the search, while we do not impose that the images have similar content and therefore we cannot use such methods.

Evaluation metric

When we calculate and sort the distances between the image patches and the query patch, ideally we expect to see all the relevant images (true positives or tp) at the beginning of the retrieved list. However, since the descriptors and distance metric are not perfect, we also see irrelevant images (false positives or f p) in the list and even some of the relevant images are not retrieved (false negatives or f n). Mean Average Precision (mAP) is a measure to evaluate the performance of a retrieval system [START_REF] Manning | Introduction to information retrieval[END_REF]. It shows how well the relevant images are ranked in the retrieved list and it is calculated based on precision and recall.

Precision is defined as the ratio of relevant retrieved images to all retrieved images or tp{ptp `f pq and it is computed for each rank position in the retrieved list. In the example of Figure 2 if baseball field is the relevant class, the precisions for first seven ranks will be t 1 1 , 1 2 , 2 3 , 2 4 , 2 5 , 3 6 , 3 7 u. Recall is defined as the ratio of retrieved relevant images to all relevant images in the dataset or tp{ptp `f nq. Assume that there are four relevant baseball fields in Figure 2 and that three of them are retrieved in the first seven ranks. Then, the recalls for these ranks will be t 1 4 , 1 4 , 2 4 , 2 4 , 2 4 , 3 4 , 3 4 u. Average Precision (AP) is the area under precision-recall curve. Practically, it is the average of precisions at the positions of relevant retrieved images (where the recall is changed). In the above example AP is calculated as 1 3 p 1 1 `2 3 `3 6 q " 13 18 . Finally, mAP is simply the mean value of all APs for different queries.

3 Theoretical background

Tree Representation

Usual image representation by a matrix of pixels gray levels does not include spatial relations between image components. In contrast, hierarchical image representation shows an image with a tree structure that considers spatial relations [START_REF] Bosilj | Partition and inclusion hierarchies of images: A comprehensive survey[END_REF]. Each node in this tree structure represents a connected component in the original image. Each pair of components in the image are either disjoint or one of them is a subset of the other one making a child-parent relationship. The leaves of the tree consist of the smallest components while the tree root represents the entire image. Such a hierarchical structure lets us quickly process the image in different scales. Depending on how to define the image components, various tree types can be built including max-tree, min-tree, Tree-of-Shapes, α-tree, ω-tree and Binary Partition Tree [START_REF] Bosilj | Partition and inclusion hierarchies of images: A comprehensive survey[END_REF].

Figure 3 shows a simple gray scale image on the left with its max-tree representation in middle. Max-tree is built using the connected components of upper level sets of the image [START_REF] Jones | Component trees for image filtering and segmentation[END_REF] that is defined as:

H max " tCpL k q|v min ď k ď v max u ( 2
)
where k is a threshold that is selected between minimum gray level v min and maximum gray level v max . Level set L k is an image that is obtained by thresholding the original image with k. It only contains all the components with intensity values equal or greater than k and CpL k q is the set of those components. Therefore, H max is the set of all connected components in different levels.

Each circle in Figure 3 represents a node in the tree structure and the corresponding connected component in the image is shown beside it. Also the gray level k is written next to each tree node. The leaves of the max-tree consists of the brightest image components or local maxima while the tree root represents the whole image. Fig. 3: Computation of pattern spectra using area attribute of a max-tree nodes. Tree nodes are shown by circles in the tree structure. The number inside each circle shows the node area and the number next to it represents the node altitude. Each tree node contributes to a histogram bin by its area multiplied by the gray level difference with its parent.

Pattern Spectra

Tree representation is an efficient structure to process the image at multliple scales. Pattern Spectra is a morphological image descriptor that can be computed using the tree representation [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF]. It is a histogram-like feature vector that shows the distribution of components in an image. After building the tree structure, various geometrical or statistical attributes can be computed for each tree node and PS shows the frequencies of those attributes in the image. In the example of Figure 3, area attribute for each node is written in the relevant circle in the tree structure and PS is shown on the right of the figure. In this example, PS consists of five bins for the area attribute and each tree node contributes to one relevant bin. Each node contributes by its volume which is defined as the node area multiplied by the difference between node gray level and its parent gray level. Contributions of all nodes to the relevant PS bin are summed and the final descriptor will be r8, 14, 0, 0, 48s. Generally, bin j of pattern spectra P S is computed based on tree representation of the image as:

P S j " ÿ @i,ApciqPrbj´1,bj r apc i q ˇˇkpc i q ´k`p pc i q ˘ˇ( 3)
where apc i q, kpc i q and ppc i q represent the area, the gray level and the parent of the component c i respectively. Therefore, the area of each component is weighted by the grey level difference with its parent. Apc i q denotes an attribute of the component c i . All the components that their Apc i q attributes are located in an interval between b j´1 and b j contribute to bin j of the PS. The histogram bin edges tb 0 , ..., b max u of the desired attribute can be divided linearly or logarithmic.

In the above example, the attribute Apc i q is defined as the area of the component c i and it is distributed linearly between 0 and 50 producing a size pattern spectra. However, it can be any other type of attribute such as compactness value that produces shape pattern spectra.

Pattern spectra describes the probability of presence of a component with a certain size or shape in the image. It gives an estimation for the amount of details that is removed from the image after filtering its components. PS is a translation and scale invariant image descriptor that can be used for image classification and retrieval. PS is also rotation invariant if it is computed based on rotation invariant attributes. Several node attributes may be used to compute a multidimensional PS. The number of bins for each attribute is a parameter that should be selected for each dimension. Two dimensional PS based on elongation and entropy was used in [START_REF] Bosilj | Retrieval of remote sensing images with pattern spectra descriptors[END_REF] for satellite image retrieval. The authors also proposed local pattern spectra that computes PS in overlapped patches in different scales of the image. The final image descriptor is the aggregation of all computed PS.

Proposed Method

In conventional image retrieval, similar images to a query image are selected according to the distances between their descriptors. In this paper, we aim to extend that concept to patch retrieval where we look for some patches in the dataset that are similar to a query patch. The similar patches can be found in the same image of the query patch or in other images as shown in the example of Figure 1.

The naive method is to divide each image to a regular grid of patches with the same size of the query patch. Then, a sliding window is used and the similarity of each window with the query patch is computed. The computation time of this approach is proportional to the number of sliding windows that is too much for a large image. Notice that the size (and so the number) of the sliding windows depends on the user selected query patch. Therefore, it is not possible to compute and store the descriptors of the sliding windows in advance. In this section, we exploit the tree representation of the image and introduce the notion of fast pattern spectra descriptor to reduce the searching time. Finally the computation time is related to the number of tree nodes instead of number of the sliding windows.

Sub-Pattern Spectra

As explained in Section 3.1, pattern spectra is a histogram-like image descriptor that represents the distribution of connected components in the image. Each image component contributes to one of the histogram bins and all image components build the whole histogram. Since a small patch in the image contains a subset of image components, we expect that it contributes to only a few number of histogram bins. This intuition is shown in Figure 4 where all the image components make the full histogram and the components in the small patch make a sparse histogram called here "sub-pattern spectra". We also notice that each image patch consists of a subset of nodes in the tree structure of the whole image and this subset of nodes does not necessarily form a single integrated sub-tree structure. Having all related sub-pattern spectra computed, we can calculate the distance between the query patch and the sliding windows in the image to find the similar windows as depicted in Figure 5. The straightforward approach to compute sub-pattern spectra for each image patch is to directly compute pattern spectra on each window that are called here "local pattern spectra". Two issues arise from this approach. First, the components in the patch are not the same components in the whole image, because some of the image components are partially located in the sliding window and create new local components (Figure 6 Top). Second, there are many sliding windows in a large image and directly computing a pattern spectra for each of them will be time consuming. We address these problems by introducing "fast pattern spectra" in next section. Fig. 5: Similarity of a query patch with a sliding window is computed based on the distance between their sub-pattern spectra.

Fast Pattern Spectra

Our motivation in this paper is to avoid directly computing a descriptor for each image patch. For this purpose we rely on the histogram nature of pattern spectra and propose an indirect approach to compute the histogram in each patch. In this case, sub-pattern spectra consists of a subset of tree nodes in the whole image. Each node in the tree structure comes with some information such as center and bounding box of the corresponding component in the image. Therefore, when we select a patch from the image, we know which tree nodes are located inside it. We use these nodes to create our sub-pattern spectra and refer to it as fast pattern spectra (Figure 6 Bottom). Fig. 6: Sub-pattern spectra of a given patch can be computed in two different ways. (Top) local pattern spectra: A small tree is built directly on the patch and its pattern spectra is computed. The nodes of this small tree may be different from the big tree of the whole image. (Bottom) fast pattern spectra: Tree representation of the whole image is constructed and only the nodes that are located inside the patch are selected to compute the histogram.

The proposed fast pattern spectra is an estimation of local pattern spectra. Instead of computing PS on the patch separately, we use the tree representation of the whole image to compute descriptors for all image patches. The computation time of the proposed descriptor is independent from the number of the patches that is determined by user selection. Instead, the processing time depends on the number of tree nodes that are fixed for each image.

Figure 7 shows the process of computing fast pattern spectra. We build a tree representation of a given image and compute the desired attribute for its nodes. According to Equation 3, pattern spectra is the sum of components volumes which are computed by multiplying the nodes areas by their gray level differences with their parents vpc i q " apc i q ˇˇkpc i q ´k`p pc i q ˘ˇ. We calculate the volumes of all tree nodes and store them in a "tree table" along with nodes locations in the original image. Then, we compute a full pattern spectra and use it to assign a bin number to each tree node. The bin number is also stored in the tree table. When a user selects a query patch, we divide the image into sliding windows based on that selection and create an empty histogram for each sliding window. Therefore, we have a three dimensional matrix P Srx, y, bs spanned by image locations and histogram bins. Each tree node c i contributes by its volume vpc i q to the histogram bin bpc i q in a specific location pxpc i q, ypc i qq. In order to compute fast pattern spectra, we traverse the tree table and simply add the volume of each node to the corresponding histogram bin. It is formulated in Equation 4that is computed with complexity Opnq where n is number of the tree nodes. In this way, we get the sub-pattern spectra for all sliding windows of the image in one pass of the tree structure. We point out that the tree table is pre-computed and fixed. Every time that the user selects a new query patch, only Equation 4 needs to be evaluated again. P Srpxpc i q, ypc i q, bpc i qs += vpc i q , i " 1, ..., n

Experimental Results

In this section, we experimented our proposed fast pattern spectra for the application of patch retrieval on selected images and compare it with baseline local pattern spectra. All the experiments were done using Python with Higra package [START_REF] Perret | Higra: Hierarchical graph analysis[END_REF] in Windows 7 64-bit Operating System running on a laptop with Intel ® Core™ i7 CPU @ 2 GHz and 8 GB RAM.

Dataset

The ideal dataset for patch retrieval, as we addressed in this paper, consists of a large aerial image with annotations for various landmarks on it. The available datasets for content based image retrieval such as "UC Merced Land Use Dataset" [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF] are not suitable for this purpose since they consist of separate image patches with different classes. As far as we know, there is no available dataset with the desired annotations. Therefore, in order to conduct the reported experiments, we selected two satellite RGB images and manually annotated them with 6 classes including farm, forest, lake, airport, road and building as shown in Figure 8. We used Figure 8(a) as query image where the user selects its query patches. Three user selected patches with different sizes and contents are marked in this image. Figure 8(c) was used to search for the similar patches. Each image patch may contain various pixel labels. We assigned the majority of pixel labels inside each patch as its class.

Experiment Settings

We compare the proposed fast pattern (fast-PS) spectra with the baseline local pattern spectra (local-PS). In both cases, the user selects an arbitrary patch from the query image. Then, the search image is divided to regular overlapped windows with the same size of the query patch and the stride of half size of the query patch horizontally and vertically. In the case of baseline local-PS, pattern spectra is computed independently on each search window. However, in the proposed fast-PS approach, a big tree is built on the whole search image and the pattern spectra is computed for all the search windows in parallel as detailed in Section 4.2.

In these experiments, pattern spectra is computed based on the max-tree structure that is built on the gray scale image. We use one dimensional pattern spectra based on Apc i q " compactness attribute. The compactness of a tree node is defined as its area divided by the square of its perimeter. We divided nodes compactness linearly into 30 bins resulting to 30-dimensional descriptors for image patches. The similarity between descriptors are computed using χ 2 distance.

Results and Discussions

We selected three query patches with different sizes and contents as marked in Figure 8(a) and retrieved similar patches from Figure 8(c). For each case, we computed the average precision as explained in Section 2.2 using baseline and proposed methods. The average precision and computation time are reported in Table 1. The second column of table shows the number of sliding windows in the search image that depends on the size of the selected query patch. As we expected, local-PS needs more computation time which depends on number of sliding windows. However, the computation time of fast-PS is almost constant since it depends only on the number of tree nodes. As we explained in Section 4.2, in proposed fast-PS we need to store a table of tree nodes attributes for search image. Last column of Table 1 reports the required time to compute this tree table. Notice that the tree table is independent from query patch. It is computed only one time for each search image and is stored along the image. The stored tree table can be used later for retrieval. Therefore, its computation time does not affect the retrieval time. However, it adds an extra cost for the storage. In our experiment the search image is 1.09 MB and the stored tree table is 4.42 MB. We basically proposed fast-PS (as its name suggests) as an approximation of local pattern spectra that is computed quickly over sliding windows in a large image. Table 1 shows that this novel descriptor actually works better than the original local-PS for the application of patch retrieval. A possible reason can be explained referring to Figure 6. When we crop a patch from a large image, it contains additional components that do not exist in the original image and may be irrelevant to its contents. In this case, fast-PS helps to rely on relevant components from original image and achieves better results. Figure 9 shows the precision-recall curves for three selected query patches. We observe that fast-PS performs better than local-PS especially at first ranks.

Conclusion and Future Lines

We extended the problem of content-based image retrieval and introduced a new application that we call patch retrieval where a user selects a query patch with arbitrary size from a large image and we look for similar patches in the same or other image. Available solutions from conventional image retrieval need to divide the image into many sliding windows and compute an image descriptor on each window separately. This approach takes a lot of time especially for large images. We relied on pattern spectra that is an image descriptor based on tree representation of the image. We exploited the histogram nature of pattern spectra and proposed fast pattern spectra. It is computed for all sliding windows of the image in parallel and its processing time is independent from the number of the sliding windows. Interestingly, it is not only faster than local pattern spectra, but also achieves a higher average precision. However, the proposed method comes with a cost that needs to store an extra structure along with the original image. We showed that the fast pattern spectra is very promising on a simple dataset. More experiments on larger datasets are necessary to explore its full potential.

Our main goal in this work was to show the capability of fast-PS to offer interactive search in a large image when the size of the query patch is not known in advance. Although we found out that the average precision of fast-PS is higher than local-PS, we recall that our objective is not to surpass previous descriptors in term of accuracy. Rather, we relied on the histogram nature of pattern spectra to propose a novel descriptor that is computed quickly on a large image while achieving acceptable accuracy. We did not compare our proposed descriptor with state-of-the-art descriptors such as deep features and left it for future work. While other descriptors possibly lead to a higher accuracy, they cannot adapt to our framework for parallel computation since they have to be computed locally in each patch. We assume that the retrieval time is more important in a large-scale dataset where our proposed method is more efficient.

There is a lot of rooms for improvement of the proposed method. We used a max-tree structure and computed one dimensional pattern spectra based on the compactness attribute. However, our method is not restricted to the tree types and their attributes. We designed an efficient framework that can be realised with various types of pattern spectra to find the best combination of tree structures and attributes. Also, we simply calculated the distance between descriptors using χ 2 . So, we did not rely on a learnt distance metric. We expect that learning a weighting matrix to compute the distances will lead to a higher accuracy.

Fig. 4 :

 4 Fig. 4: (Left) Pattern spectra is a histogram-like distribution of image components. (Right) The components of each image patch contribute to a few numbers of histogram bins creating a sub-pattern spectra.

Fig. 7 :

 7 Fig. 7: Computation of fast pattern spectra consists of populating the local histograms for sliding windows. The algorithm traverses the tree structure and adds the contribution of each node to the corresponding histogram bin.

Fig. 8 :

 8 Fig. 8: Selected satellite images for experiments. (a) 617 ˆ647 RGB image that user selects query patches from it. Three selected query patches are marked with red rectangles. (b) Ground truth labels for query image pixels showing 6 different classes. (c) 633 ˆ668 RGB image to search for similar patches (d) ground truth labels for search image pixels.

Fig. 9 :

 9 Fig. 9: Precision-recall curve and average precision (AP) for three selected query patches with different sizes using baseline local-PS and proposed fast-PS.

Table 1 :

 1 Average precision and computation time for retrieval of three selected query patches using local pattern spectra and fast pattern spectra.

			Average Precision (%) Computation Time (seconds)
	query size (pixels) total patches local-PS	fast-PS local-PS fast-PS tree table
	(a) 80 ˆ80	210	30.83	38.89	3.91	1.97
	(b) 70 ˆ90	221	41.03	46.92	4.12	1.98	2.86
	(c) 90 ˆ70	234	30.54	41.39	4.38	2.01

Acknowledgment

This work was funded by DAJ-AR-NO-2018.0010814 project from CNES.