Gradient and Hessian for polynomial chaos metamodels

Python usage
1

Hi,

i have constructed a polynomial chaos metamodel of a function with input dimension dx and output dimension dy, and i would like to evaluate the gradient of this metamodel in a vectorized way, i.e. for an input array with size (n,dx), i would like to have the corresponding metamodel-based gradients in a single call, which would be an array with size (n, dy, dx). I also would like to have hessian evaluations in a similar way, i.e. to have a (n, dy, dx, dx) array from an input (n, dx) array.
Is there a way to do this in OT ?

I only found a way to do this in a “loopy” fashion, by using the getMetaModel() and then the gradient() method, which only takes a single Point object in argument (then, for the gradient, i could obtain the (n,dy,dx) output array by looping n times on the lines of the input (n,dx) array, by converting each line to a Point object, then computing the gradient with the gradient() method at this point, and so on… however, this is obviously very slow ! :smile: )

Best regards,
Fab

1 Like
2

Hi,
Unfortunately it is not possible to perform the computation in a vectorized way. I profiled the following code:

from openturns import *
from math import pi

# Problem parameters
dimension = 3
a = 7.0
b = 0.1
# Create the Ishigami function
inputVariables = ["xi1", "xi2", "xi3"]
formula = ["sin(xi1) / (1+xi2^2) + (" + str(a) + \
    ") * (sin(xi2)) ^ 2 + (" + str(
        b) + ") * xi3^4 * sin(xi1)"]
model = SymbolicFunction(inputVariables, formula)

# Create the input distribution
distribution = ComposedDistribution([Uniform(-pi, pi)] * dimension)

# Create the orthogonal basis
enumerateFunction = LinearEnumerateFunction(dimension)
productBasis = OrthogonalProductPolynomialFactory(
    [LegendreFactory()] * dimension, enumerateFunction)

# Create the adaptive strategy
degree = 20
adaptiveStrategy = FixedStrategy(productBasis, enumerateFunction.getStrataCumulatedCardinal(degree))

# Create the projection strategy
samplingSize = 10000
projectionStrategy = LeastSquaresStrategy(MonteCarloExperiment(samplingSize))
# Create the polynomial chaos algorithm
algo = FunctionalChaosAlgorithm(
    model, distribution, adaptiveStrategy, projectionStrategy)
print("run")
algo.run()
result = FunctionalChaosResult(algo.getResult())

metaModel = result.getMetaModel()
print("done")
from time import time
import numpy as np
benchSize = 100000
inSample = distribution.getSample(benchSize)
resGradient = np.zeros((benchSize, metaModel.getInputDimension(), metaModel.getOutputDimension()))
t0 = time()
for i in range(benchSize):
    resGradient[i,:,:] = metaModel.gradient(inSample[i, :])
t1 = time()
print("Loop, t=", t1 - t0, "s, speed=", benchSize / (t1 - t0), "grads/s")

and using FlameGraph I get:

perf_gradient
perf_gradient2835Ă—646 75.6 KB

which shows that the evaluation of the underlying univariate polynomials takes only 26% of the time spent in the evaluation of the gradient, which is bad! A lot of time is lost in the creation and destruction of many small objects (Point, Matrix). It should be possible to gain a significant factor here (like 3x speed-up), and if we can add parallelism the gain should scale linearly with the number of cores.
Unfortunately I am already working hard on other parts of the library. Don’t hesitate to open an issue on github, you will have a better chance to see some improvement in a near future this way!

Cheers
RĂ©gis

3

I worked a little bit on the implementation of OrthogonalUniVariatePolynomial and the net result is a 10% reduction on the benchmark I posted. Not that much but still good to take! See Improved OrthogonalUniVariatePolynomial by regislebrun · Pull Request #1961 · openturns/openturns · GitHub.