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)) 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
v = zeros(size(X)) 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?