Re: FindRoot repeatedly evaluating function

```On Sat, 27 Aug 2011, DrMajorBob wrote:

>> 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[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
>>
>> 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]};
>>

Apparently it was not entirely obvious what to do here so here is the
code:

ClearAll[f]
With[{eqns = seqns, vars = svars}, Clear[f, J, ysol];
f[y_?NumericQ] := Part[ysol[y], 1];
J[y_?NumericQ] := (Print["Jackpot"]; {{Part[ysol[y], 2]}});
ysol[y_?NumericQ] := Module[{sol, yvars}, yvars = vars /. p -> y;
sol = First[NDSolve[eqns /. p -> y, yvars, {S, 10, 10}]];
sol = (yvars /. sol) /. S -> 10;
Print[InputForm[{y, sol}]];
sol]]

;-)

Oliver

>> 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: decoding inbuilt function
• Previous by thread: Re: FindRoot repeatedly evaluating function
• Next by thread: Re: FindRoot repeatedly evaluating function