<table><caption>**Table 3.2:** Borehole function: Input random variables and their distributions</caption><thead><tr><th>Variable</th><th>Distribution</th><th>Description</th></tr></thead><tbody><tr><td>r<sub>w</sub></td><td>N(0.10, 0.0161812)</td><td>borehole radius</td></tr><tr><td>L</td><td>U([1120, 1680])</td><td>borehole length</td></tr><tr><td>K<sub>w</sub></td><td>U([9855, 12045])</td><td>borehole hydraulic conductivity</td></tr><tr><td>T<sub>u</sub></td><td>U([63070, 115600])</td><td>transmissivity of upper aquifer</td></tr><tr><td>T<sub>l</sub></td><td>U([63.1, 116])</td><td>transmissivity of lower aquifer</td></tr><tr><td>H<sub>u</sub></td><td>U([990, 1110])</td><td>potentiometric head of upper aquifer</td></tr><tr><td>H<sub>l</sub></td><td>U([700, 820])</td><td>potentiometric head of lower aquifer</td></tr><tr><td>r</td><td>Lognormal([7.71, 1.0056])</td><td>radius of influence</td></tr></tbody></table>

the methods in detail on a subset of four *spotlight models*, each of which possesses characteristic properties that might influence the approximation quality of sparse PCE methods: the *Ishigami function* is low-dimensional and highly compressible in the PCE basis but requires a high-degree basis to be approximated accurately. The *borehole function* is smooth and nonlinear and therefore is an example for a well-behaved engineering model. A *two-dimensional diffusion model*, a stochastic heat diffusion PDE in two physical dimensions, is high-dimensional, not analytical, and the magnitude of its expansion coefficients decays only slowly. Finally, the *100D function* is high-dimensional, analytical, and compressible.

The Ishigami model is the well-known three-dimensional, highly nonlinear, smooth analytical function

$$f(X_1, X_2, X_3) = \sin(X_1) + a \sin^2(X_2) + b X_3^4 \sin(X_1) \quad (3.10)$$

taking uniform input $\mathbf{X} \sim U([-\pi, \pi]^3)$. A typical choice is $a = 7, b = 0.1$. For this function, any sparse solver should be able to find a sparse solution.

The borehole function simulates the water flow through a borehole between two aquifers (Harper and Gupta, 1983). It is an eight-dimensional nonlinear function which, despite having an analytical form, is not trivial to approximate. It is defined by

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

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

The two-dimensional heat diffusion model (Konakli and Sudret, 2016) is defined by the partial differential equation (PDE)

$$-\nabla \cdot (\kappa(\mathbf{x})\nabla T(\mathbf{x})) = Q1_A(\mathbf{x}) \quad \text{in } \Omega = [-0.5, 0.5]^2 \quad (3.12)$$

with boundary conditions $T = 0$ on the top boundary and $\nabla T \cdot \mathbf{n} = 0$ on the left, lower, and right boundaries of the square domain $\Omega$, where $\mathbf{n}$ denotes the outer unit normal (see Konakli and Sudret (2016) for an illustration of the setup). Here, the source is in $A = [0.2, 0.3]^2$ with strength $Q = 500$. The output quantity of interest is the average temperature in $B = [-0.3, -0.2]^2$. The diffusion coefficient $\kappa(\mathbf{x})$ is modelled by a lognormal random field with mean $\mu_\kappa = 1$ and standard deviation $\sigma_\kappa = 0.3$. The autocorrelation function of the underlying Gaussian random field is an isotropic squared-exponential with length scale $l = 0.2$. The random field $\kappa(\mathbf{x})$ is discretized using the EOLE method (Li and Der Kiureghian, 1993) with $d = 53$ terms, which
comprises 99% of its variance. The solution to an individual heat diffusion problem is computed using an in-house finite element code (Konakli and Sudret, 2016). The input comprises $d = 53$ independent standard normal random variables.

Finally, the so-called 100D function is an analytical model of the form

$$f(\mathbf{X}) = 3 - \frac{5}{d} \sum_{i=1}^{d} iX_i + \frac{1}{d} \sum_{i=1}^{d} iX_i^3 + \frac{1}{3d} \sum_{i=1}^{d} i \ln (X_i^2 + X_i^4) + X_1X_2^2 + X_2X_4 - X_3X_5 + X_{51} + X_{50}X_{54}^2 \quad (3.13)$$

taking uniform inputs $X_i \sim U([1, 2])$, $i \neq 20$, and $X_{20} \sim U([1, 3])$. We use $d = 100$. This function was designed for sensitivity analysis: the first-order sensitivity indices of the input variables are generally nonlinearly increasing with their index, with certain variables having especially high sensitivity. The model also contains four interaction terms. It is an example from UQLab.<fnref n="8" />

For each of the models, we use a fixed basis for the benchmark. In general, the best total degree $p$ and the hyperbolic truncation $q$ are a priori unknown. We heuristically choose $q = 1$ for low-dimensional models ($d \le 10$) and $q = 0.5$ for high-dimensional models ($d \ge 20$). The degree $p$ is chosen so that the number of basis functions $P$ is approximately $\frac{10}{3}N_{max}$, where $N_{max}$ is the largest number of experimental design points for the specific benchmark. This choice is based on the reasoning that for an experimental design of size $N$, sparse solvers like LARS often select an active basis of size $\approx \frac{N}{3}$, and that the candidate basis might be 10 times larger than the final active basis to be sufficiently rich. We focus on rather small experimental designs, since our goal is not to investigate the convergence of the methods as $N \to \infty$ (which has been demonstrated elsewhere), but to decide which methods are most efficient for small $N$. This results in the choice of values for $p$, $q$, and $N_{max}$ displayed in Table 3.1.

### 3.3.4 Results: Comparison of solvers

First, we use a fixed sampling scheme (LHS) to compare the performance of the six solvers LARS, OMP, subspace pursuit (SP) using $k$-fold cross-validation, Subspace Pursuit using LOO cross-validation (SPLOO), FastLaplace (BCS), and SPGL1 on all 11 benchmark models described above.

In Figure 3.2 we display boxplots (50 replications) of relative MSE against experimental design size for all six solvers for the four spotlight models. For the remaining seven models, the corresponding boxplots (30 replications) of relative MSE against experimental design size are provided in Figure 3.C.1 in Appendix 3.C. In the plots, the lines as well as the dot inside the white circle denote the median of the relative MSE. We make the following observations:

* For the smallest experimental designs, all solvers perform similarly poorly: there is not enough information in the ED to construct an accurate surrogate model. For larger experimental designs, there can be considerable differences between the solvers' generalization errors of up to several orders of magnitude, which demonstrates that the solvers do not use the available information in identical ways.

<footnote n="8">⁸https://www.uqlab.com/sensitivity-high-dimension</footnote>
* BCS, $SP_{LOO}$, and OMP are often among the best solvers. BCS performs especially well for smaller experimental design sizes. In the case of the Ishigami model, it seems to plateau earlier than the other solvers. It also tends to find sparser solutions than the other solvers (not shown in plots), which might explain these observations: sparse solutions are advantageous when only limited data is available, but at the same time they carry the risk of ignoring important terms. For large ED sizes and the two highly compressible models (Ishigami and 100D function), BCS has a larger spread than the other solvers, possibly because the sparsity-enforcing procedure does not always include all of the important terms. In contrast, the greedy solver OMP returns rather dense solutions, which seem to generalize well. $SP_{LOO}$ performs well in general for low-dimensional models.

* SP does not perform well for small ED sizes, but for large ED sizes it sometimes outperforms the other solvers. Together with BCS, it tends to find sparser solutions than the other solvers (not shown in plots).

* LARS and SPGL1 often achieve a similar generalization error, which is often larger, sometimes significantly, than that of the other solvers. SPGL1 tends to return rather dense solutions (not shown in plots), which might not generalize as well as other solutions.

* Some models characterized by relatively poor compressibility (e.g. diffusion and frame) show comparable performance among all solvers. This is expected, as in such cases the sparsity assumption is a rather weak proxy for solution quality.

* The exceptions to the general observations outlined above are the damped oscillator and the Morris function, for which the solver performance is reversed, with LARS and SPGL1 among the best solvers, and OMP and $SP_{LOO}$ among the worst. In these cases, however, none of the methods achieves satisfactory accuracy within the available computational budget.

To objectively assess the performance of the methods, we now aggregate the results across models. Since all solvers are tested on the same set of experimental designs, we can determine the ranking of solvers for each experimental design (50 replications × 6 – 7 ED sizes for the four spotlight models, and 30 replications × 5 – 7 ED sizes for each of the seven additional models, resulting in 2620 PCEs) and count how often each solver achieved each rank. This is displayed in Figure 3.3 in the form of stacked bar plots, where the counts are given as percentages. The counts have been normalized by the number of replications and ED sizes used for each model, so that each of the models contributes equally to the final percentages.

This ranking alone, however, does not provide a complete picture; e.g., a solver ranked last can be off by orders of magnitude or barely worse than the best-performing one. Therefore, we added an additional set of triangle markers detailing for how many of the EDs the respective solver returned a result that was within two, five, or 10 times the smallest relative MSE attained by any of the six solvers on the same ED. For example, the red triangle in the top row of Figure 3.3a indicates that in ca. 30% of the runs, the solution returned by LARS had a validation error that was at most twice as large as that of the best solution on the same ED.

We have grouped the analysis results separately for the 6 low-dimensional ($d \le 10$) and the 5 high-dimensional ($d \ge 20$) models, because we observed that dimension had a significant impact on the rankings, and this information is readily available even for black-box models. We also
<figure><img src="image_8.png" /><figcaption>**Figure 3.2:** Boxplots of relative MSE against ED size from the solver benchmark for the four spotlight models. Results for six sparse solvers and LHS design. For the corresponding plots for the seven remaining models, see Figure 3.C.1 in Appendix 3.C.</figcaption></figure>

analyze small and large ED sizes separately, where the first half of considered ED sizes (3–4) are regarded small, and the second half large.

We make the following observations:

* *Low-dimensional models, small ED sizes:* BCS is the best solver most often (31% of runs) and also most often within two times the smallest error (72% of runs). However, it comes within one order of magnitude of the smallest relative MSE (a property which we here call *robustness*) in only 87% of runs, while $SP_{LOO}$ achieves this in 94% of runs. In the two former metrics, $SP_{LOO}$ comes second. LARS and SP perform worst, followed by OMP and SPGL1. Here and in the following, we observe that OMP and $SP_{LOO}$ often result in quite robust solutions.

* *Low-dimensional models, large ED sizes:* $SP_{LOO}$ outperforms the other solvers in all respects. It provides the smallest relative MSE of all solvers in 42% of runs, is in 73% of runs within two times of the smallest relative MSE, and is even in 99% of runs within one order of magnitude of the smallest relative MSE. For the other solvers, the statistics confirm the observations highlighted by the spotlight models: LARS and SPGL1 overall do not perform well. SP and BCS come close to the best solution quite often, but are less robust, whereas OMP is robust but often not as close to the best solution.

* *High-dimensional models:* For the small as well as the large ED sizes, we see that BCS
<figure><img src="image_9.png" /><figcaption>(a) low-dim models ($d \le 10$), small ED sizes</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) low-dim models ($d \le 10$), large ED sizes</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) high-dim models ($d \ge 20$), small ED sizes</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) high-dim models ($d \ge 20$), large ED sizes</figcaption></figure>

**Figure 3.3:** Aggregated results for the solver benchmark (Section 3.3.4), separately for low-dimensional (a),(b) and high-dimensional models (c),(d). Left: small ED sizes. Right: large ED sizes. For each model and experimental design, the ranking of the six solvers is determined. The stacked bar plots visualize how often each solver achieved the respective rank. The triangle markers in hues of red additionally demonstrate in how many runs the obtained relative MSE was within a factor of {2, 5, 10} of the smallest relative MSE achieved on this experimental design.

It is the best solver in 45% (53%) of runs and comes within two times the smallest relative MSE in even 96% (97%) of runs. Both BCS and OMP attain in all cases a relative MSE within 10 times the smallest relative MSE. SP performs better for large rather than small ED sizes. While SP<sub>LOO</sub> performed best for the low-dimensional models, here it shows poor performance. Note that all solvers come within 10 times the smallest relative MSE in more than 89% of all runs, showing that the choice of solver has a smaller impact for high-dimensional than for low-dimensional models.

### 3.3.5 Results: Comparison of sampling schemes together with solvers

We pair the five solvers LARS, OMP, SP, SP<sub>LOO</sub>, and BCS<sup>9</sup> with the sampling schemes MC, LHS, coherence-optimal, and D-optimal based on a coherence-optimal candidate set. We use the abbreviations coh-opt and D-opt(coh-opt) for the latter two. Since SP<sub>LOO</sub> performed poorly

<footnote n="9">⁹Since SPGL1 did not perform well in the previous section, and is quite slow, we do not include it further in this benchmark.</footnote>
for high-dimensional models in the solver benchmark of the previous section, we do not consider
it here for the high-dimensional models.

We run the benchmark for the low-dimensional models Ishigami, undamped oscillator, borehole,
damped oscillator, and wingweight function, and for the high-dimensional models Morris func-
tion, two-dimensional diffusion, one-dimensional diffusion, and 100D function. The truss model
has Gumbel input, for which we (as of now) cannot construct a coherence-optimal sample. The
same holds for the structural frame model with its dependent input.

Detailed boxplots of the relative MSE against ED size for the spotlight models, showing how
each solver performed when paired with the sampling schemes, can be found in Appendix 3.C,
Figures 3.C.2–3.C.5. For the sake of readability, in this section we only show results that are
aggregated over models, separately for the low- and high-dimensional cases. For every model
and repetition index, we determine the relative ranking of the 20 (16) combinations (5 (4) solvers
× 4 sampling strategies). We also determine which of the combinations came within a factor
of {2, 5, 10} of the smallest relative MSE among this set (robustness). Then, we count how
often each combination achieved each rank, and how often each combination achieved a relative
MSE within a factor of the smallest relative MSE. The results are displayed in Figure 3.4 in
the form of stacked bar plots for the ranks, with triangle markers denoting the percentage of
robust runs. The combinations are sorted by the percentage of runs in which they achieved a
relative MSE within two times the smallest relative MSE, because we find that this metric is a
good compromise between performance and robustness. We analyze small and large ED sizes
separately, where the first half of considered ED sizes (3–4) are regarded as small and the second
half as large.

Note that, as opposed to the aggregated results in Figure 3.3, where the solvers are compared on
the same experimental designs, here the comparison is done on *different* experimental designs,
which are matched randomly.<fnref n="10" /> Our results in Figure 3.4 are bootstrapped four times using
random permutations of the replication index, corresponding in total to 250 replications, to
minimize the influence of this randomness (which is in any case not large, as can be seen from
permutation tests).

From Figure 3.4 and Figures 3.C.2–3.C.5, we make the following observations:

* There can be considerable differences in the performance of different combinations of solvers and sampling schemes. The differences are larger for low-dimensional models, visible in the spread of triangle markers in Figures 3.4a and 3.4b. For high-dimensional models, many combinations find an error that is close to the smallest error, which can be seen from the red triangle markers at high percentages in Figures 3.4c and 3.4d, and from the clustered boxplots in Figures 3.C.4 and 3.C.5, (e) and (f).

* MC and LHS perform comparably, and for the high-dimensional models almost identically.

<footnote n="10">¹⁰The relative MSE of each combination is interpreted as a random variable $E_{\text{solver}}^{\text{sampling}}$, where the randomness is induced by the randomness in the experimental design. E.g., $E_{\text{BCS}}^{\text{LHS}}$ is the random variable of relative MSE attained by BCS applied to an LHS design of specified size. The reference error ("smallest relative MSE") $E^*$ is a random variable as well, defined as the minimum over *one realization* of each combination of methods: $E^* = \min_{s \in \text{sampling}, t \in \text{solvers}} E_s^t$. The plots in Figure 3.4 are therefore read as follows: e.g., in the low-dimensional, small ED size case (Figure 3.4a): $P(E_{\text{BCS}}^{\text{D-opt}} = E^*) = 0.16$, $P(E_{\text{BCS}}^{\text{D-opt}} \le 10E^*) = 0.78$.</footnote>
<figure><img src="image_13.png" /><figcaption>(a) $d \le 10$, small ED sizes</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) $d \le 10$, large ED sizes</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) $d \ge 20$, small ED sizes</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) $d \ge 20$, large ED sizes</figcaption></figure>

**Figure 3.4:** Aggregated results for the five low-dimensional models Ishigami, undamped oscillator, borehole, damped oscillator, and wingweight (top), and for the four high-dimensional models Morris function, structural frame, two-dimensional diffusion, and 100D function (bottom). Separately for small (a),(c) and large (b),(d) experimental designs. For the low-dimensional (high-dimensional) case, we investigate five (four) solvers and four sampling schemes, resulting in 20 (16) combinations. For each model and repetition, the ranking of all the combinations is determined (note that as opposed to Figure 3.3, here the comparison is done on different EDs, which are matched randomly. Results are bootstrapped four times by random permutations to increase robustness). The stacked bar plots visualize how often each combination achieved the respective rank. The triangle markers in hues of red additionally demonstrate in how many runs the obtained relative MSE was within a factor of {2, 5, 10} of the smallest relative MSE achieved in this comparison. The combinations are sorted by the percentage of runs in which they achieved a relative MSE within two times the smallest relative MSE of the respective random pairing (red triangle marker). Plots of relative MSE against ED size for the spotlight models can be found in Appendix 3.C, Figures 3.C.2–3.C.5.
For the low-dimensional models, LHS sampling has, in most cases, median error and variability that are the same as or smaller than MC. This is consistent with the literature (Shields and Zhang, 2016). These observations are confirmed by the plots in Figure 3.4: for almost every solver, the combination with LHS is slightly better than the corresponding one with MC in each of the metrics.

* The advanced sampling schemes coh-opt and D-opt(coh-opt) show a clear advantage over MC and LHS sampling for low-dimensional models and large experimental designs (consistent with theoretical considerations and numerical experiments (Hampton and Doostan, 2015b)). For low-dimensional models and small experimental designs, they show mixed performance; for high-dimensional models, they perform the same as or worse than LHS and MC.

It is known that coh-opt sampling leads to a greater improvement over MC sampling for low-dimensional, high-degree expansions than for high-dimensional, low-degree expansions (Hampton and Doostan, 2015b; Alemazkoor and Meidani, 2018a). Note also that all numerical experiments in the literature testing coh-opt sampling were performed in $d \le 30$ dimensions, using only models with uniform input, or manufactured sparse PCE, i.e., polynomial models with an exactly sparse representation (Hampton and Doostan, 2015b,a; Alemazkoor and Meidani, 2018a; Diaz et al., 2018).

Both coh-opt and D-opt are sampling methods that aim to improve properties of the regression matrix. They are adapted to the candidate basis. If the candidate basis is large and contains many regressors that are not needed in the final sparse expansion, this adaptation might even deteriorate the solution.

* BCS is one of the best-performing solvers, almost regardless of sampling scheme. The exceptions are low-dimensional models with large experimental designs, where $SP_{LOO}$ with coh-opt sampling outperforms all other solvers. This might be related to BCS plateauing earlier than other solvers (see Figure 3.2a and 3.2b). It seems BCS is preferable whenever the information content is low (small ED sizes or high-dimensional models).

* OMP and $SP_{LOO}$ are generally quite robust (within one order of magnitude of the best solution). However, OMP often does not come close to the best solution, especially when paired with LHS or MC. BCS is more robust for high-dimensional models than for low-dimensional models. LARS and SP show mixed performance, with LARS being one of the least robust solvers.

* Aggregating the results for each sampling scheme separately (not shown), we observe that the behavior of the solvers is very similar in terms of ranking and robustness to the behavior observed on LHS (Figure 3.3), suggesting that the ranking of solvers is mostly independent of the sampling scheme.

Note that the results in Figures 3.C.2, 3.C.3, and 3.C.5 exhibit plateauing for larger sample sizes. This indicates that the maximal accuracy achievable with this set of basis functions has been reached. Using a larger basis might lead to more accurate solutions, if it contains an important regressor that was previously missing. However, note that a larger basis can also lead to less accurate solutions: when the experimental design size is held fixed while a larger basis is used, the ratio of experimental design points to basis functions is smaller, and the properties of the regression matrix might deteriorate.
### 3.3.6 Results: Comparison of sampling schemes together with solvers, using a smaller candidate basis

We repeat the experiments from the previous section for the Ishigami and borehole models, using a smaller candidate basis for which near-optimal sampling is feasible. The tested solvers are LARS, OMP, SP, $SP_{LOO}$, and BCS. We use the sampling schemes MC, LHS, coh-opt, D-opt(coh-opt) and near-opt(coh-opt). Boxplots of relative MSE against ED size are shown in Figures 3.5 and 3.6. For the sake of conciseness, we only show the combinations involving OMP and $SP_{LOO}$. The remaining plots are provided in Appendix 3.C, Figure 3.C.6.

We observe the following:

* Since the basis is smaller, the relative MSE reaches a plateau already for smaller experimental design sizes.

* Most qualitative observations regarding solver and sampling performance are the same as in the previous section, where a larger basis was used.

* Near-optimal sampling often achieves the same or a slightly smaller error than coh-opt sampling, which is consistent with (Alemazkoor and Meidani, 2018a). In many cases, near-optimal sampling achieves the smallest median error. For the Ishigami model, near-optimal sampling additionally exhibits small variability, while for the borehole model, it has a rather large spread, i.e., several outliers.

<figure><img src="image_17.png" /><figcaption>**Figure 3.5:** Results for the Ishigami model with a smaller basis ($d = 3, p = 12, q = 1$). Results for two sparse solvers and five experimental design schemes. 50 replications. For the remaining plots, see Figure 3.C.6 in Appendix 3.C.</figcaption></figure>

## 3.4 Discussion and conclusions

In this paper, we investigated sparse PCE methods with the goal of computing accurate surrogate models based on a few model evaluations.

We presented a literature survey and a framework describing the general computation procedure for sparse PCE. We have seen that the existing literature on sparse PCE can be fit into this framework and that methods developed for different components of the framework can naturally be combined.
<figure><img src="image_18.png" /><figcaption>**Figure 3.6:** Results for the borehole model with a smaller basis ($d = 8, p = 4, q = 1$). Results for two sparse solvers and five experimental design schemes. 50 replications. For the remaining plots, see Figure 3.C.6 in Appendix 3.C.</figcaption></figure>

In order to give recommendations to practitioners who want to use sparse PCE surrogates for their applications, we performed a numerical benchmark based on 11 example functions which are intended to resemble real-world engineering problems presenting different challenges. We tested several popular sparse solvers and sampling schemes on a fixed set of basis functions, using a range of experimental design sizes and 30–50 replications, and made the following observations:

* The choice of sampling scheme and sparse regression solver can make a difference of up to several orders of magnitude for the relative MSE. Mostly, the rankings of solvers and sampling schemes seem to be independent of one another: an experimental design that works best for one solver will also perform well with other solvers, and the ranking of solvers looks similar independent of which sampling scheme is used. Both solvers and sampling schemes make a greater difference for low-dimensional models.

* For low-dimensional models and small ED sizes, the solver BCS performs best, regardless of sampling scheme (with D-opt(coh-opt) being slightly preferable), while the solver SP$_{LOO}$ (a variant of SP) appears to be especially robust.

* For low-dimensional models and large ED sizes, SP$_{LOO}$ together with coherence-optimal sampling outperforms all other combinations.

* For low-dimensional models, and when the basis is small enough to make it feasible, near-optimal sampling outperforms all other sampling schemes, regardless of the solver.

* For high-dimensional models, BCS is by far the best solver. All solvers perform better when paired with LHS; in other words, no advanced sampling scheme appears competitive compared to LHS for such problems, whatever the solver used.

The benchmark results demonstrate that in costly, real-world applications it is worth choosing the sparse PCE training strategy carefully, since the methods can make a substantial difference in the quality of the resulting surrogate. While a more accurate surrogate model is generally desirable, in industrial applications it might have a higher impact for purposes such as optimization, rather than, e.g., sensitivity analysis.

Our conclusions are based on a number of benchmark models, which we consider representative of engineering models in terms of dimensionality and complexity. Naturally, however, no selection
of models can cover the whole space of engineering models. Further work would be required to
understand the connection among model properties, basis choice, experimental design size, and
sparse PCE techniques like solvers and sampling schemes.

All results were obtained using a fixed basis based on a heuristic choice (see Section 3.3.3). Generally, when the optimal degree of the basis is unknown, degree adaptivity (based on a cross-validation error) can be a useful strategy. Due to time and space constraints, this was not investigated in the present work. Adaptivity critically depends on the availability of an accurate error estimator. Some of the best solvers in this study (i.e., OMP and BCS) tend to underestimate the generalization error (not shown in the plots), which might be a drawback in the setting of adaptive degree selection and might change the effect and ranking of solvers and sampling methods. For a detailed discussion and benchmark of basis-adaptive schemes, we refer the reader to Lüthen et al. (2022a).

As evident from the extensive literature on the topic, sparse PCE is an already well-established
technique, as well as an active field of methodological research. Recent innovations include
Bayesian techniques for sparse PCE and the identification of suitable rotated coordinates for
the expansion. Such innovative ideas are expected to lead to further improvements in the
computation of sparse PCE, which will in turn benefit applications as well as all advanced
methods that use sparse PCE as one of their building blocks (see, e.g., Schöbi et al. (2015);
Chatterjee et al. (2019); Zhu and Sudret (2020); Marelli et al. (2021b)).

PCE is a popular metamodelling tool in the engineering community, and many different methods
are available. Up to now, the choice of which of the many PCE methods to apply was mostly
left to chance or the personal experience of the practitioner. In our benchmark, we explored
a significant set of methods that have received attention in the past few years. We hope that
this work can serve as a basis for further benchmarking efforts, in order to identify which of
the many available methods are most suitable for real-world problems. These might include
sequential enrichment of experimental design, Gaussian adaptation of the input space, stepwise
regression algorithms and many other ideas for sparse solvers, as well as methods for extremely
high-dimensional problems.

Our benchmark code is available on request. The solvers BCS and SP<sub>LOO</sub> will be made available in the 1.4.0 release of UQLab. For a description of how to add custom sampling schemes and sparse solvers for PCE to UQLab, we refer the reader to the supplementary material accompanying this paper. To facilitate easier benchmarking of PCE techniques on a large number of examples in a standardized setup, we are actively engaging in designing and developing a benchmarking platform for surrogate modelling methods similar to the UCI machine learning repository<sup>11</sup> or the structural reliability platform RPrepo<sup>12</sup> where data sets, models, and methods can be made available for testing and benchmarking.

Acknowledgments

We thank Dr. Emiliano Torre for helpful discussions.

<footnote n="11">¹¹https://archive.ics.uci.edu</footnote>

<footnote n="12">¹²https://rprepo.readthedocs.io/en/latest/reliability_problems.html</footnote>
References

Abraham, S., Raisee, M., Ghorbaniasl, G., Contino, F., and Lacor, F. (2017). A robust and efficient stepwise regression method for building sparse polynomial chaos expansions. *J. Comput. Phys.*, 332:461–474.

Alemazkoor, N. and Meidani, H. (2017). Divide and conquer: An incremental sparsity promoting compressive sampling approach for polynomial chaos expansions. *Comput. Methods Appl. Mech. Engrg.*, 318:937–956.

Alemazkoor, N. and Meidani, H. (2018a). A near-optimal sampling strategy for sparse recovery of polynomial chaos expansions. *J. Comput. Phys.*, 371:137–151.

Alemazkoor, N. and Meidani, H. (2018b). A preconditioning approach for improved estimation of sparse polynomial chaos expansions. *Comput. Methods Appl. Mech. Engrg.*, 342:474–489.

Arjoune, Y., Kaabouch, N., El Ghazi, H., and Tamtaoui, A. (2017). Compressive sensing: Performance comparison of sparse recovery algorithms. In *2017 IEEE CCWC*, pages 1–7. IEEE.

Babacan, D. (2011). MATLAB code for the TIP 2010 paper "Bayesian Compressive Sensing using Laplace Priors". http://www.dbabacan.info/software.html. [Online; accessed 28-August-2019].

Babacan, S. D., Molina, R., and Katsaggelos, A. K. (2010). Bayesian compressive sensing using Laplace priors. *IEEE Trans. Image Process.*, 19(1):53–63.

Baptista, R., Stolbunov, V., and Nair, P. B. (2019). Some greedy algorithms for sparse polynomial chaos expansions. *J. Comput. Phys.*, 387:303–325.

Baumann, D. and Baumann, K. (2014). Reliable estimation of prediction errors for QSAR models under model uncertainty using double cross-validation. *J. Cheminf.*, 6(1):47.

Berchier, M. (2015). Orthogonal matching pursuit for sparse polynomial chaos expansions. ETH Zürich. Semester project.

Berveiller, M. (2005). *Stochastic finite elements : intrusive and non intrusive methods for reliability analysis*. PhD thesis, Université Blaise Pascal, Clermont-Ferrand.

Berveiller, M., Sudret, B., and Lemaire, M. (2006). Stochastic finite elements: a non intrusive approach by regression. *Eur. J. Comput. Mech.*, 15(1–3):81–92.

Bhattacharyya, B. (2020). Global sensitivity analysis: A Bayesian learning based polynomial chaos approach. *J. Comput. Phys.*, 415:109539.

Birgin, E. G., Martínez, J. M., and Raydan, M. (2000). Nonmonotone spectral projected gradient methods on convex sets. *SIAM J. Optim.*, 10(4):1196–1211.

Blatman, G. and Sudret, B. (2008). Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach. *C. R. Mécanique*, 336(6):518–523.
Blatman, G. and Sudret, B. (2010a). An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. *Prob. Eng. Mech.*, 25:183-197.

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

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

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

Candès, E. and Romberg, J. (2005). *ℓ₁-MAGIC*. https://statweb.stanford.edu/~candes/software/11magic/. [Online; accessed 22-January-2020].

Candès, E., Romberg, J., and Tao, T. (2006). Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. *IEEE Trans. Inform. Theory*, 52(2):489-509.

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

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

Candès, E. J., Wakin, M. B., and Boyd, S. P. (2008). Enhancing sparsity by reweighted $ℓ_1$ minimization. *J. Fourier Anal. Appl.*, 14(5-6):877-905.

Carron, I. (2013). *Compressive Sensing: The Big Picture*. https://sites.google.com/site/igorcarron2/cs. [Online; accessed 22-January-2020].

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

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

Cheng, K. and Lu, Z. (2018a). Adaptive sparse polynomial chaos expansions for global sensitivity analysis based on support vector regression. *Comput. Struct.*, 194:86-96.

Cheng, K. and Lu, Z. (2018b). Sparse polynomial chaos expansion based on D-MORPH regression. *Appl. Math. Comput.*, 323:17-30.

Cohen, A. and Migliorati, G. (2017). Optimal weighted least-squares methods. *SMAI J. Comp. Math.*, 3:181-203.

Cook, R. D. and Nachtsheim, C. J. (1980). A comparison of algorithms for constructing exact D-optimal designs. *Technometrics*, 22(3):315-324.
Dai, W. and Milenkovic, O. (2009). Subspace pursuit for compressive sensing signal reconstruc-
tion. *IEEE Trans. Inform. Theory*, 55(5):2230–2249.

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

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

Donoho, D., Drori, I., Stodden, V., Tsaig, Y., and Shahram, M. (2007). SparseLab - Seeking Sparse Solutions to Linear Systems of Equations. http://sparselab.stanford.edu/. [Online; accessed 22-January-2020].

Donoho, D. L. (2006). Compressed sensing. *IEEE Trans. Inform. Theory*, 52(4):1289–1306.

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

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

Dutta, S. and Gandomi, A. H. (2020). Design of experiments for uncertainty quantification based
on polynomial chaos expansion metamodels. In *Handbook of Probabilistic Models*, pages 369–
381. Elsevier.

Dykstra, O. (1971). The augmentation of experimental data to maximize [$X'X$]. *Technometrics*, 13(3):682–688.

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

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

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

Fajraoui, N., Marelli, S., and Sudret, B. (2017). Sequential design of experiment for sparse polynomial chaos expansions. *SIAM/ASA J. Uncertain. Quantif.*, 5(1):1061–1085.

Faul, A. C. and Tipping, M. E. (2002). Analysis of sparse Bayesian learning. In *Advances in neural information processing systems*, pages 383–389.

Fedorov, V. V. (2013). *Theory of optimal experiments*. Elsevier.

Figueiredo, M. A. (2003). Adaptive sparseness for supervised learning. *IEEE Trans. Pattern Anal. Mach. Intell.*, 25(9):1150–1159.

Figueiredo, M. A. and Nowak, R. D. (2001). Wavelet-based image estimation: an empirical Bayes approach using Jeffrey's noninformative prior. *IEEE Trans. Image Process.*, 10(9):1322–1331.
Forrester, A., Sobester, A., and Keane, A. (2008). *Engineering design via surrogate modelling: a practical guide*. Wiley.

Ghanem, R. G. and Spanos, P. (1991). *Stochastic finite elements – A spectral approach*. Springer Verlag, New York. (Reedited by Dover Publications, Mineola, 2003).

Gu, M. and Eisenstat, S. C. (1996). Efficient algorithms for computing a strong rank-revealing QR factorization. *SIAM J. Sci. Comput.*, 17(4):848–869.

Guo, L., Narayan, A., Liu, Y., and Zhou, T. (2020). Sparse approximation of data-driven polynomial chaos expansions: An induced sampling approach. *Commun. Math. Res.*, 36(2):128–153.

Guo, L., Narayan, A., Zhou, T., and Chen, Y. (2017). Stochastic collocation methods via $ℓ_1$ minimization using randomized quadratures. *SIAM J. Sci. Comput.*, 39(1):A333–A359.

Hadigol, M. and Doostan, A. (2018). Least squares polynomial chaos expansion: A review of sampling strategies. *Comput. Methods Appl. Mech. Engrg.*, 332:382–407.

Halton, J. H. (1960). On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. *Numer. Math.*, 2(1):84–90.

Hampton, J. and Doostan, A. (2015a). Coherence motivated sampling and convergence analysis of least squares polynomial chaos regression. *Comput. Methods Appl. Mech. Engrg.*, 290:73–97.

Hampton, J. and Doostan, A. (2015b). Compressive sampling of polynomial chaos expansions: Convergence analysis and sampling strategies. *J. Comput. Phys.*, 280:363–386.

Hampton, J. and Doostan, A. (2017). COH-OPT. https://github.com/CU-UQ/COH-OPT. [Online; accessed 22-January-2020].

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

Harper, W. V. and Gupta, S. K. (1983). Sensitivity/uncertainty analysis of a borehole scenario comparing Latin hypercube sampling and deterministic sensitivity approaches. Technical Report No. BMI/ONWI-516, Battelle Memorial Inst. - Office of Nuclear Waste Isolation, Columbus, OH (USA).

Hastie, T., Tibshirani, R., and Friedman, J. (2001). *The elements of statistical learning: Data mining, inference and prediction*. Springer, New York.

Hong, Y. P. and Pan, C.-T. (1992). Rank-revealing QR factorizations and the singular value decomposition. *Math. Comput.*, 58(197):213–232.

Hosder, S., Walters, R., and Balch, M. (2007). Efficient sampling for non-intrusive polynomial chaos applications with multiple uncertain input variables. In *48th AIAA/ASME/ASCE/AHS/ASC Struct. Struct. Dyn. Mater. Conf.*, page 1939.

Hu, C. and Youn, B. D. (2011). Adaptive-sparse polynomial chaos expansion for reliability analysis and design of complex engineering systems. *Struct. Multidisc. Optim.*, 43:419–442.
Hu, R. and Ludkovski, M. (2017). Sequential design for ranking response surfaces. *SIAM/ASA J. Uncertain. Quantif.*, 5(1):212–239.

Huan, X., Safta, C., Sargsyan, K., Vane, Z. P., Lacaze, G., Oefelein, J. C., and Najm, H. N. (2018). Compressive sensing with cross-validation and stop-sampling for sparse polynomial chaos expansions. *SIAM/ASA J. Uncertain. Quantif.*, 6(2):907–936.

Isukapalli, S. S. (1999). *Uncertainty analysis of transport-transformation models*. PhD thesis, The State University of New Jersey.

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

Jakeman, J. D., Narayan, A., and Zhou, T. (2017). A generalized sampling and preconditioning scheme for sparse approximation of polynomial chaos expansions. *SIAM J. Sci. Comput.*, 39(3):A1114–A1144.

Ji, S., Xue, Y., and Carin, L. (2008). Bayesian compressive sensing. *IEEE Trans. Signal Process.*, 56(6):2346–2356.

Kiefer, J. and Wolfowitz, J. (1959). Optimum designs in regression problems. *Ann. Math. Stat.*, 30(2):271–294.

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

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

Lataniotis, C. (2019). *Data-driven uncertainty quantification for high-dimensional engineering problems*. PhD thesis, ETH Zurich.

Li, C. C. and Der Kiureghian, A. (1993). Optimal discretization of random fields. *J. Eng. Mech.*, 119(6):1136–1154.

Li, G. and Rabitz, H. (2010). D-MORPH regression: application to modeling with unknown parameters more than observation data. *J. Math. Chem.*, 48(4):1010–1035.

Liu, Z., Lesselier, D., Sudret, B., and Wiart, J. (2020a). Surrogate modeling based on resampled polynomial chaos expansions. *Reliab. Eng. Sys. Safety*, 202:107008.

Liu, Z., Lesselier, D., Sudret, B., and Wiart, J. (2020b). Surrogate modeling of indoor down-link human exposure based on sparse polynomial chaos expansion. *Int. J. Uncertainty Quantification*, 10(2):145–163.

Lüthen, N., Marelli, S., and Sudret, B. (2022). Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications. *Int. J. Uncertainty Quantification*, 12(3):49–74.
Marelli, S. and Sudret, B. (2014). UQLab: A framework for uncertainty quantification in Matlab. In *Vulnerability, Uncertainty, and Risk (Proc. 2nd Int. Conf. on Vulnerability, Risk Analysis and Management (ICVRAM2014), Liverpool, United Kingdom)*, pages 2554–2563.

Marelli, S. and Sudret, B. (2019). UQLab user manual – Polynomial chaos expansions. Technical report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich. Report# UQLab-V1.3-104.

Marelli, S., Wagner, P.-R., Lataniotis, C., and Sudret, B. (2021). Stochastic spectral embedding. *Int. J. Uncertainty Quantification*, 11(2):25–47.

Mathelin, L. and Gallivan, K. (2012). A compressed sensing approach for partial differential equations with random input data. *Commun. Comput. Phys.*, 12(4):919–954.

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

Mikhalev, A. and Oseledets, I. V. (2018). Rectangular maximum-volume submatrices and their applications. *Linear Algebra Appl.*, 538:187–211.

Montgomery, D. C. (2004). *Design and analysis of experiments*. John Wiley and Sons, New York.

Narayan, A., Jakeman, J., and Zhou, T. (2017). A Christoffel function weighted least squares algorithm for collocation approximations. *Math. Comput.*, 86(306):1913–1947.

Needell, D. and Tropp, J. A. (2009). CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. *Appl. Comput. Harmon. A.*, 26(3):301–321.

Nguyen, N.-K. and Miller, A. J. (1992). A review of some exchange algorithms for constructing discrete D-optimal designs. *Comput. Stat. Data An.*, 14(4):489–498.

Owen, A. B. (1994). Controlling correlations in Latin hypercube samples. *J. Am. Stat. Assoc.*, 89(428):1517–1522.

Papaioannou, I., Ehre, M., and Straub, D. (2019). PLS-based adaptation for efficient PCE representation in high dimensions. *J. Comput. Phys.*, 387:186–204.

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

Peng, J., Hampton, J., and Doostan, A. (2014). A weighted $l_1$-minimization approach for sparse polynomial chaos expansions. *J. Comput. Phys.*, 267:92–111.

Perkó, Z., Gilli, L., Lathouwers, D., and Kloosterman, J. L. (2014). Grid and basis adaptive polynomial chaos techniques for sensitivity and uncertainty analysis. *J. Comput. Phys.*, 260:54–84.
Pronzato, L. and Müller, W. G. (2012). Design of computer experiments: space filling and beyond. *Stat. Comput.*, 22(3):681–701.

Qaisar, S., Bilal, R. M., Iqbal, W., Naureen, M., and Lee, S. (2013). Compressive sensing: From theory to applications, a survey. *J. Commun. Netw.*, 15(5):443–456.

Rauhut, H. and Ward, R. (2012). Sparse Legendre expansions via $ℓ_1$-minimization. *J. Approx. Theory*, 164(5):517–533.

Sargsyan, K., Safta, C., Najm, H., Debusschere, B., Ricciuto, D., and Thornton, P. (2014). Dimensionality reduction for complex models via Bayesian compressive sensing. *Int. J. Uncertainty Quantification*, 4(1):63–93.

Schöbi, R., Sudret, B., and Wiart, J. (2015). Polynomial-chaos-based Kriging. *Int. J. Uncertainty Quantification*, 5(2):171–193.

Seeger, M. W. and Nickisch, H. (2008). Compressed sensing and Bayesian experimental design. In *Proc. of the 25th Int. Conf. on Machine Learning*, pages 912–919. ACM.

Seshadri, P., Narayan, A., and Mahadevan, S. (2017). Effectively subsampled quadratures for least squares polynomial approximations. *SIAM/ASA J. Uncertain. Quantif.*, 5(1):1003–1023.

Settles, B. (2012). Active learning. *Synthesis Lectures on Artificial Intelligence and Machine Learning*, 6(1):1–114.

Shao, Q., Younes, A., Fahs, M., and Mara, T. (2017). Bayesian sparse polynomial chaos expansion for global sensitivity analysis. *Comput. Meth. Appl. Mech. Eng.*, 318:474–496.

Shen, D., Wu, H., Xia, B., and Gan, D. (2020). Polynomial chaos expansion for parametric problems in engineering systems: A review. *IEEE Syst. J.*

Shields, M. D. and Zhang, J. (2016). The generalization of Latin hypercube sampling. *Reliab. Eng. Syst. Saf.*, 148:96–108.

Shin, Y. and Xiu, D. (2016a). Nonadaptive quasi-optimal points selection for least squares linear regression. *SIAM J. Sci. Comput.*, 38(1):A385–A411.

Shin, Y. and Xiu, D. (2016b). On a near optimal sampling strategy for least squares polynomial regression. *J. Comput. Phys.*, 326:931–946.

Sobol', I. M. (1967). Distribution of points in a cube and approximate evaluation of integrals. *U.S.S.R Comput. Maths. Math. Phys.*, 7:86–112.

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

Tang, G. and Iaccarino, G. (2014). Subsampled Gauss quadrature nodes for estimating polynomial chaos expansions. *SIAM/ASA J. Uncertain. Quantif.*, 2(1):423–443.

Tarakanov, A. and Elsheikh, A. H. (2019). Regression-based sparse polynomial chaos for uncertainty quantification of subsurface flow models. *J. Comput. Phys.*, 399:108909.
Tipireddy, R. and Ghanem, R. (2014). Basis adaptation in homogeneous chaos spaces. *J. Comput. Phys.*, 259:304–317.

Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. *J. Mach. Learn. Res.*, 1(Jun):211–244.

Tipping, M. E. and Faul, A. C. (2003). Fast marginal likelihood maximisation for sparse Bayesian models. In *Proc. 9th Int. Workshop on Artificial Intelligence and Statistics*.

Tropp, J. A. and Gilbert, A. C. (2007). Signal recovery from random measurements via orthogonal matching pursuit. *IEEE Trans. Inform. Theory*, 53(12):4655–4666.

Tsilifis, P., Huan, X., Safta, C., Sargsyan, K., Lacaze, G., Oefelein, J. C., Najm, H. N., and Ghanem, R. G. (2019). Compressive sensing adaptation for polynomial chaos expansions. *J. Comput. Phys.*, 380:29–47.

Tsilifis, P., Papaioannou, I., Straub, D., and Nobile, F. (2020). Sparse Polynomial Chaos expansions using variational relevance vector machines. *J. Comput. Phys.*, 416:109498.

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

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

Wipf, D. P., Palmer, J., and Rao, B. D. (2004). Perspectives on sparse Bayesian learning. In *Advances in neural information processing systems*, pages 249–256.

Wipf, D. P. and Rao, B. D. (2004). Sparse Bayesian learning for basis selection. *IEEE Trans. Signal Process.*, 52(8):2153–2164.

Xiu, D. and Hesthaven, J. S. (2005). High-order collocation methods for differential equations with random inputs. *SIAM J. Sci. Comput.*, 27(3):1118–1139.

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

Yan, L., Guo, L., and Xiu, D. (2012). Stochastic collocation algorithms using $\ell_1$-minimization. *Int. J. Uncertainty Quantification*, 2(3).

Yang, X. and Karniadakis, G. E. (2013). Reweighted $\ell_1$ minimization method for stochastic elliptic differential equations. *J. Comput. Phys.*, 248:87–108.

Yang, X., Li, W., and Tartakovsky, A. (2018). Sliced-inverse-regression-aided rotated compressive sensing method for uncertainty quantification. *SIAM/ASA J. Uncertain. Quantif.*, 6(4):1532–1554.

Yin, P., Lou, Y., He, Q., and Xin, J. (2015). Minimization of $\ell_{1-2}$ for compressed sensing. *SIAM J. Sci. Comput.*, 37(1):A536–A563.
Zankin, V. P., Ryzhakov, G. V., and Oseledets, I. V. (2018). Gradient descent-based D-optimal design for the least-squares polynomial approximation. *arXiv preprint arXiv:1806.06631*.

Zein, S., Colson, B., and Glineur, F. (2013). An efficient sampling method for regression-based polynomial chaos expansion. *Commun. Comput. Phys.*, 13(4):1173–1188.

Zhang, Z., Xu, Y., Yang, J., Li, X., and Zhang, D. (2015). A survey of sparse representation: algorithms and applications. *IEEE Access*, 3:490–530.

Zhao, H., Gao, Z., Xu, F., Zhang, Y., and Huang, J. (2019). An efficient adaptive forward-backward selection method for sparse polynomial chaos expansion. *Comput. Methods Appl. Mech. Engrg.*, 355:456–491.

Zhao, W. and Bu, L. (2019). Global sensitivity analysis with a hierarchical sparse metamodeling method. *Mech. Syst. Signal Pr.*, 115:769–781.

Zhou, T., Narayan, A., and Xu, Z. (2014). Multivariate discrete least-squares approximations with a new type of collocation grid. *SIAM J. Sci. Comput.*, 36(5):A2401–A2422.

Zhou, Y., Lu, Z., and Cheng, K. (2019a). Sparse polynomial chaos expansions for global sensitivity analysis with partial least squares and distance correlation. *Struct. Multidiscip. O.*, 59(1):229–247.

Zhou, Y., Lu, Z., Cheng, K., and Ling, C. (2019b). An efficient and robust adaptive sampling method for polynomial chaos expansion in sparse Bayesian learning framework. *Comput. Methods Appl. Mech. Engrg.*, 352:654–674.

Zhou, Y., Lu, Z., Cheng, K., and Shi, Y. (2019c). An expanded sparse Bayesian learning method for polynomial chaos expansion. *Mech. Syst. Signal Pr.*, 128:153–171.

Zhou, Y., Lu, Z., Hu, J., and Hu, Y. (2020). Surrogate modeling of high-dimensional problems via data-driven polynomial chaos expansions and sparse partial least square. *Comput. Methods Appl. Mech. Engrg.*, 364:112906.

Zhu, X. and Sudret, B. (2020). Replication-based emulation of the response distribution of stochastic simulators using generalized lambda distributions. *Int. J. Uncertainty Quantification*, 10(3):249–275.
# Appendix

## 3.A Details on experimental design sampling techniques

It depends on certain properties of the regression matrix $\Psi$ whether or not sparse regression techniques are able to find the true sparse solution of a linear system of equations (assuming that it exists). In the context of polynomial chaos, the entries of the regression matrix are the basis polynomials evaluated at the design points. The basis polynomials are determined by the distribution of the input random variables and the choice of the index set $A$, while the design points $\{x^{(j)}\}_{j=1}^N$ can be chosen freely from the input space to optimize properties of the resulting regression matrix.

In the following, we present sampling schemes that have been proposed in the literature for the computation of sparse PCE. Some of the schemes come with theoretical results about their performance for sparse PCE, others have heuristic justification or have guarantees for least-squares regression. They can be broadly grouped into three categories:

* Sampling according to the input distribution

- MC (Doostan and Owhadi, 2011; Hampton and Doostan, 2015b)

- LHS (McKay et al., 1979)

* Sampling from a modified distribution (induced sampling)

- asymptotic (Hampton and Doostan, 2015b)

- coherence-optimal (Hampton and Doostan, 2015b)

- Christoffel sparse approximation (Jakeman et al., 2017)

* Optimizing matrix properties

- D-optimal (Diaz et al., 2018)

- S-optimal (Shin and Xiu, 2016a; Fajraoui et al., 2017)

- near-optimal (Alemazkoor and Meidani, 2018a)

Some of the sampling schemes are nontrivial or costly to evaluate, or even not available for all input distributions. However, the bottleneck in surrogate modelling for practical applications is typically the repeated evaluation of the model, which justifies the use of a complex sampling scheme if it allows better approximation with fewer samples.
### 3.A.1 Sampling according to the input distribution

This class of sampling methods consists of all methods that are oblivious to the choice of truncation set $\mathcal{A}$ and whose main objective is to distribute design points evenly in the quantile space. Heuristically, the more uniformly the points are distributed in the quantile space, the more information about the model is captured in the model evaluations, since no region of the input domain is forgotten.

LHS is one technique for achieving a space-filling design. For each component of the input random vector $\mathbf{X}$, the corresponding quantile space is divided into $N$ intervals. In each interval, one point is sampled uniformly at random. Then, the points for each dimension are combined randomly into vectors and finally transformed into the input space using an isoprobabilistic transform. LHS can be shown to reduce the variance of linear regression estimates when the main effects are dominant, i.e., when the most important terms have interaction order one (Shields and Zhang, 2016). LHS can be combined with heuristic criteria such as the maximin distance strategy, where among several random LHS designs the one with the largest minimal pairwise distance between points is chosen, to further improve on the space-filling property.

Stratified sampling is a related sampling technique in which the input space is divided into disjoint regions, called strata, from which points are sampled and weighted according to the probability mass of their stratum. Stratified sampling reduces the variance of statistical estimators (McKay et al., 1979). There exists a range of methods between stratified sampling and LHS, called partially stratified sampling, which are able to reduce the variance of statistical estimators when interaction terms are dominant (Shields and Zhang, 2016). The authors propose an additional method called Latinized partially stratified sampling (LPSS) which combines LHS and stratified sampling with the aim of minimizing the variance of the resulting estimator. It is especially beneficial when there is prior knowledge about which variable groups interact, and it has been used for several problems with input dimension $d = 100$.

MC sampling, i.e., sampling from the input distribution, is a special case of the coherence-based theory detailed in Section 3.A.2.1 below, and bounds on the coherence and the associated number of points needed for recovery can be derived (Doostan and Owhadi, 2011; Rauhut and Ward, 2012; Yan et al., 2012; Hampton and Doostan, 2015b).

### 3.A.2 Sampling from a different distribution

The ability of sparse regression to recover the true sparse solution (if it exists; otherwise it recovers the best sparse approximation) largely depends on the regression matrix. In the case of PCE, the entries of this matrix are the evaluations of the basis polynomials at the experimental design points. The points can be chosen in a way that improves the recovery properties of the matrix.

Several approaches exist in which the $l^1$-minimization problem is modified into a weighted problem and samples are drawn not from the input distribution, but from a suitable modified distribution. The idea of these approaches is as follows. Define a weight function $w(x) : \Omega \to \mathbb{R}$ in a suitable way, which will be explained later. For an ED $\{\mathbf{x}^{(k)}\}_{k=1}^N$, define the diagonal matrix
$W = \text{diag}(w(\boldsymbol{x}^{(1)}), \dots, w(\boldsymbol{x}^{(N)}))$. Then the following modified system is solved:

$$ \min_{c} \|c\|_{1} \quad \text{s.t.} \quad \|W\Psi c - Wy\|_{2} \le \epsilon. \qquad (3.14) $$

Depending on $w(\boldsymbol{x})$, this modification can improve or deteriorate the solution $c$. Of course, the weight function is chosen to improve the solution. The matrix $W\Psi$ can also be interpreted as the evaluation of a modified basis $\{\tilde{\psi}_{\alpha}(\boldsymbol{x}) = w(\boldsymbol{x})\psi_{\alpha}(\boldsymbol{x})\}_{\alpha \in \mathcal{A}}$. To ensure orthonormality of the columns of $W\Psi$, the design points are drawn from a suitably modified input distribution $f_{\tilde{X}}$.

### 3.A.2.1 Coherence, isotropy, and weighted orthonormal systems

In this section, we define concepts that are the basis for guarantees on accuracy and stability for different sampling distributions. We mainly follow the exposition in Hampton and Doostan (2015b).

In the setting of PCE, the coherence of an orthonormal system $\{\psi_{\alpha}\}_{\alpha \in \mathcal{A}}$ is defined by

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

For distributions for which this quantity would be $\infty$, such as a Gaussian distribution, see the remark below.

A second important concept is isotropy (Candès and Plan, 2011): a random matrix, whose rows are chosen randomly following some distribution $a \sim F_a$, is isotropic if it holds that $E[a^T a] = 1$. In the case of PCE, $F_a$ is induced by propagating the input distribution $F_X$ through the basis functions. By construction, the regression matrix of standard PCE is isotropic if the ED is sampled from the input distribution. Under the assumption that the regression matrix $\Psi$ is isotropic, the number of samples needed for perfect recovery of sparse solutions in the noiseless case is proportional to $\mu(\mathcal{A}, \{\psi_{\alpha}\})s \log(P)$ with high probability (Candès and Plan, 2011), where $s$ is the sparsity of the solution vector and $P = |\mathcal{A}|$ is the number of basis functions. A similar result holds in the noisy case.

Thus, an orthonormal system $\{\psi_{\alpha}\}_{\alpha \in \mathcal{A}}$ with low coherence $\mu(\mathcal{A}, \{\psi_{\alpha}\})$ requires fewer samples for perfect recovery. The goal of coherence-optimal sampling is to find a weighted system $\{\tilde{\psi}_{\alpha}(\boldsymbol{x}) = w(\boldsymbol{x})\psi_{\alpha}(\boldsymbol{x})\}_{\alpha \in \mathcal{A}}$ that achieves $\mu(\mathcal{A}, \{\tilde{\psi}_{\alpha}\}) < \mu(\mathcal{A}, \{\psi_{\alpha}\})$ and is orthonormal with respect to some distribution $\tilde{f}_X$.

The ideas of isotropy and coherence were applied to PCE by Hampton and Doostan (2015b), who construct an isotropic regression matrix with improved coherence as follows. Let $B: \mathcal{D} \to \mathbb{R}$ be the tight upper bound for the polynomial basis,

$$ B(\boldsymbol{x}) = \max_{\alpha \in \mathcal{A}} |\psi_{\alpha}(\boldsymbol{x})|. \qquad (3.16) $$

Let $G: \mathcal{D} \to \mathbb{R}$ be a loose upper bound with $G(\boldsymbol{x}) \ge B(\boldsymbol{x}) \quad \forall \boldsymbol{x} \in \mathcal{D}$. (G is useful because using a simple expression for the upper bound can in some cases result in $\tilde{f}_X$ being a well-known distribution that can be sampled from easily.) Define a new probability distribution $\tilde{f}_X(\boldsymbol{x})$ by

$$ \tilde{f}_X(\boldsymbol{x}) = c^2 G(\boldsymbol{x})^2 f_X(\boldsymbol{x}), \qquad (3.17) $$
where $c = (\int_{\Omega} f_X(x)G(x)^2dx)^{-\frac{1}{2}}$ is the normalizing constant. Then, with the weight function

$$w(x) = \frac{1}{cG(x)}, \qquad (3.18)$$

the set of functions {$\tilde{\psi}_{\alpha}(x) = w(x)\psi_{\alpha}(x)}_{\alpha \in \mathcal{A}}$ is an orthonormal system with respect to the distribution $\tilde{f}_{\mathcal{X}}$. This follows directly from the orthonormality of {$\psi_{\alpha}}_{\alpha \in \mathcal{A}}$ with respect to $f_{\mathcal{X}}$. Furthermore, if $G = B$, the coherence $\mu(\mathcal{A}, \{\tilde{\psi}_{\alpha}\})$ is minimal.

**Remark** Some polynomial bases (e.g. Hermite polynomials) do not have a finite upper bound. It is still possible to obtain similar results by considering a smaller domain $S \subset \mathcal{D}$ on which the upper bound is finite and the isotropy is still approximately fulfilled. The modified probability distribution is then $\tilde{f}_{\mathcal{X}}(x) = c^2G(x)^2f_{\mathcal{X}}(x)1_S(x)$.

### 3.A.2.2 Sampling using a loose upper bound ("asymptotic sampling")

In the case of Legendre and Hermite polynomials, and using a certain loose upper bound $G(x) \ge \max_{\alpha \in \mathcal{A}} |\psi_{\alpha}(x)|$, analytical expressions for distributions with improved coherence can be obtained (Hampton and Doostan, 2015b).

In the case of Legendre polynomials on $[-1, 1]^d$, a loose upper bound on the polynomials is given by $G(x) \propto \prod_{i=1}^d (1 - x_i^2)^{-\frac{1}{4}}$, which leads to the Chebyshev distribution $\tilde{f}_{\mathcal{X}}(x) = \prod_{i=1}^d \frac{1}{\pi\sqrt{1-x_i^2}}$ and to the weight function $w(x) = \prod_{i=1}^d (1 - x_i^2)^{\frac{1}{4}}$.

In the case of Hermite polynomials for standard Gaussian variables, a loose upper bound on the polynomials is given by $G(x) \propto \exp(\frac{1}{4} \|x\|_2^2)$, and the subset $S$ is chosen to be the $d$-dimensional ball with radius $\sqrt{2\sqrt{2p+1}}$. This leads to a uniform distribution $\tilde{f}_{\mathcal{X}}$ on $S$ and to the weight function $w(x) = \exp(-\frac{1}{4} \|x\|_2^2)$.

Additionally, asymptotic distributions for Laguerre polynomials (corresponding to the Gamma distribution) and for Jacobi polynomials (Beta distribution) have been implemented in the software package COH-OPT (Hampton and Doostan, 2017).

For Legendre polynomials, asymptotic sampling has a smaller coherence than standard sampling in the case $d < p$ (asymptotically). In the case $d > p$, which is more common in applications, standard sampling has (asymptotically) a smaller coherence. According to theory, the sampling scheme with smaller coherence should exhibit better recovery rates. This is confirmed numerically (Hampton and Doostan, 2015b, section 5.1). For Hermite polynomials, the same observation is made.

### 3.A.2.3 Coherence-optimal sampling

The choice $G = B$ leads to the minimum possible coherence $\mu(\mathcal{A}, \{\tilde{\psi}_{\alpha}\})$ (Hampton and Doostan, 2015b, Theorem 4.5). $B$ is simple to evaluate for a single point $x \in \mathcal{D}$, but its functional form is in general not known. Therefore, Hampton and Doostan (2015b) suggest sampling $\tilde{f}_{\mathcal{X}} \propto B^2 f_{\mathcal{X}}$ using Markov chain Monte Carlo (MCMC) sampling with proposal distribution equal to the input distribution in the case $d \ge p$ and equal to the asymptotic distribution in the case $d < p$.
The resulting (unnormalized) weights are $w(\boldsymbol{x}) = \frac{1}{B(\boldsymbol{x})}$. As expected from theory, numerical examples indicate that coherence-optimal sampling achieves better recovery and a smaller error in various norms than both standard and asymptotic sampling (Hampton and Doostan, 2015b). Coherence-optimal sampling can be shown to have good properties also when used as a sampling scheme for least-squares regression (Hampton and Doostan, 2015a).

MATLAB code for MCMC-based coherence-optimal sampling is available (Hampton and Doostan, 2015a,b). However, MCMC-based coherence-optimal sampling can be very slow for high-dimensional input. An alternative is rejection-based coherence-optimal sampling. Here, samples $\boldsymbol{x}_{cand}$ are generated from a proposal distribution $f_{prop}$, which has the property that there is a $\gamma \in \mathbb{R}$ such that $\gamma f_{prop}(\boldsymbol{x}) \ge \tilde{f}_{X}(\boldsymbol{x})$ for all $\boldsymbol{x} \in \mathcal{D}$. Uniform random numbers $u \sim \text{i.i.d. } \mathcal{U}([0, 1])$ are generated. A proposed point $\boldsymbol{x}_{cand}$ is accepted if $u \le \frac{\tilde{f}_{X}(\boldsymbol{x}_{cand})}{\gamma f_{prop}(\boldsymbol{x}_{cand})}$. This is the implementation used in this benchmark. We use a product proposal density whose marginals are determined by the input marginals, the dimension of the problem, and the degree of the expansion: we choose uniform proposal marginals for uniform input marginals. For Gaussian input marginals, we use Gaussian proposal marginals if $d \ge p$; otherwise, we use the corresponding asymptotic distribution. As usual, lognormal input is mapped to Gaussian random variables before sampling (Blatman and Sudret, 2011).

Note that for Gaussian input, coherence-optimal and asymptotic sampling have a significantly larger spread than input sampling, as can be seen from Figure 3.A.2. Their support is the ball of radius $r = \sqrt{2\sqrt{2p} + 2}$ (as implemented in (Hampton and Doostan, 2017)). This can potentially cause problems in engineering applications, for which simulations may be less accurate when the input parameters are far from typical operating conditions.

### 3.A.2.4 Christoffel sparse approximation

A similar weighted sampling scheme is Christoffel sparse approximation (Narayan et al., 2017; Jakeman et al., 2017; Cohen and Migliorati, 2017). Those authors propose to use the weight function

$$w(\boldsymbol{x}) = \left( \frac{1}{|\mathcal{A}|} \sum_{\boldsymbol{\alpha} \in \mathcal{A}} |\psi_{\boldsymbol{\alpha}}(\boldsymbol{x})|^2 \right)^{-\frac{1}{2}} \quad (3.19)$$

which leads to a modified basis that has pointwise minimal average squared basis magnitude (compare to (3.18) with $G = B$). This quantity (3.6) is a measure similar to coherence (3.5) and is used by Hampton and Doostan (2015a) and Cohen and Migliorati (2017) together with the induced probability measure to obtain convergence results for weighted least-squares regression. Narayan et al. (2017) and Jakeman et al. (2017) choose as probability distribution the so-called weighted pluripotential equilibrium measure (possibly degree-dependent), which asymptotically coincides with $\tilde{f}(\boldsymbol{x}) = c^2 w(\boldsymbol{x})^2 f(\boldsymbol{x})$ when the total degree of the truncated basis $p \to \infty$. However, the modified basis is not orthonormal with respect to this measure, which leads to weaker theoretical recovery results. Theoretical results are available only for the univariate case. In numerical examples, the method performs well for low-dimensional high-degree cases and often very similarly to asymptotic sampling. In high dimensions, it performs worse than input sampling (i.e., MC). It has not been compared to coherence-optimal sampling.
### 3.A.3 Choosing points according to an optimality criterion from a candidate set

The following methods aim to improve the properties of the regression matrix by choosing the “best” design points from a large candidate set. The methods differ in the criterion defining what are the “best” points. Most presented algorithms are greedy or heuristic and are actually only able to find a suboptimal design (local optimum). The choice of the candidate set obviously influences the quality of the resulting design. In the literature, candidate sets were sampled from MC (Diaz et al., 2018), LHS (Fajraoui et al., 2017), coherence-optimal sampling (Diaz et al., 2018; Alemazkoor and Meidani, 2018a), or Christoffel sparse approximation (Shin and Xiu, 2016). In the case of coherence-optimal sampling or Christoffel sparse approximation, the resulting optimized sample inherits the weights. It also often preserves the spread of the candidate set, as can be seen in Figure 3.A.2.

#### 3.A.3.1 D-optimal sampling

D-optimal design of experiments (Kiefer and Wolfowitz, 1959; Dykstra, 1971) aims at maximizing the determinant of the so-called information matrix $\frac{1}{N} \Psi^T \Psi \in \mathbb{R}^{P \times P}$. The D-value is defined as

$$D(\Psi) = \det(\Psi^T \Psi). \qquad (3.20)$$

Sometimes the determinant of the inverse information matrix is minimized (Nguyen and Miller, 1992), or the Pth root is taken for normalization purposes (Diaz et al., 2018). The maximization of this determinant is connected to the minimization of the variance of the PCE coefficient estimate (Nguyen and Miller, 1992; Zein et al., 2013). Note that $D(\Psi) = 0$ if $N < P$.

There exists a large selection of methods for constructing D-optimal experimental designs. For an overview of methods for constructing designs following alphabetic optimality criteria (such as A-, D-, or E-optimality), see (Hadigol and Doostan, 2018, Section 4.5).

Here, we only discuss D-optimal sampling based on rank-revealing QR decomposition (RRQR) (Diaz et al., 2018), since this is the technique used in our benchmark. We decided to use RRQR-based D-optimal sampling because it can be used even in the case $N < P$ when other D-optimal methods fail due to singularity of the information matrix. Note that RRQR is not guaranteed to find a design with maximal D-value but only a local optimum (Diaz et al., 2018, Section 3.4).

Let $\Psi_{cand} \in \mathbb{R}^{M \times P}$ be the regression matrix evaluated at a set of $M$ candidate points. The goal is to select $N \le M$ points from this candidate set with the property that the D-value of the resulting regression matrix $\Psi \in \mathbb{R}^{N \times P}$ is as large as possible. Since in the case of sparse PCE often $N < P$, which leads to $D(\Psi) = 0$, another strategy is necessary. The RRQR decomposition, also known as pivoted QR, aims at permuting the columns of the original matrix in a way that ensures the R-matrix of the associated QR decomposition is as well-behaved<fnref n="13" /> as possible. This is useful for inexpensively determining the numerical rank of a matrix (Hong and Pan, 1992; Gu and Eisenstat, 1996). RRQR has a strong connection to SVD and to the

<footnote n="13">¹³$R = \begin{pmatrix} A_k & B_k \\ 0 & C_k \end{pmatrix}$ where $A_k$ is well-conditioned and $||C_k||_2$ is small</footnote>
selection of submatrices of maximal determinant (Hong and Pan, 1992). Gu and Eisenstat (1996) propose a pivoted QR decomposition where pivots are chosen to maximize the determinant of the resulting quadratic submatrix of R. The exchange of rows is based on a formula relating the determinant of the quadratic submatrix of R before and after the row exchange by a simple factor (Gu and Eisenstat, 1996, Lemma 3.1). This algorithm can be used together with SVD to perform subset selection to construct an initial experimental design from a large set of candidate samples (Seshadri et al., 2017; Diaz et al., 2018). Here, first an SVD of the matrix $\Psi_{cand}^T$ is computed. Then RRQR is applied to the transpose of the matrix consisting of the first N right singular vectors. The resulting permutation matrix is used to determine the points to be chosen from the candidate set.

### 3.A.3.2 Quasi-optimal sampling based on the S-value

Here, the idea is to select samples from a pool of candidate points so that the PCE coefficients obtained using the selected set are as close as possible to the coefficients that would be obtained if the whole set of candidate points was used (Shin and Xiu, 2016a). Under the assumption that the columns of the matrix $\Psi_{cand}$ are mutually orthogonal, the S-value is defined by

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

where $\Psi_i$ denotes the ith column of the regression matrix $\Psi$. Its maximization has the (heuristic) 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). It holds that $S(\Psi) \in [0, 1]$ due to Hadamard's inequality. If $N < P$, $S(\Psi) = 0$. If $N \ge P$, $S(\Psi) = 1$ if and only if the columns of $\Psi$ are mutually orthogonal. There exists an update formula for the S-value when the regression matrix is augmented by one row, which thus avoids the repeated calculation of determinants.

Shin and Xiu (2016a) suggest a greedy algorithm that in every iteration augments the current matrix by an additional row which maximizes the S-value of the resulting matrix among all candidate rows. When the current number of rows in the matrix $\Psi$ is smaller than the number of columns, the procedure can be adapted to avoid $S(\Psi) = 0$. We do not include it in our benchmark because it is not well suited for situations where there are more basis polynomials than design points, which is the case in sparse PCE without experimental design enrichment. However, in a sequential enrichment context (Fajraoui et al., 2017) and for least-squares regression (Shin and Xiu, 2016), this algorithm performs well.

### 3.A.3.3 Near-optimal sampling

The coherence parameter (3.15) gives a bound on the recovery rate, but it is not the only criterion that has been studied with respect to recovery accuracy. Two other matrix properties related to recovery accuracy are mutual coherence and average cross-correlation. Both of them consider the correlation between normalized columns of the regression matrix, i.e., their scalar product. They are scalar measures of how “orthonormal” the columns of a rectangular matrix
$\Psi \in \mathbb{R}^{N \times P}$ with $N < P$ are. The heuristic idea is that columns should point in as different directions as possible, so that the multiplication with sparse coefficient vectors, which results in a linear combination of a subset of the columns, is “as unique as possible”. This facilitates the recovery of the true sparse solution (assuming that it exists).

The mutual coherence is defined by

$$ \mu(\Psi) = \max_{i \neq j} \frac{|\Psi_i^T \Psi_j|}{\|\Psi_i\|_2 \|\Psi_j\|_2}, \qquad (3.22) $$

where $\Psi_i$ denotes the $i$th column of the regression matrix $\Psi \in \mathbb{R}^{N \times P}$. The mutual coherence is the worst-case cross-correlation between any two columns of the matrix. It is zero for orthonormal matrices and positive for $N < P$.

The average (squared) cross-correlation is defined by

$$ \gamma(\Psi) = \frac{1}{P(P-1)} \|1_P - \tilde{\Psi}^T \tilde{\Psi}\|_F^2 = \frac{1}{P(P-1)} \sum_{i \neq j} \frac{|\Psi_i^T \Psi_j|^2}{\|\Psi_i\|_2^2 \|\Psi_j\|_2^2} \qquad (3.23) $$

where $\tilde{\Psi}$ is the column-normalized version of $\Psi$, and $\Psi_i$ denotes the $i$th column of the regression matrix. The norm is the Frobenius-norm, taking the sum of squares of all matrix entries, and the factor $P(P-1)$ is the number of column pairs.

Alemazkoor and Meidani (2018a) suggest simultaneously optimizing mutual coherence and average cross-correlation by using the greedy procedure described in Algorithm 1 below: In each iteration, the current regression matrix is augmented by one row. This row corresponds to that point $x_j$ from the large pool of candidate points which minimizes the (normalized) distance of $(\mu'_j, \gamma'_j) \in \mathbb{R}^2$ to the “utopia point” $(\min(\mu'), \min(\gamma'))$ among all candidate points.

<caption>**Algorithm 1** Near-optimal sampling (Alemazkoor and Meidani, 2018a).</caption>

1: Sample a large number $M$ of candidate points from the coherence-optimal distribution and compute candidate rows arranged in a matrix $\Psi_{cand}$
2: Initialize $\Psi_{opt(1)}$ to be a random row from $\Psi_{cand}$
3: **for** $i = 2 \dots N$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1 \dots M$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \Psi_{\text{temp}} = \text{row-concatenate}(\Psi_{\text{opt}(i-1)}, \Psi_{\text{cand}}^{(j)}) $$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \mu'_{j} = \mu(\Psi_{\text{temp}}) \text{ and } \gamma'_{j} = \gamma(\Psi_{\text{temp}}) $$
7: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
8: &nbsp;&nbsp;&nbsp;&nbsp;$$ \mu' = (\mu'_{1}, \dots, \mu'_{M}) \text{ and } \gamma' = (\gamma'_{1}, \dots, \gamma'_{M}) $$
9: &nbsp;&nbsp;&nbsp;&nbsp;$$ j^* = \arg\min_j \left( \frac{\mu'_{j} - \min(\mu')}{\max(\mu') - \min(\mu')} \right)^2 + \left( \frac{\gamma'_{j} - \min(\gamma')}{\max(\gamma') - \min(\gamma')} \right)^2 $$
10: $$ \Psi_{\text{opt}(i)} = \text{row-concatenate}(\Psi_{\text{opt}(i-1)}, \Psi_{\text{cand}}^{(j^*)}) $$
11: **end for**

The algorithm is called near-optimal because it is a greedy algorithm, finding only a local optimum, and because optimized mutual coherence and average cross-correlation are only hinting at, but not guaranteeing, good recovery accuracy (Alemazkoor and Meidani, 2018a). Its computational complexity is $O(NMP^2)$, where $N$ is the size of the final experimental design, $M$ is the
number of candidate samples (chosen to be, e.g., proportional to *P* (Diaz et al., 2018)), and *P* is the number of regressors. This makes the algorithm prohibitively expensive in the case of large bases (*P* in the order of thousands), which is why we do not use it for some of the benchmark examples.

### 3.A.4 Illustration of sampling schemes

In Figures 3.A.1 and 3.A.2, we show illustrations of experimental designs in *d* = 2 dimensions with *N* = 100 and *p* = 12 for selected sampling techniques. The candidate set has a size of *M* = 1000. Figure 3.A.1 presents experimental designs for uniform input in the interval [-1, 1], while Figure 3.A.2 presents experimental designs for standard Gaussian input.

Note that in the standard Gaussian case, the asymptotic distribution, the coherence-optimal distribution, and the matrix-optimal distributions based on a coherence-optimal candidate set all have a very large spread that grows with the total degree of the basis. For degree *p* = 12, some points are seven standard deviations away from the mean. Engineering models are typically calibrated only for a certain region of the input domain corresponding to nonnegligible probability, and they may be less accurate (or even fail) outside of this region.
<figure><img src="image_19.png" /><figcaption>**Figure 3.A.1:** Visualization of experimental designs constructed for uniform input in [-1, 1]² for degree p = 12. Red filled points denote the chosen experimental design, while blue circles denote the candidate set. Size of the ED: N = 100, size of the candidate set: M = 1000.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 3.A.2:** Visualization of the experimental design constructed for standard Gaussian input in $d = 2$ dimensions for degree $p = 12$. The gray surface plot illustrates the Gaussian probability density function. Red filled points denote the chosen experimental design, while blue circles denote the candidate set. Size of the ED: $N = 100$; size of the candidate set: $M = 1000$. The support of the asymptotic and the coherence-optimal distribution is the ball of radius $r = \sqrt{2\sqrt{2p+2}} \approx 7.2$. Note that engineering models may be less accurate in regions where the input distribution has negligible mass.</figcaption></figure>
## 3.B Details on sparse regression solvers

In this appendix, we describe the sparse solvers used in our benchmark in more detail: LARS, OMP, subspace pursuit, SPGL1 and FastLaplace (BCS). In addition, we present an overview of greedy stepwise regression solvers for sparse PCE.

There exist various formulations for the sparse regression problem. The typical form minimizes the $ℓ^2$-norm of the empirical error under an additional constraint that is designed to enforce sparsity.

Sparsity is measured by the number of nonzero entries in a vector, formally denoted by $||c||_0 = \sum_i 1_{\{c_i \neq 0\}}$ (even though this expression is not a norm). This results in the sparse regression problem

$$ \hat{c} = \arg \min_{c} ||\Psi c - y||_2^2 + \lambda ||c||_0 \quad (3.24) $$

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

The convex relaxation of this problem is $ℓ^1$-minimization, where $||c||_0$ is replaced by $||c||_1 = \sum_i |c_i|$. There are several equivalent formulations of the relaxed problem, namely

$$ \hat{c} = \arg \min_{c} ||\Psi c - y||_2^2 + \lambda ||c||_1 \quad (3.25) $$

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

$$ \hat{c} = \arg \min_{c} ||\Psi c - y||_2 \quad \text{s.t.} \quad ||c||_1 \le \tau \quad (3.27) $$

called Lagrangian formulation, basis pursuit denoising (BPDN), and least absolute shrinkage and selection operator (LASSO), respectively. It has been shown that under certain conditions, the solutions to $ℓ^0$-minimization and $ℓ^1$-minimization coincide (Bruckstein et al., 2009). However contrary to (3.24), formulations (3.25)–(3.27) are convex problems and allow a numerical solution with considerably smaller cost. The three formulations (3.25), (3.26), and (3.27) are equivalent in the sense that if $\hat{c}$ is solution to one of the formulations, there exists a value of constraint parameter $\sigma, \tau$, or $\lambda$ so that $\hat{c}$ is also a solution to the other formulations. However, the relationship between the parameters $\sigma, \tau, \lambda$ that makes the problems equivalent depends on $\Psi$ and $y$ and is not known in advance (van den Berg and Friedlander, 2008).

There exist other sparsity-enforcing formulations, such as $ℓ^p$-norms (Bruckstein et al., 2009), $ℓ^1$ – $ℓ^2$-minimization (Yin et al., 2015), or elastic net (Tarakanov and Elsheikh, 2019). One example that we will describe is Bayesian compressive sensing, where a sparsity-enforcing prior is used for the coefficients of the PCE, resulting in a formulation that is equivalent to a sparse regression problem with a different kind of sparsity constraint, e.g., one related to the Student-t distribution (Tipping, 2001).

In the following descriptions of the algorithms, $A \subset \mathbb{N}^d$ with various sub- or superscripts denotes a set of multi-indices, which by definition of PCE can be identified with a set of basis polynomials. With the notation of Section 3.2.1, $y \in \mathbb{R}^N$ denotes the vector of model responses, $\Psi \in \mathbb{R}^{N \times P}$ denotes the regression matrix of basis polynomials evaluated at the $N$ experimental design points, and $c \in \mathbb{R}^P$ denotes the coefficients of a PCE. The residual is defined by $r = y - \Psi c$, so that the norm of the residual is the empirical error.
The polynomials used for building the PCE are sometimes also called basis functions, regressors, or predictors. A sparse PCE is a PCE for which only some of the basis functions have nonzero coefficients: these basis functions are called active. We assume the regressors are normalized, so that the correlation between two regressors is equivalent to their inner product (and to the cosine of the angle between them).

### 3.B.1 Orthogonal matching pursuit (OMP)

Orthogonal matching pursuit (OMP), also called forward stepwise regression, is a classical greedy technique for finding approximate solutions to the $l^0$-minimization problem (3.24) (Pati et al., 1993; Bruckstein et al., 2009). Despite its being a heuristic method, under certain assumptions there are theoretical guarantees for the solutions returned by OMP (Tropp and Gilbert, 2007; Bruckstein et al., 2009). OMP is an iterative algorithm that starts out with an empty model and adds the regressors one by one to the set of active regressors. In each iteration, OMP selects the regressor that is most correlated with the current residual, adds it to the set of active regressors, and then updates the coefficients of all active regressors to make sure the new residual is orthogonal to all of them and has smallest possible norm. The updating of the coefficients can be done through an update formula (Berchier, 2015) or by computing the least-squares solution to the system of equations involving only the active regressors (Marelli and Sudret, 2019).

The technique is presented in Algorithm 2. The iterations are continued until $\min\{N, P\}$ basis functions are in the active set (then either all polynomials are selected, or there are not enough points in the experimental design to use least-squares anymore).

**Algorithm 2** Orthogonal matching pursuit (OMP) (Pati et al., 1993; Tropp and Gilbert, 2007; Marelli and Sudret, 2019).

1: Given a set of candidate basis functions $A_{cand}$

2: Initialize all coefficients to zero: $c_0 = 0$. Set $A_0 = \emptyset$

3: Set the residual vector $r := y$

4: **for** $i = 1, \dots, m$ **do**
$$m \le \min\{N, P\}$$

▷ OMP can be stopped early when the error did not decrease anymore for a while

5: Find $\alpha^* \in A_{cand} \setminus A_{i-1}$ with maximal correlation with the residual by solving

$$\alpha^* = \arg \max_{\alpha \in A_{cand} \setminus A_{i-1}} |\mathbf{r}^T \psi_{\alpha}|$$

▷ The entries of vector $\psi_{\alpha}$ are evaluations of the basis function $\psi_{\alpha}$ at the ED

6: $A_i = A_{i-1} \cup \{\alpha\}$
$$\text{▷ Current set of active predictors}$$

7: Compute the coefficients $c_i$ by least-squares using only the active indices $A_i$

▷ This can be done in $O(iN)$ when maintaining a QR factorization (Tropp and Gilbert, 2007)

8: Update the residual $\mathbf{r} = \mathbf{y} - \mathbf{\Psi}_{A_i} c_i$

9: end for

OMP does not, per se, return a sparse solution. If a desired level of sparsity $K$ is known a priori, the algorithm can be stopped after $K$ iterations. Another possibility is to stop the algorithm as
soon as the residual norm is smaller than some error threshold (Bruckstein et al., 2009; Doostan and Owhadi, 2011; Jakeman et al., 2015), where the best error threshold is determined through cross-validation. A third possibility is to determine the best number of active basis functions through a model selection criterion, e.g. the LOO error (Marelli and Sudret, 2019). Since the coefficients are computed by OLS on the active basis, the LOO can be computed cheaply (Chapelle et al., 2002; Blatman and Sudret, 2011). Typically, for an increasing sequence of basis functions the LOO error first decreases (reduction of underfitting), then increases (overfitting). This can be utilized to terminate the algorithm early once the LOO error starts rising (early-stop criterion) (Marelli and Sudret, 2019).

The computational complexity of OMP is $O(mNP)$ (Tropp and Gilbert, 2007; Dai and Milenkovic, 2009), where $m \le \min\{N, P\}$ is the number of iterations. The computation of the correlations of the current residual with all regressors is $O(NP)$ and has to be performed $m$ times. The computation of the least-squares solution in step $i$ can be done in $O(iN)$, e.g., by maintaining a QR factorization of the information matrix (Tropp and Gilbert, 2007), or by using Schur's complement to update the information matrix inverse whenever a new regressor is added.

From the authors' experience, OMP often suffers from overfitting and can produce an unreliable LOO error estimate, which can be detrimental in basis-adaptive settings (see also (Lüthen et al., 2022a)).

OMP is available in many software packages, among them UQLab (Marelli and Sudret, 2014).

### 3.B.2 Least angle regression (LARS)

Least-angle regression (LARS, sometimes also abbreviated LAR) is a greedy technique that finds an approximate solution to the $l^1$-minimization problem (Efron et al., 2004). It is similar to OMP in that the algorithm starts out with an empty model and adds regressors one by one based on their correlation with the residual. However, unlike OMP, which updates the coefficients using least-squares (making the residual orthogonal to all active regressors in each step), LARS updates the coefficients in such a way that all active regressors have equal correlation with the residual. LARS can be interpreted as producing a path of solutions to (3.27), corresponding to increasing $\tau$. The coefficients are increased in the equiangular direction until a nonactive regressor has as much correlation with the residual as all the active regressors. This regressor is then added to the set of active regressors and the new equiangular direction is computed. The optimal stepsize between the addition of subsequent regressors can be computed analytically (Efron et al., 2004). This algorithm solves (3.27) approximately. A slightly modified version of LARS, called LARS-LASSO, removes regressors whenever the sign of their coefficient changes, and it has been proven to solve (3.27) (or its noiseless counterpart) exactly under certain conditions (Efron et al., 2004; Bruckstein et al., 2009).

The LARS technique is presented in Algorithm 3. It returns a sequence $A_1 \subset A_2 \subset \dots \subset A_m$ of sets containing indices of active basis functions, with $m = \min\{N, P\}$. Just like OMP, LARS can be stopped when a predefined sparsity $K$ is reached or when the norm of the residual $\|r\|_2$ falls below a predefined error threshold.
**Algorithm 3** Least angle regression (LARS) (Efron et al., 2004; Blatman and Sudret, 2011).

1: Given a set of candidate basis functions $\mathcal{A}_{\text{cand}}$
2: Initialize all coefficients to zero: $c_0 = 0$. Set $\mathcal{A}_0 = \emptyset$
3: Set the residual vector $\mathbf{r} := \mathbf{y}$
4: **for** $i = 1, \dots, m$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;Find $\alpha \in \mathcal{A}_{\text{cand}} \setminus \mathcal{A}_{i-1}$ with maximal correlation with the residual
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ For $i > 1$, $\alpha_i$ is the element "responsible for" $\gamma_i$
6: &nbsp;&nbsp;&nbsp;&nbsp;7: Compute $\mathbf{c}_i$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Current set of active predictors
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Equiangular direction for all $\alpha \in \mathcal{A}_i$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ $\mathbf{c}_1$ is equal to the first selected predictor
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ (Efron et al., 2004, Eq. 2.6)
8: Compute $\gamma_i$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ Optimal stepsize: using this, there is a new regressor that is as much correlated with $\mathbf{r}$ as all regressors in $\mathcal{A}_i$ are
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;▷ (Efron et al., 2004, Eq. 2.13)
9: Compute the new coefficients $\mathbf{c}_i = \mathbf{c}_{i-1} + \gamma_i \mathbf{c}_i$ &nbsp;&nbsp;&nbsp;&nbsp;▷ Move the coefficients jointly into the direction of the least-squares solution until one of the other predictors in $\mathcal{A}_{\text{cand}} \setminus \mathcal{A}_i$ has as much correlation with the residual as the predictors in $\mathcal{A}_i$ (ensured by choice of $\gamma_i$ and $\mathbf{c}_i$)
10: Update the residual $\mathbf{r} = \mathbf{y} - \mathbf{\Psi}_i \mathbf{c}_i$
11: **end for**

A modified version of LARS, called *hybrid LARS*, uses the equicorrelated approach to select the predictors, but computes the coefficients of the metamodel by least-squares (Efron et al., 2004; Blatman and Sudret, 2011). Once the LARS algorithm has finished and returned the sequence of basis sets $\mathcal{A}_1, \dots, \mathcal{A}_m$, the corresponding coefficients are recomputed by least squares, $\mathbf{c}_i = \mathbf{c}_i^{\text{LSQ}}$, which ensures minimal empirical error for every metamodel $(\mathcal{A}_i, \mathbf{c}_i)$. Hybrid LARS facilitates another way to choose the best sparsity level: as for OMP, a model selection criterion (e.g. LOO) for each metamodel can be evaluated, and the best one is chosen. This procedure is detailed in Algorithm 4. Cheap OLS-based computation of LOO (Chapelle et al., 2002; Blatman and Sudret, 2011) and the early-stop criterion (Marelli and Sudret, 2019) can be applied as well.

**Algorithm 4** Hybrid LARS with LOO-CV (Blatman and Sudret, 2011; Marelli and Sudret, 2019)

1: Initialization as in LARS (Algorithm 3)
2: **for** $i = 1, \dots, m$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Run one step of LARS and obtain $(\mathcal{A}_i, \mathbf{c}_i)$ &nbsp;&nbsp;&nbsp;&nbsp;▷ Algorithm 3
4: &nbsp;&nbsp;&nbsp;&nbsp;Recompute the coefficient vector using least-squares on the selected basis $\mathcal{A}_i$ only, obtaining $\mathbf{c}_i^{\text{OLS}}$ (the coefficients corresponding to $\mathcal{A} \setminus \mathcal{A}_i$ are set to zero) &nbsp;&nbsp;&nbsp;&nbsp;▷ Hybrid LARS
5: &nbsp;&nbsp;&nbsp;&nbsp;Compute the LOO error $\epsilon_{\text{LOO}}(i)$ for $\mathbf{c}_i^{\text{OLS}}$ &nbsp;&nbsp;&nbsp;&nbsp;▷ OLS-based LOO computation (Chapelle et al., 2002; Blatman and Sudret, 2011)
6: **end for** &nbsp;&nbsp;&nbsp;&nbsp;▷ early stopping possible by monitoring the LOO error (Marelli and Sudret, 2019)
7: Return the metamodel $(\mathcal{A}_{i^*}, \mathbf{c}_{i^*})$ with $i^* = \arg\min_i \epsilon_{\text{LOO}}(i)$
As for OMP, the computational complexity of LARS (in the case $N < P$) is $O(mNP)$, where $m \le \min\{N, P\}$ is the number of iterations. This is due to matrix-vector multiplication and matrix inversion which have to be performed in every iteration. The latter can be computed in $O(mN)$, when using techniques such as Schur's complement to update the information matrix inverse whenever a new regressor is added.

LARS is available in many software packages, e.g., as MATLAB implementation in UQLab (Marelli and Sudret, 2014).

### 3.B.3 Subspace pursuit (SP)

Another formulation of the $\ell^0$-minimization problem is

$$ \min_{c \in \mathbb{R}^P} \|\Psi c - u\|_2 \quad \text{s.t.} \quad \|c\|_0 = K \qquad (3.28) $$

which is equivalent to (3.24) for a certain choice of $\lambda$.

Subspace pursuit (SP) seeks to identify a solution to (3.28) by iteratively and greedily enlarging and shrinking the set of active basis functions (Dai and Milenkovic, 2009). As with LARS and OMP, regressors are added to the set of active basis functions according to their correlation with the residual. However, the regressors are not added one by one, but batchwise. More precisely, SP maintains at all times an active basis of size $K$, where $K$ denotes the desired sparsity. In each iteration, it adds $K$ regressors at once and computes the coefficients of the active regressors by OLS. Then, it removes the $K$ regressors with the smallest-in-magnitude coefficients. This is continued until convergence. Under certain assumptions, there are theoretical guarantees for the solution that SP returns (Dai and Milenkovic, 2009). To make the augmentation of the basis and the OLS regression feasible, it must hold that $2K \le \min\{N, P\}$.

The technique is described in Algorithm 5 for a fixed value of sparsity $K$. The residual of a vector and a regression matrix is defined as

$$ \text{residual}(y, \Psi) = y - \Psi \Psi^{\dagger} y \qquad (3.29) $$

where $\Psi^{\dagger}$ denotes the pseudoinverse of $\Psi$ and $\Psi^{\dagger} y = c$ is the least-squares solution to $\Psi c \approx y$ (the case of an overdetermined system). The algorithm returns a set $A$ containing $K$ multi-indices.

For arbitrary sparse vectors, the computational complexity is $O(N(P+K^2)K)$ (Dai and Milenkovic, 2009). For very sparse vectors with $K^2 \in O(P)$, the complexity thus becomes $O(NPK)$, comparable to the runtime of OMP. The number of iterations that the SP algorithm performs can be shown to be $O(K)$ in general and even $O(\log K)$ in certain cases (Dai and Milenkovic, 2009).

When the optimal sparsity level $K$ is unknown, it can be determined e.g. by cross-validation: Diaz et al. (2018) suggest running Algorithm 5 for a range of $N_K = 10$ different values for $K$ and choosing the one with the smallest 4-fold cross-validation error. In this paper, we propose to use leave-one-out cross-validation instead of 4-fold cross-validation, resulting in the SP variant $SP_{LOO}$.
**Algorithm 5** Subspace pursuit (SP) (Dai and Milenkovic, 2009).

1: Given desired sparsity $K \le \min\{\frac{N}{2}, \frac{P}{2}\}$

2: Given the experimental design and the candidate basis $\mathcal{A}_{\text{cand}}$, compute the associated regression matrix $\mathbf{\Psi}$ and the right-hand-side $\mathbf{y}$

3: $\mathcal{A}^0 = \{K \text{ indices corresponding to the largest magnitude entries in } \mathbf{\Psi}^T \mathbf{y}\}$

▷ Scalar product of columns of $\mathbf{\Psi}$ with $\mathbf{y}$

4: $\mathbf{y}_{\text{res}}^0 = \text{residual}(\mathbf{y}, \mathbf{\Psi}_{\mathcal{A}^0})$

▷ Residual of least-squares solution based on active basis

5: **for** $l = 1, 2, \dots$ **do**

6: $\mathcal{S}^l = \mathcal{A}^{l-1} \cup \{K \text{ indices corresponding to the largest magnitude entries in } \mathbf{\Psi}^T \mathbf{y}_{\text{res}}^{l-1}\}$

▷ Augment by indices of full basis that correlate best with the residual

7: $\mathbf{c} = \mathbf{\Psi}_{\mathcal{S}^l}^\dagger \mathbf{y}$

▷ Least-squares solution based on set $\mathcal{S}^l$ of size $2K$

8: $\mathcal{A}^l = \{K \text{ indices corresponding to the largest magnitude entries in } \mathbf{c}\}$

9: $\mathbf{y}_{\text{res}}^l = \text{residual}(\mathbf{y}, \mathbf{\Psi}_{\mathcal{A}^l})$

▷ Residual of least-squares solution based on $\mathcal{A}^l$

10: **if** $\|\mathbf{y}_{\text{res}}^l\|_2 \ge \|\mathbf{y}_{\text{res}}^{l-1}\|_2$ **then**

▷ if new $K$-sparse approx. is worse than the previous one

11: STOP iteration and **return** $\mathcal{A}^{l-1}$.

12: **end if**

13: **end for**

**Remark:** In line 10, the original publication (Dai and Milenkovic, 2009) uses ">" instead of ">=", but we also want to stop when the set has converged.

A related algorithm is CoSAMP (Needell and Tropp, 2009), which differs from SP mainly in the number of regressors added in each iteration.

Subspace pursuit is available as MATLAB implementation in the software package `DOPT_PCE` (Diaz et al., 2018; Diaz, 2018).

### 3.B.4 SPGL1

$l^1$-minimization is a convex problem, since both the objective function and the constraint are convex functions. Therefore, convex optimization methods can be used to find a solution. In this section, we describe the algorithm SPGL1 (van den Berg and Friedlander, 2008).

For a given value of $\tau$, formulation (3.27) (LASSO) can be solved by spectral projected gradient (SPG) descent (Birgin et al., 2000; van den Berg and Friedlander, 2008).<fnref n="14" /> However, for real-world problems, we often do not know a priori an appropriate value for $\tau$. On the other hand, a sensible range of values for $\sigma$ in formulation (3.26) (BPDN) can typically be estimated based on the noise level in the data and the expected model fit. In the case of PCE metamodelling, $\sigma$ can be related to an estimate of the relative MSE through $\text{RelMSE} = \frac{\sigma^2}{N\text{Var}[\mathbf{y}]}$, whose values are for engineering models typically between $10^{-10}$ and $10^0 = 1$.

The main idea of the solver SPGL1 is to solve BPDN through a detour over LASSO. Let $c_\tau$ be

<footnote n="14">¹⁴SPG is a gradient-based optimization algorithm with several enhancements (Barzilai-Borwein spectral step length and the Grippo-Lampariello-Lucidi scheme of nonmonotone line search) and projection onto the feasible set $\Omega_\tau = \{c \in \mathbb{R}^P : \|c\|_1 \le \tau\}$</footnote>
the solution to LASSO for a given $\tau$. Define a function $\phi: \mathbb{R}_+ \to \mathbb{R}_+$ by

$$ \phi(\tau) := \| \Psi c_{\tau} - y \|_2. \qquad (3.30) $$

Then the solution to BPDN with $\sigma := \phi(\tau)$ is $c_{\tau}$. In other words, $\phi$ is the functional relationship between $\sigma$ and $\tau$ that makes the two formulations BPDN and LASSO equivalent for given $\Psi$ and $y$. $\phi$ is the Pareto front of LASSO and BPDN and shows the trade-off between the minimal achievable $\ell^1$-norm of the coefficients and the minimal $\ell^2$-norm of the corresponding residual. The Pareto front is convex, nonincreasing and differentiable with an analytically computable derivative (van den Berg and Friedlander, 2008).

To find a solution to BPDN with a given $\sigma$, LASSO is solved with SPG several times for a sequence of $\tau$ until one is found with $\phi(\tau) = \sigma$. The sequence of $\tau$ is created by performing Newton's root finding algorithm on the function $f(\tau) = \sigma - \phi(\tau)$.

Each SPG iteration has a computational complexity of $O(NP + P \log P)$ (from matrix-vector multiplication and $\ell^1$-projection). Multiplying this with the number of SPG steps and the number of Newton steps yields the computational complexity of SPGL1.

This algorithm is available as MATLAB package SPGL1 (van den Berg and Friedlander, 2008; Van den Berg and Friedlander, 2015).

In our numerical benchmarks computing sparse PCE for compressible models, SPGL1 was among the slowest solvers and often returned rather dense solutions.

### 3.B.5 Sparse Bayesian learning

Methods from the class of Bayesian compressive sensing (BCS), also known as sparse Bayesian learning (SBL), embed the regression problem in a probabilistic framework (Tipping, 2001; Ji et al., 2008; Babacan et al., 2010; Sargsyan et al., 2014; Tsilifis et al., 2020). The goal is to compute, for a given model response vector $y$ and a regression matrix $\Psi$, the coefficient vector $c^{\text{MAP}}$ which maximizes the posterior distribution $p(c|y)$. Another quantity of interest could be the most probable value $y^*$ at a new point $x^*$ maximizing $p(y^*|y)$.

In BCS, it is assumed that the “measurements” $y$ are generated by adding zero-mean, finite-variance noise to the evaluations of the true model. This noise is often assumed to be Gaussian white noise with standard deviation $\sigma$, which, for a given input $x$, results in a Gaussian distribution for its output $y$ with mean $\Psi c$ and covariance matrix $\sigma^2\mathbf{1}$, i.e., $y|c, x, \sigma \sim N(\Psi c, \sigma^2\mathbf{1})$. Note that in the case of PCE, this is generally not a valid assumption: when an important term is missing from the PCE model, the discrepancy between measurements and PCE model evaluations can be highly correlated, heteroscedastic, and non-Gaussian, and have nonzero mean. However, even though the assumptions might not be fulfilled, this framework can still be useful for finding sparse solutions.

The class of BCS algorithms comprises several methods that differ in the assumptions on the distributions of the various hyperparameters and in the (usually iterative, approximate) techniques for computing the posterior quantities.
In Figure 3.B.1a we present the general setup of sparse Bayesian learning. The measurements $\mathbf{y}$ are assumed to follow a Gaussian distribution as described above. The noise variance $\sigma^2$ is assumed to be a random variable whose distribution has to be specified (e.g. fixed, uniform, or inverse-Gamma). The coefficients $c_i$ are assumed to be random variables as well, drawn from a normal distribution with mean zero and variance $\gamma_i$, i.e., each weight has its own variance. $\gamma$ is a so-called hyperparameter, parametrizing the distribution of a parameter.

<figure><img src="image_21.png" /><figcaption>**Figure 3.B.1:** Illustration of the general setup of BCS (a) and the hierarchical generalization of Babacan et al. (2010) (b). (a): The likelihood and the prior on the coefficients are usually Gaussian, but the choice of $p(\sigma^2)$ and $p(\gamma)$ differs between publications, as well as the resulting solution algorithm. (b): Babacan et al. (2010) makes a specific choice for $p(\gamma)$ and includes an additional layer of hyperparameters. Shaded variables are held fixed.</figcaption></figure>

So far, the described setup with fixed $\sigma^2$ and $\gamma_i$ would yield (weighted) ridge regression. The sparsity comes into play through an assumption on the distribution of the hyperparameter $\gamma$. For specific choices of $p(\gamma_i)$, it can be shown that the resulting effective prior on the coefficients $p(\mathbf{c}) = \int p(\mathbf{c}|\gamma)p(\gamma)d\gamma$ is a sparsity-encouraging distribution, i.e., one that has a sharp peak at zero, encouraging zero values, while at the same time a heavy tail, allowing for large coefficient values as well. Examples are the Laplace distribution and the generalized Student-t distribution (Wipf et al., 2004; Babacan et al., 2010; Figueiredo and Nowak, 2001).

Such sparsity-encouraging distributions are often intractable to use, because they do not allow for analytical computation of the desired values (such as the most likely coefficients given the data, or the prediction of the measurement value at a new point). However, feasible algorithms can be developed based on a suitable approximation step. Various frameworks for sparse Bayesian learning have been proposed whose setup follows the general structure of Figure 3.B.1a, but which differ in the choice of priors for the hyperparameters and employ different solution algorithms for the MAP estimate of these hyperparameters (Tipping, 2001; Faul and Tipping, 2002; Figueiredo, 2003; Tipping and Faul, 2003; Wipf and Rao, 2004; Ji et al., 2008; Seeger and Nickisch, 2008; Babacan et al., 2010; Sargsyan et al., 2014; Tsilifis et al., 2020; Bhattacharyya, 2020). The MAP estimate of the hyperparameters is inserted into the distribution for $\mathbf{c}$. Because $p(\mathbf{c}|\gamma)$ and $p(\mathbf{y}|\mathbf{c}, \sigma^2)$ are normal distributions, all subsequent computations can be carried out analytically (Tipping, 2001; Wipf et al., 2004; Wipf and Rao, 2004). The sparsity of $\mathbf{c}$ is
enforced because by the choice of $p(\gamma)$ and the other distributions, many of the components $\gamma_i^{\text{MAP}}$ of $\gamma^{\text{MAP}}$ will actually be zero, forcing the corresponding $c_i$ to be zero as well.

BCS in the implementation of Babacan et al. (2010) was suggested for sparse PCE by Sargsyan et al. (2014). This approach employs an additional layer of hyperparameters as displayed in Figure 3.B.1b. The prior on the coefficient variances is an exponential distribution $p(\gamma_i|\lambda) = \text{Exp}(\gamma_i|\frac{\lambda}{2})$ with shared hyperparameter $\lambda$. The hyperparameter $\lambda$ follows a Gamma distribution $p(\lambda|\nu) = \Gamma(\lambda|\frac{\nu}{2}, \frac{\nu}{2})$ with hyperparameter $\nu$. $\nu \to 0$ implies $p(\lambda) \propto \frac{1}{|\lambda|}$ (improper prior) and $\nu \to \infty$ implies the certain value $\lambda = 1$. In practice, Babacan et al. (2010) find that $\nu = 0$ gives the best results. The prior on $\beta = \sigma^{-2}$ is a Gamma distribution $p(\beta) = \Gamma(\beta|a, b)$ with hyperparameters $a, b$. In practice, the algorithm does not estimate $\beta$ well, which is, however, crucial; therefore, it is set to a fixed value (e.g. $\beta^{-1} = 0.01 \|y\|_2^2$ in Babacan et al. (2010); in our benchmark, we use cross-validation to determine the best value for this parameter, similarly to the strategy for SPGL1). The objective function is the logarithm of the joint distribution $\mathcal{L}(\gamma, \lambda, \beta) = \log p(y, \gamma, \lambda, \beta)$, which is an analytical expression. To maximize it, Babacan et al. (2010) adapt the fast approximate algorithm of Tipping and Faul (2003); Faul and Tipping (2002) to their generalized hierarchical setting. Here, the derivatives of the objective function with respect to the hyperparameters $\lambda, \beta$, and $\gamma_i$, $i = 1, \dots, P$ are computed. This results in an iterative scheme where these parameters are optimized one at a time while the other ones are held fixed. The algorithm is explained in detail in Babacan et al. (2010, Algorithm 1) and has been implemented in MATLAB under the name FastLaplace (Babacan, 2011).

### 3.B.6 Greedy stepwise regression solvers

Many of the sparse regression solvers that have been proposed for computing sparse PCE belong to the class of greedy stepwise regression. Here, starting from an empty model, the regressors are added one by one according to a selection criterion (forward selection). Some methods also include a backward elimination step. Then the coefficients of the selected regressors are computed. The procedure is iterated until a stopping criterion is reached. Alternatively, several models are built and one is selected in the end using a model selection criterion. We summarize some greedy stepwise regression techniques proposed for sparse PCE, together with their choices for selection criterion, coefficient computation, and stopping criterion, in Table 3.B.1, including the well-known methods OMP and LARS.

New greedy methods in the fashion of Table 3.B.1 can easily be derived by pairing other methods for the regressor selection, the coefficient computation method, and the model selection criterion. Note that except for LARS and OMP, these greedy methods are heuristic (no theoretical guarantee of convergence) and often depend on a number of tuning parameters.
<table><caption>**Table 3.B.1:** A selection of greedy stepwise regression algorithms proposed for sparse PCE. F = forward selection, FB = forward selection and backward elimination.</caption><thead><tr><th>Ref.</th><th>F/B</th><th>Regressor selection</th><th>Computation of the coefficients</th><th>Model selection / stopping criterion</th><th>Comment</th></tr></thead><tbody><tr><td>OMP (Tropp and Gilbert, 2007)</td><td>F</td><td>correlation with residual</td><td>OLS</td><td>several choices: given number of regressors (Tropp and Gilbert, 2007); (modified) LOO (Marelli and Sudret, 2019); norm of residual (Baptista et al., 2019); threshold on moving average of LOO (Baptista et al., 2019)</td><td>theoretical guarantees exist (Tropp and Gilbert, 2007)</td></tr><tr><td>LARS (Efron et al., 2004; Blatman and Sudret, 2011)</td><td>F</td><td>correlation with residual</td><td>least angle strategy</td><td>several choices: coefficient of determination (Efron et al., 2004); (modified) LOO (Blatman and Sudret, 2011)</td><td>theoretical guarantees exist (Efron et al., 2004; Bruck-stein et al., 2009)</td></tr><tr><td>Blatman and Sudret (2008)</td><td>FB</td><td>coefficient of determination</td><td>OLS</td><td>coefficient of determination</td><td>degree-adaptive</td></tr><tr><td>Blatman and Sudret (2010)</td><td>FB</td><td>coefficient of determination</td><td>OLS</td><td>LOO</td><td>degree-adaptive; with ED enrichment</td></tr><tr><td>Hu and Youn (2011)</td><td>F</td><td>testing bivariate interaction</td><td>stepwise moving least-squares</td><td>coefficient of determination</td><td>Only for interaction order ≤ 2. ED-adaptive. Degree of interaction terms determined in inner loop.</td></tr><tr><td>Abraham et al. (2017)</td><td>FB</td><td>F: one-predictor regression criterion B: confidence intervals of coefficients</td><td>OLS</td><td>given number of regressors or iterations</td><td>-</td></tr><tr><td>Shao et al. (2017)</td><td>F</td><td>partial correlation coefficient</td><td>variant of BCS (fixed priors)</td><td>KIC</td><td>prior depends on interaction order and degree of the regressor; degree-adaptive</td></tr><tr><td>Cheng and La (2018a)</td><td>FB</td><td>variance contribution of coefficient</td><td>support vector regression (SVR)</td><td>coefficient of determination</td><td>-</td></tr><tr><td>AFBS (Zhao et al., 2019)</td><td>FB</td><td>correlation with residual</td><td>OLS</td><td>threshold for residual norm</td><td>Main difference with OMP: backward selection. Employs several elimination checks.</td></tr><tr><td>HSPLSR-PCE (Zhao and Bu, 2019)</td><td>FB</td><td>F: partial least-squares (PLS) (regressors partitioned into blocks)</td><td>linear PLS on the regression matrix</td><td>modified (pseudo) LOO</td><td>Note that the input is assumed to be Gaussian, but the (nonlinear) regressors do not follow a Gaussian distribution.</td></tr><tr><td>Zhou et al. (2019a)</td><td>F</td><td>(partial) distance correlation</td><td>PLS</td><td>coefficient of determination</td><td>adaptive in degree and interaction order</td></tr><tr><td>Zhou et al. (2019c)</td><td>FB</td><td>correlation with residual</td><td>BCS</td><td>KIC</td><td>nesting iterative BCS procedure with forward-backward-selection scheme</td></tr></tbody></table>
3.C Additional results

In this appendix, we display additional results that complement the results shown in Sec-
tions 3.3.4–3.3.6. For a detailed description of the setup, we refer the reader to Section 3.3.

3.C.1 Comparison of sparse solvers

In Figure 3.C.1 we display the boxplots of relative MSE for the seven additional models presented in Table 3.1.

3.C.2 Comparison of sampling schemes together with solvers

In Section 3.3.5, Figure 3.4 we showed aggregated results for the benchmark of solvers and sampling schemes. To give a more tangible impression of the data, in Figures 3.C.2–3.C.5 we display the boxplots of relative MSE against ED size for the four models Ishigami, borehole, two-dimensional diffusion, and 100D function. We show all combinations of solvers and sampling schemes, resulting in 16–20 combinations. Solvers are denoted by different colors. Sampling schemes are shown in varying shades and line styles. We also show the same results sliced at small and large ED sizes to compare the performance between solvers.

3.C.3 Comparison of sampling schemes together with solvers, using a smaller
candidate basis

Due to space limitations, in Section 3.3.6 (Figures 3.5 and 3.6) we only showed results for two of
the five solvers (OMP and SP<sub>LOO</sub>). In Figure 3.C.6, we show boxplots of relative MSE against
ED size for the three remaining solvers LARS, SP, and BCS.

**3.D Benchmark studies**

An overview of articles and benchmark studies comparing sparse PCE methods, including their
main results, is given in Table 3.D.1.
<figure><img src="image_22.png" /><figcaption>(a) Undamped oscillator ($d = 6$)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Damped oscillator ($d = 8$)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) Wingweight function ($d = 10$)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) Truss model ($d = 10$)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Morris function ($d = 20$)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(f) Structural frame model ($d = 21$)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(g) 1-dim diffusion model ($d = 62$)</figcaption></figure>

**Figure 3.C.1:** Results for seven additional models (see Table 3.1 in Section 3.3.4 for more details), complementing the results in Figure 3.2. Boxplots of relative MSE against experimental design for six sparse solvers and LHS design. Thirty replications. Note that the damped oscillator and the Morris function are very challenging for PCE: no solver achieves a relative MSE significantly smaller than 0.1, even when large EDs are used.
<figure><img src="image_29.png" /><figcaption>(a) LARS</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) OMP</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(c) SP</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(d) SP<sub>LOO</sub></figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(e) BCS</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(f) Small ED (70 points)</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(g) Large ED (150 points)</figcaption></figure>

**Figure 3.C.2:** Boxplots of relative MSE from the benchmark of five solvers and four sampling schemes for the Ishigami model ($d = 3, p = 14, q = 1$). Solvers are coded by colors. Sampling schemes are shown in varying shades and line styles. In (f) and (g), we show the relative MSE of each of the solvers combined with each sampling scheme in the order MC-LHS-coh-opt-D-opt(coh-opt).
<figure><img src="image_36.png" /><figcaption>**Figure 3.C.3:** Boxplots of relative MSE from the benchmark of five solvers and four sampling schemes for the borehole model ($d = 8, p = 4, q = 1$). Solvers are coded by colors. Sampling schemes are shown in varying shades and line styles. In (f) and (g), we show the relative MSE of each of the solvers combined with each sampling scheme in the order MC-LHS-coh-opt-D-opt(coh-opt).</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>**Figure 3.C.4:** Boxplots of relative MSE from the benchmark of four solvers and four sampling schemes for the two-dimensional diffusion model ($d = 53, p = 4, q = 0.5$). Solvers are coded by colors. Sampling schemes are shown in varying shades and line styles. In (e) and (f), we show the relative MSE of each of the solvers combined with each sampling scheme in the order MC-LHS-coh-opt-D-opt(coh-opt).</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 3.C.5:** Boxplots of relative MSE from the benchmark of four solvers and four sampling schemes for the 100D function ($d = 100, p = 4, q = 0.5$). Solvers are coded by colors. Sampling schemes are shown in varying shades and line styles. In (e) and (f), we show the relative MSE of each of the solvers combined with each sampling scheme in the order MC-LHS-coh-opt-D-opt(coh-opt).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>(a) Ishigami model, LARS</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(b) Borehole model, LARS</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(c) Ishigami model, SP</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(d) Borehole model, SP</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(e) Ishigami model, BCS</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>(f) Borehole model, BCS</figcaption></figure>

**Figure 3.C.6:** Left column: results for the Ishigami model with a smaller basis ($d = 3, p = 12, q = 1$), complementing the plots in Figure 3.5. Right column: results for the borehole model with a smaller basis ($d = 8, p = 4, q = 1$), complementing the plots in Figure 3.6. Results for three sparse solvers and five experimental design schemes. Fifty replications.
<table><caption>**Table 3.D.1:** Overview of some articles and benchmark studies comparing sparse PCE methods. The “best method” listed in the last column is the one delivering the smallest target error, as reported in the respective publications. Target quantities can be moments, Sobol' indices or the generalization error of the PCE surrogate. Abbreviations: ‘>’ stands for ‘better than’; for the acronyms of solvers and sampling schemes see Sections 3.2.5 and 3.2.6 or the respective publications; for the sampling schemes, the method given in parentheses indicates how the corresponding candidate set is created (e.g., D-opt(coh-opt) stands for D-optimal sampling based on a coherence-optimal candidate set). If the cited paper proposed a new method, this method is marked by a star (*).</caption><thead><tr><th>Ref.</th><th>Type</th><th>Methods compared</th><th>Result</th></tr></thead><tbody><tr><td>Hampton and Doostan (2015b,a)</td><td>sampling</td><td>MC, asymptotic*, and coh-opt*</td><td>coh-opt best</td></tr><tr><td>Fajraoui et al. (2017)</td><td>sampling</td><td>LHS, Sobol, D-opt(LHS), S-opt(LHS)*; sequential sampling*</td><td>sequential S-opt(LHS) best; D-opt(LHS) worst</td></tr><tr><td>Hadigol and Doostan (2018)</td><td>sampling (OLS)</td><td>MC, LHS, coh-opt, A-opt(coh-opt)*, D-opt(coh-opt)*, E-opt(coh-opt)*</td><td>D-opt(coh-opt) best for p > d; MC, LHS best for p < d</td></tr><tr><td>Jakeman et al. (2017)</td><td>sampling</td><td>MC, asymptotic, CSA*</td><td>for p > d: CSA better than MC and asymptotic</td></tr><tr><td>Alemazkoor and Meidani (2018a)</td><td>sampling</td><td>MC, coh-opt, near-opt(coh-opt)*</td><td>near-opt(coh-opt) &gt; coh-opt &gt; MC</td></tr><tr><td>Diaz et al. (2018)</td><td>sampling</td><td>coh-opt, D-opt(coh-opt), sequential D-opt(coh-opt)*</td><td>seq. D-opt(coh-opt) &gt; D-opt(coh-opt) &gt; coh-opt</td></tr><tr><td>Dutta and Gandomi (2020)</td><td>sampling</td><td>MC, LHS, Sobol, Importance Sampling</td><td>LHS best</td></tr><tr><td>Hu and Ludkovski (2017)</td><td>solvers</td><td>OMP, SPGL1, BCS (Babacan et al., 2010)</td><td>BCS &gt; OMP &gt; SPGL1</td></tr><tr><td>Huan et al. (2018)</td><td>solvers</td><td>l1_1s, SpaRSA, CGIST, FPC_AS, ADMM with default parameters</td><td>all showed similar performance; ADMM slightly advantageous</td></tr><tr><td>Liu et al. (2020b)</td><td>solvers</td><td>OMP, LARS, rPCE*</td><td>rPCE &gt; LARS &gt; OMP</td></tr><tr><td>Baptista et al. (2019)</td><td>solvers</td><td>OMP, SPGL1, and two variants* of OMP (modified regressor selection, randomization)</td><td>best: OMP, and OMP with modified regressor selection</td></tr><tr><td>Tarakanov and Elsheikh (2019)</td><td>solvers</td><td>OMP, LARS, Rank-PCE*</td><td>best: Rank-PCE</td></tr><tr><td>Zhou et al. (2019c)</td><td>solvers</td><td>LARS, BCS (Wipf and Rao, 2004), BCS (Ji et al., 2008), D-MORPH-reweighted Cheng and Lu (2018b), stepwise regression based on Bayesian ideas*</td><td>best: stepwise regression (based on 1 Sobol' design)</td></tr></tbody></table>
CHAPTER 4

# Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications

This chapter contains the postprint of

Lüthen, N., Marelli, S., Sudret, B. (2022). Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications. *Int. J. Uncertainty Quantification* 12(3), 49–74.<fnref n="1" />

published in the *International Journal for Uncertainty Quantification*, differing from the published paper only in terms of layout and formatting.

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

## Abstract

Sparse polynomial chaos expansions (PCE) are an efficient and widely used surrogate modeling method in uncertainty quantification for engineering problems with computationally expensive models. To make use of the available information in the most efficient way, several approaches for so-called basis-adaptive sparse PCE have been proposed to determine the set of polynomial regressors (“basis”) for PCE adaptively.

<footnote n="1">¹First published in *International Journal for Uncertainty Quantification* in Volume 12, Number 3, 2022, published by Begell House, Inc. Copyright © 2022 by Begell House, Inc.</footnote>
The goal of this paper is to help practitioners identify the most suitable methods for constructing a surrogate PCE for their model. We describe three state-of-the-art basis-adaptive approaches from the recent sparse PCE literature and conduct an extensive benchmark in terms of global approximation accuracy on a large set of computational models. Investigating the synergies between sparse regression solvers and basis adaptivity schemes, we find that the choice of the proper solver and basis-adaptive scheme is very important, as it can result in more than one order of magnitude difference in performance. No single method significantly outperforms the others, but dividing the analysis into classes (regarding input dimension and experimental design size), we are able to identify specific sparse solver and basis adaptivity combinations for each class that show comparatively good performance.

To further improve on these findings, we introduce a novel solver and basis adaptivity selection scheme guided by cross-validation error. We demonstrate that this automatic selection procedure provides close-to-optimal results in terms of accuracy, and significantly more robust solutions, while being more general than the case-by-case recommendations obtained by the benchmark.

## 4.1 Introduction

Surrogate modeling techniques are a popular tool in applied sciences and engineering, because they can significantly reduce the computational cost of uncertainty quantification analysis for costly real-world computational models. Here, the computational model is approximated by a cheaper-to-evaluate function, which is created based on a small number of model evaluations, the so-called experimental design. One well-known and popular surrogate modeling technique is polynomial chaos expansion (PCE), which approximates the output of a computational model by a spectral expansion in terms of an orthonormal polynomial basis in the input random variables (Xiu and Karniadakis, 2002). PCE is particularly well suited for surrogating smooth models in low to medium dimension, and for the efficient computation of moments and Sobol' indices (Sudret, 2008; Le Gratiet et al., 2017). Engineering models are often challenging due to their computational cost: complex models often depend on a large number of input parameters, but we can only afford a few tens or hundreds of model evaluations. This is the so-called low-data regime. Sparse PCE techniques, which aim to compute an expansion involving only few terms, have proven especially powerful and cost-efficient for real-world engineering problems such as, among many others, surrogate-assisted robust design optimization (Chatterjee et al., 2019), hybrid simulation for earthquake engineering (Abbiati et al., 2021), dam engineering (Guo et al., 2019; Hariri-Ardebili and Sudret, 2020), and wind turbine design (Slot et al., 2020). Note that real-world applications are typically not exactly sparse; however, sparse regression-based PCE is a useful tool to find good approximations at low computational cost. The aim in this context is not sparsity, but the accurate approximation of the computational model. In particular, we are interested in global approximation accuracy, as measured by the relative mean-squared error. Applications with specific requirements on accuracy, such as optimization or reliability analysis, might rely instead on more specialized surrogate techniques.

In the last decade, a large amount of literature on sparse PCE has been published, proposing methods that make sparse PCE more accurate, efficient and applicable to high-dimensional
problems. However, it is often not obvious how well these methods perform when compared to and combined with each other, especially on real-world engineering problems. In an attempt to fill this gap, the authors recently conducted a literature survey on sparse PCE and a classification of the available methods into a general framework, as well as a benchmark of selected methods on a broad class of computational models (Lüthen et al., 2021). This benchmark extensively compared sparse regression solvers and experimental design sampling schemes, using a fixed polynomial basis to focus on the effect of those two classes of methods. We found that the choice of sparse solver and sampling scheme can make a difference of up to several orders of magnitude in the resulting approximation error, and that different methods are preferable in different regimes of ED size. The performance of solvers and sampling schemes seemed to be mostly independent from one another.

The goal of the present paper is to build on and complement this earlier benchmark by exploring the promising field of *basis-adaptive sparse PCE*, in which the basis is iteratively augmented and refined. So far, novel basis-adaptive methods have been proposed in isolation and not been compared to one another. We want to help practitioners choose the most suitable methods for constructing a PCE surrogate for their applications by answering the following questions: (1) Is there a significant difference between different combinations of sparse solvers and basis adaptivity strategies, and does the proper choice matter in actual applications? (2) In the case that no combination clearly emerges as superior, is there a smart strategy to automatically select a good combination? To answer the first question, we describe three basis-adaptive schemes from the sparse PCE literature in detail, namely, degree and q-norm adaptivity (Blatman and Sudret, 2011; Marelli et al., 2021a), forward neighbor degree adaptivity (Jakeman et al., 2015), and anisotropic degree basis adaptivity (Hampton and Doostan, 2018). We then evaluate the performance and synergies of combinations of these basis-adaptive schemes with several sparse regression solvers in terms of validation error<fnref n="2" /> on a library of 11 computational models of varying complexity, representative of a broad class of engineering models. These range from three- to 100-dimensional and include both analytical and numerical models. Since no combination of solver and basis adaptivity significantly outperforms the others, we address the second question by introducing an additional selection step, choosing one among several candidate PCEs computed by different combinations of methods on the same experimental design, using a cross-validation estimate for the generalization error.

The paper is structured as follows. In Section 4.2, we recall the definition of (sparse) PCE and the computing framework introduced in Lüthen et al. (2021). We discuss various estimators for the generalization error, sparse regression solvers, and basis adaptivity, and introduce automatic selection. The associated benchmark results for basis adaptivity and automatic selection are presented in Section 4.3. Finally, a discussion and summary of our results is provided in Section 4.4. Additional information and results can be found in the Appendix.

<footnote n="2">²We focus on global approximation accuracy instead of the computational cost of training, since the main cost in surrogate modeling for engineering applications is typically incurred by the model evaluations themselves. Indeed, the training of each of the surrogates presented in this work requires between a few seconds and a few minutes on a standard business laptop. Even orders of magnitude difference in the training costs of each surrogate do not matter as long as the resulting accuracy is appreciably better.</footnote>
## 4.2 Sparse polynomial chaos expansions

Let $\mathbf{X}$ be a $d$-dimensional random vector with mutually independent components and probability density function $f_X$. Denote by $\mathcal{D}_X$ the domain of the random vector $\mathbf{X}$. Define the space $L^2_{f_X}(\mathcal{D}) = \{h : \mathcal{D}_X \to \mathbb{R} \mid \mathrm{Var}_X[h(\mathbf{X})] < +\infty\}$ of all scalar valued models with finite variance under $f_X$. Under very general conditions on the input distribution $f_X$, there exists a polynomial basis $\{\psi_\alpha : \alpha \in \mathbb{N}^d\}$ of $L^2_{f_X}(\mathcal{D})$ (Xiu and Karniadakis, 2002; Ernst et al., 2012). Since the components of $\mathbf{X}$ are mutually independent, each polynomial basis function can be built as a product of univariate polynomials in $X_1, \dots, X_d$ and characterized by the multi-index $\alpha \in \mathbb{N}^d$ whose entries are equal to the degrees of the univariate terms.

For a computational model $\mathcal{M} \in L^2_{f_X}(\mathcal{D})$, let $Y = \mathcal{M}(\mathbf{X})$ denote the output random variable. $Y$ can be cast as the following spectral expansion:

$$Y = \mathcal{M}(\mathbf{X}) = \sum_{\alpha \in \mathbb{N}^d} c_{\alpha} \psi_{\alpha}(\mathbf{X}). \qquad (4.1)$$

In practice, a finite, truncated polynomial chaos expansion

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

is computed, where $\mathcal{A} \subset \mathbb{N}^d$ is the truncation set defining the basis elements used in the expansion. The accuracy of the expansion depends on $\mathcal{A}$ and the coefficient values $(c_\alpha)_{\alpha \in \mathcal{A}} =: c \in \mathbb{R}^P$, with $P = \mathrm{card}(\mathcal{A})$.

Among several methods for computing the coefficient vector $c$ for a given truncation set, sparse regression is a particularly powerful and efficient method (Doostan and Owhadi, 2011; Blatman and Sudret, 2011). In this approach, the model is evaluated at a number of points $\mathbf{X} = \{\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(N)}\} \subset \mathcal{D}_X$ called the experimental design (ED), yielding the vector of model responses $\mathbf{y} = (\mathcal{M}(\mathbf{x}^{(1)}), \dots, \mathcal{M}(\mathbf{x}^{(N)}))^T$. Let $\{\alpha_j\}_{j=1}^P$ be an arbitrary ordering of the multi-indices in the truncation set and define the regression matrix $\mathbf{\Psi} \in \mathbb{R}^{N \times P}$ with entries $\Psi_{ij} = \psi_{\alpha_j}(\mathbf{x}^{(i)})$. Sparse regression methods determine a vector $c$ that minimizes the residual norm $\|\mathbf{\Psi}c - \mathbf{y}\|_2$ under the constraint that it has only few nonzero entries i.e., it is sparse. This is usually achieved by regularized regression, resulting e.g. in the LASSO formulation

$$\hat{c} = \min_{c \in \mathbb{R}^P} \|\mathbf{\Psi}c - \mathbf{y}\|_2^2 \quad \text{s.t.} \quad \|c\|_1 \le \tau, \qquad (4.3)$$

where $\tau$ is a parameter regulating the sparsity of $c$. A PCE with a sparse coefficient vector $c$ is called sparse PCE. Provided that the regression matrix fulfills certain properties (Candès and Wakin, 2008; Bruckstein et al., 2009; Candès and Plan, 2011), sparse regression can find robust solutions to underdetermined systems of linear equations, which means that the experimental design can be smaller than the number of unknown coefficients.

The quality of the solution depends on the choice of the basis $\mathcal{A}$, on the experimental design $\mathbf{X}$, and on the method used for computing the coefficients. For each of these components, many different methods have been proposed in recent years, including iterative methods which adaptively determine $\mathcal{A}$, or the experimental design $\mathbf{X}$. These methods were recently surveyed
<figure><img src="image_45.png" /><figcaption>**Figure 4.1:** Framework for computing sparse PCE introduced by Lüthen et al. (2021), who conducted a literature survey and a benchmark for the central components “Experimental design” and “Sparse regression” (adapted from Lüthen et al. (2021)). In the present work, we discuss and benchmark the components marked in orange. In particular, we discuss the component “Basis adaptivity” and explore its relationship to sparse regression solvers.</figcaption></figure>

and classified into the framework shown in Fig. 4.1 (modified from Lüthen et al. (2021)). In the present contribution, we focus on the question of how to determine a suitable basis $\mathcal{A}$. To this end, we benchmark several basis-adaptive approaches and explore their interplay with sparse regression solvers.

### 4.2.1 Error estimation and model selection

Model selection is applied on several levels in the sparse PCE framework:

* Within the sparse solver to select its hyperparameter (see Section 4.2.2)

* Within the basis adaptivity scheme to select a basis (see Section 4.2.3)

* Finally, between solvers and basis adaptivity schemes, to automatically select a combination that is close to best (see Section 4.2.4).

Our main quantity of interest is the generalization error, in the form of the relative mean squared error normalized by the model variance

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

It can be approximated by the validation error in the form of the relative mean squared error (MSE)

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

computed on a large validation set $\{\mathbf{x}_{\text{val}}^{(i)}\}_{i=1}^{N_{\text{val}}} \sim \text{i.i.d. } f_{\mathbf{X}}$. We only consider model selection criteria that are estimates of the generalization error.

To get an accurate estimate of the generalization error, without using an additional validation set, a widely used method is cross-validation (Vapnik, 1995). Here, the available data is repeat-
edly divided into two disjoint sets, one for computing the solution (training) and the other for evaluating the error (validation). Aggregating the error estimates from the different splits, we get a cross-validation estimate of the generalization error.

One cross-validation strategy is to divide the data randomly into $k$ disjoint, roughly equally large parts, and use each of the parts in turn as validation set. This is called *k-fold cross-validation*. If $k$ is chosen to be equal to the size of the experimental design, the strategy is known as *leave-one-out cross-validation* (LOO). This is closest to using all data points to compute a solution, since in each iteration, only one point is left out from the training set.

In general, LOO can be quite costly, since for an experimental design of size $N$, the method under consideration has to be applied $N$ times. A cheap approximation to the LOO error, which requires only one application of the method, is available for certain sparse regression solvers, namely for those which in their last step recompute the solution with ordinary least-squares (OLS) on the set of regressors with nonzero coefficient (called *active basis*). (Blatman and Sudret, 2010). In particular, this is the case for the solvers hybrid least angle regression (LARS) (Blatman and Sudret, 2011), orthogonal matching pursuit (OMP) (Pati et al., 1993; Jakeman et al., 2015), and subspace pursuit (SP) (Diaz et al., 2018) (see Section 4.2.2). For these solvers, the OLS-based LOO estimate coincides with the true LOO error if the following holds: regardless of which experimental design point is left out, the sparse regression solver consistently yields the same active basis.

Since the repeated use of LOO for model selection often results in the generalization error being underestimated, especially on small experimental designs, Blatman and Sudret (2011) proposed to use a modification factor originally developed for the empirical error (Chapelle et al., 2002), defined by

$$T = \frac{N}{N - P_{\text{active}}} \left( 1 + \text{tr}((\Psi_{\text{active}}^{T} \Psi_{\text{active}})^{-1}) \right), \quad (4.6)$$

where $P_{\text{active}}$ denotes the number of nonzero coefficients (the corresponding basis functions are called *active*), and $\Psi_{\text{active}}$ denotes the regression matrix containing only the active regressors. The product of the modification factor $T$ with the LOO error is called *modified LOO error*.

### 4.2.2 Sparse regression solvers

Various sparse regression solvers available for solving sparse PCE were described in detail in Lüthen et al. (2021). We give here only a short overview of the solvers used in our benchmark, and refer to Lüthen et al. (2021) for further details. These solvers are common choices in the sparse PCE literature.

* Hybrid least angle regression (LARS) (Blatman and Sudret, 2011): adding regressors one-by-one, following a least-angle strategy. The hybrid approach recomputes the final coefficient values by ordinary least squares (OLS) on the selected regressors.

* Orthogonal matching pursuit (OMP) (Pati et al., 1993; Jakeman et al., 2015): greedily adding orthogonal regressors one-by-one, computing their coefficients by OLS.

* Subspace pursuit (SP) (Diaz et al., 2018): searching iteratively for a solution with a certain $l^0$-norm by adding and removing regressors from the active set. Coefficients are computed
by OLS. We include two variants of SP in this benchmark, one which determines its
hyperparameter by 5-fold cross-validation similar to the implementation in (Diaz, 2018),
which we denote by SP<sub>k=5</sub>, and one where it is determined by OLS-based LOO, introduced
in (Lüthen et al., 2021) and denoted by SP<sub>LOO</sub>.

* Bayesian compressive sensing (BCS) (Babacan et al., 2010): using a Bayesian framework to enforce sparsity of the coefficients.

* SPGL1 (van den Berg and Friedlander, 2008): a convex optimization solver following the Pareto front of the residual-sparsity trade-off.

Each of the solvers features at least one hyperparameter whose value needs to be determined
via cross-validation in order to get a good solution. For LARS, OMP, SP<sub>k=5</sub>, and SP<sub>LOO</sub>,
this hyperparameter is the number of active basis functions (nonzero coefficients) of the final
solution. For BCS and SPGL1, it is the bound on the residual norm in the sparse regression
formulation.

The benchmark in Lüthen et al. (2021) of sparse regression solvers on a non-adaptive polynomial basis showed that BCS and SP<sub>LOO</sub> generally outperform other sparse solvers for low-dimensional models, while for high-dimensional models, BCS is by far the best sparse solver.

4.2.3 Basis adaptivity

The sparse solver and the experimental design, which were benchmarked in (Lüthen et al., 2021),
are not the only ingredients to a sparse PCE. The choice of the *candidate basis*, from which the
sparse solver determines the *active basis* (i.e., the set of regressors with nonzero coefficient), is
another important ingredient: if the candidate basis is chosen too small, important terms might
be missing, which might lead to a large model error. On the other hand, if the candidate basis
is too large, the ratio of the number of experimental design points to the number of coefficients
is small, which causes some properties of the regression matrix to deteriorate and can result in
a large approximation error.

Of course, it is not possible to know a-priori the best choice of the candidate basis. Basis-adaptive
schemes start with an initial candidate basis and adapt it iteratively, adding or removing basis
functions in each iteration according to various heuristics. The performance of the bases in the
different iterations is evaluated using a model selection criterion, typically an estimate of the
generalization error.

Several procedures for basis-adaptive sparse PCE have been proposed in the literature. We
discuss three approaches in detail, namely

* degree and q-norm ("p&q") basis adaptivity, as implemented in UQLab (Marelli et al., 2021a), see Section 4.2.3.1;

* forward neighbor basis adaptivity (Jakeman et al., 2015), see Section 4.2.3.2; and

* anisotropic degree basis adaptivity (Hampton and Doostan, 2018), from the software BASE_PC (Hampton and Doostan, 2017), see Section 4.2.3.3.

We also briefly mention several other approaches found in the literature.
### 4.2.3.1 Degree and q-norm (“p&q”) adaptivity

A typical choice for a PCE basis is the *basis of total degree p* defined by the set of multi-indices

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

Furthermore, *hyperbolic truncation* (Blatman and Sudret, 2011) uses the *q-(quasi-)norm*

$$ \|\boldsymbol{x}\|_q = \left( \sum_{i=1}^{d} |x_i|^q \right)^{\frac{1}{q}} \qquad (4.8) $$

with $q \in (0, 1]$ to truncate a total-degree basis further:

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

Hyperbolic truncation has the effect of excluding some of the basis functions with high degree and high interaction order. This is particularly effective for high-dimensional problems.

A simple basis-adaptive scheme is *degree adaptivity* (Blatman and Sudret, 2011), which computes a number of PCEs on a sequence of total-degree candidate bases of increasing size, and returns the PCE minimizing a certain error estimate as final solution. Analogously, a q-norm adaptive scheme can be developed, and easily be combined with degree adaptivity (Marelli et al., 2021a), yielding *degree and q-norm (p&q) adaptivity*. Degree and q-norm adaptivity is solution-agnostic, i.e., it does not take the solution computed in the previous iteration into account.

### 4.2.3.2 Forward neighbor basis adaptivity

Jakeman et al. (2015) suggest a basis-adaptive algorithm based on a graph view of the PCE regressors (see also Gerstner and Griebel (2003); Narayan et al. (2014)). Since the input random vector is assumed to have independent components, the basis functions have tensor product structure. The basis functions can be considered nodes of a directed graph constructed as follows (see also Fig. 4.2): two regressors are considered neighbors if their multi-index of degrees differs only in one dimension by 1, i.e., there is a directed edge from $\psi_\alpha$ to $\psi_\beta$ iff $\gamma := \beta - \alpha$ is a multi-index with $\gamma_i = 1$ for one $i \in \{1, \dots, d\}$ and $\gamma_j = 0, j \neq i$. Forward neighbors of a regressor are regressors reachable by an outgoing edge, and backwards neighbors are regressors connected by an incoming edge.

In the context of basis-adaptive PCE, this construction is used to determine a number of candidate regressors to be added to the currently active basis, starting from an initial basis of small total degree. Assume that an important high-degree or high-order regressor is not yet part of the candidate basis but some of its backwards neighbors are. The fact that it is missing should be visible in the coefficients of its backwards neighbors, which can be expected to have a significant nonzero coefficient to compensate for the fact that the important high-degree regressor is not yet part of the basis.

This heuristic is the foundation of the algorithm whose pseudocode is summarized in Algorithm 6. In each iteration, the current set of active basis functions is determined (restriction
<figure><img src="image_46.png" /><figcaption>**Figure 4.2:** Illustration of forward neighbor basis adaptivity (Jakeman et al., 2015) in $d = 3$ dimensions. The regressors are seen as nodes in a graph, where a directed edge connects regressors whose multi-index differs by 1 in exactly one dimension. Assume that the set of active regressors $\mathcal{A}^{\text{active}}$ initially consists of the gray nodes. All admissible nodes are highlighted in green. Together with the gray nodes, they constitute the set $\mathcal{A}^{(1)}$ and would be added to the model in the next expansion step of forward neighbor basis adaptivity (Step 3 in Algorithm 6).</figcaption></figure>

step). All forward neighbors of these active basis functions are surveyed and added to the candidate basis if they are admissible, i.e., if all of their backwards neighbors are in the active set (expansion step). Jakeman et al. (2015) employ several expansion steps to prevent premature termination, and use cross-validation error estimates. We call this algorithm forward neighbor basis-adaptive.

<caption>**Algorithm 6** Forward neighbor basis adaptivity (Jakeman et al., 2015)</caption>

1: Initial PCE (basis chosen a-priori, typically small total-degree basis)
2: Restriction: retain only the active regressors $\mathcal{A}^{\text{active}}$
3: Expansion: let $\mathcal{A}^{(0)} = \mathcal{A}^{\text{active}}$. For $t = 1, \dots, T$, obtain the set $\mathcal{A}^{(t)}$ by augmenting $\mathcal{A}^{(t-1)}$ by all its admissible forward neighbors.
4: Compute a PCE and its error estimate for each augmented basis $\mathcal{A}^{(1)}, \dots, \mathcal{A}^{(T)}$.
5: Choose the PCE with the lowest error estimate among the $T$ candidates. Stop if this error estimate is larger than the previously obtained best error estimate. Else, continue iteration with restriction step (Step 2)

This algorithm is implemented in the software Dakota (Adams et al., 2014). We use our own Matlab-based implementation of the algorithm. Consistently with Jakeman et al. (2015), we set $T = 3$.

### 4.2.3.3 Anisotropic degree basis adaptivity

Hampton and Doostan (2018) propose an algorithm called BASE-PC, which combines a basis-adaptive scheme with sequential enrichment of a coherence-optimal experimental design. The two steps are executed alternatingly: based on the current ED, a suitable basis is chosen; and according to the currently active basis functions, the ED is enriched and the weights are updated.

In the present work, we solely consider the basis adaptivity part of the BASE-PC algorithm. The main idea of this algorithm is to use an *anisotropic degree basis*, defined by a *degree vector* $\mathbf{p} \in \mathbb{N}^d$. A related idea was explored earlier by Blatman and Sudret (2009). Similarly to a
total-order basis, an anisotropic degree basis is defined by

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

If all entries of $\mathbf{p}$ are the same, i.e., $p_1 = p_2 = \dots = p_d = p$, this definition reduces to a total-order basis of degree $p$. The equation $\sum_{i=1}^{d} \frac{\alpha_i}{p_i} = 1$ defines a hyperplane that cuts the $i$-th coordinate axis at $\alpha_i = p_i$.

In each iteration, the algorithm determines the current anisotropic degree vector based on the currently active basis. A new, larger candidate basis is then constructed by considering the anisotropic degree basis corresponding to a uniformly increased anisotropic degree vector.

We use our own, slightly modified implementation of the basis adaptivity part of BASE-PC, as summarized in Algorithm 7. The most costly operation is the computation of the anisotropic-degree basis (line 6). Hampton and Doostan (2018) developed a specialized efficient algorithm to generate the multi-indices of an anisotropic degree basis, which we utilize in our implementation. We call Algorithm 7 anisotropic degree basis-adaptive.

<caption>**Algorithm 7** Anisotropic degree basis adaptivity (Hampton and Doostan, 2018)</caption>

1: Initial PCE (fixed basis of low order)
2: **for** $o = 1, \dots, 10$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Restriction: denote by $\mathcal{A}^{\text{active}}$ the set of active regressors of the last selected PCE
4: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1, \dots, 10$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Additional restriction: remove $\frac{i-1}{10}$ regressors from the set $\mathcal{A}^{\text{active}}$ starting with the ones with smallest-in-magnitude coefficients, obtaining a set $\mathcal{A}^i \subseteq \mathcal{A}^{\text{active}}$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Expansion: compute the dimension-wise maximal degree of the regressors in $\mathcal{A}^i$, denoted by degree vector $\mathbf{p}^{\text{max}}$. Compute $\mathbf{p}^{\text{new}} = (p_1^{\text{max}} + 1, p_2^{\text{max}} + 1, \dots, p_d^{\text{max}} + 1)$ and the associated anisotropic-degree basis $\mathcal{A}^{i, \text{new}}$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute a PCE on the basis $\mathcal{A}^{i, \text{new}}$ and its error estimator $e^i$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;If $e^i \ge e^{i-1}$, increase the so-called strike counter by 1. Break from the inner loop if the strike counter is $\ge 3$
9: **end for**
10: From the 10 candidate PCEs, select the PCE with the lowest error estimate
11: Break from the outer loop if the error estimate has increased three times in subsequent iterations (ExpansionEarlyStop; same idea as inner loop strike counter)
12: **end for**
13: Return the PCE with the lowest error estimate among all PCEs selected in the outer loop.

### 4.2.3.4 Other basis adaptivity algorithms

We briefly summarize other approaches for basis adaptivity for sparse PCE. These approaches will not be investigated in our benchmark.

Many algorithms which use stepwise regression to build up a sparse basis regressor-by-regressor can be classified both as sparse regression solvers (as done in Lüthen et al. (2021)) and as basis-adaptive algorithms. As an example, the approach by Blatman and Sudret (2010) adds and
removes regressors one-by-one based on the induced change in LOO error, thereby building up
a small set of active basis functions which is sparse in a larger total-degree basis. In this case,
the candidate basis coincides with the active basis. Mai and Sudret (2015) use the “principle of
heredity” together with LARS to determine additional bivariate interaction terms to be added to
the model, once a univariate term is identified as relevant. We do not consider these approaches,
since we are interested in algorithms modifying the candidate basis not only one regressor at a
time, but at a larger scale.

Alemazkoor and Meidani (2017) propose a basis-adaptive algorithm relying on sparsity and step-by-step augmentation of the basis. In each step, the candidate basis is a total-order basis for a subset of input random variables, while the remaining input random variables are considered constant. The initial basis consists only of the constant term. In each iteration, either one of the constant dimensions is activated, or the total degree of the candidate basis is increased by one (active dimensions only), depending on the resulting error estimate or sparsity of the solution. The coefficients are solved by OLS until a pre-defined threshold of residual norm is reached. Then, the sparse regression solver SPGL1 is used, which identifies the sparsest solution whose residual norm is smaller than the given threshold. We do not consider this method due to its high computational cost and its strong tie to the solver SPGL1, making it less effective when paired with other sparse regression solvers.

Loukrezis et al. (2019) propose a basis-adaptive algorithm for interpolating PCE on Leja-based
sparse grids. Their algorithm is a more cautious version of forward neighbor basis adaptivity
(Section 4.2.3.2): after adding all admissible forward neighbors to the basis and computing
their coefficients, all of the recently added regressors are removed again, except for the one that
achieved the largest-in-magnitude coefficient. Thapa et al. (2020) suggest a basis-adaptive pro-
cedure that relies on total-order bases of increasing degree. In contrast with degree adaptivity
(Section 4.2.3.1), the basis functions of degree $p+1$ are not added all at once, but in chunks of
a certain size dependent on the dimension and the degree $p$. After adding a chunk of basis func-
tions, the PCE is recomputed and regressors with a coefficient smaller than a certain threshold
are removed from the basis. We do not consider these two approaches because they are similar
to the previously presented approaches while being more costly.

### 4.2.4 Automatic selection of a sparse PCE solution from several candidate solutions

For realistic simulators used in engineering and applied sciences, evaluating the computational
model is the expensive part of the surrogate modeling process. Once the model evaluations are
obtained, all further computations are post-processing steps, and are computationally cheap in
comparison to the model evaluations. Thus, it is feasible to apply several adaptive sparse PCE
methods and choose the best one.

We therefore propose the use of an additional layer of selection which we call here *automatic selection*. It can be seen as a simple case of *ensemble modeling* (Sagi and Rokach, 2018) where a single metamodel (“strong learner”) is chosen based on its cross-validation performance – or, in other words, model selection. For a given experimental design, we compute several sparse
PCEs through different combinations of sparse solvers and basis adaptivity schemes. From the resulting PCE solutions, we choose one based on the value of a suitable estimate of generalization error as model selection criterion. There are several possibilities. One possible choice is the estimator already used for selecting the hyperparameter of the solver and for basis adaptivity, i.e., (modified) LOO for LARS, OMP and $SP_{LOO}$, and $k$-fold cross-validation for $SP_{k=5}$, BCS and SPGL1. However, this estimator might not be consistent between solvers, which is however necessary for such an automatic selection. Furthermore, this estimate might be biased due to its repeated use.

A second class of estimators is given by the so-called *hybrid cross-validation error estimators*. The word *hybrid* is a reference to Efron et al. (2004), who created the hybrid version of least-angle regression (LARS) which uses LARS only to identify the set of active basis functions, and then recomputes the coefficients by ordinary least-squares (OLS) on this active set. To compute the hybrid leave-one-out (LOO) or hybrid $k$-fold cross-validation error estimate for any PCE solution, the same procedure is used: first, the active basis functions are identified using the whole experimental design. Then, the coefficients are recomputed several times using OLS, following the chosen cross-validation framework. This requires solving a linear system of equations $k$ times in the case of $k$-fold cross-validation. In case of LOO, only one linear system of equations needs to be solved (Blatman and Sudret, 2010). Furthermore, we can make use of the modification factor of Eq. (4.6) to compute the hybrid modified LOO error estimate.

As baseline cases, we will also select a solution 1) randomly from a set of generally well-performing methods, and 2) corresponding to the best combination identified in the benchmark of basis adaptivity schemes (Section 4.3.1).

## 4.3 Numerical results

In this benchmark, we compare the performance of various combinations of sparse regression solvers with basis-adaptive schemes. We use the following methods and associated implementations. The sparse solvers (wrapped to fit into our benchmark framework) are:

* Hybrid-LARS: UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* OMP: UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* $SP_{k=5}$: own implementation of Diaz et al. (2018), available in UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* $SP_{LOO}$: own adaptation of Diaz et al. (2018), available in UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* BCS: own implementation of FastLaplace (Babacan et al., 2010), available in UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* SPGL1: SPGL1 v1.9³ (van den Berg and Friedlander, 2008; Van den Berg and Friedlander, 2015)

The basis adaptivity schemes are:

<footnote n="3">³Our benchmark is performed with SPGL1 v1.9. A new version of SPGL1, v2.1, is available since June 2020. In our tests, the new version (with “hybrid” mode) did not perform significantly better than the old version. The numerical results show therefore results for SPGL1 v1.9 with default parameters.</footnote>
* p&q adaptivity: UQLab (Marelli et al., 2021a; Marelli and Sudret, 2014)

* forward neighbor basis adaptivity: own implementation of the algorithm based on the description in Jakeman et al. (2015)

* anisotropic degree basis adaptivity: own implementation of an algorithm adapted from the basis-adaptive part of BASE-PC (v1) (Hampton and Doostan, 2018), using the function `basis_anisotropic_total_order.m` from Hampton and Doostan (2017) to generate the multi-indices of an anisotropic degree basis (see Section 4.2.3.3)

To reduce the complexity of our benchmark, we choose Latin hypercube sampling with max-
imin distance optimization to generate the experimental design (ED) since Lüthen et al. (2021)
demonstrated that the choices of sparse solver and sampling scheme are mostly independent
from each other.

We use the following model selection methods:

* For the selection of the hyperparameters of the sparse regression solvers, we use
  - modified OLS-based LOO for the solvers LARS, OMP and SP<sub>LOO</sub>
  - k-fold CV for the solvers SP<sub>k=5</sub> (k = 5), BCS (k = 10) and SPGL1 (k = 10)

* The basis adaptivity schemes use the same criterion as the respective solver uses.

* We investigate in Section 4.3.2 which criterion is suited best for the final model selection.

I.e., instead of prescribing fixed values for the hyperparameters, we let the parameters be determined adaptively. In this sense, we give each method equal opportunity to produce the best possible solution (assuming that the CV error is a good proxy for solution quality).

For our benchmark, we use 11 benchmark models ranging from low-dimensional analytical mod-
els to high-dimensional differential equations. All of these models have previously been used
as numerical examples for surrogate modeling or reliability analysis. None of these models has
an exactly sparse representation. Note that we do not aim at benchmarking the ability of the
methods to recover “true” sparse solutions, but instead their approximation capabilities on en-
gineering models, which are typically not exactly sparse but *compressible*, meaning that the
magnitude of their PC coefficients decays rapidly. Table 4.1 provides an overview of the bench-
mark models. For a more detailed presentation, we refer the interested reader to the respective
publications (see column "Reference" of Table 4.1).

4.3.1 Basis adaptivity

We benchmark the sparse solvers LARS, OMP, SP<sub>k=5</sub>, SP<sub>LOO</sub>, BCS, and SPGL1 combined with
the basis-adaptive schemes described in Section 4.2.3:

1. degree and q-norm (p&q) adaptivity (abbreviation: PQ)

2. forward neighbor basis adaptivity (FN)

3. anisotropic degree basis adaptivity (AD)

As a base case, we include a *static basis* following the rule $P \approx \frac{10}{3}N$ (abbreviation: ST), where we use hyperbolic truncation with $q = 0.5$ for high-dimensional models ($d \ge 20$).

<footnote n="4">⁴https://www.uqlab.com/sensitivity-high-dimension</footnote>
<table><caption>**Table 4.1:** Overview of the 11 computational models used in our benchmark. Finite element models are marked in italic font, all other models are analytical. The column “Reference” provides the literature in which the models and their input are described in detail. The column “ED sizes” contains the two sizes of experimental design (small and large) used in the basis adaptivity benchmark.</caption><thead><tr><th>Model</th><th>Dimension</th><th>Distributions</th><th>Reference</th><th>ED sizes (small, large)</th></tr></thead><tbody><tr><td>Ishigami function</td><td>3</td><td>uniform</td><td>Blatman and Sudret (2011)</td><td>50, 150</td></tr><tr><td>Undamped oscillator</td><td>6</td><td>Gaussian</td><td>Echard et al. (2013)</td><td>60, 120</td></tr><tr><td>Borehole function</td><td>8</td><td>Gaussian, lognormal, uniform</td><td>Morris et al. (1993)</td><td>100, 250</td></tr><tr><td>Damped oscillator</td><td>8</td><td>lognormal</td><td>Dubourg (2011)</td><td>150, 350</td></tr><tr><td>Wingweight function</td><td>10</td><td>uniform</td><td>Forrester et al. (2008)</td><td>100, 250</td></tr><tr><td>Truss model</td><td>10</td><td>lognormal, Gumbel</td><td>Blatman and Sudret (2011)</td><td>100, 250</td></tr><tr><td>Morris function</td><td>20</td><td>uniform</td><td>Blatman and Sudret (2010b)</td><td>150, 350</td></tr><tr><td>Structural frame model</td><td>21</td><td>lognormal, Gaussian; dependent input variables</td><td>Blatman and Sudret (2010)</td><td>150, 350</td></tr><tr><td>2-dim diffusion model</td><td>53</td><td>Gaussian</td><td>Konakli and Sudret (2016)</td><td>100, 400</td></tr><tr><td>1-dim diffusion model</td><td>62</td><td>Gaussian</td><td>Fajraoui et al. (2017)</td><td>100, 400</td></tr><tr><td>100D function</td><td>100</td><td>uniform</td><td>UQLab example<sup>4</sup></td><td>400, 1200</td></tr></tbody></table>

The benchmark is performed on all 11 models presented in Table 4.1. The experimental design (ED) is created by Latin hypercube sampling (LHS) with optimized maximin distance. We investigate one “small” and one “large” ED size per model (see last column of Table 4.1), which correspond to the second-smallest and second-largest experimental design size, respectively, of an earlier benchmark (Lüthen et al., 2021) dedicated to investigating solvers and sampling schemes on a static basis. The small ED size is chosen to be at the lower end of the range of reasonable ED sizes. The large ED size represents the “highly informative” regime (in the engineering sense), which for costly engineering models is in the order of a few hundred model evaluations. Here we restrict ourselves to two ED sizes to control the complexity of the results; however, two ED sizes are needed since the earlier benchmark showed that the solver behavior is different in the two regimes. Since almost no engineering model is exactly sparse in the PCE basis, adaptivity can be expected to help identify relevant higher-order terms in all regimes of ED size – as soon as the ED is large enough to contain some information about the model. For each model and experimental design size, we repeat the analysis $R = 30$ times to account
for the stochasticity of the sampling method. Due to their excessive computational cost<fnref n="5" />, we omit SPGL1 and anisotropic degree basis adaptivity from the benchmark for high-dimensional models ($d \ge 20$). More details on the settings for the basis adaptivity schemes (e.g., initial basis and investigated degree ranges) can be found in 4.A.

### 4.3.1.1 Boxplots of results for the Ishigami function

The results from this benchmark for the Ishigami function are displayed in Fig. 4.1 in the form of boxplots. Results for the remaining models are shown in 4.B.1. The boxplots visualize the results for all combinations on 30 independent ED realizations. The star-shaped markers denote for one (arbitrarily selected) ED realization the respective validation errors of each of the combinations, highlighting one set of data points which is also contained in the boxplots. We observe that the ranking based on statistics (e.g., median as denoted by the white circle) and the ranking based on the actually attained error on a specific ED can be quite different (see e.g. SP$_{k=5}$ & FN vs. BCS & FN in Fig. 4.1a). From the star-shaped markers, it is obvious that the solvers and basis-adaptive schemes do not exploit the available information in the same ways: while some combinations show their best performance on the selected ED, others perform average or worse on the same ED. We want to find the “best” method, i.e., a solver-basis adaptivity combination that consistently, on each different model and ED realization, achieves a close-to-optimal error. In other words, in this comparative study we are less interested in the absolute values of the error (since we assume that the tested methods have been validated before in the associated literature), but rather in the *relative performance* of the methods. The boxplots alone do not give the full picture, since they do not show which errors correspond to the same ED realization. Also, looking at the results for this and other benchmark models in 4.B.1, it is difficult to visually extract information about the overall performance of the methods.

### 4.3.1.2 Bar plots – results aggregated over models and replications, separately per model dimensionality and ED size

Therefore, we propose an aggregated visualization of the results as displayed in Fig. 4.2 and described in the following. For every model and ED size, we determine on each unique ED realization which combination attained the smallest error $\epsilon^*$ (for example, for the ED realization visualized in Fig. 4.1a by the star-shaped markers, SP$_{LOO}$ & PQ attains rank 1). We also determine on the same ED for each combination whether it came within a factor of 2, 5, or 10 of $\epsilon^*$. For example, for the ED realization visualized in Fig. 4.1a by the star-shaped markers, the best error is attained by SP$_{LOO}$ & PQ. All other combinations achieve an error that is more than a magnitude larger, and would therefore not get a count for this ED realization.

We then aggregate the counts over all repetitions and all models in four different classes: low dimensionality ($d \le 10$) and small experimental design, low dimensionality and large experi-

<footnote n="5">⁵Anisotropic degree basis adaptivity increases the degree of the basis uniformly by 1 in each dimension, which for high-dimensional models often results in infeasibly large bases. SPGL1, or the implementation we are using, takes considerably longer than the other solvers for an increasing number of basis elements $P$. Since our benchmark setup requires several thousand runs of the sparse solver for the case of high-dimensional models, the computational cost of the full benchmark is infeasible given our computational resources.</footnote>
<figure><img src="image_47.png" /><figcaption>**Figure 4.1:** Comparison of all possible combinations of solvers and basis adaptivity schemes for the Ishigami function. Abbreviations of basis adaptivity schemes: ST – static basis, PQ – degree- and q-norm, FN – forward neighbor, AD – anisotropic degree. We show validation errors attained by all 24 combinations of methods on 30 realizations of experimental design. The boxplots visualize the attained errors across all 30 realizations. The star-shaped markers denote for one arbitrarily selected ED realization the errors that the different combinations attain (highlighting one set of data points which is also contained in the boxplots). Left: small ED; right: large ED.</figcaption></figure>

mental design, and their corresponding high dimensionality counterparts ($d \ge 20$). In other words, for each of the four classes we count how often each combination achieved the best error or came within a factor of 2, 5, or 10 of the best error, for each unique experimental design. These counts, scaled as percentages of all runs$^6$, are then visualized in the form of overlapping horizontal bars in Figs. 4.2a–4.2e as follows: The dark red bar shows the percentage of runs in which the particular combination of sparse solver and basis-adaptive scheme found the smallest relative MSE $\epsilon^*$. The other three bars in hues of red illustrate how often the respective combination was within a factor of 2, 5, or 10 of $\epsilon^*$ on the same ED.

These bars are therefore a measure of the dispersion of the accuracy of each method in terms of distance to the best attainable accuracy. This measure is more interesting than the absolute scatter of a method, since a large variability does not matter as long as the method always provides a close-to-optimal solution; and small variability is of no value if the method is always orders of magnitude away from the best possible solution. We only show the six combinations of solver and basis adaptivity scheme whose relative MSE was most often within two times the best relative MSE (denoted by the bright red bar), and sort the combinations according to this criterion. (For the full list of combinations, see Figs. 4.B.3–4.B.6 in the appendix.)

Comparing the subplots 4.2a–4.2d, we observe that the results for small and large ED sizes and low- and high-dimensional models are indeed quite different, which justifies analyzing the results separately. As detailed in the previous paragraph, the plots show both which combinations of solver and basis adaptivity scheme attain the smallest relative MSE how often (length of the dark red bar), and how robust the combination is (lengths of the other three bars), i.e., how often the solution was within a small factor of the best solution. We observe that depending on the criterion considered to assess the combinations of solver and basis adaptivity scheme, different combinations turn out to perform best. In the following, we summarize the performance of solver-

<footnote>$^6$There are 6 × 30 realizations for the two low-dimensional cases (Figs. 4.2a and 4.2b), and 5 × 30 realizations for the two high-dimensional cases (Figs. 4.2c and 4.2d). Fig. 4.2e aggregates over all 660 ED realizations.</footnote>
<figure><img src="image_48.png" /><figcaption>(a) low-dimensional models, small ED (6 × 30 runs)</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(b) low-dimensional models, large ED (6 × 30 runs)</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(c) high-dimensional models, small ED (5 × 30 runs)</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(d) high-dimensional models, large ED (5 × 30 runs)</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(e) Aggregated results for all 11 models and both ED sizes (11 × 2 × 30 runs)</figcaption></figure>

**Figure 4.2:** Aggregated results of our comparative study. For each of the 11 models, 2 ED sizes, and 30 replications (i.e., ED realizations), we have run each of the combinations of sparse solver and basis adaptivity scheme and computed the corresponding relative mean-squared errors (MSE) (see Fig. 4.1 for results for the Ishigami function). We determine the relative performance of the method combinations separately on each ED realization: denoting the best attained relative MSE on a specific ED by $\epsilon^*$, we record which of the method combinations reached this error or came within a factor of $\{2, 5, 10\}$ of $\epsilon^*$ on the same ED. We then aggregate the counts over all replications and models, separating by low/high dimensionality and small/large ED size (panels 4.2a-4.2d). This results in 6 × 30 runs for the low-dimensional cases, and 5 × 30 runs for the high-dimensional cases (per method combination). For readability, in panels 4.2a-4.2d only the six best combinations are shown (for the full results, see 4.B.2). The results are displayed in the form of overlapping bars, representing percentages. The dark red bar represents the percentage of runs where the method combination reached the best relative MSE on the specific ED realization. The other bars in shades of red denote the percentages of runs where the error of the particular method combination was within a factor of $\{2, 5, 10\}$ of the best error on the same ED realization. The method combinations are sorted by the length of the bright red bar (number of runs with an error $\le 2\epsilon^*$).

In panel 4.2e, we display the relative performance of the method combinations aggregated over all models, ED sizes and replications (i.e., 11 × 2 × 30 runs).
basis adaptivity combinations with three numbers X-Y-Z, where X denotes the percentage of
runs where the respective combination turned out best (i.e., the length of the dark red bar), Y
denotes the percentage of runs that were within a factor of two of the smallest error on that
ED (length of the bright red bar), and Z denotes the percentage of runs that were within 1
order of magnitude of the smallest error (light red bar). The numbers are rounded to integer
values. Italic numbers indicate that this is the best value achieved among all combinations. The
enumeration tags refer to the panels in Figure 4.2.

(4.2a) *Low-dimensional models, small ED:* BCS together with forward neighbor basis adaptivity (22-56-75) achieves the smallest error more often than any other combination. However, this combination is not the most robust: in 100-75 = 25% of runs, its error is more than one magnitude worse than the best error. SP<sub>LOO</sub> together with p&q adaptivity (11-32-89) is the most robust instead.

(4.2b) *Low-dimensional models, large ED:* by far the best combination in all three categories is SP<sub>LOO</sub> together with forward neighbor basis adaptivity (42-69-96).

(4.2c) *High-dimensional models, small ED:* there are two combinations that outperform the others: SP<sub>k=5</sub> with forward neighbor basis adaptivity (29-74-80) and BCS with a static basis (3-73-100). A close third is BCS with p&q adaptivity (27-69-100).

(4.2d) *High-dimensional models, large ED:* the two best combinations are SP<sub>k=5</sub> with forward neighbor basis adaptivity (53-77-81) and SP<sub>LOO</sub> with forward neighbor basis adaptivity (21-62-99).

We observe from this enumeration that only the solvers BCS, SP<sub>k=5</sub> and SP<sub>LOO</sub> are found among the best combinations. Considering the best six combinations (based on the second number Y) as in Fig. 4.2, only LARS is joining the list. OMP does not perform well in any of the cases (see also Appendix, Figs. 4.B.3–4.B.6). This is likely due to its tendency to severely underestimate the generalization error, which might make the comparison between error estimates of different bases unreliable. Likewise, SPGL1 is never among the best six combinations for low-dimensional models.

Regarding basis adaptivity schemes, the static basis is not among the best six combinations (based on the second number Y) for low-dimensional models. However for high-dimensional models, especially for the small ED case, it is among the six best combinations several times and performs well.

Note that in the high-dimensional, small ED case (Figs. 4.2c and 4.B.4), the static basis shows the most robust behavior (in terms of the third number Z) while any solver with forward neighbor basis adaptivity is off by more than a magnitude in at least 20% of runs. This is due to its bad performance for the 100D function (see also Fig. 4.B.2u). The 100D function is also responsible for the peculiar results of the high-dimensional, large ED case (Figs. 4.2d and 4.B.6): only SP<sub>LOO</sub> together with forward neighbor basis adaptivity is able to find a very accurate solution, reaching a validation error of 10<sup>-10</sup> which is several orders of magnitude smaller than what the other combinations reach. However, several other combinations are within one order of magnitude of the smallest error for all other high-dimensional models, reaching overall 80% for Z, as can be seen from Fig. 4.B.6.
### 4.3.1.3 Bar plots – results aggregated over all models, ED sizes, and replications

To get a complete picture of solver and basis adaptivity scheme performance, we display in Fig. 4.2e the relative performance of all combinations of solvers and basis adaptivity schemes, averaged over all models and experimental designs regardless of their dimensionality and size. We exclude combinations involving SPGL1 and anisotropic degree basis adaptivity (AD), which were only tested for low-dimensional models. Therefore, Fig. 4.2e contains 3 × 5 = 15 combinations. We see that SP<sub>LOO</sub> together with forward neighbor basis adaptivity (26-59-89) performs best in all three categories. SP<sub>k=5</sub> together with forward neighbor basis adaptivity (24-56-81) is on the second place in terms of the first two criteria (X and Y), while SP<sub>LOO</sub> together with p&q basis adaptivity (6-29-86) is the second-best solver in terms of achieving an error within one order of magnitude of the best error most often (Z). BCS together with any basis adaptivity scheme performs well, while all combinations involving OMP are on the bottom of the list. Combinations with a static basis are found more towards the end of the list, and those with forward neighbor basis adaptivity are found more towards the beginning of the list.

### 4.3.1.4 Conclusion of the comparative study

From these plots, we see clearly that there is no single combination of sparse solver and basis adaptivity scheme that always outperforms all others. While SP<sub>LOO</sub> together with forward neighbor basis adaptivity shows superior performance when averaged over all models and ED sizes, we identify better-performing combinations when we differentiate by model dimension and ED size. In some cases, we are faced with a trade-off between accuracy and robustness: e.g., for low-dimensional models and small EDs, the choice of BCS & FN has a higher probability of yielding a near-optimal solution (within a factor of 2 to the best), while the choice of SP<sub>LOO</sub> & PQ has a higher probability of not being more than a magnitude off from the optimal one.

It seems to be easier to identify trends regarding solvers than regarding basis adaptivity schemes. For example, BCS performs well when information is scarce (small ED sizes or high-dimensional models) while SP<sub>LOO</sub> performs well for large ED sizes, as already observed in (Lüthen et al., 2021), and OMP generally underperforms. This might be because sparse regression solvers are based on quite different principles, from gradient descent over greedy stepwise regression to Bayesian reasoning, while the basis adaptivity schemes all work with variations of the same concept – namely, gradually increasing the degree of the basis. Basis-adaptive schemes generally outperform the static basis, most likely because they offer more basis elements to choose from. However, it seems to not matter as much how the additional basis elements are generated, since none of the three basis-adaptive methods always finds the best basis. Note that the selection of the bases is guided by the cross-validation error. If the cross-validation error does not correlate well with the real validation error, the scheme will select a suboptimal basis. OMP generally has an unreliable cross-validation error, which might explain its bad performance with basis-adaptive schemes.
4.3.2 Automatic selection of sparse solver and basis adaptivity scheme

As we saw in the previous section, there is no single best-performing combination of sparse
solver and basis adaptivity scheme. In this section, we investigate the question: is there any
criterion which can help us choose the best combination of solver and basis adaptivity scheme
for a given problem? One such criterion could be a deterministic rule based on dimensionality of
the model and ED size, relying on the results from the previous section. Another option could
be to use cross-validation-based model selection criteria as described in Section 4.2.4 (see also
Section 4.2.1). We call this process of automatically choosing a well-performing combination
from a number of candidate combinations using a model selection criterion automatic selection.

Due to the performance of the methods in the benchmark in the previous sections, we restrict
our investigation to the solvers SP<sub>k=5</sub>, SP<sub>LOO</sub> and BCS, and to the basis adaptivity schemes

* p&q adaptivity (PQ), forward neighbor adaptivity (FN), and anisotropic degree adaptivity (AD) for low-dimensional models

* static basis (static), p&q adaptivity (PQ), and forward neighbor adaptivity (FN) for high-dimensional models

resulting in 9 possible solutions for each model, ED size and repetition.

We use the following model selection criteria (see also Section 4.2.4):

1. the oracle solution, i.e., the smallest relative MSE attained among *all* 24 or 15 combinations of methods (for low-dimensional or high-dimensional models, respectively) – i.e., among all methods tested in Section 4.3.1, not only among the 9 candidate solutions considered here – on each ED realization, as an ideal lower bound. Of course, this information is not available in practice.

2. the criterion used for basis and hyperparameter selection by the respective solver (see Section 4.2.1)

3. hybrid LOO, computed by OLS on the active basis functions only

4. hybrid modified LOO, computed by OLS on the active basis functions only, and using the correction factor from Eq. (4.6)

5. hybrid 10-fold cross-validation error, computed by OLS on the active basis functions only

6. A fixed rule dependent on dimensionality and ED size, according to the findings in Section 4.3.1, choosing the solver that achieved the smallest error most often (i.e., having the longest dark red bar, and coincidentally also the longest bright red bar):

* low-dimensional models, small ED: BCS & FN

• low-dimensional models, large ED: SP<sub>LOO</sub> & FN

* high-dimensional models, small ED: SP<sub>k=5</sub> & FN

* high-dimensional models, large ED: SP<sub>k=5</sub> & FN

This implies that by design, the lengths of the bars for this selection criterion are identical
with the results for the respective combination in Fig. 4.2.

7. A randomly picked combination of solver and basis adaptivity scheme from the 9 available options as upper bound: any sensible model selection criterion must perform better than this.

The results are presented in Figure 4.3. Note that we do not display which sparse solver and
<figure><img src="image_53.png" /><figcaption>(a) low-dimensional models, small ED</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(b) low-dimensional models, large ED</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>(c) high-dimensional models, small ED</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>(d) high-dimensional models, large ED</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>(e) Aggregated over all 11 models and both ED sizes</figcaption></figure>

**Figure 4.3:** Testing different automatic selection strategies. As before, the overlapping bar plot visualizes the percentage of runs where the respective strategy achieved the best error $\epsilon^*$ or was within a factor of {2, 5, 10} of the best error $\epsilon^*$ among the seven selection strategies. The bar labels are explained in the enumeration of model selection criteria in the text.
basis adaptivity scheme was chosen – we only show how close the chosen solution comes to the best possible solution.

* By construction, the oracle selection performs best in all cases. (We only include it into the plot to emphasize that all other criteria are measured against the best attained error on each ED.)

* The random selection is by far the worst selection criterion, which shows that automatic selection provides a statistically significant improvement over random selection.

* All three hybrid selection criteria are in general more robust than the fixed rule, attaining an error within 2 times of the best solution (bright red bar) similarly or more often than the fixed rule, and being within one order of magnitude of the best error (light red bar) even more than 90% of runs (low-dimensional models) or 99% of runs (high-dimensional models). Averaging over all runs (Fig. 4.3e), the hybrid criteria attain the best error almost as often as the fixed rule, which was chosen for this reason, i.e., because it attained the best error most often. This shows that automatic selection is able to return a solution with almost the same accuracy as the best-performing method, while being more robust than a fixed rule or a random choice of methods.

* The criterion used for basis and hyperparameter selection, i.e., k-fold CV for $SP_{k=5}$ and BCS, and LOO for $SP_{LOO}$, performs slightly worse than the three hybrid selection criteria in all cases shown in Fig. 4.3. The slightly worse performance of this selection criterion might be explained by selection bias: since the criterion was already used twice for selection, it is likely that it underestimates the generalization error. Another reason might be that the criteria used by different solvers (OLS-based LOO vs. k-fold CV) might not be comparable among each other (inconsistency).

We conclude that automatic selection, i.e., using cross-validation to choose a PCE solution from a number of candidate solutions computed with different methods on the same experimental design, can achieve results that are as accurate as and more robust than results from a fixed rule for which combination of solver and basis adaptivity scheme to use, even if this rule is based on a thorough benchmark such as the one in Section 4.3.1. While automatic selection leads to increased computational cost due to the need of training the surrogate with different solvers and basis adaptivity schemes on the same experimental design, it is a case-independent strategy applicable to models of any dimensionality and experimental designs of any size, and it finally results in a PCE that is both accurate and more reliable.

## 4.4 Conclusion and discussion

Our goal was to provide guidance for the choice of sparse PCE methods in engineering applications, by answering the two questions: (1) Is there a significant difference between different combinations of sparse solvers and basis adaptivity strategies, and does the proper choice matter in actual applications? (2) Is there a smart strategy to automatically select a good combination?

To answer these questions, we performed an extensive comparative study investigating several approaches for computing sparse PCE with the goal of surrogate modeling, using the relative
mean squared error on a validation set as the main performance metric. In particular, we studied the performance of combinations of basis adaptivity schemes and sparse solvers in order to identify combinations that yield the smallest generalization error. Our investigations are based on 11 analytical and numerical benchmark models of varying input dimension and complexity, representative of a wide range of engineering problems. We considered the sparse solvers least angle regression (LARS), orthogonal matching pursuit (OMP), subspace pursuit based on k-fold cross-validation ($SP_{k=5}$), subspace pursuit based on LOO ($SP_{LOO}$), Bayesian compressive sensing/FastLaplace (BCS), and spectral projected gradient-$\ell^1$ (SPGL1). The basis adaptivity schemes we compared were a fixed basis truncation scheme following the rule $N \approx \frac{10}{3}P$, degree- and q-norm adaptivity, forward neighbor basis adaptivity, and anisotropic degree basis adaptivity. We made a distinction between four cases, namely low- and high-dimensional models as well as small and large experimental design sizes.

The comparative study revealed that it is important to carefully select the strategy, since the difference in generalization error between different combinations of methods can be large, even more than an order of magnitude. No single solver or basis-adaptive scheme significantly outperformed the others. However, by dividing the analysis into classes (low-/high-dimensional models and small/large ED sizes), some significant patterns can be identified. The combinations that performed well always involved the solvers $SP_{LOO}$, $SP_{k=5}$, or BCS, but never SPGL1 or OMP. For the basis-adaptive schemes, the picture is less clear, except that combinations involving forward neighbor basis adaptivity achieved the best accuracy most often, and that the static basis (i.e., no basis adaptivity) was in nearly all cases outperformed by basis-adaptive schemes. This might imply that basis-adaptive schemes provide the opportunity to improve the solution by offering more regressors to choose from, but that none of them consistently generates the most suitable basis. The choice of the final basis is guided by the cross-validation error, which is computed based on the solution returned by the respective sparse solver. A good correlation of cross-validation and validation error is crucial for the scheme to make good choices – which for example is not the case for the solver OMP. Overall, there is no combination of solver and basis-adaptive scheme that was consistently best across all models and experimental design sizes, although averaging over all four cases, $SP_{LOO}$ together with forward neighbor basis adaptivity outperformed all other combinations.

Since no clear best combination of solver and basis adaptivity scheme emerged, we introduced an automatic selection step. Based on a suitable error estimate, automatic selection chooses one PCE solution out of several candidate solutions computed by various combinations of solvers and basis adaptivity schemes. We found that automatic selection using any hybrid cross-validation error estimate performs better than the fixed rules obtained from the basis adaptivity benchmark: automatic selection attains the best possible relative MSE almost as often as the best solver-basis adaptivity combination, while being significantly more robust. An additional advantage of automatic selection is that it is independent of model dimension or size of the available experimental design, unlike the proposed fixed rules, which rely on the somewhat arbitrary, albeit simple, classification we applied (low/high dimension and small/large experimental design).

These findings demonstrate that when building a sparse PCE for an expensive black-box computational model, it is worth it to carefully select a sparse solver, and to apply a basis-adaptive
scheme, because the difference in relative MSE between different combinations of methods on
the same experimental design can be larger than one order of magnitude. While we could iden-
tify a number of methods that generally perform well, and others that should be avoided, as we
described above, a superior strategy is to compute several PCE solutions and perform a final
model selection using one of the presented hybrid cross-validation error estimators.

Further research could investigate the use of true cross-validation for automatic selection instead
of the hybrid estimates which we used here. Also, it might be possible to identify other problem
characteristics besides model dimension and the size of the available experimental design to
guide the choice of methods in the sparse PCE framework. A promising class of methods
combines basis adaptivity with the sequential enrichment of the experimental design, adapted
to the current candidate or active basis (Fajraoui et al., 2017; Diaz et al., 2018; Hampton and
Doostan, 2018), which might be able to further improve on the results obtained here.

Acknowledgements

We thank John Jakeman (Sandia National Laboratories), Negin Alemazkoor (University of Virginia), Hadi Meidani (University of Illinois at Urbana-Champaign), and Jerrad Hampton (University of Colorado Boulder) for generously providing their code and explanations. We thank John Jakeman for his useful hints to relevant literature on basis adaptivity.

References

Abbiati, G., Marelli, S., Tsokanas, N., Sudret, B., and Stojadinovic, B. (2021). A global sensitivity analysis framework for hybrid simulation. *Mech. Syst. Signal Pr.*, 146(106997).

Adams, B. M., Bauman, L. E., Bohnhoff, W. J., Dalbey, K. R., Ebeida, M. S., Eddy, J. P., Eldred, M. S., Hough, P. D., Hu, K. T., Jakeman, J. D., Stephens, J. A., Swiler, L. P., Vigil, D. M., and Wildey, T. M. (2014). *Dakota, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.0 User's Manual*. Sandia National Laboratories. Technical Report SAND2014-4633 (Updated November 2015 (Version 6.3)).

Alemazkoor, N. and Meidani, H. (2017). Divide and conquer: An incremental sparsity promoting
compressive sampling approach for polynomial chaos expansions. *Comput. Methods Appl.
Mech. Engrg.*, 318:937–956.

Babacan, S. D., Molina, R., and Katsaggelos, A. K. (2010). Bayesian compressive sensing using Laplace priors. *IEEE Trans. Image Process.*, 19(1):53–63.

Blatman, G. and Sudret, B. (2009). Anisotropic parcimonious polynomial chaos expansions based on the sparsity-of-effects principle. In Furuta, H., Frangopol, D., and Shinozuka, M., editors, *Proc. 10th Int. Conf. Struct. Safety and Reliability (ICOSSAR’2009), Osaka, Japan.*