Add matrix inversion (#582)

* adds matrix inversion for n x n matrices

* adds link in readme and edge cases in file

* fixes typo

Author: abertsch72

README.md

@@ -228,6 +228,7 @@ If you want to uninstall algorithms, it is as simple as:
     - [copy_transform](algorithms/matrix/copy_transform.py)
     - [count_paths](algorithms/matrix/count_paths.py)
     - [matrix_rotation.txt](algorithms/matrix/matrix_rotation.txt)
+    - [matrix_inversion](algorithms/matrix/matrix_inversion.py)
     - [matrix_multiplication](algorithms/matrix/multiply.py)
     - [rotate_image](algorithms/matrix/rotate_image.py)
     - [search_in_sorted_matrix](algorithms/matrix/search_in_sorted_matrix.py)

algorithms/matrix/matrix_inversion.py

@@ -0,0 +1,124 @@
+"""
+Inverts an invertible n x n matrix -- i.e., given an n x n matrix A, returns
+an n x n matrix B such that AB = BA = In, the n x n identity matrix.
+
+For a 2 x 2 matrix, inversion is simple using the cofactor equation. For
+larger matrices, this is a four step process:
+1. calculate the matrix of minors: create an n x n matrix by considering each
+position in the original matrix in turn. Exclude the current row and column
+and calculate the determinant of the remaining matrix, then place that value
+in the current position's equivalent in the matrix of minors.
+2. create the matrix of cofactors: take the matrix of minors and multiply
+alternate values by -1 in a checkerboard pattern.
+3. adjugate: hold the top left to bottom right diagonal constant, but swap all
+other values over it.
+4. multiply the adjugated matrix by 1 / the determinant of the original matrix
+
+This code combines steps 1 and 2 into one method to reduce traversals of the
+matrix.
+
+Possible edge cases: will not work for 0x0 or 1x1 matrix, though these are
+trivial to calculate without use of this file.
+"""
+import fractions
+
+
+def invert_matrix(m):
+    """invert an n x n matrix"""
+    # Error conditions
+    if not array_is_matrix(m):
+        print("Invalid matrix: array is not a matrix")
+        return [[-1]];
+    elif len(m) != len(m[0]):
+        print("Invalid matrix: matrix is not square")
+        return [[-2]];
+    elif len(m) < 2:
+        print("Invalid matrix: matrix is too small")
+        return [[-3]];
+    elif get_determinant(m) == 0:
+        print("Invalid matrix: matrix is square, but singular (determinant = 0)")
+        return [[-4]];
+
+    # Calculation
+    elif len(m) == 2:
+        # simple case
+        multiplier = 1 / get_determinant(m)
+        inverted = [[multiplier] * len(m) for n in range(len(m))]
+        inverted[0][1] = inverted[0][1] * -1 * m[0][1]
+        inverted[1][0] = inverted[1][0] * -1 * m[1][0]
+        inverted[0][0] = multiplier * m[1][1]
+        inverted[1][1] = multiplier * m[0][0]
+        return inverted
+    else:
+        """some steps combined in helpers to reduce traversals"""
+        # get matrix of minors w/ "checkerboard" signs
+        m_of_minors = get_matrix_of_minors(m)
+
+        # calculate determinant (we need to know 1/det)
+        multiplier = fractions.Fraction(1, get_determinant(m))
+
+        # adjugate (swap on diagonals) and multiply by 1/det
+        inverted = transpose_and_multiply(m_of_minors, multiplier)
+
+        return inverted
+
+
+def get_determinant(m):
+    """recursively calculate the determinant of an n x n matrix, n >= 2"""
+    if len(m) == 2:
+        # trivial case
+        return (m[0][0] * m[1][1]) - (m[0][1] * m[1][0])
+    else:
+        sign = 1
+        det = 0
+        for i in range(len(m)):
+            det += sign * m[0][i] * get_determinant(get_minor(m, 0, i))
+            sign *= -1
+        return det
+
+
+def get_matrix_of_minors(m):
+    """get the matrix of minors and alternate signs"""
+    matrix_of_minors = [[0 for i in range(len(m))] for j in range(len(m))]
+    for row in range(len(m)):
+        for col in range(len(m[0])):
+            if (row + col) % 2 == 0:
+                sign = 1
+            else:
+                sign = -1
+            matrix_of_minors[row][col] = sign * get_determinant(get_minor(m, row, col))
+    return matrix_of_minors
+
+
+def get_minor(m, row, col):
+    """
+    get the minor of the matrix position m[row][col]
+    (all values m[r][c] where r != row and c != col)
+    """
+    minors = []
+    for i in range(len(m)):
+        if i != row:
+            new_row = m[i][:col]
+            new_row.extend(m[i][col + 1:])
+            minors.append(new_row)
+    return minors
+
+
+def transpose_and_multiply(m, multiplier=1):
+    """swap values along diagonal, optionally adding multiplier"""
+    for row in range(len(m)):
+        for col in range(row + 1):
+            temp = m[row][col] * multiplier
+            m[row][col] = m[col][row] * multiplier
+            m[col][row] = temp
+    return m
+
+
+def array_is_matrix(m):
+    if len(m) == 0:
+        return False
+    first_col = len(m[0])
+    for row in m:
+        if len(row) != first_col:
+            return False
+    return True

tests/test_matrix.py

@@ -3,6 +3,7 @@
     copy_transform,
     crout_matrix_decomposition,
     cholesky_matrix_decomposition,
+    matrix_inversion,
     multiply,
     rotate_image,
     sparse_dot_vector,
@@ -127,6 +128,37 @@ def test_cholesky_matrix_decomposition(self):
                              [[5, 1.2, 0.3, -0.6], [1.2, 6, -0.4, 0.9],
                               [0.3, -0.4, 8, 1.7], [-0.6, 0.9, 1.7, 10]]))
 
+class TestInversion(unittest.TestCase):
+    """[summary]
+    Test for the file matrix_inversion.py
+
+    Arguments:
+        unittest {[type]} -- [description]
+    """
+    def test_inversion(self):
+        from fractions import Fraction
+
+        m1 = [[1, 1], [1, 2]]
+        self.assertEqual(matrix_inversion.invert_matrix(m1), [[2, -1], [-1, 1]])
+
+        m2 = [[1, 2], [3, 4, 5]]
+        self.assertEqual(matrix_inversion.invert_matrix(m2), [[-1]])
+
+        m3 = [[1, 1, 1, 1], [2, 2, 2, 2]]
+        self.assertEqual(matrix_inversion.invert_matrix(m3), [[-2]])
+
+        m4 = [[1]]
+        self.assertEqual(matrix_inversion.invert_matrix(m4), [[-3]])
+
+        m5 = [[1, 2, 3]   , [4, 5, 6], [7, 8, 9]]
+        self.assertEqual(matrix_inversion.invert_matrix(m5), [[-4]])
+
+        m6 = [[3, 5, 1], [2, 5, 0], [1, 9, 8]]
+        self.assertEqual(matrix_inversion.invert_matrix(m6), [[Fraction(40, 53), Fraction(-31, 53), Fraction(-5, 53)],
+                                                              [Fraction(-16, 53), Fraction(23, 53), Fraction(2, 53)],
+                                                              [Fraction(13, 53), Fraction(-22, 53), Fraction(5, 53)]])
+
+
 
 class TestMultiply(unittest.TestCase):
     """[summary]