Problem 24, Section 6.4

This is supposed to serve as an illustration for the kinds of techniques that would be useful to you in problem 25 from the same section. I'll be following along the textbook's hints, with computations and comments.

In [1]:
from sympy import *
init_printing(use_latex='mathjax')

I'm supposed to obtain an orthogonal basis for a matrix $A$ via Gram-Schmidt, as in example 2 of the same section. First, let's set down $A$.

In [3]:
A=Matrix([[-10,13,7,-11],
          [2,1,-5,3],
          [-6,3,13,-3],
          [16,-16,-2,5],       
          [2,1,-5,-7]])
A
Out[3]:
$\displaystyle \left[\begin{matrix}-10 & 13 & 7 & -11\\2 & 1 & -5 & 3\\-6 & 3 & 13 & -3\\16 & -16 & -2 & 5\\2 & 1 & -5 & -7\end{matrix}\right]$

Luckily for us, Gram-Schmidt is built right into SymPy! It works by passing a list of column vectors to the GramSchmidt function. Since the columnspace method produces exactly such a list out of the columns of our matrix, let's chain the two functions.

In [5]:
GramSchmidt(A.columnspace())
Out[5]:
$\displaystyle \left[ \left[\begin{matrix}-10\\2\\-6\\16\\2\end{matrix}\right], \ \left[\begin{matrix}3\\3\\-3\\0\\3\end{matrix}\right], \ \left[\begin{matrix}6\\0\\6\\6\\0\end{matrix}\right], \ \left[\begin{matrix}0\\5\\0\\0\\-5\end{matrix}\right]\right]$

That should be it: the vectors displayed above are supposed to be orthogonal, if GramSchmidt has done its job. Note though that the vectors are not orthonormal, just orthogonal: they're not normalized so that their lengths are $1$. Although the problem doesn't ask for this, you can normalize them by passing an extra option to GramSchmidt:

In [6]:
GramSchmidt(A.columnspace(), orthonormal=True)
Out[6]:
$\displaystyle \left[ \left[\begin{matrix}- \frac{1}{2}\\\frac{1}{10}\\- \frac{3}{10}\\\frac{4}{5}\\\frac{1}{10}\end{matrix}\right], \ \left[\begin{matrix}\frac{1}{2}\\\frac{1}{2}\\- \frac{1}{2}\\0\\\frac{1}{2}\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{3}}{3}\\0\\\frac{\sqrt{3}}{3}\\\frac{\sqrt{3}}{3}\\0\end{matrix}\right], \ \left[\begin{matrix}0\\\frac{\sqrt{2}}{2}\\0\\0\\- \frac{\sqrt{2}}{2}\end{matrix}\right]\right]$

See what happened there? You got the same vectors as before, but rescaled: the first one was scaled down by $20$, the second one got scaled down by $6$, etc.

You know what though? For good measure, why don't we check that these guys are indeed orthonormal? The way to do this is to

  • collect them as columns into a single matrix $M$,
  • and then simply check that $M^T \cdot M$ is the identity matrix $I_4$.

I'll use NumPy's hstack function to make my list of column vectors into a matrix. First, let's name that list of columns something:

In [8]:
cols=GramSchmidt(A.columnspace(), orthonormal=True)
cols
Out[8]:
$\displaystyle \left[ \left[\begin{matrix}- \frac{1}{2}\\\frac{1}{10}\\- \frac{3}{10}\\\frac{4}{5}\\\frac{1}{10}\end{matrix}\right], \ \left[\begin{matrix}\frac{1}{2}\\\frac{1}{2}\\- \frac{1}{2}\\0\\\frac{1}{2}\end{matrix}\right], \ \left[\begin{matrix}\frac{\sqrt{3}}{3}\\0\\\frac{\sqrt{3}}{3}\\\frac{\sqrt{3}}{3}\\0\end{matrix}\right], \ \left[\begin{matrix}0\\\frac{\sqrt{2}}{2}\\0\\0\\- \frac{\sqrt{2}}{2}\end{matrix}\right]\right]$

And then pass it off to NumPy:

In [9]:
import numpy as np
M=Matrix(np.hstack(cols))
M
Out[9]:
$\displaystyle \left[\begin{matrix}- \frac{1}{2} & \frac{1}{2} & \frac{\sqrt{3}}{3} & 0\\\frac{1}{10} & \frac{1}{2} & 0 & \frac{\sqrt{2}}{2}\\- \frac{3}{10} & - \frac{1}{2} & \frac{\sqrt{3}}{3} & 0\\\frac{4}{5} & 0 & \frac{\sqrt{3}}{3} & 0\\\frac{1}{10} & \frac{1}{2} & 0 & - \frac{\sqrt{2}}{2}\end{matrix}\right]$

All right! Now that it's a matrix, for the moment of truth:

In [12]:
M.T*M
Out[12]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$

One final thought: note that since $M$ is $5\times 4$ and hence not square, the fact that it is left invertible with left inverse $M^T$ means it must by necessity not be right-invertible. In particular, there's no way $M\cdot M^T$ can be the identity $I_5$:

In [13]:
M*M.T
Out[13]:
$\displaystyle \left[\begin{matrix}\frac{5}{6} & \frac{1}{5} & \frac{7}{30} & - \frac{1}{15} & \frac{1}{5}\\\frac{1}{5} & \frac{19}{25} & - \frac{7}{25} & \frac{2}{25} & - \frac{6}{25}\\\frac{7}{30} & - \frac{7}{25} & \frac{101}{150} & \frac{7}{75} & - \frac{7}{25}\\- \frac{1}{15} & \frac{2}{25} & \frac{7}{75} & \frac{73}{75} & \frac{2}{25}\\\frac{1}{5} & - \frac{6}{25} & - \frac{7}{25} & \frac{2}{25} & \frac{19}{25}\end{matrix}\right]$

In other words, for non-square matrices the fact that their columns are orthonormal doesn't mean their rows are orthonormal as well (in fact they definitely won't be).