Матричная инверсия без Numpy
Я хочу инвертировать матрицу без использования numpy.linalg.inv .
Причина в том, что я использую Numba для ускорения кода, но numpy.linalg.inv не поддерживается, поэтому мне интересно, могу ли я инвертировать матрицу с «классическим» кодом Python.
- Как совместить два массива
- изменение имени импорта в python
- Ошибка при использовании `struct.pack` для записи данных в файл
- Как пропустить поддеревья в обход AST, используя привязки python для libclang
- вернуть значение из одного скрипта python в другое
- Не-ASCII-символ '\ x90', выполняющий pserve на окнах внутри virtualenv
С numpy.linalg.inv пример кода будет выглядеть так:
import numpy as np M = np.array([[1,0,0],[0,1,0],[0,0,1]]) Minv = np.linalg.inv(M) 3 Solutions collect form web for “Матричная инверсия без Numpy”
Вот более элегантное и масштабируемое решение, imo. Он будет работать для любой матрицы nxn, и вы можете использовать его для других методов. Обратите внимание, что getMatrixInverse (m) принимает в качестве массива массив массивов. Пожалуйста, не стесняйтесь задавать любые вопросы.
def transposeMatrix(m): t = [] for r in range(len(m)): tRow = [] for c in range(len(m[r])): if c == r: tRow.append(m[r][c]) else: tRow.append(m[c][r]) t.append(tRow) return t def getMatrixMinor(m,i,j): return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])] def getMatrixDeternminant(m): #base case for 2x2 matrix if len(m) == 2: return m[0][0]*m[1][1]-m[0][1]*m[1][0] determinant = 0 for c in range(len(m)): determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c)) return determinant def getMatrixInverse(m): determinant = getMatrixDeternminant(m) #special case for 2x2 matrix: if len(m) == 2: return [[m[1][1]/determinant, -1*m[0][1]/determinant], [-1*m[1][0]/determinant, m[0][0]/determinant]] #find matrix of cofactors cofactors = [] for r in range(len(m)): cofactorRow = [] for c in range(len(m)): minor = getMatrixMinor(m,r,c) cofactorRow.append(((-1)**(r+c)) * getMatrixDeternminant(minor)) cofactors.append(cofactorRow) cofactors = transposeMatrix(cofactors) for r in range(len(cofactors)): for c in range(len(cofactors)): cofactors[r][c] = cofactors[r][c]/determinant return cofactors
Для матрицы размером 4 x 4, вероятно, просто ОК, чтобы использовать математическую формулу, которую вы можете найти с помощью формулы «Google for for 4 to 4». Например, здесь (я не могу ручаться за его точность):
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
В общем, инвертирование общей матрицы не для слабонервных. Вы должны знать все математически сложные случаи и знать, почему они не будут применяться к вашему использованию, и поймать их, когда вы получаете математически патологические входы (которые или возвращают результаты с низкой точностью или числовым мусором в знании, что это не имеет значения в вашем случае использования, если вы на самом деле не делите на нуль или переполняете MAXFLOAT … который вы можете поймать с обработчиком исключений и представить как «Ошибка: матрица является единственной или очень близкой к ней»).
Как правило, лучше, как программист использовать библиотечный код, написанный специалистами по математической математике, если вы не хотите тратить время на понимание физической и математической природы конкретной проблемы, с которой вы обращаетесь, и стать вашим собственным экспертом по математике в своей собственной специализации.
Я использовал формулу из http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html, чтобы написать функцию, которая выполняет инверсию матрицы 4×4:
import numpy as np def myInverse(A): detA = np.linalg.det(A) b00 = A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] - A[1,3]*A[2,2]*A[3,1] b01 = A[0,1]*A[2,3]*A[3,2] + A[0,2]*A[2,1]*A[3,3] + A[0,3]*A[2,2]*A[3,1] - A[0,1]*A[2,2]*A[3,3] - A[0,2]*A[2,3]*A[3,1] - A[0,3]*A[2,1]*A[3,2] b02 = A[0,1]*A[1,2]*A[3,3] + A[0,2]*A[1,3]*A[3,1] + A[0,3]*A[1,1]*A[3,2] - A[0,1]*A[1,3]*A[3,2] - A[0,2]*A[1,1]*A[3,3] - A[0,3]*A[1,2]*A[3,1] b03 = A[0,1]*A[1,3]*A[2,2] + A[0,2]*A[1,1]*A[2,3] + A[0,3]*A[1,2]*A[2,1] - A[0,1]*A[1,2]*A[2,3] - A[0,2]*A[1,3]*A[2,1] - A[0,3]*A[1,1]*A[2,2] b10 = A[1,0]*A[2,3]*A[3,2] + A[1,2]*A[2,0]*A[3,3] + A[1,3]*A[2,2]*A[3,0] - A[1,0]*A[2,2]*A[3,3] - A[1,2]*A[2,3]*A[3,0] - A[1,3]*A[2,0]*A[3,2] b11 = A[0,0]*A[2,2]*A[3,3] + A[0,2]*A[2,3]*A[3,0] + A[0,3]*A[2,0]*A[3,2] - A[0,0]*A[2,3]*A[3,2] - A[0,2]*A[2,0]*A[3,3] - A[0,3]*A[2,2]*A[3,0] b12 = A[0,0]*A[1,3]*A[3,2] + A[0,2]*A[1,0]*A[3,3] + A[0,3]*A[1,2]*A[3,0] - A[0,0]*A[1,2]*A[3,3] - A[0,2]*A[1,3]*A[3,0] - A[0,3]*A[1,0]*A[3,2] b13 = A[0,0]*A[1,2]*A[2,3] + A[0,2]*A[1,3]*A[2,0] + A[0,3]*A[1,0]*A[2,2] - A[0,0]*A[1,3]*A[2,2] - A[0,2]*A[1,0]*A[2,3] - A[0,3]*A[1,2]*A[2,0] b20 = A[1,0]*A[2,1]*A[3,3] + A[1,1]*A[2,3]*A[3,0] + A[1,3]*A[2,0]*A[3,1] - A[1,0]*A[2,3]*A[3,1] - A[1,1]*A[2,0]*A[3,3] - A[1,3]*A[2,1]*A[3,0] b21 = A[0,0]*A[2,3]*A[3,1] + A[0,1]*A[2,0]*A[3,3] + A[0,3]*A[2,1]*A[3,0] - A[0,0]*A[2,1]*A[3,3] - A[0,1]*A[2,3]*A[3,0] - A[0,3]*A[2,0]*A[3,1] b22 = A[0,0]*A[1,1]*A[3,3] + A[0,1]*A[1,3]*A[3,0] + A[0,3]*A[1,0]*A[3,1] - A[0,0]*A[1,3]*A[3,1] - A[0,1]*A[1,0]*A[3,3] - A[0,3]*A[1,1]*A[3,0] b23 = A[0,0]*A[1,3]*A[2,1] + A[0,1]*A[1,0]*A[2,3] + A[0,3]*A[1,1]*A[2,0] - A[0,0]*A[1,1]*A[2,3] - A[0,1]*A[1,3]*A[2,0] - A[0,3]*A[1,0]*A[2,1] b30 = A[1,0]*A[2,2]*A[3,1] + A[1,1]*A[2,0]*A[3,2] + A[1,2]*A[2,1]*A[3,0] - A[1,0]*A[2,1]*A[3,2] - A[1,1]*A[2,2]*A[3,0] - A[1,2]*A[2,0]*A[3,1] b31 = A[0,0]*A[2,1]*A[3,2] + A[0,1]*A[2,2]*A[3,0] + A[0,2]*A[2,0]*A[3,1] - A[0,0]*A[2,2]*A[3,1] - A[0,1]*A[2,0]*A[3,2] - A[0,2]*A[2,1]*A[3,0] b32 = A[0,0]*A[1,2]*A[3,1] + A[0,1]*A[1,0]*A[3,2] + A[0,2]*A[1,1]*A[3,0] - A[0,0]*A[1,1]*A[3,2] - A[0,1]*A[1,2]*A[3,0] - A[0,2]*A[1,0]*A[3,1] b33 = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,0] + A[0,2]*A[1,0]*A[2,1] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2] - A[0,2]*A[1,1]*A[2,0] Ainv = np.array([[b00, b01, b02, b03], [b10, b11, b12, b13], [b20, b21, b22, b23], [b30, b31, b32, b33]]) / detA return Ainv