PaperlayerView exact revision on arXiv

Source-linked AI summary

Sparse Polynomial Chaos Expansions: Literature Survey and Benchmark

Nora Lüthen, Stefano Marelli, Bruno Sudret

arXiv:2002.01290v4math.NAstat.COstat.ME

TL;DR

Computer models can involve many input parameters and expensive evaluations, while the choice among numerous sparse PCE methods has lacked clear guidance. The paper surveys and classifies these methods, benchmarks selected approaches, and finds that solver and sampling choices can change relative MSE by several orders of magnitude, with different methods favored across dimensionality and experimental-design regimes.

  • Problem

    Computer models can involve many input parameters and expensive evaluations, and the growing number of sparse PCE methods makes choosing an approach difficult.

  • Method

    The paper presents a literature survey, develops a framework for describing sparse PCE computation, and benchmarks selected methods on engineering models varying in dimensionality and complexity.

  • Results

    Solver and sampling-scheme choices can change relative MSE by several orders of magnitude; BCS performs best for high-dimensional models, while SPLOO with coherence-optimal sampling suits low-dimensional models with large experimental-design sizes.

  • Takeaways & Limitations

    Sparse PCE methods are not uniformly superior: practical performance depends on model dimensionality and experimental-design size.

  • Takeaways & Limitations

    The benchmark does not support a universally best method because no selection is optimal across all engineering models, and large validation sets may be unavailable under computational constraints.

Abstract

from arXiv · show

Sparse polynomial chaos expansions (PCE) are a popular surrogate modelling method that takes advantage of the properties of PCE, the sparsity-of-effects principle, and powerful sparse regression solvers to approximate computer models with many input parameters, relying on only few model evaluations. Within the last decade, a large number of algorithms for the computation of sparse PCE have been published in the applied math and engineering literature. We present an extensive review of the existing methods and develop a framework for classifying the algorithms. Furthermore, we conduct a unique benchmark on a selection of methods to identify which approaches work best in practical applications. Comparing their accuracy on several benchmark models of varying dimensionality and complexity, we find that the choice of sparse regression solver and sampling scheme for the computation of a sparse PCE surrogate can make a significant difference, of up to several orders of magnitude in the resulting mean-squared error. Different methods seem to be superior in different regimes of model dimensionality and experimental design size.

1 Introduction

Sparse PCE offers a surrogate-modelling approach for expensive, high-dimensional computer models when only a few evaluations are affordable. This paper surveys the rapidly growing sparse-PCE literature, organizes its methods in a common framework, and benchmarks selected approaches.

  • Motivation: Expensive computer models with many input parameters can make the model evaluations required for uncertainty quantification infeasible.Surrogate models address this constraint by approximating the original model while requiring relatively few evaluations.
  • Sparse polynomial chaos expansions: Sparse PCE represents computational models with polynomial bases while setting most coefficients to zero.The sparsity-of-effects principle and compressibility motivate this representation.
  • Sparse polynomial chaos expansions: Nonintrusive sparse PCE treats the computational model as a black box and can support uncertainty propagation and sensitivity analysis.Moments and Sobol’ sensitivity indices can be computed analytically from PCE.
  • Research gap: The literature contains many promising sparse-PCE algorithms but often lacks thorough comparisons with previously published methods.This limits clear guidance on which approaches perform best.
  • Contributions: The paper develops a general classification framework and benchmarks selected methods on representative realistic models.The benchmark assesses which methods perform best across the selected models.

2 Framework and literature survey for sparse polynomial chaos expansions

The framework decomposes sparse-PCE computation into basis selection, experimental design, sparse regression, and model selection, while organizing the literature around these choices. It highlights why sparse regression can work with fewer evaluations than ordinary least squares and shows that method performance depends on the problem regime.

  • Basis and coefficients: A truncated PCE selects a finite basis from an orthonormal polynomial expansion, with accuracy depending on the chosen multi-index set and coefficients.The input components are assumed independent, and basis functions are products of univariate orthonormal polynomials.
  • Regression: Ordinary least squares uses an N × P regression matrix of basis-function evaluations, with a heuristic requirement of N ≈ 2P–3P evaluations.This requirement can be infeasible for high-dimensional or high-degree models.
  • Regression: Sparse regression imposes coefficient sparsity, enabling recovery or approximation with fewer design points than basis functions.Sparsity also contributes to robustness to noise and makes sparse PCE preferable when OLS-based evaluation costs are infeasible.
  • Framework: The proposed framework classifies available choices and enables structured comparison of their impact on sparse-PCE performance.It also reveals combinations of methods not previously considered in the literature.
  • Framework: The computation pipeline chooses candidate polynomials and an experimental design, obtains model evaluations, applies a sparse solver, and evaluates a model-selection criterion.If the solution is unsatisfactory, the basis or experimental design may be adapted before repeating the process.
  • Scope: The benchmark does not consider iterative basis adaptation or experimental-design enrichment and is not exhaustive in its experimental-design review.These boundaries are left for future benchmarks, while the review acknowledges many additional design approaches.
  • Experimental design: Sampling methods show regime dependence: coherence-optimal sampling is favored for low dimension and high degree, whereas Monte Carlo often performs better for high dimension and low degree.A sequential design strategy generally produces smaller validation error than the compared alternative.

3 Numerical results

The benchmark finds that solver performance depends strongly on model dimensionality and experimental-design size, with substantial differences in generalization error. BCS is strongest for small designs and high-dimensional models, while SPLOO excels for low-dimensional models with large designs.

  • The benchmark addresses a gap in comprehensive comparisons of sparse PCE methods, whose published evaluations were usually limited.
  • Differences between solver generalization errors can reach several orders of magnitude, showing that solvers exploit available information differently.
  • Low-dimensional models, small ED sizes: 31% of runs identified BCS as best for low-dimensional models with small experimental designs, while SPLOO was within two times the smallest error in 72%.SPLOO was within one order of magnitude in 94% of runs, compared with 87% for BCS.
  • Low-dimensional models, large ED sizes: SPLOO was best for low-dimensional models with large experimental designs, achieving the smallest relative MSE in 42% of runs and staying within two times the smallest in 73%.It was within one order of magnitude of the smallest relative MSE in 99% of runs.
  • High-dimensional models: 45% (53%) of high-dimensional runs made BCS the best solver for small (large) designs, and it was within two times the smallest relative MSE in 96% (97%).Solver choice had less impact for high-dimensional than low-dimensional models.
  • Sampling schemes: Sampling effects were largest for low-dimensional models: coh-opt and D-opt(coh-opt) helped with large designs but were mixed for small designs and no better than LHS for high-dimensional models.MC and LHS performed comparably, with LHS usually slightly better for low-dimensional models.

4 Discussion and conclusions

The paper surveys sparse PCE methods, organizes them in a general computation framework, and benchmarks solver–sampling combinations across models and experimental-design sizes. Results show that method choice can substantially affect surrogate quality, with preferred combinations depending on dimensionality and design size.

  • The paper presents a literature survey and framework for describing the general computation procedure for sparse PCE.
  • The benchmark uses 11 example functions representing engineering problems with different dimensionality and complexity, testing sparse solvers and sampling schemes on fixed basis functions.
  • Up to several orders of magnitude separate solver–sampling choices in relative MSE, with rankings mostly independent of one another and effects larger for low-dimensional models.
  • For low-dimensional models with small designs, BCS performs best most often, while SPLOO is especially robust and dominates with large designs.
  • For low-dimensional models with feasible basis sizes, near-optimal sampling outperforms other sampling schemes regardless of solver.
  • For high-dimensional models, BCS is by far the best solver, and LHS outperforms advanced sampling schemes across solvers.
  • The conclusions are limited because the benchmark models cannot cover the full space of engineering models and basis adaptivity was not investigated.

A.1 Sampling according to the input distribution

Input-distribution sampling methods select experimental-design points directly from, or systematically across, the input space. The section covers Monte Carlo, Latin hypercube, stratified, and weighted approaches.

  • Latin hypercube sampling: Latin hypercube sampling divides each marginal quantile space into N intervals and randomly combines one sample from each interval into input vectors.The vectors are transformed into the input space using an isoprobabilistic transform.
  • Stratified and partially stratified sampling: LHS can reduce the variance of linear regression estimates when main effects dominate, while partially stratified methods target dominant interaction terms.LPSS combines LHS and stratified sampling and is especially beneficial when prior knowledge identifies interacting variable groups, including problems with d = 100.
  • Stratified and partially stratified sampling: Stratified sampling divides the input space into disjoint strata and weights samples by each stratum’s probability mass, reducing statistical-estimator variance.Partially stratified methods extend this idea toward interaction-dominated settings.
  • Monte Carlo sampling: Monte Carlo sampling draws directly from the input distribution and is a special case of coherence-based sampling theory.Coherence bounds yield associated sample requirements for sparse recovery.
  • Weighted sampling: Weighted sampling modifies both the regression system and sampling distribution to improve sparse-recovery properties.The weighted matrix can be interpreted as evaluating a modified basis, with design points drawn from a corresponding modified distribution to preserve orthonormality.

A.2.1 Coherence, isotropy, and weighted orthonormal systems

Coherence and isotropy govern sparse-recovery guarantees for PCE regression matrices. Weighted basis systems and distributions can lower coherence while retaining orthonormality, although some bases and distributions create scope limitations.

  • Isotropy: PCE regression matrices are isotropic when experimental-design points are sampled from the input distribution.The induced row distribution arises by propagating the input distribution through the basis functions.
  • Recovery concepts: For an isotropic regression matrix, noiseless perfect recovery of an s-sparse solution requires a sample count proportional to µ(A, {ψα})s log(P) with high probability.Here s is solution sparsity and P = |A| is the number of basis functions; analogous results hold in noise.
  • Coherence: Lower coherence in an orthonormal system requires fewer samples for perfect recovery, motivating coherence-optimal weighted systems.The target system is orthonormal under a modified distribution and has lower coherence than the original system.
  • Weighted orthonormal systems: Weighting with w(x) modifies the regression system into a basis {˜ψα(x) = w(x)ψα(x)} and requires sampling from a suitably modified input distribution.The resulting construction preserves orthonormality with respect to the modified distribution.
  • Scope conditions: If G = B, the weighted system achieves minimal coherence, but some polynomial bases lack a finite upper bound needed by this construction.For Gaussian distributions, the relevant quantity can be infinite; restricting the domain can recover an approximately isotropic construction.

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

Asymptotic sampling uses loose polynomial upper bounds to construct tractable modified distributions. Its coherence advantage depends on the relationship between dimension and polynomial degree.

  • Construction: Asymptotic sampling uses a loose upper bound G(x) ≥ maxα∈A |ψα(x)| to construct improved-coherence sampling distributions for Legendre and Hermite bases.Related asymptotic distributions have also been implemented for Laguerre and Jacobi polynomial families.
  • Hermite polynomials: For Hermite polynomials with standard Gaussian variables, the construction restricts sampling to a bounded subset and yields a uniform distribution there.The associated support is selected using a degree-dependent radius.
  • Legendre polynomials: For Legendre polynomials, asymptotic sampling has smaller coherence when d < p, whereas standard sampling has smaller coherence when d > p.The latter regime is described as more common in applications.
  • Recovery behavior: The sampling scheme with smaller coherence exhibits better recovery rates, and this theoretical pattern is confirmed numerically.The comparison is stated for standard versus asymptotic sampling in the Legendre setting.

A.2.3 Coherence-optimal sampling

Coherence-optimal and candidate-based designs improve regression-matrix properties through weighting or point selection. Their theoretical and practical advantages depend on basis, dimension, candidate sets, and design size.

  • Coherence-optimal sampling: Coherence-optimal sampling chooses a modified distribution and weights so the transformed basis is orthonormal and, when G = B, has minimum possible coherence.The benchmark implements rejection-based sampling; MCMC alternatives can be very slow in high dimensions.
  • Coherence-optimal sampling: Numerical examples report better recovery and smaller errors in various norms for coherence-optimal sampling than for standard and asymptotic sampling.This comparison is attributed to prior numerical studies.
  • Coherence-optimal sampling: For Gaussian inputs, coherence-optimal and asymptotic sampling have a larger spread than input sampling, potentially reducing simulation accuracy far from typical operating conditions.Their support is described as a ball, so the practical concern is especially relevant to engineering simulations.
  • Candidate-based designs: Candidate-based designs select points from a pool to improve regression-matrix criteria, but greedy methods generally find only local optima and depend on candidate-set quality.Candidate sets in the literature include MC, LHS, coherence-optimal, and Christoffel sparse-approximation samples.
  • D-optimal sampling: RRQR-based D-optimal sampling remains usable when N < P, where conventional D-value criteria yield D(Ψ) = 0.RRQR is not guaranteed to find the maximal D-value and may return only a local optimum.
  • S-optimal sampling: S-value selection heuristically improves column orthogonality and the determinant of the information matrix, but S(Ψ) = 0 when N < P.The criterion is defined under a mutual-orthogonality assumption for candidate-matrix columns.

A.4 Illustration of sampling schemes

The appendix illustrates experimental designs for uniform and standard Gaussian inputs and introduces sparse regression solvers used in the benchmark. It describes greedy stepwise methods, including OMP, and their sparsity-selection and stopping strategies.

  • Sampling illustrations: For d = 2, degree p = 12, and N = 100, the illustrated designs use candidate sets of size M = 1000.The examples cover both uniform and standard Gaussian inputs.
  • Sampling illustrations: Standard Gaussian sampling schemes can spread points up to seven standard deviations from the mean at degree p = 12.The asymptotic, coherence-optimal, and matrix-optimal distributions have spreads that grow with basis degree.
  • Sampling illustrations: Engineering models may be less accurate or fail outside the input region where their calibration distribution has nonnegligible mass.This bounds the practical usefulness of designs with substantial probability spread in low-mass regions.
  • Sparse solvers: Sparse regression commonly replaces infeasible ℓ0-minimization with convex ℓ1 formulations such as BPDN and LASSO.Under certain conditions, ℓ1- and ℓ0-minimization solutions coincide, while convex formulations have considerably smaller numerical cost.
  • Sparse solvers: OMP adds the regressor most correlated with the residual, refits active coefficients by least squares, and iterates until a stopping rule is reached.Stopping can use a known sparsity K, a cross-validated residual threshold, or an early-stop criterion based on LOO error.
  • Sparse solvers: OMP has complexity O(mNP), but may overfit, produce unreliable LOO estimates, and fail to return a sparse solution without an explicit stopping rule.The active set can grow to min{N, P} basis functions unless iterations are stopped earlier.

B.3 Subspace pursuit (SP)

Subspace pursuit is a fixed-sparsity greedy solver that alternates between augmenting and pruning an active basis. The section also describes related sparse-regression formulations and SPGL1’s benchmark behavior.

  • Subspace pursuit: SP maintains an active basis of size K, adds K residual-correlated regressors per iteration, refits by OLS, and removes the K smallest coefficients.Iterations continue until convergence, producing a K-sparse approximation.
  • Subspace pursuit: SP requires 2K ≤ min{N, P} for basis augmentation and OLS regression to remain feasible.Its complexity is O(N(P+K^2)K), reducing to O(NPK) for very sparse vectors with K^2 ∈ O(P).
  • Subspace pursuit: SP performs O(K) iterations in general and O(log K) in certain cases.When K is unknown, cross-validation can select it over a range of candidate sparsity levels.
  • Subspace pursuit: The paper proposes leave-one-out cross-validation instead of 4-fold cross-validation for an SP variant.The proposed variant selects among runs with different K values using leave-one-out cross-validation.
  • SPGL1: In benchmarks on compressible models, SPGL1 was among the slowest solvers and often returned rather dense solutions.This is an empirical benchmark observation rather than a general theoretical property of SPGL1.

B.5 Sparse Bayesian learning

Sparse Bayesian learning embeds PCE regression in a probabilistic model with hyperparameterized coefficient variances and sparsity-inducing priors. Its assumptions can mismatch PCE errors, but the framework can still produce sparse solutions.

  • Assumptions and limitations: BCS generally assumes Gaussian measurement noise, but omitted PCE terms can produce correlated, heteroscedastic, non-Gaussian, and nonzero-mean discrepancies.Thus the probabilistic noise assumption is not generally valid for PCE, although the framework can remain useful for finding sparse solutions.
  • General framework: Bayesian compressive sensing models coefficients and noise probabilistically, with hyperparameters governing coefficient variances and sparsity.Different BCS methods vary in their hyperparameter priors and approximate posterior-computation algorithms.
  • Hierarchical BCS: The hierarchical Babacan et al. formulation adds hyperparameters for coefficient variances and uses exponential and Gamma priors.Its algorithm iteratively optimizes λ, β, and γ_i while holding the other parameters fixed.
  • General framework: Sparsity arises from hyperparameter priors whose effective coefficient distributions peak sharply at zero while retaining heavy tails.This combination encourages zero coefficients while allowing large coefficients.
  • Computation: Sparse Bayesian methods use approximate algorithms because sparsity-encouraging priors are generally intractable for direct analytical computation.The resulting procedures estimate posterior quantities or MAP hyperparameters through iterative approximations.
  • Hierarchical BCS: The MAP hyperparameters induce sparsity because zero variance components force their corresponding PCE coefficients to zero.Gaussian conditional distributions allow subsequent computations to be carried out analytically.
  • Hierarchical BCS: In practice, the hierarchical algorithm does not estimate β well, so the benchmark determines its value by cross-validation.The paper identifies β estimation as crucial despite the algorithm’s poor performance for that parameter.

C Additional results

The appendix reports additional benchmark results complementing Sections 3.4–3.6. Figure 10 presents relative-MSE boxplots for seven additional models.

  • Additional benchmark results: Figure 10 shows boxplots of relative MSE for seven additional benchmark models.The appendix refers readers to Section 3 for the detailed experimental setup.

C.2 Comparison of sampling schemes together with solvers

The benchmark compares solver–sampling-scheme combinations across four models using relative MSE and varying experimental-design sizes. It also examines performance at small and large design sizes.

  • Benchmark setup: Figures 11–14 display boxplots of relative MSE against experimental-design size for Ishigami, borehole, two-dimensional diffusion, and 100D models.The benchmark includes all solver–sampling combinations, yielding 16–20 combinations per model.
  • Benchmark setup: Solvers are encoded by colors, while sampling schemes use varying shades and line styles.The visual encoding distinguishes solver effects from sampling-scheme effects.
  • Benchmark setup: The figures also slice results at small and large experimental-design sizes to compare solver performance.These slices provide direct comparisons between solvers under different design-size regimes.

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

The smaller-basis comparison extends relative-MSE boxplots across additional models and solver–sampling combinations. It includes separate visualizations for Ishigami, borehole, diffusion, and 100D settings.

  • Smaller candidate basis: Figure 15 compares relative MSE against experimental-design size for Ishigami and borehole using a smaller basis.The remaining solvers shown are LARS, SP, and BCS.
  • Additional models: Figure 10 adds seven models and shows boxplots for six sparse solvers with Latin hypercube sampling.The figure complements the results in Figure 2 and uses thirty replications.
  • Additional models: The damped oscillator and Morris function remain challenging: no solver achieves relative MSE significantly below 0.1, even with large experimental designs.This limitation is reported for the additional-model results.
  • Solver–sampling comparisons: Figures 11–14 compare solver–sampling combinations for Ishigami, borehole, two-dimensional diffusion, and the 100D function.The figures encode solvers by colors and sampling schemes by shades and line styles.
  • Solver–sampling comparisons: The combined sampling schemes are ordered MC–LHS–coh-opt–D-opt(coh-opt) in the per-combination displays.This ordering is used for the relative-MSE comparisons of solvers paired with sampling schemes.

D Benchmark studies

The paper surveys articles and benchmark studies comparing sparse PCE methods and summarizes their reported main results. It also explains how methods, target quantities, and sampling schemes are represented in the overview.

  • Survey overview: Table 4 provides an overview of articles and benchmark studies comparing sparse PCE methods, including their main results.The table is presented as a synthesis of prior comparisons.
  • Comparison criteria: The table’s “best method” is the method delivering the smallest target error reported in each publication.The target may concern moments, Sobol’ indices, or PCE-surrogate generalization error.
  • Notation: The overview defines abbreviations for solvers and sampling schemes and explains parenthetical notation for candidate-set construction.For example, D-opt(coh-opt) denotes D-optimal sampling based on a coherence-optimal candidate set.
  • Notation: A star marks articles that proposed a new method.The notation is part of the table’s explanatory conventions.
Loading 2002.01290v4…