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 experiment. 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: [3.25754e-05]
  • residuals: [0.010775]
Index Multi-index Coeff.
0 [0,0,0,0] 26.10131
1 [1,0,0,0] 0.01833501
2 [0,1,0,0] 0.02643858
3 [0,0,1,0] 24.89616
4 [0,0,0,1] 3.998141
5 [2,0,0,0] -0.007498848
6 [1,1,0,0] 0.5836605
7 [1,0,1,0] 0.0008685963
8 [1,0,0,1] -0.01786958
9 [0,2,0,0] 0.005226953
10 [0,1,1,0] 0.02195566
11 [0,1,0,1] -0.06375109
12 [0,0,2,0] 9.063002
13 [0,0,1,1] -0.02439593
14 [0,0,0,2] 2.366967
15 [3,0,0,0] 0.003033619
16 [2,1,0,0] -0.01237667
17 [2,0,1,0] -0.01342518
18 [2,0,0,1] 0.01335615
19 [1,2,0,0] -0.009671841
20 [1,1,1,0] 0.004816446
21 [1,1,0,1] -0.02295916
22 [1,0,2,0] -0.01663346
23 [1,0,1,1] 0.007606394
24 [1,0,0,2] 0.03900244
25 [0,3,0,0] 0.02087256
26 [0,2,1,0] -0.01402968
27 [0,2,0,1] -0.04129633
28 [0,1,2,0] -0.01182831
29 [0,1,1,1] -0.01354742
30 [0,1,0,2] 0.09136447
31 [0,0,3,0] 0.02163227
32 [0,0,2,1] 0.02992845
33 [0,0,1,2] 0.04678639
34 [0,0,0,3] 1.057842


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

Print residuals.

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


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: 24
  • relative errors: [3.16988e-07]
  • residuals: [0.0102016]
Index Multi-index Coeff.
0 [0,0,0,0] 26.07564
1 [1,0,0,0] -0.003290232
2 [0,0,1,0] 24.86277
3 [0,0,0,1] 4.00878
4 [1,1,0,0] 0.5800043
5 [1,0,1,0] 0.005553137
6 [0,1,0,1] -0.03536237
7 [0,0,2,0] 9.006102
8 [0,0,0,2] 2.32328
9 [2,0,1,0] -0.01830508
10 [1,1,1,0] -0.001079013
11 [1,0,2,0] -0.01280951
12 [0,1,0,2] 0.07909247
13 [0,0,0,3] 1.0497
14 [0,2,0,2] 0.03390499
15 [0,0,2,2] -0.09301932
16 [2,3,0,0] 0.006792632
17 [2,2,1,0] -0.001824762
18 [1,0,1,3] -0.00731192
19 [0,1,0,4] 0.04903297
20 [0,4,2,0] 0.01100148
21 [0,2,2,2] 0.04061781
22 [0,1,1,4] -0.01867403
23 [0,0,2,4] -0.09888105


Get the number of coefficients in the PCE.

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

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.