Re: FindRoot repeatedly evaluating function
- To: mathgroup at smc.vnet.net
- Subject: [mg121103] Re: FindRoot repeatedly evaluating function
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Sun, 28 Aug 2011 04:05:39 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201108251105.HAA24905@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
> 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++]
The first half of that DOES look clever, but the second has nothing to do
with it.
What did you leave out?
Bobby
On Sat, 27 Aug 2011 07:18:24 -0500, Oliver Ruebenkoenig
<ruebenko at wolfram.com> wrote:
>
>
> 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
>
--
DrMajorBob at yahoo.com
- References:
- FindRoot repeatedly evaluating function
- From: Simon Pearce <Simon.Pearce@nottingham.ac.uk>
- FindRoot repeatedly evaluating function