가우스 소거법

미지수가 $n$개인 연립일차방정식

$\displaystyle \begin{cases} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 \\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2 \\ \quad \vdots \\a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n \end{cases}$

은 다음과 같은 행렬로 나타낼 때 $Ax=b$이다. 행렬 곱셈을 해 보면 같은 방정식을 얻을 수 있다.

$\displaystyle \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n\end{bmatrix}=\begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n\end{bmatrix}$

이 방정식의 해 $x=(x_1,\ x_2,\ \cdots,\ x_n)$를 얻기 위해서 다음과 같은 과정을 따를 수 있다.

• 두 번째 방정식에서, 첫 번째 방정식의 양변에 $a_{21}/a_{11}$을 곱한 것을 뺀다. 즉 두 번째 방정식을 이 결과로 대체한다. 그렇다면 두 번째 방정식에 있는 $x_1$ 항이 사라진다.
• 세 번째 방정식에서, 첫 번째 방정식의 양변에 $a_{31}/a_{11}$을 곱한 것을 뺀다. 즉 세 번째 방정식을 이 결과로 대체한다. 그렇다면 세 번째 방정식에 있는 $x_1$ 항이 사라진다.
• 이 과정을 마지막 방정식까지 반복한다. 그러면 첫 번째 방정식을 제외하고 $x_1$ 항이 사라진다.
• 세 번째 방정식에서, 두 번째 방정식의 양변에 $a'_{32}/a'_{22}$를 곱한 것을 뺀다. 즉 세 번째 방정식을 이 결과로 대체한다. 그렇다면 세 번째 방정식에 있는 $x_2$ 항이 사라진다.
• 네 번째 방정식에서, 두 번째 방정식의 양변에 $a'_{42}/a'_{22}$를 곱한 것을 뺀다. 즉 네 번째 방정식을 이 결과로 대체한다. 그렇다면 네 번째 방정식에 있는 $x_2$ 항이 사라진다.
• 이 과정을 마지막 방정식까지 반복한다. 그러면 두 번째 방정식을 제외하고 $x_2$ 항이 사라진다.
• 위와 같은 과정을 $x_n$ 항까지 반복한다.

이를 forward elimination(전진 소거)이라고 한다. 이제 행렬은 다음과 같은 upper triangular matrix가 된다:

$\displaystyle \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ 0 & a'_{22} & \cdots & a'_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a^{'\cdots'}_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n\end{bmatrix}=\begin{bmatrix} b_1 \\ b'_2 \\ \vdots \\ b^{'\cdots'}_n\end{bmatrix}$

가장 아래에서 $x_n=b^{'\cdots'}_{n}/a^{'\cdots'}_{nn}$를 바로 얻고, 이를 위로 대입해 나가면서 차례로 $x_{n-1},\ \cdots,\ x_1$까지 구하는 back substitution(후진 대입)을 통해 연립일차방정식의 해를 구할 수 있다. 이 알고리즘을 Gaussian elimination(가우스 소거법)이라고 한다.

예외 처리

$a_{11},\ a'_{22}$와 같은 나누는 수, 즉 pivot에서 $0$이 등장한다면, 뒤의 행에서 해당 열에 $0$이 아닌 경우가 있을 경우에 그 행과 치환하여 계산할 수 있다. 그러나 뒤의 행에서도 해당 열이 모두 $0$인 경우에, 이 연립방정식은 적당한 상수를 곱하면 좌변이 똑같아지는 행을 두 개 이상 가진다. 예를 들어 $2x+3y$와 $4x+6y$ 같은 두 좌변이 있는 것인데, 우변을 보아서 $4$와 $8$처럼 두 식이 같게 되면 연립방정식은 해가 무수히 많은 것이며, $4$와 $6$처럼 두 식이 다르게 되면 연립방정식은 해가 없는 것이다.

편의상 성분 $a_{ij}$에서 소거한 결과를 $a^{'\cdots'}_{ij},\ b^{'\cdots'}_{i}$와 같이 작성하였는데, 이를 복잡하게 $a_{ij},\ b_i$로 전개해 놓으면 pivot들의 곱이 되는 분모는 행렬 $A$의 determinant(행렬식)이고, 그 전개식은 $n\times n$ 행렬의 determinant를 구하는 일반적인 공식이 된다. 따라서 pivot들 중에 $0$이 있으면 $\det A=0$이며 이때 행렬 $A$를 singular matrix(특이 행렬)라고 한다.

행렬 $A$에 대해서만 생각할 때 첫 번째 방정식을 제외하고 $x_1$항이 사라지게 하려면 첫 번째 방정식을 제외한 행렬의 모든 성분, 즉 $n^2-n$개의 값을 바꿔야 한다. 이와 같은 과정을 마지막 방정식까지 반복하므로 Gaussian elimination의 연산 비용은 좌변에서 $\displaystyle n^2-n+(n-1)^2-(n-1)+\cdots+1^2-1=\frac{n^3-n}{3}$이다. 이를 최적화할 때는 다음과 같은 numerical instability를 고려하는 것이 좋다.

• $Ax=b$의 $A$가 singular matrix에 근접할 때, $b$의 작은 변화가 $x$의 큰 변화를 유도할 수 있다.
• 각 행의 같은 열에서 큰 pivot을 고르기 위한 행 교환을 partial pivoting이라 하고, 열 교환까지 고려하면 complete pivoting이라고 한다.
• positive definiteness, diagonal dominance

LU 분해

가우스 소거법의 각 과정을 주어진 행렬에 다른 행렬들이 작용한 결과로 묘사해 보자. 예를 들어 행렬의 곱셈

$\displaystyle \begin{bmatrix} 1 & 0 & 0 \\ -k & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33}\end{bmatrix}=\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ -ka_{11}+a_{21} & -ka_{12}+a_{22} & -ka_{13}+a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}$

은 미지수가 세 개인 연립일차방정식 $Ax=b$가 주어졌을 때 첫 번째 방정식에 $k$를 곱해 두 번째 방정식에서 빼는 경우를 좌변에 한해서 보여 주고 있다. 즉 $i\neq j$일 때 $E_{ij}=k$인 행렬 $E$를 곱한다는 것은 $j$번째 행에 $k$를 곱한 것을 $i$번째 행에 더하겠다는 뜻이다. 이런 종류의 행렬 $E$를 elementary matrix(기본 행렬)라고 한다.

따라서 전진 소거로 $Ax=b$를 $Ux=c$로 만들었다면, 적당한 기본 행렬 $E_1,\ E_2,\ E_3,\ \cdots,\ E_{(n^2-n)/2}$들이 있어서 $U=E_{(n^2-n)/2}\cdots E_3E_2E_1A$와 같이 쓸 수 있다. 이때 같은 열의 항을 소거하는 두 행렬은 교환 법칙이 성립하는 두 행렬이지만, 다른 열의 항을 소거하는 두 행렬은 교환 법칙이 성립하지 않는 두 행렬이 된다. 이제 $Ux=c$에서 $Ax=b$를 얻는 방법을 생각해 보자. $j$번째 행에 $k$를 곱한 것을 $i$번째 행에 더하는 행위를 상쇄하려면 $j$번째 행에 $k$를 곱한 것을 $i$번째 행에서 빼야 한다. 즉 $k$의 부호만 바꾸면 된다. 이렇게 만든 역행렬들 $E_1^{-1},\ E_2^{-1},\ \cdots E_{(n^2-n)/2}^{-1}$에 대해서 $A=E_1^{-1}E_2^{-1}\cdots \cdots E_{(n^2-n)/2}^{-1}U$이다. 이 식을 다음과 같이 쓴다:

$A=LU$

$L$은 lower triangular matrix이고 $U$는 upper triangular matrix이다. 행 치환이 필요하지 않은 행렬 $A$에 대해서 LU decomposition(LU 분해)을 할 수 있다. LU 분해가 주어지면 $Lc=b$에서 전진 대입으로 $c$를 구하고 $Ux=c$에서 $x$를 구할 수 있다. 그러면 $n^2$의 연산 비용으로 $Ax=b$의 해를 구할 수 있고, 좌변과 우변을 함께 생각할 필요가 없다.

치환 행렬

$k\neq 0$일 때 $E_{ii}=1/k$인 행렬은 $i$번째 행을 $k$로 나눌 때 쓴다. $E_{ii}=0$은 단위 행렬에서 치환할 두 행을 바꾼 것만 허용한다. 이렇게 세 가지의 elementary matrix가 있으며 $n\times n$ 행렬마다 $n!$개의 permutation matrix(치환 행렬)가 있다.

$\displaystyle I=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix},\ P_{12}=\begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix},\ P_{13}=\begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \\ \end{bmatrix},\ P_{23}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{bmatrix},\ P_{23}P_{12}=\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ \end{bmatrix},\ P_{12}P_{23}=\begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}$

cycle notation으로 나타내면 $(1)(2)(3),\ (2\ 1)(3),\ (2)(3\ 1),\ (1)(3\ 2),\ (3\ 1\ 2),\ (3\ 2\ 1)$이다. 행 치환이 필요한 행렬 $A$가 singular matrix가 아니면 치환 행렬 $P$에 대해서 LU decompostion을 $PA=LU$로 쓸 수 있다. 행 교환의 횟수는 $\det A$의 부호를 결정하며, 홀수일 때 음수이다.

대칭 행렬

$U$는 $L$에 비해 무거우므로 $U_{ii}$를 모두 $1$로 만들어 줄 수 있다. $D_{ii}$를 각 pivot으로 놓은 대각 행렬 $D$와, $U$의 각 행을 pivot으로 나눈 행렬 $U'$에 대해서 $U=DU'$이므로 $A=LDU'$로 쓸 수 있다. 따라서 하삼각행렬을 $LD$로 하여 $Ax=b$의 해를 구할 수도 있다.

행 치환이 필요하지 않은 행렬에 대해서 LDU 분해는 유일하므로 그러한 행렬 $A$가 대칭 행렬이면, 전치 행렬의 성질에 의해서 $A=LDL^{T}$이다. 대체로 소거법을 사용하는 거대한 행렬은 symmetric sparse matrix이다. 예를 들어 미분 방정식 $f''(x)=g(x)$의 이산적인 해 $f(0),\ f(h),\ \cdots,\ f(nh)=f(1)$를 구하려면 $f(h),\ \cdots,\ f((n-1)h)$를 해로 가지는 연립일차방정식을 구성하여야 한다. 우선 $f''(x)$로부터 일차항과 상수항을 알 수 없으므로 boundary condition $f(0)=0,\ f(1)=0$을 지정하자. 이계 도함수를

$\displaystyle g(x)=\lim_{h\to 0}\frac{1}{h}\left(\frac{f(x+h)-f(x)}{h}-\frac{f(x)-f(x-h)}{h}\right)$

로 쓸 수 있으므로 $-h^2g(ih)\approx -f((i+1)h)+2f(ih)-f((i-1)h)$로 근사할 수 있다. $h$와 $g(ih)$는 주어져 있으므로

$\displaystyle \begin{bmatrix} 2 & -1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & \cdots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} f(h) \\ f(2h) \\ f(3h) \\ \vdots \\ f((n-1)h)\end{bmatrix}=\begin{bmatrix} -h^2g(h) \\ -h^2g(2h) \\ -h^2g(3h) \\ \vdots \\ -h^2g((n-1)h)\end{bmatrix}$

로 나타낼 수 있고 LDU 분해하면 $\displaystyle L_{(i+1)i}=-\frac{i}{i+1},\ D_{ii}=\frac{i+1}{i},\ U_{i(i+1)}=-\frac{i}{i+1}$이다. 연산 비용은 $2n$이 된다.

가우스-요르단 소거법

$Ax=b$를 풀면 $x=A^{-1}b$이며, 이는 $AA^{-1}=I$인 $A^{-1}$을 구하는 것과 유사하다. 예를 들어 세 개의 연립일차방정식

$\displaystyle \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3\end{bmatrix}=\begin{bmatrix} 1 \\ 0 \\ 0\end{bmatrix},\ \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3\end{bmatrix}=\begin{bmatrix} 0 \\ 1 \\ 0\end{bmatrix},\ \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3\end{bmatrix}=\begin{bmatrix} 0 \\ 0 \\ 1\end{bmatrix}$

을 풀어 각 해가 각 열을 이루도록 붙이면 역행렬을 얻는다. 이들을 한꺼번에 풀기 위해서, 방정식 $AX=I$를 두고 다음과 같은 과정을 따를 수 있다.

• 전진 소거 과정을 가우스 소거법과 동일하게 적용한다. 그러면 우변의 항등 행렬은 기본 행렬들의 곱 $E_{(n^2-n)/2}\cdots E_3E_2E_1I$이 되고, 이 식은 $UX=L^{-1}$이다.
• $n-1$번째 방정식에서, $n$번째 방정식의 양변에 적당한 상수를 곱한 것을 뺀다. 그렇다면 $n-1$번째 방정식의 마지막 항이 사라진다.
• 이 과정을 첫 번째 방정식까지 반복한다. 그러면 마지막 방정식을 제외하고 마지막 항이 사라진다.
• 마지막 항에 적용한 과정을 첫 번째 항까지 반복한다.
• $A$의 대각 성분으로 우변의 각 행을 나누면 $IX=U^{-1}L^{-1}$을 얻는다.

이 알고리즘은 Gauss–Jordan elimination(가우스-요르단 소거법)이다. 역행렬을 구하는 연산 비용은 $A^2$를 계산하는 연산 비용과 같다.

참고 자료

• Gilbert Strang. Linear Algebra and Its Applications.