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

Probabilistic learning constrained by realizations using a weak formulation of Fourier transform of probability measures

Keywords: Probabilistic learning, realizations as targets, statistical inverse problem, Kullback-Leibler divergence, uncertainty quantification

This paper deals with the taking into account a given target set of realizations as constraints in the Kullback-Leibler divergence minimum principle (KLDMP). We present a novel probabilistic learning algorithm that makes it possible to use the KLDMP when the constraints are not defined by a target set of statistical moments for the quantity of interest (QoI) of an uncertain/stochastic computational model, but are defined by a target set of realizations for the QoI for which the statistical moments associated with these realizations are not or cannot be estimated. The method consists in defining a functional constraint, as the equality of the Fourier transforms of the posterior probability measure and the target probability measure, and in constructing a finite representation of the weak formulation of this functional constraint. The proposed approach allows for estimating the posterior probability measure of the QoI (unsupervised case) or of the posterior joint probability measure of the QoI with the control parameter (supervised case). The existence and the uniqueness of the posterior probability measure is analyzed for the two cases. The numerical aspects are detailed in order to facilitate the implementation of the proposed method. The presented application in high dimension demonstrates the efficiency and the robustness of the proposed algorithm.

Introduction

The Bayesian inference is a fundamental method of computational statistics, which permits the effective integration of data (target set) into a predictive uncertain/stochastic computational model (SCM) in physics and engineering sciences. This target set is constituted of realizations/samples of a quantity of interest (QoI) of the SCM and allows thus the posterior probability measure to be estimated. Concerning the Bayesian updating, one can refer the reader, for example, to [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] for general aspects, [START_REF] Kaipio | [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] 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] Picchini | Coupling stochastic em and approximate Bayesian computation for parameter inference in state-space models[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, or [START_REF] Fearnhead | Exact and efficient Bayesian inference for multiple changepoint problems[END_REF] for Bayesian inference for changepoint problems. 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]. It should be noted that Bayesian inference can be tricky to use [START_REF] Owhadi | On the brittleness of Bayesian inference[END_REF], in particular for the high dimension.

On the other hand, when the target set does not consist of realizations, but is made up of statistical moments of the QoI, the Bayesian inference cannot easily be used and the Kullback-Leibler divergence minimum principle (KLDMP) [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][START_REF] Gentle | Computational statistics[END_REF] allows for estimating the posterior probability measure given its prior probability measure and the constraints related to the statistical moments. This principle has extensively been used over the last three decades for imposing constraints in the framework of learning with statistical models (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]).

As the statistical moments cannot easily be taken into account with the Bayesian approach, a natural way (or at least the one that is proposed in this paper) is that of the use of the KLDMP. However, there is a difficulty in this case, because the constraints associated with the target set must be expressed as mathematical expectations, which is the case for the target statistical moments, but which is not the case for the target realizations. It is this difficulty that we propose to solve by presenting a method to take into account target realizations within the framework of the KLDMP. We propose a novel functional approach that consists in constructing and analyzing a weak formulation of the Fourier transform of the probability measure and to derive from it a finite representation of the functional constraint. This functional constraint is then expressed as a mathematical expectation that is thus adapted to the KLDMP formulation and that depends only on the target realizations. We then obtain a novel probabilistic learning algorithm that allows for integrating a target set of realizations of the QoI into a stochastic computational model. It should be noted that the constraint defined by the target realizations, which is imposed by this novel functional representation, can also be seen as a method that can replace the constraints defined by the statistical moments associated with these target realizations, statistical moments that are not estimated for the reasons above (see also Framework 2 below in which additional explanations are given).

Two main frameworks can a priori be identified to use the method presented in this paper. Framework 1. For the first one, the random vector of the training set is in high dimension and the target set is simultaneously made up of (i) a relatively small number of target realizations for certain components of the QoI and (ii) of a target set of statistical moments (also called target statistical moments) for other components. In this paper, it is assumed that the statistical moments associated with the target realizations are not or cannot be well estimated (see Framework 2 below). For instance, the following case can be considered.

(i) Concerning the target realizations whose number is relatively small and knowing that the QoI and the control parameters are in high dimension, one does not want to degrade the information by replacing the constraint defined by the target realizations with associated estimated statistical moments that are not sufficiently converged (note that the high-order statistical moments that allow for better representing the probability measure cannot be estimated). These target realizations can come from experimental measurements or from numerical simulations carried out using a high-fidelity uncertain/stochastic computational model.

(ii) Concerning the target statistical moments, the associated realizations are not available because these realizations have been destroyed or lost, but the target statistical moments were published or can be found in scientific/technical reports (this is often the case for some experimental data relative to complex systems carried out in past decades). These target statistical moments can also correspond to the mean-square minimization of the random residue of the stochastic equation of the computational model, which is imposed as constraint for estimating the posterior probability measure in the KLDMP-based learning procedure (see for instance [START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]). Framework 2. For the second one, as previously sketched, the method presented is an alternative to the moment problem when the underlying realizations are known and when we want to best represent a non-Gaussian posterior probability measure. In such a case, several statistical moments have to be considered such as the mean vector, the second-, third-, and fourth-order tensor-valued moments (and possibly some higher-order tensors) for a non-Gaussian random vector that is in high dimension. This is not easy and, as previously mentioned, a large number of realizations are required to obtain a reasonable convergence of these tensor-valued moments. In addition, these tensor-valued moments must satisfy a condition to be the moments of a probability measure [START_REF] Shohat | The problem of moments[END_REF]. In high dimension and for a set of tensor-valued moments, it is difficult enough. We thus see that imposing realizations of a high-dimensional random vector as a constraint for the KLDMP is an interesting alternative to impose a constraint defined by a set of tensor-valued moments, and the convergence analysis is greatly simplified.

In order not to complicate the presentation of the second part of this paper, dedicated to the KLDMP for estimating the posterior model, we only consider the case for which the target set of realizations is applied as constraints (we do not consider the mixed case whose some components of the QoI are constrained by statistical moments, but there is no difficulty in considering simultaneously the two types of constraints).

1.1. Framework of the considered problem, objectives of the paper, and methodology proposed (i) First case referred as the unsupervised case. The quantity of interest is a R n q -valued random variable Q, defined on a probability space (Θ, T , P), whose prior probability measure is P Q (dq) on R n q . This prior probability measure is unknown but is the underlying probability measure that has been used to generate the training set D d = {q 1 d , . . . , q N d d } constituted of N d independent realizations {q j d ∈ R n q , j = 1, . . . , N d } of Q (the subscript "d" is introduced to reference the "data" of the training set). It is assumed that n q is big (high-dimension problem). For instance, q j d can be the realizations of the discretization of a random field indexed by a bounded part of R d with d ≥ 2. Related to Q, a target set D targ = {q 1 targ , . . . , q N r targ } is given, constituted of N r given points q r targ in R n q , which are N r independent realizations of a R n q -valued random variable Q targ defined on (Θ, T , P), independent of Q, whose probability measure P targ Q of Q targ is assumed to be unknown. Giving the training set D d = {q 1 d , . . . , q N d d } of Q and the target set D targ = {q 1 targ , . . . , q N r targ } of Q targ , the KLDMP will allow for identifying the probability measure P post Q that is closest to P Q (dq) while satisfying the constraint defined by D targ . The measure P post Q , which is the measure updated with the constraint, will be called the posterior probability measure of the R n q -valued random variable Q post defined on (Θ, T , P). The probabilistic learning thus consists in using a MCMC algorithm for generating N realizations {q ℓ post , ℓ = 1, . . . , N} of Q post . Regarding the resampling of a probability measure with MCMC algorithms, it should also be noted that, when the available training set is composed of a small number of points, suitable algorithms should be used like those which have been specifically developed to deal with the case of small data (see [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Perrin | Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints[END_REF][START_REF] Farhat | Feasible probabilistic learning method for model-form uncertainty quantification in vibration analysis[END_REF][START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF][START_REF] Guilleminot | Data-driven enhancement of fracture paths in random composites[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Arnst | Computation of sobol indices in global sensitivity analysis from small data sets by probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning on manifolds (plom) with partition[END_REF] for data-driven problems and [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF][START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF][START_REF] Capiez-Lernout | Nonlinear stochastic dynamics of detuned bladed disks with uncertain mistuning and detuning optimization using a probabilistic machine learning tool[END_REF] for optimization problems).

(ii) Second case referred as the supervised case. The quantity of interest is the above R n q -valued random variable Q and there is a control parameter that is a R n w -valued random variable W. The random variables Q and W are defined on the probability space (Θ, T , P), whose prior joint probability measure is P Q,W (dq, dw) on R n q × R n w . As for the unsupervised case, this prior joint probability measure is unknown but is the underlying probability measure that has been used to generate the training set

D d = {x 1 d , . . . , x N d d } constituted of N d independent realizations, {x j d = (q j d , w j d ), j = 1, . . . , N d } of the R n x -valued random variable X = (Q, W) with n x = n q + n w . The probability measure of X is P X (dx) = P Q,W (dq, dw).
It is assumed that n q and n w are big (high dimension problem). This supervised case can correspond to Q = f (W) in which f is an unknown measurable mapping from R n w into R n q or to Q = f (W, U) in which f is also an unknown measurable mapping from R n w × R n u into R n q and where U is an uncontrolled R νvalued random variable defined on (Θ, T , P). In the first case, q j d = f (w j d ) and P Q,W (dq, dw) has no density with respect to dq ⊗ dw, while in the second case, q j d = f (w j d , u j d ) and P Q,W (dq, dw) can have a density. As for the unsupervised case, we consider a given target set D targ = {q 1 targ , . . . , q N r targ } for the quantity of interest, constituted of N r independent realizations of the R n q -valued random variable Q targ that is independent of Q. Note that no target realization is given for the control variable W. If we gave target realizations for W, which would amount to giving ourselves a target set of realizations for X, then in terms of the methodology presented in this paper, we would be in a situation similar to that of the unsupervised case. In the supervised case that we consider here, the considered system is under-observed with respect to the given target set of realizations. Similarly to the unsupervised case, giving the training set D d = {x 1 d , . . . , x N d d } of X and the target set D targ = {q 1 targ , . . . , q N r targ } of Q targ , the KLDMP will allow for identifying the probability measure P post X that is closest to P X (dx) while satisfying the constraint defined by D targ . The measure P post X , which is the measure updated with the constraint on Q, will be called the posterior probability measure of the R n q × R n w -valued random variable (Q post , W post ) defined on (Θ, T , P). The probabilistic learning thus consists in using a MCMC algorithm for generating N realizations {x ℓ post , ℓ = 1, . . . , N} of X post , that is to say, N realizations {(q ℓ post , w ℓ post ), ℓ = 1, . . . , N} of (Q post , W post ).

Novelty of the paper

In this paper, we present a novel method that makes it possible to use the KLDMP when the constraints are not defined by a target set of statistical moments for the QoI, but are defined by a target set of realizations for the QoI for which the statistical moments associated with these realizations are not estimated. We obtain a novel probabilistic learning algorithm that allows for integrating a target set of realizations of the QoI into an uncertain/stochastic computational model. The proposed methodology can be used for both unsupervised and supervised cases. The functional constraint, defined by the target set of realizations, is formulated in a Hilbertian framework as the equality of the Fourier transforms of the posterior probability measure and the target probability measure. We present the construction and the mathematical analysis of a finite representation of the weak formulation of this functional constraint. We prove an important property of the finite representation, which allows for analyzing the existence and the uniqueness of the posterior probability measure.

Organization of the paper

All the developments given in this paper will be presented within the framework of the supervised case. Given the proposed approach, the unsupervised case follows immediately. This paper is organized in three parts.

The first part (Sections 2 to 5) is devoted to the formulation and the construction of a finite representation of the functional constraint. Section 2 deals with the scaling and the reduced representation of random vector X for which the realizations are the points of the training set, which allows for constructing a normalized random variable H with values in R ν with ν < N d ≤ n x . In Section 3, Definition 1 and Lemma 1 defined the functional constraint as an equality of the Fourier transform of the probability measures of Q and Q targ and then of the random variable H and H targ , in which H targ is the "projection" of Q targ on the model. Section 4 is devoted to the weak formulation of the functional constraint imposed to random variable H giving realizations of H targ .

It should be noted that the explicitness of such a functional constraint requires to sample the probability measures in R ν , what is not possible in high dimension (ν ≫ 1). If such a sampling method was used, then the number of constraints that should be considered in the KLDMP would be huge or even unrealistic. In order to circumvent this difficulty, a weak formulation of this functional constraint is constructed using the fundamental properties of the Fourier transform of the probability measures. Then a finite representation of the weak formulation of this functional constraint is obtained by projection on a vector subspace spanned by a countable family of functions that is explicitly constructed.We thus obtain a finite representation of the functional constraint whose dimension is the number N r of target realizations, which is relatively small. Under adapted mathematical hypotheses for covering a large part of applications, Theorem 1 (proven with the help of three Lemmas) gives the required mathematical results that are necessary to construct the weak formulation (Definition 2) of the functional constraint defined on the space

H 1 = L 1 (R ν , C)∩L 2 (R ν , C).
In Section 5, we present the construction and the analysis of a finite representation of the functional constraint derived from the weak formulation, which is restricted to a Hilbert space H 1, µ that is a subset of H 1 in which µ = p ν (v) dv is a Gaussian probability measure on R ν . Theorem 2 studies the Fourier transform φ of a function φ in H 1, µ , which is a C-valued analytic function on R ν and which belongs to

H 0 = C 0 ∩ L 2 (R ν , C) (in
which C 0 is the space of all the C-valued continuous functions on R ν , which go to zero at infinity). While the considered weak formulation of the functional constraint is posed on a Hilbert space, the Hilbertian structure leads naturally to introduce a Hilbertian basis used to construct a finite representation of the weak formulation of the functional constraint. Given the fact that the measure µ related to H 1, µ is Gaussian, the multi-dimensional Hermite polynomials could be used. However, the multi-index is in high dimension ν and consequently, the curse to dimensionality prevents using this type of finite representation. Based on Theorem 2, Lemmas 5 and 6 give sought construction of the countable family of functions in H 1, µ , which allows for constructing the finite representation that is explicitly described in Definition 4 and that uses the realizations of Q targ (the points of target set D targ ). Lemma 7 gives an important property of the constructed finite representation, which will allow for analyzing the existence and the uniqueness of the posterior probability measure. The first part ends with a numerical illustration of the behavior of the finite representation of functional constraint that is proposed.

The second part of this paper corresponds to Section 6 in which we present the methodology to construct the posterior probability measure based on the use of the KLDMP with the prior model and the target set. This methodology is similar to the one we have used in [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF], but for which the constraints are now the one presented in Section 5. Thus the mathematical proofs are adapted and modified because the hypotheses are no longer the same. The finite representation of the weak formulation of the functional constraint is taken into account by introducing a vector-valued Lagrange multiplier λ. The posterior probability measure is constructed as the limit of a sequence of random variables {H λ } λ indexed by λ. Theorems 3 and 4 give the explicit construction of the probability measure of H λ and its MCMC generator based on the nonlinear stochastic dissipative Hamiltonian system studied in [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF]. This second part ends with the iterative algorithm for computing the optimal value of λ and gives elements for its numerical implementation.

The last part, Section 7, is devoted to a numerical illustration of the supervised case for which the training set D d is made up of N d independent realizations x j = (q j d , w j d ) ∈ R n x = R n q × R n w of random variable X = (Q, W) for which n x = 430 098, n q = 10 098, n w = 420 000, and N d ∈ {100, 200, 300, 400}. The target set D targ is made up of N r independent realizations q r targ ∈ R n q of random variable Q targ for which N r ∈ [50 , N targ ] with N targ ∈ {100, 200, 300, 400}. As we will see, we will also give a lighting on the associated unsupervised case to this supervised case.

Convention for the variables, vectors, and matrices

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.

[X]: boldface upper-case letters between brackets are matrix-valued random variables.

]: identity matrix in M n . ∥ [x] ∥ = sup ∥ y ∥ =1 ∥ [x] y ∥. ∥ [x] ∥ F : Frobenius norm of matrix [x].

Convention used for random variables

In this paper, for any finite integer m ≥ 1, the Euclidean space R m is equipped with the σ-algebra B R m . If Y is a R m -valued random variable defined on the probability space (Θ,

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

Scaling and reduced representation

Before performing the construction of the reduced representation that is performed by using a principal component analysis (PCA) of X, it is assumed that training set D d is scaled using the formulation presented in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. The target set D targ is also scaled using the same transformation that the one used for obtaining the scaled training set D d .

Let x j = x j d -x be the realization of

X with x = (1/N d ) N d j=1 x j d ∈ R n x . Let [ x] = [ x1 . . .
Φ] is in M n x ,N d and [Φ] T [Φ] = [I N d ]. For ν ≤ N d -1, let X (ν) be the representation of X defined by X (ν) = x + [Φ ν ] [κ] 1/2 H , (2.1) in which [Φ ν ] ∈ M n x ,ν is made up of the first ν columns of matrix [Φ], where [κ] is the diagonal matrix in M + ν such that κ α = [κ] αα = S 2 α /(N d -1)
, and where H = (H 1 , . . . , H ν ) is the R ν -valued random variable whose N d independent realizations are Throughout the rest of the paper, in order to simplify the notations, the superscript "(ν)" will be omitted and the random variable X (ν) = (Q (ν) , W (ν) ) will simply be denoted by X = (Q, W). From Eq. (2.1), it can be deduced that

η j d = [κ] -1/2 [Φ ν ] T (x j d -x) , j = 1, . . . , N d . ( 2
d ) = E{∥ X -X (ν) ∥ 2 } E{∥ X ∥ 2 } ≃ 1 - ν α=1 κ α tr{[ C X ]} ≤ ε PCA , ν < N d -1 , (2.3
Q = q + [Φ q ] [κ] 1/2 H , W = w + [Φ w ] [κ] 1/2 H , (2.4)
in which x = (q, w) ∈ R n x = R n q ×R n w , and where [Φ q ] ∈ M n q ,ν and [Φ w ] ∈ M n w ,ν are the corresponding block extraction of [Φ ν ] with respect to Q and W. The training set related to H is

D d = {η 1 d , . . . , η N d d } , η j d ∈ R ν , j = 1, . . . , N d , (2.5)
in which η j d is given by Eq. (2.2). Using D d , the estimates η ∈ R ν and [ C H ] ∈ M + ν of the mean value and the covariance matrix of H are such that

η = 0 ν , [ C H ] = [I ν ] . (2.6)
The first Eq. (2.4) allows for defining the R n q -valued random variable Q such that

Q = Q -q , Q = [Φ q ] [κ] 1/2 H . (2.7)
Note that we have

[Φ q ] T [Φ q ] ̸ = [I ν ].
We also introduce the R n q -valued random variable

Q targ = Q targ -q , (2.8)
whose N r realizations are qr targ = q r targ -q , r ∈ {1, . . . , N r } .

(2.9)

Note that Q targ is generally not centered because q is not the mean value of Q targ , but H is a centered one (see Eq. (2.6)).

Definition of the functional constraint for estimating the posterior probability measure

The objective is to construct the posterior probability measure P post Q (dq) = P post Q,W (dq, R n w ) that is closest to P targ Q (dq), which is equivalent (see Eqs. (2.7) and (2.8)) to construct the posterior probability measure P post Q (dq) that is closest to the probability measure P targ Q (dq) of Q targ for which N r independent realizations { qr targ , r = 1, . . . , N r } are given. For using the KLDMP, we need to express the constraint as a mathematical expectation of a random variable. As previously explained, we propose to define the constraint as the equality of the Fourier transforms of the probability measures (characteristic functions) instead of the probability measures themselves for the following reasons. Such a formulation makes to write the constraint as a mathematical expectation, which is required for the use of the KLDMP. On the other hand, the explicitness of the equality of probability measures on R ν would require to sample the probability measures in R ν , what is not possible in high dimension (ν ≫ 1). This is another reason to consider the equality of the Fourier transform as a functional constraint which makes it possible to construct its weak formulation and then to derive a finite representation from it.

Definition 1 (Constraint defined by the equality of the Fourier transform of the probability measures). The constraint is defined as follows,

∀y ∈ R n q , f(y) = f targ (y) , (3.1)
in which the complex-valued functions f and f targ are the characteristic functions defined on R n q of the R n q -valued random variables Q and Q targ ,

f(y) = E{e i ⟨ y, Q ⟩ } , f targ (y) = E{e i ⟨ y, Q targ ⟩ } , (3.2)
and where i is the imaginary unit and ⟨ •, • ⟩ is the usual Euclidean inner product.

The constraint defined by Eq. (3.1) is in high dimension because y ∈ R n q . We then propose to reduce the dimension using the representation of Q given by Eq. (2.7).

Lemma 1 (Functional constraint using the representation of Q). It is chosen to project the target on the prior model. Therefore, let H targ be the R ν -valued random variable defined by

H targ = [V] T Q targ , (3.3
)

in which the matrix [V] ∈ M n q ,ν is written as [V] = [Φ q ] ([Φ q ] T [Φ q ]) -1 [κ] -1/2
. Using the representation of Q defined by Eq. (2.7), the functional constraint associated with Eq. (3.1) is written as

∀v ∈ R ν , f (v) = f targ (v) , (3.4)
in which the complex-valued functions f and f targ are the characteristic functions, defined on R ν , of the R ν -valued random variables H and H targ ,

f (v) = E{e i ⟨ v, H ⟩ } , f targ (v) = E{e i ⟨ v, H targ ⟩ } . (3.5)
PROOF. (Lemma 1). Using Eq. (2.7), for all y in R n q , we have

⟨ y , Q ⟩ = ⟨ y , [Φ q ] [κ] 1/2 H ⟩ = ⟨ [κ] 1/2 [Φ q ] T y , H ⟩ = ⟨ v , H ⟩ in which v = [κ] 1/2 [Φ q ] T y ∈ R ν .
In the other hand, we perform the projection of Q targ on the prior model.

Since the matrix [Φ

q ] T [Φ q ] ∈ M + ν is invertible and since the diagonal matrix [κ] ∈ M + ν is also invertible, we introduce the pseudo-inverse [V] of [κ] 1/2 [Φ q ] T ) (projection) such that ([κ] 1/2 [Φ q ] T ) [V] = [I ν ]. Therefore, taking y = [V] v for all v in R ν , we have ⟨ y , Q targ ⟩ = ⟨ [V] v , Q targ ⟩ = ⟨ v , [V] T Q targ ⟩ = ⟨ v , H targ ⟩.

2).

Remark 1. There are two difficulties.

(i) The constraint does not concern all the variables, that is to say X = (Q, W), but only the quantity of interest Q. We are therefore in an under-observed case with respect to the applied constraint. This choice of the developments framework imposes to project the target Q targ on the prior model in order to obtain a representation of the target H targ that only depends on Q targ and not on Q targ and W targ because W targ is not given as a constraint.

(ii) As previously explained, the explicitness of the constraint defined by Eq. (3.4) requires to sample v in R ν , what is not easy, not efficient, and not accurate in high dimension (ν ≫ 1). If such a sampling method was used, then the number of constraints that should be considered in the KLDMP would be huge or even unrealistic. In the next section, we thus propose to construct a weak formulation of the functional equation defined by Eq. (3.4) using the fundamental properties of the Fourier transform of the probability measures (see for instance [START_REF] Gelfand | [END_REF]).

Weak formulation of the functional constraint

Notation 1 (Defining C b , C 0 , H 0 and H 1 complex vector spaces). (i) Let C 0 (R ν , C) (resp. L ∞ (R ν , C)) be the complex vector space of continuous (resp. bounded) functions on R ν with values in C. The norm ∥ g ∥ L ∞ in L ∞ (R ν , C) of a function g ∈ C 0 (R ν , C) ∩ L ∞ (R ν , C) is ∥ g ∥ L ∞ = ess. sup v | g(v) | = sup v | g(v) | . (4.1)
The norm

∥ φ ∥ L q in L q (R ν , C) of the complex-valued functions on R ν is ∥ φ ∥ L q = R ν | φ(v) | q dv 1/q , 1 ≤ q < +∞ . (4.2)
(ii) We define the vector spaces of complex-valued functions, C b , C 0 , H 0 and H 1 , such that

C b = C 0 (R ν , C) ∩ L ∞ (R ν , C) , (4.3)
C 0 = g ∈ C 0 (R ν , C) , | g(v) |→ 0 as ∥ v ∥→ +∞ ⊂ C b ⊂ L ∞ (R ν , C) , (4.4)
H 0 = C 0 ∩ L 2 (R ν , C) , (4.5
)

H 1 = L 1 (R ν , C) ∩ L 2 (R ν , C) . ( 4
(resp. p H targ ) ∈ C 0 (R ν , R) ∩ L 1 (R ν , R) ∩ L 2 (R ν , R) . (4.7) Remark 2. (i) A probability density function on R ν is always in L 1 (R ν , R). The unusual hypothesis is the belonging to C 0 (R ν , R) ∩ L 2 (R ν , R).
As part of the method we propose, the prior probability density p H of H will be estimated using the kernel density estimation (KDE) method of the nonparametric statistics with the training set (see Section 6.1). In this situation, the hypothesis will be verified. This hypothesis will allow us an efficient finite representation of the weak formulation to be constructed.

(ii) The prior pdf p H will effectively be used to construct the posterior pdf p post H by using the KLDMP under the constraint defined by the target set D targ as we have previously explained. The pdf p H targ is not used in the methodology proposed and moreover, if it were to be used, there would be a difficulty because it is assumed that ν is large enough and that N r is not sufficiently large for obtaining a converged estimate of p H targ using the KDE method with the target set D targ . The hypothesis defined by Eq. (4.7) for p H targ will strongly be used and is coherent with the one introduced for p H . Lemma 2 (Fourier transform of the probability measures P H and P H targ ). The Fourier transform of the probability measures P H and P H targ (characteristic functions

) v → f (v) = E{e i ⟨ v, H ⟩ } and v → f targ (v) = E{e i ⟨ v, H targ ⟩ } from R ν into C are such that f ∈ H 0 , f targ ∈ H 0 . (4.8) PROOF. (Lemma 2). Since ∀ v ∈ R ν , f (v) = R ν e i ⟨ v, η ⟩ p H (η)
dη and taking into account Eq. (4.7), it is deduced that f ∈ C 0 because p H ∈ L 1 (R ν , C), and on the other hand, since

p H ∈ L 2 (R ν , C) then f ∈ L 2 (R ν , C).
Consequently, f ∈ H 0 and the proof is similar for f targ .

Remark 3. Since f is the Fourier transform of probability measure P H (dη) = p H (η) dη (positive bounded measure), then it is known (Bochner's theorem) that f (same properties for f targ ) is a nonnegative-type function, that is to say, for all integer m ≥ 1, for all complex numbers z 1 , . . . , z m , and for all vectors v 1 , . . . , v m in R ν , we have

m k=1 m k ′ =1 f (v k -v k ′ ) z k ′ z k ≥ 0 , (4.9)
in which z k is the conjugate of complex number z k . Note that Eq. (4.9) can simply be deduced without evoking the

Bochner theorem because k,k ′ f (v k -v k ′ ) z k ′ z k = k,k ′ E{e i ⟨ v k -v k ′ , H ⟩ } z k ′ z k = E | m k ′ =1 e -i ⟨ v k ′ , H ⟩ } z k ′ | ≥ 0.
We have a similar property to Eq. (4.9) for f targ , that is to say,

m k=1 m k ′ =1 f targ (v k -v k ′ ) z k ′ z k ≥ 0 . (4.10)
Lemma 3 (Convolution operator A). Let f be defined by Eq. (3.5), which belongs to H 0 (see Lemma 2). Let A be the convolution operator defined, for all v ∈ R ν , by

(Aφ)(v) = R ν f (v -v ′ ) φ(v ′ ) dv ′ , (4.11) in which φ is a function from R ν in C. Since f ∈ H 0 , then f ∈ L ∞ (R ν , C) and A is a continuous linear operator from L 1 (R ν , C) into C b and we have ∥ Aφ ∥ L ∞ ≤ c 1 ∥ φ ∥ L 1 in which c 1 is such that ∥ f ∥ L ∞ = c 1 < +∞.
This Lemma holds by replacing f by f targ .

PROOF. (Lemma 3). For the proof of this usual result, we refer the reader, for instance, to [START_REF] Dieudonné | Treatise on Analysis[END_REF] or to Propositions 1 and 2, Pages 164-165 of [START_REF] Soize | Méthodes Mathématiques en Analyse du Signal)[END_REF].

Remark 4. Since f ∈ H 0 , then f also belongs to L 2 (R ν , C) and consequently, A is also a continuous linear operator from L 1 (R ν , C) into L 2 (R ν , C), and also from L 2 (R ν , C) into C 0 , but we do not need to use these properties. This remark holds by replacing f by f targ .

Lemma 4 (Hermitian form F on H 1 associated with f ). Let φ → F(φ) be the functional defined on H 1 with values in C, such that

F(φ) = R ν R ν f (v -v ′ ) φ(v ′ ) φ(v) dv ′ dv . (4.12)
Then F is a positive Hermitian form on H 1 . There is a finite positive constant 0 < c 1 < +∞, such that for all φ in H 1 ,

| F(φ) | ≤ c 1 ∥ φ ∥ 2 L 1 < +∞ , F(φ) ≥ 0 . (4.13)
This Lemma holds by replacing f by f targ .

PROOF. (Lemma 4). (i) Using Eq. (4.11), we have C). The property F(φ) ≥ 0 is similar to the one defined by Eq. (4.9) or (4.10). Since | F(φ) | < +∞, F(φ) exists and

| F(φ) | = | R ν (Aφ)(v) φ(v) dv | ≤ R ν | (Aφ)(v) | | φ(v) | dv. As φ ∈ H 1 , we thus have φ ∈ L 1 (R ν , C). Using Lemma 3, since f ∈ H 0 , ∥ Aφ ∥ L ∞ ≤ c 1 ∥ φ ∥ L 1 . Consequently, | F(φ) | ≤ ∥ Aφ ∥ L ∞ R ν | φ(v) | dv = ∥ Aφ ∥ L ∞ ∥ φ ∥ L 1 , and | F(φ) | ≤ c 1 ∥ φ ∥ 2 L 1 < +∞, which proves the first part of Eq. (4.13). (ii) Let us proof that F(φ) ∈ R. We have F(φ) = R ν R ν f (v -v ′ ) φ(v ′ ) φ(v) dv ′ dv = R ν R ν f (v ′ -v) φ(v) φ(v ′ ) dv dv ′ and f (v ′ -v) = f (v -v ′ ). We then have F(φ) = F(φ). (iii) Since φ ∈ H 1 , we have φ ∈ L 1 (R ν ,
F(φ) = R ν R ν f (v -v ′ ) φ(v ′ ) φ(v) dv ′ dv = E | R ν e -i ⟨ v ′ , H ⟩ φ(v ′ ) dv ′ | 2 ≥ 0.
Theorem 1 (Representation of Hermitian forms F and F targ ). (i) The positive Hermitian form φ → F(φ) : H 1 → R + , defined by Eq. (4.12), can be rewritten as

F(φ) = R ν f (v) ψ(v) dv , (4.14)
in which v → ψ(v) : R ν → C is such that ψ = φ ∨ * φ ∈ L 1 (R ν , C), with φ ∨ (v) = φ(-v) and where (φ ∨ * φ)(v) = R ν φ ∨ (v -v ′ ) φ(v ′ ) dv ′ is the convolution product of φ ∨ with φ. (ii) The Fourier transform η → φ(η) = R ν e -i ⟨ η , v ⟩ φ(v) dv of φ ∈ H 1 is such that φ ∈ H 0 . For all η in R ν , the Fourier transform ψ(η) = R ν e -i ⟨ η , v ⟩ ψ(v) dv of ψ is such that ψ(η) = | φ(η) | 2 ∈ R + and ψ is a positive-valued function that belongs to C 0 ∩ L 1 (R ν , C). (iii) We have the following representation of F(φ), ∀ φ ∈ H 1 , F(φ) = E{ | φ(H) | 2 } = E{ ψ(H) } < +∞ . (4.15)
Results (i) and (ii) hold by replacing f by f targ and F by F targ and we thus have the following representation of F targ (φ),

∀ φ ∈ H 1 , F targ (φ) = E{ | φ([V] T Q targ ) | 2 } = E{ ψ([V] T Q targ )} < +∞ . ( 4
F(φ) = R ν R ν f (u ′ ) φ(v -u ′ ) φ(v) dv du ′ . Using the notation φ ∨ , F(φ) can be rewritten as F(φ) = R ν f (u ′ ) R ν φ ∨ (u ′ - v) φ(v) dv du ′ = R ν f (v) ψ(v) dv in which ψ(v) = R ν φ ∨ (v -v ′ ) φ(v ′ ) dv ′ . As φ ∈ H 1 , we thus have φ and thus φ ∨ in L 1 (R ν , C). The convolution product of two functions in L 1 (R ν , C) is a function in L 1 (R ν , C).
(ii) Function φ belongs to H 1 . Therefore, φ belongs to H 0 . Since ψ ∈ L 1 (R ν , C), its Fourier transform η → ψ(η) on R ν belongs to C 0 and is written as ψ

(η) = φ ∨ (η) × φ(η) = | φ(η) | 2 ,
which shows that ψ is a positive-valued function. Since φ ∈ H 0 , this means that φ ∈ C 0 ∩ L 2 (R ν , C) and then ψ is a positive-valued function that belongs to C 0 ∩ L 1 (R ν , C). (iii) Using Eqs. (3.5), (4.14), and Hypothesis 1 yield

F(φ) = R ν R ν e i ⟨ η , v ⟩ ψ(v) dv p H (η) dη = R ν ψ(η)p H (η) dη = R ν | φ(η) | 2 p H (η) dη. Since φ is a continuous function on R ν , φ(H) is a C-valued random variable such that F(φ) = E{ | φ(H) | 2 },
Remark 5 (Functional F(ψ) on L 1 (R ν , C)). Due to Theorem 1-(i), functional φ → F(φ) on H 1 can also be viewed as a functional ψ → F(ψ) on L 1 (R ν , C), such that F(ψ) = F(φ) and F(ψ) = R ν f (v) ψ(v) dv , (4.17)
with | F(ψ) | < +∞ (due to Eq. (4.15)). Since C 0 ⊂ L ∞ (R ν , C), the right-hand side of Eq. (4.17) can be seen as the duality bracket of

L ∞ (R ν , C) and L 1 (R ν , C). Similarly, φ → F targ (φ) can be viewed as a functional ψ → F targ (ψ) = F targ (φ) on L 1 (R ν , C).
Definition 2 (Weak formulation of the constraint). Lemma 2 and Theorem 1 can be applied to function

v → f (v) = E{ exp(i ⟨ v, H ⟩) } and to function v → f targ (v) = E{ exp(i ⟨v, [V] T Q targ ⟩) } on R ν . For all φ in H 1 , we have F(φ) = R ν f (v) ψ(v) dv ≥ 0 , F targ (φ) = R ν f targ (v) ψ(v) dv ≥ 0 , (4.18)
in which ψ = φ ∨ * φ ∈ L 1 (R ν , C) (see Theorem 1-(i)). We thus define a weak formulation of the constraint { f (v) = f targ (v), ∀ v ∈ R ν }, as follows ∀ φ ∈ H 1 , F(φ) = F targ (φ) . (4.19)
Using Eqs. (4.15) and (4.16) for F and F targ yields, for all φ in H 1 and

ψ = φ ∨ * φ in L 1 (R ν , C), E{ ψ(H) } = B c ( ψ) , (4.20
)

in which ψ(η) = | φ(η) | 2 ∈ C 0 ∩ L 1 (R ν , C) , {η → φ(η) = R ν e -i ⟨ η , v ⟩ φ(v) dv} ∈ H 0 (see Theorem 1-(ii))
, and

B c ( ψ) = E{ ψ([V] T Q targ ) } < +∞ . (4.21)
For N r sufficiently large, the right-hand side of Eq. (4.21) can be estimated with the realizations of Q targ , yielding

B c ( ψ) = 1 N r N r r ′ =1 ψ([V] T qr ′ targ ) < +∞ . (4.22)
Remark 6 (Rewriting the weak formulation). It can easily be seen that the weak formulation defined by Eq. (4. [START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF]) can be written as follows

φ ∈ H 1 , ψ = φ ∨ * φ ∈ L 1 (R ν , C) , F(ψ) = F targ (ψ) , (4.23)
in which using Remark 5 yields

F(ψ) = R ν f (v) ψ(v) dv = E{ ψ(H)} , (4.24)
F targ (ψ) = R ν f targ (v) ψ(v) dv = B c ( ψ) . (4.25)

Construction and analysis of a finite representation of the constraint derived from the weak formulation

We have seen that the weak formulation of constraint was defined for φ ∈ H 1 (see Definition 2). If the Gaussian KDE of p H targ was used with a finite number of realizations ηr targ = [V] T qr targ of H targ , then p H targ would be decreasing as η → exp(-∥ η ∥ 2 /(2 s2 )) in which s would be, for instance, the Silverman bandwidth [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations[END_REF]. Then the Fourier transform would be decreasing as v → exp(-( s2 /2) ∥ v ∥ 2 ). This remark leads us to restrict the weak formulation defined by Eq. (4.19) to a subspace H 1, µ ⊂ H 1 defined as follows.

Definition 3 (Definition of vector space H 1, µ ). Let µ(dv) be the probability measure on R ν defined by

µ(dv) = p ν (v) dv , p ν (v) = s ν ν 2π ν/2 exp - νs 2 2 ∥ v ∥ 2 , (5.1)
in which s is written as

s = 4 N r (2 + ν) 1/(ν+4)
.

(5.

2)

The subspace H 1, µ of H 1 is then defined by

H 1, µ = { φ : R ν → C , φ(v) = ∼ φ(v) p ν (v) , ∼ φ ∈ L 2 µ (R ν , C) } , (5.3)
in which the Hilbert space L 2 µ (R ν , C) is equipped with the inner product and the associated norm,

( ∼ φ , ∼ φ ′ ) L 2 µ = R ν ∼ φ(v) ∼ φ ′ (v) p ν (v) , ∥ ∼ φ ∥ L 2 µ = R ν | ∼ φ(v) | 2 p ν (v) 1/2 . ( 5.4)
Remark 7. (i) The choice of the probability measure defined by Eq. ( 5.1) will appear later. But already now, it can be seen that, for all η in R ν ,

pν (η) = R ν e -i ⟨ η , v ⟩ p ν (v) dv = exp - 1 2νs 2 ∥ η ∥ 2 .
(5.5)

Since s → 0 as N r → +∞, the sequence of measures ( √ 2π s √ ν) -ν pν (η) dη converges to the Dirac measure δ 0 (η) on R ν in the space of the bounded measures on R ν . (ii) Parameter ν has been introduced in the exponential of p ν for numerical conditioning. It can be seen that, if all the components of η are of order 1, then ∥ η ∥ 2 ∼ ν and consequently, ∥ η ∥ 2 /ν ∼ 1. (iii) Note also that s defined by Eq. (5.2) is the Silveman bandwidth corresponding to N r realizations of Q targ and not to the N d realizations of Q. As we explained, we have to construct a finite representation of the constraint, which is consistent with a "projection on the model" of the target set of the realizations, that is to say, of the realizations of the random variable H targ = [V] T Q targ . We recall that the pdf of H targ is assumed to be unknown and will not be estimated with the Gaussian KDE from the training set D targ = {q 1 targ , . . . , q N r targ }. Indeed, we have only assumed that the unknown probability measure P H targ of H targ , which is unknown, admits a density p H targ with respect to dη, which belongs to

C 0 (R ν , C) ∩ H 1 (see Eq. (4.7)). Theorem 2 (Properties of the Fourier transform of H 1, µ ). (i) For all ∼ φ ∈ L 2 µ (R ν , C), the complex-valued function v → φ(v) = ∼ φ(v) p ν (v) on R ν belongs to H 1, µ ⊂ H 1 , in which p ν is defined by Eq. (5.1). Let η → φ(η) = R ν e -i ⟨ η , v ⟩ φ(v) dv be the Fourier transform of φ on R ν . Then φ belongs to H 0 . (ii) The complex-valued function φ is analytic on R ν . (iii) Let H targ = [V] T Q targ be the R ν -valued random variable defined in Lemma 1, whose N r realizations are {η r targ = [V] T qr targ , r = 1, . . . , N r }. If φ(η r targ ) = 0 for all r in N * = N\0, then ∼ φ = 0, dv-almost everywhere. PROOF. (Theorem 2). (i) ∥ φ ∥ L 1 = R ν | ∼ φ(v) | p ν (v) dv = R ν | ∼ φ(v) | p ν (v) 1/2 p ν (v) 1/2 dv ≤ R ν | ∼ φ(v) | 2 p ν (v) dv 1/2 × R ν p ν (v) dv 1/2 = ∥ ∼ φ ∥ L 2 µ < +∞ because ∼ φ ∈ L 2 µ (R ν , C). In addition, it can be seen that ∥ φ ∥ L 2 = R ν | ∼ φ(v) | 2 p ν (v) 2 dv ≤ (sup v p ν (v)) R ν | ∼ φ(v) | 2 p ν (v) dv = s ν ((ν/(2π)) ν/2 ∥ ∼ φ ∥ 2 L 2 µ < +∞. Consequently, φ ∈ L 1 (R ν , C) ∩ L 2 (R ν , C) = H 1 .
Since φ ∈ H 1 , its Fourier transform φ belongs to H 0 .

(ii) We have now to prove that the complex-valued function φ in analytic on R ν using a proof similar to the one of Proposition II.2.36 of [START_REF] Krée | Mathematics of Random Phenomena[END_REF]). For j = 1, . . . , ν, let z j = u j +i η j ∈ C with u j and η j in R. Let ε > 0 be a given real number.

Let u = (u 1 , . . . , u ν ) be in R ν such that max j | u j |≤ ε. Consequently, for all v = (v 1 , . . . , v ν ) in R ν , | e j z j v j |= e j u j v j ≤ e ε | v | , with | v |= |v 1 |+ . . . + |v ν |. Let us consider the Laplace transform φ(z) = R ν e j z j v j ∼ φ(v) p ν (v) dv for z ∈ C ν of function v → φ(v) = ∼ φ(v) p ν (v). We have | φ(z) | ≤ R ν e ε | v | | ∼ φ(v) | p ν (v) dv = R ν | ∼ φ(v) | p ν (v) 1/2 e ε | v | p ν (v) 1/2 dv ≤ R ν | ∼ φ(v) | 2 p ν (v) dv 1/2 R ν e 2ε | v | p ν (v) dv 1/2 < +∞ because we have R ν | ∼ φ(v) | 2 p ν (v) dv = ∥ ∼ φ ∥ 2 L 2 µ < +∞ and R ν e 2ε | v | p ν (v) dv = s ν (ν/(2π)) ν/2 R ν exp(-(νs 2 /2) ∥ v ∥ 2 + 2ε | v |) dv < +∞. Consequently, z → φ(z) exits in the domain D φ = {z ∈ C ν , u j ∈] -ε , ε [ , j = 1, . . . , ν} and is a holomorphic function in D φ ⊂ C ν . Therefore, the conjugate φ(η) of the Fourier transform φ(η) = R ν e -i ⟨ η , v ⟩ φ(v) dv can be written as φ(η) = φ(0 + i η), which shows that η → φ(η) is a C-valued analytic function on R ν .
(iii) Finally, we have to prove the last assertion of Theorem 2. It should be noted that, although φ is an analytic function on R ν , the conditions φ(η r targ ) = 0 for all r in N * do not imply, a priori, that ∼ φ = 0 because the independent realizations {η r targ , r ∈ N * } constitute a countable number of zeros of φ. However, φ(η r targ ) = 0 for r ∈ {1, . . . , N r } implies that (1/N r ) N r r=1 | φ(η r targ ) | 2 = 0 and therefore, for N r → +∞, implies that E{| φ(H targ ) | 2 } = 0. Since φ ∈ H 0 (and is also analytic), φ is a continuous function on R ν and H targ has a probability measure p H targ (η) dη on R ν for which the pdf belongs to C 0 (R ν , C) ∩ H 1 (see Eq. (4.7)). We can then conclude that E{| φ(H targ )

| 2 } = 0 implies φ = 0. Since v → φ(v) = ∼ φ(v) p ν (v) ∈ H 1, µ ⊂ H 1 , then the Plancherel equality, ∥ φ ∥ L 2 = (2π) -ν/2 ∥ φ ∥ L 2
shows that φ = 0 implies φ = 0, dvalmost everywhere, and since p ν (v) > 0 for all v in R ν , then this implies that ∼ φ = 0, dvalmost everywhere.

About a possible use of a polynomial representation

The weak formulation defined by Eq. (4.19), restricted to subspace H 1, µ of H 1 , is written as

∀ φ ∈ H 1, µ , F(φ) = F targ (φ) . (5.6)
For constructing a finite representation of Eq. (5.6), a classical method consists in performing the expansion of φ ∈ H 1, µ with respect to the orthogonal polynomials in L 2 µ (R ν , C). Let α = (α 1 , . . . , α ν ) be the multi-index in N ν . We introduce the classical notations:

|α|= α 1 + . . . + α ν , α! = α 1 ! × . . . × α ν !, i α = i |α| , and for η = (η 1 , . . . , η ν ) ∈ R ν , η α = η α 1 1 ×. . .×η α ν ν . For v = (v 1 , . . . , v ν ) ∈ R ν , let H α (v) = H α 1 (v 1 )×. . .×H α ν (v ν
) be the multi-index Hermite polynomial on R ν of degree |α| such that the real Hermite polynomials H k (y) on R are H 0 (y) = 1, H 1 (y) = y, H 2 (y) = y 2 -1, H 3 (y) = y 3 -3y, etc. It is known that the countable family

{ ∼ ϕ α , α ∈ N ν } such that ∼ ϕ α (v) = (α! ) -1/2 H α (s √ ν v) is a Hilbertian basis of L 2 µ (R ν , R)
and is also a Hilbertian basis of L 2 µ (R ν , C) considered as the complexified space of L 2 µ (R ν , R). We then have (

∼ ϕ α , ∼ ϕ β ) L 2
µ = δ αβ for α and β in N ν in which δ αβ is the Kronecker symbol. Therefore, any function

∼ φ in L 2 µ (R ν , C) can be written as ∼ φ(v) = +∞ α, |α| = 0 ξ α ∼ ϕ α (v) in which ξ α = ( ∼ φ, ∼ ϕ α ) L 2 µ = R ν ∼ φ(v) ∼ ϕ α p ν (v) dv.
The series in the right-hand side member of the expansion of

∼ φ is convergent in L 2 µ (R ν , C) and we have ∥ ∼ φ ∥ 2 L 2 µ = +∞ α, |α| = 0 | ξ α | 2 < +∞. It can be deduced that all function φ in H 1, µ ⊂ H 1 can be written as φ(v) = +∞ α, |α| = 0 ξ α ∼ ϕ α (v) p ν (v)
. For all v ∈ R ν , from a classical formula, we can deduce the following one,

R ν e -i ⟨ η ,v ⟩ H α (s √ ν v) p ν (v) dv = -i s √ ν | α | η α exp - 1 2νs 2 ∥ η ∥ 2 .
Thus, the Fourier transform φ of φ ∈ H 1, µ belongs to H 0 and can be written, for all η in R ν , as φ(η)

= +∞ α, |α| = 0 ξ α ϕ α (η) in which ϕ α (η) = 1 √ α! -i s √ ν | α | η α exp - 1 2νs 2 ∥ η ∥ 2 , α ∈ N ν .
(5.7)

Using the family { ϕ α } α defined by Eq. (5.7) and Definition 2, the finite representation of the weak formulation defined by Eq. (5.6), can be written as,

E{h c α (H)} = b c α , α ∈ {α (1) , . . . , α (N r ) } , (5.8)
in which

h c α (η)} = | ϕ α (η) | 2 = Ψ α (η) , α ∈ N ν , (5.9)
Ψ α (η) = (νs 2 ) -| α | α! η 2α exp - 1 νs 2 ∥ η ∥ 2 , (5.10)
and where, for N r sufficiently large,

b c α ≃ 1 N r N r r ′ =1 Ψ α ([V] T qr targ ) , α ∈ N ν .
(5.11)

The finite representation defined by Eq. (5.8) with Eqs. (5.9) and (5.10) is not efficient as soon as ν is large. Indeed, since α = (α 1 , . . . , α ν ) is a multi-index in N ν and since N r multi-indices α (1) , . . . , α (N r ) have to be considered, even for N r not too large, the number of components of b c can become enormous and unrealistic as soon as ν is large (a few tens or a few hundreds). In addition, Ψ α defined by Eq. (5.10) does not depend on the sampling defined by the points of the target set D targ .

Construction of an adapted finite representation of the functional constraint

The following Lemma gives the construction of a family in H 1, µ , which is based on the sampling points of the target set D targ .

Lemma 5 (Construction of a family

{ ψ r } r in C 0 ∩ L 1 (R ν , C)). Let r be fixed in {1, . . . , N r }. (i) Let v → φ r (v) ∈ H 1, µ defined, for all v in R ν , by φ r (v) = ∼ φ r (v) p ν (v) , ∼ φ r (v) = exp( i ⟨v, [V] T qr targ ⟩) , (5.12)
in which p ν is defined by Eq. (5.1) with Eq. (5.2) and where

∼ φ r ∈ L 2 µ (R ν , C) such that ∥ ∼ φ r ∥ L 2 µ = 1.
Then the Fourier transform η → φr (η) = R ν e -i ⟨ η , v ⟩ φ r (v) dv belongs to H 0 and is written as

φr (η) = exp - 1 2νs 2 ∥ η -[V] T qr targ ∥ 2 , ∀ η ∈ R ν .
(5.13)

Eq. (5.13) shows that φr , which is analytic (see Theorem 2), also belongs to H 0 ∩ L q (R ν , C) for all 3 ≤ q < +∞.

(ii) Let ψ r ∈ L 1 (R ν , C) be the function defined by ψ r = φ ∨ r * φ r . Its Fourier transform is such that (see Theorem 1) ψr ∈ C 0 ∩ L 1 (R ν , C) and is written as

ψr (η) = | φr (η) | 2 = exp - 1 νs 2 ∥ η -[V] T qr targ ∥ 2 , ∀ η ∈ R ν . (5.14)
Note that ψr is also in H 1 ∩ L q (R ν , C) for all 3 ≤ q < +∞ and function ψ r is written as

ψ r (v) = s 2 ν 2π ν/2 p ν (v) 1/2 exp i ⟨ v, [V] T qr targ ⟩ , ∀v ∈ R ν . (5.15) PROOF. (Lemma 5). (i) We have ∥ ∼ φ r ∥ L 2 µ = R ν p ν (v) dv = 1.
From Eqs. (5.1) and (5.12), and introducing

σ = 1/(s √ ν) yield φr (η) = R ν e -i ⟨ η , v ⟩ exp( i ⟨ v, [V] T qr targ ⟩) s ν 2π ν exp -νs 2 2 ∥ v ∥ 2 dv = R ν exp( i ⟨ v, [V] T qr targ - η ⟩) ( √ 2π σ) -ν exp -(2σ 2 ) -1 ∥ v ∥ 2 dv
, which gives Eq. (5.13). (ii) Eq. (5.14) is obtained by substituting Eq. (5.13)

in ψr (η) = | φr (η) | 2 . For all v ∈ R ν , we have ψ r (v) = (2π) -ν R ν e i ⟨ v, η ⟩ ψr (η) dη = (2π) -ν R ν e i ⟨ v, η ⟩ exp -1 νs 2 ∥ η - [V] T qr targ ∥ 2 dv = (s/2) ν (ν/π) ν/2 exp i ⟨ v, [V] T qr targ ⟩ -νs 2 4
∥ v∥ 2 , which can be rewritten as Eq. (5.15) by using Eq. (5.1).

Lemma 6 (Orthonormal family {

∼ φ r , r = 1, . . . , N r } in L 2 µ (R ν , C) for N r → +∞ ). For r ∈ {1, . . . , N r }, let ∼ φ r be the function in L 2
µ (R ν , C) defined by Eq. (5.12). For N r → +∞, the family

{ ∼ φ r , r = 1, . . . , N r } in L 2 µ (R ν , C) goes to an orthonormal family in L 2 µ (R ν , C): lim N r →+∞ ( ∼ φ r , ∼ φ r ′ ) L 2 µ = δ rr ′ . Let L 2,(N r ) µ = span { ∼ φ 1 , . . . , ∼ φ N r
} be the subspace of

L 2 µ (R ν , C) spanned by { ∼ φ r , r = 1, . . . , N r }.
For N r → +∞, the sequence of subspaces L 2,(N r ) µ goes to a subspace that is dense in L 2 µ (R ν , C).

PROOF. (Lemma 6). We have (

∼ φ r , ∼ φ r ′ ) L 2 µ = R ν ∼ φ r (v) ∼ φ r ′ (v) p ν (v) dv = R ν exp( i ⟨ v, [V] T ( qr targ -qr ′ targ ⟩) p ν (v) dv.
Using Eq. (5.5) allows for writing (

∼ φ r , ∼ φ r ′ ) L 2 µ = exp(-(2νs 2 ) -1 ∥ [V] T ( qr targ -qr ′ targ ∥ 2 ). For r = r ′ , we have ∥ ∼ φ r ∥ L 2 µ = 1.
Let us now consider the case r ̸ = r ′ . For N r → +∞, we have s → 0, and consequently, (

∼ φ r , ∼ φ r ′ ) L 2 µ → 0.
We then have proven the first part of the Lemma. We have now to prove that, for any

∼ φ in L 2 µ (R ν , C) and for N r → +∞ if ( ∼ φ, ∼ φ r ) L 2 µ = 0 , ∀ r, then ∼ φ = 0. We have ( ∼ φ, ∼ φ r ) L 2 µ = R ν ∼ φ(v) ∼ φ r (v) p ν (v) dv = R ν exp(-i ⟨ v, [V] T qr targ ⟩) φ(v) dv, in which φ(v) = ∼ φ(v) p ν (v). We then obtain ( ∼ φ, ∼ φ r ) L 2 µ = φ([V] T qr targ ). Using Theorem 2, the condition φ([V] T qr targ ) = 0, ∀ r ∈ N * implies that
c (H)} = b c on R N r , (5.16
)

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

h c r (η) = exp - 1 νs 2 ∥ η -[V] T qr targ ∥ 2 , ∀ η ∈ R ν , (5.17)
b c r = E exp - 1 νs 2 ∥ [V] T ( Q targ -qr targ ∥ 2 , (5.18)
which can be estimated, for N r sufficiently large, by

b c r = 1 N r N r r ′ =1 exp - 1 νs 2 ∥ [V] T ( qr ′ targ -qr targ ∥ 2 .
(5.19)

Remark 8. (i) Definition 4 shows that if the random variables H and H targ = [V] T Q targ are isonomic, then the constraint defined by Eq. (5.16) is exactly satisfied. Consequently, the use of the KLDMP for estimating the posterior probability measure under this constraint will be well posed.

(ii) Let ψ(N r ) be the function in

C 0 ∩ L 1 (R ν , C), such that for all η ∈ R ν , ψ(N r ) (η) = 1 N r N r r=1 (s √ νπ) -ν ψr (η). Using
Eq. (5.14), we have the following equality in the space of the bounded measures, ψ(N r ) (η

) dη = 1 N r N r r=1 ( √ 2π ∼ σ) -ν exp -(2 ∼ σ 2 ) -1 ∥ η -[V] T qr targ ∥ 2 dη, in which ∼ σ = s √ ν/2.
For N r → +∞, since ∼ σ → 0 because s → 0, it can be seen that the right-hand side of this last equality goes to the probability measure p H targ (η) dη in which

{ η r targ = [V] T qr targ , r = 1, . . . , N r } are N r independent realizations of H targ = [V] T Q targ .
We then have the following convergence property in the space of the bounded measures, lim N r →+∞ ψ(N r ) (η) dη = p H targ (η) dη. This result contributes to justify the construction presented in Definition 4.

The following Lemma will be used in the next section for analyzing the existence and uniqueness of the posterior probability measure constructed by using the KLDMP.

Lemma 7 (Positive definiteness of matrix E{h c (H) ⊗ h c (H)}). Let us consider any finite positive fixed value of integer N r . Let H be the R ν -valued random variable whose probability measure P H (dη) is such that the following

M +0 N r -valued matrix exits, E{h c (H) ⊗ h c (H)} = R ν h c (η) ⊗ h c (η) P H (dη) , (5.20)
in which, h c r (η) = exp -(νs 2 ) -1 ∥ ηη r targ ∥ 2 with η r targ = [V] T qr targ for r = 1, . . . , N r . Therefore, this matrix is positive definite, E{h c (H) ⊗ h c (H)} ∈ M + N r .

(5.21)

PROOF. (Lemma 7). For all w = (w 1 , . . . ,

w N r ) ∈ R N r , we have ⟨ E{h c (H) ⊗ h c (H)} w , w ⟩ = E{⟨ h c (H) , w ⟩ 2 } ≥ 0.
We then have to prove that, for all w in R N r with ∥ w ∥ ̸ = 0, we have E{⟨ h c (H) , w ⟩ 2 } > 0 or equivalently, that

E{⟨ h c (H) , w ⟩ 2 } = 0 ⇔ w = 0. We have E{⟨ h c (H) , w ⟩ 2 } = R ν ⟨ h c (η) , w ⟩ 2 P H (dη) = 0 ⇐⇒ R ν ⟨ h c (η) , w ⟩ 2 dη = 0 ⇐⇒ N r r=1 N r r ′ =1 w r w r ′ R ν exp ( -(νs 2 ) -1 γ rr ′ (η)) dη = 0 in which for all η in R ν , γ rr ′ (η) = ∥ η -η r targ ∥ 2 +∥ η - η r ′ targ ∥ 2 = 2 ∥ η -1 2 (η r targ + η r ′ targ )∥ 2 + 1 2 ∥ η r targ -η r ′ targ ∥ 2 . Consequently, E{⟨ h c (H) , w ⟩ 2 } = 0 ⇐⇒ N r r=1 N r r ′ =1 w r w r ′ exp ( - (2νs 2 ) -1 ∥ η r targ -η r ′ targ ∥ 2 ) R ν exp ( -2(νs 2 ) -1 ∥ η -1 2 (η r targ + η r ′ targ )∥ 2 ) dη = 0 ⇐⇒ ((s/2) √ 2πν) ν N r r=1 N r r ′ =1 w r w r ′ exp ( - (2νs 2 ) -1 ∥ η r targ -η r ′ targ ∥ 2 ) = 0 ⇐⇒ w = (w 1 , . . . , w N r ) = 0 N r .

Illustration of the numerical behavior of the functional constraint

In order to illustrate the numerical behavior of the finite representation of the weak formulation of the functional constraint, we consider the following simple numerical case. We assume that H and H targ = [V] T Q targ , which are statistically independent, are Gaussian random vectors. Random vector H is centered and with an identity covariance matrix [I ν ] (see Eq. (2.6)). The mean value of random vector H targ is written as η targ = m targ a ∈ R ν in which m targ is given in R and where a ∈ R ν is any given realization of a uniform random vector on [0 , 1] ν with independent components. The covariance matrix of H targ is written as σ targ [I ν ] in which σ targ is given in ]0 , +∞[. For analyzing the numerical behavior of Eqs.

Kullback-Leibler minimum principle for estimating the posterior model

In this section we reuse part of the developments that we presented in paper [START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]. We do not want to limit ourselves to referring the reader to this reference, because the hypotheses are not the same, the Lemmas and Theorems must be reformulated, and their proofs must be adapted and modified. In addition, the presentation chosen makes it easier to read and understand, thus avoiding going back and forth with this reference.

Prior probability measure of H

Let P H (dη) = p H (η) dη be the prior probability measure on R ν of H, whose probability density function η → p H (η) : R ν → R + is estimated by using the Gaussian kernel-density estimation (KDE) with the training set D d = {η 1 d , . . . , η N 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 SB is the Silverman bandwidth, ν+4) , and where η → ζ(η) : R ν → R + is written as

p H (η) = c ν ζ(η) , ∀η ∈ R ν , c ν = ( √ 2π ŝ) -ν , (6.1) in which ŝ = s SB s 2 SB + (N d -1)/N d -1/2 with s SB = (4/(N d (2 + ν))) 1/(
ζ(η) = 1 N d N d j=1 exp - 1 2 ŝ2 ∥ ŝ s SB η j d -η ∥ 2 . (6.2)
We define the potential function η → ϕ(η) : R ν → R, related to p H , which will be used in Lemma 8 and such that

ζ(η) = exp{-ϕ(η)} . (6.3)
With such a modification and using Eq. (2.6), the normalization of H is preserved for any value of N d , that is to say,

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

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

Posterior estimate using the Kullback-Leibler divergence minimum principle

The posterior probability density function η → p post H (η) on R ν of the R ν -valued random variable H post = (H post,1 , . . . , H post,ν ) is estimated by using the KLDMP [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][START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]. This estimation of p post H is performed by using the prior pdf η → p H (η) on R ν in which p H is defined by Eqs. (6.1) and (6.2), and by using the constraint defined by Eq. (5.16). The pdf p post H on R ν , which satisfies the constraint defined by Eq. (5.16) and which is closest to p H defined by Eq. (6.1), is thus the solution of the following optimization problem,

p post H = arg min p∈C ad,p R ν p(η) log p(η) p H (η) dη , (6.6)
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 . (6.7)

Methodology for solving the optimization problem

The constraints defined in admissible set C ad,p are taken into account by introducing two Lagrange multipliers, λ 0 -1 with λ 0 ∈ R + associated with the normalization condition and λ ∈ C ad,λ ⊂ R N r associated with the functional constraint. The admissible set C ad,λ of λ is, a priori, a subset of R N r , which will be defined in Section 6.4 (in fact, we will see that C ad,λ = R N r ). The Lagrange multiplier λ 0 is eliminated as a function of λ. In Eq. (6.6), the posterior pdf p post H is constructed as the limit of a sequence {p H λ } λ of probability density functions of a sequence {H λ } λ of R ν -valued random variables H λ = (H λ,1 , . . . , H λ,ν ) that depend on λ. For λ fixed in C ad,λ , a MCMC algorithm is used for generating the 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 post = {η 1 post , . . . , η N post } is generated. This set is made up of N independent realizations {η ℓ post , ℓ = 1, . . . , N} of H post whose probability measure is p post H (η) dη. The MCMC generator will be 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 removing the transient part to rapidly reach the stationary response associated with the invariant measure for which measure p post H (η) dη is the marginal measure. The ISDE is solved by using the Störmer-Verlet algorithm, which yields an efficient and accurate MCMC algorithm. This algorithm can then easily be parallelized for strongly decreasing the elapsed time on a multicore computer (See Algorithm 1 in Section 6.7-(ii)). Note that this MCMC generator can be considered to belong to the class of Hamiltonian Monte Carlo methods [START_REF] Neal | MCMC using hamiltonian dynamics[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF] but is not similar due to the dissipative term, and is a MCMC algorithm [START_REF] Kaipio | [END_REF][START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization: Estimation[END_REF].

Let us assumed that the optimization problem defined by Eq. ( 6.6) has one solution p post H and that p = p post H 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 ⟩ .
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 ν , (6.8)
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 }). Since Lemma 7 holds for any probability measure P H on R ν with support R ν , we can conclude that the N r constraints defined by the components of Eq. ( 5.16) are algebraically independent. Consequently, 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 post H for λ = λ sol and λ 0 =log(c 0 (λ sol )/c ν ). Consequently, p post H = p H λ sol (. ; λ sol ) and Eq. (6.8) yield

p post H (η) = c 0 (λ sol ) ζ(η) exp{-⟨λ sol , h c (η)⟩} , ∀ η ∈ R ν . (6.9)
Therefore, p post H is the unique solution of the optimization problem defined by Eq. (6.6), in which λ sol will be the unique solution of a convex optimization problem that will be defined by Theorem 3 in Section 6.4) and which will be the solution of the following nonlinear algebraic equation in λ, R ν h c (η) p H λ (η ; λ) dη = b c .

Analysis of the optimization problem

In this section, we study the admissible set of the Lagrange multiplier, we analyze the integrability properties of the probability density function p H λ of H λ , and we give an explicit construction of p H λ . Lemma 8 (Admissible set C ad,λ of Lagrange's multiplier and integrability properties). Let N r be fixed. Let h c (η) = (h c 1 (η), . . . , h c N r (η)) be the function on R ν with values in R N r , defined by Eq. (5.17), and let p H be the prior probability density function on R ν of H, defined by Eq. (6.1). (a) The admissible set C ad,λ of the Lagrange multiplier λ , which is defined by

C ad,λ = λ ∈ R N r | 0 < E{ exp ( -⟨ λ , h c (H) ⟩ ) } < +∞ , (6.10
)

is such that C ad,λ = R N r . (b) For all λ in C ad,λ , let η → V λ (η) be the R-valued function on R ν such that V λ (η) = ϕ(η) + ⟨λ , h c (η)⟩ , (6.11)
in which ϕ(η) =log ζ(η) (see Eq. (6.3)). We then have,

0 < R ν exp{-V λ (η)} dη < +∞ . (6.12)
(c) The pdf η → p H λ (η ; λ), defined by Eq. (6.8), which can be written as

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

0 < c 0 (λ) < +∞ , ∀λ ∈ C ad,λ . (6.14)
(d) For all λ in C ad,λ , we have V λ (η) → +∞ if ∥ η ∥→ +∞, and

R ν ∥ h c (η) ∥ 2 exp{-V λ (η)} dη < +∞ , R ν ∥ [∇ η h c (η)] ∥ F exp{-V λ (η)} dη < +∞ . (6.15)
PROOF. (Lemma 8).

(a) For all η ∈ R ν and for all r ∈ {1, . . . , N r } , Eq. (5.17) shows that 0 < h c r (η) ≤ 1. It can then be deduced that, for all λ ∈ R N r , we have 0 < E{ exp ( -⟨ λ , h c (H) ⟩ ) } < +∞, which proves that C ad,λ = R N r . (b) Using Eqs. (6.1), (6.3), and (6.11), yields

R ν exp{-V λ (η)} dη = c -1 ν R ν exp{-⟨λ , h c (η)⟩} p H (η) dη = c -1
ν E{ exp{-⟨λ , h c (H)⟩} }, which is positive and finite due to Eq. (6.10) and to 0 < c ν < +∞. We have thus proven Eq. (6.12). (c) Using Eqs. (6.12) and (6.13), and since we need to have R ν p H λ (η) dη = 1, we deduce Eq. (6.14). (d) As h c is continuous on R ν , (see Eq. (5.17)), ∀λ ∈ C ad,λ , η → exp{-V λ (η)} is continuous on R ν and then is locally integrable on R ν . Eq. (6.12) implies the integrability at infinity of η → exp{-V λ (η)}. Since η → V λ (η) is continuous on R ν , it can be deduced that V λ (η) → +∞ if ∥ η ∥→ +∞. Using Eq. (5.17), Eq. (6.15) can easily be proven.

Theorem 3 (Construction of the probability measure of H λ ). For all λ in C ad,λ , let

p H λ (η ; λ) = c 0 (λ) ζ(η) exp ( -⟨ λ , h c (η) ⟩ ) (6.16)
be the pdf of H λ (see Eq. (6.8)) with c 0 (λ) satisfying Eq. (6.14)).

(a) The R N r -valued random variable h c (H λ ) is a second-order random variable,

E{∥ h c (H λ ) ∥ 2 } < +∞ . (6.17)
(b) Let λ → Γ(λ) be the real-valued function defined on C ad,λ such that

Γ(λ) = ⟨ λ , b c ⟩ -log c 0 (λ) , (6.18)
in which b c is given in R N r . For all λ in C ad,λ , we have

∇ λ Γ(λ) = b c -E{h c (H λ )} ∈ R N r , (6.19)
[Γ ′′ (λ)] = [cov{h c (H λ )}] ∈ M + N r , (6.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,λ Γ(λ) , (6.21)
which is the unique solution in λ of the following equation,

∇ λ Γ(λ) = 0 N r . (6.22)
The pdf p post H of H post , which satisfies the constraint E{h c (H post )} = b c , is written (see Eq. (6.16)) as

p post H (η) = p H λ sol (η ; λ sol ) , ∀η ∈ R ν . (6.23)
PROOF. (Theorem 3). (a) Using Eq. (6.13), Eq. (6.14), and the first equation Eq. (6.15) yield

E{∥ h c (H λ ) ∥ 2 } = R ν ∥ h c (η) ∥ 2 c 0 (λ) exp{-V λ (η)} dη < +∞ .
(b) Let us prove Eqs. (6. [START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF]) and (6.20) using a similar proof to the one introduced in the discrete case for finding the maximum entropy probability measure [START_REF] Agmon | An algorithm for finding the distribution of maximal entropy[END_REF][START_REF] Kapur | Entropy Optimization Principles with Applications[END_REF]. Eq. (6.11) yields ∇ λ V λ (η) = h c (η) and from Eq. ( 6.13), it can be deduced that

∇ λ p H λ (η ; λ) = c 0 (λ) -1 ∇ λ c 0 (λ) -h c (η) p H λ (η ; λ) . (6.24)
By integrating the two members of Eq. ( 6.24) with respect to η on R ν , we obtain

c 0 (λ) -1 ∇ λ c 0 (λ) = R ν h c (η) p H λ (η ; λ) dη = E{h c (H λ )} . (6.25)
Eq. ( 6.18) yields ∇ λ Γ(λ) = b cc 0 (λ) -1 ∇ λ c 0 (λ), which proves Eq. (6.19) by using Eq. (6.25). Note that Eq. (6.17) implies the existence of the mean value E{h c (H λ )}. Taking the derivative of Eq. ( 6. [START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF]) with respect to λ yields

[Γ ′′ (λ)] = - R ν h c (η) ⊗ ∇ λ p H λ (η ; λ) dη . (6.26)
Substituting Eq. ( 6.25) into Eq. ( 6.24) yields

∇ λ p H λ (η ; λ) = (E{h c (H λ )} -h c (η) ) p H λ (η ; λ), which with Eq. (6.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 r -valued random variable h c (H λ ). Again Eq. (6.17) proves the existence of matrix [Γ ′′ (λ)] as a covariance matrix, which is semi-positive definite. Using Lemma 7, this matrix is positive definite. (c) Since [Γ ′′ (λ)] is a positive-definite matrix for all λ in C ad,λ , it can then be deduced that λ → Γ(λ) is strictly convex on C ad,λ . Therefore, Eq. (6.21) holds, λ sol is unique, and Eq. ( 6. [START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF] shows that E{h c (H λ sol )} = b c . Taking into account Eq. (6.9), the solution is given by Eq. (6.23) and is unique due to the uniqueness of solution λ sol of ∇ λ Γ(λ) = 0 N r .

Dissipative stochastic Hamiltonian system as a MCMC generator of H λ

For the reasons given in Section 6.3, the chosen MCMC generator is 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].

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 ν almost surely (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 (t

-t ′ ) [I ν ].
Theorem 4 (MCMC generator of H λ ). Let h c = (h c 1 , . . . , h c N r ) be the function whose component h c r is defined by Eq. (5.17). Let λ be fixed in C ad,λ . Consequently, Eq. (6.15) of Lemma 8 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

(u 0 , v 0 ) given in R ν × R ν , dU λ (t) = V λ (t) dt , (6.27)
dV λ (t) = L λ (U λ (t)) dt - 1 2 f 0 V λ (t) dt + f 0 dW wien (t) , (6.28)
U λ (0) = u 0 , V λ (0) = v 0 a.s. (6.29)
(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 6.7-(i)) 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), which can be written as

L λ (u) = 1 ζ(u) ∇ u ζ(u) -[∇ u h c (u)] λ . (6.30)
(d) The stochastic solution {(U λ (t), V λ (t)), t ≥ 0} of the ISDE defined by Eqs. (6.27) to (6.29) is unique, has almostsurely continuous trajectories, and is a second-order diffusion stochastic process. For t → +∞, this diffusion process converges to a stationary second-order diffusion stochastic process {(U st λ (τ), 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) , (6.31)
in which p H λ (η ; λ) is the pdf defined by Eq. (6.13).

(e) For t s sufficiently large, H λ is chosen 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 density function is p H λ (η ; λ), consists in solving Eqs. (6.27) to (6.29) for t ∈ [0 , t s ] and then using the realizations of U λ (t s ) (see the numerical aspects in Section 6.7).

PROOF. (Theorem 4). For r ∈ {1, . . . , N r }, function η → h c r (η) defined be Eq. (5.17) is twice continuously differentiable. Since ϕ(u) =log ζ(u) with ζ(u) given by Eq. ( 6.2), it can be deduced that function u → V λ (u) defined by Eq. (6.11) is also twice continuously differentiable. Consequently, u → ∥ ∇ u V λ (u) ∥ is locally bounded on R ν . Using Eqs. (6.11) and (6.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. (6.2), (6.3), and (6.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 , (6.32)
because

∥ [∇ u h c (u)] λ ∥ ≤ ∥ [∇ u h c (u)] ∥ ∥ λ ∥ and ∥ [∇ u h c (u)] ∥ ≤ ∥ [∇ u h c (u)] ∥ F .
From Eqs. (6.16) and (6.2), the first term in the right-hand side member of Eq. (6.32) is finite, while from the second equation (6.15), the second term is also finite. It can then be deduced that the left-hand side member of Eq. (6.32) 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. (6.27) to (6.29) 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. (6.31). For any τ > 0, U st λ (τ) = lim t→+∞ U λ (t + τ) in probability measure.

Iterative algorithm for calculating λ sol

Let us consider Theorem 3. For λ fixed in C ad,λ , the value of Γ(λ) cannot be evaluated in high dimension using Eq. (6.18) due to the presence of the normalization constant c 0 (λ). Consequently, λ sol cannot directly be estimated using, for instance, the gradient descent algorithm applied to the convex optimization problem defined by Eq. (6.21). We will then calculate λ sol as the unique solution in λ of equation ∇ λ Γ(λ) = 0 N r (see Eq. (6.22)), that is to say (see Eq. (6.19)), solving the following equation in λ on R N r , E{h c (H λ )} -b c = 0 N r . (6.33) This equation is solved by 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 -α relax (i) [Γ ′′ (λ i )] -1 ∇ λ Γ(λ i ) , i = 0, 1, . . . , i max , (6.34)
with λ 0 = 0 N r , in which α relax (i) ∈]0 , 1] is a relaxation factor, where ∇ λ Γ(λ) and [Γ ′′ (λ)] are defined by Eqs. (6. [START_REF] Spantini | Goal-oriented optimal approximations of Bayesian linear inverse problems[END_REF]) and (6.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) , (6.35)
in which the error function i → err(i) from {1, . . . , i max } into R + is defined by

err(i) = 1 ∥ b c ∥ ∥ b c -E{h c (H λ i )}∥ . (6.36)
Although Lemma 7 proves that matrix [Γ ′′ (λ i )] is invertible, when N r is large, this matrix can be ill conditioned. In such a case the spectral decomposition of this positive-definite matrix can be introduced as [Γ ′′ (λ

i )] = [φ Γ ] [λ Γ ] [φ Γ ] T with λ Γ,1 ≥ . . . ≥ λ Γ,m ≥ . . . ≥ λ Γ,N r > 0.
Integer m is chosen such that, for all r = 1, . . . , m, λ Γ,r /λ Γ,1 > ε Γ in which ε Γ ≪ 1 (for instance, ε Γ = 10 -6 ). Introducing the diagonal matrix [λ -1 Γ ] with [λ -1 Γ ] rr = 1/λ Γ,r for r = 1, . . . , m and = 0 for r = m + 1, . . . , N r , the following approximation can be used for the inverse, [Γ ′′ (λ

i )] -1 ≃ [φ Γ ] [λ -1 Γ ] [φ Γ ] T .

Numerical implementation

A time-discretization scheme (see for instance [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 to solve the ISDE defined by Eqs. (6.27) to (6.29) for t ∈ [0 , t s ] with the initial condition at t = 0 defined in Theorem 4, in order to generate the constrained learned set 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. 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. Presently we propose to use 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), which is an 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 ave 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][START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning on manifolds (plom) with partition[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF].

(i) 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 the Lagrange multipliers that are computed using Eq. (6.34) with λ 0 = 0 N r . For m = 0, 1, . . . , M s (with M s > 1 an integer), let t m = m ∆t be the time sampling, which is such that t M s = t s with t s = M s ∆t. 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 Theorem 4, let u 1 0 , . . . , u N 0 in R ν 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 Theorem 4. For ℓ = 1, . . . , N, the realizations ∆W ℓ m+1 , u ℓ 0 , and v ℓ 0 are independent of {λ i } i . For i ∈ {0, 1, . . . , i max } and for ℓ ∈ {1, . . . , N}, we introduce the realizations U i,ℓ m = U λ i (t m ; θ ℓ ) and V i,ℓ m = V λ i (t m ; θ ℓ ). For m ∈ {0, 1, . . . , M s -1}, the Störmer-Verlet scheme applied to realization θ ℓ of Eqs. (6.27) to (6.29) yields the following recurrence,

U i,ℓ m+1/2 = U i,ℓ m + ∆t 2 V i,ℓ m , (6.37)
V i,ℓ m+1 = 1 -γ 1 + γ V i,ℓ m + ∆t 1 + γ L λ i-1 (U i,ℓ m+1/2 ) + f 0 1 + γ ∆W ℓ m+1 , (6.38)
U i,ℓ m+1 = U i,ℓ m+1/2 + ∆t 2 V i,ℓ m+1 , (6.39)
with the initial condition

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

L λ i-1 (u) = 1 ζ(u) ∇ u ζ(u) -[∇ u h c (u)] λ i-1 (6.41)
Using Eq. (5.17), for α ∈ {1, . . . , ν} and r ∈ {1, . . . , N r }, the entry [∇ u h c (u)] αr of matrix [∇ u h c (u)] ∈ M ν,N r is written as

[∇ u h c (u)] αr = 2 νs 2 ({[V] T qr targ } α -u α ) exp - 1 νs 2 ∥ [V] T qr targ -u ∥ 2 . (6.42)
The recurrence defined by Eqs. (6.37) to (6.42) allows D H λ i to be calculated as

D H λ i = {η 1 λ i , . . . , η N λ i } , η ℓ λ i = U λ i (t s ; θ ℓ ) = U i,ℓ M s . (6.43)
(ii) Summary of the algorithm. The algorithm for calculating λ sol and the D H post = {η 1 post , . . . , η N post } with η ℓ post = η ℓ λ sol for ℓ = 1, . . . , N is summarized in Algorithm 1. We then obtained the N independent realizations η 1 post , . . . , η N post of the posterior R ν -valued random variable H post . Then, the realizations q 1 post , . . . , q N post of the posterior observations Q post and the realizations w 1 post , . . . , w N post of the posterior control parameter W post are computed using Eq. (2.4), that is to say, err(i) using Eq. (6.36)

Q post = q + [Φ q ] [κ] 1/2 H post , W post = w + [Φ w ] [κ]
9:

λ i+1 = λ i -α relax (i) [Γ ′′ (λ i )] -1 ∇ λ Γ(λ i )

Numerical illustration

We consider a supervised case. The training set D d = {x 1 , . . . , x N d } with x j = (q j d , w j d ) ∈ R n x = R n q ×R n w is made up of N d independent realizations of random variable X = (Q, W), which are generated as explained in Appendix A for which n x = 430 098, n q = 10 098, n w = 420 000, and N d ∈ {100, 200, 300, 400}. The target set D targ is made up of N r independent realizations q r targ ∈ R n q of random variable Q targ for which N r ∈ [50 , N targ ] with N targ ∈ {100, 200, 300, 400}. The R n q -valued random variable Q corresponds to the finite element discretization of a R 3valued random field {U(ω) = (U 1 (ω), U 2 (ω), U 3 (ω) , ω ∈ Ω} and the R n w -valued random variable W is a nonlinear transformation of the finite element discretization of a M + 6 -valued random field {[G(ω)] , ω ∈ Ω}, in which Ω is the open bounded set of R 3 defined in Appendix A.1, and where Q and W are constructed in Appendix A.2. Regarding the presentation of the results, and having to limit the number of figures, the probability density functions and the convergence curves when they will be relative to Q, will be limited to 3 components, denoted by Q obs,1 , Q obs,2 , and Q obs,3 that are also defined at the end of Appendix A.2.

Training set

The training set is generated as explained in Appendix A.2 with the stochastic boundary value problem defined in Appendix A.1. For illustration, Fig. 2 shows one realization θ ∈ Θ of the components (1, 1), (1, 2), and (4, 4) of the M + 6 -valued random field (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 )] in the plane ω 3 = 0.095774 and Fig. 3 shows the corresponding realization of the components k = 1, 2, 3 of the real-valued random field (ω 1 , ω 2 ) → U k (ω 1 , ω 2 , ω 3 ) in the plane In all the calculations and for any value of λ, the number N of the independent realizations {η 1 λ , . . . , η N λ } of the constrained learned set D H λ (see Section 6.7) is fixed to the value N = 1 000 (this value has been obtained from a convergence analysis with respect to N). The parameters of the Störmer-Verlet scheme are M s = 30, ∆t = 0.2188, and f 0 = 4. Due to the choice f 0 = 4, the stationary regime of the ISDE is obtained from instant 21 × ∆t and the realizations is extracted at 30 × ∆t. to N r = 100, display the probability density functions of the random variables Q obs,1 , Q obs,2 , and Q obs,3 , estimated with the training set, with the constrained learned set (the posterior), and the reference. Figs. 7 and8 show that the posteriors are close to the targets (this good result holds for all the components of Q). For this supervised framework, it can be seen that the proposed method performs very well for the quantity of interest Q for which a target has been given, which means that there would also be very good behavior of the method if used in an unsupervised setting. a) pdf q → p Q obs,1 (q) for N r = 20 (b) pdf q → p Q obs,2 (q) for N r = 20 (c) pdf q → p Q obs,3 (q) for N r = 20 (d) pdf q → p Q obs,1 (q) for N r = 100 (e) pdf q → p Q obs,2 (q) for N r = 100 (f) pdf q → p Q obs,3 (q) for N r = 100 a target. In the case of the numerical illustration that we present, the map f such that Q = f (W) is known numerically, that is, for w j d given, we have calculated q j d = f (w j d ) with the boundary value problem to generate the training set. This situation is particular and is not that of the general framework that we have given ourselves. Nevertheless, this particular situation allows us to use another method to qualify the quality of the probability measure of W post estimated with the constrained learned set, as follows. The proposed constrained-learned-set algorithm allows for computing the realizations {η 1 post , . . . , η N post } of the posterior random variable H post , which constitute the points of the constrained learned set. Using Eq. (6.44), the realizations q 1 post , . . . , q N post and w 1 post , . . . , w N post are computed by the equations Q post = q + [Φ q ] [κ] 1/2 H post and W post = w + [Φ w ] [κ] 1/2 H post . We can then compare Q post = q + [Φ q ] [κ] 1/2 H post with Q QA = f (W post ) in which W post = w + [Φ w ] [κ] 1/2 H post and where mapping f is evaluated with the computational model.

ω 3 = 0.095774. (a) (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 ; θ)] 11 (b) (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 ; θ)] 12 (c) (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 ; θ)] 44
a) (ω 1 , ω 2 ) → U 1 (ω 1 , ω 2 , ω 3 ; θ) (b) (ω 1 , ω 2 ) → U 2 (ω 1 , ω 2 , ω 3 ; θ) (c) (ω 1 , ω 2 ) → U 3 (ω 1 , ω 2 , ω 3 ; θ)
For N r = 100, the convergence of Q QA = f (W post ) with respect to N d has been analyzed by studying, for k = 1, 2, 3, the mean-square norm |||Q QA obs,k |||= {E{(Q QA obs,k ) 2 }} 1/2 of random component Q QA obs,k of Q QA (which depends on N d ). For k = 1, 2, 3, Table 1 yields the values of |||Q QA obs,k |||. The expected convergence can be viewed with respect to N d (this result is consistent with the convergence of the pdf's shown in Fig. 10). In addition to the convergence analysis of the mean-square norm with respect to N d , Fig. 10 displays the probability density functions of the random variables Q obs,1 , Q obs,2 , and Q obs,3 , for Q estimated with the training set, for Q QA = f (W post ) estimated with W post , and for Q targ corresponding to the reference. Similarly to the convergence of the mean-square norm, this figure shows a clear convergence with respect to N d . For N d = 100 and N r = 20 or even N r = 100, compared to the reference, the posterior Q QA evaluated with W post thanks to the knowledge of f , is less good than the prediction of the posterior Q post . This is mainly due to the use of the reduced representation for a problem in high dimension. Fig. 10 shows that the prediction can be improved by increasing the value of N d , that is to say, by increasing the value of ν, which requires to increase the number of points in the training set and consequently, which can induce potential difficulties if the numerical cost for constructing each point of the training set is high. However, the presented numerical illustration shows that the proposed method allows for integrating,with a good quality, a target set of realizations (i.e. data) in a supervised model, which is defined only by a small number of points in a training set and for which the target set of realizations are specified only for the quantities of interest (output) and not for the control variable (input). Finally, it should be noted that, when the training set is generated using a stochastic boundary value problem, there is also another method as we have proposed and validated in [START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]. It consists, for the generation of the constrained learned set, to introduce an additional scalar constraint to minimize the norm of the residue of the partial differential equations of the boundary value problem. This procedure can be implemented without difficulty in the methodology presented in this paper, involving only one additional component in the vector-valued function h C and the vector b c .

Conclusion

When the target set does not consist of realizations of the QoI of an uncertain/stochastic computational model but is made up of statistical moments of this QoI, the Bayesian inference cannot easily be used and the KLDMP allows for estimating the posterior probability measure given its prior probability measure and the constraints related to the statistical moments. Until now there were no general methods to impose constraints defined by realizations in the KLDMP, without estimating the statistical moments from these realizations. In this paper, we have presented a novel functional approach that makes it possible to take into account a target set of realizations in the KLDMP for constructing a posterior probability measure from a prior probability measure defined by a given training set of realizations. The framework is the one for which the statistical moments associated with these realizations are not estimated. The method consists in constructing a functional constraint, defined by the target set of realizations, which is formulated in a Hilbertian framework as the equality of the Fourier transforms of the posterior probability measure and the target probability measure. We have presented the construction and the mathematical analysis of the finite representation of the weak formulation of this functional constraint. The proven mathematical properties have and n w = 420 000. The target set D targ = {q 1 targ , . . . , q N r targ } is generated using the stochastic computational model with another values of the parameters (see Appendix A.3).

1. 5 .

 5 Notations for sets of matrices and their norm M n,m : set of the (n × m) real matrices. M n : set of the square (n × n) real matrices. M + n : set of the positive-definite (n × n) real matrices. M +0 n : set of the positive (n × n) real matrices. [I n

  xN d ] be the matrix in M n x ,N d and let [Φ] [S ] [Φ] T = [ x] be the economy size SVD (thin SVD [54]) of matrix [ x] in which [Φ] T is the transpose of matrix [Φ]. The diagonal entries of [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 [

  Eqs. (3.4) and (3.5) are then deduced from Eqs. (3.1) and (3.

  .16) PROOF. (Theorem 1). (i) The change of variable vv ′ = u ′ in the right-hand side of Eq. (4.12) yields

  and due to Eq. (4.13), we have F(φ) = | F(φ) | < +∞. The proof of Eq. (4.[START_REF] Bilionis | Bayesian uncertainty propagation using Gaussian processes[END_REF]) is similar to the proof of Eq. (4.15) by using also Eq. (3.3) with the second equation in Eq.(3.5)

∼φ

  = 0 dv-almost everywhere. We then have proven the Lemma. Definition 4 (Finite representation of the functional constraint). Using the family { ψr , r = 1, . . . , N r } in C 0 ∩ L 1 (R ν , C) defined in Lemma 5 (see Eq. (5.14), taking into account Lemma 6 and using Definition 2 of the weak formulation of the constraint (see Eqs.

  (4.20) to (4.22)), restricted to φ ∈ H 1, µ ⊂ H 1 , the finite representation of the constraint is written as E{h

  (5.16) to (5.19), we introduce the function (m targ , σ targ ) → J(m targ , σ targ ) = ∥ E{h c (H)} -b c ∥ for m targ ∈ [-3 , 3] and σ targ ∈ [0.1 , 2.3]. For ν = 100, the values of function J are estimated using the realizations {η j d , j = 1, . . . , N d } of H and the realizations {η r targ , r = 1, . . . , N r } of H targ with N d = 1000 and N r = 100. Fig. 1 displays the graph of function (m targ , σ targ ) → J(m targ , σ targ ). It can be seen that J is effectively minimum in the region centered at point (m targ = 0, σ targ = 1) (what was expected because H and H targ are isonomic when m targ = 0 and σ targ = 1) while J is larger when H targ is not isonomic to H, that is to say for m targ ̸ = 0 and/or σ targ ̸ = 1.

Figure 1 :

 1 Figure 1: Graph of function (m targ , σ targ ) → J(m targ , σ targ ) = ∥E{h c (H)} -b c ∥ that illustrates the numerical behavior of the finite representation of the weak formulation of the functional constraint used to identify the posterior probability measure of H from the target set D targ .

Figure 2 :

 2 Figure 2: For the training set, example of one realization θ ∈ Θ of the components (1, 1), (1, 2), and (4, 4) of the random field (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 )] in the plane ω 3 = 0.095774.

(

Figure 3 :

 3 Figure 3: For the training set, example of one realization θ ∈ Θ of the components 1, 2, and 3 of the random field (ω 1 , ω 2 ) → U(ω 1 , ω 2 , ω 3 ) in the plane ω 3 = 0.095774.

Figure 4 :

 4 Figure 4: Principal component analysis of the training set performed for N d = 100. Graph of the eigenvalues α → κ α of [ C X ] and error function ν → err PCA (ν ; N d ) defined by Eq. (2.3).

7. 3 .

 3 Numerical values of the algorithm parameters for computing the constrained learned set.

7. 4 . 7 . 5 .

 475 Iterative algorithm for computing λ sol .For computing the solution λ sol of the Lagrange multiplier, Algorithm 1 is used. For N d = 100 and N r = 20, Fig.5displays the graph of the relaxation factor i → α relax (i) and the graph of the error function i → err(i) defined by Eq. (6.36). It can be seen a fast convergence of the algorithm as a function of the iteration number. Fig. 6 shows the graph of function r → b c r defined by Eq. (5.19), representing the components of vector b c = (b c 1 , . . . , b c N r ), and the graph of function r → λ sol r , representing the components of vector λ sol = (λ sol 1 , . . . , λ sol N r ) defined by Eq. (6.35). Posterior probability measure of Q post estimated with the constrained learned setFor N d = 100, the convergence with respect to N r of the posterior probability measure of Q post (N r ) estimated with the constrained learned set has been analyzed by studying, for k = 1, 2, 3, the mean-square norm|||Q obs,k (N r )|||= {E{Q obs,k (N r ) 2 }} 1/2 of random component Q obs,k (N r ) of Q post (N r ) (which depends on N r ). For k = 1, 2, 3,Fig. 7 shows the graph of function N r → |||Q obs,k (N r )||| as well as the corresponding value for the training set and for the reference, which are both independent of N r . This figure shows the good convergence with respect to N r , which is reached for N r = 20. Always for N d = 100, Fig. 8 (a), (b), and (c) related to N r = 20, and Fig. 8 (d), (e), and (f) related

Figure 5 :

 5 Figure 5: Iterative algorithm for calculating λ sol with the constrained learned set for N d = 100 and N r = 20. Graph of the relaxation factor i → α relax (i) defined in Algorithm 1 and graph of the error function i → err(i) defined by Eq. (6.36).

Figure 6 :

 6 Figure 6: Constraint and optimal Lagrange multiplier estimated with the constrained learned set for N d = 100 and N r = 20. Graph of function r → b c r defined by Eq. (5.19) and graph of function r → λ sol r with λ sol defined by Eq. (6.35).

Figure 7 :

 7 Figure 7: Convergence of the mean-square norm with respect to N r for N d = 100. Graph of N r → |||Q obs,1 (N r )|||, |||Q obs,2 (N r )|||, and |||Q obs,3 (N r )||| estimated with the training set (top black line), with the constrained learned set (the posterior) (blue line), and the reference (red thick line).

7. 6 .

 6 Posterior probability measure of W post estimated with the constrained learned set

Fig. 9

 9 Fig. 9 is related to the standard deviation fields (ω 1 , ω 2 ) → σ G11 (ω 1 , ω 2 , ω 3 ), σ G12 (ω 1 , ω 2 , ω 3 ), and σ G44 (ω 1 , ω 2 , ω 3 ) in the plane ω 3 = 0.095774 of the components (1, 1), (1, 2), and (4, 4) of the random field (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 )] for the training set with N d = 100, for the transformation of the posterior W post (see Appendix A.2) computed with the constrained learned set for which N d = 100 and N r = 20, and finally, for the reference. Since there is no target for the control variable, we cannot directly compare W post (that is to say the random field {[G post (ω)], ω ∈ Ω}) with

(

Figure 8 :

 8 Figure 8: For N d = 100, for N r = 20 (Figs. (a), (b), and (c)) and for N r = 100 (Figs. (d), (e), and (f)), probability density functions of the random variables Q obs,1 , Q obs,2 , and Q obs,3 , estimated with the training set (black line), with the constrained learned set (the posterior) (blue line), and the reference (red thick line).

Figure 9 :

 9 Figure 9: For the training set with N d = 100 (Figs. (a, b, and c), for the posterior estimated with the constrained learned set with N d = 100 and N r = 20 (Figs. (d, e, and f), and for the reference (Figs. (g,h, and i), standard deviation fields (ω 1 , ω 2 ) → σ G11 (ω 1 , ω 2 , ω 3 ), σ G12 (ω 1 , ω 2 , ω 3 ), and σ G44 (ω 1 , ω 2 , ω 3 ) in the plane ω 3 = 0.095774 of the components (1, 1), (1, 2), and (4, 4) of the random field (ω 1 , ω 2 ) → [G(ω 1 , ω 2 , ω 3 )].

  .6) Hypothesis 1 (Existence, regularity, and integrability of the density of P H (resp. P H targ )). It is assumed that the probability measure P H (resp. P H targ ) on R ν of the R ν -valued random variable H (resp. H targ ) is defined by a density η → p H (η) (resp. η → p H targ (η)) with respect to the Lebesgue measure dη, such that p H

  Algorithm for calculating λ sol and DH post = {η 1 post , . . . , η N 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 r 3: for i = 1 : i max doD H λ i = {η 1 λ i , .. . , η N λ i } from Eq. (6.43), using Eqs. (6.37) to (6.42) and D H λ i-1 (D H λ 0 not used for i = 1)

	4:	for ℓ = 1 : N (loop in parallel computation) do
	5:
	6:	end for
	7:	∇ λ Γ(λ i ) and [Γ ′′ (λ i )] using Eqs. (6.19) and (6.20), and D H λ i
	8:

1/2 H post . (6.44) Algorithm 1 post }. 1:

  using Eq. (6.34) with a relaxation factor α relax (i) ∈]0 , 1]D H λ i-1 ← D H λ i 12:end for 13: λ sol = λ i sol , i sol = arg min i err(i) from Eq. (6.35) 14: D H post ← D H λ sol

	10:	λ i ← λ i+1
	11:

Table 1 :

 1 For N r = 100, convergence of the mean-square norm of |||Q QA obs,1 |||, |||Q QA obs,2 |||, and |||Q QA obs,3 ||| as a function of N d . ||| 7.59 × 10 -3 8.15 × 10 -3 8.12 × 10 -3 8.20 × 10 -3 |||Q QA obs,2 ||| 2.44 × 10 -2 2.60 × 10 -2 2.63 × 10 -2 2.67 × 10 -2 |||Q QA obs,3 ||| 5.11 × 10 -3 5.46 × 10 -3 6.35 × 10 -3 5.56 × 10 -3

	N d	100	200	300	400
	|||Q QA obs,1

Preprint submitted to Computational Statistics, SpringerDecember 23, 2022

(a) pdf q → p Q obs,1 (q) for N d = 100 (b) pdf q → p Q obs,2 (q) for N d = 100 (c) pdf q → p Q obs,3 (q) for N d = 100 (d) pdf q → p Q obs,1 (q) for N d = 200 (e) pdf q → p Q obs,2 (q) for N d = 200 (f) pdf q → p Q obs,3 (q) for N d = 200

(g) pdf q → p Q obs,1 (q) for N d = 300 (h) pdf q → p Q obs,2 (q) for N d = 300 (i) pdf q → p Q obs,3 (q) for N d = 300 (j) pdf q → p Q obs,1 (q) for N d = 400 (k) pdf q → p Q obs,2 (q) for N d = 400 (l) pdf q → p Q obs,3 (q) for N d = 400 allowed us to show the existence and the uniqueness of the posterior probability measure constructed by using the KLDMP. The numerical aspects have been detailed in order to facilitate the implementation of the algorithms. This method is a novel theoretical contribution for which a nontrivial and in high-dimensional numerical illustration has been presented in the field of engineering sciences and which demonstrates the efficiency and the robustness of the proposed approach.

Appendix A. Generation of the training set, target set, and numerical values of the parameters

The training set D d = {x 1 , . . . , x N d } with x j = (q j d , w j d ) ∈ R n x = R n q ×R n w is made up of N d independent realizations of random variable X = (Q, W), which are generated by using a stochastic computational model corresponding to the finite element discretization of a stochastic elliptic boundary value problem for which n x = 430 098, n q = 10 098,

The outward unit normal to ∂Ω is denoted by n = (n 1 , n 2 , n 3 ). We use the usual convention of summation on repeated Latin indices. Domain Ω is occupied by a heterogeneous and anisotropic elastic random medium for which the elastic properties are defined by the fourthorder tensor-valued non-Gaussian random field

The stochastic boundary value problem is written, for k = 1, 2, 3 and almost surely, as

)

in which the stress tensor σ is related to the strain tensor ε by ε nq = (∂U n /∂ω q + ∂U q /∂ω n )/2 by the constitutive equation, σ km (ω) = A kmnq (ω) ε nq (U(ω)). For k = 1, 2, 3, the applied stresses p k on Γ 1 are defined as follows:

1}. Using the matrix representation in Voigt notation, the random elasticity field is rewritten, for k, m, n, and q in {1, 2, 3}, as [ A(ω)] ij = A kmnq (ω) with i = (k, m) with 1 ≤ k ≤ m ≤ 3 and j = (n, q) with 1 ≤ n ≤ q ≤ 3 in which indices i and j belong to {1, . . . , 6}. The M + 6 -valued random field {[ A(ω)], ω ∈ Ω} is a non-Gaussian, second order, and statistically homogeneous. Its mean function is the given ω-independent matrix

6 corresponding to a homogeneous isotropic elastic material whose Young modulus is 10 10 N/m 2 and Poisson coefficient 0.15 (note that the fluctuations around the mean are those of a heterogeneous anisotropic elastic material). The non-Gaussian M + 6valued random field {[ A(ω)] , ω ∈ Ω} is constructed using the stochastic model [START_REF] Soize | Non Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[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] Soize | Uncertainty Quantification[END_REF] of random elasticity fields for heterogeneous anisotropic elastic media that are isotropic in statistical mean and exhibit anisotropic statistical fluctuations, for which the parameterization consists of spatial-correlation lengths and of a positive-definite lower bound. The random field

where ϵ is a given positive number (which can be chosen arbitrarily small), and where {[G(ω)], ω ∈ R 3 } is a M + 6 -valued random field (by construction), defined on (Θ, T , P), indexed by R 3 . Then [G] is homogeneous, mean-square continuous, and such that E{[G(ω))]} = [I 6 ] for all ω ∈ R 3 . Note that the lower bound ϵ [ A ]/(1 + ϵ) used in Eq. (A.5) could be replaced by a more general lower bound [A b ] in M + 6 as proposed 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 | Uncertainty Quantification[END_REF]. For all ω fixed in R 3 , the M + 6 -valued random variable [G(ω)] has been constructed by using the Maximum Entropy Principle under the following available information, E{[G(ω)]} = [I 6 ] and E{log(det[G(ω)])} = b G with |b G | < +∞, which has been introduced in order that the random matrix [G(ω)] -1 (that exists almost surely) be such that )] of 6 × (6 + 1)/2 = 21 independent normalized Gaussian real-valued random variables denoted by {G mn (ω), 1 ≤ m ≤ n ≤ 6} and such that E{G mn (ω)} = 0 and E{G mn (ω) 2 } = 1. The spatial correlation structure of random field {[G(ω)], ω ∈ R 3 } is introduced by considering 21 independent real-valued random fields {G mn (ω), ω ∈ R 3 } for 1 ≤ m ≤ n ≤ 6, corresponding to 21 independent copies of a unique normalized Gaussian homogeneous mean-square continuous real-valued random field {G(ω), ω ∈ R 3 } whose normalized spectral measure is given and has a support that is controlled by three spatial correlation lengths L c1 = L c2 = L c3 = 0.4. Note that this Gaussian field G can be replaced by a non-Gaussian field for taking into account uncertainties in the spectral measure [START_REF] Soize | Stochastic elliptic operators defined by non-Gaussian random fields with uncertain spectrum[END_REF]. The constant b G is eliminated in favor of a hyperparameter δ G > 0, which allows for controlling the level of statistical fluctuations of [G(ω)], defined by

, which is independent of ω and such that δ G = 0.6.

Appendix A.2. Stochastic computational model for generating the training set D d and observed quantities of interest

The stochastic boundary value problem defined by Eqs. (A.1) to (A.4) is discretized by the finite element method. Domain Ω is meshed with 50×10×5 = 2 500 finite elements using 8-nodes finite elements. There are 3 366 nodes and 10 098 dofs (degrees of freedom). The displacements are locked at all the 66 nodes belonging to surface Γ 0 and therefore, there are 198 zero Dirichlet conditions. There are 8 integration points in each finite element. Consequently, there are N p = 20 000 integration points ω 1 , . . . , ω N p . The R n w -valued random variable W is generated as follows. For all p = 1, . . . , N p , let [G log p ] = log M ([G(ω p )]) ∈ M 6 in which log M is the logarithm of positive-definite matrices. The R n w -valued random variable W is then defined as the vector that is the reshaping of the upper triangular part of the N p matrices { [G log p ], p = 1, . . . , N p }.We then have n w = 21 × N p = 420 000. The finite element discretization of random field {U(ω), ω ∈ Ω } is the R n q -valued random variable Q with n q = 10 098. Consequently X = (Q, W) is a random variable with values in R n x with n x = n q + n w = 430 098. The stochastic computational model is then represented by a stochastic linear matrix equation that is solved by using the Monte Carlo numerical simulation method yielding the training set

), the computed realizations being independent. For studying the convergence properties, the considered values of N d are N d ∈ {100, 200, 300, 400}.

The components of the quantity of interest Q, which will be observed for presenting the results, are the 3 components denoted by Q obs,1 , Q obs,2 , and Q obs,3 that correspond to the 3 dofs along directions ω 1 , ω 2 , and ω 3 of the finite element node of coordinates (0, 0, 0.1) (located at top corner in which the displacements are significant and result from tension, torsion, and two bendings contributions).

Appendix A.3. Target set of realizations

The target set D targ = {q 1 targ , . . . q N r targ } is generated using the stochastic boundary value problem defined in Section Appendix A.1 for which the elasticity matrix [ A targ ] is the one of a homogeneous and isotropic elastic material with a Young modulus 9 × 10 9 N/m 2 and a Poisson coefficient ν = 0.15. The level of statistical fluctuations of the random field {G targ (ω), ω ∈ R 3 } is δ targ G = 0.3. In order to analyze the convergence with respect to N r , we have considered, in consistency with the values of N d , the intervals N r ∈ [50 , N targ ] with N targ ∈ {100, 200, 300, 400}.

Conflict of interest

The author declares that he has no conflict of interest.