Polynomial Chaos Expansions: Fast Surrogates for Uncertainty Propagation

Polynomial chaos expansion is a surrogate modeling technique that replaces an expensive simulation with a series of orthogonal polynomials whose coefficients can be used to estimate uncertainty directly. Once fitted, a PCE surrogate gives you analytical access to output means, variances, confidence intervals, and Sobol sensitivity indices without running a single additional simulation.

This article covers PCE from the theoretical foundation through practical implementation. You’ll see exactly how polynomial bases are selected, how coefficients are estimated by regression, quadrature, or projection, how sparse techniques extend PCE to moderate dimensions, and which failure modes break standard PCE. The included Python examples use Chaospy, UQpy, and scikit-learn to show complete, validated workflows.

Key Takeaways

• PCE replaces random sampling with deterministic approximation: Instead of drawing thousands of Monte Carlo samples, PCE fits a polynomial to the model response. The polynomial coefficients then provide analytical moments and sensitivity indices.
• The polynomial basis must match the input distribution: Legendre for uniform, Hermite for Gaussian, Laguerre for exponential. A mismatched basis destroys orthogonality and invalidates all coefficient estimates.
• Non-intrusive PCE is the default choice for most workflows: You keep your simulator as a black box, evaluate it at training points, and estimate coefficients via regression. This avoids modifying solver equations.
• Sparse PCE methods handle ~50 parameters: Hyperbolic truncation, LASSO, and Smolyak grids reduce the polynomial basis size dramatically. They work when the model response has dominant low-order interactions, but the sparsity assumption must be verified.
• PCE is not universal: Discontinuous responses, threshold events, contact changes, and solver failures produce non-physical oscillations (the Gibbs phenomenon). When the response is discontinuous, use Multi-Element PCE, rational polynomial chaos, or a hybrid surrogate.

How Polynomial Chaos Expansion Works

Polynomial chaos expansion represents a model output as a weighted sum of orthogonal polynomials. Consider a scalar output Y that depends on an input vector X:

Y = M(X)

The PCE approximation replaces M(X) with a polynomial:

Y(X) ≈ Σ cₖ Ψₖ(X)

Here, Ψₖ are multivariate orthogonal polynomials and cₖ are deterministic coefficients. The set of polynomials must be orthogonal with respect to the input probability distribution, and the coefficients are determined by the chosen estimation method.

Orthogonal Polynomials and Distribution Matching

The polynomial family is tied directly to the input distribution. This is not a design choice — it is a mathematical requirement for orthogonality.

The orthonormality condition states:

∫ Ψᵢ(ξ) Ψⱼ(ξ) w(ξ) dξ = δᵢⱼ

where w(ξ) is the probability density of the input distribution and δᵢⱼ is the Kronecker delta. When this condition holds, the PCE coefficients relate directly to output statistics through Parseval’s identity.

Multivariate Polynomial Construction

When inputs are independent, the multivariate polynomial basis is a tensor product of univariate polynomials:

Ψₖ(X) = Π ψₖᵢ(Xᵢ)

For a total-degree basis of order p in d uncertain parameters, the number of terms is:

N = (d + p)! / (d! · p!)

For 20 parameters and order three, this produces 1,771 terms. A stable regression typically requires more training points than coefficients, so the computational demand grows quickly.

Non-Intrusive vs. Intrusive PCE

There are two fundamentally different approaches to constructing a PCE. The non-intrusive approach treats the simulator as a black box. The intrusive approach reformulates the governing equations.

Non-Intrusive PCE

Non-intrusive PCE evaluates the deterministic simulator at selected points in uncertainty space and fits a polynomial regression. The workflow:

1. Select a polynomial basis aligned with the input distributions.
2. Generate training points (quadrature or sampling).
3. Run the simulator at those points.
4. Estimate expansion coefficients.
5. Validate on independent test points.

This approach adds no solver modifications. You can wrap an existing CFD, FEM, phase-field, or molecular model and get a PCE surrogate without changing any internal code.

Intrusive PCE (Stochastic Galerkin)

Intrusive PCE substitutes the polynomial expansion directly into the governing equations and derives a coupled deterministic system for the coefficients. For a simple model:

∂u/∂t = L(u) + f → Σ ∂cₖΨₖ/∂t = L(Σ cₖΨₖ) + f

After projecting onto each basis function, the coupled system solves for all coefficients simultaneously. This can be more efficient because it uses the full structure of the equations, but it requires modifying the solver to handle the coupled system. For most production workflows, the intrusive approach is impractical — the non-intrusive regression approach is preferred.

How to Estimate PCE Coefficients

The coefficient estimation method determines both accuracy and computational cost. There are three principal approaches.

Regression-Based Estimation

The non-intrusive regression approach is the most widely used. It samples input points, evaluates the model, and fits the polynomial by linear regression:

y ≈ Φ c

where Φ is the polynomial evaluation matrix and y is the vector of model evaluations. If N_samples ≥ N_terms, the least-squares solution is:

c = (ΦᵀΦ)⁻¹ Φᵀ y

This is straightforward to implement with scikit-learn’s PolynomialFeatures:

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import Ridge

# Training data: N_samples × N_inputs
X_train = np.random.uniform(-1, 1, (100, 3))

# Polynomial features up to degree 2
poly = PolynomialFeatures(degree=2, include_bias=False)
Phi = poly.fit_transform(X_train)  # 100 × 10 matrix

# Ridge regression (regularization stabilizes coefficient estimation)
model = Ridge(alpha=1e-4)
c = model.fit(Phi, y_train).coef_

Regression with regularization (Ridge, LASSO) is numerically stable and works well for moderate dimensions. The downside is that the regression matrix becomes ill-conditioned when the training set is small relative to the basis size.

Quadrature-Based (Gauss-Formula) Estimation

When the coefficient index k is one-dimensional (a single polynomial term), the coefficient cₖ can be computed by numerical quadrature:

cₖ = ∫ M(X) Ψₖ(X) w(X) dX ≈ Σᵢ M(Xᵢ) Ψₖ(Xᵢ) w(Xᵢ) / ||Ψₖ||²

For multivariate PCE, sparse quadrature grids are needed. This approach is exact when the model is a polynomial of degree ≤ 2p, but it becomes impractical in high dimensions because the quadrature grid size grows exponentially.

Projection-Based (Stochastic Collocation)

Stochastic collocation constructs the PCE by interpolating the model at selected collocation points. The coefficient for term k is:

cₖ = Σᵢ M(Xᵢ) Ψₖ(Xᵢ) / Ψₖ(Xᵢ)²  (diagonal weight)

Collocation avoids the integration step of quadrature and is computationally efficient when sparse grids or adaptive refinement strategies are used. It is the preferred approach when the simulator is cheap enough to evaluate at all quadrature points.

Sparse PCE Methods

The full polynomial basis grows combinatorially with dimension. Sparse PCE methods reduce the basis size while retaining accuracy when the model response is dominated by low-order interactions.

Hyperbolic Truncation

The hyperbolic basis restriction (HBR) with parameter q < 1 limits which multi-indices are included:

∑ αᵢ^(1/q) ≤ p

When q → 1, the HBR approaches the total-degree basis. When q → 0, it approaches the generalized Pareto basis. The parameter q controls the trade-off between basis size and accuracy. For most engineering models, q = 0.5 to 0.7 provides good compression while retaining important interaction terms.

The hyperbolic truncation reduces a 20-parameter, order-3 basis from 1,771 terms to roughly 200-400 terms depending on q.

LASSO and Regularized Regression

LASSO (Least Absolute Shrinkage and Selection Operator) performs regression with an L1 penalty:

min (||Φc||² + λ ||c||₁)

The L1 penalty drives many coefficients exactly to zero, producing a sparse solution. Elastic Net combines L1 and L2 penalties to handle correlated predictors:

min (||Φc||² + λ₁ ||c||₁ + λ₂ ||c||₂²)

LASSO-based sparse PCE is implemented in UQpy as HyperbolicTruncLARS (hyperbolic truncation with Least-Angle Regression). The workflow:

from uqpy.surrogates import HyperbolicTruncLARS
from uqpy.distributions import MultivariateGaussian, Uniform

# Define distributions
dist = MultivariateGaussian(np.zeros(5), np.eye(5))

# Generate samples and run model
samples = dist.sample(n=500)
y_train = model_function(samples)

# Build sparse PCE with q=0.5
pce = HyperbolicTruncLARS(distribution=dist, q=0.5, p=3)
pce.fit(samples, y_train)

# Validate
validation_samples = dist.sample(n=100)
y_pred = pce.predict(validation_samples)
rmse = np.sqrt(np.mean((y_pred - y_test)**2))

Smolyak Sparse Grids

The Smolyak construction combines univariate quadrature rules using the sparse grid formula:

Qₘ = Σ_{|k|=m} (-1)^(m-|k|) · Πᵢ Q_{kᵢ}

where Q_d is the d-dimensional quadrature operator and |k| = Σ kᵢ. Smolyak grids reduce the number of multidimensional points from exponential to polynomial growth in the dimension. They are most effective when the response is smooth and low-order interactions dominate.

When PCE Is Effective

PCE converges exponentially fast for smooth model responses. The convergence rate depends on:

• Smoothness: Analytic responses converge geometrically. Piecewise smooth responses converge algebraically.
• Effective dimension: Even with 50 nominal parameters, if only 5 control the response, PCE converges on 5 dimensions.
• Parameter interactions: Strong higher-order interactions increase the required basis size.
• Experimental design quality: Sobol sequences or sparse quadrature grids perform better than uniform random sampling.

Baseline Context

• PCE provides analytical moments and Sobol indices once fitted, with no additional model evaluations.
• The regression approach treats the simulator as a black box.
• Sparse methods can extend PCE to ~50 uncertain parameters when the model exhibits sparse interactions.
• The Python UQ ecosystem includes Chaospy, UQpy, SALib, SMT, and Uncertainpy.

PCE vs Kriging vs Neural Networks

PCE is one surrogate among many. Understanding the trade-offs helps you choose the right method for your problem.

PCE is the only surrogate that gives exact analytical Sobol indices from the coefficients. Kriging gives prediction intervals but requires numerical differentiation for sensitivity. Neural networks require perturbation-based sensitivity. Reduced-order models are deterministic projections.

When to Choose PCE

PCE is the best choice when:

1. You need analytical moments and sensitivity indices: The coefficients give you E[Y], Var(Y), and Sobol indices without additional computation.
2. The response is smooth and low-dimensional: Analytic responses converge geometrically in PCE.
3. You need fast Monte Carlo on the surrogate: Once fitted, evaluating the polynomial is O(1). This is invaluable for reliability analysis or threshold probability estimation.
4. The simulator is moderately expensive: If each evaluation takes minutes to hours, a 100-500 point PCE is feasible.

When to Avoid PCE

Avoid PCE when:

1. The response is discontinuous: The Gibbs phenomenon produces oscillations near discontinuities. Use Multi-Element PCE or a hybrid surrogate.
2. You need a high-dimensional model: Even sparse PCE becomes unstable beyond ~50 parameters without verified sparsity.
3. The response is highly nonlinear with strong threshold behavior: PCE smooths thresholds. Use Monte Carlo or randomized QMC instead.
4. Data is extremely scarce: With fewer than 20 training points, PCE is unreliable. Kriging or polynomial response surfaces may be more robust.

The Gibbs Phenomenon

The Gibbs phenomenon is the most fundamental failure mode of standard PCE. For discontinuous responses, threshold events, contact changes, or solver failures, the polynomial approximation produces non-physical oscillations near the discontinuity. This is mathematically identical to the oscillations seen when approximating a step function with a Fourier series.

For a smooth response, PCE converges rapidly. For a response with a single jump discontinuity, the polynomial oscillates with amplitude proportional to the jump height and decaying only as 1/√N. Increasing the polynomial order makes the oscillation worse, not better.

Mitigation Strategies

Multi-Element PCE (ME-PCE): The input domain is partitioned into subregions. A separate polynomial is fitted within each subregion, with continuity enforced at boundaries. This is the most reliable approach for discontinuous responses.

Polynomial Chaos-Kriging Hybrid: A PCE captures the smooth component of the response, while a Kriging model captures the residual discontinuity. This leverages the strengths of both methods.

Rational Polynomial Chaos: Padé-Legendre approximations use a ratio of polynomials instead of a single polynomial. The rational form can represent discontinuous jumps more accurately than a pure polynomial.

Active Learning Domain Decomposition: Novák et al. (2023) introduced a sequential approach where the input domain is iteratively decomposed based on the model’s behavior in each region. This is the most flexible approach for complex discontinuous responses.

Python Implementation Guide

Chaospy

Chaospy is a focused library for PCE and probability distributions. It does not aim to be an all-in-one UQ framework.

import chaospy as cp
import numpy as np

# Independent uncertain parameters
distribution = cp.J(
    cp.Uniform(0.8, 1.2),
    cp.Normal(2.0, 0.1)
)

# Total-degree polynomial basis
order = 3
expansion = cp.generate_expansion(order, distribution)

# Experimental design (Sobol)
samples = distribution.sample(100, rule='sobol')

# Model evaluations
def model(sample):
    conductivity = sample[0]
    source = sample[1]
    return source / conductivity + 0.2 * conductivity**2

evaluations = np.array([model(samples[:, i]) for i in range(samples.shape[1])])

# Fit PCE by regression
approximation = cp.fit_regression(expansion, samples.T, evaluations)

# Extract statistics (analytical, no additional model evaluations)
mean = cp.E(approximation, distribution)
variance = cp.Var(approximation, distribution)
first_order = cp.Sens_m(approximation, distribution)
total_order = cp.Sens_t(approximation, distribution)

print(f"Mean: {mean}")
print(f"Standard deviation: {np.sqrt(variance)}")
print(f"First-order Sobol: {first_order}")
print(f"Total-order Sobol: {total_order}")

UQpy

UQpy provides sparse PCE, LARS, and extensive validation utilities.

from uqpy.distributions import MultivariateGaussian
from uqpy.surrogates import HyperbolicTruncLARS, MonteCarlo
import numpy as np

# Define distribution
dist = MultivariateGaussian(np.zeros(5), np.eye(5))

# Generate samples and run model
samples = dist.sample(n=500)
y_train = np.array([model_function(x) for x in samples])

# Build sparse PCE
pce = HyperbolicTruncLARS(distribution=dist, q=0.5, p=3)
pce.fit(samples, y_train)

# Validate on independent test points
test_samples = dist.sample(n=100)
y_pred = pce.predict(test_samples)
y_test = np.array([model_function(x) for x in test_samples])

# Root mean squared error
rmse = np.sqrt(np.mean((y_pred - y_test)**2))
print(f"RMSE: {rmse}")

# Monte Carlo uncertainty propagation on the surrogate
n_mc = 10000
mc_samples = dist.sample(n_mc)
mc_predictions = pce.predict(mc_samples)

mc_mean = np.mean(mc_predictions)
mc_std = np.std(mc_predictions)
print(f"MC propagation mean: {mc_mean}")
print(f"MC propagation std: {mc_std}")

Scikit-learn Regularized PCE

For projects that already use scikit-learn, Ridge or LASSO regression is straightforward.

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import Lasso

# Generate Sobol samples
N = 300
d = 4
from scipy.stats import qmc
sampler = qmc.Sobol(d)
sample = sampler.random(base2=np.ceil(np.log2(N)).astype(int))

# Model evaluations
y = np.array([model_function(row) for row in sample])

# Polynomial features with hyperbolic truncation (manual basis selection)
poly = PolynomialFeatures(degree=2, include_bias=False)
Phi = poly.fit_transform(sample)

# LASSO regression
lasso = Lasso(alpha=1e-3)
c = lasso.fit(Phi, y).coef_

print(f"Non-zero coefficients: {np.sum(c != 0)} out of {len(c)}")
print(f"Coefficients: {c}")

Validating a PCE Surrogate

A PCE surrogate should never be trusted without validation. The polynomial coefficients provide exact moments, but those moments are meaningless if the surrogate does not approximate the model accurately.

Validation Protocol

1. Generate validation points: Independent from training. Use Sobol sequences or Latin hypercube sampling.
2. Compute RMSE and maximum error: Compare PCE predictions against true model evaluations at validation points.
3. Check moment accuracy: Compare PCE-extracted mean and variance against Monte Carlo estimates on the validation set.
4. Monitor sensitivity stability: If Sobol indices change significantly between training and validation, the PCE is not converged.
5. Check for Gibbs oscillations: Plot PCE predictions against a grid. Non-smooth oscillations indicate discontinuous response failure.

Convergence Monitoring

Track the polynomial order and number of training points. As both increase, the RMSE should decrease monotonically. If RMSE stagnates or increases with higher order, the basis is overfitting the training data.

# Monitoring convergence as polynomial order increases
for order in [2, 3, 4, 5]:
    expansion = cp.generate_expansion(order, distribution)
    samples = distribution.sample(N)
    evaluations = np.array([model(s) for s in samples])
    pce = cp.fit_regression(expansion, samples.T, evaluations)
    
    # Validate
    test = distribution.sample(N_test)
    y_pred = np.array([pce(s) for s in test])
    y_true = np.array([model(s) for s in test])
    
    rmse = np.sqrt(np.mean((y_pred - y_true)**2))
    print(f"Order {order}: RMSE = {rmse:.6f}")

Common Mistakes and How to Avoid Them

1. Mismatched Polynomial Basis

Using Legendre polynomials for Gaussian inputs destroys orthogonality. The coefficients become biased and the moment extraction is invalid. Always match the polynomial family to the input distribution.

2. Overfitting the Polynomial Basis

A degree-5 basis in 10 parameters has 2,002 terms. With 2,000 training points, the regression is barely stable. Use information criteria (AIC, BIC) or cross-validation to select the polynomial order.

3. Trusting PCE Moments Without Validation

The mean and variance extracted from PCE coefficients are exact mathematical expressions. But they are only valid if the polynomial approximation is accurate. Always validate against independent test points before reporting moments.

4. Ignoring the Gibbs Phenomenon

For discontinuous responses, the polynomial oscillates near discontinuities. Report the RMSE on a grid that includes discontinuity regions. If RMSE is dominated by a few oscillatory points, the PCE is not valid.

5. Assuming Sparsity Without Verifying

LASSO or hyperbolic truncation assumes the model response is dominated by low-order interactions. If the true model has strong high-order interactions, the sparse PCE will miss them. Verify sparsity by comparing the coefficient norms of full vs. sparse bases.

Summary and Next Steps

Polynomial chaos expansion converts a deterministic simulator into a stochastic surrogate that provides analytical access to output statistics, uncertainty intervals, and sensitivity indices. The method excels when the response is smooth, the dimensionality is moderate, and computational cost is high enough to justify surrogate construction.

Recommended Workflow

1. Verify numerical convergence of the deterministic simulator before UQ. Discretization error and numerical uncertainty should be established before propagation.
2. Assign distributions from experimental data or documented expert judgment. Do not default to uniform distributions.
3. Select the polynomial basis by matching distributions to orthogonal polynomials. Verify orthonormality conditions.
4. Choose a coefficient estimation method: regression for non-intrusive PCE, quadrature for efficient low-dimensional cases, collocation for cheap models.
5. Validate the surrogate against independent test points and monitor convergence with increasing polynomial order.
6. Report moments and sensitivity indices only after validation. Include confidence intervals for estimated statistics.

What to Avoid

• Using PCE for discontinuous responses without mitigation
• Assuming sparse PCE works for every high-dimensional problem
• Reporting PCE moments without validation
• Confusing non-intrusive PCE with intrusive stochastic Galerkin
• Using full polynomial bases in problems where the basis size exceeds the training set

Related Guides

• Uncertainty Quantification and Sensitivity Analysis in Scientific Simulation — Broader UQ workflow, Sobol sensitivity analysis, Monte Carlo methods, and Python UQ tools.
• Monte Carlo Methods for Scientific Simulations — Compare random sampling, QMC, LHS, and convergence behavior against PCE.
• Machine Learning Surrogates for Scientific Simulations — Build and validate neural network and Gaussian process surrogates alongside PCE.
• Verification vs Validation in Scientific Simulations — Establish numerical convergence before starting a UQ campaign.

References

• Caflisch, R. E., & Mosolova, N. (1998). Monte Carlo and Quasi-Monte Carlo methods. UCLA CAM Report.
• Novák et al. (2023). Active Learning Domain Adaptive Localized Surrogates. Computer Methods in Applied Mechanics and Engineering. DOI
• Schneider (2023). Rational Polynomial Chaos Expansion. TUM dissertation. PDF
• Helmholtz UQ Dictionary. Polynomial Chaos Expansion entry. URL
• COMSOL Learning Center. Surrogate Modeling Theory and PCE. URL

Need Help Building PCE Surrogates for Your Simulations?

If you are struggling with polynomial chaos implementation, basis selection, sparse PCE regularization, or surrogate validation, we can help. We specialize in Python-based uncertainty quantification workflows and can guide you from basic PCE scripts to production-grade surrogate pipelines.