# Non-linear latency of sparse-dense matrix multiplication

**URL:** <https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424>\
**Category:** Numerics\
**Created:** [March 3, 2021, 5:08pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424 "2021-03-03T17:08:11Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:08pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/1 "2021-03-03T17:08:11Z")

</div>

Hello,

I’m measuring the latency of `mul!(C, A, B)`, where `C` and `B` are dense and `A` is sparse, and I’m seeing something surprising: the latency isn’t linear in the number of operations required. Does anyone know why this is?

The number of columns of `A` is `m=1000`, the number of columns of `B` is `k=3`, and I’m varying the number of rows of `A`, denoted by `n`. I’m on an Intel system and I’m including `MKLSparse` (but I see the same behaviour regardless of if I use `MKLSparse`).

My experiment looks like the following code, except that it’s in a function, I’m running it many times, and I discard the first sample (since it includes compilation time)

```julia
A = sprand(n, m, 0.05)
B = randn(m, k)
C = zeros(n, k)
latency_mul = @elapsed mul!(C, A, B)
latency_mymul! = @elapsed mymul!(C, A, B) # mymul! defined below

```

Here’s the definition of `mymul!`:

```julia
function mymul!(C::Matrix, A::SparseMatrixCSC, B::Matrix)
    @boundscheck size(C, 1) == size(A, 1) || throw(DimensionMismatch("C has dimensions $(size(C)), but A has dimensions $(size(A))"))    
    @boundscheck size(C, 2) == size(B, 2) || throw(DimensionMismatch("C has dimensions $(size(C)), but B has dimensions $(size(B))"))        
    @boundscheck size(A, 2) == size(B, 1) || throw(DimensionMismatch("A has dimensions $(size(A)), but B has dimensions $(size(B))"))
    m, n = size(A)
    k = size(B, 2)
    rows = rowvals(A)
    vals = nonzeros(A)
    C .= 0
    for col = 1:n # column of A
        for i in nzrange(A, col)
            @inbounds row = rows[i] # row of A
            @inbounds val = vals[i] # A[row, col]
            @simd for j = 1:k # columns indices of B
                @inbounds C[row, j] += val * B[col, j]
            end
        end
    end
    C
end

```

This is the latency I’m seeing. I’d have expected that latency would increase linearly with the number of rows, which is what I’m seeing for `mymul!` (ish, for large enough matrices). However, for `mul!` I see something else. The dashed lines are quadratic functions fit to the data. The latency looks linear in the number of rows when running the same experiment for dense matrices so it’s got something to do with the sparse matrix multiplication.

![sparse_mul](https://global.discourse-cdn.com/julialang/original/3X/7/2/7247bcf9756b1fd0d31b33293c6a918d1fcc7c83.png)

---

<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:** [March 3, 2021, 5:14pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/2 "2021-03-03T17:14:39Z")

</div>

What a coincidence! I just need an ad hoc implementation of subspace iteration, where one of the major costs is the multiplication of a beefy sparse matrix with a huge dense matrix of vectors as columns.  
I guess my system is not as fast as yours, since in my case the execution time (which I think you call “latency”) is on the order of seconds. (Of course, it also depends on how sparse the matrix is. I know.)

I haven’t actually tried what you describe as an in-house solution. I rather tried using threads to calculate the individual columns. I did not get any speed up, unfortunately.

---

<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:** [March 3, 2021, 5:18pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/3 "2021-03-03T17:18:54Z")

</div>

Isn’t the access to `C` row-major in your code?

---

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:20pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/4 "2021-03-03T17:20:13Z")

</div>

Threading should ideally be handled by the matrix multiplication routines without you having to think about it. If you’re on an Intel system I’d recommend trying out `MKLSparse`. There’s some text suggesting it can speed up sparse matrix operations significantly when using multiple threads, see [https://github.com/JuliaSparse/MKLSparse.jl/issues/24#issuecomment-658086827](https://github.com/JuliaSparse/MKLSparse.jl/issues/24#issuecomment-658086827)

---

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:22pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/5 "2021-03-03T17:22:09Z")

</div>

The access to `C` is based on what I believe is the fastest way to iterate over the non-zero entries of `A`. I may be wrong though. On the other hand, I don’t really care about the `mymul!`. I’ve only implemented it as a reference.

---

<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:** [March 3, 2021, 5:23pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/6 "2021-03-03T17:23:35Z")

</div>

I understand. Still, one wonders if it could be done faster.

---

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:24pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/7 "2021-03-03T17:24:51Z")

</div>

I’m more interested in the non-linearity. There’s something I’m not understanding about how `mul!` is computed.

---

<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:** [March 3, 2021, 5:25pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/8 "2021-03-03T17:25:24Z")

</div>

Is the sparsity changing as you increase the size?

---

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:30pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/9 "2021-03-03T17:30:06Z")

</div>

No.  
Everything is kept constant except the number of rows.

---

<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:** [March 3, 2021, 5:33pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/10 "2021-03-03T17:33:22Z")

</div>

Does that mean the “sparsity”, or the average number of nonzeros per column?

---

<div class="post-metadata">

**Author:** ![albin](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@albin](https://discourse.julialang.org/u/albin)\
**Post date:** [March 3, 2021, 5:35pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/11 "2021-03-03T17:35:45Z")

</div>

Does it make a difference? 🙂

I generate `A` like this:

```julia
A = sprand(n, 1000, 0.05)

```

---

<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:** [March 3, 2021, 5:37pm UTC](https://discourse.julialang.org/t/non-linear-latency-of-sparse-dense-matrix-multiplication/56424/12 "2021-03-03T17:37:52Z")

</div>

I suppose not that much. However the clustering of the non-zeros might.