[BUG] `*(block sparse, dense)` (and vice versa) is slow

[xuanzhaogao]

I found that even with an empty BlockSparseArray, doing a matrix vector mul is slow and causes large allocations. Is this expected behavior?

julia> using BlockArrays, BlockSparseArrays

julia> using BenchmarkTools

julia> cnt = fill(16, 128);

julia> bsa = BlockSparseArray{T}(undef, cnt, cnt);

julia> o = ones(2048);

julia> @benchmark $(bsa) * $(o)
BenchmarkTools.Trial: 4 samples with 1 evaluation per sample.
 Range (min … max):  1.603 s …   1.627 s  ┊ GC (min … max): 21.85% … 21.13%
 Time  (median):     1.613 s              ┊ GC (median):    21.64%
 Time  (mean ± σ):   1.614 s ± 10.349 ms  ┊ GC (mean ± σ):  21.63% ±  0.36%

  █          █                      █                     █
  █▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.6 s          Histogram: frequency by time        1.63 s <

 Memory estimate: 17.06 GiB, allocs estimate: 29376515.

[mtfishman]

Thanks for the report. We haven't implemented dense * block sparse yet (we've only implemented block sparse * block sparse, and specifically the case when the blocking on the shared axis is the same), probably it is falling back to a slow generic base implementation.

[xuanzhaogao]

Thanks, it seems this is the only way to make it work

julia> bsa = BlockSparseArray{Float64}(undef, cnt, cnt);

julia> o = BlockSparseArray{Float64}(undef, cnt, [1]);

julia> @benchmark $(bsa) * $(o)
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
 Range (min … max):  5.819 μs …  1.266 ms  ┊ GC (min … max): 0.00% … 98.76%
 Time  (median):     6.306 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   6.551 μs ± 12.607 μs  ┊ GC (mean ± σ):  1.91% ±  0.99%

      ▃█▆▇▃▁▃▁ ▄ ▂          ▁ ▁
  ▂▃▆██████████████▇▇▅▅▅▄▅▆▆████▆▆▅▅▆▅▅▄▄▄▃▃▃▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁ ▄
  5.82 μs        Histogram: frequency by time        7.58 μs <

 Memory estimate: 8.73 KiB, allocs estimate: 136.

[mtfishman]

If you wanted to perform the original operation (where one of the arrays is filled entirely with ones) you could construct it as a block sparse array:

julia> using BlockSparseArrays, BenchmarkTools

julia> r = fill(16, 128);

julia> a1 = BlockSparseArray{Float64}(undef, r, r);

julia> a2 = BlockSparseArray{Float64}(undef, r, [1]);

julia> a2 .= 1;

julia> @benchmark $a1 * $a2
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):   9.625 μs …  51.875 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     10.375 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   10.425 μs ± 559.138 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                          ▂ ▆ ▇▃ █ █ ▇ ▅▃ ▁                     
  ▂▂▁▂▁▂▁▂▂▁▂▁▃▁▃▃▁▄▁▄▁▆█▁█▁█▁██▆█▁█▁█▁██▁█▁█▁▆▅▁▅▁▄▁▃▃▁▃▁▂▁▂▂ ▄
  9.62 μs         Histogram: frequency by time           11 μs <

 Memory estimate: 9.77 KiB, allocs estimate: 136.

Maybe that could be optimized more but at least it is a lot faster than the original code you ran. Probably in practice to implement block sparse * dense I would implement it that way, i.e. convert the dense matrix to the equivalent block sparse one.

I see there is a bug with block sparse mat-vecs (i.e. the code above errors if I set a2 = BlockSparseArray{Float64}(undef, r)), I'll have to look into that issue, that should be easy to fix.

[mtfishman]

Out of curiosity, do you have a use case for block sparse * dense operations? Generally we would use one or the other (so for example if we had some block sparse/symmetric MPS and wanted to operate it with a dense MPO to study a case where the symmetry is broken, we would just convert it to dense).

[xuanzhaogao]

Thanks. I just created a PR for block sparse * dense #198, which simply apply the blocks on views of the dense vector, and returns a BlockSparseVector so that it behavior is unchanged.

[xuanzhaogao]

Out of curiosity, do you have a use case for block sparse * dense operations? Generally we would use one or the other (so for example if we had some block sparse/symmetric MPS and wanted to operate it with a dense MPO to study a case where the symmetry is broken, we would just convert it to dense).

Yes, I do have such a use case, but not for quantum many-body simulations. I use it in my project on boundary integral equations, where I need a correction term that can be represented by such a BlockSparseArray, which is applied onto a dense vector representing the density on boundary.