Chaos Vs Saltelli

Methodology
,
1

Hi,

I have performed a bunch of 100k simulation using a satellite reentry software and I’m trying to analyse its sensitivity to perturbations in the input parameters.

Using the PCE approach I get results that are consistent with the expectations, however I can only use smaller sample size (max 50k samples) on my PC as the computation time gets quickly prohibitive in time.

The sampling methods such as Saltelli or Martinez seems to be faster on my PC and use less resources. Also these methods provide confidence intervals. However they yield very inconsistent results exactly on the same data. All the indices are close to zero or negative.

With 100k simulations and a 5-dimensionnal input space I can use a size of 14000 which I though was enough.

The only difference with the example given on the OpenTURNS site ( Sobol’ indices for the Ishigami function) is that I do not use OT to generate the input sequence.

Does this detail matter ?

thank you for your help,

JM

2

Hi @jmp, you need to generate the input design with the OpenTURNS class SobolIndicesExperiment because both SaltelliSensitivityAlgorithm and MartinezSensitivityAlgorithm assume that you did that. To quote the doc pages of these 2 classes:

The class constructor Saltelli [or] MartinezSensitivityAlgorithm(inputDesign, outputDesign, N) requires an outputDesign produced by SobolIndicesExperiment (see example below). Otherwise, results will be worthless. [Actually, there is a typo in this sentence: it is the inputDesign that must be produced by SobolIndicesExperiment . The outputDesign of course contains the model output corresponding to the input design.]

You could, in theory, construct your inputDesign outside OpenTURNS, but you would need to follow the procedure detailed in the doc page of SobolIndicesExperiment. It is easier to let OpenTURNS handle that.

Basically, in order to estimate Sobol’ indices, you need two independent samples from the distribution of the input variables. Then you need to create additional samples based on these two by switching columns. The inputDesign expected by SaltelliSensitivityAlgorithm and MartinezSensitivityAlgorithm is the concatenation of all these samples.

3

Hi,
I am very surprised by the fact that you cannot use the full 100k sample with PCE in a reasonable amount of time. I suspect that you used the simplified interface, i.e. FunctionalChaosAlgorithm(inputSample, outputSample), and the algorithm is probably spending a huge amount of time trying to build the input distribution (which is done when the algorithm is created), not so much during the computation of the approximation.
You can check this by increasing the verbosity: add Log.Show(Log.ALL) just before the creation of the algorithm, and it should pause on one of the most expansive distribution factories.
A better way to use this algorithm is to go deeper into its elements. If you already know the input distribution there is no need to estimate it. If you know the smoothness of your function there is no need to go to high degrees. Here is a full script allowing to play with the different elements:

from openturns import *

# Log.Show(Log.ALL)
# Generate some data. In your case, you just have to load it from file
# in_sample = Sample.ImportFromCSVFile("in_sample.csv")
# in_dim = in_sample.getDimension()
# out_sample = Sample.ImportFromCSVFile("out_sample.csv")
N = 100000
in_dim = 5
d_ref = Normal(in_dim)
in_sample = d_ref.getSample(N)
f = SymbolicFunction(["x1","x2","x3","x4","x5"], ["x1*sin(x2+x3)+sin(x4+sin(x5))"])
out_sample = f(in_sample)
####################################
# PCE using the detailed interface #
####################################
# Let's assume we don't know d_ref. Otherwise, use the actual input distribution
marginals = [HistogramFactory().build(in_sample.getMarginal(i)) for i in range(in_dim)]
d = ComposedDistribution(marginals)
# Build the polynomial basis
polynomials = [StandardDistributionPolynomialFactory(m) for m in marginals]
basis = OrthogonalProductPolynomialFactory(polynomials)
# Build the projection strategy
projection = LeastSquaresStrategy(in_sample, out_sample)
# Build the adaptive strategy
max_degree = 6
total_size = basis.getEnumerateFunction().getStrataCumulatedCardinal(max_degree)
adaptive = FixedStrategy(basis, total_size)
# Create the PCE algo
from time import time
t0 = time()
algo = FunctionalChaosAlgorithm(in_sample, out_sample, d, adaptive, projection)
# If you want to get the exact same results using the simplified interface replace all the preceding lines by
# ResourceMap.SetAsUnsignedInteger("FunctionalChaosAlgorithm-MaximumTotalDegree", 6)
# ResourceMap.SetAsScalar("FunctionalChaosAlgorithm-QNorm", 1.0)
# algo = FunctionalChaosAlgorithm(in_sample, out_sample)
t1 = time()
print("Creation t=", t1 - t0, "s")
t0 = time()
algo.run()
t1 = time()
result = algo.getResult()
print("Run t=", t1 - t0, "s")
print(FunctionalChaosSobolIndices(result).summary())
##############
# Validation #
##############
in_sample_valid = d_ref.getSample(N)
out_sample_valid = f(in_sample_valid)
validation = MetaModelValidation(in_sample_valid, out_sample_valid, result.getMetaModel())
print("Q2=", validation.computePredictivityFactor())

Using the simplified interface I get the following results:

Creation t= 146.34749007225037 s
Run t= 0.2926785945892334 s
 input dimension: 5
 output dimension: 1
 basis size: 101
 mean: [-0.00257423]
 std-dev: [0.878122]
------------------------------------------------------------
Index   | Multi-indice                  | Part of variance  
------------------------------------------------------------
      4 | [0,0,0,1,0]                   | 0.302016
     17 | [1,0,1,0,0]                   | 0.176462
     16 | [1,1,0,0,0]                   | 0.170647
      5 | [0,0,0,0,1]                   | 0.135441
     51 | [0,0,0,2,1]                   | 0.0700867
     14 | [0,0,0,3,0]                   | 0.0496075
     70 | [1,3,0,0,0]                   | 0.0331452
     74 | [1,0,3,0,0]                   | 0.0287432
     15 | [0,0,0,0,3]                   | 0.0235836
------------------------------------------------------------


------------------------------------------------------------
Component | Sobol index            | Sobol total index      
------------------------------------------------------------
        0 | 0.000308946            | 0.409668              
        1 | 0.000163768            | 0.204437              
        2 | 0.000341013            | 0.205788              
        3 | 0.354638               | 0.428833              
        4 | 0.160817               | 0.235004              
------------------------------------------------------------


Q2= [0.792232]

and using the detailed interface:

Creation t= 0.13089418411254883 s
Run t= 0.9613080024719238 s
 input dimension: 5
 output dimension: 1
 basis size: 462
 mean: [-0.00169308]
 std-dev: [0.984487]
------------------------------------------------------------
Index   | Multi-indice                  | Part of variance  
------------------------------------------------------------
      4 | [0,0,0,1,0]                   | 0.241599
      7 | [1,1,0,0,0]                   | 0.141194
      8 | [1,0,1,0,0]                   | 0.140077
      5 | [0,0,0,0,1]                   | 0.107964
     75 | [1,1,2,0,0]                   | 0.0709469
     72 | [1,2,1,0,0]                   | 0.0691765
     53 | [0,0,0,2,1]                   | 0.0526255
     52 | [0,0,0,3,0]                   | 0.0403512
     81 | [1,0,3,0,0]                   | 0.0238355
     71 | [1,3,0,0,0]                   | 0.0218094
     55 | [0,0,0,0,3]                   | 0.0184457
    326 | [1,3,2,0,0]                   | 0.0127656
    332 | [1,2,3,0,0]                   | 0.0125651
    249 | [0,0,0,2,3]                   | 0.0115486
------------------------------------------------------------


------------------------------------------------------------
Component | Sobol index            | Sobol total index      
------------------------------------------------------------
        0 | 4.75502e-06            | 0.508491              
        1 | 8.76558e-06            | 0.343322              
        2 | 1.65735e-05            | 0.344294              
        3 | 0.284018               | 0.363659              
        4 | 0.127906               | 0.207842              
------------------------------------------------------------


Q2= [0.990007]

You can see the significant difference in terms of CPU time (147s vs 1s, on a quite powerful workstation) and the improved quality on this specific example (but you can get the exact same results by tweaking the behavior of PCE using ResourceMap, see the script).

Cheers

Régis

4

Hi @josephmure
thank you for your explanations, I think I do not feed the algorithm with the correct data.
I can always split my data into 2 independant sets but I do not really understand how to get through the right steps in order to achieve the sensitivity analysis using Saltelli algorithm. It could be I’m missing some theory…

Regards,
JM

5

Hi @regislebrun,

OK I get thet point.
Indeed when providing the data only, the algorithm attempts to estimate the distributions and it takes a while for a big sample size.
I’ll try to provide the distribution to the algo it should speed up things.

Thanks,

JM

6

Hi @regislebrun,

Thank you for the tips you provided, I could indeed provide the prefined distribution to the PCE algorithm and benefit from some speed up. The speed is rather moderate on my machine, I have the feeling that the most time consuming phase was not the estimation of the distribution but rather the optimisation phase. It often ends with “Terminating: Residual and Jacobian evaluation failed” probably meaning that the convergence for the Jacobian evaluation criteria was not met (that’s only an assumption)

I’m currently using the 100000 x 5 dataset and unfortunately perfoming the same evaluation on a 100000 x 28 dataset takes much much longer (I already tried with 9), even on my 16 cores machine with 32GB ram.

I was wondering whether we could take benefit from the hardware architecture of the machine using e.g. process parallelisation ?

Thank you for your help,
Cheers,
JM

7

Hi @jmp,
The message you get is coming from the Ceres solver, which is strange as no optimization algorithm is involved in the FunctionalChaosAlgorithm class. To be honest I don’t understand what you mean when talking about an “optimization phase”. If you use the sparsity inducing technics (LARS(), KFold()) then there is an exploratory phase in the algorithm, which can take a significant amount of time, but cannot produce the message you mention as it is not based at all on optimization solvers.

I adapted the script I already published to check the memory consumption in the case of 28 input variables and an input sample of size 1e5. Depending on how you parametrized your polynomial basis, you may have too much polynomials to consider. For example, if you keep all the polynomials of total degree less or equal to 4, you already have 35960 polynomials, resulting into a design matrix of 27Go, huge. It will take a significant amount of time to solve the resulting least-squares problem, and you will probably have to swap data on disk, which is terrible for the performance. If your installation of OT is based on parallel blas (it is the case if you installed OT using conda) you will fully benefit from your 32 cores.

What you can do to improve things:

  • compute a first meta-model using a moderate total degree (2 or 3), and use it to compute Sobol total indices. Based on these indices, reduce the input dimension by selecting the input variables with the largest total indices
  • if you don’t want to throw out input variables, use these indices to weight the univariate polynomials in your basis:
    enum_function=HyperbolicAnisotropicEnumerateFunction([1.0/s for s in sobolTotalIndices])
    or another formula giving more penalization weight to the variables of low total indices. Then you use this function to create a new polynomial basis:
    basis = OrthogonalProductPolynomialFactory(polynomials, enum_function)
  • re-run the algorithm using a maximum number of functions in the basis instead of a maximum total degree

You can also play with the sparsity-inducing selection of functions:
projection = LeastSquaresStrategy(LeastSquaresMetaModelSelectionFactory(LARS(), CorrectedLeaveOneOut()))
but it can take significantly more time.

Of course it is a kind of chicken and egg problem: you select variables based on a poor approximation of the Sobol indices, then you refine the meta-model in order to get a better approximation of the Sobol indices of the remaining variables. The key point is to check the quality of the last meta-model using its Q_2 predictivity factor. If it is close to 1 this approach is valid, as you kept most of the variance of the output so its decomposition through Sobol indices is valid. Otherwise you are in trouble. I have other ideas but it would be better to first try what I proposed.

Cheers

Régis