BandedMatrices.jl
BandedMatrices.jl copied to clipboard
Feature request: convenience functions to construct types equivalent to those from LinearAlgebra
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.
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).