What is the correct way to perform the Gram-Schmidt process?

Question

I'm trying to do a simple processing in Mathematica. My original code was in MATLAB. Part of the computation is the Gram-Schmidt process, which is an iterative calculation. Every vector is changed according to its predecessors. I hold all vectors as rows of a matrix.

1. $r1' = r1 / |r1|$
2. $r2' = r2 - (r2.r1')r1'/|r2 - (r2.r1')r1'|$
3. $r3' = r3 - (r3.r2')r2' - (r3.r1')r1' / \text{norm-of-the-numerator}$
4. ...

Now, it is clear that, for example

$$r2' = r2 - (r2.r1 / |r1|)r1 / |r1|/|r2 - (r2.r1 / |r1|)r1 / |r1||$$

since I replaced $r1'$ according to (1). If I understand Mathematica correctly, it remembers all this stages. So when I get to (6) $r6' = \ldots$ the computations is taking too long. I aborted after a few minutes, I have to get up to 11.

How can I deal with it? I read about memoization and tried it. This is what I wrote

gsstep[mat_, v_, rr_] := 
  (lp  = v.mat[[rr]]; 
   lv = v - lp*mat[[rr]];
   lv)

  For[r = 1, r <= dw1, r++,
    Print["r = ", r]
    Clear[v];
    v = mW[[r]];
    For[rr = 1, rr <= r - 1, rr++, v = gsstep[mW, v, rr]];
    Print["done"];
    norm = Norm[v];
    mW[[r]] = v/norm;];

Unfortunately, this didn't do the trick.

Can someone please give me some direction about what the correct way to do this in Mathematica?

Answer 1

Well, as others have said, you could use Orthogonalize to solve this problem, and it will likely be the fastest way to do so, but let's assume that (for whatever reason) you don't want to use Orthogonalize. I'll assume it's something of a learning exercise, so we'll create a function (which I end up calling gsOrthognalize), and we want to pass to it a list of vectors (i.e. a matrix), and we want it to output the Gram-Schmidt orthogonalization. We can write the function as follows:

gsOrthogonalize[mat__] :=
 Quiet@Module[{input = mat, rr, rtemp},
   rr = {input[[1, All]]/Norm[input[[1, All]]]};
   For[i = 2, i <= Length[input], i++,
    Clear@rtemp;
    rtemp = input[[i, All]];
    For[j = 1, j < i, j++,
     rtemp = rtemp - (rr[[j, All]].input[[i, All]]) rr[[j, All]];
     ];
    rtemp = rtemp/Norm@rtemp;
    rr = rr~Join~{rtemp};
    ];
   rr // FullSimplify
   ]

Now, when we call gsOrthoganlize on a matrix, we will get an orthonormal set out.

As an example:

tmat = {{1, 0, 1}, {2, 6, 3}, {1, 1, 1}, {2, 3, 5}};
gsOrthogonalize[tmat]//MatrixForm
Orthogonalize[tmat]//MatrixForm

we get

$$ \left( \begin{array}{ccc} \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{146}} & 6 \sqrt{\frac{2}{73}} & \frac{1}{\sqrt{146}} \\ \frac{6}{\sqrt{73}} & \frac{1}{\sqrt{73}} & -\frac{6}{\sqrt{73}} \\ 0 & 0 & 0 \\ \end{array} \right)$$

as the output from both our function gsOrthogonalize as well as the Orthogonalize command.

Answer 2

As Michael Witt has mentioned this is mainly interesting as a learning exercise as there is Orthogonalize. As such, I think it is worth mentioning that it is quite easy to do this a lot more efficient by simplifying each row before calculating the next row which will keep the underlying expressions a lot simpler and smaller. Here is code which does do that (and is also a bit simplified and optimized):

gsOrthogonalize[mat_?MatrixQ] := Module[{
    rr = mat, rtemp
  },
  rr[[1]] = mat[[1, All]]/Norm[mat[[1, All]]];
  Do[
    rtemp = mat[[i, All]];
    Do[
      rtemp = rtemp - (rr[[j, All]].mat[[i, All]]) rr[[j, All]],
      {j, 1, i - 1}
    ];
    rtemp = rtemp/Quiet[Norm[rtemp], N::meprec];
    rr[[i]] = FullSimplify[rtemp];
    ,
    {i, 2, Length[mat]}
  ];
  rr
];

Some notes: I have removed some unnecessary local variables and the Clear which also seems not to be necessary. I also am "preallocating" rr by just setting it to mat initially and then overwrite the entries within the outer loop, which is somewhat more efficient than the Join within the loop. But what really makes a difference is that I put the FullSimplify into the outer loop, so that consecutive rows will be calculated with an already simplified version of the previous rows. Then I think it is good practice to avoid the use of Quiet but where not possible limit it to that part of the code which gives the unwanted messages. I also have changed the argument pattern to check that what we get is actually a matrix. And finally I have changed For loop with a Do loop, as I think a For loop always makes code more complicated than necessary.

The net result is that this version will solve the integer 5x5 case almost immediately, where the original code will take several minutes on my machine.

Answer 3

Here is a functional programming approach which is very fast.

First we define a function that computes the projection of v onto u.

proj[u_, v_] := v.u/u.u u

Next is the definition of the Gram-Schmidt orthonormal function.

gsOrthonormal[vectorSet_?MatrixQ] :=
 Map[
  With[
    {
     norm = Norm[#]
     },
    If[norm > 0,
     #/Norm[#],
     #]
    ] &,
  Fold[
   Function[{spanVec, v},
      Join[
       spanVec,
       {v - Total@Map[proj[#, v] &, spanVec]}
       ]
      ][#1, #2] &,
   {vectorSet[[1]]}, Drop[vectorSet, 1]
   ]]

The gsOrthonormal function collects a list called spanVec using Fold. A key feature of Fold is the ability to use the previous result in subsequent iterations.

Below spanVec is shown for the first few steps assuming the input matrix is called vs.

1. spanVec = {vs[[1]]}
2. spanVec = {vs[[1]], vs[[2]] - proj[vs[[1]],vs[[2]]}

Note that after the second step, spanVec[[2]] is orthoganal to spanVec[[1]].

3. spanVec = {vs[[1]], spanVec[[2]], vs[[3]]-proj[vs[[1]],vs[[3]] - proj[spanVec[[2]],vs[[3]]}

spanVec is the output of Fold. Then it is wrapped in Map and normalized.

Validation

The function is tested using

vectorSet = {{1, 1, 0}, {0, 1, 1}, {1, 0, 1}};

a simple three dimensional set of vectors.

gsOrthonormal[vectorSet]

produces

This matches the output of Orthoganlize.

SameQ[Orthogonalize[vectorSet], gsOrthonormal[vectorSet]]
(* True *)

Timing

We will use a vector set where the dimension is eleven.

vectorSet11 = RandomInteger[{-9, 9}, {11, 11}];

Orthogonalize[vectorSet11]; // Timing
(* {0.0936006, Null} *)

The result using gsOrthonormal is actually faster than the built-in

gsOrthonormal[vectorSet11]; // Timing
(* {0.0780005, Null} *)

Here is the timing for Albert Retey's procedural version.

gsOrthogonalize[mat_?MatrixQ] := Module[{
    rr = mat, rtemp
  },
  rr[[1]] = mat[[1, All]]/Norm[mat[[1, All]]];
  Do[
    rtemp = mat[[i, All]];
    Do[
      rtemp = rtemp - (rr[[j, All]].mat[[i, All]]) rr[[j, All]],
      {j, 1, i - 1}
    ];
    rtemp = rtemp/Quiet[Norm[rtemp], N::meprec];
    rr[[i]] = FullSimplify[rtemp];
    ,
    {i, 2, Length[mat]}
  ];
  rr
];

takes 18 seconds on my machine

gsOrthogonalize[vectorSet11]; // Timing

{18.2521, Null}

Michael Witt's procedural version doesn't finish on my machine after waiting over 30 minutes.

Answer 4

Just for reference, here is a very compact implementation of the modified Gram-Schmidt algorithm (MGS, see this and this as well):

mgs[mat_?MatrixQ] := Module[{a = mat, n},
    n = Length[mat];
    Do[a[[k]] = Normalize[a[[k]]];
       Do[a[[j]] -= (a[[j]].a[[k]]) a[[k]],
          {j, k + 1, n}], {k, n}];
    a]

As noted by Golub and Van Loan, the modified version is much more stable in inexact arithmetic than the classical version of the algorithm.

Some tests:

mgs[{{1, 0, 1}, {2, 6, 3}, {1, 1, 1}, {2, 3, 5}}] // FullSimplify // MatrixForm

$$\begin{pmatrix} \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{146}} & 6 \sqrt{\frac{2}{73}} & \frac{1}{\sqrt{146}} \\ \frac{6}{\sqrt{73}} & \frac{1}{\sqrt{73}} & -\frac{6}{\sqrt{73}} \\ 0 & 0 & 0 \\ \end{pmatrix}$$

With[{ε = 2^-52}, 
     mgs[{{1, ε, 0, 0}, {1, 0, ε, 0}, {1, 0, 0, ε}}] // Simplify] // N
   {{1., 2.22045*10^-16, 0., 0.}, {1.57009*10^-16, -0.707107, 0.707107, 0.},
    {9.06493*10^-17, -0.408248, -0.408248, 0.816497}}

With[{ε = $MachineEpsilon}, 
     mgs[{{1., ε, 0., 0.}, {1., 0., ε, 0.}, {1., 0., 0., ε}}]]
   {{1., 2.22045*10^-16, 0., 0.}, {0., -0.707107, 0.707107, 0.},
    {0., -0.408248, -0.408248, 0.816497}}

With[{ε = $MachineEpsilon}, 
     Orthogonalize[{{1., ε, 0., 0.}, {1., 0., ε, 0.}, {1., 0., 0., ε}}]]
   {{1., 2.22045*10^-16, 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}}

where we see that mgs[] is in fact able to return better results than Orthogonalize[] for this case.

Answer 5

The answers already given obscure the simple geometric idea behind the Gram-Schmidt process. To see the geometric idea, one could start with:

(* project a vector onto the span of some orthonormal vectors *)
project[u_List, x_?VectorQ] := Total[(x.#) # & /@ u]

(* do one step of the Gram-Schmidt:
   given a list u of vectors that are already orthonormal...
   result is what will be the next unit vector. *) 
gsStep[u_List, vi_?VectorQ] := Normalize[vi - project[u, vi]]

Now one may use gsStep iteratively, recursively, or functionally (with Fold) to accomplish the complete Gram-Schmidt process.

No claims for efficiency are made here!