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:

```
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:

- Using the
`*`

operator, i.e. computing`v = X * β`

. - Using a for loop:

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

- Using a for loop, but iterating only over those indices for which
`β[i]`

is non-zero (stored in an integer array named`nonzero_indices`

).

```
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:

- 28.182 ms
- 1.272 μs
- 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!