[PYTHON] Numpy없이 행렬 반전

Numpy없이 행렬 반전

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)

해결법

1. ==============================

   1.여기에 더 우아하고 확장 가능한 솔루션 인 imo가 있습니다. 어떤 nxn 행렬에서도 작동 할 것이고 다른 방법에 대한 사용법을 찾을 수 있습니다. getMatrixInverse (m)는 배열의 배열을 입력으로받습니다. 언제든지 질문하십시오.

   여기에 더 우아하고 확장 가능한 솔루션 인 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
   

2. ==============================

   2.4 x 4 행렬의 경우 수학 공식을 사용하는 것이 좋습니다. Googling "4 x 4 행렬 역 분식의 수식"을 사용하여 찾을 수 있습니다. 예를 들면 다음과 같습니다 (정확성을 보장 할 수는 없습니다).

   4 x 4 행렬의 경우 수학 공식을 사용하는 것이 좋습니다. Googling "4 x 4 행렬 역 분식의 수식"을 사용하여 찾을 수 있습니다. 예를 들면 다음과 같습니다 (정확성을 보장 할 수는 없습니다).

   http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html

   일반적으로 일반적인 행렬을 뒤집는 것은 희미한 것이 아닙니다. 수학적으로 어려운 모든 경우에 대해 알고 있어야하며 사용법에 적용되지 않는 이유를 알아야하며 수학적으로 병적 인 입력을 받거나 (즉, 정확도가 낮거나 숫자가 많은 결과를 당신이 실제로 0으로 나누거나 MAXFLOAT를 넘치게하지 않는다면, 예외 처리기로 잡아서 "오류 : 행렬은 단수 또는 매우 가깝습니다"라고 제시하면, 사용 사례에서는 중요하지 않습니다.

   일반적으로 수학 전문가가 작성한 라이브러리 코드를 사용하는 프로그래머로서 더 낫습니다. 특정 분야의 문제를 수학적으로 이해하고 자신의 전문 분야에서 수학 전문가가되기까지는 시간을 투자하지 않는 것이 좋습니다.

3. ==============================

   3.적어도 7 월 16 일, 2018 Numba는 빠른 매트릭스 역행렬을 가지고 있습니다. (표준 NumPy 역수 및 기타 연산을 여기에서 오버로드하는 방법을 볼 수 있습니다.)

   적어도 7 월 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)
   

   스피드 업은 SciPy가 아니라 NumPy 역행렬에서만 작동합니다 (예상대로).

   약간 큰 매트릭스 :

   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는 따라 잡고 있습니다.

4. ==============================

   4.http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html의 공식을 사용하여 4x4 행렬의 역변환을 수행하는 함수를 작성했습니다.

   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
   

from https://stackoverflow.com/questions/32114054/matrix-inversion-without-numpy by cc-by-sa and MIT license