Gaussian Process Regression: choose an arbitrary trend

The goal of this example is to show how to configure an arbitrary trend in a Kriging metamodel. In the Gaussian Process Regression: choose a polynomial trend and Gaussian Process Regression: choose a polynomial trend on the beam model examples, we show how to configure a polynomial trend.

In general, any collection of multivariate functions can be used as the basis argument of a GaussianProcessFitter or a GaussianProcessRegression. In practice, it might not be convenient to create a multivariate basis and this is why we sometimes create it by tensorization of univariate functions. In this example, we first use Legendre polynomials as our univariate functions, then we create an orthogonal polynomial basis corresponding to the input marginals.

For this purpose, we use the cantilever beam example.

from openturns.usecases import cantilever_beam
import openturns as ot
import openturns.experimental as otexp

Definition of the model

We load the cantilever beam use case

cb = cantilever_beam.CantileverBeam()

We load the function (model) which evaluates the output Y depending on the inputs.

model = cb.model

Then we define the distribution of the input random vector.

dimension = cb.dim  # number of inputs
input_dist = cb.distribution

Create the design of experiments

We consider a simple Monte-Carlo sampling as a design of experiments. This is why we generate an input sample using the getSample method of the distribution. Then we evaluate the output using the model function.

sampleSize_train = 20
X_train = input_dist.getSample(sampleSize_train)
Y_train = model(X_train)

Create the Legendre basis

We first create a Legendre basis of univariate polynomials.

The LegendreFactory class creates Legendre polynomials.

univariateFactory = ot.LegendreFactory()

These polynomials are orthogonal with respect to the Uniform distribution on , as we can see.

univariateFactory.getMeasure()

Uniform

• name=Uniform
• dimension=1
• weight=1
• range=[-1, 1]
• description=[X0]
• isParallel=true
• isCopula=false

Even if the measure Uniform([-1, 1]) does not correspond to the input marginal distributions, these polynomials will, anyway, create a consistent trend for the Gaussian process regression metamodel.

We use the same polynomial family for each input. The multivariate polynomials are created by tensorization of the univariate polynomials. The linear enumerate function specifies the order followed to create the multivariate polynomials.

polyColl = [univariateFactory] * dimension
enumerateFunction = ot.LinearEnumerateFunction(dimension)
productBasis = ot.OrthogonalProductPolynomialFactory(polyColl, enumerateFunction)

The basis contains the first multivariate polynomials ordered according to the linear enumerate function.

functions = []
numberOfTrendCoefficients = 12
for i in range(numberOfTrendCoefficients):
    multivariatepolynomial = productBasis.build(i)
    print(multivariatepolynomial)
    functions.append(multivariatepolynomial)

1
1.73205 * x0
1.73205 * x1
1.73205 * x2
1.73205 * x3
-1.11803 + 3.3541 * x0^2
1.73205 * x0 * 1.73205 * x1
1.73205 * x0 * 1.73205 * x2
1.73205 * x0 * 1.73205 * x3
-1.11803 + 3.3541 * x1^2
1.73205 * x1 * 1.73205 * x2
1.73205 * x1 * 1.73205 * x3

basis = ot.Basis(functions)

Create the Gaussian Process Regression metamodel

In order to create the Gaussian process regression metamodel, we choose the function basis created previously and a squared exponential covariance model.

covariance_model = ot.SquaredExponential([1.0] * dimension, [1.0])

First, we estimate a Gaussian process approximating the model with the class GaussianProcessFitter.

algo_fit = otexp.GaussianProcessFitter(X_train, Y_train, covariance_model, basis)
print("First run: algo GPFitter = ", algo_fit.getOptimizationAlgorithm())
algo_fit.setOptimizationAlgorithm(ot.TNC())
algo_fit.run()
fitter_result = algo_fit.getResult()

First run: algo GPFitter =  class=OptimizationAlgorithm implementation=class=Cobyla class=OptimizationAlgorithmImplementation problem=class=OptimizationProblem implementation=class=OptimizationProblemImplementation objective=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[] evaluationImplementation=class=NoEvaluation name=Unnamed gradientImplementation=class=NoGradient name=Unnamed hessianImplementation=class=NoHessian name=Unnamed equality constraint=none inequality constraint=none bounds=none minimization=true dimension=0 startingPoint=class=Point name=Unnamed dimension=0 values=[] maximumIterationNumber=100 maximumCallsNumber=1000 maximumAbsoluteError=1e-05 maximumRelativeError=1e-05 maximumResidualError=1e-05 maximumConstraintError=1e-05 rhoBeg=0.1

Then, we condition the estimated gaussian process to make it interpolate the data set using the class GaussianProcessRegression.

gpr_algo = otexp.GaussianProcessRegression(fitter_result)
gpr_algo.run()
gpr_result = gpr_algo.getResult()

We get the Gaussian process regression metamodel.

gpr_metamodel = gpr_result.getMetaModel()

The method getTrendCoefficients() returns the coefficients of the trend. We see that the number of coefficients in the trend corresponds to the number of functions in the basis.

print(gpr_result.getTrendCoefficients())

[7.28015e-46,8.51624e-35,3.75184e-43,3.21875e-45,1.81292e-52,1.11494e-23,4.38873e-32,3.76531e-34,2.12187e-41,2.17629e-40,1.65929e-42,9.3468e-50]#12

We get the updated covariance model. The SquaredExponential model has one amplitude coefficient and 4 scale coefficients. This is because this covariance model is anisotropic : each of the 4 input variables is associated with its own scale coefficient.

print(gpr_result.getCovarianceModel())

SquaredExponential(scale=[1,1,0.17976,0.01], amplitude=[0.0407412])

Create an orthogonal multivariate polynomial factory

We suggest here to create a multivariate polynomial basis where each marginal polynomial family corresponds to the polynomial family orthogonal to the marginal input distribution. We use the class OrthogonalProductPolynomialFactory created from the input marginal distributions. This corresponds to the basis we usually use in the polynomial chaos expansion.

We first create the mutlivariate polynomial basis as a tensorized product of the univariate polynomial basis orthogonal to the marginal input distributions.

multivariateBasis = ot.OrthogonalProductPolynomialFactory([cb.E, cb.F, cb.L, cb.II])

Then we create the multivariate basis which has maximum degree equal to 2. There are 15 functions contained in the basis.

totalDegree = 2
enumerateFunction = multivariateBasis.getEnumerateFunction()
numberOfTrendCoefficients = enumerateFunction.getStrataCumulatedCardinal(totalDegree)
print("Number of functions in the basis: ", numberOfTrendCoefficients)

Number of functions in the basis:  15

functions = []
for i in range(numberOfTrendCoefficients):
    multivariatepolynomial = productBasis.build(i)
    print(multivariatepolynomial)
    functions.append(multivariatepolynomial)

1
1.73205 * x0
1.73205 * x1
1.73205 * x2
1.73205 * x3
-1.11803 + 3.3541 * x0^2
1.73205 * x0 * 1.73205 * x1
1.73205 * x0 * 1.73205 * x2
1.73205 * x0 * 1.73205 * x3
-1.11803 + 3.3541 * x1^2
1.73205 * x1 * 1.73205 * x2
1.73205 * x1 * 1.73205 * x3
-1.11803 + 3.3541 * x2^2
1.73205 * x2 * 1.73205 * x3
-1.11803 + 3.3541 * x3^2

basis = ot.Basis(functions)

algo_fit = otexp.GaussianProcessFitter(X_train, Y_train, covariance_model, basis)
print("Second run: algo GPFitter = ", algo_fit.getOptimizationAlgorithm())
algo_fit.setOptimizationAlgorithm(ot.TNC())
algo_fit.run()
fitter_result = algo_fit.getResult()
gpr_algo = otexp.GaussianProcessRegression(fitter_result)
gpr_algo.run()
gpr_result = gpr_algo.getResult()
print(gpr_result.getTrendCoefficients())

Second run: algo GPFitter =  class=OptimizationAlgorithm implementation=class=Cobyla class=OptimizationAlgorithmImplementation problem=class=OptimizationProblem implementation=class=OptimizationProblemImplementation objective=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[] evaluationImplementation=class=NoEvaluation name=Unnamed gradientImplementation=class=NoGradient name=Unnamed hessianImplementation=class=NoHessian name=Unnamed equality constraint=none inequality constraint=none bounds=none minimization=true dimension=0 startingPoint=class=Point name=Unnamed dimension=0 values=[] maximumIterationNumber=100 maximumCallsNumber=1000 maximumAbsoluteError=1e-05 maximumRelativeError=1e-05 maximumResidualError=1e-05 maximumConstraintError=1e-05 rhoBeg=0.1
[7.28015e-46,8.51624e-35,3.75184e-43,3.21875e-45,1.81292e-52,1.11494e-23,4.38873e-32,3.76531e-34,2.12187e-41,2.17629e-40,1.65929e-42,9.3468e-50,1.50988e-44,8.01436e-52,-8.13946e-46]#15

Conclusion

The trend that we have configured corresponds to the basis that we would have used in a full polynomial chaos expansioncomputed with least squares.

Other extensions of this work would be:

• to use a Fourier basis with FourierSeriesFactory,

• wavelets with HaarWaveletFactory,

or any other univariate factory.