Create a full or sparse polynomial chaos expansion

In this example we create a global approximation of a model using polynomial chaos expansion based on a design of experiments. The goal of this example is to show how we can create a full or sparse polynomial chaos expansion depending on our needs and depending on the number of observations we have. In general, we should have more observations than parameters to estimate. This is why a sparse polynomial chaos may be interesting: by carefully selecting the coefficients we estimate, we may reduce overfitting and increase the predictions of the metamodel.

import openturns as ot

Define the model

Create the function.

myModel = ot.SymbolicFunction(
    ["x1", "x2", "x3", "x4"], ["1 + x1 * x2 + 2 * x3^2 + x4^4"]
)

Create a multivariate distribution.

distribution = ot.JointDistribution(
    [ot.Normal(), ot.Uniform(), ot.Gamma(2.75, 1.0), ot.Beta(2.5, 1.0, -1.0, 2.0)]
)

In order to create the PCE, we can specify the distribution of the input parameters. If not known, statistical inference can be used to select a possible candidate, and fitting tests can validate such an hypothesis. Please read Fit a distribution from an input sample for an example of this method.

Create a training sample

Create a pair of input and output samples.

sampleSize = 250
inputSample = distribution.getSample(sampleSize)
outputSample = myModel(inputSample)

Build the orthogonal basis

In the next cell, we create the univariate orthogonal polynomial basis for each marginal.

inputDimension = inputSample.getDimension()
coll = [
    ot.StandardDistributionPolynomialFactory(distribution.getMarginal(i))
    for i in range(inputDimension)
]
enumerateFunction = ot.LinearEnumerateFunction(inputDimension)
productBasis = ot.OrthogonalProductPolynomialFactory(coll, enumerateFunction)

We can achieve the same result using OrthogonalProductPolynomialFactory.

marginalDistributionCollection = [
    distribution.getMarginal(i) for i in range(inputDimension)
]
multivariateBasis = ot.OrthogonalProductPolynomialFactory(
    marginalDistributionCollection
)
multivariateBasis
  • dimension: 4
  • enumerate function: class=LinearEnumerateFunction dimension=4
Index Name Distribution Univariate polynomial
0 X0 Normal HermiteFactory
1 X1 Uniform LegendreFactory
2 X2 Gamma LaguerreFactory
3 X3 Beta JacobiFactory


Create a full PCE

Create the algorithm. We compute the basis size from the total degree. The next lines use the LeastSquaresStrategy class with default parameters (the default is the PenalizedLeastSquaresAlgorithmFactory class). This creates a full polynomial chaos expansion, i.e. we keep all the candidate coefficients produced by the enumeration rule. In order to create a sparse polynomial chaos expansion, we must use the LeastSquaresMetaModelSelectionFactory class instead.

totalDegree = 3
candidateBasisSize = enumerateFunction.getBasisSizeFromTotalDegree(totalDegree)
print("Candidate basis size = ", candidateBasisSize)
adaptiveStrategy = ot.FixedStrategy(productBasis, candidateBasisSize)
projectionStrategy = ot.LeastSquaresStrategy()
algo = ot.FunctionalChaosAlgorithm(
    inputSample, outputSample, distribution, adaptiveStrategy, projectionStrategy
)
algo.run()
result = algo.getResult()
result
Candidate basis size =  35
FunctionalChaosResult
  • input dimension: 4
  • output dimension: 1
  • distribution dimension: 4
  • transformation: 4 -> 4
  • inverse transformation: 4 -> 4
  • orthogonal basis dimension: 4
  • indices size: 35
Index Multi-index Coeff.
0 [0,0,0,0] 26.12842
1 [1,0,0,0] -0.01980392
2 [0,1,0,0] 0.02165394
3 [0,0,1,0] 24.87216
4 [0,0,0,1] 3.993464
5 [2,0,0,0] -0.01137476
6 [1,1,0,0] 0.5769773
7 [1,0,1,0] -0.00666321
8 [1,0,0,1] -0.04540182
9 [0,2,0,0] 0.002660498
10 [0,1,1,0] 0.02179118
11 [0,1,0,1] -0.05160322
12 [0,0,2,0] 9.070564
13 [0,0,1,1] 0.03938575
14 [0,0,0,2] 2.426618
15 [3,0,0,0] -0.001311593
16 [2,1,0,0] 0.01813314
17 [2,0,1,0] -0.001881889
18 [2,0,0,1] 0.02305903
19 [1,2,0,0] 0.01089606
20 [1,1,1,0] 0.01388923
21 [1,1,0,1] -0.03398956
22 [1,0,2,0] -0.004318601
23 [1,0,1,1] 0.0480162
24 [1,0,0,2] 0.01735828
25 [0,3,0,0] -0.00550616
26 [0,2,1,0] 0.02198017
27 [0,2,0,1] 0.0003605928
28 [0,1,2,0] 0.05869172
29 [0,1,1,1] -0.004482646
30 [0,1,0,2] 0.05797735
31 [0,0,3,0] -0.01249474
32 [0,0,2,1] 0.01957749
33 [0,0,1,2] -0.001761463
34 [0,0,0,3] 0.9833603


Get the number of coefficients in the PCE.

selectedBasisSizeFull = result.getIndices().getSize()
print("Selected basis size = ", selectedBasisSizeFull)
Selected basis size =  35

We see that the number of coefficients in the selected basis is equal to the number of coefficients in the candidate basis. This is, indeed, a full PCE.

Use the PCE

Get the metamodel function.

metamodel = result.getMetaModel()

In order to evaluate the metamodel on a single point, we just use it as any other openturns.Function.

xPoint = distribution.getMean()
yPoint = metamodel(xPoint)
print("Value at ", xPoint, " is ", yPoint)
Value at  [0,0,2.75,1.14286]  is  [17.7498]

Based on these results, we may want to validate our metamodel. More details on this topic are presented in Validate a polynomial chaos.

Create a sparse PCE

In order to create a sparse polynomial chaos expansion, we use the LeastSquaresMetaModelSelectionFactory class instead.

totalDegree = 6
candidateBasisSize = enumerateFunction.getBasisSizeFromTotalDegree(totalDegree)
print("Candidate basis size = ", candidateBasisSize)
adaptiveStrategy = ot.FixedStrategy(productBasis, candidateBasisSize)
selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory()
projectionStrategy = ot.LeastSquaresStrategy(selectionAlgorithm)
algo = ot.FunctionalChaosAlgorithm(
    inputSample, outputSample, distribution, adaptiveStrategy, projectionStrategy
)
algo.run()
result = algo.getResult()
result
Candidate basis size =  210
FunctionalChaosResult
  • input dimension: 4
  • output dimension: 1
  • distribution dimension: 4
  • transformation: 4 -> 4
  • inverse transformation: 4 -> 4
  • orthogonal basis dimension: 4
  • indices size: 10
Index Multi-index Coeff.
0 [0,0,0,0] 26.11651
1 [0,0,1,0] 24.87469
2 [0,0,0,1] 3.974438
3 [1,1,0,0] 0.5773503
4 [0,0,2,0] 9.082951
5 [0,0,0,2] 2.419672
6 [2,0,1,0] 1.724218e-15
7 [0,0,0,3] 0.9657677
8 [1,0,0,3] -2.727748e-15
9 [0,0,0,4] 0.2409653


Get the number of coefficients in the PCE.

selectedBasisSizeSparse = result.getIndices().getSize()
print("Selected basis size = ", selectedBasisSizeSparse)
Selected basis size =  10

We see that the number of selected coefficients is lower than the number of candidate coefficients. This may reduce overfitting and can produce a PCE with more accurate predictions.