LU Decomposition
정칙 정방 행렬 (Non-singular Square Matrix) 은 분해가 가능하다.
분해는
형태의 Linear System 을 직접적으로 (directly) 푸는 효율적인 방법이다. 일단
분해를 해놓으면
의 값에 상관없이 두번의 후 대입 (back substitution) 만으로 해를 구할 수 있다.
분해는 유일한 형태가 아니며, 그 형태에 따라 Doolittle, Crout, Cholesky 분해 알고리즘으로 나눠진다. 먼저 일반적인 형태를 알아보고 각각의 알고리즘을 자세히 살펴보자.
가 정칙 정방 행렬이고 행 교환이 필요없다고 하면 아래와 같은 일반적인 형태로
분해가 가능하다.
이것을 상삼각 성분과 하삼각 성분으로 나누면
약간 고치면,
(1)
(2)
여기서 의 대각성분이 모두
이면 Doolittle 분해,
의 대각성분이 모두
이면 Crout 분해라고 한다. 쉽게 생각하면
로 분해가 되는데, 대각행렬
가 왼쪽에 붙는지 오른쪽에 붙는지에 따라 알고리즘이 달라진다고 보면 된다. 그리고 만약
가 대칭행렬 (또는 hermitian 행렬) 이라면
로 분해가 가능하다. 이때
의 대각 성분이 양수라면
형태로 분해가 가능하고 이것을 Cholesky 분해라고 부른다. 세가지 알고리즘 모두 in-place 로 구할 수 있으므로 메모리를 절약할 수 있다.
Doolittle Algorithm
Doolittle 알고리즘은 식 (1), (2) 에서 이므로,
를 구하기 위해서는
성분들이 필요하다. 또
를 구하려면
성분들이 필요하다. 그림으로 나타내면,
파란색으로 표시된 성분을 구하려면 빨간색 성분들이 필요하다. 잘 보면 좌상단부터 차례로 계산해 나갈 수 있다.
따라서 in-place 로 들을 계산해 나가려면 좌상단 성분부터 차례로 한 행씩 순서대로 구하면 된다. 의사 코드는 다음과 같다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
decomp_doolittle(A) for i = 1, ..., n % Compute a line of L for j = 1, ..., i - 1 sum = 0 for k = 1, ..., j - 1 sum += a(i, k) * a(k, j) end a(i, j) = (a(i, j) - sum) / a(j, j) end % Compute a line of U for j = i, ..., n sum = 0 for k = 1, ..., i - 1 sum += a(i, k) * a(k, j) end a(i, j) = a(i, j) - sum end end end |
한 행씩 계산하므로 row major 행렬에 적합하다.
다음은 를 푸는 의사 코드다. in-place 로 계산되었기 때문에
에
가 모두 포함되어있다고 가정한다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
solve_doolittle(A, b) % Solve Ly = b for i = 1, ..., n sum = 0 for j = 1, ..., i - 1 sum += a(i, j) * y(j) end y(i) = b(i) - sum end % Solve Ux = y for i = n, ..., 1 sum = 0 for j = i + 1, ..., n sum += a(i, j) * x(j) end x(i) = (y(i) - sum) / a(i, i) end end |
Crout Algorithm
Crout 알고리즘은 식 (1), (2) 에서 이므로,
를 구하기 위해서는
성분들이 필요하다. 또
를 구하려면
성분들이 필요하다. 그림으로 나타내면,
따라서 in-place 로 들을 계산해 나가려면 좌상단 성분부터 한 열씩 구하면 된다. 의사 코드는 다음과 같다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
decomp_crout(A) for j = 1, ..., n % Compute a column of U for i = 1, ..., j - 1 sum = 0 for k = 1, ..., i - 1 sum += a(i, k) * a(k, j) end a(i, j) = (a(i, j) - sum) / a(i, i) end % Compute a column of L for i = j, ..., n sum = 0 for k = 1, ..., j - 1 sum += a(i, k) * a(k, j) end a(i, j) = a(i, j) - sum end end end |
한 열씩 계산하므로 column major 행렬에 적합하다.
다음은 를 푸는 의사 코드다. in-place 로 계산되었기 때문에
에
가 모두 포함되어있다고 가정한다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
solve_crout(A, b) % Solve Ly = b for i = 1, ..., n sum = 0 for j = 1, ..., i - 1 sum += a(i, j) * y(j) end y(i) = (b(i) - sum) / a(i, i) end % Solve Ux = y for i = n, ..., 1 sum = 0 for j = i + 1, ..., n sum += a(i, j) * x(j) end x(i) = y(i) - sum end end |
Cholesky Algorithm
Cholesky 알고리즘은 설명했듯이 가 대칭행렬인 것을 가정한다.
이므로 식 (1), (2) 는 다음과 같이 바꿀 수 있다.
in-place 로 을 구할 것인지
를 구할 것인지에 따라 두가지 버젼으로 나눌 수 있다.
을 구할 때는 열 우선으로,
를 구할 때는 행 우선으로 구 할 수 있다.
먼저 을 구하는 Cholesky 알고리즘은 다음과 같다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
decomp_lcholesky(A) for j = 1, ..., n % Compute a diagonal element of L sum = 0 for k = 1, ..., j - 1 sum += a(j, k) * a(j, k) end a(j, j) = sqrt(a(j, j) - sum) % Compute a column which is below the diagonal of L for i = j + 1, ..., n sum = 0 for k = 1, ..., j - 1 sum += a(i, k) * a(j, k) end a(i, j) = (a(i, j) - sum) / a(j, j) end end end |
에
성분만 in-place 로 구했기 때문에 대각성분을 제외한 상삼각성분은 그대로 남아있다.
다음은 를 푸는 의사 코드다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
solve_lcholesky(A, b) % Solve Ly = b for i = 1, ..., n sum = 0 for j = 1, ..., i - 1 sum += a(i, j) * y(j) end y(i) = (b(i) - sum) / a(i, i) end % Solve L^Tx = y for i = n, ..., 1 sum = 0 for j = i + 1, ..., n sum += a(j, i) * x(j) end x(i) = (y(i) - sum) / a(i, i) end end |
다음으로 를 구하는 Cholesky 알고리즘은 다음과 같다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
decomp_ucholesky(A) for i = 1, ..., n % Compute a diagonal element of U sum = 0 for k = 1, ..., i - 1 sum += a(k, i) * a(k, i) end a(i, i) = sqrt(a(i, i) - sum) % Compute a row which is right the diagonal of U for j = i + 1, ..., n sum = 0 for k = 1, ..., i - 1 sum += a(k, i) * a(k, j) end a(i, j) = (a(i j) - sum) / a(i, i) end end |
마찬가지로 에
성분만 in-place 로 구했기 때문에 대각성분을 제외한 하삼각성분은 그대로 남아있다.
다음은 를 푸는 의사 코드다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
solve_ucholesky(A, b) % Solve U^Ty = b for i = 1, ..., n sum = 0 for j = 1, ..., i - 1 sum += a(j, i) * y(j) end y(i) = (b(i) - sum) / a(i, i) end % Solve Ux = y for i = n, ..., 1 sum = 0 for j = i + 1, ..., n sum += a(i, j) * x(j) end x(i) = (y(i) - sum) / a(i, i) end end |