Multiplication of a dense matrix and a symmetric sparse matrix is slow

[MqCreaple]

An example code

A = Symmetric(sprand(20000, 20000, 0.001))
B = rand(100, 20000)
C = B * A

After executing this code, the type of C is a sparse matrix but should be a dense matrix. This results in really slow matrix multiplication in comparison to directly multiplying dense matrix with sparse matrix.

[MqCreaple]

Probably this issue has to be fixed in LinearAlgebra package

[ViralBShah]

Pretty sure it is in the SparseArrays package. Probably just need to pull out the parent and pass it to the multiplication (and ignore the Symmetric nature) as a quick fix.

[KlausC]

pull out the parent and pass it to the multiplication (and ignore the Symmetric nature) as a quick fix

Not seriously??? We won't return incorrect results. In this case it would be optimal to transform B to sparse(B).

[ViralBShah]

I meant convert the Symmetric matrix to a regular matrix and use the fallback.

[KlausC]

To demonstrate the run times:

julia> @time sparse(B) * A
  0.285723 seconds (60 allocations: 165.698 MiB, 8.80% gc time)
100×20000 SparseMatrixCSC{Float64, Int64} with 2000000 stored entries:
⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛

julia> @time B * Matrix(A)
  6.020736 seconds (126.85 k allocations: 3.007 GiB, 0.11% gc time, 1.76% compilation time)
100×20000 Matrix{Float64}:

[ViralBShah]

Right I meant to remove the symmetric tag not the sparsity ...

[KlausC]

You are right, this is even better:

julia> @time B * sparse(A);
  0.052730 seconds (9 allocations: 21.517 MiB)

julia> @time sparse(B) * A;
  0.271246 seconds (62 allocations: 165.721 MiB, 5.44% gc time)

[KlausC]

I propose to change this

https://github.com/JuliaLang/julia/blob/6412a56223e38824ce6eff78bf34662637971e1c/stdlib/LinearAlgebra/src/matmul.jl#L139-L142

method definition of * and make it dispatch to SparseArrays if one of the arguments is sparse (in the sense of issparse). To avoid dependency of SparseArrays, that PR would make use of specialized traits. The alternative of providing a new version of mult! does not work, because the result type shall be dense, not sparse as the current code enforces.

[KlausC]

The proposed change in LinearAlgebra is a precondition for additional changes in SparseArrays in order to get this issue resolved. I am going to start PRs for both libraries, if there are no objections.

[ViralBShah]

Cc @dkarrasch @Wimmerer