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]:
#clear

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]:
#clear

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