Eigenvalue problem - simple iterations

Preliminaries

Clear all symbols from previous evaluations to avoid conflicts

Clear["Global`*"]

Problem

[Taken from Trefethen, Bau: Numerical Linear Algebra, SIAM 1997, p. 315]

n=1000;
A=SparseArray[{{i_,i_}->0.5+Sqrt[i],{i_,j_}/;Abs[i-j]==1->1.0,{i_,j_}/;Abs[i-j]==100->1.0},{n,n}];
b=ConstantArray[1.0,n];

The matrix is symmetric positive definite - all eigenvalues are positive

eigenA=Sort[Eigenvalues[A]];
maxλ=Max[eigenA]; minλ=Min[eigenA];
Print["κ(A) = ",Max[eigenA]/Min[eigenA]]
Print["max(λ) = ",maxλ]
Print["min(λ) = ",minλ]
ListPlot[{eigenA},PlotRange->All,AxesLabel->{"n","eigenvalue"}]

Maximal eigenvalue

Power iteration - quite slow convergence

niter=2000;
error=ConstantArray[0.0,niter];
x=b/Norm[b];
Do[
w=A.x;
λ=x.w;
error[[i]]=Abs[maxλ-λ];
x=w/Norm[w],
{i,1,niter}
];
Print["Exact maximum   λ: ",maxλ]
Print["Power iteration λ: ",λ, ", error after ", niter," iterations: ",error[[niter]]]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Inverse iteration - reasonable convergence if we set the shift μ close to actual eigenvalue

niter=10;
error=ConstantArray[0.0,niter];
x=b/Norm[b];
μ=34.5;
Do[
w=LinearSolve[A-μ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[i]]=Abs[maxλ-λ],
{i,1,niter}
];
Print["Exact maximal λ: ",maxλ]
Print["Inverse iteration λ: ",λ, ", error after ", niter," iterations: ",error[[niter]]]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Rayleigh quotient iteration - fast convergence  if we set initial guess of λ close to actual eigenvalue

niter=3;
error=ConstantArray[0.0,niter];
x=b/Norm[b];
(* x=xInvIter; *)
λ=34.5;
Do[
w=LinearSolve[A-λ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[i]]=Abs[maxλ-λ],
{i,1,niter}
];
Print["Exact maximal λ: ",maxλ,", Rayleigh quotient iteration λ: ",λ]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Minimal eigenvalue

Inverse iteration - again quite fast convergence

niter=10;
error=ConstantArray[0.0,niter];
x=b/Norm[b];
μ=0.0;
Do[
w=LinearSolve[A-μ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[i]]=Abs[minλ-λ],
{i,1,niter}
];
Print["Exact minimal λ: ",minλ,", inverse iteration λ: ",λ]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Rayleigh quotient iteration - fast convergence but to an incorrect eigenvalue 12.57265157916

niter=4;
error=ConstantArray[0.0,niter];
x=b/Norm[b]; (* not a good initial guess for the minimal eigenvalue *)
λ=0.0;
Do[
w=LinearSolve[A-λ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[i]]=Abs[Min[eigenA]-λ],
{i,1,niter}
];
Print["Exact minimal λ: ",minλ,", inverse iteration λ: ",λ]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Rayleigh quotient iteration - fast convergence if we get close enough using inverse iteration first

niterInverse=2;
niterRayleigh=3;
error=ConstantArray[0.0,niterInverse+niterRayleigh];
(* inverse iteration *)
x=b/Norm[b];
μ=0.0;
Do[
w=LinearSolve[A-μ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[i]]=Abs[Min[eigenA]-λ],
{i,1,niterInverse}
];
(* Rayleigh quotient iteration *)
Do[
w=LinearSolve[A-λ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
error[[niterInverse+i]]=Abs[Min[eigenA]-λ],
{i,1,niterRayleigh}
];
Print["Exact minimal λ: ",minλ,", inverse iteration λ: ",λ]
ListLogPlot[{error},PlotRange->All,AxesLabel->{"niter","error"}]

Other eigenvalues

eigenA[[380;;400]]

Inverse iteration - finds an eigenvalue closest to the initial quess

x=b/Norm[b];
μ=20.0; (* initial guess *)
Do[
w=LinearSolve[A-μ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
Print[λ],
{i,1,10}
];

Rayleigh quotient iteration

λ=20.0;
x=b/Norm[b];
Do[
w=LinearSolve[A-λ IdentityMatrix[n],x];
x=w/Norm[w];
λ=x.(A.x);
Print[λ],
{i,1,4}
];

Basic QR algorithm

n=300;
A=SparseArray[{{i_,i_}->0.5+Sqrt[i],{i_,j_}/;Abs[i-j]==1->1.0,
{i_,j_}/;Abs[i-j]==10->1.0},{n,n}];
(* A=SparseArray[{{i_,i_}->0.5+Sqrt[i],{i_,j_}/;Abs[i-j]==1->1.0},{n,n}]; *)
niter=3000;

AA=A;
SumOffDiagQR=ConstantArray[0.0,niter+1];
SumOfLastRowQR=ConstantArray[0.0,niter+1];
SumOffDiagQR[[1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowQR[[1]]=Total[Abs[AA[[n,1;;n-1]]]^2];
Do[
{Q,R}=QRDecomposition[AA];
AA=R.ConjugateTranspose[Q];
SumOffDiagQR[[i+1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowQR[[i+1]]=Total[Abs[AA[[n,1;;n-1]]]^2],
{i,1,niter}
];
ListLogPlot[SumOffDiagQR,AxesLabel->{"niter","sum of off-diagonal elements"}]
Norm[Sort[Diagonal[AA]]-Sort[Eigenvalues[A]]]

ListLogPlot[SumOfLastRowQR[[1;;100]],AxesLabel->{"niter","sum of last row"},PlotRange->{10^-30,10}]

QR algorithm with shifts

Rayleigh quotient shift

AA=A;
SumOffDiagRay=ConstantArray[0.0,niter+1];
SumOfLastRowRay=ConstantArray[0.0,niter+1];
SumOffDiagRay[[1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowRay[[1]]=Total[Abs[AA[[n,1;;n-1]]]^2];
μ=AA[[n,n]];
Do[
{Q,R}=QRDecomposition[AA-μ IdentityMatrix[n]];
AA=R.ConjugateTranspose[Q]+μ IdentityMatrix[n];
μ=AA[[n,n]];
SumOffDiagRay[[i+1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowRay[[i+1]]=Total[Abs[AA[[n,1;;n-1]]]^2],
{i,1,niter}
];
ListLogPlot[{SumOffDiagQR,SumOffDiagRay},AxesLabel->{"niter","sum off-diag"},PlotLegends->{"Basic","Rayleigh"}]
Norm[Sort[Diagonal[AA]]-Sort[Eigenvalues[A]]]

ListLogPlot[{SumOfLastRowQR[[1;;40]],SumOfLastRowRay[[1;;20]]},AxesLabel->{"niter","sum of last row"},PlotRange->{10^-40,10}]

Wilkinson shift

AA=A;
SumOffDiagW=ConstantArray[0.0,niter+1];
SumOfLastRowW=ConstantArray[0.0,niter+1];
SumOffDiagW[[1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowW[[1]]=Total[Abs[AA[[n,1;;n-1]]]^2];
Do[
δ=(AA[[n-1,n-1]]-AA[[n,n]])/2;
μ=AA[[n,n]]-Sign[δ]*AA[[n,n-1]]^2/(Abs[δ]+Sqrt[δ^2+AA[[n,n-1]]^2]);
{Q,R}=QRDecomposition[AA-μ IdentityMatrix[n]];
AA=R.ConjugateTranspose[Q]+μ IdentityMatrix[n];
SumOffDiagW[[i+1]]=Total[Abs[LowerTriangularize[AA,-1]]^2,2];
SumOfLastRowW[[i+1]]=Total[Abs[AA[[n,1;;n-1]]]^2],
{i,1,niter}
];
ListLogPlot[{SumOffDiagQR,SumOffDiagRay,SumOffDiagW},AxesLabel->{"niter","sum off-diag"},PlotLegends->{"Basic","Rayleigh","Wilkinson"}]
Norm[Sort[Diagonal[AA]]-Sort[Eigenvalues[A]]]

ListLogPlot[{SumOfLastRowQR[[1;;40]],SumOfLastRowRay[[1;;20]],SumOfLastRowW[[1;;20]]},AxesLabel->{"niter","sum of last row"},PlotRange->{10^-40,10}]