MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: FindRoot repeatedly evaluating function

  • To: mathgroup at smc.vnet.net
  • Subject: [mg121104] Re: FindRoot repeatedly evaluating function
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Sat, 27 Aug 2011 08:18:24 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <201108251105.HAA24905@smc.vnet.net>


Hello Simon,

On Thu, 25 Aug 2011, Simon Pearce wrote:

> Hi Mathgroup,
>
> When I use FindRoot[f[y],y] I find that the inner function f is evaluated 3 or 4 times at each value of y (or at least very similar values), even if y is far from the root. This has obvious implications to the speed of my code.
> Can anyone explain why this is the case, and tell me any way to stop it from repeatedly evaluating f? If I use f[a]:=f[a]=... then it uses the stored result, but I don't want to store thousands of such real valued expressions.
>
> The following simple code shows the essence of the problem, using Print to show where the function is evaluated and its value there.
>
> f[a_?NumericQ]:=Module[{sol},
>  sol=NDSolve[{x''[S]-x'[S]+x[S]==0,x[0]==1,x'[0]==a},x,{S,0,10}][[1]];
>  Print[{a,x[10]/.sol}]; x[10]/.sol ]
> FindRoot[f[y],{y,6}]
>
> Thanks,
> Simon Pearce
> Postdoctoral Researcher
> The Centre for Plant Integrative Biology
> School of Biosciences
> University of Nottingham
>
>
> This message and any attachment are intended solely for the addressee and may contain confidential information. If you have received this message in error, please send it back to me, and immediately delete it.   Please do not use, copy or disclose the information contained in this message or in any attachment.  Any views or opinions expressed by the author of this email do not necessarily reflect the views of the University of Nottingham.  This message has been checked for viruses but the contents of an attachment may still contain software viruses which could damage your computer system: you are advised to perform your own checks. Email communications with the University of Nottingham may be monitored as permitted by UK legislation.
>

here are some ideas you may find useful.

Here is a slight modification of your example:

f[a_?NumericQ] :=
  Module[{sol},
   sol = NDSolve[{x''[S] - x'[S] + x[S] == 0, x[0] == 1, x'[0] == a},
      x, {S, 0, 10}][[1]];
   Print[InputForm[{a, x[10] /. sol}]]; x[10] /. sol]

e = s = 0;
FindRoot[f[y], {y, 6}
  , StepMonitor :> s++, EvaluationMonitor :> e++
  ]

on my machine this gives e==17 and s==5.

Sometimes the Jacobian that needs to be computed does not need to be 
updated every step.

s = e = 0;
FindRoot[f[y], {y, 6}, Method -> {"Newton", "UpdateJacobian" -> 2}
  , StepMonitor :> s++, EvaluationMonitor :> e++]

This gives me e==11 and s==4.

Another thing to try is to use a finite difference approximation to the 
Jacobian that is of lover order.

s = e = 0;
FindRoot[f[y], {y, 6}
  , Jacobian -> {"FiniteDifference", "DifferenceOrder" -> 2}
  , StepMonitor :> s++, EvaluationMonitor :> e++]

this gives me e==10 and s==3.

Combining the two last approaches:

s = e = 0;
FindRoot[f[y], {y, 6}
  , Method -> {"Newton", "UpdateJacobian" -> 2}
  , Jacobian -> {"FiniteDifference", "DifferenceOrder" -> 2}
  , StepMonitor :> s++, EvaluationMonitor :> e++]

e==9 and s==3.

Since you seem concerned about performance in both speed and memory, here
is another idea.

First you cloud call NDSolve like this

if = First[
   x /. NDSolve[{x''[S] - x'[S] + x[S] == 0, x[0] == 1, x'[0] == 1},
     x, {S, 10, 10}]]

{S,10,10} stores just the last value in the interpolation function.

if["ValuesOnGrid"]

Compare that to {S,0,10}.

Next, you could separate the equation processing from the solution like 
so:

nds = First@
    NDSolve`ProcessEquations[{x''[S] - x'[S] + x[S] == 0, x[0] == 1,
      x'[0] == 0}, {}, {S, 10, 10}];

ClearAll[ff]
With[{state = nds},
  ff[a_?NumericQ] := (Module[{sol, tmp, tnds},
     tnds = First@
       NDSolve`Reinitialize[state, {x[0] == 1, x'[0] == a}];
     NDSolve`Iterate[tnds];
     sol = NDSolve`ProcessSolutions[tnds, 1];
     tmp = x[10.] /. sol;
     Print[InputForm[{a, tmp}]]; tmp])]

s = 0; e = 0
FindRoot[ff[y], {y, 6}
  , StepMonitor :> s++, EvaluationMonitor :> e++]

gives e==17 and s==5

If you allow caching (May I ask why you do not?)

ClearAll[ff]
With[{state = nds},
  ff[a_?NumericQ] := (ff[a] =
     Module[{sol, tmp, tnds},
      tnds = First@NDSolve`Reinitialize[state, {x[0] == 1, x'[0] == a}];
      NDSolve`Iterate[tnds];
      (*Print["sv: ",tnds@"SolutionVector"];*)
      sol = NDSolve`ProcessSolutions[tnds, 1];
      tmp = x[10.] /. sol;
      Print[InputForm[{a, tmp}]]; tmp])]

s = 0; e = 0
FindRoot[ff[y], {y, 6}
  , StepMonitor :> s++, EvaluationMonitor :> e++]

s and e stay the same but as can be seen from the print statement those 
values are fetched. To be quite honest I have not looked at the code and 
do not know how the algorithm works exactly.

Possibly another speed improvement is possible if you do not construct the 
interpolation function but use the value from the tnds@"SolutionVector" 
directly.

A colleague of mine came up with this reformulation. This is very clever 
since then you have access to the Jacobian.

eqns = {Derivative[2, 0][x][S, p] - Derivative[1, 0][x][S, p] +
      x[S, p] == 0, x[0, p] == 1, Derivative[1, 0][x][0, p] == p};
seqns = {eqns, D[eqns, p]};
svars = {x[S, p], Derivative[0, 1][x][S, p]};

s = 0; e = 0;
FindRoot[f[y], {y, 6}, Jacobian -> J[y]
  , StepMonitor :> s++, EvaluationMonitor :> e++]

gives

s==3 and e==4

So, if you have the Jacobian it is good to use it. Even if you do not have 
the Jacobian you can sometimes say something about it's structure and then 
use a SparseArray to convey that to FindRoot.

Here are some links to the documentation

(* FindRoot *)
(*tutorial/UnconstrainedOptimizationNewtonsMethodRoot*)
(*tutorial/UnconstrainedOptimizationIntroductionNonlinearEquations*)
(*tutorial/UnconstrainedOptimizationSpecifyingDerivatives*)


(* NDSolve`XY *)
(* tutorial/NDSolveStateData *)

Hope this helps.

Oliver




  • Prev by Date: Re: FindRoot repeatedly evaluating function
  • Next by Date: Re: Problem with NIntegrating a compiled function
  • Previous by thread: Re: FindRoot repeatedly evaluating function
  • Next by thread: Re: FindRoot repeatedly evaluating function