Shadowing Chaotic Motion

Stan Wagon

Macalester College, St. Paul, Minnesota

Rob Knapp

Wolfram Research, Inc., Champaign, Illinois

How the optimization works:

To see how FindMinimum is working to get the best value, let's dissect what we did for the logistic map:

Shadow[noisy_, n_, start_] := (
residual[x0_?NumericQ] :=Block[{$MinPrecision = n,$MaxPrecision = n},
fTraj[x0,n]-Take[noisy,n+1]];
FindMinimum[0,{x0,SetPrecision[start, n]}, WorkingPrecision->n,
Method->{"LevenbergMarquardt","Residual"->residual[x0],"Jacobian"->sensitivity[x0,n]}])

Note that the function to FindMinimum is just given as 0 here, which, of course would not be terribly interesting to optimize.  When Method->"LevenbergMarquardt" is used with a specified residual, r, the objective function is taken to be (1/2)r·r.

Some definitions

Here are some definitions from before just to make sure everything works as I expect:

r=4; fQuad[x_]:=r  x (1-x);
der[x_]:=Evaluate[fQuad'[x]]
fIter[x_,n_]:=Nest[fQuad,x,n];
fTraj[x_,n_]:=NestList[fQuad,x,n];
sensitivity[x_,n_]:=Most[FoldList[Times,{1},der[fTraj[x,n]]]];

init = 1/10;
noisy=SetPrecision[fTraj[N[init],100],100];

SetAttributes[FixedPrecisionEvaluate, HoldRest];
FixedPrecisionEvaluate[p_, toeval_] := Block[{$MinPrecision = p, $MaxPrecision = p}, toeval];

Extent[points_] := Module[{xpts, ypts}, 
{xpts, ypts} = Transpose[points];
 {{Min[xpts],Max[xpts]},{Min[ypts],Max[ypts]}}];

PlotExtent[plot_]:= Extent[Apply[Join,Cases[plot, Line[pts_]:>pts,Infinity]]];

qminpt[qp_, x_] := Module[{cl = CoefficientList[qp,x],m},m=-cl[[2]]/(2 cl[[3]]);{m, qp /. x->m}]

The function

Define the residual as a function of (Subscript[x, 0]) for a given segment of the noisy trajectory (computed to the precision of the noisy trajectory):

rfun[ntraj_List][x0_?NumberQ] := Module[{p = Precision[ntraj], n = Length[ntraj] - 1},
FixedPrecisionEvaluate[p, fTraj[SetPrecision[x0, p], n] - ntraj]];

and f = (1/2)r.r

fobj[ntraj_List][x0_?NumberQ]:=Module[{res=rfun[ntraj][x0],p=Precision[ntraj]},FixedPrecisionEvaluate[p,((res.res)/2)]]

It is informative to look at the function for different lengths n:

Manipulate[Module[{ϵ = (10^le), p = Precision[noisy]},
Plot[fobj[Take[noisy, n]][init + δ],{δ, -ϵ,ϵ}, WorkingPrecision->p, PlotRange->All]],{n,2,Length[noisy],1},{{le,0,"(Subscript[log, 10])(ϵ)"},-20,0,1}]

Quadratic Model

Most of the methods that FIndMinimum uses develop a quadratic model at a given point and determine the next step by the minimizer of that quadratic model.

The quadratic model for the LevenbergMarquardt (GaussNewton) method is formed by approximating the Hessian matrix using (J^T)J where J is the Jacobian matrix ((∂r)/(∂x)):

f = (1/2)(r^T).r
((∂f)/(∂x))=(1/2)(((∂(r^T))/(∂x)).r+(r^T)((∂r)/(∂x)))=(1/2)((J^T)r+(r^T)J) = (J^T)r
(((∂^2)f)/(∂(x^2)))= ((∂(J^T))/(∂x)).r+(J^T)((∂r)/(∂x))=((∂(J^T))/(∂x)).r+(J^T)J 

if the residual is small near the minimum, generally the second derivative term ((∂(J^T))/(∂x))·r is negligible enough so that you still get the same quadratic convergence for smooth functions as you do with Newton's method.

First define the Jacobian:

Jfun[ntraj_][x0_?NumberQ] :=Module[{p = Precision[ntraj], n = Length[ntraj] - 1},
FixedPrecisionEvaluate[p, sensitivity[SetPrecision[x0, p], n]]]

This defines a the quadratic model function for a given point (Subscript[x, 0]):

Clear[qmodel];
qmodel[ntraj_][x0_] := Function[Evaluate[Module[{p = Precision[ntraj],x0p, res, J, Jt},
x0p = SetPrecision[x0, p];
res = rfun[ntraj][x0];
 J = Jfun[ntraj][x0];
Jt = Transpose[J];
FixedPrecisionEvaluate[p, Part[res.res/2 +Jt.res (#-x0) + Jt.J (#-x0)^2/2,1,1]]]]]

qmodel[Take[noisy, 50]][init]

Manipulate[
Module[{
ntraj = Take[noisy, n], p = Precision[noisy], x0p, f0, pf, prange},
x0p = SetPrecision[x0, p];
f0 = fobj[ntraj][x0p];
pf = Plot[fobj[ntraj][x],{x,0,1}, PlotRange->All, WorkingPrecision->p];
prange = PlotExtent[pf];
Show[
Plot[Evaluate[qmodel[ntraj][x0p][x]],{x,0,1}, PlotStyle->Red, Epilog->{Red,PointSize[0.025],Point[{x0p, f0}]}, WorkingPrecision->p],
pf, PlotRange->prange]],
{n, 2, 5,1},{{x0,init},0,1}]

Manipulate[
Module[{ϵ = (10^le),
ntraj = Take[noisy, n], p = Precision[noisy], x0p, f0, pf, prange},
x0p = SetPrecision[init+SetPrecision[x0r,p]*ϵ, p];
f0 = fobj[ntraj][x0p];
pf = Plot[fobj[ntraj][init+δ],{δ,-ϵ,ϵ}, PlotRange->All, WorkingPrecision->p];
prange = PlotExtent[pf];
Plot[{qmodel[ntraj][x0p][init+δ],fobj[ntraj][init+δ]},{δ,-ϵ,ϵ}, PlotStyle->{Red,Blue}, Epilog->{Red,PointSize[0.025],Point[{x0p-init, f0}]}, WorkingPrecision->p, PlotRange->prange]],
{n, 2,70,1},{{x0r,0},-1,1},{{le,0,"(Subscript[log, 10])(ϵ)"},-20,0,1}]

Sequence of Steps

LMSteps[noisy_, start_]:=
Module[{p = Precision[noisy],c = 0, g = 0, s= 0},Quiet[Join[{start}, Reap[FindMinimum[fobj[noisy][x],{{x,SetPrecision[start, p]}},Method->{"LevenbergMarquardt", "Residual"->rfun[noisy][x], "Jacobian"->{Jfun[noisy][x], EvaluationMonitor:>g++}}, WorkingPrecision->p, StepMonitor:>(s++;Sow[x]), EvaluationMonitor:> c++];If[TrueQ[stats],Print[c, " residual evaluations and ",g, " Jacobian evaluations and ",s," steps"]]][[2,1]]]]]

Manipulate[
Module[{
ntraj = Take[noisy, n], p = Precision[noisy], x0p, f0, pf, prange,steps, colors,x,models, minpoints, tangentpoints, showpoints,xmin, xmax},
x0p = SetPrecision[x0, p];
f0 = fobj[ntraj][x0p];
steps = LMSteps[ntraj, x0p];
models = Map[qmodel[ntraj][#][x]&, steps];
minpoints = Map[qminpt[#,x]&, models];
tangentpoints = Map[{#, fobj[ntraj][#]}&, steps];
{xmin, xmax}= First[Extent[Join[minpoints, tangentpoints]]];
xmin = Min[0, xmin - .1];
xmax = Max[xmax+.1, 1];
colors = ColorData["Rainbow"];
colors =Table[colors[c], {c,0,1,1./(Length[steps] - 1)}];
pf = Plot[fobj[ntraj][x],{x,xmin, xmax}, PlotRange->All, WorkingPrecision->p, PlotStyle->{Thickness[0.01], Green}];
prange = Extent[Join[Transpose[PlotExtent[pf]], minpoints, tangentpoints]];
showpoints = MapThread[{#1, {PointSize[0.015],Point[#2]}, {PointSize[0.025],Point[#3]}}&, {colors,  minpoints, tangentpoints}];
Show[
Plot[Evaluate[models],{x,xmin, xmax}, PlotStyle->colors, WorkingPrecision->p, Epilog->showpoints],
pf, PlotRange->prange]],
{n, 2, 5,1},{{x0,init},0,1}]

Step restriction

As you saw above, the quadratic model is not just blindly used -- there is some sort of step restriction.

A line search considers values along the step direction and ensures that the new value has sufficiently smaller function value than the prevous one -- if not, the step is reduced until this is satisfied.

In a trust region method, a region where the quadratic model is "trusted" be to a reasonably good model for the function is used to limit the step size. If the model function value and the actual function value are too different, the region size is reduced. The "LevenbergMarquardt" method uses a trust region approach.

Why not another method?

Why did Stan use the LevenbergMarquardt method?

For one, the default method is a disaster here:

FMSteps[method_][ntraj_, start_] :=Block[{p = Precision[noisy],c=0,s=0,g = 0},Join[{start}, Reap[FindMinimum[fobj[ntraj][x],{{x,SetPrecision[init, p]}},Gradient:> {Dot[Transpose[Jfun[ntraj][x]], rfun[ntraj][x]],EvaluationMonitor:>g++},Method->method, WorkingPrecision->p, StepMonitor:>(s++;Sow[x]), EvaluationMonitor:>c++];If[TrueQ[stats],Print[c, " function evaluations and ",g, " gradient evaluations and ",s," steps"]]][[2,1]]]]

stats = True;
n50 = Take[noisy, 50];

FMSteps[Automatic][n50, init]

The method runs into trouble because the sensitivity is great -- the default method builds up an approximation to the Hessian matrix, but at the beginning, it is not good enough to give a reasonable step.  If you decrease n, it becomes less sensitive and works to some extent.

FMSteps[Automatic][Take[noisy, 20], init]

Newton's method works quite well:

FMSteps["Newton"][n50, init]

FMSteps[{"Newton", "StepControl"->"TrustRegion"}][n50, init]

LMSteps[n50, init]

The results are very close:

Last[%] - Last[%%]

The big difference is in the number of gradient evaluations. The reason for this is that Newton's method requires a Hessian matrix. When it is not supplied and cannot be computed symbolically as is the case here, FindMinimum uses finite differences (of the gradient) to approximate the Hessian.

Now recall that the LevenbergMarquardt method also uses an approximate Hessian, but in this case, it comes from the Jacobian and no second derivatives need to be approximated directly.

Summary

There are a few points I hope that I have illustrated:

• To get a good global result on a function with lots of local minima, FindMinimum needs a decent starting value. NMinimize could be used here instead, but when you can give some reasonable starting information, FindMinimum is much faster.

• If FindMinimum does not work with the default method, don't give up -- the Newton and LevenbergMarquardt methods are both very robust.  There is lots of information on details of FIndMinimum and its methods in the help browser Unconstrained Optimization documentation.

• Having derivative information makes a big difference, especially for problems with high senstivity like this.