Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Diagonality and bandedness checks with InfiniteArrays.jl? #56

Open
TSGut opened this issue Jan 6, 2021 · 5 comments
Open

Diagonality and bandedness checks with InfiniteArrays.jl? #56

TSGut opened this issue Jan 6, 2021 · 5 comments

Comments

@TSGut
Copy link
Contributor

TSGut commented Jan 6, 2021

The methods in InfiniteArrays.jl and InfiniteLinearAlgebra.jl don't currently seem to agree on what is a diagonal and what is a banded matrix. To give some concrete examples:

A minimal working example for Diagonal is

julia> LinearAlgebra.isdiag(InfiniteArrays.Diagonal(0:∞))
true

but the following never finishes computing, i.e. does not even return false:

julia> LinearAlgebra.isdiag(InfiniteLinearAlgebra.BandedMatrix(0=>(0:∞)))

Furthermore, we have a similar situation for banded matrix types:

julia> InfiniteLinearAlgebra.BandedMatrix(0 => (1:∞)) isa InfiniteLinearAlgebra.BandedMatrix
true

julia> InfiniteArrays.Diagonal(1:∞) isa InfiniteLinearAlgebra.BandedMatrix
false

What one wants these to evaluate to might be up to the design philosophy between the packages to some degree but the current behavior does seem odd.

@dlfivefifty
Copy link
Member

dlfivefifty commented Jan 6, 2021

julia> LinearAlgebra.isdiag(InfiniteLinearAlgebra.BandedMatrix(0=>(0:∞)))

This should be added to BandedMatrices.jl, but note we need to check off-diagonal bands too:

julia> isdiag(zeros(5,5))
true

something like

function isdiag(A::AbstractBandedMatrix)
   bandwidths(A) == (0,0) && return true
   Base.invoke(isdiag, Tuple{AbstractMatrix}, A)
end
julia> InfiniteArrays.Diagonal(1:∞) isa InfiniteLinearAlgebra.BandedMatrix
false 

This isn't changeble since Diagonal is a StdLib type.

@dlfivefifty
Copy link
Member

Actually, note isdiag calls isbanded which calls istriu and istril so its really the last two that need overloading. So we can modify the code from StdLib, say in ArrayLayouts.jl. Something like:

function istriu(A::LayoutMatrix, k::Integer = 0)
    require_one_based_indexing(A)
    m, n = size(A)
    for j in 1:min(n, m + k - 1)
        for i in (max(1, j - k + 1):m)  colsupport(A,j)
            iszero(A[i, j]) || return false
        end
    end
    return true
end

and similar for istril.

@dlfivefifty
Copy link
Member

@TSGut do you want to have a go making a PR for this?

@TSGut
Copy link
Contributor Author

TSGut commented Jan 7, 2021

Sure. You said you think it would fit into BandedMatrices.jl?

EDIT: Actually just saw you also mentioned ArrayLayouts.jl. I suppose I'll try to get it working and passing tests for now. I'll submit the PR to the package you think it would most fit into, so just let me know.

@dlfivefifty
Copy link
Member

Hmm actually you could just put it in BandedMatrices.jl since we can do it a bit smarter by adding:

function istriu(A::AbstractBandedMatrix, k::Integer)
    l,_ = bandwidths(A)
    -l  k && return true
    ...
end

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants