Implement finite section based adaptive QL

[TSGut]

Another PR related to JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl#70.

Here, I am working towards implementing infinite dimensional QL based on taking the QL decomposition of finite sections and checking for convergence in the desired block compared to a much larger finite section. The basic functionality works. Some stuff I am currently thinking about:

1) How to expose this functionality, if at all. Overloading ql() would be weird considering the method inherently wants a tolerance and the way the factors are updated does not really lend itself to the MatrixFactorizations type. Perhaps it should stay its own thing that a user only calls when they know what they are doing? I am very open to suggestions on this.

2) sort a bunch of optimizations including getting the MemoryLayout of the adaptive cached L array to be LowerTriangular etc.

[TSGut]

I believe tests are failing bc of an issue possibly (?) unrelated to my additions, see #126

I will have a closer look later.

[dlfivefifty]

There’s no reason to do QL on the full finite section. One can check convergence by checking just a piece of the bottom right

though the details might be tricky

[TSGut]

There’s no reason to do QL on the full finite section.

Do you mean that as in when updating/expanding an existing finite section?

[dlfivefifty]

Right. You can look at Q,L = ql(A[m:n, m:n]). If L[1:b,1:b] is unchanged as n -> ∞ n is sufficiently large. But we can also let m grow with n.

[TSGut]

That would make the update computation a lot cheaper for large n, I will look into it.

[TSGut]

Alright, I think functionality-wise this can do everything I need for the SemiclassicalOPs application we have in mind. @dlfivefifty Any further comments? Also looking for opinions on if and if yes how to expose this to users or to just have people explicitly import this.

Still intend to add some more tests.

[dlfivefifty]

This isn’t a good way to do things. One should follow the Adaptive QR example and have a single data structure to represent the underlying ql(A).factors (which is a banded matrix that contains both the Q and L). Then one can deduce the L via LowerTriangular(ql(A).factors) and the Q similarly.

It is odd in this version that you are using a single type for storing both Q and L, which shouldn’t be constructed at all.

[TSGut]

I can match the QL structure but I will have to see how that works for Q without recomputing it

EDIT: Actually, I think it should be fine.

[TSGut]

So @dlfivefifty just to confirm, you don't want me to mimic QLProduct as implemented in this package here:

InfiniteLinearAlgebra.jl/src/banded/infqltoeplitz.jl

Line 148 in 8b122c8

struct QLProduct{T,QQ<:Tuple,LL} <: Factorization{T}

which seems to be called in the tests when running ql(A) for A an infinite BandedMatrix:

julia> A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5,∞), 1 => Fill(0.5,∞))
ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (1, 1) with data (3-element Vector{Float64}) * (1×ℵ₀ Ones{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(3)×OneToInf() with indices OneToInf()×OneToInf():
 5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
 2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  5.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
 ⋮                        ⋮                        ⋮                        ⋮                   ⋱  

julia> F = ql(A)
InfiniteLinearAlgebra.QLProduct{Float64, Tuple{LowerHessenbergQ{Float64, Fill{MatrixFactorizations.QLPackedQ{Float64, Matrix{Float64}, Vector{Float64}}, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, LowerTriangular{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(hcat), Tuple{Vector{Float64}, ApplyArray{Float64, 2, typeof(*), Tuple{Vector{Float64}, Ones{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}}}}}, InfiniteArrays.OneToInf{Int64}}}}
Q factor:
ℵ₀×ℵ₀ InfiniteLinearAlgebra.ProductQ{Float64, Tuple{LowerHessenbergQ{Float64, Fill{MatrixFactorizations.QLPackedQ{Float64, Matrix{Float64}, Vector{Float64}}, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
  0.99454      -0.104356      0.0           0.0           0.0           0.0          0.0        …  
 -0.103786     -0.98911      -0.104356      0.0           0.0           0.0          0.0           
  0.0108307     0.10322      -0.98911      -0.104356      0.0           0.0          0.0           
 -0.00113025   -0.0107716     0.10322      -0.98911      -0.104356      0.0          0.0           
  0.000117949   0.00112408   -0.0107716     0.10322      -0.98911      -0.104356     0.0           
 -1.23087e-5   -0.000117305   0.00112408   -0.0107716     0.10322      -0.98911     -0.104356   …  
  1.28448e-6    1.22415e-5   -0.000117305   0.00112408   -0.0107716     0.10322     -0.98911       
 -1.34044e-7   -1.27747e-6    1.22415e-5   -0.000117305   0.00112408   -0.0107716    0.10322       
  1.39883e-8    1.33312e-7   -1.27747e-6    1.22415e-5   -0.000117305   0.00112408  -0.0107716     
  ⋮                                                                     ⋮                       ⋱  
L factor:
ℵ₀×ℵ₀ LowerTriangular{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(hcat), Tuple{Vector{Float64}, ApplyArray{Float64, 2, typeof(*), Tuple{Vector{Float64}, Ones{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
  4.76513     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅       …  
 -2.5       -4.79129     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅          
 -0.208712  -2.5       -4.79129     ⋅          ⋅          ⋅          ⋅          ⋅         ⋅          
   ⋅        -0.208712  -2.5       -4.79129     ⋅          ⋅          ⋅          ⋅         ⋅          
   ⋅          ⋅        -0.208712  -2.5       -4.79129     ⋅          ⋅          ⋅         ⋅          
   ⋅          ⋅          ⋅        -0.208712  -2.5       -4.79129     ⋅          ⋅         ⋅       …  
   ⋅          ⋅          ⋅          ⋅        -0.208712  -2.5       -4.79129     ⋅         ⋅          
   ⋅          ⋅          ⋅          ⋅          ⋅        -0.208712  -2.5       -4.79129    ⋅          
   ⋅          ⋅          ⋅          ⋅          ⋅          ⋅        -0.208712  -2.5      -4.79129     
  ⋮                                                      ⋮                                        ⋱  

You want it to instead follow this struct: https://github.com/JuliaLinearAlgebra/MatrixFactorizations.jl/blob/d72eb0626794736c7096f7b65c28494625020153/src/ql.jl#L37

Is that correct? It should be fine either way, just need to know what to work towards

[TSGut]

Okay, now featuring .factors and .τ syntax - probably still some things to tidy up but it all works fine and we can now choose not to construct Q if it isn't needed in the usual QL syntax way. Still not overloading ql() at all for now as I am unsure on what the best way to do that would be (specifically how one would choose which infinite ql to use without a whole bunch of type inference?).

Note that this relies on the current master branch of MatrixFactorizations, so we should version tag that before we merge this.

Intend to still add more tests.

[codecov]

Codecov Report

Patch coverage: 94.78% and project coverage change: +1.88 🎉

Comparison is base (3ae935a) 78.00% compared to head (be7b705) 79.88%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #125      +/-   ##
==========================================
+ Coverage   78.00%   79.88%   +1.88%     
==========================================
  Files          11       11              
  Lines        1314     1407      +93     
==========================================
+ Hits         1025     1124      +99     
+ Misses        289      283       -6     

Impacted Files	Coverage Δ
src/InfiniteLinearAlgebra.jl	92.30% <ø> (ø)
src/infql.jl	62.32% <94.69%> (+15.78%)	⬆️
src/banded/hessenbergq.jl	71.84% <100.00%> (ø)

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

[TSGut]

No idea what happened with CI here but this is doing everything I want now and all tests pass when julia actually builds (?).

@dlfivefifty any further comments or is this good to merge? I am trying to wrap up the semiclassical decomposition stuff and reduce the number of loose branches floating around before I give that the final push

[TSGut]

alright, the CI demons are gone and it all passes.

[dlfivefifty]

What is this type for?

[TSGut]

the factors field should really only be of the adaptive factors type, I will add that explicit declaration then T is used

[TSGut]

unless you mean what is the struct for in general? In that case: this is the main struct that handles the new QL variant. as far as I know without a custom struct how we can tell julia to construct Q and L for us from .factors information when prompted for Q and L or tau without storing redundant data. not sure what the question means in that case basically, this is just the type that implements the adaptive storing and expansion ruleset

Let me know what exactly you meant there, I made the other minor changes

[dlfivefifty]

I meant why do we need AdaptiveQLFiniteSection when it just wraps QLFiniteSectionFactors?

[dlfivefifty]

Also it only works for tridiagonal but doesn't actually error for larger bandwidths

[TSGut]

The tridiagonal part is certainly not intentional. I think it should be an easy fix but I have to check where I assumed tridiagonal (probably in the adjoint lmul? I'll check and push a bugfix soon). It used to work for non tridiagonal so a recent change must have broken it. I'll also add a test to ensure it.

For the degenerate Q that was intentional but I guess it is up for debate whether it should have that behaviour? If we don't want it to do this then we'd have to think about how the algorithm is supposed to detect if it runs into that scenario.

[dlfivefifty]

I think the only possible case of a degenerate QL is the shift operator so it could just check if the first column of Q is zero (which probably has a simple condition on the τ but I'd have to think about it)

[TSGut]

I believe this update fixed the tridiagonal assumption - the test you added passes. Once you have MatrixFactorizations.jl updated we can remove some of the then-redundant overloads here.

Next I will have a look what we can do about detecting the degenerate case via tau. Just to clarify, you think the correct behavior should be to issue an error when this occurs?

[dlfivefifty]

Ideally it would error as it is no long orthogonal, just an isometry

[dlfivefifty]

The tests should "pass" once JuliaArrays/LazyArrays.jl#250 is tagged but you'll notice I found another couple broken cases due to assuming tridiagonal

[TSGut]

Yep, I saw the broken tests. Will investigate

[TSGut]

I guess the problem with multiplication by Q as in the broken tests is that this will result in an infinite vector with all non-zero entries even if the vector we are multiplying with has only finitely many non-zeros.

Alternatives:

1. I can implement it to give a custom cached infinite output vector with it knowing how to adaptively resize. Guess this would require a new struct <:AbstractCachedVector.
2. alternatively we can just return ApplyArray(*,Q,B) which already works (but I assume this is slower since it basically just computes all the entries of Q).

Which do you prefer?

[marcusdavidwebb]

For tridiagonal A the only degenerate case for Q is the shift operator (since Q must be lower Hessenberg and there are limited forms this can take), but for higher bandwidths I'm not sure that there can't be other cases where Q is merely a transposed isometry instead of unitary.

[dlfivefifty]

The best thing to do the same thing as QR and apply the householder reflections until convergence:

InfiniteLinearAlgebra.jl/src/infqr.jl

Line 282 in 3ae935a

function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:AbstractBlockBandedLayout}},<:PaddedLayout}; tolerance=1E-30)

But since Q * b is not really needed for solving systems just leaving it lazy for now is fine

[dlfivefifty]

but for higher bandwidths I'm not sure that there can't be other cases where Q is merely a transposed isometry instead of unitary.

Hmm maybe we just ignore this issue for the time being and return the isometry. We can return to the question how to detect if Q is merely an isometry later. (Probably there's a simple check of linear independence in the first u-1 columns where u is the upper bandwidth but that'll take some thought)

[marcusdavidwebb]

To clarify: The upper triangular shift isn't an isometry, but its transpose is.

Also, now I think about it properly, there might be other possible "transpose isometries" that are Hessenberg... let me think some more

Is it possible to find multiple QL factorisation for a given operator by using this finite section approach?

[marcusdavidwebb]

Consider:

Q =
[1 0 ..... ]
[0 0 1 ... ]
[0 0 0 1 ... ]
[0 0 0 0 1 ... ]
......

This operator satisfies QQ' == I but not Q'Q == I. The tail would always need to be a shift if it is orthogonal and Hessenberg.

[TSGut]

Okay, I replaced multiplication by Q with lazy multiplication. We can always make this more sophisticated in the future if we want to multiply by Q in some application. Tests will continue to fail until LazyArrays is tagged.

Is it possible to find multiple QL factorisation for a given operator by using this finite section approach?

Maybe.. I hadn't thought about it yet. Perhaps one could use an infinite orthogonal matrix U as a left-preconditioner to change what Q we converge to, then apply U' to the resulting Q? Would this necessarily converge to the same Q as just doing QL if there are multiple solutions - maybe.. Do we have a constructed example where we know two Qs exist and what they look like?

[dlfivefifty]

Do we have a constructed example where we know two Qs exist and what they look like?

Yes, any toeplitz operator where the winding number of the symbol is > 0 (or maybe it’s < 0)

[dlfivefifty]

The “branch” optional argument examples on the tests should provide examples

[TSGut]

I will have a look and see if some modification can change what we converge to.

[dlfivefifty]

Deflation?

[dlfivefifty]

But this is a research question so not needed to merge this PR

[TSGut]

The only open question for this PR is how we handle the degenerate case. If we leave it as is, then Julia will incorrectly overload ldiv for the shift operator special case with lmul by Q', causing a silent bug. But I guess we said we leave it like this for now?

[dlfivefifty]

Don't add this dependency

[dlfivefifty]

Don't add this dependency

[TSGut]

Yeah those were just used in dev, not meant to be added