Polynomial chaos

Polynomial chaos, also called a polynomial chaos expansion, is a representation of a random quantity as a series of orthogonal polynomials in standardized random variables. Despite its name, the method is not a theory of chaotic dynamics. The word “chaos” derives from the historical description of Gaussian functionals as homogeneous chaoses within probability theory.

Polynomial chaos provides a spectral description of uncertainty. Random dependence is resolved through an orthogonal polynomial basis, while spatial or temporal dependence remains in deterministic coefficients. The resulting representation is used in uncertainty quantification, particularly for mathematical models whose inputs contain specified probability distributions.

Mathematical formulation

Let (Y) be a square-integrable random quantity defined on a probability space ((\Omega,\mathcal F,\mathbb P)). A polynomial chaos expansion expresses (Y) in the form

[ Y(\omega)=\sum_{\boldsymbol{\alpha}\in\mathbb N_0^d} y_{\boldsymbol{\alpha}}, \Psi_{\boldsymbol{\alpha}}!\left(\boldsymbol{\xi}(\omega)\right), ]

where (\boldsymbol{\xi}=(\xi_1,\ldots,\xi_d)) is a finite-dimensional random vector and (\boldsymbol{\alpha}) is a multi-index. The functions (\Psi_{\boldsymbol{\alpha}}) are multivariate polynomials orthogonal with respect to the probability measure of (\boldsymbol{\xi}). The quantities (y_{\boldsymbol{\alpha}}) are deterministic coefficients.

Orthogonality means that

[ \mathbb E!\left[ \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi}) \Psi_{\boldsymbol{\beta}}(\boldsymbol{\xi}) \right]

\gamma_{\boldsymbol{\alpha}} \delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}, ]

where (\gamma_{\boldsymbol{\alpha}}) is the squared norm of the basis polynomial and (\delta_{\boldsymbol{\alpha}\boldsymbol{\beta}}) is the multi-index form of the Kronecker delta. When the basis is normalized, every (\gamma_{\boldsymbol{\alpha}}) equals one.

For an orthogonal basis, the coefficient associated with (\boldsymbol{\alpha}) is

[ y_{\boldsymbol{\alpha}}

\frac{ \mathbb E!\left[ Y\Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi}) \right] }{ \mathbb E!\left[ \Psi_{\boldsymbol{\alpha}}^2(\boldsymbol{\xi}) \right] }. ]

This identity is the orthogonal projection of (Y) onto the polynomial subspace. If the polynomial system is complete in the relevant Hilbert space, the infinite expansion converges to (Y) in the mean-square sense.

For a random field (u(\mathbf x,\omega)), the coefficients depend on the deterministic coordinate:

[ u(\mathbf x,\omega)

\sum_{\boldsymbol{\alpha}} u_{\boldsymbol{\alpha}}(\mathbf x) \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi}(\omega)). ]

The expansion therefore separates stochastic coordinates from physical coordinates without asserting statistical independence between the output coefficients.

Historical development

Norbert Wiener introduced homogeneous polynomial chaoses in 1938 while studying functionals of Brownian motion. His construction used Hermite polynomials in Gaussian random variables and decomposed square-integrable Gaussian functionals into mutually orthogonal subspaces. The term “polynomial chaos” developed from this Gaussian setting rather than from the later mathematical study of deterministic chaos.

In 1946, You Watanabe formulated the finite-dimensional coefficient projection identity for Hermite-chaos truncations and separated the mean-square residual into contributions from omitted polynomial degrees. Watanabe’s formulation treated the truncated expansion as an orthogonal projection rather than as an algebraic approximation, establishing the error interpretation subsequently used in computational implementations.

R. H. Cameron and W. T. Martin proved in 1947 that suitably constructed Hermite functionals form a complete orthogonal system for square-integrable functionals on Wiener space. The resulting Cameron–Martin expansion supplied the functional-analytic foundation for Wiener’s Gaussian chaos decomposition.

The computational use of polynomial chaos in stochastic engineering developed much later. Roger Ghanem and Pol Spanos applied Hermite expansions and stochastic finite-element formulations to systems with Gaussian uncertainty. Dongbin Xiu and George Karniadakis subsequently established the generalized polynomial chaos framework, in which the polynomial family is matched to the probability measure of the input rather than being restricted to the Gaussian case.

Basis construction

The classical Wiener expansion uses Hermite polynomials because they are orthogonal under the Gaussian measure. If the components of (\boldsymbol{\xi}) are independent standard Gaussian variables, a multivariate basis is formed by tensor products:

[ \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi})

\prod_{j=1}^{d} H_{\alpha_j}(\xi_j). ]

Generalized polynomial chaos replaces the Gaussian measure with the distribution of the standardized input. A uniformly distributed input is associated with Legendre polynomials, whose orthogonality measure is constant on a bounded interval. A beta-distributed input leads to Jacobi polynomials, with parameters determined by the exponents in the beta density. A gamma-distributed input leads to Laguerre polynomials, which are orthogonal under a weighted measure on the nonnegative half-line.

These correspondences belong to the broader theory of orthogonal polynomials. For distributions outside the classical families, an orthogonal basis can be defined through the input measure itself. Completeness requires more than the formal existence of moments; it depends on whether polynomials are dense in the corresponding (L^2) space.

Dependence among input variables changes the construction. Tensor-product polynomials remain directly orthogonal only when the joint measure factors into marginal measures. A dependent joint distribution instead requires polynomials orthogonal under the complete joint measure, or an equivalent transformation to independent reference variables.

Truncation and approximation spaces

Computations use a finite index set (\mathcal A):

[ Y_{\mathcal A}

\sum_{\boldsymbol{\alpha}\in\mathcal A} y_{\boldsymbol{\alpha}} \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi}). ]

A total-degree truncation retains all multi-indices satisfying

[ |\boldsymbol{\alpha}|_1

\alpha_1+\cdots+\alpha_d \leq p. ]

For (d) random dimensions and maximum total degree (p), the number of retained terms is

[ \binom{d+p}{p}. ]

This combinatorial growth is a form of the curse of dimensionality. It reflects the size of the approximation space rather than a property unique to any coefficient-computation method. Restricted multi-index sets reduce the number of mixed high-degree terms when the model response depends weakly on high-order interactions.

For an orthogonal projection, the truncation error satisfies

[ \left|Y-Y_{\mathcal A}\right|_{L^2}^2

\sum_{\boldsymbol{\alpha}\notin\mathcal A} |y_{\boldsymbol{\alpha}}|^2 \gamma_{\boldsymbol{\alpha}}. ]

This expression follows from Parseval's identity. It distinguishes approximation error caused by the finite polynomial space from numerical error in estimated coefficients.

Computation of coefficients

An intrusive stochastic Galerkin method substitutes the truncated expansion into the governing equation and projects the residual onto each basis function. For differential equations, this operation converts the stochastic problem into a coupled deterministic system for the coefficient fields. The coupling structure is determined by inner products involving products of basis polynomials and by the manner in which uncertain parameters enter the original operator.

Non-intrusive projection evaluates the deterministic model at selected realizations of the random inputs and approximates the expectation defining each coefficient. Its mathematical basis is numerical integration, including quadrature rules adapted to the orthogonality measure. The deterministic model remains separate from the spectral reconstruction.

Regression instead determines a finite coefficient vector from sampled model outputs. The corresponding design matrix contains evaluations of the basis polynomials at the sample points. Least-squares regression approximates the projection when the sampling measure, weighting, and sample population adequately represent the polynomial space. Regularized variants impose additional structure when the retained basis is large relative to the available model evaluations.

Statistical interpretation

When the constant basis function is (\Psi_{\boldsymbol 0}=1), the mean is the constant coefficient:

[ \mathbb E[Y]=y_{\boldsymbol 0}. ]

For an orthogonal expansion, the variance is

[ \operatorname{Var}(Y)

\sum_{\boldsymbol{\alpha}\neq\boldsymbol 0} |y_{\boldsymbol{\alpha}}|^2 \gamma_{\boldsymbol{\alpha}}. ]

The multi-index structure also partitions variance according to dependence on particular input variables. Terms whose multi-indices have nonzero entries only in one coordinate represent contributions associated with that coordinate alone. Terms involving several nonzero entries represent interactions among the corresponding inputs. This decomposition connects polynomial chaos coefficients with Sobol indices when the input variables are independent.

Higher moments can be derived from products of expansion terms. Their evaluation involves integrals of three or more basis polynomials, so orthogonality alone does not eliminate every cross-term. The polynomial representation therefore provides direct formulas for low-order statistics while retaining the full approximated random response.

Convergence and limitations

Polynomial chaos convergence depends on the regularity of the model response as a function of its random inputs. Analytic dependence can produce rapid decay of polynomial coefficients, whereas discontinuities in random space commonly yield slower convergence and oscillatory finite expansions. These properties parallel those of deterministic spectral methods.

A finite polynomial expansion also depends on the selected stochastic coordinates. If the random input is originally a field, a representation such as the Karhunen–Loève theorem can reduce it to a sequence of random variables. Truncating that representation introduces an input-model error distinct from the polynomial truncation error.

Polynomial chaos does not make a stochastic model accurate independently of its probability law. The expansion represents uncertainty under the specified input measure, and a change in that measure changes both the orthogonal basis and the resulting coefficients. Numerical convergence within one measure therefore does not establish the validity of the underlying probabilistic assumptions.