Reducing the computation time used by parallel table when finding a fixed point

Question

I am using Intel Core(™) processor, 32 bit operating system, 4GB RAM, Windows 7. I am running the following code and it is taking hours to obtain the result. Can any one help me in improving the codes so that it takes less time.

μ = 0.01214398779977472`

μ1 = 328823.8592475095`

Subscript[e, 1] = 0.0167`

A = 389.1779396462019`

ϕ = 0

Subscript[e, 2] = 0.0549`

Vx = 
  x + μ/(y^2 + (-1 + x + μ)^2)^(3/2) - x μ/(y^2 + (-1 + x + μ)^2)^(3/2) - 
  μ^2/(y^2 + (-1 + x + μ)^2)^(3/2) - x/(y^2 + (x + μ)^2)^(3/2) - 
  μ/(y^2 + (x + μ)^2)^(3/2) + x μ/(y^2 + (x + μ)^2)^(3/2) + 
  μ^2/(y^2 + (x + μ)^2)^(3/2) + x μ1/(2 A^3) + (3 x μ1 Cos[2 ϕ])/A^3 + 
  (3 y μ1 Sin[2 ϕ])/A^3 - (μ1 Cos[ϕ] Subscript[e, 1]^2)/(2 A^2); 

Vy =
  y - (y μ)/(y^2 + (-1 + x + μ)^2)^(3/2) - y/(y^2 + (x + μ)^2)^(3/2) + 
  (y*μ)/(y^2 + (x + μ)^2)^(3/2) + (y*μ1)/(2 A^3) - (3 y μ1 Cos[2 ϕ])/A^3 + 
  (μ1 Sin[ϕ])/A^2 + (3 x μ1 Sin[2 ϕ])/A^3 - (μ1 Sin[ϕ] Subscript[e, 1]^2)/(2 A^2);

Vxx = D[Vx, x];
Vyy = D[Vy, y];
Vxy = D[Vx, y];
Vyx = D[Vy, x];

newton[{x_, y_}] := 
  {x, y} + 
    {Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)], 
     Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]};

 sol = 
   {{ 0.188277, -0.9824026}, { 0.188277, 0.9824026}, { 1.152706, 0}, 
    {-0.998550, 0}, { 0.987856, 0}, { 0.835518, 0}, { -0.012144, 0}};

 tab = 
   ParallelTable[
     N[FixedPoint[newton, {i, j}, 10]], 
     {j, -2, 2, 0.003},
     {i, -2, 2, 0.003}]

 rules = Rule @@@ Transpose[{sol, Range[Length[sol]]}];

 newtab = Map[First @ Nearest[rules, #] &, tab, {2}];

 P1 = 
   ArrayPlot[newtab, 
     ColorFunction -> "Rainbow",  DataReversed -> True]

Answer 1

The principle issue with the code in the question is that the definition of newton does not replace x and y in Vx, etc. by the values of the arguments of newton. There are a number of ways to resolve this issue, including defining

newton[{x1_, y1_}] := {x1, y1} + {Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)], 
    Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]}/.{{x -> x1, y -> y1}}

With this resolved, we can focus on increasing the computational speed. Since almost all the time goes into repeatedly evaluating newton, our effort should go there. A factor of three or more is obtained by redefining newton as

newton[{x1_, y1_}] := Block[{}, {vx, vy, vxx, vyy, vxy} = {Vx, Vy, Vxx, Vyy, Vxy} /. 
  {x -> x1, y -> y1}; {x1, y1} + {-(vx vyy - vy vxy), (vx vxy - vy vxx)}/(vyy vxx - vxy^2)]

which takes advantage of the fact that Vxy == Vyx and also eliminates duplicated evaluations of intermediate results. Eliminating Simplify also saves a bit of time, because it is ineffective in the present calculation. (N in the definition of tab also can be eliminated, because it has no effect, but doing so saves negligible time.) With these changes, the computation completes in about eleven minutes on my four-processor PC, yielding

Perhaps, restructuring Vx, Vy, and their derivatives could save a bit more time.

Answer 2

There are several issues:

1. newton does not replace x and y. (See bbgrodfrey's answer.)

2. Every call to newton calls also a Simplify which is pretty slow. Simplifying it once at the time of definition should suffice.

You can resolve point 1 and 2 by defining:

 newton[{x_, y_}] := 
  Evaluate[{x, y} + {Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)], 
     Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]}];

Even that is probably still too slow. Compiling the code helps and it makes parallelization easier.

With[{xnew = N@newton[{x, y}][[1]], ynew = N@newton[{x, y}][[2]]}, 
  newtonize = Compile[{{x0, _Real, 1}, {maxiter, _Integer}},
     Block[{x, y, xbuffer},
      x = Compile`GetElement[x0, 1];
      y = Compile`GetElement[x0, 2];
      Do[
       xbuffer = xnew;
       y = ynew;
       x = xbuffer;
       ,
       {maxiter}];
      {x, y}],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"];
  ];

This CompiledFunction performs maxiter iterations on each input vector. Note that we use the option RuntimeAttributes -> {Listable} which creates a CompiledFunction that threads automatically over the matrix x0data of starting vectors below:

sol = Developer`ToPackedArray[N[{{0.188277, -0.9824026}, {0.188277, 0.9824026}, {1.152706, 0.}, {-0.998550, 0.}, {0.987856, 0.}, {0.835518, 0.}, {-0.012144, 0.}}]];
cols = Developer`ToPackedArray[
   List @@@ ColorData["Rainbow"] /@ Rescale[Range[Length[sol]]]];
colorize = Compile[{{a, _Integer}, {cols, _Real, 2}},
   Compile`GetElement[cols, a],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

δx = 0.003;
maxiter = 10;
tab = newtonize[
     Outer[List, Table[N[j], {j, -2, 2, δx}], Table[N[j], {j, -2, 2, δx}]], 
     maxiter
     ]; // AbsoluteTiming // First
newtab = ArrayReshape[
     Nearest[sol -> Automatic, Flatten[tab, 1]],
     Dimensions[tab][[1 ;; 2]]
     ]; // AbsoluteTiming // First
size = 600;

P1 = ImageRotate[
     ImageResize[Image[colorize[newtab, cols]], {size, size}]
     ]; // AbsoluteTiming // First
P1

0.578838

0.048854

0.183145

0.1231

Performing that with δx = 0.0005 and maxiter = 30 leads after 2 minutes to this picture:

Note also the following tweaks:

• We do no longer use a nested table for the construction of the initial points but Outer (among other things, this ensures that the generated array is a PackedArray).

• We do not Map a NearestFunction onto tab but rather apply Nearest with a syntax that enables internal low level parallelism.

• We ensure that sol is also a PackedArray so that Nearest can profit from this.

• We use Image in conjunction with a hand written method colorize instead of ArrayPlot in order to improve performance and stability (my FrontEnd had the urgent need to shut down several times).

Update

Corrected a bug in the CompiledFunction: The former version changed the value of x before the new value of y was computed. This was somewhat hidden by the way the code for x and y is injected into the compiled function. Buffering the new value for x (in the variable xbuffer) remedies this.