Я хочу инвертировать матрицу без использования numpy.linalg.inv.
Причина в том, что я использую Numba для ускорения кода, но numpy.linalg.inv не поддерживается, поэтому мне интересно, могу ли я инвертировать матрицу с "классическим" кодом Python.
С numpy.linalg.inv пример кода будет выглядеть следующим образом:
import numpy as np
M = np.array([[1,0,0],[0,1,0],[0,0,1]])
Minv = np.linalg.inv(M)
Я использовал формулу http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html для записи функции, которая выполняет инверсию матрицы 4x4:
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
Вот более элегантное и масштабируемое решение, imo. Он будет работать для любой матрицы nxn, и вы можете использовать его для других методов. Обратите внимание, что getMatrixInverse (m) принимает в качестве массива массив массивов. Пожалуйста, не стесняйтесь задавать любые вопросы.
def transposeMatrix(m):
return map(list,zip(*m))
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 х 4, вероятно, всего лишь ОК, чтобы использовать математическую формулу, которую вы можете найти с помощью формулы "Google" для формулы 4 на 4 матрицы. Например, здесь (я не могу ручаться за его точность):
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
В общем, инвертирование общей матрицы не для слабонервных. Вы должны знать все математически сложные случаи и знать, почему они не будут применяться к вашему использованию, и поймать их, когда вы получаете математически патологические входы (которые или возвращают результаты с низкой точностью или числовым мусором в знании, что это не имеет значения в вашем случае использования, если вы на самом деле не делите на нуль или переполняете MAXFLOAT... который вы можете поймать с обработчиком исключений и представить как "Ошибка: матрица является единственной или очень близкой к ней" ).
Как правило, лучше программист использовать библиотечный код, написанный специалистами по математической математике, если вы не хотите потратить время на понимание физической и математической природы конкретной проблемы, к которой вы обращаетесь, и стать вашим собственным экспертом по математике в своем собственном специалисте поле.
По крайней мере, 16 июля 2018 года Numba имеет быструю обратную матрицу. (Вы можете увидеть, как они перегружают стандартные обратные NumPy и другие операции здесь.)
Вот результаты моего бенчмаркинга:
import numpy as np
from scipy import linalg as sla
from scipy import linalg as nla
import numba
def gen_ex(d0):
x = np.random.randn(d0,d0)
return x.T + x
@numba.jit
def inv_nla_jit(A):
return np.linalg.inv(A)
@numba.jit
def inv_sla_jit(A):
return sla.inv(A)
Для маленьких матриц это особенно быстро:
ex1 = gen_ex(4)
%timeit inv_nla_jit(ex1) # NumPy + Numba
%timeit inv_sla_jit(ex1) # SciPy + Numba
%timeit nla.inv(ex1) # NumPy
%timeit sla.inv(ex1) # SciPy
[Из]
2.54 µs ± 467 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
67.3 µs ± 9.18 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
63.5 µs ± 7.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
56.6 µs ± 5.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Обратите внимание, что ускорение работает только для инверсии NumPy, а не SciPy (как и ожидалось).
Чуть крупнее матрица:
ex2 = gen_ex(40)
%timeit inv_nla_jit(ex2) # NumPy + Numba
%timeit inv_sla_jit(ex2) # SciPy + Numba
%timeit nla.inv(ex2) # NumPy
%timeit sla.inv(ex2) # SciPy
[Из]
131 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
278 µs ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
231 µs ± 24.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
189 µs ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Так что здесь все еще ускорение, но SciPy наверстывает упущенное.