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

R Ghanem
email: ghanem@usc.edu

Probabilistic-learning-based stochastic surrogate model from small incomplete datasets for nonlinear dynamical systems

Keywords: Probabilistic learning, realizations as targets, statistical inverse problem, Kullback-Leibler divergence, uncertainty quantification degree of freedom. ISDE: Itô stochastic differential equation. KDE: kernel density estimation. KLDMP: Kullback-Leibler divergence minimum principle. PCA: principal component analysis. pdf: probability density function. POD: proper orthogonal decomposition. ROB: reduced order basis. SCM: stochastic computational model. SSM: statistical surrogate model

We consider a high-dimensional nonlinear computational model of a dynamical system, parameterized by a vectorvalued control parameter, in the presence of uncertainties represented by an uncontrolled parameter modeled by a vector-valued random variable, and possibly with stochastic excitation. The objective is to construct a statistical surrogate model where the input is any deterministic value of the control parameter, and the output is a vector-valued observation of the computational model, which is a random vector whose probability measure is updated using a target dataset. To construct this statistical surrogate model, the stochastic response of the computational model must be built, which is a vector-valued time-discretized stochastic process in high dimension, depending on the control parameter. It is assumed that the computational cost of a single evaluation of the deterministic model is high. For the probabilistic updating, we consider a subset of the components of the observation of the computational model, defined as the "identification observation" of the computational model, for which a small target dataset is available. Therefore, the target dataset is associated with partial observability, corresponding to an incomplete data case. Given a prior probability model of the random control and uncontrolled parameters, a training dataset is constructed, consisting of realizations of the random triplet composed of the stochastic response, the random identification observation, and the random control parameter. Since the computational cost of a single evaluation of the deterministic model is assumed to be large, the training dataset is also of small size. The main challenges in this problem are the high dimensionality, partial observability leading to incomplete data in the target dataset for the identification observation of the computational model (which is not sufficient to identify the computational stochastic responses), and the availability of a small training dataset. To address these challenges, we propose a methodology based on statistical methods for constructing necessary reduced representations, direct probabilistic learning under constraints using probabilistic learning on manifolds (PLoM) constrained by the target dataset, and the use of a weak formulation of the Fourier transform of probability measures. Statistical conditioning is also employed to explore the learned dataset. The constructed predictive statistical surrogate model can be implemented in the context of online computation. We apply this approach to a problem of nonlinear stochastic dynamics in high dimensions within the framework of deformable solids mechanics.

Introduction

Context of the paper

The development of surrogate models for parameterized large-scale computational models is challenging and an extremely active research topic, resulting in a huge number of publications. The scientific community has proposed several classes of methods, including deterministic representations, probabilistic/statistical-based approaches, and more recently, Machine Learning tools, with or without probabilistic/statistical formulations. It is not possible here, nor is it the purpose, to attempt to review this vast field. Instead, we will focus on aspects directly related to the method proposed in this paper. After referencing some general works on surrogate modeling, we will specifically discuss methods for handling incomplete data, approaches for model reduction in the presence of uncertainties, and probabilistic approaches developed for small data. These three components are utilized in the proposed probabilisticlearning-based stochastic surrogate model from small incomplete datasets.

(i) Concerning the general aspects of the surrogate models, their construction can be carried out using a parametric and nonparametric approaches in a deterministic or a probabilistic framework (see for instance [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF][START_REF] Alizadeh | Managing computational complexity using surrogate models: a critical review[END_REF]).

-Typically the parametric approaches can be based on deterministic polynomial regressions [START_REF] Peixoto | A property of well-formulated polynomial regression models[END_REF][START_REF] Ostertagová | Modelling using polynomial regression[END_REF][START_REF] Kleijnen | Regression and kriging metamodels with their experimental designs in simulation: a review[END_REF], kriging approaches [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF][START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF][START_REF] Kersaudy | A new surrogate modeling technique combining kriging and polynomial chaos expansions-application to uncertainty analysis in computational dosimetry[END_REF][START_REF] Qian | A sequential constraints updating approach for kriging surrogate model-assisted engineering optimization design problem[END_REF][START_REF] Zhou | Surrogate modeling of high-dimensional problems via data-driven polynomial chaos expansions and sparse partial least square[END_REF], polynomial chaos representations (see for instance [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF][START_REF] Doostan | Red-Horse, Stochastic model reduction for chaos representations[END_REF][START_REF] Das | Asymptotic sampling distribution for polynomial chaos representation from data: a maximum entropy and fisher information approach[END_REF][START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF][START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF][START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF][START_REF] Tipireddy | Basis adaptation in homogeneous chaos spaces[END_REF][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF][START_REF] Abraham | A robust and efficient stepwise regression method for building sparse polynomial chaos expansions[END_REF][START_REF] Thimmisetty | Homogeneous chaos basis adaptation for design optimization under uncertainty: Application to the oil well placement problem[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF], for general methods, [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF][START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansions[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Rosić | Sampling-free linear Bayesian update of polynomial chaos representations[END_REF][START_REF] Madankan | Polynomial-chaos-based Bayesian approach for state and parameter estimations[END_REF][START_REF] Chen-Charpentier | Parameter estimation using polynomial chaos and maximum likelihood[END_REF][START_REF] Elsheikh | Efficient bayesian inference of subsurface flow models using nested sampling and sparse polynomial chaos surrogates[END_REF][START_REF] Giraldi | Bayesian inference of earthquake parameters from buoy data using a polynomial chaos-based surrogate[END_REF][START_REF] Tsilifis | Bayesian adaptation of chaos representations using variational inference and sampling on geodesics[END_REF] for integrating data using the maximum likelihood or the Bayesian method, [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Shao | Bayesian sparse polynomial chaos expansion for global sensitivity analysis[END_REF][START_REF] Luthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF]) for sparse modeling, and [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF][START_REF] Tsilifis | Reduced Wiener chaos representation of random fields via basis adaptation and projection[END_REF][START_REF] Tsilifis | Bayesian adaptation of chaos representations using variational inference and sampling on geodesics[END_REF][START_REF] Liu | Surrogate modeling based on resampled polynomial chaos expansions[END_REF][START_REF] Kontolati | Manifold learning-based polynomial chaos expansions for high-dimensional surrogate models[END_REF] for polynomial chaos on manifolds.

-The nonparametric approaches can be based on kernel-based regression (see for instance [START_REF] Akian | Learning best kernels from data in gaussian process regression. with application to aerodynamics[END_REF]), and more generally on spectral approaches, including deterministic projection methods based on the choice or the construction of a deterministic reduced-order basis (ROB) (such as a Hilbert basis adapted to the operators of the problem under consideration; see for instance [START_REF] Nouy | Generalized spectral decomposition for stochastic nonlinear problems[END_REF][START_REF] Nouy | Low-rank tensor methods for model order reduction[END_REF]) or ROB constructed in the probabilistic context, such as principal component analysis (PCA), Proper Orthogonal Decomposition (POD), and Karhunen-Loéve expansion (these three methods belonging to the same mathematical class of the spectral methods). A direct construction of a random ROB has also been proposed (see below).

-We must obviously mention the most emerging approaches, those of Machine Learning (ML), both parametric and nonparametric, that is to say the deep learning for artificial neural networks, which are effective and validated only for "complex mapping" when big data are available. In this ML context, many efforts are performed to develop probabilistic/statistical surrogate models for the small dataset cases. This last approach is the one followed in the present work.

(ii) The developments of methods in Machine Learning for the case of incomplete datasets are not recent and have given rise to a huge number of publications with different fields of application. It is therefore not a question here of offering a review, even a very partial one, on this vast subject. We only give a few aspects and references to situate the method that is proposed to take into account incomplete data. Two main classes can be identified, statistical methods and methods based on reduced-order representation. For the statistical methods, we can refer, for instance, to methods based on the likelihood and Bayesian approaches [START_REF] Ghahramani | Learning from incomplete data, MIT A.I. Memo 1509[END_REF], to regression analyses [START_REF] Harel | The estimation of r 2 and adjusted r 2 in incomplete data sets using multiple imputation[END_REF], to regression neural network ensembles for multiple imputation [START_REF] Gheyas | A neural network-based framework for the reconstruction of incomplete data sets[END_REF][START_REF] Hittawe | Efficient sst prediction in the red sea using hybrid deep learningbased approach[END_REF], to Bayesian-based density-based clustering approaches [START_REF] Xue | Effective density-based clustering algorithms for incomplete data[END_REF], or to an adaptive Bayesian SLOPE with missing values [START_REF] Jiang | Adaptive bayesian SLOPE: model selection with incomplete data[END_REF]. Concerning the reduced-order representation methods, see, for instance, [START_REF] Stacklies | pcamethods -a bioconductor package providing pca methods for incomplete data[END_REF][START_REF] Podani | Principal component analysis of incomplete data-a simple solution to an old problem[END_REF] for PCA-based methods, [START_REF] Caiafa | Decomposition methods for machine learning with small, incomplete or noisy datasets[END_REF] for decomposition methods, [START_REF] Cai | Subspace estimation from unbalanced and incomplete data matrices: l2,8 statistical guarantees[END_REF] for spectral methods (in fact, these three methods belong to the large class of the spectral methods). In this paper, the context of incomplete data differs from that which is generally considered in the very large number of published papers relating to this subject. Indeed, we consider dynamical systems that are under observed and we do not seek to complete the "incomplete experimental data", but to take these "experimental data" into account to update the probability measure of all the observations of the dynamical system using the Kullback-Leibler divergence minimum principle with respect to the prior probability measure constructed using only the points of the training dataset. To take into account this framework of incomplete data, the proposed method is based both on a projection derived from a PCA representation and on the use of probabilistic learning based on the Kullback minimum principle. The method is therefore, straddling statistical methods and reduced-order representation methods.

(iii) The projection-based reduced-order computational model (ROM) that has a very small dimension with respect to the one of the large-scale computational model is a very attractive and efficient method in nonlinear computational dynamics (see for instance [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Carlberg | Efficient non-linear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF][START_REF] Carlberg | The gnat method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF][START_REF] Farhat | Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency[END_REF][START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Paul-Dubois-Taine | An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models[END_REF]). It should be noted that a hyperreduction method allows for achieving computational efficiency in arbitrarily nonlinear parameterized settings [START_REF] Farhat | Computational bottlenecks for PROMS: precomputation and hyperreduction, Model order reduction[END_REF]. However, in order to obtain a robust ROM against uncertainties, which is required, for example, to transform such a ROM into a digital twin of the physical system [START_REF] Jones | Characterising the digital twin: A systematic literature review[END_REF][START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF], it is necessary to take into account the model uncertainties induced by modeling errors so that the computational model becomes predictive. The nonparametric probabilistic method (NPM) developed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF][START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Farhat | Feasible probabilistic learning method for model-form uncertainty quantification in vibration analysis[END_REF][START_REF] Wang | Modeling uncertainties in molecular dynamics simulations using a stochastic reduced-order basis[END_REF][START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF][START_REF] Azzi | Acceleration of a physics-based machine learning approach for modeling and quantifying model-form uncertainties and performing model updating[END_REF][START_REF] Zhang | A riemannian stochastic representation for quantifying model uncertainties in molecular dynamics simulations[END_REF] is a way for taking into account model-form uncertainties in nonlinear computational model. The NPM consists in substituting the reduced-order basis by a random basis whose hyperparameters of its prior probability model can be identified from targets, and therefore, allows the model to be enriched using data. We have explained above this approach that couples the parameterized reduced-order models with the NPM formulation of the model-form uncertainty, enriched by target datasets, because the method proposed is an alternative way which makes it possible to build a surrogate predictive statistical model that can also be used online. The proposed method will couple the reduced-order representation and the probabilistic learning under the constraints defined by the target datasets (see below).

(iv) The probabilistic learning is also a very active domain of research for constructing surrogate models (see for instance, [START_REF] Talwalkar | Large-scale manifold learning[END_REF][START_REF] Öztireli | Spectral sampling of manifolds[END_REF][START_REF] Marzouk | Sampling via measure transport: An introduction[END_REF][START_REF] Parno | Transport map accelerated markov chain Monte Carlo[END_REF][START_REF] Perrin | Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints[END_REF][START_REF] Kevrekidis | Manifold learning for parameter reduction[END_REF][START_REF] Pan | Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability[END_REF][START_REF] Kalogeris | Diffusion maps-based surrogate modeling: An alternative machine learning approach[END_REF]). In this context, the probabilistic learning on manifolds (PLoM) method has specifically been developed for the small dataset cases [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Entropy-based closure for probabilistic learning on manifolds[END_REF][START_REF] Soize | Sampling of Bayesian posteriors with a non-Gaussian probabilistic learning on manifolds from a small dataset[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF] for which several extensions have been proposed to take into account implicit constraints induced by physics, computational models, and measurements [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF], to reduce the stochastic dimension using a statistical partition approach [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF], and to update the prior probability measure by a target dataset whose points are, for instance, experimental realisations of the system observations [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF]. This last capability of PLoM can also be viewed as an alternative method to the Bayesian inference for the high dimension [START_REF] Kennedy | Bayesian calibration of computer models[END_REF][START_REF] Marzouk | Stochastic spectral methods for efficient Bayesian solution of inverse problems[END_REF][START_REF] Gentle | Computational statistics[END_REF][START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF][START_REF] Owhadi | On the brittleness of Bayesian inference[END_REF][START_REF] Matthies | Inverse problems in a Bayesian setting[END_REF][START_REF] Dashti | The Bayesian approach to inverse problems[END_REF][START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF][START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF][START_REF] Perrin | Adaptive method for indirect identification of the statistical properties of random fields in a Bayesian framework[END_REF] and is a complementary approach to existing methods in machine learning for sampling distributions on manifolds under constraints. It allows for solving unsupervised and supervised problems under uncertainty for which the training sets are small. This situation is encountered in many problems of physics and engineering sciences with expensive function evaluations. The exploration of the admissible solution space in these situations is thus hampered by available computational resources. The PLoM was successfully adapted to tackle these challenges for several related problems including nonconvex optimization under uncertainty [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF][START_REF] Soize | Design optimization under uncertainties of a mesoscale implant in biological tissues using a probabilistic learning algorithm[END_REF][START_REF] Ghanem | Optimal well-placement using probabilistic learning[END_REF][START_REF] Farhat | Feasible probabilistic learning method for model-form uncertainty quantification in vibration analysis[END_REF][START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF][START_REF] Capiez-Lernout | Nonlinear stochastic dynamics of detuned bladed disks with uncertain mistuning and detuning optimization using a probabilistic machine learning tool[END_REF][START_REF] Almeida | A probabilistic learning approach applied to the optimization of wake steering in wind farms[END_REF], fracture paths in random composites [START_REF] Guilleminot | Data-driven enhancement of fracture paths in random composites[END_REF], updating digital twins under uncertainties [START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF], calculation of Sobol's indices [START_REF] Arnst | Computation of sobol indices in global sensitivity analysis from small data sets by probabilistic learning on manifolds[END_REF], dynamic monitoring [START_REF] Soize | Machine learning for detecting structural changes from dynamic monitoring using the probabilistic learning on manifolds[END_REF], surrogate modeling of structural seismic response [START_REF] Zhong | Surrogate modeling of structural seismic response using Probabilistic Learning on Manifolds[END_REF], for the waterflooding in oil reservoir [START_REF] Almeida | Uncertainty quantification of waterflooding in oil reservoirs computational simulations using a probabilistic learning approach[END_REF]. As we indicated previously, the third ingredient of the proposed method to build the predictive surrogate model is a probabilistic learning method based on the Kullback-Leibler divergence minimum principle that allows the target dataset to be integrated for updating the prior probability measure that is built with the training set. In this context, the PLoM algorithm is also used to avoid the scattering of learned realizations associated with the updated probability measure in order to preserve its concentration in the neighborhood of the random manifold defined by the stochastic computational model.

Novelty of the paper

The novelty of this paper consists in the development of a methodology for constructing a predictive statistical surrogate model to represent any parameterized, uncertain, stochastic, nonlinear computational model in cases where the dimension is high, there is partial observability leading to incomplete data for the small target dataset of identification observations, the nonlinear mapping that computes the identification observations from the computational stochastic responses is not injective, and only a small training dataset is available. Additionally, the stochasticity of the random responses is caused by uncontrolled random parameters and stochastic excitations in the nonlinear computational model. Machine learning formulations based on artificial neural networks are not well-suited for such cases due to the high dimension, small training dataset, small target dataset, incomplete data, random uncontrolled parameters, and stochastic excitations in the nonlinear computational model. The proposed approach and algorithms are based on statistical methods, adapted reduced representations, direct probabilistic learning under constraints using PLoM constrained by the target dataset (utilizing a weak formulation of the Fourier transform of probability measures), and an effective description of a predictive statistical surrogate model using conditional statistics. These statistics explore the learned dataset and can be carried out online for any given value of the control parameter without invoking the stochastic computational model.

Organization of the paper

The framework of the problem under consideration and a summary of the methodology have been presented in Section 2. Additional developments primarily focused on convergence analyses and algorithms are provided in Section 3. An application is then carried out in Section 4 and involves the mechanical system described in [START_REF] Soize | Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics[END_REF], which is used to validate the proposed methodology. This application pertains to a three-dimensional MEMS device in which the nonlinear stochastic dynamics are studied. It is particularly interesting and challenging for the proposed approach as the effects of the nonlinearities on the nonlinear stochastic responses are highly sensitive to the values of the control parameters.

Terminology and notations

(i) Convention used for the variables, vectors, and matrices. A lower-case Latin or Greek letter, such as x or η, is a deterministic real variable. A boldface lower-case Latin or Greek letter, such as x or η, is a deterministic vector. An upper-case Latin letter, such as X, is a real-valued random variable. A boldface upper-case Latin letter, such as X, is a vector-valued random variable. A lower-or upper-case Latin letter between brackets, such as [x] or [X], is a deterministic matrix. A boldface upper-case letter between brackets, such as [X], is a matrix-valued random variable.

(ii) Probability space, random variable, probability measure, and probability density function. For any finite integer m ≥ 1, the Euclidean space R m is equipped with the σ-algebra B R m . If Y is a R m -valued random variable defined on the probability space (Θ,

T , P), Y is a mapping θ → Y(θ) from Θ into R m , measurable from (Θ, T ) into (R m , B R m ), and Y(θ) is a realization (sample) of Y for θ ∈ Θ.
The probability distribution of Y is the probability measure P Y (dy) on the measurable set (R m , B R m ) (we will simply say on R m ). The Lebesgue measure on R m is noted dy and when P Y (dy) is written as p Y (y) dy, p Y is the probability density function (pdf) on R m of P Y (dy) with respect to dy. Finally, E denotes the mathematical expectation operator and a.s. means "almost surely".

(iii) Algebraic notations. N, R: set of all the integers {0, 1, 2, . . .}, set of all the real numbers. R n : Euclidean vector space on R of dimension n. M n,m : set of all the (n × m) real matrices. M n : set of all the square (n × n) real matrices. M + n : set of all the positive-definite symmetric (n × n) real matrices. Control parameter: the control parameter, w, of the computational model corresponds to the "input" of the "statistical surrogate model" (SSM). It is a deterministic parameter that belongs to an admissible set C w ⊂ R n w . For the construction of the SSM, the control parameter w is modeled by a random variable W whose realizations are generated with its given prior probability measure for which its support is the admissible set C w . These realizations allow the training dataset to be constructed. Uncontrolled parameter: there is an internal vector-valued parameter u in the computational model that is not used as an input or an output of the statistical surrogate model. For the construction of the SSM, the uncontrolled control parameter u is modeled by a random variable U whose realizations are generated with its given prior probability measure. These realizations allow the training dataset to be constructed. Observation: the vector-valued "observation" of the computational model corresponds to the "output" of the statistical surrogate model (the quantity of interest). Concerning the prediction performed with the SSM, for a deterministic value w of the control parameter given in its admissible set, the SSM estimates the probability measure of the "updated random observation" (the conditional random output given w).

[I n ]: identity matrix in M n . x = (x 1 , . . . , x n ): point in R n . ⟨x, y⟩ = x 1 y 1 + . . . + x n y n : inner product in R n . ∥ x ∥: norm in R n such that ∥ x ∥ = ⟨x, x⟩. ∥ [x] ∥ = sup ∥ y ∥ =1 ∥ [x]
Identification observation and target dataset: For the construction of the SSM, a specific vector-valued observation of the computational model is introduced and is defined as the "identification observation". The "target dataset" is constituted of a set of values of the identification observation, which are either experimental measures or other simulated data.

Training dataset: given a prior probability model of the random control and uncontrolled parameters, a training dataset is constructed, consisting of realizations of the random triplet composed of the stochastic response of the computational model, the random identification observation, and the random control parameter.

Learned datasets: there are several learned datasets, which are related to the different random variables introduced in the construction of the SSM. Each "learned dataset" is built from the corresponding "training dataset" associated with each random variable. In particular, the learned dataset for updating the random observation of the computational model, is obtained by constraining the "identification observation" by the "target dataset".

Reference dataset for validation: To validate the SSM, a specific set of values of the control parameter w, denoted w ref , is introduced. This set is defined as the "reference dataset" and is made up of different points from those in the training dataset.

Summary of the proposed methodology

In this section, in order to facilitate the reading of this paper, we present a summary of the proposed method, indicating the framework of the developments and giving the main hypotheses used and the reasons for which they are introduced.

Parameterized nonlinear uncertain stochastic computational model and its observation

We consider a high-dimension nonlinear uncertain stochastic computational model of a dynamical system, parameterized by a vector-valued control parameter w, in presence of uncertainties represented by an uncontrolled parameter that is a R n u -valued random variable U = (U 1 , . . . , U n u ) whose support of its probability distribution is C u ⊂ R n u , and possibly with a stochastic excitation. The control parameter w = (w 1 , . . . , w n w ) is with values in an admissible subset C w of R n w . The time-discretized nonlinear stochastic computational model, simply called stochastic computational model (SCM), corresponds to the time-discretized nonlinear differential equation of order m diff in t, which is written as, N Y(t; w), Y (1)

Identification observation and small target dataset

For the probabilistic updating using a target dataset, a subset of components of O(w) is considered, which is represented by the random variable O id (w) with values in R n o with n o ≪ N o , and which will be called the "identification observation". The corresponding restriction of mapping O is denoted by O id . The random identification observation is then written as O id (w) = O id (Y(. ; w)). Relatively to O id (w), a deterministic target dataset D targ (o id ) = {o 1 targ , . . . , o N r targ } is given. For each r in {1, . . . , N r }, the vector o r targ ∈ R n o corresponds to a "measurement" performed on the dynamical system for a given value w r targ ∈ C w of the control parameter w. Let W targ = {w 1 targ , . . . , w N r targ } be the set of these control parameter values. The target dataset is used to update the probability measure of the random surrogate model whose statistical fluctuations are induced by the modeling errors existing in the SCM and that are simulated thanks to the presence of random vector U that models uncertainties. It should be noted that the target dataset is associated with a partial observability, which thus corresponds to an incomplete data case, that is an important difficulty.

Small training dataset

For constructing the training dataset, a prior probability model is introduced for the random variables W, U, and for the time series F that are assumed to be statistically independent. Let {(w j , u j , f j ), j = 1, . . . , n d } be n d independent realizations of (W, U, F) and let W train = {w 1 , . . . , w n d } be the set of sampled values of W. The training dataset, D train (y, o id , w), is constituted of n d points constructed using the SCM. For every j in {1, . . . , n d }, let y j = y(. ; w j , u j , f j ) be the solution of the deterministic computational model. Once y j is known, the corresponding realization o j id = O id (y j ) of O id (w j ) is computed. The training dataset is thus defined by

D train (y, o id , w) = {(y j , o j id , w j ), j = 1, . . . , n d } with y j = y(. ; w j , u j , f j ) and o j id = O id (y j ) . (2.
2)

It should be noted that y j = {y j (t), t ∈ J} is a function defined on J with values in R n y while (o j id , w j ) is in R n o × R n w . It is assumed that the numerical cost of a single evaluation of the deterministic computational model is large. Therefore the number n d is assumed to be small, a few tens or a few hundreds. Under these conditions the number of points in the training dataset is small and we are in the case of a small training dataset (as opposed to big data). Under these conditions, machine learning formulations based on the learning of artificial neural networks are not suitable for the considered framework of the high dimension and of small training and target datasets. For this reason, we propose a methodology based on statistical methods for constructing reduced representation and also on a direct probabilistic learning under constraints of the probability measures.

Predictive statistical surrogate model

The objective of this work is the construction of a predictive statistical surrogate model defined by the family {O ud (w), w ∈ C w } of random variables with values in R N o , which is the probabilistic updating of the family {O(w), w ∈ C w }. This means that we have to estimate the w-dependent probability measure P O ud (w) (do; w) of the R N o -valued random variable O ud (w) and to develop a generator of independent realizations. This problem is difficult due to the high dimension, to the partial observability inducing incomplete data for the target dataset of the identification observations, to the non-injectivity of nonlinear mapping O, and to small training and target datasets.

2.5. Large learned dataset from the small training dataset, built with a probabilistic learning constrained by the small target dataset Below, we summarize the main steps of the proposed methodology that is developed for circumventing the identified difficulties and we explain the reasons for the choices made.

(i) Reduced representation of Y with a (t, w)-independent ROB in R n y represented by a matrix [V] in M n y ,n q . The R N o - valued random observation O(w)
is expressed as a nonlinear mapping of the finite family {Y(t; w), t ∈ J} of R n y -valued random variables. The probabilistic learning constrained by the target dataset required to perform the learning for a random vector of length n y × n time , for which n d realizations are available, and for which n d × n MC learned realizations should be generated. Therefore, learned data would be represented by a matrix containing (n y × n time ) × (n d × n MC ) 64-bit words. For instance, for n y = 10 5 , n time = 10 4 , n d = 10 2 , and n MC = 10 4 , this matrix requires 10 15 64-bit words! It is then necessary to construct a reduced representation of {Y(t; w), t ∈ J}. However, we want to perform a global reduced representation that yields a good approximation for all t in J and for all w in C w . We then use the existing POD methodology for constructing the reduced representation Y (n q ) (• ; w) of Y(• ; w), which is written as,

Y (n q ) (t; w) = [V] Q(t; w) , ∀t ∈ J , ∀w ∈ C w .
(2.

3)

The matrix [V] ∈ M n y ,n q is independent of t and w, whose columns constitute an orthonormal reduced-order basis (ROB) in R n y (for the usual Euclidean inner product), and where n q ≪ n y is estimated by a convergence analysis.

Let J snps = {τ 1 , . . . , τ n snps } ⊂ J be the subset of snapshots time with n snps < n time (possibly with n snps = n time ) and let [y j ] = [y j (τ 1 ) . . . y j (τ n snps )] be the matrix in M n y ,n snps . Then matrix [V] is constructed by compression of the matrix

[y] = [y 1 . . . y n d ] ∈ M n y ,N d , N d = n d × n snps , [y j ] = [y j (τ 1 ) . . . y j (τ n snps )] ∈ M n y ,n snps , (2.4
) and[y n d ]. For all w in C w , {Q(t; w), t ∈ J} is a finite family of R n q -valued random variables (sometimes called the random vector of the generalized coordinates). The finite family of functions {y 1 , . . . , y n d } with y j : J → R n y is replaced by the finite family {q 1 , . . . , q n d } with q j : J → R n q such that, for all t in J,

and is such that [V] T [V] = [I n q ]. It should be noted that [y] is the collection of matrices [y 1 ], [y 2 ], . . .,
q j (t) = [V] T y j (t) , ∀t ∈ J , ∀ j ∈ {1, . . . , n d } . (2.5)
This means that q j (t) is the projection of y j (t) on the subspace of R n y spanned by the ROB represented by [V].

(ii) Definition of the training dataset D train (x). For all j ∈ {1, . . . , n d }, we define x j ∈ R n x such that

x j = (q j , o j id , w j ) , q j = (q j (t 1 ), . . . , q j (t n time )) ∈ R n time ×n q , n x = n time × n q + n o + n w .

(2.6)

The training dataset relative to the n d points x j is thus D train (x) = {x 1 , . . . , x n d } and is represented by the matrix

[x d ] ∈ M n x ,n d defined by [x d ] = [x 1 . . . x n d ] ∈ M n x ,n d . (2.7)
There are three main difficulties for using a probabilistic learning from the training dataset D train (x) under the constraint defined by the target dataset D targ (o id ). The first difficulty, which is the most important, is that we find ourselves in a case of incomplete data induced by a partial observability. Indeed, the target dataset D targ (o id ) concerns only the part o j id ∈ R n o of the vector x j ∈ R n x . There is not a target value for the q j ∈ R n time ×n q part. The second difficulty is due to the fact that there does not exist a mapping o j id → q j that is the inverse of the non-injective nonlinear mapping q j → o j id . Therefore, we cannot complete the data by this way. Finally, the third difficulty is related to the fact that n d is not equal to N r , and therefore, the possibility to replace x j = (q j , o j id , w j ) by a vector such as (q j , o j id , w j , o r targ ) cannot be done. We thus propose the following approach to bypass these three difficulties.

(iii) Reduced representation of X using a PCA. Let X = (Q, O id , W) be the R n x -valued random variable for which x 1 , . . . , x n d , grouped in matrix [x d ], are n d independent realizations. A principal component analysis (PCA) of X is performed using [x d ].
The reduced representation X (ν) of X is thus obtained for which the mean-square convergence (in the vector space of all second-order R n x -valued random variables) is controlled with respect to the dimension ν of the reduction, which is written as

X (ν) = x + [Φ] [λ] 1/2 H , ν < n x , (2.8)
in which H is a centered R ν -valued random variable with covariance matrix equal to [I ν ]. In general, n o is large and consequently, the reduction will be important yielding ν ≪ n o . The vector x is the mean value of X estimated with x 1 , . . . , x n d . The matrix [Φ] ∈ M n x ,ν and the positive-definite diagonal matrix [λ] ∈ M ν are the eigenvectors and the eigenvalues of the covariance matrix of X, which is estimated with x 1 , . . . ,

x n d . Matrix [Φ] is such that [Φ] T [Φ] = [I ν ]. By construction, we have [x d ] = [x] + [Φ] [λ] 1/2 [η d ] in which [x] = [x . . . x] ∈ M n x ,n d . The columns of the matrix [η d ] ∈ M ν,n d are the n d realizations η 1 , . . . , η n d of random vector H. At mean-square convergence, matrix [η d ] is computed by [η d ] = [λ] -1/2 [Φ] T ([x d ] -[x]
) .

(2.9)

The training dataset for the R ν -valued random variable H is then defined as D train (η) = {η 1 , . . . , η n d }.

(iv) "Projection" of the target on the model in the context of incomplete data due to the partial observability. Now we have to associate a vector η r targ ∈ R n o to o r targ for r = 1, . . . , N r . For that we need to build a mapping that associates a vector η ∈ R ν to each o id ∈ R n o . We now remove superscript (ν) for simplifying the writing. For any realization η in R ν of H, the corresponding realization x of X is given by Eq. (2.8),

x = x + [Φ] [λ] 1/2 η. Since x = (q, o id , w), the extraction of o id ∈ R n o from x ∈ R n x yields o id = o id + [A o ] η , o id ∈ R n o , [A o ] = [Φ o ] [λ] 1/2 ∈ M n o ,ν , [Φ o ] ∈ M n o ,ν .
(2.10)

The matrix

[A o ] ∈ M n o ,ν admits a unique left pseudo-inverse [A inv o ] ∈ M ν,n o The desired mapping is constructed by solving the equation [A o ] η = o id -o id in the linear least-squares sense, which admits the unique solution η = [A inv o ] (o id -o id ). We then have η r targ = [A inv o ] (o r targ -o id ) , r ∈ {1, . . . , N r } . (2.11)
Note that this construction can be viewed as a "projection" of the target onto the model in the context of incomplete data due to partial observability. In order to simplify the following explanation, let us consider that F = f is a deterministic time function. Therefore, if the trace of the covariance matrix of the random uncontrolled parameter U tends to zero, then for each given value w of the control parameter, the level of statistical fluctuations of the stochastic response Y(•, ; w) = y(•, ; w, U, f) will also tend to zero (assuming continuity with respect to U). This means that for small statistical fluctuations of U, if the "distance" (between two clusters) from the training dataset D train (η) to the target dataset D targ (η targ ) is "significant", then the projection of the target dataset onto the model we have presented will only give a weak "contribution," and the constraint will have little effect in the probabilistic updating. To remedy this, the prior probability model of U has to be carefully defined in order to generate enough statistical fluctuations in the stochastic response y(•, ; w, U, F). Thus, the level of statistical fluctuations of U makes it possible to control the "diameter" of the generated family of computational models. This situation is similar (although very different methodologically and also in terms of assumptions and objectives) to that of Gaussian Bayesian inference: for a given level of Gaussian noise, the support of the prior probability measure must be adapted to the domain where the Gaussian likelihood function has significant contributions (which can be a difficult problem in high dimensions).

(v) Probabilistic learning on manifolds (PLoM) constrained by the target dataset. Using the small training dataset D train (η) = {η j , j = 1, . . . , n d } relative to H and using the associated small target dataset D targ (η targ ) = {η r targ , r = 1, . . . , N r }, the probabilistic learning on manifolds (PLoM) is carried out under the constraints defined by the target dataset and based on the use of a weak formulation of Fourier transform of probability measures. This step allows for generating a large learned dataset D learn (η ud ) = {η ℓ ud , ℓ = 1, . . . , N ud } of the R ν -valued random variable H ud that is the updating of H under the constraint defined by D targ (η targ ). The number of points in the learned dataset is

N ud = n d ×n MC in which n MC ≪ 1 is given. (vi)
D learn (o ud , w ud ) = {(o ℓ ud , w ℓ ud ) ∈ R N o × R n w , ℓ = 1, . . . , N ud } .
(2.12)

Predictive statistical surrogate model available in an online computational context

The updated statistical surrogate model is defined by the w-dependent R N o -valued random variable O ud (w). Let us assume that all the considered probability measures involving in the conditional statistical estimations admit densities with respect to the Lebesgue measures do and dw on R N o and R n w . For any w given in C w , the w-dependent probability density function on R N o of the random variable O ud (w) is given by The joint pdf p O ud , W ud (o, w) is estimated using the Gaussian KDE method with the learned dataset D learn (o ud , w ud ) defined by Eq. (2.12). The estimation of the pdf p W ud (w) is deduced by an explicit integration with respect to o of the KDE estimate of p O ud , W ud (o, w). Once steps (v-1) to (v-4) have been performed (offline computation), step (v-5) corresponds to the effective description of the predictive statistical surrogate model for which the output statistics can rapidly be computed in an online context using only conditional statistics on the learned dataset. For any value of the control parameter w given in C w , this predictive statistical surrogate model yields the statistics of the random observations, such as mean values, variances, probability density functions, confidence domains, etc.

p O ud (w) (o; w) = p O ud | W ud (o|w) , o ∈ R N o , w ∈ C w ⊂ R n w , ( 2

Algorithmic complements and convergence analysis

In this section we present additional developments that allow us to specify the convergence criteria and the algorithms to implement the proposed methodology for building the predictive statistical surrogate model. All notations, assumptions, developments introduced in Section 2 are used without repeating them and without systematically referring to this section.

Illustration of a parameterized nonlinear uncertain computational model

As an illustration of the dynamic computational model whose time discretization yields Eq. (2.1), we consider the w-parametric, nonlinear, uncertain, computational model

[M(w, U)] Ÿ(t; w) + g Y(t, w), Ẏ(t; w); w, U = F(t; w) , t ∈ ]t 0 , T ] , (3.1)
arising from a large finite element semi-discretization of a boundary value problem governing the dynamic equilibrium of a mechanical structure, where Y(t; w) = y(t; w, U, F) represents the n y displacement dofs at time t, Ẏ(t; w) = dY(t; w)/dt is the velocity vector and Ÿ(t; w) = d 2 Y(t; w)/dt 2 the acceleration vector at time t. The initial conditions associated with Eq. (3.1) are written as

Y(t 0 ) = y 0 , Ẏ(t 0 ) = y 1 ,
where y 0 and y 1 are two given vectors in R n y . For all w ∈ C w and u ∈ C u , [M(w, u)] is the mass matrix belonging to M + n y , g (y(t; w, u), dy(t; w, u)/dt; w, u) is the R n y vector representing the internal nonlinear forces at time t, and F(• ; w) is a given R n y -valued stochastic process modeling the external forces.

Construction of the ROB represented by matrix [V] and convergence criterion

Let ε comp be the given relative tolerance (for instance 10 -5 ) to perform the data compression of [y] ∈ M n y ,N d defined by Eq. (2.4). Let n max q be a given integer whose value is of the order of the rank of matrix [y] ∈ M n y ,N d with N d = n d ×n snps and such that n max q < n y . The singular value decomposition (SVD) of [y], restricted to the only calculation of the n max q largest singular values σ 1 ≥ . . . , ≥ σ n max q (ordered by descending values), allows for computing the n max q associated leftsingular vectors, represented by the matrix

[V max ] ∈ M n y ,n max q such that [V max ] T [V max ] = [I n max q ]. If (σ n max q /σ 1 ) 2 > ε comp , then n max
q has been chosen too small and the computation has to be restarted with a larger value of n max q . The optimal value of n q ≤ n max q is then calculated such that

err comp (n q ; N d ) ≤ ε comp < err comp (n q + 1; N d ) with err comp (n q ; N d ) = 1 - n q α=1 σ 2 α n max q α=1 σ 2 α .
(3.2)

Finally, matrix [V] ∈ M n y ,n q is made up of the first n q columns of [V max ] and is such that [V] T [V] = [I n q ]. For this fixed value of n q , matrix [V] ∈ M n y ,n q is such that [V] = arg min [v]∈M ny ,nq ,[v] T [v]=[I nq ] ∥ [y] -[v] [v] T [y] ∥ 2 F ,
in which ∥ • ∥ F is the Frobenius norm.

Reduced representation of X using PCA and convergence criterion

The reduced representation X (ν) defined by Eq. (2.8) of the R n x -valued random variable X = (Q, O id , W) is constructed as follows, using the training dataset D train (x) = {x 1 , . . . , x n d }. Beforehand the PCA computation, the training dataset D train (x) is scaled using the formulation presented in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] (at the end of the numerical procedure, a back scaling must be carried out). Let x j c = x j dx be the realization of

X with x = (1/n d ) n d j=1 x j d ∈ R n x . Let [x c ] = [x 1 c . . . x n d
c ] be the matrix in M n x ,n d . Since n x ≫ n d , the economy size SVD (thin SVD [START_REF] Golub | Matrix Computations, Second Edition[END_REF] of matrix [x c ] is carried out, which allows for obtaining the left-singular vectors represented by the matrix

[Φ c ] ∈ M n x ,n d such that [Φ c ] T [Φ c ] = [I n d ]. The corresponding singular values S 1 ≥ . . . ≥ S n d -1 > S n d = 0 are in decreasing order. For ν ≤ n d -1, the reduced representation X (ν) of X is given by Eq. (2.8), in which [Φ] ∈ M n x ,ν is made up of the first ν columns of matrix [Φ c ] and where [λ] is the diagonal matrix in M + ν such that [λ] αα = λ α = S 2 α /(n d -1)
. Note that the positive real numbers λ 1 ≥ . . . ≥ λ ν > 0 are the ν largest positive eigenvalues of the estimated covariance matrix [ C X ] of the covariance matrix [C X ] of X, performed using the training dataset. Therefore, [λ] and [Φ] depend on n d . As it can be seen, these eigenvalues and the associated eigenvectors are computed without computing [ C X ] because n x can be very large. It should also be noted that, if n d = n x and if ν < n d -1, then the sequence of random variables X (ν) is mean-square convergent to X when ν goes to n d -1, and if ν = n d -1 = n x -1, then Eq. (2.8) is not an approximation and corresponds to a change of basis. In general, for the high-dimension problems and small training datasets, n x is large and n d ≪ n x . Therefore, Eq. (2.8) corresponds to a reduced representation, which is an approximation whose accuracy depends on ν and n d and which is classically controlled as follows. For n d ≪ n x and for ν < n d -1, parameter ν is chosen such that

err PCA (ν ; n d ) = E{∥ X -X (ν) ∥ 2 } E{∥ X ∥ 2 } ≃ 1 - ν α=1 λ α tr[ C X ] ≤ ε PCA , ν < n d -1 , (3.3)
in which ε PCA is a given positive real number sufficiently small, where ∥ • ∥ is the usual Euclidean norm, and where E is the mathematical expectation operator. The trace, tr[ C X ], of matrix [ C X ] is calculated by estimating the diagonal entries of [ C X ] using the training dataset. Note that ν → err PCA (ν ; n d ) defined by Eq. (3.3) gives the relative error as a function of ν < n d -1 for a fixed value of n d .

"Projection" of the target on the model

We have to specify the meaning of Eq. (2.11) that we have defined as a "projection" of the target dataset onto the model in a context of incomplete data due to the partial observability. For all r ∈ {1, . . . , N r }, let η r targ ∈ R n o be the linear least squares solution of the equation

[A o ] η r targ = b r with b r = o r targ -o id ∈ R n o , which is the unique solution of the optimization problem η r targ = min η ∈ R ν ∥ [A o ] η -b r ∥ . ( 3
[A inv o ] = ([A o ] T [A o ]) -1 [A o ] T . For all η ∈ R ν , it is known that ∥ [A o ] η -b r ∥ ≥ ∥ [A o ] η r targ -b r ∥ in which η r targ is given by η r targ = [A inv o ] b r .
Therefore, the unique solution of the optimization problem defined by Eq. (3.4) is η r targ = [A inv o ] b r , (see Eq. (2.11)). The relative error in the "projection" of the target dataset on the model can be quantified, for fixed N r , by the mapping,

ν → err targ (ν; N r ) = N r r=1 ∥([I n o ] -[A o ] [A inv o ])(o r targ -o id )∥ 2 N r r=1 ∥o r targ -o id ∥ 2 , (3.5

Probabilistic learning constrained by the target dataset using a weak formulation of Fourier transform of probability measures

The construction of the learned dataset D learn (η ud ) = {η ℓ ud , ℓ = 1, . . . , N ud } of the R ν -valued random variable H ud is carried out using the Kullback-Leibler divergence minimum principle (KLDMP) based on the prior probability measure of H constructed with the training dataset D train (η) = {η j , j = 1, . . . , n d } and constrained by the target dataset D targ (η targ ) = {η r targ , r = 1, . . . , N r }. It should be noted that the constraints imposed for the KLDMP must be described by statistical moments, i.e., taking the form of a mathematical expectation. In the present case, the constraints are described by realizations that constitute the points of the target dataset. Therefore, to be able to impose the constraints using realizations, we use the extension [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF] of the constrained PLoM [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]. With this extension, the constraint, which is defined by the realizations of the target dataset, is transformed into a constraint expressed in the form of a mathematical expectation, which is necessary to implement the KLDMP. This transformation of constraint formulation can be done thanks to the use of a weak formulation of the Fourier transform of the probability measures. We refer the reader to [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF] for details of this method and the associated algorithms. However, we give below a very brief summary of the essential points of this method in order to facilitate the reading.

Prior probability measure of H. Let P H (dη) = p H (η) dη be the prior probability measure on R ν of H, whose probability density function η → p H (η) : R ν → R + is estimated by using the Gaussian kernel-density estimation (KDE) with the training dataset D train (η) = {η j , j = 1, . . . , n d }, involving 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] for which s SB = (4/(n d (2 + ν))) 1/(ν+4) is the Silverman bandwidth,

p H (η) = 1 n d n d j=1 1 ( √ 2π ŝ) ν exp - 1 2 ŝ2 ∥ ŝ s SB η j -η ∥ 2 , ∀η ∈ R ν , (3.6)
in which ŝ = s SB s 2 SB + (n d -1)/n d -1/2 . With such a modification, the normalization of H is preserved for any value of n d , that is to say,

E{H} = R ν η p H (η) dη = 1 2 ŝ2 η = 0 ν , (3.7
)

E{H ⊗ H} = R ν (η ⊗ η) p H (η) dη = ŝ2 [I ν ] + ŝ2 s 2 SB (n d -1) n d [ C H ] = [I ν ] , (3.8)
in which η ∈ R ν and [ C H ] ∈ M + ν are the estimates of the mean value and the covariance matrix of H, performed with D train (η). Theorem 3.1 in [START_REF] Soize | Probabilistic learning on manifolds[END_REF] proves that, for all η fixed in R ν , Eq. (3.6) is a consistent estimation of the sequence {p H } n d for n d → +∞.

Representation of the constraint defined by the target dataset D targ (η targ ). In [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF], it is proven that the constraint defined by the target dataset D targ (η targ ) = {η r targ , r = 1, . . . , N r } can be written as

E{h c (H)} = b c on R N r , (3.9)
in which h c (η) = (h c 1 (η), . . . , h c N r (η)) and b c = (b c 1 , . . . , b c N r ) are the vectors in R N r , which are written, for r ∈ {1, . . . , N r } and η ∈ R ν , as

h c r (η) = exp - 1 νs 2 ∥ η -η r targ ∥ 2 , b c r = 1 N r N r r ′ =1 exp - 1 νs 2 ∥ η r ′ targ -η r targ ∥ 2 , (3.10)
in which s = 4 (N r (2 + ν)) -1 1/(ν+4) .

Updated estimate using the Kullback-Leibler divergence minimum principle under the constraint.

C ad,p = η → p(η) : R ν → R + , R ν p(η) dη = 1 , R ν h c (η) p(η) dη = b c . (3.12)
It has been proven that there exists a unique solution to the optimization problem defined by Eqs. (3.11) and (3.12), which is reformulated using Lagrange multipliers to account for the constraints in the admissible set (refer to Theorem 3 in [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF] for the construction of the probability measure of H ud and the proof of its existence and uniqueness). To generate the learned dataset D learn (η ud ) = {η ℓ ud , ℓ = 1, . . . , N ud }, an MCMC algorithm is required [START_REF] Kaipio | [END_REF][START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization: Estimation[END_REF]. In this work, the MCMC generator used is a nonlinear Itô stochastic differential equation (ISDE) associated with the nonlinear stochastic dissipative Hamiltonian dynamical system proposed in [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle. applications to stochastic processes, random fields and random matrices[END_REF] and based on [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF]. This MCMC generator allows for the removal of the transient part to rapidly reach the stationary response associated with the invariant measure, for which the measure p H ud (η) dη is the marginal measure (refer to Theorem 4 in [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF] for the mathematical analysis of this MCMC generator). The ISDE is solved using the Störmer-Verlet algorithm, which provides an efficient and accurate MCMC algorithm. This algorithm can be easily parallelized to significantly reduce the elapsed time on a multicore computer. It should be noted that this MCMC generator can be considered as belonging to the class of Hamiltonian Monte Carlo methods [START_REF] Neal | MCMC using hamiltonian dynamics[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF], but it is distinct due to the presence of the dissipative term. Finally, the algorithm presented in [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF] can be readily extended to integrate the diffusion maps tool, which forms the basis of the PLoM algorithm ( [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF]). It is this extended PLoM algorithm that will be used in the application presented in Section 4.

Predictive statistical surrogate model

Once the learned dataset D learn (o ud , w ud ) = {(o ℓ ud , w ℓ ud ) ∈ R N o × R n w , ℓ = 1, . . . , N ud }, which are realizations of the updated R N o × R n w -valued random variable (O ud , W ud ), has been constructed (see Eq. (2.12)), the predictive statistical surrogate model defined by Eq. (2.13) can be implemented as explained at the end of Section 2.

(i) In order to prepare the developments that will be used to present the results of the application, we are going to reshape the updated R N o -random variable O ud into a M n obs ,n freq -valued random variable [O ud ] such that N o = n obs × n freq , which makes it possible to explain the dependence according to a parameter ω ∈ B = {ω 1 , . . . , ω n freq } ⊂ R + (for instance, ω will be the frequency and n obs will be the number of considered frequency-dependent observations). The Consequently, for any w 0 given in W ref , there is an index r 0 ∈ {1, . . . , N r } such that o r 0 targ ∈ R n o corresponds to the "measurement" performed on the dynamical system for the given value w 0 = w r 0 targ of the control parameter w.

(iii) For the given value w 0 of the control parameter, and for every observation defined by a fixed value of the observation index i in {1, . . . , n obs } and by a fixed value of the frequency index k in {1, . . . , n freq }, the statistical surrogate model consists (see Eq. (2.13)) in estimating the probability density function p Z(w 0 ) (z; w 0 ) of the real-valued random variable Z(w 0 ) defined as the conditional random variable

Z ud = [O ud ] ik given w 0 , p Z(w 0 ) (z; w 0 ) = p Z ud ,W ud (z, w 0 )/p W ud (w 0 ) , p W ud (w 0 ) = R p Z ud ,W ud (z, w 0 ) dz , (3.13
, i o ) 2 = n freq k=1 E [O ud (w 0 )] 2 i o k -(E{[O ud (w 0 )] i o k }) 2 n freq k=1 E{[O ud (w 0 )] i o k } 2 , (3.15) in which [O ud (w 0 )] i o k is the conditional random variable given W ud = w 0 such that, E{[O ud (w 0 )] α i o k } = R + o α k p [O ud (w 0 )] io k (o k ; w 0 , i o ) do k , α = 1, 2 , (3.16)
in which the conditional pdf (with support R + ) of [O ud (w 0 )] i o k is written as

p [O ud (w 0 )] io k (o k ; w 0 , i o ) = p [O ud (w 0 )] iok ,W ud (o k , w 0 ; i o )/p W ud (w 0 ) , ∀o k ∈ R + . (3.17)
The formulas, which allow for numerically estimating the conditional pdf p [O ud (w 0 )] io k (• ; w 0 , i o ) and E{[O ud (w 0 )] α i o k } using the N ud learned realizations [O ud (w 0 )] ℓ i o k and w ℓ ud for ℓ = 1, . . . , N ud , are given in Appendix A.2. For a fixed value of n d , for given w 0 and i o , the convergence analysis, with respect to n d , N r , and n MC is performed by studying the values of cv(n d , N r , n MC ; w 0 , i o ).

Distance of the statistical-surrogate-model predictions to the target

In order to evaluate the quality of the predictions carried out by the statistical surrogate model, we introduce the mean-square distance between the predictions and the target for each given value w 0 of the control parameter and for each observation i o ∈ {1, . . . , N o } for which there is a corresponding target i ∈ {1, . . . , n o } (as in Section 3.7). We thus define the real-valued random variable R i o by

R i o = ( n freq k=1 ([O ud ] i o k -[o r 0 targ ] ik ) 2 )/( n freq k=1 [o r 0 targ ] 2 ik ) , (3.18)
in which [o r 0 targ ] is the reshaping of o r 0 targ ∈ R n o . Note that R i o is not the statistical conditioning for given W ud = w 0 (we have just introduced the translation [o r 0 targ ] ik that depends on w 0 ). The N ud learned realizations of R i o are computed by

r ℓ i o = n freq k=1 ([o ℓ ud ] i o k -[o r 0 targ ] ik ) 2 n freq k=1 [o r 0 targ ] 2 ik , ℓ ∈ {1, . . . , N ud } . (3.19)
We now introduce the statistical conditioning R i o (w 0 ) of R i o given W ud = w 0 whose conditional probability density function r → p R io (w 0 ) (r; w 0 ) on R with support R + is written as p R io (w 0 ) (r; w 0 ) = p R io ,W ud (r, w 0 )/p W ud (w 0 ) , (3.20) in which p R io ,W ud is the joint probability density function of R i o and W ud . For given w 0 and i o , the mean-square distance between the prediction and the target is defined by d ms (w 0 , i o ) = (E{R i o (w 0 )}) 1/2 , and consequently,

d ms (w 0 , i o ) 2 = R + r p R i (w 0 ) (r; w 0 ) dr . (3.21)
The formulas, which allow for numerically estimating the pdf p R i (w 0 ) (• ; w 0 ) and the mean-square distance d ms (w 0 , i o ) using the N ud learned realizations r ℓ i and w ℓ ud for ℓ = 1, . . . , N ud , are given in Appendix A.3.

Application: nonlinear stochastic dynamics of a three-dimensional MEMS device

In this application, we partially reuse the nonlinear dynamical system presented in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF], but with different parameter values and the introduction of a vector-valued control parameter and a vector-valued random uncontrolled parameter. This system is particularly interesting because its response is highly sensitive to the nonlinearities considered, and there is a significant transfer of energy in the stochastic response outside the frequency band of the external excitation applied to the system. A 3D view of the system is shown in Figure 1-mid. We provide all the dimensions and mechanical constants values to allow for reproducibility of this application. Additionally, some values have been modified (if the system definition were based on [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF], it would result in numerous back-and-forth references and make the text less readable).

4.1. Nonlinear dynamical system (i) Definition of the nonlinear dynamical system. This 3D MEMS device considered is constituted of a mobile part made up of a square frame with a vertical beam attached to it. Its suspended part is constituted of a parallelepipedic The nominal mass density of the silicon material is 2 330 Kg/m 3 . A nonlinear elastic material is inserted between the aforementioned beams. Its constitutive equation corresponds to a cubic, elastic, restoring force with elastic constant k b whose nominal value is k b = 2 × 10 12 N/m. A zero x 1 -, x 2 -, and x 3 -displacement boundary conditions with respect to the reference frame Ox 1 x 2 x 3 are prescribed at the base of the mobile part of the device. In the reference frame, the following time-dependent, square integrable, and real-valued x 1 -acceleration is applied to the base,

Γ(t) = Γ 0 {sin(t(ω c + ∆ω c /2)) -sin(t(ω c -∆ω c /2))}/(π t) , ∀t ∈ [t 0 , T ] ,
where Γ 0 is the amplitude whose nominal value is Γ 0 = 120 m/s 2 , where ω c = 2 π×13×10 6 rad/s is the central angular frequency, and where ∆ω c = 2 π×10×10 6 rad/s is the angular frequency bandwidth. The energy of the excitation signal is mainly concentrated in the frequency band [-ω e , -ω min ]∪[ω min , ω e ], where ω min = ω c -∆ω c /2 = 2π×8×10 6 rad/s and ω e = ω c + ∆ω c /2 = 2π × 18 × 10 6 rad/s. At time t 0 , the device is at rest (its displacement and velocity fields are zero). In all analyses, the time-interval of analysis is [t 0 , T ] with t 0 = -2.7778 × 10 -5 s and T = 4.6403 × 10 -5 s. For this value of T , the device is returned to its zero equilibrium with a small relative error. 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 ] in which ω max has been identified as ω max = 2π × 72 × 10 6 rad/s (see figure 1-right. The time discretization is performed with a sampling time-step ∆t = π/ω max = 6.9444 × 10 -9 s. This yields n time = 10 682 time points in the time-interval [t 0 , T ]. For computing the Fourier transform of the observations (x 1 -acceleration at observed points) from the time-discretized responses, the sampling frequency step is set to ∆ω = 2π × 13 481 rad/s, yielding also n freq = 10 682 frequency points in the frequency band of analysis [-ω max , ω max ]. The results will be presented in the frequency domain for the frequency band B a = [0, ω a ], where ω a = 2 π × 70 × 10 6 rad/s.

(ii) Control parameter w. There are two control parameters and consequently, w = (w 1 , w 2 ) ∈ C w ⊂ R n w with n w = 2 and where the admissible set of w is C w = [0.5 , 1] × [0.5 , 0.9]. Component w 1 allows for controlling the amplitude Γ 0 of the x 1 -acceleration applied to the base, such that Γ 0 = w 1 Γ 0 . Component w 2 allows for controlling the elastic constant k b of the nonlinear elastic material, such that k b = w 2 k b .

(iii) Uncontrolled random parameter U. As we have explained in Section 2, the level of statistical fluctuations of U makes it possible to control the "diameter" of the generated family of computational models. We choose this uncontrolled random parameter U as the components of the random elasticity matrix [C elas ] with values in M + 6 whose mean value is [C sym elas ] ∈ M + 6 and whose level of statistical fluctuations is controlled by a hyperparameter δ train that will be fixed to 0.2. This value has been identified as large enough to generate a training dataset with sufficiently large statistical fluctuations to correctly update the probability measure of H ud under the constraint of the target dataset. The mean model is orthotropic and we could used a random orthotropic model as presented in [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF]. Nevertheless, for the reasons given above, it is better adapted to use for the statistical fluctuations an anisotropic model, the one introduced in [START_REF] Soize | Non Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF]

] yielding [C elas ] = [L sym elas ] T [G] [L sym elas ] in which [C sym elas ] = [L sym elas ] T [L sym elas ]
and where the probability measure, the hyperparameter δ train , and the random generator of random matrix [G] with values in M + 6 is given Page 103 of [START_REF] Soize | Uncertainty Quantification[END_REF], and is such that E{[G]} = [I 6 ]. This random generator allows for calculating n d independent realizations {u j , j = 1, . . . n d } of U.

(iv) Observation O(w) and identification observation O id (w). The observation and the identification observation are relative to the x 1 -acceleration in the frequency domain of n obs = 744 observed nodes of the finite element mesh (see Paragraph (v) below) distributed across various locations of the MEMS device, in particular on the boundaries of the vertical beams of the suspended and mobile parts.

(iv-1) Observation O(w). The frequency sampling is n obs freq = 5 193 frequency points ω obs k = ω 1 + (k -1) ∆ω for k ∈ {1, . . . , n obs freq } and with ω 1 = 2π × 13 481 rad/s and ∆ω = 2π × 13 481 rad/s. We defined the matrix-valued observation [O(w)] with values in M n obs ,n obs freq such that, for m ∈ {1, . . . , n obs } and k ∈ {1, . . . ,

n obs freq }, [O(w)] mk = log((ω obs k ) 2 | Y i m (ω obs k )|) with Y(ω obs k ) = T t 0 exp(-iω obs k t) Y(t) dt.
In this formula, i m is the dof number corresponding to the x 1 -displacement of the m-th observed node. Observation O(w), which is obtained by reshaping [O(w)], is with values in R N o with N o = n obs × n obs freq = 3 863 592 and is used for presenting the predictions obtained with the statistical surrogate model. Nevertheless, in order to keep the number of figures within a reasonable limit, two observations are selected among all 744 possible ones. These observations are relative to two observed nodes in the suspended part, the first one having coordinates (14, 6, 0) microns, located at the down left corner of the elastic beam of the suspended part and the second one having coordinates (7, 17, 0) microns, located at the down left corner of the massive suspended part (see Figure 1-left). These two observations are the R n obs freq -valued random variable denoted by O 1 obs (w) and O 2 obs (w). (iv-2) Identification observation O id (w). There are n id freq = 519 frequency points ω id k = ω 1 + (k -1) δω id for k ∈ {1, . . . , n id freq } and with δω id = 2π × 134 810 rad/s (there are 10 times less frequency points). We defined the matrixvalued identification observation [O id (w)] with values in M n obs ,n id freq such that, for m ∈ {1, . . . , n obs } and k ∈ {1, . . . , n id freq },

[O id (w)] mk = log((ω id k ) 2 | Y i m (ω id k )|) with Y(ω id k ) = T t 0 exp(-iω id k t) Y(t)

Computational model

The finite element mesh is shown in Figure 1-mid. There are 7 328 eight-nodes solid elements, 10 675 nodes, and n y = 32 025 dofs. There are 205 of these nodes, which belong to the base of the mobile part of the device: at each of these nodes, all displacement dofs are constrained to zero in the moving reference frame Ox 1 x 2 x 3 , due to the boundary conditions. Hence, there are 615 zero Dirichlet conditions. The governing equation for the relative displacement vector, Y(t; w), is of the type of Eq. (3.1) and written as

[M] Ÿ(t; w) + [D] Ẏ(t; w) + [K(U)] Y(t; w) + f NL (Y(t; w); w 2 ) = f(t; w 1 ) , t ∈ ]t 0 , T ] , (4.1)
in which the external force f is deterministic and depends on the control parameter w 1 via the amplitude Γ 0 = w 1 Γ 0 of the imposed acceleration. The internal nonlinear elastic force f NL , which are generated by the nonlinear elastic material inserted between the elastic beams, depends on the control parameter w 2 via the elastic constant k b = w 2 k b . The stiffness matrix [K(U)] of the suspended and mobile parts depends on the uncontrolled random parameter U that is the reshaping of the random elasticity matrix [C elas ] of the linear elastic silicon. The mass matrix [M] and the damping matrix [D] is independent of w and U. Damping matrix [D] is constructed using the global damping model described in Appendix A of [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF], which is applied to the nominal model and for which the damping rate is ξ d = 0.02. The nominal computational model is given by Eq. (4.1) in which U is replaced by the deterministic vector u corresponding to the reshaping of the nominal elasticity matrix [C sym elas ].

Target dataset and reference dataset for validation

The set W targ = {w

Constructing the training dataset

The prior probability model of the control parameter w is the R n w -valued random variable W for which its probability measure is uniform on the admissible set C w defined in Section 4.1(ii). Consequently, for given n d , the set W targ = {w 1 , . . . . The reduced representation of {Y(t; w), t ∈ J, w ∈ C w }, defined in Section 2.5(i), is constructed with the (t, w)-independent ROB in R n y represented by matrix [V] ∈ M n y ,n q that is computed as explained in Section 3.2. The computation of [V] is carried out with n snps = n time = 10 682 and n max q = 20. For n d ∈ {100, 200, 300}, the value of N d = n d × n snps are 1 068 200, 2 136 400, and 3 204 600, respectively, and Figure 2-left displays the graph of function n q → err comp (n q ; N d ) defined by Eq. (3.2). The three curves are almost superimposed. Using the criterion defined by Eq. (3.2) with err comp = 10 -5 yields n q = 11 that is independent of the three considered values of n d . This small value of n q shows the efficiency of the POD method for reducing data.

(ii) Constructing training dataset D train (x) and PCA-based reduced representation of X. The training dataset D train (x) is then constructed using Section 2.5(ii) and for j ∈ {1, . . . , n d }, its points x j ∈ R n x are defined by Eq. (2.6) with n x = 10 682 × 12 + 386 136 + 2 = 514 322 that is, in this case, independent of n d . The reduced representation of X using PCA is performed as explained in Section 3.3. Note that, for this application, n x is independent of N d , but the points of D train (x) depend on n d , and consequently, for a fixed value of the tolerance ε PCA , dimension ν is n d -dependent. For n d ∈ {100, 200, 300}, Figure 2-right displays the graph of function ν → err PCA (ν ; n d ) defined by Eq. (3.3). Choosing ε PCA = 0.01 yields for ν the values 44, 57, and 62, respectively. It should be noted that the chosen value of ε PCA = 10 -2 might appear to large. In fact, calculations were made with smaller values (up to 10 -6 ) and showed that, for this application, there was no significant impact on the predictions made by the statistical surrogate model. This choice makes it possible to reduce the computational costs.

Construction of the learned dataset D learn (η ud ) using PLoM constrained by the target dataset

The methodology presented in Sections 2.5(v) and 3.5 is used for constructing the learned dataset D learn (η ud ). In addition the PLoM method is used, i.e. the MCMC generator, which is based on an Itô equation derived from a dissipative Hamiltonian formulation, is projected on the diffusion maps basis (see the methodology in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF] and the last version of the used algorithm in [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF]). Let ε opt be the parameter of the kernel used for constructing the transition matrix defined in Section 5.2 of [START_REF] Soize | Probabilistic learning on manifolds[END_REF] and let m opt = ν + 1 be the dimension of the diffusion maps basis (see [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF]). Then for ε PCA = 10 -2 , for n d = 100, 200, 300, yielding ν = 44, 57, 62 respectively, we have ε opt = 106.5, 150.5, 192.5 respectively, and the distributions of the eigenvalues α → Λ α of the transition matrix are shown in Figure 4-left. As explained in Section 3.5, the optimization problem defined by Eqs. (3.11) and (3.12) is solve with the algorithm presented in Sections 6.3 to 6.6 of [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF]. For N r = 208, for the three considered values of n d , and for n MC = 10 000, Figure 4-right shows the graphs of the error function i → err(i) defined by Eq. (6.36) of [START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF], which allows for controlling the convergence of the iteration algorithm as a function of the iteration number i in order to estimate the optimal value of the vector-valued Lagrange multiplier. The results show in Figure 5 are obtained for ε PCA = 10 -2 , The construction of the large learned dataset D learn (o ud , w ud ) defined by Eq. (2.12) is carried out by using Sections 2.5(vi) and 3.6.

The convergence analysis with respect to N r is carried out using the coefficient of variation cv(w j targ , i; N r ) defined by Eq. (3.15), for the points of the reference dataset

D ref (o ud ) = {{O i ud,obs (w 1 ref ), . . . , O i ud,obs (w 4 ref )}, i = 1, 2} defined in Section 4.3.
For n d = 300, ε PCA = 10 -2 , and for n MC = 10 000, Figure 6-left for observation i = 1 and Figure 6-right for observation i = 2 display the graph of function N r → cv(w j ref , i; N r ) for the four values j = 1, 2, 3, 4 of the reference control parameters. In Figure 6, it can be seen that the coefficient of variation is decreasing with respect to N r when N r ≥ 60 for observation 1 and N r ≥ 100 for observation 2. Fluctuations in the convergence for the small values of N r can be observed. This is due to the fact that the coefficient of variation is sensitive to the points chosen to define the target dataset. Nevertheless, such a sensitivity is lost when N r increases because the distribution of points in the target dataset tends to become more homogeneous. Clearly, if the distribution of points is homogeneous for a given value of N r , then it would also be homogeneous for a large value of N r . Similarly, for N r = 208, the convergence analysis with respect to n d ∈ {100, 200, 300} for ε PCA = 10 -2 and n MC = 10 000 is carried out using the coefficient of variation cv(w j ref , i; n d ). Figure 7-left for observation i = 1 and Figure 7-right for observation i = 2 display the graph of function n d → cv(w j ref , i; n d ) for the four values j = 1, 2, 3, 4 of the reference control parameters. Figure 7 shows that, for n d ≥ 100, the coefficient of variation is decreasing with respect to n d (except for the first point w 1 ref ) for observation 1 while is decreasing for observation 2 for j equal to 2 and 4. It is clear that there is no reason for the convergence with respect to n d to be monotonous for n d ≥ 1. Since ε PCA = 10 -2 is fixed for all values of n d , the reduced dimension ν varies according to n d .

Distance between the surrogate-model prediction and the target

For n d = 300, ε PCA = 10 -2 , and n MC = 10 000, we have estimated the square d ms (w j ref , i; N r ) 2 of the distance defined by Eq. (3.21) between the surrogate-model prediction and the target, for the two observations {O i ud,obs (w j ref ), i = 1, 2} and for the four values {w j ref , j = 1, 2, 3, 4} of the reference control parameters. Figure 8-left for observation 1 and Figure 8-right for observation 2 displays the set of the values d ms (w j ref , i; N r ) 2 for N r ∈ {4, 16, 30, 63, 104, 208}. These figures show that, for N r = 208, the values belong to the interval [0.012 , 0.075] for observation 1 and to the interval [0.023 , 0.12] for obseravtion 2. These figures quantify the quality of the predictive statistical surrogate model and it can be seen that the prediction is good enough. For ε PCA = 10 -2 , n d = 300, n MC = 10 000, N ud = 3 000 000, and N r = 208, the graph of the conditional pdf r → p R i (w j ref ) (r; w j ref ), defined Eq. (3.20), of the random variable R i given W ud = w j ref , for the two observations {O i ud,obs (w j ref ), i = 1, 2} and for the four values {w j ref , j = 1, 2, 3, 4} of the reference control parameters, is displayed in Figure 9-left for observation 1 (with a bimodal pdf for j = 3) and in Figure 9-right for observation 2 (with a multimodal pdf for j = 4).

Prediction of the statistical surrogate model

The prediction of the statistical surrogate model is carried out by using Section 3.6(iii) with ε PCA = 10 -2 , n d = 300, n MC = 10 000, N ud = 3 000 000, and N r = 208. For these values of the main parameters of the algorithms, the projection error of the target dataset on the model stays sufficiently small and is converged with respect to ν (see Section 4.6), the convergence is reasonnably reached with respect to n d and N r (see Section 4.8), and the quality of the predictive surrogate model is good enough (see Section 4.9). The prediction of the statistical surrogate model consists in estimating the conditional confidence regions for the two observations {O i ud,obs (w j ref ), i = 1, 2} and for the four values {w j ref , j = 1, 2, 3, 4} of the reference control parameters. Taking into account the definition of these two observations given in Section 4.1(iv-1), we introduce the frequency dependent function ω → dB 1 j (ω) and ω → dB 2 j (ω) such that dB i j (ω k ) = {O i ud,obs (w j ref )} k for k = 1, . . . , n obs freq . For each observation i ∈ {1, 2} and for each value w j ref with j ∈ {1, 2, 3, 4} of the reference control parameter, the conditional confidence region is estimated with a probability level p c = 0.98. Figures 10 and11 display, for observations 1 and 2 respectively, the conditional confidence region of f → dB i j (2πf) and also the deterministic functions for the training and for the target corresponding to w = w j ref .

These figures show the the prediction is good. In particular, it can be seen in Figures 10-(c) and (d) and in Figures 11-(c) and -(d. The effects of the learning under the constraints defined by the target dataset are very visible: the target line is inside the conditional confidence region whereas the training line is outside this region.

Conclusion

We have presented a novel methodology that addresses the challenging problem of constructing predictive statistical surrogate models for parameterized uncertain nonlinear computational models. These difficulties arise from the chosen framework, which corresponds to real-world situations involving large stochastic computational models of complex systems encountered in engineering sciences. The main challenges are primarily attributed to the high dimensionality of the considered under-observed nonlinear uncertain computational model, its partial observability leading to incomplete data in the target dataset of the identification observations, the non-injectivity of the nonlinear mapping used to compute the identification observations from the stochastic responses of the computational model, and, above all, the limited availability of a small training dataset. The proposed approach is purely probabilistic. The surrogate model is not directly represented by an algebraic model but indirectly represented by a probability measure, with its generator facilitating the construction of a large learned dataset. We have provided an effective description of the statistical surrogate model, which allows for rapid computation of output statistics in an online context using only conditional statistics that explore the learned dataset. For any given value of the vector-valued control parameter, this statistical surrogate model provides the statistics of any observation of the stochastic computational model, including mean values, variances, probability density functions, and confidence intervals. The presented developments have been illustrated through a representative application that highlight all the aforementioned difficulties. The obtained results contribute to the validation of the proposed approach.

  y ∥ the operator norm. [x] T : transpose of matrix [x]. tr{[x]}: trace of the square matrix [x]. ∥ [x] ∥ F : Frobenius norm of matrix [x]. δ kk ′ : Kronecker's symbol. δ x 0 : Dirac measure at point x 0 . (iv) List of symbols. C: fourth-order tensor-valued random field. C u ⊂ R n u : support of the probability distribution of U. C w ⊂ R n w : admissible set of the control parameter w, which is the support of the probability distribution of W. D learn (η ud ): learned dataset for the updated random vector H ud . D learn (o ud , w ud ): learned dataset for random vector (O ud , W ud ). D targ (η targ ): set of the N r vectors, η r targ . D targ (o id ): target dataset constituted of the N r vectors o r targ , associated with the O id . D learn (x ud ): learned dataset for the updated random vector X ud . D train (η): training dataset for random vector H. D train (x): training dataset constituted of the n d realizations x j of X. D train (y, o id , w): training dataset constituted of the n d realizations (y j , o j id , w j ).

  .13) in which the conditional pdf p O ud | W ud (o|w) of O ud given W ud = w is written as p O ud | W ud (o|w) = p O ud , W ud (o, w)/p W ud (w).

. 4 )

 4 Since ν < n o , the rank of matrix [A o ] ∈ M n o ,ν is less than or equal to ν. There is always a unique left-pseudo inverse [A inv o ] ∈ M ν,n o such that on the one hand [A inv o ] [A o ] and [A o ] [A inv o ] are symmetric matrices and on the other hand such that [A o ] [A inv o ] [A o ] = [A o ] and [A o ] [A inv o ] [A o ] = [A o ]. The left pseudo-inverse can be computed using the SVD of [A o ] and then deducing [A inv o ] by inverting the non-zero singular values. If the rank of [A o ] is equal to ν, then all the singular values are positive and [A o ] T [A o ] is invertible. In such a case the left-pseudo inverse can be written as

  ) in which [A o ] [A invo ] is the orthogonal projector onto the range of [A o ] and consequently, where[I n o ] -[A o ] [A inv o ]is the orthogonal projector onto the null space of [A o ] T .

  corresponding learned realizations for [O ud ] are {[o ℓ ud ] ∈ M n obs ,n freq , ℓ = 1, . . . , N ud }. Since the statistical conditioning of O ud given W ud = w has been written as O ud (w), the statistical conditioning of [O ud ] given W ud = w is written as [O ud (w)]. (ii) In addition, the results will be presented for a given set W ref = {w 1 ref , . . . , w n ref ref } of values of the control parameter, whose n ref points are defined as the "reference dataset for validation" for which the reference responses are known what allows for validating the predictions. For this reason, this subset is chosen such that W ref ⊂ W targ but W ref W train .

  ) in which p Z ud ,W ud is the joint probability density function of Z ud and W ud . It can then deduce the estimate of the confidence interval [z -(w 0 ), z + (w 0 )] of Z(w 0 ) such that z -(w 0 ) :z -(w 0 ) -∞ p Z(w 0 ) (z; w 0 ) dz = 1p c , z + (w 0 ) : z + (w 0 ) -∞ p Z(w 0 ) (z; w 0 ) dz = p c ,(3.14)in which p c is a given probability level (for instance p c = 0.98). The formulas to numerically estimate the confidence interval using the N ud learned realizations z ℓ ud = [o ℓ ud ] ik and w ℓ ud for ℓ = 1, . . . , N ud , are given in Appendix A.1. 3.7. Mean-square convergence criterion of the surrogate-model prediction with respect to n d ,N r , and N ud In order to analyze the convergence of the surrogate-model prediction with respect to the number n d of points in training dataset D train (y, o id , w), the number N r of points in target dataset D targ (o id ) = {o 1 targ , . . . , o N r targ }, and the number N ud = n d × n MC of points in learned dataset D learn (o ud , w ud ), we introduce a mean-square convergence criterion of the surrogate-model prediction. For each given value w 0 of the control parameter and for each observation i o ∈ {1, . . . , N o }, for which the corresponding component of the target is i ∈ {1, . . . , n o } and which corresponds to the "measurement" superscript r = r o associated with w 0 , this criterion is defined as the square of a coefficient of variation, cv(w 0

Figure 1 :

 1 Figure 1: Left figure: 2D view of the MEMS device with lengths in 10 -6 m (from[START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF]); observation node of coordinates (14, 6, 0) microns: marked by symbol x; observation node of coordinates (7, 17, 0) microns: marked by symbol o. Mid figure: 3D view of the finite element mesh of the nonlinear computational stochastic dynamical model (from[START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-dimensional and high-dimensional nonlinear models[END_REF]) for which the reference frame is the Cartesian coordinate system Ox 1 x 2 x 3 attached to the mobile part, with origin O located at the bottom left corner of the device, axis Ox 1 horizontal and oriented positively from left to right, axis Ox 2 vertical and oriented positively from bottom to top, axis Ox 3 perpendicular to the plane Ox 1 x 2 oriented positively from bottom to top. Right figure: illustration of the linear response (dashed line) and of the nonlinear response (solid line) in the frequency domain and in dB of the nominal dynamical system for a given observation and for a given value of the control parameter.

  w n d } is made up of n d realizations drawn from this uniform distribution on C w . The considered values of n d are 100, 200, and 300. For each fixed value of n d , the training dataset D train (y, o id , w) is generated using the computational model defined by Eq. (3.20) and the methodology presented in Section 2.

4. 5 .

 5 Reduced representation of Y, training dataset D train (x), and reduced representation of X (i) SVD-based data compression for computing [V]

Figure 2 :

 2 Figure 2: Left figure: for n d ∈ {100, 200, 300} yielding N d ∈ {1 068 200, 2 136 400, 3 204 600}, graph of function n → err comp (n; N d ) (the three curves are almost superimposed). Right figure: for n d = 100 (thin line), 200 (mid thickness line), 300 (thick line), graph of function ν → err PCA (ν ; n d ).

4. 6 .

 6 Convergence analysis of the "projection" of the target on the model (i) For n d = 300, ε PCA = 10 -2 , and N r = 208, Figure3-left shows the graph of function ν → err targ (ν ; N r ) defined by Eq. (3.5) for which ν ∈ {1, . . . , 299}, and which allows for analyzing the convergence of the "projection" of the target on the model. The dimension ν of the reduced representation of X is constructed with the PCA of X and the truncation error is given by the function ν → err PCA (ν ; n d ) defined by Eq. (3.3) (see Figure2-right). It can be seen that for ν = 62 the PCA error is err PCA (62 ; 300) = 0.00997 and the corresponding value of the projection error of the target on the model is err targ (62 ; 208) = 0.1126. It should be noted that, for ν = 62, this projection error of the target on the model is small enough not to significantly impact the updating of the learned probability measure of H ud obtained by the probabilistic learning under the constraint defined by the projected target. Note also that for this case (n d = 300), the asymptotic value of the error, obtained for ν = 299 is err targ (299 ; 208) = 0.094. Since there is a partial observability (744 components of the 32 025 components of Y(• ; w) are only used to construct the identification observation), which induces incomplete data. This level, 0.094, of the projection error of the target could be reduced in considerably increasing the number n d of points in the training dataset, that would not be coherent with the framework of the proposed methodology for which a small training dataset is assumed (small value of n d ).(ii) In addition, for ε PCA fixed to the value 10 -2 and for N r = 208, Figure3-right displays the graph of function n d → err targ (ν(n d ) ; N r ) for n d ∈ {100, 200, 300} whose corresponding values of ν(n d ) are {44, 57, 62}. It can be seen that for these fixed values of N r and ε PCA , the projection error stays sufficiently small.

Figure 3 :

 3 Figure 3: For n d = 300, ε PCA = 10 -6 , and N r = 208, graph of function ν → err targ (ν ; N r ) (left figure). For ε PCA = 10 -2 and N r = 208, graph of function n d → err targ (ν(n d ) ; N r ) with n d ∈ {100, 200, 300} and ν(n d ) ∈ {44, 57, 62} (right figure).

Figure 4 :

 4 Figure 4: For ε PCA = 10 -2 and n d = 100 (thin line), 200 (medium line), and (thick line); graph of function α → Λ α (left figure). For N r = 208 and for n MC = 10 000, graph of function i → err(i) (right figure).

  n d = 300, n MC = 10 000, N ud = n d × n MC = 3 000 000 learned realizations, and N r = 208. Figure 5-left shows the cloud of the learned realizations {(η ℓ ud,1 , η ℓ ud,2 , η ℓ ud,3 ), ℓ = 1, . . . , N ud } (red points) and also the cloud of the n d = 300 points {(η j 1 , η j 2 , η j 3 ), j = 1, . . . , n d } (blue points) of the training dataset. It can be seen that the constraints defined by the target dataset has a strong effect; the prior probability measure of H is strongly modified by the constraints yielding the updated probability measure of H ud . As an illustration, Figure 5-right displays the pdf η 3 → p H 3 (η 3 ) estimated with the n d points of the training set and the pdf η 3 → p H ud,3 (η 3 ) estimated with the N ud points of the learned dataset (solid line). It can be seen the strong effect of the constraints defined by the target set. 4.8. Construction of the large learned dataset D learn (o ud , w ud ) and convergence analysis

Figure 5 :

 5 Figure 5: Learned dataset for ε PCA = 10 -2 , n d = 300, n MC = 10 000, N ud = 3 000 000, and N r = 208. Left figure: cloud of the learned realizations {(η ℓ ud,1 , η ℓ ud,2 , η ℓ ud,3 ), ℓ = 1, . . . , N ud } (red points) and cloud of the n d = 300 points {(η j 1 , η j 2 , η j 3 ), j = 1, . . . , n d } (blue points) of the training dataset. Right figure: pdf η 3 → p H 3 (η 3 ) estimated with the training set (dashed line) and pdf η 3 → p H ud,3 (η 3 ) estimated with the learned dataset (solid line).

Figure 6 :

 6 Figure 6: For n d = 300, ε PCA = 10 -2 , and n MC = 10 000, convergence analysis with respect to N r : graph of function N r → cv(w j ref , i; N r ) for observations O 1 ud,obs (w j ref ) (left figure) and O 2 ud,obs (w j ref ) (right figure), for the four values j = 1, 2, 3, 4 of the reference control parameters.

Figure 7 :

 7 Figure 7: For N r = 208, for ε PCA = 10 -2 , and for n MC = 10 000, convergence analysis with respect to n d : graph of function n d → cv(w j ref , i; n d ) for observations O 1 ud,obs (w j ref ) (left figure) and O 2 ud,obs (w j ref ) (right figure), for the four values j = 1, 2, 3, 4 of the reference control parameters.

Figure 8 :

 8 Figure 8: For n d = 300, ε PCA = 10 -2 , n MC = 10 000, and for N r ∈ {4, 16, 30, 63, 104, 208}, set of the values d ms (w j ref , i; N r ) 2 of the square of the distance between the surrogate-model prediction and the target for observations O 1 ud,obs (w j ref ) (left figure) and O 2 ud,obs (w j ref ) (right figure), for the four values j = 1, 2, 3, 4 of the reference control parameters.

Figure 9 :Figure 10 :

 910 Figure 9: For ε PCA = 10 -2 , n d = 300, n MC = 10 000, N ud = 3 000 000, and N r = 208, graph of the conditional pdf r → p R i (w j ref ) (r; w j ref ) for observation O 1 ud,obs (w j ref ) (left figure) and for observation O 2 ud,obs (w j ref ) (right figure), for the four values {w j ref , j = 1, 2, 3, 4} of the reference control parameters.

Figure 11 :

 11 Figure 11: Observation 2: conditional confidence region of f → dB 2 j (2πf) for W = w j ref (yellow region with orange edges), training (blue dashed line) and target (red solid thick line) corresponding to w = w j ref .

  (t; w), . . . , Y (m diff ) (t; w); t, w, U = F(t; w) , ∀ t ∈ J = {t 1 , . . . , t n time } ⊂ [t 0 , T ] , (2.1) in which N is a nonlinear operator and where F = {F(t; w), t ∈ J} with values in R n y is a given time-discretized stochastic process (time series) that depends on w. A deterministic initial condition is given at t = t 0 , where t n = t 0 + n ∆t for n = 1, . . . , n time . The random response is the time-discretized stochastic process Y(• ; w) = {Y(t; w), t ∈ J} with values in R n y , which depends on the control parameter w but also on U and F, and can be written as Y(t; w) = The SCM is assumed to be in high dimension, that is to say, n y is large. For all w ∈ C w and for a given prior probability measure of U and F, it is assumed that Y(• ; w) is the unique solution of the nonlinear SCM, and that Y(• ; w), Y(1) (• ; w), . . . , Y (m diff ) (• ; w) are second-order stochastic processes indexed by J with values in R n y . We are interested in the prediction of an observation that is represented by a w-dependent random variable O(w) with values in R N o , defined by a non-injective nonlinear mapping O of Y(. ; w) such that, O(w) = O(Y(. ; w)).

	y(t; w, U, F) in which y is a deterministic vector-valued function, (t, w, u, f) → y(t, w, u, f). Finally, {Y (m) (t; w) =
	d m Y(t; w)/dt For
	example, the random vector O(w) can be related to the logarithm of the modulus of the frequency-sampled Fourier
	transform in time of a subset of components {Y k (t; w), t ∈ J} of {Y(t; w), t ∈ J}. The objective of the predictions is not
	that of the prediction of Y(• ; w) but is that of the prediction of O(w). This remark is important within the framework
	of the proposed methodology.

m , 1 ≤ m ≤ m diff }.

  Construction of the large learned dataset D learn (o ud , w ud ). Eq. (2.8) allows for constructing the learned dataset D learn (x ud ) = {x ℓ N o of the updated w-independent random observation O ud are computed using nonlinear operator O such that o ℓ

ud , ℓ = 1, . . . , N ud }. From Eqs. (2.3) and (2.6) and for all ℓ ∈ {1, . . . , N ud }, we deduce the learned realizations y ℓ ud = (y ℓ ud (t 1 ), . . . , y ℓ ud (t n time )) ∈ R n time ×n y and w ℓ ud ∈ R n w . The learned realizations o ℓ ud ∈ R ud = O(y ℓ ud ). Finally, we obtain the learned dataset

  dt. Identification observation O id (w), which is obtained by reshaping [O id (w)], is with values in R n o with n o = n obs × n id freq = 386 136 and is used to construct the learned dataset.

  1 targ , . . . , w N r targ } of the considered values of the control parameters for generating the target dataset D targ (o id ) = {o 1 targ , . . . , o N r targ } defined in Section 2.2, results from an experimental plan of points chosen on a Cartesian grid of N r points in the admissible set C w defined in Section 4.1(ii). For each value of N r , the nodes of this 2Dgrid is approximatively uniform and lightly modified in order to contain four points w ref ∈ C w ⊂ R n w defined in the reference dataset W ref = {w 1 ref = (0.6, 0.6), w 2 ref = (0.6, 0.8), w 3 ref = (0.9, 0.6), w 4 ref = (0.9, 0.8)} of the reference control parameters. The considered values of N r are 4, 16, 30, 63, 104, and 208. It should be noted that it is normal to choose the set W ref ⊂ W targ because the reference responses that are used to validate the predictions made with the statistical surrogate model must correspond to values of the control parameter for which the response "measurements" were made. For each fixed value of N r , target dataset D targ (o id ) is generated using the methodology presented in Section 2.2 and the deterministic nominal computational model defined by Eq. (3.20) in which the nominal value U = u has been replaced by another value U = u ref corresponding to a modification of the elasticy constants in elasticity matrix [C sym elas ]. For the validation, the reference dataset is then defined by D ref (o) = {O i

	obs (w 1 ref ), . . . , O i obs (w 4 ref )}, i = 1, 2 ,

Acknowledgments

The authors acknowledge partial funding from DOE SciDAC FASTMath institute.

Appendix A. Formulas for conditional statistics

In this Appendix, we give the explicit formulas for computing the conditional statistics related to the predictive statistical surrogate model defined in Section 3.6(iii), related to the mean-square convergence criterion of the surrogate-model prediction defined in Section 3.7, and related to the distance of the statistical-surrogate-model predictions to the target defined in Section 3.8. We reuse the notations of these three sections, but for simplifying the writing subscript "ud" is removed. Consequently, for k fixed in {1, . . . , n freq }, the N ud learned realizations w ℓ ud of W ud , for ℓ = 1, . . . , N ud , are simply rewritten as w ℓ for W with ℓ = 1, . . . , N. Considering W = (W 1 , . . . , W n w ), w = (w 1 , . . . , w n w ), and w 0 = (w 0,1 , . . . , w 0,n w ), for j = 1, . . . , n w , let w j = 1 N N ℓ=1 w ℓ j and let σ j = 1

N-1 N ℓ=1 (w ℓ jw j ) 2 (if σ j = 0, then σ j is set to 1). Let W j = (W jw j )/σ j whose realizations are wℓ j = (w ℓ jw j )/σ j for ℓ = 1, . . . , N, and let w0, j = (w 0, jw j )/σ j .

Appendix A.1. Formulas for conditional statistics of the predictive statistical surrogate model

We reuse the notations of Section 3.6(iii), but simplifying the writing, for k fixed in {1, . . . , n freq }, the N ud learned realizations z ℓ ud = [o ℓ ud ] ik of Z ud for ℓ = 1, . . . , N ud , are simply rewritten as z ℓ for Z. The associated conditional random variable Z ud (w 0 ) given W = w 0 is also simply rewritten as

on R be the conditional cumulative distribution function of Z(w 0 ) given W = w 0 , which is defined by

We then have,

) is written as

y 0 e -t 2 dt. For each k fixed in {1, . . . , n freq }, the conditional upper bound z + (w 0 ) = Proba{Z ≤ z + (w 0 ) | w 0 } and the conditional lower bound z -(w 0 ) = Proba{Z ≤ z -(w 0 ) | w 0 } of the confidence interval [z -(w 0 ) , z + (w 0 )] of the random variable Z(w 0 ) given W = w 0 , for a given probability level p c ∈]0, 1[ (for instance, p c = 0.98) are computed by solving the equations,

. Formulas for mean-square convergence criterion of the surrogate-model prediction

We reuse the notations of Section 3.7 but again as in Appendix Appendix A.1, for simplifying the writing we define

Therefore, the N learned realizations of Q are denoted by q ℓ for ℓ = 1, . . . , N. The associated conditional random variable of Q given W = w 0 is written as

whose realizations are qℓ k = (q ℓ kq k )/σ q,k for ℓ = 1, . . . , N and k = 1, . . . , n freq . The square of the coefficient of variation cv(w 0 , i o ) 2 defined by Eq. (3.15) is written as

and where for α = 1 or 2, n+4) with n = n freq + n w and where γ 1 = 0 and γ 2 = s 2 .

Appendix A.3. Formulas for the distance of the statistical-surrogate-model predictions to the target

We reuse the notations of Section 3.8, removing subscript i o when no confusion is possible. .

Conflict of interest

The authors declare that he has no conflict of interest.