(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



Clear[HessianUpdate]
SetAttributes[HessianUpdate, HoldFirst]


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]
];
]


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


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


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];
]


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


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


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



