Skip to content

Gram-Schmidt Orthogonalization and Projection

from sympy import *

Gram-Schmidt Procedure

Suppose we have a subspace \(\mathbb{S}\) of \(\mathbb{R}^n\) whose basis consists of \(k\) vectors \(\vec{v}_1,\vec{v}_2, \ldots , \vec{v}_k\).

\[ \mathbb{S} = \operatorname{span}\left\{ \vec{v}_1,\vec{v}_2, \ldots , \vec{v}_k \right\}. \]

The Gram-Schmidt procedure computes from this basis a new basis for \(\mathbb{S}\) that is orthonormal (that is, the vectors of of unit length and mutually orthogonal). That is,

\[ \operatorname{gram\_schmidt}\left(\left\{ \vec{v}_1,\vec{v}_2, \ldots , \vec{v}_k \right\}\right) \rightarrow \left\{ \vec{w}_1,\vec{w}_2, \ldots , \vec{w}_k \right\} \]

where

  • \(|\vec{w}_i| = 1\), for all \(i\),
  • \(\vec{w}_i \cdot \vec{w}_j = 0\) for all \(i\ne j\).

The procedure usually constructs the vectors by, setting \(\vec{w}_1 = \vec{v}_1/|\vec{v}_1|\), and then recursively building the rest of the vectors by

\[ \vec{w}_j = \operatorname{perp}_{\operatorname{span}\{\vec{w}_1, \ldots, \vec{w}_{j-1}\}}(\vec{v}_j), \quad 1< j \le k. \]

It follows that, for each \(1\le i \le k\)

\[ \operatorname{span} \left\{ \vec{v}_1, \vec{v}_2 \ldots, \vec{v}_i \right\} = \operatorname{span} \left\{ \vec{w}_1, \vec{w}_2 \ldots, \vec{w}_i \right\} \]

There is a nice Gram-Schmidt orthogonalizer which will take a list of vectors and orthogonalize them with respect to another. The vectors in the list must be linearly independent. There is an optional argument which specifies whether or not the output should also be normalized, it defaults to False.

Example

Use the Gram-Schmidt procedure on the set

\[\left\{ \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} , \quad \begin{bmatrix} -1 \\ 2 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix} \right\} \]
A = [Matrix([1, 1, 0]), Matrix([-1, 2, 1]), Matrix([0, 1, 1])]
pprint(GramSchmidt(A))

\(\left[ \left[\begin{matrix} 1\\1\\0 \end{matrix} \right], \left[\begin{matrix} -3/2\\3/2\\1 \end{matrix} \right], \left[\begin{matrix} 2/11\\-2/11\\6/11 \end{matrix} \right] \right]\)

Use optional argument True to have a normalized basis returned.

A = [Matrix([1, 1, 0]), Matrix([-1, 2, 1]), Matrix([0, 1, 1])]
pprint(GramSchmidt(A, True))

\(\left[ \left[\begin{matrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ 0 \end{matrix} \right], \left[\begin{matrix} \frac{-3}{\sqrt{22}}\\ \frac{3}{\sqrt{22}}\\ \frac{2}{\sqrt{22}} \end{matrix} \right], \left[\begin{matrix} \frac{1}{\sqrt{11}}\\ \frac{-1}{\sqrt{11}}\\ \frac{3}{\sqrt{11}} \end{matrix} \right] \right]\)

If we have a matrix of column vectors we wish to orthogonalize, first extract them to a list.

A = Matrix([[1, 1, 0], [-1, 2, 1], [0, 1, 1]]).T
vectorList = [A.col(i) for i in range(A.shape[1])]
GramSchmidt(vectorList)

Note: we can check that the vectors returned by GramSchmidt are orthonormal. This would mean that the matrix of these vectors must be orthogonal. Recall, a matrix \(Q\) is orthogonal if \(Q^TQ = I\).

A = [Matrix([1, 1, 0]), Matrix([-1, 2, 1]), Matrix([0, 1, 1])]
L = GramSchmidt(A,True)
B = L[0]
for col in L[1:]:
    B = B.row_join(col)
B.T*B

\(\left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right]\)

Example

Use the Gram-Schmidt procedure on the set

\[\left\{ \begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 1 \\ 0 \\ 1 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 2 \\ 1 \\ 1 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \end{bmatrix} \right\} \]

Note: the third vector is in the span of the first two. Since GramSchmidt requires the list of vectors to be linearly independent we'll use rref to extract the linearly independent columns.

A = Matrix([[1, 1, 0, 1],[1, 0, 1, 1],[2, 1, 1, 2],[1, 0, 0, 1]]).T
A

\(\left[\begin{matrix} 1 & 1 & 2 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 1 & 1 & 2 & 1 \\ \end{matrix} \right]\)

pivots = A.rref()[1]
L = []
for i in pivots:
    L.append(A[:,i])
pprint(L)

\(\left[ \left[\begin{matrix} 1 \\ 1 \\ 0 \\ 1 \end{matrix}\right], \left[\begin{matrix} 1 \\ 0 \\ 1 \\ 1 \end{matrix}\right], \left[\begin{matrix} 1 \\ 0 \\ 0 \\ 1 \end{matrix}\right] \right]\)

Now we have a list of linearly independent vectors from the original list. We can apply GramSchmidt to this list.

pprint(GramSchmidt(L)) # (1)
  1. use optional second argument True to return a set of normalize vectors

\(\left[ \left[ \begin{matrix} 1 \\ 1 \\ 0 \\ 1 \end{matrix} \right], \left[ \begin{matrix} 1/3 \\ -2/3 \\ 1 \\ 1/3 \end{matrix} \right], \left[ \begin{matrix} 1/5 \\ -2/5 \\ -2/5 \\ 1/5 \end{matrix} \right] \right]\)

Gram-Schmidt via QR Decomposition

The Gram-Schmidt can be found hidden in the SymPy's QR decomposition. We'll cover a number of matrix decompositions methods in the next section, but for now we note that QR decomposition is a factorization of the form

\[A = QR\]

where \(Q\) is orthogonal (i.e. columns are orthonormal), and \(R\) is upper triangular. This means that column \(i\) of \(A\) is a linear combination of the first \(i\) columns of \(Q\). Column \(i\) of \(R\) tells us exactly what this linear combination is. This is true for each column \(i\).

The function sympy.QRdecomposition returns the two matrices \(Q\) and \(R\) associated with matrix \(A\).

Let's redo the two examples above.

Example

Use sympy.QRdecomposition to find the Gram-Schmidt basis for the set

\[\left\{ \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} , \quad \begin{bmatrix} -1 \\ 2 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix} \right\} \]

First form a column matrix, then return the matrix \(Q\) from QR decomposition.

A = Matrix([[1, 1, 0],[-1, 2, 1],[0, 1, 1]]).T
Q, R = A.QRdecomposition()
Q

\(\left[\begin{matrix}\frac{\sqrt{2}}{2} & - \frac{3 \sqrt{22}}{22} & \frac{\sqrt{11}}{11}\\\frac{\sqrt{2}}{2} & \frac{3 \sqrt{22}}{22} & - \frac{\sqrt{11}}{11}\\0 & \frac{\sqrt{22}}{11} & \frac{3 \sqrt{11}}{11}\end{matrix}\right]\)

This is the same collection of orthonormal vectors we found above (up to scaling by \(-1\)).

Here is the example we looked at above, where the original set of vectors is linearly dependent.

Example

Use sympy.QRdecomposition to find the Gram-Schmidt basis for the set

\[\left\{ \begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 1 \\ 0 \\ 1 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 2 \\ 1 \\ 1 \\ 1 \end{bmatrix} , \quad \begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \end{bmatrix} \right\} \]

First form a column matrix, then return the matrix \(Q\) from QR decomposition.

A = Matrix([[1, 1, 0, 1],[1, 0, 1, 1],[2, 1, 1, 2],[1, 0, 0, 1]]).T
Q, R = A.QRdecomposition()
Q

\(\left[\begin{matrix}\frac{\sqrt{3}}{3} & \frac{\sqrt{15}}{15} & \frac{\sqrt{10}}{10}\\\frac{\sqrt{3}}{3} & - \frac{2 \sqrt{15}}{15} & - \frac{\sqrt{10}}{5}\\0 & \frac{\sqrt{15}}{5} & - \frac{\sqrt{10}}{5}\\\frac{\sqrt{3}}{3} & \frac{\sqrt{15}}{15} & \frac{\sqrt{10}}{10}\end{matrix}\right]\)

This is the same collection of vectors we found above, they are just normalized here..

Projection onto a Subspace

Suppose we have an orthogonal basis \(\{ \vec{w}_1, \vec{w}_2, \ldots, \vec{w}_k \}\) for a subspace \(\mathbb{S}\). The the projection of a vector \(\vec{v}\) onto \(\mathbb{S}\) is given by:

\[ \operatorname{proj}_{\mathbb{S}}(\vec{v}) = \frac{\vec{v}\cdot\vec{w}_1}{\vec{w}_1\cdot\vec{w}_1}\vec{w}_1 + \frac{\vec{v}\cdot\vec{w}_2}{\vec{w}_2\cdot\vec{w}_2}\vec{w}_2 + \cdots + \frac{\vec{v}\cdot\vec{w}_k}{\vec{w}_k\cdot\vec{w}_k}\vec{w}_k \]

This is used in the Gram-Schmidt procedure. Here we just implement the projection and perpendicular as functions in python so we can use them as needed.

Projection onto a Subspace
def proj(v, A):
    '''input : v is a vector, and A is a matrix whose
               column vectors are form an orthogonal basis.
       output: the projection of v onto the subspace spanned
               columns of A.'''
    num_rows = A.shape[0]
    num_cols = A.shape[1]
    w = zeros(num_rows)
    for i in range(num_cols):
        coli = A[:,i]
        w += (v.dot(coli)/coli.dot(coli))*coli
    return w

def perp(v,A):
    '''input : v is a vector, and A is a matrix whose
               column vectors are form an orthogonal basis.
       output: the perpendicular of the projection of v onto
               the subspace spanned columns of A.'''
    return v - proj(v,A)

Example

Let \(\mathbb{S} = \operatorname{span}\left\{ \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}, \begin{bmatrix} 1 \\ -1 \\ 1 \\ -1\end{bmatrix} \right\}\) and let \(\vec{x} = \begin{bmatrix} 2\\5\\-7\\3 \end{bmatrix}\). Determine \(\operatorname{proj}_{\mathbb{S}}(\vec{x})\) and \(\operatorname{perp}_{\mathbb{S}}(\vec{x})\).

Note: Basis is already orthogonal, so we can project directly onto it.

A = Matrix([[1.,1,1,1],[1,-1,1,-1]]).T
v = Matrix([2,5,-7,3])
print('projection of v onto columnspace of A: ')
pprint(proj(v,A))
print('perpendicular of v onto columnspace of A: ')
pprint(perp(v,A))

\(\text{projection of v onto columnspace of A:}\)

\(\left[\begin{matrix} -2.5 \\ 4 \\ -2.5 \\ 4 \end{matrix}\right]\)

\(\text{perpendicular of v onto columnspace of A:}\)

\(\left[\begin{matrix} 4.5 \\ 1 \\ -4.5 \\ -1\end{matrix}\right]\)

In the next example, we use the Gram-Schmidt procedure to first find an orthogonal basis for the subspace, then we project the vector onto this.

Example

Project \(\vec{v} = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\) onto the space with basis:

\[ \vec{v}_1 = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} , \quad \vec{v}_2= \begin{bmatrix} -1 \\ 2 \\ 1 \end{bmatrix} \]

The first thing to do would be to find an orthogonal basis, we could use Gram-Schmidt to do this (see above), but here we didn't normalize:

\[ \vec{w}_1 = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} , \quad \vec{w}_2= \begin{bmatrix} -3 \\ 3 \\ 2 \end{bmatrix} \]

Then projecting \(\vec{x}\) onto this basis we get

\[ \operatorname{proj}_{\mathbb{S}}(\vec{v}) = \frac{\vec{v}\cdot\vec{w}_1}{\vec{w}_1\cdot\vec{w}_1}\vec{w}_1 + \frac{\vec{v}\cdot\vec{w}_2}{\vec{w}_2\cdot\vec{w}_2}\vec{w}_2 \]
\[ = \begin{bmatrix} 2/2 - 6/22 \\ 2/2 + 6/22 \\ 2/11 \end{bmatrix} = \begin{bmatrix} 8/11 \\ 14/11 \\ 2/11 \end{bmatrix} \]

Also,

\[ \operatorname{perp}_{\mathbb{S}}(\vec{v}) = \vec{v} - \operatorname{proj}_{\mathbb{S}}(\vec{v}) = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} - \begin{bmatrix} 8/11 \\ 14/11 \\ 2/11 \end{bmatrix} = \begin{bmatrix} 3/11 \\ -3/11 \\ 9/11 \end{bmatrix} \]

Let's check this in python:

v1 = Matrix([1,1,0])
v2 = Matrix([-1,2,1])
A = [v1,v2]
Q = GramSchmidt(A)
pprint(Q)

\(\left[ \left[\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}\right], \left[\begin{matrix} -3/2 \\ 3/2 \\ 1 \end{matrix}\right] \right]\)

Q = Matrix(Q).reshape(2,3).T
Q

\(\left[\begin{matrix} 1 & -3/2 \\ 1 & 3/2 \\ 0 & 1 \\ \end{matrix} \right]\)

v = Matrix([1,1,1])
print('projection of v onto columnspace of Q:')
pprint(proj(v,Q))

\(\text{projection of v onto columnspace of Q:}\)

\(\left[\begin{matrix} 8/11 \\ 14/11 \\ 2/11 \end{matrix}\right]\)

print('perpendicular of v onto columnspace of Q:')
pprint(perp(v,Q))

\(\text{perpendicular of v onto columnspace of Q:}\)

\(\left[\begin{matrix} 3/11 \\ -3/11 \\ 9/11 \end{matrix}\right]\)