Gram-Schmidt Process

GramSchmidt[V_]:= Module[
{n,m,Q=0*V,qi,qj},
{n,m}=Dimensions[V];
Do[
qi=V[[All,i]];
Do[
qj=Q[[All,j]];
qi=qi-(qi.qj)qj,
{j, 1,i-1}];
Q[[All,i]]=qi/Norm[qi],
{i,1,m}];
Q
]

A small test.

{n,m}={4,3};
V=RandomReal[{-1,1},{n,m}];
Q=GramSchmidt[V];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
MatrixPlot[OrthCheck]
MatrixForm[OrthCheck]

Testing

A larger test.

{n,m}={700,689};
V=RandomReal[{-1,1},{n,m}];
Timing[Q=GramSchmidt[V];]
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
MatrixPlot[Log[10,Abs[OrthCheck]]]

There seems to be a pattern!  With larger errors in the far corner.

Histogram[ 
Log[10,Abs[OrthCheck[[{1,123,680}]]]],
ChartStyle->{Red,Green,Blue}
]

Histogram[ 
Transpose[Log[10,Abs[OrthCheck[[All,{1,123,680}]]]]],
ChartStyle->{Red,Green,Blue}
]

Histogram[ 
Log[10,Abs[Diagonal[OrthCheck]]],
ChartStyle->{Pink}
]

A torture test.

Instead of just taking larger matrices we can take harder matrices.  The vandermonde matrix is a classical ill-behaved matrix.

Vandermonde[n_,m_]:= Module[
{p},
Transpose[
Table[
p=N[j/m];
Table[(p^(i-1)),{i,1,n}],
{j,m}]
]
]

Lets see how the poor thing works

{n,m}={8,7};
V=Vandermonde[n,m];
Q=GramSchmidt[V];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]

Lets see if we can fix it.  Lets try doing it twice.

{n,m}={18,17};
V=Vandermonde[n,m];
Q=GramSchmidt[V];
Q=GramSchmidt[Q];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]

Is three times better?

{n,m}={80,79};
V=Vandermonde[n,m];
Q=GramSchmidt[V];
Q=GramSchmidt[Q];
Q=GramSchmidt[Q];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]

A slightly fancier implementation

This one just deos the orthogonalization twice inside the loop.

GramSchmidtDouble[V_]:= Module[
{n,m,Q=0*V,qi,qj},
{n,m}=Dimensions[V];
Do[
qi=V[[All,i]];
Do[
Do[
qj=Q[[All,j]];
qi=qi-(qi.qj)qj,
{j, 1,i-1}],
{2}];
Q[[All,i]]=qi/Norm[qi],
{i,1,m}];
Q
]

It does better!

{n,m}={80,79};
V=Vandermonde[n,m];
Q=GramSchmidtDouble[V];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]

You can reorder the operations but this one is not so good!

GramSchmidtDoubleReorder[V_]:= Module[
{n,m,Q=0*V,qi,qj},
{n,m}=Dimensions[V];
Do[
qi=V[[All,i]];
Do[
qj=Q[[All,j]];
Do[
qi=qi-(qi.qj)qj,
{2}],
{j, 1,i-1}];
Q[[All,i]]=qi/Norm[qi],
{i,1,m}];
Q
]

{n,m}={80,79};
V=Vandermonde[n,m];
Q=GramSchmidtDoubleReorder[V];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]

In fact it is just about the same as the original!

{n,m}={80,79};
V=Vandermonde[n,m];
Q=GramSchmidt[V];
OrthCheck=IdentityMatrix[m]-Transpose[Q].Q;
Norm[OrthCheck]
ListContourPlot[Log[10,Abs[OrthCheck]],
ContourLabels->All]