Householder Reflectors

HouseU[Col_]:=Module[{x,γ,τ=Max[Abs[Col]]},
(* Watkins 3rd 3.2.35*)
(* returns γ and u *)
If[τ==0.0,Return[{0,Col}]];
x=Col/τ;
τ=Sign[x[[1]]]Norm[x];
x[[1]]+=τ;
γ=x[[1]]/τ;
x=x/x[[1]];
{γ,x}]

testing

n=34;
A0=RandomReal[{-10,10},{n,n}];
A=A0;
Id=IdentityMatrix[n];
i=2;
{γ,u}=HouseU[A[[All,i]]];
Q=IdentityMatrix[n]-γ Outer[Times,u,u];
(*A-=Outer[Times,γ u,Transpose[A].u];*)
MatrixPlot[Chop[Q.A]]

QR Decomposition

Sweeping through the submatrices starting from the top corner reduces A to an upper triangular matrix.

n=11;
A0=RandomReal[{-10,10},{n,n}];
b0=RandomReal[{-1,1},n];
{A,b}={A0,b0};

Id=IdentityMatrix[n];
Q=Id;
Pics=Table[{γ,u}=HouseU[A[[i;;-1,i]]];
A[[i;;-1,i;;-1]]-=Outer[Times,γ u,Transpose[A[[i;;-1,i;;-1]]].u];b[[i;;-1]]-=γ u (b[[i;;-1]].u);
{QOld,AOld}={Q,A};
Q[[i;;-1]]-=γ Outer[Times,u,Transpose[Q[[i;;-1]]].u];(*Plot QOld-Q or AOld-a to see the updating*)
MatrixPlot[Chop[A],PlotLabel->StringForm["i=``",i]],{i,1,n-1}]

We have Q.(Subscript[A, 0])=A=R.

Norm[Q.Transpose[Q]-Id]
Norm[Q.A0-A]

What is it good for?

LinearSolve[A0,b0]
LinearSolve[A,b]

Histogram[{Flatten[A],Flatten[A0]},n]

Hessenberg reduction

Similarity transformations leave eigenvalues unchanged i.e (Subscript[A, new]) and (Subscript[A, old]) have the same eigenvalues for  any invertible matrix H where

	(Subscript[A, new])=H.(Subscript[A, old]).(H^(-1))

It is easy to compute the eigenvalues of an upper triangular matrix.  

We can not reduce directly to an upper traingular matrix using Householder reflections because the multiplication by the (Q^T) destroys any zeros we just made!

A=RandomReal[{-10,10},{n,n}];

Id=IdentityMatrix[n];
Q=Id;
i=1;
{γ,u}=HouseU[A[[All,i]]];
Q=Id-Outer[Times,γ u,u];
{MatrixPlot[Chop[Q.A],PlotLabel->"Q.A"],MatrixPlot[Chop[Q.A.Transpose[Q]],PlotLabel->"Q.A.(Q^T)"]}
{Eigenvalues[A],Eigenvalues[Q.A.Transpose[Q]]}

Although it preserves eigenvalues it does not make zeros!  The Hessenberg trick is to be less greedy: instead of making all but the top of the first column zero we make all but the top two entries of the first column zero.

n=11;
A=RandomReal[{-10,10},{n,n}];

Id=IdentityMatrix[n];
i=1;
{γ,u}=HouseU[A[[2;;-1,i]]];
U=Join[{0},u]
Q=Id-Outer[Times,γ U,U];
{MatrixPlot[Chop[Q.A],PlotLabel->"Q.A"],MatrixPlot[Chop[Q.A.Transpose[Q]],PlotLabel->"Q.A.(Q^T)"]}
{Eigenvalues[A],Eigenvalues[Q.A.Transpose[Q]]}

Here is the sweep all the way through the matrix.

n=16;
A0=RandomReal[{-10,10},{n,n}];
(*A0=A0+Transpose[A0];*)
A=A0;

Id=IdentityMatrix[n];
Table[
{γ,u}=HouseU[A[[i+1;;-1,i]]];
U=Join[ConstantArray[0,i],u];
Q=Id-Outer[Times,γ U,U];
A=Chop[Q.A.Transpose[Q]];
MatrixPlot[A,PlotLabel->"Q.A.(Q^T)"],
{i,1,n-2}]
Eigenvalues[A]-Eigenvalues[A0]

This preserves the eigenvalues and does not break the zeros that were made by Q.

Symmetric Tridiagonal Eigenvalues:
Bulge Chasing

First a simple command to make a Hessenberg reduced test matrix!

HessenbergReducedSymmetric[n_]:=Module[{A=RandomReal[{-1,1},{n,n}],Id=IdentityMatrix[n],γ,u,i},A=A+Transpose[A];
Do[
{γ,u}=HouseU[A[[i+1;;-1,i]]];
U=Join[ConstantArray[0,i],u];
Q=Id-Outer[Times,γ U,U];
A=Chop[Q.A.Transpose[Q]],{i,1,n-2}];
A]

n=25;
Id=IdentityMatrix[n];
A0=HessenbergReducedSymmetric[n];
MatrixPlot[A0]

ListPlot[Abs[Map[Diagonal[A0,#]&,{0,-1}]],
PlotRange->All,
Joined->True,
Mesh->All,
PlotStyle->{Red, Blue}]

Now we try to make the first off diagonal entries zero and make a "bulge". Traditionally at the top.

A=A0;
i=1;
{γ,u}=HouseU[A[[i;;-1,i]]];
U=Join[ConstantArray[0,i-1],u];
Q=Id-Outer[Times,γ U,U];
A=Q.A.Transpose[Q];
MatrixPlot[A,PlotLabel->StringForm["i=``",i]]

And chase the bulge down through the matrix to the bottom.

Table[
{γ,u}=HouseU[A[[i+j;;i+j+1,j]]];
U=Join[ConstantArray[0,j],u,ConstantArray[0,n-j-2]];Q=Id-γ Outer[Times,U,U];
A=Chop[Q.A.Transpose[Q]];
MatrixPlot[A,
PlotLabel->StringForm["j=``",j]
],
{j,1,n-2}]
ListPlot[
Abs[Map[Diagonal[A,#]&,{0,-1}]],
PlotRange->All,
Joined->True,
Mesh->All
]

This got better but we can make it get better faster!  How can we tell how nearly diagonal each bit is?

ListLogPlot[
Abs[Diagonal[A,1]/Drop[Diagonal[A,0],-1]]
]

HessenbergDiagonalRatio[A_]:=Module[{},
Abs[Diagonal[A,-1]/Drop[Diagonal[A,0],-1]]
]
CreateBulge[AIn_,ρ_]:=Module[
{A=AIn,i=1,γ,v,u,U,Q,Id =IdentityMatrix[Length[AIn]]},
v=A[[i;;-1,i]];
v[[1]]-=ρ;
{γ,u}=HouseU[v];
U=Join[ConstantArray[0,i-1],u];
Q=Id-Outer[Times,γ U,U];
A=Q.A.Transpose[Q];
A]
ChaseBulge[AIn_]:=Module[
{A=AIn,γ,u,U,Q,j,Id =IdentityMatrix[Length[AIn]],i=1},
Do[
{γ,u}=HouseU[A[[i+j;;i+j+1,j]]];
U=Join[ConstantArray[0,j],u,ConstantArray[0,n-j-2]];
Q=Id-γ Outer[Times,U,U];
A=Chop[Q.A.Transpose[Q]],
{j,1,n-2}
];
A
]
FrancisIteration[AIn_,ρ_]:=Module[
{A=AIn},
A=CreateBulge[A,ρ];
ChaseBulge[A]
]

Testing these

n=25;
Id=IdentityMatrix[n];
A0=HessenbergReducedSymmetric[n];
A=A0;
A=CreateBulge[A,0];
A=ChaseBulge[A];
Eigenvalues[A]-Eigenvalues[A0]
ListLogPlot[
Map[HessenbergDiagonalRatio,{A,A0}],
Joined->True,
Mesh->All,
PlotStyle->{Red, Blue}]

Lets try to understand the convergence.

BarChart[Abs[Transpose[{Drop[Diagonal[A,0],-1],Diagonal[A,-1]}]]]

The large entries in the bulge are always at the corners. Does this need to be this way?

Behavior of Shifted and unshifted algorithm.

How well does this work.

WilkinsonShift[ASub_]:=Sort[Eigenvalues[ASub],Abs[#1-A[[-1,-1]]]<Abs[#2-A[[-1,-1]]]&][[1]]
WilkinsonShift[{{2,1},{1,1.8}}]

A=A0;
MaxIter=20;
ListContourPlot[
Table[
ρ=0; (* unshifted *)
A=FrancisIteration[A,ρ];
Log[10,HessenbergDiagonalRatio[A]],
{MaxIter}],
ContourLabels->All,
FrameLabel->{"mat","iter"},
PlotRange->All
]

A=A0;
MaxIter=20;
A=A0;
MaxIter=20;
ListContourPlot[
Table[
ρ=A[[-1,-1]]; (* Rayleigh unshifted *)
A=FrancisIteration[A,ρ];
Log[10,HessenbergDiagonalRatio[A]],
{MaxIter}],
ContourLabels->All,
FrameLabel->{"mat","iter"},
PlotRange->All
]

The funny stuff in the corner is because the bottom subdiagonal is zero after a couple of steps.

A=A0;
MaxIter=20;
A=A0;
MaxIter=20;
ListContourPlot[
Table[
ρ=WilkinsonShift[A[[-2;;-1,-2;;-1]]]; (* Wilkinson shift *)
A=FrancisIteration[A,ρ];
Log[10,HessenbergDiagonalRatio[A]],
{MaxIter}],
ContourLabels->All,
FrameLabel->{"mat","iter"},
PlotRange->All
]

The funny stuff in the corner is because the bottom subdiagonal is zero after a couple of steps.

MatrixPlot[A]

Once the entry is small enough we should "deflate" the problem and restart with a smaller matrix.

Parallelization

Multiple bulges can run down the matrix together. This is a standard process.  I do not think this lends itself to the GPU.

Bulges can be bigger and for CUDA bigger is better. How much bigger?  We would like to get up to 16 or so

LAPACK implements m=6 for large eigenvalue problems bigger than 50×50.

Watkins "Evalues" has a discussion about error, initiation and propagation in larger bulges and recommends no larger than 8.

Several questions about the numerical experiments underlying this recomendation:

Is it general for any initial reduction?  Or is it specific for starting with a Hessenberg reduced from a dense.  As we saw this results in large sub-diagonal entries.  Is it different for a more accurate "jacobi" initial reduction? How is this precision dependent?

If bulges can only have m=4 or 8 can we do anything with AD derivatives wrt the shifts ρ.

Possibly attempt to improve accuracy.

Possibly attempt to improve convergence.

Possibly attempt to create deflation in the middle of the matrix.

Possibly attempt to create deflation say 16 up from the bottom of the matrix.

Reduction to Hessenberg form is time consuming it seems counter intuitive to reduce to one sub-diagonal and then immediately start computing with more subdiagonals.  Can one expose more parallelization in the initial reduction if your target is to get to say a 8 subdiagonal hessenberg matrix.  If one does does the first step of the bulge chase from this form do you avoid the initial loss of accuracy that plagues high-order bulge chasing algorithms.