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

@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! — ABCDEMMM, Dec 29 '21 at 04: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. — xzczd, Dec 29 '21 at 04:40
For a fair comparison, you should include H = zeros(s,s) inside tic..toc, imo. — Michael E2, Dec 29 '21 at 04:44
@MichaelE2 in such case: >> untitled Elapsed time is 0.125228 seconds.

untitled

Elapsed time is 0.122755 seconds.

untitled

Elapsed time is 0.128294 seconds. — ABCDEMMM, Dec 29 '21 at 04:46
Another hint: Do[Do[Hmm[[r, c]] = 1/(r + c - 1), {r, 1}], {c, 2}] // AbsoluteTiming. (Basically, preallocation was pointless.) — Michael E2, Dec 29 '21 at 04:47

Michael E2 · Accepted Answer · 2021-12-29T05:26:24.523

28

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.

edited Dec 29 '21 at 05:26

answered Dec 29 '21 at 05:20

Michael E2

235,386
17
334
747

While I'm typing! :D – Silvia Dec 29 '21 at 05:30
@Silvia Yes, I just saw yours. :) – Michael E2 Dec 29 '21 at 05: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. – Silvia Dec 29 '21 at 05:34
2

@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. – Michael E2 Dec 29 '21 at 06:11
@MichaelE2 why is "DO" loop in Matlab so much faster (nearly 360 times)? Here: for <-> Fo[] – ABCDEMMM Dec 29 '21 at 18:49
1

@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. – Michael E2 Dec 29 '21 at 18:56

xzczd · Answer 2 · 2021-12-30T13:56:21.967

21

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

edited Dec 30 '21 at 13:56

answered Dec 29 '21 at 01:57

xzczd

65,995
9
163
468

1

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. – Greg Hurst Dec 29 '21 at 04: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. – xzczd Dec 29 '21 at 16:41
@xzczd matlab solver does not run in fortran... it runs in c/c++... also with MKL .... – ABCDEMMM Dec 29 '21 at 19:20
@a I use fortran for estimating the limit of timing because fortran is famous for its speed. – xzczd Dec 30 '21 at 00:10
@xzczd but C/C++ is also famous for its speed... e.g. Eigen (C++ library) – ABCDEMMM Dec 30 '21 at 01:28
@a AFAIK, the traditional view is, they're almost equally fast. – xzczd Dec 30 '21 at 01:32
@xzczd is there a good way to call Eigen (C++ library) in MMA? – ABCDEMMM Dec 30 '21 at 01:39
Let us continue this discussion in chat. – xzczd Dec 30 '21 at 01:40
1

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 – chyanog Dec 30 '21 at 02:55
@chyanog Wow, why not posting it as a separate answer? – xzczd Dec 30 '21 at 03:06
@xzczd see the new update! – ABCDEMMM Dec 30 '21 at 10:09
@ABCDEMMM Once again, I'm afraid this is related to your compiler. On my laptop the timing is around 0.50517. – xzczd Dec 30 '21 at 10:21
@xzczd if we directly run his code, it takes such time, but only run "Hmm=cf[Range[s]]" , it is very fast, see the new update! – ABCDEMMM Dec 30 '21 at 10:34
1

@xzczd the Julia code might be faster as Hmm = 1 ./ ((1:s) .+ transpose(1:s) .- 1) – Simon Woods Dec 31 '21 at 14:15
@Simon Not obvious, but it's indeed slightly faster. (This is a bit suprising to me actually, I remember the so-called devectorization is encouraged around julia v0.2! ) – xzczd Dec 31 '21 at 15:16
@xzczd how is the simple loop-performance in MMA-GPU? also very slow? – ABCDEMMM Jan 04 '22 at 21:26
@a I don't have experience with GPU, but AFAIK Do loop is always not fast if not compiled. – xzczd Jan 05 '22 at 02:13

score 19 · Answer 3 · answered Dec 29 '21 at 05:29

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} *)

score 13 · Answer 4 · answered Dec 29 '21 at 02:05

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 https://mathematica.stackexchange.com/q/40927 for more details. — Greg Hurst, Dec 29 '21 at 04:28

score 11 · Answer 5 · answered Dec 29 '21 at 01:54

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

michalkvasnicka · Answer 6 · 2021-12-30T08:16:02.477

10

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.

edited Dec 30 '21 at 08:16

answered Dec 29 '21 at 15:49

michalkvasnicka

101
3

I believe s=5000 is too small, what if you test with s=15000? (See my post for timing of fortran and julia. ) – xzczd Dec 30 '21 at 02:12
Ok, let's talk only about Mathematica and Matlab. Matlab vectorized code is about 25x faster than its for-loop version. I think for benchmark purposes should be taken the vectorized version. – michalkvasnicka Dec 30 '21 at 08:20
Once again, what if you test with s=15000? – xzczd Dec 30 '21 at 08:30
Elapsed time is 1.359625 seconds.... For-loop Elapsed time is 0.047916 seconds.... Vectorized – michalkvasnicka Dec 30 '21 at 10:02
@xzczd can we have such performance (multi-threading) for this kind of simulation in MMA? – ABCDEMMM Dec 30 '21 at 10:12
Then again I feel the timing of MATLAB suspicious, because the timing of fortran and julia seems to suggest the limit of s=15000 is around 0.7 second for an average PC, and a 10x speed-up is just incredible in my eyes, but who knows. To @a : I don't think parallelization will help in this case, the communication is too expensive, have you read this?: https://mathematica.stackexchange.com/q/48295/1871 – xzczd Dec 30 '21 at 10:29
@xzczd I added paralleltable results, see the update! – ABCDEMMM Dec 30 '21 at 10:47
@a ParallelTable cannot be compiled, see this post for more info: https://mathematica.stackexchange.com/q/1096/1871 – xzczd Dec 30 '21 at 10:52
1

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. – michalkvasnicka Dec 30 '21 at 11:48
For s=20000, still relevant Matlab timing results corresponding to expexted complexity. Elapsed time is 2.526585 seconds (for-loop) Elapsed time is 0.100607 seconds (vectorized) – michalkvasnicka Dec 30 '21 at 11:55
Can you add the hardware parameters of your PC to the answer? BTW, since this is no longer a pure Mathematica question, I started a new question in stackoverflow: https://stackoverflow.com/q/70531385/1551513 – xzczd Dec 30 '21 at 12:03
1

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. – Wolfie Dec 30 '21 at 12:36
@Wolfie so we're dealing with some kind of expression cache here? – Pavel Perikov Dec 30 '21 at 13:27
I am afraid that the situation is not so clear as I expected. All previous timings I performed on Cloud Matlab service Matlab OnLine (virtualized Matlab instances). On regular 8core PC (INTEL i9 CPU) is the situation different: (s = 15000) Elapsed time is 0.477092 seconds. ... (for-loop) Elapsed time is 0.106886 seconds. ... (vectorized) – michalkvasnicka Dec 30 '21 at 13:34
1

@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. – Wolfie Dec 30 '21 at 14:02
@Wolfie Intuitively I don't think any kind of memory pre allocation can give such speedup but I may be wrong – Pavel Perikov Dec 30 '21 at 14:47
1

@Pavel The MATLAB license fee lets you leave intuition at the door and does a lot of heavy lifting to justify its cost ;) – Wolfie Dec 30 '21 at 16:05

score 9 · Answer 7 · answered Dec 31 '21 at 05:44

9

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}

answered Dec 31 '21 at 05:44

Greg Hurst

35,921
1
90
136

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! – SHuisman Jan 01 '22 at 20:28

Pavel Perikov · Answer 8 · 2021-12-30T14:49:09.920

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

Speed up Do loops in MMA 13.0

8 Answers8

Linked