# Gram-Schmidt and Modified Gram-Schmidt¶

In [1]:
#keep
import numpy as np
import numpy.linalg as la

In [2]:
#keep
A = np.random.randn(3, 3)

In [3]:
#keep
def test_orthogonality(Q):
print("Q:")
print(Q)

print("Q^T Q:")
QtQ = np.dot(Q.T, Q)
QtQ[np.abs(QtQ) < 1e-15] = 0
print(QtQ)

In [4]:
#keep
Q = np.zeros(A.shape)


Now let us generalize the process we used for three vectors earlier:

In [5]:
for k in range(A.shape[1]):
avec = A[:, k]
q = avec
for j in range(k):
q = q - np.dot(avec, Q[:,j])*Q[:,j]

Q[:, k] = q/la.norm(q)


This procedure is called Gram-Schmidt Orthonormalization.

In [6]:
#keep
test_orthogonality(Q)

Q:
[[-0.65609385 -0.75361761 -0.04001703]
[ 0.65802183 -0.59722164  0.4586214 ]
[ 0.36952419 -0.2745666  -0.88773028]]
Q^T Q:
[[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
[  0.00000000e+00   1.00000000e+00  -1.05471187e-15]
[  0.00000000e+00  -1.05471187e-15   1.00000000e+00]]


Now let us try a different example (Source):

In [7]:
#keep

np.set_printoptions(precision=13)

eps = 1e-8

A = np.array([
[1,  1,  1],
[eps,eps,0],
[eps,0,  eps]
])

A

Out[7]:
array([[  1.0000000000000e+00,   1.0000000000000e+00,   1.0000000000000e+00],
[  1.0000000000000e-08,   1.0000000000000e-08,   0.0000000000000e+00],
[  1.0000000000000e-08,   0.0000000000000e+00,   1.0000000000000e-08]])
In [8]:
#keep
Q = np.zeros(A.shape)

In [9]:
#keep
for k in range(A.shape[1]):
avec = A[:, k]
q = avec
for j in range(k):
print(q)
q = q - np.dot(avec, Q[:,j])*Q[:,j]

print(q)
q = q/la.norm(q)
Q[:, k] = q
print("norm -->", q)
print("-------")

[  1.0000000000000e+00   1.0000000000000e-08   1.0000000000000e-08]
norm --> [  1.0000000000000e+00   1.0000000000000e-08   1.0000000000000e-08]
-------
[  1.0000000000000e+00   1.0000000000000e-08   0.0000000000000e+00]
[  0.0000000000000e+00   0.0000000000000e+00  -1.0000000000000e-08]
norm --> [ 0.  0. -1.]
-------
[  1.0000000000000e+00   0.0000000000000e+00   1.0000000000000e-08]
[  0.0000000000000e+00  -1.0000000000000e-08   0.0000000000000e+00]
[  0.0000000000000e+00  -1.0000000000000e-08  -1.0000000000000e-08]
norm --> [ 0.              -0.7071067811865 -0.7071067811865]
-------

In [10]:
#keep
test_orthogonality(Q)

Q:
[[  1.0000000000000e+00   0.0000000000000e+00   0.0000000000000e+00]
[  1.0000000000000e-08   0.0000000000000e+00  -7.0710678118655e-01]
[  1.0000000000000e-08  -1.0000000000000e+00  -7.0710678118655e-01]]
Q^T Q:
[[  1.0000000000000e+00  -1.0000000000000e-08  -1.4142135623731e-08]
[ -1.0000000000000e-08   1.0000000000000e+00   7.0710678118655e-01]
[ -1.4142135623731e-08   7.0710678118655e-01   1.0000000000000e+00]]


Questions:

• What happened?
• How do we fix it?
In [11]:
#keep
Q = np.zeros(A.shape)

In [12]:
for k in range(A.shape[1]):
q = A[:, k]
for j in range(k):
q = q - np.dot(q, Q[:,j])*Q[:,j]

Q[:, k] = q/la.norm(q)

In [13]:
#keep
test_orthogonality(Q)

Q:
[[  1.0000000000000e+00   0.0000000000000e+00   0.0000000000000e+00]
[  1.0000000000000e-08   0.0000000000000e+00  -1.0000000000000e+00]
[  1.0000000000000e-08  -1.0000000000000e+00   0.0000000000000e+00]]
Q^T Q:
[[  1.0000000000000e+00  -1.0000000000000e-08  -1.0000000000000e-08]
[ -1.0000000000000e-08   1.0000000000000e+00   0.0000000000000e+00]
[ -1.0000000000000e-08   0.0000000000000e+00   1.0000000000000e+00]]


This procedure is called Modified Gram-Schmidt Orthogonalization.

Questions:

• Is there a difference mathematically between modified and unmodified?
• Why are there $10^{-8}$ values left in $Q^TQ$?
In [16]:
A = np.random.rand(10,10)
print(np.linalg.cond(A))
print(np.linalg.cond(A)**2)
print(np.linalg.cond(A.T.dot(A)))

53.4739304307
2859.4612357
2859.4612357

In [ ]:


In [ ]: