JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

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

WIP: Replace Lanzos with Cholesky #70

Closed TSGut closed 1 year ago

TSGut commented 1 year ago

@ioannisPApapadopoulos @dlfivefifty

I have started the process replacing Lancos with Cholesky. There are some workarounds present in the current version which can be removed once these two issues are resolved: https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/issues/124 https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/issues/123

The package currently requires ClassicalOrthogonalPolynomials.jl from my branch here. That branch has an adhoc workaround for the isposdef try/catch bug mentioned above so it is prone to break until we fix that.

Furthermore, I currently have to make use of and include SelectInfiniteBand in this repo, which I described here. I expect that this functionality will end up being moved into e.g. InfiniteLinearAlgebra.jl or something to that effect. For now I keep it here so I can test the functionality. No longer needed due to simpler but still temporary workaround mentioned here.

Please bear with me as what I am trying to do here is basically ripping out the spine of this package and replacing it with a new one -- there will be a lot of bugs to fix; I am happy to get bug reports!

All these caveats aside, here is a demo of something that already works which might be relevant for @ioannisPApapadopoulos :

julia> ρ=2/100; t=inv(1-ρ^2);

julia> @time T = SemiclassicalJacobi(t,1,0,100)
  0.446175 seconds (12.08 k allocations: 10.336 MiB)
SemiclassicalJacobi with weight x^1.0 * (1-x)^0.0 * (1.0004001600640255-x)^100.0

julia> @time T2 = SemiclassicalJacobi(t,0,0,99)
  0.333058 seconds (11.37 k allocations: 9.088 MiB)
SemiclassicalJacobi with weight x^0.0 * (1-x)^0.0 * (1.0004001600640255-x)^99.0

julia> @time T \ T2
  0.000432 seconds (244 allocations: 18.938 KiB)
(ℵ₀×ℵ₀ Diagonal{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiBands{Float64}, Tuple{Int64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiBands{Float64}, Tuple{Int64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}}, LazyBandedMatrices.SymTridiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiBands{Float64}, Tuple{Int64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}}}, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiBands{Float64}, Tuple{Int64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}}}}}}, Diagonal{Float64, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 1.0  0.970729  -0.0192308    ⋅           ⋅           ⋅          ⋅        ⋅    ⋅    ⋅   …  
  ⋅   1.3802     1.32704    -0.0394673    ⋅           ⋅          ⋅        ⋅    ⋅    ⋅      
  ⋅    ⋅         1.65075     1.5716     -0.0757516    ⋅          ⋅        ⋅    ⋅    ⋅      
  ⋅    ⋅          ⋅          1.86393     1.49156    -0.93884     ⋅        ⋅    ⋅    ⋅      
  ⋅    ⋅          ⋅           ⋅          2.2975      0.742187  -1.98953   ⋅    ⋅    ⋅      
 ⋮                                                   ⋮                                  ⋱  

(PS: I haven't yet done any sanity testing on whether the above is a reasonably accurate, too much still to do before I get to that. In low weight params where I can test easily it all seems correct.)

TSGut commented 1 year ago

The last push (temporarily) adds SingularIntegrals.jl via main branch and the same for ClassicalOPs (since the QR/Cholesky jacobi matrices are not tagged yet) and updates the requirements on several package dependencies. See https://github.com/JuliaApproximation/SingularIntegrals.jl/issues/9

Before merging (when it stops being WIP), we should replace this with registered versions.

dlfivefifty commented 1 year ago

Please remove the "WIP" in the PR name before merging