Skip to content

Instantly share code, notes, and snippets.

@ad-1
Last active April 13, 2023 20:56
Show Gist options
  • Select an option

  • Save ad-1/2e600d0f17f1dd3cac09e4f5d3870abc to your computer and use it in GitHub Desktop.

Select an option

Save ad-1/2e600d0f17f1dd3cac09e4f5d3870abc to your computer and use it in GitHub Desktop.
Find the Inverse of a Matrix using Python
# form the augmented matrix by concatenating A and I
M = np.concatenate((M, I), axis=1)
# identity matrix with same shape as A
I = np.identity(n=n)
# matrix to find inverse
A = np.array([[1, 0, 1], [0, 1, 1], [1, 1, 1]])
import time
# ===========
# start timer
start = time.time()
# run custom method
result_1 = invert_matrix(A)
# record time
time_1 = time.time() - start
# ===========
# start timer
start = time.time()
# run numpy method
result_2 = np.linalg.inv(A)
# record time
time_2 = time.time() - start
# assert numpy arrays are equal
np.testing.assert_array_almost_equal(result_1, result_2, decimal=6)
# compare times
print(f'numpy is {time_1 - time_2} seconds faster than the manual implementation')
import time
import numpy as np
# write rows in reduced row echelon (rref) form
def invert_matrix(M):
# store dimension
n = M.shape[0]
# A must be square with non-zero determinant
# assert np.linalg.det(M) != 0
# identity matrix with same shape as A
I = np.identity(n=n)
# form the augmented matrix by concatenating A and I
M = np.concatenate((M, I), axis=1)
# move all zeros to buttom of matrix
M = np.concatenate((M[np.any(M != 0, axis=1)], M[np.all(M == 0, axis=1)]), axis=0)
# iterate over matrix rows
for i in range(0, n):
# initialize row-swap iterator
j = 1
# select pivot value
pivot = M[i][i]
# find next non-zero leading coefficient
while pivot == 0 and i + j < n:
# perform row swap operation
M[[i, i + j]] = M[[i + j, i]]
# incrememnt row-swap iterator
j += 1
# get new pivot
pivot = M[i][i]
# if pivot is zero, remaining rows are all zeros
if pivot == 0:
# return inverse matrix
return M[:, n:]
# extract row
row = M[i]
# get 1 along the diagonal
M[i] = row / pivot
# iterate over all rows except pivot to get augmented matrix into reduced row echelon form
for j in [k for k in range(0, n) if k != i]:
# subtract current row from remaining rows
M[j] = M[j] - M[i] * M[j][i]
# return inverse matrix
return M[:, n:]
if __name__ == '__main__':
# matrix to find inverse
A = np.array([[1, 0, 1], [0, 1, 1], [1, 1, 1]])
# A = np.random.rand(1000, 1000)
start = time.time()
result_1 = invert_matrix(A)
time_1 = time.time() - start
start = time.time()
result_2 = np.linalg.inv(A)
time_2 = time.time() - start
print(f'numpy is {time_1 - time_2} seconds faster than the manual implementation')
np.testing.assert_array_almost_equal(result_1, result_2, decimal=6)
# define large number of rows
size = 500
# create a square array of random floating point values
A = np.random.rand(size, size)
# define large number of rows
size = 500
# create a square array of random floating point values
A = np.random.rand(size, size)
# iterate over matrix rows
for i in range(0, n):
# select pivot value
pivot = M[i][i]
# extract row
row = M[i]
# get 1 along the diagonal
M[i] = row / pivot
# iterate over all rows except pivot to get augmented matrix into reduced row echelon form
for j in [k for k in range(0, n) if k != i]:
# subtract current row from remaining rows
M[j] = M[j] - M[i] * M[j][i]
# extract inverse matrix
inverse = M[:, n:]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment