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>

**Re: FindRoot repeatedly evaluating function**

**Re: FindRoot repeatedly evaluating function**

**Re: FindRoot repeatedly evaluating function**

**Re: FindRoot repeatedly evaluating function**