# Speeding up elementwise Vector-SparseMatrixCSC multiplication broadcasting

**URL:** <https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437>\
**Category:** Performance\
**Tags:** sparse\
**Created:** [April 13, 2022, 2:24pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437 "2022-04-13T14:24:04Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 2:24pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/1 "2022-04-13T14:24:04Z")

</div>

I noticed today a bottleneck in my code which was causing significant performance loss.

Given a vector b and a sparse matrix A, let b \odot A be the elementwise multiplication of b on each column of A. An example in Julia could look like this:

```julia
using SparseArrays
n = 3
b = collect(1:n)
A = sparse(collect(1:n), collect(1:n), rand(1:n))
b .* A

```

It seems however, that this is not natively handled very well as seen by cranking n up.

```julia
@time begin
   n = 3
   b = collect(1:n)
   A = sparse(collect(1:n), collect(1:n), rand(1:n))
   b .* A
end
0.000028 seconds (28 allocations: 2.375 KiB)

@time begin
   n = 10000
   b = collect(1:n)
   A = sparse(collect(1:n), collect(1:n), rand(1:n))
   b .* A
end
0.247870 seconds (51 allocations: 1.491 GiB, 25.08% gc time)

@time begin
   n = 100000
   b = collect(1:n)
   A = sparse(collect(1:n), collect(1:n), rand(1:n))
   b .* A
end
ERROR: OutOfMemoryError()

```

However, if I define the following function, which computes the same thing, I can do much better. It’s significantly faster and uses much less memory.

```julia
function fun(b::Vector, A::SparseMatrixCSC)
   B = copy(A)
   B.nzval .= B.nzval .* b[B.rowval]
   return B
end

@time begin
   n = 100000
   b = collect(1:n)
   A = sparse(collect(1:n), collect(1:n), rand(1:n))
   fun(b, A)
end
0.005198 seconds (43 allocations: 11.446 MiB)

@time begin
   n = Int(1e6)
   b = collect(1:n)
   A = sparse(collect(1:n), collect(1:n), rand(1:n))
   fun(b, A)
end
0.068245 seconds (43 allocations: 114.443 MiB, 7.60% gc time)

```

Is there a reason `.*` does not behave internally in this way? Can I overload it myself to make it behave in this way? I looked a bit at `Base.broadcast` and I have made custom broadcasting for my own callable struct types before, in places where I could gain significant speed-ups by computing the elementwise operations differently than a naive elementwise function call. But I couldn’t seem to figure out how overloading broadcasting for `*` is supposed to be expressed.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 13, 2022, 3:51pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/2 "2022-04-13T15:51:19Z")

</div>

It looks like this is one of the more sane use cases for sparse matrix broadcasts, so shouldn’t it be supported by [SparseArrays](https://github.com/JuliaSparse/SparseArrays.jl)?

---

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 4:19pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/3 "2022-04-13T16:19:24Z")

</div>

Thanks, I’m always happy to be affirmed in my sanity! 😅

That aside, yes, SparseArrays is of course where I’d put it too idealy.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 13, 2022, 4:29pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/4 "2022-04-13T16:29:52Z")

</div>

> [@bcsj](#):
>
> Thanks, I’m always happy to be affirmed in my sanity!

Sorry, I didn’t mean you: I saw some issues on github discussing broadcast between sparse and dense matrices. So it is surprising to me that the basic operation you checked isn’t working appropriately. Ah, wait a second: the proposed optimization doesn’t work for addition? Then you probably should really special case it in your code?

---

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 4:36pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/5 "2022-04-13T16:36:42Z")

</div>

Haha, all good.

Yes, you are right, this optimization is exclusively useful for multiplication. Should probably be in the post title, I’ll edit it in.

Broadcasting a vector elementwise addition to each column would result in a mostly dense matrix, unless the vector itself was very sparse.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 13, 2022, 5:04pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/6 "2022-04-13T17:04:38Z")

</div>

> [@bcsj](#):
>
> elementwise multiplication of b on each column of A

You could also express this as `Diagonal(b) * A` (using the `Diagonal` type from `LinearAlgebra`), since SparseArrays has a specialized method for this product.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [April 13, 2022, 5:33pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/7 "2022-04-13T17:33:15Z")

</div>

The difference lies in `b .* A` vs. `B.nzval .= B.nzval .* b[B.rowval]` - the former manifests the full, dense matrix (with all zero elements as well), while the latter only saves what’s required - the nonzeros:

```julia
julia> b .* A |> typeof
Matrix{Int64} (alias for Array{Int64, 2})

julia> fun(b, A) |> typeof
SparseMatrixCSC{Int64, Int64}

```

You can probably avoid that full allocation if you pass in a sparse array as a receiving buffer and do `B .= b .* A`, since that may specialize on the `setindex!` during broadcasting.

By the way, there’s no need to `collect` your ranges - they are already `<: AbstractVector` and behave like their “full” counterparts, while taking up a miniscule fraction of their memory.

---

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 5:40pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/8 "2022-04-13T17:40:07Z")

</div>

I’m not sure that is correct. On my machine the output type remains sparse.

```julia
julia> using SparseArrays
julia> n = 3
julia> b = collect(1:3)
3-element Array{Int64,1}:
 1
 2
 3
julia> a = sparse(collect(1:n),collect(1:n),rand(1:n))
3×3 SparseMatrixCSC{Int64,Int64} with 3 stored entries:
  [1, 1] = 2
  [2, 2] = 2
  [3, 3] = 2
julia> b .* a
3×3 SparseMatrixCSC{Int64,Int64} with 3 stored entries:
  [1, 1] = 2
  [2, 2] = 4
  [3, 3] = 6
julia> b .* a |> typeof
SparseMatrixCSC{Int64,Int64}

```

---

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 5:50pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/9 "2022-04-13T17:50:38Z")

</div>

That does seem like a good work around solution. Also has marginally better performance compared to my little function.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [April 13, 2022, 6:34pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/10 "2022-04-13T18:34:42Z")

</div>

> [@bcsj](#):
>
> I’m not sure that is correct. On my machine the output type remains sparse.

That’s very interesting! I’m on a recent master, so either this is a regression or intentional 🤔

---

<div class="post-metadata">

**Author:** ![bcsj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcsj/32/18526_2.png) [@bcsj](https://discourse.julialang.org/u/bcsj)\
**Post date:** [April 13, 2022, 7:08pm UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/11 "2022-04-13T19:08:51Z")

</div>

Very curious indeed. The code I ran in response to you comment was on my private laptop which runs an older `1.5.2 (2020-09-23)`. But I’m 99% confident my work laptop, which runs a later version, though not the very latest, also remained sparse in output.

---

<div class="post-metadata">

**Author:** ![axsk](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axsk/32/1952_2.png) [@axsk](https://discourse.julialang.org/u/axsk)\
**Post date:** [December 12, 2024, 9:45am UTC](https://discourse.julialang.org/t/speeding-up-elementwise-vector-sparsematrixcsc-multiplication-broadcasting/79437/12 "2024-12-12T09:45:23Z")

</div>

This still is not working out of the box, right?