Update Strategies for (Subscript[B, k])

Note: 
1) (Subscript[B, k])=Id in (Subscript[p, k])=-((Subscript[B, k])^(-1))∇(Subscript[f, k]) gives 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 (because (∇^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

SR1 is the simplest
	(Subscript[B, k+1])=(Subscript[B, k])+(((δf-(Subscript[B, k])δx)(((δf-(Subscript[B, k])δx))^T))/((((δf-(Subscript[B, k])δx))^T)δx))=(Subscript[B, k])+(((δf-(Subscript[B, k])δx)⊗(δf-(Subscript[B, k])δx))/((δf-(Subscript[B, k])δx).δx))
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])+(((δx-(Subscript[H, k])δf)⊗(δx-(Subscript[H, k])δf))/((δx-(Subscript[H, k])δf).δf))
Here δx=(Subscript[x, k+1])-(Subscript[x, k]) and δf=∇(Subscript[f, k+1])-∇(Subscript[f, k]).

Note: As always it is possible for this formula to have problems.  What if the inner product on the bottom is zero.

SR1Update[B_,df_,dx_]:= Module[
{v=df-B.dx},
B+ (1/(v.dx))Outer[Times,v,v]
]

n=2;
A= RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
A=A.A;
f[x_]:= x.A.x;
M=10;
B=IdentityMatrix[n];
{ContourPlot[f[{x1,x2}],
{x1,-M,M},
{x2,-M,M},
Contours->40],
ContourPlot[{x1,x2}.B.{x1,x2},
{x1,-M,M},
{x2,-M,M},
Contours->40]
}
{p11,p12}=RandomReal[{-1,1},{2,n}];
B=SR1Update[B, A.p11 -A.p12, p11-p12]
{ContourPlot[f[{x1,x2}],
{x1,-M,M},
{x2,-M,M},
Contours->40],
ContourPlot[{x1,x2}.B.{x1,x2},
{x1,-M,M},
{x2,-M,M},
Contours->40,
Epilog->{Red, Map[Point,{p11,p12}]}]
}
{p21,p22}=RandomReal[{-1,1},{2,n}];
B=SR1Update[B, A.p21 -A.p22, p21-p22]
{ContourPlot[f[{x1,x2}],
{x1,-M,M},
{x2,-M,M},
Contours->40],
ContourPlot[{x1,x2}.B.{x1,x2},
{x1,-M,M},
{x2,-M,M},
Contours->40,
Epilog->{Red, Map[Point,{p11,p12}], Green, Map[Point,{p21,p22}]}]
}

The right way to think about this is probably not with random points and certainly not by looking at pictures (because we are restricted to two dimensions).

n=134;
A= RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
A=A.A;

One way to see it work is by watching the norm converge!

B=IdentityMatrix[n];
NIter=145;
ListLogPlot[
Table[
{p1,p2}=RandomReal[{-1,1},{2,n}];
B=SR1Update[B,A.p1-A.p2,p1-p2];
Norm[B-A],
{NIter}],
PlotRange->All
]

Another is to see the eigenvalues converge.

B=IdentityMatrix[n];
NIter=145;
ListLogPlot[
Transpose[Table[
{p1,p2}=RandomReal[{-1,1},{2,n}];
B=SR1Update[B,A.p1-A.p2,p1-p2];
((Eigenvalues[A])/(Eigenvalues[B])),
{NIter}]],
PlotRange->All
]

ListLogPlot[ {Table[1,{145}],Eigenvalues[A]}]

B=IdentityMatrix[n];
NIter=145;
ListLogPlot[
Transpose[Table[
{p1,p2}=RandomReal[{-1,1},{2,n}];
B=SR1Update[B,A.p1-A.p2,p1-p2];
((Eigenvalues[A])/(Eigenvalues[B]))+((Eigenvalues[B])/(Eigenvalues[A])),
{NIter}]],
PlotRange->All
]

BFGS

BFGS is a little more complicated
	(Subscript[B, k+1])=(Subscript[B, k])+((δf⊗δf)/(δf.δx))-((((Subscript[B, k])δx)⊗((Subscript[B, k])δx))/(δx.((Subscript[B, k])δx)))
which is a rank two update becuase it has two extra piece composed of the outer product of vectors.  There are many rank two updates satisfying the secant condition.  This is the one that minimizes the change in B in a specific manner. It also has an explicit update formula for the inverse
	(Subscript[H, k+1])=(Subscript[H, k])+(((δf.δx + δf.(Subscript[H, k])δf))/(((δf.δx))^2))δx⊗δx-(1/(δx.δf))((Subscript[H, k])δf⊗δx +δx⊗(Subscript[H, k])δf )
As before δx=(Subscript[x, k+1])-(Subscript[x, k]) and δf=∇(Subscript[f, k+1])-∇(Subscript[f, k]).

Note: As always it is possible for this formula to have problems.  What if either inner product on the bottom is zero.

DFP

DFP is a second rank two update formula
	(Subscript[B, k+1])=(Id - γ  δf⊗δx)(Subscript[B, k])(Id - γ  δf⊗δx)+γ δf⊗δf
where γ=1/(δx.δf) It also has an inverse formulation
	(Subscript[H, k+1])=(Subscript[H, k])-(1/(δf.(Subscript[H, k])δf))((Subscript[H, k])δf)⊗((Subscript[H, k])δf)+(1/(δx.δf))δx⊗δx
DFP is the update that minimizes the change in H.