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: Sun, 28 Aug 2011 04:05:50 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201108251105.HAA24905@smc.vnet.net> <201108271218.IAA17810@smc.vnet.net>
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 >> > > >
- References:
- FindRoot repeatedly evaluating function
- From: Simon Pearce <Simon.Pearce@nottingham.ac.uk>
- Re: FindRoot repeatedly evaluating function
- From: Oliver Ruebenkoenig <ruebenko@wolfram.com>
- FindRoot repeatedly evaluating function