Bruno Bouchard
email: bouchard@ceremade.dauphine.fr

Xiaolu Tan
email: tan@ceremade.dauphine.fr

Xavier Warin
email: warin@edf.fr

Numerical approximation of general Lipschitz BSDEs with branching processes

Keywords: BSDE, Monte-Carlo methods, branching process Primary 65C05, 60J60; Secondary 60J85, 60H35

, this requires a localization procedure that uses a priori estimates on the true solution, so as to ensure the well-posedness of the involved Picard iteration scheme, and the global convergence of the algorithm. When, the nonlinearity depends on the gradient, the later needs to be controlled as well. This is done by using a face-lifting procedure. Convergence of our algorithm is proved without any limitation on the time horizon. We also provide numerical simulations to illustrate the performance of the algorithm.

Introduction

The aim of this paper is to extend the branching process based numerical algorithm proposed in Bouchard et al. [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF] to general BSDEs in form:

Y • = g(X T ) + T • f (X s , Y s , Z s ) ds - T • Z s dW s , (1.1)
where W is a standard d-dimensional Brownian motion, f :

R d × R × R d → R
is the driver function, g : R d → R is the terminal condition, and X is the solution of

X = X 0 + • 0 µ(X s ) ds + • 0 σ(X s ) dW s , (1.2)
with constant initial condition X 0 ∈ R d and coefficients (µ, σ) : R d → R d × R d×d , that are assumed to be Lipschitz1 . From the PDE point of view, this amounts to solving the parabolic equation

∂ t u + µ • Du + 1 2 Tr[σσ D 2 u] + f •, u, Duσ = 0, u(T, •) = g.
The main idea of [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF] was to approximate the driver function by local polynomials and use a Picard iteration argument so as to reduce the problem to solving BSDE's with (stochastic) global polynomial drivers, see Section 2, to which the branching process based pure forward Monte-Carlo algorithm of [START_REF] Henry-Labordère | Cutting CVA's Complexity[END_REF][START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF][START_REF] Henry-Labordere | A numerical algorithm for a class of BSDEs via the branching process[END_REF]] can be applied. See for instance [START_REF] Mckean | Application of Brownian motion to the equation of Kolmogorov-Petrovskii-Piskunov[END_REF][START_REF]Skorokhod Branching diffusion processes[END_REF][START_REF] Watanabe | On the branching process for Brownian particles with an absorbing boundary[END_REF] for the related Feynman-Kac representation of the KPP (Kolmogorov-Petrovskii-Piskunov) equation. This algorithm seems to be very adapted to situations where the original driver can be well approximated by polynomials with rather small coefficients on quite large domains. The reason is that, in such a situation, it is basically a pure forward Monte-Carlo method, see in particular [3, Remark 2.10(ii)], which can be expected to be less costly than the classical schemes, see e.g. [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF][START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF][START_REF] Bouchard | Monte-Carlo valuation of American options: facts and new algorithms to improve existing methods[END_REF][START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF][START_REF] Zhang | A numerical scheme for backward stochastic differential equations[END_REF] and the references therein. However, the numerical scheme of [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF] only works when the driver function (x, y, z) → f (x, y, z) is independent of z, i.e. the nonlinearity in the above equation does not depend on the gradient of the solution.

Importantly, the algorithm proposed in [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF] requires the truncation of the approximation of the Y -component at some given time steps. The reason is that BSDEs with polynomial drivers may only be defined up to an explosion time. This truncation is based on a priori estimates of the true solution. It ensures the well-posedness of the algorithm on an arbitrary time horizon, its stability, and global convergence.

In the case where the driver also depends on the Z component of the BSDE, a similar truncation has to be performed on the gradient itself. It can however not be done by simply projecting Z on a suitable compact set at certain time steps, since Z only maters up to an equivalent class of ([0, T ] × Ω, dt × dP).

Alternatively, we propose to use a face-lift procedure at certain time steps, see (2.10). Again this time steps depend on the explosion times of the corresponding BSDEs with polynomial drivers. Note that a similar face-lift procedure is used in Chassagneux, Elie and Kharroubi [START_REF] Chassagneux | A numerical probabilistic scheme for super-replication with Delta constraints[END_REF] 2 in the context of the discrete time approximation of BSDEs with contraint on the Z-component.

We prove the convergence of the scheme. The very good performance of this approach is illustrated in Section 4 by a numerical test case.

Notations: All over this paper, we view elements of R d , d ≥ 1, as column vectors. Transposition is denoted by the superscript . We consider a complete filtered probability space (Ω, F, F = (F t ) t≤T , P) supporting a ddimensional Brownian motion W . We simply write

E t [•] for E[•|F t ], t ≤ T .
We use the standard notations S 2 (resp. L 2 ) for the class of progressively measurable processes ξ such that ξ

S 2 := E[sup [0,T ] ξ 2 ] 1 2 (resp. ξ L 2 := E[ T 0 ξ s 2 ds] 1
2 ) is finite. The dimension of the process ξ is given by the context. For a map (t, x) → ψ(t, x), we denote by ∂ t ψ is derivative with respect to its first variable and by Dψ and D 2 ψ its Jacobian and Hessian matrix with respect to its second component.

Approximation of BSDE using local polynomial drivers and the Picard iteration

For the rest of this paper, let us consider the (decoupled) forward-backward system (1.1)-(1.2) in which f and g are both bounded and Lipschitz contin-uous, and σ is non-degenerate such that there is a constant a 0 > 0 satisfying

σσ (x) ≥ |a 0 | 2 I d , ∀x ∈ R d . (2.1)
We also assume that µ, σ, Dµ and Dσ are all bounded and continuous. In particular, (1.1)

-(1.2) has a unique solution (X, Y, Z) ∈ S 2 × S 2 × L 2 . The above conditions indeed imply that |Y |+ σ(X) -1 Z ≤ M on [0, T ],
for some M > 0.

Remark 2.1. The above assumptions can be relaxed by using standard localization or mollification arguments. For instance, one could simply assume that g has polynomial growth and is locally Lipschitz. In this case, it can be truncated outside a compacted set so as to reduce to the above. Then, standard estimates and stability results for SDEs and BSDEs can be used to estimate the additional error in a very standard way. See e.g. [START_REF] Karoui | Backward stochastic differential equations in finance[END_REF].

Local polynomial approximation of the generator

As in [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF], our first step is to approximate the driver f by a driver f • that has a local polynomial structure. The difference is that it now depends on both components of the solution of the BSDE. Namely, let

f • (x, y, z, y , z ) := j• j=1 ∈L c j, (x)y 0 q• q=1 (b q (x) z) q ϕ j (y , z ), (2.2)
in which (x, y, z, y , z

) ∈ R × R × R d × R × R d , L := {0, • • • , L • } q•+1 (with L • , q • ∈ N)
, the functions (b q ) 0≤q≤q• and (c j, , ϕ j ) ∈L,1≤j≤j• (with j • ∈ N) are continuous and satisfy

|c j, | ≤ C L , b q ≤ 1, |ϕ j | ≤ 1, and Dϕ j ≤ L ϕ , (2.3)
for all 1 ≤ j ≤ j • , 0 ≤ q ≤ q • and ∈ L, for some constants C L , L ϕ ≥ 0.

For a good choice of the local polynomial f • , we can assume that

(x, y, z) → f• x, y, z := f • x, y, z, y, z
is globally bounded and Lipschitz. Then, the BSDE

Ȳ• = g(X T ) + T • f• (X s , Ȳs , Zs ) ds- T • Z s dW s , (2.4)
has a unique solution ( Ȳ , Z) ∈ S 2 × L 2 , and standard estimates imply that ( Ȳ , Z) provides a good approximation of (Y, Z) whenever f• is a good approximation of f :

Y -Ȳ S 2 + Z -Z L 2 ≤ C • (f -f• )(X, Y, Z) L 2 , (2.5)
for some C • > 0 that depends on the global Lipschitz constant of f• (but not on the precise expression of f• ), see e.g. [START_REF] Karoui | Backward stochastic differential equations in finance[END_REF].

One can think at the (c j, ) ∈L as the coefficients of a polynomial approximation of f in terms of (y, (b q (x) z)) q≤q• on a subset A j , the A j 's forming a partition of [-M, M ] 1+d . Then, the ϕ j 's have to be considered as smoothing kernels that allow one to pass in a Lipschitz way from one part of the partition to another one. The choice of the basis functions (b q ) q≤q• as well as (ϕ j ) 1≤j≤j• will obviously depend on the application, but it should in practice typically be constructed such that the sets

A j := {(y, z) ∈ R × R d : ϕ j (y, z) = 1} (2.6)
are large and the intersection between the supports of the ϕ j 's are small. See [3, Remark 2.10(ii)] and below. Finally, since the function f• is chosen to be globally bounded and Lipschitz, by possibly adjusting the constant M , we can assume without loss of generality that

| Ȳ | + Z σ(X) -1 ≤ M. (2.7)
For later use, let us recall that Ȳ is related to the unique bounded and continuous viscosity solution ū of

∂ t ū + µ • Dū + 1 2 Tr[σσ D 2 ū] + f• •, ū, Dūσ = 0, ū(T, •) = g, though Ȳ = ū(•, X). (2.8) Moreover,
ū is bounded by M and M -Lipschitz.

Picard iteration with truncation and face-lifting

Our next step is to introduce a Picard iteration scheme to approximate the solution Ȳ of (2.4) so as to be able to apply the branching process based forward Monte-Carlo approach of [START_REF] Henry-Labordère | Cutting CVA's Complexity[END_REF][START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF][START_REF] Henry-Labordere | A numerical algorithm for a class of BSDEs via the branching process[END_REF] to each iteration: given ( Ȳ m-1 , Zm-1 ), use the representation of the BSDE with driver f • (X, •, •, Ȳ m-1 , Zm-1 ). However, although the map (y, z) → f• (x, y, z) = f • (x, y, z, y, z) is globally Lipschitz, the map (y, z) → f • (x, y, z, y , z ) is a polynomial, given fixed (x, y , z ), and hence only locally Lipschitz in general. In order to reduce to a Lipschitz driver, we need to truncate the solution at certain time steps, that are smaller than the natural explosion time of the corresponding BSDE with (stochastic) polynomial driver. As in [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF], it can be performed by a simple truncation at the level of the first component of the solution. As for the second component, that should be interpreted as a gradient, a simple truncation does not make sense, the gradient needs to be modified by modifying the function itself. Moreover, from the BSDE viewpoint, Z is only defined up to an equivalent class on ([0, T ] × Ω, dt × dP), so that changing its value at a finite number of given times does not change the corresponding BSDE. We instead use a face-lifting procedure, as in [START_REF] Chassagneux | A numerical probabilistic scheme for super-replication with Delta constraints[END_REF].

More precisely, let us define the operator R on the space of bounded functions ψ : R d → R by

R[ψ] := (-M ) ∨ sup p∈R d ψ(• + p) -δ(p) ∧ M,
where

δ(p) := sup q∈[-M,M ] d (p q) = M d i=1 |p i |.
The operation ψ → sup p∈R d ψ(• + p) -δ(p) maps ψ into the smallest M -Lipschitz function above ψ. This is the so-called face-lifting procedure, which has to be understood as a (form of) projection on the family of M -Lipschitz functions, see e. Let now (h • , M h• ) be such that (A.3) and (A.4) in the Appendix hold. The constant h • is a lower bound for the explosion time of the BSDE with driver (y, z) → f • (x, y, z, y , z ) for any fixed (x, y , z ). Let us then fix h ∈ (0, h • ) such that N h := T /h ∈ N, and define

t i = ih and T := {t i , i = 0, • • • , N h }.
Our algorithm consists in using a Picard iteration scheme to solve (2.4), which re-localize the solution at each time step of the grid T by applying operator R.

Namely, using the notation X t,x to denote the solution of (1.2) on [t, T ] such that X t,x t = x ∈ R d , we initialize our Picard scheme by setting

Y x,0 T = Ȳ x,0 T = g(x) (Y x,0 , Z x,0 ) = ( Ȳ x,0 , Zx,0 ) = (y, Dy)(•, X t i ,x ) on [t i , t i+1 ), i ≤ N h -1, in which y is a continuous function, M -Lipschitz in space, continuously dif- ferentiable in space on [0, T ) × R d and such that |y| ≤ M and y(T, •) = g. Then, given ( Ȳ x,m-1 , Zx,m-1 ), for m ≥ 1, we define ( Ȳ x,m , Zx,m ) as follows: 1. For i = N h , set ūm t i = ūm T := g 2. For i < N h , given ( Ȳ x,m-1 , Zx,m-1 ): (a) Let (Y x,m • , Z x,m • ) be the unique solution on [t i , t i+1 ) of Y x,m • = ūm t i+1 (X t i ,x t i+1 ) + t i+1 • f • X t i ,x s , Y x,m s , Z x,m s , Ȳ x,m-1 s , Zx,m-1 s ds - t i+1 • (Z x,m s ) dW s . (2.9) (b) Let u m t i : x ∈ R d → Y x,m t i
, and set

ūm t i := R[u m t i ]. (2.10) (c) Set Ȳ x,m := Y x,m on (t i , t i+1 ], Ȳ x,m t i := ūm t i (x), and Zx,m := Z x,m on [t i , t i+1 ), for x ∈ R d .
3. We finally define Ȳ m T = g(X T ) and

( Ȳ m , Zm ) := ( Ȳ Xt i ,m , ZXt i ,m ) on [t i , t i+1 ), i ≤ N h . (2.11)
In above, the existence and uniqueness of the solution (Y x,m , Z x,m ) to (2.9) is ensured by Proposition A.1. The projection operation in (2.10) is consistent with the behavior of the solution of (2.4), recall (2.7), and it is crucial to control the explosion of ( Ȳ m , Zm ) and therefore to ensure both the stability and the convergence of the scheme. This procedure is non-expansive, as explained in the following Remark, and therefore can not alter the convergence of the scheme.

Remark 2.2. Let ψ, ψ be two measurable and bounded maps on R d . Then,

sup p∈R d |ψ(• + p) -δ(p) -ψ (• + p) + δ(p)| = sup x∈R d |ψ(x) -ψ (x)|, and therefore R[ψ] -R[ψ ] ∞ ≤ ψ -ψ ∞ .
In particular, since ū defined through (2.8) is M -Lipschitz in its space variable and bounded by M , we have ū(t,

•) = R[ū(t, •)] for t ≤ T and therefore R[ψ] -ū(t, •) ∞ ≤ ψ -ū(t, •) ∞
for all t ≤ T and all measurable and bounded map ψ.

Also note that, if we had

( Ȳ m-1 t , Zm-1 t ) ∈ A j if and only if ( Ȳt , Zt ) ∈ A j , for all j ≤ j • , then we would have ( Ȳ m-1 , Zm-1 ) = ( Ȳ , Z), recall (2.6
) and the definition of f• in terms of f • . This means that we do not need to be very precise on the original prior, whenever the sets A j can be chosen to be large.

From the theoretical viewpoint, the error due to the above Picard iteration scheme can be deduced from classical arguments. Recall that (h • , M h• ) is such that (A.3) and (A.4) in the Appendix hold.

Theorem 2.1. For each m ≥ 0, the algorithm defined in 1.-2.-3. above provides the unique solution

( Ȳ m , Zm ) ∈ S 2 × L 2 . Moreover, it satisfies | Ȳ m | ∨ Zm σ(X) -1 ≤ M h• ,
and there exists a measurable map (ū m , vm ) :

[0, T ] × R d → R 1+d , that is continuous on ∪ i<N h (t i , t i+1 ) × R d , such that ūm (t i , •) is continuous on R d for all i ≤ N h , and
Ȳ m = ūm (•, X) on [0, T ] P-a.s.
(2.12)

Zm = vm (•, X) dt × dP-a.e. on [0, T ] × Ω. vm = Dū m σ on ∪ i<N h (t i , t i+1 ) × R d .
Moreover, for any constant ρ ∈ (0, 1), there is some constant

C ρ > 0 such that | Ȳ m t -Ȳt | 2 + E t [ T t Zm s -Zs 2 ds] ≤ C ρ ρ m , for all t ≤ T .
Proof. i) Recall from Remark 2.2 that R maps bounded functions into M -Lipschitz functions that are bounded by M . Then, by Proposition A.1 in the Appendix, the solutions (Y x,m , Z x,m ) as well as ( Ȳ x,m , Zx,m ) are uniquely defined in and below (2.9). Moreover, one has

| Ȳ x,m t |∨ ( Zx,m t ) σ(X t i ,x t ) -1 ≤ M h• , for all x ∈ R d , i < N h and t ∈ [t i , t i+1 ). As a consequence, ( Ȳ m , Zm ) is also uniquely defined and satisfies | Ȳ m | ∨ Zm σ(X) -1 ≤ M h• .
Using again Proposition A.1, one has the existence of (ū m , vm ) satisfying the condition in the statement. ii) We next prove the convergence of the sequence ( Ȳ m , Zm ) m≥0 to ( Ȳ , Z). Since {( Ȳ x,m , Zx,m ), x ∈ R d } is uniformly bounded, the generator f • (x, y, z, y , z ) in (2.9) can be considered to be uniformly Lipschitz in (y, z) and (y , z ). Assume that the corresponding Lipschitz constants are L 1 and L 2 . Let us set Θ x,m := (Y x,m , Z x,m ) and define (∆Y x,m , ∆Z x,m ) := (Y x,m -Ȳ x , Z x,m -Zx ) where Θx := ( Ȳ x , Zx ) denotes the solution of

Ȳ x • = ū(t i+1 , X t i ,x t i+1 ) + t i+1 • f• X t i ,x s , Ȳ x s , Zx s ds - t i+1 • ( Zx s ) dW s ,
on each [t i , t i+1 ], recall (2.8). In the following, we fix β > 0 and use the notation

ξ β,t := E t [ t i+1 t e βs |ξ s | 2 ds] 1 2 for ξ ∈ L 2 , t ∈ [t i , t i+1 ). Fix t ∈ [t i , t i+1 )
. By applying Itô's formula to (e βs (∆Y x,m+1 s ) 2 ) s∈[t,t i+1 ] and then taking expectation, we obtain

e βt |∆Y x,m+1 t | 2 + β ∆Y x,m+1 2 β,t + 2 ∆Z x,m+1 2
β,t ≤ E t e βt i+1 (∆Y x,m+1 t i+1 -) 2 + 2E t t i+1 t e βs ∆Y x,m+1 s f • (X t i ,x s , Θ x,m+1 s , Θ x,m s ) -f • (X t i ,x s , Θx s , Θx s ) ds .
Using the Lipschitz property of f • and the inequality λ + 1 λ ≥ 2 for all λ > 0, it follows that, for all λ 1 , λ 2 > 0,

e βt |∆Y x,m+1 t | 2 + (β -(2L 1 + λ 1 L 1 + λ 2 L 2 )) ∆Y x,m+1 2 β,t (2.13) + (2 - L 1 λ 1 ) ∆Z x,m+1 2
β,t ≤ E t e βt i+1 |∆Y x,m+1 t i+1 -| 2 + L 2 λ 2 ∆Y x,m 2 β,t + L 2 λ 2 ∆Z x,m 2 β,t . iii) Let us now choose 1 > ρ = ρ 0 > ρ 1 > • • • > ρ N h > 0 such that (m + 1)e βh ≤ ρ k ρ k+1 m+1
for all m ≥ 0.

(2.14)

For i = N h -1, we have ∆Y x,m t i+1 -= 0 for all m ≥ 1. Choosing λ 1 , λ 2 and β > 0 in (2.13) such that L 2 λ 2 1 β -(2L 1 + λ 1 L 1 + λ 2 L 2 ) ≤ ρ N h and L 2 λ 2 1 2 -L 1 /λ 1 ≤ ρ N h , it follows from (2.13) that, for t ∈ [t N h -1 , t N h ), ∆Y x,m+1 2 β,t + ∆Z x,m+1 2 β,t ≤ C(ρ i ) m+1
, for m ≥ 0, where

C := esssup sup (s,x )∈[0,T ]×R d e βT |∆Y x ,0 s | 2 + ∆Z x ,0 s 2 < ∞,
and then, by (2.13) again,

|∆Y x,m+1 t | 2 ≤ C(ρ i ) m+1 , for t ∈ [t i , t i+1 ), i = N h -1, m ≥ 0.
Recalling Remark 2.2, this shows that

| Ȳ x,m+1 t -Ȳ x t | 2 ≤ C i (ρ i ) m+1 , for t ∈ [t i , t i+1 ), i = N h -1, m ≥ 0, (2.15)
in which

C N h -1 := C.
Assume now that (2.15) holds true for i + 1 ≤ N h and some given C i+1 > 0.

Recall that ρ i ≥ ρ N h . Applying (2.13) with the above choice of λ 1 , λ 2 and β, we obtain

∆Y x,m+1 2 β,t + ∆Z x,m+1 2 β,t ≤ e βh C i+1 (ρ i+1 ) m+1 + ρ i ∆Y x,m 2 β,t + ∆Z x,m 2 β,t ,
which, by (2.14) and the fact that ρ i < 1, induces that

∆Y x,m+1 2 β,t + ∆Z x,m+1 2 β,t ≤ (m + 1)e βh C i+1 (ρ i+1 ) m+1 + (ρ i ) m+1 ∆Y x,0 2 β,t + ∆Z x,0 2 β,t ≤ C i (ρ i ) m+1 (2.16)
where

C i := C i+1 + C.
Let us further choose λ 2 > 0 such that L 2 /λ 2 ≤ ρ N h , and recall that ρ i ≥ ρ N h . Then, using again (2.13), (2.14), (2.15) applied to i + 1, we obtain, for t ∈ [t i , t i+1 ),

|∆Y x,m+1 t | 2 ≤ e βh C i+1 (ρ i+1 ) m+1 + C i (ρ i ) m+1 ≤ C i (ρ i ) m+1 ,
so that it follows from Remark 2.2 that

| Ȳ x,m+1 t -Ȳ x t | 2 ≤ C i (ρ i ) m+1 , (2.17)
where

C i := e βh C i+1 + C i . Since ( Ȳ , Z) = ( Ȳ Xt i , ZXt i ) and ( Ȳ m , Zm ) = ( Ȳ Xt i ,m , ZXt i ,m ) on each [t i , t i+1 ),
this concludes the proof.

A branching diffusion representation for Ȳ m

We now explain how the solution of (2.9) on [t i , t i+1 ) can be represented by means of a branching diffusion system. We slightly adapt the arguments of [START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF].

Let us consider an element (p

) ∈L ⊂ R + such that ∈L p = 1, set K n := {(1, k 2 , . . . , k n ) : (k 2 , . . . , k n ) ∈ {1, . . . , (q • + 1)L • } n-1 } for n ≥ 1, and K := ∪ n≥1 K n . Let (W k
) k∈K be a sequence of independent d-dimensional Brownian motions, (ξ k = (ξ k,q ) 0≤q≤q• ) k∈K and (δ k ) k∈K be two sequences of independent random variables, such that

P[ξ k = ] = p , ∈ L, k ∈ K, and
F (t) := P[δ k > t] = ∞ t ρ(s)ds, t ≥ 0, k ∈ K,
for some continuous strictly positive map ρ : R + → R + . We assume that (W k ) k∈K , (ξ k ) k∈K , (δ k ) k∈K and W are independent.

Given the above, we construct particles X (k) that have the dynamics (1.2) up to a killing time T k at which they split in ξ k 1 := ξ k,0 + • • • ξ k,q• different (conditionally) independent particles with dynamics (1.2) up to their own killing time. The construction is done as follows. First, we set T (1) := δ 1 , and, given k

= (1, k 2 , . . . , k n ) ∈ K n with n ≥ 2, we let T k := δ k + T k-in which k-:= (1, k 2 , . . . , k n-1 ) ∈ K n-1 .
We can then define the Brownian particles (W (k) ) k∈K by using the following induction: we first set

W ((1)) := W 1 1 [0,T (1) ] , K 1 t := {(1)}1 [0,T (1) ] (t) + ∅1 [0,T (1) ] c (t), t ≥ 0,
then, given n ≥ 2 and k ∈ Kn-1

T := ∪ t≤T K n-1 t
, we let

W (k⊕j) := W (k) •∧T k + W k⊕j •∨T k -W k⊕j T k 1 [0,T k⊕j ] , 1 ≤ j ≤ ξ k 1 ,
in which we use the notation

(1, k 1 , . . . , k n-1 ) ⊕ j := (1, k 1 , . . . , k n-1 , j), and
Kn t := {k ⊕ j : k ∈ Kn-1 T , 1 ≤ j ≤ ξ k 1 s.t. t ∈ (0, T k⊕j ]}, Kt := ∪ n≥1 Kn t , K n t := {k ⊕ j : k ∈ Kn-1 T , 1 ≤ j ≤ ξ k 1 s.t. t ∈ (T k , T k⊕j ]}, K t := ∪ n≥1 K n t .
Now observe that the solution X x of (1.2) on [0, T ] with initial condition X x 0 = x ∈ R d can be identified in law on the canonical space as a process of the form Φ[x](•, W ) in which the deterministic map (x, s, ω) → Φ[x](s, ω) is B(R d ) ⊗ P-measurable, where P is the predictable σ-filed on [0, T ] × Ω. We then define the corresponding particles (X x,(k) ) k∈K by X x,(k) := Φ[x](•, W (k) ). Moreover, we define the d × d-dimensional matrix valued tangent process

∇X x,(k) defined on [[T k-, T k ]] by ∇X x,(k) T k-= I d (3.1) d∇X x,(k) t = Dµ(X x,(k) t )∇X x,(k) t dt + d i=1 Dσ i (X x,(k) t )∇X x,(k) t dW (k),i t ,
where I d denotes the d × d-dimensional identity matrix, and σ i denotes the i-th column of matrix σ.

Finally, we give a mark 0 to the initial particle (1), and, for every particle k, knowing ξ k = (ξ k,0 , • • • , ξ k,q• ), we consider its offspring particles (k ⊕ j) j=1,••• , ξ k 1 and give the first ξ k,0 particles the mark 0, the next ξ k,1 particles the mark 1, the next ξ k,2 particles the mark 2, etc. Thus, every particles k carries a mark θ k taking values in 0 to q • .

Given the above construction, we can now provide the branching process based representation of ( Ȳ m ) m≥0 . We assume here that (ū m-1 , vm-1 ) defined in (2.12) are given for some m ≥ 1, recall that (ū 0 , v0 ) = (y, Dy) by construction. We set ũm (T,

•) = g. Then, for i = N h -1, • • • , 0, we define (ũ m , ṽm ) on each interval [t i , t i+1 ) recursively by ũm (t, x) := E U m t,x 1 t =t i + R[E U m t,• ](x)1 t=t i ṽm (t, x) := E V m t,x σ(x), (3.2) for (t, x) ∈ [t i , t i+1 ) × R d , in which U m t,x := k∈K t i+1 -t G m t,x (k)W t,x (k) k∈ Kt i+1 -t \K t i+1 -t A m t,x (k)W t,x (k) V m t,x := U m t,x 1
T (1)

T (1) 0 σ -1 (X x, (1) s
)∇X x, (1) s dW (1) r where

G m t,x (k) := ũm t i+1 , X x,(k) t i+1 -t -ũm t i+1 , X x,(k-) T k- 1 {θ(k) =0} F (t i+1 -t -T k-) , A m t,x (k) := j• j=1 c j,ξ k X x,(k) T k ϕ j (ū m-1 , vm-1 )(t + T k , X x,(k) T k ) p ξ k ρ(δ k ) , W t,x (k) := 1 {θ k =0} + 1 {θ k =0} T k -T k- b θ k (X x,(k) T k ) • T k T k- σ -1 (X x,(k) s )∇X x,(k) s dW (k) r .
Compare with [13, (3.4) and (3.10)].

The next proposition shows that (ũ m , ṽm ) actually coincides with (ū m , vm ) in (2.12), a result that follows essentially from [START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF]. Nevertheless, to be more precise on the square integrability of U m t,x and V m t,x , one will fix a special density function ρ as well as probability weights (p ) ∈L . Recall again that (ū m , vm ) are defined in Theorem 2.1 and satisfy (2.12), and that (h • , M • ) are chosen such that (A.3)-(A.4) in the Appendix hold.

Proposition 3.1. Let us choose ρ(t) = 1 3 t -2/3 1 {t∈[0,1]} and p = c ∞ c 1,∞ with c 1,∞ := ∈L c ∞ .
Let h • and M h • be given by (A.8) and (A.9). Assume that h ∈ (0,

h • ∧ h • ). Then, E[|U m t,x | 2 ] ∨ E[ V m t,x 2 ] ≤ (M h • ) 2 , for all m ≥ 1 and (t, x) ∈ [0, T ] × R d . Moreover, ũm = ūm on [0, T ] × R d and ṽm = vm on ∪ i≤N h -1 (t i , t i+1 ) × R d .
The proof of the above mimics the arguments of [START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF]Theorem 3.12] and is postponed to the Appendix A.2.

Remark 3.1. The integrability and representation results in Proposition 3.1 hold true for a large class of parameters ρ, (p ) ∈L and (c ) ∈L (see e.g. [13, Section 3.2] for more details). We restrict to a special choice of parameters in Proposition 3.1 in order to compute explicitly the lower bound h • for the explosion time as well as the upper bound M h • for the variance of the estimators.

Remark 3.2. a. The above scheme requires the computation of conditional expectations that involve the all path (ũ m , ṽm )(•, X t i ,x ) on each (t i , t i+1 ), for all x ∈ R d , and therefore the use of an additional time space grid on which ũm and ṽm are estimated by Monte-Carlo. For t i < t < t i+1 , the precision does not need to be important because the corresponding values are only used for the localization procedure, see the discussion just after Remark 2.2, and the grid does not need to be fine. The space grid should be finer at the times in T because each ũm (t i , •) is used per se as a terminal condition on (t i , t i+1 ],

and not only for the localization of the polynomial. This corresponds to the Method B in [3, Section 3]. One can also consider a very fine time grid T and avoid the use of a sub-grid. This is Method A in [3, Section 3]. The numerical tests performed in [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF] suggest that Method A is more efficient. b. From a numerical viewpoint, ṽm can also be estimated by using a finite difference scheme based on the estimation of ũm . It seems to be indeed more stable. c. Obviously, one can also adapt the algorithm to make profit of the ghost particle or of the normalization techniques described in [START_REF] Warin | Variations on branching methods for non linear PDEs[END_REF], which seem to reduce efficiently the variance of the estimations, even when ρ is the exponential distribution.

Numerical example

This section is dedicated to a simple example in dimension one showing the efficiency of the proposed methodology. We use the Method A of [3, Section 3] together with the normalization technique described in [START_REF] Warin | Variations on branching methods for non linear PDEs[END_REF] and a finite difference scheme to compute ṽm . In particular, this normalization technique allows us to take ρ as an exponential density rather than that in Proposition 3.1. See Remark 3.2.

We consider the SDE with coefficients

µ(x) = -0.5(x + 0.2) and σ(x) = 0.1 (1 + 0.5((x ∨ 0) ∧ 1)) .
The maturity is T = 1 and the non linearity in (1.1) is taken as

f (t, x, y, z) = f (y, z) + 1 2 e t-T 2 cos(x) σ(x) 2 2 - 1 2 (1 + cos(x)) + µ(x) sin(x) - 2 4 + | sin(x)(1 + cos(x))|e t-T where f (y, z) = 1 2(1 + |yz|) . (4.1)
It is complemented by the choice of the terminal condition g(x) = 1+cos(x) 2 , so that an analytic solution is available :

u(t, x) = 1 + cos(x) 2 e t-T 2 .
We use the algorithm to compute an estimation of u(0, •) on X := [-1., 1.].

To construct our local polynomials approximation of f , we use a linear spline interpolation in each direction, obtained by tensorization, and leading to a local approximation on the basis 1, y, z, yz on each mesh of a domain The parameters affecting the convergence of the algorithm are:

[0, 1] × [-1, 1].
• The couple of meshes (n y , n z ) used in the spline representation of (4.1), where n y (resp. n z ) is the number of linear spline meshes for the y (resp. z) discretization.

• The number of time steps N h .

• The grid and the interpolation used on X at t 0 = 0 and for all dates t i , i = 1, . . . , N h . Note that the size of the grid has to be adapted to the value of T , because of the diffusive feature of (1.2). All interpolations are achieved with the StOpt library (see [START_REF] Gevret | The StOpt library[END_REF][START_REF] Gevret | STochastic OPTimization library in C++[END_REF]) using a modified quadratic interpolator as in [START_REF] Warin | Some Non-monotone Schemes for Time Dependent Hamilton-Jacobi-Bellman Equations in Stochastic Control[END_REF]. In the following, ∆x denotes the mesh of the space discretization.

• The time step is set to 0.002 and we use an Euler scheme to approximate (1.2).

• The accuracy of the estimation of the expectations appearing in our algorithm. We compute the empirical standard deviation θ associated to each Monte Carlo estimation of the expectation in (3.2). We try to fix the number M of samples such that θ/ M does not exceed a certain level, fixed at 0.000125, at each point of our grid. We cap this number of simulations at 510 5 .

• The intensity, set to 0.4, of the exponential distribution used to define the random variables (δ k ) k∈K .

Finally, we take M = 1 in the definition of R.

We only perform one Picard iteration with initial prior (ũ 0 , ṽ0 ) = (g, Dgσ).

On the different figures below, we plot the errors obtained on X for different values of N h , (n y , n z ) and ∆x. We first use 20 time steps and an interpolation step of 0.1 In figure 2, we display the error as a function of the number of spline meshes. We provide two plots:

• On the left-hand side, n y varies above 5 and n z varies above 10,

• On the right-hand side, we add (n y , n z ) = [START_REF] Bouchard | Monte-Carlo valuation of American options: facts and new algorithms to improve existing methods[END_REF][START_REF] Bouchard | Monte-Carlo valuation of American options: facts and new algorithms to improve existing methods[END_REF]. It leads to a maximal error of 0.11, showing that a quite good accuracy in the spline representation in z is necessary. In Figure 4, we finally let the number of time steps N h vary. Once again we give two plots:

• one with N h above or equal to 20,

• one with small values of N h .

The results clearly show that the algorithm produces bad results when N h is too small: the time steps are too large for the branching method. In this case, it exhibits a large variance. When N h is too large, then interpolation errors propagate leading also to a deterioration of our estimations. Numerically, it can be checked that the face-lifting procedure is in most of the cases useless when only one Picard iteration is used:

• When the variance of the branching scheme is small, the face-lifting and truncation procedure has no effect,

• When the variance becomes too large, the face-lifting procedure is regularizing the solution and this permits to reduce the error due to our interpolations.

In figure 5, we provide the estimation with and without face-lifting, obtained with N h = 10, (n y , n z ) = (20, 10) and a space discretization ∆x = 0.1. For N h = 20 the computational time is equal to 320 seconds on a 3 year old laptop.

A Appendix

A.1 A priori estimates for the Picard iteration scheme

In this section, we let ∇X be the tangent process associated to X on [0, T ] by

∇X 0 = I d , d∇X t = Dµ(X t )∇X t dt + d i=1 Dσ i (X t )∇X t dW i t ,
and we define

N t s := s t 1 s -t (σ(X r ) -1 ∇X r ) dW r ∇X -1 t
Proof. We construct the required solution by using Picard iterations. We set ( Ỹ 0 , Z0 ) = (y, Dy)(•, X), and define recursively on [0, T ] the couple ( Ỹ n+1 , Zn+1 ) as the unique solution of

Ỹ n+1 • = g(X T ) + T • f • (X s , Ỹ n s , Zn s , ( Ũ , Ṽ )(s, X s )) ds- T • ( Zn+1 s ) dW s ,
whenever it is well-defined. It is the case for n = 0. We now assume that ( Ỹ n , Zn ) is well-defined and such that

| Ỹ n | ∨ σ(X) -1 Zn ≤ M h• on [0, T ] for some n ≥ 0. Then, | Ỹ n+1 • | ≤ g ∞ + sup D M h• |f • | h • ≤ M h• ,
in which we used (A.3) for the second inequality. On the other hand, up to using a mollifying argument, one can assume that g is C 1 b and that (U, V ) is Lipschitz. Then, it follows from the same arguments as in [15, Theorem 3.1, Theorem 4.2] that Zn+1 admits the representation

( Zn+1 t ) =E t Dg(X T )∇X T ∇X -1 t σ(X t ) + E t T t f • (X s , Ỹ n s , Zn s , (U, V )(s, X s ))N t s ds σ(X t ).
By combining the above together with (A.1) and (A.2), we obtain that

( Zn+1 t ) σ(X t ) -1 ≤ M e Cµ,σ(T -t) + C µ,σ sup D M h• |f • | (h • ) 1 2 ≤ M h• ,
in which we used (A.4) for the second inequality. The above proves that the sequence ( Ỹ n , Zn ) n≥0 is uniformly bounded on [0, T ]. Therefore, we can consider f • as a Lipschitz generator, and hence ( Ỹ n , Zn ) n≥0 is in fact a Picard iteration that converges to a solution of (A.5) with the same bound. The existence of the maps ũ and ṽ such that ṽ = Dũσ follows from [15, Theorem 3.1] applied to (A.5) when (g, f [START_REF] Ma | Representation theorems for backward stochastic differential equations[END_REF]Theorem 4.2] combined with a simple approximation argument, see e.g. [8, (ii) of the proof of Proposition 3.2], then shows that the same holds on [0, T ) under our conditions.

• (•, y, z, (U, V )(t, •))) is C 1 b , uniformly in (t, y, z) ∈ [0, T ] × R d+1 . The representation result of

A.2 Proof of the representation formula

We adapt the proof of [13, Theorem 3.12] to our context. We proceed by induction. In all this section, we fix

(t, x) ∈ [t i , t i+1 ) × R d ,
and assume that the result of Proposition 3.1 holds up to rank m -1 ≥ 0 on [0, T ] (with the convention U 0 • = y, V 0 • := Dy), and up to rank m on [t i+1 , T ]. In particular, we assume that ũm (t i+1 , •) = ūm (t i+1 , •). We fix = 4 and define

C 1, := M ∨ sup t≤s, x∈R d , q=1,••• ,q•, ξ ≤M E ξ • (X t,x s -x) b q (x) • W t,x,s , C 2, := sup t≤s≤t i+1 , x∈R d , q=1,••• ,q• E √ s -t b q (x) • W t,x,s where W t,x,s := 1 s -t s t σ -1 (X t,x r )∇X t,x r dW r
in which ∇X t,x is the tangent process of X t,x with initial condition I d at t. We then set

Ĉ1, := C 1, F (T ) -1 , Ĉ2, := C 2, j • sup j≤j•, ∈L,t∈(0,h•] c j, ∞ p t -/(2( -1)) ρ(t) -1
.

Since F is non-increasing and ũm (t i+1 , •) = ūm (t i+1 , •) is bounded by M and M -Lipschitz, direct computations imply that Proof. Let b ∈ R d be a fixed vector. Set Q t,x := ∇X t,x b. Then, it follows from direct computations that

E[|U m t,x | ] ∨ E V m t,x ] ≤ E
E max [t,t i+1 ] Q t,x ≤ C 0 b .
Further, remember that each b q is assumed to be bounded by 1, so that b q σ -1 2 is uniformly bounded by λ(σσ ) -1 . Then, direct computations lead to We can now conclude the proof of Proposition 3.1.

C 1, ≤ 1 ∨ M ∨ 2 -1 C 0 + µ q T 2 C bdg C 0 λ(σσ ) -
Proof of Proposition 3.1. In view of (A.7), Lemma A.2 implies that {|U m t,x | 2 + V m t,x 2 , (t, x) ∈ [t i , t i+1 )×R d } is uniformly integrable (with a bound that does not depend on i < N h for 0 < h ≤ h • ). Then, arguing exactly as in [START_REF] Bouchard | Numerical approximation of BSDEs using local polynomial drivers and branching processes[END_REF]Proposition A.2] leads to ũm = ūm on [t i , t i+1 ). Combined with [START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF]Proposition 3.7], the uniform integrability also implies that Dũ m = ṽm σ on (t i , t i+1 ) × R d , and one can conclude from Theorem 2.1 that ṽm = ũm on (t i , t i+1 ) × R d . By the induction hypothesis of the beginning of this section, this proves that the statements of Proposition 3.1 hold.

Remark A.1. The constants Ĉ1 and Ĉ2 (and hence h • and M h • ) are clearly not optimal for applications. For instance, if σ ≡ σ • , for some non-degenerate constant matrix, the constants C 1, and C 2, can be significantly simplified as shown in [START_REF] Henry-Labordere | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF]Remark 3.9].

g. [ 2 ,

 2 Exercise 5.2], see also Remark 2.2 below. The outer operations in the definition of R are just a natural projection on [-M, M ].

Figure 1

 1 displays the function f and the error obtained by a discretization of 20 × 20 meshes. The driver f Error on the driver due to the linear spline representation.

Figure 1 :

 1 Figure 1: The driver f and its linear spline representation error for 20 × 20 splines.

Figure 2 :

 2 Figure 2: Error plot depending on the couple (n y , n z ) for 20 time steps, a space discretization ∆x = 0.1.

Figure 3 :

 3 Figure 3: Error plot depending on ∆x for (n y , n z ) = (20, 10), N h = 20.

Figure 4 :

 4 Figure 4: Error plot depending on N h for (n y , n z ) = (20, 10), a space discretization ∆x = 0.1

Figure 5 :

 5 Figure 5: Error plot for N h = 10 for (n y , n z ) = (20, 10), a space discretization ∆x = 0.1 with and without face-lifting.

  -T k-) k∈ Kt\Kt Ĉ2, p ξ k ρ(δ k ) . (A.7)We first estimate the right-hand side, see (A.10) below.Let us denote by C bdg the constant in the Burkholder-Davis-Gundy inequality such thatE sup 0≤t≤T |M t | ] ≤ C bdg E[( M T ) 2 ]for any continuous martingale M with M 0 = 0. Denote furtherC 0 := (3 × 3) -1 1 + ( /( -1)) 1 + (1 + | λDµ T | e C Q T ) ,where λDµ the largest eigenvalue of the matrix Dµ, λDσ the largest eigenvalue of the matrix Dσ i , i ≤ d, and C Q := λDµ + d ( -1) λDσ /2. Define also λ(σσ ) -1 as the largest eigenvalue of matrix (σσ ) -1 .Lemma A.1. Under the Assumptions of Proposition 3.1,Ĉ1, ≤ Ĉ1 := 2 1 ∨ M ∨ 2 -1 (M √ d) C 0 + µ T 2 C BDG C 0 λ(σσ ) -1 2 ,and Ĉ2, ≤ Ĉ2 := C bdg C 0 λ(σσ ) -1 2 .

1 2 , and C 2 , 2 . 1 ( 1 - 1 . (A. 9 )Lemma A. 2 .

 2211192 ≤ C bdg C 0 λ(σσ ) -1It remains to use our specific choice of ρ and (p ) ∈L in Proposition 3.1 to conclude.Let us now choose h • and Mh • such that h • < 1 ∧ Ĉ-(|L|-1) |L| + 1)(|L| -1) Ĉ2 , (A.8) and (M h • ) 4 := Ĉ1-|L| h • (|L| + 1)(|L| -1) Ĉ2(1-|L|) -Let the conditions of Proposition 3.1 hold. Then, the ordinary differential equation η (t) = ∈L Ĉ2 η(t) 1 with initial condition η(0) = Ĉ1 ≥ 1 has a unique solution on [0, h • ], and it is bounded by(M h • ) 2 . Moreover, E k∈Kt Ĉ1 F (t -T k-) k∈ Kt\Kt Ĉ2 p ξ k ρ(δ k ) ≤ η(t) ≤ (M h • ) 4 , (A.10) for all t ∈ [0, h • ].Proof. The result follows from exactly the same arguments as in [3, Lemma A.1].

As usual, we could add a time dependency in the coefficients f , µ and σ without any additional difficulty.

We are grateful to the authors for several discussions on this subject.

Acknowledgements

This work has benefited from the financial support of the Initiative de Recherche "Méthodes non-linéaires pour la gestion des risques financiers" sponsored by AXA Research Fund. Bruno Bouchard and Xavier Warin acknowledge the financial support of ANR project CAESARS (ANR-15-CE05-0024).. Xiaolu Tan acknowledges the financial support of the ERC 321111 Rofirm, the ANR Isotace (ANR-12-MONU-0013), and the Chairs Financial Risks (Risk Foundation, sponsored by Société Générale) and Finance and Sustainable Development (IEF sponsored by EDF and CA).

for t ≤ s ≤ T . Standard estimates lead to

for some C µ,σ > 0 that only depends on µ ∞ , σ ∞ , Dµ ∞ , Dσ ∞ , a 0 in (2.1) and T . In particular, it does not depend on X 0 . Up to changing this constant, we may assume that

and let M h• ≥ M and h • ∈ (0, T ] be such that

and

(A.4)

The existence of h • and M h• follows from (2.2) and (2.3). Note that they do not depend on X 0 .

Proposition A.1. Let g : R d → R be bounded by M and M -Lipschitz. Fix

Then, there exists a unique bounded solution on

Moreover, there exists a bounded continuous map (ũ, ṽ) : [0, T ]×R d → R×R d such that Ỹ = ũ(•, X) on [0, T ] P-a.s. and Z = ṽ(•, X) dt × dP-a.e. on [0, T ] × Ω. It satisfies ṽ = Dũσ on [0, T ).


Amelia Carolina Sparavigna

Thermal Conductivity of the Crystalline Silicon

Keywords: Solid State physics, Materials Science, Thermal Conductivity, Crystalline Silicon, Isotopic Mass Defects, Phonon Scattering

In this paper we discuss and compare theoretical and experimental data of thermal conductivity in crystalline Silicon (natural and enriched). We focus in particular on the role that the isotopic mass defects have in determining the value of this conductivity.

Introduction

In the theoretical study of the thermal conductivity in non-metallic solids, a widely used approach is that of solving the Boltzmann equation of the phonon thermal transport to obtain the related conductivity [START_REF] Ziman | Electrons and phonons: The theory of transport phenomena in solids[END_REF][START_REF] Gurevich | Transport in phonon systems[END_REF][START_REF] Srivastava | The physics of phonons[END_REF]. Many researchers are using methods of solution based on the relaxation time approximation, of which we proposed a simplified discussion in [START_REF] Sparavigna | The Boltzmann equation of phonon thermal transport solved in the relaxation time approximation-I-Theory[END_REF][START_REF] Sparavigna | The Boltzmann Equation Of Phonon Thermal Transport Solved In the Relaxation Time Approximation -II -Data Analysis[END_REF]. Other researchers prefer the use of methods which are involving more specific atomic models [START_REF] Omini | Beyond the isotropic-model approximation in the theory of thermal conductivity[END_REF][START_REF] Omini | Heat transport in dielectric solids with diamond structure[END_REF][START_REF] Broido | Lattice thermal conductivity of silicon from empirical interatomic potentials[END_REF][START_REF] Kaburaki | Dynamical thermal conductivity of argon crystal[END_REF][START_REF] Turney | Predicting phonon properties and thermal conductivity from anharmonic lattice dynamics calculations and molecular dynamics simulations[END_REF][START_REF] Mcgaughey | Quantitative validation of the Boltzmann transport equation phonon thermal conductivity model under the single-mode relaxation time approximation[END_REF][START_REF] Li | ShengBTE: A solver of the Boltzmann transport equation for phonons[END_REF]. In general, all the used methods are able of describing the behavior of the thermal conductivity, when they are involving a proper formulation of the normal and umklapp phonon-phonon scattering and of the other scattering mechanisms due to the presence of defects in the crystal. If we consider a solid where only phonons are acting in the thermal transport, the sample turns out to be like a box containing them. At very low temperatures, we have a few phonons interacting only with the walls of the box. As the temperature increases, the population of phonons increases too and consequently their thermal conductivity. However, with a further increase of temperature, also the interactions among phonons become relevant. Therefore, after the thermal conductivity has reached a maximum, its value starts decreasing with a further increase of the temperature. Among the processes which are influencing the values of the thermal conductivity we find the scattering of phonons due to the presence of isotopic mass defects. Actually, the natural crystals of several elements are composed by different isotopes; that is, in these crystals with natural composition we have the sites of the crystalline lattice occupied by atoms of the same element having different masses. This disorder of the lattice reduces the thermal conductivity. For an experimental study of the thermal conductivity without the isotope effect, the sample needs to be made of an isotopically enriched material. Among the materials which have received specific enrichment for the analysis of thermal conductivity we find diamond-like Carbon, Silicon and Germanium [START_REF] Anthony | Thermal diffusivity of isotopically enriched C 12 diamond[END_REF][START_REF] Wei | Thermal conductivity of isotopically modified single crystal diamond[END_REF][START_REF] Onn | Some aspects of the thermal conductivity of isotopically enriched diamond single crystals[END_REF][START_REF] Olson | Thermal conductivity of diamond between 170 and 1200 K and the isotope effect[END_REF][START_REF] Asen-Palmer | Thermal conductivity of germanium crystals with different isotopic compositions[END_REF]. As resulting from experiments, the isotopic purity of the solid is strongly influencing the thermal conductivity. This is important for technological applications of these materials, such as for a better comprehension of the thermal transport.

Here we discuss the case of crystalline Silicon. Specifically, we discuss some of the literature on it, which is providing experimental data. Theoretical works had been proposed in the past by the author, who solved the phononic Boltzmann equation in crystals in which different isotopes are present. Let us start for these theoretical works on the thermal conductivity of Silicon.

An iterative approach

In [START_REF] Omini | An iterative approach to the phonon Boltzmann equation in the theory of thermal conductivity[END_REF][START_REF] Omini | Effect of phonon scattering by isotope impurities on the thermal conductivity of dielectric solids[END_REF] we proposed an iterative approach to the solution of the Boltzmann equation. In the calculation of the thermal conductivity, we inserted competing scattering mechanism due to isotope impurities and the three-phonon collisions, without resorting to the systematic use of the relaxation time approximation. The model that we used was not so specific as in the following works [START_REF] Omini | Beyond the isotropic-model approximation in the theory of thermal conductivity[END_REF][START_REF] Omini | Heat transport in dielectric solids with diamond structure[END_REF], but it was able to give a good agreement with experimental data for silicon, here shown in the Figure 1.

Figure 1: Thermal conductivity of Silicon as a function of temperature at room pressure. The curves A and B are the results of calculations proposed in [START_REF] Omini | Effect of phonon scattering by isotope impurities on the thermal conductivity of dielectric solids[END_REF], with and without isotopic mass defects, in the case of absence of the boundary scattering. This scattering is considered for curve B'. The experimental points are from Holland and Neuringer [START_REF] Holland | The effect of impurities on the lattice thermal conductivity of silicon[END_REF].

Let us note that the results given in the Figure 1 had not been obtained by using adjustable parameters. Actually, the use of adjustable parameters is rather common in the methods involving the relaxation time approximation. A criticism to this use was made by Srivastava too [START_REF] Srivastava | Calculation of lattice thermal conductivity of Ge from 4 to 900 K[END_REF].

In [START_REF] Omini | Effect of phonon scattering by isotope impurities on the thermal conductivity of dielectric solids[END_REF] we have also considered the isotope effect in Germanium. And for it too we observed the relevant role of the purity of the sample. For our further discussion on Silicon, let us extrapolate from the Figure 1 the thermal conductivity at 200 K. For the natural Silicon it is of 2.2 W/cm-K and for the pure sample, it is of 3.8 W/cm-K.

Following the methods given in [START_REF] Omini | Beyond the isotropic-model approximation in the theory of thermal conductivity[END_REF][START_REF] Omini | Heat transport in dielectric solids with diamond structure[END_REF], the effect of isotopes on the thermal conductivity of Silicon was investigated again in [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF]. Actually, in [START_REF] Omini | Beyond the isotropic-model approximation in the theory of thermal conductivity[END_REF], [START_REF] Omini | Heat transport in dielectric solids with diamond structure[END_REF] and [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF], we applied the iterative method on a revised model of the solid: we used a crystal described by its real Brillouin zone and Hamiltonian to determine properly the three-phonon normal and umklapp processes. Again, for the further discussion, let us extrapolate from the Figure 2 the thermal conductivity at 200 K. For the natural Silicon it is of 2. W/cm-K and for the pure sample, it is of 3. W/cm-K .

As we can see from the values obtained in [START_REF] Omini | Effect of phonon scattering by isotope impurities on the thermal conductivity of dielectric solids[END_REF] and [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF], we have that the ratio of the thermal conductivity of enriched and natural Silicon, at 200 K, is 1.72 in the case of [START_REF] Omini | Effect of phonon scattering by isotope impurities on the thermal conductivity of dielectric solids[END_REF], and it is 1.5 in the case of [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF]. It means that, using a more realistic mode, the effect of the isotopic mass defects is reduced. From [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF], we estimate that the thermal conductivity at 200 K is reduced of 30% at the most.

Experimental results

At the time we investigated the use of the iterative approach to solve the Boltzmann equation, we had not the opportunity of having a service, such as Google Scholar or other web supports, to find easily the literature on the subject. As the reader can see, we had at our disposal only the access to a few experimental data for Silicon. Now, let us examine the further data that we find on the web today. First, we consider the paper by Capinski, Maris, Bauser, et al. [START_REF] Capinski | Thermal conductivity of isotopically enriched Si[END_REF]. The experimental data are given in the Figure 3.

Fig. 3: Thermal conductivity of isotopically pure 28 Si (squares) and Si of natural isotopic abundance (triangles) [START_REF] Capinski | Thermal conductivity of isotopically enriched Si[END_REF]. The circles are the data from Ref. [START_REF] Holland | Analysis of lattice thermal conductivity[END_REF].

The data from Ref.24 tell us that, at 200 K, the thermal conductivity of the natural Silicon is of about 2.8 W/cm-K and that of the isotopically pure Silicon of about 4 W/cm-K . The ratio of the thermal conductivity is 1.42. The isotope scattering is reducing the thermal conductivity of 30%. These results seem in good agreement with the previously given theoretical data.

In 2000, new data on the thermal conductivity of Silicon had been published [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF]. Here in the Figure 4, we give some of the experimental points.

Fig. 4: Thermal conductivity of isotopically pure 28 Si and Si of natural isotopic abundance. Some data from Ref. [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF]. The open dots represent the measurements on sample SIN1 of natural silicon. As told in [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF] the "plus" symbols denote the "standard" curve for natural Si, as given by [START_REF] Touloukian | TPRC data Series[END_REF].

Let us consider the data from Ref.26 at 200 K. The thermal conductivity of the natural Silicon is of about 2.3 W/cm-K and that of the isotopically pure Silicon of about 4 W/cm-K (Figure 3). In [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF], the authors discussed the result in the following manner. "In principle, the thermal conductivity enhancement, which is most prominent at the maximum of kappa, where umklapp processes are frozen out, should vanish rapidly at higher temperatures. However, we find an almost constant enhancement in kappa of about 60% between 100 and 300 K and estimate this effect to persist up to at least 400 K. This observation also contradicts the predictions of a recent calculation for the thermal conductivity enhancement in isotopically pure silicon". The authors are referring to article [START_REF] Capinski | Analysis of the effect of isotope scattering on the thermal conductivity of crystalline silicon[END_REF].

The behavior at high temperature Data published in [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF] had been discussed in the work by Inyushkin in [START_REF] Inyushkin | Thermal Conductivity of Isotopically Modified Silicon[END_REF]. He noted, about the phonon free path at low temperatures, that a "tremendous discrepancy" between theoretical and experimental values exists, and it "cannot be rationalized in the framework of the current views of heat conduction in crystals in the boundary scattering regime". Actually, also the lack of results for sample SIN1 (natural Silicon) at low temperatures was stressed. The authors of [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF] considered Inyushkin's analysis and reconsider their results. "In view of the interest, in this problem" [START_REF] Inyushkin | Thermal Conductivity of Isotopically Modified Silicon[END_REF], "we have recently remeasured the same and other similar samples with a different apparatus around T = 300 K. We consistently found an increase in the k of 28 Si with respect of that of nat Si of about 10 %, whereas the values of this enhancement measured below 100 K are the same as those in the original data" [START_REF] Ruf | Erratum to "Thermal conductivity of isotopically enriched silicon[END_REF].

Fig. 5: The same as in the Figure 4 with the red squares representing the correction reported in [START_REF] Ruf | Erratum to "Thermal conductivity of isotopically enriched silicon[END_REF] for the data of pure 28 Si.

The authors, in the new measurements, found the thermal conductivity of the natural Silicon the same as before, but the "values reported earlier for 28 Si are systematically 90 W/m-K higher than the present ones". "This difference ... below 100 K … reduces the enhancement … reported earlier for 28 Si from 60% to about 10%. The reason for the discrepancy is not known, but we believe the new data to be more reliable". Therefore we can revise accordingly the data of the Figure 4, obtaining Figure 5.

Let us consider the value of the thermal conductivity of the enriched Silicon at 200 K. From the Figure 4, we have a value of 400 W/m-K: if we reduce it of 90 W/m-K, we obtain a value of 3 W/cm-K (Figure 5). Let us note that this value is in good agreement with the result of the theoretical model in [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF]. From the plot, the conductivity of natural Silicon seems being 2.2 W/cm-K. The ratio of the conductivity of enriched and natural Silicon seems being 1.36. Then, the enrichment increases, at 200 K, the conductivity of 36% .

Discussion

Let us consider again the results given in [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF]. In the Figure 6, we can see side by side the theoretical results given in the Ref. [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF] and the experimental data given in the Figure 5, with the corrections for the enriched sample (red squares).

Figure 6: On the left, the theoretical curves from [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF], on the right the experimental data [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF][START_REF] Ruf | Erratum to "Thermal conductivity of isotopically enriched silicon[END_REF]. The red lines had been added to help the reader in the comparison.

We can discuss the comparison. First of all, I have removed the experimental data of the enriched Silicon below 3 K. These data depend on the boundary scattering and therefore on the size of the sample. For the enriched sample, the specific value of the thermal conductivity at the peak is different, however the thermal conductivity obtained from theory between 3 K and 10 K and between 30 K and 200 K is in good agreement with experimental data. For the natural sample, above 10 K, the overall behaviors are in agreement. Moreover, the theoretical curve is in agreement with the open circles from [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF].

New data

In [START_REF] Kremer | Thermal conductivity of isotopically enriched 28 Si: revisited[END_REF], the researchers measured the thermal conductivity of isotopically enriched 28 Si "independently in three laboratories by high precision experiments on a total of four samples of different shape and degree of isotope enrichment in the range from 5 to 300 K with particular emphasis on the range near room temperature. The results obtained in the different laboratories are in good agreement with each other. They indicate that at room temperature the thermal conductivity of isotopically enriched 28 Si exceeds the thermal conductivity of Si with a natural, unmodified isotope mixture by 10±2%. This finding is in disagreement with an earlier report by Ruf et al.". The results given in [START_REF] Kremer | Thermal conductivity of isotopically enriched 28 Si: revisited[END_REF] are in agreement with the data from [START_REF] Gusev | Thermal Conductivity of 28 Si from 80 to 300 K[END_REF]. Let us note that the results for natural Silicon from [START_REF] Kremer | Thermal conductivity of isotopically enriched 28 Si: revisited[END_REF], at 200 K, are rather different from data given in [START_REF] Capinski | Thermal conductivity of isotopically enriched Si[END_REF], and also from [START_REF] Ruf | Thermal conductivity of isotopically enriched silicon[END_REF] and [START_REF] Ruf | Erratum to "Thermal conductivity of isotopically enriched silicon[END_REF]. Fig. 7: Data from [START_REF] Kremer | Thermal conductivity of isotopically enriched 28 Si: revisited[END_REF].

From the Figure 7, at 200 K, the conductivity is comprised between 2.5 and 2.9 W/cm-K. The ratio of the conductivity is 1.16. Actually, the new measurements had increased the value of the thermal conductivity in the natural sample of 16%. Let us stress this fact, that it is a good set of measurements on the natural sample which is relevant for any model of the phonon scattering by the isotopic mass defects. For what concerns the calculation in [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF], the given model based on a real representation of the crystalline lattice and of the scattering processes produces values of the thermal conductivity in good agreement for the pure crystal. For the crystal with natural isotopic mass defects, the formula used for the scattering produced by these defects is giving, at 200 K, a conductivity 25% lower than the new experimental value. In the case of Diamond and Germanium, that were studied in [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF] too, the agreement between theory and experiments is good on all the range of temperature.

Conclusion

In this paper we have presented some theoretical curves and experimental data on the thermal conductivity of crystalline Silicon. Let us stress that a good overall agreement of the theoretical curves with experimental values exists, in particular for enriched Silicon [START_REF] Sparavigna | Influence of isotope scattering on the thermal conductivity of diamond[END_REF]. The theoretical curves had been obtained from a realistic model, without the use of adjustable parameters. For what concerns the natural Silicon, the most recent data here considered give a thermal conductivity at 200 K which is about 10% larger than the previously determined value. It would be interesting to understand the origin of this difference in the experimental values.

Figure 2 :

 2 Figure 2: Thermal conductivity of Silicon as a function of temperature at room pressure. The curves A and B are the results of calculations proposed in [22], with and without isotope scattering. The experimental points are from Refs. 20 and 23.


Pierre Chalimbaud
email: chalimba@univ-bpclermont.fr

François Berry
email: berry@univ-bpclermont.fr

Embedded active vision system based on an FPGA architecture

In computer vision and more particularly in vision processing, the impressive evolution of algorithms and the emergence of new techniques dramatically increase algorithm complexity. In this paper, a novel FPGA-based architecture dedicated to active vision (and more precisely early vision) is proposed. Active vision appears as an alternative approach to deal with artificial vision problems. The central idea is to take into account the perceptual aspects of visual tasks, inspired by biological vision systems. For this reason, we propose an original approach based on a System On Programmable Chip implemented in an FPGA connected to a CMOS imager and an inertial set. With such a structure based on reprogrammable devices, this system admits a high degree of versatility and allows the implementation of parallel image processing algorithms.

Introduction

FPDs and in particular FPGAs have achieved rapid acceptance and growth over the past decade because they can be applied to a very wide range of applications [START_REF] Dehon | The density advantage of configurable computing[END_REF]. One of the most interesting applications of FPGAs is the prototyping of designs to be implemented as gate arrays. Another is the emulation of entire large hardware systems. Apart from prototyping, an emerging topic for FPGA applications is their use in custom computing machines. This involves using the programmable parts to "execute" software, rather than compiling the software for execution on a regular CPU. For the latter, the notion of such Soft-core CPU or hardware overload of the instruction set becomes crucial. Such approaches offer a good trade-off between the performance of fixedfunctionality hardware and the flexibility of softwareprogrammable substrates. These different aspects are a great advantage in the design of an embedded sensing system, in particular when there are several data flows. Like ASICs, the main benefit of these systems is their ability to implement specialized circuitry directly in hardware. However, fast prototyping is easier for FPGAs. Conse-quently, in the design of a versatile embedded system dedicated to image processing, the FPGA solution proves to be the better way.

In computer vision and especially in vision processing, the impressive evolution of algorithms and the emergence of new techniques drastically increase the complexity of algorithms. This computational aspect is crucial for the majority of real-time applications and in most cases programmable devices are the best option. For example, FP-GAs have already been used to accelerate real-time point tracking [START_REF] Benedetti | Real-time 2-d feature detection on a reconfigurable computer[END_REF], stereo-vision computing [START_REF] Woodfill | Real-time stereo vision on the parts reconfigurable computer[END_REF], color-based object detection [START_REF] Benitez | Reactive computer vision system with reconfigurable architecture[END_REF], and video and image compression [START_REF] Böhm | Mapping a single assignment programming language to reconfigurable systems[END_REF].

In this paper, an architecture dedicated to computer vision is proposed. Our approach towards a smart camera consists in performing most of the early vision processing at the sensor level, before transmitting the information to the main processing unit. This behavior is inspired by the human vision system, where eyes are responsible for attention and fixation tasks, sending to the brain only pertinent information about the observed scene. As a matter of fact, the amount of visual data to be transmitted and analyzed is strongly reduced and communication bottlenecks can be avoided. The adaptation of perceptual aspects from biological vision to artificial systems, which is known as active vision and active perception, is briefly explained in section 2 as the principal motivation of this work. Consequently, the main originality of this work is to use the concepts developed in active vision and more generally in bio-inspired computer vision in order to design suitable hardware. In the third section, the hardware of the smart camera is described. The technological choices are argued according to the objectives given in the previous section. The different modules are fully described and the different data flows are explained. The fourth section presents the core of the FPGA design, in particular the specific modules like the address generation unit or the Fixed Pattern Noise (FPN) correction unit. Finally, we present the results of two image processing algorithms (motion detection and high-speed template tracking implementation).

p. 1

One of the numerous objectives in artificial vision research is to build computer systems that analyze images automatically, determining what the computer "sees" or "recognizes" and "understands" from the environment.

In what follows, the problem is to perform the process of interpretation of sensorial data within an environmental model. The first ways of treating the "vision problem" used passive vision and dynamic vision approaches. Passive vision comprises the classical analysis of images. The approach that David Marr explicitly advocated [START_REF] Marr | San Francisco[END_REF] and to which many others subscribe, it has led to a thriving research field that has been dominant in visual science in recent years. From David Marr, "Vision is a process that produces from images of the external world a description that is useful to the viewer and not cluttered with irrelevant information". David Marr proposes a model of visual processing that begins by identifying the "zero-crossings" (edges) in the image, uses this edge information to provide a crude segmentation of surfaces called the 2-D sketch, and finally extracts from this sketch the three-dimensional spatial information. That spatial interpretation is expressed in terms of geometrical primitives such as generalized cylinders or cones, so that the only data which must be explicitly stored are the x, y, z locations, alpha, beta, gamma orientations, and aspect ratios of each of the cylinders and a symbolic code of the relations between them. In this way the complex scene is reduced to a highly compressed set of meaningful numbers. The problem with this model is that nobody has ever been able to define how such spatial information can be reliably extracted from the scene. Moreover, the visual world contains far too many ambiguities to be handled successfully. Dynamic vision is a complementary approach which corresponds to the study of visual information, but in an unbounded sequence of views. This approach introduces time into the image processing, while movement (measured by optical flow) is used in the perception process. Some classical approaches using these strategies revolve around recovering structure from motion.

In contrast to these two approaches, [START_REF] Aloimonos | Active vision[END_REF], [START_REF] Bajcsy | Active perception[END_REF] and [START_REF] Ballard | Animate vision[END_REF] have proposed the active vision approach. Active vision techniques are derived from attempts to simulate the human visual system. In human vision, head motion, saccadic eye movement, and the eye's adaptation to lighting variations are important in the perception process. Active vision therefore aims to simulate the power of this adaptation. In other words, active vision is an alternative approach to dealing with artificial vision problems. The central idea, also known as the task-driven paradigm, is to take into account the perceptual aspect of visual tasks. Therefore, instead of a full 3D representation of the observed scene, the system is supposed to extract only the information useful for solving a given problem, through a task-driven observation strategy (Figure 1). An artificial active vision system uses observer-controlled input sensors. Its main goal must be actively to extract the requested information in order to solve a given task. A wide literature proposes many systems built around the active vision paradigm. The majority of these systems have been driven by the "robotic" approach and are based on a robotic head. A large survey up to 1996 can be found in [START_REF] Andersen | A Framework for control of a camera head[END_REF] and in [START_REF] Vieville | A few steps towards 3d active vision[END_REF].

Another trend considers algorithmic aspects and focuses on gaze control using foveated sensors with a log-polar mapping. This method can be applied at the sensor level (imager), at the image processing level or both. At the sensor level, some dedicated imagers based on logarithmic-structured space-variant pixel geometry have been implemented. The main advantage of these methods is the ability to cover wide work-spaces with high acuity and a small number of pixels. Several descriptions of the advantages of using space-variant, or logmap, architectures for video sensors have been proposed [START_REF] Rojer | Design considerations for a space-variant visual sensor with complexlogarithmic geometry[END_REF][2] [START_REF] Tistarelli | On the advantage of polar and log-polar mapping for direct estimation of time-to-impact from optical flow[END_REF]. Another logmap device consists of an emulated sensor based on a conventional CCD and an image warp algorithm embedded on a microcontroller [START_REF] Bederson | A Miniature Space-Variant Active Vision System: Cortex-I[END_REF]. More recently, a new trend towards smaller active vision systems comparable in size to the human head is pushing the limit of motor, gearbox and camera design [START_REF] Kuniyoshi | Active stereo vision system with foveated wide angle lenses[END_REF][START_REF] Sharkey | On the kinematics of robot heads[END_REF][START_REF] Truong | A cabledriven pan-tilt mechanism for active vision[END_REF].

However, as mentioned above, most work dedicated to active vision systems is concentrated in the robotic field. In contrast, the main motivation for the work presented here, is to propose a system truly resulting from the hup. 2 man visual system. Consequently, our approach needs a dedicated architecture for which the FPGA proves to be essential. This architecture is presented and discussed in the next section.

Architectural features

The main purpose of our architecture is to allow the implementation of early vision processes as in the human or primate visual system. In these systems it is well known that the first neural layers (in the retina) pre-filter the visual data flow in order to select only the conspicuous information. From this pre-filtered information, an attentional processing allows focusing on the selected target. In the literature, several computational models of visual attention can be found. The first representative model was proposed by Koch and Ullman in [START_REF] Koch | Shifts in selective visual attention: Towards the underlying neural circuitry[END_REF] and has been recently revised by Itti [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF]. In these models, the purpose of the saliency map is to combine the "salient" or "conspicuous" location information from each of the lower feature maps into a global measure to determine how different a given location is from its surroundings. This technique is used to guide selective attention. The design of our active vision system is based on this kind of approach where we assume that the strategy of visual processes can be divided into three successive tasks: Attention: This is the initializing step of the process. Whole images are grabbed while waiting for the building of the saliency maps. These maps are built in parallel and represent/code conspicuity within the visual field along particular dimensions (e.g. color, orientation, or motion). The result of this step is a set of ROIs (Regions Of Interest).

Focusing: This step allows the generation of the geometry of an ROI (rectangular, tilted, foveal, circular, . . . ) and the optimization of the Signal/Noise Ratio: Contrast optimization in an ROI [START_REF] Chalimbaud | Contrast optimization in a multi-windowing image processing architecture[END_REF], Tracking of an ROI in motion, . . .

High-level processing: This last step comprises different kinds of tasks such as identification and classification.

The attention stage needs strong parallelization, on the one hand to respect real-time constraints, and on the other hand because of the intrinsic characteristics of the algorithms. As examples, some classical algorithms in an attention task used to build an efficient saliency map are motion detection, Gabor filters, and color segmentation. However, the characteristics of particular visual tasks may require dedicated image processing and only an FPGA approach allows such flexibility. For architectures such as these, a Stratix EP1S60 from Altera has been chosen. This choice is detailed below. The need for strong parallelization was what led us to connect 5×2MB SRAM synchronous memory blocks. Each 2MB memory has private data and address buses. Consequently, in the FPGA, 5 attention processes (using 2MB each) can address all the memory at the same time and an SDRAM module socket provides an extension of the memory to 64 MB (Figure 2).

The focusing stage must control the imaging devices in order to address only the ROI and to optimize the analog signal conversion. That is the reason why the sensing board has been designed around a CMOS imager and a set of 4 Digital/Analog Converters. A set of inertial sensors has been added in order to estimate the movements of the camera and improve the perception (stabilization and depth estimation [START_REF] Chalimbaud | Design of a hybrid visuo-inertial smart sensor[END_REF]).

In our approach, the high-level processing has to be performed on a host computer rather than on the embedded system. In order to send the data, the smart camera is connected via a high speed communication (USB2.0 or FireWire). The embedded system is integrated into a modular architecture consisting of three boards: the sensing board, the processing board and the communication board. An overview of the smart camera is shown in Figure 3 and a structural description presents the stacked structure with 3 boards.

System On Programmable Chip Features

As described in the previous section, the sensor was designed around a Stratix EP1S60 manufactured by Altera. This component enables a high density of integration (57120 Logic Elements). It also has three further main advantages which guided our choice.

Firstly, the Stratix is optimized to maximize the performance benefits of SOPC integration based on a NIOS embedded processor. A NIOS processor is a userconfigurable soft core processor, allowing many implementation and optimization options. The NIOS CPU is a pipelined general-purpose RISC microprocessor which supports both 32-bit and 16-bit architectural variants. Both 16 and 32-bit variants use 16-bit instructions. For p. 3 our sensor, the main advantage of this soft core processor is its extensibility and adaptability. Indeed, users can incorporate custom logic directly into the NIOS arithmetic logic unit (ALU). Furthermore, thanks to a dedicated bus (Avalon Bus), users can also connect into the SOPC onchip processor and custom peripherals. They can thus define their own instructions and processor peripherals to optimize the system for a specific application.

Secondly, the Stratix integrates DSP Blocks. These embedded DSP Blocks have been optimized to implement several DSP functions with maximum performance and minimum logic resource utilization. Each DSP block offers multipliers, adders, subtractors, accumulators and a summation unit, functions that are frequently required in typical DSP algorithms. Each DSP block can also support a variety of multiplier bit sizes (9×9, 18×18, 36×36) and operation modes (multiplication, complex multiplication, multiply-accumulation and multiply-addition) and can offer a DSP throughput of 2.8 GMACS per DSP block. The EP1S160 device has 18 DSP Blocks that can support up to 144 9×9 multipliers. These embedded DSP Blocks can also be used to create DSP algorithms and complex math routines in high-performance hardware. These can then be accessed as regular software routines or implemented as custom instructions on the NIOS CPU. For example, a cumbersome algorithm can be implemented in hardware and directly executed in software using a custom instruction. This gives designers the flexibility and portability of high-level software design, while maintaining the performance benefits of parallel hardware operations in FPGAs.

Lastly, the Stratix device incorporates a configurable internal memory called TriMatrix memory. The TriMatrix memory is composed of three sizes of embedded RAM blocks. The Stratix EP1S60 TriMatrix memory includes 574 M512 blocks (32×18 bits), 292 M4K blocks (128×36 bits), and 6 M-RAM blocks (4K×144 bits). Each of these blocks can be configured to support a wide range of features and to synthesize a wide variety of RAM (FIFO, double ports). With up to 5 Mbits of fast RAM, the TriMatrix memory structure is therefore appropriate for handling the bottlenecks arising in sensor vision algorithms.

Sensing device

This module consists of a CMOS imager manufactured by Neuricam, two 2D-accelerometers from Analog Devices and three 1D-gyrometers from Murata. The imager allows full 2D addressing with a column bus and a row bus. It has a resolution of 640×480 (VGA) and provides a broad dynamic range (120db) due to the logarithmic response of its pixel structure. Four digital-to-analog converters allow modification of the four analog voltages of the imager: analog signal offset, digital conversion range, voltage reference and a pixel precharge voltage. These four converters are used to optimize the conversion range. In effect, the CMOS imager has a logarithmic curve that enables the broad dynamic range (120dB). The major difference between biological and artificial vision systems most probably lies in their flexibility. In order to develop adaptive capacities, the hardware architecture previously described is designed to implement some specific low-level processing dedicated to early vision. This low-level processing attempts to establish an efficient interface between sensitive elements and high-level perception systems. As explained in section 3, our strategy for efficient visual perception is based on three layers: attention, focusing and high-level processing. This approach adopts a pyramidal method which reduces the amount of data flow. Typically, we can consider a simple system built around an attentional module based only on color segmentation and a focusing module based on template tracking (Figure 5) In the following, a detailed de-Figure 5: Example of perception strategy scription of the FPGA organization is presented. All the "standard modules" (FPN correction, Addressing module) are described and designs for focusing and attentional modules are proposed.

Implementation approach

The implementation of such an approach requires the management, sequentially and concurrently, of the execution of the routines previously described. Indeed, all task-oriented execution (attention, focusing and identification) is controlled by supplied results and these three layers possibly have to share areas of interest. Moreover, the information bottleneck located in the imager level should be continuously optimized to ensure high performance. In our hardware architecture, these functions are carried out by what we term a "Sequencer" (6-M0) and are performed on the NIOS soft core processor. This solution has two main advantages. Firstly, we benefit from software flexibility to define the routines' interactions, and secondly, the soft core processor allows an efficient architectural matching with the other parts of the supervision unit.

An internal RAM (Figure 6-R0) is used to store the instruction sequences which define the sequencer behavior according to the task under consideration. The host computer which uses our embedded system communicates it through a communication bus (USB 2.0 protocol) and sends requests in order indicate the sequencer the relevant behavior to adopt. More precisely, according to (and potentially set of meters) are passed in a dedicated stack, the sequencer chooses pre-established between the modules (Figure 6-P4) which constitute a dedicated processing chain. This architectural module implements the previously described routines of environmental adaptation, attention, focusing and low-level identification (Figure 6-P6). A number of these modules, due to environmental adaptation (Processing No. 1 to i), modify the pixel flow which is going to be used by the attention, focusing and identification modules (Processing No. j to j+1). The different data flows (corrected windows of interest and inertial measurements) can be used by these modules to perform computing. The set of results that are provided by these processing modules are collected in a buffer. This is how the sequencer selects results to send to the host computer. The sequencer is going to use a part of these results to perform visual feedback on sensing devices (Figure 6-S0) using dedicated control modules (Figure 6-P0, P1 and P3). We note in Figure 6 the module P2 which works the external RAM R1. This module performs the Fixed Pattern Noise Correction which is absolutely essential with the image sensor technology we use (described in section 4.3). Lastly, the dedicated communication module P5 is a multiplexer that synchronizes the corrected pixels flows and the sequencer results flows for sending to the host computer. The sequencer constitutes an active interface between the sensing device, processing chain and the host computer. The modular processing chain is synchronized with raw pixel flow control provided by the CMOS imager. Finally, this control allows dynamic control of the global sensor state according to global visual data coherence.

Addressing module (Figure 6-P0)

The goal of the address generator device is to compute line and column addresses of the current window of interest. The shape of the window is actually rectangular and is set by 5 parameters: position (X,Y), size (H,W) and orientation (α).

The address generation is based on the well-known "Bresenham" graphical algorithm [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF]. For the computation of the tilted rectangular window addresses, we have implemented an architecture based on a recent method for drawing straight lines suitable for raster-scan displays and plotters developed by C. Bond [START_REF] Bond | A New Line Drawing Algorithm : Based on Sample Rate Conversion[END_REF].

The approach proposed by Bond is based on signal processing concepts related to resampling, multirate processing and sample rate conversion. The x-coordinates of each pixel can be viewed as a uniform sample set, and the y-coordinates represent another sample set. Assuming the slope of the line is within the first octant, the y-coordinate is generated by resampling the set of x-coordinates. To control this resampling, the fractional part of the y-coordinate is stored in a control variable as an integer. The algorithm can be summarized by these few lines of code:

• Initialization step Incr=((Y2-Y1)/(X2-X1))*2^n X=X1 Y=Y1 Cvar=1/2*2^n • Loop: repeat{ Cvar+= Incr Y+= Carry X++ } until {X = X2}
(X1, Y 1) and (X2, Y 2) are the coordinates of the segment represented in Figure 7. n is the number of bits which are used to store the fractional part of the y-coordinates. Cvar is the fractional part of the ycoordinate. Incr is an integer variable used to store slope value. Carry is an overflow indicator for the operation p. 6 The extension of the algorithm to other octants is performed by interchanging the roles of x and y, and changing the signs of the coordinates x and y. The internal architecture of the address generator device is illustrated in Figure 8. The range of window tilt is encoded in a natural binary-coded variable named Angle. The first three MSB bits of this variable define the octant (Figure 8). The other bits of Angle are used by two functions, based on the Bond algorithm, to generate the fixed-sample variable F ixCoor and the resampled variable ResCoord for each window of interest dimension. According to the quadrant, a decoder defines the sign of the F ixCoor and ResCoord variables. When the line is located in octants 2, 3, 6 or 7, the decoder causes the use of the complementary angle and the inversion of F ixCoor and ResCoord. Finally, the sum of the line coordinates, column coordinates and position vector of the window on the imager give the iterative X and Y address of each pixel. The implementation of this module in the FPGA is characterized by the following parameters:

Total logic elements 150 < 1% Total memory bits 0 0% DSP block 9-bit elements 0 0 % System Clock Frequency 190 MHz Due to the technological limits, a classical CMOS imager (without an embedded CDS correction) provides a raw pixel flow with a high FPN. Indeed, the non-uniformity of the electrical characteristics of each pixel involves an additional stage of correction. In order to even out the electrical response of each pixel, the module called FPN correction module (Figure 6) subtracts the reference values (of the FPN) from the pixel flow. These referent values represent offset differences between each pixel for the same illumination (Figure 9). The set of these values constitutes a reference image. In order to carry out this correction, the sensor integrates a module in order to load it from an external RAM (Figure 6-R1). The implementation of this module in the FPGA is characterized by the following parameters: Total logic elements 65 < 1% Total memory bits 85 < 1% DSP block 9-bit elements 0 0 % System Clock Frequency 241 MHz

Example of attention module implementation: Motion detection

Based on an image difference method (Figure 10), this algorithm looks for moving objects in a scene. In the image plan, motion is transduced to temporal and spatial graylevel changes. This module detects temporal changes and defines a rectangular window around the moving object. In the first step, a difference image is obtained through subtraction of images i and i-1. Then the difference image is thresholded, and its vertical projection is calculated (line by line pixel sum in each column). A peak detector is applied to the vertical projection, giving the horizontal position of the moving objects found in the scene. The horizontal projection inside each vertical zone detected previously is then calculated, and a second peak detection is applied to define the vertical position of the moving objects. In this way, it is possible to define the position of several moving objects in the image simultaneously. This information can be used as a parameter for another algorithm, such as a tracker.

Example of focusing module: Template tracking

As explained in paragraph 2.2, we assume that an active vision process can be split into three layers: attention, focusing and interpretation. The custom module proposed This module (Figure 11) comprises two parts: a memory to store the reference template denoted I * and a dedicated architecture for the displacement estimation. The architecture adopted is based on a derivation of the Kanade-Lucas-Tomasi algorithm [START_REF] Tomasi | Detection and tracking of point features[END_REF]. This algorithm is an iterative method to estimate displacement between two frames (I * and I). The proposed method is based on the calculation of the dissimilarity between two images as follows: where N is the maximum number of iterations allowed by the processing. In our case, because of NIOS control, the iterative process is limited to 80 MHz and the pixel flow runs at 8 MHz (corresponding to N = 80 8 = 10 iterations max).

ε = W I(x + d 2 ) -I * (x - d 2 )
The architecture developed to implement this algorithm is presented in Figure 12. The first module called "storage device" allows the storage of the reference frame I * and swapping of the current pixel flow between two doubleport memories. This module performs the two functions simultaneously in order to ensure the pixel flow rate. The displacement between the reference template I * and current window I of interest is simultaneously computed during the storage of the next window of interest.

D.Gx

According to the signs and the comparison of Λ, Λ x and Λ y , the displacement counters CT x , CT y are updated as p. 9 follows :

if Λ x > Λ 4 then if sign(Λ x ) ⊕ sign(Λ) then CT x -- else CT x + + and if Λ y > Λ 4 then if sign(Λ y ) ⊕ sign(Λ) then CT y -- else CT y + +
The updating of the reference template position is carried out according to the counter values. This process is iteratively repeated and allows detection of the correct translation vector between the two frames. Lastly, the estimated translation vector is used to update the position of the window of interest in the CMOS imager in order to track the reference template.

The implementation of this architecture on a Stratix EP1S60 leads to the following parameters: To illustrate the algorithm, several images and simulation results are presented. This simulation was performed using ModelSim software with the VHDL description of our system. To evaluate the robustness of the approach, the image is artificially moved with a given displacement indicated under each image.

Conclusion

The computation of low-level vision tasks in real-time is the first and fundamental step in building an interactive vision system. This paper proposes an alternative to classical architectures with a highly versatile architecture dedicated to early image processing. The proposed embedded system attempts to define a global coordination between sensitive elements, low-level processing and visual tasks.

The approach, based on FPGA Technology and a CMOS imager, reduces the classical bottleneck between sensor and processing. The FPGA component ensures a high interaction rate between the CMOS imager and lowlevel processing. This interaction is used to select useful information earlier in the acquisition chain than for more traditional systems. It then focuses processing resources. This capacity is used to control the sensor state according to the visual task and the environment evolution. Our implementation of the FPGA and CMOS imager technologies results in high-speed vision, realenvironment vision and the efficient design of embedded systems. Among prospective algorithm candidates, we can cite the works performed on dynamically reconfigurable components such as the ARDOISE 1 project [START_REF] Demigny | Vers une méthodologie pour la programmation des architectures à reconfiguration dynamique[END_REF]. This evolution of FPGA technology seems to be attractive for performing dynamic control of the acquisition chain. Rather than having a control system state, the system itself can be physically changed and giving a higher level of suitability for many algorithms. We have also developed a DSP board in order to improve the computation capabilities of our system. With this device, our embedded system will evolve into a heterogeneous architecture and new research into co-design between the FPGA and the DSP will be necessary. Moreover, to test the validity of our approach, several visual tasks will be implemented. Our objective is to identify elementary functions in order to define a library of architectural modules. Of course, this library will provide efficient solutions to attention resolution, focusing and identification of sub-tasks according to specific applications. Finally, we plan to work on the development of software tools to facilitate the implementation of complex vision tasks

Figure 1 :

 1 Figure 1: Saccadic eye movements and task-driven strategy: Examples of eye-scanning records obtained by Yarbus [1]. Observers were given different instructions while viewing the picture "They did not expect him" by Ilya Repin. Each of the traces shows a three minute record of eye scanning with the following instructions: a-Free examination, b-Following request to give the ages of the people, c-Remember the position of the people and objects in the room.

Figure 2 :

 2 Figure 2: Architecture of the sensor

Figure 3 :

 3 Figure 3: a -Front-view of the camera, b-Back-view of the camera, c -Structural description

Figure 4 :

 4 Figure 4: Sensing board (a) Front view (b) back view-The white circles show the 4 Digital-Analog Converters. (c) Orthogonal board. (d) Global synoptic of the sensing boards.

Figure 6 :

 6 Figure 6: Block diagram of architecture adopted

Figure 7 :

 7 Figure 7: Sampled tilted Line and Octants encoding

Figure 8 :

 8 Figure 8: Synoptic diagram of Addressing module

Figure 9 :

 9 Figure 9: Images (a) without and (b) with the subtraction of the image reference

Figure 10 :

 10 Figure 10: Motion detection processus

Figure 11 :

 11 Figure 11: Template tracking module

2 w 2 W

 22 (x)dx where x = [x y] T , the displacement d = [d x d ] T ∂x , and the weighting function w(x) is usually set to the constant 1. To find the d, set ∂d to zero. consider Taylor series expansion I and I * respectively about x -d 2 and x + d 2 , we obtain: ∂ε ∂d = W I(x) -I * (x) + 1 2 g T (x)d g(x)w(x)dx = 0 where g = ∂ ∂x (I + I * ) ∂ ∂y (I + I * ) . Finally, the displacement d can be estimated by solving the equation: Zd = e where Z is the following 2x2 matrix Z = W g(x)g T (x)w(x)dx and e is the following 2x1 vector: e = [I * (x) -I(x)] g(x)w(x)dx Due to the Taylor expansion, this algorithm is not exact and needs iterations to find the correct displacement. The estimated displacement at each iteration i is denoted s i = (s x s y ) T and is calculated by: s = 2.Z -1 .e where e = W D.Gx W D.Gy t and Z = W Gx.Gx W Gx.Gy W Gx.Gy W Gy.Gy . where G x = ∂(I+I * ) ∂x , G y = ∂(I+I * ) ∂y and D = (I * -I) is the interframe difference. The displacement d = (dx dy) T is the result of the iterative process such that d = i=N i=1 s i

Figure 12 :

 12 Figure 12: Implemented architecture

W

Figure 13 :

 13 Figure 13: Results of the template tracking architecture.