# Gram-Schmidt and Modified Gram-Schmidt¶

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

In [35]:
A = np.random.randn(3, 3)

In [36]:
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 [37]:
Q = np.zeros(A.shape)


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

In [38]:
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 [39]:
test_orthogonality(Q)

Q:
[[-0.6932589320501 -0.6855758663147 -0.2222111263183]
[-0.7199408381809  0.6447564063972  0.2568547564853]
[ 0.032821374928  -0.3380457187077  0.9405572015626]]
Q^T Q:
[[ 1.  0.  0.]
[ 0.  1.  0.]
[ 0.  0.  1.]]


Now let us try a different example (Source):

In [2]:
np.set_printoptions(precision=13)

eps = 1e-8

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

A

Out[2]:
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 [3]:
Q = np.zeros(A.shape)

In [17]:
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 [43]:
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 [44]:
Q = np.zeros(A.shape)

In [47]:
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 [48]:
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$?