Lanczos -> Cholesky

[TSGut]

Going to try and finish this now since a lot of the dependencies have been merged recently. Continued from #70

[TSGut]

Still need to fix a bunch of stuff but the core functionality works now and we can do

julia> P = SemiclassicalJacobi.(1.1,0:100,0:100,0:100)
101-element SemiclassicalOrthogonalPolynomials.SemiclassicalJacobiFamily{Float64, UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}:
 SemiclassicalJacobi with weight x^0.0 * (1-x)^0.0 * (1.1-x)^0.0
 SemiclassicalJacobi with weight x^1.0 * (1-x)^1.0 * (1.1-x)^1.0
 SemiclassicalJacobi with weight x^2.0 * (1-x)^2.0 * (1.1-x)^2.0
 SemiclassicalJacobi with weight x^3.0 * (1-x)^3.0 * (1.1-x)^3.0
 SemiclassicalJacobi with weight x^4.0 * (1-x)^4.0 * (1.1-x)^4.0
 SemiclassicalJacobi with weight x^5.0 * (1-x)^5.0 * (1.1-x)^5.0
 SemiclassicalJacobi with weight x^6.0 * (1-x)^6.0 * (1.1-x)^6.0
 SemiclassicalJacobi with weight x^7.0 * (1-x)^7.0 * (1.1-x)^7.0
 ⋮
 SemiclassicalJacobi with weight x^93.0 * (1-x)^93.0 * (1.1-x)^93.0
 SemiclassicalJacobi with weight x^94.0 * (1-x)^94.0 * (1.1-x)^94.0
 SemiclassicalJacobi with weight x^95.0 * (1-x)^95.0 * (1.1-x)^95.0
 SemiclassicalJacobi with weight x^96.0 * (1-x)^96.0 * (1.1-x)^96.0
 SemiclassicalJacobi with weight x^97.0 * (1-x)^97.0 * (1.1-x)^97.0
 SemiclassicalJacobi with weight x^98.0 * (1-x)^98.0 * (1.1-x)^98.0
 SemiclassicalJacobi with weight x^99.0 * (1-x)^99.0 * (1.1-x)^99.0
 SemiclassicalJacobi with weight x^100.0 * (1-x)^100.0 * (1.1-x)^100.0

julia> X = jacobimatrix.(P)
jacobimatrix.(101-element SemiclassicalOrthogonalPolynomials.SemiclassicalJacobiFamily{Float64, UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}):
 [0.5 0.28867513459481287 …  … 0.28867513459481287 0.5 …  … … ; … ; ]
 [0.41666666666666663 0.20749832663314555 …  … 0.20749832663314555 0.4727342549923196 …  … … ; … ; ]
 [0.39169675090252704 0.1681738384210662 …  … 0.1681738384210662 0.44241100225595364 …  … … ; … ; ]
 [0.38009797220323516 0.14490005502018063 …  … 0.14490005502018063 0.4230522045733082 …  … … ; … ; ]
 [0.373428048624142 0.12916692265397328 …  … 0.12916692265397328 0.4101374615066672 …  … … ; … ; ]
 [0.36909988941748256 0.117636942113484 …  … 0.117636942113484 0.40099223002606 …  … … ; … ; ]
 [0.3660653124171152 0.10872422803705963 …  … 0.10872422803705963 0.3941981988492811 …  … … ; … ; ]
 [0.3638200630909342 0.10156959856743515 …  … 0.10156959856743515 0.38895938015609216 …  … … ; … ; ]
 ⋮
 [0.35004220804369834 0.029374771188954582 …  … 0.029374771188954582 0.3524682844710059 …  … … ; … ; ]
 [0.35002920982345737 0.02921950301665129 …  … 0.02921950301665129 0.35243002816644325 …  … … ; … ; ]
 [0.35001648330579516 0.02906667121792179 …  … 0.02906667121792179 0.3523925640156556 …  … … ; … ; ]
 [0.3500040200597182 0.028916212735689447 …  … 0.028916212735689447 0.35235586766735966 …  … … ; … ; ]
 [0.3499918119994847 0.028768066774281436 …  … 0.028768066774281436 0.35231991575821453 …  … … ; … ; ]
 [0.349979851367106 0.028622174696214853 …  … 0.028622174696214853 0.35228468586323247 …  … … ; … ; ]
 [0.3499681307159118 0.028478479924682394 …  … 0.028478479924682394 0.3522501564491314 …  … … ; … ; ]
 [0.3499566428950835 0.02833692785136984 …  … 0.02833692785136984 0.3522163068304576 …  … … ; … ; ]

[codecov]

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.29 🎉

Comparison is base (f0f7afe) 93.03% compared to head (7d17bcd) 93.33%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #75      +/-   ##
==========================================
+ Coverage   93.03%   93.33%   +0.29%     
==========================================
  Files           4        4              
  Lines         704      585     -119     
==========================================
- Hits          655      546     -109     
+ Misses         49       39      -10     

Impacted Files	Coverage Δ
src/SemiclassicalOrthogonalPolynomials.jl	90.87% <100.00%> (+2.59%)	⬆️
src/neg1c.jl	100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[TSGut]

Okay, with all tests passing (except for coverage which I will keep pushing up as I go) this is a good point to pause and look at what works and what still needs to be done.

What works:

• Jacobi matrices
• hierarchical construction (as seen in previous comment)
• decomposition based raising operators via ldiv
• expansions and evaluations
• uses QR or Cholesky depending on even and odd integer shifts etc.

Basically 95% of the functionality of the old version is reproduced now, but...

What is missing:

• There is a new bug somewhere (not quite pinned down) that is causing powers of jacobi matrices to not compute under certain circumstances (i.e. P.X^3 or something like that will just get stuck in an infinite loop) and a possibly unrelated adaptive cholesky bug where finite cholesky is called instead and we get stuck in an infinite loop there as well. This is currently solved by an ad hoc workaround where I decompose the Jacobi matrix via something like Symmetric(BandedMatrix(0=>Q.X.dv, 1=>Q.X.ev). This is costly and should not be needed of course, so we need to fix those bugs before this can be replaced with just Q.X in all places. Possibly something strange is going on with SymTridiagonal infinite matrices. This is not using the new LazyArrays, so it's unrelated to that.
• Lowering via QL. This should be next since the QL PR is pretty much ready to go once the LazyArrays PR is merged

[TSGut]

@ioannisPApapadopoulos With the last change this now will always use the step-by-step hierarchy, even if you just call

P = SemiclassicalJacobi(t, 1, 1, 100)

ensuring that it will not fail. Do keep in mind that it discards the steps if you call it in this way, so if you intend to use the individual steps later, it's still better to call

P = SemiclassicalJacobi.(t, 1, 1, 0:100)

[TSGut]

All seems good here in terms of compatibility with the Q variant in ClassicalOPs. But to make this work right now requires the 0.0.2 branch of SingularIntegrals.jl which is not currently registered, so for now it has to be manually added.

I will now use this implementation to generate the plots for @ioannisPApapadopoulos

[dlfivefifty]

tagged SingularIntegrals v0.0.2

Can you remove the WIP if its ready to be reviewed?

[dlfivefifty]

Shouldn't the neg1c case be superseded by doing a reverse Cholesky?

[TSGut]

Ah, sorry I missed your comments here.

I guess it's a matter of design philosophy for the neg1c stuff. In any case, reverse Cholesky isn't implemented anywhere yet though (just QL) so we can't do that yet. Should be pretty similar to QL in terms of implementation.

There is still a workaround to a bug in this PR (marked by todo comments) that I want to fix (needs changes in other packages) before merging this, unless you want to use this code then temporarily this is fine.

[TSGut]

The PR branch JuliaLinearAlgebra/InfiniteLinearAlgebra.jl#140 will allow me to remove the workaround

[dlfivefifty]

@TSGut if this is ready to merge remove the WIP from the title

[TSGut]

just need to make sure all the dependencies are correct and tests pass first

[TSGut]

@dlfivefifty Ok looking good. Tests pass without the nasty workarounds and all the hierarchy stuff is fast and passes tests. Had to lock to Julia 1.9+ though, hope that's not an issue.

The one obvious thing that is missing is lowering via reverse Cholesky / QL but I think that can be done separately once we have implemented adaptive reverse Cholesky.

This is a huge change to the package obviously, so probably we want to make a "breaking changes" sort of version change. 0.4 maybe? at the moment it's just 0.3.5