Gram-Schmidt Orthogonization using Numpy

gram_schmidt.py

import numpy as np

def gs_cofficient(v1, v2):
    return np.dot(v2, v1) / np.dot(v1, v1)

def multiply(cofficient, v):
    return map((lambda x : x * cofficient), v)

def proj(v1, v2):
    return multiply(gs_cofficient(v1, v2) , v1)

def gs(X):
    Y = []
    for i in range(len(X)):
        temp_vec = X[i]
        for inY in Y :
            proj_vec = proj(inY, X[i])
            temp_vec = map(lambda x, y : x - y, temp_vec, proj_vec)
        Y.append(temp_vec)
    return Y

#Test data
test = np.array([[3.0, 1.0], [2.0, 2.0]])
test2 = np.array([[1.0, 1.0, 0.0], [1.0, 3.0, 1.0], [2.0, -1.0, 1.0]])
print np.array(gs(test))
print np.array(gs(test2))

Hi .
Thanks for this code.
I need some thing very similar to this .
Can you please help me .
i have a matrix A which i have to write down in the form
A=UDR

where U is a well-conditioned orthogonal matrix, D is
a diagonal matrix and R is a right tri-
angular matrix .

hi, thx for sharing. it's not normalized and use array form, so i wrote another one: https://github.com/whille/mylab/blob/master/algrithm/gram_schmidt.py

Hey,
A vectorized (faster) version:

def gs(X, row_vecs=True, norm = True):
    if not row_vecs:
        X = X.T
    Y = X[0:1,:].copy()
    for i in range(1, X.shape[0]):
        proj = np.diag((X[i,:].dot(Y.T)/np.linalg.norm(Y,axis=1)**2).flat).dot(Y)
        Y = np.vstack((Y, X[i,:] - proj.sum(0)))
    if norm:
        Y = np.diag(1/np.linalg.norm(Y,axis=1)).dot(Y)
    if row_vecs:
        return Y
    else:
        return Y.T

I think the fastest & easiest way to do this with NumPy is to use its built-in QR factorization:

def gram_schmidt_columns(X):
    Q, R = np.linalg.qr(X)
    return Q

After checking ingmarschuster implementation I was not satisfied with speed so I wrote it using generators. It normalizes by default and assumes the vectors are in rows. Morover, if a vector is already there it is omitted from the normalization

def gram_schmidt(vectors):
    basis = []
    for v in vectors:
        w = v - np.sum( np.dot(v,b)*b  for b in basis )
        if (w > 1e-10).any():  
            basis.append(w/np.linalg.norm(w))
    return np.array(basis)

@jakevdp

Every page I look at people keep pointing towards the qr function. However I cannot seem to get it to give me sane results. Let's take an example:

# A semi-interesting set of vectors
vectors = np.array(
    [[0, 0, 1],
     [0, 0, 1],
     [1, 1, 1],])

print gram_schmidt(vectors) # From JJGO's answer
#[[ 0.          0.          1.        ]
# [ 0.70710678  0.70710678  0.        ]]

print gram_schmidt_columns(vectors)
#[[ 0.  0. -1.]
# [-0.  1.  0.]
# [-1. -0.  0.]]

If the input requires some preprocessing step before applying qr then I wish people would mention this. Anyway, thanks @JJGO for posting some useful code.

@AnnanFay

Gramm-Schmidt process only works for non-singular matrices, i.e. your column and row vectors need to be linearly independent.

@NarekMargaryan that's true, but it is often difficult to know what the rank of a large matrix is. Gram Schmidt can be modified to allow singular matrices, where you discard the projections of a previously-calculated linearly dependent vector. In other words, the vectors calculated after finding a linear dependent vector can be assumed to be zeros. I wrote a similar script with vectorized code: https://gist.github.com/anmolkabra/b95b8e7fb7a6ff12ba5d120b6d9d1937.

The script takes care of singular matrices @AnnanFay. Although I don't stop when I hit a linearly dependent vector, it's not hard to do so.

I tested both codes with OP's matrices:

test = np.array([[3.0, 1.0], [2.0, 2.0]])
test2 = np.array([[1.0, 1.0, 0.0], [1.0, 3.0, 1.0], [2.0, -1.0, 1.0]])

In [41]: jjgo_numpy(test)
Out[41]:
array([[ 0.83205029, -0.5547002 ],
       [ 0.5547002 ,  0.83205029]])

In [42]: gs_by_qr(test)
Out[42]:
array([[-0.83205029, -0.5547002 ],
       [-0.5547002 ,  0.83205029]])

In [43]: jjgo_numpy(test2)
Out[43]:
array([[ 0.40824829,  0.20739034, -0.88900089],
       [ 0.40824829,  0.82956136,  0.38100038],
       [ 0.81649658, -0.51847585,  0.25400025]])

In [44]: gs_by_qr(test2)
Out[44]:
array([[-0.40824829, -0.20739034, -0.88900089],
       [-0.40824829, -0.82956136,  0.38100038],
       [-0.81649658,  0.51847585,  0.25400025]])

QR factorization (@jakevdp) and @JJGO implementation gives the same results, except for the sign.

Here is another example of how to implement this using numpy. The previous versions obscured the details of the algorithm.

import numpy as np


def gram_schmidt(A):
    """Orthogonalize a set of vectors stored as the columns of matrix A."""
    # Get the number of vectors.
    n = A.shape[1]
    for j in range(n):
        # To orthogonalize the vector in column j with respect to the
        # previous vectors, subtract from it its projection onto
        # each of the previous vectors.
        for k in range(j):
            A[:, j] -= np.dot(A[:, k], A[:, j]) * A[:, k]
        A[:, j] = A[:, j] / np.linalg.norm(A[:, j])
    return A

if __name__ == '__main__':
    A = np.array([[1.0, 1.0, 0.0], [1.0, 3.0, 1.0], [2.0, -1.0, 1.0]])
    print(gram_schmidt(A))

@mretegan your implementation would not work for a non-square matrix A. n = A.shape[0] gets the number of rows not columns

@seralouk

Why would you need that? The matrix is only used to store the vectors that have to be orthogonalized; the function assumes that they are the columns of matrix A. You can easily write the function to accept a list of vectors as input instead of a matrix.

Assume that I have A = np.array([[1.0, 1.0, 0.0,1.], [1.0, 3.0, 1.0,0.], [2.0, -1.0, 1.0,1.]]) and that I want to orthogonalize the columns of A., so that I have an orthogonalized version of A. How can I do this?

There is indeed a problem with the code. n should be number of columns and not rows: n = A.shape[1]. I've edited my previous answer.

There is indeed a problem with the code. n should be number of columns and not rows: n = A.shape[1].

Exactly. I just wanted to check that I am not missing something here. n = A.shape[1] should be fine

Maybe using list comprehension would make the code a bit less sluggsih for bigger matrices...so this version uses taht, You can even use it to check if a suspected "orthogonal" matric is actually othogonal by automatic_check=True. To check if you had two or more linearly dependent vectors used in the process, simply set orthogonality_check=True, and if the fucntion return False, then you had a linearly dependent vector in your set of vectors.

def Grahm_Schmidt(matrix, orthogonality_check=False, automatic_check=False, error_tol=1.e-10):
     """
     matrix is a matrix whose rows are  vectors to be turned into an ON-basis
     """
     veclist = list(matrix)
     newbasis = []
     def orth_check(Matrix):
          """
          This fucntion check for the pairwise orthogonality of the new basis
          """
          list_ = list(Matrix)
          dot_matrix = np.array([[m(item1, item2) for item1 in list_] for item2 in list_])
          if (dot_matrix - np.eye(dot_matrix.shape[0]) < error_tol).all():
               return True
          else:
               error = dot_matrix - np.eye(dot_matrix.shape[0])
               return False, np.max(error), np.min(error)

     def m(a,b):
               return np.matmul(a,b)
     def n(a):
          return np.linalg.norm(a)

     def motor(vector, ind):
          if ind == 0:
               new = vector/n(vector)
          else:
               L = np.array([newbasis[i]*m(newbasis[i],vector) for i in range(len(newbasis))])
               projections = np.sum(L, axis=0)
               NEW = vector - projections
               new = NEW/n(NEW)
          newbasis.append(new)
          
     [motor(vector, ind) for ind,vector in enumerate(veclist)]
     newbasis_matrix = np.array(newbasis)

     if orthogonality_check and automatic_check == False:
          return orth_check(newbasis_matrix)

     elif automatic_check:
          return orth_check(matrix)

     else:
          return newbasis_matrix

Pytorch version :

import torch    
def orthonormalize(vectors):    
    """    
        Orthonormalizes the vectors using gram schmidt procedure.    
    
        Parameters:    
            vectors: torch tensor, size (dimension, n_vectors)    
                    they must be linearly independant    
        Returns:    
            orthonormalized_vectors: torch tensor, size (dimension, n_vectors)    
    """    
    assert (vectors.size(1) <= vectors.size(0)), 'number of vectors must be smaller or equal to the dimension'    
    orthonormalized_vectors = torch.zeros_like(vectors)    
    orthonormalized_vectors[:, 0] = vectors[:, 0] / torch.norm(vectors[:, 0], p=2)    
    
    for i in range(1, orthonormalized_vectors.size(1)):    
        vector = vectors[:, i]    
        V = orthonormalized_vectors[:, :i]    
        PV_vector= torch.mv(V, torch.mv(V.t(), vector))    
        orthonormalized_vectors[:, i] = (vector - PV_vector) / torch.norm(vector - PV_vector, p=2)    
    
    return orthonormalized_vectors

Similar to @ingmarschuster 's vectorized version

def gram_schmidt(A,norm=True,row_vect=False):
    """Orthonormalizes vectors by gram-schmidt process
    
    Parameters
    -----------
    A : ndarray,
    Matrix having vectors in its columns
    
    norm : bool,
    Do you need Normalized vectors?
    
    row_vect: bool,
    Does Matrix A has vectors in its rows?
    
    Returns 
    -------
    G : ndarray,
    Matrix of orthogonal vectors 
    
    """
    if row_vect :
        # if true, transpose it to make column vector matrix
        A = A.T
    
    no_of_vectors = A.shape[1]
    G = A[:,0:1].copy() # copy the first vector in matrix
    # 0:1 is done to to be consistent with dimensions - [[1,2,3]]
    
    # iterate from 2nd vector to number of vectors
    for i in range(1,no_of_vectors):
        
        # calculates weights(coefficents) for every vector in G
        numerator = A[:,i].dot(G)
        denominator = np.diag(np.dot(G.T,G)) #to get elements in diagonal
        weights = np.squeeze(numerator/denominator)
        
        # projected vector onto subspace G 
        projected_vector = np.sum(weights * G,
                                  axis=1,
                                  keepdims=True)
        
        # orthogonal vector to subspace G
        orthogonalized_vector = A[:,i:i+1] - projected_vector
        
        # now add the orthogonal vector to our set 
        G = np.hstack((G,orthogonalized_vector))
        
    if norm :
        # to get orthoNORMAL vectors (unit orthogonal vectors)
        # replace zero to 1 to deal with division by 0 if matrix has 0 vector
        G = G/replace_zero(np.linalg.norm(G,axis=0))
    
    if row_vect:
        return G.T
    
    return G

def replace_zero(array): 
    
    for i in range(len(array)) :
        if array[i] == 0 : 
            array[i] = 1
    return array

G = np.array([[1,0,0],[1,1,0],[1,1,1],[1,1,1]])
gram_schmidt(G)
>
array([[ 0.5       , -0.8660254 ,  0.        ],
       [ 0.5       ,  0.28867513, -0.81649658],
       [ 0.5       ,  0.28867513,  0.40824829],
       [ 0.5       ,  0.28867513,  0.40824829]])

B = np.array([[1,0],[-2,0],[2,0]])
gram_schmidt(B)
>
array([[ 0.33333333,  0.        ],
       [-0.66666667,  0.        ],
       [ 0.66666667,  0.        ]])

Is there a 3D implementation of this procedure? I have n image in a non-orthogonal 3D system and I am trying to convert it to an orthogonal 3D image.

I simply implemented formula from a book (Maths for ML) and it is working fine with any matrix dimension (eg., 4x5).

import numpy as np
from numpy.linalg import inv

def _from_basis_to_project_matrix(B):
    """Formula page 87 book 'Maths For ML' """
    inv_Bt_B = inv(np.dot(B.T, B))
    proj_matrix = np.dot(B, np.dot(inv_Bt_B, B.T))
    return proj_matrix

def _projec(x, basis):
    """project `x` in the new vector space defined by `basis` """
    proj_matrix = _from_basis_to_project_matrix(basis)
    proj_x = np.dot(proj_matrix, x)
    return proj_x

def _get_col(x, col_id):
    """return column `col_id` from matrix `x` """
    raw_col_val = x[:, col_id]
    col_as_row = np.array([raw_col_val]).T
    return col_as_row

def gram_schmidt(B):
    nb_col = B.shape[1]
    first_col=_get_col(B, 0)
    
    U_vecs = [first_col]
    
    for i in range(1, nb_col):
        B_i = _get_col(B, i)
        U_i_1 = np.concatenate(U_vecs, axis=1)
        p = _projec(B_i, U_i_1)
        U_i = B_i - p
        U_vecs.append(U_i)

    return U_vecs

if __name__ == '__main__':
    # B=np.array([[2,1],[0,1]])
    # B=np.array([[1,0],[1,1],[1,2]])
    B = np.array([[0, 1, -3, -1],
                  [-1, -3, 4, -3],
                  [2, 1, 1, 5],
                  [0, -1, 2, 0],
                  [2, 2, 1, 7]
                  ])

    U = gram_schmidt(B)

    print(np.round(U, 2))

Is there a 3D implementation of this procedure? I have n image in a non-orthogonal 3D system and I am trying to convert it to an orthogonal 3D image.

If you are working solely with vectors, then there is no problem with any number of dimensions (assuming you have sufficiently many linearly independent vectors to form an orthonormal basis).

This should be fine even if you are using an MxNxP (matrix/tensor, whatever). In that case you'd be able to use vectors of length MP or NP. Just have some mapping from the interval [1,MP] to a "2D" slice of the "3D" object.
For instance, take phi: x -> (x // M, x % M) [or something like it] to map x, the index into the vector, to i.j the indices in the slice.
phi would then just fill rows or columns (depending on whether you're using a row-major convention or not).

#  this maps a vector of length 25 to a 5x5 matrix
for x in range(1,25):
    print(f'({x // 5},{x % 5})')

Linearity shouldn't be compromised by the method you use to traverse the "lower-dimensional" slice.