Legendre, Hermite, Laguerre, Chebyshev polynomials

[PaulXiCao]

Legendre polynomial is implemented via the Bonnet recursion formula.

Documentation + testing is already inlcuded.

This is in relation to issue #124 .

[ludoro]

Nice!
I will soon be in need of Legendre, Hermite, Laguerre and Jacobi polynomials for my package.
I guess that a recursion algorithm might be slow but I do not know how other packages compute those! I have been looking around for some papers for the computational aspects but haven't found any, do you have by any chance some recursive formula for the other? I'd love to open a PR to help out!

[PaulXiCao]

Nice!
I will soon be in need of Legendre, Hermite, Laguerre and Jacobi polynomials for my package.
I guess that a recursion algorithm might be slow but I do not know how other packages compute those! I have been looking around for some papers for the computational aspects but haven't found any, do you have by any chance some recursive formula for the other? I'd love to open a PR to help out!

I havent looked into those yet but I'd recommend to look at the implementation of python's numpy. I found their Hermite evaluation at github.com/numpy/....

[PaulXiCao]

@simonbyrne I dont get why the checks fail. Is the build process somehow broken?

[PaulXiCao]

Nice!
I will soon be in need of Legendre, Hermite, Laguerre and Jacobi polynomials for my package.
I guess that a recursion algorithm might be slow but I do not know how other packages compute those! I have been looking around for some papers for the computational aspects but haven't found any, do you have by any chance some recursive formula for the other? I'd love to open a PR to help out!

@ludoro You can easily implement Hermite, Laguerre, Legendre, Chebyshev polynomials via recurrence relations. They can all be described by the same recurrence relation and I implemented Hermite and Legendre already. Check out my last commits.

[PaulXiCao]

@simonbyrne I dont plan on implementing any new features for this branch. Let me know what you think of it.

I need this branch to start working on implementing spherical harmonics..

[simonbyrne]

You need to add Polynomials to the test dependencies in the Project.toml (similar to the Test package: https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/Project.toml)

[PaulXiCao]

You need to add Polynomials to the test dependencies in the Project.toml (similar to the Test package: https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/Project.toml)

Done.

[PaulXiCao]

continuous-integration/drone/pr still fails because it only runs julia 1.0.
thus, range doesnt accept stop parameter as positional one..

[PaulXiCao]

Quote from @MikaelSlevinsky

For evaluating degree-n Chebyshev polynomials, the recurrence relations cost O(n) per point. Why not use

* the explicit formulas, say, `T_n(x) = cos(n*acos(x))` (I think the answer is low accuracy); or,

* the O(log n) algorithm based on reducing the degree such as `T_{2n}(x) = 2 T_n(x)^2 - 1` in the even case or `T_{2n+1}(x) = 2 T_{n+1}(x) T_n(x) - x` in the odd case? 

I liked the formula to reduce the degree and implemented that in the most current commit. Do you know similar formulas for the other polynomials?

As to the explicit formula: Usually recurrence formulas should be more efficient (in terms of speed as well as accuracy) but maybe there are parameter ranges where it is better to use the explicit formula? Maybe we could benchmark that at some point.. that could also be done in another issue...

[MikaelSlevinsky]

Thanks for taking this on, @PaulXiCao; it's a big task.

Most of the functions in this library are specialized to a set of floating-point types, with some fallbacks for other numbers. In this sense, what makes including orthogonal polynomials challenging is not so much evaluating the recurrences, but that doing so defines the scope of orthogonal polynomials for the entire community.

According to Wikipedia, Chebyshev polynomials are also defined for non-natural degree by the formula T_n(x) = cos(acos(n*x)), provided -1 <= x <= 1, and other similar-looking formulas exist for other values of x. There is simply no comparison between this formula with the degree doubling formulas or the three-term recurrences for evaluating chebyshevt(0.123, 0.456).

On the other hand, the O(log n) complexity for the degree doubling formulas is significantly faster than the three-term recurrences when x is a square matrix, but it's not clear when comparing with the analytical formula, which has essentially fixed cost independent of degree. For low degree, the degree doubling seems better, but that changes when the degree gets large:

julia> chebyshevt(n, x) = cos(n*acos(x))
chebyshevt (generic function with 1 method)

julia> function chebyshevt(n::Integer, x)
           if n == 0
               return one(x)
           elseif n == 1
               return x
           elseif iseven(n)
               Tn2 = chebyshevt(n÷2, x)
               return 2*Tn2*Tn2-one(x)
           else
               Tn2 = chebyshevt(n÷2, x)
               Tn2p1 = chebyshevt(n÷2+1, x)
               return 2*Tn2*Tn2p1-x
           end
       end
chebyshevt (generic function with 2 methods)

julia> x = randn(1000, 1000)/sqrt(2000); x = (x+x')/2;

julia> @time chebyshevt(100, x);
  0.321130 seconds (120 allocations: 442.509 MiB, 9.06% gc time)

julia> @time chebyshevt(100.0, x);
  1.351499 seconds (84 allocations: 420.002 MiB, 15.79% gc time)

julia> @time chebyshevt(10000, x);
  2.151449 seconds (728 allocations: 2.697 GiB, 14.03% gc time)

julia> @time chebyshevt(10000.0, x);
  1.510139 seconds (96 allocations: 511.555 MiB, 2.31% gc time)

julia> 

Another tradeoff is that it appears that the degree-doubling formulas for Chebyshev polynomials have a higher error growth rate than the traditional three-term recurrences.

All four types of Chebyshev polynomials, T, U, V, W, would satisfy a degree-doubling-like formula because they're essentially trig functions with a variable transformation.

The classical orthogonal polynomials have so-called interior and exterior asymptotic approximations, some of which are mathematically rigorous; see the DLMF https://dlmf.nist.gov/18.15. Implementation of these formulas would naturally come with branching in the code but would make evaluating e.g. Legendre polynomials of extreme degree, say P_{1e24}(0.248), feasible.

[MikaelSlevinsky]

I'm also beginning to wonder whether this is the appropriate package for classical orthogonal polynomials instead of something like JuliaMath/ClassicalOrthogonalPolynomials.jl

All classical OPs are of hypergeometric type. So they could be evaluated for argument, parameters, and degree that are ::Number by calling https://github.com/JuliaMath/HypergeometricFunctions.jl (which depends on https://github.com/JuliaMath/SpecialFunctions.jl).

[giordano]

I'm also beginning to wonder whether this is the appropriate package for classical orthogonal polynomials instead of something like JuliaMath/ClassicalOrthogonalPolynomials.jl

A package like https://github.com/jverzani/SpecialPolynomials.jl? CC: @jverzani

[MikaelSlevinsky]

That repository should use FastGaussQuadrature.jl for faster algorithms than GLR and FastTransforms.jl for the connection problems.

Most of these variations on the same ideas are already unified under one tent: ApproxFun or more generally the JuliaApproximation org.

[dlfivefifty]

Also, it has its own implementation of pFq instead of calling HypergeometricFunctions.jl...

The argument for these being in SpecialFunctions.jl is it makes it a “canonical” implementation that can be refined, where putting it in yet another package will mean that they will continue to be reimplemented.

[jverzani]

Yes, agree with all of that. I put a notice in the docs for SpecialPolynomials to use FastGaussQuadratures for any real work for finding quadrature nodes and weights, but really should put in notices elsewhere to redirect if performance is important. The only advantage here, is perhaps a bit more generality in how the orthogonal polynomials are defined, but that is a stretch. The one in SpecialPolynomials is more for pedagogical purposes. Thanks for the pointer about FastTransforms, I hadn't stumbled across that package. Will have a look. Finally, the pFq function is really just used for testing, and I wanted to be able to pass a variable from Polynomials, rather than a number, and that didn't seem to work with the one from HypergeometricFunctions. The one in SpecialPolynomials doesn't do anything very reasonable when the series isn't finite..

[MikaelSlevinsky]

SpecialPolynomials.jl a pretty cool looking package!

There's a lot to do to make HypergeometricFunctions.jl more user-friendly, so if you have any use cases, please file issues.

The package is more than "just the Maclaurin series," though. As an example, if one defines Chebyshev polynomials by Gauss' hypergeometric function, then

julia> using HypergeometricFunctions

julia> chebyshevt(n, x) = _₂F₁(-n, n, 1//2, (1-x)/2)

julia> @time chebyshevt(0.123, 0.456)
  0.000003 seconds (5 allocations: 176 bytes)
0.9909056375600692

the actual branch that's called is https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/62ca7c6704cc27c96bbcc76e3573bd780b9ba1e8/src/gauss.jl#L16 which is pretty efficient (further defined here https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/62ca7c6704cc27c96bbcc76e3573bd780b9ba1e8/src/specialfunctions.jl#L81).

[jverzani]

Hi @MikaelSlevinsky that is impressive performance. SpecialPolynomials can evaluate the polynomial in that time, but if you include construction is about 5 times slower :( The only thing I was missing with these implementations was being able to use x=Polynomials.variable(); chebyt(5,x), say. For that path needs something with sqrt defined. I'll see how to add methods to call into HypergeometricFunctions for evaluation.

[PaulXiCao]

I guess someone needs to decide if this pr is going forward or if it should rather go to a different repo..

[dlfivefifty]

Let's do a poll: Thumbs up for this repos and Thumbs down for it in a different repos

[PaulXiCao]

I will merge the latest changes from master and resolve the conflicts in the hope that this pr will be merged at some point. Hope dies last ;)

[codecov]

Codecov Report

Merging #175 (67ded57) into master (281292f) will increase coverage by 0.52%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #175      +/-   ##
==========================================
+ Coverage   88.17%   88.69%   +0.52%     
==========================================
  Files          11       12       +1     
  Lines        2630     2751     +121     
==========================================
+ Hits         2319     2440     +121     
  Misses        311      311              

Flag	Coverage Δ
unittests	88.69% <100.00%> (+0.52%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/SpecialFunctions.jl	100.00% <ø> (ø)
src/14_legendre.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 281292f...67ded57. Read the comment docs.

[stevengj]

As I've commented on several of these issues, I tend to feel that orthogonal polynomials belong in a separate package.

[dlfivefifty]

Note I just registered the package:

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl

This supports some of these functions, but not all. (It was previously called OrthogonalPolynomialsQuasi.jl)

[ViralBShah]

It seems like the consensus here is that these belong in a separate package and that @dlfivefifty has registered ClassicalOrthogonalPolynomials.jl. I'm closing this issue for that reason.