Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2011

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

Search the Archive

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

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