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
- Follow-Ups:
- Re: FindRoot repeatedly evaluating function
- From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
- Re: FindRoot repeatedly evaluating function
- References:
- FindRoot repeatedly evaluating function
- From: Simon Pearce <Simon.Pearce@nottingham.ac.uk>
- FindRoot repeatedly evaluating function