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

Quy-Dong To
email: quy-dong.to@univ-eiffel.fr

Polynomial-chaos-based conditional statistics for probabilistic learning with heterogeneous data applied to atomic collisions of Helium on graphite substrate

Keywords: Polynomial chaos expansion, Statistical surrogate model, Probabilistic learning, Heterogeneous data, Uncertainty quantification, Atomic collisions

A formulation and an algorithm are presented to construct a truncated polynomial chaos representation of a vectorvalued random output. This representation depends on a vector-valued random input with a known probability measure and a vector-valued random latent variable with an unknown probability measure. The construction of this PCE representation relies solely on a training set comprising a small number of independent realizations of the non-Gaussian dependent random output and input vectors. The training set consists of heterogeneous data, which poses challenges in accurately estimating the chaos coefficients. Despite the heterogeneity of the data, the proposed formulation and algorithm allow for the construction of a highly accurate global surrogate model. Additionally, we propose an alternative approach by constructing a surrogate model based on prior separation of the heterogeneous dataset into subsets, each containing "quasi-homogeneous" data. The separation method is designed to account for a partial overlap of the probability measure supports associated with the subsets. The identification of the PCE is performed offline. By utilizing the PCE, a fast online surrogate model is obtained, enabling analysis of large dynamical systems beyond the computational capabilities currently available. An application to atomic collisions of Helium on a graphite substrate is presented, where the training set was generated by Molecular Dynamics simulations done in a previous paper. The obtained results demonstrate accuracy of the proposed approach.

Introduction

Physics problem addressed in the paper and its statistical surrogate model. To solve gas flow problems in engineering applications, a variety of simulation methods have been developed in the literature. They can be classified into two main types: continuum-based methods (Navier Stokes, moment equations, Burnette etc) and particle-based methods (Direct Simulation Monte Carlo, Lattice Boltzmann, Molecular Dynamics, etc.). In addition to the bulk behavior representing fluid-fluid interaction, it must be completed by the interaction between the fluid and the solid boundary. To avoid the huge computation cost relating to the modeling the solid phase, the interaction fluid-solid is usually substituted by statistical surrogate models. The construction of the latter can be done by studying separately the gas-wall collisions using the Molecular Dynamics. The presented application is the case of Helium and graphite at low temperature where complex phenomenon like adsorption and surface diffusion are present and dominant. However, the developed methodology is general and can be used for any complex dynamical systems and in particular, in the context of the application presented, can be applied to any gas-wall couple. The input vector to this statistical surrogate model consists of a realization (sample) w j = v j in of the random velocity vector W (the control variable) representing the velocity of the incident particle on the layer. The output vector is the corresponding realization q j of the random vector Q = (V out , ∆ t , D x , D y ). This vector includes the random reflected velocity vector V out (output) of the particle, the random absorption duration ∆ t (residence time), and the two displacement components D x and D y in the (oxy) plane within the layer. There exist two primary physical regimes depending on the realization w j of W. In one regime, the particle undergoes quasi-reflection by the layer, resulting in a short absorption time. In the other regime, the particle is absorbed within the layer and then emerges after a more or less extended random duration and with random displacements in the plane of the layer. On the other hand, the dynamical system is extremely complex, and the state Q cannot be solely explained by W. This implies that a deterministic mapping such as Q = f(W) does not exist. Instead, there exists a random mapping F such that Q = F(W). To construct the statistical surrogate model, it is necessary to introduce a hidden explanatory vector-valued random variable (latent variable). We can introduce a random vector U of unknown dimension, which is assumed to be statistically independent of W. Thus, the relationship for Q can be expressed as Q = f(W, U). The only available information consists of the realizations {(q j d , w j d ), j = 1, . . . , n d } of (Q, W) obtained from the MD simulations, where n d is relatively small. The objective is to construct a truncated polynomial chaos representation, Q chaos = f chaos (W, U), of Q = F(W), which defines the statistical surrogate model for Q. A probability measure P U (du) for U is then constructed to make the probability measure of Q chaos as close as possible to the probability measure of Q. This representation allows us to generate a realization q 0 of Q chaos ≃ Q corresponding to a given realization w 0 of W, ensuring that (q 0 , w 0 ) is a consistent realization with respect to the probability measure of (Q, W). Therefore, computing the realization q 0 of Q chaos , given W = w 0 , is done quickly using q 0 = f chaos (w 0 , u 0 ), where u 0 is any realization of U (which is independent of w 0 ).

Main difficulties related to the construction of the statistical surrogate model. There are three main difficulties. The first is that only a small training dataset {(q j d , w j d ), j = 1, . . . , n d } of dimension n d is available. The number of points, n d , is too small to construct the chaos representation Q chaos of Q. Therefore, it is necessary to generate a large learned dataset {(q ℓ , w ℓ ), ℓ = 1, . . . , N} consisting of N ≫ n d learned realizations that follow the probability measure of (Q, W). This dataset will be created using only the information provided by the training dataset. The second problem is associated with the presence of the hidden random vector U, of which we do not know the dimension and its probability measure. The third is related to the training dataset, which comprises heterogeneous data. Within this dataset, there is a combination of data related to two different regimes of physics processes. As we have explained, one regime corresponds to the scenario where the particle reflects off the layer, meaning there is no absorption of the particle. In this case, the duration ∆ t and the displacements (D x , D y ) are small. The other regime corresponds to the possibility of particle absorption by the layer. Here, the particle moves inside the layer and emerges after a random duration ∆ t and random displacements (D x , D y ), which can be large. Separating these two probabilistic phenomena (and therefore separating the points of the training dataset into distinct clusters) is challenging due to partial overlap in the supports of the probability measures associated with these two regimes. Consequently, the training dataset consists of heterogeneous data.

Rewriting the addressed physics problem within a broader context. We will present a general methodology that can be applied to other situations. To do this, we reformulate in a broader context, the problem we presented earlier for the specific physics problem in question. We consider a large-scale stochastic computational model that depends on an unknown (and therefore uncontrolled) random parameter, which is modeled by a random variable with value R n u , denoted U (the latent variable). The dimension n u and the probability measure P U (du) of U are both unknown. The control parameter is the random variable W with values in R n w , whose probability measure is P W (dw). We assume that W and U are independent, so we have P W,U (dw, du) = P W (dw) ⊗ P U (du). The quantity of interest is the R n qvalued random variable Q = F(W). Let C w ⊂ R n w be the support of P W (dw). The random mapping F is unknown. The joint probability measure of Q and W is denoted by P Q,W (dq, dw). In such a scenario, the random mapping F can be rewritten as w → F(w) = f(w, U) in which the deterministic mapping f is also unknown. Let X = (Q, W) be the R n x -valued random variable, where n x = n q + n w , and its probability measure is denoted by P X (dx). Random vector X is associated with the random manifold defined by the random graph {(F(w), w), w ∈ C w } and is fully characterized by its probability measure P X (dx), which represents the joint probability measure P Q,W (dq, dw) of Q and W. For this problem, the only available information is the training set

D train (X) = {x j d , j = 1, . . . , n d } consisting of n d independent realizations x j d = (q j d , w j d ) ∈ R n x of the R n x -valued random variable X = (Q, W).
It is assumed that n d is small. For each given realization w j d of W, the corresponding realization q j d of Q is obtained from the computational model. However, there are no corresponding realizations u ℓ d of U available for the given x j d . In addition, as we have explained, it will be necessary to construct a large learned dataset D learn (X) = {x ℓ , ℓ = 1, . . . , N} with N ≫ n d in order to build f chaos . The estimation of the probability measure P X will be performed using the learned dataset. Note that the realization u ℓ of U corresponding to x ℓ is also not usable.

Short overview of the related works on polynomial chaos expansion (PCE) for developing a methodology to solve the problem. Given the statistical dependence between Q and W, the joint probability measure P Q,W (dq, dw) = P X (dx) can be accurately estimated using the large learned dataset D learn (X) as well as the probability measure P W of W. If Q did not depend on U, then the problem would easily be obtained by a projection on the chaos polynomials constructed with P W . When Q depends on U, the coefficients obtained with such a projection become random coefficients yielding a truncated PCE Q chaos with random coefficients for which a detailed analysis has been proposed in [START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF].

The concept of PCE for stochastic processes was first introduced by Wiener and Cameron in their seminal works [START_REF] Wiener | The homogeneous chaos[END_REF][START_REF] Cameron | The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals[END_REF], while Ghanem and colleagues pioneered an effective Karhunen-Loève-based construction for random fields [START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Ghanem | Stochastic Finite a Spectral Approach[END_REF]. The Wiener-Askey PCE was employed by Xiu [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF], and the development of random fields in polynomial chaos for arbitrary probability measures was introduced by Soize [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]. The PCE with random coefficients were explored by [START_REF] Lucor | Generalized polynomial chaos and random oscillators[END_REF][START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF][START_REF] Dolgov | Polynomial chaos expansion of random coefficients and the solution of stochastic partial differential equations in the tensor train format[END_REF], while Tipireddy presented a basis adaptation in homogeneous chaos spaces [START_REF] Tipireddy | Basis adaptation in homogeneous chaos spaces[END_REF]. A compressed principal component analysis of non-Gaussian vectors using symmetric polynomial chaos was proposed by Mignolet [START_REF] Mignolet | Compressed principal component analysis of non-Gaussian vectors[END_REF]. Significant works are also devoted to the acceleration of stochastic convergence of 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]. Polynomial chaos expansions have been and continue to be intensively used in both finite and infinite dimensions for uncertainty modeling and propagation [START_REF] Debusschere | Numerical challenges in the use of polynomial chaos representations for stochastic processes[END_REF][START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF][START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF][START_REF] Das | Polynomial chaos representation of spatio-temporal random fields from experimental measurements[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] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] (see also hereinafter the stochastic solvers and the stochastic finite elements).

After constructing probabilistic models of uncertainties, it becomes essential to investigate how these uncertainties propagate within systems. This requires the use of methods for solving stochastic equations. The initial set of methods is grounded on Monte Carlo numerical simulation techniques. The second set is based on spectral projection methods [START_REF] Babuska | Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation[END_REF][START_REF] Matthies | Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations[END_REF][START_REF] Babuska | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: Invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | Generalized spectral decomposition for stochastic nonlinear problems[END_REF][START_REF] Nouy | Proper generalized decomposition and separated representations for the numerical solution of high dimensional stochastic problems[END_REF], such as those based on polynomial chaos expansions [START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF][START_REF] Ghanem | Stochastic Finite a Spectral Approach[END_REF][START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Najm | Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics[END_REF] and called stochastic finite element method when the discretization method of the boundary value problems are performed using the finite element method [START_REF] Ghanem | Stochastic Finite a Spectral Approach[END_REF][START_REF] Ghanem | Numerical solution of spectral stochastic finite element systems[END_REF][START_REF] Ghanem | Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach[END_REF][START_REF] Ghanem | Ingredients for a general purpose stochastic finite elements formulation[END_REF][START_REF] Pellissetti | Iterative solution of systems of linear equations arising in the context of stochastic finite elements[END_REF][START_REF] Deb | Solution of stochastic partial differential equations using galerkin finite element techniques[END_REF][START_REF] Ghanem | Stochastic Finite Elements: A spectral Approach[END_REF] and also [START_REF] Frauenfelder | Finite elements for elliptic problems with stochastic coefficients[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF][START_REF] Xu | A multiscale stochastic finite element method on elliptic problems involving uncertainties[END_REF][START_REF] Matthies | Stochastic finite elements: Computational approaches to stochastic partial differential equations[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Luthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF].

Most often, the uncertainties probability model is a prior model. If targets are available for observations of the system, coming from experimental measurements or from more precise numerical simulations, a posterior probability model of uncertainties can be estimated by solving inverse statistical problems based on the maximum likelihood, the Bayesian inference, and machine learning. For general overviews on statistical inverse methods, see [43,[START_REF] Kennedy | Bayesian calibration of computer models[END_REF][START_REF] Tarantola | Inverse Problem Theory And Methods For Model Parameter Estimation[END_REF][START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF][START_REF] Owhadi | On the brittleness of Bayesian inference[END_REF][START_REF] Matthies | Inverse problems in a Bayesian setting[END_REF][START_REF] Dashti | The Bayesian approach to inverse problems[END_REF][START_REF] Soize | Sampling of Bayesian posteriors with a non-Gaussian probabilistic learning on manifolds from a small dataset[END_REF], and for complements, see [START_REF] Fearnhead | Exact and efficient Bayesian inference for multiple changepoint problems[END_REF][START_REF] Golightly | Bayesian sequential inference for nonlinear multivariate diffusions[END_REF][START_REF] Zabaras | A scalable framework for the solution of stochastic inverse problems using a sparse grid collocation approach[END_REF][START_REF] Ma | An efficient Bayesian inference approach to inverse problems based on an adaptive sparse grid collocation method[END_REF][START_REF] Flath | Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial hessian approximations[END_REF][START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF][START_REF] Perrin | Karhunen-loève expansion revisited for vector-valued random fields: Scaling, errors and optimal basis[END_REF][START_REF] Najm | Inference given summary statistics[END_REF][START_REF] Tsilifis | Bayesian adaptation of chaos representations using variational inference and sampling on geodesics[END_REF][START_REF] Zhou | An approximate empirical Bayesian method for large-scale linear-Gaussian inverse problems[END_REF][START_REF] Perrin | Adaptive method for indirect identification of the statistical properties of random fields in a Bayesian framework[END_REF]. The statistical identification of the coefficients of polynomial chaos representations of random fields can be found in [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Das | Asymptotic sampling distribution for polynomial chaos representation from data: a maximum entropy and fisher information approach[END_REF][START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansions[END_REF], in particular in [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] Soize | A computational inverse method for identification of non-Gaussian random fields using the Bayesian approach in very high dimension[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF] for high dimension, and for representations of random vectors in [START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF][START_REF] Madankan | Polynomial-chaos-based Bayesian approach for state and parameter estimations[END_REF][START_REF] Chen-Charpentier | Parameter estimation using polynomial chaos and maximum likelihood[END_REF][START_REF] Elsheikh | Efficient bayesian inference of subsurface flow models using nested sampling and sparse polynomial chaos surrogates[END_REF][START_REF] Nagel | Spectral likelihood expansions for Bayesian inference[END_REF][START_REF] Sraj | Coordinate transformation and polynomial chaos for the Bayesian inference of a Gaussian process with parametrized prior covariance function[END_REF][START_REF] Shao | Bayesian sparse polynomial chaos expansion for global sensitivity analysis[END_REF]. The inverse identification of random matrices have been proposed in [START_REF] Soize | Probabilistic model identification of uncertainties in computational models for dynamical systems and experimental validation[END_REF][START_REF] Desceliers | Probabilistic model of the human cortical bone with mechanical alterations in ultrasonic range[END_REF][START_REF] Arnst | Identification and sampling of Bayesian posteriors of high-dimensional symmetric positive-definite matrices for data-driven updating of computational models[END_REF]. Statistical inverse methods are also used to perform model updating [START_REF] Beck | Updating models and their uncertainties. i: Bayesian statistical framework[END_REF][START_REF] Beck | Bayesian updating of structural models and reliability using markov chain monte carlo simulation[END_REF][START_REF] Ching | Bayesian state and parameter estimation of uncertain dynamical systems[END_REF][START_REF] Cheung | Calculation of posterior probabilities for bayesian model class assessment and averaging from posterior samples based on dynamic system data[END_REF], model selection [START_REF] Parussini | Multi-fidelity Gaussian process regression for prediction of random fields[END_REF][START_REF] Jiang | Adaptive bayesian SLOPE: model selection with incomplete data[END_REF], and to construct surrogate models (or metamodels) [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF][START_REF] Constantine | Active subspace methods in theory and practice: applications to kriging surfaces[END_REF][START_REF] Kersaudy | A new surrogate modeling technique combining kriging and polynomial chaos expansions-application to uncertainty analysis in computational dosimetry[END_REF][START_REF] Kleijnen | Regression and kriging metamodels with their experimental designs in simulation: a review[END_REF][START_REF] Giovanis | Data-driven surrogates for high dimensional models using Gaussian process regression on the Grassmann manifold[END_REF][START_REF] Liu | Surrogate modeling based on resampled polynomial chaos expansions[END_REF][START_REF] Zhou | Surrogate modeling of high-dimensional problems via data-driven polynomial chaos expansions and sparse partial least square[END_REF].

The machine learning tools and artificial intelligence [START_REF] Korb | Bayesian artificial intelligence[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] Ghahramani | Probabilistic machine learning and artificial intelligence[END_REF][START_REF] Russel | Artifical Intelligence, A Modern Approach[END_REF] provide methods that make it possible to solve problems in UQ in the field of physics and engineering sciences. These problems could not be solved without these learning methods because the use of the usual methods would require computer resources, which are not available. Regarding these learning methods, let us cite, for example, learning with kernels [START_REF] Hofmann | Kernel methods in machine learning[END_REF][START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF][START_REF] Akian | Learning best kernels from data in gaussian process regression. with application to aerodynamics[END_REF], probabilistic and statistical learning [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] James | An Introduction to Statistical Learning[END_REF][START_REF] Taylor | Statistical learning and selective inference[END_REF][START_REF] Swischuk | Projection-based model reduction: Formulations for physics-based machine learning[END_REF], learning on the manifolds [START_REF] Öztireli | Spectral sampling of manifolds[END_REF][START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Perrin | Nested polynomial trends for the improvement of Gaussian process-based predictors[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF][START_REF] Perrin | Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints[END_REF][START_REF] Soize | Entropy-based closure for probabilistic learning on manifolds[END_REF][START_REF] Kevrekidis | Manifold learning for parameter reduction[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Kontolati | Manifold learning-based polynomial chaos expansions for high-dimensional surrogate models[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF][START_REF] Almeida | A probabilistic learning approach applied to the optimization of wake steering in wind farms[END_REF][START_REF] Zhong | Surrogate modeling of structural seismic response using Probabilistic Learning on Manifolds[END_REF][START_REF] Almeida | Uncertainty quantification of waterflooding in oil reservoirs computational simulations using a probabilistic learning approach[END_REF], and probabilistic physics-based learning [START_REF] Pan | Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability[END_REF][START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF][START_REF] Soize | Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures[END_REF].

Heterogeneous data are ubiquitous in many scientific fields, and their importance continues to grow. As a result, research on heterogeneous data is important, leading to the emergence of new methods and techniques for analysis and modeling. Common methods for analyzing heterogeneous data at various levels of abstraction include regression methods that use polynomial chaos representations, kriging methods that are statistical interpolation techniques, and neural network and deep learning methods that can model complex relationships between input and output variables but require large training datasets. More advanced techniques comprise Hidden Markov Models (HMM) [START_REF] Rabiner | An introduction to hidden Markov models[END_REF] that are particularly useful for modeling time series data with different regimes, Dirichlet Processes (DP) [START_REF] Rodriguez | The nested Dirichlet process[END_REF] that are nonparametric probabilistic models capable of clustering, hierarchical clustering methods that group similar data into clusters [START_REF] Jung | A decision criterion for the optimal number of clusters in hierarchical clustering[END_REF], Probabilistic Graphical Models (PGM) [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF] that are effective for modeling systems with complex variable interactions, Hidden Markov Networks (HMN) [START_REF] Ghahramani | An introduction to hidden Markov models and Bayesian networks[END_REF] that are suitable for modeling complex systems with unknown graph structures, and kernel algorithm [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF].

Novelty. The novelty of this work is directly linked to the main challenges that we have identified. First, we present a formulation and an algorithm to build a surrogate statistical model in the form of a truncated polynomial chaos representation, Q chaos = f chaos (W, U), of the vector random variable Q = F(W). In this context, F represents an unknown random mapping, and the available information solely consists of n d independent realizations of the non-Gaussian dependent random vector (Q, W), where the value of n d is small and the realizations correspond to heterogeneous data. Due to this limited heterogeneous dataset, accurately estimating the chaos coefficients becomes challenging. Furthermore, the vector-valued random variable U represents the latent random variables within F, and it is necessary to determine its dimension and to chose its probability measure. Although the data is heterogeneous, the proposed formulation and algorithm enable the construction of a global surrogate model with excellent accuracy. Although the global surrogate model is highly accurate, we also propose an alternative approach by constructing a surrogate model based on a prior separation of the heterogeneous dataset into subsets, each consisting of "quasi-homogeneous" data. The proposed separation method is developed in the context of the existence of a partial overlap of the supports of probability measures associated with the subsets. It should be noted that in this case, no method can achieve an "exact" separation based solely on the available information.

Organization of the paper. It should be noted that Section 1 presents the physics problem addressed in the paper. It provides a definition of the statistical surrogate model of interest and presents the broader context in which the physics problem is posed. Section 2 focuses on the generation of a large learned dataset from a given small training dataset. To circumvent the numerical difficulties during the construction of the polynomial chaos representations, a normalization and a scaling of the learned dataset are carried out. This involves transforming Q and W into normalized/scaled random variables, denoted as Y and Ξ respectively. Section 3 deals with the construction of the non-separated multivariate polynomial chaos for Ξ, while Section 4 addresses the construction of separated multivariate polynomial chaos for the latent random variable U. In Section 5, we present the polynomial chaos expansion of random vector Y. This expansion is obtained in tensorizing the two Hilbert bases related to Ξ and U. Subsequently, we can construct the truncated polynomial chaos expansion of random vector Y, which serves as the basis for constructing the statistical surrogate model. Section 7 is devoted to the identification of the unknown coefficients in the truncated PCE of Y. The optimization problem involved is nonconvex, with a constraint, and its cost function is based on the Overlap criterion. To simplify the search for an optimal solution, we transform this problem into an unconstrained optimization problem. Section 8 finalizes the first part dedicated to the construction of the global statistical surrogate model, which is based on a representation of Q in polynomial chaos of W and U. In section 9, we present an alternative approach, which consists in using the separation into clusters of the learned dataset, to construct the conditional statistics based on the polynomial chaos. The last section deals with the application to atomic collisions of Helium on graphite substrate.

Notations

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

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

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

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

C: set of all the complex numbers. M n,m : set of the (n × m) real matrices. Convention used for random variables. In this paper, for any finite integer m ≥ 1, the Euclidean space R m is equipped with the σ-algebra B R m . If Y is a R m -valued random variable defined on the probability space (Θ,

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

Generation of a large learned dataset, its normalization and scaling

Generation of a large learned dataset. As described in Section 1, a large learned dataset D learn (X) = {x ℓ , ℓ = 1, . . . , N} is generated from the training dataset D train (X) = {x j q , j = 1, . . . , n d }, where N ≫ n d . The learned dataset is generated using the PLoM algorithm under constraints [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF] to enforce the learned probability measure to match the given mean value and covariance matrix. From this, we obtain the learned realizations {q ℓ , ℓ = 1, . . . , N} and {w ℓ , ℓ = 1, . . . , N} for the random vectors W and Q, respectively, where (q ℓ , w

ℓ ) = x ℓ ∈ R n q × R n w = R n x .
Scaling the control random parameter W to obtain the random vector Ξ. Since W = (W 1 , . . . , W n w ) is a R n w -valued random variable, in order to avoid numerical difficulties while constructing its polynomial chaos, we introduce a scaled

R n w -valued random variable Ξ = (Ξ 1 , . . . , Ξ n w ) such that W = s w (Ξ) in which s w = (s w,1 , . . . , s w,n w ) is the mapping from [-1 , +1] n w into R n w such that, for all k ∈ {1, . . . , n w }, we have w k = s w,k (ξ k ) = a k ξ k + b k in which a k and b k are such that s w,k (-1) = min ℓ w ℓ k and s w,k (+1) = max ℓ w ℓ k . The support of the probability measure P Ξ (dξ) of Ξ on R n w is chosen as the compact subset C ξ = [-1 , +1] n w ⊂ R n w . The N independent realizations {ξ ℓ , ℓ = 1, . . . , N} of Ξ are given by ξ ℓ = s -1 w (w ℓ ).
With such a scaling, the support of the probability measure P W (dw) is defined as the compact subset C w = s w (C ξ ). We then have introduced the mapping,

ξ → w = s w (ξ) : C ξ → R n w such that W = s w (Ξ) .
(2.1)

Normalization and scaling of the random vector Q to obtain random vector Y. To address potential numerical difficulties in constructing the polynomial chaos expansion of Q, we employ a normalization technique that involves a principal component analysis (PCA) followed by scaling. This process yields the R n q -valued random variable Y. Let q be the empirical mean value of Q and [ C Q ] its empirical covariance matrix, which are estimated with the independent realizations {q ℓ , ℓ = 1, . . . , N}. Note that n q ≪ N and we assume that

[ C Q ] belongs to M + n q . Let [V] ∈ M n q be the or- thogonal matrix, such that [V] T [V] = [V] [V] T = [I n q ], constituted of the eigenvectors of matrix [ C Q ] and let [ζ] ∈ M + n q
be the diagonal matrix of the eigenvalues that are all positive. We thus have [

C Q ] = [V] [ζ] [V] T . We then define the normalized R n q -valued random variable R such that Q = q + [V] [ζ] 1/2 R and consequently, R = [ζ] -1/2 [V] T (Q -q). The realizations {r ℓ , ℓ = 1, . . . , N} of R are therefore computed by r ℓ = [ζ] -1/2 [V] T (q ℓ -q).
Consequently, the empirical mean value r and the empirical covariance matrix

[ C R ] estimated with {r ℓ , ℓ = 1, . . . , N} are such that r = 0 n q and [ C R ] = [I n q ]. The random variable R is now scaled in a random variable Y. Let Y = (Y 1 , . . . , Y n q ) be the scaled R n q -valued random variable R = (R 1 , . . . , R n q ) such that for all i ∈ {1, . . . , n q }, we have R i = s i × Y i in which s i = (max ℓ |r ℓ i |) > 0, yielding Y i = R i /s i . The N independent realizations {y ℓ , ℓ = 1, . . . , N} of Y are such that y ℓ i = r ℓ i /s i .
The composition of these two transformations allows the bijective mapping q to be defined,

y → q = q(y) : R n q → R n q such that Q = q(Y) and Y = q -1 (Q) . (2.2)
The empirical estimate of the covariance matrix [

C Y ] ∈ M + n q of Y is such that [ C Y ] ii ′ = (s i ) -2 δ ii ′ . Since R is centered, Y will be centered, and the second-order moment-matrix [M Y ] = E{Y Y T } ∈ M + n q estimated by [M Y ] = (N - 1) -1 N ℓ=1 y ℓ (y ℓ ) T is equal to [ C Y ]. Matrix [M Y ] can then be rewritten as [M Y ] = [L Y ] 2 with [L Y ] ii ′ = (s i ) -1 δ ii ′ . The mapping f from R n w × R n u into R n q , which have defined in Section 1 and which is such that Q = f(W, U), is then transformed in a mapping y from R n w × R n u into R n q such that Y = y(Ξ, U) = q -1 (f(s w (Ξ), U)) .
(2.3)

Non-separated multivariate polynomial chaos for Ξ

Hilbert space H associated with the random vector Ξ. Let Ξ = (Ξ 1 , . . . , Ξ n w ) be the R n w -valued random variable defined in Section 2. Let H = L 2 P Ξ (R n w , R) be the Hilbert space of all the functions from R n w into R, equipped with the inner product ⟨h , h⟩ H = R nw h(ξ) h(ξ) P Ξ (dξ), and the associated norm

∥ h ∥ H = ⟨h , h⟩ 1/2 H . For any h in H, H = h(Ξ) is a second-order real-valued random variable such that E{H 2 } = E{h(Ξ) 2 } = R nw h(ξ) 2 P Ξ (dξ) < +∞.
Non-separated multivariate polynomial chaos in H. Let α = (α 1 , . . . , α n w ) ∈ N n w be the multi-index, which includes α (1) = (0, . . . , 0). For all α in N n w , let Ψ α 1 ,...,α nw (ξ 1 , . . . , ξ n w ), rewritten as Ψ α (ξ), be the multivariate polynomials (nonseparated with respect to ξ 1 , . . . , ξ n w ), which are orthonormal in H, ⟨Ψ α , Ψ β ⟩ H = δ αβ , and such that Ψ α (1)

(ξ) = 1. It is known that {Ψ α , α ∈ N n w } is a Hilbert basis of H.

Polynomial chaos expansion of H = h(Ξ).

For any h in H, the PCE of the second-order random variable H is written as

H = α∈N nw h α Ψ α (Ξ)
, where the series of the PCE converges with respect to the norm of H. The real coefficients h α are such that ∥H∥ 2 H = α∈N nw (h α ) 2 < +∞ and can be calculated by

h α = ⟨H , Ψ α ⟩ H = E{H Ψ α (Ξ)} = E{h(Ξ) Ψ α (Ξ)} = R nw h(ξ) Ψ α (ξ) P Ξ (dξ).
Truncated PCE H chaos = h chaos (Ξ) for representing H = h(Ξ). Let N g be the maximum degree of the considered truncated polynomial chaos expansion. We then have

|α| = α 1 + . . . + α n w ≤ N g . The set {α ∈ N n w , |α| ≤ N g } of all the multi-indices in N n w such that |α| ≤ N g is rewritten as {α (k) = (α (k) 1 , . . . , α (k) n w ) ∈ N n w , k = 1, .
. . , κ} with α (1) = (0, . . . , 0) and where κ = (N g + n w )!/(N g ! n w !). To simplify the notation, the polynomial chaos of multi-index α (k) is rewritten as ψ k (Ξ) = Ψ α (k) (Ξ). We then have,

ψ 1 (Ξ) = 1 , ⟨ψ k , ψ k ′ ⟩ H = δ kk ′ , E{ψ k (Ξ)} = δ 1k . (3.1)
The truncated PCE of H = h(Ξ) is thus given by

H chaos = h chaos (Ξ) = κ k=1 h k ψ k (Ξ) , h k = ⟨H , ψ k ⟩ H , (3.2)
in which h k is the rewriting of h α (k) .

Separated multivariate polynomial chaos for U

As previously mentioned, the vector-valued random variable U is a latent variable for which no information is available. Consequently, in the context of PCE, we choose U as a normalized R n u -valued Gaussian random variable whose probability measure is written as P U (du) = (2π) -n u /2 exp(-∥u∥ 2 /2) du. The dimension n u ≥ 1 is unknown and needs to be determined. We use notations similar to those introduced in Section 3.

Hilbert space G associated with the random vector U. Let G = L 2 P U (R n u , R) be the Hilbert space of all the functions from R n u into R, equipped with the inner product ⟨g , g⟩ G = R nu g(u) g(u) P U (du), and the associated norm

∥ g ∥ G = ⟨g , g⟩ 1/2 G . For any g in G, G = g(U) is a second-order real-valued random variable such that E{G 2 } = E{g(U) 2 } = R nu g(u) 2 P U (du) < +∞.
Normalized multivariate Hermite polynomials as the separated polynomial chaos in G. Let a = (a 1 , . . . , a n u ) ∈ N n u be the multi-index, which includes the zero multi-index a (1) = (0, . . . , 0). For all a in N n u , let Φ a 1 ,...,a nu (u 1 , . . . , u n u ) = Φ a 1 (u 1 ) × . . . × Φ a nu (u n u ), rewritten as Φ a (u), such that Φ a (1) (u) = 1, and where Φ a i are the normalized Hermite polynomials on R. Therefore, {Φ a , a ∈ N n u } is an orthonormal family in G, ⟨Φ a , Φ b ⟩ G = δ ab , and constitutes a Hilbert basis of G.

Polynomial chaos expansion of G = g(U). For any g in G, the PCE of the second-order random variable G is written as G = a∈N nu g a Φ a (U), where the series of the PCE converges with respect to the norm of G. The real coefficients g a are such that ∥G∥ 2 G = a∈N nu (g a ) 2 < +∞ and can be calculated by

g a = ⟨G , Φ a ⟩ G = E{G Φ a (U)} = E{g(U) Φ a (U)} = R nu g(u) Φ a (u) P U (du).
Truncated PCE G chaos = g chaos (U) for representing G = g(U). Let n g be the maximum degree of the considered truncated polynomial chaos expansion. We then have |a| = a 1 + . . .

+ a n u ≤ n g . The set {a ∈ N n u , |a| ≤ n g } of all the multi-indices in N n u such that |a| ≤ n g is rewritten as {a (m) = (a (m) 1 , . . . , a (m) n u ) ∈ N n u , m = 1, .
. . , µ} with a (1) = (0, . . . , 0) and where µ = (n g + n u )!/(n g ! n u !). To simplify the notation, the polynomial chaos of multi-index a (m) is rewritten as φ m (U) = Φ a (m) (U). As in Section 3, we have,

φ 1 (U) = 1 , ⟨φ m , φ m ′ ⟩ G = δ mm ′ , E{φ m (U)} = δ 1m . (4.1)
The truncated PCE of G = g(U) is thus given by

G chaos = g chaos (U) = µ m=1 g m φ m (U) , g m = ⟨G , φ m ⟩ G , (4.2)
in which g m is the rewriting of g a (m) .

Polynomial chaos expansion of random vector Y

Hilbert space F associated with the random variable (Ξ, U). Let F = H ⊗ G be the Hilbert space defined according to the universal property of the tensor product of H and G, which has to be understood as the completion H ⊗G of space H ⊗ G. Hilbert space F is equipped with the inner product

⟨f , f⟩ F = R nw R nu f(ξ, u) f(ξ, u) P Ξ (dξ) ⊗ P U (du) , (5.1)
and the associated norm

∥ f ∥ G = ⟨f , f⟩ 1/2 F .
Multivariate polynomial chaos in F. The family of functions {Γ αa = Ψ α ⊗ Φ a , α ∈ N n w , a ∈ N n u }, in which Ψ α and Φ a are defined in Sections 3 and 4, is a Hilbert basis of F. We then have

⟨Γ αa , Γ βb ⟩ F = ⟨Ψ α , Ψ β ⟩ H × ⟨Φ a , Φ b ⟩ G = δ αβ δ ab , (5.2)
and

Γ α (1) a (1) (ξ, u) = Ψ α (1) (ξ) × Φ a (1) (u) = 1. Polynomial chaos expansion of Y = y(Ξ, U). Let F n q = L 2 P Ξ ⊗P U (R n w × R n u , R n q
) denote the Hilbert space of the square-integrable functions on R n w × R n u with values in R n q , with respect to the probability measure P Ξ ⊗ P U . Since

F n q = L 2 P Ξ ⊗P U (R n w × R n u ) ⊗ R n q , it can be deduced that F n q = F ⊗ R n q .
Considering the introduced hypotheses, the mapping (ξ, u) → y(ξ, u) defined by Eq. (2.3) belongs to F n q . Hence, the second-order R n q -valued random variable Y = y(Ξ, U) admits the polynomial chaos expansion,

Y = α∈N nw a∈N nu z αa Γ αa (Ξ, U) , (5.3)
whose coefficients z αa in R n q are such that α∈N nw a∈N nu ∥z αa ∥ 2 < +∞. The series of the PCE is convergent for the norm of F n q .

Remark on the impossibility of using projection to compute the coefficients. It may seem that coefficients z αa in R n q could be computed by the projection z αa = E{Y Γ αa (Ξ, U)}. However, it is not possible due to the fact that the random vector U, introduced as a latent random variable, is independent of Y despite the statistical dependence between Y and Ξ. Consequently, we have E{Y Ψ α (Ξ) Φ a (U)} = E{Y Ψ α (Ξ)} E{Φ a (U)} = 0 n q for all a ∈ N n u , except for a = a (1) = (0, . . . , 0).

Truncated polynomial chaos expansion of random vector Y

Using the index renumbering introduced in Eqs. (3.2) and (4.2), the truncated PCE of Y = y(Ξ, U)) is written as,

Y chaos = y chaos (Ξ, U)) = κ k=1 µ m=1 z km γ km (Ξ, U) , (6.1)
in which, for k = 1, . . . , κ and m = 1, . . . , µ, the vector-valued coefficient z km is a rewriting of z α (k) a (m) and where

γ km (Ξ, U) = Γ α (k) a (m) (Ξ, U) = ψ k (Ξ) φ m (U) , ⟨γ km , γ k ′ m ′ ⟩ F = ⟨ψ k , ψ k ′ ⟩ H × ⟨φ m , φ m ′ ⟩ G δ kk ′ δ mm ′ . ( 6
(Ξ, U) = 1, γ k1 (Ξ, U) = ψ k (Ξ), γ 1m (Ξ, U) = φ m (U)
, and Eq. ( 6.2) yields E{γ km (Ξ, U)} = δ 1k δ 1m . In addition, it can easily be seen that E{∥Y chaos ∥ 2 } = κ k=1 µ m=1 ∥z km ∥ 2 < +∞. From Eqs. (6.1) and ( 6.2), it can be deduced that

E{Y chaos } = z 11 , E{Y chaos (Y chaos ) T } = κ k=1 µ m=1 z km (z km ) T . (6.3)
As explained at the end of Section 5, the projection of Y yields z km = E{Y γ km (Ξ, U)}. However, only the vector-valued coefficients {z k,1 , k = 1, . . . , κ} can be calculated through this projection. By rewriting z k,1 as z k , these coefficients are given by z k = E{Y ψ k (Ξ)} . (6.4)

We then obtain the PCE Y chaos proj of Y through a projection without considering the latent random variable U. In other words,

Y chaos proj = κ k=1 z k ψ k (Ξ) . (6.5)
The error, ∥Y -Y chaos proj ∥ F nq , between Y = y(Ξ, U) (see Eq. (2.3)) and Y chaos proj , is significant and can only be reduced by including the latent random vector U.

Matrix representation of the realizations of the PCE Y chaos of Y. The N realizations of Y chaos , as defined by Eq. (6.1), are given by In the following, it is assumed that N ≫ J. By employing the global index j, Eq. (6.6) can be rewritten in the following matrix form,

y chaos,ℓ = κ k=1 µ m=1 z km γ km (ξ ℓ , u ℓ ) . ( 6
[y chaos ] = [z] [γ] ∈ M n q ,N , [z] ∈ M n q ,J , [γ] ∈ M J,N , (6.7)
in which the entries of matrices [y chaos ], [z], and [γ] are

[y chaos ] iℓ = y chaos,ℓ i , [z] i j = z km i , [γ] jℓ = γ km (ξ ℓ , u ℓ ) . (6.8)
Using Eqs. (6.6) to (6.8), Eq. ( 6.2) can be rewritten as

1 N -1 [γ] [γ] T = [I J ] . (6.9)
It should be noted that the factor 1 N-1 , which is used instead of 1 N , originates from the statistical estimator employed to compute the realizations of polynomial chaos (refer to Page 122 of [START_REF] Soize | Uncertainty Quantification[END_REF]).

Computation of matrix

[γ]. Since [γ] jℓ = γ km (ξ ℓ , u ℓ ) = ψ k (ξ ℓ ) φ m (u ℓ ) (see Eq. (6.2)), by introducing the matrices [ψ] ∈ M κ,N and [φ] ∈ M µ,N such that [ψ] kℓ = ψ k (ξ ℓ ) and [φ mℓ ] = φ m (u ℓ ), the entries of matrix [γ] ∈ M J,N can be written as [γ] jℓ = [ψ] kℓ [φ mℓ ] with j = (k, m). Matrix [ψ]
, which is associated with the N realizations of the nonseparated Hilbert basis constructed with P Ξ (dξ) on R n w , and matrix [φ], related to the N realizations of the separated normalized-Hermite-based Hilbert basis constructed with P U (du) on R n u , are computed using the algorithm detailed on Page 122 of [START_REF] Soize | Uncertainty Quantification[END_REF]) (see also [START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF]).

Constraint on matrix

[z] [z] T = [M Y ] , [z] ∈ M n q ,J .
(6.10)

Relationship between the matrix [y chaos ] representing the realizations of Y chaos and its counterpart [q chaos ] representing the realizations of Q chaos . Using the mapping y → q = q(y) defined by Eq. (2.2), we have for the random variables,

Q = q(Y) , Q chaos = q(Y chaos ) . (6.11)
From the N realizations of Y chaos , represented by the matrix [y chaos ] ∈ M n q ,N , we can derive the N corresponding realizations {q chaos,ℓ , ℓ = 1, . . . , N} as expressed in the matrix [q chaos ] ∈ M n q ,N . In this matrix, each entry [q chaos ] iℓ represents q chaos,ℓ i = q i (y chaos,ℓ ). This relationship can be written in matrix form as:

[q chaos ] = [Q([y chaos ])] , Q : M n q ,N → M n q ,N .
(6.12)

Optimization problem for estimating the coefficients of the truncated PCE

Optimal value [z opt ] of matrix [z]. To estimate the matrix [z] ∈ M n q ,J that contains the coefficients of the truncated PCE of Y, as defined by Eq. (6.7), various approaches can be employed, in particular the maximum likelihood method (see the references given in Section 1 regarding the statistical identification of PCE coefficients). Among all the methods, we propose using the OVL i (Overlap) indicator, which quantifies the overlap between the PDF of the components Q i of Q and the PDF of the components

Q chaos i of Q chaos . The OVL i indicator is associated with the L 1 -norm of functions q i → p Q i (q i ) -p Q chaos i (q i ; [z]) on R, where p Q i and p Q chaos i (• ; [z]
) are the probability density functions of the real-valued random variables Q i and Q chaos i , respectively. These PDFs are estimated using Gaussian KDE applied to the realizations {q ℓ i , ℓ = 1, . . . , N} for Q i and to the realizations {q chaos,ℓ i , ℓ = 1, . . . , N} for Q chaos i , which correspond to the columns of matrix [q chaos ] = [Q([z] [γ])] (see Eq. (6.12) with Eq. (6.7)). For i ∈ {1, . . . , n q }, OVL i ([z]) is written as

OVL i ([z]) = 1 - 1 2 R |p Q i (q i ) -p Q chaos i (q i ; [z])| dq i , (7.1)
and the cost function is defined by

J([z]) = 1 n q n q i=1 OVL i ([z]) . (7.
2)

It can be seen that 0 ≤ J([z]) ≤ 1 and the upper bound is reached when

p Q i = p Q chaos i (• ; [z]
) for all i. Hence, the optimization problem is written as

[z opt ] = arg max [z]∈C ad J([z]) , (7.3)
in which the admissible set C ad ⊂ M n q ,J allows the constraint defined by Eq. (6.10) to be taken into account,

C ad = { [z] ∈ M n q ,J , [z] [z] T = [M Y ] } . (7.4)
Transforming the constrained optimization problem into an unconstrained optimization problem. The optimization problem defined by Eq. ( 7.3) is nonconvex. We propose to solve it using an algorithm designed for unconstrained optimization problem. Due to the nonconvex nature of the problem, the estimated solution will strongly depend on the initial point [z 0 ] chosen for the initialization of the optimization algorithm. We then need to transform the optimization problem on C ad into an unconstrained optimization problem on M n q ,J , and also carefully choose the initial point. To do this, we need to introduce a [ẑ] → [z] transformation to eliminate the constraint, which will then be automatically satisfied. Additionally, we have to introduce a second transformation [z] → [ẑ] to search for an optimal solution in the vicinity of [z 0 ] = [z] ∈ M n q ,J . The columns of [z], denoted z 1 , . . . , z J , are defined as the projection on the polynomial chaos (see Section 6).

(a) Transformation [ẑ] → [z] from M n q ,J into M n q ,J . Let [ẑ] be any unconstrained matrix given in M n q ,J . Let [c] in M n q be the upper triangular matrix resulting from the Cholesky factorization of the matrix

[ẑ] [ẑ] T ∈ M + n q . We have [c] T [c] = [ẑ] [ẑ] T , which implies the existence of [c] -1 . Utilizing the decomposition [M Y ] = [L Y ] 2 introduced in Section 2, the desired transformation is expressed as [ẑ] → [z] = [L Y ] [c] -T [ẑ] : M n q ,J → M n q ,J . (7.5)
It can easily be verified that for any [ẑ] in M n q ,J , we have

[z] [z] T = [M Y ]. (b) Transformation [z] → [ẑ] from M n q ,J into M n q ,J . Matrix [z]
being constructed as the projection of [y] onto the subspace spanned by [γ], using Eqs. (6.7) and (6.9) yields

[z] = 1 N -1 [y] [γ] T . (7.6) It is important to note that [z] [z] T = [y] [χ] [y] T where [χ] = [γ] T [γ]/(N -1) [I N ]. Consequently, [z] [z] T [M Y ].
For [z] ∈ M n q ,J defined by Eq. (7.6), the transformation [z] → [ẑ] from M n q ,J into M n q ,J is defined by

[ẑ] i j = [z] i j (1 + [z] i j ) , i ∈ {1, . . . , n q } , j ∈ {1, . . . , J} . (7.7)
This transformation shows that as [z] explores M n q ,J in the vicinity of [z 0 ] = [0], [ẑ] also explores M n q ,J in the vicinity of [z], and (d) Reformulation of the optimization problem and algorithm. The optimization problem defined by Eq. ( 7.3) can be rewritten as,

[z] = [L Y ] [c] -T [ẑ] satisfies the constraint [z] [z] T = [M Y ]. (c) Transformation [z] → [z] = [z([z])] from M n q ,J into M n q ,J .
[z opt ] = [z([z opt ])] , [z opt ] = arg max [z]∈M nq,J J([z]) , J([z]) = J([z([z])]) . (7.8)
This nonconvex optimization problem can be solved using various algorithms. As we are searching for a solution [z opt ] in the vicinity of [0 n q ,J ], an appropriate choice for the algorithm may be the unconstrained quasi-Newton algorithm, initialized with [z 0 ] = [0 n q ,J ], where the gradient is not explicitly provided.

Optimal truncated PCE of Y. The Gaussian probability measure P U (du) of the latent random variable U depends on its dimension n u and the polynomial chaos expansion in U depends on the maximum degree n g (see Section 4).

Consequently, the PCE Y chaos = y chaos (Ξ, U) of Y = y(Ξ, U), which depends on N g , also depends on n u and n g . The presented methodology allows for estimating an optimal value [z opt (N g , n u , n g )] that depends on N g , n u , and n g . For a given value of each of these three integers, the error between Y and its chaos representation Y chaos can be quantified by evaluating

J opt (N g , n u , n g ) = J([z opt (N g , n u , n g )]) , (7.9)
in which J([z]) is defined by Eq. (7.2). Consequently, the optimal truncated PCE of Y is obtained by using the optimal values N opt g , n opt u , and n opt g of N g , n u , and n g , respectively, such that (N opt g , n opt u , n opt g ) = arg max

N g ≥2 , n u ≥1 , n g ≥1 J opt (N g , n u , n g ) , (7.10)
Comments about the algorithm for solving the optimal values N opt g , n opt u , and n opt g . These optimal values are obtained by solving the optimization problem defined by Eq. (7.10). A simple and direct approach, although it can be computationally expensive, is to compute the value of J opt (N g , n u , n g ) at each point in a three-dimensional grid. This grid corresponds to a discretization of the domain [2,

N max g ] × [1, n max u ] × [1, n max g ] ⊂ N 3 , in which N max g , n max
u , and n max g are set to sufficiently large values. Another approach, which is less computationally expensive, is to use the assumptions introduced in section 1, which are related to the underlying physical problem that generated the data. In this context, if it were not necessary to introduce the latent variable U, then the solution would be [z opt (N g )] = [z(N g )] given by Eq. (7.6). The corresponding value of the overlap can be expressed as follows,

J(N g ) = J([z(N g )]) . (7.11)
Since the latent variable U is essential, at convergence with respect to N g , the overlap J(N g ) will be less than 1.

The difference 1 -J(N g ) makes it possible to quantify the error induced by the projection method. We can then quickly estimate an optimal value N opt g for N g using a one-dimensional grid over [2, N max g ] ⊂ N. In the presence of the latent variable U, we set N g to this optimal value N opt g and then we search for the optimal values n opt u and n opt g on a two-dimensional grid over [1,

n max u ] × [1, n max g ] ⊂ N 2 .

Polynomial-chaos-based statistical surrogate model

Presenting the problem to be solved. The problem at hand involves computing the realization q chaos 0 of Q chaos given Ξ = ξ 0 , where ξ 0 represents a realization of Ξ following the probability measure P Ξ (dξ), and u 0 represents a realization of U following the probability measure is P U (du). In fact, a realization w 0 of W is given, then the realization ξ 0 is calculated by ξ 0 = s -1 w (w 0 ) (see Eq. (2.1)).

Polynomial chaos-based algorithm for computing the conditional realization q chaos 0 of Q chaos given Ξ = ξ 0 . Using the optimal truncated PCE Q chaos of Q, we compute the corresponding realization q chaos 0 of Q chaos according to Eqs. (6.11), (6.7), and (6.8), by

q chaos 0 = q(y chaos 0 ) , y chaos 0 = [z opt (N opt g , n opt u , n opt g )] γ opt 0 , γ opt 0 ∈ R J . (8.1)
Validation. The verification of Eq. (8.1) is performed in a probability framework as follows. We generate N v ∼ N additional realizations ξ 0 of Ξ and u 0 of U, which are distinct from the realizations {(ξ ℓ , u ℓ ), ℓ = 1, . . . , N}. For each realization (ξ 0 , u 0 ), we compute the realization q chaos 0 of Q chaos using Eq. (8.1). We quantify the error between the reference Q (defined by the N realizations from the learned dataset) and Q chaos (defined by the N v realizations generated using Eq. (8.1)), by computing J opt using Eq. (7.9) for the optimal values N opt g , n opt u , and n opt g . Additionally, we compare the probability density functions of Q i and Q chaos i for i = 1, . . . , n q .

Cluster separation of the learned dataset constituted of heterogeneous data

As we explained in Section 1, in the case of heterogeneous data, it can be interesting to use an adapted approach based on the formulation presented in Sections 3 to 8. This approach involves performing a prior cluster separation of the learned dataset into distinct clusters, each consisting of "quasi-homogeneous" data. It is important to note that the objective remains the same as before, which is to construct a global polynomial chaos representation

Q chaos = f chaos (W, U) of Q = f(W, U
) based on the polynomial chaos representation constructed for each distinct cluster. This is particularly useful when the training dataset is generated by physical processes that exhibit multiple regimes simultaneously. Although the numerical offline computation may increase, the online prediction with the PCE Q chaos given Ξ = ξ 0 remains unaffected. It should be noted that the cluster separation is performed on the learned dataset rather than the training dataset, which is too small. This choice offers an advantage because the PLoM algorithm used to generate the learned dataset from the training dataset preserves the concentration of the probability measure in the vicinity of the manifolds while enhancing the available information contained in the points of the training dataset. The use of prior separation into distinct clusters, allows us to better understand how the learned dataset is structured and facilitates the construction of the representation in polynomial chaos for each cluster.

Hypothesis used for building the cluster separation of the heterogeneous learned dataset. It is assumed that there are multiple physical regimes, which may involve a cluster separation of the learned dataset D learn (X) into K distinct clusters based on their proximity or similarity.

Comments on existing algorithms for constructing cluster separation into distinct clusters of points. There are various algorithms available for constructing a cluster separation of a set of points. The principal methods include K-means [START_REF] Hartigan | A K-means clustering algorithm[END_REF][START_REF] Sinaga | Unsupervised K-means clustering algorithm[END_REF], Hierachical Agglomerative Clusturing [START_REF] Lukasová | Hierarchical agglomerative clustering procedure[END_REF], Divisive Clusturing [START_REF] Savaresi | Cluster selection in divisive clustering algorithms[END_REF], Density-Based Spectral Clusturing with its variants [START_REF] Ester | Density-based spatial clustering of applications with noise[END_REF][START_REF] Sander | Density-based clustering in spatial databases: The algorithm GDBSCAN and its applications[END_REF][START_REF] Bäcklund | A density-based spatial clustering of application with noise[END_REF], Mean-Shift Clusturing [START_REF] Wu | Mean shift-based clustering[END_REF], Spectral Clusturing [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], Gaussian Mixture Model [START_REF] Maugis | Variable selection for clustering with gaussian mixture models[END_REF], and more. In the context of the previously introduced hypothesis, our objective is to determine the number K of clusters while considering the existence of K physical regimes. Therefore, algorithms based on Density-Based Spectral Clustering are not well suited for this purpose. For large datasets, algorithms such as Agglomerative Clustering, Mean-Shift Clustering, or Spectral Clustering (which is well adapted to image segmentation) require significant computational resources in terms of RAM and CPU. The Gaussian mixture model could be a candidate approach. However, this method is sensitive to dimensionality and outliers, and gives for each identified cluster a Gaussian measure, which may not be suitable for the non-Gaussian case. For the application presented in this paper, three clustering algorithms have been tested: Divisive Clustering, Hierarchical Agglomerative Clustering, and K-means with the squared Euclidean distance. In K-means, each centroid is calculated as the mean of the points in its respective cluster.

Notation and formulation related to the cluster separation. The chosen algorithm is applied to perform the partition into K clusters of D learn (Q) = {q ℓ , ℓ = 1, . . . , N}. For σ = 1, . . . , K, let I σ = {ℓ (σ) 1 , . . . , ℓ (σ) N σ } ⊂ {1, . . . , N} be the subset of indices resulting from this cluster separation. We then have N = K σ=1 N σ , and

∪ K σ=1 I σ = {1, . . . , N} , ∩ K σ=1 I σ = {∅} . (9.1)
Consequently, the resulting partition of D learn (X) can be expressed as,

D learn (X) = ∪ K σ=1 D σ (X (σ) ) , ∩ K σ=1 D σ (X (σ) ) = {∅} , D σ (X (σ) ) = {x ℓ = (q ℓ , w ℓ ) , ℓ ∈ I σ } , (9.2)
where

X (σ) = (Q (σ) , W (σ)
) is the R n x -valued random variable for which {x ℓ , ℓ ∈ I σ } are N σ independent realizations. Furthermore, we introduce a discrete-valued random variable B d with values in {1, 2, . . . , K}, which is statistically dependent on W. The realizations {b ℓ d , ℓ = 1, . . . , N} of B d , corresponding to {w ℓ , ℓ = 1, . . . , N}, are defined as follows,

if

ℓ ∈ I σ then b ℓ d = 2σ -(K + 1) , σ ∈ {1, . . . , K} . (9.3)
For instance, for

K = 2, if ℓ ∈ I 1 , then b ℓ d = -1 while if ℓ ∈ I 2 , then b ℓ d = +1.
It is important to note that it is impossible to obtain a "perfect" separation, regardless of the algorithm used within the framework of probability theory, because of the potential partial overlap of the supports of the probability measures for the random variables X (1) , . . . , X (K) . For a given realization w 0 of W and u 0 of U, we need to identify the σ-index in {1, . . . , K} of the associated cluster, as explained earlier. Considering the generation of realizations {b ℓ d , ℓ = 1, . . . , N} defined by Eq. ( 9.3), along with the construction of the smoothed random variable B, which will be represented by B chaos , we propose the following method for identifying the σ-index. Let B σ be the subset of R defined by

B 1 =] -∞, 2 -K], B K =]K -2, +∞[, and if K ≥ 3, B σ =]2σ -K -2, 2σ -K] for 2 ≤ σ < K. Let b 0 = b(
w 0 , u 0 ) be the corresponding realization of B, that will be approximated by b chaos 0 = b chaos (w 0 , u 0 ). We identify the value of the σ-index in {1, . . . , K} as the index σ of interval B σ to which b 0 belongs.

Truncated PCE of B and of {X (σ) , σ = 1, . . . , K}, computation of the PCE-based statistical surrogate model, realizations, and validation. The truncated PCE of B and X (σ) for all σ = 1, . . . , K are constructed using the methodology outlined in Sections 3 to 8. Following the approach described in Section 8, we consider any realization ξ 0 of Ξ from P Ξ (dξ) and any realization u 0 of U from P U (du). By utilizing the surrogate model based on cluster separation, we compute the conditional realization q 0 of Q chaos given Ξ = ξ 0 . The validation process follows the methodology presented in Section 8.

Application to atomic collisions of Helium on graphite substrate

Description of the physical system and its Molecular Dynamics simulations for generating the training dataset

The details of the dynamical system and its analysis can be found in [START_REF] Magnico | Collisions, adsorption and self diffusion of gas in nanometric channels by molecular dynamics and stochastic simulation and the case of helium gas in graphitic slit pore[END_REF]. Here we provide a brief overview. The wall model consists of three layers of graphene with dimensions 17.04 × 17.22 Å 2 in the xOy plane and is periodically replicated in the x and y directions. The z coordinates are orthogonal to the xOy plane. The wall contains 336 carbon atoms and has a width of 6.8, AA. All molecular dynamics simulations involving graphite, which is composed of carbon atoms (C), and helium (He) gas atoms, were performed using the LAMMPS (Large-scale Atomic Molecular Massively Parallel Simulator) package. The adaptive intermolecular reactive bond order potential was used to model the interactions between the carbon atoms in graphite. The 12-6 Lennard Jones potential was employed to describe the interactions between carbon and helium atoms. A cutoff distance r cut of 12 Å has been chosen for these interactions. A graphitic wall served as the lower boundary, and a reflective plane was located at a distance of 18.8 Å from the surface, serving as the upper boundary. A control plane was placed at a distance of z = r cut from the carbon surface, where information about the crossing atoms was recorded, including their incoming/outgoing times, velocities, and displacements during their residence time. During the simulation, the bottom sheet was fixed, while the other two layers were kept at a constant temperature of 50 • K using the Nosé-Hoover thermostat in the NVT ensemble, with a relaxation temperature parameter set to 40 time steps. The layers were free to interact with the helium atoms. Approximately 1 × 10 6 time steps of 1 f s were spent equilibrating the system at the desired temperature before proceeding with subsequent statistical analysis. The velocity-Verlet algorithm was used for time integration.

In the following, the three components of the control parameter W correspond to the components (V x,in , V y,in , V z,in ) of the velocity V in of the incident particle on the layer. The random output vector Q = (V out , log(∆ t ), D x , D y ) consists of the random output velocity vector V out = (V x,out , V y,out , V z,out ) (reflected velocity) of the particle, the random logarithm duration of absorption log(∆ t ) (residence time), and the two displacement components D x and D y in the xOy plane within the layer. In this application, all distances are measured in Angstroms (Å), and time is measured in picoseconds (ps). Figure 1: For N = 200 000 learned realizations, PDF of the components W 1 , W 2 , and W 3 of the random control parameter W, which represents the incident velocity vector V in = (V x,in , V y,in , V z,in ). (f) PDF q 6 → p Q 6 (q 6 ) Figure 2: For N = 200 000 learned realizations, PDF of the components Q 1 to Q 6 of the random output vector Q, which represents (V x,out , V y,out , V z,out , log(∆ t ), D x , D y ).

plays the PDF of the logarithm of the residence time (duration of absorption) Q 4 = log(∆ t ), reveals the existence of two distinct physical regimes. The first regime corresponds to a short residence time, with a PDF peak at approximatively exp(1.38) ≃ 4 ps. Conversely, the peak of the long residence-time regime is around exp(4.5) ≃ 90 ps. It should be noted that these two residence-time regimes are not entirely separated, as the PDF exhibits a local minimum at approximately exp(3) ≃ 20 ps, where the PDF has a significant value of 0.11. The presence of these two physical regimes is also evident in the PDF of Q 5 = D x (Figure 2e) and Q 6 = D y (Figure 2f). These PDFs exhibit strong non-Gaussian behavior and display two distinct patterns: one for the distance in the range of [-120, 120], Å, and the other for distances in the intervals ] -∞, -120[∪]120, +∞[ Å. For a more comprehensive analysis of the underlying physics, we refer the reader to [START_REF] Magnico | Collisions, adsorption and self diffusion of gas in nanometric channels by molecular dynamics and stochastic simulation and the case of helium gas in graphitic slit pore[END_REF].

PCE constructed through projection and quantification of the induced error. As explained in Section 7, the projection-based truncated PCE of Y is constructed in estimating the coefficients-matrix [z] in Eq. (6.7) by using the projection method yielding [z] = [z] defined by Eq. (7.6). Figure 3 displays the graph of function N g → J(N g ) defined by Eq. (7.11). Since, with a such projection, the random latent variable U is not taken into account (corresponding to (n u , n g ) = (0, 0)), the projection method induces a significant error of 1 -0.7425 for N g = 12, which is approximately a 25% error. We intentionally limited the degree N g to 12, yielding κ = 455 polynomial chaos. In fact, increasing the degree leads to a higher computational cost with little improvement in convergence. Introducing the latent variable U is necessary to reduce the error. Figures 4a to 4f show the PDFs of the 6 components of Q. These PDFs correspond to the reference values (estimated using the learned dataset) and the values estimated using the truncated PCE with the projection method and N g = 12. All computation are performed using N = 200 000 learned realizations. As expected, there is a significant error between the reference values and the PCE constructed through projection. This comparison is interesting because it demonstrates that the presence of the two physical regimes can only be reproduced with the PCE by introducing the latent variable U. The control variable W, which represents the incident velocity vector alone, insufficient to explain the phenomena. (f) PDF q 6 → p Q 6 (q 6 ) Figure 4: For N = 200 000 learned realizations, PDF of the components Q 1 to Q 6 of the random output Q, corresponding to the reference (blue thick line) and estimated with the truncated PCE by the projection method with N g = 12 (black thin line).

Optimization problem for estimating the coefficients of the truncated PCE and convergence analysis with respect to n u and n g . For a total of 200,000 learned realizations, the convergence with respect to N g , n u , and n g , has been analyzed by studying the function (N g , n u , n g ) → J opt (N g , n u , n g ) defined by Eq. (7.9). As previously explained, N d = 12 is a suitable value for N opt g , which will be confirmed below. We have studied the function (n u , n g ) → J opt (12, n u , n g ) for N g = 12. The calculation gives J opt (12, 0, 0) = 0.7425, and for the points (1, 1), (1, 2), and (2, 2), we have got 0.996. This indicates that convergence is reached when n u = 1 and n g = 1. However, looking at the PDFs of the components of Q chaos , we see a slightly better match to the PDF of Q when n u = 2 and n g = 2, corresponding to µ = 6 polynomial chaos. Figure 5a displays the graph of the function N g → J opt (N g , 2, 2). This figure shows that N opt g = 12 is an optimal choice, and gives excellent convergence. As previously mentioned, the optimization problem defined by Eq. (7.8) is solved using the quasi-Newton algorithm. For N g = 12, n u = 2, and n g = 2, yielding J = κ × µ = 455 × 6 = 2 730 coefficients with values in R 6 for the truncated PCE, Figure 5b displays the graph of the function ι → J([z ι ]), where ι represents the iteration number in the quasi-Newton algorithm. This graph illustrates the rapid convergence rate. Validation of the polynomial-chaos based conditional statistics. Once the optimal coefficients of the truncated PCE have been estimated as described above, the validation process is carried out following the explanation in Section 8.To do so, N v = 200 000 new realizations of random control variable W are generated according to the probability measure P W (dw). Then, N v corresponding realizations of Q chaos are generated using Eq. (8.1). Figures 6a to 6f display the PDFs of the 6 components of Q. These PDFs correspond to the reference values (estimated with the learned dataset) and the values estimated using the optimal truncated PCE with N g = 12, n u = 2, and n g = 2. The predictions obtained with this statistical substitution model based on the truncated PCE are very good.

Polynomial-chaos-based surrogate model using a cluster separation

As we have explained in Section 10.2 (also see Figure 2d, the heterogeneous data is generated by two distinct physical regimes that cannot be perfectly separated. To identify a good cluster-separation algorithm adapted to the considered heterogeneous data, we tested three algorithms: K-means, Divisive Clusturing, and Hierachical Agglomerative Clusturing, with K = 2 and K = 3 clusters. The best separation was achieved with K-means for K = 2 clusters. All computations were performed using the N = 200 000 points of D learn (Q). Using the notation introduced in Eqs. (9.1) to (9.3), we obtained a first cluster σ = 1 with I 1 consisting of N 1 = 21 059 points, and the second cluster σ = 2 with I 2 comprising N 2 = 178 941 points. Cluster D 1 (X (1) ) corresponds to a physical regime with a long residence time, while cluster D 2 (X (2) ) corresponds to a mixture of the two physical regimes, primarily containing short residence time but also medium and long residence time. Certainly, a better separation could be obtained "manually", based on physics expertise without using a clustering algorithm and based on the analysis of the single component Q 4 (see [START_REF] Magnico | Collisions, adsorption and self diffusion of gas in nanometric channels by molecular dynamics and stochastic simulation and the case of helium gas in graphitic slit pore[END_REF] in the context of this application). However, we wish to present here a cluster separation using a separation algorithm which acts "simultaneously" on all the components of the random vector Q (because these components are statistically dependent). The goal is to illustrate the algorithm we propose to construct realizations using a statistical surrogate model constructed from the PCE of each cluster and the PCE of the random variable B. In all the figures of Section 10.3 displaying PDFs, the thin blue lines correspond to the reference estimated by the Gaussian KDE using N 1 = 21 059 realizations for Q (1) , N 2 = 178 941 realizations for Q (2) , N = 200 000 realizations for both B and Q. (f) PDF q 6 → p Q 6 (q 6 ) Figure 6: PDF of the components Q 1 to Q 6 of the random output vector Q, corresponding to the reference (blue thin line) and estimated with the optimal truncated PCE with N g = 12, n u = 2, and n g = 2 (red thick line). In figures (a), (b), and (c), the curves are almost identical.

These realizations are obtained from the learned dataset.

Optimal PCE of Q (1) corresponding to cluster σ = 1. Convergence is reached for N opt g = 4, n opt u = 2, and n opt g = 2. The corresponding value of J opt (N opt g , n opt u , n opt g ) is 0.990, indicating a very good convergence. Figures 7a to 7f display the PDFs of the 6 components of Q (1) (the reference) compared to the estimated PDFs of Q (1),chaos using the optimal truncated PCE with N 1 points. The predictions obtained through this approach for the first cluster are accurate. Figure 7d shows that this cluster, D 1 (X (1) ), corresponds to a long residence time, as the peak of the PDF of Q (1) Optimal PCE of Q (2) corresponding to cluster σ = 2. Convergence is reached for N opt g = 8, n opt u = 2, and n opt g = 2. The corresponding value of J opt (N opt g , n opt u , n opt g ) is 0.987, indicating a good convergence. Figures 8a to 8f display the PDFs of the 6 components of Q (2) (the reference) compared to the estimated PDFs of Q (2),chaos using the optimal truncated PCE with N 2 points. The predictions obtained through this approach for the second cluster are accurate. Figure 8d shows a mixture of the two physical regimes in this cluster, D 2 (X (2) ), with the dominant regime of short residence time, as the first peak of the PDF of Q (2) 4 occurs at ∆ t = exp(1.33) ≃ 3.8 ps.

Optimal PCE of random variable B for identifying the cluster number σ. Convergence is reached for N opt g = 14, n opt u = 2, and n opt g = 2. The corresponding value of J opt (N opt g , n opt u , n opt g ) is 0.900. It should be noted that the rate of convergence is not critical as we are dealing with two clusters (K = 2), where B 1 =] -∞, 0] and B 2 =]0, +∞[. Thus, we only need to determine if, for each realization b chaos 0 = b chaos (w 0 , u 0 ) of B chaos , we have b chaos 0 ≤ 0 or b chaos 0 > 0. Figure 9 displays the PDF of B (the reference) compared to the estimated PDF of B chaos using the optimal truncated PCE with N points. The prediction obtained through this approach for the second cluster is sufficiently accurate.

Validation of the PCE-based statistical surrogate model using the cluster separation. Once the optimal coefficients of the truncated PCEs for Q (1),chaos , Q (2),chaos , and B chaos have been estimated as described above, the validation process is carried out following the explanation in Section 9. To do so, N v = 200 000 new realizations of random control (f) PDF q 6 → p Q

(q 6 ) Figure 7: PDF of the components Q (1) 1 to Q (1) 6 of the random output vector Q (1) , corresponding to the reference (blue thin line) and Q (1,chaos) 1 to Q (1,chaos) 6 of the random output vector Q (1,chaos) , estimated with the optimal truncated PCE with N g = 4, n u = 2, and n g = 2 (red thick line). (q 6 ) Figure 8: PDF of the components Q (2) 1 to Q (2) 6 of the random output vector Q (2) , corresponding to the reference (blue thin line) and Q (2,chaos) Figure 9: Graph of the function N g → J(N g ) representing the relationship between the degree N g of the PCE constructed using the projection method, enabling quantification of the induced error.

1 to Q (2,chaos)
variable W are generated according to the probability measure P W (dw). Then, N v corresponding realizations of Q chaos are generated. Figures 10a to 10f display the PDFs of the 6 components of Q. These PDFs correspond to the reference values (estimated with the learned dataset) and the values estimated using the algorithm based on the three PCEs. The predictions obtained through this PCE approach are highly accurate. (f) PDF q 6 → p Q 6 (q 6 ) Figure 10: PDF of the components Q 1 to Q 6 of the random output vector Q, corresponding to the reference (blue thin line) and estimated on the base of the optimal truncated PCEs of the clusters (red thick line).

Conclusion

We have presented a formulation and an algorithm developed a polynomial chaos of a vectorvalued random quantity of interest (the output) as a function of an input composed of a part of a vector-valued random control parameter with known probability measure and another part of a vector-valued random latent variable with unknown probability measure. The training dataset consists of heterogeneous data, which poses challenges for accurately estimating the chaos coefficients. The proposed approach enables the construction of a global, accurate, and highly efficient statistical surrogate model for evaluating an output realization given an input realization. As an alternative, we have also proposed a construction based on the clustering of the learned dataset, which can facilitate offline construction in the case of heterogeneous data. The application presented, which concerns the collisions of helium atoms on a graphite substrate, demonstrates the accuracy and efficiency of the proposed approach. This approach, which allows the construction of a rapid online surrogate model, is general and should allow the analysis of large complex systems beyond the computational capabilities currently available.

  M n : set of the square (n × n) real matrices. M + n : set of the positive-definite (n × n) real matrices. M +0 n : set of the positive (n × n) real matrices. N: number of points in the learned dataset. n d : number of points in the training dataset. N: set of all the integers {0, 1, 2, . . .}. R: set of all the real number.

R

  n : Euclidean space of dimension n. [I n ]: identity matrix in M n . x = (x 1 , . . . , x n ): point in R n . [x] T : transpose of matrix [x]. tr{[x]}: trace of the square matrix [x]. α: multi-index α = (α 1 , . . . , α n ) in N n . δ 0 : Dirac measure on R at point 0. δ kk ′ : Kronecker's symbol. δ αβ : δ α 1 β 1 × . . . × δ α n β n . E: mathematical expectation operator. KDE: kernel density estimation. PDF: probability density function. PCE: polynomial chaos expansion.

. 6 )

 6 Instead of using indices k and m, we introduce the global index j such that j = (k, m) ∈ {1, . . . , J} for (k, m) ∈ {1, . . . , κ} × {1, . . . , µ} with J = κ × µ .

  [z] defined by the second-order moment-matrix[M Y ] of Y. Similarly to the introduced notation [y chaos ], we define [y] ∈ M n q ,N such that [y] iℓ = y ℓ i . The estimate of [M Y ] = E{YY T } is then expressed using the same notation: [M Y ] = [y] [y] T /(N -1) (note that [M Y ]is a given data derived from the learned dataset, as explained in Section 2). Hence, by imposing the equation E{Y chaos (Y chaos ) T } = [M Y ], we derive the following constraint,

  The composition of transformation [z] → [ẑ] defined by Eq. (7.5) and transformation [ẑ] → [z] defined by Eq. (7.7) is a well-defined transformation [z] → [z] = [z([z])].

1 N N ℓ=1 1 √

 11 Smoothing the discrete random variable B d into a real-valued random variable B. In order to use the same algorithm as the one presented in Sections 3 to 8 for constructing the truncated PCE B chaos = b chaos (W, U) of B d = b d (W, U), we introduce a real-valued random variable B that corresponds to smoothing B d on R. The probability density function of B = b(W, U) is constructed using the Gaussian KDE with the realizations {b ℓ d , ℓ = 1, . . . , N} defined by Eq. (9.3). It can be expressed as p B (b ; β) = 2π β exp{-1 2β 2 (bb ℓ d ) 2 } where β represents the bandwidth. As β approaches zero, p B (b, ; β), db converges to the discrete probability measure P B d (db) = 1 N N ℓ=1 δ 0 (bb ℓ d ). The smoothing of B d with B involves reducing the Silverman bandwidth, which is accomplished by selecting β = 1 2 ( 4 3N ) 1/5 σ B d , where σ B d is the standard deviation of B d estimated using the realizations {b ℓ d , ℓ = 1, . . . , N}. Identifying the cluster number σ based on a given realization of the real-valued random variable B. The presented formulation for cluster separation enables the construction of the truncated PCE, B chaos = b chaos (W, U), of B = b(W, U).

10. 2 .

 2 Polynomial-chaos-based surrogate modelGeneration of the large learned dataset and probability density functions of W and Q. For the application under consideration, all the presented results correspond to a temperature of 50 • K. The large learned dataset is generated as explained in Section 2, with N = 200 000 learned realizations on the basis of the training dataset consisting of n d = 2 000 points. The training dataset was generated as explained in Section 10.1. The probability density functions of the 3 components of W, estimated using the Gaussian KDE with the N learned realizations, are shown in Figures1a to 1c. The corresponding probability density functions of the 6 components of Q are shown in Figures2a to 2f. These figures define the reference and will be used for comparing with the PDFs given by the PCE. Figure2d, which dis-PDF w 3 → p W 3 (w 3 )

4 (

 4 d) PDF q 4 → p Q 4 (q 4 ) PDF q 5 → p Q 5 (q 5 )

Figure 3 :

 3 Figure3: Graph of the function N g → J(N g ) allowing the error induced by the projection method to be quantified.

3 ( 3 (

 33 PDF q 1 → p Q 1 (q 1 ) b) PDF q 2 → p Q 2 (q 2 ) c) PDF q 3 → p Q 3 (q 3 )

  Graph of function N g → J opt (N g , 2, 2).

  Graph of function ι → J([z ι ]).

Figure 5 :

 5 Figure 5: Convergence analysis of the truncated PCE with respect to N g for n u = n g = 2 (a) and convergence of the quasi-Newton algorithm as a function of the iteration number ι (b).

4 (

 4 d) PDF q 4 → p Q 4 (q 4 ) PDF q 5 → p Q 5 (q 5 )

4 occurs

 4 at ∆ t = exp(4.8) ≃ 120 ps. The peak of the PDF of Q 5 = D x is reached at D x = -123 Å, and the peak of the PDF of Q 6 = D y is reached at D y = 190 Å.

of the random output vector Q(2,chaos) , estimated with the optimal truncated PCE with N g = 8, n u = 2, and n g = 2 (red thick line). In figures (b) and (c), the curves are almost identical.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work the authors do not use Generative AI and AI-assisted technologies in the writing process and take full responsibility for the content of the publication.