BiCGSTAB

1. 개요

BiCGSTAB Bi-Conjugate Gradient Stabilized
푸는 것	일반(비대칭) 희소 선형계 Ax = b
반복당 비용	행렬-벡터 곱 2회, 내적 4회 (전치 곱 불필요)
저장	벡터 약 7~8개, 반복 수와 무관
제안	헨크 반 데르 포르스트, 1992
계보	BiCG(플레처 1976) → CGS(소네벨트 1989) → BiCGSTAB → BiCGStab(ℓ), IDR(s)
약점	고유값의 허수부가 크면 정체, 붕괴(breakdown) 가능

GMRES 는 메모리가 무섭고, BiCG 는 들쭉날쭉하고, CGS 는 더 들쭉날쭉하다. 그럼 CGS 의 제곱을 버리고 매 스텝 한 번씩 다림질하자.

BiCGSTAB(Bi-Conjugate Gradient Stabilized)은 비대칭 선형계를 푸는 크리로프 부분공간법으로, 쌍켤레 기울기법(BiCG)의 짧은 점화식을 유지하면서 매 반복 1차원 잔차 최소화 단계를 하나씩 곱해 수렴 곡선을 매끄럽게 만든 방법이다. 헨크 반 데르 포르스트가 1992년 발표했고, 저장량이 반복 수와 무관하며 전치 행렬 곱이 필요 없어서 OpenFOAM의 PBiCGStab, PETSc의 KSPBCGS 같은 이름으로 비대칭 문제의 기본 솔버 자리를 차지했다. 이 논문은 1990년대 응용수학 분야에서 가장 많이 인용된 논문으로 꼽힌다.

GMRES가 “잔차를 최소화한다”는 최적성을 사는 대가로 반복마다 벡터를 하나씩 쌓아야 한다면, BiCGSTAB 은 최적성을 포기하고 고정 메모리를 산다. 이 문서는 그 거래가 무엇인지를 다항식으로 읽고, 알고리즘 각 줄의 역할, 수렴이 갑자기 멈추는 두 가지 이유, 그리고 실측 비교로 “언제 GMRES 대신 이걸 쓰나”를 정리한다.

2. 계보 — 잔차 다항식의 이야기

크리로프 방법의 $k$ 번째 잔차는 $\mathbf r_k = P_k(A)\,\mathbf r_0$, $P_k(0) = 1$ 인 어떤 $k$ 차 다항식으로 쓸 수 있다. 각 방법은 이 다항식을 어떻게 고르느냐로 갈린다.

• BiCG(플레처 1976, 뿌리는 란초스 1952): 비대칭 란초스 알고리즘을 써서 $A$ 와 $A^{\mathsf T}$ 로 두 크리로프 공간을 동시에 키우고, 잔차를 그림자 공간(shadow space)에 직교시켜 $P_k$ 를 정한다. 3항 점화식이라 저장이 싸지만 $A^{\mathsf T}$ 곱이 필요하고, 잔차가 들쭉날쭉하다.
• CGS(소네벨트 1989): $A^{\mathsf T}$ 곱이 하는 일은 그림자 쪽 계수 계산뿐이라는 걸 간파하고, 잔차를 $\mathbf r_k = P_k(A)^2\,\mathbf r_0$ 로 제곱해 전치 곱을 없앴다. 수렴하는 구간에서는 BiCG 의 두 배 빠르지만, 튀는 구간에서도 두 배로 튄다. 제곱한 잔차가 $10^{10}$ 까지 솟구쳤다 내려오면 그 사이 반올림 오차가 해를 망친다.¹
• BiCGSTAB: 두 개의 $P_k$ 중 하나를 다른 다항식으로 바꾼다.

$$\mathbf r_k = Q_k(A)\,P_k(A)\,\mathbf r_0, \qquad Q_k(z) = \prod_{j=1}^{k} (1 - \omega_j z)$$

$P_k$ 는 여전히 BiCG 다항식이고(그래서 BiCG 의 계수를 그대로 쓴다), $Q_k$ 는 매 스텝 $\omega_j$ 를 현재 잔차 노름이 최소가 되도록 하나씩 고른 1차 인자들의 곱이다. 즉 BiCG 한 스텝 + GMRES(1)(=1차원 최소 잔차) 한 스텝을 번갈아 곱하는 구조다. 제곱처럼 진폭을 키우지 않고 오히려 매 스텝 깎아 내려서 “STAB(ilized)“이다.

3. 알고리즘 — 줄마다 하는 일

$\mathbf r_0 = \mathbf b - A\mathbf x_0$, 그림자 벡터 $\hat{\mathbf r}_0$ 를 고른다(보통 $\hat{\mathbf r}_0 = \mathbf r_0$). $\rho_0 = \alpha = \omega_0 = 1$, $\mathbf v_0 = \mathbf p_0 = 0$. $i = 1, 2, \dots$ 에 대해:

$$\begin{aligned} \rho_i &= (\hat{\mathbf r}_0,\ \mathbf r_{i-1}), \qquad \beta = \frac{\rho_i}{\rho_{i-1}}\cdot\frac{\alpha}{\omega_{i-1}} \\ \mathbf p_i &= \mathbf r_{i-1} + \beta\,(\mathbf p_{i-1} - \omega_{i-1}\mathbf v_{i-1}) \\ \mathbf v_i &= A\mathbf p_i, \qquad \alpha = \frac{\rho_i}{(\hat{\mathbf r}_0,\ \mathbf v_i)} \\ \mathbf s &= \mathbf r_{i-1} - \alpha\,\mathbf v_i \\ \mathbf t &= A\mathbf s, \qquad \omega_i = \frac{(\mathbf t,\ \mathbf s)}{(\mathbf t,\ \mathbf t)} \\ \mathbf x_i &= \mathbf x_{i-1} + \alpha\,\mathbf p_i + \omega_i\,\mathbf s, \qquad \mathbf r_i = \mathbf s - \omega_i\,\mathbf t \end{aligned}$$

• $\mathbf s$ 까지가 BiCG 반 스텝이다. $\alpha$ 는 BiCG 의 그 $\alpha$ 다.
• $\omega_i$ 는 $\|\mathbf s - \omega A\mathbf s\|_2$ 를 최소화하는 스칼라 — 한 줄짜리 최소제곱이다. 이게 안정화 스텝이다.
• 행렬-벡터 곱은 $A\mathbf p_i$, $A\mathbf s$ 두 번. 내적은 $\rho$, $(\hat{\mathbf r}_0,\mathbf v)$, $(\mathbf t,\mathbf s)$, $(\mathbf t,\mathbf t)$ 네 번(수렴 판정용 노름은 별도).
• $\mathbf s$ 가 이미 충분히 작으면 $\mathbf x_{i-1} + \alpha\mathbf p_i$ 로 끝내는 조기 탈출을 넣는 게 일반적이다.

전처리기 $M$ 은 오른쪽 전처리로 넣는 게 관례다. $\mathbf p_i$, $\mathbf s$ 에 $A$ 를 곱하기 직전 $M^{-1}$ 을 적용하고 해 갱신도 전처리된 벡터로 하면, 잔차가 원래 문제의 잔차 그대로라 정지 판정이 정직해진다. 반복당 전처리 적용도 2회다.

4. 수렴이 멈추는 두 가지 이유

① 란초스 붕괴. $\rho_i = (\hat{\mathbf r}_0, \mathbf r_{i-1}) = 0$ 이나 $(\hat{\mathbf r}_0, \mathbf v_i) = 0$ 이면 나눗셈이 터진다. BiCG 에서 물려받은 문제로, 정확히 0 은 드물지만 거의 0 이면 계수가 커져 잔차가 튄다. 그림자 벡터를 다시 잡고 재시작하는 것이 실무 처방이다.

② $\omega$ 붕괴 — 허수 고유값의 저주. $\omega_i = (\mathbf t,\mathbf s)/(\mathbf t,\mathbf t) = (A\mathbf s,\mathbf s)/\|A\mathbf s\|^2$ 이다. $A$ 의 대칭부가 약하면 $(A\mathbf s, \mathbf s) \approx 0$ 이 되어 $\omega \approx 0$, 그러면 $\mathbf x$ 갱신이 사라지고 다음 스텝의 $\beta$ 는 $1/\omega$ 로 폭발한다.

극단적 예가 반대칭 행렬이다. $A^{\mathsf T} = -A$ 면 모든 $\mathbf s$ 에 대해 $\mathbf s^{\mathsf T}A\mathbf s = 0$ 이다. $200\times200$ 무작위 반대칭 행렬로 첫 스텝을 계산하면 $\omega_1 \approx 2.8\times10^{-19}$ — 반올림 오차 수준의 0 이다. SciPy bicgstab 은 5000 회를 돌고도 수렴하지 못했고, 같은 행렬에서 GMRES 는 200 회(= 행렬 크기)에 $10^{-15}$ 까지 갔다.

근본 원인은 $Q_k$ 의 인자 $1 - \omega z$ 가 실수 계수 1차식이라는 것이다. 이 인자는 실축 근처 고유값은 잘 누르지만, 허수부가 큰 고유값 쌍은 하나의 실수 근으로 동시에 누를 수가 없다. 대류가 지배적인 유동 문제가 정확히 이 경우다. 그래서 나온 것이 BiCGStab($\ell$)(슬레이펜·포켐마 1993)으로, $\ell$ 스텝마다 $\ell$ 차 최소 잔차 다항식을 붙인다. $\ell = 2$ 면 복소 켤레 쌍을 실수 계수 2차식으로 누를 수 있다. PETSc KSPBCGSL, MATLAB bicgstabl 이 이것이다.

5. 숫자로 — 대류-확산 방정식

$[0,1]^2$ 에서 $-\Delta u + \mathrm{Pe}\,u_x = 1$ 을 $63\times63$ 내부 격자, 중심차분으로 이산화하고($n = 3969$) 전처리 없이 상대잔차 $10^{-8}$ 까지 풀었다(SciPy, $\mathbf x_0 = 0$). 셀 페클레 수 $\mathrm{Pe}\,h/2$ 가 1 을 넘으면 중심차분이 대각 우세를 잃고 고유값이 복소평면으로 크게 퍼진다.

Pe	셀 Pe	BiCGSTAB 반복 (곱)	GMRES(30) 곱	완전 GMRES 곱
0	0	82 (164)	614	117
50	0.39	97 (194)	308	124
200	1.56	150 (300)	329	107
1000	7.8	1085 (2170)	355	157

• 약한 대류(Pe ≤ 200)에서 BiCGSTAB 은 재시작 GMRES(30) 보다 행렬-벡터 곱이 적고, 벡터 8개 남짓으로 버틴다. 완전 GMRES 는 곱 수가 가장 적지만 벡터를 100개 넘게 저장한다 — 3D 백만 셀이면 그게 수 GB 다.
• 강한 대류(셀 Pe 7.8)에서 BiCGSTAB 은 곱 2170 회로 GMRES(30) 의 6 배를 썼다. 위의 “허수 고유값의 저주”가 그대로 숫자로 나온 것이다.
• 참고로 Pe = 0 은 대칭 양정치라 켤레기울기법이 정답이다. 같은 문제에서 CG 는 118 회(곱 118)로, BiCGSTAB 의 164 곱보다 싸다. 대칭 문제에 BiCGSTAB 을 쓰는 건 손해다.

물론 실전에서는 전처리 없이 돌리지 않는다. ILU 류 전처리를 붙이면 고유값이 1 근처로 모여 허수부 문제가 크게 줄어드는데, BiCGSTAB 이 CFD 코드에서 잘 버티는 이유의 상당 부분이 이것이다.²

6. 고르는 요령

• 대칭 양정치 → CG. 대칭 부정치 → MINRES. BiCGSTAB 은 이 둘이 안 될 때다.
• 비대칭, 메모리 빠듯, 전처리가 괜찮음 → BiCGSTAB. 대부분의 유한체적법 압력·운동량 방정식이 여기다.
• 비대칭, 대류 강함·전처리 약함, 잔차가 계단처럼 멈춤 → BiCGStab(2) 나 IDR(s), 아니면 메모리를 내고 GMRES.
• 수렴이 너무 들쭉날쭉해서 정지 판정이 불안 → 잔차 노름의 이동 최솟값을 기록하거나 QMR 계열(QMRCGSTAB)로 매끄럽게 한다.

IDR($s$)(소네벨트·반 기젠 2008)는 BiCGSTAB 을 $s = 1$ 특수 경우로 포함하는 일반화로, 그림자 공간을 $s$ 차원으로 늘려 고정 메모리로 GMRES 에 더 가까운 수렴을 낸다. 베설링과 소네벨트가 1980년에 낸 원래 IDR 이 CGS·BiCGSTAB 의 조상이었다는 점에서, 30년 만에 계보가 한 바퀴 돈 셈이다.³

7. 관련 문서

• 크리로프 부분공간법 · GMRES · 켤레기울기법 · MINRES
• 전처리기 · 불완전 촐레스키 · 란초스 알고리즘
• 대류-확산 방정식 · 페클레 수 · 유한체적법
• OpenFOAM · PETSc · IDR(s) · 쌍켤레 기울기법

8. Footnotes

1. CGS 의 잔차 로그를 처음 보면 심전도인 줄 안다. 반 데르 포르스트 논문 제목의 「smoothly converging」은 바로 이 심전도를 겨냥한 정면 저격이었다. ↩

2. OpenFOAM 사용자 사이에서 “수렴 안 되면 일단 PBiCGStab + DILU” 가 국룰이 된 이유다. 수렴이 안 되는 원인이 메쉬 품질이면 솔버를 백 번 바꿔도 안 된다는 것까지가 국룰의 후반부다. ↩

3. 원래 IDR(1980)은 다중격자 학회 논문집 한구석에 묻혀 있다가, 소네벨트 본인이 2008년에 일반화하면서 재조명됐다. 소네벨트는 CGS 로 한 번, IDR(s) 로 또 한 번 계보를 바꿨으니 크리로프 업계의 장기 흥행 감독이다. ↩