Line Search: Steepest Descent and Newton Direction

f[{x_,y_}]:=Exp[Sin[50 x]]+Sin[60 Exp[y]]+Sin[70 Sin[x]]+Sin[Sin[80 y]]-Sin[10 (x+y)]+1/4 (x^2+y^2);
Df[{x_,y_}]=D[f[{x,y}],{{x,y}}];
DDf[{x_,y_}]=D[Df[{x,y}],{{x,y}}];
M=6;
ϵ=(10^(-1));
{x0,y0}= RandomReal[{-M,M},{2}];
pSD=-Df[{x0,y0}];
Hess = DDf[{x0,y0}];
pN = LinearSolve[Hess, pSD];
NormpSD =Norm[pSD]; 
ContourPlot[f[{x,y}],{x, -M,M},{y,-M,M},
MaxRecursion->2,
Epilog->{
Red, Point[{x0,y0}], 
Pink, Arrow[ {{x0,y0},{x0,y0}+M  (pSD/NormpSD)}],
Green,  Arrow[ {{x0,y0},{x0,y0}+pN}]
}
]
ContourPlot[f[{x,y}],{x, x0-ϵ,x0+ϵ},{y,y0-ϵ,y0+ϵ},
MaxRecursion->2,
Epilog->{
Red, Point[{x0,y0}], 
Pink, Arrow[ {{x0,y0},{x0,y0}+ϵ (pSD/NormpSD)}],
Green,  Arrow[ {{x0,y0},{x0,y0}+pN}]
},
PlotLabel->StringForm["Norm pSD = ``",NormpSD]
]

c1=(10^(-3));
c2=0.9;
p=pSD/NormpSD; (* p = pSD/NormpSD or p=pN*)
p=pN; 
Wolfe2[α_,c2_]:=p. Df[{x0,y0}+ α p ]≥c2 p.Df[{x0,y0}] ;
StrongWolfe2[α_,c2_]:=Abs[p. Df[{x0,y0}+ α p ]]≤Abs[c2 p.Df[{x0,y0}]]; 
αMax=2.2;
Show[
Plot[ {
f[{x0,y0}+ α p ],
f[{x0,y0}]- α,
f[{x0,y0}]- c1 α
}  , {α,0,αMax}
],
Plot[
f[{x0,y0}+ α p ], {α,0, αMax},
RegionFunction:>Function[α,StrongWolfe2[α,c2]],
AxesLabel->{"α","f"},
PlotStyle->{Thickness[0.01],Green}
],
Plot[
f[{x0,y0}+ α p ], {α,0, αMax},
RegionFunction:>Function[α,Wolfe2[α,c2]],
AxesLabel->{"α","f"},
PlotStyle->{Thickness[0.005],Red}
]
]

There are other similar conditions.  They all have the same flavor:

At least some (usually small) fraction of the decrease the gradient predicts you should be able to attain.  It is not enough to say (Subscript[f, new])<(Subscript[f, old]). These always permit tiny steps because a sufficiently small step does enjoy the tangent line decrease.

At least some change in slope so that the slope at the new point is (at least) significantly less negative than the old point.

Adjustment Strategies Search

If we are using the Wolfe Conditions we need to try a shorter step if it fails the first condition and a longer step if it passes the first and fails the second.  So we need two adjustement startegies. There are lots of choices but a common one is:  if you need a shorter step multiply by a constant less than one; if you need a longer step multiply by a constant bigger than one.

Here is a very simple implementation suitable for a Quasi Newton scheme.  Note it assumes that p is a descent direction and starts with 1 because that is a suitable length for a Newton step.

Step[x0_,p_,{f_,Df_}]:=Module[
{αMin=(10^(-6)),α=1,c1=(10^(-4)),c2=0.9,cBig=1.5,cSmall=0.5,fOld=f[x0],fVal,DfOld = Df[x0],DfVal},
fVal=f[x0+α p];
DfVal=Df[x0+α p];
Do[
fVal=f[x0+α p];DfVal=Df[x0+α p];
Print[{fVal≤fOld + c1 p.DfOld,p.DfVal>=c2 p.DfOld}];
If[fVal≤fOld + c1 p.DfOld,If[p.DfVal>=c2 p.DfOld,Break[],α=cBig*α],α=cSmall*α],
{10}
];
α
]

Step[{x0,y0},p,{f, Df}]

Efficiency issues.  This way of doing things has unneccessary function and gradient evaluations.

Parrallelization Strategies

Suppose you could compute a number of functions or gradients at the same time.  This means you could for instance compute a bunch of points and take the lowest.  Here is a simple proposal possibly suitable for a QuasiNewton iteration.  It assumes that p is a descent direction.
Compute a bunch of points along the line.  Select the α with the smallest fvalue.  This is your step provided it is not the last point.  If there are no points satisfying Wolfe1 take a smaller step using a standard strategy.  If the smallest value satisfies Wolfe1 but is the last value you computed then take a longer step by recentering your computation on the last point. Continue as usual.

ϵ=0.1;
p=(pSD/pSD);
pts = Table[ α=(((1+ϵ))^n);{α, f[{x0,y0}+α p]},{n,-7,8}];
Plot[{
f[{x0,y0}+α p],
f[{x0,y0}]+α c1  p. Df[{x0,y0}]
},{α, 0,αMax},
Epilog:>Map[Point,pts],
PlotRange->All
]

Update Strategies for (Subscript[B, k])

Note: 
1) (Subscript[B, k])=Id in (Subscript[p, k])=-((Subscript[B, k])^(-1))∇(Subscript[f, k]) is steepest descent.  
2) Everytime we compute a ∇f we get some information about (∇^2)f near the point (Subscript[x, k])
3) If we assume that (∇^2)f does not change fast we could build up an approximation of (∇^2)f from these evaluations.
4) We NEED to make sure our approximation is always symmetric ((∇^2)f is) and positive definate so that -((Subscript[B, k])^(-1))∇(Subscript[f, k]) is a descent direction.
5) Info we have is that ∇f((OverVector[b]))-∇f((OverVector[a]))=(∇^2)f(s(OverVector[a])+(1-s)(OverVector[b])).((OverVector[b])-(OverVector[a])) for some s satisfying 0≤s≤1.  Assuming that (∇^2) is approximately constant gives the "secant equation"  ∇(Subscript[f, k+1])-∇(Subscript[f, k])=(Subscript[B, k+1])((Subscript[x, k+1])-(Subscript[x, k])) for the updated approximation to (∇^2)f to be used for the next step k+1->k+2.  This is one vector equation in (ℝ^n) for the ((n(n-1))/2) entries in a symmetric matrix. There is clearly more than one answer.  There are three common (SR1, BFGS, and DFP) update formulas and others that are less common.

SR1 is the simplest
	(Subscript[B, k+1])=(Subscript[B, k])+((((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k]))((((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k])))^T))/(((((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k])))^T)(Subscript[δx, k])))=(Subscript[B, k])+((((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k]))⊗((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k])))/(((Subscript[δf, k])-(Subscript[B, k])(Subscript[δx, k])).(Subscript[δx, k])))
which is a rank one update becuase it has one extra piece composed of the outer product of vectors.  It is the only symmetric rank one update formula. I has a conveniently easy inverse update formula
	(Subscript[H, k+1])=(Subscript[H, k])+((((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k]))((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k])))/(((((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k])))^T)(Subscript[δf, k])))=(Subscript[H, k])+((((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k]))⊗((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k])))/(((Subscript[δx, k])-(Subscript[H, k])(Subscript[δf, k])).(Subscript[δf, k])))
BFGS and DFP are rank two updates and have two extra terms in them.  We will look at them Wednesday.