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

Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints

Keywords: Probabilistic learning, statistical inverse problem, Kullback-Leibler divergence, implicit constraints, stochastic homogenization, uncertainty quantification surely. BVP: boundary value problem. det: determinant. dof: degree of freedom. E: mathematical expectation operator. err: error function. ISDE: Itô stochastic differential equation. KDE: kernel density estimation. pdf: probability density function. PDE: partial differential equation

This paper deals with a probabilistic learning inference that allows for integrating data (target set) into predictive models for which the target set is constituted of statistical moments of the quantity of interest (QoI) and for which the training set is constituted of a small number of points. In a first part, we present a mathematical analysis of a general methodology that allows for estimating a posterior probability model for a stochastic boundary value problem from a prior probability model of a control parameter. The given targets are statistical moments of the QoI for which the underlying realizations are not available. Under these conditions, the Kullback-Leibler minimum principle is used for estimating the posterior probability measure. The constraints are represented by an implicit nonlinear mapping for which a statistical surrogate model is introduced. The MCMC generator and the necessary numerical elements are given to facilitate the implementation of the methodology. In a second part, an application is presented to illustrate the proposed theory and is also, as such, a contribution to the three-dimensional stochastic homogenization of heterogeneous linear elastic media in the case of a non-separation of the microscale and macroscale. Consequently, the macroscale is another mesoscale at larger scale with random effective/apparent elastic properties. For the construction of the posterior probability measure by using the probabilistic learning inference, in addition to the constraints defined by given statistical moments of the random effective/apparent elasticity tensor, the second-order moment of the random normalized residue of the stochastic partial differential equation has been added as a constraint. This constraint guarantees that the algorithm seeks to bring the statistical moments closer to their targets while preserving a small residue.

Introduction

(i) Objective of the paper. This paper deals with a very important question related to probabilistic learning algorithms that allow for integrating data (target set) into predictive models for which the training set is constituted of a small number of points and for which the target set is made up of statistical moments of some quantities of interest (QoI). The considered constraints are thus implicit. It is assumed that these statistical moments such as mean values, second-order moments, have been estimated with realizations (samples) that are no longer available. This situation occurs quite frequently when the data (the realizations/samples) have been lost, or deteriorated, or not commented on, or no longer exist, or are not accessible, but for which the values of some statistical moments have been published or are given in technical reports. As the experimental measurements remain very expensive, and many have been done in the past decades in many fields of physics and engineering sciences, it is very interesting to integrate into the models these experiments to which we know only statistical moments. Additional physics-based constraints can simultaneously considered. This is the case if we want the learning process to be also controlled by the model, that is to say, that the learned probability measure minimizes, for instance, the mean-square norm of the random normalized residue of the stochastic partial differential equation of the boundary value problem (BVP); this type of constraint is thus also implicit. We therefore consider a model driven by a stochastic BVP for which there are a vector-valued control parameter and a vector-valued quantity of interest (QoI) that is defined by a given transformation of the solution of the stochastic BVP. The graph of the functional dependence between the control parameter and the QoI defines a manifold that is not explicitly described. We consider the statistical inverse problem consisting in identifying the control parameter by giving the target set related to the QoI. A prior probabilistic model of the control parameter is introduced, which induces a random QoI. A training set, constituted of independent realizations belonging to the manifold, is constructed by solving the BVP. The framework of the developments presented is that for which the probabilistic learning inference is made using a computational model of which a single evaluation is assumed very high. Thus it is assumed that the training set consists of a small number of points as opposed to big data. The prior probability measure on the manifold is estimated with the points of the training set. The probabilistic learning inference consists in estimating the posterior probability model on the manifold, which is closest to the prior probability measure while satisfying the constraints defined by an implicit function h c for which an algebraic representation is not available.

(ii) Ways for addressing the integration of data into the predictive models and methodology proposed in the paper. The literature related to the field treated is so vast that it would not be reasonable to attempt a review of the state of the art. We will therefore limit ourselves to mentioning the main statistical methods, which make it possible to integrate data into the predictive models in order to allow the positioning of the method that is proposed in this paper.

(ii-1) The Bayesian method provides a rational and efficient framework for integrating data into the predictive models. This statistical tool allows for estimating the posterior probability measure of the control parameter giving its prior probability measure and a target set of realizations for the QoI. Consequently, the Bayesian inference is certainly one of the most popular statistical tools to solve the statistical inverse problems in the context of parametric or nonparametric statistics (see, for instance, [START_REF] Bernardo | Bayesian Theory[END_REF][START_REF] Kennedy | Bayesian calibration of computer models[END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization: Estimation[END_REF][START_REF] Congdon | Bayesian Statistical Modelling[END_REF][START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Gentle | Computational Statistics[END_REF][START_REF] Marin | Approximate Bayesian computational methods[END_REF][START_REF] Givens | Computational Statistics[END_REF][START_REF] Scott | Bayes and big data: The consensus monte carlo algorithm[END_REF][START_REF] Ghanem | Handbook of Uncertainty Quantification[END_REF][START_REF] Soize | Sampling of Bayesian posteriors with a non-Gaussian probabilistic learning on manifolds from a small dataset[END_REF] for general aspects, [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF][START_REF] Marzouk | Stochastic spectral methods for efficient Bayesian solution of inverse problems[END_REF][START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF][START_REF] Soize | A computational inverse method for identification of non-Gaussian random fields using the Bayesian approach in very high dimension[END_REF][START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF][START_REF] Matthies | Inverse problems in a Bayesian setting[END_REF][START_REF] Bilionis | Bayesian uncertainty propagation using Gaussian processes[END_REF][START_REF] Dashti | The Bayesian approach to inverse problems[END_REF][START_REF] Arnst | Itô-SDE MCMC method for Bayesian characterization of errors associated with data limitations in stochastic expansion methods for 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] for specific aspects related to statistical inverse problems, [START_REF] Shen | Variational Markov chain Monte Carlo for Bayesian smoothing of non-linear diffusions[END_REF][START_REF] Depraetere | A comparison of variational approximations for fast inference in mixed logit models[END_REF] for variational Bayesian methods, [START_REF] Golightly | Bayesian sequential inference for nonlinear multivariate diffusions[END_REF] for Bayesian sequential inference, [START_REF] Fearnhead | Exact and efficient Bayesian inference for multiple changepoint problems[END_REF] for Bayesian inference for changepoint problems or [START_REF] Owhadi | Bayesian numerical homogenization[END_REF] for Bayesian numerical homogenization). Bayesian inferences have also been considered in the framework of machine learning [START_REF] Neil | Inference in hybrid Bayesian networks using dynamic discretization[END_REF][START_REF] Sambasivan | A Bayesian perspective of statistical machine learning for big data[END_REF] and probabilistic learning for small data sets and in high dimension [START_REF] Soize | Sampling of Bayesian posteriors with a non-Gaussian probabilistic learning on manifolds from a small dataset[END_REF]. Bayesian inference is therefore a powerful statistical tool for integrating raw data but requires that targets be given in the form of realizations, which is not the hypothesis introduced in this paper. Note also that Bayesian inference can remains tricky to use [START_REF] Owhadi | On the brittleness of Bayesian inference[END_REF], in particular for the high dimension.

The maximum-likelihood method (see [START_REF] Congdon | Bayesian Statistical Modelling[END_REF][START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization: Estimation[END_REF]) is also used for estimating the hyperparameters of probability measure from a target set of realizations (see also [START_REF] Murshudov | Refinement of macromolecular structures by the maximum-likelihood method[END_REF][START_REF] Zhang | Optimization of identification structure parameters based on recursive maximum likelihood iteration[END_REF][START_REF] Vinayak | Maximum likelihood estimation for learning populations of parameters[END_REF][START_REF] Sur | A modern maximum-likelihood theory for high-dimensional logistic regression[END_REF][START_REF] Lio | Uncertain maximum likelihood estimation with application to uncertain regression analysis[END_REF] for improvements of the basis methodology). This statistical tool also allows for solving statistical inverse problem in high dimension for non-Gaussian random field (see [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Desceliers | Identification of chaos representations of elastic properties of random media using experimental vibration tests[END_REF]). Similarly to the Bayesian inference, the maximum-likelihood method is also a very effective statistical method for estimating hyperparameters of probability measures, but requires that targets be given in the form of realizations.

All the parameterized functional representations of random quantities in finite or in infinite dimension, such as the spectral approaches, allow for integrating data into predictive models by using the Bayesian inference, the maximum likelihood, or the least-square approach for identifying the parameters of the representations. In this context, it should be noted that a general parameterized representation of non-Gaussian second-order random quantities (vector-valued random variables, stochastic processes, and random fields) is given by the polynomial chaos expansions (see [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach[END_REF][START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF][START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[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] 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] Mignolet | Compressed principal component analysis of non-Gaussian vectors[END_REF][START_REF] Luthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF]).

(ii-2) In Artificial Intelligence (see [START_REF] Russel | Artifical Intelligence, A Modern Approach[END_REF] for a remarkable synthesis of this field), Machine learning has taken a considerable place in all fields of sciences and engineering and today provides an essential and very effective tool for many applications. This tool is mainly based on the use of Artificial Neural Networks (ANNs), which allow building an algebraic representation from data for both unsupervised and supervised cases. By construction, this method makes it possible to integrate data. In this context, physics-informed machine learning has seen a growing interest in recent years. Numerous works have been carried out in computational science and engineering [START_REF] Saha | Hierarchical deep learning neural network (HiDeNN): An artificial intelligence (AI) framework for computational science and engineering[END_REF], for instance in the field of computational fluid dynamics [START_REF] Wang | Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data[END_REF][START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF][START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework[END_REF][START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF], for projection-based model reduction [START_REF] Swischuk | Projection-based model reduction: Formulations for physics-based machine learning[END_REF], for hierarchical materials and their performance using multiresolution analysis and clustering discretization [START_REF] Li | Clustering discretization methods for generation of material performance databases in machine learning and design optimization[END_REF][START_REF] Yu | Multiresolution clustering analysis for efficient modeling of hierarchical material systems[END_REF]. Likewise, physics-informed machine learning has also been used in the context of linear [START_REF] Raissi | Machine learning of linear differential equations using Gaussian processes[END_REF], quasilinear [START_REF] Sirignano | DGM: A deep learning algorithm for solving partial differential equations[END_REF], and nonlinear partial differential equations [START_REF] Raissi | Hidden physics models: Machine learning of nonlinear partial differential equations[END_REF][START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Kats | A physics-informed machine learning method for predicting grain structure characteristics in directed energy deposition[END_REF]. It should be noted that the use of ANNs requires the availability of big data in order to be able to identify the network parameters by learning algorithms such as the deep learning [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF] and therefore, to obtain a predictive ANN. As we have previously explained, for the probabilistic learning that is presented in this paper, it is assumed that the training set is made up of a small number of points and consequently, these types of approaches are, in general, not available.

The statistical and probabilistic learning approaches deals with the problem of inferring a function based on data. In particular for supervised learning, it consists in training from a set of data consisting of input-output pairs. These methods have also been extensively developed (see for instance, [START_REF] Schölkopf | Kernel principal component analysis[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] Aggarwal | Mining Text Data[END_REF][START_REF] Dalalyan | Sparse regression learning by aggregation and Langevin Monte-Carlo[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] Balcan | Statistical active learning algorithms[END_REF][START_REF] James | An Introduction to Statistical Learning[END_REF][START_REF] Dong | Knowledge vault: A web-scale approach to probabilistic knowledge fusion[END_REF][START_REF] Ghahramani | Probabilistic machine learning and artificial intelligence[END_REF][START_REF] Taylor | Statistical learning and selective inference[END_REF]) and play an increasingly important role in computational physics and engineering science, in particular for problems using large scale computational models. In this framework, statistical methods for supervised learning have been developed in the form of surrogate models from which approximations of the expensive functions can be easily evaluated [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF][START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Eldred | Design under uncertainty employing stochastic expansion methods[END_REF][START_REF] Yao | Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles[END_REF][START_REF] Homem-De Mello | Monte Carlo sampling-based methods for stochastic optimization[END_REF][START_REF] Papadopoulos | A neural network-based surrogate model for carbon nanotubes with geometric nonlinearities[END_REF][START_REF] Giovanis | Data-driven surrogates for high dimensional models using Gaussian process regression on the Grassmann manifold[END_REF]. The probabilistic learning is the framework of this paper as we have explained in paragraph (i).

(ii-3) In the context of probabilistic learning and Information Theory, the Kullback-Leibler minimum principle [START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF][START_REF] Cover | Elements of Information Theory[END_REF] is a method that allows for estimating the posterior probability measure on the manifold given its prior probability measure and the constraints related to the statistical moments for which targets are given. This method is the one that is used in the methodology presented in this paper. It should be noted that this principle has widely been used in many fields (see for instance, [START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF][START_REF] Vasconcelos | The Kullback-Leibler kernel as a framework for discriminant and localized representations for visual recognition[END_REF][START_REF] Zhang | Local Gabor binary patterns based on Kullback-Leibler divergence for partially occluded face recognition[END_REF][START_REF] Cappé | Kullback-Leibler upper confidence bounds for optimal sequential allocation[END_REF][START_REF] Saleem | Low rank sparse decomposition model based speech enhancement using gammatone filterbank and Kullback-Leibler divergence[END_REF]) in particular, for reinforcement learning [START_REF] Filippi | Optimism in reinforcement learning and Kullback-Leibler divergence[END_REF] and for probabilistic learning [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]). The posterior probability measure is represented by an algebraic expression of the prior probability measure and of a vector-valued Lagrange multiplier λ associated with the vector-valued function h c . The optimal value λ sol of the Lagrange multipliers is obtained as the limit of a sequence {λ i } i of Lagrange multipliers allowing for constructing a sequence of probability measures whose limit, for λ = λ sol is the searched posterior probability measure. In this paper, we will call "constrained learned set", the set of realizations that results from the learning process in taking into account the constraints. As the framework of the developments presented is that of the high dimension with a small training set, for each value λ i of λ, the constrained learned set must be generated with a MCMC algorithm. The MCMC generator is based on a nonlinear Itô stochastic differential equation (ISDE) associated with a nonlinear stochastic dissipative Hamiltonian dynamical system. The presence of a dissipative term makes it possible to delete the transient part in order to quickly reach the stationary response associated with the invariant measure. The evaluation of the drift vector of the ISDE requires to evaluate the gradient of function h c a large number of times. In the framework of this work, there is not an available algebraic expression of this gradient. Furthermore, it is not possible to do a direct numerical calculation of it, taking into account the high dimension (function h c is implicit and for each evaluation of it, the BVP must be used). For instance, components of h c can be related to the norm of the random normalized residue of the partial differential equation (including the boundary conditions) of the BVP and the statistical moments of random observations. The construction of a surrogate model of implicit function h c by using a deterministic approach, such as the meshless methods [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF][START_REF] Duarte | H-p clouds, an h-p meshless method, Numerical Methods for Partial Differential Equations[END_REF][START_REF] Breitkopf | Explicit form and efficient computation of MLS shape functions and their derivatives[END_REF][START_REF] Rassineux | Surface remeshing by local Hermite diffuse interpolation[END_REF][START_REF] Zhang | Meshless methods based on collocation with radial basis functions[END_REF], is not adapted taking into account a possible high dimension of the space on which h c is defined. To circumvent this difficulty, we generalize the approach proposed in [START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF], which consists in constructing a statistical surrogate model ĥN of h c , depending on the number N of points generated in the constrained learned set, for which its gradient has an explicit algebraic representation.

(iii) Organization and novelties of the paper. First of all, let us point out that a similar problem has been tackled in [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF] devoted to take into account constraints in the PLoM (probabilistic learning on manifolds) method [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 with partition[END_REF]. However, in [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF], function h c is explicit. In this work, h c is implicit (that is to say the use of the BVP is required for evaluating its value in any given point). The presented methodology is novel and general (note that PLoM is not used, but could be implemented if necessary, but this would be prejudicial to the clarity of the developments). In addition, a mathematical analysis of the methodology is presented. This analysis is necessary because, due to the use of the statistical surrogate model ĥN of h c , the constrained learned realizations of the posterior probability measure on the manifold are generated with the MCMC generator that is the limit of a sequence of MCMC generators, each generator of the sequence depending on ĥN instead of h c . This means that convergence properties with respect to N must be studied. The organization of the paper is the following.

In Section 2, we set the problem and summarize the methodology of the probabilistic learning inference that is proposed. The used hypotheses are discussed in order to expose the difficulties involved by the choice of a general framework for the developments: small number of points in the training set, targets defined by statistical moments of the observations, implicit description of function h c and the necessity to construct a statistical surrogate model. Section 3 deals with the mathematical analysis of the proposed methodology. Lemma 1 proves the convexity of the admissible set C ad,λ ⊂ R n c of the Lagrange multipliers λ and the integrability properties related to the sequence of posterior probability measures indexed by λ. Proposition 1 gives the construction of λ sol as the unique solution of a convex optimization problem posed in C ad,λ . After building the statistical surrogate model ĥN of h c (Definition 3), Proposition 2 proves the convergence of the sequence { ĥN } N towards h c and of its gradient for λ fixed in C ad,λ . Proposition 3 proves the existence and uniqueness of the invariant measure of the ISDE for each fixed value of λ and gives the generator of the constrained learned set depending on λ. Proposition 4 proves the mean-square convergence, as N goes to infinity, of the stationary stochastic solution of the ISDE calculated with ĥN towards the one calculated with h c . After having given the construction of the iteration algorithm for calculating the sequence {λ i } i that converges to λ sol , Proposition 5 gives the rate of convergence of the sequence of posterior probability measures indexed by λ.

In Section 4, a few numerical elements are given for implementing the methodology. The Störmer-Verlet scheme that allows for solving the ISDE is detailed and the explicit expression of the gradient of ĥN is given. We present the algorithm for calculating λ sol and for generating the constrained learned set D H c . Section 5 is devoted to stochastic homogenization that has given rise to a large number of works (see for instance [START_REF] Papanicolaou | Boundary Value Problems with Rapidly Oscillating Random Coefficients[END_REF][START_REF] Torquato | Microstructure of two-phase random media. v. the n-point matrix probability functions for impenetrable spheres[END_REF][START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Sab | On the homogenization and the simulation of random materials[END_REF][START_REF] Andrews | Stochastic homogenization of elliptic boundary-value problems with Lp-data[END_REF][START_REF] Ostoja-Starzewski | Random field models of heterogeneous materials[END_REF][START_REF] Jikov | Homogenization of Differential Operators and Integral Functionals[END_REF]) and for which the analysis of the representative volume element (RVE) size has received a particular attention (see [START_REF] Drugan | A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites[END_REF][START_REF] Ren | Effects of grain sizes, shapes, and distribution on minimum sizes of representative volume elements of cubic polycrystals[END_REF][START_REF] Sab | Periodization of random media and representative volume element size for linear composites[END_REF][START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF][START_REF] Ostoja-Starzewski | Comparisons of the size of the representative volume element in elastic, plastic, thermoelastic, and permeable random microstructures[END_REF][START_REF] Yin | Statistical volume element method for predicting microstructure-constitutive property relations[END_REF][START_REF] Soize | Tensor-valued random fields for meso-scale stochastic model of anisotropic elastic microstructure and probabilistic analysis of representative volume element size[END_REF][START_REF] Zhang | Robust multiscale identification of apparent elastic properties at mesoscale for random heterogeneous materials with multiscale field measurements[END_REF][START_REF] Soize | Computational stochastic homogenization of heterogeneous media from an elasticity random field having an uncertain spectral measure[END_REF]). In the context of the stochastic homogenization of the three-dimensional elastic microstructures of materials, the effective quantities are sought at the macroscopic scale by making a spatial average over a representative volume element (RVE). The effective quantities are those of the macroscopic scale when the dimensions of the RVE in the three directions are taken sufficiently large so that the statistical fluctuations of the effective quantities are negligible compared to the statistical mean value. In this paper, we present an application related to the stochastic homogenization of a random linear elastic medium at mesoscale posed on the domain Ω ⊂ R 3 of the microstructure, which is not a RVE, which means that there is no scale separation between the mesoscale and the "macroscale". Since there is no scale separation, the "macroscale" is in fact another mesoscale at a larger scale than the initial mesoscale, which is described by random effective/apparent quantities. Consequently, the effective/apparent elasticity tensor at macroscale is not deterministic and has statistical fluctuations. This case is obtained when the spatial correlation lengths of the random apparent elasticity field at mesoscale are not sufficiently small with respect to the size of domain Ω. For instance, such a situation is encountered when the size of an experimental specimen, which is measured with the objective to perform an inverse identification of its effective properties, is not sufficiently large compared to the size of the heterogeneities. It is then interesting to identify the probabilistic model of the random apparent elasticity field, without having a scale separation. One can then estimate the effective properties using the identified random apparent elasticity field with a stochastic computational model for which domain Ω is an RVE, that is to say, has a size that is largest than the one of the specimen. In order to analyze the scale separation, we will consider several values of the spatial correlation lengths to cover three cases: SC1: partial separation (separated for two directions but not in the third one). SC2: not separated in the three directions. SC3: strongly separated in the three directions.

Notations

x, η: lower-case Latin or Greek letters are deterministic real variables.

x, η: boldface lower-case Latin or Greek letters are deterministic vectors. X: upper-case Latin letters are real-valued random variables. X: boldface upper-case Latin letters are vector-valued random variables.

[x]: lower-case Latin letters between brackets are deterministic matrices.

Framework of the developments presented in the paper

For instance, we consider a stochastic elliptic BVP on an open bounded domain Ω ⊂ R d with d ≥ 1, whose partial differential equation (PDE) is written as N(Y, G, W) = 0 a.s. The unknown is the non-Gaussian vector-valued field {Y(ξ), ξ ∈ Ω} defined on (Θ, T , P) and which satisfies the boundary conditions. The coefficients of the stochastic elliptic operator depend on a non-Gaussian second-order vector-valued random field G defined on (Θ, T , P) and on a random vector-valued control parameter W also defined on (Θ, T , P). It is assumed that the weak formulation of this stochastic BVP admits a unique strong stochastic solution Y = f (G, W), which is a second-order random field. The observation (quantity of interest) is, for instance, a second-order vector-valued random variable

Q = O(Y, G, W) in which O is a given measurable mapping.
The problem under consideration belongs to the class of the statistical inverse problems. A prior probability model of {G, W} is given and we are interested in estimating a posterior model {G c , W c } of {G, W} in order that some statistical moments of the posterior observations

Q c = O(Y c , G c , W c ) with Y c = f (G c , W c
), are equal to some given targets (the superscript "c" is introduced to designate the solution with the constraints, which corresponds to the posterior model). The statistical moments of Q c are globally written as E{M c (Q c )} = b c in which b c ∈ R n c is the target, q → M c (q) is a given measurable mapping, and E is the mathematical expectation operator. As explained in Section 1, the realizations that have allowed b c to be estimated are not available. Only b c is known (it can be the case when one or several components of b c have been estimated with experimental realizations that are not available). However, as we will see in Section 5, a component of the vector of statistical moments can also be related to the random normalized residue of the PDE of the stochastic boundary value problem.

j d = O(y j d , g j d , w j d ).
The training set is then made up of a small number N d of points x j d = {y j d , g j d , w j d } for j = 1, . . . , N d , which are N d independent realizations of X = {Y, G, W}. Note that a strong hypothesis in the present work is that the BVP can only be solved a small number of times. This means that the training set is a small data set (as opposed to a big data set). Consequently, for constructing the posterior probability measure on the manifold, a learning tool must be used for generating the constrained learned realizations of X without solving the BVP, but using only the training set.

(ii) Statistical inverse problem. The available information consists of the targets of statistical moments related to observation Q c and represented by b c . As we have explained in Section 1, the realizations that have been used for estimating b c are not available. This is the second strong hypothesis used in this work. Consequently the classical statistical tools such as the Bayesian inference or the maximum likelihood method cannot easily be used for solving the statistical inverse problem under consideration.

(iii) Finite reduced-order representation. The second-order random variable X = {Y, G, W}, defined on (Θ, T , P), is assumed to be with values in a real Hilbert space X equipped with the inner product X , X X and its associated norm X X = X , X 1/2 X . Consequently, X belongs the Hilbert space L 2 (Θ, X) of the equivalent class of all the second-order random variables with values in X, equipped with the inner product X , X = E{ X , X X } for which the square of the associated norm is

||| X ||| 2 = E{ X 2 X } = Θ X(θ) 2 X d P(θ).
Since the problem is in infinite dimension, in order to implement the probabilistic learning inference that we propose, we need to introduce a finite representation X (ν) of dimension ν of random variable X in L 2 (Θ, X). Assuming that the covariance operator is a Hilbert-Schmidt [START_REF] Gelfand | Les Distributions. Tome 4. Application de l'Analyse Harmonique[END_REF], symmetric, positive operator in X, X (ν) can be represented using the truncated Karhunen-Loève expansion [START_REF] Karhunen | Ueber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Functions aleatoires du second ordre[END_REF] of X,

X (ν) = x + ν α=1 √ κ α ϕ α H α , (2.1)
in which the eigenvalues of the covariance operator are κ 1 ≥ . . . ≥ κ ν ≥ . . . = 0 with +∞ α=1 κ 2 α < +∞, where the family of the eigenfunctions {ϕ α } α is a Hilbertian basis of X, where x = E{X}, and where H = (H 1 , . . . , H ν ) is a second-order, centered, R ν -valued random variable whose covariance matrix is the identity matrix [I ν ] in M ν . It should be noted that the covariance operator of X is not positive definite but only positive (its kernel is not reduced to zero) and therefore, there is a zero eigenvalue with a finite multiplicity, which is not taken into account in the truncated representation of X, defined by Eq. (2.1), in which κ ν > 0. The following heuristic argument can be given for explaining why the covariance operator of X is only positive. Since the deterministic mapping N relates the random quantities Y, G, and W, the "functional components" of X are not "independent", and so the covariance operator is not positive definite. For α ∈ {1, . . . , ν}, the component H α is written as

H α = κ -1/2 α X -x , ϕ α X . The training set D d related to H is made up of the N d independent realizations {η j d , j = 1, . . . , N d } such that η j d = κ -1/2 α x j d -x , ϕ α X .
If the kernel of the covariance operator was explicitly known, then ν would be chosen in order that ||| X -X (ν) ||| ≤ ε for a sufficiently small value of ε. In this paper, we assume that the kernel of the covariance operator is unknown. Therefore, we can only obtain an approximation of the covariance operator using an empirical estimator built with the N d points {x j d , j = 1, . . . , N d } (for instance, see Section 5.6). Under these conditions the largest value of ν will be N d -1 and the discretization of Eq. (2.1) will simply correspond to a normalization of the points that constitute the training set D d .

Formulation using the Kullback-Leibler divergence minimum principle

Taking into account Section 2.2-(ii), we use the Kullback-Leibler divergence minimum principle [START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF][START_REF] Cover | Elements of Information Theory[END_REF] for estimating the posterior probability measure

P H c (dη) = p H c (η) dη on R ν of the R ν -valued random variable H c = (H c 1 , . . . , H c ν )
. This estimation of P H c is performed (1) using the prior probability measure P H (dη) = p H (η) dη on R ν in which p H is estimated with the Gaussian KDE method using the points {η 1 d , . . . , η N d d } of the training set D d , and (2) using the constraint defined by the given statistical moments,

E{h c (H c )} = b c ∈ R n c . (2.2) In Eq. (2.2), η → h c (η) is a function from R ν into R n c with n c > 1, such that h c (H c ) is a random variable equal to M c (Q c ) that
p H c = arg min p∈C ad,p R ν p(η) log p(η) p H (η) dη , (2.3)
in which the admissible set C ad,p is defined by

C ad,p = η → p(η) : R ν → R + , R ν p(η) dη = 1 , R ν h c (η) p(η) dη = b c . (2.4)
It should be noted that, since we are interested in divergence (from p H (η) dη) of probability measure p(η) dη that satisfies the constraints expressed in C ad,p , we are in fact interested in the divergence and therefore, the symmetry condition is irrelevant in this case.

Methodology used for solving the optimization problem and MCMC generator

The constraints, which are defined in admissible set C ad,p , are taken into account by introducing the Lagrange multipliers λ 0 -1 with λ 0 ∈ R + , which is associated with the normalization condition, and λ ∈ C ad,λ ⊂ R n c , which is associated with the moments constraints. The admissible set C ad,λ of λ is a subset of R n c , which is completely defined by Definition 2 in Section 3. The Lagrange multiplier λ 0 is eliminated as a function of λ. In Eq. (2.3), the posterior pdf p H c is constructed as the limit of a sequence {p H λ } λ of probability density functions of a R ν -valued random variable H λ = (H λ,1 , . . . , H λ,ν ) that depends on λ. The construction of {p H λ } λ requires to generate with a MCMC algorithm a constrained learned set

D H λ = {η 1 λ , . . . η N λ } constituted of N N d independent realizations {η λ , = 1, . . . , N} of H λ .
When the convergence is reached with respect to λ, the constrained learned set D H c = {η 1 c , . . . , η N c } is generated. This set is made up of N independent realizations {η c , = 1, . . . , N} of H c whose probability measure is p H c (η) dη (subscript or superscript "c" is introduced to designate the quantities related to the constrained learned set (posterior model)). For λ fixed in C ad,λ , D H λ is generated using a MCMC generator based on 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 deleting the transient part to rapidly reach the stationary response associated with the invariant measure for which measure p H c (η) dη is a marginal measure of this invariant measure. The ISDE is solved by using the Störmer-Verlet algorithm, which yields an efficient and accurate MCMC algorithm. Using an ergodicity property, obtaining a portion of the realization of the stationary process of the ISDE, which corresponds to its invariant measure, is computed by integrating the ISDE over a long period of time. Snapshots are taken on this portion. The successive snapshots are taken sufficiently far apart to guarantee a quasi statistical independence. This approach is efficient but is numerically "recursive" and therefore cannot be parallelized. The proposed MCMC algorithm, which is detailed in Section 4.1, makes it possible to reach the stationary solution very quickly thanks to the introduction of the dissipation term. For each realization of the initial condition and of the Wiener process, we very quickly obtain a realization associated with the invariant measure. By generating independent realizations of the initial condition and of the Wiener process, independent realizations are thus calculated and these calculations can therefore be parallelized for strongly decreasing the elapsed time on a multicore computer. The independent realizations of the initial condition are generated by scanning all the points of the training set. The parallel loops are explicitly defined in Algorithm 1, which is detailed in Section 4.4).

Formal construction of the optimal solution using the sequence {p H λ } λ

Let us assumed that the optimization problem defined by Eq. ( 2.

3) has almost one solution p H c and that p = p H c is a regular point of the continuously differentiable functional p → R ν h c (η) p(η) dη -b c . For λ 0 ∈ R + and λ ∈ C ad,λ , we define the Lagrangian,

Lag(p, λ 0 , λ) = R ν p(η) log p(η) p H (η) dη + (λ 0 -1) ( R ν p(η) dη -1) + λ , R ν h c (η) p(η) dη -b c . For all η in R ν , the pdf of H is written as p H (η) = c ν ζ(η) in which the positive-valued function ζ is integrable on R ν , is such that supp ζ = R ν , and where the positive constant c ν is such that c -1 ν = R ν ζ(η) dη. We define the sequence {p H λ } λ of pdf η → p H λ (η ; λ) on R ν , indexed by λ, such that p H λ (. ; λ) is an extremum of functional p → Lag(p, λ 0 , λ).
Using the calculus of variations yields

p H λ (η ; λ) = c 0 (λ) ζ(η) exp{-λ , h c (η) } , ∀ η ∈ R ν , (2.5)
in which c 0 (λ) is the constant of normalization that depends on λ (note that λ 0 is eliminated and we have c 0 (λ) = c ν exp{-λ 0 }). The existence of a unique solution requires that the constraints be algebraically independent in the following sense: given any bounded positive measure P(dη) on R ν with support R ν , there exists a bounded set B in R ν with P(B) > 0 such that

∀ v ∈ R n c , v = 0 , B h c (η) , v 2 P(dη) > 0 . (2.6)
Under the condition defined by Eq. (2.6), there exists (see [START_REF] Luenberger | Optimization by Vector Space Methods[END_REF])

λ sol in C ad,λ such that the functional (p, λ 0 , λ) → Lag(p, λ 0 , λ) is stationary at point p = p H c for λ = λ sol and λ 0 = -log(c 0 (λ sol )/c ν ). Consequently, p H c = p H λ sol (. ; λ sol ) and Eq. (2.5) yield p H c (η) = c 0 (λ sol ) ζ(η) exp{-λ sol , h c (η) } , ∀ η ∈ R ν . (2.7)
Taking into account the introduced hypotheses, p H c is the unique solution of the optimization problem defined by Eq. (2.3), in which λ sol is the unique solution of a convex optimization problem that will be defined by Proposition 1 in Section 3) and which will be, under required hypotheses, the solution of the following nonlinear algebraic equation in λ, R ν h c (η) p H λ (η ; λ) dη = b c .

Implicit definition of function h c , resulting difficulties, and necessity to construct a statistical surrogate model

Using the methodology presented in Section 2.4, the drift vector of the ISDE will involve the matrix [∇ η h c (η)] ∈ M ν,n c (the transpose of the Jacobian matrix of h c ). As we have explained in Section 1, an additional strong assumption used in this paper is that function η → h c (η) from R ν into R n c is not explicitly defined by an algebraic expression. It is assumed that we can only compute a = h c (η ) ∈ R n c for any point η given in R ν (for instance and as previously underlined, a component of h c can be related to the square of a norm of the random normalized residue of the stochastic PDE). The MCMC generator requires the evaluation of [∇ η h c (η)] for a large number of values of η. Consequently, a statistical surrogate model ĥN of h c is constructed and allows for deducing an algebraic representation η

→ [∇ η ĥN η)] of function η → [∇ η h c (η)] from R ν into M ν,n c
. Such statistical surrogate model is an approximation whose convergence with respect to N will be given by Proposition 2 in Section 3.

Mathematical analysis of the proposed methodology

In this section, there are some repetitions with respect to Section 2, but we have preferred to do them so that Section 3 be mathematically coherent and self contained.

, . . . , η N d d given in R ν such that η = 1 N d N d j=1 η j d = 0 ν , [ C H ] = 1 N d -1 N d j=1 (η j d -η) ⊗ (η j d -η) = [I ν ] . (3.1)
Let H = (H 1 , . . . , H ν ) be the R ν -valued random variable defined on the probability space (Θ, T , P) whose probability measure P H (dη) = p H (η) dη is defined by the probability density function η → p H (η) : R ν → R + with respect to the Lebesgue measure dη on R ν , ν+4) , and where η → ζ(η) : R ν → R + is written as

p H (η) = c ν ζ(η) , ∀η ∈ R ν , c ν = ( √ 2π ŝ) -ν , (3.2) in which ŝ = s s 2 + (N d -1)/N d -1/2 with s = (4/(N d (2 + ν))) 1/(
ζ(η) = 1 N d N d j=1 exp - 1 2 ŝ2 ŝ s η j d -η 2 . (3.3)
We define the potential function η → φ(η) : R ν → R, related to p H , such that

ζ(η) = exp{-φ(η)} . (3.4)
Remark 1. Definition 1 of p H corresponds to a Gaussian kernel-density estimation (KDE) using the training set D 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 is the Silverman bandwidth. With such a modification, the normalization of H is preserved for any value of N d ,

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

E{H ⊗ H} = R ν η ⊗ η p H (η) dη = ŝ2 [I ν ] + ŝ2 s 2 (N d -1) N d [ C H ] = [I ν ] . (3.6)
Theorem 3.1 in [START_REF] Soize | Probabilistic learning on manifolds[END_REF] proves that, for all η fixed in R ν , Eq. (3.2) with Eq. (3.3) is a consistent estimation of the sequence {p H } N d for N d → +∞.

Hypothesis 1 (Concerning function h c ). Let n c be the integer such that 1 ≤ n c ≤ ν. It is assumed that η → h c (η) verifies the property defined by Eq. (2.6), is continuously differentiable from R ν into R n c ,

h c ∈ C 1 (R ν , R n c ) , (3.7)
and there exist constants α > 0, β > 0, c α > 0, and c β > 0, independent of η, such that for η → +∞,

h c (η) ≤ c α η α , [∇ η h c (η)] F ≤ c β η β , (3.8)
in which

[∇ η h c (η)] ∈ M ν,n c with [∇ η h c (η)] αk = ∂h c k (η)/∂η α ,
C ad,λ = λ ∈ R n c | 0 < E{ exp{-λ , h c (H) } < +∞ , (3.9)
in which the pdf of the R ν -valued random variable H is defined by Eq. (3.2). It is also assumed that h c is such that C ad,λ is not reduced to the empty set, C ad,λ = ∅ .

V λ (η) = φ(η) + λ , h c (η) , (3.11) in which φ(η) = -log ζ(η) (see Eq. (3.4)). One then has 0 < R ν exp{-V λ (η)} dη < +∞ . (3.12)
(c) The pdf η → p H λ (η ; λ) with respect to dη, defined by Eq. (2.5), which can be written as

p H λ (η ; λ) = c 0 (λ) exp{-V λ (η)} , ∀η ∈ R ν , (3.13)
is such that the constant c 0 (λ) of normalization verifies

0 < c 0 (λ) < +∞ , ∀λ ∈ C ad,λ . (3.14)
(d) ∀λ ∈ C ad,λ , we have V λ (η) → +∞ if η → +∞, and we have,

R ν h c (η) 2 exp{-V λ (η)} dη < +∞ , R ν [∇ η h c (η)] F exp{-V λ (η)} dη < +∞ . (3.15) PROOF. (Lemma 1). (a) C ad,λ is a convex subset if ∀µ ∈ [0, 1], ∀λ and λ in C ad,λ , µλ + (1 -µ)λ ∈ C ad,λ . Let A = exp{-λ , h c (H) } and B = exp{-λ , h c (H) } be R + -
R ν exp{-V λ (η)} dη = 1 c ν R ν exp{-λ , h c (η) } p H (η) dη = 1 c ν E{ exp{-λ , h c (η) } },
which is positive and finite due to Eq. (3.9) and to 0 < c ν < +∞. We have thus proven Eq. (3.12). (c) Using Eqs. (3.12) and (3.13), and since we must have (3.17)

R ν p H λ (η) dη = 1, we deduce Eq. (3.14). (d) Since h c is continuous on R ν , (see Eq. (3.7)), ∀λ ∈ C ad,λ , η → exp{-V λ (η)} is continuous on R ν and then is locally integrable on R ν . Eq. (3.12) implies the integrability at infinity of η → exp{-V λ (η)}. Since η → V λ (η) is continuous on R ν , it can be deduced that V λ (η) → +∞ if η → +∞. Using Eq. (3.8), for η → +∞, one has h c (η) 2 exp{-V λ (η)} ≤ c 2 α η 2α exp{-V λ (η)} and [∇ η h c (η)] F exp{-V λ (η)} ≤ c β η β exp{-V λ (η)},
(b) Let λ → Γ(λ) : C ad,λ → R be defined by Γ(λ) = λ , b c -log c 0 (λ) , (3.18)
in which b c is given in R n c . For all λ in C ad,λ , we have

∇ λ Γ(λ) = b c -E{h c (H λ )} ∈ R n c , (3.19) [Γ (λ)] = [cov{h c (H λ )}] ∈ M + n c , (3.20)
where the positive-definite covariance matrix

[Γ (λ)] of h c (H λ ) is such that [Γ (λ)] kk = ∂ 2 Γ(λ)/∂λ k ∂λ k . (c) Γ is a strictly convex function on C ad,λ .
There is a unique solution λ sol in C ad,λ of the convex optimization problem,

λ sol = arg min λ∈C ad,λ Γ(λ) . (3.21) If the following equation in λ, ∇ λ Γ(λ) = 0 n c , (3.22
E{ h c (H λ ) 2 } = R ν h c (η) 2 c 0 (λ) exp{-V λ (η)} dη < +∞ .
(b) The definition of function Γ given by Eq. (3.18) is similar to the one introduced in the discrete case for finding the probability measure of maximal entropy [START_REF] Agmon | An algorithm for finding the distribution of maximal entropy[END_REF][START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF]. Let us prove Eqs. (3.19) and (3.20). Eq. (3.11) yields ∇ λ V λ (η) = h c (η) and from Eq. (3.13), it can be deduced that

∇ λ p H λ (η ; λ) = c 0 (λ) -1 ∇ λ c 0 (λ) -h c (η) p H λ (η ; λ) . (3.24)
Integrating Eq. (3.13) on R ν and taking the logarithm yields log c 0 (λ) =log R ν exp{-V λ (η)} dη and consequently,

c 0 (λ) -1 ∇ λ c 0 (λ) = R ν h c (η) p H λ (η ; λ) dη = E{h c (H λ )} . (3.25) Eq. (3.18) yields ∇ λ Γ(λ) = b c -c 0 (λ) -1 ∇ λ c 0 (λ),
[Γ (λ)] = - R ν h c (η) ⊗ ∇ λ p H λ (η ; λ) dη . (3.26) Substituting Eq. (3.25) into Eq. (3.24) yields ∇ λ p H λ (η ; λ) = (E{h c (H λ )} -h c (η) ) p H λ (η ; λ), which with Eq. (3.26), gives [Γ (λ)] = E{h c (H λ ) ⊗ h c (H λ )} -(E{h c (H λ )}) ⊗ (E{h c (H λ )}
) that is the covariance matrix of the R n c -valued random variable h c (H λ ). Again Eq. (3.17) proves the existence of matrix [Γ (λ)] as a covariance matrix, which is semi-positive definite. We have to prove that this matrix is positive definite, which will be true if the matrix

[M λ ] = E{h c (H λ ) ⊗ h c (H λ )} belongs to M + n c , that is to say if [M λ ] v , v 2 > 0 for all v in R n c with v = 0.
Since p H λ (η ; λ) dη is a probability measure, this will be true if η → h c (η) , v 2 is not zero on a set B such that B p H λ (η ; λ) dη > 0 that is the case due to the hypothesis defined by Eq. (2.6). (c) From Lemma 1-(a), C ad,λ is a convex set and from Eq. (3.20), [Γ (λ)] is a positive-definite matrix for all λ in C ad,λ . It can then be deduced that λ → Γ(λ) is strictly convex on C ad,λ and therefore, Eq. (3.21) holds and λ sol is unique. Note that the existence of a solution of Eq. (3.22), which would then be a global minimum, could be not in C ad,λ . However, if the equation ∇ λ Γ(λ) = 0 n c admits a solution λ = λ ∈ R n c that belongs to C ad,λ , then this solution is unique and we have λ sol = λ, which is the solution of the convex optimization problem defined by Eq. (3.21), and then Eq. (3.22) holds for λ = λ sol . Finally, under the condition that λ belongs to C ad,λ , Eq. (3.19) shows that E{h c (H λ sol )} = b c . Taking into account Eq. (2.7), the solution is given by Eq. (3.23) and is unique due to the uniqueness of solution λ sol of ∇ λ Γ(λ) = 0 n c . Definition 3 (Surrogate model ĥN of h c ). Let λ be fixed in C ad,λ and let D H λ = {η 1 λ , . . . , η N λ } be the constrained learned set whose points are N N d independent realizations of the R ν -valued random variable H λ for which the pdf η → p H λ (η ; λ) is defined by Eq. (3.13). Let A λ = h c (H λ ) be the R n c -valued random variable defined on (Θ, T , P) whose N independent realizations a 1 λ , . . . , a N λ are such that

a λ = h c (η λ ) ∈ R n c for = 1, . . . , N. The surrogate model η → ĥN (η ; λ) : R ν → R n c of h c is defined, for all η in R ν , by ĥN (η ; λ) = N =1 a λ β N η (η λ ) N =1 β N η (η λ ) , (3.27)
in which for all η and η in R ν ,

β N η ( η) = exp{- 1 2s 2 SB η -η 2 H } , (3.28) η -η 2 H = [σ H λ ] -2 ( η -η) , η -η , (3.29)
in which [σ H λ ] is the diagonal positive-definite matrix in M + ν such that [σ H λ ] αα is the standard deviation of the realvalued random variable H λ,α , estimated using D H λ , and where s SB is the Silverman bandwidth that depends on N and written as

s SB = 4 N(2 + n c + ν) 1/(n c +ν+4)
.

(3.30)

Remark 2 (Rationale of Definition 3). Let us assume that λ is fixed in C ad,λ . For all η in R ν , we have the following identity,

h c (η) = E{h c (H λ ) | H λ = η} , (3.31) in which E{h c (H λ ) | H λ = η} is the conditional mathematical expectation of the R n c -valued second-order random variable h c (H λ ) given H λ = η. Let P A λ ,H λ (da, dη ; λ) be the probability measure on R n c × R ν of the R n c × R ν - valued random variable ( A λ , H λ ). Since A λ = h c (H λ ), the support of P A λ ,H λ is the manifold defined by the graph {(a, η) ∈ R n c × R ν ; a = h c (η)}. Let P H λ (dη ; λ) = p H λ (η ; λ) dη
h c (η) = a∈R nc a P A λ | H λ (da | η ; λ) , (3.32)
in which

P A λ | H λ (da | η ; λ) is the conditional probability measure of A λ given H λ = η in R ν
, which could also be written as δ 0 (a -h c (η)) in which δ 0 is the Dirac measure on R n c at point a = 0 n c . Note that if the joint probability measure P A λ ,H λ (da, dη ; λ) of ( A λ , H λ ) had a density p A λ ,H λ (a, η ; λ) with respect to da ⊗ dη (that is not the case), then Eq. (3.32) could be written as

h c (η) = a∈R nc a p A λ | H λ (a | η ; λ) da with p A λ | H λ (a | η ; λ) = p A λ ,H λ (a | η ; λ)/p H λ (η ; λ).
The statistical surrogate model η → ĥN (η ; λ) : R ν → R n c of h c is then defined by approximating (regularizing)

P A λ | H λ (da | η ; λ) by the conditional probability measure, pN A λ | H λ (a | η ; λ) da = ( pN H λ (η ; λ)) -1 pN A λ ,H λ (a, η ; λ) da , (3.33)
in which the pdf pN A λ ,H λ (a, η ; λ) on R n c × R ν with respect to da ⊗ dη is defined by the following Gaussian kernel density representation, based on the N independent realizations {(a λ , η λ ), = 1, . . . , N} of ( A λ , H λ ), and where pN

H λ (η ; λ) = R nc pN A λ ,H λ (a, η ; λ) da. Therefore, we have pN A λ ,H λ (a, η ; λ) = 1 N N =1 ( √ 2π s SB ) n c +ν det[σ A λ ] det[σ H λ ] -1 exp - 1 2s 2 SB a λ -a 2 A + η λ -η 2 H , (3.34)
in which s SB is defined by Eq. (3.30), where [σ A λ ] is the diagonal positive-definite matrix in M + n c such that [σ A λ ] kk is the standard deviation of the real-valued random variable A λ,k , estimated using the realizations {a λ , = 1, . . . , N}, and where for all ã and a in R n c , ã -

a 2 A = [σ A λ ] -2 ( ã -a)
, ã -a . From Eq. (3.32) and using the approximation of P A λ | H λ (da | η ; λ) defined by Eq. (3.33), we have ĥN (η ; λ) = ( pN ĥN (η ; λ)]} N ). Let λ be fixed in C ad,λ and let η be fixed in R ν . Under the hypothesis h c ∈ C 1 (R ν , R n c ) (see Eq. (3.7)), ∀ε > 0, there exists a finite integer N ε (η, λ) depending on ε, η, and λ, such that ∀N ≥ N ε (η, λ),

H λ (η ; λ)) -1 a∈R nc a pN A λ ,H λ (a, η ; λ) da,
ĥN (η ; λ) -h c (η) ≤ ε , [∇ η ĥN (η ; λ)] -[∇ η h c (η)] F ≤ ε . ( 3
> 0, let η → κ η ( η ; s, λ) = ( √ 2π s) -ν (det[σ H λ ]) -1 exp{-1 2s 2 η -η 2
H } be the function defined on R ν with values in R + . Since ∀α ∈ {1, . . . , ν}, lim s→0 + {s [σ H λ ] αα } = 0, it can be seen that lim

s→0 + κ η ( η ; s, λ) d η = δ 0 ( η -η) , (3.36)
in the vector space of bounded measures, in which δ 0 ( η) is the Dirac measure on R ν at point η = 0 ν . Eq. (3.27) with Eqs. (3.28) and (3.29) can be rewritten as ĥN (η ; λ) = h N (η ; s SB , λ) in which for all s > 0, h N (η ; s, λ) =

(1/N) N =1 h c (η λ ) κ η (η λ ; s, λ) (1/N) N =1 κ η (η λ ; s, λ) -1
. Since h c , p H λ ( . ; λ), and κ η ( . ; s, λ) are continuous on R ν , and since for any value of N, η 1 λ , . . . ,

η N λ are N independent realizations of R ν -valued random variable H λ , for s > 0 fixed, lim N→+∞ h N (η ; s, λ) = E{h c (H λ ) κ η (H λ ; s, λ)} E{κ η (H λ ; s, λ)} -1
, which can be rewritten as

lim N→+∞ h N (η ; s, λ) = R ν h c ( η) p H λ ( η ; λ) κ η ( η ; s, λ) R ν p H λ ( η ; λ) κ η ( η ; s, λ) -1 . Since h c belongs to C 1 (R ν , R n c )
and p H λ ( . ; λ) to C 0 (R ν , R + ), using Eq. (3.36) yields, for η fixed in R ν , lim s→0

+ lim N→+∞ h N (η ; s, λ) = h c (η) ∈ R ν and lim s→0 + lim N→+∞ [∇ η h N (η ; s, λ)] = [∇ η h c (η)] ∈ M ν,n c . Consequently, ∀λ ∈ C ad,λ , ∀η ∈ R ν , ∀ε > 0,
there exists s ε > 0 and a finite integer N ε (η, λ) depending on ε, η, and λ, such that

h N ε (η ; s ε , λ) -h c (η) ≤ ε and [∇ η h N ε (η ; s ε , λ)] -[∇ η h c (η)] ≤ ε. The Silverman bandwidth s SB = s SB (N) defined by Eq. (3.30), goes to 0 when N → +∞. Therefore, choosing N ≥ N ε (η, λ) such that s SB (N) < s ε (that is always possible) yields h N (η ; s SB (N), λ) -h c (η) ≤ ε and [∇ η h N (η ; s SB (N), λ)] -[∇ η h c (η)] ≤ ε, which proves Eq. (3.35).
Notation (Normalized R ν -valued Wiener stochastic process). Let {W wien (t) = (W wien 1 (t), . . . , W wien ν (t)), t ≥ 0} be the Wiener process, defined on (Θ, T , P), indexed by R + , with values in R ν , such that W wien 1 , . . . , W wien ν are mutually independent, W wien (0) = 0 ν a.s., W wien is a process with independent increments such that, for all 0 ≤ t < t < +∞, the increment W wien (t) -W wien (t ) is a R ν -valued second-order, Gaussian, centered random variable whose covariance matrix is (tt ) [I ν ].

Proposition 3 (MCMC generator of H λ ). Let h c be the function satisfying Hypothesis 1. Let λ be fixed in C ad,λ . Consequently, Lemma 1-(d) holds. Let {(U λ (t), V λ (t)), t ≥ 0} be the stochastic process, defined on (Θ, T , P), indexed by R + , with values in R ν × R ν , which verifies the following ISDE for t > 0, with the initial condition 2) while the initial condition v 0 is chosen as any realization of a normalized Gaussian R ν -valued random variable V G , independent of W wien , whose probability density function with respect to dv is p V G (v) = (2π) -ν/2 exp{-v 2 /2}. (b) The parameter f 0 > 0 allows the dissipation term in the dissipative Hamiltonian system to be controlled and to rapidly reach the stationary response associated with the invariant measure (the value f 0 = 4 is generally a good choice). (c) For all u = (u 1 , . . . , u ν ) in R ν , the vector L λ (u) in R ν is defined by L λ (u) = -∇ u V λ (u) that can be written as

(u 0 , v 0 ) given in R ν × R ν , dU λ (t) = V λ (t) dt , (3.37)
dV λ (t) = L λ (U λ (t)) dt - 1 2 f 0 V λ (t) dt + f 0 dW wien (t) , (3.38)
U λ (0) = u 0 , V λ (0) = v 0 a.s. ( 3
L λ (u) = 1 ζ(u) ∇ u ζ(u) -[∇ u h c (u)] λ . ( 3
(τ), V st λ (τ)), τ ≥ 0} associated with the unique invariant probability measure on R ν × R ν , p H λ ,V G (η, v ; λ) dη ⊗ dv = (p H λ (η ; λ) dη) ⊗ (p V G (v) dv) , (3.41)
in which p H λ (η ; λ) is the pdf defined by Eq. (3.13).

(e) For t s sufficiently large, we can choose H λ as U λ (t s ). The generation of the constrained learned set D H λ = {η 1 λ , . . . , η N λ }, made up of N N d independent realizations of H λ whose probability measure is p H λ (η ; λ) dη, consists in solving Eqs. (3.37) to (3.39) for t ∈ [0 , t s ] and then using the realizations of U λ (t s ) (see the numerical aspects in Section 4).

PROOF. (Proposition 3). Since h

c ∈ C 1 (R ν , R n c ) and φ(u) = -log ζ(u) with ζ(u)
given by Eq. (3.3), function u → V λ (u) defined by Eq. (3.11) belongs to C 1 (R ν , R). Therefore, u → ∇ u V λ (u) is locally bounded on R ν . Using Eqs. (3.11) and (3.12), it can be seen that, for all λ ∈ C ad,λ , inf u >R V λ (u) → +∞ if R → +∞, and inf u∈R ν V λ (u) is a finite real number. Using Eqs. (3.3), (3.4), and (3.11) yields

R ν ∇ u V λ (u) p H λ (u ; λ) du ≤ R ν 1 ζ(u) ∇ u ζ(u) p H λ (u ; λ) du + R ν [∇ u h c (u)] F λ p H λ (u ; λ) du , (3.42) because [∇ u h c (u)] λ ≤ [∇ u h c (u)] λ and [∇ u h c (u)] ≤ [∇ u h c (u)] F .
From Eqs. (3.12) and (3.3), the first term in the right-hand side member of Eq. (3.42) is finite, while from the second equation (3.15), the second term is also finite. It can then be deduced that the left-hand side member of Eq. (3.42) is finite. Consequently, Theorems 6, 7, and 9 in Pages 214 to 216 of [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF], and the expression of the invariant measure given by Theorem 4 in Page 211 of the same reference, for which the Hamiltonian is H(u, v) = v 2 /2 + V λ (u), prove that the solution of Eqs. (3.37) to (3.39) is unique and is a second-order diffusion stochastic process with almost-surely continuous trajectories, which converges for t → +∞ to a second-order stationary diffusion process with almost surely continuous trajectories {(U st λ (τ), V st λ (τ)), τ ≥ 0} associated with the invariant probability measure given by Eq. (3.41). For any τ > 0, U st λ (τ) = lim t→+∞ U λ (t + τ) in probability measure.

Proposition 4 (Convergence of the sequence of MCMC generator using the statistical surrogate model). Let λ be fixed in C ad,λ and let us use Proposition 3. For all η in R ν , let ĥN (η ; λ) be the approximation of h c (η) defined by Eq. (3.27) and let u → Lλ (u) be the twice continuously differentiable function on R ν with values in R ν such that, for all

u in R ν , LN λ (u) = 1 ζ(u) ∇ u ζ(u) -[∇ u ĥN (u ; λ)] λ , (3.43
dV N λ (t) = LN λ (U N λ (t)) dt - 1 2 f 0 V N λ (t) dt + f 0 dW wien (t) , (3.45)
U N λ (0) = u 0 , V N λ (0) = v 0 a.s. , (3.46)
and where u 0 , v 0 , f 0 , and W wien are the quantities defined in Proposition 3. Then the stochastic solution {(U N λ (t), V N λ (t)), t ≥ 0} of Eqs. (3.44) to (3.46) is unique, has almost-surely continuous trajectories, and is a second-order diffusion stochastic process, which converges to a stationary second-order diffusion stochastic process for t → +∞, associated with the unique invariant probability measure on

R ν × R ν , pN H λ ,V G (η, v ; λ) dη ⊗ dv = ( pN H λ (η ; λ) dη) ⊗ (p V G (v) dv), in which pN H λ (η ; λ) = ĉN 0 (λ) exp{-V N λ (η)} with V N λ (η) = -log ζ(ζ) + λ , ĥN (η ; λ)
. Then for all t ∈ [0 , t s ] with t s < +∞, the sequence {U N λ (t)} N of second-order R ν -valued random variables converges in mean-square to the second-order R ν -valued random variable U λ (t) of Proposition 3,

lim N→+∞ E{ U N λ (t) -U λ (t) 2 } = 0 , ∀t ∈ [0 , t s ] . (3.47)
PROOF. (Proposition 4). The classical theorem, such as Theorem 5.1 Page 118 of [START_REF] Friedman | Stochastic Differential Equations and Applications[END_REF], cannot directly be used because, for η → +∞, the growth of the drift vector must be at most linear, which is not verified taking into account Eq. (3.40) and the fact that [∇ η h c (η)] F can have any polynomial growth at infinity (see Eq. (3.8)). Consequently, an adapted proof of this Proposition 4 must be done. For λ be fixed in C ad,λ , the unique second-order stochastic process with almost-surely continuous trajectories {(U λ (t), V λ (t)), t ≥ 0} of Proposition 3 can be written as

Z λ (t) = z 0 + t 0 a λ (Z λ (τ)) dτ + t 0 [b] dW wien (τ) , (3.48)
in which z 0 = (u 0 , v 0 ), a λ (z) = (v , L λ (u) -(1/2) f 0 v) and where z = (u, v) with u and v in R ν , where z 0 , z, and a λ (z) are in R 2ν = R ν × R ν , and where

[b] = [ [0 ν ] f 0 [I ν ] ] T ∈ M 2ν,ν .
Reusing the proof of Proposition 3, it can be seen that Eqs. (3.44) to (3.46) admits a unique solution {(U N λ (t), V N λ (t)), t ≥ 0} (with the properties given in Proposition 3), which can be written as,

Z N λ (t) = z 0 + t 0 a N λ (Z N λ (τ)) dτ + t 0 [b] dW wien (τ) , (3.49) in which a N λ (z) = (v , LN λ (u) -(1/2) f 0 v).
Subtracting the two equations Eqs. (3.48) and (3.49) yields

Z N λ (t) -Z λ (t) = χ N λ (t) + t 0 a N λ (Z N λ (τ)) -a N λ (Z λ (τ)) dτ , χ N λ (t) = t 0 a N λ (Z λ (τ)) -a λ (Z λ (τ)) dτ . (3.50)
Let L 2 (Θ, R m ) be the Hilbert space of the equivalent classes of second-order R m -valued random variables defined on (Θ, T , P), equipped with the inner product A , A = E{ A , A } and the associated norm

||| A |||= (E{ A 2 }) 1/2 . Let us define G(τ) = a N λ (Z N λ (τ)) -a N λ (Z λ (τ)) and 1(τ) = 1. The Schwarz inequality t 0 1(τ) × ||| G(τ) ||| dτ ≤ ( t 0 1(τ) 2 dτ) 1/2 ( t 0 ||| G(τ) ||| 2 dτ) 1/2 yields t 0 ||| G(τ) ||| dτ 2 ≤ t t 0 ||| G(τ) ||| 2 dτ. Using the inequality 2αβ ≤ 2α 2 + β 2 /2 for all α > 0 and β > 0 and since ||| a N λ (Z N λ (τ)) -a N λ (Z λ (τ)) ||| 2 = (1 + f 2 0 /4) ||| V N λ (τ) -V λ (τ) ||| 2 +||| LN λ (U N λ (τ)) - LN λ (U λ (τ)) ||| 2 , we obtain ||| Z N λ (t) -Z λ (t) ||| 2 ≤ 3 ||| χ N λ (t) ||| 2 + 3 2 t t 0       (1+ f 2 0 4 ) ||| V N λ (τ) -V λ (τ) ||| 2 +||| LN λ (U N λ (τ)) -LN λ (U λ (τ)) ||| 2       dτ. (3.51)
Let t be fixed such that 0 < t ≤ t s < +∞. Since {U N λ (τ), τ ≥ 0} and {U λ (τ), τ ≥ 0} are dependent second-order R ν -valued stochastic process with almost-surely continuous trajectories, there exists a finite positive constant r λ (t s ) depending on λ and t s , 0

< r λ (t s ) < +∞, such that sup 0≤τ≤t s ||| U N λ (τ) ||| < r λ (t s ) and sup 0≤τ≤t s ||| U λ (τ) ||| < r λ (t s ). Let U t s be the open ball of L 2 (Θ, R ν ) such that U t s = U ∈ L 2 (Θ, R ν ) ; ||| U ||| < r λ (t s ) . Due to the convexity of the open ball in a normed vector space, U t s is a convex open set of L 2 (Θ, R ν ) (∀U, U ∈ U t s , ∀µ ∈ [0 , 1], we have ||| (1 -µ) U + µ U ||| ≤ (1 -µ) ||| U |||+µ |||U ||| ≤ (1 -µ) r λ (t s ) + µ r λ (t s ) = r λ (t s ), which shows that (1 -µ) U + µ U ∈ U t s ). Let U → LN λ (U) be the mapping from L 2 (Θ, R ν ) into L 2 (Θ, R ν ),

in which LN

λ is defined by Eq. (3.43). Since u → ĥN (u ; λ) is twice continuously differentiable on R ν , then function u

→ LN λ (u) is continuously differentiable on R ν . It can easily be verified that ∀U ∈ U t s , E{ [∇ u LN λ (U)] 2 F } 1/2
≤ k λ (t s ), in which k λ (t s ) is a finite positive constant depending on λ and t s . Consequently, using Theorem 3.3.2 Page 45 of [START_REF] Cartan | Cours de Calcul Différentiel[END_REF] for Banach spaces, for all 0 ≤ τ ≤ t ≤ t s , we have

||| LN λ (U N λ (τ)) -LN λ (U λ (τ)) ||| ≤ k λ (t s ) ||| U N λ (τ) -U λ (τ) ||| . (3.52)
From Eqs. (3.51) and (3.52), it can be deduced using the Gronwall Lemma (see [START_REF] Gronwall | Note on the derivatives with respect to a parameter of the solutions of a system of differential equations[END_REF] or [START_REF] Dieudonné | Calcul Infinitésimal[END_REF], Page 362) that, for all 0

≤ t ≤ t s < +∞, ||| Z N λ (t) -Z λ (t) ||| 2 ≤ c λ (t s ) ||| χ N λ (t) ||| 2 , (3.53) with c λ (t s ) = (9/2) t 2 s max{(1 + f 2 0 /4) , k λ (t s ) 2 } < +∞. Eq. (3.50) yields ||| χ N λ (t) ||| ≤ t 0 ||| LN λ (U λ (τ)) -L λ (U λ (τ)) ||| dτ. Using Eqs. (3.40) and (3.43) yields ||| χ N λ (t) ||| ≤ λ t 0 E{ [∇ u ĥN (U λ (τ) ; λ)] -[∇ u h c (U λ (τ))] 2 F 1/2 dτ . (3.54)
From Proposition 2, it can be deduced that, for N → +∞, the right-hand side member of Eq. (3.54) goes to 0 and consequently, Eq. (3.53) shows that Z N λ (t) → Z λ (t) for the mean-square convergence.

Iterative algorithm for calculating λ sol . Under Proposition 1, for λ ∈ C ad,λ , since Γ(λ) cannot be evaluated in high dimension using Eq. (3.18) due to the presence of constant c 0 (λ) (the normalization constant), λ sol cannot directly be estimated using the gradient descent algorithm applied to the convex optimization problem defined by Eq. (3.21). We will then assumed that λ sol can be calculated as the unique solution in λ of equation ∇ λ Γ(λ) = 0 n c (see Proposition 1-(c)) and in particular Eq. (3.22)), that is to say (see Eq. (3.19)), solving the following equation in λ on R n c ,

E{h c (H λ )} -b c = 0 n c . (3.55)
This equation is solved using the Newton iterative method [START_REF] Kelley | Solving Nonlinear Equations With Newton's Method[END_REF] applied to function λ → ∇ λ Γ(λ) as proposed in [START_REF] Batou | Calculation of Lagrange multipliers in the construction of maximum entropy distributions in high stochastic dimension[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF], that is to say,

λ i+1 = λ i -[Γ (λ i )] -1 ∇ λ Γ(λ i ) , i = 0, 1, . . . , i max , (3.56)
with λ 0 = 0 n c , in which ∇ λ Γ(λ) and [Γ (λ)] are defined by Eqs. (3.19) and (3.20), and where i max is a given integer sufficiently large. An estimation of λ sol is chosen as

λ sol = λ i sol , i sol = arg min i=1,...,i max err(i) , (3.57)
in which the error function i → err(i) : {1, . . . , i max } → R + is defined as follows for taking into account the possible types of algebraic quantities (scalar, vectors, tensors) that are used for defining function h c . Therefore, let M be an integer such that 1 ≤ M ≤ n c and for which h c and b c are written as

h c (η) = (h c,k 1 (η), . . . , h c,k M (η)) and b c = (b c,k 1 , . . . , b c,k M ) with M m=1 k m = n c
.The error function is then defined by

err(i) =        M m=1 w m err m (i) err m (1) 2        1/2 , ( 3
(i) = 1 b c,k m b c,k m -E{h c,k m (H λ i )} . (3.59)
Proposition 5 (Rate of convergence of the sequence {p H λ } λ ). Let λ i+1 and λ i be given values of λ in C ad,λ and let p H λ ( . ; λ) be the pdf of H λ defined by Eq. (3.13). For λ i+1 -λ i sufficiently small, we have

p H λ i+1 ( . ; λ i+1 ) -p H λ i ( . ; λ i ) L 1 (R ν ,R) ≤ λ i+1 -λ i tr [Γ (λ i )] 1/2 + o ( λ i+1 -λ i ) , (3.60)
in which [Γ (λ i )] ∈ M + n c is defined by Eq. (3.20) for λ = λ i .

PROOF. (Proposition 5

). In this proof, for simplifying the writing, λ i is simply written as λ. Proposition 1 shows that Γ is twice differentiable in C ad,λ . For all η fixed in R ν , the Taylor expansion of λ i+1 → p H λ i+1 (η ; λ i+1 ) around λ, truncated at the first order, is written as,

p H λ i+1 (η ; λ i+1 ) = p H λ (η ; λ) + ∇ λ p H λ (η ; λ) , λ i+1 -λ + . . . (3.61)
Eq. (3.18) is written as c 0 (λ) = exp{ λ , b c -Γ(λ) } and is substituted in Eq. (3.13) (or equivalently, in Eq. (2.5)) of

p H λ (η ; λ) yielding p H λ (η ; λ) = ζ(η) exp{-Γ(λ) -λ , h c (η) -b c }.
The gradient with respect to λ can be written as

∇ λ p H λ (η ; λ) = -∇ λ Γ(λ) + h c (η) -b c p H λ (η ; λ) . (3.62)
Let us introduce the score variable

v(η ; λ) = ∇ λ log p H λ (η ; λ) = p H λ (η ; λ) -1 ∇ λ p H λ (η ; λ) , (3.63)
which can be rewritten, using Eq. (3.62), as v(η ; λ) = -∇ λ Γ(λ) + h c (η) -b c yielding with the use of Eq. (3.19),

v(η ; λ) = -h c (η) -E{h c (H λ )} . (3.64) Eq. (3.64) shows that E{ v(H λ ; λ) 2 } = E h c (H λ ) -E{h c (H λ )} 2 = tr [cov{h c (H λ )}
] and using Eq. (3.20) yields,

E{ v(H λ ; λ) 2 } = R ν v(η; λ) 2 p H λ (η ; λ) dη = tr [Γ (λ)] . (3.65)
Eq. (3.63) can be written as

∇ λ p H λ (η ; λ) = p H λ (η ; λ) v(η ; λ). Consequently, R ν | ∇ λ p H λ (η ; λ) , λ i+1 -λ | dη = R ν | v(η ; λ) , λ i+1 -λ | p H λ (η ; λ) dη ≤ λ i+1 -λ R ν v(η ; λ) p H λ (η ; λ) dη ≤ λ i+1 -λ R ν p H λ (η ; λ) dη 1/2 R ν v(η ; λ) 2 p H λ (η ; λ) dη 1/2 . (3.66)
Since R ν p H λ (η ; λ) dη = 1 and using Eq. (3.65), we obtain

R ν | ∇ λ p H λ (η ; λ) , λ i+1 -λ | dη ≤ tr [Γ (λ)] 1/2 λ i+1 -λ . (3.67)
From Eq. (3.61), it can be deduced that

R ν | p H λ i+1 (η ; λ i+1 ) -p H λ (η ; λ) | dη ≤ R ν | ∇ λ p H λ (η ; λ) , λ i+1 -λ | dη + o ( λ i+1 -λ i ) that
D H λ = {η 1 λ , . . . , η N λ } with N N d . It is assumed that N is written as N = N d × n MC with n MC
1. Therefore, a discretization scheme [START_REF] Kloeden | Numerical Solution of Stochastic Differentials Equations[END_REF][START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] must be used. The case of Hamiltonian dynamical systems has been analyzed in [START_REF] Talay | Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF] by using an implicit Euler scheme. The Störmer-Verlet scheme (see [START_REF] Hairer | Geometric numerical integration illustrated by the Störmer-Verlet method[END_REF] for the deterministic case and [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF] for the stochastic case) is a very efficient scheme that allows for having a long-time energy conservation for non-dissipative Hamiltonian dynamical systems. In [START_REF] Soize | Time-domain formulation in computational dynamics for linear viscoelastic media with model uncertainties and stochastic excitation[END_REF], we have proposed to use an extension of the Störmer-Verlet scheme for stochastic dissipative Hamiltonian systems, that we have also used in [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF][START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF].

Störmer-Verlet scheme and computation of the constrained learned set D H λ i

Let i be the index of the sequence {λ i , i = 0, 1, . . . , i max } of Lagrange multipliers computed using Eq. (3.56) with λ 0 = 0 n c . Let t m = m ∆t for m = 0, 1, . . . , M s (with M s > 1 an integer) be the time sampling in which t s = M s ∆t (and thus t M s = t s ). Let ∆W wien m+1 = W wien (t m+1 ) -W wien (t m ) be the Gaussian, second-order, centered, R ν -valued random variable such that E{∆W wien m+1 ⊗ ∆W wien m+1 } = ∆t [I ν ]. Let {θ , = 1, . . . , N} be N independent realizations in Θ. For m = 0, 1, . . . , M s -1, let ∆W m+1 = ∆W wien m+1 (θ ) be the realization θ of ∆W wien m+1 . Following the choice of (u 0 , v 0 ) defined in Proposition 3-(a), let u 1 0 , . . . , u N 0 in R ν be such that for k = 1, . . . , n MC and for j = 1, . . . , N d , we take u 0 = η j d with = j + (k -1) × N d . Let v 1 0 , . . . , v N 0 in R ν be N independent realizations of the R ν -valued random variable V G also defined in . Note that the realizations ∆W m+1 , u 0 , and v 0 , for = 1, . . . , N are independent of {λ i } i . For i ∈ {0, 1, . . . , i max } and for ∈ {1, . . . , N}, we introduce the realizations

U i, m = U N λ i (t m ; θ ) and V i, m = V N λ i (t m ; θ ). For m ∈ {0, 1, . . . , M s -1},
U i, m+1/2 = U i, m + ∆t 2 V i, m , (4.1)
V i, m+1 = 1 -γ 1 + γ V i, m + ∆t 1 + γ LN λ i-1 (U i, m+1/2 ) + f 0 1 + γ ∆W m+1 , (4.2)
U i, m+1 = U i, m+1/2 + ∆t 2 V i, m+1 , (4.3)
with the initial condition

U i, 0 = u 0 , V i, 0 = v 0 , (4.4)
in which γ = f 0 ∆t/4 and where, using Eq. (3.43),

LN λ i-1 (u) = 1 ζ(u) ∇ u ζ(u) -[∇ u ĥN (u ; λ i-1 )] λ i-1 (4.5)
The matrix [∇ u ĥN (u ; λ i-1 )] is given by Eq. (4.8), which depends on

D H λ i-1 = {η 1 λ i-1 , . . . , η N λ i-1 } . (4.6)
It should be noted that, in Eq. (4.2), LN λ i-1 has been used instead of LN λ i because LN λ i depends on D H λ i that is unknown, recalling that the aim of the recurrence defined by Eqs. (4.1) to (4.4) is precisely to calculate D H λ i that is written as

D H λ i = {η 1 λ i , . . . , η N λ i } , η λ i = U N λ i (t s ; θ ) = U i, M s . (4.7)

Explicit expression of the gradient of the statistical surrogate model of h c

Using Definition 3, for fixed value of λ ∈ C ad,λ , the gradient [∇ η ĥN (η ; λ)] ∈ M ν,n c at point η ∈ R ν of the statistical surrogate model ĥN of h c can be written as

[∇ η ĥN (η ; λ)] = N =1 γ λ ⊗ a λ , γ λ = ∇ η        β N η (η λ ) N =1 β N η (η λ )        , (4.8)
in which a λ = h c (η λ ) ∈ R n c and where γ λ ∈ R ν is explicitly calculated using Eqs. (3.28) and (3.29).

Summary of the complete algorithm

The algorithm for calculating λ sol and

D H c = {η 1 c , . . . , η N c } with η c = η λ sol for = 1, . . . , N is summarized in Algorithm 1.
λ i+1 = λ i -α relax [Γ (λ i )] -1 ∇ λ Γ(λ i
λ i ← λ i+1 14: D H λ i-1 ← D H λ i 15: end for 16: λ sol = λ i sol , i sol = arg min i err(i) from Eq. (3.57) 17: D H c ← D H λ sol

Application to stochastic homogenization without scale separation

In this section, we consider the stochastic boundary value problem associated with the stochastic homogenization of a random elastic medium without scale separation, which has been presented in Section 1. The physical space R 3 is referred to a Cartesian reference system whose the generic point is ξ = (ξ 1 , ξ 2 , ξ 3 ). We consider the stochastic homogenization of a heterogeneous linear elastic microstructure occupying the 3D bounded open domain 3 (square thick plate) with boundary ∂Ω. The homogenization method on Ω in the one proposed in [START_REF] Bornert | Homogenization in Mechanics of Materials[END_REF] that we have already used in [START_REF] Soize | Tensor-valued random fields for meso-scale stochastic model of anisotropic elastic microstructure and probabilistic analysis of representative volume element size[END_REF][START_REF] Soize | Computational stochastic homogenization of heterogeneous media from an elasticity random field having an uncertain spectral measure[END_REF]. In this section, we use the convention for summation over repeated Latin indices taking values in {1, 2, 3}.

Ω =] 0 , 1 [×] 0 , 1 [×] 0 , 0.1 [⊂ R

Stochastic elliptic boundary value problem

For all m and r in {1, 2, 3} the unknown field is the R 3 -valued random field

{Y(ξ) = (Y 1 (ξ), Y 2 (ξ), Y 3 (ξ))
, ξ ∈ Ω} defined on (Θ, T , P), indexed by Ω, such that for i = 1, 2, 3, and almost surely,

- ∂ ∂ξ j C i jpq (ξ) ε pq (Y mr (ξ)) = 0 3 , ∀ξ ∈ Ω , (5.1)
Y mr (ξ) = y mr 0 , ∀ξ ∈ ∂Ω , (5.2)
in which the strain tensor is ε pq (y) = (∂y p /∂ξ q +∂y q /∂ξ p )/2 for all y = (y 1 , y 2 , y 3 ). For all ξ ∈ ∂Ω, y mr 0 = (y mr 0,1 , y mr 0,2 , y mr 0,3 ) is defined by

y mr 0, j = (δ jm ξ r + δ jr ξ m )/2 , j ∈ {1, 2, 3} , (5.3)
in which δ jm is the Kronecker symbol. At mesoscale, the linear elastic heterogeneous medium is described by the random apparent elasticity field {C(ξ), ξ ∈ R 3 }, which is a non-Gaussian fourth-order tensor-valued random field C = {C i jpq } i jpq with i, j, p, and q in {1, 2, 3}, defined on (Θ, T , P). The stochastic homogenization consists, for i, j, m, and r in {1, 2, 3}, in analyzing at macroscale the component C eff i jmr of the random effective/apparent elasticity tensor {C eff i jmr } i jmr , which is defined by

C eff i jmr = 1 |Ω| Ω C i jpq (ξ) ε pq (Y mr (ξ)) dξ , (5.4)
in which Y mr is the R 3 -valued random field that satisfies Eqs. (5.1) to (5.3) and where |Ω|= Ω dξ. The random effective/apparent elasticity tensor C eff is symmetric and positive definite almost surely. If there was a scale separation, then the statistical fluctuations of this tensor would be negligible.

Prior probability model of random field C

The prior probability model of C used for generating the training set is the one presented in [START_REF] Soize | Computational stochastic homogenization of heterogeneous media from an elasticity random field having an uncertain spectral measure[END_REF][START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF]. This is a second-order, non-Gaussian, positive-definite fourth-order tensor-valued homogeneous random field, indexed by R 3 , defined on (Θ, T , P), with a spectral random measure. This random field is parameterized as

C(ξ) = c(G(ξ), z) , ξ ∈ Ω , (5.5)
in which {G(ξ), ξ ∈ R 3 } is a non-Gaussian second-order, homogeneous, R 21 -valued random field indexed by R 3 , defined on (Θ, T , P), with random spectral measure, where z = (z 1 , z 2 , z 3 ) is the nominal value of a R 3 -valued control parameter, and where c is a given mapping from R 21 × R 3 into the fourth-order tensor on R 3 . (ii) Anisotropic statistical fluctuations at mesoscale. The statistical fluctuations of the random medium are assumed to be anisotropic, which means that, for all ξ fixed in R 3 , the random apparent elasticity tensor C(ξ) is a full anisotropic tensor. The hyperparameters that control the anisotropic statistical fluctuations of the random apparent elasticity field (see [START_REF] Soize | Computational stochastic homogenization of heterogeneous media from an elasticity random field having an uncertain spectral measure[END_REF]) are:

(1) the dispersion coefficient δ C = 0.3 that controls the level of statistical fluctuations of the random medium.

(2) the spatial correlation lengths L c1 , L c2 , and L c3 (for directions ξ 1 , ξ 2 , and ξ 3 ) of the random field {G(ξ), ξ ∈ R 3 } and the dispersion coefficient δ s = 0.1 that controls the level of uncertainties of its spectral measure (these spectral-measure uncertainties will not be controlled, which means that δ s will not be a control parameter and its value is fixed). As explained at the end of Section 1, three cases, SC1, SC2, and SC3, of the correlation lengths are considered for analyzing the level of scale separation and are defined in Table 1. Taking into account the definition of domain Ω and the values of the spatial correlation lengths, there will not have a scale separation and consequently, the effective/apparent elasticity tensor will exhibit statistical fluctuations. (iv) Random control parameter for the probabilistic learning inference. For estimating the posterior model using the probabilistic learning inference methodology, in addition to the prior probability model of random field G, we introduce a R n w -valued random control parameter W = (W 1 , . . . , W n w ) defined on (Θ, T , P) and independent of G such that, n w = 3 and

W 1 = log C bulk , W 2 = log C shear , W 3 = log δ C , (5.6)
in which (a) C bulk and C shear are Gamma independent random variables (see [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF]) whose mean values are the nominal values C bulk and C shear previously defined and for which the coefficient of variation of C bulk is chosen as δ bulk = 0.5 yielding δ shear = 0.25 (note that, with the model proposed in [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF], δ shear is deduced from δ bulk and cannot arbitrarily be chosen); (b) δ C is chosen as a uniform random variable on [0.1 , 0.5] whose mean values is δ C .

It should be noted that the statistical fluctuations of C bulk , C shear , and δ C are chosen sufficiently large in order that the range of fluctuations of the random effective/apparent elasticity tensor covers the experimental target (see Section 5.5) in order to be able to improve the prior probabilistic model with the posterior probabilistic model by solving the inverse statistical problem with the proposed probabilistic learning inference approach.

Solution of the stochastic BVP

Under the hypotheses introduced for constructing random field C, Proposition 5.1 of [START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF] proves that for 1 ≤ m ≤ r ≤ 3, the strong stochastic solution {Y mr (ξ), ξ ∈ Ω} of the weak formulation of the stochastic elliptic BVP defined by Eqs. (5.1) to (5.3) exists, is unique, and is a second-order random field, E{ Y mr (ξ) 2 } < +∞ , ∀ξ ∈ Ω .

(5.7)

Due to Corollary 5.1 of [START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF] and its proof, the random effective/apparent elasticity tensor C eff is a second-order random variable,

X = x + [Φ] [κ] 1/2 H , (5.20)
in which [κ] is the diagonal matrix such that κ α = [κ] αα = S 2 α /(N d -1), and where H is the R ν -valued random variable whose N d independent realizations are computed by

η j d = [κ] -1/2 [Φ] T (x j d -x) , j = 1, . . . , N d . (5.21)
Note that {κ α } α are the eigenvalues of the covariance matrix of X estimated with {x 1 d , . . . , x N d d }. Random vector H is then normalized. The empirical estimation of its mean value and its covariance matrix are given by Eq. (3.1) (centered and identity matrix). From Eq. (5.20), it can be deduced that, for 1 ≤ m ≤ r ≤ 3, we have

Y mr = y mr + [Φ mr y ] [κ] 1/2 H , y mr ∈ R n y , [Φ mr y ] ∈ M n y ,ν , (5.22)
G = g + [Φ g ] [κ] 1/2 H , g ∈ R n g , [Φ g ] ∈ M n g ,ν , (5.23)
W = w + [Φ w ] [κ] 1/2 H , w ∈ R n w , [Φ w ] ∈ M n w ,ν .
(5.24)

The training set D d introduced in Definition 1 is written as

D d = {η 1 d , . . . , η N d d } , η j d ∈ R ν .
(5.25)

Figure 1 displays the distribution of the eigenvalues κ α for the 3 cases, SC1, SC2, and SC3 of scale separation, computed with the training set and the prior probability model.

Definition of the random normalized residue induced by the use of the constrained learned set

The construction of the random normalized residue is based on a similar approach of the one that has been used in [START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF]. Let us consider the constrained learned set D H λ i = {η 1 λ i , . . . , η N λ i } generated with Algorithm 1 for iteration i. Using Eqs. (5.22) to (5.24), for ∈ {1, . . . , N}, the corresponding realizations y mr, λ i ∈ R n y for 1 ≤ m ≤ r ≤ 3, g λ i ∈ R n g , and w λ i ∈ R n w , are computed by

y mr, λ i = y mr + [Φ mr y ] [κ] 1/2 η λ i , (5.26)
g λ i = g + [Φ g ] [κ] 1/2 η λ i , (5.27)
w λ i = w + [Φ w ] [κ] 1/2 η λ i .
(5.28)

For 1 ≤ m ≤ r ≤ 3 and for ∈ {1, . . . , N} the realization R mr,

λ i
of the R n y -valued random residue are computed using Eq. (5.9) and Eqs. (5.26) to (5.28), R mr,

λ i = N mr (y mr, λ i , g λ i , w λ i ).
We define the realization ρ λ i of the random residue ρλ i by

ρ λ i = 1 6 n y         1≤m≤r≤3 R mr, λ i 2         1/2 .
(5.29)

Finally, we define the realization ρ λ i of the random normalized residue ρ λ i by

ρ λ i = ρ λ i ρ 0 , (5.30)
in which ρ 0 is the estimation of E{ ρ λ 1 } using the constrained learned set D H λ 1 of the first iteration i = 1. Since for i = 1, λ i-1 = λ 0 = 0 n c , ρ λ 1 is the random normalized residue of the constrained learned set without taking into account the constraints and we have, for i = 1, E{ ρ λ 1 } = E{ ρλ 1 }/ρ 0 = 1.

Defining function h c related to the constraints and defining the targets represented by b c

We define the function η → h c (η) : R ν → R n c related to the constraints (see Eq. (2.2)) and we define the target represented by vector b c given in R n c . Three constraints are introduced, the second-order moment of the random normalized residue (see Section 5.7) and two statistical moments: the normalized mean value of the random effective/apparent elasticity matrix [C eff ] and its coefficient of dispersion (see Section 5.5). Below we consider iteration i, and then λ 1 , . . . , λ i-1 , λ i are known.

(i) The random normalized residue ρ λ i whose realizations are defined by Eq. (5.30), is an implicit function of H λ i , that we can rewrite as ρ λ i (H λ i ). The second-order moment of the normalized random residue is then written as (iii) The last constraint is given by Eq. (5.17) that is rewritten, using Eq. (3.55), as E{ h c δ (H λ i )} = b c δ , in which b c δ = (µ eff δ exp /µ exp ) 2 and where η → h c δ (η) : R ν → R is a positive-valued implicit function.

E{(ρ λ i (H λ i )) 2 } yielding h c ρ (H λ i ) = (ρ λ i (H λ i )) 2 .
(iv) Finally, for given λ, and in particular for λ = λ sol yielding H c = H λ , the constraint is defined by Eq. (3.55) with

η → h c (η) = (h c ρ (η), h c C (η), h c δ (η)) : R ν → R n c = R × R 21 × R , (5.31) b c = (b c ρ , b c C , b c δ ) ∈ R n c = R × R 21 × R , (5.32)
in which n c = 23. Using the mathematical developments presented in [START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF], it can be verified that h c satisfies Hypothesis 1. We will also assume that Eq. (2.6) holds. Indeed, verifying a priori that the components of h c are algebraically independent is very difficult for the considered problem. This hypothesis will indirectly be verified by checking that the solution λ sol is well identified (see Proposition 1-(c) and Section 5.9).

Error function and convergence analysis of the sequence of MCMC generator

The error function is defined by Eqs. (3.58) and (3.59). The constraint on the random normalized residue (see Section 5.8-(i)) is introduced with b c ρ = 1 in order to avoid the increasing of E{ h c ρ (H c )} with respect to E{ h c ρ (H λ 1 )} (first iteration without constraints effects). This constraint, which is taken into account in Algorithm 1 for computing {λ i , i = 1, . . . , i max }, is not taken into account in the error function (see Eq. (3.58)) to identify the index i sol (see Eq. (3.57)) of the optimal value λ sol = λ i sol of λ i . Consequently, Eq. (3.58) is written, for i ∈ {1, . . . , i max },

err(i) = (err C (i) / err C (1)) 2 + (err δ (i) / err δ (1)) 2 , (5.33) err C (i) = b c C -E{ h c C (H λ i )} / b c C , err δ (i) = | b c δ -E{ h c δ (H λ i )} | / b c δ .
(5.34)

For the three cases, SC1, SC2,and SC3, Fig. 2 displays the error function i → err(i) defined by Eq. (5.33), computed with the constrained learned set D H λ i for N = 1000, 2000, 6000, and 10 000. It can be seen that convergence is reached for N = 10 000 and that, at convergence, function i → err(i) is relatively smooth (that is not the case for N = 1000). These graphs show a good illustration of the convergence of the sequence in N of the MCMC generator using the statistical surrogate model ĥN of h c (see Proposition 4). When the convergence is reached for N = 10 000, Table 3 gives the value of i sol such that λ sol = λ i sol (see Eq. is the pdf of the random variable ∆ eff 2,c defined by Eq. (5.18), at convergence λ = λ sol (introduction of subscript "c"). For the three cases, SC1, SC2, and SC3, Fig. 3 (5.35). Fig. 5 shows the graph of the posterior pdf r → p ρ c (r) of the random normalized residue ρ c , estimated with the constrained learned set for N = 10 000 and its counterpart for the estimation performed with the training set (constructed using the prior probability model). (iv) For cases SC1, SC2, and SC3, Table 4 gives the posterior statistics computed with the constrained learned set for N = 10 000 (subscript "c"), the prior statistics computed with the training set (subscript "d"), and the targets (superscript "exp"), for the second-order moment of the random normalized residue, for the Frobenius norm of the mean value of the random effective/apparent elasticity matrix, and for the coefficient of dispersion of this random matrix. For the same three cases, Table 5 gives the values of the entries of the mean matrices [ C eff d ] computed with the training set, [ C eff c ] computed with the constrained learned set for N = 10 000, and [ C exp ] for the targets. Note that entries (4, 5), [START_REF] Congdon | Bayesian Statistical Modelling[END_REF][START_REF] Gentle | Computational Statistics[END_REF], and [START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Gentle | Computational Statistics[END_REF], which are small with respect to the other entries, are not given. (thick blue line) estimated with the constrained learned set for N = 10 000 and its prior counterpart (thin black line) corresponding to an estimation with the training set (constructed using the prior probability model).

Posterior probability model of parameters

The prior probability model concerns the R n g -valued random variable G that corresponds to the spatial discretization of the R 21 -valued random field G, and the R 3 -valued random variable W that is related (see Eq. , estimated with the constrained learned set for N = 10 000, and their prior counterparts estimated with the training set constructed using the prior probability model. It should be noted that for each one of the random variables C bulk , C shear , and δ C , its prior probability model is the same for the three cases and consequently, does not depend on the case contrary to its posterior probability model that depends on it.

5.12. Discussion about the presented results (i) The results obtained with the posterior model (see Tables 4 and5) show that the constrained learned set significantly improves the prior probability model used for generating the training set. The comparison of the posterior statistics with the targets are good.

(ii) As explained in Section 5.9 the residue constraint is taken into account for the generation of the constrained learned set, but does not intervene in the estimation of the optimal value λ sol of λ using the error function. Nevertheless and as expected, Fig. 5 and Table 4 show that the residue is controlled and stayed small with respect to the reference (the training) for the optimal solution. (iii) The convergence of the sequence of MCMC generators with respect to the number of points generated in the constrained learned set is good as shown by Figs. 0.6702 0.6909 0.6958 0.6785 0.6985 0.7003 0.6835 0.6976 0.7027 [START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF][START_REF] Carlin | Bayesian Methods for Data Analysis[END_REF] 0.6714 0.6903 0.6949 0.6732 0.6950 0.6960 0.6726 0.6872 0.6980 [START_REF] Gentle | Computational Statistics[END_REF][START_REF] Gentle | Computational Statistics[END_REF] 0.6713 0.6924 0.6960 0.6727 0.6958 0.6970 0.6749 0.6933 0.6991 (iv) The dispersion of the target, measured by the value of δ exp , is smaller than the one of the prior probability model (with which the training set has been constructed). So, it was expected that the posterior random effective/apparent elasticity matrix be less dispersed than the one exhibited by the prior probability model. The results confirm this point as it can be seen in Tables 4 and5 and also in Figs. 3 and6.

(v) Comment (ii) is extended as follows. It can be seen that the target is "well" reached for the mean value of the random effective/apparent elasticity matrix (see Table 5), while it is "less well" reached for its coefficient of dispersion computed using the maximum likelihood. This can be explained by the fact that the second-order moment E{ρ 2 c } of the random normalized residue of the equation is kept small during the probabilistic learning of the posterior probability measure, what prevents the dispersion coefficient from reaching its target. We have carried out numerical tests without imposing the constraint related to E{ρ 2 c }, which should stay close to 1. We have observed that the dispersion coefficient "reasonably" reached its target but that E{ρ 2 c } did not stay close to 1 but took on significant values greater than 1. There is indeed a choice of objective between (1) correctly satisfying the constraint on the dispersion coefficient while degrading the value of E{ρ 2 c } or (2) preserving a small value of E{ρ 2 c } to the detriment of perfectly reaching the target for the dispersion coefficient. We have chosen to present the compromise consisting in taking into account the constraint on the residue during the probabilistic learning process, but the error function that we have chosen to identify the optimal value λ sol of λ does not take into account the constraint on the residue. (vi) A last comment concerns the effects of no scale separation. As expected, for the three cases SC1, SC2, and SC3, Table 4 shows that the coefficient of dispersion is significant and increases with the spatial correlation lengths of the random apparent elasticity field at mesoscale, inducing statistical fluctuations of the effective/apparent elasticity tensor at macroscale. It should be noted that, even for the case SC1, for which homogenization in the plane of the plate (domain Ω) is guaranteed (the correlation lengths L C1 and L C2 being much lower than 1), this is not the case for the correlation length L C3 that is equal to the thickness of the plate. Consequently, there is no homogenization at the macroscopic scale and the effective/apparent elasticity tensor remains random and is not deterministic.

Conclusions

In this paper, we have presented general methodology to estimate a posterior probability model for a stochastic value problem from a prior probability model. The given targets are statistical moments for which the underlying realizations are not available. Under these conditions, it has been proposed to use the Kullback-Leibler divergence minimum principle for estimating the posterior probability measure, given the prior probability measure and the constraints related to the targets of the statistical moments. We have proposed the construction of a statistical surrogate model of the implicit mapping that represents the constraints. The constrained learned set, which defines the posterior model, is constructed using only a training set constituted of a small number of points. A mathematical analysis of the proposed methodology has been presented. We have defined the required mathematical hypotheses, which have allowed us to prove the convergence of introduced approximations. We have also given all the necessary numerical elements, which facilitate the implementation of the methodology in a parallel computing framework.

The application presented to illustrate the proposed theory is also, as such, a contribution to the three-dimensional stochastic homogenization of heterogeneous linear elastic media in the case of a non-separation of the microscale and macroscale, that is to say, when there are significant statistical fluctuations in the effective/apparent elasticity tensor at macroscale. The prior stochastic model of the elasticity tensor field at the mesoscopic scale is an advanced model, recently proposed, which takes into account uncertainties on its spectral measure. In addition to the statistical moments of the random effective/apparent elasticity tensor, for the construction of the posterior probability measure by probabilistic learning inference, the second-order moment of the random normalized residue of the stochastic partial differential equation has been added as a constraint. This constraint guarantees that the algorithm seeks to bring the statistical moments closer to their targets while preserving a small residue. The results obtained are those which were expected and give a very good illustration of the theory developed for a non-trivial application.

Definition 1 (

 1 Training set D d and probability measure of H). Let ν and N d be integers such that N d > ν. Let D d be the set of N d points η 1 d

(3. 10 )

 10 Lemma 1 (Convexity of C ad,λ and integrability properties). Under Hypothesis 1 and with Definition 2, (a) C ad,λ defined by Eq. (3.9) is a convex open subset of R n c . (b) ∀λ ∈ C ad,λ , let η → V λ (η) : R ν → R be the function defined by

  which allow for proving the integrability at infinity and then proving Eq.(3.15). Proposition 1 (Construction of the probability measure of H λ ). We consider Hypothesis 1 and Definition 2. For all λ in C ad,λ , let p H λ (η ; λ) = c 0 (λ) ζ(λ) exp{-λ , h c (η) } (3.16) be the pdf of H λ (see Eq. (2.5)) with c 0 (λ) satisfying Eq. (3.14)). (a) The R n c -valued random variable h c (H λ ) is of second-order, E{ h c (H λ ) 2 } < +∞ .

  ) has a solution λ that belongs to C ad,λ , then this solution is unique and we have λ sol = λ. The pdf p H c of H c , which satisfies the constraint E{h c (H c )} = b c is written (see Eq. (3.13) or (3.16)) as p H c (η) = p H λ sol (η ; λ sol ) , ∀η ∈ R ν . (3.23) PROOF. (Proposition 1). (a) Using Eq. (3.13), Eq. (3.14), and the first equation Eq. (3.15) yield

  which using Eq. (3.34) yields Eq. (3.27). Proposition 2 (Convergence of sequences { ĥN (η ; λ)} N and {[∇ η

. 39 )

 39 (a) The initial condition u 0 ∈ R ν is chosen from the points of the training set D d = {η 1 d , . . . , η N d d } (see Section 4.

. 58 )

 58 in which {w m ≥ 0 , m = 1, . . . M} are given real numbers and where err m

  the Störmer-Verlet scheme applied to realization θ of Eqs. (3.44) to (3.46) yields the following recurrence,

Algorithm 1 for = 1 :: end for 10 :

 1110 Algorithm for calculating λ sol and D H c = {η 1 c , . . . , η N c }. 1: Data:N d , D d = {η 1 d , . . . , η N d d }, N, i max , M s , t s , ∆t, f 0 , γ = f 0 ∆t/4 2: Init: ∆W m+1 , ∈ {1, . . . , N}, m ∈ {1, . . . , M s -1}, u 0 and v 0 for ∈ {1, . . . , N}, λ 0 = 0 n c 3: for i = 1 : i max do 4: for = 1 : N (loop in parallel computation) do 5:D H λ i = {η 1 λ i , . . . ,η N λ i } from Eq. (4.7), using Eqs. (4.1) to (4.4) and D H λ i-1 (D H λ 0 not used for i = 1) N (loop in parallel computation) do 8: h c (η λ i ), = 1, . . . , N using the BVP 9∇ λ Γ(λ i ) and [Γ (λ i )] using Eqs. (3.19) to (3.20) and D H λ i 11: err(i) using Eqs. (3.58) with (3.59) 12:

  ) using Eq. (3.56) with a relaxation factor α relax ∈]0 , 1] 13:

( i )

 i Isotropic nominal model at mesoscale. At mesoscale, the nominal model is a linear, elastic, homogeneous, isotropic medium whose elasticity tensor C depends only on the Young modulus E = 1.7 × 10 11 N/m 2 and Poisson coefficient ν P = 0.24. The corresponding bulk modulus C bulk = E /(3(1 -2ν)) and shear modulus C shear = E /(2(1 + ν)) are 1.08974 × 10 11 N/m 2 and 6.85484 × 10 10 N/m 2 .

L

  c1 L c2 L c3 SC1 0.1 0.1 0.1 SC2 0.3 0.3 0.1 SC3 0.5 0.5 0.2 (iii) Nominal value z of the control parameter. It is defined by z 1 = C bulk , z 2 = C shear , and z 3 = δ C .

  . The Monte Carlo numerical simulation method is used with N d = 50 independent realizations and the prior probability model of G and W. We then obtain the training set {x j d , j = 1, . . . , N d } with x j d = ( {y mr, j d , 1 ≤ m ≤ r ≤ 3}, g j d , w j d ) in which y mr, j d ∈ R n y is the solution (see Eq. (5.9)) ofN mr (y mr, j d , g j d , w j d ) = 0 n y ,(5.19)that is to say of the linear equation [a mr ( g j d , w j d )] y mr, j d = b mr ( g j d , w j d ) (see Section 5.4). The statistical moments defined in Section 5.5 can then be computed using theN d independent realizations { [C eff, j d ], j = 1, . . . , N d } of the M + 6 -valued random matrix [C effd ] defined by Eq. (5.10) (a subscript "d" is introduced to designate the computation done with the training set, based on the prior probability model). Since N d n x , the normalization of X is performed using the principal component analysis as follows. Letx j ctr = x j d -x with x = (1/N d ) N d j=1 x j d ∈ R n x . Let [x ctr ] = [x 1 ctr . . . x N d ctr ] be the matrix in M n x ,N d and let [Φ] [S ] [Ψ] T = [x ctr ] be the thin SVD [151] (economy size SVD) of matrix [x ctr ]. The diagonal entries of diagonal matrix [S ] are the singular values S 1 ≥ . . . ≥ S N d -1 > S N d = 0 that are in decreasing order and we have S N d = 0. The matrix [Φ] is in M n x ,ν with ν = N d -1 and [Φ] T [Φ] = [I ν ]. Then random vector X can be written as

Figure 1 :

 1 Figure 1: For cases, SC1, SC2, and SC3, distribution of the eigenvalues κ α .

  Therefore, using the notation of Eq. (3.55), we haveE{ h c ρ (H λ i )} = b cρ , in which we choose b c ρ = 1 (this value is close to the value of the second-order moment of the random normalized residue without constraint) and where η → h c ρ (η) : R ν → R is a positive-valued implicit function.(ii) The second constraint is given by Eq. (5.13). Transforming the upper triangular matrix of [C eff n ] ∈ M + 6 in a R21 vector, Eq. (5.13) is rewritten, similarly to Eq. (3.55), asE{ h c C (H λ i )} = b c C , in which η → h c C (η) : R ν → R21 is an implicit function and where b c C ∈ R21 is the reshaping of the upper triangular matrix of [ C exp n ] ∈ M + 6 defined by Eq. (5.12).

  (3.57)) and the corresponding value err(i sol ) of the error.

Figure 2 :

 2 Figure 2: For scale separation, SC1 (left figure), SC2 (central figure), and SC3 (right figure), graph of error function i → err(i) for N = 1000, 2000, 6000, and 10 000.

  displays the graph of the posterior pdf δ 2 → p ∆ eff 2,C (δ2) of random variable ∆ eff 2,c , estimated with the constrained learned set for N = 10 000 and its prior counterpart estimated with the training set (constructed using the prior probability model). Fig. 4 (left figure) displays the graph of the Frobenius norm N → [ C eff c (N)] F of the posterior mean value [ C eff c (N)] = E{ [C eff c (N)]} of the random effective/apparent elasticity matrix [C eff c (N)] as a function of N and estimated using the constrained learned set, while Fig. 4 (right figure) shows the graph of the maximum likelihood N → δ eff c,ML (N) of the coefficient of dispersion of [C eff c (N)] defined by Eq.

Figure 3 :

 3 Figure 3: For scale separation, SC1 (left figure), SC2 (central figure), and SC3 (right figure), graph of the posterior pdf of random variable ∆ eff 2,c

Figure 4 :

 4 Figure 4: For scale separation, SC1, SC2, and SC3, graph of N → [ C eff c (N)] F (left figure) and N → δ eff c,ML (N) (right figure), computed with the constrained learned set.

Figure 5 :

 5 Figure 5: For scale separation, SC1 (left figure), SC2 (central figure), and SC3 (right figure), graph of the posterior pdf of random normalized residue ρ c (solid line) estimated with the constrained learned set for N = 10 000 and its prior counterpart (dashed line) corresponding to an estimation with the training set (constructed using the prior probability model).

  2 and 4 in accordance to Proposition 4.

Figure 6 :

 6 Figure 6: Graphs of the posterior pdf c → p C bulk (c) (left figure), c → p C shear (c) (central figure), and δ → p δ C (δ) (right figure), estimated with the constrained learned set for N = 10 000, for cases, SC1 (thin solid line), SC2 (med solid line), and SC3 (thick solid line), and the corresponding prior pdf estimated with the training set constructed using the prior probability model (dashed line).

  H ) is not a distance because the symmetry property and the triangle inequality are not verified. It can easily be seen that the cross entropy S (p, p H ) = -R ν p(η) log(p H (η)) dη and the entropyS (p) = -R ν p(η) log(p(η)) dη are related to D(p, p H ) by S (p, p H ) = S (p) + D(p, p H ).Finally, it can be proven (see for instance[START_REF] Cover | Elements of Information Theory[END_REF]) that (p, p H ) → D(p, p H ) is a convex function in the pair (p, p H ). The probability density function p H c on R ν , which satisfies the constraint defined by Eq. (2.2) and which is closest to p H , is the solution of the optimization problem (see for instance[START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF][START_REF] Cover | Elements of Information Theory[END_REF][START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF]),

is expressed as a function of H c . We recall that the Kullback-Leibler divergence between two probability measures p(η) dη and p H (η) dη on R ν is defined by D(p, p H ) = R ν p(η) log(p(η)/p H (η)) dη, and is such that D(p, p H ) ≥ 0 (that can be proven by applying the Jensen inequality

[START_REF] Jensen | Sur les fonctions convexes et les inégalités entre les valeurs moyennes[END_REF][START_REF] Durrett | Probability, Theory and Examples[END_REF]

) and D(p, p H ) = 0 if and only if p = p H . Note that (p, p H ) → D(p, p

  and where . F is the Frobenius norm. Definition 2 (Admissible subset C ad,λ of R n c ). Under Hypothesis 1, the admissible set C ad,λ of Lagrange multiplier λ is defined as the open subset of R n c such that

  valued random variables. Since λ and λ are in C ad,λ , we have E{A} < +∞ and E{B} < +∞. We have to prove that E{ exp{-µλ + (1 -µ)λ , h c (H) } } < +∞, that is to say, E{A

µ B 1-µ } < +∞. For µ = 0 and µ = 1, it is verified. For µ ∈]0, 1[, since A and B are almost-surely positive and using the Hölder inequality yield E{A µ B 1-µ } ≤ (E{A}) µ × (E{B}) 1-µ < +∞, which finishes the proof of the convexity of C ad,λ . (b) Using Eqs. (3.2), (3.4), and (3.11), yields

  be the probability measure of H λ , for which the density p H λ is defined by Eq. (3.13) with supp p H λ = R ν . Eq. (3.31) can then be rewritten, for all η in R ν , as

  .[START_REF] Lio | Uncertain maximum likelihood estimation with application to uncertain regression analysis[END_REF] PROOF. (Proposition 2). The probability measure P H λ (dη ; λ) admits a continuous density p H λ ( . ; λ) with respect to dη on R ν (see Eq. (3.13)). Using the notation of Definition 3, for s

  yields Eq. (3.60) by using Eq. (3.67).4.A few numerical elements for implementation of the methodology 4.1. Choice of the integration scheme for solving the ISDE introduced in Propositions 3 and 4As we have previously explained, for λ ∈ C ad,λ , the ISDE defined by Eqs. (3.44) to(3.46), must be solved for t ∈ [0 , t s ] (see Proposition 3-(e)) with the initial condition at t = 0 defined in Proposition 3-(a), in order to generate the constrained learned set

Table 1 :

 1 Values of the spatial correlation lengths L c1 , L c2 , and L c3 for cases SC1, SC2, and SC3 of scale separation.

Table 3 :

 3 For cases SC1, SC2, and SC3 of scale separation, value err(i sol ) of the error function for the solution λ sol = λ i sol computed with the constrained learned set for N = 10 000. Residue and posterior second-order moments of the random effective/apparent elasticity matrix estimated with the constrained learned set The posterior statistics of the random normalized residue and the random effective/apparent elasticity matrix are estimated with the constrained learned set D H c = D H λ sol as a function of N for the three cases, SC1, SC2, and SC3. (i) The posterior second-order moment of the random normalized residue is E{ ρ 2 c } and compared to 1, see Section 5.8-(i). (ii) The posterior mean value is [ C eff c ] = E{[C eff c ]} and is compared to [ C exp ] (see Eq. (5.13) for the normalized version). We also consider the Frobenius norm [ C eff c ] F that is compared to [ C exp ] F . (iii) For the posterior dispersion coefficient of [C eff c ], instead of comparing δ eff c to δ exp (see Eq. (5.16)), we use the maximum likelihood consisting in comparing δ eff c,ML to δ exp in which

	Case	SC1	SC2	SC3
	i max	30	30	20
	i sol	27	30	17
	err(i sol ) 0.5924 0.5495 0.5364

Table 4 :

 4 For cases SC1, SC2, and SC3, posterior statistics computed with the constrained learned set for N = 10 000 (subscript "c"), prior statistics (subscript "d") computed with the training set, and targets. ×1011 4.2106 4.1925 4.1943 [ C eff c ] F ×10 11 4.6294 4.6923 4.6816 [ C exp ] F ×10 11 4.6317 4.6549 4.6706 δ eff

		SC1	SC2	SC3
	E{ρ 2 c } b c ρ	1.2938 1.2687 1.2413 1 1 1
	[ C eff d ] F d,ML	0.2257 0.2469 0.2701
	δ eff c,ML	0.1329 0.1476 0.1671
	δ exp	0.0946 0.1374 0.1825

Table 5 :

 5 For cases SC1, SC2, and SC3, values of [ C eff d ] i j computed with the training set, [ C eff c ] i j computed with the constrained learned set for N = 10 000, and [ C exp ] i j for the targets.

			SC1			SC2			SC3
	Entries of
	(6 × 6)	[ C eff d ]	[ C eff c ]	[ C exp ]	[ C eff d ]	[ C eff c ]	[ C exp ]	[ C eff d ]	[ C eff c ]	[ C exp ]
	matrix
	(1, 1)	2.0904 2.2600 2.2652 2.0792 2.2914 2.2810 2.0465 2.2751 2.2946
	(1, 2)	0.7427 0.8809 0.8753 0.7269 0.8982 0.8809 0.7140 0.8874 0.8824
	(1, 3)	0.7458 0.8804 0.8745 0.7324 0.8983 0.8800 0.7381 0.8999 0.8826
	(2, 2)	2.0832 2.2603 2.2668 2.0786 2.2946 2.2846 2.0917 2.2830 2.2822
	(2, 3)	0.7451 0.8802 0.8734 0.7486 0.8950 0.8754 0.7471 0.9045 0.8808
	(3, 3)	2.0839 2.2647 2.2680 2.0777 2.2841 2.2697 2.1038 2.2928 2.2812
	(4, 4)

Stochastic computational model and random effective/apparent elasticity matrix

The finite element method is used for discretizing the weak formulation of the stochastic BVP. The finite element mesh of domain Ω is made up of 60 × 60 × 6 = 21 600 eight-nodes solid elements, 26 047 nodes, and 78 141 dofs (25 926 Dirichlet conditions on ∂Ω and n y = 52 215 the remaining dofs). There are 2 3 integration points in each finite element, which yields n p = 172 800 integration points for the spatial discretization of the random fourth-order tensorvalued elasticity field {C(ξ), ξ ∈ Ω}. The discretization of random field C is expressed as a function of a R n g -valued random variable G corresponding to the spatial discretization of random field G with n g = 21 × n p = 3 628 800. For 1 ≤ m ≤ r ≤ 3, let Y mr be the R n y -valued random variable of the free dofs (that corresponds to the dofs of the nodes inside Ω of the finite element discretization of random field {Y mr (ξ), ξ ∈ Ω}). Therefore, the stochastic computational model can be written as the R n y -valued stochastic equation,

which is a stochastic linear equation that can be rewritten as [a mr (G, W)] Y mrb mr (G, W) = 0 n y in which, for g d ∈ R n g and w d ∈ R n w , [a mr ( g d , w d )] is a matrix in M + n y (thus, invertible) and where b mr ( g d , w d ) is a vector in R n y , which depends on the Dirichlet condition defined by Eq. (5.2). For 1 ≤ m ≤ r ≤ 3, the finite element discretization of the right-hand side member of Eq. (5.4) yields the M + 6 -valued random effective/apparent elasticity matrix [C eff ] such that [C eff ] ij = C eff i jmr in which the indices i = (i, j) with 1 ≤ i ≤ j ≤ 3 and j = (m, r) with 1 ≤ m ≤ r ≤ 3 are with values in {1, . . . , 6}. This random matrix can be written as

(5.10)

Note that, for 1 ≤ m ≤ r ≤ 3, Y mr satisfying Eq. (5.9) corresponds to the strong stochastic solution of the finite element discretization of the weak formulation of the stochastic BVP (see Section 5.3). Due to Eqs. (5.7) and (5.8), we have,

which proves that Y mr and [C eff ] are second-order random variables.

Definition of the statistical moments and their targets

From Eq. (5.11), the M + 6 -valued random variable [C eff ] is of second-order. We can then define its first two moments. (i) The first statistical moment of interest is the mean value [

. Introducing the subscript "n" to designate a normalization, we define the normalized quantities with respect to µ exp as,

(5.12)

The corresponding constraint equation will then be written as

(ii) The second statistical moment of interest is the coefficient of dispersion δ eff of random matrix [C eff ] and its target counterpart δ exp . Let ∆ eff 2 be the positive-valued random variable defined by

Consequently, δ eff that is defined by

can be rewritten as

Defining µ eff = [ C eff ] F , the constraint equation that is defined by

can be rewritten, using Eqs. (5.12) to (5.15), as

(5.17) while Eq. (5.14) yields

It should be noted that, if δ eff goes to zero, then the statistical fluctuations of [C eff ] goes to zero because, due to the Tchebychev inequality, [C eff ] goes in probability to its mean value [ C eff ] (this would be the case of a scale separation).

(iii) As we have explained and as we will see in Sections 5.7 and 5.8, the constrained learned set is constructed by also taking as a constraint the second-order moment of the normalized random residue of the stochastic PDE. On the other hand, as shown by equations Eqs. (5.19) and (5.29), the spatial discretization G of the random field {G(ξ), ξ ∈ Ω} is obviously a fundamental quantity that intervenes in the calculation of this normalized random residue during the constrained learning process. Furthermore, it is essential that the mathematical properties of {C(ξ), ξ ∈ Ω} (random field with values in the positive-definite symmetric fourth-order tensors, lower bound, statistical homogeneity, control of the random spectral measure), defined by Eq. (5.5), remain verified. This implies strong mathematical properties on the random field {G(ξ), ξ ∈ Ω}. This class of random fields with a random spectral measure, which was constructed and analyzed in [START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF] and that we are using in this paper, already corresponds to a large class, given all the algebraic and probabilistic properties that this class must verify. The targets have been simulated by using the stochastic computational model while keeping the same class of random fields for {C(ξ) ∈ Ω} but for which the values of its hyperparameters are different of those used for the prior probabilistic model. It should be noted that changing the class by keeping this type of algebraic representation would be very difficult. A way for performing an extension of the class would then be to use a polynomial Chaos representations as proposed in [START_REF] Nouy | Random field representations for stochastic elliptic boundary value problems and statistical inverse problems[END_REF][START_REF] Soize | Uncertainty Quantification[END_REF] to properly take into account all the required properties. As the objective of this paper is not to build an even larger class for generating the targets, we have preferred to keep the same class. Finally, note that this choice does not facilitate the analysis, because the methodology and its proven mathematical properties are independent of the choice of the random fields class provide that the class be correctly constructed for the stochastic elliptic operator. Nevertheless, the proximity of the constructed posterior probability measure to the prior probability measure is certainly dependent on the chosen class.

(iv) For the three considered cases of scale separation, the numerical values of µ eff for the training set (and thus denoted by µ eff,d ) and µ exp for the target are given in Table 2. The stochastic computational model defined in Section 5.4 is used for generating the training set related to the random variable X = ( {Y mr , 1 ≤ m ≤ r ≤ 3}, G, W) with values in R n x = R 6×n y × R n g × R n w with n x = 6 n y + n g + n w =