Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Replace Lanzos with Cholesky #70

Merged
merged 7 commits into from
Apr 27, 2023
Merged

WIP: Replace Lanzos with Cholesky #70

merged 7 commits into from
Apr 27, 2023

Conversation

TSGut
Copy link
Member

@TSGut TSGut commented Jan 23, 2023

@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:
JuliaLinearAlgebra/InfiniteLinearAlgebra.jl#124
JuliaLinearAlgebra/InfiniteLinearAlgebra.jl#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
Copy link
Member Author

TSGut commented Apr 19, 2023

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 JuliaApproximation/SingularIntegrals.jl#9

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

@TSGut TSGut merged commit 0bb78e2 into JuliaApproximation:master Apr 27, 2023
@TSGut TSGut mentioned this pull request Apr 27, 2023
@dlfivefifty
Copy link
Member

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants