Marie Billaud-Friess

Arthur Macherey

Anthony Nouy

Clémentine Prieur

A probabilistic reduced basis method for parameter-dependent problems

Keywords: Reduced basis method, probabilistic greedy algorithm, parameter-dependent partial differential equation, Feynman-Kac formula 2010 AMS Subject Classifications: 65N75, 65D15

Probabilistic variants of Model Order Reduction (MOR) methods have recently emerged for improving stability and computational performance of classical approaches. In this paper, we propose a probabilistic Reduced Basis Method (RBM) for the approximation of a family of parameterdependent functions. It relies on a probabilistic greedy algorithm with an error indicator that can be written as an expectation of some parameter-dependent random variable. Practical algorithms relying on Monte Carlo estimates of this error indicator are discussed. In particular, when using Probably Approximately Correct (PAC) bandit algorithm, the resulting procedure is proven to be a weak greedy algorithm with high probability. Intended applications concern the approximation of a parameter-dependent family of functions for which we only have access to (noisy) pointwise evaluations. As a particular application, we consider the approximation of solution manifolds of linear parameter-dependent partial differential equations with a probabilistic interpretation through the Feynman-Kac formula.

Introduction

This article focuses on the approximation of a family of functions M = {u(ξ) : ξ ∈ Ξ} indexed by a parameter ξ, each function u(ξ) being an element of some high-dimensional vector space V . The functions u(ξ) can be known a priori, or implicitly given through parameter-dependent equations. In multi-query contexts such as optimization, control or uncertainty quantification, one is interested in computing u(ξ) for many instances of the parameter. For complex numerical models, this can be computationally intractable. Model order reduction (MOR) methods aim at providing an approximation u n (ξ) of u(ξ) which can be evaluated efficiently for any ξ in the parameter set Ξ. For linear approximation methods, an approximation u n (ξ) is obtained by means of a projection onto a low-dimensional subspace V n which is chosen to approximate at best M, uniformly over Ξ for empirical interpolation method (EIM) or reduced basis method (RBM), or in mean-square sense for proper orthogonal decomposition (POD) or proper generalized decomposition (PGD) methods (see, e.g., the survey [START_REF] Nouy | Low-Rank Tensor Methods for Model Order Reduction[END_REF]) .

Probabilistic variants of MOR methods have been recently proposed for improving stability and computational performance of classical MOR methods. In [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], the authors introduced a probabilistic greedy algorithm for the construction of reduced spaces V n , which uses different training sets in Ξ with moderate cardinality, randomly chosen at each iteration, that allows a sparse exploration of a possibly high-dimensional parameter set. In [START_REF] Cai | A Stochastic Discrete Empirical Interpolation Approach for Parameterized Systems[END_REF], the authors derive a similar probabilistic EIM using sequential sampling in Ξ, which provides an interpolation with a prescribed precision with high probability. Let us also mention that a control variate method using a reduced basis paradigm has been proposed in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF] for Monte Carlo (MC) estimation of the expectation of a collection of random variables u(ξ) in a space V of second-order random variables. A greedy algorithm is introduced to select a subspace V n of random variables, that relies on a statistical estimation of the projection error. This algorithm has been analyzed in [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF] and proven to be a weak greedy algorithm with high probability. Probabilistic approaches have also been introduced for providing efficient and numerically stable error estimates for reduced order models [START_REF] Homescu | Error Estimation for Reduced-Order Models of Dynamical Systems[END_REF][START_REF] Janon | Goal-oriented error estimation for parameterdependent nonlinear problems[END_REF][START_REF] Smetana | Randomized residual-based error estimators for the proper generalized decomposition approximation of parametrized problems[END_REF][START_REF] Smetana | Randomized residual-based error estimators for parametrized equations[END_REF]. In [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF][START_REF] Balabanov | Randomized linear algebra for model reduction. Part II: minimal residual methods and dictionary-based approximation[END_REF][START_REF] Balabanov | Preconditioners for model order reduction by interpolation and random sketching of operators[END_REF][START_REF] Saibaba | Randomized Discrete Empirical Interpolation Method for Nonlinear Model Reduction[END_REF][START_REF] Zahm | Interpolation of inverse operators for preconditioning parameterdependent equations[END_REF], random sketching methods have been systematically used in different tasks of projection-based model order reduction, including the construction of reduced spaces or libraries of reduced spaces, the projection onto these spaces, the error estimation and preconditioning.

Here, we consider the problem of computing an approximation u n of u within a reduced basis framework. The reduced basis method performs in two steps, offline and online. During the offline stage, a reduced space V n is generated from snapshots u(ξ i ) at parameter values ξ i greedily selected by maximizing over Ξ (or some subset of Ξ) an error indicator ∆(u n-1 (ξ), ξ) which provides a measure of the discrepancy between u(ξ) and u n-1 (ξ). Then, during the online step, u n (ξ) is obtained by some projection onto V n . In this paper, we propose a probabilistic greedy algorithm for which ∆(u n (ξ), ξ) is the square error norm u(ξ)u n (ξ) 2 V , expressed as the expectation of some parameter-dependent random variable Z n (ξ), ∆(u n (ξ), ξ) = E(Z n (ξ)).

(

) 1
For maximizing E(Z n (ξ)), we use a PAC (Probably Approximately Correct) bandit algorithm proposed by the authors in [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF], which relies on adaptive Monte-Carlo estimations of E(Z n (ξ)). The algorithm only requires a limited number of samples by preferably sampling random variables associated with a probable maximizer ξ. It is particularly suitable for applications where the random variable Z n (ξ) is costly to sample. Under suitable assumptions on the distribution of Z n (ξ), it provides a PAC maximizer in relative precision, meaning that with high probability the parameter ξ is a quasi-optimal solution of the optimization problem. We prove in this work that the resulting greedy algorithm is a weak-greedy algorithm with high probability.

Intended applications concern the approximation of a parameter-dependent family of functions u(ξ) defined on a bounded domain D for which we have access to (possibly noisy) pointwise evaluations u(ξ)(x) := u(x, ξ) for any x ∈ D. The proposed probabilistic greedy algorithm can be used to generate a sequence of spaces V n and corresponding interpolations u n of u onto V n . Assuming u(ξ) ∈ L 2 (D) and we have a direct access to pointwise evaluations, the square error norm u(ξ)-u n (ξ) 2

L 2 (D)
used to select the parameter ξ can be estimated from samples of Z n (ξ) = |D||u(Y, ξ)u n (Y, ξ)| 2 with Y a uniform random variable over D. It results in a probabilistic EIM in the spirit of [START_REF] Cai | A Stochastic Discrete Empirical Interpolation Approach for Parameterized Systems[END_REF]. In a fully discrete setting where Ξ and D are finite sets, u can be identified with a matrix and the proposed algorithm is a probabilistic version of adaptive cross approximation for low-rank matrix approximation [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF][START_REF] Tyrtyshnikov | Incomplete cross approximation in the mosaic-skeleton method[END_REF], with a particular column-selection strategy. Another context is the solution of a linear parameter-dependent partial differential equation (PDE) defined on a bounded domain D and whose solution u(ξ) admits a probabilistic representation through the Feynman-Kac formula. This allows to express a pointwise evaluation u(x, ξ) as the expectation of a functional of some stochastic process.

The problem being linear, the error u(ξ)u n (ξ) also admits a Feynman-Kac representation, which again allows to express the square error norm ∆(u n (ξ), ξ) = u(ξ)u n (ξ) 2 L 2 (D) as the expectation of some random variable Z n (ξ) and to estimate it through Monte-Carlo simulations of stochastic processes. This is a natural framework to apply the proposed probabilistic greedy algorithm, which allows a direct estimation of the targeted error norm and avoids the use of possibly highly biased residual based error estimates. In practice, as the exact solution of the PDE is not available, the snapshots used for generating the reduced space V n are numerical approximations computed from pointwise evaluations of the exact solution u(ξ) by some interpolation or learning procedure. This results in a fully probabilistic setting which opens the route for the solution of high-dimensional PDEs (see, e.g. [START_REF] Billaud-Friess | Stochastic methods for solving highdimensional partial differential equations[END_REF] where the authors rely on interpolation on sparse polynomial spaces). This paper is structured as follows. In Section 2 we recall basic facts concerning reduced basis method. Then in Section 3 we present and analyze our new probabilistic greedy algorithm. Based on this algorithm, we derive in Section 4 a new reduced basis method for parameter-dependent PDEs with probabilistic interpretation. Numerical results illustrating the performance of the proposed approaches are presented in Section 5.

Reduced basis greedy algorithms

As discussed in the introduction, reduced basis method relies on two steps. We mainly focus on the offline stage during which the reduced subspace V n ⊂ V is constructed. In particular, we recall in this section some basic facts concerning greedy algorithms usually considered in that context. For detailed overview on that topic see, e.g., surveys [START_REF] Haasdonk | Reduced Basis Methods for Parametrized PDEs-A Tutorial Introduction for Stationary and Instationary Problems, chapter 2[END_REF][START_REF] Nouy | Low-Rank Methods for High-Dimensional Approximation and Model Order Reduction[END_REF]. Throughout this paper, V is some Hilbert space equipped with a norm • V . We seek an approximation u n (ξ) of u(ξ) in a low-dimensional space V n which is designed to well approximate the solution manifold M = {u(ξ) : ξ ∈ Ξ}.

A benchmark for optimal linear approximation is given by the Kolmogorov n-width

d n (M) V := inf dim Vn=n sup u∈M u -P Vn u V ,
where the infimum is taken over all n-dimensional subspaces V n of V and where P Vn stands for the orthogonal projection onto V n . However, an optimal space V n is in general out of reach. A prominent approach is to rely on a greedy algorithm for generating a sequence of spaces from suitably selected parameter values. Starting from V 0 = {0}, the n-th step of this algorithm reads as follows. Given {ξ 1 , . . . , ξ n-1 } ⊂ Ξ and the corresponding subspace

V n-1 = span{u(ξ 1 ), . . . , u(ξ n-1 )}, a new parameter value ξ n is selected as u(ξ n ) -u n-1 (ξ n ) V = sup ξ∈Ξ u(ξ) -u n-1 (ξ) V , (2)
where u n-1 stands for an approximation of u(ξ) in V n-1 . However, this ideal algorithm is still unfeasible in practice, at least for the two following reasons:

1) computing the error u(ξ)u n-1 (ξ) V for all ξ ∈ Ξ may be unfeasible in practice (e.g. when u(ξ) is only given by some parameter-dependent equation), and

2) maximizing this error over Ξ is a non trivial optimization problem.

Point 1) is usually tackled by selecting a parameter ξ n which maximizes some surrogate error indicator ∆(u n-1 (ξ), ξ) that can be easily estimated. Assuming u n (ξ) is a quasi-optimal projection of u(ξ) onto V n and assuming ∆(u n-1 (ξ), ξ) is equivalent to u(ξ)u n-1 (ξ) V , there exists γ ∈ (0, 1] such that u(ξ n ) -

P V n-1 u(ξ n ) V ≥ γ sup ξ∈Ξ u(ξ) -P V n-1 u(ξ) V , (3)
which yields a weak-greedy algorithm. Quasi-optimality means that the approximation

u n of u in V n satisfies u(ξ) -u n (ξ) V ≤ C u(ξ) -P Vn u(ξ) V (4)
for some constant C independent from V n and ξ. Although the generated sequence V n is not optimal, it has been proven in [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF][START_REF] Devore | Greedy Algorithms for Reduced Bases in Banach Spaces[END_REF] that the approximation error

σ n (M) V := sup u∈M u -P Vn u V
has the same type of decay as the benchmark d n (M) V for algebraic or exponential convergence.

Remark 2.1. In the case of parameter-dependent linear equation arising e.g., from the discretization of some parameter-dependent linear PDE of the form r(u(ξ), ξ) = 0 with u(ξ) ∈ V = R N , the approximation u n (ξ) is typically obtained through some (Petrov-)Galerkin projection onto V n , with a complexity depending on n N . In such a context, a weak greedy algorithm classically involves a certified residual based error estimate ∆(u n (ξ), ξ), that is an upper bound of the true error. However, for some applications, such an error estimate can be pessimistic (when the underlying discrete operator is badly conditionned) so that the generated sequence V n is far from being optimal. A possible strategy to improve such an estimate is to consider a preconditioned residual [START_REF] Cohen | Adaptivity and variational stabilization for convectiondiffusion equations[END_REF][START_REF] Zahm | Interpolation of inverse operators for preconditioning parameterdependent equations[END_REF][START_REF] Balabanov | Preconditioners for model order reduction by interpolation and random sketching of operators[END_REF]. In Section 4, we overcome this limitation by considering for ∆(u n (ξ), ξ) the targeted square error norm u(ξ)u n (ξ) 2 V , which is evaluated using adaptive Monte-Carlo estimations. Point 2) is addressed by transforming the continuous optimization problem over Ξ into a discrete optimization over a finite subset Ξ ⊂ Ξ. Choosing the training set Ξ is a delicate task. As pointed out in [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]Section 2], if Ξ is an ε-net of Ξ, then a greedy algorithm for the approximation of the discrete solution manifold M = {u(ξ) : ξ ∈ Ξ} generates a sequence of spaces that are able to achieve a precision in O( ) with similar performance as the ideal greedy algorithm. However, the cardinal of Ξ may be very large for a parameter set Ξ in a high-dimensional space R p and when a low precision ε is required. In [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], the authors propose a greedy algorithm which uses different training sets randomly chosen at each step. Under suitable assumptions on the approximability of the solution map ξ → u(ξ) by sparse polynomial expansions, training sets can be chosen of moderate size independent of the parametric dimension p.

To conclude this section, we give a practical deterministic (weak)-greedy algorithm that can be summarized as follows.

Algorithm 2.2 (Deterministic greedy algorithm). Let Ξ ⊂ Ξ be a discrete training set and V 0 = {0}. For n ≥ 1 proceed as follows.

(Step 1.) Select ξ n ∈ arg max ξ∈ Ξ ∆(u n-1 (ξ), ξ). (Step 2.) Compute u(ξ n ) and update V n = span{u(ξ 1 ), . . . , u(ξ n )}.
Usually, Algorithm 2.2 is stopped when ∆(u n (ξ), ξ) is below some target precision ε > 0 or for a given dimension n.

A probabilistic greedy algorithm

In this section, we motivate and present a probabilistic variant of Algorithm 2.2. Such an algorithm relies on the concept of Probably Approximately Correct (PAC) maximum. It is proven to be a weak greedy algorithm with high probability.

As a starting point for our work, we assume that the error estimator required at each step of Algorithm 2.2 admits the following form

∆(u n (ξ), ξ) := u(ξ) -u n (ξ) 2 V = E(Z n (ξ)), (5)
where Z n (ξ) is some parameter-dependent real valued random variable, defined on the probability space (Ω, F, P). Here, E(Z n (ξ)) is a probabilistic representation of the current square error u(ξ)u n (ξ) 2 V depending on the targeted applications as discussed in what follows.

Example 3.1 (Estimate of the norm of approximation error). Let suppose that u(ξ) belongs to V = L 2 (D) the Lebesgue space of square integrable functions defined on a bounded set

D ⊂ R d . If ∆(u n (ξ), ξ) = u(ξ) -u n (ξ) 2 L 2 , then (5) holds with Z n (ξ) = |D||u(Y, ξ) -u n (Y, ξ)| 2 where Y ∼ U(D) is a random variable with uniform distribution over D.
Example 3.2 (Greedy algorithm for control variate [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF][START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF]). Let suppose that we want to compute a MC estimate of the expectation of a parameter-dependent family of random variables u(ξ) belonging to a Hilbert space of centered second-order random variables. MC estimate is known to slowly converge with respect to the number of samples of u(ξ). Variance reduction techniques based on control variates are usually used to improve MC estimates. In [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF] the authors propose a RB paradigm to compute a control variate with a greedy algorithm of the form of Algorithm 2.2 where

∆(u n (ξ), ξ) = E(Z n (ξ)) with Z n (ξ) = |u(ξ) -u n (ξ)| 2 in (5).

Main algorithm

Solving the following optimization problem

ξ n ∈ arg max ξ∈ Ξ E(Z n-1 (ξ)) (6)
is in general out of reach, since E(Z n-1 (ξ)) is unknown a priori or too costly to compute. Then, we propose a greedy algorithm with an approximate solution of (6).

Algorithm 3.3 (Probabilistic greedy algorithm). Let Ξ ⊂ Ξ be a discrete training set and (λ n ) n≥1 ⊂ (0, 1) N + . Starting from V 0 = {0}, proceed, for n ≥ 1, as follows.

(Step 1.) Select ξ n ∈ S(Z n-1 (ξ), Ξ). (Step 2.) Compute u(ξ n ) and update V n = span{u(ξ 1 ), . . . , u(ξ n )}.
The question is now how to choose properly the set of candidate parameter values S(Z n-1 (ξ), Ξ)? In view of numerical applications, a first practical and naive approach is to seek ξ n maximizing the empirical mean, i.e.

S(Z

n-1 (ξ), Ξ) := arg max ξ∈ Ξ Z n-1 (ξ) K where Z n-1 (ξ) K = 1 K K i=1 (Z n-1 (ξ)) i with K i.i.d. copies of Z n-1 (ξ).
Despite its simplicity, it is well known that such an estimate for the expectation suffers from low convergence with respect to the number of samples leading to possible high computational costs especially if Z n (ξ) is expensive to evaluate. Moreover, nothing ensures that the returned (random) parameter ξ n is a (quasi-)optimum for [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF], almost surely or at least with high probability.

Instead, the so-called bandit algorithms (see, e.g., monograph [START_REF] Lattimore | Bandit Algorithms[END_REF]) are good candidates to address [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF]. Here, we particularly focus on PAC bandit algorithms that for each n return a parameter value ξ n which is a probably approximately correct (PAC) maximum in relative precision for E(Z n (ξ)) over Ξ (see [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF]). For a given ε ∈ (0, 1) and probability λ n ∈ (0, 1), letting ξ n ∈ arg max ξ∈ Ξ E[Z n-1 (ξ)], such an algorithm returns ξ n satisfying

P (E(Z n-1 (ξ n )) -E(Z n-1 (ξ n )) ≤ εE(Z n-1 (ξ n ))) ≥ 1 -λ n . (7)
We use the notation S(Z n-1 (ξ), Ξ) := PAC λn,ε (Z n-1 , Ξ) when ξ n satisfies [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF]. In practice, the adaptive bandit algorithm in relative precision introduced in [6, §3.2] is particularly interesting in the case where Z n (ξ) is costly to evaluate since it preferentially samples the random variable Z n (ξ) for the parameter values for which it is more likely to find a maximum. Hence, it outperforms the mean complexity of a naive approach in terms of number of generated samples. Appendix A gives a detailed presentation of such a PAC adaptive bandit algorithm. As stated in Proposition A.3, such an algorithm provides a PAC maximum in relative precision, that fulfills [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF], in the particular case where {Z n (ξ), ξ ∈ Ξ} are random variables satisfying some concentration inequality. The interested reader can refer to [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF] and included references for more details.

Analysis of a greedy algorithm with PAC maximum

Now, we propose and analyze a probabilistic greedy algorithm where the parameter ξ n is a PAC maximum in relative precision for [START_REF] Billaud-Friess | Stochastic methods for solving highdimensional partial differential equations[END_REF] i.e. satisfying [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF], at each step n.

At a step n of the Algorithm 3.3, the reduced space V n = span{u(ξ 1 ), . . . , u(ξ n )}, as well as the approximation u n (ξ) are no longer deterministic. Indeed, they are related to the selected parameters ξ 1 , . . . , ξ n depending themselves on the errors at the previous steps through i.i.d. samples of the random variables Z i (ξ) for all ξ ∈ Ξ and i < n (required during PAC selection of ξ n ). Now, we prove that Algorithm 3.3 is a weak greedy algorithm with high probability.

Theorem 3.4. Take (λ n ) n≥1 ⊂ (0, 1) N + such that n≥1 λ n = λ < 1, ε ∈ (0, 1) and Ξ ⊂ Ξ a discrete training set. Moreover, suppose that for n ≥ 1 the approximation u n of u in V n is quasi-optimal in the sense that it satisfies (3) with Ξ replaced by Ξ. Then, Algorithm 3.3 is a weak-greedy algorithm of parameter

√ 1-ε C , with probability at least 1 -λ, i.e. P u(ξ n ) -P V n-1 u(ξ n ) V ≥ √ 1 -ε C max ξ∈ Ξ u(ξ) -P V n-1 u(ξ) V , ∀n ≥ 1 ≥ 1 -λ. (8)
Proof. Let first introduce some useful notation. We denote by P n-1 (•) := P(•|Z <n ) the conditional probability measure with respect to Z <n that denotes the collection of random variables Z i (ξ) k for all ξ ∈ Ξ and i < n, where Z i (ξ) k are i.i.d. copies of Z i (ξ). The related conditional expectation is

E n-1 [•] = E [•|Z <n ]. Now, let A = ∩ n≥1 A n , each event A n being defined as A n := E n-1 [Z n-1 (ξ n )] -E n-1 [Z n-1 (ξ n )] ≤ εE n-1 [Z n-1 (ξ n )] , with ξ n ∈ arg max ξ∈ Ξ E n-1 [Z n-1 (ξ)]
. Then, at each step n of Algorithm 3.3, the parameter ξ n is a PAC maximum knowing Z <n i.e.

P n-1 (A n ) ≥ 1 -λ n . (9)
Finally, as u n-1 (ξ) is completely determined by all the steps before n (i.e. depending only on Z <n ), we have

u(ξ) -u n-1 (ξ) 2 V = E n-1 (Z n-1 (ξ)). ( 10
)
For all n ≥ 1, the quasi-optimality condition (4) and probabilistic representation (10) lead to

u(ξ n ) -P V n-1 u(ξ n ) 2 V ≥ 1 C 2 u(ξ n ) -u n-1 (ξ n ) 2 V = 1 C 2 E n-1 [Z n-1 (ξ n )] (11)
Moreover, if A holds we have for n ≥ 1

E n-1 [Z n-1 (ξ n )] ≥ (1 -ε)E n-1 [Z n-1 (ξ n )] ≥ (1 -ε) max ξ∈ Ξ u(ξ) -P V n-1 u(ξ) 2 V ( 12
)
by definition of ξ n . Thus, by combining ( 11) and ( 12) we have that A implies for all n ≥ 1

u(ξ n ) -P V n-1 u(ξ n ) V ≥ √ 1 -ε C max ξ∈ Ξ u(ξ) -P V n-1 u(ξ) V .
We now estimate P(A)

P(A) = 1 -P(A) ≥ 1 - n≥1 P(A n ) = 1 - n≥1 E 1 An = 1 - n≥1 E[E 1 An |Z <n P n-1 (An) ] ≥ 1 - n≥1 λ n ,
where the last inequality derives from [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF], which concludes the proof.

Remark 3.5. Theorem 3.4 proves that Algorithm 3.3 is a weak greedy algorithm, with probability 1λ, for the approximation of the discrete solution manifold M. Thus, the approximation error σ n ( M) has the same decay rate as d n ( M) for algebraic or exponential convergence. In the lines of [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], it is possible to consider also a fully probabilistic variant of Algorithm 3.3, in which a training set Ξ n randomly chosen is used at each step n of Algorithm 3.3 instead of Ξ. For a particular class of functions that can be approximated by polynomials with a certain algebraic rate, it can be proven that, for suitable chosen size of random training set Ξ n , the resulting algorithm is a weak greedy algorithm with high probability with respect to the continuous solution manifold M.

Reduced basis method for parameter-dependent PDEs with probabilistic interpretation

We recall that this work is motivated by the approximation, in a reduced basis framework, of a costly function u(ξ) : D → R defined on the spatial domain D ⊂ R d depending on the parameters ξ lying in Ξ ⊂ R p . Here we consider the problem where u is the solution of a parameter-dependent PDE with probabilistic interpretation.

Let D be an open bounded domain in R d . For any parameter ξ ∈ Ξ, we seek u(ξ) : D → R the solution of the following boundary value problem,

-A(ξ)u(ξ) = g(ξ) in D, u(ξ) = f (ξ) on ∂D, (13)
where f (ξ) : ∂D → R, g(ξ) : D → R are respectively the boundary condition and source term, and A(ξ) is a linear and elliptic partial differential operator.

Since the exact solution of ( 13) is not computable in general, it is classical to consider instead u h (ξ) an approximation in some finite dimensional space V h ⊂ V deduced from some numerical discretization of the PDE. Classical RBM applies in that context, relying on some variational principles to provide an approximation u n (ξ) of u h (ξ) in a reduced space V n ⊂ V h , of small dimension, obtained through greedy algorithm (see Section 2). Here, we overcome such a priori discretization of the PDE and directly address the approximation of the solution u(ξ) of ( 13). The key idea is to use the socalled Feynman-Kac representation formula that allows to compute pointwise estimates of u(ξ) for any x ∈ D. This particular framework rises the following practical questions. During the offline step, how to choose a computable error estimator ∆(u n (ξ), ξ) required in greedy algorithm and compute the snapshots required for generating the reduced basis and related reduced space V n ? During the online step, how to compute the approximation u n ?

To that goal, in this section, a probabilistic RBM using only (noisy) pointwise evaluations is presented. We first recall the Feynman-Kac formula in Section 4.1. In Section 4.2, we detail a probabilistic greedy algorithm for construction of the reduced space V n in this setting. Finally, in Section 4.3, we discuss possible approaches for computing the approximation u n (ξ).

Feynman-Kac representation formula for an elliptic PDE

In what follows, W = (W t ) t≥0 denotes a standard d-dimensional Brownian motion defined on the probability space (Ω, F, P) endowed with its natural filtration (F t ) t≥0 . For the sake of simplicity, the dependence to parameter ξ is omitted in the presentation of the Feynman-Kac formula.

Let us consider the boundary problem [START_REF] Cohen | Adaptivity and variational stabilization for convectiondiffusion equations[END_REF], where the partial differential operator A is given as

A = 1 2 d i,j=1 (σσ T ) ij ∂ 2 ∂x i ∂x j + d i=1 b i ∂ ∂x i . (14)
It is the infinitesimal generator associated to the parameter-dependent d-dimensional diffusion process X x = (X x t ) t≥0 , adapted to (F t ) t≥0 , solution of the following stochastic differential equation (SDE) As problem ( 13) is defined on a bounded domain, we define the first exit time of D for the process X x as τ x = inf {s > 0 :

dX x t = b(X x t )dt + σ(X x t )dW t , X x 0 = x ∈ D, (15)
X x s / ∈ D} . (17)
Also, we assume some regularity property on the spatial domain D and data.

(x) = E f (X x τ x ) + τ x 0 g(X x t )dt := E(F (x, X x )), ( 18
)
where X x is the stopped diffusion process solution of (15).

Theorem 4.5 allows to derive a probabilistic numerical method for the computation of pointwise MC estimate of u, see Appendix B.

Offline step

During the offline step, the probabilistic greedy algorithm presented in Section 3 is considered to construct the reduced space V n . The keystone of such an algorithm is the probabilistic reinterpretation of the error estimate ∆(u n (ξ), ξ) as in Equation [START_REF] Billaud-Friess | Stochastic methods for solving highdimensional partial differential equations[END_REF]. Using the Feynman-Kac representation formula, we show in Section 4.2.1 that it is possible to rewrite the square of the approximation error as an expectation. Then, in Section 4.2.2, we discuss possible strategies for practical implementation of such an algorithm.

Probabilistic error estimate

Let us assume that u n (ξ) is a linear approximation of u(ξ) in a given reduced space V n ⊂ V (e.g., obtained using Algorithm 3.3). We recall that u(ξ) ∈ C(D) ∩ C 2 (D) is the unique solution of [START_REF] Cohen | Adaptivity and variational stabilization for convectiondiffusion equations[END_REF] with the following probabilistic representation

u(x, ξ) = E F (x, X x,ξ , ξ) , x ∈ D, (19)
with X x,ξ the parameter-dependent stopped diffusion process solution of [START_REF] Friedman | Stochastic Differential Equations and Applications[END_REF]. In classical RB methods, the error estimate ∆(u n (ξ), ξ) used in Algorithm 2.2 is usually related to some suitable norm of the equation residual. Here, we follow another path by considering the L 2 -norm of the current approximation error e n (ξ

) = u(ξ) -u n (ξ), i.e. ∆(u n (ξ), ξ) = e n (ξ) 2 L 2 .
In what follows, we give a possible probabilistic reinterpretation of this error. Assuming that u n (ξ) is regular enough (the regularity being inherited from the snapshots), the error e n (ξ) := u(ξ)u n (ξ) is the unique solution, for all ξ in Ξ, of

-A(ξ)e n (ξ) = g n (ξ) on D, e n (ξ) = f n (ξ) on ∂D, (20)
where f n (ξ) := f (ξ)u n (ξ) and g n (ξ) = g(ξ) + A(ξ)u n (ξ). By Feynman-Kac representation theorem, for all ξ in Ξ, e n (ξ) is the unique solution of (20) in C(D) ∩ C 2 (D) and satisfies for all

x ∈ D e n (x, ξ) = E f n (X x,ξ τ x,ξ , ξ) + τ x,ξ 0 g n (X x,ξ t , ξ)dt := E F n (x, X x,ξ , ξ) , (21)
with X x,ξ the stopped diffusion process solution of [START_REF] Friedman | Stochastic Differential Equations and Applications[END_REF]. Then, we have the following probabilistic reinterpretation for e n (ξ) 2 L 2 .

Theorem 4.6. Let Y ∼ U (D) be uniformly distributed on D. Let W and W be two independent standard d-dimensional Brownian motions defined on (Ω, F, P) and independent of Y . For any x ∈ D, let X x,ξ and Xx,ξ be solutions of (15) with W , W respectively. Then we have for any ξ in Ξ Proof. We first recall

e n (ξ) 2 L 2 = |D|E (Z n (ξ)) , ( 22
)
with Z n (ξ) = F n (Y, X Y,ξ , ξ)F n (Y,
e n (ξ) 2 L 2 = D e n (x, ξ) 2 dx = |D|E e n (Y, ξ) 2 .
Since, for any x, X x,ξ and Xx,ξ are i.i.d. random processes, we have

E e n (Y, ξ) 2 = E E F n (Y, X Y,ξ , ξ))|Y 2 = E E F n (Y, X Y,ξ , ξ)|Y E F n (Y, XY,ξ , ξ)|Y , and e n (ξ) 2 L 2 = |D|E E F n (Y, X Y,ξ , ξ)F n (Y, XY,ξ , ξ)|Y .
Then by the law of iterated expectation we obtain [START_REF] Lattimore | Bandit Algorithms[END_REF].

A probabilistic greedy algorithm using pointwise evaluations

For the purpose of numerical applications, we can apply Algorithm 3.3 together with the error estimate [START_REF] Lattimore | Bandit Algorithms[END_REF] for the construction of the reduced space V n .

Sample computation. The samples of Z n (ξ) (as defined in Theorem 4.6) are generated from the functional F n and independent trajectories of the discrete diffusion process X Y,∆t , with Y ∼ U(D). The discrete diffusion process X Y,∆t is computed using a suitable time integration scheme (see Appendix B) for the stochastic ODE [START_REF] Friedman | Stochastic Differential Equations and Applications[END_REF].

Snapshot computation. Within Algorithm 3.3, the reduced space corresponds to V n = span{u(ξ 1 ), . . . , u(ξ n )}. However, the snapshots {u(ξ 1 ), . . . , u(ξ n )} are generally not available since it requires to compute the exact solution of (13) for parameter instances {ξ 1 , . . . , ξ n }. From Feynman-Kac formula [START_REF] Gobet | Sequential Control Variates for Functionals of Markov Processes on JSTOR[END_REF], it is possible to compute MC estimates u ∆t,M (x, ξ) of u(ξ) from independent realizations of the diffusion process X x,∆t starting from x ∈ D (as detailed in Appendix B). Then, a global numerical approximation can be computed in some finite dimensional linear space, e.g., by interpolation or least-square method, from these MC pointwise estimates. To compensate possible slow convergence of MC estimates, one can consider a sequential approach which uses the approximation error at each step as control variate in order to reduce the variance of MC estimates. Such a strategy has been initially proposed in [START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF][START_REF] Gobet | Sequential Control Variates for Functionals of Markov Processes on JSTOR[END_REF] for interpolation, and recently extended for high dimensional problems in [START_REF] Billaud-Friess | Stochastic methods for solving highdimensional partial differential equations[END_REF].

Projection computation. The projection u n (ξ) is obtained following the same procedure as for online step as described in forthcoming section.

Online step

Given the reduced space V n obtained during the offline stage, we now turn to the computation of the approximation u n (ξ) during the online step. For RBM to be of practical interest, the approximation u n (ξ) in V n obtained at a low complexity depending only on n. Here, following the same path as for the snapshots computation during offline stage, u n (ξ) can be computed by interpolation or a least-square method using MC estimates u ∆t,M (x, ξ) given by (36).

Numerical applications

The aim of this section is two folds. We first illustrate the feasibility of a greedy algorithm with probabilistic error estimate for the approximation of a parameter-dependent function from its pointwise evaluations. Then, we present some numerical experiment concerning the probabilistic RBM, discussed in Section 4, for the solution of parameter-dependent PDEs with probabilistic interpretation.

Approximation of parameter-dependent functions

Let us consider the problem of computing an approximation u n (ξ) of u(ξ), from its pointwise evaluations at given points in D. Particularly, we seek u n (ξ) as the interpolation of u(ξ) in the finite dimensional space V n ⊂ V , such that

u n (x i , ξ) = u(x i , ξ), x i ∈ Γ,
with Γ = {x 1 , . . . , x n } an unisolvant grid of suitably chosen interpolation points in D.

Procedures for the construction of V n

For constructing the space V n = span{u(ξ 1 ), . . . , u(ξ n )}, we compare different greedy procedures for the selection of the snapshots u(ξ i ), i = 1, . . . , n. First, we use the deterministic greedy Algorithm 2.2, for which ∆(u n-1 (ξ), ξ) is a numerical estimate of the L 2 -norm of the approximation error u(ξ)u n (ξ) L 2 (D) using some integration rule. This approach is confronted to probabilistic alternatives relying on probabilistic reinterpretation of the approximation error

u(ξ) -u n (ξ) 2 L 2 (D) = E(Z n (ξ)).
where p m -p , p = 2 and λ n = λ n . We also consider a variant relying on asymptotic concentration inequality as Central Limit Theorem (CLT), which overcomes the necessity of computing any bound for Z n (ξ) and defining the sequence (d m ) m≥1 . These probabilistic approaches are also compared to another naive approach, in which ξ n is chosen at random in Ξ (without replacement) at each step n of Algorithm 3. 3. In what follows, the deterministic approach is called D-Greedy, whereas the probabilistic ones using MC estimate and bandit algorithms are named MC-greedy and PAC-greedy relying on non asymptotic (Bounded) or asymptotic concentration inequalities (CLT). The last one is simply referred as to Random.

Z n (ξ) = |D||u n (X, ξ) -u(X, ξ)| 2 , X ∼ U(D),

Numerical setting

We perform some numerical tests with the methods discussed in the previous section for the approximation of the two subsequent functions

u(x, ξ) = 10x sin(2πxξ), x ∈ [0, 1], ξ ∈ [2, 4]
and, following [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF],

u(x, ξ) = (x + 0.1) 1/2 1 [0,ξ] (x) + 1 2 (ξ + 0.1) -1/2 (x -ξ) + (ξ + 0.1) 1/2 1 [ξ,1] (x), x ∈ [0, 1], ξ ∈ [0, 1].
The test case related to each function will be designated by (TC1) and (TC2), respectively.

For the numerical experiments, the training set Ξ is obtained using # Ξ = 300 equally spaced points in Ξ, similarly for D made from equally spaced points in D (10000 for (TC1), and 1000 for (TC2)). Then, the L 2 -norm of the approximation error is estimated by trapezium rule. Both deterministic and probabilistic greedy algorithms are stopped for given n = 20 for (TC1) and n = 30 for (TC2). The interpolation grid Γ is set to be the sequence of magic points [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF], with respect to the basis of V n . For the probabilistic procedure with naive MC estimate, we set K ∈ {1, 50}. Finally, for bandit algorithms, the stopping criterion is ε = 0.9 and the probability of failure is set to λ = 0.1.

Numerical results

Let us first study the quality of the approximations provided by the different approaches. Figures 1 and2 represent the evolution of the estimated expectation E ξ and maximum, with respect to ξ, of the approximation error u n (ξ)u(ξ) L 2 ( D) for (TC1) and (TC2). These estimates have been computed using 100 independent realizations of u n (ξ) obtained from uniform draws of ξ in Ξ. In that case, only one realisation of the probabilistic algorithms is performed for the comparison. For both test cases, D-greedy, MC-greedy and PAC-greedy procedures behave similarly with the same error decay with respect to n reaching a precision around 10 -14 for (TC1) and 10 -5 for (TC2). Let us mention, that for (TC2), the function to approximate has a slow decay of its Kolmorov n-width (see e.g. discussion in [8, Section 4.3.2]) which explains higher error for larger n. For the random approach, the selection of interpolation points is less optimal. For first iterates it behaves similarly as other approaches, but we observe that the approximation u n (ξ) is less accurate with n, from around 15 for (TC1) and 10 for (TC2) respectively. However, despite no guarantee on the optimality of the returned parameter ξ n , the PAC algorithm with asymptotic concentration inequality and especially the MC greedy algorithm, either with a single random evaluation of the error estimate (K = 1), lead to very satisfactory results with an error close to the deterministic interpolation approach for both test cases. Greedy procedure. We now turn to the study of the greedy procedures used for the selection of the snapshots. Figures 34represent the error estimate ∆(u n-1 (ξ), ξ) as well as the number of samples m n (ξ) required during the greedy selection of ξ n for deterministic and probabilistic greedy algorithms based on bandit algorithms for (TC1) and (TC2). These curves corresponds to one realisation of the probabilistic algorithms. First we observe that the parameters selected (indicated with the symbol * on the curves) by bandit algorithms do not necessary coincide with the maximum of EIM. Second, as expected the number of samples m n (ξ) is adapted for both algorithms resulting in higher sampling in the region where it is likely to find maximum. Globally, we observe that PAC-greedy (CLT) works quite similarly as PAC-greedy (Bounded). But the two approaches differ in terms of required number of samples. Indeed, CLT based approach only requires around a maximum of 10 -10 2 samples whereas the one based on concentration inequalities requires between 10 3 -10 5 samples. Note that, for the proposed test cases the PAC-greedy (CLT) is quite competitive in comparison to MC approach with fixed number of samples K. CPU times. Now, we briefly discuss the CPU computational costs of compared methods which are summarized in Table 1. The more costly approach is the PAC-Greedy (Bounded) method since it requires a high number of samples for constructing suitable confidence intervals used during adaptive bandit algorithm. At the opposite, the Random approach is the cheapest. Let mention, that the PAC-Greedy (CLT) method, is quite competitive with MC-Greedy with K = 50 and underlines the interest of using some adaptive procedure for snapshot selection. However, the more interesting trade-off between efficiency and accuracy is the MC-Greedy approach whose computational time remains close to the D-Greedy method when K = 1.

Method

In regard to these observations, in the next section the MC-greedy approach with few samples K for the error estimation will be retained in practice to reduce computational costs.

Parameter-dependent PDE

Now, we focus on the solution of parameter-dependent PDEs, as introduced in Section 4. Given Ξ = [0.005, 1], we seek u(ξ), ξ ∈ Ξ, solution on D =]0, 1[ of the following boundary problem

-A(ξ)u(ξ) := -a(ξ)u (ξ) -b(ξ)u (ξ) = g(ξ) in ]0, 1[, u(ξ) = f (ξ) at x ∈ {0, 1}, (23)
for given boundary values f (ξ) : {0, 1} → R and source term g(ξ) : [0, 1] → R. Moreover, we denote a(ξ) ∈ R * + and b(ξ) ∈ R the diffusion and advection coefficients respectively. We assume that (23) admits a unique solution in C 2 ([0, 1]) whose probabilistic representation is given by

u(x, ξ) = E(F (x, X x,ξ , ξ)) := E f (X x,ξ τ x,ξ ) + τ x,ξ 0 g(X x,ξ t , ξ)dt .
The associated parameter-dependent stopped diffusion process X x,ξ is solution of

dX x,ξ t = b(ξ)dt + 2a(ξ)dW t , X x,ξ 0 = x. (24)
In the following we take a(ξ) = ξ and b(ξ) = -10. The source term as well as Dirichlet boundary conditions are set such that the exact solution to Equation ( 23) is

u(x, ξ) = exp(x/ξ) -1 exp(1/ξ) -1 • (25)

Compared procedures

In what follows, we test the probabilistic RBM discussed in Section 4. Since the exact solution is known, we use it for the snapshots. The projection u n is obtained from evaluations of the exact solution through interpolation (Interp) using magic points or a least-squares (LS) approximation using a set of points D in D, with # D ≥ n. It is also compared to the minimal residual based (MinRes) approximation defined by

u n (ξ) = arg min v∈Vn x i ∈ D |A(x i , ξ)v(x i , ξ) + g(x i , ξ)| 2 + |u(0, ξ) -v(0, ξ)| 2 + |u(1, ξ) -v(1, ξ)| 2 . ( 26
)
During the offline stage, we consider different greedy algorithms for the generation of the reduced spaces. First, we perform the proposed probabilistic greedy Algorithm 3.3 for the construction of the reduced space V n using MC estimates with Z n (ξ) defined as in Theorem 4.6. This approach is compared to an alternative RBM in deterministic setting. Since u(ξ) is implicitly known through the boundary value problem [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF], the greedy selection of V n is performed using Algorithm 2.2 where ∆(u n (ξ), ξ) is an estimate (using integration rules) of the residual based error estimate A(ξ)u n (ξ) + g(ξ) 2 L 2 (D) during the offline stage.

In the presented results, the residual based greedy algorithm is referred to as Residual. The MC estimate using Feynman-Kac representation is named FK-MC. These approaches are also compared to a naive one named Random in which the parameters ξ 1 , . . . , ξ n that define V n are chosen at random.

Numerical results

For the numerical experiments, the training set Ξ is made of 200 samples from a log uniform distribution over Ξ. This distribution is chosen as the solution strongly varies with respect to the viscosity ξ, in particular we want to reach small viscosities. Realizations of Z n (ξ), given by (36), are computed using M = 500 realizations of the approximate stochastic diffusion process solution of Equation ( 24), obtained by Euler-Maruyama scheme (see Appendix B) with ∆t = 10 -3 . For computing u n (ξ), magic points are used for interpolation while for LS and MinRes approaches, we choose for D a regular grid of 100 points in D. Here, greedy algorithms are stopped when n = 30. The MC-FK greedy algorithm is performed for K ∈ {1, 10}.

Figure 5 represents the estimated expectation and maximum, with ξ, of the approximation error for the compared procedures, with respect to n. The obtained results underline the importance of both projection and reduced basis construction. When MinRes method is used, the approximation is less accurate than for other deterministic approaches (dashed yellow curve). Especially, when it is compared to approaches using a residual based error estimate (blue curves) during the greedy procedure but using interpolation or LS for projection, the mean approximation error is of order 10 -11 for MinRes against 10 -15 for the latter (for the maximum error, we have 10 -10 against 10 -14 ). Second, let us comment the impact of the probabilistic reduced basis selection. As for function approximation, picking at random the reduced basis for the considered problem is far from optimal since the error expectation tends to stagnate around 10 -10 for n ≥ 20 (around 10 -9 for the maximum). However, when considering residual based or FK-based error estimates (even with K = 1), with interpolation or least-square ofr projection, the error behaves quite similarly and reaches 10 -15 for n = 30. This shows that the proposed probabilistic based error procedure performs well in practice.

Conclusion

In this work we have proposed a probabilistic greedy algorithm for the approximation of a family of parameter-dependent functions for which we only have access to (noisy) pointwise evaluations. It relies on an error indicator that can be written as an expectation of some parameter-dependent random variable. Different variants of this algorithm have been proposed using either naive MC estimates or PAC bandit algorithms, the latter leading to a weak greedy algorithm with high probability. Several test cases have demonstrated the performances of the proposed procedures.

For parameter-dependent PDEs whose solution admits a probabilistic representation, through the Feynman-Kac formula, such an algorithm can be embedded within a probabilistic RBM using only (noisy) pointwise evaluations. Numerical results have shown the main relevance of considering Feynman-Kac error based estimate for greedy basis selection. We have also illustrated the influence of the projection during online and offline step. Following the discussion of Section 4.2.2, using a sequential procedure as proposed in [START_REF] Billaud-Friess | Stochastic methods for solving highdimensional partial differential equations[END_REF][START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF][START_REF] Gobet | Sequential Control Variates for Functionals of Markov Processes on JSTOR[END_REF] should be an interesting alternative to avoid limitations of residual based projections. However, further work should be conducted to provide a projection with controlled error and at low cost, which is crucial for efficient model reduction.

A Adaptive bandit algorithms

We present an adaptive bandit algorithm to find a PAC maximum in relative precision of E(Z(ξ)) over the discrete set Ξ. Here {Z(ξ) : ξ ∈ Ξ} is a finite collection of random variables satisfying E[Z(ξ)] = 0, defined on the probability space (Ω, F, P). After introducing some required notation, we present a practical adaptive bandit algorithm which returns a PAC maximum in relative precision for [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF] when assuming suitable assumptions on the distribution of Z(ξ).

A.1 Notation and assumptions

We denote by Z(ξ) m the empirical mean of Z(ξ) and V (ξ) m its empirical variance, respectively defined by

Z(ξ) m = 1 m m k=1 Z(ξ) k and V (ξ) m = 1 m m k=1 Z(ξ) k -Z(ξ) m 2 , (27)
where Z(ξ) 1 , . . . , Z(ξ) m are m independent copies of Z(ξ). Moreover, the random variable Z(ξ) is assumed to satisfy the following concentration inequality

P |Z(ξ) m -E[Z(ξ)]| ≤ c(m, x, ξ) ≥ 1 -x, (28)
for each ξ ∈ Ξ, 0 ≤ x ≤ 1 and m ≥ 1. The bound c(m, x, ξ) depends on the probability distribution of Z(ξ).

For example, for second order random variables, the central limit theorem provides such a property with c(m, x, ξ) = γ x V (ξ) m m , and γ x the (x/2)-quantile of the normal distribution.

A.2 Algorithm

Now, let us define a sequence (d m ) m≥1 ⊂ (0, 1) N , independent from ξ, and such that m≥1 d m < ∞.

Then we introduce c ξ,m = c(m, d m , ξ), and β ± ξ,m(ξ) = Z(ξ) m(ξ) ± c ξ,m(ξ) . Note that the concentration inequality [START_REF] Smetana | Randomized residual-based error estimators for parametrized equations[END_REF]

Then, the adaptive bandit algorithm proposed in [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF] is as follows.

Algorithm A.2 (Adaptive bandit algorithm). for all ξ ∈ Ξ do Sample Z(ξ), increment m(ξ) and update ε ξ,m(ξ) .

5:

Compute the estimate Êm(ξ) [Z(ξ)] using (30).

6:

end for

7:

Increment and put in Ξ every ξ ∈ Ξ such that At each step , the principle of Algorithm A.2 is to successively increase the number of samples m(ξ) of the random variables Z(ξ) in the subset Ξ ⊂ Ξ, obtained using confidence intervals [β - m(ξ) (ξ), β + m(ξ) (ξ)] of E(Z(ξ)) according to (31). The idea behind is to use those confidence intervals to find regions of Ξ where one has a high chance to find a maximum. Then, ξ is returned as a maximizer over Ξ of the expectation estimate defined by [START_REF] Zahm | Interpolation of inverse operators for preconditioning parameterdependent equations[END_REF]. Under suitable assumptions, and when using certified non-asymptotic concentration inequalities, this algorithm returns a PAC maximum in relative precision of E(Z(ξ)) over Ξ, as stated in [6, Proposition 3.2], recalled below.

β + ξ,m(ξ) ≥ max ν∈ Ξ β - ν,m(ν) . (31)
which satisfies (32) for any p > 1.

In general, confidence intervals based on asymptotic theorems are much smaller than those obtained with non-asymptotic concentration inequalities, and yield a selection of much smaller sets Ξ of candidate maximizers, hence a much faster convergence of the algorithm. However, when using asymptotic theorems, we can not guarantee to obtain a PAC maximizer.

Assumption 4 . 1 .Assumption 4 . 2 .

 4142 where b(•) : R d → R d and σ(•) : R d → R d×d are the drift and diffusion coefficients, respectively. Before recalling Feynman-Kac formula, we introduce additional assumptions and notation. Denoting by • both euclidean norm on R d and Frobenius norm on R d×d , we first introduce the assumption that b and σ are Lipschitz continuous. There exists a constant 0 < M < +∞ such that for all x, y ∈ D we have b(x)b(y) + σ(x)σ(y) ≤ M xy . (16) Under Assumption 4.1, there exists a unique strong solution to Equation (15) (see e.g. [15, §5, Theorem 1.1.]).Denoting a = σσ T , we introduce the following uniform ellipticity assumption. There exists c > 0 such that y T a(x)y ≥ c y 2 , for all y ∈ R d , x ∈ D.

Assumption 4 . 3 .Assumption 4 . 4 .

 4344 The domain D is an open connected bounded domain of R d , regular in the sense that it satisfies P(τ x = 0) = 1, x ∈ ∂D. We assume that f is continuous on ∂D, g is Hölder-continuous on D.The following probabilistic representation theorem [15, §6, Theorem 2.4] holds. Theorem 4.5 (Feynman-Kac formula). Under Assumptions 4.1-4.4 there exists a unique solution of (13) in C(D) ∩ C 2 (D), which satisfies for all x ∈ D u

  XY,ξ , ξ) and |D| the Lebesgue measure of D.

Remark 4 . 7 .

 47 Assuming the existence of probabilistic representations for the gradient of u(ξ) and u n (ξ), it would be possible to consider probabilistic interpretation of other norms of the approximation error, such as the H 1 -norm. Such probabilistic representations have been derived in simple cases, see e.g.[START_REF] Gobet | Monte-Carlo Methods and Stochastic Processes: From Linear to Non-linear[END_REF] Corollary IV.5.2].

as discussed in Example 3 . 1 .

 31 In this setting, the set S(Z n-1 (ξ), Ξ) within Algorithm 3.3 is obtained using either a crude MC estimate of the expectation or adaptive bandit algorithms discussed in Appendix A. When non asymptotic concentration inequalities are used, the parameter ξ n returned by Algorithm 3.3 is a PAC maximum under suitable assumptions on the distribution of Z n (ξ). In particular, for any ξ ∈ Ξ, if there exist a n (ξ), b n (ξ) ∈ R such that a n (ξ) ≤ Z n (ξ) ≤ b n (ξ) a.s., concentration inequalities under the form (28) hold. Having such a knowledge a priori of the distribution of Z n (ξ) and finding the bounds a n (ξ), b n (ξ) is not an easy task. Here, since u(ξ) is known, we set the following heuristic bounds a n (ξ) = min x∈ D |u(x, ξ)u n (x, ξ)| 2 and b n (ξ) = max x∈ D |u(x, ξ)u n (x, ξ)| 2 , to perform our computations with D ⊂ D a finite subset. Moreover, by Remark A.4, we have to define the sequence (d m ) m≥1 with d m = λn # Ξ (p-1)

Figure 1 :

 1 Figure 1: (TC1) Evolution with respect to n, of the estimated expectation and maximum of the approximation error in L 2 -norm, computed for 100 instances of ξ, for one realisation of the probabilistic greedy algorithms compared to the deterministic one.

Figure 2 :

 2 Figure 2: (TC2) Evolution with respect to n, of the estimated expectation and maximum of the approximation error in L 2 -norm, computed for 100 instances of ξ, for one realisation of the probabilistic greedy algorithms compared to the deterministic one.

2 Figure 3 : 2 Figure 4 :

 2324 Figure 3: (TC1) Evolution of the error during greedy procedures based on PAC bandit algorithms (top), together with required samples m n (ξ) for selecting ξ n ∈ Ξ.

  K = 1 (Interp) FK-MC K = 1 (LS) FK-MC K = 10 (Interp) FK-MC K = 10 (LS)

Figure 5 :

 5 Figure 5: Parameter-dependent equation: evolution with respect to n, of the estimated expectation and maximum of the approximation error in L 2 -norm, computed for 100 instances of ξ, for one realisation of the probabilistic greedy algorithms compared to the deterministic one.

Remark A. 1 .

 1 In view of numerical experiments, we can consider standard concentration inequalities for sub-Gaussian or bounded random variables, see e.g.[6, §2]. In particular, if for any ξ ∈ Ξ, there exists a(ξ), b(ξ) ∈ R such that almost surely we have a(ξ) ≤ Z(ξ) ≤ b(ξ), then (28) holds withc(m, x, ξ) = 2V (ξ) m log(3/x) m + 3 (b(ξ)a(ξ)) log(3/x) m .An alternative to (28) is to rely on asymptotic confidence intervals for E(Z(ξ)) based on limit theorems of the form P |Z(ξ) m -E[Z(ξ)]| ≤ c(m, x, ξ) → 1x, as m → ∞.

  implies that [β - m(ξ) (ξ), β + m(ξ) (ξ)] is a confidence interval for E(Z(ξ)) with level 1-d m(ξ) . Letting s(ξ) := sign(Z(ξ) m(ξ) ) and ε ξ,m(ξ) = c ξ,m(ξ) |Z(ξ) m(ξ) | if Z(ξ) m(ξ) = 0, +∞ otherwise., we define the following estimate for E[Z(ξ)] given byÊm(ξ) [Z(ξ)] = Z(ξ) m(ξ)ε ξ,m(ξ) s(ξ)c ξ,m(ξ) if ε ξ,m(ξ) < 1, Z(ξ) m(ξ)otherwise.

1 :

 1 Let ε, λ ∈ (0, 1) and K ∈ N. Set = 0, Ξ 0 = Ξ, m(ξ) = K and ε ξ,m(ξ) = +∞ for all ξ ∈ Ξ. 2: while #Ξ > 1 and max ξ∈Ξ ε ξ,m(ξ) > ε 2+ε do 3:

8: end while 9 :

 9 Select ξ such that ξ ∈ arg max ξ∈Ξ Êm(ξ) [Z(ξ)] .

Proposition A. 3 .

 3 Let ε, λ ∈ (0, 1) and (d m ) m≥1 ⊂ (0, 1) be a sequence that satisfies m≥1 d m ≤ λ # Ξ and log(1/d m )/m → c satisfies (28). If for all ξ in Ξ, Z(ξ) is a random variable with E[Z(ξ)] = 0, then Algorithm A.2 almost surely stops and returns a PAC maximum in relative precision, i.e. ξ = PAC λ,ε (Z, Ξ). Remark A.4. In practice a possible choice for the sequence (d m ) m≥1 is to take d m = δcm -p with δ = λ # Ξ and c = p -1 p ,

Table 1 :

 1 CPU times for compared approaches.

	(TC1)	(TC2)

B Probabilistic approximation of the solution of a PDE Here we discuss the numerical computation of an estimate of u(x) for any x ∈ D. To that goal, we use a suitable integration scheme to get an approximation of the diffusion process X x and a MC method to evaluate the expectation in formula [START_REF] Gobet | Sequential Control Variates for Functionals of Markov Processes on JSTOR[END_REF].

An approximation of the diffusion process is obtained using a Euler-Maruyama scheme. More precisely, setting t n = n∆t, n ∈ N, X x is approximated by a piecewise constant process X x,∆t , where

where ∆W n = W n+1 -W n is an increment of the standard Brownian motion.

Numerical computation of u(x) for all x ∈ D requires the computation of a stopped process X x,∆t at time τ x,∆t , an estimation of the first exit time of D. Here, we consider the simplest way to define this discrete exit time τ x,∆t = min t n > 0 :

Such a discretization choice may lead to over-estimation of the exit time with an error in O(∆t 1/2 ). More sophisticated approaches are possible to improve the order of convergence, as Brownian bridge, boundary shifting or Walk On Sphere (WOS) methods, see e.g., [START_REF] Gobet | Monte-Carlo Methods and Stochastic Processes: From Linear to Non-linear[END_REF]Chapter 6]. These are not considered here.

Letting {X x,∆t (ω m )} M m=1 be M independent samples of X x,∆t , we obtain a MC estimate noted u ∆t,M (x) for u(x) defined as