LDL 분해

1. 개요

LDL 분해 LDLT Factorization
대상	대칭 행렬 (양정치·부정치 모두)
형태	PAPT = LDLT, L 단위 하삼각, D 대각 또는 1×1·2×2 블록 대각
연산량	밀집 n³/3 flops (촐레스키와 같음)
덤	D 의 부호로 행렬의 관성(양·음·영 고유값 개수)을 공짜로 얻음
대표 구현	LAPACK xSYTRF · MA57 · MUMPS · PARDISO · QDLDL

촐레스키는 양정치 행렬만 받아 준다. 세상의 대칭 행렬 절반은 양정치가 아니다. 그 나머지 절반을 위해 존재하는 분해.

LDL 분해(LDLᵀ factorization)는 대칭 행렬 $A$ 를 단위 하삼각 행렬 $L$ 과 (블록) 대각 행렬 $D$ 의 곱 $A = LDL^{\mathsf T}$ 로 쓰는 분해다. 양정치 행렬이라면 $D$ 가 양의 대각이고 $L D^{1/2}$ 가 곧 촐레스키 분해의 인자라서, 제곱근만 뺀 촐레스키 변형에 불과하다. LDL 이 독자적인 의미를 갖는 건 대칭이지만 부정치(indefinite)인 행렬, 즉 고유값에 음수가 섞인 행렬에서다.

부정치 대칭 행렬은 시뮬레이션에서 흔하다. 비압축 유동이나 혼합유한요소법의 속도-압력 계, 구속조건이 붙은 최적화의 카루시-쿤-터커 조건 행렬, 이동(shift)을 건 고유값 문제의 $K - \sigma M$ 이 전부 그렇다. 이런 행렬에 LU 분해를 쓰면 대칭성을 버리느라 메모리와 연산이 두 배가 되고, 촐레스키는 아예 깨진다. LDL 은 대칭을 지키면서 부정치를 견디는 거의 유일한 직접법이다. 다만 그러려면 $D$ 에 2×2 블록을 허용해야 한다는 게 이 문서의 핵심이다.

2. 양정치일 때 — 제곱근 없는 촐레스키

$A$ 가 양정치면 피벗팅 없이 $j = 1,\dots,n$ 순서로

$$d_j = a_{jj} - \sum_{k<j} \ell_{jk}^2 d_k,\qquad \ell_{ij} = \frac{1}{d_j}\Big(a_{ij} - \sum_{k<j} \ell_{ik}\ell_{jk} d_k\Big)\quad (i>j)$$

를 계산하면 된다. 촐레스키의 교과서 예제

$$A = \begin{pmatrix} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \end{pmatrix}$$

에 넣으면 $D = \mathrm{diag}(4,\ 1,\ 9)$, $L$ 의 아래쪽 원소가 $\ell_{21}=3$, $\ell_{31}=-4$, $\ell_{32}=5$ 다. 촐레스키 인자 $\mathrm{diag}(2,1,3)$ 을 열마다 곱하면 정확히 촐레스키의 $L$ 이 된다. 정수 행렬이 정수로 떨어지는 건 제곱근이 없어서다.

이 형태가 쓸모 있는 곳은 스카이라인 저장 같은 열 단위 유한요소 솔버(능동 열 LDLᵀ)와, 대각이 0에 가까워지는 걸 부호 검사로 잡고 싶은 경우다. 촐레스키는 음수 제곱근에서 그냥 멈추지만, LDL 은 $d_j$ 의 부호를 보고 “몇 번째 피벗에서 양정치성이 깨졌는지”를 계속 셀 수 있다.

3. 부정치일 때 — 대각 피벗만으로는 안 된다

가장 작은 반례가 이것이다.

$$A = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$$

고유값이 $\pm1$ 이라 조건수는 1, 더할 나위 없이 착한 행렬이다. 그런데 $d_1 = a_{11} = 0$ 이라 첫 단계에서 0으로 나눈다. 대칭 치환 $PAP^{\mathsf T}$ 는 대각 원소를 대각 원소끼리만 섞으니, 대각이 전부 0인 이 행렬은 어떤 대칭 치환으로도 대각 $D$ 를 갖는 LDL 이 존재하지 않는다.

0이 아니어도 작으면 문제다. $A = \begin{pmatrix} \varepsilon & 1 \\ 1 & 1 \end{pmatrix}$ 을 피벗 없이 분해하면 $\ell_{21} = 1/\varepsilon$, $d_2 = 1 - 1/\varepsilon$ 이다. $\varepsilon = 10^{-10}$ 이면 원래 원소가 1 이하인 행렬에서 인자 원소가 $10^{10}$ 으로 자란다. 이 성장 인자가 크면 후진 오차 상계가 같이 커져서, 원래 행렬과 상관없는 반올림 오차가 해에 실린다. 가우스 소거법이 부분 피벗팅을 하는 이유와 똑같다.

해법은 $D$ 를 1×1 과 2×2 대각 블록의 혼합으로 허용하는 것이다. 위 $\begin{pmatrix} 0&1\\1&0 \end{pmatrix}$ 은 그 자체로 2×2 피벗 하나가 되고, $L = I$ 로 끝난다. 2×2 블록 $E$ 는 대칭이고 역행렬을 닫힌 꼴로 구할 수 있으니, 그 블록으로 아래쪽을 소거하는 건 1×1 피벗과 거의 같은 비용이다.

4. 피벗 전략

어느 순간에 1×1 을 쓰고 언제 2×2 를 쓸지, 어떤 행·열을 피벗 자리로 끌어올지가 알고리즘의 전부다.

• 번치-파렛(Bunch–Parlett, 1971). 남은 부분행렬 전체를 뒤져 가장 큰 원소를 찾는 완전 피벗팅. 안정성은 좋지만 비교 연산이 $O(n^3)$ 이라 소거 자체만큼 비싸다.
• 번치-카우프만 피벗(1977). 한 단계에 열 두 개만 본다. 현재 열의 대각 $|a_{11}|$ 과 그 열의 최대 비대각 $\lambda$, 그 최대 원소가 있는 행 $r$ 의 최대 비대각 $\sigma$ 를 비교해 1×1(원래 자리 또는 $r$ 자리)과 2×2 중 하나를 고른다. 판정 상수

$$\alpha = \frac{1+\sqrt{17}}{8} \approx 0.6404$$

는 원소 성장의 상계를 최소화하도록 고른 값이고,¹ 그때 성장은 단계당 최대 $1 + 1/\alpha \approx 2.57$ 배로 묶인다. 비교 비용이 $O(n^2)$ 이라 LAPACK xSYTRF 의 기본값이 됐다.

• 룩(rook) 피벗팅. 번치-카우프만은 $D$ 의 성장은 묶지만 $L$ 의 원소 크기는 묶지 못해서, 애슈크래프트·그라임스·루이스(1998)가 $\lVert L\rVert$ 이 터지는 예를 보였다. 체스의 룩처럼 행과 열을 번갈아 더듬어 “자기 행·열에서 최대인” 원소를 찾으면 $L$ 의 원소가 상수로 묶인다. LAPACK 에는 xSYTRF_ROOK 으로 들어가 있다.
• 아센(Aasen, 1971). 블록 대각 대신 삼중대각 $T$ 로 $PAP^{\mathsf T} = LTL^{\mathsf T}$ 를 만든다. 피벗이 부분 피벗팅처럼 단순하고 블록화가 잘 돼서 최근 다시 주목받았다(xSYTRF_AA).

희소 행렬에서는 사정이 더 복잡하다. 최소차수 순서화·중첩절단으로 정한 순서를 수치 피벗이 흐트러뜨리기 때문이다. 다중전방 솔버(MA57·MUMPS)는 전방행렬 안에서만 피벗을 고르고, 안 되면 그 피벗을 부모 노드로 미룬다(delayed pivot). PARDISO 는 미루는 대신 작은 피벗을 $\pm\epsilon\lVert A\rVert$ 크기(대칭 부정치 기본값 $\epsilon = 10^{-8}$)로 교란하고 나중에 반복 정련으로 보정한다. 어느 쪽이든 희소 직접 솔버에서 부정치 행렬이 양정치보다 몇 배 다루기 어려운 이유가 이 피벗 문제다.

5. 관성 — D 의 부호가 고유값을 센다

실베스터 관성 법칙에 따르면 합동 변환 $A \mapsto SAS^{\mathsf T}$ (가역 $S$)은 양·음·영 고유값의 개수를 보존한다. $A = LDL^{\mathsf T}$ 는 바로 그런 합동이니, $D$ 의 고유값 부호를 세면 $A$ 의 관성이 나온다(2×2 블록은 행렬식이 음수면 양·음 하나씩). 고유값을 하나도 계산하지 않고 얻는 정보다.

5.1. 검산 — 이동을 건 개수 세기

$A$ 를 $3\times3$ 이산 라플라시안 $\mathrm{tridiag}(-1,\,2,\,-1)$ 이라 하자. 고유값은 $2-\sqrt2 \approx 0.586$, $2$, $2+\sqrt2 \approx 3.414$ 다. $\sigma = 2.5$ 로 이동한 $A - 2.5I = \mathrm{tridiag}(-1,\,-0.5,\,-1)$ 을 분해하면

$$d_1 = -0.5,\qquad d_2 = -0.5 - \frac{1}{-0.5} = 1.5,\qquad d_3 = -0.5 - \frac{1}{1.5} \approx -1.167$$

로 음수가 둘이다. 실제로 2.5 아래 고유값은 0.586과 2, 두 개다. 이 셈이 대칭 삼중대각에서는 스텀 수열의 부호 변화 개수와 같은 것이고, 이분법으로 $k$ 번째 고유값을 가두는 근거가 된다.

실무에서는 이게 안전장치로 쓰인다.

• 고유값 누락 검사. 유한요소 모드 해석에서 란초스 알고리즘이 구간 $[\sigma_1, \sigma_2]$ 의 모드를 다 찾았는지, $K - \sigma_1 M$ 과 $K - \sigma_2 M$ 을 LDL 분해해 음수 피벗 개수 차이와 맞춰 본다. 구조해석 업계에서 “스텀 검사(Sturm check)“라 부르는 것이다. 시프트-인버트에 어차피 그 분해가 필요하니 검사는 공짜다.
• 최적화의 관성 보정. 내점법 솔버 IPOPT 는 뉴턴 스텝마다 KKT 행렬 $\begin{pmatrix} W & J^{\mathsf T} \\ J & 0 \end{pmatrix}$ 을 LDL 분해하고, 관성이 $(n, m, 0)$ — 양수 $n$ 개, 음수 $m$ 개 — 인지 확인한다. 아니면 헤시안 블록에 $\delta I$ 를 더해 다시 분해한다. 관성이 맞아야 그 스텝이 제약 위에서 하강 방향이기 때문이다. 그래서 IPOPT 의 선형 솔버(MA27·MA57·MUMPS 등)는 반드시 관성을 보고할 수 있어야 한다.²

6. 준정치 행렬 — 피벗 없이도 되는 특수한 부정치

부정치라고 다 2×2 피벗이 필요한 건 아니다. 블록 구조가

$$K = \begin{pmatrix} H & A^{\mathsf T} \\ A & -G \end{pmatrix},\qquad H \succ 0,\ G \succ 0$$

인 준정치(quasi-definite) 행렬은 반더베이(1995)가 보였듯 어떤 대칭 치환에 대해서도 대각 $D$ 를 갖는 LDLᵀ 가 존재한다. 그러니 순서는 순전히 채움을 줄이는 기준으로만 고르고(기호 분해를 한 번만 하면 된다), 수치 단계는 피벗 없이 민다.

정규화된 내점법이나 이차계획법이 일부러 오른쪽 아래에 $-\delta I$ 를 넣어 KKT 행렬을 준정치로 만드는 이유가 이것이다. 이차계획 솔버 OSQP 가 쓰는 QDLDL 은 이름 그대로 준정치 LDL 전용의 수백 줄짜리 구현이다. 대가는 원래 문제에서 $\delta$ 만큼 벗어난다는 것이고, 반복 정련 몇 번으로 그 차이를 메운다.³

7. 관련 문서

• 촐레스키 분해 · LU 분해 · 가우스 소거법 · 스카이라인 저장
• 희소 직접 솔버 · 최소차수 순서화 · 중첩절단 · MUMPS
• 번치-카우프만 피벗 · 실베스터 관성 법칙 · 스텀 수열 · 란초스 알고리즘
• 내점법 · 카루시-쿤-터커 조건 · 안장점 · 혼합유한요소법

8. Footnotes

1. LAPACK 소스를 열면 이 값이 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 로 박혀 있다. 17 이라는 상수를 위해 SEVTEN 이라는 변수까지 따로 만들어 둔 정성이 포트란답다. ↩

2. IPOPT 로그에 찍히는 lg(rg) 열이 바로 이 헤시안 보정량 $\delta$ 의 로그다. 그 열이 계속 차 있으면 문제가 비볼록하다는 뜻이고, 대개 사람이 모델링을 다시 봐야 한다. ↩

3. 정확한 문제를 못 풀면 살짝 다른 쉬운 문제를 풀고 차이를 메운다. 수치해석 전체를 한 줄로 요약하면 대략 이렇다. ↩