Gram-Schmidt Orthonormalisation

Turn the columns of a matrix into an orthonormal basis for the same column space.

Process the columns left to right. For each one, subtract off its projection onto every basis vector already accepted, then normalise what is left:

v = a_j - sum over accepted q of (q . a_j) q
q_new = v / ||v||

If ||v|| is below 1e-10 the column added nothing new, so skip it rather than dividing by a near-zero. That is what makes the output a basis for the column space rather than a set of d vectors, so a matrix with dependent columns returns fewer columns than it started with.

Return an (n, r) array whose columns are the basis, in the order the original columns were accepted.

Use the modified form: subtract each projection from the running vector one at a time, updating v as you go, rather than computing all the projections against the original column and subtracting them at the end. The two are algebraically identical, but the classical version loses orthogonality badly in floating point once the columns are nearly dependent.

Function signature

def gram_schmidt(A):

Example

Input

A = np.array([[1.0, 1.0], [1.0, 0.0], [0.0, 0.0]])
gram_schmidt(A)

Output

array([[ 0.70710678,  0.70710678],
       [ 0.70710678, -0.70710678],
       [ 0.        ,  0.        ]])

Reasoning

The first column (1, 1, 0) normalises to (1, 1, 0) / sqrt(2). The second column (1, 0, 0) minus its projection on that is (0.5, -0.5, 0), which normalises to (1, -1, 0) / sqrt(2).

Constraints and notes

• Skip a column when what is left has norm below 1e-10.
• Use the modified form, which subtracts one projection at a time.
• The result can have fewer columns than the input.