# How to speed up nested loop? I have used Do loop and Table as well

I am trying to calculate mean for example if list is {1,2,3,4,5} then i have to calculate mean of Cos(1+2)+Cos(1+3)+Cos(1+4)+Cos(1+5)+Cos(2+3)+Cos(2+4)+Cos(2+5)+Cos(3+4)+Cos(3+5)+Cos(4+5)

My simple code is

     SetDirectory["/Users/macos/"]

phia = Flatten[Import["p180.dat", "Table"]];(* for negative charges*)
Length[phia]
res = 0;
l = 0;
Do[
Do[
res = res + Cos[phia[[m]] + phia[[n]]];
l++;
, {m, 2, Length[phia]}];
, {n, Length[phia]-1}]; // AbsoluteTiming
mean=res/l;


it takes 5943 seconds. it is very slow. how to speed this process?[![enter image description here][2]][2] thanking in advance i also used Table command that also does not helped much.

• Iterator specs to Table and Do can be done all in a single call. Your performance issue probably comes at least in part from that l++ though. You can do this much more functionally, by first building your list of phia[[n]]+phia[[m]] using Table, calling Cos on that built list (Cos is Listable), then calling Total on that. – b3m2a1 Oct 16 '17 at 15:13

I do not agree with @b3m2a1: functional programming does not necessarily shine here, because things involving e.g., Subsets[phia,2] require $O(n^2)$ flops and $O(n^2)$ memory, while your implementation requires $O(n^2)$ flops but only $O(n)$ memory. So, in principle, your procedural code is fine (apart from the negligible fact that it does not seem to do what you expect); the issue is that the Mathematica kernel is not good at executing procedural code without preparation. That's what Compile is good for!

So, let's compile your code. I removed the ++l, because it is not necessary.

cf = Compile[{{x, _Real, 1}}, Block[{res = 0.},
Do[
Do[
If[m != n, res += Cos[x[[m]] + x[[n]]]],
{m, 1, Length[x]}],
{n, 1, Length[x]}];
res/(Length[x] (Length[x] - 1))
],
CompilationTarget -> "C"];


Now let's generate some data:

phia = RandomReal[{-1, 1}, 44359];


And let it run:

mean1 = cf[phia]; // AbsoluteTiming // First

(* 19.3747 *)


Still, cf requires $O(n^2)$ evaluations of trigonometric functions. We can exploit the identity $\cos(a+b) = \cos(a) \, \cos(b) - \sin(a) \, \sin(b)$ in order to cut that down to $O(n)$ trigonometric function evaluations as follows:

cf2 = Compile[{{cosx, _Real, 1}, {sinx, _Real, 1}},
Block[{res = 0., diagres = 0.},
Do[
Do[
res += (cosx[[m]] cosx[[n]] - sinx[[m]] sinx[[n]]),
{m, 1, Length[cosx]}],
{n, 1, Length[cosx]}
];
Do[diagres += (cosx[[m]]^2 - sinx[[m]]^2), {m, 1, Length[cosx]}];
(res - diagres)/(Length[cosx] (Length[cosx] - 1))
],
CompilationTarget -> "C"];

mean2 = cf2[Cos[phia], Sin[phia]]; // AbsoluteTiming// First
mean1 == mean2

(* 2.302 *)
(* True *)


Finally, a version that also employs parallelization; it includes some rounding error, but I think it's still okay.

cf3 = Compile[{{cosx, _Real}, {sinx, _Real}, {cosy, _Real, 1}, {siny, _Real, 1}},
cosx Total[cosy] - sinx Total[siny],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
f3 = x \[Function] With[{cosx = Cos[x], sinx = Sin[x]},
(Total[cf3[cosx, sinx, cosx, sinx]] - Total[cosx^2] + Total[sinx^2])/(Length[x] (Length[x] - 1))
];

mean3 = f3[phia]// AbsoluteTiming// First
Abs[mean2 - mean3]

(* 0.804182 *)
(* 4.05276*10^-12 *)


Running b3m2a1's implementation on this data set takes considerably longer, consumes over 50 GB(!) RAM and finally crashes the kernel. Thus, I give an example with a considerably smaller data set:

phia = RandomReal[{-1, 1}, 4435];
mean1 = cf[phia]; // AbsoluteTiming // First
mean2 = cf2[Cos[phia], Sin[phia]]; // AbsoluteTiming // First
mean3 = f3[phia]; // AbsoluteTiming // First
mean4 = Mean@Cos@Map[Total, Subsets[phia, {2}]]; // AbsoluteTiming // First

(* 0.188377 *)
(* 0.022097 *)
(* 0.008892 *)
(* 2.141 *)


Testing the accuracy:

mean2 - mean1
mean3 - mean1
mean4 - mean1

(* -2.66454*10^-15 *)
(* -7.14984*10^-14 *)
(* -5.88418*10^-14 *)


Fun fact

Due to the simplification by trigonometrics, we can obtain the same result with $O(n)$ complexity and in a completely vectorized way:

f5 = x \[Function] With[{cosx = Cos[x], sinx = Sin[x]},
(Total[cosx]^2 - Total[sinx]^2 - Total[cosx^2] +
Total[sinx^2])/(Length[cosx] (Length[cosx] - 1))
];
mean5 = f5[phia]; // AbsoluteTiming // First
mean1 - mean5

(* 0.000195 *)
(* -1.00253*10^-13 *)


In the large example above, f5[phia] takes 0.00089 seconds; this is several million times faster than the OP's implementation.

• You should compare the timing of your compiled version with b3m2a1's functional method (on the same machine). – bill s Oct 16 '17 at 16:36
• Good point. Done. – Henrik Schumacher Oct 16 '17 at 17:13
• wow. Thanks a lot Henry Schumacher, yes this is what i am doing . but when i run your code it gives me error CreateLibrary::cmperr: Compile error: LINK : fatal error LNK1104: cannot open file 'LIBCMT.lib' >> Compile::nogen: A library could not be generated from the compiled function. >> – IrfanS Oct 17 '17 at 5:01
• Hm. I googled the error message LNK1104 and it seems to be related to Visual Studio. This sort of tells me you are running Mathematica on Windows, right? Do you have a C compiler installed on your system? Could you please run the following code on Mathematica and send me the outcome? Needs["CCompilerDriver"]; CCompilers[] – Henrik Schumacher Oct 17 '17 at 5:45
• You need to enclose anything with backticks in double backticks: Needs["CCompilerDriver"]. – J. M. will be back soon Oct 25 '17 at 9:11

phia = RandomReal[{-1, 1}, 4985];