아놀디 알고리즘

수치해석 시뮬레이션 마지막 수정: 2026-07-29 04:43:22

1. 개요[편집]

대칭이면 란초스. 아니면 아놀디. 그리고 세상의 행렬은 대부분 대칭이 아니다.

아놀디 알고리즘(Arnoldi iteration)은 대규모 비대칭(비에르미트) 행렬 $A$ 의 극단 고유값 몇 개를 구하기 위해, 크릴로프 부분공간에 $A$ 를 정사영해 훨씬 작은 상헤센베르크 행렬로 압축하는 반복법이다. 1951년 W. E. 아놀디가 “최소화 반복의 원리”라는 제목으로 발표했고, 오늘날 MATLAB eigs, SciPy eigs, 그리고 그 밑의 ARPACK까지 이어지는 계보의 출발점이다.

핵심 아이디어는 크리로프 부분공간법 전반과 같다. 시작 벡터 $v_1$ 에 대해

\mathcal{K}_m(A, v_1) = \mathrm{span}\{v_1, Av_1, A^2 v_1, \dots, A^{m-1}v_1\}

을 만들고, $A$ 의 작용을 이 $m$ 차원 공간으로 제한한다. 필요한 것은 행렬-벡터 곱 하나뿐 — $A$ 를 명시적으로 저장할 필요조차 없어서(matrix-free), 희소행렬이나 함수로만 주어지는 연산자에 그대로 쓸 수 있다.

2. 아놀디 관계와 헤센베르크 압축[편집]

크릴로프 기저 $\{v_1, Av_1, \dots\}$ 를 그대로 쓰면 벡터들이 급격히 최대 고유벡터 방향으로 쏠려 수치적으로 완전히 망가진다(사실상 거듭제곱법을 여러 번 한 것과 같아진다). 그래서 매 단계 정규직교화를 한다. 수정 그람-슈미트(MGS)로 $Av_j$ 에서 이전 $v_1,\dots,v_j$ 성분을 하나씩 빼고 정규화하면, 그 과정에서 얻은 계수 $h_{ij}$ 가 상헤센베르크 행렬을 채운다.

A V_m = V_m H_m + h_{m+1,m}\, v_{m+1} e_m^{T} = V_{m+1} \tilde{H}_m

여기서 $V_m = [v_1 \dots v_m]$ 은 $n \times m$ 정규직교 기저, $H_m = V_m^{*} A V_m$ 은 $m \times m$ 상헤센베르크(부대각선 아래가 전부 0), $\tilde{H}_m$ 은 마지막 행에 $h_{m+1,m}$ 하나만 더 붙은 $(m+1) \times m$ 행렬이다. 이 한 줄이 알고리즘 전부라고 해도 된다.

만약 어느 단계에서 $h_{j+1,j} = 0$ 이 되면 그건 고장이 아니라 행운의 붕괴(lucky breakdown)다. $\mathcal{K}_j$ 가 $A$ -불변 부분공간이 됐다는 뜻이고, $H_j$ 의 고유값은 전부 $A$ 의 정확한 고유값이다.

3. 리츠값 — 언제 믿을 것인가[편집]

$H_m$ 의 고유쌍 $(\theta_i, y_i)$ 에서 $\theta_i$ 를 리츠값, $\tilde{x}_i = V_m y_i$ 를 리츠벡터라 한다. 이것이 $A$ 의 고유쌍 근사다. 좋은 점은 잔차를 공짜로 알 수 있다는 것 — 위 아놀디 관계에 $y_i$ 를 곱하면

\| A\tilde{x}_i - \theta_i \tilde{x}_i \|_2 = |h_{m+1,m}| \cdot |e_m^{T} y_i|

가 되어, $n$ 차원 벡터를 만들지 않고 작은 고유벡터의 마지막 성분 하나만 보면 수렴 여부를 판정할 수 있다. 큰 행렬을 다시 곱해 볼 필요가 없다.

수렴은 스펙트럼의 바깥쪽부터 온다. 모듈러스가 큰 고유값, 실부가 최대인 고유값, 스펙트럼 볼록껍질의 꼭짓점 근처가 먼저 잡히고 내부는 한참 뒤다. 이는 크릴로프 부분공간이 본질적으로 $p(A)v_1$ 꼴의 다항식 필터이기 때문이고, 고유값 문제에서 거듭제곱법·부분공간 반복이 갖는 편향의 세련된 버전이다.

4. 란초스와의 관계, 그리고 비용[편집]

$A$ 가 에르미트이면 $H_m = V_m^{*} A V_m$ 도 에르미트다. 그런데 헤센베르크이면서 에르미트인 행렬은 삼중대각이다. 즉 직교화 계수 $h_{ij}$ 가 $|i-j| \le 1$ 일 때만 살아남고, 아놀디는 자동으로 3항 점화식 — 즉 란초스 알고리즘으로 축약된다. 아놀디는 란초스의 비대칭 일반화이고, 반대로 란초스는 아놀디의 대칭 특수화다.

문제는 대칭성을 잃으면 이 축약이 사라진다는 것. $m$ 스텝 아놀디의 비용은 다음과 같다.

항목	아놀디(비대칭)	란초스(대칭)
저장	기저 벡터 $m$ 개, $O(mn)$	원리상 벡터 3개, $O(n)$
직교화 연산	$O(m^2 n)$	$O(mn)$
스텝당 내적	$j$ 개 (누적)	2개 (고정)

$m$ 이 커질수록 스텝당 비용이 선형으로 늘고 메모리도 같이 늘어난다. 그래서 $m$ 을 무한정 키우는 것은 불가능하고, 재시작이 선택이 아니라 필수가 된다.¹

5. 재시작 — 명시적, 그리고 암시적(IRAM)[편집]

명시적 재시작은 소박하다. $m$ 스텝 돌려 원하는 리츠벡터들을 얻고, 그 조합을 새 시작 벡터 $v_1$ 으로 삼아 처음부터 다시 돌린다. 구현은 쉬운데 정보 손실이 크고, 축퇴·복소 켤레쌍 처리가 지저분하다.

암시적 재시작 아놀디(Implicitly Restarted Arnoldi Method, IRAM)는 소렌슨(1992)의 작품으로, 이 분야의 판을 바꿨다. 원하는 고유값 $k$ 개 + 여유분 $p$ 개, 즉 $m = k+p$ 스텝 분해를 만든 다음, 버리고 싶은 리츠값 $p$ 개를 시프트로 삼아 $H_m$ 에 이동 QR 스텝을 $p$ 번 적용한다. 이 조작은 시작 벡터에 필터 다항식

v_1^{+} \propto \prod_{j=1}^{p} (A - \mu_j I)\, v_1

을 곱한 것과 정확히 동등하면서도, $A$ 를 한 번도 다시 곱하지 않는다. 결과로 남는 앞쪽 $k$ 열은 그 자체로 유효한 $k$ 스텝 아놀디 분해라 거기서 이어 붙이면 된다. 원하지 않는 방향은 다항식이 눌러 죽이고 원하는 방향만 증폭되므로, 메모리를 $m$ 개로 고정한 채 수렴을 밀어붙일 수 있다. 이것을 포트란으로 구현한 것이 ARPACK이고, eigs(A, k, 'LM')을 칠 때마다 뒤에서 도는 물건이 바로 이것이다.

후속작인 Krylov–Schur(Stewart, 2001)는 헤센베르크 대신 슈어 형식을 유지한다. IRAM의 벌지 체이싱 과정에서 수렴한 리츠값을 걷어내고(purging) 잠그는(locking) 조작이 수치적으로 까다로웠던 문제를, 슈어 형식의 대각 블록 재배열로 깔끔하게 해결했다. SLEPc 같은 현대 라이브러리의 기본값이 Krylov–Schur인 이유다.

6. GMRES — 같은 기저 위의 자매 알고리즘[편집]

아놀디 기저를 고유값이 아니라 **선형계 $Ax=b$ **에 쓰면 GMRES(Generalized Minimal RESidual, Saad & Schultz 1986)가 된다. $x_0 + \mathcal{K}_m$ 안에서 잔차 2-노름을 최소화하는데, $r = b - Ax = V_{m+1}(\beta e_1 - \tilde{H}_m y)$ 이므로 문제가

\min_{y} \| \beta e_1 - \tilde{H}_m y \|_2

라는 $(m+1) \times m$ 작은 최소제곱으로 줄어든다. $\tilde{H}_m$ 이 헤센베르크라 기븐스 회전으로 한 열씩 QR을 갱신할 수 있고, 그 부산물로 $x$ 를 만들지 않고도 잔차 노름이 매 반복 공짜로 나온다. 같은 아놀디 관계에서 고유값을 보면 아놀디, 잔차를 최소화하면 GMRES — 둘은 한 뿌리다. 저장·연산이 $m$ 과 함께 자라는 문제도 똑같이 물려받아, 실무에서는 GMRES( $m$ ) 재시작을 쓴다.

7. 실무 이슈[편집]

재직교화. 고전 그람-슈미트는 직교성을 유한 정밀도에서 빠르게 잃는다. MGS가 기본이고, 그마저 부족하면 재직교화를 한 번 더 한다. 전량 재직교화는 비싸므로 DGKS 판정(직교화 후 노름이 이전의 $1/\sqrt{2}$ 배 미만으로 줄면 상쇄가 심했다고 보고 한 번 더 돌린다)으로 선택적으로 건다. 병렬 환경에서는 통신을 줄이려 고전 GS + 재직교화(CGS2)를 일부러 쓰기도 한다.
이동-역변환(shift-invert). 내부 고유값을 노릴 때는 $A$ 대신 $(A - \sigma I)^{-1}$ 에 아놀디를 돌린다. 고유값이 $1/(\lambda - \sigma)$ 로 사상되어 $\sigma$ 에 가장 가까운 $\lambda$ 가 지배적이 되므로, 극단 수렴 편향이 오히려 무기가 된다. 대가는 매 반복 선형계 풀이 — 보통 희소 LU 분해를 한 번 해두고 재사용한다. 일반화 문제 $Ax = \lambda Bx$ 에는 케일리 변환을 쓴다. 이 발상 자체는 역반복법과 같다.
비정규 행렬의 함정. $A$ 가 정규가 아니면 고유벡터 기저의 조건수가 폭발할 수 있고, 잔차가 작아도 리츠값이 참 고유값에서 멀 수 있다. 실제로 유동 안정성 해석(오어-조머펠트 방정식)처럼 강한 비정규 연산자에서는 고유값이 전부 안정 영역에 있어도 유한 시간 동안 에너지가 수백 배 증폭되는 과도 성장이 일어난다. 이럴 때 봐야 하는 것은 스펙트럼이 아니라 의사스펙트럼 $\Lambda_\epsilon(A) = \{ z : \|(zI-A)^{-1}\| > 1/\epsilon \}$ 이다.²
행렬 함수. $f(A)b \approx \|b\|_2 V_m f(H_m) e_1$ 이라는 근사가 성립해, 지수 적분기의 $e^{\tau A}b$ 나 $\varphi$ -함수 계산에도 그대로 쓰인다. 아놀디를 “고유값 알고리즘”으로만 기억하면 절반만 아는 셈이다.
블록 아놀디. 축퇴 고유값이나 여러 시작 벡터가 있을 때는 블록 버전을 쓴다. 단일 벡터 크릴로프 공간은 원리적으로 축퇴 고유값을 하나밖에 못 잡는다.³

8. 관련 문서[편집]

“메모리가 남는데 왜 재시작하냐”고 물으면 답은 연산량이다. $m=500$ 쯤 가면 스텝당 내적이 500개라, 행렬-벡터 곱보다 직교화가 더 비싼 지경이 된다. 고유값 몇 개 구하려다 CPU를 직교화에 갈아 넣는 참사. ↩
트레페던의 표현을 빌리면 “비정규 행렬에게 고유값만 물어보는 것은 사람에게 주민등록번호만 물어보는 것과 비슷하다”. 계산은 되는데 알고 싶은 걸 못 알아낸다. ↩
기하적 중복도가 2 이상인 고유값의 고유공간에서, $\mathcal{K}_m(A,v_1)$ 은 $v_1$ 의 사영 방향 딱 하나만 담는다. 대칭 문제에서 란초스가 축퇴 모드를 놓치는 고전적 사고와 같은 원인이다. 실제로는 반올림 오차 덕에 유령처럼 늦게 나타나기도 하는데, 그걸 믿고 설계하면 안 된다. ↩