LU 분해

1. 개요

한 번 분해해두면, 우변만 바뀌는 문제는 공짜로 푼다. 게으른 자를 위한 선형대수의 축복.

LU 분해(LU decomposition)는 정사각행렬 $A$를 하삼각행렬(lower triangular) $L$과 상삼각행렬(upper triangular) $U$의 곱으로 인수분해하는 것이다.

$$A = LU$$

여기서 $L$은 대각 성분이 모두 1인 단위 하삼각행렬, $U$는 상삼각행렬이다. 본질적으로 가우스 소거법의 전진 소거 과정에서 나온 배수(multiplier)들을 버리지 않고 $L$에 저장한 것이 LU 분해다. 즉 둘은 같은 계산의 두 얼굴인데, LU 분해가 훨씬 실용적인 이유는 한 번 분해해두면 재활용이 가능하기 때문이다.

이 재활용성 덕에 LU 분해는 수치해석의 밥줄이 되었다. 구조해석에서 하중 조건(우변)만 바꿔가며 여러 케이스를 돌리거나, 반복법의 전처리기로 불완전 LU를 쓸 때 모두 이 분해가 주인공이다.

2. 왜 분해하는가

$A\mathbf{x} = \mathbf{b}$를 그냥 가우스 소거법으로 풀면 우변 $\mathbf{b}$가 하나 바뀔 때마다 $\tfrac{2}{3}n^3$ flops짜리 소거를 처음부터 다시 해야 한다. 미련한 짓이다.

반면 $A = LU$로 한 번 분해($\tfrac{2}{3}n^3$ flops)해두면, 이후 임의의 우변에 대해 다음 두 단계로 해를 구할 수 있다.

$$L\mathbf{y} = \mathbf{b} \quad(\text{전진대입}), \qquad U\mathbf{x} = \mathbf{y} \quad(\text{후진대입})$$

이 두 삼각 시스템 풀이는 각각 $2n^2$ flops, 합쳐 $4n^2$밖에 안 든다. $n = 1000$이고 우변이 100개인 상황을 실제로 세어 보면 차이가 확연하다.

방식	분해	우변 100개 풀이	합계
매번 소거	—	$100 \times 6.7 \times 10^8$	약 $6.7 \times 10^{10}$
한 번 LU	$6.7 \times 10^8$	$100 \times 4 \times 10^6$	약 $1.1 \times 10^9$

60배다. 우변이 많아질수록 비율은 $\tfrac{n}{6}$에 가까워진다. 이게 LU 분해가 단순 소거법을 이긴 결정적 한 방이고, 하중 케이스가 수십~수백 개인 구조해석에서 직접법이 여전히 사랑받는 이유이기도 하다.

3. 전진대입과 후진대입

삼각행렬을 푸는 것은 놀랍도록 간단하다. 전진대입(forward substitution)은 $L\mathbf{y} = \mathbf{b}$에서 위에서부터 차례로 풀어 내려간다.

$$y_i = b_i - \sum_{j=1}^{i-1} \ell_{ij}\, y_j, \qquad i = 1, 2, \dots, n$$

($L$의 대각이 1이라 나눗셈이 없다.) 이렇게 얻은 $\mathbf{y}$로 후진대입(back substitution)을 수행한다.

$$x_i = \frac{1}{u_{ii}}\left( y_i - \sum_{j=i+1}^{n} u_{ij}\, x_j \right), \qquad i = n, n-1, \dots, 1$$

두 대입 모두 미지수를 하나씩 순차적으로 확정하므로 재귀적 의존성이 있다. 이 때문에 삼각 풀이는 병렬 컴퓨팅에서 병렬화가 까다로운 대표적 골칫거리로 꼽힌다.¹ 게다가 연산량 $4n^2$에 대해 읽어야 할 데이터도 $O(n^2)$이라 산술강도(arithmetic intensity)가 $O(1)$이다 — 즉 삼각 풀이는 메모리 대역폭에 묶인(memory-bound) 연산이고, 분해가 연산에 묶인(compute-bound) 연산인 것과 대조된다.

4. 피벗팅과 존재성

모든 행렬이 곱게 LU로 쪼개지는 건 아니다. 소거 도중 대각 피벗이 0이 되면 그 자리에서 알고리즘이 멈춘다. 그래서 실무에서는 가우스 소거법과 마찬가지로 부분 피벗팅을 곁들여 행을 교환하며, 그 결과는 치환행렬 $P$를 동원한 형태가 된다. 피벗팅 없는 순수 $A=LU$는 대각 우세(diagonally dominant)하거나 대칭 양의 정부호인 얌전한 행렬에만 안전하게 통한다.

여기서 규약이 라이브러리마다 다르다는 점이 초보를 잡는다. $P$가 치환행렬이라 $P^{-1} = P^{\mathsf{T}}$이므로 다음 두 표현은 수학적으로 같은데, 함수가 돌려주는 $P$는 서로 전치 관계다.

규약	쓰는 곳
$A = P L U$	LAPACK dgetrf 문서, SciPy scipy.linalg.lu
$P A = L U$	MATLAB [L,U,P] = lu(A), 대부분의 교과서

받은 $P$를 어느 쪽에 곱해야 하는지 헷갈리면 $\|PLU - A\|$와 $\|LU - PA\|$를 둘 다 재 보면 끝난다. 0에 가까운 쪽이 그 라이브러리의 규약이다.² LAPACK은 실제로는 행렬 $P$가 아니라 교환 인덱스 배열 ipiv를 돌려주고, dlaswp로 그걸 적용한다 — 치환행렬을 실제로 만들어 곱하는 건 $O(n^3)$ 낭비라서 아무도 하지 않는다.

5. 블록 알고리즘과 BLAS 3

교과서의 스칼라 삼중 루프를 그대로 구현하면 현대 CPU에서 이론 성능의 몇 퍼센트밖에 못 낸다. 갱신 한 번에 값 하나만 쓰고 버려서 캐시를 못 쓰기 때문이다. 그래서 LAPACK의 dgetrf는 블록 right-looking 알고리즘을 쓴다. 열을 블록 폭 $n_b$(보통 32~256)만큼 묶어 그 좁은 패널만 비블록 루틴으로 분해하고, 남은 부분행렬은 dtrsm(삼각 풀이)과 dgemm(행렬 곱)으로 한꺼번에 갱신한다.

핵심은 전체 $\tfrac{2}{3}n^3$ flops의 대부분이 dgemm 한 곳으로 몰린다는 것이다. 행렬-행렬 곱은 데이터 $O(n^2)$당 연산 $O(n^3)$이라 캐시 블로킹이 잘 먹히는 유일한 연산이고, 그래서 블록 LU는 피크 성능의 상당 부분을 실제로 뽑아낸다. “LU 분해를 직접 구현하지 마라”는 조언은 알고리즘이 어려워서가 아니라, 이 메모리 계층 튜닝을 혼자 이길 수 없기 때문이다.

6. 촐레스키 분해

행렬이 대칭 양의 정부호(symmetric positive definite, SPD)라는 특권을 누린다면, 촐레스키 분해 $A = LL^{\mathsf{T}}$를 쓸 수 있다. 저장 공간이 절반, 연산량도 절반($\tfrac{1}{3}n^3$ flops)으로 줄고, 피벗팅 없이도 무조건 안정해서 코드가 깔끔하다. SPD 행렬은 구조해석의 강성행렬이나 유한요소법 이산화에서 흔하게 튀어나오므로, 촐레스키는 이 바닥에서 사랑받는 국룰 분해다. 제곱근 계산이 부담스럽거나 대칭이지만 부정부호(indefinite)인 행렬까지 다뤄야 하면 $A = LDL^{\mathsf{T}}$ 변형을 쓴다.³

7. 희소성과 채움현상

LU 분해가 희소행렬을 만나면 이야기가 복잡해진다. 원래 0이던 자리가 소거 과정에서 0이 아닌 값으로 채워지는 채움현상(fill-in)이 발생하기 때문이다.⁴ 이걸 방치하면 $L$과 $U$가 원래 $A$보다 훨씬 빽빽해져서 메모리와 연산이 폭발한다.

그래서 희소 직접 솔버는 분해 전에 행과 열을 재배열(reordering)해 채움현상을 최소화한다. AMD(approximate minimum degree)나 중첩절단(nested dissection) 같은 순서화 알고리즘이 여기 동원되고, PARDISO·MUMPS·UMFPACK 같은 라이브러리가 이 마법을 대신 부려준다. 다만 순서화로 이길 수 있는 한계가 차원에 따라 갈린다 — 2차원 정규격자는 중첩절단이 $O(N \log N)$ 메모리로 막아 주지만, 3차원에서는 분리자가 면이라 $O(N^{4/3})$ 메모리·$O(N^2)$ 연산이 하한이다. 그래서 3차원 대규모 문제는 결국 반복법과 다중격자법에 자리를 내준다. 자세한 비교는 희소 직접 솔버 문서에 있다.

8. 관련 문서

• 가우스 소거법
• 촐레스키 분해
• 희소행렬 · 희소 직접 솔버
• LAPACK
• 반복법
• 전처리기
• 조건수
• QR 분해
• 강성행렬
• 고유값 문제
• 수치해석
• 크리로프 부분공간법

9. Footnotes

1. 순차 의존성 때문에 삼각 풀이(triangular solve)는 GPU에서 특히 애를 먹인다. 레벨 스케줄링(level scheduling)으로 병렬 가능한 미지수들을 묶어 처리하지만, 밀집 행렬 곱셈만큼 시원하게 뻥 뚫리진 않는다. ↩

2. 이 규약 차이는 MATLAB 코드를 NumPy로 옮기다 “답이 행만 섞여서 나온다”는 형태로 발견되는 게 보통이다. 부호도 값도 다 맞는데 순서만 틀린 답을 보면 대체로 $P$ 방향을 뒤집은 것이다. ↩

3. 앙드레루이 촐레스키(André-Louis Cholesky, 1875~1918)는 프랑스 육군 포병 장교이면서 측지·지형 부서에서 삼각측량 계산을 담당했다. 최소제곱 정규방정식을 손으로 빨리 풀려고 만든 방법인데, 정작 본인은 1918년 8월 31일 서부전선에서 입은 부상으로 전사했고, 방법은 동료 장교 브누아(Benoît)가 1924년 Bulletin Géodésique에 발표해 세상에 알려졌다. 측량하다 만든 게 지금은 분자동역학까지 굴린다. ↩

4. 채움현상이 얼마나 얄미운가 하면, 아주 성긴 행렬도 순서를 잘못 잡으면 분해 결과가 거의 밀집 행렬이 되어버린다. “0을 잘 정렬하는 것”이 곧 성능이라는, 직관에 반하는 세계다. ↩