Polynomial Chaos Expansion

Resources

Main Idea

For random variables y:ΩRd and a quantity of interest u=f(y), polynomial chaos expansions approximate f. Originally, this was with normally distributed y, where the optimal polynomials are Hermite, but this has been generalized to other distributions (Wiener-Askey schemes).

Intrusive

For solving Partial Differential Equations, this idea can be introduced intrusively as a Stochastic Galerkin method, where x are just more variables making up the domain:

u^(x,t,y)=j=0Puj(x,t)unknownψj(y)knownu^tN(x,t,y,u^,u^x,),ψi=0,i=0,1,,P

This condition of Galerkin Projection (bottom) is a system of P+1 coupled, deterministic PDEs for the unknowns uj(x,t) (for the original scalar PDE).

Non-Intrusive

Alternatively, PCE can be done non-intrusively, using a block box model for f (rather than specifying it implicitly as the solution of a PDE above). With u=f(y), and y following some distribution, we solve for cj such that

u(y)=j=0Pcjψj(y),

where ψj are the associated Wiener-Askey polynomials. For instance, if yU([1,1]), ψj are the Legendre Polynomials. Note that these polynomials are orthogonal with respect to the measure associated with the density of y (p). If also normalized, this means

ψi(y)ψj(y)dp=δij.

Non-intrusively, we sample N instances of y and compute the associated u. We want to choose the coefficients of our expansion so that it matches this data {(yi,ui)}i=1N. For instance, we may wish to minimize the square of the error (Least Squares):

minci=1N(uiu(yi;c))2.

This is solved through the normal equations when N>P+1, creating the standard N×(P+1) measurement matrix Ψ and N×1 solution vector u,

\Psi { #T} \Psi \, \mathbf{c} = \Psi ^ T \mathbf{u} .

This system should be solved with factorization that makes use of the positive definite structure (Singular Value Decomposition or QR Decomposition). Note that in the limit of N, \Psi { #T} \Psi \to \mathbf{I}, due to the aforementioned orthonormality.

Compressive Sensing

The above least squares requires N>P+1, but for a d-dimensional y, and polynomials of maximum total order p,

P+1=(p+dd)=(p+d)!p!d!.

The total number of terms grows rapidly in both p and d. Thus, we may need very many samples in order to even be determined (let alone sufficiently overdetermined).