JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

A Julia repository for semiclassical orthogonal polynomials
MIT License
7 stars 3 forks source link

Constructing with `t = 1` and `c` non-integer takes an insanely long time #119

Open DanielVandH opened 3 months ago

DanielVandH commented 3 months ago
julia> @time SemiclassicalJacobi(1.0, 1/2, 1.0, 1/2)
284.716841 seconds (218.66 k allocations: 4.023 GiB, 0.34% gc time)
SemiclassicalJacobi with weight x^0.5 * (1-x)^1.0 * (1.0-x)^0.5 on 0..1

I can of course just get around this by checking first if t=1 and then replacing SemiclassicalJacobi(t, a, b, c) with SemiclassicalJacobi(t, a, b + c, 0.0), but could be nice to fix this at some point

TSGut commented 3 months ago

I'm surprised this even works. Is it using Lanczos fallback here?

We could always assert t>1 somewhere.

DanielVandH commented 3 months ago

It's evaluating this

jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), jacobi(b, a, UnitInterval{T}())))

Would it be more appropriate to assert t >= 1 instead and then, in semiclassical_jacobimatrix, check isone(t)? Do these polynomials technically still make sense with t = 1? I assumed $P^{1, (a, b, c)} \equiv P^{1, (a, b + c, 0)}$ which just becomes an affine Jacobi.

TSGut commented 3 months ago

What you wrote is definitely correct mathematically.

Since for t=1 the weight can have a singularity on the domain boundary I suspect there will be convergence issues unless one is careful which may be why it's so slow as is. I think redirecting it to b+c could make sense but whether to do that or just ban t=1 is a design philosophy question

If you think redirecting would be useful you can definitely try implementing it, probably first would want to replace the Lanczos fallback though if that hasn't been done yet.

DanielVandH commented 3 months ago

Since it might be a bit questionable, for now I'll probably just work around it my own code and if it comes up again I'll think about it some more. Working with it in my own code might reveal some more reasons not to use it as well.

dlfivefifty commented 3 months ago

in an ideal world we come up with an algorithm whose cost is uniform as t -> 1 but in the meantime probably best to just throw an error