Reference

“Block Householder transformation for parallel QR factorisation”, Applied Mathematics Letters 12, pp.29-34, 1999, F. Rotella, I. Zambettakis

Trial Implementation

First a blocked definition of a Householder reflector.

BlockHouse[V_,0]:= Module[
{n=Length[V]},
SparseArray[Band[{1,1}]->1,{n,n}]-2V.Inverse[Transpose[V].V].Transpose[V]
]

BlockHouse[V_,1]:= Module[
{n=Length[V]},
SparseArray[Band[{1,1}]->1,{n,n}]-2V.PseudoInverse[V]
]

Testing for orthogonality.

{m,r}={323,24};
V=RandomReal[{-1,1},{m,r}];
H= BlockHouse[V,0];
Norm[IdentityMatrix[m]-H.Transpose[H]]

Speed trials. The zero version seems to sin.

{m,r}={323,126};
V=RandomReal[{-1,1},{m,r}];
Timing[BlockHouse[V,0];]
Timing[BlockHouse[V,1];]

Now a Discussion of how to kill the bottom bit of some columns.

HouseholderV[A1A2_,0]:= Module[
{A1,A2, m,r,P,λ,X},
{m,r}=Dimensions[A1A2];
A1=Take[A1A2,r];A2=Drop[A1A2,r];
P=A2.Inverse[A1]; POut=P;
{λ,P}=Eigensystem[IdentityMatrix[r]+ Transpose[P].P];
X=Transpose[P].DiagonalMatrix[(Sqrt[λ])].P.A1;
ArrayFlatten[{{A1+X},{A2}}]
]

HouseholderV[A1A2_,1]:= Module[
{A1,A2, m,r,P,λ,X},
{m,r}=Dimensions[A1A2];
A1=Take[A1A2,r];A2=Drop[A1A2,r];
P=A2.Inverse[A1];
X=MatrixPower[IdentityMatrix[r]+ Transpose[P].P,1/2].A1;
ArrayFlatten[{{A1+X},{A2}}]
]

Testing the zeroing property. Looks like both work.

{m,r}={25,12};
A1A2=RandomReal[{-1,1},{m,r}];
V=HouseholderV[A1A2,0];
H=BlockHouse[V,0];
Norm[ IdentityMatrix[m]-H.Transpose[H]]
MatrixPlot[Chop[H.A1A2]]

Time trials.  Looks as though the 1th version is faster!

{m,r}={223,24};
A1A2=RandomReal[{-1,1},{m,r}];
{Timing[HouseholderV[A1A2,0];] ,Timing[HouseholderV[A1A2,1];]}
{m,r}={223,56};
A1A2=RandomReal[{-1,1},{m,r}];
{Timing[HouseholderV[A1A2,0];] ,Timing[HouseholderV[A1A2,1];]}

Note that for r=1 this is the same as the Householder guy we already looks at

{m,r}={25,1};
A1A2=RandomReal[{-1,1},{m,r}];
V=HouseholderV[A1A2,0];
H=BlockHouse[V,0];
Norm[ IdentityMatrix[m]-H.Transpose[H]]
MatrixPlot[Chop[H.A1A2]]

QR Decomposition

The goal of this is to use orthogonal matrices to reduce a matrix to triangular. Here is a quick review of the block version of the QR decomposition.

{NumberOfBlocks,BlockSize}={5,8};
n=NumberOfBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;

U=A[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
A=Chop[H.A];
Q=H;
MatrixPlot[A,
PlotLabel->{Norm[Q.A-AOriginal],Norm[IdentityMatrix[n]-Q.Transpose[Q]]}
]

i=2;
ASub=A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]];
U=ASub[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
HBig=ArrayFlatten[{
{IdentityMatrix[(i-1)*BlockSize],0},
{0, H}
}];
Q=Q.HBig;
ASub=Chop[H.ASub];
A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]]=ASub;
MatrixPlot[A,
PlotLabel->{Norm[Q.A-AOriginal],Norm[IdentityMatrix[n]-Q.Transpose[Q]]}
]

i=3;
ASub=A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]];
U=ASub[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
HBig=ArrayFlatten[{
{IdentityMatrix[(i-1)*BlockSize],0},
{0, H}
}];
Q=Q.HBig;
ASub=Chop[H.ASub];
A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]]=ASub;
MatrixPlot[A,
PlotLabel->{Norm[Q.A-AOriginal],Norm[IdentityMatrix[n]-Q.Transpose[Q]]}
]

Making a loop. No thought process here about efficiency.

BlockHouseQRreduce[AIn_,BlockSize_]:= Module[
{U,V,H,A,Q,HBig,n=Length[AIn]},

U=AIn[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
A=Chop[H.AIn];
Q=H;

Do[
ASub=A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]];
U=ASub[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
HBig=ArrayFlatten[{
{IdentityMatrix[(i-1)*BlockSize],0},
{0, H}
}];
Q=Q.HBig;
ASub=Chop[H.ASub];
A[[1+(i-1)*BlockSize;;n,1+(i-1)*BlockSize;;n]]=ASub,
{i,2,Floor[n/BlockSize]-1}];

{Q,A}
]

Testing the code.

BlockSize=5;
BlockCount=5;
n=BlockSize*BlockCount;
A=RandomReal[{-1,1},{n,n}];
{Q,R}=BlockHouseQRreduce[A,BlockSize];
MatrixPlot[R,
PlotLabel->{
Norm[IdentityMatrix[n]-Q.Transpose[Q]],
Norm[A-Q.R]}
]

Hessenberg Reduction

The goal of this is to use orthogonal matrices to reduce a matrix to almost triangular form using similarity transforms. We can not do this exactly the same way we did the QR decomposition because post multiplying by Transpose[Q] trashes the zeros we just made.

{NumberOfBlocks,BlockSize}={5,8};
n=NumberOfBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;

U=A[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
A=Chop[H.A.Transpose[H]];
MatrixPlot[A,
PlotLabel->Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]
]

The trick is to get a little bit less greedy and leave the top corner alone and make some zeroes further down.  Here is an example where we start further down by the block size.

{NumberOfBlocks,BlockSize}={5,8};
n=NumberOfBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;

U=A[[1+BlockSize;;n,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
HBig=ArrayFlatten[{
{IdentityMatrix[BlockSize],0},
{0,H}}];
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,
PlotLabel->Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]
]

There are other choices. Here is an example that shifts down by two rows. As a result it only preserves two of the columns of zeros.

{NumberOfBlocks,BlockSize}={5,5};
Shift=2;
n=NumberOfBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;

U=A[[1+Shift;;n,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
HBig=ArrayFlatten[{
{IdentityMatrix[Shift],0},
{0,H}}];
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,
PlotLabel->Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]
]

Built In Hessenberg Reduction

The built-in Mathematica Hessenberg decompostion does the non-blocked version.

n=23;
A= RandomReal[{-1,1},{n,n}];
(*A=A.Transpose[A];*)
{Q,H}=HessenbergDecomposition[A];
MatrixPlot[H,
PlotLabel->{Norm[IdentityMatrix[n]-Q.Transpose[Q]],Norm[Transpose[Q].A.Q-H]}]

Discussion

One thing that is clear about the non-blocked version is that it produces large sub-diagonal entries. One thing that is less obvious is that Rotations do not change the “Frobenius” norm of a matrix.

A= RandomReal[{-1,1},{n,n}];
(*A=A.Transpose[A];*)
{Q,H}=HessenbergDecomposition[A];
Map[Norm[#,"Frobenius"]&,{A,H}]

We can compute a scaled measure of how upper triangular a matrix is by the ratio of the Frobenius norm of the bottom triangular bit (not including the diagonal) to the Frobenius norm of the whole thing. We can mess around with the Squares and Square roots that are in the defn of the Frobenius norm.

n=1230;
A=RandomReal[{-1,1},{n,n}];
MeasureUpperT[A_]:= (Sqrt[(((Norm[LowerTriangularize[A,-1],"Frobenius"])^2)/((Norm[A,"Frobenius"])^2))])
{MeasureUpperT[A],MeasureUpperT[H]}