Marie Billaud

Marie Billaud-Friess

Fabio Nobile

Keywords: greedy algorithm

come

Introduction (version française)

Mes travaux de recherche relèvent principalement de l'analyse numérique selon deux thématiques distinctes. Les travaux menés durant mes doctorat et post-doctorats ont porté sur le développe-ment de méthodes numériques pour la simulation d'écoulements à plusieurs composants. A mon arrivée à Centrale Nantes, j'ai été amenée à changer de thématique pour me focaliser sur le développement et l'analyse de méthodes d'approximation et réduction de modèle pour les équations aux dérivées partielles. Dans ce manuscrit, j'ai choisi de me concentrer sur mes plus récentes contributions portant sur le second thème. Un bref aperçu de mes travaux de thèse et post-doctorats [MBF1, MBF2, MBF4, MBF3, MBF5, MBF6, MBF7] est néanmoins donné en annexe (voir page 63).

Le présent mémoire se constitue de trois chapitres. Cependant, on peut voir se dessiner deux axes suivant lesquels les contributions présentées peuvent être rangées avec d'une part des méthodes déterministes et d'autre part des approches probabilistes. Par souci de concision, les concepts utilisés seront motivés et introduits au fur et à mesure des enjeux et des objectifs soulignés.

Approximation de faible rang pour les problèmes dépendant de paramètres

La première direction concerne le développement de méthodes d'approximation de faible rang pour la réduction de modèle de problèmes dépendant de paramètres.

Dans le Chapitre 1, des contributions dédiées à la résolution de tels problèmes motivés par des applications en quantification d'incertitude sont présentées. Le chapitre est divisé en deux sections selon les deux points de vue suivants: d'une part les méthodes d'approximation de faible rang sous format de tenseur et d'autre part les méthodes basées sur des projections. En particulier, une méthode d'approximation sous format de tenseur, basée sur des formulations idéales en minimum de résidu pour des problèmes en grande dimension, est d'abord proposée [START_REF] Billaud-Friess | Méthode de réduction de modèle a priori basée sur des formulations idéales en minimum de résidu[END_REF][START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF]. Puis, des méthodes basées sur des projections pour pour l'approximation de quantités d'intérêt à valeur fonctionnelle ou vectorielle [START_REF] Zahm | Projection-Based Model Order Reduction Methods for the Estimation of Vector-Valued Variables of Interest[END_REF] sont discutées. Ces travaux ont été réalisés en collaboration avec A. Nouy et O. Zahm.

Dans le Chapitre 2, les méthodes d'approximation sous format de faible rang, discutées au chapitre précédent, sont étendues pour des problèmes paramétrés dépendant du temps. En particulier, des approches de faible rang dites dynamiques sont discutées. Premièrement, une méthode de type bases réduites dynamiques est présentée [START_REF] Billaud-Friess | Dynamical Model Reduction Method for Solving Parameter-Dependent Dynamical Systems[END_REF]. La méthode proposée surpassant les approches classiques pour des problèmes de type transport, des pistes d'extension pour les systèmes de loi de conservation hyperboliques ont été explorées. Notamment, une méthode de reconstruction dans un cadre volume fini a été proposée et validée pour des problèmes de transport linéaires et dépendant de paramètres [START_REF] Billaud-Friess | Reconstruction of finite volume solution for parameter-dependent linear hyperbolic conservation laws[END_REF]. Dans un second temps, une méthode de faible rang dynamique travaillant dans la variété des matrices de rang fixé est considérée. Cette approche repose sur la paramétrisation géométrique de cette variété, proposée dans [START_REF] Billaud-Friess | Principal bundle structure of matrix manifolds[END_REF], à partir de laquelle un nouveau schéma d'intégration de type splitting a été dérivé pour l'approximation de faible rang dynamique de solutions de systèmes dynamiques matriciels [START_REF] Billaud-Friess | A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds[END_REF]. Ces diverses contributions ont été le fruit d'une collaboration avec T. Heuzé, A. Falcó et A. Nouy.

Méthodes d'approximation probabilistes pour les EDP

Le deuxième axe de ce manuscrit concerne le développement de méthodes numériques pour l'approximation de la solution d'équations aux dérivées partielles (éventuellement dépendant de paramètres) en utilisant des estimations ponctuelles de la fonction recherchée.

Dans ce but, des approches probabilistes sont discutées dans le Chapitre 3. La clef de voûte des contributions présentées est le théorème de représentation de Feynman-Kac qui permet d'expliciter la solution d'une équation aux dérivées partielles, évaluée ponctuellement, comme l'espérance d'une fonctionnelle stochastique. En utilisant cet outil, une méthode probabiliste d'interpolation parcimonieuse a été proposée pour des équations en grande dimension [START_REF] Billaud-Friess | A probabilistic reduced basis method for parameter-dependent problems[END_REF]. Par ailleurs, une méthode de bases réduites probabiliste basée sur des estimations ponctuelles est proposée [START_REF] Billaud-Friess | A probabilistic reduced basis method for parameter-dependent problems[END_REF] pour la résolution d'EDP paramétrées. En particulier, une telle méthode s'appuie sur une procédure glouton avec un estimateur d'erreur probabiliste combiné avec un algorithme de type bandit adaptatif "probably approximately correct" [START_REF] Billaud-Friess | A PAC algorithm in relative precision for bandit problem with costly sampling[END_REF]. Ces travaux ont été réalisés en collaboration avec A. Macherey, A. Nouy et C. Prieur.

Introduction

My research fits in the domain of numerical analysis. The work done during my PhD and post-doctorates concerned the development of numerical methods for the simulation of multicomponent flows. Since I arrived in Centrale Nantes, my research work has mainly focused on the development of approximation and model order reduction methods for partial differential equations. In this manuscript, I have made the choice to not present in detail the work [MBF1, MBF2, MBF4, MBF3, MBF5, MBF6, MBF7] done before arriving in Centrale Nantes. Only a concise summary is provided in appendix (see page [START_REF] Dick | High-dimensional integration: The quasi-Monte Carlo way Acta Numer[END_REF]. This document is mainly composed of three chapters, following two principal directions with deterministic and probabilistic approximation methods.

Low-rank approximation for parameter-dependent problems

The first direction concerns the development of suitable low-rank approximation methods for model order reduction of parameter-dependent problems.

In Chapter 1, contributions devoted to the resolution of such problems motivated by uncertainty quantification are presented. The chapter is divided in two parts following a natural presentation of the proposed approaches within the framework of low-rank approximation mainly : low-rank approximation in tensor format and projection based low-rank approximation for parameter-dependent problems. More particularly tensor approximation methods based on ideal minimal residual formulations for the solution of high-dimensional problems [START_REF] Billaud-Friess | Méthode de réduction de modèle a priori basée sur des formulations idéales en minimum de résidu[END_REF][START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF] are first proposed. Then, goal oriented projection based model order reduction methods for the estimation of vector-valued variables of interest [START_REF] Zahm | Projection-Based Model Order Reduction Methods for the Estimation of Vector-Valued Variables of Interest[END_REF] are discussed. These contributions were done in collaboration together with A. Nouy and O. Zahm. In Chapter 2, low-rank approximation methods discussed in Chapter 1 are extended for the approximation of time and parameter-dependent problems. Here, dynamical low-rank approximation methods are discussed following two subsequent directions. First, dynamical reduced basis method for parameter-dependent dynamical systems that can be interpreted as a dynamical low-rank approximation approach with a subspace point of view is presented [START_REF] Billaud-Friess | Dynamical Model Reduction Method for Solving Parameter-Dependent Dynamical Systems[END_REF]. As the proposed method is shown to perform better than usual reduced basis for transport (dominated) problems, possible extension to parameter-dependent hyperbolic conservation laws is discussed. In this direction a reconstruction method in finite volume setting has been developed and validated for linear parameter-dependent transport problems [START_REF] Billaud-Friess | Reconstruction of finite volume solution for parameter-dependent linear hyperbolic conservation laws[END_REF]. Secondly, we consider a dynamical low-rank approximation method that works directly in the set of fixed rank matrices. Based on a suitable geometric parametrization of this set [START_REF] Billaud-Friess | Principal bundle structure of matrix manifolds[END_REF], a new splitting integration scheme for the approximation of the solution of matrix dynamical systems [START_REF] Billaud-Friess | A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds[END_REF], arising from discretization of parameter-dependent dynamical systems, has been proposed. These works were the result of a joint work with T. Heuzé, A. Falcó and A. Nouy.

Probabilistic approximation methods for PDEs

The second direction deals with the problem of constructing an approximation of the solution of a partial differential equation (possibly parameter-dependent) using approaches that rely on sampled estimates of the function to approximate.

To that goal, probabilistic approaches are presented in Chapter 3. The main corner stone of the proposed contributions is the Feynman-Kac representation formula of the solution of a partial differential equation. Using this key ingredient, a probabilistic sparse polynomial interpolation method has been proposed to deal with high dimensional problem [START_REF] Billaud-Friess | Stochastic Methods for Solving High-Dimensional Partial Differential Equations[END_REF]. Then reduced basis method using pointwise estimates is discussed [START_REF] Billaud-Friess | A probabilistic reduced basis method for parameter-dependent problems[END_REF] for solving parameterdependent partial differential equations. Especially, it relies on a greedy procedure with a probabilistic error estimate gathered with a probably approximately correct bandit algorithm proposed in [START_REF] Billaud-Friess | A PAC algorithm in relative precision for bandit problem with costly sampling[END_REF]. These contributions were done in collaboration with A. Macherey, A. Nouy and C. Prieur.

Chapter 1

Low-rank approximation for parameter-dependent equations Parameter-dependent problems are ubiquitous in large variety of applications such as uncertainty quantification, sensitivity, inverse problems, control or optimization. This chapter presents contributions devoted to the resolution of such problems with application in uncertainty quantification using low-rank approximation methods.

The chapter is divided in two parts following a natural presentation of the proposed approaches within the framework of low-rank approximation : low-rank approximation in tensor format and projection based low-rank approximation for parameter-dependent problems. After a brief introduction on low-rank tensor approximation methods for model order reduction of parameterdependent problems in Section 1.1, we present in Section 1.2 a tensor approximation method based on ideal minimal residual formulation for the solution of high-dimensional problems. Then, in Section 1.3, we propose projection based methods for the estimation of vector-valued variables of interest.

General context

In this chapter, we are concerned with the resolution of high-dimensional partial differential equations (PDEs) arising in physics or stochastic calculus, or the solution of parameterdependent or stochastic equations. The following general linear parameter-dependent equation is considered A(ξ)u(ξ) = b(ξ).

(1.1.1)

Here, A is a linear operator, depending on the (random) parameters ξ = (ξ 1 , . . . , ξ p ) T that take their values in some set Ξ ⊂ R p (with p possibly large). The parameter-dependent solution is u(ξ) which belongs to some vector space V .

Low-rank approximation for model reduction

Computing a solution u to the problem (1.1.1) rises various difficulties, related to the complexity and cost of the computations, when one wants to access to pointwise evaluations of u over Ξ or is interested in having a global approximation of the solution map u : Ξ → V , in particular in high-dimension (when p is large). As discussed in what follows, low-rank approximation methods are nowadays widely used for the solution of problem (1.1.1) in these contexts. We refer e.g. the reader to the surveys [START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF] or [START_REF] Nouy | Low-rank tensor methods for model order reduction[END_REF][START_REF] Nouy | Chapter IV : Low-rank methods for high-dimensional approximation and model order reduction[END_REF].

In a multiquery context, one has to evaluate the solution u of (1.1.1) for many instances of the parameter ξ. For complex numerical model, arising for example from the standard linear approximation of a PDEs with dim(V ) = n (n large), the computation of the solution for one instance of the parameter can be very costly. To tackle this problem, Model Order Reduction (MOR) methods have emerged these last decades. The idea of MOR is to provide an accurate surrogate of the solution map u : Ξ → V cheaper to compute, and allowing rapid evaluation of the solution or quantities of interest. Such methods usually rely on linear approximation since there are aimed at approximating u as follows

u(ξ) ≈ r i=1 α i (ξ)v i , (1.1.2)
where {v 1 , . . . , v r } are elements in V and {α 1 , . . . , α r } are elements of a some vector space S of functions defined on Ξ. Such an approximation can be interpreted as a rank-r approximation in the tensor space X := V ⊗ S. Standard MOR methods aim at finding an approximation under the form (1.1.2) with a small number of terms r << n. To compute an approximation under the form (1.1.2) we can distinguish two kinds of approaches depending on either they rely on constructing a basis {v i } r i=1 ⊂ V , or directly computing an approximation under the form (1.1.2) in a low-rank subset of X. The approaches of the first kind are called projection based low-rank approximation methods, since the approximation is obtained by a projection (e.g. Galerkin, Least-Square, interpolation . . . ) of u in the subspace of V spanned by a basis selected a priori. Classical projection based low-rank methods include the Proper Orthogonal Decomposition (POD) method, the Empirical Interpolation Method (EIM) or the Reduced Basis (RB) method. The approaches of the second kind include low-rank tensor methods that rely on direct optimization in low-rank manifold. Classical examples are the best r-term approximation that coincides with the rank-r truncated Singular Value Decomposition (SVD), as well as the Proper Generalized Decomposition (PGD) methods. All, these methods differ by their selection of the reduced basis for different controls of the error over the parameter set (uniform control for EIM and RB or control in mean-square sense for POD and PGD see, e.g., [START_REF] Nouy | Low-rank tensor methods for model order reduction[END_REF]).

Computing a global approximation of the solution u as a function of multiple parameters ξ = (ξ 1 , . . . , ξ p ) ∈ Ξ, is challenging in particular for high-dimensional problems. A widespread idea is to construct some suitable approximation of the solution map using a separated representation under the form

u(ξ) ≈ r 1 ν 1 =1 • • • rp νp=1 a ν 1 ,...,νp φ 1 ν 1 (ξ 1 ) . . . φ p νp (ξ p ). (1.1.3)
Seeking an approximation in a tensorized basis is not fully satisfactory because of the so-called curse of dimensionality [START_REF] Bellman | Dynamic Programming[END_REF]. In brief, the storage or computational complexity for solving (1.1.1) grows exponentially in the parametric dimension p, i.e., it has the form ρ p with ρ > 1.

The last decades significant research effort has been done in order to develop tensor formats whose complexity increases slowly with the number of parameters p. In this manuscript, we focus on low-rank approximations in tensor format that aims at reducing the complexity by exploiting (high-order) low-rank structures of multivariate functions, considered as elements of tensor product spaces. Classical low-rank tensor formats include canonical tensor formats, Tucker formats, Tensor Train formats [START_REF] Holtz | On manifolds of tensors with fixed TT rank[END_REF][START_REF] Oseledets | Tensor-train decomposition[END_REF], Hierarchical Tucker formats [START_REF] Hackbusch | A New Scheme for the Tensor Representation[END_REF] or more general tree-based Hierarchical Tucker formats [START_REF] Falcó | Tree-based tensor formats[END_REF]. Other approaches are also available for the construction of structured approximation of the solution map. Among them, let mention sparse tensor methods that consist in searching an approximation of the form (1.1.3) with only a few non-zero terms. To this end, the idea is to properly select a small number of functions in a set of functions as e.g. multivariate polynomial [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF][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].

Contributions

The contributions presented in this chapter are mainly related to the PhD of O. Zahm (IN-RIA, Grenoble-Rhône Alpes) [START_REF] Zahm | Reduced Order Modeling for parameter dependent equations -Application to uncertainty quantification[END_REF] supervised together with A. Nouy (LMJL, ECN). In this work, we have dealt with the problem of low-rank approximation in the context of uncertainty quantification addressing the following points.

• Generally, the solution u is unknown and only given implicitly through the equation (1.1.1). As a consequence practical low-rank approximation methods often use the residual of (1.1.1) measured with a certain norm. Such methods are not able to provide optimal approximations, but only quasi-optimal approximations in the best case. This loss of accuracy can be problematic for the efficiency of the methods.

• In many applications (e.g. motivated by uncertainty quantification), we only need to compute a partial information, called a quantity of interest (QoI), which is a function of the solution map. In particular, the main challenge is to develop goal-oriented low-rank approximation methods which provide accurate approximation of QoI while reducing the computational complexity.

Low-rank approximation in tensor format

We recall briefly in Section 1.

v µ = v 1 ⊗ . . . ⊗ v p , v µ ∈ X µ
, an elementary tensor. A Hilbert tensor space X (see [START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF] for more details) equipped with the norm • X is obtained by the completion with respect to • X of the algebraic tensor space

X = • X p µ=1 X µ = span{⊗ p µ=1 v µ : v µ ∈ X µ , 1 ≤ µ ≤ p} • X . A natural inner product on X is induced by inner products •, • µ in X µ , 1 ≤ µ ≤ d.
It is defined by v, w X = p µ=1 v µ , w µ µ for v, w two elementary tensors. Finally, we denote by X the topological dual of X and by •, • X ,X the duality pairing between X and X . We denote by R X : X → X the Riesz isomorphism.

Classical tensor formats

Here, we give some brief summary of the most classical tensor formats used in the litterature. The reader is referred to the textbook [START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF] and survey [START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF] for a more exhaustive review of the different tensor formats and their mathematical properties. In this section, for the sake of presentation, we assume that dim(X µ ) = n, µ = 1, . . . , p.

• The set of tensors in canonical format with rank lower than r ∈ N is defined as the set

T can r (X) = v = r i=1 ⊗ p µ=1 v µ i : v µ i ∈ X µ .
Here, r is the canonical rank of the tensor v in X, i.e. the smallest integer such that v can be written under the form r i=1 ⊗ p µ=1 v µ i . The canonical tensor format has a storage complexity in O(rnp).

• The set of tensors in Tucker format with multilinear rank bounded by r = (r 1 , . . . , r p ) ∈ N p is given by

T tuc r (X) =    v = r 1 i 1 =1 . . . rp ip=1 α i 1 ,...,ip ⊗ p µ=1 v µ iµ : v µ iµ ∈ X µ , α i 1 ,...,ip ∈ R    In this format α ∈ R r 1 ו•
•×rp is a tensor of order p. An approximation v in Tucker format, can be interpreted as an approximation in the tensor space

U 1 ⊗ • • • ⊗ U p where U µ = span(v µ 1 , . . . , v µ rp ) is a r p -dimensional subspace of X µ .
Actually, if the multilinear rank is r = (r, . . . , r) (with r = max µ (r µ )) the storage complexity is in O(r p + nrp), which limits the applicability of the Tucker format for very large values of p.

• To overcome such a limitation due to storage complexity, Tree-based (or Hierarchical)

Tucker formats [START_REF] Falcó | Geometric Structures in Tensor Representations (Final Release)[END_REF][START_REF] Falcó | Tree-based tensor formats[END_REF][START_REF] Hackbusch | A New Scheme for the Tensor Representation[END_REF] based on a notion of rank associated with a dimension partition tree have been introduced. As a particular example, let mention the set of tensor in Tensor Train (TT) format [START_REF] Oseledets | Tensor-train decomposition[END_REF] with TT-rank bounded by r = (r 1 , . . . , r p-1 )

∈ N p-1 T T T r (X) =    v = r 1 i 1 =1 . . . rp ip=1 v 1 1,i 1 ⊗ v 2 i 1 ,i 2 ⊗ • • • ⊗ v µ iµ,1 : v µ iµ-1,iµ ∈ X µ    Let r = (r, . . . , r) (with r = max µ (r µ
)) be the the TT-rank. In comparison to the Tucker format, TT format has a storage of lower complexity in O(pr 2 n) (which does not grow exponentially with p).

Best low-rank approximation

Low-rank tensor approximation methods consist in searching for approximations in a subset of tensors of X noted

M ≤r (X) = {v ∈ X : rank(v) ≤ r} ,
where the notion of rank(v) depends of the tensor format used for the approximation v, as previously discussed (i.e. canonical rank, multilinear rank or TT-rank). Then, the best approximation of u ∈ X in a given set of tensor format M ≤r with respect to the norm • X is

min v∈M ≤r (X) u -v X . (1.2.1)

Existence of a solution

Low-rank tensor subsets M ≤r are neither linear subspaces nor convex sets. However, they usually satisfy topological and geometrical properties that make the best approximation problem (1.2.1) meaningful M ≤r (X) is weakly closed (or simply closed in finite dimension), (1.2.2)

M ≤r (X) ⊂ γM ≤r (X) for all γ ∈ R. (1.2.3) Property (1.2.
3) is satisfied by all the classical tensor subsets (i.e. T can r (X), T tuc r (X) or T T T r (X)). Property (1.2.2) ensures the existence of minimum in the set M ≤r (X) to problem (1.2.1). This property, under some suitable conditions on the norm • X (which is naturally satisfied in finite dimension), is verified by most tensor subsets used for approximation. For the the set M ≤r (X) with the notion of canonical rank for p ≥ 3 the property (1.2.2) is only fullfiled for the set of elementary tensors M r≤1 (X) [START_REF] Falcó | A Proper Generalized Decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF]. Subsets of tensors associated with the notions of multilinear rank or TT-rank have better properties. They are closed sets and they are differentiable manifolds [START_REF] Falcó | Geometric Structures in Tensor Representations (Final Release)[END_REF][START_REF] Holtz | On manifolds of tensors with fixed TT rank[END_REF][START_REF] Uschmajew | Geometric Methods on Low-Rank Matrix and Tensor Manifolds[END_REF]. This has useful consequences in practice for optimization [START_REF] Uschmajew | Geometric Methods on Low-Rank Matrix and Tensor Manifolds[END_REF] or for the projection of dynamical systems in these manifolds (see, e.g., [START_REF] Lubich | Dynamical approximation by hierarchical tucker and tensor-train tensors[END_REF] and Chapter 2).

Practical aspects

For the case p = 2, all classical low-rank tensor formats coincide with the canonical format so that M ≤r (X) = T can r (X) for some rank r. When the norm • X is the canonical norm, then u r coincides with a rank-r SVD of u. Efficient algorithms for computing the SVD can therefore be applied. For d > 2 and when the norm • X is the canonical norm, different algorithms based on optimization methods have been proposed for the different tensor formats (see, e.g., [START_REF] Espig | A regularized newton method for the efficient approximation of tensors represented in the canonical tensor format[END_REF][START_REF] Holtz | The Alternating Linear Scheme for Tensor Optimisation in the TT format[END_REF] or [START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF] for a recent review). Very efficient algorithms based on higher order SVD have also been proposed in [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF], [START_REF] Grasedyck | Hierarchical singular value decomposition of tensors[END_REF] and [START_REF] Oseledets | Breaking the curse of dimensionality, or how to use SVD in many dimensions[END_REF], respectively for Tucker, Hierarchical Tucker and Tensor Train tensors. Note that these algorithms provide quasi-best approximations with a quasi-optimality factor bounded by a function of the dimension d. For a general norm • X , the computation of a global optimum to the best approximation problem is still an open problem for all tensor subsets, and methods based on SVD cannot be applied anymore. However, classical optimization methods can still be applied (such as Alternating Least Square (ALS)) in order to provide an approximation of the best approximation [START_REF] Espig | A regularized newton method for the efficient approximation of tensors represented in the canonical tensor format[END_REF][START_REF] Rohwedder | On local convergence of alternating schemes for optimization of convex problems in the tensor train format[END_REF][START_REF] Uschmajew | Local Convergence of the Alternating Least Squares Algorithm for Canonical Tensor Approximation[END_REF].

Approximation of the solution of high-dimensional problems

Given X (resp. Y ) a Hilbert tensor space equipped with inner product •, • X (resp. •, • Y ) and associated norm • X (resp. • Y ), we consider the following general problem Au = b, u ∈ X.

(1.2.4)

with A ∈ L(X, Y ). We assume that A is weakly coercive operator satisfying the assumptions of the Nečas' theorem (see, e.g." [71, Chapter 2] for details). In particular, we have for all v ∈ X,

α v X ≤ Av Y ≤ β v X , (1.2.5)
The constants α, β stand for the weakly coercive and continuity constants of the operator A defined by

inf 0 =v∈X sup 0 =w∈Y Av, w Y ,Y v X w Y := α > 0 and sup 0 =v∈X sup 0 =w∈Y Av, w Y ,Y v X w Y := β < ∞,
Under these assumptions the problem (1.2.4) is well-posed, and the condition number of the operator A is

κ(A) = β α ≥ 1.

Usual algorithms

Here, we want to compute an approximation in M ≤r (X) of the solution u ∈ X of the problem (1.2.4). As u is not explictly given, the best approximation problem (1.2.1) cannot be solved directly. Low-rank approximation methods in tensor format typically rely on the definition of approximations based on the residual Aub, which is a computable quantity. Different strategies have been proposed for the construction of low-rank approximations of the solution of equations in tensor format.

The first family of methods rely on classical iterative algorithms for linear or nonlinear systems of equations with low-rank tensor algebra (using low-rank tensor compression) for standard algebraic operations [START_REF] Ballani | A projection method to solve linear systems in tensor format[END_REF][START_REF] Khoromskij | Tensor-structured Galerkin approximation of parametric and stochastic elliptic PDEs[END_REF][START_REF] Kressner | Low-rank tensor krylov subspace methods for parametrized linear systems[END_REF][START_REF] Matthies | Solving stochastic systems with low-rank tensor compression[END_REF]. The second family of methods consists in directly computing an approximation of u in M ≤r (X) by minimizing some residual norm [START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF][START_REF] Doostan | A least-squares approximation of partial differential equations with high-dimensional random inputs[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF] min

v∈M ≤r (X) Av -b Y , (1.2.6)
for some choice of the norm • Y . To that goal, constructive approaches well known as greedy algorithms have been proposed (see [START_REF] Temlyakov | Greedy approximation[END_REF] and monograph [START_REF] Temlyakov | Greedy Approximation[END_REF]). Such algorithms, often referred as the PGD [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF][START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF] consist in computing successive corrections in a small low-rank tensor subset, typically the set of rank-1 canonical tensors. A step r of this algorithm can be summarized as follows:

w r ∈ arg min w∈T can 1 (X) A(u r-1 + w) -b Y , (1.2.7
)

u r = u r-1 + w r . (1.2.8)
where u r stands for a rank-r approximation in T can r (X). In practice, the rank-one approximation that defines the correction w r can be solved using ALS algorithm. These greedy algorithms have been analyzed in several papers [START_REF] Ammar | On the convergence of a greedy rank-one update algorithm for a class of linear systems[END_REF][START_REF] Cances | Convergence of a greedy algorithm for highdimensional convex nonlinear problems[END_REF][START_REF] Cances | Greedy algorithms for high-dimensional nonsymmetric linear problems[END_REF][START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor banach spaces[END_REF][START_REF] Falcó | A Proper Generalized Decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF][START_REF] Figueroa | Greedy approximation of high-dimensional Ornstein-Uhlenbeck operators[END_REF] and a series of improved algorithms have been introduced in order to increase the quality of suboptimal greedy constructions [START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor banach spaces[END_REF][START_REF] Giraldi | Tensor-based methods for numerical homogenization from high-resolution images[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF][START_REF] Nouy | Proper Generalized Decompositions and Separated Representations for the Numerical Solution of High Dimensional Stochastic Problems[END_REF][START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF].

Good choice of the residual norm

Although minimal residual based approaches are well founded, they generally provide lowrank approximations that can be very far from optimal approximations (1.2.1), with respect to the natural norm • X , when using usual measures of the residual. Indeed, for such a choice, the operator A is often badly conditioned. If we are interested in obtaining an optimal approximation with respect to the norm • X , it implies to choose the residual norm [START_REF] Cohen | Adaptivity and variational stabilization for convection-diffusion equations[END_REF][START_REF] Dahmen | Adaptive petrov-galerkin methods for first order transport equations[END_REF] such that

Av Y = v X , v ∈ X, (1.2.9)
in (1.2.6) to recover the best approximation problem (1.2.1). This implies

I X = R -1 X A * R -1 Y A ⇔ R Y = AR -1 X A * ⇔ R X = A * R -1 Y A, (1.2.10) thus w Y = A * w X , (1.2

.11)

where A * : Y → X stands for the adjoint operator of A. This choice implies that α = β = 1, and therefore κ(A) = 1, meaning that problem (1.2.4) is ideally conditioned. Imposing (1.2.9) can be interpreted as applying implictely some preconditioning operator R -1

X A * R -1 Y to A.
Let us mention, that it is also possible to construct a priori an explicit preconditioner. See, e.g., [START_REF] Desceliers | Polynomial chaos representation of a stochastic preconditioner[END_REF] for stochastic problems and more recently [START_REF] Zahm | Interpolation of Inverse Operators for Preconditioning Parameter-Dependent Equations[END_REF] for parameter dependent algebraic problems. In the latter, the construction of an explicit parameter dependent preconditioner is obtained by interpolation of the inverse operator.

Ideal minimal residual method

For solving (1.2.6), we have proposed in [START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF] the following gradient-type algorithm. Letting u 0 = 0, we construct a sequence {u k } k≥0 in M ≤r (X) and a sequence {y k } k≥0 in Y defined for k ≥ 0 by

y k = R -1 Y (Au k -b) u k+1 ∈ Π r (u k -ρR -1 X A * y k ) (1.2.12)
with ρ > 0. Here, Π r : X → M ≤r (X) denotes the set-valued map that associates to an element u ∈ X the set of best approximations of u in M ≤r :

Π r (u) = arg min v∈M ≤r (X) u -v X . (1.2.13)
The algorithm defined by (1.2.12) can be interpreted as a preconditioned gradient algorithm with an implicit preconditioner. Let us observe at step k we have

u k+1 ∈ Π r (u k -ρR -1 X A * R -1 Y (Au k -b))
where the residual is ideally preconditioned. The resulting method is then refered as Ideal Minimal Residual (IMR) method. We have the following convergence result [MBF9, Proposition 4

.3]. Proposition 1.2.1. Assuming that γ(ρ) := max{|1-ρβ 2 |, |1-ρα 2 |} < 1/2, the sequence {u k } k≥1 defined by (1.2.12) is such that u k -u X ≤ (2γ) k u 0 -u X + 1 1 -2γ u -Π r (u) X (1.2.14)
and

lim sup k→∞ u k -u X ≤ 1 1 -2γ u -Π r (u) X (1.2.15)
The condition γ(ρ) < 1/2 imposes β α < √ 3 and ρ ∈ ( 1 2α 2 , 3 2β 2 ) which is a very restrictive condition, generally not satisfied without an excellent preconditioning of the operator A. Thus, for the ideal choice of norms (1.2.11), we have α = β = 1. Then , we have convergence for all ρ ∈ [ 1 2 , 3 2 ] towards a neighborhood of Π r (u) of size 2γ 1-2γ u-Π r (u) X with γ = |1-ρ|. Moreover, if ρ = 1, then u 1 ∈ Π r (u), the algorithm converges in one iteration for any initialization u 0 .

Approximate Ideal Minimal Residual method

The practical implementation of the IMR method rises some difficulties discussed as follows.

1. The computation of the iterate y k requires to solve an auxiliary problem that involves the Riesz map inverse

R -1 Y such that R Y = AR -1 X A * .
It is equivalent to compute the initial operator inverse A -1 . To overcome this difficulty, following [START_REF] Cohen | Adaptivity and variational stabilization for convection-diffusion equations[END_REF], we use instead Λ δ (r k ) an approximation of the residual r k = R -1 Y (Au kb) ∈ Y with a controlled relative precision δ > 0. Here, we have introduced Λ δ : Y → Y the approximation map such that

Λ δ (y) -y Y ≤ δ y Y , y ∈ R -1 Y (Av -b); v ∈ M ≤r (X) . (1.2.16)
2. Moreover, we may not be able to solve the best approximation problem in M ≤r (X) exactly, because there is no available algorithm for computing a global optimum, or because the algorithm has been stopped at a finite precision. Here, we introduce a set of quasi-optimal approximations

Π η r (u) ⊂ M ≤r (X) such that u -Π η r (u) X ≤ η u -Π r (u) X (η ≥ 1).
(1.2.17)

Using the two approximation maps defined by (1.2.16)-(1.2.17), we have proposed a perturbated version of the ideal gradient type algorithm amenable for numerical use. Letting u 0 = 0, we construct two sequences {u k } k≥0 ⊂ M ≤r and {y k } k≥0 ⊂ Y defined for k ≥ 0 by

y k = Λ δ (R -1 Y (Au k -b)) u k+1 ∈ Π η r (u k -R -1 X A * y k ) (1.2.18)
The resulting algorithm can be interpreted as a preconditioned gradient algorithm with an implicit preconditioner that approximates the ideal preconditioner. The algorithm works into steps. At the first step, an auxiliary problem provides an approximation of the residual R -1 Y (Au kb) up to precision δ. Then u k is computed as a quasi-optimal approximation in M ≤r (X) of the update u k -R -1 X A * y k . The resulting low-rank approximation approach is called Approximate Ideal Minimal Residual (A-IMR) method. In practice, these two steps are performed using practical low-rank approximation algorithms discussed previously as e.g. SVD, ALS or PGD. For a detailed presentation of the method see [MBF9, Section 5]. Again we have the following convergence result [MBF9, Proposition 5.3] for the perturbed gradient type algorithm.

Proposition 1.2.2. Assume (1.2.9), (1.2.16), and (1.2.17), with δ(1+η) < 1. Then, the sequence {u k } k≥1 defined by (1.2.18) is such that

u k -u X ≤ ((1 + η)δ) k u 0 -u X + η 1 -δ(1 + η) u -Π r (u) X . (1.2.19) Note that (1.2.19) implies that u -Π r (u) X ≤ u -u k X ≤ (1 + γ k ) u -Π r (u) X , (1.2.20)
This proposition states that the sequence {u k } k≥0 converges towards a neighborhood of the best approximation of u in M ≤r (X). The error depends on γ k such that lim sup k→∞ γ k ≤ c(η, δ) where c(η, δ) depends on the parameters δ and η that control the quality of this approximation.

Remark 1.2.3. We have proposed in a [MBF9, Section 6] a (weak) greedy algorithm for the adaptive construction of an approximation of the solution of problem (1.2.4). In particular, we use the A-IMR method for the computation of the greedy corrections in rank-1 tensor subset (see Equation (1.2.7)). A convergence proof of this algorithm is provided under some conditions on the parameters δ, η.

Numerical applications

We illustrate the behavior of the A-IMR method for the numerical solution of steady reactionadvection-diffusion stochastic problem on a two-dimensional unit square domain

D = [0, 1] 2 -∇(ξ) • (κ(ξ)∇u(ξ)) + c(ξ) • ∇u(ξ) + a(ξ)u(ξ) = f in D, (1.2.21) u = 0 on ∂D.
The solution u(ξ) : D → R depends on ξ = (ξ 1 , . . . , ξ p ) a random vector. Here a, κ are the reaction and diffusion coefficients respectively and advection field c. Combining a finite element discretization in space and stochastic Galerkin method for the parameters, the resulting Galerkin approximation still denoted by u, identified with its set of coefficients, is a tensor

u ∈ X = R n ⊗ (⊗ p µ=1 R s ).
Here, n and s are the dimensions of the finite element and piecewise polynomial space respectively. The tensor u is solution of the following algebraic tensor equation

Au = b, (1.2.22)
where the tensor operator A and tensor b are the algebraic representations of the bilinear and linear forms arising in the Galerkin formulation of the considered problem.

Two dimensional case with p = 1

We present numerical results concerning the A-IMR method applied to equation (1.2.22). This method provides an approximation ũ of the best approximation of u in M ≤r (X) with respect to • 2 the natural canonical norm on X. This choice corresponds to an operator R X = I n ⊗ I s , where I n (resp. I s ) is the identity in R n (resp. R s ). We take M ≤r (X) = T can r (X). For considered numerical applications, we have dim(X) ≥ 10 5 . The exact solution u of (1.2.22) is computed and used as a reference solution. For that choice, the • X -best approximation of u in M ≤r (X) is a rank-r SVD that can be computed exactly using classical algorithms (see Section 1.2.2). For the construction of an approximation in T can r (X) using A-IMR, we consider two strategies: the direct approximation in T can r (X) using the gradient type algorithm with M ≤r (X) = T can r (X), and a greedy algorithm that consists in a series of r corrections in T can 1 (X) computed using with the gradient type algorithm, with M ≤r (X) = T can 1 (X) and with an updated residual b at each correction. The A-IMR is compared to a standard approach, denoted CMR, which consists in minimizing the canonical norm of the residual. The latter approach has been introduced and analyzed in different papers, using either direct minimization or greedy rank-one algorithms [START_REF] Ammar | On the convergence of a greedy rank-one update algorithm for a class of linear systems[END_REF][START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF][START_REF] Doostan | A least-squares approximation of partial differential equations with high-dimensional random inputs[END_REF], and is known to suffer from ill-conditioning of the operator A.

The convergence curves with respect to the rank are shown in Figure 1.2.1, where the error is measured in the • 2 norm. Concerning the direct approach, we observe that the different algorithms have roughly the same rate of convergence. The A-IMR convergence curves are close to the optimal SVD for a wide range of values of δ. Concerning the greedy approach, we observe a significant difference between A-IMR and CMR. We note that A-IMR is close to the optimal SVD up to a certain rank (depending on δ) after which the convergence rate decreases but remains better than the one of CMR.

Higher dimensional case with p = 9

Here, we compute low-rank approximations of the solution of u in the canonical tensor subset T can r (X) with r ≥ 1. We rely on the greedy algorithm with successive corrections in M ≤r = T can 1 (X) computed with the gradient type algorithm. For considered numerical applications, we have dim(X) ≥ 10 11 . As a reference solution, we consider a low-rank approximation u ref of u computed using a greedy rank-one algorithm based on a canonical minimal residual formulation.

On Figure 1.2.2, we plot the convergence with the rank r of the approximations computed by both A-IMR and CMR algorithms and of the greedy approximations u of the reference solution u ref . We observe (as for the lower-dimensional example) that, with different values of the parameter δ (up to 0.9), the A-IMR method provides a better approximation of the solution in comparison to the CMR method. When decreasing δ, the proposed algorithm seems to provide approximations that tend to present the same convergence as the greedy approximations.

Summary

We have proposed a new algorithm for the construction of low-rank approximations of the solutions of high-dimensional weakly coercive problems formulated in a tensor space X. This algorithm is based on the approximate minimization of a particular residual norm on given low-rank tensor subsets M ≤r (X), the residual norm coinciding with some measure of the error in solution. Therefore, the algorithm is able to provide a quasi-best low-rank approximation with respect to a norm • X that can be designed for a certain objective. In view of goal oriented approximation, it is possible to design the norm • X such that it takes into account some particular quantity of interest, see first results in [START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF]. This topic has been studied more deeply in [START_REF] Zahm | Reduced Order Modeling for parameter dependent equations -Application to uncertainty quantification[END_REF]Chapter 5], but not considered for publication. However, it is discussed in Section 1.3 within the context of projection based low-rank approximation.

Let mention that the method has been applied with encouraging results for the solution of transient elastodynamics equations in space-time domain [START_REF] Boucinha | Ideal minimal residual-based proper generalized decomposition for non-symmetric multi-field models -Application to transient elastodynamics in space-time domain[END_REF]. Some extension to saddle points problems [START_REF] Billaud-Friess | A tensor-based algorithm for the model reduction of high dimensional problems[END_REF] with possible perturbed Uzawa algorithm have also been investigated.

Projection based goal oriented low-rank approximation

We recall briefly in Section 1.3.1 projection based low-rank approximation methods for the solution of parameter dependent linear problem under the form (1.1.1) formulated in Hilbert space. Then, in Section 1.3.2, we discuss possible goal oriented projection methods when the quantity to compute is not the solution itself but some Quantity of Interest (QoI) expressed as s = Lu, with L a linear map. From these two ingredients we derive a goal oriented reduced basis method. This work is related to the publication [START_REF] Zahm | Projection-Based Model Order Reduction Methods for the Estimation of Vector-Valued Variables of Interest[END_REF]. Some illustrative examples of the proposed approach are given in the last section.

Notations

Projection based low-rank approximation

Let V and W be two Hilbert spaces. We consider the linear parameter-dependent equation

A(ξ)u(ξ) = b(ξ) (1.3.1)
with A(ξ) ∈ L(V, W ) and b(ξ) ∈ W , where the operator A(ξ) and right-hand side b(ξ) depend on a parameter ξ which takes values in some parameter set Ξ ⊂ R p . We assume that A(ξ) is a norm-isomorphism such that for all u ∈ V ,

α(ξ) u V ≤ A(ξ)u W ≤ β(ξ) u V ,
where

inf 0 =v∈V sup 0 =w∈W A(ξ)v, w v V w W := α(ξ) > 0 and sup 0 =v∈V sup 0 =w∈W A(ξ)v, w v V w W := β(ξ) < ∞,
which ensures the well-posedness of (1.3.1).

We adopt here a subspace point of view for the low-rank approximation of the solution map u : Ξ → V . It relies on two key ingredients : properly chosen (small) subspaces of V , and appropriate projections on this subspaces. Here, we will focus in particularly on the well known Reduced Basis (RB) method, see monograph [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF], the surveys [START_REF] Haasdonk | Chapter II : Reduced basis methods for parametrized PDEs -A tutorial introduction for stationary and instationary problems[END_REF], and [START_REF] Nouy | Low-rank tensor methods for model order reduction[END_REF][START_REF] Nouy | Chapter IV : Low-rank methods for high-dimensional approximation and model order reduction[END_REF] for a presentation within low-rank approximation framework. RB method performs in two steps.

During an offline step, the subspaces are constructed from snapshots of the exact solution u(ξ) for a small number of values of the parameter ξ. Then the low-rank approximation defined as some projection of the solution u(ξ) in this subspace, is computed for any new instance of the parameter ξ ∈ Ξ during the online stage.

Online step: Galerkin type projections

In what follows, we consider that we are given a finite-dimensional subspace V r ⊂ V with dim(V r ) = r , the reduced space together with the associated reduced basis noted {v 1 , . . . , v r } ⊂ V .

Classical methods used for problems under the form (1.3.1) are Galerkin-type projections.

They provide an approximation u r (ξ) of u(ξ), in the dimensional subspace V r deduced from some orthogonality condition on the residual r(u(ξ), ξ) = A(ξ)u(ξ)b(ξ) with respect to a rdimensional space or by minimizing some residual norm. The resulting approximation u r ∈ V r is quasi-optimal that for all ξ ∈ Ξ i.e. there exists γ(ξ) ≥ 1

u(ξ) -u r (ξ) V ≤ γ(ξ) u(ξ) -P Vr u(ξ) V , (1.3.2)
were P Vr is the orthogonal projection of u in V r , defined as

u(ξ) -P Vr u(ξ) V = inf v∈Vr u(ξ) -v V .
We give here two typical examples of Galerkin projection methods.

Minimal residual based projection. The minimal residual projection u r is defined through

min v∈Vr r(v, ξ) W . (1.3.3)
In that case (1.3.2) holds with γ(ξ) = β(ξ) α(ξ) the condition number of the operator A(ξ). For problems badly conditioned, the minimal residual approximation method may lead to an approximation u r far from the best approximation in V r for the norm • V . In that context, it should be relevant to use parameter dependent preconditioner or ideal minimal residual based projection as discussed in Section 1.2.3.

Petrov-Galerkin projection. Let assume that the approximation u r ∈ V r is defined as a Petrov-Galerkin projection of u characterized by

r(u r (ξ), ξ), y = 0, ∀y ∈ W r , (1.3.4)
where W r ⊂ W is a test space of dimension r. Under the assumption that

α Vr,Wr (ξ) = inf 0 =v∈Vr sup 0 =y∈Wr A(ξ)v, y v V y W > 0, (1.3.5)
Céa's lemma [START_REF] Ern | Theory and practice of finite elements[END_REF]Lemma 2.8] states that (1.3.2) holds with γ(ξ) = 1 + β(ξ) α Vr ,Wr (ξ) .

The resulting low-rank approximation is a Galerkin projection of u(ξ) in V r given by

u r (ξ) = r i=1 α i (ξ)v i ,
where {α 1 , . . . , α r } are the coefficients of u r (ξ) in the reduced basis. If the residual r(u r (ξ), ξ) (through the operator A(ξ) and the right-hand side b(ξ)) admits low-rank representations, also called affine representation, it can be computed with a complexity only depending on r. If not, such representations can be obtained by the Empirical Interpolation Method (EIM) [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF]. It allows to compute the Galerkin projection u r (ξ) by solving a reduced system of r equations on the unknown coefficients {α i (ξ)} r i=1 with a computational complexity only depending on r.

Offline step: Construction of subspaces

At this step, the question is how to properly select the reduced space V r so that the solution manifold K := {u(ξ) : ξ ∈ Ξ} can be well approximated. To answer these question, we can use an error measure that gives a good indication on how well K can be approximated by a r-linear subspace. Assuming that the manifold K is compact, let us define

d ∞ r (K) V = min dim(Vr)=r sup ξ∈Ξ u(ξ) -P Vr u(ξ) V
which is the Kolmogorov r-width of K. This measure of error is a good indicator for uniformly accurate approximation over the parameter set Ξ. In particular, when the Kolmogorov r-width has a rapid decay with r, one can expect to obtain a good low-rank approximation. It is also relevant to consider subspaces optimal with respect to the L 2 -norm (or more generally L p -norm). In that case, we could consider

d 2 r (K) V = min dim(Vr)=r Ξ u(ξ) -P Vr u(ξ) 2 V dµ(ξ) 1/2
, where µ is the probability measure of ξ.

A simple choice is to take the subspaces V r generated from evaluations (snapshots) of the solution u(ξ) at for some a finite number of instances (possibly randomly chosen) Ξ = {ξ 1 , . . . , ξ K } ⊂ Ξ of the parameter ξ. But nothing guarantees that this choice provides reduced space close to the optimal r-dimensional subspaces. Thus, in practice V r is sought as the best r-dimensional subspace satisfying some optimality criterion related the measures

d 2 r (K) V or d ∞ r (K) V .
From samples. For the norm L 2 -norm, V r is sought such that

min dim(Vr)=r 1 K K k=1 u(ξ k ) -P Vr u(ξ k ) 2 V
where Ξ is a K-sample according to the probability measure µ. This choice is at the basis of POD methods for parameter-dependent equations [START_REF] Kahlbacher | Galerkin proper orthogonal decomposition methods for parameter dependent elliptic systems[END_REF][START_REF] Volkwein | Model reduction using proper orthogonal decomposition[END_REF]. For L ∞ -norm, V r is such that

min dim(Vr)=r max ξ∈ Ξ u(ξ) -P Vr u(ξ) V .
where Ξ is a set of K points in Ξ (called training set). In view of computation, a tractable definition of the subspaces relies on finding an appropriate set of parameters {ξ 1 , . . . , ξ r } ⊂ Ξ and impose V r = span{u(ξ 1 ), . . . , u(ξ r )}. Practically, greedy algorithm (e.g. [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Bebendorf | Comparison of some reduced representation approximations[END_REF] for (EIM) ) are prominent techniques for constructing such a nested sequence of reduced space from snapshots of the solution computed for adaptively chosen parameter values. It works as follows : given {ξ 1 , . . . , ξ r-1 } and the corresponding subspace

V r-1 = span{u(ξ 1 ), . . . , u(ξ r-1 )}, a new parameter value ξ r is selected such that u(ξ r ) -P Vr u(ξ r ) V = max ξ∈ Ξ u(ξ) -P V r-1 u(ξ) V ,
Then the reduced space is V r = span{u(ξ 1 ), . . . , u(ξ r )}.

From equations. As in the present case, u(ξ) is solution of the parameter-dependent equation (1.3.1), the quantities u(ξ) -P Vr u(ξ) V are not longer computable. Instead, a posteriori error estimates ∆(u r (ξ), ξ), computable without requiring u(ξ) and with low cost, are used. Typical choice consists in taking for ∆(u r (ξ), ξ) a bound of a certain norm of the residual r(u r (ξ), ξ) which is behind the PGD and RB methods. In the particular case of RB method, when ∆(u r (ξ), ξ) is certified and under the quasi-optimality of the Galerkin projection (1.3.2), the resulting algorithm is the so-called weak greedy algorithm. Convergence results on greedy algorithms can be found in [START_REF] Binev | Convergence Rates for Greedy Algorithms in Reduced Basis Methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF][START_REF] Devore | Greedy algorithms for reduced bases in Banach spaces[END_REF], especially they have proven that the approximation error in V r has the same type of decay as the Kolmogorov r-width for algebraic and exponential convergence.

Goal oriented projections

In many applications one is not interested in the solution u(ξ) itself, but only in a QoI s(ξ) which is a functional of u(ξ). Let Z be an Hilbert space, we consider

s(ξ) = L(ξ)u(ξ),
with L(ξ) ∈ L(V, Z). For example, for boundary value problems, L(ξ) can be defined as the trace operator providing the restriction of the solution to the boundary of the domain. In this case the QoI belongs to an infinite dimensional space or, after discretization, to a finite but possibly high dimensional space.

MOR for evaluation of QoI

Efficient goal-oriented methods have been proposed for the estimation of a scalar-valued variable of interest s(ξ). A standard method consists in computing an approximation of the solution of the so-called dual problem associated to (1.3.1) which is used to correct the estimation of s(ξ). We refer to [START_REF] Pierce | Adjoint Recovery of Superconvergent Functionals from PDE Approximations[END_REF] for a general survey on primal-dual methods and to [START_REF] Chen | Certified Reduced Basis Methods and Output Bounds for the Harmonic Maxwell's Equations[END_REF][START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF][START_REF] Haasdonk | Chapter II : Reduced basis methods for parametrized PDEs -A tutorial introduction for stationary and instationary problems[END_REF][START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF] for the application in the context of RB methods. Difficulties arise when the variable of interest s(ξ) takes its values in a vector space of finite or infinite dimension. The standard approach, which consists in treating s(ξ) as a collection of scalar-valued variables of interest and in building one reduced dual space for each of them, has a complexity which grows proportionally to the dimension of s(ξ). To circumvents this issue, one could construct a single reduced dual space, thus allowing to handle variables of interest with high and potentially infinite dimension. Such an approach can be found for parametric dynamical systems, see the monograph [START_REF] Benner | Dimension reduction of large-scale systems[END_REF] for a general introduction. In this framework, projection-based model order reduction methods are used for the approximation of s(ξ) which is an output of the dynamical system. Petrov-Galerkin methods have been proposed with different ways of constructing the reduced basis for the test and trial space, such as the balanced truncation methods, (balanced) POD method, moment matching methods, etc. We refer to [START_REF] Benner | Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems[END_REF] for a recent review on these methods.

In what follows, we present some work [START_REF] Zahm | Projection-Based Model Order Reduction Methods for the Estimation of Vector-Valued Variables of Interest[END_REF] concerning goal-oriented low-rank approximation for parameter dependent PDEs following the subsequent road map.

1. Develop an appropriate method for the approximation of variable of interest s(ξ) taking values in a single vector space of finite or infinite dimension.

2. Choose suitable projection principle (Petrov-Galerkin, Primal-Dual approach, Saddle point formulation) that guarantees the best approximation error on the QoI while using a single reduced space.

3. Using computable error estimate, propose greedy approaches to select the reduced spaces associated to the solution and the QoI.

Goal oriented projections

To alleviate the notations, we have omitted the dependence to ξ in this section. The approximation of any quantity a is noted ã.

A posteriori computation of QoI. We assume that ũ = u r ∈ V r , the Petrov-Galerkin projection of u by Equation (1.3.4) with W r satisfying the assumption (1.3.5). Then, we consider s as the standard approximation of s, computed a posteriori from ũ as

s = Lũ. (1.3.6)
The next proposition [MBF10, Proposition 2.1] provides an error bound for the standard approximation s.

Proposition 1.3.1. Let u = u r be the Petrov-Galerkin projection defined by (1.3.4), the approximation s of s defined by (1.3.10) satisfies

s -s Z ≤ δ L Wr 1 -(δ Vr,Wr ) 2 min v∈Vr u -v V , (1.3.7)
with δ Vr,Wr = max

0 =v∈Vr min y∈Wr v -R -1 V A * y V v V < 1, (1.3.8)
and

δ L Wr = sup 0 =z ∈Z min y∈Wr L * z -A * y V z Z . (1.3.9)
Primal-dual approach. Let us introduce the dual variable Q ∈ L(Z , W ) defined by the dual problem A * Q = L * . For given approximations u and Q, we define the approximation s by

s = L u + Q * (b -A u), (1.3.10)
where Q * (b -A u) is a correction using the approximation of the dual variable. The following proposition [MBF10, Proposition 2.5] provides an error bound on the QoI, which is a generalization of the classical error bound for scalar-valued variables of interest (see [START_REF] Pierce | Adjoint Recovery of Superconvergent Functionals from PDE Approximations[END_REF]) to vector-valued variables of interest.

Proposition 1.3.2. Let u = u r be the Petrov-Galerkin projection defined by (1.3.4), the approximation s of s defined by (1.3.10) satisfies

s -s Z ≤ L * -A * Q Z →V 1 -(δ Vr,Wr ) 2 min v∈Vr u -v V , (1.3.11)
where

L * -A * Q Z →V = sup 0 =z ∈Z (L * -A * Q)z V z Z . (1.3.12)
For the approximation

Q of Q ∈ L(Z , W ), the bound (1.3.11) suggests that L * -A * Q Z →V should be small. We then propose to choose Q as a solution of min Q∈L(Z ,W Q k ) L * -A * Q Z →V , (1.3.13)
where W Q k ⊂ W is a given approximation space (different from W r ). The next proposition [MBF10, Proposition 2.6] shows how to construct a solution of (1.3.13).

Proposition 1.3.3. The operator Q k : Z → W Q k defined for z ∈ Z by Q k z = arg min y k ∈W Q k L * z -A * y k V (1.3.14)
is linear and is a solution of (1.3.13).

Moreover Q k z ∈ W Q k is characterized by L * z -A * Q k z , R -1 V A * y k = 0, ∀y k ∈ W Q k . (1.3.15)
We give now a new bound [MBF10, Proposition 2.8] of the error on the QoI.

Proposition 1.3.4. Let u = u r be the Petrov-Galerkin projection defined by (1.3.4) and let Q = Q k be defined by (1.3.15). Then the approximation s defined by (1.3.10) satisfies

s -s Z ≤ δ L W Q k 1 -(δ Vr,Wr ) 2 min v∈Vr u -v V . (1.3.16)
where

δ L W Q k = sup 0 =z ∈Z min y∈W Q k L * z -A * y V z Z . (1.3.17)
For such choice of Q, the Proposition 1. Projection based on a saddle point problem In the line of [START_REF] Dahmen | Double greedy algorithms: Reduced basis methods for transport dominated problems[END_REF] for the approximation of (vectorvalued) variables of interest, we define the projection of u on the reduced space V r by means of a saddle point problem.

Let us equip W with a norm • W such that the relation y W = A * y V holds for any y ∈ W , which is equivalent to the following relation between the Riesz maps R W and R V :

R W = AR -1 V A * . (1.3.18)
We introduce a subspace T p ⊂ W of dimension p and we define the projection u r,p in V r as the solution of the saddle point problem

min v∈Vr max w∈Tp w W =1 Av -b, w . (1.3.19)
Under the condition (discrete inf-sup condition)

inf 0 =v∈Vr sup 0 =y∈Tp Av, y v V y W =: α Vr,Tp > 0, (1.3.20)
we have proven [MBF10, Proposition 2.10] that there exists a unique solution (u r,p , y r,p ) in

V r × T p to R W y r,p , y + Au r,p , y = b, y ∀y ∈ T p , (1.3.21a) A * y r,p , v = 0 ∀v ∈ V r . (1.3.21b)
Now, we define the approximation s as

s = Lũ + LR -1 V A * ỹ, (1.3.22)
where (ũ, ỹ) := (u r,p , y r,p ) ∈ V r × T p is the solution of the saddle point problem (1.

s -s Z ≤ δ L Tp 1 -(δ Vr,Tp ) 2 min v∈Vr u -v V . (1.3.23) with δ L Tp = sup 0 =z ∈Z min y∈Tp L * z -A * y V z Z (1.3.24)
and δ Vr,Tp = max

0 =v∈Vr min y∈Tp v -R -1 V A * y V v V . (1.3.25)
In practice, we will consider for T p spaces of the form

T p = W r + W Q k , (1.3.26) with dim(W r ) = r which implies δ L Tp ≤ δ L W Q k
and δ Vr,Tp ≤ δ Vr,Wr , so that the error bound (1.3.23) for the QoI is better than the error bound (1.3.16) of the primaldual method with primal approximation space V r , primal test space W r and dual approximation space W Q k . Therefore, we expect the approximation u r,p to be closer to the solution u than the Petrov-Galerkin projection u r . Also, the approximation of the QoI is expected to be improved.

Practical aspects

The projection methods introduced for the estimation of the QoI rely on the introduction of three spaces: the primal approximation space V r , the primal test space W r and the dual approximation space W Q k and for the saddle point approach the space T p = W r +W Q k . We adopt an offline/online strategy. Reduced (low-dimensional) spaces V r , W r and W Q k are constructed during the offline phase. Then, the projections on these reduced spaces and the evaluations of the QoI are rapidly computed for any parameter value ξ ∈ Ξ during the online phase.

Reduced space W r . Given the primal approximation space V r , we know that W r should be chosen such that δ Vr,Wr is as close to zero as possible (see Propositions 1.3.1, 1.3.4 and 1.3.5).

To that goal, we define the (parameter-dependent) test space as

W r (ξ) = P m (ξ) * R V V r , (1.3.27)
where P m (ξ) is an interpolation of the inverse of A(ξ) using m interpolation points in the parameter set Ξ [START_REF] Zahm | Interpolation of Inverse Operators for Preconditioning Parameter-Dependent Equations[END_REF]. The interpolation points are chosen as the points where solutions (primal and dual) have already been computed, i.e. the points given by (1.3.28).

Reduced spaces V r and W Q k . From suitable error estimates ∆(ξ) for the vector-valued QoI (see for details [MBF10, Section 3] and Section 1.3.3), different greedy algorithms for the construction of the reduced spaces V r and W Q k can be proposed. Given Ξ ⊂ Ξ a training set, they broadly consist in searching for a parameter value ξ * ∈ Ξ, a finite set where the error estimate ∆(ξ) is maximum, i.e.

ξ * ∈ arg max

ξ∈ Ξ ∆(ξ). (1.3.28)
A first strategy is to simultaneously enrich both the primal approximation space

V r+1 = V r + span(u(ξ * )) (1.3.29)
and the dual approximation space

W Q k+l = W Q k + range(Q(ξ * )) (1.3.30)
at each iteration. This strategy is referred as the simultaneous construction, as opposed to the alternate construction which consists in enriching

W Q k (resp. V r ) if V r (resp. W Q k )
were enriched at the previous greedy iteration step.

Numerical applications

We consider the benchmark problem of the cooling of electronic components [START_REF] Zahm | Interpolation of Inverse Operators for Preconditioning Parameter-Dependent Equations[END_REF]. The equation to solve is an advection-diffusion equation over the domain

D ⊂ R 2 -∇ • (κ(ξ)∇T (ξ)) + c(ξ) • ∇T (ξ) = f, (1.3.31)
whose solution T (ξ) : D → R is the temperature field. Here κ and c denote respectively the parameter-dependent diffusion coefficient and the advection field. Here, we focus on the resulting algebraic parameter-dependent equation coming from stabilized finite element discretization of (1.3.31), that is A(ξ)u(ξ) = b(ξ), where u(ξ) ∈ R n are the coefficients of the finite element approximation T h = n i=1 u i ϕ i of T , and where ξ = (ξ 1 , . . . , ξ 4 ) is a 4-dimensional random vector. The space V = R n with n = 2.8 × 10 4 is endowed with the norm

• V = • V 0 which corresponds to the H 1 (D)-norm 1 . The variable of interest s(ξ) = (s 1 (ξ), s 2 (ξ))
is the mean temperature of both electronic components, with

s 1 (ξ) = 1 |D IC 1 | D IC 1 T h (ξ) , s 2 (ξ) = 1 |D IC 2 | D IC 2 T h (ξ), (1.3.32)
where D IC i (i = 1, 2) are two subdomains of D. Then we can write s(ξ) = Lu(ξ) for an appropriate L ∈ R l×n , with l = 2. Here we have Z = R 2 , which we equip with the canonical norm on R 2 .

Comparison of the projections methods Here the approximation spaces V r , W Q k and the test space W r are given. We denote by V r , W Q k and W r the matrices containing the basis vectors of the corresponding subspaces. The matrix V r = (u(ξ 1 ), . . . , u(ξ 50 )) contains 50 snapshots of the primal solution (r = 50), and

W Q k = (Q(ξ 1 ) . . . Q(ξ 25
)) contains 25 snapshots of the dual solution so that the dimension of

W Q k is k = 25l = 50.
The test space W r is defined according to (1.3.27). The matrix associated to the test space is given by W r (ξ) = P T m (ξ)R V V r . We consider a samples set Ξ t of size t = 10 4 . For any ξ ∈ Ξ t we compute the exact quantity of interest s(ξ) and the approximation s(ξ) by the following methods.

• Primal only: solve the linear system

W T r (ξ)A(ξ)V r U r (ξ) = W r (ξ) T b of size r and compute s(ξ) = LV r U r (ξ).
• Dual only: solve the linear system

(W Q k ) T A(ξ)R -1 V A(ξ) * W Q k Y k (ξ) = (W Q k ) T b of size k and compute s(ξ) = LR -1 V A(ξ) * W Q k Y k (ξ).
• Primal-dual: solve the linear system of the Primal only method, solve the linear system

(W Q k ) T A(ξ)R -1 V A(ξ) * W Q k Y k (ξ) = (W Q k ) T b -(W Q k ) T A(ξ)V r U r (ξ)
of size k, and compute

s(ξ) = LV r U r (ξ) + LR -1 V A(ξ) * W Q k Y k (ξ).
• Saddle point: solve the linear system of size p + r

T T p (ξ)A(ξ)R -1 V A(ξ) * T p (ξ) T T p (ξ)A(ξ)V r T T p (ξ)A(ξ)V r T 0 Y r,p (ξ) U r,p (ξ) = T p (ξ) T b 0 with T p (ξ) = W r (ξ), W Q k , and compute s(ξ) = LV r U r,p (ξ) + LR -1 V A(ξ) * T p (ξ) Y r,p (ξ).
The numerical results are given in Figure 1.3.1. The saddle point method leads to the lowest error on the variable of interest. Also, we see that a good preconditioner (for example with m = 30) improves the accuracy for the saddle point method, the primal only method and the primal-dual method. However, this improvement is not really significant for the considered application: the errors are barely divided by 2 compared to the non preconditioned Galerkin projection (m = 0). In fact, the preconditioner improves the quality of the test space, and the choice W r = V r (yielding the standard Galerkin projection) is sufficiently accurate for this example and for the chosen norm on V .

Quality of error estimate

We discuss now the quality of the error estimate ∆(ξ) for the variable of interest. We consider

∆(ξ) = P m (ξ)(A(ξ)u r (ξ) -b(ξ)) V 0 L(ξ) * -A(ξ) * Q k (ξ) Z →V 0 (1.3.33)
for the primal-dual method, and

∆(ξ) = P m (ξ)(A(ξ)t r,p (ξ) -b(ξ)) V 0 sup 0 =z ∈Z inf y∈Tp L(ξ) * z -A(ξ) * y V 0 z Z (1.3.34)
for the saddle point method. We see that the max-min ratio and the normalized standard deviation are decreasing with m: this indicates an improvement of the error estimate. Furthermore, the mean value of η(ξ) seems to converge (with m) to 19.5 for the primal-dual method, and to 13.8 for the saddle point method. In fact, with a good preconditioner,

P m (ξ)(A(ξ)u r (ξ) -b(ξ)) V 0 (or P m (ξ)(A(ξ)t r,p (ξ) -b(ξ)) V 0
) is expected to be a good approximation of the primal error u(ξ)u r (ξ) V 0 (or u(ξ)t r,p (ξ) V 0 ), but this does not ensure that the effectivity index η(ξ) will converge to 1.

Summary

We have proposed and analyzed projection based methods for the estimation of vector-valued variables of interest in the context of parameter-dependent equations. This includes a generalization of the classical primal-dual method to the case of vector-valued variables of interest, and also a Petrov-Galerkin method based on a saddle point problem.

The error estimates proposed and used during the greedy enrichment of the reduced spaces, involve the use of Cauchy-Schwarz inequalities, which are clearly not optimal. To overcome this limitation, randomized a posteriori error estimator for low-rank approximations, reliable and efficient at given high probability have been proposed recently in [START_REF] Smetana | Randomized Residual-Based Error Estimators for Parametrized Equations[END_REF]. Chapter 2

E(η(ξ)) max η(ξ) min η(ξ) Var(η(ξ)) 1/2 E(η(ξ)) Primal-dual m = 0 5.

Dynamical low-rank approximation

In this chapter, model order reduction methods introduced previously are extended for the solution of parameter and time-dependent problems. Among these methods, we particularly focus on dynamical low-rank approximation approaches.

After a brief discussion on possible low-rank approximation methods for model order reduction of time-dependent and parameter-dependent problems in Section 2.1, we present in Section 2.2 dynamical reduced basis method that can be interpreted as a dynamical low-rank approximation method with a subspace point of view. Then, in Section 2.3, we present a new splitting integration scheme suitable for dynamical low-rank approximation in the set of fixed rank matrices.

General context

Here, we focus on parameter and time-dependent problems. Typical problems could either be ordinary differential equations (ODEs), or partial differential equations (PDEs) depending on some (random) parameters. When PDEs are considered, we apply a method of lines approach, that is, we first discretize in space such that the problem is approximated by a system of ODEs. Thus, we are concerned with the solution of parameter-dependent non-autonomous dynamical systems of the form

u (t, ξ) = f (u(t, ξ), t, ξ), t ∈ (0, T ], u(0, ξ) = u 0 (ξ), (2.1.1)
where the flux f and initial condition u 0 depend on some (random) parameters ξ = (ξ 1 , . . . , ξ p ) T with values in a parameter set Ξ ⊂ R p . We denote T > 0 the final time. We assume that the solution u(t, ξ) belongs to some vector space V for all t ∈ [0, T ].

Low-rank approximation for parameter and time-dependent problems

In Chapter 1, we have seen that low-rank approximation methods are effective approaches for solving time-independent parameter-dependent problems (1.1.1). Seeking a solution u to the problem (2.1.1) leads to similar difficulties in term of complexity and computation costs. Here, we discuss how to extend such low-rank approximation methods to compute an approximation of the solution map ξ → u(t, ξ) for all t ∈ [0, T ].

Low-rank projection based approaches

In the context of RB methods (see, e.g., [START_REF] Haasdonk | Chapter II : Reduced basis methods for parametrized PDEs -A tutorial introduction for stationary and instationary problems[END_REF][START_REF] Ohlberger | Reduced Basis Methods: Success, Limitations and Future Challenges[END_REF]), the reduced approximation is classically obtained as the (Petrov-)Galerkin projection of the solution u of (2.1.1) onto a time-independent sub-space V r of V . It writes under following form

u(t, ξ) ≈ r i=1 v i α i (t, ξ), (2.1.2)
where {v 1 , . . . , v r } is a time-independent basis of V r and {α 1 (t), . . . , α r (t)} are elements of some vector space S of functions defined on Ξ. Such an approximation, can be seen as a rank-r approximation in the tensor space V ⊗ S at each time t. If the flux f admits a low-rank representation, the approximation (2.1.2) can be computed efficiently by solving a reduced dynamical system of r equations satisfied by the coefficients {α i (ξ)} r i=1 , using e.g. a time-stepping scheme. Sample-based methods have been proposed for the construction of V r (see, e.g., [START_REF] Baur | Comparison of Methods for Parametric Model Order Reduction of Instationary Problems[END_REF][START_REF] Bui | Model reduction for large-scale systems with highdimensional parametric input space[END_REF][START_REF] Janon | Certified reduced-basis solutions of viscous Burgers equation parametrized by initial and boundary values[END_REF][START_REF] Nguyen | Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers' equation[END_REF]). In [START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF], V r is obtained as the span of snapshots u(t k , ξ i ) (in both time and parameter) of the solution of the full-order model. However, for a high-dimensional space V , it is not feasible (and far from optimal) to retain a large number of snapshots. Then, one can rely on a POD of the snapshots matrix in order to extract subspaces which are optimal in a mean-square sense. In particular, adaptive construction of V r can be performed through a POD-greedy algorithm using a posteriori error estimates [START_REF] Drohmann | Reduced Basis Approximation for Nonlinear Parametrized Evolution Equations based on Empirical Operator Interpolation[END_REF][START_REF] Eftang | An hp certified reduced basis method for parametrized parabolic partial differential equations[END_REF][START_REF] Haasdonk | Reduced basis method for finite volume approximations of parametrized linear evolution equations[END_REF][START_REF] Haasdonk | Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition[END_REF][START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF][START_REF] Wirtz | A-posteriori error estimation for DEIM reduced nonlinear dynamical systems[END_REF]. As pointed in [START_REF] Ohlberger | Reduced Basis Methods: Success, Limitations and Future Challenges[END_REF], efficiency of such methods is limited by the effectivity of error estimators that may show bad long-time effectivity. Moreover, for parameter and time dependent problems with complex dynamics (e.g. advection dominated or hyperbolic parameter-dependent equations), the optimal subspace on which to approximate the solution at each time instant can considerably change over time. In that context, classical RB methods may require very high-dimensional (time-independent) reduced spaces for computing an approximation given by (2.1.2) with prescribed accuracy. PGD method has also been considered in [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Proper Generalized Decompositions and Separated Representations for the Numerical Solution of High Dimensional Stochastic Problems[END_REF] for the low-rank approximation of the solution of stochastic evolution equations. In this context, the approximation is sought under the form

u(t, ξ) ≈ r i=1 v i (t)α i (ξ), (2.1.3)
where {v 1 (t), . . . , v r (t)} is a time-dependent basis of V r (t) and {α 1 , . . . , α r } are elements of S. Approximation (2.1.3) is a rank-r element in the tensor space V [0,T ] ⊗ S. The resulting approximation can be seen as a projection of u(t, ξ) ∈ V onto a time-dependent sub-space V r (t) = span{v 1 (t), . . . , v r (t)} ⊂ V which allows to well capture transient phenomena. Contrary to RB basis method previously discussed, the approximation in (2.1.3) is obtained with a global in time variational principle that is not optimal at each instant t. Let us mention, that reduced basis method based on Petrov-Galerkin space-time (PG-ST) formulations has also been introduced, see for example [START_REF] Steih | Space-time Reduced Basis methods for time-periodic parametric partial differential equations[END_REF][START_REF] Urban | A new error bound for reduced basis approximation of parabolic partial differential equations[END_REF]. Such an approach provides a low-rank approximation of the form (2.1.3). At the discrete level, previous approaches differ from RB approaches since they do not rely on time-stepping scheme except in very particular cases [START_REF] Glas | Two Ways to Treat Time in Reduced Basis Methods[END_REF] .

Dynamical low-rank approximation methods

Dynamical low-rank approximation (DLRA) methods (see survey [161, Section 9.5]) aim at computing an approximation of u under the form

u(t, ξ) ≈ r i=1 v i (t)α i (t, ξ). (2.1.4)
The approximation (2.1.4) can be interpreted as a projection onto a time-dependent reduced space V r (t) obtained through principles which are local in time (e.g. Dirac Frenkel principle). DLRA methods are particular low-rank tensor methods since they rely on computing a rank-r approximation in the tensor space V ⊗ S at each time t. They first appeared in [START_REF] Koch | Dynamical low-rank approximation[END_REF][START_REF] Nonnenmacher | Dynamical low-rank approximation: applications and numerical experiments[END_REF] in the case of low-rank matrices and then extended to low-rank approximation in tensor format (see, e.g., [START_REF] Lubich | Dynamical approximation by hierarchical tucker and tensor-train tensors[END_REF][START_REF] Uschmajew | The geometry of algorithms using hierarchical tensors[END_REF]) for the solution of high dimensional time dependent problems. Let us mention that dynamically (bi-)orthogonal methods that provide approximations under the form (2.1.4), have been also proposed [START_REF] Cheng | A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations I: Derivation and algorithms[END_REF][START_REF] Cheng | A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations II: Adaptivity and generalizations[END_REF][START_REF] Feppon | Dynamically Orthogonal Numerical Schemes for Efficient Stochastic Advection and Lagrangian Transport[END_REF][START_REF] Sapsis | Dynamically orthogonal field equations for continuous stochastic dynamical systems[END_REF] and analysed [START_REF] Musharbash | On the dynamically orthogonal approximation of time-dependent random PDEs[END_REF] for the numerical solution of timedependent stochastic equations.

Contributions

The contributions [MBF11,MBF13,MBF14,MBF16] presented in the next sections are the result of collaborations with A. Falcó (Univ. CEU Cardenal Herrera) and A. Nouy (ECN, LMJL) on dynamical low-rank approximation methods for the solution of parameter dependent dynamical systems and T. Heuzé (ECN,GeM) for possible extension to hyperbolic parameter dependent conservation laws.

• We have proposed a dynamical RB method for the solution of parameter-dependent dynamical systems of the form (2.1.1). It relies on the greedy construction of time-dependent sub-spaces V r (t) ⊂ V generated from samples of the solution of the full-order model at some selected parameter values. The resulting method can be interpreted as a DLRA method with a subspace point of view and a uniform control of the error over the parameter set. Some prospects for extension of such dynamical RB to parameter dependent hyperbolic conservation laws, within finite volume framework, are also discussed.

• Regarding to matrix ODE (e.g., arising from (2.1.1) after discretization of the parameter domain), the DLRA relies on approximation in the set of fixed rank matrices. However, the representation of any element of this set is not unique. To dodge this undesirable property, we have proposed a geometric description of the set of fixed rank matrices endowed with the structure of analytic principal bundle, relying on an explicit description of local charts. From this description, we have derived a new splitting integration scheme suitable for DLRA.

The work concerning dynamical RB methods has been partially supported by the GdR MOMAS through the project REDYN (2015). The work concerning application to hyperbolic conservation laws has been founded through the project PEPS : DROME by the Cellule Energie du CNRS (2019), co-driven with T. Heuzé.

Projection based dynamical low-rank approximation

Notations. We consider in this section the resolution of a dynamical system under the form (2.1.1) in algebraic setting i.e. V = R N . In what follows, V is equipped with canonical inner product u, v := u T v and associated norm u .

The dynamical RB approach developed in [START_REF] Billaud-Friess | Dynamical Model Reduction Method for Solving Parameter-Dependent Dynamical Systems[END_REF] performs in two steps. During an offline step, the time-dependent subspaces are constructed in greedy fashion from the trajectories t → u(t, ξ) of the exact solution for the selected values of the parameter ξ. Then the lowrank approximation is defined as the Galerkin projection of the solution u(t, ξ) in the time dependent subspace, and computed for any new instance of the parameter ξ ∈ Ξ during the online stage. For the sake of brevity, we only detail the projection step and greedy construction of the time-dependent reduced spaces.

Construction of the approximation

We consider that for any t > 0, the reduced space V r (t) ⊂ V is given and satisfies dim V r (t) = r for all t > 0.

(2.2.1)

We denote by {v 1 (t), . . . , v r (t)} the an orthonormal basis of V r (t) and V (t) = [v 1 (t), . . . , v r (t)] ∈ R N ×r the associated orthogonal matrix. Finally Π Vr (t) = V (t)V T (t) ∈ R d×d is the orthogonal projector onto V r (t). We seek an approximation u r (t, ξ) ∈ V r (t) of u(t, ξ) under the form

u r (t, ξ) := r i=1 α i (t, ξ)v i (t) for all t > 0. (2.2.2)
where α r = (α 1 , . . . , α r ) T is solution of the small r-order system obtained through the Galerkin projection of (2.1.1) onto V r (t)

α r (t, ξ) = f r (α r (t, ξ), t, ξ), α r (0, ξ) = α 0 r (ξ), (2.2.3)
with α 0 r (ξ) = V T (0)u 0 (ξ) and a reduced flux f r defined by for all t ≥ 0, where ∆ r (t, ξ) is the solution of the ordinary differential equation

f r (α r , t, ξ) = V (t) T f (V (t)α r , t, ξ) -V (t) T V (t)α r . ( 2
∆ r (t, ξ) = L[f ](u r (t, ξ))∆ r (t, ξ) + r(t, ξ) , ∆ r (0, ξ) = u 0 (ξ) -Π Vr (0)u 0 (ξ) , (2.2.6) with r(t, ξ) = Π V ⊥ r (t)(V (t)V (t) T u r (t, ξ) -f (u r (t, ξ), t, ξ))
, the log-lipschitz constant defined as

L[f ](v) = sup u∈V,u =v u -v, f (u) -f (v) u -v 2 .
Remark 2.2.2. Note that the last term in (2.2.3) takes into account the time dependency of the subspace V r , otherwise V (t) = 0 and we recover the usual projected dynamical system within the framework of classical RB method (see, e.g., [START_REF] Haasdonk | Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition[END_REF]). Remark 2.2.3. In practice, the coefficients α r as well as the error bound ∆ r are computed using an appropriate integration scheme. The resulting numerical method can be implemented in an efficient time-stepping offline/online fashion provided suitable assumptions on the flux (see [MBF11, Section 4] for details).

Selection of sub-spaces

We present greedy algorithms for the construction of an increasing sequence of reduced spaces {V r } r>0 . These spaces are generated from successive samples of the solution u of the full-order model at parameters values {ξ r } r>0 which are selected adaptively using the an a posteriori estimate ∆ r among a finite training set Ξ ⊂ Ξ.

T-Greedy We have proposed in [START_REF] Billaud-Friess | Dynamical Model Reduction Method for Solving Parameter-Dependent Dynamical Systems[END_REF], an adaptive greedy strategy called T-greedy algorithm. A step r of this algorithm is as follows:

ξ r+1 ∈ arg max ξ∈ Ξ ∆ (0,T ) r (ξ), (2.2.7) V r+1 (t) = V r (t) + span{u(t, ξ r+1 )}. (2.2.8)
with the following global error bound ∆ (0,T ) r (ξ) := ∆ r (•, ξ) (0,T ),p , where • (0,T ),q denotes the natural norm in L q (0, T ) with q = 2, ∞.

POD-Greedy For the adaptive construction of time-independent reduced spaces V r , we recall the POD-greedy algorithm [START_REF] Haasdonk | Reduced basis method for finite volume approximations of parametrized linear evolution equations[END_REF]. It selects a new parameters value

ξ r+1 ∈ arg max ξ∈ Ξ ∆ (0,T ) r (ξ).
then the trajectory u(•, ξ r+1 ) of the full-order model is computed and the space V r is enriched by the -dimensional subspace generated by the first POD modes of u(•, ξ r+1 ) -Π Vr u(•, ξ r+1 ).

Remark 2.2.4. If the error estimate is such that

c r ∆ r (t, •) ≤ u(t, •) -u r (t, •) ≤ C r ∆ r (t, •),
(2.2.9)

with 0 < c r ≤ C r < ∞ some positive constants, the POD-greedy algorithm has been proved to converge with quasi-optimal rates [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF].

Advantages of the proposed approach

• For time-dependent subspace V r (t), the Galerkin projection u r of the reduced dynamical system interpolates the solution u for the instances of the parameters ξ 1 , . . . , ξ r selected with the T-greedy algorithm. Whereas, for time-independent subspace V r obtained through the POD-greedy algorithm, the approximation u r does not in general interpolate the function ξ → u(•, ξ) for ξ 1 , . . . , ξ r .

• POD-greedy algorithm may lead to high-dimensional reduced spaces for reaching a desired accuracy. A typical example is the advection problem for which a very high-dimensional time-independent reduced space may be required to approximate the solution even for one instance of the parameter. We will see on numerical examples (see Section 2.2.3), that using time-dependent reduced spaces allows to overcome such a limitation.

Numerical applications

The proposed approach is applied for the numerical approximation of the solution of the following nonlinear partial differential equation defined on a open interval D ⊂ R and a time interval

I = (0, T ) ∂ ∂t u(ξ) + ∂ ∂x (c(u(ξ), ξ)u(ξ)) + a(ξ) ∂ ∂x u(ξ) -µ(ξ) ∂ 2 ∂x 2 u(ξ) = g(ξ), on D × I, (2.2.10)
with appropriate boundary conditions and a parameter-independent initial condition u 0 (x).

Here ξ = (ξ 1 , . . . , ξ d ) denotes a random vector with values in Ξ and with independent components. The functions c : X × Ξ → R and a : D × Ξ → R will be specified in each test case. Finally, µ : Ξ → R is a parameter-dependent coefficient and g : D × I × Ξ → R is a given source term. We consider an approximation of the solution of (2.2.10), still noted u, obtained with an appropriate scheme yielding a system of N ordinary differential equations of the form (2.1.

1).

In what follows, we denote by MTI (resp. MTD) the method using time-independent (resp. time-dependent) reduced spaces.

Advection equation

Let D = (0, 1) and I = (0, 0.2). The first application is an advection equation µ(ξ) = 0, g(x, t, ξ) = 0, c(u, ξ) = 0 and a(x, ξ) = a(ξ) = a 0 + a 1 ξ 1 , with ξ 1 ∼ U (-1, 1), a 0 = 1 and a 1 = 0.5. We impose periodic boundary conditions and consider a finite difference upwind scheme with N = 2000, together with an explicit Euler time integration scheme. We consider two functions as initial conditions. One is a smooth function given by u

0 cont (x) = 1 √ 2π exp(-x-0.6 0.05 2 )
whereas the second is a discontinuous function given by u 0 disc (x) = 1 [0.1,09] (x)( 3x + sin(10x)).

Deterministic case

We first consider a deterministic problem with ξ = 0.65 and compare the approximations obtained by projections on reduced spaces constructed in two different ways.

In the first method (MTI), the reduced space V r is time-independent and generated by the first r modes of the POD of the trajectory t → u(t, 0.65), with 1 ≤ r ≤ 200. In the second method (MTD), we consider the one-dimensional time-dependent space V r (t) = span{u(t, 0.65)}, r = 1.

For obtaining a very accurate precision, we have shown that the MTI requires a reduced space with rather high dimension (i.e. 50) when u 0 = u 0 cont . This is worst when u 0 = u 0 disc where a reduced space with higher dimension (i.e. 200) is required to well approximate the discontinuous solution. Concerning the MTD, we always obtain relative errors up to the machine precision, in both cases, with only one-dimensional time dependent reduced space.

General case

We now consider the parameter-dependent problem and compare the approximations obtained with MTI and MTD for a subspace of dimension r = 20 generated from a training set of size 30. We estimate the expectation E(E q (ξ)) and maximum max ξ∈Ξ (E q (ξ)) of the relative errors E q = u ru I,q / u I,q , with q = 2, ∞, by the empirical mean and by the maximum of the values of E q taken at 50 randomly chosen values of the parameters. These quantities are depicted on Figure 2.2.1 for different values of r and for both continuous and discontinuous initial conditions. For the same dimension of the reduced spaces, MTD clearly provides a more accurate approximation than MTI in particular when considering discontinuous initial condition.

10 0 rank r max(E 2 ) max(E ∞ ) E(E 2 ) E(E ∞ )

Viscous Burger's equation

We study a nonlinear viscous Burger's equation with uncertain parameters [START_REF] Wirtz | A-posteriori error estimation for DEIM reduced nonlinear dynamical systems[END_REF]. We consider a spatial domain D = (0, 1), a time interval I = (0, 1), homogeneous Dirichlet boundary conditions and an initial condition u 0 = 0. We consider c(u, ξ) = 1 2 u, a(x, ξ) = 0 and a diffusion coefficient µ(ξ) = ξ, with ξ ∼ U (0.01, 0.06) a uniform random variable. The source term is defined by g(x, t, ξ) = g 1 (x, t) + g 2 (x, t), with g 1 (x, t) = 4e -( x-0.2 0.03 ) 2 1 [0.1,0.3] (x) sin(4πt) and g 2 (x, t) = 4 • 1 [0.6,0.7] (x)1 [0.2,0.4] (t). We use a finite difference scheme in space (N = 298) and a semi-implicit Euler scheme in time (with K = 200 time steps). Reduced spaces are constructed with POD-greedy and T-greedy algorithms with a training set of size 60. The superiority of MTD over MTI is confirmed by Figure 2.2.2 where we have plotted statistical estimations (with a random sample of size 50) of the expectation and maximum of the relative errors. As we can see, MTD provides with r = 15 an approximation with a relative error of 10 -10 , whereas MTI only provides an approximation with a relative error of 10 -2 for the same dimension of the reduced space. For r = 50, MTI provides an approximation with relative error 10 -6 , which is still higher than MTD with r = 15.

10 0 rank r max(E 2 ) max(E ∞ ) E(E 2 ) E(E ∞ )

Prospects for parameter-dependent conservation laws

In this section, we give some elements of discussion concerning the extension of the proposed dynamical RB for parameter-dependent hyperbolic conservation laws. Let us consider the following one dimensional parameter-dependent scalar conservation equation. The space domain is the bounded open interval D ⊂ R, and I = [0, T ] is the time domain. The function u(•, t; ξ), belonging to the space V , is solution of the equation

∂ t u(x, t, ξ) + ∂ x f (u(x, t, ξ), ξ) = 0, (x, t) ∈ D × I, (2.2.11)
with initial condition u 0 : D → R and suitable boundary conditions. The considered initial boundary value problem (2.2.11) depends on parameters ξ ∈ Ξ ⊂ R p through the conservative flux f (•, ξ).

Context and challenges

It is now well understood that classical low-rank approximation methods are linear approximation methods which are not suitable for solving (parameter-dependent) hyperbolic problems because the solution manifold can not be well approximated with a linear space of small dimension r. Indeed, for such problems, the Kolmogorov r-width decreases slowly with r (see, e.g., [START_REF] Greif | Decay of the Kolmogorov N-width for wave problems[END_REF][START_REF] Welper | Interpolation of Functions with Parameter Dependent Jumps by Transformed Snapshots[END_REF]). To overcome this issue, new ROM approaches have emerged this past years. A seminal work for dealing with transport dominated problem has been considered in [START_REF] Ohlberger | Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing[END_REF] using a freezing method. It opened the route to various MOR approaches most of them relying on constructing suitable approximations from transformed snapshots see, e.g., [START_REF] Black | Projection-based model reduction with dynamically transformed modes[END_REF][START_REF] Cagniart | Model Order Reduction for Hyperbolic Problems: a new framework[END_REF][START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Reiss | The Shifted Proper Orthogonal Decomposition: A Mode Decomposition for Multiple Transport Phenomena[END_REF][START_REF] Rim | Manifold Approximations via Transported Subspaces: Model reduction for transport-dominated problems[END_REF][START_REF] Welper | Interpolation of Functions with Parameter Dependent Jumps by Transformed Snapshots[END_REF].

Moreover, defining suitable projection based ROM methods providing approximations that preserve as much as possible the mathematical features of solutions (conservative, entropic, monotonic, total variation dimininishing...) of hyperbolic PDEs remain an opened question. In this direction, it was proposed in [START_REF] Abgrall | Crisovan Robust model reduction by L 1 -norm minimization and approximation via dictionaries: application to nonlinear hyperbolic problems[END_REF][START_REF] Abgrall | Crisovan Model reduction using L 1 -norm minimization as an application to nonlinear hyperbolic problems[END_REF] an approximation problem based on the minimization of the residual of the discretized equations in L 1 -norm for ROM of hyperbolic conservation laws. It has the advantage to provide non oscillatory approximation, especially in presence of shocks. Projection-based hyper-reduced models of nonlinear conservation laws globally conservative that inherits a semi-discrete entropy inequality has also been proposed in [START_REF] Chan | Entropy stable reduced order modeling of nonlinear conservation laws[END_REF].

Transformed dynamical reduced basis

As illustrated in Section 2.2.3, the dynamic RB method discussed previously has shown to be interesting candidate when dealing with parameter-dependent transport dominated problems. We have considered, within the project PEPS DROME, to combine this approach with adapted local basis, with respect to the parameter, for solving Equation (2.2.11). The idea is to design reduced spaces spanned by local basis functions {φ i (x, t; ξ)} r i=1 deduced from {u(ϕ i (x, t, ξ), t), t; ξ i )} r i=1 , which correspond to transformed snapshots of the solution for given parameter instances with ϕ(•, t, ξ, ξ i ) : D → D a parameter-dependent space transformation at the instant t. Then, the solution u(x, t, ξ) of Equation(2.2.11) is approximated by the following rank-r approximation u(x, t, ξ) ≈ r i=1 α i (t, ξ)φ i (x, t; ξ).

(2.2.12)

From a physical viewpoint, a relevant choice is to derive the parameter-dependent transformations ϕ(ξ, ξ i ) from the characteristic associated with the hyperbolic system (2.2.11), such that the approximation (2.2.12) captures well the features of the original solution. Such approximation is exact with r = 1 for parameter-dependent linear transport equation when the space transformation ϕ i (ξ, ξ i ) is a suitable chosen parameter and time dependent space shift (see, e.g., [START_REF] Rim | Manifold Approximations via Transported Subspaces: Model reduction for transport-dominated problems[END_REF]Example 3.6]). Thus, efficiently computing the transformed snapshots of the trajectory of the true solution for given instances of the parameter is crucial in that context.

Reconstruction of finite volume solution from snapshots

In the paper [START_REF] Billaud-Friess | Reconstruction of finite volume solution for parameter-dependent linear hyperbolic conservation laws[END_REF], we have focused on the design of a robust approximation for discontinuous solutions generated by parameter-dependent linear hyperbolic systems. The proposed approach is derived in Finite Volume (FV) framework where the snapshots consist of trajectories of the numerical approximation provided by a known FV scheme. To this end, the Reconstruct-Translate-Average (RTA) method have been introduced, which is inspired from Godunov's method also interpreted as Reconstruct-Evolve-Average (REA) [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. The proposed method is detailed and analyzed for one dimensional parameter-dependent linear scalar transport equation, and an application to system of linear hyperbolic conservation laws is shown. In that case, one snapshot is sufficient to approximate the FV solution for any parameter values. This leads to a simple and efficient method that does not require any projection step neither time stepping procedure. The proposed approach has been applied for the reconstruction of solutions for the parameter-dependent transport problem and the wave equation.

Summary and discussion for future research

We have discussed previously projection-based model order reduction approach for the solution of parameter-dependent dynamical systems. This generalizes classical reduced basis methods with time-independent reduced spaces and performs better in practice. In particular, the proposed T-Greedy provides smaller reduced spaces V r (t) allowing reduced computational cost during the online stage.

Here are some some perspectives and questions that follow from the presented contributions.

Towards nonlinear approximation.

• Obviously the RTA approach proposed in [START_REF] Billaud-Friess | Reconstruction of finite volume solution for parameter-dependent linear hyperbolic conservation laws[END_REF] is only a tool to compute efficiently a robust approximation of the finite volume solution of linear conservation laws, for any parameter values, from precomputed snapshots. It is a first step toward efficient nonlinear ROM strategies, in particular to design dynamical RB method with adapted local basis in finite volume framework e.g. from some precomputed dictionnary. In the lines of [START_REF] Cagniart | Model Order Reduction for Hyperbolic Problems: a new framework[END_REF][START_REF] Rim | Manifold Approximations via Transported Subspaces: Model reduction for transport-dominated problems[END_REF], using RTA approach together with efficient strategies to compute both transformations and reduced basis would be interesting for attacking general hyperbolic conservation laws.

• More generally, designing approximation methods for solving parameter-dependent evolution equations has become an active field of research in the MOR community leading to the emergence of nonlinear approximation methods. Beyond approximation relying on transformed snapshots, many effort has also been put on finding suitable approximation format for attacking the approximation of the solution itself but the question of constructing such an approximation, when only the equation is available, remains a challenging task. In this direction, recent approaches involving machine learning techniques such as e.g. Neural Galerkin schemes [START_REF] Bruna | Eijnden Neural Galerkin Scheme with Active Learning for High-Dimensional Evolution Equations[END_REF], that are dynamical nonlinear approximation methods, seem appealing.

Randomized linear algebra for dynamical reduced basis. Despite the proposed time-dependent RB approach yields a significant reduction of online computation costs, it requires during the offline step the storage of the whole trajectory of the snapshots as well as related reduced quantities (matrices, error estimates etc.) which can be limiting for very large complex numerical model. Recently, novel projection methods based on randomized linear algebra have been proposed in [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF][START_REF] Balabanov | Randomized linear algebra for model reduction. Part II: minimal residual methods and dictionary-based approximation[END_REF]. The idea behind is to approximate the ROM from its random sketch, which is a set of low-dimensional random projections of the reduced approximation space and the spaces of associated residuals. An interesting fact is that such an approach can be used for reducing both complexity and memory requirements of classical ROM approaches. Moreover, it allows to provide sharp error estimates that can be computed efficiently [START_REF] Homescu | Error Estimation for Reduced-Order Models of Dynamical Systems[END_REF][START_REF] Smetana | Randomized Residual-Based Error Estimators for Parametrized Equations[END_REF][START_REF] Smetana | Randomized residual-based error estimators for the proper generalized decomposition approximation of parametrized problems[END_REF] and improve projection stability [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF]. It would be of interest of extending such approaches within the framework of dynamical RB proposed in this section for parameter and time dependent problem.

Geometry based dynamical low-rank approximation

We now focus on low-rank approximation of time-dependent matrices A(t) ∈ R n×m . Introducing Ȧ(t) = d dt A(t) the time derivative, the matrix A(t) is defined as the solution of the following matrix ODE

Ȧ(t) = F (A(t), t), A(0) = A 0 , (2.3.1)
given A 0 ∈ R n×m and F : R n×m × [0, T ] → R n×m . In view of MOR, the goal of DLRA methods methods is to approximate the solution A(t) of (2.3.1) with a matrix Z(t), cheaper to compute, which belongs to the nonlinear manifold of fixed rank matrices

M r (R n×m ) = {Z ∈ R n×m : rank(Z) = r},
where r min(n, m) stands for the rank. When A(t) is known, Z(t) can be defined as the best rank-r approximation solution of

Z(t) = arg min W ∈Mr(R n×m ) A(t) -W , (2.3.2)
with • the Frobenius norm. In that case, Z is obtained through a truncated SVD of A(t) for each instant t. Nevertheless, as A is implicitly given by (2.3.1), it is more relevant to introduce low-rank approximation using Ȧ. To that goal, the approximation Z is classically obtained through its derivative Ż which satisfies the Dirac-Frenkel variational principle

Ż(t) = arg min δW ∈T Z(t) Mr(R n×m ) δW -F (Z(t), t) , (2.3.3)
given Z(0) = Z 0 ∈ M r (R n×m ) the best rank-r approximation of A(0) and T Z(t) M r (R n×m ) the tangent space to M r (R n×m ) at Z(t). Equivalently, Ż(t) corresponds to the orthogonal projection of F (Z(t), t) on the solution dependent tangent space, i.e.

Ż(t) = P T Z (t) F (Z(t), t), Z(0) = Z 0 , (2.3.4)
where P T Z denotes the projection onto T Z(t) M r (R n×m ).

Towards numerics. One way to deal with the numerical integration of (2.3.4) is to use projection based methods (see, e.g., [START_REF] Kieri | Projection methods for dynamical low-rank approximation of high-dimensional problems[END_REF]), in the lines of Riemaniann optimization. Such methods work as follows. Perform one step of the numerical scheme (e.g. Runge Kutta) leaving M r (R n×m ), and then project back to it by means of retraction. This latter step is usually performed using a r-terms truncated SVD. As these approaches work directly in the ambiant space R n×m , they do not exploit the geometry of the set of fixed rank matrices M r (R n×m ). Following [START_REF] Koch | Dynamical low-rank approximation[END_REF], we adopt a second point of view and propose numerical integration of the differential equation (2.3.4) using suitable geometric description of M r (R n×m ) which requires to deals with the subsequent difficulties.

1. The first difficulty relies on the proper geometric description of the set of fixed rank matrices M r (R n×m ). In practice, a way to compute the rank-r matrix Z(t) is done through its parametrization

Z(t) = U (t)G(t)V (t) T , (2.3.5)
with U (t) ∈ R n×r , V (t) ∈ R m×r and G(t) ∈ R r×r . Such a parametrization of the matrix Z(t) is not unique. A way to dodge this undesirable property is to properly define the tangent space T Z(t) M r (R n×m ) imposing the so-called gauge conditions on U, V . In that case, the matrix Z(t) admits a unique decomposition (see [START_REF] Koch | Dynamical low-rank approximation[END_REF]Proposition 2.1]) of the form (2.3.5). In addition, the system (2.3.4) results in a system ODEs driving the evolution of the parameters U, G, V .

2. The second difficulty appears when numerical integration is performed for solving the resulting system of ODEs governing the evolution of parameters U, V and G [START_REF] Koch | Dynamical low-rank approximation[END_REF][START_REF] Nonnenmacher | Dynamical low-rank approximation: applications and numerical experiments[END_REF].

In presence of small singular values for Z(t), the matrix G(t) may be ill-conditioned. As consequence, classical integration schemes may be unstable for large rank and time steps (see, e.g., [START_REF] Kieri | Discretized dynamical low-rank approximation in the presence of small singular values[END_REF]Section 2.1]). Moreover, in case of overapproximation, i.e. when the approximation Z(t) has a rank r greater than the rank of the exact solution A(t), these methods fail since G(t) becomes singular.

Using a chart based geometric description of M r (R n×m ) relying on a natural parametrization of matrices, we have proposed numerical integrator in local coordinates for solving (2.3.4), suitable for MOR.

Geometry of the set of fixed-rank matrices

In what follows, we present briefly the geometric description of the fixed-rank matrix manifold and the associated tangent space proposed in [START_REF] Billaud-Friess | Principal bundle structure of matrix manifolds[END_REF].

Any matrix Z ∈ M r (R n×m ) can be represented (in non-unique way) by means of the factorization

Z = U GV T , (2.3.6)
with U ∈ M r (R n×r ), V ∈ M r (R m×r ) and G ∈ GL r with GL r the Lie group of r × r invertible matrices. We consider

U ⊥ ∈ M n-r (R n×(n-r) ), V ⊥ ∈ M m-r (R n×(m-r) ), the matrices such that U T ⊥ U = 0 and V T ⊥ V = 0. The neighborhood U Z of Z in M r (R n×m
) is defined as the set

U Z = {(U + U ⊥ X)H(V + V ⊥ Y ) T : (X, Y, H) ∈ R (n-r)×r × R (m-r)×r × GL r }. M r (R n×m ) Manifold U Z • × W = (U + U ⊥ X)H(V + V ⊥ Y ) T Z = U GV T Parameters (X, Y ) H • (0, 0) G× (R (n-r)×r × R (m-r)×r ) × GL r θ Z Figure 2.3.1: Representation of the local chart θ Z that associates to W = (U + U ⊥ X)H(V + V ⊥ Y ) T in U Z ⊂ M r (R n×m ) the parameters (X, Y, H) ∈ R (n-r)×r × R (m-r)×r × GL r .
We associate to the neighborhood U Z of Z the local chart θ Z : U Z → R (n-r)×r × R (m-r)×r × GL r (see Figure 2.3.1) which is given by

θ Z (W ) = (U + ⊥ W (V + ) T (U + W (V + ) T ) -1 , V + ⊥ W T (U + ) T (V + W T (U + ) T ) -1 , U + W (V + ) T )
for any W ∈ U Z . Here U + and V + stand for the Moore-Penrose pseudo-inverses1 of U and V respectively. This means that any matrix W belonging to the neighborhood U Z admits a unique parametrization

W = (U + U ⊥ X)H(V + V ⊥ Y ) T , with parameters (X, Y, H) ∈ R (n-r)×r × R (m-r)×r × GL r .
In this description, the parameters are not longer U, V, G but X, Y, H.

Such geometric description confers the set M r (R n×m ) the structure of an analytic r(n + mr)dimensional manifold. Moreover, M r (R n×m ) is an analytic principal bundle with typical fiber GL r and base

2 G r (R n ) × G r (R m ) (see [MBF13, Theorem 7]).
We now define the tangent map at Z ∈ M r (R n×m ) noted T Z i by

T Z i : R (n-r)×r × R (m-r)×r × R r×r → R n×m , (δX, δY, δH) → U ⊥ δXGV T + U G(V ⊥ δY ) T + U δHV T .
Then, the tangent space to M r (R n×m ) at Z is defined as the image through T Z i of the tangent space in the local coordinates in R (n-r)×r × R (m-r)×r × R r×r

T Z M r (R n×m ) = {U ⊥ δXGV T + U G(V ⊥ δY ) T + U δGV T : δX ∈ R (n-r)×r , δY ∈ R (m-r)×r , δG ∈ R r×r }.
As stated in [MBF13, Proposition 6], T Z i is an isomorphism between T Z M r (R n×m ) and R (n-r)×r × R (m-r)×r ×R r×r . Thus, any tangent matrix δZ ∈ T Z M r (R n×m ) admits a unique parametrization of the form

δZ = T Z i(δX, δY, δH) = U ⊥ δXGV T + U G(V ⊥ δY ) T + U δHV T , (2.3.7)
where (δX, δY, δH) ∈ R (n-r)×r × R (m-r)×r × R r×r are uniquely given through

δX = U + ⊥ δZ(V + ) T G -1 , δY = V + ⊥ δZ T (U + ) T G -T , δH = U + δZ(V + ) T .
(2.3.8)

Dynamical low-rank approximation

In the context of dynamical low-rank approximation, we recall that Z is given through the projected Equation (2.3.4). By definition of the tangent map (2.3.7), the tangent matrix is given by Ż = T Z i( Ẋ, Ẏ , Ḣ) where the parameters satisfies

( Ẋ, Ẏ , Ḣ) = T Z i -1 (P T Z F (Z))
(2.3.9) with

P T Z F (Z) = P ⊥ U F (Z)P T V + P T U F (Z)(P ⊥ V ) T + P U F (Z)P T V .
Here P U = U U + , P V = V V + denote the projections associated to U, V respectively, and their related orthogonal projections P ⊥ U = I -P U , P ⊥ V = I -P V . Equation (2.3.9) yields equivalently to the following system of ODEs on the parameters:

Ẋ = U + ⊥ F (Z)(V + ) T G -1 , Ẏ = V + ⊥ F (Z) T (U + ) T G -T , Ḣ = U + F (Z)(V + ) T .
(2.3.10) Unlike [114, Proposition 2.1], the proposed geometrical description allows to avoid the non uniqueness of the parametrization (2.3.5) without imposing some gauge conditions a posteriori trough the tangent space. More precisely, if we impose ( U , V , Ġ) = (U ⊥ Ẋ, V ⊥ Ż, Ḣ) and assuming that U, V orthogonal, our parametrization naturally fulfills the gauge conditions U T U = 0 and V T V = 0 and we recover the following system of ODEs [START_REF] Koch | Dynamical low-rank approximation[END_REF] governing the evolution of the factors U, G, V

U = (I -U U T )F (Z)V G -1 , V = (I -V V T )F (Z) T U G -T , Ġ = U T F (Z)V.
(2.3.11)

Remark 2.3.1. Let us remark that if F : R n×m → R n×m is a C p is vector field then (2.3.1) admits a unique maximal solution. With the proposed chart-based based description, the projected dynamical system (2.3.9) inherits such properties. Indeed, T Z i -1 P T Z F : R n×m → R (n-r)×r × R (m-r)×r × R r×r is also a C p vector field [78, Theorem 2.6, Theorem 3.5].

Integration schemes

In seminal work [START_REF] Lubich | A projection-splitting integrator for dynamical low-rank approximation[END_REF], and more recently in [START_REF] Ceruti | Time integration of symmetric and anti-symmetric low-rank matrices and Tucker tensors[END_REF][START_REF] Ceruti | An unconventional robust integrator for dynamical low-rank approximation[END_REF], projection-splitting integrators have been proposed to deal with numerical integration of (2.3.4). They have the advantage to be accurate, robust and stable independently of the rank r and do not need any regularization in case of over-approximation. The integrator initially proposed in [START_REF] Lubich | A projection-splitting integrator for dynamical low-rank approximation[END_REF], is based on an explicit scheme using a Lie-Trotter splitting of the projection operator P T Z(t) . In addition to its simplicity, it has the advantage to remains robust in case of small singular values and especially for overapproximation as it avoids the inversion of G(t). Following the same idea, we have derived in [START_REF] Billaud-Friess | A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds[END_REF] a first order numerical integrator in local coordinates for solving (2.3.4) that can be interpreted as a splitting integrator.

Symmetric splitting method

We first consider the setting of the classical description [START_REF] Koch | Dynamical low-rank approximation[END_REF]. To perform time integration, a symmetric Lie-Trotter splitting method [START_REF] Lubich | A projection-splitting integrator for dynamical low-rank approximation[END_REF] is applied to Equation (2.3.3). This integration scheme relies on a decomposition of the projection P T Z (t) as follows

P T Z (t) = Q 1 -Q 2 + Q 3 , (2.3.12)
where Q 1 , Q 2 and Q 3 are three projections respectively defined by

Q 1 A = AP T V , Q 2 A = P U AP T V , Q 3 A = P U A, (2.3.13)
for any A ∈ R n×m . Using this splitting, one integration step from t 0 to t 1 starting from the factorized rank-r matrix Z 0 = Z(t 0 ) under the form Z 0 = U 0 G 0 V T 0 reads as follows.

1. Integrate on [t 0 , t 1 ] the n × r matrix differential equation

d dt (U G) = F ((U G)V T )V, V = 0, with initial conditions (U G)(t 0 ) = U 0 G 0 , V (t 0 ) = V 0 . Then set U 1 and Ĝ1 such that U 1 Ĝ1 = (U G)(t 1 ).
2. Integrate on [t 0 , t 1 ] the r × r matrix differential equation

Ġ = -U T F (U GV T )V, U = 0, V = 0, with initial conditions G(t 0 ) = Ĝ1 , U (t 0 ) = U 1 , V (t 0 ) = V 0 . Set G1 = G(t 1 ). 3. Integrate on [t 0 , t 1 ] the m × r matrix differential equation d dt (V G T ) = F (U (V G T ) T ) T U, U = 0, with initial conditions (V G) T (t 0 ) = V 0 GT 1 , U (t 0 ) = U 1 . Set V 1 and G 1 such that V 1 G T 1 = (V G T )(t 1 ).
After these steps, the obtained approximation is

Z(t 1 ) = U 1 G 1 V T 1 .
Each step corresponds to the integration of the right hand side of (2.3.4) associated to the projections Q 1 , Q 2 and Q 3 respectively (for details see [122, Lemma 3.1]).

Chart based splitting method

In [START_REF] Billaud-Friess | A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds[END_REF], we have proposed a numerical integrator relying on the fibre bundle structure of the manifold of fixed rank matrices. Thus, the guiding idea is to perform some update of the parameters (X, Y, H), instead of (U, G, V ) directly, whose dynamic is governed by the system of ODEs (2.3.10). Working in a fixed neighborhood U Z of Z, the matrices U, V, G are fixed and Equation (2.3.10) writes equivalently

Ḣ = U + F (Z)(V + ) T , ẊG = U + ⊥ F (Z)(V + ) T , Ẏ G T = V + ⊥ F (Z) T (U + ) T .
(2.3.14)

Then, we integrate the system (2.3.14) from t 0 to t 1 in three steps. Letting

U (t) = U (t 0 ) + U (t 0 ) ⊥ X(t), V (t) = V (t 0 )+V (t 0 ) ⊥ Y (t) and Z(t) = U (t)H(t)V (t) T , we start from (X(t 0 ), Y (t 0 ), H(t 0 )) = (0, 0, G(t 0
)) and we proceed as follows.

1. Integrate on [t 0 , t 1 ] the r × r matrix differential equation

Ḣ = U + F (Z)(V + ) T , Ẋ = 0, Ẏ = 0, with initial conditions X(t 0 ) = 0, H(t 0 ) = G(t 0 ) and Y (t 0 ) = 0. Set H 1 = H(t 1 ).
2. Integrate on [t 0 , t 1 ] the n × r matrix differential equation

ẊH = U + ⊥ F (Z)(V + ) T , Ḣ = 0, Ẏ = 0,
with initial conditions X(t 0 ) = 0, H(t 0 ) = H 1 and Y (t 0 ) = 0. Then set X 1 = X(t 1 ).

3.

Integrate on [t 0 , t 1 ] the m × r matrix differential equation

Ẏ H T = V + ⊥ F (Z) T (U + ) T , Ḣ = 0, Ẋ = 0. with initial conditions X(t 0 ) = X 1 , H(t 0 ) = H 1 and Y (t 0 ) = 0. Then set Y 1 = Y (t 1 ).
After these three steps, we obtain an approximation Z(t

1 ) := U 1 H 1 V T 1 with U 1 = U (t 0 ) + U (t 0 ) ⊥ X 1 and V 1 = V (t 0 ) + V (t 0 ) ⊥ Y 1 .
The chart based method can be interpreted as a Lie-Trotter splitting that relies on the following decomposition of the projection

P T Z (t) = P 1 + P 2 + P 3 , (2.3.15)
where P 1 , P 2 and P 3 are three projections respectively defined by

P 1 A = P U AP T V , P 2 A = P ⊥ U AP T V , P 3 A = P U A(P ⊥ V ) T , ( 2

.3.16)

for any A ∈ R n×m . Here each term P i of the projection is associated to the ODE solved at Step i (see [MBF14, Appendix 1]).

Comparaison of the two approaches

• The inversion of G in the symmetric splitting method (and H in the chart based splitting method) is avoided. This convenient choice allows to deal with the case of overapproximation.

• The symmetric splitting and the chart based methods differ by the update order. Indeed, the chart based method first updates H and then X, Y (or equivalently G and then U, V ) whereas the symmetric splitting updates U G, then G and finally V G. Moreover, Step 2 of the symmetric splitting method described in Section 2.3.3 can be interpreted as a backward evolution problem that can be ill-conditionned, as pointed out in [9, Section 5].

In the chart based method, the update of G at Step 1 is still a forward evolution problem due to our splitting choice.

• In practice, we provide practical formulation of those methods amenable for numerical use. Explicit approximation of the flux is performed leading first order accurate schemes.

The two method are proven to coincide and are exact for matrix approximation [MBF14, Lemma 3.11].

Numerical applications

We consider the approximation of the parameter-dependent Burger's viscous equation in one dimension. To that goal, let Ω × I = (0, 1) × [0, 1] be a space time domain. We seek u(ξ) the solution of The problem (2.3.17) is semi-discretized in space by means of finite difference (FD) schemes with n nodes and m instances of the parameter ξ such that we get the following dynamical system

∂ t u(ξ) -µ(ξ)∂ 2 xx u(ξ) + u(ξ)∂ x u(ξ) = f (ξ), on Ω × I, (2.3
Ẋ(t) = LX(t) + h(X, t), X(0) = 0, (2.3.18)
with solution X(t) ∈ R n×m . We denote L the operator obtained by tensorization of discrete Laplacian and diagonal matrix obtained from m instances of ξ 1 . Moreover, we define the matrix valued function h obtained by the discrete version of the first derivative.

In what follows we denote KSL the symmetric approach and chart based approach our approach.

Single parameter Burger's problem

In this section, the parameter ξ 1 takes its values in [0.01, 0.06] while the others are fixed to ξ 2 = 4 and ξ 3 = 0. We chose the following initial condition u 0 (x; ξ) = sin(xξ 1 )e -100(x-10ξ 1 ) 2 .

and set n = 100 and m = 60. For solving the matrix ODE given by Equation (2.3.18), we first confront the chart based and KSL algorithms for various ranks, fixing ∆t = 10 -4 . The approximations obtained are compared to a reference solution noted {X k ref } K k=0 computed with an explicit Euler scheme with ∆t = 5.10 -6 .

On Figure 2.3.2, the evolution of e k = X k ref -Z k for both algorithms is studied for different ranks. As we can observe, both methods seem to provide an approximation with similar accuracy, except for ranks 15 and 20 where the chart based algorithm provides a slightly more accurate approximation than the KSL algorithm. This observation is confirmed by Figure 2.3.3 where we analyze the convergence with respect to rank and time step. As we can see, the error e K = X K ref -Z K decreases with respect to the rank and time step, and is the smallest (by up to one order of magnitude) for the chart based algorithm for larger rank.

Multiple varying parameters Burger's problem

To conclude this section, we illustrate the behavior of the two methods for the case where ξ is a vector of independent random parameters uniformly distributed on [0.01, 0.06]× [START_REF] Abgrall | Crisovan Model reduction using L 1 -norm minimization as an application to nonlinear hyperbolic problems[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF]×[0.01, 0.1] with u 0 (x; ξ) = ξ 3 e -100(x-10ξ 2 ) 2 . The numerical simulations are performed for m = 60 and n = 100. Here, the two KSL and Chart methods are run with ∆t = 10 -4 and compared to true numerical reference solution {X k ref } K k=0 obtained with the explicit Euler scheme with the same time step. We represent on Figure 2.3.4 the approximation error at final time for the three approaches with respect to r ∈ {5, 10, 12, 15, 20, 25, 30}. We clearly observe that the chart method provides a better approximation than KSL for which the error seems to stagnate after r ≥ 12.

Summary and discussion for future research

We have introduced and compared some geometry based algorithms for dynamical low-rank approximation. Using a different geometry description of the set of fixed rank matrices relying on charts, we generalized the description of [START_REF] Koch | Dynamical low-rank approximation[END_REF]. Then, from this description we derive a new splitting algorithm motivated by fibre bundle structure of the manifold of fixed rank matrices.

Let us now outline some possible perspectives and applications.

Analysis of convergence order and stability of the splitting scheme. The presented results demonstrate the feasibility of an alternative splitting naturally provided by the proposed geometric description. It would clearly be interesting to conduct further numerical experiments and theoretical analysis, to better understand the differences with KSL splitting algorithm, especially in term of accuracy with respect to time step and rank following [START_REF] Kieri | Discretized dynamical low-rank approximation in the presence of small singular values[END_REF] as well as a study of the influence of backward/forward integration for G on the scheme stability in the lines of [START_REF] Bachmayr | Existence of dynamical lowrank approximations to parabolic problems[END_REF][START_REF] Kazashi | Stability properties of a projector-splitting scheme for dynamical low rank approximation of random parabolic equations[END_REF].

Alternating scheme relying on the fibre bundle structure of the set of fixed rank matrices.

The proposed splitting scheme is a first step towards designing new algorithms integrating the geometric structure of the manifold of fixed rank matrices, by working in neighborhoods. Especially, deriving a numerical scheme working alternatively with the fiber and base, (with possible parallelization as in [START_REF] Ceruti | An unconventional robust integrator for dynamical low-rank approximation[END_REF] to enhance computation efficiency) should be interesting.

Toward adaptive DLRA. Presented methods compute low-rank approximations with a fixed rank r. This choice, done a priori, may not be optimal at each time t. Especially if the actual rank of X(t) varies with time (e.g. for applications arising from the discretization of a particular PDE). This is a first reason why allowing the rank of the approximation to evolve with time is interesting. Moreover, adapting the rank through the iterative process enables some computational time saving. DLRA with rank adaptivity has been the object of recent prospects. In [START_REF] Ceruti | A rank-adaptive robust integrator for dynamical lowrank approximation[END_REF] the authors propose a variant of the splitting integrator introduced in [START_REF] Ceruti | An unconventional robust integrator for dynamical low-rank approximation[END_REF] with an additional truncation step, to allow-rank adaptivity.

• To our knowledge, rank adaptivity is introduced at numerical level trough efficient algorithms in view of preserving some particular features as : exactness, robustness or symmetry properties. An interesting question concerns the possibility of formulating an adaptive DRLA method based on a global optimization in a "dynamical" low-rank manifold with suitable geometry description (e.g. integrating the rank).

• Finally, when considering the general case of tensors, as in [START_REF] Ceruti | Rank-adaptive time integration of tree tensor networks[END_REF][START_REF] Dektor | Rank-adaptive tensor methods for high-dimensional nonlinear PDEs[END_REF][START_REF] Ehrlacher | Adaptive Hierarchical Subtensor Partitioning for Tensor Compression[END_REF][START_REF] Rodgers | Adaptive integration of nonlinear evolution equations on tensor manifolds[END_REF], one could imagine the development of adaptive algorithms for the construction of DLRA that adapts dynamically not only the rank but also the tensor format, e.g., with dynamical choice of tree-based formats.

Revisiting classical Riemannian optimization algorithm. To finish, let us mention that the geometric description of the low-rank matrix manifold, proposed in [START_REF] Billaud-Friess | Principal bundle structure of matrix manifolds[END_REF], was motivated within the context of DLRA for parameter-dependent dynamical systems. However, it can be carried over to optimization in the manifold of fixed rank matrices, in particular to revisit classical Riemannian optimization algorithm. This has been object of investigation during the internship of M. M. Abedin Nejad [START_REF] Abedin Nejad | Optimization Algorithms on Tensor Manifolds using Charts Intership report[END_REF].

Chapter 3

Approximation for PDEs with probabilistic interpretation

In the previous chapters, we focused on low-rank approximation methods for the MOR of complex numerical models typically arising from the discretization of parameter-dependent PDEs.

Here, we consider another direction by addressing directly the problem of constructing an approximation of the solution of the considered PDEs (possibly parameter-dependent) using approaches that rely on pointwise evaluations (possibly noisy) of the function to approximate.

This chapter is organized as follows. After giving a brief discussion on probabilistic approximation methods in Section 3.1, we give required notations and results concerning Monte Carlo method for estimating the solution of elliptic PDEs in Section 3.2. In section 3.3, we introduce a probabilistic sparse polynomial interpolation methods to deal with high dimensional equations. Then, in section 3.4, we discuss some probabilistic RB method for parameter dependent equations.

General context

In this chapter, we focus on computing an approximation of the solution u : D → R of the following continuous boundary problem (possibly parameter dependent)

-Au = g in D, u = f on ∂D, (3.1.1)
where D is an open bounded domain in R d with boundary ∂D. Here, A is some linear elliptic partial differential operator and the functions f : ∂D → R and g : D → R are the source and boundary terms respectively. Here, we assume that u lies in V , some high dimensional vector space.

Probabilistic approaches for approximation

When solving a problem of the form (3.1.1), standard numerical approaches (e.g. finite difference, finite element ...) are usually considered to provide a global approximation u n of u, for a given mesh of D, in a finite dimensional space V n ⊂ V . Such methods meet their limit when addressing parameter dependent problems, when many evaluations of the complex numerical solution u n is required, or are no longer suitable in high-dimension. As discussed in previous chapters, approximations in low-rank tensor formats are good candidates to tackle those limitations. Classically, the presented methods are deterministic and rely on projection principles or optimization problems, that provide an approximation u r of the complex numerical solution u n in linear subspace V r ⊂ V n of small dimension. Here, we follow another path, by considering probabilistic approaches, in the sense that they integrate some randomness, which only use pointwise evaluations of u.

Monte Carlo methods for approximation

Monte Carlo (MC) methods are well known to be efficient probabilistic approaches when computing numerically an integral or an expectation (see, e.g., monographs [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF][START_REF] Graham | Stochastic Simulation and Monte Carlo Methods[END_REF]), especially in high dimension.

When turning to the computation of the solution of a PDE under the form (3.1.1), MC method allows to provide an estimate of the solution u(x) at any point x ∈ D without suffering from the curse of dimensionality. Such an estimate relies on the so-called Feynman-Kac probabilistic representation formula. In brief, for any x ∈ D, u(x) can be interpreted as the expectation of some given stochastic functional depending on the operator A and the data f, g. The main limitations of such approaches are that they only provide pointwise estimates of u and may suffer from slow convergence of MC estimate. To tackle the first limitation, it is possible to compute approximations in a given approximation format through classical interpolation or regression (see, e.g., [START_REF] Beck | Machine learning approximation algorithms for highdimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations[END_REF][START_REF] Beck | Solving the Kolmogorov PDE by means of deep learning[END_REF][START_REF] Weinan | Deep learning-based numerical methods for highdimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]) within deep learning context. Let us mention, that it is also possible to consider some general approximation setting relying on the minimization of an empirical risk on a selected model class [START_REF] Eigel | Variational Monte Carlo-bridging concepts of machine learning and high-dimensional partial differential equations[END_REF][START_REF] Eigel | Convergence bounds for empirical nonlinear least-squares[END_REF]. The resulting method called variational Monte Carlo method is shown to be applicable to a broad range of problems including e.g. approximation using neural networks, or approximation in low-rank tensor format. Concerning the second difficulty, variance reduction techniques (see Section 3.2) can be considered. In particular let us mention that 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 have proposed to combine some polynomial interpolation for reconstructing a global approximation u n of u in some polynomial space V n ⊂ V from MC estimates, such that the approximation error e n = uu n can be used as control variate during some iterative procedure.

Probabilistic approaches within MOR

In MOR, various probabilistic approaches have recently emerged in view of improving the efficiency of the underlying methods either in term of their approximation properties or their numerical complexity and stability.

Within the framework of RB method for parameter-dependent problems, the authors in [START_REF] Cohen | Reduced Basis Greedy Selection Using Random Training Sets[END_REF] have proposed a probabilistic greedy algorithm that uses, at each step, different training sets with moderate cardinality, randomly chosen. In practice, the resulting method allows to deal efficiently with high dimensional parametric problems. In [START_REF] Cai | A Stochastic Discrete Empirical Interpolation Approach for Parameterized Systems[END_REF], the authors derive a similar probabilistic EIM using sequential sampling in the parameter set, which provides an interpolation with a prescribed precision with high probability. Probabilistic approaches can be also used to provide numerically stable and efficient error estimates for low-rank approximations, see, e.g., [START_REF] Smetana | Randomized Residual-Based Error Estimators for Parametrized Equations[END_REF][START_REF] Smetana | Randomized residual-based error estimators for the proper generalized decomposition approximation of parametrized problems[END_REF] for contributions using randomized (a posteriori) error estimates. In [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF][START_REF] Balabanov | Randomized linear algebra for model reduction. Part II: minimal residual methods and dictionary-based approximation[END_REF], this idea has been generalized within the framework of randomized linear algebra to enhance the efficiency and stability of projection-based model order reduction methods for solving parameter-dependent equations. In particular, it allows drastic computational savings. Finally, let us also mention that a control variate using a reduced basis paradigm, has been proposed in [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF][START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF] for MC estimate of the expectation of some parameter dependent random variable.

Contributions

Within this chapter, our main concern is the development of probabilistic approximation methods, as previously discussed, that use on pointwise estimates, over D, of the solution u of (3.1.1).

To that goal, probabilistic interpretation of PDEs using the Feynman-Kac formula was the com-mon thread in contributions [MBF12, MBF15, MBF17] within the PhD of A. Macherey [START_REF] Macherey | Approximation and model reduction for partial differential equa-tions with probabilistic interpretation[END_REF] supervised together with A. Nouy (LMJL, ECN) and C. Prieur (LJK, UGA). Especially, we have followed two subsequent directions.

• From pointwise estimates of the solution u over D, we propose to compute a global approximation through polynomial interpolation overcoming possible limitations due to high dimension and slow convergence rates of Monte Carlo estimates.

• The idea is to propose a RB method in full sampling setting, i.e. that uses pointwise estimates of the solution u over D for given instances of the parameter. In particular, we propose a greedy algorithm using some probabilistic MC error estimate. In that case, the difficulty is to perform such algorithm keeping a reasonable budget (for MC sampling) in particular during the selection of the parameter and ensuring that the resulting procedure is a weak greedy algorithm with high probability.

Monte Carlo methods for PDEs

In this section, we introduce essential ingredients and notations concerning numerical probabilistic methods for PDEs that will be considered in the next sections. A deeper discussion on that topic can be found e.g. in [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF].

Notations. In what follows, (Ω, F, P) denotes a probability space, with Ω is the sample space, F a σ-algebra and P a probability measure.

Monte Carlo methods

Let Z be a random variable with finite expectation E(Z) and variance V(Z). MC methods are simple methods that compute an estimate of E(Z) through the empirical mean

ZM := 1 M M m=1 Z m , (3.2.1)
where {Z 1 , . . . , Z M } are M independent random variables having the same distribution as Z.

It is well known, that ZM is an unbiased estimate that converges almost surely to E(Z), by the strong law of large numbers. At this point, a natural question is to understand how good is the estimate ZM for the expectation E(Z), in which sense (e.g. error bounds via confidence interval, root mean square...) and as M → ∞.

Error bounds via confidence intervals

It is possible to derive some asymptotic error bounds, through confidence intervals, of the error |E(Z) -ZM |, as M → ∞. Using Central Limit theorem (other asymptotic controls are also possible see [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF]Chapter 2.3]), we get the following result. Let x ∈ (0, 1) and c(M, x) = z x VM M with, we have

P |Z M -E[Z]| ≤ c(M, x) → 1 -x,
as M → ∞. Here z x ∈ R satisfies F (z x ) = 1x/2 with F the cumulative distribution function of the centered reduced normal distribution and

V M = 1 M M m=1 Z m -Z M 2
is the empirical variance. This result means that for M large enough, the error is bounded by c(M, x) with a probability 1x. In practice, the question of how to choose M , to use these error bounds, remains a delicate task. Alternatively, non-asymptotic error bounds, such that

P |Z M -E[Z]| ≤ c(M, x) ≥ 1 -x, (3.2.2)
is valid for any value of M , are preferable. Some of this non-asymptotic error bounds are related to concentration inequalities (e.g. Heoffding, Bernstein). We refer the interested reader to [START_REF] Boucheron | Concentration Inequalities: A Non asymptotic Theory of Independence[END_REF] for a detailed overview.

Root mean square error and variance reduction techniques

Another relevant measure of the quality of the estimate (3.2.1) is given by

E(|E(Z) -ZM | 2 ) = V( ZM ) = V(Z) M . (3.2.3)
This result says that ZM is an estimate of the expectation E(Z) for which the root mean square error behaves like O(M -1/2 ). This show the limitations of the MC method that requires a large number of samples of Z to converge. In practice, this may induce high computation cost when the samples of Z are costly to compute (e.g. when Z is related to the evaluation of the solution of a PDEs).

The so-called variance reduction methods [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF]Chapter 3], are aimed at improving the variance in (3.2.3). Among these approaches, let mention the antithetic sampling, conditioning or stratification and importance sampling methods. Some variance reduction methods use control variate. They consist in using an additional centered square integrable random variable Y that is simulated jointly with Z such that

Zc M := 1 M M m=1 (Z m -Y m ), (3.2.4)
with {Y 1 , . . . , Y M } is a M -sample of Y , is a new unbiased estimate of E(Z), that converges to E(Z). In practice, Y is chosen such that the new estimate variance

V( Zc M ) = E(|E(Z) -Zc M | 2 ) = V(Z -Y ) M ,
is smaller than V(Z). A possible optimal choice, is to seek Y as a linear combinaison of finite number of random variables, that minimize E((Z -Y ) 2 ). This approach has been extended to multilevel approaches that construct sequences of control variates via Multilevel Monte-Carlo method (MLMC) (see survey [START_REF] Giles | Multilevel Monte Carlo methods Acta Numer[END_REF]). Let us mention that other reduction variance approaches, using control variates, have been also developed in the context of the numerical simulation of PDEs, possibly parameter-dependent (see, e.g., [START_REF] Blel | Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables[END_REF][START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF][START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF]).

Probabilistic representation of solution of PDEs

Notations. Let W = (W t ) t≥0 be a d-dimensional brownian motion, i.e. W = (W 1 , . . . , W d ) T where each W i , i = 1, . . . , d are independent standard brownian motions defined on (Ω, F, P) endowed with its natural filtration (F t ) t≥0 .

Given D be an open bounded domain in R d , we consider the following boundary problem

-Au = g in D, u = f on ∂D, (3.2.5)
with solution u : D → R and f : ∂D → R, g : D → R the boundary condition and source terms, respectively. Here, we assume that A is linear and elliptic partial differential operator given by

A = 1 2 d i,j=1 (σσ T ) ij ∂ 2 x i x j + d i=1 b i ∂ x i . (3.2.6)
It corresponds to the infinitesimal generator associated to the d-dimensional diffusion process X x = (X x t ) t≥0 which is the solution of the stochastic differential equation Denoting a = σσ T , we introduce the following uniform ellipticity assumption.

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

Assumption 2 (A2).

There exists c > 0 such that

y T a(x)y ≥ c y 2 , y ∈ R d , x ∈ D.
As problem (3.2.5) is defined on a bounded domain, we have to define the first exit time of D for the process X x as τ x = inf {s > 0 : X x s / ∈ D} .

(3.2.9) Also, we assume some regularity property on the spatial domain D and data.

Assumption 3 (A3). The domain D is an open connected bounded domain of R d , regular according to the two following assumptions:

i) P(τ x = 0) = 1, x ∈ ∂D,
ii) 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}.

Assumption 4 (A4). We assume that f is continuous on ∂D, g is Hölder-continuous on D.

The following probabilistic representation theorem [57, Theorem 2.4] holds.

(x) = E f (X x τ x ) + τ x 0 g(X x t )dt , (3.2.10)
with X x the stopped diffusion process solution of (3.2.7).

[87]. 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 , 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 ). Variance reduction methods, discussed in Section 3.2, can be used to further improve the convergence.

High dimensional PDEs

Keeping the same notations as in Section 3.2.2, we assume that u : D → R is solution of (3.2.5).

We recall that under the assumptions (A1)-(A4), u admits, by Theorem 3.2.1, the following probabilistic representation

u(x) = E[F (u, X x )] := E u(X x τ x ) - τ x 0 Au(X x t )dt , x ∈ D
where F depends on the data of the problem and X x is the stopped diffusion process solution of Equation (3.2.7). An MC estimate of u(x) is

u ∆t,M (x) := 1 M M m=1 F u, X x,∆t (ω m ) .
with X x,∆t the approximated stopped diffusion process obtained through Euler-Maruyama scheme.

Within this section, we consider the problem of reconstructing a global approximation of the solution u from pointwise estimates u ∆t,M (x) for any x ∈ D. It is structured as follows. In the Section 3.3.1, we discuss possible strategies for global reconstruction of a function from pointwise evaluations. Then, we present the control variate method for high dimensional PDEs proposed in [START_REF] Billaud-Friess | Stochastic Methods for Solving High-Dimensional Partial Differential Equations[END_REF]. Some numerical experiments illustrating the behavior of the proposed approach are given in Section 3.3.3.

Approximation from pointwise evaluations

We consider the problem of constructing a global approximation u n of u, in some finite dimensional linear space V n ⊂ V , from pointwise evaluations of u at given points Γ m = {x 1 , . . . , x m } ⊂ D with m ≥ n. Two projection based strategies are classically used to construct u n namely : the interpolation and least-squares methods.

Interpolation method

Well known interpolation methods include polynomial and spline interpolation as well as Kriging. Let u be continuous and Γ n be an unisolvent grid of interpolation points (here m = n) for V n , Lagrange interpolation methods aim at constructing the approximation u n of u as the unique function in V n such that

u n (x i ) = u(x i ), x i ∈ Γ n .
Such an approximation satisfies the following error bound in L ∞ -norm

u -u n ∞ ≤ (1 + L n ) inf v∈Vn u -v ∞ ,
where L n stands for the Lebesgue constant. The choice of the points Λ n is crucial to control the growth of the Lebesgue constant with respect to n. In dimension one, interpolation is quite well understood for polynomial approximation (see, e.g., [START_REF] Ibrahimoglu | Lebesgue functions and Lebesgue constants in polynomial interpolation[END_REF][START_REF] Smith | Lebesgue constants in polynomial interpolation[END_REF]). In that context, it is well known that taking uniformly distributed interpolation points leads to a Lebesgue constant L n that grows exponentialy with n. A suitable choice of interpolation points, is the sequence of Chebyshev points for which L n grows only logarithmically with n. When turning to general multivariate domains with arbitrary linear approximation space (i.e. V n is not necessary polynomial space), a general procedure consists in selecting the interpolation points by maximizing over D n of the determinant of the collocation matrix. This provides the sequence of Fekete points, for which the Lebesgue constant grows linearly with n. However, it remains computationally untractable in high dimension as it requires to solve a multivariate non convex optimization problem. Another procedure, first proposed in the context of reduced basis [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF], is to select interpolation points in greedy fashion providing sequence of nested interpolations points Γ n ⊂ Γ n+1 . This method has a reasonable computational cost for moderate d and is independent on n. The resulting interpolation points, called magic points [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF] coincide, in the univariate case for polynomial approximation, with Leja point. In that case, an upper bound of L n behaves quadratically with n. When turning to high dimensional problem, sparse polynomial interpolation methods [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF][START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs[END_REF] are suitable approaches that allow to construct adaptively the sequence of interpolation point {Γ n } n≥0 for which the Lebesgue constant behaves well and do not suffer from the curse of dimensionality.

Least-square method

Let u be a fonction in L 2 µ (D) with some probability measure µ on D. Least square fitting methods construct the approximation u n of u in V n as solution of the following least square regression problem

min v∈Vn m i=1 |u(x i ) -v(x i )| 2 (3.3.1)
where the {x 1 , . . . , x m } is m-sample for the probability measure µ. The least squares method is stable and accurate for large m large enough [START_REF] Cohen | On the Stability and Accuracy of Least Squares Approximations[END_REF]. Weighted least squares have been proposed to overcome this limitation. In that case, the approximation u n is sought as the solution of

min v∈Vn m i=1 w(x i )|u(x i ) -v(x i )| 2 (3.3.2)
where the w(x i ) are some weights related to {x 1 , . . . , x m }. Here the x i are independent and identically distributed copies of a new measure ρ. It is shown in [START_REF] Cohen | Optimal weighted least-squares methods[END_REF], that for an optimal choice of weight w and sampling measure ρ, stability and optimal accuracy are achieved if m scales linearly with n up to an additional logarithmic factor. To that goal the weight w and sampling measure ρ are given through the inverse Christoffel function that depends on the space V n and the measure µ. It can be boosted when the measure ρ is chosen so as to ensure the stability of the weighted least-squares projection [START_REF] Haberstich | Boosted optimal weighted least-squares[END_REF]. In that case, it has the advantage to requires a sample of size close to the ideal m = n, as for interpolation.

Control variate method for high dimensional PDEs

We seek a global approximation of u from pointwise estimates u ∆t,M (x) for any x ∈ D. In that context, we have proposed in [START_REF] Billaud-Friess | Stochastic Methods for Solving High-Dimensional Partial Differential Equations[END_REF] algorithms that combine the variance reduction method proposed in [START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF][START_REF] Gobet | Sequential control variates for functionals of Markov processes[END_REF] with adaptive sparse polynomial interpolation [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF][START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs[END_REF] for approximation of u for high-dimensional problem.

Notations. Let Λ ⊂ N d be a finite set of multi-indices and V Λ ⊂ C 2 (D) an approximation space of finite dimension #Λ and let I Λ : R D → V Λ be the interpolation operator associated with an unisolvent grid Γ Λ = {x ν : ν ∈ Λ}. We let {ϕ ν } ν∈Λ denote the (unique) basis of V Λ that satisfies the interpolation property ϕ ν (x µ ) = δ νµ for all ν, µ ∈ Λ. The interpolation I Λ (w) = ν∈Λ w(x ν )ϕ ν (x) of function w is then the unique function in V Λ such that

I Λ (w)(x ν ) = w(x ν ), ν ∈ Λ.
A sequential control variate method

The sequential control variate procedure of [START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF] provides a sequence of approximations {ũ k } k≥1 of u in V Λ , which are defined for all k ≥ 1 by ũk = ũk-1 + ẽk , where ẽk is an approximation of e k , solution of

-A(e k )(x) = g(x) + A(u k-1 )(x), x ∈ D, e k (x) = f (x) -u 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(u k-1 ) on D and fu k-1 on ∂D. Here, the approximation ẽk = I Λ (e k ∆t,M ) is defined as the interpolation 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 ). Now, let us introduce the time integration error at point x for a function h

e ∆t (h, x) = E[F (h, X ∆t,x )] -E[F (h, X x )].
Then the following theorem [MBF12, Theorem 2] gives a control of the error ũk+1 (x ν )u(x ν ) in expectation at step k of the control variate algorithm.

sup ν∈Λ E ũk+1 (x ν ) -u(x ν ) C(∆t, Λ) sup ν∈Λ E ũk (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 lim sup n→∞ sup ν∈Λ E ũk (x ν ) -u(x ν ) C 1 (∆t, Λ) 1 -C(∆t, Λ) . (3.3.3)
For ∆t sufficiently small (i.e. under the condition C(∆t, Λ) < 1), 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). This section concerns the approximation of the solution manifold M = {u(ξ) : ξ ∈ Ξ}. First, we motivate a probabilistic RB method in full sampling setting in Section 3.4.1. It leads us to consider some probabilistic greedy procedure with a bandit algorithm for the selection of a new snapshot detailed in Section 3.4.2.

E ũk -u ∞ C 1 (∆t, Λ) 1 -C(∆t, Λ) L Λ + u -I Λ (u) ∞ , (3.3
Here, we present some working directions and preliminary results obtained in [START_REF] Macherey | Approximation and model reduction for partial differential equa-tions with probabilistic interpretation[END_REF]. Practical implementation and numerical illustrations are under consideration for future publication [START_REF] Billaud-Friess | A probabilistic reduced basis method for parameter-dependent problems[END_REF].

Sample based probabilistic reduced basis method

Notations. Throughout this section, we assume that u(ξ) belongs to V ⊂ L 2 (D) a Hilbert space equipped with its inner product •, • V and the related norm • V .

Classical framework

We recall that RB method (see Section 1.3) is aimed at providing a low-rank approximation u r of u in a low-dimensional subspace V r ⊂ V , of small dimension r, following two steps: the offline and online stages.

During the offline step, the reduced basis space V r is deduced from snapshots of u computed for r instances of the parameter ξ selected in a training set Ξ ⊂ Ξ using a greedy procedure relying on some computable error estimate ∆(u r (ξ), ξ). It is summarized in Algorithm 3. (Step 2.) Compute u(ξ r ) and update V r = span{u(ξ 1 ), . . . , u(ξ r )}.

Usually, the Algorithm 3.4.1 is stopped when ∆(u r-1 (ξ), ξ) is below some target precision ε > 0 or for a given space dimension r. Then, during the offline stage, u r is computed as an approximation of u through suitable projection in V r at low complexity depending only on r.

We recall that, if the projection u r ∈ V r satisfies the following quasi-optimality condition

u(ξ) -u r (ξ) V ≤ C u(ξ) -P Vr u(ξ) V , ξ ∈ Ξ (3.4.5)
for some constant C independent from V r and ξ, and if the error estimate

∆(u r (ξ), ξ) is chosen such that u(ξ r ) -P Vr u(ξ r ) V ≥ γ sup ξ∈ Ξ u(ξ) -P Vr u(ξ) V , γ ∈ (0, 1], ( 3

Probabilistic and sample based framework

As the exact solution u(ξ) of (3.4.1) is not computable in general, a numerical approximation u n (ξ) in a finite dimensional subspace V n ⊂ V is usually considered. It is usually related to some discretization a priori for a given mesh of D of the boundary problem (3.4.1). Instead, we assume that we can only compute pointwise estimates u ∆t,M (x, ξ) of u(x, ξ) for any (x, ξ) ∈ Ξ × D. In view of practical applications, it rises several questions listed below.

1) How to compute the snapshots required for generating the reduced spaces V r and the approximation u r ?

2) How to choose a computable error estimator ∆(u r (ξ), ξ) in Algorithm 3.4.1 ?

The first point could be addressed by computing a snapshot as a global approximation u n (ξ) over D obtained through suitable projection in some finite dimensional vector space V n ⊂ V , from pointwise estimate u ∆t,M (x, ξ). To that goal we can, for example, consider a least-square or interpolation procedure in polynomial space as discussed in Section 3.3. Following the same path, it is possible to compute u r (ξ) the approximation of u(ξ) in V r . At this step, the goal is to compute such a projection using a small number of samples in comparison to seeking directly an approximation u n (ξ) in some approximation space V n .

Concerning the second point, the error estimate ∆(u r (ξ), ξ) required in Algorithm 3.4.1 is usually related to some norm of the equation residual in classical RB methods (see Section 1.3.1). Here, we follows another path by considering the L 2 -norm of the current approximation error e r (ξ) = u(ξ)u r (ξ) i.e. ∆(u r (ξ), ξ) = e r (ξ) 2 L 2 . By construction, u r (ξ) is a linear approximation in V r , regular enough (inherited from the snapshots), so that e r (ξ) can be written as the solution, for all ξ in Ξ, of

-A(ξ)e r (ξ) = g r (ξ) on D, e r (ξ) = f r (ξ) on ∂D, (3.4.7)
where f r (ξ) := f (ξ)u r (ξ) and g r (ξ) = g(ξ) + A(ξ)u r (ξ). By Theorem 3.

Probabilistic greedy algorithm

In this section, we present and analyze a probabilistic variant of Algorithm 3.4.1, in general setting, motivated by the applications discussed in the previous section. As a starting point, we assume that the error estimator required at the r-th iteration is under the form

∆(u r (ξ), ξ) := u(ξ) -u r (ξ) 2 V = E(Z r (ξ)), (3.4.10)
where Z r (ξ) is some parameter-dependent real valued random variable, defined on the probability space (Ω, F, P). Solving the following optimization problem,

ξ r ∈ arg max ξ∈ Ξ E(Z r-1 (ξ)), (3.4.11)
at step r of Algorithm 3.4.1, is in general out of reach, since E(Z r-1 (ξ)) is unknown a priori or too costly to compute.

In view of numerical applications, a first practical and naive approach is to seek ξ r as the maximum of a crude MC estimate (empirical mean) of E(Z r-1 (ξ)) , computed using i.i.d. samples of Z r-1 (ξ), ξ ∈ Ξ. Despite its simplicity, we have seen that such an estimate for the expectation suffers from low convergence with respect to the number of samples leading to possible high computational costs especially here where Z r-1 (ξ) can be expensive to evaluate. Moreover, nothing ensures that returned (random) parameter ξ r is a (quasi-)optimum for (3.4.11), almost surely or at least with high probability.

Probably Approximately Correct maximum

Instead, the so-called bandit algorithms (see, e.g., monograph [START_REF] Lattimore | Bandit Algorithms[END_REF]) are good candidates to address (3.4.11). In particular probably approximately correct (PAC) bandit algorithms are designed not to return an optimum almost surely but a quasi-optimum, either in absolute or relative precision, with high probability. Here, we particularly focus on PAC bandit algorithms that return ξ r (random) as a PAC maximum in relative precision for E(Z r (ξ)) over Ξ, i.e.,

P (E(Z r-1 (ξ r )) -E(Z r-1 (ξ r )) ≤ εE(Z r-1 (ξ r ))) ≥ 1 -λ, (3.4.12)
for any prescribed threshold ε in (0, 1), probability of failure λ in (0, 1) and

ξ r ∈ arg max ξ∈ Ξ E[Z r-1 (ξ)].
In what follows, we denote ξ r := PAC λ,ε (Z r-1 , Ξ) a PAC maximum in relative precision for (3.4.11), if it satisfies (3.4.12).

Practical algorithms Classical PAC bandit algorithms are aimed at returning a PAC maximum in absolute precision [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF]. Here, we are rather interested in having a PAC maximum in relative precision as proposed in [START_REF] Billaud-Friess | A PAC algorithm in relative precision for bandit problem with costly sampling[END_REF]. In this paper, we studied two variants : a naive and an adaptive one. The interest of the adaptive algorithm is that it outperforms the mean complexity of the naive algorithm in terms of number of required samples. Such an algorithm preferentially samples the random variable Z r (ξ) for the parameter values for which it is more likely to find a maximum. Thus it is particularly well suited for applications with high sampling cost. Both algorithms have been proven [MBF15, Propositions 3.1 & 3.2] to provide PAC solution in relative precision for a family of random variables {Z r (ξ), ξ ∈ Ξ}, indexed by ξ, satisfying suitable concentration inequalities under the form (3.2.2).

A weak greedy algorithm with high probability

Now, we present a probabilistic weak greedy algorithm where the parameter ξ r is a PAC maximum in relative precision at each step r. Applications of the probabilistic reduced basis framework. The approach, discussed in Section 3.4.1, has been motivated by the solution of PDEs with probabilistic representation. It is none the less usable for the general problem of approximating any family of functions (possibly random) {u(ξ) : ξ ∈ Ξ}, indexed by the parameter ξ. In particular, if u(ξ) is known a priori but costly to evaluate, the proposed probabilistic greedy algorithm can be used to generate a sequence of spaces V r and corresponding interpolant u r of u onto V r . Then, the proposed method leads to a probabilistic EIM as in [START_REF] Cai | A Stochastic Discrete Empirical Interpolation Approach for Parameterized Systems[END_REF]. In particular, in a fully discrete setting, u can be identified with a matrix and the resulting algorithm is a probabilistic version of adaptive cross approximation for low-rank matrix approximation [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF][START_REF] Tyrtyshnikov | Incomplete Cross Approximation in the Mosaic-Skeleton Method[END_REF], where the proposed PAC bandit algorithm provides a particular column-selection strategy that sparsely explores the matrix entries.

Time dependent problems. The proposed approaches were mainly applied for the approximation of linear (parameter-dependent) elliptic problems. However, Feynman-Kac representation formula still holds for parabolic problems (see e.g. [START_REF] Comets | Calcul stochastique et modèles de diffusions[END_REF][START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF]). Thus, applying for example the proposed interpolation method, or a learning approach with low-rank tensor format as previously discussed, for dealing with time dependent high-dimensional and/or parameter-dependent problem approaches is straight-forward. However, extending the presented probabilistic and sampling based RB paradigm for time-dependent and parameter-dependent problems (as discussed in Chapter 2) is a challenging and open question.

Past research contributions

In this chapter we give a concise summary of the contributions [MBF2, MBF4, MBF3, MBF5, MBF6, MBF7] concerning the development of numerical methods for the simulation of multicomponent or multi-phase flows.

Compressible-incompressible interface flows

In the literature, several approaches are available to treat diphasic flows constituted of immiscible gaz and liquid. Some of them consider both fluids as e.g. incompressible [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows[END_REF] or compressible [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF]. However, in one approach the one does not ensure the volume variation of the gaz whereas the second does not respect the volume conservation of the liquid. As [START_REF] Caiden | A numerical method for two phase flow consisting of separate compressible and incompressible regions[END_REF], we consider an approach that consider the gaz as compressible fluid and the liquid as an incompressible fluid. To deal numerically with such interface flows, one faces at least to the following difficulties: (1) transition from a gas density linked to the local temperature and pressure by an equation of state to a liquid density mainly constant in space, (2) proper approximation of the divergence constraint in incompressible regions and (3) wave transmission at the interface.

We proposed in [MBF2, MBF3] a global numerical approach (i.e. the same for each phase) that allow to easily simulate compressible-incompressible on unstructured meshes (allowing complex geometry) which possible high order precision. This work was related to my PhD thesis

Global model

The proposed numerical approach relies on the formulation of a global model using the same set of primitive unknowns and equations used everywhere in the flow, but with a dynamic parameterization that changes from compressible to incompressible regions.

Flows under consideration are viscous, Newtonian and consequently described by the Navier-Stokes equations. On one hand, the compressible Navier-Stokes equations are considered under weakly compressibility assumption so that a non-conservative formulation can be used. On the other hand, the incompressible non-isothermal model is retained. In addition, the level set transport equation is used to capture the interface position needed to identify the local characteristics of the fluid and to recover the adequate local modelling.

The first step of unification consists in choosing a same unknown suitable for both compressible and incompressible models. This is possible for both primitive and entropic unknowns. For flows without shocks, we consider quasi-linear unified equations [START_REF] Hauke | A unified approach to compressible and in-compressible flows[END_REF][START_REF] Hauke | A comparative study of different sets of variables for solving compressible and incompressible flows[END_REF], for primitives unknowns Y Y Y = (p, u t , T) t (with p the pressure, u the velocity and T the temperature)

L(Y, χ, ρ, )Y = S(Y, ρ), (1)
∂ x i K ij (Y, µ, κ)∂ x j .
Here, A 0 , A i correspond the jacobian matrices derived from Navier-Stokes equations under their conservative form. The K ij are the diffusion matrices (see [MBF3, Section 2.2]). Here, the jacobian matrices A 0 , A i depend on three quantities namely the density ρ, and two thermodynamical coefficients

χ = (α p , β T ) = - 1 ρ ∂ρ ∂T p , 1 ρ ∂ρ ∂p T
that represent respectively the volumic dilation effects (α p ) and isothermal compressibility (β T ).

In the particular case of perfect gaz (indicated with subscript g ) we have

χ g = 1 T , 1 p and ρ g = p RT , (2)
whereas for an incompressible (indicated with subscript l ) fluid it holds

χ l = (0, 0) et ρ l = constant. (3)
The diffusion matrices also depend on the viscosity coefficient µ and, κ the thermal conductivity which are different for each fluid. For the liquid (gas), they are noted µ l , κ l (respectively µ g , κ g ).

The interface is captured through Level Set method [START_REF] Osher | Fronts propagating with curvature dependant speed: algorithms based on based on Hamilton-Jacobi formulations[END_REF]. It is given as the zero isovalue of a continuous function φ, called Level Set function that satisfies φ < 0 in the liquid and φ > 0 in the gaz. Here, we consider diphasic flow without any state change, so that the Level Set φ satisfies the following transport equation

∂ t φ + u • ∇φ = 0, ( 4
)
where u is the fluid flow velocity.

Combining (1)-( 4), we derive the following global system of equations

L(Y Y Y , φ)Y Y Y = S S S(Y Y Y , φ), (5)
∂ t φ + u • ∇φ = 0, (6)
where L is defined by In the previous definition H(φ) = 1 φ≥0 is the characteristic function equals to 1 in the gas.

L(Y Y Y , φ) = A 0 (Y Y Y , φ)∂ t + d i=1 A i (Y Y Y , φ)∂ x i - d i,j=1 ∂ x i K ij (Y Y Y , φ) ∂ x j ,

Global finite element approximation

The resolution of ( 5)-( 6) relies on a time splitting. A each time step, we first solve the global equation ( 5) then the interface position is updated through [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF]. The two relates subproblems are solved independently with different numerical schemes but with the same mesh in time and space (thus the time step fulfills a global CFL condition). As they are adapted to treat complex geometry and provide high order approximation, Stabilized Finite Element (SFE) have been used for the space discretization of (5) together with explicite Euler scheme in time. A Runge Kutta / Discontinous Galerkin (RKDG) [START_REF] Di Pietro | Mass preserving finite element implementations of the level set method[END_REF][START_REF] Marchandise | A quadrature-free discontinuous Galerkin method for the level set equation[END_REF] has been considered for solving the transport equation on the level set [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF]. The combination, a posteriori, of these two solvers has been done following the procedure proposed in [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows[END_REF] for incompressible and diphasique flows (see details in [START_REF] Billaud | A simple stabilized finite element method for solving two phase compressible-incompressible interface flows[END_REF]). Here, we only present the "coupled" and "global" numerical methop used for the numerical solution of (5).

For single fluid flow Lagrange FE method are no suitable for numerical solution of ( 1) are they provides unstable approximations in case of large advection, and violate the inf-sup condition (for the pressure). To overcome these limitations, it is possible to add to the Galerkin formulation some additional term called stabilization. To that goal, a Galerkin Least Square (GLS) [START_REF] Hauke | A unified approach to compressible and in-compressible flows[END_REF][START_REF] Hauke | A comparative study of different sets of variables for solving compressible and incompressible flows[END_REF] method is retained since it allows to deal with those problems for both incompressible and compressible equations. GLS method is a Petrov-Galerkin method, where the test fonction is W W W h plus an additional term LW W W h multiplied by a matrix containing coefficients τ e . These coefficients defined locally on each element e of the mesh act as stabilization term. This matrix is defined differently depending on the compressible [START_REF] Hauke | Simple stabilizing matrices for the computation of compressible flows in primitive variables[END_REF] or incompressible case [START_REF] Franca | Stabilized finite element methods : II. The incompressible Navier-Stokes equations[END_REF]. Thus, the resulting numerical scheme is said to be unified, as it allows to solve numerically both compressible and incompressible Navier-Stokes equations.

This unified FE scheme is a good candidate for numerical solution of (5), however it requires some particular care for the treatment of the elements crossed by the interface (thus containing incompressible and compressible fluid). Indeed, it requires the computation of "discontinuous" integrals (i.e. where the quantity to integrate is discontinuous due to density, viscosity and conductivity) that allows correct wave transmission at the interface and the derivation of a suitable stabilization term adapted to each phase. To this end, we use averaged quantities in each element e relying on the volume fraction of the compressible fluid in each element noted ε n e (obtained using φ). We take • arithmetical average for the unknown Y Y Y n e , the matrix of stabilization τ n e , and for the the thermodynamical coefficients χ χ χ n e , • while harmonic average is retained for the density ρ n e , viscosity µ n e and conductivity κ n e . We finally get a general numerical method relying on the following global discrete variational formulation. Let be given Y Y Y n h at t n and the interface position through φ n , Equation ( 5) is solved from t n to t n+1 for 1 ≤ n ≤ N -1 by seeking Y Y Y n+1 h ∈ S p h × S u u u h × S T h such that for all W W W h ∈ T p h × T u u u h × T T h we have with τ e,g , τ e,l the stabilization matrices for a compressible and incompressible fluide respectively.

                      

Numerical results

First numerical results have concerned numerical simulation of smooth flow for which the global behavior is driven by the gaz compressibility. For the considered cases, comparisons with analytical computation of relevant quantity related to the flow permitted the validation of the proposed approach as well as convergence study, in both one and two dimensions. Then more general test cases have been considered as rising bubble or oscillation water column (see [START_REF] Billaud | A simple stabilized finite element method for solving two phase compressible-incompressible interface flows[END_REF] for details) .

Multi-material compressible flows

My post-doctorate was dedicated to the simulation of multi-material compressible flows arising in astrophysics ou inertial confinement fusion (ICF). Several approaches are available for numerical simulation of such flows, using either interface capturing or tracking methods. Here we have adopted subsequently each point of view leading to several publications [MBF4, MBF5, MBF6, MBF7]. These contributions are related to several collaborations with the CEA with Edouard Audit (Maison de la simulation, CEA), Jérôme Breil (CEA-CESTA, Celia), Guy Schurtz (CEA-CESTA, Celia, retired) and Samuel Kokh (CEA, Saclay) but also during research stay during CEMRACS 2010 with the team SIMCAPIAD.

Eulerian methods: We consider problems restricted to multi-material flows with m component on cartesian mesh. In this framework, some methods [169,[START_REF] Allaire | A five-equation model for the numerical simulation of interfaces in two-phase flows[END_REF] consider a global model to represent the flow with particular Equation of States (EOS). Chosing suitable EOS allows to treat the diffusion region (it is only numerical i.e. not related to any physics) appearing near the interface when solving numerically those equations (e.g. using finite volume (FV)). In particular, it allows the numerical solution to converge to the correct solution. Some approaches have been developped to limitate the interface diffusion [START_REF] Després | Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire[END_REF][START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] using anti-diffusive numerical scheme. In [169,[START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluides by means of a five-equation model[END_REF], a 5 equation model (quasi-conservative) is proposed to modelize the evolution of m = 2 compressible and immiscible fluides separated by an interface. It corresponds to Euler equations combined to partial mass equations, and the transport of a color function.

In [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluides by means of a five-equation model[END_REF], this system is solved using anti-diffusive FV scheme, in order to capture accurately the interface (less numerical diffusion). In that goal, a splitting is proposed 1) the acoustic part of the equations is solved using a classical FV scheme (e.g. relaxation, acoustic scheme), 2) the transport is solved using the anti-diffusive scheme [START_REF] Després | Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire[END_REF][START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF]. This approach is only proposed to two materials compressible flows and is only first order accurate.

1. In [START_REF] Billaud-Friess | A second order anti-diffusive Lagrange-remap scheme for two-component flows[END_REF] we have proposed to improve the accuracy of this method to second order using MUSCL resonstruction scheme in space (for acoustic fluxes) together with Heun scheme in time. The approach has been validated in 1D, and 2D dimensional case have been treated using directional splitting.

) 8
where ρ is the density, u the velocity, p the pressure, E total energy and the Y k are mass fractions for each fluid k, 1 ≤ k ≤ m. This system, together with properly chosen equation of states is hyperbolic. Following [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluides by means of a five-equation model[END_REF], we proposed to apply a lagrange-projection scheme with an anti-diffusive scheme for the projection step when solving [START_REF] Bachmayr | Sparse polynomial approximation of parametric elliptic pdes. Part ii: lognormal coefficients[END_REF]. In that case, the main difficulty was to generalize the proposed scheme to m > 2 color function Z k preserving the maximum principle and such that the sum is equal to 1. To that goal, a recursive construction of the fluxes for the color function has been proposed in the line of [START_REF] Jaouen | Numerical transport of an arbitrary number of components[END_REF] . The resulting method has been tested and validated on several test cases in 1D and 2D for flows with 3 and 5 materials [START_REF] Billaud-Friess | An anti-diffusive Lagrange-Remap scheme for multimaterial compressible flows with an arbitrary number of components[END_REF][START_REF] Billaud-Friess | Simulation of sharp interface multi-material flows involving an arbitrary number of components through an extended five-equation model[END_REF] for different equations of state.

Contribution on ALE method: Here we focus on the simulation of multi-material compressible flows on unstructured meshes. In that context, we consider an interface capturing approach that consists in explicit reconstruction of interfaces from color functions see e.g. VOF (volume of fluid) method (m = 2) or MOF (moment of fluid) method (m ≥ 2) which are very efficient when combined with ALE (Arbitrary Lagrangian Eulerian) method to track the interface. The latter allows to combine the advantages of an eulerian approach adapted for complex flows with distortion and a lagrangian approach that allows to exactly capture material fronts as interfaces. Here, we retained a CCALE-MOF (Cell-Centered Arbitrary Lagrangian Eulerian MOF), where the color functions are the fluid moment (e.g. mass fraction) [START_REF] Galera | A two-dimensional unstructured cell-centered multmaterial ALE scheme using VOF interface reconstruction[END_REF][START_REF] Galera | A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction[END_REF]. This algorithm works according to the following steps.

1. Lagrangian step: First, the mesh as well as conservative unknowns (ρ, ρu, ρE) are updated during a Lagrangian step using a centered VF: Explicit Unstructured Cell-Centered Lagrangian HYDrodynamics (EUC-CLHYD). The color functions (e.g. mass fraction) remain unchanged during this step.

2. Reconstruction step: Interfaces are reconstructed from fluid moment allowing perfect recovering of interfaces.

3. Rezoning step: During this step the mesh is "reshaped", i.e, if the mesh is very distorted during lagrangian step it is rectifed.

4. Projection step: Finally, a projection step is done. It constists in interpolating in conservative way the physical quantities on the new mesh.

The CCALE-MOF approach, works well for cartesian coordinates and can be extended quite easily for cylindrical coordinates. The objective of the work [START_REF] Billaud-Friess | A Multi-Material CCALE-MOF Approach in Cylindrical Geometry[END_REF] was to develop some "hybrid" approach valid for both coordinates with axisymmetric geometries. The method was tested with success for the numerical simulation of complex multimaterial flows with axisymmetric geometries for cartesian, polar or unstructured meshes.

Figure 1 . 2 . 1 :

 121 Figure 1.2.1: Comparison of minimal residual methods for M ≤r (X) = T can r (X) and • X = • 2 . Convergence with the rank r of the approximations obtained with CMR and A-IMR with different precisions δ, and with direct (left) or greedy rank-one (right) approaches.

Figure 1 . 2 . 2 :

 122 Figure 1.2.2: Convergence with the rank of approximations obtained with the greedy CMR and A-IMR algorithms for different precisions δ. Convergence is plotted with respect to the norm • 2

  3.4 provides a sharper bound than Proposition 1.3.2 for the error on the QoI by taking advantage of the orthogonality property (1.3.15).

Proposition 1 . 3 . 5 .

 135 3.21). The following proposition [MBF10, Proposition 2.14] provides an error bound for the approximation of the QoI. Let (ũ, ỹ) := (u r,p , y r,p ) ∈ V r × T p be the solution of the saddle point problem (1.3.21). The approximation s defined by (1.3.22) satisfies

10 - 6

 6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 PDF of the error. Three different preconditioners P m (ξ) are used: m = 0 (dotted lines), m = 10 (dashed lines) and m = 30 (continuous lines). Primal only L ∞ -norm L 2 -norm m = 0 1.284 × 10 0 1.245 × 10 -1 m = 5 1.203 × 10 0 9.637 × 10 -2 m = 10 1.458 × 10 0 1.064 × 10 -1 m = 20 1.068 × 10 0 8.386 × 10 -2 m = 30 1.066 × 10 0 7.955 × 10 -2 Primal-dual L ∞ -norm L 2 -norm m = 0 2.751 × 10 -1 1.085 × 10 -2 m = 5 1.308 × 10 -1 5.708 × 10 -3 m = 10 1.333 × 10 -1 5.807 × 10 -3 m = 20 1.232 × 10 -1 5.465 × 10 -3 m = 30 1.224 × 10 -1 5.408 × 10 -3 Saddle point L ∞ -norm L 2 -norm m = 0 1.023 × 10 -1 4.347 × 10 -3 m = 5 9.715 × 10 -2 3.389 × 10 -3 m = 10 9.573 × 10 -2 3.867 × 10 -3 m = 20 6.022 × 10 -2 2.996 × 10 -3 m = 30 5.705 × 10 -2 2.896 × 10 -3 (b) L ∞ and L 2 norm of the error.

Figure 1 . 3 . 1 :

 131 Figure 1.3.1: Application 2: Probability density function, L ∞ norm and L 2 norm of the error s(ξ)s(ξ) Z estimated from a samples set of size 10 4 .

  Figure 1.3.2 shows statistics of the effectivity index η(ξ) = ∆(ξ)/ s(ξ)s(ξ) Z for different numbers m of interpolation points for the preconditioner.

  PDF of η(ξ) for the primal-dual methods and the saddle point methods. Three different preconditioners P m (ξ) are used: m = 0 (dotted lines), m = 10 (dashed lines) and m = 30 (continuous lines)

Figure 2 . 2 . 1 :

 221 Figure 2.2.1: Advection equation: Statistical estimations of the expectation and maximum of the relative errors E 2 and E ∞ with respect to the reduced dimension r, for MTD (left) and MTI (right) for both continuous (top) and discontinuous (bottom) initial conditions.

Figure 2 . 2 . 2 :

 222 Figure 2.2.2: Burger's equation: Statistical estimations of the expectations and maxima of the relative errors E 2 and E ∞ with respect to the reduced dimension r, for MTD (left) and MTI (right).

. 17 )

 17 with the initial data u 0 (ξ) : Ω → R and supplemented with homogeneous Dirichlet boundary conditions. The solution u(•; ξ) depends on the parameter ξ ∈ R 3 through the viscosity µ(ξ) = ξ 1 , the initial condition and the source term defined by means of the function f (ξ) : Ω × I → R given by f (x, t; ξ) = ξ 2 exp(-(x -0.2) 2 /0.03 2 ) sin(ξ 2 πt)1 [0.1,0.3] (x).

Figure 2 . 3 . 2 :

 232 Figure 2.3.2: Burgers's equation: evolution of e k for both algorithms for r ∈ {5, 8, 10, 12, 15, 20}.

Figure 2 . 3 . 3 :

 233 Figure 2.3.3: Burgers's equation: error e K at final time for different (r, ∆t).

Figure 2 . 3 . 4 :

 234 Figure 2.3.4: Burgers's equation (multiple parameter case): error e K at final time for different r.

. 7 )

 7 where b : R d → R d and σ : R d → R d×d are the drift and diffusion coefficients, respectively.Before recalling Feynman-Kac formula, we first need to introduce additional assumptions and notations. Denoting by • both euclidean norm on R d and Frobenius norm on R d×d , we introduce the following assumption on b and σ.Assumption 1 (A1).There exists a constant 0 < M < +∞ such that for all x, y ∈ D we haveb(x)b(y) + σ(x)σ(y) ≤ M xy , (3.2.8) Under Assumption 1, by [57, Theorem 5.3], there exists an unique strong solution to the Equation (3.2.7).

Theorem 3 . 2 . 1 (

 321 Feynman-Kac formula). Under assumptions (A1)-(A4) there exists a unique solution of (3.2.5) in C D ∩ C 2 (D), which satisfies for all x ∈ D u

Theorem 3 . 3 . 1 .

 331 Assuming (A2)-(A4), it holds

Theorem 3 . 3 .

 33 1 can be extended [MBF12, Corollary 1] to provide a convergence result in L ∞ (D) using Lebesgue constants in L ∞ -norm associated to I Λ . Corollary 3.3.2 (Convergence in L ∞ ). Assuming (A2)-(A4), it holds lim sup n→∞

. 4 )

 4 with the Lebesgue constant L Λ = sup x∈D ν∈Λ |ϕ ν (x)|.

4 . 1 . 3 . 4 . 1 (

 41341 Algorithm Deterministic greedy algorithm). Let Ξ ⊂ Ξ be a finite training set and V 0 = {0}. For r ≥ 1 proceed as follows.(Step 1.) Select ξ r ∈ arg max ξ∈ Ξ ∆(u r-1 (ξ), ξ).

.4. 6 )

 6 thus the Algorithm 3.4.1 is a weak greedy algorithm with respect to the discrete solution manifold K = {u(ξ) : ξ ∈ Ξ}.

2 . 1

 21 for all ξ in Ξ, e r (ξ) is the unique solution of(3.4.7) in C(D) ∩ C 2 (D) and satisfies for allx ∈ D e r (x, ξ) = E F r (x, X x,ξ , ξ) = E f r (X x,ξ τ x,ξ , ξ) + τ x,ξ 0 g r (X x,ξ t , ξ)dt . (3.4.8) Then, we have [123, Theorem 5.7] the following probabilistic reinterpretation of e r (ξ) 2 L 2 . Theorem 3.4.2. Taking Y ∼ U(D) uniformly distributed on D, we have for any ξ in Ξ e r (ξ) 2 L 2 = |D|E (Z r (ξ)) .(3.4.9)Here Z r (ξ) = F r (Y, X Y,ξ , ξ)F r (Y, XY,ξ , ξ) where X x,ξ and Xx,ξ are two i.i.d random stopped diffusion processes solution of (3.4.3), for any x ∈ D, and |D| is the Lebesgue measure of D.

Algorithm 3 . 4 . 3 (

 343 Probabilistic greedy algorithm). Let Ξ ⊂ Ξ be a discrete training set and (λ r ) r≥1 ⊂ (0, 1). Starting from V 0 = {0}, proceed, for r ≥ 1, as follows.(Step 1.) Select ξ r := PAC λr,ε (Z r-1 , Ξ) (Step 2.) Compute u(ξ r ) and update V r = span{u(ξ 1 ), . . . , u(ξ r )}.

  [START_REF] Billaud | Stabilized finite element method for incompressible-compressible two-phase flows[END_REF] made at CEA-CESTA (Comissariat à l'énergie atomique et aux énergies alternatives, Centre d'études scientifiques et techniques d'Aquitaine), under the supervision of Jean-Paul Caltagirone (Trèfle, I2M, Université de Bordeaux), Gérard Gallice (CEA-CESTA, retired), Boniface Nkonga (Laboratoire J.A. Dieudonné, Université de Nice Sophia Antipolis).

  withL(Y, χ, ρ, µ, κ) = A 0 (Y, χ, ρ)∂ t + d i=1 A i (Y, χ, ρ)∂ x i -

A 0

 0 (Y Y Y , φ) = A 0 (Y Y Y , χ χ χ(φ), ρ(φ)), A i (Y Y Y , φ) = A i (Y Y Y , χ χ χ(φ), ρ(φ)), K ij (Y Y Y , φ) = K ij (Y Y Y , µ(φ), κ(φ)), S(Y Y Y , φ) = S S S(Y Y Y , ρ(φ)), (φ) = (1 -H(φ))χ χ χ + H(φ)χ χ χ g , ρ(φ) = (1 -H(φ))ρ + H(φ)ρ g , µ(φ) = (1 -H(φ))µ + H(φ)µ g , κ(φ) = (1 -H(φ))κ + H(φ)κ g ,

h

  is the FE approximation at t n and W W W h the test fonction. Finally the average matrices A

2 .∂

 2 Then, to deal with flows with a arbitrary number of component, we have proposed an extension of the five equation model[169], leading to the following m-component modelt (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, moment conservation ∂ t (ρE) + ∇ • ((ρE + p)u) = 0, energie conservation ∂ t (ρY k ) + ∇ • (ρuY k ) = 0, partial mass conservation ∂ t Z k + u • ∇Z k = 0, color function transport, (

  . For a Hilbert space H, we denote by R H ∈ L(H, H ) the Riesz map such thatv 2 H = R H v, v , where •, • denotes the duality pairing. The dual norm • H on H is such that R H = R -1 H . Then v H = R H v H and | v, w | ≤ v H w H holdfor any v ∈ H and w ∈ H . For any operator C ∈ L(H 1 , H 2 ), with H 1 and H 2 two Hilbert spaces, C * ∈ L(H 2 , H 1 ) denotes the adjoint of C, such that Cv 1 , v 2 = v 1 , C * v 2 for any v 1 ∈ H 1 and v 2 ∈ H 2 .

It means that v V0 = v h H 1 (D) for all v ∈ V , where v h = n i=1 v i ϕ i .

For any A ∈ R n×m , the Moore-Penrose pseudo inverse is given by A = (A T A) -1 A T .

Here G r (R p ) = {V r ⊂ R p : dim(V r ) = r} denotes the Grassmann manifold.

Numerical estimation of the solution of a PDEs

We now discuss the numerical computation of an estimate of u(x). We perform a suitable integration scheme, to get an approximation of the stopped diffusion process X x with a MC method to evaluate the expectation in formula (3.2.10).

Integration of diffusion process

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

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

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

Such a discretization choice may leads to over-estimation of the exit time with an error in O(∆t 1/2 ). More sophisticated approaches [START_REF] Gobet | Monte-Carlo methods and stochastic processes: from linear to non-linear[END_REF]Chapter 6] exist to improve the order of converge, as Brownian bridge, boundary shifting or Walk On Sphere (WOS) methods.

Monte Carlo estimate

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

Approximation error

Introducing the following pointwise estimate,

the error is the sum of two contributions,

where ε ∆t is the time integration error and ε M C is the Monte-Carlo estimation error.

Under assumptions (A1)-(A4), in particular assumption (A3-ii) which ensures that D does not have singular points, it can be proven [86, §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

Adaptive sparse (polynomial) interpolation

When D admits some cartesian product structure, multivariate interpolation can be performed on tensor product of univariate bases. For high dimensional problems, there can be some advantage in considering approximation spaces that are sparse and anisotropic with respect to the different variables. To that goal, adaptive selection of the interpolation points can be envisaged. In such a procedure the set of interpolation points is successively enriched, as the polynomial dimension increases, at different stages of the computation. In what follows, we recall in brief the main elements of sparse polynomial interpolation following [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF][START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs[END_REF].

Notations Here, we assume that the function to interpolate u is defined on

k } k∈N where the ϕ

k are univariate polynomial of degree k defined on D i . For a multi-index ν = (ν 1 , . . . , ν d ) ∈ N d , we introduce the multivariate function

Let Λ ⊂ N d downward closed i.e. for all ν ∈ Λ, such that µ ≤ ν then µ ∈ Λ. 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

Main idea The goal of sparse interpolation is to provide a sequence of approximations {I Λn (u)} n≥1 , in increasing polynomial spaces {P Λn } n≥1 , associated with a nested sequence of downward closed subsets {Λ n } n≥1 . The sequence of downward closed subsets (Λ n ) n≥1 is selected such that the best approximation error and the Lebesgue constant satisfy

In practice, 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] provides a practical algorithm that constructs a good sequence of sets Λ n .

Proposed algorithms

We have proposed (see [START_REF] Billaud-Friess | Stochastic Methods for Solving High-Dimensional Partial Differential Equations[END_REF] for details) two algorithms combining sparse adaptive interpolation and sequential control variate procedure discussed previously. Here we give the general guidelines of the proposed algorithms.

1. The first algorithm called Perturbed Adaptive Sparse Interpolation algorithm. It is a perturbed version of the sparse adaptive algorithm where at each step n the computation of the exact interpolant I Λn (u) from u(x ν ) is replaced by an approximation denoted u n in V n = P Λn based on MC estimate u ∆t,M (x ν ). In particular, this approximation approximation is computed with the sequential control variate procedure stopped for a given number of iterations or according to some stopping criterion.

2. As a second algorithm, we consider an Adaptive Sequential Control Variates algorithm where an approximation e k n k of e k is obtained by applying the adaptive interpolation algorithm 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, e k n k therefore belongs to a different approximation space V n k = P Λn k .

Numerical applications

In this section, we give a short illustration of the behavior of the perturbed adaptive sparse algorithm for polynomial interpolation in high dimension. A detailed study of the different algorithms discussed previously is available in [MBF12, Section 4.3].

Let us consider the following simple diffusion equation

were

The source terms and boundary conditions will be specified later for each test case. The stochastic differential equation associated to (3.3.5) is the following

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.

We consider two test cases noted respectively (TC1) and (TC2).

(TC1) As first test case in dimension d = 10, we define (3.3.5) such that its solution is the Henon-Heiles potential

(TC2) We also consider the problem (3.3.5) whose exact solution is a sum of non-polynomial functions in dimension d = 20, given by u(x) = x 2 1 + sin(x 12 ) + exp(x 5 ) + sin(x 15 )(x 8 + 1).

MC estimate are computed for (∆t, M ) = 10 -4 , 1000) and use set to 30 the maximal number of iterates in the controle variate procedure. Since for both test cases the exact solution is known, we propose to compare the behavior exact and perturbed sparse adaptive algorithms. In Table 3.1 and Table 3.2, we summarize the results associated to both algorithms for (TC1) and (TC2) respectively. We observe that the perturbed sparse adaptive algorithm performs well in comparison to sparse adaptive algorithm, for (TC1). Indeed, we get an approximation with a precision below the prescribed value ε for both algorithms. Similar observation holds for (TC2) in Table 3.2 and this despite the fact that the test case involves higher dimensional problem.

Parameter-dependent PDEs

Now, we consider a parameter-dependent variant of the model problem (3.2.5). Let Ξ ⊂ R p be a set of parameter. Keeping the same notations as in Section 3.2, the problem of interest is

where the boundary and source terms f and g as well as the linear elliptic operator A (through the coefficients b, σ in (3.2.6)), depend on ξ. Also, we assume that, for any ξ ∈ Ξ, there is a unique solution u(ξ) to (3.4.1). Moreover, if for all parameter ξ ∈ Ξ, the assumptions (A1)-(A4) are fulfilled, by Theorem 3.2.1, the solution u(ξ) admits a probabilistic representation under the form

Here, X x,ξ is a parameter-dependent diffusion process, defined for all t ≥ 0 by

where W stands for the standard brownian motion. A MC estimate of u(x, ξ) is then

for m-independent realisations of X x,∆,ξ the approximated diffusion process obtained through Euler-Maruyama scheme.

We have proven in [START_REF] Macherey | Approximation and model reduction for partial differential equa-tions with probabilistic interpretation[END_REF]Chapter 5] that Algorithm 3.4.3 is a weak greedy algorithm with high probability.

Theorem 3.4.4. Take (λ r ) r≥1 ∈ (0, 1) such that r≥1 λ r = λ < 1, ε ∈ (0, 1) and Ξ ⊂ Ξ a discrete training set. Moreover, suppose that for r ≥ 1, the approximation u r of u in V r satisfies the following quasi-optimality condition (3.4.5). Then, Algorithm 3.4.3 is a weakgreedy algorithm of parameter √ 1-ε C , with probability at least 1λ i.e.

P u(ξ r ) -

General comments

• At a step r of the Algorithm 3.4.3, the reduced space V r = span{u(ξ 1 ), . . . , u(ξ r )}, as well as the approximation u r (ξ) are no longer deterministic. Indeed, they are related to the selected parameters ξ 1 , . . . , ξ r depending themselves on the errors at the previous steps through i.i.d. samples of the random variables Z i (ξ) for all ξ ∈ Ξ and i < r (required during PAC selection of ξ r ). Thus, the result stated in Theorem 3.4.4 holds with high probability instead of with certainty, which is the price to pay when introducing some randomness in the greedy procedure. Especially, it demonstrates that Algorithm 3.4.3 is a weak greedy algorithm, with high probability, with respect to the discrete solution manifold K = {u(ξ) : ξ ∈ Ξ}.

• In the lines of [START_REF] Cohen | Reduced Basis Greedy Selection Using Random Training Sets[END_REF], it is also possible to consider a fully probabilistic variant of Algorithm 3.4.3, in which a training set Ξ r randomly chosen is used at each step r of Algorithm 3.4.3 instead of Ξ. For a particular class of functions, of ξ, that can be approximated by polynomials with a certain algebraic rate, it can be proven that, for suitable chosen size of random training set Ξ r , the resulting algorithm is a weak greedy algorithm with high probability with respect to the continuous solution manifold K = {u(ξ) : ξ ∈ Ξ} (some elements of proof are given in [123, Chapter 5]).

Summary and discussion for future research

In this chapter, we have addressed the problem of approximating the solution of PDEs, whose solution admits a probabilistic representation. Within this framework, we have proposed a probabilistic sparse interpolation method for high dimensional problems.

Ongoing work concerns some probabilistic reduced basis method using samples of the function to approximate. Beyond the validation of the proposed approach, one pending and challenging question that originally motivated this work concerns the feasibility of such an approach for parameter-dependent problems in high dimension.

In what follows we give additional elements concerning possible extensions to these contributions.

Learning approaches with tensor using F-K representation for PDEs. Several recent works on developing algorithms that are learning approaches in tensor format i.e. that use random and unstructured evaluation of the function to approximate have been proposed see e.g. [START_REF] Grelier | Learning with tree-based tensor formats[END_REF][START_REF] Haberstich | Active learning of tree tensor networks using optimal least-squares[END_REF][START_REF] Nouy | Higher-order principal component analysis for the approximation of tensors in tree-based low rank formats[END_REF] or [START_REF] Eigel | Variational Monte Carlo-bridging concepts of machine learning and high-dimensional partial differential equations[END_REF]. Within this framework, an interesting idea would be to combine it with Feynman-Kac pointwise estimations of the solution to perform suitable approximation of high dimensional PDEs.