크리로프 부분공간법

1. 개요

크리로프 부분공간법 Krylov Subspace Methods
분류	대규모 선형계 반복 해법
대표 알고리즘	CG, GMRES, BiCGSTAB, MINRES
핵심 도구	부분공간, 직교화, 아르놀디/란초스
단짝	전처리기

역행렬은 만드는 게 아니다. 곱하는 것이다.

크리로프 부분공간법(Krylov Subspace Methods)은 커다란 선형계 $A\mathbf{x} = \mathbf{b}$를, 행렬 $A$를 벡터에 곱하는 연산만으로 반복적으로 풀어내는 수치 해법의 총칭이다. 이름은 러시아의 해군 조선공학자이자 응용수학자인 알렉세이 크릴로프(Alexei N. Krylov, 1863~1945)가 1931년 발표한, 특성다항식 계산에 이 벡터열을 쓴 논문에서 왔다. 정작 그는 오늘날의 반복 해법을 만든 적이 없다.¹ 크리로프 반복은 20세기 10대 알고리즘 목록에 이름을 올린 몇 안 되는 선형대수 기법이기도 하다.

핵심 아이디어는 이렇다. $A^{-1}$을 직접 구하는 건 미친 짓이지만², $A$를 벡터에 한 번 곱하는 건 싸다. 특히 유한체적법이나 유한요소법에서 튀어나오는 행렬은 대부분 희소행렬(sparse matrix)이라, 0이 아닌 성분만 곱하면 연산량이 격자 크기에 거의 선형이다. 그래서 “곱셈만으로 승부 본다”는 이 접근이 대규모 시뮬레이션의 사실상 표준 솔버가 되었다.

2. 크리로프 부분공간이란

시작점은 잔차 $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$ 하나다. 여기에 $A$를 계속 곱해서 만든 벡터들이 펼치는 공간이 바로 크리로프 부분공간이다.

$$\mathcal{K}_m(A, \mathbf{r}_0) = \mathrm{span}\{\mathbf{r}_0,\ A\mathbf{r}_0,\ A^2\mathbf{r}_0,\ \dots,\ A^{m-1}\mathbf{r}_0\}$$

이 $m$차원 공간 안에서 “가장 그럴듯한 근사해”를 골라내는 게 게임의 전부다. 왜 하필 이 공간이냐면, 케일리-해밀턴 정리에 따라 정칙 행렬 $A$의 역행렬이 $A$의 거듭제곱들의 선형결합으로 표현되고, 따라서 $A^{-1}\mathbf{b}$가 $A$에 대한 어떤 다항식 $p(A)\mathbf{b}$와 정확히 같아지기 때문이다. 크리로프 부분공간은 바로 그 다항식들이 사는 동네다.

더 날카롭게 말하면 종료 시점을 결정하는 것은 행렬 크기 $n$이 아니라 $\mathbf{r}_0$에 대한 $A$의 최소다항식 차수 $d$다. 크리로프 부분공간의 차원은 $d$에서 성장을 멈추므로(더 곱해 봐야 새 방향이 안 나온다) 정확해는 늦어도 $d$번째 반복에서 나온다. 고유값이 $d$종류밖에 없는 행렬이면 $n$이 100만이어도 $d$번에 끝난다는 뜻이고, 전처리기가 노리는 것이 정확히 이 $d$(정확히는 고유값 뭉치의 개수)를 줄이는 일이다.³

문제는 순진하게 $\mathbf{r}_0, A\mathbf{r}_0, A^2\mathbf{r}_0, \dots$을 그냥 쌓으면 이 벡터들이 전부 최대 고유값 방향으로 쏠려서 거의 평행해진다는 것. 수치적으로 서로 구분이 안 되는 쓰레기가 된다. 그래서 반드시 직교화가 필요하다.

3. 직교화 — 아르놀디와 란초스

크리로프 벡터들을 그람-슈미트로 하나씩 직교정규화하며 정규직교 기저 $\{\mathbf{q}_1, \dots, \mathbf{q}_m\}$를 만드는 절차가 아르놀디 반복(Arnoldi iteration, 1951)이다. 이 과정은 다음 관계로 요약된다.

$$A Q_m = Q_{m+1} \tilde{H}_m$$

여기서 $\tilde{H}_m$은 상부 헤센베르크(Hessenberg) 형태의 작은 행렬이다. 즉 거대한 $A$의 작용을, 크리로프 부분공간 위에서는 조그만 $(m+1) \times m$ 행렬로 압축해 놓은 셈. 실제 반복법은 이 작은 행렬만 상대하면 되니 계산이 감당 가능해진다.

행렬 $A$가 대칭이면 이 헤센베르크 행렬이 삼중대각(tridiagonal)으로 붕괴하는데, 이 특수한 경우가 란초스 반복(Lanczos iteration, 1950)이다. 삼중대각이면 새 벡터를 직전 두 개하고만 직교화하면 되므로, 저장 비용과 연산량이 극적으로 줄어든다. 대칭 여부가 알고리즘 선택을 가르는 첫 번째 갈림길인 이유다.

4. 3대 알고리즘

크리로프 계열은 행렬의 성질에 따라 쓰는 물건이 다르다. 국룰은 대략 이렇다.

알고리즘	적용 행렬	특징
CG (켤레기울기법)	대칭 양정치(SPD)	짧은 점화식, 메모리 저렴, 가장 우아함
MINRES	대칭 (부정치 허용)	잔차 노름 최소화, 3항 점화식
GMRES	일반 비대칭	강력하지만 벡터를 다 저장해야 함
BiCGSTAB	일반 비대칭	짧은 점화식, 저장 저렴, 수렴이 요동칠 수 있음

CG(Conjugate Gradient, 켤레기울기법, Hestenes & Stiefel 1952)는 구조해석의 강성행렬처럼 대칭 양정치 행렬을 만났을 때의 왕도다. 삼중대각 구조 덕에 벡터 몇 개만 들고 있으면 되고, 수렴 과정이 에너지 노름을 단조 감소시키는 최적화로 해석되어 이론도 깔끔하다.

GMRES(Generalized Minimal RESidual, Saad & Schultz 1986)는 대칭성이 깨진 일반 행렬 — 대류항이 지배적인 나비에-스토크스 방정식 이산화 같은 경우 — 의 표준이다. 매 스텝 잔차 노름을 최소화하므로 수렴이 흔들리지 않고 안정적이지만, 아르놀디 벡터를 전부 저장해야 해서 반복이 길어지면 메모리가 폭발한다. 그래서 실무에서는 $m$번마다 다시 시작하는 GMRES(m) 재시작 방식을 쓴다. 재시작 주기를 너무 짧게 잡으면 수렴이 정체되는(stagnation) 함정이 있으니 주의.⁴

BiCGSTAB(Bi-Conjugate Gradient Stabilized, van der Vorst 1992)는 비대칭인데 GMRES의 메모리는 감당 못 하겠다 싶을 때의 현실적 타협안. 짧은 점화식으로 저장은 싸지만 잔차가 널을 뛰기도 한다. OpenFOAM을 돌려본 사람이라면 로그에 찍히는 PBiCGStab, PCG가 바로 이 친구들이다.

5. 수렴과 조건수

크리로프 해법의 수렴 속도는 결국 행렬 $A$의 스펙트럼, 즉 고유값 분포에 달렸다. CG의 경우 오차가 다음처럼 조건수 $\kappa(A)$에 지배된다.

$$\frac{\|\mathbf{e}_k\|_A}{\|\mathbf{e}_0\|_A} \le 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k$$

조건수 $\kappa = \lambda_{\max}/\lambda_{\min}$가 커지면 $\sqrt{\kappa}$에 비해 분자·분모의 $\mp 1$이 무의미해져 괄호 안의 수축률이 1에 딱 붙고, 그만큼 수렴이 기어간다. 이 감쇠율은 $1 - 2/\sqrt{\kappa}$ 정도라서, 반복 횟수는 대략 $\sqrt{\kappa}$에 비례한다.

5.1. 숫자로 확인해 보기

상대오차 $10^{-6}$을 목표로 필요한 반복 횟수 $k \approx \ln(2/10^{-6}) / \ln\frac{\sqrt{\kappa}+1}{\sqrt{\kappa}-1}$ 를 계산해 보자.

$\kappa$	$\sqrt{\kappa}$	수축률	$10^{-6}$까지 반복
$10^2$	10	0.818	약 73회
$10^4$	100	0.980	약 726회
$10^6$	1000	0.998	약 7,250회

조건수가 100배 뛸 때마다 반복이 딱 10배씩 는다 — $\sqrt{\kappa}$ 비례가 눈에 보인다. 3차원 푸아송 문제의 이산화 행렬은 $\kappa = O(h^{-2})$이므로 격자를 한 방향으로 2배 촘촘히 하면 조건수가 4배, 반복은 2배가 되고, 거기에 미지수가 8배 늘어난 것까지 곱해져 총 비용이 16배가 된다. 순수 크리로프 해법만으로 큰 문제를 풀 수 없는 이유가 이 곱셈이다.⁵

그래서 실전에서 크리로프 해법은 절대로 혼자 다니지 않는다. 반드시 전처리기와 짝을 이뤄, 원래 문제 대신 조건수가 확 개선된 등가 문제를 푼다. 심지어 다중격자법을 전처리기로 끼워 넣기도 하는데, 이 조합은 대규모 타원형 문제에서 격자 크기와 무관하게 일정한 반복 횟수로 수렴하는 성능을 낸다. “크리로프 + 전처리기”는 사실상 한 단어로 취급해야 한다.⁶

6. 관련 문서

• 전처리기
• 다중격자법
• 켤레기울기법 · GMRES
• 란초스 알고리즘
• 조건수 · 희소행렬
• 유한체적법 · 유한요소법
• 구조해석
• 뉴턴-랩슨법
• 나비에-스토크스 방정식
• OpenFOAM

7. Footnotes

1. 크릴로프의 1931년 논문은 “기계계의 미소진동을 결정하는 문제의 수치해법”에 관한 것으로, 특성다항식 계수를 구하려고 $\mathbf{b}, A\mathbf{b}, A^2\mathbf{b}, \dots$ 를 쌓았다. 이 접근은 정확히 위에서 말한 “전부 최대 고유값 방향으로 쏠린다”는 이유로 수치적으로 악명 높았다. 이름만 남고 방법은 바뀐 셈이다. ↩

2. $n \times n$ 밀집 행렬의 역행렬을 가우스 소거로 구하면 $O(n^3)$이다. 격자 100만 개짜리 문제에서 $n = 10^6$이면 $10^{18}$번 연산인데, 이건 초당 페타플롭 기계로도 며칠이 걸린다. 게다가 성긴 행렬의 역행렬은 대부분 꽉 찬 밀집 행렬이라 저장조차 불가능하다. 역행렬은 칠판에서만 아름답다. ↩

3. 어디까지나 이론상이다. 부동소수점 반올림 오차로 직교성이 슬금슬금 무너져서 실제로는 $d$번째에 딱 떨어지지 않는다. 애초에 $d$가 100만이면 100만 번 돌릴 거고, 그럴 바엔 반복법을 왜 쓰나. 크리로프 해법의 진짜 가치는 “유한 종료”가 아니라 “몇십 번 만에 충분히 정확한 근사해”에 있다. ↩

4. 재시작 GMRES가 정체되면 로그의 잔차가 어느 값에서 딱 멈춰서 안 내려간다. 이때 초심자는 “수렴은 신에게 맡긴다”며 반복만 늘리다 밤을 새우는데, 정답은 전처리기를 갈아 끼우거나 재시작 주기 $m$을 키우는 것이다. ↩

5. 그래서 “우리 솔버가 $O(N)$이다”라는 자랑은 반복 횟수가 격자와 무관하다는 주장이며, 그걸 실제로 해내는 물건이 다중격자법 전처리다. 나머지는 대개 $O(N^{4/3})$ 언저리에서 정직하게 고생한다. ↩

6. 그래서 논문이나 코드에서 이들을 묶어 “전처리된 크리로프 해법(Preconditioned Krylov Solver)“이라 부른다. PCG의 P, PGMRES의 P가 다 Preconditioned의 P다. 전처리기 없는 크리로프는 헬스장 안 가는 헬창 같은 거라 보면 된다. ↩