켤레기울기법

1. 개요

켤레기울기법 Conjugate Gradient Method (CG)
제안	헤스테네스 & 슈티펠 (1952), 미국 표준국
대상	대칭 양정치(SPD) 선형계 Ax = b
정체	크리로프 부분공간법 — A-노름 최소화
반복당 비용	행렬-벡터곱 1회 + 내적 2회, 벡터 4개 저장
수렴 상계	2((√κ−1)/(√κ+1))k
실무	전처리기가 전부 — PCG, AMG-PCG
비선형 확장	플레처-리브스, 폴락-리비에르

켤레기울기법(CG)은 대칭 양정치 행렬에 대한 선형계 $A\mathbf{x}=\mathbf{b}$ 를, $A$ 에 대해 서로 켤레인 방향들을 차례로 따라가며 푸는 반복 해법이다. 헤스테네스와 슈티펠이 1952년 미국 표준국에서 (서로 독립적으로 같은 것을 발견하고는 합쳐서) 발표했다.

$A$ 가 SPD면 선형계 풀이는 이차형식의 최소화와 완전히 같은 문제다.

$$\phi(\mathbf{x}) = \tfrac12 \mathbf{x}^\top A \mathbf{x} - \mathbf{b}^\top \mathbf{x}, \qquad \nabla\phi(\mathbf{x}) = A\mathbf{x} - \mathbf{b} = -\mathbf{r}$$

$\phi$ 는 강볼록이므로 유일한 최소점이 $A^{-1}\mathbf{b}$ 다. 그러니 경사하강법을 그냥 돌리면 되지 않느냐 — 여기서 CG의 존재 이유가 나온다. 최급강하법은 같은 방향을 몇 번이고 다시 밟는다. CG는 “이미 처리한 방향은 두 번 다시 건드리지 않도록” 방향을 고르고, 그 대가로 이론상 $n$ 스텝 안에 정확한 해에 도달한다.

CG가 구조해석의 강성행렬, 전산유체역학의 압력 포아송 방정식, 영상 복원의 정규화 최소제곱에서 사실상 기본값인 이유는 세 가지다. 행렬을 만들 필요가 없고($A\mathbf{v}$ 만 계산하면 된다), 메모리가 벡터 몇 개뿐이며, 전처리기를 끼울 자리가 깔끔하게 열려 있다.¹

2. A-켤레 방향

두 방향 $\mathbf{p}_i, \mathbf{p}_j$ 가 $\mathbf{p}_i^\top A \mathbf{p}_j = 0$ 을 만족하면 $A$-켤레(A-conjugate) 또는 $A$-직교라 한다. $A=I$ 면 그냥 직교다.

왜 이 조건이 특별한가. $\{\mathbf{p}_0,\dots,\mathbf{p}_{n-1}\}$ 이 서로 $A$-켤레이면 이 방향들을 좌표축 삼아 $\phi$ 를 다시 쓸 때 교차항이 전부 사라진다. 즉 $\mathbf{x} = \sum_k \alpha_k \mathbf{p}_k$ 로 두면

$$\phi\Bigl(\sum_k \alpha_k \mathbf{p}_k\Bigr) = \sum_k \Bigl( \tfrac12 \alpha_k^2\, \mathbf{p}_k^\top A \mathbf{p}_k - \alpha_k\, \mathbf{b}^\top \mathbf{p}_k \Bigr)$$

로 $n$ 개의 독립적인 1변수 문제로 분리된다. 각 방향에서 정확한 1차원 최소화를 한 번씩만 하면 끝이고, 한 번 최적화한 방향은 다른 방향의 갱신에 영향받지 않는다. 이것이 “이미 처리한 방향을 다시 밟지 않는다”의 정확한 의미다.

대조군이 최급강하법이다. 잔차 방향 $\mathbf{r}_k$ 로 정확한 라인서치를 하면 다음 잔차가 직전 방향에 유클리드 직교($\mathbf{r}_{k+1}^\top\mathbf{r}_k=0$)해지지만, $A$-직교하지는 않는다. 그래서 다음 스텝이 직전 스텝이 이룬 성과를 일부 되돌리고, 등고선이 길쭉한 타원일수록 골짜기를 가로지르며 지그재그로 기어간다. 수렴 인자가 $\frac{\kappa-1}{\kappa+1}$ 로 조건수에 그대로 지배되는 반면, CG는 $\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$ 로 제곱근이 붙는다. 조건수가 $10^4$ 이면 실효 조건수가 $10^2$ 로 줄어드는 셈이다.

3. 알고리즘 — 짧은 점화식이라는 기적

$A$-켤레 방향 $n$ 개를 미리 만들려면 그람-슈미트를 $A$-내적으로 돌려야 하고, 그러면 방향을 전부 저장해야 해서 $O(n^2)$ 메모리가 든다. CG의 진짜 마법은 그럴 필요가 없다는 것이다.

$$\begin{aligned} &\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0, \quad \mathbf{p}_0 = \mathbf{r}_0 \\ &\alpha_k = \frac{\mathbf{r}_k^\top \mathbf{r}_k}{\mathbf{p}_k^\top A \mathbf{p}_k}, \qquad \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k, \qquad \mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A\mathbf{p}_k \\ &\beta_k = \frac{\mathbf{r}_{k+1}^\top \mathbf{r}_{k+1}}{\mathbf{r}_k^\top \mathbf{r}_k}, \qquad \mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \end{aligned}$$

새 방향은 새 잔차와 직전 방향, 딱 둘만 조합해서 만든다. 그런데도 $\mathbf{p}_{k+1}$ 은 이전의 모든 방향에 대해 자동으로 $A$-켤레이고, 잔차들은 서로 직교한다. 저장할 것은 $\mathbf{x},\mathbf{r},\mathbf{p}$ 와 임시 $A\mathbf{p}$ 넷뿐이고, 반복당 비용은 행렬-벡터곱 1회와 내적 2회다.

왜 이게 가능한가. CG는 크리로프 부분공간법이고, 대칭 행렬에서 아놀디 반복이 삼중대각으로 붕괴하는 것(란초스 알고리즘)이 그 근거다. 삼중대각이란 곧 “새 기저벡터는 직전 두 개하고만 직교화하면 된다”는 뜻이며, CG의 3항 점화식이 정확히 그 그림자다. 그리고 이건 대칭에만 있는 특권이다 — 페이버-만토이펠 정리는 “짧은 점화식 + 최적성”을 동시에 갖는 크리로프 방법이 본질적으로 에르미트 행렬에만 존재함을 말한다. GMRES가 벡터를 전부 저장해야 하는 것은 구현이 게을러서가 아니라 정리 때문이다.

기하적으로 CG가 하는 일은 이렇게 요약된다.

$$\mathbf{x}_k = \arg\min_{\mathbf{x} \in \mathbf{x}_0 + \mathcal{K}_k(A,\mathbf{r}_0)}\ \|\mathbf{x} - \mathbf{x}^\star\|_A, \qquad \|\mathbf{e}\|_A^2 = \mathbf{e}^\top A \mathbf{e}$$

$k$ 차원 크리로프 부분공간 안에서 $A$-노름 오차를 최소화하는 점을 매 스텝 정확히 찾아 준다. $\|\mathbf{e}\|_A^2 = 2(\phi(\mathbf{x})-\phi(\mathbf{x}^\star))$ 이므로 이 노름은 에너지 오차이며, 구조 문제에서는 문자 그대로 변형에너지 오차다.

4. 수렴 — 상계보다 스펙트럼

최적성 성질을 다항식으로 번역하면 수렴 이론이 나온다. $\mathbf{x}_k$ 는 $\mathbf{e}_k = q_k(A)\mathbf{e}_0$ 꼴이고 $q_k$ 는 $q_k(0)=1$ 인 $k$ 차 다항식이므로,

$$\frac{\|\mathbf{e}_k\|_A}{\|\mathbf{e}_0\|_A} \;\le\; \min_{\substack{q \in \mathcal{P}_k \\ q(0)=1}}\ \max_{\lambda \in \sigma(A)} |q(\lambda)|$$

**”$A$ 의 고유값 전부에서 작고 원점에서 1인 다항식이 존재하느냐”**가 수렴 속도의 전부다. 여기에 $\sigma(A) \subset [\lambda_{\min},\lambda_{\max}]$ 만 알고 체비쇼프 다항식을 대입하면 그 유명한 상계가 나온다.

$$\frac{\|\mathbf{e}_k\|_A}{\|\mathbf{e}_0\|_A} \;\le\; 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{k}, \qquad \kappa = \frac{\lambda_{\max}}{\lambda_{\min}}$$

그런데 실무에서 CG의 거동은 이 상계보다 훨씬 좋은 경우가 흔하다. 상계가 스펙트럼 정보를 “구간 하나”로 뭉개 버렸기 때문이다.

• 고유값이 뭉쳐 있으면 빠르다. $A$ 의 서로 다른 고유값이 $m$ 개뿐이면, 그 $m$ 개 점에서 0이 되는 다항식을 만들 수 있으므로 CG는 정확히 $m$ 스텝 만에 끝난다(정확 연산에서). 뭉치(cluster)가 $m$ 개면 $m$ 스텝쯤에 오차가 뚝 떨어진다. 전처리기의 목표가 “조건수를 줄인다”보다 “고유값을 뭉친다”로 서술되는 이유가 이것이다.
• 초선형 거동. 반복이 진행되면 CG는 스펙트럼 양 끝의 고유값에 대응하는 성분을 먼저 처리하고, 그 뒤로는 남은 좁은 구간에 대한 실효 조건수로 움직인다. 그래서 수렴 곡선이 로그 스케일에서 아래로 휘는 모습이 자주 보인다.
• 극단적 고유값 몇 개가 붙어 있으면 느리다. 조건수가 같아도 스펙트럼이 $[\lambda_{\min},\lambda_{\max}]$ 에 고르게 깔린 경우와 양 끝에 몇 개만 떨어져 있는 경우의 거동은 완전히 다르다. 조건수 하나로 CG를 예측하려는 시도는 대개 실망으로 끝난다.

4.1. 검산 — 반복 수는 대체 몇 번인가

상계를 숫자로 풀면 CG를 쓰는 이유가 한눈에 보인다. 상대오차 $\varepsilon = 10^{-6}$ 까지 줄이는 데 필요한 반복 수는 각각

$$k_{\mathrm{CG}} \approx \tfrac12\sqrt{\kappa}\,\ln\frac{2}{\varepsilon}, \qquad k_{\mathrm{SD}} \approx \tfrac12\kappa\,\ln\frac{1}{\varepsilon}$$

이고, 상계 부등식을 그대로 풀어 얻은 값은 다음과 같다.

$\kappa(A)$	CG 상계	최급강하법 상계	배수
$10^2$	73	691	10배
$10^4$	726	69,078	95배
$10^6$	7,255	6,907,756	952배

조건수가 100배 나빠질 때 CG는 10배만 느려지고 최급강하법은 100배 느려진다. 제곱근 하나가 만드는 격차가 이것이다. 그리고 $\kappa = 10^6$ 짜리 문제에서 최급강하법이 요구하는 690만 번은 그냥 “안 돌아간다”는 뜻이다.

같은 표를 반대로 읽는 것이 더 실용적이다. 반복 수를 절반으로 줄이고 싶으면 조건수를 4분의 1로 만들어야 한다. 잔차 로그가 안 내려간다고 최대 반복 수를 2배로 키우는 것과, 전처리기로 $\kappa$ 를 4분의 1로 줄이는 것이 같은 값이라는 얘기다. 전자는 시간을 2배 쓰고 후자는 안 쓴다.

5. 유한 종료 성질이 왜 안 지켜지나

이론상 CG는 $n$ 번 안에 정확한 해에 도달한다. 잔차들이 서로 직교하고 $\mathbb{R}^n$ 에는 직교벡터가 $n$ 개까지만 있으니 당연하다. 그러나 부동소수점 연산에서는 이 성질이 거의 반드시 깨진다.

원인은 잔차 직교성의 상실이다. $\mathbf{r}_{k+1}=\mathbf{r}_k-\alpha_k A\mathbf{p}_k$ 로 잔차를 재귀적으로 갱신하는데, 반올림 오차가 매 스텝 조금씩 섞이면서 오래된 잔차들과의 직교성이 서서히 무너진다. 결과는 두 가지다.

• 수렴 지연. 이미 처리했어야 할 고유방향 성분이 되살아나(고스트 고유값) 같은 일을 다시 한다. 페이지의 란초스 해석과 그린바움의 결과가 이 현상을 정량화했는데, 유한정밀도 CG는 원래 $A$ 의 고유값 주위에 작은 뭉치를 붙인 더 큰 행렬에 대한 정확 CG처럼 행동한다는 것이 그 요지다. 뭉치 안에서 여러 번 반복하는 것이 지연의 정체다.
• 참잔차와 재귀잔차의 괴리. 반복이 길어지면 코드가 들고 있는 $\mathbf{r}_k$ 와 실제 $\mathbf{b}-A\mathbf{x}_k$ 가 어긋난다. 로그의 잔차는 $10^{-12}$ 인데 다시 계산해 보면 $10^{-7}$ 인 상황이 여기서 나온다. 수렴 판정 직전에 참잔차를 한 번 계산해 보는 습관이 정신건강에 좋다.

그래서 실무에서 CG는 “정확히 푸는 직접법”이 아니라 **“몇십 번 만에 충분히 정확한 답을 주는 반복법”**으로 쓴다. $n=10^7$ 인 문제를 $10^7$ 번 돌릴 생각은 애초에 아무도 하지 않는다.²

6. 전처리 — 사실상 실무의 전부

$M$ 이 SPD이고 $M\mathbf{z}=\mathbf{r}$ 을 푸는 것이 싸면, CG를 $M$-내적 위에서 돌리는 것으로 조건수를 $\kappa(M^{-1}A)$ 로 바꿀 수 있다. 알고리즘 변경은 매 반복 $\mathbf{z}_k = M^{-1}\mathbf{r}_k$ 를 한 번 풀고 $\mathbf{r}^\top\mathbf{r}$ 을 $\mathbf{r}^\top\mathbf{z}$ 로 바꾸는 것뿐이다. 이것이 PCG이고, 현실의 CG는 100% 이 형태다.

• 야코비(대각) 전처리. $M=\mathrm{diag}(A)$. 공짜에 가깝고 병렬화가 완벽하지만 효과도 딱 그만큼. 물성이 셀마다 크게 다른 문제(재료가 섞인 구조, 밀도비 큰 다상 유동)에서는 스케일만 맞춰 줘도 의외로 값을 한다.
• 불완전 촐레스키 분해(IC). 촐레스키 분해를 하되 채움(fill-in)을 버려 희소성을 유지한다. SPD 문제의 국민 전처리기지만, $A$ 가 M-행렬이 아니면 분해 도중 음의 피벗이 나와 실패할 수 있다. 대각을 $\alpha$ 만큼 밀어 주는 시프트 IC가 표준 처방이다.
• 대수적 다중격자(AMG). 격자 정보 없이 행렬 성분만 보고 성긴 레벨을 만든다. 타원형 문제에서 격자를 아무리 촘촘히 해도 반복 횟수가 거의 일정해지는 $h$-독립성을 준다. 셋업 비용이 비싸지만 같은 행렬로 여러 번 풀거나 시간전진 스텝마다 재사용하면 압도적이다. 대규모 압력 포아송 방정식의 표준 답안이 AMG-PCG인 이유.

한 가지 함정. 전처리기는 SPD여야 하고, 반복 중에 바뀌면 안 된다. 내부에서 반복 해법을 근사로 돌리는 전처리기(예: 안쪽 CG 몇 번)는 매 반복 다른 연산자가 되어 CG의 켤레성 가정을 깨뜨린다. 이 경우 유연 CG(FCG)나 GMRES 같은 유연 계열로 갈아타야 한다.

7. 비선형 CG

일반 비이차 함수 $f(\mathbf{x})$ 의 최소화로 확장한 것이 비선형 CG다. $A\mathbf{p}$ 를 계산할 수 없으니 $\alpha_k$ 는 라인서치(강 울프 조건)로 정하고, $\beta_k$ 는 잔차 대신 기울기로 쓴다.

$$\beta_k^{\mathrm{FR}} = \frac{\|\mathbf{g}_{k+1}\|^2}{\|\mathbf{g}_k\|^2}, \qquad \beta_k^{\mathrm{PR}} = \frac{\mathbf{g}_{k+1}^\top(\mathbf{g}_{k+1}-\mathbf{g}_k)}{\|\mathbf{g}_k\|^2}$$

앞이 플레처-리브스(1964), 뒤가 폴락-리비에르(1969)다. 이차함수에서는 둘이 같은 값이지만 일반 함수에서는 성격이 갈린다.

• FR은 강 울프 라인서치와 함께 전역수렴이 증명되지만, 나쁜 방향에 걸리면 아주 작은 스텝을 오래 반복하며 회복이 느리다.
• PR은 실전에서 대체로 빠르다. 기울기가 거의 안 변하면 $\beta\approx0$ 이 되어 자동으로 최급강하로 리셋되기 때문이다. 대신 파월의 반례처럼 수렴하지 않고 순환할 수 있다. 표준 처방은 $\beta^{\mathrm{PR+}} = \max\{\beta^{\mathrm{PR}},0\}$ 으로 음수를 잘라 내는 것.
• 재시작은 선택이 아니라 필수다. 이차성이 깨지면 켤레성 보증이 사라지므로 $n$ 반복마다, 또는 $|\mathbf{g}_{k+1}^\top\mathbf{g}_k|$ 가 $\|\mathbf{g}_{k+1}\|^2$ 에 비해 커지면(파월 재시작) $\mathbf{p}=-\mathbf{g}$ 로 초기화한다.

비선형 CG의 존재 이유는 메모리다. 준-뉴턴법의 L-BFGS가 대체로 더 빠르지만 벡터 쌍을 $m$ 개 들고 있어야 하고, 비선형 CG는 서너 개면 된다. 변수가 수억 개인 위상 최적화나 대규모 학습에서 아직 자리가 있는 이유다.

8. 행렬 없는 구현과 응용

CG가 $A$ 에 요구하는 것은 $\mathbf{v} \mapsto A\mathbf{v}$ 하나뿐이다. 행렬을 조립할 필요도, 저장할 필요도 없다. 이 성질이 실무에서 갖는 무게는 생각보다 크다.

• 유한요소법의 요소별 곱. 전역 강성행렬을 만들지 않고 요소별 기여를 그때그때 더해 $A\mathbf{v}$ 를 만든다. 고차 요소일수록 행렬 저장이 비싸지므로 이 방식이 이득이고, GPU에서는 메모리 대역폭이 병목이라 더욱 그렇다.
• 압력 포아송 방정식. 비압축성 유동의 투영법이나 SIMPLE 알고리즘의 압력 보정 단계는 대칭 양정치(정확히는 노이만 경계에서 상수 모드만큼 준정치)라 PCG의 정확한 사냥감이다. 실제 CFD 해석 시간의 상당 부분이 이 한 방정식에 들어가며, OpenFOAM 로그의 PCG/GAMG가 바로 그 자리다. 순수 노이만 문제는 상수 모드를 제거하거나 한 점의 압력을 고정해 특이성을 처리해야 한다.
• 뉴턴-CG(절단 뉴턴법). 뉴턴 방향을 구하는 내부 선형계를 CG로 대충 푼다. 헤세 행렬을 만들지 않고 헤세-벡터곱만 자동 미분이나 유한차분으로 계산하면 되고, 초기 반복에서는 정확도를 낮게 잡아 낭비를 막는다. 신뢰영역과 결합할 때는 CG 도중 음의 곡률($\mathbf{p}^\top A\mathbf{p}\le0$)을 만나면 그 방향으로 신뢰영역 경계까지 가고 멈추는 슈타이하우크 규칙을 쓴다 — 비볼록에서 CG가 터지는 지점을 오히려 정보로 재활용하는 영리한 처리다.

9. 함정

• $A$ 가 SPD가 아니면 쓰면 안 된다. 부정치면 $\mathbf{p}^\top A\mathbf{p}$ 가 0이나 음수가 되어 그냥 터진다. 대칭 부정치는 MINRES, 비대칭은 GMRES나 BiCGSTAB이다. “대칭인 것 같은데”는 근거가 아니다 — 경계조건 처리에서 대칭성이 깨지는 일이 흔하다.
• 정규방정식으로 밀어 넣지 마라. 최소제곱을 $A^\top A\mathbf{x}=A^\top\mathbf{b}$ 로 바꿔 CG를 돌리면 조건수가 제곱된다. LSQR·LSMR처럼 정규방정식을 명시적으로 만들지 않는 알고리즘을 쓰는 게 맞다.
• 수렴 판정은 상대잔차로. $\|\mathbf{r}_k\|\le\texttt{tol}\cdot\|\mathbf{b}\|$ 처럼 스케일에 무관한 기준을 쓰고, 좌전처리를 쓰면 로그에 찍히는 잔차가 전처리된 잔차라는 점을 기억해야 한다.
• 전처리기 없이 반복 수만 늘리는 것은 답이 아니다. 격자를 절반으로 줄이면 조건수는 네 배가 되고 반복 수는 두 배가 된다. 잔차 로그가 안 내려간다고 최대 반복 수만 키우는 것은 크리로프 부분공간법의 가장 흔한 오답이다.³

10. 관련 문서

• 크리로프 부분공간법 · GMRES · 란초스 알고리즘 · 아놀디 알고리즘
• 전처리기 · 다중격자법 · 촐레스키 분해 · 영역 분할법
• 조건수 · 희소행렬 · 고유값 문제 · 그람-슈미트
• 경사하강법 · 준-뉴턴법 · 라인서치 · 신뢰 영역 방법
• 유한요소법 · 유한체적법 · SIMPLE 알고리즘 · 포아송 방정식
• 구조해석 · 전산유체역학 · PETSc · OpenFOAM

11. Footnotes

1. 재미있는 역사가 하나 있다. 1952년 당시 CG는 ”$n$ 스텝에 끝나는 직접법”으로 소개됐고, 반올림 오차 때문에 그 약속이 안 지켜지자 20년 가까이 찬밥 신세였다. 반복법으로 다시 보고 전처리기를 붙이자는 관점(리드 1971, 마이어링크-판 데르 포르스트의 불완전 촐레스키 1977)이 자리잡고 나서야 현대적 지위를 얻었다. 알고리즘이 아니라 알고리즘을 보는 각도가 바뀐 사례. ↩

2. 그래서 “CG는 $n$ 스텝 만에 정확히 끝난다”는 문장은 시험에서는 맞고 현장에서는 무의미하다. 실제로 이 성질에 의존하는 코드는 없다. 굳이 쓸모를 찾자면 유닛 테스트에서 $5\times5$ 짜리 문제를 5스텝 안에 푸는지 확인하는 용도 정도. ↩

3. 반복 수 상한을 200에서 2000으로 올려 놓고 “돌아간다”고 보고하는 것은 문제를 해결한 게 아니라 청구서를 미룬 것이다. 격자를 한 번 더 조이는 순간 20000이 필요해진다. 반복 수가 격자와 함께 늘어나는 것이 보이면 그건 전처리기를 바꾸라는 신호이지 인내심을 시험하라는 신호가 아니다. ↩