Krylov Subspaces and Arnoldi Iteration

Pick a random vector (OverVector[b]).  Hit it with a matrix A the chances are that you have caught some information about the largest eigenvalue. Hit it again and you have amplified whatever information you got about the largest eigenvalue and picked up some information about the second alrgest eigenvalue, …

The Arnoldi iteration filters and cleans this information at each n for the Krylov subspaces
	(Subscript[𝒦, n])=⟨b,Ab,(A^2)b,…,(A^(n-1))b⟩
by maintaining a list of orthogonal vectors (generated using modified gram-schmidt) that span (Subscript[𝒦, n]) i.e.
	(Subscript[𝒦, n])=⟨(Subscript[q, 1]),(Subscript[q, 2]),…,(Subscript[q, n])⟩.
The code is fairly easy!  Here is a non-function version that does not store the H array and leaves space for m  vectors in the Q.

m=12;
nMax=4;
A = RandomReal[{-1,1},{m,m}];
b = RandomReal[{-1,1},{m}];
Q=0*A;
Q[[All,1]]=b/Norm[b];
Do[
v=A.Q[[All,n]];
Do[
h = Q[[All,j]].v;
v=v-h Q[[All,j]],
{j,1,n}
];
h=Norm[v];
Q[[All,n+1]]=v/h,
{n,1, nMax-1}
]
MatrixForm[Q]
MatrixForm[Transpose[Q].Q]

Here s a function that stores Q and H in appropriate;y sized matrices.

Arnoldi[A_,b_,nMax_]:=Module[
{m=Length[A],Q,H,n,j,h,v},
Q=ConstantArray[0,{m,nMax+1}];
H=ConstantArray[0,{nMax+1,nMax}];
Q[[All, 1]] = b/Norm[b];
Do[
v = A.Q[[All, n]];
Do[
H[[j,n]] = Q[[All, j]].v;
v = v - H[[j,n]] Q[[All, j]],
{j, 1, n}
];
H[[n+1,n]] = Norm[v];
Q[[All, n+ 1]] = v/H[[n+1,n]],
{n, 1, nMax}
];
{Q,H}
]

TestMatrix=Import["/home/faculty/struther/Desktop/Link to 4610/TestMatricesETC/cavity03.rua"];
MatrixPlot[TestMatrix];
m=Length[TestMatrix];
nMax=12;
b=RandomReal[{-1,1},m];
{Qn,Hn}=Arnoldi[TestMatrix,b,nMax];
{QQn,HHn}=Arnoldi[TestMatrix,b,nMax+1];
Max[Abs[Hn-HHn[[1;;nMax+1,1;;nMax]]]]
Max[Abs[Qn-QQn[[All,1;;nMax+1]]]]

Eigenvalues

The idea is to run the Arnoldi iteration to generate a Krylov sapce of dimension n.  Compute eginevalues of the much smalle Hessenberg matroix (Subscript[H, n]) and hope that it approximates some of the eigen values of the larger matrix.  I will do this with one of the matrices from matrix market 
	http://math.nist.gov/MatrixMarket/.  
The idea works on general matrices.  We are interested in something sparse.

As a small example I downloaded a matrix cavityrua.org.  It is relatively small and non-symmetric.  It has somethin to do with a cavity driven flow!

TestMatrix=Import["/home/faculty/struther/Desktop/Link to 4610/TestMatricesETC/cavity03.rua"];
MatrixPlot[TestMatrix]

Mathematica's built in command takes a quarter of a second to find all the evals (using a dense technique).  It takes a lot less time to find the 10 largest.

Timing[evals=Eigenvalues[TestMatrix];]
Timing[Eigenvalues[TestMatrix,10];]

Lets make a picture

SpecPic=ListPlot[
Table[
{Re[evals[[i]]],Im[evals[[i]]]},
{i, Length[evals]}],
PlotRange->All,
Frame->True,
Axes->None,
GridLines->Automatic,
PlotLabel->"Spectra of cavity03.rua",
FrameLabel->{"real", "imag"}
]

Now we need to see what happens with the eigenvalues of the (Subscript[H, n]) from arnoldi iteration.

m=Length[TestMatrix];
b =RandomReal[{-1,1},{m}];
nMax=5;
{Q,H}=Arnoldi[TestMatrix,b,nMax];
MatrixForm[H]
MatrixForm[Drop[H,-1]]

nMax=124;
b =RandomReal[{-1,1},{m}];
{Q,H}=Arnoldi[TestMatrix,b,nMax];
ArnolidEvals=Eigenvalues[Drop[H,-1]];
ListPlot[{
Table[
{Re[evals[[i]]],Im[evals[[i]]]},
{i, Length[evals]}],
Table[
{Re[ArnolidEvals[[i]]],Im[ArnolidEvals[[i]]]},
{i, Length[ArnolidEvals]}]
},
PlotRange->All,
Frame->True,
Axes->None,
GridLines->Automatic,
PlotLabel->"Spectra of cavity03.rua",
FrameLabel->{"real", "imag"},
PlotStyle->{{PointSize[0.005],Red}, {PointSize[0.03],Opacity[0.3],Blue}}
]

GMRES

The idea is to run the Arnoldi iteration to generate a Krylov space of dimension n and then approximate the solution of 
	A.x=d
in the krylov space.  There are a couple of key observations.

From the theory of the Arnoldi iteration
	A.(Subscript[Q, n])=(Subscript[Q, n+1]).(Subscript[H, n])
The Qs are orthogonal and the columns of (Subscript[Q, n]) span the Krylov space (Subscript[𝒦, n]). Moreover, there is good reason to suspect that the Kryov spaces contain good (and improving as we increase n)  information about the eigen spaces associated with large eigen values of A.

The plan is to minimize the residual
	(Subscript[r, n])=(min_(x∈(Subscript[𝒦, n])))||A.x-b||
for larger and larger values of n untill the residual is sufficiently small.  The key to doing this efficiently is the basis for (Subscript[𝒦, n]) provided by the columns of (Subscript[Q, n]). The algorithm uses the invertible linear mapping y∈(ℝ^n)⟶x=(Subscript[Q, n]).y∈(Subscript[𝒦, n]) to provide an explicit representation of the least squares problem in the simpler space (ℝ^n). 
	||A.x-d||=||A.(Subscript[Q, n]).y-d||=||(Subscript[Q, n+1]).(Subscript[H, n]).y-d||=||(Subscript[H, n]).y-(Subscript[Q, n+1]^T).d||
If we take b=d the last term becomes cleaner since all the columns of (Subscript[Q, n+1]) (and hence rows of (Subscript[Q, n+1]^T)) except the first are orthogonal to b by construction.  So if b=d then we have 
	||A.x-d||=||(Subscript[H, n]).y-(Subscript[Q, n+1]^T).d||=||(Subscript[H, n]).y-||b||(Subscript[e, 1])||
The key to how well this works is how rapidly the Krylov subspaces provides good estimates and how well the iteration preserves orhtogonality.  Since
	(Subscript[𝒦, 0])⊂(Subscript[𝒦, 1])⊂(Subscript[𝒦, 2])⊂…⊂(Subscript[𝒦, m-1])⊂(Subscript[𝒦, m])=(ℝ^m)
the residuals (in exact arithmetic) satisfy
	(Subscript[r, 0])≥(Subscript[r, 1])≥(Subscript[r, 2])≥…≥(Subscript[r, m-1])≥(Subscript[r, m])=0
ans so we should have convergence.

A=Import["/home/faculty/struther/Desktop/Link to 4610/TestMatricesETC/fidap001.mtx"];
b=Import["/home/faculty/struther/Desktop/Link to 4610/TestMatricesETC/fidap001_rhs1.mtx"][[All,1]];
m=Length[A]; 
Timing[MmaSol=LinearSolve[A,b];]
ListPlot[MmaSol,
PlotRange->All]

Lets try the iteration!

nMax=183;
{Q,H}=Arnoldi[A,b,nMax];
Qn=Q[[All,1;;nMax]];
yn=LeastSquares[A.Qn, b];
Yn=LeastSquares[H,Transpose[Q].b];
xn =Qn.yn;
Xn=Qn.Yn;
ListPlot[ {xn,Xn,MmaSol},
PlotRange->All]
{((Norm[xn-Xn])/(Norm[MmaSol])),((Norm[xn-MmaSol])/(Norm[MmaSol]))}

These are hard test matrices.  We can make them easier by making them more diagonal.

ListPlot[Diagonal[A],
PlotRange->All
]

Of course, it is not the same matric but with stuff added to the diagonal it becomes much easier to solve accurately.

ANew=A +  IdentityMatrix[Length[b]];
MmaSol=LinearSolve[ANew,b];
nMax=23;
{Q,H}=Arnoldi[ANew,b,nMax];
Qn=Q[[All,1;;nMax]];
yn=LeastSquares[ANew.Qn, b];
Yn=LeastSquares[H,Transpose[Q].b];
xn =Qn.yn;
Xn=Qn.Yn;
ListPlot[ {xn,Xn,MmaSol},
PlotRange->All]
{((Norm[xn-Xn])/(Norm[MmaSol])),((Norm[xn-MmaSol])/(Norm[MmaSol]))}