촐레스키 분해

1. 개요

촐레스키 분해 Cholesky Decomposition
분해 형태	A = LLᵀ
적용 조건	대칭 양정치 행렬 (SPD)
연산량	LU의 약 1/2 (n³/3)
관련 분해	LU 분해, QR 분해
주요 응용	유한요소 강성행렬, 정규방정식

대칭이고 양정치라면, 굳이 두 번 일할 필요 없다.

촐레스키 분해(Cholesky Decomposition)는 대칭 양정치(symmetric positive-definite, SPD) 행렬 $A$를 하삼각행렬 $L$과 그 전치의 곱 $A = LL^{\mathsf T}$로 나타내는 분해다. LU 분해의 특수한 경우이지만, 행렬이 대칭이고 양정치라는 강한 조건을 만족하면 위쪽 삼각형을 따로 저장·계산할 필요가 없어져 연산량과 메모리가 정확히 절반으로 줄어든다.¹ 그래서 SPD 행렬을 만났을 때 촐레스키를 안 쓰는 건 그냥 손해다.

프랑스 육군 측지 장교였던 앙드레루이 촐레스키(André-Louis Cholesky)가 삼각측량 문제를 풀려고 고안했는데, 정작 본인은 제1차 세계대전에서 전사하고 유고로 발표되며 세상에 알려졌다.² 오늘날엔 구조해석의 강성행렬, 최소자승의 정규방정식, 몬테카를로 방법의 상관 샘플링까지, SPD가 나오는 곳이라면 어디서나 등장하는 일꾼이 됐다.

2. 알고리즘과 절반의 연산량

분해식은 $A = LL^{\mathsf T}$를 성분별로 풀어 재귀적으로 얻는다. 대각 성분과 비대각 성분이 각각

$$L_{jj} = \sqrt{A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2}, \qquad L_{ij} = \frac{1}{L_{jj}}\left(A_{ij} - \sum_{k=1}^{j-1} L_{ik}L_{jk}\right)$$

로 계산된다. 총 연산량은 대략 $n^3/3$으로, LU 분해($2n^3/3$)의 딱 절반이다. 대칭성 덕에 절반만 훑으면 되기 때문. 구현에서 곱셈·나눗셈·뺄셈을 직접 세어 보면 선행항 근사가 얼마나 잘 맞는지 바로 보인다.

$n$	촐레스키 실측	$n^3/3$	LU 실측	$2n^3/3$
50	42,925	41,667	82,075	83,333
200	2,686,700	2,666,667	5,313,300	5,333,333
1000	333,833,500	333,333,333	666,166,500	666,666,667

$n$ = 1000에서 두 실측값의 비는 1.9955 — “딱 절반”이라는 말이 선행항 수준에서 정확하다. 작은 $n$ 에서 실측이 약간 큰 것은 $O(n^2)$ 항과 제곱근 $n$ 회 때문이고, $n$ 이 커지면 묻힌다.

게다가 SPD 행렬은 피벗팅 없이도 수치적으로 안정해서, 가우스 소거법에서 골머리를 앓던 부분 피벗팅을 생략할 수 있다는 게 큰 매력이다. 공짜로 되는 게 아니라 SPD라는 조건이 성장 인자(growth factor)를 묶어 주기 때문이고, 후진 오차가 $\kappa(A)$ 에 비례하는 선에서 끝난다는 것이 증명돼 있다. 피벗팅을 생략하면 메모리 접근 패턴이 규칙적이 되어 희소행렬 구조를 분해 전에 미리 확정할 수 있다는 부수 효과도 크다 — 이게 희소 촐레스키가 희소 LU보다 다루기 쉬운 근본 이유다.

다만 대각 성분에 제곱근이 들어가는 게 거슬리면, 제곱근을 피한 $A = LDL^{\mathsf T}$ 변형을 쓰기도 한다. 부동소수점 연산에서 제곱근은 나눗셈급으로 비싸고, 무엇보다 준정치 경계에서 음수의 제곱근을 밟을 위험을 $D$ 의 부호 검사로 바꿀 수 있다.

3. 양정치 판정기

촐레스키 분해에는 숨은 보너스가 있다. 양정치성 판정 도구로 쓸 수 있다는 것. 분해를 진행하다 대각 성분 $L_{jj}$의 제곱근 안이 0 이하가 되면(음수의 제곱근이 나오면), 그 순간 “이 행렬은 양정치가 아니다”라고 확정할 수 있다. 고유값을 다 구해서 부호를 확인하는 것보다 훨씬 싸고 빠르다.³ LAPACK의 dpotrf가 반환하는 INFO 값이 바로 이 판정이다 — 0이면 성공, $j > 0$ 이면 $j$ 번째 선행 주부분행렬에서 양정치성이 깨졌다는 뜻이다.

3.1. 판정에는 수치적 단서가 붙는다

주의할 점은 이 판정이 배정도 산술에서의 판정이라는 것이다. 수학적으로 멀쩡한 SPD 행렬도 조건수가 충분히 나쁘면 촐레스키가 터진다. 고전적인 병적 행렬인 힐베르트 행렬 $H_{ij} = 1/(i+j-1)$ 로 확인해 보면 이렇다. 이 행렬은 모든 $n$ 에서 수학적으로 SPD다.

$n$	8	10	12	13	14
$\kappa(H_n)$	1.5×10¹⁰	1.6×10¹³	1.8×10¹⁶	—	—
배정도 촐레스키	성공	성공	성공	성공	실패

$n$ = 14에서 13번째 대각 성분의 근호 안이 $-5.3\times10^{-16}$ 으로 음수가 되어 멈춘다. 조건수가 배정도 기계 epsilon의 역수($1/\varepsilon \approx 4.5\times10^{15}$)를 넘어서는 지점과 대체로 겹친다. 그러니 “촐레스키가 실패했다”는 신호를 받았을 때 결론은 두 갈래다 — 진짜로 양정치가 아니거나, 양정치이긴 한데 조건수가 배정도로 감당이 안 되거나. 둘을 구분하지 않으면 멀쩡한 모델을 의심하며 며칠을 날린다.⁴

최적화에서 실패를 “수정”하는 수정 촐레스키(modified Cholesky)도 이 지점에서 설계가 갈린다. 대각에 $\tau I$ 를 더해 억지로 밀어 올리는 단순한 방식은 $\tau$ 를 너무 크게 잡으면 뉴턴 방향이 최급강하 방향으로 퇴화한다. 대칭 부정치 분해를 먼저 하고 필요한 최소 보정만 계산하는 쪽이 이 절충을 더 잘 다룬다.

이 성질은 최적화에서 특히 요긴하다. 뉴턴-랩슨법 기반 최적화에서 헤시안 행렬이 양정치인지(즉 현재 방향이 하강 방향인지) 확인하는 데, 촐레스키가 성공하느냐 실패하느냐가 곧 판정이 된다. 실패하면 헤시안을 수정(modified Cholesky)해 강제로 양정치로 만들어 밀고 나가기도 한다.

4. 응용

4.1. 유한요소 강성행렬

유한요소법으로 조립한 선형 정적 구조 문제의 강성행렬 $K$는 (적절한 경계 조건이 걸리면) 대칭 양정치가 된다. 그래서 $Ku = f$를 푸는 표준 직접법이 바로 촐레스키다. 실제 상용 구조해석 솔버(예: NASTRAN)의 직접 솔버 상당수가 희소행렬용 촐레스키(sparse Cholesky)를 기반으로 한다. 여기서는 채우기(fill-in)를 줄이는 재정렬(reordering)과 결합해 대규모 문제를 다룬다. 어떤 순서로 소거하느냐가 연산량을 수십 배 가르므로, 최소차수 순서화(AMD)나 중첩절단(nested dissection) 같은 순서화 알고리즘이 사실상 솔버 성능의 본체다. 전체 구조와 다른 분해와의 비교는 희소 직접 솔버가 따로 다룬다.

4.2. 정규방정식

최소자승법의 정규방정식 $A^{\mathsf T}A\,x = A^{\mathsf T}b$에서 $A^{\mathsf T}A$는 (풀랭크면) SPD다. 촐레스키로 빠르게 풀 수 있지만, 조건수가 제곱으로 나빠진다는 함정이 있어 병적 문제에선 QR 분해나 특이값 분해에 자리를 내준다. “싸고 빠르지만 예민하다”가 요약.

이 함정의 크기는 재 보면 바로 나온다. 40개 점에서 다항식 차수를 올려 가며 반데르몽드 최소자승을 풀고(참해를 전부 1로 심어 두고) 정규방정식+촐레스키와 하우스홀더 QR을 비교하면 이렇다.

미지수 $n$	$\kappa(A)$	$\kappa(A^{\mathsf T}A)$	정규방정식 오차	QR 오차
4	1.2×10²	1.4×10⁴	1.2×10⁻¹²	1.8×10⁻¹⁴
8	1.1×10⁵	1.2×10¹⁰	7.4×10⁻⁷	3.6×10⁻¹²
10	3.5×10⁶	1.2×10¹³	4.2×10⁻⁴	2.6×10⁻¹⁰
12	1.2×10⁸	1.4×10¹⁶	2.6×10⁻¹	3.9×10⁻⁹

$n$ = 12에서 정규방정식의 답은 상대오차 26% 로 사실상 쓸 수 없는데, 같은 문제를 QR로 풀면 아직 $10^{-9}$ 다. $\kappa(A)$ 가 $10^8$ 밖에 안 되는데도 그렇다 — 제곱되어 $10^{16}$ 이 되는 순간 배정도가 바닥나기 때문이다. “조건수 제곱”은 비유가 아니라 유효숫자 절반을 그 자리에서 잃는다는 뜻이다. 반대로 조건수가 온순한 문제($n$ = 4)에서는 정규방정식도 $10^{-12}$ 라 아무 문제가 없으니, 촐레스키를 쓰지 말라는 얘기가 아니라 $\kappa(A)$ 를 먼저 보라는 얘기다.

4.3. 몬테카를로 상관 샘플링

주어진 공분산행렬 $\Sigma$를 따르는 상관 있는 난수를 만드는 표준 레시피가 촐레스키다. $\Sigma = LL^{\mathsf T}$로 분해한 뒤, 독립 표준정규 벡터 $z$에 $L$을 곱하면 $Lz$가 정확히 원하는 공분산 구조를 갖는다.

$$\mathrm{Cov}(Lz) = L\,\mathrm{Cov}(z)\,L^{\mathsf T} = LL^{\mathsf T} = \Sigma$$

몬테카를로 방법 기반 금융 시뮬레이션, 불확실성 정량화, 가우시안 프로세스 회귀가 모두 이 한 줄에 기대고 있다. 가우시안 프로세스에서는 공분산행렬이 커널 길이 척도에 따라 심하게 병적이 되는 일이 흔해서, 대각에 작은 지터(jitter)를 더해 촐레스키를 통과시키는 게 관행이다 — 위에서 본 “양정치인데 조건수가 감당 안 되는” 경우의 가장 흔한 실전 사례다.

5. 현업에서의 촐레스키

• SPD임이 보장되는 순간, 다른 분해는 고민할 이유가 없다. 반값에 안정성까지 공짜다.
• 대각 성분 제곱근 안이 음수로 나왔다면? 축하한다, 당신의 행렬은 양정치가 아니거나 반올림 오차로 준정치 경계에 걸린 것이다. 물성치나 경계 조건, 특이한 구속을 의심하라.⁵
• 대규모 희소 문제에서는 순수 직접법인 촐레스키와 반복법인 크리로프 부분공간법(특히 켤레기울기법, CG) 사이에서 저울질한다. CG의 전처리기로 불완전 촐레스키(incomplete Cholesky)를 쓰는 하이브리드도 흔하다.
• 분해를 혼합 정밀도로 돌리는 것도 선택지다. 단정도로 분해해 두고 반복 정련(iterative refinement)으로 배정도 정확도를 회복하는 패턴인데, 위 힐베르트·반데르몽드 실측이 보여 주듯 조건수가 나쁜 문제에서는 이 회복이 보장되지 않는다. 정밀도를 깎기 전에 $\kappa(A)$ 를 먼저 보는 순서는 여기서도 같다.

6. 관련 문서

• LU 분해 · QR 분해 · 특이값 분해
• 가우스 소거법
• 강성행렬 · 유한요소법 · 구조해석
• 몬테카를로 방법
• 크리로프 부분공간법 · 켤레기울기법 · 전처리기
• 희소행렬 · 희소 직접 솔버 · LAPACK
• 조건수 · 부동소수점 연산
• 가우시안 프로세스 · 불확실성 정량화

7. Footnotes

1. “절반”이 나오는 이유는 순전히 대칭성 때문이다. 위삼각과 아래삼각이 서로 전치라 한쪽만 저장·계산하면 나머지는 공짜로 따라온다. 공짜라는 말에 엔지니어가 약한 건 국룰이다. ↩

2. 촐레스키는 1918년 전사했고, 그의 방법은 동료 브누아(Commandant Benoît)가 1924년에 소개하며 세상에 나왔다. 본인은 자기 이름이 SPD 행렬의 대명사가 될 줄 몰랐을 것이다. ↩

3. 물론 이론적으로는 “모든 선행 주부분행렬식이 양수” 같은 실베스터 판정도 있지만, 실전에서 행렬식을 계산하는 건 촐레스키를 돌리는 것보다 훨씬 비싸고 불안정하다. 그냥 촐레스키를 돌려보고 터지나 안 터지나 보는 게 최고다. ↩

4. 구분하는 실용적인 방법이 있다. 실패한 위치 $j$ 가 어디냐를 보는 것이다. 경계 조건 누락이나 물성 입력 실수는 보통 특정 절점·요소에 대응하는 자리에서 터지고 근호 안이 크게 음수다. 조건수 문제는 분해가 거의 끝까지 간 뒤 근호 안이 $-10^{-16}$ 수준으로 아슬아슬하게 음수가 된다. 위 힐베르트 $n$ = 14가 정확히 후자다 — 13번째, 즉 마지막 직전에서 터졌다. ↩

5. 유한요소에서 강체 모드(rigid body mode)를 구속하지 않으면 강성행렬이 특이해져 촐레스키가 대각에서 0을 만나 멈춘다. “왜 안 돌아가지”의 8할은 경계 조건 누락이다. ↩