|
| 1 | +""" |
| 2 | +Inverts an invertible n x n matrix -- i.e., given an n x n matrix A, returns |
| 3 | +an n x n matrix B such that AB = BA = In, the n x n identity matrix. |
| 4 | +
|
| 5 | +For a 2 x 2 matrix, inversion is simple using the cofactor equation. For |
| 6 | +larger matrices, this is a four step process: |
| 7 | +1. calculate the matrix of minors: create an n x n matrix by considering each |
| 8 | +position in the original matrix in turn. Exclude the current row and column |
| 9 | +and calculate the determinant of the remaining matrix, then place that value |
| 10 | +in the current position's equivalent in the matrix of minors. |
| 11 | +2. create the matrix of cofactors: take the matrix of minors and multiply |
| 12 | +alternate values by -1 in a checkerboard pattern. |
| 13 | +3. adjugate: hold the top left to bottom right diagonal constant, but swap all |
| 14 | +other values over it. |
| 15 | +4. multiply the adjugated matrix by 1 / the determinant of the original matrix |
| 16 | +
|
| 17 | +This code combines steps 1 and 2 into one method to reduce traversals of the |
| 18 | +matrix. |
| 19 | +
|
| 20 | +Possible edge cases: will not work for 0x0 or 1x1 matrix, though these are |
| 21 | +trivial to calculate without use of this file. |
| 22 | +""" |
| 23 | +importfractions |
| 24 | + |
| 25 | + |
| 26 | +definvert_matrix(m): |
| 27 | +"""invert an n x n matrix""" |
| 28 | +# Error conditions |
| 29 | +ifnotarray_is_matrix(m): |
| 30 | +print("Invalid matrix: array is not a matrix") |
| 31 | +return [[-1]]; |
| 32 | +eliflen(m)!=len(m[0]): |
| 33 | +print("Invalid matrix: matrix is not square") |
| 34 | +return [[-2]]; |
| 35 | +eliflen(m)<2: |
| 36 | +print("Invalid matrix: matrix is too small") |
| 37 | +return [[-3]]; |
| 38 | +elifget_determinant(m)==0: |
| 39 | +print("Invalid matrix: matrix is square, but singular (determinant = 0)") |
| 40 | +return [[-4]]; |
| 41 | + |
| 42 | +# Calculation |
| 43 | +eliflen(m)==2: |
| 44 | +# simple case |
| 45 | +multiplier=1/get_determinant(m) |
| 46 | +inverted= [[multiplier]*len(m)forninrange(len(m))] |
| 47 | +inverted[0][1]=inverted[0][1]*-1*m[0][1] |
| 48 | +inverted[1][0]=inverted[1][0]*-1*m[1][0] |
| 49 | +inverted[0][0]=multiplier*m[1][1] |
| 50 | +inverted[1][1]=multiplier*m[0][0] |
| 51 | +returninverted |
| 52 | +else: |
| 53 | +"""some steps combined in helpers to reduce traversals""" |
| 54 | +# get matrix of minors w/ "checkerboard" signs |
| 55 | +m_of_minors=get_matrix_of_minors(m) |
| 56 | + |
| 57 | +# calculate determinant (we need to know 1/det) |
| 58 | +multiplier=fractions.Fraction(1,get_determinant(m)) |
| 59 | + |
| 60 | +# adjugate (swap on diagonals) and multiply by 1/det |
| 61 | +inverted=transpose_and_multiply(m_of_minors,multiplier) |
| 62 | + |
| 63 | +returninverted |
| 64 | + |
| 65 | + |
| 66 | +defget_determinant(m): |
| 67 | +"""recursively calculate the determinant of an n x n matrix, n >= 2""" |
| 68 | +iflen(m)==2: |
| 69 | +# trivial case |
| 70 | +return (m[0][0]*m[1][1])- (m[0][1]*m[1][0]) |
| 71 | +else: |
| 72 | +sign=1 |
| 73 | +det=0 |
| 74 | +foriinrange(len(m)): |
| 75 | +det+=sign*m[0][i]*get_determinant(get_minor(m,0,i)) |
| 76 | +sign*=-1 |
| 77 | +returndet |
| 78 | + |
| 79 | + |
| 80 | +defget_matrix_of_minors(m): |
| 81 | +"""get the matrix of minors and alternate signs""" |
| 82 | +matrix_of_minors= [[0foriinrange(len(m))]forjinrange(len(m))] |
| 83 | +forrowinrange(len(m)): |
| 84 | +forcolinrange(len(m[0])): |
| 85 | +if (row+col)%2==0: |
| 86 | +sign=1 |
| 87 | +else: |
| 88 | +sign=-1 |
| 89 | +matrix_of_minors[row][col]=sign*get_determinant(get_minor(m,row,col)) |
| 90 | +returnmatrix_of_minors |
| 91 | + |
| 92 | + |
| 93 | +defget_minor(m,row,col): |
| 94 | +""" |
| 95 | + get the minor of the matrix position m[row][col] |
| 96 | + (all values m[r][c] where r != row and c != col) |
| 97 | + """ |
| 98 | +minors= [] |
| 99 | +foriinrange(len(m)): |
| 100 | +ifi!=row: |
| 101 | +new_row=m[i][:col] |
| 102 | +new_row.extend(m[i][col+1:]) |
| 103 | +minors.append(new_row) |
| 104 | +returnminors |
| 105 | + |
| 106 | + |
| 107 | +deftranspose_and_multiply(m,multiplier=1): |
| 108 | +"""swap values along diagonal, optionally adding multiplier""" |
| 109 | +forrowinrange(len(m)): |
| 110 | +forcolinrange(row+1): |
| 111 | +temp=m[row][col]*multiplier |
| 112 | +m[row][col]=m[col][row]*multiplier |
| 113 | +m[col][row]=temp |
| 114 | +returnm |
| 115 | + |
| 116 | + |
| 117 | +defarray_is_matrix(m): |
| 118 | +iflen(m)==0: |
| 119 | +returnFalse |
| 120 | +first_col=len(m[0]) |
| 121 | +forrowinm: |
| 122 | +iflen(row)!=first_col: |
| 123 | +returnFalse |
| 124 | +returnTrue |