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

R Ghanem
email: ghanem@usc.edu

Polynomial chaos representation of databases on manifolds

Keywords: Polynomial chaos expansion, Arbitrary probability measure, Concentration of probability, Measure concentration, Generator, Probability distribution on manifolds, Random sampling generator, MCMC generator, Diffusion maps, Statistics on manifolds

Characterizing the polynomial chaos expansion (PCE) of a vector-valued random variable with probability distribution concentrated on a manifold is a relevant problem in data-driven settings. The probability distribution of such random vectors is multimodal in general, leading to potentially very slow convergence of the PCE. In this paper, we build on a recent development for estimating and sampling from probabilities concentrated on a diffusion manifold. The proposed methodology constructs a PCE of the random vector together with an associated generator that samples from the target probability distribution which is estimated from data concentrated in the neighborhood of the manifold. The method is robust and remains efficient for high dimension and large datasets. The resulting polynomial chaos construction on manifolds permits the adaptation of many uncertainty quantification and statistical tools to emerging questions motivated by data-driven queries.

Introduction

This work is a continuation of a recent paper [1] in which a methodology has been proposed for identifying, from a database made up of N samples of a R ν -valued random variable (possibly for a high value of ν), its non-Gaussian probability distribution that is assumed to be unknown and that is concentrated on an unknown subset S ν of R ν . The method proposes an identification of subset S ν and the construction of generator of samples, based on an Itô stochastic differential equation (ISDE), for the unknown probability distribution that allows for preserving the concentration on S ν and consequently, avoiding the scattering of the generated samples. The main contribution of the present paper is to propose a mathematical formulation and an algorithm for constructing an intrinsic analytical representation based on a polynomial chaos representation (PCE) of the database that is concentrated on subset S ν . The proposed PCE and its associated generator of samples preserves the concentration over subset S ν . In order to obtain an efficient algorithm for constructing the PCE in high dimension, a non classical construction is proposed for which the random germ of the chaos is not constituted of independent normalized Gaussian random variables but depends on the correlation structure of the germ used by the ISDE-based generator of samples. Consequently, the orthonormal multivariate polynomials are constructed with an efficient algorithm for a multivariate probability density function that is not separable with respect to the coordinates.

Constructing the high-dimensional polynomial chaos expansion (PCE) of a random vector H with values in the Euclidean space R ν remains a challenging problem when the probability distribution of the random vector is concentrated on a subset (a manifold) S ν of R ν . In such cases, the probability distribution on R ν is in general multimodal resulting in a potentially very slow convergence of the PCE.

In a recent paper [1], a new methodology was proposed for constructing a generator of samples from a given normalized dataset related to a second-order R ν -valued random variable H = (H 1 , . . . , H ν ) for which the probability density function (pdf), η → p H (η), with respect to the Lebesgue measure dη on R ν is unknown and is concentrated on an unknown subset S ν of R ν . The method consists of first delineating a diffusion manifold for the available data [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF], and then developing a reduced-order Itô stochastic differential equation (ISDE) by the projection on this manifold of an ISDE that admits p H (η) dη as a unique invariant measure. Specifically, we first introduce a random matrix [H] = [H 1 . . . H N ] with values in M ν,N , whose columns are independent copies H 1 , . . . , H N of random vector H and for which the number N of columns is equal to the number of independent data points in an initial dataset. This dataset, suitably normalized, is represented by a matrix [η d ] given in M ν,N , which is taken as a sample of random matrix [H]. The approach presented in [1] is then summarized using the following steps, with an expanded summary provided in section 2 below:

• Multidimensional kernel-density estimation methods [3,[START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] are used to construct the pdf [η] → p [H] ([η]) with respect to the volume element d[η] on M ν,N of random matrix [H], which is assumed to be a second-order random variable. This estimate of the pdf is related to the data and can be used for generating samples of [H]. It should be noted that the method proposed in [1] differs from the MCMC methods on Riemann manifolds that have recently been presented in the very good paper [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF], for which the manifold is the locus of density functions and not of the data itself.

• A generator for random matrix [H] is constructed using the approach 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][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] belonging to the class of Hamiltonian Monte Carlo methods [START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF], which is an MCMC algorithm [START_REF] Kaipio | Statistical and 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]. These samples are obtained by solving an ISDE corresponding to a stochastic nonlinear dissipative Hamiltonian dynamical system, for which p H (η) dη is the unique invariant measure.

• A diffusion-map approach [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF][START_REF] Coifman | Diffusion maps[END_REF][START_REF] Talmon | Intrinsic modeling of stochastic dynamical systems using empirical geometry[END_REF] is then used to discover and characterize the local geometry structure of the normalized dataset concentrated in the neighborhood of S ν . The diffusion-map vectors [g] = [g 1 . . . g m ] ∈ M N,m are thus constructed. They are associated with the first m eigenvalues of the transition matrix relative to the local geometric structure of the given normalized dataset.

• A reduced-order representation [H] = [Z] [g] T is constructed on the manifold in which [Z] is a random matrix with values in M ν,m (with m N ).

• The reduced-order ISDE is obtained by projecting the ISDE from the second step above onto the diffusion manifold by using the reduced-order basis represented by matrix [g] T . The constructed reduced-order ISDE is then used for generating additional samples [z A numerical validation of this procedure was obtained through a number of applications ranging from standard templates to real experimental datasets. These results show the efficiency and the robustness of the method proposed, which allows for concentrating the additional generated samples in the neighborhood of subset S ν .

In the present paper, we propose the construction of a polynomial chaos expansion (PCE) of the second-order random matrix [Z], with the associated generator of samples, which preserves the concentration of the probability measure of the column H 1 , . . . , H N of [H] = [Z] [g] T around S ν . Such a construction will give an intrinsic analytical representation of random matrix [H] through the polynomial chaos expansion of random matrix [Z], which is very useful, for instance, in uncertainty quantification and statistical data analysis, in stochastic modeling and associated statistical inverse problems for boundary value problems. In particular, such a representation could be used for response surface methodology, for computing extreme value statistics, for constructing parameterized reduced-order models, for accelerating robust updating, robust optimization, and robust design.

A systematic construction of PCE of second-order random fields and their use for analyzing the uncertainty propagation in boundary value problems was initiated in [START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A spectral Approach[END_REF]. Since then, PCE and their use in spectral approaches for solving linear and nonlinear stochastic boundary value problems, and some associated statistical inverse problems has rapidly grown [START_REF] Ghanem | Numerical solution of spectral stochastic finite element systems[END_REF][START_REF] Xiu | Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Frauenfelder | Finite elements for elliptic problems with stochastic coefficients[END_REF][START_REF] Debusschere | Numerical challenges in the use of polynomial chaos representations for stochastic processes[END_REF][START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Babuska | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF][START_REF] Doostan | Red-Horse, Stochastic model reduction for chaos representations[END_REF][START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Das | Asymptotic sampling distribution for polynomial chaos representation of data: A maximum-entropy and fisher information approach[END_REF][START_REF] Ghanem | A probability construction of model validation[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Najm | Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics[END_REF][START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Le Maitre | Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics[END_REF][START_REF] Doostan | A non-adapted sparse approximation of PDEs with stochastic inputs[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Nouy | Random fields representations for stochastic elliptic boundary value problems and statistical inverse problems[END_REF]). Several extensions have been proposed concerning the generalized chaos expansions, the PCE for an arbitrary probability measure, and the PCE with random coefficients [START_REF] Lucor | Generalized polynomial chaos and random oscillators[END_REF][START_REF] Soize | Physical systems with random uncertainties: Chaos representation with arbitrary probability measure[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF][START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF][START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF], and recently, the construction of a basis adaptation in homogeneous chaos spaces [START_REF] Tipireddy | Basis adaptation in homogeneous chaos spaces[END_REF].

Significant work has also been devoted to the acceleration of stochastic convergence of the PCE [START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF][START_REF] Keshavarzzadeh | Convergence acceleration of polynomial chaos solutions via sequence transformation[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF][START_REF] Tipireddy | Basis adaptation in homogeneous chaos spaces[END_REF][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF]. Nevertheless, the construction of a PCE of a second-order random vector for which the probability measure is concentrated on a subset S ν of R ν , remains a difficult and challenging problem.

In [1], it has been shown that a direct sampling of H using, for instance, a nonparametric estimation of the pdf of H constructed with the given normalized dataset and an MCMC method for generating samples, yields samples that are not concentrated around subset S ν , but that are scattered. This was indeed the motivation for the method proposed in [1] using the diffusion maps and summarized above. In order to preserve the concentration of the probability measure for the analytical representation of the second-order random matrix [H], constraining associated samples of its columns H 1 , . . . , H N to be concentrated in the neighborhood of S ν , the following methodology is proposed.

• The PCE of random M ν,m -valued random variable [Z] is performed with respect to a M Ng,µ -valued random germ [Ξ], with 1 ≤ N g ≤ ν and 1 ≤ µ ≤ m. The random germ [Ξ]
is constructed as a linear mapping of the matrix-valued Wiener process that is used in the reduced-order ISDE that allows the samples of random matrix [Z] to be generated on the manifold. • The R ν -valued random variable H is an application-specific normalization of a R n -valued random variable X (with ν ≤ n) for which N samples constitute the given dataset. According to our formulation, these N samples of X constitute one sample of random matrix [X] with values in M n,N .

Remark on the methodology proposed. Another approach could be envisaged for generating additional samples without recourse to the reduced-order ISDE. Such an approach would consist in performing a direct generation of samples of random matrix [H] from the multidimensional Gaussian kernel estimation method for which the concentration would be lost due to unavoidable scattering of the generated samples. These scattered samples would then be projected on the subspace spanned by the m N diffusion-maps vectors represented by matrix [g] in order to generate the corresponding samples of [Z]. With such an approach, the stochastic germ, which would be used by this Gaussian kernel generator, would not live on the "manifold" identified by the diffusion maps. The corresponding generator would produce samples that would belong to the big set M ν,N , before projecting them on the identified "manifold". Such an approach has not been investigated in this work for the following theoretical reason. The proposed method has been developed with two constraints in mind. The first one concerns the stochastic germ of the reduced-order ISDE that is used for generating the samples of [Z], which must live on the "manifold" identified by the diffusion maps, that is to say, which must live on a subset of the set M ν,m , which has a small dimension because m N . The second criterion requires the samples of [Z] to be directly generated by the reduced-order ISDE on the "manifold" which is a subset of M ν,m with a small dimension. The proposed generator based on a reduced-order ISDE, which allows these two constraints to be easily implemented. These two constraints are significant since they allows for constructing the PCE on the "manifold" without ambiguity and with a very efficient algorithm, as proposed in the manuscript.

Organization of the paper. The paper is organized as follows. Section 2 is devoted to a brief review of the methodology and algorithm proposed in [1] for generating samples that follow the pdf of the dataset, in the neighborhood of S ν . Section 3 deals with the PCE of random matrix [Z], which capitalizes on the results summarized in Section 2. In Section 4, using the PCE of random matrix [Z], we present the generation of additional samples of random matrix [X] for which the dataset constitutes one sample. Finally, Section 5 is devoted to numerical examples.

2. Short summary of the methodology and algorithm proposed in [1] for a concentrated probability measure

In this section, we give a brief summary of the methodology presented in [1], which underlies the proposed PCE construction. Let X = (X 1 , . . . , X n ) be a random vector, defined on a probability space (Θ, T , P), with values in R n , for which the probability density function (pdf) p X on R n is unknown but is concentrated on an unknown subset S n of R n , and for which N data points x d,1 , . . . , x d,N in R n correspond to N statistically independent samples of X, which are given and which are represented by the matrix

[x d ] = [x d,1 . . . x d,N ] in M n,N . Only using [x d ]
as available information, this method allows for generating samples of random vector X that are concentrated in the neighborhood of S n . The steps detailed in [1] are summarized hereinafter.

Scaling and normalizing the given dataset [x d ]

From a given unscaled dataset represented by a matrix

[x uns d ] in M n,N , the (scaled) dataset [x d ] in M n,N is easily deduced from [x uns d ].
The normalization of X is obtained by a principal component analysis. Introducing the random matrix [X] = [X 1 , . . . , X N ] with values in M n,N , whose columns are N independent copies of random vector X, we introduce the corresponding normalized random matrix [H] = [H 1 , . . . , H N ] with values in M ν,N , whose columns are N independent copies of random vector H, with ν ≤ n, such that

[X] = [x] + [ϕ] [λ] 1/2 [H] , (1)
in which [λ] is the (ν × ν) diagonal matrix of the ν positive eigenvalues of the empirical estimation of the covariance matrix of X, where [ϕ] is the (n×ν) matrix of the associated eigenvectors such

[ϕ] T [ϕ] = [I ν ],
and where [x] is the matrix in M n,N for which each column is the empirical estimation of the mean value of random vector X. The sample

[η d ] = [η d,1 . . . η d,N ] ∈ M ν,N of [H] (associated with the sample [x d ] of [X]) is thus computed by [η d ] = [λ] -1/2 [ϕ] T ([x d ] -[x]) . (2)
Consequently, the empirical estimates of the mean value and of the covariance matrix of random vector H are 0 ν and [I ν ], respectively.

Constructing the diffusion-maps vectors for H

Let k ε (η, η ) = exp(-1 4ε η -η 2
) be the kernel defined on R ν × R ν , depending on a real smoothing parameter ε > 0. It should be noted that this kernel could be replaced by another one satisfying the symmetry, the positivity preserving, and the positive semi-definiteness properties. For m ≤ N , let [g] = [g 1 . . . g m ] ∈ M N,m be the "diffusion-maps basis" associated with kernel k ε , which is defined and constructed in Appendix A (for m = N , [g] is an algebraic basis of R N ). For α = 1, . . . , m, the diffusion-maps vector g α ∈ R N is defined by Eq. (A.3). The subspace of R N spanned by the vector basis {g α } α allows for characterizing the local geometry structure of the dataset concentrated in the neighborhood of S ν .

Introducing the reduced-order representation of random matrix [H]

The reduced-order representation is obtained in projecting each column of the M N,ν -valued random matrix [H] T on the subspace of R N , spanned by {g 1 . . . g m }.

Introducing the random matrix [Z] with values in M ν,m , the following reducedorder representation of [H] is introduced,

[H] = [Z] [g] T .
( 3) yields the least squares approximation to Z in the form,

) Since the matrix [g] T [g] ∈ M m is invertible, Eq. ( 3
[Z] = [H] [a] , [a] = [g] ([g] T [g]) -1 ∈ M N,m . (4)
In particular, matrix

[η d ] ∈ M ν,N can be written as [η d ] = [z d ] [g] T in which the matrix [z d ] ∈ M ν,m is given by [z d ] = [η d ] [a] ∈ M ν,m . (5)
Consequently, from Eqs. ( 1) and ( 4), it can be deduced the following representation of random matrix [X] as a function of random matrix [Z],

[X] = [x] + [ϕ] [λ] 1/2 [Z] [g] T . (6)

Estimating dimension m of the reduced-order representation of random matrix [Z]

For a given value of integer ζ related to the analysis scale of the local geometric structure of the dataset (see Eq. (A.3) in Appendix A) and for a given value of smoothing parameter ε > 0, the decreasing of the graph α → Λ α of the positive eigenvalues of transition matrix [P] (see Appendix A) yields a criterion for choosing the value of m that allows the local geometric structure of the dataset represented by [η d ] to be discovered. Nevertheless, this criterion may not be sufficient, and the L 2 -convergence may need to be enforced by increasing, as required, the value of m. However, if the value of m is chosen too large, the localization of the geometric structure of the dataset is lost. Consequently, a compromise must be reached between the very small value of m given by the decreasing criteria of the eigenvalues of matrix [P] ∈ M N and a larger value of m which is necessary for obtaining a reasonable mean-square convergence. A criterion for estimating an optimal value of m is given in Appendix B.

Reduced-order ISDE for generating additional samples

{([Z(r)], [Y(r)]), r ∈ R + }, defined on (Θ, T , P), indexed by R + = [0 , +∞[, with values in M ν,m ×M ν,m
, which is the unique second-order stationary (for the shift semi-group on R + ) and ergodic diffusion stochastic process, of the following reduced-order ISDE, for r > 0,

d[Z(r)] = [Y(r)] dr , (7)
d[Y(r)] = [L([Z(r)])] dr - 1 2 f 0 [Y(r)] dr + f 0 [dW(r)] , (8)
with the initial condition

[Z(0)] = [H d ] [a] , [Y(0)] = [N ] [a] a.s , (9)
in which the random matrices

[L([Z(r)])] and [dW(r)] with values in M ν,m are such that [L([Z(r)])] = [L([Z(r)] [g] T )] [a] , [dW(r)] = [dW(r)] [a] . ( 10
) (i) For all [u] = [u 1 . . . u N ] in M ν,N with u = (u 1 , . . . , u ν ) in R ν , the matrix [L([u])] in M ν,N
is defined, for all k = 1, . . . , ν and for all = 1, . . . , N , by

[L([u])] k = 1 q(u ) {∇ u q(u )} k , (11)
q(u ) = 1 N N j=1 exp{- 1 2 s 2 ν s ν s ν η d,j -u 2 } , (12)
∇ u q(u ) = 1 s 2 ν 1 N N j=1 ( s ν s ν η d,j -u ) exp{- 1 2 s 2 ν s ν s ν η d,j -u 2 } , (13)
s ν = 4 N (2 + ν) 1/(ν+4) , s ν = s ν s 2 ν + N -1 N . ( 14
) (ii) The stochastic process {[dW(r)], r ≥ 0} with values in M ν,N is such that [dW(r)] = [dW 1 (r) . . . dW N (r)] in which the columns W 1 , . . . , W N are N in- dependent copies of the normalized Wiener process W = (W 1 , . . . , W ν ) de- fined on (Θ, T , P), indexed by R + with values in R ν . The matrix-valued au- tocorrelation 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 ν ].
(iii) The probability distribution of the random matrix (iv) The free parameter f 0 > 0 allows the dissipation term of the nonlinear secondorder dynamical system (dissipative Hamiltonian system) to be controlled.

[H d ] with values in M ν,N is p [H] ([η]) d[η]. A known sample of [H d ] is matrix [η d ]
(v) The algorithm for solving Eqs. ( 7) to ( 9) is detailed in [1] and is summarized in Appendix C. 7) to ( 9)). For ρ = M 0 ∆r, in which ∆r is the sampling step of the continuous index parameter r used in the integration scheme (see Appendix C), and where M 0 is a positive integer greater or equal to 1 such that M = M 0 × n MC , the additional samples [z

Generating additional samples

[x 1 ar ], . . . , [x nMC ar ] of random matrix [X] For θ fixed in Θ, the deterministic quantities {[W(r; θ)], r ≥ 0}, [Z(0; θ)] = [η d ] [a]
[z ar ] = [Z( ×ρ ; θ)] , [η ar ] = [z ar ] [g] T , = 1, . . . , n MC . (15)
• If M 0 = 1, then ρ = ∆r and the n MC additional samples are dependent, but the ergodic property of {[Z(r)], r ∈ R + } can be used for obtaining the convergence of statistics constructed using [z

m k = 1 N×n MC N j=1 nMC =1 [η ar ] kj , (16)
[ c ] kk = 1 N×n MC -1 N j=1 nMC =1 ([η ar ] kj -m k )([η ar ] k j -m k ) . (17)
[η ar ] kj = ν k =1 {[ L] -T } kk ([η ar ] k j -m k ) . (18)
Consequently, the empirical estimates of the mean value and covariance matrix performed with the re-normalized samples {[η ar ], = 1, . . . n MC } yield the zero vector and the identity matrix for random vector H, respectively.

• Using Eqs. (4), it is deduced that the samples [z 18), can be computed by

[z ar ] = [η ar ] [a] , = 1, . . . n MC . (19)
Using Eq. ( 6), the additional samples [x 1 ar ], . . . , [x nMC ar ] of random matrix [X] can be generated by

[x ar ] = [x] + [ϕ] [λ] 1/2 [z ar ] [g] T , = 1, . . . n MC . (20)

Polynomial chaos expansion of random matrix [Z]

3.1. Construction of an adapted Gaussian germ for the PCE on the manifold In general, N g ≤ ν and µ ≤ m. However, for explaining the construction of the Gaussian germ [Ξ] with values in M Ng,µ , we construct the germ [Ξ(ν, m)] with values in M ν,m and then, for N g < ν and µ < m, the germ [Ξ] with values in M Ng,µ will be obtained by extracting the block

(N g × µ) from [Ξ(ν, m)].
The idea is to construct a Gaussian random matrix [Ξ(ν, m)], defined on (Θ, T , P) with values in M ν,m that is statistically dependent of the M ν,m -valued stochastic process {[W(r)], r ≥ 0}, in order to perform the computation of the coefficients of the PCE of random matrix [Z] by a direct use of the projection formula. Nevertheless, in order to avoid potential numerical difficulties for computing the samples of the PCE with respect to [Ξ], we need to introduce a germ whose statistical fluctuations have unit variance. Consequently, we will perform the construction by normalizing {[W(r)], r ≥ 0} for which the statistical fluctuations can have a variance greater than 1 due to the presence of matrix

Taking into account Eq. (C.3) of Appendix C, it can easily be seen that, for all fixed in {1, . . . , M }, the fourth-order covariance tensor of random matrix [∆W ], is independent of , and is such that, for all α and α in {1, . . . , m}, and for all k and k in {1, . . . , ν},

E{ [∆W ] kα [∆W ] k α } = δ kk [A] αα , (22)
in which the positive-definite matrix [A] is written as

[A] = ∆r [a] T [a] ∈ M + m . (23)
We then defined the random matrix [Ξ(ν, m)] as the Gaussian, centered, random matrix with values in M ν,m , for which the fourth-order covariance tensor is such that, for all k and k in {1, . . . , ν} and for all α and α in {1, . . . , m},

E{ [Ξ(ν, m)] kα [Ξ(ν, m)] k α } = δ kk [A] αα [A] αα [A] α α , (24)
in which the matrix [A] is defined by Eq. [START_REF] Ganapathysubramanian | Sparse grid collocation schemes for stochastic natural convection problems[END_REF]. By comparing Eq. ( 22) with Eq. ( 24), and in taking ρ and the sample θ in Θ used in Eq. ( 15), it can be concluded that n MC samples, {[ξ (ν, m)], = 1, . . . , n MC } of random matrix [Ξ(ν, m)] can be constructed, for = 1, . . . , n MC , as

[ξ (ν, m)] = [∆W ×ρ (θ)] [d ] , (25)
in which [d ] is the (m × m) real diagonal positive-definite matrix which is written as,

[d ] αα = δ αα [A] αα . ( 26
)
For fixed N g ≤ ν and µ ≤ m, the n MC samples, {[ξ ], = 1, . . . , n MC } of random matrix [Ξ] are then given by extracting, for = 1, . . . , n MC , the block

(N g × µ) from [ξ (ν, m)], [ξ ] kα = [ξ (ν, m)] kα , k = 1, . . . , N g , α = 1, . . . , µ . (27) 3.2
B N d ,Ng,µ = {β ∈ N µNg | 0 ≤ |β| ≤ N d } . (28)
The K elements of set B N d ,Ng,µ , which depend on N d , N d , and µ (and also written as K(N d , N g , µ)), are denoted by β (1) , . . . β (K) in which β (1) is the multi-index (0, . . . , 0), and where the integer

K(N d , N g , µ) is such that K(N d , N g , µ) = (µN g + N d )! (µN g )! N d ! . ( 29
)
Orthonormal multivariate polynomials. Let p [Ξ] ([ξ]) d[ξ] be the Gaussian, centered, probability measure of the random matrix [Ξ] with values in M Ng,µ , for which the fourth-order covariance tensor, which is derived from Eq. ( 24), is written, for all k and k in {1, . . . , N g }, and for all α and α in {1, . . . , µ}, as

E{ [Ξ] kα [Ξ] k α } = δ kk [A] αα [A] αα [A] α α . ( 30
)
Note that the probability measure

p [Ξ] ([ξ]) d[ξ] is such that p [Ξ] ([ξ]) d[ξ] = p Ξ 1 ,...,Ξ µ (ξ 1 , . . . , ξ µ ) dξ 1 . . . dξ µ . ( 31
)
in which ξ α ∈ R Ng and where dξ α is the Lebesgue measure on R Ng . For all multiindices β ∈ N µNg , we introduce the real-valued multivariate polynomials Ψ β ([ξ]) on M Ng,µ , which is defined by

Ψ β ([ξ]) = Ψ β 1 ,...,β µ (ξ 1 , . . . , ξ µ ) , (32)
in which

[ξ] = [ξ 1 . . . ξ µ ]. Let {Ψ β ([ξ]), β ∈ N µNg }
Ψ β ([ξ]) Ψ β ([ξ]) p [Ξ] ([ξ]) d[ξ] = E{Ψ β ([Ξ]) Ψ β ([Ξ])} = δ ββ , ( 33
)
where δ ββ is the Kronecker symbol. By convention, for β = β (1) = (0, . . . , 0), Ψ β (1) ([ξ]) = 1 is the constant normalized multivariate polynomial. Since the random vectors Ξ 1 , . . . , Ξ µ are mutually dependent, the probability density function p Ξ 1 ,...,Ξ µ (ξ 1 , . . . , ξ µ ) is not separable in ξ 1 , . . . , ξ µ and consequently, the orthonormal multivariate polynomials Ψ β ([ξ]) cannot be written as

Ψ β 1 (ξ 1 ) × . . . × Ψ β µ (ξ µ ).
From a theoretical point of view, the orthonormal multivariate polynomials can be viewed as the result of a Gram-Schmidt orthonormalization algorithm of the multivariate monomials (see Section 3.5) defined, for

β ∈ N µNg , by M β ([ξ]) = M β 1 (ξ 1 ) × . . . × M β µ (ξ µ ) , (34)
in which for all α in {1, . . . , µ},

M β α (ξ α ) = (ξ α 1 ) β α 1 × . . . × (ξ α Ng ) β α Ng . (35)

Polynomial chaos expansion of random matrix [Z]

In this section, we construct the PCE of random matrix [Z]. Since [Z] is a second-order random variable with values in M ν,m , random matrix [Z] can be expanded in polynomial chaos by using the classical theory. What is different in our case, it is the fact that we propose to construct the truncated PCE of random matrix [Z] with respect to the germ [Ξ] defined in Section 3.2, which is related to the correlation structure of the germ on the manifold. We then have to analyze the problem related to the L 2 convergence of the truncated PCE of

[Z(N d , N g , µ)] = K(N d ,Ng,µ) κ=1 [y (κ) ] Ψ β (κ) ([Ξ]) , [y (κ) ] ∈ M ν,m , (36)
in which the integer K(N d , N g , µ) and the multi-indices β (κ) are defined in Section 3.2 and where Ψ β ([ξ]) is defined by Eq. ( 32). The family {[y (κ) ]} κ of the PCE coefficients are matrices in M ν,m , which are computed by using the formula,

[y (κ) ] = E{[Z] Ψ β (κ) ([Ξ])} . ( 37
)
The L 2 -convergence,

[Z] = lim µ→m Ng →ν N d →+∞ [Z(N d , N g , µ)] , (38)
is analyzed in paragraph (ii) below.

It should be noted that the statistical dependence between random matrix [Ξ] and random matrix [Z] allows the development of very efficient algorithms for computing the coefficients [y (κ) ] of the PCE of [Z] using Eq. [START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF]. Indeed, if [Ξ] had been chosen as an independent random matrix of [Z], then Eq. ( 37) would give [y (κ) ] = [0] for all κ, because E{[Z]} = [0], and consequently, could not be used. The maximum likelihood method should then be used for identifying the K(N d , N g , µ) matrix-valued coefficients. This would entail an optimization problem with ν × µ × K(N d , N g , µ) variables, which would quickly become computationally prohibitive.

Z] = f (ν,m) ([Ξ(ν, m)]) in which [ξ] → f (ν,m) ([ξ]) is a measurable mapping from M ν,m into M ν,m . As [Z] is a second-order random variable, we have Mν,m f (ν,m) ([ξ]) 2 F p [Ξ] ([ξ]) d[ξ] < +∞ ,
in which • F is the Frobenius norm and where p [Ξ] is defined by Eq. ( 31). On the other hand, it can easily be seen that, for any multi-index β (κ) with κ in N, the multivariate monomial M β (κ) ([ξ]) defined by Eq. ( 34) is such that

Mν,m |M β (κ) ([ξ])| p [Ξ] ([ξ]) d[ξ] < +∞ .
Consequently, the mathematical results proved in [START_REF] Soize | Physical systems with random uncertainties: Chaos representation with arbitrary probability measure[END_REF] allow us to conclude that the sequence

{[Z(N d , ν, m)]} N d is convergent in L 2 to [Z] = f (ν,m) ([Ξ(ν, m)]) when N d → +∞, that is to say, [Z] = lim N d →+∞ [Z(N d , ν, m)] . (39)
(ii-2) The correlation structure that we have retained (see Eq. ( 30) deduced from Eq. ( 24)) for [Ξ] shows that the components (Ξ α 1 , . . . , Ξ α Ng ) of each column Ξ α with values in R Ng of random matrix [Ξ] are mutually independent (Gaussian, centered, and not correlated), but the columns Ξ 1 . . . Ξ µ are mutually dependent Gaussian vectors. Consequently, for µ = m, the truncated PCE

[Z(N d , N g , m)] of [Z] with respect to [Ξ] yields a corresponding truncated PCE Z α (N d , N g , m) of each column Z α of [Z] that depends on [Ξ] = [Ξ 1 . . . Ξ m ]. In particular, the truncated PCE Z α (N d , N g , m) of column Z α with values in R ν depends on Ξ α
with values in R Ng , for which all the components (Ξ α 1 , . . . , Ξ α Ng ) are independent (because Z α depends on Ξ 1 , . . . , Ξ µ and consequently, depends on Ξ α ). This property and Eq. ( 39) guarantee that such a truncated PCE Z

α (N d , N g , m) of Z α with values in R ν is convergent in L 2 when N g → ν and N d → +∞. Conse- quently, if [Z(N d , N g , µ)
] denotes the truncated PCE of random matrix [Z], it can be deduced that, for µ = m, we have the following L 2 -convergence,

[Z] = lim Ng →ν N d →+∞ [Z(N d , N g , m)] . (40)
(ii-3) However, for µ fixed such that µ < m, due to the dependence of the Gaussian random vectors Ξ 1 , . . . , Ξ µ , Eqs. ( 39) and ( 40) do not prove the convergence in

L 2 of the truncated PCE [Z(N d , N g , µ)] of [Z] when N g → ν and N d → +∞.
As the explicit expression of the probability distribution of random matrix [Z] is unknown, it seems difficult to prove such a convergence. In this case we settle for a numerical exploration using an adapted criterion such that the one defined by Eq. ( 50)) below. Nevertheless, due to Eq. ( 40), for µ = m, the given theoretical proof ensures convergence.

Estimation of the matrix-valued coefficients

In this section, integers n MC , N g , µ, and N d are fixed. The construction of the truncated PCE [Z(N d , N g , µ)] of random matrix [Z] requires two steps. The first one is the numerical calculation of the samples of the polynomial chaos. The second one is the computation of the PCE coefficients of random matrix [Z]. These steps are detailed next and for simplifying the notation, K(N d , N g , µ) is sometimes simply written as K.

(i) Samples of the polynomial chaos for the computation of the PCE coefficients.

Let Ψ([Ξ]) = (Ψ β (1) ([Ξ]), . . . , Ψ β (K) ([Ξ]
)) be the random vector with values in R K . The n MC samples of the K multivariate polynomials are represented by the matrix [Ψ] ∈ M K,nMC such that

[Ψ] κ = Ψ β (κ) ([ξ ]) , κ = 1, . . . , K , = 1, . . . , n MC , (41)
where [ξ ] is defined by Eq. ( 27) (because we are computing the coefficients {[y (κ) ]} κ ). Due to Eq. ( 33), matrix [Ψ] must be such that

lim nMC→+∞ 1 (n MC -1) [Ψ] [Ψ] T = [I K ] . ( 42
)
(ii) Computation of the coefficients. From Eqs. ( 37), [START_REF] Ghosh | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF], and ( 42), it can be deduced that, for n MC sufficiently large (greater than K), the entries [y (κ) ] kα of an estimation of matrix-valued coefficient [y [κ) ], can be written, for κ = 1, . . . , K, for k = 1, . . . , ν, and for α = 1, . . . , m, as

[y (κ) ] kα = [ y α ] kκ , (43)
[ y α ] 1 (n MC -1) [ z α ar ] [Ψ] T . ( 44
)
In Eq. ( 44), the matrix [ z α ar ] ∈ M ν,nMC is defined by

[ z α ar ] k = [z ar ] kα , = 1, . . . , n MC , (45)
in which [z 1 ar ], . . . , [z nMC ar ] are the samples of random matrix [Z], computed with Eq. ( 19).

Computation of matrix [Ψ]

For the reasons detailed in [START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF], matrix [Ψ] is computed with the algorithm based on the Cholesky factorization proposed in [START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF] (we refer the reader to [START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF] for a method involving the singular value decomposition and to [START_REF] Arnst | Hybrid sampling/spectral method for solving stochastic coupled problems[END_REF] for a method involving the QR factorization). This algorithm preserves the orthogonality property defined by Eq. ( 42) for the high dimensions and requires that n MC > K.

Let [M] be the (K × n MC ) real matrix of samples of the monomials defined by Eq. ( 34) such that

[M] κ = M β (κ) ([ξ ]) , κ = 1 . . . , K , = 1, . . . , n MC , (46)
in which [ξ ] is defined by Eq. ( 27) and where M β (κ) ([ξ ]) is computed by using Eqs. [START_REF] Lucor | Generalized polynomial chaos and random oscillators[END_REF] and [START_REF] Soize | Physical systems with random uncertainties: Chaos representation with arbitrary probability measure[END_REF]. The algorithm is then as follows:

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

(K ×K) real matrix [F] = 1 nMC-1 [M] [M] T .
Since n MC > K, then matrix [F] to be positive definite.

• Compute the lower triangular

(K × K) real matrix [L] from the Cholesky decomposition [L] [L] T of positive-definite symmetric matrix [F].
• Compute the (K × n MC ) real matrix [Ψ] as the solution of the linear matrix equation

[L] [Ψ] = [M].
It should be noted that, when the rank of matrix [M] is greater or equal to K(N d , N g , µ), which corresponds to n MC > K, the proposed algorithm shows that the constructed polynomials are exactly orthonormal.

Computation of n MC samples of [Z] using its PCE

Let N max d ≥ 1 be the greatest value of N d , which is considered for performing the L 2 -convergence analysis of the PCE of [Z]. We thus consider N d , N g , and

µ fixed such that 1 ≤ N d ≤ N max d , 1 ≤ N g ≤ ν, and 1 ≤ µ ≤ m (see Sec- tion 3.3-(i)). We are thus interested in computing n MC samples [z 1 PCE ], . . . , [z nMC PCE ] of the random matrix [Z(N d , N g , µ)] such that, for = 1, . . . , n MC , [z PCE ] kα = [ z α PCE ] k , k = 1, . . . , ν , α = 1, . . . , m , (47)
in which, for α = 1, . . . , m, the matrix

[ z α PCE ] ∈ M ν,nMC is computed by [ z α PCE ] = [ y α ] [Ψ] . (48)
The matrix [ y α ] ∈ M m,K is given by

[ y α ] kκ = [y (κ) ] kα , k = 1, . . . , ν , κ = 1, . . . , K , (49)
in which matrix [y (κ) ] ∈ M ν,m has been computed by using Eqs. ( 43) to ( 45) and where The matrix [Ψ] ∈ M K,nMC has been computed by using Eq. ( 41). 38)) is performed by constructing the error function

(N d , N g , µ) → error [Z] (N d , N g , µ) defined on [1, N max d ]×[1, ν]× [1, m] by error [Z] (N d , N g , µ) = nMC =1 [z ar ] -[z PCE ] 2 F nMC =1 [z ar ] 2 F . ( 50
)
4. Generating additional samples of random vector X using its analytical representation and L 2 -error

For N d , N g , and µ fixed to the values identified by the analysis of the error function (N d , N g , µ) → error [Z] (N d , N g , µ) defined by Eq. ( 50), the number K = K(N d , N g , µ) of coefficients is known. For a given number n sim of additional samples, the methodology consists of • computing n sim samples of the analytical representation X(N d , N g , µ) of random vector X by using Eq. ( 6) and the PCE [Z(N d , N g , µ)] of random matrix [Z];

• estimating the L 2 -error for X(N d , N g , µ) with respect to random vector X.

4.2. Estimation of the L 2 -error for X(N d , N g , µ) with respect to random vector X For given N d , N g , µ, and n sim , the L 2 -error of the random vector X(N d , N g , µ) with respect to random vector X is estimated by

error X (N d , N g , µ) = [cov X(N d ,Ng,µ) ] -[cov X ] F [cov X ] F , (51)
in which the covariance matrix [cov X(N d ,Ng,µ) ] ∈ M n of the random vector X(N d , N g , µ) is estimated with the n sim × N samples x 1 PCE,sim , . . . , x nsim×N PCE,sim , and where the covariance matrix [cov X ] ∈ M n of the random vector X is estimated with the n MC × N samples x 1 ar , . . . , x nMC×N ar that are the columns of the matrices [x 1 ar ], . . . [x nMC ar ], that are computed using Section 2 (see Eq. ( 20)). It should be noted that the number of samples for random vectors X and X(N d , N g , µ) are not the same and, in addition, X and X(N d , N g , µ) are statistically independent. Consequently, the classical L 2 -norm cannot be introduced for quantifying the convergence, and since X and X(N d , N g , µ) have the same mean vectors, the distance between the covariance matrices is used.

Applications

Three applications presented in [1] are partially reused for illustrating the methodology proposed. The three examples consist of obtaining samples of random vector X with values in R n such that:

1. in example 1, the dimension is n = 2 and there are N = 230 given data points in subset S n , for which the mean value is made up of two circles in the plane). 2. in example 2, the dimension is n = 3 and there are N = 400 given data points in subset S n , for which the mean value is made up of a helix in threedimensional space). 3. the third example corresponds to a petro-physics database that is made up of experimental measurements (downloaded from [46]) and detailed in [START_REF] Thimmisetty | Multiscale stochastic representation in high-dimensional data using Gaussian processes with implicit diffusion metrics[END_REF], for which the dimension is n = 35 and for which N = 13, 056 given data points are concentrated in an unknown "complex" subset S n of R n , which cannot be easily described once discovered. For this first application, the statistical fluctuations are around two circles. The number of given data points is N = 230, no scaling of data is performed, but the normalization defined in Section 2.1 is done. Fig. 1 (left) displays the 230 given data points for random vector X = (X 1 , X 2 ) of the dataset represented by matrix [x d ] in M 2,230 , and shows that the given data points are concentrated in the neighborhood of two circles. The kernel is defined in Section 2.2, the value used for the smoothing parameter is ε = 100 × 2π/N = 2.73, parameter ζ defined in Appendix A is chosen to 1, and the graph of the eigenvalues of the transition matrix for random vector H is displayed in Fig. 1 (right). This graph shows that dimension m can be chosen to be 3, for which the value of e red (m) (defined by Eq. (B.4)) is 6.96 × 10 -4 (it can thus be considered that a reasonable mean-square convergence has been reached). Fig. 2 (left) displays the pdf for random variables X 1 and X 2 computed with a nonparametric estimation from the data points. For generating 9, 200 additional samples, the numerical values of the parameters are f 0 = 1.5, ∆r = 0.1179, M 0 = 110, and n MC = 40, yielding 4, 400 for the parameter M defined in Appendix C. The additional samples [x 1 ar ], . . . , [x nMC ar ] (with n MC = 40) are computed by using Section 2.6 and are displayed in Fig. 2 (right), which shows the 230 given data points and the 9, 200 additional samples generated using the reduced-order ISDE. It can be seen that the additional samples are effectively concentrated in subset S n (as it has be shown in [1], if a direct simulation was used instead of the method proposed, then the 9, 200 additional samples would not be concentrated in subset S n , but would be scattered). The number of given data points is N = 400, no scaling of data is performed, but the normalization defined in Section 2.1 is done. Fig. 6 (left) displays the 400 given data points for random vector X = (X 1 , X 2 , X 3 ) of the dataset represented by matrix [x d ] in M 3,400 , and shows that the given data points are concentrated in the neighborhood of a helical. The kernel is defined in Section 2.2, the value of the smoothing parameter that is retained is ε = 100×2π/N = 1.57, parameter ζ defined in Appendix A is set to 1, and the graph of the eigenvalues of the transition matrix for random vector H is displayed in Fig. 6 (right). This graph shows that dimension m can be chosen to be 4 for which the value of e red (m) (defined by Eq. (B.4)) is 5.48 × 10 -4 (it can thus be considered that a reasonable mean-square convergence has been reached). Fig. 7 (left) displays the pdf for random variables 7 (right), which shows the 400 given data points and the 8, 000 additional samples generated using the reduced-order ISDE. It can be seen that the additional samples are effectively concentrated in subset S n (as it has be shown in [1], if a direct simulation was used instead of the method proposed, then the 8, 000 additional samples would not be concentrated in subset S n , but would be scattered).

The analytical representation of [X] is constructed using the methodology presented in Sections 3 and 4. All the computations are performed with f 0 = 1.5, ∆r = 0.1197, M 0 = 110, and n MC = 2, 000. The value of µ is fixed such that µ = m = 4. ples generated using the analytical representation. It can be seen that the samples are effectively concentrated in subset S n .

Application 3: Dimension n = 35 with N = 13, 056 given data points

The database used corresponds to petro-physics field observations. The dimension of random vector X is n = 35 and the number of given data points is N = 13, 056. The scaling defined in Section 2.1 is necessary and has been done. The scaled data points are then represented by the matrix [x d ] in M 35,13056 . The normalization is performed using Eq. (1) for which ν = 32. Figure 11 displays the 13, 056 given data points viewed from coordinates x 16 and x 28 , from coordinates x 27 and x 28 , and from coordinates x 30 , x 32 , and x 33 . Although only a partial representation of the data points are shown in these three figures, it can be seen that S n is certainly a complex subset of R n . The kernel is defined in Section 2.2, the value of the The analytical representation is constructed using the methodology presented in Section 3. All the computations are performed with f 0 = 1.5, ∆r = 0.06142, M 0 = 1, and n MC = 25, 000 (the ergodic property is used for estimating the coefficients of the PCE). 51)) is error X (N d , N g , µ) = 6.6 × 10 -3 . The same components X 16 , X 27 , X 28 , X 30 , X 32 , and X 33 that those used in Fig. [START_REF] Talmon | Intrinsic modeling of stochastic dynamical systems using empirical geometry[END_REF] The algorithm for solving the reduced-order ISDE defined by Eqs. [START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF] to ( 9) is detailed in [1] and is summarized hereinafter. The Störmer-Verlet scheme is used. Let M = n MC × M 0 be the positive integer in which n MC and M 0 have been introduced in Section 2.6. The reduced-order ISDE is solved on the finite interval R = [0 , M ∆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 + 1 sampling points r such that r = ∆r for = 0, . . . , M . The following notations are introduced:

[z 1

 1 ar ], . . . , [z nMC ar ] of random matrix [Z] For m, ε, and ζ fixed, we introduce the Markov stochastic process

  defined by Eq. (2). The random matrix[N ] with values in M ν,N is written as [N ] = [N 1 . . . N N ] in which the columns N 1 , . . . , N N are N independent copies of the normalized Gaussian vector N with values in R ν (this means that E{N } = 0 and E{N N T } = [I ν ]). The random matrices [H d ] and [N ], and the normalized Wiener process {W(r), r ≥ 0} are assumed to be independent.

  , and [Y(0; θ)] = [N (θ)] [a] are independent samples of the stochastic process {[W(r)], r ≥ 0}, the random matrix [Z(0)], and the random matrix [Y(0)]. Let {([Z(r; θ)], [Y(r; θ)]), r ∈ R + } be the corresponding sample of the unique stationary diffusion process {([Z(r)], [Y(r)]), r ∈ R + } of the problem defined by Eqs. (

(

  ii) Assuming that matrix [ c ] is positive definite, computing the Cholesky factorization [ c ] = [ L] T [ L] and computing the re-normalized samples {[η ar ], = 1, . . . n MC } such that, for all = 1, ..., n MC , j = 1, . . . , N , and k = 1, . . . , ν,

  [a] related to the manifold. Let {[∆W ], = 1, . . . , M } be the sequence of the independent Gaussian centered random matrices with values in M ν,m , which are defined by Eq. (C.2) of Appendix C, [∆W ] = [∆W ] [a] .

  . Orthonormal multivariate polynomials Definition of the multi-indices for the PCE of [Z] with respect to [Ξ]. Let N g and µ be fixed such that 1 ≤ N g ≤ ν and 1 ≤ µ ≤ m. For α ∈ {1, . . . , µ}, let β α = (β α 1 , . . . , β α Ng ) be the multi-index such that β α ∈ N Ng and let be |β α | = β α 1 + . . . + β α Ng . Relatively to the columns Ξ 1 , . . . , Ξ µ of matrix [Ξ], we introduce the multi-index β = (β 1 , . . . , β µ ) that belongs to N µNg = N Ng × . . . × N Ng such that |β| = |β 1 | + . . . + |β µ |. For a fixed value of N d that is the maximum degree of the orthonormal multivariate polynomials, such that N d ≥ 1, the following set B N d ,Ng,µ of multi-indices is introduced,

  [Z]. (i) Truncated PCE and convergence. For N g ≤ ν, µ ≤ m, and for N d > 1 fixed, the truncated PCE denoted by [Z(N d , N g , µ)] of random matrix [Z] with values in M ν,m , with respect to random matrix [Ξ] with values in M Ng,µ , is written as

(

  ii) Comments about the L 2 -convergence of the truncated PCE of [Z] when the germ is chosen as [Ξ]. (ii-1) As ν and m are finite, for N g = ν and µ = m, the L 2 -convergence of the sequence {[Z(N d , ν, m)]} N d of random matrices can be proved as follows. Let [Ξ(ν, m)] be the Gaussian random matrix with values in M ν,m defined in Section 3.2. From the construction of the M ν,m -valued random matrix [Z] presented in Section 2.5, [Z] can be written as [

3. 7 .

 7 Quantification of the L 2 -convergence The quantification of the L 2 -convergence of the sequence of random matrices {[Z(N d , N g , µ)]} N d ,Ng,µ towards [Z] (see Eq.(

Figure 1 :

 1 Figure 1: 230 given data points (left). Eigenvalues of the transition matrix for random vector H (right).

5. 1 .

 1 Application 1: Dimension n = 2 with N = 230 given data points

Figure 2 :

 2 Figure 2: pdf for random variables X 1 (solid line) and X 2 (dashed line) obtained by a nonparametric estimation from data points (left). 230 given data points (blue symbols) and 9, 200 additional samples (red symbols) generated using the reduced-order ISDE with m = 3 (right).

Figure 3 :

 3 Figure 3: Graph of the error function N d → error [Z] (N d , N g , µ) for N g = 1 (dashed line) and for N g = 2 (solid line), and for µ = m = 3.

Figure 4 :

 4 Figure 4: For N d = 5, N g = 2, and µ = m = 3, graphs of the reference pdf (dashed line) and of the pdf computed with the analytical representation (solid line) for random variables X 1 (left figure) and X 2 (right figure). In each figure, the dashed line and the solid line are superimposed.

Figure 5 :

 5 Figure 5: 230 given data points (blue symbols) and 9, 200 additional samples (red symbols) generated using the analytical representation.

Figure 6 :

 6 Figure 6: 400 given data points (left). Eigenvalues of the transition matrix for random vector H(right).

X 1 ,

 1 X 2 , and X 3 computed with a nonparametric estimation from the data points. For generating 8, 000 additional samples, the numerical values of the parameters are f 0 = 1.5, ∆r = 0.1179, M 0 = 110, and n MC = 20, yielding 2, 200 for the parameter M defined in Appendix C. The additional samples [x 1 ar ], . . . , [x nMC ar ] (with n MC = 20) are computed by using Section 2.6 and are displayed in Fig.

Fig. 8

 8 displays the graph of the error function (N d , N g ) → error [Z] (N d , N g , µ) defined by Eq. (50) for N d = 1, . . . , 4, for N g = 1, 2, 3, and

Figure 7 :

 7 Figure7: pdf for random variables X 1 (solid line), X 2 (dashed line), and X 3 (dotted line) obtained by a nonparametric estimation from data points (left). 400 given data points (blue symbols) and 8, 000 additional samples (red symbols) generated using the reduced-order ISDE with m = 4 (right).

Figure 8 :

 8 Figure 8: Graph of the error function N d → error [Z] (N d , N g , µ) for N g = 1 (dotted line), N g = 2 (dashed line), N g = 3 (solid line), and for µ = m = 4. In each figure, the dashed line and the solid line are superimposed.

Figure 9 :

 9 Figure 9: For N d = 4, N g = 3, and µ = m = 4, graphs of the reference pdf (dashed line) and of the pdf computed with the analytical representation (solid line) for random variables X 1 (left figure), X 2 (central figure), and X 3 (right figure).

Figure 10 :

 10 Figure 10: 400 given data points (blue symbols) and 8, 000 additional samples (red symbols) generated using the analytical representation.

Figure 11 :

 11 Figure 11: 13, 056 given data points viewed from coordinates x 16 and x 28 (left), viewed from coordinates x 27 and x 28 (center), and viewed from coordinates x 30 , x 32 , and x 33 (right).

Figure 12 :

 12 Figure 12: Eigenvalues of the transition matrix for random vector H (left). Graph m → e red (m) in log 10 scale (right).

Figure 13 :

 13 Figure 13: Left figures: pdf for some components of random vector X estimated from the data points [x d ] (dashed lines) and estimated with the simulated samples [x 1 ar ], . . . , [x nMC ar ] (solid lines). Right figures: 13, 056 given data points (blue symbols) and 39, 168 simulated samples (red symbols).

N

  g and N d is achieved for µ = 8. No additional computation are thus performed for 8 < µ ≤ 50. In order to limit the volume of data presented in the figures, only n MC = 3 samples of random matrix [X] with values in M 35,13056 are displayed, yielding 39, 168 samples of random vector X. The simulated samples [x 1 ar ], . . . , [x nMC ar ] of [X] are

Figure 14 :

 14 Figure 14: Graph of the error function N d → error [Z] (N d , N g , µ) for N g = 1 (dotted line), N g = 2 (dashed line), N g = 3 (thin solid line), and for µ = 8.

  Fig. 14 displays the graph of the error function (N d , N g ) → error [Z] (N d , N g , µ) defined by Eq. (50) for N d = 1, . . . , 4, for N g = 1, . . . , 3, and for µ = 8 (for the reason given before). For N d = 4, N g = 3, and µ = 8, the number of vector-valued coefficients is K(N d , N g , µ) = 20, 475 and the value of the error function (defined by Eq. (

Figure 15 :

 15 Figure 15: Left figures: pdf for some components of random vector X estimated without the analytical representation (dashed lines) and estimated with the simulated samples [x 1 ar ], . . . , [x nMC ar ] (solid lines). Right figures: 13, 056 given data points (blue symbols) and 39, 168 simulated samples (red symbols).

  [Z ] = [Z(r )], [Y ] = [Y(r )], and [W ] = [W(r )], for = 0, . . . , M , with[Z 0 ] = [H d ] [a] , [Y 0 ] = [N ] [a] , [W 0 ] = [0 ν,m ] a.s . (C.1) For = 0, . . . , M -1, let [∆W +1 ] = [∆W +1 ] [a] ,(C.2) be the sequence of random matrices with values in M ν,m , in which [∆W +1 ] = [W +1 ] -[W ]. The increments [∆W 1 ], . . . , [∆W M ] are M independent random matrices with values in M ν,N . For all k = 1, . . . , ν and for all j = 1, . . . , N , the

  diffusion-map vectors related to the manifold. The quantity Ψ β ([Ξ]) is a shortand notation for the non-separable orthonormal multivariate polynomials Ψ β 1 ,...β µ (Ξ 1 , . . . , Ξ µ ), in which Ξ 1 , . . . , Ξ µ are the µ columns of random matrix [Ξ] and are mutually dependent Gaussian R Ng -valued random variables, β α = {β α 1 , . . . , β α Ng }. Denoting the maximum degree of the polynomials as N d (i.e. |β 1 | + . . . + |β µ | ≤ N d with |β α | = β α 1 + . . . + β α Ng ), we denote the PCE of random matrix [Z] by [Z(N d , N g , µ)].• Since the germ of the reduced-order ISDE (which allows for generating the T , = 1, . . . , n MC } of [H], whose columns will be concentrated on S ν .

	• The PCE [Z(N d , N g , µ)] of [Z] allows for generating the samples [z 1 ar ], . . . , [z nMC ar ] of random matrix [Z] and therefore, for deducing the samples {[η ar ] =
	[z ar ] [g]

• The PCE of random matrix [Z] is constructed relative to a basis of multivariate polynomials Ψ β ([Ξ]) orthonormal with respect to the Gaussian centered measure of random matrix [Ξ] for which the given covariance tensor depends on the samples for [Z]) is statistically dependent on the germ [Ξ] that is used for constructing [Z(N d , N g , µ)], the coefficients of [Z(N d , N g , µ)] can easily be computed. • The analytical representation of random matrix [H] with values in M ν,N is then obtained by replacing [Z] by [Z(N d , N g , µ)] in the reduced-order representation [H] = [Z] [g] T , for which the germ is the random matrix [Ξ] with values in M Ng,µ .

  The representation defined by Eq. (1) has been constructed in order that the empirical estimates of the mean value and the covariance matrix of random vector H are 0 ν and [I ν ], respectively. As the generator of random matrix [H] defined by Eqs. (7) to (14) can introduce a bias induced by the integration scheme defined by Eqs. (C.4) to (C.7) in Appendix C, the samples {[η ar ], = 1, . . . n MC } defined by Eq. (15) are re-normalized using the following algorithm:(i) Computation of the empirical estimate m = ( m 1 , . . . , m

	•
	1 ar ], . . . , [z nMC ar ] for random
	matrix [Z] .
	• If integer M 0 is chosen sufficiently large (such that ρ is much larger than
	the relaxation time of the dissipative Hamiltonian dynamical system), then
	[z 1 ar ], . . . , [z nMC ar ] can approximatively be considered as independent samples
	of random matrix [Z].

ν ) ∈ R ν of the mean vector of H and of the empirical estimate [ c ] ∈ M ν of the covariance matrix of H, such that, for all k and k in {1, . . . , ν},

sentation for n MC = 25, 000). It can be seen that the convergence (in measure) of the analytical representation of X is good (the dashed lines are superimposed to the solid lines). In Fig. 15, each one of the three right figures displays the 13, 056 given data points [x d ] and the 39, 168 samples corresponding to the 3-first samples [x 1 PCE,sim ], . . . , [x 3 PCE,sim ] generated with the analytical representation. It can be seen that the samples computed with the analytical representation are effectively concentrated in subset S n .

Conclusions

Starting with a dataset concentrated around a manifold, a new methodology has been presented and validated for constructing a probabilistic characterization of the dataset in the form of a polynomial chaos expansion. Mathematically, this takes the form of a known affine transformation of a matrix-valued random variable for which its polynomial chaos expansion is constructed and is concentrated on the manifold that is identified from the dataset. The proposed methodology is robust and can be used for high dimension and for large given datasets. In the first article, we have proposed a method for constructing a generator of samples on the manifold on which the dataset lives. In this paper, we have proposed an analytical characterization of the dataset on the manifold, which also allows for generating samples on the manifold, which is easy to implement. In addition, a germ is characterized on the manifold, which allows other processes to be constructed on the manifold and which also allows for implementing the spectral methods for any nonlinear transformations defined on the manifold.

Acknowledgments

This research was supported by the ScramJet-UQ project funded under DARPA's EQUIPS Program.

Appendix A. Construction of the diffusion-maps basis

In this appendix, we summarize the construction of the diffusion map-basis based on [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF][START_REF] Coifman | Diffusion maps[END_REF] and detailed in [1]. Let [b] be the positive-definite diagonal real matrix in

Appendix B. Criterion for estimating an optimal value of m

In this Appendix, we recall the criterion introduced in [1] for estimating a value of dimension m. Let [x d ] ∈ M n,N be the matrix of the dataset introduced in Section 2.1 and let [η d ] ∈ M ν,N be the matrix computed with Eq. ( 2). We then introduced the matrix [x red (m)] ∈ M n,N such that (see Eqs. ( 4) and ( 6)), For = 0, . . . , M -1, the Störmer-Verlet scheme applied to Eqs. ( 7) and ( 8) yields

with the initial condition defined by [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF], where b = f 0 ∆r /4, and where [L + 1 2 ] is the M ν,m -valued random variable such that


Shahrokh Ghavamian

Gilles Pijaudier-Cabot

Jacky Mazars

Error Indicator to Assess Quality of Finite-Element Frame Analyses

Keywords: Finite elements, Frames, Analysis, Computer software, Errors

The purpose of this article is to describe a tool developed in order to equip a finite-element software based on nonlinear beam analysis with an error indicator aimed at measuring discretization errors. The technique is an extension of the error estimate devised by Ladeve `ze et al. in 1991 in which the finite-element solution is compared to a statically admissible distribution of the generalized stresses. This last distribution is built element by element which is fast and easy to implement. A nonlinear analysis of a single bay frame in static's is presented as an example.

Introduction

Advances achieved in construction engineering and more restrictive safety requirements, especially with respect to severe and exceptional loads such as earthquakes, have induced major needs in understanding and predicting the response of civil engineering structures. In most instances, the validity of computational results and their efficiency are of great concern. Quality is, however, quite intricate to evaluate in a broad sense because several factors are important and interact. A finite-element model is based on a series of assumptions, with respect to its geometry, to the boundary conditions, to the representation of the applied loads, to the construction process, and also to the constitutive relations of the different materials. The purpose of this technical note is to present a simple tool for evaluating one of these components: the quality of the discretisation of the finite-element ͑FE͒ model.

Many techniques can be found in the literature, related to adaptive meshing and mesh optimization. A posteriori error indicators are mostly elaborated from four different approaches. The first one relies on the analysis of the stress field in the discretised structure ͑see, e.g., Zienckiewicz and Zhu 1987͒. Two sets of results are compared: one being discrete, the other being smoothened. The method devised by Babuska and Rheinholdt ͑1982͒ is based on the measure of the residual forces in the equations of equilibrium which are never totally balanced in nonlinear FE analyses. Huerta and Diez ͑2000͒ use a comparison between two finite-element models. A very fine discretization is considered to be the reference ͑i.e., quasi-exact͒ solution, against which the user's mesh is compared. Finally, Ladeve `ze et al. ͑1991͒ define an error indicator as the distance between the finite-element solution and a statically admissible solution.

Most of these indicators provide an upper bound of the true error in elastic analyses. It is not rigorously demonstrated that they provide an upper bound of the error in general nonlinear analyses but they are convergent, i.e., equal to zero when the exact solution is found.

In this short technical note, we show that the technique due to Ladeve `ze et al. can be easily implemented in finite-element analyses of frames.

Finite-Element Model

The finite-element model used here is a nonlinear static and dynamic computer program based on a layered finite-element description of reinforced concrete beams and frames ͑Bazant et al. 1987͒. This program uses Euler-Bernoulli beam elements. The two-dimensional beam elements possess two nodes. The horizontal displacement u x ͑parallel to the neutral axis of the beam͒ is interpolated with a linear function over the element. The vertical displacement u y ͑perpendicular to the neutral axis of the beam͒ is interpolated with a classical third order hermitian function. These displacements are computed at the centroid of the cross section:

ͫ u x ͑ x ͒ u y ͑ x ͒ ͬϭN ͫ u i u j i j i j ͬ (1)
where N denotes the matrix of shape functions, expressed in term of the element neutral axis coordinate x. u i , v i , and i are the horizontal displacement, the vertical displacement and the rotation at node i, respectively. According to the Euler-Bernoulli theory, the horizontal displacement at any point of the cross section of coordinates ͑x,y͒ is u(x,y)ϭu x (x)Ϫyץu y (x)/ץx and the nonzero component of the strain ͑in direction x͒ is (x,y) ϭץu(x,y)/ץx. The beam cross section is decomposed into layers made of material with nonlinear uniaxial stress-strain laws. The constitutive model used for concrete is a rate independent damage model. A complete description of the constitutive equations can be found

Error Indicator

The finite-element solution provides the values of the generalized stresses, i.e., the normal and shear forces (N k ,V k ) and the bending moment M k at each nodal point. The upper-script k indicates that these expressions are obtained from a displacement-based ͑finite element͒ approximation.

The error estimation consists in measuring the difference between two generalized stress fields. The first one is the finiteelement ͑FE͒ solution. It is kinematically admissible, which means that it satisfies the displacement boundary conditions, compatibility conditions, but that the generalized stresses field derived with the help of the constitutive relations do not satisfy the equilibrium equations pointwise. The second one is statically admissible. It satisfies the equilibrium equations pointwise, and the force boundary conditions as well. With the help of the constitutive relations, strains can be derived from this field but they are, in general, not compatible with a displacement field and with the displacement boundary conditions. Theoretically, if those two solutions are identical, all the governing equations of the problem are satisfied and the exact solution has been found. The distance between the kinematically ͑FE͒ and statically admissible fields is an indicator of the error due to mesh discretisation. It does not mean that the FE solution is better than the statically admissible solution or conversely, any solution can be used for design purposes, but both are approximates.

Statically admissible generalized stress field

The statically admissible solution ͑generalized stress͒ denoted (N s ,V s ,M s ) should verify exactly the differential equations of equilibrium/motion:

dN s ͑ x,t ͒ dx ϩ f x ϭ d 2 ͓u x ͑ x,t ͔͒ dt 2 (2a) dV s ͑ x,t ͒ dx ϩ f y ϭ d 2 ͓u y ͑ x,t ͔͒ dt 2 (2b) dM s ͑ x,t ͒ dx ϩV s ͑ x,t ͒ϭ0 ( 2 c)
where ( f x , f y ) are, for most applications, constant distributed forces, ϭmass per unit length of beam, and rotational inertia terms have been neglected.

For the sake of simplicity, let us omit the inertia terms and focus on statics. Equation ͑2͒ may be integrated over each beam element ͑recall that distributed forces are constant͒:

N s ͑ x ͒ϭaϩ f x x (3a) M s ͑ x ͒ϭ ͑ x ͒ 2 2 ϩbxϩc (3b)
Constants ͑a, b, c͒ need to be computed. Since each finite element is in equilibrium, according to the FE solution, nodal values in Eq. ͑3͒ are equated with the FE solution. Constant a is computed from Eq. ͑3a͒ by expressing this relation at either nodes of the element with the same result. Constants ͑b, c͒ are computed from Eq. ͑3b͒ expressed at each node:

N s ͑ 0 ͒ϭaϭN k ͑ 0 ͒ M s ͑ 0 ͒ϭcϭM k ͑ 0 ͒ (4) M s ͑ 1 ͒ϭ L 2 2 ϩbLϩcϭM k ͑ 1 ͒
where Lϭelement length. Within each beam element, Eq. ͑4͒ provides a statically admissible solution. It is entirely local, i.e., expressed at the element level. Hence its resolution is a post processing operation which is performed element per element successively. The same method applies to more sophisticated load distributions without any modifications.

It is important also to satisfy equilibrium across each beam element. As opposed to the case of two-dimensional finite elements where stress vectors are not continuous across each element, this is automatically satisfied in the beam formulation ͑because of the order of continuity of the interpolation functions͒. This is the major reason why the statically admissible solution (N s ,V s ,M s ) is so easy to construct. Note that another technique due to Carol and Murcia ͑1989͒ could have been used for obtaining statically admissible fields. It would have been more computer time consuming because a global set of equations of size equal to the number of degree of freedom of the problem needs to be solved.

In dynamics, inertial forces must be taken into account in the derivation of the «statically» admissible field of generalized forces. With a classical Newmark integration scheme, which assumes a constant acceleration over each time step, these inertia forces can be viewed as additional distributed loads which are constant during each time step and follow the interpolation of the displacements in space. Hence the «statically» admissible field of normal force becomes a polynomial of order 2 and the «statically» admissible distribution of bending moment becomes a polynomial of order 5. The derivation of the statically admissible solution is exactly the same as in static's, except that the external loads are slightly more complex.

Definition of the error indicator

We have now two solutions, which can be compared. The distance E between the two solutions is defined as an L 2 norm according to Ladeve `ze et al. ͑1991͒. E is further denoted as the global error indicator. It is an approximation of the true error

͐ elements P ͱ͑N s ͑ x ͒ϪN k ͑ x ͒͒ 2 dx ͐ structure ͱ͑N k ͑ x ͒ϩN s ͑ x ͒͒ 2 dx ϩ ͐ element p ͱ͑M s ͑ x ͒ϪM k ͑ x ͒͒ 2 dx ͐ structure ͱ͑M k ͑ x ͒ϩM s ͑ x ͒͒ 2 dx (5)
where e p is the local error indicator in element p ͑distance between the two solutions in the element͒. Numerically, the local error indicator is computed over the element according to the quadrature of the finite-element solution. We may retain, in each element, the maximum value of the indicator over the entire calculation. Consequently, the local and global error indicators are either constant or monotonically increasing in the course of the calculation. The above definitions can be used in statics and in dynamics.

Application

As an illustrative example, consider the single bay plane structure shown in Fig. 2. It is made of reinforced concrete. The reinforcement is constant over the frame; it is incorporated in the beam element assuming a perfect bond between steel and concrete. The loading is displacement controlled: a horizontal displacement is applied to the left corner node. The maximum displacement is 0.5 m, decomposed into a loading history of 100 increments of equal size. Four different meshes shown in Fig. 3 are considered. The number of finite elements ranges from 3 to 25. The distributions of tensile damage, of the plastic strain in the reinforcement, and of the local error indicator over the structure at the end of the loading process are also shown in this figure .  Fig. 4 shows the load-deflection curves ͑horizontal force versus horizontal displacement at the top left corner͒ for each finite- element mesh. It can be observed that they are quite similar for the two finest meshes. In fact the FE solution seems to have converged toward the exact one, in terms of mesh refinement. In terms of local quantities ͑damage, plastic strain, and local errors͒, there are still differences as shown in Fig. 3. The minimum and maximum local errors keep on diminishing upon refinement. The maximum local error, however, concentrates in the joints where material nonlinearity is more intense, as expected. According to this analysis, the global response can be obtained quite accurately with a moderately dense mesh, local quantities such as damage or plastic strain are much more demanding from the viewpoint of the finite-element discretization.

Conclusion

An extension of the error indicator devised by Ladeve `ze et al. ͑1991͒ to the simple case of beam analysis has been presented.

The method relies on the comparison between the finite-element solution and a statically admissible field of generalized stresses. It is quite simple to implement because statically admissible fields of generalized stresses are constructed element per element in a post-processing course of the FE calculation. It can be easily extended to a three-dimensional configuration, or to cases where distributed forces are arbitrary ͑known͒ functions, such as in dynamics.

Two forms of error estimate can be considered: the global error which is an overall measure of the quality of the finiteelement discretization and the distribution of the local errors which pin-points locations where the discretization is too coarse or could be coarser. This latter information could serve as a point of entry in the design of optimized meshes.

1

  LMT-Cachan, ENS Cachan, 61 avenue du Président Wilson F-94235 Cachan cedex, France 2 R&DO, Laboratoire de Génie Civil de Nantes-Saint Nazaire, Ecole Centrale de Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France in Dube ´et al. ͑1996͒. Fig.1shows the concrete response for tension-compression cycles. For the reinforcements, a onedimensional elastoplastic model with linear hardening is implemented. It should be pointed out that any other constitutive relations could be implemented as well, without any modification of the calculation of the error indicator described in the next section.

Fig. 1 .

 1 Fig. 1. Response of concrete constitutive law to tension-compression loading

Fig. 2 .Fig. 3 .

 23 Fig. 2. Single bay case study