`stack` not specialized for SparseVectors as it is for`hcat`

[VPetukhov]

This is a follow up on #49887 . Replacing calls to hcat with stack leads to 500+ fold drop in performance:

Code:

using SparseArrays
using BenchmarkTools

vs = [sprand(Bool, 100, 0.1) for i in 1:20000];

Benchmarks:

@benchmark stack(vs, dims=2)

 Range (min … max):  698.288 ms … 750.043 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     703.131 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   708.654 ms ±  17.297 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%
 Memory estimate: 3.79 MiB, allocs estimate: 27.

@benchmark hcat(vs...)

BenchmarkTools.Trial: 2524 samples with 1 evaluation.
 Range (min … max):  1.112 ms … 31.250 ms  ┊ GC (min … max): 0.00% … 88.59%
 Time  (median):     1.814 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.968 ms ±  1.940 ms  ┊ GC (mean ± σ):  6.47% ±  6.20%
 Memory estimate: 2.64 MiB, allocs estimate: 20008.

The test was done on julia 1.8.5, as hcat(vs...) fails on 1.9.*. However, testing stack on 1.9.2 produces almost the same results.

versioninfo()

For 1.8.5:

Julia Version 1.8.5
Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × Intel(R) Xeon(R) CPU E5-2697A v4 @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, broadwell)
  Threads: 1 on 64 virtual cores

For 1.9.2:

Julia Version 1.9.2
Commit e4ee485e909 (2023-07-05 09:39 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × Intel(R) Xeon(R) CPU E5-2697A v4 @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, broadwell)
  Threads: 1 on 64 virtual cores

[mcabbott]

As noted, the splat of 20_000 items gives a StackOverflowError on Julia master (and other recent versions). That's issue #49887 closed on the grounds that this is a bad idea.

Maybe worth adding that the specialised hcat(Xin::AbstractSparseVector...) method seems to have become slower, both for a vector of 200 items (small enough to avoid the overflow) and for a tuple of 20. On master this method appears to be slower than stack / reduce(hvat, vs), neither of which have any specialisations for sparse arrays. (Their speeds have not changed.)

This slow-down is an orthogonal question to whether stack should gain methods for SparseArrays. Perhaps they should be different issues.

julia> vs_tup = Tuple([sprand(Float32, 100, 0.1) for i in 1:20]);  # a reasonable thing to splat

julia> m1t = @btime stack($vs_tup);
  min 39.959 μs, mean 49.163 μs (11 allocations, 11.67 KiB)  # similar on 1.8 (via Compat) & master

julia> m2t = @btime reduce(hcat, $vs_tup);  # pairwise
  min 78.708 μs, mean 113.434 μs (1080 allocations, 200.83 KiB)

julia> m3t = @btime hcat($vs_tup...);
  min 105.625 μs, mean 112.608 μs (227 allocations, 17.98 KiB) # Julia 1.11.0-DEV.137
  min 643.293 ns, mean 1.047 μs (24 allocations, 3.86 KiB)     # Julia 1.8.2, note the units!

julia> using SparseArrays, BenchmarkTools, Compat

julia> vs_float = [sprand(Float32, 100, 0.1) for i in 1:200];  # fewer vectors, avoids stack overflow

julia> m1 = @btime stack($vs_float);
  min 445.459 μs, mean 472.581 μs (50 allocations, 85.13 KiB)

julia> m2 = @btime reduce(hcat, $vs_float);  # fast path AbstractVector{<:AbstractVecOrMat}
  min 705.417 μs, mean 1.201 ms (2322 allocations, 2.90 MiB)

julia> m3 = @btime hcat($vs_float...);  # method for AbstractSparseVector...
  min 6.488 ms, mean 6.678 ms (2225 allocations, 763.93 KiB) # Julia 1.11.0-DEV.137
  min 8.819 μs, mean 13.111 μs (204 allocations, 31.38 KiB)  # Julia 1.8.2

julia> @which hcat(vs_float...)
hcat(Xin::AbstractSparseVector...)
     @ SparseArrays ~/.julia/dev/julia/usr/share/julia/stdlib/v1.11/SparseArrays/src/sparsevector.jl:1107

julia> @which reduce(hcat, vs_float)  # no specialisation for sparse
reduce(::typeof(hcat), A::AbstractVector{<:AbstractVecOrMat})
     @ Base abstractarray.jl:1708

julia> @which reduce(hcat, vs_tup)  # worse, pairwise
reduce(op, itr; kw...)
     @ Base reduce.jl:490