Relatively slow evaluation of moments and sensitivities

[utsekaj42]

I am wondering if there are any quick improvements to make the evaluation of cp.Std, and cp.Sens_m, cp.Sens_t more performant.

In principle once the expansion coefficients are determined, these should be computed from linear combinations of these coefficients, thus it is surprising that he bulk of evaluation time is used for these computations and not for cp.orth_ttr or cp.fit_regression

I am using this script to test things

"""basic_chaospy.py"""
import chaospy as cp
import numpy as np

@profile
def test():
    dim = 8
    dist = cp.Iid(cp.Uniform(), dim)
    order = 4
    poly = cp.orth_ttr(order, dist)
    Ns = len(poly)*2
    print(Ns)
    samples = dist.sample(Ns)
    y_sample =  np.sum(samples, axis=0)

    approx = cp.fit_regression(poly, samples, y_sample)
    exp_pc = cp.E(approx, dist)
    std_pc = cp.Std(approx, dist)
    S_pc = cp.Sens_m(approx, dist) 
    S_t_pc = cp.Sens_m(approx, dist)

if __name__ == "__main__":
    test()

Some results , using line_profiler are

kernprof -v -l basic_chaospy.py

168
Wrote profile results to basic_chaospy.py.lprof
Timer unit: 1e-06 s

Total time: 110.117 s
File: basic_chaospy.py
Function: test at line 4

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     4                                           @profile
     5                                           def test():
     6         1          2.0      2.0      0.0      dim = 6
     7         1         32.0     32.0      0.0      dist = cp.Iid(cp.Uniform(), dim)
     8         1          1.0      1.0      0.0      order = 3
     9         1      33499.0  33499.0      0.0      poly = cp.orth_ttr(order, dist)
    10         1          3.0      3.0      0.0      Ns = len(poly)*2
    11         1         22.0     22.0      0.0      print(Ns)
    12         1       1513.0   1513.0      0.0      samples = dist.sample(Ns)
    13         1         16.0     16.0      0.0      y_sample =  np.sum(samples, axis=0)
    14
    15         1      36276.0  36276.0      0.0      approx = cp.fit_regression(poly, samples, y_sam
ple)
    16         1     136131.0 136131.0      0.1      exp_pc = cp.E(approx, dist)
    17         1   20470060.0 20470060.0     18.6      std_pc = cp.Std(approx, dist)
    18         1   20871824.0 20871824.0     19.0      S_pc = cp.Sens_m(approx, dist) #Si from chaos
py
    19         1   68567731.0 68567731.0     62.3      S_t_pc = cp.Sens_t(approx, dist) #Si from cha
ospy

I don't include results currently using higher order approximants or more dimensions, but the computational (and memory) demands for these estimates are prohibitive in some cases when seeking to verify convergence of the PCE.

[jonathf]

Yes, you are right.

The issue has been on my radar a while now. In fact I am working on switching the backend for the polynomials these days, which hopefully will fix the problem. But it will take a little time to get all the pieces together for this refactor. I can keep you posted if you like.

[utsekaj42]

Staying posted sounds great! If there are any specific portions of refactoring that could be delegated, I might be able to put some time trying to implement some improvements, but I suppose if the solution is just to replace the polynomial backend it may be hard to delegate.

[utsekaj42]

I have implemented a workaround for my use case, but perhaps this can be integrated in a more general manner. Currently, this only works if one uses normed polynomials for the expansion, and then use retall in the fitting method to return the coefficients of the PCE as uhat:

import chaospy as cp
import numpy as np

@profile
def test():
    dim = 5
    dist = cp.Iid(cp.Uniform(), dim)
    order = 2
    poly = cp.orth_ttr(order, dist, normed=True)
    Ns = len(poly)*2
    print(Ns)
    samples = dist.sample(Ns)
    y_sample =  np.sum(samples, axis=0)

    approx, uhat = cp.fit_regression(poly, samples, y_sample, retall=1)
    exp_pc = cp.E(approx, dist)
    print(exp_pc, uhat[0])
    std_pc = cp.Std(approx, dist)
    variance = np.sum(uhat[1:]**2)
    std_coefs = np.sqrt(np.sum(uhat[1:]**2))
    print(std_pc, std_coefs)
    S_pc = cp.Sens_m(approx, dist) #Si from chaospy
    S_t_pc = cp.Sens_t(approx, dist) #Si from chaospy
    main, total = calc_sens(dim, poly, uhat)
    print(S_pc, main)
    print(S_t_pc, total)

def calc_sens(dim, poly, uhat):
    variance = np.sum(uhat[1:]**2)
    main = np.zeros(dim)
    total = np.zeros(dim)
    #Compute sens
    for idx, pol in enumerate(poly[1:]):
        add_total = [False,]*dim
        add_main = [True,]*dim
        for term in pol.A:
            for var in range(dim):
                if term[var] > 0:
                    add_total[var] = True
                    add_main[var] = add_main[var] and True
                    add_main[0:var] = [False,]*var
                    try:
                        add_main[var+1::] = [False,]*(dim-var)
                    except IndexError:
                        pass
        for var in range(dim):
            if add_main[var]:
                main[var] += uhat[idx+1]**2
            if add_total[var]:
                total[var] += uhat[idx+1]**2

    main = main/variance
    total = total/variance
    return main, total

if __name__ == "__main__":
    test()

which is much faster:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     4                                           @profile
     5                                           def test():
     6         1          2.0      2.0      0.0      dim = 5
     7         1         33.0     33.0      0.0      dist = cp.Iid(cp.Uniform(), dim)
     8         1          0.0      0.0      0.0      order = 2
     9         1      17789.0  17789.0      0.4      poly = cp.orth_ttr(order, dist, normed=True)
    10         1          3.0      3.0      0.0      Ns = len(poly)*2
    11         1         20.0     20.0      0.0      print(Ns)
    12         1       1953.0   1953.0      0.0      samples = dist.sample(Ns)
    13         1         15.0     15.0      0.0      y_sample =  np.sum(samples, axis=0)
    14                                           
    15         1      19713.0  19713.0      0.4      approx, uhat = cp.fit_regression(poly, samples, y_sample, retall=1)
    16         1      23543.0  23543.0      0.5      exp_pc = cp.E(approx, dist)
    17         1        250.0    250.0      0.0      print(exp_pc, uhat[0])
    18         1     824198.0 824198.0     17.3      std_pc = cp.Std(approx, dist)
    19         1         15.0     15.0      0.0      variance = np.sum(uhat[1:]**2)
    20         1         10.0     10.0      0.0      std_coefs = np.sqrt(np.sum(uhat[1:]**2))
    21         1         20.0     20.0      0.0      print(std_pc, std_coefs)
    22         1    1002124.0 1002124.0     21.0      S_pc = cp.Sens_m(approx, dist) #Si from chaospy
    23         1    2878358.0 2878358.0     60.3      S_t_pc = cp.Sens_t(approx, dist) #Si from chaospy
    24         1       5134.0   5134.0      0.1      main, total = calc_sens(dim, poly, uhat)
    25         1        503.0    503.0      0.0      print(S_pc, main)
    26         1        479.0    479.0      0.0      print(S_t_pc, total)

[jonathf]

Cool.

The major components of the refactoring is already in place. It consist of redefining the polynomial class as a numpy data type. I have split it out into a new repo which can be found here: https://github.com/jonathf/numpoly

I have now also made an experimental branch in chaospy that switches out the backend: https://github.com/jonathf/chaospy/pull/171

I will say that it is beta-ready, but a lot of stuff needs polish. I don't imagine it is any faster at this point. Bot the docs and the testing needs an overhaul.

What it really needed at this stage is testing and profiling. What works, what doesn't. Where are the bottlenecks. So if you have the time, do your analysis there. And if you patch works in the new environment, I definitely would like to see a patch.

Things to note at this stage:

• lots of functionality is removed/replaced/updated. Mostly backend stuff.
• old polynomial names were reference-by-index only, fixed to q0,...,qn. In the new system labels are used. they are stored in poly.names, which gives a reference-by-index behavior. How that plays out in practice is still a bit undecided.
• Sensetivity analysis doesn't seem to work for your example. Not sure why, I'll have to take a close look.
• A major bottleneck is in chaospy.distributions.evaluations.moments. I plan to add caching there soon, hopefully improving the performance a bit there.

[utsekaj42]

Excellent. I have updated calc_sens to work with numpoly, though I do have some questions.

What do you mean by "reference-by-index"? Now I rely on the fact that the suffix "q0", "q2" corresponds to the ordering of the joint distribution, but perhaps this assumption will be invalidated. (Otherwise, perhaps sensitivities should be returned as a dict like struct indexed on variable instead of an array.)

I haven't looked into why cp.sens_m isn't working, but as a general comment, it seems that the general case of computing sensitivities of an arbitrary polynomial wrt. an arbitrary distribution will always be less efficient than computing sensitivities with respect to a polynomial already expressed in the orthogonal basis of the distribution. Maybe a special case will exist where the approx is represented in the basis of a specific distribution and thus cp.sens_m(approx) would exploit that unless explicitly told to use a different distribution or approx is not expressed in the orthogonal basis.

[jonathf]

You are understanding the problem correctly. numpoly allows the indeterminats to have any name, not only q0, q1, .... The system is still positional in the sense that if you do xyz = numpoly.symbols("x y z"), x, y, z is index 0, 1 and 2. The problem lies in how to enforce that distribution index matches with polynomial index. Either chaospy has to enforce the qN convension, or as you suggest return values have to be change some how. Or another option, which I am playing with, is to some how tag the positional distributions against the polynomial name.

[tostenzel]

I would also love to see an option to calculate the sobol indices by post-processing the PCE coefficients instead of Monte-Carlo!

[jonathf]

So update time, sine a couple of months have passed.

Looking through the different option, I have come to the conclusion that though numpoly can do things by keyword, chaospy will not. To get around this, chaospy will import all it needs from numpoly except the constructors (like symbols, monomial, ...). These will have custom wrappers in chaospy.

This will ensure that the backward compatibility is maintained, while still moving the majority of the code canges into the dependency.

[jonathf]

@tostenzel, the current implementation have two sets of definitions: Sens_* and Sens_*_sample. The former bypasses the sampling and calculates sensitivity without Monte-Carlo.

Though implementation is faster than Monte-Carlo, it is still slower than it needs to be. This is what needs improving.

[tostenzel]

@jonathf Ah OK, got you. And great thanks for the update!

[jonathf]

@tostenzel, @utsekaj42, after the decision was made regarding approach, then adding change went fast. There is now a beta version out: v3.1.0b0.

@utsekaj42, I haven't run you example in full, but for the cp.Std part of the example, the execution time is down from over 400 sec down to under 35 sec.

Please test it out and let me know what works and what doesn't.

[jonathf]

Version 3.1.0 is now out.

I've also added a few tasks to https://github.com/jonathf/chaospy/blob/master/tasks.rst that could also improve speed with quite a bit. On my todo-list.

[gaoutham92]

I have implemented a workaround for my use case, but perhaps this can be integrated in a more general manner. Currently, this only works if one uses normed polynomials for the expansion, and then use retall in the fitting method to return the coefficients of the PCE as uhat:

import chaospy as cp
import numpy as np

@profile
def test():
    dim = 5
    dist = cp.Iid(cp.Uniform(), dim)
    order = 2
    poly = cp.orth_ttr(order, dist, normed=True)
    Ns = len(poly)*2
    print(Ns)
    samples = dist.sample(Ns)
    y_sample =  np.sum(samples, axis=0)

    approx, uhat = cp.fit_regression(poly, samples, y_sample, retall=1)
    exp_pc = cp.E(approx, dist)
    print(exp_pc, uhat[0])
    std_pc = cp.Std(approx, dist)
    variance = np.sum(uhat[1:]**2)
    std_coefs = np.sqrt(np.sum(uhat[1:]**2))
    print(std_pc, std_coefs)
    S_pc = cp.Sens_m(approx, dist) #Si from chaospy
    S_t_pc = cp.Sens_t(approx, dist) #Si from chaospy
    main, total = calc_sens(dim, poly, uhat)
    print(S_pc, main)
    print(S_t_pc, total)

def calc_sens(dim, poly, uhat):
    variance = np.sum(uhat[1:]**2)
    main = np.zeros(dim)
    total = np.zeros(dim)
    #Compute sens
    for idx, pol in enumerate(poly[1:]):
        add_total = [False,]*dim
        add_main = [True,]*dim
        for term in pol.A:
            for var in range(dim):
                if term[var] > 0:
                    add_total[var] = True
                    add_main[var] = add_main[var] and True
                    add_main[0:var] = [False,]*var
                    try:
                        add_main[var+1::] = [False,]*(dim-var)
                    except IndexError:
                        pass
        for var in range(dim):
            if add_main[var]:
                main[var] += uhat[idx+1]**2
            if add_total[var]:
                total[var] += uhat[idx+1]**2

    main = main/variance
    total = total/variance
    return main, total

if __name__ == "__main__":
    test()

which is much faster:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     4                                           @profile
     5                                           def test():
     6         1          2.0      2.0      0.0      dim = 5
     7         1         33.0     33.0      0.0      dist = cp.Iid(cp.Uniform(), dim)
     8         1          0.0      0.0      0.0      order = 2
     9         1      17789.0  17789.0      0.4      poly = cp.orth_ttr(order, dist, normed=True)
    10         1          3.0      3.0      0.0      Ns = len(poly)*2
    11         1         20.0     20.0      0.0      print(Ns)
    12         1       1953.0   1953.0      0.0      samples = dist.sample(Ns)
    13         1         15.0     15.0      0.0      y_sample =  np.sum(samples, axis=0)
    14                                           
    15         1      19713.0  19713.0      0.4      approx, uhat = cp.fit_regression(poly, samples, y_sample, retall=1)
    16         1      23543.0  23543.0      0.5      exp_pc = cp.E(approx, dist)
    17         1        250.0    250.0      0.0      print(exp_pc, uhat[0])
    18         1     824198.0 824198.0     17.3      std_pc = cp.Std(approx, dist)
    19         1         15.0     15.0      0.0      variance = np.sum(uhat[1:]**2)
    20         1         10.0     10.0      0.0      std_coefs = np.sqrt(np.sum(uhat[1:]**2))
    21         1         20.0     20.0      0.0      print(std_pc, std_coefs)
    22         1    1002124.0 1002124.0     21.0      S_pc = cp.Sens_m(approx, dist) #Si from chaospy
    23         1    2878358.0 2878358.0     60.3      S_t_pc = cp.Sens_t(approx, dist) #Si from chaospy
    24         1       5134.0   5134.0      0.1      main, total = calc_sens(dim, poly, uhat)
    25         1        503.0    503.0      0.0      print(S_pc, main)
    26         1        479.0    479.0      0.0      print(S_t_pc, total)

Hi,

The sensitivity indices calculation from Chaospy requires too much time, when I was checking the cause for that I found that the code is designed to take the conditional expectation of the response surrogate with respect to the input variables and taking the variance of the conditional expectation(that is why it is taking a lot of time). This might be the usual procedure for Monte Carlo-based Sobol' indices calculation, however, for PCE, we can just simply do basic post-processing on coefficients to get variance corresponding to each input parameter and divide it with the total variance.

Actually, this issue has been raised earlier and @utsekaj42 has given the code for the same. When I try to use the code it says numpoly has no attribute A.

I also attached a reference paper from ETH Zurich, they are actually computing the sensitivity indices by simply operating the coefficients. (Pg 161-162 sec. 3.3 in the reference)

SPC_Hydrogen_ReliabilityEngineering&_System_safety.pdf

I also tried to check the results from Chaospy and UQLab (Matlab) for 10 D problem, the condition expectation code for the sensitivity index you made give similar results to UQLab up to 10 variables (but took a much longer time). However, Chaospy sensitivity indices gave wrong values for more than 10 variables.

So can you please tell whether your approach of taking conditional expectation on polynomial surrogates is the right approach for calculating sensitivity indices from PC?. I think when we are dealing with pc coefficients, this is not needed.

                                                      Thankyou

Regards Gowtham Radhakrishnan

[jonathf]

Hello Gowtham!

Amazing that you have taken the time to look into this. I will give it attention, but likely not this week as I have other things on my plate.

Last month I added some major improvements to the way expected values are handled iternally. To utilize these improvements, you need python >= 3.7, chaospy >= 4.3.2 and numpoly >= 1.2.3. I believe it solves a lot of the issues address by both you and @utsekaj42, but some diagnostics is required to know for sure.

Can you report which version you are using, and if you are using the older version, upgrade and run again?

Oh, and pol.A was deprecated late 2019. The correct syntax is now: pol.exponents, if you want to recreate the original code example.

Jonathan

[gaoutham92]

Hi,

Thanks for the reply. I think I am using Chaospy version 3.3.8, I will upgrade and run it

Regards Gowtham R