Gram-Schmidt and Modified Gram-Schmidt

import numpy as np
import numpy.linalg as la

A = np.random.randn(3, 3)

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)

Q = np.zeros(A.shape)

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

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.

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):

np.set_printoptions(precision=13)

eps = 1e-8

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

A

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]])

Q = np.zeros(A.shape)

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]
-------

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?

Q = np.zeros(A.shape)

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)

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$?