Skip to content

Instantly share code, notes, and snippets.

@iizukak
Created October 14, 2011 18:18
Show Gist options
  • Save iizukak/1287876 to your computer and use it in GitHub Desktop.
Save iizukak/1287876 to your computer and use it in GitHub Desktop.
Gram-Schmidt Orthogonization using Numpy
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))
@PierrickPochelu
Copy link

PierrickPochelu commented Nov 16, 2022

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))

@egoughnour
Copy link

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment