E Gobet

J G López-Salas

C Vázquez

Quasi-Regression Monte-Carlo scheme for semi-linear PDEs and BSDEs with large scale parallelization on GPUs

Keywords: BSDEs, semi-linear PDEs, dynamic programming equations, empirical regressions, parallel computing, GPUs, CUDA Mathematics Subject Classification (2010) 49L20, 62Jxx, 65C30, 93E24, 68W10

In this article we design a novel quasi-regression Monte Carlo algorithm in order to approximate the solution of discrete time backward stochastic differential equations (BSDEs), and we analyze the convergence of the proposed method. The algorithm also approximates the solution to the related semi-linear parabolic partial differential equation (PDE) obtained through the well known Feynman-Kac representation. For the sake of enriching the algorithm with high order convergence a weighted approximation of the solution is computed and appropriate conditions on the parameters of the method are inferred. With the challenge of tackling problems in high dimensions we propose suitable projections of the solution and efficient parallelizations of the algorithm taking advantage of powerful many core processors such as graphics processing units (GPUs).

Introduction

The problem. In this work we are interested in numerically approximating the solution (X, Y, Z) of a decoupled forward-backward stochastic differential equation (FBSDE in short)

Y t = g(X T ) + T t f (s, X s , Y s , Z s )ds - T t Z s dW s , (1.1)
X t = x + t 0 b(s, X s )ds + t 0 σ(s, X s )dW s . (1.
2)

The terminal T > 0 is fixed. These equations are considered in a filtered probability space (Ω, F, P, (F t ) 0≤t≤T ) supporting a q ≥ 1 dimensional Brownian motion W . The filtration is assumed to be the natural filtration of W augmented with the P-null sets. In this representation, X is a d-dimensional adapted continuous process (called the forward component), Y is a scalar adapted continuous process and Z is a q-dimensional progressively measurable process. Consideration of multidimensional Y is quite straightforward and is discarded for the sake of simplicity. Regarding terminology, g(X T ) is called terminal condition and f the driver. Under standard assumptions (like L 2 -integrability and Lipschitz regularity), there is a unique solution (X, Y, Z) (see [START_REF] Pardoux | Stochastic Differential Equations, Backward SDEs, Partial Differential Equations[END_REF] for broad account on existence/uniqueness assumptions). Detailed hypotheses for the current study are exposed later. Our aim in this article is to design a new numerical scheme for approximating the solution to the above BSDE. When f is linear in (Y, Z), explicit representations are available [START_REF] El Karoui | Backward stochastic differential equations in finance[END_REF]Theorem 1.1], which can be easily turned into simple and efficient Monte-Carlo algorithms. Our focus is thus more on the case of non-linear driver yielding non-linear equations.

Feynman-Kac representations. Such equations, first studied by Pardoux and Peng [START_REF] Pardoux | Adapted solution of a backward stochastic differential equation[END_REF] in a general form, are important tools in semi-linear PDEs and stochastic control theory. Let us explain. Denote by A the infinitesimal generator of X:

A = i b i (t, x)∂ x i + 1 2 i,j
[σσ ] i,j (t, x)∂ 2 x i ,x j , where denotes the transpose operator. Consider the parabolic partial differential equation (PDE) of the form ∂ t u(t, x) + Au(t, x) + f (t, x, u(t, x), (∇ x uσ)(t, x)) = 0, t < T, u(T, .) = g(.).

Assuming that u is smooth enough, the Itô formula leads to u(t, X t ) = u(0, x) + t 0 (∂ t u(s, X s ) + Au(s, X s ))ds + t 0 (∇ x uσ)(s, X s )dW s .

Using the PDE (1.3) and noting u(T, X T ) = g(X T ), we get u(t, X t ) = g(X T )+ T t f (s, X s , u(s, X s ), (∇ x uσ)(s, X s ))ds-T t (∇ x uσ)(s, X s )dW s , therefore proving that the solution (Y, Z) can be represented as (Y t , Z t ) = (u(t, X t ), (∇ x uσ)(t, X t )). This is the well-known Feynman-Kac representation, making connection between semi-linear PDEs and BSDEs. One may prefer one or the other representations depending on the application at hand, on the numerical scheme to use, and on the theoretical estimates available (of probabilistic or PDE nature). As a consequence, solving numerically (Y, Z) for any realization of X or computing (u, ∇u) in the full space R d gives, in both cases, essentially access to the same solution. The above reasoning holds for Cauchy boundary conditions and parabolic PDEs as presented above, however extensions in the directions of various boundary conditions (Dirichlet, Neumann, free-boundary) or stationary problems (elliptic PDEs) are also available, see [39, Sections 3.8, 5.4, 5.6, 5.7 and 5.8] and references therein.

In view of (1.3), the link with stochastic control appears clearly, by analogy with the Hamilton-Jacobi-Bellman equation. Since the non-linearity in f is related only to the first derivative of u (the Z component), the underlying control problem is associated to a control in the drift only (not in the diffusion coefficient). See [START_REF] El Karoui | Backward stochastic differential equations in finance[END_REF][START_REF] El Karoui | Backward stochastic differential equations and applications[END_REF] for some examples in finance.

Applications. We now conclude the presentation of connection between BSDEs and PDEs by giving a broader overview of applications where these semi-linear PDEs come into play. We follow the presentation of [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF]Section 7.1], some examples can be found in [25, Chapter 2] and [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF]Chapter 14].

Ecology. In a region of the plane (d = 2), suppose the existence of 2 interacting species in the same environment and denote u k (t, x) (k = 1, 2) the population density of the individuals of each species at a point of the space x and at a time t. These population densities obey to the system      ∂ t u 1 (t, x) = α 1 ∆u 1 (t, x) + f 1 (t, x, u 1 , u 2 ), ∂ t u 2 (t, x) = α 2 ∆u 2 (t, x) + f 2 (t, x, u 1 , u 2 ), t > 0, x ∈ R 2 , u i (0, •) given, (1.4) where the two Laplacian operators model the diffusive movement (as a Brownian motion) of each species. Note that the initial condition at t = 0 can be turned into a terminal condition at t = T as in our previous setting, this is just a simple time reversal t ↔ T -t.

The function f = (f 1 , f 2 ) plays the role of a growth rate of the population, describing the available local resources and interactions between the species: the case ∂ u 2 f 1 < 0 and ∂ u 1 f 2 > 0 corresponds to the predator-prey model, where the growth rate of the species 1 (prey) decreases in the case of high density of the species 2 (predator) and conversely for the growth rate of predators; the case ∂ u 2 f 1 > 0 and ∂ u 1 f 2 > 0 models symbiosis, where each species benefits from the other; ∂ u 2 f 1 < 0 and ∂ u 1 f 2 < 0 describes competition between species. This kind of model, introduced in the 30s by Fisher, Kolmogorov, Petrovsky and Piskunov, is discussed with extra references in [START_REF]Biological invasions : theory and practice[END_REF]. Chemistry. Suppose that a container contains N chemical compounds taking part in R independent reactions. Denote c i the concentration of the i-th compound and θ the temperature. Then their evolution follows the N + 1 equations (in

R + × R 3 )      ε p ∂ t c i = D i ∆c i + R
j=1 ν ij g j (c 1 , . . . , c N , θ), i = 1, . . . , N, ρc p ∂ t θ = k∆θ -R j=1 N i=1 ν ij H i g j (c 1 , . . . , c N , θ), c i (0, •) and θ(0, •) given, (1.5) where g j is the speed of the j-th reaction and H i is the partial molar enthalpy of the i-th compound. For more details, see [START_REF] Gavalas | Nonlinear differential equations of chemically reacting systems[END_REF]. For reaction-diffusion equations, see [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF].

Neuroscience. The famous model of Hodgkin and Huxley is a set of equations describing the psychological phenomenon of signal transmission in the axon (nerve fiber), showing the dependences between electrical excitability and various chemical ion concentrations. The PDE system is of size 4, the first unknown u represents the electric potential and the three other unknowns (v 1 , v 2 , v 3 ) are chemical concentrations: the system writes for (t, x) ∈ R + × R (linear neuron)

               c 0 ∂ t u = 1 R ∂ 2 xx u + κ 1 v 3 1 v 2 (c 1 -u) + κ 2 v 4 3 (c 2 -u) + κ 3 (c 3 -u), ∂ t v 1 = ε 1 ∂ 2 xx v 1 + g 1 (u)(h 1 (u) -v 1 ), ∂ t v 2 = ε 2 ∂ 2 xx v 2 + g 2 (u)(h 2 (u) -v 2 ), ∂ t v 3 = ε 3 ∂ 2 xx v 3 + g 3 (u)(h 3 (u) -v 3
), u(0, •) and v i (0, •) given, (1.6) for different positive constants c i , R, κ i , ε i and different functions g i , h i . There also exists a simplification of this model, known as the FitzHugh-Nagumo model. Materials physics. The Allen-Cahn equation is a prototype model of phase transition with a diffusive interface, used to model, for example, a solid/liquid phase transition [START_REF] Visintin | Models of phase transitions[END_REF]. The system is one-dimensional, and takes the form (in R + ×R 3 ) ∂ t u = ε∆u + u(1 -u 2 ), u(0, •) given.

(1.7)

The solution u thus represents an order parameter defining the arrangement of atoms in a crystal lattice. Economy and finance. There are numerous applications of BSDEs (and semilinear PDEs) in finance, we refer the reader to [START_REF] Cvitanic | Hedging options for a large investor and Forward-Backward SDEs[END_REF][START_REF] El Karoui | Backward stochastic differential equations in finance[END_REF][START_REF] Crepey | Financial modeling: a Backward Stochastic Differential Equations perspective[END_REF]. Let us mention: valuation/hedging of contingent claim in complete markets; some market imperfections, such as higher interest rate for borrowing [START_REF] Bergman | Option pricing with differential interest rates[END_REF]; incomplete markets [START_REF] Pham | Continuous-time stochastic control and optimization with financial applications[END_REF] or super-replication [START_REF] El Karoui | Dynamic programming and pricing of contingent claims in an incomplete market[END_REF]; connections with recursive utilities in [START_REF] Duffie | Stochastic differential utility[END_REF]; non-linear pricing rule [START_REF] Peng | Nonlinear expectations, nonlinear evaluations and risk measures[END_REF], etc. See [START_REF] El Karoui | Backward stochastic differential equations and applications[END_REF][START_REF] Crepey | Financial modeling: a Backward Stochastic Differential Equations perspective[END_REF] for extra references.

Numerical approximations. In the two past decades, there have been several attempts to design efficient approximation schemes for BSDEs. It has been initiated by [START_REF] Bally | Approximation scheme for solutions of BSDE[END_REF][START_REF] Chevance | Numerical methods for backward stochastic differential equations[END_REF] and at that time, only time-discretization issue was considered. Since then, optimal discretization strategies have been derived: see [START_REF] Geiss | Generalized fractional smoothness and Lp-variation of BSDEs with non-Lipschitz terminal condition[END_REF] and references therein.

If the time grid is {t 0 := 0, . . . , t i , . . . , t N =: T }, it leads to a nested expectation problem of the form (also called Dynamic Programming Equation, or DPE for short)

Z t i ≈ E   g(X T ) + N -1 j=i+1 f (t j , X t j , Y t j+1 , Z t j )(t j+1 -t j ) (W t i+1 -W t i ) t i+1 -t i | F t i   , (1.8)
Y t i ≈ E   g(X T ) + N -1 j=i f (t j , X t j , Y t j+1 , Z t j )(t j+1 -t j ) | F t i   (1.9)
for i = N -1, . . . , 0. We use the symbol ≈ to emphasize the heuristics, a rigorous justification can be done using the aforementioned references. The above backward dynamic programming equation makes appear clearly that at each of N grid times, one has to compute a conditional expectation function, that is a function from R d to R, this is in a way equivalent to solving the PDE at each time point t i . This is a difficult problem, with increasing difficulty as the dimension gets large and/or as the unknown function gets less smooth (curse of dimensionality), like for finite differences schemes for PDEs. Among the numerous approaches, let us mention the quantization method of Bally and Pagès [START_REF] Bally | A quantization algorithm for solving multi-dimensional optimal stopping problems[END_REF], which is an optimal space discretization of the underlying dynamic programming equation, suitable to compute easily the conditional expectations on the space grid; in [START_REF] Bouchard | Discrete time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF] the authors compute the regression function with Malliavin calculus integration by parts combined with Monte Carlo simulations; the use of empirical regression methods and machine learning techniques is developed in [START_REF] Gobet | A regression-based Monte-Carlo method to solve backward stochastic differential equations[END_REF][START_REF] Lemor | Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations[END_REF] (also called regression Monte-Carlo method ); branching processes for small time or small linearity are used in [START_REF] Henry-Labordere | A numerical algorithm for a class of BSDEs via the branching process[END_REF]. Optimal convergence rates with regression methods are derived in [START_REF] Gobet | Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions[END_REF].

Parallel computations. Most of the above methods do not allow parallel computations because of the non-linearity of the equations. Recently in [START_REF] Gobet | Stratified Regression Monte-Carlo Scheme for Semilinear PDEs and BSDEs with Large Scale Parallelization on GPUs[END_REF], the authors have proposed a novel approach inspired from regression Monte-Carlo schemes suitable for massive parallelization, with excellent results on GPUs. It takes advantage of local approximations of the regression function in different strata of the state space R d (partitioning estimates) and of a new stratified sampling scheme to guarantee that enough simulations are available in each stratum to compute accurately the regression functions. The parallelization step is made with respect to the strata. Our current work lies in this vein and is aimed at overcoming some drawbacks of the stratified regression scheme of [START_REF] Gobet | Stratified Regression Monte-Carlo Scheme for Semilinear PDEs and BSDEs with Large Scale Parallelization on GPUs[END_REF]. In the previous scheme, the solution u to the PDE is represented as a piecewise polynomial function (like in a discontinuous Galerkin method [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]) and therefore, even if the exact solution is continuous, the approximation may be not, especially at the interface of the different strata. In the current work, we propose a different scheme, where the approximated solution is decomposed on basis functions that are smooth and that have natural orthonormal properties with respect to the sampling distribution, thus avoiding the possible discontinuity at the interfaces. Moreover, the sampling distribution is required to have some nice stability properties, useful for the propagation of error in the DPE , see Theorem 2; we prove that the Student's t-distribution meets this requirement, see Theorem 1. One possible set of smooth and orthonormal functions is given by the trigonometric basis composed with a space transformation related to the Cumulative Distribution Function (CDF for short) of the Student's t-distribution. We call such a basis the Student-cosine basis. Instead of a regression Monte Carlo scheme, we design a quasi-regression Monte-Carlo scheme: we take advantage of the fact that the basis functions enjoy orthonormal properties, in order to compute the basis coefficients as simple Monte-Carlo averages, instead of solutions of least-squares problems. On the one hand, it leads to a simplified and quicker procedure (because we avoid inverting linear systems). On the other hand, we obtain a globally smooth approximation which is mainly suitable for smooth solutions u [START_REF] Funaro | Polynomial Approximation of Differential Equations[END_REF][START_REF] Canuto | Spectral Methods, Fundamentals in Single Domains[END_REF]: this is a significant difference with partitioning estimates which can approximate well locally (ir-)regular functions. In this new scheme, the notion of strata is lost and the parallelization is now made with respect to the numerous required simulations and basis coefficients. Moreover, we establish tight non-asymptotic error estimates (as a function of the size of the function basis and the number of Monte-Carlo simulations).

In recent years, GPUs have become increasingly attractive for general purpose parallel computations. GPUs are not only impressive graphics tools but also highly parallel programmable processors that can deliver order-of-magnitude performance gains over optimized CPU applications. CPU hardware is designed to handle complex tasks, instead, GPUs only do one thing well, they handle millions of repetitive low level tasks. These high performance computing devices offer an arithmetic capability and a memory bandwidth that substantially outpace their CPU counterpart. In fact, a single NVIDIA Titan V GPU can sustain over 13.8 tera-floating-point operations per second (Tflops) and a memory bandwidth of 652.8 GB/s. In the last years GPUs have been heavily used in scientific computation in order to achieve dramatic runtime reductions over traditional CPU codes. Particularly successful applications of GPUs include Monte Carlo methods like the one we target in the present work. Indeed, Monte Carlo algorithms follow all design principles of GPU computing, they give rise to embarrassingly parallel not divergent codes that can exploit the GPU arithmetic intensity and its streaming memory bandwidth. In the last part of this article our goal is to offer an efficient implementation of the proposed algorithm taking advantage of several GPUs, which will allow us to address the solution of semi-linear PDEs in high dimensions.

Organization of the paper. In Section 2, we set the framework of our study, define some notations used throughout the paper and describe our algorithm based on the approximation of conditional expectations by a projection on a finite orthonormal basis of functions. Compared to the usual DPE to solve (2.1), we propose an equivalent modification that includes a weighted factor: this weight is related to the sampling distribution and its crucial purpose is to lead to higher order approximation of the unknown function under the condition it is smooth (see the Fourier-based analysis in Section 2.4). Proofs of the propagation of error are postponed to Section 3.2. This new numerical scheme has been tested on GPU, which is exposed in Section 4. Numerical results are collected in Section 5.

Notations of the paper.

(i) log(x) stands for the natural logarithm of x ∈ R + .

(ii) |x| stands for the Euclidean norm of the vector x = (x 1 , . . . , x n ) ∈ R n . (iii) x stands for the l 1 -norm of the vector x = (x 1 , . . . , x n ) ∈ R n , i.e.

x = x 1 = n i=1 |x i |. It will be mainly used when x is a multi-index of integers used for differentiation. (iv) (x) + means the positive part of x ∈ R defined by (x) + = max(x, 0). (v) Re (z) denotes the real part of the complex z ∈ C; we set i := √ -1. (vi) For any finite L > 0 and x = (x 1 , . . . , x n ) ∈ R n , define the soft truncation function

T L (x) := (-L ∨ x 1 ∧ L, . . . , -L ∨ x n ∧ L).
Moreover, if L : R n → R + is a measurable function, then the above definition of T L is extended by replacing L by L(x):

T L (x) := (-L(x) ∨ x 1 ∧ L(x), . . . , -L(x) ∨ x n ∧ L(x)). (1.10) (vii) Let j, k ∈ N d be two multi-indices. j ≤ k means that j l ≤ k l for all l = 1, . . . , d. (viii) The finite multi-indices set Λ ⊂ N d is downward closed if (k ∈ Λ and j ≤ k) ⇒ j ∈ Λ.
(ix) Let Λ ⊂ N d be an arbitrary downward closed multi-index set and f : R d → R.

f is Λ-smooth if ∂ n x f exist and are continuous for all n ∈ Λ.

2 Mathematical framework and basic properties 2.1 Model

Dynamic Programming Equation

From now on, we aim at solving the dynamic programming equation of the form (1.8)-(1.9). First, for the sake of simplicity, we intentionally consider the case where the non-linearity does not depend on the gradient (the Z component). Second, we consider the case of N equidistant grid times of [0, T ], i.e. we set t k = k∆ with ∆ = T /N . Third, the forward component X is approximated by an Euler scheme with time-step ∆, see (2.6) below. These two last approximations are known to be accurate enough when data (f, g, b, σ) are sufficiently regular (see [START_REF] Geiss | Generalized fractional smoothness and Lp-variation of BSDEs with non-Lipschitz terminal condition[END_REF] for a precise statement). To simplify notation, we write X i and Y i for an approximation of X t i and Y t i . All in all, our aim now is to solve

Y i = E   g(X N ) + N -1 j=i f j (X j , Y j+1 )∆ | F t i   for i ∈ {N -1, . . . , 0}, (2.1)
where f j (x, y) := f (t j , x, y), f being the driver in (1.1). In other words, our subsequent scheme will approximate the solutions to

X t = x + t 0 b(s, X s )ds + t 0 σ(s, X s )dW s , Y t = E g(X T ) + T t f (s, X s , Y s )ds | F t , (2.2)
and ∂ t u(t, x) + Au(t, x) + f (t, x, u(t, x)) = 0 for t < T and u(T, .) = g(.).

(2.3)

Our assumptions on the functions b, σ, f, g are as follows.

(A f ,g ) g is a polynomially bounded measurable function from R d to R: there exist two finite non-negative constants C g and η g such that

|g(x)| ≤ C g (1 + |x| 2 ) η g /2 , ∀x ∈ R d .
For every i < N , f i (x, y) is a measurable function R d × R to R, and there exist three finite non-negative constants L f , C f , η f such that, for every i < N ,

|f i (x, y) -f i (x, y )| ≤ L f |y -y |, ∀(x, y, y ) ∈ R d × (R) 2 , |f i (x, 0)| ≤ C f (1 + |x| 2 ) η f /2 , ∀x ∈ R d . (A b,σ ) The drift coefficient b : [0, T ] × R d → R d and the diffusion coefficient σ : [0, T ] × R d → R d ⊗ R q are bounded measurable functions.
One important observation is that, due to the Markov property of the Euler scheme, for every i, there exist measurable deterministic functions

y i : R d → R, such that Y i = y i (X i ),
almost surely. See the justification in [START_REF] Gobet | Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions[END_REF]Section 4.1]. Moreover, we have an uniform absolute bound for the function y i (•), the proof is postponed in Section 3.3.

Proposition 1 Assume (A f ,g ) and (A b,σ ). Then, for any x ∈ R d and any 0 ≤ i ≤ N , we have

|y i (x)| ≤ C η (C g + T C f )e C η L f T (1 + |x| 2 ) ηg ∨η f 2 =: L (x) (2.4)
with

C η := sup 0≤i≤j≤N,x∈R d E (1 + |X j | 2 ) ηg ∨η f 2 | X i = x (1 + |x| 2 ) ηg ∨η f 2 < +∞. ( 2

.5)

A second crucial observation is that the value functions y i (•) are independent of how we initialize the forward component. Our subsequent algorithm takes advantage of this observation. For instance, let X i i be a random variable in R d with some distribution ν (more details on this to follow) and let X i j be the Euler scheme evolution of X j starting from X i ; it writes

X i j+1 = X i j + b(t j , X i j )∆ + σ(t j , X i j )(W t j+1 -W t j ), j ≥ i. (2.6)
This flexibility property w.r.t. the initialization then writes

y i (X i i ) := E   g(X i N ) + N -1 j=i f j X i j , y j+1 (X i j+1 ) ∆ | X i i   . (2.7)
Approximating the solution to (2.1) is actually achieved by approximating the functions y i (•). In this way, we are directly approximating the solution to the semi-linear PDE (2.3). For some reasons that will appear clear later, we define a weighted modification of y i by

y (q) i (x) = y i (x) (1 + |x| 2 ) q/2
(2.8)

for a damping exponent q ≥ 0. For q = 0, y (q) i and y i coincide. The previous DPE (2.7) becomes

y (q) i (X i i ) := E g(X i N ) (1 + |X i i | 2 ) q/2
(2.9)

+ N -1 j=i f j X i j , y (q) j+1 (X i j+1 )(1 + |X i j+1 | 2 ) q/2 (1 + |X i i | 2 ) q/2 ∆ | X i i .
The introduction of the polynomial factor (1 + |X i i | 2 ) q/2 gives higher flexibility in the error analysis, it ensures that y (q) i decreases faster at infinity, which will provide nicer estimates on the approximation error when dealing with Fourier-basis, see Section 2.4.

Randomizing the initial value X i i

For the subsequent algorithm, X i will be sampled randomly according to the distribution ν, which we list the properties below.

(A ν ) The distribution ν(dx) of X i i has a density with respect to the Lebesgue measure on R d , of the form

ν(x) = d l=1 ν l (x l ), x = (x 1 , . . . , x d ) ∈ R d , (2.10)
for some probability density functions ν l , which are assumed to be strictly positive.

The function

x ∈ R d → L (x) (1 + |x| 2 ) q/2
(where L is defined in (2.4)) is squareintegrable on R d with respect to ν(x)dx. We denote by F ν l the cumulative distribution function (CDF) of the onedimensional marginal of ν l ,

F ν l (x) := x -∞ ν l (x )dx ,
and by F -1 ν l its inverse (which is standardly defined since F ν l is continuous and strictly increasing). Note that the above integrability condition ensures that y (q) i (for any i) is square integrable on R d with respect to ν(x)dx. In this work, we deal essentially with one type of distribution.

Example 1 (Student's t-distribution) We deal with the non-standardized Student's t-distribution on R d with parameter µ > 0 defined by its probability density function of the form (2.10) with

ν l (x l ) = c µ (1 + x 2 l ) (µ+1)/2 , (2.11)
where c µ ∈ R + . This is the density of X/ √ µ where X is a standard Student's t-distribution with µ degrees of freedom. It leads to numerous ways to simulate ν l , as acceptance rejection or ratio-of-uniforms methods, for example. From [27, p. 507] we know that

c µ = Γ ((µ + 1)/2) √ µ Γ (µ/2) √ πµ = Γ ((µ + 1)/2) Γ (µ/2) √ π ,
where Γ is the Gamma function defined in [27, p. 501]. Note that the marginal CDF and its inverse are not explicit in general for the Student's t-distribution, however their values (for usual parameter µ) are tabulated in most of books in statistics and computable using standard software libraries. Moreover, a priori estimates on F ν l , F -1 ν l and d n u F -1 ν l (u), which will be used in the error analysis of the algorithm, are computed in the Appendix A.

Note that the square integrability of

x ∈ R d → L (x) (1 + |x| 2 ) q/2 (see Assumption (A ν )) is guaranteed by the condition µ > η g ∨ η f -q.
The distribution we have considered so far has mean zero, this is suitable for problems where 0 has a natural interpretation as the center in the application at hand. In order to solve the problem far away from the neighborhood of zero we should allow non zero means (say x ∈ R d ): this is easy to handle by shifting the above distribution by x. For instance, in the case of the Student's t-distribution (Example 1), we would take

ν l (x l ) = c µ (1 + (x l -xl ) 2 ) -(µ+1)/2 .
(2.12)

Since this centering adjustment is quite obvious, we mainly discuss the case of Example 1 in what follows.

More importantly, the probability measure ν(dx) is expected to satisfy another assumption stated below.

H ρ (q) ∃c ν > 0 (depending on q, b, σ, ν, T ) such that for any measurable function

γ : R d → R with γ ∈ L 2 ν (R d ), any N ≥ j ≥ i ≥ 0, we have E |γ(X i j )| 2 (1 + |X i j | 2 ) q (1 + |X i i | 2 ) q ≤ c ν E |γ(X i i )| 2 . (2.13)
This is a stability property about expectation of weighted functionals of the Markov chain (X i j ) i≤j≤N . This property will be crucial for analyzing the propagation of error in the quasi-regression scheme (see the proof in Section 3.2). We now illustrate this property in relation with previous Example 1, the proof is postponed to Section 3.1.

Theorem 1 (Norm-stability) Assume that X is the Euler scheme

(2.6) with time step ∆ of a SDE (2.2) which drift coefficient b : [0, T ] × R d → R d is bounded measurable, and which diffusion coefficient σ : [0, T ] × R d → R d ⊗ R q
is bounded measurable and η-Hölder continuous in space uniformly in time (for some η ∈ (0, 1]),

sup 0≤t≤T sup x =x |σ(t, x) -σ(t, x )| |x -x | η < +∞.
In addition, we assume that σ fulfills a uniform ellipticity condition:

inf 0≤t≤T,x∈R d inf ζ∈R d s.t.|ζ|=1 ζ.σσ (t, x)ζ > 0.
Let X i i be sampled from distribution ν in the case of Student's t-distribution (Example 1): then H ρ (q) holds for any q ≥ 0.

Approximation spaces and basis functions

The aim of this section is to define some proper basis functions which satisfy orthogonality properties in R d and which span some L 2 space, see Proposition 2. It turns out that the choice of measure for defining the L 2 space has to coincide with the sampling measure of X i i defined in (A ν ): this plays an important role for designing Algorithm 1. Our strategy for defining such basis functions is to start from trigonometric basis on [0, 1] d and then to apply appropriate changes of variable: later, this transform will allow to easily quantify the approximation error when truncating the basis, see Propositions 3 and 4.

Notation. We start by setting notations to deal with appropriate L 2 space on R d or [0, 1] d .

(i) Let ν : R d → R be a probability density satisfying (A ν ). We will consider the space of measurable functions

L 2 ν (R d ) := v : R d → R such that R d v 2 νdx < ∞ , L 2 U ([0, 1] d ) := w : [0, 1] d → R such that [0,1] d w 2 du < ∞ .
In these spaces we define the inner products,

(v 1 , v 2 ) L 2 ν (R d ) = R d v 1 v 2 νdx, ∀v 1 , v 2 ∈ L 2 ν (R d ), (w 1 , w 2 ) L 2 U ([0,1] d ) = [0,1] d w 1 w 2 du, ∀w 1 , w 2 ∈ L 2 U ([0, 1] d ),
and the norms,

v L 2 ν (R d ) = (v, v) L 2 ν (R d ) = R d v 2 νdx 1 2 , ∀v ∈ L 2 ν (R d ), w L 2 U ([0,1] d ) = (w, w) L 2 U ([0,1] d ) = [0,1] d w 2 du 1 2 , ∀w ∈ L 2 U ([0, 1] d ).
Before the general multidimensional case, we start with an easy computation in dimension 1.

Lemma 1 (Orthonormal basis functions in dimension 1) Assume that ν satisfies

(A ν ) in dimension d = 1. Let x ∈ R and k ∈ N. The basis functions φ k (x) :=    1, if k = 0, √ 2 cos (kπF ν 1 (x)) , if k > 0, (2.14)
are orthonormal in R with respect to the probability density weight ν = ν 1 .

Proof Let i, j ∈ N. It is sufficient to see that ∞ -∞ cos (iπF ν 1 (x)) cos (jπF ν 1 (x)) ν 1 (x)dx = 1 0 cos(iπu) cos(jπu)du,
where we have used the change of variable u = F ν 1 (x). Clearly the cosine basis functions

C k (u) =    1, if k = 0, √ 2 cos(kπu), if k > 0, (2.15)
are orthonormal in [0, 1] with respect to the Lebesgue measure, therefore

(φ i , φ j ) L 2 ν (R) = 1 i=j . Note that (C k • F ν 1 )(x) = φ k (x) and (φ k • F -1 ν 1 )(u) = C k (u)
. This change of variable will be useful in order to analyze the truncation error of ν 1 series in terms of the truncation error of Fourier cosine series (see Section 2.4). Now that the reader gets acquainted with the one-dimensional case, we move to the more general result for any dimension d ≥ 1. First, the natural way to build an orthonormal multidimensional set of basis functions is to take tensor products of orthonormal one-dimensional basis functions. To this end we need to add the superscript (l) to the notations of the basis functions in dimension 1. Therefore, from now on, φ (l)

k l and C (l)
k l stand for the one-dimensional basis functions (2.14) and (2.15) in the l-th coordinate, respectively. For any multiindex k = (k 1 , . . . , k d ), set

φ k (x) := d l=1 φ (l) k l (x l ), x := (x 1 , . . . , x d ) ∈ R d , (2.16
)

C k (u) := d l=1 C (l) k l (u l ), u := (u 1 , . . . , u d ) ∈ [0, 1] d .
(2.17)

The transformation function F ν : R d → [0, 1] d and its inverse are defined by

F ν (x) := (F ν 1 (x 1 ), . . . , F ν 1 (x d )), F -1 ν (u) := (F -1 ν 1 (u 1 ), . . . , F -1 ν d (u d )).
Therefore, from the above and from (2.16)-(2.17), we readily check the relations

(C k • F ν )(x) := d l=1 (C (l) k l • F ν l )(x l ) = φ k (x), (2.18) (φ k • F -1 ν )(u) := d l=1 (φ (l) k l • F -1 ν l )(u l ) = C k (u). (2.19)
To analyze orthogonality properties, we repeatedly use the following lemma, which proof is quite straightforward and left to the reader (use the same change-ofvariable arguments as in the proof of Lemma 1).

Lemma 2 Assume ν satisfies (A ν ). Then, for all functions v 1 , v 2 ∈ L 2 ν (R d ), we have the following identities

(v 1 , v 2 ) L 2 ν (R d ) = (v 1 •F -1 ν , v 2 •F -1 ν ) L 2 U ([0,1] d ) , v 1 L 2 ν (R d ) = v 1 • F -1 ν L 2 U ([0,1] d )
.

As a consequence of Lemma 2 and Equality (2.19), since the tensored cosine functions (C k (.) : k ∈ N d ) are orthonormal in [0, 1] d with respect to the Lebesgue measure, it implies that (φ k (.) : k ∈ N d ) are orthonormal in R d with respect to the probability measure ν(dx). Moreover, the cosine functions (C k (.) :

k ∈ N d ) form a complete family in L 2 U ([0, 1] d ):
this is a standard result in Fourier analysis (observe that here, the sine functions are not needed to get a complete basis since the scalar product (., .)

L 2 U ([0,1] d ) is defined on [0, 1] d and not [-1, 1] d ).
We summarize our findings in a proposition ready to be used.

Proposition 2 Assume ν satisfies (A ν ). The family of functions (φ k (.) : k ∈ N d ) forms a complete orthonomal basis of R d with respect to the measure ν.

Global Quasi-Regression algorithm and convergence results

We now describe our main algorithm which is able to compute accurately the function y (q) i defined in (2.9). For this, define

S (q) i (x i i:N ) := g(x i N ) (1 + |x i i | 2 ) q/2 + N -1 j=i f j x i j , y (q) j+1 (x i j+1 )(1 + |x i j+1 | 2 ) q/2 (1 + |x i i | 2 ) q/2
∆,

(2.20)

x i i:N := (x i i , . . . , x i N ) ∈ (R d ) N -i+1 ,
which allows us to rewrite the exact solution as

y (q) i (x) = E S (q) i X i i:N |X i i = x , x ∈ R d .
Note that it is easy to justify that S (q)

i (X i i:N ) is square-integrable under (A f ,g ), (A ν ) and H ρ (q) (remind of the bound (2.4) and the integrability condition in (A ν )). Therefore, the exact solution y

(q) i is in L 2 ν (R d )
, and consequently (Proposition 2) it can be decomposed in its φ series:

y (q) i (x) = k∈N d α (q) i,k φ k (x) (2.21)
for some coefficients (α

(q) i,k : k ∈ N d ).
In practice, the coefficients α

(q) i,k are computed only in a finite set of multi-indices Γ ⊂ N d . Different choices of Γ are discussed later in Section 2.4.2.
The orthonormality property of the φ k 's (Proposition 2) combined with (2.21) yields α

(q) i,k = (y (q) i , φ k ) L 2 ν (R d ) . (2.22)
Thus, using the tower property of conditional expectation, we obtain, ∀k ∈ N d ,

α (q) i,k = E y (q) i (X i i )φ k (X i i ) = E E S (q) i X i i:N |X i i φ k (X i i ) = E S (q) i X i i:N φ k (X i i ) . (2.23)
The above expectations in (2.23) can be computed using Monte-Carlo simulations and evaluations of S (q) i along different sampled paths x i i:N . Since S (q) i depends on the unknown solution, we have to replace it by approximations, using previous steps in the dynamic programming equation (2.9). It gives rise to the following Algorithm 1, which we will call Global Quasi-Regression Multistep-forward Dynamical Programming (GQRMDP) algorithm. To maintain good integrability properties along the iteration, we use the truncation function T L (1.10) with the upper bound (2.4).

Algorithm 1 (GQRMDP algorithm)

Initialization. Set ȳ(q,M)

N (x N ) := g(x N ) (1+|x N | 2 ) q/2 . Backward iteration for i = N -1 to i = 0, ȳ(q,M) i (•) := k∈Γ ᾱ(q,M) i,k φ k (•), (2.24)
where for all k ∈ Γ ,

ᾱ(q,M) i,k := 1 M M m=1 S (q,M ) i (X i,m i:N )φ k (X i,m i ), (2.25)
and

S (q,M ) i (x i i:N ) := g(x i N ) (1 + |x i i | 2 ) q/2
(2.26)

+ N -1 j=i f j x i j , T L ȳ(q,M) j+1 (x i j+1 )(1 + |x i j+1 | 2 ) q/2 (1 + |x i i | 2 ) q/2
∆.

In the algorithm, in order to compute ȳ(q,M) i (•) we have considered M independent copies of

X i i Euler scheme ----------→ X i i+1 Euler scheme ----------→ . . . Euler scheme ----------→ X i N ,
which form a cloud of simulations C i = X i,m i:N : m = 1, . . . , M . We assume that the clouds of simulations (C i : 0 ≤ i < N ) are independently generated.

Our main result is the following error propagation theorem, which proof is postponed to Section 3.2.

Theorem 2 Assume (A f ,g ), (A b,σ ), (A ν ) and H ρ (q). Let k ∈ N d be a multiindex and let Γ ⊂ N d be a finite set of multi-indices. Let (φ k (.) : k ∈ N d ) be the basis functions defined in (2.18), as an orthonormal basis in L 2 ν (R d ), and let L Γ be the related Christoffel number [START_REF] Nevai | Géza Freud, orthogonal polynomials and Christoffel functions. a case study[END_REF] restricted to Γ :

L Γ = k∈Γ φ k (•) 2 ∞ ≤ 2 d #Γ. (2.27)
For each i ∈ {0, . . . , N -1}, define the following local error term

E i := k∈(N d -Γ ) (α (q) i,k ) 2 + L Γ M E S (q) i (X i i:N ) 2 .
(

2.28)

There is a finite positive constant C (depending on model parameters but not on N, M, Γ ) such that

E ȳ(q,M) i (•) -y (q) i (•) 2 ν ≤ C   E i + N -1 j=i E j+1 ∆   .
Observe that provided that L Γ M remains bounded, the constant C can be taken as uniform w.r.t. the algorithm parameters. This is the situation in practice since for getting accurate approximations, we need L Γ M to be small.

The convergence rate of the term

k∈(N d -Γ ) (α (q) i,k ) 2 is discussed in Section 2.4
considering ν as the Student's t-distribution (Example 1). This analysis will be undertaken assuming that y i (x) has a polynomial growth in x at infinity (coherently with the a priori estimate of Proposition 1) as well as its derivatives.

Truncation error for the Student-cosine basis

From Equations (2.21)-(2.22) and Lemma 2 we recall that y

(q) i = k∈N d α (q) i,k φ k , where α (q)
i,k = (y

(q) i , φ k ) L 2 ν (R d ) = (y (q) i • F -1 ν , C k ) L 2 U ([0,1] d ) . (2.29)
Define the truncation of the ν-series of y (q) i by

P Γ y (q) i := k∈Γ α (q) i,k φ k .
Therefore, from orthonormal properties of the basis functions, we have

y (q) i -P Γ y (q) i 2 L 2 ν (R d ) = y (q) i • F -1 ν - k∈Γ α (q) i,k C k 2 L 2 U ([0,1] d ) = k∈(N d -Γ ) |α (q) i,k | 2 .
In other words, α

i,k are the Fourier coefficients of the projection of

h (q) i := y (q) i • F -1 ν ∈ L 2 U ([0, 1] d ) (2.30)
over cosine basis functions. In the following section we will take advantage of this property to derive a tight bound for the truncation error k∈(N d -Γ ) |α

(q)
i,k | 2 , using assumptions on y (q) i and its derivatives only (and not on y

(q) i • F -1 ν ).
The specific choice of the Student's t-distribution measure ν of Example 1 plays an important role to get tractable conditions. Now our aim is twofold. Firstly we prove that if y i is r times (r ≥ 1) continuously differentiable on R d , with (for some p ∈ N),

lim sup |x|→+∞ |∂ n x l y i (x)| |x| p-n < +∞, for n = 0, . . . , r, l = 1, . . . , d,
and if the parameters are such that

p + r - 1 2 µ < q,
then the truncation error is controlled as

k∈N d :∃k l >K l |α (q) i,k | 2 ≤ d l=1 C (K l + 1) 2r
for some constant C > 0. In the following, we first establish the above result in dimension 1 (Proposition 3) and then in dimension d (Proposition 4) .

Our second goal is improving this result when mixed derivatives of y i are also supposed smooth. Thus, we prove that if y i is Λ-smooth, where Λ(r) = {0, . . . , r} d , with (for some p ∈ N),

lim sup |x|→+∞ |∂ n x y i (x)| |x| p-n < +∞, for n ∈ Λ,
and if the parameters are such that

p + d r - 1 2 µ < q,
then the truncation error satisfies the tail inequality

k∈ ΓH (DEG) |α (q) i,k | 2 ≤ C (DEG + 1) 2r ,
for some constant C > 0 (Proposition 5).

One-dimensional case

In this section, we let ν = ν 1 be the Student density (2.11). Lemma 3 Assume that y i is r times (r ∈ N) continuously differentiable on R, with the existence of a polynomial growth degree p ∈ N such that lim sup x→±∞ |d n

x y i (x)| |x| p-n < +∞, for n = 0, . . . , r.

(2.31)

Assume the relation p + rµ < q (2.32) between r, p, the Student parameter µ > 0 and the damping exponent q ≥ 0. Then h (q) i is r times continuously differentiable on (0, 1), and its first r derivatives vanish at u = 0 and u = 1:

lim u→0 + ,1 -d n u h (q) i (u) = 0, n = 0, . . . , r.
Proof The regularities of y (q) i and h

(q) i are straightforward in view of (2.8)-(2.30) and Lemma 15. The estimate lim u→0 + ,1 -h (q) i (u) = 0 is clearly true since it is equivalent to 0 = lim x→±∞ y (q) i (x), which holds owing to the condition q > p.

From now on, we may consider the case n ≥ 1 and r ≥ 1. We start showing that for 1 ≤ n ≤ r,

d n x y (q) i (x) = O(x p-n-q ) as x → ±∞. (2.33)
In fact, using Leibniz rule,

d n x y (q) i (x) = n m=0 n m d n-m x y i (x) d m x 1 (1 + x 2 ) q/2 . Besides, d n-m x y i (x) = O(x p-n+m ), d m x 1 (1 + x 2 ) q/2 = O(x -q-m ), as x → ±∞,
which leads directly to the desired result (2.33).

We now estimate d n u h (q) i (u). In order to compute the n-th derivative of h (q) i we take advantage of the following Faà di Bruno's formula [29, pp. 224-226],

d n u h (q) i (u) = d n u y (q) i (F -1 ν 1 (u)) = n m=1 1 m! d m x y (q) i (x) x=F -1 ν 1 (u) j∈J m,n n! j 1 !j 2 ! • • • j m ! m i=1 d j i u F -1 ν 1 (u), J m,n = {j = (j 1 , . . . , j m ) ∈ N m + : j 1 + . . . + j m = n}. (2.

34)

As u → 0 + , using also (A.7),

d n u h (q) i (u) = O max 1≤m≤n, j∈J m,n F -1 ν 1 (u) p-m-q d j 1 u F -1 ν 1 (u) • • • d j m u F -1 ν 1 (u) = O max 1≤m≤n, j∈J m,n u -p-m-q µ u - 1+j 1 µ µ • • • u -1+jm µ µ = O u -p+nµ-q µ . (2.35)
The above bound goes to 0 as u → 0 + since p + nµ -q ≤ p + rµ -q < 0. Analogously, as u → 1 -, we have

d n u h (q) i (u) = O (1 -u) -p+nµ-q µ (2.36)
which converges to 0 as u → 1 -since p + nµ < q.

Lemma 4 Assume the hypotheses of Lemma 3 and replace the condition (2.32) by the weaker condition

p + r - 1 2 µ < q.
(2.37)

Then the r first derivatives of h (q) i are square integrable on (0, 1):

d n u h (q) i 2 L 2 U ([0,1])
< +∞, n = 0, . . . , r.

Proof The integrability of d n u h

(q) i on any compact set of the form [ε, 1 -ε] (ε > 0) is obvious owing to the continuity of the function to integrate. As a consequence one has only to consider what happens around 0 and 1. As we have just seen in the proof of Lemma 3 (see (2.35)

-(2.36)), d n u h (q) i (u) = O u→0 + u -p+nµ-q µ , d n u h (q) i (u) = O u→1 -(1 -u) -p+nµ-q µ ;
therefore it is square-integrable around u = 0 and u = 1 provided that -p+nµ-q µ > -1 2 . Since n ≤ r the condition (2.37) is sufficient.

Proposition 3 Assume that y i is r times (r ≥ 1) continuously differentiable on R, with the existence of a polynomial growth degree p ∈ N such that (2.31) holds. Assume the relation (2.37) holds. Then, the coefficients (α

(q) i,k , k ∈ N) defined in (2.29) satisfy the tail inequality (for any K ≥ 0) ∞ k=K+1 |α (q) i,k | 2 ≤ C (r) (2.38) (K + 1) 2r , where C (r) (2.38) := 1 π 2r d r u h (q) i 2 L 2 U ([0,1]) < +∞.
(2.38)

Proof Under (2.37), we can invoke Lemma 3 (with r -1 derivatives) and Lemma 4 (with r derivatives), to get

lim u→0 + ,1 -d n u h (q) i (u) = 0, for n = 0, . . . , r -1, d n u h (q) i 2 L 2 U ([0,1])
< +∞, for n = 0, . . . , r.

(2.39)

The coefficients (α

(q)
i,k , k ∈ N) coincide with the Fourier coefficients of the decomposition of h (q) i with respect to the orthonormal basis functions

{C k , k ∈ N} in L 2 U ([0, 1]) (see the relations (2.29)-(2.30)). The C k -cosine series of h (q) i writes h (q) i (u) = ∞ k=0 α (q) i,k C k (u) = α (q) i,0 + ∞ k=1 α (q) i,k C k (u),
where

α (q) i,0 = 1 0 h (q)
i (u)du and α

(q) i,k = √ 2 1 0 h (q)
i (u) cos(kπu)du for k > 0.

Applying integration by parts (h (q) i is continuously differentiable) and using that lim u→0 + ,1 -h (q) i (u) = 0 we get (for k = 0)

α (q) i,k = Re √ 2 lim ε→0 1-ε ε h (q) i (u)e ikπu du = Re √ 2 ikπ h (q) i (u)e ikπu 1 - 0 + - √ 2 ikπ 1 0 d u h (q) i (u)e ikπu du = Re √ 2 -ikπ 1 0 d u h (q) i (u)e ikπu du .
(2.40)

Thanks to (2.62), we can repeat r -1 times the integration by parts as above, and we obtain i,k : k ≥ 1} are the projection coefficients of ±d r u h

α (q) i,k = Re √ 2 (-ikπ) r 1 0 d r u h (q) i (u)e ikπu du = 1 (kπ) r 1 0 d r u h (q) i (u)Re √ 2 i r e ikπu du. ( 2
(q) i on these orthonormal functions. As such, the Bessel inequality gives

∞ k=1 (kπ) 2r |α (q) i,k | 2 ≤ d r u h (q) i 2 L 2 U ([0,1]) = π 2r C (r) (2.38) .
It readily follows

∞ k=K+1 |α (q) i,k | 2 ≤ 1 π 2r ∞ k=K+1 (kπ) 2r (K + 1) 2r |α (q) i,k | 2 ≤ C (r) (2.38) (K + 1) 2r .

Multi-indices sets

In this work only arbitrary downward closed multi-indices sets are considered. As already stated in the introduction, the finite multi-indices set

Γ ⊂ N d is downward closed if (k ∈ Γ and j ≤ k) ⇒ j ∈ Γ,
where j ≤ k means that j l ≤ k l for all l = 1, . . . , d. The most well-known of such multi-indices sets is the full or tensor product multi-indices set. For any K = (K 1 , . . . , K d ) where K l ∈ N for all l = 1, . . . , d, the full anisotropic multiindices set and its complement are defined by

Γ F (K 1 , . . . , K d ) := {k ∈ N d : k l ≤ K l , ∀l = 1 . . . d}, ΓF (K 1 , . . . , K d ) := {k ∈ N d : k / ∈ Γ F (K 1 , . . . , K d )}.
In high dimensions the full multi-indices set suffers from the curse of dimensionality. In order to effectively compute the projection of y (q) i over a multidimensional set of of basis functions we propose to consider the following two sparse sets of downward closed multi-indices sets, where DEG denotes a nonnegative integer:

-Total degree set defined by Secondly we divide Γ d 1 (DEG) in the following disjoint sets,

Γ T (DEG) := k ∈ N d : d l=1 k l ≤ DEG ,
Γ d 1 (DEG) = DEG g=2 Υ d (g) ∪ {1}.
For any integer g ∈ {2, . . . , DEG} having the prime factorization

g = i p v i i ,
#Υ d (g) is equal to the number of possible ways to write g as a product of exactly d of its factors,

#Υ d (g) = i v i + d -1 d -1 . Therefore, #Γ d 1 (DEG) = DEG g=2 #Υ d (g) + 1.
On the other hand, let us compute #Γ 2 (DEG). Firstly we define Υ c (DEG) as the subset of Γ 2 (DEG) containing the multi-indices of Γ 2 (DEG) with exactly c of the

d indexes equal to zero, 0 < c < d. Note that #Υ c (DEG) = d c #Γ d-c 1 (DEG). Therefore, since Γ 2 (DEG) = d-1 c=1 Υ c (DEG) ∪ {0}, #Γ 2 (DEG) = d-1 c=1 d c #Γ d-c 1 (DEG) + 1.
Finally, we conclude using that #Γ

H (DEG) = #Γ d 1 (DEG) + #Γ 2 (DEG).
Note that ∀g ∈ N, g ≤ (DEG -d + 1) + , Γ H (g) ⊂ Γ T (DEG). See Figure 1 and Table 1. Table 1 Cardinalities of Γ T and Γ H . between r, p, the Student parameter µ > 0 and the damping exponent q ≥ 0. Then h (q) i is r times continuously differentiable on (0, 1) d , and its first r coordinate-wise derivatives vanish at ∂[0, 1] d :

d = 3 d = 4 DEG #Γ T (DEG) #Γ H (DEG -2) DEG #Γ T (DEG) #Γ H (DEG -
lim u→∂[0,1] d ∂ n u l h (q) i (u) = 0, n = 0, . . . , r, l = 1, . . . , d.
Proof As in dimension 1, the regularities of y (q) i and h

(q) i are straightforward to obtain. The statement lim u→∂[0,1] d h (q) i (u) = 0 is clear since the limit coincides with that of y (q) i (x) as |x| → +∞, which equals to 0 in view of q > p. Therefore, from now on, we may consider n ≥ 1 and let l = 1, . . . , d be fixed. In the following, to deal with comparison of functions as |x| tends to +∞, we use the same notation C for any generic constants, that do not depend on x, but that may depend on the order n, r, . . . of derivatives under consideration. We will keep the same C from line to line, in order the alleviate the computations, although its value changes.

We start by showing that for n = 1, . . . , r,

|∂ n x l y (q) i (x)| ≤ C(1 + |x| 2 ) p/2-n/2-q/2 . (2.44)
By applying the Leibniz rule to the product (2.8) (as in dimension 1, since we differentiate n times with respect to the same variable x l ) we get

∂ n x l y (q) i (x) = n m=0 n m ∂ n-m x l y i (x) ∂ m x l 1 (1 + |x| 2 ) q/2 .
(2.45)

Our standing assumptions write

|∂ n-m x l y i (x)| ≤ C(1 + |x|) p-(n-m) . ( 2

.46)

From (A.2) we know that

∂ m x l 1 (1 + |x| 2 ) q/2 ≤ C(1 + |x| 2 ) -q/2-m/2 .
By plugging the previous upper bound and (2.46) into the Leibniz formula (2.45), we get

∂ n x l y (q) i (x) ≤ C max 0≤m≤n (1 + |x|) p-(n-m) (1 + |x| 2 ) -q/2-m/2 ≤ C(1 + |x| 2 ) p/2-q/2-n/2 ,
we are done with (2.44).

In order to compute the n-th derivative of h

(q) i (u) = y (q) i (F -1 ν 1 (u 1 ), . . . , F -1 ν d (u d ))
we again take advantage of the Faà di Bruno's formula (as in (2.34))

∂ n u l h (q) i (u) = n m=1 1 m! ∂ m x l y (q) i (x) x=(F -1 ν 1 (u 1 ),...,F -1 ν d (u d )) j∈J m,n n! j 1 !j 2 ! • • • j m ! m i=1 d j i u l F -1 ν l (u l ).
Combining the above with Lemma 11 and the estimates (A.7), we obtain

∂ n u l h (q) i (u) ≤ C max 1≤m≤n, j∈J m,n 1 + d i=1 |F -1 ν i (u i )| 2 p/2-q/2-m/2 m i=1 d j i u l F -1 ν l (u l ) ≤ C max 1≤m≤n, j∈J m,n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m 1 u 1/µ l (1 -u l ) 1/µ m+µ m i=1 j i ≤ C max 1≤m≤n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ m+nµ ≤ C 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q+nµ . (2.47)
Since p + nµ < q, we get that lim u→∂

[0,1] d ∂ n u l h (q)
i (u) = 0. Lemma 6 Assume the hypotheses of Lemma 5 and replace the condition (2.43) by the weaker condition p + r -1 2 µ < q.

(2.48)

Then the r first coordinate-wise derivatives of h (q) i are square integrable on (0, 1) d :

∂ n u l h (q) i (u) 2 L 2 U ([0,1] d )
< +∞, n = 0, . . . , r, l = 1, . . . , d.

Proof In view of the bound (2.47), it is sufficient to check the integrability on (0, 1)

d of l(u) := 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ 2(p-q+rµ)
.

Observe that the function is symmetric with respect to each hyperplan {u ∈ (0, 1) d : u l = 1 2 }, l = 1, . . . , d. Therefore, only the integrability on (0, 1/2) d has to be established, which is equivalent to prove

(0,1/2) d 1 + d i=1 1 u 1/µ i 2(p-q+rµ) du < +∞
Use the change of variables z i = u -1/2µ i , then the above integral writes (up to a non-zero factor)

I := (2 1/(2µ) ,+∞) d 1 + d i=1 z 2 i µ-ε d i=1 z -2µ-1 i dz,
where ε > 0 is given by 2(p -q + rµ) = µ -ε. To show that I < +∞, we may decrease ε ↓ 0 as much as needed to get µ -ε > 0. Once done, using the inequality

(x 1 + x 2 ) α ≤ 2 (α-1) + (x α 1 + x α 2 ),
available for any x 1 ≥ 0, x 2 ≥ 0, α > 0, it yields (by direct induction) the existence of a constant C α,d such that

1 + d i=1 x i α ≤ C α,d 1 + d i=1
x α i for any x 1 ≥ 0, . . . , x d ≥ 0. Consequently,

I ≤ C µ-ε,d (2 1/(2µ) ,+∞) d 1 + d i=1 z 2(µ-ε) i d i=1 z -2µ-1 i dz < +∞.
We are done.

Proposition 4

Assume that y i is r times (r ≥ 1) coordinate-wise continuously differentiable on R d with the existence of a polynomial growth degree p ∈ N such that (2.42) holds. Assume the relation (2.48) holds. Then, the coefficients (α

(q) i,k , k ∈ N d ) defined in (2.29) satisfy the tail inequality k∈ ΓF (K 1 ,...,K d ) |α (q) i,k | 2 ≤ d l=1 C (l,r) (2.49) (K l + 1) 2r (2.49)
where

C (l,r) (2.49) := 1 π 2r ∂ r u l h (q) i 2 L 2 U ([0,1] d ) < +∞.
Proof By definition of Γ F (K 1 , . . . , K d ), we have the union bound

{k ∈ ΓF (K 1 , . . . , K d )} = d l=1 {k ∈ N d : k l > K l },
therefore the advertised inequality (2.49) directly follows from

k∈N d :k l >K l |α (q) i,k | 2 ≤ C (l,r) (2.49) (K l + 1) 2r , l = 1, . . . , d. (2.50)
We only prove the case l = 1, the arguments apply similarly for other l.

Under (2.48), we can invoke Lemma 5 (with r -1 derivatives) and Lemma 6 (with r derivatives), to get lim

u→∂[0,1] d ∂ n u 1 h (q)
i (u) = 0, for n = 0, . . . , r -1,

∂ n u 1 h (q) i 2 L 2 U ([0,1] d )
< +∞, for n = 0, . . . , r.

(2.51)

In view of (2.29)-(2.30), the Fourier coefficients of the decomposition of h (q) i with respect to the orthonormal basis functions

{C k , k ∈ N d } in L 2 U ([0, 1] d ) are given by (α (q) i,k , k ∈ N d )n and the C k -cosine series of h (q) i writes h (q) i (u) = k∈N d α (q) i,k C k (u), where α (q) i,k = β k [0,1] d h (q) i (u) cos(k 1 πu 1 ) cos(k 2 πu 2 ) • • • cos(k d πu d )du = β k Re [0,1] d h (q) i (u)e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d )du , β k := d l=1 (1 k l =0 + 1 k l >0 √ 2).
Applying integration by parts in u 1 (h (q) i is uniformly bounded on (0, 1) d , continuously differentiable and the derivative with respect to u 1 is integrable) and using that lim

u 1 →0 + ,1 -h (q) i (u) = 0, we get for k 1 > 0 α (q) i,k = Re β k ik 1 π [0,1] d-1 h (q) i (u)e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d )du 2 . . . du d u 1 =1 - u 1 =0 + - β k ik 1 π [0,1] d ∂ u 1 h (q) i (u)e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d )du = Re β k -ik 1 π [0,1] d ∂ u 1 h (q) i (u)e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d )du .
Thanks to (2.51), we can repeat r -1 times the integration by parts as above, and we obtain

α (q) i,k = Re β k (-ik 1 π) r [0,1] d ∂ r u 1 h (q) i (u)e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d )du = 1 (k 1 π) r [0,1] d ∂ r u 1 h (q) i (u)Re β k i r e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d ) du.
(2.52)

Observe that the sets

{u → β k cos(k 1 πu 1 ) cos(k 2 πu 2 ) • • • cos(k d πu d ), k ∈ N d }, {u → β k sin(k 1 πu 1 ) cos(k 2 πu 2 ) • • • cos(k d πu d ), k ∈ N d }, are made of orthonormal functions in L 2 U ([0, 1] d ). Consequently, Re β k i r e ik 1 πu 1 cos(k 2 πu 2 ) • • • cos(k d πu d ) , k ∈ N d
is (for any value of r) a set of orthonormal functions and the equality (2.52) means that {(k 1 π) r α (q) i,k : k ∈ N d } are the projection coefficients of ±∂ r u 1 h (q) i on these orthonormal functions. As such, the Bessel inequality gives

k∈N d (k 1 π) 2r |α (q) i,k | 2 ≤ ∂ r u 1 h (q) i 2 L 2 U ([0,1] d ) = π 2r C (1,r) (2.49) .
It readily follows

k∈N d :k 1 >K 1 |α (q) i,k | 2 ≤ 1 π 2r k∈N d :k 1 >K 1 (k 1 π) 2r (K 1 + 1) 2r |α (q) i,k | 2 ≤ C (1,r) (2.49) (K 1 + 1) 2r .
The proof of (2.50) is now finished and Proposition 4 is proved. Now we consider that not only y i and its coordinate-wise derivatives are smooth but also its mixed derivatives.

Lemma 7 Let Λ ⊂ N d be an arbitrary downward closed multi-index set and set r = max r∈Λ (r). Assume that y i is Λ-smooth with the existence of a polynomial growth degree p ∈ N such that

lim sup |x|→+∞ |∂ n x y i (x)| |x| p-n < +∞, for n ∈ Λ. (2.53)
Assume the relation p + rµ < q (2.54) between r, p, the Student parameter µ > 0 and the damping exponent q ≥ 0. Then h (q) i is Λ-smooth on (0, 1) d and its derivatives vanish at ∂[0, 1] d :

lim u→∂[0,1] d ∂ n u h (q) i (u) = 0, n ∈ Λ.
Proof We start by showing that for n > 0, n ∈ Λ,

|∂ n x y (q) i (x)| ≤ C(1 + |x| 2 ) p/2-n/2-q/2 . (2.55)
By applying the multivariate Leibniz rule to the product (2.8) we get

∂ n x y (q) i (x) = {m∈N d :m≤n} C n,m ∂ n-m x y i (x) ∂ m x 1 (1 + |x| 2 ) q/2 .
(2.56)

Our standing assumptions write

|∂ n-m x y i (x)| ≤ C(1 + |x|) p-(n-m) ≤C(1 + |x| 2 ) p/2-n/2+m/2 .
(2.57) From (A.4) we know that

∂ m x 1 (1 + |x| 2 ) q/2 ≤ C(1 + |x| 2 ) -q/2-m/2 .
By plugging the previous upper bound and (2.57) into the Leibniz formula (2.56), we get

∂ n x y (q) i (x) ≤ C max 0≤m≤n (1 + |x| 2 ) p/2-n/2+m/2 (1 + |x| 2 ) -q/2-m/2 ≤ C(1 + |x| 2 ) p/2-q/2-n/2 ,
we are done with (2.55).

Let [ * , . . . , * ] denotes a list. In order to compute the n

(n > 0, n ∈ Λ ) derivative of h (q) i (u) = y (q) i (F -1 ν 1 (u 1 ), . . . , F -1 ν d (u d ))
we take advantage of the multivariate Faà di Bruno's formula derived from [7, Theorem 2.1]:

∂ n u h (q) i (u) = m∈M n ∂ m x y (q) i (x) x=(F -1 ν 1 (u 1 ),...,F -1 ν d (u d )) J∈A m,n C m,J d i=1 m i k=1 d j ik u i F -1 ν i (u i ), M n = {m = (m 1 , . . . , m d ) ∈ N d + : ∀i = 1, . . . , d, 0 ≤ m i ≤ n i , m = 0}, A m,n = {J = [j 1 ∈ N m 1 + , . . . , j d ∈ N m d + ] : j i = n i ∀i = 1, . . . , d}.
Combining the above with Lemma 11 and the estimates (A.7), we obtain

∂ n u h (q) i (u) ≤ C max m∈M n , J∈A m,n 1 + d i=1 |F -1 ν i (u i )| 2 p/2-q/2-m/2 d i=1 m i k=1 d j ik u i F -1 ν i (u i ) ≤ C max m∈M n , J∈A m,n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m d i=1 1 u 1/µ i (1 -u i ) 1/µ m i +µ m i l=1 j il ≤ C max m =0,m≤n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m d i=1 1 u 1/µ i (1 -u i ) 1/µ m i +µn i ≤ C max m =0,m≤n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m d i=1   1 + d j=1 1 u 1/µ j (1 -u j ) 1/µ   m i +µn i ≤ C max m =0,m≤n 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q-m+m+µn (2.58) ≤ C 1 + d i=1 1 u 1/µ i (1 -u i ) 1/µ p-q+nµ . (2.59) Since p + nµ < q, we get that lim u→∂[0,1] d ∂ n u h (q) i (u) = 0.
Lemma 8 Assume the hypotheses of Lemma 7 and replace the condition (2.54) by the weaker condition

p + r - 1 2 µ < q.
(2.60)

Then the Λ derivatives of h (q) i are square integrable on (0, 1) d :

∂ n u h (q) i (u) 2 L 2 U ([0,1] d ) < +∞, n ∈ Λ.
Proof The proof follows readily from the proof of Lemma 6 by replacing r by r = max r∈Λ (r).

Proposition 5 Let Λ(r) = {0, . . . , r} d . Assume that y i is Λ-smooth with the existence of a polynomial growth degree p ∈ N such that (2.53) holds. Assume the relation (2.60) holds, i.e. p + (d r -1/2)µ < q.

Then, the coefficients (α

(q) i,k , k ∈ N d ) defined in (2.29) satisfy the tail inequality k∈ ΓH (DEG) |α (q) i,k | 2 ≤ C (2.61) (DEG + 1) 2r , (2.61)
where

C (2.61) := 1 π 2dr ∂ (r,...,r) u h (q) i 2 L 2 U ([0,1] d ) < +∞.
Proof Under (2.60), we can invoke Lemma 7 (with #Λ(r) -1 derivatives) (condition p + (dr -1)µ < q) and Lemma 8 (with #Λ(r) derivatives) (condition p + (dr -1/2)µ < q), to get lim

u→∂[0,1] d ∂ n u h (q) i (u) = 0, n ∈ Λ(r), n = (r, . . . , r), ∂ n u h (q) i (u) 2 L 2 U ([0,1] d ) < +∞, n ∈ Λ(r).
(2.62)

By doing integration by parts on the Fourier coefficients of h (q) i as we have already done in Proposition 4 we readily obtain

α (q) i,k = [0,1] d h (q) i (u)β k cos(k 1 πu 1 ) • • • cos(k d πu d )du, = 1 π dr (k 1 • • • k d ) r [0,1] d ∂ (r,...,r) u h (q) i (u)β k cos ,r (k 1 πu 1 ) • • • cos ,r (k d πu d )du, (2.63)
where cos ,r (.) is the r-th antiderivative of cos(•) such that

cos ,4m (•) = cos(•), cos ,4m+1 (•) = sin(•), cos ,4m+2 (•) = -cos(•), cos ,4m+3 (•) = -sin(•), m ∈ N.

Equation (2.63) tells us that h

(q) i is a function in a Korobov space [30, Definition 6.5.6, p. 225]. Observe that the sets

{u → β k cos(k 1 πu 1 ) cos(k 2 πu 2 ) • • • cos(k d πu d ), k ∈ N d }, {u → β k sin(k 1 πu 1 ) cos(k 2 πu 2 ) • • • cos(k d πu d ), k ∈ N d }, . . . {u → β k sin(k 1 πu 1 ) sin(k 2 πu 2 ) • • • sin(k d πu d ), k ∈ N d }, are made of orthonormal functions in L 2 U ([0, 1] d ). Consequently, β k cos ,r (k 1 πu 1 ) • • • cos ,r (k d πu d ), k ∈ N d
is (for any value of r) a set of orthonormal functions and the equality (2.63) means that {(k

1 π) r • • • (k d π) r α (q) i,k : k ∈ N d } are the projection coefficients of ∂ (r,...,r) u h (q) i
on these orthonormal functions. As such, the Bessel inequality gives

k∈N d π 2dr (k 1 • • • k d ) 2r |α (q) i,k | 2 ≤ ∂ (r,...,r) u h (q) i 2 L 2 U ([0,1] d ) = π 2dr C (2.61) .
It readily follows

k∈N d , k / ∈Γ H (DEG) |α (q) i,k | 2 ≤ 1 π 2dr k∈N d , k / ∈Γ H (DEG) π 2dr (k 1 • • • k d ) 2r (DEG + 1) 2r |α (q) i,k | 2 ≤ C (2.61) (DEG + 1) 2r .
Equation (2.63) means that the most significant coefficients α (q) i,k are the ones for which the product k 1 • • • k d is as small as possible. Under this setting we prefer to consider the hyperbolic cross index set (Proposition 5) rather than the total degree set (Proposition 4). In fact, although the hyperbolic truncation error with Γ H ((DEG -d + 1) + ) will be greater than the total truncation error with Γ T (DEG), if h (q) i belongs to a Korobov space, the hyperbolic truncation error will be near to the total truncation error (here we are assuming that DEG is big enough to take into account all the significant Fourier coefficients). Besides, having in mind that #Γ H (DEG -d + 1) #Γ T (DEG) (specially in high dimensions), the statistical error, driven by 2 d #Γ M , will be smaller for the hyperbolic multi-indices set, considering the same number of simulations M . Last but not least, the algorithm executed with the hyperbolic multi-indices set will be less time consuming than with the total degree set. This behavior is shown in the numerical experiments of Section 5.

3 Proofs related to error analysis

Proof of Theorem 1

First observe that we can restrict in the proof of (2.13) to the case of bounded measurable functions γ: indeed, for γ ∈ L 2 ν (R d ) (possibly unbounded), apply the inequality (2.13) to the bounded function x → T L (γ(x)) where L > 0 to get (for any 0

≤ i ≤ j ≤ N ) E |T L (γ(X i j ))| 2 (1 + |X i j | 2 ) q (1 + |X i i | 2 ) q ≤ c ν E |T L (γ(X i i ))| 2 , (3.1)
and take the limit L → +∞. On the one hand, the dominated convergence theorem (because γ ∈ L 2 ν ) ensures that the r.h.s. of (3.1) converges to c ν E |γ(X i i )| 2 . On the other hand, by Fatou's lemma, the liminf of the l.h.s. of (3.1)

is bounded from below by E |γ(X i j )| 2 (1 + |X i j | 2 ) q (1 + |X i i | 2 ) q
, which all in all, gives the result (2.13) for the square integrable function γ. Now we establish (2.13) for bounded functions γ. We follow and adapt the arguments of [24, Proposition 3.1 and Proposition 3.3 item (b)] which paves the way for q = 0 (although the distribution ν is different). Under the current assumptions on b and σ, the Euler scheme X is a non-homogeneous Markov chain which transition kernel has a probability density w.r.t. the Lebesgue measure on R d . This density is denoted by p i,j (x, y) where i < j stands for the time indices and x (resp. y) stands for the initial (resp. final) space point. This probability density satisfies to a Gaussian upper bound of Aronson type (see [START_REF] Menozzi | On some non asymptotic bounds for the Euler scheme[END_REF]Theorem 2.1]): there is a finite constant C (3.2) > 0 (depending on d, q, T , b, σ but not on N ) such that

p i,j (x, y) ≤ C (3.2) exp -|y-x| 2 2C (3.2) (t j -t i ) (2π(t j -t i )) d/2 , (3.2)
for any 0 ≤ i < j ≤ N and any x, y ∈ R d . Assume that ν is such that, for any Λ ≥ 0, there is a finite constant C ν (Λ) > 0 (depending on µ and q) such that, for any λ ∈ [0, Λ] and any y ∈ R d , we have

R d ν(y + z √ λ) (1 + |y| 2 ) q (1 + |y + z √ λ| 2 ) q e -|z| 2 /2 (2π) d/2 dz ≤ C ν (Λ)ν(y). (3.3)
Then we claim that the inequality (2.13) holds, still considering bounded γ (which clearly ensures that both sides of (2.13) are finite). Indeed, combining (3.2) and (3.3), we get

E |γ(X i j )| 2 (1 + |X i j | 2 ) q (1 + |X i i | 2 ) q = R d R d ν(x)p i,j (x, y) (1 + |y| 2 ) q (1 + |x| 2 ) q γ 2 (y)dxdy (use (3.
2) and the change of variable (x, y) Therefore, it remains to establish (3.3); since

→ (y + z C (3.2) (t j -t i ), y)) ≤ C (3.2) R d R d ν(x) e - |y-x| 2 2C (3.2) (t j -t i ) (2π(t j -t i )) d/2 (1 + |y| 2 ) q (1 + |x| 2 ) q γ 2 (y)dxdy = C 1+d/2 (3.2) R d R d ν y + z C (3.2) (t j -t i ) × (1 + |y| 2 ) q (1 + |y + z C (3.2) (t j -t i )| 2 ) q e -|z| 2 /2 (2π) d/2 γ 2 (y)dzdy ≤ C 1+d/2 (3.2) C ν (C (3.2) T ) R d
(1 + z 2 ) 1/2 ≤ (1 + z) ≤ √ 2(1 + z 2 ) 1/2 , ∀z ≥ 0,
for getting (3.3) it is enough to prove (for any y = (y 1 , . . . , y d ))

R d d l=1 1 (1 + |y l + z l √ λ|) µ+1 (1 + |y| 2 ) q (1 + |y + z √ λ| 2 ) q e -|z| 2 /2 (2π) d/2 dz ≤ C ν (Λ) 2 (µ+1)d/2 d l=1 1 (1 + |y l |) µ+1 . (3.4) Set A λ := {z ∈ R d : |z| √ λ ≤ |y|/2}.
Denote by I 1 (resp. I 2 ) the above integral restricted to A λ (resp. A c λ ). On A λ , use that 1+|y +z

√ λ| 2 ≥ 1+ 1 2 |y| 2 -|z √ λ| 2 ≥ 1 4 + 1 4 |y| 2 to obtain I 1 ≤ 4 q A λ d l=1 1 (1 + |y l + z l √ λ|) µ+1 e -|z| 2 /2 (2π) d/2 dz ≤ 4 q d l=1 R 1 (1 + |y l + z l √ λ|) µ+1 e -|z l | 2 /2 √ 2π dz l .
Split each integral according to the set A λ,l := {z l ∈ R :

|z l | √ λ ≤ |y l |/2} and use the bound 1 + |y l + z l √ λ| ≥ 1 A c λ,l + 1 2 (1 + |y l |)1 A λ,l : it gives R 1 (1 + |y l + z l √ λ|) µ+1 e -|z l | 2 /2 √ 2π dz l ≤ A c λ,l e -|z l | 2 /2 √ 2π dz l + 2 µ+1 A c λ,l 1 (1 + |y l ) µ+1 e -|z l | 2 /2 √ 2π dz l ≤ 2N - |y l | 2 √ Λ + 2 µ+1 1 (1 + |y l ) µ+1 ≤ C µ,Λ (1 + |y l ) µ+1
where N (•) is the cumulative density function of the standard normal distribution, where C µ,Λ is a finite positive constant depending only on µ and Λ, and where at the last line, we have used that N (-x) is decreasing when x → +∞ faster than any polynomial. Therefore, we have proved

I 1 ≤ 4 q d l=1 C µ,Λ (1 + |y l ) µ+1 . (3.5)
Now we deal with I 2 , i.e. the part of the integral of (3.4) restricted to A c λ . On this set, |z| 2 /2 ≥ |z| 2 /4 + |y| 2 /(16Λ), therefore we obtain

I 2 ≤ (1 + |y| 2 ) q A c λ e -|z| 2 /2 (2π) d/2 dz ≤ (1 + |y| 2 ) q e -|y| 2 /(16Λ) A c λ e -|z| 2 /4 (2π) d/2 dz ≤ 2 d/2 (1 + |y| 2 ) q e -|y| 2 /(16Λ) ≤ C µ,q,d,Λ d l=1 1 (1 + |y l |) µ+1 (3.6)
where the last inequality follows from usual comparisons between exponential and polynomial growth. Gathering (3.5) and (3.6) gives the announced estimate (3.4) for some constant C ν (Λ). Theorem 1 is proved.

Proof of Theorem 2

In all this section, we assume the assumptions of Theorem 2 are in force, without further reference.

In addition to the coefficients ᾱ(q,M)

i,k calculated with the responses computed using the erroneous functions ȳ(q,M)

j i<j<N
(see the definitions (2.24) and (2.25)), we define the following coefficients ᾱ(q) i,k calculated with the responses computed with the exact functions y

(q) j i<j<N , ᾱ(q) i,k := 1 M M m=1 S (q) i X i,m i:N φ k (X i,m i ), (3.7) and thus ȳ(q)
i (•) := k∈Γ ᾱ(q) i,k φ k (•). (3.8)

Preliminary results

The next result states a decomposition of the local error as a superposition of an approximation error (truncation over the functions basis) and a statistical error.

Lemma 9 (Bias/variance decomposition) For each i ∈ {0, . . . , N -1} we have

E ȳ(q) i (•) -y (q) i (•) 2 ν ≤ k∈(N d -Γ ) (α (q) i,k ) 2 + L Γ M E S (q) i X i i:N 2 =: E i
where L Γ is the Christoffel number of the functions basis defined in (2.27) and the definition of E i is stated in (2.28).

Proof By definition of ȳ(q) i and y

(q) i , E ȳ(q) i (•) -y (q) i (•) 2 ν = E   - k∈(N d -Γ ) α (q) i,k φ k (•) + k∈Γ ᾱ(q) i,k -α (q) i,k φ k (•) 2 ν   .
By orthonormality of the {φ k (•)} we obtain

E ȳ(q) i (•) -y (q) i (•) 2 ν = k∈(N d -Γ ) (α (q) i,k ) 2 + k∈Γ E ᾱ(q) i,k -α (q) i,k 2 . For k ∈ Γ , observe that E ᾱ(q) i,k = α (q) i,k , therefore E ᾱ(q) i,k -α (q) i,k 2 = Var ᾱ(q) i,k = 1 M Var S (q) i X i i:N φ k (X i i ) ≤ 1 M φ k (•) 2 ∞ E S (q) i X i i:N 2 ,
which leads to the desired result.

Before stating the next preliminary result we define the following σ-algebras generated by the cloud of simulations,

G i := σ(C i ) and G i:N := σ(C i , . . . , C N ). Lemma 10 We have k∈Γ E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) 2 ≤ E S (q,M ) i X i i:N -S (q) i X i i:N 2 , (3.9
)

k∈Γ E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) | G i+1:N 2 ≤ E S (q,M ) i X i i:N -S (q) i X i i:N 2 | G i+1:N . (3.10)
Proof We start with the proof of (3.9). Let β i,k be the projection coefficients in

L 2 ν (R d ) of the function h i (x) := E S (q,M ) i X i i:N -S (q) i X i i:N | X i i = x
over the functions {φ k (•)}. Since the latter form an orthonormal basis of L 2 ν (R d ) (Proposition 2) and owing to the tower property of conditional expectation, we have

β i,k =E h i (X i i:N )φ k (X i i ) = E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) . (3.11)
By the Parseval identity and the Jensen inequality, we get

k∈N d β 2 i,k = E h 2 i (X i i ) ≤ E E S (q,M ) i X i i:N -S (q) i X i i:N 2 | X i i = E S (q,M ) i X i i:N -S (q) i X i i:N 2 .
By recalling the representation (3.11) for β i,k , the above inequality readily leads to (3.9). We now prove (3.10). Similarly, denote by

β i,k (G i+1:N ) the L 2 ν (R d )-projection coefficients of the (random) function h i (x, G i+1:N ) = E S (q,M ) i X i i:N -S (q) i X i i:N | X i i = x, G i+1:N .
Remind of the independence of G i+1:N and X i i:N . Therefore, using similar arguments as before, we obtain

β i,k (G i+1:N ) = E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) | G i+1:N and k∈N d |β i,k (G i+1:N )| 2 ≤ E S (q,M ) i X i i:N -S (q) i X i i:N 2 | G i+1:N ,
from which (3.10) readily follows. We are done.

Proof of Theorem 2

Let

ε i := E ȳ(q,M) i (•) -y (q) i (•) 2 ν
. By definition of ȳ(q,M)

i (•) and ȳ(q) i (•),

ε i = E   k∈Γ ᾱ(q,M) i,k φ k (•) -y (q) i (•) 2 ν   = E   ȳ(q) i (•) -y (q) i (•) + k∈Γ ᾱ(q,M) i,k - ᾱ(q) i,k φ k (•) 2 ν   .
Then, using the triangular inequality for • ν norm and the inequality (a + b) 2 ≤ 2a 2 + 2b 2 , we obtain

ε i ≤ 2E ȳ(q) i (•) -y (q) i (•) 2 ν + 2E   k∈Γ ᾱ(q,M) i,k - ᾱ(q) i,k φ k (•) 2 ν   := ε (a) i + ε (b) i .
(3.12)

The first term ε (a) i is bounded taking into account Lemma 9:

ε (a) i ≤ 2E i , (3.13)
Regarding the second term ε

(b) i , the orthonormality of {φ k (•)} gives ε (b) i = 2 k∈Γ E ᾱ(q,M) i,k - ᾱ(q) i,k 2 = 2 k∈Γ Var ᾱ(q,M) i,k - ᾱ(q) i,k + 2 k∈Γ E ᾱ(q,M) i,k - ᾱ(q) i,k 2 := ε (b,c) i + ε (b,d) i . (3.14)
To handle ε

(b,d) i write ᾱ(q,M) i,k - ᾱ(q) i,k = 1 M M m=1 S (q,M ) i X i,m i:N -S (q) i X i,m i:N φ k (X i,m i ),
and apply Lemma 10, thus obtaining

ε (b,d) i = 2 k∈Γ E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) 2 ≤ 2E S (q,M ) i X i i:N -S (q) i X i i:N 2 . (3.15)
Now we focus on ε (b,c) i

. We recall the conditional variance formula Var

[X] = E [Var [X | G]] + Var [E [X | G]]
, available for any scalar random variable X and any sigma-algebra G: then,

ε (b,c) i = 2 k∈Γ E Var ᾱ(q,M) i,k - ᾱ(q) i,k | G i+1:N + 2 k∈Γ Var E ᾱ(q,M) i,k - ᾱ(q) i,k | G i+1:N := ε (b,c,e) i + ε (b,c,f ) i .
We start with ε

(b,c,e) i : ε (b,c,e) i = 2 k∈Γ E Var 1 M M m=1 S (q,M ) i X i,m i:N -S (q) i X i,m i:N φ k (X i,m i ) | G i+1:N = 2 k∈Γ E 1 M 2 M m=1 Var S (q,M ) i X i,m i:N -S (q) i X i,m i:N φ k (X i,m i ) | G i+1:N ≤ 2 k∈Γ E 1 M 2 M m=1 φ k (•) 2 ∞ E S (q,M ) i X i,m i:N -S (q) i X i,m i:N 2 | G i+1:N ≤ 2L Γ M E S (q,M ) i X i i:N -S (q) i X i i:N 2 .
We continue with ε (b,c,f ) i

: use Lemma 10 to write

ε (b,c,f ) i = 2 k∈Γ Var E ᾱ(q,M) i,k - ᾱ(q) i,k | G i+1:N ≤ 2 k∈Γ E E ᾱ(q,M) i,k - ᾱ(q) i,k | G i+1:N 2 = 2E k∈Γ E S (q,M ) i X i i:N -S (q) i X i i:N φ k (X i i ) | G i+1:N 2 ≤ 2E E S (q,M ) i X i i:N -S (q) i X i i:N 2 | G i+1:N = 2E S (q,M ) i X i i:N -S (q) i X i i:N 2 .
Summarizing the above computations, we get a bound for ε

(b,c) i : ε (b,c) i ≤ 2 L Γ M + 1 E S (q,M ) i X i i:N -S (q) i X i i:N 2 .
All in all, in view of (3.14)-(3.15) and the above, we have derived a bound on ε

(b) i : ε (b) i ≤ 2 L Γ M + 2 E S (q,M ) i X i i:N -S (q) i X i i:N 2 . (3.16)
It remains to estimate the above term with the difference of the responses. Using the Cauchy-Schwarz inequality for the summation in j, and taking into account that f is Lipschitz in y, we obtain S (q,M ) i

X i i:N -S (q) i X i i:N 2 = 1 (1 + |X i i | 2 ) q N -1 j=i f j X i j , T L (ȳ (q,M ) j+1 (X i j+1 )(1 + |X i j+1 | 2 ) q/2 ) -f j X i j , y (q) j+1 (X i j+1 )(1 + |X i j+1 | 2 ) q/2 ∆ 2 ≤ (T -t i ) (1 + |X i i | 2 ) q N -1 j=i f j X i j , T L (ȳ (q,M ) j+1 (X i j+1 )(1 + |X i j+1 | 2 ) q/2 ) -f j X i j , y (q) j+1 (X i j+1 )(1 + |X i j+1 | 2 ) q/2 2 ∆ using that T L (y (q) j+1 (x)(1 + |x| 2 ) q/2 ) = y (q) j+1 (x)(1 + |x| 2 ) q/2
and that T L is 1-Lipschitz (see Proposition 1 and Definition (1.10))

≤ (T -t i ) (1 + |X i i | 2 ) q N -1 j=i L 2 f ȳ(q,M) j+1 (X i j+1 ) -y (q) j+1 (X i j+1 ) 2 (1 + |X i j+1 | 2 ) q ∆.
Combine this with (3.16) to obtain,

ε (b) i ≤ 2 L Γ M + 2 T L 2 f N -1 j=i E (1 + |X i j+1 | 2 ) q (1 + |X i i | 2 ) q ȳ(q,M) j+1 (X i j+1 ) -y (q) j+1 (X i j+1 ) 2 ∆.
Finally, taking advantage of the norm-stability property of the distribution ν (see (2.13)) we get,

ε (b) i ≤ 2 L Γ M + 2 T L 2 f c ν N -1 j=i ε j+1 ∆. (3.17) Let us set C ε := 2( L Γ M + 2)T L 2 f c ν .
At this point, in view of (3.12), (3.13) and (3.17) we have obtained

ε i ≤2E i + C ε N -1 j=i ε j+1 ∆.
Apply now [START_REF] Gobet | Approximation of backward stochastic differential equations using Malliavin weights and least-squares regression[END_REF]Proposition 2.4] with α = β = 1/2 to get the existence of a constant C(C ε , T ) (depending only on C ε , T ) such that

ε i ≤ 2E i + C(C ε , T )C ε N -1 j=i E j+1 ∆.
We are done.

Proof of Proposition 1

Without loss of generality, we can assume from now that η f = η g because assumptions on f and g are still valid by choosing the maximum parameter between η f and η g . Set η := η f = η g . First note that the constant C η in (2.5) is finite since b and σ are bounded, this follows from standard manipulations which are left to the reader.

Second we prove (2.4) by induction on i = N, . . . , 0. It is obvious for i = N since C y ≥ C g (note that C η ≥ 1). Now assume that sup

x∈R d |y j (x)| (1 + |x| 2 ) η 2 ≤ C j (3.18)
for j = N, . . . , i + 1 for some i ∈ {N, . . . , 0} and for some finite C j 's, and let us prove that the inequality holds also for j = i. In the proof below, we will get a relation between the C j 's leading to the statement of Proposition 1. We start from (2.7) and make use of the growth assumptions on f , g and the Lipschitz property of f : we have

|y i (x)| ≤ E C g (1 + |X N | 2 ) η 2 + N -1 j=i C f (1 + |X j | 2 ) η 2 + L f C j+1 (1 + |X j+1 | 2 ) η 2 ∆ | X i = x ≤ (1 + |x| 2 ) η 2   C g C η + T C f C η + C η L f N j=i+1 C j ∆   .
Therefore, (3.18) holds for j = i too, and thus owing to the induction principle, it holds for any j ∈ 0, . . . , N . Moreover, we can take

C i := C g C η + T C f C η + C η L f N j=i+1 C j ∆.
The Gronwall lemma implies that

C i ≤ C η (C g + T C f )e C η L f T .

GPU implementations

Two implementations on GPUs of Algorithm 1 are proposed in this section. The first one is a parallelization of the algorithm on a single GPU, we will refer to it as the (mono-)GPU version. This version includes two kernels, one simulates the paths of the forward process and computes the associated responses, the other one computes the regression coefficients (α

(q) i,k , k ∈ Γ ).
In the first kernel the initial value of each simulated path of the forward process is stored in a device vector in global memory, it will be read later in the second kernel. In order to minimize the number of memory transactions and therefore maximize performance, all accesses to global memory have been implemented in a coalesced way. The random numbers needed for the path generation of the forward process were generated on the fly (inline generation) taking advantage of the NVIDIA cuRAND library [9] and the generator MRG32k3a proposed by L'Ecuyer in [START_REF] L'ecuyer | Good parameters and implementations for combined multiple recursive random number generators[END_REF]. Therefore, inside this kernel the random number generator is called as needed. Another approach would be the pre-generation of the random numbers in a separate previous kernel, storing them in GPU global memory and reading them back from this device memory in the next kernel. Both alternatives have advantages and drawbacks. In this work we have chosen inline generation having in mind that this option is faster and saves global memory. Besides, register swapping was not observed on the implementation and the quality of the obtained solutions is similar to the accuracy of pure sequential traditional CPU solutions achieved employing more complex random number generators. Additionally, the pre-generation of random numbers is not feasible when dealing with a huge number of simulations M , as required for high dimensional problems, since it will require more memory than the available in the GPU. In the second kernel, in order to compute the regression coefficients, a parallelization not only over the multi-indices k ∈ Γ but also over the simulations 1 ≤ m ≤ M was proposed. Thus, blocks of threads parallelize the outer for loop ∀k ∈ Γ , whilst the threads inside each block carry out in parallel the inner loop traversing the vectors of the responses and the simulations. We stress that the parallel sum reduction (2.25) has been implemented in shared memory avoiding bank conflicts in pursuance of getting the most of the GPU.

A host can have several GPUs attached (in our machine, 4 GPUs TITAN BLACK). When dealing with high dimensional problems demanding a huge number of simulations M it is compelling to harness the computing power of these extra GPUs to further reduce execution times. The approach we have followed is to use OpenMP [37] to control all GPUs inside the host node. This strategy involves launching as many CPU threads as GPUs available on the node, so that each CPU thread handles one GPU. We will allude this version as M(ulti)GPU. In order to simulate the forward process and compute the responses, each GPU runs in parallel the first previous mono-GPU kernel, storing in its global memory the initial Student random variable of each path and the computed responses of the paths. Having done that, in order to compute the regression coefficients, the outer loop ∀k ∈ Γ is executed sequentially in CPU, whereas the inner loop on the simulations 1 ≤ m ≤ M is performed in parallel across all GPUs. More precisely, the parallel reduction operation is executed in two levels. The first reduction level is carried out inside each one of the GPUs, it has been performed using the reduce primitive of the NVIDIA Thrust library [START_REF]NVIDIA Thrust[END_REF]. The second level of reduction was made in the plane of CPU threads. This approach of parallelizing only over the simulations of the forward process behaves optimally when the number of simulations is large enough to fully occupy the GPU. In fact, this is the scenario we will be dealing with when tackling high dimensional problems. Actually, as sketched in Theorem 2, in order to get accurate approximations we need 2 d #Γ M to be small, therefore one need to choose M much bigger than #Γ . This fact motivates a small variant of this multi-GPU version, specifically designed for problems where M is so large that we cannot store the initial values of the simulated paths in the set of global memories of the GPUs. In the first kernel in charge of simulating the paths of the forward process and computing the associated responses, now the initial values of the forward process at each time step are not stored, thus saving d × M × 4 bytes of global memory1 . Instead, we just store the initial random seeds of the launched threads, which will be used in the next kernel in order to resimulate the same initial values of the forward process and allow the computation of the regression coefficients. This strategy is even as efficient as the previous one, on the grounds that GPUs are (most often) faster recomputing data than getting this data back from global memory. This is due to the fact that the version recomputing data is less memory bounded than the first one.

Numerical experiments

We introduce the following functions

g(x) := 1 + κ + sin λ d l=1 x l , f (t, x, y) := min   1, y -κ -1 -sin λ d l=1 x l exp λ 2 d(t -T ) 2 2   .
We aim to solve the PDE:

∂ t u(t, x) + 1 2 (∆ x u)(t, x) + f (t, x, u(t, x)) = 0, t ≤ T, u(T, x) = g(x).
(5.1)

Note that

u(t, x) = 1 + κ + sin λ d l=1 x l exp λ 2 d(t -T ) 2 ,
is solution to (5.1). The probabilistic formulation goes as follows:

dX t = dW t , -dY t = f (t, X t , Y t )dt -Z t dW t , Y T = g(X T ).
Then, the processes (Y, Z) satisfy Y t = u(t, X t ) and Z t = ∇ x u(t, X t ). For the numerical results the final time T is fixed to 1, κ = 6 10 and λ

= 1 √ d .
The µ parameter of the ν l Student's t-distribution is chosen as µ = 2. Under this choice ν l (x l ) = 0.5(1 + x 2 l ) -3/2 . The marginal CDF and its inverse are given by

F ν l (x) = 1 2 x √ x 2 + 1 + 1 , F -1 ν l (u) =                - -u 2 + u -0.25 u(u -1) if 0 < u ≤ 0.5, -u 2 + u -0.25 u(u -1) if 0.5 < u < 1.
(5.2)

In order to sample according to ν, we draw d independent random variables (U 1 , . . . , U d ) with uniform distribution on (0, 1), and we compute

X i i := F -1 ν 1 (U 1 ), . . . , F -1 ν d (U d ) d ∼ ν.
First we present some numerical results in dimension 1. Taking into account that

d r x y i (x) = ±λ exp t -1 2 [1 r even sin(x) + 1 r odd cos(x)] ,
in order to fit the assumptions of Lemma 3 we select p = r. The values of q are selected using relation (2.37) under p = r for r = 0, 1, 2 and 3. Therefore, for q = 0 only the function h (q) i is square integrable; for q = 2.1 the function h (q) i and its first spatial derivative are square integrable and the function h (q) i vanish for u → 0 + , 1 -; finally, for q = 8.1 the function h (q) i and its first two spatial derivatives are square integrable and vanish for u → 0 + , 1 -, its third spatial derivative is square integrable, as well.

In Figure 2, the analytical solution y (q) i versus the computed approximation ȳ(q,M) i are plotted for t = 0.5 and the values of q = 0, 2.1, 5.1, 8.1. The function was approximated on 1000 independent and identically ν distributed simulations. As long as q increases the convergence of the algorithm is improved (the truncation error shrinks to 0, coherently with Proposition 3). In Figure 3 the corresponding results are presented for h (q) i and h(q,M) i ; for q large the functions are very flat at the boundaries 0 and 1, as expected. The parameters of the GQRMDP algorithm were chosen ∆ = 0.01, K = 200 meaning that y (q) i was projected over 201 basis functions, and M = 2 × 10 6 . This selection makes 2 d #Γ M small, in fact it is 2.01 × 10 -4 , thus we control well the statistical error of the algorithm.

In order to assess the performance of the algorithm, we compute the average mean squared error (M SE) over 10 3 independent runs of the algorithm for two error indicators:

M SE max := ln    10 -3 max 0≤i≤N -1 10 3 m=1 y (q) i (R i,m ) - ȳ(q,M) i (R i,m ) 2    , M SE av := ln    10 -3 N -1 10 3 m=1 N -1 i=0 |y (q) i (R i,m ) - ȳ(q,M) i (R i,m )| 2    , (5.3)
Fig. 3 For t = 0.5, h (q) i versus h(q,M) i considering q = 0, 2.1, 5.1, 8.1. h(q,M) i were computed using ∆ = 0.01, K = 200, M = 2 × 10 6 . exponent q = 0, 2.1, 5.1, 8.1. As we have just seen in Figure 2, we expect (see Proposition 3) that the convergence of the algorithm is improved for large values of the damping exponent q. Besides, as long as the time step decreases, our parameter choice is such that K increases and 2K

M decreases so that the algorithm is more and more accurate. The GPU implementation provides a significant reduction in the computational time: the GPU speed-up reaches the value 12.2. One usually expects a GPU speed-up factor between 3x-10x vs two sockets of CPU [START_REF] Lee | Debunking the 100X GPU vs. CPU Myth: An Evaluation of Throughput Computing on CPU and GPU[END_REF], therefore both CPU and GPU codes were thoroughly optimized. In order to assess the accuracy of the proposed algorithm, in the following Table 3, for t = 0, q = 0 and different values of N , K and M , 99% confidence intervals for ȳ(q,M) i (0) are shown in the knowledge that u(0, 0) = 1.6. These confidence intervals were computed using 50 independent runs of the algorithm. Since we are computing the solution at the origin, similar results are obtained for values of q > 0.

q = 0 q = 2.1 q = 5.1 ∆ K M M SEmax M SEav M SEmax M SEav M SEmax M SEav 0.
Table 3 Now we present results in dimensions greater than one, for which the total degree and the hyperbolic cross index sets are different. Coherently with Section 2.4.3, we start showing the behaviour of the algorithm when the discussed conditions are enforced only on coordinate-wise derivatives of h (q) i . Table 4 shows results for d = 2. Convergence is clearly improved by increasing q, considering q = 5.1 (p = r = 2), not only h (q) i and its first two coordinate-wise derivatives are square integrable but also h (q) i and its first coordinate-wise derivatives vanish at the boundaries of [0, 1]2 . Additionally, convergence is enhanced as long as ∆ decreases: DEG and the number of simulations M increase in a way such that decreases. Finally, the hyperbolic cross index set offers better performance than the total degree set, not only an improved rate of convergence but also much faster computational times. In fact, for ∆ = 0.01, the hyperbolic cross index set is not only almost 3.4 times faster than the total degree set (on GPU), but also more accurate achieving M SE av = -10.137 instead of M SE av = -9.107. The GPU speed-up with respect to a fully optimized CPU parallel version reaches the value 11.21. Tables 5 and6 show results for d = 3 and 5, respectively. As before, the convergence of the algorithm improves for large values of the damping exponent, so that we benefit from small approximation error (see Proposition 4 with the conditions (2.43)-(2.48)). Furthermore, the hyperbolic cross index set outperforms the total degree set (which can be expected from the equation (2.63)). Except for the cases with ∆ = 0.05, where there are not enough simulations to fully take advantage of four GPUs, the multi-GPU implementation provides a significant reduction in the computational time: the speed-up with respect to the mono-GPU version reaches the value 3.75. The ideal speed-up would be 4, although the multi-GPU parallelization introduces some overhead which makes this value unattainable. Finally, results for d = 7, 9 and 11 are shown in Table 7 considering only the scenario performing better, i.e. the multi-GPU version considering a large value of the damping exponent along with the hyperbolic cross index set for the selection of the basis functions on which to project y (q) i . Before the end of this section, the behavior of the algorithm is exhibited using larger values of the damping exponent q that account for the case where not only coordinate-wise derivatives of h (q) i but also its mixed derivatives are required to vanish at the boundaries of [0, 1] d and to be square integrable. Taking into account that ∂ (r,...,r) x y i (x) = ±λ exp t -1 2 [1 dr even sin(x) + 1 dr odd cos(x)] ,

d = 1, i = 0, 99% confidence intervals for ȳ(q,M) i (0). y (q) i (0) = 1.6. ∆ K M q = 0 0.
and under Λ = {0, . . . , r} d , in order to fit the assumptions of Lemma 7 we select p = dr. In the following numerical experiments r = 2 or 1 will be selected. In order to avoid too large values of q, which will hinder the computation of the damping factor (1 + |x| 2 ) q/2 when single and even double precision are considered, the µ parameter of the ν l Student's t-distribution is chosen as µ = 1 from now on. Under this choice we deal with the Cauchy distribution ν l (x l ) = π -1 (1 + x 2 ) -1 , whose CDF and its inverse are given by F ν l (x) = 0.5 + π -1 arctan(x) and F -1 ν l (u) = tan (πx -0.5π), respectively.

Tables 8 and9 show the behavior of the algorithm using a value of q satisfying the condition (2.60) for r = 2, which boils down to q > 2dr -0.5 in our benchmark scenario. As expected from equation (2.63) for the Fourier coefficients of h (q) i , the hyperbolic cross index set clearly outperforms the total degree set. In order to take the most of our runtime environment, the codes were executed taking advantage of the whole available GPUs and the computations were performed using single precision. Ultimately, the semi-linear PDE (5.1) is solved for dimension 5 and 7 in Table 10, and for dimension 9 and 11 in Table 11. The GQRMDP algorithm was executed using double precision and considering only the hyperbolic cross index set for the projection of y (q) i . For d = 5, 7, the so-discussed conditions are imposed for Λ = {0, 1, 2} d , while for d = 9, 11, in order to avoid larger values of q, these smoothness conditions are only forced for Λ = {0, 1} d . The error indicators and the computational time for d = 7 inside Table 10 are to be compared with those of Table 7: as it was anticipated, for q = 27.6 the algorithm obtains much better convergence, whereas the execution time is larger due to more accurate computations on double precision. Proof We only have to consider the case n ≥ 1. In order to compute ∂ n x l 1 (1 + |x| 2 ) q/2 we use the following Faà di Bruno's formula (as long as the partial derivative is only with respect to one variable we can extend in the following way the one-dimensional Faà di Bruno's formula (2.34)):

∂ n x l f (g(x)) = n m=1 1 m! d m u f (u) u=g(x) j∈Jm,n n! j 1 !j 2 ! • • • jm! m i=1 ∂ j i
x l g(x), Jm,n = {j = (j 1 , . . . , jm) ∈ N m + : j 1 + . . . + jm = n}.

(A.3)
Here, consider f (u) = u -q/2 and g(x) = 1 + |x| 2 : d m u f (u) = Cm,qu -q/2-m and ∂ j x l g(x) = 2x l 1 j=1 + 21 j=2 . Observe that in the Faà di Bruno's formula (A.3), the sum over j ∈ Jm,n will be made only for j i = 1 with m 1 := #{j i : j ∈ Jm,n, j i = 1} terms, and for j i = 2 with m 2 := #{j i : j ∈ Jm,n, j i = 2} terms. Owing to this property and by definition of Jm,n, we have m 1 + 2m 2 = n and m 1 + m 2 = m, which results in m 1 = 2m -n and m 2 = n -m. Invoking the form of the derivatives ∂ j i

x l g(x), it readily follows that m i=1

∂ j i x l g(x) = 2 m x m 1 l .
All in all, we deduce

∂ n x l 1 (1 + |x| 2 ) q/2 ≤ C n m= n/2 (1 + |x| 2 ) -q/2-m (1 + |x l |) 2m-n ≤ C n m= n/2
(1 + |x| 2 ) -q/2-m (1 + |x| 2 ) m-n/2 ≤ C(1 + |x| 2 ) -q/2-n/2 . Lemma 13 (Estimates on derivatives of where the positive constants C m,J depend on m and the matrix J under consideration. Note that non coordinate-wise derivatives of g(x) = 1+|x| 2 are all zero. Besides, coordinatewise derivatives of g are ∂ j

x l g(x) = 2x l 1 j=1 + 21 j=2 . For a given J, let a l ∈ N ∀l = 1, . . . , d denote the number of vectors in J with 1 in the l-th coordinate and zero anywhere else. It holds that m i=1

∂ j i x g(x) = 1 2m≥n Cx a 1 1 • • • x a d d with a 1 + . . . + a d = 2m -n, m i=1 ∂ j i x g(x) ≤ C(1 + |x| 2 ) a 1 /2 • • • (1 + |x| 2 ) a d /2 1 d l=1 a l =2m-n ≤ C(1 + |x| 2 ) m-n/2

  where #Γ T (DEG) = DEG+d d . -Hyperbolic cross index set defined by Γ H (DEG) := k ∈ N d : d l=1 max(k l , 1) ≤ DEG . In order to compute the number of elements in Γ H (DEG) we split this set in the following two disjoint sets: -Γ d 1 (DEG) := k ∈ N d , ∀l, k l = 0 : d l=1 k l ≤ DEG , -Γ 2 (DEG) := k ∈ N d , ∃l, k l = 0 : d l=1 max(k l , 1) ≤ DEG . On the one hand, let us compute #Γ d 1 (DEG). Firstly we define Υ d (DEG) := k ∈ N d , ∀l, k l = 0 : d l=1 k l = DEG .

Fig. 1

 1 Fig. 1 For d = 2 examples of total and hyperbolic sets. Note that Γ H (9) ⊂ Γ H (15) ⊂ Γ T (16).

γ 2 2 ( 3 . 2 )

 2232 (y)ν(y)dy which gives the announced result (2.13) with c ν := C 1+d/C ν (C (3.2) T ).

1 (1+|x| 2 ) q/ 2 )n x 1 ( 1 +

 12211 Let x ∈ R d , q ∈ R + , n = (n 1 , . . . , n d ) ∈ N d + .The following upper bound holds∂ |x| 2 ) q/2 ≤ C(1 + |x| 2 ) -q/2-n/2 . (A.4)Proof We only have to consider the case n ≥ 1. For such n, we use the following multivariate Faà di Bruno's formula derived from [7, Corollary 2.10] m × N d : j 1 + . . . + jm = n, j 1 , . . . , jm = 0}, (A.5)

  In this section, we let ν be the Student density (2.10)-(2.11) in dimension d. In the first part of the section only y i and its coordinate-wise derivatives are assumed to be smooth. In the second part, on top of that mixed derivatives of y i will be supposed smooth. Let us start with the first case. Lemma 5 Assume that y i is r times (r ∈ N) coordinate-wise continuously differentiable on R d , with the existence of a polynomial growth degree p ∈ N such that

	2.4.3 Multidimensional case
	lim sup |x|→+∞	|∂ n x l y i (x)| |x| p-n < +∞, for	n = 0, . . . , r, l = 1, . . . , d.	(2.42)
	Assume the relation
			p + rµ < q			(2.43)
						3)
	6	84	50	5	126	48
	8	165	86	7	330	136
	10	286	123	9	715	248
	12	455	165	11	1365	368

Table 2 d

 2 = 1.

Table 4 d

 4 = 2.

						q = 0	q = 5.1
		∆	DEG #Γ (•)	M	M SEmax	M SEav	M SEmax	M SEav	CPU	GPU
	Total	0.05 0.02 0.01	20 25 31	231 351 528	2 × 10 4 2 × 10 5 2 × 10 6	-2.322 -2.869 -2.978	-2.773 -3.688 -3.966	-4.646 -6.761 -8.793	-5.239 -7.186 -9.107 25144.70 2243.06 3.24 0.33 362.55 35.51
	Hyper.	0.05 0.02 0.01	19 24 30	99 2 × 10 4 133 2 × 10 5 172 2 × 10 6	-2.526 -2.745 -2.898	-3.062 -3.550 -3.770	-5.588 -7.487 -9.389 -10.137 -6.196 -8.191	1.43 129.99 7221.68	0.16 12.38 662.54

Table 5 d

 5 = 3.

						q = 0	q = 5.1
		∆	DEG #Γ (•)	M	M SEmax	M SEav	M SEmax	M SEav	GPU	MGPU
	Total	0.05 0.02 0.01	22 24 26	2300 2925 3654	2 × 10 5 2 × 10 6 2 × 10 7	-2.164 -2.423 -2.484	-2.785 -3.580 -3.763	-5.699 -7.693 -9.836 -10.487 220687.47 58783.53 -6.141 53.89 121.75 -8.329 4361.22 1351.36
	Hyper.	0.05 0.02 0.01	20 22 24	411 459 528	2 × 10 5 2 × 10 6 2 × 10 7	-2.255 -2.249 -2.248	-3.191 -3.356 -3.380	-6.883 -7.778 -8.387	-7.662 -8.875 -9.576	8.37 592.19 27317.65	22.99 191.59 7306.78

Table 6 d

 6 = 5.

						q = 0	q = 5.1
		∆	DEG #Γ (•)	M	M SEmax	M SEav	M SEmax	M SEav	MGPU
	Tot.	0.2 0.1	20 21	53130 2 × 10 6 65780 2 × 10 7	-1.906 -2.221	-2.049 -3.356	-6.537 -8.482	-6.994 -9.103 19163.72 874.36
	Hyper.	0.2 0.1 0.1	16 17 36	3042 3122 10503 2 × 10 7 2 × 10 6 2 × 10 7	-2.353 -1.955 -2.216	-3.259 -3.278 -3.494	-7.312 -7.438 -8.609	-7.835 -8.125 -9.407	49.27 773.70 2718.07
		0.1	66	24893 2 × 10 7	-2.381	-3.603	-8.819	-9.703	6563.77

Table 7 d

 7 = 7, 9, 11, hyperbolic cross index set.

	q = 5.1

Table 8 d

 8 = 2, Λ = {0, 1, 2} 2 .

						q = 7.6
		∆	DEG #Γ (DEG)	M	M SEmax	M SEav	MGPU
	Total	0.05 0.02 0.01	31 41 51	528 903 1378 2 × 10 7 2 × 10 5 2 × 10 6	-7.453 -8.916 -10.528 -11.416 13220.48 -7.871 27.98 -9.598 297.99
	Hyper.	0.05 0.02 0.01	30 40 50	172 239 308	2 × 10 5 2 × 10 6 2 × 10 7	-7.710 -9.019 -10.152 -8.494 -10.425 -11.964	9.94 73.66 2646.07

Table 9 d

 9 = 3, Λ = {0, 1, 2} 3 .

						q = 11.6
		∆	DEG #Γ (DEG)	M	M SEmax	M SEav	MGPU
	Tot.	0.2 0.2	142 152	497640 608685	2 × 10 6 2 × 10 7	-5.009 -6.848	-5.423 -7.085 27347.05 6950.81
	Hyper.	0.2 0.2 0.1	140 150 160	4881 2 × 10 6 5322 2 × 10 7 5763 2 × 10 8	-7.955 -8.173 -8.319	-8.693 -9.039 -9.494	68.73 205.92 8030.68

Table 10 d

 10 = 5, 7, hyperbolic cross index set, Λ = {0, 1, 2} d , double precision.

	d	∆	DEG #Γ H (DEG)	M	q	M SEmax	M SEav	MGPU
	5 0.1 100	45443 2 × 10 7	19.6	-9.065 -10.198 16889.41
	7 0.1	36	94460 2 × 10 7	27.6	-15.054 -16.381 48021.18

Table 11 d

 11 = 9, 11, hyperbolic cross index set, Λ = {0, 1} d , double precision.

	d	∆	DEG #Γ H (DEG)	M	q	M SEmax	M SEav	MGPU
	9 0.2	5	14336 2 × 10 7	17.6	-21.016 -21.852 2466.21
	11 0.2	2	13312 2 × 10 7	21.6	-25.559 -27.161 2648.61

Single precision data type is assumed.

[START_REF] Bally | A quantization algorithm for solving multi-dimensional optimal stopping problems[END_REF] #Γ M

The first author research is part of the Finance for Energy Markets (FiME) lab, of the Chair Financial Risks of the Risk Foundation and of the ANR project CAESARS (ANR-15-CE05-0024). The second author has been financially supported by the Chair Financial Risks of the Risk Foundation, the Spanish grant MTM2016-76497-R and the Xunta de Galicia 2018 postdoctoral grant. The third author was partially supported by Spanish grant MTM2016-76497-R.

where the simulations {R i,m ; i = 0, . . . , N -1, m = 1, . . . , 10 3 } are independent and identically ν-distributed, and independently drawn from the simulations used for the GQRMDP scheme.

In order to test the theoretical results we compare the performance of the algorithm according to the choice of the damping exponent q and the multi-indices set. We will show the impact of these choices on the convergence of the approximation of the semi-linear PDE and on the computational performance in terms of computational time and memory consumption.

The numerical experiments have been performed with the following hardware and software configurations: four GPUs GeForce GTX TITAN Black with 6 GBytes of global memory per device (see [START_REF]NVIDIA's Next Generation CUDA Compute Architecture[END_REF] for details in the architecture), two multicore Intel Xeon CPUs E5-2620 v2 clocked at 2.10 GHz (6 cores per socket) with 62 GBytes of RAM, CentOS Linux, NVIDIA CUDA SDK 7.5 and INTEL C compiler 15.0.6. The CPU programs were optimized and parallelized using OpenMP [37]. The multicore CPUs time (CPU) and the GPU time (GPU) will all be measured in seconds in the forthcoming tables. CPU times correspond to executions using 24 threads so as to take advantage of Intel Hyperthreading. Unless otherwise indicated, the results are obtained using single precision, both in CPU and GPU.

All examples will be run using 64 thread blocks, each with 256 threads. In the following Table 2 we show results for d = 1 and the choices of the damping

A Appendix. Technical estimates on Student's t-distribution

In this appendix, estimates on the CDF of the Student's t-distribution, its inverse and the derivatives of the inverse of the CDF are computed. As we have been doing in the article, we use the same notation C for any generic constants, we will keep the same C from line to line, in order the alleviate the computations, although its value changes. In this appendix ν l stands for Student's t-distribution density (2.11) with parameter µ, Fν l its marginal CDF and F -1 ν its inverse.

Lemma 11 (Estimates on Fν l and F -1 ν l ) Let x ∈ R, u ∈ (0, 1). The following estimates hold

Proof In view of (2.11) we know that ν l (x l ) ∼ x l →±∞ cµ |x l | µ+1 . Besides,

The distribution ν l being symmetric, we have Fν l (x) + Fν l (-x) = 1 and F -1 ν l (u) = -F -1 ν l (1u). Therefore, the estimate 1 -Fν l (x) for x → +∞ follows from the case x → -∞.

In order to compute F -1 ν l we use that Fν l (F -1

. Therefore,

, which is the advertised result. The case u → 1 -follows from the case for u → 0 + using again the symmetry of the distribution ν l .

Lemma 12 (Estimates on derivatives of

The following upper bound holds

for 2m ≥ n. All in all, we deduce

Lemma 14 (Estimates on derivatives of Student's t-distribution) Let n ∈ N. The following upper bound holds

Proof This result follows readily using the previous estimate (A.2).

Lemma 15 (Estimates on derivatives of F -1 ν l ) Let n ∈ N. The following upper bounds hold

where C is a constant depending on µ and n.

Proof Using the combinatorial formula for higher derivatives of inverses [START_REF] Johnson | Combinatorics of Higher Derivatives of Inverses[END_REF] and the fact that

Therefore, using (A.6),

.

Finally, combining the previous upper bound with (A.1) yields

where the constant C may change value along computations.