Non-Sequitor: MachineEpsilon

So I can find out what Mathematica (on my machine) thinks machine expsilon is.

({{, $MachineEpsilon
gives the difference between 1.0 and the next-nearest number representable as a machine-precision number. }})

$MachineEpsilon
x=1.0 +$MachineEpsilon
x-1.0

explain the foillowing plot

Plot[(1+x)-1,
{x,-3$MachineEpsilon,3 $MachineEpsilon},PlotRange->All]

I can test for machine esilon by using the definition.  Machine ϵ is the largest number for which 
	1.0+ϵ=1.0
Note:  This depends on exactly how you test (matlab and mathematica understand that floating point tests are funky)
and how you change the numbers.

ϵ=0.1;
While[
(* Print[ϵ];*)
1.0+ϵ>1.0,
ϵ=0.9 ϵ
]
ϵ
1.0 +ϵ>1.0
1.0 +(1/0.9)ϵ>1.0

ϵ=1.0 (10^(-20));
While[
(* Print[ϵ];*)
1.0+ϵ==1.0,
ϵ=1.1 ϵ
]
ϵ
1.0+ϵ==1.0
1.0+(1/1.1)ϵ==1.0

Under and overflow are separate ideas.

$MinMachineNumber
$MaxMachineNumber

Loss of precision

Clear[f,g,x]
f[x_]:= 1/(x^2) ;
g[x_]:=1/(x^2) -1/(x^3);
fTakeAwayg[x_]=f[x]-g[x]
LogLogPlot[{
f[x]-g[x],
fTakeAwayg[x]},
{x, 1,(10^152)},
PlotStyle->{{Thickness[0.02],Red}, Blue},
Frame->True,
GridLines->Automatic]

Non-Sequitor: 0.1 in binary

RealDigits[x]
gives a list of the digits in the approximate real number x, together with the number of digits that are to the left of the decimal point. 

RealDigits[x,b]
gives a list of base-b digits in x.

{digits,power}=RealDigits[0.1,2]

The general pattern is 1100110011001100110011001100…  The approximation to 0.1 is 
	(1+1/2+0/4+0/8+1/16+1/32+0/64+0/128+…+ 0/(2^47)+1/(2^48)+1/(2^49)+0/(2^50)+1/(2^51)+0/(2^53))/(2^(3+1))

0.1`123 -Sum[ 1/(2^(4 i))+1/(2^(4 i + 1)),{i, 0, 10}]/(2^(3+1))
0.1`123 -Sum[ 1/(2^(4 i))+1/(2^(4 i + 1)),{i, 0, 11}]/(2^(3+1))
0.1`123 -Sum[ 1/(2^(4 i))+1/(2^(4 i + 1)),{i, 0, 12}]/(2^(3+1))
0.1`123 -Sum[ 1/(2^(4 i))+1/(2^(4 i + 1)),{i, 0, 13}]/(2^(3+1))
0.1`123 -Sum[ 1/(2^(4 i))+1/(2^(4 i + 1)),{i, 0, 14}]/(2^(3+1))

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

Timings: HouseholderV

A perfectly good question is “How does this scale with m and r?”

m=800;
RunSize=10;
DataFixedm=Table[
r=Floor[20*(1.5^p)];
A1A2=RandomReal[{-1,1},{m,r}];
Data=Table[Timing[HouseholderV[A1A2,1];][[1]],{RunSize}];
{r,Mean[Data]},
{p,0,9,1}]

ListLogLogPlot[{
DataFixedm,
{{50, (0.01/(50^2)) (50^2)},{700, (0.01/(50^2)) (700^2)}}
},
Frame->True,
Joined->True,
Mesh->All,
GridLines->Automatic,
FrameLabel->{"r","Sec"},
PlotLabel->StringForm["m=``",m]
]

r=16;
RunSize=10;
DataFixedr=Table[
m=Floor[20*(1.5^p)];
A1A2=RandomReal[{-1,1},{m,r}];
Data=Table[Timing[HouseholderV[A1A2,1];][[1]],{RunSize}];
{m,Mean[Data]},
{p,3,20,1}]

ListLogLogPlot[{
DataFixedr,
{{50, (0.0001/50) 50},{70000, (0.0001/50) 70000}}
},
Frame->True,
Joined->True,
Mesh->All,
GridLines->Automatic,
FrameLabel->{"m","Sec"},
PlotLabel->StringForm["r=``",r]
]

So as a guess.  Forming the householder “V” reducer  is linear in m and quadratic in r.

Timings: BlockHouse

A perfectly good question is “How does this scale with m and r?”

m=2000;
RunSize=5;
DataFixedm=Table[
r=Floor[20*(1.5^p)];
A1A2=RandomReal[{-1,1},{m,r}];
Data=Table[Timing[BlockHouse[A1A2,0];][[1]],{RunSize}];
{r,Mean[Data]},
{p,0,11,1}]
ListLogLogPlot[{
DataFixedm,
{{400, (0.5/400) 400},{1700, (0.5/400) 1700}}
},
Frame->True,
Joined->True,
Mesh->All,
GridLines->Automatic,
FrameLabel->{"r","Sec"},
PlotLabel->StringForm["m=``",m]
]

Not sure about you but this may not look like a straight line on the log-log plot.  On the oher hand I guess the one below does.

r=16;
RunSize=10;
DataFixedr=Table[
m=Floor[20*(1.5^p)];
A1A2=RandomReal[{-1,1},{m,r}];
Data=Table[Timing[BlockHouse[A1A2,0];][[1]],{RunSize}];
{m,Mean[Data]},
{p,3,12,1}]

ListLogLogPlot[{
DataFixedr,
{{50, (0.0001/(50^2)) (50^2)},{7000, (0.0001/(50^2)) (7000^2)}}
},
Frame->True,
Joined->True,
Mesh->All,
GridLines->Automatic,
FrameLabel->{"m","Sec"},
PlotLabel->StringForm["r=``",r]
]

So as a guess.  Forming the householder “V” reducer  is linear in m and quadratic in r.

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

Lets see if the zero stripes match.

Map[Norm, R-ROut]
Map[Norm, 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}]]
]
]
]
SignIndifferentMatrixDiff[R, ROut]

QR Pivoting

By default the built-in QR decomposition does not pivot.  It does tend to reveal the pivots in decreasing order.  Turning on the pivoting makes it run slower.

n=345
A=RandomReal[{-1,1},{n,n}];
Timing[{Q,R}=QRDecomposition[A];]
Timing[{Qp,Rp,P}=QRDecomposition[A,
Pivoting->True];]
Map[ListPlot[Abs[Diagonal[#]]]&, {R,Rp}]

The matrix P is one of the permutation matrices. It is automatically stored as a sparse array and does the same thing as for the LU decomposition.

Head[P]
Norm[A.P-Transpose[Qp].Rp]
MatrixPlot[P]

Usually you do not need to pivot in a QR unless the matrix is deficient.

Gram -Schmidt

Thm: 3.4.2:   A matrix is orthogonal iff its cols (rows) are orthonormal.

Def: 3.4.2:   A matrix Q is an isometry if  its cols are orthonormal.

Thm: 3.4.2:   A matrix Q∈(ℝ^(n×m)) (n>m) is an isometry iff (Q^T).Q=Id∈(ℝ^(m×m))

Someone gives us a set of vectors
	{(Subscript[v, 1]),(Subscript[v, 2]),…,(Subscript[v, m])}⊂(ℝ^n)
Lets try to think how we might go about building an orthonormal set 
	{(Subscript[q, 1]),(Subscript[q, 2]),…,(Subscript[q, m])}⊂(ℝ^n)
that spans the same subspace of (ℝ^n).  We can easily build a set so that at each stage we 
span the same set.

The first coilumn just needs scaled.

{n,m}={3,4};
V=RandomReal[{-1,1},{m,n}];
(*MatrixForm[A]*)
q1=V[[All,1]];
q1=q1/Norm[q1]

The second column needs to be made orthogonal to the first and scaled.

q2=V[[All,2]];
q2=q2-(q1.q2)q1;
q2=q2/Norm[q2]

The third column needs to be made orthogonal to the first and second then scaled.

q3=V[[All,3]];
q3=q3-(q1.q3)q1;
q3=q3-(q2.q3)q2;
q3=q3/Norm[q3]