Gram-Schmidt and Modified Gram-Schmidt¶

In [3]:
import numpy as np
import numpy.linalg as la
In [4]:
A = np.random.randn(3, 3)
In [5]:
def test_orthogonality(Q):
    print("Q:")
    print(Q.round(9))
    
    print("Q^T Q:")
    QtQ = np.dot(Q.T, Q)
    QtQ[np.abs(QtQ) < 1e-15] = 0
    print(QtQ.round(9))
In [6]:
Q = np.zeros(A.shape)

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

In [7]:
for k in range(A.shape[1]):
    avec = A[:, k]
    q = avec
    for j in range(k):
        coeff = avec @ Q[:,j]
        q = q - coeff*Q[:,j]
    
    Q[:, k] = q/la.norm(q)

This procedure is called Gram-Schmidt Orthonormalization.

In [8]:
test_orthogonality(Q)

Q:
[[-0.58228521  0.04686309 -0.81163279]
 [-0.70728599  0.4630539   0.53416067]
 [ 0.40086215  0.88509035 -0.23648384]]
Q^T Q:
[[ 1.  0. -0.]
 [ 0.  1.  0.]
 [-0.  0.  1.]]

Now let us try a different example (Source):

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

eps = 1e-8

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

A
Out[15]:
array([[1.e+00, 1.e+00, 1.e+00],
       [1.e-08, 1.e-08, 0.e+00],
       [1.e-08, 0.e+00, 1.e-08]])
In [16]:
Q = np.zeros(A.shape)
In [22]:
for k in range(A.shape[1]):
    avec = A[:, k]
    q = avec
    for j in range(k):
        coeff = avec @ Q[:, j]
        print(q, coeff)
        q = q - coeff * Q[:,j]
    
    print(q)
    Q[:, k] = q/la.norm(q)
    print("norm -->", Q[:, k])
    print("-------")

[1.e+00 1.e-08 1.e-08]
norm --> [1.e+00 1.e-08 1.e-08]
-------
[1.e+00 1.e-08 0.e+00] 1.0
[ 0.e+00  0.e+00 -1.e-08]
norm --> [ 0.  0. -1.]
-------
[1.e+00 0.e+00 1.e-08] 1.0
[ 0.e+00 -1.e-08  0.e+00] -1e-08
[ 0.e+00 -1.e-08 -1.e-08]
norm --> [ 0.              -0.7071067811865 -0.7071067811865]
-------
In [18]:
test_orthogonality(Q)

Q:
[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 1.00000000e-08  0.00000000e+00 -7.07106781e-01]
 [ 1.00000000e-08 -1.00000000e+00 -7.07106781e-01]]
Q^T Q:
[[ 1.00000000e+00 -1.00000000e-08 -1.40000000e-08]
 [-1.00000000e-08  1.00000000e+00  7.07106781e-01]
 [-1.40000000e-08  7.07106781e-01  1.00000000e+00]]

Questions:

  • What happened?
  • How do we fix it?
In [19]:
Q = np.zeros(A.shape)
In [21]:
for k in range(A.shape[1]):
    q = A[:, k]
    for j in range(k):
        coeff = q @ Q[:,j]
        print(q, coeff)
        q = q - np.dot(q, Q[:,j])*Q[:,j]

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

[1.e+00 1.e-08 1.e-08]
norm --> [1.e+00 1.e-08 1.e-08]
-------
[1.e+00 1.e-08 0.e+00] 1.0
[ 0.e+00  0.e+00 -1.e-08]
norm --> [ 0.  0. -1.]
-------
[1.e+00 0.e+00 1.e-08] 1.0
[ 0.e+00 -1.e-08  0.e+00] 0.0
[ 0.e+00 -1.e-08  0.e+00]
norm --> [ 0. -1.  0.]
-------
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$?