Orthogonalization

by Rolf Sulanke

Institut für Mathematik der Humboldt Universität zu Berlin
http://www2.mathematik.hu-berlin.de/~sulanke
sulanke@mathematik.hu-berlin.de

Sept  7, 2003

Summary

In this notebook two orthogonalization procedures: esorthonorm, contained in our package euvec.m, and GramSchmidt, contained in the standard Mathematca package LinearAlgebra/Orthogonalization.m, are comparad. It results that esorthonorm is  10 times faster than GramSchmidt in orthonormalizing the first 25 polynomial to construct the Legendre polynomials. Besides of this, GramSchmidt produces its output in an unsimplified manner, what gives a very long and complicated result, and, may be, is the reason of the longer time needed to proceed. For the test the package  euvec.m must not be loaded; the definition of esorthonorm is reproduced in this notebook.  The test has been carried out on a Dell Inspiron 3200 laptop.

1. The Module esorthonorm

\(\(\(esorthonorm::usage = "esorthonorm[b_List, opts] is the Erhard Schmidt
orthogonalization for the list b of vectors with respect to a positive
semi-definite scalar product which can be specified as an option; the
default is Dot."\)\(
\)
    \)\)
\(\(\(innerprod::usage = "innerprod is an option for esorthonorm and orthonorm;the
default is Dot (esorthonorm), resp. pssp (orthonorm)."\)\(
\)
    \)\)
\(\(\(normed::usage = "normed is an option for esorthonorm and orthonorm which
 specifies whether
the vectors produced by esorthonorm or orthonorm are normed, or only
 orthogonal;the default
is True."\)\(
\)
    \)\)
\(\(\(print::usage = "print is an option for esorthonorm and orthonorm which specifies
 whether
the (unnormed) vectors produced by esorthonorm or orthonorm are printed during
 the
orthogonalization procedure;the default is False."\)\(
\)
    \)\)
neglect::usage = "neglect is an option for esorthonorm and \
orthonorm which
 specifies whether
vectors with components or scalar square less than 10^neglect
produced by esorthonorm  or orthonorm during the orthogonalization procedure \
are
 set to
zero;the default is neglect -> -10.For neglect -> 0 no neglecting takes \
place."

\(\(\(Options[esorthonorm] = {innerprod -> Dot, normed -> True, 
        print -> False, neglect -> \(-10\)}\)\(
\)
    \)\)
esorthonorm[m_List, opts___] := 
      Module[{ma = m, mb = {}, k = Length[m], rd = Length[m], j = 0, i, ne, 
          prod, pr, bn, sm, aus, h, 
          f}, {prod, nor, pr, 
            bn} = {innerprod, normed, print, neglect} /. {opts} /. 
            Options[esorthonorm]; 
If[
          Head[m[[1]]] == List, 
          rd = Length[m[[1]]]]; 
f[h_, bn_: - 10] := 
          If[bn == 0, h, Chop[h, 10^bn]]; 
sm[aus_] := 
          MapAll[f, aus]; 
Label[ne]; 
If[
          Length[mb] ≥ Min[rd, k], 
          Return[mb]]; 
j = j + 1; 
If[
          j > k, Return[mb]]; 
ma[[j]] = 
          sm[Simplify[
              ma[[j]] - 
                Sum[mb[[i]]*
                    prod[ma[[j]], mb[[i]]]/
                      prod[mb[[i]], mb[[i]]], {i, 1, 
                    Length[mb]}]]]; 
If[pr == True, 
          Print[ma[[j]]]]; 
If[
          prod[ma[[j]], ma[[j]]] == 0, Goto[ne], 
          If[nor == True, 
            AppendTo[mb, 
              sm[Simplify[
                  ma[[j]]/Sqrt[prod[ma[[j]], ma[[j]]]]]]], 
            AppendTo[mb, Simplify[ma[[j]]]]], 
          If[nor == True, 
            AppendTo[mb, 
              sm[Simplify[
                  ma[[j]]/Sqrt[prod[ma[[j]], ma[[j]]]]]]], 
            AppendTo[mb, Simplify[ma[[j]]]]]]; 
Goto[
          ne]]

2. The Legendre polynomials with esorthonorm and GramSchmidt

The procedure esorthonorm  may be applied also for function spaces, for which a positive semidefinite scalar product must be defined in terms of Mathematica.

prodint[f_, g_]  :=  Integrate[f*g, {x,  -1, 1}]

prodint[1, x^2]

seq[h_]  :=  Table[x^i, {i, 0, h}]

seq[3]

esorthonorm[seq[3],  innerprod  ->  prodint]

esorthonorm needed 26.27 seconds for the next orthonormalization:

Timing[esorthonorm[seq[25],  innerprod  ->  prodint]]

legesortho = %;

Now we load the package containing the GramSchmidt orthogonalization:

<< "LinearAlgebra/Orthogonalization.m"

?GramSchmidt

Options[GramSchmidt]

GramSchmidt needed 268.26 s for the same orthonormalization, and the result looks very unsimplified::

Timing[GramSchmidt[seq[25],  InnerProduct  ->  prodint]]

gsortho  =  %;

The results of both procedures coincide:

Simplify[gsortho[[2]]  -  legesortho[[2]]]

Length[%]