Code Pieces

Test Functions

Quadratic Test Functions

Clear[x,y,QuadSPD,QuadNonSPD]
QuadSPD[0][{x_,y_}]:= 12.3 (((x- 2 y))^2)+ 1.2 (((x+y))^2)+ 3 x - 4 y + 1.1
QuadSPD[1][{x_,y_}]:={3+24.6` (x-2 y)+2.4` (x+y),-4-49.2` (x-2 y)+2.4` (x+y)}
QuadSPD[2][{x_,y_}]:= {{27.0`,-46.8`},{-46.8`,100.8`}}
QuadNonSPD[0][{x_,y_}]:= 1.1`+3 x+12.3` (((x-2 y))^2)-4 y-0.1` (((x+y))^2)
QuadNonSPD[1][{x_,y_}]:={3+24.6` (x-2 y)-0.2` (x+y),-4-49.2` (x-2 y)-0.2` (x+y)}
QuadNonSPD[2][{x_,y_}]:={{24.4`,-49.4`},{-49.4`,98.2`}}

Trust region control

The following is a direct implementation of Alg 4.1 on p69 of Nocedal

Clear[TrustRadiusAdjust41]
TrustRadiusAdjust41[ {f_, m_},{xIn_,pIn_}, ΔIn_,ΔMax_:10,η_:0.125]:= Module[
{ρ, mFun,pp,xOut,ΔOut},
mFun[pp_]:=m[[1]]+m[[2]].pp + (1/2)pp.m[[3]].pp;
ρ=((f[0][xIn]-f[0][xIn+pIn])/(mFun[0*pIn]-mFun[pIn]));
Which[
ρ<0.25, ΔOut=0.25ΔIn,
And[ ρ>0.75, Norm[pIn]==Δ0],  ΔOut=Min[2 ΔIn,ΔMax],
True, ΔOut=ΔIn
];
If[ρ>η, xOut=xIn+pIn, xOut=xIn];
{xOut, ΔOut}
]

f=QuadSPD;
y0={1.2,3.4};
p0={0.1,0.1};
Δ0=0.1;
f[0][y0]
m={f[0][y0],f[1][y0],f[2][y0]};
TrustRadiusAdjust41[{f,m},{y0,p0},Δ0]

ϵ=10;
M=25;
Show[
ContourPlot[
f[0][y0+{p1,p2}],
{p1,-M,M},
{p2,-M,M},
Contours->30,
ColorFunction->Hue
],
ContourPlot[
mFun[{p1,p2}],
{p1,-ϵ,ϵ},
{p2,-ϵ,ϵ},
Contours->30 (ϵ/M)
]
]

Hessian Updates and Approximations

I need a number of Hessian updates etc. I need BFGS, DFP, SR1, Newton, and Indefinate Newton.  If I get fancy I should worry about efficiency.  The general fomr for all of these is going to be
	HessianUpdate[HessianAddress, updateinfo, Tag]
This should be sufficiently flexible.

The following makes all the variants of the Hessian update hold the first argument.

Clear[HessianUpdate]
SetAttributes[HessianUpdate, HoldFirst]

BFGS

HessianUpdate[B_, {y_,s_},"BFGS_B"]:= Module[
{sy=s.y,BS=B.s,Tol=(10^(-6))},
If[sy>Tol , 
B=B -(1/(s.BS)) Outer[Times, BS, BS]+(1/sy)Outer[Times, y,y]
];
]

DFP

HessianUpdate[B_, {y_,s_},"DFP_B"]:= Module[
{sy=s.y,Tol=(10^(-6))},
If[ sy>Tol,  
ρ=(1/sy);
M=IdentityMatrix[Length[y]]-ρ Outer[Times, y,s];
B=M.B.Transpose[M]+ ρ Outer[Times, y,y]
]
]

SR1

HessianUpdate[B_, {y_,s_},"SR1_B"]:= Module[
{yMinusBs=y-B.s, sDotyMinusBs,Tol=(10^(-9))},
sDotyMinusBs = s.yMinusBs;
If[ Abs[sDotyMinusBs]>Tol, 
B= B + (1/sDotyMinusBs)Outer[Times, yMinusBs, yMinusBs]
]
]

Newton and ModNewton

Clear[HessianComp]
SetAttributes[HessianComp, HoldFirst]
HessianComp[B_, {f_,x_},"Newton_B"]:= Module[
{},
B=f[2][x];
]

HessianComp[H_, {f_,x_},"Newton_H"]:= Module[
{B},
B=f[2][x];
H=LinearSolve[B];
]

Testing Hessian Computations

B=IdentityMatrix[2];
x0={1.2,3.4};
f=QuadSPD;
HessianComp[B,{f,x0},"Newton_B"]
{B,f[2][x0]}

H=IdentityMatrix[2];
x0={1.2,3.4};
u={1.2,3.4};
f=QuadSPD;
HessianComp[H,{f,x0},"Newton_H"]
{H[u],LinearSolve[f[2][x0], u]}

Testing Hessian Updates

First BFGS etc. on a Quadratic function.

x0={1.2,3.4};
f=QuadSPD;
MaxIter=120;
Table[
B=IdentityMatrix[2];
xOld=RandomReal[{-1,1},{2}];
DfOld=f[1][xOld];
ListLogLogPlot[
Table[
xNew=RandomReal[{-1,1},{2}];
DfNew=f[1][xNew];
s=xNew-xOld;
y=DfNew-DfOld;
Sow[s.y,"dotprods"];
HessianUpdate[B,{y,s},QNewtScheme];
{xOld, DfOld}={xNew,DfNew};
Norm[B-f[2][x0]],
{MaxIter}],
PlotLabel->QNewtScheme,
AxesLabel->{"iter","MatrixNorm"}
],
{QNewtScheme,{"BFGS_B","DFP_B","SR1_B"}}
]

Second BFGS etc. on a Quadratic non SPD function.

x0={1.2,3.4};
f=QuadNonSPD;
MaxIter=120;
Table[
B=IdentityMatrix[2];
xOld=RandomReal[{-1,1},{2}];
DfOld=f[1][xOld];
ListLogLogPlot[
Table[
xNew=RandomReal[{-1,1},{2}];
DfNew=f[1][xNew];
s=xNew-xOld;
y=DfNew-DfOld;
Sow[s.y,"dotprods"];
HessianUpdate[B,{y,s},QNewtScheme];
{xOld, DfOld}={xNew,DfNew};
Norm[B-f[2][x0]],
{MaxIter}],
PlotLabel->QNewtScheme,
AxesLabel->{"iter","MatrixNorm"}
],
{QNewtScheme,{"BFGS_B","DFP_B","SR1_B"}}
]

x0={1.2,23.4};
f=QuadNonSPD;
MaxIter=12;
Table[
B=IdentityMatrix[2];
xOld=RandomReal[{-1,1},{2}];
DfOld=f[1][xOld];
ListLogPlot[
Table[
xNew=RandomReal[{-1,1},{2}];
DfNew=f[1][xNew];
s=xNew-xOld;
y=DfNew-DfOld;
HessianUpdate[B,{y,s},QNewtScheme];
{xOld, DfOld}={xNew,DfNew};
Norm[B-f[2][x0]],
{MaxIter}],
PlotLabel->QNewtScheme,
AxesLabel->{"iter","MatrixNorm"}
],
{QNewtScheme,{"BFGS_B","DFP_B","SR1_B"}}
]

Main Drivers

The following is a fairly direct implementation of Alg 4.1 on p69 of Nocedal showing how I would identify subtasks.  The first thing to do is to identify the data for these algorithms.  I think the primary data at any stage is
	{(Subscript[x, k]),{(Subscript[f, k]),∇(Subscript[f, k]),(Subscript[B, k])},(Subscript[Δ, k])} 
where (Subscript[B, k])≈(∇^2)(Subscript[f, k])  with model function 
	mk[p]=fk+Dfk.p+(1/2)p.Bk.p 
and Δk is the current trust region radius. In the code the model function is described by the triple m={f,Df,B}

TrustRegionDriver41[ f_, x0_,
{ModelFunctionInitializationAlgorithm_,TrustRegionSubSolverAlgorithm_, TrustRegionUpdateAlgorithm_},
Δ0_:1.0]:= Module[
{n=Length[x0],m,x,Δ,p,ρ,
MaxIter=23},
{x,m,Δ}={x0, ModelFunctionInitialize[f,x0, ModelFunctionInitializationAlgorithm],Δ0};
(* Initializes triple describing model function *)
Do[
p = TrustRegionSubSolver[m,Δ,TrustRegionSubSolverAlgorithm];
{x, m,Δ}=ModelFunctionUpdate[{x,m, Δ},p,f,ModelFunctionUpdateAlgorithm],
{MaxIter}];
{x,m,Δ}
]

Sub Problem Solvers

TrustRegionSubSolver[m_,Δ_,"Cauchy"]:= Module[
{g=m[[2]], gBg,τ}, (* m = {f,g,B} *)
gBg=g.m[[3]].g;         (* m = {f,g,B} *)
(* See 4.11 and 4.12 p72 Nocedal *)
τ= Which[
gBg≤0,1,
True, Min[1, ((((Norm[g])^3) )/(Δ gBg))]
];
- τ (Δ/(Norm[g]))g
]

TrustRegionSubSolver[m_,Δ_,"Newton"]:= Module[
{pNewt,NormpNewt}, (* m = {f,g,B} *)
pNewt=-LinearSolve[m[[3]],m[[2]]];         (* m = {f,g,B} *)
NormpNewt=Norm[pNewt];
If[NormpNewt<Δ,pNewt,Δ (pNewt/NormpNewt)]
]

TrustRegionSubSolver[m_,Δ_,"DogLeg"]:= Module[
{g=m[[2]], gBg,τ,pNewton,pCauchy,NormpCauchy,pCauchyToNewton}, (* m = {f,g,B} *)
gBg=g.m[[3]].g;         (* m = {f,g,B} *)
pCauchy=-(((Norm[g])^2)/gBg)g;
NormpCauchy=Norm[pCauchy];
If[ NormpCauchy>Δ, Return[Δ (pCauchy/NormpCauchy)]];
pNewton=-LinearSolve[m[[3]],m[[2]]];
pCauchyToNewton=pNewton-pCauchy;
If[Norm[pNewton]≤Δ, Return[pNewton]];
{a,b,c}= {pCauchyToNewton.pCauchyToNewton,2 pCauchyToNewton.pCauchy, pCauchy.pCauchy-(Δ^2)} ;
t= ((-b+(Sqrt[(b^2)-4 a c]))/(2a));
pCauchy+t pCauchyToNewton
]

TrustRegionSubSolver[m_,Δ_,"2Dv1"]:= Module[
{pNewt,pSD,α1,α2},
pSD=-m[[2]];                                               (* m = {f,g,B} *)
pSD=pSD/Norm[pSD];
pNewt=LinearSolve[m[[3]],pSD];     (* m = {f,g,B} *)
pNewt= pNewt - (pNewt.pSD)pSD;
pNewt=pNewt/Norm[pNewt];
(α1 pNewt+ α2 pSD)/.FindMinimum[ {
(1/2)(α1 pNewt+ α2 pSD).m[[3]].(α1 pNewt+ α2 pSD)+m[[2]].(α1 pNewt+ α2 pSD),
((α1)^2)+((α2)^2)<= (Δ^2)},
{α1,α2}][[2]]
]

TrustRegionSubSolver[m_,Δ_,"2Dv2"]:= Module[
{pNewt,pSD,α1,α2, M=m[[3]],M2,g2,p1,p2,λ},
pSD=-m[[2]];                                               (* m = {f,g,B} *)
pSD=pSD/Norm[pSD];
pNewt=LinearSolve[m[[3]],pSD];     (* m = {f,g,B} *)
pNewt= pNewt - (pNewt.pSD)pSD;
pNewt=pNewt/Norm[pNewt];
(α1 pNewt+ α2 pSD)/.FindMinimum[ {
(1/2)(α1 pNewt+ α2 pSD).m[[3]].(α1 pNewt+ α2 pSD)+m[[2]].(α1 pNewt+ α2 pSD),
((α1)^2)+((α2)^2)<= (Δ^2)},
{α1,α2}][[2]];
M2={
{pNewt.M.pNewt, pSD.M.pNewt},
{pSD.M.pNewt,pSD.M.pSD}
};
g2={m[[2]].pNewt,m[[2]].pSD};
Data={{α1 pNewt+α2 pSD},λ}/.Solve[{((α1)^2)+((α2)^2)==(Δ^2), M2.{α1,α2}+g2 - λ {α1,α2}=={0,0}}, {α1,α2,λ}];
Select[Data,And[Im[#[[2]]]==0.0,#[[2]]<0]&][[All,1]]
]

m={1.2, {1.2,3.4}, {{1.2,0.1},{0.1, 23.4}}};
Δ=0.7; 
TrustRegionSubSolver[m, Δ, "2Dv1"]
Fred=TrustRegionSubSolver[m, Δ, "2Dv2"]

model[p_]:= m[[1]]+m[[2]].p + (1/2)p.m[[3]].p
ContourPlot[ 
model[{p1,p2}],
{p1,-Δ,Δ},
{p2, -Δ,Δ},
Contours->40,
Epilog->{Thickness[0.01],Circle[{0,0},Δ],
PointSize[0.04],Point[{0,0}],
Green, Point[TrustRegionSubSolver[m, Δ, "Cauchy"]],
Blue, Point[TrustRegionSubSolver[m, Δ, "Newton"]],
PointSize[0.03],
Pink,Map[Point,TrustRegionSubSolver[m, Δ, "2Dv2"]],
PointSize[0.02],
Cyan,Point[TrustRegionSubSolver[m, Δ, "DogLeg"]]
},
ColorFunction->Hue
]

M2

Testing the sub problem solvers.

On a 2D SPD function

m={1.2, {1.2,3.4}, {{1.2,0.1},{0.1, 23.4}}};
Δ=0.7; 
model[p_]:= m[[1]]+m[[2]].p + (1/2)p.m[[3]].p
Manipulate[
ContourPlot[ 
model[{p1,p2}],
{p1,-Δ,Δ},
{p2, -Δ,Δ},
Contours->40,
Epilog->{Thickness[0.01],Circle[{0,0},Δ],
PointSize[0.04],Point[{0,0}],
Green, Point[TrustRegionSubSolver[m, Δ, "Cauchy"]],
Blue, Point[TrustRegionSubSolver[m, Δ, "Newton"]],
PointSize[0.03],
Pink,Point[TrustRegionSubSolver[m, Δ, "2Dv1"]],
PointSize[0.02],
Cyan,Point[TrustRegionSubSolver[m, Δ, "DogLeg"]]
},
ColorFunction->Hue
],
{Δ, 0.1, 1.4,0.1}]

On a 2D  non PD function

m={1.2, {1.2,3.4}, {{-1.2,0.1},{0.1, 23.4}}};
Δ=0.7; 
model[p_]:= m[[1]]+m[[2]].p + (1/2)p.m[[3]].p
Manipulate[
ContourPlot[ 
model[{p1,p2}],
{p1,-Δ,Δ},
{p2, -Δ,Δ},
Contours->40,
Epilog->{Thickness[0.01],Circle[{0,0},Δ],
PointSize[0.04],Point[{0,0}],
Green, Point[TrustRegionSubSolver[m, Δ, "Cauchy"]],
Blue, Point[TrustRegionSubSolver[m, Δ, "Newton"]],
PointSize[0.03],
Pink,Point[TrustRegionSubSolver[m, Δ, "2Dv1"]],
PointSize[0.02],
Cyan,Point[TrustRegionSubSolver[m, Δ, "DogLeg"]]
},
ColorFunction->Hue
],
{Δ, 0.1, 1.4,0.1}]

In 3D

It is much harder to see what is going on in 3D.

m={1.2, {1.2,3.4,5.6},
 {
{1.2,0.1,0.1},
{0.1, 3.4,0.2},
{0.1, 0.2, 6.5}
}
};
Δ=1.3; 
model[p_]:= m[[1]]+m[[2]].p + (1/2)p.m[[3]].p

pN=LinearSolve[ m[[3]],-m[[2]]]
Norm[pN]
Show[
SphericalPlot3D[
Δ,{θ,0,Pi},{ϕ,0,2Pi},
ColorFunction->Function[{x,y,z,θ,ϕ,r},Hue[model[{x,y,z}]/10]],
PlotStyle->Opacity[0.5]
],
Graphics3D[{
PointSize[0.04],Point[{0,0,0}],
Green, Point[TrustRegionSubSolver[m, Δ, "Cauchy"]],
Blue, Point[TrustRegionSubSolver[m, Δ, "Newton"]],
PointSize[0.03],
Pink,Point[TrustRegionSubSolver[m, Δ, "2Dv1"]],
PointSize[0.02],
Cyan,Point[TrustRegionSubSolver[m, Δ, "DogLeg"]]
}]
]

Runs