LU Decomposition Algorithm

LU Decomposition

정칙 정방 행렬 (Non-singular Square Matrix) 은 LU 분해가 가능하다. LU 분해는 Ax = b 형태의 Linear System 을 직접적으로 (directly) 푸는 효율적인 방법이다. 일단 LU 분해를 해놓으면 b 의 값에 상관없이 두번의 후 대입 (back substitution) 만으로 해를 구할 수 있다. LU 분해는 유일한 형태가 아니며, 그 형태에 따라 Doolittle, Crout, Cholesky 분해 알고리즘으로 나눠진다. 먼저 일반적인 형태를 알아보고 각각의 알고리즘을 자세히 살펴보자.

A 가 정칙 정방 행렬이고 행 교환이 필요없다고 하면 아래와 같은 일반적인 형태로 LU 분해가 가능하다.

    \[A = LU = \begin{bmatrix} l_{11} & 0 & \dots & 0 \\ l_{21} & l_{22} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \dots & l_{nn} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n} \\ 0 & u_{22} & \dots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & u_{nn} \end{bmatrix}\]

    \[a_{ij} = \sum_{k=1}^{\text{min}(i, j)} l_{ik} u_{kj}\]

이것을 상삼각 성분과 하삼각 성분으로 나누면

    \[\begin{aligned} a_{ij} &= \sum_{k=1}^i l_{ik} u_{kj}\quad \text{s.t. } i \le j \\ a_{ij} &= \sum_{k=1}^j l_{ik} u_{kj}\quad \text{s.t. } i \ge j \\ \end{aligned}\]

약간 고치면,

(1)   \[\begin{aligned}  a_{ij} &= \sum_{k=1}^{i-1} l_{ik} u_{kj} + l_{ii} u_{ij} \quad \text{s.t. } i \le j \\ u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}}{l_{ii}} \\ \end{aligned}\]

(2)   \[\begin{aligned}  a_{ij} &= \sum_{k=1}^{j-1} l_{ik} u_{kj} + l_{ij} u_{jj} \quad \text{s.t. } i \ge j \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}} \\ \end{aligned}\]

여기서 L 의 대각성분이 모두 1 이면 Doolittle 분해, U 의 대각성분이 모두 1 이면 Crout 분해라고 한다. 쉽게 생각하면 A = LDU 로 분해가 되는데, 대각행렬 D 가 왼쪽에 붙는지 오른쪽에 붙는지에 따라 알고리즘이 달라진다고 보면 된다. 그리고 만약 A 가 대칭행렬 (또는 hermitian 행렬) 이라면 A = LDL^T 로 분해가 가능하다. 이때 D 의 대각 성분이 양수라면 A = LL^T 형태로 분해가 가능하고 이것을 Cholesky 분해라고 부른다. 세가지 알고리즘 모두 in-place 로 구할 수 있으므로 메모리를 절약할 수 있다.

Doolittle Algorithm

Doolittle 알고리즘은 식 (1), (2) 에서 l_{ii} = 1 이므로,

    \[\begin{aligned} u_{ij} &= a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj} \quad &\text{s.t. } i \le j \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}} \quad &\text{s.t. } i > j \\ \end{aligned}\]

u_{ij} 를 구하기 위해서는 l_{i1}, l_{i2}, \dots, l_{i(i-1)} 성분들이 필요하다. 또 l_{ij} 를 구하려면 u_{1j}, u_{2j}, \dots, u_{jj} 성분들이 필요하다. 그림으로 나타내면,

doolittle-lu-order

파란색으로 표시된 성분을 구하려면 빨간색 성분들이 필요하다. 잘 보면 좌상단부터 차례로 계산해 나갈 수 있다.
따라서 in-place 로 l_{ij}, u_{ij} 들을 계산해 나가려면 좌상단 성분부터 차례로 한 행씩 순서대로 구하면 된다. 의사 코드는 다음과 같다.

한 행씩 계산하므로 row major 행렬에 적합하다.

다음은 LUx = b 를 푸는 의사 코드다. in-place 로 계산되었기 때문에 AL, U 가 모두 포함되어있다고 가정한다.

Crout Algorithm

Crout 알고리즘은 식 (1), (2) 에서 u_{jj} = 1 이므로,

    \[\begin{aligned} u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}}{l_{ii}} \quad &\text{s.t. } i < j \\ l_{ij} &= a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \quad &\text{s.t. } i \ge j \\ \end{aligned}\]

u_{ij} 를 구하기 위해서는 l_{i1}, l_{i2}, \dots, l_{ii} 성분들이 필요하다. 또 l_{ij} 를 구하려면 u_{1j}, u_{2j}, \dots, u_{j(j-1)} 성분들이 필요하다. 그림으로 나타내면,

crout-lu-order

따라서 in-place 로 l_{ij}, u_{ij} 들을 계산해 나가려면 좌상단 성분부터 한 열씩 구하면 된다. 의사 코드는 다음과 같다.

한 열씩 계산하므로 column major 행렬에 적합하다.

다음은 LUx = b 를 푸는 의사 코드다. in-place 로 계산되었기 때문에 AL, U 가 모두 포함되어있다고 가정한다.

Cholesky Algorithm

Cholesky 알고리즘은 설명했듯이 A 가 대칭행렬인 것을 가정한다. l_{ij} = u_{ji} 이므로 식 (1), (2) 는 다음과 같이 바꿀 수 있다.

    \[\begin{aligned} u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} u_{ki} u_{kj}}{u_{ii}} \quad &\text{s.t. } i \le j  \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk}}{l_{jj}} \quad &\text{s.t. } i \ge j  \\ \end{aligned}\]

in-place 로 L 을 구할 것인지 U 를 구할 것인지에 따라 두가지 버젼으로 나눌 수 있다. L 을 구할 때는 열 우선으로, U 를 구할 때는 행 우선으로 구 할 수 있다.

먼저 L 을 구하는 Cholesky 알고리즘은 다음과 같다.

AL 성분만 in-place 로 구했기 때문에 대각성분을 제외한 상삼각성분은 그대로 남아있다.

다음은 LL^Tx = b 를 푸는 의사 코드다.

다음으로 U 를 구하는 Cholesky 알고리즘은 다음과 같다.

마찬가지로 AU 성분만 in-place 로 구했기 때문에 대각성분을 제외한 하삼각성분은 그대로 남아있다.

다음은 U^TUx = b 를 푸는 의사 코드다.

Leave a Reply

Your email address will not be published. Required fields are marked *