JuliaSparse
sparse \ sparse is dense
The expression "sparse \ sparse" returns a dense matrix, while the equivalent "sparse / sparse" returns a sparse matrix:
julia> using LinearAlgebra
julia> using SparseArrays
julia> A = sprandn(2,2, 0.5) + 2I
julia> B = sprandn(2,2, 0.5) + 2I
julia> A / B
2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
1.91532 ⋅
0.411181 0.560145
julia> A \ B
2×2 Matrix{Float64}:
0.522107 0.0
-0.112086 1.78525
Note that A \ B is the same as (B' / A')', i.e. there is a straightforward way to calculate the respective sparse result.
I think I was hitting a special case here, and that sparse \ sparse is not supported in general.
In general the inverse of a sparse matrix is dense, so defaulting to a dense matrix (if we support this operation at all) eems like the right thing…
The issue is that I saw that A/B was sparse, while B\A was dense. This doesn't make sense, these two operations should return equivalent types.
However, the operation isn't supported by Julia except in special cases, so the point is moot.
Generally, I want to argue that "sparse * sparse = sparse", and thus "sparse / sparse = sparse" makes sense in many cases. It is an operation that should be supported in some way. (The case is different from inv.)
"sparse * sparse = sparse", and thus "sparse / sparse = sparse"
I don't think this follows.
sparse * sparse = sparse is not true in general. It will be fully dense in a lot of problems, and have pretty significant fill in for a lot of others. Tim Davis uses a bitmap format for the latter case, and it's far more common to see that in results than I thought at first.
I'll second on that. This is one of the reasons why it's very much worth considering representing products of sparse matrices lazily as, for instance, in LinearMaps.jl.
I meant to say that "sparse * sparse = sparse" in many cases that are interesting to me (and hopefully to others), and similarly for "sparse / sparse". Apologies for the bad notation. (For example, the matrices might be diagonal except for small dense blocks near the beginning/end of the diagonal.) That is, I would like to be able to perform these operations in some (once Julia supports them). This doesn't have to be the default.
Anyway, since this is currently not supported the point is currently moot. Apologies for the confusion.
It's an interesting problem though. I'll have to look at what LinearMaps does, but there's enough different possible result types that union splitting won't be enough if we really optimized sparse products as much as possible.