# Speed up Do loops in MMA 13.0

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)

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


• Then what's the first-time cost of your MATLAB program? Dec 29 '21 at 4:22
• @xzczd in this post. Dec 29 '21 at 4:24
• 1. Notice I'm compiling with -Ofast option of TDM-GCC, have you set the compiler properly? 2. By first-time cost I mean timing with a fresh kernel of MATLAB 3. s=5000 is small, what's the timing of larger s? Dec 29 '21 at 4:35
• @xzczd if your conclusion is "Now MATLAB and Mathematica are equally fast." please show us the run time from MATLAB in your answer! It can be very helpful! Dec 29 '21 at 4:37
• I don't have access to MATLAB, but I've shown the timing of uncompiled Table and gfortran, which should be enough for comparison. Dec 29 '21 at 4:40

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.

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.

• While I'm typing! :D Dec 29 '21 at 5:30
• @Silvia Yes, I just saw yours. :) Dec 29 '21 at 5:32
• For Do way, I thought of pre-allocating as a NumericArray, but couldn't find handy way to mutate it other than C interface in WolframNumericArrayLibrary.h. Dec 29 '21 at 5:34
• @Silvia If you use Do[], then this seems a lower limit, Do[r, {r, 5000}, {c, 5000}] // AbsoluteTiming, which takes almost 1.8 sec. on my computer. Inside Compile[], it takes only 0.1 sec. Unless there's a trick to speed up Do. Dec 29 '21 at 6:11
• @ABCDEMMM I could say, "Because the loop is so much slower in Mathematica," and, although a non-answer, it expresses the limit of my knowledge. :) Do[] is has some generality, so I would expect some overhead, but this seems a lot more than necessary to me at present. Dec 29 '21 at 18:56
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

• Perhaps OP could print UNIX time before and after the tic/toc and subtract to see if tic/toc does indeed ignore compile time or not. Dec 29 '21 at 4:23
• default real is 32bit in fortran
– I.M.
Dec 29 '21 at 16:22
• @I.M. Thx for pointing out. (I know little about fortran, too. ) The sample is revised, and the timing increases to 0.1160831 second now. Dec 29 '21 at 16:41
• This version is a little faster. s=15000;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 Dec 30 '21 at 2:55
• @xzczd the Julia code might be faster as Hmm = 1 ./ ((1:s) .+ transpose(1:s) .- 1) Dec 31 '21 at 14:15

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 // DeveloperPackedArrayQ
(* 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} *)


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.

• Taking 1.0/x to Divide[1.0, x] gives a 5.5x speedup on my machine. See mathematica.stackexchange.com/q/40927 for more details. Dec 29 '21 at 4:28

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).

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.

• I believe s=5000 is too small, what if you test with s=15000? (See my post for timing of fortran and julia. ) Dec 30 '21 at 2:12
• Vectorization is in Matlab very well supported by multi-threaded functions. But, I am sure that Matlab timing is very reliable and precise. So, this Matlab performance is definitely not any kind of artefact. Dec 30 '21 at 11:48
• Ref the stackoverflow question by xzczd, it looks like MATLAB is inadvertently "cheating" here, and the majority of the vectorisation time savings here are actually due to avoiding a cold-start on memory allocation - running this from a clear workspace shows only a 50% speed up (approx) for vectorisation. Still worth having but nowhere near the stated advantage shown here. Dec 30 '21 at 12:36
• @Pavel As noted in my answer to the linked question, I think it's just memory pre-allocation. The expression/calculation itself isn't cached, but when re-running the same command without clearing the already-assigned variables the required memory for the output array is already allocated in one of the cases, meaning this (time consuming) step can be skipped. Dec 30 '21 at 14:02
• @Pavel The MATLAB license fee lets you leave intuition at the door and does a lot of heavy lifting to justify its cost ;) Dec 30 '21 at 16:05

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}

• Literally translating code is indeed always a bad indeed. This is a very nice way of doing it in a very smart way, and very fast! Jan 1 at 20:28

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