DLA Gram-Schmidt orthogonalization

Section B.6 Gram-Schmidt orthogonalization

Warning B.6.1.

This page contains several redefinitions of a Python function inprod. You will get incorrect results or errors if you Evaluate a Sage cell defining that function in one subsection below, and then Evaluate Sage cells that use a function by that same name in a different subsection below without evaluating the appropriate Sage cell near the beginning of that different subsection that redefines that inprod function to the appropriate formula used in the example in that particular subsection.

The best way to avoid this problem is to Evaluate every Sage cell in a particular subsection, from the beginning, in order.

🔗

Subsection B.6.1 An example in $\R^4$

🔗

Let's begin simply with the standard inner product on

$\R^4\text{.}$ Suppose we want to produce an orthogonal basis for the “hyperplane” in

$\R^4$ defined by the equation

$\begin{equation*} w + 2 x + 3 y - 4 z = 0 \text{.} \end{equation*}$

🔗

First, we can obtain a basis for this space by treating this single equation as a homogeneous linear system. In other words, we want a basis for the null space of the matrix

$\begin{equation*} \begin{bmatrix} 1 \amp 2 \amp 3 \amp -4 \end{bmatrix} \text{.} \end{equation*}$

🔗

Assigning parameters

$\begin{align*} x_2 \amp = r, \amp x_3 \amp = s, \amp x_4 \amp = t \end{align*}$

🔗

leads to basis

$\begin{equation*} \basisfont{B} = \{ (-2,1,0,0), (-3,0,1,0), (4,0,0,1) \} \text{.} \end{equation*}$

🔗

Let's load this basis into Sage.

xxxxxxxxxx
 
v1 = vector((-2,1,0,0))
v2 = vector((-3,0,1,0))
v3 = vector((4,0,0,1))
print(v1,',',v2,',',v3)

🔗

Is this an orthogonal basis? We can realize the dot product just by multiplying vectors — Sage will interpret the first vector as a row vector and the second as a column matrix.

xxxxxxxxxx
 
v1 * v2

🔗

Nope, not orthogonal.

🔗

So let's carry out Gram-Schmidt. Begin by setting

$\uvec{e}_1$ to be our first basis vector.

xxxxxxxxxx
 
e1 = v1
print('e1 =', e1)

🔗

Next we can compute

$\uvec{e}_2\text{.}$

xxxxxxxxxx
 
e2 = v2 - ( (v2*e1) / (e1*e1) ) * e1
print('e2 =',e2)

🔗

And finally

$\uvec{e}_3\text{.}$

xxxxxxxxxx
 
e3 = v3 - ( (v3*e1) / (e1*e1) ) * e1 - ( (v3*e2) / (e2*e2) ) * e2
print('e3 =',e3)

🔗

The fractions are kind of ugly, and scaling doesn't affect orthogonality, so let's clear them from

$\uvec{e}_2$ and

$\uvec{e}_3\text{.}$

xxxxxxxxxx
 
e2 = 5 * e2
e3 = 7 * e3
print(e1,',',e2,',',e3)

🔗

Let's make sure that these are orthogonal vectors.

xxxxxxxxxx
 
print(e1 * e2)
print(e1 * e3)
print(e2 * e3)

🔗

Yep, all orthogonal.

🔗

If we want an orthonormal basis, we can go one step further and normalize each of these vectors. This will re-introduce some fractions, but that can't be helped.

xxxxxxxxxx
 
e1 = e1 / sqrt(e1*e1)
e2 = e2 / sqrt(e2*e2)
e3 = e3 / sqrt(e3*e3)
print(e1)
print(e2)
print(e3)

🔗

Not pretty, but they're orthogonal unit vectors.

🔗

We know that these vectors are independent because they're orthogonal, and there are three of them, so as long as they all lie in our

$3$ -dimensional hyperplane, they will for a basis for it. We can test that by making sure they satisfy the hyperplane equation,

$\begin{equation*} w + 2 x + 3 y - 4 z = 0 \text{.} \end{equation*}$

xxxxxxxxxx
 
print(e1[0] + 2*e1[1] + 3*e1[2] - 4*e1[3])
print(e2[0] + 2*e2[1] + 3*e2[2] - 4*e2[3])
print(e3[0] + 2*e3[1] + 3*e3[2] - 4*e3[3])

🔗

(Remember that Sage, like Python, uses

$0$ -based indexing, so e1[0] means the first coordinate of vector e1.)

🔗

Subsection B.6.2 An example in $\C^4$

🔗

Here we will carry out the Gram-Schmidt portion of Example 37.4.7. It will get annoying typing conjugate for every complex dot product calculation, so let's start by creating a Python function to return inner product values. Remember that the order of multiplication matters for the complex dot product.

xxxxxxxxxx
 
def inprod(u,v):
    return v.conjugate() * u
print("inprod function loaded")

🔗

Now let's load our initial basis into Sage. Remember that Sage uses I for the imaginary root of

$x^2 + 1\text{.}$

xxxxxxxxxx
 
v1 = vector((1,1,I,I))
v2 = vector((1,-1,I,I))
v3 = vector((0,0,1,0))
v4 = vector((0,0,0,1))
print(v1,',',v2,',',v3,',',v4)

🔗

Let's double-check that these vectors are linearly independent, and hence a basis for

$\C^4\text{.}$

xxxxxxxxxx
 
column_matrix([v1,v2,v3,v4]).rank()

🔗

Since the rank is equal to

$\dim \C^4\text{,}$ these vectors do indeend form a basis, just not an orthogonal one.

🔗

To get an orthogonal basis, it's time to Gram-Schmidt. The last input line gets a bit long, so we have used the Python line-continuation character \ to break it up over multiple input lines.

xxxxxxxxxx
 
e1 = v1
e2 = v2 - (inprod(v2,e1) / inprod(e1,e1)) * e1
e3 = v3 - (inprod(v3,e1) / inprod(e1,e1)) * e1 - (inprod(v3,e2) / inprod(e2,e2)) * e2
e4 = v4 \
- (inprod(v4,e1) / inprod(e1,e1)) * e1 \
- (inprod(v4,e2) / inprod(e2,e2)) * e2 \
- (inprod(v4,e3) / inprod(e3,e3)) * e3
print(e1,',',e2,',',e3,',',e4)

🔗

As in Example 37.4.7, we'll clear fractions, which doesn't affect orthogonality.

xxxxxxxxxx
 
e2 = 2 * e2
e3 = 3 * e3
e4 = 2 * e4
print(e1,',',e2,',',e3,',',e4)

🔗

That's better.

🔗

Let's double-check orthogonality. This time we'll use some Python list magic.

xxxxxxxxxx
 
E = [e1,e2,e3,e4]
for e in E:
    print("checking",e,"against all")
    for f in E:
        print(inprod(e,f))

🔗

Each block of four results contains three zeros and one non-zero value. The three zeros represents the orthogonality of that particular vector with the other three vectors in the basis, while the non-zero value is the result of the inner product of that vector paired with itself.

🔗

Subsection B.6.3 A polynomial example with an integral inner product

🔗

Here we will carry out Example 37.4.4.

🔗

First let's enter our initial (non-orthogonal) basis vectors into Sage.

xxxxxxxxxx
 
v1 = 1
v2 = x
v3 = x^2
print(v1,',',v2,',',v3)

🔗

To make life easier, let's create a Python function to compute our inner product values for us. The last line is to test it out.

xxxxxxxxxx
 
def inprod(f,g):
    return integrate(f*g,x,0,1)
inprod(1,x)

🔗

Our calculation inprod(1,x) should be computing the integral

$\begin{equation*} \int_0^1 1 x \, dx \text{,} \end{equation*}$

🔗

which should evaluate to the area of a triangle of base and height both

$1\text{.}$ Looks like the computation got it correct.

🔗

Now we carry out the Gram-Schmidt process.

xxxxxxxxxx
 
e1 = v1
e2 = v2 - (inprod(v2,e1) / inprod(e1,e1)) * e1
e3 = v3 - (inprod(v3,e1) / inprod(e1,e1)) * e1 - (inprod(v3,e2) / inprod(e2,e2)) * e2
print('e1 =',e1)
print('e2 =',e2)
print('e3 =',e3)

🔗

We now have the orthogonal basis

$\begin{equation*} \left\{ 1, x - \frac{1}{2}, x^2 - x + \frac{1}{6} \right\} \text{,} \end{equation*}$

🔗

just as in Example 37.4.4. Let's double-check:

xxxxxxxxxx
 
print(inprod(e1,e2))
print(inprod(e1,e3))
print(inprod(e2,e3))

🔗

Yep, all zeros.

🔗

If we want an orthonormal basis, just need to normalize.

xxxxxxxxxx
 
f1 = e1 / sqrt(inprod(e1,e1))
f2 = e2 / sqrt(inprod(e2,e2))
f3 = e3 / sqrt(inprod(e3,e3))
print(f1,',',f2,',',f3)

🔗

Subsection B.6.4 A polynomial example with a “sampling” inner product

🔗

Now we'll repeat the Gram-Schmidt process on the standard basis of

$\poly_2(\R)\text{,}$ but using the inner product

$\begin{equation*} \inprod{p}{q} = p(-1) q(-1) + p(0) q(0) + p(1) q(1) \text{.} \end{equation*}$

🔗

Once again we create a Python function to return inner product values, but this time we'll use slightly more sophisticated Python/Sage programming techniques.

xxxxxxxxxx
 
def inprod(p,q):
    inprod_value = 0
    for x_value in [-1, 0, 1]:
        inprod_value += p.subs(x = x_value) * q.subs(x = x_value)
    return inprod_value
print("inprod function loaded")

🔗

Now let's enter our initial (non-orthogonal) basis vectors into Sage.

xxxxxxxxxx
 
v1 = 1
v2 = x
v3 = x^2
print(v1,',',v2,',',v3)

🔗

Now we carry out the Gram-Schmidt process.

xxxxxxxxxx
 
e1 = v1
e2 = v2 - (inprod(v2,e1) / inprod(e1,e1)) * e1
e3 = v3 - (inprod(v3,e1) / inprod(e1,e1)) * e1 - (inprod(v3,e2) / inprod(e2,e2)) * e2
print('e1 =',e1)
print('e2 =',e2)
print('e3 =',e3)

🔗

Even though we are using the same initial basis as in Subsection B.6.3, the Gram-Schmidt process has produced a different result because we are using a different inner product.

🔗

Notice that this time it turns out that

$1$ and

$x$ were already orthogonal, and we just needed to “straighten out”

$x^2\text{.}$

🔗

Finally, let's double-check orthogonality.

xxxxxxxxxx
 
print(inprod(e1,e2))
print(inprod(e1,e3))
print(inprod(e2,e3))

🔗

Yep, all zeros.

🔗

Subsection B.6.5 Bonus Fun

🔗

If you've taken a Python programming class, see if you can write a procedure (using def) that takes a list B and a function inprod as parameters, and returns a list containing the results of performing the Gram-Schmidt process on the objects in B using the input parameter function inprod as the inner product function. Fun!