1
$\begingroup$

The following code is my attempt of employing Broyden's method for root finding on the function $f(x,y)=(e^{xy}-y^2-2,\cos(x+y)+\frac{1}{2})$. Where the first matrix is the Jacobian, then it gets updated using an estimate. So $k[l,k]$ computes the first 12 iterations. However, this program is extremely slow and I am not sure why. Could anyone point out why it is so slow and how I could improve it?

k[l_, k_] := Module[{m}, 
  f[{x_, y_}] := {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  h = {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  b = {x, y};
  
  m[0] := {{l}, {k}};
  
  J[0] := D[h, {b}];
  m[n_] := 
   m[n] = N[m[n - 1]] - 
      Inverse[N[J[n - 1]]].N[f[m[n - 1]]] /. {x -> 
       Part[Part[N[m[n - 1]], 1], 1], 
      y -> Part[Part[N[m[n - 1]], 2], 1]};
  J[n_] := 
   J[n] = N[
      J[n - 1]] + ((N[f[m[n]]] - N[f[m[n - 1]]] - 
          N[J[n - 1]].(N[m[n] - m[n - 1]]))/
        Norm[N[m[n] - m[n - 1]]]^2).Transpose[N[m[n]] - N[m[n - 1]]];
  Table[m[a], {a, 1, 12}]]

I would appreciate any input.

EDIT: I found some optimization to my code. I removed the redundant $N$ as per a suggestion in the comments. The optimization I found is to deal with numerical matrices form the get go. The original code did everything for arbitrary $x,y$

k[l_, k_] := Module[{m}, 
  f[{x_, y_}] := {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  h = {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  b = {x, y};
  
  m[0] := N[{{l}, {k}}];
  
  J[0] := 
   N[D[h, {b}]] /. {x -> Part[Part[m[0], 1], 1], 
     y -> Part[Part[m[0], 2], 1]};
  m[n_] := m[n] = m[n - 1] - LinearSolve[J[n - 1], f[m[n - 1]]];
  J[n_] := 
   J[n] = J[
      n - 1] + ((f[m[n]] - f[m[n - 1]] - J[n - 1].(m[n] - m[n - 1]))/
        Norm[m[n] - m[n - 1]]^2).Transpose[m[n] - m[n - 1]];
  Table[m[a], {a, 1, 15}]]
```
$\endgroup$
9
  • $\begingroup$ Your code is very hard to follow because of all the brackets, parts, unnecessary Ns, and in-body function definitions. An iterative algorithm might be better with Nest. I'm not familiar with this numerical method but I've rewritten it here: pastebin.com/Xhu4n2v8 . It's quite sensitive to the initial value and converges a bit slowly for this function, so I may not have got it right, but you get the gist of it. $\endgroup$
    – flinty
    Dec 8, 2020 at 17:42
  • 1
    $\begingroup$ @flinty Since you went to the trouble of writing code and posting it, please post it AS AN ANSWER so it gets preserved with its attached question and everybody can see it. $\endgroup$
    – MarcoB
    Dec 8, 2020 at 17:52
  • $\begingroup$ @MarcoB the question is about improving the existing code, not rewriting it, so I didn't think it constituted an answer as to why OPs code is slow. $\endgroup$
    – flinty
    Dec 8, 2020 at 18:08
  • 1
    $\begingroup$ The idea of Broyden's method is to invert/factorize the initial Jacobian only once in the beginning (e.g. with LinearSolve) and to use the Sherman-Morrison formula to deal with the rank-1-updates. At the moment OP's code looks as if the a linear solve were performed in each iteration. (The many linear solves are typically the reason why Newton's method takes so long. Not in dimension, 2 though.) $\endgroup$ Dec 8, 2020 at 19:07
  • 2
    $\begingroup$ I am just saying that you did not implement Broyden's method. $\endgroup$ Dec 8, 2020 at 19:39

1 Answer 1

3
$\begingroup$

The problem with your implementation is the fact that the values of the Jacobian are kept symbolic at every iteration, so their complexity grows exponentially.

For each iteration J[n] should have the value of {x, y} updated from m[n-1]. I am also unravelling the internal definitions of all those functions, removing the unnecessary N calls, and somewhat simplifying the rest of the code:

ClearAll[f, k, J, m]
f[{x_, y_}] = {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
J[0] = D[f[{x, y}], {{x, y}}];

(* note the replacement of the values of x and y at the end *)
(* that is critical to obtain numerical values for J at each step *)
(* that was the source of your slowdown *)
J[n_] := J[
    n] = (J[n - 
        1] + ((f@m[n] - f@m[n - 1] - J[n - 1].(m[n] - m[n - 1]))/
         Norm[m[n] - m[n - 1]]^2).(m[n] - m[n - 1])) /. 
    Thread[{x, y} -> m[n - 1]];


m[n_] := m[n] = 
  m[n - 1] - (Inverse[J[n - 1]].f@m[n - 1] /. 
     Thread[{x, y} -> m[n - 1]])

Let's see a trial run:

m[0] = N@{1, 4};
results = Table[m[n], {n, 0, 15}]
ListPlot@Transpose@results

{{1., 4.}, {1.00863, 3.17414}, {0.934639, 3.25276}, {0.886161, 3.30122}, {0.852723, 3.33465}, {0.829105, 3.35826}, {0.812153, 3.3752}, {0.799845, 3.38751}, {0.790838, 3.39651}, {0.784208, 3.40314}, {0.779307, 3.40804}, {0.775673, 3.41167}, {0.772972, 3.41437}, {0.770961, 3.41638}, {0.769463, 3.41787}, {0.768345, 3.41899}}

convergence plot


The problem still remains that this approach based on matrix inversion at every step is not very efficient, and I am also not sure that it is in the spirit of the method proposed (although I am not familiar with it).

Finally,

FindRoot[f[{x, y}], {x, 1}, {y, 4}]

(* Out: {x -> 0.764945, y -> 3.42385} *)

is built-in and orders of magnitude faster, of course.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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