Classical Gram-Schmidt algorithm

gistfile1.py

import numpy as np
from scipy import linalg

def cgs(A):
    """Classical Gram-Schmidt (CGS) algorithm"""
    m, n = A.shape
    R = np.zeros((n, n))
    Q = np.empty((m, n))
    R[0, 0] = linalg.norm(A[:, 0])
    Q[:, 0] = A[:, 0] / R[0, 0]
    for k in range(1, n):
        R[:k-1, k] = np.dot(Q[:m, :k-1].T, A[:m, k])
        z = A[:m, k] - np.dot(Q[:m, :k-1], R[:k-1, k])
        R[k, k] = linalg.norm(z) ** 2
        Q[:m, k] = z / R[k, k]
    return Q, R

if __name__ == '__main__':

    n = 5
    X = np.random.random((n, n))
    import rogues
#    X = rogues.hilb(n)
    Q, R = cgs(X)
    assert np.allclose(np.dot(Q, R), X)
    print linalg.norm(np.dot(Q.T, Q) - np.eye(5), np.inf)