TSGut
commented
1 year ago Now also featuring the weighted conversion speedup promised to @ioannisPApapadopoulos. Hopefully this getting faster and faster eventually reaches a point where it is fast enough (I suspect now it would just be code optimization rather than better formulas tho).
Here is the example from https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/issues/78 reran in current implementation:
julia> using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials
julia> t = 1.1
1.1
julia> Q₀₀ = SemiclassicalJacobi.(t, 0, 0, 0:∞);
julia> Q₁₁ = SemiclassicalJacobi.(t, 1, 1, 0:∞);
julia> L = (Weighted.(Q₀₀) .\ Weighted.(Q₁₁));
julia> @time L[20]; # before: 0.183541 seconds
0.060309 seconds (191.13 k allocations: 8.075 MiB)
julia> @time L[30]; # before: 0.396986 seconds
0.100929 seconds (402.88 k allocations: 15.625 MiB)
julia> @time L[40]; # before: 0.861748 seconds
0.212241 seconds (758.47 k allocations: 27.565 MiB)
julia> @time L[50]; # before: 1.565395 seconds
0.291458 seconds (1.26 M allocations: 43.963 MiB)
julia> @time L[60]; # before: 6.235148 seconds (with compile), 2.743544 seconds (without compile)
0.441863 seconds (1.95 M allocations: 64.999 MiB)and we also get
julia> @time L[250]
37.124558 seconds (142.60 M allocations: 3.696 GiB, 2.71% gc time)
(((ℵ₀×ℵ₀ LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Adjoint{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{Bool, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ LinearAlgebra.Diagonal{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * ((ℵ₀×ℵ₀ LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Adjoint{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ LinearAlgebra.Diagonal{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ LinearAlgebra.Diagonal{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
0.00434421 ⋅ ⋅ ⋅ …
0.00428913 0.00608072 ⋅ ⋅
-3.76572e-5 0.0059776 0.00737181 ⋅
⋅ -7.88228e-5 0.00721544 0.00842674
⋅ ⋅ -0.000127024 0.00821237
⋅ ⋅ ⋅ -0.000180951 …
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋮ ⋱Granted that is not exactly fast but there are a LOT of cholesky decompositions that need to be done to construct this. I am sure there is a lot of code level optimizations left here but let's see if this is good enough first.
codecov[bot]
ioannisPApapadopoulos
dlfivefifty
@ioannisPApapadopoulos @dlfivefifty
Basically just porting over @MikaelSlevinsky's contiguous recurrence implementation in a way that should be reasonably natural for this package. This should be both much faster and hopefully much more accurate at high c.
Previously we had
Now we get
At no cost to accuracy (the old method using bigfloat, the new one just Float64):
It's implemented to work whenever sum is broadcast over a crange vector of weights as seen above, as it is it wont be used if you just call sum on a single weight (mainly because doing that in a clean way would involve a lot of clauses about a, b and c and other contiguous recurrences), maybe we implement that at some point in the future. If you must use it on a single weight just do e.g.
sum.(SemiclassicalJacobiWeight.(1.1,2,3,3:3)).Still need to add and run tests.