C Soize
email: christian.soize@univ-paris-est.fr

Polynomial chaos expansion of a multimodal random vector

Keywords: multimodal probability distribution, polynomial chaos, random vector, high dimension AMS subject classifications. 62H10, 62H12, 62H20, 60E10, 60E05, 60G35, 60G60

A methodology and algorithms are proposed for constructing the polynomial chaos expansion (PCE) of multimodal random vectors. An algorithm is developed for generating independent realizations of any multimodal multivariate probability measure that is constructed from a set of independent realizations using the Gaussian kernel-density estimation method. The PCE is then performed with respect to this multimodal probability measure, for which the realizations of the polynomial chaos are computed with an adapted algorithm. Finally, a numerical application is presented for analyzing the convergence properties.

analyzing the convergence properties. This new class of algorithms for the multimodal case can be useful in the context of uncertainty quantification for direct and inverse problems, and in particular, for the approaches devoted to dimension reduction in chaos expansions for nonlinear coupled problems, when an iterative solver is used (see for instance [START_REF] Arnst | Dimension reduction in stochastic modeling of coupled problems[END_REF][START_REF] Arnst | Hybrid sampling/spectral method for solving stochastic coupled problems[END_REF][START_REF] Arnst | Reduced chaos expansions with random coefficients in reduced-dimensional stochastic modeling of coupled problems[END_REF]).

2. Construction of a representation of a multimodal random vector in high dimension from a set of realizations. In the first part of this paper, we propose a construction of a stochastic model of a multimodal random vector X with values in R N , defined on a probability space (Θ, T , P) using only a set of ν ≫ 1 independent realizations of X. The construction consists (i) in introducing the usual empirical estimations of the mean vector and the covariance matrix of X, (ii) in constructing a reduced-order statistical model X (n) (with values in R N ) of X (with dimension n < N ) using the classical principal component analysis, and yielding a multimodal reduced random vector H with values in R n (which is assumed to be in high dimension, that is to say, with N > n ≫ 1), (iii) in constructing a multimodal probability density function η → p H (η) on R n of H introducing an adapted nonparametric statistical estimator, (iv) in developing a generator of independent realizations of H that follows the multimodal probability distribution p H (η) dη on R n , (v) for fixed n, in computing a sequence { w M } M ≥1 } of statistical estimations of {w(H)} using this generator, in which w is a given measurable mapping defined on R n . At the end of this section, a numerical investigation is presented in order to illustrate the construction and the algorithms.

Data description and usual empirical estimations of second-order moments.

Let X = (X 1 , . . . , X N ) be a R N -valued second-order random vector defined on a probability space (Θ, T , P), whose probability distribution is represented by an unknown probability density function x → p X (x) with respect to the Lebesgue measure dx on R N . Let E be the mathematical expectation. It is assumed that ν (with ν ≫ 1) independent realizations x exp,1 , . . . , x exp,ν of X are known (coming from experimental data or from numerical simulations). Let m X and [ C X ] be the empirical estimations of the mean vector m X = E{X} and the covariance matrix [C X ] = E{(Xm X ) (Xm X ) T }, such that

m X = 1 ν ν ℓ=1
x exp,ℓ , [ C X ] = 1 ν -1 ν ℓ=1 (x exp,ℓm X ) (x exp,ℓm X ) T .

(2.1)

Note that [ C X ] can be written as

[ C X ] = [ R X ] - ν + 1 ν -1 m X m T X , [ R X ] = 1 ν -1 ν ℓ=1 x exp,ℓ (x exp,ℓ ) T , (2.2)
in which [ R X ] is the estimation of the second-order moment matrix, [R X ] = E{X X T }.

2.2.

Reduced-order statistical model of X. The eigenvalues λ X,1 ≥ λ X,2 ≥ . . . ≥ λ X,N ≥ 0, and the associated orthonormal eigenvectors φ 1 X , . . . , φ N X , such that (φ i X ) T φ j X = δ ij (the Kronecker symbol), are such that [ C X ] φ i X = λ X,i φ i X . The principal component analysis allows a reduced-order statistical model, X (n) , of X to be constructed, such that

X (n) = m X + n i=1 λ X,i H i φ i X .
(2.3)

Consequently, random vector X is such that X = X (N ) , and the value of n is fixed in [1, N ] for that err(n) ≤ ε in which ε is such that 0 ≤ ε ≪ 1, and where err is the error function defined by

err(n) = 1 - n i=1 λ X,i tr[ C X ] . (2.4)
Let H = (H 1 , . . . , H n ) be the R n -valued second-order random variable. The ν independent realizations η exp,1 , . . . , η exp,ν , with η exp,ℓ = (η exp,ℓ 1 , . . . , η exp,ℓ n ) ∈ R n , of the R n -valued random variable H are computed, for ℓ = 1, . . . , ν and i = 1, . . . , n, η exp,ℓ i = 1 λ X,i (x exp,ℓm X ) T φ i X .

(2.5)

By construction, it can easily be verified that

m H = 1 ν ν ℓ=1 η exp,ℓ = 0 , [ R H ] = 1 ν -1 ν ℓ=1 η exp,ℓ (η exp,ℓ ) T = [ I n ] .
(2.6)

Construction of the multimodal probability density function of random vector H.

As explained in the introduction of Section 2, the probability density function p H with respect to the Lebesgue measure dη on R n of the R n -valued second-order random variable H = (H 1 , . . . , H n ), must be constructed. The unknown probability density function p X has been assumed to be multimodal. Due to the reduced representation introduced by Eq. (2.3), probability density function p H (that differs from p X ), could be not multimodal. The objective of the present work is to develop a methodology adapted to the case for which p H is multimodal, and consequently, it is assumed that p H can be multimodal. It should be noted that the method proposed is very general and can be used for a probability density function p H that is or that is not multimodal. We propose to choose for p H its estimation carried out by using the Gaussian kernel-density estimation method on the basis of the knowledge of the ν independent realizations, η exp,1 , . . . , η exp,ν computed with Eq. (2.5). A modification of the classical Gaussian kernel-density estimation method [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF] is used in order that Eq. (2.6), [ R H ] = [ I n ], be preserved. The positive valued function p H on R n is then defined, for all η in R n , by

p H (η) = 1 ν ν ℓ=1 µ n, sn ( s n s n η exp,ℓ -η) , (2.7)
in which µ n, sn is the positive function from R n into ]0 , +∞[ defined, for all η in R n , by

µ n, sn (η) = 1 ( √ 2π s n ) n exp{- 1 2 s 2 n η 2 n } , (2.8)
in which η 2 n = η 2 1 + . . . + η 2 n , and where the positive parameters s n and s n are defined by

s n = 4 ν(2 + n) 1/(n+4) , s n = s n s 2 n + ν-1 ν . (2.9)
Parameter s n is the usual multidimensional optimal Silverman bandwidth (in taking into account that the empirical estimation of the standard deviation of each component is unity, and parameter s n has been introduced in order that the second equation in Eq. (2.6) holds. It should be noted that, for n fixed, parameters s n and s n go to 0 + , and s n /s n goes to 1 -, when ν goes to +∞. Using Eqs. (2.7) to (2.9), it can easily be verified that

E{H} = R n η p H (η) dη = s n s n m H = 0 , (2.10
) (iii) In the framework of the proposed developments, no assumptions (in particular concerning the support) are introduced concerning the unknown multivariate probability distribution that has to be estimated from a set of data by the nonparametric statistics using the kernel-density estimation method (kernel smoothing). Consequently, a first arbitrary choice must be done concerning the kernel smoothing, and the kernel proposed is the multivariate Gaussian kernel for centered and uncorrelated random variables (see Eq. (2.6)) (see Eq. (2.6)). The second arbitrary choice concerns the bandwidth. Again, since no information is available, the multidimensional Silverman bandwidth that is optimal for a Gaussian distribution is chosen. It should be noted that this choice for the bandwidth is consistent with the usual choice of the empirical estimators used for estimating the mean vector and the covariance matrix, which are optimal (unbiased, efficient and consistent estimators) for a Gaussian distribution. In addition, the methodology proposed in the paper is general and is independent of the choice of the kernel smoothing, provided that the assumptions used by the generator of realizations are satisfied. Therefore, if additional information are available concerning the unknown multivariate probability distribution that has to be estimated from the set of data, the Gaussian kernel µ n, sn (η) and the Silverman bandwidth s n can be replaced by a kernel and bandwidth that are better adapted. Nevertheless, it should be noted that the nonparametric estimation yields a probability density function p H that is only used for constructing a random vector H (whose probability distribution is defined by this pdf), which is used as a germ for performing the PCE in order to increase the speed of convergence of the truncated PCE.

E{H H T } = R n η η T p H (η) dη = s 2 n [ I n ] + ( s n s n ) 2 (ν -1) ν [ R H ] = [ I n ] . ( 2

Remark 2.

For i fixed in {1, . . . , n}, the probability density function p H i on R of the random variable H i is calculated in integrating Eq. (2.7) over R n-1 , and yields

p H i (η i ) = 1 ν ν ℓ=1 µ 1, sn ( s n s n η exp,ℓ i -η i ) .
(2.12)

For i and j fixed in {1, . . . , n}, the joint probability density function p H i H j on R 2 of the random variables H i and H j , is calculated in integrating Eq. (2.7) over R n-2 , and yields

p H i ,H j (η i , η j ) = 1 ν ν ℓ=1 µ 2, sn ( s n s n η exp,ℓ,ij -η ij ) , (2.13)
in which η exp,ℓ,ij = (η exp,ℓ i , η exp,ℓ j

) and η ij = (η i , η j ) belong to R 2 .

Generator for multimodal random vector H.

Let w be a mapping from R n into an Euclidean space such that w(H) is a second-order random variable. The estimation of E{w(H)} = R n w(η) p H (η) dη requires a generator of independent realizations of random vector H for which the multimodal probability distribution p H (η) dη is defined by Eq. (2.7) with Eqs. (2.8) to (2.9). Such a generator can be performed using the Markov Chain Monte Carlo method (MCMC) [START_REF] Kaipio | Statistical ans Computational Inverse Problems[END_REF][START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF]. The transition kernel of the homogeneous Markov chain of the MCMC method can be constructed using the Metropolis-Hastings algorithm [START_REF] Metropolis | The Monte Carlo method[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] (that requires the definition of a good proposal distribution), the Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distribution and the Bayesian distribution of images[END_REF] (that requires the knowledge of the conditional distribution) or the slice sampling [START_REF] Neal | Slice sampling[END_REF] (that can exhibit difficulties related to the general shape of the probability distribution, in particular for multimodal distributions). In general, these algorithms are efficient, but can also be not efficient if there exist attraction regions which do not correspond to the invariant measure under consideration and tricky even in high dimension. These cases cannot easily be detected and are time consuming. The method 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] is very robust, has recently been applied with success for complex problems in high dimension [START_REF] Batou | Calculation of Lagrange multipliers in the construction of maximum entropy distributions in high stochastic dimension[END_REF][START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF] and is reused hereinafter. It looks similar to the Gibbs approach but corresponds to a more direct construction of a random generator of realizations for random variable H whose probability distribution is p H (η) dη and is multimodal. The difference between the Gibbs algorithm and the proposed algorithm is that the convergence in the proposed method can be studied with all the mathematical results concerning the existence and uniqueness of Itô stochastic differential equation. In addition, a parameter is introduced which allows the transient part of the response to be killed in order to get more rapidly the stationary solution corresponding to the invariant measure. Thus, following [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], the construction of the transition kernel by using the detailed balance equation is replaced by the construction of an Itô Stochastic Differential Equation (ISDE), which admits p H (η) dη (defined by Eqs. (2.7)to (2.9)) as a unique invariant measure. The ergodic method or the Monte Carlo method can be used for estimating E{w(H)}.

It should be noted that the main ideas presented in this paper are not related to a specific MCMC algorithm for constructing a set of realizations. The alternative MCMC algorithm proposed hereinafter can be replaced by any traditional MCMC algorithm. Nevertheless, this alternative algorithm is very robust and very rich in terms of control based on the use of mathematical results for the Itô stochastic differential equations.

Interpretation of the multimodal probability distribution p H as the invariant measure of an It ô stochastic differential equation (ISDE

c n = 1 ( √ 2π s n ) n , Φ(η) = -log{q(η)} , (2.15)
in which η → q(η) is the continuously differentiable function from R n into R + * =]0 , +∞[ defined by

q(η) = 1 ν ν ℓ=1 exp{- 1 2 s 2 n s n s n η exp,ℓ -η 2 } , (2.16)
and where s n and s n are given by Eq. (2.9). It can then be deduced that η → Φ(η) is a continuously differentiable function on R n .

Let {(U(r), V(r)), r ∈ R + } be the Markov stochastic process defined on the probability space (Θ, T , P), indexed by R + = [0 , +∞[, with values in R n × R n , satisfying, for all r > 0, the following Itô stochastic differential equation dU(r) = V(r) dr , (2.17)

dV(r) = -∇ u Φ(U(r)) dr - 1 2 f 0 V(r) dr + f 0 dW(r) , (2.18)
with the initial condition U(0) = u 0 , V(0) = v 0 a.s. , (2.19) in which u 0 and v 0 are given vectors in R n (that will be taken as zero in the application presented later), and where W = (W 1 , . . . , W n ) is the normalized Wiener process defined on (Θ, T , P) indexed by R + with values in R n . The matrix-valued autocorrelation function [R W (r, r ′ )] = E{W(r) W(r ′ ) T } of W is then written as [R W (r, r ′ )] = min(r, r ′ ) [I n ]. In Eq. (2.18), the free parameter f 0 > 0 will allow a dissipation term to be introduced in the nonlinear second-order dynamical system (formulated in the Hamiltonian form with an additional dissipative term) in order to kill the transient part of the response and consequently, to get more rapidly the stationary solution corresponding to the invariant measure. It can easily be proved that function u → Φ(u):

(i) is continuous on R n , (ii) is such that u → ∇ u Φ(u) is a locally bounded function on R n (i.e. is bounded on all compact sets in R n that is the case because u → Φ(u) is a continuously differentiable function on R n ), and (iii) is such that, inf u >R Φ(u) → +∞ if R → +∞ , (2.20
)

inf u∈R n Φ(u) = Φ min with Φ min ∈ R , (2.21) R n ∇ u Φ(u) p H (u) du < +∞ with p H (u) = c n e -Φ(u) . (2.22)
Under hypotheses (i) to (iii), and using Theorems 4 to 7 in pages 211 to 216 of Ref. [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF], in which the Hamiltonian is taken as H(u, v) = v 2 /2 + Φ(u), and using [START_REF] Doob | Stochastic Processes[END_REF][START_REF] Khasminskii | Stochastic Stability of Differential Equations[END_REF] for the ergodic property, it can be deduced that the problem defined by Eqs. (2.17) to (2.19) admits a unique solution. This solution is a second-order diffusion stochastic process {(U(r), V(r)), r ∈ R + }, which converges to a stationary and ergodic diffusion stochastic process {(U st (r st ), V st (r st )), r st ≥ 0}, when r goes to infinity, associated with the invariant probability measure P st (du, dv) = ρ st (u, v) du dv. The probability density function (u, v) → ρ st (u, v) on R n × R n is the unique solution of the steady-state Fokker-Planck equation with the normalization condition, associated with Eqs. (2.17) and (2.18), and is written (see Propositions 8 and 9 in pages 120 to 123 of Ref. [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF]), as

ρ st (u, v) = c 0 exp{- 1 2 v 2 -Φ(u)} , (2.23)
in which c 0 is the constant of normalization. From Eqs. (2.14) and (2.23), it can be deduced that

p H (η) = R n ρ st (η, v) dv , ∀ η ∈ R n . (2.24)
It can therefore be concluded that random variable H for which the multimodal probability density function is p H , can be defined, for any r st > 0 as

H = U st (r st ) = lim r→+∞ U(r) in probability distribution . (2.25)
As explained above, the free parameter f 0 > 0 introduced in Eq. (2.18), allows a dissipation term to be introduced in the nonlinear dynamical system and consequently, allows the transient response generated by the initial conditions (u 0 , v 0 ) to be rapidly killed in order to get more rapidly the asymptotic behavior corresponding to the stationary and ergodic solution associated with the invariant measure.

Remark 3. Instead of Eq. (2.18), the following equation

dV(r) = -∇ u Φ(U(r)) dr-1 2 f 0 [D 0 ]V(r) dr + √ f 0 [S 0 ] dW(r) could be used, in which [S 0 ] would belong to M n (R)
and where [D 0 ] would be a positive symmetric matrix such that

[D 0 ] = [S 0 ] [S 0 ] T with 1 ≤ rank[D 0 ] ≤ n.
If such an equation were used, then the invariant measure would always be given by Eq. (2.23) (see page 244 of Ref. [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF]). In particular, a diagonal positive-definite damping matrix 1 2 f 0 [D 0 ] could be chosen in order trying to increase the speed of convergence towards the stationary and ergodic solution of the Itô equation. However, in order not to complicate too much setting data parameters of the algorithm, while maintaining good control of the speed of convergence towards the stationary and ergodic solution, the simpler form defined by Eq. (2.18) is proposed.

Discretization scheme of the ISDE.

A discretization scheme must be used for numerically solving the ISDE defined by Eqs. (2.17) to (2.19). For general surveys on discretization schemes for Itô stochastic differential equations, we refer the reader to [START_REF] Kloeden | Numerical Solution of Stochastic Differentials Equations[END_REF][START_REF] Tubaro | Expansion of the global error for numerical schemes solving stochastic differential equation[END_REF][START_REF] Talay | Simulation and numerical analysis of stochastic differential systems[END_REF]. Concerning the particular cases related to Hamiltonian dynamical systems (which have also been analyzed in [START_REF] Talay | Stochastic Hamiltonian system: exponential convergence to the invariant measure and discretization by the implicit Euler scheme[END_REF] using an implicit Euler scheme), we propose to use the Störmer-Verlet scheme, which is a very efficient scheme that preserves energy for nondissipative Hamiltonian dynamical systems (see [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] for reviews about this scheme in the deterministic case, and see [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF] and the references therein for the stochastic case). The Störmer-Verlet scheme has been validated 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], for solving an ISDE of the type defined by Eqs. (2.17) to (2.19), and corresponding to a weakly dissipative Hamiltonian dynamical system. We then propose to reuse hereinafter the Störmer-Verlet scheme 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]. Let M ≥ 1 be an integer. The Itô stochastic differential equation defined by Eqs. (2.17) and (2.18) with the initial condition defined by Eq. (2.19), is solved on the finite interval R = [0 , (M -1) ∆r], in which ∆r is the sampling step of the continuous index parameter r. The integration scheme is based on the use of the M sampling points r k such that r k = (k -1) ∆r for k = 1, . . . , M . The following notations are introduced:

U k = U(r k ), V k = V(r k ), and W k = W(r k ), for k = 1, . . . , M , with U 1 = u 0 , V 1 = v 0 , W 1 = 0 a.s .
(2.26)

Let {∆W k+1 = W k+1 -W k , k = 1, . . . , M -1} be the family of independent Gaussian second- order centered R n -valued random variables such that E{∆W k+1 (∆W k+1 ) T } = ∆r [I n ]. For k = 1, .
. . , M -1, the Störmer-Verlet scheme applied to Eqs. (2.17) and (2.18) is written as

U k+ 1 2 = U k + ∆r 2 V k , V k+1 = 1-b 1+b V k + ∆r 1+b L k+ 1 2 + √ f 0 1+b ∆W k+1 , U k+1 = U k+ 1 2 + ∆r 2 V k+1 , (2.27)
with the initial condition defined by (2.26), where b = f 0 ∆r /4, and where

L k+ 1 2 is the R n -valued random variable such that L k+ 1 2 = -{∇ u Φ(u)} u=U k+ 1 2 .
(2.28)

From Eqs. (2.15) and (2.16), it can be deduced that,

∇ u Φ(u) = - 1 q(u) ∇ u q(u) , (2.29)
∇ u q(u) = 1 s 2 n 1 ν ν ℓ=1 ( s n s n η exp,ℓ -u) exp{- 1 2 s 2 n s n s n η exp,ℓ -u 2 } .
(2.30)

Choosing the parameters for numerical integration of the ISDE.

In this section, we construct the values of the parameters ∆r, M 0 , M, f 0 , u 0 and v 0 , which are used in the discretization scheme of the ISDE, presented in Section 2.4.2. First, we associate with the nonlinear Hamiltonian dynamical system, a linearized diagonal second-order dynamical system in U lin = (U lin 1 , . . . , U lin n ) (the components are not coupled) such that, for all i in {1, . . . , n},

Ü lin i (r) + 1 2 f 0 U lin i (r) + K i U lin i (r) = √ f 0 Ẇi (r)
, in which Ẇ = ( Ẇ1 , . . . , Ẇn ) is the generalized Gaussian white stochastic process (the generalized derivative of W). The behavior of the nonlinear stiffness force F(u) = -L(u), with u = (u 1 , . . . , u n ) and F(u) = (F 1 (u), . . . , F n (u)), can have some fluctuations in the neighborhood of u = 0, and is such that F(0) = 0. Consequently, we cannot calculate K i in writing K i u i = {∂F i (u)/∂u i } u=0 . We propose to replace this equation by an incremental equation on a symmetric interval [-∆ , ∆] for a sufficiently large increment ∆ > 0 (typically, the problem under consideration being normalized to 1, ∆ can, for instance, be chosen as 5), which yields

K i = (F i (u i )-F i (-u i ))/(2∆) in which u i = (u i 1 , . . . , u i n ) with u i j = ∆ δ ij . Let 0 < λ 1 ≤ . . . ≤ λ n be the eigenvalues of matrix [ K ],
and let be ω min = √ λ 1 and ω max = √ λ n .

(i) A first estimation of ∆r is chosen as ∆r 0 = π/(10 ω max ). An oversampling m overs > 1 is introduced to get a sufficient accuracy of the Störmer-Verlet scheme (for instance, m overs = 10), and yields ∆r = ∆r 0 /m overs (a convergence analysis with respect to m overs must be carried out).

(ii) The minimum damping rate, ζ min of the linear second-order dynamical system is such that 2 ζ min ω min = f 0 /2 that yields f 0 = 4 ζ min ω min . The damping rate ζ min is chosen, for instance, as 0.7 to rapidly kill the transient response induced by the initial conditions (that are not distributed following the invariant measure).

(iii) The larger relaxation "time" of the linear second-order dynamical system can be defined as r 0 such that exp{-ζ min ω min r 0 } = ε 0 with ε 0 ≪ 1, which yields r 0 = -log(ε 0 )/(ζ min ω min ) (for instance, the value of ε 0 can be chosen as 1/200). The parameter M 0 is then defined as r 0 = M 0 ∆r that yields M 0 = 1 + fix(r 0 /∆r) in which fix(x) rounds x to the nearest integer towards zero. Value r 0 (and then integer M 0 ) corresponds to reaching the stationary response. Integer M 0 will be used for calculating independent realizations of H.

(iv) The initial conditions are chosen as u 0 = v 0 = 0.

(v-1) Ergodic method. The integer M is defined as M = m ergo M 0 where m ergo ≫ 1 is an integer which has to be chosen in order to reach a reasonable convergence for estimating E{w(H)} using the ergodic method (see Eq. (2.33) that we introduce later). For instance, an initial value for m ergo can be chosen as 200 or 400, but a convergence analysis must be carried out with respect to m ergo . (v-2) Monte Carlo method. An integer M 0 is introduced such that M 0 = i 0 M 0 . The integer i 0 is chosen in order that the sequence {U k , k ≥ M 0 } corresponds to the stationary solution. Taking into account the construction of integer M 0 , the integer i 0 could be chosen to 1 or 2. We then introduce the sequence of integers

{ M ℓ , ℓ = 1, . . . , ν s } such that M ℓ = (1 + ℓ) M 0 . Integer M is then defined by M = M νs = (1 + ν s ) M 0 .
By construction of M 0 , the vectors U M ℓ (θ) and U M ℓ+1 (θ) can be considered as two independent realizations of the random vector U(r) for any fixed r such that r > r 0 (stationary part of the response). The integer ν s is chosen in order to reach a reasonable convergence for estimating E{w(H)} using the Monte Carlo method (see Eq. (2.34) that we introduce later). For instance, an initial value for ν s can be chosen as 200, but a convergence analysis with respect to ν s must be carried out.

Random generator of independent realizations of H.

A random generator of ν s independent realizations, H(θ 1 ), . . . , H(θ νs ), of random vector H whose multimodal probability distribution is p H (η) dη, is constructed as follows. Let {U k (θ), k = 1, . . . , M } be constructed using the algorithm presented in Section 2.4.2. Using the sequence { M ℓ , ℓ = 1, . . . , ν s } defined in Section 2.4.3-(v-2), each independent realization H(θ ℓ ) can then be obtained, for ℓ in {1, . . . , ν s }, as

η sim,ℓ = H(θ ℓ ) = U M ℓ (θ) , M ℓ = (1 + ℓ) M 0 .
(2.31)

Estimating E{w(H)}.

In this section, two methods are proposed for estimating E{w(H)}: (i) and ergodic method, and (ii) a Monte Carlo method. (i) Estimating E{w(H)} with the ergodic method. For any realization θ in Θ, let {(U(r, θ), V(r, θ)), r ≥ 0} be the solution of Eqs. (2.17) to (2.19). Then using the ergodic property, E{w(H)} = R n w(η) p H (η) dη can be estimated (see [START_REF] Doob | Stochastic Processes[END_REF][START_REF] Kloeden | Numerical Solution of Stochastic Differentials Equations[END_REF][START_REF] Talay | Simulation and numerical analysis of stochastic differential systems[END_REF]) by

E{w(H)} = lim R→+∞ 1 R r M 0 +R r M 0 w(U(r, θ)) dr with probability 1 , (2.32)
in which r M 0 = (M 0 -1) ∆r with M 0 a fixed integer greater than 1. The parameter M 0 (estimated as explained in Section 2.4.3-(iii)) allows us to remove the transient part of the response induced by the initial condition. Let {U k (θ), k = 1, . . . , M } be the corresponding realization of the discretized solution constructed as explained in Section 2.4.2, in which M has been estimated in Section 2.4.3-(iv). Thus, the numerical approximation of Eq. (2.32) is written as

E{w(H)} = lim M →+∞ w ER M , w ER M = 1 M -M 0 + 1 M k=M 0 w(U k (θ)) .
(2.33)

(
w(U M ℓ (θ)) , M ℓ = (1 + ℓ) M 0 , (2.34)
2.6. Numerical application.

2.6.1. Data generation for the numerical application. We consider the data description introduced in Section 2.1 for N = 100 and ν = 500. The algorithm that has been used for generating the ν independent realizations x exp,1 , . . . , x exp,ν of random vector X with values in R N , is described in Appendix A. The reader can then simulate the "data description" used in this numerical application.

Defining the optimal values of parameters.

The estimations m X and [ C X ] are computed using Eqs. (2.1) and (2.2). Concerning the construction of the reduced-order statistical model, Fig. 2.1 displays the graph of the error function defined by Eq. (2.4). We choose n = 5 corresponding to an error of 0.089. With the values of the main parameters, N = 100, ν = 500 and n = 5, for the numerical integration of the ISDE (see Section 2.4.3), the optimal value of the damping rate has been found to ζ min = 0.7 that yields f 0 = 5.0565. The optimal value of the relaxation "time" parameter has been found to ε 0 = 1/200 (see Section 2.4.3-(iii)), and ∆r 0 = 0.1670. The others values of the parameters are analyzed in the next section. 2.6.3. Convergence analysis. The convergence analysis is carried out with the ergodic method and with the Monte Carlo method. (i) Ergodic method. The parameters defined in Section 2.6.2 are fixed. In order to analyze the convergence of Eq. (2.33) with respect to M for the ergodic estimation, using the integration scheme defined in Section 2.4.2, two error functions, err 1 (M ) and err 2 (M ) are introduced. The first one is related to the estimation of the mean value E{H} = m H that must be equal to 0 (see Eq. (2.10)). In Eq. (2.33), we choose w(H) = H, and then

err ER 1 (M ) = 1 n n i=1 |{ m H -w ER M } i | , w ER M = 1 M -M 0 + 1 M k=M 0 U k (θ) .
(2.35)

The second one is related to the estimation of the second-order moment matrix E{H H T } = [ R H ] that must be equal to [ I n ] (see Eq. (2.11)). In Eq. (2.33), we choose w(H) = H H T , and then

err ER 2 (M ) = [ R H ] -[ w ER M ] F [ R H ] F , [ w ER M ] = 1 M -M 0 + 1 M k=M 0 U k (θ) U k (θ) T , (2.36)
in which • F is the Frobenius norm. For m overs = 10; 100; and 1, 000 (with ∆r = ∆r 0 /m overs )), we have M 0 = 260; 2, 600; and 26, 000. Figures. 2. A reasonably good accuracy is obtained in using m ergo = 400 that yields M = 104, 000 for which err ER 1 (M ) = 0.060 and err ER 2 (M ) = 0.089. The value m ergo = 10, 000 or 100, 000 (M = 2, 600, 000 or 26, 000, 000) yields err ER 1 (M ) = 0.0080 and err ER 2 (M ) = 0.0264, or err ER 1 (M ) = 0.0021 and err ER 2 (M ) = 0.0061. equal to 0 (see Eq. (2.10)), and yields

err MC 1 (ν s ) = 1 n n i=1 |{ m H -w MC νs } i | , w MC νs = 1 ν s νs ℓ=1 U M ℓ (θ) , M ℓ = (1 + ℓ) M 0 . (2.37)
The second one is related to the estimation of the second-order moment matrix E{H H T } = [ R H ] that must be equal to [ I n ] (see Eq. (2.11)), and yields, for

M ℓ = (1 + ℓ) M 0 , err MC 2 (ν s ) = [ R H ] -[ w MC νs ] F [ R H ] F , [ w MC νs ] = 1 ν s -1 νs ℓ=1 U M ℓ (θ) (U M ℓ (θ)) T .
(2.38)

A numerical study has been performed with f 0 = 5.0565, ε 0 = 1/200, m overs = 10 that yields M 0 = 260, ∆r = 0.01638 and r 0 = 4.26. A sensitivity analysis has be carried out with respect to the value of i 0 that is such that M 0 = i 0 M 0 . The values 1, 2, and 3 for i 0 yield very close results and the differences are not significant. Consequently, i 0 is fixed to the value 1. For these values of the parameters, the graphs of the error functions ν s → err MC 1 (ν s ) and ν s → err MC 2 (ν s ) have been constructed and yield similar results to those shown in Figs. 2.2 and 2.3. A reasonably good accuracy is obtained in using ν s = 400 that yields M = 104, 000 for which err MC 1 (ν s ) = 0.0692 and err MC 2 (ν s ) = 0.1118. The value ν s = 10, 000 (M = 2, 600, 000) yields err MC 1 (ν s ) = 0.0081 and err MC 2 (ν s ) = 0.034, while for ν s = 100, 000 (M = 26, 000, 000) yields err MC 1 (ν s ) = 0.0016 and err MC 2 (ν s ) = 0.0079. In comparing these errors with those given by the ergodic method, it can be concluded that, for a same value of the CPU time, the ergodic method is slightly more efficient. Nevertheless, if we consider the estimation of statistical quantities related to the random variable Ξ that results from a nonlinear transformation h of H such that Ξ = h(H), transformation h has to be evaluated M times with the ergodic method, but only ν s times with the Monte Carlo method. Since M/ν s = M 0 , the gain is large enough when using the Monte Carlo method with respect to the ergodic method.

Estimation of the probability density function of H.

In order to numerically validate the generator (see Sections 2.4.2 and 2.5-(i)) of random vector H whose probability density function p H is defined by Eqs. (2.7) to (2.9), we propose to compare p H with its estimation p H constructed with the generator. This construction used the Gaussian kernel-density estimation method presented in Section 2.3. For all η in R n , we have p H (η) = E{δ 0 (Hη)} that can be rewritten as

p H (η) = lim σn→0 + E{µ n,σn (H -η)} , (2.39)
in which η → µ n,σn (η) is the following Gaussian function from R n into ]0 , +∞[, defined, for all η in R n , by

µ n,σn (η) = 1 ( √ 2π σ n ) n exp{- 1 2σ 2 n η 2 n } , (2.40
p H (η) = 1 M -M 0 + 1 M k=M 0 µ n, σ ER n ( σ ER n σ ER n U k (θ) -η) , (2.41)
in which the positive parameters σ ER n and σ ER n are defined by

σ ER n = 4 (M -M 0 + 1)(2 + n) 1/(n+4) , σ ER n = σ ER n (σ ER n ) 2 + M -M 0 M -M 0 +1 . (2.42)
For i fixed in {1, . . . , n}, the probability density function p H i on R of random variable H i is calculated in using Eq. (2.12). Integrating Eq. (2.41) over R n-1 yields the estimation η i → p H i (η i ) of the probability density function p H i on R,

p H i (η i ) = 1 M -M 0 + 1 M k=M 0 µ 1, σ ER n ( σ ER n σ ER n U k i (θ) -η i ) . ( 2

.43)

For i in {1, . . . , n}, Fig. 2.4 displays the graph of the probability density function p H i calculated with Eq. (2.12), which is compared with p H i calculated with Eq. (2.43) for which m ergo = 400 and m ergo = 10, 000. It can be seen that convergence is reached for m ergo = 10, 000, and that a good approximation is obtained for m ergo = 400. For i and j fixed in {1, . . . , n}, the joint probability density function p H i H j of random variables H i and H j is calculated in using again Eq. (2.13). Using the ergodic method, and in integrating Eq. (2.41) over R n-2 yields the estimation p H i H j of the joint probability density function p H i H j on R 2 , which is written as

p H i H j (η i , η j ) = 1 M -M 0 + 1 M k=M 0 µ 2, σ ER n ( σ ER n σ ER n U k,ij (θ) -η ij ) , (2.44) in which U k,ij (θ) = (U k i (θ), U k j (θ)) and η ij = (η i , η j ) belong to R 2 .
The computation has been carried out for the ten couples of indices, and all the obtained results have the same quality. Fig. 2.5 displays the graph of the joint probability density function p H i ,H j calculated with Eq. (2.13), which is compared with p H i ,H j calculated with Eq. (2.44) for which m ergo = 10, 000. (ii) Monte Carlo method. From Eqs. (2.34) and (2.40), the estimation p H of p H is written, for all η in R n , as

p H (η) = 1 ν s νs ℓ=1 µ n, σ MC n ( σ MC n σ MC n U M ℓ (θ) -η) , M ℓ = (1 + ℓ) M 0 , (2.45)
in which the positive parameters σ MC n and σ MC n are defined by

σ MC n = 4 ν s (2 + n) 1/(n+4) , σ MC n = σ MC n (σ MC n ) 2 + (ν s -1)/ν s . ( 2

.46)

For i fixed in {1, . . . , n}, the estimation p H i of the probability density function on R of R-valued random variable H i is written as

p H i (η i ) = 1 ν s νs ℓ=1 µ 1, σ MC n ( σ MC n σ MC n U M ℓ i (θ) -η i ) (2.47)
For i in {1, . . . , n}, for ν s = 400 and ν s = 10, 000, we obtain a similar result that the one displayed in Fig. 2.4, the convergence being reached for ν s = 10, 000. For i and j fixed in {1, . . . , n}, the estimation (η i , η j ) → p H i H j (η i , η j ) of the joint probability density function on R 2 of the R 2 -valued random variable (H i , H j ) is given by

p H i H j (η i , η j ) = 1 ν s νs ℓ=1 µ 2, σ MC n ( σ MC n σ MC n U M ℓ ,ij (θ) -η ij ) , (2.48) in which U M ℓ ,ij (θ) = (U M ℓ i (θ), U M ℓ j (θ)) and η ij = (η i , η j ) belong to R 2 .
For ν s = 10, 000, a similar result to the one shown in Fig. 2.5 is obtained.

Polynomial chaos expansion of a multimodal random vector.

Let h be a given measurable mapping from R N into R M , and let

Q = (Q 1 , . . . , Q M ) be the R M -valued random variable such that Q = h(X (n) ) , (3.1)
in which X (n) is the multimodal random vector defined by Eq. (2.3). Transformation h is assumed to be such that Q is a second-order random vector. Possibly, the deterministic mapping h transfers the multimodal character of random vector X (n) to random vector Q = h(X (n) ). We are interested in constructing the polynomial chaos expansion (PCE) of random vector Q. If such a PCE of Q was carried out with respect to the polynomial chaos associated with a unimodal random variable, which is usually the choice done (uniform, Gaussian, etc.), then the speed of convergence could be low. Since the multimodal probability distribution of Q is induced by the multimodal probability distribution of H, we propose an alternative approach consisting in using the multimodal random vector H (for which the probability density function p H is defined by Eqs. (2.7) to (2.9)) as the stochastic germ for the polynomial chaos expansion of random vector Q.

3.1.

Constructing independent realizations of Q and reduced-order statistical model. Although the contents of this section is very classic and well known, a brief presentation is given in order to define the notations and the mapping h. Let q sim,1 , . . . , q sim,νs be ν s independent realizations of random vector Q , which are computed, for all ℓ in {1, . . . , ν s }, by

q sim,ℓ = h( m X + n i=1 λ X,i η sim,ℓ i φ i X ) , (3.2)
η sim,ℓ = U M ℓ (θ) , M ℓ = (1 + ℓ) M 0 , (3.3)
in which η sim,1 , . . . , η sim,νs are the ν s independent realizations of the second-order R n -valued random vector H whose probability density function is defined by Eqs. (2.7) to (2.9), and for which the generator of independent realizations is detailed in Section 2.4.4 (see Eq. (2.31)). Similarly to Section 2.1, let m Q and [ C Q ] be the empirical estimations of the mean vector m Q = E{Q} and the covariance matrix

[C Q ] = E{(Q -m Q ) (Q -m Q ) T }, such that m Q = 1 ν s νs ℓ=1 q sim,ℓ , [ C Q ] = 1 ν s -1 νs ℓ=1 (q sim,ℓ -m Q ) (q sim,ℓ -m Q ) T . (3.4) Let λ Q,1 ≥ λ Q,2 ≥ . . . ≥ λ Q,M
≥ 0 be the eigenvalues, and let φ 1 Q , . . . , φ M Q , be the associated orthonormal eigenvectors ((

φ i Q ) T φ j Q = δ ij ) of the eigenvalue problem [ C Q ] φ j Q = λ Q,j φ j Q . The reduced-order statistical model Q (m) of Q is then written as Q (m) = m Q + m j=1 λ Q,j Ξ j φ j Q . (3.5)
The random vector Q is such that Q = Q (M ) , and the value of m is fixed in {1, . . . , M } for that err Q (m) ≤ ε in which ε is any positive real number, and where err Q is the error function defined by

err Q (m) = 1 - m j=1 λ Q,j tr[ C Q ] . (3.6)
Let ξ sim,1 , . . . , ξ sim,νs be the ν s independent realizations of the second-order random vector Ξ = (Ξ 1 , . . . , Ξ m ), computed, for ℓ = 1, . . . , ν s and j = 1, . . . , m, by

ξ sim,ℓ j = 1 λ Q,j (q sim,ℓ -m Q ) T φ j Q . (3.7)
The estimation m Ξ of E{Ξ} and the estimation [ R Ξ ] of E{Ξ Ξ T } are written as

m Ξ = 1 ν s νs ℓ=1 ξ sim,ℓ , [ R Ξ ] = 1 ν s -1 νs ℓ=1 ξ sim,ℓ (ξ sim,ℓ ) T . (3.8)
We then have,

m Ξ = 0 , [ R Ξ ] = [ I m ] . (3.9)
Finally, Eqs. (3.2) and (3.7) define a mapping h = (h 1 , . . . , h m ) from R n into R m such that Ξ = h(H), and for all ℓ = 1, . . . , ν s and for all j = 1, . . . , m, we have ξ sim,ℓ j = h j (η sim,ℓ ) that can globally be rewritten as ξ sim,ℓ = h(η sim,ℓ ) , ℓ = 1, . . . , ν s .

(3.10)

3.2.

Polynomial chaos expansion with respect to the arbitrary measure defined by the multimodal probability measure of H. In order to define the algebraic notations for the polynomial chaos expansions in a finite approximation, the following parameters and multi-indexes are introduced. Let N = {0, 1, 2, . . .} be the set of all the null and positive integers. For N d ≥ 1 (N d will be the maximum degree of the polynomials), the following set A N d of multi-indexes are introduced,

A N d = {α = (α 1 , . . . , α n ) ∈ N n | 0 ≤ α 1 + . . . + α n ≤ N d } . (3.11)
The 1 + K elements of set A N d are denoted by α (0) , . . . α (K) in which α (0) is the multi-index (0, . . . , 0), and where the integer K is such that

K = (n + N d )! n! N d ! -1 . (3.12)
It is assumed that K < ν s . Let {Ψ α (η) , α ∈ = (α 1 , . . . , α n ) ∈ N n } be the family of multivariate orthonormal polynomials with respect to the multimodal probability measure p H (η) dη on R n defined by Eqs. (2.7) to (2.9), such that, for all α and β in N n ,

R n Ψ α (η) Ψ β (η) p H (η) dη = E{Ψ α (H) Ψ β (H)} = δ αβ . (3.13)
In Eq. (3.13), δ αβ is the kronecker symbol and by convention, for α = 0, Ψ 0 (η) = 1 is the constant normalized polynomial. It should be noted that the probability measure p H (η) dη is non separable (the components of random vector H are not statistically independent). Consequently, the integral on R n defined by Eq. (3.13) cannot be written as a product of n integrals over R, and each multivariate orthonormal polynomial

Ψ α (η) in η = (η 1 , . . . , η n ) cannot be written as Ψ α (η) = ψ α 1 (η 1 ) × . . . × ψ αn (η n )
, but this type of tensorial product of univariate polynomials for the multivariate polynomials (tensorial product that does not hold in the present case) is not used. From a theoretical point of view, the multivariate orthonormal polynomials can be viewed as the result of a Gram-Schmidt orthonormalization algorithm of the multivariate monomials defined, for α = (α 1 , . . . , α n ) belonging to A N d (defined by Eq. (3.11)), and for η = (η 1 , . . . , η n ) in R n , by

M α (η) = η α 1 1 × . . . × η αn n . (3.14)
It should be noted that an explicit algebraic formula (constructed, for instance, with a symbolic toolbox) for each multivariate orthonormal polynomial Ψ α (η), cannot be obtained due to the non separable expression of p H (η) defined by Eqs. (2.7) to (2.9). Nevertheless, even if such an explicit expression could be obtained, this expression would not be used for the reasons detailed in [START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF] and summarized in Section 3.4, and is effectively not used in the methodology proposed for constructing the polynomial chaos expansion of the random vector Ξ (see Sections 3.3 and 3.4). Finally, it can be seen that the degree of multivariate polynomial Ψ α (η) is α 1 + . . . + α n , and that for all α in A N d , the maximum degree of the polynomials is N d .

The second-order random variable Ξ can then be expanded in polynomial chaos Ψ α as,

Ξ = lim N d →+∞ Ξ (N d ) , (3.15)
with convergence in the space L 2 (Θ, R m ) of all the R m -valued second-order random variables defined on (Θ, T , P), and where

Ξ (N d ) = K k=0 y k Ψ α (k) (H) , y 0 , . . . , y K ∈ R m . (3.16)
Taking into account Eq. (3.9), we choose to identify the coefficients y 0 , . . . , y K in R m , such that for each fixed value of N d ,

E{Ξ (N d ) } = 0 , E{Ξ (N d ) (Ξ (N d ) ) T } = [I m ] . (3.17)
From Eqs. (3.13), (3.16), and (3.17), it can be deduced that coefficients y 0 , . . . , y K must satisfy the constraints y 0 = 0 and K k=1 y k (y k ) T = [I m ]. Introducing the rectangular matrix

[ y ] = [y 1 . . . y K ] ∈ M m,K (R) , (3.18)
and the R K -valued random variable Ψ(H) = (Ψ α (1) (H), . . . , Ψ α (K) (H)), then Eq. (3.16) can be rewritten as

Ξ (N d ) = [ y ] Ψ(H) , (3.19)
in which the coefficients

[ y ] verifies [ y ] [ y ] T = [I m ].
Finally, due to the orthonormality of the polynomial chaos (see the second equation in Eq. (3.17)), it can easily be deduced the classical equation,

[ y ] = E{Ξ Ψ(H) T } = E{h(H) Ψ(H) T } , (3.20)
in which Ξ = h(H) is the R m -random variable whose ν s independent realizations ξ sim,1 , . . . , ξ sim,νs have been computed with Eq. (3.7) (see also Eq. (3.10)).

Computing the coefficients of the PCE.

For N d fixed, we propose to compute the coefficients [ y ] of the PCE in using the Monte Carlo method. The use of the ergodic method would lead us to numerically evaluate function h for the M values η k = U k (θ) for k = M 0 , . . . , M of η, while the Monte Carlo method only requires ν s evaluations, and we have seen that M ≫ ν s . In addition, the ν s evaluations of function h have already be carried out for constructing the statistical reduced-order model of random vector Q (see Eqs. (3.2), (3.3) and (3.7)). Therefore, there is a weak interest in developing an approach based on the ergodic method for computing the coefficients of the PCE of random vector Ξ. We then use the ν s independent realizations η sim,1 , . . . , η sim,νs of random vector H, defined by Eq. (3.3), and the ν s independent realizations ξ sim,1 , . . . , ξ sim,νs of random vector Ξ computed in using Eq. (3.7) (i.e. Eq. (3.10)). Taking into account Eqs. (3.8) and (3.9), the estimation [ y N d ] of [ y ] given by Eq. (3.20) is written as

[ y N d ] = 1 ν s -1 [ξ sim ] [Ψ] T . (3.21)
in which the (m × ν s ) real matrix [ξ sim ] and the (K × ν s ) real matrix [Ψ] are defined, for all j = 1, . . . , m, k = 1, . . . , K and ℓ = 1, . . . , ν s , as

[ξ sim ] jℓ = ξ sim,ℓ j , [Ψ] kℓ = Ψ α (k) (η sim,ℓ ) . (3.22)
From Eq. (3.19), random vector Ξ (N d ) , for which [ y ] is estimated by the Monte Carlo method, is written as

Ξ (N d ) = [ y N d ] Ψ(H) . (3.23)
The ν s independent realizations ξ (N d ),1 , . . . , ξ (N d ),νs of random vector Ξ (N d ) are thus given, for all ℓ = 1, . . . , ν s , by ξ

(N d ),ℓ = [ y N d ] Ψ(η sim,ℓ ). Let [ξ N d ]
be the (m × ν s ) real matrix defined, for all j = 1, . . . , m and ℓ = 1, . . . , ν s , by

[ξ N d ] jℓ = ξ (N d ),ℓ j . (3.24)
Consequently, matrix [ξ N d ] can be computed by the following equation,

[ξ N d ] = [ y N d ] [Ψ] . (3.25)
3.4. Computing the realizations of the polynomial chaos. The objective of this section is to compute the (K × ν s ) real matrix [Ψ] defined in Eq. (3.22) and constituted of the ν s independent realizations Ψ α (k) (η sim,ℓ ) of the polynomial chaos Ψ α (k) (H) for the arbitrary multimodal probability measure p H (η) dη defined by Eqs. (2.7) to (2.9). It is assumed that K < ν s . Such a calculation must preserve the orthogonality property of the PCE in high dimension which is written as

lim νs→+∞ 1 ν s -1 [Ψ] [Ψ] T = [I K ] . (3.26)
As it has been proved in [START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF], this problem is not trivial at all. For instance, for classical measures for which the polynomial chaos are explicitly known (Hermite, Legendre, and so on), if the explicit algebraic formulas (constructed with a symbolic Toolbox) or the computational recurrence relations with respect to the degree are used, then important numerical noise is induced, and the orthogonality property is lost. If a global orthogonalization was done to correct this problem, then the independence of the realizations would be lost. In addition, in the present case, we have an arbitrary multimodal measure, and the explicit expression of the polynomial chaos are unknown. We then propose to use the new algorithm proposed in [START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF], that we summarize below. It consists (i) in constructing the ν s independent realizations of the multimodal multivariate monomials using the generator of the multimodal random vector H presented in Section 2.4.4; (ii) in performing an orthogonalization of the realizations of the multimodal multivariate monomials with an algorithm different from the Gram-Schmidt orthogonalization algorithm as the latter is not stable in high dimension. For α = (α 1 , . . . , α n ) belonging to A N d (defined by Eq. (3.11)), and for η = (η 1 , . . . , η n ) in R n , let M α (η) = η α 1 1 × . . . × η αn n be the multivariate monomials. Let K 0 = 1 + K. Let [M] be the (K 0 × ν s ) real matrix of independent realizations of the monomials such that

[M] kℓ = M α (k) (η sim,ℓ ) , k = 0, 1, . . . , K , (3.27)
in which α (0) = (0, . . . , 0), and where η sim,1 , . . . , η sim,νs are the ν s independent realizations of multimodal random vector H, computed with Eq. (3.3). Let [Ψ 0 ] be the (K 0 × ν s ) real matrix such that [Ψ 0 ] kℓ = Ψ α (k) (η sim,ℓ ) for k = 0, 1, . . . , K and ℓ = 1, . . . , ν s . Matrix [Ψ 0 ] can be written as

[Ψ 0 ] = [A] [M] in which [A] is a (K 0 × K 0 ) real matrix, and [M] is a (K 0 × ν s ) real matrix. Let [F] be the positive-definite symmetric (K 0 ×K 0 ) real matrix defined by [F] = 1 νs-1 [M] [M] T = [A] -1 [A] -T .
The algorithm is then the following:

• Compute matrix [M] defined by Eq. (3.27), and then compute the

(K 0 × K 0 ) real matrix [F] = 1 νs-1 [M] [M] T . • Compute the lower triangular (K 0 × K 0 ) real matrix [L] from the Cholesky decomposition [L] [L] T of positive-definite symmetric matrix [F]. • Compute the lower triangular (K 0 × K 0 ) real matrix [A] = [L] -1 . • Compute the (K 0 × ν s ) real matrix [Ψ 0 ] = [A] [M]. • Deduce the (K × ν s ) real matrix [Ψ] obtained in deleting the first row of matrix [Ψ 0 ].

3.5.

Convergence analysis with respect to N d . The numerical convergence with respect to parameter N d is classically analyzed. For n fixed, let K = K(N d ) be defined by Eq. (3.12). Taking into account the convergence in L 2 (Θ, R m ) of the sequence {Ξ (N d ) } N d of random vectors (defined by Eq. (3.19)) towards random vector Ξ, the L 2 -convergence of the polynomial chaos expansion can be performed in analyzing the graph

N d → Ξ -Ξ (N d ) L 2 = (E{ Ξ -Ξ (N d ) 2 }) 1/2 , (3.28)
in which ξ 2 = ξ 2 1 + . . . + ξ 2 n . Introducing the estimation of the second-order moment of real valued random variable Ξ -Ξ (N d ) in using the ν s independent realizations η sim,1 , . . . , η sim,νs of H, yields Ξ -Ξ

(N d ) 2 L 2 ≃ 1 νs-1 [ξ sim ] -[ξ N d ] 2
F in which • F is the Frobenius matrix norm, and where matrices [ξ sim ] and [ξ N d ] are defined by Eqs. (3.22) and (3.25), respectively. For the numerical convergence analysis, we therefore introduce the relative error function

N d → err Ξ (N d ) such that err Ξ (N d ) = [ξ sim ] -[ξ N d ] F [ξ sim ] F . (3.29)
3.6. Numerical application.

3.6.1. Definition of the data transformation. We reuse the numerical application presented in Section 2.6 for which the values of the parameters are N = 100, ν = 500, n = 5, f 0 = 5.0565, ε 0 = 1/200, m overs = 10 that yields M 0 = 260, ∆r = 0.01638, r 0 = 4.26, i 0 = 1. The transformation h from R N in R M , introduced in Section 3 (see Eq. (3.1)), is described in Appendix B. The reader can then simulate the "data transformation" used in this numerical application. Taking into account the convergence analysis presented in Section 2.6.3-(i), and since the number of evaluation of h must not be too large for the practical applications, we choose ν s = 1, 000. In this section, we compare the graphs of the probability density functions relative to random vector Ξ (for which the ν s independent realizations are ξ sim,1 , . . . , ξ sim,νs ), with the graphs of the probability density functions relative to random vector Ξ (N d ) (for which the ν s independent realizations are ξ (N d ),1 , . . . , ξ (N d ),νs ). The computations are performed for ν s = 1, 000, N d = 3 that corresponds to K = 55 with an error err Ξ (N d ) = 0.061. For j in {1, . . . , m}, Fig. 3.4 displays the graphs of the probability density ) and ξ ij = (ξ i , ξ j ) belong to R 2 . Fig. 3.5 shows that there is good agreement between the joint densities.

Conclusions.

A methodology and algorithms have been proposed for constructing the polynomial chaos expansion (PCE) of any R m -valued random variable with respect to an arbitrary multimodal probability measure on R n . The algorithms proposed are efficient and robust, and can be used in high dimension. This new class of algorithms, which is adapted to the multimodal case, can be useful in the context of uncertainty quantification for direct and inverse problems.

Appendix A. Data generation for the numerical application. In this Appendix, the algorithm used for generating the ν independent realizations x exp,1 , . . . , x exp,ν of the random vector X with values in R N , is described. We consider, N = 100, ν = 500, n = 5, n r = 15, c sym = 0.2, and c sep = 2. We introduce a deterministic (n r × n r ) real matrix [a] which is generated by randn( ′ state ′ , 0), and then, [a] = [I nr ] + [randn(n r , n r )]/2, in which we have used the Matlab function, "randn". Then

Figure 2 . 1 .

 21 Figure 2.1. Reduced-order statistical model: graph of the error function n → err(n)

  2 and 2.3 display the graphs of function m ergo → err ER 1 (m ergo M 0 ) and m ergo → err ER 2 (m ergo M 0 ) (with M = m ergo M 0 ). In all the paper, the following parameters are fixed to the values: f 0 = 5.0565, ε 0 = 1/200, m overs = 10 that yields M 0 = 260, ∆r = 0.01638 and r 0 = 4.26.

Figure 2 . 2 .Figure 2 . 3 .

 2223 Figure 2.2. Graph of error function mergo → err ER 1 (mergo M0) for movers = 10 (thick solid line), movers = 100 (mid solid line), movers = 1000 (thin solid line)

Figure 2 . 4 .

 24 Figure 2.4. For i in {1, . . . , n}, graphs of the probability density functions ηi → pH i (ηi) (thick dashed line), and ηi → pH i (ηi) for mergo = 400 (thin solid line) and for mergo = 10, 000 (thick solid line).

Figure 2 . 5 .

 25 Figure 2.5. For (i, j) = (1, 3), (3, 4) and (4, 5), graphs of the joint probability density functions (ηi, ηj ) → pH i H j (ηi, ηj) (left figure), and (ηi, ηj ) → pH i H j (ηi, ηj ) for mergo = 10, 000 (right figure).

3. 6 . 2 .Figure 3 . 1 .Figure 3 . 2 . 3 . 6 . 4 .

 623132364 Figure 3.1. Reduced-order statistical model: graph of the error function m → errQ(m)

Figure 3 . 3 .ξ

 33 Figure 3.3. Graph N d → errΞ(N d ) of the relative error of the polynomial-chaos coefficients versus the maximum degree N d of the polynomial chaos.

pdf Ξ 1 Figure 3 . 5 .

 135 Figure 3.5.For (i, j) = (1, 2),[START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Arnst | Hybrid sampling/spectral method for solving stochastic coupled problems[END_REF] and[START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Arnst | Reduced chaos expansions with random coefficients in reduced-dimensional stochastic modeling of coupled problems[END_REF], graphs of the probability density functions (ξi, ξj) → pΞ i Ξ j (ξi, ξj) and (ξi, ξj) → p Ξ (N d ) i Ξ (N d ) j(ξi, ξj).

  ). Let η → Φ(η) be the function from R n into R such that p H (η) = c n e -Φ(η) .

	(2.14)
	From Eqs. (2.7), (2.8) and (2.14), it can be deduced that

  ii) Estimating E{w(H)} with the Monte Carlo method. Using Eq. (2.31), the Monte Carlo method for estimating E{w(H)} yields

	E{w(H)} = lim νs→+∞	w MC νs	, w MC νs =	1 ν s	νs ℓ=1

  From Eqs. (2.33) and (2.40), the estimation p H of p H is written, for all η in R n , as

	)
	in which η 2 n = η 2 1 + . . . + η 2 n .
	(i) Ergodic method.

ν independent realizations (Y 1ℓ , . . . , Y nrℓ ) (for ℓ = 1, . . . , ν) of a R nr -valued random vector Y = (Y 1 , . . . , Y nr ) are generated such that, randn( ′ state ′ , 0), and then, Y iℓ = nr j=1 [a] ij [randn(n r , ν)] jℓ . The ν independent realizations η exp,1 , . . . , η exp,ν of the R nr -valued random vector η are generated as follows. For i = 1, . . . , n r and for ℓ = 1, . . . , ν, η exp,

] iℓ where z i = 2 i/n r -0.5, and where σ i = 1 + 0.2(2 i/n r -1). Let m = (m 1 , . . . , m nr ) and λ = (λ 1 , . . . , λ nr ) be the vectors in R nr such that m j = 1.5 for all j, and where λ is the vector (1.000, 0.2459, 0.1082, 0.0605, 0.0385, 0.0266, 0.0195, 0.0149, 0.01117, 0.0095, 0.0078, 0.0065, 0.0056, 0.0048, 0.0042). A family of orthonormal vectors {φ 1 , . . . , φ nr } is constructed as the normalized eigenvectors associated with the n r largest eigenvalues of the positive symmetric Appendix B. Data transformation for the numerical application. In this Appendix, the function x → q = h(x) from R N into R M , introduced in Section 3, is defined. It is assumed that M = N . Then, for all α in {1, . . . , N }, q α = z α + { m X } α in which z α = y α /(1 + 6 y 2 α ) with y α = x α -{ m X } α .