Pablo & Milena

Francine & Yann

Marielle Félix

REMERCIEMENTS

Je tiens tout d'abord à remercier vivement Mireille Bossy et Benjamin Jourdain d'avoir accepté de rapporter ma thèse. Merci pour vos relectures attentives de mon manuscrit. Merci également à Pierre Etoré et Tony Lelièvre pour avoir accepter d'être les examinateurs de ma thèse. Merci en particulier à ce dernier d'avoir accepté de présider le jury. Merci à tou.te.s pour vos remarques et vos questions suite à la soutenance, j'ai pris un grand plaisir à échanger avec vous sur mes sujets d'étude. Espérons juste que la prochaine fois ces échanges se feront en face à face et pas derrière nos écrans.

Ensuite j'aimerai remercier mes encadrant.e.s de thèse. Tout d'abord Clémentine. Merci de ton accueil à Grenoble, dans un cadre vraiment propice au travail ... et aux activités de montagne diverses et variées que j'ai eu l'opportunité de pratiquer dans les massifs alentours. Merci ensuite à Marie pour ton encadrement bienveillant et ta grande disponibilité, que ce soit pour parler de différents aspects pratiques et théoriques de ma thèse ou pour m'aider à soumettre mon premier papier un 1er janvier. Enfin merci à Anthony. La liste des choses que je te dois est trop longue pour figurer intégralement ici, donc je vais aller à l'essentiel. Tout d'abord merci pour toutes les connaissances scientifiques que tu m'as partagées tout au long de ces trois années. Merci du temps que tu m'as accordé, pour écouter et répondre à mes questions, qu'elles soient scientifiques ou plus personnelles. Merci de m'avoir montré ce qu'est un chercheur passionné et impliqué, notamment auprès de ses doctorants. Encore merci à vous trois, vous avez su, chacun.e à votre manière m'apporter ce qui m'a permis d'aller au bout de ce travail.

Enfin si mes encadrant ;e.s ont largement contribué à l'aboutissement de cette thèse, ils ne sont pas les seul.e.s à m'avoir soutenu. Et sans ces autres personnes je ne serai pas arrivé à bon port aussi "facilement". Merci donc à toutes les personnes qui m'ont soutenu, de près ou de loin, au cours de ma thèse.

Je pense tout d'abord à la frange grenobloise des soutiens. Merci à l'équipe Airsea pour son accueil au top, moi qui débarquais à Grenoble sans y connaître personne. Merci 1 en particulier à Annie pour tout le soutien administratif. Réaliser une thèse à cheval entre deux lieux géographiques n'est pas toujours évident mais avec ton aide ça n'a jamais paru très compliqué. Bien sûr merci à mes co-bureau. Que ce soit Reda, fidèle au poste depuis la première heure et toujours disponible pour un "bon café", un match de foot le midi ou une bière en soirée. Ou Philo qui nous a rejoint en cours de route et avec qui les discussions, scientifiques évidemment, auront été au moins aussi animées que les soirées. Ou Mathieu pour quelques mois, lui qui a bien participé aux discussions "scientifiques" du bureau. Merci à tou.te.s les trois pour votre soutien, je sais qu'on se recroisera très bientôt pour fêter ma thèse, celle de Philo ou, je l'espère, la chute du système capitaliste actuel. Merci également aux autres doctorant.e.s du labo et de l'équipe, mention toute spéciale pour Sophie et Victor le mangeur de fleur ... votre bureau était bien mais le nôtre était mieux. Merci aussi à Olivier pour avoir encore resserré les liens entre permanent.e.s et doctorant.e.s dans l'équipe, liens déjà savamment entretenus par nos séminaires d'équipe très intéressants à tous les points de vue.

Merci bien sûr à tou.te.s les membres du LJK que j'ai côtoyé au cours de mes quelques mois grenoblois. Merci en particulier à Laurence, Fred et Cathy pour les coups de main (informatiques, administratifs, mais pas que). Dédicace spéciale aux 7 nains pour les séances de course à pied intenses, suivies de soirées encore plus intenses et de lendemain beaucoup plus compliqués (sauf pour Francky, ce qui restera toujours un mystère). J'en termine par les ami.e.s de Grenoble, toujours là pour un ravitaillement oenologique ou zytologique en soirée ou un aller-retour en montagne pour se changer les idées (ou se perdre pour les plus chanceux). Je pense bien sûr à

le monde académique, s'il y a bien quelqu'un qui le mérite, c'est toi. Merci à Mazen de m'avoir rejoint ces derniers mois, c'est toujours agréable une personne dans le bureau pour ne pas péter un câble tout seul, surtout en période de rédaction. Et je ne désespère pas que tu trouves un jour un boulot à la hauteur de tes capacités en recherche ... et que tu élèves ta sensibilité politico-socio-environnementales à de telles hauteurs.

Merci ensuite à tou.te.s les doctorant.e.s du LMJL. Les ancien.ne.s d'abord avec Caro, Solène, Matthieu, Germain, Maha, Zeinab, Fabien et Alexandre. Merci pour les bons moments passés en conférence, à une après-midi des doctorants ou au labo. Les derniers n'ont pas été aussi nombreux que je voulais, à cause de mes nombreux aller-retours à Centrale, mais je garderai en tête votre superbe accueil et j'espère qu'on se recroisera à d'autres occasions. Merci aussi aux autres, rencontré.e.s plus récemment, et notamment Adrien, Adrian, Antoine, Honorine et Jérôme. Merci pour les échanges à distance en période de covid et pour les quelques bières qu'on a eu le temps d'aller boire entre les différents confinements. Les deux dernier.e.s, vous récupérez le bureau 013, faites-en bon usage.

Merci aussi aux centralien.ne.s du labo et au département inof-maths dans son ensemble. Plus particulièrement, merci à Erwan et Cécile mon grand frère et ma grande soeur de thèse. Merci pour votre accueil, vos conseils et vos disponibilités quand je suis revenu sur Nantes et après. Merci aussi à Flavien, Martin et Léa pour les bons moments passés quand j'étais à Centrale. Mention spéciale à Anthony, rencontre récente mais heureuse. Bon courage pour la fin de ta thèse, je sais qu'on va vite se recroiser.

Enfin les copains et copines de Nantes. Bon je vous aurai saoulé avec ma thèse ces dernières années, maintenant c'est bon c'est fini, promis on en parle plus. Merci à Étienne pour le soutien actif et les coups de pied aux fesses quand il le fallait. Merci aussi de m'avoir montré la voie, quand j'ai vu que tu pouvais avoir une thèse, mon syndrome de l'imposteur a disparu. Merci aux autres pour les nombreux verres (d'eau) de soutien et séance intensives (mais pas très régulières pour moi) de course à pied et de vélo, en particulier Raphaëlle, Raphaël, Vincent, Constance et Ludovic. Petite dédicace à Gildas aussi, maintenant que je vais avoir du temps pour m'entraîner, plus question de perdre au squash. Et puis vu qu'on approche de la fin, il faut mentionner les meilleur.e.s. Thomas. Merci pour tout ... les dépannages d'ordinateur quand le mien m'a lâché, les (trop) nombreux verres de soutien et les coups de pied aux fesses salvateurs. Clairement si j'ai fini c'est grâce à ta phrase : "tu pourras aller dire à tous ceux qui disent de la merde dans les médias que tu viens faire pareil mais qu'en plus tu as une thèse". Enfin Marion. La collocataire parfaite avec qui faire une thèse. Franchement malgré les périodes pas faciles qu'on a traversées tous les deux (deux thésard.e.s en rédaction dans le même appartement et en temps de covid ça en fait des raisons de péter des câbles) et malgré les confinements multiples on ne s'est jamais entretué, et ça c'est grâce à ta tolérance (sisi) à toute épreuve. Merci pour toutes les discussions tard le soir (pas trop tôt le matin, bizarre) qu'on a pu avoir, sur tout et n'importe quoi et pour ton soutien aussi bien mental que logistique (on m'a soufflé que soutien physique pouvait être mal interprété donc j'ai changé). Pas sûr que j'y sois arrivé sans toi et le (brillant) exemple que tu m'as montré (en plus du soutien important sus-mentionné).

Une dernière catégorie de remerciements : les inclassables. A ma famille tout d'abord, merci pour votre soutien tout au long de ces trois années. En particulier à ma maman pour m'avoir poussé sans rien comprendre à mon sujet. A mon papa aussi pour avoir été persuadé toutes ces années qu'il pourrait comprendre ce que je fais en enseignant les maths dans le secondaire ... je ne pense pas avoir vulgarisé mon sujet aussi bien que tu le fais avec les maths en classe. Mais clairement si je me suis orienté dans les maths et vers l'enseignement je te le dois (au moins en partie). Merci aussi à mes soeurs qui ont eu des échéances importantes en même temps que moi. Je n'ai pas été trop là pour vous ces derniers temps, mais promis je vais me rattraper. Et finalement merci à mes grand-mères pour leur soutien même si "faire de la recherche" c'est pas tous les jours très clair et puis après tout "c'est quand que tu fais un vrai métier ?". Les oncles, tantes et cousin.e.s ne sont pas oublié.e.s, merci pour votre soutien, quelle que soit la forme qu'il ait pris.

Il y a aussi les copains qui ont suivi ça de près ou de loin. Hélène tout d'abord pour les multiples lettres de soutien échangées. Tes prises de recul m'auront été salutaires (enfin elles l'auront été pour que je finisse ma thèse). Merci pour les discussions sur l'après-thèse aussi, qui m'ont donné envie de finir vite pour me concentrer pleinement à l'après. Merci aux Sportings aussi, je l'ai finie ma thèse Roland-Garros ... et tout ça presque sans louper un match. Merci aux zouzous de Besançon et de Neuville, les (trop rares) petites parties de pêche, séances de jardinage ou soirées au cours de ces 3 dernières années ont toujours été un agréable bol d'air et un bon moyen de couper avec la thèse. Il s'en est passé des choses ces trois dernières années, maintenant on va avoir le temps de fêter ça.

Et une Petit(e) dernière pour la fin, toi qui a su m'apporter ton soutien et tes encouragements au cours des derniers mois passés ensemble et ce malgré tes échéances synchronisées avec les miennes. Merci Camille. Je suis content d'avoir désormais un peu de temps rien qu'à nous, libéré du poids de nos projets professionnels respectifs.

List of Figures

1. 1 Comparison between the estimation of the exit time τ t,x,∆t of D and the estimation of the exit time τ t,x,∆t 1 of D ∆t [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF], illustrated for a particular trajectory of X t,x,∆t . The domain D ∆t (delimited by ---) allows to estimate an exit position (and exit time) with a bias lower than the one using D (delimited by -). . . . . . . . . . . . . . . . .

List of Tables

Introduction (version française)

Cette thèse est consacrée à l'approximation de la solution d'équations aux dérivées partielles (EDPs) en grande dimension, éventuellement paramétrées. La résolution d'EDPs est un domaine très vaste et très ancien des mathématiques : les EDPs ont été analysées depuis la fin du 17 ème siècle par Newton et Leibniz. Alors qu'à cette époque, seuls des modèles simples avec des solutions analytiques étaient étudiés, les récents développements en informatique ont ouvert la possibilité d'aborder des EDPs plus compliquées par le biais de méthodes numériques. Plus récemment, les EDPs dépendant des paramètres ont été considérées pour étudier des problèmes plus complexes en optimisation, contrôle, analyse de sensibilité ou quantification d'incertitude.

La solution d'une EDP est une fonction multivariée u(x) dépendant de certaines variables physiques (spatiales, temporelles) x dans R d . Pour les EDPs paramétrées, la solution u(x; y) dépend de variables supplémentaires y ∈ Y ⊂ R p , appelées paramètres.

Les problèmes concrets pour lesquels les EDPs sont utilisées proviennent de différents domaines, par exemple en physique, en biologie ou en finance. Ces applications nécessitent généralement de pouvoir approcher avec une certaine précision la solution u d'une EDP. Pour les EDPs paramétrées, des approximations précises de u(•; y) pour de nombreuses valeurs y des paramètres sont généralement requises. Des approximations précises signifient généralement des modèles numériques fins, donc coûteux.

De plus, la complexité numérique de l'approximation des fonctions u(x) ou u(x; y) peut augmenter de manière exponentielle avec les dimensions d ou p. C'est ce qu'on appelle la malédiction de la dimension [START_REF] Oseledets | Breaking the curse of dimensionality, or how to use svd in many dimensions[END_REF][START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric pdes[END_REF]. Cette malédiction peut parfois être évitée, par exemple en utilisant les méthodes de Monte-Carlo pour l'estimation d'intégrales [START_REF] Novak | Tractability of Multivariate Problems: Standard information for functionals[END_REF]. Lorsqu'il s'agit d'approximation de fonctions, qui est un problème plus difficle qu'un problème d'intégration, il n'est pas toujours possible d'éviter cette malédiction en utilisant des méthodes aléatoires [START_REF] Kunsch | High-dimensional function approximation: Breaking the curse with monte carlo methods[END_REF].

Les EDPs en grande dimension donnent donc lieu à des problèmes d'approximation difficiles, aussi bien dans des cas indépendants des paramètres que dans des cas paramétrés. Pour les EDPs paramétrées, une solution consiste à utiliser un métamodèle, également appelé modèle de substitution, qui approche u et dont l'évaluation est peu coûteuse par rapport au coût d'évaluation du modèle numérique initial. Il est particulièrement pertinent d'utiliser des modèles de substitution lorsque des éval-uations de u(•; y) pour de nombreuses valeurs de y sont demandées. La construction d'un métamodèle pour les fonctions en grande dimension est aujourd'hui un sujet de grand intérêt, par exemple en quantification d'incertitude (UQ) [START_REF] Sudret | Surrogate models for uncertainty quantification: An overview[END_REF], pour des problèmes de contrôle ou d'optimisation.

De ce qui précède, au moins deux questions émergent, que nous aborderons dans ce manuscrit :

• Comment approcher la solution u(x) d'une EDP lorsque d est grand ?

• Comment approcher la variété des solutions M := {u(•; y) : y ∈ Y } par des sous-espaces de faible dimension lorsque p est grand ?

Pour aborder ces questions, nous supposerons que la solution u admet une certaine représentation probabiliste, connue sous le nom de représentation de Feynman-Kac, qui trouve son origine dans les travaux de Richard Feynman et Mark Kac dans les années 1960 [START_REF] Kac | On distributions of certain wiener functionals[END_REF]. Cette représentation est au coeur des contributions de cette thèse.

EDPs en grande dimension

Nous nous concentrons d'abord sur l'approximation de la solution u(x) d'une EDP en grande dimension. Nous sommes donc dans le cas où u : O ⊂ R d → R est une fonction en grande dimension (d 1), avec O un domaine physique (spatial, ou spatio-temporel).

Les méthodes numériques standards pour résoudre les EDPs comprennent notamment la méthode des éléments finis (FEM) [START_REF] Zienkiewicz | The finite element method[END_REF], les différences finies [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF] ou encore les volumes finis. Ces méthodes permettent d'approcher la solution globale u : O → R de l'EDP. Elles sont très efficaces pour résoudre des problèmes en dimension faible (disons d ≤ 3). Cependant, elles nécessitent souvent l'utilisation d'une discrétisation fine du domaine physique O, ce qui rend impossible leur utilisation pour les problèmes en grande dimension même pour les ordinateurs actuels. D'autres approches sans maillage ont été introduites plus récemment, comme [START_REF] Flyer | Radial basis function-generated finite differences: A mesh-free method for computational geosciences[END_REF]. Cependant, à notre connaissance, elles n'ont pas encore été adaptées pour des problèmes en grande dimension.

Pour lutter contre la malédiction de la dimension, des méthodes d'approximation par réseaux de tenseurs ou réseaux de neurones ont été proposées [4,[START_REF] Khoromskij | Tensor numerical methods for multidimensional pdes: theoretical analysis and initial applications[END_REF][START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF]. Des méthodes d'approximation parcimonieuses ont également été développées [START_REF] Bungartz | Sparse grids[END_REF]. Ces techniques ont été exploitées pour la résolutions d'EDPs en grande dimension, voir par exemple [START_REF] Shen | Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems[END_REF][START_REF] Shen | Efficient spectral sparse grid methods and applications to high-dimensional elliptic equations ii. unbounded domains[END_REF].

Pour en venir à des techniques approximation plus locales, par opposition aux techniques globales, nous rappelons qu'il est possible d'obtenir des évaluations ponctuelles de la solution en utilisant une représentation probabiliste de celle-ci. La représentation dite de Feynman-Kac mentionnée précédemment fait partie de ces représentations probabilistes. De plus, il n'y a pas de malédiction de la dimension à affronter puisque nous nous concentrons uniquement sur l'obtention d'évaluations ponctuelles. Cette représentation existe pour plusieurs classes d'EDP linéaires et non linéaires [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]. Elle permet d'écrire la solution u sous la forme

u(x) = E [Q(X x )] , ∀x ∈ O, (1)
où Q est une fonction dépendant des données de l'EDP et d'un certain processus de diffusion stochastique X x dépendant de l'EDP et du point auquel la solution u est évaluée. En utilisant une estimation de Monte-Carlo traditionnelle basée sur M échantillons i.i.d. de X x , nous avons accès à des évaluations ponctuelles bruitées de u en différents points x ∈ O. Une analyse standard montre que la précision de ces évaluations ponctuelles se comporte généralement en O( √ M -1 + ∆t α ) pour un certain α > 0, où ∆t est un paramètre d'intégration temporelle pour l'approximation du processus stochastique X x [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Lorsque l'on souhaite approcher précisément u, par exemple jusqu'à la précision machine, les coûts de calcul deviennent dès lors trop élevés, même pour les ordinateurs actuels.

Différentes méthodes permettent de réduire la variance des évaluations ponctuelles de la solution, parmi lesquelles on peut citer [START_REF] Giles | Multilevel monte carlo methods[END_REF][START_REF] Mascagni | -shell error analysis for "walk on spheres[END_REF]. Les auteurs dans [START_REF] Gobet | A spectral monte carlo method for the poisson equation[END_REF] ont introduit une méthode de réduction de variance par variables de contrôle globales en utilisant la représentation de Feynman-Kac. Dans [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF] ils ont décrit et analysé plus précisément l'utilisation de ces variables de contrôle. Ils construisent d'abord une approximation globale ũ de u dans un espace d'approximation de dimension finie P. Ils utilisent par exemple pour cela des méthodes d'interpolation reposant sur des estimations d'évaluations ponctuelles de u sur une grille unisolvante Γ P pour P, obtenues par des méthodes de Monte-Carlo pour la représentation probabiliste (1). Montrant ensuite que l'erreur e = u -ũ admet également une représentation probabiliste, les auteurs construisent une approximation globale ẽ de e basée sur des estimations d'évaluations ponctuelles obtenues par des méthodes de Monte-Carlo. Le processus peut ensuite être itéré plusieurs fois pour ũ + ẽ. Il en résulte des évaluations ponctuelles approchées de u sur la grille Γ P avec une erreur qui dépend de l'erreur d'intégration temporelle de l'erreur d'approximation dans P de la solution exacte u [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF].

Première Contribution

Lorsque d 1, l'espace d'approximation P dans lequel vivent ũ et ẽ, et la grille associée Γ P doivent être bien choisis. Ici, nous combinons la méthode des variables de contrôle de [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF] avec la stratégie d'interpolation parcimonieuse adaptative de [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF] qui construit de manière adaptative l'espace d'approximation P et en même temps la grille unisolvante associée Γ P . Nous proposons deux versions adaptatives de l'algorithme de [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. Étant donné une approximation globale ũk de u, notre premier algorithme construit au pas k une interpolation adaptative parcimonieuse ẽk de l'erreur u-ũk dans l'espace d'approximation P k (en utilisant sa grille d'interpolation associée Γ P k ). Notre deuxième contribution est un algorithme d'interpolation adaptative perturbée basé sur [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF], où à chaque étape k de l'algorithme, nous calculons une approximation ũ de u en utilisant l'algorithme utilisant des variables de contrôle de [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF], plutôt que l'interpolation sur P k de la solution exacte. Ce travail a été publié dans [START_REF] Billaud-Friess | Stochastic methods for solving high-dimensional partial differential equations[END_REF].

EDPs paramétrées et Optimisation Discrète

Nous nous concentrons désormais sur l'approximation de y ∈ Y → u(•; y) ∈ H, qui est la solution d'une certaine EDP paramétrée, où H est un espace de Hilbert de fonctions définies sur le domaine O. Les méthodes de réduction de modèle visent à réduire la complexité des modèles numériques complexes. Parmi ces méthodes, on retrouve les méthodes de base réduite qui sont très pertinentes pour certaines classes d'EDPs paramétrées et qui exploitent une bonne approximabilité de la variété des solutions M := {u(y) : y ∈ Y } par des espaces de faible dimension [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF][START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic pde's[END_REF].

Les méthodes de base réduite construisent des sous-espaces H n de dimension n dans H tels que la distance entre H n et M est proche de la n-largeur de Kolmogorov [START_REF] Pinkus | N-widths in Approximation Theory[END_REF] qui caractérise la plus petite distance entre un sous-espace de dimension n de H et M, c'est-à-dire la meilleure erreur d'approximation qui peut être atteinte par les outils d'approximation linéaire. Nous nous concentrons en particulier sur les algorithmes gloutons pour les méthodes de base réduite [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]. Dans un cadre idéalisé, les algorithmes gloutons construisent de manière adaptative une séquence imbriquée d'espaces réduits (H n ) n≥1 selon la procédure suivante. À une étape n ≥ 1, l'espace réduit courant H n-1 est enrichi avec un snapshot u(y n ) qui est sélectionné parmi toutes les éléments de M dont la distance à H n-1 est maximale. Cette distance est en pratique estimée à l'aide d'une estimation d'erreur a posteriori [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] peu coûteuse à calculer. Pour les EDPs paramétrées, une estimation d'erreur standard utilise la norme du résidu [START_REF] Haasdonk | Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition[END_REF][START_REF] Haasdonk | Reduced basis methods for parametrized pdes-a tutorial introduction for stationary and instationary problems[END_REF]. De telles estimations d'erreur dégradent généralement la performance de l'algorithme glouton, qui se transforme alors en un algorithme glouton sous-optimal [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]. Pour certains problèmes paramétrés, les méthodes d'estimation d'erreur basées sur le résidu donnent de mauvais résultats dans la construction de l'espace réduit. Des techniques alternatives ont donc été proposées, comme le préconditionnement [START_REF] Zahm | Interpolation of inverse operators for preconditioning parameter-dependent equations[END_REF] ou l'estimation d'erreur sans résidu [START_REF] Chen | A robust error estimator and a residualfree error indicator for reduced basis methods[END_REF].

Cependant, pour des applications pratiques, la maximisation de l'erreur d'approximation sur un ensemble continu Y reste impossible. C'est pourquoi on considère généralement des ensembles d'apprentissage discrets Ξ ⊂ Y et la sélection des snapshots est faite parmi les élements de la variété discrète des solutions M Ξ := {u(ξ) : ξ ∈ Ξ}. Le choix de ces ensembles d'apprentissage discrets est particulièrement difficile pour les espaces de paramètres de grande dimension, c'est-à-dire lorsque p 1. La sélection aléatoire de ces ensembles d'apprentissage semble être une technique intéressante lorsque la variété des solutions M peut être bien approchée par des espaces de faible dimension [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]. Il en résulte un algorithme glouton probabiliste.

La norme de l'erreur dans les algorithmes gloutons peut être écrite comme une espérance utilisant une représentation de Feynman-Kac de ladite erreur. Dans ce contexte, afin d'enrichir l'espace réduit avec un snapshot sélectionné à l'aide d'un algorithme glouton, nous devons résoudre un problème d'optimisation discrète du type suivant

max ξ∈Ξ E[Z(ξ)], (2)
où, dans le contexte des méthodes de base réduite, les (Z(ξ)) ξ∈Ξ sont des variables aléatoires dépendant de l'espace réduit. De tels problèmes d'optimisation discrète ont déjà été largement étudiés depuis les années 1950 [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF][START_REF] Lattimore | Bandit algorithms[END_REF], initialement pour des applications médicales. Les algorithmes de bandit utilisent des estimations de Monte-Carlo et des inégalités de concentration afin de fournir une solution quasioptimale de (2) avec une probabilité élevée. Le résultat de tels algorithmes est un paramètre ξ dans Ξ qui est probably approximately correct (PAC) [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF] dans un sens absolu, ce qui signifie que

P E[Z(ξ )] -E Z( ξ) ≤ ≥ 1 -λ,
pour des valeurs prescrites de dans (0, 1) et λ dans (0, 1), avec ξ ∈ arg max ξ∈Ξ E[Z(ξ)].

Seconde contribution

Les algorithmes de bandits existants ont besoin d'estimer intelligement afin de pouvoir être utilisés pour la maximization de l'erreur à posteriori dans des algorithmes gloutons. Pour contourner l'estimation de , la deuxième contribution de cette thèse est un algorithme probablement approximativement correct (PAC) en précision relative qui combine des estimations de Monte-Carlo avec une précision relative garantie [START_REF] Mnih | Empirical bernstein stopping[END_REF] avec un algorithme PAC [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF]. Cet algorithme fournit ξ dans Ξ qui satisfait

P E [Z(ξ )] -E Z( ξ) ≤ |E[Z(ξ )| ≥ 1 -λ, ( 3
)
pour des valeurs prescrites de dans (0, 1) et λ dans (0, 1), avec ξ ∈ arg max ξ∈Ξ E[Z(ξ)].

Troisième contribution

La troisième contribution de cette thèse combine les méthodes de base réduite) avec une sélection de snapshots basée sur le nouvel algorithme PAC en précision relative. En réinterprétant l'erreur d'approximation dans un algorithme glouton comme une espérance, nous pouvons sélectionner un snapshot pour enrichir la base réduite en utilisant l'algorithme PAC en précision relative. Il en résulte un algorithme glouton probabiliste. En combinant cet algorithme avec l'utilisation d'ensembles d'apprentissage aléatoires [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], qui sont adaptés au cas p 1, nous proposons un deuxième algorithme glouton probabiliste qui consiste à choisir aléatoirement un snapshot quasi-optimale parmi des ensembles d'apprentissage aléatoires en utilisant l'algorithme PAC en précision relative que nous avons proposé.

Plan du Manuscrit

Cette thèse est composée de cinq chapitres.

Dans un premier temps, nous nous concentrons sur la résolution d'EDPs en grande dimension à l'aide de méthodes probabilistes, ce qui correspond à notre première contribution, publiée dans [START_REF] Billaud-Friess | Stochastic methods for solving high-dimensional partial differential equations[END_REF]. Un aperçu des méthodes probabilistes pour les EDPs linéaires est présenté dans le Chapitre 1 tandis que notre contribution est présentée dans le Chapitre 2.

Le Chapitre 3 traite d'un problème d'optimisation discrète pour une fonction définie comme une espérance, et introduit notre deuxième contribution. Ce travail a été soumis comme [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF].

Enfin, nous étudions les EDPs paramétrées. Dans le Chapitre 4, nous donnons un aperçu des méthodes de base réduite pour ces modèles. Le Chapitre 5 introduit deux algorithmes gloutons probabilistes pour les méthode de base réduite, qui utilisent l'algorithme d'optimisation discrète introduit dans le Chapitre 3.

Introduction (english version)

This thesis is devoted to the approximation of the solution of high-dimensional partial differential equations (PDEs), possibly parameter-dependent. Approximating the solution of PDEs is a vast and quite old field of mathematics: PDEs have been analysed since the late 17 th century by Newton and Leibniz. While at this time only simple models with analytical solutions were studied, the recent developments in computer science have opened up the possibility to address more complicated PDEs through numerical methods. More recently parameter-dependent PDEs have been considered to study more complex problems in optimization, control, sensitivity analysis or uncertainty quantification.

The solution of a PDE is formally defined as a multi-variable function u(x) depending on some physical (spatial, temporal) variables x in R d . For parameter-dependent PDEs, the solution u(x; y) depends on additional variables y ∈ Y ⊂ R p .

The concrete problems for which PDEs are used arise in various fields, for example in physics, biology or finance. These applications usually require to be able to approximate up to some precision the solution u of a PDE. For parameter-dependent PDEs, precise approximations of u(•; y) for many instances y of the parameters are usually required. Precise approximations usually means fine, thus costly, numerical models.

Moreover, the numerical complexity to approximate functions u(x) or u(x; y) might increase exponentially with the dimensions d or p. This is known as the curse of dimensionality [START_REF] Oseledets | Breaking the curse of dimensionality, or how to use svd in many dimensions[END_REF][START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric pdes[END_REF]. This curse can sometimes be avoided, for example using Monte-Carlo methods for the estimation of integrals [START_REF] Novak | Tractability of Multivariate Problems: Standard information for functionals[END_REF]. When it comes to function approximation, which is a more complex problem than integration, avoiding this curse is not always possible using random methods [START_REF] Kunsch | High-dimensional function approximation: Breaking the curse with monte carlo methods[END_REF].

High-dimensional PDEs thus result in difficult approximation problems, in both parameter-independent and parameter-dependent cases. For parameter-dependent PDEs, a solution is to use a metamodel, also known as surrogate model, which approximates u and is cheap to evaluate compared to the initial numerical model. It is particularly relevant to use surrogates when evaluations of u(•; y) for many values of y are required. The building of metamodel for high-dimensional functions is nowadays a matter of high interest, e.g. in uncertainty quantification (UQ) [START_REF] Sudret | Surrogate models for uncertainty quantification: An overview[END_REF], control or optimisation.

From the foregoing, at least two questions arise, that we will address in this manuscript:

• How to approximate the solution x ∈ R d → u(x) of a PDE when d is large ?

• How to approximate the solution manifold M := {u(•; y) : y ∈ Y ⊂ R p } by low-dimensional subspaces of H when p is large ?

To tackle these questions, we will assume that the solution u admits some probabilistic representation, known as the Feynman-Kac representation, that originates in the work from Richard Feynman and Mark Kac in the 1960s [START_REF] Kac | On distributions of certain wiener functionals[END_REF]. This representation is at the core of the contributions of this thesis.

High-Dimensional PDEs

First we focus on the approximation of the solution u(x) to some high-dimensional PDE. We are thus in the case where u :

O ⊂ R d → R is a high-dimensional function (d 1)
, with O a physical domain (spatial or spatio-temporal).

Standard numerical methods for solving PDEs include Finite Element Method (FEM) [START_REF] Zienkiewicz | The finite element method[END_REF], finite differences [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF] or finite volumes. Such methods approximate the global solution u : O → R. They are really efficient for solving low-dimensional problems (say d ≤ 3). However they often require the use of a fine discretization of the physical domain O, which makes them intractable for high-dimensional problems, even for today's computers. Other meshless approaches have been introduced more recently, such as [START_REF] Flyer | Radial basis function-generated finite differences: A mesh-free method for computational geosciences[END_REF]. However they still have not been considered, up to our knowledge, for high-dimensional problems.

To tackle the curse of dimensionality, methods based on tensor networks or neural networks have been proposed [4,[START_REF] Khoromskij | Tensor numerical methods for multidimensional pdes: theoretical analysis and initial applications[END_REF][START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF]. Sparse approximation methods have also been developped [START_REF] Bungartz | Sparse grids[END_REF] and exploited to solve high-dimensional PDEs, see e.g. [START_REF] Shen | Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems[END_REF][START_REF] Shen | Efficient spectral sparse grid methods and applications to high-dimensional elliptic equations ii. unbounded domains[END_REF].

Moving to local approximation methods, we recall that it is possible to get pointwise evaluations of the solution using a probabilistic representation. The so-called Feynman-Kac representation mentioned earlier is such a probabilistic representation. There is no curse of dimensionality to tackle since we only focus on getting a few pointwise evaluations. This representation exists for several classes of linear and non-linear PDEs [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF] and yields an expression of the solution of the form

u(x) = E [Q(X x )] , ∀x ∈ O, ( 4
)
where Q is a functional depending on the data of the PDE and some stochastic diffusion process X x depending on the PDE and on the point at which the solution u is evaluated. Using a traditional Monte-Carlo estimate based on M i.i.d. samples of X x , we have access to noisy pointwise evaluations of u in different points x ∈ O.

A standard analysis shows that the precision of these pointwise evaluations generally behaves as O √ M -1 + ∆t α for some α > 0, where ∆t is a time-integration parameter for the approximation of the stochastic process [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. When we wish to approximate precisely u, e.g. up to machine precision, the computational costs are too high, even for today's computers.

Different methods allow to reduce the variance of the pointwise evaluations of the solution, among them we mention [START_REF] Giles | Multilevel monte carlo methods[END_REF][START_REF] Mascagni | -shell error analysis for "walk on spheres[END_REF]. The authors in [START_REF] Gobet | A spectral monte carlo method for the poisson equation[END_REF] introduced a variance reduction method with global control variates using the Feynman-Kac representation. In [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF], they described and analysed more precisely the use of these control variates. First they build a global approximation ũ of u in some finite-dimensional approximation space P. They use for example interpolation methods relying on estimations of pointwise evaluations of u on an unisolvent grid Γ P for P, obtained by Monte-Carlo methods for the probabilistic representation (4). Showing then that the error e = u -ũ admits a probabilistic representation, the authors build a global approximation ẽ of e based on estimations of pointwise evaluations obtained by Monte-Carlo methods. Updating ũ to ũ + ẽ, the process can be then iterated multiple times. It results in approximate pointwise evaluations of u on the grid Γ P with an error depending on the time integration error of the approximation error in P of the exact solution u [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF].

First contribution

When d 1, the approximation space P in which lives ũ and ẽ, and the associated grid Γ P must be well-chosen. Here we combine the control variate method from [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF] with the adaptive interpolation strategy from [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF] which builds adaptively and at the same time the approximation space P and the adapted unisolvent grid Γ P . We propose two adaptive versions of the algorithm of [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. Being given a global approximation ũk of u, our first algorithm builds at step k a sparse adaptive interpolation ẽk of the error u -ũk in the approximation space P k (and its associated grid Γ P k ). Our second contribution is a perturbed adaptive interpolation algorithm based on [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF], where at each step k of the algorithm, we compute an approximation ũ of u using the control variate algorithm from [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF], rather than the exact interpolation on P k . This work can be found in [START_REF] Billaud-Friess | Stochastic methods for solving high-dimensional partial differential equations[END_REF].

Parameter-Dependent PDEs and Discrete Optimization

We focus on approximating the map y ∈ Y → u(•; y) ∈ H, which is the solution of some parameter-dependent PDE, with H a Hilbert space of functions defined on O. Model order reduction (MOR) methods aim at reducing the complexity of complex numerical models. Among MOR we can find reduced basis methods (RBM) that are highly relevant for some classes of parameter-dependent PDEs and exploit a good approximability of the solution manifold M := {u(y) : y ∈ Y } by low-dimensional spaces [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF][START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic pde's[END_REF]. RBM build n-dimensional subspaces H n of H such that the distance between H n and M is close to the Kolmogorov n-width [START_REF] Pinkus | N-widths in Approximation Theory[END_REF] that characterizes the smallest distance between a n-dimensional subspace of H and M, that is the best approximation error that can be achieved by linear approximation tools. We focus particularly on greedy algorithms for RBM [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].

In an idealized framework, greedy algorithms adaptively build a nested sequence of reduced space (H n ) n≥1 according to the following procedure. At each step n ≥ 1, the current reduced space H n-1 is enriched with a snapshot u(y n ) which is selected among all the elements from M whose distance to H n-1 is maximal. This distance is in practice estimated using an a posteriori error estimate [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] inexpensive to compute. For parameter-dependent PDEs, a traditional error estimate is based on the norm of the residual [START_REF] Haasdonk | Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition[END_REF][START_REF] Haasdonk | Reduced basis methods for parametrized pdes-a tutorial introduction for stationary and instationary problems[END_REF]. Such error estimates usually deteriorate the performance of the greedy algorithm, that turns into a weak-greedy algorithm [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]. For some parameter-dependent problems (with low decay of the Kolmogorov n-width of their solution manifold) residual-based error estimates may yield a poor quality of the reduced space. Thus alternative techniques have been proposed, such as preconditioning [START_REF] Zahm | Interpolation of inverse operators for preconditioning parameter-dependent equations[END_REF] or residual-free error estimate [START_REF] Chen | A robust error estimator and a residualfree error indicator for reduced basis methods[END_REF].

However, in practical applications, maximizing the approximation error over a continuous set Y remains impossible. This is why discrete training sets Ξ ⊂ Y are usually considered and the selection of the snapshots is made among elements of the discretized solution manifold M Ξ := {u(ξ) : ξ ∈ Ξ}. The choice of these discrete training sets is particularly challenging for high-dimensional parameter spaces, i.e. when p 1. A random selection of these training sets seems an interesting technique when the solution manifold M can be well-approximated by low dimensional spaces [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]. This results in a probabilistic greedy algorithm.

The norm of the a posteriori error estimate in greedy algorithms can be written as an expectation using a Feynman-Kac representation of the error. In this context, in order to enrich the reduced space with a snapshot selected using a greedy procedure, we need to solve a discrete optimization problem of the following type max ξ∈Ξ E[Z(ξ)], (5) where, in the context of reduced basis methods, the (Z(ξ)) ξ∈Ξ are random variables depending on the reduced space. Such discrete optimization problems have already been widely studied since the 1950s [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF][START_REF] Lattimore | Bandit algorithms[END_REF] and were originally mentioned in medical applications. Bandit algorithms use Monte-Carlo estimates and concentration inequalities in order to provide a quasi-optimal solution of (5) with high probability. The result of such algorithms is a parameter ξ in Ξ that is probably approximately correct (PAC) [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF] in an absolute sense, which means

P E[Z(ξ )] -E Z( ξ) ≤ ≥ 1 -λ,
for some prescribed values of in (0, 1) and λ in (0, 1), with ξ ∈ arg max ξ∈Ξ E[Z(ξ)].

Second contribution

Existing bandit algorithms would need a careful estimation of in order to preserve guarantees of weak-greedy algorithms. To circumvent the estimation of , the second contribution of this thesis is a probably approximately correct (PAC) algorithm in relative precision which combines Monte-Carlo estimates with guaranteed relative precision [START_REF] Mnih | Empirical bernstein stopping[END_REF] with probably approximately correct (PAC) algorithm [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF]. It returns ξ in Ξ which satisfies

P E [Z(ξ )] -E Z( ξ) ≤ |E[Z(ξ )| ≥ 1 -λ, ( 6
)
for some prescribed values of in (0, 1) and λ in (0, 1), with ξ ∈ arg max ξ∈Ξ E[Z(ξ)].

Third contribution

The third contribution of this thesis combines reduced basis methods (RBM) with a snapshot selection based on the new PAC algorithm in relative precision. By reinterpretating the approximation error in a greedy algorithm as an expectation, we can select a snapshot to enrich the reduced basis using the PAC algorithm in relative precision. This results in a probabilistic greedy algorithm. Combining this algorithm with the use of random training sets [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], that are adapted to the case p 1, we propose a second probabilistic greedy algorithm that consists in randomly choosing a quasi-optimal snapshot among random training sets using the PAC algorithm in relative precision we proposed.

Thesis' Outline

This thesis is composed of five chapters.

First we focus on solving high-dimensional PDEs using probabilistic methods, which corresponds to our first contribution, published in [START_REF] Billaud-Friess | Stochastic methods for solving high-dimensional partial differential equations[END_REF]. A survey of probabilistic methods for linear parameter-independent PDEs is presented in Chapter 1 while the contribution is presented in Chapter 2.

In Chapter 3 we consider a discrete optimization problem for a function defined as an expectation, and introduces to our second contribution. This work has been submitted as [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF].

Finally we study parameter-dependent PDEs. In Chapter 4, we survey reduced basis methods for parameter-dependent PDEs. Chapter 5 introduces two probabilistic greedy algorithms for reduced basis methods based on the discrete optimization algorithm from Chapter 3.

Chapter 1 Probabilistic Methods for PDEs: a Survey

This chapter first introduces probabilistic representation for parameter-independent partial differential equations (PDEs) on bounded domains. Then numerical methods to simulate pointwise evaluations of the solution of parameter-independent PDEs using their probabilistic representation are surveyed.

Contents

T ] × D ⊂ R × R d .
The space variable is denoted by x = (x 1 , . . . , x d ). In the parabolic case, the time variable is denoted by t. We denote by L a linear differential operator taken under the non-divergence form Remark 1.1. The backward formulation in the parabolic case arises from problems in finance. In physics, forward problems are more common, where the terminal condition (1.3) is replaced by an initial condition. The references in the literature on theoretical results for PDEs usually treat the forward parabolic case. However, using the change of variable t → T -t, it is easy to derive theoretical results for the backward case. Moreover, both elliptic and parabolic problems can also be considered on unbounded domains. In order to simplify our study we will not consider PDEs on unbounded domains here and we refer the reader to the following non-exhaustive literature [START_REF] Comets | Calcul stochastique et modeles de diffusions[END_REF][START_REF] Friedman | Stochastic differential equations and applications[END_REF][START_REF] Karatzas | Graduate texts in mathematics[END_REF][START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF] for more details.

L(u)(t, x) = 1 2 d i,j=1 (σ(t, x)σ(t, x) T ) ij ∂ 2 x i x j u(t, x) + d i=1 b i (t, x)∂ x i u(t, x),
: [0, T ] × D → R ∂ t u + L(u) -ku + g = 0 in [0, T ) × D, (1.2) u(T, x) = f (T, x) for x ∈ D, (1.3) u(t, x) = f (t, x) for (t, x) ∈ [0, T ) × ∂D. (1.4) Equation (1.
No assumption were given yet on the different functions involved in the definition of the problems (1.2)-(1.4) and (1.5)-(1.6). For both parabolic and elliptic problems, uniqueness of the solution can be obtained using the so-called maximum-principle. We refer to [START_REF] Evans | Graduate studies in mathematics[END_REF]Chapter 2.3,Theorem 5] for the parabolic problem and to [36, Chapter 2.2, Theorem 5] for the elliptic problem. In these results, the uniqueness relates to classical solutions that are solutions in C([0, T ] × D) ∩ C 1,2 ([0, T ] × D), which means continuous over its domain of definition with continuous derivatives. It ensures that equation (1.2) is well-defined. Existence of classical solutions can be obtained too and we refer to [41, Chapter 6, Theorem 3.6] and references therein for the forward parabolic problem and to [41, Chapter 6, Theorem 2.4] and references therein for the elliptic problem.

Remark 1.2. A weaker notion of solution exists. It is related to the weak formulation of PDEs. We mention this aspect in Section 4.3.1

This chapter is a survey of probabilistic methods for both elliptic and parabolic problems defined above. We recall that such methods give access to pointwise evaluations of the solution u, i.e. to u(t, x) for chosen values of (t, x) in [0, T ] × D in the parabolic case, and to u(x) for chosen values of x in D in the elliptic case. Stochastic Differential Equations (SDEs) are at the core of the probabilistic methods introduced hereafter. We first introduce in Section 1.1 the notion of SDE. Then we show in Section 1.2 how to express the solution of a PDE as the expectation, depending on a diffusion process solution of a well-chosen SDE, related to data of the problem. This being done, we present numerical methods allowing to estimate pointwise evaluations of the solution u. Since the numerical methods are similar for both elliptic and parabolic problems, we present them in the parabolic case in Section 1.3 and adapt this presentation for the elliptic case in Section 2.2.

Stochastic Differential Equations

Stochastic calculus originates in the description of Brownian motion in the early 19 th century by Robert Brown, which characterizes a particle with erratic movements. These observations have led to the introduction of the notions of stochastic process and stochastic integrals that are at the core of the Itô calculus, introduced in the 1940s. Such notions give a rigorous framework for the definition of an SDE. We recall the most important ones here and refer the reader to [START_REF] Gobet | Introduction to stochastic calculus and to the resolution of pdes using monte carlo simulations[END_REF] and references therein for more details. First, we introduce a probability space (Ω, F, P), where Ω is a sample space, F a σ-algebra on Ω and P a probability function on F.

Definition 1.3 (Stochastic process with values in R d

). Let T be a subset of R. A stochastic process indexed by T is a family (X t ) t∈T of random variables defined on the probability space (Ω, F, P). The random variables X t may be valued in any measurable space (R d , G), where G is a σ-algebra on R d . Then we define a measurable (stochastic) process.

Definition 1.4 (Measurable and adapted process). A real-valued stochastic process

(X t ) t∈T is measurable if X : (t, ω) ∈ (T × Ω, B(T) ⊗ F) → X t (ω) ∈ (R d , B(R d )) is (B(T) ⊗ F, B(R d ))-measurable.
Recalling that a filtration is a family (F t ) t∈T ⊂ F such that

t 1 ≤ t 2 ⇒ F t 1 ⊂ F t 2 , (X t ) t∈T is said to be adapted to the filtration (F t ) t∈T if, for all t ≥ 0, X t : (Ω, F t ) → (R d , B(R d )) is (F t , B(R d ))-measurable.
Now we give a definition of the (standard) Brownian motion in dimension one, which is a particular stochastic process indexed by T = R + .

Definition 1.5 (One-dimensional Brownian motion). A (standard one-dimensional)

Brownian motion is a real-valued stochastic process (W t ) t∈R + , with continuous trajectories, such that 1. W 0 = 0 ; 2. any time increment W t -W s (0 ≤ s < t) has the distribution N (0, t -s) ; 3. for any 0 = t 0 < t 1 < . . . < t n , the increments {W t i+1 -W t i ; 0 ≤ i ≤ n -1} are independent. From this definition, we can deduce the main properties of the Brownian motion. It includes the symmetry property (if W is a Brownian motion, so is -W ). We refer the reader to [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]Chapter IV.1] for additional properties.

The extension to obtain a d-dimensional Brownian motion consists simply in considering W = (W 1 , . . . , W d ), where each W i is a one-dimensional Brownian motion, each one being independent from the others. In what follows, we will denote by (F t ) t∈R + the natural filtration of the Brownian motion. Unless otherwise specified, any adapted process will now be meant adapted to (F t ) t∈R + . We can then introduce the set of adapted and square-integrable processes indexed by T = [0, T ], which is defined by

H 2 T := H 2 [0,T ] = {φ adapted to (F t ) t∈[0,T ] such that φ 2 H 2 T := E T 0 |φ t | 2 dt < ∞}.
We also define the set of elementary processes, H elem , which contains the processes (φ t ) 0≤t≤T , such that φ t = φ t i for all t i < t ≤ t i+1 , where t 1 < t 2 < . . . is any time grid. We now introduce the notion of Stochastic Differential Equation (SDE) and conditions for existence and uniqueness of a solution to this SDE.

∀(t, x, y) ∈[s, S] × R d × R d |b(t, x) -b(t, y)| + |σ(t, x) -σ(t, y)| ≤ C b,σ |x -y|; (A1) sup s≤t≤S |b(t, 0)| + |σ(t, 0)| ≤ C b,σ . (A2)
For a given x s ∈ R d and for any s ≤ t ≤ S we consider the following SDE

X t = x s + t s b(r, X r )dr + t s σ(r, X r )dW r .
(

1.7)

There exists a unique process (X ) for some constant C depending on S -s and C b,σ . The last integral in (1.7) is defined in the sense of Ito. In particular for processes φ in H 2 T ∩ H elem we have

T 0 φ r dW r = t i ≤T φ t i W T ∧t i+1 -W t i . (1.8)
This definition can be extended for processes in H 2 T . We refer the reader to [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]Chapter IV.2] for this extension and we recommend [START_REF] Karatzas | Graduate texts in mathematics[END_REF] for a complete presentation on stochastic integration.

Feynman-Kac Representation

Our goal here is to present the links between SDEs and PDEs. To do so, we introduce the notion of Ito process and the Ito formula. Definition 1.7 (Ito process). Let (b t ) t≥0 : R + → R d and (σ t ) t≥0 = ((σ 1,t , . . . , σ d,t )) t≥0 : R + → R d×d be two adapted processes satisfying s 0 (|b r |+|σ r |) 2 dr < +∞ almost surely for all s ≥ 0. The Ito process (X t ) t≥0 , starting at x 0 , with drift coefficient b and diffusion coefficient σ is the d-dimensional stochastic process defined by

X t = x 0 + t 0 b r dr + t 0 σ r dW r . (1.9)
Remark 1.8. The process (σ t ) t≥0 can be taken more generally with values in R d×q , for some q ∈ N possibly different from d, the space dimension. We focus in this manuscript on the case q = d. Moreover, the definition of the Ito process has been taken using s = 0 as a time reference even if it is still possible to make it start at any time s, just like the stochastic process solution of (1.7). However, we stick to this time origin to avoid introducing too much notation. Theorem 1.9 (Ito formula). Let G be a real-valued function in C 1,2 (R + × R d ) and (X t ) t≥0 the d-dimensional Ito process defined by (1.9). Thus the process

(Y t = G(t, X t )) t≥0 is a 1-dimensional Ito process given by G(t, X t ) =G(0, x 0 ) + t 0 ∂ t G(s, X s )ds + t 0 ∇ x G(s, X s )b s ds + t 0 ∇ x G(s, X s )σ s dW s + 1 2 t 0 d k,l=1 ∂ 2 x k ,x l G(s, Xs)[σ s σ T s ] kl ds.
(1.10)

Ito formula allows to introduce probabilistic representations for PDEs, also known as Feynman-Kac representations. It originates in the 1940s, with the PhD work of Richard Feynman, then developed by Mark Kac, in particular in [START_REF] Kac | On distributions of certain wiener functionals[END_REF]. Such a representation is obtained by seeing the differential operator (1.1) as an infinitesimal generator of the d-dimensional diffusion process X t,x solution of the following SDE

X t,x s = X t,x t + s t b(r, X t,x r )dr + s t σ(r, X t,x r )dW r , X t,x t = x ∈ D, (1.11)
whose coefficients are defined using the coefficients of L. Here W is a d-dimensional Brownian motion. The stochastic process X t,x is uniquely defined, under suitable assumptions such as (A1) and (A2), by its value x at time t and the functions b and σ. Remark 1.10. In the elliptic case, the different functions involved are not timedependent and we take the time origin at t = 0 in (1.11). When there is no ambiguity, we simplify our notation with the convention X 0,x = X x . Applying now the Ito formula to a well-chosen function of X t,x , it is possible to obtain the Feynman-Kac representation of the solution of the considered problem (either (1.2)-(1.4) or (1.5)-(1.6)).

Backward Parabolic Boundary Value Problem

The following Feynman-Kac representation comes from [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]Theorem IV.4.5].

Theorem 1.11 (Feynman-Kac formula for backward parabolic PDE on bounded domain). Set T > 0 and assume that (P1) assumptions (A1) and (A2) are satisified,

(P2) D ⊂ R d is a bounded domain and we define τ t,x = inf{s > t : X t,x s / ∈ D} as the first exit time of X t,x from D, (P3) the boundary ∂D is regular for X t,x in the sense that ∀t ∈ [0, T ], ∀x ∈ ∂D, P(τ t,x = 0) = 1, (1.12) (P4) the functions f, g, k : [0, T ] × D → R are continuous, (P5) there exists a continuous function u : [0, T ] × D → R, of class C 1,2 on all open subset of [0, T [×D verifying (1.2)-(1.4).
Then u admits the following probabilistic representation

u(t, x) = E[G(u, X t,x )], (1.13)
where, setting τ t,x T = τ t,x ∧ T , we have

G(u, X t,x ) = f (τ t,x T , X t,x τ t,x T ) exp - τ t,x T t k(r, X t,x r )dr + τ t,x T t g(s, X t,x s ) exp - s t k(r, X t,x r )dr ds.
Assumption (P5) could be replaced by suitable assumptions on the coefficients of the PDE and the linear operator L, see e.g. [41, Chapter 6, Theorem 3.6], that ensure the existence of a classical solution to (1.2)- (1.4). An essential regularity assumption is needed on the boundary. In some Feynman-Kac representations (see e.g. [41, Chapter 6, Theorem 5.2]), a high regularity is required for ∂D in order to get rid of assumption (P5). However, in many practical examples (e.g. on polygonal domains) this regularity for the boundary is not verified. In order to avoid these assumptions, it is possible to refer to existence and uniqueness theorem [41, Chapter 6, Theorem 3.6]. The only needed assumption on the boundary ∂D in this theorem is the existence of a barrier at every point of (0, T ) × ∂D. This condition ensures a certain regularity of the boundary and thus the existence of a classical solution.

Other conditions exist, such as the exterior cone condition, also called Zaremba's cone condition, see [START_REF] Karatzas | Graduate texts in mathematics[END_REF]Chapter 4.2.C] for a discussion on regularity conditions of ∂D. Other works have been conducted linking these conditions [START_REF] Michael | Barriers for uniformly elliptic equations and the exterior cone condition[END_REF], but we will not develop them further in this manuscript.

Elliptic Boundary Value Problem

For the elliptic problem, we recall the following notation X x = X 0,x (see Remark 1.10). The following result is taken from [START_REF] Comets | Calcul stochastique et modeles de diffusions[END_REF]Theorem 2.4].

Theorem 1.12 (Feynman-Kac formula for elliptic PDE on bounded domain). Assume that

(E1) D is an open connected bounded domain of R d , regular in the sense that, if τ x = inf {s > 0 : X x s /
∈ D} is the first exit time of D for the process X x , solution of (1.11), we have

P(τ x = 0) = 1, x ∈ ∂D, (E2) b, σ are Lipschitz functions over D, (E3) f is continuous on ∂D, g and k ≥ 0 are Hölder-continuous functions on D, (E4) (uniform ellipticity assumption) there exists c > 0 such that d i,j=1 σ(x)σ(x) T ij Υ i Υ j ≥ c d i=1 Υ 2 i , Υ ∈ R d , x ∈ D.
Then, there exists a unique solution of

(1.5) in C D ∩ C 2 (D), which satisfies for all x ∈ D u(x) = E [F (u, X x )] (1.14)
where

F (u, X x ) = u(X x τ x ) exp - τ x 0 k(X x t )dt + τ x 0 A(u)(X x t ) exp - t 0 k(X x s )ds dt, with u(X x τ x ) = f (X x τ x ) and A(u)(X x t ) = g(X x t ).
Proof. We refer the reader to [START_REF] Friedman | Stochastic differential equations and applications[END_REF] and references therein for a proof of the existence of a classical solution.

Using now the

Y s = Y 0 + s 0 b Y r dr + s 0 σ Y r dW r Z s = e s 0 crdr
we obtain

Y s Z s = Y 0 + s 0 e r 0 cs 1 ds 1 (c r Y r + b Y r )dr + s 0 e r 0 cs 1 ds 1 σ Y r dW r
Taking Y s = u(X x s ) with X x s the solution process, starting from x at time 0, of the SDE (1.11) whose infinitesimal generator is L (whose drift coefficient is b and diffusion coefficient is σ) and c s = -k(X x s )1 s 0 and using the fact that u is the solution of (1.5)-(1.6), we get

u(x) =u(X x s )e -s 0 k(X x r )dr + s 0 e -r 0 k(X x s 1 )ds 1 g(X x r )dr - s 0 e -r 0 k(X x s 1 )ds 1 ∇ x u(X x r )σ(X x r )dW r .
Choosing the stopping time as an exit time

τ x n = inf s > 0 : d(X x s , ∂D) 1 n ,
we have for s = τ x n , taking at the same time the expectation in the previous formula (the second integral has a zero-valued expectation)

u(x) = E u(X x τ x n )e -τ x n 0 k(X x r )dr + τ x n 0 e -r 0 k(X x s )ds g(X x r )dr
With τ x n → τ x we take the limit and we get (1.14) by continuity on D of u, g and k.

Compared to Theorem 1.11, Theorem 1.12 does not assume the existence of a classical solution: this existence is already ensured by the assumptions of the theorem. Among the other assumptions, (E2) ensure the uniqueness of the solution to the SDE, while (E4) ensures that the exit time τ x is almost surely finite for every x in D (see [START_REF] Freidlin | Functional integration and partial differential equations[END_REF]Lemma 3.3.1]).

Numerical Methods

In this section, we first detail the numerical tools we will use to estimate pointwise evaluations of the solution of (1.2)-(1.4) using the Feynman-Kac formula (1.13). It includes time-integration schemes in order to simulate the diffusion process X t,x and Monte-Carlo methods for expectation estimation. Finally, after a brief analysis of the error resulting from the numerical tools we use. We present in Section 1.3.3 a complexity analysis of these numerical tools together with a survey on techniques reducing this complexity.

Approximation of Diffusion Processes

The process X t,x is traditionally approximated by a time-integration scheme S designed according to the SDE satisfied by X t,x . Letting t n = t + n∆t, n ∈ N, X t,x is approximated by a piecewise constant process X t,x,∆t , where

X t,x,∆t s = X t,x,∆t t n for s ∈ [t n , t n+1 [ and X t,x,∆t t n+1 = S(X t,x,∆t t n , b, σ, t n , ∆W n , ∆t). (1.15)
Here ∆W n = W t n+1 -W t n is an increment of the standard Brownian motion. In particular the Euler-Maruyama (see Section 2.2.1) time-integration schemes is the most basic way to simulate trajectories of X t,x,∆t . Remark 1.13. Piecewise constant is not mandatory in the definition of X t,x,∆t . For example in [5] the authors use piecewise affine trajectories (which are continuous). In our analysis we chose piecewise constant approximations for a sake of simplicity. Using a time-integration scheme, we can also estimate the exit time τ t,x by

τ t,x ≈ τ t,x,∆t := min t n > 0 : X t,x,∆t t n / ∈ D . (1.16)
To quantify the accuracy of the stochastic process X t,x,∆t (one can show that X t,x,∆t is a Ito process, see e.g. [47, Theorem V.1.2]) as an approximation of X t,x defined by (1.11), we will introduce in Section 1.3.3 the notion of time-integration error.

Monte-Carlo Methods

We can use Monte-Carlo approach to approximate the time-integrated solution

u ∆t (t, x) = E G(u, X t,x,∆t ) .
(1.17) Thus, denoting by (X t,x,∆t (ω m )) m=1...M a family of M independent trajectories of the stochastic process X t,x,∆t (simulated using time-integration schemes introduced in section 1.3.1) we have

u ∆t (t, x) ≈ u ∆t,M (t, x) = 1 M M m=1 G(u, X t,x,∆t (ω m )), (1.18)
where G(u, X t,x,∆t (ω m )) can be computed by the rectangle quadrature formula with left point (which is exact here since X t,x,∆t is a piecewise constant stochastic process),

G(u, X t,x,∆t (ω m )) =f (t L(ωm) , X t,x,∆t t L(ωm) (ω m )) L(ωm)-1 l=0 exp -k(t l , X t,x,∆t t l )∆t + L(ωm)-1 l=0 g(t l , X t,x,∆t t l ) l j=0 exp -k(t j , X t,x,∆t t j )∆t , (1.19)
with L a random variable such that t L = τ t,x,∆t ∧ T .

Remark 1.14. Using the rectangle quadrature formula with right point is possible but when t L (ω m ) = τ t,x,∆t (ω m ), some quantities in (1.19) are not well-defined, since X t,x,∆ t L (ωm) is almost surely not in D, the space domain where f, g, k are defined.

Statistical Error for the Simulation of Pointwise Evaluations

The estimation error ∆t,M is defined as

∆t,M (t, x) := u(t, x) -u ∆t,M (t, x). (1.20)
It naturally depends on two parameters, the time-integration discretization parameter ∆t and the number of Monte-Carlo samples M . The error can be decomposed in two terms

∆t,M (t, x) = E [ ∆t,M (t, x)] + ∆t,M (t, x) -E [ ∆t,M (t, x)] . (1.21) We notice that E [ ∆t,M (t, x)] := u(t, x) -u ∆t (t, x)
. This error will be referred to as the time-integration error and be denoted ∆t (t, x). The remaining term is due to the Monte-Carlo method and will be noted M C (t, x). We notice that it corresponds to u ∆t (t, x)-u ∆t,M (t, x) and that its variance is equal to the variance of the estimate. On one hand, the time-integration error ∆t is due to the time-integration scheme (1.15) we use. To quantify the precision of the time-integration scheme, we study its strong convergence.

Definition 1.15 (Strong convergence of time-integrated stochastic process).

A discrete time approximation X t,x,∆t converges strongly with order γ, on [t, T ], towards X t,x when there exists a positive constant C, independent of ∆t, such that, for all 0 < ∆t ≤ ∆t lim , we have

E sup t≤s≤T X t,x s -X t,x,∆t s ≤ C∆t γ . (1.22)
Under assumptions slightly stronger than (A1) and (A2) on the coefficients b and σ of (1.11), it has been proven in [47, Theorem V.2.1] that the Euler-Maruyama scheme (presented in Section 1.2.2 in the stationary case) has a strong order of convergence γ = 1 2 . Other numerical schemes, assuming stronger conditions on the coefficients of the SDE defining X t,x , have been proposed to improve the order of convergence. We refer to [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] for explicit time-integration schemes of higher order. Remark 1.16. An interesting point is that, when the coefficients b and σ involved in the definition of the infinitesimal generator of X t,x , see (1.11), are constant, the Euler-Maruyama scheme does not differ from these high-order schemes: the corrective terms added to increase the order of convergence depend on the derivatives of b and σ. In this situation, other techniques exist to improve the convergence order. We will present them later. However, the strong convergence order does not allow to study the behaviour of ∆t directly. In particular, assumptions on the functional involved in the Feynman-Kac representation have to be made. For example, recalling and assuming G(u, •) is Lipschitz, we can show that | ∆t (t, x)| ∆t 1/2 when we use the Euler-Maruyama time-integration scheme. However this is a really naive way to study ∆t . A more advanced way is often referred to as the weak error study and can be made independently from the strong error.

On the other hand, the Monte-Carlo error M C can be controlled in probability non-asymptotically using concentration inequalities (see [START_REF] Frikha | Concentration bounds for stochastic approximations[END_REF] in the parabolic case on unbounded domain where G(u, X t,x,∆t ) = f (T, X t,x,∆t T ) for a deterministic value of T ) which means that for all > 0, it exists C, not depending on or M , such that

P (|u ∆t (t, x) -u ∆t,M (t, x)| > ) ≤ 2 exp -M 2 C .
From this, we can deduce confidence intervals with length behaving in M -1/2 as expected for Monte-Carlo methods.

In light of these results and the error decomposition (1.21), a balance has to be found between a fine time-discretization parameter ∆t and a high number of Monte-Carlo samples M . This balance depends naturally on the order of convergence of ∆t with respect to ∆t (equal to 1 2 here) and also on the order of convergence of M C with 1/M (also equal to 1 2 here). If we denote by the target accuracy we aim at reaching for our pointwise evaluations, we have to take M = O( -2 ) and ∆t = O( 2), due to the order of convergence of the time-integration error and the Monte-Carlo error that are in both cases equal to 1 2 . The complexity being O(M ∆t -1 ), as function of the wanted precision it behaves like O( -4 ). Thus, for very precise pointwise evaluations, the complexity becomes too high. A solution to overcome these accuracy limitations is to improve the order of convergence of both the time-integration technique and the Monte-Carlo method. We recall hereafter a few ways to enhance these orders. In Section 1.3.3.1 we survey high-order time-integration methods while Section 1.3.3.2 is devoted to enhanced Monte-Carlo methods.

Time-Integration Error ∆t

It has been proven in [47, Theorem V.3.1] that for a parabolic Cauchy problem on unbounded domains, i.e. when no exit-time has to be estimated, the timeintegration error with an Euler-Maruyama scheme behaved as O(∆t). Combined with classical Monte-Carlo methods this would allow an overall complexity in O( -3 ). However, the assumptions on the coefficients of the PDE were very strict to obtain this result (C 1,4 ([0, T ] × D), plus the existence of a classical solution with the same regularity). For a simpler but classical study case (parabolic on unbounded domain with g ≡ k ≡ 0), the time-integration error E

[f (X t,x T )] -E[f (X t,x,∆t T
)] has been studied in [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]5] while using the Euler-Maruyama time-integration scheme. In [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], taking f and the coefficients of the SDE driving X t,x to be smooth functions, the authors give an expansion of the error in term of powers of ∆t

E[f (X t,x T )] -E[f (X t,x,∆t T )] = R r=1 a r ∆t r + O(∆t R+1 ), (1.23)
where the coefficients (a r ) R+1 r=1 do not depend on ∆t. In [5], a similar expansion is obtained when f is only measurable and under non-degeneracy assumption on the infinitesimal generator of X t,x . Such expansions are then useful to propose new approximations (of higher order) for E[f (X t,x T )] by the principle of Romberg extrapolation, that we present here when R = 1 in (1.23). In this case we can notice that

E[f (X t,x T )] -2E[f (X t,x,∆t/2 T )] -E[f (X t,x,∆t T )] = O(∆t 2 ), (1.24)
and then we have a time-integrated approximation whose order of convergence is 2. This technique can be extended to any order R + 1 if we have an expansion of order R. For more generalities on the subject, see for example [99, Section 2.3] or [START_REF] Pagès | Multi-step richardson-romberg extrapolation: remarks on variance control and complexity[END_REF].

Remark 1.17. Such a technique can be used to propose high-order approximation in the case where the coefficients b and σ are constant (see Remark 1.16). These techniques, even if interesting in practice, cannot be used for any function f (for which the existence of an expansion might be harder to prove), and we still have to use time-integration schemes with time-integration parameter lower than ∆t (∆t/2 for example when using the Romberg extrapolation technique to get an order 2 approximation), which is numerically costly.

Another major concern, not treated by the already presented literature, is the exittime estimation which introduces also a bias: we easily notice that τ t,x,∆t ≥ τ t,x which introduces a bias in O( √ ∆t) according to [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF]Theorem 4]. To deal with this overestimation it is proposed in [50, Theorem 5] to consider a slightly modified exit-time to regain a time-integration error in o( √ ∆t): this exit time is

τ t,x,∆t 1 := inf{s > t : X t,x,∆t / ∈ D ∆t },
where

D ∆t = {x ∈ D : d(x, ∂D) > C √ ∆t} is a smaller domain of D, see Figure 1.1.
This technique is very useful in practice since it is easy to implement. However the resulting convergence has only been proved to behave like o(∆t 1/2 ), the exact order has not been explicited yet.

In the general parabolic case (1.2)-(1.4), high order of convergence are difficult to reach using traditional techniques without very restrictive assumptions on the PDE. This is due to the complicated form of the Feynman-Kac functional in (1.13).

However, for easier study cases, interesting alternatives can be proposed. In fact one can notice that the simulation of stochastic processes' trajectories is complicated. If we were considering problems for which the simulation of a whole trajectory is not needed, we could investigate better precision without higher numerical complexity.

Here we consider two different elliptic problems (1.5)-(1.6) for which we do not need the simulation of the entire trajectory. A first one where k ≡ b ≡ g ≡ 0, f ≡ σ ≡ 1.

In this simplified case, known as the Poisson problem, the Feynman-Kac formula becomes

u(x) = E[f (X x τ x )], (1.25)
Figure 1.1: Comparison between the estimation of the exit time τ t,x,∆t of D and the estimation of the exit time τ t,x,∆t 1 of D ∆t [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF], illustrated for a particular trajectory of X t,x,∆t . The domain D ∆t (delimited by ---) allows to estimate an exit position (and exit time) with a bias lower than the one using D (delimited by -).

where τ x = τ 0,x by convention. This problem consist here in estimating the exit point of the domain D. A second one where k ≡ b ≡ 0, f ≡ g ≡ σ ≡ 1. The Feynman-Kac formula thus becomes

u(x) = E[f (X x τ x ) + τ x ]. (1.26)
Here we have to estimate both the exit point and the exit time of the domain D. [START_REF] Hwang | Analysis and comparison of green's function first-passage algorithms with "walk on spheres[END_REF]. Both procedures work the same, until the selection of x 5 . When x 4 has been drawn, it lays at a distance to the boundary smaller than δ GF F P . So we can use the GFFP method and draw the next point on the intersection of a sphere and the boundary (in bold on the left). The WOS draws the next sample still on a sphere (in bold on the right) and stops since the next sample x 5 lays at a distance lower than δ W OS of the boundary.

A groundbreaking technique has been theoretically introduced in [START_REF] Muller | Some continuous monte carlo methods for the dirichlet problem[END_REF] to solve (1.25): the Walk-On-Sphere (WOS) procedure. It is based on the observation that the exit position of a sphere, for a Brownian motion starting in the center of this sphere, is uniformly distributed over the sphere. Then, consider procedure described here 1. Start in x 0 = x and set k = 0.

2. Find the largest sphere, centered in x k , contained in D. Draw according to the uniform distribution a point x k+1 on this sphere. If this point is on ∂D, stop.

3. Otherwise, increment k and repeat the procedure from step 2.

It returns a point on ∂D and this point has the same distribution as X x τ x . However, this procedure almost surely does not stop. Simply modifying it by stopping when the point x k+1 is closer than some parameter δ W OS from the boundary, [START_REF] Motoo | Some evaluations for continuous monte carlo method by using brownian hitting process[END_REF] analyzes the numerical error resulting from this technique. It gives then a result on the average number N (δ W OS ) of spheres before stopping: it is bounded, up to a constant, by | log(δ W OS )|. Under regularity assumptions on the boundary, the error of such method is in O(δ W OS ) [START_REF] Mascagni | -shell error analysis for "walk on spheres[END_REF]. Thus, the complexity (to get one sample) as a function of the precision behaves in O(log( -1 )) which is better than techniques whose complexity is a power of -1 . The WOS procedure has then been adapted to rectangular domains in [START_REF] Deaconu | A random walk on rectangles algorithm[END_REF]. It has been later improved under the name of Green Function First Passage (GFFP) [START_REF] Given | A first-passage algorithm for the hydrodynamic friction and diffusion-limited reaction rate of macromolecules[END_REF] in order to avoid the absorption layer of size δ W OS and then to simulate exactly τ x . An illustrated comparison of GFFP and WOS methods can be found on Figure 1.2, see [START_REF] Hwang | Analysis and comparison of green's function first-passage algorithms with "walk on spheres[END_REF] for a comparison in term of complexity of the two methods. This GFFP method thus requires the exact simulation of the exit time over domains intersected with spheres (see Figure 1.2). An improvement of GFFP, called the Simulation-Tabulation method has even been proposed in [START_REF] Hwang | The simulation-tabulation method for classical diffusion monte carlo[END_REF]. These techniques are however difficult to extend to non-classical domain D (or to high dimensional problems on complicated domain) since it is necessary to get samples according to the distribution of the exit time of a Brownian motion from domain more complicated than spheres (see [START_REF] Given | A first-passage algorithm for the hydrodynamic friction and diffusion-limited reaction rate of macromolecules[END_REF]Appendix] among other references). Moreover, it is only available for the Poisson equation, when the diffusion process is in fact the Brownian motion.

All these techniques to estimate the exit position from a bounded domain of the Brownian motion have been completed more recently by the study of the exit time of other one-dimensional stochastic processes [START_REF] Herrmann | The first-passage time of the brownian motion to a curved boundary: an algorithmic approach[END_REF][START_REF] Herrmann | Exact simulation of first exit times for onedimensional diffusion processes[END_REF][START_REF] Herrmann | Exact simulation of diffusion first exit times: algorithm acceleration[END_REF]. This new class of algorithms, made to estimate (1.26), are called the Walk-On-Moving-Spheres (WOMS) algorithms [START_REF] Deaconu | The walk on moving spheres: a new tool for simulating brownian motion's exit time from a domain[END_REF].

Another field of interest consists in studying PDEs with non-regular coefficients, which differs from the classical presentation adopted in Theorem 1.11 or Theorem 1.12. We refer the reader to [START_REF] Bossy | Probabilistic interpretation and random walk on spheres algorithms for the poisson-boltzmann equation in molecular dynamics[END_REF] and references therein for an overview for estimating (1.26) with such coefficients.

Promising numerical results allow to exactly simulate trajectories of diffusion processes in the case d = 1 [START_REF] Beskos | Exact simulation of diffusions[END_REF][START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF][START_REF] Etore | Exact simulation for solutions of one-dimensional stochastic differential equations with discontinuous drift[END_REF]. They are presented as allowing to reach higher precision for the same numerical complexity. However, up to our knowledge, such methods have not yet been employed for solving PDEs.

As a conclusion for this part, we summarize all the literature presented above in Table 1 O( -γ ), γ ∈ (0, 1] depending on the order of the expansion (1.25), small d WOS and GFFP [START_REF] Hwang | Analysis and comparison of green's function first-passage algorithms with "walk on spheres[END_REF] High convergence order but difficult to extend to complex or high-dimensional space domain or to other diffusion processes O(log( -1 ))

(1.26), small d WOMS [START_REF] Deaconu | The walk on moving spheres: a new tool for simulating brownian motion's exit time from a domain[END_REF] or [START_REF] Herrmann | Exact simulation of first exit times for onedimensional diffusion processes[END_REF][START_REF] Herrmann | Exact simulation of diffusion first exit times: algorithm acceleration[END_REF] High convergence order but difficult to extend to more complex or high-dimensional space domain or to other diffusion processes [START_REF] Étoré | Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process[END_REF][START_REF] Etore | Exact simulation for solutions of one-dimensional stochastic differential equations with discontinuous drift[END_REF] Exact simulation of trajectories but the diffusion coefficient has to be equal to 1

O(log( -1 )) (X t,x s ) 0≤s≤T , d = 1
Lower complexity at given precision compared to standard simulation of trajectories.

Table 1.1: Summary of the advantages and drawbacks of each numerical timeintegration technique. The complexity column corresponds to the complexity to reach a precision of for ∆t .

Monte-Carlo Error M C

Regarding the Monte-Carlo error M C , the order of convergence of the error, with respect to the number of Monte-Carlo samples M , behaves in O(M -1/2 ). This order is difficult to reduce, except using multilevel Monte-Carlo methods that we will present later. However it is possible to improve the estimate u ∆t,M (t, x) by using variance reduction techniques. In order to see the interest of these techniques, we introduce

V [u ∆t,M (t, x)] = V[G(u, X t,x,∆t )] M . (1.27)
When the variance V[G(u, X t,x,∆t )] is large, the estimation using the classic Monte-Carlo estimate is not adapted for precise estimation of u(t, x). It is thus possible to use variance reduction techniques, e.g. antithetic sampling. Suppose that X t,x,∆t and Xt,x,∆t have the same law. In particular since their expectations are equal, we have

u ∆t (t, x) = E[G(u, X t,x,∆t )] = E[G(u, Xt,x,∆t )] = E[G(u, X t,x,∆t ) + G(u, Xt,x,∆t )] 2
Using the last equality, it is possible to propose the following estimate

ũ∆t,M (t, x) = 2 M M/2 m=1 G(u, X t,x,∆t (ω m )) + G(u, Xt,x,∆t (ω m )) 2 , (1.28)
where the X t,x,∆t (ω m ) are i.i.d. samples of X t,x,∆t and the Xt,x,∆t (ω m ) are i.i.d. samples of Xt,x,∆t . This estimate is better than u ∆t,M (t, x) when we can prove that the covariance Cov G(u, X t,x,∆t ), G(u, Xt,x,∆t ) < 0. This variance reduction technique is particularly relevant and easy to implement in some cases here since when X t,x,∆ is the standard discretized Brownian motion, then so is -X t,x,∆ . Other variance reduction techniques exist and we refer the reader to [47, Chapter III] for a complete presentation.

We focus now on multilevel methods [START_REF] Giles | Multilevel monte carlo methods[END_REF]. In some cases, these techniques also allow to improve the convergence rate of M C with reference to M . If we want to estimate E[Z L ] but Z L is costly to sample, it might me better, if we have access to an approximation Z L-1 of Z L which is less costly to sample but also less precise, to use the following equality

E[Z L ] = E[Z L-1 ] + E[Z L -Z L-1 ],
in order to propose the following estimate

Ê[Z L ] = M 1 m=1 Z L-1 (ω 1 m ) M 1 + M 2 m=1 Z L (ω 2 m ) -Z L-1 (ω 2 m ) M 2 , (1.29)
where the Z j (ω i ) are independent samples drawn according to the law of Z j . It is important to make sure that the samples used in the difference Z L (ω 2 m ) -Z L-1 (ω 2 m ) are drawn according the same underlying elementary event ω 2 m . This multilevel Monte-Carlo method with two levels can then be generalized to a L level multilevel Monte-Carlo method using the following equality

E[Z L ] = L-1 l=0 E[Z L-l -Z L-l-1 ], Z 0 = 0, (1.30)
just assuming that we have access to a sequence of approximation with decreasing cost and precision (Z L-l ) L l=1 . Multilevel methods are particularly relevant for the simulation of Feynman-Kac representation, e.g. [START_REF] Higham | Mean exit times and the multilevel monte carlo method[END_REF], since it is easy to have access to the different levels required to use them: each level l corresponds to a Euler-Maruyama time-integration scheme of parameter ∆t l , with ∆t l that decreases with l. Of course an optimization has to be made for the choice of the number of samples M l to be used at each level l, as well as the choice of the different Z l , i.e. of ∆t l , when this choice can be driven. Improved non-asymptotic confidence intervals have been proposed for the multilevel method [START_REF] Jourdain | Non-asymptotic error bounds for the multilevel monte carlo euler method applied to sdes with constant diffusion coefficient[END_REF]. The Monte-Carlo error M C using these multilevel methods behaves at most and under additional assumptions, as a function of the reached precision , in O( log( -1 )) (see among others [47, Theorem IV.3.1]).

Chapter 2 Probabilistic Methods for Solving High-Dimensional Elliptic PDEs

This chapter introduced our first contribution. It consists in an algorithm to approximate up to any precision the solution of a high-dimensional elliptic partial differential equation (PDE) on a bounded domain. As we have seen in Section 1.3.3, the complexity as a function of the reached precision behaves polynomially in -1 . Taking a closer look at this complexity, reaching low precision, e.g. machine precision, through these methods would cost a high number of Monte-Carlo samples together with time-integration schemes with very fine parameter ∆t. While it may be possible to conduct such techniques for a single pointwise evaluation, such an approach are not relevant in a multi-query context, which is when several pointwise evaluations are required. Particularly since deterministic methods are relevant in this context for low-dimensional problems, i.e. when d is small.

Contents

Moreover the most efficient time-integration techniques presented in Section 1.3.3 mainly for low-dimensional problems are difficult to implement or to extend to highdimensional problems. Thus it is particularly costly to reach machine precision for pointwise evaluations of high-dimensional problems.

In this chapter, we present work that has already been published in [START_REF] Billaud-Friess | Stochastic methods for solving high-dimensional partial differential equations[END_REF]. We consider high-dimensional PDEs, i.e. d 1, and in the multi-query context. We propose two algorithms combining global approximation techniques with the probabilistic methods for PDEs seen in Chapter 1. These algorithms provide heuristics allowing to reach machine precision while approximating the solution of a high-dimensional PDE on bounded domain.

Title: Probabilistic methods for solving high-dimensional PDEs. Authors: Marie Billaud-Friess ♣ , Arthur Macherey ♣,♠ , Anthony Nouy ♣ and Clémentine Prieur ♠ .

♣ Centrale Nantes, LMJL, UMR CNRS 6629, 1 rue de la Noë, 44321 Nantes ♠ Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP*, LJK, 38000 Grenoble, France *Institute of Engineering Univ. Grenoble Alpes Abstract: We propose algorithms for solving high-dimensional partial differential equations (PDEs) that combine a probabilistic interpretation of PDEs, through Feynman-Kac representation, with sparse interpolation. Monte-Carlo methods and time-integration schemes are used to estimate pointwise evaluations of the solution of a PDE. We use a sequential control variates algorithm, where control variates are constructed based on successive approximations of the solution of the PDE. Two different algorithms are proposed, combining in different ways the sequential control variates algorithm and adaptive sparse interpolation. Numerical examples will illustrate the behavior of these algorithms.

Introduction

We consider the solution of an elliptic partial differential equation We are interested in approximating the solution of (2.1) up to a given precision. For high dimensional PDEs (d 1), this requires suitable approximation formats such as sparse tensors [START_REF] Bungartz | Sparse grids[END_REF][START_REF] Shen | Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems[END_REF] or low-rank tensors [START_REF] Oseledets | Tensor-train decomposition[END_REF][START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF][START_REF] Hackbusch | Numerical tensor calculus[END_REF]4,[START_REF] Nouy | Low-rank methods for high-dimensional approximation and model order reduction. Model reduction and approximation[END_REF]. Also, this requires algorithms that provide approximations in a given approximation format. Approximations are typically provided by Galerkin projections using variational formulations of PDEs. Another path consists in using a probabilistic representation of the solution u through Feynman-Kac formula, and Monte-Carlo methods to provide estimations of pointwise evaluations of u (see e.g., [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]). This allows to compute approximations in a given approximation format through classical interpolation or regression [START_REF] Beck | Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations[END_REF][START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward Bibliography stochastic differential equations[END_REF][START_REF] Beck | Solving stochastic differential equations and kolmogorov equations by means of deep learning[END_REF]. 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[END_REF], the authors consider interpolations on fixed polynomial spaces and propose a sequential control variates method for improving the performance of Monte-Carlo estimation. In this paper, we propose algorithms that combine this variance reduction method with adaptive sparse interpolation [START_REF] Chkifa | Sparse adaptive taylor approximation algorithms for parametric and stochastic elliptic pdes[END_REF][START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF].

A(u) = g in D, u = f on ∂D, ( 2
The outline is as follows. In Section 2.2, we recall the theoretical and numerical aspects associated to probabilistic tools for estimating the solution of (2.1). We also present the sequential control variates algorithm introduced 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[END_REF]. In Section 2.3 we introduce sparse polynomial interpolation methods and present a classical adaptive algorithm. In Section 2.4, we present two algorithms combining the sequential control variates algorithm from Section 2.2 and adaptive sparse polynomial interpolation. Finally, numerical results are presented in Section 2.4.

Probabilistic Tools for Solving PDEs

We consider the problem (2.1) with a linear partial differential operator defined by

A(u) = -L(u) + ku,
where k is a real valued function defined on D, and where

(2.2) L(u)(x) = 1 2 d i,j=1 (σ(x)σ(x) T ) ij ∂ 2 x i x j u(x) + d i=1 b i (x)∂ x i u(x)
is the infinitesimal generator associated to the d-dimensional diffusion process X x solution of the stochastic differential equation

(2.3) dX x t = b(X x t )dt + σ(X x t )dW t , X x 0 = x ∈ D,

Pointwise Evaluations of the Solution

The following theorem recalls the Feynman-Kac formula (see [START_REF] Comets | Calcul stochastique et modeles de diffusions[END_REF]Theorem 2.4] or [START_REF] Friedman | Stochastic differential equations and applications[END_REF]Theorem 2.4] and the references therein) that provides a probabilistic representation of u(x), the solution of (2.1) evaluated at x ∈ D.

σ(x)σ(x) T ij Υ i Υ j ≥ c d i=1 Υ 2 i , Υ ∈ R d , x ∈ D.
Then, there exists a unique solution of (2.1) in C D ∩ C 2 (D), which satisfies for all x ∈ D

u(x) = E [F (u, X x )] (2.4)
where

F (u, X) = u(X x τ x ) exp - τ x 0 k(X x t )dt + τ x 0 A(u)(X x t ) exp - t 0 k(X x s )ds dt, with u(X x τ x ) = f (X x τ x ) and A(u)(X x t ) = g(X x t
). Note that F (u, X x ) in (2.4) only depends on the values of u on ∂D and A(u) on D, which are the given data f and g respectively. A Monte-Carlo method can then be used to estimate u(x) using (2.4), which relies on the simulation of independent samples of an approximation of the stochastic process X x . This process is here approximated by an Euler-Maruyama scheme. More precisely, letting t n = n∆t, n ∈ N, X x is approximated by a piecewise constant process X x,∆t , where

X x,∆t t = X x,∆t t n for t ∈ [t n , t n+1 [ and X x,∆t t n+1 = X x,∆t t n + ∆t b(X x,∆t t n ) + σ(X x,∆t t n ) ∆W n , X x,∆t 0 = x. (2.5)
Here ∆W n = W t n+1 -W t n is an increment of the standard Brownian motion. For details on time-integration schemes, the reader can refer to [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Letting {X x,∆t (ω m )} M m=1 be independent samples of X x,∆t , we obtain an estimation u ∆t,M (x) of u(x) defined as

u ∆t,M (x) := 1 M M m=1 F u, X x,∆t (ω m ) = 1 M M m=1 f (X x,∆t τ x,∆t (ω m )) exp - τ x,∆t 0 k(X x,∆t t (ω m ))dt + τ x,∆t 0 g(X x,∆t t (ω m )) exp - t 0 k(X x,∆t s (ω m ))ds dt (2.6)
where τ x,∆t is the first exit time of D for the process X x,∆t (ω m ), given by

τ x,∆t = inf t > 0 : X x,∆t t / ∈ D = min t n > 0 : X x,∆t tn / ∈ D .
Remark 2.2. In practice, f has to be defined over R d and not only on the boundary ∂D. Indeed, although X x τ x ∈ ∂D with probability one, X x,∆t τ x,∆t ∈ R d \D with probability one. The error can be decomposed in two terms

u(x) -u ∆t,M (x) = ε ∆t u(x) -E F u, X x,∆t + E F u, X x,∆t - 1 M M m=1 F u, X x,∆t (ω m ) ε M C , (2.7)
where ε ∆t is the time integration error and ε M C is the Monte-Carlo estimation error. Before discussing the contribution of each of both terms to the error, let us introduce the following additional assumption, which ensures that D does not have singular points [1] .

(E5) Each point of ∂D satisfies the exterior cone condition which means that, for all x ∈ ∂D, there exists a finite right circular cone K, with vertex x, such that K ∩ D = {x}. Under assumptions (E1)-(E5), it can be proven [49, §4.1] that the time integration error ε ∆t converges to zero. It can be improved to O(∆t 1/2 ) by adding differentiability assumptions on the boundary [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF]. The estimation error ε M C is a random variable with zero mean and standard deviation converging as O(M -1/2 ). The computational complexity for computing a pointwise evaluation of u ∆t,M (x) is in O (M ∆t -1 ) in expectation for ∆t sufficiently small [2] , so that the computational complexity for achieving a precision (root mean squared error) behaves as O( -4 ). This does not allow to obtain a very high accuracy in a reasonable computational time. The convergence with ∆t can be improved to O(∆t) by suitable boundary corrections [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF], therefore yielding a convergence in O( -3 ). To further improve the convergence, high-order integration schemes could be considered (see [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] for a survey). Also, variance reduction methods can be used to further improve the convergence, such as antithetic variables, importance sampling, control variates (see [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to nonlinear[END_REF]). Multilevel Monte-Carlo [START_REF] Giles | Multilevel monte carlo methods[END_REF] can be considered as a variance reduction method using several control variates (associated with processes X x,∆t k using different time discretizations). Here, we rely on the sequential control variates algorithm proposed in [START_REF] Gobet | A spectral monte carlo method for the poisson equation[END_REF] and analyzed in [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. This algorithm constructs a sequence of approximations of u. At each iteration of the algorithm, the current approximation is used as a control variate for the estimation of u through Feynman-Kac formula.

A Sequential Control Variates Algorithm

Here we recall the sequential control variates algorithm introduced in [START_REF] Gobet | A spectral monte carlo method for the poisson equation[END_REF] in a general interpolation framework. We let V Λ ⊂ C 2 (D) be an approximation space of finite dimension #Λ and let I Λ : R D → V Λ be the interpolation operator associated with a unisolvent grid Γ Λ = {x ν : ν ∈ Λ}. We let (l ν ) ν∈Λ denote the (unique) basis of V Λ that satisfies the interpolation property l ν (x µ ) = δ νµ for all ν, µ ∈ Λ. The interpolation I Λ (w) = ν∈Λ w(x ν )l ν (x) of function w is then the unique function in

V Λ such that I Λ (w)(x ν ) = w(x ν ), ν ∈ Λ.
The following algorithm provides a sequence of approximations (ũ k ) k≥1 of u in V Λ , which are defined by ũk = ũk-1 + ẽk , where ẽk is an approximation of e k , solution of

A(e k )(x) = g(x) -A(ũ k-1 )(x), x ∈ D, e k (x) = f (x) -ũk-1 (x), x ∈ ∂D.
Note that e k admits a Feyman-Kac representation e k (x) = E(F (e k , X x )), where F (e k , X x ) depends on the residuals g -A(ũ k-1 ) on D and f -ũk-1 on ∂D. The approximation ẽk is then defined as the interpolation I Λ (e k ∆t,M ) of the Monte-Carlo estimate e k ∆t,M (x) of e k ∆t (x) = E(F (e k , X x,∆t )) (using M samples of X x,∆t ). For Algorithm 2.1 Sequential Control Variates Algorithm [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF] 1: Set ũ0 = 0, k = 1 and S = 0. 2: while k ≤ K and S < n s do 3:

Compute e k ∆t,M (x ν ) for x ν ∈ Γ Λ .

4:

Compute ẽk = I Λ (e k ∆t,M ) = ν∈Λ e k ∆t,M (x ν )l ν (x).

5:

Update ũk = ũk-1 + ẽk .

6:

If ũkũk-1 2 ≤ tol ũk-1 2 then S = S + 1 else S = 0.

7:

Set k = k + 1.

8: end while practical reasons, Algorithm 2.1 is stopped using an heuristic error criterion based on stagnation. This criterion is satisfied when the desired tolerance tol is reached for n s successive iterations (in practice we chose n s = 5). Now let us provide some convergence results for Algorithm 2.1. To that goal, we introduce the time integration error at point x for a function h

(2.8) e ∆t (h, x) = E[F (h, X x,∆t )] -E[F (h, X x )].
Then the following theorem [49, Theorem 3.1] gives a control of the error in expectation.

Theorem 2.3. Assuming (E2)-(E5), it holds

sup ν∈Λ E ũn+1 (x ν ) -u(x ν ) C(∆t, Λ) sup ν∈Λ |E [ũ n (x ν ) -u(x ν )]| + C 1 (∆t, Λ) with C(∆t, Λ) = sup ν∈Λ µ∈Λ |e ∆t (l µ , x ν )| and C 1 (∆t, Λ) = sup ν∈Λ e ∆t (u -I Λ (u), x ν ) . Moreover if C(∆t, Λ) < 1, it holds (2.9) lim sup n→∞ sup ν∈Λ |E [ũ n (x ν ) -u(x ν )]| C 1 (∆t, Λ) 1 -C(∆t, Λ)
.

The condition C(∆t, Λ) < 1 implies that in practice ∆t should be chosen sufficiently small [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]Theorem 4.2]. Under this condition, the error at interpolation points uniformly converges geometrically up to a threshold term depending on time integration errors for interpolation functions l ν and the interpolation error u-I Λ (u).

Theorem 2.3 provides a convergence result at interpolation points. Below, we provide a corollary to this theorem that provides a convergence result in L ∞ (D). This result involves the Lebesgue constants in L ∞ -norm associated to I Λ , defined by (2.10)

L Λ = sup v∈C 0 (D) I Λ (v) ∞ v ∞ ,
and such that for any v ∈ C 0 (D),

(2.11) v -I Λ (v) ∞ ≤ (1 + L Λ ) inf w∈V Λ v -w ∞ .
Throughout this article, we adopt the convention that supremum exclude elements with norm 0. We recall also that the L ∞ Lebesgue constant can be expressed as

L Λ = sup x∈D ν∈Λ |l ν (x)|. Corollary 2.4 (Convergence in L ∞ ). Assuming (E2)-(E5), it holds (2.12) lim sup n→∞ E [ũ n -u] ∞ C 1 (∆t, Λ) 1 -C(∆t, Λ) L Λ + u -I Λ (u) ∞ .
Proof. By triangular inequality, we have

E [ũ n -u] ∞ E [ũ n -I Λ (u)] ∞ + I Λ (u) -u ∞ . We can build a continuous function w such that w(x ν ) = E [ũ n (x ν ) -u(x ν )] for all ν ∈ Λ, and such that w ∞ = sup ν∈Λ |w(x ν )| = sup ν∈Λ |E [ũ n (x ν ) -u(x ν )]| .
We have then

E [ũ n -I Λ (u)] ∞ = I Λ (w) ∞ ≤ L Λ w ∞ .
The result follows from the definition of the function w and Theorem 2.3.

Remark 2.5. Since for bounded domains D, we have

v 2 ≤ |D| 1/2 v ∞ ,
for all v in C 0 (D), where |D| denotes the Lebesgue measure of D, we can deduce the convergence results in L 2 norm from those in L ∞ norm.

Adaptive Sparse Interpolation

We here present sparse interpolation methods following [START_REF] Chkifa | Sparse adaptive taylor approximation algorithms for parametric and stochastic elliptic pdes[END_REF][START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF].

Sparse Interpolation

For 1 ≤ i ≤ d, we let {ϕ (i) k } k∈N 0 be a univariate polynomial basis, where ϕ

(i) k (x i ) is a polynomial of degree k. For a multi-index ν = (ν 1 , . . . , ν d ) ∈ N d 0 , we introduce the multivariate polynomial ϕ ν (x) = d i=1 ϕ (i) ν i (x i ).
For a subset Λ ⊂ N d , we let

P Λ = span{ϕ ν : ν ∈ Λ}. A subset Λ is said to be downward closed if ∀ν ∈ Λ, µ ≤ ν ⇒ µ ∈ Λ.
If Λ is downward closed, then the polynomial space P Λ does not depend on the choice of univariate polynomial bases and is such that P Λ = span{x ν : ν ∈ Λ}, with

x ν = x ν 1 1 . . . x ν d d .
In the case where D = D 1 × . . . × D d , we can choose for {ϕ

(i) k } k∈N 0 an orthonormal basis in L 2 (D i ) (i.
e. a rescaled and shifted Legendre basis). Then {ϕ ν } ν∈N d 0 is an orthonormal basis of L 2 (D). To define a set of points Γ Λ unisolvent for P Λ , we can proceed as follows. For each dimension 1 ≤ i ≤ d, we introduce a sequence of points {z

(i) k } k∈N 0 in D i such that for any p ≥ 0, Γ (i) p = {z (i) k } p k=0 is unisolvent for P p = span{ϕ (i) k : 0 ≤ k ≤ p}, therefore defining an interpolation operator I (i) p . Then we let Γ Λ = {z ν = (z (1) ν 1 , . . . , z (d) ν d ) : ν ∈ Λ} ⊂ D.
This construction is interesting for adaptive sparse algorithms since for an increasing sequence of subsets Λ n , we obtain an increasing sequence of sets Γ Λn , and the computation of the interpolation on P Λn only requires the evaluation of the function on the new set of points Γ Λn \ Γ Λ n-1 . Also, with such a construction, we have the following property of the Lebesgue constant of I Λ in L ∞ -norm. This result is directly taken from [22, Section 3].

Proposition 2.6. If for each dimension

1 ≤ i ≤ d, the sequence of points {z (i) k } k∈N 0 is such that the interpolation operator I (i)
p has a Lebesgue constant L p ≤ (p + 1) s for some s > 0, then for any downward closed set Λ, the Lebesgue constant L Λ satisfies

(2.13) L Λ ≤ (#Λ) s+1 .
Leja points or magic points [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF] are examples of sequences of points such that the interpolation operators I (i) p have Lebesgue constants not growing too fast with p. For a given Λ with ρ i := max ν∈Λ ν i , it is possible to construct univariate interpolation grids Γ (i) ρ i with better properties (e.g., Chebychev points), therefore resulting in better properties for the associated interpolation operator I Λ . However for Chebychev points, e.g., ρ i ≤ ρ i does not ensure Γ (i)

ρ i ⊂ Γ (i) ρ i
. Thus with such univariate grids, an increasing sequence of sets Λ n will not be associated with an increasing sequence of sets Γ Λn , and the evaluations of the function will not be completely recycled in adaptive algorithms. However, for some of the algorithms described in Section 2.4, this is not an issue as evaluations can not be recycled anyway. Note that for general domains D which are not the product of intervals, the above constructions of grids Γ Λ are not viable since it may yield to grids not contained in the domain D. For such general domains, magic points obtained through greedy algorithms could be considered.

Adaptive Algorithm for Sparse Interpolation

An adaptive sparse interpolation algorithm consists in constructing a sequence of approximations (u n ) n≥1 associated with an increasing sequence of downward closed subsets (Λ n ) n≥1 . According to (2.11), we have to construct a sequence such that the best approximation error and the Lebesgue constant are such that

L Λn inf w∈P Λn u -w ∞ -→ 0 as n → ∞
for obtaining a convergent algorithm. For example, if (2.14) inf

w∈P Λn u -w ∞ = O((#Λ n ) -r )
holds [3] for some r > 1 and if

L Λn = O((#Λ n ) k ) for k < r, then the error u-u n ∞ = O(n -r
) tends to zero with an algebraic rate of convergence r = r-k > 0. Of course, the challenge is to propose a practical algorithm that constructs a good sequence of sets Λ m . We now present the adaptive sparse interpolation algorithm with bulk chasing procedure introduced in [START_REF] Chkifa | Sparse adaptive taylor approximation algorithms for parametric and stochastic elliptic pdes[END_REF]. Let θ be a fixed bulk chasing parameter in (0, 1) and let E Λ (v) = P Λ (v) Compute M Λn .

4:

Set Λ n = Λ n ∪ M Λn and compute I Λ n (u).

5:

Select N n ⊂ M Λn the smallest such that E Nn (I Λ n (u)) ≥ θE M Λn (I Λ n (u)) 6:
Update Λ n+1 = Λ n ∪ N n .

7:

Compute u n+1 = I Λ n+1 (u) (this step is not necessary in practice).

8:

Compute ε n .

9:

Update n = n + 1.

10: end while reduced margin of Λ n defined by

M Λn = {ν ∈ N d \ Λ n : ∀j s.t. ν j > 0, ν -e j ∈ Λ n },
where (e j ) k = δ kj . The reduced margin is such that for any subset S ⊂ M Λn , Λ n ∪ S is downward closed. This ensures that the sequence (Λ n ) n≥1 generated by [3] see e.g. [START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF] for conditions on u ensuring such a behavior of the approximation error.

the algorithm is an increasing sequence of downward closed sets. Finally, Algorithm 2.2 is stopped using a criterion based on

ε n = E Mn (I Λ n (u)) E Λ n (I Λ n (u)) .

Combining Sparse Adaptive Interpolation with Sequential Control Variates Algorithm

We present in this section two ways of combining Algorithm 2.1 and Algorithm 2.2. First we introduce a perturbed version of Algorithm 2.2 and then an adaptive version of Algorithm 2.1. At the end of the section, numerical results will illustrate the behavior of the proposed algorithms.

Perturbed Version of Sparse Interpolation

As we do not have access to exact evaluations of the solution u of (2.1), Algorithm 2.2 can not be used for interpolating u. So we introduce a perturbed version of this algorithm, where the computation of the exact interpolant I Λ (u) is replaced by an approximation denoted ũΛ , which can be computed for example with Algorithm 2.1 stopped for a given tolerance tol or at step k. This brings the following algorithm. Compute M Λn .

4:

Set Λ n = Λ n ∪ M Λn and compute ũΛ n .

5:

Select N n as the smallest subset of M Λn such that E Nn (ũ

Λ n ) ≥ θE M Λn (ũ Λ n ) 6:
Update Λ n+1 = Λ n ∪ N n .

7:

Compute ũΛ n+1 .

8:

Compute εn .

9:

Update n = n + 1.

10: end while

Adaptive Version of Sequential Control Variate Algorithm

As a second algorithm, we consider the sequential control variates algorithm (Algorithm 2.1) where at step 4, an approximation ẽk of e k is obtained by applying the adaptive interpolation algorithm (Algorithm 2.3) to the function e k ∆t,M , which uses Monte-Carlo estimations e k ∆t,M (x ν ) of e k (x ν ) at interpolation points. At each iteration, ẽk therefore belongs to a different approximation space P Λ k . In the numerical section, we will call this algorithm adaptive Algorithm 2.1.

Numerical Results

In this section, we illustrate the behavior of algorithms previously introduced on different test cases. We consider the simple diffusion equation

(2.15) -u(x) = g(x), x ∈ D, u(x) = f (x), x ∈ ∂D, were D =] -1, 1[ d .
The source terms and boundary conditions will be specified later for each test case. The stochastic differential equation associated to (2.15) is the following (2. [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF])

dX x t = √ 2dW t , X x 0 = x, where (W t ) t≥0 is a d-dimensional Brownian motion.
We use tensorized grids of magic points for the selection of interpolation points evolved in adaptive algorithms.

Small dimensional test case. We consider a first test case (TC1) in dimension d = 5.

Here the source term and the boundary conditions in problem (2.15) are chosen such that the solution is given by Step k u -ũk Step k u -ũk We first test the influence of ∆t and M on the convergence of Algorithm 2.1 when Λ is fixed. In that case, Λ is selected a priori with Algorithm 2.2 using samples of the exact solution u for (TC1), stopped for ε ∈ {10 -6 , 10 -8 , 10 -10 }. In what follows, the notation Λ i stands for the set obtained for ε = 10 -i , i ∈ {6, 8, 10}. We represent on Figure 2.1 the evolution of the absolute error in L 2 -norm (similar results hold for the L ∞ -norm) between the approximation and the true solution with respect to step k for Λ = Λ 6 . As claimed in Corollary 1, we recover the geometric convergence up to a threshold value that depends on ∆t. We also notice faster convergence as M increases and when ∆t decreases. We fix M = 1000 in the next simulations. We study the impact of the choice of Λ i on the convergence of Algorithm 2.1. Again we observe on Figure 2.2 that the convergence rate gets better as ∆t decreases. Moreover as #Λ increases the threshold value decreases. This is justified by the fact that interpolation error decreases as #Λ i increases (see Table 2.1). Nevertheless, we observe that it may also deteriorate the convergence rate if it is chosen too large together with ∆t not sufficiently small. Indeed for the same number of iterations k = 10 and the same time-step ∆t = 2.5 • 10 -3 , we have an approximate absolute error equal to 10 -7 for Λ 8 against 10 -4 for Λ 10 . We present then results provided with the adaptive Algorithm 2.1. The parameters chosen for the adaptive interpolation are ε = 5 • 10 -2 , θ = 0.5. K = 30 ensures the stopping of Algorithm 2.1. As illustrated by Figure 2.4, we recover globally the same behavior as for Algorithm 2.1 without adaptive interpolation. Indeed as k increases, both absolute errors in L 2 -norm and L ∞ -norm decrease and then stagnate. Again, we notice the influence of ∆t on the stagnation level. Nevertheless, the convergence rates are deteriorated and the algorithm provides less accurate approximations than Algorithm 2.3. This might be due to the sparse adaptive interpolation procedure, which uses here pointwise evaluations based on Monte-Carlo estimates, unlike Algorithm 2.3 which relies on pointwise evaluations resulting from Algorithm 2.1 stopping for a given tolerance.

(TC1) u(x) = x 2 1 + sin(x 2 ) + exp(x 3 ) + sin(x 4 )(x 5 + 1), x ∈ D.
Λ n #Λ n ε n ||u -u n || 2 ||u -u n || ∞ 1 6.183372e
Finally in Table 2.2, we compare the algorithmic complexity of these algorithms to reach a precision of 3 • 10 -5 for (∆t, M ) = (10 -4 , 1000). For adaptive Algorithm 2.1, Λ k refers to the set of multi-indices considered at step k of Algorithm 2.1. For Algorithm 2.3, N n stands for the number of iteration required by Algorithm 2.1 to reach tolerance ε tol at step n. Finally, Algorithm 2.1 is run with full-grid Λ = Λ max where Λ max = {ν ∈ N d : ν i ≤ 10} is the set of multi-indices allowing to reach the machine precision. In this case, N stands for the number of steps for this algorithm to converge. Higher-dimensional test cases. Now, we consider two other test cases noted respectively (TC2) and (TC3) in higher dimension.

(∆t) -1 ( k #Λ k ) M (∆t) -1 ( n #Λ n N n ) M (∆t) -1 #Λ max N Est.
(TC2) As second test case in dimension d = 10, we define (2.15) such that its solution is the Henon-Heiles potential

u(x) = 1 2 d i=1 x 2 i + 0.2 d-1 i=1 x i x 2 i+1 -x 3 i + 2.5 10 -3 d-1 i=1 x 2 i + x 2 i+1 2 , x ∈ D.
We set (∆t, M ) = (10 -4 , 1000) and K = 30 for Algorithm 2.1.

(TC3) We also consider the problem (2.15) whose exact solution is a sum of nonpolynomial functions, like (TC1) but now in dimension d = 20, given by u(x) = x 2 1 + sin(x 12 ) + exp(x 5 ) + sin(x 15 )(x 8 + 1).

Here, the Monte-Carlo simulations are performed for (∆t, M ) = (10 -4 , 1000) and K = 30.

Since for both test cases the exact solution is known, we propose to compare the behavior of Algorithm 2.3 and Algorithm 2.2. Again, the approximations ũn , at each step n of Algorithm 2.3, are provided by Algorithm 2.1 stopped when a stagnation is detected. In both cases, the parameters for Algorithm 2.3 are set to θ = 0.5 and

ε = 10 -15 .
In Table 2.3 and Table 2.4, we summarize the results associated to the exact and perturbed sparse adaptive algorithms for (TC2) and (TC3) respectively. We observe that Algorithm 2.3 performs well in comparison to Algorithm 2.2, for (TC2). Indeed, we get an approximation with a precision below the prescribed value ε for both algorithms. Similar observation holds for (TC3) in Table 2.4 and this despite the fact that the test case involves higher dimensional problem.

#Λ n ε n u -u n ∞ u -u n 2 #Λ n εn u -ũΛn ∞ u -
#Λ n ε n u -u n ∞ u -u n 2 #Λ n εn u -ũΛn ∞ u -

Conclusion

In this paper we have introduced a probabilistic approach to approximate the solution of high-dimensional elliptic PDEs. This approach relies on adaptive sparse polynomial interpolation using pointwise evaluations of the solution estimated using a Monte-Carlo method with control variates. Especially, we have proposed and compared different algorithms. First we proposed Algorithm 2.1 which combines the sequential algorithm proposed in [START_REF] Gobet | A spectral monte carlo method for the poisson equation[END_REF] and sparse interpolation. For the non-adaptive version of this algorithm we recover the convergence up to a threshold as the original sequential algorithm [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. Nevertheless it remains limited to small-dimensional test cases, since its algorithmic complexity remains high. Hence, for practical use, the adaptive Algorithm 2.1 should be preferred. Adaptive Algorithm 2.1 converges but it does not allow to reach low precision with reasonable number of Monte-Carlo samples or time-step in the Euler-Maruyama scheme. Secondly, we proposed Algorithm 2.3. It is a perturbed sparse adaptive interpolation algorithm relying on inexact pointwise evaluations of the function to approximate. Numerical experiments have shown that the perturbed algorithm (Algorithm 2.3) performs well in comparison to the ideal one (Algorithm 2.2) and better than the adapted Algorithm 2.1 with a similar algorithmic complexity. Here, since only heuristic tools have been provided to justify the convergence of this algorithm, the proof of its convergence, under assumptions on the class of functions to be approximated, should be addressed in a future work.

Chapter 3 A PAC Algorithm in Relative Precision for Bandit Problem with Costly Sampling

This chapter introduces our second contribution. It consists in a discrete optimization algorithm whose output is probably quasi-optimal in relative precision, with controlled probability and relative quasi-optimality. We present in this chapter a contribution that has been submitted as [START_REF] Billaud-Friess | A pac algorithm in relative precision for bandit problem with costly sampling[END_REF] in Mathematical Methods of Operations Research and is currently under review.

This contribution is about discrete optimization for functions defined using a expectation. The probably approximately correct (PAC) algorithm in relative precision presented and analysed in this chapter will then be used in Chapter 5 as a tool to select a quasi-optimal snapshot to enrich the reduced space.

Title: A PAC algorithm in relative precision for bandit problem with costly sampling

Authors: Marie Billaud-Friess ♣ , Arthur Macherey ♣,♠ , Anthony Nouy ♣ and Clémentine Prieur ♠ .

Introduction

We consider an optimization problem

(3.1) max ξ∈Ξ E[Z(ξ)],
where E[Z(ξ)] is the expectation of a random variable Z(ξ), and where we assume that the set Ξ is finite. Such a problem is encountered in different fields such as reinforcement learning [START_REF] Sutton | Reinforcement learning: An introduction[END_REF] or robust optimization [START_REF] Beyer | Robust optimization-a comprehensive survey[END_REF].

To solve (3.1), classical optimization methods include random search algorithms [START_REF] Gong | Stochastic comparison algorithm for discrete optimization with estimation[END_REF][START_REF] Yan | Stochastic discrete optimization[END_REF], stochastic approximation methods [START_REF] Dupa | Stochastic approximation on a discrete set and the multi-armed[END_REF][START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF] and bandit algorithms [START_REF] Lattimore | Bandit algorithms[END_REF]2,3,[START_REF] Garivier | The kl-ucb algorithm for bounded stochastic bandits and beyond[END_REF].

In this paper, we focus on unstructured stochastic bandit problems with a finite number of arms where "arms" stands for "random variables" and corresponds here to the Z(ξ), ξ ∈ Ξ (see, e.g., [START_REF] Lattimore | Bandit algorithms[END_REF]Section 4]). Stochastic means that the only way to learn about the probability distribution of arms Z(ξ), ξ ∈ Ξ is to generate i.i.d. samples from it. Unstructured means that knowledge about the probability distribution of one arm Z(ξ) does not restrict the range of possibilities for other arms

Z(ξ ), ξ = ξ.
Additionally, we suppose here it is numerically costly to sample the random variables Z(ξ), ξ ∈ Ξ. Our aim is thus to solve (3.1) by sampling as few as possible the random variables Z(ξ), ξ ∈ Ξ. However, it is not feasible to solve (3.1) almost surely using only a finite number of samples from the random variables Z(ξ), ξ ∈ Ξ. Thus, it is relevant to adopt a probably approximately correct (PAC) approach (see e.g. [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF]). For a precision τ abs and a probability λ ∈ (0, 1), a (τ abs , λ)-PAC algorithm returns ξ such that

(3.2) P E[Z(ξ )] -E[Z( ξ)] ≤ τ abs ≥ 1 -λ, ξ ∈ arg max ξ∈Ξ E[Z(ξ)].
Until recently, one of the main focus of bandit algorithms was the best arm (random variable) identification [2], through the use of Successive Rejects algorithm or Up-per Confidence Bounds algorithms. Such algorithms are (0, λ)-PAC algorithms, as stated in [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF]. Racing algorithms [3] were designed to solve the best arm identification problem too and are mainly analyzed in a finite budget setting, which consists in fixing a maximum number of samples that can be used. While trying to identify the best arm, bandit algorithms also aim at minimizing the regret [1,[START_REF] Bubeck | Pure exploration in finitely-armed and continuous-armed bandits[END_REF][START_REF] Garivier | The kl-ucb algorithm for bounded stochastic bandits and beyond[END_REF]. More recently, other focuses have emerged, such as the identification of the subset of Ξ containing the m best arms [START_REF] Kalyanakrishnan | Pac subset selection in stochastic multi-armed bandits[END_REF][START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF] or the identification of "good arms" (also known as thresholding bandit problem) that are random variables whose expectation is greater or equal to a given threshold [START_REF] Kano | Good arm identification via bandit feedback[END_REF][START_REF] Tao | Thresholding bandit with optimal aggregate regret[END_REF][START_REF] Mukherjee | Thresholding bandits with augmented ucb[END_REF][START_REF] Locatelli | An optimal algorithm for the thresholding bandit problem[END_REF].

The (τ abs , λ)-PAC algorithms mentioned above measure the error in absolute precision. However, without knowing E[Z(ξ )], providing in advance a relevant value for τ abs is not an easy task. In this work, we rather consider (τ, λ)-PAC algorithms in relative precision that return ξ ∈ Ξ such that

(3.3) P E[Z(ξ )] -E[Z( ξ)] ≤ τ |E[Z(ξ )]| ≥ 1 -λ,
where τ and λ are set in advance in (0, 1). We introduce two algorithms that yield a solution ξ satisfying (3.3). The first algorithm builds an estimate precise enough for each expectation E[Z(ξ)]. This naive approach drives a majority of the budget on the random variables with the lowest expectations in absolute value. In order to avoid this drawback and thus to reduce the number of samples required to reach the prescribed relative precision, we propose a second algorithm which adaptively samples random variables exploiting confidence intervals obtained from an empirical Berstein concentration inequality.

The outline of the paper is as follows. In section 3.2, we present a Monte-Carlo estimate for the expectation of a single random variable that has been proposed in [START_REF] Mnih | Efficient stopping rules[END_REF]. It provides an estimation of the expectation with guaranteed relative precision, with high probability. In section 3.3, we introduce two new algorithms that rely on these Monte-Carlo estimates and yield a solution to (3.3). Then, we study numerically the performance of our algorithms and compare them to algorithms from the literature, possibly adapted to solve (3.3).

Monte-Carlo Estimate with Guaranteed Relative Precision

In what follows, we consider a random variable Z defined on probability space (Ω, F, P). We denote by Z m the empirical mean of Z and by V m its empirical variance, respectively defined by

Z m = 1 m m i=1 Z i and V m = 1 m m i=1 Z i -Z m 2 ,
where (Z i ) i≥1 is a sequence of i.i.d. copies of Z. The aim is to provide an estimate

Ê[Z] of E[Z] which satisfies (3.4) P | Ê[Z] -E[Z]| ≤ |E[Z]| ≥ 1 -δ,
Chapter 3. PAC Algorithm in Relative Precision with ( , δ) ∈ (0, 1) 2 given a priori. For that, we will rely on Theorem 3.1 hereafter.

Theorem 3.1. If Z takes its values in a bounded interval [a, b],

for any m ∈ N and x ∈ (0, 1), we have

(3.5) P   Z m -E[Z] ≤ 2V m log(3/x) m + 3(b -a) log(3/x) m   ≥ 1 -x.
Proof. We simply apply [1, Theorem 1] to Z -a which is a positive random variable whose values are lower than b -a.

Based on Theorem 3.1, several estimates for E(Z) have been proposed in [START_REF] Mnih | Efficient stopping rules[END_REF][START_REF] Mnih | Empirical bernstein stopping[END_REF].

We focus in this paper on the estimate introduced in [78, Equation (3.7)].

Monte-Carlo Estimate

Considering a sequence (d m ) m≥1 in (0, 1), we introduce the sequence (c m ) m≥1 defined, for all m ≥ 1, by Let M be an integer-valued random variable on (Ω, F, P), such that

c m = 2V m log(3/d m ) m + 3(b -a) log(3/d m ) m . ( 3
(3.8) c M ≤ |Z M |,
with ∈ (0, 1). Then, we define the following estimate

ÊM [Z] = Z M -sign(Z M )c M . (3.9) Proposition 3.2. Let , δ ∈ (0, 1). Assume that Z takes its values in a bounded interval [a, b] and that (d m ) m≥1 satisfies (3.10) m≥1 d m ≤ δ.
Then the estimate ÊM [Z] defined by (3.9), with M satisfying (3.8), is such that

(3.11) P ÊM [Z] -E[Z] ≤ |E[Z]| ≥ 1 -δ.

Complexity Analysis

In this section, we state a complexity result. Following [START_REF] Mnih | Efficient stopping rules[END_REF], we focus on a particular sequence (d m ) m≥1 defined by

ν = min max(σ 2 , 2 µ 2 ) (b -a) 2 , 2 µ 2 (1 + ) 2 max(σ 2 , 2 µ 2 )γ , with γ = ( 2 + 2 √ 2 + 2/3 + 3) 2 . Moreover, E(M ) ≤ 2 ν p log 2p ν + log 3 cδ + 4δ/3.
Proof. See Appendix.

Remark 3.5. The result from Proposition 3.4 helps in understanding the influence of parameters (ε, δ) appearing in (3.4) on M . Indeed, we deduce from this result that for δ < 1/2,

E(M ) ν -1 log(ν -1 ) + (ν -1 + 1) log(δ -1 ).
We first observe a weak impact of δ on the average complexity. When → 0, we have ν ∼ 2 µ 2 σ 2 γ . Then for fixed δ and → 0, the bound for

E(M ) is in O( -2 σ 2 µ 2 ).
As expected, the relative precision has a much stronger impact on the average complexity.

Optimization Algorithms with Guaranteed Relative Precision

In this section we consider a finite collection of bounded random variables Z(ξ) on (Ω, F, P), indexed by ξ ∈ Ξ, and such that E[Z(ξ)] = 0. Each random variable Z(ξ) takes its values in a bounded interval [a(ξ), b(ξ)], which is assumed to be known.

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

(3.18) Z(ξ) m = 1 m m i=1 Z(ξ) i and V (ξ) m = 1 m m i=1 Z(ξ) i -Z(ξ) m 2 ,
where {(Z(ξ) i ) i≥1 : ξ ∈ Ξ} are independent i.i.d. copies of {Z(ξ) : ξ ∈ Ξ}. We also introduce #Ξ different sequences

c ξ,m = 2V (ξ) m log(3/d m ) m + 3 (b(ξ) -a(ξ)) log(3/d m ) m ,
where (d m ) m≥1 is a positive sequence, independent from ξ, such that m≥1 d m ≤ δ.

Taking in (0, 1), for each ξ in Ξ, we define, as in (3.14),

(3.19) m(ξ) = min m ∈ N : c ξ,m ≤ |Z(ξ) m | .
Then defining s(ξ) := sign(Z(ξ) m(ξ) ), we propose the following estimate for E[Z(ξ)]:

(3.20) Êm(ξ) [Z(ξ)] = Z(ξ) M (ξ) -s(ξ)c ξ,m(ξ) .
These notation being introduced, we propose below two algorithms returning ξ in Ξ such that

(3.21) P E[Z(ξ )] -E[Z( ξ)] ≤ τ |E[Z(ξ )]| ≥ 1 -λ, ξ ∈ arg max ξ∈Ξ E[Z(ξ)],
for given (τ, λ) in (0, 1) 2 .

Non-Adaptive Algorithm

We first propose a non-adaptive algorithm that provides a parameter ξ satisfying (3.21), by selecting the maximizer of independent estimates Êm(ξ) [Z(ξ)] of E[Z(ξ)] over Ξ.

Algorithm 3.1 Non-Adaptive PAC Algorithm in Relative Precision

Require: τ , λ, {Z(ξ)} ξ∈Ξ . Ensure: ξ Then, for all ξ in Ξ, the estimate Êm(ξ) [Z(ξ)] is well defined and satisfies (3.23) ,P(A(ξ)) ≥ 1λ #Ξ and by a union bound argument, P(∩ ξ∈Ξ A(ξ)) ≥ 1 -λ. To prove that ξ satisfies (3.21), it remains to prove that

(3.23) P E[Z(ξ)] -Êm(ξ) [Z(ξ)] ≤ |E[Z(ξ)]| ≥ 1 -δ, with = τ 2+τ .
A(ξ) = E [Z(ξ)] -Êm(ξ) [Z(ξ)] ≤ |E [Z(ξ)]| . By
∩ ξ∈Ξ A(ξ) implies E [Z(ξ )] -E[Z( ξ)] ≤ τ |E [Z(ξ )]|.
In what follows we suppose that ∩ ξ∈Ξ A(ξ) holds. Since < 1, E[Z(ξ)], Êm(ξ) [Z(ξ)] and Z(ξ) m(ξ) have the same sign, that we denote by s(ξ). Since A(ξ ) ∩ A( ξ) holds, we have

E [Z(ξ )] -E[Z( ξ)] ≤ E [Z(ξ )] - Êm( ξ) [Z( ξ)] 1 + s( ξ) ≤ E [Z(ξ )] - Êm(ξ ) [Z(ξ )] 1 + s( ξ) ≤ E [Z(ξ )] - 1 -s(ξ ) 1 + s( ξ) E [Z(ξ )] = (s(ξ ) + s( ξ)) 1 + s( ξ) E [Z(ξ )] .
Then we deduce

(3.24) E [Z(ξ )] -E[Z( ξ)] ≤ 2 1 - |E[Z(ξ )]| = τ |E[Z(ξ )]|,
which ends the proof.

Remark 3.7. If E[Z(ξ )] > 0, we can prove that s( ξ) = s(ξ ) = 1, so that the inequality (3.24) becomes E [Z(ξ )] -E[Z( ξ)] ≤ 2 1 + |E[Z(ξ )]|.
Therefore, we can set = τ 2-τ in Algorithm 3.1 to lower the complexity and still guarantee that ξ satisfies (3.21). Algorithm 3.1 provides for each random variable an estimate Êm(ξ) [Z(ξ)] that satisfies (3.23). However, as will be illustrated later, this algorithm tends to use many samples for variables with a low expectation in absolute value. We propose in the next subsection an adaptive algorithm avoiding this drawback by using confidence intervals, which results in a lower overall complexity.

Adaptive Algorithm

The idea of the adaptive algorithm is to successively increase the number of samples m(ξ) of a subset of random variables Z(ξ) that are selected based on confidence intervals of E[Z(ξ)] deduced from the concentration inequality of Theorem 3.1. This algorithm follows the main lines of the racing algorithms [79, section 4]. However racing algorithms do not allow to sample again a random variable discarded in an earlier step of the algorithm. The adaptive algorithm presented hereafter allows it.

In order to present this adaptive algorithm, for each ξ, we introduce the confidence interval

[β - m(ξ) (ξ), β + m(ξ) (ξ)], with (3.25) β - ξ,m(ξ) = Z(ξ) m(ξ) -c ξ,m(ξ) and β + ξ,m(ξ) = Z(ξ) m(ξ) + c ξ,m(ξ) . From (3.7), we have that (3.26) P β - ξ,m(ξ) ≤ E(Z(ξ)) ≤ β + ξ,m(ξ) ≥ 1 -d m(ξ) . We define ξ,m(ξ) by ξ,m(ξ) = c ξ,m(ξ) |Z(ξ) m(ξ) | if Z(ξ) m(ξ) = 0, or ξ,m(ξ) = +∞ otherwise. Letting s(ξ) := sign(Z(ξ) m(ξ) ), we use as an estimate for E[Z(ξ)] (3.27) Êm(ξ) [Z(ξ)] = Z(ξ) m(ξ) -ξ,m(ξ) s(ξ)c ξ,m(ξ) if ξ,m(ξ) < 1, Z(ξ) m(ξ) otherwise.
If ξ,m(ξ) < 1, we note that

Êm(ξ) [Z(ξ)] = (Z(ξ) m(ξ) ∓ c ξ,m(ξ) ) 1 ± s(ξ) ξ,m(ξ) , so that (3.28) Êm(ξ) [Z(ξ)] 1 ± s(ξ) ξ,m(ξ) = β ∓ ξ,m(ξ) .
The adaptive algorithm is described in Algorithm 3.2. At each iteration n, one sample of Z(ξ) is drawn for each ξ in a subset Ξ n selected according to (3.29). In the next proposition, we prove that the algorithm returns a solution to (3.21) under suitable assumptions.

Proposition 3.8. Let (τ, λ) ∈ (0, 1) 2 . We assume that (d m ) m≥1 is a positive sequence that satisfies

(3.30) +∞ m=1 d m ≤ λ #Ξ and log(3/d m )/m → m→+∞ 0.
Moreover, we assume that, for all ξ in Ξ, Z(ξ) is a bounded random variable with E[Z(ξ)] = 0. Then, it holds almost surely that Algorithm 3.2 stops and ξ satisfies (3.21).

Algorithm 3.2 Adaptive PAC Algorithm in Relative Precision

Require: τ , λ, {Z(ξ)} ξ∈Ξ . Ensure: ξ. 1: Set n = 0, Ξ 0 = Ξ, ξ,0 = +∞ and m(ξ) = 0 for all ξ ∈ Ξ. for all ξ ∈ Ξ n do 4:

Sample Z(ξ), increment m(ξ) and update ξ,m(ξ) .

5:

Build the estimate Êm(ξ) [Z(ξ)] using (3.27).

6:

end for 7:

Increment n and put in Ξ n every ξ ∈ Ξ such that Proof. Let m n (ξ) denote the number of samples of Z(ξ) at iteration n of the algorithm. We first prove by contradiction that Algorithm 3.2 stops almost surely. Let us suppose that Algorithm 3.2 does not stop with probability η > 0, that means

(3.29) β + ξ,m(ξ) ≥ max ν∈Ξ β - ν,m(ν) (ν).
(3.31) P ∀n > 0, #Ξ n > 1 and max ξ∈Ξn ξ,mn(ξ) > τ 2 + τ = η > 0.
Since (Ξ n ) n≥1 is a sequence from a finite set, we can extract a constant subsequence, still denoted (Ξ n ) n≥1 , equal to Ξ ⊂ Ξ, with Ξ = ∅ such that we increase m n (ξ) for all ξ in Ξ , we have that m n (ξ) → ∞ as n → ∞ for all ξ ∈ Ξ . Therefore, lim n→+∞ max ξ∈Ξ ξ,mn(ξ) = 0 holds almost surely, which contradicts (3.32).

We now prove that ξ satisfies (3.21). For clarity, we remove the index n from m n (ξ).

Defining

A(ξ) = |Z(ξ) m(ξ) -E[Z(ξ)]| ≤ c ξ,m(ξ)
for all ξ in Ξ, we proceed as in (3.13) to obtain

P (A(ξ)) ≥ 1 -λ/#Ξ.

Thus, by a union bound argument

P (∩ ξ∈Ξ A(ξ)) ≥ 1 -λ. It remains to prove that ∩ ξ∈Ξ A(ξ) implies E[Z(ξ )] -E[Z( ξ)] ≤ τ |E[Z(ξ )]|
in order to prove that ξ satisfies (3.21). In the rest of the proof, we suppose that ∩ ξ∈Ξ A(ξ) holds. First, for all ξ / ∈ Ξ n , using (3.25), we have

(3.33) E[Z(ξ)] ≤ β + ξ,m(ξ) < max ν∈Ξ β - ν,m(ν) ≤ max ν∈Ξ E[Z(ν)] = E[Z(ξ )],
that implies ξ ∈ Ξ n . If the stopping condition is #Ξ n = 1, we then have ξ = ξ . If the stopping condition is max ξ∈Ξn ξ,m(ξ) ≤ τ 2+τ < 1, it means that, for all ξ in Ξ n , ξ,m(ξ) ≤ τ 2+τ < 1. Then for all ξ ∈ Ξ n , using Proposition 3.2 with = ξ,m(ξ) < 1 and δ = λ/#Ξ < 1 and the fact that ∩ ξ∈Ξ A(ξ) holds, we obtain that the estimate Êm(ξ) [Z(ξ)] satisfies

(3.34) Êm(ξ) [Z(ξ)] -E[Z(ξ)] ≤ ξ,m(ξ) |E[Z(ξ)]|.
We have that ξ,m(ξ) < 1 and (3.28) hold for all ξ ∈ Ξ n . In particular, since ξ, ξ ∈ Ξ n we get

E [Z(ξ )] -E[Z( ξ)] ≤ E [Z(ξ )] - Êm( ξ) [Z( ξ)] 1 + s( ξ) ξ,m( ξ) ≤ E [Z(ξ )] - Êm(ξ ) [Z(ξ )] 1 + s( ξ) ξ,m( ξ) ≤ E [Z(ξ )] - 1 -s(ξ ) ξ ,m(ξ ) 1 + s( ξ) ξ,m( ξ) E [Z(ξ )] = s(ξ ) ξ ,m(ξ ) + s( ξ) ξ,m( ξ) 1 + s( ξ) ξ,m( ξ) E [Z(ξ )] .
Then we deduce

(3.35) E [Z(ξ )] -E[Z( ξ)] ≤ 2τ /(2 + τ ) 1 -τ /(2 + τ ) |E[Z(ξ )]| = τ |E[Z(ξ )]|,
which ends the proof.

Remark 3.9. As for Algorithm 3.1 (see Remark 3.7), if E[Z(ξ )] > 0, we can set = τ 2-τ in Algorithm 3.2 to lower the complexity and still guarantee that ξ satisfies (3.21). Remark 3.10. A variant of Algorithm 3.2 using batch sampling would generate several samples of Z(ξ) at step 4. The result of Proposition 3.8 also holds for the algorithm with batch sampling. An optimal choice of the number of samples should depend on sampling costs.

Numerical Results

In this section, we propose a numerical study of the behaviour of our algorithms on a toy example. We consider the set of random variables Z(ξ) = f (ξ) + U (ξ), ξ ∈ Ξ, where f (ξ) = sin(ξ) + sin(10ξ/3), the U (ξ) are i.i.d. uniform random variables over (-1/20, 1/20), and Ξ = {3 + 4i/100 : 0 ≤ i ≤ 100}. The numerical results are obtained with the sequence (d m ) m≥1 defined by (3.16) with p = 2. We set τ = 0.1 and λ = 0.1.

We first compare our algorithms with two existing ones. The first one is the Median Elimination (ME) algorithm (see [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF] for a description of the algorithm), that solves problem (3.2). We take τ abs = τ |E[Z(ξ )]| to ensure ME algorithm provides a solution that also guarantees (3.3). Of course, this is not feasible in practice without knowing the solution of the optimization problem or at least a bound of |E[Z(ξ )]|. The second algorithm which we compare to our algorithms is the UCB-V Algorithm (see [1, section 3.1]). It consists in only resampling the random variable whose confidence interval has the highest upper bound. To do so, we replace Steps 3 to 6 of Algorithm 3.2 by: Compute ξ + = arg max ξ∈Ξn β + ξ,m(ξ) , Sample Z(ξ + ), increment m(ξ + ) and update ξ + ,m(ξ + ) .

We choose these algorithms to perform the comparison because i) ME Algorithm ensures theoretical guarantees similar to ours (although in absolute precision) and ii) the UCB-V Algorithm is optimal, in a sense that we will define later, for solving the optimization problem (3.1).

]

m(ξ) E[Z(ξ)] Êm(ξ) [Z(ξ)]
(a) Algorithm 3.1

] We illustrate on Figure 3.1 the behavior of algorithms. The results that we show on Figure 3.1 are the ones of a single run of each algorithm. On the left scale, we plot the estimates Êm(ξ) [Z(ξ)] as defined in (3.27) and the associated confidence intervals

m(ξ) E[Z(ξ)] Êm(ξ) [Z(ξ)], ξ ∈ Ξn Êm(ξ) [Z(ξ)], ξ / ∈ Ξn (b) Algorithm 3.
m(ξ) E[Z(ξ)] Êm(ξ) [Z(ξ)] (c) ME Algorithm
m(ξ) E[Z(ξ)] Êm(ξ) [Z(ξ)], ξ ∈ Ξn Êm(ξ) [Z(ξ)], ξ / ∈ Ξn (d) UCB-V Algorithm
[β - ξ,m(ξ) , β + ξ,m(ξ) ] of level 1 -d m(ξ)
given by (3.26). The estimates and confidence intervals for ξ ∈ Ξ n are drawn in blue, while the ones for ξ / ∈ Ξ n are drawn in red. On the right scale, we plot the number of samples m(ξ) generated for each ξ ∈ Ξ. We observe that Algorithm 3.1 samples too much the random variables with low expectation in absolute value. This is responsible for the three peaks on m(ξ) observed on Figure 3.1a. Algorithm 3.2 avoids this drawback as it does not try to reach the condition ξ,m(ξ) < 1 for all random variables. The UCB-V algorithm samples mostly the two random variables with highest expectations (more than 99% of the samples are drawn from these random variables). Other random variables are not sufficiently often sampled for reaching rapidly the stopping condition based on confidence intervals. The Median Elimination Algorithm oversamples all random variables in comparison with other algorithms.

Complexity.

To perform a quantitative comparison with existing algorithms in the case of costly sampling, a relevant complexity measure is the total number of samples generated after a single run of the algorithm

M = ξ∈Ξ m(ξ).
Table 3.1 shows the average complexity E(M) estimated using 30 independent runs of each algorithm. We observe that the expected complexity of Algorithm 3.2 is far below the one of the other algorithms. It means that, for the complexity measure E(M), the adaptive algorithm we have proposed performs the best. ME Alg. Alg. 3.1 Alg. 3.2 UCB-V Alg. 2.0 • 10 7 1.4 • 10 8 1.9 • 10 3 1.9 • 10 8

Table 3.1: Average complexity E(M), estimated using 30 runs for each algorithm, with τ = 0.1, λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME algorithm.

We now compare the four algorithms in terms of expected runtime, that is a measure of complexity taking into account the sampling cost and the cost of all other operations performed by the algorithms. Denoting by t the time (assumed constant) for generating one sample from a distribution, the runtime of an algorithm is a random variable T = Mt + N , where Mt is the sampling time, and N is the (random) time taken by all other operations. The expected runtime is then

E(T ) = E(M)t + E(N ).
From the values of E(N ) and E(M), estimated over 30 runs of the algorihms, we deduce Table 3.2, which shows the average runtime E(T ) for different values of t . We observe that Algorithm 3.2 has the smallest average runtime whatever the sampling cost. The first line corresponds to E(N ) and shows that Algorithm 3.2 performs the best when sampling cost t = 0 (or negligible). The impressive gain for large sampling costs t is due to the small value of the average number of samples E(M) required by the algorithm.

[Z(ξ)]

m1(ξ) E[Z(ξ)] Êm1(ξ) [Z(ξ)], ξ ∈ Ξ1 (a) n = 1 3 3.5 4 4.5 5 5.5 6 6.5 7 0 10 20 30 40 50 ξ m10(ξ) data4 -4 -3 -2 -1 0 1 2 Êm10(ξ)
[Z(ξ)]

m10(ξ) E[Z(ξ)] Êm10(ξ) [Z(ξ)], ξ ∈ Ξ10 Êm10(ξ) [Z(ξ)], ξ / ∈ Ξ10 (b) n = 10
[Z(ξ)]

m21(ξ) E[Z(ξ)] Êm21(ξ) [Z(ξ)], ξ ∈ Ξ21 Êm21(ξ) [Z(ξ)], ξ / ∈ Ξ21 (c) n = 21
[Z(ξ)] Behavior of Algorithm 3.2. Now, we illustrate the behavior of Algorithm 3.2 on Figure 3.2 and show the evolution with n of Ξ n and m n (ξ) for a single run of Algorithm 3.2, where m n (ξ) denotes the total number of samples from Z(ξ) generated from iteration 1 to iteration n. When n = 1, the algorithm has sampled every random variable once, which is not enough to distinguish some confidence intervals. So Ξ 1 is equal to Ξ. When n = 10, some confidence intervals can be distinguished and the algorithm has identified two groups of values where a quasi-maximum could be. These two groups correspond to the two groups of random variables in Ξ 10 .

m214(ξ) E[Z(ξ)] Êm214(ξ) [Z(ξ)], ξ ∈ Ξ214 Êm214(ξ) [Z(ξ)], ξ / ∈ Ξ214 (d) n = 214
When n = 21, the algorithm has identified the main peak of the function. However, the values of ξ,m(ξ) for ξ in Ξ 21 are not small enough for the algorithm to stop. Then the algorithm continues sampling the random variables in Ξ n , updating Ξ n when it is necessary. ξ,m(ξ) for ξ in Ξ n decreases since m(ξ) is increasing for these values of ξ and the algorithm stops at n = 214 when max

ξ∈Ξ 211 ξ,m(ξ) < τ 2+τ .
Figure 3.3 shows the influence of τ and λ on the average complexity E(M) of Algorithm 3.2. We observe that τ has a much bigger impact than λ. This observation is consistant with the impact of = τ /(2 + τ ) and δ = λ/#Ξ on the expected number of sampling E(M ) to build an estimate ÊM [Z] of E[Z(ξ)] with relative precision with probability 1 -δ (see Remark 3.5).

3.A Intermediate Results

Here we provide intermediate results used thereafter for the proof of Proposition 3.4 in section 3.B. We first recall a version of Bennett's inequality from [1, Lemma 5].

Lemma 3.11. Let U be a random variable defined on (Ω, F, P) such that U ≤ b almost surely, with b ∈ R. Let U 1 , . . . , U m be i.i.d. copies of U and U = 1 i=1 U i . For any x > 0, it holds, with probability at least 1 -exp(-x), simultaneously for all

1 ≤ ≤ m (3.36) U -E [U ] ≤ 2mE [U 2 ] x + b + x/3, with b + = max(0, b).
Now, the following result provides a bound with high probability for the estimated variance of an i.i.d. sequence of bounded random variables.

Lemma 3.12. Let X be a bounded random variable defined on (Ω, F, P), such that a ≤ X ≤ b almost surely, with a < b two real numbers. Let X 1 , . . . , X m be i.i.d.

copies of X and V m = 1 m m i=1 (X i -X m ) 2 where X m = 1 m m i=1 X i . Then, for any x > 0 (3.37) P   V m ≤ V[X] + 2V[X] (b -a) 2 x m + x(b -a) 2 3m   ≥ 1 -exp(-x).
Proof. Let us define U = (X -E[X]) 2 which satisfies U ≤ (b -a) 2 almost surely. Applying Lemma 3.11 with U defined previously with = m gives for any x > 0

P m U m -E[U ] ≤ 2mE[U 2 ]x + x(b -a) 2 3 ≥ 1 -exp(-x).
Moreover, as 2 and using the boundedness of U we get

U m = V m + (X m -E[X])
P   V m ≤ E[U ] + 2E[U ] (b -a) 2 x m + x(b -a) 2 3m   ≥ 1 -exp(-x),
which ends the proof since

E[U ] = V[X].
We recall a second result in the line of [START_REF] Mnih | Efficient stopping rules[END_REF]Lemma 3].

Lemma 3.13. Let q, k be positive real numbers. If t > 0 is a solution of

(3.38) log qt t = k, then (3.39) t ≤ 2 k log 2q k .
Moreover, if t is such that

(3.40) t ≥ 2 k log 2q k , then (3.41) log qt t ≤ k.
Proof. Let t > 0 be a solution of (3.38). Since the function log is concave, it holds for all s > 0 kt = log(qt) ≤ log(qs) + t -s s .

In particular, for s = 2 k > 0 we get

(3.42) t ≤ 2 k log 2q k -1 ≤ 2 k log 2q k ,
which yields (3.39). Now, let ϕ : s → log(qs) s defined for s > 0. This function is continuous, strictly increasing on (0, e q ] and strictly decreasing on [ e q , ∞) so it admits a maximum at t = e q . The existence of a solution t > 0 of (3.38) implies k ≤ q e . If k = q e then t = e q and ϕ(t) is the maximum of ϕ. For any t > 0, in particular satisfying (3.40), we have ϕ(t ) ≤ ϕ(t) = k which is (3.41). If 0 < k < q e , there are two solutions t 1 , t 2 to (3.38) such that 0 < t 1 < e q < t 2 . By (3.39) and (3.40) we have t ≥ t 2 > e q and since ϕ is stricly discreasing on [ e q , ∞) it holds ϕ(t ) ≤ ϕ(t 2 ) = k , that is (3.41).

3.B Proof of Proposition 3.4

Let us define the two events A = m≥1 A m and B = m≥1 B m with

A m = V m ≤ σ 2 + 2σ 2 (b -a) 2 log(3/d m )/m + log(3/d m )(b -a) 2 /
d m ≤ σ 2 ,
where

M σ 2 = 2(b -a) 2 σ 2 p log 2p(b -a) 2 σ 2 + log 3 cδ .
Again, Lemma 3.13 with k =

d m ≤ 2 µ 2 ,
where

M 2 µ 2 = 2(b -a)
V m ≤ (1 + √ 2 + 1/3) max(σ 2 , 2 µ 2 ).
In what follows, we define M = min (M σ 2 , M 2 µ 2 ). Now, we deduce from (3. .

If M ≥ M + then M ≥ M + .
Otherwise, M < M + and we have

M + = min m ≥ M : c m ≤ |µ| 1 +
. Moreover, as (3.49) holds for all m ≥ M , we get the inclusion difficulties when it comes to approximation. We have seen in Chapter 2 that approximation of parameter-independent high-dimensional problems is challenging. When it comes to high-dimensional parameter-dependent problems, we can imagine that other problems arise.

   m ∈ N * : m ≥ M , α log(3/d m ) m ≤ |µ| 1 +    ⊂ m ∈ N * : m ≥ M , c m ≤ |µ| 1 + .
We focus here on a particular context to introduce these solutions. It is usually referred as the multi-query context, see e.g. [START_REF] Haasdonk | Reduced basis methods for parametrized pdes-a tutorial introduction for stationary and instationary problems[END_REF], where many snapshots, i.e. u(y) for many different values of y in Y , are required. In such a context, one cannot only approximate the solution, e.g. using an algorithm such as Algorithm 2.3, as many times as the number of required snapshots: it would be too costly. One has to try to use already computed snapshots in order to reduce computational cost. This is what reduced basis methods (RBM) and more generally model order reduction (MOR) methods do. Our particular focus is on RBM that build a reduced space H n using snapshots u(y 1 ), . . . , u(y n ) (offline phase) and then project at low computational cost any desired snapshot u(y) onto H n (online phase). This offline-online decomposition is particularly relevant in many applications, see e.g. [START_REF] Salmoiraghi | Advances in geometrical parametrization and reduced order models and methods for computational fluid dynamics problems in applied sciences and engineering: overview and perspectives[END_REF], where calculus has to be made in real-time with a possible pre-computation of some needed quantities (here the snapshots). Of course, such techniques are particularly relevant when the solution manifold M := {u(y), y ∈ Y } can be well-approximated by lowdimensional subspaces of H.

We discuss theoretical aspects for reduced basis methods in Section 4.1, recalling an existing ideal algorithm to build a space H n with good approximation properties. Then in Section 4.2 we discuss practical implementation of the ideal algorithm. We finish our survey by considering in Section 4.3 the particular case of parameterdependent PDEs.

Idealized RBM

We present here basic notions related to RBM from a theoretical point of view. Then we recall a first algorithm to build reduced spaces together with first approximation results.

Solution Manifold and Kolmogorov Width

As we mentioned earlier, RBM are particularly relevant for problems where the solution manifold M defined below

M = {u(y) : y ∈ Y } ⊂ H, (4.2)
can be well-approximated by low-dimensional spaces H n .

We assume here that M is compact. Compactness can be obtained under some assumptions on Y and P, see e.g. [START_REF] Cohen | Kolmogorov widths under holomorphic mappings[END_REF]. Since we will focus on linear approximation tools, we recall the notion of Kolmogorov n-width for a compact set M noted d n (M) H and defined by

d n (M) H := inf dim(W )=n dist(M, W ) H , (4.3)
where dist(M, W ) H := sup u∈M d(u, W ) H and d(u, W ) H := inf w∈W u -w H . Other widths [START_REF] Pinkus | N-widths in Approximation Theory[END_REF] can be considered when it comes to non-linear approximation, which is not our concern here, but for which we refer to [START_REF] Devore | Nonlinear approximation[END_REF].

Problems for which it can be interesting to use RBM are thus the ones where the Kolmogorov n-width of M decreases quickly with n [START_REF] Pinkus | N-widths in Approximation Theory[END_REF]. However, these theoretical decays can be impacted by the dimension d, which can cause troubles when d 1, even for smooth functions [START_REF] Novak | Approximation of infinitely differentiable multivariate functions is intractable[END_REF]. As an example of the impact of the dimension, we recall a result from [START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF], that gives an estimate of the n-width of the unit ball K of

C m ([-1, 1] d ) in L ∞ -norm c d n -m/d ≤ d n (K) L ∞ ≤ C d n -m/d
where c d and C d are constants that grow exponentially with the dimension d. We thus observe that the dimension d impacts negatively the decay of the width, a phenomenon known as the curse of dimensionality, while the regularity improves this decay, phenomenon known as the blessing of regularity. However, for some other classes of more structured functions it is possible to circumvent the curse of dimensionality [START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF].

Greedy Algorithm

We give a first ideal algorithm for building reduced spaces. To that end, we first recall that in a Hilbert space we have inf w∈W u -w H = u -P W u H (4.4) where P W stands for the H-orthogonal projection onto W . We can now describe the greedy algorithm for RBM in Algorithm 4.1. We fix a finite number of steps n > 0 for this algorithm. Other stopping criteria exist, such as stopping when the desired approximation precision is reached [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]. Algorithm 4.1 Greedy Algorithm for Reduced Basis Methods [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF] 1: Set H 0 = {0}. 2: for i = 1 . . . n do 3:

Select y i ∈ Y such that u(y i ) -P H i-1 u(y i ) H ≥ γ sup y∈Y u(y) -P H i-1 u(y) H . (4.5) 4:
Compute u(y i ) and update H i = span{u(y 1 ), . . . , u(y i )}.

5: end for

Readers should notice that when γ = 1, Algorithm 4.1 corresponds to the strong greedy algorithm described for example in [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF]. Otherwise, when 0 < γ < 1, Algorithm 4.1 corresponds to a γ-weak-greedy algorithm such as the one described in [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF]. We introduce the following notation

σ n (M) H = dist(M, H n ) H , (4.6)
where H n has been constructed with Algorithm 4.1. We distinguish different results for the study of σ n (M) H . First of all, indirect comparison results assume a particular kind of upper bound for d n (M) H in order to derive an upper bound for σ n (M) H . For a γ-weak-greedy algorithm when H is a Hilbert space, we have the following existing result. Theorem 4.1. [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]Corollary 3.3] Assume that M is a compact set included in the unit ball of H. Assume moreover that, for all n ≥ 1,

d n (M) H ≤ C exp(-cn α )
with α > 0. Then, for all n ≥ 1,

σ n (M) H ≤ C 1 exp(-c 1 n α ), C 1 = √ 2Cγ -1 and c 1 = 2 -1-2α c.
Second of all, direct comparisons compare directly d n (M) H and σ n (M) H for any compact manifold M. Those direct comparison results are usually more pessimistic than indirect comparisons results, due to less strong assumptions. Still in the case of a γ-weak-greedy algorithm when H is a Hilbert space, we have the following result.

Theorem 4.2. [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]Corollary 3.3] For any compact set M included in the unit ball of H and for any n ≥ 1, we have

σ n (M) H ≤ √ 2γ -1 min 1≤m≤n (d m (M) H ) n-m n . (4.7)
It gives in particular, for all n ≥ 1

σ 2n (M) H ≤ √ 2γ -1 d n (M) H . (4.8)

Snapshot-Based Reduced Spaces

Algorithm 4.1 builds a reduced basis using snapshots u(y) in M, that are usually ortho-normalized for practical reasons. However, nothing indicates that optimal subspace for the Kolmogorov n-width can be generated by snapshots. We recall a result here from [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF] showing a relatively small loss when building of the reduced spaces using snapshots. Moreover, for any n > 1 and > 0, there exists a set M 1 such that

d n (M 1 ) H ≥ (n -1 -)d n (M 1 ) H .
When

Feasible Version of RBM

Algorithm 4.1 is impossible to implement in practice for several reasons. First, the exact solution u is in general inaccessible since it resides in an infinite-dimensional Hilbert space H. From this arise some limitations: 1) the impossibility to compute exactly a snapshot 2) the impossibility to provide an error estimate for u(y) -

P H i-1 u(y)
H that takes as a reference the exact solution u and thus the impossibility to solve the optimization problem max y∈Y u(y) -P H i-1 u(y) H . Other problems arise. Among others 3) the impossibility to compute an exact orthogonal projection

P H i-1 (

Truth Approximation

Since the exact solution u is inaccessible, we rather use an accurate approximation u N of u, called the truth approximation. This approximation lives in some highdimensional space H N , called the truth space, and can be computed by traditional approximation methods for the approximation of the solution of high-dimensional PDEs (see Chapter 2). We assume that it approximates u at some precision ∈ (0, 1), which means that for any y in Y we have u(y) -u N (y) H ≤ . (4.9) Two different possibilities exist for the space H N . Either the space H N ⊂ H of dimension N is fixed and large, i.e. N 1. Either the space H N (y) ⊂ H depends on y, e.g. through a dependence of its dimension N = N (y). In the first case H N ensures that, for any y in Y , the solution u(y) in H can be well represented by elements in H N . An example of such space is a finite element approximation space with N degrees of freedom. In the second case, the dependence on y allows important computational savings when it comes to the computation of the truth approximation of some snapshot up to some given precision , see (4.9). While more flexible, the second case is more difficult to describe in a few lines. Thus, for the sake of clarity, we stick to the first case in our presentation.

In the rest of our presentation, we assume that the exact and approximated solution coincide, which means that we take 1 of the order of magnitude of the machine precision for example. It allows us to neglect the differences between the high-dimensional and high-fidelity truth approximation u N and the exact solution u. Considering the truth approximation rather than the exact solution is not problematic when 1. In this case we can prove in a similar way as we will prove Theorem 4.4 that the behavior of both Kolmogorov n-width if M and M N := {u N (y) : y ∈ Y } are similar up to a small threshold 2 .

Feasible Projection

The truth approximation answers problem 1). In order to answer 3), we avoid using the exact orthogonal projection P H i-1 for the benefit of an approximated projection u i-1 := u H i-1 of u, that is easy to compute and that presents interesting properties, such as a quasi-optimality property. Quasi-optimality means that there exists a constant C > 0, such that for any y in Y we have

u(y) -u i-1 (y) H ≤ C u(y) -P H i-1 u(y) H . (4.10)
The independence of C on y or n and the existence of such projection are important points that we will discuss in Section 4.2.4 for parameter-dependent PDEs.

Discrete Training Set

In order to answer 4), we have to come up with an optimization strategy to solve max y∈Y u(y) -P H i-1 u(y) H . A commonly used technique is to discretize the solution manifold, which means considering a finite discrete set Ξ ⊂ Y and solve the optimization problem over this so-called training set. Thus the optimization problem becomes a discrete optimization problem which is possible to solve if we know how to compute u(y) -P H i-1 u(y) H for a finite number of values of y in Y which corresponds to Ξ. Then, denoting by M Ξ the following discrete solution manifold M Ξ = {u(ξ) : ξ ∈ Ξ}, (4.11) we can derive the following result.

Practical Greedy Algorithm

Using the tools introduced in Section 4.2.1, we can now propose a realistic version of Algorithm 4.1 in the sense that it is possible to implement. This is a weak-greedy algorithm for (4.1) over Y which means that, at each step i of the algorithm, we add a snapshot u(y i ) such that y i is the training set and satisfies

u(y i ) -P H i-1 u(y i ) H ≥ γ sup y∈Y u(y) -P H i-1 u(y) H (4.12)
for some fixed γ ∈ (0, 1). An algorithm for which the parameter γ does not change from one step to another is called a γ-weak-greedy algorithm over Y .

In order to make Algorithm 4.2 entirely practical, and thus to tackle problem 3), we assume that we have access to some error estimate ∆(u i-1 (y)) for u(y) -u i-1 (y) H for all y in Y , which is necessary since the H-norm might be hard to compute exactly. We moreover assume that this estimate satisfies

c 1 u(y) -u i-1 (y) H ≤ ∆(u i-1 (y)) ≤ C 1 u(y) -u i-1 (y) H , (4.13)
for some constants c 1 and C 1 independent from y and i. Compute u N (ξ i ) and update H i = span{u N (ξ 1 ), . . . , u N (ξ i )}.

5: end for

The following theorem shows that Algorithm 4.2 is a weak-greedy algorithm under suitable assumptions.

u(ξ i ) -P H i-1 u(ξ i ) H ≥ γ max ξ∈Ξ u(ξ) -P H i-1 u(ξ) H , (4.15)
where γ = c 1 C 1 C . Thus Algorithm 4.2 is a γ-weak-greedy algorithm for (4.1) over Ξ.

Proof. We take 1 ≤ i ≤ n. We have using (4.14) and(4.13)

∆(u i-1 (ξ i )) = max ξ∈Ξ ∆(u i-1 (ξ)) ≥ c 1 max ξ∈Ξ u(ξ) -u i-1 (ξ) H ≥ c 1 max ξ∈Ξ u(ξ) -P H i-1 u(ξ) H .
We have also from (4.10) and (4.13)

∆(u i-1 (ξ i )) ≤ C 1 u(ξ i ) -u i-1 (ξ i ) H ≤ C 1 C u(ξ i ) -P H i-1 u(ξ i ) H which concludes the proof.
Using this result and if we further satisfy

max ξ∈Ξ u(ξ) -P H i-1 u(ξ) H ≥ γ 1 max y∈Y u(y) -P H i-1 u(y) H , (4.16)
we obtain a γ 1 γ-weak-greedy algorithm over Y . In this case, all results seen in Section 4.1.2 apply to study the consequences of approximating M using the reduced space H n built by Algorithm 4.2.

RBM for Parameter-Dependent PDEs

We consider the particular case of parameter-dependent PDEs. In order to directly consider practical methods, we consider the discrete solution manifold defined in (4.11) and work with Ξ rather than Y . Such parameter-dependent problems are particularly interesting since, for several types of parameter-dependent PDEs, it is possible to circumvent the curse of dimensionality mentioned before [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric pdes[END_REF][START_REF] Cohen | Kolmogorov widths under holomorphic mappings[END_REF].

Weak-Formulation

We consider the following parameter-dependent weak-formulation of a PDE: For any ξ in Ξ, the goal is to find u(ξ) in H such that a(u(ξ), v; ξ) = f (v; ξ), ∀v ∈ H, (4.17) with H some function space. For all ξ in Ξ, a(•, •; ξ) is a bilinear form while f (•; ξ) is a linear form. We assume the existence of positive constants α

(ξ) ≤ α max < ∞ such that a(w, v; ξ) ≤ α(ξ) w H v H , ∀(w, v, ξ) ∈ H × H × Ξ. (4.18)
We also assume the existence of constants β(ξ) ≥ β min > 0 such that, for all ξ in Ξ

a(v, v; ξ) ≥ β(ξ) v 2 H ≥ β min v 2 H . (4.19)
In this context, the traditional truth space consists in Galerkin projections of u(ξ) onto some high-dimensional space H N . The Galerkin projection of u(ξ) onto a space H N , denoted by u N (ξ) is defined for each ξ in Ξ as the unique solution of

a(u N (ξ), v; ξ) = f (v; ξ), ∀v ∈ H N . (4.20)
Galerkin projection has the following property (Céa's lemma)

u(ξ) -u N (ξ) H ≤ α max β min u(ξ) -P H N u(ξ) H . (4.21)
Having a basis for H N , we can solve (4.20) for a fixed value of ξ by solving a highdimensional linear system (see e.g. discrete reduced basis in [START_REF] Haasdonk | Reduced basis methods for parametrized pdes-a tutorial introduction for stationary and instationary problems[END_REF]). We show now how to use the Galerkin projection to propose an efficient offline-online decomposition.

Affine Decomposition and Empirical Interpolation Method

We assume now that the forms a and f admit an affine decomposition of the following type

a(w, v; ξ) = Qa i=1 Θ a i (ξ)a i (w, v), (4.22) f (v; ξ) = Q f i=1 Θ f i (ξ)f i (v), (4.23)
where Θ a i , Θ f i are parameter-dependent functions, a i , f i are parameter-independent forms and Q a , Q f are strictly positive integers.

Once the snapshots (u N (ξ j )) n i=1 are computed, it is possible to compute the Galerkin projection u N n (ξ) of u(ξ) onto H n = span u N (ξ 1 ), . . . , u N (ξ n ) at low cost by precomputing the matrices A j = [a j (u N (ξ i ), u N (ξ k ))] j,k for all j, the vectors f j = [f j (u N (ξ k ))] n k=1 . Using (4.22) and (4.23), we end up with a linear system

A n (ξ)u n (ξ) = f n (ξ), where A n (ξ) = Qa i=1 Θ a i (ξ)A i and f n (ξ) = Q f i=1 Θ f i (ξ)f i ,
where the vector u n (ξ) contains the coefficients of the Galerkin projection onto H n expressed in the basis of the snapshots (u N (ξ j )) n i=1 . Since the system size does not depend on N , the computational cost for solving these systems (in O(n κ )) is much smaller than the computational cost to obtain a truth snapshot (in O(N κ ) for some κ > 1). It justifies the interest of such methods when α and β are of the same order of magnitude. The values of α and β depend on the nature (sparse, tridiagonal . . . ) of the linear system to solve.

When the affine decomposition of a or f is not obvious, the empirical interpolation method (EIM) [START_REF] Barrault | An 'empirical interpolation'method: application to efficient reduced-basis discretization of partial differential equations[END_REF] provides such a decomposition. It allows then an efficient offlineonline decomposition as we have seen before.

Error Estimation

There is still one missing piece to propose a practical version of Algorithm 4.1 for parameter-dependent PDEs: the estimation of the error u(ξ) -P H i-1 u(ξ) H , where H i-1 is the already built reduced space. A common approach in the literature is to propose an estimator based on the residual. We define the residual for all v ∈ H as

r i-1 (v; ξ) := a(u(ξ) -u i-1 (ξ), v; ξ). We obtain that r i-1 (v; ξ) = f (v; ξ) -a(u i-1 (ξ), v; ξ), ∀v ∈ H.
Then we have the following inequalities

β(ξ) e i-1 (ξ) H ≤ r i-1 (•; ξ) H (O) ≤ α(ξ) e i-1 (ξ) H , (4.24) where r i-1 (•; ξ) H (O) := sup v∈H r i-1 (v;ξ) v H and e i-1 (ξ) = u(ξ) -u i-1 (ξ). The quantity ∆(u i-1 (ξ)) := r i-1 (•; ξ) H (O) β min (4.25)
is an upper bound of the true error e i-1 (ξ) H and is usually taken as an error estimate. With this choice we have the following result. Theorem 4.6. We assume that (4.18), (4. [START_REF] Bungartz | Sparse grids[END_REF]) and (4.24) hold and that the error estimate is as in (4.25). Then Algorithm 4.2 adapted to parameter-dependent PDEs is a γ-weak-greedy algorithm over any training set Ξ with γ = (β min ) 2 (α max ) 2 . The use of an error estimate based on the residual degrades the performance of the practical weak-greedy algorithm for parameter-dependent PDEs. It would be more interesting to directly deal with the approximation error u(ξ) -u i-1 (ξ) H , particularly for problems with small β min . We present a new technique to avoid using residual-based error estimate in Chapter 5.

Our problems of interest are parameter-dependent PDEs whose solution u(y) admits a probabilistic representation, which means that there exists some parameterdependent functional Q y : R d → R, depending on the data from (5.1), and some random variable X y defined on a probability space (Ω, F, P) such that

u(y) = E [Q y (X y )] . (5.2)
This probabilistic representation of the solution is motivated by the Feynman-Kac representation of solutions of partial differential equations (PDEs). We give in Section 5.1 two guiding examples of such representations for solutions of elliptic and parabolic parameter-dependent PDEs. Those are extensions to parameterdependent PDEs of the probabilistic representations given in Chapter 1. In Section 5.2 we use this probabilistic representation of the solution to derive a probabilistic representation of the approximation error u(y) -u n (y) for any y in Y , where u n : Y → H n is an approximation of u onto the finite dimensional (reduced) space H n ⊂ H. In particular we will reinterpret the H-norm of the approximation error as an expectation. It allows to propose a first probabilistic greedy algorithm to construct the reduced space H n that is analyzed in Section 5.3.1. Then our second probabilistic greedy algorithm combines the first algorithm with random training sets from [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] which allows us to treat parameter-dependent problems with highdimensional parametric dimension.

Probabilistic Representation for Parameter-Dependent PDEs

We extend the Feynman-Kac representation from Chapter 1 to parameter-dependent PDEs.

Elliptic Parameter-Dependent PDEs

We consider an elliptic parameter-dependent PDE on a open bounded domain D ⊂ R d . For any y in Y , let -L y u(x; y) = g(x; y) x ∈ D, u(x; y) = f (x; y) x ∈ ∂D, (5.3) where the partial differential operator is defined as follow

L y = 1 2 d i,j=1 (σ(x; y)σ(x; y) T ) ij ∂ 2 x i x j + d i=1 b i (x; y)∂ x i while g : D × Y → R
= x ∈ D. (5.4)
For a fixed value of y in Y , we can then state a similar result as Theorem 1.12 that applies to (5.3).

Υ i Υ j ≥ c d i=1 Υ 2 i , Υ ∈ R d , x ∈ D.
Then, for all y in Y , there exists a unique solution u(•; y) of (5.3) in C(D) ∩ C 2 (D).

It satisfies for all x ∈ D u(x; y) = E f (X x;y τ x;y ; y) + τ x;y 0 g(X x;y t ; y)dt . (5.5)

Parabolic Parameter-Dependent PDEs

We consider a parabolic parameter-dependent PDE. For any y in Y (5.6) where A y = ∂ t + L y , with the differential operator L y is defined as follow

A y u(t, x; y) + g(t, x; y) = 0 (t, x) ∈ [0, T ) × D, u(T, x; y) = f (T, x; y) x ∈ D, u(t, x; y) = f (t, x; y) (t, x) ∈ [0, T ) × ∂D,
L y = 1 2 d i,j=1 (σ(t, x; y)σ(t, x; y) T ) ij ∂ 2 x i x j + d i=1 b i (t, x; y)∂ x i , while g : [0, T ] × D × Y → R is the source term and f : [0, T ] × D × Y → R represents
both the terminal and boundary condition. In the parabolic case, L y can also been seen as the infinitesimal generator of the stochastic process X ; y)ds (5.9) where τ t,x;y T = τ t,x;y ∧ T .

Error Estimation using Probabilistic Representation

In the context of RBM, a greedy algorithm consists, at each step i + 1 ≥ 1, in seeking y that maximizes the approximation error to enrich the current reduced space H i ⊂ H. This optimization problem usually takes the following form sup y∈Y u(y) -u i (y) H , (5.10) where u i (y) : O → R is some approximation of u(y) in some space H i (e.g. obtained by interpolation, least-squares or Galerkin projection). We assume all along this section that both the approximation u i and the space H i are deterministic. In this section, we propose to interpret the squared norm of the approximation error u(y) -u i (y) 2 H as an expectation using probabilistic representations introduced in Section 5.1. We assume here that, for any y in Y , we have u i (y) = i j=1 α j (y)φ j , (5.11) where the family of functions (φ j ) j=1...i is a basis of H i , obtained from snapshots u(y j ) selected for example by a greedy algorithm.

Parameter-Dependent Elliptic PDEs

Assuming that u i (y) is regular enough (see Theorem 5.3) for all y in Y (through the choice of the reduced basis functions φ j ), we have that e i (y) := u(y) -u i (y) is the unique solution, for all y in Y , of -L y e i (y) = g i (y) on D, e i (y) = f i (y) on ∂D, (5.12) where f i (y) := f (y) -u i (y) and g i (y) = g(y) + L y u i (y). Applying Theorem 5.1, we obtain the following probabilistic representation for e i (y) for any y in Y : Then, for all y in Y , e i (y) is the unique solution of (5.12) in C(D) ∩ C 2 (D) and satisfies for all x ∈ D e i (x; y) = E [F i (x, X x;y ; y)] := E f i (X x;y τ x;y ; y) + τ x;y 0 g i (X x;y t ; y)dt .

(5. [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF] In what follows we suppose that H ⊂ L 2 (D) is equipped with the L 2 -norm over D.

Using the probabilistic representation of e i (y), we prove the following reinterpretation of e i (y) 2 L 2 .

Theorem 5.4 (Reinterpretation of e i (y) 2 L 2 ). Taking V ∼ U(D) uniformly distributed on D, we have for any y in Y e i (y) 2 L 2 = |D|E F i (V, X V ;y ; y)F i (V, XV ;y ; y) , (5.14) where X x;y and Xx;y are two i.i.d random processes for any x ∈ D and |D| is the Lebesgue measure of D.

Proof. We first recall

e i (y) 2 L 2 = D (e i (x; y)) 2 dx = |D|E (e i (V ; y)) 2 .
Then, by the tower property

E (e i (V ; y)) 2 = E E F i (V, X V ;y ; y))|V 2 = E E F i (V, X V ;y ; y)|V E F i (V, XV ;y ; y)|V .
Then, by independence of X V ;y |V and XV ;y |V we can write

E (e i (V ; y)) 2 = E E F i (V, X V ;y ; y)F i (V, XV ;y ; y)|V .
Finally, by the law of iterated expectation, we get (5.14

Parameter-Dependent Parabolic PDEs

Here the functions u(y), u i (y), f (y), g(y) and e i (y) depend on t and x, which we will not precise when there is no ambiguity. In this section, the L 2 -norm is taken over [0, T ] × D.

Assuming that u i (y) is regular enough (see Theorem 5.6) for all y in Y (through the choice of the reduced basis functions φ i ), we have that the approximation error e i (y) := u(y) -u i (y) satisfies, for all y in Y , the following equations

A y e i (y) + g i (y) = 0 on [0, T ) × D, e i (y) = f i (y) on {T } × D, e i (y) = f i (y) on [0, T ) × ∂D, ( 5.15)
where f i (y) := f (y) -u i (y) and g i (y) = g(y) + A y u i (y). Applying Theorem 5.2, we obtain the following probabilistic representation for the unknown e i (y) = u(y)u i (y). Then, for all y in Y , e i (•, y) admits the following probabilistic representation e i (t, x; y) = E[F i (t, x, X t,x;y ; y)] := E f i (τ x;y , X t,x;y τ x;y ; y) + τ x;y t g i (s, X t,x;y s ; y)ds . (5.16) Using the above probabilistic representation, we prove the following expression of the squared L 2 -norm of the approximation error, following the proof of Theorem 5.4.

e i (y) 2 L 2 = T |D|E F i (T, V, X T,V ;y ; y)F i (T, V, XT,V ;y ; y) ,
where X t,x;y and Xt,x;y are two i.i.d random processes for any (t, x) ∈ [0, T ] × D.

Probabilistic Greedy Algorithms

We present in Section 5.3.1 a first greedy algorithm using probabilistic error estimate. Then we recall in Section 5.3.2 another greedy algorithm from [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] using random training sets. Then these algorithms are combined in Algorithm 5.3.

A similar contribution can be found in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF], where the authors proposed a variance reduction algorithm using the reduced basis paradigm for parameter-dependent stochastic differential equations. The analysis of this algorithm can be found in [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF].

One has to be careful in this section, since neither the approximation u i nor the reduced space H i are deterministic. They depend on the previous (probabilistic) steps of the greedy algorithms. We have thus to make a distinction between what is simulated before stage i and at stage i in order to analyse these algorithms. We will define appropriate conditional expectations before presenting each one of them.

In this section, we consider the L 2 -norm in order to reuse the probabilistic interpretation of the error proposed in Section 5.2.

With Probabilistic Error Estimate

We describe a first probabilistic greedy algorithm that selects randomly the snapshot to enrich the reduced space using the probabilistic representation of the error introduced in Section 5.2. Select ξ i using PAC algorithm to solve max ξ∈Ξ E[Z i (ξ)] with relative precision and probability greater than 1 -λ i .

4:

Compute u(ξ i ) and update H i = span{u(ξ 1 ), . . . , u(ξ i )}.

5: end for Theorem 5.10. Take (λ i ) 1≤i≤n ∈ (0, 1) such that n i=1 λ i < 1, ∈ (0, 1) and Ξ a discrete training set. Assume moreover that for all i = 1, . . . , n u(y) -u i-1 (y) L 2 ≤ C u(y) -P H i-1 u(y) L 2 , (5. [START_REF] Bungartz | Sparse grids[END_REF] for some constant C independent from H i-1 and y. Then, with probability at least 1 -n i=1 λ i , Algorithm 5.1 is a weak-greedy algorithm of parameter Proof. In order to analyse Algorithm 5.1, we define here the following conditional expectation, for any i ≥ 1 and any random variable on (Ω, F, P) (5.21) where Z <i is the random variables Z j,k (ξ) for ξ in Ξ and j < i, where Z j,k (ξ) are i.i.d. copies of Z j (ξ) generated by PAC algorithm at step j of the greedy algorithm. Then for any element A in F, we let

E i [•] = E •|Z <i ,
P i (A) = E i [1 A ] . (5.22)
Using these notations we can write, in the case where u i is randomly determined by the steps before i, that

e i (y) 2 L 2 = u(y) -u i (y) 2 L 2 = E i+1 [Z i+1 (ξ)] . (5.23)
These notations being introduced, we define the following event A = ∩ n i=1 A i , each A i being defined as

A i := E i [Z i (ξ i, )] -E i [Z i (ξ i )] ≤ E i [Z i (ξ i, )] .
We have P i (A i ) ≥ 1 -λ i almost surely. Now, for all ω ∈ Ω such that A is verified, we have for any i ≥ 1

u(ξ i ) -P H i-1 u(ξ i ) 2 L 2 ≥ 1 C 2 u(ξ i ) -u i-1 (ξ i ) 2 L 2 = 1 C 2 E i [Z i (ξ i )] ≥ 1 - C 2 E i [Z i (ξ i, )] = 1 - C 2 max ξ∈Ξ E i [Z i (ξ)] = 1 - C 2 max ξ∈Ξ u(ξ) -u i-1 (ξ) 2 L 2 ≥ 1 - C 2 max ξ∈Ξ u(ξ) -P H i-1 u(ξ) 2 L 2
Taking the square root of the obtained inequality, we obtain (5.20). We now estimate

P(A)

P(A) = 1 -P(A) ≥ 1 -

n i=1 P(A i ) = 1 - n i=1 E 1 A i = 1 - n i=1 E[E 1 A i |Z <i P i (A i ) ] ≥ 1 - n i=1 λ i ,
where the last inequality derives from the selection procedure of ξ i , see step 3 of Algorithm 5.1. It gives us the desired result.

With Random Training Sets

We consider now a Hilbert space H. We have not mentioned how we choose Ξ so far. For low-dimensional parameter spaces Y ⊂ R p , i.e. when p is relatively small, it is possible to take for Ξ some -net of Y with a small . However when p is too high, this choice for Ξ is not reasonable since the number of parameters in Ξ grows exponentially with p. The authors in [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] proposed a solution when p is large, by using random training sets. We recall in Algorithm 5.2 the greedy algorithm using these random training sets whose analysis in made in Theorem 5.11. Increment i ← i + 1.

4:

Select ξ i ∈ Ξ i such that

ξ i ∈ arg max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) H ,
and σi = max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) H .

5:

Compute u(ξ i ) and update H i = span{u(ξ 1 ), . . . , u(ξ i )}.

6: end while

Before presenting the analysis of this algorithm, we recall a few notations. We introduce the polynomial space P Λ := span{y → y ν := p j=1 y ν j j : ν ∈ Λ} where ν is a multi-index and Λ a downward closed set of multi-indices. Then Q Λ := H ⊗ P Λ is the set of functions from Y to H that can be written u(y) = ν∈Λ c ν y ν with c ν ∈ H. We also define Σ m := ∪ #Λ=m Q Λ the set of functions from Y to H that admits a m-term expansion on the monomial basis {y ν : ν ∈ N p }. With these notations, the class of functions A r := A r ((Σ m ) m≥1 ) is composed of all the functions v : Y → H in L ∞ (Y, H) such that δ m (u) := inf The norm • A r can be defined on A r by v A r := max{ v L ∞ (Y,H) , |v| A r }, where |v| A r is defined as the smallest C such that δ m (v) ≤ Cm -r . Theorem 5.11. [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]Theorem 2] Take η, in (0, 1). Assume that

• y → u(y) belongs to A r for some r > 2 ;

• u A r ≤ M 0 for some M 0 > 0 ;

• m is the smallest integer such that 32M 0 m -r+2 ≤ and 2 4r+2 m -r ≤ 1 ; • N is the smallest integer such that (1 - Then, with probability greater than 1-η, Algorithm 5.2 is a γ-weak-greedy algorithm with γ = 1/8m over Y , which means that at each step before the stopping of the algorithm u(ξ i ) -P H i-1 u(ξ i ) H ≥ γ sup y∈Y u(y) -P H i-1 u(y) H . (5.25) Moreover, denoting by H f the reduced space built at the end of Algorithm 5.2, we have dist(M, H f ) := sup y∈Y u(y) -P H f u(y) H ≤ with probability greater than 1 -η.

Proof. We only repeat here the proof of [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]Theorem 4.1] for the first part of the proof, with some details useful for the next results.

We define the following event B = ∩ f i=1 B i where f is when Algorithm 5.2 stops and B i is defined by

B i = {σ i ≥ γσ i (M) H } ,
where σ i (M) H = dist(M, H i-1 ) H . We recall σi = max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) H . We observe B ⊃ ∩ m 2 i=1 B i . Thus we can write

P(B) = 1 -P(B) ≥ 1 - m 2 i=1 P(B i ) = 1 - m 2 i=1 E 1 B i = 1 - m 2 i=1 E E 1 B i |Ξ 1 , . . . , Ξ i-1 .
For fixed training sets Ξ 1 , . . . , Ξ i-1 , H i-1 is fixed and [27, Lemma 3.1] gives E 1 B i |Ξ 1 , . . . , Ξ i-1 = P(B i |Ξ 1 , . . . , Ξ i-1 ) ≤ η m 2 , which concludes our proof. Algorithm 5.1 and Algorithm 5.2 can be combined to give Algorithm 5.3 that we analyze in Theorem 5.12. In order to do so, we take H = L 2 (O) and suppose that we consider only the L 2 -norm.

Algorithm 5.3 Full Probabilistic Greedy Algorithm

Require: (λ i ) 1≤i≤m 2 ∈ (0, 1), (Ξ i ) i=1...m 2 a collection of discrete training sets.

1: Set i = 0, H 0 = {0} and σ0 = 8m . 2: while σi ≥ 8m and i < m 2 do 3:

Increment i ← i + 1.

4:

Select ξ i using PAC algorithm to solve max ξ∈Ξ i E[Z i (ξ)] with relative precision and probability greater than 1 -λ i .

5:

Set σi = max ξ∈Ξ i E[Z i (ξ)].

6:

Compute u(ξ i ) and update H i = span{u(ξ 1 ), . . . , u(ξ i )}. Proof. In order to analyse Algorithm 5.3, we define here the following conditional expectation, for any i ≥ 1 and any random variable on (Ω, F, P) (5.26) where Z <i is the Ξ j for j < i and the random variables Z j,k (ξ) for ξ in Ξ j and j < i, where Z k,j (ξ) are i.i.d. copies of Z j (ξ) generated by PAC algorithm at step j of the greedy algorithm. Using these notations we can write, in the case where u i is randomly determined by the steps before i, that e i (y) 2

E i [•] = E •|Z <i ,
A i := E i [Z i (ξ i, )] -E i [Z i (ξ i )] ≤ E i [Z i (ξ i, )] , B i = {σ i ≥ γσ i (M) L 2 } ,
where σ i (M) H = dist(M, H i-1 ) H . We recall σi = max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) H . We first observe, that for all ω in Ω such that A ∩ B is verified, we have for any i u(ξ i ) -

P H i-1 u(ξ i ) L 2 ≥ √ 1 - C max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) L 2 = √ 1 - C σi ≥ √ 1 - 8mC σ i (M) L 2 = √ 1 - 8mC max y∈Y u(y) -P H i-1 u(y) L 2 .
Using A ⊃ ∩ m 2 i=1 A i and B ⊃ ∩ m 2 i=1 B i , we can write

P(A ∩ B) = 1 -P(A ∪ B) ≥ 1 - m 2 i=1 P(A i ) - m 2 i=1 P(B i ) = 1 - m 2 i=1 E 1 A i - m 2 i=1 E 1 B i = 1 - m 2 i=1 E E 1 A i |Z <i , Ξ i - m 2 i=1 E E 1 B i |Z <i .
Then we proceed as in the proofs of Theorem 5.10 and Theorem 5.11 to conclude the proof.

Remark 5.13. We keep the same stopping criterion as in Algorithm 5.2 for Algorithm 5.3. However this stopping criterion would not be available for a practical implementation of Algorithm 5.3 since we cannot compute or give any precise estimation of max ξ∈Ξ i u(ξ) -P H i-1 u(ξ) L 2 = max ξ∈Ξ i E[Z i (ξ)]. However, e.g. using the bounds of the confidence intervals provided for the estimates of the expectations by PAC algorithm, it is possible to give a practical stopping criterion that can be used for Algorithm 5.3 while ensuring the reaching of a given precision.

Conclusion and Outlooks

In Chapter 1 we have surveyed standard probabilistic methods for partial differential equations (PDEs), including the probabilistic representation for the solution of a PDE and Monte-Carlo methods coupled with time-integration techniques to get pointwise estimations of this solution. Then, in Chapter 2, we have proposed two algorithms to approximate the global solution of high-dimensional PDEs whose solution admits a probabilistic interpretation. These algorithms combine an adaptive sparse interpolation method from [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF] with a control variates technique from [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. Our first algorithm uses the adaptive sparse interpolation method from [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF] to build adaptively the control variate from [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF]. This first algorithm is not very stable and the approximation it builds is not precise compared to the one resulting from our second algorithm. This second algorithm is a perturbed version of the adaptive sparse interpolation method from [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric pdes[END_REF]. At each step of the algorithm, we compute the approximation resulting from the control variates technique from [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF] on the interpolation basis rather than the interpolation of the exact solution on this interpolation basis. Both algorithms were illustrated theoretically and numerically. We have in particular showed how the second algorithm does not only allow to get precise estimations of the solution on a fixed interpolation grid, such as it has been done in [START_REF] Gobet | Sequential control variates for functionals of markov processes[END_REF], but on the entire domain of definition of the solution. An approximation up to machine precision was built by the second algorithm for an elliptic PDE in dimension 20.

Other adaptive approximation techniques, such as adaptive least-squares, can be used instead of adaptive interpolation techniques. Together with a precise study of the time-integration error for some functional of a stochastic process, it might lead to more advanced results of the study of the convergence of our algorithms.

Then we have studied in Chapter 3 an optimization problem over a discrete set of random variables where the quantity to optimize is their expectation. We have combined expectation estimation in relative precision [START_REF] Mnih | Empirical bernstein stopping[END_REF] with a probably approximately correct algorithm in absolute precision [START_REF] Even-Dar | Pac bounds for multi-armed bandit and markov decision processes[END_REF] to propose two probably approximately correct algorithms in relative precision. The first one naively builds precise estimates for each expectation to finally select the random variable with the highest estimate. The second one adaptively samples each random variable and relies on confidence intervals for these estimates to decide which random variable to continue sampling and which random variable to stop sampling. Both algorithms provide the same theoretical guarantees for their output. However the second algorithm allows important savings in terms of number of samples used compared to the first In a first part, we consider partial differential equations in high dimension. Based on a probabilistic interpretation of the solution which allows to obtain pointwise evaluations of the solution using Monte-Carlo methods, we propose an algorithm combining an adaptive interpolation method and a variance reduction method to approximate the global solution.

In a second part, we focus on reduced basis methods for parametric partial differential equations. We propose two greedy algorithms based on a probabilistic interpretation of the error. We also propose a discrete optimization algorithm probably approximately correct in relative precision which allows us, for these two greedy algorithms, to judiciously select a snapshot to add to the reduced basis based on the probabilistic representation of the approximation error.

1 . 2

 12 Comparison bewteen Walk on Sphere (left) and Green Function First Passage (right) techniques [62]. . . . . . . . . . . . . . . . . . . . . . 2.1 (TC1) Algorithm 2.1 for fixed Λ : evolution of u -ũk Λ 6 with respect to k for various M (left figure), and various ∆t (right figure). . . . . 2.2 (TC1) Algorithm 2.1 for fixed Λ i : evolution of u-ũk Λ i 2 with respect to k for i = 8 (left figure), and i = 10 (right figure). . . . . . . . . . . 2.3 (TC1) Comparison of Algorithm 2.2 applied to exact solution and Algorithm 2.3 : (left) absolute error in L 2 -norm (right) evolution of ε n and εn with respect to #Λ n . . . . . . . . . . . . . . . . . . . . . . 2.4 (TC1) Adaptive Algorithm 2.1: evolution of u -u k Λ k 2 (continuous line) and u -u k Λ k ∞ (dashed line) with respect to step k and ∆t. . . 3.1 Final state of each algorithm after one run with τ = 0.1, λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME Algorithm. Left scale : values of the estimates Êm(ξ) [Z(ξ)] together with the associated confidence intervals of level 1 -d m(ξ) . Right scale : values of m(ξ). . . . . . . . . . . . . . 3.2 Evolution of Ξ n and number of samples m n (ξ) with n for Algorithm 3.2 with τ = λ = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Average complexity E(M) of Algorithm 3.2 with respect to τ and λ (in log-log scale). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 1

 1 Average complexity E(M), estimated using 30 runs for each algorithm, with τ = 0.1, λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME algorithm. 3.2 Estimated runtime T (in seconds) for different values of t , with τ = λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME algorithm. . . . . . . . . . . . .

(1. 1 )

 1 where b = (b 1 , . . . , b d ) : [0, T ] × D → R d and σ : [0, T ] × D → R d×d . Chapter 1. Probabilistic Methods for PDEs With these notations, the backward parabolic boundary-value problem is composed of the following equations satisfied by the unknown solution u

  2) is a parabolic PDE. We have k : [0, T ]×D → R, and f : [0, T ]×D → R represents both the terminal and boundary conditions (respectively (1.3) and (1.4)). The function g : [0, T ] × D → R is the source term of the PDE. In the literature, the set of equations (1.2)-(1.4) is also referred to as the Cauchy-Dirichlet problem. We study also elliptic problems, where the time variable does not appear. It is composed of the following equations satisfied by the unknown solution u : D → R -L(u) + ku = g in D, (1.5) u = f on ∂D. (1.6) Equation (1.6) stands for a Cauchy boundary condition. The functions k : D → R, f : ∂D → R and g : D → R are not time dependent in this case. Neither are the functions b : D → R d and σ : D → R d×d defining the coefficients of the linear differential operator (1.1).

Theorem 1 . 6 (

 16 Existence and uniqueness of a solution to a SDE). Let W be a ddimensional Brownian motion. Take s ∈ R + and s < S < +∞. Let b : [s, S] × R d → R d and σ : [s, S] × R d → R d×d be two continuous functions satisfying the following conditions for a constant C b,σ

  Ito formula for the product of the Ito process Y s and the function Z s defined by (where b Y and σ Y are scalar drift and diffusion coefficients of the SDE Chapter 1. Probabilistic Methods for PDEs satisfied by Y )

Figure 1 . 2 :

 12 Figure 1.2: Comparison between Walk on Sphere (left) and Green Function First Passage (right) techniques [62]. Both procedures work the same, until the selection of x 5 . When x 4 has been drawn, it lays at a distance to the boundary smaller than

. 1 )

 1 where u : D → R is a real-valued function, and D is an open bounded domain in R d . A is an elliptic linear differential operator and f : ∂D → R, g : D → R are respectively the boundary condition and the source term of the PDE.

  where W is a d-dimensional Brownian motion and b := (b 1 , . . . , b d ) T : R d → R d and σ : R d → R d×d stand for the drift and the diffusion respectively.

Theorem 2 . 1 (

 21 Feynman-Kac formula for elliptic PDE on bounded domain). Assume that(E1) D is an open connected bounded domain of R d , regular in the sense that, if τ x = inf {s > 0 : X x s /∈ D} is the first exit time of D for the process X x , solution of (2.3), we haveP(τ x = 0) = 1, x ∈ ∂D,(E2) b, σ are Lipschitz functions over D, (E3) f is continuous on ∂D, g and k ≥ 0 are Hölder-continuous functions on D, (E4) (uniform ellipticity assumption) there exists c > 0 such that d i,j=1

Algorithm 2 . 3 2 :

 232 Perturbed adaptive sparse interpolation algorithm 1: Set Λ 1 = {0 d } and n = 1. while n ≤ N and εn-1 > ε do 3:

5 •

 5 10 -3 ∆t = 6.25 • 10 -4

Figure 2 . 1 :

 21 Figure 2.1: (TC1) Algorithm 2.1 for fixed Λ : evolution of u -ũk Λ 6 with respect to k for various M (left figure), and various ∆t (right figure).

2 Figure 2 . 2 :

 222 Figure 2.2: (TC1) Algorithm 2.1 for fixed Λ i : evolution of u -ũk Λ i 2 with respect to k for i = 8 (left figure), and i = 10 (right figure).

Figure 2 . 3 :

 23 Figure 2.3: (TC1) Comparison of Algorithm 2.2 applied to exact solution and Algorithm 2.3 : (left) absolute error in L 2 -norm (right) evolution of ε n and εn with respect to #Λ n .

Figure 2 . 4 :

 24 Figure 2.4: (TC1) Adaptive Algorithm 2.1: evolution of u-u k Λ k 2 (continuous line) and u -u k Λ k ∞ (dashed line) with respect to step k and ∆t.

. 6 )

 6 Using Theorem 3.1, we see that c m stands for the half-length of a confidence interval of level 1 -d m for E[Z], i.e. (3.7) P(|Z m -E[Z]| ≤ c m ) ≥ 1 -d m .

(3. 16 )Proposition 3 . 4 .

 1634 d m = δcm -p , c = p -1 p , which satisfies (3.10) and (3.15), for any p > 1. The following result extends the result of [78, Theorem 2] stated for random variables Z with range in [0, 1]. Let 0 < δ ≤ 3/4 and let Z be a random variable taking values in a bounded interval [a, b], with expectation µ = E[Z] and variance σ 2 = V[Z]. If µ = 0 and (d m ) m≥1 satisfies (3.16), then M defined by (3.14) satisfies M < +∞ almost surely and (3.17) where • denotes the ceil function and

1 : 3 :

 13 Set = τ 2+τ and δ = λ/#Ξ. 2: for all ξ ∈ Ξ do Build an estimate Êm(ξ) [Z(ξ)] of E[Z(ξ)] using (3.20) with and δ as above.

4: end for 5 :Proposition 3 . 6 .

 536 Select ξ such that ξ ∈ arg max ξ∈Ξ Êm(ξ) [Z(ξ)] . Let (τ, λ) ∈ (0, 1) 2 . We assume that, for all ξ ∈ Ξ, Z(ξ) is a bounded random variable with E[Z(ξ)] = 0. Moreover we assume that the sequence (d m ) m≥1 is such that and log(3/d m )/m → m→+∞ 0.

2 :

 2 while #Ξ n > 1 and max

8: end while 9 :

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

  log(3/d m )/m → 0 as m → +∞ and E[Z(ξ)] = 0 for all ξ, we have that ξ,m a.s → m→+∞ 0 for all ξ in Ξ. Yet, since at each iteration n (from the subsequence),

2

 2

Figure 3 . 1 :

 31 Figure 3.1: Final state of each algorithm after one run with τ = 0.1, λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME Algorithm. Left scale : values of the estimates Êm(ξ) [Z(ξ)] together with the associated confidence intervals of level 1 -d m(ξ) . Right scale : values of m(ξ).

• 10 1 3 • 1 . 4 • 10 9 Table 3 . 2 :

 314932 10 -3 1.2 • 10 3 10 -6 2.0 • 10 1 1.9 • 10 2 3.8 • 10 -3 1.4 • 10 3 10 -4 2.0 • 10 3 1.9 • 10 4 3.8 • 10 -1 1.4 • 10 5 10 -2 2.0 • 10 5 1.9 • 10 6 3.8 • 10 1 1.4 • 10 7 1 2.0 • 10 7 1.9 • 10 8 3.8 • 10 3 Estimated runtime T (in seconds) for different values of t , with τ = λ = 0.1 and τ abs = τ |E[Z(ξ )]| for ME algorithm.

Figure 3 . 2 :

 32 Figure 3.2: Evolution of Ξ n and number of samples m n (ξ) with n for Algorithm 3.2 with τ = λ = 0.1.

5 τFigure 3 . 3 :

 533 Figure 3.3: Average complexity E(M) of Algorithm 3.2 with respect to τ and λ (in log-log scale).

3 , 3 3 δc 1 /p gives for any integer m ≥ M σ 2 ( 3

 333123 3m , and B m = |Z m -µ| ≤ c m . Applying Lemma 3.12 with x = log(3/d m ) for A m , m ≥ 1 together with a union bound argument leads to P(A) ≥ 1 -δ/3. Similarly, using a union bound argument and Theorem 3.1 with x = log(3/d m ), for B m , m ≥ 1, gives P(B) ≥ 1 -δ. By gathering these two results we have (3.43) P (A ∩ B) ≥ 1 -(P(A c ) + P(B c )) ≥ 1 -4δ where A c and B c correspond respectively to the complementary events of A and B. It remains to prove that A ∩ B implies (will prove (3.17). In what follows, we suppose that A ∩ B holds. First we derive an upper bound for V m . Since A holds, we have (3.45) V m ≤ σ 2 + 2σ 2 (b -a) 2 log(3/d m )/m + log(3/d m )(b -a) 2 /3m. Lemma 3.13 with k = σ 2 p(b-a) 2 and q =

2 µ 2 p

 22 (b-a) 2 and q = 3 δc 1/p gives for any integer m ≥ M 2 µ 2 (3.47) (b -a) 2 m log 3

2 ,

 2 48) an upper bound for c m . By definition, c m = 2V m log(3/d m ) m + 9(b -a) 2 log(3/d m ) 2 m then for all integer m ≥ M and using either (3.46), or (3.47), we have (3.49) c m ≤ α log(3/d m ) m , with α := ( 2 + 2 √ 2 + 2/3 + 3) 2 max(σ 2 , 2 µ 2 ). Now, using (3.49), we seek a bound for M , the smallest integer such that c M ≤ |Z M |. To that aim, let us introduce the integer M , valued random variable M + (3.51) M + = min m ∈ N * : c m ≤ |µ| 1 +

2 (

 2 Taking the min leads again to M ≥ M + . Moreover, since B holds, |µ|-c M + ≤ |Z M + | and using (3.51) it implies that c M + ≤ |Z M + |. By definition of M we get M + ≥ M . Hence, we have M ≥ M . To conclude the proof, it remains to find an upper bound for M . Applying again Lemma 3.13 with k = 2 µ 1+ ) 2 αp and q = 3 δc 1/p gives for any integer m ≥ M f (3.52) α log(3/d m ) m ≤ 2 µ 2 (1 + ) 2

Theorem 4 . 3 .

 43 [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF] Theorem 4.1] Let d n (M) H = inf W ∈Wn(M) sup y∈Y u(y) -P W u(y) H , where W n (M) is composed of all n-dimensional spaces spanned by n snapshots from M. For any compact set M and any n ≥ 0 d n (M) H ≤ (n + 1)d n (M) H .

Algorithm 4 . 2 3 :

 423 Practical Greedy Algorithm for RBM1: Set X 0 = {0}. 2: for i = 1 . . . n do Select ξ i ∈ Ξ such that ξ i ∈ arg max ξ∈Ξ ∆(u i-1 (ξ)). (4.14) 4:

Theorem 4 . 5 .

 45 Assume that Ξ is a training set and u = u N is the finite-dimensional solution to (4.1). Moreover assume that (4.13) and (4.10) hold. Then we have for all i = 1, . . . , n that

Theorem 5 . 1 (

 51 Feynman-Kac formula for an elliptic parameter-dependent PDE on bounded domain). Assume that, for all y in Y (EP1) D is an open connected bounded domain of R d , regular in the sense that, if τ x;y = inf {s > 0 : X x;y s / ∈ D} is the first exit time of D for the process X x;y , solution of (5.4), we have P(τ x;y = 0) = 1, x ∈ ∂D, (EP2) b(•; y), σ(•; y) are Lipschitz functions over D, (EP3) f (•; y) is continuous on ∂D, g(•; y) is Hölder-continuous on D, (EP4) (uniform ellipticity assumption) there exists c > 0 such that d i,j=1 σ(x; y)σ(x; y) T ij

Theorem 5 . 2 (

 52 Feynman-Kac formula for a backward parabolic parameter-dependent PDE on bounded domain). Let T > 0 and assume that, for all y in Y (PP1) assumptions (A1) and (A2) are satisfied for b(•; y) and σ(•; y), (PP2) D ⊂ R d is a bounded domain, (PP3) the boundary ∂D is regular for X t,x;y in the sense that ∀t ∈ [0, T ], ∀x ∈ ∂D, P(τ t,x;y = 0) = 1, (5.8) with τ t,x;y = inf{s > t : X t,x;y s / ∈ D} the first exit time of X t,x;y from D, (PP4) the functions f (•; y), g(•; y) : [0, T ] × D → R are continuous, (PP5) there exists a continuous function u(•; y) : [0, T ] × D → R, of class C 1,2 on all open subsets of [0, T ) × D verifying (5.6).

Theorem 5 . 3 (

 53 Feynman-Kac formula for e i (y) = u(y) -u i (y) for elliptic parameter-dependent PDEs). Assume that, for all y in Y , (EP1)-(EP3) hold and u i (y) is of class C(D) ∩ C 2 (D) such that L y u i (y) is Hölder-continuous on D.

Theorem 5 . 6 (

 56 Feynman-Kac formula for e i (y) = u(y) -u i (y) for parabolic parameter-dependent PDEs). Assume that, for all y in Y , (PP1)-(PP5) hold and assume moreover that u i (y) is continuous on [0, T ] × D, of class C 1,2 on all open sets of [0, T ) × D and such that A y u i is continuous on [0, T ] × D.

Theorem 5 . 7 (

 57 Reinterpretation of e i (y) 2 L 2 ). Taking V ∼ U(D) uniformly distributed on D and T ∼ U(0, T ) uniformly distributed on [0, T ], we have for any y in Y

Algorithm 5 . 1 1 :

 511 Probabilistic Greedy Algorithm for RBM Require: (λ i ) 1≤i≤n ∈ (0, 1), ∈ (0, 1), Ξ a discrete training set. Set H 0 = {0}.2: for i = 1 . . . n do 3:

  which means that for all i = 1, . . . , n u(ξ i ) -P H i-1 u(ξ i ) L 2 ≥ √ 1 -C max ξ∈Ξ u(ξ) -P H i-1 u(ξ) L 2 .(5.20)

Algorithm 5 . 2 1 :

 521 Greedy Algorithm with Random Training Sets[START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] Require: (Ξ i ) i=1...m 2 a collection of discrete training sets. Set i = 0, H 0 = {0} and σ0 = 8m .2: while σi ≥ 8m and i < m 2 do 3:

  ) -P (y) H ≤ Cm -r , m ≥ 1.(5.24)

3 4m 2 )

 2 N ≤ η m 2 ; • (Ξ i ) i=1...m 2 ⊂ Y is a collection of discretetraining sets obtained from N samples drawn according to the uniform distribution over Y = [-1, 1] d .

7: end while Theorem 5 . 12 .

 512 Under assumptions of Theorem 5.10 and Theorem 5.11, we have that Algorithm 5.3 is a γ-weak-greedy algorithm over Y with γ = √ 1-8mC , and with probability greater than 1 -η -m 2 i=1 λ i .

L 2 =

 2 u(y) -u i (y) 2 L 2 = E i+1 [Z i+1 (ξ)] . (5.27)Taking f the final step of the algorithm, letA = ∩ f i=1 A i and B ∩ f i=1 B i with

Titre:

  Approximation et réduction de modèle pour les équations aux dérivées partielles avec interprétation probabiliste Mot clés : équations aux dérivées partielles, approximation en grande dimension, méthodes de bases réduites, stratégies adaptatives, méthodes de Monte-Carlo, optimisation discrète Résumé : Nous nous intéressons dans cette thèse à la résolution numérique de modèles régis par des équations aux dérivées partielles admettant une interprétation probabiliste. Dans un premier temps, nous considérons des équations aux dérivées partielles en grande dimension. En nous basant sur une interprétation probabiliste de la solution qui permet d'obtenir des évaluations ponctuelles de celle-ci via des méthodes de Monte-Carlo, nous proposons un algorithme combinant une méthode d'interpolation adaptative et une méthode de réduction de variance pour approcher la solution sur tout son domaine de définition.Dans un deuxième temps, nous nous intéressons aux méthodes de bases réduites pour les équations aux dérivées partielles paramétrées. Nous proposons deux algorithmes gloutons reposant sur une interprétation probabiliste de l'erreur. Nous proposons également un algorithme d'optimisation discrète probably approximately correct en précision relative qui nous permet, pour ces deux algorithmes gloutons, de sélectionner judicieusement un snapshot à ajouter à la base réduite en se basant sur la représentation probabiliste de l'erreur d'approximation.Title: Approximation and model reduction for partial differential equations with probabilistic interpretationKeywords: partial differential equations, high-dimensional approximation, reduced basis methods, adaptive strategies, Monte-Carlo methods, discrete optimization Abstract: In this thesis, we are interested in the numerical solution of models governed by partial differential equations that admit a probabilistic interpretation.

  1.1 Summary of the advantages and drawbacks of each numerical timeintegration technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Algorithm 2.2 computed on the exact solution of (TC1): evolution of #Λ n , error criterion ε n and interpolation errors in norms L 2 and L ∞ at each step n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 (TC1) Comparison of the algorithmic complexity to reach the precision 3 • 10 -5 , with (∆t, M ) = (10 -4 , 1000). . . . . . . . . . . . . . . . 2.3 (TC2) Comparison of Algorithm 2.2 (first four columns) and Algorithm 2.3 (last four columns). . . . . . . . . . . . . . . . . . . . . . . 2.4 (TC3) Comparison of Algorithm 2.2 (first four columns) and Algorithm 2.3 (last four columns). . . . . . . . . . . . . . . . . . . . . . .

  .1.

	Problem studied Technique used Advantages/Drawbacks	Complexity
	(1.14) or (1.13)	Euler-Maruyama	Easy to implement and to apply	O( -1 ) or O( -2 )
		scheme	to any diffusion process but low	depending on the
			convergence order	regularity for the
						functions in (1.14)
		High-order	High convergence order com-	O( -γ ), γ ∈ (0, 1]
		schemes [71]	pared	to	Euler-Maruyama	depending on the
			scheme but not available when	order of the scheme
			the SDE has constant coeffi-
			cients, less easy to implement in
			high-dimension
		Boundary shifting	Easy to implement	O( -1 )
		[50]
	Poisson equation on	Romberg extrapo-	Allow high-order of convergence
	unbounded domain	lation [99, section	even when the SDE has constant
		2.3]	coefficients but difficult to ex-
			tend to complex function of X t,x

  2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.2 Probabilistic Tools for Solving PDEs . . . . . . . . . . . . 47 2.2.1 Pointwise Evaluations of the Solution . . . . . . . . Main Contributions . . . . . . . . . . . . . . . . . . . . . . 54 2.4.1 Perturbed Version of Sparse Interpolation . . . . . 54 2.4.2 Adaptive Version of Sequential Control Variate Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 54 2.4.3 Numerical Results . . . . . . . . . . . . . . . . . . 55 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

	2.2.2 A Sequential Control Variates Algorithm . . . . . . 2.3 Adaptive Sparse Interpolation . . . . . . . . . . . . . . . . 2.3.1 Sparse Interpolation . . . . . . . . . . . . . . . . . 2.3.2 Adaptive Algorithm for Sparse Interpolation . . . . 2.4	47 50 52 52 53

Table 2 .

 2 1: Algorithm 2.2 computed on the exact solution of (TC1): evolution of #Λ n , error criterion ε n and interpolation errors in norms L 2 and L ∞ at each step n.

	Λ 6 Λ 8 Λ 10	10 20 26 30 33 40 41 50 55	-01 1.261601e+00 4.213566e+00 2.792486e-02 1.204421e-01 3.602629e-01 2.178450e-05 9.394419e-04 3.393999e-03 9.632815e-07 4.270457e-06 1.585129e-05 9.699704e-08 2.447475e-06 8.316435e-06 4.114730e-09 2.189518e-08 9.880306e-08 1.936050e-10 6.135776e-10 1.739848e-09 1.008412e-11 9.535433e-11 4.781375e-10 1.900248e-14 1.004230e-13 4.223288e-13 7.453467e-15 2.905404e-14 1.254552e-13
	We present now the behavior of Algorithm 2.3. Simulations are performed with a bulk-chasing parameter θ = 0.5. At each step n of Algorithm 2.3, we use Algorithm 2.1 with (∆t, M ) = (10 -4 , 1000), stopped when a stagnation is detected. As shown on the left plot of Figure 2.3, for #Λ n = 55 we reach approximately a precision of 10

-14

as for Algorithm 2.2 performed on the exact solution (see Table

2

.1). According to the right plot of Figure

2

.3, we also observe that the enrichment procedure behaves similarly for both algorithms (ε n and ε n are almost the same). Here using the approximation provided by Algorithm 2.1 has a low impact on the behavior of

  Moreover, we observed that while adaptive version of Algorithm 2.1 stagnates at a precision of 3 • 10 -5 , Algorithm 2.3, with the same parameters ∆t and M , converges almost up to the machine precision. This is why the high-dimensional test cases will be run only with Algorithm 2.3.

	Complexity	4 • 10 9 operations	16 • 10 9 operations	10 12 N operations
	Table 2.2: (TC1) Comparison of the algorithmic complexity to reach the precision 3 • 10

-5

, with (∆t, M ) = (10 -4 , 1000).

We observe that both the adaptive version of Algorithm 2.1 and Algorithm 2.3 have a similar complexity, which is better than for the full-grid version of Algorithm 2.1.

Table 2 .

 2

	ũΛn 2

4: (TC3) Comparison of Algorithm 2.2 (first four columns) and Algorithm 2.3 (last four columns).

  Adaptive Algorithm . . . . . . . . . . . . . . . . . . 69 3.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . 71 3.A Intermediate Results . . . . . . . . . . . . . . . . . . . . . 75 3.B Proof of Proposition 3.4 . . . . . . . . . . . . . . . . . . . 77

	. 3.2 Monte-Carlo Estimate with Guaranteed Relative Precision 3.2.1 Monte-Carlo Estimate . . . . . . . . . . . . . . . . 3.2.2 Complexity Analysis . . . . . . . . . . . . . . . . . 3.3 Optimization Algorithms with Guaranteed Relative Precision 66 62 63 64 66 3.3.1 Non-Adaptive Algorithm . . . . . . . . . . . . . . . 67 3.3.2

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . .

  Moreover, the value ξ returned by Algorithm 3.1 satisfies (3.21).

	Proof. The assumptions on (d m ) m≥1 in (3.22) combined with E[Z(ξ)] = 0 ensure that for all ξ in Ξ , M (ξ) is almost surely finite. Then, for all ξ in Ξ, the estimate Êm(ξ) [Z(ξ)] is well defined. Applying Proposition 3.2 for each Z(ξ) with δ = λ/#Ξ and = τ 2+τ , we obtain (3.23). Now let

  For all m ≥ min (M σ 2 , M 2 µ 2 ), i.e. m ≥ M σ 2 or m ≥ M 2 µ 2 , we obtain from(3.45)

	2 µ 2	2	p log	2p(b -a) 2 2 µ 2	+ log	3 cδ	.
	and (3.46), or (3.45) and (3.47), that
	(3.48)

  the decay of d n (M) H is too slow, e.g. polynomial in n with a degree lower or equal to 1, building a good reduced space H n with a greedy algorithm is not necessarily feasible. However, due to Theorem 4.2, the loss in the decay of the Kolmogorov n-width resulting from the use of snapshots to build the reduced spaces is much smaller than the loss resulting from the use of greedy algorithm. Moreover since we have d n (M) H ≤ σ n (M) H , any upper bound for σ n (M) H (see e.g. Theorem 4.1) holds also for d n (M) H .

  which is an even more important limitation when u is not accessible) and 4) the difficulty to solve the optimization problem over a continuous set Y . The last point is due to the complete absence of knowledge on properties (convexity, continuity ...) of y → u(y) -P H i-1 u(y) H . To circumvent these problems, we recall a realistic framework for RBM in Section 4.2.1 before showing how it is possible to derive a weaker version of Algorithm 4.1 in Section 4.2.4.

  is the source term and f : ∂D × Y → R is the boundary condition. The differential operator L y is also the infinitesimal generator of the stochastic process X x;y which is solution of the following Stochastic Differential Equation (SDE)

	X x;y t	= X x;y 0 +	0	t	b(s, X x;y s ; y)ds
		+	0	t	σ(s, X x;y

s ; y)dW s , X x;y 0

  ). Assuming the existence of probabilistic representations for the gradient of u(y) and u i (y), it would be possible to consider probabilistic interpretation of other squared 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 nonlinear[END_REF] Corollary IV.5.2].

	Remark 5.5.

∆t (t, x) = E[G(u, X t,x )] -E[G(u, X t,x,∆t )]

List of Algorithms

Proof. We have

(3.12) Then using (3.7) and (3.10), we deduce that It suffices to consider the case E[Z] > 0 and we have

Also

which concludes the proof.

In practice, the computation of the estimate given by (3.9) requires a particular choice for the random variable M and for the sequence (d m ) m≥1 . A natural choice for M which satisfies (3.8) is

we have that c m converges to 0 almost surely. Moreover if E[Z] = 0, it is sufficient to ensure that M < +∞ almost surely.

Remark 3.3. When choosing M as in (3.14), the estimate defined by (3.9) is the one proposed in [78, equation (3.7)]. A variant of this estimate can be found in [START_REF] Mnih | Empirical bernstein stopping[END_REF].

with

If M f ≤ M , (3.50) and (3.52) imply M = M , where • denotes the ceil function.

Otherwise M f > M and we obtain M ≤ M f . Thus, it provides the following upper bound

Since M ≥ M and A ∩ B implies (3.53), we deduce that A ∩ B implies (3.44), which concludes the proof of the first result.

Let us now prove the result in expectation. Let K := 2

which ends the proof.

Chapter 4

Reduced Basis Methods

This chapter surveys reduced basis methods (RBM) for parameter-dependent problems with a particular focus on parameter-dependent partial differential equations (PDEs). We assume that the problem (4.1) is well-posed in H, which means that for all instance y ∈ Y , there exists a unique solution u(y) ∈ H that satisfies (4.1). The map u is defined on a high-dimensional set O × Y ⊂ R d × R p , which causes some Theorem 4.4. Assume that the map u : Y → H is Hölder-continuous with reference to y, which means

Contents

a constant L independent of y 1 and y 2 . Moreover, assume that Ξ is an -net of Y , which means Ξ is a finite subset of Y such that for all y in Y there exists a ξ in Ξ satisfying

Then we have

Proof. We have

Then we notice that for all y in Y and all space W ⊂ H

for any ξ in Ξ. In particular, for ξ y such that ξ y -y R D ≤ (that exists by assumption), it brings

Taking the supremum over Y and the infimum over W , we obtain the second inequality.

Theorem 4.4 ensures that if an -net Ξ is dense enough in Y ⊂ R p we can expect similar behavior between d n (M) H and d n (M Ξ ) H , up to a threshold depending on the parameter . However when the parametric dimension p is large, i.e. p 1, the number of points of -net of Y is too large since it grows exponentially with the dimension p: this is another phenomenon illustrating the curse of dimensionality.

To avoid using an -net, the authors of [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] have proposed to use random training sets, introducing at the same time the notion of probabilistic greedy algorithm. This contribution will be discussed further in Chapter 5. Another contribution [START_REF] Hesthaven | Efficient greedy algorithms for highdimensional parameter spaces with applications to empirical interpolation and reduced basis methods[END_REF] has first investigated ways to avoid computing at each step of a greedy algorithm the error u(ξ) -P H i-1 u(ξ) H over the whole training set Ξ: this technique reduces the computational cost of using large training sets but not necessarily the storage problem. The authors proposed in the same paper a heuristic method to avoid considering too big training sets by safety checking the approximation error of the reduced space they had built on random points outside of the training set. They came up with many numerical experiments illustrating the interest of their method. They also allowed their algorithm to remove and add points to the training set, which thus changes with the number of steps.

Chapter 5 Probabilistic Reduced Basis Methods

In this chapter we present our third contribution. This contribution consists in two probabilistic greedy algorithms for reduced basis methods (RBM). In order to present these algorithms, we first provide a Feynman-Kac representation of the approximation error which allows us interpret its squared norm as an expectation and thus to use Monte-Carlo methods to estimate this squared norm. Then combining this Monte-Carlo estimation with the probably approximately correct (PAC) algorithm from Chapter 3 we come up with a new technique to adaptively select snapshots for the construction of the reduced basis. where the unknown is a real-valued function u(y) := u(•; y) belongs to some Hilbert space H.

In Chapter 4 traditional deterministic reduced basis methods (RBM) were presented for (5.1). Here we present probabilistic reduced basis methods where probabilistic methods are considered to select the snapshot that enriches the reduced space.

Probably Approximately Correct Algorithm

As explained in Chapter 4, solving the global optimization problem (5.10) is not realistic when Y is not a discrete set. We rather consider the following discretized optimization problem max ξ∈Ξ u(ξ) -u i (ξ) H , (5.17) where Ξ ⊂ Y is a finite set of parameters. Thus, using Theorem 5.4 or Theorem 5.7 when • H = • L 2 , (5.17) becomes a discrete optimization problem where the quantity to optimize is the expectation of some parameter-dependent random variable.

We describe the procedure to sample X W ;ξ (with W = V for the elliptic case or W = (T, V ) for the parabolic case):

1. draw a sample w according to the distribution of W ; 2. simulate using a time-integration scheme an approximate trajectory of the diffusion process X w;ξ that will be denoted X ∆t,w;ξ , see (1.15).

These samples are used to provide Monte-Carlo estimates of u(ξ) -u i (ξ) 2 L 2 for different values of ξ in Ξ. However our goal is not to provide precise estimates for u(ξ) -u i (ξ) 2 L 2 for every value ξ in Ξ, which would require a high number of samples from X ∆t,W ;ξ . We rather use the PAC algorithm from Chapter 3 to solve (5.17) in high probability with a guaranteed relative quasi-optimality, sampling less the random variables Z i+1 (ξ) = |O|F i (W, X W,ξ ; ξ)F i (W, XW,ξ ; ξ) (|O| is Lebesgue measure of the physical domain O and is equal to|D| in the elliptic case and to T |D| in the parabolic case). This algorithm returns ξ i+1 in Ξ such that (5.18) for any prescribed threshold in (0, 1) and probability of failure λ in (0, 1), where

Remark 5.8. For parabolic parameter-dependent PDEs, the random variables Z i+1 (ξ) are bounded for all ξ in Ξ since τ t,x;ξ 1 ≤ T , see (5.9). This is not the case for elliptic parameter-dependent PDEs since we usually cannot provide a deterministic bound for the exit time τ x;ξ (see (5.5)). Remark 5.9. We chose here not to consider

instead of Z i+1 (ξ) in order to avoid dealing with the time-integration error.

Conclusion and Outlook

algorithm.

These discrete optimization techniques could have been considered as remote from the original interest of this thesis, which is partial differential equations. However, such optimization techniques could be adapted to find extrema of the solution of a PDE. Indeed, when a solution admits a probabilistic representation, it can be written as an expectation of some random process, as we have seen in Chapter 1.

Moreover this probably approximately correct algorithm in relative precision can be used in reduced basis methods. We have surveyed first in Chapter 4 standard reduced basis methods with a particular focus on parameter-dependent PDEs. Then we have proposed two probabilistic greedy algorithms based on a probabilistic interpretation of the squared norm of the a posteriori approximation error. This squared norm written as an expectation allows to select the next snapshot to enrich the current reduced space using the adaptive probably approximately correct algorithm in relative precision presented in Chapter 3. Our first algorithm used PAC algorithm for snapshots' selection on a fixed training set. Our second algorithm used PAC algorithm for snapshots' selection on random training sets [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF]. Both algorithms offer a probabilistic selection of the snapshots and the second one is promising for high-dimensional parameter-dependent PDEs.

Further outlooks include goal-oriented reduced basis methods, such as it can be encountered in Uncertainty Quantification (UQ) when we look into some quantity of interest depending on the solution and not the solution itself. Using probabilistic interpretation of the solution to derive probabilistic goal-oriented reduced basis methods could be of high interest when the quantity of interest can be estimated using Monte-Carlo methods. It includes in particular quantities of interest such as integrals of the solution in a deterministic setting or statistical moments of some functional of the solution in a stochastic setting.