Number of points for sparse grid

[samira-m94]

Hi, I am trying to use the sparse grid feature and unfortunately the number of points I am getting is not in agreement with the literature like:

https://link.springer.com/content/pdf/10.1007%2Fs00158-009-0441-x.pdf

for example using Gaussian quadrature for a 5 dimensional case with accuracy of 2 the code suggests 61 points for each dimension whereas the paper mentions 66. However, for the accuracy of 1 they are actually the same.

I was wondering where the difference is coming from?

[jonathf]

I am having a little bit of an issue accessing journal articles at the moment, so I can not answer concretely, but I still think I get what is going on here from what you are telling me.

First of, please take look at #139 and what I am saying about nested samples in Smolyak sparse grids.

Next I want to note that Gaussian quadrature do not really nest, but there is a few exceptions. For e.g. Hermite-Gauss quadrature the origin point repeats every second order, as the only point that nests.

So here is my educated guess:

66 quadrature points are correct, however 5 lower level grids (one for each dimension) nests into a higher order one in the origin point. In practice when you evaluate on grid points x0 ... x4 that are all the same (here origin), you can simplify:

f(x0)*w0 + ... + f(x4)*w4 =
f(x)*w0 + ... + f(x)*w4 =
f(x)*(w0 + ... + w4)

In other words, when samples are the same, one can in practice reduce the number of evaluations, by collapsing the weights. This is supported in chaospy and done automatically.

[samira-m94]

Thanks for the response. Another question I have is that as I try to run the code below for evaluating moments using sparse grid there are some jumps in the values of the moments as I increase the order.

import chaospy as cp
import numpy as np

def CPCX7_func(x):
    b = [0.5630, 0.0513, 0.4498, 0.7311, 0.3501, 0.1530, 0.8833]
    y = (b[0]+b[1]*x[0]+b[2]*x[1]+b[3]*x[0]*x[1])/(1+b[4]*x[0]+b[5]*x[1]+b[6]*x[0]*x[1])
    return y

Func = 'CPCX7'
Dist = 'Uniform'
dist = cp.Iid(cp.Uniform(0,10),2)
fcall = 0
m = 2
while fcall <= 1000:
    print('m',m)
    file = open("SG2.txt","a")
    n = 0
    M = []
    X, W = cp.generate_quadrature(m, dist,sparse=True, rule="G")
    fcall = len(W)
    while n < 4:
        b = sum((eval(Func+'_func')(X)**(n+1))*W)
        M.append(b)
        n += 1
    miu = M[0]
    std = abs(M[1] - miu ** 2) ** 0.5
    skw = (M[2] - 3 * miu * std ** 2 - miu ** 3) / std ** 3
    kurt = (M[3] - 4 * miu * M[2] + 6 * (miu ** 2) * M[1] - 3 * miu ** 4) / std ** 4
    N = [miu, std, skw, kurt]

Mean	Std	Skew	Kurt	Fcalls	order
0.833066	0.105857	0.106937	2.056261	14	2
0.833386	0.127966	0.246703	3.202096	30	3
0.834014	0.140938	0.406833	4.393493	55	4
0.834601	0.148342	0.574505	5.574581	91	5
0.66074	0.372296	-0.80092	2.613906	139	6
1.009509	0.390086	-12.1715	-9.04258	203	7
0.835366	0.156231	0.948395	8.381025	285	8
0.835429	0.15699	1.015105	8.969121	385	9
0.835462	0.157419	1.060128	9.400309	506	10
0.835478	0.157661	1.089496	9.704907	650	11
0.835487	0.157795	1.108158	9.913529	819	12
0.835492	0.157869	1.119773	10.05273	1015	13
0.75103	0.2975	-1.20039	5.252516	1239	14
0.91996	0.229992	-17.7255	-23.7969	1495	15
0.835496	0.157945	1.133709	10.23828	1785	16
0.835497	0.157951	1.135207	10.26102	2109	17
0.768346	0.276862	-1.26016	6.187896	2469	18
0.902648	0.189133	-20.011	-37.1708	2869	19
0.798205	0.234046	-1.20288	8.477458	3306	20
0.872789	0.087506	-13.0299	-341.133	3790	21
0.745034	0.304129	-1.17119	4.995767	4321	22
0.92596	0.242816	-17.0821	-21.1906	4897	23
0.75858	0.288829	-1.22885	5.639868	5522	24
0.912414	0.212962	-18.6361	-28.1846	6198	25
0.822769	0.188266	-0.45927	10.87087	6925	26
0.848226	0.118886	7.642926	-8.48902	7709	27
0.781937	0.258684	-1.27504	7.105616	8551	28
0.889057	0.15084	-22.123	-66.1418	9451	29
0.816453	0.201322	-0.79842	10.32046	10410	30

Generally when we increase number of Function evaluations we should get closer to the real values of the moments but as we increase the order and consequently number function evaluations here we get unstable numbers

[jonathf]

The good thing about sparse grid is that it scales a lot better in higher dimensions (when nestedness works as expected, which your example does not). However, it also has a much stronger assumption that the underlying function is well approximated using polynomials. By commenting out the part of your expression starting with the division sign, you will see a much higher precision. You also see the same if you switch the flag to sparse=False. I suspect that your problem isn't a good candidate for sparse grid.

Another aspect, which I don't know matter as much, is that the formula for the higher order moments uses the expanded form, which is known to be less precise than there collapsed counterpart. E.g.

f = vars()[f"{Func}_func"]
miu = cp.sum(f(X)*W)
std = np.sqrt(np.abs(cp.sum((f(X)-miu)**2*W)))
skw = cp.sum(((f(X)-miu)/std)**3*W)
kurt = cp.sum(((f(X)-miu)/std)**4*W)

might help a little.

(The fact that np.abs to get std also hints to the numerical instability of the polynomial assumption.)

[samira-m94]

Thanks for the detailed answer. I have already tried the sparse=False and I call it full factorial numerical integration (FFNI) but I wanted to assess the sparse grid performance in compare to FFNI that's why I have insist on doing this using sparse grid. But still I did not get the answer for my question which is the fact that the moments are unstable although we are increasing number of function calls. I know it could be natural to have all these instabilities according to all of the factors you mentioned above but as function evaluations in different locations are added I would expect more of a stable performance at least in lower order moments. I have tried bunch of different functions with different formulations and dimensions but it is the case for most of them

[jonathf]

As sparse grid is sensitive to noise, remove potential noise sources becomes important to maintain numerical stability. And finding out what that often require a bit of diagnostics.

I have taken the liberty to explore your problem here. The solution lies in the fact that Uniform(0, 10) is a lot less well-behaved than the canonical Uniform(-1, 1). Simply explained, moments involving subtracting almost identical large numbers really get high rounding errors for floating point numbers. Combined with sparse grid, and it all falls apart. The trick to getting around it is to model your problem in Uniform(-1, 1) space and convert that to Uniform(0, 10) yourself.

Here is a possible solution:

import chaospy as cp
import numpy as np

def CPCX7_func(x):
    b = [0.5630, 0.0513, 0.4498, 0.7311, 0.3501, 0.1530, 0.8833]
    x = (x+1)*5  # map from [-1, 1] -> [0, 10]
    y = (b[0]+b[1]*x[0]+b[2]*x[1]+b[3]*x[0]*x[1])/(1+b[4]*x[0]+b[5]*x[1]+b[6]*x[0]*x[1])
    return y

Func = 'CPCX7'
dist = cp.J(cp.Uniform(-1, 1), cp.Uniform(-1, 1))

for m in range(2, 31):
    X, W = cp.generate_quadrature(m, dist, sparse=True, rule="G")
    fcall = len(W)
    f = vars()[f"{Func}_func"]
    miu = cp.sum(f(X)*W)
    std = np.sqrt(np.abs(cp.sum((f(X)-miu)**2*W)))
    skw = cp.sum(((f(X)-miu)/std)**3*W)
    kurt = cp.sum(((f(X)-miu)/std)**4*W)
    N = [miu, std, skw, kurt]
    print(f"{miu:.4f} {std:.4f} {skw:.4f} {kurt:.4f} {fcall:4d} {m:2d}")

Here is the output when I run this code:

0.4624 0.2954 1.5800 2.6392   13  2
1.2040 0.4633 -0.9166 0.8857   29  3
0.8340 0.1409 0.4068 4.3935   53  4
0.8346 0.1483 0.5745 5.5746   89  5
0.8350 0.1525 0.7282 6.6688  137  6
0.8352 0.1549 0.8539 7.6145  201  7
0.8354 0.1562 0.9484 8.3810  281  8
0.8354 0.1570 1.0151 8.9691  381  9
0.8355 0.1574 1.0601 9.4003  501 10
0.8355 0.1577 1.0895 9.7049  645 11
0.8355 0.1578 1.1082 9.9135  813 12
0.8355 0.1579 1.1198 10.0527 1009 13
0.8355 0.1579 1.1269 10.1435 1233 14
0.8355 0.1579 1.1312 10.2017 1489 15
0.8355 0.1579 1.1337 10.2383 1777 16
0.8355 0.1580 1.1352 10.2610 2101 17
0.8355 0.1580 1.1361 10.2750 2461 18
0.8355 0.1580 1.1366 10.2834 2861 19
0.8355 0.1580 1.1369 10.2885 3301 20
0.8355 0.1580 1.1370 10.2915 3785 21
0.8355 0.1580 1.1371 10.2933 4313 22
0.8355 0.1580 1.1372 10.2944 4889 23
0.8355 0.1580 1.1372 10.2950 5513 24
0.8355 0.1580 1.1372 10.2954 6189 25
0.8355 0.1580 1.1372 10.2956 6917 26
0.8355 0.1580 1.1372 10.2957 7701 27
0.8355 0.1580 1.1372 10.2958 8541 28
0.8355 0.1580 1.1373 10.2958 9441 29
0.8355 0.1580 1.1373 10.2958 10401 30

Let me know if this answers your question.

Side note: I am wondering if there is a way to achieve this automatically inside chaospy. I will have to explore what is possible.