Lowering/Raising for TwoBand

[ioannisPApapadopoulos]

I want to implement raising and lowering for TwoBand. I am able to get a hold of the necessary coefficient but I am not sure how to put them into the correct matrix format. E.g.

using SemiclassicalOrthogonalPolynomials

P, Q = TwoBandJacobi(ρ,0,0,0),TwoBandJacobi(ρ,1,1,0)
A = Q.P \ P.P   # Coefficients for even degree
B = Q.Q \ P.Q   # Coefficients for odd degree

# Want to interlace A and B together

The issue is actually getting ahold of the diagonal and the super diagonals is proving tricky. A and B are made from a series of nested calls, e.g.A.args[1].args[1].args[2].dv gives some diagonal. But to form A that's multiplied by other matrices and coefficients. Is there a smart way to interlace A and B? @dlfivefifty ?

[TSGut]

If I am understanding your intent correctly then by chance I actually just needed this exact type of construction when working on the other issue you opened. I made a non-cached version that allows one to select bands from pretty much any type of infinite array, even if the data structure does not make it "nice" to do so:

using InfiniteArrays
import Base: size, getindex
mutable struct SelectInfiniteBand{T} <: AbstractVector{T}
    data::AbstractMatrix{T}
    band::Integer
    SelectInfiniteBand{T}(M::AbstractMatrix{T}, band::Integer) where T = new{T}(M, band)
end
function SelectInfiniteBand(M::AbstractMatrix{T}, band::Integer) where T
    # only for infinite square matrices
    @assert length(axes(M,1)) == length(axes(M,2))
    @assert length(axes(M,1)) == ℵ₀
    return SelectInfiniteBand{T}(M, band)
end
Base.size(::SelectInfiniteBand) = (ℵ₀,)
function Base.getindex(M::SelectInfiniteBand, nm::Int)
    if M.band <= 0
        return M.data[nm-M.band,nm]
    else
        return M.data[nm,nm+M.band]
    end
end

Some examples:

julia> using LinearAlgebra, LazyArrays

julia> X = ApplyArray(*,jacobimatrix(Legendre()) + π*I,Diagonal(Fill(2,∞)))
(ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ Diagonal{Int64, Fill{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 6.28319  0.666667   ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        …  
 2.0      6.28319   0.8       ⋅         ⋅         ⋅         ⋅         ⋅         ⋅           
  ⋅       1.33333   6.28319  0.857143   ⋅         ⋅         ⋅         ⋅         ⋅           
  ⋅        ⋅        1.2      6.28319   0.888889   ⋅         ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅       1.14286   6.28319   0.909091   ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅        ⋅        1.11111   6.28319   0.923077   ⋅         ⋅        …  
  ⋅        ⋅         ⋅        ⋅         ⋅        1.09091   6.28319   0.933333   ⋅           
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅        1.07692   6.28319   0.941176     
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅        1.06667   6.28319      
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅         ⋅        1.05882      
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅           
 ⋮                                               ⋮                                       ⋱  

julia> SelectInfiniteBand(X,0)
ℵ₀-element SelectInfiniteBand{Float64} with indices OneToInf():
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 6.283185307179586
 ⋮

julia> SelectInfiniteBand(X,1)
ℵ₀-element SelectInfiniteBand{Float64} with indices OneToInf():
 0.6666666666666666
 0.8
 0.8571428571428571
 0.8888888888888888
 0.9090909090909091
 0.9230769230769231
 0.9333333333333333
 0.9411764705882353
 0.9473684210526315
 0.9523809523809523
 0.9565217391304348
 0.96
 0.9629629629629629
 0.9655172413793104
 ⋮

julia> SelectInfiniteBand(X,-1)
ℵ₀-element SelectInfiniteBand{Float64} with indices OneToInf():
 2.0
 1.3333333333333333
 1.2
 1.1428571428571428
 1.1111111111111112
 1.0909090909090908
 1.0769230769230769
 1.0666666666666667
 1.0588235294117647
 1.0526315789473684
 1.0476190476190477
 1.0434782608695652
 1.04
 1.037037037037037
 ⋮

It would be straightforward to implement a cached version too if we want it. Let me know if you find any bugs with this. Might make a PR of this into... I don't know where it would belong @dlfivefifty - maybe InfiniteLinearAlgebra? Using something like this we could also overload the ArrayLayouts.jl functions diagonaldata() and supdiagonaldata() etc. for all infinite arrays which would make a lot of things easier for that other high c issue that I am working on.

[dlfivefifty]

Why don't you just use cache(view(X,band(0)))?

[dlfivefifty]

And I don't think want to overload those functions. Those are meant for TridiagonalLayout and doing what you propose will be very slow.

What you probably want is just cache(X).

[TSGut]

Why don't you just use cache(view(X,band(0)))?

That would be because I didn't know the functionality existed and I think it isn't exported by default by BandedMatrices.

For @ioannisPApapadopoulos, what Sheehan said seems to live here: https://github.com/JuliaLinearAlgebra/BandedMatrices.jl/blob/afdd8a24c0a6d9562314e4d2ad5e3c9c6f0a92da/src/generic/Band.jl

Thanks for the pointer, that simplifies things! In what sense is it slow, though? Both what I did and your version take nanoseconds to complete lazily for infinite matrices. Obviously how long it takes to fill for large indices will depend on the matrix structure and can be slow but a user would likely be aware of that? In any case, not a big concern considering I can just replace diagonaldata(M) with view(M,band(0)),

[TSGut]

Nevermind actually, it is explicitly exported by BandedMatrices.jl, I just didn't have it as an explicit dep in the repo I was testing in.

[ioannisPApapadopoulos]

Exactly what I need, thank you both!