Building a Polynomial Chaos Surrogate

PyApprox Tutorial Library

Approximate an expensive model with a polynomial chaos expansion.

Learning Objectives

After completing this tutorial, you will be able to:

  • Construct a polynomial chaos expansion (PCE) and connect basis polynomials to the input distribution
  • Fit a PCE from training data using least squares
  • Assess surrogate accuracy on an independent test set
  • Visualize how the surrogate captures the model’s input–output relationship

Surrogate models replace expensive simulators with cheap approximations — see The Surrogate Modeling Workflow for an introduction to the four-step surrogate workflow.

Setup: The Beam Benchmark

We use a 2D cantilever beam with \(d = 3\) uncertain parameters (one KLE mode coefficient per subdomain) parameterizing a spatially varying Young’s modulus field, and \(K = 2\) quantities of interest (QoIs): tip deflection and total von Mises stress.

import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.benchmarks.instances import cantilever_beam_2d_linear
from pyapprox.benchmarks.instances.pde.cantilever_beam import MESH_PATHS

bkd = NumpyBkd()

benchmark = cantilever_beam_2d_linear(bkd, num_kle_terms=1, mesh_path=MESH_PATHS[4])
model = benchmark.function()
prior = benchmark.prior()
marginals = prior.marginals()

nvars = model.nvars()
nqoi = model.nqoi()

The Polynomial Chaos Expansion

A PCE approximates the model output as a finite sum of multivariate polynomials:

\[ f(\boldsymbol{\theta}) \approx \hat{f}(\boldsymbol{\theta}) = \sum_{j=0}^{P-1} c_j \, \Psi_j(\boldsymbol{\theta}) \]

where:

  • \(c_j \in \mathbb{R}^K\) are the expansion coefficients (one per QoI)
  • \(\Psi_j(\boldsymbol{\theta})\) are multivariate polynomial basis functions
  • \(P\) is the number of terms, determined by the polynomial degree and the number of input variables

Basis Polynomials and the Input Distribution

The basis functions are products of univariate orthogonal polynomials:

\[ \Psi_j(\boldsymbol{\theta}) = \prod_{i=1}^d \psi_{j_i}^{(i)}(\theta_i) \]

where \(\psi_n^{(i)}\) is the degree-\(n\) polynomial orthogonal with respect to the marginal distribution of \(\theta_i\). The key correspondence is:

Input distribution Orthogonal polynomials
Gaussian (\(\mathcal{N}(0,1)\)) Hermite
Uniform (\(\mathcal{U}(-1,1)\)) Legendre
Beta Jacobi
Gamma Laguerre

For our beam benchmark, the KLE coefficients are standard normal, so the basis is a product of Hermite polynomials. PyApprox selects the correct polynomial family automatically from the marginal distributions via create_pce_from_marginals.

Index Sets

For \(d = 3\) variables and polynomial degree \(p\), the number of multivariate basis functions grows rapidly. A total-degree index set keeps the expansion manageable by excluding high-order interaction terms:

\[ \Lambda_p = \left\{ \mathbf{j} \in \mathbb{N}_0^d : \|\mathbf{j}\|_1 \leq p \right\} \]

from pyapprox.surrogates.affine.indices import (
    compute_hyperbolic_indices,
)

for deg in [1, 2, 3, 4, 5]:
    indices = compute_hyperbolic_indices(nvars, deg, 1.0, bkd)

Collecting Training Data

To fit the PCE, we need training data: pairs \(\{(\boldsymbol{\theta}^{(i)}, \mathbf{q}^{(i)})\}_{i=1}^N\) where \(\mathbf{q}^{(i)} = f(\boldsymbol{\theta}^{(i)})\). The simplest strategy is to draw random samples from the prior.

N_train = 100
np.random.seed(42)
samples_train = prior.rvs(N_train)       # (nvars, N_train)
values_train = model(samples_train)       # (nqoi, N_train)
NoteBetter sampling strategies exist

Random sampling from the prior is the simplest approach, but not the most efficient. The accuracy of a least-squares fit depends on the sample locations relative to the polynomial basis. Induced sampling draws samples from a distribution weighted by the basis functions, producing better-conditioned design matrices and more accurate surrogates for the same number of samples. Choosing the right number of samples and polynomial degree is covered in Overfitting and Cross-Validation.

We also set aside an independent test set to evaluate surrogate accuracy:

N_test = 500
np.random.seed(123)
samples_test = prior.rvs(N_test)
values_test = model(samples_test)     # (nqoi, N_test)

Fitting the PCE

Building the Expansion

We construct a PCE using create_pce_from_marginals, which automatically selects Hermite polynomials for Gaussian marginals and builds the orthonormal basis with a total-degree index set at the specified degree.

from pyapprox.surrogates.affine.expansions.pce import (
    create_pce_from_marginals,
)

degree = 3
pce = create_pce_from_marginals(marginals, degree, bkd, nqoi=nqoi)

nterms = pce.nterms()

Least-Squares Fitting

The coefficients are found by solving a least-squares problem:

\[ \hat{\mathbf{c}} = \arg\min_{\mathbf{c}} \sum_{i=1}^N \left\| \mathbf{q}^{(i)} - \sum_{j=0}^{P-1} c_j\, \Psi_j(\boldsymbol{\theta}^{(i)}) \right\|^2 \]

In matrix form: \(\hat{\mathbf{C}} = (\boldsymbol{\Phi}^\top \boldsymbol{\Phi})^{-1} \boldsymbol{\Phi}^\top \mathbf{V}\), where \(\Phi_{ij} = \Psi_j(\boldsymbol{\theta}^{(i)})\) is the basis matrix and \(\mathbf{V}\) contains the training values.

We use LeastSquaresFitter to solve this:

from pyapprox.surrogates.affine.expansions.fitters.least_squares import (
    LeastSquaresFitter,
)

ls_fitter = LeastSquaresFitter(bkd)
ls_result = ls_fitter.fit(pce, samples_train, values_train)
pce = ls_result.surrogate()
NoteAlternative fitters

Least squares is the simplest fitting method, but others exist. For high-dimensional problems where most coefficients are near zero, Orthogonal Matching Pursuit (OMP) or BPDN (L1-regularized least squares) exploit sparsity by selecting only the most important basis terms. These are covered in Sparse Polynomial Approximation.

Assessing Surrogate Accuracy

We evaluate the fitted PCE on the independent test set and compute the relative \(L_2\) error for each QoI:

\[ \epsilon_k = \frac{\| f_k(\boldsymbol{\Theta}_{\text{test}}) - \hat{f}_k(\boldsymbol{\Theta}_{\text{test}}) \|_2}{\| f_k(\boldsymbol{\Theta}_{\text{test}}) \|_2}, \qquad k = 1, \ldots, K \]

predictions_test = pce(samples_test)     # (nqoi, N_test)

qoi_labels = [r"$\delta_{\mathrm{tip}}$", r"$\sigma_{\mathrm{tot}}$"]

Figure 1 shows the predicted vs. true QoI values on the test set. Points near the diagonal indicate a good fit.

Figure 1: PCE predictions vs. true model values on the independent test set. Left: tip deflection. Right: total von Mises stress. Points near the diagonal indicate good surrogate accuracy.

Exploring the Surrogate

With the expensive model, generating a dense cross-section plot — fixing some inputs and sweeping others — would require hundreds of additional FEM solves. With the fitted PCE, we can produce these plots instantaneously.

PyApprox provides CrossSectionReducer and PairPlotter to automate this. The reducer fixes non-plotted variables at nominal values (here, zero for standard-normal inputs); the pair plotter creates a lower-triangular grid of 1D cross sections (diagonal) and 2D contours (off-diagonal).

from pyapprox.interface.functions.marginalize import CrossSectionReducer
from pyapprox.interface.functions.plot.pair_plot import PairPlotter

nominal = bkd.zeros(nvars)
plot_limits = bkd.full((nvars, 2), 0.0)
plot_limits[:, 0] = -3.0
plot_limits[:, 1] = 3.0

input_labels = [rf"$\theta_{{{i+1}}}$" for i in range(nvars)]

reducer = CrossSectionReducer(pce, nominal, bkd)
plotter = PairPlotter(reducer, plot_limits, bkd, variable_names=input_labels)

fig, axes = plotter.plot(npts_1d=61,
                         contour_kwargs={"qoi": 0, "cmap": "cividis", "levels": 20},
                         line_kwargs={"qoi": 0})
plt.tight_layout()
plt.show()
Figure 2: PCE cross sections for tip deflection. Diagonal: 1D slices through each input. Lower triangle: 2D contour plots for each pair of inputs. Non-plotted variables are fixed at zero.

These plots reveal how each input affects tip deflection and where interactions between inputs are strong (curved contours) or weak (parallel lines). Producing them required zero additional model evaluations — only cheap polynomial evaluations through the surrogate.

Key Takeaways

  • Surrogate models replace expensive simulators with cheap approximations, enabling analyses — like dense cross-section sweeps — that would be infeasible with the original model
  • A polynomial chaos expansion represents the QoI as a polynomial in the uncertain inputs, with basis functions chosen to be orthogonal with respect to the input distribution
  • PyApprox automatically selects the correct polynomial family from the marginal distributions via create_pce_from_marginals
  • A degree-3 PCE fitted with 100 training samples already captures the beam model’s behavior well, as verified by test-set error

Exercises

  1. Increase the polynomial degree from 3 to 5. How does the test-set error change? What happens to the oversampling ratio?

  2. In the cross-section plots, change the nominal values from zero to \(\pm 1\) using CrossSectionReducer. How do the cross sections change?

  3. Reduce the training set to \(N = 30\). How does the scatter plot change? At what oversampling ratio does the fit start to degrade?

Next Steps