가우스 소거법

미지수가 $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를 고려하는 것이 좋다.

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$의 해를 구할 수 있고, 좌변과 우변을 함께 생각할 필요가 없다.

LDU 분해

$U$는 $L$에 비해 무거우므로, $U_{ii}$를 모두 $1$로 만들어 줄 수 있다. $D_{ii}$를 각 pivot으로 놓은 대각 행렬 $D$와, $U$의 각 행을 pivot으로 나눈 행렬 $U'$에 대해서 $U=DU'$이므로 $A=LDU'$로 쓰면 이러한 분해는 유일하다. 이때 하삼각행렬을 무겁게 하여 $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_{(1\ 2)}=\begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix},\ P_{(1\ 3)}=\begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \\ \end{bmatrix},\ P_{(2\ 3)}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{bmatrix},\ P_{(2\ 3)}P_{(1\ 2)}=\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ \end{bmatrix},\ P_{(1\ 2)}P_{(2\ 3)}=\begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}$

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

가우스-요르단 소거법

$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$를 계산하는 연산 비용과 같다.