These notes are still quite rough - I’m still editing them at the moment, and will update with a revised version in a couple days.

Table of Contents

• A Short Protocol for Polynomial Chaos Expansion (PCE)
  • Setting up a surrogate PCE model
  • Using the PCE model for Bayesian inference
• Background
  • Hilbert Space and Inner Products
  • Multivariate Polynomial Basis
    • Tensor product construction
  • Computing PCE coefficients
    • Procedure for PCE analysis
  • Evaluating PCE accuracy
    • Curse of dimensionality
  • Compressive sensing approaches
    • Least Angle Regression (LAR)
  • Sobol decomposition for PCE
    • General Sobol method
    • Sobol method for PCEs

A Short Protocol for Polynomial Chaos Expansion (PCE)

The following are methods to solve a forward problem: Monte Carlo integration, point collocation, pseudo-spectral projection, and intrusive Galerkin

The approach of PCE is to model a deterministic forward problem $u(q)$ by an expansion: \(u(q) \approx \hat{u}(q) = \sum_{i=1}^{N} c_i P_i (\xi_i)\) Where, $q$ are parameters and $P_i$ are polynomials that are orthogonal relative to the probability density function of the random variable $\hat{u}(q)$. $c_i$ are called the Fourier coefficients.

Setting up a surrogate PCE model

1. Declare priors for parameters in model-space ($q_i$) and PCE-space ($\xi_i$)
2. Set up transform from model-space to PCE-space. For example,

Let $q_i \sim \text{Uniform}(a, b)$ and $\xi_i \sim \text{Uniform}(-1, 1)$. Then, \(q_i = F^{-1}_i (T(\xi_i))\) where $F_{i}^{-1}$ is the inverse cumulative distribution function (cdf) of $q_i$, and $T$ is the cdf of $\xi_i$.

1. Construct an $N_P$-dimensional basis of orthogonal polynomials, where

\(N_P = \frac{(N+P)!}{N! P!} - 1\) Where, $N$ is the total number of random variables $\xi_i$ and $P$ the maximum order of the polynomials.

1. Generate $M$ parameter sets sampled from the priors of $\mathbf{\xi}$, and arrange them into a $(M \times (N_P + 1))$-dimensional matrix $\mathbf{H}$.

1. Evaluate the PCE model at the selected points to yield a $M \times 1$ output vector

1. Estimate the Fourier coefficients based on

$\mathbf{c}$ can be found by minimizing the 2-norm: \(|| \mathbf{z} - \mathbf{H}\mathbf{c} ||\)

for which the least-squares solution is: \(\mathbf{\hat{c}} = (\mathbf{H^T}\mathbf{H})^{-1}\mathbf{H^T} \mathbf{z}\)

Using the PCE model for Bayesian inference

We are interested in estimating the posterior distribution \(\pi(\mathbf{q} \mid \mathbf{d}) \propto \mathcal{L}(\mathbf{d} - \mathbf{z}(t; \mathbf{\xi})) \cdot \pi_0(\mathbf{q}) \\\)

Where, $\mathcal{L}$ is the likelihood, $\mathbf{d}$ is the data, $\mathbf{q}$ are parameters in model-space, and $\pi_0$ are prior distributions, also in model-space.

1. Set model-space priors $\pi_0(\mathbf{q})$
2. Construct PCE model $\mathbf{z}(\mathbf{\xi})$ over selected parameter space.
3. Draw samples $\mathbf{ (q^{k}) }^{K}_{k=1}$ from the priors;

• for k = 1:K
  • Evaluate prior and likelihood
  • For importance sampling:
    • Evaluate posterior/prior ratios for each sample, then
    • Normalize resulting weights

1. Transform samples $\mathbf{xi}$ back to model-space to recover posterior distributions.
2. Find MAP estimates by: \(\hat{q}_{\text{MAP}} = \argmin_\mathbf{q} (-\pi(\mathbf{q} \mid \mathbf{d}))\)

Background

Hilbert Space and Inner Products

PCE is a Hilbert space technique equivalent to Fourier series expansions for periodic signals. Hilbert space is a generalization of Euclidean spaces to spaces with possibly infinite dimensions. Hilbert spaces have a defined inner product operator, which thus defines distances between objects in that space, as well as orthogonality (i.e. if the inner product between objects is zero). The dot product gives both length and angles:

More generally, a vector space that has a qualifying inner product is called an inner product space. Every such space is also a Hilbert space. The dot product is one type of inner product, but more are possible. The criteria for inner products are:

1. the operation is commutative, 
2. linear, 
3. positive definite 

More on Hilbert spaces

The coefficients $c_i$ are found by \(c_i = \frac{<u_i , P_i>}{<P_i , P_i>}\)

One method to compute these is by quadrature integration, which involves approximating the integral of a function by that of a weighted sequence of polynomials.

Multivariate Polynomial Basis

Consider the following PCE expansion for random output $Y = M(\mathbf{X})$,

where

• $Y$ are the observed data,
• $M$ is our model,
• $\mathbf{X}$ are model parameters (aka ‘input variables’)
• $m$ is the dimensionality of $\mathbf{X}$
• $\Phi_i$ is a PCE basis, comprised of multivariate orthogonal/orthonormal polynomials
• $c_i$ are Fourier coefficients/coordinates

The prior probability density functions (PDFs) of (random) input parameters $\mathbf{X}$ is given by

where $\mathbf{x}$ is a single realization (draw/sample).

A basis of orthogonal polynomials satisfies the property that the inner product of any two polynomials in the basis is zero. In the univariate case, let $P_{k}^i, \ k \in \mathcal{N}$ be the $k$-th polynomial component for the input variable $X_i$. Note $\mathcal{N}$ is used here (and below) to denote the set of positive integers, or ‘natural numbers.’ Then,

When $X_i \sim \text{Uniform}(-1, 1)$, $P_i$ are called Legendre polynomials, and Hermite polynomials if $X_i \sim \text{Normal}(0, 1)$.

An orthonormal basis can be constructed by normalizing the $P^{i}_j$:

Tensor product construction

The above tells us how to construct a basis for a single input parameter, but how do we construct a basis for the full set of $\mathbf{X}$?

To do this, we define a multivariate polynomial as the product of univariate polynomials corresponding to each input parameter $x_i$.

where, $\alpha = (\alpha_1, \dots \alpha_m)$ are multi-indices (degree of the polynomial in each variable, i.e. partial degree in each dimension).

Computing PCE coefficients

PCEs can be written as the sum of a truncated series and a residual term:

where:

• $\mathbf{C} = { c_\alpha , \ \alpha \in \mathcal{A} } = { c_0, \dots c_{P-1} }$ contains the $P$ unknown coefficients. See below for a definition of $\mathcal{A}$.
• $\mathbf{\Phi}(x) = { \Phi_0 (x) , \dots \Phi_{P-1} (x) }$ holds the $P$ corresponding polynomial bases

Procedure for PCE analysis

The idea is to find coefficients $\mathbf{\hat{C}}$ that the minimizes the difference between the outputs of the true model $M$ and a (sparse) PCE model.

1. Select a ‘truncation scheme’ - that is, after computing a ‘full order’ model (see equation for ‘cardinality’), we only keep a subset of terms to obtain a ‘sparse’ PCE.

A common truncation scheme is to set some maximal order. For instance, with an $m$-dimensional input space and maximum polynomial order $p$, we can write:

1. Simulate the model response for $n$ parameter sets sampled from the input space.

1. Construct and evaluate the polynomial basis over the selected parameters and hold the results in an “experimental matrix” $\mathbf{A}$:

1. Finally, we solve the linear system for coefficients: \(\mathbf{\hat{C}} = (\mathbf{A^T}\mathbf{A})^{-1}\mathbf{A^T}\mathbf{M}\)

Evaluating PCE accuracy

We also want to optimize the polynomial order, or ‘truncation scheme.’ We could do this by:

• Minimizing the MSE (on a different set of data) between the true model and PCEs of different orders, but this would require running the true model many times, which is computationally intensive.
• Or, we can just re-use the same data used for calibration (i.e. compute the empirical error), but this would result in overfitting.

Alternatively, we can use Leave-one-out (LOO) cross-validation:

1. Select an ‘experimental design’ $\mathcal{X} = { x^j ,\ j=1, \dots n}$, i.e. predictions of the true model with a low-discrepancy set as input
2. Build PCEs using all points but one, e.g. leaving out the $i$-th point: $\mathcal{X}_{i \neq j}$.
3. Compute the error at the held-out $i$-th point, which is the local error at that point.
4. Repeat this for all points $x^i$.

To get the LOO error without re-computing $n$ PCEs (for each validation point), we can instead run a single computation via a corrected LOO error:

Where, $h_i$ is the $i$-th diagonal term of the matrix \(\mathbf{A}(\mathbf{A^T}\mathbf{A})^{-1}\mathbf{A^T}\)

This effectively bypasses the last step of LOO. That is, we get the same LOO error after just a single PCE analysis (i.e. select a truncation order, construct a PCE basis, fit Fourier coefficients).

Curse of dimensionality

For $m$ model parameters and polynomial order $p$, the cardinality of the truncation scheme (number of terms) is:

Compressive sensing approaches

Least Angle Regression (LAR)

Least angle regression is one approach to perform this regularized minimization, and does so by solving the LASSO problem for different values of $\lambda$ simultaneously. A Python implementation in scikit-learn is available.

LAR is an algorithm that finds coefficients for linear regressions:

1. Initially, all coefficients are set to zero.
2. For the coefficient whose value is most correlated to the norm of the residual, move it in the direction of correlation until the correlation of some other coefficient catches up (i.e. since the correlation of the first coefficient will decrease as it approaches its least-square value).
3. Extend (i.e. modify the value of) the two coefficients in the direction that is equiangular to (i.e. ‘in the middle’ between) their individual correlations.
4. Repeat for all coefficients.

Sobol decomposition for PCE

General Sobol method

The Sobol method looks at the contribution of individual parameters $x_i \in \mathbf{X}$ to the total variance of the model output $M(\mathbf{X})$. The central quantity are the Sobol indices, or first-order effects, which describe the partial variance due to a single parameter. The $u$-th first-order Sobol index is given by:

Where $D_u$ and $D$ are the partial and total variances, respectively. The so-called ‘total’ order Sobol index of the $u$-th parameter is given by one minus the partial variance due to all parameters besides the $u$-th parameter.

Sobol method for PCEs

For PCE, the first-order Sobol index for the $i$-th parameter can be computed as the sum of products between all of the $i$-th coefficients and polynomial bases.

The formula for total-order indices is analogous to the general method: we sum products for all but the $i$-th terms. \(\widetilde{\mathcal{A}}_u = \{ \alpha \in \mathcal{A} \mid \alpha_i = 0 \ \forall \ i \neq \mathbf{u} \}\)