       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>

```> A colleague of mine came up with this reformulation. This is very clever
>
> 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==1,x'==a},x,{S,0,10}][];
>>  Print[{a,x/.sol}]; x/.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 == 1, x' == a},
>       x, {S, 0, 10}][];
>    Print[InputForm[{a, x /. sol}]]; x /. 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 == 1, x' == 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 == 1,
>       x' == 0}, {}, {S, 10, 10}];
>
> ClearAll[ff]
> With[{state = nds},
>   ff[a_?NumericQ] := (Module[{sol, tmp, tnds},
>      tnds = First@
>        NDSolve`Reinitialize[state, {x == 1, x' == 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 == 1, x' == 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
>
> 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

```

• Prev by Date: Re: FindRoot repeatedly evaluating function
• Next by Date: Re: FindRoot repeatedly evaluating function
• Previous by thread: Re: FindRoot repeatedly evaluating function
• Next by thread: Re: FindRoot repeatedly evaluating function