# Performance of sparse matrix multiplied with a reshaped view

**URL:** <https://discourse.julialang.org/t/performance-of-sparse-matrix-multiplied-with-a-reshaped-view/44869>\
**Category:** New to Julia\
**Created:** [August 13, 2020, 12:13pm UTC](https://discourse.julialang.org/t/performance-of-sparse-matrix-multiplied-with-a-reshaped-view/44869 "2020-08-13T12:13:54Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ckoe-bccms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ckoe-bccms/32/1816_2.png) [@ckoe-bccms](https://discourse.julialang.org/u/ckoe-bccms)\
**Post date:** [August 13, 2020, 12:13pm UTC](https://discourse.julialang.org/t/performance-of-sparse-matrix-multiplied-with-a-reshaped-view/44869/1 "2020-08-13T12:13:54Z")

</div>

Hello everybody,

this question is very similar to other performance questions with sparse matrix multiplications where performance hits boil down to julia falling back to unoptimized methods. Still, I would like to ask for confirmation and eventually some hint how to regain performance.

First I would like to establish a baseline timing for the MWE:

```julia
julia> using LinearAlgebra, SparseArrays, BenchmarkTools
julia> function mymul(A,x)
         z=A*x
         return z
       end
mymul (generic function with 1 method)

julia> A=rand(100,100);
julia> B=sparse(rand(100,100));

julia> vec=rand(100);

julia> @btime Z=mymul($A,$vec);
  1.059 μs (1 allocation: 896 bytes)
julia> @btime Z=mymul($B,$vec);
  6.660 μs (1 allocation: 896 bytes)

```

Multiplying a sparse matrix with a dense vector is about six to seven times slower than multiplying a dense matrix with the same vector. That appears to be reasonable in this example.

Changing the vector to a view of a vector does not change the timings

```julia
julia> vec2=rand(200);
julia> vecview=@view vec2[1:100];

julia> @btime Z=mymul($A,$vecview);
  1.058 μs (1 allocation: 896 bytes)
julia> @btime Z=mymul($B,$vecview);
  6.882 μs (1 allocation: 896 bytes)

```

which is very nice. I was surprised how well the view works for either case, the overhead of a view is gone in the current version.

However, when using `reshape()` to make a vector form a view inside an array

```julia
julia> ar=rand(102,102);
julia> arview=@view ar[5:14,5:14];
julia> vec3=reshape(arview,:);

julia> @btime Z=mymul($A,$vec3);
  1.894 μs (1 allocation: 896 bytes)
julia> @btime Z=mymul($B,$vec3);
  141.315 μs (1 allocation: 896 bytes)

julia> @which A*vec3
*(A::AbstractArray{T,2}, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /home/ckoe/bin/julia-1.5.0/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:49
julia> @which B*vec3
*(A::AbstractArray{T,2}, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /home/ckoe/bin/julia-1.5.0/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:49

```

the performance hit in the dense matrix case is moderate while in the sparse case it is quite big.  
I would like to ask some questions about this.

Is the loss of performance in the dense case from the elements in `vec3` being no longer adjacent in memory ? Is this a case where for the given types the fallback method for the sparse matrix case is just slow ? This is not obvious to me from the output of `@which`.

Finally, is there a better way to construct the vector `vec3` ? Of course a `collect()` will remove the performance penalty but that creates a copy which I would like to avoid.

Best Regards

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [August 13, 2020, 12:32pm UTC](https://discourse.julialang.org/t/performance-of-sparse-matrix-multiplied-with-a-reshaped-view/44869/2 "2020-08-13T12:32:04Z")

</div>

> [@ckoe-bccms](#):
>
> Is the loss of performance in the dense case from the elements in `vec3` being no longer adjacent in memory ?

Yes, I think so:

```julia
julia> arview isa StridedArray
true

julia> vec3 isa StridedArray
false

```

The dispatch to `generic_matvecmul!` happens a few steps later, I think here: `@which mul!(similar(vec3), A, vec3, true, false) ` which you can get to with `@less` a lot, or `Dubugger.@enter`.

Instead of `collect(vec3)`, perhaps you can re-use some vectors (both for this and for the output) allocated outside the fast bit?

Why `B*vec3` is slower I’m not sure.

---

<div class="post-metadata">

**Author:** ![ckoe-bccms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ckoe-bccms/32/1816_2.png) [@ckoe-bccms](https://discourse.julialang.org/u/ckoe-bccms)\
**Post date:** [August 13, 2020, 1:16pm UTC](https://discourse.julialang.org/t/performance-of-sparse-matrix-multiplied-with-a-reshaped-view/44869/3 "2020-08-13T13:16:44Z")

</div>

> [@mcabbott](#):
>
> Yes, I think so:
> 
> ```julia
> julia> arview isa StridedArray
> true
> 
> julia> vec3 isa StridedArray
> false
> 
> ```
> 
> The dispatch to `generic_matvecmul!` happens a few steps later, I think here: `@which mul!(similar(vec3), A, vec3, true, false) ` which you can get to with `@less` a lot, or `Dubugger.@enter` .

Thank you for the confirmation. I will have to remember that it is possible to check if something is of a given type (if you now the type). I might try to drill deeper to understand what happens in the sparse case.

> [@mcabbott](#):
>
> Instead of `collect(vec3)` , perhaps you can re-use some vectors (both for this and for the output) allocated outside the fast bit?

Yes, passing a “scratch vector” around would be a pragmatic solution.