domain: spi phys

text:
Nicolas Vauchel
email: nicolas.vauchel@onera.fr

Éric Garnier

Thomas Gomez

A multi-element non-intrusive Polynomial Chaos method using agglomerative clustering based on the derivatives to study irregular and discontinuous Quantities of Interest

Keywords: Uncertainty Quantification, Multi-element Polynomial Chaos, Machine learning

A non-intrusive method to get a multi-element Polynomial Chaos model is developed. This method is called ME-ACD, for Multi-Element based on Agglomerative Clustering on Derivatives. It aims at approximating a Quantity of Interest which presents discontinuities or irregularities making it difficult to be accurately approximated by standard Polynomial Chaos models. The method permits to efficiently split the parameter space and to train local polynomial models of lower degrees on every element where the local pieces of the Quantity of Interest are smoother. The algorithm is based on agglomerative clustering of the observations in a well-chosen abstract space taking into account the value of the Quantity of Interest and of its derivatives with respect to the stochastic input parameters. The same observations are used for both partitioning the space and training the local models. Several partitions of the parameter space are tested, and the one leading to local models minimizing a crossvalidation error is selected. Once the training observations are labelled with a class number indicating the element they are located in, a neural network classifier is trained to determine which local model to use for further evaluations. The method has proven to efficiently split the parameter space for a set of applications of moderate dimension. The piecewise chaos model is compared with the ones of a standard Polynomial Chaos non-intrusive method and of a Gradient Boosted Trees method in terms of accuracy.

Introduction

Uncertainty Quantification (UQ) has become paramount in physical simulations in a wide range of fields. Deterministic simulations, with fixed parameters, have been the core of numerical analysis during many years, but they do not take into account the potential uncertainties of the parameters, notably the ones coming from data. UQ focuses on the impact of these uncertainties on an output, which is often denoted the "output" or the "Quantity of Interest" (QOI). The studied uncertainties are called the "input parameters" and are supposed to follow a stochastic distribution law. The determination of the distribution law of an input parameter is one step of UQ which is not discussed in the context of this article. Here is assumed that the input parameters are independent with a known probability distribution function (pdf). Some transforms are known in the literature to handle the case of independent input parameters and map them to independent ones [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF][START_REF] Nataf | Determination des distribution don't les marges sont donnees[END_REF]. The classical Monte-Carlo (MC) method is known to be robust and simple to implement, but has a low convergence rate and only gives information on the statistic moments of the QOI. The use of surrogate models overcomes these drawbacks. Surrogate models approximate the relation between the QOI and the input parameters by an analytical formula. In this article, we are interested in the Polynomial Chaos expansion (PC), first introduced by Wiener in [START_REF] Wiener | The homogeneous chaos[END_REF]. Originally, this expansion aimed at building an approximation of a stochastic process as a series of Hermite polynomials of gaussian random variables. PC was later extended in [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] for other stochastic distributions using other orhogonal families, introducing the gPC (generalized polynomial chaos) framework. In the article, PC denotes an expansion in the gPC framework. Constructing a PC model for an Ordinary Differential Equation or a Partial Differential Equation is classically done by Stochastic Galerkin (SG) methods, also called "Intrusive methods" [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] in the field of chaos polynomials. SG methods extend both input parameters and QOI into a chaos series and use orthogonalisation to obtain a system of coupled differential equations with every coefficient of the polynomial approximation series of the QOI. These methods can therefore result in large systems and need to develop new simulation codes, which is extremely difficult for complex systems. In the context of this article, "Non-intrusive" (NI) methods are considered [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF][START_REF] Xiu | Efficient collocational approach for parametric uncertainty analysis[END_REF][START_REF] Lin | Stochastic computational fluid mechanics[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Guo | Constructing least-squares polynomial approximations[END_REF]. NI methods consist in fixing a value of the input parameters and then running the deterministic associated system, which is suitable for QOIs with complex computation, seeing it as a "black box". NI methods are then close to Machine learning (ML) problems. A set of fixed values of the input parameters and the corresponding value of the QOI with these parameters is called an "observation" (or in some references a "node"). Observations can be represented in an abstract space. The parameter space is one possible abstract space where every coordinate of a point corresponds to the value taken by the input parameters. The number of input parameters is then representing the dimension of this space. The PC model is obtained with the values of the QOI on a set of observations, called "training observations". There are numerous NI methods. Some of them are based on classical interpolation or linear regression methods but with observations representing stochastic uncertainties. For more details on "Intrusive" or "Non-intrusive" methods, the reader is referred to [START_REF] Xiu | Fast numerical methods for stochastic computations: a review[END_REF]. PC models (both obtained with SG and NI) have been proven to be very efficient in many cases when the QOI is smooth with respect to the input parameters but its accuracy highly deteriorates in presence of irregularities and discontinuities. In this article, the term of irregularities refers to a continuous behaviour needing a high truncation degree of the chaos expansion to be accurately approximated. As shown below, the number of terms of the expansion grows in factorial with the truncation degree and with the dimension of the parameter space (the number of input parameters), which can become a burden both in term of computational time and memory to approximate such quantities. Moreover, a discontinuous QOI suffers from Gibbs phenomenon, which does not vanish when the truncation degree is increased. Among other methods, Multi-Element (ME) methods were developed to address these issues [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): Error analysis and applications[END_REF][START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF][START_REF] Sargsyan | Uncertainty quantification given discontinuous model response and a limited number of model runs[END_REF]. The first NI ME method, ME-PCM [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): Error analysis and applications[END_REF] (the non-intrusive equivalent of the multi-element intrusive method ME-PC [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF]), consists in refining the parameter space into smaller and smaller hypercubes, and applying stochastic collocation methods on these elements to find the chaos coefficients of the local PC models. In that way, a QOI can be approximated locally by models with lower degrees, which can reduce the cost of the algorithm in time and in memory. In the case of a discontinuous QOI, the Gibbs phenomenon is confined in the hypercubes which contain the discontinuity so that its impact on the variance of the QOI is reduced. However, if the discontinuity or the irregularity is not aligned with the Cartesian axis of the parameter space, the number of hypercubes can grow dramatically. In [START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF], the minimal multi-element method has been developed to study a discontinuous QOI. Relying on polynomial annihilation edge detection [START_REF] Archibald | Polynomial fitting for edge detection in irregularly sampled signals and images[END_REF], this algorithm aims at splitting the parameter space along the discontinuities, which totally removes the Gibbs phenomenon as the QOI is smooth on each element. The same observations are used to divide the parameter space and to train the surrogate models, using the least orthogonal interpolation [START_REF] Narayan | Stochastic collocation methods on unstructured grids in high dimensions via interpolation[END_REF] on every element. Similar methods split the parameter space along the discontinuity using Bayesian inference to detect the discontinuity and to approximate its curve equation [START_REF] Sargsyan | Uncertainty quantification given discontinuous model response and a limited number of model runs[END_REF]. Nevertheless, these two methods only focus on discontinuous behaviours but do not simplify the study of irregular continuous behaviours as ME-PCM did. Moreover, to the author knowledge, these methods do not permit to split the parameter space in presence of a discontinuity which runs partially through it. In this article, a new NI ME method is developed. This method aims at partitioning the parameter space according to a potential discontinuity or efficiently partitioning it to use local models of lower degree in presence of irregularities. This method is called "ME-ACD", which stands for Multi Element method based on Agglomerative Clustering on Derivatives. To split the parameter space, the used method is agglomerative clustering on the training observations represented in an well chosen abstract space taking into account the value of the QOI and the values of its derivatives with respect to the input parameters. This method also permits to split the parameter space in presence of a discontinuity which runs partially through it. The same training observations are used for both partitioning the parameter space and for training the local surrogate models with linear regression. Once the piecewise model is obtained, the observations are labelled with an index stating in what element they are in and a supervised classifier is trained to determine which local polynomial model is used for a given evaluation. The chosen classifier is the Multi-Layer Perceptron (MLP) classifier, based on artificial neural network. The use of clustering to study discontinuous UQ has already been done in [START_REF] Rai | Sparse low rank approximation of multivariate functions-Applications in uncertainty quantification[END_REF] and [START_REF] Sargsyan | Dimensionality reduction for complex models via bayesian compressive sensing[END_REF] but both are only driven by the QOI values on the observations. To the author knowledge, the development of a method which selects the most suitable partition of the parameter space using clustering on a abstract space containing the derivatives of the QOI had never been done in the field of UQ.

The article is outlined as following : PC framework and its extension to multi-element is developed in Section 2. Section 3 presents in details the different steps of ME-ACD. Section 4 shows its performances on a set of applications. The accuracy of ME-ACD is compared with the one of an algorithm of Gradient Boosted Trees (GBT) [START_REF] Jerome H Friedman | Greedy function approximation: a gradient boosting machine[END_REF] and with the accuracy of classical NI PC method on a set of discontinuous and irregular continuous QOIs. GBT is not a PC method, but a ML one which can be used for regression. As NI PC methods are classical regression or interpolation methods used in the context of UQ with a polynomial surrogate model, GBT can be used in this same context, except that it does not give orthogonal polynomial series. GBT are adapted to the study as they are said to be able to approximate even irregular and discontinuous QOIs. In a similar way, the method developed in this article in the context of UQ could be used in other approaches as a new regression method.

from PC expansion to ME model

In this section, notations and the standard PC framework are introduced and its limits are underlined. Then, some NI methods are recalled. Finally, the ME approach is introduced.

QOI and input parameters

The studied QOI is denoted f . f can be a given function, the solution of an Ordinary Differential Equation system at a given instant or the solution of a Partial Differential Equation system at a given location and instant. In this work, the possible dependencies in time and in space of the QOI are ignored to only focus on its dependency in the input parameters, without loss of generality. f depends on d stochastic input parameters p = (p 1 , ..., p d ) T which are supposed stochastic and independent. With this hypothesis, p follows the following multivariate pdf

ρ(p) = d m=1 ρ m (p m ) (1)
where ρ m is the pdf of p m . All the pdfs are supposed known. f is therefore stochastic and is supposed having a finite variance.

PC expansion

PC methods aims at approximating the dependence of the QOI f with respect to the stochastic input parameters p as a sum of multivariate orthogonal polynomials. As orthogonal polynomials take input parameters in a standard form, a bijective mapping is needed. In [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF], every input parameter is mapped with a standard Gaussian variable and multivariate Hermite polynomials are selected as orthogonal basis for the polynomial expansion. Isoprobabilisitic transformations can also be used to map every input parameter to an input parameter following the same kind of stochastic distribution but which is standardized. In that way, the used multivariate orthogonal polynomials are the ones being products of univariate polynomials corresponding to the standard distribution law of the input parameter through the Wiener-Askey scheme [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]. In this article, another mapping, developed below, is considered. For NI methods, as for ML methods, a sampling of the input parameters is effectuated. In the case of UQ, the observations are representing stochastic parameters, so the sampling is usually done in respect to their distribution. Observations are then mapped. The orthogonal properties are then used, notably to estimate the moments of the QOI and the Sobol sensitivity indices [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF].

With the ME decomposition, the orthogonal properties of the polynomial family cannot be used as local models only fits data on its element and is ignored in the rest of the parameter space, whereas the orthogonal properties need information of the model on the whole parameter space. A method permitting to use the orthogonal properties of the polynomials and then estimating Sobol indices with a ME model is the subject of further researches. Here, as orthogonal properties are not required for the learning of the model, the canonical multivariate polynomials basis

(Ψ k ) is used. For k ∈ N d Ψ k = d m=1 X km (2)
Moreover, a bijective mapping from p to x, where every x m follows an uniform law on [0, 1], is used. To do this, another hypothesis is considered on top of the others. Every ρ m is strictly positive (instead of being only positive), so that its cumulative distribution function C m is bijective.

x m = C m (p m ) (3)
From now on, the input parameters refer to x, the parameter space refers to X = [0, 1] d , the unit ddimensional hypercube which contains all possibles x and the QOI refers to f

f (x) = f (C 1 (p 1 ), ..., C d (p d )) = f (p) (4)
Once the model is trained, the same mapping is necessary for evaluating it with new values.

The truncation degree of the chaos series is denoted n ∈ N. There are several way to truncate a series of multidimensional polynomials, which leads to a finite number of terms M . The classical truncation scheme consists on selecting all the polynomials with

||k|| 1 = d m=1 k m ≤ n. This set of polynomials has a size of M (n, d) = (n + d)! n!d! (5)
M highlights the famous "curse of dimensionality", which underlines the dramatical increase of the number of terms of the series with d and n. For high values of n and d, the computational cost in time and in memory can become problematic. From this moment, the multi-index k is not used anymore and is replaced by a scalar index j. The polynomials are indexed from 1 to M (n, d) by increasing truncation degree. The order between polynomials with the same truncation degree is not important. Moreover, the dependence on n and d is sometimes omitted to simply denote the number of chaos terms M . In the approach, the standard computational PC model f is then written

f (x) = M j=1 c j Ψ j (x) (6)
PC model f converges fast when the QOI f is smooth. However, as stated above, the model is not accurate when the QOI is irregular or discontinuous. Figure 1 underlines this loss of accuracy for a 1D piecewise constant QOI f (1)

f (1) 1 (x) = -1 if x < 1 2 1 otherwise (7)
where x follows an uniform law on [0, 1].

The Gibbs phenomenon induces that even by increasing the truncation degree, the oscillations close to the discontinuity are not vanishing. In this example, a piecewise surrogate model would have been very accurate, approximating only the two constant parts of the function independently with 0 th -degree local polynomial models.

(dashed lines) and of its piecewise polynomials model (called Multi-element model below). Pseudo-spectral projection [START_REF] Xiu | Efficient collocational approach for parametric uncertainty analysis[END_REF][START_REF] Lin | Stochastic computational fluid mechanics[END_REF], one of the stochastic collocation method, uses the orthogonality of the chosen chaos polynomials basis to find the expression of the coefficients with an integral. Quadrature rules are then used to evaluate the integral. This implies that the training observations are structured. The structure can correspond to the Gauss-Legendre grid or the Smolyak sparse grids [START_REF] Abramovich | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF] for example. In the context of this article, the ME approach based on clustering implies that the observations are unstructured, and the use of Pseudo-spectral projection method is not possible.

Non-intrusive methods

Two NI methods are able to get the chaos coefficients with an unstructured observation grid : interpolation on unstructured grid and linear regression. Interpolation on unstructured grid is developed in [START_REF] Narayan | Stochastic collocation methods on unstructured grids in high dimensions via interpolation[END_REF]. The surrogate model is forced to be equal to the real value of the QOI at each training observation. This method is used in [START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF][START_REF] Yous | An adaptive minimum spanning tree multielement method for uncertainty quantification of smooth and discontinuous responses[END_REF].

In the present study, Ordinary Least Square (OLS) method, a linear regression method, is used to find the coefficients. Contrary to interpolation on unstructured grid, regression methods [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Guo | Constructing least-squares polynomial approximations[END_REF] do not imply that the surrogate model is equal to the real value of the QOI on the training observations. The aim of OLS is to minimize the euclidean distance between the model value and the real value on the N training observations

C opt = argmin C ∥ AC -B ∥ 2 (8)
where C ∈ R M is the vector containing the unknown chaos coefficients, A ∈ M N,M is the design matrix and B ∈ R N is the vector of predicted values. For the particular case of finding coefficients of a polynomials model :

A ij = Ψ (d) j (x i ), B i = f (x i ) and C j = c j
The problem is equivalent to solve this matrix equation, usually referred to as the normal equations

A T AC = A T B ( 9
)
The solution is therefore

C opt = A + B (10)
where

A + = (A T A) -1 A T .
The number of terms of the series M (determined in the case of finding the chaos coefficients by the chosen truncation degree n and the dimension d through equation ( 5)) is a crucial variable for regression. Compared to the number of training observations N , if M is too low, the regression is poorly accurate (usually called "underfitting" or "bias" in the ML field). If M is too high, what is called "overfitting" or "variance" in the ML field can occur. Overfitting is the fact that the model is accurate on the training observation but do not generalize well on other observations, thus leading to a poorly accurate model in between the training observations. These aspects can be taken into account by testing all regressions with the truncation degrees between the chosen n min and n max and selecting the optimal model as the one minimizing a cross validation error.

The considered error has to be an error between the surrogate model values and the real values of the QOI evaluated on the training observations. In that way, this error can be used when the algorithm is running to select the most suitable degree of a model and more generally the best hyperparameters to use in an range of values. Using the Mean Square Error (MSE) between the local surrogate model values and the values of the QOI on the training observations would give very poor information of the model accuracy if overfitting occurs. In that way, cross-validation is used. The Predicted Residual Sum of Squares (PRESS) error (based on the Leave-One-Out error) used in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] is really efficient in the context of PC models learned with linear regression. PRESS error is defined by

δ = 1 N N i=1 f (i) (x i ) -f (x i ) 2 (11)
where f (i) is the surrogate model learned on the set T \{x i }.

Fortunately, in the context of OLS, a method exists to express analytically every term of the sum without having to train a new model on the set T \{x i } for every observation x i . This formula [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF], which only needs the model trained on the observations of T , leads to

δ = 1 N N i=1 f (x i ) -f (x i ) 1 -h i 2 (12)
where h i is the i th diagonal term of the matrix H = AA + known as the "Hat matrix" or the "projection matrix". This error correctly gives information on the accuracy of the local surrogate models, even when overfitting occurs.

To compute A + and the PRESS error, the Singular Value Decomposition (SVD) of A is used [START_REF] Kalman | A singularly valuable decomposition: the svd of a matrix[END_REF]. If A is a N × M matrix, finding the SVD of A consists in finding a N × N orthogonal matrix U, a N × M diagonal rectangular matrix Σ and a M × M orthogonal matrix V so that

A = UΣV T (13)
The columns of U are called "left singular vectors" of A, they are the eigenvectors of AA T . The columns of V are called "right singular vectors" of A, they are the eigenvectors of A T A. The diagonal terms of Σ are the so-called "singular values" of A, which are the eigenvalues of AA T . SVD decomposition is not unique, but what is usually done is to take the one which leads to the singular values ordered by decreasing value in the diagonal of Σ. If the rank of A is equal to r, A has r non-zero singular values.

SVD is really useful to compute efficiently A + . It is proven in [START_REF] Kalman | A singularly valuable decomposition: the svd of a matrix[END_REF] that

A + = VΣ + U T (14)
where Σ + is the diagonal rectangular M × N matrix with the inverse of the non-zero singular values on the diagonal and zeros elsewhere.

SVD is also useful to compute the "Hat matrix" and then the PRESS error. Indeed, we have

H = AA + = UΣV T VΣ + U T = UΣΣ + U T = U r U T r (15)
where U r is the square matrix containing the first r rows and columns of U. The diagonal terms of H, needed for the computation of the PRESS error, are then equals to

h i,i = r j=1 (U r ) i,j 2 (16)
Not all the rows and the columns of U and V have to be computed to get A + and H. Taking only the r first columns of U, Σ and the first r rows of V T is known as "compact SVD".

Algorithm 1 underlines the pseudo-code to compute OLS and the PRESS error.

Note : in all the pseudo-codes of this article, X ← value means that variable X gets the given value. Moreover, a ":" used as an index of a vector or a matrix denotes all the rows or all the columns of the item.

If A is ill-conditioned, the user can ignore a set of singular value less than a given threshold, thus considering r ≤ r singular values instead of r.

Algorithm 1 Pseudo-code of the computation of OLS and PRESS error

A ← N × M design matrix B ← vector of size N containing the values of the QOI on the training observations U, Σ, V ← SVD(A) (compact form) Σ + ← min(N, M ) × min(N, M ) matrix with 1 σ i,i on its diagonal and zeros elsewhere

X ← VΣ + U T B H ← r j=1 (U ) i,j
2 (Here H is the vector containing the diagonal terms of the "Hat matrix" and not the matrix itself)

Y = AX δ ← 1 N N i=1 Y i -B i 1-H i 2 Return X, δ
Let's consider that N ≥ M without loss of generality (if M > N the operations can be conducted on the transpose of the matrix). In the most time expensive case, r = M . Time complexity of a matrix product between a N × M matrix A and a M × P matrix B is in O(N M P ). In that way, time complexity of calculating X is in

O(M 3 + M 2 N + M N ) = O(M 2 N ).
To the author knowledge, the complexity of compact SVD is O(M 2 N ). Time complexity of the algorithm of OLS and PRESS Error is then O(M 2 N ).

In [START_REF] Guo | Constructing least-squares polynomial approximations[END_REF], the authors state that a standard random sampling often lead to suboptimal approximations, and show the existence of a special sampling leading to nearly optimal solutions. However, as it is described in section 2.4, in the current paper, regression is used on subsets of the sampling, which induces that for every approximation, the user does not have the control on the position of the training observations which are used in the training of the local surrogate models.

Multi-element approach

As stated above in section 2, to approximate an irregular or a discontinuous QOI, one solution is to use a ME method. ME methods consist into partitioning the parameter space X into N E non-overlapping subsets E e called elements

X = N E e=1 E e (17)
E e 1 E e 2 = ∅ if e 1 ̸ = e 2 (18)
Once X is partitioned, a local surrogate model fe is obtained in every element E e with the training observations located in this element. The global surrogate model is then a piecewise model composed of every local model. The global piecewise model is written

f (x) = N E e=1 fe (x)I Ee (x) = N E e=1 Me j=1 c ej Ψ (d) j (x) I Ee (x) (19)
with M e the number of multivariate monomials of the local PC model fe . I Ee is the indicator function of the element E e satisfying

I Ee (x) = 1 if x ∈ E e 0 otherwise (20)
In numerical applications, the classifier plays the role of the indicator function. Once the global model is trained, to evaluate the model on a new point, the classifier detects in which E e the observation lies to know which local model to use. Every coefficient of the local polynomials models fe is obtained with the OLS algorithm on the training observations located in E e .

ME-ACD algorithm

Abstract spaces

In the following algorithm, different abstract spaces are considered:

• P : the original d-dimensional parameter space of points of coordinates (p 1 , ..., p d ).

• X : the mapped d-dimensional parameter space containing points of coordinates (x 1 , ..., x d ). This space is the one referred to as the parameter space.

• F : the d + 1-dimensional space containing points of coordinates (x 1 , ..., x d , f (x)).

• D : the d + 1-dimensional space containing points of coordinates ( f (x), ∂ f ∂x 1 (x), ..., ∂ f ∂x d (x)
). In that way, an observation can be considered in those different spaces. The objective is to split the parameter space X into a partition of N E elements E e . To do this, agglomerative clustering is used. The results of the clustering are different according to the space the observations areconsidered in. The space has to be chosen so that the clustering efficiently splits the parameter space both in the case of a discontinuous QOI and in the case of a continuous but irregular QOI.

Step 1 -Sampling of the parameter space

The first step of the algorithm consists in properly sampling the parameter space X to get a training set T = {x i } N i=1 of N training observations. One aim of ME-ACD is to work on every sampling, so that it could be used with experimental data. In that way, an adaptive sampling is not used, conversely to other methods [START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF][START_REF] Yous | An adaptive minimum spanning tree multielement method for uncertainty quantification of smooth and discontinuous responses[END_REF]. In this article, Sobol quasi-random sequences [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] are used. Sobol sequences, as Latin Hypercube Sampling, aims to sample a d-dimensional space leaving less sparse zones in the hypercube than purely random Monte-Carlo sampling. More information on the Sobol sequences can be found in [START_REF] Joe | Notes on generating sobol sequences[END_REF]. The initial direction numbers used in the article are the one suggested on the following website : https://web.maths.unsw.edu.au/ fkuo/sobol/ (last accessed : April 2022).

With this strategy, the unit hypercube is sampled and the corresponding p i of input parameters of training observation x i can be obtained by inverting formula [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF]. In some cases, a slight modification on the sampling can be performed for the observations located on the hyperfaces of the unit hypercube. They can be translated of a distance ϵ = 10 -3 perpendicularly to the hyperface (ϵ is added to every components of every x i which are equals to 0 and is subtracted to every components equals to 1). This avoids having an infinitely far point for normal distribution for example. Figure 2 shows a Monte-Carlo sampling and a Sobol sequence for a 2D example. 500 observations of p with p 1 following an uniform law on [1, 1.5] and p 2 following a normal law of mean 10 and of standard deviation 5 sampled with the two strategies are represented in X and in P. Monte Carlo and Sobol samplings are very close to each other, but the latter is elaborated to guarantee a lower discrepancy.

Once the sampling is performed, the QOI is evaluated on every training observations of T to get a vector B ∈ R N with :

B i = f (x i ) = f (p i ) (21)
Time complexity of the computation of a Sobol sequence is in O(N ×d). Time complexity of evaluating the QOI is independent of the algorithm.

Step 2 -Computation of the global design matrix

In Step 4, the parameter space is partitioned and OLS with a multidimensional polynomial model of degree varying from n min to n max is used on every element to find the chaos coefficients of the local models. A major advantage of using OLS in this approach is that a single design matrix has to be computed once. This matrix will be referred to as the global design matrix. When OLS is called to learn the coefficients of a local model, only a submatrix of this global design matrix is used. The rows of the submatrix are the rows corresponding to training observations located in the element. Its columns are the first M (n, d) columns corresponding to the chosen truncature degree n. This global matrix, denoted A, is a N × M (n max , d) matrix, M (n max , d) being the number of terms of the polynomial series of degree n max (see formula ( 5)). The global design matrix terms are written

A i,j = Ψ j (x i ) (22)
As a reminder, j is a scalar index used to denote the multidimensional canonical polynomial Ψ j . For theses polynomials :

Ψ j = d m=1 X km m (23) n = d m=1 k m (24)
with n ∈ [0, n max ] being the degree of Ψ j .

Computing this global matrix in a naive way would result in recomputing every components (x i ) k l l for the calculus of every Ψ j (x i ) thus leading to a complexity in O(d 2 ) for the computation of every Ψ j (x i ).

Computing the global matrix would then have a complexity in

O N × M (n max , d) × d 2 .
Instead of using this naive computation, for every training observations x i , we can store the values of every

(x i ) km m in a d × n max matrix D (i) with D (i) m,l = (x i ) l m .
Every term A i,j of the global design matrix is then the product of terms of D (i) :

A i,j = d m=1 D (i) m,km (25)
With this computation scheme, the computation of a matrix

D (i) is in O(d × n max ). The computation of the term Ψ j (x i ) with these matrices is in O(d). Time complexity of computing the global matrix is then in O N × d × n max + N × M (n max , d) × d = O N × M (n max , d) × d .
M ← (nmax+d)! nmax!d! Initiate a N × M matrix A with zeros for i ∈ [1, N ] do x ← T i Initiate a d × n max matrix D with zeros for l ∈ [1, d] do for m ∈ [1, n max ] do D l,m ← x m l end for end for for j ∈ [1, M ] do k ← multidimensional degree corresponding to the mapped scalar degree j a = 1 for l ∈ [1, d] do a ← a × D l,k l end for A i,j ← a end for end for

Step 3 -Estimation of the derivatives

In Step 4, agglomerative clustering is used on the observations considered as points of D. Therefore, the values of the derivatives of the QOI on the unstructured training points grid need to be estimated. To estimate these values, the first order Taylor's theory is employed. For two different observations x and x

′ f (x) ≈ f (x ′ ) + ∂ f ∂x 1 (x ′ ) × (x 1 -x ′ 1 ) + ... + ∂ f ∂x d (x ′ ) × (x d -x ′ d ) (26)
In that way, the values of the derivatives of the QOI on a given point x ′ can be obtained by using OLS algorithm, solving the previous formula with its nearest neighbours.

{x i } N neigh i=1
denoting the set of the N neigh nearest neighbours of x ′ , the following system has to be solved

   x 1 1 -x ′ 1 • • • x 1 d -x ′ d . . . . . . x N neigh 1 -x ′ 1 • • • x N neigh d -x ′ d        ∂ f ∂x 1 (x ′ ) . . . ∂ f ∂x d (x ′ )     =    f (x 1 ) -f (x ′ ) . . . f (x N neigh ) -f (x ′ )    (27)
The observations considered as neighbours of x ′ have to be determined. A first idea would be to consider as neighbours the observations corresponding to the nearest points in X . Nevertheless, if the QOI is discontinuous, the calculus of the gradient between two observations which are from different sides of the discontinuity is going to be highly altered as the QOI is not derivable in this zone. In that way, the N neigh closest neighbours are going to be determined as the nearest neighbours of x ′ in the space F, so that in this calculus, only points on the same side of the discontinuity are going to be considered as neighbours. The Ball tree algorithm [START_REF] Stephen | Five balltree construction algorithms[END_REF] is used to search the neighbours. To determine the best number of neighbours, OLS algorithm is effectuated for N neigh varying from k min to k max and the one minimizing the PRESS error is selected. Now the observations can be represented in the space D and the agglomerative clustering to test the different partitions of the space can be used with these points. Algorithm 3 is the pseudo-code for Step 3. This step contains the construction of a Ball tree on the training set considered as N points of a d + 1-dimensional space. Construction time of the ball tree is in O(N log(N )) [START_REF] Stephen | Five balltree construction algorithms[END_REF]. According to author knowledge, query time in this ball tree is in O(d log(N )), so the time of all the queries for nearest neighbours is in O(dN log(N )). This step also contains the construction of a design matrix for every observations for a total time complexity in O(N × dk max ), and k maxk min + 1 regressions with submatrices of size N neigh × d of that matrix. If k min and k max are proportional to the dimension d, time complexity of the k maxk min + 1 regressions for a an observation is in O(d 4 ). Time complexity of this step is more likely to be dominated by the regressions and is in O(N × d 4 ). Agglomerative clustering follows a "bottom-up" approach : each point of the considered space begins with its own cluster and at each step, the 2 closest clusters are merged together, until the wanted number

Step 4 -Partition of the space and training of the local models

for i ∈ [1, N ] do x ← T i I ← J i,:
Initiate a k max × d matrix A and a vector F of length k max with zeros

for m ∈ [1, k max ] do F m ← B i -B Im for l ∈ [1, d] do A m,l ← x l -T Im,l end for end for Initiate a d × (k max -k min + 1) matrix C and a vector E of length k max -k min + 1 with zeros for m ∈ [k min , k max ] do C :,m-k min +1 , δ m = OLS (A 1:m,: , F 1:m ) end for k opt ← argmin m∈[k min ,kmax] δ m D i,: ← C :,kopt end for return(D)
of clusters is reached. Ward clustering is used. In this clustering, the two clusters having the minimal sum of square of the distance between pairs of points (one from one cluster and the other from the other cluster) are merged at each step. To take advantage of the bottom-up approach, only the partition with N E min elements is made. The partitions from N E min + 1 to N Emax corresponds to the last steps of the agglomerative clustering with N E min elements, and are also accessible. From the partition with N E = N Emax elements, two of the elements are merged to get the partition with N Emax -1 elements and so on until the partition with N E = N E min elements is reached. In that way, elements of partitions with N E < N Emax elements are unions of the elements of the partition with N Emax elements. A compact way to store this information is to give a label to every training observations from 1 to N Emax corresponding to the index of the element E Np,e the corresponding points belongs (partition with N Emax elements). As elements E p,e of a given partition are union of elements E Np,e , an element can be stored as a list of the indices corresponding to elements E Np,e which compose it.

The information of all the partitions obtained by clustering can be stored in :

• a vector L of size N containing the labels for every training observations, labels corresponding to the index e of the element E Np,e (maximal partition) they belong.

• N p lists of elements E p , one by tested partition, containing N E elements E p,e . As E p,e is the union of some elements E Np,e ′ , one element is represented as a list of labels, labels corresponding to the index e ′ of the elements E Np,e ′ which compose it.

Figure 3 illustrates this idea in a tree like structure diagram. In this example, N E min = 2 and N Emax = 4. Elements are represented as circles and every step of agglomerative clustering is a depth of the tree. In the example, L is a vector of size N containing labels from 1 to 4 and they are 3 lists :

• E 1 is a list containing the elements E 1,e being stored as

E 1,1 = [1] and E 1,2 = [2, 3, 4]
• E 2 is a list containing the elements E 2,e being stored as

E 2,1 = [1], E 2,2 = [2, 3] and E 2,3 = [4]
• E 3 is a list containing the elements E 3,e being stored as With this data structure, for a given partition, the observations located in an element E p,e are the ones with a label L equals to one of the indices in the list which store the information of E p,e . For example, an observation located in E 1,2 in the above example has a label of 2, 3 or 4 (corresponding to the elements of the partition with N E = N Emax elements).

E 3,1 = [1], E 3,2 = [2], E 3,3 = [3] and E 2,4 = [4] N E = N Emax N E = N Emin E 3,1 E 3,2 E 3,3 E 3,4 E 2,1 E 2,2 E 2,3 E 1,1 E 1,2 Ω
Applying the clustering step on observations considered in different spaces give different results. Making a clustering in X would not have any interest. Making a clustering of the observations considered as points in F would efficiently split the parameter space for a discontinuous QOI, but would not be efficient in the case of a continuous irregular QOI. An irregular behaviour of a QOI is due to a high variation of its derivatives. The idea is then to make clusters in which the QOI has "close" derivatives' values. In that way, agglomerative clustering on the observations considered in D is used. Therefore, the space is split according to discontinuities, as the value of the QOI is present in the coordinates of a point in D, and according to irregular behaviours, as the values of the derivatives are also present in the coordinates. Figure 4 represents the plots of observations of two 1D-QOIs f (1) 1 (see equation ( 7)) and f (1) 2 (see equation ( 28) below) considered in the spaces D and F.

f

2 is an irregular QOI defined by

f (1) 2 (x) = 10 tanh 10(x - 1 2 ) (28)
As d = 1 for the examples of the figure, the considered sampling is an uniform 1D grid with a constant step. Figure 4 One can add a connectivity diagram between the training observations represented in X before the clustering (as an optional step between step 3 and step 4). With a connectivity diagram, only two clusters which have at least one observation connected together can be merged. This connectivity reduces the chances to gather two clusters which are far from one another in X . If the user decide to use connectivity, one way to construct the diagram is to connect an observation with its k connect nearest neighbours in X , k connect being a hyperparameter to choose. The choice of using or not a connectivity diagram is modelled by the boolean isConnect.

A low accuracy of the gradient calculus of step 3 can alter the partition but if a point is miscalculated with this method, it can be seen as an outlier in D and is more likely to be isolated by the partition approach, however creating unnecessary elements. Some action could then be taken to remove the observation from the training set without counting a supplementary element. In the scope of this article, no outlier removing actions are taken.

Once the N p partitions have been found, the one leading to the most accurate model has to be selected. For every partitions, local models are trained on the elements. The selection criterion of the optimal partition is based on the PRESS errors of the local models. The partition of N E ∈ [N E min , N Emax ] elements is considered. For each of its elements E p,e , a local surrogate model is trained using the OLS on the observations located in the element to find the coefficients of the chaos series. As stated in 2.3, the optimal truncation degree of the polynomials model is obtained minimising PRESS errors. For every local surrogate model, regressions are carried out with a truncation degree from n min to n max . The selected truncation degree is the one minimizing the PRESS error δ p,e,n of the considered local surrogate model fe trained with the observations of E p,e of the considered partition p with N E elements with a truncation degree n :

δ p,e = min n∈[n min ,nmax] δ p,e,n (29)
This error δ p,e represents the accuracy of the local surrogate model fe trained on the observations located in E p,e and having an optimal degree between n min and n max . The error δ p , image of the accuracy of the considered partition p with N E elements, is defined as the maximal value of the errors δ p,e of all the local models trained on its elements with an optimal degree :

δ p = max e∈[1,N E ] δ p,e (30)
To get the optimal partition, we first consider the one with N E min elements. Local models for every of its elements are trained with optimal degree. From one partition to the next one, only one element is split into two, going down the agglomerative clustering tree. In that way, only the local models on the two new elements can be trained at each step, the other ones having already be trained for previous partitions. During the different partitions of the space, a test can be added : if the selection error δ p of the current tested partition is lower than a given threshold δ min , a suitable partition is already found and there is no need to look for others with a lowest δ p error. This test is referred as the early stop.

As stated above in section 2.3, one major advantage of using OLS in this approach is that the global design matrix A, of size (N × M (n max , d)), is only computed once, at Step 2. When a regression on a element is wanted, only a sub-matrix composed of the lines corresponding to the observations located in the element and of the columns corresponding to the considered truncation degree (the first M (n, d) columns for a given degree n) is used. Algorithm 4 is the pseudo-code for Step 4. Time complexity of this step is given by the construction of the agglomerative clustering with N E min elements and by all the regressions. Ward clustering time complexity is in O(d × N 3 ). In a case without early stop, N E min + 2 × N Emax -N E min local models are trained. Every training needs n maxn min + 1 regressions to determine the local models with optimal degree, dominated by the one with n = n max . Time complexity of all these regressions is then in

O N E min + 2 × N Emax -N E min × max N, M (n max , d) × min N, M (n max , d) 2 .
Algorithm

Step 5 -Training of a classifier

Once the most suitable partition is obtained, with its set of N E elements and N E local surrogate models, a classifier is needed. Indeed, for further model evaluations, the first step is to determine in which element E e the evaluation is, to know which local model to use to perform the evaluation. The use of a local model corresponding to a wrong element leads to huge errors, polynomial models being potentially highly divergent far from their training observations. To model the action of the indicator function denoted I Ee in formula [START_REF] Rai | Sparse low rank approximation of multivariate functions-Applications in uncertainty quantification[END_REF], a classifier is used.

The user is free to choose any supervised classification method. In this article, Multi Layer Perceptron (MLP) classifier is chosen. MLP classifier is an artificial neural network used in a purpose of classifying observations labelled with class numbers. In Step 4, a label corresponding to the element in the partition with N E = N Emax is given to every training observation. In Step 5, training observations are relabelled with this time a class number equal to the index of the element in the partition with N E elements, N E being the optimal number of elements found at the end of Step 4.

In all the below applications, 2 hidden layers of 50 neurons are taken. For the hidden layers, Rectified Linear Unit activation function is used. For the output layer, Softmax function is used. Cross-entropy cost function is minimised. The training of the classifier is done using backpropagation with Adam method [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] (with the classical decay values of β 1 = 0.9 and β 2 = 0.999, a mini-bacth size of 200 observations and a max number of iterations over the entire training set of n epochs = 200 with an early stop if the weights have converged). Several classifiers are trained with different learning rates. The training which maximises the number of well classified training observations, giving the best accuracy being the fraction of well classified observations. In the article, a naive method testing all N l learning rates is used. No regularisation term is added as classification overfitting is not a problem in the approach (all the training observations have to be well-classified to avoid using a wrong local model).

To the author knowledge, time complexity of one optimisation step with feedforward, backpropagation and Adam bias momentum correction is in O (d × n neurons + n 2 neurons + n neurons × N E ) × N , with n neurons being the number of neurons in the hidden layers (50 in our case). In that case, the learning of one classifier is in in O(n epochs × n 2 neurons × N train ), where n epochs = 200 in the worst case is the number of steps of optimisation going on all the training set. So the learning of all classifiers is in O(N l ×n epochs ×n 2 neurons ×N ).

The input of the classifier is a given point in X (a set of value of input parameters) which need to be evaluated by the piecewise model. This point is denoted the evaluation point. The output is an array of N E components containing the probabilities that this evaluation point is located in the element E e . Therefore, a test can be added to check if the evaluation point lies in an indecision zone : if the highest probability is lower than θ, θ ∈ [0, 1] being a threshold to choose, the evaluation point is considered as located in an indecision zone. In this case, the point is not evaluated by the model and is ignored. θ is a trade-off parameter. The greater θ is, the more prudence is taken, but the more evaluation points close to the boundaries of the elements are ignored. With θ = 0, no evaluation point is ignored but the ones located in the indecision zone are going to alter the accuracy of the model. By increasing θ, an evaluation point is more likely to be ignored if it is potentially affecting the accuracy. After a given value, rising θ induces that some evaluation points are ignored even if they are not really altering the accuracy. In [START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF], the minimal multi-element method presents a similar procedure, using the resolution of the sampling strategy to ignore points lying within a distance of the boundary lower than this resolution. θ is a parameter which can be tuned after the training of the model.

Summary of ME-ACD algorithm

The numerical parameters and hyperparameters of the algorithm are :

• N : the number of observations generated in the Sobol sequence. It determines the size of the training set.

• N E min : the minimal number of wanted elements.

• N Emax : the maximal number of wanted elements.

• N neigh-min : the minimal number of neighbours in F used in the estimation of the derivatives.

• k max : the maximal number of neighbours in F used in the estimation of the derivatives.

• n min : the minimum chaos degree a local surrogate model can have

• n max : the maximum chaos degree a local surrogate model can have.

• δ min : an optional threshold to early stop the search of partition if one suitable partition has already been found.

• A set of N l learning rates to get the optimal classifier.

• isConnect : a boolean indicating if a connectivity diagram is used.

• (Optional) k connect : number of connected neighbours an observation have in the optional connectivity diagram

For the classifier, hyperparameters are given in section 3.6. The caution threshold θ ∈ [0, 1] of the classifier is a post-processing parameter. Once the model is obtained with its attached classifier, the user is still free to tune θ without launching the training once again.

For agglomerative clustering, research of the nearest neighbours and MLP classifiers, the Scikit-learn toolkit has been used [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].

In term of time complexity :

• Step 1 : O(N × d) + N × time complexity for computing a value of f • Step 2 : O N × M (n max , d) × d • Step 3 : O(N × d 4 ) (if k min
and k max are set to be proportional to d)

• Step 4 : O N E min + 2 × N Emax -N E min × max N, M (n max , d) × min N, M (n max , d) 2 • Step 5 : O(N l × n epochs × n 2 neurons × N )
The dominant step is dependent of the choice of the hyperparameters and of the dimension of the QOI.

Applications

In this section, performances of ME-ACD are demonstrated on a set of test functions. To evaluate these performances, the possible causes of errors between the piecewise surrogate model and the QOI it approximates have to be identified. These errors can come from 3 phenomenons :

• Underfitting (or bias) : the optimal degrees of the local polynomials models are too low to approximate the QOI on the elements. Increasing N , n max and/or N Emax could solve the issue.

• Overfitting (or variance) : as the partitions of the parameter space and the optimal degrees of the local models are selected with the PRESS error, overfitting is not likely to occur. If it still occurs, lowering n min and/or increasing N Emax could solve the issue.

• Misclassification : as the polynomial models are highly divergent far from their training observations, if the evaluation point is misclassified in a wrong element, the wrong local model is used and the accuracy is highly altered. As a reminder, the caution threshold θ is here to avoid these errors, ignoring the point when the probability of being in the element given by the classifier is strictly lower than it. Increasing N or adding training observations close to the boundaries of the elements and/or modifying the classifier parameters could solve the issue. Increasing θ makes the algorithm stricter and makes it ignoring more evaluation points which are more likely to be responsible of misclassification errors.

For the classifier, the learning rate is the one within [10 -m , 2 × 10 -m , 5 × 10 -m ], ∀m ∈ [START_REF] Archibald | Polynomial fitting for edge detection in irregularly sampled signals and images[END_REF][START_REF] Genz | Testing Multidimensional Integration Routines[END_REF] given the classifier maximising the accuracy on the training set itself (the optimal classifier between N l = 21 is therefore selected). The threshold is set to θ = 0.85, which means that the evaluation point is not ignored by the caution procedure if the classifier is sure at more than 85% that it is located in a given element.

The performances of the algorithm are compared to the ones of the standard PC model of optimal degrees between a given n min-P C and n max-P C and of Gradient Boosted Trees (GBT) algorithm trained on the same test observations. In fact, in the applications, the standard PC model is obtained with the ME-ACD method with N E min = N Emax = 1. To set the 3 hyperparameters of GBT : the learning rate, the number of trees and the depth of the trees ; a validation set of N val = 10 5 observations for every training set is sampled. A brute force search is realised for every training set : every combination of hyperparameters is tried to train a GBT model and the one minimizing the MSE between the true value of the QOI and the value of the given GBT model on the validation observations is chosen. In the applications, the learning rates is located in [10 -m , 2 × 10 -m , 5 × 10 -m ], ∀m ∈ [START_REF] Archibald | Polynomial fitting for edge detection in irregularly sampled signals and images[END_REF][START_REF] Genz | Testing Multidimensional Integration Routines[END_REF], the number of trees in [100, 200, 300, 400, 500] and the maximum depth of the trees in [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF][START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): Error analysis and applications[END_REF].

To evaluate the precision of the algorithms, a test set of N test = 10 5 observations is randomly sampled. This test set is different from the validation ones used to tune the GBT models, but remains the same for a given convergence study. Conversely, the validation sets used to tune the hyperparameters of GBT models is different for every different training set (for every different N of the convergence study).

The precision of each model is given by E, the MSE between the true value of the QOI and the value given by the predictions of the model on the test set. To see the impact of the caution parameter, E θ is also studied, which is the same error but evaluated with the N θ test observations kept by the caution procedure of the ME-ACD model with θ = 85%. The value V θ=0.85 = N θ=0. 85 Ntest , fraction of the kept test observations by the caution procedure, is given for each N of the convergence study.

Discontinuous functions

In this section, a set of 3 different discontinuous functions taking d = 2 input parameters are studied :

f hyperplane = 1 if x 1 + 0.3x 2 -0.5 > 0 -1 otherwise (31)
This function presents a discontinuity which is a line of equation x 1 + 0.3x 2 -0.5 = 0.

f hypersphere = 1 if (x 1 -0.5) 2 + (x 2 -0.5) 2 < (0.25) 2 -1 otherwise (32)
This function presents a discontinuity which is a circle of center (0.5, 0.5) and of radius 0.25.

f multiple =        3 if (x 1 ) 2 + (x 2 ) 2 < (0.25) 2 2 else if x 1 < 0.5 -0.3x 2 1 else if (x 1 -1) 2 + (x 1 -1) 2 > 0.3 2 0 otherwise (33)
This function presents 3 discontinuities, two fractions of circles in the corners of the unit hypercube and the same line as in f hyperplane .

For these test functions, optimal elements are clearly defined. In that way, a true label can be provided indicating in which element a test observation is located. With these labels, the accuracy of the classifier, denoted F and representing the fraction of well-classified test observations, can be computed. For the next sections, as the Multi-element procedure is also used to study continuous irregular behaviour with polynomials series of fewer degrees, the optimal partition is not clearly defined and the computation of this accuracy is not possible.

Figure 5 shows 1000 training observations labelled by the clustering step of ME-ACD and represented in the parameter space for the 3 test functions. This first visualisation highlights that in that dimension, ME-ACD method is able to detect discontinuities and to partition the parameter space according to these discontinuities. functions on the test set. Table 1 contains the hyperparameters and the parameters selected for PC and for ME-ACD for these convergence studies. Hyperparameters for GBT are the best from the brute-force search described above in section 4 for every N of the convergence studies.

The convergence plots of E (Figures 6a, Figure 6c and Figure 6e) highlight that ME-ACD method is able to give a more accurate model than the 2 other methods, but error is still relatively high. As the functions are piecewise constant and the model is also piecewise constant(as n min = n max = 0), the cause of the errors for the ME model are missclassification errors. As F is close to 99% for every N , around 1% of the observations are responsible for all the error. They are the observations located close to the discontinuities.

Figure 7 shows this time a convergence plot of the MSE on the N θ test observations kept by the caution procedure for every N . With θ = 85%, V θ is around 90% for every N (shown by Figure 6b, Figure 6d and Figure 6f). In most of the case, all the problematic observations are removed by the caution procedure, leading to errors of magnitude 10 -30 , stating that no remaining observation is missclassified, or leading to errors of magnitude 10 -5 , highlighting than only around 10 of the N θ ≈ 10 4 kept test observations are missclassified. With θ = 85%, around 9% of the test observations are ignored despite the fact they are not leading to any error. The caution procedure is then able to detect the problematic observations and the user is free to act consequently to avoid these missclassification errors.

Even better classifier results (F and V θ ) could have been obtained with another sampling. Here, Sobol sampling is used to model a general case, leading to a homogeneous sampling no matter the potential discontinuities. In [START_REF] Yous | An adaptive minimum spanning tree multielement method for uncertainty quantification of smooth and discontinuous responses[END_REF], an adaptive sampling is developed. This sampling is able to get more observations close to discontinuities. If the user is free to sample the parameter space, more observations close to the discontinuities should lead to an even better accuracy of the classifier and to an even more efficient caution procedure.

Figure 8 shows the optimal degrees for the PC models in the convergence studies for the 3 test functions. The corresponding number of terms of the series is also shown on the left in parenthesis. With increasing number of training observations, PC is able to select higher optimal degrees without suffering from overfitting. Compared to these PC models, local polynomials models of the ME method have only one term and lead to a better accuracy. The economy in number of terms of the PC series induced by using lower degrees local model would be even more marked at higher dimensions.

Discontinuity which runs partially through the parameter space

In this subsection, the following QOI is studied :

f partial = 0.5x 2 -0.25 if x 1 >= 0.5 -0.5x 2 + 0.25 otherwise (34)
This piecewise function presents a discontinuity which "runs partially" through the parameter space. The function is linear on each side of the hyperplane x 1 = 0.5 and so has the geometry of an hyperplane. However, when x 2 becomes close to 0.5, the gap between the two surfaces is so thin that numerically, it becomes more and more regular (although the discontinuity is still here mathematically). To the author knowledge, none of the previous multi-element methods permits to split the parameter space in presence of a discontinuity which runs partially through it. As stated in [START_REF] Yous | An adaptive minimum spanning tree multielement method for uncertainty quantification of smooth and discontinuous responses[END_REF], this is a problematic for some of the modern methods in the multi-element domain.

Figure 9, containing results obtained with the parameters and hyperparameters gathered in Table 2, shows similar results for the convergence study of this function as for the discontinuous functions presented at the previous subsection. ME-ACD is then able to handle discontinuities running partially through the parameter space in the same way as any other discontinuity. Conversely to polynomial annihilation used in [START_REF] Archibald | Polynomial fitting for edge detection in irregularly sampled signals and images[END_REF][START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF], agglomerative clustering of the observations in space D is not affected by thinner gap close to the boundaries of elements.

Bifurcation function

One major advantage of ME-ACD is that it is able to accurately model continuous regular functions, continuous irregular functions and discontinuous functions. The adaptation of ME-ACD to model solutions of Ordinary Differential Equation (ODE) is under development. This adaptation will then be suitable to approximate solutions of ODE with the apparition of bifurcations (thus leading to irregularities and discontinuities with respect to the input parameters).

In this subsection, the following QOI is studied :

f bif -t (x 1 , x 2 ) = (1 + x 2 ) tanh( 1 2 (3x 1 -1)(1 + x 2 )t) ( 35
)
this function is the solution at time t of the following ODE :

ẏ = -p 1 (y 2 -p 2 2 ) (36)
with y(t = 0) = 0, p 1 following an uniform distribution law on [-1 2 , 1] and p 2 following an uniform distribution law on [START_REF] Archibald | Polynomial fitting for edge detection in irregularly sampled signals and images[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF]. The transformation to get the corresponding x appears in the form of the mapped function. This system models the apparition of a bifurcation in a system. The two equilibria of the system (if we ignore the punctual case p 1 = 0) are -p 2 , stable if p 1 < 0 and p 2 , stable if p 1 > 0.

At the first instant, the solution is regular in function of x. When t rises, the solution is more and more irregular with respect to x. Asymptotically (t → ∞), the solution presents a discontinuity which is a line of equation p 1 = 0 (x 1 = 1 3 ), symbol of a bifurcation with respect to p 1 (x 1 ).

The 3 considered instants are : t = 1 (regular), t = 10 (irregular) and t = 1000 (discontinuous). Figure 10 shows the variation of the function with x 1 for x 2 = 1 (the value of x 2 leading to the widest gap).

Similarly to previous sections, Figure 11 shows the results of convergence studies of the MSE on all the test set and of the kept fraction of test observations by the caution procedure of ME-ACD. Figure 12 underlines the convergence of the MSE E θ=0.85 on the kept observations. The hyperparameters and the parameters used for the functions in these studies are given by Table 3.

At t = 1, the early stop of ME-ACD makes the algorithm not using the ME procedure, thus being equivalent to PC, but with a lower n max . PC is then outmatching ME-ACD is this case. With a similar n max , both methods would have been identical. As N Eopt = 1 for all N , 100% of the test set is always kept by the caution procedure. At t = 10, ME-ACD split the parameter space into 3 or more elements to use local models of lower degrees. ME-ACD is able to accurately give an approximation of a continuous irregular QOI using lower degrees local models. As shown in Figure 12b, the caution procedure permits to have an even more accurate model, but as the QOI is not discontinuous, missclassification errors are less altering the accuracy. is located in this zone, and a 2 elements model with two hyperplanes is given by ME-ACD, explaining why the accuracy is better with a low N . Generating more observations close to the quasi-discontinuity with an adaptive sampling or ignoring these points seeing them at outliers should solve the issue.

f bif -1000 (x 1 , x 2 = 1) (c)
ME-ACD is then able to give an accurate model in the three cases of a regular QOI, a continuous irregular QOI and a discontinuous QOI, which makes it a model able to approximate a lot of QOIs and a good candidate to approximate solutions of ODE with potential bifurcations.

Genz function in dimension 5

In this section, a 5-dimensional function a QOI from the Genz family [START_REF] Genz | Testing Multidimensional Integration Routines[END_REF], depending on d = 5 input parameters is studied : f genz (x) = 0 if one of the x i > 0.5 exp( 5 i=1 x i ) otherwise (37) design matrix having M (n max = 10, d = 5) = 3003 columns, but ME-ACD leads to a far better accuracy with its ME procedure. Optimal degree of PC is equal to 9.

Conclusion and discussions

A new NI-ME method, ME-ACD, is developed in this paper. This method is able to efficiently partition the parameter space in order to construct a piecewise chaos polynomial model in the case of a discontinuous QOI, as the methods relying on Bayesian inference [START_REF] Sargsyan | Uncertainty quantification given discontinuous model response and a limited number of model runs[END_REF] or polynomial annihilation edge detection [START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF]. ME-ACD is also able to partition the space in presence of a discontinuity which runs partially through the parameter space, which was not possible with previous methods to the author knowledge. Moreover, the method permits to study a continuous but irregular QOI, which needs a lot of chaos modes to be accurately approximate by the standard PC models. ME-ACD model is composed of local models of lower degrees, taking the same logic as ME-PCM [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): Error analysis and applications[END_REF] but with a partition following the irregular behaviours. The method has been successfully tested on several discontinuous and continuous irregular QOIs from d = 2 to d = 5. From the results of the paper, the method is general and is very promising to approximate a large set of QOIs. The method will be tested on other kind of theoretical QOIs and on practical applications in the near future. On the majority of the discontinuous and irregular applications presented in the paper, ME-ACD models are more accurate than standard PC and GBT models. Moreover, ME-ACD can lead to physical interpretation, partitioning the space into elements which could be interpreted. With the use of a neural network classifier, the model is also able to determine which evaluation lie in a zone in between the given elements, and is so able to be cautious with these points, ignoring them instead of giving a value which risks to be inaccurate due to misclassification errors. ME-ACD has been developed in the context of UQ, but as other NI PC methods, it can be used in another context as a new regression method, with inputs which are not representing uncertaintites. ME-ACD presents limits which will be the subject of future researches. The partition performed by the ME-ACD algorithm is sensitive to the chosen hyperparameters and to the sampling strategy. In that way, the partition and the classification may lack of robustness according to the studied QOI and to the chosen sampling. One of the aims of ME-ACD is to work on every sampling. However, ME method may fail to give a suitable partition leading to an accurate model with the chosen sampling and/or hyperparameters. According to the QOI, some zones of the parameter space may need more observations than others. With a general sampling, sparsity can lead to selection of nearest neighbours from different sides of the discontinuity, leading to poorly accurate gradient estimation. Moreover, sparsity in some zones can lead to observations which are isolated from the others in D, not being enough to constitute a proper element with an accurate local model, like in the case of the QOI f bif -1000 studied in section 4.3, and therefore spoiling the accuracy of another element. Finally, this sparsity can alter the accuracy of the classification step, leading to missclassification errors undetected by the caution procedure. For now, even if these drawbacks were present, ME-ACD have been tested successfully on theoretical QOIs depending on a low or moderate number of input parameters, up to d = 5. For higher dimensions, the notion of nearest neighbours which is used for the gradient calculus and for the connectivity diagram, might become less "meaningful", as highlighted in [START_REF] Beyer | When is "nearest neighbor" meaningful?[END_REF]. ME-ACD is therefore expected to be even less robust as the dimension increases. Moreover, as the number of input parameters rises, the value of the QOI, which is one of the d + 1 considered coordinate in F or in D, is expected to have less impact in the selection of the nearest neighbours and in the agglomerative clustering step. The risk of selecting nearest neighbours from different sides of a discontinuity is therefore expected to increase, altering the calculus of the derivatives and then of the partition. First tests have been made with QOIs depending on d = 9 input parameters. The partition step worked on a simple constant piecewise QOI with d = 9, but failed with a more complex one where an element represented only 0.2% of the parameter space. These first tests showed that the method is promising with dimensions d ∈ [START_REF] Foo | The multi-element probabilistic collocation method (me-pcm): Error analysis and applications[END_REF][START_REF] John D Jakeman | Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions[END_REF], even if the lack of robustness persists and seems exacerbated. Further works will therefore consist in adapting the method for higher dimensions and finding a sampling strategy suitable to make the algorithm more robust, if the user has control over sampling. With such improvements, the use of local models with lower degrees would be promising to approximate higher dimensionality examples (d > 5). Finally, ME-ACD is also going to be adapted to approximate time-dependent functions being the solution of an ODE.

Figure 1 : 1 ( 1 (

 111 Figure 1: (a) Plot of f (1) 1 (dashed lines) and of its PC surrogate model f (1) 1 (plain lines) with n = 150. (b) Plot of f (1)

  Before introducing ME techniques, three NI methods are recalled. NI methods aim at finding the chaos coefficients of the model by running simulations on a set of observations called the training set or the training observations. In the whole article, the number of training observations is denoted N and the set of input parameters taken in the training observations is denoted T = {x i } N i=1 .

Figure 2 :

 2 Figure 2: N train = 500 observations generated by 2 sampling strategies. (a) Monte-Carlo sampling in X . (b) Mapping of the Monte-Carlo sampling in P. (c) Sobol sampling in X . (d) Mapping of the Sobol sampling in P.

Algorithm 2

 2 is the pseudo-code of the global design matrix computation with this method. Algorithm 2 Pseudo-code of the computation of the global design matrix d ← dimension of the parameter space T ← input parameters of training observations sampled in Step 1 n max ← maximal wanted degree of the local models

  This step is the core of ME-ACD. It consists in testing different partitions of the parameter space X using agglomerative clustering and to select the one providing the most accurate local models. Partitions withN E ∈ [|N E min , N Emax |]elements are considered leading to N p = N Emax -N E min + 1 tested partitions. In this section, as several partitions are considered, an element is denoted E p,e , where p ∈ [|1, N p |] is the index of the partition with N E = N E min + p -1 elements and e ∈ [|1, N E |] is the index of the element in this partition. E 1,e correspond to the elements of the partition with N E = N E min elements and E Np,e to the elements of the partition with N E = N Emax elements.

Algorithm 3

 3 Pseudo-code of Step 3 d ← dimension of the parameter space T ← input parameters of training observations sampled in Step 1 B ← values of the QOI of the training observations computed in Step 1 Initiate a N × d matrix D with zeros Construct the Ball tree with training data considered in F (using information of T and of B) J ← N × k max matrix where every row contains the indices of the k max nearest neighbours of the observation of index the index of the row sorted from the nearest to the furthest.

Figure 3 :

 3 Figure 3: Tree-like diagram giving an example of agglomerative clustering. Clusters are represented by circles and steps from N E = 1 to N Emax = 4 are displayed ; N E min = 2 in this example.

  highlights how the clustering is performed in space D (subfigures (b) and (d)), and how it permits the use of local models of lower degrees (as shown by subfigures (a) and (c)) in both examples.

Figure 4 :

 4 Figure 4: Plots of 50 observations of two 1D-QOIs f (1) 1 and f (1) 2 considered in the spaces F and D. The different colors of the observations stand for different clusters. (a) Observations of f (1) 1 represented in F. (b) Observations of f (1) 1 represented in D. (c) Observations of f (1) 2 represented in F. (d) Observations of f (1) 2 represented in D.

Figure 5 :

 5 Figure 5: Visualisation of 1000 labelled training observations represented in the 2D parameter space X partitioned for the 3 test functions. Different colors stands for different elements obtained with the clustering. (a) Visualisation of f hyperplane . (a) Visualisation of f hypersphere . (a) Visualisation of f multiple .

Figure 6

 6 Figure 6 gathers all convergence plots of the mean square error E, of the accuracy F of the ME-ACD classifier and of the fraction of kept observations by the caution procedure of ME-ACD on the 3 test

Figure 6 :

 6 Figure 6: On the left, convergence plots of the MSE between the values of the models and the true value of the 3 discontinuous functions on the test set. (a) f hyperplane . (c) f hypersphere . (e) f multiple . On the right, convergence plots of the accuracy of the classifier of ME-ACD and of the kept fraction of observations by the caution procedure of the ME-ACD, on the test set. (b) f hyperplane . (d) f hypersphere . (f) f multiple .

Figure 7 :

 7 Figure 7: Convergence plots of the MSE between the values of the models and the true value, on the N θ=0.85 of the of N test = 10 5 tests observations kept by the caution procedure set observations at each N . (a) Convergence plot of the MSE for f hyperplane on the kept observations. (b) Convergence plot of the MSE for f hypersphere on the kept observations. (c) Convergence plot of the MSE for f multiple on the kept observations.

Figure 8 :

 8 Figure 8: Evolution of the optimal degrees of PC models in the convergence studies for the 3 disconitnuous functions f hyperplane , f hypersphere and f multiple . The corresponding number of terms of the series M (n, d = 2) is shown on the left in parenthesis.

Figure 9 :

 9 Figure 9: (a) Convergence plot of the MSE for f partial on its test set. (b) Convergence plot of the accuracy of the ME-ACD classifier and fraction of kept observations by the caution procedure of ME-ACD for f partial on its test set. (c) Convergence plot of the MSE between the values of the models and the true value on the N θ=0.85 tests observations kept by the caution procedure at each N for f partial .

Figure 10 :

 10 Figure 10: Plot of f bif -t (x 1 , x 2 = 1) vs x 1 at 3 different instants : t = 1 (regular), t = 10 (irregular) and t = 1000 (discontinuous).

Figure 11 :Figure 12 :

 1112 Figure 11: On the left, convergence plots of the MSE between the values of the models and the true value of the 3 discontinuous functions on the test set. (a) f bif -1 . (c) f bif -10 . (e) f bif -1000 . On right, convergence plots of the kept fraction of observations by the caution procedure of the ME-ACD on the test set. (b) f bif -1 . (d) f bif -10 . (f) f bif -1000 .

4

  Pseudo-code of Step 4 d ← dimension of the parameter space N E min ← minimal number of elements N Emax ← maximal number of elements N p = N Emax -N E min + 1 δ min ← minimal partition error for a potential early stop of the search of partitions B ← values of the QOI on the training observations computed in Step 1 A ← global design matrix computed at Step 2 D ← derivatives of the QOI on the training observations computed at Step 3 L, E 1 , ..., E Np ← Results of ward clustering with N E min elements on the observations considered as points of D (using information of B and D) Initiate a variable δ opt to a huge value Initiate three void lists C, δ and done_elements for p ∈ [1, N p ] do N E = N E min + p -1 for e ∈ [1, N E ] do if the list e of E ← maximum of elements of δ corresponding to the elements of the partition p if δ p ≤ δ opt then N Eopt ← N E δ opt = δ p if δ p ≤ δ min then

	end if
	end for
	δ p Break the loops
	end if
	end if
	end for
	the coefficients of the local models are the coefficients in C corresponding to the elements of the partition
	with N Eopt elements

p , E p,e , is not already in done_elements then E p,e is appended to done_elements I ← indices of L where the items are equals to one of the labels contained in E p,e b ← subvector of B with the rows of indices in I Initiate a void list c and initiate a vector δ of size n maxn min + 1 with zeros for n ∈ [n min , n max ] do a ← submatrix of A with the M (n, d) first column and the rows of indices in I the result of the OLS with matrices a and b is appended to c and the related PRESS error is stored as δn-n min +1 end for n opt ← argmin n∈[n min ,nmax] δ c nopt-n min +1 is appended to C δ nopt-n min +1 is appended to δ

Table 1 :

 1 Hyperparameters and parameters selected for ME-ACD and PC for the convergence study of f hyperplane , f hypersphere and f multiple .

	3 test

Table 2 :

 2 Hyperparameters and parameters selected for ME-ACD and PC for the convergence study of f partial .

		PC	ME-ACD
	N	[100, 2000]
	N E min	1	1
	N Emax	1	10
	k min	not relevant	6
	k max	not relevant	10
	n min	0	0
	n max	20	10
	δ min	10 -5	10 -5
	learning rate	see section 4
	isConnect	False	False
	θ	not relevant	0.85

Table 3 :

 3 ). Hyperparameters and parameters selected for ME-ACD and PC for the convergence study of f bif -1 , f bif -10 and f bif -1000 .

		PC	ME-ACD
	N	[100, 2000]
	N E min	1	1
	N Emax	1	10
	k min	not relevant	6
	k max	not relevant	10
	n min	0	0
	n max	20	10
	δ min	10 -5	10 -5
	learning rate	see section 4
	isConnect	False	True
	k connect	not relevant	10
	θ	not relevant	0.85

Table 4

 4 gathers all selected hyperparameters and parameters, and Table5shows the performances of the 3 methods with N = 50000 training observations. Once again, ME-ACD is able to give a relatively accurate model and is able to reduce missclassification errors with the caution procedure, ignoring test observations located in the indecision zone. For this application, ME-ACD and PC use the same global

		PC	ME-ACD
	N	50000
	N E min	1	1
	N Emax	1	5
	k min	not relevant	15
	k max	not relevant	25
	n min	0	0
	n max	10	10
	δ min	10 -5	10 -5
	learning rate	see section 4
	isConnect	not relevant	True
	k connect	not relevant	30
	θ	not relevant	0.85

Table 4 :

 4 Hyperparameters and parameters selected for ME-ACD and PC for the convergence study of f genz

		PC	GBT	ME-ACD
	E	0.15415751 0.01737417 0.02062906
	E θ=0.85	not relevant	0.00182077
	V θ=0.85	not relevant	99.651%

Table 5 :

 5 Performances of the 3 methods for the 5d QOI f genz trained with N = 50000 observations. Performances are evaluated on a test set of N test = 10 5 test observations.

CRediT authorship contribution statement

Nicolas Vauchel: Writing, investigation, reviewing and editing. Éric Garnier: Supervision, Rewiewing. Thomas Gomez: Supervision, Rewiewing.

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.

Acknowledgements

The authors wish to thank ONERA and Hauts-de-France region for their funding.