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!