Operation on transpose sparse array is slow

General Usage Performance
2

Transposing a sparse matrix with Y = X' doesn’t actually copy any data, it wraps it with a special Adjoint type that gives a (conjugate-)transposed view of the underlying X matrix.

Adjoint views of sparse matrices have several optimized methods, e.g. Y * vector is fast, but not all of the optimized methods for sparse matrices are implemented for the Adjoint wrapper. In particular, “logical indexing” with a vector of booleans has a fast implementation for sparse matrices, but not for Adjoint wrappers thereof, and in consequence the corresponding operation on Y falls back to the generic method for dense matrices, which is very slow because it does not exploit sparsity.

Three possible workarounds:

  1. Set Y = SparseMatrixCSC(X') — this makes an explicit sparse-matrix copy of the Adjoint wrapper, which uses up some extra memory but then subsequent operations will be fast.
  2. Re-organize your code so that you do your calculations directly on the un-transposed matrix X
  3. Define the missing specialized methods for Adjoint wrappers of sparse matrices. (This is “type piracy” but may be a reasonable stopgap measure in the short run.)

For example, your code above should become fast again if we define:

using LinearAlgebra, SparseArrays
Base.getindex(A::Adjoint{<:Any,<:SparseArrays.AbstractSparseMatrixCSC}, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = parent(A)[findall(J),findall(I)]'

If someone wants to volunteer to make a PR to SparseArrays.jl with this and similar enhancements, that would be great. See also getindex is very slow for adjoints of sparse arrays · Issue #36 · JuliaSparse/SparseArrays.jl · GitHub and Strange performance with Adjoint structures

3 Likes