BandedMatrices.jl icon indicating copy to clipboard operation
BandedMatrices.jl copied to clipboard

Feature request: convenience functions to construct types equivalent to those from LinearAlgebra

Open jishnub opened this issue 3 years ago • 1 comments

Proposal:

Add functions diagonal, tridiagonal etc that construct appropriate banded matrices.

Eg:

julia> diagonal(d) = BandedMatrix(0=>d);

julia> diagonal(ones(4))
4×4 BandedMatrix{Float64} with bandwidths (0, 0):
 1.0   ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0

julia> tridiagonal(l, d, u) = BandedMatrix(-1=>l, 0=>d, 1=>u);

julia> tridiagonal(ones(3), -2ones(4), ones(3))
4×4 BandedMatrix{Float64} with bandwidths (1, 1):
 -2.0   1.0    ⋅     ⋅ 
  1.0  -2.0   1.0    ⋅ 
   ⋅    1.0  -2.0   1.0
   ⋅     ⋅    1.0  -2.0

Reasoning:

Many of the structured matrices in LinearAlgebra may be represented as BandedMatrices. The advantage of a BandedMatrix type is that the type information is preserved under certain multiplications:

julia> BT = BandedMatrix(1=>ones(3), -1=>-ones(3))
4×4 BandedMatrix{Float64} with bandwidths (1, 1):
  0.0   1.0    ⋅    ⋅ 
 -1.0   0.0   1.0   ⋅ 
   ⋅   -1.0   0.0  1.0
   ⋅     ⋅   -1.0  0.0

julia> BT * BT
4×4 BandedMatrix{Float64} with bandwidths (2, 2):
 -1.0   0.0   1.0    ⋅ 
  0.0  -2.0   0.0   1.0
  1.0   0.0  -2.0   0.0
   ⋅    1.0   0.0  -1.0

julia> T = Tridiagonal(-ones(3), zeros(4), ones(3))
4×4 Tridiagonal{Float64, Vector{Float64}}:
  0.0   1.0    ⋅    ⋅ 
 -1.0   0.0   1.0   ⋅ 
   ⋅   -1.0   0.0  1.0
   ⋅     ⋅   -1.0  0.0

julia> T * T
4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries:
 -1.0    ⋅    1.0    ⋅ 
   ⋅   -2.0    ⋅    1.0
  1.0    ⋅   -2.0    ⋅ 
   ⋅    1.0    ⋅   -1.0

The structured matrix types in LinearAlgebra lose the banded-ness on multiplications, whereas a BandedMatrix preserves this information. As a consequence, in certain applications, it might make sense to have an easy way to construct a drop-in replacement for structured matrix types from LinearAlgebra. This would help a user to easily switch between the two packages without having to change the syntax.

jishnub avatar Jan 13 '22 10:01 jishnub

I don't like the proposed names. In fact we already have ArrayLayouts.diagonal (which is used to support both Diagonal and QuasiDiagonal in QuasiArrays.jl).

dlfivefifty avatar Jan 13 '22 10:01 dlfivefifty