# Matrix-vector multiplication slower than a 'naive' for loop?

**URL:** <https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910>\
**Category:** Performance\
**Tags:** vector\
**Created:** [July 29, 2020, 7:30pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910 "2020-07-29T19:30:00Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![barankarakus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barankarakus/32/15938_2.png) [@barankarakus](https://discourse.julialang.org/u/barankarakus)\
**Post date:** [July 29, 2020, 7:30pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/1 "2020-07-29T19:30:00Z")

</div>

Hi,

I’m new to Julia and have been playing around with it to test its performance claims.

I am looking to multiply a matrix `X` by a sparse vector `β`, generated via the following code:

```julia
n = 500
p = 100000

β = rand(p); 
zero_indices = Int[]
nonzero_indices = Int[]
for i in 1:p
    if rand() > 0.3
        push!(zero_indices, i)
    else
        push!(nonzero_indices, i)
    end
end
for i in zero_indices
    β[i] = 0 
end 
X = rand(n, p);

```

I tried three ways to do this:

1. Using the `*` operator, i.e. computing `v = X * β`.
2. Using a for loop:

```julia
v = zeros(size(X)[1])
for i in 1:length(β)
    v += X[:, i] * β[i]
end

```

1. Using a for loop, but iterating only over those indices for which `β[i]` is non-zero (stored in an integer array named `nonzero_indices`).

```julia
v = zeros(size(X)[1])
for i in nonzero_indices
    v += X[:, i] * β[i]
end

```

I used the `@benchmark ` macro to time each of these approaches. Here are the mean times:

1. 28.182 ms
2. 1.272 μs
3. 1.132 μs

It makes sense to me that method 3 is faster than method 2: we’re looping over fewer things. I suspected it may even be faster than method 1, since we’re leveraging the sparsity of the vector β. But why on earth is method 2 so much faster than method 1? That is, why is a for loop implementation of matrix-vector multiplication so much faster than an operator dedicated to this computation?

Thanks!

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [July 29, 2020, 7:52pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/2 "2020-07-29T19:52:53Z")

</div>

Double-check your benchmarking code. I’m getting:

```julia
julia> @btime f1($X, $β);
  19.875 ms (1 allocation: 4.06 KiB)

julia> @btime f2($X, $β);
  186.214 ms (300001 allocations: 1.16 GiB)

julia> @btime f3($X, $β, $nonzero_indices);
  56.557 ms (89530 allocations: 355.19 MiB)

```

which makes sense given your observations. My function definitions are:

```julia
julia> function f1(X, β)
         X * β
       end
f1 (generic function with 1 method)

julia> function f2(X, β)
         v = zeros(size(X)[1])
         for i in 1:length(β)
           v += X[:, i] * β[i]
         end
         v
       end
f2 (generic function with 1 method)

julia> function f3(X, β, nonzero_indices)
         v = zeros(size(X)[1])
         for i in nonzero_indices
           v += X[:, i] * β[i]
         end
         v
       end

```

Now if you want a _faster_ implementation, you just need to follow the [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html). Here’s an easy improvement of `f3` which makes it faster than `f1` by more than a factor of 2:

```julia
julia> function f4(X, β, nonzero_indices)
         v = zeros(size(X, 1))
         for i in nonzero_indices
           v .+= @view(X[:, i]) .* β[i]
         end
         v
       end
f4 (generic function with 1 method)

julia> @btime f4($X, $β, $nonzero_indices);
  7.871 ms (29844 allocations: 1.37 MiB)

```

By the way, Julia 1.5 (released soon) will make `f4` even faster by avoiding memory allocation for the views:

```julia
julia-1.5> @btime f4($X, $β, $nonzero_indices);
  7.289 ms (1 allocation: 4.06 KiB)

```

---

<div class="post-metadata">

**Author:** ![barankarakus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barankarakus/32/15938_2.png) [@barankarakus](https://discourse.julialang.org/u/barankarakus)\
**Post date:** [July 29, 2020, 8:03pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/3 "2020-07-29T20:03:11Z")

</div>

Thanks!

Your timing results sure make more sense. I re-ran my experiment but unfortunately I’m still getting similar results. I generated the matrix and vector exactly as above, and ran the following cells in a Jupyter notebook:

 ![Screenshot 2020-07-29 at 21.00.57](https://global.discourse-cdn.com/julialang/original/3X/e/4/e46b7b07ee3b83745a715a20f32424e946be1c73.png)

Do you have any idea what might be going wrong here?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [July 29, 2020, 8:10pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/4 "2020-07-29T20:10:01Z")

</div>

`quote` creates an _expression_ when metaprogramming. It doesn’t actually _do_ anything with that expression. You’re benchmarking essentially the time it takes to parse the code down into an expression object, not how long it takes to _run_ it.

Get rid of the `quote` and use `let` or `begin` to introduce a new block. Even better, use an actual function (see the [very first performance tip](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Avoid-global-variables-1)).

---

<div class="post-metadata">

**Author:** ![barankarakus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barankarakus/32/15938_2.png) [@barankarakus](https://discourse.julialang.org/u/barankarakus)\
**Post date:** [July 29, 2020, 8:10pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/5 "2020-07-29T20:10:38Z")

</div>

Ah I see… Very silly of me 🤣! Thank you.

---

<div class="post-metadata">

**Author:** ![barankarakus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barankarakus/32/15938_2.png) [@barankarakus](https://discourse.julialang.org/u/barankarakus)\
**Post date:** [July 29, 2020, 9:43pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/6 "2020-07-29T21:43:21Z")

</div>

> [@rdeits](#):
>
> Now if you want a _faster_ implementation, you just need to follow the [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html). Here’s an easy improvement of `f3` which makes it faster than `f1` by more than a factor of 2:
> 
> ```julia
> julia> function f4(X, β, nonzero_indices)
> v = zeros(size(X, 1))
> for i in nonzero_indices
> v .+= @view(X[:, i]) .* β[i]
> end
> v
> end
> f4 (generic function with 1 method)
> 
> julia> @btime f4($X, $β, $nonzero_indices);
> 7.871 ms (29844 allocations: 1.37 MiB)
> 
> ```
> 
> By the way, Julia 1.5 (released soon) will make `f4` even faster by avoiding memory allocation for the views:
> 
> ```julia
> julia-1.5> @btime f4($X, $β, $nonzero_indices);
> 7.289 ms (1 allocation: 4.06 KiB)
> 
> ```

This is a huge reduction in memory. From the [documentation](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Consider-using-views-for-slices-1) I understand that the point of views is to avoid copying the data, instead referencing it in-place. I’m curious to know why these views involve so much memory allocation - and why they don’t in Julia 1.5. Would you be able to enlighten? 🙂

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [July 29, 2020, 10:02pm UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/7 "2020-07-29T22:02:05Z")

</div>

A view is a lightweight reference to another array, so when you create a view you don’t have to copy the _data_ in that array, but you still have to construct the view itself. In Julia versions before 1.5, the view itself was often allocated on the heap, so code using views would still show a large number of allocations, but each allocation would be very small (just a pointer to the original array and some extra information about the size). That’s what happened here. Compare `f3` (no views) which had:

```julia
89530 allocations: 355.19 MiB

```

with `f4` (using views) which had:

```julia
29844 allocations: 1.37 MiB

```

The number of allocations is only reduced by about 3X, but the total amount of allocated memory is reduced by 300X because we’re no longer copying the data for each slice of `X`. The average allocation size (1.37 MiB / 29844) is less than 50 bytes, which is tiny (just a few `Int64`s or pointers). That makes sense since we’re allocating a large number of lightweight views.

The improvement in Julia 1.5 is that the view itself can now be constructed on the stack and therefore does not require memory allocation. The details are pretty low-level, but the relevant PR is here: [https://github.com/JuliaLang/julia/pull/33886](https://github.com/JuliaLang/julia/pull/33886)

---

<div class="post-metadata">

**Author:** ![nhavt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nhavt/32/16356_2.png) [@nhavt](https://discourse.julialang.org/u/nhavt)\
**Post date:** [July 30, 2020, 8:29am UTC](https://discourse.julialang.org/t/matrix-vector-multiplication-slower-than-a-naive-for-loop/43910/8 "2020-07-30T08:29:20Z")

</div>

Hi, I propose another version which is as fast as f4 and requires fewer allocations.

```julia
function f5(X, β, nonzero_indices)
    v = zeros(size(X, 1))
    @inbounds @simd for i in nonzero_indices
                      for j in 1:size(X,1)
                        v[j] += X[j,i] * β[i]
                      end
                    end
    v
end

```

Below is the computing time with my computer:

```julia
@btime f4($X, $β, $nonzero_indices);
8.376 ms (29681 allocations: 1.36 MiB)
@btime f5($X, $β, $nonzero_indices);
8.288 ms (1 allocation: 4.06 KiB)

```