Conjugate Gradients

Krylov Space

Conjugate Gradients (and a bunch of similar algorithms) are based on Krylov spaces.  For A∈(ℝ^(n×n)) and q∈(ℝ^n)
	K(A,z,k)=Span[{q,A q,…,(A^(k-1))q}
All of the Krylov methods systematically generate least squares like solutions to A.x=b in Krylov spaces of increasing dimension.

Conjugate Gradient Ideas

Conjugate gradients only works for SPD matrices.  It is based on the non-standard inner product/norm pair
	(Subscript[(x,y), A])=x.A.y  and ||x(Subscript[||, A])=(Sqrt[x.A.x])
Vectors are conjugate if x.A.y=0.

Any set of n conjugate directions (Subscript[{(Subscript[p, i])}, i=1]^n) spans (ℝ^n) so we can write x=(∑_(i=1)^n)(Subscript[α, i])(Subscript[p, i]) for some coefficients α.

Conjugacy gives (Subscript[α, i])=(Subscript[(x,(Subscript[p, i])), A])/(Subscript[((Subscript[p, i]),(Subscript[p, i])), A])

The conjugate gradient “trick” is to build a set of conjugate directions as you go that rapidly build a good approximation

(Subscript[x, k])=(∑_(i=1)^k)(Subscript[α, i])(Subscript[p, i])

The result of the algorithm is a sequence of aproximation (Subscript[x, k]) that systematically decrease the residual

x.A.x/2-b.x

Defining F[x]= x.A.x/2-b.x  gives ∇F[x]=A.x-b.

Note: minimizing F solves A.x=b

The direction b-A.x is the “straight downhill direction” for F.

Conjugate Gradient Algorithm

So here is the plan.

(Subscript[x, 0])=0

Compute (Subscript[p, 1])=-∇F[(Subscript[x, 0])]=b  and normalize so that (Subscript[p, 1]).A.(Subscript[p, 1])=1.

Compute (Subscript[x, 1])=(Subscript[α, 1])(Subscript[p, 1])  so that F[(Subscript[x, 0])+(Subscript[α, 1])(Subscript[p, 1])]≤F[(Subscript[x, 0])+α (Subscript[p, 1])] for all α.

Differentiating gives A.((Subscript[x, 0])+(Subscript[α, 1])(Subscript[p, 1]))-b=0 ⟹ (Subscript[α, 1])=(Subscript[p, 1]).(b-A.(Subscript[x, 0])).

Compute (Subscript[p, 2])=-∇F[(Subscript[x, 1])]=b-A.(Subscript[x, 1])  and orth-normalize so that (Subscript[p, 2]).A.(Subscript[p, 1])=0 and (Subscript[p, 2]).A.(Subscript[p, 2])=1

Compute (Subscript[x, 2])=(Subscript[x, 1])+(Subscript[α, 2])(Subscript[p, 2]) so that F[(Subscript[x, 1])+(Subscript[α, 2])(Subscript[p, 2])]≤F[(Subscript[x, 0])+α (Subscript[p, 2])] for all α

Differentiating gives A.((Subscript[x, 1])+(Subscript[α, 2])(Subscript[p, 2]))-b=0 ⟹ (Subscript[α, 2])=(Subscript[p, 2]).(b-A.(Subscript[x, 1])).

…

Stop when ∇F[(Subscript[x, k])] is small enough.

Our first version is simply going to implement this directly.  The residual r=b-A.(Subscript[x, i]) is used repeatedly

CG[A_,b_]:=Module[
{Tol=(10^(-8)),n=Length[b],MaxIter=(10^3),P,X,p,r,x,i=1},
P= ConstantArray[0,{MaxIter,n}];
X= ConstantArray[0,{MaxIter,n}];
(* Initializing the Conjugate Directions using b*)
r=b;
p=r/(Sqrt[r.A.r]);
x=(p.r)p;
{X[[i]],P[[i]]}={x,p};
While[
And[Norm[r]>Tol,i++<MaxIter],
(* Computing residual *)
r=b-A.x;
(* Conjugating *)
p=r;
Do[p=p-(p.A.P[[j]]) P[[j]],{j,1,i-1}];
(*Do[p=p-(p.A.P[[j]]) P[[j]],{j,Max[1,i-12,i-1}];*)
p=p/(Sqrt[p.A.p]);
(* Updating x *)
x+=(p.r)p;
(* Storing Data *)
{X[[i]],P[[i]]}={x,p}
];
i=Min[i,MaxIter];
{Take[X,i],Take[P,i]}
]

Testing

A Random Test

n=124;
A=RandomReal[{0,1},{n,n}];
A=A.Transpose[A]+4.9IdentityMatrix[n];
xSol=RandomReal[{-1,1},n];
b=A.xSol;
{X,P}=CG[A,b];
ErrDat=Table[Norm[X[[i]]-xSol],{i, Length[X]}];
AErrDat=Table[res=X[[i]]-xSol;(Sqrt[res.A.res]),{i, Length[X]}];
GraphicsGrid[{{MatrixPlot[X],MatrixPlot[P]},
{ListLogPlot[ErrDat],ListLogPlot[AErrDat]}}]

Lets see what is happening with the “conjugacy”.

ConjTest=Table[
P[[i]].A.P[[j]],
{i,Length[P]},
{j, Length[P]}
]-IdentityMatrix[Length[P]];
Log10ConjErr=Log[10,Abs[ConjTest]];
{MatrixPlot[ConjTest],ListContourPlot[Log10ConjErr]}

Testing

A 1D Laplacian Test

n=123;
A=SparseArray[{Band[{1,1}]->2.0,Band[{1,2}]->1.0, Band[{2,1}]->1.0},{n,n}];
xSol=RandomReal[{-1,1},n];
b=A.xSol;
{X,P}=CG[A,b];
ErrDat=Table[Norm[X[[i]]-xSol],{i, Length[X]}];
AErrDat=Table[res=X[[i]]-xSol;(Sqrt[res.A.res]),{i, Length[X]}];
GraphicsGrid[{{MatrixPlot[X],MatrixPlot[P]},
{ListLogPlot[ErrDat,PlotRange->All],
ListLogPlot[AErrDat,PlotRange->All]}}]

Lets see what is happening with the “conjugacy”.

ConjTest=Table[
P[[i]].A.P[[j]],
{i,Length[P]},
{j, Length[P]}
]-IdentityMatrix[Length[P]];
Log10ConjErr=Log[10,Abs[ConjTest]];
{MatrixPlot[ConjTest],ListContourPlot[Log10ConjErr]}