6.154 반복법과 크릴로프 부분 공간 방법

1. 직접법과 반복법의 비교

선형 시스템 \mathbf{A}\mathbf{x} = \mathbf{b}를 풀기 위한 수치적 방법은 직접법(direct method)과 반복법(iterative method)으로 크게 나뉜다. 직접법은 유한한 단계에서 정확한 해(부동소수점 오차 내에서)를 구하며, LU 분해, Cholesky 분해, QR 분해 등이 이에 해당한다. 직접법의 계산량은 \mathcal{O}(n^3)이므로, n이 수천 이상인 대규모 시스템에서는 계산 비용이 과도해진다.

반복법은 초기 추정값 \mathbf{x}_0에서 출발하여 점진적으로 정확한 해에 수렴하는 일련의 근사해 \mathbf{x}_1, \mathbf{x}_2, \ldots를 생성한다. 반복법의 주요 장점은 다음과 같다.

1. 각 반복 단계에서 행렬-벡터 곱 \mathbf{A}\mathbf{v}만을 필요로 하므로, 행렬 \mathbf{A}를 명시적으로 저장하지 않아도 된다.
2. 희소 행렬에 대해 행렬-벡터 곱의 계산량이 \mathcal{O}(\text{nnz})이므로, 대규모 희소 시스템에서 효율적이다.
3. 높은 정확도가 불필요한 경우 적은 반복 횟수로 충분한 근사해를 얻을 수 있다.

2. 고전적 반복법

2.1 야코비 반복법

행렬 \mathbf{A}를 대각 성분 \mathbf{D}, 하삼각 성분 \mathbf{L}, 상삼각 성분 \mathbf{U}로 분해하면 \mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}이다. 야코비(Jacobi) 반복법은 다음과 같이 정의된다.

\mathbf{x}_{k+1} = \mathbf{D}^{-1}(\mathbf{b} - (\mathbf{L} + \mathbf{U})\mathbf{x}_k)

성분별로 표현하면 다음과 같다.

x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{\substack{j=1 \\ j \neq i}}^{n} a_{ij} x_j^{(k)}\right), \quad i = 1, 2, \ldots, n

야코비 반복법의 수렴 조건은 반복 행렬 \mathbf{G}_J = -\mathbf{D}^{-1}(\mathbf{L} + \mathbf{U})의 스펙트럼 반경(spectral radius)이 1보다 작은 것이다.

\rho(\mathbf{G}_J) = \max_i \vert \lambda_i(\mathbf{G}_J) \vert < 1

행렬 \mathbf{A}가 엄격한 대각 우위(strictly diagonally dominant), 즉 \vert a_{ii} \vert > \sum_{j \neq i} \vert a_{ij} \vert이면 야코비 반복법은 반드시 수렴한다.

가우스-자이델 반복법

가우스-자이델(Gauss-Seidel) 반복법은 야코비 반복법을 개선한 것으로, 각 성분을 갱신할 때 이미 계산된 최신 값을 즉시 사용한다.

x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)}\right)

행렬 형태로는 다음과 같다.

\mathbf{x}_{k+1} = (\mathbf{D} + \mathbf{L})^{-1}(\mathbf{b} - \mathbf{U}\mathbf{x}_k)

가우스-자이델 반복법의 반복 행렬은 \mathbf{G}_{GS} = -(\mathbf{D} + \mathbf{L})^{-1}\mathbf{U}이며, \mathbf{A}가 대칭 양의 정부호(symmetric positive definite, SPD)이면 항상 수렴한다.

연속 과완화법(SOR)

연속 과완화법(Successive Over-Relaxation, SOR)은 가우스-자이델 반복법에 완화 매개변수 \omega를 도입한 것이다.

x_i^{(k+1)} = (1 - \omega)x_i^{(k)} + \frac{\omega}{a_{ii}}\left(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)}\right)

0 < \omega < 2일 때 수렴이 보장되며(SPD 행렬의 경우), \omega = 1이면 가우스-자이델 반복법과 동일하다. 최적 \omega 값은 반복 행렬의 스펙트럼 반경을 최소화하도록 선택되며, 많은 경우 이론적 최적값의 해석적 유도가 가능하다.

3. 크릴로프 부분 공간의 정의

고전적 반복법의 수렴 속도가 느린 경우가 많으므로, 현대 수치 선형대수학에서는 크릴로프 부분 공간(Krylov subspace) 방법을 주로 사용한다. 행렬 \mathbf{A} \in \mathbb{R}^{n \times n}와 벡터 \mathbf{r} \in \mathbb{R}^n에 대하여, m차 크릴로프 부분 공간은 다음과 같이 정의된다.

\mathcal{K}_m(\mathbf{A}, \mathbf{r}) = \text{span}\{\mathbf{r}, \mathbf{A}\mathbf{r}, \mathbf{A}^2\mathbf{r}, \ldots, \mathbf{A}^{m-1}\mathbf{r}\}

이 부분 공간은 행렬-벡터 곱의 반복적 적용으로 자연스럽게 생성되며, 차원이 m까지 증가한다. 크릴로프 부분 공간 방법의 핵심 아이디어는 전체 n차원 공간 대신 저차원의 크릴로프 부분 공간에서 최적 근사해를 탐색하는 것이다.

크릴로프 부분 공간이 효과적인 이유는 다음과 같다. 임의의 다항식 p(\mathbf{A})를 \mathbf{r}에 적용한 결과는 크릴로프 부분 공간에 속한다.

p(\mathbf{A})\mathbf{r} = \sum_{j=0}^{m-1} c_j \mathbf{A}^j \mathbf{r} \in \mathcal{K}_m(\mathbf{A}, \mathbf{r})

따라서 크릴로프 부분 공간에서의 최적화는 행렬 다항식에 의한 최적 필터링과 동치이다.

4. 켤레 기울기법(CG)

4.1 알고리즘의 유도

켤레 기울기법(Conjugate Gradient, CG)은 SPD 행렬에 대한 가장 대표적인 크릴로프 부분 공간 방법이다. CG는 \mathbf{x}_k \in \mathbf{x}_0 + \mathcal{K}_k(\mathbf{A}, \mathbf{r}_0)에서 \mathbf{A}-노름 오차를 최소화하는 해를 구한다.

\mathbf{x}_k = \arg\min_{\mathbf{x} \in \mathbf{x}_0 + \mathcal{K}_k} \|\mathbf{x} - \mathbf{x}^*\|_\mathbf{A}

여기서 \|\mathbf{v}\|_\mathbf{A} = \sqrt{\mathbf{v}^\top \mathbf{A} \mathbf{v}}이고 \mathbf{x}^* = \mathbf{A}^{-1}\mathbf{b}는 정확한 해이다. 이는 이차 함수의 최소화와 동치이다.

\phi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^\top \mathbf{A}\mathbf{x} - \mathbf{b}^\top \mathbf{x}

4.2 CG 알고리즘

CG 알고리즘의 절차는 다음과 같다.

초기화: \mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0, \mathbf{p}_0 = \mathbf{r}_0

k = 0, 1, 2, \ldots에 대해 반복한다.

\alpha_k = \frac{\mathbf{r}_k^\top \mathbf{r}_k}{\mathbf{p}_k^\top \mathbf{A}\mathbf{p}_k}

\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k

\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k \mathbf{A}\mathbf{p}_k

\beta_k = \frac{\mathbf{r}_{k+1}^\top \mathbf{r}_{k+1}}{\mathbf{r}_k^\top \mathbf{r}_k}

\mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k

수렴 판정: \|\mathbf{r}_{k+1}\| / \|\mathbf{b}\| < \epsilon이면 종료한다.

CG의 계산량

각 CG 반복 단계에서 필요한 주요 연산은 다음과 같다.

희소 행렬에서 \text{nnz} = \mathcal{O}(n)이면 각 반복의 비용은 \mathcal{O}(n)이며, m번 반복 시 총 비용은 \mathcal{O}(mn)이다. m \ll n이면 직접법의 \mathcal{O}(n^3)에 비해 극적인 절감이 이루어진다.

CG의 수렴 해석

CG의 수렴 속도는 \mathbf{A}의 조건수(condition number) \kappa(\mathbf{A}) = \lambda_{\max} / \lambda_{\min}에 의해 결정된다. k번째 반복에서의 오차 상한은 다음과 같다.

\frac{\|\mathbf{x}_k - \mathbf{x}^*\|_\mathbf{A}}{\|\mathbf{x}_0 - \mathbf{x}^*\|_\mathbf{A}} \leq 2\left(\frac{\sqrt{\kappa(\mathbf{A})} - 1}{\sqrt{\kappa(\mathbf{A})} + 1}\right)^k

잔차를 \epsilon 이하로 줄이기 위해 필요한 반복 횟수는 대략 다음과 같다.

k \approx \frac{1}{2}\sqrt{\kappa(\mathbf{A})} \ln\left(\frac{2}{\epsilon}\right)

따라서 조건수가 클수록 수렴이 느려지며, 이를 개선하기 위해 전처리(preconditioning) 기법이 필수적으로 사용된다.

GMRES

비대칭 시스템에 대한 일반화

행렬 \mathbf{A}가 대칭이 아닌 경우 CG를 직접 적용할 수 없다. 일반화된 최소 잔차법(Generalized Minimal Residual, GMRES)은 비대칭 행렬에 대한 크릴로프 부분 공간 방법으로, \mathcal{K}_k(\mathbf{A}, \mathbf{r}_0) 내에서 잔차 노름을 최소화하는 해를 구한다.

\mathbf{x}_k = \arg\min_{\mathbf{x} \in \mathbf{x}_0 + \mathcal{K}_k} \|\mathbf{b} - \mathbf{A}\mathbf{x}\|_2

4.3 아르놀디 과정

GMRES는 아르놀디(Arnoldi) 과정을 통해 크릴로프 부분 공간의 정규직교 기저 \{\mathbf{q}_1, \mathbf{q}_2, \ldots, \mathbf{q}_k\}를 구축한다. 아르놀디 과정의 j번째 단계는 다음과 같다.

\tilde{\mathbf{q}}_{j+1} = \mathbf{A}\mathbf{q}_j - \sum_{i=1}^{j} h_{ij}\mathbf{q}_i

h_{ij} = \mathbf{q}_i^\top \mathbf{A}\mathbf{q}_j, \quad i = 1, \ldots, j

h_{j+1,j} = \|\tilde{\mathbf{q}}_{j+1}\|_2, \quad \mathbf{q}_{j+1} = \frac{\tilde{\mathbf{q}}_{j+1}}{h_{j+1,j}}

이 과정은 다음의 아르놀디 분해 관계를 만족한다.

\mathbf{A}\mathbf{Q}_k = \mathbf{Q}_{k+1}\tilde{\mathbf{H}}_k

여기서 \mathbf{Q}_k = [\mathbf{q}_1, \ldots, \mathbf{q}_k] \in \mathbb{R}^{n \times k}이고 \tilde{\mathbf{H}}_k \in \mathbb{R}^{(k+1) \times k}는 상 헤센베르그(upper Hessenberg) 행렬이다.

4.4 GMRES의 최소화 문제

\mathbf{x}_k = \mathbf{x}_0 + \mathbf{Q}_k \mathbf{y}_k로 표현하면, 잔차 최소화 문제는 소규모 최소자승 문제로 변환된다.

\mathbf{y}_k = \arg\min_{\mathbf{y} \in \mathbb{R}^k} \|\beta \mathbf{e}_1 - \tilde{\mathbf{H}}_k \mathbf{y}\|_2

여기서 \beta = \|\mathbf{r}_0\|_2이다. 이 (k+1) \times k 최소자승 문제는 기븐스 회전(Givens rotation)을 이용하여 \mathcal{O}(k^2)에 풀 수 있다.

GMRES의 계산량과 재시작

GMRES의 k번째 반복에서의 주요 비용은 아르놀디 과정의 직교화이다. 수정된 그람-슈미트(modified Gram-Schmidt) 과정을 사용하면 각 반복의 비용은 \mathcal{O}(kn)이며, k번 반복까지의 총 비용은 \mathcal{O}(k^2 n)이다. 또한 k개의 기저 벡터를 저장해야 하므로 메모리 요구량은 \mathcal{O}(kn)이다.

k가 커지면 비용과 메모리가 급증하므로, 실무에서는 재시작 GMRES (GMRES(m))를 사용한다. m번 반복 후 현재 근사해를 초기값으로 하여 아르놀디 과정을 재시작한다.

쌍켤레 기울기법(BiCG)과 변형

비대칭 행렬에 대한 또 다른 접근으로 쌍켤레 기울기법(Biconjugate Gradient, BiCG)이 있다. BiCG는 \mathbf{A}와 \mathbf{A}^\top 양쪽에 대해 크릴로프 부분 공간을 동시에 구축하며, 란초스(Lanczos) 쌍직교화 과정에 기반한다.

BiCG의 안정화된 변형인 BiCGSTAB(BiCG Stabilized)은 \mathbf{A}^\top과의 곱을 피하면서 잔차의 진동을 억제한다. 각 반복 단계에서 행렬-벡터 곱 2회를 필요로 하며, 수렴 행동이 GMRES보다 불규칙할 수 있으나 메모리 요구량이 \mathcal{O}(n)으로 일정하다는 장점이 있다.

크릴로프 방법의 비교

로봇공학에서의 적용

접촉 역학 시뮬레이션

다수의 접촉점이 존재하는 로봇 시뮬레이션에서는 선형 상보성 문제(Linear Complementarity Problem, LCP)를 풀어야 하며, 이 과정에서 대규모 선형 시스템이 반복적으로 등장한다. 접촉 자코비안 \mathbf{J}_c와 관성 행렬 \mathbf{M}으로부터 형성되는 슈어 보수(Schur complement) 행렬 \mathbf{S} = \mathbf{J}_c \mathbf{M}^{-1} \mathbf{J}_c^\top에 대한 시스템은 CG 또는 GMRES로 효율적으로 풀 수 있다.

유한 요소 기반 변형 해석

소프트 로봇이나 유연 링크 로봇의 유한 요소(FE) 해석에서는 강성 행렬(stiffness matrix) \mathbf{K}와 질량 행렬이 매우 큰 희소 SPD 행렬이 된다. 이러한 시스템에서 전처리된 CG(Preconditioned CG, PCG)는 가장 표준적인 풀이 방법이다.

최적 제어와 궤적 최적화

모델 예측 제어(Model Predictive Control, MPC)에서는 매 제어 주기마다 이차 프로그래밍(Quadratic Programming, QP) 문제를 풀어야 하며, 내부적으로 KKT 시스템의 선형 풀이가 요구된다. KKT 행렬은 대칭이지만 부정치(indefinite)일 수 있으므로 MINRES가 적합하며, 문제의 구조를 이용한 전처리가 수렴 속도를 크게 향상시킨다.

참고문헌

• Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM.
• Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
• Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. SIAM.
• Shewchuk, J. R. (1994). “An Introduction to the Conjugate Gradient Method Without the Agonizing Pain.” Technical Report, Carnegie Mellon University.
• Featherstone, R. (2008). Rigid Body Dynamics Algorithms. Springer.

v 0.1