In [1]:
from numpy import *
from mth309 import *

QR factorization via Gram-Schmidt

Example of Gram-Schmidt orthogonalization from last time

Three random (non-orthogonal, but LI) vectors in $\mathbb{R}^4$ (columns of X):

In [2]:
A = Matrix([[7,10,12,2],[-1,4,3,6],[1,2,3,0]]).T
A
Out[2]:
  7 -1 1
 10  4 2
 12  3 3
  2  6 0

Give names to the columns:

In [3]:
x1,x2,x3 = A.T
In [4]:
x1
Out[4]:
7 10 12 2
In [5]:
x2
Out[5]:
-1 4 3 6
In [6]:
x3
Out[6]:
1 2 3 0

The first vector in our orthogonal basis:

In [7]:
u1 = x1

The second:

In [8]:
u2 = x2 - dot(u1,x2)/dot(u1,u1)*u1
u2
Out[8]:
-32/11 14/11 -3/11 60/11

Check orthogonality of u1,u2:

In [9]:
dot(u1,u2)
Out[9]:
Fraction(0, 1)

The third orthogonal basis vector:

In [10]:
u3 = x3 - dot(u1,x3)/dot(u1,u1)*u1 - dot(u2,x3)/dot(u2,u2)*u2
u3
Out[10]:
-752/1317 -110/1317 196/439 -346/1317

Check orthogonalities:

In [11]:
dot(u1,u3)
Out[11]:
Fraction(0, 1)
In [12]:
dot(u2,u3)
Out[12]:
Fraction(0, 1)
In [13]:
Q = vstack([u1,u2,u3]).astype(float).T
Q
Out[13]:
array([[  7.        ,  -2.90909091,  -0.57099468],
       [ 10.        ,   1.27272727,  -0.08352316],
       [ 12.        ,  -0.27272727,   0.44646925],
       [  2.        ,   5.45454545,  -0.2627183 ]])

Not orthonormal.

In [14]:
Q /= linalg.norm(Q,axis=0)
Q
Out[14]:
array([[ 0.4061812 , -0.46049124, -0.73631282],
       [ 0.58025885,  0.20146492, -0.10770533],
       [ 0.69631062, -0.04317105,  0.57573396],
       [ 0.11605177,  0.86342108, -0.33878223]])
In [15]:
dot(Q.T,Q).round(15)
Out[15]:
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
In [20]:
R = dot(Q.T,A.astype(float))
R.round(15)
Out[20]:
array([[ 17.23368794,   4.70009671,   3.65563078],
       [ -0.        ,   6.31736424,  -0.18707457],
       [ -0.        ,  -0.        ,   0.7754784 ]])

Verify that QR = A

In [22]:
QR =dot(Q,R)
QR.round(15)
Out[22]:
array([[  7.,  -1.,   1.],
       [ 10.,   4.,   2.],
       [ 12.,   3.,   3.],
       [  2.,   6.,   0.]])

A recovered!