Operations with Eye do not preserve sparsity

[amilsted]

Hi, I've noticed that SquareEye, since it is a LinearAlgebra.Diagonal preserves sparsity in many operations. For example:

julia> using FillArrays

julia> using SparseArrays

julia> E = Eye(2)
2×2 Eye{Float64}

julia> M = sprand(2,2, 0.5)
2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 0.906128   0.408236
 0.0285949   ⋅ 

julia> E + M
2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.90613    0.408236
 0.0285949  1.0

julia> kron(E, E)  # although it would be nice is this were another SquareEye!
4×4 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0

However, general Eye typically converts to a dense matrix:

julia> E = Eye(2,2)
2×2 Eye{Float64}

julia> E + M
2×2 Matrix{Float64}:
 1.90613    0.408236
 0.0285949  1.0

julia> kron(E,E)
4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

This was quite surprising to me. Could we add some operator overloads to keep things sparse in the Eye case?

[jishnub]

The broadcast is actually a bit involved, and defining everything here might not be the best idea. The interplay between the structured matrix broadcast style and the sparse matrix style is handled by SparseArrays, so perhaps that is the package that should provide a hook to access that functionality for custom structured matrix types.

Regarding kron, note that the result is not an Eye in general for rectangular matrices, although perhaps we may express this as a sparse matrix

julia> E = Eye(2,3)
2×3 Eye{Float64}

julia> kron(E, E)
4×9 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0

The Diagonal case should be fixed in #268

[amilsted]

Right, I realize you can't make kron of rectangular Eye still be an Eye, but that result is very sparse, so having kron(Eye(m,n), Eye(x,y)) return a sparse matrix seems reasonable. I don't think this case can be handles in SparseArrays.

I guess another alternative is to return a SciMLOperators lazy tensor product...

[jishnub]

Now the rectangular kron is also addressed in #268, and it'll return a sparse array

[putianyi889]

getindex also doesn't preserve sparsity

julia> Eye(5)[:,:]
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

I'm not sure if it's worth adding BandedMatrices to the dependency.

[dlfivefifty]

Adding banded matrices as a dependency won’t work

we could add a new type called OffDiagonal

[dlfivefifty]

Note layout_getindex(D, :, :) will produce a BandedMatrix