Compute the inverse of a matrix
Posted on Tue 14 June 2022 in math
Introduction
Recently, I wrote a C++ Library [1] for linear algebra, and wanted to wrote a method to inverse a matrix. I got stuck for a while, but finally, I found a solution. Here it is.
The algorithm
The algorithm is rather simple. It firstly check if the matrix is squared, as is it required for the inverse.
Then it computes the determinant of the matrix, and check if it is not zero.
If not, it computes for each element of the matrix, the submatrix of the transposed matrix without the row and column of the element.
Finally it set the element of the matrix as the formula goes:
$$
A^{-1}_{i, j} = (-1)^{i+j} \times \det(subA.T) / \det(A)
$$
Python implementation
For the sake of simplicity, here we describe the algorithm in Python, rather than C++.
#!/usr/bin/env python3
"""Get the submatrix of the matrix without the row and column of the element at position (a, b)."""
def submatrix(a, b, matrix):
submatrix = [[0 for i in range(len(matrix[j]) - 1)] for j in range(len(matrix) - 1)]
for j in range(len(submatrix)):
y = j if j < b else j + 1
for i in range(len(submatrix[j])):
x = i if i < a else i + 1
submatrix[j][i] = matrix[y][x]
return submatrix
"""Return (-1)^i"""
def toggle(i):
return 1 if i % 2 == 0 else -1
"""Compute the determinant of the matrix"""
def determinant(matrix):
if len(matrix) != len(matrix[0]):
raise ValueError("The matrix must be squared.")
if len(matrix) == 1:
return matrix[0][0]
elif len(matrix) == 2:
return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
else:
det = 0
for i in range(len(matrix)):
det += toggle(i) * matrix[0][i] * determinant(submatrix(i, 0, matrix))
return det
"""Get the transposed matrix"""
def transpose(matrix):
transposed = [[0 for i in range(len(matrix[j]))] for j in range(len(matrix))]
for i in range(len(matrix)):
for j in range(len(matrix[i])):
transposed[j][i] = matrix[i][j]
return transposed
"""Invert the matrix if possible"""
def invert(matrix):
if len(matrix) != len(matrix[0]):
raise ValueError("The matrix must be squared.")
det = determinant(matrix)
if det == 0:
raise ValueError("The matrix is not invertible.")
else :
inverse = [[0 for i in range(len(matrix[j]))] for j in range(len(matrix))]
transposed = transpose(matrix)
for j in range(len(matrix)):
for i in range(len(matrix[j])):
sm = submatrix(i, j, transposed)
sm_det = determinant(sm)
inverse[j][i] = toggle(i + j) * sm_det / det
return inverse
def main():
# Test submatrix calculation
matrix = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
minor = [[5, 6], [8, 9]]
assert submatrix(0, 0, matrix) == minor, "Minor does not work properly."
# Test submatrix calculation 2
matrix = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
minor = [[1, 3], [7, 9]]
assert submatrix(1, 1, matrix) == minor, "Minor does not work properly."
# Test tranposition of matrix
matrix = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
transposed = [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
assert transposed == transpose(matrix), "Transpose does not work properly."
# Test determinant calculation
matrix = [[8, 2, 3], [4, 7, 6], [7, 8, 9]]
assert determinant(matrix) == 81, "Determinant does not work properly."
# Test inversion of a matrix
matrix = [[1, 4],[3, 2]]
assert invert(matrix) == [[-1/5, 2/5], [3/10, -1/10]], "Inversion does not work properly."
print("All is fine !")
if __name__ == "__main__":
main()
This little piece of code is available here.