The hyperbolic enumeration function can be slow when q = 1

Python usage
1

Hi!
I created a polynomial chaos expansion of a function in relatively large dimension using regression and the LARS selection method. This may require to use the hyperbolic enumeration function. It appears, however, that this can be much slower when the quasi-norm parameter is equal to 1, even in the particular case where the weights of the “anisotropic” part are all equal. This is surprising, since the hyperbolic (isotropic) enumeration rule is mathematically equivalent to the linear enumeration.

The linear enumeration rule for a total polynomial degree equal to d in dimension p is:

\left\{\boldsymbol{\alpha} \in \mathbb{N}^p \; \textrm{s.t.} \; \sum_{i = 1}^p \alpha_i \leq d\right\}

The hyperbolic (isotropic) enumeration rule with quasi-norm parameter q for a total polynomial degree equal to d in dimension p is:

\left\{\boldsymbol{\alpha} \in \mathbb{N}^p \; \textrm{s.t.} \; \left( \sum_{i = 1}^p \alpha_i^q\right)^{\frac{1}{q}} \leq d\right\}

So, if q = 1, both rules are mathematically equivalent. I understand that the implementation may not reflect this, because the enumeration algorithm have to take into account that the parameter q can be different from 1. Hence, producing the multi-indices can be different, even if q = 1. What I do not understand, however, is the amount of the loss of speed.

In the next script, I create a function with input dimension 20. Then I create a sample of this function with size 200. I split that sample into a training set (75%) and a test set (25%). The key part is where I select either the linear enumeration rule or the hyperbolic enumeration rule with quasi-norm parameter equal to 1.

import openturns as ot
import time

t1 = time.time()

# Parameters
totalDegree = 3
useLinear = False  # Set to True to use linear rule

# Define the input distribution
inputDimension = 20
distributionList = [ot.Uniform(0, 1)] * inputDimension
myDistribution = ot.ComposedDistribution(distributionList)

# Create a sample
inputLabels = ["x%d" % (i) for i in range(inputDimension)]
g = ot.SymbolicFunction(inputLabels, ["1 + x1 + 2 * x3 + 3 * x5 - 5 * x13"])
inputTrain = myDistribution.getSample(100)
outputTrain = g(inputTrain)

# Split
split_fraction = 0.75
training_sample_size = inputTrain.getSize()
split_index = int(split_fraction * training_sample_size)
inputTest = inputTrain.split(split_index)
outputTest = outputTrain.split(split_index)

# Create enumeration function
if useLinear:
    # Trick: this is much faster than a Hyperbolic rule in this case
    enumerateFunction = ot.LinearEnumerateFunction(inputDimension)
else:
    quasiNorm = 1.0
    enumerateFunction = ot.HyperbolicAnisotropicEnumerateFunction(
        inputDimension, quasiNorm)

# Create the polynomial collection and the multivariate basis

polyColl = [
    ot.StandardDistributionPolynomialFactory(distributionList[i]) for i in range(inputDimension)
]

multivariateBasis = ot.OrthogonalProductPolynomialFactory(
    polyColl, enumerateFunction
)
print("Run!")
selectionAlgorithm = ot.LeastSquaresMetaModelSelectionFactory()
projectionStrategy = ot.LeastSquaresStrategy(
    inputTrain, outputTrain, selectionAlgorithm
)
enumerateFunction = multivariateBasis.getEnumerateFunction()
basisSize = enumerateFunction.getBasisSizeFromTotalDegree(totalDegree)
print("Basis size = ", basisSize)
adaptiveStrategy = ot.FixedStrategy(multivariateBasis, basisSize)
chaosAlgorithm = ot.FunctionalChaosAlgorithm(
    inputTrain, outputTrain, myDistribution, adaptiveStrategy, projectionStrategy
)
chaosAlgorithm.run()
chaosResult = chaosAlgorithm.getResult()
# Get number of coefficients in the selection
indices = chaosResult.getIndices()
number_of_coefficients = indices.getSize()
print("Selected number of coefficients =", number_of_coefficients)
coefficients = chaosResult.getCoefficients()
for i in range(number_of_coefficients):
    index_in_the_basis = indices[i]
    multiindex = enumerateFunction(index_in_the_basis)
    print("#", i, ", index = ", index_in_the_basis, " : ", 
          multiindex, ", coeff = ", coefficients[i, 0])

# Test
metamodel = chaosResult.getMetaModel()
val = ot.MetaModelValidation(inputTest, outputTest, metamodel)
q2score = val.computePredictivityFactor()[0]
print("Q2 = %.2f%%" % (100 * q2score))

t2 = time.time()
elapsed = t2 - t1
print("elapsed = %.1f (s)" % (elapsed))

Here is the output of this code.

  • With linear: Basis size = 1771, Selected number of coefficients = 6, Q2 = 100.00%, elapsed = 0.1 (s) 
# 0 , index =  0  :  [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  1.5000000000000007
# 1 , index =  2  :  [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.288675134594813
# 2 , index =  4  :  [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.577350269189626
# 3 , index =  6  :  [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.8660254037844387
# 4 , index =  14  :  [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0]#20 , coeff =  -1.443375672974065
# 5 , index =  1443  :  [0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0]#20 , coeff =  -3.5016511832433925e-16
  • With hyperbolic and q = 1: Basis size = 1771, Selected number of coefficients = 5, Q2 = 100.00%, elapsed = 26.6 (s) 
# 0 , index =  0  :  [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  1.5
# 1 , index =  2  :  [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.2886751345948131
# 2 , index =  4  :  [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.5773502691896257
# 3 , index =  6  :  [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#20 , coeff =  0.866025403784439
# 4 , index =  14  :  [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0]#20 , coeff =  -1.443375672974065

We see that the linear enumeration is much faster, although it produces the same output. I noticed that the same trick is used in the FunctionalChaosAlgorithm source code in the case where the enumeration function is not given by the user. So this problem seems to be known to some developers. I closely inspected the multi-indices produced by both rules: they are the same and the layers in which they appear is the same. Moreover, I benchmarked the speed of the production of the multi-indices.

image
image731×378 32.8 KB

The hyperbolic enumeration rule is slightly slower, but that does not explain that drastic speed difference: notice that the different in speed is seen for basis sizes larger than 10^4. In my case,
the basis size is only 1771 and the calculation of the multi-indices is (almost) instantaneous. For example, producing 1000 multi-indices in dimension 20 requires 0.4 seconds with the hyperbolic enumeration rule with q = 1.

Why is there such a difference in speed between the hyperbolic (with q= 1) and linear enumeration rule when used in the contexte of the FunctionalChaosAlgorithm class?

Regards,
Michaël