파이썬으로 행렬 구현하기 (3) - 가우스 조던 소거법을 이용한 역행렬

[알고리즘 & 코딩테스트/코딩테스트] - 파이썬으로 행렬 구현하기 (1)

[알고리즘 & 코딩테스트/코딩테스트] - 파이썬으로 행렬 구현하기 (2)

---

지난 포스트에서의 행렬 덧셈, 뺄셈, 곱셈에 이어 이번 포스트에서는 행렬의 역변환, 역행렬을 구현해보도록 하겠습니다. 행렬의 행과 열의 크기가 같은 정방행렬의 역행렬을 구하는 방법으로는 대표적으로 가우스-조던 소거법이 있습니다. 

$AX=I$에서 A의 역행렬 $X$를 구하는 가우스-조던 소거법의 기본은 $AX$를 $IX$로 바꿔주기 위해 1) 행렬의 한 행을 상수배, 2) 행렬의 두 행을 바꿈, 3) 한 행을 상수배하여 다른 행에 더하는 기본 행연산을 수행한다는 것입니다. 행렬의 각 행은 일차방정식이라 볼 수 있고 각 일차방정식을 identity matrix에 맞는 계수를 가지도록 행연산을 수행한다는 것이죠. 우리가 예전에 연립일차방정식을 풀 때, 상수배를 하거나 그것을 더하거나 뺀다고 해도 해 자체가 바뀌지 않았던 것을 기억하실 겁니다. 기본 행연산들을 $E$라고 할 때, 기본 행연산의 조합으로 $A$의 역행렬을 구할 수 있습니다. 자세한 내용은 이 포스트를 참고 바랍니다.

$AX=I$

$E_1 E_2 ... E_n AX = I E_1 E_2 ... E_n$

$IX = E_1 E_2 ... E_n = A^{-1}$

가우스-조던 소거법의 디테일은 위의 링크를 참고하고 여기에서는 가우스-조던 소거법을 지난 포스트의 행렬 클래스에 덧붙여 파이썬으로 구현해보고자 합니다.

 

역행렬 구현

    def __rtruediv__(self, x):
        '''Computes x/A'''

        import copy
        if self.ncols != self.nrows:
            raise ArithmeticError('matrix not squred')

        A = copy.deepcopy(self)
        B = Matrix.identity(self.ncols, x)
        for c in range(A.ncols):
            ## 1
            for r in range(c+1, A.ncols):
                if abs(A[r,c]) > abs(A[c,c]):
                    A.swap_rows(r,c)
                    B.swap_rows(r,c)
            p = 0. + A[c,c]  # trick to make sure it is float
            ## 2
            for k in range(A.ncols):
                A[c,k] = A[c,k] / p
                B[c,k] = B[c,k] / p
            ## 3
            for r in list(range(0, c)) + list(range(c+1, A.ncols)):  # except 'c' column
                p = 0. + A[r,c]
                for k in range(A.ncols):
                    A[r,k] -= A[c,k]*p
                    B[r,k] -= B[c,k]*p

        return B

    def __truediv__(self, A):
        if isinstance(A, Matrix):
            return self*(1./A)  # matrix / matrix
        else:
            return (1./A)*self  # matrix / scalar
            
            
A = Matrix([[1,2],[4,9]])
print(1 / A)
>>> [[9.0, -2.0], [-4.0, 1.0]]

print(A / A)
>>> [[1.0, 0.0], [0.0, 1.0]]

print(A / 2)
>>> [[0.5, 1.0], [2.0, 4.5]]

• 파이썬 3부터는 __div__가 __truediv__로 변경되었습니다.
• 가우스-조던 알고리즘대로 컬럼 별로 순회하면서 1번, 2번, 3번 과정을 거칩니다.
• 1번 과정은 "c" 칼럼의 행을 돌면서 밑의 행의 칼럼이 값이 크면 바꿔주는 로직입니다. 없어도 상관 없습니다.
• 2번 과정은 A의 대각 성분이 1이 되도록 해당 행의 값을 대각 성분으로 나눠주는 과정입니다.
• 3번 과정은 "c" 칼럼의 대각 성분이 아닌 성분의 값을 0으로 만들도록 소거해주는 과정입니다.

 

전치 행렬 구현

추가적으로 전치 행렬을 구현합니다. 전치 행렬은 행렬의 행과 열의 개수가 서로 바뀌어 정의됩니다. 해당 행렬 클래스의 속성으로 호출되도록 @property를 이용해 구현합니다.

    @property
    def T(self):                    
        return Matrix(self.ncols, self.nrows, fill=lambda r, c: A[c,r])
        
        
A = Matrix([[1,2],[4,9]])
print(A.T)
>>> [[1, 4], [2, 9]]

 

마지막으로 주어진 정방행렬이 대칭행렬인지 아닌지 판단하는 함수를 구현합니다.

def is_almost_symmetric(A, ap=1e-6, rp=1e-4):
    if A.nrows != A.ncols: return False

    for r in range(A.nrows):
        for c in range(r):
            delta = abs(A[r,c] - A[c,r])
            if delta > ap and delta > max(abs(A[r,c]), abs(A[c,r]))*rp:
                return False

    return True

 

선형방정식 풀기

이제는 우리가 만든 행렬 클래스를 이용해 간단한 선형방정식을 풀 수 있습니다. 다음 방정식에 대하여 $Ax=b$ 꼴로 행렬을 나타낸 뒤에 $x=A^{-1}*b$ 수식으로 $x=(x_0,x_1,x_2)$ 값을 구할 수 있습니다.

$x_0 + 2x_1 + 2x_2=3$

$4x_0+4x_1+4x_2=6$

$4x_0+6x_1+4x_2=10$

A = Matrix([[1,2,2],[4,4,2],[4,6,4]])
b = Matrix([[3],[6],[10]])
x = (1/A)*b
print(x)
>>> [[-1.0], [3.0], [-1.0]]

 

---

[알고리즘 & 코딩테스트/코딩테스트] - 파이썬으로 행렬 구현하기 (4) - 마무리