희소 직접 솔버

수치해석 시뮬레이션 마지막 수정: 2026-09-01 04:51:44

1. 개요[편집]

희소 직접 솔버 Sparse Direct Solver
하는 일	$PAQ = LU$ 또는 $PAP^{\mathsf T} = LL^{\mathsf T}$ 를 명시적으로 만든다
최대 적	채움현상(fill-in) — 0이던 자리가 채워진다
승부처	순서화. 나머지는 다 구현 디테일
3단계	심볼릭 분해 → 수치 분해 → 전진·후진 대입
2D 정규격자	중첩절단으로 $\mathrm{nnz}(L)=O(N\log N)$, 연산 $O(N^{1.5})$ — 최적
3D 정규격자	$\mathrm{nnz}(L)=O(N^{4/3})$, 연산 $O(N^{2})$ — 여기서 죽는다
대표 구현	MUMPS · UMFPACK · CHOLMOD · PARDISO · SuperLU

반복법은 “수렴은 신에게 맡긴다”이고, 직접법은 “메모리는 신에게 맡긴다”이다. 둘 중 뭘 맡길지가 3차원에서 갈린다.

희소 직접 솔버는 희소행렬 $A$ 를 유한 번의 소거로 삼각인수 $LU$ 또는 $LL^{\mathsf T}$ 로 분해한 뒤, 전진·후진 대입으로 $Ax=b$ 를 정확히(반올림 오차 범위에서) 푸는 알고리즘군이다. 반복법과 달리 수렴 여부를 걱정할 필요가 없고, 초기 추정값도 전처리기 튜닝도 필요 없다. 한 번 분해해 두면 우변이 몇 개든 대입만으로 풀린다. MATLAB에서 무심코 치는 A\b 가 희소행렬을 만나면 뒤에서 돌아가는 물건이 정확히 이것이다.

대가는 하나뿐인데, 그게 치명적이다. 소거를 하면 원래 0이던 자리에 0이 아닌 값이 새로 생긴다. 이 채움현상(fill-in) 때문에 성기던 행렬이 분해 후에는 거의 밀집이 되고, 애써 아낀 메모리가 도로 폭발한다. 이 문서의 8할은 채움을 어떻게 줄이느냐의 이야기다. LU 분해·촐레스키 분해 자체의 수식은 각 문서에 있으니, 여기서는 희소성이 붙었을 때 무엇이 달라지는가만 본다.

2. 채움현상과 소거 그래프[편집]

대칭 양정치 행렬 $A$ 의 희소 패턴을 그래프 $G(A)$ 로 보자 — 정점은 미지수, 간선은 $a_{ij} \ne 0$ . 촐레스키 소거의 $k$ 단계에서 $x_k$ 를 소거할 때 일어나는 갱신은

a_{ij} \;\leftarrow\; a_{ij} - \frac{a_{ik}a_{kj}}{a_{kk}}

인데, $a_{ik}$ 와 $a_{kj}$ 가 모두 0이 아니면 $a_{ij}$ 는 원래 0이었더라도 0이 아니게 된다. 그래프 언어로 옮기면 규칙이 딱 한 줄이다.

정점 $k$ 를 제거하고, $k$ 의 이웃들을 서로 전부 연결해 클리크로 만든다.

이것이 소거 게임(elimination game)이고, 새로 생긴 간선이 곧 채움이다. 여기서 두 가지가 바로 따라 나온다.

차수가 큰 정점을 먼저 소거하면 재앙이다. 이웃이 $d$ 개면 $\binom{d}{2}$ 개의 간선이 한 번에 생긴다.
채움이 0이 되는 순서가 존재할 조건은 그래프가 현(chordal) 그래프인 것이다. 길이 4 이상의 사이클에 지름길 간선이 항상 있으면 완전 소거 순서가 존재한다. 삼중대각 행렬이나 트리 구조 행렬이 이 경우이고, 그래서 토머스 알고리즘은 채움이 아예 없다.

문제는 일반 행렬이다. 채움을 최소화하는 순서를 찾는 문제는 NP-완전이다(야나카키스, 1981). 그래서 실무는 전부 휴리스틱이고, 놀랍게도 그 휴리스틱들이 대단히 잘 작동한다.

20×20 격자 5점 라플라시안(N=400)을 소거 순서만 바꿔 심볼릭 촐레스키로 분해한다. 빨간 점이 fill-in — 자연순서 nnz(L)=8019(fill 6859), 중첩절단 5062, 최소차수 3729로 같은 행렬의 flops가 8.63e4 대 2.62e4로 갈린다. 자연순서 L은 밴드 촐레스키로 실제로 풀어 ‖Ax−b‖∞=4.0e−15 까지 확인한 자리다.

3. 순서화가 전부다[편집]

같은 행렬을 어떤 순서로 소거하느냐에 따라 메모리와 연산량이 자릿수 단위로 달라진다. 2차원 $k\times k$ 정규격자( $N = k^2$ 미지수, 5점 스텐실)를 벤치마크로 삼아 비교하면 이렇다.

순서화	발상	2D 격자 $\mathrm{nnz}(L)$	2D 격자 연산량
자연순서 (행 우선)	아무 생각 없음	$O(N^{3/2})$	$O(N^{2})$
역 커스힐-맥키 (RCM)	대역폭·프로파일 축소	$O(N^{3/2})$	$O(N^{2})$
최소차수 (MD / AMD)	매 단계 차수 최소 정점 선택	실측상 ND에 근접	실측상 ND에 근접
중첩절단 (ND)	분리자로 재귀 이분	$O(N \log N)$	$O(N^{3/2})$

자연순서는 대역폭이 $k=\sqrt N$ 인 띠행렬을 만든다. 띠 안은 전부 채워지므로 $\mathrm{nnz}(L) \approx N\sqrt N$ . 100만 미지수면 $L$ 만 10억 개 — 이미 끝났다.

역 커스힐-맥키(Cuthill–McKee 1969, 역순 취하기는 조지 1971)는 너비 우선 탐색으로 레벨 집합을 만들어 번호를 매긴 뒤 뒤집는다. 대역폭·프로파일을 크게 줄여서 띠 솔버에는 큰 이득이고, 캐시 지역성이 좋아져 반복법의 행렬-벡터 곱도 빨라진다. 다만 채움 최소화 관점에서는 AMD·ND에 밀린다.

최소차수는 마코비츠(1957)의 비대칭판과 티니-워커(1967)의 대칭판이 원조다. 매 단계 소거 그래프에서 차수가 가장 작은 정점을 고른다 — 완벽한 그리디 알고리즘이다. 문제는 소거할 때마다 그래프가 바뀌어 차수 갱신 비용이 크다는 것이고, 이를 근사 최소차수(AMD; 아메스토이·데이비스·더프 1996)가 차수 대신 상한을 쓰고 원소(element) 개념으로 클리크를 압축 표현해 해결했다. 오늘날 중소 규모 문제의 기본값이며, 이론적 보장은 거의 없는데 실측 성능이 대단히 좋은 전형적인 “일단 돌려” 알고리즘이다.¹

중첩절단(nested dissection)은 조지(1973)가 정규 유한요소 격자에 대해 제안했다.² 그래프를 대략 균등한 두 조각 $\Omega_1, \Omega_2$ 와 이들을 갈라놓는 분리자 $S$ 로 나눈 뒤, $\Omega_1 \to \Omega_2 \to S$ 순으로 번호를 매기고 각 조각에 재귀한다. $\Omega_1$ 을 소거할 때 생기는 채움은 $\Omega_1 \cup S$ 안에 갇히고, $\Omega_2$ 의 정점과는 영원히 연결되지 않는다 — 사이에 $S$ 가 있어 경로가 없기 때문이다. 이 “차단” 효과가 채움을 구조적으로 봉쇄한다.

2D $k \times k$ 격자면 최상위 분리자가 정점 $k = \sqrt N$ 개짜리 십자선이고, 재귀적으로 분리자 크기가 절반씩 줄어든다. 계산해 보면

\mathrm{nnz}(L) = \Theta(N \log N), \qquad \text{연산량} = \Theta(N^{3/2})

이 나온다. 그리고 호프만·마틴·로즈(1973)가 같은 격자에 대해 어떤 순서를 써도 이보다 잘할 수 없다는 하한을 증명했다. 즉 중첩절단은 상수배 안에서 최적이다. 립턴·로즈·타잔(1979)이 이를 일반화해, $O(\sqrt n)$ 크기 분리자를 갖는 그래프(평면 그래프 전부 포함)에 같은 결과가 성립함을 보였다.

3차원으로 가면 이야기가 달라진다. $k\times k\times k$ 격자의 분리자는 면이라 $k^2 = N^{2/3}$ 개 정점이고, 최상위 분리자 하나를 소거하는 것만으로 밀집 $N^{2/3} \times N^{2/3}$ 행렬 분해가 필요하다.

\mathrm{nnz}(L) = \Theta(N^{4/3}), \qquad \text{연산량} = \Theta(N^{2})

$N = 10^6$ (격자 $100^3$ )이면 $L$ 만 대략 $10^8$ 개 원소, 배정밀도로 800 MB에 인덱스까지 얹힌다. $N = 10^7$ 이면 연산량이 $10^{14}$ — 이 지점에서 3D 직접법은 사실상 사망 선고를 받는다. 직접법이 3D에서 죽는 이유는 시간이 아니라 메모리다.³

실무 순서화는 하이브리드다. 그래프 상위 몇 단계는 METIS·SCOTCH 같은 그래프 분할 도구로 중첩절단을 하고(피들러 벡터 기반 스펙트럴 이분이나 다단계 KL 정련이 쓰인다), 조각이 충분히 작아지면 AMD로 마무리한다. 이 조합이 거의 모든 현대 솔버의 기본 설정이다. 비대칭 $LU$ 는 순서화를 $A + A^{\mathsf T}$ 나 $A^{\mathsf T}A$ 의 패턴에 적용하며(COLAMD 계열), 이 미묘한 차이가 성능을 꽤 흔든다.

4. 세 단계[편집]

희소 직접 솔버는 반드시 세 단계로 쪼개져 있고, 이 분리가 성능의 핵심이다.

1. 심볼릭 분해(해석 단계). 수치값은 전혀 보지 않고 패턴만으로 $L$ 의 0 아닌 자리를 전부 예측한다. 순서화도 여기 포함된다. 촐레스키에서는 이게 정확히 맞아떨어지는데, 대칭 양정치 행렬은 피벗팅이 필요 없어서 소거 순서가 값과 무관하게 미리 확정되기 때문이다. 이 단계가 산출하는 것이 소거 트리(elimination tree)로, 정의는 간단하다.

\mathrm{parent}(j) \;=\; \min \{\, i > j \;:\; \ell_{ij} \ne 0 \,\}

소거 트리는 열 $j$ 의 갱신이 어느 열로 전파되는지를 말해 주고, 따라서 형제 서브트리는 서로 독립이라 병렬로 처리할 수 있다. 병렬 컴퓨팅 스케줄링이 통째로 이 트리 위에서 이뤄진다. 열별 nonzero 개수도 소거 트리로 거의 선형 시간에 셀 수 있어서, 심볼릭 단계는 대개 수치 단계보다 훨씬 싸다.

2. 수치 분해. 예측된 자리에 실제 값을 채운다. 전체 시간의 대부분이 여기 있고, 뒤에 나올 초노드·다중전방 기법이 전부 이 단계를 BLAS-3로 돌리기 위한 장치다.

3. 전진·후진 대입. $Ly = b$ , $L^{\mathsf T}x = y$ . 비용이 $O(\mathrm{nnz}(L))$ 로 분해에 비해 압도적으로 싸다. 이 비대칭이 직접법의 존재 이유다 — 우변이 100개면 분해 한 번에 대입 100번이고, 반복법으로 100번 푸는 것보다 훨씬 빠르다.

행렬이 그대로이고 값만 바뀌는 상황(예: 뉴턴 반복에서 패턴이 같은 자코비안, 고정 시간간격 음해법)에서는 심볼릭 결과를 재사용하고 수치 분해만 다시 돌린다. 상용 솔버 API가 analyze / factorize / solve 로 나뉘어 있는 이유가 이것이다.

5. 초노드와 다중전방법[편집]

1980년대까지의 희소 솔버는 열 단위로 인덱스를 따라다니며 값을 갱신했고, 그래서 간접 참조 지옥이었다. 부동소수점 유닛은 놀고 메모리 지연만 먹는 전형적인 상황. 돌파구는 “희소 안에서 밀집 덩어리를 찾아내자”였다.

초노드(supernode)는 대각 아래 희소 패턴이 동일한 연속된 열들의 묶음이다. 소거 트리에서 자식이 하나뿐인 경로가 대개 초노드가 된다. 한 초노드는 밀집 사다리꼴 블록이므로, 갱신을 dgemm 한 번으로 처리할 수 있다. BLAS-1 수준이던 연산이 BLAS-3로 올라가면서 캐시 재사용률이 뛰고, 실측 성능이 수 배에서 열 배까지 개선된다. 행렬 곱셈이 왜 그렇게 최적화되어 있는지가 여기서 보답받는다.

다중전방법(multifrontal method, 더프·리드 1983)은 한 발 더 나간다. 소거 트리의 각 노드마다 작은 밀집 전방행렬(frontal matrix)을 만들고, 거기서 몇 개 변수를 소거한 뒤 남은 슈어 보수를 갱신행렬로 부모에게 넘긴다. 부모는 자기 원소와 자식들의 갱신행렬을 인덱스 맞춰 더해(extend-add) 자기 전방행렬을 만든다. 결과적으로

모든 실수 연산이 밀집 행렬 커널 안에서 일어나고,
트리의 형제 노드는 완전히 독립이라 자연스러운 병렬 분해 단위가 되며,
활성 데이터가 스택 형태로 관리되어 메모리 국소성이 좋다.

MUMPS의 이름 자체가 MUltifrontal Massively Parallel Solver다. 트리 아래쪽의 잔가지들은 노드별 병렬(트리 병렬)로, 위쪽의 거대한 분리자 노드는 노드 안 밀집 분해를 여러 프로세스로 쪼개(노드 병렬) 처리하는 2단 전략이 표준이다.

최근에는 전방행렬의 비대각 블록이 수치적으로 저랭크라는 성질을 이용해 $O(N^2)$ 를 깎는 방향이 활발하다. MUMPS의 BLR(block low-rank), STRUMPACK의 HSS 등이 그것으로, 저랭크 근사를 허용 오차 안에서 도입해 3D의 메모리 벽을 미루는 시도다. 정확한 분해를 포기하고 근사 분해 + 반복 정련으로 가는 셈이라, 직접법과 반복법의 경계가 흐려지는 지점이기도 하다.

6. 피벗팅과 안정성[편집]

여기서 대칭과 비대칭의 운명이 갈린다.

대칭 양정치면 촐레스키는 어떤 순서로 소거해도 후진 안정이다. 그래서 순서를 순전히 채움 최소화 기준으로 고를 수 있고, 심볼릭 분해가 정확히 예측된다. 희소 솔버 세계에서 SPD 문제가 편애받는 이유.

비대칭이면 얘기가 완전히 달라진다. 안정성을 위해 부분 피벗팅을 하면 행 교환이 일어나고, 행 교환은 희소 패턴을 바꿔 심볼릭 예측을 무효화한다. 채움 최소화와 수치 안정성이 정면충돌하는 것이다. 처방은 셋이다.

임계값 부분 피벗팅. 열 최대값의 $\tau$ 배 이상이기만 하면 대각 원소를 그대로 쓴다( $\tau = 0.1$ 정도가 국룰). 안정성을 조금 양보하고 채움을 크게 아낀다. UMFPACK 계열의 기본 전략.
정적 피벗팅. 아예 피벗을 미리 고정한다. 분해 도중 피벗이 너무 작으면 $\varepsilon\|A\|$ 정도로 살짝 흔들어 진행하고, 그렇게 생긴 오차는 나중에 반복 정련(iterative refinement)으로 회수한다. SuperLU_DIST가 대규모 분산 환경에서 통신을 줄이려고 택한 길이다.
MC64 스케일링·순열. 큰 원소가 대각에 오도록 최대 가중치 이분 매칭을 푸는 전처리(더프·코스터 1999). 열마다 스케일을 곱해 대각을 1, 나머지를 1 이하로 만든 뒤 순열을 적용하면 정적 피벗팅이 훨씬 잘 먹는다. 이분 매칭이 수치선형대수 한복판에 등장하는 재미있는 사례.

대칭이지만 부정치인 경우(안장점 문제, 혼합유한요소법의 KKT 행렬 등)는 $LDL^{\mathsf T}$ 에 $1\times1$ / $2\times2$ 블록 피벗(번치-카우프만)을 섞어 대칭을 유지한 채 안정성을 확보한다.

무엇을 하든 마지막에는 잔차 $\|b - A\hat x\|$ 를 확인해야 한다. 직접법은 “정확한 해”를 준다는 인상을 주지만 실제로 보장하는 것은 후진 안정성이고, 전진 오차는 조건수에 비례해 부푼다. 자세한 것은 후진 오차 해석 문서 참조.

7. 구현체와 손익분기[편집]

실무에서 쓰는 물건들은 대략 이렇게 나뉜다.

이름	특기	비고
CHOLMOD	대칭 양정치, 초노드 촐레스키	SuiteSparse. MATLAB `chol` 의 엔진
UMFPACK	비대칭 다중전방 LU	SuiteSparse. MATLAB 백슬래시의 비대칭 경로
MUMPS	다중전방, MPI 분산, 대칭·비대칭 모두	학계·오픈소스 진영의 사실상 표준
PARDISO	초노드 LU, 공유메모리 병렬	Intel MKL에 포함된 버전이 널리 쓰임
SuperLU / SuperLU_DIST	초노드 LU, 정적 피벗팅 분산판	대규모 병렬 쪽
SPQR	희소 QR (다중전방)	최소자승 문제용

이제 본론. 언제 직접법을 쓰고 언제 크리로프 부분공간법이나 다중격자법으로 가야 하는가?

직접법이 이기는 판.

2차원 문제. 수백만 미지수까지도 중첩절단이 $O(N\log N)$ 으로 막아 준다. 사실상 무적.
우변이 많은 경우. 주파수 스윕, 다중 하중 케이스, 감도 해석의 수반 방정식. 분해 한 번을 수백 번 재사용한다.
이동-반전 고유값 해석. ARPACK·SLEPc의 shift-and-invert 모드는 매 반복마다 $(A-\sigma B)^{-1}$ 을 곱해야 하고, 이건 사실상 희소 직접 솔버 호출이다. 모드 해석이 직접 솔버에 목매는 이유.
반복법이 손드는 행렬. 강한 이방성·불연속 계수, 안장점 구조, 고주파 헬름홀츠 방정식, 미분대수방정식에서 나오는 극도로 나쁜 조건수. 전처리기를 아무리 만져도 GMRES가 정체하는 상황에서 직접법은 그냥 푼다.
회로·네트워크 문제. SPICE 계열은 행렬이 격자가 아니라 매우 성기고 불규칙해서 반복법의 전처리 이론이 잘 안 먹히고, 전통적으로 KLU 같은 전용 직접 솔버를 쓴다.

반복법이 이기는 판.

3차원 대규모 타원형 문제. 위에서 본 $O(N^{4/3})$ 메모리· $O(N^2)$ 연산이 벽이다. 대수적 다중격자를 전처리기로 쓴 CG는 좋은 조건에서 거의 $O(N)$ 이고 메모리도 $O(N)$ 이다. 수백만 이상 미지수의 3D 확산·구조 문제는 논쟁의 여지가 없다.
행렬을 명시적으로 만들 수 없을 때. 행렬-벡터 곱만 정의된 무행렬(matrix-free) 방식이면 직접법은 애초에 선택지가 아니다.
정확도가 덜 필요할 때. 시간 전진 한 스텝의 중간 해를 $10^{-6}$ 까지 풀 이유가 없다면, 반복법은 필요한 만큼만 풀고 멈출 수 있다. 직접법에는 그런 손잡이가 없다.

현실적인 절충안이 널리 쓰인다. 불완전 LU/촐레스키(ILU/IC)는 채움을 임계값이나 레벨로 잘라 버린 “고의로 부정확한” 직접법이고, 그 결과를 전처리기로 크리로프 법에 먹인다. 영역 분할법은 부분영역마다 직접 솔버를 돌리고 계면만 반복으로 처리한다. 즉 현대의 대규모 솔버는 대부분 직접법과 반복법의 혼합물이며, “직접이냐 반복이냐”는 이미 20년 전에 끝난 논쟁이다. 남은 질문은 어느 계층에서 자를 것인가뿐이다.

8. 실무 체크리스트[편집]

직접 솔버가 느리거나 죽었을 때 순서대로 볼 것들.

순서화를 확인했나. 라이브러리 기본값이 AMD인 경우가 많은데, 3D 큰 문제면 METIS 기반 중첩절단으로 바꾸는 것만으로 메모리가 절반 이하로 떨어지는 일이 흔하다. 대부분의 솔버가 순서화 옵션을 노출한다.
심볼릭을 재사용하고 있나. 시간 전진이나 뉴턴 반복에서 패턴이 안 바뀌는데 매번 analyze 를 다시 부르고 있다면, 그것만으로 20~40%를 버리는 중이다.
분해 후 메모리를 예측했나. 대부분의 솔버는 해석 단계가 끝나면 예상 $\mathrm{nnz}(L)$ 과 필요 메모리를 보고한다. 이걸 안 보고 잡을 던지는 것이 3D OOM의 표준 경로다.
행렬이 정말 비대칭인가. 패턴만 대칭이고 값이 비대칭이면 symmetric pattern 옵션으로 순서화 비용과 저장을 크게 줄일 수 있다. 대칭 양정치인데 비대칭 LU로 풀고 있으면 연산량이 두 배다.
스케일링을 했나. 물성치 단위가 뒤섞여 행마다 크기가 $10^{12}$ 배 차이 나는 행렬(전자기-열-구조 연성이 대표적)은 피벗 선택을 망가뜨린다. 행·열 균형 스케일링이나 MC64를 먼저 걸어야 한다.
잔차를 봤나. 분해가 성공했다고 답이 맞는 게 아니다. $\|b-A\hat x\| / \|b\|$ 를 찍고, 크면 반복 정련 한두 번을 돌린다. 대개 그걸로 해결되고, 안 되면 행렬 자체가 특이에 가까운 것이다.
정말 직접법이어야 하나. 우변이 하나뿐이고 3D 타원형이고 미지수가 수백만이면, 답은 거의 항상 “아니오”다.

9. 관련 문서[편집]

희소행렬 · LU 분해 · 촐레스키 분해 · 가우스 소거법
중첩 분할 · 최소 차수 순서화 · 그래프 분할 · METIS · 피들러 벡터
반복법 · 크리로프 부분공간법 · 켤레기울기법 · 다중격자법 · 전처리기
슈어 보수 · 저랭크 근사 · 영역 분할법 · 행렬 곱셈
조건수 · 후진 오차 해석 · 부동소수점 연산
MUMPS · PETSc · LAPACK · SPICE
깊이 우선 탐색 · 너비 우선 탐색 · 이분 매칭 · 강한 연결 요소

AMD가 대단히 잘 작동하는 이유에 대한 만족스러운 이론은 아직 없다. 최악의 경우를 만드는 반례 그래프는 얼마든지 만들 수 있는데, 실제 공학 문제의 행렬에서는 거의 항상 중첩절단에 근접한 결과를 낸다. “왜 되는지는 모르겠는데 잘 된다”는 점에서 난류 모델링의 $k$ - $\varepsilon$ 과 정신적으로 형제 사이다. ↩
George, A. (1973). “Nested dissection of a regular finite element mesh”. SIAM J. Numer. Anal. 10(2), 345–363. 앨런 조지는 이 논문 하나로 희소행렬 순서화 이론의 방향을 통째로 바꿨다. 재미있는 것은 발상이 순전히 기하학적이라는 점 — “격자를 십자로 자르면 양쪽이 서로 못 본다”는 그림 한 장이 이론의 전부이고, 대수는 그 그림을 따라간다. ↩
3D 직접법의 메모리 폭발은 실무에서 대개 “잘 돌던 해석이 격자를 두 배 촘촘히 했더니 스왑을 치기 시작했다”는 형태로 찾아온다. 격자 간격 절반이면 $N$ 이 8배, $\mathrm{nnz}(L)$ 은 $8^{4/3} = 16$ 배, 연산량은 64배다. 반복법이었다면 8배와 8배로 끝났을 일이다. 이 계산을 미리 해 두지 않으면 금요일 저녁에 잡을 던지고 월요일에 OOM 로그를 보게 된다. ↩