3중 대각행렬의 풀이

이번 글에서는 여러 형태의 행렬 중 3중 대각 행렬(tridiagonal matrix)에 대해 다뤄보겠습니다.

 

행렬의 원소 대부분의 값이 0인 행렬을 희소행렬(sparse matrix)이라 합니다. 원소 몇 개를 제외하고는 모두 0이기 때문에 일반적인 행렬보다 다루기가 쉽고,  희소 행렬이 쓰인 연립방정식의 풀이를 좀 더 효율적으로 할 수 있습니다.

 

 

삼중 대각행렬이란?

희소 행렬 중 널리 쓰이는 형태 중 하나가 바로 3중 대각 행렬입니다. 이름에서도 알 수 있듯이 행렬의 대각선과 그 대각선 바로 위, 아래에 위치한 대각선을 제외하고 나머지 원소가 0 인 행렬을 뜻합니다. 수식으로 써보면

 

$$
\mathbf{T}= \begin{pmatrix}
b_1 & c_1 &  & & \\
a_2 & b_2 & c_2 & & \\
    & a_3 & b_3 &\ddots & \\
  & & \ddots&\ddots &c_{n-1} \\
 & & & a_n & b_n
\end{pmatrix}
$$

 

와 같은 형태입니다. 위 행렬의 비어있는 세 대각선에 위치한 원소를 제외하고는 모두 0입니다.

$$ \mathbf{x} = \pmatrix{x_1\\x_2\\ \vdots \\x_n} , \mathbf{d} = \pmatrix{d_1\\d_2\\ \vdots \\d_n}$$

라 할 때, 행렬 방정식

$$ \mathbf{Tx} = \mathbf{d}\tag{1}$$

를 푸는 경우가 종종 등장합니다. 

 

삼중대각행렬 방정식의 풀이

식(1)의  $\mathbf{T}$ 가 삼중대각행렬인지 여부에 상관없이 $\mathbf{T}^{-1}$이 존재한다면

$$ \mathbf{x } = \mathbf {T}^{-1} \mathbf{d} \tag{2}$$

로 풀 수 있습니다. 하지만 역행렬을 계산하는 과정은 복잡할 뿐 아니라 계산시간도 많이 소요되지요. 더군다나 $\mathbf{T}$처럼 거의 모든 원소가 0인 경우는 굳이 역행렬을 구하지 않더라도 간단한 알고리즘이 있을 듯합니다.

그 알고리즘이 바로 Thomas Algorithm입니다.

 

자세한 수학적 원리는 여기를 참고하시기 바랍니다. 이 곳에 python code도 준비되어 있군요! 

 

def thomas(a, b, c, d):
    """ A is the tridiagnonal coefficient matrix and d is the RHS matrix"""
    """
    a is lower diagonal a2,a3,..,a_N, meaning
    b is diagonal b1,b2,b3,..,b_N meaning
    c is upper diagonal c1,c2,c3,.. c_{N-1} meaning
    """
    N = len(a)
    cp = np.zeros(N, dtype='float64')  # store tranformed c or c'
    dp = np.zeros(N, dtype='float64')  # store transformed d or d'
    X = np.zeros(N, dtype='float64')  # store unknown coefficients

    # Perform Forward Sweep
    # Equation 1 indexed as 0 in python
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    # Equation 2, ..., N (indexed 1 - N-1 in Python)
    for i in np.arange(1, (N), 1):
        dnum = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / dnum
        dp[i] = (d[i] - a[i] * dp[i - 1]) / dnum

    # Perform Back Substitution
    X[(N - 1)] = dp[N - 1]  # Obtain last xn

    for i in np.arange((N - 2), -1, -1):  # use x[i+1] to obtain x[i]
        X[i] = (dp[i]) - (cp[i]) * (X[i + 1])

    return (X)

 

첨자만 주의해서 보면 나머지는 어려울 것이 없습니다. 즉. $\mathbf{T}$ 행렬의 $i$번째 행에는 

$$a_i , b_i , c_i$$

의 원소가 있는 것입니다. 따라서 주 대각선보다 하나 밑인 $\mathbf{a}$ 대각선은 $a_1$은 무시하고

$a_2, a_3, \cdots, a_n$으로 이루어져 있고, 주 대각선보다 위에 위치한 $\mathbf{c}$ 대각선은 $c_n$은 무시하고

$c_1, c_2, \cdots, c_{n-1}$로 이루어집니다.

 

 

역행렬 풀이와 Thomas Algorithm 풀이 비교

 

식 (1)을 역행렬로 푸는 방법, 즉 (2)의 방법을 이용하는 것과 Thomas Algorithm을 이용하는 방법이

• 결과는 같게 나오는지
• 계산시간은 각각 어떻게 되는지

를 알아보겠습니다.

 

import numpy as np
import time

def inv_vs_thomas_algorithm():
    nSize = 10
    a = np.random.randn(nSize)
    b = np.random.randn(nSize)
    c = np.random.randn(nSize)
    d = np.random.randn(nSize)

    trid= np.zeros((nSize,nSize))

    for i in range(nSize):
        trid[i,i] = b[i]
        if i>0:
            trid[i,i-1]=a[i]
            trid[i-1,i]=c[i-1]

    x_inv = np.linalg.inv(trid) @ d
    x_thomas = thomas(a,b,c,d)

    print(x_inv- x_thomas)
    print(np.allclose(x_inv, x_thomas))

    nIteration = 100000
    start_time = time.time()
    for i in range(nIteration):
        x_inv = np.linalg.inv(trid) @ d
    time1= time.time() - start_time

    start_time = time.time()
    for i in range(nIteration):
        x_inv = thomas(a,b,c,d)
    time2 = time.time() - start_time

    print('elapsed time using inverse matrix : {:.2f}'.format(time1))
    print('elapse time using Thomas algorithem : {:.2f}'.format(time2))

 

code를 간단히 살펴보겠습니다.

 

import time	# 계산시간을 측정하기 위해 필요한 library

 

    nSize = 10					# 10 × 10 tridiagonal을 만들 예정
    a = np.random.randn(nSize)	# 주대각선 바로 아래 대각선 원소 배열
    b = np.random.randn(nSize)	# 주대각선 원소 배열
    c = np.random.randn(nSize)	# 주대각선 바로 위 대각선 원소 배열
    d = np.random.randn(nSize)	# Tx = d, 즉 행렬방정식의 right hand side 배열 

    trid= np.zeros((nSize,nSize))	# 10 × 10 영행렬을 만든 후

    for i in range(nSize):
        trid[i,i] = b[i]			# 주대각선에 b배열 대입
        if i>0:
            trid[i,i-1]=a[i]		# 주대각선 바로 아래 a 배열 대입
            trid[i-1,i]=c[i-1]		# 주대각선 바로 위 c 배열 대입

 

    x_inv = np.linalg.inv(trid) @ d		# inverse matrix을 이용한 풀이
    x_thomas = thomas(a,b,c,d)			# Thomas algorithm 사숑한 풀이

 

 

    print(x_inv- x_thomas)					# 두 결과 차이 출력
    print(np.allclose(x_inv, x_thomas))		# 주어진 tolerance 하에서 두 행렬이 같은지 true/false 판정

○ 부동 소수 계산으로 인해 같은 결과가 나와야 함에도 불구하고 극미하게 결과가 다른 경우가 있습니다.

○ 이때를 위해 주어진 tolerance 하에서 차이를 보는 numpy.allclose라는 함수가 있습니다.

 

 

    nIteration = 100000					# 역행렬이용법, Thomas algorithm 각각 10만번씩 계산예정
    start_time = time.time()			# 계산시간 측정을 위한 start time 저장
    for i in range(nIteration):
        x_inv = np.linalg.inv(trid) @ d	# 역행렬을 계산하여 행렬방정식을 품
    time1= time.time() - start_time		# 계산종료시점과 start time간의 차이로 10만번 계산시간 측정

    start_time = time.time()
    for i in range(nIteration):
        x_inv = thomas(a,b,c,d)			# Thomas algorithm의 10만번 계산소요시간 테스트
    time2 = time.time() - start_time

 

 

결과

 

[-3.33066907e-16 -5.55111512e-16  2.22044605e-16 -4.44089210e-16
 -1.11022302e-16  1.11022302e-16  5.55111512e-17  1.11022302e-16
 -3.88578059e-16  2.22044605e-16]
 # a,b,c,d가 랜덤으로 추출한 숫자이긴 하지만 두 방법 결과의 차이는 극미하다.(원래 차이가 0이어야 정상)
True	#allclose로 두 결과가 같은지 여부를 판단
elapsed time using inverse matrix : 25.10	
elapse time using Thomas algorithem : 6.53

Process finished with exit code 0

○ 계산 차이가 확연히 나죠. 역행렬 구하는 시간이 오래 걸리기 때문입니다. Thomas algorighm을 뜨면  계산 시간이 거의 1/4로 줄어드는 느낌입니다.

 

3중 대각 행렬의 응용

미분방정식을 푸는 방법 중 FDM이라는 방법이 있습니다. FDM에 대해서는

2022.07.22 - [수학의 재미/아름다운 이론] - 미분방정식을 연립방정식으로, FDM #1

미분방정식을 연립방정식으로, FDM #1

의 FDM 시리즈를 보시면 되겠습니다.  이 방법을 사용하여 문제를 풀 때, 3중 대각 행렬이 많이 등장합니다. 전의 글에서는 3중 대각 행렬의 inverse matrix를 구해서 결과를 얻었는데, 이렇게 되면 시간 소모적입니다.

FDM을 풀 때, 여러 번의 3중 대각 행렬이 나오거든요. 따라서 Thomas algorithm을 사용하는 편이 훨씬 효과적입니다.

앞으로 FDM을 풀 때는 위에 만들어놓은

 

def thomas(a, b, c, d):

을 적극 이용해야겠습니다.