Polynomial chaos expansions

A recent mathematical framework for estimating uncertainty is that of polynomial chaos expansions (Xiu and Hesthaven, 2005). Polynomial chaos expansions can be seen as a subset of polynomial approximation methods. For a review of polynomial chaos expansions see (Xiu, (2010)). Polynomial chaos expansions are much faster than (quasi-)Monte Carlo methods as long as the number of uncertain parameters is relatively low, typically smaller than about twenty (Crestaux et al.,2009). This is the case for many neuroscience models, and even for models with a higher number of uncertain parameters, the analysis could be performed for selected subsets of the parameters.

The general idea behind polynomial chaos expansions is to approximate the model \(U\) with a polynomial expansion \(\hat{U}\):

\[U \approx \hat{U}(\boldsymbol{x}, t, \boldsymbol{Q}) = \sum_{n=0}^{N_p - 1} c_n(\boldsymbol{x}, t) \boldsymbol{\phi}_n (\boldsymbol{Q}),\]

where \(\boldsymbol{\phi}_n\) denote polynomials and \(c_n\) denote expansion coefficients. The number of expansion factors \(N_p\) is given by

\[N_p = \binom{D+p}{p},\]

where \(p\) is the polynomial order. The number of expansion coefficients in the multivariate case (\(D>1\)) is greater than the polynomial order. This is because the multivariate polynomial is created by multiplying univariate polynomials together. The polynomials \(\phi_n(\boldsymbol{Q})\) are chosen so they are orthogonal with respect to the probability density function \(\rho_{\boldsymbol{Q}}\), which ensures useful statistical properties.

When creating the polynomial chaos expansion, the first step is to find the orthogonal polynomials \(\boldsymbol{\phi}_n\), which in Uncertainpy is done using the so called three-term recurrence relation (Xiu, 2010). The next step is to estimate the expansion coefficients \(c_n\). The non-intrusive methods for doing this can be divided into two classes, point-collocation methods and pseudo-spectral projection methods, both of which are implemented in Uncertainpy.

Point collocation is the default method used in Uncertainpy. This method is based on demanding that the polynomial approximation is equal to the model output evaluated at a set of collocation nodes drawn from the joint probability density function \(\rho_{\boldsymbol{Q}}\). This demand results in a set of linear equations for the polynomial coefficients \(c_n\), which can be solved by the use of regression methods. The regression method used in Uncertainpy is Tikhonov regularization (Rifkin and Lippert, 2007).

Pseudo-spectral projection methods are based on least squares minimization in the orthogonal polynomial space, and finds the expansion coefficients \(c_n\) through numerical integration. The integration uses a quadrature scheme with weights and nodes, and the model is evaluated at these nodes. The quadrature method used in Uncertainpy is Leja quadrature, with Smolyak sparse grids to reduce the number of nodes required (Narayan and Jakeman, 2014; Smolyak, 1963). Pseudo-spectral projection methods are only used in Uncertainpy when requested by the user.

Several of the statistical metrics of interest can be obtained directly from the polynomial chaos expansion \(\hat{U}\). The mean is simply

\[\mathbb{E}[U] \approx \mathbb{E}[\hat{U}] = c_0,\]

and the variance is

\[\mathbb{V}[U] \approx \mathbb{V}[\hat{U}] = \sum_{n=1}^{N_p - 1} \gamma_n c_n^2,\]

where \(\gamma_n\) is a normalization factor defined as

\[\gamma_n = \mathbb{E}\left[\boldsymbol{\phi}_n^2(\boldsymbol{Q})\right].\]

The first and total order Sobol indices can also be calculated directly from the polynomial chaos expansion (Sudret, 2008; Crestaux et al.,2009). On the other hand, the percentiles and prediction interval must be estimated using \(\hat{U}\) as a surrogate model, and then perform the same procedure as for the (quasi-)Monte Carlo methods.