Faster blockdiag for uniform input value-index types

[irhum]

The current blockdiag implementation in SparseArrays uses type promotion to identify what types the values and index should be for the output SparseArray.

This can be quite slow when the input to blockdiag is many small, sparse arrays (for the sake of example, over 2500+). If all the input sparse arrays have the same value and index types, then it's not necessary to compute them via type promotion and we can directly proceed to creating the new matrix.

This can provide significant speedups; for M generated with

M = [sprand(rand(15:50), rand(15:50), 0.02) for i in 1:2500]

and benchmarking the existing implementation with BenchmarkTools

@benchmark output = blockdiag(M...)

we get

BenchmarkTools.Trial: 
  memory estimate:  192.35 MiB
  allocs estimate:  37956
  --------------
  minimum time:     1.404 s (1.30% GC)
  median time:      1.434 s (1.31% GC)
  mean time:        1.438 s (1.24% GC)
  maximum time:     1.482 s (0.89% GC)
  --------------
  samples:          4
  evals/sample:     1

running the same benchmark, with same M, with the new implementation, we get a speedup of nearly 1000x

@benchmark output = blockdiag(M...)

  memory estimate:  3.00 MiB
  allocs estimate:  5019
  --------------
  minimum time:     1.072 ms (0.00% GC)
  median time:      1.162 ms (0.00% GC)
  mean time:        1.386 ms (11.86% GC)
  maximum time:     5.059 ms (60.41% GC)
  --------------
  samples:          3591
  evals/sample:     1

the current blockdiag implementation is hence split into 3 possible methods

blockdiag(X::AbstractSparseMatrixCSC...) # for the general case
blockdiag(X::AbstractSparseMatrixCSC{Tv, Ti}...) where {Tv, Ti <: Integer} # for increased speed when all input X have the same index and value types
blockdiag(::Type{Tv}, ::Type{Ti}, X::AbstractSparseMatrixCSC...) where {Tv, Ti <: Integer} # internally called by the two above methods, and available to the user when they're dealing with heterogenous inputs, and wish to specifiy the output index and value types manually for the increased speed benefits

[irhum]

added another method to deal with the empty input case, blockdiag(), otherwise it defaults to

blockdiag(X::AbstractSparseMatrixCSC{Tv, Ti}...) where {Tv, Ti <: Integer}

instead of

blockdiag(X::AbstractSparseMatrixCSC...)

resulting in a Tv not defined error, making an explicit definition necessary.

The pull request, I would suggest is important given the case where all inputs have the same types for the value and index would be more common (and would benefit from increased performance) than the heterogenous case

[dkarrasch]

Very nice contribution, @irhum, and welcome! I have one comment that I'm curious what other people think about.

[dkarrasch]

You'll need to rename the function at the call sites as well.

[irhum]

The function call sites have been changed to make the appropriate _blockdiag call.

Since we also need to account for the case where no inputs are given, a clean solution seems to be to pass it as an optional argument to one of the methods, specifically:

blockdiag(X::AbstractSparseMatrixCSC{Tv, Ti}...=spzeros(0,0)) where {Tv, Ti <: Integer}

Since we're explicitly handing it, isempty(X) in

Ti = isempty(X) ? Int : promote_type(map(x->eltype(rowvals(x)), X)...)

is always false, and hence can be replaced with

Ti = promote_type(map(x->eltype(rowvals(x)), X)...)

Opinions?

[dkarrasch]

Looks good in principle, but there is this subtle issue with empty arguments.

[dkarrasch]

This default value makes the function change behavior for empty arguments. Currently, we have

julia> A = blockdiag()
0×0 SparseMatrixCSC{Union{},Int64} with 0 stored entries

julia> B = spzeros(0,0)
0×0 SparseMatrixCSC{Float64,Int64} with 0 stored entries

Changing that behavior may be for the better, but I'm not sure. This definitely requires a broader discussion, I'd say.

[irhum]

I would rather not change the behavior inside this pull request anyway; the goal here is to provide increased speed when dealing with CSC arrays with the same types for the index and values; if needed, a separate issue can be opened whether the default should be changed.

Created a new commit that fixes this (the output is now SparseMatrixCSC{Union{},Int64}) by explicitly spelling out the blockdiag() case, to make it clear what it does today, to avoid any confusions in the future.

[dkarrasch]

I'd suggest to have

Suggested change

	blockdiag() = spzeros(Union{}, Int, 0, 0)
	blockdiag() = spzeros(promote_type(), Int, 0, 0)

to avoid hardcoding some "non-sense" type and let the ecosystem decide. Pinging @tkf for this nasty "reduction over empty tuples" issue.

[irhum]

Thanks for pointing this out, it is more sensible this way. I believe that resolves the issue for now, of speeding up this function while still making sure blockdiag() output remains unchanged. Would appreciate what @tkf has to say about the reduction on empty tuple issue

[tkf]

I think Union{} makes sense (i.e., it's the identity element of promote_type as a monoid). Though promote_type() sounds also fine. It's just a LISPy way of saying that we want the identity of promote_type.

[dkarrasch]

Thanks, @tkf. I just had no idea what might be subject to potential future design changes and what is "set in stone". I'll go ahead and merge this as is then.