Lanczos -> Cholesky - Githubissues

TSGut commented 1 year ago

Going to try and finish this now since a lot of the dependencies have been merged recently. Continued from https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/70

TSGut commented 1 year ago

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[bot] commented 1 year ago

Codecov Report

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

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

Additional details and impacted files

```diff @@ 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](https://app.codecov.io/gh/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/75?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation) | Coverage Δ | | |---|---|---| | [src/SemiclassicalOrthogonalPolynomials.jl](https://app.codecov.io/gh/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/75?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation#diff-c3JjL1NlbWljbGFzc2ljYWxPcnRob2dvbmFsUG9seW5vbWlhbHMuamw=) | `90.87% <100.00%> (+2.59%)` | :arrow_up: | | [src/neg1c.jl](https://app.codecov.io/gh/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/75?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation#diff-c3JjL25lZzFjLmps) | `100.00% <100.00%> (ø)` | |

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

TSGut commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

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 commented 1 year ago

tagged SingularIntegrals v0.0.2

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

dlfivefifty commented 1 year ago

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

TSGut commented 1 year ago

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 commented 1 year ago

The PR branch https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/pull/140 will allow me to remove the workaround

dlfivefifty commented 1 year ago

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

TSGut commented 1 year ago

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

TSGut commented 1 year ago

@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

JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

Lanczos -> Cholesky #75

Codecov Report