Is the timing of MATLAB reliable? If yes, can we reproduce the performance with julia, fortran, etc.?

Just to add on the Julia side - make sure you use BenchmarkToolsto benchmark, wrap the code you want to benchmark in functions so as not to benchmark in global scope, and interpolate any variables you pass to @btime.

Here's how I would do it:

julia> s = 15_000;

julia> function f_loop!(H)
           for c ∈ 1:size(H, 1)
               for r ∈ 1:size(H, 1)
                   H[r, c] = 1 / (r + c - 1)
               end
           end
       end
f_loop! (generic function with 1 method)

julia> function f_vec!(H)
           c = 1:size(H, 1)
           r = c'
           H .= 1 ./ (r .+ c .- 1)
       end
f_vec! (generic function with 1 method)

julia> H = zeros(s, s);

julia> using BenchmarkTools

julia> @btime f_loop!($H);
  625.891 ms (0 allocations: 0 bytes)

julia> H = zeros(s, s);

julia> @btime f_vec!($H);
  625.248 ms (0 allocations: 0 bytes)

So both versions come in at the same time, which is what I'd expect for such a straightforward operation where a properly type-inferred code should compile down to roughly the same machine code.

tic/toc should be fine, but it looks like the timing is being skewed by memory pre-allocation.

I can reproduce similar timings to your MATLAB example, however

• On first run (clear workspace)

  • Loop approach takes 2.08 sec
  • Vectorised approach takes 1.04 sec
  • Vectorisation saves 50% execution time

• On second run (workspace not cleared)

  • Loop approach takes 2.55 sec
  • Vectorised approach takes 0.065 sec
  • Vectorisation "saves" 97.5% execution time

My guess would be that since the loop approach explicitly creates a new matrix via zeros, the memory is reallocated from scratch on every run and you don't see the speed improvement on subsequent runs.

However, when HH remains in memory and the HH=___ line outputs a matrix of the same size, I suspect MATLAB is doing some clever memory allocation to speed up the operation.

We can prove this theory with the following test:

Test Num  |  Workspace cleared  |    s    |  Loop (sec)  |  Vectorised (sec) 
    1     |       Yes           |  15000  |    2.10      |       1.41
    2     |        No           |  15000  |    2.73      |       0.07
    3     |        No           |  15000  |    2.50      |       0.07
    4     |        No           |  15001  |    2.74      |       1.73

See the variation between tests 2 and 3, this is why timeit would have been helpful for an average runtime (see footnote). The difference in output sizes between tests 3 and 4 are pretty small, but the execution time returns to a similar magnitude of that in test 1 for the vectorised approach, suggesting that the re-allocation to create HH costs most of the time.

Footnote: tic/toc timings in MATLAB can be improved by using the in-built timeit function, which essentially takes an average over several runs. One interesting thing to observe from the workings of timeit though is that it explicitly "warms up" (quoting a comment) the tic/toc function by calling it a couple of times. You can see when running tic/toc a few times from a clear workspace (with no intermediate code) that the first call takes longer than subsequent calls, as there must be some overhead for getting the timer initialised.

I hope that the following modified benchmark could bring some new light to the problem:

s = 15000;

tic
% for-loop version
H = zeros(s,s);
for i =1:10
  for c = 1:s
    for r = 1:s
        H(r,c) = H(r,c) + 1/(r+c-1+i);
    end
  end
end
toc

tic;
% vectorized version
HH = zeros(s,s);
c = 1:s;
r = c';
for i=1:10
   HH= HH + 1./(r+c-1+i);
end
toc

isequal(H,HH)

In this case any kind of "cashing" is avoided by changing of matrix H (HH) at each for-loop (over "i") iteration.

In this case we get:

Elapsed time is 3.737275 seconds. (for-loop)
Elapsed time is 1.143387 seconds. (vectorized)

So, there is still performance improvement (~ 3x) due to the vectorization, which is probably done by implicit multi-threading implementation of vectorized Matlab commands.

Yes, tic/toc vs timeit is not strictly consistent, but the overall timing functionality is very similar.

To add to this, here is a simple python script which does the vectorized operation with numpy:

from timeit import default_timer
import numpy as np

s = 15000

start = default_timer()
# for-loop
H = np.zeros([s, s])
for c in range(1, s):
    for r in range(1, s):
        H[r, c] = 1 / (r + c - 1)
end = default_timer()
print(end - start)

start = default_timer()
# vectorized
c = np.arange(1, s).reshape([1, -1])
r = c.T
HH = 1 / (c + r - 1)
end = default_timer()
print(end - start)

for-loop: 32.94566780002788 seconds
vectorized: 0.494859800033737 seconds

While the for-loop version is terribly slow, the vectorized version is faster than the posted fortran/julia times. Numpy internally tries to use special SIMD hardware instructions to speed up arithmetic on vectors, which can make a significant difference. It's possible that the fortran/julia compilers weren't able to generate those instructions from the provided code, but numpy/matlab were able to. However, Matlab is still about 10x faster than the numpy code, which I don't think would be explained by better use of SIMD instructions. Instead, they may also be using multiple threads to parallelize the computation, since the matrix is fairly large.

Ultimately, I think the matlab numbers are plausible, but I'm not sure exactly how they're getting their speedup.