C Soize
email: c.soize@univ-eiffel.fr

R Ghanem
email: ghanem@usc.edu

Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets

Keywords: probabilistic learning, PLoM, partial differential equations, small dataset, non-Gaussian, statistical constraints, stochastic computational models, Navier-Stokes equations, nonlinear solid dynamics, machine learning, data driven, data science, Kullback-Leibler, minimum cross-entropy, uncertainty quantification

A novel extension of the Probabilistic Learning on Manifolds (PLoM) is presented. It makes it possible to synthesize solutions to a wide range of nonlinear stochastic boundary value problems described by partial differential equations (PDEs) for which a stochastic computational model (SCM) is available and which depend on a vector-valued random control parameter. The cost of a single numerical evaluation of this SCM is assumed to be such that only a limited number of points can be computed for constructing the training dataset (small data). Each point of the training dataset is made up of realizations from a vector-valued stochastic process (the stochastic solution) and the associated random control parameter on which it depends. The presented PLoM constrained by PDE allows for generating a large number of learned realizations of the stochastic process and its corresponding random control parameter. These learned realizations are generated so as to minimize the vector-valued random residual of the PDE in the mean-square sense. Appropriate novel methods are developed to solve this challenging problem. Three applications are presented. The first one is a simple uncertain nonlinear dynamical system with a nonstationary stochastic excitation. The second one concerns the 2D nonlinear unsteady Navier-Stokes equations for incompressible flows in which the Reynolds number is the random control parameter. The last one deals with the nonlinear dynamics of a 3D elastic structure with uncertainties. The results obtained make it possible to validate the PLoM constrained by stochastic PDE but also provide further validation of the PLoM without constraint.

Introduction

(i) Context. Physics-informed machine learning (eg. neural networks and deep learning) has seen a growing interest in recent years. Numerous works in this area of machine learning have been carried out for Schrödinger equation [START_REF] Mills | Deep learning and the schrödinger equation[END_REF], for reconstructing Reynolds stress model in turbulent flow [START_REF] Wang | Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on DNS data[END_REF] and for augmenting turbulence model [START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF][START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework[END_REF], for predicting large eddy simulation [START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF], for projection-based model reduction [START_REF] Swischuk | Projection-based model reduction: Formulations for physics-based machine learning[END_REF]. Physics-informed machine learning has also been used in the context of linear [START_REF] Raissi | Machine learning of linear differential equations using gaussian processes[END_REF], quasilinear [START_REF] Sirignano | DGM: A deep learning algorithm for solving partial differential equations[END_REF], and nonlinear partial differential equations [START_REF] Raissi | Hidden physics models: Machine learning of nonlinear partial differential equations[END_REF][START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. All these approaches are based on the use of Neural Networks, which, in general, require big data for training. In machine learning, probabilistic or statistical learning deals with the problem of inferring a function based on data. In particular for supervised learning, it consists in training from a set of data consisting of input-output pairs. The statistical and probabilistic learning approaches have been extensively developed (see for instance, [START_REF] Schölkopf | Kernel principal component analysis[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] Aggarwal | Mining Text Data[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] Balcan | Statistical active learning algorithms[END_REF][START_REF] James | An Introduction to Statistical Learning[END_REF][START_REF] Dong | Knowledge vault: A web-scale approach to probabilistic knowledge fusion[END_REF][START_REF] Ghahramani | Probabilistic machine learning and artificial intelligence[END_REF][START_REF] Taylor | Statistical learning and selective inference[END_REF]) and play an increasingly important role in computational physics and engineering science, in particular for problems using large scale computational models. In this framework, statistical methods for supervised learning have been developed in the form of surrogate models from which approximations of the expensive functions can be easily evaluated [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF][START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Eldred | Design under uncertainty employing stochastic expansion methods[END_REF][START_REF] Yao | Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles[END_REF][START_REF] Homem-De Mello | Monte Carlo sampling-based methods for stochastic optimization[END_REF][START_REF] Papadopoulos | A neural network-based surrogate for carbon nanotubes with geometric nonlinearities[END_REF][START_REF] Giovanis | Data-driven surrogates for high dimensional models using gaussian process regression on the grassmann manifold[END_REF]. Among the many methodologies available on this topic, we are interested here in Probabilistic Learning on Manifolds (PLoM), which was initially introduced in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] and specially developed for solving unsupervised and supervised problems in the case of small data, as opposed to big data. The PLoM is not based on the use of neural networks. The first initial work published in 2016 was completed from a theoretical point of view in [START_REF] Soize | Entropy-based closure for probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF], extended to the case of the polynomial chaos representation [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF], to the case of the learning in presence of physics constraints [START_REF] Soize | Physics-constrained non-gaussian probabilistic learning on manifolds[END_REF], and to the sampling of Bayesian posteriors [START_REF] Soize | Sampling of bayesian posteriors with a non-gaussian probabilistic learning on manifolds from a small dataset[END_REF]. The PLoM has been used with success for probabilistic nonconvex optimization under constraints and uncertainties [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF] that has allowed for analyzing very complex engineering systems such as hypersonic combustion flows [START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF][START_REF] Soize | Enhancing model predictability for a scramjet using probabilistic learning on manifolds[END_REF], optimal placement of wells [START_REF] Ghanem | Optimal well-placement using probabilistic learning[END_REF], the design optimization under uncertainties of mesoscale implants [START_REF] Soize | Design optimization under uncertainties of a mesoscale implant in biological tissues using a probabilistic learning algorithm[END_REF], the quantification of model uncertainties in nonlinear computational dynamics [START_REF] Farhat | Feasible probabilistic learning method for model-form uncertainty quantification in vibration analysis[END_REF][START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF], the fracture paths in random composites [START_REF] Guilleminot | Data-driven enhancement of fracture paths in random composites[END_REF], the digital twin for composite material systems [START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF], the computation of Sobol indices for complex stochastic system [START_REF] Arnst | Computation of sobol indices in global sensitivity analysis from samll data sets by probabilistic learning on manifolds[END_REF], and the detection of structural changes from dynamic monitoring [START_REF] Soize | Machine learning for detecting structural changes from dynamic monitoring using the probabilistic learning on manifolds[END_REF].

(ii) Objective of the paper. As already mentioned, the PLoM is not based on the use of neural networks and was developed for the unsupervised or supervised cases of small data. For the supervised case, these small data constitute the training dataset that is constructed as N d input-output pairs. The PLoM then allows, using only the training dataset, and by a relying on purely probabilistic arguments, to generate a large number n MC N d of additional realizations (the learned dataset) while preserving the concentration of the probability measure on an intrinsic manifold delineated using diffusion maps [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF]. We are then entitled to ask the following fundamental question: do the n MC additional realizations generated by the PLoM satisfy, in a "reasonable" sense, the stochastic Partial Differential Equation (PDE) that was used to build the training dataset. As previously explained, the PLoM has been validated for many situations and the convergence of the learning with respect to N d has theoretically and numerically been studied. This paper will give further confirmation of this point. However, for certain problems, we want to be able to control that the additional realizations satisfy, with a quantifiable metric, the stochastic PDE underlying the stochastic computational model used to generate N d points of the training dataset. This is the objective of this paper: to develop and validate the PLoM method constrained by the stochastic PDE that has been used to build the training dataset. It should be noted that the method requires evaluations of the residual of the PDE. As the proposed method is placed in a general framework, these evaluations cannot be obtained explicitly by an algebraic calculation. They can therefore only be calculated using the computational model associated with the stochastic boundary value problem. Under these conditions, the intrusive character of the numerical evaluations of the residue can arise. However, in the framework of the algorithm proposed in this paper, the commercial and in-house software allow, in general, such residual evaluations to be made in a nonintrusive manner with an almost negligible numerical cost (this is a possibility offered by the software used in computational sciences and engineering ).

(iii) Novelty of the paper. Two novel extensions of the PLoM are presented, which allow to give a constructive answer to our objective. The first is to extend the initial PLoM for nonlinear stochastic dynamical systems (that is to say, physics system depending on time) such as those of the unsteady computational fluid dynamics (CFD) or the nonlinear solid dynamics (SD) in presence of a vector-valued random control parameter. Indeed the PLoM was developed for random vectors and it is therefore necessary to extend it to the case of vector-valued stochastic processes for analyzing these types of stochastic dynamical systems. The second is to extend the PLoM method to take into account the constraint that the generated additional realizations verify the stochastic PDE in a certain sense. These two novel extensions constitute a real challenge and correspond to a difficult problem, which will necessitate the introduction of appropriate methods.

(iv) Organization of the paper. In order to facilitate the reading of this work, we have chosen to present an abstract formulation in the first part in order to better introduce the concepts and methods. These developments are presented for a generic nonlinear stochastic computational dynamical model, which results from the spatial discretization of a nonlinear stochastic boundary value problem. We delay the introduction of the complexity in the stochastic boundary value problem until the sections where the respective applications are introduced. We can thus focus our attention initially on setting up the PLoM with PDE constraints. Later in the paper, we consider nonlinear PDEs with random coefficients. We intersperse in that presentation additional new developments required to efficiently implement our proposed procedure. We specifically introduce two applications that are representative of broad classes of problems, namely a first order (CFD case) and a second order (SD case) problems, which involve one field (SD case) or several fields (CFD case) whose partial time derivatives are different according to the field in question (eg. velocity field and pressure field).

(v) Brief discussion on the hypotheses on which the PLoM method has been built [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF]. The PLoM approach starts from a training dataset made up of a relatively small number N d of points (initial realizations). For the supervised case, it is assumed that the training dataset is related to an underlying stochastic manifold related to a R n x -valued random variable X = (Q, W) with Q = f(U, W) = F(W). The measurable mapping f is unknown and F is also an unknown stochastic mapping. The probability distributions of the vector-valued random variables W (control parameter) and U (non-controlled parameter) are given. The stochastic manifold is defined by the unknown graph {w, F(w)} for w belonging to an admissible set C w . In the PLoM construction, its is not assumed that this stochastic manifold can directly be described; for instance, it is not assumed that there exist properties of local differentiability (moreover, the manifold is stochastic). Under these conditions, the probability measure of X is concentrated in a region of R n x for which the only available information is the cloud of the points of the training dataset. The PLoM method makes it possible to generate the learned dataset whose n MC N d points (additional realizations) are generated by the probability measure that is estimated from the training dataset. The concentration of the probability measure is preserved thanks to the use of the diffusion-maps basis that allows to enrich the available information from the training dataset. It should also be noted that the estimate of this unknown probability measure cannot be performed from the training dataset by using an arbitrarily estimator. It must be parameterized in a manner that permits convergence to any probability measure as its number of points goes towards infinity. The PLoM method therefore does not only consist in generating points that belong to the region in which the measure is concentrated, but also allows these additional points to be realizations of the estimate probability measure with the convergence properties evoked above. The choice of the kernel estimation method for estimating the probability measure from the training dataset guarantees that this required fundamental property is satisfied [START_REF] Soize | Probabilistic learning on manifolds[END_REF]. Considering all of these assumptions in the construction and in the methodology of the PLoM makes PLoM as a general approach for small datasets for arbitrary non-Gaussian measures.

Notations

Lower-case letters such as q or η are deterministic real variables. Boldface lower-case letters such as q or η are deterministic vectors. Upper-case letters such as X or H are real-valued random variables. Boldface upper-case letters such as X or H are vector-valued random variables. Lower-case letters between brackets such as [x] or [η] are deterministic matrices. Boldface upper-case letters between brackets such as [X] or [H] are matrix-valued random variables.

[I n ]: identity matrix in M n . M n,N : set of all the (n × N) real matrices. M N : set of all the square (n × n) real matrices.

M +0

n : set of all the positive symmetric (n × n) real matrices. M + n : set of all the positive-definite symmetric (n × n) real matrices. R: set of all the real numbers. R n : Euclidean vector space on R of dimension n.

[x] k j : entry of matrix [x].

[x] T : transpose of matrix [x].

• : usual Euclidean norm on R n . < x , y >: usual Euclidean inner product in R n . δ αβ : Kronecker's symbol.

In this paper, any Euclidean space E (such as R N ) is equipped with its Borel field B E , which means that (E, B E ) is a measurable space on which a probability measure can be defined. The mathematical expectation is denoted by E.

2.

Abstract nonlinear stochastic computational model, training dataset, and time sampling

Nonlinear stochastic equation with its initial condition

We consider the abstract nonlinear stochastic computational model (SCM) on R N written as

N t Y(t), t, W = 0 a.s. , t ∈]t 0 , T ] , (1)
L IC t 0 Y = y 0 a.s. , t = t 0 , (2)
in which [t 0 , T ] is a given compact interval of R and Y = Y(t), t ∈ [t 0 , T ] is the unknown R N -valued stochastic process. The vector W is the R n w -valued random control parameter of the system, which can correspond either to a set of system parameters, or to the set of generalized coordinates of a reduced-order representations based on the Karhunen-Loève or chaos expansions, (1) of random fields that are the coefficients of partial differential operators of the problem and/or (2) of stochastic processes modeling the time-dependent excitations of the system. In Eq. ( 1), N t is a given nonlinear differential operator in time t of order m d ≥ 1 with values in R N . Equation ( 2) is the deterministic initial condition at t = t 0 . The operator L IC t 0 Y is a linear differential operator with values in R N 0 in which N 0 ≥ N depends on m d . The deterministic vector y 0 is given in

R N 0 . For instance, L IC t 0 Y = Y(t 0 ) and N 0 = N for m d = 1 (first-order nonlinear differential operator N t ) or L IC t 0 Y = Y(t 0 ), Ẏ(t 0 ) with N 0 = 2N for m d = 2
(second-order nonlinear differential operator). In addition, in order to simplify the presentation, we will assume that the first m d derivatives of Y with respect to t at initial time t = t 0 take deterministic values (this simplified hypothesis can easily be removed without difficulty and does not call into question the presented developments). These operators will be detailed for each of the applications presented in the sequel.

Vector-valued random control parameter

The random vector W = (W 1 , . . . , W n w ) is defined on a probability space (Θ, T , P). Its prior probability distribution P W (dw) = p W (w) dw is defined by a given probability density function w → p W (w) on R n w with respect to the Lebesgue measure dw on R n w . Its support, C w ⊂ R n w , is the admissible set of W, which will be specified for each of the applications.

Uniqueness of the stochastic solution of the SCM

Let m d and y 0 be as defined in Section 2.1.

Let C m d ad = y ∈ C m d ([t 0 , T ], R N ) ; L IC t 0 (y) = y 0 in which C m d ([t 0 , T ], R N
) is the set of all the m d -times continuously differentiable functions t → y(t) = (y 1 (t), . . . , y N (t)) defined on [t 0 , T ] with values in R N . For m d = 1, the first time derivative is denoted by ẏ(t) = dy(t)/dt and for m d = 2, the second time derivative is denoted by ÿ(t) = d 2 y(t)/dt 2 , which are continuous functions. Let w = W(θ ) ∈ C w ⊂ R n w be any realization of W for θ ∈ Θ. It is assumed that the problem defined by

N t (y (t), t, w ) = 0 , t ∈]t 0 , T ] , (3)
L IC t 0 (y ) = y 0 , t = t 0 . (4)
admits a unique solution

y = t → y (t) from [t 0 , T ] into R N , which is m d -times differentiable in time, that is to say, (y , w ) ∈ C m d ad × C w .
Consistent with this hypothesis, it is assumed that Eqs. ( 1)-( 2) admit a unique stochastic solution Y sol = Y sol (t), t ∈ [t 0 , T ] that depends on W and y 0 , defined on (Θ, T , P), indexed by [t 0 , T ], with values in R N , which is assumed to be a second-order stochastic process whose trajectories (sample paths) are thus m d -times differentiable in time. Consequently, for all θ in Θ,

Y sol (•; θ ), W(θ ) ∈ C m d ad × C w .
It should be noted that the probability measure of (Y sol , W) (defined for instance by its system of marginal (or cylindrical) distributions) cannot be exactly constructed in the considered general case. We will limit ourselves to building an approximation of this measure using the training dataset and the PLoM without or under constraints.

Training dataset Y d N d

We consider N d independent realizations w ∈ C w of W and the corresponding N d independent realizations y ∈ C m d ad of Y sol (see Section 2.3). The training dataset Y d N d is then defined by

Y d N d = (y 1 , w 1 ), . . . , (y N d , w N d ) ⊂ (C m d ad × C w ) N d . (5)
Thus the construction of the training dataset Y d N d requires solving N d samples of the SCM defined by Eqs. ( 3)-( 4). As we have explained in Section 1, it is assumed that the numerical cost for one computation is very high and consequently, only a small value of N d can be considered (small data case).

Time sampling and time integration scheme

The construction of the training dataset requires the numerical solution of Eqs. ( 3)-( 4). Therefore, we have to introduce a time sampling and a time integration scheme. For all ∈ 1, . . . , N d , the time sampling is performed (thus, is independent from ) with a constant time step ∆t and n time sampling times t 1 , . . . , t n time constituting a finite part J of [t 0 , T ], such that

t n = t 0 + n ∆t , n = 1, . . . n time , J = t 1 , . . . , t n time ⊂ [t 0 , T ] , (6)
in which ∆t and n time verify Tt 0 = n time ∆t. For all in 1, . . . , N d , Eqs. ( 3)-( 4) are solved using a time integration scheme (such as central difference scheme or Newmark scheme, etc), which must be adapted to the type of nonlinear operator N t , that is to say, which depends on each application (see Sections 10,[START_REF] Schölkopf | Kernel principal component analysis[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. We then obtain the time sampling of function t → y (t) ∈ C m d ad ,

y (t), t ∈ J = y (t 1 ), . . . , y (t n time ) , y (t n ) ∈ R N . (7)
Using the time integration scheme, the time derivative up to order m d of y can be computed at all sampling times belonging to J. It should be noted that for each pair (y , w ) in Y d N d , function y is only known at the sampling times of J.

Definition of the stochastic process Y and its second-order moments

Let us consider the couple (Y, W) in which Y = Y(t), t ∈ [t 0 , T ] is a mean-square continuous second-order stochastic process, defined on (Θ, T , P), indexed by [t 0 , T ], with values in R N , whose trajectories are assumed to be m d -times differentiable in time, and where the N d independent realizations of (Y, W) are those defined in Eq. [START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF]. Therefore, (Y, W) is an approximation of (Y sol , W) which depends on

N d . Let t → y(t) = E Y(t) be the mean function of Y, which is a continuous function from [t 0 , T ] into R N and let (t, t ) → [C Y (t, t )] = E Y(t)-y(t) Y(t )-

y(t )

T be its matrix-valued covariance function that is a continuous function from [t 0 , T ] × [t 0 , T ] into M N . These two functions are estimated at the sampling times as follows: for n and n in 1, . . . , n time ,

y(t n ) = 1 N d N d =1 y (t n ) , [C Y (t n , t n )] = 1 N d -1 N d =1 y c (t n )y c (t n ) T (8)
in which, for ∈ 1, . . . , N d , y c is the centered realization such that

y c (t n ) = y (t n ) -y(t n ) , n ∈ 1, . . . , n time . (9)

Reduced-order representation of the nonstationary stochastic process Y

The construction of the reduced-order representation of stochastic process Y = Y(t), t ∈ [t 0 , T ] is based on its classical Karhunen-Loève expansion that will effectively be calculated at sampling times t 0 , t 1 , . . . , t n time .

Eigenvalue problem for the Karhunen-Loève (KL) expansion of stochastic process

{Y(t), t ∈ [t 0 , T ]}
The KL expansion of stochastic process Y(t), t ∈ [t 0 , T ] is based on the use of the spectral analysis of its covariance operator. Let H = L 2 µ 0 ([t 0 , T ], R N ) be the Hilbert space of all the square integrable functions t → y(t) from [t 0 , T ] into R N , equipped with the inner product < y , y > H and its associated norm

y H =< y , y > 1/2 H such that < y , y > H = T t 0 < y(t) , y (t) > µ 0 (dt), in which the measure µ 0 (dt) on R is defined by µ 0 (dt) = (T -t 0 ) -1 dt with dt the Lebesgue measure on R. Since the covariance function (t, t ) → [C Y (t, t )] is continuous on the compact set [t 0 , T ] × [t 0 , T ], it can be seen that T t 0 T t 0 [C Y (t, t )] 2 F µ 0 (dt) µ 0 (dt ) < +∞ in which • F is the Froebenius norm. Therefore, the covariance linear operator C Y such that (C Y y)(t) = T t 0 [C Y (t, t )]
y(t ) µ 0 (dt ) belongs to L HS (H) (space of the Hilbert-Schmidt operators in H) [START_REF] Gelfand | Les Distributions. Tome 4. Application de l'Analyse Harmonique[END_REF]. Consequently, the eigenvalue problem

C Y ϕ = Λ ϕ has a countable spec- trum Λ 1 ≥ Λ 2 ≥ . . . such that α Λ 2 α < +∞.
3.2. Reduced-order representation of discrete stochastic process {Y(t), t ∈ J} Instead of discretizing the continuous eigenvalue problem C Y ϕ = Λ ϕ, we choose to approximate it by the eigenvalue problem related to the covariance matrix of the random vector Y J = Y(t 1 ), . . . , Y(t n time ) with values in R n time N , which corresponds to the time sampling of stochastic process Y(t), t ∈ [t 0 , T ] . This choice is justified by Section 3.1 (the spectrum of C Y is countable) and by the fact that for applications devoted to nonlinear dynamical system, number n time of sampling times is generally very high, and consequently, the constructed approximation is excellent. Consequently, the KL expansion of stochastic process Y is replaced by the Principal Component Analysis (PCA) of random vector Y J .

(i) Eigenvalue problem for the PCA expansion. The eigenvalue problem related to the covariance matrix of random vector Y J is written, using Eq. ( 8), as

1 n time n time n =1 [C Y (t n , t n )] ϕ α (t n ) = Λ α ϕ α (t n ) , n ∈ 1, . . . , n time . ( 10
)
The eigenvalues are ordered as Λ 1 ≥ Λ 2 ≥ . . . ≥ 0 and the family of eigenvectors ϕ α (t 1 ), . . . , ϕ α (t n time ) α in R n time N satisfies the orthogonality property, (1/n time )

n time n=1 < ϕ α (t n ) , ϕ β (t n ) >= δ αβ .
(ii) Reduced-order representation of discrete stochastic process Y(t), t ∈ J . The reduced representation Y (n q ) (t), t ∈ J of Y(t), t ∈ J of order n q is constructed using the eigenvectors associated with the first n q n time N strictly positive eigenvalues, and is written as

Y (n q ) (t n ) = y(t n ) + [V(t n )] Q , n ∈ 1, . . . , n time , (11)
in which y(t n ) ∈ R N is defined by Eq. ( 8) and where [V(t n )] is the matrix defined by

[V(t n )] = [Λ 1/2 1 ϕ 1 (t n ) . . . Λ 1/2 n q ϕ n q (t n )] ∈ M N,n q . ( 12
)
In Eq. ( 11), Q = (Q 1 , . . . , Q n q ) is a second-order centered and orthogonal random variable defined on probability space (Θ, T , P) with values in R n q , which is written as

Q = [Λ] -1 1 n time n time n=1 [V(t n )] T Y(t n ) -y(t n ) , (13)
in which [Λ] is the positive-definite diagonal (n q × n q ) matrix, [Λ] αβ = Λ α δ αβ . The N d independent realizations q in R n q of Q, associated with the N d points of the training dataset Y d N d defined by Eq. ( 5), are computed by

q = [Λ] -1 1 n time n time n=1 [V(t n )] T y (t n ) -y(t n ) . (14)
The empirical estimates of the mean vector q ∈ R n q and the covariance matrix

[C Q ] ∈ M + n q of random vector Q are such that q = (1/N d ) N d =1 q = 0 n q and [C Q ] = 1/(N d -1) N d =1 q (q ) T = [I n q ]
. Using a time integration scheme and the initial condition at t 0 yields the m d time derivatives of y(t) and [V(t)] at t ∈ t 0 , t 1 , . . . , t n time and consequently, allows for calculating the m d time derivatives of Y (n q ) (t n ) for n ∈ 0, 1, . . . , n time . For this type of calculation, an adapted methodology must be developed for each application because N and n time can be very large, inducing possible difficulties due to Random Access Memory (RAM) limitations.

(iii) Calculation of the reduction order n q . Let n q → err KL (n q ) be the error function related to the relative mean-square error between Y(t), t ∈ J and Y (n q ) (t), t ∈ J , defined by err KL (n q ) = 1 -

n q α=1 Λ α (1/n time ) n time n=1 tr [C Y (t n , t n )] . (15)
The value of n q n time N is calculated such that Λ n q > 0 and err KL (n q ) ≤ ε KL where ε KL is a given small positive real number used for error control.

(iv) Algorithm for solving the eigenvalue problem. For ∈ 1, . . . , N d and for α ∈ 1, . . . , n q , we introduce the vectors y J = y c (t 1 ), . . . , y c (t n time ) and ϕ α J = ϕ α (t 1 ), . . . , ϕ α (t n time ) that belong to R n time N , in which y c (t n ) is defined by Eq. ( 9) and ϕ α (t n ) by Eq. [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. To these vectors, we associate the matrices [

y J ] ∈ M n time N,N d and [ϕ J ] ∈ M n time N,n q such that [y J ] = [y 1 J . . . y N d J ] and [ϕ J ] = [ϕ 1 J . . . ϕ n q J ]
. Using these notations, the eigenvalue problem defined by Eq. ( 10) can be rewritten as [C J ] ϕ α J = Λ α ϕ α J in which [C J ] is a symmetric real (n time N ×n time N) matrix whose rank is n q ≤ N d , and which is written as

[C J ] = 1/(N d -1) (n -1/2 time [y J ]) (n -1/2 time [y J ]) T .
In general, for the considered nonlinear SCM, N is large and the number n time of sampling times is also large. Consequently, the dimension (n time N × n time N) of matrix [C J ] of this eigenvalue problem is a full matrix that cannot be assembled (the number of entries is, for a medium problem, of order 10 18 ). In the framework of the present work, N d is small (small size of the training dataset, for instance, of order 100). Therefore, n time N N d , which allows for solving the eigenvalue problem, without assembling [C J ], by using a "thin SVD" [START_REF] Golub | Matrix Computations, Second Edition[END_REF] applied to rectangular matrix n -1/2 time [y J ]. Keeping only the n q largest singular values s α = [s] αα ordered as s 1 ≥ s 2 ≥ . . . ≥ s n q > 0, one obtains the following factorization,

n -1/2 time [y J ] = [b J ] [s] [v] in which [v] ∈ M n q verifies [v] [v] T = [I n q ] and where [b J ] ∈ M n time N,n q is such that [b J ] T [b J ] = [I n q ]. Finally, eigenvalues Λ 1 , . . . , Λ n q and matrix [ϕ J ] ∈ M n time N,n q are given by Λ α = s 2 α /(N d -1) and [ϕ J ] = n 1/2 time [b J ].

Random residual associated with the nonlinear SCM and time subsampling

In this section, we consider any mean-square continuous second-order stochastic process Y(t), t ∈ [t 0 , T ] with values in R N , defined on (Θ, T , P), whose trajectories belong to C m d ad , and which is statistically dependent on a random vector W defined on (Θ, T , P), with values in C w ⊂ R n w . It is assumed that n r independent realizations ( y , w ) of ( Y, W) are known, that is to say, for all ∈ 1, . . . , n r and θ ∈ Θ,

( y , w ) ∈ C m d ad × C w , y = Y(•; θ ) , w = W(θ ) . (16)
As previously, Y(t), t ∈ J and y (t), t ∈ J will denote the time sampling of Y and y . In this section, we will formulate the residual calculation for the couple ( Y, W). We will use it for computing the random residual when ( Y, W) will represent either (Y, W) with n r = N d , W defined in Section 2.2, and Y defined in Section 2.6 (the reference) or will represent (Y (n q ) , W) with n r = n MC , Y (n q ) defined by Eq. [START_REF] Schölkopf | Kernel principal component analysis[END_REF], and for which the probability measure of random vector X = (Q, W) will be estimated either with the PLoM without constraint or with the PLoM under constraints.

Random residual

(i) Definition of the time dependent vector-valued random residual R(t). For all t fixed in [t 0 , T ], the time dependent R N -valued residual R(t) = (R 1 (t), . . . , R N (t)), associated with the couple ( Y, W), is defined by

R(t) = N t Y(t), t, W , t ∈ [t 0 , T ] . (17)
It is assumed that R(t), t ∈ [t 0 , T ] is a mean-square continuous second-order R N -valued stochastic process defined on (Θ, T , P), whose trajectories are almost surely continuous. For ∈ 1, . . . , n r , and using Eqs. ( 16) and ( 17), the realization r (t) of R(t) is written as

r (t) = N t ( y (t), t, w ) . ( 18
)
The R N -valued function t → r (t) is continuous on [t 0 , T ] and is evaluated at the sampling times of J. It should be noted that the approximations of the time derivatives involved in the differential operator for the calculation of the residual are exactly the same as those used for the calculation of each realization of the random solution of the stochastic boundary problem, using the stochastic computational model (see Sections 11 and 12 devoted to Applications 2 and 3).

(ii) Definition of the time dependent real-valued random residual R norm (t). For all t fixed in [t 0 , T ], the time dependent real-valued random residual R norm (t) is the positive-valued random variable is defined by

R norm (t) = 1 √ N R(t) , (19)
in which • is the usual Euclidean norm.

(iii) Definition of the random residual ρ. The random residual associated with the couple ( Y, W) is then defined as follows,

ρ = R norm µ , R norm 2 µ = 1 T -t 0 T t 0 R norm (t) 2 dµ(t) , (20)
in which dµ is a given bounded positive measure on R, which will be defined in Section 4.2 when introducing the time subsampling for computing the residual. From Eqs. [START_REF] Ghahramani | Probabilistic machine learning and artificial intelligence[END_REF] and [START_REF] Taylor | Statistical learning and selective inference[END_REF], it can be seen that

ρ 2 = 1 N(T -t 0 ) T t 0 R(t) 2 dµ(t) . (21)
It should be noted that, if dµ(t) was chosen as dt on bounded interval [t 0 , T ], then for ∈ 1, . . . , n r , the realization ρ of random variable ρ could only be computed using the time sampling J, that is to say, by the formula ( ρ

) 2 = (N n time ) -1 n time n=1 r (t n ) 2 .
(iv) Definition of the L 2 -norm of random residual ρ. Since R is mean-square continuous on compact set [t 0 , T ], it can be deduced [START_REF] Doob | Stochastic processes[END_REF][START_REF] Krée | Mathematics of Random Phenomena[END_REF] that ρ is a real-valued (in fact a positive-valued) second-order random variable. Consequently, ρ belongs to L 2 (Θ, R) and its norm is

ρ L 2 = E ρ 2 . If ( Y, W) is the couple (Y sol , W) defined in Section 2.3, then ρ L 2 = 0. Using the n r independent realizations ρ of ρ, the empirical estimate of ρ L 2 is computed as (1/n r ) n r =1 ( ρ )
Y (n q ) (t), t ∈ [t 0 , T ] of Y(t), t ∈ [t 0 , T ] ,
that is to say, ( Y, W) is the couple (Y (n q ) , W). The positive-valued normalized random residual ρ is then defined by

ρ = ρ ρ ref , (22)
in which ρ is the random residual defined by Eq. [START_REF] Taylor | Statistical learning and selective inference[END_REF]. Consequently, if n r = N d,ref and ( Y, W) is the couple (Y, W) defined in Section 2.6 with N d = N d,ref , then the empirical mean value of ρ is equal to 1 (because, in this case ρ is ρ ref ).

Time subsampling for computing the realizations of the random residual

The role played by measure dµ introduced in Eq. ( 20) is the following. The computation of the n r realizations r (t), t ∈ J of R(t), t ∈ J using Eq. ( 18) requires n r × n time numerical evaluations of the nonlinear operator N t and requires a RAM of n r × n time × N 64-bit words. The values of N, n time , and n r can be large. In general the CPU time is not a limitation because a parallel computation can be done for the loop on the number n r of realizations, but difficulties can occur with possible RAM limitations. On the other hand, such a large value of n time is not necessary for constructing the constraint on ρ and a reasonable approximation of the integral on [t 0 , T ] in the Eq. ( 21) can be used to compute the n r realizations ρ of random residual ρ. For such a computation, we introduce the following subsampling times,

τ n s , n s = 1, . . . , n sp ⊂ t n , n = 1, . . . , n time , n sp n time , (23)
in which τ 1 , . . . , τ n sp are defined as the subsampling times. The measure dµ is then defined by

dµ(t) = (T - t 0 )/n sp n sp n s =1 δ 0 (t -τ n s )
, in which δ 0 is the Dirac measure on R at 0. With this approximation, there is a big gain with respect to the RAM limitations and also a gain for the CPU time. Using this measure, Eq. ( 21) yields

ρ 2 = 1 Nn sp n sp n s =1 R(τ n s ) 2 . ( 24
)
The n r independent realizations ρ of random residual ρ are such that

( ρ ) 2 = 1 Nn sp n sp n s =1 r (τ n s ) 2 , ∈ 1, . . . , n r , (25)
in which r (τ n s ) ∈ R N is computed with Eq. ( 18) for t = τ n s . The time subsampling has to be defined for each application.

Definition of random vectors X and H

The formulation of PLoM under constraints for stochastic processes requires the definition of the initial dataset D d N d of the random vector X = (Q, W) and its PCA constructed using D d N d .

(i) Definition of the initial dataset D d N d . For the probabilistic learning, we introduce the initial dataset

D d N d that is associated with the training dataset Y d N d = (y 1 , w 1 ), . . . , (y N d , w N d
) defined by Eq. ( 5) and with the reduced-order representation defined in Section 3.2-(ii). Each realization y of Y is associated with the realization w of W (see Section 2.3). The reduced-order representation Y (n q ) , defined by Eq. ( 11), allows the N d independent realizations q , = 1, . . . , N d of the R n q -valued random variable Q to be computed by using Eq. ( 14). Let X = (Q, W) be the random variable defined on (Θ, T , P) with values in R n x = R n q × R n w . Each realization q depends on realization w , and x = (q , w ) is a realization of X. The initial dataset D d N d is made up of the N d independent realizations x ,

D d N d = x = (q , w ) ∈ R n x = R n q × R n w , = 1, . . . , N d , (26)
(ii) PCA of random vector X. The first step consists of performing a PCA of X. This step is mandatory to formulate the constraints and to extend the PLoM under constraints. It should be noted that the reduced-order representation of stochastic process Y has already introduced the random vector Q, which is centered and whose covariance matrix is the identity matrix. However, random vector X is not centered and its covariance matrix is not the identity matrix. The purpose of the PCA is to improve the statistical computation of X through de-correlation and normalization. Let x ∈ R n x and [C X ] ∈ M +0 n x be the empirical estimates of the mean vector and the covariance matrix of

X such that x = (1/N d ) N d =1 x and [C X ] = 1/(N d -1) N d =1 (x -x) (x -x) T . Let ξ 1 ≥ ξ 2 ≥ . . . ≥ ξ n x ≥
0 be the eigenvalues and let ψ 1 , . . . , ψ n x be the orthonormal eigenvectors of the eigenvalue problem [C X ] ψ α = ξ α ψ α . For ε PCA > 0 fixed, let ν be the integer such that n q + 1 ≤ ν ≤ n x and xi ν > 0, verifying,

err PCA (ν) = 1 - ν α=1 ξ α tr [C X ] ≤ ε PCA . ( 27
)
The PCA of X yields its representation X (ν) that is written as

X (ν) = x + [ψ] [ξ] 1/2 H , E X -X (ν) 2 ≤ ε PCA E X 2 , ( 28
) in which [ψ] = [ψ 1 . . . ψ ν ] ∈ M n x ,ν such that [ψ] T [ψ] = [I ν ] and where [ξ] is the positive-definite diagonal (ν × ν) matrix such that [ξ] αβ = ξ α δ αβ . The R ν -valued random variable H is obtained by projection, H = [ξ] -1/2 [ψ] T (X -x) , (29)
and its N d independent realizations are such that

η d = [ξ] -1/2 [ψ] T (x -x) ∈ R ν , ∈ 1, . . . , N d . (30)
Using

η 1 d , . . . , η N d d
, the estimates of the mean vector η and the covariance matrix

[C H ] of random vector H are such that η = 0 ν and [C H ] = [I ν ].
(iii) Probability density function of H estimated with the initial dataset. Let P H (dη) = p H (η) dη be the probability distribution of the R ν -valued random variable H defined by Eq. ( 29). The pdf η → p H (η) on R ν with respect to dη is estimated by using the Gaussian kernel-density estimation method (KDE) with the N d independent realizations η d , = 1, . . . , N d defined by Eq. [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. Using the modification proposed in [START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] of the classical formulation [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations[END_REF], p H (η) is estimated, for all η in R ν , as

p H (η) c ν ζ(η) , c ν = 1 ( √ 2π s ν ) ν , (31)
ζ(η) = 1 N d N d =1 exp - 1 2 s 2 ν s ν s ν η d -η 2 , (32)
in which

s ν = 4/ N d (2 + ν) 1/(ν+4)
is the usual multidimensional optimal Silverman bandwidth (taking into account that the covariance matrix of H is [I ν ]), and where

s ν = s ν s 2 ν + (N d -1)/N d -1/2 has been introduced to ensure that E H = ( s ν /s ν ) η = 0 ν and E HH T = s 2 ν [I ν ] + s 2 ν s -2 ν (N d -1)N -1 d [C H ] = [I ν ].
The pdf of H, defined by Eqs. ( 31) and [START_REF] Soize | Probabilistic learning on manifolds[END_REF], is rewritten as

p H (η) = c ν e -φ(η) , φ(η) = -log ζ(η) . (33)
(iv) Remark. By construction of the PCA of X, we have ν ≤ N d . For the applications, ν is such that ν < N d and even can be such that ν N d when the dataset is small. Consequently, there is no overfitting of the KDE with the data. given by Eq. ( 30), the PLoM method [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF] allows for generating n MC N d additional realizations η 1 ar , . . . , η n MC ar of H. Using Eq. ( 28)

allows for computing n MC additional realizations x 1 ar , . . . , x n MC ar of X (ν) such that x ar = x + [ψ] [ξ] 1/2 η ar for all in 1, . . . , n MC . For random vector X (ν) , the learned dataset is defined by D ar n MC = x ar ∈ R n x , = 1, . . . , n MC . From the knowledge of the n MC realizations x ar , we can deduce the n MC additional realizations (q ar , w ar ) in

R n q × R n w of random variable Q (ν) , W (ν) = X (ν) such that (q ar , w ar ) = x ar , ∈ 1, . . . , n MC . (34)
For each additional realization q ar given by Eq. ( 34), and using Eq. ( 11), the corresponding additional realization y ar = y ar (t), t ∈ J of stochastic process Y (n q ) (t), t ∈ J is computed by y ar (t) = y(t) + [V(t)] q ar . The learned dataset for the couple (Y (n q ) , W) is then defined by Y ar n MC = (y 1 ar , w 1 ar ), . . . , (y n MC ar , w n MC ar ) , in which y ar = y ar (t) t ∈ J and where w ar ∈ C w ⊂ R n w . As explained in Section 2.5, the time derivatives up to order m d of y ar are computed at all the sampling times belonging to J by using the time integration scheme. It should be noted that the PLoM method introduces a M ν,N d -valued random variable [H] related to random vector H whose probability distribution is defined by Eq. ( 33). The method allows for generating n MC additional realizations

[η 1 ar ], . . . , [η n MC ar ] in M ν,N d of ran- dom matrix [H].
For each , the additional realization η ar is obtained by a random extraction of one column of matrix [η ar ], such that, for k ∈ 1, . . . , ν and for ∈ 1, . . . , n MC , we have η ar k = [η ar ] k j 0 in which j 1 0 , . . . , j n MC 0 are n MC independent realizations of a uniform discrete random variable J 0 with values in 1, . . . , N d , which are constructed at the beginning of the algorithms and are never modified.

(ii) Computing the random residual of the stochastic equations for the learned realizations. The n MC realizations ρ ar of random residual ρ, corresponding to the learned dataset Y ar n MC , are computed (see Eq. ( 25)), for ∈ 1, . . . , n MC , by ( ρ ar ) 2 = 1/(N n sp ) 18)), for t = τ n s , by r ar (τ n s ) = N t y ar (τ n s ), τ n s , w ar . The corresponding additional realizations ρ 1 ar , . . . , ρ n MC ar of the normalized random residual ρ can then be computed using Eq. ( 22) and written as,

n sp n s =1 r ar (τ n s ) 2 , in which r ar (τ n s ) ∈ R N is computed (see Eq. (
ρ ar = ρ ar / ρ ref , = 1, . . . , n MC . (35)

Definition of the constraints for the PLoM under constraints

In a first formulation of the proposed theory, only the constraint on the normalized random residual ρ of the stochastic equations was introduced in the PLoM under constraints (it will be algo 1 defined in Section 8.3-(i)). Some applications revealed a small difficulty when the constraint on ρ alone was taken into account. As the learning is done on the random vector X = (Q, W), the calculation of the probability measure of X, which is based on the minimization of the Kullback-Leibler cross-entropy does not guarantee the preservation of the prior probability distribution of W. Indeed, we obtained good results, but the prior probability distribution of W is altered during the minimization process. Since the objective is to find the optimal probability measure of the couple (Y (n q ) , W), that is to say, the optimal probability measure of X, which preserves the prior probability distribution of W while minimizing the L 2 -norm of ρ, we have introduced constraints on W in order to best preserve its prior probability distribution (it will be algo 2 or algo 3 defined in Section 8.3-(ii-iii)). In this section, we define the constraints on ρ and W. The PLoM under constraints is formulated using the Kullback-Leibler principle in order to satisfy the constraints on ρ and W. However, as this principle is formulated for the random vector H, the constraints must be expressed in terms of a mathematical expectation of a function of H. For this, we must therefore use the representations constructed previously to express these constraints as a function of H. We introduce the block decomposition of x ∈ R n x and [ψ] ∈ M n x ,ν such that

x = x q x w , [ψ] = [ψ q ] [ψ w ] , ( 36
) in which x q ∈ R n q , x w ∈ R n w , [ψ q ] ∈ M n q ,ν , and [ψ w ] ∈ M n w ,ν . Therefore, X (ν) = (Q (ν) , W (ν)
) defined by Eq. ( 28) can be rewritten as

Q (ν) = x q + [ψ q ] [ξ] 1/2 H , W (ν) = x w + [ψ w ] [ξ] 1/2 H . (37)
Equation ( 11) is ν) . Substituting in it the representation of Q (ν) given by Eq. (37) yields

Y (n q ) (t n ) = y(t n ) + [V(t n )] Q (
Y (n q ) (t n ) = y H (t n ) + [V H (t n )] H , (38)
in which y

H (t n ) = y(t n ) + [V(t n )] x q and [V H (t n )] = [V(t n )] [ψ q ] [ξ] 1/2 .
(i) Constraint associated with the normalized random residual of the stochastic equation. Equation ( 17) is written as

R(t n ) = N t (Y (n q ) (t n ), t n , W (ν) )
. Substituting in its right-hand side the representations Y (n q ) (t n ) given by Eq. ( 38) and W (ν) given by Eq. [START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF] show that R(t n ) depends only on t n and H. Consequently, using Eqs. ( 22) with ( 21),

ρ 2 is a nonlinear mapping η → h c ρ (η) of H such that ρ 2 = h c ρ (H). The constraint consists of imposing a reference value to ρ 2 L 2 = E ρ 2 that is chosen as b c ρ = E ρ 2 ref with ρ ref = ρ ref / ρ ref ,
for which the mathematical expectation is computed using N d = N d,ref realizations of the reduced-order representation Y (n q ) and where N d,ref is the given reference value of N d introduced in Section 4.1-(v). The constraint on ρ is then written as

E h c ρ (H) = b c ρ . (39)
Note that the nonlinear application h c ρ cannot be constructed explicitly by an algebraic calculation. Only, a numerical evaluation of h c ρ (η) for each given η in R ν can be performed using the computational model. In addition we will need to calculate the gradient of h c ρ (η) with respect to η. This important difficulty will be resolved by statistical conditioning (see Section 9).

(ii) Constraints associated with the prior probability model of W (ν) . In order to avoid a significant modification of the prior probability distribution of W (ν) , we impose the mean value and the second-order moment of each component W (ν) j of W (ν) . For all j ∈ 1, . . . , n w , we thus write E W (ν) 39) and ( 40) can be rewritten as

E h c (H) = b c . (41)

Probabilistic learning on manifolds under constraints for stochastic process Y c

For the PLoM under constraints, stochastic process Y (n q ) is denoted by Y c .

Methodology

For the different random variables of the problem, when they will be relative to the PLoM under constraints, we will use an exponent c. Thus H relates to the training dataset and to the PLoM without constraint, while H will be noted H c for the PLoM under constraints. In these conditions, the random variable X (ν) = (Q (ν) , W (ν) ) defined by Eq. ( 37) is denoted by X c = (Q c , W c ) when the constraints are applied. We have to find among all the probability density functions that satisfy the constraint E h c (H c ) = b c , the pdf p H c that is closest to the pdf p H . The use of the Kullback-Leibler minimum cross-entropy principle [START_REF] Bhattacharyya | On the measures of divergence between two statistical populations defined by their probability distributions[END_REF][START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF] yields

p H c = arg min p ∈ C ad,H c R ν p(η) log p(η) p H (η) dη , (42)
in which the admissible set is defined by

C ad,H c = p : R ν → R + , R ν p(η) dη = 1, R ν h c (η) p(η) dη = b c .

Expressing p H c as the solution of an optimization problem

The optimization problem defined by Eq. ( 42) can only be solved by using an adapted numerical method. A major difficulty is due to the normalization constant of p H c that induces considerable difficulties when ν is large, which is bypassed by using the following method. Let λ 0 > 0 be the real Lagrange multiplier that allows for imposing the normalization constraint R ν p(η) dη = 1. Let λ ∈ C ad,λ ⊂ R 1+2n w be the one introduced for imposing the constraint

R ν h c (η) p(η) dη = b c with C ad,λ = λ ∈ R 1+2n w , R ν e -φ(η) -<λ ,h c (η)> dη < +∞ , in
which φ is defined by Eq. [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF]. For λ given in C ad,λ , let H λ be the R ν -valued random variable whose pdf with respect to dη on R ν is written as

p H λ (η) = c 0 (λ) exp -φ(η) -< λ , h c (η) > , (43)
in which c 0 (λ) is the normalization constant that depends on λ. Using the normalization condition of p [START_REF] Agmon | An algorithm for finding the distribution of maximal entropy[END_REF]). The gradient of Γ at point λ is

H λ yields c 0 (λ) = ( R ν exp{-φ(η) -< λ , h c (η) >} dη) -1 . Consequently, (∇c 0 (λ))/c 0 (λ) = -E{h c (H λ )}. Let λ → Γ(λ) be the real function on C ad,λ such that Γ(λ) =< λ , b c > -log c 0 (λ) (see
∇Γ(λ) = b c -E h c (H λ ) , (44)
while its Hessian matrix [Γ (λ)] = [cov h c (H λ ) ] is the covariance matrix of random vector h c (H λ ). It will be assume that [Γ (λ)] belongs to M + 1+2n w for all λ ∈ C ad,λ . This hypothesis is related to the algebraic independence of the constraints and will be checked in the algorithm for each application. Consequently, λ → Γ(λ) is a strictly convex function in C ad,λ . Nevertheless, it cannot be proven that admissible set C ad,λ is convex for the general case because the mapping h c ρ depends on the application that is performed. So, we will assume that the following optimization problem

λ sol = arg min λ∈C ad,λ Γ(λ) , (45)
has a unique solution λ sol that is such that ∇Γ(λ sol ) = 0 1+2n w and that

p H c (η) = p H λ (η) λ=λ sol , p H (η) = p H λ (η) λ=0 1+2nw . (46)
It should be noted that Γ(λ) cannot be evaluated for any given λ because function Γ depends on normalization constant c 0 (λ) that involves an integral on R ν , which cannot be calculated when ν is large.

Iterative algorithm for computing the Lagrange multipliers λ sol

Taking into account Eqs. ( 41) and ( 44), the solution λ sol of the optimization problem defined by Eq. ( 45) is computed by solving the equation ∇Γ(λ) = 0 1+2n w in λ using a Newton iterative method, which allows for constructing a solution H λ sol that satisfies the constraint defined by Eq. ( 41). The Lagrange multiplier λ is rewritten as λ = (λ ρ , λ w ) in which the real-valued Lagrange multiplier λ ρ is associated with the constraint E h c ρ (H) = b c ρ (see Eq. ( 39)) and where the vector-valued multiplier λ w is associated with the constraint E h c w (H) = b c w (see Eq. ( 40)). On this basis, three iteration algorithms, referred as algo 1 to algo 3, have been thoroughly tested through the applications presented in Sections 10, 11, and 12.

-The first algorithm, algo 1, does not take into account the constraints associated with the prior probability distribution of W (ν) . As we have explained in Section 7, this algorithm can be not suitable.

-The two other algorithms, algo 2 and algo 3, take into account not only the constraint on ρ but also the constraints on W (ν) . Algorithm algo 2 is the classical approach for which the constraints vector is globally taken. We have been able to verify that this algo 2 is not always efficient, involving convergence difficulties. This phenomenon is linked to the fact that the two constraints (one on ρ and the other on W (ν) ) play in opposite directions (which does not prevent obtaining the minimum, but with potential numerical difficulties). It is for this reason that we have introduced algorithm algo 3 for which the constraint on ρ and the constraints on W (ν) are decoupled for the computation of the Lagrange multipliers. The numerous numerical tests, which have been carried out, have shown that algo 3 is often the most efficient algorithm.

(i) Algorithm algo 1. The only constraint is on the normalized random residual ρ. The iteration algorithm is written, for i ≥ 1, as

λ i+1 ρ = λ i ρ -Γ ρ (λ i ρ ) -1 Γ ρ (λ i ρ ) with λ 1 ρ = 0, where Γ ρ (λ i ρ ) = b c ρ -E h c ρ (H λ i ρ ) and where Γ ρ (λ i ρ ) = var(h c ρ (H λ i ρ ))
is the variance. At each iteration i, the convergence of the iteration algorithm is controlled by calculating the error err R (i) for R = ρ 2 and the corresponding error err W (i) for W (ν) , such that err

R (i) = | b c ρ -E h c ρ (H λ i ρ ) | / b c ρ , err W (i) = b c w -E h c w (H λ i ρ ) / b c w . (47)
(ii) Algorithm algo 2. The constraints are on the normalized random residual ρ and on the random control parameter W (ν) . The iteration algorithm is written as

λ i+1 = λ i -[Γ (λ i )] -1 Γ (λ i ) , i ≥ 1 , (48)
with

λ 1 = 0 1+2n w , where Γ (λ i ) = b c -E h c (H λ i ) , and where [Γ (λ i )] = [cov h c (H λ i ) ]
is the covariance matrix. At each iteration i, by using λ i = (λ i ρ , λ i w ), the convergence of the iteration algorithm is controlled by calculating the error

err RW (i) = err R (i) 2 + err W (i) 2 , ( 49
)
in which err R (i) and err W (i) are defined by err

R (i) = | b c ρ -E h c ρ (H λ i ) | / b c ρ , err W (i) = b c w -E h c w (H λ i ) / b c w . (50)
(iii) Algorithm algo 3. It is algorithm algo 2 for which λ ρ and λ w are decoupled in Eq. ( 48). The iteration algorithm is then written, for i ≥ 1, as

(λ i ρ , λ i w ) = λ i , λ i+1 ρ = λ i ρ -Γ ρ (λ i ) -1 Γ ρ (λ i ) , λ i+1 w = λ i w -[Γ w (λ i )] -1 Γ w (λ i ) , λ i+1 = (λ i+1 ρ , λ i+1 w ) , (51)
with λ 1 = 0 1+2n w , where (Γ ρ (λ i ), Γ w (λ i )) = Γ (λ i ), and where

Γ ρ (λ i ) = [Γ (λ i )] 11 and [Γ w (λ i )] j j = [Γ (λ i )]
1+ j,1+ j for j and j in 1, . . . , n w . At each iteration i, the convergence of the iteration algorithm is controlled by calculating the errors err R (i) and err W (i) by Eq. ( 50), and err RW (i) by Eq. ( 49).

Additional realizations of H λ i and learned dataset for

(Y c , W c ) (i) Learned dataset Y c,ar n MC for (Y c , W c ).
For each λ i of the iteration algorithm, n MC N d addition realizations of H λ i are generated for estimating E h c (H λ i ) and [cov h c (H λ i ) ]. At convergence for i = i opt , we have λ sol λ i opt . Taking into account Eq. ( 46), we have p H c p H λ for λ = λ i opt . Consequently, at the last iteration, the generator allows for computing n MC additional realizations η 1 c,ar , . . . , η MC c,ar of H c for which the generator is described in Paragraphs (iv) to (vi) below. The learned dataset Y c,ar n MC under constraints is defined by Y c,ar n MC = (y 1 c,ar , w 1 c,ar ), . . . , (y n MC c,ar , w n MC c,ar ) , in which w c,ar and y c,ar = y c,ar (t n ) ∈ R N , n = 1, . . . , n time are given, using Eqs. [START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF] and [START_REF] Soize | Enhancing model predictability for a scramjet using probabilistic learning on manifolds[END_REF], by w c,ar = x w + [ψ w ] [ξ] 1/2 η c,ar and for n ∈ 1, . . . , n time , by y c,ar

(t n ) = y H (t n ) + [V H (t n )] η c,ar .
(ii) Computing the random residual of the stochastic equation for the learned realizations under constraints. The n MC realizations ρ c,ar of random residual ρ, corresponding to the learned dataset Y c,ar n MC under constraints, are computed (see Eq. ( 25)), for ∈ 1, . . . , n MC , by ( ρ c,ar 2 , in which r c,ar (τ n s ) ∈ R N is computed (see Eq. ( 18)), for t = τ n s , by r c,ar (τ n s ) = N t y c,ar (τ n s ), τ n s , w c,ar .

) 2 = 1/(N n sp ) n sp n s =1 r c,ar (τ n s )
(iii) Construction of the diffusion-maps basis. The PLoM [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF] uses the diffusion-maps basis [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF] for projecting the nonlinear ISDE constructed by the PLoM method. This projection allows for preserving the concentration of the probability measure of X c = (Q c , W c ) on the manifold (see the definition of X c at the beginning of Section 8.1). The diffusion-maps basis depends only on η 1 d , . . . , η N d d defined by Eq. ( 30) (initial dataset) and is independent of the iteration algorithm. This basis is represented by the matrix

[g] = [g 1 . . . g m ] ∈ M N d ,m with 1 < m ≤ N d . Let [K] be the symmetric (N × N) real matrix such that [K] = exp -(4 ε diff ) -1 η d -η d 2 ,
b] -1/2 [K] [b] -1/2 γ α = κ α γ α related
to a positive-definite symmetric real matrix with the orthonormalization < γ α , γ β >= δ αβ . Therefore, g α can be deduced from γ α by g α = [b] -1/2 γ α . This construction introduces two hyperparameters: the dimension m ≤ N d and the smoothing parameter ε diff > 0. We refer the reader to [START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Entropy-based closure for probabilistic learning on manifolds[END_REF] for estimating their values. We introduce the matrix [a] defined by

[a] = [g] [g] T [g] -1 ∈ M N d ,m . (52)
(iv) Generator of random vector H λ i by using PLoM algorithm under constraints. Now and throughout the rest of Section 8.4, we remove the exponent i on λ. Let [Z λ (r)], [X λ (r)] , r ≥ 0 be the stochastic process with values in M ν,m × M ν,m , which is the solution of the following ISDE,

d[Z λ (r)] =[X λ (r)] dr , r > 0 , (53)
d[X λ (r)] =[L λ ([Z λ (r)])] dr - 1 2 f 0 [X λ (r)] dr + f 0 d[W wien (r)] [a] , r > 0 , (54)
[Z λ (0)] =[η d ] [a] , [X λ (0)] = [v 0 ] [a] a.s. , r = 0 , (55)
in which:

(1) [a] is the matrix defined by Eq. ( 52),

[η d ] = [η 1 d . . . η N d d ] ∈ M ν,N d , and [v 0 ] ∈ M ν,N d is such that [v 0 ] k k are ν × N d independent
normalized Gaussian real-valued random variables that are fixed and independent of iteration i.

(2) f 0 > 0 is a free parameter allowing the dissipation to be controlled in the stochastic dynamical system and therefore to quickly reach the stationary solution corresponding to the invariant measure. This parameter is chosen such that f 0 ≤ 4/ s ν . The value 4 is generally a good choice. (3) [W wien (r)], r ∈ R + is the stochastic process, defined on (Θ, T , P), indexed by R + , with values in M ν,N d , for which the columns of [W wien ] are N d independent copies of the R ν -valued normalized Wiener process whose matrix-valued autocorrelation function is min(r, r ) [I ν ]. This stochastic process is independent of iteration i (this is the same for all the iterations) and is statistically independent of the random control parameter W.

[L λ ([u])] k = 1 ζ(u ) ∂ζ(u ) ∂u k -< λ , ∂h c (u ) ∂u k > , (56)
in which ζ is defined by Eq. ( 32) and where the gradient of ζ is such that

∇ u ζ(u ) = 1 s 2 ν 1 N d N d j=1 s ν s ν η j d -u exp - 1 2 s 2 ν s ν s ν η j d -u 2 . ( 57
)
For f 0 fixed, let r 0 be a sufficiently large positive real number in order that the stochastic process ([Z λ (r)], [X λ (r)]), r ≥ r 0 is stationary. For given r st such that r st ≥ r 0 , the probability distribution of random variable [Z λ (r st )], [X λ (r st )] is the invariant measure of Eqs. ( 53) and [START_REF] Bhattacharyya | On the measures of divergence between two statistical populations defined by their probability distributions[END_REF]. As proven in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF], the R ν -valued random variable H λ is a column

(see Paragraph (v) below) of random matrix [H λ ] with values in M ν,N d such that [H λ ] = [Z λ (r st )] [g] T . ( 58
) (v)
λ k = [η λ ] k j 0 in which [η λ ] = [z λ ] [g] T ∈ M ν,N
and where j 0 is defined in Section 6-(i). 9. Computational statistics method for estimating the derivative of h c Equation [START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF] shows that we have to compute, for any given u = (u 1 , . . . , u ν ) in R ν , the partial derivative

∂h c (u )/∂u k for k = 1, . . . , ν in which h c (u ) = h c ρ (u ), h c w (u ) ∈ R 1+2n w (see Section 7-(iii)). The first component h c
ρ of h c involves the nonlinear stochastic equation and since ν can be large, the computation of the partial derivatives is a very difficult problem. The construction of a surrogate model of function u → h c ρ (u ) by using a deterministic approach, such as the meshless methods [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Duarte | H-p clouds, an h-p meshless method, Numerical Methods for Partial Differential Equations[END_REF][START_REF] Breitkopf | Explicit form and efficient computation of mls shape functions and their derivatives[END_REF][START_REF] Rassineux | remeshing by local hermite diffuse interpolation[END_REF][START_REF] Zhang | Meshless methods based on collocation with radial basis functions[END_REF], is not adapted to this difficult and numerically expensive problem taking into account the dimension ν of u and the possible complexity of function h c ρ . We thus propose to develop a computational statistics method based on statistical conditioning for the first component h c ρ of h c . For the other components, that is to say for function h c w , there is no difficulty. We then have to construct, for all η in R ν , an estimation of the gradient ∇h c ρ (η) of h c ρ (η) and, for all i in 1, . . . , 2n w , the gradient ∇h c w,i (η) of h c w,i (η).

9.1. Estimation of the gradient of h c ρ (η) For all η in R ν , the following statistical conditioning is used for calculating h c ρ (η),

h c ρ (η) = E ρ 2 | H = η , (59)
in which the conditional mathematical expectation of ρ 2 given H = η is estimated by using the n MC realizations ρ 1 ar , . . . , ρ n MC ar of ρ (see Eq. ( 35)) and the n MC realizations η 1 ar , . . . , η MC ar of random vector H (see Section 6-(i)), computed by the PLoM without constraint. These realizations have been chosen because the PLoM method without constraint already gives a learned dataset of dimension n MC N d , which already verifies well the stochastic equation.

(i) Normalization and gradient of h c ρ . The scaling of H and of its realizations η ar is necessary. Let η ar = (η ar,1 , . . . , η ar,ν ) and σ ar = (σ ar,1 , . . . , σ ar,ν ) be the empirical estimates of the mean value and the standard deviation of H, computed with η 1 ar , . . . , η MC ar . Let H = ( H 1 , . . . , H ν ) and η = ( η 1 , . . . , η ν ) such that, for all η = (η 1 , . . . , η ν ) in R ν and for all k ∈ 1, . . . , ν , H k = (H k -η ar,k )/σ ar,k and η k = (η k -η ar,k )/σ ar,k . Consequently, Eq. ( 59) can be rewritten as

h c ρ (η) = h c ρ ( η) , η k = η ar,k + σ ar,k η k , (60)
h c ρ ( η) = E ρ 2 | H = η , (61)
in which η = η(η) has to be seen as a function of η. The realizations of H are η 1 ar , . . . , η MC ar such that, for = 1, . . . , n MC and k = 1, . . . , ν, we have η ar,k = (η ar,k -η ar,k )/σ ar,k . From Eqs. [START_REF] Duarte | H-p clouds, an h-p meshless method, Numerical Methods for Partial Differential Equations[END_REF], it can be deduced that, for all k ∈ 1, . . . , ν , we have

∂h c ρ (η)/∂η k = σ -1 ar,k ∂ h c ρ ( η)/∂ η k .
(ii) Estimation of h c ρ ( η). Let (r, η) → p R, H (r, η) be the joint pdf on R × R ν of the positive-valued random variable R = ρ 2 and H with respect to the Lebesgue measure on R × R ν . We then have

p H ( η) = R p R, H (r, η) dr . ( 62
)
From Eq. ( 61), it can be deduced that

h c ρ ( η) = 1 p H ( η) R r p R, H (r, η) dr . (63)
Let r 1 ar , . . . r n MC ar be the realizations of R such that r ar = (ρ ar ) 2 , and let σ R be its standard deviation estimated using r 1 ar , . . . r n MC ar . For n MC sufficiently large, the Gaussian kernel-density estimation [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations[END_REF] of p R, H (r, η) can be written as

p R, H (r, η) 1 n MC n MC =1 1 σ R ( √ 2πs SB ) ν+1 exp - 1 2s 2 SB ( r ar -r σ R ) 2 + η ar -η 2 , (64)
in which s SB = 4/ n MC (2 + ν + 1) 1/(4+ν+1) is the usual multidimensional optimal Silverman bandwidth. Using Eqs. ( 62) and ( 64) yields

p H ( η) 1 n MC n MC =1 1 ( √ 2πs SB ) ν exp - 1 2s 2 SB η ar -η 2 . ( 65
)
From Eqs. ( 63) to ( 65), it can be deduced that

h c ρ ( η) = A( η) B( η) , (66)
in which A( η) and B( η) are defined by

A( η) = n MC =1 r ar exp - 1 2s 2 SB η ar -η 2 , B( η) = n MC =1 exp - 1 2s 2 SB η ar -η 2 . ( 67
)
(iii) Estimation of ∇ h c ρ ( η). From Eqs. ( 66) and ( 67), it can be seen that

∇ η h c ρ ( η) = 1 B( η) ∇ η A( η) -h c ρ ( η) ∇ η B( η) , (68)
in which ∇ η A( η) and ∇ η B( η) are given by

∇ η A( η) = 1 s 2 SB n MC =1 r ar ( η ar -η) exp - 1 2s 2 SB η ar -η 2 , ∇ η B( η) = 1 s 2 SB n MC =1 ( η ar -η) exp - 1 2s 2 SB η ar -η 2 . ( 69
)
(iv) Numerical conditioning. Equations ( 66) to [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF] show the presence of exponentials of -o ( η) with o ( η) = 1 2s 2 SB η ar -η 2 . If, for a given η, o ( η) is large enough for all = 1, . . . , n MC , then we can get a value close to 0/0 or c/0 for the computation of ∇ η h c ρ ( η), which is "almost" indeterminate or "almost" in overflow, since expo ( η) can be very small for all . It can be seen that if expo ( η) is replaced by exp -(o ( η)e( η) in Eqs. ( 66) to [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF] in which e( η) is independent of , then h c ρ ( η) defined by Eq. ( 66) and ∇ η h c ρ ( η) defined by Eq. ( 68) are not modified. Therefore, for improving the numerical conditioning, in Eqs. [START_REF] Griebel | Numerical Simulation in Fluid Dynamics: a Practical Introduction[END_REF] to [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF]

, exp -1 2s 2 SB η ar -η 2 is replaced by exp -1 2s 2 SB η ar -η 2 -e( η) in which e( η) = (1/n MC ) n MC =1 η ar -η 2 /(2s 2 SB ).

Estimation of the gradient of h c

w,i (η) for i = 1, . . . , 2n w From Section 7-(ii) and for all j ∈ 1, . . . , n w , we have h c w, j (H) = W (ν) j and h c w, j+n w (H) = (W (ν) j ) 2 . Using Eq. ( 37) yields, for all η = (η 1 , . . . , η ν ) ∈ R ν ,

h c w, j (η) = x w +[ψ w ] [ξ] 1/2 η j , h c w, j+n w (η) = x w +[ψ w ] [ξ] 1/2 η j 2 . ( 70
)
It can easily be deduced that, for all j ∈ 1, . . . , n w and for all k ∈ 1, . . . , ν ,

∂h c w, j (η) ∂η k = [ψ w ] [ξ] 1/2 jk , ∂h c w, j+n w (η) ∂η k = 2 x w +[ψ w ] [ξ] 1/2 η j [ψ w ] [ξ] 1/2 jk . ( 71
)

Application 1: nonlinear stochastic dynamical system

This application deals with the nonstationary stochastic response of a nonlinear dynamic system with one degree of freedom (Duffing oscillator) excited by a nonstationary stochastic process and with a random control parameter. This example was deliberately chosen to be very simple in order to quickly show to the reader the role played by the constraints on the stochastic equation during the learning carried out by the PLoM. For this application, we will not give the algebraic complements necessary for its implementation and we will not do the detailed analysis. These aspects will be presented for Applications 2 and 3.

Ÿ(t) + 2 χ d g 1 (W 1 ) Ẏ(t) + g 1 (W 1 ) 2 1 + k b Y(t) 2 Y(t) = γ(t, W 2 ) , t ∈]0, T ] , Y(0) = Ẏ(0) = 0 a.s. , in which W = (W 1 , W 2 )
is the R 2 -valued random control parameter (n w = 2) defined on the probability space (Θ, T , P) such that W 1 and W 2 are independent normalized Gaussian random variables, where χ d = 0.05, k b = 5 × 10 8 , T = 0.7325 s, and where Y = Y(t), t ∈ [0, T ] is the unknown real-valued stochastic process defined on (Θ, T , P).

(ii) Nonstationary stochastic excitation and random control parameters. The stochastic excitation is defined by

γ(t, W 2 ) = t T g 2 (W 2 ) 2 1 + 0.05 sin(ω b t) exp - 1 1 -( 2t T -1) 4 , ∀ t ≥ 0 , in which ω b = 2π × 100 rad/s. For j = 1, 2, we have g j (W j ) = g j 1 + √ 3 δ j (2 U j -1) in which U j = 1 + erf(W j / √
2) /2, with δ j = 0.2 and where g j = E g j (W j ) with g 1 = ω b and g 2 = 6.

(iii) Values of the main parameters. The time sampling is done with n time = 2 931 sampling times for which the time step is ∆t = 2.5×10 -4 s. There is no time subsampling, that is to say, n sp = n time . The number of points in the training dataset is N d = 80, the KL-expansion of stochastic process Y has been performed with ε KL = 10 -6 yielding n q = 25. The PCA of random vector X of dimension n x = n q + n w = 27 has been performed with ε PCA = 10 -6 yielding ν = 27.

For ε diff = 63, the dimension of the diffusion-maps basis is m = 28. In the ISDE, f 0 = 4 and the parameters for the Störmer-Verlet algorithm (Appendix A) are ∆r = 0.19914, l 0 = 100, and M 0 = 20.

Reference stochastic solution

Figure 1 is related to the reference solution that is constructed using a training dataset with N d = 1 000 points. Figure 1-(a) displays the time evolution of the pdf y → p Y (y, t) of random variable Y(t) for t = 0.1785, 0.2977, 0.3960, and 0.5535 s, which are estimated with N d = 1 000 and 10 000 points in the training dataset. It can be seen that convergence with respect to N d is obtained for N d = 1 000 (the curves for 10 000 are very close to those for 1 000). The convergence of the reference solution constructed with N d = 1 000 can also be viewed in Fig. 1-(

Analysis of the role played by theconstraints on the probability distribution of W

For simplifying the writing, W (ν) is simply noted W. Figure 2 allows for analyzing the role played by the constraints on the probability distribution of W during the minimization of the normalized random residual. The computation is performed with N d = 80 points in the training dataset and n MC = 1 000 additional realizations for the PLoM under constraints. Figure 2-(a) displays the graph of error function i → err W (i) (see Eq. ( 47)) computed with algo 1 of the PLoM under the constraint on ρ but without the constraints on W. At iteration i opt = 10, the relative error err R (i opt ) (see Eq. ( 47)) is minimum while i → err W (i) (see Eq. ( 47)) is increasing from i = 5 and does not reach a minimum for i opt = 10. This phenomenon can also be seen of Figure 2-(b) that displays the pdf w 2 → p W 2 (w 2 ) of W 2 estimated with the training dataset, with algo 1 of the PLoM under the constraint on ρ but without the constraints on W, and with algo 3 of the PLoM under the constraints on ρ and W. It can be seen that algo 3 allows for obtaining a pdf of W 2 , which is close to the pdf that corresponds to the training dataset.

Results of the PLoM without and under constraints for a given dimension of the training dataset

In this section, the computation is performed with N d = 80 points in the training dataset and n MC = 1 000 additional realizations for the PLoM under constraints. In order to visualize the type of sample path of stochastic process Y sol , Fig. 3-(a) displays the graph of 5 realizations y = Y sol (•, θ ) of stochastic process Y sol . Fig. 3-(

b) displays n MC = 1 000 additional realizations of components (W c 1 , W c 2 , Q c 2 ) of random vector X c = (Q c , W c ).
A concentration of these realizations can be observed. Concerning the pdf y → p Y (y, t) of random variable Y(t), Fig. 4 shows, for the two instants t = 0.2977 s and t = 0.5535 s, the reference constructed with N d = 1 000 and, for N d = 80 points and n MC = 1 000 additional realizations, the estimation constructed with the PLoM without constraint and the one with the PLoM under constraints using algo 3 for which the minimum is obtained for i opt = 8. This figure already shows that the PLoM without constraint gives good results but that the PLoM under constraints significantly improves the result. In particular, it can be seen in figures (c) and (d), which are in a logarithmic scale for the ordinates, that the PLoM under constraints makes it possible to correctly predict the tails of the probability density function of the reference.

Application 2: unsteady 2D Navier-Stokes equation for incompressible fluid

Boundary value problem

We consider a dimensionless 2D formulation in a Cartesian coordinate system (O, x 1 , x 2 ), whose generic point is x = (x 1 , x 2 ) . The 2D open domain containing the fluid is Ω =]0, x 1 [×]0, x 2 [ with x 1 = 0.9 and x 2 = 1. The boundary is written as

∂Ω = Γ 0 ∪ Γ u in which Γ u = x 1 ∈ [0, x 1 ], x 2 =
x 2 and where Γ 0 = ∂Ω\Γ u . Let u(x, t) = (u 1 (x, t), u 2 (x, t)) be the dimensionless velocity field and p(x, t) be the dimensionless pressure field, in which t is the dimensionless time. These quantities are related to the non dimensionless quantities (denoted with a tilde) as follows: It is assumed that ν CFD is spatially uniform (independent of x and t) and consequently, the energy equation is not required. The boundary conditions on the non dimensionless boundary

x = x/ ref , t = t V ref / ref in
∂ Ω = Γ 0 ∪ Γ u are written, for all t in [t 0 , T ], as u(•, t) = (0, 0) on Γ 0 and u(•, t) = ( V CFD , 0) on Γ u where V CFD in m/s is written as V CFD = V CFD × V ref in which V CFD
is the given dimensionless velocity that is assumed independent of x and t. The dimensionless boundary value problem for the unsteady incompressible flow in Ω is defined as follows. For all t ∈]t 0 , T ] and for all x ∈ Ω, we have the following dimensionless 2D Navier-Stokes equations for an incompressible flow,

∂u ∂t + u • ∇u + ∇p - 1 Re ∇ 2 u = 0 , (72)
∇ • u = 0 , (73)
This staggered grid allows the boundary conditions to be easily implemented. Let be N u = N u 1 + N u 2 and N = N u + N p . At time t, the spatial discretization of random fields u(•, t) = u 1 (•, t), u 2 (•, t) and p(•, t) are the R N u -and R N p -valued random variables U(t) and P(t), and we define the R N -valued random variable Y(t) by

Y(t) = U(t), P(t) . ( 77
)
(iii) Computational model. The time discretization scheme and the spatial discretization allow for computing the time-spatial-discretization Y(t n ) = U(t n ), P(t n ) , and the time derivative U(t n ), of the solution of the boundary value problem, at time t n , for all n in 1, . . . , n time . In that condition, Eq. ( 1) is rewritten, at sampling time t n , as

N NS U(t n ), U(t n ), P(t n ), W = 0 N u a.s. , (78)
N DIV U(t n ) = 0 N p a.s. , (79)
in which N NS is a nonlinear operator and where N DIV is a linear operator.

(iv) Values of the parameters for the Navier-Stokes solver. The staggered grid is such that n x 1 = 100 and n x 2 = 111.

We then have N u 1 = 10 989, N u 2 = 11 000, N u = 21 989, N p = 11 100, and N = 33 089. The time step has been chosen to ∆t = 0.03 and leads to a convergent solution up to the steady (stationary) stochastic flow with a L 2 -norm of the residual (level of that is acceptable (see Section 11.4). Clearly, this approximation could be improved in decreasing the value of ∆t but there is no point in unnecessarily increasing the computational costs, because that does not change anything on the validation of the PLoM under constraints. (i) For computation of the steady (timestationary) stochastic flow (see Section 11.4-(i)), T = 300 and n time = 10 000. For computing the random residual, n sp = 500. (ii) For computation of the reference unsteady (time-nonstationary) stochastic flow (see Section 11.4-(ii)), t = 30 is chosen for analyzing the performance of the proposed PLoM method under constraints, which is an instant corresponding to the unsteady flow regime. We then have T = 30 and n time = 1 000.

(v) Reduced-order representations (KL expansions) of the nonstationary stochastic processes Y = (U, P) and U. The KL expansion of the nonstationary R N -valued stochastic process Y(t), t ∈ [t 0 , T ] is performed at sampling times t 1 , . . . , t n time , yielding the representation

Y(t n ) = y(t n ) + [V(t n )] Q. Since Y = (U, P), we have U(t n ) = u(t n ) + [V u (t n )] Q , ( 80
) P(t n ) = p(t n ) + [V p (t n )] Q , ( 81
) in which u(t n ) ∈ R N u , [V u (t n )] ∈ M N u ,n q , p(t n ) ∈ R N p , and [V p (t n )] ∈ M N p ,n q . The time derivative of U(t) at t n is calculated by U(t n ) = (U(t n ) -U(t n-1
))/∆t, which allows for deducing u(t n ) and [ Vu (t n )]. The KL expansion of U can then be written as,

U(t n ) = u(t n ) + [ Vu (t n )] Q . ( 82
)
For N d = 100 and T = 30, Figs. 5-(a) and-(b) show the distribution of the eigenvalues α → Λ α of the KL-expansion of stochastic process Y(t) = U(t), P(t) , t ∈ [t 0 , T ] and the error function n q → err KL (n q ) defined by Eq. ( 15). Fixing the tolerance ε KL to the value 10 -6 yields the optimal value n q = 24.

(vi) Calculation of the random residual ρ using the subsampling times. Let R NS u 1 (τ n s ) and R NS u 2 (τ n s ) with values in R N u 1 and R N u 2 be the block decomposition of R NS (τ n s ) associated with the components relative to the x 1 -and x 2 -axes. We then have

R NS (τ n s ) = R NS u 1 (τ n s ), R NS u 2 (τ n s ) with values in R N u = R N u 1 × R N u 2 .
The positive-valued random variable R norm (τ n s ) is defined by Eq. ( 19), but the Euclidean norm is adapted to the structure of random vector R(τ n s ) as

R norm (τ n s ) = R NS u 1 (τ n s ) 2 / N 2 u 1 + R NS u 2 (τ n s ) 2 / N 2 u 2 + R DIV (τ n s ) 2 / N 2 p 1/2 . ( 85
)
The random residual ρ is then calculated by using Eq. ( 24). For t ≥ 300, the stochastic solution is stationary in time. The mean value m P (x, t) and the standard deviation σ P (x, t) of the pressure random field P(x, t), x ∈ Ω, t ≥ 300 are independent of t. For N d = 50 realizations of the training dataset, Fig. 6 shows the estimates of the mean value and standard deviation of the pressure random field at t = 300.

(ii)-Reference unsteady (time-nonstationary) stochastic flow. Instant t = 30 is selected for analyzing the performance of the proposed method in the unsteady flow regime. A reference unsteady solution has been computed with N d = 4 000 realizations from the training dataset requiring a significant CPU effort. This reference solution is used for analyzing all the convergence aspects of the PLoM without and under constraints. Figures 7-(a) and (b) show the estimates of the mean value and standard deviation of the pressure random field at t = 30. The estimate of E ρ is 1.53 × 10 -3 while the estimate of ρ L 2 is 1.64 × 10 -3 . Figure 7-(c) displays the graphs of the pdf r → p ρ N d (r) of the random variable respectively. Figure 7-(d) displays the graph N d → ρ L 2 (N d ). These results show that the convergence with respect to N d is reasonably obtained for N d = 4 000.

ρ N d = ρ/ ρ N d in which ρ N d = E ρ is estimated using N d realizations. Therefore, E ρ N d = 1. For N d = 400,
11.5. Analysis of the role played by the constraints on the probability distribution of W As in Application 1, for simplifying the writing, W (ν) is simply noted W. We consider the training dataset of dimension N d = 100 and n MC = 400 additional realizations (learning) are computed by using the PLoM under constraints. Figures 8-(a) and (b) correspond to algo 1 of the PLoM under the constraint on ρ but without the constraints on W. The mean value ρ ref = E ρ ref is estimated using the N d,ref = 100 realizations. For N d = 100, n MC = 400, t = 30, and using algo 1, Figs. 8-(a) and (b) show the error functions i → err R (i) and i → err W (i) defined by Eq. ( 47). As expected, Fig. 8-(a) shows the decrease in the error function of ρ, which shows that the optimization process decreases the norm of ρ. On the other hand, since there is no constraint on W, the optimization process to decrease the random residual of the stochastic PDE, modifies the prior probability distribution of W, which is seen by the growth of function i → err W (i) that is plotted in Fig. 8-(b). Figures 8-(c) and (d) quantify this modification of the prior probability density function of each component of random vector W = (W 1 , W 2 ). Figure 8-(c) shows the pdf w → p W 1 (w) of W 1 while Fig. 8-(d) the pdf w → p W 2 (w) of W 2 . On each of these two figures, it is plotted the estimated pdf, on the one hand with the N d = 100 realizations (prior pdf estimated with the training), on the other hand with the n MC = 400 additional realizations (learning) calculated using algo 1 (converged for i opt = 12) and using algo 3 (converged for i opt = 5) of the PLoM under the constraints on ρ and W. It can be seen notable improvement of the probability density function of W = (W 1 , W 2 ) when the constraints on W are taken into account. 11.6. Analysis of the behavior of the iterative algorithms as a function of n MC and convergence criteria The three algorithms, algo 1 to algo 3 have deeply been tested. For the reasons given in Section 8.3, we have used and we present the results obtained with algo 3 for the iterative algorithm. For unsteady stochastic flow at t = 30 and for N d = 100, Fig. 9 shows the evolution of the error functions i → err R (i) , i → err W (i), and i → err RW (i) defined in algo 3, for n MC = 400, 800, and 1 600. For n MC = 400, a minimum of err RW is reached for i opt = 5 with err RW (i opt ) = 2.92, while, for n MC = 1 600, a minimum of err RW is reached for i opt = 7 with err RW (i opt ) = 2.99.

Although the two values are close, one might think that the solution obtained for n MC = 400 is better than that obtained for n MC = 1 600. In fact it is not the case because the larger n MC the better the convergence of the statistical estimators of the solution. The effects of the small difference between 2.92 and 2.99 are negligible.

Results of the PLoM without and under constraints for a given dimension of the training dataset

We consider the training dataset with N d = 100. A large number n MC N d of additional realizations is generated with the PLoM under the constraints on ρ and W, by using algorithm algo 3. The unsteady stochastic flow is observed at t = 30. The challenge is therefore, for this small value N d of points in the training dataset, to show that, for a large value of n MC (with convergence with respect to n MC ), a good approximation of the reference is obtained whereas the use of the PLoM without the constraints gives an acceptable approximation, which can therefore be improved. For this application, the PLoM method under constraints of the stochastic PDE allows an improvement of the PLoM without constraint [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF].

(i) Figure 10 shows the probability density function r → p ρ (r) of ρ as a function of the number n MC of additional realizations generated with the PLoM without constraint and under constraints. Figure 10-(c) shows, for n MC = 4 000, the comparison of this pdf estimated with the PLoM without constraint and under constraints. The analysis of these figures shows a good convergence of the pdf of ρ with respect to n MC . (ii) For the unsteady stochastic flow at t = 30, for N d = 100 and n MC = 4 000, Fig. 11 shows the evolution of the error functions i → err R (i), i → err W (i), and i → err RW (i) computed with algo 3. The iterative algorithm is stopped at iteration 14 for which there is a local minimum that is err RW (14) = 3.65 and for which err R (14) = 5.16 and err W (14) = 0.0015. Figure 12 shows the prior probability density function w → p W 1 (w) and w → p W 2 (w) of To illustrate the distribution of the realizations generated by the PLoM under constraints, Fig. 13 displays the 4 000 additional realizations of components of X c = (Q c , W c ). A concentration of these realizations can be observed.

Figure 14 displays the 2D-plot of x → 1 000 × m P (x, t) and x → 1 000 × σ P (x, t) estimated with n MC = 4 000 additional realizations generated by the PLoM without constraint and are compared to the reference (see Fig. 7). These results show that the standard PLoM method (that is to say without constraint) is already good. It allows for generating a large number n MC of additional realizations (learning) from a small number N d of points of the training dataset. These realizations allow for estimating the second-order quantities of the random field, which are good for the mean and less good for the standard deviation (see the zoom and the associated comment of Fig. 15-(c)). Similarly, Fig. 15 displays the 2D-plot of x → 1 000×m P (x, t) and x → 1 000×σ P (x, t) estimated with n MC = 4 000 additional realizations computed with the PLoM under constraints and shows the comparison with the reference (see Fig. 7). These results show that the PLoM under constraints give an excellent prediction of the second-order quantities, for both the mean and the standard deviation of the random fields. Figure 15-(c) corresponds to a zoom of Fig. 14-(b). At the center of the vortex, it can be seen that the value of the innermost level line (in the center of the vortex) is 0.95772 for the reference solution (solid line) while that which corresponds for the PLoM without constraint, is not the proximal dashed line, but the second deviating from the center. The value of the proximal dashed line is 1.1173 instead of 0.95772. In contrast, Fig. 15-(d) that is a zoom of Fig. 15-(b) shows an excellent coincidence of the level lines for the same region of the fluid domain.

Discussion

The results obtained for this application devoted to unsteady nonlinear stochastic computational fluid dynamics show that the PLoM without constraint is already giving good results (as it was expected, taking into account the previous works). However the novel proposed extension, consisting in constraining the PLoM method by the stochastic PDE allows to improve the results, which are then quite excellent (see zoom in Fig. 15-(d)).

12. Application 3: nonlinear 3D solid dynamics 12.1. Description of the nonlinear 3D dynamical system (i) System definition and geometry. It is a 3D MEMS accelerometer whose scheme is shown in Fig. 1-(A) of Ref. [START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF]. The mobile part is constituted of a deformable rectangular frame in silicon with one vertical deformable central beam attached to it. The suspended part is made up of a deformable parallelepipedic solid with two attached vertical deformable beams that are around the central beam, and is attached to the mobile part by a 3D suspension constituted of 20 linear springs. The MEMS is represented by an 3D-domain Ω ⊂ R 3 and is described in a Cartesian coordinate system Ox 1 x 2 x 3 that is attached to its mobile part and adopted as the reference coordinate system. The origin O is located at the bottom left corner of the device. The axis Ox 1 is horizontal and oriented positively from left to right. The axis Ox 2 is vertical and oriented positively from bottom to top. The axis Ox 3 is perpendicular to the plane Ox 1 x 2 . The overall dimensions in microns ( ×10 -6 m ) of the rectangular frame are 30 × 31 × 4.

(ii) Materials. The suspended and mobile parts are made of a homogeneous, orthotropic, linear elastic material, which is a standard (100) silicon wafer [START_REF] Hopcroft | What is the young's modulus of silicon?[END_REF], whose Young moduli are (86)

12.3. Computational model, KL expansion, random residual, and PLoM parameters (i) Discrete nonlinear dynamical equation. The semidiscretization of the boundary value problem of the nonlinear 3D dynamical system is performed using the Finite Element (FE) method. The FE mesh is shown in Fig. 1-(B) of Ref. [START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF]. There are 7 328 eight-nodes solid elements, 10 675 nodes, and N = 31 410 free dofs. In the reference coordinate system, the semidiscretized stochastic dynamical system is written as

[M] Ÿ(t) + [D] Ẏ(t) + [K(W)] Y(t) + fNL(Y(t)) = γ(t) f M , t ∈]t 0 , T ] , (87)
Y(t 0 ) = 0 N , Ẏ(t 0 ) = 0 N . (88)
In Eq. (87), Y(t) = (Y 1 (t), . . . , Y N (t)) is the vector of the N free degrees of freedom of the displacements in the reference coordinate system. The corresponding velocity and acceleration vectors are Ẏ(t) and Ÿ(t). The matrices [M],

[D], and [K(W)] are the FE mass, damping, and stiffness matrices with values in M + N . The stiffness matrix depends on W via the use of Eq. ( 86) and fNL(Y(t)) is the R N -vector of the nonlinear internal forces at time t induced by the nonlinear material. Finally, f M is the R N vector independent of time t, which represents the spatial part of the external forces in the reference coordinate system, due to the imposed acceleration of boundary Γ 0 . Note that, for the dynamical system under consideration, f M is independent of the damping and stiffness matrices. The time-dependent load vector is such that f M = 9.64 × 10 -12 , max t | γ(t) | = 2 × 10 7 , and consequently, max t γ(t)f M = 5.59 × 10 -6 .

(ii) Damping modeling. Let ϕ 1 , . . . , ϕ n mode be the first n mode = 40 elastic modes, normalized with respect to the mass matrix, associated with the first n mode lowest eigenfrequencies 0 < ω 1 ≤ . . . ≤ ω n mode such that [K(w)] ϕ β = ω (iii) Time sampling and discretization scheme. Let ∆t = 6.9444 × 10 -9 s be the time step. The sampling times are t 1 , . . . , t n time such that t n = t 0 + n ∆t for n = 1, . . . n time in which n time = 6 466 and Tt 0 = n time ∆t = 4.4903 × 10 -5 s. Using Eqs. (87) and (88) yields, because f NL (0

N ) = 0 N , Y(t 0 ) = 0 N , Ẏ(t 0 ) = 0 N , Ÿ(t 0 ) = ÿ(t 0 ) , ÿ(t 0 ) = γ(t 0 ) [M] -1 f M . (89)
The centered Newmark integration method (implicit scheme) is used. At each sampling time, the nonlinear algebraic equation is solved by using a fixed point method with a relative precision of 10 -6 . In order to guarantee the convergence of the fixed point method, a local adaptive time step is used. Knowing Y(t n-1 ), Ẏ(t n-1 ), and Ÿ(t n-1 ) at sampling time t n-1 , and knowing Y(t n ), the following equations are used for calculating Ÿ(t n ) and Ẏ(t n ),

Ÿ(t n ) = a 0 (Y(t n ) -Y(t n-1 )) -a 1 Ẏ(t n-1 ) -a 2 Ÿ(t n-1 ) , (90)
Ẏ(t n ) = Ẏ(t n-1 ) + a 3 ( Ÿ(t n-1 ) + Ÿ(t n )) , (91)
in which a 0 = 4/∆t 2 , a 1 = 4/∆t, a 2 = 1, and a 3 = ∆t/2.

(iv) Computational model. The FE discretization and time discretization scheme previously defined allow for computing the R N -valued time series (Y(t n ), Ẏ(t n ), Ÿ(t n )), n = 1, . . . , n time . In that condition, Eq. ( 1) is rewritten, at sampling time t n , as

N SD ( Ÿ(t n ), Ẏ(t n ), Y(t n ), t n , W) = 0 N a.s. , (92)
in which N SD is a nonlinear operator with values in R N .

(v) KL expansion of the nonstationary stochastic processes Y. The KL expansion of the nonstationary R N -valued stochastic process Y(t), t ∈ [t 0 , T ] is performed at the sampling times, t 1 , . . . , t n time , yielding the representation

Y(t n ) = y(t n ) + [V(t n )] Q. Knowing y(t n ) and [V(t n )]
for n ∈ 1, . . . , n time , their first and second time derivatives are computed by the following recurrences based on the centered Newmark integration method used in Eqs. ( 90) and ( 91).

(a) Initialization of the recurrences. Y(t 0 ) = Ẏ(t 0 ) = 0 N and Ÿ(t 0 ) = ÿ(t 0 ) in which ÿ(t 0 ) is defined in Eq. ( 89). Choosing y(t 0 ) = ẏ(t 0 ) = 0 N and ÿ(t 0 ) = ÿ(t 0 ), it can then be deduced that [V(t

0 )] = [ V(t 0 )] = [ V(t 0 )] = [0 N,n q ].
(b) Recurrences. For computing the mean function t → y(t), the recurrence is written, for n ∈ 1, . . . , n time , as ÿ

(t n ) = a 0 (y(t n ) -y(t n-1 )) -a 1 ẏ(t n-1 ) -a 2 ÿ(t n-1 ) with ẏ(t n ) = ẏ(t n-1 ) + a 3 (ÿ(t n-1 ) + ÿ(t n )), while for the function t → [V(t)], the recurrence is [ V(t n )] = a 0 ([V(t n )] -[V(t n-1 )]) -a 1 [ V(t n-1 )] -a 2 [ V(t n-1 )] with [ V(t n )] = [ V(t n-1 )] + a 3 ([ V(t n-1 )] + [ V(t n )]
). Consequently, the KL expansions of stochastic processes Y, Ẏ, and Ÿ at the sampling times t 1 , . . . , t n time are

Y(t n ) = y(t n ) + [V(t n )] Q , Ẏ(t n ) = ẏ(t n ) + [ V(t n )] Q , Ÿ(t n ) = ÿ(t n ) + [ V(t n )] Q . ( 93
)
For N d = 120, Figs. 16-(a) and -(b) show the distribution of the eigenvalues α → Λ α of the KL-expansion of stochastic process Y(t), t ∈ [t 0 , T ] and the error function n q → err KL (n q ) defined by Eq. ( 15). Fixing the tolerance ε KL to the value 10 -6 yields the optimal value n q = 107. b(t n ) . For calculating the random residual with the KL expansions, the family Y(τ n s ), Ẏ(τ n s ), Ÿ(τ n s ), n s = 1, . . . , n sp is extracted from the family {Y(t n ), Ẏ(t n ), Ÿ(t n ), n = 1, . . . , n time } computed with Eq. (93). Taking into account Eq. (92) at subsampling time τ n s , the R N -valued random residual R(τ n s ), defined by Eq. ( 17), is rewritten as

R(τ n s ) = N SD ( Ÿ(τ n s ), Ẏ(τ n s ), Y(τ n s ), τ n s , W) a.s. (94)
The positive-valued random variable R norm (τ n s ) is defined by Eq. ( 19) and writes R norm (τ n s ) = N -1/2 R(τ n s ) . The random residual ρ is then calculated using Eq. ( 24) with n sp = 60. The mean value ρ ref is calculated for the fixed value N d,ref = 120. For N d = 120 and using the KL expansion, the mean value of ρ is 5.48 × 10 -7 and ρ L 2 = 6.93 × 10 -7 .

(vii) Values of the parameters for the PLoM without and under constraints. We consider always the case N d = 120 for which the KL expansion of stochastic process Y has been performed. The dimension of X = (Q, W) is n x = n q + n w = 107 + 24 = 131. Concerning the PCA of X, Fig. 16-(c) shows the error function ν → err PCA (ν) defined by Eq. ( 27).

Fixing the tolerance ε PCA to the value 10 -3 yields the optimal value ν = 112. For the PLoM without constraint or under constraints, the diffusion-maps basis is computed with ε diff = 240, which yields the optimal dimension m = 113. Concerning the Störmer-Verlet algorithm for solving the nonlinear ISDE: f 0 = 4 and ∆r = 0.2147.

12.4. Stochastic solution of reference and predictions by the PLoM without constraint (i) Solution of reference. A reference solution has been computed using N d = 1 000 realizations from the training dataset requiring a significant CPU effort. Figure 17 shows, for each one of the two observations, two types of graphs.

The first type corresponds to the deterministic function ν freq → dB lin, j (ν freq ) computed by using the deterministic linear computational model for which the control parameters are fixed to w and the nonlinear material between the beams are removed. The second type corresponds, for the probability level p c = 0.98, to the confidence region of the random function ν freq → dB j (ν freq ) computed using with the nonlinear stochastic computational model estimated using the N d = 1 000 realizations. It can be observed that the effect of the structural nonlinearity is very important outside the frequency band [ω min , ω e ] of the excitation. which ρ N d = E ρ is estimated using N d realizations. We then have E ρ N d = 1 and the standard deviations are, for N d = 30, 60, 120, 440, and 1 000, respectively, 0.080, 0.094, 0.091, 0.087, and 0.087. Figure 18-(b) shows the graph N d → ρ L 2 (N d ). These results show that the convergence with respect to N d is reasonably obtained for N d = 1 000, that is the value selected for the reference.

(ii) Predictions using the PLoM without constraint. For N d = 120 and using the nonlinear stochastic computational model, for each one of the two observations and for the probability level p c = 0.98, Fig. 19 shows the confidence region of the random function ν freq → dB j (ν freq ) estimated using n MC = 1 000 additional realizations generated with the PLoM without constraint. These confidence regions are almost identical to those of the reference stochastic solution shown in Fig. 17 (see also the discussion given in Section 12.7).

Analysis of the role played by constraints on the probability distribution of random control parameter W

As previously, for simplifying the writing, W (ν) is simply noted W. For N d = 120, n MC = 1 000, and using algo 1, Fig. 20-(a) shows the decrease in the error function i → err R (i) while Fig. 20-(b) shows the growth of function i → err W (i) for i > 24, proving that the constraint on W must be taken into account. (i) Figure 21 shows the pdf r → p ρ (r) of ρ as a function of the number n MC of additional realizations generated with the PLoM without any constraint and under the constraints. The analysis of these figures shows a good convergence of the pdf of ρ with respect to n MC .

(ii) The case N d = 120 is analyzed with n MC = 1 000 additional realizations generated with the PLoM under constraints. Figure 22 shows the evolution of the error functions i → err R (i), i → err W (i), and i → err RW (i). There is a local minimum for i opt = 72 such that err RW (72) = 0.0032, for which err R (72) = 8.9 × 10 -6 and err W (72) = 0.0023. one hand with the N d = 120 realizations from the training dataset, on the other hand with the n MC = 1 000 additional realizations (learning) estimated with the PLoM without constraint and with the PLoM under the constraints on ρ and W by using algo 3. The figures show that the estimation performed with the PLoM without constraint is already good. Indeed, Fig. 22-(b) shows that the error err W (1) = 0.013 at iteration i = 1, which corresponds to PLoM without constraint, is already small. So, for this case, the introduction of the constraint on W is not decisive, although it allows for improving the little that there is to improve. The mean value of ρ is 6.35 × 10 -7 , its standard deviation is 2.74 × 10 -7 , and ρ L 2 = 6.95 × 10 -7 . In order to illustrate the distribution of the additional realizations generated by the PLoM under constraints, Fig. 24 displays the n MC = 1 000 additional realizations of components of random vector X c = (Q c , W c ). The probability distribution of (W c 1 , W c 2 , Q c 1 ) is concentrated while the one for (W c 4 , W c 10 , Q c 1 ) is not. For N d = 120 and using the nonlinear stochastic computational model, for each one of the two observations, and for the probability level p c = 0.98, Fig. 25 shows the confidence region of the random function ν freq → dB j (ν freq ) estimated using n MC = 1 000 additional realizations generated with the PLoM under constraints. These confidence regions are almost identical to those of the reference stochastic solution shown in Fig. 17.

Discussion

Similarly to the CFD application, the results obtained for this application devoted to the nonlinear computational stochastic dynamics of elastic solids show that the PLoM method without constraints gives good results, which was expected for but not at this level of quality. The use of the proposed extension, consisting in constraining the PLoM method by the stochastic PDE allows to verify the convergence and the robustness of the algorithms. The examination of the confidence domains of the stochastic responses in the frequency domain does not show any significant differences (compare Fig. 19 with Fig. 25). This application therefore also provides an additional validation of the PLoM method without constraint.

Conclusion and discussion

A novel extension of the PLoM method has been It makes it possible to generate a large number of time-dependent additional realizations that verify the stochastic PDE used to generate the training dataset. This is paramount for the case where only a small number of solutions can be secured from the numerical model due to computational cost. The method is general and can be applied to linear and nonlinear time-dependent stochastic boundary value problems. Its implementation requires specialized algorithms related to the Karhunen-Loève expansions of high-dimension nonstationary stochastic processes, in order to circumvent potential difficulties associated with RAM limitations, thus providing an efficient alternative to I/O (Input-Output) on external memory. Taking into account the scalar constraint related to the random residual of the stochastic PDE, a novel statistical formulation and its implementation have been developed for estimating the gradient of the potential that appears in the drift vector of the reduced nonlinear Itô stochastic differential equation, which is the MCMC generator. The prior probability model of the vector-valued random control parameter of the stochastic boundary value problem is given and must not be modified by the PLoM under the constraint related to the random residual. In order to guarantee this property, we have had to introduce a vector constraint on the random control parameters. The results obtained through the presented applications in the field of nonlinear computational fluid and solid stochastic dynamics, allow the following conclusions to be drawn.

-The PLoM constrained by PDE provides further validation of the PLoM method without constraint, proposed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] for which a complete mathematical analysis has been proposed in [START_REF] Soize | Probabilistic learning on manifolds[END_REF]. Although the obtained results concerning the PLoM without constraint were expected, they are particularly good and interesting in terms of statistical methods of learning with small data.

-The PLoM constrained by stochastic PDE allows a direct control of the random residual of the stochastic PDE and thus one is assured that the probability measure constructed by PLoM under constraint, checks better the stochastic boundary value problem that is considered. This novel approach makes it possible to improve learning from small data when the random residual of the stochastic PDE computed with the PLoM without constraint is not sufficiently small. We observe this for the CFD application for which the PLoM under constraints gives excellent results and improve the PLoM without constraint. It should also be noted that the robustness of the PLoM algorithms under constraints allows for maintaining the quality of the results given by the PLoM without constraint when these results are already good, as can be seen in application 3 devoted to nonlinear solid dynamics.

6 .

 6 Probabilistic learning on manifolds (PLoM) without constraint (i) Generation of the learned datasets D ar n MC and Y ar n MC . Using only the realizations η 1 d , . . . , η N d d

j

  = b c w, j and E (W (ν) j ) 2 = b c w, j+n w in which b c w, j and b c w, j+n w are the mean value and the second-order moment of W j (and not of W (ν) j ), which are computed with the training dataset for N d = N d,ref . The use of the representation W (ν) = x w + [ψ w ] [ξ] 1/2 H allows for defining a nonlinear mapping η → h c w (η) from R ν into R 2n w such that E h c w, j (H) = b c w, j and E h c w, j+n w (H) = b c w, j+n w . Introducing the vector b c w = (b c w,1 , . . . , b c w,n w , b c w,1+n w , . . . , b c w,2n w ) ∈ R 2n w , the constraint on W (ν) is written as E h c w (H) = b c w . (40) (iii) Global constraint imposed to random vector H. Defining the vector b c = (b c ρ , b c w ) ∈ R 1+2n w and the mapping h c = (h c ρ , h c w ) with values in R 1+2n w , Eqs. (

  which depends on a real smoothing parameter ε diff > 0. Let [P] be the transition matrix in M N d of a Markov chain such that[P] = [b] -1 [K] in which [b] is the positive-definite diagonal real matrix such that [b] i = δ i N d =1 [K] .For m fixed in 1, . . . , N d , let g 1 , . . . , g m be the right eigenvectors in R N d of matrix [P] such that [P] g α = κ α g α in which the eigenvalues are sorted such that 1 = κ 1 > κ 2 > . . . > κ m . It can easily be proven that the eigenvalues are positive and the largest is κ 1 = 1 for which all the components of the corresponding eigenvector g 1 are equal. Considering the normalization [g] T [b] [g] = [I m ], the right-eigenvalue problem of the nonsymmetric matrix [P] can then be done solving the eigenvalue problem [

( 4 )

 4 the random matrix [L λ ([Z λ (r)])] with values in M ν,m is written as [L λ ([Z λ (r)])] = [L λ ([Z λ (r)] [g] T )] [a],in which, for k = 1, . . . , ν, for = 1, . . . , N d , and for [u] = [u 1 . . . u N d ] ∈ M ν,N d with u = (u 1 , . . . , u ν ), we have

10. 1 .

 1 Nonlinear stochastic equation (i) Nonlinear stochastic dynamics equation with initial conditions. We consider the Duffing oscillator with N = 1,

  b) that shows the graph of function N d → ρ L 2 (N d ) of the random residual ρ. The values of ρ L 2 (N d ) are close for N d = 1 000 and N d = 10 000.

Figure 1 :

 1 Figure 1: (a) time evolution of the pdf y → p Y (y, t) for t = 0.1785 s (green), t = 0.2977 s (black), t = 0.3960 s (blue), and t = 0.5535 s (red) estimated with N d = 1 000 points (thick line) and with N d = 10 000 points (thin line) in the training dataset; (b) graph of N d → ρ L 2 (N d ).

Figure 2 :

 2 Figure 2: For N d = 80 and n MC = 1 000: (a) relative error function i → err W (i) computed with algo 1 of the PLoM under the constraint on ρ but without the constraints on W; (b) pdf w 2 → p W 2 (w 2 ) of W 2 estimated with the training dataset (dashed line), with algo 1 of the PLoM under the constraint on ρ but without the constraints on W (thin solid line), and with algo 3 of the PLoM under the constraints on ρ and W (thick solid line).

Figure 3 :Figure 4 :

 34 Figure 3: (a) Example of 5 realizations y of Y sol . (b) 3D-plot of realizations (w c ar,1 , w c ar,2 , q c ar,2 ) of (W c 1 , W c 2 , Q c 2 ) computed with the PLoM under the constraints: n MC = 1 000 additional realizations (red symbols) and N d = 80 points from the training dataset (blue symbols).

  which ref is a given reference length and where V ref is the reference velocity in m/s (independent of x and t), u = u/ V ref , and p = p/( ρ ref V 2 ref ) with ρ ref the mass density in kg/m 3 (independent of x and t). Let Re = V ref ref /ν CFD be the Reynolds number in which ν CFD is the kinematic viscosity such that ν CFD = η ref / ρ ref where η ref is the dynamic viscosity in Pa × s.

Figure 5 :

 5 Figure 5: For N d = 100, T = 30, and n time = 1 000, KL expansion of nonstationary stochastic process Y: (a) distribution of the eigenvalues α → Λ α and (b) error function n q → err KL (n q ). PCA of random vector X: (c) error function ν → err PCA (ν).

  The mean value ρ ref of the reference random residual ρ ref is calculated for the fixed value N d,ref = 100 of the number N d of training points in the training dataset.(vii) Values of the parameters for the PLoM without and under constraints. We consider the case N d = 100 and T = 30 for which the KL expansions of stochastic processes U, U, and P have been performed. Random vector X = (Q, W) has dimension n x = n q + n w = 26. Concerning the PCA of X, Figure5-(c) shows the error function ν → err PCA (ν) defined by Eq.[START_REF] Homem-De Mello | Monte Carlo sampling-based methods for stochastic optimization[END_REF]. Fixing the tolerance ε PCA to the value 10 -6 yields the optimal value ν = 26. For the PLoM without constraint or under constraints, the diffusion-maps basis is computed with ε diff = 95, which yields the optimal dimension m = 27. Concerning the Störmer-Verlet algorithm for solving the nonlinear ISDE: f 0 = 4 and ∆r = 0.1974. The convergence analysis with respect to n MC is performed by considering the following values 100, 400, 800, 1 200, 1 600, 2 000, and 4 000 of n MC . 11.4. Steady and unsteady stochastic flow (i)-Steady (time-stationary) stochastic flow. For all the possible values of random control parameters ν CFD ∈ [1 , 3] × 10 -5 and V ref ∈ [0.1 , 0.3] m/s of the prior probability model, the steady flow (time-stationary regime of the stochastic flow) has been obtained at t = 300 (almost-sure stochastic solution of the Navier-Stokes equations).

Figure 7 :

 7 Figure 7: Reference unsteady stochastic flow at t = 30: (a) mean value x → 1 000 × m P (x, t); (b) standard deviation x → 1 000 × σ P (x, t); (c) pdf r → p ρ N d (r) of the random variable ρ N d ; (d) L 2 -norm N d → ρ L 2 (N d ) of the random residual ρ.

Figure 8 :

 8 Figure 8: For N d = 100, n MC = 400, unsteady stochastic flow at t = 30, and for the PLoM computation under the constraint on ρ but without the constraints on W: (a) function i → err R (i); (b) function i → err W (i); (c) pdf w → p W 1 (w) and (d) pdf w → p W 2 (w) estimated with N d = 100realizations (prior model) (dashed line), with n MC = 400 additional realizations (learning) calculated using algo 1 of the PLoM under the constraint on ρ but without the constraints on W (thin solid line), and with n MC = 400 additional realizations (learning) calculated using algo 3 of the PLoM under the constraints on ρ and W (thick solid line).

Figure 9 :

 9 Figure9: For unsteady stochastic flow at t = 30, for N d = 100, and using algo 3 for the iteration algorithm, evolution of the error functions as a function of n MC : i → err R (i) (a), i → err W (i) (b), and i → err RW (i) (c).

Figure 10 :

 10 Figure 10: For unsteady stochastic flow at t = 30 and for N d = 100, graph of pdf r → p ρ (r) of ρ as a function of n MC , estimated (a) with the PLoM without constraint, (b) with the PLoM under constraints with algo 3, and (c) for n MC = 4 000, the PLoM without and under constraints.

Figure 11 :

 11 Figure 11: For unsteady stochastic flow at t = 30, N d = 100, n MC = 4 000, and using algo 3, evolution of the error functions i → err R (i) (a), i → err W (i) (b), and i → err RW (i) (c).each component W 1 and W 2 of random vector W. In each figure is plotted the estimated pdf, on the one hand with the N d = 100 realizations from the training dataset (prior pdf estimated with the training), on the other hand with the n MC = 4 000 additional realizations (learning) estimated with the PLoM without constraint and with the PLoM under the constraints on ρ and W. The mean value E ρ of ρ is 3.23 × 10 -3 , its standard deviation is 2.74 × 10 -3 , and ρ L 2 is 4.23 × 10 -3 . These values must be compared with those of the reference given in.To illustrate the distribution of the realizations generated by the PLoM under constraints, Fig.13displays the 4 000 additional realizations of components of X c = (Q c , W c ). A concentration of these realizations can be observed.

Figure 12 :

 12 Figure 12: For N d = 100 and unsteady stochastic flow at t = 30, pdf w → p W 1 (w) (a) and pdf w → p W 2 (w) (b), estimated with the N d = 100 realizations (prior model) from the training dataset (solid line), with n MC = 4 000 additional realizations (learning) estimated with the PLoM without constraint (dashed line) and under the constraints (thick solid line).

Figure 13 :

 13 Figure 13: 3D-plot of realizations (w c ar,1 , w c ar,2 , q c ar,2 ) of (W c 1 , W c 2 , Q c 2 ) (left figure) and (w c ar,1 , w c ar,2 , q c ar,7 ) of (W c 1 , W c 2 , Q c 7 ) (right figure), computed with the PLoM under constraints: n MC = 4 000 additional realizations symbols) and N d = 100 points from the training dataset (blue symbols).

Figure 14 :

 14 Figure14: For unsteady stochastic flow at t = 30 and for N d = 100, 2D-plot of mean value x → 1 000 × m P (x, t) (a) and of standard deviation x → 1 000 × σ P (x, t) (b) of the pressure random field corresponding to the reference (black thin solid lines) and estimated with n MC = 4 000 additional realizations computed with the PLoM without constraint (blue thick dashed lines).

Figure 15 :

 15 Figure 15: For unsteady stochastic flow at t = 30 and for N d = 100, 2D-plot of mean value x → 1 000 × m P (x, t) (a) and of standard deviation x → 1 000 × σ P (x, t) (b) of the pressure random field corresponding to the reference (black thin solid lines) and estimated with n MC = 4 000 additional realizations computed with the PLoM under constraints (red thick solid lines). Fig. (c) is a zoom of Fig. 14-(b). Fig. (d) is a zoom of Fig. (b).

E 11 =

 11 E 22 = 169 × 10 9 N/m 2 , E 33 = 130 × 10 9 N/m 2 , with Poisson ratios ν 23 = 0.36, ν 31 = 0.28, ν 12 = 0.064, and with shear moduli G 23 = G 31 = 79.6 × 10 9 N/m 2 and G 12 = 50.9 × 10 9 N/m 2 . Let [C] ∈ M + 6 be the Kelvin representation of the fourth-order elasticity tensor of this medium, which belongs to the class of orthotropic symmetry. The mass density is 2 330 Kg/m 3 . There are three types (iii) Random vector W. Let W = (W 1 , . . . , W n w ) be the R n w -valued random variable defined by W = (W s , W C ) with n w = 3 + 21 = 24. For constituting the training dataset, the N d realizations w d , = 1, . . . , N d of W are generated with its prior probability model. The random stiffnesses K s and the random elasticity matrix [C] are defined by a given R n w -valued mapping f param on R n w such that K s , [C] = f param (W) .

  2 β [M] ϕ β in which w = E W . We have ω 1 = 2π × 4.1341 × 10 5 rad/s and ω mode = 2π × 1.2400 × 10 8 rad/s. There are 25 elastic modes in the frequency band B o . The computational model of matrix [D] ∈ M +N is derived from the one proposed Page 93 in[START_REF] Ohayon | Structural Acoustics and Vibration: Mechanical Models, Variational Formulations and Discretization[END_REF], and is written as[D] = n mode β=1 2 χ d (ω β -ω n mode ) [M] ϕ β ([M] ϕ β ) T + 2 χ d ω n mode [M], in which χ d = 0.02. This damping model yields a constant damping rate χ d for the first n mode elastic modes and yields the damping rates χ d ω n mode /ω n mode +1 , . . . , χ d ω n mode /ω N for the Nn mode higher elastic modes.

Figure 16 :

 16 Figure16: For N d = 120, for Tt 0 = 4.4903 × 10 -5 s, and for n time = 6 466, KL expansion of nonstationary stochastic process Y: (a) distribution of the eigenvalues α → Λ α and (b) error function n q → err KL (n q ). PCA of random vector X: (c) error function ν → err PCA (ν).

(

  vi) Calculation of random residual ρ using the subsampling times . Let b(t n ) = (1/N d ) N d =1 y (t n ) be the estimate of the mean value of Y(t n ) computed with the training dataset. The subsampling times τ n s , n s = 1, . . . , n sp ⊂ t n , n = 1, . . . , n time are chosen as the n sp sampling times for which b(τ 1 ), . . . , b(τ n sp ) correspond to the n sp largest values of the set b(t 1 ), . . . ,

Figure 18 -Figure 17 :

 1817 Figure 17: Reference stochastic solution for N d = 1 000: for Obs j , j = 1 (left figure) and j = 6 (right figure), graph of the deterministic functionν freq → dB lin, j (ν freq ) (blue thin line), and graph of the confidence region (yellow region with black thin lines for the upper and lower envelopes) of the random function ν freq → dB j (ν freq ). The black thick line inside the yellow region is the estimate of mean value ν freq → E{dB j (ν freq )}.

Figure 18 :Figure 19 :

 1819 Figure 18: Reference stochastic (a) evolution of the pdf r → p ρ N d (r) of ρ N d as a function of the number N d of realizations in the training dataset; (b) N d → ρ L 2 (N d ) of the random residual ρ.

Figure 20 :

 20 Figure 20: For N d = 120 and with n MC = 1 000 for the PLoM computation under the constraint on ρ but without constraint on W (algo 1): (a) function i → err R (i) and (b) function i → err W (i).

12. 6 .

 6 Results of the PLoM without and under constraints for a given dimension of the training datasetThe training dataset has N d = 120 realizations and n MC N d additional realizations are generated with the PLoM under constraints using algorithm algo 3.

Figure 21 :

 21 Figure 21: For N d = 120 and using algo 3, graph of pdf r → p ρ (r) of ρ as a function of n MC , estimated with the PLoM without constraint (a), under constraints (b), and for n MC = 1 000, comparison of the estimates obtained with the PLoM without and under constraints (c).

Figure 23 Figure 22 :

 2322 Figure23shows the pdf w → p W 2 (w) and the pdf w → p W 5 (w). In each figure is plotted the estimated pdf, on the

Figure 23 :

 23 Figure 23: For N d = 120, pdf w → p W 2 (w) (a) and pdf w → p W 5 (w) (b), estimated with the N d = 120 realizations (prior model) from the training dataset (black thin solid line), with n MC = 1 000 additional realizations (learning) estimated with the PLoM without constraint (blue dashed line) and under the constraints by using algo 3 (blue thick solid line).

Figure 24 :

 24 Figure 24: 3D-plot of realizations (w c ar,1 , w c ar,2 , q c ar,1 ) of (W c 1 , W c 2 , Q c 1 ) (left figure) and (w c ar,4 , w c ar,10 , q c ar,1 ) of (W c 4 , W c 10 , Q c 1 ) (right figure), computed with the PLoM under constraints: n MC = 1 000 additional realizations (red symbols) and N d = 120 points from the training dataset (blue symbols).

Figure 25 :

 25 Figure 25: For N d = 120, for n MC = 1 000 additional realizations generated by the PLoM under constraints, for Obs j , j = 1 (left figure) and j = 6 (right figure), graph of the confidence region (yellow region with thin lines for the upper and lower envelopes) of the random function ν freq → dB j (ν freq ). The thick line inside the yellow region is the estimate of mean value ν freq → E{dB j (ν freq )}.

  Generation of n MC realizations of H λ . The Störmer-Verlet algorithm (see Appendix A) is used for solving Eqs. (53) to (55) as proposed in [52, 30], which allows for generating n MC realizations [z 1 λ ], . . . , [z n MC λ ] in M ν,m of random matrix [Z λ ]. It can then be deduced n MC additional realizations η 1 λ , . . . , η n MC

λ in R ν of random vector H λ such that, for k ∈ 1, . . . , ν and for ∈ 1, . . . , n MC , we have η

  The subsampling times τ n s , n s = 1, . . . , n sp (see Section 4.2) is built with a constant step ∆τ = (Tt 0 )/n sp and n sp = 50. The family U(τ n s ), U(τ n s ), P(τ n s ), n s = 1, . . . , n sp is extracted from the family U(t n ), U(t n ), P(t n ), n = 1, . . . , n time computed with Eqs. (80) to (82). Taking into account Eqs. (78) and (79) at subsampling time τ n s , the R N -valued random residual R(τ n s ), defined by Eq. (17), is rewritten as R(τ n s ) = R NS (τ n s ), R DIV (τ n s ) in which the R N u -valued random variable R NS (τ n s ) and the R N p -valued random variable R DIV (τ n s ) are written, for all n s ∈ 1, . . . , n sp , as R NS (τ n s ) =N NS ( U(τ n s ), U(τ n s ), P(τ n s ), W) a.s. , (83) R DIV (τ n s ) =N DIV (U(τ n s )) a.s. .

	(84)

  1 000, 2 000, and 4 000, the standard deviation of ρ N d is equal to 0.396, 0.383, 0.379, and 0.377, Figure 6: Mean value x → 1 000 × m P (x, t) (a) and standard deviation x → 1 000 × σ P (x, t) (b) estimated with N d = 50 realizations of the training dataset, at t = 300 of the time-stationary pressure random field {P(x, t), x ∈ Ω, t ≥ 300}.

								1											1
								0.8 0.9		9 0 5 -0 .48 71 -0 .7 54 89 -1 .0 22 7 1					-1 .2 90 5 -1 .0 22 7 -0 .7 54 89 -0 .4 87 11 -0 .2 19 32	0.8 0.9		0.3 42 46 0. 17 12 3			36 8 0. 51	0.3 42 46
										.2
								0.7		-1									0.7
								0.5 0.6					-2 .6 29 4	-2.3616	-2.0938 -1. 82 6 -1. 55 82	0.5 0.6		0.68491	0.85614	1.0274
								0.4											0.4
										-1
								0.2 0.3		.0 22 7 -0 .7 54 89 11 -0 .2 19 32 -0 .4 87								0.2 0.3					0. 17 12 3 8 0. 51 36 0. 34 24 6
								0.1											0.1
								0	0	0.2	0.4		0.6		0.8	0	0	0.2	0.4	0.6	0.8
											(a) Mean value									(b) Standard deviation
	1								1										1							1.9	10 -3
	0.8		-0.526 05	-0.7 769 2			-0.7 769 -0. 526 05 2 -0. 27 51 8 -0. 02 43 07 0.2 26 56 0.4 77 43 0.7 28 3 1.2 3	0.8	0.15	962					0.1 59	62	0.3 192 0.4 788 6 4	0.8						N d = 2 000 N d = 400 N d = 1 000	1.85
	0.6		-0.2 751 8		-2. 03 13 -1.7804	-1.5295	-1.2787	-1.02 78	0.6			0.4788 6 0.63848	0.7981 0.95772	1.1173					0.6						N d = 4 000	1.75 1.8
	0.4								0.4							0.15 962 0.319 24			0.4							1.7
	0.2		-0.024 307						0.2										0.2							1.65
	0	0	0.2	0.4	0.6		0.8	0	0	0.2	0.4		0.6	0.8	0	0	0.5		1	1.5	2	2.5	1.6	0	1000	2000	3000	4000

Acknowledgments

This material is based upon work supported by the Defense Advanced Research Projects Agency (DARPA) under Agreement No. HR00111990032 and by FASTMATH, a SciDAC Institute (Scientific Discovery through Advanced Computing) funded by the U.S. Department of Energy (DOE), Office of Science, Advanced Scientific Computing Research (ASCR).

with the boundary conditions that are written for all t ∈ [t 0 , T ], u(x, t) = (0, 0) , x ∈ Γ 0 , u(x, t) = (V CFD , 0) , x ∈ Γ u ,

and with the following initial condition at time t 0 , u(x, t 0 ) = (0, 0) for all x ∈ Ω. We choose ref = 0.9 m and V ref = v V CFD with v = 1 m/s. Consequently, the Reynolds can be rewritten as Re = 0.9V CFD /ν CFD .

Prior probability model of uncertainties

Parametric probabilistic uncertainties and their quantification in computational fluid dynamics have extensively been studied (see for instance [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Bijl | Uncertainty Quantification in Computational Fluid Dynamics[END_REF]). The present application is developed to validate the PLoM under constraints. We must therefore introduce uncertainties and we therefore choose a parametric probabilistic approach. The uncertain control parameters are V CFD and ν CFD . They are modeled by two independent positive-valued random variables, defined on the probability space (Θ, T , P), with values in [V min , V max ] and [ν min , ν max ], written as Consequently, the support of the probability distribution of the random Reynolds number Re is the interval [3 000, 27 000]. The uniform random variables are classically written as

2) /2, in which W 1 and W 2 are independent normalized Gaussian random variables. Let W = (W 1 , W 2 ) be the R n w -valued random variable with n w = 2. For j = 1, 2, let w d, j , = 1, . . . , N d be N d realizations generated with the prior probability model for constituting the training dataset. The random Reynolds number Re is defined by a given real-valued mapping f Re on R 2 such that Re = f Re (W) .

11.3. Computational model, KL expansion, random residual, and PLoM parameters (i) Time discretization scheme. We use the time discretization scheme proposed in [START_REF] Griebel | Numerical Simulation in Fluid Dynamics: a Practical Introduction[END_REF]. choose t 0 = 0. Let ∆t be the time step, t n = n ∆t, n = 1, . . . n time be the time sampling, and T = n time ∆t be the final time. Let us assumed that the velocity field x → u (n) (x) = u(x, t n ) and the pressure field x → p (n) (x) = p(x, t n ) are known at the n-th sampling time and satisfy Eqs. [START_REF] Ohayon | Structural Acoustics and Vibration: Mechanical Models, Variational Formulations and Discretization[END_REF] and (73). The following scheme allows for calculating the approximation of the solution at time t n+1 (below, the x dependence is removed). In Eq. ( 72), the partial time derivative of u at time t n+1 , which is discretized as (u (n+1)u (n) )/∆t, is written as the sum of three terms: the nonlinear force that is explicitly treated, the viscosity force that is implicitly treated, and the pressure correction term that is implicitly treated for enforcing incompressibility defined by to Eq. ( 73),

The implicit equation for the calculation of p (n+1) is obtained by substituting u (n+1) given by the third Eq. ( 75) into the fourth and gives the elliptic equation,

(ii) Spatial discretization. The spatial discretization, the discretization of the spatial derivatives, and the implementation of the boundary conditions, are performed as explained in [START_REF] Griebel | Numerical Simulation in Fluid Dynamics: a Practical Introduction[END_REF]. A grid is defined using n x 1 grid points in the x 1 -axis with a spatial step h x 1 = x 1 /n x 1 and n x 2 grid points in the x 2 -axis with

The spatial discretization is performed on a staggered grid for which u 1 is placed on the x 2 (vertical) cell interfaces, u 2 is placed on the x 1 (horizontal) cell interfaces, and p is placed in the cell midpoints. We then have, for the interior points, (iii) Time-dependent excitation. The base Γ 0 of the mobile part is submitted to a given acceleration in the x 1 -direction. The time-dependent modulation is γ(t) = γ 0 sin t(ω c +∆ω c /2) -sin t(ω c -∆ω c /2) /(π t), in which γ 0 = 120 m/s 2 , ω c = 2 π × 13 × 10 6 rad/s is the central angular frequency and where ∆ω c = 2 π × 10 × 10 6 rad/s is the bandwidth. In the positive-frequency domain, the energy of signal γ(t) is mainly concentrated in the frequency band [ω min , ω e ], with ω min = 2π × 8 × 10 6 rad/s and ω e = 2π × 18 × 10 6 rad/s.

(iv) Time interval of analysis, frequency band of analysis, and observations. This nonlinear dynamical system is analyzed for t ∈ [t 0 , T ] with t 0 = -1.3889 ×10 -5 s and T = 3.1014 × 10 -5 s. At time t 0 , the device is at rest. The frequency band of analysis is [0 , ω max ] with ω max = 2π × 72 × 10 6 rad/s. The frequency band of observation is

where ω o = 2 π × 70 × 10 6 rad/s. If the dynamical system was linear, the energy of the response signal would be concentrated in the same frequency band [ω min , ω e ] as that of the energy of the excitation signal. Due to the nonlinearity however, part of the energy of the excitation signal is transferred outside its frequency band and consequently, the frequency band of the response is not [ω min , ω e ] but [0, ω max ], where ω max > ω e . The results will be presented for two selected observations Obs 1 and Obs 6 . The first one is the x 1 -dof, denoted as dof 1 , of the node of coordinates (14.5, 8.0, 0.0) × 10 -6 m located on the vertical beam of the mobile part. The second one is the x 1 -dof, denoted as dof 6 , of the node of coordinates (23.0, 17.0, 0.0) × 10 -6 m located at the down right corner of the suspended part. The results will be presented in the frequency domain (in Hz) for each one of the two observations and will be the random function ν freq → dB j (ν freq ) = log 10 (| Y dof j (2πν freq ) |) for j = 1 and j = 6, in which ω → Y dof j (ω) is the Fourier transform of the displacement t → Y dof j (t).

12.2. Prior probability model of uncertainties (i) Prior probability model of the springs stiffnesses. Let K s = (K s 1 , K s 2 , K s 3 ) be the R 3 -valued random variable defined on a probability space (Θ, T , P). For j ∈ {1, 2, 3}, K s j is the random stiffness of the spring of type j, with mean value k s j and with a uniform probability distribution on [k s j (1 -

, in which δ s = 0.2 for all j. Random variable K s j can then be written as

2) /2, where W s j is a normalized Gaussian random variable. The random variables W s 1 , W s 2 , W s 3 are assumed to be independent and we introduce the R 3 -valued random variable W s = (W s 1 , W s 2 , W s 3 ).

(ii) Prior probability model of the elasticity tensor of the silicon. The material in Ω is a random homogeneous medium. The Kelvin matrix representation [C] of the random fourth-order elasticity tensor is defined on (Θ, T , P), independent of W s , with values in M + 6 . Its mean value is matrix

can be constructed either in the class of orthotropic symmetry [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF] or both in the orthotropic symmetry class and the anisotropic symmetry class [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF]. Here, we consider the statistical fluctuations in the anisotropy symmetry class [START_REF] Soize | Non gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF] and consequently,

is a random (6 × 6) upper triangular matrix defined as follows [START_REF] Soize | Uncertainty Quantification[END_REF]. Let {W i j } 1≤i≤ j≤6 be 21 independent normalized Gaussian real-valued random variables defined on (Θ, T , P), independent of W s . The Finally, we introduce the random vector

Appendix A. Störmer-Verlet algorithm for solving Eqs. [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations[END_REF] to [START_REF] Kullback | On information and sufficiency[END_REF] The ISDE defined by Eqs. ( 53) to ( 55) is solved for r ∈ [0, r max ] with r max = ( 0 + n MC × M 0 ) ∆r in which ∆r is the sampling step, where 0 is chosen in order that the solution of Eqs. ( 53) to (55) has reached the stationary regime, and where M 0 allows for controlling the number of sampling steps between two consecutive realizations for obtaining a quasi-independence of the computed realizations. Taking as the origin of , the value 0 , for = 0, 1, . . . , n MC × M 0 , we consider the sampling points r = ∆r and the following notations:

, and [W wien ] = [W wien (r )]. The Störmer-Verlet scheme is written, for = 0, 1, . . . , n MC × M 0 , as

with the initial condition defined by Eq. ( 55), where β = f 0 ∆r /4, and where [L + 1 2 ] is the M ν,m -valued random variable such that [L + The increments are statistically independent. For all α = 1, . . . , ν and for all j = 1, . . . , N d , the real-valued random variables [∆W wien +1 ] α j α j are independent, Gaussian, second-order, and centered random variables such that E [∆W wien +1 ] α j [∆W wien +1 ] α j = ∆r δ αα δ j j .