Linear Conjugate Gradient: Convergence etc.

We implement the algorithm on p108 directly.

CGp108[A_,b_,x0_] := Module[
{x=x0,r,p,α,β,k, n=Length[A]+1},
r=A.x-b;
p=-r;
Table[
α =  - ((r.p)/(p.(A.p)));
x=x + α p;
r = A.x-b;
β=((r.(A.p))/(p.(A.p)));
p=β p - r;
x,
{k,1, n}]
]

Now we are going to test it.  Remember it only works for Symmetric Positive Definate matrices. First we need to make SPD matrices.

n=10;
A= RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
{λs,Qs}=Eigensystem[A];
A-Transpose[Qs].DiagonalMatrix[λs].Qs;
APlus=Transpose[Qs].DiagonalMatrix[Abs[λs]].Qs;
b=RandomReal[{-1,1},n];
xInf = LinearSolve[APlus, b];
x0=RandomReal[{-1,1},n];
Data=CGp108[APlus,b,x0];
ListLogPlot[
Table[ Norm[Data[[k]]-xInf],{k,1,n+1}]
]

A slightly larger problem.

n=300;
A= RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
{λs,Qs}=Eigensystem[A];
A-Transpose[Qs].DiagonalMatrix[λs].Qs;
APlus=Transpose[Qs].DiagonalMatrix[Abs[λs]].Qs;
b=RandomReal[{-1,1},n];
xInf = LinearSolve[APlus, b];
x0=RandomReal[{-1,1},n];
Data=CGp108[APlus,b,x0];
ListLogPlot[
Table[ Norm[Data[[k]]-xInf],{k,1,n}]
]

n=500;
A= RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
{λs,Qs}=Eigensystem[A];
A-Transpose[Qs].DiagonalMatrix[λs].Qs;
λs=RandomChoice[Range[1, 40,1],n];
APlus=Transpose[Qs].DiagonalMatrix[Abs[λs]].Qs;
b=RandomReal[{-1,1},n];
xInf = LinearSolve[APlus, b];
x0=RandomReal[{-1,1},n];
Data=CGp108[APlus,b,x0];
ListLogPlot[
Table[ Norm[Data[[k]]-xInf],{k,1,n}],
PlotRange->All
]

Note this is not the correct fancy error norm.  But the drop off still happens.

Linear Conjugate Gradient: p112 Implementation

We implement the algorithm on p112 directly.

Clear[RandomSPDMat]
RandomSPDMat[n_, Dist_:UniformDistribution[{1,10}]]:= Module[
{Q=RandomReal[{-1,1},{n,n}],λs=Abs[RandomReal[Dist,n]]},
Q=Eigenvectors[Q+Transpose[Q]];
Transpose[Q].DiagonalMatrix[λs].Q
]

Clear[CGp112Step,CGp112]
CGp112Step[{xIn_,rIn_,pIn_},{A_,b_}]:= Module[
{Ap=A.pIn,xOut,rOut,pOut,α,β},
α=((rIn.rIn)/(pIn.Ap));         (* value of α giving the min in direction p *)
xOut=xIn+α pIn;(* New x with this α in deirection p *)
rOut=rIn+α Ap;   (* New gradient *)
β=((rOut.rOut)/(rIn.rIn));
pOut=-rOut+ β pIn;
{xOut,rOut,pOut}
];

CGp112[x0_,{A_,b_},Tol_] := Module[
{r0,p},
r0=A.x0-b;
NestWhile[
CGp112Step[#,{A,b}]&, 
{x0,r0,-r0}, 
(Sow[Norm[#[[2]]],"CGp112"]>Tol)&][[1]]
]

Now we are going to test it.  Remember it only works for Symmetric Positive Definate matrices. So we make an SPD matrix and test the thing.

n=30;
A= RandomSPDMat[n];
b=RandomReal[{-1,1},n];
x0=RandomReal[{-1,1},n];
Tol=(10^(-12));
ListLogPlot[
Reap[xCG=CGp112[x0,{A,b},Tol];][[2,1]]
]
xSol=LinearSolve[A,b];
xCG-xSol

A slightly larger problem.

n=300;
A= RandomSPDMat[n];
b=RandomReal[{-1,1},n];
x0=RandomReal[{-1,1},n];
Tol=(10^(-12));
ListLogPlot[
Reap[xCG=CGp112[x0,{A,b},Tol];][[2,1]]
]

NonLinear Conjugate Gradient: FR p121 Preliminary Implementation

We implement the algorithm on p121 directly.  It needs a line search algorithm.  This could be any of our earlier line search calls. For clarity we are going to use the built-in 1D mathematica FindMinimum command.

Clear[f, Df, x1,x2]
f2D[{x1_,x2_}]:= (x1^2)+ (x2^2)- 2Sin[x1 x2] +  2x1 - 3 x2 ;
Df2D[{x1_,x2_}]= D[ f2D[{x1,x2}], {{x1,x2}}];
M=5;
m=0.5;
x0 =RandomReal[m{-M,M},{2}];
p0= -Df2D[x0];
α =α/. FindMinimum[{
f2D[x0+α p0],
0≤α≤1},
{α,0.0}
] [[2]];

xNew =x0+α p0;
ContourPlot[ f2D[{x1,x2}],{x1, -M, M},{x2, -M, M},
Contours->60,
Epilog->{
{Red,Arrow[ {x0, x0 + p0}]},
{Green,Point[x0]},
{Yellow, Point[xNew]}
}
]

Clear[CGFRStep,CGFR]
CGFRStep[{xIn_,rIn_,pIn_},{f_,Df_}]:= Module[
{f0=f[xIn],xOut,rOut,pOut,α,β},
α=α/.FindMinimum[{  (* Mma Line Search *)
f[xIn+α pIn],
0≤α≤1},
{α,0.0}
] [[2]];                    (* OK α in direction p *)
xOut =xIn+α pIn;(* New x with OK α in direction p *)
rOut= Df[xOut];    (* New gradient *)
β=((rOut.rOut)/(rIn.rIn));
pOut=-rOut+ β pIn;
{xOut,rOut,pOut}
];

CGFR[x0_,{f_,Df_},Tol_] := Module[
{r0,p},
r0=A.x0-b;
NestWhile[
CGFRStep[#,{f,Df}, ]&, 
{x0,r0,-r0}, 
(Norm[#[[2]]]>Tol)&][[1]]
]

First we are going to test on a small sample problem.

Clear[f, Df, x1,x2]
f2D[{x1_,x2_}]:= (x1^2)+ (x2^2)- 2Sin[x1 x2] +  2x1 - 3 x2 ;
Df2D[{x1_,x2_}]= D[ f2D[{x1,x2}], {{x1,x2}}];
M=5;
m=0.5;
x0 =RandomReal[m{-M,M},{2}];
{x,r,p}={x0, Df2D[x0], -Df2D[x0]};
Data=Table[
{x,r,p}=CGFRStep[{x,r,p},{f2D, Df2D}],
{24}];
xs=Join[{x0},Data[[All,1]]]

ContourPlot[ f2D[{x1,x2}],{x1, -M, M},{x2, -M, M},
Contours->60,
Epilog->{
{Red,Line[ xs]},
{Green,Point[x0]}
}
]

We need a bigger test problem.

n=30;
A= RandomSPDMat[n];
b=RandomReal[{-1,1},n];
fBig[x_]:= Module[
{x1,x2},
{x1,x2}=Take[x,2];
(1/2)x.A.x - b.x + f2D[{x1,x2}]
]
DfBig[x_]:=Module[
{x1,x2,Der, Der2D},
{x1,x2}=Take[x,2];
Der2D=Df2D[{x1,x2}];
Der=A.x-b;
Der[[{1,2}]]+=Der2D;
Der
]
M=3;
x0 =RandomReal[{-M,M},{n}];
{x,r,p}={x0, DfBig[x0], -DfBig[x0]};
Data=Table[
{x,r,p}=CGFRStep[{x,r,p},{fBig, DfBig}],
{15}];
xs=Data[[All,1]];
rs=Data[[All,2]];
ListLogPlot[Map[Norm,rs],
PlotRange->All]