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

Question

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.

Answer 1

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.

Answer 2

So I left a comment about this, but here's where functional programming shines. Consider this:

phia = RandomReal[{-1, 1}, 4985];
Mean@Cos@Map[Total, Subsets[phia, {2}]] // AbsoluteTiming

{3.20715, 0.70717}

If I understand correctly, that's all you want out of your code. And it's short and fast.