# Gram-Schmidt Process for Polynomials

I want to implement the Gram Schimdt procedure to the vector space of polynomials of degree up to 5, i.e. I want to find an orthogonal basis from the set of vectors $v=(1,x,x^2,x^3,x^4,x^5)$. The inner product is defined as $$\langle f,g\rangle=\int_{-\pi}^{\pi}f(x)g(x)dx.$$ The Gram Schimdt process is to first find $u$ such that \begin{align} &u_1=v_1\\ &u_i=v_i -\sum_{j=1}^{i-1}\frac{\langle u_j,v_i\rangle}{\langle u_j,u_j\rangle}u_j,i\ge 1 \end{align} and then the orthogonal basis is found by letting $e_i=u_i/\|u_i\|$. I wrote the following code to find $u$, but it went to endless loop. I am an amatur in mathematica and don't know where I am wrong.

v = Function[x, Evaluate[x^#]] & /@ Range[0, 5];
u = v;
For[i = 2, i < 7, i++,
u[[i]] = Function[x,
v[[i]][x] -
Sum[Integrate[v[[i]][x]*u[[j]][x], {x, -Pi, Pi}]/
Integrate[u[[j]][x]^2, {x, -Pi, Pi}]*u[[j]][x], {j, 1, i - 1}]
]
]


I am reading Linear Algebra Done Right P108 and want to implement this by hand to deepen my understanding of this topic and I also want to learn more about mathematica functions. Thanks for any help.

• You could do this with Orthogonalize[x^Range, Integrate[#1*#2, {x, -\[Pi], \[Pi]}] &] – acl Mar 25 '13 at 15:05
• Thanks this works perfect! But I am reading a linear algebra textbook and want to implement this by hand to deepen my understanding of this topic. – user6193 Mar 25 '13 at 15:07
• I happen to be teaching exactly this topic in my PDE class today. My Legendre polynomial notebook is the last link on our class webpage: sites.google.com/a/unca.edu/mark-mcclure/home/classes/… – Mark McClure Mar 25 '13 at 16:43

You can do this using built-in functions using

Orthogonalize[x^Range, Integrate[#1*#2, {x, -\[Pi], \[Pi]}] &]


By hand, a literal implementation of your equations is

inn = Integrate[#1*#2, {x, -\[Pi], \[Pi]}] &;
v = x^Range;
u = ConstantArray[0, Length@v];
u[] = v[];
Do[
u[[i]] = v[[i]] - Total@Table[
u[[j]]*inn[v[[i]], u[[j]]]/inn[u[[j]], u[[j]]],
{j, 1, i - 1}],
{i, 2, Length@v}
]

u

(*{x, x^2, -((3 \[Pi]^2 x)/5) + x^3, -(5/7) \[Pi]^2 x^2 +
x^4, -((3 \[Pi]^4 x)/7) + x^5 -
10/9 \[Pi]^2 (-((3 \[Pi]^2 x)/5) + x^3)}*)


You can see that they are orthogonal but not orthonormal:

Outer[inn, u, u] // MatrixForm • The basis consists of six vectors so I guess in my code Range should be replaced by Range[0,5]. – user6193 Mar 25 '13 at 17:39
• Yes, if you want {1,x,x^2} for instance, that is what you need. I did not pay much attention to that sorry. – acl Mar 25 '13 at 17:40
• Never mind. Thanks for your brilliant code. I googled but found nothing about the @ you used in your code(It worked perfectly). Is it an alternative to [] ? It makes the code much handy. – user6193 Mar 25 '13 at 17:46
• if you put the cursor after @ and press F1, you'll be taken to the doc page. In brief, yes, f@x is the same as f[x] – acl Mar 25 '13 at 17:52

This will not answer your question directly, but seems generally on-topic. I have replied to a similar question before on StackOverflow, producing a complete and general Gram-Schmidt process implelentation which I believe to be rather idiomatic for Mathematica. I then specialized it to polynomials. Since it is not going to be migrated, I will simply reproduce my answer here verbatim.

I will throw in a complete code for Gram - Schmidt and an example for function addition etc, since I happened to have that code written about 4 years ago. Did not test extensively though. I did not change a single line of it now, so a disclaimer (I was a lot worse at mma at the time). That said, here is a Gram - Schmidt procedure implementation, which is a slightly generalized version of the code I discussed here:

oneStepOrtogonalizeGen[vec_, {}, _, _, _] := vec;

oneStepOrtogonalizeGen[vec_, vecmat_List, dotF_, plusF_, timesF_] :=
Fold[plusF[#1, timesF[-dotF[vec, #2]/dotF[#2, #2], #2]] &, vec,  vecmat];

GSOrthogonalizeGen[startvecs_List, dotF_, plusF_, timesF_] :=
Fold[Append[#1,oneStepOrtogonalizeGen[#2, #1, dotF, plusF, timesF]] &, {},  startvecs];

normalizeGen[vec_, dotF_, timesF_] := timesF[1/Sqrt[dotF[vec, vec]], vec];

GSOrthoNormalizeGen[startvecs_List, dotF_, plusF_, timesF_] :=
Map[normalizeGen[#, dotF, timesF] &, GSOrthogonalizeGen[startvecs, dotF, plusF, timesF]];


The functions above are parametrized by 3 functions, realizing addition, multiplication by a number, and the dot product in a given vector space. The example to illustrate will be to find Hermite polynomials by orthonormalizing monomials. These are possible implementations for the 3 functions we need:

hermiteDot[f_Function, g_Function] :=
Module[{x}, Integrate[f[x]*g[x]*Exp[-x^2], {x, -Infinity, Infinity}]];

SetAttributes[functionPlus, {Flat, Orderless, OneIdentity}];
functionPlus[f__Function] :=   With[{expr = Plus @@ Through[{f}[#]]}, expr &];

SetAttributes[functionTimes, {Flat, Orderless, OneIdentity}];
functionTimes[a___, f_Function] /; FreeQ[{a}, # | Function] :=
With[{expr = Times[a, f[#]]}, expr &];


These functions may be a bit naive, but they will illustrate the idea (and yes, I also used Through). Here are some examples to illustrate their use:

In:= hermiteDot[#^2 &, #^4 &]
Out= (15 Sqrt[\[Pi]])/8

In:= functionPlus[# &, #^2 &, Sin[#] &]
Out= Sin[#1] + #1 + #1^2 &

In:= functionTimes[z, #^2 &, x, 5]
Out= 5 x z #1^2 &


Now, the main test:

In:=
results =
GSOrthoNormalizeGen[{1 &, # &, #^2 &, #^3 &, #^4 &}, hermiteDot,
functionPlus, functionTimes]

Out= {1/\[Pi]^(1/4) &, (Sqrt #1)/\[Pi]^(1/4) &, (
Sqrt (-(1/2) + #1^2))/\[Pi]^(1/4) &, (2 (-((3 #1)/2) + #1^3))/(
Sqrt \[Pi]^(1/4)) &, (Sqrt[2/3] (-(3/4) + #1^4 -
3 (-(1/2) + #1^2)))/\[Pi]^(1/4) &}


These are indeed the properly normalized Hermite polynomials, as is easy to verify. The normalization of built-in HermiteH is different. Our results are normalized as one would normalize the wave functions of a harmonic oscillator, say. It is trivial to obtain a list of polynomials as expressions depending on a variable, say x:

In:= Through[results[x]]
Out= {1/\[Pi]^(1/4),(Sqrt x)/\[Pi]^(1/4),(Sqrt (-(1/2)+x^2))/\[Pi]^(1/4),
(2 (-((3 x)/2)+x^3))/(Sqrt \[Pi]^(1/4)),(Sqrt[2/3] (-(3/4)+x^4-3 (-(1/2)+x^2)))/\[Pi]^(1/4)}