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

ot.Log.Show(ot.Log.NONE)

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
  • relative errors: [4.8855e-05]
  • residuals: [0.0142158]
Index Multi-index Coeff.
0 [0,0,0,0] 26.12489
1 [1,0,0,0] -0.008183713
2 [0,1,0,0] 0.03242933
3 [0,0,1,0] 24.89812
4 [0,0,0,1] 3.977448
5 [2,0,0,0] -0.007980056
6 [1,1,0,0] 0.566605
7 [1,0,1,0] -0.01328741
8 [1,0,0,1] -0.01765562
9 [0,2,0,0] 0.003943191
10 [0,1,1,0] 0.02739853
11 [0,1,0,1] -0.0491181
12 [0,0,2,0] 9.098184
13 [0,0,1,1] -0.01695863
14 [0,0,0,2] 2.463476
15 [3,0,0,0] -0.01985544
16 [2,1,0,0] -0.007527326
17 [2,0,1,0] 0.02265436
18 [2,0,0,1] 0.02087004
19 [1,2,0,0] 0.02155156
20 [1,1,1,0] -0.01319882
21 [1,1,0,1] -0.01864534
22 [1,0,2,0] 0.0136746
23 [1,0,1,1] 0.02377915
24 [1,0,0,2] 0.05699986
25 [0,3,0,0] -0.006699248
26 [0,2,1,0] -0.02076579
27 [0,2,0,1] -0.0100742
28 [0,1,2,0] -0.02103417
29 [0,1,1,1] -0.008148804
30 [0,1,0,2] 0.0206051
31 [0,0,3,0] -0.01129148
32 [0,0,2,1] -0.0109149
33 [0,0,1,2] -0.003731515
34 [0,0,0,3] 0.946916


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.7186]

Print residuals.

result.getResiduals()
class=Point name=Unnamed dimension=1 values=[0.0142158]


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: 9
  • relative errors: [8.37134e-33]
  • residuals: [1.83709e-15]
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 [0,0,0,3] 0.9657677
7 [0,0,0,4] 0.2409653
8 [0,0,2,4] -1.030422e-14


Get the number of coefficients in the PCE.

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

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.