Classical

A direct implementation of algorithm 7.1

CGS[A_]:= Module[
{m,n,Q,R,i,j,v,r},
{m,n}=Dimensions[A];
Q=Table[0,{m},{n}];
R=Table[0,{n},{n}];
Do[
v=A[[All,j]];
Do[
r=Q[[All,i]].A[[All,j]];
v=v-r Q[[All,i]];
R[[i,j]]=r,
{i,j-1}];
r=Norm[v];
Q[[All,j]]=v/r;
R[[j,j]]=r,
{j,1,n}
];
{Q,R}
]

Testing the code.  The output Q should be approximately orthogonal and have the same column space as the input. It is easy to check small test problems for Orthogonality and that the code can easily genearte an R matrix

Testing: Small

{m,n}={9,9};
A= RandomReal[{-1,1},{m,n}];
{Q,R}=CGS[A];
RT= Transpose[Q].A;
MatrixForm[Transpose[Q].Q]
MatrixForm[RT]
MatrixForm[R-RT]

It is worth noting the following: the diagonal values in R tend to come out in decreasing numerical order and the largest "zero" entries tend to be in high indices.  This is because the algorithm tends to pull out big bits first and errors build up through the computation and pollute the little bits that are left until the end.

Testing: Large

The same pattern holds true for larger examples.

{m,n}={400,400};
A= RandomReal[{-1,1},{m,n}];
{Q,R}=CGS[A];
MatrixPlot[
ΔI=Abs[Transpose[Q].Q-IdentityMatrix[m]],
PlotLabel->StringForm["Max[ΔI]=`` at ``",MaxΔ=Max[ΔI], Position[ΔI,MaxΔ]]
]
ListPlot[Diagonal[R]]

Timing

TimingData=Table[
m=Floor[(2^p)];
n=m;
A= RandomReal[{-1,1},{m,n}];
{m, Timing[CGS[A];][[1]]},
{p,4,9,0.5}];
ListLogLogPlot[TimingData]

Modified

A direct implementation of algorithm 8.1

Clear[MGS,A]
MGS[A_]:= Module[
{m,n,V,Q,R,i,j,q, v,r},
{m,n}=Dimensions[A];
V=A;
Q=Table[0,{m},{n}];
R=Table[0,{n},{n}];
Do[
v=V[[All,i]];
r = Norm[v];
q=v/r;
Q[[All,i]]=q;
R[[i,i]]=r;
Do[
v=V[[All,j]];
r=q.v;
v=v-r q;
V[[All,j]]=v;
R[[i,j]]=r,
{j,i+1,n}
],
{i,1,n}
];
{Q,R}
]

Testing the code.  The output Q should be approximately orthogonal and have the same column space as the input. It is easy to check small test problems for Orthogonality and that the code can easily genearte an R matrix

Testing: Small

{m,n}={4,4};
A= RandomReal[{-1,1},{m,n}];
{Q,R}=MGS[A];
MatrixForm[Transpose[Q].Q]
MatrixForm[R]
MatrixForm[Transpose[Q].A]

Testing: Large

The same pattern holds true for larger examples.

{m,n}={400,400};
A= RandomReal[{-1,1},{m,n}];
{Q,R}=MGS[A];
MatrixPlot[
ΔI=Abs[Transpose[Q].Q-IdentityMatrix[m]],
PlotLabel->StringForm["Max[ΔI]=`` at ``",MaxΔ=Max[ΔI], Position[ΔI,MaxΔ]]
]
ListPlot[Diagonal[R],
PlotLabel->"Diagonal Entries of R"]

Timing

TimingData=Table[
m=Floor[(2^p)];
n=m;
A= RandomReal[{-1,1},{m,n}];
{m, Timing[MGS[A];][[1]]},
{p,4,9,0.5}];
ListLogLogPlot[TimingData]

Fit[Log[TimingData],{1,t},t]

Recreation of Plot

To make the picture in the text I need to create a hard matrix with a wide range of SVs.

dim=80;
{U,X}=QRDecomposition[RandomReal[{-1,1},{dim,dim}]];
{V,X}=QRDecomposition[RandomReal[{-1,1},{dim,dim}]];
InVals=Table[(2^i),{i,-1.0,-dim,-1}];
S=DiagonalMatrix[InVals];
A=U.S.V;
ListLogPlot[{
InVals,
Fred=Diagonal[CGS[A][[2]]],
Barney=Diagonal[MGS[A][[2]]]
},
PlotRange->All,
PlotStyle->{Red,Green, Blue}]

Table[
ϵ=(2.0^(-i));
{i,ϵ,1.0==1.0-ϵ,1.0==1.0+ϵ},
{i,45,48}]