Underflow in stieltjes() when trying to compute recursion coefficients for very high degree

[HumpyBlumpy]

Hello,

I am trying to compute a high number of recursion coefficients (on the order of 500).

Using the OrthoPoly constructor as in the tutorial returns the error

ERROR: DomainError with 2.2250738585072014e-307:
Underflow in stieltjes() for k=255; try using `removeZeroWeights`

I do not see a way to use `removeZeroWeights'.
I would be grateful for any help or tips in general about obtaining such high number of coefficients.

[timueh]

Hi there, would you be willing to post an MWE? I'll take a look then. :)

[HumpyBlumpy]

Hi,

Here is an example, almost directly pulled from the documentation:

using PolyChaos

supp = (0, 1)
w(t) = t
my_meas = Measure("my_meas", w, supp, false, Dict())
my_op = OrthoPoly("my_op", 500, my_meas; Nquad=1000);
show(my_op)

This is of course trivial (these are just Jacobi polynomials) but it already breaks down. I am interested in more complicated weight functions ( x*cos(x a)^2 or even more complicated).

[timueh]

I'll take a look. Sorry for the delay. :)

I'll try the following (but go ahead and do it yourself): just call stieltjes directly, setting the removeZeroWeights parameter to false. Also, note that the underflow limit is hard-coded (currently)

https://github.com/timueh/PolyChaos.jl/blob/88e49689ca1511b5e99ba87177674322d1c03935/src/stieltjes.jl#L21

I'm not sure how the nature of the weight function influences the underflow --> stieltjes is a numerical procedure after all that will never attain the same exactness as analytical solution. Sure -- that doesn't help you :D

[timueh]

just out of curiosity -- what's your use case anyway?

[HumpyBlumpy]

Thanks for your help!

I figured out how to call stieltjes directly. I tried setting removeZeroWeights=false but it did not change anything.
In fact, it seems that problem occurs exactly at n=255. So n<255 works, but n>= 255 does not.
Seems to be independent of the weight function (at least for the few that I checked).

The use case is basically performing a change of basis to simplify numerical simulations of quantum many-body hamiltonians (see https://arxiv.org/abs/1006.4507).

[timueh]

Mh, I guess reaching the numerical limit of the method. Have you tried the lanczos method? But I doubt it's going to be better.

Do you really need so many coefficients anyway? Perhaps there is a way you can resort to a weight function to which there is an analytical solution to the three-term recurrence coefficients? But I'm sure you have thought about that.

[HumpyBlumpy]

Hi,

Sorry for the long delay.
I tried the lanczos method and it seems to work, in that it returns me the alpha, beta coefficients which seem reasonable. I have checked them against exact solution(when it exists) and they matched.

I stumbled upon another, potentially related issue.

I need to be able to evaluate the orthonormal polynomials, so I wrote a function that uses evaluate and simply divides by the norm of the monic polynomial:

function eval_pol(x,n,alphas,betas)
    return evaluate(n,x,alphas,betas)/sqrt(computeSP2(n,betas)[end])
end

The problem is that the norm of the monic polynomials appears to decrease rapidly with n. For n around 250, the norms are about 10^-150, and then shortly after my function just returns NaNs. The range of the monic polynomials themselves is extremely small, so we just end up dividing extremely small numbers.
What I ultimately need is to evaluate integrals of the form where are the orthonormal polynomials.

The orthonormal polynomials are much better behaved than the monic ones. For instance, for n=200, the range is around [-50,50], except at k=0 where it sharply increases to 400, which are very reasonable numbers (and I probably dont even care about the "divergence" at 0 since my f(0)=0). So I think what I am trying to do should be well behaved, even for large n.

It seems the issue is dealing with monic polynomials rather than the orthonormal ones.
Is there a way of directly evaluating the orthonormal polynomials, given the recursion coefficients, without first evaluating the monic one and dividing by the norm?
Could that also be related to why the stieltjes method fails?

[timueh]

Hi there,

isn't there an elegant way to compute the recurrence coefficients of the orthonormal polynomials directly from the recurrence coefficients of the monic orthogonal polynomials? I am almost positive. It must be in Gautschi's book on orthogonal polynomials.

Here's the thing: if you had the recurrence coefficients of the orthonormal polynomials instead of the orthogonal monic ones, then you're fine.