JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

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

Why do we still use Lanczos for c ≠ integer? #114

Open dlfivefifty opened 1 month ago

dlfivefifty commented 1 month ago

Cholesky should always be faster. otherwise we've done something seriously wrong!

@tsgut @ioannisPApapadopoulos @DanielVandH

ioannisPApapadopoulos commented 1 month ago

I do not think we do. Looking at https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/blob/a1e4c3e65f164bc55f97f3574eb0da3ca62468ff/src/SemiclassicalOrthogonalPolynomials.jl#L133 it seems it only uses Lanczos if c is not an integer?

dlfivefifty commented 1 month ago

That's what I meant

TSGut commented 1 month ago

It's just a matter of a more subtle implementation. I guess we'd want to go to the integer floor of the given parameters via a hierarchy, then go the final non integer amount as a last step. It should be fine to do that, though it won't be a polynomial modification at the end, so it won't be a banded connection. Not that big of an issue though.

The answer to your question is just that it's a bit more subtle and we didn't need it. But I can get that working pretty quickly if we do.

dlfivefifty commented 1 month ago

you just approximate (t-x)^c by a polynomial

TSGut commented 1 month ago

Well, as long as what you mean is to reach $P^{a,b,c}(x)$ with $c$ a general real, we ladder to $P^{a,b,\lfloor c \rfloor}(x)$ as now and then from there use decomposition methods on the Clenshaw-computed remaining weight $(tI-X)^{c-\lfloor c \rfloor}$. Otherwise it will be unstable for high $c$.

But both of those parts are implemented. The first is just what SemiclassicalOPs.jl does and the second already works in ClassicalOPs.jl. So it's just a well-placed if condition in SemiclassicalOPs.jl to do this.

dlfivefifty commented 1 month ago

what makes you think lanczos is stable for high $c$?

TSGut commented 1 month ago

Fair point but we moved away from Lanczos for specifically that reason. You are likely right that for low c where Lanczos is reasonable you can just call the ClassicalOP.jl convertedop directly. I really would prefer we do it properly though. 😂