Support Jacobi matrix computation via Cholesky and QR #120

TSGut · 2023-01-23T23:07:09Z

After the unnecessary episode of the premature merger, here is the current status of this PR. Still a lot left to sort out but the functionality is there and works. The old PR with chat history is here.

The big bug to fix for this PR is here: JuliaLinearAlgebra/InfiniteLinearAlgebra.jl#123
Due to that bug I currently have to rely on a very adhoc fix for loss of positive definiteness.

This branch is intended to be used with JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl#70, but both are still WIP.

Get PR uptodate

codecov · 2023-01-25T23:32:54Z

Codecov Report

Patch coverage: 97.97% and project coverage change: +0.96 🎉

Comparison is base (60abdde) 88.74% compared to head (4496079) 89.70%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #120      +/-   ##
==========================================
+ Coverage   88.74%   89.70%   +0.96%     
==========================================
  Files          16       17       +1     
  Lines        1652     1749      +97     
==========================================
+ Hits         1466     1569     +103     
+ Misses        186      180       -6

Impacted Files	Coverage Δ
src/ClassicalOrthogonalPolynomials.jl	`86.60% <0.00%> (ø)`
src/choleskyQR.jl	`100.00% <100.00%> (ø)`

... and 1 file with indirect coverage changes

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

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

get up to date

src/choleskyQR.jl

TSGut · 2023-04-17T15:58:34Z

ok @dlfivefifty, this should be better. let me know if you spot anything else. Tests are passing and I am pretty much happy with it as is for now.

I think we also want to implement sqrtweighted but I would rather do that as a separate PR once this is merged

src/choleskyQR.jl

*) split storage 2xN matrix into 2 vectors. *) remove an instance of Float64 hardcoding *) remove an unnecessary infinite tail of zeros *) pre-fill cached arrays to avoid step-by-step expansion

TSGut · 2023-04-19T05:41:18Z

@dlfivefifty: Another round of minor improvements, with band storage now styled after the band storage in SemiclassicalOPs, i.e. like in https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/blob/5e4d53d5d875d54a8b459839237a7c537d9d859d/src/SemiclassicalOrthogonalPolynomials.jl#L110

Also some minor typing improvements here and there.

Not sure how to share profiler results in a sensible manner here since I mostly rely on ProfileView.jl, but here are some benchmarks on my laptop:

julia> @benchmark qr_jacobimatrix(sqrtW, P, false)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  16.051 μs …   8.576 ms  ┊ GC (min … max):  0.00% … 99.45%
 Time  (median):     17.477 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   21.229 μs ± 136.794 μs  ┊ GC (mean ± σ):  11.10% ±  1.72%

    ▄█▆▃                                                        
  ▁▅████▇▅▄▃▃▂▃▂▃▃▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  16.1 μs         Histogram: frequency by time         29.5 μs <

 Memory estimate: 28.17 KiB, allocs estimate: 296.

julia> @benchmark qr_jacobimatrix(sqrtW, P, false)[1000,1000]
BenchmarkTools.Trial: 3487 samples with 1 evaluation.
 Range (min … max):  1.212 ms …  10.888 ms  ┊ GC (min … max): 0.00% … 84.50%
 Time  (median):     1.328 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.430 ms ± 884.072 μs  ┊ GC (mean ± σ):  6.12% ±  8.55%

       ▃▃ ▂▁▁█▅▂▅▁                                             
  ▂▂▂▂▄███████████▅▅▅▅▅▄▃▃▄▄▄▄▃▃▃▂▃▃▂▂▂▂▂▂▂▁▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▂▂ ▃
  1.21 ms         Histogram: frequency by time        1.77 ms <

 Memory estimate: 914.77 KiB, allocs estimate: 11400.

Compared with Lanczos for the same problem:

julia> @benchmark jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(legendre(0..1))))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  290.774 μs …  15.843 ms  ┊ GC (min … max): 0.00% … 90.28%
 Time  (median):     311.731 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   331.059 μs ± 409.251 μs  ┊ GC (mean ± σ):  3.26% ±  2.60%

      ▂▆███▅▅▄▂▂                                                 
  ▁▁▃▆███████████▇▆▅▅▅▄▄▄▃▃▃▃▃▃▃▃▂▂▃▃▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  291 μs           Histogram: frequency by time          404 μs <

 Memory estimate: 97.06 KiB, allocs estimate: 833.

julia> @benchmark jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(legendre(0..1))))[1000,1000]
BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  1.907 s …    3.121 s  ┊ GC (min … max): 1.39% … 1.56%
 Time  (median):     2.104 s               ┊ GC (median):    2.02%
 Time  (mean ± σ):   2.377 s ± 651.585 ms  ┊ GC (mean ± σ):  1.65% ± 0.33%

  █        █                                               █  
  █▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.91 s         Histogram: frequency by time         3.12 s <

 Memory estimate: 307.94 MiB, allocs estimate: 12968641.

It looks good to me but if there are still allocations you don't expect let me know and I can have another look.

dlfivefifty · 2023-04-19T07:32:58Z

OK if its faster than the existing Lanczos then its fine, don't worry about making it faster.

TSGut · 2023-04-19T08:42:40Z

Yes, I have been benchmarking against Lanczos and it is significantly faster.

I think this is ready to be merged if you are happy with the changes. Then I can focus on wrapping up the PR in SemiclassicalOPs as well.

dlfivefifty · 2023-04-19T09:44:01Z

Don't we want a CholeskyPolynomial to replace LanczosPolynomial that uses this?

Or perhaps we just need a more general type:

struct GeneralOrthogonalPolynomial{T, WW, XX} <: OrthogonalPolynomial{T}
    weight::WW
    X::XX
end

TSGut · 2023-04-19T10:11:40Z

I think the simplicity of the three weight terms in the hierarchy of semiclassical Jacobi means we still want to implement that separately like I have in the PR in that repo since a lot of simplifications happen there for computing the hierarchy step by step.

But in general yes I don't see why we can't have a GeneralOrthogonalPolynomial. Raises some design philosophy questions though, since I don't know how to generically "detect" that an arbitrary supplied weight function is the square of a polynomial, so this wouldn't be able to leverage the QR variant? We always could do CholeskyOrthogonalPolynomial and QROrthogonalPolynomial separately but the latter isn't really a nice name. 😅

I think we should definitely keep the Lanzcos approach around, though. It's an extremely helpful sanity check on these things.

TSGut · 2023-04-19T11:29:28Z

BTW, it would be nice to tag a new version of this for use in SemiclassicalOPs.

dlfivefifty · 2023-04-19T14:07:52Z

src/choleskyQR.jl

+end
+
+# The generated Jacobi operators are symmetric tridiagonal, so we store their data in cached bands
+mutable struct CholeskyJacobiBands{dv,T} <: AbstractCachedVector{T}


Wait I think we should rename this CholeskyJacobiBand

Should take no time at all, should I just open a new PR?

I'll just do it

dlfivefifty · 2023-04-19T14:09:39Z

src/choleskyQR.jl

+
+An optional bool can be supplied, i.e. `cholesky_jacobimatrix(sqrtw, P, false)` to disable checks of symmetry for the weight multiplication matrix and orthonormality for the basis (use with caution).
+"""
+function cholesky_jacobimatrix(w::Function, P::OrthogonalPolynomial, checks::Bool = true)


Need to remove ::OrthogonalPolyomial as not all orthogonal polynomials are <: OrthogonalPolynomial, eg, any mapped polynomials

Is there a sensible replacement? I guess I can just hope the user only gives it sensible input? I can try to compute P's jacobi matrix and then catch an error to return something like "Failed to compute Jacobi matrix, check the supplied Polynomials" or something?

Who cares if the user doesn't give sensible input? It will just error anyways so why even check??

And who's "the user"??? Me??

We are the user (for now but others are using this package). 😄 It's nice to get sensible error messages but tbf something like "no method matching jacobimatrix(::Type)" will be fine even for those standards and that will be automatic

dlfivefifty · 2023-04-19T14:17:09Z

You could in theory just check if it is of the form w .^ 2 but there's no reason to add that. I think the interface would be:

OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthonormalpolynomial(singularities(w)))
OrthogonalPolynomial(w::AbstractQuasiVector, P::AbstractQuasiMatrix) = GeneralOrthogonalPolynomial(w, cholesky_jacobimatrix(w, P))

But this does have a catch: we lose reference to the P which we know we can convert to-from.

So perhaps a better option is something like:

# represent an Orthogonal polynomial which has a conversion operator to P
struct ConvertedOrthogonalPolynomial{T, WW, XX, PP} <: OrthogonalPolynomial{T}
    weight::WW
    X::XX
    P::PP
end

abstract type AbstractConvertedOPLayout <: AbstractOPLayout end
struct ConvertedOPLayout <: AbstractConvertedOPLayout end
# also change NormalizedOPLayout <: ConvertedOPLayout

# transform to P * U if needed for differentiation, etc.
arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, Q.X.U

# also change all the NormalizedOPLayout
@inline copy(L::Ldiv{Lay,<: AbstractConvertedOPLayout}) where Lay = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))

Maybe I'll do the PR for this

dlfivefifty · 2023-04-19T14:19:37Z

I think we should definitely keep the Lanzcos approach around, though

In general better to delete unused code: keeping it around just creates admin, and Git will always keep the code around. Though in this case we are using in the tests.

dlfivefifty · 2023-04-19T14:29:37Z

BTW, it would be nice to tag a new version of this for use in SemiclassicalOPs.

Let me fix it first.

TSGut · 2023-04-19T14:30:40Z

Maybe I'll do the PR for this

It's up to you! Just assign an issue to me if you want me to do something related to this. I assume that whoever does this new PR will also make the small changes mentioned here, i.e. change "Bands" to "Band" in the QR and Cholesky case and remove the ::OrthogonalPolynomial type from the function calls?

But this does have a catch: we lose reference to the P which we know we can convert to-from.

Remembering this is probably better. we have to carry the conversion operator around either way for expanding the cached values so imo it makes sense to remember the original basis somewhere

TSGut added 4 commits January 23, 2023 11:09

Merge pull request #1 from JuliaApproximation/main

ccea7f9

Get PR uptodate

workaround for bug

21bd437

Update README.md

ef8c485

Merge remote-tracking branch 'origin/main' into main

ebd9e47

TSGut mentioned this pull request Jan 23, 2023

WIP: Replace Lanzos with Cholesky JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl#70

Merged

TSGut added 4 commits January 25, 2023 17:28

begin also implementing QR variant

4310453

rm wrong type test

b60c890

rename files

ad2ebdb

more name changing

eb0fb68

TSGut added 7 commits January 25, 2023 23:53

minor changes

2c7eeb8

Merge branch 'main' into main

74334d1

Merge pull request #7 from JuliaApproximation/main

b8f8118

get up to date

minor changes to tests

0c2db8d

make some safety checks optional via bool

5a64f79

Merge branch 'main' into main

a8a747b

fix minor typo in main branch comments

10088b6

dlfivefifty requested changes Apr 17, 2023

View reviewed changes

src/choleskyQR.jl Outdated Show resolved Hide resolved

src/choleskyQR.jl Outdated Show resolved Hide resolved

src/choleskyQR.jl Outdated Show resolved Hide resolved

TSGut added 6 commits April 17, 2023 15:14

change entry computation + remove redundant field

c74ef5e

clarify a comment

d78f256

add more tests

d0cf141

Merge remote-tracking branch 'origin/main' into main

ce99a40

use known lengths instead of padding

51cff6d

remove redundant line in test

b5d622d

TSGut changed the title ~~WIP: Cholesky Polynomials revival~~ Support Jacobi matrix computation via Cholesky and QR Apr 17, 2023

TSGut requested a review from dlfivefifty April 18, 2023 14:22

dlfivefifty requested changes Apr 18, 2023

View reviewed changes

TSGut added 2 commits April 19, 2023 05:50

lots of small changes, see description

7269e53

*) split storage 2xN matrix into 2 vectors. *) remove an instance of Float64 hardcoding *) remove an unnecessary infinite tail of zeros *) pre-fill cached arrays to avoid step-by-step expansion

minor type adjustments

6693844

remove now unnecessary methods

f2c7d16

simplify resizing since we are using vectors now

4496079

dlfivefifty approved these changes Apr 19, 2023

View reviewed changes

dlfivefifty merged commit a424701 into JuliaApproximation:main Apr 19, 2023

dlfivefifty reviewed Apr 19, 2023

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Support Jacobi matrix computation via Cholesky and QR #120

Support Jacobi matrix computation via Cholesky and QR #120

TSGut commented Jan 23, 2023 •

edited

Loading

codecov bot commented Jan 25, 2023 •

edited

Loading

TSGut commented Apr 17, 2023

TSGut commented Apr 19, 2023 •

edited

Loading

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023

TSGut commented Apr 19, 2023

dlfivefifty Apr 19, 2023

TSGut Apr 19, 2023

dlfivefifty Apr 19, 2023

dlfivefifty Apr 19, 2023

TSGut Apr 19, 2023

dlfivefifty Apr 19, 2023

dlfivefifty Apr 19, 2023

TSGut Apr 19, 2023 •

edited

Loading

dlfivefifty commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023 •

edited

Loading

Support Jacobi matrix computation via Cholesky and QR #120

Support Jacobi matrix computation via Cholesky and QR #120

Conversation

TSGut commented Jan 23, 2023 • edited Loading

codecov bot commented Jan 25, 2023 • edited Loading

Codecov Report

TSGut commented Apr 17, 2023

TSGut commented Apr 19, 2023 • edited Loading

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023

TSGut commented Apr 19, 2023

dlfivefifty Apr 19, 2023

Choose a reason for hiding this comment

TSGut Apr 19, 2023

Choose a reason for hiding this comment

dlfivefifty Apr 19, 2023

Choose a reason for hiding this comment

dlfivefifty Apr 19, 2023

Choose a reason for hiding this comment

TSGut Apr 19, 2023

Choose a reason for hiding this comment

dlfivefifty Apr 19, 2023

Choose a reason for hiding this comment

dlfivefifty Apr 19, 2023

Choose a reason for hiding this comment

TSGut Apr 19, 2023 • edited Loading

Choose a reason for hiding this comment

dlfivefifty commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

dlfivefifty commented Apr 19, 2023

TSGut commented Apr 19, 2023 • edited Loading

TSGut commented Jan 23, 2023 •

edited

Loading

codecov bot commented Jan 25, 2023 •

edited

Loading

TSGut commented Apr 19, 2023 •

edited

Loading

TSGut Apr 19, 2023 •

edited

Loading

TSGut commented Apr 19, 2023 •

edited

Loading