Math 236  •  12 April 2001

Iterative and recursive programming 
of the Gram - Schmidt algorithm

This notebook explores both iterative (explicitly looping) and recursive ways to express the Gram-Schmidt algorithm in Mathematica.  For further discussion and examples of iteration and recursion, see the separate notebook Factorial.nb.

Auxiliary functions: project, norm, normalize

Orthogonal projection project onto the span of orthnormal vectors

The function project

The first task is to write a function to calculate the orthogonal projection of a vector onto the span of a given list of orthonormal vectors.  (So this is entirely different from the method you are using for your function proj, which projects onto the span of an arbitrary list of linearly independent vectors.)  The definition following is a direct expression of our mathematical formula for this:

project[u_List,  x_?VectorQ]  :=  
      Sum[(x  .  
              u[[
                j]]) u[[
            j]],  {j,  1,  Length[u]}]

If you wish, instead of Sum, you may use summation notation with the aid of the BasicInput palette:

projectWithSigma[u_List,  
        x_?VectorQ]  :=  (∑_(j = 1)^(Length[u]))(x  .  
              u[[
                j]]) u[[
            j]]

A set of orthonormal vectors

Here is a test set of three vectors in (ℝ^4):

onVectors  =  {{(2/(Sqrt[21])), -(4/(Sqrt[21])), 0, (1/(Sqrt[21]))}, {(1/(Sqrt[15])), 
            (1/(Sqrt[15])), -(Sqrt[(3/5)]), (2/(Sqrt[15]))}, {(1/(Sqrt[10])), 
            (1/(Sqrt[10])), (Sqrt[(2/5)]), (Sqrt[(2/5)])}};

These vectors really are orthonormal, as you an check by forming the dot products of all their pairs.

n  =  Length[onVectors];  
    Table[onVectors[[i]]  .  
          onVectors[[j]],  {i,  1,  
          n},  {j,  1,  n}]  //  TableForm

A more concise way of obtaining the same result is functionally, by using the function Outer:

Outer[Dot,  onVectors,  onVectors,  1]   //  TableForm

How do you think the list onVectors was concocted?  (The secret: start with a list of independent vectors and apply the Gram-Schmidt process to it!)

Testing the function project

project[
      onVectors,   {2,  0,  -1,  1}]

Ugh!  An attempt at simplification is surely called for:

project[onVectors,  {2,  0,  -1,  1}]  //  
      Simplify

The summation counter is local!

Notice that the variable j, used as a counter in the Sum definition of project, does not now have a value:

j

Now we assign a value to the name j to make a point about Mathematica.

j  =  "Hello, world!";

Let's calculate with project again.

project[onVectors,  {2,  0,  -1,  1}]  //  
      Simplify

j

Even though the counter j used in the definition of project (and projectWithSigma) seems to take successive values from 1 through Length[u], nonetheless after Sum has done its work the variable j still has its original value!  Why?  Because the counter j in a Sum expression automatically is made a local variable: evaluating the Sum expression actually creates, behind the scenes, a new variable and uses that in place of j.  Hence when Sum is done it has not assigned any original or new value to your variable j.

Expressing project functionally

If you wish, you can form the needed sum in a functional way through use of Apply, Map, and the pure function x . #&.  This pure function is just the dot product in which the first of the two variables is fixed to be a given vector x; you can read as "dot product with x".

projectFunctional[u_List,  x_?VectorQ]  :=  
      Apply[Plus,  Map[(x  .  #) # &,  u]]

projectFunctional[onVectors,  {2,  0,  -1,  1}]  //  
      Simplify

Normalizing a vector: norm and normalize

norm[v_?VectorQ]  :=  (Sqrt[v  .  v])

normalize[
        v_?VectorQ]   :=  (1/norm)[v] 
          v  /;  ! (v   ==  0 v)

normalize[{3,  4,  0}]

Observe the condition expression
	/; !(v  == 0 v)
in the definition of normalize.  The expression !(v  == 0 v) says that v is not a zero vector.  Then the entire expression /; !(v  == 0 v) may be read as "provided that v is not a zero vector".  The expression normalize[vec] for an argument value vec is actually evaluated only in the case that vec is not a zero vector; otherwise, when vec is a zero vector, the entire expression normalize[vec] is returned without evalation of normalize.  This is the way we guard against accidentally trying to divide by the length, 0, of a zero vector.

normalize[{0, 0, 0}]

normalize[{1, 2, 3}  -  {1, 2, 3}]

Test vectors to orthonormalize

Here is a test set of three vectors in (ℝ^4):

vecs  =  { {2,  -4,  0,  1},  {0,  2,  -2,  
            1},  {0,  0,  2,  0} };

Exercise:  Use Mathematica to check that these vectors really are linearly independent.

General step in Gram-Schmidt algorithm: function gsStep

The function gsStep carries out one step of the Gram-Schmidt process.  Given a list u of vectors that are already orthonormal and a "new" vector called here vi, that is not already in the span of the vectors in u, then the result of gsStep[u, vi] is what will become the next unit vector, orthogonal to all the vectors already in u.  (The actual work of appending this newly obtained  vector gsStep[u, vi] to the list u is what will be in the main function gs (see the next section).

gsStep[u_List,  vi_?VectorQ]  :=  
      normalize[vi  -  project[u,  vi]]

Iterative implementation of Gram-Schmidt algorithm using a Do loop

gs[v_List]  :=  
      Module[{i,  
          u}, 
u  =  {normalize[
              v[[1]]]}; 
Do[
          u  =  Append[u,  
              gsStep[u,  
                v[[i]]]],  {i,  2,  
            Length[v]}]; 
u 
]

Now test gs with these vectors:

gs[vecs]

The result is, of course, precisely the test set of orthonormal vectors used earlier to test the function project.  That they are orthonormal was verified earlier.

Exercise:  Check that these orthonormal vectors have the same span as the vectors vecs.

Note that the following would not work:

	gs[v_List] := Module[{i, u},
	    u = {normalize[v[[1]]]};
 	    Do[u = Append[u, gsStep[u, v[[i]]]], {i, 2, Length[v]}]
 	 ]
 	 
The reason it would not work is explained in the subsection "A few words abot Do loops".

A few words about Do loops

Note that Do is not any special kind of "keyword" in Mathematica; it is just another function!  Unlike most functions, however, it returns a Null result, which therefore normally displays nothing in an output cell.  For example:

Do[(i^2),  {i,  1,  10}]

No output!

Do[(i^2),  {i,  1,  10}]  ==  Null

If you had really wanted to write a Do loop that would form a list of squares, you could do it in several ways.  Here are three:

Do[Print[(i^2)],  {i,  1,  5}]

tbl  =  {};
Do[
        tbl  =  Append[tbl,  (i^2)],  {i,  1,  
          5}];
tbl

Clear[sqs];
Do[sqs[i]  =  (i^2),  {i,  1,  5}];
Table[sqs[i],  {i,  1,  5}]

Of course this is just a "toy" example:  If you really wanted to produce a list of the first 10 perfect squares, you could simply evaluate the following:

Range[5]^2

The next thing to know is that Do, like Sum, treats the counter—iteration variable i (or whatever name you used for it) as a local variable: it doesn't leave i with any value when it's done, and it doesn't change any value that may already have been assigned to i.  For example:

Clear[i];
Do[(i^2),  {i,  1,  5}];
i

i  =  "Today is the beginning of the rest of your life.";
Do[(i^2),  {i,  1,  5}];
i

Iteration through functional programming

The simplest way to implement the Gram-Schmidt iteratively is undoubtedly the one shown earlier, using a Do loop.  Mathematica also provides the functions Nest (along with the related function NestList) and Fold (along with the related function FoldList) that allow you to accomplish the same iteration in a terser fashion—provided you can describe what you are doing in terms of starting from some initial value (Subscript[x, 0]) and repeatedly evaluating the same function f some number n of times so as to obtain the last of the iterates:
	(Subscript[x, 0]),  (Subscript[x, 1])  =  f((Subscript[x, 0]) ),  
      (Subscript[x, 2])  =  f((Subscript[x, 1]))  =  f(f((Subscript[x, 0]))),  
      (Subscript[x, 3]) =  f((Subscript[x, 2]))  =  
          f(f(f((Subscript[x, 0])))),  … (Subscript[x, n])  =  f(
            (Subscript[x, n  -  1]))  =  
          f(f(f( … f((Subscript[x, 0])) …)))
The catch is to be able to cook up the necessary function f.

In our situation of the Gram-Schmidt process, it's pretty clear what the initial value should be: the normalized form of the first of the original list of vectors, that is:
	normalize[First[v]]
And it seems that the function that needs to be iterated should be gsStep.  Unfortunately, gsStep takes two arguments, not one.  That's why Fold is a more promising construct to use.  We leave it as an exercise for you to learn about Fold and to use it to implement the Gram-Schmidt algorithm functionally.

Other forms of iteration

Besides Do, the functions While and For also allow you to perform iteration explicitly step-by-step.  In a situation such as the one here, where you know ahead oftime that here are a fixed number of steps to the iteration, it is almost always preferable to use Do instead of While or For; the reason is that you specifying the initial and final value of the counter is more complicated with While and For, and you have to explicitly increment the counter with an argument such as i = i + 1 (or the abbreviated form i++).  Of course, if you need delicate control over the iteration—because, for example, you do not know ahead of time the number of iteration steps, or because you need to change the values of the counter in some irregular way, then by all means consider a While or For.

Recursive version of Gram-Schmidt algorithm

First write a rule for the base step—for a list consisting of a single (nonzero) vector.  Note that the result needs to be again a list consisting of a single vector.

gsRecur[{v_?VectorQ}]  :=  {normalize[
            v]}  /;  !  (v  ==  0 v)

Try gsRecur with the first of the test set vecs of vectors:

vecs[[1]]

{2,  0,  -1,  1}

gsRecur[ {vecs[[1]] }]

Exercise:  Check that the preceding result is a unit vector that has the same span as vecs[[1]].

Now write the recurrence relation, which tells how to effect the Gram-Schmidt process on an arbitrary list of some number k (linearly independent) vectors:  Carry out the Gram-Schmidt process on the first k  -  1 of the vectors—by means of the instructions described in the function gsRecur—to obtain a list u of  k  -  1 orthonormal vectors.  Then append to those the kth orthonormal vector as obtained from the kth of the original vectors (and the k  -  1 orthonormal vectors) by using gsStep.

gsRecur[v_List]  :=  
      Module[{u},  
u  =  
          gsRecur[Drop[v,  -1]]; 
Append[u,  
          gsStep[u,  Last @ v]]
]

Try this recursive function on the same test set of linearly independent vectors as before:

gsRecur[vecs]

Notice that you don't have to tell Mathematica in any further detail that in order to effect the Gram-Schmidt process on the first k  -  1 of the original vectors, it now has to in turn first effect the process on the first k  -  2 of those, etc.  All that is taken care of automatically by the recurrence relation.  This recurrence relation becomes one of the two rules that Mathematica knows about gsRecur:

?gsRecur

As is the case with any rules entered into Mathematica's database of rules at a session, when you ask Mathematica to evaluate an expression involving gsRecur, it first tries to use the special rule (for a list of one nonzero vector).  It that is not relevant, then Mathematica tries to use the general rule.  But when it tries to use the general rule, Mathematica finds it has either to use the special rule (in case the original list consisted of only two vectors) or else it has to use the general rule again.  And when evaluating any expression, Mathematica tries to use all rules it knows, as many times as it needs to use them, until the value of the expression undergoes no further changes.  That is why the "recusrive descent" downwards works automatically.