Parallelize for loop in mathematica

Question

I have the following code which essentially consists of a first for loop that iterates through a very large number of rows and a small for loop that iterates through a small number of columns. The algorithm computes a sum and I wish to parallelize the outer for loop to obtain a better performance. Every time the outer for loop gets executed a new function gets defined and then inside the inner for loop new terms get added to the sum. I have little experience with parallelization toolkits in Mathematica and wonder if anybody has suggestions to increase the performance of my algorithm? It seems that parallelsum might not be right thing for me, as the summand gets changed in every iteration?

 Matrix = Import["/tilde_residues.txt", "CSV"];  
 Dims := Dimensions[Ytilde]; Rows := Dims[[1]]; Columns := Dims[[2]] - 1; 
 l := Columns/4; 
 w := E^((2*Pi*I)/3); 
 sum = 0; 
For[i = 1, i <= Rows, i++, m := Matrix[[i, 4*l + 1]]; 
  f[x_] := 1/Product[(-(w*Matrix[[i, k + l]]) + w^2*Matrix[[i, k]] + x)*
          (w^2*Matrix[[i, k + 2*l]] - w*Matrix[[i, k + 3*l]] + 
    x), {k, 1, l}]; 
   For[j = 1, j <= l, j++, 
    sum = Residue[f[x], 
                 {x, -(w^2*Matrix[[i, j]] - 
           w*Matrix[[i, j + l]])}] + sum]; ]; 
     Print[Simplify[sum]];

Answer 1

It is tempting to simply parallelize your approach. One seemingly obvious way would be to replace For and Sum loops by Do loops (For and Sum loops in Mathematica are super inefficient). And then to replace the outer Do loop by a ParallelDo. However, since quite a lot data sharing is going on, speeding up the code this way is not as easy as it sounds.

Instead I try to address the main bottleneck of your code: The Residue function. It has to perform some symbolic computation and this is always slower than pure machine number crunching. In fact, the residue in your case can be computed exactly if we assume that the data in each row of Matrix is sufficiently generic so that poles can only occur for $j = k$.

Moreover I assume about Matrix that it is a packed array of double precision complex numbers. In particular, no symbols are involved. Here a random matrix that I use for test purposes.

SeedRandom[123];
Matrix = RandomComplex[{-1. - 1. I, 1. + 1. I}, {100, 20}];

This is your code:

RepeatedTiming[
 Module[{Dims, l, w, sum, x, f},
  Dims := Dimensions[Matrix]; Rows := Dims[[1]]; 
  Columns := Dims[[2]] - 1;
  l := Quotient[Columns, 4];
  w := E^((2*Pi*I)/3);
  sum = 0;
  For[i = 1, i <= Rows, i++,
   m := Matrix[[i, 4*l + 1]];
   f[x_] := 1/Product[(-(w*Matrix[[i, k + l]]) + w^2*Matrix[[i, k]] + x)*(w^2*Matrix[[i, k + 2*l]] - w*Matrix[[i, k + 3*l]] + x), {k, 1, l}];
   For[j = 1, j <= l, j++,
    sum = Residue[f[x], {x, -(w^2*Matrix[[i, j]] - w*Matrix[[i, j + l]])}] + sum];
   ];
  sum
  ]
 ]

{0.267779, -144.106 + 98.4052 I}

This is what I get if I apply some simplifications and compute the residues by hand:

RepeatedTiming[
 Module[{rows, cols, l, w, Y, Z, sum2, y, z},
  {rows, cols} = Dimensions[Matrix];
  l = Quotient[cols - 1, 4];
  w = E^((2. Pi I)/3);
  
  Y = w^2 Matrix[[All, 1 ;; l]] - w Matrix[[All, 1 + l ;; 2 l]];
  Z = w^2 Matrix[[All, 2 l + 1 ;; 3 l]] - w Matrix[[All, 3 l + 1 ;; 4 l]];
  
  sum2 = 0. I;
  
  Do[
   y = Y[[i]];
   z = Z[[i]];
   Do[
    sum2 += 1./Times[
       Product[(y[[k]] - y[[j]]) (z[[k]] - y[[j]]), {k, 1, j - 1}],
       (z[[j]] - y[[j]]),
       Product[(y[[k]] - y[[j]]) (z[[k]] - y[[j]]), {k, j + 1, l}]
       ]
    , {j, 1, l}];
   , {i, 1, rows}];
  sum2
  ]
 ]

{0.0033698, -144.106 + 98.4052 I}

This is already 75 times faster - without any parallelization involved!

A somewhat cleaner version using Sum and Product along with Boole[i==j] looks like this:

RepeatedTiming[
 Module[{rows, cols, l, w, Y, Z, sum2, y, z, yj},
  {rows, cols} = Dimensions[Matrix];
  l = Quotient[cols - 1, 4];
  w = E^((2. Pi I)/3);
  Y = w^2 Matrix[[All, 1 ;; l]] - w Matrix[[All, 1 + l ;; 2 l]];
  Z = w^2 Matrix[[All, 2 l + 1 ;; 3 l]] - 
    w Matrix[[All, 3 l + 1 ;; 4 l]];
  sum2 = 0. I;
  Sum[
   y = Y[[i]];
   z = Z[[i]];
   Sum[
    yj = y[[j]];
    1./Product[(y[[k]] - yj + Boole[j == k]) (z[[k]] - yj), {k, 1, l}]
    , {j, 1, l}]
   , {i, 1, rows}]
  ]
 ]

Next, I compile the two do loops into a Listable CompiledFunction, that can be mapped parallely in a list of vectors (a.k.a. a matrix):

cRowSum = Compile[{{y, _Complex, 1}, {z, _Complex, 1}},
   Block[{prod, l, yj, yk, zk, sum},
    l = Min[Length[y], Length[z]];
    sum = 0. I;
    Do[
     prod = 1. + 0. I;
     yj = Compile`GetElement[y, j];
     Do[
      yk = Compile`GetElement[y, k];
      zk = Compile`GetElement[z, k];
      prod *= (yk - yj + Boole[j == k]) (zk - yj);
      , {k, 1, l}];
     sum += (1. + 0. I)/prod;
     , {j, 1, l}];
    
    sum
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Then we can get the result simply by evaluating:

RepeatedTiming[
 sum3 = Total[cRowSum[
    w^2 Matrix[[All, 1 ;; l]] - w Matrix[[All, 1 + l ;; 2 l]],
    w^2 Matrix[[All, 2 l + 1 ;; 3 l]] - w Matrix[[All, 3 l + 1 ;; 4 l]]
    ]]
 ]

{0.0000526625, -144.106 + 98.4052 I}

This says it takes just 53 microseconds. (But we have to be aware that timing statistics for such small timings are not very precise - you get the idea anyways.) This is about 5000 times faster than the original code!

Btw.: Parallelization is counterproductive for such short runtimes. The overhead of parallelization is typically too large. When I turn off parallelization with Parallelization -> False, then I get twice as fast (0.0000273632 sec) on my machine.