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

Code

CG

CG[A_,b_]:=Module[
{Tol=(10^(-10)),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]}
]

Test Matrices

The plan is to make matrices with specified distributions of eigenvalues.

RandomEvalSymPDMat[ EvalDist_]:= Module[
{evals,n,Q},
evals=Abs[Flatten[
Table[ 
RandomVariate[EvalDist[[i,1]],EvalDist[[i,2]]],
{i, Length[EvalDist]}]
]
];
n=Length[evals];
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]];
Q.DiagonalMatrix[evals].Transpose[Q]
]

A=RandomEvalSymMat[{
{NormalDistribution[10,0.00001],9},
{NormalDistribution[10,0.00001],45}
}];
n=Length[A];
xSol=RandomReal[{-1,1},n];
b=A.xSol;
{X,P}=CG[A,b];
ResData=Table[Norm[A.X[[i]]-b],{i, Length[X]}];
ErrData=Table[(1/2)X[[i]].A.X[[i]]-b.X[[i]]-((1/2)xSol.A.xSol-b.xSol),{i, Length[X]}];
{ListLogPlot[ResData],ListLogPlot[ErrData]}