5
$\begingroup$

I tried to implement the fastest way of root finding of a system of nonlinear equations but still it is almost two times slower than the embedded FindRoot function in the Mathematica. In the following functions, f is the evaluation of functions and fs is the gradient of the function at the given initial points. The difference of the the following functions lies in the implementation of f and fs. In the first function, Experimental`CreateNumericalFunction were used to evaluate Jacobian of the function numerically. In the second function f was compiled and the third one inverse of the Jacobian function was first evaluated symbolically and then used later using injac function. I was wondering if someone can help to speed up the functions? possibly with NestWhile?

List of variables used in the functions

  • vars= List of variables
  • expr= List of equations to be solved
  • init= List of initial points
  • nit= Number of iterations
  • Prints= Printing the norms(norm[x_i+1-x_i] and norm[f[x_i+1]] ) if set to 1
  • tol = Tolerence

Function 1

Clear[solve];
solve[vars_List, expr_List, init_List, nit_Integer, Prints_Integer, 
  tol_Real] := 
 Block[{$MinPrecision = 16, $MaxPrecision = 16, ax, res, it = init, 
   f = Function[Evaluate[vars], #] & /@ expr, 
   fs = Experimental`CreateNumericalFunction[vars, 
     expr, {Length@vars}]},
  Unprotect[In, Out];
  Clear[In, Out];
  Protect[In, Out];
  ClearSystemCache[];
  Table[
   res = N[
     it - Inverse[
         fs["Jacobian"[#]] &@
          N[it, 16]].Transpose[{f[[#]] @@ N[it, 16] & /@ 
           Range[Length@vars]}] // Flatten, 16];
   If[Prints == 1,
    Print["norm[x_i+1-x_i]=", Norm[(it - res), 2]];
    Print["norm[f[x_i+1]]=", 
     Norm[(f[[#]] @@ N[res, 16] & /@ Range[Length@vars]), 2]], 
    Unevaluated[Sequence[]]];
   If[N[Norm[(it - res), 2], 16] < tol, Return["Exit", Table]];
   If[i == nit, Print["Solution has not converged"]];
   it = res, {i, nit}];
  res
  ]

Function 2

Clear[solve1];
solve1[vars_List, expr_List, init_List, nit_Integer, Prints_Integer, 
  tol_Real] := 
 Block[{$MinPrecision = 16, $MaxPrecision = 16, ax, res, it = init, 
   f = Compile[##] &[vars, #, RuntimeAttributes -> {Listable}, 
       "Parallelization" -> True] & /@ expr, 
   fs = Experimental`CreateNumericalFunction[vars, 
     expr, {Length@vars}]},
  Unprotect[In, Out];
  Clear[In, Out];
  Protect[In, Out];
  ClearSystemCache[];
  Table[
   res = N[
     it - Inverse[
         N[fs["Jacobian"[#]], 16] &@
          N[it, 16]].Transpose[{N[f[[#]], 16] @@ N[it, 16] & /@ 
           Range[Length@vars]}] // Flatten, 16];
   If[Prints == 1,
    Print["norm[x_i+1-x_i]=", Norm[(it - res), 2]];
    Print["norm[f[x_i+1]]=", 
     Norm[(f[[#]] @@ N[res, 16] & /@ Range[Length@vars]), 2]], 
    Unevaluated[Sequence[]]];
   If[N[Norm[(it - res), 2], 16] < tol, Return["Exit", Table]];
   If[i == nit, Print["Solution has not converged"]];
   it = res, {i, nit}];
  res
  ]

Function 3

injac[f_List?VectorQ, x_List] := 
     Inverse@Outer[D, f, x] /; Equal @@ (Dimensions /@ {f, x});
Clear[solve2];
solve2[vars_List, expr_List, init_List, inv_, nit_Integer, 
  Prints_Integer, tol_Real] := 
 Block[{$MinPrecision = 16, $MaxPrecision = 16, ax, res, it = init, 
   f = Function[Evaluate[vars], #] & /@ expr, fs = inv[expr, vars]},
  Unprotect[In, Out];
  Clear[In, Out];
  Protect[In, Out];
  ClearSystemCache[];
  Table[
   res = N[
     it - (fs /. 
          MapThread[
           Rule, {vars, 
            it}]).Transpose[{N[f[[#]], 16] @@ N[it, 16] & /@ 
           Range[Length@vars]}] // Flatten, 16];
   If[Prints == 1,
    Print["norm[x_i+1-x_i]=", Norm[(it - res), 2]];
    Print["norm[f[x_i+1]]=", 
     Norm[(f[[#]] @@ N[res, 16] & /@ Range[Length@vars]), 2]], 
    Unevaluated[Sequence[]]];
   If[N[Norm[(it - res), 2], 16] < tol, Return["Exit", Table]];
   If[i == nit, Print["Solution has not converged"]];
   it = res, {i, nit}];
  res
  ] 

Using the functions and results:

enter image description here

$\endgroup$
3
$\begingroup$

FindRoot does some fancy checking at each step. You could avoid that to save time. I doubt the basic Newton step can be done any faster than it is in FindRoot (however it is done). For an arbitrary n-dimensional solver, you would probably want to use LinearSolve than Inverse. Your input shows machine precision numbers, so I don't know if you really wanted to control the precision of the computation with $MinPrecision and $MaxPrecision. Those system parameters affect only arbitrary precision numbers but not machine reals or exact numbers. In other words, they have no effect on the examples shown in the OP.

Here's a stripped down solver that's almost twice as fast on a similar 2D example:

newtonStep[f_, vars_] := newtonStep[f, D[f, {vars}], vars];
newtonStep[f_, jac_, vars_][x_] := Block[vars,
   vars = x;
   x - LinearSolve[jac, f]];
FixedPoint[
  newtonStep[{x^2 + y^2 - 2, (x - 1)^2 + (y - 2)^2 - 5}, {x, y}],
  {1., -2.},
  100,
  SameTest -> Equal] // RepeatedTiming
(*  {0.00023, {1.4, -0.2}}  *)

FindRoot[{x^2 + y^2 - 2, (x - 1)^2 + (y - 2)^2 - 5}, {{x, 1.}, {y, -2.}}] // RepeatedTiming
(*  {0.00040, {x -> 1.4, y -> -0.2}}  *)

Note that there is some tolerance in Equal when comparing numbers. On numbers that differ only in their last seven bits, Equal will return True.

Now, if you want all of the robustness of FindRoot, then you have to code up some checks and so forth, and you find the gap in timing starts to close.

$\endgroup$
  • $\begingroup$ Computing the symbolic derivative only once in the beginning and creating a function fromt hat should also improve the speed a bit. $\endgroup$ – Henrik Schumacher Apr 7 '18 at 17:38
  • $\begingroup$ @HenrikSchumacher It does compute the derivative only once. Maybe having two forms of newtonStep is confusing? $\endgroup$ – Michael E2 Apr 7 '18 at 17:39
  • 1
    $\begingroup$ Ah, now I come to see it. The "double" construction and the Block-trick are indeed a bit confusing but very elegant! $\endgroup$ – Henrik Schumacher Apr 7 '18 at 17:52
  • 1
    $\begingroup$ @HenrikSchumacher It also depends on newtonStep evaluating to the version with the Jacobian before being passed to FixedPoint. If FixedPoint were HoldAll, it might not work efficiently. So one has to think carefully about all those evaluation rules. Every time I have such problem, I think, how neat you can do it this way in Mathematica; but every time I have to look at someone else's code, or my old code, I think, now what in the world is going on?!? $\endgroup$ – Michael E2 Apr 7 '18 at 17:55
  • 1
    $\begingroup$ Hm. I think the confusion is mostly caused by the naming of things. I would have found createNewtonIterator[f_, vars_] := With[{jac = D[f, {vars, 1}]}, Function[x, Block[vars, vars = x; x - LinearSolve[jac, f]]] ]; FixedPoint[ createNewtonIterator[{x^2 + y^2 - 2, (x - 1)^2 + (y - 2)^2 - 5}, {x, y}], {1., -2.}, 100, SameTest -> Equal] a bit more transparent. But that's a matter of taste, isn't it? (Although, the pure function created seems to be a tiny bit faster.) Anyway, I learnt something new! Many thanks! $\endgroup$ – Henrik Schumacher Apr 7 '18 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.