Add normalized polynomials

[dlfivefifty]

I was going to propose trying to stick to linear algebra language as much as possible, in which case normalisation is really just a QR. That is:

P = Legendre()
Q,D = qr(P)
qr(P) isa OrthogonalPolynomialQR{Float64}
@test Q isa Normalized{Float64,Legendre{Float64}}
@test D isa Diagonal
x,n = 0.1,5
@test P[x,n] == (Q*D)[x,n]

For other OPs we would need to weight them first, to avoid dealing with inner products:

T = Chebyshev()
wT = SqrtWeighted(T)
@test wT isa SqrtWeighted{Float64,Chebyshev{Float64}}
wQ,D = qr(wT)
@test wQ isa SqrtWeighted{Float64,Normalized{Float64,Chebyshev{Float64}}}
@test wQ[x,n] == (T*D)[x,n]/sqrt(1-x^2)

But this seems to suggest that Normalized{T,<:OrthogonalPolynomial} <: OrthogonalPolynomial is the more fundamental type.

[guilgautier]

Hi @dlfivefifty,

First, thanks for your nice work on this project!

I've recently started playing with Julia; I like the philosophy and flexibility of the language.
I am now looking to evaluate classical orthonormal polynomials (Jacobi, Laguerre, Hermite, Fourier).

As I understand, the package is based on QuasiArrays.jl and ContinuumArrays.jl, which enables fast evaluation of the orthogonal polynomials because of the special indexing features of such arrays.

Is there is a minimal version of the package (with the least dependencies possible) to simply evaluate to evaluate orthogonal, monic orthogonal and orthonormal polynomials ?

Btw, when I install the package, I can't use Normalized nor qr.

Cheers

[dlfivefifty]

I am now looking to evaluate classical orthonormal polynomials (Jacobi, Laguerre, Hermite, Fourier).

Already works in v0.4 (which I'm about to tag):

julia> P = Normalized(Legendre())
Normalized(Legendre{Float64})

julia> @time P[range(-1,1; length=1000), 1:1000] # 1000x1000 Vandermonde matrix
  0.119438 seconds (4.00 M allocations: 373.499 MiB, 20.28% gc time)
1000×1000 Array{Float64,2}:
 0.707107  -1.22474  1.58114  -1.87083  2.12132  -2.34521  2.54951  -2.73861  2.91548  -3.08221  …  31.504      -31.5198    31.5357    -31.5515     31.5674    -31.5832    31.5991   -31.6149
 0.707107  -1.22229  1.57165  -1.84841  2.07904  -2.27527  2.44339  -2.58716  2.70901  -2.81059      2.19638     -2.33678    2.46783    -2.58899     2.69979    -2.79978    2.88855   -2.96576
 0.707107  -1.21984  1.56218  -1.82611  2.03715  -2.20632  2.33939  -2.43979  2.50974  -2.55086      2.65661     -2.62344    2.56926    -2.49451     2.39978    -2.28583    2.15358   -2.00408
 0.707107  -1.21739  1.55274  -1.80392  1.99563  -2.13833  2.23749  -2.29643  2.31751  -2.30269      0.803525    -0.549829   0.289527   -0.0257482  -0.23834     0.499566  -0.75479    1.00095
 0.707107  -1.21494  1.54331  -1.78184  1.95449  -2.07131  2.13766  -2.15701  2.13217  -2.06576      1.62597     -1.80848    1.96202    -2.08414     2.17288    -2.22682    2.2451    -2.22741
 0.707107  -1.21249  1.53389  -1.75987  1.91372  -2.00525  2.03987  -2.02147  1.95356  -1.83973  …   0.0824217    0.21793   -0.513919    0.799619   -1.06931     1.31759   -1.5395     1.73058
 0.707107  -1.21003  1.5245   -1.73801  1.87333  -1.94013  1.9441   -1.88975  1.78153  -1.6243      -1.50506      1.69742   -1.849       1.95616    -2.01633     2.02805   -1.99105    1.90622
 0.707107  -1.20758  1.51513  -1.71626  1.83331  -1.87595  1.85032  -1.76177  1.61594  -1.41917     -1.15878      1.40491   -1.61165     1.77323    -1.8851      1.94414   -1.94869    1.89862
 0.707107  -1.20513  1.50578  -1.69463  1.79366  -1.81271  1.75852  -1.63748  1.45663  -1.22402      0.762819    -0.442394   0.107797    0.230252   -0.560926    0.873632  -1.15835    1.40597
 ⋮                                                ⋮                                              ⋱                                       ⋮                                           
 0.707107   1.20513  1.50578   1.69463  1.79366   1.81271  1.75852   1.63748  1.45663   1.22402      0.762819     0.442394   0.107797   -0.230252   -0.560926   -0.873632  -1.15835   -1.40597
 0.707107   1.20758  1.51513   1.71626  1.83331   1.87595  1.85032   1.76177  1.61594   1.41917     -1.15878     -1.40491   -1.61165    -1.77323    -1.8851     -1.94414   -1.94869   -1.89862
 0.707107   1.21003  1.5245    1.73801  1.87333   1.94013  1.9441    1.88975  1.78153   1.6243      -1.50506     -1.69742   -1.849      -1.95616    -2.01633    -2.02805   -1.99105   -1.90622
 0.707107   1.21249  1.53389   1.75987  1.91372   2.00525  2.03987   2.02147  1.95356   1.83973      0.0824217   -0.21793   -0.513919   -0.799619   -1.06931    -1.31759   -1.5395    -1.73058
 0.707107   1.21494  1.54331   1.78184  1.95449   2.07131  2.13766   2.15701  2.13217   2.06576  …   1.62597      1.80848    1.96202     2.08414     2.17288     2.22682    2.2451     2.22741
 0.707107   1.21739  1.55274   1.80392  1.99563   2.13833  2.23749   2.29643  2.31751   2.30269      0.803525     0.549829   0.289527    0.0257482  -0.23834    -0.499566  -0.75479   -1.00095
 0.707107   1.21984  1.56218   1.82611  2.03715   2.20632  2.33939   2.43979  2.50974   2.55086      2.65661      2.62344    2.56926     2.49451     2.39978     2.28583    2.15358    2.00408
 0.707107   1.22229  1.57165   1.84841  2.07904   2.27527  2.44339   2.58716  2.70901   2.81059      2.19638      2.33678    2.46783     2.58899     2.69979     2.79978    2.88855    2.96576
 0.707107   1.22474  1.58114   1.87083  2.12132   2.34521  2.54951   2.73861  2.91548   3.08221     31.504       31.5198    31.5357     31.5515     31.5674     31.5832    31.5991    31.6149

Is there is a minimal version of the package (with the least dependencies possible) to simply evaluate to evaluate orthogonal, monic orthogonal and orthonormal polynomials ?

No. Note this package uses LazyArrays.jl and InfiniteArrays.jl to efficiently deal with recurrence relationships; so its not even the case that you could pull out the current code to make it "lightweight": you'd need to re-implement forwardrecurrence for each class of polynomials.

[guilgautier]

I am now looking to evaluate classical orthonormal polynomials (Jacobi, Laguerre, Hermite, Fourier).

Already works in v0.4 (which I'm about to tag):

julia> P = Normalized(Legendre())
Normalized(Legendre{Float64})

Ok great, I'll check out v0.4!

Is there is a minimal version of the package (with the least dependencies possible) to simply evaluate to evaluate orthogonal, monic orthogonal and orthonormal polynomials ?

No. Note this package uses LazyArrays.jl and InfiniteArrays.jl to efficiently deal with recurrence relationships; so its not even the case that you could pull out the current code to make it "lightweight": you'd need to re-implement forwardrecurrence for each class of polynomials.

Using LazyArrays.jl and InfiniteArrays.jl seem to be the right thing to do.
It would be nice to encapsulate key features like evaluation of orthogonal polynomials in a lightweight environnement.
This would definitely benefit other projects.

[dlfivefifty]

QuasiArrays.jl and ContinuumArrays.jl are pretty lightweight so I don't think there's much benefit in dropping the dependency. Dropping FastTransforms.jl would be bigger.

But you could do it and move out the forwardrecurrence code from FastTransforms.jl, then something like this could work:

julia> using FillArrays, LazyArrays

julia> import OrthogonalPolynomialsQuasi: initiateforwardrecurrence # TODO: move to new package

julia> chebyshevtrecurrence() = ([1; Fill(2,∞)], Zeros{Int}(∞), Ones{Int}(∞))
chebyshevtrecurrence (generic function with 1 method)

julia> function chebyshevt(x::Number, n::Integer)
           n == 0 && return one(x)
           initiateforwardrecurrence(n, chebyshevtrecurrence()..., x, one(x))[2]
       end
chebyshevt (generic function with 1 method

julia> chebyshevt(0.1, 10)
-0.5388927487999999

julia> ChebyshevT()[0.1,11]
-0.5388927487999999

Adding other OPs is straightforward: just add the recurrences. You could even make things like chebyshevt.(0.1, 1:100_000) fast by overloading broadcasting to call forwardrecurrence!.

I don't have time to do this though but give you pointers if you feel like taking it up. The relevant code is:

https://github.com/JuliaApproximation/FastTransforms.jl/blob/8a000df9e2381f414fece33924bc7ac00b059570/src/clenshaw.jl#L1

OrthogonalPolynomialsQuasi.jl/src/clenshaw.jl

Line 2 in 7ee21f8

##

OrthogonalPolynomialsQuasi.jl/src/chebyshev.jl

Line 109 in 7ee21f8

recurrencecoefficients(C::ChebyshevT) = (Vcat(1, Fill(2,∞)), Zeros{Int}(∞), Ones{Int}(∞))

OrthogonalPolynomialsQuasi.jl/src/jacobi.jl

Line 192 in 7ee21f8

function recurrencecoefficients(::Legendre{T}) where T

OrthogonalPolynomialsQuasi.jl/src/jacobi.jl

Line 207 in 7ee21f8

function recurrencecoefficients(P::Jacobi)

[dlfivefifty]

FYI there was an effort to add this to SpecialFunctions.jl, w/o any of the LazyArrays.jl machinery:

JuliaMath/SpecialFunctions.jl#175

[guilgautier]

Thanks for the pointers @dlfivefifty!
I might give it a try when I will have more time and feel more comfortable with Julia.