PolynomialChaos

Computes and evaluates polynomial chaos surrogate model.

Overview

Polynomial chaos is a surrogate modeling technique where a quantity of interest (QoI) that is dependent on input parameters is expanded as a sum of orthogonal polynomials. Given a QoI QQ dependent on a set of parameters ξ\vec{\xi}, the polynomial chaos expansion (PCE) is:

Q(ξ)=i=1PqiΦi(ξ),Q(\vec{\xi}) = \sum_{i=1}^{P}q_i\Phi_i(\vec{\xi}) ,(1)

where PP is the multidimensional polynomial order and qiq_i are coefficients that are to be computed. These coefficients can be found using intrusive and non intrusive techniques. The intrusive technique is quite difficult to generalize and very computationally demanding. Since the polynomial basis is orthogonal, a non intrusive technique is developed where the coefficients are found by performing a Galerkin projection and integrating:

qi=Q(ξ)Φi(ξ)Φi(ξ),Φi(ξ),q_i = \frac{\left\langle Q(\vec{\xi})\Phi_i(\vec{\xi})\right\rangle}{\left\langle\Phi_i(\vec{\xi}),\Phi_i(\vec{\xi})\right\rangle},(2)

where,

a(ξ)b(ξ)=a(ξ)b(ξ)f(ξ)dξ.\left\langle a(\vec{\xi}) b(\vec{\xi}) \right\rangle = \int_{-\infty}^{\infty}a(\vec{\xi}) b(\vec{\xi}) f(\vec{\xi}) d\vec{\xi} .

The weight function (f(ξ)f(\vec{\xi})) and bases (Φi(ξ)\Phi_i(\vec{\xi})) are typically products of one-dimensional functions:

f(ξ)=d=0Dfd(ξd),f(\vec{\xi}) = \prod_{d=0}^{D}f_d(\xi_d) ,Q(ξ)=i=1Pqid=0Dϕkd,id(ξd).Q(\vec{\xi}) = \sum_{i=1}^{P}q_i\prod_{d=0}^{D}\phi^d_{k_{d,i}}(\xi_d) .

The weighting functions are defined by the probability density function of the parameter and the polynomials are based on these distributions, Table 1 is a list of commonly used distributions and their corresponding orthogonal polynomials.

Table 1: Common probability density functions and their corresponding orthogonal polynomials

DistributionDensity Function (fd(ξ)f_d(\xi))Polynomial (ϕi(ξ)\phi_i(\xi))Support
Normal12πeξ2/2\frac{1}{2\pi}e^{-\xi^2/2}Hermite[,][-\infty, \infty]
Uniform12\frac{1}{2}Legendre[1,1][-1, 1]
Beta(1ξ)α(1+ξ)β2α+β+1B(α+1,β+1)\frac{(1-\xi)^{\alpha}(1+\xi)^{\beta}}{2^{\alpha+\beta+1}B(\alpha+1,\beta+1)}Jacobi[1,1][-1,1]
Exponentialeξe^{-\xi}Laguerre[0,][0,\infty]
GammaξαeξΓ(α+1)\frac{\xi^{\alpha}e^{-\xi}}{\Gamma(\alpha+1)}Generalized Laguerre[0,][0,\infty]

Computing Coefficients

There are currently two techniques to compute the coefficients in Eq. (1), the first performs the integration described by Eq. (2). For Monte Carlo sampling, the integration is in the form

qi=1Φi(ξ),Φi(ξ)1Nmcn=1NmcQ(ξn)Φi(ξn).q_i = \frac{1}{\left\langle\Phi_i(\vec{\xi}),\Phi_i(\vec{\xi})\right\rangle} \frac{1}{N_{\mathrm{mc}}}\sum_{n=1}^{N_{\mathrm{mc}}} Q(\vec{\xi}_n)\Phi_i(\vec{\xi}_n) .

And for using a QuadratureSampler

qi=1Φi(ξ),Φi(ξ)n=1NqwnQ(ξn)Φi(ξn).q_i = \frac{1}{\left\langle\Phi_i(\vec{\xi}),\Phi_i(\vec{\xi})\right\rangle}\sum_{n=1}^{N_q} w_n Q(\vec{\xi}_n)\Phi_i(\vec{\xi}_n) .

The numerical quadrature method is typically much more efficient that than the Monte Carlo method and has the added benefit of exactly integrating the polynomial basis. However, the quadrature suffers from the curse of dimensionality. The naive approach uses a Cartesian product of one-dimensional quadratures, which results in (max(kid)+1)D(\max(k^d_i) + 1)^D quadrature points to be sampled. Sparse grids can help mitigate the curse of dimensionality significantly.

The other technique is using ordinary least-squares (OLS) regression, like in PolynomialRegressionTrainer. In the majority of cases, OLS is more accurate than integration for Monte Carlo sampling, as shown in the figure below.

251002510002510k51251025100
OLSOLS with Regularization (penalty = 1)IntegrationNumber of SamplesCoefficient RRMSE

Comparison between OLS regression and integration methods for computing polynomial chaos coefficients

Generating a Tuple

In polynomial chaos, a tuple describes the combination of polynomial orders representing the expansion basis (kd,ik_{d,i}). Again, the naive approach would be to do a tensor product of highest polynomial order, but this is often wasteful since generating a complete monomial basis is usually optimal. Below demonstrates the difference between a tensor basis and a complete monomial basis:

D=2,kmax=2,Tensor product:Φ0=ϕ01ϕ02,Φ1=ϕ11ϕ02,Φ2=ϕ01ϕ12,Φ3=ϕ11ϕ12Φ4=ϕ21ϕ02,Φ5=ϕ01ϕ22,Φ6=ϕ11ϕ22,Φ7=ϕ21ϕ12,Φ8=ϕ21ϕ22,\begin{aligned} &D=2, \, k_{\max}=2, \, \text{Tensor product:} \\ &\Phi_0 = \phi^1_0\phi^2_0,\, \Phi_1 = \phi^1_1\phi^2_0,\, \Phi_2 = \phi^1_0\phi^2_1,\, \Phi_3 = \phi^1_1\phi^2_1 \\ &\Phi_4 = \phi^1_2\phi^2_0,\, \Phi_5 = \phi^1_0\phi^2_2,\, \Phi_6 = \phi^1_1\phi^2_2,\, \Phi_7 = \phi^1_2\phi^2_1,\, \Phi_8 = \phi^1_2\phi^2_2 , \end{aligned}D=2,kmax=2,Complete monomial:Φ0=ϕ01ϕ02,Φ1=ϕ11ϕ02,Φ2=ϕ01ϕ12,Φ3=ϕ11ϕ12Φ4=ϕ21ϕ02,Φ5=ϕ01ϕ22.\begin{aligned} &D=2, \, k_{\max}=2, \, \text{Complete monomial:} \\ &\Phi_0 = \phi^1_0\phi^2_0,\, \Phi_1 = \phi^1_1\phi^2_0,\, \Phi_2 = \phi^1_0\phi^2_1,\, \Phi_3 = \phi^1_1\phi^2_1 \Phi_4 = \phi^1_2\phi^2_0,\, \Phi_5 = \phi^1_0\phi^2_2 . \end{aligned}

The tuple is generated and stored as matrix in the userobject, below is an example of this matrix with D=3D=3 and kmax=3k_{\max}=3:

kd,i=[010021100032211100000010010210010210321000010010120010120123]k_{d,i} = \left[ \begin{array}{c|ccc|cccccc|cccccccccc} 0 & 1 & 0 & 0 & 2 & 1 & 1 & 0 & 0 & 0 & 3 & 2 & 2 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 2 & 1 & 0 & 0 & 1 & 0 & 2 & 1 & 0 & 3 & 2 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 2 & 0 & 0 & 1 & 0 & 1 & 2 & 0 & 1 & 2 & 3 \end{array} \right]

Computing Statistical Moments

Statistical moments are based on the expectation of a function of the quantity of interest:

E[g(Q(ξ))]=g(Q(ξ))f(ξ)dξ.E\left[g\left(Q(\vec{\xi})\right)\right] = \int_{-\infty}^{\infty}g\left(Q(\vec{\xi})\right)f(\vec{\xi})d\vec{\xi} .

The first four statistical moments, and the most common ones, are defined as:

Mean: μ=E[Q],\text{Mean: } \mu = E[Q] ,Variance: σ2=E[(Qμ)2]=E[Q2]μ2,\text{Variance: } \, \sigma^2 = E\left[\left(Q-\mu\right)^2\right] = E[Q^2] - \mu^2 ,Skewness: Skew=E[(Qμ)3]σ3=E[Q3]3σ2μμ3σ3,\text{Skewness: } \mathrm{Skew} = \frac{E\left[\left(Q-\mu\right)^3\right]}{\sigma^3} = \frac{E[Q^3] - 3\sigma^2\mu - \mu^3}{\sigma^3} ,Kurtosis: Kurt=E[(Qμ)4]σ4=E[Q4]4E[Q3]μ+6σ2μ2+3μ4σ4.\text{Kurtosis: } \mathrm{Kurt} = \frac{E\left[\left(Q-\mu\right)^4\right]}{\sigma^4} = \frac{E[Q^4] - 4E[Q^3]\mu + 6\sigma^2\mu^2 + 3\mu^4}{\sigma^4} .

Because of the orthogonality of the polynomials, mean and variance are trivial to compute:

E[Qpc]=q0,E\left[Q_{\mathrm{pc}}\right] = q_0 ,E[Qpc2]=i=1Pqi2d=1Dϕkd,i2,E\left[Q^2_{\mathrm{pc}}\right] = \sum_{i=1}^{P}q_i^2 \prod_{d=1}^{D}\left\langle\phi_{k_{d,i}}^2\right\rangle ,

where Φi2\left\langle\Phi_i^2\right\rangle is known analytically. The higher order moments are significantly more taxing to compute since it does not take advantage of orthogonality:

E[Qpc3]=i=1Pj=1Pk=1Pqiqjqkd=1Dϕkd,iϕkd,jϕkd,k,E\left[Q^3_{\mathrm{pc}}\right] = \sum_{i=1}^{P}\sum_{j=1}^{P}\sum_{k=1}^{P}q_iq_jq_k\prod_{d=1}^{D}\left\langle\phi_{k_{d,i}}\phi_{k_{d,j}}\phi_{k_{d,k}}\right\rangle ,E[Qpc4]=i=1Pj=1Pk=1P=1Pqiqjqkqd=1Dϕkd,iϕkd,jϕkd,kϕkd,,E\left[Q^4_{\mathrm{pc}}\right] = \sum_{i=1}^{P}\sum_{j=1}^{P}\sum_{k=1}^{P}\sum_{\ell=1}^{P}q_iq_jq_kq_{\ell}\prod_{d=1}^{D}\left\langle\phi_{k_{d,i}}\phi_{k_{d,j}}\phi_{k_{d,k}}\phi_{k_{d,\ell}}\right\rangle ,

where the polynomial norms are computed using one-dimensional quadrature. We see here the number of operations required to compute Kurtosis is approximately DN4DN^{4}. If the number of coefficients is sufficiently high, these moments would probably be best computed inexactly by sampling the PCE surrogate.

Implementation

The PolynomialChaos user object takes in a list of distributions and constructs a polynomial class based on their type. Given a sampler and a vectorpostprocessor of results from sampling, it then loops through the MC or quadrature points to compute the coefficients. The statistical moments are then computed based on user preferences and the model can be evaluated using the evaluate function from any other moose object that has the reference. The algorithm uses the parallelization of the sampler to compute the coefficients, no other part of the algorithm is parallelized.

Example Input File Syntax

The example involves a homogeneous, one-dimensional diffusion-reaction problem, where the diffusion coefficient (DD) and reaction coefficient (Σ\Sigma) are uncertain with a uniform probability:

Dd2udx2+Σu=Q,-D\frac{d^2u}{dx^2} + \Sigma u = Q ,

where the QoI is the integrated average of u(x)u(x).

After the sub-app is set with the diffusion-reaction problem, distributions are set for each uncertain parameter:

[Distributions]
  [D_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
  [S_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

A sampler is then defined, either using Monte Carlo,

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 100
    distributions = 'D_dist S_dist'
    execute_on = initial
  []
[]

[Trainers]
  [poly_chaos]
    type = PolynomialChaosTrainer
    execute_on = timestep_end
    order = 5
    distributions = 'D_dist S_dist'
    sampler = sample
    response = storage/data:avg:value
    regression_type = integration
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_mc.i)

or quadrature,

[Samplers]
  [quadrature]
    type = Quadrature
    distributions = 'D_dist S_dist'
    execute_on = INITIAL
    order = 5
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

It is important that the order in the quadrature sampler input matches the order in the PolynomialChaos input. The sampler is then used by the MultiApp and Transfers to sample the sub-app, the QoI from the app is then put in a reporter:

[Reporters]
  [storage]
    type = StochasticReporter
  []
  [pc_samp]
    type = EvaluateSurrogate
    model = poly_chaos
    sampler = sample
    parallel_type = ROOT
    execute_on = final
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

All this information is ready to be sent to the PolynomialChaos trainer:

[Surrogates]
  [poly_chaos]
    type = PolynomialChaos
    trainer = poly_chaos
  []
[]

[Trainers]
  [poly_chaos]
    type = PolynomialChaosTrainer
    execute_on = timestep_end
    order = 5
    distributions = 'D_dist S_dist'
    sampler = quadrature
    response = storage/data:avg:value
  []
[]
(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_quad.i)

Input Parameters

  • filenameThe name of the file which will be associated with the saved/loaded data.

    C++ Type:FileName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the file which will be associated with the saved/loaded data.

  • trainerThe SurrogateTrainer object. If this is specified the trainer data is automatically gathered and available in this SurrogateModel object.

    C++ Type:UserObjectName

    Unit:(no unit assumed)

    Controllable:No

    Description:The SurrogateTrainer object. If this is specified the trainer data is automatically gathered and available in this SurrogateModel object.

Optional Parameters

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Set the enabled status of the MooseObject.

Advanced Parameters