Products and inverses of symmetric banded matrices

[miromarszal]

Symmetric banded matrices seem very fragile and fall back to a generic dense matrix on occasions. Consider the following:

B = BandedMatrix(Ones(5,5), (0,0))
B * B

5×5 BandedMatrix{Float64} with bandwidths (0, 0):
 1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅    ⋅   1.0

Product of banded matrices is banded, as expected. However, when the matrix is additionally specified as Symmetric, a product turns out to be dense:

S = Symmetric(B)
S * S

5×5 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

Moreover, inverse of a symmetric matrix should be symmetric. This works with regular matrices, but when it is a BandedMatrix, the result is once again a generic matrix:

inv(S)

5×5 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

EDIT: I discovered more issues:

• Symmetric{BandedMatrix} x Diagonal becomes dense, while BandedMatrix x Diagonal remains BandedMatrix.
• Symmetric{BandedMatrix} x Symmetric{Matrix} throws an error: MethodError: *(::Symmetric{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}, ::Symmetric{Float64, Matrix{Float64}}) is ambiguous.

[jishnub]

S * S seems to call https://github.com/JuliaMatrices/ArrayLayouts.jl/blob/ade9957631d9153734938bfe4ccf42e38fb7b3be/src/muladd.jl#L41, which materializes the array. Perhaps similar may be extended here so that the output is banded.

[dlfivefifty]

I think probably we have a few options:

1. Overload similar(::MulAdd{<:SymmetricLayout{<:AbstractBandedLayout},<:SymmetricLayout{<:AbstractBandedLayout}}) to create a banded matrix.
2. Overload copy(::Mul{<:SymmetricLayout{<:AbstractBandedLayout},<:SymmetricLayout{<:AbstractBandedLayout}) (or possibly mulreduce) to convert arguments to banded matrices first.