Notebook

Topic 7. Gram-Schmidt Orthogonalization¶

Author: Jaron Ellingson, jaronce@byu.edu

Introduction¶

Gram-Schmit Orthogonalization (sometimes refered to as the Gram-Schmidt Process) takes a set of vectors and finds an orthonormal basis which spans that set of vectors. The importance of the Gram-Schmidt Process lies in the fact that sometimes it is convienent to have an orthonormal basis for a set of vectors. Once you have an orthonormal basis, you can find the fourier decomposition, and other series representations such as spherical harmonics and Bessel functions. In addition, this process can be applied to other norm and inner product spaces such as $\mathcal{L}_{2}$ , which allows for finding an orthonormal set of functions.

Explanation of the theory¶

Given a set of of vectors $\textit{T}$ = $\{\textbf{p}_1,\textbf{p}_2,\dotso,\textbf{p}_{n}\}$ we want to find set of vectos $\textit{T'}$ = $\{\textbf{e}_1,\textbf{e}_2,\dotso,\textbf{e}_{n'}\}$ with $n'\leq n$ so that

$\begin{equation} span\{\textbf{e}_1,\textbf{e}_2,\dotso,\textbf{e}_{n'}\} = span\{\textbf{p}_1,\textbf{p}_2,\dotso,\textbf{p}_{n}\} \end{equation}$

and

$\begin{equation} \left\langle \textbf{e}_i,\textbf{e}_j \right\rangle = \delta_{i,j} \end{equation}$ Assume that none of the

$\textbf{p}_i$ vectors are zero vectors.

The process will be developed stepwise. The norm $\left\lVert \cdot \right\Vert$ is the induced norm.

Normalize the first vector

$\begin{equation} \textbf{e}_1 = \frac{\textbf{p}_1}{\left\lVert \textbf{p}_1 \right\Vert} \end{equation}$

2. Compute the difference between the projection of $\textbf{p}_2$ onto $\textbf{q}_1$ and $\textbf{p}_2$ by the orthogonality theorem, this is orthogonal to $\textbf{p}_1$ :

$\begin{equation} \textbf{q}_2 = \textbf{p}_2 - \left\langle \textbf{p}_2,\textbf{e}_1 \right\rangle \textbf{e}_1 \end{equation}$ If

$\textbf{q}_2 = 0$ , then

$\textbf{e}_2 \in span(\textbf{e})_1$ and can be discarded; we well asume that such discards are done as necessary in what follows. If

$\textbf{q}_2 \neq 0$ ,then normalize

$\begin{equation} \textbf{e}_2 = \frac{\textbf{q}_2}{\left\lVert \textbf{q}_2 \right\Vert} \end{equation}$

These steps are shown in the figure below.

In [2]:

%matplotlib inline
import matplotlib.pyplot as plt

plt.grid(False)
plt.arrow(-1, -1, 2, 1.5, color="Red", head_width=.1, length_includes_head=True)
plt.arrow(-1, -1, 3, 0, color="Red", head_width=.1, length_includes_head=True)
plt.arrow(-1, -1, 2, 0, color="Green", head_width=.1, length_includes_head=True)
plt.arrow(1, -1, 0, 1.5, color="Black", head_width=.1, length_includes_head=True)
plt.arrow(-1, -1, 0.75, 0, color="Blue", head_width=.1, length_includes_head=True)
plt.arrow(-1, -1, 0, 0.75, color="Blue", head_width=.1, length_includes_head=True)
plt.text(-.5, -1.25,"$e_1$", color="Blue",size=14)
plt.text(-1.14, 0,"$e_2$", color="Blue",size=14)
plt.text(1.5, -1.25,"$p_1$", color="Red",size=14)
plt.text(.5, 0.5,"$p_2$", color="Red",size=14)
plt.text(1.1, 0.25,"$q_2$", color="Black",size=14)
#plt.text(1.1, 0.25,"$q_2 = p_2 - proj_{p_1}(p_2)$", color="Black",size=14)
plt.axis((-3, 3, -2, 2))
#plt.axis('off')

Out[2]:

(-3, 3, -2, 2)

At the next stage, a vector orthongonal to $\textbf{e}_1$ and $\textbf{e}_2$ obained from the error between $\textbf{p}_3$ and its projection onto span( $\textbf{e}_1,\textbf{e}_2$ ):

$\begin{equation} \textbf{q}_3 = \textbf{p}_3 - \left\langle \textbf{p}_3,\textbf{e}_1 \right\rangle \textbf{e}_1 - \left\langle \textbf{p}_3,\textbf{e}_2 \right\rangle \textbf{e}_2 \end{equation}$

This is normalized to produce $\textbf{e}_3$

$\begin{equation} \textbf{e}_3 = \frac{\textbf{q}_3}{\left\lVert \textbf{q}_3 \right\Vert} \end{equation}$

In [3]:

plt.grid(False)
plt.arrow(-1, 0, 1.5, 1, color="Red", head_width=.1, length_includes_head=True)
#plt.arrow(-1, -1, 3, 0, color="Red", head_width=.1, length_includes_head=True)
plt.arrow(-1, 0, 1.5, -1, color="Green", head_width=.1, length_includes_head=True)
plt.arrow(0.5, -1, 0, 2, color="Black", head_width=.1, length_includes_head=True)
plt.arrow(-1, 0, 0.75, 0, color="Blue", head_width=.1, length_includes_head=True)
plt.arrow(-1, 0, -0.5, -0.5, color="Blue", head_width=.1, length_includes_head=True)
plt.arrow(-1, 0, 0, 0.75, color="Blue", head_width=.1, length_includes_head=True)
plt.arrow(0.5, -1, 0, 0.75, color="Blue", head_width=.1, length_includes_head=True)
plt.text(-.5, .15,"$e_1$", color="Blue",size=14)
plt.text(-1.14, 1,"$e_2$", color="Blue",size=14)
plt.text(-1.8, -.7,"$e_3$", color="Blue",size=14)
#plt.text(1.5, -1.25,"$p_1$", color="Red",size=14)
plt.text(0, 1.1,"$p_3$", color="Red",size=14)
plt.text(0.6, 0.7,"$q_2$", color="Black",size=14)
#plt.text(1.1, 0.25,"$q_2 = p_2 - proj_{p_1}(p_2)$", color="Black",size=14)
plt.axis((-3, 3, -2, 2))
#plt.axis('off')

Out[3]:

(-3, 3, -2, 2)

Now proceed inductively: To form the next orthongoanl vecotr using $\textbf{p}_k$ , determine the component orthongonal to all previously found vectors.

$\begin{equation} \textbf{q}_k = \textbf{p}_k - \sum_{i=1}^{k-1} \left\langle \textbf{p}_k,\textbf{e}_i \right\rangle \textbf{e}_i \end{equation}$

and normalize

$\begin{equation} \textbf{e}_k = \frac{\textbf{q}_k}{\left\lVert \textbf{q}_k \right\Vert} \end{equation}$

Simple Numerical Examples¶

Example #1¶

Given a set of vectors (shown here as test and test2), this function calculates the gram schmidt orthonormalization of the set. This example prints out the final results from the function.

In [4]:

import numpy as np

def gram_schmidt(vectors):
    basis = []
    for v in vectors:
        w = v - np.sum( np.dot(v,b)*b  for b in basis )
        if (w > 1e-10).any():  
            basis.append(w/np.linalg.norm(w))
    return np.array(basis)

test = np.array([[3.0, 1.0], [2.0, 2.0]])
test2 = np.array([[1.0, 1.0, 0.0], [1.0, 3.0, 1.0], [2.0, -1.0, 1.0]])

print(np.array(gram_schmidt(test)))
print(np.array(gram_schmidt(test2)))

[[ 0.9486833   0.31622777]
 [-0.31622777  0.9486833 ]]
[[ 0.70710678  0.70710678  0.        ]
 [-0.57735027  0.57735027  0.57735027]
 [ 0.40824829 -0.40824829  0.81649658]]

Example #2¶

The gram schmidt method can also be used for QR decomposition (A = Q*R). Where Q (Q*transpose(Q) = Idenity) is a unitary matrix and R is an upper triangular matrix.

In [5]:

def gram_schmidt(A,verbose=False):
    m = A.shape[0]
    n = A.shape[1]
    R=np.zeros([m,n])
    Q=np.zeros([m,m])
    R[0,0]=np.linalg.norm(A[:,0])
    if R[0, 0] == 0:
        print ("cannot complete Gram_Schimidt decomposition")
    else:
        Q[:,0]=A[:,0]/R[0,0]
    for i in range(1,n):
        Q[:,i]=A[:,i]
        for j in range(i):
            R[j,i]=np.dot(Q[:,j].T,Q[:,i])
            Q[:,i]-=R[j,i]*Q[:,j]
        R[i,i]=np.linalg.norm(Q[:,i])
        if R[i,i]==0:
            print("cannot complete Gram_Schimidt decomposition")
        else:
            Q[:, i] = Q[:, i] / R[i, i]
    return  Q,R

A = np.random.randint(0,10000,(4,4))
print("A = ", A)

Q, R  = gram_schmidt(A)
print ("Q * R = ", np.dot(Q, R))

print("Q = ", Q)
print("R = ", R)
print("Q * Tranpose(Q) = ", np.dot(Q, np.transpose(Q)))

A =  [[7452 1050 8608 1157]
 [3652 3175 1384 6425]
 [8434  204 7797 3573]
 [7642 7533  566 9293]]
Q * R =  [[7452. 1050. 8608. 1157.]
 [3652. 3175. 1384. 6425.]
 [8434.  204. 7797. 3573.]
 [7642. 7533.  566. 9293.]]
Q =  [[ 0.52905388 -0.31446399  0.74287648 -0.26334183]
 [ 0.25927332  0.28476492  0.25073724  0.88815377]
 [ 0.59877085 -0.51828775 -0.58984596  0.15790194]
 [ 0.54254291  0.74256256 -0.19325401 -0.34190778]]
R =  [[14085.52192856  5587.82432055  9888.62569001  9459.20589068]
 [    0.          6061.93445714 -5933.59056943  6514.57147057]
 [    0.             0.          2033.29041635 -1432.93427866]
 [    0.             0.             0.          2788.53614384]]
Q * Tranpose(Q) =  [[ 1.00000000e+00  3.60822483e-16  1.17961196e-16 -1.52655666e-16]
 [ 3.60822483e-16  1.00000000e+00  4.44089210e-16 -2.77555756e-16]
 [ 1.17961196e-16  4.44089210e-16  1.00000000e+00 -2.01227923e-16]
 [-1.52655666e-16 -2.77555756e-16 -2.01227923e-16  1.00000000e+00]]

An Engineering Application¶

This example gives the Chebyshev polynomials as a function of cos. These polynomials are simular to the Legendre polynomials but start with the set of functions: $\{cos(x), cos(2x), ... cos(nx)\}$ . The Gram-Schmidt Orthonormalization will find the orthonormal functions which span this set.

In [6]:

def norm(vector,t):
    return np.sqrt(np.trapz(vector*vector,t))

def inner(vector1,vector2,t):
    return np.trapz(vector1*vector2,t)

N = 100000
t = np.linspace(-1,1,N)
p = np.zeros((N,5))
for kk in range(5):
    p[:,kk] = np.cos(t*kk)
    
e1 = p[:,0]/norm(p[:,0],t)
q2 = p[:,1]-inner(p[:,1],e1,t)*e1
e2 = q2/norm(q2,t)
q3 = p[:,2]-inner(p[:,2],e2,t)*e2-inner(p[:,2],e1,t)*e1
e3 = q3/norm(q3,t)
q4 = p[:,3]-inner(p[:,3],e3,t)*e3-inner(p[:,3],e2,t)*e2-inner(p[:,3],e1,t)*e1
e4 = q4/norm(q4,t)
q5 = p[:,4]-inner(p[:,4],e4,t)*e4-inner(p[:,4],e3,t)*e3-inner(p[:,4],e2,t)*e2-inner(p[:,4],e1,t)*e1
e5 = q5/norm(q5,t)

plt.plot(t,e1,label='e_1')
plt.plot(t,e2,label='e_2')
plt.plot(t,e3,label='e_3')
plt.plot(t,e4,label='e_4')
plt.plot(t,e5,label='e_5')
plt.grid()
plt.legend()

Out[6]:

<matplotlib.legend.Legend at 0x7f34d7526ac8>

The Chebyshev and Legendre polynomials are important in many physics applications. One major example is shown in quantum phyisics where these polynomials appear as solutions to the Schrodinger equation. These polynomials allow for calculations of the wavefunctions for the hydrogen atom and for determining electrical potential in spherical coordinates.

In [ ]: