# Sparse matrix-vector product: much more slow than Matlab

**URL:** https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840
**Category:** Performance
**Tags:** matlab, optimization
**Created:** [December 18, 2017, 9:14pm UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840 "2017-12-18T21:14:14Z")
**Posts on this page:** 5
**Page:** 2

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 20, 2017, 2:19am UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840/21 "2017-12-20T02:19:07Z")

</div>

Are you sure it’s still computing the same thing? I would do the changes one at a time making sure the result is the same each time. If you change it to non-allocating iteratively, it’ll be clear what change caused the problem. I cannot see what it is right off the bat though. What’s `sparseM,sparseR,sparseS,q`? Can you make it an MWE, i.e. should be able to copy/paste to run it.

---

<div class="post-metadata">

### Author: ![hongyan](https://avatars.discourse-cdn.com/v4/letter/h/4da419/32.png) [@hongyan](https://discourse.julialang.org/u/hongyan)
#### Post date: [December 20, 2017, 4:07am UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840/22 "2017-12-20T04:07:18Z")

</div>

Hi @ChrisRackauckas, by your suggestions, I do the changes one at a time and find what caused the difference. Here are the codes:

Version 1: (You can copy, paste and run)

```julia
numGrids = 799;
q = rand(numGrids);
rows1 = zeros(Int64,numGrids*numGrids);
cols1 = zeros(Int64,numGrids*numGrids);
vals1 = zeros(numGrids*numGrids);
rows2 = zeros(Int64,numGrids*numGrids);
cols2 = zeros(Int64,numGrids*numGrids);
vals2 = zeros(numGrids*numGrids);
rows3 = zeros(Int64,numGrids*numGrids*numGrids);
cols3 = zeros(Int64,numGrids*numGrids*numGrids);
vals3 = zeros(numGrids*numGrids*numGrids);

count1 = 0;
count2 = 0;
count3 = 0;

for i = 1:numGrids
    for j = 1:numGrids
        r = rand();
        if (r > 0.965)
            count1 += 1;
            rows1[count1] = i;
            cols1[count1] = j;
            vals1[count1] = r;
        end

        r = rand();
        if (r > 0.965)
            count2 += 1;
            rows2[count2] = i;
            cols2[count2] = j;
            vals2[count2] = r;
        end
    end

    for k = 1:(numGrids*numGrids)
        r = rand();
        if (r > 0.977)
            count3 += 1;
            rows3[count3] = i;
            cols3[count3] = k;
            vals3[count3] = r;
        end
    end
end

rows1 = rows1[1:count1];
cols1 = cols1[1:count1];
vals1 = vals1[1:count1];
rows2 = rows2[1:count2];
cols2 = cols2[1:count2];
vals2 = vals2[1:count2];
rows3 = rows3[1:count3];
cols3 = cols3[1:count3];
vals3 = vals3[1:count3];

sparseM = sparse(rows1,cols1,vals1,numGrids,numGrids);
sparseR = sparse(rows2,cols2,vals2,numGrids,numGrids);
sparseS = sparse(rows3,cols3,vals3,numGrids,numGrids*numGrids);

function testFun(sparseM,sparseR,sparseS,q,numIters,numGrids)
    tmp = Array{Float64}(numGrids, numGrids);
    tmp2 = Array{Float64}(numGrids, numGrids);
    tmp3 = Array{Float64}(numGrids, numGrids);
    for time = 1:numIters
        tmp .= 0.5 .* 0.05 .* (sparseR .+ reshape(q'*sparseS, numGrids, numGrids));
        tmp2 .= sparseM .+ tmp;
        tmp3 .= sparseM .- tmp;
        rhs = tmp2 * q;
        forwardQ = tmp3 \ rhs;
    end
end

@time testFun(sparseM,sparseR,sparseS,q,10,799)

```

It shows **1.535764 seconds (113.92 k allocations: 122.168 MiB, 10.65% gc time)**

Version 2: (you can also copy, paste and run)

```julia
numGrids = 799;
q = rand(numGrids);
rows1 = zeros(Int64,numGrids*numGrids);
cols1 = zeros(Int64,numGrids*numGrids);
vals1 = zeros(numGrids*numGrids);
rows2 = zeros(Int64,numGrids*numGrids);
cols2 = zeros(Int64,numGrids*numGrids);
vals2 = zeros(numGrids*numGrids);
rows3 = zeros(Int64,numGrids*numGrids*numGrids);
cols3 = zeros(Int64,numGrids*numGrids*numGrids);
vals3 = zeros(numGrids*numGrids*numGrids);

count1 = 0;
count2 = 0;
count3 = 0;

for i = 1:numGrids
    for j = 1:numGrids
        r = rand();
        if (r > 0.965)
            count1 += 1;
            rows1[count1] = i;
            cols1[count1] = j;
            vals1[count1] = r;
        end

        r = rand();
        if (r > 0.965)
            count2 += 1;
            rows2[count2] = i;
            cols2[count2] = j;
            vals2[count2] = r;
        end
    end

    for k = 1:(numGrids*numGrids)
        r = rand();
        if (r > 0.977)
            count3 += 1;
            rows3[count3] = i;
            cols3[count3] = k;
            vals3[count3] = r;
        end
    end
end

rows1 = rows1[1:count1];
cols1 = cols1[1:count1];
vals1 = vals1[1:count1];
rows2 = rows2[1:count2];
cols2 = cols2[1:count2];
vals2 = vals2[1:count2];
rows3 = rows3[1:count3];
cols3 = cols3[1:count3];
vals3 = vals3[1:count3];

sparseM = sparse(rows1,cols1,vals1,numGrids,numGrids);
sparseR = sparse(rows2,cols2,vals2,numGrids,numGrids);
sparseS = sparse(rows3,cols3,vals3,numGrids,numGrids*numGrids);

function testFun(sparseM,sparseR,sparseS,q,numIters,numGrids)
    tmp = Array{Float64}(numGrids, numGrids);
    tmp2 = Array{Float64}(numGrids, numGrids);
    tmp3 = Array{Float64}(numGrids, numGrids);
    longVector = Array{Float64}(1,numGrids * numGrids);
    for time = 1:numIters
        At_mul_B!(longVector, q, sparseS);
        tmp .= 0.5 .* 0.05 .* (sparseR .+ reshape(longVector, numGrids, numGrids));
        tmp2 .= sparseM .+ tmp;
        tmp3 .= sparseM .- tmp;
        rhs = tmp2 * q;
        forwardQ = tmp3 \ rhs;
    end
end

@time testFun(sparseM,sparseR,sparseS,q,10,799)

```

It shows **18.268159 seconds (55.34 k allocations: 31.397 MiB)**

Version 2 is much more slow than version 1 but it has less memory allocations. The only difference between version 1 and version 2 is that in version 2, I used **At\_mul\_B!(longVector, q, sparseS);** do you know why it makes the speed so slow?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 20, 2017, 4:41am UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840/23 "2017-12-20T04:41:47Z")

</div>

I ran

```julia
function testFun2(sparseM,sparseR,sparseS,q,numIters,numGrids)
    tmp = Array{Float64}(numGrids, numGrids);
    tmp2 = Array{Float64}(numGrids, numGrids);
    tmp3 = Array{Float64}(numGrids, numGrids);
    longVector = Array{Float64}(1,numGrids * numGrids);
    for time = 1:numIters
        At_mul_B!(longVector, q, sparseS);
        @show @which q'*sparseS
        tmp .= 0.5 .* 0.05 .* (sparseR .+ reshape(q'*sparseS, numGrids, numGrids));
        tmp2 .= sparseM .+ tmp;
        tmp3 .= sparseM .- tmp;
        rhs = tmp2 * q;
        forwardQ = tmp3 \ rhs;
    end
end

@time testFun2(sparseM,sparseR,sparseS,q,10,799)

```

which tells me two things. First I get the speed slowdown, so `At_mul_B!` is seemingly not what I want (it must just be slow for sparse matrices? That makes sense given how sparse matrices are stored). Then it shows

```julia
@which(q' * sparseS) = *(rowvec::RowVector, mat::AbstractArray{T,2} where T) in Base.LinAlg at linalg\rowvector.jl:181

```

which points me to

```julia
@inline *(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat.' * transpose(rowvec))

```

and checking

```julia
@which rand(3,3).' * transpose(rand(3)')

```

points me to

```julia
function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S}
    TS = promote_op(matprod, T, S)
    A_mul_B!(similar(x, TS, size(A,1)), A, convert(AbstractVector{TS}, x))
end

```

so that tells us how the fast version is lowering (and it shows you where the temporary is made). But I think I missed something. This brings me to two separate versions:

```julia
function testFun2(sparseM,sparseR,sparseS_transp,q,numIters,numGrids)
    forwardQ = Array{Float64}(numGrids,1);
    tmp = Array{Float64}(numGrids, numGrids);
    tmp2 = Array{Float64}(numGrids, numGrids);
    tmp3 = Array{Float64}(numGrids, numGrids);
    longVector = Array{Float64}(numGrids * numGrids);
    for time = 1:numIters
        A_mul_B!(longVector, sparseS_transp, q);
        tmp .= 0.5 .* 0.05 .* (sparseR .+ reshape(longVector, numGrids, numGrids));
        tmp2 .= sparseM .+ tmp;
        tmp3 .= sparseM .- tmp;
        A_mul_B!(forwardQ,tmp2,q);
        A_ldiv_B!(lufact!(tmp3),forwardQ);
    end
end

@time testFun2(sparseM,sparseR,sparseS.',q,10,799)

function testFun(sparseM,sparseR,sparseS,q,numIters,numGrids)
    forwardQ = Array{Float64}(numGrids,1);
    tmp = Array{Float64}(numGrids, numGrids);
    tmp2 = Array{Float64}(numGrids, numGrids);
    tmp3 = Array{Float64}(numGrids, numGrids);
    for time = 1:numIters
        tmp .= 0.5 .* 0.05 .* (sparseR .+ reshape(q'*sparseS, numGrids, numGrids));
        tmp2 .= sparseM .+ tmp;
        tmp3 .= sparseM .- tmp;
        A_mul_B!(forwardQ,tmp2,q);
        A_ldiv_B!(lufact!(tmp3),forwardQ);
    end
end

@time testFun(sparseM,sparseR,sparseS,q,10,799)

```

I might have misinterpreted what Julia is internally doing a little bit since `testFun2` is about 2x slower than `testFun` so there’s something going a little wrong. But those are the tools you’d use to find out so I’ll leave it at that. But what this is pointing to is that you might want to just be storing the transpose of that sparse matrix for this algorithm.

---

<div class="post-metadata">

### Author: ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [December 20, 2017, 8:10am UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840/24 "2017-12-20T08:10:05Z")

</div>

It is possible to use [https://github.com/JuliaSparse/MKLSparse.jl](https://github.com/JuliaSparse/MKLSparse.jl) to use MKL for sparse operations.

---

<div class="post-metadata">

### Author: ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [December 20, 2017, 10:13am UTC](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840/25 "2017-12-20T10:13:56Z")

</div>

Have you considered transposing S (before, not in the loop)?

```julia
n = 799; S = sprand(n*n,n,1.0-0.977); q=rand(n); qo=zeros(n*n);
@time A_mul_B!(qo,S,q);
  0.156302 seconds (4 allocations: 160 bytes)

n = 799; St = sprand(n,n*n,1.0-0.977); q=rand(n); qo=zeros(n*n);
@time At_mul_B!(qo,St,q);
  0.032524 seconds (4 allocations: 160 bytes)

```

Reason: q fits into L1, qo does not. I don’t know a lot about how julia handles sparse matrices, but we clearly want to operate with sequential writes and random reads here.

Edit: I overlooked it in the thread, you did consider it. Alternatively, you don’t really need to store S as a matrix. Just sort first by output-index, then by input-index, and implement the matrix multiplication by hand (trivial loop).

Edit2: At\_mul\_B! is faster than the trivial loop. So use that, I guess?

[Previous page](https://discourse.julialang.org/t/sparse-matrix-vector-product-much-more-slow-than-matlab/7840.md?page=1)