I was benchmarking different ways of summing a matrix by its second dimension and I got some surprising results (for me) when doing it through a for loop.
I had thought that in general the loop would be faster if the inner loop was going through the rows, so it loads the matrix one column at a time in memory. But given the example below, function f, which loops through columns first, is faster than its counterpart g.
I’m guessing this is because it’s faster to iterate through the vector res over and over again, than calculating the sum for each entry one at a time. Any intuition on why this is the case? Does the CPU guess to keep res in cache if you’re looping through it repeatedly?
julia> using BenchmarkTools
julia> X = randn(10000, 100); Xt = Matrix(X');
julia> function f(x)
n, m = size(x)
res = zeros(m)
for i = 1:n, j = 1:m
res[j] += x[i, j]
end
return res
end;
julia> function g(x)
n, m = size(x)
res = zeros(m)
for j = 1:m, i = 1:n
res[j] += x[i, j]
end
return res
end;
julia> function h(x)
m, n = size(x)
res = zeros(m)
for i = 1:n, j = 1:m
res[j] += x[j, i]
end
return res
end;
julia> function p(x)
m, n = size(x)
res = zeros(m)
for j = 1:m, i = 1:n
res[j] += x[j, i]
end
return res
end;
julia> @benchmark f($X)
BenchmarkTools.Trial: 7641 samples with 1 evaluation.
Range (min … max): 535.404 μs … 1.974 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 612.767 μs ┊ GC (median): 0.00%
Time (mean ± σ): 651.517 μs ± 124.330 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▂██▇▆▄▃▃▂▁ ▁
▂▃▃▄█████████████▇█▇█▇▇▇▇▇█▇▇▇▇▇▇███▇▇▇▇▇▅▆▆▆▇▆▄▅▃▅▃▃▄▅▅▂▄▄▅▄ █
535 μs Histogram: log(frequency) by time 1.22 ms <
Memory estimate: 896 bytes, allocs estimate: 1.
julia> @benchmark g($X)
BenchmarkTools.Trial: 5624 samples with 1 evaluation.
Range (min … max): 852.255 μs … 1.672 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 855.119 μs ┊ GC (median): 0.00%
Time (mean ± σ): 886.767 μs ± 66.944 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▂▁▃▁▁ ▁ ▂▁▃▄▁
████████████████████▇▇▆▆▆▆▆▅▅▅▅▆▅▅▄▅▇▅▅▄▄▄▂▅▄▄▄▄▄▅▂▄▄▄▃▄▃▂▂▄ █
852 μs Histogram: log(frequency) by time 1.19 ms <
Memory estimate: 896 bytes, allocs estimate: 1.
julia> @benchmark h($Xt)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 389.913 μs … 3.518 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 449.312 μs ┊ GC (median): 0.00%
Time (mean ± σ): 492.996 μs ± 204.491 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▃█▇▄▃▂▁▁▁▁▁ ▂
████████████████▇▆▅▆▆▆▅▃▅▄▄▅▄▅▄▅▆▄▅▄▄▅▅▄▅▅▄▅▄▄▄▅▅▃▆▅▅▁▄▅▄▃▅▄▃ █
390 μs Histogram: log(frequency) by time 1.65 ms <
Memory estimate: 896 bytes, allocs estimate: 1.
julia> @benchmark p($Xt)
BenchmarkTools.Trial: 3554 samples with 1 evaluation.
Range (min … max): 1.236 ms … 6.085 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.278 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.402 ms ± 428.580 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▇▆▄▃▂▂▁▂▁▁▁ ▁
█████████████▇█▇▇▆▇▆▆▇▆▆▅▅▆▅▅▅▁▆▆▅▅▅▅▅▆▁▄▅▆▅▆▅▄▅▄▃▆▃▃▁▃▃▁▁▅ █
1.24 ms Histogram: log(frequency) by time 3.29 ms <
Memory estimate: 896 bytes, allocs estimate: 1.
I find that g is faster on my computer.
The difference between the performance might be that your CPU is able to run fewer instructions in parallel if they all write to the same index of f.
Modern CPUs are superscalar, running multiple instructions at a time per core. In g, you write to the same res[j] multiple iterations after each other. This prevents the CPU from doing multiple iterations at a time, since each iteration has to wait for the last one. This is not the case for f or h.
You can get even faster than this if you enable SIMD.
It’s so odd that f is faster than h. The two are the same except for the order in which they access x, which should be much favorable in h since it is column-major order.