Systems of equations

In general optimization techniques are sometimes used to solve systems of equations.  If we want to solve
	(Subscript[g, 1])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]==0
	(Subscript[g, 2])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]==0
	(Subscript[g, 3])[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])]==0
then we can try to minimize
	f[(OverVector[x])]= (Subscript[w, 1])(((Subscript[g, 1])[(OverVector[x])])^2)+(Subscript[w, 2])(((Subscript[g, 2])[(OverVector[x])])^2)+(Subscript[w, 3])(((Subscript[g, 3])[(OverVector[x])])^2)=(OverVector[g])[(OverVector[x])].W.(OverVector[g])[(OverVector[x])]
and if we find a minimum of zero for any weights (Subscript[w, 1]),(Subscript[w, 2]) etc we have a solution of the equations.

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

f[x_:10]:= x
f[3]
f[]

TrustRadiusAdjust41[ {f_, m_}, Δ0_]:= Module[
{ΔMax=10.0,η=0.125},
ΔMax
]

TrustRadiusAdjust41[{f,m},s]

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_,Δ0_:1.0]:= Module[
{n=Length[x0],m,x,Δ,p,ρ,
MaxIter=23},
m={f[0][x0],f[1][x0],IdentityMatrix[n]};
(* Initializes triple describing model function *)
Do[
p = TrustRegionSubSolver[m,Δ,TrustRegionAlgorithm];
{x, m,Δ}=TrustRegionUpdate[{x,m, Δ},p,f,UpdateAlgorithm],
{MaxIter}]
]

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"]]
}]
]

Contour Pics and Constrained Optimization

f[x_,y_]:= (((x-y))^2)+ 4((( y+ 2 x))^2) + 2 x
ContourPlot[f[x,y],{x,-0.6,0.6},{y,-0.6,0.6},
AxesLabel->{"x","y"},
Contours->30,
Epilog->{Red,Circle[{0,0},0.5]}
]

f[x_,y_]:= (((x-y))^2)- 4((( y+ 2 x))^2) + 3 x
ContourPlot[f[x,y],{x,-0.6,0.6},{y,-0.6,0.6},
AxesLabel->{"x","y"},
Contours->30,
Epilog->{Red,Circle[{0,0},0.5]}
]

Clear[x,y,λ]
f[x_,y_]:= (((x-5y))^2)- 0.01((( y+ 2 x))^2) +12 x
{x0,y0} = {0.2,0.1};
Δ=0.7;
γ=2.6;
Hess=D[f[x,y],{{x,y},2}];
Gradf[x_,y_]=D[f[x,y],{{x,y}}];
g=-Gradf[x0,y0];
Temp=Simplify[LinearSolve[ Hess-2 λ(({{1, 0}, {0, 1}})), g]]
αSol[λ_]=Temp.Temp;
ConstrainedPointsData={{x0,y0}+Temp,λ}/.Solve[αSol[λ]==(Δ^2),λ];
ConstrainedPointsData =Select[ConstrainedPointsData,Im[#[[2]]]==0.0&];
MinPoints=Select[ConstrainedPointsData,#[[2]]<=0 &];
NewtonPoint={x0,y0}+LinearSolve[Hess, g];

{((149.698`-5.6899999999999995` λ)/(1.209999999999997`+25.95` λ-(λ^2))),((29.999000000000002`-1.495` λ)/(1.209999999999997`+25.95` λ-(λ^2)))}

ContourPlot[ f[x, y],
 {x, x0 - γ Δ, x0 + γ Δ},
 {y, y0 - γ Δ, y0 + γ Δ} ,
 Contours -> 20,
 Epilog -> {Red, 
   Thickness[0.01], Circle[{x0, y0}, Δ],
   PointSize[0.02], Point[{x0, y0}],
   Blue, Point[NewtonPoint],
   Green, Map[Point, ConstrainedPointsData[[All, 1]]],
   Cyan, Map[Point, MinPoints[[All, 1]]],
   Black, Table[
    {{x, y}, λ} = ConstrainedPointsData[[i]];
    Text[{λ, f[x, y]}, {x, y}, {-1, -1},
     Background -> White], {i, Length[ConstrainedPointsData]} ]
   }
 ]

More dimensions get more solutions.

Clear[x,y,z,f,λ]
f[x_,y_,z_]:= (((x-y))^2)- 4((( y+ 2 x))^2) + 3 x + (((z-x))^2)
{x0,y0,z0} = {0.2,0.1,0.2};
Δ=0.4;
γ=1.6;
Hess=D[f[x,y,z],{{x,y,z},2}];
Gradf[x_,y_,z_]=D[f[x,y,z],{{x,y,z}}];
g=-Gradf[x0,y0,z0];
Temp=LinearSolve[ Hess-2 λ(({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})), g];
αSol[λ_]=Simplify[Temp.Temp]
ConstrainedPointsData={{x0,y0,z0}+Temp,λ}/.Solve[αSol[λ]==(Δ^2),λ]
NewtonPoint={x0,y0,z0}+LinearSolve[Hess,g]

((10.170000000000002` (2.0052941376563047`-2.492257955039293` λ+(λ^2)) (18.230503663976485`-1.8086270007129175` λ+(λ^2)))/((((-2.2453251484312795`+λ))^2) (((-0.8400822807887431`+λ))^2) (((19.085407429220023`+λ))^2)))

{{{0.537149, 0.314906, 0.212088}, -26.8911}, {{-0.169841, -0.04934, 0.169751}, -11.2267}, {{-0.0878079, 0.290086, -0.00256866}, -0.420792}, {{-0.0986001 + 0.179271 I, 0.277048 - 0.349084 I, 0.697275 + 0.231934 I}, 1.35508 - 0.526121 I}, {{-0.0986001 - 0.179271 I, 0.277048 + 0.349084 I, 0.697275 - 0.231934 I}, 1.35508 + 0.526121 I}, {{0.265536, -0.293914, 0.17683}, 3.82844}}

{-0.125, 0.375, -0.125}

MinData = Select[ConstrainedPointsData, (z = #[[2]]; And[Im[z] == 0, z < 0]) &]

{{{0.537149, 0.314906, 0.212088}, -26.8911}, {{-0.169841, -0.04934, 0.169751}, -11.2267}, {{-0.0878079, 0.290086, -0.00256866}, -0.420792}}

TableForm[MinData]
NewtonPoint = {x0, y0, z0} + LinearSolve[Hess, g]

({{({{0.5371490311252809`}, {0.31490558473794217`}, {0.21208802960484746`}}), -26.891148693917824`}, {({{-0.16984127722752318`}, {-0.04934002868660661`}, {0.169751272909462`}}), -11.22667241899285`}, {({{-0.08780785229208715`}, {0.2900857167945954`}, {-0.0025686560889938703`}}), -0.42079163602509856`}})

{-0.125, 0.375, -0.125}