# Sparse spectral surrogate models for deterministic and stochastic computer simulations

Nora Lüthen

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nora Lüthen. Sparse spectral surrogate models for deterministic and stochastic computer simulations. Statistics [math.ST]. ETH Zürich, 2022. English. NNT: tel-04284700

HAL Id: tel-04284700

https://hal.science/tel-04284700

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Sparse spectral surrogate models for
deterministic and stochastic
computer simulations

A thesis submitted to attain the degree of

DOCTOR OF SCIENCES of ETH ZURICH

(Dr. sc. ETH Zurich)

presented by

Nora Lüthen

M.Sc. Mathematics, University of Bonn

born on 26.03.1991

accepted on the recommendation of

Prof. Dr. Bruno Sudret, examiner

Dr. Sebastián Da Veiga, co-examiner

Dr. Stefano Marelli, co-examiner

Prof. Dr. Olivier Roustant, co-examiner

15.09.2022
<figure><img src="image_3.png" /></figure>
Abstract

Computer simulations are an invaluable tool for modeling and investigating real-world phenomena and processes. However, as any model, simulations are affected by uncertainty caused by imperfect knowledge or natural variability of their parameters, initial conditions, or input values. This leads to uncertainty in the model response, which needs to be quantified to make subsequent conclusions and decisions trustworthy. To alleviate the considerable cost of uncertainty analyses for expensive computational models, the latter are often replaced by surrogates, i.e., by approximations with an explicit functional form that can be created based on a rather small number of model evaluations, and can be evaluated at low cost.

Some computer models are affected by uncertainty only through their input parameters: for fixed values of the inputs, they always produce the same response. These models are called deterministic simulators. In contrast, models that feature inherent stochasticity are called stochastic simulators. The latter generate a different result each time they are run even if their input parameters are held at fixed values. In other words, they behave like random fields whose index set is the space of input parameters.

In this thesis, we investigate spectral surrogate models, which are a class of global non-intrusive methods that expand the computational model onto an orthonormal basis of a suitable function space. We focus on sparse expansions, i.e., representations that only include a small finite subset of the basis elements. Sparse representations are typically computed by regression with sparsity-encouraging constraints, often using ideas originating from the field of compressed sensing.

In particular, for deterministic simulators we explore the popular sparse polynomial chaos expansions (PCE) method, which utilizes a polynomial basis that is orthonormal with respect to the distribution of the input variables. We conduct an extensive literature survey as well as a benchmark of several promising methods on multiple models of varying dimensionality and complexity. The benchmark results are aggregated and visualized in a novel way to extract reliable recommendations about which methods should be used in practice.

We also investigate the recently proposed Poincaré chaos expansions, which rely on a generally non-polynomial basis consisting of eigenfunctions of a specific differential operator connected to the Poincaré inequality. By construction, this basis is well suited for derivative-based global sensitivity analysis, which we explore both analytically and numerically.

Furthermore, we propose a new surrogate model for stochastic simulators. Taking the random function view of a stochastic simulator, we approximate its trajectories by sparse PCE and perform Karhunen-Loève expansion on them. The latter is a well-known spectral representation for a random field which separately characterizes its spatial and stochastic variation. The joint distribution of the random coefficients is inferred using the marginal-copula framework. The resulting surrogate model is able to approximate marginal distributions, mean, and covariance function of the stochastic simulator, and can generate new trajectories.
Zusammenfassung

Computersimulationen sind ein wertvolles Werkzeug für die Modellierung und Erforschung realer
Phänomene und Prozesse. Allerdings sind Simulationen, wie alle Modelle, von Ungewissheit
betroffen, die durch fehlende Informationen oder durch die natürliche Variabilität ihrer Param-
eter, Anfangsbedigungen oder Eingabewerte verursacht wird. Dies führt zu Ungewissheit in der
Modellantwort, welche quantifiziert werden muss, damit anschliessende Schlussfolgerungen und
Entscheidungen verlässlich sind. Um die beträchtlichen Kosten von Ungewissheitsanalysen für
teure Computermodelle zu verringern, werden letztere häufig durch Ersatzmodelle (engl. surro-
gates) angenähert, d.h. durch Funktionen mit expliziter analytischer Form, die mithilfe weniger
Modellauswertungen erstellt und kostengünstig ausgewertet werden können.

Einige Computermodelle sind nur durch ihre Eingabeparameter von Ungewissheit betroffen: für
eine gegebene Kombination von Eingabewerten erzeugen sie immer dieselbe Ausgabe. Solche
Modelle heissen deterministische Simulatoren. Im Gegensatz dazu werden Modelle mit in-
härenter Stochastizität stochastische Simulatoren genannt. Diese erzeugen bei jeder Auswertung
eine andere Ausgabe, selbst wenn die Eingabewerte immer gleich sind. Mit anderen Worten,
sie verhalten sich wie Zufallsfelder (engl. random fields), wobei der Raum der Eingabewerte als
Indexmenge fungiert.

In dieser Arbeit untersuchen wir spektrale Ersatzmodelle, eine Klasse globaler nicht-intrusiver
Methoden, welche das Computermodell mithilfe einer orthonormalen Basis eines geeigneten
Funktionenraums darstellen. Wir konzentrieren uns auf dünnbesetzte (engl. sparse) Darstellun-
gen, d.h. solche, die nur eine geringe endliche Anzahl Basisfunktionen verwenden. Dünnbesetzte
Darstellungen erhält man üblicherweise mithilfe von Regression mit Nebenbedingungen, welche
die Dünnbesetztheit begünstigen, wobei oft Ideen aus dem Gebiet der komprimierten Erfassung
(engl. compressed sensing) verwendet werden.

Für deterministische Simulatoren haben wir uns speziell mit der beliebten Methode der dünnbe-
setzten polynomiellen Chaos-Expansion (PCE) beschäftigt, welche eine polynomielle Basis ver-
wendet, die orthonormal bezüglich der Eingabeverteilung ist. Wir erstellen eine umfassende
Literaturübersicht und führen eine Vergleichsstudie vielversprechender Methoden auf mehreren
Modellen unterschiedlicher Dimension und Komplexität durch. Die Ergebnisse werden auf eine
neue Weise zusammengeführt und visualisiert, um zuverlässige Empfehlungen zur Methodenwahl
in der Praxis zu erhalten.

Wir untersuchen auch die relativ neue Methode der Poincaré-Chaos-Expansion, die eine im
Allgemeinen nicht polynomielle Basis verwendet, welche aus Eigenfunktionen eines Differenzial-
operators mit Verbindungen zur Poincaré-Ungleichung besteht. Diese Basis ist per Definition
gut für die gradientenbasierte globale Sensitivitätsanalyse geeignet, was wir sowohl analytisch
als auch numerisch untersuchen.

Ausserdem entwickeln wir ein neues Ersatzmodell für stochastische Simulatoren. Unter Verwen-
dung der Zufallsfunktionen-Sicht eines stochastischen Simulators nähern wir dessen Trajektorien
mit dünnbesetzten PCE an und berechnen die Karhunen-Loève-Expansion derselben. Letztere
ist eine bekannte spektrale Darstellung für Zufallsfelder, welche die räumliche und stochastische
Variabilität getrennt modelliert. Die multivariate Verteilung der zufälligen Koeffizienten wird
mithilfe von Marginalverteilungen und Copulas bestimmt. Das resultierende Ersatzmodell kann Marginalverteilungen sowie Mittelwert- und Kovarianzfunktion des stochastischen Simulators annähern und neue Trajektorien erzeugen.
Acknowledgments

Life is random and uncertain. Whether or not someone has the opportunity to achieve something is influenced by many random events and decisions. Of course, my decision to pursue a PhD in the group of Prof. Bruno Sudret at ETH Zürich was a conscious and determined one, but there were lots of chances and lucky coincidences leading up to it. Certainly, many people had (and are having) a great and deterministic impact on my life, changing it to the better. I feel extraordinarily lucky and grateful to have them in my life.

First and foremost, I want to express my sincere gratitude to my supervisor Prof. Bruno Sudret. Thank you for accepting me into your group and for giving me the opportunity to work with you on these exciting problems. Thank you for your guidance, encouragement, and trust throughout these years, and for all the opportunities (writing, teaching, presentations, collaborations) and feedback you provided that allowed me to develop my skills and grow professionally. Working with you and in your group has been an exceptionally positive experience, and I am very grateful for it.

I was very lucky to have Dr. Stefano Marelli as my second supervisor. Thank you Stefano for the patient support and sound advice in all matters of PhD life, for the (sometimes heated, always super interesting) scientific discussions and valuable feedback, and for always having an “offenes Ohr” and being extremely generous and helpful whenever needed. I feel very lucky that I had the opportunity to work with you and learn from you during these years.

Furthermore, I want to thank Prof. Sebastián Da Veiga and Prof. Olivier Roustant for serving on my doctoral committee. Thank you very much for the time and energy you put into reading my thesis and into the very interesting, insightful discussion during my defense. To be honest, I would not have minded discussing with you even longer!

I would like to thank Prof. Olivier Roustant and Prof. Fabrice Gamboa for our inspiring collaboration. Thank you for giving me the opportunity to stay a month with you in Toulouse, and for the time and energy you took for our discussions. It is an honor to work with you!

I had the pleasure of spending my PhD years with a wonderful group of smart, fun and friendly colleagues. I immensely enjoyed our lunch discussions, coffee breaks, dinners, and all our other shared adventures. Thank you Adéla, Anderson, Andréia, Biswarup, Chiara, Christos, Damar, Diego, Emiliano, Fritz, Katerina, Marco, Margaretha, Max, Moustapha, Nicole, Paul, Philippe, René, Roberto, Rui, Stefano, Styfen, Tong, and Xujia! I especially want to thank Xujia, my PhD sibling: I am so happy that we were hired for the same project. Thank you for our many interesting discussions and shared highs and lows. I very much enjoyed going on this journey with you! Paul, my officemate for more than three years, thank you so much for integrating me into the group and to the Civil Engineering world in general, and for all the fun times and useful advice! Chiara, thank you for your moral and culinary support during the last stressful months—it was a pleasure to share the office with you!

Furthermore, I would like to thank Dr. Hélène Schernberg for taking the time to share her valuable experience and thoughtful advice regarding science and life choices with me. This has helped and continues helping me immensely.
During my studies, I had the honor of working with Prof. Martin Rumpf and his group who guided me on my first small steps in mathematics research—thank you for this enjoyable experience that inspired me to pursue a PhD! Furthermore, I would like to thank my teacher Ralf Wambach, who nearly 20 years ago patiently started encouraging me to join mathematics camps and competitions. Without him, I would almost surely have followed a different trajectory.

Of course, work needs to be balanced with free time, and I am grateful to all my wonderful friends in Switzerland and abroad for providing this balance. Thank you for distracting me with dinners, hikes, gardening, skiing, Ultimate Frisbee, or phone calls whenever the work thoughts became heavy, for keeping in touch and sharing the ups and downs, and generally for all the good times we had during all the last years!

I feel very grateful to have the best family one could wish for. Thank you to my parents and sisters for believing in me and always providing the right loving kind of understanding, encouragement, and advice.

Finally, thank you Jerri for your constant and patient support, the good food you cook and your wonderful presence in my life.
# Contents

<table><tr><td>1</td><td>Overview</td><td>1</td></tr><tr><td>1.1</td><td>Motivation and context</td><td>1</td></tr><tr><td>1.2</td><td>Research objectives</td><td>4</td></tr><tr><td>1.3</td><td>Outline</td><td>6</td></tr><tr><td>2</td><td>Review of the state of the art</td><td>7</td></tr><tr><td>2.1</td><td>Probabilistic uncertainty quantification</td><td>7</td></tr><tr><td>2.2</td><td>Surrogate modeling for deterministic simulators</td><td>9</td></tr><tr><td>2.3</td><td>Surrogate modeling for stochastic simulators</td><td>14</td></tr><tr><td>2.4</td><td>Spectral methods</td><td>21</td></tr><tr><td>2.5</td><td>Sparsity and compressed sensing</td><td>30</td></tr><tr><td>2.6</td><td>Sensitivity analysis</td><td>36</td></tr><tr><td>2.7</td><td>Software for uncertainty quantification</td><td>40</td></tr><tr><td>3</td><td>Sparse polynomial chaos expansions: Literature survey and benchmark</td><td>43</td></tr><tr><td>3.1</td><td>Introduction</td><td>44</td></tr><tr><td>3.2</td><td>Framework and literature survey for sparse polynomial chaos expansions</td><td>45</td></tr><tr><td>3.3</td><td>Numerical results</td><td>56</td></tr><tr><td>3.4</td><td>Discussion and conclusions</td><td>67</td></tr><tr><td>3.A</td><td>Details on experimental design sampling techniques</td><td>79</td></tr><tr><td>3.B</td><td>Details on sparse regression solvers</td><td>90</td></tr><tr><td>3.C</td><td>Additional results</td><td>100</td></tr><tr><td>3.D</td><td>Benchmark studies</td><td>100</td></tr></table>
4 Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications 109

<table><tr><td>4.1 Introduction</td><td>110</td></tr><tr><td>4.2 Sparse polynomial chaos expansions</td><td>112</td></tr><tr><td>4.3 Numerical results</td><td>120</td></tr><tr><td>4.4 Conclusion and discussion</td><td>130</td></tr><tr><td>4.A Details on the settings for the basis adaptivity schemes in the benchmark</td><td>137</td></tr><tr><td>4.B Additional results</td><td>137</td></tr></table>

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

<table><tr><td>5.1 Introduction</td><td>146</td></tr><tr><td>5.2 Mathematical background</td><td>148</td></tr><tr><td>5.3 Computation of sparse Poincaré expansions</td><td>156</td></tr><tr><td>5.4 Numerical results</td><td>161</td></tr><tr><td>5.5 Conclusion</td><td>170</td></tr><tr><td>5.A Additional results</td><td>177</td></tr></table>

6 A spectral surrogate model for stochastic simulators computed from trajectory samples 181

<table><tr><td>6.1 Introduction</td><td>182</td></tr><tr><td>6.2 Theoretical foundation</td><td>186</td></tr><tr><td>6.3 Surrogating a stochastic simulator from a set of samples</td><td>191</td></tr><tr><td>6.4 Numerical experiments</td><td>195</td></tr><tr><td>6.5 Considerations on the number of modes</td><td>210</td></tr><tr><td>6.6 Discussion and conclusions</td><td>211</td></tr><tr><td>6.A Analytical derivations for extended KLE on PCE trajectories</td><td>219</td></tr></table>

7 Discussion and conclusions 223

<table><tr><td>7.1 Sparse polynomial chaos expansions</td><td>224</td></tr><tr><td>7.2 Poincaré chaos expansions</td><td>229</td></tr><tr><td>7.3 Spectral surrogate for stochastic simulators</td><td>230</td></tr><tr><td>7.4 Final conclusion</td><td>231</td></tr></table>

Bibliography 233
# CHAPTER 1

## Overview

### 1.1 Motivation and context

Computer technology has made incredible progress. Just 80 years ago, only a handful of people in the world were concerned with electronic computing machines. The first computers of that time were big devices filling whole rooms, and to enter a program into such a machine, the operators had to physically manipulate its wires and switches. One generation ago, when my parents wrote their dissertations using typewriters, computer technology had already advanced a lot. Powerful operating systems and programming languages had been developed, the first personal computers had become available, and around the time I was born, the world wide web started out as a project of a few especially tech-savvy people. Nowadays, the majority of people in our society owns at least a smartphone or a personal computer, not speaking of the various other types of smart devices from watches to cars that are becoming increasingly prevalent. Many of us use the internet every day, even depend on it to do their work and stay connected with friends and family, and a large number of companies and government institutions would not be able to operate anymore without computers. Computers are everywhere.

In addition to the well-known developments that came with computers, such as automation, big data, and anything related to the internet from Wikipedia to social media, there is another domain in which computers critically influence our everyday lives. It is perhaps less familiar to the lay-person, but no less important: namely, the domain of *computer simulation*. This broad term describes the replication of real-world processes and phenomena by mathematical models solved by computer calculations, exploiting the fact that computers can perform billions of calculation steps and process huge amounts of data per second. To name some prominent examples, computer simulations play an essential role in weather forecasting, climate science, and the aerospace industry. Recently, computer simulations of the spreading of infectious diseases have attracted media attention due to the COVID-19 pandemic.

Yet computer simulations, also called *computer models*, are used in many more fields and for
many more purposes. Nowadays, there is hardly a scientific field that does not enrich its theoret-
ical and experimental research with insights and predictions gained from computer models. In
industry, products and processes are optimized using computer simulation. Car manufacturers
optimize the safety and cost of critical components. Insurance companies use computer models
to assess risks and adjust their contract pricing. Even in politics, computer simulations can
be used to aid policy decisions. Much of our modern, comfortable standard of living, coun-
ting on affordable products and services and reliable infrastructure, is facilitated by computer
simulations that have helped explore and optimize the respective system of interest.

In all of the cases described above, computer simulations are crucial because they can give us
information that could not otherwise be obtained: real-world experiments would be unethical,
too expensive or impossible. At the same time, while it would theoretically be possible to
perform the required calculations by hand, this would take many (hundreds of) years.

Of course, to be able to rely on the predictions and responses of a computer simulation, we need
to be able to trust that the values it returns are reasonable. This concerns the mathematical
equations used for the model, the algorithm used to solve them, and the implementation in
computer code: the equations should describe the phenomenon well, the algorithm should be
suitable, and the implementation should be free of bugs.

However, even if the equations, the algorithm, and the implementation are confirmed to be
adequate, there are further sources of inaccuracy or uncertainty that can influence the predictions
of a computer simulation. Every model depends on some parameters and input values that
describe the properties and characteristics of the specific phenomenon of interest: e.g., in a
space flight simulation, such parameters would include the weight, speed, and dimensions of the
spacecraft. In an epidemics simulation, relevant parameters would be the infection and recovery
rates as well as the size of the potentially affected population. However, the natural world is
so complex that we can hardly ever know these parameter values exactly. In other words, the
parameters are affected by uncertainty.

To illustrate this, let us consider an example from structural engineering, namely the compressive
strength of concrete. Concrete is one of the most widely used building materials of our time, and
its strength is a crucial parameter in any simulation of a concrete structure. Perhaps surprisingly
(at least to non-engineers), it turns out that it is virtually impossible to produce two batches
of concrete that have exactly the same strength. Concrete is made of cement, water, and a
granular component called aggregate, and its strength depends on the properties and amounts
of all these ingredients. The same company can produce concrete with the same recipe and still
end up with concrete having a widely varying strength – an uncertainty that is usually dealt
with by using generous safety margins in structural design. Likewise, if the concrete strength
was an input parameter to a computational model, this uncertainty would need to be taken into
account. The uncertainty in the input parameters to a computational model is called *parametric
uncertainty*.

Does the uncertainty in the input parameters affect the output of the computer model, and if
so, how severely? In what way are our subsequent conclusions and decisions affected by this un-
certainty? Answering these questions is the main focus of the field of *uncertainty quantification*
(UQ). More precisely, typical questions tackled in UQ include the following:
* How uncertain is the simulation prediction as a consequence of the uncertainty in the input? (*Uncertainty propagation*)

* Which input uncertainties have the largest influence on the output uncertainty, and which ones are unimportant? (*Sensitivity analysis*)

* What is the probability that the system attains a state that is not considered safe anymore (damaged or failed)? (*Reliability analysis*)

* Which parameters are optimal for a certain objective, if the system is affected by uncertainty? (*Robust optimization*)

* Which (range of) input parameters makes the model output fit best to given observations? (*Bayesian calibration*)

* Which set of parameter combinations provides the most information about the behavior of the model? (*Optimal experimental design*)

Clearly, questions from this list arise in almost any research field that uses computational models. To answer these questions, we need to combine knowledge and ideas from various scientific areas, which makes UQ an interdisciplinary science at the intersection of fields such as computational science, numerical and functional analysis, statistics, computer science, and whichever engineering or applied science field that is interested in answering the above-mentioned questions.

To understand the behavior of a computational model under uncertainty, a single model run is not enough – it must be run many times in different configurations to explore its behavior. The most widely-known uncertainty propagation method, *Monte Carlo simulation*, does exactly this: the computer model is run for a large number of different randomly generated parameter configurations, and the resulting spread of responses is interpreted as the output uncertainty. This method is universal and can be applied readily for computer simulations that only involve a few calculations and can be run in a fraction of a second. However, the typical computer simulation in current applied science and engineering research or in industrial applications is much more resource-intensive, which often makes the Monte Carlo approach overly expensive or even infeasible. Such an expensive computer model often involves the simulation of several components, each of which requires solving complex mathematical equations over a long time scale or with a fine spatial resolution. Performing a single simulation might take hours or days. Even though the computational power of modern computing facilities is increasing steadily due to technological progress, this is compensated for by the increasingly complex state-of-the-art simulations and the associated ever-growing computational demands. Examples of such complex simulations from civil engineering are the response of high-rise buildings to earthquakes, or the power production of an off-shore wind turbine that is subject to varying wind and wave conditions.

In these situations, where it is infeasible to run the computational model sufficiently many times to perform a thorough uncertainty analysis, a convenient solution is to take a detour through a *surrogate model*. A surrogate model is a mathematical function that approximates the input-output relationship of the original computational model. Being a model of a model, it is also sometimes called *metamodel*; other synonyms are *response surface* or *emulator*. It is often trained based on a relatively small number of samples from the original model, and can
be evaluated at a fraction of its cost. Subsequently, the inexpensive surrogate can be used for
uncertainty analysis in place of the expensive original model.

If no knowledge of the original computational model except for a number of input-output pairs for training is required, we call the surrogate *non-intrusive* and the computational model a *black-box*. Non-intrusive surrogates for black-box models are the topic of this thesis. More specifically, we will investigate several instances of a specific class of surrogate models called *spectral expansions*, which combine a strong mathematical foundation with good performance in practical problems. The specific objectives of our research are explained in the following section.

1.2 Research objectives

Before we present the research questions investigated in this thesis, we must briefly introduce the two classes of computational models that we will consider, namely *deterministic* and *stochastic simulators*.

When investigating parametric uncertainty, it is convenient to treat the computer model as a
deterministic simulator: there is no uncertainty within the computational model. Only the input
parameters are assumed to be affected by uncertainty, and this uncertainty is modeled separately.
Surrogating such a deterministic simulator can also be viewed as a function approximation task:
an expensive deterministic function is to be approximated by another cheap-to-evaluate one,
based on a few evaluations of the former. An introduction to surrogates for deterministic models
is provided in Section 2.2.

However, there also exist computer experiments for which it is infeasible or undesirable to com-
pletely separate all the uncertainty from the computational model. For example, for stochastic
differential equation models that involve a driving white noise term, it is usually not meaningful
to treat the latter as a (high-dimensional) input variable, which would be necessary to make
the model deterministic. Such models that are affected by additional stochasticity on top of the
explicitly modeled randomness in the input variables are called *stochastic simulators*. Running
such a simulator repeatedly with identical input parameter values can still result in varying re-
sponses because of the additional stochasticity. We will describe several examples of stochastic
simulators in Section 2.3.

The overarching goal of this thesis is to develop cost-effective surrogate models for deterministic and stochastic simulators, and to investigate their performance and limitations.

Our methodology includes two main components: *spectral expansions* and *sparsity*. *Spectral expansions* are representations of a function in terms of an orthonormal basis of a suitable function space, and we will consider three instances: *polynomial chaos expansion* (PCE), *Poincaré chaos expansion*, and *Karhunen-Loève expansion*. The first and the last ones are well-known since the seminal work by Ghanem and Spanos (1991) establishing stochastic computation in engineering. The second one is a recent development (Roustant et al., 2017, 2020a). All details about these spectral expansions can be found in Section 2.4 below.
To represent a function in terms of a basis, it is expressed as a weighted sum of basis functions. If most of the weights are zero, i.e., if only a few basis functions are involved in the sum, we call the representation *sparse*. Sparsity is a useful concept that can achieve considerable savings in computational cost. In the domain of signal processing, it has led to the development of the fruitful field of *compressed sensing* (Donoho, 2006; Candès et al., 2006). It has also successfully been applied in the field of spectral expansions, notably for polynomial chaos expansions, resulting in *sparse PCE* (Blatman and Sudret, 2008, 2010, 2011; Doostan and Owhadi, 2011; Mathelin and Gallivan, 2012). Sparsity and compressed sensing are discussed in Section 2.5.

Conceptually, the work presented in this thesis can be divided into three parts, each part addressing a different topic in the area of surrogate modeling with sparse spectral methods:

## 1. Benchmark and literature review of methods for computing sparse polynomial chaos expansions (PCE).

(a) The sparse PCE literature is vast, and for practitioners and newcomers it can be difficult to identify the most recent and relevant methodological developments, often resulting in suboptimal choices for their applications. Can we provide a comprehensive and accessible overview of ideas and concepts?

(b) Which of the many sparse PCE computation techniques results in the best surrogate? How does this depend on the characteristics of the model that is to be surrogated?

These questions are addressed in the two journal papers Lüthen et al. (2021) and Lüthen et al. (2022a) reported here in Chapters 3 and 4:

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

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

## 2. Regression-based sparse Poincaré chaos expansions.

What are the advantages of using a different, non-polynomial basis of the function space for the spectral expansion? Which properties does this basis have? How can derivative information be utilized for spectral surrogate modeling and sensitivity analysis?
These questions are addressed in the submitted journal paper Lüthen et al. (2022c) provided in Chapter 5:

* N. Lüthen, O. Roustant, F. Gamboa, B. Iooss, S. Marelli, B. Sudret (2022). Global sensitivity analysis using derivative-based sparse Poincaré chaos expansions.
DOI: 10.48550/arXiv.2107.00394
3. Surrogate modeling for stochastic simulators.

Sparse PCEs show good performance for deterministic simulators. How can their efficiency
be exploited to build a surrogate model for stochastic simulators?

This is addressed in the submitted journal paper Lüthen et al. (2022b) provided in Chap-
ter 6:

* N. Lüthen, S. Marelli, B. Sudret (2022). A spectral surrogate model for stochastic simulators computed from trajectory samples.

DOI: 10.48550/arXiv.2207.05653

1.3 Outline

The thesis is organized as follows. Following this introductory chapter, we describe the current state of the art of surrogate modeling for deterministic and stochastic simulators with a focus on spectral methods and sparsity in Chapter 2. Chapters 3 to 6 contain the journal papers and preprints that resulted from our PhD research. The resulting findings and associated insights as well as connections to other research areas are discussed in Chapter 7.
# CHAPTER 2

## Review of the state of the art

In this chapter, we describe the current state of the art of the most relevant concepts and methods in uncertainty quantification that constitute the foundation for the research developed during this PhD (Chapters 3 to 6). We first give an overview of probabilistic uncertainty quantification and the associated notation in Section 2.1. Then, we review the relevant literature for surrogate modeling for deterministic (Section 2.2) and stochastic simulators (Section 2.3). In both cases, our focus lies on a certain class of methods called *spectral*, which we explain in detail in Section 2.4. Another useful concept in surrogate modeling, complementing spectral methods well, is *sparsity*, which we discuss in Section 2.5. We then give an overview of sensitivity analysis methods in Section 2.6. Finally, in Section 2.7 we present available software implementations of uncertainty quantification methods.

### 2.1 Probabilistic uncertainty quantification

Mathematically, a deterministic computational model can be expressed as a function $\mathcal{M} : \mathcal{D} \to \mathcal{Y}$ that maps from the set $\mathcal{D}$ of admissible inputs to the set $\mathcal{Y}$ of possible outputs (or *responses*). Usually, it is assumed that the input domain $\mathcal{D}$ is a subset of $\mathbb{R}^d$, where $d$ is the dimension of the input space, i.e., the number of scalar parameters. In this thesis, we assume that the output range $\mathcal{Y}$ is a subset of $\mathbb{R}$, i.e., the model is *real-valued*.<fnref n="1" />

When performing uncertainty quantification, the first crucial step is to define and implement the computational model. In the workflow sketched in Figure 2.1 (adapted from Sudret (2007)), this is displayed as Step A. In particular, we need to define the output quantity of interest for the uncertainty analysis: it is often advantageous not to use the full set of output values the computer simulation provides, but instead a meaningful aggregate value or a set of values related to the specific question of interest.

<footnote n="1">¹In general, the response of the model can also be vector-valued or even functional in space or time, which is often addressed by using linear (Nagel et al., 2020; Perrin et al., 2021) or nonlinear (Giovanis and Shields, 2018; Kontolati et al., 2022a) dimension reduction methods in the output space.</footnote>
The second step (Step B in Figure 2.1) in the UQ workflow is to find a suitable mathematical representation of the uncertainty in the input parameters. While there exist several approaches to do so, by far the most common framework and the one used in this thesis is probability theory (Grigoriu, 2002).<fnref n="2" /> We model the input of a deterministic simulator as a d-dimensional random vector X with finite variance and joint probability density function f_X, where d is typically in the order of ten to a hundred.<fnref n="3" /> Here and in the following, uppercase letters denote random variables, whereas lowercase letters denote real numbers (e.g., realizations of a random variable). Boldface letters denote vectors or matrices.

Since the input is a random vector, the output $Y = \mathcal{M}(X)$ is a random variable. We assume that the model output $Y$ has finite variance, which is equivalent to the model $\mathcal{M}$ being in the function space

$$L_{f_X}^2(\mathcal{D}) = \left\{ h : \mathcal{D} \to \mathbb{R} \mid \int_{\mathcal{D}} h(x)^2 f_X(x) dx < \infty \right\}, \quad (2.1)$$

i.e., being square-integrable with respect to the density $f_X$. The uncertainty analysis depicted in Figure 2.1 as Step C is then conducted for the random variable $Y$. The goal of uncertainty propagation is to characterize the moments or even the full distribution of $Y$. Sensitivity analysis, introduced in more detail in Section 2.6, aims at ranking the input parameters by their importance for the variability of the model output. Further analyses could be performed such as reliability analysis, reliability-based design optimization, or Bayesian calibration of parameters from data (see also the list of UQ questions in Section 1.1); however, as they are not discussed in this thesis, they are omitted from Figure 2.1.

To reduce the computational cost of the uncertainty analysis for the original model $\mathcal{M}$, a surrogate model $\tilde{\mathcal{M}}$ might be constructed, which can provide an approximate answer to the questions of Step C. This is illustrated in Figure 2.1 in the form of an additional purple box. The construction and use of a surrogate model can be seen as a detour which aims at reducing the overall computational cost, potentially at the expense of accuracy. We discuss surrogate models for deterministic models in Section 2.2.

So far, we described the mathematical setup for deterministic simulators. In contrast, a stochastic simulator is a computer model $\mathcal{M}$ for which the mapping from the input variables to the response is not deterministic. As before, we assume that the d-dimensional input random vector $X$ has finite variance. However, due to residual stochasticity in the model, now the output $Y_x = \mathcal{M}(x)$ is a random variable even for fixed $x \in \mathcal{D}$. In other words, $\mathcal{M}$ can be modeled as a random field, i.e., a collection of random variables indexed by the set of inputs:

$$\mathcal{M} : \mathcal{D} \times \Omega \to \mathbb{R}$$

- For each $x \in \mathcal{D}$, $\mathcal{M}(x, \cdot) : \Omega \to \mathbb{R}$ is a random variable

- For each $\omega \in \Omega$, $\mathcal{M}(\cdot, \omega) : \mathcal{D} \to \mathbb{R}$ is a function (deterministic simulator)

<footnote n="2">²The uncertainty affecting a system is sometimes further classified as epistemic, i.e., due to incomplete knowledge, and aleatory, i.e., due to natural variability of the parameter in question. While some researchers apply different methods to the two types of uncertainty, using e.g. imprecise probability for epistemic uncertainty (Faes et al., 2021), we model both types probabilistically. Note also that the classification can depend on the point of view, as illustrated by the example of the compressive strength of concrete in Section 1.1.</footnote>

<footnote n="3">³For an overview of the state-of-the-art treatment of very high-dimensional input in UQ, we refer to Kontolati et al. (2022b).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.1:** A diagram illustrating the three steps of an uncertainty analysis, adapted from Sudret (2007). An additional purple box shows the surrogate model (here illustrated with the functional form of PCE), which after its construction can be used in place of the original model to perform the UQ analysis. While there are more areas of UQ that could be included here, such as reliability analysis or Bayesian inference, this diagram focuses on the topics discussed in this thesis.</figcaption></figure>

Here, $\omega \in \Omega$ is a random event responsible for the stochasticity of the simulator. It is not explicitly modeled and instead considered to be an unobservable, latent parameter whose value is not known. In the case of computer simulations, it is usually possible to fix $\omega$ by fixing the seed of the pseudo-random number generator, which results in the same sequence of random numbers being generated. However, the seed is in general not a useful parametrization of the latent stochasticity and should not be treated as an input variable. We will discuss stochastic simulators, random fields, and the associated surrogate models in Section 2.3.

## 2.2 Surrogate modeling for deterministic simulators

We consider the case of real-valued deterministic simulators $\mathcal{M} : \mathcal{D} \to \mathcal{Y} \subset \mathbb{R}$. Uncertainty analysis is concerned with the properties of the induced random variable $Y = \mathcal{M}(\mathbf{X})$. There are two possible approaches to constructing a surrogate model for $Y$:

* The response random variable $Y$ could be approximated, aiming for equality in distribution: $\tilde{Y} \stackrel{d}{=} Y$.

* The mapping $\mathcal{M}$ could be approximated, aiming for almost sure equality $\tilde{\mathcal{M}}(\mathbf{X}) \stackrel{\text{a.s.}}{=} Y = \mathcal{M}(\mathbf{X})$.

The second approach is nowadays much more widespread. One reason for this might be that when constructing surrogate models for black-box models, we usually assume to have at our disposal a number of input-output pairs generated by the model. It is therefore natural to use this information to approximate the mapping itself instead of only the resulting distribution.

In the context of deterministic simulators, a surrogate model $\tilde{\mathcal{M}}$ is a mathematical function that approximates the input-output relationship of the computational model, which is considered a black-box. The surrogate model is computed based on a set of *training data* consisting of input-output pairs, also called *experimental design*, of the original deterministic simulator.<fnref n="4" />

<footnote n="4">⁴Due to their construction from data, such surrogate models are sometimes also called *data-fit* or *data-driven* surrogates (Robinson et al., 2008; Asher et al., 2015) to distinguish them from both reduced-order models, which</footnote>
As opposed to real-world experiments, deterministic computer simulations are not affected by noise. Still, instead of interpolating the training data, most surrogate modeling methods approximate it, i.e., they determine a function that minimizes a suitable distance to the data points under some assumptions and constraints. Depending on the functional form of the surrogate, approximation can have a smoothing effect, i.e., it might capture the low-frequency components of the data and disregard the high-frequency ones. It is well-known that polynomial interpolation can be unstable unless specific experimental designs are used, which is known as Runge's phenomenon (Xiu, 2010). More generally, interpolating the training data can lead to overfitting, i.e., the surrogate might lose the ability to generalize to unseen points, whereas approximation is often more stable.

The various families of surrogate models differ in their functional form and in the number of parameters they involve. Determining the exact functional form and the parameter values fitting the provided data best is often called training, in reference to machine learning terminology. As usual, there is a trade-off between flexibility of the model and requirements on the amount of data: a function with many free parameters and a flexible functional form is able to approximate a wide variety of models, but it requires a lot of training data to avoid overfitting. On the other hand, a function making many structural assumptions and involving only a few parameters can be trained based on a small set of data, but it might be inexact if the structural assumptions do not match the true shape of the model. When knowledge about model characteristics is available, it can be used to select a surrogate model with suitable structure.

We give an overview of common surrogate models for deterministic simulators in Section 2.2.1. Then, in Section 2.2.2 we highlight some connections between surrogate modeling and (supervised) machine learning techniques. Finally, we discuss in Section 2.2.3 how the performance of a surrogate model can be assessed.

## 2.2.1 Overview of common surrogate models

In this section we give a short overview of the families of surrogate models that are typically used in UQ for deterministic simulators.<fnref n="5" />

Polynomial chaos expansions (PCE) approximate the original model by a weighted sum of multivariate polynomials, which are orthonormal with respect to the density $f_X$ of the input. Under certain assumptions on $f_X$, the orthonormal polynomials constitute a basis for $L^2_{f_X}(\mathcal{D})$ defined in Eq. (2.1) and are therefore suitable for approximating square-integrable models. One of the advantages of PCE is that moments and certain kinds of sensitivity indices can be computed analytically from its functional form (see Section 2.6). PCE is one of the main surrogate modeling techniques used and investigated in this thesis, and is explained in further detail in Section 2.4.1.1.

<footnote>are simplified models utilizing knowledge about the physics and equations of the original model (Benner et al., 2015), and hierarchical (multifidelity) surrogates. However, as we explain in more detail in Section 2.2.2, the typical setup in UQ does not rely on data alone but assumes to have access to both the input distribution and the computational model from which the data is constructed (see also Section 2.1).</footnote>

<footnote n="5">⁵The use of surrogates is not restricted to the field of UQ. Surrogate models are used also for other purposes such as optimization (Robinson et al., 2008; Forrester and Keane, 2009; Razavi et al., 2012).</footnote>
Kriging, or Gaussian process modeling, is another popular surrogate technique in UQ. It has its origins in mining and geostatistics (Krige, 1951; Matheron, 1963) and was introduced for general computational models by Sacks et al. (1989). In the machine learning community, it is also known as Gaussian process regression (Rasmussen and Williams, 2006). In Kriging, we construct a Gaussian stochastic process that is conditioned on the given experimental design. Depending on the formulation, the resulting process can be interpolating or approximating. The choices of covariance function and associated hyperparameters govern the regularity of the realizations. The final prediction, i.e. the surrogate, is usually the mean of the inferred process. Being a stochastic process, Kriging provides confidence bounds for the prediction, which can be used for goal-oriented active learning, i.e., for the sequential enrichment of the experimental design with additional informative points. Kriging is popular because of its suitability for active learning and because its covariance function is a kernel, which is a powerful mathematical concept allowing to work even with non-numeric data. A drawback of Kriging is its unfavorable scaling with the number of experimental design points: to compute a Kriging metamodel, a dense covariance matrix needs to be inverted, which becomes prohibitively expensive for large experimental designs. Many of the approaches for surrogate modeling of stochastic simulators in the literature rely on Kriging in some way (see Section 2.3.3).

As a hybrid between the two types of approaches, polynomial chaos-Kriging or PC-Kriging has been developed (Kersaudy et al., 2015; Schöbi et al., 2015). Here, a Kriging metamodel is constructed whose trend function is a PCE, thus combining the global behavior of PCE with the local correction capabilities of Kriging.

Several other types of surrogate models have been proposed, including low-rank tensor approximations, also known as canonical decompositions (Nouy, 2010; Chevreuil et al., 2015; Konakli and Sudret, 2016), support vector machines (Smola and Schölkopf, 2004; Bourinet et al., 2011), regression trees (Breiman et al., 1984) and other localized regression methods (Marelli et al., 2021b), radial basis functions (Buhmann, 2000), and artificial neural networks (Hajela and Berke, 1992; Tripathy and Bilionis, 2018).

There are numerous papers reviewing and comparing surrogate modeling methods for different purposes (sensitivity analysis, optimization, etc.) and different fields of applications (Forrester and Keane, 2009; Razavi et al., 2012; Asher et al., 2015; Le Gratiet et al., 2017; Chatterjee et al., 2019; Becker, 2020; Cheng et al., 2020). Unsurprisingly, no surrogate modeling method is consistently superior to all others. The decision for a specific type of surrogate model depends on the characteristics of the computational model such as dimensionality and regularity, on the amount and properties of the available data, and on the purpose of the analysis (uncertainty propagation, reliability, sensitivity, optimization, etc.).

Since no single method performs best in all cases, several authors propose combining the predictions of several surrogates, or choosing one out of many based on model selection criteria (Viana et al., 2009; Sagi and Rokach, 2018; Liu et al., 2020b; Cheng and Lu, 2020; Parisi et al., 2021). This concept of aggregating surrogate models is known as ensemble modeling. For an overview of common model selection criteria, see Section 2.2.3.

Some surrogate models provide more than just an approximation to the model: for example, they might additionally yield analytical estimates of the response moments of the model, sensitivity
indices, or confidence bounds for their predictions. However, since evaluations of surrogate
models are cheap, it is always possible to compute moments numerically using Monte Carlo
simulation. Bootstrap can often be used to obtain confidence bounds (Hastie et al., 2001;
Marelli and Sudret, 2018).

2.2.2 Surrogate modeling versus supervised machine learning

Readers familiar with machine learning (Bishop, 2006) will have recognized many methods
and concepts mentioned in the previous section, which gives rise to the question: what is the
difference between surrogate modeling in the field of UQ and *supervised learning* in the field
of machine learning? The goal of supervised learning is, for a given set of input-output data
pairs, to learn an suitable map that is able to make predictions for unseen data points. If the
domain of the output data is a discrete set (e.g., when the outputs are categorical), the learning
task is commonly called *classification*, while in the continuous case it is called *regression*. In
this section, we summarize some of the relevant differences and similarities between black-
box surrogate modeling and supervised machine learning regression, and point out interesting
connections between the fields of uncertainty quantification and machine learning.

Mathematically, the setup of the problems looks similar. Both black-box surrogate modeling and supervised machine learning regression rely on a set of training data, and make a number of structural assumptions in order to approximate the input-output (or in machine learning terminology, feature-response) relationship using a mathematical model.

However, they differ in the data-generating process and the size of the training set. In UQ, we usually assume that both the computational model and a probabilistic model of its input are available (as explained in Section 2.1). Due to the cost of the computational model, only a few model evaluations can be afforded; however, it is possible to choose the experimental design points in order to collect especially informative model evaluations. The evaluations of the deterministic computational model are not affected by noise. To compensate for the relatively small data sets, UQ methods typically make rather strict structural assumptions.

In contrast, in machine learning the data-generating process is usually not available, but only
the data itself is given, and it might be noisy. The amount of data is typically large (a situation
described with the well-known term “big data”), which permits the use of very expressive models
such as neural networks that have many free parameters. A challenge in machine learning is the
choice of network architecture as well as the training process, both of which involve numerous
hyperparameters that need to be chosen carefully, an operation known as *tuning*.

Machine learning methods typically need to deal with thousands (if not millions) of features, for
example, when the features are the pixels of an image in a classification task. Also, the features
are often not real numbers, but might be categorical variables, images or even text or speech.
While many classical UQ techniques were initially developed for rather low-dimensional numer-
ical input in the order of a few tens of continuous variables, recent UQ research is increasingly
concerned with input spaces of very high dimension (Lataniotis et al., 2020; Kontolati et al.,
2022b) and with the treatment of categorical or discrete inputs (Lauvernet and Helbert, 2020;
Roustant et al., 2020b; Zhu and Sudret, 2021a; Rohmer et al., 2022).
Some methods, like Gaussian process regression (Kriging) (Rasmussen and Williams, 2006), support vector machines (Vapnik, 1995) and principal component analysis (Pearson, 1901; Jolliffe, 2002) have been used for years by both the UQ and the machine learning communities. Artificial neural networks have been used for structural engineering and for uncertainty quantification (Hajela and Berke, 1992; Tripathy and Bilionis, 2018), while Torre et al. (2019) and Timpe et al. (2020) have explored the use of PCE for machine learning applications. As a surrogate modeling method that is not purely black-box, the so-called physics-informed neural networks have received a lot of attention in recent years for the solution of (stochastic) differential equations (Raissi et al., 2019; Zhu et al., 2019; Yang et al., 2020; Karniadakis et al., 2021). Here, the solution is represented by a neural network, and the differential equations are embedded into the loss function used for training.

The connections between UQ and machine learning go beyond the data-driven identification of input-output maps. While neural networks are under certain assumptions universal approximators (Cheng and Titterington, 1994) and perform well for many practical problems when enough data is available, in their pure form they also have some disadvantages: it is often not obvious to understand why a neural network reaches a certain conclusion or what factors influence its response, nor how confident the network is in its prediction, all questions that by definition are central to the field of UQ. Therefore, recent machine learning research is increasingly concerned with questions of fairness, explainability and confidence (Sundararajan et al., 2017; Abdar et al., 2021; Bénesse et al., 2022) and also uses traditional UQ techniques for this purpose (Liu et al., 2020a; Fel et al., 2021).

### 2.2.3 Surrogate model assessment and model selection

In order to be useful, a surrogate model must accurately approximate the input-output relationship of the computational model. Therefore, we need to be able to evaluate how well a given surrogate is performing. Furthermore, surrogate models often depend on a number of tunable hyperparameters that need to be chosen according to such a performance criterion during training.

A global measure of the approximation quality is the mean-squared error between the model $\mathcal{M}$ and the surrogate model $\tilde{\mathcal{M}}$ defined by

$$ \epsilon = \mathbb{E}_{\mathbf{X}} \left[ \left( \mathcal{M}(\mathbf{X}) - \tilde{\mathcal{M}}(\mathbf{X}) \right)^2 \right]. \qquad (2.2) $$

Here the expectation is with respect to $\mathbf{X}$, so that the squared difference is weighted by the input density $f_{\mathbf{X}}$.

Denote by $\mathcal{X} = \{\mathbf{x}^{(i)}\}_{i=1}^N \subset \mathcal{D}$ the experimental design used to train the surrogate model $\tilde{\mathcal{M}}$, usually sampled from $f_{\mathbf{X}}$. A simple approximation to $\epsilon$ is the empirical error

$$ \epsilon_{\text{emp}} = \frac{1}{N} \sum_{i=1}^{N} (\mathcal{M}(\mathbf{x}^{(i)}) - \tilde{\mathcal{M}}(\mathbf{x}^{(i)}))^2. \qquad (2.3) $$

However, judging the approximation quality of the surrogate from the empirical error, i.e., from its performance on the same points that were used to train it, will in general underestimate the true prediction error.
A better approximation to $\epsilon$ is the validation error

$$ \epsilon_{\text{val}} = \frac{1}{N_{\text{val}}} \sum_{i=1}^{N_{\text{val}}} (\mathcal{M}(\boldsymbol{x}_{\text{val}}^{(i)}) - \tilde{\mathcal{M}}(\boldsymbol{x}_{\text{val}}^{(i)}))^2, \quad (2.4) $$

where $\mathcal{X}_{\text{val}} = \{\boldsymbol{x}_{\text{val}}^{(i)}\}_{i=1}^{N_{\text{val}}}$ is a large set of points sampled from $f_X$ independently of $\mathcal{X}$.

If such a large set of additional model evaluations is not available, which is usually the case for all but toy problems, the error can be approximated by cross-validation (Hastie et al., 2001) based on $\mathcal{X}$ as follows: the experimental design is split into several parts, called folds. Each fold is in turn treated as validation set, while the remaining points of the experimental design are used to construct the surrogate model. Finally, the errors computed on the left-out sets are averaged, and serve as an estimate for the validation error. If the number of folds is equal to the size of the experimental design, the method is known as leave-one-out (LOO) cross-validation (Hastie et al., 2001).

Cross-validation is a general error assessment technique that is not only used to judge the performance of the final surrogate, but is also often used already during the training process for selecting hyperparameters, or – in the case of iterative methods – for deciding when to stop the iterations. These decisions are also called *model selection*, because we have to decide between different surrogate models (of the same type) arising from different choices of hyperparameters.

For methods which aim at fitting a random quantity to given data, such as Kriging which infers a stochastic process (Rasmussen and Williams, 2006), or weak polynomial chaos approximation which fits a distribution (Xiu, 2010), an alternative discrepancy measure is the *likelihood*, which is the probability of observing the training data given a certain choice of parameters. A model with many free parameters is typically more flexible and can achieve a larger likelihood than one with fewer parameters. Therefore, model selection based on the likelihood alone can lead to overfitting. To prevent this, model selection criteria such as Akaike's information criterion (AIC) (Akaike, 1974) and the Bayesian information criterion (BIC) (Schwarz, 1978) additionally take the number of free parameters into account. It has been shown that AIC is asymptotically equivalent to cross-validation (Stone, 1977).

It should be noted that the repeated application of the same criterion for training and selection among a large number of models can also lead to overfitting (Rao et al., 2008; Arlot and Celisse, 2010; Cawley and Talbot, 2010).

## 2.3 Surrogate modeling for stochastic simulators

The classical framework of UQ shown in Figure 2.1 assumes that all stochasticity influencing the system can be modeled probabilistically in the form of random input variables. The computational model itself is then deterministic and can be surrogated by the deterministic function approximation methods described in Section 2.2. However, there are cases in which it is not desired or even not possible to completely separate the stochasticity from the computational model. The model contains residual randomness which causes the response to be a random quantity, even for fixed input parameters. In this case, we call the model a *stochastic simulator*. There are two main cases from which stochastic simulators arise.
### 2.3.1 Examples of stochastic simulators

In some instances, the stochasticity is an intrinsic part of the computational model. For example, consider the stochastic susceptible-infected-recovered (SIR) model developed in epidemiology to study the spread of an infectious disease (Britton, 2010). Here, the population consists of individuals that are either susceptible to a disease, infected by it, or recovered and therefore immune. Two numbers, the rates of infection and recovery, govern the stochastic evolution of the disease with respect to time. Due to the stochasticity of infection and recovery, the evolution might look different each time the simulation is run, even if the population size and the rate parameters are kept at the same values. More generally, any stochastic differential equation with a random forcing term is a stochastic simulator of this class. In these models, the stochasticity drives the simulation and therefore cannot be modeled separately: the model is intrinsically a stochastic simulator.

Stochastic simulators also arise in cases when the complete characterization of the stochasticity would in principle be possible but is not desired, e.g., if the latter is too high-dimensional or too complicated to be modeled explicitly. For example, when designing a wind turbine for a new geographical site, the structural safety of the turbine over a period of many years must be assessed for the specific wind conditions at the respective location. This is done by employing a wind generation model together with a deterministic aero-servo-elastic simulation of the wind turbine structure (Dimitrov et al., 2018; Slot et al., 2020).

The wind generation model generates wind boxes, i.e., realizations of a vector-valued random field defined in both time and space. A wind box is often summarized by a few statistics such as mean wind speed, turbulence, inflow angle etc., which together are called wind climate. For each wind climate, an arbitrary number of wind boxes can be generated that match these values. If we completely separated the wind box generation from the multi-physics simulation of the turbine, the latter would be a deterministic simulator; however, the wind fields are too high-dimensional to be modeled explicitly using the classical methodology described in Section 2.2.<fnref n="6" /> Therefore, it can be desirable to treat the generation of wind boxes from wind conditions as a part of the computer model, resulting in a stochastic simulator which, for given wind climate, returns a different response each time the simulation is run. The same reasoning shows that the simulation of earthquake damage to a structure can be regarded as a stochastic simulator (Zhu et al., 2022). In these cases, the stochasticity comes from latent variables that are not explicitly modeled probabilistically as input variables.

### 2.3.2 Stochastic simulators as random fields

As explained in Section 2.1, stochastic simulators can be seen as random fields for which the space of explicitly modeled parameters $\mathcal{D}$ assumes the role of the index set, and the latent stochasticity is represented as an event from an abstract space $\Omega$.

<footnote n="6">⁶A conceptually quite different approach, namely autoregressive models, has recently been proposed for constructing custom surrogates for such high-dimensional time-dependent random input (Dimitrov et al., 2022).</footnote>
The random field $\mathcal{M}: \mathcal{D} \times \Omega \to \mathbb{R}$ has the properties that

$$ - \text{ For each } \boldsymbol{x} \in \mathcal{D}, Y_x = \mathcal{M}(\boldsymbol{x}, \cdot) : \Omega \to \mathbb{R} \text{ is a random variable;} \quad (2.5) $$

$$ - \text{ For each } \omega \in \Omega, \mathcal{M}(\cdot, \omega) : \mathcal{D} \to \mathbb{R} \text{ is a deterministic function, called trajectory.} \quad (2.6) $$

We call Eq. (2.5) the random variable view of the random field, and Eq. (2.6) the random function view. These two perspectives lead to different types of surrogate models, as we will see below.

We assume that $\mathcal{M} \in L_{f_x}^2(\mathcal{D}) \times L^2(\Omega, \mathcal{F}, P)$ where $(\Omega, \mathcal{F}, P)$ is the probability space of the stochastic event $\omega \in \Omega$.

Just like deterministic computer models, stochastic simulators can only be evaluated pointwise, i.e., for a finite number of experimental design points $\boldsymbol{x}^{(i)} \in \mathcal{D}$ and stochastic events $\omega^{(j)} \in \Omega$. Repeated evaluations for the same location $\boldsymbol{x} \in \mathcal{X}$, each time using a different $\omega^{(j)}, j = 1, \dots, R$, are called replications. For many computational models, it is possible to fix the sequence of random numbers used for the computations by fixing the seed of the random number generator, which allows the (pointwise) evaluation of the simulator along trajectories even if the underlying stochastic event $\omega$ is not known.

To fully characterize a general random field, it is necessary to specify the collection of all its finite-dimensional distributions (Grigoriu, 2002, Chapter 3.5)

$$ \{F_{Y_{x_1}, \dots, Y_{x_n}} : n \ge 1, \boldsymbol{x}_1, \dots, \boldsymbol{x}_n \in \mathcal{D}\}, \quad (2.7) $$

where $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)$. The univariate distributions $F_{Y_x}$ are called marginals or marginal distributions. Extending the concept of moments of random variables to random fields, the deterministic mean function of $\mathcal{M}$ is given by $\mu(\boldsymbol{x}) = \mathbb{E}[Y_x]$. If $\mu(\boldsymbol{x}) = 0$, the random field is called centered. The (auto-)covariance function of $\mathcal{M}$ is defined by

$$ C(\boldsymbol{x}, \boldsymbol{x}') = \text{Cov}[Y_x, Y_{x'}]. \quad (2.8) $$

To enable the theory presented in Section 2.4.2 below, we assume that $C$ is continuous on $\mathcal{D} \times \mathcal{D}$.

While Eq. (2.7) is needed for general random fields, there are a few special cases that can be characterized by a smaller number of properties, most notably the family of Gaussian processes, for which all finite-dimensional joint distributions are multivariate Gaussian distributions. A Gaussian process is uniquely defined by its mean and covariance function. Its conditional distributions are again multivariate Gaussians, which motivates the popular surrogate modeling technique Kriging a.k.a. Gaussian process modeling (Section 2.2.1). While Gaussian random fields are computationally convenient, random fields encountered in real-world problems (and in particular, stochastic simulators) are often not Gaussian. One obvious argument is that Gaussian variables are unbounded, whereas physical quantities are almost always bounded (Grigoriu, 2002, p. 118).

### 2.3.3 Overview of surrogate modeling approaches

Most surrogate modeling approaches for stochastic simulators developed to date have focused on the random variable view (Eq. (2.5)), i.e., on their interpretation as a collection of random
variables indexed by the set of inputs. The goal of these approaches, which we describe in Section 2.3.3.1, is to approximate the random variables $\mathcal{M}(\boldsymbol{x}, \cdot)$ in distribution, usually relying on the regularity of the mapping with respect to $\boldsymbol{x}$. In other words, they approximate the marginal distributions $F_{Y_{\boldsymbol{x}}}$, without attempting to match any of the higher-order joint distributions (i.e., $n > 1$ in Eq. (2.7)). The resulting surrogate can be used to generate pointwise samples, which will approximately follow the correct marginal distributions.

By construction, this class of surrogates is in general not able to generate trajectories. To obtain a surrogate model that has this feature, we must take the random function view (Eq. (2.6)), and data sampled from trajectories must be available. In the context of stochastic simulators, the literature on such surrogate models is scarce. However, in the broader context of inferring general random fields from given samples, there are a number of methods, which we discuss in Section 2.3.3.2.

### 2.3.3.1 Surrogates in distribution

We first give an overview of surrogate models that characterize the response random variable in distribution. This can be seen as conditional distribution estimation, aiming at characterizing the quantity $P(Y | \boldsymbol{X} = \boldsymbol{x})$.

To obtain an explicit expression for the statistics or the full distribution of $\mathcal{M}(\boldsymbol{x}, \cdot)$ at each location in the input domain $\mathcal{D} \ni \boldsymbol{x}$, various parametric approaches have been developed, which often rely on deterministic surrogate models to emulate these statistics or other auxiliary quantities (such as expansion coefficients) over $\mathcal{D}$. (Joint) Gaussian process models have been proposed to address scalar statistics of the random response variables, such as mean and variance (Iooss and Ribatet, 2009; Ankenman et al., 2010; Marrel et al., 2012; Binois et al., 2018). Other approaches emulate quantiles using quantile Kriging (Plumlee and Tuo, 2014) or other quantile regression methods (Torossian et al., 2020). Assuming that the full marginal distributions at a set of experimental design points $\boldsymbol{X} \subset \mathcal{D}$ are available, Moutoussamy et al. (2015) propose to represent the probability density in terms of a suitable functional basis. Similarly, Browne et al. (2016) propose to represent the quantile function of the stochastic simulator output in terms of a suitable basis, and emulate the parameters throughout the domain using Gaussian processes.

Another parametric choice is the generalized lambda model (Zhu and Sudret, 2020) that utilizes the generalized lambda distribution (Freimer et al., 1988), a parametric family of distributions with only four parameters that is able to closely approximate many of the classical unimodal parametric families (Gaussian, uniform, exponential, lognormal etc.). Zhu and Sudret (2020) apply this distribution to emulate the response densities based on replications, and extrapolate to the whole domain using PCE, using a joint computation scheme that simultaneously fits the response densities and the PCE models to the given data. The generalized lambda model has further been refined to avoid the need for replications by employing maximum conditional likelihood estimation (Zhu and Sudret, 2021b,c).

Replacing the generalized lambda distribution by a PCE approximation of the conditional random variables, Zhu and Sudret (2022) recently developed the replication-free stochastic PCE method. This approach uses PCE with two additional variables (a latent and a noise variable)
to emulate the stochastic simulator, inspired by the weak PCE methodology based on maximum likelihood estimation (Xiu, 2010).

A classical nonparametric method for conditional density estimation, likewise not requiring replications, is the well-known kernel density estimation (Fan and Gijbels, 1996; Hall et al., 2004).

Finally, in the field of machine learning, generative models like variational autoencoders (Kingma and Welling, 2014) and generative adversarial networks (Goodfellow et al., 2014) are designed to learn a target distribution from data and to sample from it, thus being a kind of surrogate in distribution. Similarly, Bayesian neural networks (MacKay, 1992; Goan and Fookes, 2020), whose weights are independent Gaussian random variables, are able to approximate conditional distributions. Bayesian methods such as Markov Chain Monte Carlo or variational inference are used to infer the parameters of the weight densities from the given data.

### 2.3.3.2 Surrogates using the random function view

Whereas the methods presented in the previous section are able to accurately characterize the marginal response distributions, they do not model the joint distributions of higher order. Of course, unless we are in the specific case of a Gaussian random field, modeling all the finite-dimensional distributions in Eq. (2.7) is infeasible. If we still want to capture some of the higher-order information about the random field, we need to make use of its structure, i.e., of the properties of its covariance function (Eq. (2.8)) and the induced regularity of its trajectories. In other words, we need to utilize the random function view of random fields shown in Eq. (2.6). The overall goal is to find a representation of the random field that can be sampled from, both pointwise and trajectory-wise.

In the literature, two classes of methods for the representation of random fields can be distinguished. The first is the construction of a representation for given properties such as marginal distributions, covariance function, or spectral density.<fnref n="7" /> These properties are assumed to be known analytically, and the goal is to generate samples from a random field with these properties. Note that general non-Gaussian random fields are not completely characterized by their marginals and covariance function (Grigoriu, 2002). We will give a short overview of such methods in the following.

The second class, which is most relevant in the context of surrogating stochastic simulators, is the inference of a representation from data, i.e., from pointwise samples of the random field along trajectories. This data is generated by a real-world process (such as an earthquake) or by a stochastic computational model, both of which behave like a random field. However, it might be impossible or too expensive to use these processes or models directly for generating a large amount of samples, and therefore a surrogate needs to be constructed. These methods are discussed below.

<footnote n="7">⁷The spectral density exists only for weakly stationary random processes, i.e., for processes with constant mean and an autocorrelation function $r(\tau) = \mathbb{E}[Y_x Y_{x+\tau}]$ that only depends on the difference between two locations. For a one-dimensional centered process, it is defined by $s(\nu) = \frac{1}{2\pi} \int_{-\infty}^{\infty} r(\tau) \exp(-i\nu\tau) d\tau$, i.e., it is the Fourier transform of the autocorrelation function (Grigoriu, 2002).</footnote>
**Simulation of non-Gaussian random fields from properties.** A simple way to discretize a random field is to represent it by a random vector (i.e., a finite number of random variables) on a mesh defined in the index set. The value of the random field at other locations is then computed from the random vector. Sudret (2007) gives an overview of such methods, including the midpoint method (Der Kiureghian and Ke, 1988) and spatial averaging (Vanmarcke and Grigoriu, 1983). However, Sudret and Der Kiureghian (2000) point out that these methods are rather inefficient.

Alternatively, series expansion methods exploit the structure of the random field in order to represent it using “spatial” functions (i.e., functions defined over the index set) and random coefficients (Sudret, 2007). Grigoriu (2010) classifies the series expansion methods into linear models with dependent parameters and nonlinear models, also known as translation processes.

Linear series expansion methods represent the random field as a series of spatial functions $h_k(\boldsymbol{x})$ weighted by random coefficients $\xi_k(\omega)$, with the general form (Grigoriu, 2006, 2010)

$$ \mathcal{M}(\boldsymbol{x}, \omega) = \mu(\boldsymbol{x}) + \sum_{k=0}^{\infty} \xi_k(\omega)h_k(\boldsymbol{x}). \qquad (2.9) $$

The spectral representation method expands a weakly stationary random process in terms of trigonometric polynomials (Fourier expansion) (Deodatis and Shinozuka, 1991; Grigoriu, 2006).<fnref n="8" /> Orthogonal series expansion (OSE) expands a random field onto any orthonormal basis for the spatial component (Zhang and Ellingwood, 1994). Finally, Karhunen-Loève expansion (KLE) expands a general random field onto an optimal orthonormal spatial basis, which is computed as the solution of an integral eigenvalue problem involving the covariance function (Karhunen, 1946; Loève, 1978). We will explain KLE in detail in Section 2.4.2. In all these cases, the random coefficients $\xi_k$ of the resulting expansion are in general dependent, and failing to account for this dependence can result in improper characterizations of the random field (Grigoriu, 2010).

Note that the linear series expansion in Eq. (2.9) transforms the random field, which (taking the random variable view of Eq. (2.5)) is an uncountably infinite collection of correlated random variables, into a countable collection of random variables $\{\xi_k\}_{k=0}^{\infty}$ together with spatial functions $\{h_k\}_{k=0}^{\infty}$ capturing the correlation.

For the special case of Gaussian random fields, the expansion optimal linear estimation (EOLE) method (Li and Der Kiureghian, 1993) provides a compromise between pointwise and series representation, by representing the random field by a random vector defined on a mesh, and interpolating between the mesh points using optimal linear estimation, relying on the Gaussian assumption just like Kriging. The random vector can be further compressed by using a form of dimension reduction closely related to KLE (see Section 2.4.2).

Translation processes make use of the fact that simulating a Gaussian process G is comparatively easy, e.g. by EOLE or KLE, since such a process is fully characterized by its mean and covariance function, and all conditional distributions are Gaussian (Grigoriu, 2006). For given marginals and covariance function or spectral density of the general non-Gaussian random field of interest,

<footnote n="8">⁸Strictly speaking, this interpretation only holds for periodic processes. For general processes, the analysis is slightly more involved (Grigoriu, 2006).</footnote>
the goal of translation process algorithms is to find a mapping $\mathcal{T}$ such that the induced ran-
dom field $\mathcal{T}(G)$ has the requested properties (Grigoriu, 1984; Yamazaki and Shinozuka, 1988;
Grigoriu, 1998; Deodatis and Micaletti, 2001; Sakamoto and Ghanem, 2002; Poirion and Puig,
2010; Shields et al., 2011; Kim and Shields, 2015). Here, the difficulty is that the marginals and
covariance function might be inconsistent, i.e., there might be no such mapping $\mathcal{T}$ that achieves
the desired properties (Grigoriu, 2010).

**Inference of random fields from data.** The methods described above all assume that the marginals and the covariance function or the spectral density of the random field of interest are known analytically. Now, we turn to methods that aim to characterize a random field solely based on a set of samples. Many of these methods rely on the spectral method *Karhunen-Loève expansion* (KLE), which is described in detail in Section 2.4.2 below.

One popular approach is based on the repeated evaluation of the random field at the points of a given mesh. In other words, the random field is reduced to a high-dimensional random vector with one component per mesh point, similar to the pointwise representation described above. First, discrete KLE (i.e., principal component analysis (PCA); see Section 2.4.2) is performed on the realizations of the high-dimensional random vector to reduce its dimensionality. Then, the random coefficients are modeled, often by PCE (Ghanem and Spanos, 1991; Soize, 2010). The PCE coefficients are determined, e.g., by maximum likelihood (weak PCE) (Desceliers et al., 2006; Perrin et al., 2012) or Bayesian inference (Ghanem and Doostan, 2006). Das et al. (2009) propose two methods to approximate the experimentally determined multivariate joint distributions of the random coefficients by PCE based on conditional distributions and the Rosenblatt transform (Rosenblatt, 1952). To further model the uncertainty in the stochastic model induced by the finite size of the sample set, a number of papers explores the inference of random instead of deterministic PCE coefficients (Ghanem and Doostan, 2006; Das et al., 2008; Arnst et al., 2010; Soize, 2010).

The extension of the inferred high-dimensional random vector (defined on the mesh) to the
continuous index set (where the original random field is defined) is often not discussed, because
it is considered sufficient to know the random field on a fine mesh. Among the few publications
explicitly addressing this topic, Desceliers et al. (2006) suggest to use the interpolation functions
of the finite element method. Das et al. (2009) generally refer to this aspect as “a task of
interpolation or/and approximation technique”. In the specific context of stochastic simulators,
Azzi et al. (2019) propose to use PCE or Kriging for this extension.

A variant of the KLE-PCE methodology described above is the intrusive method proposed by
Doostan et al. (2007), which was later extended to the non-intrusive case by Raisee et al. (2015)
and Abraham et al. (2018). Starting from a PCE expansion for each of the random variables
of a coarse mesh computed from given replications, this method first computes a KLE on the
coarse mesh (“optimal expansion”). The resulting random variables of the KLE, modeled as
PCE, are then used to derive spatial functions on a finer mesh.

In the context of earthquake engineering, in a series of papers Zentner and Poirion (2012); Poirion and Zentner (2013, 2014) reverse this construction: first, they interpolate the trajectories based on the given data and express them in terms of an orthonormal basis. Then they perform KLE
on the interpolated trajectories, which – due to the orthonormality of the basis – reduces to a
discrete eigenvalue problem. The random coefficients are modeled by kernel density estimation,
taking bivariate dependence into account (Poirion and Zentner, 2014).

Azzi et al. (2019) propose another approach based on KLE, which first infers a continuous
representation of the covariance function based on the sample covariance matrix computed from
replications at mesh points. Then, the eigenfunctions of KLE are computed on a finer mesh.

Because KLE is linear in the random coefficients (see Eq. (2.9)), the complexity of non-Gaussian random fields potentially requires the random coefficients to follow intricate joint distributions (see also Section 2.4.2). To avoid having to infer complicated joint distributions, Sarma et al. (2008) and Ma and Zabaras (2011) instead suggest using kernel PCA, a nonlinear extension of KLE. The data is first mapped to the so-called feature space using a nonlinear transformation, and then KLE is applied, with the hope that the data is (close to) linear in the feature space. A difficulty of kernel PCA is the pre-image problem: the obtained realizations in the feature space must be mapped back to the original space of the data. Sarma et al. (2008) use a polynomial kernel while Ma and Zabaras (2011) use a Gaussian kernel. Both assume that in the feature space, the random coefficients of the linear expansion (Eq. (2.9)) are independent.

In the context of stochastic differential equations with stochastic forcing as well as random input parameters, Navarro Jimenez et al. (2017) suggest to model the resulting random field by PCE with random coefficients (i.e., choosing PCE basis functions for $h_k$ in Eq. (2.9)). The stochastic coefficients are determined by non-intrusive pseudospectral projection. They do not discuss resampling the random field, because the purpose of their construction is sensitivity analysis.

## 2.4 Spectral methods

The topic of this thesis is surrogate modeling with spectral methods. The word *spectral* usually hints to the connection to an eigenvalue problem, since the *spectrum* of an operator denotes the set of its eigenvalues. Solving an eigenvalue problem yields, in addition to the eigenvalues, eigenvectors or eigenfunctions that are orthogonal in the respective space where the problem is set. Under appropriate conditions, the eigenvectors or -functions constitute a basis for this space.

A *spectral expansion* is then the representation of an element from that space in terms of such an orthogonal basis. More generally, this expression is used for any representation as long as the respective basis is orthogonal.

### 2.4.1 Chaos expansions

In UQ, a spectral expansion in $L_{fx}^2(\mathcal{D})$ is also often called *chaos expansion*, in reference to an early contribution by Wiener (1938) who introduced them for stochastic processes using Hermite polynomials.<fnref n="9" />

<footnote n="9">⁹Therefore, the use of the term *chaos expansion* in UQ predates the mathematical field of chaos theory by several decades (Xiu, 2010). Chaos expansions as introduced here do not belong to the field of chaos theory.</footnote>
The general mathematical framework of chaos expansions is as follows. Together with the inner product

$$ \langle g, h \rangle = \int_{\mathcal{D}} g(\boldsymbol{x})h(\boldsymbol{x})f_{\boldsymbol{X}}(\boldsymbol{x}) \, d\boldsymbol{x}, \qquad (2.10) $$

the space of square-integrable functions $L_{f_X}^2(\mathcal{D})$ defined in Eq. (2.1) is a Hilbert space.<fnref n="10" /> Interpreted probabilistically, a model $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ has finite variance under $\boldsymbol{X}$, i.e., $\mathrm{Var}_{\boldsymbol{X}}[\mathcal{M}(\boldsymbol{X})] < \infty$.

Let $\{\psi_k\}_{k=0}^{\infty}$ be an orthonormal basis for $L_{f_X}^2(\mathcal{D})$, i.e.,

1. for all $i, j$, $\langle \psi_i, \psi_j \rangle = \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{otherwise} \end{cases}$; and

2. $\{\psi_k\}_{k=0}^{\infty}$ is dense in $L_{f_X}^2(\mathcal{D})$.

Then the chaos expansion of $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ in terms of this basis is given by

$$ \mathcal{M}(\boldsymbol{x}) = \sum_{k=0}^{\infty} c_k \psi_k(\boldsymbol{x}), \qquad (2.11) $$

where

$$ c_k = \langle \mathcal{M}, \psi_k \rangle \qquad (2.12) $$

is the projection of $\mathcal{M}$ onto $\psi_k$.

It follows that the random response vector $Y = \mathcal{M}(\boldsymbol{X})$ can be expressed as a series of random variables $Z_k = \psi_k(\boldsymbol{X})$ by

$$ Y = \sum_{k=0}^{\infty} c_k Z_k. \qquad (2.13) $$

From now on, we assume that the components of the random vector $\boldsymbol{X}$ are independent.<fnref n="11" /> Therefore, the density $f_X$ factorizes into $f_X(\boldsymbol{x}) = \prod_{i=1}^d f_{X_i}(x_i)$. A basis of $L_{f_X}^2(\mathcal{D}) = L_{f_{X_1}}^2(\mathcal{D}_1) \otimes \dots \otimes L_{f_{X_d}}^2(\mathcal{D}_d)$ can be obtained from one-dimensional bases $\{\psi_k^{(i)}\}_{k=0}^{\infty}$ of $L_{f_{X_i}}^2(\mathcal{D}_i)$ by a standard tensor product construction: an element $\psi_\alpha$ of the multivariate basis is given by

$$ \psi_{\alpha}(\boldsymbol{x}) = \prod_{i=1}^{d} \psi_{\alpha_i}^{(i)}(x_i), \qquad (2.14) $$

where the so-called multi-index $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{N}_0^d$ collects the indices of the respective univariate basis elements.

We will describe two particular types of chaos expansion below, namely polynomial chaos expansions (PCE) in Section 2.4.1.1 and Poincaré chaos expansions in Section 2.4.1.2. Beyond these, a chaos expansion in terms of a basis consisting of B-splines has recently been proposed (Rahman, 2020; Eckert et al., 2020). While rarely used in UQ, another instance that fits into the chaos expansion framework is the well-known Fourier expansion, whose orthonormal basis consists of sines and cosines.

<footnote n="10">¹⁰A Hilbert space is a complete metric space with an inner product that induces the metric. It is a powerful concept in functional analysis and applied mathematics. Under certain conditions it admits a countable orthonormal basis, which is essential for function approximation.</footnote>

<footnote n="11">¹¹We discuss the case of dependent input variables in Section 2.4.1.1 for PCE and in Section 2.6.3 for sensitivity analysis.</footnote>
2.4.1.1 Polynomial chaos expansion

Polynomial chaos expansion denotes the expansion of a model $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ in terms of a multi-variate orthonormal polynomial basis obtained as the tensor product of univariate orthonormal polynomial bases (Eq. (2.14)). In the following, we describe the construction of the latter.

Let $f_X$ be the probability density function of a random variable $X$ defined on $\mathcal{D} \subset \mathbb{R}$. An orthonormal polynomial basis for $L_{f_X}^2(\mathcal{D})$ exists if and only if the distribution function $F_X$ is uniquely defined by the sequence of its moments (Ernst et al., 2012). This is the case for most of the standard parametric distributions such as uniform, Gaussian, Beta or Gamma distributions, but not for the lognormal distribution. A sufficient condition for existence is that all moments are finite and the support of the distribution is compact (Ernst et al., 2012).

For the classical parametric distribution families, the associated orthonormal polynomials are well-known, and their properties have been studied for centuries (Szegö, 1939; Gautschi, 2004). When $X$ is Gaussian, the associated orthonormal polynomials are the (probabilists’) *Hermite polynomials*. A uniformly distributed $X$ is associated to the *Legendre polynomials*. A Beta distribution generates the *Jacobi polynomials*, and a Gamma distribution the *Laguerre polynomials*. In Figure 2.2, we illustrate the first five orthonormal polynomials associated to two different densities.

<figure><img src="image_5.png" /><figcaption>**Figure 2.2:** The first five orthonormal polynomials associated to two different densities: Left, Hermite polynomials associated to the Gaussian distribution. On the right, arbitrary polynomials associated to a custom (non-classical) multimodal density. The basis functions are plotted in color, while the densities are shown in gray.</figcaption></figure>

Originally, PCE was proposed in terms of Hermite polynomials by Ghanem and Spanos (1991) corresponding to the Gaussian distribution. The approach was generalized by Xiu and Karniadakis (2002) to distributions from the Askey scheme, i.e., the parametric families mentioned above. However, PCE can even be computed for arbitrary distributions (including discrete ones), as long as they admit an orthonormal polynomial basis, which is known under the name arbitrary PCE (Soize and Ghanem, 2004; Wan and Karniadakis, 2006; Oladyshkin and Nowak, 2012; Torre et al., 2019).
To construct an orthonormal polynomial basis for a given distribution, there are several options. For the three classical polynomial families Hermite, Jacobi (which includes Legendre as special case), and Laguerre, one way is to solve the associated eigenvalue problem which involves a differential operator (Bochner, 1929). An alternative for any type of distribution is orthonormalization of the set of monomials {1, x, x², x³, ...} using the Gram-Schmidt algorithm.

A third, often more stable, method is the *Stieltjes procedure* (Gautschi, 1982, 2004): Orthonormal polynomials obey an equation known as the three-term recurrence rule

$$ \sqrt{\beta_{n+1}} \psi_{n+1}(x) = (x - \alpha_n) \psi_n(x) - \sqrt{\beta_n} \psi_{n-1}(x), \quad (2.15) $$

starting with $\psi_{-1}(x) = 0$ and $\psi_0(x) = 1$, where $\alpha_n$ and $\beta_n$ are given by

$$ \alpha_n = \frac{\langle x \psi_n, \psi_n \rangle}{\langle \psi_n, \psi_n \rangle}, \quad (2.16) $$

$$ \beta_n = \frac{\langle \psi_n, \psi_n \rangle}{\langle \psi_{n-1}, \psi_{n-1} \rangle}. \quad (2.17) $$

The coefficients $\alpha_n$ and $\beta_n$ are known analytically for the classical orthonormal families. For other cases, they can be computed by numerical integration.

A univariate orthonormal polynomial basis constructed by Eq. (2.15) has exactly one element of each degree $k \in \mathbb{N}_0$. The degree is used as the index of the elements in the basis.

If an orthonormal polynomial basis does not exist, such as in the case of the lognormal distribution, we can first map the given input to another distribution and then use the associated polynomials for the PCE. Such a mapping is called *isoprobabilistic transform*. For example, in the case of a lognormal variable $X$ a transformation to Gaussian variables is easily achieved by taking the logarithm: $Z = \mathcal{T}(X) = \log(X)$ is Gaussian. The PCE is then computed for the transformed model $\tilde{\mathcal{M}}(Z) = \mathcal{M}(\mathcal{T}^{-1}(Z))$.

In principle, this procedure can be used for any two continuous random variables $X$ and $Z$ using the probability transform $\mathcal{T} = F_Z^{-1} \circ F_X$. However, especially when one of the variables is bounded and the other unbounded, the transform can introduce significant nonlinearity into the composite model $\tilde{\mathcal{M}} = \mathcal{M} \circ \mathcal{T}^{-1}$, which in turn requires high-degree polynomials to be approximated well (Lei et al., 2019; Torre et al., 2019).

The case of dependent input variables is an active area of research. If the components of $X$ are not independent, it is still possible to construct a multivariate polynomial basis that is orthonormal with respect to the distribution of $X$, e.g., by applying the Gram-Schmidt procedure (Torre et al., 2019; Jakeman et al., 2019; Lei et al., 2019); however, it is not unique but depends on the ordering of terms. Alternatively, the dependent random vector $X$ can be transformed to an independent one using a suitable probability transform (Lebrun and Dutfoy, 2009) before computing the PCE, which likewise is a non-unique operation depending on the ordering of the variables. Torre et al. (2019) demonstrated that for the purpose of global approximation, it can be advantageous to use a basis that is orthonormal to the marginals only, ignoring the dependence structure.
### 2.4.1.2 Poincaré chaos expansion

The Poincaré chaos expansion has recently been proposed by Roustant et al. (2020a), expanding on the results of an earlier paper by Roustant et al. (2017) with the goal of computing Poincaré constants for univariate distributions.

In one dimension, the Poincaré constant $C_P(f_X)$ associated to the probability density $f_X$ of a random variable $X$ is the smallest possible constant $C$ with

$$ \mathbb{E}_X[h(X)^2] \leq C \mathbb{E}_X[(h'(X))^2] \quad (2.18) $$

for all $h \in H_{f_X}^1$ with $\mathbb{E}_X[h(X)] = 0$, where $H_{f_X}^1$ is the Sobolev space of weakly differentiable functions, i.e., $H_{f_X}^1 = \{h \in L_{f_X}^2 : h' \in L_{f_X}^2\}$. The Poincaré inequality (2.18) bounds the variance of any function $h$ by its expected squared derivative. The constant $C_P(f_X)$ is informative about the underlying distribution: it is a measure of multimodality and can be related to the convergence rate of Markov chains (Pillaud-Vivien et al., 2020).

Roustant et al. (2017) introduced a computation procedure for Poincaré constants of univariate measures $f_X dx$ fulfilling the following assumption:

**Assumption 1 (Assumption on $f_X$).** Assume that $f_X$ is supported on a bounded interval $(a, b)$ and that $f_X(x) = e^{-V(x)}$ with $V$ continuous and piecewise $C^1$ on $[a, b]$.

Under this assumption, the Poincaré constant can be obtained by solving the following eigenvalue problem:

$$ \begin{gathered} L\psi_k = \psi_k'' - V'\psi_k' = -\lambda_k\psi_k, \quad (2.19) \\ \psi_k'(a) = \psi_k'(b) = 0, \end{gathered} $$

which can be expressed in weak form as

$$ \text{for all } h \in H_{f_X}^1 \qquad \langle h', \psi_k' \rangle = \lambda_k \langle h, \psi_k \rangle. \quad (2.20) $$

The solutions $(\lambda_k, \psi_k)$ have the following properties:

* After normalization, the eigenfunctions $(\psi_k)_{k \ge 0}$ form an orthonormal basis of $L_{f_X}^2$.

* The eigenvalues satisfy $0 = \lambda_0 < \lambda_1 < \dots \to \infty$.

* For the first eigenpair, $\lambda_0 = 0$ and $\psi_0(x) = 1$.

* For the second eigenpair, $C_P(f_X) = \frac{1}{\lambda_1}$, and $\psi_1$ attains equality in Eq. (2.18).

* If $f_X \in C^m$, then $\psi_k \in C^{m+1}$.

* $\psi_k$ has $k$ roots.

We show examples for univariate Poincaré bases in Figure 2.3. The Poincaré basis is in general not polynomial, except for the special case of the Gaussian distribution, for which the associated Poincaré basis coincides with the Hermite polynomials. The Poincaré basis associated to the uniform distribution can be shown to consist of cosine functions.

The eigenvalue problem in Eq. (2.19) can be solved numerically e.g. by one-dimensional finite elements. Since by Assumption 1 the measure is supported on a bounded interval, for most purposes it is enough to use an evenly spaced grid with sufficiently fine resolution.

As a consequence of Eq. (2.20), the Poincaré basis has several interesting properties related to derivatives. This is explored in our paper Lüthen et al. (2022c) available in Chapter 5.
<figure><img src="image_6.png" /><figcaption>Figure 2.3: The first five basis functions of Poincaré bases associated to different densities: Truncated Gaussian, uniform, triangular, and a custom density. The basis functions are plotted in colors, while the densities are shown in gray. We observe that higher-order functions oscillate more. The first eigenfunction is constant, and the functions have zero derivative at the boundary since they observe the homogeneous Neumann boundary condition (Eq. (2.19)).</figcaption></figure>

### 2.4.1.3 Truncation

In order to compute the chaos expansion in Eq. (2.11) for a given model, the coefficients of the expansion need to be determined. In general, infinitely many coefficients are needed for an exact representation, which is of course impossible in practice. Therefore, the expansion is truncated to a finite number of terms characterized by a multi-index set $A \subset \mathbb{N}^d$. The most widely used truncation for PCE is the total-degree truncation of order $p \in \mathbb{N}$, where the total degree of each retained basis element does not exceed $p$:

$$A^p = \{\alpha \in \mathbb{N}^d : \sum_{i=1}^{d} \alpha_i \le p\}. \qquad (2.21)$$

For the non-polynomial Poincaré chaos expansion, the word “degree” refers to the index of the univariate basis function or, equivalently, to the number of its roots.
More details about other possible truncation strategies and adaptive choices of truncation sets can be found in our publication Lüthen et al. (2022a), which is provided in Chapter 4.

### 2.4.1.4 Computation of the coefficients

Since the most popular chaos expansion in the field of UQ is the PCE, most of the UQ literature about coefficient computation discussed in this section is dealing with orthonormal polynomial bases, but many results can directly be applied or be generalized easily to other types of chaos expansions.

Historically, PCE were first applied in an *intrusive* context with the goal of solving stochastic differential equations, i.e., differential equations with random parameters (Ghanem and Spanos, 1991). In the *stochastic Galerkin method*, polynomial functions of one or several random variables serve as a basis for the random part of the solution, in the same way as the (finite element) functions on a spatial grid might serve as a basis for the spatial part of the solution. The stochastic Galerkin method is called *intrusive* because in order to perform an uncertainty analysis for a given deterministic differential equation, the solution scheme itself needs to be modified to incorporate the additional basis functions for the random part.

However, because the modification of existing code is inconvenient (and impossible for certain industrial applications, as the source code of the used computer models is often not available), most research efforts nowadays focus on *non-intrusive* PCE, which are computed purely from input-output data pairs of the computational model. In the following, we briefly discuss the three main non-intrusive approaches: projection, interpolation, and regression.

*Projection* relies on Eq. (2.12) to compute each coefficient of the expansion by projecting the model onto the respective basis function. The integrals can be evaluated using numerical integration, e.g. Gauss quadrature, with methods such as *non-intrusive spectral projection* (Le Maître et al., 2002) or *sparse pseudospectral projection* (Constantine et al., 2012; Conrad and Marzouk, 2013). While the number of model evaluations for naive tensor product quadrature scales as $(p+1)^d$, sparse quadrature based on the Smolyak scheme (Gerstner and Griebel, 2003; Conrad and Marzouk, 2013) can alleviate this to some degree (Eldred et al., 2008).

Whereas projection uses the pointwise model evaluations indirectly for computing approximations to the integrals in Eq. (2.12), another class of methods directly *interpolates* the model in the given pointwise evaluations. In other words, for a set of *experimental design points* $\mathcal{X} = \{x^{(1)}, \dots, x^{(N)}\} \subset \mathbb{R}^d$ the goal is to find the coefficient vector $c$ such that

$$ \mathcal{M}(x^{(i)}) = \mathcal{M}_{\mathcal{X}}^{\text{PCE}}(x^{(i)}; c) \text{ for all } i = 1, \dots, N, \qquad (2.22) $$

where the notation $\mathcal{M}_{\mathcal{X}}^{\text{PCE}}(\cdot; c)$ denotes the PCE adapted to the random vector $\mathcal{X}$ and using the coefficient vector $c$. This is sometimes called *stochastic collocation* (Xiu and Hesthaven, 2005; Babuška et al., 2007).<fnref n="12" /> The choice of points $\mathcal{X}$ crucially influences the stability of the

<footnote n="12">¹²There is some confusion around the term “stochastic collocation” in PCE because it has been used by various authors to refer to different classes of methods. Narayan and Zhou (2015) use this term for both interpolation and regression-based methods. It is used to refer to least-squares regression by e.g. Hosder et al. (2007); Eldred et al. (2008); Narayan et al. (2017). Xiu (2009) also counts pseudospectral projection as stochastic collocation method, to contrast it with stochastic Galerkin methods.</footnote>
interpolation procedure (Babuška et al., 2007; Narayan and Jakeman, 2014; Narayan and Zhou, 2015; Loukrezis et al., 2019).

Alternatively, the coefficients can be computed by regression, which aims to find an approximate solution to Eq. (2.22). One possibility is ordinary least-squares regression (OLS) (Isukapalli, 1999; Berveiller et al., 2006). Here, the model as well as the basis functions are evaluated on the experimental design $\mathcal{X}$, yielding the response vector $\mathbf{y} = (\mathcal{M}(\mathbf{x}^{(1)}), \dots, \mathcal{M}(\mathbf{x}^{(N)}))^{T} \in \mathbb{R}^{N}$ and the regression matrix $\mathbf{\Psi} \in \mathbb{R}^{N \times P}$ with entries $\Psi_{ij} = \psi_{\alpha_{j}}(\mathbf{x}^{(i)})$ for an enumeration $\mathcal{A} = \{\alpha_{j}\}_{j=1}^{P}$. To obtain the coefficients, the OLS problem

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

is solved, which has the closed-form solution $\hat{\mathbf{c}} = (\mathbf{\Psi}^{T}\mathbf{\Psi})^{-1}\mathbf{\Psi}^{T}\mathbf{y}$. The inversion requires $N \ge P$, where $P = (d+d')$ for a total-degree basis of degree $p$ in $d$ dimensions (Eq. (2.21)). To avoid overfitting, a heuristic recommendation in the engineering community is to choose $N \ge 2P$ or $\ge 3P$ (Hosder et al., 2007; Fajraoui et al., 2017).<fnref n="13" />

The continuous version of the OLS problem, using the notation of Eq. (2.22), is

$$ \hat{\mathbf{c}} = \arg \min_{\mathbf{c}} \mathbb{E} \left[ \left( \mathcal{M}(\mathbf{X}) - \mathcal{M}_{\mathbf{X}}^{\text{PCE}}(\mathbf{X}; \mathbf{c}) \right)^2 \right]. \qquad (2.24) $$

The analytical solution to this problem is given by Eq. (2.12), which makes the continuous version of OLS equivalent to the exact projection-based computation of the coefficients.

An alternative regression-based method that is often more sample-efficient than OLS is given by sparse regression, which we will explain in detail in Section 2.5.2.

## 2.4.2 Karhunen-Loève expansion

The chaos expansion described in the previous section was defined for deterministic models $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ – or, equivalently, for the induced random variables $Y = \mathcal{M}(\mathbf{X})$. Now, we turn to random fields $\{Y_x = \mathcal{M}(\mathbf{x}, \cdot)\}_{\mathbf{x} \in \mathcal{D}}$ as defined in Section 2.3.2, which are collections of random variables indexed over an index set $\mathcal{D}$.

### 2.4.2.1 The standard case

Let $\mu(\mathbf{x})$ be the mean function and $C(\mathbf{x}, \mathbf{x}')$ be the covariance function of the random field (see Section 2.3.2). Assume that $\mathcal{D}$ is closed and bounded, $C$ is continuous, and that $\mathcal{M}(\mathbf{x}, \cdot)$ has finite variance for all $\mathbf{x} \in \mathcal{D}$.

Then the Karhunen-Loève expansion (KLE) of the random field $\mathcal{M}$ is given by the following representation (Kosambi, 1943; Karhunen, 1946; Loève, 1978; Ghanem and Spanos, 1991):

$$ \mathcal{M}(\mathbf{x}, \omega) = \mu(\mathbf{x}) + \sum_{k=1}^{\infty} \sqrt{\lambda_k} \xi_k(\omega) \phi_k(\mathbf{x}) \qquad (2.25) $$

<footnote n="13">¹³In fact, rigorous theoretical analyses of least squares methods for orthonormal polynomials have resulted in the following more demanding bounds: in one dimension, $\frac{N}{\log N} \ge c_1 P^2$ for OLS (Migliorati et al., 2014); in $d$ dimensions, $N \ge c_2 P \log P$ (Hampton and Doostan, 2015a) or $\frac{N}{\log N} \ge c_3 P$ for certain weighted least-squares methods (Cohen and Migliorati, 2017); where in all cases, $c_i$ is a constant that does not depend on $N$ or $P$.</footnote>
where {$\phi_k$}$_{k=1,2,...}$ is an orthonormal basis of $L^2(D)$, $\lambda_1 \ge \lambda_2 \ge ... \to 0$ is a decreasing sequence of real numbers, and {$\xi_k$}$_{k=1,2,...}$ is a countable family of zero mean, unit variance, uncorrelated random variables. The quantities ($\lambda_k, \phi_k$) arise as solutions of the integral eigenvalue problem

$$ \int_D C(\mathbf{x}, \mathbf{x}')\phi_k(\mathbf{x}')d\mathbf{x}' = \lambda_k\phi_k(\mathbf{x}), \qquad (2.26) $$

which constitutes a homogeneous Fredholm integral equation of the second kind (Atkinson, 1996). The KL-random variables (KL-RV) $\xi_k$ are the result of the projection of $\mathcal{M}$ onto the basis:

$$ \xi_k(\omega) = \frac{1}{\sqrt{\lambda_k}} \int_D (\mathcal{M}(\mathbf{x}, \omega) - \mu(\mathbf{x})) \phi_k(\mathbf{x}) d\mathbf{x}. \qquad (2.27) $$

Note that both {$\phi_k$}$_{k=1}^\infty$ and {$\xi_k$}$_{k=1}^\infty$ are orthonormal bases, one capturing the “spatial” variation of $\mathcal{M}$ over $D$ (in terms of $\mathbf{x}$), the other capturing the stochastic variation of $\mathcal{M}$ (in terms of $\omega$). In that sense, Eq. (2.25) can be called a diagonalization of the random field (Venturi, 2011).

The covariance function can be expressed as the uniformly converging series (Mercer's theorem) (Grigoriu, 2002)

$$ C(\mathbf{x}, \mathbf{x}') = \sum_{k=1}^{\infty} \lambda_k \phi_k(\mathbf{x}) \phi_k(\mathbf{x}'). \qquad (2.28) $$

It follows that the average variance of the random field over the domain $D$ is equal to $\sum_{k=1}^\infty \lambda_k$. The series in Eq. (2.28) only depends on the spatial basis functions: the KL random variables do not enter this expression.

### 2.4.2.2 Non-Gaussian random fields

For a Gaussian random field, the KL random variables $\xi_k$ are standard Gaussian and independent. However, KLE also exists for non-Gaussian fields. The non-Gaussianity is governed by the joint distribution $f_\xi$, which might in general be complex and difficult to model or infer, leading some authors to conclude that KLE is mainly useful for Gaussian random fields (Grigoriu, 2006; Betz et al., 2014).

Still, several approaches for modeling the (joint) non-Gaussian distribution of KL random variables have been proposed. For given marginals and covariance function, Phoon et al. (2002, 2005) propose an algorithm which iteratively updates the distributions of the KL-RV in order to match the requested marginals. Other authors suggest to model the KL-RV by PCE using an iterative algorithm for the coefficients (Dai et al., 2019) or the method of fractional moments (Zhang et al., 2022). When only samples are given, the method of choice is often kernel density estimation (Grigoriu, 2010; Poirion and Zentner, 2013, 2014).

### 2.4.2.3 Extended KLE

While in Eqs. (2.25) to (2.27), KLE is defined for a closed and bounded index set $D$ and uses unweighted integrals, it can also be generalized to the weighted Hilbert space $L^2_{f_X}(D)$ with probability density $f_X$. It is then called extended KLE (Lemma et al., 2006).
### 2.4.2.4 Numerical computation

In practice, the infinite expansion in Eq. (2.25) must be truncated, resulting in

$$ \mathcal{M}(\mathbf{x}, \omega) \approx \mu(\mathbf{x}) + \sum_{k=1}^{M} \sqrt{\lambda_k} \xi_k(\omega) \phi_k(\mathbf{x}). \qquad (2.29) $$

K is usually chosen so that the fraction of explained variance satisfies

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

for a small threshold parameter $\epsilon > 0$. The basis $\{\phi_k\}_{k=1,2,...}$ is optimal in the sense that for any $M$, Eq. (2.29) has the smallest truncation error among all possible expansions in $L_{fX}^2(D)$.

In most cases, the integral eigenvalue problem in Eq. (2.26) does not have an analytical solution. The numerical solution is challenging because the integral is potentially high-dimensional and the quantities of interest are not only vectors but functions. Numerous methods for solving Eq. (2.26) have been proposed, such as Galerkin methods (Ghanem and Spanos, 1991; Gutiérrez et al., 1992; Schwab and Todor, 2006), the Nyström method (Atkinson, 1996; Betz et al., 2014), expansion optimal linear estimation (EOLE) (Li and Der Kiureghian, 1993), orthogonal series expansion (OSE) (Zhang and Ellingwood, 1994), and degenerate kernel methods (Atkinson, 1996; Betz et al., 2014). We refer to the comprehensive review by Betz et al. (2014) for more details.

### 2.4.2.5 Functional PCA

KLE can be seen as the continuous equivalent of the well-known dimension reduction method principal component analysis (PCA) (Jolliffe, 2002), which for a given random vector (or for realizations of a random vector) determines a new orthogonal coordinate system that is optimal for a truncated expansion. KLE does the same for random functions, and is therefore also known as functional PCA (Besse and Ramsay, 1986; Jolliffe, 2002; Ramsay and Silverman, 2005).

## 2.5 Sparsity and compressed sensing

Many numerical methods are affected by the so-called *curse of dimensionality*. This term denotes the explosion of the computational cost that often occurs when solving problems with a large number of parameters, which is referred to as their *dimension*. For example, the size of a PCE basis of total degree $p$ as defined in Eq. (2.21), which in $d$ dimensions is given by $P = \binom{p+d}{p}$, grows fast with increasing $d$. Solving a system of linear equations is another example: the number of required operations scales with $N^3$, where $N$ is the number of degrees of freedom (e.g., the number of nodes in a finite element mesh with increasingly fine resolution).

Due to the *curse of dimensionality*, many high-dimensional problems may exceed the capacity even of today's powerful supercomputers. While the capabilities of computers continue to grow due to technological progress, the *curse of dimensionality* will always be an issue as long as the
algorithm in question scales worse than the improvement in technology. The only solution here is to alleviate the curse of dimensionality by devising better algorithms.

Such a tool, able to mitigate the curse of dimensionality for certain applications, is given by sparsity. Let $v \in \mathbb{R}^P$ be a vector, sometimes also called signal.<fnref n="14" /> For example, $v$ could be the vector of coefficients for a representation of an object in a suitable basis. The signal $v$ is called sparse if most of its entries are equal to zero. More precisely, it is called s-sparse if it does not have more than $s$ nonzero entries. A closely related concept is compressibility: a signal is called compressible if most of its entries are close to zero and only a few of them are large in magnitude. In other words, a compressible signal is well-approximated by a sparse signal.

### 2.5.1 The setting of compressed sensing

The concept of sparsity has had a big impact on the field of signal processing with the introduction and rigorous analysis of compressed sensing (Candès et al., 2006; Donoho, 2006; Candès and Wakin, 2008), which is sometimes also known as compressive sensing or compressive sampling. Traditionally, signal processing was bound by the Nyquist rate: to reconstruct a time-dependent signal with a frequency of at most $B$ Hz from measurements, at least $2B$ samples at constant rate are necessary. The groundbreaking result of compressed sensing is the following: under the assumption that the representation of the signal in some suitable basis is s-sparse, and using a specific type of measurement acquisition, the signal can be reconstructed with much less than $2B$ samples – notably without knowing where in the sparse representation the $s$ nonzero entries occur.

To make these notions more precise, let $c \in \mathbb{R}^P$ be the sparse signal of interest.<fnref n="15" /> It is measured by the sensing matrix $\Psi \in \mathbb{R}^{N \times P}$, which results in the vector of measurements $y = \Psi c \in \mathbb{R}^N$, where usually $N < P$. In other words, $y$ is a linear combinations of the columns of $\Psi$, weighted by the entries of $c$. Since $N < P$, the resulting linear system of equations $\Psi c' = y$ is underdetermined. Equivalently, in the compressed sensing literature this situation is also often referred to as overcomplete basis, i.e., there are many more elements in the basis (columns in $\Psi$) than necessary to explain the measurements $y$.

The goal is to reconstruct the sparse signal from the measurements, i.e., to find the sparsest solution $c$ of the equation $\Psi c' = y$. Mathematically, this can be written using the $l^0$-"norm"<fnref n="16" /> $\|c\|_0 = \sum_i 1_{c_i \neq 0}$ counting the number of nonzero entries in $c$ as follows:

$$c = \arg \min_{c' \in \mathbb{R}^P} \|c'\|_0 \quad \text{subject to} \quad \Psi c' = y. \qquad (2.31)$$

The only way to solve this so-called $l^0$-minimization problem exactly is by a combinatorial search through all possible nonzero patterns for $c$. However, this is infeasible for large problem sizes.

<footnote n="14">¹⁴This term is due to the origins of this theory in the field of signal processing.</footnote>

<footnote n="15">¹⁵We assume for simplicity that the signal of interest $c$ is sparse in $\mathbb{R}^P$. More generally, we can assume that $c$ is sparse when expressed in a certain basis $\Phi \in \mathbb{R}^{P \times P}$, i.e., $c = \Phi a$ with $a$ sparse. Then, the compressed sensing methodology will be applied to find $a = \arg \min_{a' \in \mathbb{R}^P} \|a'\|_0$ subject to $\Psi \Phi a' = y$.</footnote>

<footnote n="16">¹⁶The $l^0$-"norm" is not actually a norm because the triangle inequality is not fulfilled.</footnote>
The convex relaxation of this problem is given by $\ell^1$-minimization, where $\|c\|_0$ is replaced by the $\ell^1$-norm $\|c\|_1 = \sum_i |c_i|$:

$$c = \arg \min_{c' \in \mathbb{R}^P} \|c'\|_1 \quad \text{subject to} \quad \Psi c' = y. \qquad (2.32)$$

Under certain assumptions on $\Psi$ and on the sparsity of $c$, the solution of Eq. (2.32) coincides with the solution of Eq. (2.31). At the same time, Eq. (2.32) is much easier to solve because this problem is convex. An exemplary mathematical result of compressed sensing, making precise what assumptions are needed, is given in Section 2.5.2 below. Making compressed sensing even more powerful, the same methodology can be used to compute best $s$-sparse approximations $c_s$ to compressible signals $c$, achieving a trade-off between the accuracy of the solution $c_s$ and the number $N$ of required measurements.

Considering signals that are sparse in some basis has two main advantages:

1. It improves the interpretability of the solution: it is immediately clear which basis elements are needed to explain the signal.

2. Due to sparsity, reconstruction is possible even in the case $N < P$, when the linear system of equations is underdetermined, which alleviates the curse of dimensionality.

Compressed sensing, and more generally sparsity, has had significant impact not only in the field of signal processing but also in many other fields, from biomedicine to information security (Kougioumtzoglou et al., 2020). Specifically, the field of engineering mechanics has benefitted from compressed sensing techniques, as discussed in the extensive review by Kougioumtzoglou et al. (2020). We will discuss sparsity in the context of engineering models in more detail in Section 2.5.3.

In the field of UQ, sparsity has been utilized in particular for sparse regression-based PCE, starting in 2008 with a method proposed by Blatman and Sudret (2008). We give an overview of the literature on sparse regression-based PCE in Section 2.5.4.

We now cite one exemplary result from the compressed sensing literature to illustrate the type of mathematical guarantee obtained by compressed sensing.

## 2.5.2 An exemplary result of compressed sensing

The following results are from Candès and Plan (2011), adapted to our notation. Let $c \in \mathbb{R}^P$ be an $s$-sparse signal. Let $\Psi \in \mathbb{R}^{N \times P}$ be the so-called sensing matrix (to be multiplied with the signal) and assume that each of its rows were sampled i.i.d. from a vector-valued distribution $F$ with the following two properties:

1. **Isotropy**: For a row $\Psi \sim F$, it holds that $E[\Psi^T \Psi] = \mathbf{1}_P$ (with $\mathbf{1}_P$ the $P \times P$ identity matrix).

2. **Incoherence**: For the coherence parameter $\mu(F)$ defined by

$$\mu(F) = \sup_{\Psi \sim F} \max_{k=1, \dots, P} |\Psi_k|^2, \qquad (2.33)$$

it holds that $\mu(F) < \infty$, where $\Psi_k$ is the $k$-th entry of a random row vector $\Psi \sim F$.<fnref n="17" />

<footnote n="17">¹⁷If such a bound does not exist (e.g., for the Gaussian example $F = \mathcal{N}(0, \mathbf{1}_P)$), a probabilistic bound can be defined so that $\mu(F)$ is with high probability an upper bound.</footnote>
For example, the isotropic distribution $F = N(0, 1_P)$ results in a sensing matrix $\Psi$ whose entries are i.i.d. standard Gaussian random variables. Its coherence can be shown to be $\mu(F) = 6 \log P$ (Candès and Plan, 2011).

PCE provides another example of an isotropic sensing matrix: assuming that the experimental design $\mathcal{X}$ has been sampled i.i.d. from the input distribution $f_X$, a PCE regression matrix $\Psi$ containing the basis function evaluations at these locations is isotropic, because the basis functions are orthonormal w.r.t. $f_X$: $\mathbb{E}_X[\psi_k(\mathbf{X})\psi_l(\mathbf{X})] = \delta_{kl}$. We will talk more about the application of compressed sensing theory to PCE in Section 2.5.4 below.

Now, consider the $\ell^1$-minimization problem in Eq. (2.32). Candès and Plan (2011) obtain the following result:

**Theorem 1 (Noiseless incoherent sampling (Candès and Plan, 2011)).** Let $c$ be a fixed but otherwise arbitrary $s$-sparse vector in $\mathbb{R}^P$ and pick any scalar $\beta > 0$. Then with probability at least $1 - \frac{5}{P} - e^{-\beta}$, $c$ is the unique minimizer to Eq. (2.32) with $\mathbf{y} = \mathbf{\Psi}c$ provided that [the number $N$ of rows of $\mathbf{\Psi}$, each sampled from $F$, fulfills]

$$N \geq C_{\beta} \mu(F) s \log(P). \qquad (2.34)$$

More precisely, $C_{\beta}$ may be chosen as $C_0(1 + \beta)$ for some positive numerical constant $C_0$.

If we wanted to solve $\mathbf{\Psi}c = \mathbf{y}$ by inversion of $\mathbf{\Psi}$, we would need at least $N = P$ points. In contrast, the bound in Eq. (2.34) grows only logarithmically in $P$, up to the factor $\mu(F)$ that depends on the way the sensing matrix $\mathbf{\Psi}$ is constructed. This makes the coherence parameter $\mu(F)$ an important property that crucially influences the efficiency of the procedure. A number of contributions in the PCE literature computes and optimizes the coherence for the distribution $F$ induced by the PCE regression matrix (Doostan and Owhadi, 2011; Yan et al., 2012; Hampton and Doostan, 2015a,b; Shin and Xiu, 2016; Jakeman et al., 2017; Narayan et al., 2017); see also Section 2.5.4 below.

The cited result relies on the coherence parameter $\mu(F)$ to guarantee recovery. Bounds similar to Eq. (2.34) have been obtained for other properties of the sensing matrix $\mathbf{\Psi}$, for example the spark (Donoho and Elad, 2003), the mutual coherence (Donoho et al., 2006), and the restricted isometry property (RIP) (Candès et al., 2008).

While the above result is for noiseless recovery, i.e., $\mathbf{\Psi}c = \mathbf{y}$, similar results exist for the case that $c$ is not exactly $s$-sparse but only compressible, or that $\mathbf{y}$ is corrupted by noise. For this case, Candès and Plan (2011) propose to solve the so-called LASSO problem<fnref n="18" />

$$\hat{c} = \arg \min_{c' \in \mathbb{R}^P} \|\mathbf{\Psi}c' - \mathbf{y}\|_2^2 + \lambda \|c'\|_1. \qquad (2.35)$$

This results in bounds on the norms $\|c - \hat{c}\|_1$ and $\|c - \hat{c}\|_2$ which involve the best $s$-sparse approximation to the true solution $c$, while requiring the same number of samples as the noiseless case in Eq. (2.34). $\lambda$ is a hyperparameter governing the trade-off between accuracy and sparsity, chosen as $\lambda = 10\sqrt{P}$ by Candès and Plan (2011). While Eq. (2.35) looks similar to ridge regression, which uses the $\ell^2$-norm for the constraint in the second summand and can be solved

<footnote n="18">¹⁸LASSO stands for “least absolute shrinkage and selection operator”.</footnote>
analytically, the LASSO problem does not have a closed-form solution (Hastie et al., 2001, Ch. 3.4.2).<fnref n="19" />

Besides Eqs. (2.31) and (2.35), there are many other problem formulations and algorithms that provably encourage sparse solutions. For example, a formulation equivalent to Eq. (2.35) is *basis pursuit denoising* (BPDN) given by

$$c = \arg \min_{c'} \|c'\|_1 \quad \text{s.t.} \quad \|\Psi c' - y\|_2 \le \sigma. \qquad (2.36)$$

Instead of measuring sparsity with the $ℓ^1$-norm $||c||_1 = \sum_i |c_i|$, other sparsity-enhancing formulations make use e.g. of a combination of different norms such as $ℓ^1 - ℓ^2$ minimization (Yin et al., 2015) or elastic net (Hastie et al., 2001), or of Bayesian formulations utilizing a sparsity-encouraging prior on the coefficients (Tipping, 2001; Tipping and Faul, 2003; Babacan et al., 2010). A plethora of algorithms has been developed for solving these formulations, from gradient-descent based methods such as spectral projected gradient-$ℓ^1$ (SPGL1) for (2.36) (van den Berg and Friedlander, 2008) to stepwise regression algorithms such as orthogonal matching pursuit (OMP) (Pati et al., 1993; Tropp and Gilbert, 2007) or least-angle regression (LARS) (Efron et al., 2004) to iterative algorithms for sparse Bayesian learning (Tipping and Faul, 2003; Wipf and Rao, 2004; Ji et al., 2008; Babacan et al., 2010), to name just a few. General surveys of compressed sensing algorithms are available, e.g., in Bruckstein et al. (2009); Qaisar et al. (2013); Zhang et al. (2015); Arjoune et al. (2017).

### 2.5.3 Sparsity in engineering models

In computational as well as experimental engineering, it is common to deal with sparse data, i.e., only a small number of point measurements of the system of interest in a few configurations is available (Kougioumtzoglou et al., 2020). This is the case, for example, in structural health monitoring where a small number of sensors on a structure is used to determine the state of the system (Bao et al., 2011). The sparse set of data must be used together with some structural assumptions to draw conclusions about the system of interest.

As elaborated in the extensive review by Kougioumtzoglou et al. (2020), compressed sensing has successfully been applied to several engineering subfields, from structural health monitoring to uncertainty modeling and uncertainty propagation. With our focus on surrogate modeling, this thesis belongs to the third subfield. In particular, we will use sparse PCE and its close relative, sparse PoinCE, to surrogate deterministic simulators, while for stochastic simulators we will employ a method that uses sparse PCE as one of its building blocks (see Chapters 3 to 6).

The successful application of compressed sensing requires a suitable basis in which the signal of interest is sparse. As evident from the success of sparse PCE for engineering models (see, e.g., Blatman and Sudret (2008); Chatterjee et al. (2019); Hariri-Ardebili and Sudret (2020)), the polynomial chaos basis seems to be well suited for this task. A tentative explanation for

<footnote n="19">¹⁹In the special case that $\Psi$ has orthogonal columns and is overdetermined, the solution to Eq. (2.35) can be obtained analytically by shrinking and truncating the coefficients of the OLS solution in Eq. (2.23) (soft thresholding) (Hastie et al., 2001, Ch. 3.4). However, this formula is not applicable in the underdetermined case, in which the OLS solution cannot be computed.</footnote>
this observation could be that many quantities of interest computed from engineering models
are effectively smooth functions of their input parameters, or arise from physical processes
that have (close to) linear, exponential or low-frequency oscillatory behavior and are thus well
approximated by polynomials.

Furthermore, due to independence, the orthonormal polynomial basis of PCE is constructed as a
tensor product of univariate bases (Section 2.4.1). By construction, most of the resulting terms
only involve a subset of all variables. The number of basis functions in a total-degree basis is
large, but often only a few of these functions are actually needed to represent the model. A few
heuristic arguments can be found in the literature to explain this observation:

* If the model is high-dimensional, it is often the case that only a few of its input parameters are important. Becker (2020) relates this to the well-known *Pareto principle* (also known as “80/20 rule”): a small number of variables is able to explain a large amount of variation. In other words, models often have *low effective dimensionality*.

* The response of a system is often dominated by main effects (i.e., terms involving only a single variable) and low-order interactions between variables, an observation called *sparsity-of-effects principle* (Montgomery, 2004).

However, it should be noted that not all engineering models are sparse in the PCE basis. Most
are only compressible, i.e., the magnitude of their PCE coefficients decays fast. At the same
time, there are also examples that are not well approximated by polynomials, e.g., functions
such as likelihoods with highly localized behavior, non-differentiable functions arising e.g. from
taking maxima, or discontinuous models. For such models, a possible solution are adaptive
domain partitioning schemes in which local surrogates are constructed (Le Maître et al., 2004;
Wan and Karniadakis, 2006; Moustapha and Sudret, 2019; Marelli et al., 2021b).

We will numerically explore the sparsity of engineering models in terms of the non-polynomial
Poincaré basis (Section 2.4.1.2) in Chapter 5.

2.5.4 Sparse regression for polynomial chaos expansions

Shortly after the breakthrough in compressed sensing around 2006 (Donoho, 2006; Candès et al., 2006), these ideas also arrived in the UQ community working on PCE. The first papers suggest stepwise regression algorithms to build a sparse solution step by step (Blatman and Sudret, 2008, 2010). Subsequently, the $l^1$-minimization formulation and the associated rigorous algorithms were applied to PCE, yielding *sparse PCE* (Blatman and Sudret, 2011; Doostan and Owhadi, 2011; Mathelin and Gallivan, 2012) which continues to receive considerable attention from both the engineering and the applied mathematics communities.

Regression-based PCE perfectly fits into the framework of compressed sensing. For an experimental design $\mathcal{X} = \{\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(N)}\} \subset \mathcal{D}$ sampled from the input distribution $f_{\mathbf{X}}$, the regression matrix $\Psi \in \mathbb{R}^{N \times P}$ is obtained by evaluating the orthonormal polynomial basis functions at the experimental design points, resulting in entries $\Psi_{ij} = \psi_{\alpha_j}(\mathbf{x}^{(i)})$. The columns of $\Psi$ are by design orthogonal in expectation: in the terminology of Section 2.5.2, $\Psi$ is isotropic. The sparse signal to be “recovered” is the vector of coefficients $c$ that defines the PCE through
$\mathcal{M}^{\text{PCE}}(\boldsymbol{x}) = \sum_{\alpha \in \mathcal{A}} c_{\alpha} \psi_{\alpha}(\boldsymbol{x})$. The measurements are the model evaluations collected into the response vector $\boldsymbol{y} = (\mathcal{M}(\boldsymbol{x}^{(1)}), \dots, \mathcal{M}(\boldsymbol{x}^{(N)}))^{\text{T}}$.

As introduced in Eq. (2.23), the ordinary least squares (OLS) regression problem is to find $\boldsymbol{c}$ with $\boldsymbol{c} = \arg \min_{\boldsymbol{c}' \in \mathbb{R}^P} \|\boldsymbol{y} - \boldsymbol{\Psi}\boldsymbol{c}'\|_2^2$. This problem has the unique solution $\boldsymbol{c}_{\text{OLS}} = (\boldsymbol{\Psi}^T\boldsymbol{\Psi})^{-1} \boldsymbol{\Psi}^T\boldsymbol{y}$ if the regression matrix $\boldsymbol{\Psi}$ has rank $P$, which requires $N \ge P$ (Berveiller et al., 2006).

In the case of PCE, a sparse solution is desired not because sparse solutions are better approximations but for computational reasons: a sparse solution of sufficient accuracy can often be obtained with considerably fewer model evaluations than a dense solution would require. However, if enough data was available, a denser solution might provide an overall better approximation.

To find sparse coefficient vectors $\boldsymbol{c}$ for given data, the PCE community has developed numerous approaches ranging from methods that closely follow recent developments in compressed sensing (Blatman and Sudret, 2011; Hampton and Doostan, 2015b; Jakeman et al., 2017; Alemazkoor and Meidani, 2018a) to heuristic algorithms which achieve sparsity by iteratively adding or removing nonzero coefficients from the vector (Blatman and Sudret, 2008, 2010; Shao et al., 2017; Alemazkoor and Meidani, 2017). Not all methods that are currently in use can be traced back to a specific article introducing it to the PCE community. There are a number of methods that have been used for years without an explicit introduction, e.g., orthogonal matching pursuit (OMP) (Pati et al., 1993; Tropp and Gilbert, 2007), spectral projected gradient-$\ell^1$ (SPGL1) (van den Berg and Friedlander, 2008), and subspace pursuit (Dai and Milenkovic, 2009).

Furthermore, as explained in Section 2.5.2, the properties of the regression matrix are important for the quality of the result. The regression matrix is induced by the sampling of the experimental design, whose choice can therefore influence the final result. Again, the literature on sampling schemes is vast. Currently used and investigated schemes include, among others, Latin hypercube sampling (LHS) (McKay et al., 1979; Shields and Zhang, 2016), D-optimal sampling (Kiefer and Wolfowitz, 1959; Diaz et al., 2018), and coherence-optimal sampling (Hampton and Doostan, 2015b). Some of these methods introduce weights into the regression problem.

Further contributions to the sparse PCE literature include preconditioning schemes (Alemazkoor and Meidani, 2018b), basis adaptivity (Jakeman et al., 2015; Hampton and Doostan, 2018), rotation of the input coordinates (Tipireddy and Ghanem, 2014; Yang et al., 2018; Tsilifis et al., 2019; Papaioannou et al., 2019), and sequential sampling of the experimental design (Fajraoui et al., 2017; Hampton and Doostan, 2018; Diaz et al., 2018).

For a thorough and extensive literature review, we refer to our publications Lüthen et al. (2021) and Lüthen et al. (2022a) provided in Chapters 3 and 4.

## 2.6 Sensitivity analysis

An important part of understanding the model behavior under uncertainty is the analysis of the relative importance of each of the input variables, which is known as sensitivity analysis (Saltelli et al., 2008). There are numerous approaches for sensitivity analysis from local to global methods, ranging from correlation-based methods suited mainly for linear models to the more
general but more costly variance-based or density-based approaches. For extensive reviews and comparisons of the available methods, see Iooss and Lemaître (2015); Borgonovo and Plischke (2016); Borgonovo (2017); Becker (2020); Da Veiga et al. (2021).

To assess the global sensitivity of a real-valued nonlinear model with independent input variables, one of the most widespread approaches is the variance-based Sobol' sensitivity analysis, which we describe in detail in Section 2.6.1. Sobol' indices determine the contribution of each input parameter (and of each possible combination of parameters) to the variance of the output.

In the case that partial derivatives of the computational model are available, the derivative-based global sensitivity measures (DGSM) can be computed. These indices are defined as the $L_{f_X}^2(D)$-norm of the partial derivatives over the input space. We present this method in Section 2.6.2.

Further notable developments in global sensitivity analysis include the following. Shapley indices were originally introduced in game theory to attribute the “created value” of a team to each of the team members. In UQ, they can be applied for sensitivity analysis using the “team” of input variables and a suitable definition for the created value, e.g., the unnormalized closed Sobol' index (Owen, 2014; Owen and Prieur, 2017). The density-based Borgonovo indices measure the expected $L^1$-difference between the unconditional output density and the density conditioned on the i-th variable (Borgonovo, 2007). These were generalized by Da Veiga (2015) by replacing the $L^1$-difference by other dissimilarity measures between distributions. Recently, Da Veiga (2021) proposed two types of moment-independent indices based on kernel mean embeddings and the maximum mean discrepancy measure, or the Hilbert-Schmidt independence criterion (HSIC), respectively. Both indices allow an ANOVA-like decomposition analogous to Eq. (2.38) and can be seen as density-based generalizations of Sobol' indices.

In Sections 2.6.1 and 2.6.2, we present two global methods in more detail, namely Sobol' indices and derivative-based global sensitivity indices (DGSM), which will be needed in Chapter 5. In Section 2.6.3 we discuss which sensitivity analysis methods are applicable in the more general cases when the input variables are dependent or when the output is not scalar and real-valued. The latter includes the case of stochastic simulators.

## 2.6.1 Hoeffding-Sobol' decomposition and Sobol' indices

Assuming that the input random variables $X_i, i = 1, \dots, d$ are independent, any function $\mathcal{M} \in L_{f_X}^2(D)$ can be decomposed uniquely as a sum of terms of increasing complexity

$$ \mathcal{M}(\boldsymbol{x}) = m_0 + \sum_{1 \le i \le d} m_i(x_i) + \sum_{1 \le i < j \le d} m_{i,j}(x_i, x_j) + \sum_{1 \le i < j < k \le d} m_{i,j,k}(x_i, x_j, x_k) + \dots + m_{1,\dots,d}(x_1, \dots, x_d), \quad (2.37) $$

where the terms satisfy $\int m_I(\boldsymbol{x}_I)f_{X_k}(x_k)d\boldsymbol{x}_k = 0$ for all $I \subset \{1, \dots, d\}$ and $k \in I$. This implies that the terms are mutually orthogonal, i.e., $\int m_I(\boldsymbol{x}_I)m_J(\boldsymbol{x}_J)f_X(\boldsymbol{x})d\boldsymbol{x} = 0$ for $I \neq J$. This decomposition is known as Hoeffding-Sobol' decomposition or ANOVA decomposition (Hoeffding, 1948; Efron and Stein, 1981; Sobol', 1993).
ANOVA stands for “analysis of variance”, which refers to the following decomposition resulting from taking the variance on both sides of Eq. (2.37):

$$ \mathrm{Var}[\mathcal{M}(\mathbf{X})] = \sum_{1 \le i \le d} \mathrm{Var}[m_i(X_i)] + \sum_{1 \le i < j \le d} \mathrm{Var}[m_{i,j}(X_i, X_j)] + \dots + \mathrm{Var}[m_{1, \dots, d}(X_1, \dots, X_d)]. \quad (2.38) $$

Since all these terms are non-negative, it is clearly visible which amount of the output variance is explained by variable $X_i$ alone: this quantity is called first-order Sobol' index and is given by

$$ S_i^1 = \frac{\mathrm{Var}[m_i(X_i)]}{\mathrm{Var}[\mathcal{M}(\mathbf{X})]} \quad (2.39) $$

The total influence of a variable $X_i$ on the output variance is given by the total Sobol' index

$$ S_i^T = \frac{1}{\mathrm{Var}[\mathcal{M}(\mathbf{X})]} \sum_{J:i \in J} \mathrm{Var}[m_J(\mathbf{X}_J)], \quad (2.40) $$

gathering all the contributions by terms that involve variable $X_i$. Due to the normalization, $0 \le S_i^1 \le 1$ and $0 \le S_i^T \le 1$. It holds that $S_i^1 \le S_i^T$, and the difference between these two indices quantifies the interaction effect involving $X_i$. If $S_i^1$ is large, $X_i$ is important. On the other hand, $S_i^T$ close to zero implies that $X_i$ does not have a large effect on the variance of the output; in other words, setting $X_i$ to a constant value will not have a large impact on the output uncertainty.

Sobol' indices can be computed by Monte Carlo-based methods (Homma and Saltelli, 1996; Jansen, 1999; Sobol', 2001; Janon et al., 2014; Becker, 2020). An alternative is to compute them from a surrogate model (Iooss and Lemaître, 2015; Le Gratiet et al., 2017). For this, polynomial chaos expansions are especially suitable, because the tensor product basis functions directly relate to the terms of the Hoeffding-Sobol' decomposition as follows (Sudret, 2008). Let $\mathcal{M}(\mathbf{X}) = \sum_{\alpha \in \mathbb{N}^d} c_\alpha \psi_\alpha(\mathbf{X})$. Then

$$ m_I(\mathbf{X}_I) = \sum_{\substack{\alpha: \alpha_i > 0 \text{ for } i \in I, \\ \alpha_j = 0 \text{ for } j \notin I}} c_\alpha \psi_\alpha(\mathbf{X}) \quad (2.41) $$

and it follows that the Sobol' indices can be computed from the PCE coefficients as

$$ S_i^1 = \frac{1}{D} \sum_{\substack{\alpha: \alpha_i > 0, \\ \alpha_j = 0 \text{ for } j \neq i}} c_\alpha^2, \quad (2.42) $$

$$ S_i^T = \frac{1}{D} \sum_{\alpha: \alpha_i > 0} c_\alpha^2 \quad (2.43) $$

with the total variance $D = \mathrm{Var}[\mathcal{M}(\mathbf{X})] = \sum_{\alpha \neq 0} c_\alpha^2$.

More generally, the formulas for the first-order and total Sobol' indices Eqs. (2.42) and (2.43) hold for any tensorized spectral expansion whose associated univariate bases each contain the constant function as first element.

## 2.6.2 Derivative-based global sensitivity measure (DGSM)

When the partial derivatives of the computational model are available, a global sensitivity index can be obtained by integration. The so-called derivative-based sensitivity measure (DGSM)
index of $\mathcal{M}$ with respect to $X_i$ (Sobol and Gresham, 1995; Kucherenko et al., 2009; Sobol' and Kucherenko, 2009) is defined by

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

This can be seen as a continuous version of the Morris sensitivity indices (Morris, 1991; Lamboni et al., 2013). DGSM can efficiently be computed by Monte Carlo simulation (Kucherenko et al., 2009). An alternative is the computation through PCE, which is possible analytically for certain polynomial families (Sudret and Mai, 2015).

### 2.6.3 Dependent input and non-scalar output

Most of the methods described above assume that the input variables are independent and that the output is real-valued and continuous. We now briefly discuss extensions to the more general cases of dependent input and non-scalar output. In particular, the latter includes sensitivity analysis for stochastic simulators.

When the input variables are dependent, several methods, including the popular Sobol' indices, are not applicable or give answers that are difficult to interpret. Efficiently conducting meaningful sensitivity analysis for dependent variables is an active area of research, and several solutions have been proposed from local polynomial regression (Da Veiga et al., 2009) over modified variance-based indices utilizing random variable transformations (Mara and Tarantola, 2012; Mara et al., 2015) to Shapley indices (Owen, 2014; Owen and Prieur, 2017; Da Veiga, 2021). For a comprehensive overview of the current state of the art in this field, we refer to Da Veiga et al. (2021).

For the treatment of time- or space-dependent vector-valued or functional output, Nagel et al. (2020) and Perrin et al. (2021) reduce the dimensionality of the output by principal component analysis (PCA) before applying Sobol' analysis to the principal components. Nagel et al. (2020) analyze the evolution of the resulting Sobol' indices with time, while Perrin et al. (2021) aggregate the indices using the generalized sensitivity index of Lamboni et al. (2011) and Gamboa et al. (2014). Density-based methods like the Borgonovo indices or the dissimilarity measure generalization by Da Veiga (2015) can take multivariate output into account by integrating the density measure over the whole multivariate output range. Finally, the kernel-based methods proposed by Da Veiga (2021) can be applied to output of any type (categorical, vector-valued, functional) provided that suitable kernels are used.

Sensitivity analysis can also be performed for stochastic simulators (Section 2.3). As discussed by Zhu and Sudret (2021c), Sobol' indices can be generalized to stochastic simulators in several ways. The first is to consider both the explicit and latent variables as input variables, resulting in classical Sobol' indices for the extended input space (Iooss and Ribatet, 2009). The second way is to consider a statistic of the output random variable, which effectively reduces the stochastic simulator to a deterministic one, and to perform Sobol' analysis for this statistic (Iooss and Ribatet, 2009; Azzi et al., 2020). Finally, Sobol' analysis can be applied to each of the deterministic trajectories of the stochastic simulator obtained by fixing the latent variables.
The resulting indices are random variables that depend on the latent variables (Hart et al., 2017; Navarro Jimenez et al., 2017).

Kernel-based methods provide another alternative for the sensitivity analysis of stochastic sim-
ulators. Using the random variable view of Eq. (2.5), a stochastic simulator can be seen as a
deterministic model with random variable output. Therefore, with a suitable choice of kernel
for probability distributions, the recently proposed kernel-based indices by Da Veiga (2021) can
be applied.

2.7 Software for uncertainty quantification

There is an increasing awareness in the computational science community that understanding
the influence of uncertainty on the model output is an important part of the quality assurance
of computational models, and that UQ methods can give valuable insights about the model
itself. UQ is an active field of research, and a vast amount of efficient and reliable techniques
have been and continue being developed to address the questions from Section 1.1. This also
poses a challenge: modelers and practitioners, while being experts in their respective field of
applied science or engineering, cannot always be aware of the most recent developments and
recommendations in UQ, and might perform outdated, inaccurate or inefficient analyses (Saltelli
et al., 2019).

In view of this, it is necessary to test the developed methods on practical problems, and to
benchmark them to provide recommendations. At the same time, it is important to make
the recent developments available to the public in form of software that is thoroughly tested
and easy to use, so that not every researcher has to implement the UQ method of interest
(which might be mathematically or computationally quite involved) from scratch. Fortunately,
numerous such software projects have been developed by experts in the field of UQ, written in
various programming languages and with different focus areas and target groups. Among others,
there are the Python libraries UQPy (Olivier et al., 2020), UQ[Py]Lab (Lataniotis et al., 2021),
and PyApprox (Jakeman, 2022); the C++/Python libraries OpenTURNS (Andrianov et al.,
2007) and UQTk (Sargsyan et al., 2022); the C++ toolkit Dakota (Adams et al., 2014); the R
sensitivity package (Iooss et al., 2021); the R/Shiny platform Lagun (Da Veiga et al., 2022);
and the Matlab-based software packages Cossan (Patelli et al., 2014) and UQLab (Marelli and
Sudret, 2014).

For the research conducted in the context of this thesis, we have made extensive use of UQLab.
UQLab is a general-purpose UQ software developed by the Chair of Risk, Safety and Uncertainty
Quantification at ETH Zürich. It is based on Matlab (Mat, 2020) and comprises numerous
modules covering all steps of the UQ framework (Section 2.1), notably the modeling of the
input, surrogate models like sparse PCE and Kriging, as well as sensitivity analysis, but also
further modules such as Bayesian inversion, random fields, or reliability analysis, just to name
a few. Selected results of our research related to sparse PCE solvers (Chapter 3) were already
made available in the recent release UQLab Version 1.4.0 (Marelli et al., 2021a).
Summary and outlook on the following chapters

This concludes our review of the state of the art in uncertainty quantification, in which we specifically focused on the topics relevant for this thesis, i.e., surrogate modeling by spectral methods, sparsity, and sensitivity analysis. In the subsequent Chapters 3 to 6, we report our research results which build on the presented material as follows.

In Chapters 3 and 4, we survey and classify the literature on sparse regression-based PCE, whose basics we reviewed in Section 2.4.1 (PCE) and Section 2.5 (sparse regression). Furthermore, we perform an extensive benchmark of the most promising methods using UQLab (Section 2.7).

In Chapter 5, we analytically and numerically investigate the properties of sparse Poincaré chaos expansions, relying on the theory presented in Section 2.4.1.2 (Poincaré chaos), Section 2.5 (sparse regression), and Section 2.6 (sensitivity analysis).

Finally, in Chapter 6 we propose a new sparse spectral surrogate model for stochastic simulators, which belongs to the class of surrogate models taking the random function view as explained in Section 2.3.3.2. We use sparse polynomial chaos expansion (Sections 2.4.1 and 2.5) and Karhunen-Loève expansion (Section 2.4.2) to construct a surrogate that is able to approximate non-Gaussian stochastic simulators.
CHAPTER 3

# Sparse polynomial chaos expansions: Literature survey and benchmark

This chapter contains the postprint of

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

as published in the *SIAM/ASA Journal of Uncertainty Quantification*, differing from the published paper only in terms of layout and formatting.

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

## Abstract

Sparse polynomial chaos expansions (PCE) are a popular surrogate modelling method that takes advantage of the properties of PCE, the sparsity-of-effects principle, and powerful sparse regression solvers to approximate computer models with many input parameters, relying on only few model evaluations. Within the last decade, a large number of algorithms for the computation of sparse PCE have been published in the applied math and engineering literature. We present an extensive review of the existing methods and develop a framework for classifying

<footnote n="1">¹First published in SIAM/ASA Journal of Uncertainty Quantification in Volume 9, Number 2, 2021, published by the Society for Industrial and Applied Mathematics (SIAM). Copyright © by SIAM and ASA. Unauthorized reproduction of this article is prohibited.</footnote>
the algorithms. Furthermore, we conduct a unique benchmark on a selection of methods to identify which approaches work best in practical applications. Comparing their accuracy on several benchmark models of varying dimensionality and complexity, we find that the choice of sparse regression solver and sampling scheme for the computation of a sparse PCE surrogate can make a significant difference, of up to several orders of magnitude in the resulting mean-squared error. Different methods seem to be superior in different regimes of model dimensionality and experimental design size.

## 3.1 Introduction

Computer modelling is used in nearly every field of science and engineering. Often, these computer codes model complex phenomena, have many input parameters, and are expensive to evaluate. In order to explore the behavior of the model under uncertainty (e.g., uncertainty propagation, parameter calibration from data or sensitivity analysis), many model runs are required. However, if the model is costly, only a few model evaluations can be afforded, which often do not suffice for thorough uncertainty quantification. In engineering and applied sciences, a popular work-around in this situation is to construct a surrogate model. A surrogate model is a cheap-to-evaluate proxy to the original model, which typically can be constructed from a relatively small number of model evaluations and approximates the input-output relation of the original model well. Since the surrogate model is cheap to evaluate, uncertainty quantification can be performed at a low cost by using the surrogate model instead of the original model. Therefore, surrogate modelling aims at constructing a metamodel that provides an accurate approximation to the original model while requiring as few model evaluations as possible for its construction.

In this article, we focus on *nonintrusive regression-based sparse polynomial chaos expansions* (PCE), which is a popular surrogate modelling technique, and within the last decade is has received attention from the communities of applied mathematics and engineering. PCE express the computational model in terms of a basis of polynomials orthonormal with respect to the input random variables (Xiu and Karniadakis, 2002) and work well for globally smooth problems, which are common in many engineering applications. In addition to being a surrogate model, PCE are also often used for uncertainty propagation and sensitivity analysis, since moments and Sobol' sensitivity indices can be computed analytically (Sudret, 2008). *Nonintrusive PCE* treat the model as a black box (unlike intrusive PCE commonly used for solving stochastic PDEs). It is often advantageous to compute a *sparse PCE*, which is an expansion for which most coefficients are zero. This can be justified by the *sparsity-of-effects principle* and by *compressibility*: The sparsity-of-effects principle is a heuristic stating that most models describing physical phenomena are dominated by main effects and interactions of low order (Montgomery, 2004). Furthermore, PCE of real-world models are usually either sparse or at least compressible, meaning that the PCE coefficients, sorted by magnitude, decay quickly. Additional advantages of sparse expansions are given in Section 3.2.2.

Within the last decade, a large number of articles has been published on the topic of regression-based sparse PCE, each containing promising improvements on how to perform sparse PCE but
often lacking a thorough comparison to previously published methods. In this work, we survey the state-of-the-art literature, develop a general framework into which the various approaches can be fit, and carry out a numerical benchmark of a selection of methods to assess which of the many sparse PCE methods perform best on a representative set of realistic benchmark models.

The paper is structured as follows. Section 3.2 contains the description of our framework for classifying the sparse PCE literature as well as the extensive literature review. Section 3.3 contains the benchmark description and the numerical results. Finally, conclusions are drawn in section 3.4. More detailed descriptions of selected sparse solvers and experimental design techniques are given in the Appendices using unified notation.

## 3.2 Framework and literature survey for sparse polynomial chaos expansions

### 3.2.1 Regression-based polynomial chaos expansions

Let $\mathbf{X}$ be a $d$-dimensional random vector on a domain $\mathcal{D} \subset \mathbb{R}^d$ with independent components and probability density function $f_X(\mathbf{x}) = \prod_{i=1}^d f_{X_i}(x_i)$. Let $L_{f_X}^2(\mathcal{D})$ be the space of all scalar-valued models with finite variance under $f_X$, i.e., $L_{f_X}^2(\mathcal{D}) = \{h : \mathcal{D} \to \mathbb{R} \mid \mathrm{Var}_X[h(\mathbf{X})] < +\infty\}$. Under certain assumptions on the input distribution $f_X$ (Xiu and Karniadakis, 2002; Ernst et al., 2012), there exists a polynomial orthonormal basis $\{\psi_\alpha : \alpha \in \mathbb{N}^d\}$ for $L_{f_X}^2(\mathcal{D})$. Since the components of $\mathbf{X}$ are assumed to be independent, the basis elements are products of univariate orthonormal polynomials and are characterized by the multi-index $\alpha \in \mathbb{N}^d$ of polynomial degrees in each dimension.

We consider a particular model $\mathcal{M} \in L_{f_X}^2(\mathcal{D})$ and denote by $Y = \mathcal{M}(\mathbf{X})$ the corresponding output random variable. $Y$ can be represented exactly through an infinite expansion in $\{\psi_\alpha : \alpha \in \mathbb{N}^d\}$. In practice, however, not all infinitely many coefficients can be computed, and we are interested in a truncated expansion

$$Y = \mathcal{M}(\mathbf{X}) \approx \mathcal{M}^{\mathrm{PCE}}(\mathbf{X}) = \sum_{\alpha \in \mathcal{A}} c_\alpha \psi_\alpha(\mathbf{X}) \quad (3.1)$$

whose accuracy depends on the choice of the finite set $\mathcal{A} \subset \mathbb{N}^d$ (i.e., on the basis elements used for the expansion) as well as on the coefficients $c_\alpha$. Several truncation techniques are described in Section 3.2.4.

To compute the coefficients, one well-known and practical approach is regression (Isukapalli, 1999; Berveiller et al., 2006)$^2$. The basic regression approach is ordinary least squares (OLS). Let $\{\mathbf{x}^{(k)}\}_{k=1}^N \subset \mathcal{D}$ be a sample of the input space called experimental design (ED). Let $\mathbf{y} = (y^{(1)}, \dots, y^{(N)})^T$ be the vector of model responses with $y^{(k)} = \mathcal{M}(\mathbf{x}^{(k)})$. Define the matrix

<footnote>$^2$Other, earlier approaches for computing the coefficients are stochastic Galerkin and stochastic collocation methods (Ghanem and Spanos, 1991; Xiu and Hesthaven, 2005; Shen et al., 2020). For a comparison of their performance to regression-based PCE, see, e.g., Berveiller (2005); Hosder et al. (2007); Doostan and Owhadi (2011); Mathelin and Gallivan (2012).</footnote>
of basis function evaluations $\Psi$ with entries $\Psi_{ij} = \psi_j(x^{(i)})$, where the basis functions are enumerated in an arbitrary way. Denoting the number of basis functions with $P$, we see that the regression matrix $\Psi$ is an $N \times P$-matrix. Then, the OLS regression problem can be written as

$$ \hat{c} = \arg \min_{c \in \mathbb{R}^P} \| \Psi c - y \|_2. \qquad (3.2) $$

For a unique and robust solution, a heuristic number of model evaluations is $N \approx 2P, 3P$ (Hosder et al., 2007; Fajraoui et al., 2017), which can be infeasible for high-dimensional or high-degree PCE approximations.

### 3.2.2 Sparse PCE

Sparse coefficient vectors are determined through sparse regression, which, in addition to a good regression fit, requires that the solution be sparse. This constraint on sparsity is realized e.g. by adding as a regularization term the $l^0$-"norm" or the $l^1$-norm of the coefficient vector to the OLS formulation of (3.2) (see Appendix 3.B for more details). Many sparse regression methods used in PCE were originally developed in the context of compressive sensing (Donoho, 2006; Candès et al., 2006). For an introduction to the concepts and ideas of compressive sensing, see, e.g., Candès and Wakin (2008); Bruckstein et al. (2009); Kougioumtzoglou et al. (2020).

Unlike OLS, compressive sensing methods allow one to use fewer design points $N$ than basis functions $P$ and still recover the true sparse solution, or find a good sparse approximation to it. This and its robustness to noise, both of which are induced by the sparsity constraint, are the main reasons why sparse PCE are preferred to full PCE in practical settings when the number of model evaluations necessary for OLS-based PCE would be infeasible to compute. Note that while the use of sparse PCE for engineering models can be justified by compressibility and the sparsity-of-effects heuristic (section 3.1), the main goal in sparse PCE is to compute a good surrogate model from a few model evaluations and not to find the sparsest possible expansion. The assumption of sparsity is used as a tool for finding robust solutions to underdetermined systems of linear equations.

The first publications on sparse regression-based PCE proposed greedy forward-backward selection algorithms (Blatman and Sudret, 2008, 2010) and introduced the LARS algorithm for sparse PCE (Blatman and Sudret, 2011). On the mathematical side, Doostan and Owhadi (2011) analyzes convergence properties for sparse Legendre PCE when the design points are sampled from the uniform distribution. Another early work is Mathelin and Gallivan (2012) demonstrating that sparse PCE are less costly and more accurate than PCE based on Smolyak sparse grids. Since then, a large number of articles has been published on the topic of sparse PCE suggesting new methods for specific aspects of the sparse regression procedure. In the following, we present a framework into which the existing literature can be fit. The framework provides an overview of the available choices and enables a structured comparison of their impact on the performance of the resulting sparse PCE. Naturally, some new combinations of methods arise that have not yet been considered in the literature.
<figure><img src="image_7.png" /></figure>

**Figure 3.1:** Framework for computing sparse PCE. For each component of the framework, a number of methods has been proposed in the literature. In the first part (section 3.2) of this paper, we review the literature for each of the components. Details on selected methods are given in Appendices 3.A and 3.B. In the second part (section 3.3), we conduct a benchmark of selected methods for the components marked in orange, performing a single iteration of the framework. Iterative basis adaptation and experimental design enrichment are not considered in this work and are left for future benchmarks.

### 3.2.3 Framework: Classifying the literature on sparse PCE

Here, we present the framework we developed in order to gain an overview of the extensive literature proposing new methods for computing sparse PCE. Figure 3.1 shows a sketch of this framework. To compute a sparse PCE, the first step is to choose a set $\mathcal{A}$ of candidate polynomials for the expansion (Section 3.2.4) as well as an experimental design (Section 3.2.5). The experimental design defines the locations of the model evaluations. Once the model evaluations are obtained, the sparse solution can be computed by applying a sparse regression solver (see Section 3.2.6). This solver often depends on a number of hyperparameters that have to be selected carefully in order to get good results. Then, a suitable model selection criterion is evaluated (Section 3.2.7). If the obtained solution is satisfactory, the process can be stopped. Otherwise, the basis can be adapted (usually augmented; see Section 3.2.4), and/or the experimental design can be enriched (see Section 3.2.5.4). This process is repeated until the value of the model selection criterion either is satisfactory or cannot be reduced.

In addition to the components shown in Figure 3.1, there are methods (we call them enhancements) that aim at generally improving the solution to the sparse regression problem by, e.g., adapting the input space or preconditioning the regression matrix. They are discussed in Section 3.2.8.

### 3.2.4 Choice of basis and basis adaptation

The approximation quality of a truncated PCE hinges on the polynomial functions available for building the surrogate model, which are characterized by the associated set of multi-indices $\mathcal{A} \subset \mathbb{N}^d$. We call the finite set of polynomials $\{\psi_\alpha : \alpha \in \mathcal{A}\}$ included in the current truncated PCE model basis and call its members candidate polynomials or candidate basis functions.
A sparse PCE algorithm will find nonzero coefficients only for a subset $\mathcal{A}^{\text{active}} \subset \mathcal{A}$ of the basis functions, which we call active basis functions. On the one hand, $\mathcal{A}$ should include enough candidate polynomials to facilitate a good approximation. On the other hand, unnecessary basis functions decrease the ratio $N/P$ of model evaluations to unknown coefficients and deteriorate the regression matrix. Therefore, it is beneficial to carefully select the polynomials to be included in the expansion.

The choice of basis is often motivated by the sparsity-of-effects principle, a heuristic guideline stating that most real-world models are well approximated by terms of low degree and low interaction order. The following are popular ways to construct a basis:

* Total-degree A total-degree basis of degree $p$ is defined by $\mathcal{A}^p = \{\alpha : \|\alpha\|_1 \le p\}$.

* Hyperbolic truncation Let $p$ be fixed. Define the q-norm-truncated basis

$$ \mathcal{A}^{p,q} = \{\alpha : \|\alpha\|_q \le p\} \qquad (3.3) $$

with $q \in (0,1]$ (Blatman and Sudret, 2011) and the quasi-norm $\|x\|_q = (\sum_{i=1}^d |x_i|^q)^{\frac{1}{q}}$. For $q = 1$, $\mathcal{A}^{p,1}$ is the total-degree basis of order $p$. For smaller $q$, this truncation scheme excludes terms with high interaction order while keeping univariate polynomials up to degree $p$.

* Interaction order The interaction order of the basis can be restricted by defining

$$ \mathcal{A}^{p,r} = \{\alpha \in \mathcal{A}^p : \#\{i : \alpha_i \neq 0, i = 1, \dots, d\} \le r\} \qquad (3.4) $$

(Blatman and Sudret, 2008; Marelli and Sudret, 2019). This is useful for reducing the number of basis functions especially in high dimensions and when it is known (e.g., for physical reasons) that only a certain number of variables might interact.

Instead of using a fixed basis $\mathcal{A}$, it can be beneficial to employ an iterative scheme which starts from a small set of basis functions (low-dimensional, low-order) and, after computing a sparse solution, repeatedly adapts the basis by including a set of the most promising candidate polynomials and possibly removing others. This is called basis adaptivity (not to be confused with Gaussian adaptation (Tipireddy and Ghanem, 2014); see also Section 3.2.8).

A simple instance of basis adaptivity is degree adaptivity (Blatman and Sudret, 2010), which is based on total-degree bases. The procedure starts with a basis of low total degree and iteratively increases the total degree of the basis. Finally, a model selection criterion is used to select the best basis and associated sparse solution. Similarly, q-norm and interaction order, or a combination of all three, can also be used to design a basis adaptation scheme (Blatman and Sudret, 2010, 2011). This basis adaptivity is solution-agnostic in the sense that it does not use any information from the solutions computed in previous runs for the augmentation of the basis. Another solution-agnostic method is the dimension- and order-incrementing algorithm of Alemazkoor and Meidani (2017). Two methods that adapt the basis based on the active terms of the previous sparse solution are forward neighbor basis adaptivity (Sargsyan et al., 2014; Jakeman et al., 2015) and anisotropic degree basis adaptivity (Hampton and Doostan, 2018). These approaches keep the size of the basis small by strictly controlling which functions are added to the basis, often starting with a constant surrogate model and adding dimensions
only when necessary. A discussion and benchmark of basis-adaptive methods for sparse PCE is
available in Lüthen et al. (2022a).

**3.2.5 Experimental design**

Generally, experimental design techniques aim to select points in order to achieve certain goals related to exploration of the space (space-filling design), or to achieve certain properties of the regression matrix such as (in expectation) orthonormal columns, small determinant, small condition number, etc. For regression-based PCE, there are several main classes of experimental design techniques:

* Sampling based on the input distribution. The samples are drawn from the input dis- tribution. Techniques like Latin hypercube sampling (LHS) can be used to improve the space-filling properties.

* Sampling from a different distribution (also called induced sampling (Guo et al., 2020)). A different distribution and associated basis are constructed that have better properties than the input distribution and its basis.

* *Choosing points according to an optimality criterion from a candidate set.* Certain prop- erties of the regression matrix are optimized by choosing the design points from a suitable candidate set.

Note that none of these sampling methods considers the evaluations **y** of the computational model. For model-aware sampling techniques (active and supervised learning), see Section 3.2.5.4.

The following sections contain an overview of available sampling methods for sparse PCE. Se-
lected experimental design techniques are described in more detail in Appendix 3.A. Note that
due to the large amount of literature on experimental design for sparse PCE, our review cannot
be exhaustive. There are many more approaches available, including the deterministic Weil
points (Zhou et al., 2014), sparse grids (Perkó et al., 2014), randomized or subsampled quadra-
ture points (Berveiller et al., 2006; Tang and Iaccarino, 2014; Guo et al., 2017), etc.

**3.2.5.1 Sampling based on the input distribution**

The most basic sampling method is Monte Carlo (MC) sampling, where the points are sampled independently from the input distribution (Doostan and Owhadi, 2011; Hampton and Doostan, 2015b). LHS (McKay et al., 1979) aims at distributing the design points in a more space-filling way than MC sampling, using a stratification of the input quantile space in each dimension. LHS is known to filter main effects; i.e., it reduces the variance of linear regression estimators when the quantity of interest is dominated by terms of interaction order one (Shields and Zhang, 2016). LHS with sample decorrelation can further reduce the variance (Owen, 1994). LHS can also be used together with a criterion such as maximin distance (maximize the minimal distance between the design points in quantile space) (Pronzato and Müller, 2012), where several LHS designs are generated and the one that optimizes the criterion is returned.
A generalization of LHS that combines it with stratified sampling is Latinized partially stratified sampling (Shields and Zhang, 2016), which filters both main effects and low-order interaction terms and has been shown to consistently outperform LHS in high-dimensional cases.

Other space-filling/low-discrepancy methods are Sobol' sequences (Sobol', 1967) and Halton sequences (Halton, 1960), which are deterministic but appear to be quasi-random and space-filling in low dimensions.

### 3.2.5.2 Sampling from a different distribution

Several methods consider the coherence parameter of a basis, defined by

$$ \mu(\mathcal{A}, \{\psi_{\alpha}\}) = \sup_{x \in \mathcal{D}} \max_{\alpha \in \mathcal{A}} |\psi_{\alpha}(x)|^2, \qquad (3.5) $$

which can be used to bound the number of samples needed for accurate recovery by $\ell^1$-minimization (Candès and Plan, 2011; Hampton and Doostan, 2015b). For Hermite and Legendre polynomial bases, the coherence parameter of a total-degree basis grows exponentially with the total degree $p$ (Rauhut and Ward, 2012; Yan et al., 2012; Hampton and Doostan, 2015b).

To construct a coherence-optimal design, a new probability distribution and its associated orthonormal basis are constructed that achieve minimal coherence (Hampton and Doostan, 2015b,a). The new basis can be derived from the original PCE basis by multiplying each member by a weight function. Coherence-optimal samples can be drawn by Markov Chain Monte Carlo (MCMC) (Hampton and Doostan, 2015b) or by rejection sampling (see section 3.A.2.3). A related sampling scheme is obtained by constructing a new probability distribution and associated orthonormal basis which have improved but not optimal coherence; however, the distribution is constructed to belong to some classical family and is therefore straightforward to sample. This is called asymptotic sampling (Hampton and Doostan, 2015b,a) and results in a Chebyshev distribution for uniform input, and in a uniform distribution (within a ball of degree-dependent radius) for Gaussian input. Numerical experiments confirm the expected performance gain of coherence-optimal sampling over both MC and asymptotic sampling, and of asymptotic sampling over MC in the case of low dimension $d$ and high total degree $p$. For high-dimensional problems with low degree, MC often performs better than asymptotic sampling (Hampton and Doostan, 2015b).

The so-called Christoffel sparse approximation (CSA) (Jakeman et al., 2017; Narayan et al., 2017; Cohen and Migliorati, 2017) is a related approach which constructs a new orthonormal basis that minimizes the quantity

$$ \tilde{\mu}(\mathcal{A}, \{\psi_{\alpha}\}) = \sup_{x \in \mathcal{D}} \left( \frac{1}{|\mathcal{A}|} \sum_{\alpha \in \mathcal{A}} |\psi_{\alpha}(x)|^2 \right)^{\frac{1}{2}}. \qquad (3.6) $$

As for coherence-optimal sampling, the new basis can be derived from the original basis by multiplying each member by a weight function, which results in a weighted regression problem. The corresponding probability distribution is chosen to be the so-called weighted pluripotential equilibrium measure, which for bounded distributions is the Chebyshev distribution. For
one-dimensional Gaussian input, this measure is a symmetric Beta distribution with degree-dependent bounds.

Note that all three sampling methods described in this section introduce weights and therefore modify the objective function into a weighted regression problem $W\Psi c \approx Wy$. Since the objective function belongs to the scope of the solver, these methods cannot be considered as pure sampling methods in the sense of being completely independent of the solver.

### 3.2.5.3 Choosing points according to an optimality criterion from a candidate set

The following methods choose points from a candidate set in order to optimize properties of the regression matrix. Candidate points can be sampled, e.g., using MC, LHS (Fajraoui et al., 2017), coherence-optimal sampling (Diaz et al., 2018; Alemazkoor and Meidani, 2018a), or Christoffel sampling (Shin and Xiu, 2016). Note that some of these methods introduce weights, resulting in a weighted regression problem. The candidate set can have a large influence on the resulting design.

* D-optimal sampling aims at maximizing the determinant $D(\Psi) = \det(\frac{1}{N}\Psi^T\Psi)^{\frac{1}{P}}$ of the information matrix (Kiefer and Wolfowitz, 1959). Note that $D(\Psi) = 0$ if $N < P$. Maximizing the D-value is connected to minimizing the variance of the coefficients of the PCE estimate (Zein et al., 2013). Algorithms for D-optimal designs include greedy augmentation (Dykstra, 1971), exchange techniques (Fedorov, 2013; Cook and Nachtsheim, 1980; Nguyen and Miller, 1992; Zein et al., 2013), maxvol (Mikhalev and Oseledets, 2018), gradient descent (Zankin et al., 2018), and rank-revealing QR decomposition (RRQR)/subset selection (Diaz et al., 2018; Gu and Eisenstat, 1996). The advantage of the last method is that it can also be applied for wide matrices $\Psi \in \mathbb{R}^{N \times P}$ where $N < P$.

* S-optimal sampling (also called “quasi-optimal” in Shin and Xiu (2016a)) selects samples from a large pool of candidate points so that the PCE coefficients computed using the selected set are as close as possible to the coefficients computed from the whole set of candidate points (Shin and Xiu, 2016a). The S-value is defined by<fnref n="3" />

$$S(\Psi) = \left( \frac{\sqrt{\det \Psi^T \Psi}}{\prod_{i=1}^{P} \|\Psi_i\|_2} \right)^{\frac{1}{P}} \quad (3.7)$$

where $\Psi_i$ denotes the $i$th column of the regression matrix. Its maximization has the effect of maximizing the column orthogonality of the regression matrix while at the same time maximizing the determinant of the information matrix (Shin and Xiu, 2016a). Note that $S(\Psi) = 0$ if $N < P$. An S-optimal experimental design can be computed using a greedy exchange algorithm (Shin and Xiu, 2016a,?; Fajraoui et al., 2017).

* Near-optimal sampling simultaneously minimizes the two matrix properties mutual coherence and average cross-correlation (Alemazkoor and Meidani, 2018a), both of which quantify the correlation between normalized columns of the regression matrix (see section 3.A.3.3 for the definitions of these properties). A near-optimal design can be built by

<footnote n="3">³This definition assumes that the columns of the matrix $\Psi_{cand}$, containing the evaluations of all candidate points, are mutually orthogonal.</footnote>
a greedy algorithm (Alemazkoor and Meidani, 2018a). Note that for near-optimal sampling, it is not necessary that $N \ge P$, since this method does not rely on the determinant of the information matrix.

### 3.2.5.4 Sequential enrichment of the experimental design

Instead of sampling the whole experimental design at once, it has been proposed to use sequential enrichment. Starting with a small experimental design, additional points are chosen based on the last computed sparse PCE solution or on an augmented basis. In the context of machine learning, sequential sampling is also known as active learning (Settles, 2012). Sequential enrichment has been proposed in the context of S-optimal sampling (Fajraoui et al., 2017), D-optimal sampling (Diaz et al., 2018), and coherence-optimal sampling (Hampton and Doostan, 2018). Zhou et al. (2019) suggest an enrichment strategy based on approximations to the expected quadratic loss function, i.e., the mean-squared error. Ji et al. (2008) and Seeger and Nickisch (2008) propose choosing points that minimize the differential entropy of the posterior distribution of the coefficients (using a Bayesian regression setting). In all cases, numerical examples show that the sequential strategy generally leads to solutions with a smaller validation error compared to nonsequential strategies. Due to the complexity of the topic and the already large extent of our benchmark, this strategy, albeit promising, is not explored further in this paper.

### 3.2.6 Solution of the minimization problem

There are many formulations of the regression problem that lead to a sparse solution, such as $\ell^0$-minimization, $\ell^1$-minimization (basis pursuit denoising (BPDN), LASSO), $\ell^1 - \ell^2$ minimization, Bayesian methods, etc. (see also Appendix 3.B). Based on these formulations, a vast number of sparse solvers has been proposed in the compressed sensing literature; see, e.g., Carron (2013) and the surveys of Qaisar et al. (2013); Zhang et al. (2015); Arjoune et al. (2017). We focus here on solvers that have been proposed in the context of sparse PCE. Of course, it is straightforward to use any other sparse solver to compute a sparse PCE.

The following solvers have been proposed in the sparse PCE literature:

* Convex optimization solvers. $\ell^1$-minimization in its various formulations is a (constrained) convex optimization problem. Least angle regression (LARS) (Efron et al., 2004; Blatman and Sudret, 2011; Marelli and Sudret, 2019) is an iterative method that adds regressors one by one according to their correlation with the current residual, and updates the coefficients following a least angle strategy. With the LARS-LASSO modification, which allows for backwards elimination of regressors, LARS is able to generate the whole LASSO path (Efron et al., 2004). Unmodified LARS can also be classified as a greedy method. SPGL1 (van den Berg and Friedlander, 2008; Van den Berg and Friedlander, 2015) solves the BPDN formulation by solving a succession of LASSO instances using the spectral projected gradient (SPG) method. Other solvers belonging to this class are e.g. the solvers implemented in $\ell_1$magic (Candès and Romberg, 2005) and SparseLab (Donoho et al., 2007).
* Greedy methods are variants of stepwise regression where the regressors are added to the model one by one according to some selection criterion, aiming at finding a heuristic solution to the intractable $ℓ^0$-minimization formulation. Orthogonal matching pursuit (OMP) (Tropp and Gilbert, 2007; Doostan and Owhadi, 2011; Marelli and Sudret, 2019) is a classical forward selection algorithm in which orthonormalized regressors are added to the model one by one according to their correlation with the residual, and the coefficients are computed by least-squares. Baptista et al. (2019) suggests extensions to OMP such as parallelization, randomization and a modified regressor selection procedure. Subspace pursuit (SP) (Dai and Milenkovic, 2009; Diaz et al., 2018; Diaz, 2018) is an iterative algorithm that repeatedly uses least squares on a subset of regressors. LARS (Efron et al., 2004; Blatman and Sudret, 2011) without the LASSO modification (allowing for removal of regressors) can also be classified as a greedy method. Another greedy method is ranking-based sparse PCE (Tarakanov and Elsheikh, 2019) which employs batch updating, coordinatewise gradient descent of the elastic net formulation, and a correlation- and stability-based ranking procedure for the regressors. Many more greedy stepwise regression techniques have been proposed, utilizing various selection criteria, solvers, and stopping criteria. An overview of methods following this scheme is given in section 3.B.6.

* Bayesian compressive sensing (BCS) (a.k.a. sparse Bayesian learning) is a class of methods that use a Bayesian setting to find a sparse solution. They impose a sparsity-inducing prior on the coefficients, whose parameters are again considered to be random variables with a hyperprior (Tipping, 2001; Ji et al., 2008; Sargsyan et al., 2014; Tsilifis et al., 2020). The solution is typically the maximum a posteriori estimate of the coefficients and can be computed e.g. by differentiation (Tipping, 2001), expectation-maximization (Figueiredo, 2003; Wipf and Rao, 2004), expectation-propagation (Seeger and Nickisch, 2008), variational inference (Tsilifis et al., 2020; Bhattacharyya, 2020), or a fast approximate algorithm (Faul and Tipping, 2002; Tipping and Faul, 2003). An extension called FastLaplace with an additional layer of hyperparameters has been proven to attain even sparser solutions (Babacan et al., 2010; Babacan, 2011). A greedy algorithm using the Bayesian setting to select the regressors is the greedy Bayesian Kashyap information criterion (KIC)-based algorithm (Shao et al., 2017).

* Iteratively reweighted methods. Iteratively reweighted $ℓ^1$-minimization uses the coefficients computed in a previous iteration to construct a weighted $ℓ^1$-minimization problem (Candès et al., 2008; Yang and Karniadakis, 2013). Cheng and Lu (2018b) suggest an iterative reweighted method with D-MORPH regression (Li and Rabitz, 2010) as its computational core, which is a technique that follows a certain path, defined by a quadratic objective function, on the manifold of solutions to the underdetermined system.

Each of the solvers mentioned above features one or more hyperparameters whose values must be calibrated. This is usually done by cross-validation. Popular choices are leave-one-out (LOO) cross-validation (accelerated for least-squares solutions) (Blatman and Sudret, 2010, 2011), LOO cross-validation with a modification factor for small sample sizes (Chapelle et al., 2002; Blatman and Sudret, 2011), and k-fold cross-validation (Doostan and Owhadi, 2011; Jakeman et al., 2015; Huan et al., 2018).
Selected sparse regression solvers are described in more detail in Appendix 3.B.

### 3.2.7 Model selection criterion

To decide whether to continue iterating in the framework or stop the process, we need to assess how well the current sparse solution performs. Our main quantity of interest is the generalization error, which quantifies the mean-square accuracy of the surrogate. It is given by

$$E_{\text{gen}} = \mathbb{E}_{\mathbf{X}} \left[ \left( \mathcal{M}(\mathbf{X}) - \mathcal{M}^{\text{PCE}}(\mathbf{X}) \right)^2 \right] \quad (3.8)$$

where $\mathcal{M}$ is the computational model, $\mathbf{X}$ is the random input vector, and $\mathcal{M}^{\text{PCE}}$ is the sparse PC surrogate.

The generalization error can be approximated by the validation error, which is the MC estimate of (3.8) on a validation set $\{(\mathbf{x}_{\text{val}}^{(i)}, y_{\text{val}}^{(i)}) : \mathbf{x}_{\text{val}}^{(i)} \sim \text{i.i.d.} f_{\mathbf{X}}, y_{\text{val}}^{(i)} = \mathcal{M}(\mathbf{x}_{\text{val}}^{(i)}), i = 1, \dots, N_{\text{val}}\}$. To make the validation error independent of the scaling of the model, it is convenient to use the relative mean-squared error defined by

$$\text{RelMSE} = \frac{\sum_{i=1}^{N_{\text{val}}} (y_{\text{val}}^{(i)} - \mathcal{M}^{\text{PCE}}(\mathbf{x}_{\text{val}}^{(i)}))^2}{\sum_{i=1}^{N_{\text{val}}} (y_{\text{val}}^{(i)} - \bar{y})^2} \quad (3.9)$$

where $\bar{y} = \frac{1}{N_{\text{val}}} \sum_{i=1}^{N_{\text{val}}} y_{\text{val}}^{(i)}$.

The best surrogate model $\mathcal{M}^{\text{PCE}}$ (defined by $\mathcal{A}$ and $\mathcal{c}$) is the one that has the smallest generalization error. In practical applications, the generalization error typically cannot be computed, and a large validation set is not available due to computational constraints. Instead, we define a model selection criterion that acts as a proxy for the generalization error. A typical stopping criterion in the PCE framework of Figure 3.1 is the observation that the model selection criterion no longer improves. The following model selection criteria have been proposed in the sparse PCE literature:

* k-fold cross-validation (CV) (Hastie et al., 2001; Jakeman et al., 2015; Hampton and Doostan, 2018), which approximates the validation error by building a surrogate several times on different subsets of the data, and evaluating the error on the remaining data points.

* Leave-one-out (LOO) cross-validation (Hastie et al., 2001; Blatman and Sudret, 2010, 2011), which is N-fold cross-validation (where N is the size of the experimental design). For PCE approximations computed by OLS, there exists an efficient formula to evaluate the LOO error (Blatman and Sudret, 2011, Appendix D).

* Modified LOO (Blatman and Sudret, 2011), which uses a correction factor for the LOO which was derived for the empirical error for OLS with small sample size (Chapelle et al., 2002). The correction factor depends on the experimental design and the active basis functions.

* Kashyap information criterion (KIC) (Shao et al., 2017; Zhou et al., 2019c), an approximation to the Bayesian model evidence, which is the likelihood of observations given the model.
* *Sparsity* (Alemazkoor and Meidani, 2017), which uses the idea that a larger basis should lead to a sparser solution when the necessary basis functions enter the model, unless the ratio of basis functions to model evaluations becomes too large.

A model selection criterion is also often used to determine the hyperparameter of the sparse solver (see Section 3.2.6).

If cross-validation is used to select the solver hyperparameter, this estimate of the validation error is often too optimistic due to model selection bias. Instead of reusing this estimate for model selection, it is better to perform an outer loop of *k*-fold or LOO cross-validation, a procedure called *double cross-validation* or *cross-model validation* (Baumann and Baumann, 2014; Liu et al., 2020b).

### 3.2.8 Further enhancements of sparse PCE

There are many enhancements to the simple scheme for sparse PCE presented in Figure 3.1. The following methods have been suggested to improve the accuracy of the solution and reduce the number of model evaluations needed:

* Alemazkoor and Meidani (2018b) construct a preconditioning matrix for a given regression matrix which reduces the mutual coherence while avoiding deterioration of the signal-to-noise ratio.

* Huan et al. (2018) suggest a technique called stop-sampling, which guides the decision of whether to obtain more samples (sequential ED enrichment) by observing the decrease of the CV error.

* In the case when prior information about the magnitude of the coefficients is available, Peng et al. (2014) use this information to construct a weighted regression problem which allows a more accurate solution with fewer points (similar to iteratively reweighted $\ell^1$-minimization).

* Liu et al. (2020b) use resampled PCE, which is a technique for improving the PCE solution by aggregating the results of several solver runs on different subsets of the data. Only the terms that are chosen most often by the solvers are retained in the final solution.

* Several methods exist to reduce the dimension of the input space before computing the sparse PCE. Unsupervised methods are principal component analysis (PCA) and kernel PCA (Lataniotis, 2019). “Basis adaptation” methods (referring to a basis of the input random space) determine a suitable rotation of the input space, often assumed to be independent standard Gaussian, into new coordinates which permit a sparser representation in fewer coordinates (see Tipireddy and Ghanem (2014); Yang et al. (2018); Tsilifis et al. (2019) and others). A related technique is nonlinear PCE-driven partial least squares (PLS) (Papaioannou et al., 2019; Zhou et al., 2020), which reduces the input dimension by identifying directions in the input space that are able to explain the output well in terms of a sum of one-dimensional PCEs.
## 3.3 Numerical results

### 3.3.1 Benchmark design

While the number of methods for computing sparse PCE is large, to the best of the authors' knowledge there is no comprehensive benchmark study on this topic. Most publications only compare the newly developed method to one or two baseline methods. An overview of publications containing comparisons of sparse PCE methods is presented in Appendix 3.D.

Since the number of possible combinations of sampling schemes, sparse regression solvers, basis adaptation schemes, model evaluation criteria etc. is huge (see Section 3.2.3 and thereafter), we restrict our benchmark as follows to some of the most promising and best-known methods:

* We consider the sampling schemes MC, LHS, coherence-optimal, and D-optimal. LHS is used together with a maximin criterion to improve the space-filling property (using the MATLAB function `lhsdesign`). D-optimal designs are constructed from a coherence-optimal candidate set<fnref n="4" /> using the subset selection/RRQR algorithm, which allows for the construction of D-optimal experimental designs with size *N* smaller than the number of regressors *P* (Diaz et al., 2018). There is no such algorithm for S-optimal sampling, which is why we do not consider the latter in this benchmark. We do not consider Sobol' sequences, since they have been shown to be outperformed by LHS in sparse PCE applications (Fajraoui et al., 2017). Near-optimal sampling can realistically be used only for rather small bases (*P* ∈ O(100)), since its algorithm scales as O(MP²), with *M* = 10*P* as suggested by (Diaz et al., 2018) (see below). We use it with a coherence-optimal candidate set for two models with small basis.

* We consider the sparse regression solvers LARS, OMP, subspace pursuit (SP), FastLaplace (which we call here BCS), and SPGL1. Each of these solvers involves at least one hyperparameter, whose range is chosen according to reasonable guesses. For LARS and OMP, the hyperparameter is the number *K* of selected regressors and its range is [1, min{*P*, *N* − 1}]. For SP, *K* must fulfill $2*K* \leq \min\{*P*, *N*\}$. For BCS and SPGL1, the hyperparameter *σ* is chosen from the range $\sigma^2 \in N \cdot \text{Var}[y] \cdot [10^{-16}, 10^{-1}]$ which resembles a suitable range of possible relative MSE values. The hyperparameter values of LARS and OMP are determined by modified LOO cross-validation, while the hyperparameters of SP, BCS, and SPGL1 are determined by *k*-fold cross-validation (Section 3.2.6). In addition, we consider a variant of SP which uses LOO cross-validation instead of *k*-fold cross-validation, which we name *SP*<sub>LOO</sub>.

* We only consider the nonadaptive setting, in which both the basis and the size of the experimental design are fixed before the sparse PCE is computed.

* For each model, we define a reasonable range of 5–7 experimental design sizes. Each experiment is repeated 30–50 times to account for statistical uncertainty. The experimental

<footnote n="4">⁴We have also conducted all benchmark experiments with D-optimal designs constructed from LHS candidate sets, but we do not display these results, because we found that in most cases, D-opt(LHS) sampling performs (significantly) worse than most other sampling schemes, and often worse than its candidate set LHS. This matches with the results of Fajraoui et al. (2017) who observed this in a sequential enrichment setting and with the LARS solver.</footnote>
designs are generated anew for each repetition and each ED size. All solvers are tested on
the same ED realizations.

* The coherence-optimal candidate sets from which the D-optimal designs are selected have size $M = 10P$ as in (Diaz et al., 2018). For computational reasons, they are not sampled completely anew for each replication, but are rather drawn uniformly at random without replacement from a larger set of size $2M = 20P$ as in (Diaz et al., 2018).

* Since we are interested in sparse PCE for the purpose of surrogate modelling, our main quantity of interest is the relative mean-squared error<sup>56</sup> (RelMSE) as defined in (3.9). We investigate the RelMSE for several models, sparse solvers, and experimental design techniques. Typically, the practical interest lies in small experimental designs.

* Since the experimental design is random, the resulting validation error is a random variable. We visualize the data with boxplots. When comparing the performance of different methods, we consider the median performance and the spread of the resulting validation error. However, often there can be considerable overlap of validation errors between methods.

3.3.2 Software

For the implementation of the benchmark, we use the general-purpose uncertainty quantification
software UQLab (Marelli and Sudret, 2014). UQLab supports the integration of other software
packages.<fnref n="7" /> We utilize the following code:

* UQLab for MC sampling and LHS (Marelli and Sudret, 2014).

* DOPT_PCE for D-optimal sampling (subset selection/RRQR) and subspace pursuit (Diaz et al., 2018; Diaz, 2018).

* An in-house developed rejection-based implementation of coherence-optimal sampling.

* An in-house implementation of near-optimal sampling based on the description by Ale-
  mazkoor and Meidani (2018a).

* UQLab for the solvers LARS and OMP (Marelli and Sudret, 2014).

* spgl1-1.9 for SPGL1 (van den Berg and Friedlander, 2008; Van den Berg and Friedlander, 2015).

* FastLaplace for the hierarchical implementation “FastLaplace” of BCS (Babacan et al., 2010; Babacan, 2011).

<footnote n="5"><sup>5</sup>Note that some authors such as Doostan and Owhadi (2011); Hampton and Doostan (2015b); Diaz et al. (2018); Alemazkoor and Meidani (2018a) choose to normalize instead by $\sum_{x \in \mathcal{X}_{\text{val}}} M(x)^2$ or use the unnormalized mean-squared error (Shin and Xiu, 2016). To assess the recovery of sparse vectors just as in compressed sensing, some consider the error in the coefficient vector instead of the error in the model approximation (Alemazkoor and Meidani, 2018a).</footnote>

<footnote n="6"><sup>6</sup>Since a typical application of PCE is the computation of moments and Sobol' indices, the error in these quantities is another possible performance measure. However, globally accurate prediction as considered in this paper is more challenging that the prediction of moments and Sobol' indices, which are accurate if the largest-in-magnitude coefficients are estimated accurately. If a globally accurate surrogate model can be constructed, typically also the moments and Sobol' indices are accurate.</footnote>

<footnote n="7"><sup>7</sup>A description of how to use custom sparse solvers and sampling schemes in the UQLab framework can be found in the supplementary material.</footnote>
### 3.3.3 Benchmark: Considered models

Our benchmark is performed on a selection of 11 computational models of varying complexity and input dimensionality, which are typical benchmark models in the context of sensitivity and reliability analysis. An overview of these models is given in Table 3.1. For details on the models, we refer the reader to the respective publications supplied in the last column of the table. While of course not representative of all possible classes of engineering models, we believe that this sample provides a good testing ground for the comparative performance among different approaches for computing sparse PCE.

<table><caption>**Table 3.1:** Overview of the 11 computational models used in our benchmark. Italic font denotes finite element (FE) models, all other models are analytical. For each model, a static total-degree basis with hyperbolic truncation defined by *p* and *q* is used. The values are chosen to fulfill *P* ≈ 10/3 *N*<sub>max</sub>, where *N*<sub>max</sub> is the largest tested experimental design size. The values for *p* in parentheses for the Ishigami and borehole models refer to the smaller basis used in Section 3.3.6. The column “Reference” provides the relevant literature in which the models and their probabilistic inputs are described in detail.</caption><thead><tr><th>Model</th><th>Dimension</th><th>Input distributions</th><th>Basis</th><th>N<sub>max</sub></th><th>Reference</th></tr></thead><tbody><tr><td>Ishigami function</td><td>3</td><td>uniform</td><td>p = 14 (12),<br>q = 1</td><td>200</td><td>(Blatman<br>and Sudret, 2011)</td></tr><tr><td>Undamped oscillator</td><td>6</td><td>Gaussian</td><td>p = 5,<br>q = 1</td><td>150</td><td>(Echard<br>et al., 2013)</td></tr><tr><td>Borehole function</td><td>8</td><td>Gaussian, lognormal,<br>uniform</td><td>p = 5 (4),<br>q = 1</td><td>300</td><td>(Harper and<br>Gupta, 1983)</td></tr><tr><td>Damped oscillator</td><td>8</td><td>lognormal</td><td>p = 5,<br>q = 1</td><td>400</td><td>(Dubourg,<br>2011)</td></tr><tr><td>Wingweight function</td><td>10</td><td>uniform</td><td>p = 4,<br>q = 1</td><td>300</td><td>(Forrester<br>et al., 2008)</td></tr><tr><td>Truss model</td><td>10</td><td>lognormal, Gumbel</td><td>p = 4,<br>q = 1</td><td>300</td><td>(Blatman<br>and Sudret, 2011)</td></tr><tr><td>Morris function</td><td>20</td><td>uniform</td><td>p = 8,<br>q = 0.5</td><td>400</td><td>(Blatman<br>and Sudret, 2010b)</td></tr><tr><td>Structural frame model</td><td>21</td><td>lognormal, Gaussian;<br>dependent input<br>variables</td><td>p = 8,<br>q = 0.5</td><td>400</td><td>(Blatman<br>and Sudret, 2010)</td></tr><tr><td>2-dim diffusion model</td><td>53</td><td>Gaussian</td><td>p = 4,<br>q = 0.5</td><td>500</td><td>(Konakli<br>and Sudret, 2016)</td></tr><tr><td>1-dim diffusion model</td><td>62</td><td>Gaussian</td><td>p = 4,<br>q = 0.5</td><td>500</td><td>(Fajraoui<br>et al., 2017)</td></tr><tr><td>100D function</td><td>100</td><td>uniform</td><td>p = 4,<br>q = 0.5</td><td>1400</td><td>UQLab<br>example<sup>8</sup></td></tr></tbody></table>

In addition to analyzing aggregated performance on all 11 models, we investigate the behavior of