Diagonal broadcast results in SparseMatrixCSC

[garrison]

I was surprised to notice that broadcast on a Diagonal matrix results in a sparse matrix, not a dense one:

julia> Diagonal(-1:0) .+ 1
2×2 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1, 1]  =  0
  [2, 1]  =  1
  [1, 2]  =  1
  [2, 2]  =  1

In particular, the 0 in the above example is a "stored entry", so even in the case where there are many zeros in the result, it will still take more memory than a dense matrix.

[kshyatt]

@Sacha0 seems dually-qualified to explain/fix/???? this :)

[martinholters]

The sparse matrix can cope with a broadcast resulting in a dense result, but can also exploit the common situation of the broadcast giving a sparse (in particular: diagonal) result.

[garrison]

the common situation of the broadcast giving a sparse (in particular: diagonal) result.

example, please?

[Sacha0]

example, please?

julia> Diagonal(1:4) .+ Diagonal(1:4)
4×4 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1, 1]  =  2
  [2, 2]  =  4
  [3, 3]  =  6
  [4, 4]  =  8

julia> Diagonal(1:4) .* Bidiagonal(1:4, 1:3, :U)
4×4 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1, 1]  =  1
  [2, 2]  =  4
  [3, 3]  =  9
  [4, 4]  =  16

julia> Diagonal(1:4) .+ sprand(4, 4, 0.2)
4×4 SparseMatrixCSC{Float64,Int64} with 6 stored entries:
  [1, 1]  =  1.0
  [3, 1]  =  0.110989
  [2, 2]  =  2.88061
  [3, 3]  =  3.0
  [3, 4]  =  0.0425649
  [4, 4]  =  4.0

For the design discussion in all its detail, please see #18590 and related threads. Where the result happens to be dense in sparse form, we may yet avoid storing numerical zeros. But otherwise this behavior is expected :). Best!

[garrison]

To give some examples, here's what I might have naively expected before reading #18590:

• Diagonal(1:4) .+ Diagonal(5:8) -> Diagonal{Int}
• Diagonal(1:4) .* Diagonal(5:8) -> Diagonal{Int}
• Diagonal(1:4) .+ 1 -> Matrix{Int}
• Diagonal(1:4) .* 2 -> Diagonal{Int}
• mytimes(x,y) = x * y; mytimes.(Diagonal(1:4), 2) -> SparseMatrixCSC{Int,Int}

but after reading #18590 (comment) I no longer believe the effort it would take to make broadcasting on special matrices always return the "optimal" type is worth it.

With this, I believe there are two questions that remain:

1. Should a sentence or two be added to the Broadcasting section of the manual mentioning the reasoning for returning a sparse matrix for all these operations?
2. Should dropzeros! be applied automatically to the results of these operations?

[Sacha0]

Should a sentence or two be added to the Broadcasting section of the manual mentioning the reasoning for returning a sparse matrix for all these operations?

Sounds great! :)

Should dropzeros! be applied automatically to the results of these operations?

Do you mean specifically to the results of operations that we cannot presently determine are zero-preserving? If so, yes, chances are we should do that implicitly in the future; a couple existing issues discuss this possibility. (Operations that we can presently determine are zero-preserving do this implicitly already.) Best!

[oscardssmith]

I really think this path is a bad one. If our only justification (from #18590) for breaking with user expectations in our API design is "We didn't feel like writing some trivial methods", that suggests that the API design is wrong.

[garrison]

I'm going to close this. It will be great if a sentence is added to the Broadcasting section of the manual explaining the reasoning here, but this ticket currently serves as a poor reminder to do so.