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}={123,24};
A1A2=RandomReal[{-1,1},{m,r}];
V=HouseholderV[A1A2,0];
H=BlockHouse[V,0];
{MatrixPlot[Chop[H.A1A2]],MatrixPlot[Chop[Transpose[H].H]]}

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];]}

QR Decomposition

n=16*34;
A=RandomReal[{-1,1},{n,n}];
U=A[[All,1;;16]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
A=Chop[H.A];
MatrixPlot[A]

ASub=A[[17;;-1,17;;-1]];
U=ASub[[All,1;;16]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
ASub=Chop[H.ASub];
A[[17;;-1,17;;-1]]=ASub;
MatrixPlot[A]

Making a loop. No thought process here about efficiency.

BlockHouseQRreduce[AIn_,BlockSize_]:= Module[
{A=AIn},
Do[
ASub=A[[
1+i*BlockSize;;-1,
1+i*BlockSize;;-1]];
U=ASub[[All,1;;BlockSize]];
V=HouseholderV[U,1];
H=BlockHouse[V,0];
ASub=Chop[H.ASub];
A[[
1+i*BlockSize;;-1,
1+i*BlockSize;;-1]]=ASub,
{i,0, Floor[((Length[A])/BlockSize)]-2}
];
A
]

Testing the code.

BlockSize=5;
BlockCount=15;
n=BlockSize*BlockCount;
A=RandomReal[{-1,1},{n,n}];
AOut=BlockHouseQRreduce[A,BlockSize];
MatrixPlot[AOut]

So it is block uper triangular. How am I going to test that I have not messed up! well I need to go the last few steps and compare the result to a standard QR decomposition.  It is slightly messier than you think because there are ± sign possibilities!

{Q,R}=QRDecomposition[A];
{QOut,ROut}=QRDecomposition[AOut];
Map[MatrixPlot,Chop[{R-ROut,R+ROut}]]

The zero stripes match up just like they look they do. I made a function to check

SignIndifferentMatrixDiff[R, ROut]

Looks as though they do.  If I was going to do anything real with this I would make a function.

SignIndifferentMatrixDiff[A1_,A2_]:= Module[
{},
Max[
Map[Min,
Transpose[Map[Norm,{A1+A2,A1-A2},{2}]]
]
]
]

Hessenberg Reduction

BlockSize=4;
NumBlocks=6;
n=NumBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
AOriginal=AOriginal.Transpose[AOriginal];
A=AOriginal;
Do[
U=A[[All,1+(i-1)*BlockSize;;i*BlockSize]];
ZeroBlock=0*U[[1;;i*BlockSize]];
U=Drop[U, i*BlockSize];
V=HouseholderV[U,1];
V=Join[ZeroBlock,V];
H=BlockHouse[V,0];
A=Chop[H.A.Transpose[H]];
MatrixPlot[A],
{i,1, NumBlocks-2}]

{QHBI,HBI}=HessenbergDecomposition[AOriginal];
{QHBI2,HBI2}=HessenbergDecomposition[A];
{MatrixPlot[HBI] ,MatrixPlot[HBI2]} 
Map[ ((Norm[LowerTriangularize[#,-1]])/(Norm[#,"Frobenius"]))&,{HBI,HBI2}]

Map[ Norm[UpperTriangularize[#,1],"Frobenius"]&,{HBI,HBI2}]

ListPlot[ {Abs[Diagonal[HBI,-1]],Abs[Diagonal[HBI2,-1]]},
PlotRange->All,
PlotStyle->{{PointSize[0.02],Red}, Blue}]

ListLogPlot[Abs[Eigenvalues[A]-Eigenvalues[AOriginal]],
Frame->True, PlotRange->All,GridLines->Automatic]
ListLogPlot[((Abs[Eigenvalues[A]-Eigenvalues[AOriginal]])/(Abs[Eigenvalues[A]])),
Frame->True, PlotRange->All,GridLines->Automatic]

ListLogPlot[{
Sort[Abs[Diagonal[AOriginal]]],
Sort[Abs[Diagonal[AOriginal,12]]],
Sort[Abs[Diagonal[A]]]
},
GridLines->Automatic,
Frame->True,
PlotRange->All,
PlotStyle->{Red,Green,Blue}
]

Hessenberg Reduction: Fancy

BlockSize=8;
NumBlocks=6;
n=NumBlocks*BlockSize;
AOriginal=RandomReal[{-1,1},{n,n}];
(*AOriginal=AOriginal+Transpose[AOriginal];*)
A=AOriginal;
Table[
U=A[[All,i;;i+BlockSize]];
ZeroBlock=0*U[[1;;i+BlockSize]];
U=Drop[U, i+BlockSize];
V=HouseholderV[U,0];
V=Join[ZeroBlock,V];
H=BlockHouse[V,0];
A=Chop[H.A.Transpose[H]];
MatrixPlot[A],
{i,1,12,1}]

ListPlot3D[
Abs[A],
PlotRange->All
]

A=Chop[H.A.Transpose[H]];
MatrixPlot[A]

{QHBI,HBI}=HessenbergDecomposition[AOriginal];
{QHBI2,HBI2}=HessenbergDecomposition[A];
{MatrixPlot[HBI] ,MatrixPlot[HBI2]} 
Map[ ((Norm[LowerTriangularize[#,-1]])/(Norm[#,"Frobenius"]))&,{HBI,HBI2}]

Map[ Norm[UpperTriangularize[#,2],"Frobenius"]&,{HBI,HBI2}]

ListPlot[ {Abs[Diagonal[HBI,-1]],Abs[Diagonal[HBI2,-1]]},
PlotRange->All,
PlotStyle->{{PointSize[0.02],Red}, Blue}]

ListLogPlot[Abs[Eigenvalues[A]-Eigenvalues[AOriginal]],
Frame->True, PlotRange->All,GridLines->Automatic]
ListLogPlot[((Abs[Eigenvalues[A]-Eigenvalues[AOriginal]])/(Abs[Eigenvalues[A]])),
Frame->True, PlotRange->All,GridLines->Automatic]

ListLogPlot[{
Sort[Abs[Diagonal[AOriginal]]],
Sort[Abs[Diagonal[AOriginal,12]]],
Sort[Abs[Diagonal[A]]]
},
GridLines->Automatic,
Frame->True,
PlotRange->All,
PlotStyle->{Red,Green,Blue}
]