Generalized Polynomial Chaos: A UQ Framework

• Generalized Polynomial Chaos (gPC) is a framework that approximates square-integrable outputs using orthonormal polynomials matched to the input probability law.
• It computes expansion coefficients via projection, quadrature, regression, collocation, or Galerkin methods, enabling direct access to statistical moments.
• gPC extends to handle dependent variables, sparse sampling, and diverse applications ranging from aerodynamic CFD to Bayesian inversion and control under uncertainty.

Generalized polynomial chaos (gPC) is a spectral-like representation of uncertainty in which a square-integrable random output is expanded in orthogonal or orthonormal polynomials matched to the probability law of the uncertain inputs. In its standard form, one approximates a computational model output, a stochastic state, or a parametric solution map by a finite polynomial expansion in a basis adapted to the input measure, and then computes the expansion coefficients by projection, quadrature, regression, collocation, or Galerkin projection. Across the literature, gPC appears as a framework for uncertainty quantification, stochastic Galerkin discretization, Bayesian inversion, calibration, control under parametric uncertainty, and reduced-order surrogate construction (Savin et al., 2015, Tagade et al., 2012, Sutulovic et al., 2024).

1. Core formulation

For a computational model mapping uncertain inputs $\xi=(\xi_1,\dots,\xi_D)^\top$ to an output of interest $Y=g(\xi)$, gPC seeks a finite-dimensional approximation

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$

Here $\Psi_i(\xi)$ is the $i$-th multivariate orthonormal polynomial of total degree $\le p$, $c_i=E[Y(\xi)\Psi_i(\xi)]$ are the gPC coefficients, $\xi$ is the random input vector with joint PDF $\rho(\xi)$, and the expansion projects the model output on a finite-dimensional polynomial subspace of $L^2(\rho)$ (Savin et al., 2015).

Equivalent multi-index formulations are standard. For a random output $Y=g(\xi)$0, one writes

$Y=g(\xi)$1

where $Y=g(\xi)$2, $Y=g(\xi)$3, and $Y=g(\xi)$4 (Sutulovic et al., 2024). With total-order truncation, the cardinality of the retained basis is

$Y=g(\xi)$5

a formula that recurs in forward UQ, inverse modeling, and optimal control under uncertainty (Karagiannis et al., 2019, Aoyama et al., 2021).

Because orthogonality isolates the constant mode, the first two moments are directly accessible from the coefficients. In the orthonormal setting, $Y=g(\xi)$6 and $Y=g(\xi)$7 for the truncated generalized polynomial chaos expansion (Rahman, 2018). In control-oriented formulations, the mean is the zeroth coefficient $Y=g(\xi)$8 and the covariance is assembled from the higher modes as $Y=g(\xi)$9 (Aoyama et al., 2021).

2. Orthogonal bases and measure consistency

The choice of polynomial family is determined by the input distribution. The Wiener–Askey correspondence appears repeatedly: Hermite polynomials for Gaussian random variables, Legendre polynomials for uniform random variables, Jacobi polynomials for Beta-type random variables, and Laguerre polynomials for Gamma random variables (Sutulovic et al., 2024). In aerodynamic computations with Beta-distributed uncertainties, the marginal distributions determine Gauss–Jacobi polynomials rescaled to the parameter interval, and the multivariate basis is assembled as a product

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$0

with orthonormality

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$1

This tensor-product structure is the classical construction for mutually independent inputs (Savin et al., 2015).

A major extension concerns dependent random variables. A generalized polynomial chaos expansion in dependent random variables can be built from measure-consistent multivariate orthonormal polynomials, and unlike existing PCEs, whether classical or generalized, no tensor-product structure is assumed or required (Rahman, 2018). In that setting, monomials are ordered by increasing total degree and a multivariate Gram–Schmidt procedure in $$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$2 produces orthonormal polynomials $$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$3. The resulting construction yields an orthogonal decomposition of polynomial spaces, completeness under finite-moment assumptions, and mean-square convergence of truncated expansions to the correct limit (Rahman, 2018).

This distinction between independent-input tensorization and dependent-input measure consistency is structurally important. For product measures, tensorized bases exploit separability directly. For non-product-type densities, the basis itself must encode dependence, otherwise orthogonality under the true joint density is lost (Rahman, 2018).

3. Intrusive and non-intrusive coefficient construction

The coefficients of a gPC approximation are projections of the stochastic output onto the polynomial basis. In non-intrusive formulations, these projections are approximated from deterministic model evaluations. A standard cubature formula is

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$4

where $$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$5 is a quadrature or cubature rule (Savin et al., 2015). Collocation and regression variants treat the forward model as a black box, evaluate it at selected samples $$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$6, and reconstruct the coefficients by discrete inner products or least-squares solves (Sutulovic et al., 2024). Pseudo-spectral projection is used in both classical and mapped settings, with the expansion coefficients assembled from weighted evaluations of the original deterministic solver (Georg et al., 2020).

Intrusive formulations substitute the truncated gPC expansion into the governing equations and enforce orthogonality of the residual against every basis polynomial. For quasilinear hyperbolic conservation laws,

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$7

a symmetrization based on the left-eigenvector matrix of the Jacobian yields a symmetric form, and stochastic Galerkin projection produces a deterministic coupled system for the gPC coefficients. In the one-dimensional setting, the resulting deterministic gPC Galerkin system is proved to be symmetrically hyperbolic, which permits higher-order and path-conservative finite volume WENO discretization in space together with a third-order total variation diminishing Runge–Kutta method in time (Wu et al., 2016). The same paper extends the method to two-dimensional quasilinear systems via operator splitting while carrying over symmetric hyperbolicity of the one-dimensional subsystems (Wu et al., 2016).

Intrusive projection also appears in lower-dimensional dynamical systems. For a parametric ODE

$$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$8

Galerkin projection yields a deterministic system of ODEs for the time-dependent spectral coefficients $$Y(\xi)\approx \hat Y_p(\xi)=\sum_{i=0}^{P} c_i\,\Psi_i(\xi),\qquad P+1=\frac{(p+D)!}{p!\,D!}.$$9 (Sutulovic et al., 2024). In optimal control under parametric uncertainty, the lifted coefficient vector $\Psi_i(\xi)$0 obeys a deterministic evolution law $\Psi_i(\xi)$1 obtained by the same principle (Aoyama et al., 2021).

4. Sampling, sparsity, and high-dimensional acceleration

The principal computational difficulty is the rapid growth of basis size and sampling cost with dimension and polynomial order. Direct tensor-product constructions suffer from the “curse of dimensionality,” and tensor-product quadrature in $\Psi_i(\xi)$2 dimensions uses $\Psi_i(\xi)$3 points, which is costly for $\Psi_i(\xi)$4 (Savin et al., 2015, Sutulovic et al., 2024). One standard mitigation is sparse quadrature. Starting from one-dimensional Gauss–Jacobi-Lobatto rules, Smolyak’s algorithm constructs sparse grids of level $\Psi_i(\xi)$5 with $\Psi_i(\xi)$6 points while remaining exact for all polynomials up to total degree $\Psi_i(\xi)$7 (Savin et al., 2015).

A second mitigation is coefficient sparsity. In aerodynamic computations, observed “sparsity-of-effects” means that in many CFD outputs low-order and main-effect terms dominate, high-order cross-terms are negligible, and the gPC coefficient vector is compressible (Savin et al., 2015). This motivates compressed sensing formulations. With measurement matrix $\Psi_i(\xi)$8 and data vector $\Psi_i(\xi)$9, one recovers the coefficient vector by Basis Pursuit Denoising,

$i$0

where $i$1 accounts for truncation or modeling error (Savin et al., 2015). Under incoherence or Restricted Isometry Property conditions and if $i$2, the sparse vector is stably recovered (Savin et al., 2015). A more general framework enhances sparsity by iterative rotations of the random inputs, using an alternating-direction method and PCA on surrogate gradients to identify rotated coordinates in which the expansion is sparser; the paper also notes that for non-Gaussian bases excessive coherence can negate sparsity gains (Yang et al., 2017).

Model-reduction and derivative-enhanced constructions provide additional acceleration. A weighted reduced basis method can be integrated with non-intrusive gPC so that the reduced surrogate is used at quadrature nodes, with rigorous error control and reported speed-ups of $i$3–$i$4 as parameter dimension grows from $i$5 to $i$6 (Jiang et al., 2016). Sensitivity-enhanced generalized polynomial chaos enriches least-squares regression with adjoint-computed gradients of the quantity of interest, and the reported computational cost scales as $i$7, instead of $i$8 of the standard LSQ formulation (Kantarakias et al., 2022). For time-average quantities of chaotic systems, the same se-gPC idea, coupled to the adjoint of the shadowing operator in the frequency domain, yields in its simplest form a computational cost that is independent of the number of random variables (Papadakis et al., 2023).

5. Bayesian inversion, calibration, and control

A prominent use of gPC is surrogate-assisted inference. In hydrologic parameter inversion for CLM4, the surrogate model is built in a Bayesian framework with a spike-and-slab prior on the gPC coefficients, allowing basis selection uncertainty to be quantified through marginal inclusion probabilities $i$9 (Karagiannis et al., 2019). Two surrogate summaries are emphasized: Bayesian model averaging, based on ergodic means of the sampled coefficients, and the Median Probability Model, which retains basis terms with posterior inclusion probability greater than $\le p$0 and then refits on the reduced basis (Karagiannis et al., 2019). The surrogate is then inserted into a block-Gibbs/Metropolis–Hastings inversion scheme to recover marginal posteriors and credible intervals for the hydrologic parameters (Karagiannis et al., 2019).

A related calibration framework combines Karhunen–Loève expansion and gPC. Independent stationary Gaussian processes with uncertain hyper-parameters represent parametric and structural uncertainties, the KL expansion of a prior Gaussian process is projected on a generalized Polynomial Chaos basis, and intrusive Galerkin projection is utilized to calculate the associated coefficients of the simulator output (Tagade et al., 2012). Bayesian inference then updates the prior probability distribution of the generalized Polynomial Chaos basis, which along with the chaos expansion coefficients represent the posterior probability distribution (Tagade et al., 2012).

Conditional constructions sharpen this strategy when data about the uncertain field are available. In PDE-based parameter estimation for Darcy flow, Gaussian Process regression conditions the KL expansion of the unknown space-dependent parameter field, reducing the stochastic dimension from $\le p$1 to $\le p$2, and a conditional gPC surrogate of the PDE state is then built in the reduced random space (Li et al., 2019). The resulting parameter estimate exactly matches the parameter measurements, and the conditional gPC surrogate can also be used to estimate the states variance, which can in turn guide data acquisition (Li et al., 2019).

In control, gPC is used not only to propagate uncertainty but also to express higher order moments of the uncertainty distribution without any Gaussianity assumption (Aoyama et al., 2021). Receding-horizon Differential Dynamic Programming under parametric uncertainty lifts the stochastic system to a deterministic coefficient system, computes mean and covariance from the gPC coefficients, and incorporates chance constraints through a deterministic approximation involving the state covariance (Aoyama et al., 2021). The receding-horizon implementation resets the higher-order coefficients to zero at each measurement update, a step reported to prevent unbounded growth of variance and ensure that the chance-constraints remain satisfied in closed loop (Aoyama et al., 2021).

6. Generalizations, validation, and newer computational directions

Several extensions alter the basis itself rather than only the coefficient solver. Conformally mapped polynomial chaos expansions introduce a real-to-real conformal map, define a transformed density $\le p$3, and build orthonormal polynomials with respect to the transformed weight (Georg et al., 2020). The mapped basis remains orthonormal with respect to the original PDF, so moments and Sobol sensitivity indices remain directly computable from the coefficients, while the transformed integrand can exhibit a larger Bernstein ellipse of analyticity and thus a faster convergence constant (Georg et al., 2020).

Another extension establishes rigorous a posteriori validation of gPC expansions. In a Banach-space formulation $\le p$4, one defines computable quantities $\le p$5, $\le p$6, and a nonlinear bound $\le p$7; if

$\le p$8

then a Newton–Kantorovich-type argument yields existence, uniqueness, and an explicit radius $\le p$9 for a true solution near the numerical one (Breden, 2022). This framework is applied to random invariant periodic orbits in the Lorenz system and to branches and $c_i=E[Y(\xi)\Psi_i(\xi)]$0-dimensional manifolds of steady states of the Swift–Hohenberg equation (Breden, 2022).

A still newer direction concerns countably-parametric holomorphic maps and quantum representations. For $c_i=E[Y(\xi)\Psi_i(\xi)]$1-holomorphic maps on $c_i=E[Y(\xi)\Psi_i(\xi)]$2, best $c_i=E[Y(\xi)\Psi_i(\xi)]$3-term truncations of Chebyshev-gPC expansions admit approximation rates of order $c_i=E[Y(\xi)\Psi_i(\xi)]$4 in the $c_i=E[Y(\xi)\Psi_i(\xi)]$5 norm and of order $c_i=E[Y(\xi)\Psi_i(\xi)]$6 in the $c_i=E[Y(\xi)\Psi_i(\xi)]$7 norm (Aftab et al., 2 Jun 2025). Parameterized quantum circuit encodings of these truncations are then constructed through tensorization of univariate PQCs for Chebyshev polynomials and linear combination of unitaries, with width $c_i=E[Y(\xi)\Psi_i(\xi)]$8 and depth $c_i=E[Y(\xi)\Psi_i(\xi)]$9 (Aftab et al., 2 Jun 2025). This suggests that gPC has become a structural interface not only for classical UQ algorithms but also for rigorous quantum approximation schemes.

7. Applications and reported performance

Published applications span aerodynamic CFD, land-surface modeling, Maxwell problems, neural dynamics, heterogeneous neural networks, hyperbolic conservation laws, singular-potential wave equations, and chaotic dynamics. The reported results are domain-specific, but they consistently use the same building blocks: basis adaptation to the input law, truncation in polynomial degree, and a coefficient solver chosen to match the forward model and the available data.

Domain	gPC formulation	Reported result
Steady RANS around RAE2822 airfoil in transonic flow	3-D Jacobi polynomials up to total order $\xi$0; full tensor, sparse Smolyak, and compressed sensing surrogates	Tensor grid: $\xi$1 CFD runs; sparse grid: $\xi$2 runs; CS: $\xi$3 runs; means $\xi$4 and standard deviations $\xi$5 agree within $\xi$6; $\xi$7 contributes $\xi$8 of output variance (Savin et al., 2015)
CLM4 hydrologic parameter inversion at the US-ARM site	Bayesian surrogate with gPC variable selection, Gibbs sampling, and posterior inversion	Only a handful of PC terms are important, and the method obtains tight posterior distributions for those parameters (Karagiannis et al., 2019)
Quasilinear hyperbolic conservation laws	Intrusive gPC stochastic Galerkin with symmetric hyperbolicity, WENO, and TVD Runge–Kutta	In 2D Riemann tests, gPC-SG with $\xi$9 on $\rho(\xi)$0 meshes reproduces mean and $\rho(\xi)$1 contours in close agreement with a $\rho(\xi)$2-point stochastic collocation run on $\rho(\xi)$3 (Wu et al., 2016)
Neural dynamics	Intrusive and non-intrusive gPC for Hindmarsh–Rose, Jansen–Rit, and Epileptor models	HR collocation with $\rho(\xi)$4 yields a $\rho(\xi)$5 speedup over MC; Epileptor long simulations reduce from $\rho(\xi)$6 for MC to $\rho(\xi)$7 for gPC collocation (Sutulovic et al., 2024)
Aerodynamic UQ with $\rho(\xi)$8 stochastic parameters	Sensitivity-enhanced generalized polynomial chaos	Reported cost scaling $\rho(\xi)$9 instead of $L^2(\rho)$0; for transonic inviscid flow around NACA0012, se-gPC uses $L^2(\rho)$1 evaluations for $L^2(\rho)$2, $L^2(\rho)$3 for $L^2(\rho)$4, and $L^2(\rho)$5 for $L^2(\rho)$6 (Kantarakias et al., 2022)

Other studies emphasize different strengths. In conformally mapped gPC for Maxwell’s source problem with random input data, mapped gPC converges roughly $L^2(\rho)$7 faster in polynomial degree than standard gPC in the $L^2(\rho)$8-error of $L^2(\rho)$9, and to reach a given $Y=g(\xi)$00 accuracy one needs about half as many finite-element solves (Georg et al., 2020). In heterogeneous neural networks, gPC combined with anchored ANOVA reduces the representation of four heterogeneities to $Y=g(\xi)$01 basis functions per state variable, and ANOVA with $Y=g(\xi)$02 uses $Y=g(\xi)$03 collocation points versus $Y=g(\xi)$04 sparse-grid points versus $Y=g(\xi)$05 Monte Carlo samples while yielding visually indistinguishable mean and variance (Choi et al., 2016). In wave equations with a singular potential, gPC coupled with a consistent Chebyshev collocation method is used to determine critical parameters efficiently and accurately, including a Klein–Gordon critical value $Y=g(\xi)$06 and a sine–Gordon critical value $Y=g(\xi)$07 in the reported test configurations (Chakraborty et al., 2011).

Taken together, these results show that generalized polynomial chaos is not a single algorithm but a family of measure-adapted polynomial representations that support multiple computational regimes: intrusive Galerkin discretization, non-intrusive projection and regression, sparse recovery, Bayesian variable selection, reduced-basis acceleration, sensitivity enrichment, validated numerics, and dimension-reduced inverse modeling. The choice among these regimes is controlled by the input measure, the regularity of the response map, the sparsity pattern of the coefficients, and the computational structure of the underlying forward model.