# Implementation of the Newton-Raphson method to solve a system of nonlinear equations(Speed comparison)

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, ExperimentalCreateNumericalFunction 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 = ExperimentalCreateNumericalFunction[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 = ExperimentalCreateNumericalFunction[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 /.
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: 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.

• Computing the symbolic derivative only once in the beginning and creating a function fromt hat should also improve the speed a bit. Apr 7 '18 at 17:38
• @HenrikSchumacher It does compute the derivative only once. Maybe having two forms of newtonStep is confusing? Apr 7 '18 at 17:39
• Ah, now I come to see it. The "double" construction and the Block-trick are indeed a bit confusing but very elegant! Apr 7 '18 at 17:52
• @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?!? Apr 7 '18 at 17:55
• 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! Apr 7 '18 at 18:11