Blatman, G. and Sudret, B. (2010a). An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. *Prob. Eng. Mech.*, 25:183-197.

Blatman, G. and Sudret, B. (2010b). Efficient computation of global sensitivity indices using sparse polynomial chaos expansions. *Reliab. Eng. Sys. Safety*, 95:1216-1229.

Blatman, G. and Sudret, B. (2011). Adaptive sparse polynomial chaos expansion based on Least Angle Regression. *J. Comput. Phys*, 230:2345-2367.

Bruckstein, A. M., Donoho, D. L., and Elad, M. (2009). From sparse solutions of systems of equations to sparse modeling of signals and images. *SIAM Rev.*, 51(1):34-81.

Candès, E. J. and Plan, Y. (2011). A probabilistic and RIPless theory of compressed sensing. *IEEE Trans. Inform. Theory*, 57(11):7235-7254.

Candès, E. J. and Wakin, M. B. (2008). An introduction to compressive sampling: A sensing/sampling paradigm that goes against the common knowledge in data acquisition. *IEEE Signal Process. Mag.*, 25(2):21-30.

Chapelle, O., Vapnik, V., and Bengio, Y. (2002). Model selection for small sample regression. *Machine Learning*, 48(1):9-23.

Chatterjee, T., Chakraborty, S., and Chowdhury, R. (2019). A critical review of surrogate assisted robust design optimization. *Arch. Comput. Method. E.*, 26(1):245-274.

Diaz, P. (2018). DOPT_PCE. https://github.com/CU-UQ/DOPT_PCE. [Online; accessed 14-May-2020].

Diaz, P., Doostan, A., and Hampton, J. (2018). Sparse polynomial chaos expansions via compressed sensing and D-optimal design. *Comput. Methods Appl. Mech. Engrg.*, 336:640-666.

Doostan, A. and Owhadi, H. (2011). A non-adapted sparse approximation of PDEs with stochastic inputs. *J. Comput. Phys.*, 230(8):3015-3034.

Dubourg, V. (2011). *Adaptive surrogate models for reliability analysis and reliability-based design optimization*. PhD thesis, Université Blaise Pascal, Clermont-Ferrand, France.

Echard, B., Gayton, N., Lemaire, M., and Relun, N. (2013). A combined importance sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models. *Reliab. Eng. Syst. Saf.*, 111:232-240.

Efron, B., Hastie, T., Johnstone, I., and Tibshirani, R. (2004). Least angle regression. *Ann. Stat.*, 32:407-499.

Ernst, O., Mugler, A., Starkloff, H.-J., and Ullmann, E. (2012). On the convergence of generalized polynomial chaos expansions. *ESAIM: Math. Model. Numer. Anal.*, 46(02):317-339.

Fajraoui, N., Marelli, S., and Sudret, B. (2017). Sequential design of experiment for sparse polynomial chaos expansions. *SIAM/ASA J. Uncertain. Quantif.*, 5(1):1061-1085.
Forrester, A., Sobester, A., and Keane, A. (2008). *Engineering design via surrogate modelling: a practical guide*. Wiley.

Gerstner, T. and Griebel, M. (2003). Dimension-adaptive tensor-product quadrature. *Computing*, 71:65–87.

Guo, X., Dias, D., and Pan, Q. (2019). Probabilistic stability analysis of an embankment dam considering soil spatial variability. *Comput. Geotech.*, 113:103093.

Hampton, J. and Doostan, A. (2017). BASE_PC. https://github.com/CU-UQ/BASE_PC. [Online; accessed 14-May-2020].

Hampton, J. and Doostan, A. (2018). Basis adaptive sample efficient polynomial chaos (BASE-PC). *J. Comput. Phys.*, 371:20–49.

Hariri-Ardebili, M. and Sudret, B. (2020). Polynomial chaos expansion for uncertainty quantification of dam engineering problems. *Eng. Struct.*, 203(#109631).

Jakeman, J. D., Eldred, M. S., and Sargsyan, K. (2015). Enhancing $l_1$-minimization estimates of polynomial chaos expansions using basis selection. *J. Comput. Phys.*, 289:18–34.

Konakli, K. and Sudret, B. (2016). Global sensitivity analysis using low-rank tensor approximations. *Reliab. Eng. Sys. Safety*, 156:64–83.

Le Gratiet, L., Marelli, S., and Sudret, B. (2017). *Metamodel-based sensitivity analysis: polynomial chaos expansions and Gaussian processes*, chapter 38, pages 1289–1325. Springer International Publishing. Cham, Switzerland.

Loukrezis, D., Römer, U., and De Gersem, H. (2019). Assessing the performance of Leja and Clenshaw-Curtis collocation for computational electromagnetics with random input data. *Int. J. Uncertainty Quantification*, 9(1).

Lüthen, N., Marelli, S., and Sudret, B. (2021). Sparse polynomial chaos expansions: Literature survey and benchmark. *SIAM/ASA J. Uncertain. Quantif.*, 9(2):593–649.

Mai, C. V. and Sudret, B. (2015). Hierarchical adaptive polynomial chaos expansions. In Papadrakakis, M., Papadopoulos, V., and Stefanou, G., editors, *1st Int. Conf. on Uncertainty Quantification in Computational Sciences and Engineering (UNCECOMP)*, Creta, Greece.

Marelli, S., Lüthen, N., and Sudret, B. (2021). UQLab user manual – Polynomial chaos expansions. Technical report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich, Switzerland. Report # UQLab-V1.4-104.

Marelli, S. and Sudret, B. (2014). UQLab: A framework for uncertainty quantification in Matlab. In *Vulnerability, Uncertainty, and Risk (Proc. 2nd Int. Conf. on Vulnerability, Risk Analysis and Management (ICVRAM2014)), Liverpool, United Kingdom)*, pages 2554–2563.

Morris, M. D., Mitchell, T. J., and Ylvisaker, D. (1993). Bayesian design and analysis of computer experiments: use of derivatives in surface prediction. *Technometrics*, 35:243–255.
Narayan, A., Gittelson, C., and Xiu, D. (2014). A stochastic collocation algorithm with multi-fidelity models. *SIAM J. Sci. Comput.*, 36(2):A495–A521.

Pati, Y. C., Rezaiifar, R., and Krishnaprasad, P. S. (1993). Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In *Proc. of 27th Asilomar Conf. on Signals, Systems and Computers*, pages 40–44. IEEE.

Sagi, O. and Rokach, L. (2018). Ensemble learning: A survey. *Wires. Data Min. Knowl.*, 8(4):e1249.

Slot, R. M., Sorensen, J. D., Sudret, B., Svenningsen, L., and Thogersen, M. L. (2020). Surrogate model uncertainty in wind turbine reliability assessment. *Renewable Energy*, 151:1150–1162.

Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. *Reliab. Eng. Sys. Saf.*, 93:964–979.

Thapa, M., Mulani, S. B., and Walters, R. W. (2020). Adaptive weighted least-squares polynomial chaos expansion with basis adaptivity and sequential adaptive sampling. *Comput. Meth. Appl. Mech. Eng.*, 360:112759.

van den Berg, E. and Friedlander, M. P. (2008). Probing the Pareto frontier for basis pursuit solutions. *SIAM J. Sci. Comput.*, 31(2):890–912.

Van den Berg, E. and Friedlander, M. P. (2015). SPGL1 - A Matlab solver for sparse optimization. https://friedlander.io/software/spgl1/. [Online; accessed 28-August-2019].

Vapnik, V. N. (1995). *The Nature of Statistical Learning Theory*. Springer-Verlag, New York.

Xiu, D. and Karniadakis, G. E. (2002). The Wiener-Askey polynomial chaos for stochastic differential equations. *SIAM J. Sci. Comput.*, 24(2):619–644.
# Appendix

## 4.A Details on the settings for the basis adaptivity schemes in the benchmark

In Table 4.A.1, we list the basis adaptivity settings for each of the 11 models in our benchmark.

* The rule for the static basis (ST) is to choose a total-degree basis with $p$ so that the number of basis functions $P$ is closest to $\frac{10}{3}N$, where $N$ is the number of ED points. Here, for low-dimensional models, the q-norm is chosen as $q = 1$, while for high-dimensional models, we use $q = 0.5$.

* For degree and q-norm basis adaptivity (PQ), we choose the ranges for degree and q-norm as large as possible while still keeping the number of basis functions computationally manageable (rule of thumb $\lesssim 10^4$).

* For forward neighbor degree adaptivity (FN), the degree of the initial basis is chosen so that the size of the basis is closest to $10N$ (as recommended in (Jakeman et al., 2015)), while we set the q-norm to the maximum of the q-norm-range for PQ basis adaptivity.

* Finally, for anisotropic degree basis adaptivity, which is only used for low-dimensional models, we use $q = 1$ and $p = \lceil \frac{p_{max}}{2} \rceil$, where $p_{max}$ is the maximum of the degree range for PQ basis adaptivity.

## 4.B Additional results

### 4.B.1 Basis adaptivity benchmark: raw data

Figs. 4.B.1 and 4.B.2 show the raw data of the basis adaptivity benchmark: for each model and ED size, we display the boxplots of resulting relative MSE for each combination of sparse solver and basis-adaptive scheme. The star-shaped markers visualize the attained error for one specific ED realization. They illustrate that the resulting errors are not always well correlated: on the same ED, some methods achieve one of their smallest errors while others produce one of their largest. Therefore, we continue in Section 4.3.1 by investigating the relative performance of the different methods.
<table><caption>**Table 4.A.1:** Details on the initial bases and degree and q-norm ranges for the various basis adaptivity schemes</caption><thead><tr><th>Model</th><th>dim d</th><th>static basis</th><th>PQ range</th><th>FN initial</th><th>AD initial</th></tr></thead><tbody><tr><td>Ishigami function</td><td>3</td><td>p = 8 (small ED)<br/> / p = 12 (large),<br/> q = 1</td><td>p ∈ [1, ..., 25],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 12 (small ED)<br/> / p = 19 (large)</td><td>p = 13</td></tr><tr><td>Undamped oscillator</td><td>6</td><td>p = 4 / p = 4,<br/> q = 1</td><td>p ∈ [1, ..., 10],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 5 / p = 6</td><td>p = 5</td></tr><tr><td>Borehole function</td><td>8</td><td>p = 4 / p = 4,<br/> q = 1</td><td>p ∈ [1, ..., 10],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 5 / p = 6</td><td>p = 5</td></tr><tr><td>Damped oscillator</td><td>8</td><td>p = 4 / p = 5,<br/> q = 1</td><td>p ∈ [1, ..., 7],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 5 / p = 6</td><td>p = 4</td></tr><tr><td>Wingweight function</td><td>10</td><td>p = 3 / p = 4,<br/> q = 1</td><td>p ∈ [1, ..., 7],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 4 / p = 5</td><td>p = 4</td></tr><tr><td>Truss model</td><td>10</td><td>p = 3 / p = 4,<br/> q = 1</td><td>p ∈ [1, ..., 6],<br/> q ∈ [0.5, 0.6, ..., 1]</td><td>p = 4 / p = 5</td><td>p = 3</td></tr><tr><td>Morris function</td><td>20</td><td>p = 6 / p = 8,<br/> q = 0.5</td><td>p ∈ [1, ..., 8],<br/> q ∈ [0.4, 0.5, 0.6]</td><td>p = 6 / p = 8,<br/> q = 0.6</td><td>—</td></tr><tr><td>Structural frame model</td><td>21</td><td>p = 5 / p = 8,<br/> q = 0.5</td><td>p ∈ [1, ..., 8],<br/> q ∈ [0.4, 0.5, 0.6]</td><td>p = 6 / p = 8,<br/> q = 0.6</td><td>—</td></tr><tr><td>2-dim diffusion model</td><td>53</td><td>p = 3 / p = 4,<br/> q = 0.5</td><td>p ∈ [1, ..., 6],<br/> q ∈ [0.4, 0.5, 0.6]</td><td>p = 4 / p = 5,<br/> q = 0.6</td><td>—</td></tr><tr><td>1-dim diffusion model</td><td>62</td><td>p = 3 / p = 4,<br/> q = 0.5</td><td>p ∈ [1, ..., 5],<br/> q ∈ [0.4, 0.5, 0.6]</td><td>p = 3 / p = 4,<br/> q = 0.6</td><td>—</td></tr><tr><td>100D function</td><td>100</td><td>p = 3 / p = 4,<br/> q = 0.5</td><td>p ∈ [1, ..., 5],<br/> q = 0.5</td><td>p = 4 / p = 5,<br/> q = 0.5</td><td>—</td></tr></tbody></table>
<figure><img src="image_58.png" /><figcaption>(a) Ishigami function, $N = 50$</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>(b) Ishigami function, $N = 150$</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>(c) Undamped oscillator, $N = 60$</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>(d) Undamped oscillator, $N = 120$</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>(e) Borehole function, $N = 100$</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>(f) Borehole function, $N = 250$</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>(g) Damped oscillator, $N = 150$</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(h) Damped oscillator, $N = 350$</figcaption></figure>

**Figure 4.B.1:** Comparison of different combinations of solvers and basis adaptivity schemes for **low-dimensional models**. We show validation errors attained by all 24 combinations of methods on 30 realizations of experimental design. The boxplots visualize the attained errors across all 30 realizations. The star-shaped markers denote the attained errors of all combinations for one selected ED realization. They highlight one set of data points which is also part of the larger set visualized by the boxplots. Left: small ED; right: large ED. Abbreviations of basis-adaptive schemes: ST: static basis; PQ: degree and q-norm; FN: forward neighbor; AD: anisotropic degree. Continued in next figure.
<figure><img src="image_66.png" /><figcaption>(i) Truss model, $N = 100$</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>(j) Truss model, $N = 250$</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>(k) Wingweight function, $N = 100$</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(l) Wingweight function, $N = 250$</figcaption></figure>

**Figure 4.B.1:** Continued. Comparison of different combinations of solvers and basis adaptivity schemes for **low-dimensional models**. We show validation errors attained by all 24 combinations of methods on 30 realizations of experimental design. The boxplots visualize the attained errors across all 30 realizations. The star-shaped markers denote the attained errors of all combinations for one selected ED realization. They highlight one set of data points which is also part of the larger set visualized by the boxplots. Left: small ED; right: large ED. Abbreviations of basis-adaptive schemes: ST: static basis; PQ: degree and q-norm; FN: forward neighbor; AD: anisotropic degree.

<figure><img src="image_70.png" /><figcaption>(m) Morris function, $N = 150$</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>(n) Morris function, $N = 350$</figcaption></figure>

**Figure 4.B.2:** Comparison of different combinations of solvers and basis adaptivity schemes for **high-dimensional models**. We show validation errors attained by all 15 combinations of methods on 30 realizations of experimental design. The boxplots visualize the attained errors across all 30 realizations. The star-shaped markers denote the attained errors of all combinations for one selected ED realization. They highlight one set of data points which is also part of the larger set visualized by the boxplots. Left: small ED; right: large ED. Abbreviations of basis-adaptive schemes: ST: static basis; PQ: degree and q-norm; FN: forward neighbor. Continued in next figure.
<figure><img src="image_72.png" /><figcaption>(o) Structural frame, $N = 150$</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(p) Structural frame, $N = 350$</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(q) Diffusion 2D, $N = 100$</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(r) Diffusion 2D, $N = 400$</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>(s) Diffusion 1D, $N = 100$</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(t) Diffusion 1D, $N = 400$</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(u) 100D function, $N = 400$</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(v) 100D function, $N = 1200$</figcaption></figure>

**Figure 4.B.2:** Continued. Comparison of different combinations of solvers and basis adaptivity schemes for **high-dimensional models**. We show validation errors attained by all 15 combinations of methods on 30 realizations of experimental design. The boxplots visualize the attained errors across all 30 realizations. The star-shaped markers denote the attained errors of all combinations for one selected ED realization. They highlight one set of data points which is also part of the larger set visualized by the boxplots. Left: small ED; right: large ED. Abbreviations of basis-adaptive schemes: ST: static basis; PQ: degree and q-norm; FN: forward neighbor; AD: anisotropic degree.
4.B.2 Bar plots for all combinations of methods

In Figs. 4.B.3-4.B.6, we display the aggregated results for all combinations of solvers and sam-
pling schemes (while in Section 4.3.1, we only displayed the six best combinations). We sort the
method combinations according to how often they achieved an error within two times the best
relMSE (bright red bar). For a detailed description of how to read this plot, we refer to Section
4.3.1.

<figure><img src="image_80.png" /><figcaption>**Figure 4.B.3: Small ED sizes.** Aggregated results for low-dim models.</figcaption></figure>
<figure><img src="image_81.png" /><figcaption>**Figure 4.B.4: Small ED sizes.** Aggregated results for **high-dim** models.</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>**Figure 4.B.5: Large ED sizes.** Aggregated results for **low-dim** models.</figcaption></figure>
<figure><img src="image_83.png" /><figcaption>**Figure 4.B.6: Large ED sizes.** Aggregated results for **high-dim** models.</figcaption></figure>
CHAPTER 5

Global sensitivity analysis using
derivative-based sparse Poincaré
chaos expansions

This chapter contains the preprint of

Lüthen, N., Roustant, O., Gamboa, F., Iooss, B., Marelli, S., Sudret, B. (2022).

Global sensitivity analysis using derivative-based sparse Poincaré chaos expansions.

https://arxiv.org/abs/2107.00394.

submitted for publication in the *International Journal of Uncertainty Quantification* in March 2022.

**Author contributions.** N. Lüthen: Conceptualization, Methodology, Software, Validation, Investigation, Writing - Original Draft, Visualization. O. Roustant: Conceptualization, Methodology, Formal Analysis, Writing - Original Draft, Supervision. F. Gamboa: Conceptualization, Methodology, Formal Analysis, Writing - Original Draft. B. Iooss: Conceptualization, Methodology, Formal Analysis, Writing - Original Draft. S. Marelli: Writing - Review & Editing, Supervision. B. Sudret: Conceptualization, Writing - Review & Editing, Supervision.

Abstract

Variance-based global sensitivity analysis, in particular Sobol' analysis, is widely used for de-
termining the importance of input variables to a computational model. Sobol' indices can
be computed cheaply based on spectral methods like polynomial chaos expansions (PCE). An-
other choice are the recently developed Poincaré chaos expansions (PoinCE), whose orthonormal
tensor-product basis is generated from the eigenfunctions of one-dimensional Poincaré differen-
tial operators. In this paper, we show that the Poincaré basis is the unique orthonormal basis
with the property that partial derivatives of the basis form again an orthogonal basis with re-
spect to the same measure as the original basis. This special property makes PoinCE ideally
suited for incorporating derivative information into the surrogate modelling process. Assuming
that partial derivative evaluations of the computational model are available, we compute spec-
tral expansions in terms of Poincaré basis functions or basis partial derivatives, respectively, by
sparse regression. We show on two numerical examples that the derivative-based expansions
provide accurate estimates for Sobol' indices, even outperforming PCE in terms of bias and
variance. In addition, we derive an analytical expression based on the PoinCE coefficients for a
second popular sensitivity index, the derivative-based sensitivity measure (DGSM), and explore
its performance as upper bound to the corresponding total Sobol' indices.

5.1 Introduction

Computer models simulating physical phenomena and industrial systems are commonly used
in engineering and safety studies, for prediction, validation or optimisation purposes. These
numerical models often take as inputs a high number of physical parameters, whose values
are variable or not perfectly known, creating the need for uncertainty quantification on model
computations (Smith, 2014). Uncertainty quantification typically becomes more challenging the
higher the input dimension is (curse of dimensionality). In this situation, global sensitivity
analysis (GSA) is an invaluable tool that allows the analyst to rank the relative importance
of each input of the model and to detect non-influential inputs (Borgonovo and Plischke, 2016;
Razavi et al., 2021). Most often relying on a probabilistic modeling of the model input variables,
GSA tries to explain model output uncertainties on the basis of model input uncertainties,
accounting for the full range of variation of the variables.

A well-known and widely used GSA method is Sobol' analysis (Sobol', 1993), which relies on the
functional ANOVA (analysis of variance) decomposition (Efron and Stein, 1981). For a square-
integrable model and independent input variables, Sobol' analysis determines which part of the
model output variance can be attributed to each input and to each interaction between inputs.
The overall contribution of each input, including interactions with other inputs, is provided by
the total Sobol' index (Homma and Saltelli, 1996). Sobol' indices can be estimated efficiently
using various Monte Carlo-based techniques as well as metamodel-based techniques (Prieur and
Tarantola, 2017). The latter save on expensive model evaluations by first performing a small
number of model runs, which are used to compute an accurate approximation to the original
model – the meta- or surrogate model – from which the Sobol' indices are finally computed
(Fang et al., 2006; Le Gratiet et al., 2017).

One of the most popular and powerful metamodelling methods is the polynomial chaos expansion
(PCE) (Xiu and Karniadakis, 2002). PCE represents the model in a specific basis consisting of
polynomials that are orthonormal with respect to the input distribution. Orthogonal polynomial
systems have been studied throughout the last century and they have many useful properties (see,
e.g., Szegö (1939) and Simon (2010)). One particular strength of PCE is that once it is computed,
it easily gives all the variance-based quantities defined through the ANOVA decomposition, and
in particular the Sobol' indices at all orders (Sudret, 2008). In practice, the expansion cannot
use infinitely many terms and must be truncated. Among the many approaches available to
compute the expansion coefficients, sparse regression techniques combined with adaptive basis
selection appear to be especially promising (see Lüthen et al. (2021, 2022a) for an overview).
Here, a small number of terms is selected which is able to best represent the computational
model based on the available model evaluations.

In some practical situations, partial derivatives of the model output with respect to each input
are easily accessible, for example by algorithmic differentiation of the numerical model in the
reverse (adjoint) mode (Griewank and Walther, 2008). This technique allows for computing all
partial derivatives of the model output at a cost independent of the number of input variables.
Since PCEs are such a well-established metamodelling tool, there have been many efforts to
leverage the additional information contained in model derivatives to improve the performance
of PCE. The idea of including derivative information into sparse regression problems, often
called *gradient-enhanced *ℓ*¹-*minimization*, is tested by Jakeman et al. (2015) for one numerical
example with uniform inputs, and analyzed theoretically and numerically by Peng et al. (2016)
for Hermite PCE. Both report favorable results. Roderick et al. (2010) and Li et al. (2011)
apply polynomial regression (PCE) in the context of nuclear engineering. They include deriva-
tive information into the least-squares regression formulation and observe that most polynomial
families are not orthogonal with respect to the *H*¹ inner product. This may deteriorate certain
properties of the regression matrix. To alleviate this issue, Guo et al. (2018) develop a precon-
ditioning procedure for gradient-enhanced sparse regression with certain polynomial families,
with the goal of improving the orthogonality properties of the regression matrix. In all these
approaches, the utilization of derivative information is not straightforward, but requires spe-
cific polynomial families and/or specialized sampling and preconditioning, because the partial
derivatives of a PCE basis do in general not form an orthogonal system. Gejadze et al. (2019)
have derived derivative-enhanced projection methods to compute the PCE coefficients but their
method is restricted to Hermite polynomials and low polynomial degree.

On a different note, the availability of model derivatives has implications also for GSA. The so-
called *Derivative-based Global Sensitivity Measures* (DGSM) are computed by integrating the
squared partial derivatives of the model output over the domain of the inputs. These indices have
been shown to be efficiently estimated by sampling techniques (as Monte Carlo or quasi-Monte
Carlo) as well as from PCE (Sudret and Mai, 2015), and have been proven to be an excellent
screening technique (i.e., detecting all the non-influential inputs among a large number), see
e.g. the review in Kucherenko and Iooss (2017). Indeed, the interpretation of DGSM indices
is straightforward due to their inequality relationship with Sobol' indices: multiplied with the
associated Poincaré constant, DGSM indices provide an upper bound of the total Sobol' index
(Lamboni et al., 2013), regardless of the input probability distribution.

Another way to utilize model derivatives, which solves the issues present for the polynomial
chaos formulation, and naturally provides sharp lower bounds as well as upper bounds on total
Sobol' indices, is to compute *Poincaré chaos expansions* (Roustant et al., 2020a), which we will
abbreviate by *PoinCE* in the sequel. Similar to PCE, PoinCE is a spectral expansion in terms
of an orthonormal basis whose elements are eigenfunctions of the so-called *Poincaré differential operator*. The eigenfunctions are in general non-polynomial, except for the special case of the Gaussian distribution, where they coincide with the Hermite polynomials. The key property of PoinCE is that the partial derivatives of the basis form again an orthogonal basis with respect to the input distribution. This allows to conveniently expand the derivative of the computational model in terms of partial derivatives of the basis (PoinCE-der), which yields another estimator for partial variances and Sobol' indices. If the partial derivatives of the model have smaller variability than the model itself, the estimates based on the model derivatives might be more accurate. This makes PoinCE(-der) an efficient tool for screening (Roustant et al., 2020a).

Our present contribution to the field of generalized chaos expansions and GSA is two-fold. On the theoretical side, we provide a proof that the Poincaré basis is in fact characterized uniquely as the orthonormal basis which remains an orthogonal basis (w.r.t. the same probability measure) after differentiation. Furthermore, we show how PoinCE naturally generalizes an analytical formula for DGSM originally developed for Hermite PCE (Sudret and Mai, 2015), which implies that PoinCE simultaneously and efficiently provides lower and upper bounds to all partial variances. On the computational side, we improve on Roustant et al. (2020a), which introduced projection-based Poincaré chaos and demonstrated that small Sobol' indices were approximated particularly well by the derivative expansion. In this contribution, we compute PoinCE by sparse regression, thus generalizing the powerful and cost-effective sparse PCE methodology to non-polynomial functions. We explore the performance of PoinCE as an estimator for partial variances (upper and lower bounds) and compare it to standard PCE.

This paper is organized as follows. Section 5.2 revisits the mathematical foundations of PoinCE and presents several analytical results related to Sobol' indices and DGSM. Section 5.3 explains the computation of PoinCE basis functions, and the sparse regression methodology adapted from PCE to PoinCE. The methodology is applied in Section 5.4, where two example problems are investigated to demonstrate its performance for sensitivity analysis and screening. Finally, we summarize our conclusions in Section 5.5.

## 5.2 Mathematical background

### 5.2.1 Orthonormal bases in $L^2$

In this section, we recall some important facts about orthonormal bases in $L^2(E, \mu)$ where $E$ is a subset of $\mathbb{R}^d$ and $\mu$ is a probability measure on $E$. We first outline the general theory in Section 5.2.1.1. The particular cases of polynomial and Poincaré bases in several dimensions are developed in Sections 5.2.1.2 and 5.2.1.3.

#### 5.2.1.1 General theory

To begin with, recall that $L^2(E, \mu)$ endowed with the inner product

$$ \langle f, g \rangle = \int_E f(x)g(x)\mu(dx), \quad \text{for } f, g \in L^2(E, \mu) \qquad (5.1) $$
is a Hilbert space. Recall that a sequence of functions $(\Phi_{\alpha})_{\alpha \in \mathcal{I}}$ ($\mathcal{I} \subset \mathbb{N}$) is an orthonormal system in $L^2(E, \mu)$ if it satisfies the two following assumptions:

1) For all $\alpha \neq \alpha'$, $\langle \Phi_{\alpha}, \Phi_{\alpha'} \rangle = 0$, (orthogonality)

2) For all $\alpha$, $\langle \Phi_{\alpha}, \Phi_{\alpha} \rangle = 1$ (unit norm).

An orthonormal system in $L^2(E, \mu)$ is called complete if the closure of the span generated by $(\Phi_{\alpha})$ is $L^2(E, \mu)$. In this case, the system $(\Phi_{\alpha})$ is called an Hilbertian or orthonormal basis of $L^2(E, \mu)$ and for any function $f \in L^2(E, \mu)$ the following expansion holds:

$$f = \sum_{\alpha} \langle \Phi_{\alpha}, f \rangle \Phi_{\alpha}, \quad (\mu \text{ almost surely}). \qquad (5.2)$$

When used to represent random variables in terms of a basis of uncorrelated random variables, such an expansion is often called *chaos expansion* in the uncertainty quantification literature (Wiener, 1938; Ghanem and Spanos, 1991; Ernst et al., 2012).

An archetypal example of chaos expansion is given by the so-called *Fourier expansion*. This corresponds to the case where the set $E = [0, 1]$ is endowed with the Lebesgue measure and we have for $\alpha \in \mathbb{Z}$,

$$\Phi_{\alpha}(x) = \sqrt{2} \cos(2\pi\alpha x) \text{ if } \alpha < 0, \Phi_{0}(x) = 1, \text{ and } \Phi_{\alpha}(x) = \sqrt{2} \sin(2\pi\alpha x) \text{ if } \alpha > 0.$$

In this frame, any square-integrable function $f$ may be expanded as

$$f(x) = a_0 + \sqrt{2} \sum_{\alpha>0} (a_{\alpha} \cos(2\pi\alpha x) + b_{\alpha} \sin(2\pi\alpha x)).$$

Here, for all $\alpha \in \mathbb{Z}_*$,

$$a_0 = \int_0^1 f(x) dx, \quad a_{\alpha} = \sqrt{2} \int_0^1 f(x) \cos(2\pi\alpha) dx, \quad b_{\alpha} = \sqrt{2} \int_0^1 f(x) \sin(2\pi\alpha) dx.$$

When the probability measure is a product measure $\mu = \mu_1 \otimes \dots \otimes \mu_d$ on a product space of intervals $E = E_1 \times E_2 \times \dots \times E_d$, there is a canonical way to build a Hilbertian basis from a collection of univariate Hilbertian ones. Indeed, for $i = 1, \dots, d$ assume that $(\Phi_{\alpha_i}^{(i)})$ is a Hilbertian basis of $L^2(E_i, \mu_i)$. Then, setting $\alpha := (\alpha_1, \dots, \alpha_d)$ and defining the tensor product functions $\Phi_{\alpha} := \prod_{i=1}^d \Phi_{\alpha_i}^{(i)}$, we obtain that $(\Phi_{\alpha})$ is an orthonormal basis of $L^2(E, \mu)$.

In the following we describe two particular chaos types, namely the classical polynomial chaos and the recently developed Poincaré chaos, for a probability measure $\mu$ on $E \subset \mathbb{R}$.

### 5.2.1.2 Polynomial chaos

A classical family of chaos expansions on an interval $E$ of $\mathbb{R}$ endowed with a probability measure are polynomial chaos expansions (PCE) given by orthonormal polynomial bases. A well-known example is the Hermite expansion for which the set $E$ is the whole line $\mathbb{R}$ endowed with the standard Gaussian distribution. In this example, for $\alpha \in \mathbb{N}$, $\Phi_{\alpha} = H_{\alpha}$ is the Hermite polynomial of degree $\alpha$. The first Hermite polynomials are

$$H_0(x) = 1, H_1(x) = x, H_2(x) = \frac{x^2 - 1}{2}, H_3(x) = \frac{x^3 - 3x}{6}, H_4(x) = \frac{x^4 - 6x^2 + 3}{24} \quad (x \in \mathbb{R}).$$
In general, there exists an orthonormal polynomial basis for $L^2(\mu)$ whenever the moment problem for $\mu$ is determinate. This includes the uniform, Gaussian, Beta and Gamma distributions, as well as all distributions with compact support (Ernst et al., 2012).

### 5.2.1.3 Poincaré chaos

The Poincaré basis is another example of an orthonormal basis of $L^2(\mu)$, consisting of functions that admit weak derivatives, i.e. that belong to $H^1(\mu) = \{f \in L^2(\mu) \text{ s.t. } f' \in L^2(\mu)\}$. Recall that $H^1(\mu)$, endowed with the norm $\|f\|_{H^1(\mu)}^2 = \|f\|^2 + \|f'\|^2$, is a Hilbert space. This short summary is based on Roustant et al. (2017) in which more details can be found. We assume that:

**Assumption 2.** The probability measure $\mu$ is supported on a bounded interval $(a, b)$ and admits a density of the form $\rho = e^{-V}$, where $V$ is continuous and piecewise $C^1$ on $[a, b]$ with respect to the Lebesgue measure.

This assumption is sufficient to guarantee the existence of a Poincaré basis. On the topological side, it implies that the Hilbert space $L^2(\mu)$ (resp. $H^1(\mu)$) is equal to $L^2(a, b)$ (resp. $H^1(a, b)$), with an equivalent norm. Indeed, $\mu$ is a bounded perturbation of the uniform measure on $[a, b]$, meaning that the pdf $\rho$ is bounded from below and above by strictly positive constants (by continuity of $V$ on the compact support $[a, b]$).

**Theorem 2 (1D Poincaré basis).** Under Assumption 2, there exists an orthonormal basis $(\varphi_\alpha)_{\alpha \ge 0}$ of $L^2(\mu)$ such that for all $f \in H^1(\mu)$ and for all integer $\alpha \ge 0$, we have:

$$ \langle f', \varphi'_{\alpha} \rangle = \lambda_{\alpha} \langle f, \varphi_{\alpha} \rangle, \qquad (5.3) $$

where $(\lambda_\alpha)_{\alpha \ge 0}$ is an increasing sequence that tends to infinity:

$$ 0 = \lambda_0 < \lambda_1 < \lambda_2 < \dots < \lambda_\alpha \xrightarrow{\alpha \to \infty} +\infty. $$

Here, the inner product $\langle \cdot, \cdot \rangle$ is the one on $L^2(\mu)$ as defined in (5.1). The basis functions $\varphi_\alpha$ are unique up to a sign change, and form the so-called Poincaré basis. Notice that $\varphi_0$ is a constant function equal to $\pm 1$; by convention, we choose $\varphi_0 = 1$.

Furthermore, the Poincaré basis functions are the eigenfunctions of the differential operator

$$ L(f) = f'' - V'f' $$

i.e. satisfy $L(f) = -\lambda f$, subject to Neumann conditions $f'(a) = f'(b) = 0$. The $(\lambda_\alpha)_{\alpha \ge 0}$ are the corresponding eigenvalues.

Finally, $\alpha^{-2}\lambda_\alpha \to \pi^2$ when $\alpha$ tends to infinity, and for all $\alpha \in \mathbb{N}^*$, the eigenfunction $\varphi_\alpha$ has exactly $\alpha$ zeros in $(a, b)$.

*Proof.* The main part of the Theorem can be found in Roustant et al. (2017) or Bakry et al. (2014). The two last assertions come by rewriting the differential equation $f'' - V'f' = -\lambda f$ in the Sturm-Liouville form

$$ -(pf')' + qf = \lambdawf \qquad (5.4) $$
with $p = w = \rho$ and $q = 0$. Then, by the Sturm-Liouville theory (see e.g. Zettl (2010, Theorem 4.3.1, (1), (6) and (7))), we have that $\alpha^{-2}\lambda_{\alpha} \to \pi^2$ when $\alpha$ tends to infinity, and for all $\alpha \in \mathbb{N}^*$, the eigenfunction $\varphi_{\alpha}$ has exactly $\alpha$ zeros in $(a, b)$. $\square$

The Poincaré basis shares some similarity with both the polynomial chaos and the Fourier basis in terms of oscillations: by Theorem 2, the higher the order of the eigenvalue, the more oscillating the corresponding eigenfunction.

For some specific cases, the Poincaré basis is known analytically. For instance, for the uniform distribution, the Poincaré basis is a kind of Fourier basis (see e.g. Roustant et al. (2020a, §4)). Otherwise it has to be computed numerically, e.g., by a finite element technique (see Section 5.3.1).

Note that Assumption 2 is a convenient sufficient condition which guarantees the existence of a Poincaré basis. It is satisfied for a large range of truncated parametric probability distributions. The set of probability distributions for which the Poincaré chaos exists is larger, but not well known. For instance, the Poincaré chaos is defined for the Gaussian distribution, and then coincides with polynomial chaos, corresponding to Hermite polynomials. This is the only case where Poincaré chaos and polynomial chaos coincide (Bakry et al., 2014, §2.7). On the other hand, Poincaré chaos is not defined for the Laplace distribution, since the eigenvalues of the associated operator do not form a countable set (Bakry et al., 2014, §4.4.1).

The Poincaré basis is useful for sensitivity analysis. First, it is linked to the Poincaré inequality

$$ \mathrm{Var}_{\mu}(f) \leq C_P(\mu) \int f'^2 d\mu, \qquad (5.5) $$

which holds for all functions $f \in H^1(\mu)$ under the assumptions on $\mu$. Indeed, the smallest constant $C_P(\mu)$ such that (5.5) is satisfied is equal to $C_P(\mu) = 1/\lambda_1$, and choosing $f = \varphi_1$ corresponds to the equality case (Roustant et al., 2017). Roughly speaking, the Poincaré basis function associated to the first non-zero eigenvalue is the function with the largest possible variance for a given amount of integrated squared derivative (in the $L^2$ sense). A second appealing property for the analysis of variance is that the derivatives of the Poincaré basis remain orthogonal functions:

**Proposition 1.** Under Assumption 2, the sequence $(\frac{1}{\sqrt{\lambda_{\alpha}}} \varphi'_{\alpha})_{\alpha \geq 1}$ is an orthonormal basis of $L^2(\mu)$.

*Proof.* The orthonormality of the sequence is a consequence of (5.3) by choosing $f = \varphi_{\beta}$, with $\beta \in \mathbb{N}^*$. It remains to show that the system is dense in $L^2(\mu)$, or equivalently, that its orthogonal is null. Let thus $f \in L^2(\mu)$ such that

$$ \langle f, \varphi'_{\alpha} \rangle = 0, \quad \text{for all } \alpha \geq 1. $$

As explained when stating Assumption 2, $L^2(\mu)$ (resp. $H^1(\mu)$) is equal to $L^2(a, b)$ (resp. $H^1(a, b)$), with an equivalent norm. Now, there exists $g \in H^1(\mu)$ such that $f = g'$. Indeed, let us define $g$ by $g(x) = g(a) + \int_a^x f(t)dt$. As $f \in L^2(\mu) = L^2(a, b)$, then $g$ belongs to $H^1(a, b) = H^1(\mu)$, and $g' = f$. Then we have

$$ \langle g', \varphi'_{\alpha} \rangle = 0, \quad \text{for all } \alpha \geq 1. $$
By (5.3), we obtain $\langle g, \varphi_\alpha \rangle = 0$ for all $\alpha \ge 1$ (as $\lambda_\alpha > 0$ for $\alpha \ge 1$). As the functions $\varphi_\alpha$ form an orthonormal basis of $L^2(\mu)$ with $\varphi_0 = 1$, this implies that $g$ is a constant function, and finally $f = 0$. The proof is completed. $\square$

In fact, the property in Proposition 1, i.e., that the derivatives of the Poincaré basis form again an orthogonal basis in $L^2(\mu)$, uniquely characterizes the Poincaré basis:

**Proposition 2.** Under Assumption 2, Poincaré bases are the only orthonormal bases $(\varphi_\alpha)$ of $L^2(\mu)$ in $H^1(\mu)$ such that $(\varphi'_\alpha)$ is an orthogonal basis of $L^2(\mu)$.

This result seems difficult to find in the literature. German-speaking readers can find a similar proposition in Mikolas (1955), stated in the frame of Sturm-Liouville theory for twice-differentiable functions satisfying boundary conditions. See also Kwon and Lee (2003) for a similar result under the assumption that all functions involved in the Sturm-Liouville problem Eq. (5.4) are of class $C^\infty$. We provide below a proof based on Hilbertian arguments.

As a corollary, if there exists a basis different from the Poincaré basis for which derivatives form an orthogonal system, then that system is not dense in $L^2(\mu)$. As an example, for the uniform probability measure on $[0, 2\pi]$, consider the usual Fourier basis formed by $\{\cos(nx), \sin(nx) : n \ge 0\}$ (up to multiplicative constants). Taking derivatives results in the same set of functions (up to multiplicative constants) – except for the constant function $\cos(0x) = 1$. Thus, the derivatives form an orthogonal system which covers the orthogonal of constant functions in $L^2(\mu)$, which is a strict subspace of $L^2(\mu)$. Meanwhile, the Poincaré basis for this probability measure is formed by functions proportional to $\cos(\frac{n}{2}x)$ for $n \ge 0$. Proposition 2 guarantees that all functions of $L^2(\mu)$, including the constant functions, are spanned by the derivatives. Indeed, this is explained intuitively by the presence of half-frequencies: when $n$ is odd, the functions $\sin(\frac{n}{2}x)$ are not orthogonal to 1.

*Proof of Proposition 2.* The fact that a Poincaré basis remains an orthogonal basis by derivation has been proved in Proposition 1. Conversely, let $(\varphi_\alpha)_{\alpha \ge 0}$ be a system of $H^1(\mu)$, with $\varphi_0 = 1$, such that $(\varphi_\alpha)$ is an orthonormal basis of $L^2(\mu)$ and $(\varphi'_\alpha)_{\alpha \ge 1}$ is an orthogonal basis of $L^2(\mu)$. Let us first prove that $(\varphi_\alpha)$ is an orthogonal basis of $H^1(\mu)$. The orthogonality is a direct consequence of the definition of the inner product of $H^1(\mu)$:

$$ \langle \varphi_\alpha, \varphi_\beta \rangle_{H^1(\mu)} = \langle \varphi_\alpha, \varphi_\beta \rangle_{L^2(\mu)} + \langle \varphi'_\alpha, \varphi'_\beta \rangle_{L^2(\mu)} = (1 + \|\varphi'_\alpha\|_{L^2(\mu)}^2) \delta_{\alpha,\beta}. $$

Let us prove that $(\varphi_\alpha)$ is dense in $H^1(\mu)$. As explained when stating Assumption 2, $L^2(\mu)$ (resp. $H^1(\mu)$) is equal to $L^2(a, b)$ (resp. $H^1(a, b)$), with an equivalent norm. Hence, it is equivalent to prove that $(\varphi_\alpha)$ is dense in $H^1(a, b)$. Now, let $f$ be in $H^1(a, b)$. As $(\varphi'_\alpha)$ is dense in $L^2(a, b)$ (equivalently in $L^2(\mu)$), then $f'$ expands as $f' = \sum_{\alpha \in \mathbb{N}} c_\alpha \varphi'_\alpha$. In $H^1(a, b)$ each function is equal to the primitive function of its derivative, hence we have:

$$ \begin{aligned} \left| f(t) - f(a) - \sum_{\alpha=1}^{N} c_{\alpha}(\varphi_{\alpha}(t) - \varphi_{\alpha}(a)) \right| &= \left| \int_{a}^{t} \left( f'(x) - \sum_{\alpha=1}^{N} c_{\alpha} \varphi'_{\alpha}(x) \right) dx \right| \\ &\le (b-a) \left\| f' - \sum_{\alpha=1}^{N} c_{\alpha} \varphi'_{\alpha} \right\|_{L^2(a,b)} \end{aligned} $$
where the inequality comes from the Cauchy-Schwarz inequality. We deduce that $\|f - f(a) - \sum_{\alpha=1}^{N} c_{\alpha}(\varphi_{\alpha} - \varphi_{\alpha}(a))\|_{L^2(a,b)} \to 0$ when $N$ tends to infinity. Together with $f' = \sum_{\alpha \in \mathbb{N}} c_{\alpha} \phi'_{\alpha}$, this implies that $\|f - f(a) - \sum_{\alpha=1}^{N} c_{\alpha}(\varphi_{\alpha} - \varphi_{\alpha}(a))\|_{H^1(a,b)} \to 0$. As $\varphi_0 = 1$, this proves that $(\varphi_{\alpha})$ is dense in $H^1(a, b)$, which was to be proved.

Now, let us fix $\alpha \ge 0$. Consider the linear form $L_{\alpha}$ defined on $H^1(\mu)$ by $L_{\alpha}(f) = \langle f', \varphi'_{\alpha} \rangle_{L^2(\mu)}$. The Cauchy-Schwarz inequality gives $|L_{\alpha}(f)| \le \|f'\|_{L^2(\mu)} \|\varphi'_{\alpha}\|_{L^2(\mu)} \le \|f\|_{H^1(\mu)} \|\varphi'_{\alpha}\|_{L^2(\mu)}$. This proves that $L_{\alpha}$ is continuous. Hence, by the Riesz representation theorem, there exists a unique $\zeta_{\alpha} \in H^1(\mu)$ such that for all $f \in H^1(\mu)$, $L_{\alpha}(f) = \langle f, \zeta_{\alpha} \rangle_{H^1(\mu)}$, i.e. $\langle f', \varphi'_{\alpha} \rangle_{L^2(\mu)} = \langle f, \zeta_{\alpha} \rangle_{H^1(\mu)}$. Choosing $f = \varphi_{\beta}$ with $\beta \ne \alpha$, we obtain by orthogonality of $(\varphi'_{\alpha})$ that for all $\beta \ne \alpha$, $\langle \varphi_{\beta}, \zeta_{\alpha} \rangle_{H^1(\mu)} = 0$. As $(\varphi_{\beta})_{\beta \ge 0}$ is an orthogonal basis of $H^1(\mu)$, this implies that $\zeta_{\alpha}$ is collinear to $\varphi_{\alpha}$, i.e., there exists $\tilde{\lambda}_{\alpha} \in \mathbb{R}$ such that $\zeta_{\alpha} = \tilde{\lambda}_{\alpha} \varphi_{\alpha}$. Thus, for all $f \in H^1(\mu)$, we have $\langle f', \varphi'_{\alpha} \rangle_{L^2(\mu)} = \tilde{\lambda}_{\alpha} \langle f, \varphi_{\alpha} \rangle_{H^1(\mu)}$. Choosing $f = \varphi_{\alpha}$, we get $\tilde{\lambda}_{\alpha} = \frac{\|\varphi'_{\alpha}\|_{L^2(\mu)}^2}{1+\|\varphi'_{\alpha}\|_{L^2(\mu)}^2}$, which belongs to $[0, 1)$. Finally, we obtain that $\langle f', \varphi'_{\alpha} \rangle_{L^2(\mu)} = \lambda_{\alpha} \langle f, \varphi_{\alpha} \rangle_{L^2(\mu)}$, where $\lambda_{\alpha} = \frac{\tilde{\lambda}_{\alpha}}{1-\tilde{\lambda}_{\alpha}}$ is a non-negative real number. As it is true for all $f$ in $H^1(\mu)$ and all $\alpha \in \mathbb{N}$, this implies, by uniqueness of the Poincaré basis (under Assumption 2), that $(\varphi_{\alpha})_{\alpha \ge 0}$ is a Poincaré basis. $\square$

Turning to higher dimensions, we assume that for all $i = 1, \dots, d$, the probability measure $\mu_i$ satisfies Assumption 2, and we denote by $(\varphi_{i,\alpha_i})_{\alpha_i \ge 0}$ the sequence of 1-dimensional Poincaré basis functions, and by $(\lambda_{i,\alpha_i})_{\alpha_i \ge 0}$ the sequence of associated eigenvalues. The Poincaré chaos basis is then defined by the tensor product $\Phi_{\alpha} = \varphi_{1,\alpha_1} \otimes \dots \otimes \varphi_{d,\alpha_d}$. Using the properties of $L^2$ bases, (5.3) thus implies that for all $f \in H^1(\mu)$, for all $i = 1, \dots, d$:

$$ \left\langle \frac{\partial f}{\partial x_i}, \frac{\partial \Phi_{\alpha}}{\partial x_i} \right\rangle = \lambda_{i, \alpha_i} \langle f, \Phi_{\alpha} \rangle. \qquad (5.6) $$

Similarly, applying Proposition 1, we get:

**Proposition 3.** Under Assumption 2, for all $i = 1, \dots, d$, the sequence $\left( \frac{1}{\sqrt{\lambda_{i,\alpha_i}}} \frac{\partial \Phi_{\alpha}}{\partial x_i} \right)_{\alpha, \alpha_i \ge 1}$ is an orthonormal basis of $L^2(\mu)$.

### 5.2.2 Variance-based indices, derivative-based indices

We first recall the definition of variance-based sensitivity indices, which quantify the importance of each input variable in terms of function response variability.

Let $f$ be a real-valued function defined on $E = E_1 \times \dots \times E_d \subseteq \mathbb{R}^d$. The uncertainty of the inputs is represented by a random vector $\mathbf{X} = (X_1, \dots, X_d)^T$ with probability measure $\mu$ on $E$. We further assume that the $X_i$'s are independent and that $f(\mathbf{X})$ belongs to $L^2(E, \mu)$. Denoting by $\mu_i$ the marginal distribution of $X_i$ on $E_i$ ($i = 1, \dots, d$), we then have $\mu = \mu_1 \otimes \dots \otimes \mu_d$. In this framework, $f(\mathbf{X})$ can be decomposed uniquely as a sum of terms of increasing complexity

$$ f(\mathbf{X}) = f_0 + \sum_{1 \le i \le d} f_i(X_i) + \sum_{1 \le i < j \le d} f_{i,j}(X_i, X_j) + \dots + f_{1, \dots, d}(X_1, \dots, X_d) \quad (5.7) $$

under centering conditions $\mathbb{E}[f_I(X_I)] = 0$ and non-overlapping conditions $\mathbb{E}[f_I(X_I)|X_J] = 0$, for all sets $I \subseteq \{1, \dots, d\}$ and all strict subsets $J$ of $I$. We have used the set notation $X_I$ to represent
the subvector of **X** obtained by selecting the coordinates belonging to **I**. These conditions imply that all the terms of (5.7) are orthogonal, leading to the variance decomposition

$$ \mathrm{Var}f(\mathbf{X}) = \sum_{1 \le i \le d} \mathrm{Var}f_i(X_i) + \sum_{1 \le i < j \le d} \mathrm{Var}f_{i,j}(X_i, X_j) + \cdots + \mathrm{Var}f_{1, \dots, d}(X_1, \dots, X_d) \quad (5.8) $$

Due to this property, the functional decomposition (5.7) is often called ANOVA (ANalysis Of VAriance) decomposition. Originating from Hoeffding (1948), it was revisited by Efron and Stein (1981), Antoniadis (1984), and Sobol' (1993). For a given set $I \subseteq \{1, \dots, d\}$, we call the corresponding term of (5.8) partial variance (denoted $D_I$), and call its normalized version Sobol' index (denoted $S_I$):

$$ D_I = \mathrm{Var}(f_I(X_I)), \quad S_I = \frac{D_I}{D}, $$

where $D = \mathrm{Var}f(\mathbf{X})$ is the overall variance (total variance). In particular, for $i \in \{1, \dots, d\}$, the first-order Sobol' index $S_i$ corresponds to the proportion of variance of $f(\mathbf{X})$ explained by $X_i$ only. In order to include also the interactions of $X_i$ with the other variables, the total partial variance and the total Sobol' index are defined by

$$ D_i^{\mathrm{tot}} = \sum_{I \supseteq \{i\}} \mathrm{Var}(f_I(X_I)), \quad S_i^{\mathrm{tot}} = \frac{D_i^{\mathrm{tot}}}{D}. $$

Note that practitioners also call the (total) partial variances unnormalized (total) Sobol' indices. In the sequel, we will use these two words interchangeably.

The total Sobol' index can be used for screening. Indeed, under mild conditions, if $S_i^{\mathrm{tot}} = 0$ then the function $f$ does not depend on $x_i$ over $E$ (in the pointwise sense).

When the derivatives are available, a global sensitivity index can be obtained by integration. The so-called *derivative-based sensitivity measure* (DGSM) index of $f$ with respect to $X_i$ (Sobol and Gresham, 1995; Kucherenko et al., 2009) is defined by

$$ \nu_i = \mathbb{E} \left[ \left( \frac{\partial f}{\partial x_i}(\mathbf{X}) \right)^2 \right] = \int_{\mathbb{R}^d} \left( \frac{\partial f}{\partial x_i}(\mathbf{x}) \right)^2 d\mu(\mathbf{x}) = \left\| \frac{\partial f}{\partial x_i} \right\|^2. \quad (5.9) $$

Contrarily to variance-based indices, DGSM are not associated to a variance decomposition. Nevertheless, they can be used for screening. Indeed, under mild conditions, $\nu_i = 0$ implies that $f$ does not depend on $x_i$ over $E$.

### 5.2.3 Chaos expansion serving sensitivity analysis

One main advantage of using an orthonormal basis for sensitivity analysis is that, once the expansion has been obtained, the variance-based indices can be computed in a straightforward way as a sum of squared coefficients (Sudret, 2006, 2008). More precisely, let $f$ be in $L^2(\mu)$, and let $(\Phi_\alpha)_{\alpha \in \mathbb{N}^d}$ be a multivariate orthonormal basis obtained by tensorization as described in Section 5.2.1. The expansion of $f$ in this basis is given by

$$ f = \sum_{\alpha \in \mathbb{N}^d} c_\alpha \Phi_\alpha. \quad (5.10) $$
By using the orthonormality we obtain the expression of the total variance

$$D = \sum_{\alpha \neq 0} c_{\alpha}^{2}. \qquad (5.11)$$

The expression of the total Sobol' index $S_i^{\text{tot}}$ is obtained by only considering the terms of the decomposition (5.10) that contain the variable $x_i$, i.e. such that $\alpha_i \ge 1$. Hence, we have $S_i^{\text{tot}} = \frac{D_i^{\text{tot}}}{D}$ with $D_i^{\text{tot}}$ the total partial variance

$$D_i^{\text{tot}} = \sum_{\alpha, \alpha_i \ge 1} c_{\alpha}^2. \qquad (5.12)$$

The first-order Sobol' index $S_i^1$ relies on the terms that include $x_i$ only, i.e., $S_i^1 = \frac{D_i^1}{D}$ with

$$D_i^1 = \sum_{\substack{\alpha, \alpha_i \ge 1, \\ \alpha_j=0 \text{ for } j \neq i}} c_{\alpha}^2. \qquad (5.13)$$

Let us now consider the case where the gradient of $f$ is available. The Poincaré basis is particularly suited to this situation. Indeed, we can derive in a straightforward way expressions of both variance-based and derivative-based indices, involving the derivatives of $f$. Due to orthonormality, the coefficients of the basis expansion in (5.10) are given by the projection of $f$ onto the associated basis element:

$$c_{\alpha} = \langle f, \Phi_{\alpha} \rangle. \qquad (5.14)$$

From now on, let $(\Phi_{\alpha})_{\alpha}$ denote the Poincaré basis. Combining (5.6) and (5.14), and assuming that $\alpha_1 \ge 1$, $c_{\alpha}$ can be written using the partial derivatives w.r.t variable $X_1$ (Roustant et al., 2020a):

$$c_{\alpha} = \langle f, \Phi_{\alpha} \rangle = \frac{1}{\lambda_{1,\alpha_1}} \left\langle \frac{\partial f}{\partial x_1}, \frac{\partial \Phi_{\alpha}}{\partial x_1} \right\rangle = \frac{1}{\lambda_{1,\alpha_1}} \left\langle \frac{\partial f}{\partial x_1}, \frac{\partial \varphi_{1,\alpha_1}}{\partial x_1} \otimes \varphi_{2,\alpha_2} \otimes \cdots \otimes \varphi_{d,\alpha_d} \right\rangle \quad (5.15)$$

and equivalently using partial derivatives w.r.t variable $X_i$ if $\alpha_i \ge 1$. Thus, (5.11), (5.12) and (5.13) can also be computed using the various partial derivatives of $f$. Whereas the theoretical expressions are equal, their estimators have different properties. For example, if the integral is evaluated by Monte Carlo simulation, the expression whose integrand has smaller variance will be more accurate. We describe in Section 5.3 the computation of the expansion coefficients by regression, and we empirically compare the two estimation procedures in Section 5.4.

Furthermore, DGSM can be computed directly from the Poincaré expansion. More precisely, we have the following proposition.

**Proposition 4 (DGSM formula for Poincaré chaos).** Let $f \in H^1(\mu)$. Let $f = \sum_{\alpha} c_{\alpha} \Phi_{\alpha}$ be the expansion of $f$ in the Poincaré chaos basis, with $c_{\alpha} = \langle f, \Phi_{\alpha} \rangle$. Then the DGSM index of $f$ with respect to $X_i$ is equal to:

$$\nu_i = \sum_{\alpha, \alpha_i \ge 1} \lambda_{i, \alpha_i} (c_{\alpha})^2. \qquad (5.16)$$

*Proof.* Write $f = \sum_{\alpha} c_{\alpha} \Phi_{\alpha}$. Then by Proposition 3, we get

$$\frac{\partial f}{\partial x_i} = \sum_{\alpha, \alpha_i \ge 1} c_{\alpha} \frac{\partial \Phi_{\alpha}}{\partial x_i},$$
where we can constrain the sum to multi-indices $\alpha$ such that $\alpha_i \ge 1$, since $\varphi_{i,\alpha_i} = 1$ for $\alpha_i = 0$. Now, using again the orthogonality of Poincaré basis derivatives (Proposition 3), it follows that

$$ \nu_i = \left\| \frac{\partial f}{\partial x_i} \right\|^2 = \sum_{\alpha, \alpha_i \ge 1} (c_\alpha)^2 \left\| \frac{\partial \Phi_\alpha}{\partial x_i} \right\|^2 = \sum_{\alpha, \alpha_i \ge 1} \lambda_{i, \alpha_i} (c_\alpha)^2 . $$

□

Formula (5.16) extends a previous result given by Sudret and Mai (2015) when all the $\mu_i$ are standard Gaussian. Indeed, in that case, Poincaré chaos coincides with polynomial chaos, and $\lambda_{i,\alpha_i} = \alpha_i$.

Using the expressions provided in (5.12) and (5.16) and an inequality derived by Sobol' and Kucherenko (2009) and Lamboni et al. (2013), we obtain lower and upper bounds to total partial variances as follows:

$$ \sum_{\alpha \in \mathcal{A}, \alpha_i \ge 1} (c_\alpha)^2 \le D_i^{\text{tot}} \le C_P(\mu_i) \nu_i = \sum_{\alpha \in \mathbb{N}^d, \alpha_i \ge 1} \frac{\lambda_{i, \alpha_i}}{\lambda_{i, 1}} (c_\alpha)^2 . \quad (5.17) $$

The lower bound is an obvious consequence of the truncation. The upper bound holds only for the full infinite expansion and is otherwise underestimated. Comparing the form of the right-hand side of Eq. (5.17) with the total Sobol' formula Eq. (5.12) gives insight into how tight this upper bound is: equality is attained only if the Poincaré chaos expansion does not contain terms of higher degree than 1 for $X_i$ (then, $\frac{\lambda_{i,\alpha_i}}{\lambda_{i,1}} = 1$). Else, depending on the decay behavior of $c_\alpha$ the gap can be significant, since the eigenvalues are diverging to infinity (see Theorem 2).

## 5.3 Computation of sparse Poincaré expansions

Let $f \in H^1(\mu, E)$ be a computational model defined on the input space $E \subset \mathbb{R}^d$, with independent input random variables and with the input probability measure $\mu$ admitting a probability density function $\rho$ fulfilling Assumption 2 for each of the marginals. In the remainder of this paper, we assume that $\rho$ is known. We also assume that we are provided with an i.i.d. sample from the input distribution and with the corresponding model evaluations and model gradient values at each of the points.

With Poincaré expansion (PoinCE) we denote the expansion of the computational model onto the Poincaré basis

$$ f(\boldsymbol{x}) = \sum_{\alpha} c_{\alpha} \Phi_{\alpha}(\boldsymbol{x}), \quad (5.18) $$

and with Poincaré derivative expansion in direction $i$ (PoinCE-der-$i$) the expression

$$ \frac{\partial f}{\partial x_i}(\boldsymbol{x}) = \sum_{\alpha, \alpha_i \ge 1} c_{\alpha}^{\partial, i} \frac{\partial \Phi_{\alpha}}{\partial x_i}(\boldsymbol{x}), \quad (5.19) $$

or the equivalent expansion using normalized basis derivatives that have unit norm in $L^2(\mu)$. Note that (5.19) is the partial derivative of (5.18) w.r.t. variable $X_i$. Because the zeroth order basis function of a Poincaré basis is the constant function, basis terms for which $\alpha_i = 0$ have
zero partial derivative w.r.t. $X_i$ and are not included in (5.19). While in theory by Equation (5.15), the two expressions (5.18) and (5.19) provide identical coefficients for corresponding basis elements, i.e., $c_\alpha = c_\alpha^{\partial, i}$ for $\alpha \in \{\alpha' \in \mathcal{A} : \alpha'_i \ge 1\}$, in practice they will not coincide when estimated from a data set of finite size. This will be investigated in Section 5.4 for a number of numerical examples.

In this section, we describe how such expansions are computed in practice: this concerns the computation of the Poincaré basis functions, the choice of truncation, the location of the sampled points, and the method for computing the coefficients. The implementation relies on and integrates into the UQLab framework (Marelli and Sudret, 2014).

### 5.3.1 Implementation of Poincaré basis functions

As described in Section 5.2.1.3, Poincaré basis functions are tensor products of univariate Poincaré basis functions. Each 1D basis consists of the eigenfunctions of the Poincaré differential operator associated with the respective marginal distribution (Theorem 2).

A Poincaré basis is guaranteed to exist for marginal distributions fulfilling Assumption 2 and for the Gaussian distribution. Other distributions have to be transformed or truncated to allow for a Poincaré basis. Since an isoprobabilistic transformation to standard variables can be highly nonlinear (Torre et al., 2019; Oladyshkin and Nowak, 2012), we opt for truncation: if the distribution is not Gaussian and has (one- or two-sided) unbounded support, we truncate it to its $10^{-6}$- and $(1 - 10^{-6})$-quantiles, respectively.<fnref n="1" />

We consider here only standard parametric families of probability densities (bounded and unbounded), although a Poincaré basis can be computed for any input distribution which after truncation fulfills Assumption 2. In particular, without any changes to the methodology PoinCE could be used in a data-driven framework (Torre et al., 2019) by computing the Poincaré basis for a dimensionwise kernel density estimate of the input distribution (assuming independence) given the available data.

As can be seen from applying the change-of-variables formula for a linear transformation to (5.3), the eigenvalues of the Poincaré differential operator scale with the inverse of the squared support interval length. To avoid numerical difficulties, we therefore linearly transform (i.e., shift and rescale) parametric families to standard parameters using

* their bounds in the case of uniform, beta, triangular;

* their location and scale parameter in the case of Gaussian, Gumbel, Gumbel-min, Laplace, logistic;

* their (inverse) scale parameter in the case of exponential, gamma, Weibull, lognormal.

<footnote n="1">¹One might argue that this can distort the results obtained with PoinCE, especially in the tails. It is true that this truncation introduces a small error. However, as all such methods, PoinCE by design approximates accurately mainly the bulk, not the tails (for this, specialized techniques like subset simulation shall be used). Furthermore, in practical applications it is a modelling choice how to represent the input distribution. Choosing an unbounded parametric distribution is common, but not necessarily the most sensible choice, since for virtually every quantity in the real world there is an upper bound that cannot be exceeded.</footnote>
In the current implementation, distributions not belonging to this group of families are not being rescaled.

For standard uniform ($\mathcal{U}([-0.5, 0.5])$) and standard Gaussian ($\mathcal{N}(0, 1)$) marginals, the Poincaré basis can be analytically computed and is given by the Fourier (cosine) basis and the Hermite polynomial basis, respectively (Roustant et al., 2020a). Therefore, in the special case of uniform or Gaussian marginals, we always (after rescaling) use the analytical solution.

For all other marginals, the Poincaré basis is computed numerically using linear finite elements. We use a fine uniform grid within the bounds and piecewise linear functions with local support, commonly called ‘hat’ functions. Using the weak formulation of the eigenvalue problem of the Poincaré differential operator given in (5.3), we arrive at the shifted generalized eigenvalue problem

$$ K\mathbf{a}^{(n)} = (\lambda_n + 1)M\mathbf{a}^{(n)} \quad (5.20) $$

as described in Roustant et al. (2017, section 4.3), where the eigenvector $\mathbf{a}^{(n)}$ denotes the vector of coefficients used to express eigenfunction $\varphi_n$ in terms of ‘hat’ functions. Here $M$ is the mass matrix, and $K$ is the sum of mass- and stiffness matrix. After solving this problem using Matlab’s builtin function `eigs`, we interpolate the discrete eigenvectors with piecewise cubic splines, prescribing zero derivatives at the interval boundaries. Then, the basis derivatives are computed using centered finite differences. While more sophisticated techniques (e.g., Hermitian $C^1$ elements, or Haar wavelets (Bujurke et al., 2008)) could of course be used to improve this numerical computation procedure, it is accurate enough for our purposes of demonstrating the usefulness of Poincaré. Eigenfunctions and eigenfunction derivatives are scaled to have unit norm with respect to the measure $\mu_i$.

### 5.3.2 Choice of the basis truncation

In practice, the series in (5.18) and (5.19) cannot include an infinite number of terms, but must be truncated to a finite expansion. We denote by $\mathcal{A} \subset \mathbb{N}^d$ the subset of multi-indices that are included in the expansion. For PCE, $\mathcal{A}$ is typically chosen to include terms up to a certain degree $p$, resulting in the so-called *total degree basis*

$$ \mathcal{A}^p = \{\alpha \in \mathbb{N}^d : \sum_{i=1}^{d} |\alpha_i| \le p\}. \quad (5.21) $$

To further restrict the number of terms used in the expansion, another common truncation method is *hyperbolic truncation* (Blatman and Sudret, 2011)

$$ \mathcal{A}^{p,q} = \{\alpha \in \mathbb{N}^d : \|\alpha\|_q \le p\}. \quad (5.22) $$

with the $\ell^q$-(*quasi-*)norm $\|\alpha\|_q = (\sum_{i=1}^d \alpha_i^q)^{\frac{1}{q}}$ for $q \in (0, 1]$.

Since the Poincaré basis is in general not polynomial, the concept of polynomial degree cannot be used to characterize the basis functions. Instead, we use the natural order of the basis functions corresponding to the increasing sequence of Poincaré eigenvalues, which also corresponds to an increasing number of oscillations (Theorem 2; recall that the $n$th eigenfunction has $n$ zeros).
Therefore, we use the PCE terminology “degree” also for PoinCE. In particular, a degree of $\alpha_i = 0$ denotes the constant basis function $\varphi_{i,0}(x_i) = 1$ associated to the eigenvalue $\lambda_{i,0} = 0$.

Often, in practice it is not known which degree is needed for a given problem. While in theory the expansion is more accurate the larger the total degree is, in practice accuracy is limited by the number of available sample points, since the quality of the regression solution (see Section 5.3.3) depends on the ratio of sample points to basis elements. In that case, a successful strategy consists of applying degree adaptivity, i.e., choosing the best degree for the expansion by cross-validation (Blatman and Sudret, 2011; Lüthen et al., 2022a). This procedure is computationally inexpensive, since it only requires a new surrogate model fit for each new total degree, but no additional model evaluations. We apply leave-one-out (LOO) cross-validation together with a modification factor introduced by Chapelle et al. (2002) (Blatman and Sudret, 2011).

Both hyperbolic truncation and degree adaptivity contribute to the sparsity of the resulting expansion by identifying a suitable subset of basis functions necessary for a good approximation. Sparsity is a successful concept in regression-based PCE (Lüthen et al., 2021). Denote by $P = |\mathcal{A}|$ the number of basis elements in the truncated expansion. $\mathcal{A}$, also called candidate basis, contains the basis elements available for approximation. We describe below how sparse regression further selects only a subset of $\mathcal{A}$ to be active, i.e., have a nonzero coefficient. The final expansion might (and indeed often will) have less than $P$ active terms.

### 5.3.3 Computation of the coefficients by sparse regression

For computing the coefficients of an orthogonal expansion as in (5.18) and (5.19), there exist two main approaches. One is projection: the model $f$ is projected onto the basis functions, see (5.2). The resulting integral may be evaluated by Monte Carlo (MC) simulation, as done by Roustant et al. (2020a) for Poincaré chaos, or by (sparse) quadrature methods (Le Maître et al., 2002; Matthies and Keese, 2005; Constantine et al., 2012). However, note that in general MC converges slowly, while quadrature (even when sparse) is affected by the curse of dimensionality.

The second approach is regression, introduced for PCE by Blatman and Sudret (2008). Here, after choosing an experimental design $\mathcal{X} = \{\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(N)}\}$ of input points, (5.18) is discretized as

$$ \mathbf{y} \approx \mathbf{\Psi} \mathbf{c} \qquad (5.23) $$

where $\mathbf{y} = (f(\mathbf{x}^{(1)}), \dots, f(\mathbf{x}^{(N)}))^T$ is the vector of model evaluations, $\mathbf{\Psi} \in \mathbb{R}^{N \times P}$ is the regression matrix with entries $\mathbf{\Psi}_{kj} = \Phi_j(\mathbf{x}^{(k)})$ where $j$ refers to an enumeration of the multivariate basis $(\Phi_\alpha)_{\alpha \in \mathcal{A}}$, and $\mathbf{c}$ is the vector of expansion coefficients. The discretization of (5.19) is analogous, with a vector

$$ \mathbf{y}_{\partial,i} = \left( \frac{\partial f}{\partial x_i}(\mathbf{x}^{(1)}), \dots, \frac{\partial f}{\partial x_i}(\mathbf{x}^{(N)}) \right)^T \qquad (5.24) $$

containing model partial derivatives and a regression matrix $\mathbf{\Psi}_{\partial,i}$ with entries

$$ \mathbf{\Psi}_{kj}^{\partial,i} = \frac{1}{\sqrt{\lambda_{i,\alpha_j(i)}}} \frac{\partial \Phi_j}{\partial x_i}(\mathbf{x}^{(k)}), \qquad (5.25) $$

where $\alpha_j(i)$ denotes the $i$th component of the $j$th basis element characterized by the multi-index $\alpha_j$ (see also Proposition 3).
The regression problem can be solved by ordinary least squares as

$$ \hat{c} = \arg \min_{c} \| \Psi c - y \|_{2}^{2}, \qquad (5.26) $$

provided that enough model evaluations are available – at least $N \ge P$, or better $N \ge kP$ with $k = 2,3$ to avoid overfitting. Due to the rapid growth of the total-degree basis with increasing dimension and degree, this requirement on model evaluations is often too restrictive for real-world problems.

To avoid this problem, sparse regression can be used, which regularizes the problem by encouraging solutions with few nonzero coefficients (Candès and Wakin, 2008; Kougioumtzoglou et al., 2020). An example is $\ell^1$-minimization:

$$ \hat{c} = \arg \min_{c} \| \Psi c - y \|_{2}^{2} + \lambda \| c \|_{1}. \qquad (5.27) $$

The $\ell^1$-norm penalizes the coefficient vector so that sparse solutions are preferred. The sparse regression formulation allows for accurate solutions even in the case $N < P$. There exist many sparse regression methods utilizing different formulations of the sparse regression problem, see e.g. Lüthen et al. (2021) for an overview of available sparse regression solvers in the context of PCE. In this work, we use the sparse solver Hybrid Least Angle Regression (Hybrid-LARS) (Blatman and Sudret, 2011; Marelli et al., 2021a).

A result by Candès and Plan (2011) on sparse recovery emphasizes the importance of isotropy of the row distribution of the regression matrix, i.e., the requirement that for a row $a = (\Phi_{\alpha_1}(x), \dots, \Phi_{\alpha_P}(x))$ of the regression matrix $\Psi$ it holds that $E[a^T a] = I_P$, where $I_P$ is the identity matrix of size $P$, and the expectation is with respect to the distribution of the experimental design points. If the experimental design points are chosen to follow the input distribution, the distributions of regression matrix rows for Poincaré as well as for normalized Poincaré derivative expansions are isotropic by construction due to orthonormality of the bases w.r.t. the input distribution. To improve the space-filling property of the experimental design, we use Latin Hypercube Sampling (LHS) (McKay et al., 1979) with maximin distance optimization.

### 5.3.4 Coefficients and Sobol' indices for Poincaré derivative expansions

Let $\hat{c}^{\partial,i}$ be the solution to the sparse regression problem corresponding to the $i$-th Poincaré derivative expansion (PoinCE-der-$i$) (5.19) with regression matrix (5.25) and data vector (5.24).<fnref n="2" /> By construction, this expansion only provides coefficients corresponding to the basis elements from the set $A^i := \{\alpha \in A : \alpha_i \ge 1\}$, since the partial derivatives w.r.t. $X_i$ of the basis elements $\{\alpha \in A : \alpha_i = 0\}$ are zero and therefore no coefficient value can be computed for those elements. Theoretically, for $\alpha \in A^i$ the coefficient $c_{\alpha}^{\partial,i}$ from (5.19) is equal to the PoinCE solution $c_{\alpha}$ from (5.18), however when estimated from a data set of finite size they will in general not coincide.

The coefficients from the set $A^i$ are sufficient for computing partial variances for variable $i$ as in (5.12) and (5.13), but not enough for computing the total variance (5.11), which requires

<footnote n="2">²Note that in practice, we normalize and rescale the regression matrix as described in Section 5.3.1 to improve the estimation of the coefficients.</footnote>
all coefficients $c_\alpha$, $\alpha \in \mathcal{A}$, and which is needed for normalizing the partial variances to Sobol' indices.

To compute the total variance from PoinCE-der expansions, we therefore aggregate the coefficients of all $d$ PoinCE-der-$i$ expansions into one vector $\hat{c}^{\partial,\text{avg}}$ as follows:

$$ \hat{c}_{\alpha}^{\partial, \text{avg}} = \frac{1}{\#\{i : \alpha_i \ge 1\}} \sum_{i:\alpha_i \ge 1} \hat{c}_{\alpha}^{\partial, i} \quad \text{for each } \alpha \in \mathcal{A} \setminus \mathbf{0}, \qquad (5.28) $$

i.e., every PoinCE-der-$i$ expansion which computed a coefficient value for the basis element with index $\alpha$ contributes equally to the averaged value. It follows that in theory, the averaged coefficient $\hat{c}_{\alpha}^{\partial,\text{avg}}$ is equal to the PoinCE solution $c_\alpha$ from (5.18), too. It can therefore be used to estimate the total variance according to (5.11).

The averaging procedure yields PoinCE-der estimates for all coefficients except for the coefficient $c_0$ corresponding to the constant term $\psi_0$. Let $\hat{c}_{\alpha}^{\partial,\text{avg}} = (\hat{c}_{\alpha}^{\partial,\text{avg}})_{\alpha \in \mathcal{A} \setminus \mathbf{0}}$ be in the form of a column vector in $\mathbb{R}^{(P-1) \times 1}$. In order to use the averaged PoinCE-der expansion also as a surrogate model, we estimate the remaining coefficient $\hat{c}_0^{\partial,\text{avg}}$ corresponding to the constant term by ordinary least-squares on the residual $y_{\text{res}}$:

$$ \begin{aligned} y_{\text{res}} &= y - \Psi \begin{pmatrix} 0 \\ \hat{c}_{\alpha}^{\partial, \text{avg}} \end{pmatrix}, \\ \hat{c}_{0}^{\partial, \text{avg}} &= \frac{1}{N} \sum_{k=1}^{N} y_{\text{res}}^{(k)} \end{aligned} $$

Note that the described construction uses model evaluations and partial derivatives separately. An obvious question is whether one could use these simultaneously to compute an estimate for the coefficients. Although tempting, the simple stacking of regression matrices $\Psi$ and $\Psi_{\partial,i}$ into a big regression matrix (as done by Peng et al. (2016) for Hermite PCE) is not satisfactory, since the increasing norm of the basis partial derivatives (see Proposition 3) introduces an undesired weighting into the problem. The simultaneous use of evaluation and derivative data is a topic of further research.

## 5.4 Numerical results

We investigate the performance of PoinCE (both based on model evaluations and on derivatives) on two numerical examples. The focus of our study is on Sobol' sensitivity analysis, but we also investigate DGSM-based upper bounds to partial variances, validation error (relative mean-squared error) and sparsity. Our implementation is based on UQLab (Marelli and Sudret, 2014) and integrates into its PCE module (Marelli et al., 2021a).

We use the following estimation techniques to compute the Sobol' indices of the models:

* PoinCE-LARS / PoinCE-der-LARS: Poincaré expansion and Poincaré derivative expansion computed by LARS as proposed in Section 5.3.3

* PoinCE-MC / PoinCE-der-MC: As a baseline, we compare to MC-based computation using the Poincaré basis/the Poincaré partial derivative basis as in Roustant et al. (2020a)
* PCE-LARS: As a second baseline, we compare to PCE computed by LARS (with generalized polynomial chaos adapted to the respective input) (Blatman and Sudret, 2011; Marelli et al., 2021a). Sparse PCE is a state-of-the-art method for computing Sobol' indices for real-world models (Le Gratiet et al., 2017).

We do not compare to any sample-based estimates of Sobol' indices, since it is known that ANOVA-based estimation outperforms sample-based estimation. For example, Sudret (2008) and Crestaux et al. (2009) have shown that polynomial chaos-based estimators of Sobol' indices are much more efficient than Monte Carlo or quasi-Monte Carlo-based estimators (for smooth models and dimensions up to 20). Recently, Becker (2020) has shown that certain sample-based approaches can be more efficient than metamodel-based ones for screening with total Sobol' indices. However, the screening performance metrics of Becker (2020) are only based on input ranking. In contrary, our practical purpose is to perform a so-called quantitative screening which aims at providing a correct screening and a good estimation of Sobol' indices.

We do not include a comparison to gradient-enhanced PCE (Peng et al., 2016; Guo et al., 2018) because so far these methods are developed only for Gaussian, uniform and Beta input and are not immediately usable for other input distributions. Furthermore, the code of the relatively involved sampling- and preconditioning approach is not readily available. The development and comparison of gradient-enhanced PoinCE to gradient-enhanced PCE is a topic of future research.

Partial variances are normalized to Sobol' indices using the total variance. For PCE-LARS and PoinCE-LARS, the total variance is computed from the expansion coefficients as in (5.11). For PoinCE-MC and PoinCE-der-MC, we use the sample variance as done by Roustant et al. (2020a). For PoinCE-der-LARS, the total variance is obtained by the procedure detailed in Section 5.3.4.

The DGSM-based upper bound to the total partial variances is computed from (5.17) using the coefficients of the PoinCE derivative expansions as described in Section 5.3.4. Note that the inequalities in (5.17) are analytical bounds that do not necessarily hold for the estimated quantities.

For uniform and Gaussian input variables, the analytical expression for the Poincaré basis functions is used, while for all others, the basis functions are computed numerically using a resolution of $10^3$ points for the uniform grid within the given bounds (see Section 5.3.1).<fnref n="3" />

### 5.4.1 Dyke cost model

Our first application is a simplified analytical model computing the cost associated to a dyke that is to be constructed along a stretch of river to prevent flooding (Iooss and Lemaître, 2015; Roustant et al., 2020a). Its output is the cost in million euros given by

$$Y = \mathbb{1}_{S>0} + \left[ 0.2 + 0.8 \left( 1 - \exp^{-\frac{1000}{S^4}} \right) \right] \cdot \mathbb{1}_{S \le 0} + \frac{1}{20} \left( 8 \cdot \mathbb{1}_{H_d \le 8} + H_d \cdot \mathbb{1}_{H_d > 8} \right) \quad (5.29)$$

where S is the maximal annual overflow and $H_d$ is the dyke height. Here, the first term represents the cost of the consequences of a flooding event, the second describes the maintenance costs,

<footnote n="3">³For the flood model, the change in the resulting Sobol' indices when instead using a grid with $10^2$ or $10^4$ points is in the order of $10^{-4}$ or $10^{-6}$, respectively.</footnote>
and the third is associated to the construction cost. S is computed from the river characteristics detailed in Table 5.1 via the 1D Saint-Venant equations under several simplifying assumptions as follows:

$$S = \left( \frac{Q}{BK_s \sqrt{\frac{Z_m - Z_v}{L}}} \right)^{0.6} + Z_v - H_d - C_b \quad (5.30)$$

The model Eq. (5.29) is continuous and piecewise $C^1$, and therefore in $H^1$. It has 8 input variables, of which $Q, K_s, Z_v$ and $H_d$ are important, and $C_b, Z_m, L$ and $B$ are unimportant (see also the last two columns of Table 5.1).

<table><caption>**Table 5.1:** Input variables to the dyke cost model and reference values of first-order and total Sobol' indices. The reference values were obtained to good precision by Monte-Carlo-based Sobol' index estimation using a large sample (Roustant et al., 2020a).</caption><thead><tr><th>Input</th><th>Function</th><th>Unit</th><th>Distribution</th><th>S<sub>i</sub></th><th>S<sub>i</sub><sup>tot</sup></th></tr></thead><tbody><tr><td>Q</td><td>Maximal annual flowrate</td><td>m<sup>3</sup>/s</td><td>Gumbel G(1013, 558) truncated to [500, 3000]</td><td>0.358</td><td>0.483</td></tr><tr><td>K<sub>s</sub></td><td>Strickler coefficient</td><td>—</td><td>Gaussian N(30, 8<sup>2</sup>) truncated to [15, +&infin;]</td><td>0.156</td><td>0.252</td></tr><tr><td>Z<sub>v</sub></td><td>River downstream level</td><td>m</td><td>Triangular T(49, 51)</td><td>0.167</td><td>0.223</td></tr><tr><td>Z<sub>m</sub></td><td>River upstream level</td><td>m</td><td>Triangular T(54, 56)</td><td>0.003</td><td>0.008</td></tr><tr><td>H<sub>d</sub></td><td>Dyke height</td><td>m</td><td>Uniform U([7, 9])</td><td>0.119</td><td>0.177</td></tr><tr><td>C<sub>b</sub></td><td>Bank level</td><td>m</td><td>Triangular T(55, 56)</td><td>0.029</td><td>0.040</td></tr><tr><td>L</td><td>Length of river stretch</td><td>m</td><td>Triangular T(4990, 5010)</td><td>0.000</td><td>0.000</td></tr><tr><td>B</td><td>River width</td><td>m</td><td>Triangular T(295, 305)</td><td>0.000</td><td>0.000</td></tr></tbody></table>

The dyke cost model has been used by Roustant et al. (2020a) to demonstrate the performance of projection-based PoinCE. We compare the new regression-based methods PoinCE-LARS and PoinCE-der-LARS with the projection-based counterparts PoinCE-MC and PoinCE-der-MC, and additionally with the standard PCE method PCE-LARS. The projection-based estimates use a basis of total degree 2, while the regression-based estimates use degree adaptivity with a degree of up to 5 (remember that for PoinCE, the degree corresponds to the ordering of the eigenfunctions by the magnitude of the eigenvalues). The experimental design is sampled by LHS with maximin distance optimization. Gradients are computed here by finite differences. For each size of the experimental design, we perform 50 independent repetitions. We display the resulting estimates in the form of boxplots. We show results only for three input variables: the most important variable $Q$, the low-importance variable $C_b$, and the unimportant variable $B$. The results for the remaining input variables can be found in Section 5.A.

### 5.4.1.1 Comparison of MC-based and regression-based computation of PoinCE(der)

First we investigate the two different ways to compute PoinCE: projection-based as in Roustant et al. (2020a) versus sparse regression-based as described in Section 5.3. Figure 5.1 and 5.2 show estimates for first-order and total Sobol' indices. We observe that in all cases the regression-based estimates have a smaller variance than the corresponding projection-based estimates.
Also, the median of the regression-based estimates is closer to the true Sobol' index value than the median of the projection-based estimates. Note that while the regression-based estimates use a degree-adaptive basis of $p \le 5$, the projection-based estimates use a fixed degree of only $p = 2$. While this choice introduces a certain bias to the projection-based estimates, a larger value for $p$ leads to unfeasibly large variance for those estimates. This is because the coefficients of higher-order terms cannot be estimated precisely with few experimental design points, which makes the overall estimate less precise.

We also observe that regression-based estimates are often clustered around the true Sobol' index already for very small experimental design sizes. This might be due to the generally smaller variance of sparse-regression-based coefficient estimates compared to MC-based estimates, as well as due to the choice of normalization factor for Sobol' indices (coefficient-based variance Eq. (5.12) vs. sample variance as described in the beginning of Section 5.4). Since regression generally leads to more precise estimates than projection, in the remainder of this paper we focus on regression-based PoinCE estimates.

Furthermore, as already observed by Roustant et al. (2020a), PoinCE-der estimates for Sobol' indices have a smaller variance than PoinCE estimates. In the case of projection-based estimates, this is the case if the derivative has a smaller variance than the original model. In the case of regression, the explanation might be that PoinCE-$i$-der has to compute less coefficients than PoinCE for the same number of experimental design points ($\{\alpha \in A : \alpha_i > 0\}$ vs. $A$), which can result in a more precise estimate of the true coefficient values.

<figure><img src="image_84.png" /><figcaption>**Figure 5.1:** Comparison of PoinCE estimates of first-order Sobol' indices for the dyke cost model. Degree $p = 2$ for the MC-based estimates and $p \le 5$ (degree-adaptive) for the regression-based estimates. Results for the remaining variables are displayed in Figure 5.A.1 in the appendix.</figcaption></figure>

### 5.4.1.2 Comparison of regression-based PoinCE-der with PCE and the DGSM-based upper bound

Next, we investigate the performance of regression-based PoinCE compared to state-of-the-art PCE, and the usefulness of the DGSM-based upper bound to partial variances derived in (5.17). The corresponding results, unnormalized$^4$ estimates for first-order and total Sobol' indices, are

<footnote>$^4$We show unnormalized indices because the DGSM-based upper bound is not normalized.</footnote>
<figure><img src="image_85.png" /><figcaption>**Figure 5.2:** Comparison of PoinCE estimates of total Sobol' indices for the dyke cost model. Degree $p = 2$ for the MC-based estimates and $p \le 5$ (degree-adaptive) for the regression-based estimates. Results for the remaining variables are displayed in Figure 5.A.2 in the appendix.</figcaption></figure>

displayed in Figure 5.3 and 5.4. Because PoinCE-der achieves more accurate estimates than PoinCE, we compute the DGSM-based upper bound using the PoinCE-der-$i$ coefficients. For total Sobol' indices, we also include a precise Monte Carlo estimate for the DGSM-based upper bound (using $10^7$ derivative samples) computed from (5.9) and the second inequality of (5.17).

We make the following observations: the PCE-LARS estimates are generally very similar to the PoinCE-LARS estimates, but the latter often have a slightly larger range. The similarity might be because both rely on model evaluations only. However, the respective basis functions have a very different shape (for inputs that do not follow a Gaussian distribution). In particular, the PoinCE basis functions by construction obey Neumann boundary conditions, i.e., have zero derivative on the boundary.

As observed before for normalized indices, PoinCE-der performs better than PoinCE: the median is closer to the true value, and the range is smaller. This effect is especially pronounced for low-importance variables. In 8 dimensions, a PoinCE-der expansion of degree 5 has 495 terms while the total-degree basis of PCE and PoinCE has 1287 terms. This means that a PoinCE-der expansion has to estimate less than half of the coefficients. PoinCE-der generally gives a tighter "lower bound" than PCE (but note that the estimates are not guaranteed to be a lower bound).

By construction (5.17), the DGSM-based upper bound estimate is larger than or equal to the corresponding total Sobol' index estimate. However, it would be an upper bound to the true Sobol' index value only if the full infinite expansion was used. This is visible in Figure 5.A.4 and 5.4: for some inputs, the upper bound estimate almost coincides with the Sobol' index estimate, and is smaller than the true Sobol' index.

For some inputs, such as $K_s$ and especially $H_d$ (see Figure 5.A.4), the DGSM-based upper bound is not tight. An explanation for the large gap between the Sobol' index and the upper bound for $H_d$, using (5.17), might be that $H_d$ is responsible for a kink in the model (through the last term of (5.29)), while the Poincaré basis is differentiable in $H_d$ (cosine functions). Therefore, its expansion needs high-order terms. From (5.17) it follows that even if the PoinCE-der found the true coefficients, the estimated upper bound is not tight due to the quickly-growing eigenvalues.
Finally, we see that in many cases the gap between the estimated upper bound estimate and its true value is larger than the gap between Sobol' index estimate and its true value. The reason might be that some higher-order terms are still missing from the considered expansion. Due to the eigenvalue factor involved in the estimate of the upper bound (5.17), this has a larger influence on the upper bound than on the Sobol' index estimate.

<figure><img src="image_86.png" /><figcaption>**Figure 5.3:** Estimates of unnormalized first-order Sobol' indices for the dyke cost model ($p \le 5$). Boxplots: in grey the PCE-based estimates. The dashed line ("True value") denotes a high-precision estimate for the unnormalized first-order Sobol' index. Results for the remaining input variables can be found in Figure 5.A.3 in Section 5.A.</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>**Figure 5.4:** Estimates of unnormalized total Sobol' indices for the dyke cost model ($p \le 5$). Boxplots: in grey the PCE-based estimates and in black the DGSM-based upper bound from (5.17). Lines: the dashed line ("True value") denotes a high-precision estimate for the unnormalized total Sobol' index, while the dotted line ("UB true value") is a MC-based high-precision estimate for the DGSM-based upper bound. Results for the remaining input variables can be found in Figure 5.A.4 in Section 5.A.</figcaption></figure>

### 5.4.1.3 Comparison of total variance, relative mean-squared error, and sparsity

To estimate Sobol' indices precisely, it is crucial to have a good estimate for the total variance. For PCE-LARS, PoinCE-LARS, and PoinCE-der-LARS, this value can directly be computed
from the expansion coefficients, while for PoinCE computed by projection, we are using the empirical variance, as detailed in the beginning of Section 5.4. In Figure 5.5a we display the scatter of the variance estimates (50 replications). The empirical estimate has the largest variation, while PoinCE-der-LARS has the smallest. PCE-LARS and PoinCE-LARS underestimate the total variance more than PoinCE-der-LARS. This is likely one reason for the good performance of PoinCE-der for the estimation of Sobol' indices: a more accurate total variance leads to more accurate Sobol' indices.

Interestingly, while PoinCE performs well for the estimation of Sobol' indices, this is not true for the generalization error, given by the relative mean-squared error

$$ \text{RelMSE} = \frac{\mathbb{E}_X [(f(X) - f^{\text{surr}}(X))^2]}{\text{Var}_X [f(X)]} \qquad (5.31) $$

with the surrogate model $f^{\text{surr}}$. The RelMSE is computed by Monte Carlo integration on a validation set of size $10^6$ sampled from the input distribution $\mu$. In Figure 5.5b we display box-plots of estimates for the generalization error on a validation set of size $10^6$ (mean-squared error normalized by the variance of the validation set). PoinCE-der attains a smaller relative MSE than PoinCE. PCE shows faster convergence behavior than both, and attains a smaller relative MSE than PoinCE. PoinCE-der performs better than PCE for the two small experimental design sizes, which shows that the information brought by derivatives might be especially useful when data is scarce.

Finally, we display the number of nonzero coefficients of each expansion in Figure 5.5c. PCE and PoinCE have a similar number of active coefficients, while PoinCE-der has considerably more active coefficients. This corresponds to a better validation error only for the smallest experimental design size (Figure 5.5b). A possible explanation is the following: using the derivative information to compute PoinCE-der, the same amount of data is used for a (truncated) expansion containing less terms than the PoinCE expansion has (since some terms, which do not contain the variable for which the partial derivative is taken, drop out). In this way, more coefficients can be computed.

## 5.4.2 Mascaret data set

Our second application focuses on a phenomenological and industrial simulation model, called Mascaret (Goutal et al., 2012), based on a 1D solver of the Saint Venant equations and aiming at computing water height for river flood events. The studied case, taken from Petit et al. (2016) and also studied in Roustant et al. (2017), is the French Vienne river in permanent regime whose uncertain input data concern flowrate, several physical parameters and geometrical data (transverse river profiles). 37 independent inputs have then been considered as random variables:

* 12 Strickler coefficients of the main channel $K_{s,c}^i$, uniform in [20, 40];

* 12 Strickler coefficients of the flood plain $K_{s,p}^i$, uniform in [10, 30];

* 12 slope perturbations $dZ^i$, standard Gaussian with bounds [-3, 3];

* 1 discharge value $Q$, Gaussian with zero mean and standard deviation 50, bounds [-150, 50].

The derivatives of the model output with respect to these 37 inputs have been efficiently (with a cost independent of the number of inputs) computed by using the adjoint model of Mascaret
<figure><img src="image_88.png" /><figcaption>**Figure 5.5:** Dyke cost model: Comparison of PCE and PoinCE(-der) with respect to the following metrics: estimation of total variance, relative mean-squared error, and number of nonzero coefficients in the respective expansions.</figcaption></figure>

(Demangeon et al., 2015). This adjoint model has been obtained by automatic differentiation (Griewank and Walther, 2008) using the automatic differentiation software Tapenade (Hascoët and Pascual, 2013). A large-size Monte Carlo sample ($n = 20000$) is available from the study of Petit et al. (2016). This data set contains all the values of the 37 inputs, the water height as output and the 37 partial derivatives of the output (one derivative with respect to each input). Note that this sample, which has a very large size, has been obtained during a research work for a demonstrative purpose. In industrial practice, the aim is to use the minimal possible sample size: it is expected to use methods able to deal with sample sizes of the order of one hundred.

Previous studies on this data set (Petit et al., 2016; Roustant et al., 2017) have identified 32 of the 37 inputs as noninfluential. In our study, we display results for the 5 remaining inputs ($K_{s,c}^{11}$, $K_{s,c}^{12}$, $dZ^{11}$, $dZ^{12}$, and $Q$) and for one of the noninfluential inputs ($K_{s,c}^1$). We choose a basis with hyperbolic truncation using $q = 0.5$, and degree adaptivity $p = 1, 2, \dots, 8$. We analyze several experimental design sizes ranging from 30 to 300. For each experimental design size, we run 30 replications, sampling the design randomly without replacement from the given full data set. “True” values for Sobol' indices and total variance are computed from a PCE using all 20000 points.

### 5.4.2.1 Comparison of regression-based PoinCE(der) with PCE and the DGSM-based upper bound

Estimates of first-order and total Sobol' indices are displayed in Figure 5.6 and 5.7. We display results for regression-based PCE, PoinCE, and PoinCE-der. In addition, we display the upper bound computed as in (5.17), computed based on PoinCE-der coefficients and normalized by the PoinCE-der total variance. We observe that for the non-influential variable $K_{s,c}^1$ (and indeed all other 31 non-influential variables), derivative-based PoinCE correctly identify a total and first-order Sobol' index of 0. Overall, PoinCE and PCE show very similar results, with PoinCE
having slightly larger variance in a few cases. For some variables such as $dZ^{11}$ and $Q$, the DGSM-based estimate of the upper bound almost coincides with the PoinCE-der estimate. Overall, we observe that PoinCE-der estimates have smaller variance than PCE and PoinCE for the important variables $K_{s,c}^{11}$, $dZ^{11}$, and $Q$, even already for 30 experimental design points. For low-importance variables such as $K_{s,c}^{12}$ and $dZ^{12}$, PoinCE-der correctly identifies a value away from zero already for the smallest experimental design, while half of the PCE and PoinCE estimates are zero. For small experimental design sizes, the PoinCE-der estimates also have a smaller bias than the PCE and PoinCE estimates. Sometimes the PoinCE-der estimates seem to systematically over- or underestimate the true Sobol' index by a small amount. However, note that the "true" value was computed by a PCE (based on all 20 000 points), and might therefore itself be slightly inaccurate.

<figure><img src="image_89.png" /><figcaption>Figure 5.6: First-order Sobol' indices for the Mascaret data set (30 replications). "True" values computed from a PCE using all 20 000 points.</figcaption></figure>

### 5.4.2.2 Comparison of total variance and relative mean-squared error

In Figure 5.8a we display various estimates for the total variance. We observe again that PoinCE-der yields an estimate with smaller variance and less bias than PoinCE and PCE. PoinCE and PCE both have smaller variance than the empirical estimate, but generally underestimate the total variance.

While PoinCE-der estimates Sobol' indices and total variance well, we observe in Figure 5.8b showing the relative MSE that PCE and PoinCE are performing better as global surrogate
<figure><img src="image_90.png" /><figcaption>**Figure 5.7:** Total Sobol' indices for the Mascaret data set (30 replications). "True" values computed from a PCE using all 20 000 points.</figcaption></figure>

models: their model approximation error is for large experimental designs almost an order of magnitude better than for PoinCE-der.

## 5.5 Conclusion

In this paper we studied PoinCE, an expansion in terms of the Poincaré basis, which is an orthonormal basis of $L^2(\mu)$ with the unique property that all its partial derivatives form again an orthogonal basis for the same space. We provided a proof of this property as well as a few analytical results as direct consequences. In particular, we showed how upper and lower bounds for partial variances can be obtained analytically from PoinCE coefficients.

We described the computation of PoinCE and Poincaré derivative expansions by sparse regression and applied the method to two numerical examples. We found that while PoinCE does not outperform PCE in terms of validation error, it can be advantageous for estimating Sobol' indices in the low-data regime. PoinCE is therefore a valuable tool if model derivatives are cheaply available (e.g., by automatic differentiation or as a by-product of the simulation). Taking partial derivatives reduces the size of the truncated basis especially for high-dimensional, low-order total-degree bases, which gives an advantage to derivative-based PoinCE over expansions relying on model evaluations.

Future work on the topic of PoinCE will investigate the simultaneous use of model evaluations
<figure><img src="image_91.png" /><figcaption>**Figure 5.8:** Mascaret data set. Comparison of PCE and PoinCE(-der) with respect to the following metrics: estimation of total variance, and relative mean-squared error. “True” value of total variance computed from a PCE using all 20 000 points.</figcaption></figure>

and derivatives for the computation of the coefficients, and compare to the related topic of
gradient-enhanced PCE.

Funding

This paper is a part of the project "Surrogate Modeling for Stochastic Simulators (SAMOS)"
funded by the Swiss National Science Foundation (Grant #200021_175524), whose support is
gratefully acknowledged. Part of this research was conducted within the frame of the Chair in
Applied Mathematics OQUAIDO, gathering partners in technological research (BRGM, CEA,
IFPEN, IRSN, Safran, Storengy) and academia (CNRS, Ecole Centrale de Lyon, Mines Saint-
Etienne, University of Grenoble, University of Nice, University of Toulouse) around advanced
methods for computer experiments. Support from the ANR-3IA Artificial and Natural Intelli-
gence Toulouse Institute is gratefully acknowledged.

Declaration of competing interests

No competing interests.

References

Antoniadis, A. (1984). Analysis of variance on function spaces. *Statistics*, 15(1):59-71.

Bakry, D., Gentil, I., and Ledoux, M. (2014). *Analysis and geometry of Markov diffusion operators*, volume 348 of *Fundamental Principles of Mathematical Sciences*. Springer, Cham.

Becker, W. (2020). Metafunctions for benchmarking in sensitivity analysis. *Reliab. Eng. Sys. Safety*, 204:107189.
Blatman, G. and Sudret, B. (2008). Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach. *C. R. Mécanique*, 336(6):518-523.

Blatman, G. and Sudret, B. (2011). Adaptive sparse polynomial chaos expansion based on Least Angle Regression. *J. Comput. Phys*, 230:2345-2367.

Borgonovo, E. and Plischke, E. (2016). Sensitivity analysis: A review of recent advances. *Eur. J. Oper. Res.*, 248(3):869-887.

Bujurke, N. M., Salimath, C. S., and Shiralashetti, S. C. (2008). Computation of eigenvalues and solutions of regular Sturm-Liouville problems using Haar wavelets. *J. Comput. Appl. Math.*, 219(1):90-101.

Candès, E. J. and Plan, Y. (2011). A probabilistic and RIPless theory of compressed sensing. *IEEE Trans. Inform. Theory*, 57(11):7235-7254.

Candès, E. J. and Wakin, M. B. (2008). An introduction to compressive sampling: A sensing/sampling paradigm that goes against the common knowledge in data acquisition. *IEEE Signal Process. Mag.*, 25(2):21-30.

Chapelle, O., Vapnik, V., and Bengio, Y. (2002). Model selection for small sample regression. *Machine Learning*, 48(1):9-23.

Constantine, P. G., Eldred, M. S., and Phipps, E. T. (2012). Sparse pseudospectral approximation method. *Comput. Methods Appl. Mech. Engrg.*, 229:1-12.

Crestaux, T., Le Maître, O., and Martinez, J.-M. (2009). Polynomial chaos expansion for sensitivity analysis. *Reliab. Eng. Sys. Safety*, 94(7):1161-1172.

Demangeon, F., Goeury, C., Zaoui, F., Goutal, N., Pascual, V., and Hascoët, L. (2015). Algorithmic differentiation applied to the optimal calibration of a shallow water model. *La Houille Blanche - Revue internationale de l'eau*, 102(4):57-65.

Efron, B. and Stein, C. (1981). The Jacknife estimate of variance. *Ann. Stat.*, 9(3):586-596.

Ernst, O., Mugler, A., Starkloff, H.-J., and Ullmann, E. (2012). On the convergence of generalized polynomial chaos expansions. *ESAIM: Math. Model. Numer. Anal.*, 46(02):317-339.

Fang, K.-T., Li, R., and Sudjianto, A. (2006). *Design and Modeling for Computer Experiments*. Computer Science and Data Analysis Series. Chapman & Hall/CRC, Boca Raton, Florida, USA.

Gejadze, I., Malaterre, P.-O., and Shutyavtsev, V. (2019). On the use of derivatives in the polynomial chaos based global sensitivity and uncertainty analysis applied to the distributed parameter models. *J. Comput. Phys.*, 381:218-245.

Ghanem, R. G. and Spanos, P. (1991). *Stochastic finite elements – A spectral approach*. Springer Verlag, New York. (Reedited by Dover Publications, Mineola, 2003).
Goutal, N., Lacombe, J.-M., Zaoui, F., and El-Kadi-Abderrezak, K. (2012). MASCARET: a 1-D open-source software for flow hydrodynamic and water quality in open channel networks. In Murillo Muñoz, R., editor, *River Flow 2012: Proceedings of the International Conference on Fluvial Hydraulics*, volume 2, pages 1169–1174, San José, Costa Rica. CRC Press.

Griewank, A. and Walther, A. (2008). *Evaluating derivatives: Principles and techniques of automatic differentiation*. SIAM Philadelphia.

Guo, L., Narayan, A., and Zhou, T. (2018). A gradient enhanced $ℓ^1$-minimization for sparse approximation of polynomial chaos expansions. *J. Comput. Phys.*, 367:49–64.

Hascoët, L. and Pascual, V. (2013). The Tapenade automatic differentiation tool: Principles, model and specification. *ACM T. Math. Software*, 39(3):1–43.

Hoeffding, W. (1948). A class of statistics with asymptotically normal distributions. *Ann. Math. Stat.*, 19:293–325.

Homma, T. and Saltelli, A. (1996). Importance measures in global sensitivity analysis of non linear models. *Reliab. Eng. Sys. Safety*, 52:1–17.

Iooss, B. and Lemaître, P. (2015). *Uncertainty Management in Simulation-Optimization of Complex Systems*, volume 59 of *Operations Research/Computer Science Interfaces Series (ORCS)*, chapter A review on global sensitivity analysis methods. Springer.

Jakeman, J. D., Eldred, M. S., and Sargsyan, K. (2015). Enhancing $ℓ_1$-minimization estimates of polynomial chaos expansions using basis selection. *J. Comput. Phys.*, 289:18–34.

Kougioumtzoglou, I. A., Petromichelakis, I., and Psaros, A. F. (2020). Sparse representations and compressive sampling approaches in engineering mechanics: A review of theoretical concepts and diverse applications. *Prob. Eng. Mech.*, 61:103082.

Kucherenko, S. and Iooss, B. (2017). Derivative-based global sensitivity measures. In *Handbook of Uncertainty Quantification*, pages 1241–1263. Springer International Publishing.

Kucherenko, S., Rodriguez-Fernandez, M., Pantelides, C., and Shah, N. (2009). Monte Carlo evaluation of derivative-based global sensitivity measures. *Reliab. Eng. Syst. Safety*, 94:1135–1148.

Kwon, K. H. and Lee, D. W. (2003). Orthogonal functions satisfying a second-order differential equation. *J. Comput. Appl. Math.*, 153(1-2):283–293.

Lamboni, M., Iooss, B., Popelin, A.-L., and Gamboa, F. (2013). Derivative-based global sensitivity measures: general links with Sobol' indices and numerical tests. *Math. Comput. Simul.*, 87:44–54.

Le Gratiet, L., Marelli, S., and Sudret, B. (2017). *Metamodel-based sensitivity analysis: polynomial chaos expansions and Gaussian processes*, pages 1289–1325. Springer.

Le Maître, O., Reagan, M., Najm, H., and Ghanem, R. (2002). A stochastic projection method for fluid flow - II. random process. *J. Comp. Phys.*, 181:9–44.
Li, Y., Anitescu, M., Roderick, O., and Hickernell, F. (2011). Orthogonal bases for polyno-
mial regression with derivative information in uncertainty quantification. *Int. J. Uncertainty
Quantification*, 1:297–320.

Lüthen, N., Marelli, S., and Sudret, B. (2021). Sparse polynomial chaos expansions: Literature survey and benchmark. *SIAM/ASA J. Uncertain. Quantif.*, 9(2):593–649.

Lüthen, N., Marelli, S., and Sudret, B. (2022). Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications. *Int. J. Uncertainty Quantification*, 12(3):49–74.

Marelli, S., Lüthen, N., and Sudret, B. (2021). UQLab user manual – Polynomial chaos expansions. Technical report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich, Switzerland. Report # UQLab-V1.4-104.

Marelli, S. and Sudret, B. (2014). UQLab: A framework for uncertainty quantification in Matlab. In *Vulnerability, Uncertainty, and Risk (Proc. 2nd Int. Conf. on Vulnerability, Risk Analysis and Management (ICVRAM2014), Liverpool, United Kingdom)*, pages 2554–2563.

Matthies, H. G. and Keese, A. (2005). Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. *Comput. Methods Appl. Mech. Engrg.*, 194:1295–1331.

McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics*, 2:239–245.

Mikolas, M. (1955). Über gewisse Eigenschaften orthogonaler Systeme der Klasse L² und die Eigenfunktionen Sturm-Liouvillescher Differentialgleichungen. *Acta Math. Acad. Sci. H.*, 6:147–190.

Oladyshkin, S. and Nowak, W. (2012). Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion. *Reliab. Eng. Syst. Saf.*, 106:179–190.

Peng, J., Hampton, J., and Doostan, A. (2016). On polynomial chaos expansion via gradient-enhanced $\ell_1$-minimization. *J. Comput. Phys.*, 310:440–458.

Petit, S., Zaoui, F., Popelin, A.-L., Goeury, C., and Goutal, N. (2016). Couplage entre indices à base de dérivées et mode adjoint pour l'analyse de sensibilité globale. Application sur le code Mascaret. Technical report. Preprint.

Prieur, C. and Tarantola, S. (2017). Variance-based sensitivity analysis: Theory and estimation algorithms. In *Springer Handbook on Uncertainty Quantification*, pages 1217–1239. Springer.

Razavi, S., Jakeman, A., Saltelli, A., Prieur, C., Iooss, B., Borgonovo, E., Plischke, E., Lo Piano, S., Iwanaga, T., Becker, W., Tarantola, S., Guillaume, J., Jakeman, J., Gupta, H., Melillo, N., Rabiti, G., Chabridon, V., Duan, Q., Sun, X., Smith, S., Sheikholeslami, R., Hosseini, N., Asadzadeh, M., Puy, A., Kucherenko, S., and Maier, H. (2021). The future of sensitivity analysis: An essential discipline for systems modelling and policy making. *Environ. Modell. Softw.*, 137(104954).
Roderick, O., Anitescu, M., and Fischer, P. (2010). Polynomial regression approaches using derivative information for uncertainty quantification. *Nucl. Sci. Eng.*, 164:122–139.

Roustant, O., Barthe, F., and Iooss, B. (2017). Poincaré inequalities on intervals – application to sensitivity analysis. *Electron. J. Stat.*, 11(2):3081–3119.

Roustant, O., Gamboa, F., and Iooss, B. (2020). Parseval inequalities and lower bounds for variance-based sensitivity indices. *Electron. J. Statist.*, 14(1):386–412.

Simon, B. (2010). *Szegő's theorem and its descendants: spectral theory for L2 perturbations of orthogonal polynomials*. Princeton University Press.

Smith, R. C. (2014). *Uncertainty Quantification: Theory, Implementation and Applications*. SIAM Computational Science and Engineering.

Sobol, I. and Gresham, A. (1995). On an alternative global sensitivity estimator. In *Proceedings of SAMO 1995*, pages 40–42, Belgirate.

Sobol', I. M. (1993). Sensitivity estimates for nonlinear mathematical models. *Math. Modeling & Comp. Exp.*, 1(4):407–414.

Sobol', I. M. and Kucherenko, S. (2009). Derivative based global sensitivity measures and their link with global sensitivity indices. *Math. Comput. Simul.*, 79(10):3009–3017.

Sudret, B. (2006). Global sensitivity analysis using polynomial chaos expansions. In Spanos, P. and Deodatis, G., editors, *Proc. 5th Int. Conf. on Comp. Stoch. Mech (CSM5)*, Rhodos, Greece.

Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. *Reliab. Eng. Sys. Saf.*, 93:964–979.

Sudret, B. and Mai, C. V. (2015). Computing derivative-based global sensitivity measures using polynomial chaos expansions. *Reliab. Eng. Sys. Safety*, 134:241–250.

Szegő, G. (1939). *Orthogonal polynomials*, volume 23. American Mathematical Society.

Torre, E., Marelli, S., Embrechts, P., and Sudret, B. (2019). Data-driven polynomial chaos expansion for machine learning regression. *J. Comput. Phys*, 388:601–623.

Wiener, N. (1938). The homogeneous chaos. *American Journal of Mathematics*, 60(4):897–936.

Xiu, D. and Karniadakis, G. E. (2002). The Wiener-Askey polynomial chaos for stochastic differential equations. *SIAM J. Sci. Comput.*, 24(2):619–644.

Zettl, A. (2010). *Sturm-Liouville theory*. Number 121. American Mathematical Society.
# Appendix

## 5.A Additional results

In Figure 5.A.1 to 5.A.4, we show additional results for the dyke cost model, namely Sobol' index estimates (normalized and unnormalized) for the remaining five input variables. For the corresponding discussion, see Section 5.4.1.

<figure><img src="image_92.png" /><figcaption>**Figure 5.A.1:** Comparison of PoinCE estimates of first-order Sobol' indices for the dyke cost model. Degree $p = 2$ for the MC-based estimates and $p \le 5$ (degree-adaptive) for the regression-based estimates. See also Figure 5.1 in the main part of the paper.</figcaption></figure>
<figure><img src="image_93.png" /><figcaption>**Figure 5.A.2:** Comparison of PoinCE estimates of total Sobol' indices for the dyke cost model. Degree $p = 2$ for the MC-based estimates and $p \le 5$ (degree-adaptive) for the regression-based estimates. See also Figure 5.2 in the main part of the paper.</figcaption></figure>
<figure><img src="image_94.png" /><figcaption>**Figure 5.A.3:** Estimates of unnormalized first-order Sobol' indices for the dyke cost model ($p \le 5$). Boxplots: in grey the PCE-based estimates. The dashed line ("True value") denotes a high-precision estimate for the unnormalized first-order Sobol' index. See also Figure 5.3.</figcaption></figure>
<figure><img src="image_95.png" /><figcaption>**Figure 5.A.4:** Estimates of unnormalized total Sobol' indices for the dyke cost model ($p \le 5$). Boxplots: in grey the PCE-based estimates and in black the DGSM-based upper bound from (5.17). Lines: the dashed line (“True value”) denotes a high-precision estimate for the unnormalized total Sobol' index, while the dotted line (“UB true value”) is a MC-based high-precision estimate for the DGSM-based upper bound. See also Figure 5.4.</figcaption></figure>
CHAPTER 6

# A spectral surrogate model for stochastic simulators computed from trajectory samples

This chapter contains the preprint of

Lüthen, N., Marelli, S., Sudret, B. (2022). A spectral surrogate model for stochastic simulators computed from trajectory samples. https://arxiv.org/abs/2207.05653

submitted for publication in the journal *Computer Methods in Applied Mechanics and Engineering* in July 2022.

**Author contributions.** N. Lüthen: Conceptualization, Methodology, Software, Investigation, Writing - Original Draft, Visualization. S. Marelli: Conceptualization, Methodology, Writing - Review & Editing, Supervision. B. Sudret: Conceptualization, Writing - Review & Editing, Supervision, Funding acquisition.

## Abstract

Stochastic simulators are non-deterministic computer models which provide a different response each time they are run, even when the input parameters are held at fixed values. They arise when additional sources of uncertainty are affecting the computer model, which are not explicitly modeled as input parameters. The uncertainty analysis of stochastic simulators requires their repeated evaluation for different values of the input variables, as well as for different realizations
of the underlying latent stochasticity. The computational cost of such analyses can be consid-
erable, which motivates the construction of surrogate models that can approximate the original
model and its stochastic response, but can be evaluated at much lower cost.

We propose a surrogate model for stochastic simulators based on spectral expansions. Con-
sidering a certain class of stochastic simulators that can be repeatedly evaluated for the same
underlying random event, we view the simulator as a random field indexed by the input pa-
rameter space. For a fixed realization of the latent stochasticity, the response of the simulator
is a deterministic function, called trajectory. Based on samples from several such trajectories,
we approximate the latter by sparse polynomial chaos expansion and compute analytically an
extended Karhunen-Loève expansion (KLE) to reduce its dimensionality. The uncorrelated but
dependent random variables of the KLE are modeled by advanced statistical techniques such
as parametric inference, vine copula modeling, and kernel density estimation. The resulting
surrogate model approximates the marginals and the covariance function, and allows to obtain
new realizations at low computational cost. We observe that in our numerical examples, the
first mode of the KLE is by far the most important, and investigate this phenomenon and its
implications.

6.1 Introduction

Nowadays, computer simulations are an essential ingredient of the research and development
workflow in virtually all fields of science and engineering. Typically, not all parameters and
conditions needed for the simulations are known exactly, and this uncertainty affects the output
of the simulations. This is the main focus of the field of uncertainty quantification (Smith, 2014).

Most computer simulations can be classified as *deterministic simulators*: repeatedly evaluating
the model $\mathcal{M}$ for the same set of input parameters $\boldsymbol{x}$ always yields the same deterministic
response $y = \mathcal{M}(\boldsymbol{x}) \in \mathbb{R}^1$. To perform uncertainty quantification, the uncertainty on the input
 parameters and conditions) is usually represented probabilistically, and we follow this approach
in this paper. Propagating the input uncertainty through the deterministic simulator, the overall
response of the simulation becomes a random quantity.

However, not all computer simulations can be classified as deterministic simulators. Some models
contain intrinsic stochasticity that cannot be modeled as input parameter, e.g., epidemiological
models where each transmission or recovery is a random event, governed by the respective rate
of occurrence. Other models depend on uncontrollable environmental variables such as wind
fields or earthquakes, for which it can be infeasible or undesirable to explicitly model their
uncertainty. In these cases, it is more convenient to use the notion of a *stochastic simulator*:
only some of the uncertainty is explicitly modeled as random input variables, and there is some
residual randomness affecting the computational model that causes the model response $\mathcal{M}(\boldsymbol{x})$
for a fixed set of input parameters $\boldsymbol{x}$ to still be a random variable: $Y_\boldsymbol{x} = \mathcal{M}(\boldsymbol{x})$. In other words,
evaluating the computer model several times with the same input parameters $\boldsymbol{x}$ will result in

<footnote n="1">¹We consider here only real-valued simulators. The extension to low-dimensional vector-valued simulators is straightforward. For the extension to high-dimensional vector-valued or function-valued simulators, see e.g. Nagel et al. (2020); Perrin et al. (2021).</footnote>
different realizations *y* of the random variable *Y*<sub>*x*</sub>. Of course, since there is no true randomness in a computer, every computer simulation can be made deterministic by fixing the random seed. However, the seed is in general not a useful parametrization of uncertainty.

Uncertainty quantification methods typically require many runs of the computational model, which can become costly or even infeasible for expensive engineering simulators. To save computational resources, the model is often replaced with a cheaper surrogate model (or metamodel), which provides a reasonably good approximation to the original model. The surrogate model is computed from a small number of model evaluations and can subsequently be evaluated many times with negligible computational cost. Surrogate models often treat the model as a black box, i.e., they do not use any specific knowledge about the model and rely only on the available input-output data samples (and sometimes on the characteristics of the input parameter space). Popular surrogate models for deterministic simulators include polynomial chaos expansions (Ghanem and Spanos, 1991; Xiu and Karniadakis, 2002), Kriging (Sacks et al., 1989; Rasmussen and Williams, 2006), radial basis functions (Buhmann, 2000), and support vector regression (Vapnik, 1995; Smola and Schölkopf, 2004).

Since the response of stochastic simulators is a random variable for every set of input parameters, even more runs might be required to analyze their uncertainty, making surrogate models all the more relevant in this case. Research on surrogating stochastic simulators is comparatively recent. Most available methods focus on the marginal response distribution $P(Y|\mathbf{X}=\mathbf{x})$ for $\mathbf{x} \in \mathcal{D}$ and emulate the conditional density itself or certain statistics of it. Early contributions aimed at characterizing the variation of the first two moments of the output response over the input domain using joint Gaussian process models (Iooss and Ribatet, 2009; Marrel et al., 2012). Another class of methods aims at directly modeling the variation of the marginal output probability density function (pdf) of the random variable $Y_x$ over the input domain. Assuming that the true marginal response pdf at a number of input locations is known, Moutoussamy et al. (2015) represent the marginal pdf of a new input point as a linear combination of training examples (i.e., kernel regression) or of specifically constructed basis functions. However, the true marginal pdf is rarely known or its generation might require a lot of samples. For a finite number of stochastic simulator evaluations over the input domain (with or without replications), Zhu and Sudret (2020, 2021b) model the variation of the marginal output pdf over the input domain using the so-called generalized lambda model, a parametric distribution family that is able to approximate many classical families. In fact, stochastic simulators are akin to real-world scientific experiments, which are usually stochastic due to unavoidable measurement error and environmental noise. Therefore, standard statistical methods like quantile regression (Torossian et al., 2020) and kernel conditional density estimation (Hall et al., 2004) can also be used to emulate the marginal distribution of the response of a stochastic simulator. Furthermore, Zhu and Sudret (2022) developed an approach that emulates the stochastic simulator response in distribution, inspired by the weak PCE methodology based on maximum likelihood estimation (Xiu, 2010).

A related method from machine learning are Bayesian neural networks, whose weights are modeled as independent Gaussian random variables (MacKay, 1992; Goan and Fookes, 2020). Bayesian methods such as Markov Chain Monte Carlo or variational inference are used to de-
termine the parameters of the weight densities from the given data. Furthermore, generative
models like variational autoencoders (Kingma and Welling, 2014) and generative adversarial
networks (Goodfellow et al., 2014) can be seen as surrogate models in distribution, learning a
conditional target density from data.

All the methods cited above aim at emulating only the univariate probability density functions of the response random variables of the stochastic simulator. However, they do not take into account the correlation and higher-order information between the stochastic simulator responses at different points in the input domain. This close relation between the responses at different input locations can be best illustrated by fixing the stochasticity of the simulator (e.g., by fixing the random seed)<fnref n="2" />: in this case, the stochastic simulator response over the input domain becomes a deterministic function, which we call a *trajectory*. In other words, the stochastic simulator can be seen as a random field, i.e., as a collection of random functions.

Surrogating a stochastic simulator based on few model evaluations becomes therefore the task of inferring a random field from discrete samples (often called “limited data”). Popular methods for modeling random fields include orthogonal series expansions, such as spectral representation (Shinozuka and Deodatis, 1991; Grigoriu, 1993) or Karhunen-Loève expansion (KLE) (Loève, 1978; Karhunen, 1946; Zhang and Ellingwood, 1994; Ramsay and Silverman, 2005; Grigoriu, 2006), and translation processes, which are mappings of Gaussian processes (Yamazaki and Shinozuka, 1988; Grigoriu, 1998; Sakamoto and Ghanem, 2002; Shields et al., 2011). To our knowledge, the only publication in the specific context of stochastic simulators which takes the random field point of view and aims at emulating trajectories (including the higher-order relations between responses at different input locations) is by Azzi et al. (2019), who construct a metamodel using Karhunen-Loève expansion together with the deterministic methods PCE and Kriging.

The goal of our paper is to develop a surrogate model that is able to emulate the trajectories of a
stochastic simulator, and allows insight into the dependence between the simulator responses at
different input locations. Our method of choice in this paper is Karhunen-Loève expansion, one
of the most popular methods for random field inference from limited data. The main challenges
in constructing a trajectory-based surrogate for a stochastic simulator (a *stochastic emulator*)
are explained in more detail in the following:

1. Accuracy and efficiency: the surrogate should be accurate while needing as few model evaluations as possible.

2. Continuous surrogate from discrete data: the surrogate should emulate the response over the whole (continuous) input domain, while the available data consists of trajectories sampled at a few points throughout the input domain (i.e., discrete samples).

3. The stochastic simulator is in general a non-Gaussian random field. This introduces additional complexity into the Karhunen-Loève model.

We are addressing each of these challenges by introducing a novel approach that combines several
state-of-the-art methodologies. We use Karhunen-Loève expansion in conjunction with sparse
PCE (Blatman and Sudret, 2011; Lüthen et al., 2021), which is a powerful and sample-efficient

<footnote n="2">²Note that this does not require this randomness to be modeled. In practice, fixing the seed might not be
possible for all computational models, as it depends on their implementation.</footnote>
surrogate modeling method for deterministic simulators, to address Challenge 1. This circumvents the otherwise high computational cost of solving the integral eigenvalue problem of KLE (Schwab and Todor, 2006; Betz et al., 2014) by reducing the integral eigenvalue problem to finite-dimensional discrete principal component analysis (PCA) in the truncated space of PCE coefficients. The joint distribution of the resulting sample of dependent random KLE coefficients (Challenge 3) is identified using statistical inference within the marginal-copula framework (Torre et al., 2019b). The procedure results in an analytical formula for the stochastic emulator that can be used for computing marginals and correlations, as well as for generating new trajectories that resemble trajectories of the original stochastic simulator.

In our approach, the extension from discrete data to the continuous model (Challenge 2) is achieved by approximating the sampled trajectories by sparse regression-based PCE. A similar approach has been used by Navarro Jimenez et al. (2017) in the context of stochastic differential equations with the goal of sensitivity analysis, using non-intrusive pseudospectral projection to compute the PCE coefficients. The representation by sparse PCE can be seen as a variant of orthogonal series expansion (OSE) (Zhang and Ellingwood, 1994), which expands a second-order random process in terms of an orthogonal basis of the associated Hilbert space.

Note that when random fields are approximated based on a set of samples, it is most often assumed that the latter are collected on a discrete mesh in the index set, whereas this is not a requirement for our method. In such mesh-based approximations to random fields, PCE is often used for modeling the random variables arising in dimension-reduced expansions (Desceliers et al., 2006; Doostan et al., 2007; Das et al., 2009; Raisee et al., 2015; Abraham et al., 2018; Dai et al., 2019). This is distinct from our approach, as we use PCE to approximate the trajectories in the input space. Our approach yields an emulator for the whole input space (including unseen locations), while existing approaches are mostly focused on building an emulator on the discrete mesh where the samples were collected.

KLE represents a random field using an optimal orthogonal basis of the index space, resulting in an expansion in terms of deterministic functions weighted by random coefficients. These random coefficients are by construction uncorrelated, but unless the random field is a Gaussian random field, they are in general statistically dependent. Inferring the joint distribution of dependent random variables from samples is challenging but necessary for approximating a general non-Gaussian random field by KLE. To address this challenge of inference, several approaches have been proposed. Grigoriu (2010) suggests two methods to infer the joint distribution of the random coefficients of a series expansion model, of which one amounts to kernel density estimation, and the other to the fitting of a discrete joint distribution. Poirion and Zentner (2013, 2014) use KLE for modeling seismic ground motion time series, and model the random KLE coefficients by 1D sample CDFs assuming at most pairwise dependence (Poirion and Zentner, 2013), or by kernel density estimation (Poirion and Zentner, 2014). In the present paper, we investigate the use of kernel density estimation and inference of parametric joint distributions based on marginals and vine copulas.

This paper is organized as follows: in Section 6.2 we recall the relevant theory and definitions. In Section 6.3 we present our new stochastic emulator. The proposed method is then applied in Section 6.4, where we assess its performance on several examples of varying complexity. Here we
observe that the KLE is often significantly dominated by its first mode, a phenomenon that we investigate in Section 6.5. Finally, we draw conclusions and give an outlook on possible further developments in Section 6.6.

## 6.2 Theoretical foundation

We provide a brief summary of the relevant theory and concepts needed to construct our proposed stochastic emulator for stochastic simulators: random fields, polynomial chaos expansions, Karhunen-Loève expansion, and inference of joint probability distributions.

### 6.2.1 Stochastic simulators as random fields

Let $\mathbf{X}$ be a random vector with values in $\mathcal{D} \subset \mathbb{R}^d$, with finite variance and joint probability density function (pdf) $f_{\mathbf{X}}$. Denote by $\omega \subset \Omega$ an abstract random event in a probability space $(\Omega, \mathcal{F}, P)$. A stochastic simulator is a mapping

$$ \mathcal{M}: \mathcal{D} \times \Omega \to \mathbb{R}, \qquad (6.1) $$

$$ (x, \omega) \mapsto \mathcal{M}(x, \omega). \qquad (6.2) $$

Fixing $x \in \mathcal{D}$, the quantity $Y_x = \mathcal{M}(x, \cdot) : \Omega \to \mathbb{R}$ is a random variable. Fixing $\omega \in \Omega$, $\mathcal{M}(\cdot, \omega) : \mathcal{D} \to \mathbb{R}$ is a function in the input parameters, which we call trajectory or realization of the stochastic simulator (see also Figure 6.1). We assume that $Y_x$ has finite variance for all $x$, and that $\mathcal{M}(\cdot, \omega) \in L^2_{f_{\mathbf{X}}} (\mathcal{D})$ for all $\omega \in \Omega$.

<figure><img src="image_96.png" /><figcaption>**Figure 6.1:** Visual representation of a stochastic simulator when either the input parameters $x$ or the random event $\omega$ are held fixed, resulting in a random variable (left) or a deterministic function (right). The computational model is a high-rise building parametrized by several properties $x$ (visualized in the sketch by the shape of the building) subject to random earthquake events $\omega$ (visualized by 1D time series in different colors), whose output $\mathcal{M}(x, \omega)$ is a real number (e.g., the maximal displacement at the top floor).</figcaption></figure>

These definitions imply that a stochastic simulator $\mathcal{M}$ can be seen as a random field (also: stochastic process or random process) $\{Y_x\}_{x \in \mathcal{D}}$ with index set $\mathcal{D}$, i.e., as a family of random variables $\{Y_x\}$ indexed by $x \in \mathcal{D}$. In the following, we provide a brief reminder of a few random field basics. For more details, see e.g. Grigoriu (2002).
To fully characterize a general random field, one needs to specify the collection of all its finite-dimensional distributions

$$F_{Y_{x_1}, \dots, Y_{x_n}}(y_1, \dots, y_n) = \mathbb{P}(Y_{x_1} \le y_1 \wedge \dots \wedge Y_{x_n} \le y_n) \quad (6.3)$$

for all $n \ge 1$ and any $\boldsymbol{x}_1, \dots, \boldsymbol{x}_n \in \mathcal{D}$. Extending the concept of moments of random variables to random fields, the deterministic mean function of the random field is given by $\mu(\boldsymbol{x}) = \mathbb{E}[Y_{\boldsymbol{x}}]$. If $\mu(\boldsymbol{x}) = 0$, the random field is called centered. The (auto-)covariance function is defined by

$$c(\boldsymbol{x}, \boldsymbol{x}') = \mathbb{E}[(Y_{\boldsymbol{x}} - \mu(\boldsymbol{x}))(Y_{\boldsymbol{x}'} - \mu(\boldsymbol{x}'))]. \quad (6.4)$$

In general, a random field is not uniquely defined by its mean and covariance function. The only exception is the family of Gaussian processes, for which all finite-dimensional joint distributions are multivariate Gaussian distributions. For Gaussian processes, conditional distributions are again multivariate Gaussians, which lies at the foundation of the popular surrogate modeling technique Kriging/Gaussian process modeling. While Gaussian random fields are computationally convenient, random fields encountered in real-world problems (and in particular, stochastic simulators) are often non-Gaussian. One obvious argument is that Gaussian variables are unbounded while physical quantities are almost always bounded (Grigoriu, 2002).

A special feature of a stochastic simulator $\mathcal{M}$, as opposed to classical random fields, is that its index set is not an interval or a hypercube, but a general domain $\mathcal{D} \in \mathbb{R}^d$ with weight function $f_X$. We will use this property to build an accurate surrogate model for $\mathcal{M}$ respecting the probability density $f_X$ of the input space.

### 6.2.2 Polynomial chaos expansion

Polynomial chaos expansion (PCE) is a technique for modeling random variables using a basis of polynomials that are orthonormal w.r.t. a given probability density function (Ghanem and Spanos, 1991; Xiu and Karniadakis, 2002). In our algorithm (Section 6.3), we will use PCE to approximate trajectories $\mathcal{M}(\cdot, \omega)$ of the stochastic simulator, which can be seen as random variables $\mathcal{M}(\boldsymbol{X}, \omega)$ with their randomness induced by the uncertainty in the input $\boldsymbol{X}$.

Consider a random vector $\boldsymbol{X}$ with values in $\mathcal{D} \subset \mathbb{R}^d$ and independent components. Let $f_X(\boldsymbol{x}) = \prod_{i=1}^d f_{X_i}(x_i)$ be its probability density function (pdf) and assume that $\boldsymbol{X}$ has finite variance. Let $L_{f_X}^2(\mathcal{D})$ be the space of real-valued function that are square-integrable under $f_X$, i.e., $L_{f_X}^2(\mathcal{D}) = \{g : \mathcal{D} \to \mathbb{R} \mid \mathrm{Var}_X[g(\boldsymbol{X})] < +\infty\}$. Under certain assumptions on the random vector $\boldsymbol{X}$ (Xiu and Karniadakis, 2002; Ernst et al., 2012), there exists an orthonormal polynomial basis $\{\psi_\alpha \mid \alpha \in \mathbb{N}^d\}$ of $L_{f_X}^2(\mathcal{D})$, where each element is the product of univariate polynomials characterized by the multi-index $\boldsymbol{\alpha}$.

Let $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ be a (computational) model. Its output $Y = \mathcal{M}(\boldsymbol{X})$ is a random variable, which can be represented in terms of the orthonormal polynomial basis as $\mathcal{M}(\boldsymbol{X}) = \sum_{\alpha \in \mathbb{N}^d} a_\alpha \psi_\alpha(\boldsymbol{X})$ with $a_\alpha = \mathbb{E}_X[\mathcal{M}(\boldsymbol{X})\psi_\alpha(\boldsymbol{X})]$. This representation is called polynomial chaos expansion. In practice, a truncated expansion is computed,

$$\mathcal{M}(\boldsymbol{X}) \approx \mathcal{M}^{\mathrm{PCE}}(\boldsymbol{X}) = \sum_{\boldsymbol{\alpha} \in \mathcal{A}} a_{\boldsymbol{\alpha}} \psi_{\boldsymbol{\alpha}}(\boldsymbol{X}), \quad (6.5)$$
where $\mathcal{A} \subset \mathbb{N}^d$ is a finite subset of the full basis. The accuracy of a truncated PCE depends on three ingredients: the choice of $\mathcal{A}$, the method used for computing the coefficients $\boldsymbol{a} = (a_\alpha)_{\alpha \in \mathcal{A}}$, and the choice of points $\mathcal{X} \subset \mathcal{D}$ used in the coefficient computation method. An extensive overview of the state-of-the-art methods to determine these is given in Lüthen et al. (2021, 2022a).

### 6.2.3 Karhunen-Loève expansion

Karhunen-Loève expansion (KLE) is a well-established spectral expansion technique through which a random field is represented in terms of an optimal orthogonal basis for the index space, weighted by random coefficients (Karhunen, 1946; Loève, 1978). KLE transforms the random field, which is an uncountably infinite but correlated family of random variables $\{\mathcal{M}_x\}_{x \in \mathcal{D}}$, into a countably infinite but uncorrelated family of different random variables $\{\xi_i\}_{i=1,2,...}$. Furthermore, the random variables $\xi_i$ are typically of decreasing importance. KLE is therefore well suited and often used for discretization and modeling efforts for random fields.

To make these notions more precise, let $\{\mathcal{M}_x(\omega)\}_{x \in \mathcal{D}}$ be a random field. Denote by $\mu(\boldsymbol{x}) = \mathbb{E}[\mathcal{M}_x]$ its mean function, and by $c(\boldsymbol{x}, \boldsymbol{x}') = \text{Cov}[\mathcal{M}_x, \mathcal{M}_{x'}]$ its covariance function. Let $\mathcal{D}$ be closed and bounded. Let $c$ be continuous on $\mathcal{D} \times \mathcal{D}$ and assume that $\mathcal{M}_x$ has finite variance for all $\boldsymbol{x} \in \mathcal{D}$. Then the Karhunen-Loève expansion of the random field $\mathcal{M}_x$ is given by

$$ \mathcal{M}_x(\omega) = \mu(\boldsymbol{x}) + \sum_{k=1}^{\infty} \sqrt{\lambda_k} \xi_k(\omega) \phi_k(\boldsymbol{x}) \quad (6.6) $$

where $(\phi_k)_{k=1,2,...}$ is an orthonormal basis of $L^2(\mathcal{D})$, $\lambda_1 \ge \lambda_2 \ge \dots \ge 0$ is a non-increasing sequence of non-negative real numbers, and $\{\xi_k\}_{k=1,2,...}$ is a countable family of zero mean, unit variance, uncorrelated random variables.

Here, $(\lambda_k, \phi_k)$ are solutions to the integral eigenvalue problem

$$ \int_{\mathcal{D}} c(\boldsymbol{x}, \boldsymbol{x}') \phi_k(\boldsymbol{x}') d\boldsymbol{x}' = \lambda_k \phi_k(\boldsymbol{x}), \quad (6.7) $$

and $\xi_k$ is the result of the projection of $\mathcal{M}$ onto the spatial basis

$$ \xi_k(\omega) = \frac{1}{\sqrt{\lambda_k}} \int_{\mathcal{D}} \mathcal{M}(\boldsymbol{x}, \omega) \phi_k(\boldsymbol{x}) d\boldsymbol{x}. \quad (6.8) $$

From Eq. (6.6) and the properties of $\phi_k$ and $\xi_k$ it follows immediately that the covariance function can be expressed as

$$ c(\boldsymbol{x}, \boldsymbol{x}') = \sum_{k=1}^{\infty} \lambda_k \phi_k(\boldsymbol{x}) \phi_k(\boldsymbol{x}') \quad (6.9) $$

(Mercer's theorem). Note that the KLE random variables $\{\xi_k\}$ (herein KL-RV) do not enter this expression.

KLE is especially well-suited to Gaussian random fields, since in this case the random variables $\xi_k$ are standard Gaussian and independent. However, Eq. (6.6) holds for all random fields fulfilling the assumptions, not only for Gaussian random fields. The non-Gaussianity is modeled by the (possibly complex) joint distribution $f_\xi$ of the KL-RV.
Eqs. (6.6) to (6.8) are formulated in terms of $L^2(\mathcal{D})$, but they can be generalized: let $\mathbf{X}$ be a random variable with values in $\mathcal{D} \subset \mathbb{R}^d$, density $f_{\mathbf{X}}$, and finite variance. Then KLE can be generalized to the space $L^2_{f_{\mathbf{X}}}(\mathcal{D})$ instead of $L^2(\mathcal{D})$. In that case, the index set $\mathcal{D}$ does not have to be bounded, since the volume of $\mathcal{D}$ under measure $f_{\mathbf{X}}d\mathbf{x}$ is finite. This is called extended KLE (Iemma et al., 2006). This property is crucial for our proposed stochastic emulator, which we will introduce in Section 6.3.

In practice, the infinite expansion in Eq. (6.6) must be truncated. From the orthonormality of $\{\phi_k\}$ it follows from Eq. (6.9) that the variance of the random field is equal to $\sum_{k=1}^{\infty} \lambda_k$. The sequence $\lambda_1 \ge \lambda_2 \ge \dots \ge 0$ is non-increasing, and typically (depending on the correlation length of the random field) this sequence decays rather quickly to zero. Loosely speaking, the higher the correlation between different locations in the index set, the fewer spatial basis functions are needed to approximate the trajectories, and the faster the decay of the eigenvalues. Knowing this, the KLE can be truncated at the $M$-th term with $M$ chosen so that the fraction of explained variance is sufficiently large:

$$ \frac{\sum_{k=1}^{M} \lambda_k}{\sum_{k=1}^{\infty} \lambda_k} > 1 - \epsilon \qquad (6.10) $$

for a small threshold parameter $\epsilon > 0$ (e.g., $\epsilon = 0.001$).

KLE is closely related to function principal component analysis (fPCA) (Besse and Ramsay, 1986; Ramsay and Silverman, 2005). To compute a solution to the integral eigenvalue problem in Eq. (6.7), there are several possibilities (Ramsay and Silverman, 2005, Section 8.4): the integrals can be approximated numerically; the eigenproblem can be discretized on a number of representative grid points in $\mathcal{D}$ (this is the approach chosen by the majority of modelers, including Azzi et al. (2019)); or the eigenproblem can be written in terms of a suitable (truncated) spatial basis, which transforms the problem into a (finite-dimensional) discrete eigenvalue problem. The third approach is related to orthogonal series expansion (OSE) (Zhang and Ellingwood, 1994). It is used by Poirion and Zentner (2014), who derive the explicit discrete problem for a basis consisting of interpolation functions, building on results by Besse and Ramsay (1986) and Besse (1991). We use this approach together with the orthogonal basis provided by polynomial chaos expansion (Section 6.3). Detailed calculations are provided in Section 6.A.

### 6.2.4 Inference of the joint distribution of the Karhunen-Loève random variables

Characterizing the dependent (but uncorrelated) Karhunen-Loève random variables (KL-RV) $\xi_k, k = 1, \dots, M$ correctly is important for the accurate modeling of a general non-Gaussian stochastic process (Grigoriu, 2010). However, inferring the joint distribution of a random vector is a challenging task. The main challenge is the scarcity of data: the higher the dimensionality, the more samples are needed to be able to correctly infer the dependence structure of the data. We need to construct a suitable parametric or non-parametric model to accurately describe the joint distribution. In the following, we introduce the marginal-copula framework, which is a powerful tool to represent and infer complex dependence structures between random variables (Nelsen, 2006; Torre et al., 2019).
Let $\mathcal{Z}$ be any $M$-dimensional random vector with multivariate cumulative distribution function (CDF) $F_{\mathcal{Z}}$ and marginal distributions $F_{Z_i}$. The so-called Sklar's theorem states that $F_{\mathcal{Z}}$ can be written as

$$F_{\mathcal{Z}}(z_1, \dots, z_M) = C(F_{Z_1}(z_1), \dots, F_{Z_M}(z_M)), \quad (6.11)$$

where the function $C: [0, 1]^d \to \mathbb{R}$ is called copula (Sklar, 1959; Nelsen, 2006). $C$ is a CDF with uniform marginals, which defines the dependence structure of the random vector $\mathcal{Z}$. $C$ is unique if all marginals $F_{Z_i}$ are continuous, and it holds that

$$C(u_1, \dots, u_M) = F_{\mathcal{Z}}(F_{Z_1}^{-1}(u_1), \dots, F_{Z_M}^{-1}(u_M)). \quad (6.12)$$

Let an i.i.d. sample $\mathcal{Z} = \{z^{(1)}, \dots, z^{(N)}\}$ of the random vector $\mathcal{Z}$ be given. The goal is to infer the joint distribution $F_{\mathcal{Z}}$ from this sample. For this, the copula representation of Eq. (6.11) is convenient, since it allows inferring the marginals and the dependence structure of the data separately, as briefly explained in the following.

To infer the marginal distributions, we consider two options. The first is parametric inference, where we choose from a set of parametric probability distributions with zero mean and unit standard deviation (see Table 6.1) the distribution with the smallest Akaike information criterion (AIC). If a distribution family has more than two parameters, its remaining parameters are chosen by maximum likelihood. We utilize the uncertainty quantification software UQLab (Marelli and Sudret, 2014; Torre et al., 2021) with a modification prescribing the desired moments.

<table><caption>**Table 6.1:** Considered marginal families with zero mean and unit standard deviation. The last column lists the remaining degrees of freedom $k$ after fixing the first two moments. The Akaike information criterion is then given as $AIC = 2k - 2 \log \mathcal{L}$, where $\mathcal{L}$ is the likelihood.</caption><thead><tr><th>Type</th><th>Parameter</th><th>k</th></tr></thead><tbody><tr><td>Uniform U([a, b])</td><td>a = -√3, b = √3</td><td>0</td></tr><tr><td>Gaussian N(μ, σ)</td><td>μ = 0, σ = 1</td><td>0</td></tr><tr><td>Gumbel (for maxima) G(μ, β)</td><td>μ ≈ -0.4501, β ≈ 0.7797</td><td>0</td></tr><tr><td>Gumbel (for minima) G<sub>min</sub>(μ, β)</td><td>μ ≈ 0.4501, β ≈ 0.7797</td><td>0</td></tr><tr><td>Logistic P(μ, s)</td><td>μ = 0, s ≈ 0.5513</td><td>0</td></tr><tr><td>Laplace L(μ, b)</td><td>μ = 1, b = 1/√2</td><td>0</td></tr><tr><td>Beta B(a, b, r, s)</td><td>a, b chosen according to data bounds<br>r = <sup>a(ab+1)</sup>⁄<sub>b−a</sub>, s = <sup>b(ab+1)</sup>⁄<sub>a−b</sub></td><td>2</td></tr></tbody></table>

A second popular method to represent marginal behavior non-parametrically is kernel density estimation (KDE) (Wand and Jones, 1995; Simonoff, 1996), which has also been proposed for estimating the distribution of KL-RV (Grigoriu, 2010; Poirion and Zentner, 2014). Here the distribution is modeled as a Gaussian mixture, where the Gaussian density functions are centered in the data points and share the same standard deviation, called bandwidth in the case of 1D KDE. We adopt a bandwidth estimation method optimal for data with Gaussian distribution (Bowman and Azzalini, 1997).

To characterize the dependence structure, we use a copula. While any multivariate CDF with uniform marginals $U([0, 1])$ constitutes a copula, there are a number of well-known parametric families (see, e.g., Nelsen (2006); Joe (2014)). Besides the independence copula and the families
derived from multivariate elliptical distributions, most of these parametric families are pair copulas, i.e., they couple only two variables. Constructing meaningful parametric copulas for more than two variables (other than elliptical copulas) is in general difficult (Nelsen, 2006).

A solution is to decompose the M-variate copula into a product of conditional pair copulas, which is known as vine copula construction (Bedford and Cooke, 2002). This is always possible as a consequence of the chain rule of probability. In general, a vine copula is the product of $\frac{M(M-1)}{2}$ pair copulas.<fnref n="3" /> The factorization into pair copulas is not unique but depends on the ordering and grouping of variables. Two classes of vine copulas, differing in the order in which the variables are grouped into pairs, are the drawable vine (D-vine) (Kurowicka and Cooke, 2005) and the canonical vine (C-vine) (Aas et al., 2009). For a more detailed description of the vine copula construction, we refer to Aas et al. (2009) and Torre et al. (2019b).

To infer a copula from data, we first map the multivariate data to $[0, 1]^d$ by applying element-wise the inferred marginal CDFs (see Eq. (6.11)). Then we infer the dependence structure by using Kendall's tau to determine the groupings of variables as well as their order in the vine copula (Aas et al., 2009; Torre et al., 2019b). For each pair copula, the parameters are identified by maximum likelihood. Finally, the best-fitting copula is chosen using AIC. This approach is implemented in the statistical inference module of UQLab (Torre et al., 2021). The list of available copula families can be found in Lataniotis et al. (2021, Section 1.4).

## 6.3 Surrogating a stochastic simulator from a set of samples

We are now ready to describe the construction of our spectral surrogate model for a stochastic simulator. Assume that discrete samples of the stochastic simulator $\mathcal{M}$ are available in the following form:

$$ \mathcal{T}_r = \left\{ (\mathbf{x}^{(r,i)}, \mathcal{M}(\mathbf{x}^{(r,i)}, \omega^{(r)})) : i = 1, \dots, N_r \right\}, \quad r = 1, \dots, R \quad (6.13) $$

i.e., in the form of discrete evaluations of the stochastic simulator on $R$ trajectories, where for every $r$, $\{\mathbf{x}^{(r,i)} : i = 1, \dots, N_r\}$ is an i.i.d. sample from the input distribution $f_X$, the so-called experimental design. In particular, for different trajectories the samples can be taken at different locations, i.e., for $r_1 \neq r_2$ we can have $\mathbf{x}^{(r_1,i)} \neq \mathbf{x}^{(r_2,i)}$ and in principle even different numbers of samples $N_{r_1} \neq N_{r_2}$. However, here we assume for notational simplicity that $N_r = N$ for all $r$.

Our proposed method consists of the following steps (see also Figure 6.1):

1. **Approximate each discrete trajectory $\mathcal{T}_r$ by a sparse PCE $\mathcal{M}_r^{\text{PCE}}$ in $L_{f_X}^2(\mathcal{D})$:**

$$ \mathcal{M}_r^{\text{PCE}}(\mathbf{x}) = \sum_{\alpha \in \mathcal{A}^{(r)}} a_{\alpha}^{(r)} \psi_{\alpha}(\mathbf{x}) \quad (6.14) $$

with $\mathcal{A}^{(r)}$ the set of regressors with nonzero associated coefficient $a_{\alpha}^{(r)}$. We use a total-degree basis with degree- and $q$-norm adaptivity to determine the truncation set $\mathcal{A}^{(r)}$ (Blatman and Sudret, 2011; Lüthen et al., 2022a) and apply the least-angle regression

<footnote n="3">³There are $M-1$ unconditional pair copulas; $M-2$ pair copulas conditioned on 1 other variable; $M-3$ conditioned on 2 other variables; and so on, until there is 1 pair copula conditioned on all except 2 variables.</footnote>
solver to compute the coefficients (sparse PCE) (Blatman and Sudret, 2011; Lüthen et al., 2021).

2. Determine a set $\mathcal{A}$ of regressors that jointly represents all trajectories well:

* Identify the union $\mathcal{A} = \bigcup_{r=1}^{R} \mathcal{A}^{(r)}$ of all chosen regressors.

* To keep the size of the basis manageable, discard the regressors with the smallest sum of squares of coefficients over all trajectories $(\sum_{r=1}^{R} (a_{\alpha}^{(r)})^2)$ until $P = |\mathcal{A}| \le \frac{N}{2}$ regressors or less are left in $\mathcal{A}$.

* To avoid discontinuous behavior resulting from sparse selection, recompute the coefficients of every trajectory by ordinary least squares (OLS), using the chosen set of regressors $\mathcal{A}$.

This results in $R$ PCE trajectories, where each trajectory uses the same set of $P$ PCE basis functions.

3. Center the PCE trajectories by subtracting the sample mean

$$ \hat{\mu}^{\text{PCE}}(\mathbf{x}) = \frac{1}{R} \sum_{r=1}^{R} \mathcal{M}_{r}^{\text{PCE}}(\mathbf{x}) = \sum_{\alpha \in \mathcal{A}} \left( \frac{1}{R} \sum_{r=1}^{R} a_{\alpha}^{(r)} \right) \psi_{\alpha}(\mathbf{x}) \quad (6.15) $$

which is itself a PCE. We denote by $\tilde{\mathcal{M}}_{r}^{\text{PCE}}(\mathbf{x}) = \mathcal{M}_{r}^{\text{PCE}}(\mathbf{x}) - \hat{\mu}^{\text{PCE}}(\mathbf{x})$ the centered PCE trajectories. Extract the coefficients $\tilde{a}_{\alpha}^{(r)}$ of the centered trajectories and store them in a $P \times R$ matrix $\tilde{a}$.

4. Apply extended KLE to the set of PCE trajectories. The sample covariance function has the form

$$ \hat{c}(\mathbf{x}, \mathbf{x}') = \frac{1}{R-1} \sum_{r=1}^{R} \tilde{\mathcal{M}}_{r}^{\text{PCE}}(\mathbf{x}) \tilde{\mathcal{M}}_{r}^{\text{PCE}}(\mathbf{x}'). \quad (6.16) $$

Computing the eigenfunctions $\phi(\mathbf{x})$ of the associated integral eigenvalue problem in Eq. (6.7) is equivalent to computing a PCA on the PCE coefficients, i.e., equivalent to solving the following $P$-dimensional eigenproblem for $\tilde{a}$:

$$ \Sigma b = \lambda b, \quad (6.17) $$

where $\Sigma = \frac{1}{R-1} \tilde{a} \tilde{a}^T$ (see Section 6.A for the derivation of this equivalence). The eigenvectors $b$ contain the coefficients of the eigenfunctions represented in the PCE basis: $\phi(\mathbf{x}) = \sum_{\alpha \in \mathcal{A}} b_{\alpha} \psi_{\alpha}(\mathbf{x})$.

5. Identify the truncation order $K \ll P$ for the KLE based on a given threshold for the explained variance. We use a threshold of 99.9% (see Eq. (6.10)).

6. Compute the realizations of the KL-RV $\xi_i$ from the sample trajectories by projecting onto the eigenfunctions. Due to the orthonormality of the PCE basis, this can be done analytically (see Section 6.A.2). Denote the realizations by $\xi^{(r)} \in \mathbb{R}^K$.

7. Infer the joint distribution $f_{\xi}$ of random KL coefficients from the data set $\{\xi^{(r)}\}_{r=1,...,R}$. We will test four methods consisting of the techniques described in Section 6.2.4:

(a) Option 1: assume standard Gaussian marginals, which implies independence;

(b) Option 2: parametric inference of the marginals (with moment constraints) and of the copula;

(c) Option 3: 1D kernel density estimation of each marginal, assuming independence;

(d) Option 4: 1D kernel density estimation of each marginal and parametric inference of the copula.
The resulting stochastic model for the random field $\mathcal{M}$ is

$$ \hat{\mathcal{M}}(\mathbf{x}, \cdot) = \hat{\mu}(\mathbf{x}) + \sum_{k=1}^{K} \sqrt{\lambda_k} Z_k(\cdot) \underbrace{\left( \sum_{\alpha \in \mathcal{A}} b_{\alpha}^{(k)} \psi_{\alpha}(\mathbf{x}) \right)}_{=\phi_k(\mathbf{x})} \quad (6.18) $$

where $\mathcal{Z} = (Z_1, Z_2, \dots, Z_K)$ is a random vector distributed according to the inferred joint distribution $f_{\xi}$.

The full procedure is visualized in Figure 6.1.

Using the stochastic emulator constructed in Eq. (6.18), we can easily compute the following quantities.

* The mean function $\hat{\mu}$ is given by the sample mean of the approximated trajectories (PCE trajectories), see Eq. (6.15).

* The covariance function $\hat{c}(\cdot, \cdot)$ can be computed from the KLE eigenfunctions using the truncated version of Eq. (6.9). Note that this relation does not involve the KL-RV.

* New trajectories (i.e., realizations of the random field) can be generated by drawing new samples of the KL-RV $\xi_k$, and evaluating Eq. (6.6).

* A histogram of the marginal pdf $f_{M_{x'}}$ of the random field at any input space location $x'$ can be created by generating many new trajectories and evaluating them at $x'$.

**Remark 1 (Another stochastic emulator).** A simple stochastic emulator able to model marginal distributions $f_{M_{x'}}$ can be constructed by evaluating all PCE trajectories from Step 2 above at the new location $x'$ and computing a kernel density estimate on the resulting set of predictions. This method will be used as a comparison method for marginal estimation in Section 6.4. However, unlike our stochastic emulator in Eq. (6.18), this simple emulator is not able to resample trajectories.

**Remark 2 (Alternatives to PCE).** We choose PCE to approximate the sampled trajectories because it is a powerful method for deterministic surrogate modeling. However, the choice of PCE in the above method is not crucial: without any changes to the methodology, PCE could be replaced by any other spectral expansion onto an orthonormal basis of $L_{fX}^2(\mathcal{D})$, e.g., a Poincaré basis (Lüthen et al., 2022c) or a spline basis (Rahman, 2020). From the orthonormality of the basis it follows that functional PCA in $L_{fX}^2(\mathcal{D})$ becomes traditional (unweighted) PCA in the coefficient space (see Section 6.A), which avoids the expensive numerical solution of the integral eigenvalue problem in $d$ dimensions, and instead solves an inexpensive discrete eigenvalue problem.
<figure><img src="image_97.png" /><figcaption>**Figure 6.1:** Sketch of our stochastic emulator, starting with stochastic simulator samples (discrete trajectories) at the top and resulting in the stochastic emulator at the bottom, which is a KLE that includes a probabilistic model of the KL-RV. The sketch is purely for illustration and does not display real data. Note that there are two equivalent ways to arrive at the third box: through extended KLE and through PCA on the coefficients.</figcaption></figure>
## 6.4 Numerical experiments

To analyse the performance of our stochastic emulator, we apply it to three models of increasing complexity: the three-dimensional Ishigami function with two random parameters (Section 6.4.1), the borehole model with five hidden (latent) variables (Section 6.4.2), and finally the Heston stochastic volatility model, a system of two stochastic ODEs with six inputs that has already been used by Zhu and Sudret (2021c) as a stochastic emulator benchmark model (Section 6.4.3).

We first investigate the pointwise approximation capabilities of our emulator by plotting the stochastic simulator and emulator responses at selected points throughout the input domain. Then, we investigate the convergence behavior of our stochastic emulator using the following global error measures:

* The global convergence of the marginal distributions is assessed using the *averaged normalized Wasserstein distance*. The *Wasserstein distance of order two* between two random variables $Y_1, Y_2$ with quantile functions (inverse CDF) $Q_1, Q_2$ is defined by (Villani, 2009)

$$d_{\text{WS}}(Y_1, Y_2) = \|Q_1 - Q_2\|_2 = \sqrt{\int_0^1 (Q_1(u) - Q_2(u))^2 \, du}. \quad (6.19)$$

To measure the global quality of marginal approximation, we consider the quantity

$$\epsilon_{\text{marg}} = \mathbb{E}_{\mathbf{X}} \left[ \frac{d_{\text{WS}}(\mathcal{M}(\mathbf{X}, \cdot), \hat{\mathcal{M}}(\mathbf{X}, \cdot))}{\sigma(\mathcal{M}(\mathbf{X}, \cdot))} \right], \quad (6.20)$$

computed by Monte-Carlo integration on a validation set with $N_{\text{val}} = 1,000$ points and $R_{\text{val}} = 10,000$ replications (Zhu and Sudret, 2021b).

* The global error between the true covariance function $c$ and the emulated one $\hat{c}$ is computed by

$$\epsilon_{\text{cov}} = \|c - \hat{c}\|_{L^2_{f_X}(\mathcal{D}) \times L^2_{f_X}(\mathcal{D})} \approx \frac{1}{N_{\text{val}}} \|C - \hat{C}\|_F \quad (6.21)$$

where $C$ and $\hat{C}$ denote the true and emulated covariance matrices for a validation sample $\{\mathbf{x}^{(i)} : i = 1, \dots, N_{\text{val}}\}$, and $\|\cdot\|_F$ is the Frobenius norm.

### 6.4.1 Stochastic Ishigami function

#### 6.4.1.1 Problem statement

The Ishigami function is a well-known benchmark function for deterministic surrogate models. It is highly non-linear and has significant interaction terms. It becomes a stochastic simulator by treating its parameters $a$ and $b$, which are usually fixed at $a = 0.7$ and $b = 0.1$, as additional random variables:

$$f(\mathbf{X}; A, B) = \sin(X_1) + A \sin(X_2)^2 + B X_3^4 \sin(X_1). \quad (6.22)$$

$A$ and $B$ have here the role of so-called *hidden* or *latent* random variables. In other words, we assume that they cannot be observed, and that therefore their values cannot be utilized in
the surrogate modeling process. They introduce stochasticity into the otherwise deterministic Ishigami model. Here we model *A* and *B* as lognormal random variables with mean 7 and standard deviation 0.7, and mean 0.1 and standard deviation 0.1, respectively. We assume that both variables are coupled with a Clayton pair copula with parameter 1.5. The non-hidden (explicit) input variables are as usual **X** = (*X*₁, *X*₂, *X*₃), which are independent and uniformly distributed in [-π, π]. A Sobol' analysis of *f*(**X**; *A*, *B*) in Eq. (6.22) reveals that the main effect of the group of explicit input parameters (*M*₁ in Eq. (6.28)) is approx. 75%, while the interaction effect between the explicit and the latent group is approx. 25%, and the main effect of the group of latent variables is negligible. Samples of the input space and the latent space are displayed in Figure 6.1.

<figure><img src="image_98.png" /><figcaption>**Figure 6.1:** Samples of the input space (left) and the latent space (right). The red line in Figure 6.1a is the trajectory along which the simulator/emulator response is plotted in Figure 6.7. The red dots annotated by small numbers denote the five points that are used for visualization in the following. Figure 6.1b shows a sample of the latent space (*A*, *B* in Eq. (6.22)).</figcaption></figure>

We use different experimental design sizes $N \in \{50, 100, 150\}$ and a maximum degree of $p = 14$ for the PCE trajectories (with degree-adaptivity (Blatman and Sudret, 2011)). This results in a relative mean-squared error in the order of $10^{-3}/10^{-5}/10^{-10}$, respectively. We also test different numbers of trajectories $R \in \{10, 30, 100, 300\}$, and use a different experimental design for each trajectory. For each combination of experimental design size and number of trajectories, we conduct 50 independent repetitions. All resulting stochastic emulators are evaluated on the same validation set, consisting of $R_{\text{val}} = 10,000$ trajectories of the true stochastic simulator, each evaluated on a set of $N_{\text{val}} = 1,000$ points in the input space.
### 6.4.1.2 Analysis of the KL-RV samples

To illustrate the type of result obtained with our proposed stochastic emulator described in Section 6.3, we now present scatterplots showing realizations of the following random quantities: 1) the KL-RV (compressed representation of PCE coefficients) resulting from step 7; 2) the PCE coefficients resulting from transforming the KL-RV samples to the PCE coefficient space. Detailed results for the prediction $M(x', \cdot)$ at a new location $x'$ for a number of new trajectories are presented in Section 6.4.1.3 below. The results shown here are based on parametric inference of marginals and copula (Option (b) of Step 7).

The truncation of the KLE (Step 5 of our algorithm) typically results in two modes with eigenvalues $\lambda_1 \in [3, 5]$ and $\lambda_2 \approx 0.1$ depending on the size and realization of the experimental design. We display a specific example in Figure 6.2, which is computed from 100 trajectories with 150 samples each, and has eigenvalues $\lambda_1 = 4.46$ and $\lambda_2 = 0.10$. The figure shows resampled values for the KL-RV: in black, samples computed from validation trajectories by projecting first onto the truncated PCE space and then onto the eigenfunctions; in red, new samples drawn from the input object inferred in Step 7 of our algorithm (Section 6.3). We see that their inferred joint distribution (Beta and Gumbel marginals, with a Clayton copula) visually matches the validation data well.

<figure><img src="image_99.png" /><figcaption>**Figure 6.2:** KL-RV coefficient samples computed from validation trajectories (black) and new samples from the stochastic emulator (red). This is data from one experiment with $N = 150$ and $R = 100$, max degree $p = 14$. Number of validation trajectories and resampled PCE coefficients: 200 each. Number of KL modes: $M = 2$. Inferred distribution of KL-RV: Beta and Gumbel, with Clayton copula (parameter 0.32). The corresponding eigenvalues are $\lambda_1 = 4.46$ and $\lambda_2 = 0.10$.</figcaption></figure>

The KL-RV are the compressed representation of the random PCE coefficients (which in turn encode trajectories). Mapping the realizations of the KL-RV back to the PCE coefficient space, we obtain the samples displayed in Figure 6.3 for $N = 150$ and $R = 100$. Validation samples from
the original stochastic simulator (generated by regressing them onto the truncated PCE basis) are displayed in black, while 200 resampled PCE coefficient vectors generated from the stochastic emulator are shown in red. We only show the 5 coefficients with maximal mean absolute value. We see that the validation samples have a slightly larger spread than the emulator samples, but that overall the behavior is matched well. Some parameters have linear functional dependence, e.g., $a_1, a_8$ and $a_{17}$, which is perfectly reproduced by the emulator. These parameters correspond to the basis functions $\alpha_1 = (0,0,0)$ (constant term), $\alpha_8 = (0,4,0)$ and $\alpha_{17} = (0,6,0)$ and are needed to emulate the second term of the stochastic Ishigami model in Eq. (6.22). There are no interactions with the other terms, therefore a different value of $A$ just proportionally changes the relative weighting of these terms. A similar explanation holds for $a_2$ and $a_6$ with $\alpha_2 = (1,0,0)$ and $\alpha_6 = (1,0,2)$, which are involved in emulating the first and the third term of Eq. (6.22) and change proportionally with $B$.

<figure><img src="image_100.png" /><figcaption>**Figure 6.3:** PCE coefficient samples computed from validation trajectories (black) and new samples from the stochastic emulator (red). This is data from one experiment with $N = 120$ and $R = 100$, max degree $p = 14$. Number of validation trajectories and resampled PCE coefficients: $R_{val} = 200$ each. The PCE coefficients are sorted by mean magnitude, and we only display the largest 5 out of total 75 nonzero coefficients.</figcaption></figure>

### 6.4.1.3 Marginal performance on selected validation points

We now investigate the performance of the stochastic emulator with parametric inference of KL-RV marginals and copulas (choice 7b) on a selection of out-of-sample validation points, i.e., points that were not used for training.
Figure 6.4 shows the histograms and pairwise scatterplots of samples from the output random variables $Y_i = \mathcal{M}(\boldsymbol{x}^{(i)}, \cdot)$ and $\hat{Y}_i = \hat{\mathcal{M}}(\boldsymbol{x}^{(i)}, \cdot)$ of the stochastic simulator and emulator, respectively. The five selected validation locations $\{\boldsymbol{x}^{(i)}\}_{i=1}^5$ in the input space are visualized in Figure 6.1a by red dots. Each black (resp. red) point in Figure 6.4 is a new trajectory of the stochastic simulator (resp. emulator) evaluated at the five given points. Both samples have the same size (200 new trajectories). Overall, the model behavior is captured well, but the stochastic simulator has a slightly larger spread (see e.g. $Y_2$ vs. $Y_3$).

<figure><img src="image_101.png" /><figcaption>**Figure 6.4:** Scatterplot of output $Y_i = \mathcal{M}(\boldsymbol{x}^{(i)}, \cdot)$ of the stochastic simulator (black) and of output $\hat{Y}_i = \hat{\mathcal{M}}(\boldsymbol{x}^{(i)}, \cdot)$ of the parametric stochastic emulator (red, created from training set with $N = 150$, $R = 100$) for five validation points sampled from the input space. The location of these five points is illustrated in Figure 6.1a.</figcaption></figure>

From the data in the off-diagonal scatter plots in Figure 6.4, we can compute the sample covariance matrix. However, we can also compute the covariance analytically from the KLE eigenfunctions, using Eq. (6.9). In Figure 6.5, we use the five illustrative points shown in Figure 6.1a to compare this covariance estimate to a validation covariance matrix computed empirically from 10,000 trajectories of the stochastic simulator. Qualitatively, the covariance is reproduced well, although the KLE-based covariance is slightly smaller in magnitude than the empirical covariance.
<figure><img src="image_102.png" /><figcaption>**Figure 6.5:** Covariance approximation for the Ishigami model for the $N_{\text{val}} = 5$ validation locations in the input space illustrated in Figure 6.1a. Computation based on $N = 150, R = 100$, max degree $p = 14$. KLE-based covariance: computed from eigenfunctions as in Eq. (6.9). Empirical covariance: based on the validation set comprising 10,000 trajectories of the stochastic simulator.</figcaption></figure>

In Figure 6.6, we visualize the marginal distribution $f_{Y_1}$ of $Y_1 = M(\boldsymbol{x}^{(1)}, \cdot)$ at one validation point (the point marked with “1” in Figure 6.1a) for an increasing number of trajectories in the training set, and 4 independent repetitions of each experiment. The estimates for the marginal distribution $f_{Y_1}$ are computed by KDE from 10,000 samples from the constructed stochastic emulator, while the histogram and the dashed curve represents a validation set of 10,000 samples of the original stochastic simulator. As expected, we observe that with an increasing number $R$ of trajectories, the shape of the predicted marginal becomes closer to the kernel density estimate of the validation set and shows less variation.

<figure><img src="image_103.png" /><figcaption>**Figure 6.6:** Prediction of Y-marginal at validation point $\boldsymbol{x}_{\text{val}}^{(1)} = (-\frac{\pi}{2}, \frac{\pi}{2}, -\frac{\pi}{2})$ for the Ishigami model and $R = \{10, 30, 100, 300\}$ trajectories with the parametric stochastic emulator for 4 independent repetitions. Visualization of predicted marginals by KDE using 10,000 samples. Number of experimental design points $N = 150$, max degree $p = 14$. The approximation error is in the order $O(10^{-10})$.</figcaption></figure>

Finally, to assess visually how well the resampled trajectories match the behavior of the original stochastic model, we plot in Figure 6.7 a 1D slice of 10 new trajectories generated by the stochastic simulator (left) and the stochastic emulator (middle). The slice through the input space is shown in Figure 6.1a by a red line. On the right, data for the same slice is shown, but this time we show quantiles aggregated over 10,000 new trajectories each. The trajectory slices look qualitatively similar, although there is a lot of variability between individual realizations. From the aggregated data on the left, we see that the bulk of the distribution (10%-90% quantile) is predicted quite accurately. Interestingly, in Figure 6.7c it seems that the trajectories of the
stochastic emulator (red) have a larger spread than the ones of the simulator (black), contrary to the results earlier in this section, which always showed the simulator having a larger spread than the emulator. This illustrates the difficulty of inferring global behavior from local observations. Theoretically, the emulator should have a smaller variance than the simulator, because terms are missing from Eq. (6.9) due to truncation.

<figure><img src="image_104.png" /><figcaption>**Figure 6.7:** Visualization of stochastic simulator/emulator response Y when following a 1D slice through the input space, which is illustrated with a red line in Figure 6.1a. N = 150, R = 100. The left and middle plots show 10 trajectories each. The right plot aggregates the values from 10,000 trajectories to show quantile information.</figcaption></figure>

### 6.4.1.4 Convergence with the number of trajectories

To assess the global performance of our proposed method, we now construct stochastic emulators for all combinations of input space experimental design sizes $N \in \{50, 100, 150\}$ and numbers of trajectories in the range $R \in \{10, 30, 100, 300\}$. We then evaluate each of the resulting stochastic emulators $R_{\text{val}} = 10,000$ times at $N_{\text{val}} = 1,000$ validation points in the input space (out-of-sample, i.e. not used for training) and compute the errors as described in Section 6.4. Each combination is independently repeated 50 times to account for the statistical uncertainty of the sampling of both experimental design and trajectories, which allows us to display results in the form of Tukey boxplots.

In Figure 6.8a, we display the global convergence of marginal predictions for the parametric stochastic emulator in terms of $\epsilon_{\text{marg}}$ defined in Eq. (6.20). Each boxplot represents one experiment (i.e., a specific number of experimental design points and number of trajectories), repeated independently 50 times. The value of the averaged and normalized Wasserstein distance $\epsilon_{\text{marg}}$ is, by itself, difficult to interpret. To aid the interpretation and give an idea of the quality of the approximation, we add two auxiliary quantities to the plot:

* The averaged and normalized Wasserstein distance is computed based on samples. As a lower bound, we independently sample $100 \times 2$ validation sets (each consisting of $R_{\text{val}} = 10,000$ trajectories evaluated at $N_{\text{val}} = 1,000$ points in the input domain; each pair of validation sets shares the same points in the input domain). We then compute the error in Eq. (6.20) for each of the 100 pairs. The median and quantiles (0.25–0.75 and 0.05–0.95)
of this value are displayed in Figure 6.8a in gray, indicating the best possible error that
can be achieved due to the natural variability of the sample estimates.

* A priori, it is unclear which value of the (averaged) normalized Wasserstein distance corresponds to predicted marginals that are visually close to the true marginals. To have some concrete examples on what a specific value of the normalized Wasserstein distance means, we consider the marginals predicted at one chosen validation point, shown in Figure 6.8b. We add the corresponding value of the normalized Wasserstein distance between simulator and emulator prediction as a small colored circle to the plot in Figure 6.8a.

We observe that the quality of the marginal estimates improves as we increase the size of the input parameter sample, which is expected since the PCE approximation of the trajectories becomes better with increasing experimental design size. $N = 50$ points are clearly too few to achieve a good estimate, whereas $N = 100$ and $N = 150$ show convergence with the number of trajectories, indicating that the error from statistical inference is the dominating one. While the convergence of the marginals with the number of trajectories looks slow, the improvement of the marginal shapes is actually significant (compare the values with Figure 6.8b).

<figure><img src="image_105.png" /><figcaption>**Figure 6.8:** Convergence of $\epsilon_{\text{marg}}$ and $\epsilon_{\text{cov}}$ (Eqs. (6.20) and (6.21)) for increasing number of available trajectories and parameter locations. Results for the stochastic emulator with parametric inference (choice 7b of our algorithm in Section 6.3) and 50 replications. The errors are computed based on a validation set of size $N_{\text{val}} = 1,000$, $R_{\text{val}} = 10,000$. The gray areas and the dashed line represent quantiles and the median of a lower bound estimate for the respective error measure computed from $100 \times 2$ independent MC samples of size $R_{\text{val}} = 10,000$ generated by the true stochastic simulator. The colored points in Figure 6.8a correspond to the results for a single replication and validation point as shown in Figure 6.8b and help assess the meaning of the numerical error measures.</figcaption></figure>
In addition to marginal predictions, our stochastic emulator can also emulate the covariance function, using Eq. (6.9). Since this equation relies only on the KL eigenfunctions, not on the KL-RV, the choice of inference method in Step 7 of our algorithm in Section 6.3 does not affect these results. In Figure 6.8c we display the convergence of $\epsilon_{cov}$ from Eq. (6.21). We observe that the error decreases with increasing numbers of trajectories. For the largest numbers of trajectories and experimental design points, the error is already in the range of the rough lower bound on achievable accuracy (obtained as described above by empirical sampling of the true model). Interestingly, unlike the marginal error in Figure 6.8a, an increasing number $N$ of input parameter samples does not lead to a smaller covariance error. This indicates that the covariance estimate is less sensitive to the quality of the trajectory approximation, while the inference of the distribution of the KL-RV is more sensitive to it.

So far, we showed results for the stochastic emulator with parametric inference only (Option 7b). Now, in Figure 6.9 we compare the four inference options described in Step 7 of the algorithm in Section 6.3 with the results of a fifth method described in Remark 1, which we call here PCE-KDE. We use the experimental design size $N = 100$, which yields PCE approximations with relative validation error of $10^{-5}$. Due to this close fit, the PCE-KDE estimate can be considered as near-optimal estimate given the available data.

The error $\epsilon_{marg}$ is again computed on a validation set consisting of $R_{val} = 10,000$ trajectories evaluated at $N_{val} = 1,000$ points in the input space.

We observe that for 20 trajectories, the five methods show almost identical performance. This suggests that 20 trajectories are not yet enough to infer a distribution that is able to generalize to unseen data, so that any marginal distribution with mean zero and unit standard deviation provides a reasonable approximation. Comparing with the results for PCE-KDE, we see that our emulator is similarly accurate in prediction at an unseen point as a kernel density estimate using the training set of highly accurate PCE trajectories. This suggests that we do not lose much accuracy by applying our KLE approach on top of the PCE approximation, which can be seen as a form of dimension reduction in the stochastic space.

For the larger number of trajectories, $R = 100$ and $R = 500$, we do observe a difference between the performance of the different marginal inference methods: standard Gaussians perform worst, while kernel density estimation without copula performs best of all the inference methods considered. Kernel density estimation with independence assumption performs almost on par with the PCE-KDE estimate. This suggests that the KL-RV are close to independent in this case, and that fitting a vine copula (using the available pair copulas) does not improve the overall inference, at least in the considered cases of limited data. It also demonstrates that the true KL-RV distribution is not well approximated by independent standard Gaussians nor by other currently available parametric families, but that it can be approximated well by the more flexible kernel density estimation. Parametric inference, offering a variety of standard marginal shapes, performs slightly better than Gaussian random variables.
<figure><img src="image_106.png" /><figcaption>**Figure 6.9:** Convergence of average normalized Wasserstein distance. Comparison of the four different methods for inferring the joint distribution of KL-RV described in Step 7 with the method PCE-KDE described in Remark 1. N = 100 and p = 14. Errors are computed based on N<sub>val</sub> = 1,000, R<sub>val</sub> = 10,000, for 50 replications. The gray areas and dashed line have the same meaning as in Figure 6.8.</figcaption></figure>

6.4.2 Borehole function with latent variables

As a second example, we consider the well-known borehole function, which computes the water
flow between two aquifers that are connected by a borehole (Harper and Gupta, 1983). It
depends on eight parameters and is defined by

$$
B(R_w, H_u, K_w, R, T_u, T_l, H_l, L) = \frac{2\pi T_u (H_u - H_l)}{\ln(R/R_w) \left(1 + \frac{2LT_u}{\ln(R/R_w)R_w^2 K_w} + \frac{T_u}{T_l}\right)}. \quad (6.23)
$$

Its input random variables and their distributions are provided in Table 6.1.

We consider five parameters ($\Xi = (R, T_u, T_l, H_l, L)$) of the borehole function to be latent, resulting in the three-dimensional stochastic simulator $\tilde{B}(r_w, h_u, k_w) = B(r_w, h_u, k_w; \Xi)$.

For the three-dimensional input space, we use $N = \{20, 30, 60\}$ input samples and a maximal PCE degree of $p = 6$. The accuracy of the borehole approximation in terms of relative mean-squared validation error is in the order of $10^{-3}/10^{-7}/10^{-10}$ for the different experimental design sizes. The number of trajectories is in the range $R = \{10, 30, 100, 300\}$.

This results in typically *M* = 2 KL modes for an explained variance threshold of 99.9%. The eigenvalues of the KLE are approximately $\lambda_1 \approx 170$ and $\lambda_2 \approx 0.5$. The first mode alone covers more than 99.5% of the total variance, even though five independent parameters are used as latent variables. Two of these have a significant total Sobol' index, and the sum of the first-order indices of the latent group is 19%. We will investigate this phenomenon in more detail in Section 6.5 below.

As before, we analyze the global convergence of the marginal and covariance approximation for
increasing numbers of input samples and trajectories, and we compare different methods for
inferring the distribution of the KL-RV as described in Step 7 of our algorithm (Section 6.3).
For the detailed explanation of the error measures, the setup of the convergence study, and the
interpretation of the plots, see Sections 6.4 and 6.4.1.

In Figure 6.10a, we see that our method converges in both global error metrics ($\epsilon_{\text{marg}}$ and $\epsilon_{\text{cov}}$)
towards the rough empirical lower bound indicated by the gray area and dashed line. For $\epsilon_{\text{marg}}$,
<table><caption>**Table 6.1:** Borehole function: Input random variables and their distributions. For the borehole stochastic simulator with hidden variables, five of the eight variables (marked by italic letters) are considered latent.</caption><thead><tr><th>Variable</th><th>Distribution</th><th>Description</th><th>Total Sobol' index</th></tr></thead><tbody><tr><td>R<sub>w</sub></td><td>N(0.10, 0.0161812)</td><td>borehole radius</td><td>6.94e-01</td></tr><tr><td>H<sub>u</sub></td><td>U([990, 1110])</td><td>potentiometric head of upper aquifer</td><td>1.06e-01</td></tr><tr><td>K<sub>w</sub></td><td>U([9855, 12045])</td><td>borehole hydraulic conductivity</td><td>2.51e-02</td></tr><tr><td>R</td><td>Lognormal([7.71, 1.0056])</td><td>radius of influence</td><td>2.77e-06</td></tr><tr><td>T<sub>u</sub></td><td>U([63070, 115600])</td><td>transmissivity of upper aquifer</td><td>2.10e-08</td></tr><tr><td>T<sub>l</sub></td><td>U([63.1, 116])</td><td>transmissivity of lower aquifer</td><td>8.23e-06</td></tr><tr><td>H<sub>l</sub></td><td>U([700, 820])</td><td>potentiometric head of lower aquifer</td><td>1.06e-01</td></tr><tr><td>L</td><td>U([1120, 1680])</td><td>borehole length</td><td>1.03e-01</td></tr></tbody></table>

the difference between the results for the three experimental design sizes $N = 20, 30$, and $60$ is small. This indicates that at least for this model, a validation mean-squared error smaller than $O(10^{-3})$ does not lead to significantly more accurate results, and that below this accuracy the error is dominated by the inference error.

The convergence of the emulated covariance function is displayed in Figure 6.10c. As expected, the global error becomes smaller with an increasing number of trajectories, and approaches the lower bound representing the variability due to the error being computed from samples. Again, the difference in the results for the three experimental design sizes $N = 20, 30, 60$ is small. Since the first mode accounts for more than 99.5% of the explained variance, the first KLE eigenfunction has the dominating influence on the covariance estimation (Eq. (6.9)). The results indicate that the first eigenfunction and its eigenvalue are estimated accurately already for the smallest experimental design sizes.

The comparison of the different inference methods for the distribution of the KL-RV (Step 7 of the algorithm) is displayed in Figure 6.11. Similarly as for the Ishigami function, we observe that for a small number of trajectories ($R = 10$ and $30$) the four inference methods and PCE-KDE show almost the same performance. Modeling the KL-RV with standard Gaussian distributions seems to offer a slight advantage (resulting in a slightly smaller median error and smaller variability) for small numbers of trajectories, probably because they make the strongest assumptions on the distribution shape, which is advantageous for generalizability in the case of small data. As the number of trajectories grows, a similar pattern as in Section 6.4.1 emerges: standard Gaussian inference shows the worst performance, followed by parametric inference. Inference with kernel density estimation (dependent and independent) shows the best performance, on par with the PCE-KDE estimate, which (due to the high accuracy of the PCE approximations for $N = 60$) represents the near-optimal estimate given the available training data.

Interestingly, there is no significant difference between the performance of KDE with and without the independence assumption. Here the magnitude of the eigenvalues might offer an explanation: with more than two orders of magnitude difference between $\lambda_1$ and $\lambda_2$, the dependence between the two random variables $\xi_1$ and $\xi_2$ does not influence the resulting predictions as much as the correct identification of the marginal shape of the first KL-RV $\xi_1$.
<figure><img src="image_107.png" /><figcaption>**Figure 6.10:** Convergence of the $\epsilon_{\text{marg}}$ and $\epsilon_{\text{cov}}$ (Eqs. (6.20) and (6.21)) for an increasing number of available trajectories and parameter locations. Results for the stochastic emulator with parametric inference (Option 7b) and 50 replications. The errors are computed based on a validation set of size $N_{\text{val}} = 1,000$, $R_{\text{val}} = 10,000$. The gray areas and the dashed line represent quantiles and the median of a lower bound estimate for the respective error measure computed from 100 × 2 independent MC samples of size $R_{\text{val}} = 10,000$ generated by the true stochastic simulator. The colored points in Figure 6.10a correspond to the results for a single replication and validation point as shown in Figure 6.10b and help assess the meaning of the numerical error measures.</figcaption></figure>

<figure><img src="image_108.png" /><figcaption>**Figure 6.11:** Convergence of $\epsilon_{\text{marg}}$. Comparison of the four different methods for inferring the joint distribution of KL-RV described in Step 7 with PCE-KDE described in Remark 1. $N = 60$ and $p = 6$. Errors are computed based on $N_{\text{val}} = 1,000$, $R_{\text{val}} = 10,000$, and 50 replications. The gray areas and dashed line have the same meaning as in Figure 6.10.</figcaption></figure>
### 6.4.3 Heston stochastic volatility model for a stock price

As a third example, we consider the Heston stochastic volatility model, which describes a stock price $Y_t$ (Heston, 1993) with its volatility $\nu_t$ modeled as stochastic process:

$$dU_t = \mu U_t dt + \sqrt{\nu_t} U_t dW_t^{(1)}, \quad (6.24)$$

$$d\nu_t = \kappa(\theta - \nu_t)dt + \sigma\sqrt{\nu_t}dW_t^{(2)} \quad (6.25)$$

with two Wiener processes $W_t^{(1)}$ and $W_t^{(2)}$ with correlation coefficient $\rho$. This model has six uniformly distributed parameters $\mathbf{X} = (\mu, \kappa, \theta, \sigma, \rho, \nu_0)$ detailed in Table 6.2, the bounds of which are calibrated from real data as described in Zhu and Sudret (2021c). The quantity of interest is

$$Y_x = U_1(\mathbf{X} = \mathbf{x}), \quad (6.26)$$

i.e., the stock price after 1 year. As proposed by Zhu and Sudret (2021c), we set $U_0 = 1$ and use the Euler-Maruyama method to integrate the system of stochastic differential equations (SDEs) and replace $\nu_t$ by $\max(\nu_t, 0)$ to avoid negative values of $\nu_t$. This model is a stochastic simulator due to the stochasticity induced by the two Wiener processes $W_t^{(1)}$ and $W_t^{(2)}$ driving the SDEs. A trajectory in the parameter space $\mathcal{D}$ is obtained by fixing the realizations of these processes and evaluating Eqs. (6.24) and (6.25) for $\mathbf{x} \in \mathcal{D}$.

For the six-dimensional input space, we use $N = \{50, 100, 150\}$ input samples and a maximal PCE degree of $p = 7$. The accuracy of the approximation in terms of relative mean-squared validation error is ca. $\mathcal{O}(0.03)/\mathcal{O}(0.02)/\mathcal{O}(0.006)$ for the different experimental design sizes. This means that the Heston model is not particularly well approximated by PCE, even for rather large experimental designs. We use a number of trajectories in the range $R = \{10, 30, 100, 300\}$.

This results in typically $M = 4$ to 6 KL modes for an explained variance threshold of 99.9%. The first eigenvalue is $\lambda_1 \approx 0.05$ and usually covers more than 97% of the variance.

<table><caption>Table 6.2: Parameters and their distributions for the Heston SDE model.</caption><thead><tr><th>Variable</th><th>Distribution</th><th>Description</th></tr></thead><tbody><tr><td>μ</td><td>U([0, 0.1])</td><td>Expected return rate</td></tr><tr><td>κ</td><td>U([0.3, 2])</td><td>Mean reversion speed of the volatility</td></tr><tr><td>θ</td><td>U([0.02, 0.07])</td><td>Long term mean of the volatility</td></tr><tr><td>σ</td><td>U([0.2, 0.4])</td><td>Volatility of the volatility</td></tr><tr><td>ρ</td><td>U([-1, -0.5])</td><td>Correlation coefficient between dW<sub>t</sub><sup>(1)</sup> and dW<sub>t</sub><sup>(2)</sup></td></tr><tr><td>ν<sub>0</sub></td><td>U([0.02, 0.07])</td><td>Initial volatility</td></tr></tbody></table>

Again, we analyze the global convergence of the marginal and covariance approximation in the same way as in the preceding sections. The marginal approximations of the parametric stochastic emulator converge with increasing experimental design size and number of trajectories, but slowly, as displayed in Figure 6.12a. There is no significant difference between the three experimental design sizes $N = 50, 100, 150$. This indicates that the improvement due to a better PCE approximation for an increasing number of experimental design points is overshadowed by the inaccuracy due to the inference of the KL-RV. This, in turn, could be because the PCE approximation is not yet sufficiently accurate (note that the relative validation error is in the
order of $10^{-2}$ for all ED sizes.) Even for the largest number of trajectories, the averaged and
normalized Wasserstein distance between the responses of the true model and the emulator is
still much larger than the variability resulting from sampling the true model, which is illustrated
by the gray areas and the dashed line in Figure 6.12a (quantiles and median, respectively).
Comparing the boxplots to the colored points corresponding to the marginal estimates illustrated
in Figure 6.12b, we observe that the marginal shape of the stochastic response for one validation
point $\mathbf{x}_{\text{val}}^{(1)}$ is already captured quite well for 100-300 trajectories (with the value for $R = 300$
being a bit of an outlier).

The convergence of the covariance function is shown in Figure 6.12c. As expected, the covariance
estimate becomes better with increasing number of trajectories, even approaching the lower
bound obtained by resampling the original model. However, we observe again that an increasing
number of experimental design points does not influence the estimate much, which indicates that
the covariance estimate is quite robust against the trajectory approximation quality. Since the
first mode is also dominant for this example (accounting for ca. 97% of the explained variance),
it indicates that the first eigenfunction is well estimated already for small experimental design
sizes.

Figure 6.13 shows the comparison between the different methods for KL-RV inference (Step 7 of
the algorithm in Section 6.3) as described in Section 6.4.1.4. All four methods perform compara-
bly. The independent standard Gaussian approximation performs slightly better than the other
methods in the case of few trajectories and slightly worse for the case of many trajectories, which
is consistent with the previous cases. Again, KDE with and without dependence shows almost
identical performance, indicating that either the true copula is the independence copula, or
that the existing parametric copulas are not suitable for capturing the dependence structure.
While the inference methods show a similar performance to PCE-KDE for smaller numbers of
trajectories, PCE-KDE finds a better marginal approximation when $R = 300$ trajectories are
available. This indicates that some information is lost in the KLE procedure of Section 6.3.
<figure><img src="image_109.png" /><figcaption>**Figure 6.12:** Convergence of $\epsilon_{\text{marg}}$ and $\epsilon_{\text{cov}}$ (Eqs. (6.20) and (6.21)) for an increasing number of available trajectories and parameter locations. Results for the stochastic emulator with parametric inference (choice 7b) and 50 replications. The errors are computed based on a validation set of size $N_{\text{val}} = 1,000$, $R_{\text{val}} = 10,000$. The gray areas and the dashed line represent quantiles and the median of a lower bound estimate for the respective error measure computed from 100 × 2 independent MC samples of size $R_{\text{val}} = 10,000$ generated by the true stochastic simulator. The colored points in Figure 6.12a correspond to the results for a single replication and validation point as shown in Figure 6.12b and help assess the meaning of the numerical error measures.</figcaption></figure>

<figure><img src="image_110.png" /><figcaption>**Figure 6.13:** Convergence of marginals (Wasserstein distance). Comparison of the four different methods for inferring the joint distribution of KL-RV described in Step 7 with PCE-KDE described in Remark 1. $N = 150$ and $p = 5$. Errors are computed based on $N_{\text{val}} = 1,000$, $R_{\text{val}} = 10,000$, and 50 replications. The gray areas and dashed line have the same meaning as in Figure 6.12.</figcaption></figure>
## 6.5 Considerations on the number of modes

In this section we investigate how many modes $K$ we can expect in the stochastic emulator of Eq. (6.18). We consider here a certain class of stochastic simulators that arise from a deterministic model $z \mapsto \mathcal{M}(z)$ by considering some of its variables as *hidden* (or *latent*). In other words, the stochastic simulator is $\mathcal{M}(\mathbf{X}, \Xi)$, where $\mathbf{X}$ are the explicit input variables, and $\Xi$ the latent variables inducing the stochasticity (the random events, see Section 6.2.1). Assume that all components of $\mathbf{Z} = (\mathbf{X}, \Xi)$ are independent, and denote by $f_k$ the marginal distribution of $Z_k$.

Assume further that the deterministic simulator $\mathcal{M}$ has finite variance under $f_Z$. Then it can be decomposed into the Hoeffding-Sobol' decomposition (a.k.a. ANOVA decomposition, analysis of variance) (Hoeffding, 1948; Sobol and Gresham, 1995) as

$$ \mathcal{M}(\mathbf{Z}) = m_0 + \sum_{1 \le i \le d} m_i(Z_i) + \sum_{1 \le i < j \le d} m_{i,j}(Z_i, Z_j) + \dots + m_{1, \dots, d}(Z_1, \dots, Z_d) \quad (6.27) $$

where the terms satisfy $\int m_I(\mathbf{Z}_I) f_k(z_k) dz_k = 0$ for all $k \in I \subset \{1, \dots, d\}$. $m_0$ is the mean of $\mathcal{M}(\mathbf{Z})$. The univariate terms $m_i$ are called main effects, and the remaining summands are interaction terms of increasing order.

Now we group the summands of Eq. (6.27) according to whether they involve only input variables, only latent variables, or some variables from both groups. This results in the following decomposition:

$$ \mathcal{M}(\mathbf{X}, \Xi) = m_0 + \mathcal{M}_1(\mathbf{X}) + \mathcal{M}_2(\Xi) + \mathcal{M}_{12}(\mathbf{X}, \Xi) \quad (6.28) $$

where, e.g., $\mathcal{M}_1(\mathbf{x}) = \sum_{I:Z_I \subset \mathbf{X}} m_I(z_I)$. The last summand of Eq. (6.28) contains all interaction terms from Eq. (6.27) that involve at least one input and at least one latent variable.

It is a rather common phenomenon in uncertainty quantification that engineering models actually have near-zero interaction terms. In that case, the model is essentially additive with respect to the two groups of variables $\mathbf{X}$ and $\Xi$:

$$ \mathcal{M}(\mathbf{X}, \Xi) \approx m_0 + m_1(\mathbf{X}) + m_2(\Xi). \quad (6.29) $$

This means that any realization $\xi$ of the unknown latent variables $\Xi$ results in a constant shift of $\mathcal{M}(\cdot, \xi)$ regardless of the value of the input parameters $\mathbf{x}$, a behavior that can be accurately modeled by a single KL mode: if equality holds in Eq. (6.29), the mean function is $\mu(\mathbf{x}) = \mathbb{E}_{\Xi}[\mathcal{M}(\mathbf{x}, \Xi)] = m_0 + m_1(\mathbf{X})$, the covariance function is $c(\mathbf{x}, \mathbf{x}') = \text{Var}[m_2(\Xi)]$, and the only nonzero eigenvalue of Eq. (6.7) is $\lambda_1 = \text{Var}[m_2(\Xi)]$ with eigenfunction $\phi_1(\mathbf{x}) = 1$.

We have observed this in the numerical examples in Section 6.4, e.g., for the case of the borehole model with hidden variables. The deterministic borehole model defined in Eq. (6.23) has relatively low interaction: despite its nonlinearity, under the input uncertainties in Table 6.1 its first order Sobol' indices (Sobol', 1993) sum up to $\sum_{i=1}^8 S_i^1 \approx 96.7\%$. The interaction effect between the explicit and the latent group is around 2%. Therefore, only one mode is sufficient to model the stochastic simulator that results from treating several of the model's variables as latent.
## 6.6 Discussion and conclusions

We presented a spectral surrogate model for stochastic simulators, a special class of computational models whose response for a given input is a random variable. Our method relies on several advanced techniques for modeling uncertainties, such as polynomial chaos expansion (PCE), Karhunen-Loève expansion (KLE), and statistical inference of multivariate distributions. The resulting surrogate model is not only able to emulate the marginal distributions and the covariance structure, but it can also generate new trajectories.

The form of our surrogate model provides insight into the model structure. The number of expansion modes indicates how high-dimensional the underlying stochasticity is. The eigenfunctions of the KLE, which are polynomials, give information about how the input parameters influence the stochastic simulator output. Even though our numerical examples were chosen to represent a range of cases of increasing complexity, we found that one mode was usually sufficient to explain more than 95% of the variance of the stochastic simulator. We were able to explain this phenomenon for the common case of stochastic simulators that arise from finite-dimensional deterministic models with independent inputs and finite variance by treating some of their input variables as latent. Considering the Hoeffding-Sobol' decomposition of the underlying deterministic simulator, we found that if the interaction terms between the explicit and latent variables are negligible, one single KL mode will be sufficient to emulate the behavior of the stochastic simulator. Indeed, by experience, this is a common occurrence in applications of uncertainty quantification. Interactions are rarely dominant in engineering problems, and so one KL mode might suffice in many cases.

From our numerical experiments, we found that the Gaussian (or more generally, parametric) approximation of marginals of the KL-random variables (KL-RV) can be preferable if the number of available trajectories is very small. When more trajectories are available, the better choice is kernel density estimation. Since the first mode was dominant in our numerical examples, the characterization of the dependence between KL-RV turned out not to be crucial, at least for the class of applications considered.

Our numerical tests reveal that the emulator is able to capture the true model behavior reasonably well, as long as enough input samples and trajectories are used. Due to the sequential nature of our approach, it is important to use enough points in the experimental design: if the PCE approximation is not accurate enough, also the inferred distribution for the KL-RV will be inaccurate. Interestingly, we observed in all three examples that the covariance estimate was not heavily influenced by using a larger experimental design, even though the latter typically results in more accurate PCE approximations of the trajectories. This indicates that the number of trajectories is more important for the covariance approximation than the quality of the PCE approximation. Also, it seems that (since the first mode was dominant for all investigated models) the first KLE eigenfunction can be identified accurately already from a small experimental design.

Note that our surrogate relies on the assumption that the trajectories are well approximated by sparse PCE, an assumption not fulfilled by some stochastic simulators, e.g., ones with discontinuous or non-differentiable trajectories in the input parameter space. This could be circumvented