Specialize multiplication of AbstractQ with sparse arrays

[dkarrasch]

It turns out that in a few packages (in particular optimization-related ones) there is multiplication of AbstractQ instances and SparseArrays. Due to a lack of specialization, this falls back to _generic_matmatmul!, which is horrible for at least three reasons: (1) it's computing the elements of the matrix representation of Q one-by-one; (2) it doesn't use sparsity at all; and (3) it writes the dense result into an empty sparse output array. I have a hard time imagining anything worse than that. This popped up once we drop AbstractQ <: AbstractMatrix subtyping, see JuliaLang/julia#46196.

[dkarrasch]

Any opinions? Since this occurs in the wild, I guess there is no option to disallow that multiplication?

[ViralBShah]

Won't this be bad for large-ish sparse matrices, where conversion to dense will run out of memory?

[dkarrasch]

That conversion seems still beneficial in that case, because otherwise you'll get dense SparseMatrixCSC, which requires the values and the index vectors to be stored. Note we don't currently have multiplication with sparse matrices at all, so that is covered by the most generic method in LinearAlgebra, not here.

[dkarrasch]

[Sh/c]ould we backport this to v1.9 (and even earlier)? It doesn't interfere with anything else AFAICT, because those products are passed through to _generic_matmatmul! down the pipeline. So this can be viewed a performance fix. @fredrikekre @KristofferC

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.7 (2022-07-19)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using LinearAlgebra, SparseArrays

julia> methods(*, (LinearAlgebra.AbstractQ, SparseMatrixCSC))
# 1 method for generic function "*":
[1] *(A::AbstractMatrix{T} where T, B::AbstractMatrix{T} where T) in LinearAlgebra at /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:151

              _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.8.3 (2022-11-14)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using LinearAlgebra, SparseArrays

julia> methods(*, (LinearAlgebra.AbstractQ, SparseMatrixCSC))
# 1 method for generic function "*":
[1] *(A::AbstractMatrix, B::AbstractMatrix) in LinearAlgebra at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/matmul.jl:139

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.9.0-alpha1 (2022-11-15)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using LinearAlgebra, SparseArrays

julia> methods(*, (LinearAlgebra.AbstractQ, SparseMatrixCSC))
# 2 methods for generic function "*" from Base:
 [1] *(Q::SparseArrays.SPQR.QRSparseQ, B::SparseMatrixCSC)
     @ SparseArrays.SPQR /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/solvers/spqr.jl:287
 [2] *(A::AbstractMatrix, B::AbstractMatrix)
     @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/matmul.jl:139

where [1] points to

(*)(Q::QRSparseQ, B::SparseMatrixCSC) = sparse(Q) * B
(*)(A::SparseMatrixCSC, Q::QRSparseQ) = A * sparse(Q)

[dkarrasch]

We can think about backporting this later, I guess. I'm merging this and include it in Julia master.

[dkarrasch]

@SobhanMP @Wimmerer Before I forget, let me bring this to your attention. Looking at

(*)(Q::QRSparseQ, B::SparseMatrixCSC) = sparse(Q) * B
(*)(A::SparseMatrixCSC, Q::QRSparseQ) = A * sparse(Q)

I was wondering whether for typical sparse Q's, the result of Q*B is sparse or dense? If it's typically dense, perhaps

(*)(Q::QRSparseQ, B::SparseMatrixCSC) = lmul!(Q, Matrix(B))
(*)(A::SparseMatrixCSC, Q::QRSparseQ) = lmul!(Matrix(A), Q)

could be a more performant option (including eltype promotion). I came to think about it because we only have [l/r]mul! methods for dense strided factors (probably for a good reason), so I thought that the general expectation may be that the result is rather dense anyway. With this proposal, we would only allocate the target array, and not first Q and then the result array.