Smolyak grids diverge for higher polynoamial orders

[robertsawko]

This is going to be just a quick reality check as I may be just using the method incorrectly. I would appreciate a comment on this. I am running the following code

def cp_pseudospectral(max_order=10, sparse=False):
    '''Test of the pseudospectral method'''
    errors = zeros(max_order)

    for order in range(max_order):
        P, norms = orth_ttr(order, dist, retall=True)
        nodes, weights = generate_quadrature(
            2 * order, dist, sparse=sparse, rule='G')
        solves = [u(t, s[0], s[1]) for s in nodes.T]
        u_hat = fit_quadrature(P, nodes, weights, solves, norms=norms)
        errors[order] = L2norm_in_random_space2D(u, u_hat)
    return errors

where u is the exponent function. I am getting this behavior

Smolyak can stay on the Tensor grid depending on the quadrature order I select, but it diverges in the end. Is that an expected behavior?

[jonathf]

I have experienced it before, yes.

Smolyiak sparse grid is much more dependent upon the smoothness of the function it models. If the smoothness criteria doesn't hold, Smolyak will diverge hard.

Also worth noting a couple of things:

• Smolyak is for batteling the curse of dimensionality. I doubt you will gain anything in a two dimenstional problem.
• Smolyak works better when the samples scheme is nested. Using rule="G" implies that they are not. Unless you use infinite bounded distributions (like Gaussian, gamma, etc.), I would recommend using sparsegrid together with something like Clenshaw-Curtis (rule="C").

[jonathf]

Having said that, try your code on a simple polynomial with order at most half the number of the polynomial order. If it still diverges, it might be a bug, and I will take a look.

[robertsawko]

I am getting strange results. I will fork your repository and submit my example to your doc. I may be making some more mistakes. In the worst case you will get another example to your documentation.

[jonathf]

I looked at your example, and it seems that it behaves as expected.

Since you have only a 2 dimensional problem, fixing the polynomial order (order_p) to collocation order (order_c) will cause trouble. There simply isn't enough samples for the sparse grid to function properly.

I'd suggest that you as bare minimum at least keep order_p<=order_c+1. And more if you having convergence problems.

[robertsawko]

I agree, I kind of forgot that this is a 2D case and was extending from my 1D example. Clearly I need more quadrature point and multi-dimensional theory is not so elegant, is that right? I've tried +1 and I also tried doubling but qualitatively the results look the same. They drop after I reach 4 and then slowly increase.

The point collocation schemes look more natural. They have a larger error but at least it's decreasing with polynomial order.

With my examples I am trying to follow certain minimality so with point collocation I picked the number of samples equal to the number of polynomials. Is there any theory which allows me to choose a number of quadrature points for the multi-dimensional case?

[jonathf]

It is the other way around. For the full grid, the number of samples is about order^dim, which is significantly larger than the number of polynomials which is comb(order+dim, dim).

For Smolyak, there isn't a clear cut formula (as far as I remember), but for the most part it is lower. (More so with Clenshaw-Curtis than Gaussian quadrature). And the way it is built it is by definition more prone to non-linearities.

I did order_c=order_p+1, and I got the following results from your code:

ps_t:
[  1.02942810e+00   1.34370962e-01   2.80030550e-15   9.01096130e-15
   3.04751800e-14   1.34021917e-13   8.11765366e-13   2.34057208e-12]
ps_s:
[  1.02942810e+00   1.22958640e+00   2.28656840e+00   1.41114798e-14
   4.51674805e-14   2.05472640e-13   1.19216404e-12   4.44270003e-12]

I would like to note that 10^15 is about as accurate as you could get with my software. When it gets there it has converged completely. What we seeing here are likely rounding errors from each polynomial coefficient is going to accumulate. That is why it rises again in relative numbers. In absolute numbers however, it has converged.

[robertsawko]

This is I am getting too and I realize that we're at the level of a rounding error. It just surprised me that it is consistently increasing for both tensor and Smolyak. I think we can close this issue here, as I am satisfied that I am using the software correctly.

Thanks for all your replies.