Speed up Do loops in MMA 13.0

Question

How to speed up Do loop in MMA 13.

We consider the following Benchmark test (4 kernels, i7, win10):

In MMA 13:

s = 5000;
Hmm = ConstantArray[0, {s, s}];

Do[Do[Hmm[[r, c]] = 1/(r + c - 1), {r, s}], {c, s}] // AbsoluteTiming

It takes 36.8516356 seconds.

But in Matlab 2021b

s = 5000;
H = zeros(s,s);

tic 
for c = 1:s
    for r = 1:s
        H(r,c) = 1/(r+c-1);
    end
end
toc

It only takes 0.114233 seconds......

Nearly 360 times slower than Matlab 2021b...

Update:

1.) If we use "Table"

s = 5000;
Hmm = ConstantArray[0, {s, s}];

AbsoluteTiming[
 Table[Table[Hmm[[r, c]] = 1/(r + c - 1), {r, 1, s}], {c, 1, s}]]

It takes 36.6726 seconds...

2.) If we use "For"

AbsoluteTiming[
 For[c = 1, c <= s, c++, 
  For[r = 1, r <= s, r++, Hmm[[r, c]] = 1/(r + c - 1)]]]

It takes 46.7529 seconds...

3.) Test results from Matlab 2021b

4.) If we try "Compile" (https://mathematica.stackexchange.com/a/261329/54516)

Compile[{}, Module[{s, Hmm}, s = 5000;
     Hmm = Table[0., s, s];
     Do[Do[Hmm[[r, c]] = 1/(r + c - 1), {r, s}], {c, s}]; 
     Hmm]][]; // AbsoluteTiming

It takes 1.0638 seconds...

Nearly 10 times slower than Matlab 2021b...

5.) If we try another "Compile" (https://mathematica.stackexchange.com/a/261329/54516)

cf0 = With[{s = s}, 
    Compile[{}, Table[1/(r + c - 1), {r, 1, s}, {c, 1, s}], 
     CompilationTarget -> "C", RuntimeOptions -> "Speed"]][[-1]];

Hmm = cf0[]; // AbsoluteTiming

It takes 0.181941 seconds...

Nearly 2 times slower than Matlab 2021b...

Note that, for this special case: MATLAB and Mathematica are NOT equally fast.

6.) Why is tic/toc used (@xzczd's Question)?

Because e.g. "Use a pair of tic and toc calls to report the total time required for element-by-element matrix multiplication; use another pair to report the total runtime of your program." (https://www.mathworks.com/help/matlab/ref/toc.html)

Please check: https://www.mathworks.com/help/matlab/ref/toc.html

7.) How about julia 1.6.3 Do loops speed

@time Hmm=[1. /(r+c-1) for r=1:s,c=1:s];
    #  0.107591 seconds (85.06 k allocations: 195.439 MiB, 44.46% compilation time)

from @xzczd: (https://mathematica.stackexchange.com/a/261329/54516):

It takes 0.107591 seconds... @xzczd.

8.) The computational performance of the @chyanog's MMA code (@chyanog's comments https://mathematica.stackexchange.com/a/261329/54516)

s = 5000; 
 cf = With[{s = s}, 
   Compile[{{r, _Integer}}, Table[1/(r + c - 1), {c, 1, s}], 
    CompilationTarget -> "C", RuntimeOptions -> "Speed", 
    RuntimeAttributes -> {Listable}]]; 
 Hmm = cf[Range[s]]; // AbsoluteTiming

It takes 0.0717162 seconds... @chyanog.

Nearly 1.5 times faster than Matlab 2021b...

9.) "ParallelTable"

Based on the update 8.), now we test the ParallelTable :

cf = With[{s = s}, 
   Compile[{{r, _Integer}}, ParallelTable[1/(r + c - 1), {c, 1, s}], 
    CompilationTarget -> "C", RuntimeOptions -> "Speed", 
    RuntimeAttributes -> {Listable}]];

Hmm = cf[Range[s]]; // AbsoluteTiming

Answer 1

Here's the fastest I've found:

foo = Divide[1., 
    Outer[Plus, Range[1., s], Range[0., s - 1]]]; // RepeatedTiming
(*  {0.175878, Null}  *)

foo == Hmm
(*  True  *)

For comparison, moving H = zeros(s,s) inside tic..toc, the MATLAB timing on my machine is 0.166379.

---

Addendum: Notes

The trouble with the OP's first code is that Hmm has to be unpacked when 1/(r+c-1) is a rational number and not an integer. Yes, I mean integer because the preallocation was an array of integers. Do[..] and Table are still disappointingly slow even with the proper preallocation and formula.

Answer 2

Hmm = 
   With[{s = s}, 
     Compile[{}, Table[1/(r + c - 1), {r, 1, s}, {c, 1, s}]]][]; // AbsoluteTiming
(* {0.707614, Null} *)

If you insist on Do:

Compile[{}, Module[{s, Hmm},
     s = 5000;
     Hmm = Table[0., s, s];
     Do[Do[Hmm[[r, c]] = 1/(r + c - 1), {r, s}], {c, s}]; Hmm](*,
    CompilationTarget -> "C", RuntimeOptions -> "Speed"*)][]; // AbsoluteTiming
(* {0.9933, Null} *)

For comparison, Table[Table[Hmm[[r, c]] = 1/(r + c - 1), {r, 1, s}], {c, 1, s}]; takes 37 seconds on my laptop.

---

Since the speed of MATLAB is suspiciously fast, I tested the following sample with gfortran (TDM-GCC-10.3.0-2, with compile option -Ofast):

program tst
use, intrinsic :: iso_fortran_env
implicit none
integer,parameter::s=5000
integer::r,c
real(real64)::hmn(s,s)
do r=1,s
    do c=1, s
        hmn(r,c)=1._real64/(r + c - 1)
    end do 
end do
print *, hmn(s,s)
end program

The compilation timing is 0.2152324 second and execution timing is 0.1160831 second. I know little about MATLAB, but perhaps the tic/toc of MATLAB doesn't count the compilation timing. I've also tested with julia 1.6.3:

@time (s=5000; Hmm=[1. /(r+c-1) for r=1:s,c=1:s];)
#  0.107591 seconds (85.06 k allocations: 195.439 MiB, 44.46% compilation time)

Anyway, if we take away the compilation timing and extract the LibraryFunction[…] as mentioned here:

cf0 = With[{s = s}, 
    Compile[{}, Table[1/(r + c - 1), {r, 1, s}, {c, 1, s}], CompilationTarget -> "C", 
     RuntimeOptions -> "Speed"]][[-1]];

Hmm = cf0[]; // AbsoluteTiming
(* {0.0954907, Null} *)

Now MATLAB and Mathematica are equally fast. Still, my compiler is TDM-GCC-10.3.0-2, with compile option -Ofast.

---

Aha, FunctionCompile performs slightly better than fortran, if only compilation timing doesn't need to be counted!:

cf = 
   FunctionCompile[
    Function[Typed[s, "MachineInteger"], 
     Table[1./(r + c - 1), {r, 1, s}, {c, 1, s}]]]; // AbsoluteTiming
(* {3.47833, Null} *)

cf[5000]; // RepeatedTiming
(* {0.0772538, Null} *)

All the tests are done with Mathematica 12.3.1, Windows 10, on a laptop with i7-8565U CPU.

---

I believe s=5000 is too small for accurate timing, so I turn to s=15000, the result is as follows:

fortran
compilation timing: 0.2057823
execution timing: 0.7179657

julia (Timing given by @time of julia oscillates quite a bit, the following is an average of 10 measurements. )
0.7945998

cf (FunctionCompile of Mathematica)
0.643739

Answer 3

As other answers have suggested, a better approach here is to stick to high efficient functional way, and try to take advantage of "vectorizable" operation when possible.

For OP's question, one possible way is (time measured on a laptop with Intel i7):

AbsoluteTiming[hmm = Range[5000] // N // Outer[Plus, #, #] & // 1/(# - 1) &;]
(* Out[] = {0.195207, Null} *)
hmm // Developer`PackedArrayQ
(* Out[] = True *)

Here Range[5000] // N creates a PackedArray that supports efficient storage and computing. Outer invokes auto-compilation with internal typed-version which is also very efficient. At last 1/(# - 1) & uses vectorized version of Plus, Power, Times, etc. to transform the inputted array "as a whole".

To compare, try this 1/100 size example without vectorization, yet it's slower than above code:

AbsoluteTiming[Range[500] // N // Outer[1/(#1 + #2 - 1) &, #, #] &;]
(* Out[] = {0.276201, Null} *)

Answer 4

Much of the time may be spent in formatting the output, that is not shown because it is very large. Also, it may be possible to reserve space for hmm, but creating a big null matrix and redefining it does not work. This is a simple way to improve the timing:

s=5000
line = Range[0, s - 1]
lineJ[j_] := (1./(# + j) & /@ line)
hmm = lineJ /@ Range[1, s] // AbsoluteTiming; (* this ; supresses output *)
hmm[[1]]

yields 1.28 s in my notebook (i7 on windows 11). Notice also that in lineJ[ ] the numerator is 1. (a real number), not 1 (an Integer); with 1, the matrix will be constituted of rationals and will take much longer to build (21 s), which is not the case with MatLab. A faster way is

rmm = (1.0/Range[#, # + s - 1]) & /@ Range[1, s] // AbsoluteTiming;

which gives 0.42 s.

Answer 5

As a general rule, avoid Do. In this case, Table is more efficient

Clear["Global`*"]

s = 5000;

The matrix initialization is part of defining the matrix with Do and should be included in the timing.

(Hmm = ConstantArray[0, {s, s}];
  Do[Do[Hmm[[r, c]] = 1/(r + c - 1), {r, s}], {c, s}]) // AbsoluteTiming

(* {31.1698, Null} *)

Note that it is more efficient to use to use a single Do rather than nested Do

(Hmm2 = ConstantArray[0, {s, s}];
  Do[Hmm2[[r, c]] = 1/(r + c - 1), {r, s}, {c, s}]) // AbsoluteTiming

(* {28.4179, Null} *)

When using Table, the assignment should be done only once, i.e., at the matrix level rather than the element level.

(Hmm3 = Table[1/(r + c - 1), {r, s}, {c, s}];) // AbsoluteTiming

(* {16.9393, Null} *)

Even Array is faster than Do

(Hmm4 = Array[1/(#1 + #2 - 1) &, {s, s}];) // AbsoluteTiming

(* {24.1783, Null} *)

Verifying the equivalence of the different approaches,

Hmm === Hmm2 === Hmm3 === Hmm4

(* True *)

I cannot speak to MatLab, as I do not use it (and don't have access to it).

Answer 6

Matlab code is not fully optimized!!! See following code:

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

tic;
% vectorized version
c = 1:s;
r = c';
HH=1./(r+c-1);
toc

isequal(H,HH)

Vectorized Matlab code is much more faster (cca 25x) than simple for-loop.

Elapsed time is 0.168488 seconds.... for-lop
Elapsed time is 0.007221 seconds.... vectorized

So, in general, Matlab is in floating point numerics always significantly faster than Mathematica. Moreover, Matlab multi-threading support is much more better for majority of Matlab built-in functions.

Answer 7

Another way is to precompute all quotients and then extract each row:

cf3 = FunctionCompile[
  Function[Typed[s, "MachineInteger"], 
    With[{vals = Divide[1.0, Range[2 s - 1]]},
      Table[vals[[k ;; k + s - 1]], {k, s}]
    ]
  ]
];

cf3[5000]; // RepeatedTiming

{0.0388719, Null}

Answer 8

Sorry, it's not the answer, just a long comment.

Considering the Outer version in this manner:

time WolframKernel <<EOF
With[{s = 15000}, 
 Table[AbsoluteTiming[
    Divide[1., Outer[Plus, Range[1., s], Range[0., s - 1]]]][[1]], 30]]
EOF

gives me

real    0m17.796s
user    0m56.334s
sys 1m18.325s

So the system spends more time in system calls than actually doing something. strace ing shows that sys calls are sched_yield so we're dealing with OMP usage here. I think this should be addressed by WRI