# Sparse matrix vector faster than dense matrix vector?

**URL:** <https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446>\
**Category:** Performance\
**Created:** [September 1, 2023, 2:25pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446 "2023-09-01T14:25:01Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![dgleich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dgleich/32/1494_2.png) [@dgleich](https://discourse.julialang.org/u/dgleich)\
**Post date:** [September 1, 2023, 2:25pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/1 "2023-09-01T14:25:01Z")

</div>

As I was teaching sparse matrices in class, I was humbled by a theory-practice gap.

The idea in theory: don’t use sparse matrices when you really have dense matrix.

The idea in practice: here, let me show you in Julia!

* * *

```julia
n = 10 
using BenchmarkTools, SparseArrays
@benchmark A*v setup=begin
    A = randn(n,n)
    v = randn(n)
end
using BenchmarkTools
@benchmark A*v setup=begin
    A = sparse(randn(n,n))
    v = randn(n)
end

```

And then this shows:

**Dense case**

```julia
BenchmarkTools.Trial: 10000 samples with 961 evaluations.
 Range (min … max): 86.672 ns … 798.907 ns ┊ GC (min … max): 0.00% … 85.88%
 Time (median): 89.663 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 92.006 ns ± 23.905 ns ┊ GC (mean ± σ): 1.33% ± 4.33%

```

**Sparse case**

```julia
BenchmarkTools.Trial: 10000 samples with 965 evaluations.
 Range (min … max): 81.649 ns … 804.188 ns ┊ GC (min … max): 0.00% … 86.96%
 Time (median): 86.010 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 89.411 ns ± 31.736 ns ┊ GC (mean ± σ): 1.85% ± 4.61%

```

So the sparse case is 3 ns faster than the dense case ?!?

As I said, I was humbled in lecture 🙂

* * *

Notes

- Yes, the difference goes away for larger n and sparse is about 3x slower than dense.
- I haven’t seen this on any Intel/AMD x86-64 architectures.
- This does _not_ happen on Julia 1.8, so something changed in 1.9 that made the sparse case faster.
- Using the following code shows that the dense case really should be faster…

```julia
## Copy-pasted sparse matrix vector product without allocations
using BenchmarkTools, SparseArrays, LinearAlgebra
function _spmatmul!(C, A, B, α, β)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    β != one(β) && LinearAlgebra._rmul_or_fill!(C, β)
    for k in 1:size(C, 2)
        @inbounds for col in 1:size(A, 2)
            αxj = B[col,k] * α
            for j in nzrange(A, col)
                C[rv[j], k] += nzv[j]*αxj
            end
        end
    end
    C
end
@benchmark _spmatmul!(x,A,v,1.0,0.0) setup=begin
    A = sparse(randn(n,n))
    v = randn(n)
    x = zeros(n) 
end
## Equivalent dense-matrix vector product 
using BenchmarkTools, SparseArrays, LinearAlgebra
function _matmul!(C, A, B, α, β)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    β != one(β) && LinearAlgebra._rmul_or_fill!(C, β)
    for k in 1:size(C, 2)
        @inbounds for col in 1:size(A, 2)
            αxj = B[col,k] * α
            for i in 1:size(A,1)
                C[i, k] += A[i,col]*αxj
            end
        end
    end
    C
end
@benchmark _matmul!(x,A,v,1.0,0.0) setup=begin
    A = randn(n,n)
    v = randn(n)
    x = zeros(n) 
end
@benchmark mul!(x,A,v,1.0,0.0) setup=begin
    A = randn(n,n)
    v = randn(n)
    x = zeros(n) 
end

```

Let’s do our bakeoff…

Julia 1.8.4 (MacBook Air M2)

dense \_matmul! - 44.7 ns  
dense LinearAlgebra.mul! - 73.8 ns  
sparse \_spmatmul! - 74.3 ns

Julia 1.9.2 (MacBook Air M2)

dense \_matmul! - 43.8 ns  
dense LinearAlgebra.mul! - 73.0 ns  
sparse \_spmatmul! - 68.2 ns

* * *

The question: does anyone know what might have changed in 1.9 that would have made the sparse method faster?

Also, anyone know why the dense code for small n is slower than expected?

---

<div class="post-metadata">

**Author:** ![devel-chm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devel-chm/32/3572_2.png) [@devel-chm](https://discourse.julialang.org/u/devel-chm)\
**Post date:** [September 1, 2023, 3:02pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/2 "2023-09-01T15:02:38Z")

</div>

Some observations:

- It doesn’t look like your sparse array is actually sparse.  
The fill factor looks greater than 0.95 (at least)

- Standard deviation in the measurements is \>24 nsec  
so a 3-5 nsec difference is in the noise (not significant)

- With matrix multiply routines, hardware/processor differences  
can affect the results.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 1, 2023, 3:57pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/3 "2023-09-01T15:57:52Z")

</div>

I believe what’s happening here is that the M1/M2 chips is really good at prefetching memory though pointers which makes the sparse stuff work better. Also with a completely dense sparse matrix, the processor is probably able to fully branch predict all the branches on where the next data will be. The really bad case for sparse arrays is between ~10% and 90% full since that is where it’s very unpredictable.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 1, 2023, 4:07pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/4 "2023-09-01T16:07:58Z")

</div>

@dgleich I wonder what the situation with multithreading is? The dense matrix multiplication would run on multiple threads by default, I believe? Now, that is probably moot for the small matrices, but at some point it will mix things up, won’t it?

---

<div class="post-metadata">

**Author:** ![dgleich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dgleich/32/1494_2.png) [@dgleich](https://discourse.julialang.org/u/dgleich)\
**Post date:** [September 1, 2023, 4:08pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/5 "2023-09-01T16:08:13Z")

</div>

The prefetching on the M1/M2 doesn’t explain why this wasn’t as fast on Julia 1.8 though. Also doesn’t explain why the dense case is slower than you might expect.

---

<div class="post-metadata">

**Author:** ![dgleich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dgleich/32/1494_2.png) [@dgleich](https://discourse.julialang.org/u/dgleich)\
**Post date:** [September 1, 2023, 4:14pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/6 "2023-09-01T16:14:51Z")

</div>

@PetrKryslUCSD Oh – good idea to check that!

```julia
Threads.nthreads() = 1
LinearAlgebra.BLAS.get_num_threads() = 4

```

so I re-ran after `LinearAlgebra.BLAS.set_num_threads(1)`

and got the same 73 ns result. So that doesn’t help explain it. I will have to double-check with bigger matrices if that explains the factor of 3 vs. sparse.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 1, 2023, 4:44pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/7 "2023-09-01T16:44:41Z")

</div>

The dense case is slower because there is a pretty high `O(1)` overhead for calling BLAS (due to calling conventions and the like). I’m also not that surprised about Julia 1.8 since we had pretty poor support for M1 in 1.8.

---

<div class="post-metadata">

**Author:** ![dgleich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dgleich/32/1494_2.png) [@dgleich](https://discourse.julialang.org/u/dgleich)\
**Post date:** [September 1, 2023, 6:08pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/8 "2023-09-01T18:08:14Z")

</div>

Okay thanks – I think that’s the other missing piece (the high BLAS overhead) – would it be far to say this is most likely due to the following?

- BLAS overhead is higher on AppleARM than x86-64?
- AppleARM gets better perf for simple codes than x86-64?

So these combine to give a different trade off point between simple codes/etc for small n?

e.g. I was finally able to replicate similar behavior on on Intel machine with n=6 matrix.

Thanks!

---

<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:** [September 1, 2023, 6:55pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/9 "2023-09-01T18:55:48Z")

</div>

> [@dgleich](#):
>
> So these combine to give a different trade off point between simple codes/etc for small n?

If you have such small matrices and the size is a fixed compile-time constant and you care a lot about performance, you are often better off telling the compiler the size of the array with [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) so that it just unrolls everything.

---

<div class="post-metadata">

**Author:** ![dgleich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dgleich/32/1494_2.png) [@dgleich](https://discourse.julialang.org/u/dgleich)\
**Post date:** [September 5, 2023, 2:15pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-faster-than-dense-matrix-vector/103446/10 "2023-09-05T14:15:44Z")

</div>

> [@stevengj](#):
>
> If you have such small matrices and the size is a fixed compile-time constant and you care a lot about performance, you are often better off telling the compiler the size of the array with [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) so that it just unrolls everything.

Thanks for pointing this out!

Agree that StaticArrays would be better for a fixed small size. In this case, it runs in 6 ns 🙂