S Marelli

P.-R Wagner

C Lataniotis

B Sudret

Stochastic Spectral Embedding

Keywords:

Constructing approximations that can accurately mimic the behavior of complex models at reduced computational costs is an important aspect of uncertainty quantification. Despite their flexibility and efficiency, classical surrogate models such as Kriging or polynomial chaos expansions tend to struggle with highly non-linear, localized or non-stationary computational models.

We hereby propose a novel sequential adaptive surrogate modeling method based on recursively embedding locally spectral expansions. It is achieved by means of disjoint recursive partitioning of the input domain, which consists in sequentially splitting the latter into smaller subdomains, and constructing a simpler local spectral expansions in each, exploiting the tradeoff complexity vs. locality. The resulting expansion, which we refer to as "stochastic spectral embedding" (SSE), is a piece-wise continuous approximation of the model response that shows promising approximation capabilities, and good scaling with both the problem dimension and the size of the training set.

We finally show how the method compares favorably against state-of-the-art sparse polynomial chaos expansions on a set of models with different complexity and input dimension.

Introduction

In the era of machine learning (ML) and uncertainty quantification (UQ), it is not surprising to see their boundary getting progressively blurred. Cross-fertilization between the two disci-1 plines is nowadays the norm, rather than an exception, and for good reasons. Physics-informed neural networks are reaching unprecedented approximation power in UQ applications (see, e.g., [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Pang | fPINNs: Fractional physics-informed neural networks[END_REF]), while sparse polynomial chaos expansions are used as denoising regressors in Torre et al. (2019a), as high-dimensional regression tools in [START_REF] Lataniotis | Extending classical surrogate modeling to high dimensions through supervised dimensionality reduction: a data-driven approach[END_REF], and on real-world experimental data in [START_REF] Abbiati | A global sensitivity analysis framework for hybrid simulation[END_REF]. UQ-born Gaussian process modeling [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]; [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] is now a staple tool in ML Rasmussen and Williams (2006), while support vector machines [START_REF] Vapnik | The nature of statistical learning theory[END_REF] found their way in rare event estimation [START_REF] Bourinet | Rare-event probability estimation with adaptive support vector regression surrogates[END_REF][START_REF] Moustapha | Comparative study of Kriging and support vector regression for structural engineering applications[END_REF].

More in general, the adoption of both surrogate models and ML is becoming mainstream in applied sciences and engineering, with applications in entirely different fields. A few examples from the recent literature include: macroeconomics [START_REF] Harenberg | Uncertainty quantification and global sensitivity analysis for economic models[END_REF], wind turbine modeling [START_REF] Slot | Surrogate model uncertainty in wind turbine reliability assessment[END_REF], nuclear engineering [START_REF] Radaideh | Analyzing nuclear reactor simulation data and uncertainty with the group method of data handling[END_REF], smart grid engineering [START_REF] Wang | Quantitative evaluations of uncertainties in multivariate operations of microgrids[END_REF], crash test simulations [START_REF] Moustapha | Comparative study of Kriging and support vector regression for structural engineering applications[END_REF], dam stability assessment [START_REF] Guo | Kriging based reliability and sensitivity analysis-application to the stability of an earth dam[END_REF], and many more. This list could be extended arbitrarily, as does the rich literature on these topics, but this task lies outside the scope of the current paper.

A common aspect across all of these works is the use of efficient and accurate functional approximation tools. Regardless of the specific technique, the general concept is straightforward:

given a finite set of input realizations and their corresponding model responses, known as the training set (ML) or experimental design (UQ), a suitable parametric function is calibrated such that it accurately approximates the underlying (possibly unknown) input-output map. For the sake of consistency, and a little bias towards UQ, we will refer to this process as surrogate modeling, acknowledging that it is also known as emulation, metamodeling, reduced order-or response surface-modeling, or sometimes simply regression. A variety of methods is available in the surrogate modeling literature, which we cluster here in two classes:

• Localized surrogates: this includes interpolants (e.g. Gaussian process modeling [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF], spline interpolation [START_REF] Reinsch | Smoothing spline functions[END_REF], sparse grids [START_REF] Bungartz | Sparse grids[END_REF]), but also local regression methods (e.g. Gaussian process regression [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], multivariate moving averages [START_REF] Lowry | A multivariate exponentially weighted moving average control chart[END_REF] or support vector machines [START_REF] Vapnik | The nature of statistical learning theory[END_REF]). These techniques rely on the availability of local information, e.g. through kernels on point-wise distance measures or support vectors, to provide predictions that are more accurate closer to the points in the training set. They therefore tend to perform better in interpolation, rather than extrapolation, tasks.

• Global surrogates: they provide global approximations without capitalizing on locally available information. Examples in this class include spectral methods (e.g. polynomial chaos expansions (PCE) [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]; [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] and Pointcaré expansions [START_REF] Roustant | Poincaré inequalities on intervals -application to sensitivity analysis[END_REF]), linear regression methods (e.g. compressive sensing [START_REF] Donoho | Compressed sensing[END_REF]; [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF], generalized linear models [START_REF] Nelder | Generalized linear models[END_REF]), artificial neural networks (ANNs) [START_REF] Goodfellow | Deep Learning[END_REF]. These techniques tend to achieve better global accuracy (e.g. in terms of generalization error), thus offering some degree of extrapolation capabilities, but also worse local accuracy than their localized counterparts.

Each of the two classes have advantages and disadvantages, but they both tend to perform well on models that show homogeneous complexity throughout the input parameter space. Some models of practical engineering relevance, however, can show a highly localized behavior in different regions of the parameter space. Common examples include likelihood functions used in Bayesian inference [START_REF] Nagel | Spectral likelihood expansions for Bayesian inference[END_REF], crash test simulations [START_REF] Serna | Advanced surrogate models for multidisciplinary design optimization. 8th Weimar Optimization and Stochastic Days[END_REF], snap-through models [START_REF] Hrinda | Snap-through instability patterns in truss structures[END_REF], and discontinuous models in general.

Different approaches with varying degree of complexity have been proposed in the UQ and ML literature to address this kind of behavior. Examples include regression trees [START_REF] Chipman | Bart: Bayesian additive regression trees[END_REF][START_REF] Breiman | Classification and regression trees[END_REF], multivariate adaptive regression splines (MARS, [START_REF] Friedman | Multivariate adaptive regression splines[END_REF]), various combinations of Kriging and PCE (PC-Kriging, [START_REF] Schöbi | Polynomial-chaos-based Kriging[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]), multi-resolution/multi-element polynomial chaos expansions [START_REF] Maître | Multi-resolution analysis of Wienertype uncertainty propagation schemes[END_REF]; [START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF]; [START_REF] Foo | The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications[END_REF] instead on the idea of partitioning the input parameter space into (often disjoint) subdomains, followed by the use of regression-based surrogates in each subdomain. This divide-and-conquer approach is particularly effective in reducing the complexity of the computational model in each subdomain, hence allowing relatively simple techniques to be used as local approximants. A prime example of this class of methods is given by regression trees [START_REF] Chipman | Bart: Bayesian additive regression trees[END_REF], where the local surrogates are as simple as constant values.

A common trait of most surrogate models used in a UQ context is that they rely on some form of regularity of the underlying computational model (e.g. smoothness or symmetry) to achieve an efficient representation based on an experimental design of relatively small size. It is therefore not surprising that they often show limited scalability with both the number of input dimensions (the well known curse of dimensionality) and with the size of the experimental design. Indeed, most local surrogates and interpolants rely on either kernel or clustering methods, neither of which scales linearly with the number of dimensions. Moreover, their training requires the solution of complex optimization problems that often have at least as many parameters as input dimension [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF][START_REF] Vapnik | The nature of statistical learning theory[END_REF]. Global regression methods, on the other hand, require the optimization and storage of a large number of parameters or coefficients, which also rarely scales linearly in high dimension for non-trivial models.

To step further into scalability considerations, the number of available samples in the experimental design deserves some discussion. Historically, UQ-based surrogate modeling has taken a parsimonious approach: focus on small but informative experimental designs (N ED ≈ 10 1-2 ), because of the high computational costs associated to engineering models, and to their smooth behavior. On the other hand, ML has seen its expansion in the era of big data, focusing on large experimental designs (N ED ≈ 10 5-7 ), with often noisy data and highly non-smooth behavior. Albeit the gap is closing over time, a no-man's land in between the two still exists: computational models that show a complexity that is too high for classical surrogate modeling (e.g. extremely non-linear, or highly localized), but are expensive enough to only allow for N ED ≈ 10 3-4 , regardless of the input dimension.

It is with this class of problems in mind that we propose a new surrogate modeling technique, namely stochastic spectral embedding (SSE), that combines global spectral representations and adaptive domain decompositions. We demonstrate that SSE can efficiently approximate models with varying degrees of complexity across the input space, while maintaining favorable scaling properties with both the input dimension and the size of the experimental design.

The paper is organized as follows: we first describe the general rationale and the details of the algorithm in Section 2. Then, in Section 3 we tackle the issue of constructing an SSE from an experimental design, in a regression context. In Section 4, we choose a reference spectral decomposition technique (polynomial chaos expansions, PCE) and we apply SSE to two highly complex analytical functions to showcase its capability to adapt to models with non-homogeneous complexity, and its scalability to high dimensions and large experimental designs. Finally, we also tackle two models of engineering complexity that are known to be challenging for classical surrogate modeling methods. We present concluding remarks in Section 5 and discuss extensions of the algorithm that could further improve its performance.

Stochastic spectral embedding: rationale and main algorithm

As the name suggests, stochastic spectral embedding (SSE) is a combination of two main ingredients: a stochastic spectral representation-based surrogate model and some form of embedding, which implies the sequential construction of subdomains of the full input space. In other words, SSE consists in iteratively refining a spectral surrogate model by means of embedding additional surrogate models in subdomains of the parent expansion. In a sense, SSE can be seen as an extension of regression trees [START_REF] Breiman | Classification and regression trees[END_REF]; [START_REF] Chipman | Bart: Bayesian additive regression trees[END_REF] to a much wider class of regression models, with the addition of a strong stochastic component due to the use of spectral representations.

Spectral expansions

We will consider herein the Hilbert space H of random variables of the form Y = M(X) with finite second moments (E Y 2 < ∞), where X is an M -dimensional random vector with joint distribution X ∼ f X (x). Let the space be equipped with the inner product:

g(X), h(X) H def = E [g(X)h(X)] = DX g(x)h(x)f X (x) dx, (1)
where D X ⊆ R M is the support of X. Then, every Y ∈ H admits a spectral representation M S of the form:

Y = M S (X) def = ∞ j=1 a j Ψ j (X), (2)
where the a j ∈ R are real coefficients, and the Ψ j 's form a countably infinite orthonormal basis of the space:

E [Ψ i (X)Ψ j (X)] = Ψ i (X), Ψ j (X) H = δ ij , (3)
where δ ij is the Kronecker delta. For notational simplicity, the inner product subscript H is omitted hereinafter.

Spectral decompositions of the form of Eq. ( 2) have a property that is particularly important for surrogate modelling, namely the fact that due to the orthogonality of the basis in Eq. (3), their (finite) second moment is given by:

E M(X) 2 = M S (X), M S (X) = ∞ j=1 a 2 j < +∞. (4)
The converging sum in Eq. ( 4) implies therefore that the coefficients a j must decay at least geometrically when sorted by decreasing absolute value. This property is sometimes referred to as compressibility, because it essentially means that most of the information on the model variability is contained in a finite set of coefficients/basis elements. This allows one to truncate the spectral decomposition in Eq. (2) even if in principle it has an infinite number of terms. The truncated version of Eq. ( 2) is given by:

M S (X) ≈ M S (X) = j∈A a j Ψ j (X), ( 5
)
where A is a truncation set (typically related to the complexity of the basis functions, e.g., maximum frequency in Fourier expansions, or maximum polynomal degree in PCE). The rapid decay in the coefficients of spectral expansions is the main reason why many powerful surrogate modeling techniques that belong to the so-called class of compressive sensing [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | An introduction to compressive sampling: A sensing/sampling paradigm that goes against the common knowledge in data acquisition[END_REF], have proven to be very effective in various recent applications [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]Torre et al., 2019a;[START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF]. Compressive sensing uses sparse regression tools to identify the best truncation set A based on the available information in the experimental design.

Because of the truncation introduced in Eq. ( 5), the expansion is in general not exact, hence we define the residual R(X) as:

R(X) = M(X) -M S (X). ( 6
)
Due to the convergent behavior of the truncated expansion, it follows that Var [R(X)]

Var M S (X) . By definition, spectral expansions belong to the class of global representations, i.e. the basis functions in Eq. ( 5) have support on the entire domain D X . Therefore, highly localized models, or those with inhomogeneous behavior throughout the input domain tend to require an extremely large number of terms in the truncated expansion to achieve satisfactory approximation accuracy (Gibbs phenomenon). As an example, the number of terms in the wellestablished polynomial chaos expansion [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] can grow very fast when the underlying model has strongly localized behavior, because a high polynomial degree is required for an accurate representation.

A sequential partitioning approach

To alleviate this limitation, while still capitalizing on the powerful convergence properties of spectral methods, SSE constructs a sequence of spectral expansions of manageable complexity on increasingly smaller subdomains of the original domain. Such subdomains are denoted D ,p X ⊆ D X , where is the expansion level and p is a subdomain index within the level. The expansion is performed only on the local residual from the previous level.

For illustration purposes, Figure 1 shows an example of sequential partitioning for a simple 2D bounded domain, obtained by splitting each subdomain in two equal subdomains across a random dimension. When = 0, there is only a single subdomain D 0,1 X def = D X and the residual is R 0,1 S (X) = M S (X) from Eq. ( 6).

While the idea of partitioning the input space in smaller subdomains and constructing a surrogate model in each is certainly not new in UQ (see, e.g. [START_REF] Maître | Multi-resolution analysis of Wienertype uncertainty propagation schemes[END_REF]; [START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF]), the use of a sequence of residuals in SSE sets it apart from other divide and conquer methods.

In more formal terms, SSE in its basic form can be written as a multi-level expansion of the form:

(a) = 0 (b) = 1 (c) = 2
M SSE (X) = L =0 P p=1 1 D ,p X (X) R ,p S (X), ( 7
)
where L is the total number of expansion levels considered, P is the number of subdomains at

level ∈ {1, . . . , L}, 1 D ,p X (X) is the indicator function of the subdomain D ,p X . Finally, R ,p S (X) is the truncated expansion of the residual of the SSE up to level -1, R (X), R (X) = M(X) - -1 k=0 P k p=1 1 D k,p X (X) R k,p S (X). ( 8
)
In the above equations, each residual term is expanded onto a local orthonormal basis as follows:

R k,p S (X) = j∈A k,p a k,p j Ψ k,p j (X). (9)
A local inner product is defined in the domain D k,p X :

Ψ k,p i (X), Ψ k,p j (X) k,p = D k,p X Ψ k,p i (x)Ψ k,p j (x)f k,p X (x) dx, (10)
where:

f k,p X (x) = 1 D k,p X (x) f X (x) V k,p (11)
is the joint PDF of the input parameters restricted to the subdomain D k,p X and rescaled by its probability mass V k,p :

V k,p = D k,p X f X (x) dx. ( 12
)
A crucial aspect of SSE is that by partitioning the entire input domain D X into smaller subdomains D ,p X , it trades the complexity of the single, often global expansion in Eq. ( 5) for a (possibly large) number of local expansions with much smaller truncation sets. In cases where the spectral basis is continuous in Eq. ( 5), SSE results in a final piecewise continuous approximation, but no continuity is ensured on the boundaries of the subdomains. Mean-square convergence of the procedure is guaranteed by the spectral convergence in each level, which implies that the residual local variance in each subdomain is in expectation decreasing rapidly. In other words, for each increasing level in Eq. ( 7), new discontinuity bounds are generated during the partitioning step, but the variance of the overall residual is reduced, thus resulting, in expectation, in lower amplitude discontinuities. This behavior is analogous to that of regression trees [START_REF] Friedman | Multivariate adaptive regression splines[END_REF]; [START_REF] Breiman | Classification and regression trees[END_REF].

The SSE algorithm

Algorithmically, SSE consists of a local refinement sequence of a global spectral expansion into sequentially smaller subdomains D ,p X . For notational simplicity, we introduce here a set of local random vectors distributed according to the local PDF in Eq. ( 11): X ,p ∼ f ,p X (x). We further choose a certain partitioning strategy that is discussed in Section 3.2.

Then, the SSE algorithm can be written as:

1. Initialization:

(a) = 0, p = 1 (b) D ,p X = D X (c) R (X) = M(X) 2. For each subdomain D ,p X , p = 1, • • • , P : (a) Calculate the truncated expansion R ,p S (X ,p ) of the residual R (X ,p ) in the current subdomain (b) Update the residual in the current subdomain R +1 (X ,p ) = R (X ,p ) -R ,p S (X ,p ) (c) Split the current subdomain D ,p X in N S subdomains D +1,{s1,••• ,sN S } X
based on a partitioning strategy (d) If < L, ← + 1, go back to 2a, otherwise terminate the algorithm 3. Termination (a) Return the full sequence of D ,p X and R ,p S (X ,p ) needed to compute Eq. ( 7).

Note that in steps 2a and 2b of the previous algorithm the residual R (X ,p ) is only indexed by , but not by the subdomain index p. This is because the residual is fully defined with respect to the previous level , which is independent on the particular subdomain under consideration (see Figure 1).

Building a stochastic spectral embedding from data

For it be useful in practical applications, SSE needs to be "trainable" from a finite-size experimental design. Hereinafter, we consider an experimental design N ) and its corresponding model evaluations Y = y (1) , • • • , y (N ) as the only data available for training.

X = x (1) , • • • , x (
Upon closer inspection of the algorithm in Section 2.3, the training phase of the SSE representation consists in estimating the following quantities from the available experimental design:

1. The expansion coefficients of the local residual in each level and subdomain a ,p (Eq. ( 9)).

2. A partitioning strategy at each level.

3. The total number of splitting levels, L.

In the following sections we introduce a comprehensive adaptive strategy based on sparse linear regression to perform each of these steps from a given experimental design.

Calculating the residual expansion coefficients

For a specific subdomain D ,p X , a local spectral expansion of the residual R S needs to be constructed from the available experimental design. We therefore define a so-called local experimental design X ,p ⊆ X , the subset of the original experimental design lying within the subdomain

D ,p X : X ,p def = x (j) , j = 1, • • • , N ,p , such that x (j) ∈ X ∩ D ,p X . (13)
A similar notation is used to identify the corresponding model responses, Y ,p . Using the auxiliary local random vector X ,p introduced in the previous section, the residual expansion in Eq. ( 9) reads:

R ,p S (X ,p ) = j∈A ,p a ,p j Ψ ,p j (X ,p ). ( 14
)
Given the local experimental design X ,p and a truncated local spectral basis Ψ ,p j , j ∈ A ,p , the task of identifying the coefficients a ,p def = a ,p j , j ∈ A ,p can then be cast as a linear regression problem (see, e.g., [START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF]):

a ,p ≈ a ,p = arg min a x (i) ∈X ,p   R (x (i) ) - j∈A ,p a ,p j Ψ ,p j (x (i) )   2 . ( 15
)
While in principle the regression problem in Eq. ( 15) can be solved through ordinary least squares, recent literature on the topic of compressive sensing has amply demonstrated that sparse regression approaches can provide great benefits in terms of accuracy, especially for relatively small experimental designs [START_REF] Donoho | Compressed sensing[END_REF]; [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]; [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF]. A review of the available techniques for this purpose lies outside the scope of this paper and is extensively explored for one popular class of spectral representations (polynomial chaos expansions) in [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF].

Partitioning strategy

A second step necessary to construct SSE from data is to identify a proper partitioning strategy between levels. Any strategy for the partitioning of the input domain D X can be employed for Eq. ( 7), under the sole condition that at each level :

P p=1 D ,p X = D X . ( 16
)
While a comprehensive study on different partitioning strategies would be interesting, for the sake of simplicity we adopt hereinafter a rather simple approach, similar in spirit to regression trees [START_REF] Breiman | Classification and regression trees[END_REF]. In other words, we split every subdomain in two parts of equal probability mass along one of the input directions, d ,p ∈ {1, • • • , M }. Note that the direction in itself can be different for each subdomain, even on the same level.

Under very general conditions, it is possible to bijectively map any random vector X with joint distribution F X , to the uniform independent random vector U ∼ U(0, 1) M through an appropriate isoprobabilistic transform (e.g., the Rosenblatt transform Rosenblatt (1952); Torre et al. ( 2019b)):

U = g(X) X = g -1 (U ), (17)
where g(•) denotes the isoprobabilistic transform. This mapping simplifies the proposed partitioning strategy: splitting is performed in the uniformly distributed quantile space U , and the resulting split domains D ,p U are mapped back to the input space X via the inverse transform (see Eq. ( 17)). This has several computational benefits, including proper treatment of unbounded variables. Figure 2 shows graphically a two-dimensional example of partitioning in the quantile (uniform) space U , and its corresponding mapping to unbounded random variables in the physical space X.

In the general case, a strategy is needed to choose a specific splitting direction d ,p ∈ {1, • • • , M } for each existing subdomain D ,p X . Different heuristic reasoning can be used to make this choice, including purely random splitting (as in Figure 2), using the direction of maximum residual difference, or estimates of the variability of the R ,p S in each direction (following the same rationale as in [START_REF] Shields | Adaptive monte carlo analysis for strongly nonlinear stochastic systems[END_REF]). The optimal criterion can be application-specific, because it may in principle depend on the chosen spectral representation.

Sparse tree representation and expansion truncation

In a regression context, it is difficult to choose a priori a truncation on the maximum number of levels L in the expansion in Eq. ( 7). Because the samples may be unequally distributed, some subdomains at each level may be empty, or more formally X ,p = ∅ for some combinations of the two probability spaces from Eq. ( 17). The splitting direction in each subdomain is determined randomly in this example.

and p. We therefore take a straightforward approach to obviate this issue, by initializing the residual at every level and subdomain to the null function, hence a ,p = 0. The coefficients are then updated only in the subdomains that satisfy X ,p ≥ N min , where N min is a parameter of the SSE algorithm that represents the minimum number of points in a subdomain required to justify an expansion. The SSE expansion is truncated when no new updates are possible given the current experimental design X , or more formally:

L = min : X ,p < N min , ∀ p ∈ {1, • • • , P } -1. ( 18
)
In addition to providing a suitable stopping criterion for the algorithm in Section 2.3, an added benefit of this strategy is that only the residual expansions that were effectively updated need to be stored in memory. This provides a degree of sparsity in the representation and potentially significantly reduces the memory fingerprint of the method, especially in the case of a large number of points in the experimental design.

Note that, for an experimental design of size N and a minimum number of points per expansion N min , the following holds:

2 LN min ≤ N, ( 19
)
where L is the expected value of the maximum L in (18) and therefore

L ≤ log 2 N N min . ( 20
)

Error estimation

In the context of surrogate modeling, assessing the accuracy of the approximation is an important task. Arguably the best known accuracy estimator in function approximation is the so-called generalization error E GEN , which for SSE is given by

E GEN def = E (M(X) -M SSE (X)) 2 . (21)
A direct estimation of this quantity is in general impossible, as it would require the availability of an extensive validation set. Instead, because we adopt a regression approach to calibrate the spectral decompositions in each subdomain, we estimate the generalization error through leave-one-out cross-validation [START_REF] Chapelle | Model selection for small sample regression[END_REF][START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] that is available for each of the local expansions.

For notational convenience, we introduce here the set of terminal domains

D T = D L,1 X , • • • , D L,PL X , i.e.
E ,p GEN = E R (X ,p ) -R ,p S (X ,p ) 2 , (22)
then the global generalization error is simply given by the average error in each terminal domain:

E GEN = E E L,p GEN = PL p=1 E L,p GEN • V L,p , (23)
which is equal to the sum of the terminal domain errors weighted by the corresponding probability mass in Eq. ( 12).

To provide an estimator based on the available experimental design X , we only need an estimator of E L,p GEN . Arguably the most common tool for the estimation of generalization error in regression problems is k-fold cross-validation (see, e.g., [START_REF] Vapnik | The nature of statistical learning theory[END_REF]). The special case of k = N is also known as leave-one-out error and marked E LOO . For ordinary least square regression it can be calculated analytically from the expansion coefficients and the basis functions [START_REF] Chapelle | Model selection for small sample regression[END_REF][START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF].

Given the sparse tree representation described in Section 3.3, it cannot be guaranteed that the residual R is expanded in every terminal domain. Therefore, during the splitting phase of the SSE algorithm (Step 2c of the algorithm in Section 2.3) we initialize the error of all the subdomains D +1,{s1,••• ,sN S } of the current subdomain D ,p to the leave-one-out error of the latter E ,p LOO :

E +1,{s1,••• ,sN S } LOO = E ,P LOO . (24)
Then, we update the error estimate in each subdomain during Step 2a only if the conditions for its expansion hold (see Section 3.3). As a result, every terminal domain is either assigned its own leave-one-out error if it contains a residual expansion, or inherits the leave-one-out error from the last ancestor domain that was expanded.

By using the leave-one-out error in each terminal domain as an estimator of its generalization error E L,p GEN = E L,p LOO , the empirical estimator of Eq. ( 23) reads:

E GEN = PL p=1 E L,p LOO • V L,p . (25)
In most metamodeling applications, it is customary to normalize the estimated error by the variance of the experimental design, to obtain a dimensionless error measure. The relative error is thus defined as:

GEN = 1 Var [Y] PL p=1 E L,p LOO • V L,p . ( 26
)
4 Applications

In this section we aim at showing the performance of SSE on a set of applications that can prove challenging for standard metamodeling techniques. Because of its widespread use in the uncertainty quantification of engineering models, we choose as a spectral decomposition technique polynomial chaos expansions [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]; [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] (hereinafter PCE). This choice is also quite convenient due to several specific properties of PCE, that combine well with SSE.

Synergies with polynomial chaos expansions

By using the same notation as in Section 2, and assuming that X has independent components, the truncated polynomial chaos expansion of a finite variance model can be written as Le [START_REF] Gratiet | Metamodel-based sensitivity analysis: polynomial chaos expansions and Gaussian processes, Chapter 38, Handbook on Uncertainty Quantification[END_REF]:

M PCE (X) = α∈A a α Ψ α (X), ( 27
)
where α is a multi-index that identifies the polynomial degree in each variable, A is a suitable truncation set (e.g. A = A M,d containing all multivariate polynomials with degree ≤ d), and the Ψ α (X) form an orthogonal basis of multivariate polynomials. The latter can be obtained via tensor product of univariate polynomials as follows:

Ψ α (X) = M i=1 Φ (i) αi (X i ), (28)
where

Φ (i)
αi is a polynomial of degree α i that belongs to the family of univariate polynomials orthogonal with respect to the input PDF of X i ∼ f Xi (x i ) and the inner product in Eq. ( 1).

An interesting property of the univariate polynomials that synergizes well with our proposed SSE, is that it is possible to construct polynomials orthogonal to almost any input PDF through Gram-Schmidt orthogonalization (for an extensive review, see [START_REF] Gautschi | Orthogonal polynomials: computation and approximation[END_REF]; [START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF]). In the context of SSE, this property has a powerful implication: in each subdomain D ,p the basis elements Ψ ,p α (X ,p ) in Eq. ( 9) are still polynomial functions of the original input variables X.

This property, together with the analytical integrability of polynomials, allows us to derive several statistics of interest of PCE-based SSE analytically. Let us first introduce the notion of flattened representation: because SSE is a polynomial in the original variables in every level and subdomain, this also holds for the terminal domains introduced in Section 3.4. Therefore, one can project the full SSE in Eq. ( 7) as a local PCE onto each terminal domain:

M F SSE (X L,p ) = α∈A T c p α Ψ L,p α (X L,p ), ( 29
)
where A T is a suitable truncation set for the projection to be exact, and the c p α are the corresponding coefficients. The latter can easily be computed either analytically or exactly through quadrature. Note that, while the basis elements in the PCE in Eq. ( 29) correspond to the Ψ k,P j in Eq. ( 9) (they only depend on the input PDF in Eq. ( 11)), in general the coefficients will not be the same, i.e. c p α = a L,p α in Eq. ( 9).

Because PCE contains as a basis element the constant term, it is straightforward to demonstrate that the expected value of Eq. ( 7) reads:

E [M SSE (X)] = PL p=1 c p 0 V L,p , (30)
which is the weighted mean of all the mean values of the flattened representation in Eq. ( 29).

Similarly, the variance can be calculated as:

Var [M SSE (X)] = PL p=1 V L,p α∈A T (c p α ) 2 -E [M SSE (X)] 2 . ( 31
)
A number of other quantities of engineering interest (e.g. conditional variances, Sobol' sensitivity indices, etc.) can be derived similarly from the flattened representation. A selection of those is reported in A.

From a technical perspective, the flattened representation in Eq. ( 29) contains all the information needed to evaluate Eq. ( 7) on a new point, but at a much lower storage cost, as only the final sets of coefficients c p α and basis indices A T need to be stored. This has additional advantages during the prediction of the response on new points, because it only requires the prediction of a single local expansion in the appropriate terminal domains, rather than that of all of its ancestors as in the original formulation in Eq. ( 7). More formally, for a point x 0 ∈ D X it is sufficient to find p 0 ∈ 1, • • • , P L for which x 0 ∈ D L,p0 X and evaluate the flattened SSE from Eq. ( 29) for p = p 0 .

Example applications and testing strategy

To compare the performance of SSE over sparse PCE, we choose four reference problems of increasing complexity: (i) a one-dimensional analytical function with localized non-polynomial behavior, (ii) a 100-dimensional analytical function with decreasing parametric importance in higher dimensions, (iii) an 8-dimensional engineering model describing the performance function of a damped oscillator and (iv) a three-dimensional discontinuous engineering model describing the snap-trough behavior of a truss structure.

Among the ingredients identified in Section 2.3 is a partitioning strategy, to choose the splitting direction in every subdomain. After extensive testing, we found that splitting according to the direction of highest variability of R ,p S proved to be the most effective, especially for smaller experimental designs. We therefore split each subdomain D ,p X into two subdomains with equal probability mass, i.e. N S = 2, along the direction that has the maximum first order Sobol' indexSobol' (1993), as analytically derived from the coefficients of R ,p S [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. In all applications we compare the convergence behavior of SSE vs. its spectral counterpart PCE as a function of the experimental design (ED) sizes N ED . To assess the robustness of the results, we consider 10 independent replications of each ED, and provide the results in Tukey boxplots. For each experimental design size, SSE construction is terminated for N min = min{5M, 50} (see Section 3.3).

As a spectral technique, we adopt the adaptive sparse-PCE based on LARS approach developed in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] in its numerical implementation in UQLab (Marelli andSudret, 2014, 2019). Each R ,p S is therefore a degree-and q-norm-adaptive polynomial chaos expansion. We further introduce a rank truncation of r = 2 to cope with high dimensional problems, i.e. we limit the maximum number of input interactions [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF] to 2 variables at a time. The truncation set for each spectral expansion (Eq. ( 27)) thus reads:

A M,p,q,r = {α ∈ N M : ||α|| q ≤ p, ||α|| 0 ≤ r}, (32)
where

||α|| q = M i=1 α q i 1 q , q ∈ (0, 1]; ||α|| 0 = M i=1 1 {αi>0} . (33)
The q-norm is adaptively increased between q = {0.5, • • • , 0.8} while the maximum polynomial degree is adaptively increased in the interval p = {0, 1, • • • , p SSE max }, where the maximum degree p SSE max is a parameter for each case study. In all examples the SSE performance is compared to standard polynomial chaos expansions on the same ED. These PCEs are constructed with the same adaptive approach used for the SSE expansions. Their maximum degree, however, is denoted by p PCE max . It is set to the highest value our computational budget admitted for a given dimensionality but at least to p PCE max > 2 • p SSE max . We compare the performance of SSE and PCE in terms of the relative mean squared error (MSE) η, a well known error estimator defined as

η def = E M(X) -M(X)

Application 1: one-dimensional analytical function

We first present a simple one-dimensional example that is meant to illustrate how SSE behaves with a model on which PCE is expected to fail. The model is given by:

M(X) = -X + 0.1 sin (30X) + exp (-(50(X -0.65)) 2 ), ( 35
)
where X ∼ U(0, 1) is a uniformly distributed random variable. The first two terms in Eq. ( 35)

(polynomial and sinusoidal) can be accurately approximated by a low degree PCE, while the third term (squared exponential) causes PCE to require extremely high degree due to the localized peak it introduces at x = 0.65 (see Figure 3). In the same X : [0, 0.5], because the available maximum degree p SSE max is not sufficiently high, resulting in a mostly constant polynomial correction. In subdomain D 1,2

X : [0.5, 1], the same problem is observed and the insufficient maximum degree results only in a small global improvement. In the next step (Figure c), the residual in D 2,1 X , D 2,2 X and D 2,4 X is significantly reduced to a very small oscillation around 0. After the final step (Figure d), the overall approximation is has a high accuracy.

From the residual progression it can be seen that the algorithm needs more levels to accurately approximate the target function near regions of high complexity, i.e., near the exponential peak.

While this property does not affect the convergence when an experimental design of fixed size is chosen, it can be exploited in adaptive experimental design settings [START_REF] Wagner | Bayesian model calibration with stochastic spectral embedding[END_REF].

As expected, the final SSE accuracy increases with the size of the experimental design. In the SSE relative MSE is at least one order of magnitude smaller than the PCE error. There is considerable variability in the relative MSE between individual realizations. This can be attributed to the squared exponential peak in Eq. ( 35): depending on the input realizations in the experimental design, it is captured better or worse by the available data.

Application 2: 100-dimensional analytical function

With this example we want to explore the scalability of SSE in high dimensional problems. This example uses a variant of a test function introduced in Zhou (1998). We modified the function to have a high nominal dimensionality (M = 100), and relatively low effective dimensionality meaning that the majority of the variability can be attributed to a small number of input parameters. It takes the form

M(X) = 10 M 2 [ϕ(10 • (X -1/3)) + ϕ(10 • (X -2/3))] ,
where ϕ(x)

def = (2π) -M/2 exp - 1 2 M i=1 a 2 i x 2 i . (36)
The factor a = (a 1 , • • • , a M ) modifies the original function and serves as a dimensiondependent weight that decays exponentially, as:

a i = e -(i-1) , with i ∈ {1, • • • , M }. (37)
The input random vector X is distributed according to a multivariate standard uniform distribution with independent marginals, X ∼ f X (X) = Two contour cross-section plots of this function are shown between dimensions {X 1 , X 2 } in To manage the computational complexity, in this example we limit the maximum polynomial degree in SSE to p SSE max = 2 and compute the SSE with total experimental design sizes of N ED = {1,000; 2,000; 5,000; 10,000}. For PCE, we choose a maximum degree of p SSE max = 7, which is the maximum degree we could run on a standard desktop computer with 16GB of RAM before incurring memory issues. The resulting comparison between PCE and SSE with respect to the relative MSE is plotted in Figure 6. In this scenario, the SSE algorithm outperforms PCE on all investigated experimental designs. In fact, PCE seems to benefit from increasing the experimental design only marginally, with a relative error of η ≈ 0.25 for all considered values of N ED . SSE shows instead a convincing convergence behavior, by reducing its residual by almost two orders of magnitude across the various experimental designs. Given the high dimensionality of this problem, the approximation power of sparse PCE is limited by the curse of dimensionality, rather than the lack of data. By reducing the complexity of the spectral representation at each level, SSE can better exploit informative datasets without incurring similar computational bottlenecks.

Application 3: damped oscillator

Damped oscillators are a class of engineering models that is commonly used in structural reliability problems [START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF]. This class of problems is known to be often difficult to surrogate, due its high non-linearity and often local behavior. A sketch of the oscillator considered in this example is displayed in Figure 7. It consists of a primary and secondary system with masses m p , m s , stiffnesses k p , k s and damping ratios ζ p , ζ s . The subscripts p and s denote the primary and secondary system properties, respectively.

In this example we consider the limit state function of the damped oscillator given by

M(X) = F s -p • k s E S [x 2 S ], (38)
where F s is the force capacity of the secondary spring, p is the so-called peak factor and x S is the relative displacement between the primary and secondary systems. The mean-square relative displacement of the secondary spring under a white noise base acceleration S is analytically given by:

E S x 2 S = π S 0 4ζ s ω 3 s ζ a ζ s ζ p ζ s (4ζ 2 a + θ 2 ) + γζ 2 a (ζ p ω 3 p + ζ s ω 3 s )ω p 4ζ a ω 4 a , (39)
where S 0 is the white noise intensity, ω p = k p /m p and ω s = k s /m s are the natural frequen- All variables but the peak factor (set to p = 3) are modelled as independent random variables and are summarized in the random vector X = {m p , m s , k p , k s , ζ p , ζ s , S 0 , F s }. Their marginal distributions are lognormal, with the parameters given in Table 1.

For the convergence study, we choose experimental design sizes of N ED = {1,000; 5,000; 10,000; 20,000}.

The maximum degrees are set to p SSE max = 4 and p PCE max = 10 for SSE and PCE, respectively. Figure 8 summarizes the results. This benchmark is known to be quite difficult to approximate with standard surrogate model techniques, as it is clear from the scale of the RMSE in Figure 8. For the two smaller experimental design sizes N ED = {1,000; 5,000}, both PCE and SSE perform quite poorly, with SSE showing a similar median behavior, but much higher variability. For larger experimental designs, however, the SSE performance improves significantly over that of PCE, until at N ED = 20,000 the relative MSE of SSE is half that of PCE. This behavior is in line with the previous findings: provided enough information, SSE can provide higher expressive power than its static counterpart.

Application 4: truss with discontinuous snap-through behavior

As a last example, we address another problem of engineering interest: the geometrically nonlinear two-bar truss structure shown in Figure 9. The structure itself is defined by the initial inclination α 0 and length 0 of the two bars. A peculiarity of this structure is that it exhibits the so-called snap-through behavior. At first, the vertical displacement w of such a structure typically increases linearly with an increasing load P (Figure a). Once a specific critical load is exceeded, the structure snaps through to another equilibrium point, at which the load can be increased further (Figure b). The main implication of this kind of behavior is that it is discontinuous in some critical regions of the input space, which are in general unknown a priori.

The vertical displacement w of the truss tip is related to the angle α by

M(X) = w = 0 cos α 0 (tan α 0 -tan α(X)). (40)
At the same time, α needs to satisfy the following constitutive equation that depends on the random vector X = {P, E, A}:

P = -2EA tan α(cos α 0 -cos α). (41)
For a given realization of X, this equation can be solved numerically for α, the value of which then is used in Eq. ( 40) to estimate the corresponding vertical displacement.

In this study we set the constants l 0 = 5 m and α 0 = 10 • , and treat the parameters X as independent random variables with marginals listed in Table 2 [START_REF] Moustapha | A two-stage surrogate modelling approach for the approximation of models with non-smooth outputs[END_REF].

We investigate experimental designs of sizes N ED = {100; 500; 1,000; 2,000; 5,000} and set the maximum polynomial degrees to p SSE max = 4 and p PCE max = 10 for SSE and PCE, respectively. Figure 10 summarizes the convergence behavior of PCE and SSE in this benchmark. For all experimental designs, SSE outperforms sparse PCE. The dispersion of the RMSE is also significantly improved. These observations can be explained with the well-known Gibbs phenomenon in spectral representations, that leads to large discrepancies close to discontinuities. The effect is far less severe (although still present) for SSE than for PCE because it is restricted to those subdomains that actually contain the discontinuity. This behavior is investigated more closely in Figure 11, where a cross section through M is shown. It is created by setting A to its mean value and drawing a map proportional to the point-wise discrepancy in the remaining directions.

We adjust transparency of the model response to reflect the underlying joint PDF: solid colors correspond to high probability, fading ones to low probability. LOO from Eq. ( 25). All plots clearly show the effect of the function discontinuity. For PCE the Gibbs phenomenon is clearly visible. It causes a large error near the discontinuity, with an oscillating error at large distances to the discontinuity. Naturally, SSE also suffers from the same problem, but its effect is only localized close to the discontinuity and it does not affect further regions. Furthermore, the available domain-wise local error estimator in Figure c gives a clear indication of local loss of accuracy of SSE. In practical applications this information can be crucial, as it allows one to assess the confidence of the surrogate model predictions, and could be used to adaptively enrich the experimental design close to critical regions [START_REF] Wagner | Bayesian model calibration with stochastic spectral embedding[END_REF].

Error estimation accuracy

For practical applications, it is important that surrogate models offer insights into their prediction accuracy. Most surrogate models offer global confidence bounds on their predictions that are typically computed through cross-validation techniques (e.g., leave-one-out error), while some techniques also offer point-wise confidence bounds (e.g., Kriging Santner et al. (2003) and bootstrap PCE [START_REF] Marelli | An active-learning algorithm that combines sparse polynomial chaos expansions and bootstrap for structural reliability analysis[END_REF]).

For SSE, the domain-wise error estimators of the local expansions give some insight into the local accuracy of SSE, as shown in the last case study (Section 4.6). The weighted sum of those domain-wise estimators can be used as a global estimator of the generalization error as proposed in Section 3.4, Eq. ( 26). To assess the accuracy of this global estimator, in Figure 12 we plot it against the relative MSE on a validation set for all presented case studies. The diagonal dashed line corresponds to perfect error estimation: η = ε GEN . In all the applications ε GEN significantly underestimates the true error. For comparison, we show the corresponding LOO estimator from PCE, which also exhibits a bias towards lower errors. On the other hand, it is clear that there is a strong correlation the SSE estimators and the true validation error across all applications, which still makes ε GEN a strong potential candidate for model selection in future extensions of the method to adaptively select the SSE hyperparameters (e.g., adaptively choosing p SSE max ).

Considerations on computational costs and scalability

For a metamodeling technique to be relevant in an engineering context, the associated training and prediction costs need to be negligible with respect to the costs of the underlying physicsbased computational model. In this sense, SSE performs quite well, as its computational costs scale primarily with the size of the experimental design. Indeed, because the residual expansions (Eq. ( 9)) can be chosen with very low degree or otherwise strict truncation, the driving cost is the total number of subdomains that are expanded in the SSE sparse tree (see Section 3.3).

Interestingly, the expected number of expansions N E depends only on the ratio between the total experimental design size N and the minimum number of points in each subdomain needed to perform a residual expansion, N min , as in (see also Eq. ( 20)):

N PCE = L =0 2 = 2 L+1 -1 = 2 log 2 (N/Nmin) +1 -1 ≤ 2N/N min -1. ( 42
)
In other words, the computational complexity increases at most linearly with the number of points in the experimental design. Additionally, the storage costs can be further reduced to 2 L = N/N min expansions when using the flattened representation in Eq. ( 29). Therefore, as required, the training and evaluation costs of SSE are normally negligible with respect to those needed to produce a training set for any realistic engineering application.

Conclusions

In an effort to extend the applicability of the powerful class of spectral decomposition-based metamodels, we propose a novel metamodel technique called stochastic spectral embedding (SSE), that exploits both recent advances in UQ (sparse spectral expansions) and in machine learning (regression trees). While our presentation was general in nature, we showed how well this approach synergizes with sparse polynomial chaos expansions. We also provided analytical formulas to calculate several statistical properties of the resulting model by means of the so-called flattened representation, which has additional benefits in terms of computational costs.

We tested the performance of SSE on both simple test functions and engineering-like examples of varying dimensionality and complexity, using varying experimental design sizes, and compared it to our best performing sparse PCE. Its generalization capabilities, especially for highly complex models and large experimental designs, outperform PCE in most cases.

We also demonstrated that the associated computational costs of the training of SSE scale linearly in expectation with the number of points in the experimental design. This compares favorably with most metamodeling techniques common in the UQ community (e.g. PCE or Kriging).

This performance, however, comes at the cost of trading the continuity of PCE for the piecewise continuity of SSE. This also implies the loss of the effective generalization error estimate provided by LOO in linear regression. To mitigate this issue, we proposed the error estimate GEN (Eq. ( 26)). Despite its absolute scale being biased towards lower values, it still shows high correlation with the actual generalization error for all experimental design sizes and dimensions. This is a promising property for further research into providing automatic selection of the hyperparameters of the algorithm (which at the moment are the maximum degree of the residual expansions p SSE max , as well as the minimum number of points in each subdomain N min required to expand the residual) and further enhance its performance.

Additional research is ongoing towards the use of the approximate local error measures pro-vided by E ,p GEN (Eq. ( 22)) for goal-oriented adaptive experimental design construction, a topic explored by the authors in [START_REF] Wagner | Bayesian model calibration with stochastic spectral embedding[END_REF]. in the full input space are not necessarily terminal subdomains in the lower dimensional, conditional expectation input space defined by u. However, because the basis functions are polynomials, it is once again possible to perform the flattening process (see Section 3.4) and obtain disjoint subdomains. By denoting as P ⊆ {1, • • • , P L } the set of terminal subdomains in the conditional input variables X u , we can rewrite Eq. ( 44) as

E [M SSE (X)| X u ] = p∈P α∈A P v=0 d p α Ψ L,p α,u (X L,p u ), (45)
where A P v=0 is a suitable multi-index set allowing an exact representation of the polynomials and d p α are the corresponding coefficients.

A.2 Partial variance and Sobol' indices

The Sobol' Hoeffding decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Gratiet | Metamodel-based sensitivity analysis: polynomial chaos expansions and Gaussian processes, Chapter 38, Handbook on Uncertainty Quantification[END_REF] of the SSE representation M SSE reads

M SSE (X) = M 0 SSE + u⊂{1,••• ,M } u =∅ M u SSE (X u ), ( 46
)
where M 0 SSE def = E [M SSE (X)] and the remaining terms can be computed recursively by

M i SSE (X i ) = E X∼i [M SSE (X)] -M 0 SSE , (47)
M ij SSE (X ij ) = E X∼ij [M SSE (X)] -M i SSE (X i ) -M j SSE (X j ) -M 0 SSE , (48)
• • • = • • • (49) M u SSE (X u ) = E Xv [M SSE (X)] - w⊂u w =∅ M w SSE (X w ) -M 0 SSE . ( 50
)
The decomposition of Eq. ( 46), allows the definition of the so-called partial variance, i.e., the fraction of the variance Var [M SSE (X)] that can be attributed to X u , defined by

V u def = DX u (M u SSE (x u )) 2 f Xu (x u ) dx u . (51)
Using Eq. ( 47), the so-called first order partial variance can therefore be written as

V i = DX i E X∼i [M SSE (X)] -M 0 SSE 2 f Xi (x i ) dx i . ( 52
)
With the expression for the conditional expectation from Eq. ( 45), this integral can be solved analytically as

V i = DX i   p∈P α∈A P ∼i=0 d p α Φ (i),L,p αi (X L,p i ) -M 0 SSE   2 f Xi (x i ) dx i (53) = p∈P V L,p i α∈A P ∼i=0 (d p α ) 2 -M 0 SSE 2 . ( 54
)

  and deep neural networks Goodfellow et al. (2016), among others. Such methods can be broadly classified in two macro-families: global approximations with local refinements (e.g. PC-Kriging), or domain-decomposition-based methods (regression trees, MARS, multi-element polynomial chaos expansions). The class of global approximations with local refinements rely on efficiently combining global surrogates (e.g. spectral decompositions as polynomial chaos expansions, or global regression models) with local interpolation techniques (e.g. Gaussian processes or splines), to provide surrogates with acceptable extrapolation capabilities and good local accuracy. The class of domain-decomposition-based methods relies

Figure 1 :

 1 Figure 1: Example of a possible partitioning sequence for an SSE with L = 2 and P = 2 , and bounded uniform marginal distributions. In this simple example, each subdomain is split in two equal parts across a random direction.

Figure 2 :

 2 Figure 2: Graphical representation of the partitioning strategy described in section 3.2 for a twodimensional problem with independent random variables. Upper row: partitioning in the quantile space; Lower row: partitioning in the original space. Red dots show a random sampling from the original distributions in both spaces, and serve as a visual aid to recognize the mapping between

  those domains that belong to the last expansion level L in Eq. (7). By definition D T is a complete partition of the input domain D X . Because of the sequential nature of SSE, which locally refines the previous approximation level with the expansion of the residual in the current subdomain, an accurate estimate of the local generalization error in each of the terminal domains would then suffice to provide an estimate of the overall E GEN of the full SSE. If we denote the local residual error:

  Figure we detail four SSE refinement stages, with N ED = 200 and p SSE max = 5. For every step we show on the top panel a graphical representation of the various subdomains identified by the algorithm, with the subdomains of level highlighted in orange. In the middle panel we plot the true model (orange solid line) and the current SSE approximation M SSE as a dashed blue line. In the bottom panel, we plot the corresponding residual M(X) -M SSE (X) as a solid blue line in the same vertical scale as in the middle panel, for comparison. In the first step in Figure a the main trend of the function is identified, leaving a residual that mainly consists of the sine oscillation and the exponential peak. In the following step (Figure b) the approximation is not greatly improved in the subdomain D 1,1

  Figure 4, we compare SSE and PCE on a set of experimental design sizes of N ED = {10, 50, 100, 200} in terms of their relative mean squared error (MSE, Eq. (34)). PCE is constructed with a maximum adaptive degree of p PCE max = 20. At the extremely small experimental design of N ED = 10, the SSE approach is comparable to PCE. As the available experimental design points increase, SSE exhibits faster convergence in RMSE than PCE, and from N ED = 50 onwards SSE consistently outperforms PCE in this problem. At larger experimental designs, SSE can accurately reproduce the localized behavior of this test function, while not being constrained by the global nature of PCE basis functions defined on the full domain. At the final ED size of N ED = 200

Figure 3 :

 3 Figure 3: One-dimensional analytical function: selected steps of the SSE construction and the resulting domains, residuals and total approximation. The terminal domains (Eq. (29)) are highlighted in orange. 18

Figure 4 :

 4 Figure 4: One-dimensional analytical function: comparison of RMSE convergence between PCE and SSE as a function of the number of points in the experimental design. A slight horizontal offset is added to improve readability.

Figure a and

  Figure a and between {X 1 , X 10 } in Figure b. The effect of the decay factor a is clearly visible and results in a close-to constant behavior along the parametric dimension X 10 .

Figure 6 :

 6 Figure 6: 100-dimensional analytical function: comparison of RMSE convergence between PCE and SSE as a function of the number of points in the experimental design. A slight horizontal offset is added to improve readability.

Figure 7 :

 7 Figure 7: 8-dimensional damped oscillator : model setup.

  cies of the two subsystems, and the further abbreviations are used: γ = m s /m p , ω a = (ω p +ω s )/2, ζ a = (ζ p + ζ s )/2 and θ = (ω pω s )/ω a .

Figure 8 :

 8 Figure 8: 8-dimensional damped oscillator : comparison of RMSE convergence between PCE and SSE as a function of the number of points in the experimental design. A slight horizontal offset is added to improve readability.

Figure 9 :

 9 Figure 9: 3-dimensional snap through truss: illustration of the truss structure subject to snapthrough

  Figures a and b show the relative

Figure 10 :

 10 Figure 10: 3-dimensional snap-through truss: comparison of RMSE convergence between PCE and SSE as a function of the number of points in the experimental design. A slight horizontal offset is added to improve readability.

Figure c shows

  Figure c shows instead the domain-wise error estimator E ,pLOO from Eq. (25). All plots clearly show the effect of the function discontinuity. For PCE the Gibbs phenomenon

  (a) PCE point-wise error: (M(x) -MPCE(x)) 2 /Var [M] (b) SSE point-wise error: (M(x) -MSSE(x)) 2 /Var [M] (c) SSE domain-wise error estimator: E ,p LOO /Var [M]

Figure 11 :

 11 Figure 11: 3-dimensional snap-through truss: comparison between PCE/SSE pointwise error and SSE domain-wise error estimator from Eq. (25)) normalized by the model variance. The fine lines show the SSE subdomains.

Figure 12 :

 12 Figure 12: Comparison of the relative generalization error GEN with the relative generalization error estimator GEN which for SSE is given by Eq. (26) and for PCE is the relative leave-one-out error E LOO /Var [Y ].

Table 1 :

 1 8-dimensional damped oscillator : marginal distributions.

	Variable Description	Distribution Mean C.O.V.
	m p	primary mass	Lognormal	1.50	0.1
	m s	secondary mass	Lognormal	0.01	0.1
	k p	primary spring stiffness	Lognormal	1.00	0.2
	k s	secondary spring stiffness	Lognormal	0.01	0.2
	ζ p	primary damping ratio	Lognormal	0.05	0.4
	ζ s	secondary damping ratio	Lognormal	0.02	0.5
	S 0	white noise intensity	Lognormal	100.00	0.1
	F s	secondary spring force capacity Lognormal	15.00	0.1

Var [M(X)] (34)where M is either the PCE or SSE surrogate model. This error measure was estimated with standard Monte Carlo simulations using a large sample of size N = 10 6 .

i=1 U(0, 1).

Acknowledgements

The PhD thesis of the second author is supported by ETH grant #44 17-1.

A Postprocessing PCE-based SSE

If polynomials chaos expansions are used to construct the residual expansions, several quantities of engineering interest can be computed analytically as a post-processing step of the final SSE.

In this section we derive expressions for (i) conditional expectations, (ii) partial variances and (iii) Sobol' indices.

A.1 Conditional expectations

Conditional expectations describe the expectation of a multivariate function of a random vector

where we used the flattened representation from Eq. ( 29). This corresponds to marginalizing over the parameters X v . Due to the local orthonormality of the SSE representation, an analytical expression for this integral can be found as

where

u is an auxiliary random variable that is only defined in the (L, p)-subdomain and Ψ L,p α,u (X L,p u )

) is a polynomial basis function of the non-marginalized variables.

The marginalization process in Eq. ( 43) creates an additional problem: this expression generally contains overlapping subdomains {D L,p Xu } p=1,••• ,PL due to the fact that terminal subdomains With the availability of the partial variances in Eq. ( 52) and the total variance from Eq. ( 31), one can analytically derive the first order Sobol' indices:

Higher order indices can be computed in a similar way.