Cholesky분해는 양정부호행렬을 하삼각행렬과 그 전치행렬의 곱으로 표현하는 것이다라는 내용은 기초 선형대수학 교과서에 나오는 얘기다. 수학과 담을 쌓고 사는 사람이 아니면 늦어도 학부 3학년에는 이것을 배운다. 수학적 중요성은 잠깐 잊어버리더라도 계산상으로 Cholesky분해는 꽤 중요하다. 우선 LU분해에 걸리는 시간의 절반밖에 안걸린다. 게다가, 수치상 안정적이어서 다른 알고리즘에 비해 계산과정상 누적되는 오차의 영향을 덜 받는다. 단점도 물론 있다. 자세한 것은 선형대수 교과서나 수치해석 교과서를 보면 잘 나온다.
주요한 사용 예를 들자면: 최적화 문제의 답을 수치해석적으로 구해야만 하는 경우에 목적함수가 터무니없이 이상하게 생기지 않았으면 가중최소제곱법을 반복해서 답을 구할 수 있다. 이때 가중최소제곱법의 답을 구할때 Cholesky분해를 이용한다. (물론, 다른 방법으로 가중최소제곱법의 답을 구할 수도 있다. 예를 들어 MINPACK은 QR분해를 이용한다.)
실제로 계산하는 법은 의외로 간단하다. 이미 구현되어 제공되는 함수들도 많다.
예를들어 LAPACK의
DPOTRF는
배정도실수 행렬을 다루는 Fortran서브루틴이다.
같은 기능을 하는 함수가 NAG와 IMSL에도 물론 있다.
그런데, 이론적으로
Choleksy분해는
양반정부호행렬까지 확장될 수 있다.
그러나, 대부분의 라이브러리의 Cholesky분해 함수들은 양정부호행렬만 다룬다.
양반정부호행렬을 다루는 구현이
이미 1968년에 공개되었다는 점을 고려하면 이런 상황이 좀 의아하다.
내가 아는한 IMSL의 imsl_f_lin_sol_nonnegdef()만이 양반정부호행렬을
다루는 C 구현이다.
이하의 구현은 위에 링크된 1968년에 공개된 Fortran 구현을 C로 다시 구현하면서
Farebrother and Berry의 버그수정과
그것을 고려한 Barrett and Healy의
개선을 포함한 것이다.
Healy의 원래 구현은 입력행렬과 출력행렬이 같은 주소를 가리켜도 되었는데,
Barrett and Healy의 개선점을 포함하면서 그것이 허용되지 않게 되었다.
양정부호행렬만 분해한다면 이 함수는 LAPACK의
DPPTRF와
같은 기능을 한다.
함수 원형
int cholesky(double A[], int n, double U[], int *nullity)
- A (입력)
- 양반정부호 행렬. 하삼각행렬부분만 1차원배열로 저장한다.
즉 A는
n*(n+1)/2개의 원소를 A11, A21, A22, A31, ... 순으로 저장한 1차원 배열. - n (입력)
- A의 차원. 즉, 수학적으로 A는 n×n행렬.
- U (출력)
- Cholesky분해의 결과로 나오는 하삼각행렬. 역시 1차원 배열로 저장된다.
- nullity (출력)
- A의 영공간의 차원. NULL이면 출력이 없다. NULL을 전달했으나 이후에 영공간의 차원을 알고자 한다면, 대각원소들 중에 0.0인 것들의 갯수를 세면 된다.
- 0
- 작업 성공
- -1
- 비정상적인 매개변수: n <= 0이거나 A와 U가 같은 주소를 가리키는 경우.
- 1 ... n
- A가 양반정부호행렬이 아니며, 리턴값에 해당하는 행에서 처음으로 양반정부호행렬이 아님이 발견됨
구현
수식자체와 양정부호행렬임을 검사하는 방법은 아무 선형대수 교과서나
수치해석 교과서를 봐도 잘 설명되어 있다.
위에 링크된 논문들과 그에 따른 구현에서 중요한 것은 언제 어떻게 선형종속임을
결정하는가 하는 부분이다.
(약간 줄 수가 많아 접어놓음.)
더 해볼만한 것
원 논문의 바로 다음에 이어지는 논문에는
이 함수를 이용해서 양(반)정부호행렬의 (일반화)역행렬을 구하는 함수가 공개되어 있다.
같은 기능을 하는 것을 스스로 만들어쓰려면 BLAS의
DTPSV를 이용하면 된다.
