LU Decomposition in IDL

LU 분해(LU Decomposition)는 선형대수학에서 정방형 행렬 A를 하삼각 행렬(Lower Triangular Matrix)와 상삼각 행렬(Upper Triangular Matrix) 두 성분으로 쪼개는 기법을 의미합니다. L과 U를 서로 행렬곱을 해주면 원래 행렬인 A를 얻게 됩니다. 개념을 좀 더 도식적으로 본다면 다음 그림과 같습니다.

이러한 작업은 결국 1차 연립방정식의 해를 구하는데 사용이 됩니다. 왜냐하면 미지수의 갯수와 방정식의 갯수가 같을 때 이 식을 연립방정식으로 표현하면, 미지수에 대한 계수들로 정방형의 A라는 행렬이 얻어지는데, 이 A에 대하여 LU 분해를 해주게 되면 결과적으로 미지수들의 해를 얻을 수 있기 때문입니다. IDL에서도 이러한 일련의 작업들을 수행할 수 있도록 해주는 내장 루틴들이 있습니다. 바로 LUDC와 LUSOL로서 제가 예전에 이 블로그에서 그 사용법을 한번 소개해드린 바 있습니다 (링크 참조). 여기서 LUDC라는 루틴은 말 그대로 LU 분해를 수행하는 역할을 하고, 그 결과를 LUSOL이 사용하여 최종적으로 연립방정식의 해를 얻게 됩니다.

그런데 LUDC라는 루틴에 대하여 좀 더 들여다보면, 이 루틴 자체는 LU 분해를 수행하지만 그 결과는 L, U 성분이 완벽히 분리되어서 전달되지 않습니다. 입력인자로 들어간 원본 행렬 A와 구조는 같지만 값들의 구성이 바뀐 형태로 변형된 결과가 전달되는데, 이 결과 역시 하나의 행렬입니다. 대신 이 결과 행렬안에는 L과 U의 성분값들이 한꺼번에 포함됩니다. 도식적으로 본다면 다음 그림과 같은 방식입니다.

물론 LU 분해를 수행한다는 것은 맞지만, 그 결과 자체는 우리가 금방 알아보기 쉬운 L과 U로 얻어지지 않는다는 의미입니다. 따라서 A의 형태로 얻어진 값들을 사용하여 재구성을 해야만 L, U를 따로따로 얻을 수 있다는 얘기가 됩니다. 물론 최종적으로 연립방정식의 해를 구하는 데에만 관심이 있을 경우에는 이것은 별다른 문제는 되지 않습니다. 그러나 L과 U를 따로따로 구하고 싶을 경우에는 LUDC의 결과만으로는 충분하지 않게 느껴질 가능성이 큽니다. 물론 LUDC의 결과를 L, U로 따로 재구성하는 기본적인 방법에 대해서는 수치계산 분야에서 유명한 'Numerical Recipes in C'라는 책의 2장 3절에 나와있긴 합니다만, 유저 입장에서는 이런 내용을 참조하여 직접 해석을 하는 작업이 상당히 귀찮을 수 밖에 없습니다.

그래서 혹시나 L, U를 쉽게 얻을 수 있는 방법에 대하여 IDL 본사인 ExelisVis의 기술지원팀에 문의를 넣어봤는데요. 몇번 문답이 오고간 끝에, 이러한 기능을 수행하는 루틴을 고맙게도 새로 하나 만들어서 보내왔습니다. 테스트해본 결과 잘 작동을 합니다. 루틴의 이름은 LU_DECOMP라는 함수인데, 소스 파일은 이 게시물에 첨부해놓았으니 받아서 사용하시면 됩니다. 사용법 및 결과의 형태에 대한 설명은 다음 그림을 참조하시면 됩니다.

이 루틴을 사용한 예제를 하나 본다면, 다음과 같은 A가 있을 때 이를 LU_DECOMP 함수에 입력인자로 넣으면 됩니다.

A = [[1, 1, 1], $

    [2, 3, 5], $

    [4, 0, 5]]

A_LU = LU_DECOMP(A)

이 때 얻어지는 결과는 구조체(Structure)의 형태로 얻어집니다. 이 구조체안에는 L, U가 독립된 필드로서 저장이 되어 있기 때문에 다음과 같이 개별적으로 조회해볼 수 있습니다.

PRINT, 'Value of Lower:', A_LU.L, /IMPLIED

PRINT, ''

PRINT, 'Value of Upper:', A_LU.U, /IMPLIED

Value of Lower:

       1.0000000       0.0000000       0.0000000

       2.0000000       1.0000000       0.0000000

       4.0000000      -4.0000000       1.0000000

Value of Upper:

       1.0000000       1.0000000       1.0000000

       0.0000000       1.0000000       3.0000000

       0.0000000       0.0000000       13.000000

그리고 다음과 같이 이 두 행렬을 서로 행렬곱으로 곱해주면 원래의 A가 다시 얻어집니다.

PRINT, 'Value of L*U to check that is correct:', A_LU.L ## A_LU.U, /IMPLIED

Value of L*U to check that is correct:

       1.0000000       1.0000000       1.0000000

       2.0000000       3.0000000       5.0000000

       4.0000000       0.0000000       5.0000000

따라서 이 LU_DECOMP라는 루틴을 사용하면 정방형 행렬 A에 대하여 LU 분해를 수행하여 얻어진 L, U 성분을 명확하게 얻을 수있습니다. 단, 유의할 점은 이 LU_DECOMP는 L, U를 따로 볼 수 있도록 결과를 산출해주는 역할만 수행합니다. 만약 연립방정식의 해를 구하는 작업까지 염두에 둘 경우에는 여전히 LUDC, LUSOL을 사용하는 것이 맞다는 점만 염두에 두면 됩니다.

끝으로 이 루틴을 제공해준 ExelisVis 기술팀의 Zachary Norman에게 감사를 표합니다. 제가 이 글 작성하면 이 분에게 링크를 알려드리기로 했기 때문에 아마 이 게시물을 보시게 될 겁니다. 그래서 영문으로도 인사말을 남겨두도록 하겠습니다.

Thank you so much for your kindness and effort, Mr. Norman!

lu_decomp.pro

lu_decomp.pro

0.0MB