분할 정복 고유값 알고리즘

1. 개요

분할 정복 고유값 알고리즘 Divide-and-conquer eigenvalue algorithm
대상	대칭 삼중대각행렬 T (밀집 대칭은 삼중대각화 후)
착상	Cuppen (1981)
분해	T = diag(T1, T2) + ρ vvT — 랭크-1 수정
병합	세큘러 방정식의 근 n 개를 극점 사이에서 격리
안정화	Gu–Eisenstat (1995) 뢰브너 재구성
LAPACK	dstedc (드라이버 dsyevd), SVD 판은 dbdsdc

고유값 문제를 반으로 자를 수 있을 리가 없다. 자른 자국이 랭크-1이라면 이야기가 다르다.

분할 정복 고유값 알고리즘(divide-and-conquer eigenvalue algorithm)은 대칭 삼중대각행렬의 비대각 원소 하나를 떼어 내 문제를 두 개의 독립된 부분문제로 쪼개고, 재귀로 얻은 두 스펙트럼을 랭크-1 수정의 세큘러 방정식을 풀어 다시 합치는 알고리즘이다. 1981년 카우펀(J. J. M. Cuppen)이 제시했다.¹

$$T \;=\; \begin{bmatrix} T_1 & \\ & T_2 \end{bmatrix} + \rho\, v v^{\top}$$

밀집 대칭 고유값 문제는 삼중대각화로 $T$ 까지 온 다음이 진짜 승부인데, 고유벡터까지 전부 필요한 큰 문제에서는 이 알고리즘이 사실상 표준이다. LAPACK dsyevd 가 부르는 것이 바로 이것이고, 같은 일을 하는 시프트 QR 반복(QR 알고리즘)보다 수십 배 빠른 경우가 드물지 않다. 대신 $O(n^2)$ 짜리 작업공간을 요구한다 — 공짜는 없다.

2. 자르기 — 왜 랭크-1인가

$T$ 의 대각을 $\alpha_1,\dots,\alpha_n$, 비대각을 $\beta_1,\dots,\beta_{n-1}$ 이라 하고 중간 지점 $m$ 을 고른다. $\beta_m$ 이 만드는 결합은 위치 $m, m+1$ 두 곳에만 걸쳐 있으므로, $v = e_m + e_{m+1}$ 로 두면

$$T = \begin{bmatrix} T_1 & \\ & T_2\end{bmatrix} + \beta_m\, v v^{\top}$$

로 정확히 쪼개진다. 여기서 $T_1$ 은 앞쪽 $m \times m$ 블록에서 $\alpha_m \to \alpha_m - \beta_m$, $T_2$ 는 뒤쪽 블록에서 $\alpha_{m+1} \to \alpha_{m+1} - \beta_m$ 으로 고친 삼중대각행렬이다. 삼중대각이라는 구조 덕분에 절단면의 랭크가 1이라는 것이 이 알고리즘 전체의 출발점이다. 헤센베르크였다면 절단면이 랭크 $O(n)$ 이 되어 아무 이득이 없다.

두 부분문제를 재귀로 풀어 $T_i = Q_i D_i Q_i^{\top}$ 를 얻고 $Q = \mathrm{diag}(Q_1,Q_2)$, $z = Q^{\top}v$ 로 두면

$$T = Q\left(D + \rho\, z z^{\top}\right)Q^{\top}, \qquad D = \mathrm{diag}(d_1,\dots,d_n),\ \ \rho = \beta_m$$

이 되고, 남은 일은 대각행렬 + 랭크-1 의 고유분해 하나뿐이다. $z$ 는 $Q_1$ 의 마지막 행과 $Q_2$ 의 첫 행을 이어 붙인 벡터라 따로 계산할 것도 없다.

3. 세큘러 방정식

$(D + \rho zz^{\top})u = \lambda u$ 를 정리하면 $u = -\rho(z^{\top}u)(D - \lambda I)^{-1}z$ 이고, 여기에 다시 $z^{\top}$ 를 곱하면 $\lambda$ 만 남는다.

$$w(\lambda) \;=\; 1 + \rho\sum_{i=1}^{n}\frac{z_i^{2}}{d_i - \lambda} \;=\; 0$$

이것이 세큘러 방정식(secular equation, 세속방정식)이다.² 이름은 천체역학의 장주기(secular) 섭동에서 왔다. 근 하나를 구하면 고유벡터는 공짜다.

$$u_i \;=\; \frac{(D-\lambda_i I)^{-1}z}{\left\lVert (D-\lambda_i I)^{-1}z\right\rVert_2}, \qquad \text{$T$ 의 고유벡터} = Q u_i$$

$d_i$ 가 서로 다르고 모든 $z_i \ne 0$ 이면 $w$ 는 각 극점 사이 구간에서 단조증가하고 $\pm\infty$ 로 발산하므로, 극점 사이마다 근이 정확히 하나 있다. $\rho > 0$ 일 때

$$d_1 < \lambda_1 < d_2 < \lambda_2 < \cdots < d_n < \lambda_n < d_n + \rho\lVert z\rVert_2^{2}$$

로 완벽하게 격리된다. 근을 찾는 데 초기 추정값 문제가 없다는 뜻이고, 이 교대(interlacing) 구조가 알고리즘을 신뢰할 수 있게 만드는 두 번째 기둥이다.

4. 근을 안전하게 뽑기

구간이 주어졌다고 뉴턴-랩슨법을 그냥 돌리면 안 된다. $w$ 는 구간 양 끝에서 극점을 향해 수직으로 치솟는 함수라 다항식 근사인 뉴턴 접선은 번번이 구간 밖으로 튀어나간다. 실무 처방은 두 가지다.

• 이동 변수. 근을 $\lambda = d_i + \mu$ 로 놓고 $\mu$ 를 직접 푼다. $d_i - \lambda$ 를 뺄셈으로 계산하면 가까운 두 수의 상쇄로 유효숫자가 통째로 날아가는데, 이동 변수를 쓰면 그 항이 $-\mu$ 로 정확히 나온다. 가까운 고유값에서 정확도가 살아남는 유일한 이유가 이것이다.
• 유리 보간. $w$ 를 좌우 극점 두 개를 그대로 살린 $\dfrac{c_1}{d_i - \lambda} + \dfrac{c_2}{d_{i+1} - \lambda} + c_3$ 꼴로 근사해 그 방정식을 닫힌 형으로 푼다(리(Li)의 “middle way”). 근사가 원 함수의 특이점 구조를 공유하므로 반복이 구간을 벗어나지 않고, 수렴 차수도 3차다. LAPACK dlaed4 가 이 방식이며 고유값 하나당 반복 두세 번이면 끝난다.

5. 디플레이션 — 실측이 이론보다 빠른 이유

두 경우에 근을 아예 구하지 않고 통과시킬 수 있다.

• $z_i$ 가 무시할 만큼 작으면 $d_i$ 자체가 이미 고유값이고, 대응 고유벡터는 $Q$ 의 그 열 그대로다.
• $d_i \approx d_j$ 면 $(i,j)$ 평면의 기븐스 회전 하나로 $z_i, z_j$ 중 하나를 0으로 만들 수 있고, 그러면 위 경우로 환원된다.

이 디플레이션(deflation)은 실전 행렬에서 대량으로 일어난다. 문제 크기가 매 병합마다 줄어드는 셈이라 실측 복잡도가 최악 경계 $O(n^3)$ 보다 한참 아래(경험적으로 $O(n^{2.3})$ 언저리)로 내려가며, 이것이 분할 정복이 이론상의 우월성보다 훨씬 큰 실측 우위를 갖는 진짜 이유다. 물론 뒤집어 말하면 성능이 데이터에 의존한다는 뜻이라, 벤치마크 행렬을 고르는 사람이 결론도 고를 수 있다.³

6. 직교성 상실과 Gu–Eisenstat

카우펀의 원 논문이 20년 가까이 라이브러리에 들어가지 못한 이유가 있다. 고유벡터 공식 $u_i \propto (D - \lambda_i I)^{-1}z$ 의 $i$ 번째 성분은 $z_i/(d_i - \lambda_i)$ 인데, $d_i$ 와 $d_{i+1}$ 이 가까우면 $d_i - \lambda_i$ 는 거의 상쇄된 작은 수다. 계산된 $\hat\lambda_i$ 의 절대오차가 아무리 작아도 이 분모의 상대오차는 1에 육박할 수 있고, 그러면 두 고유벡터가 서로 수직이기는커녕 거의 평행하게 나온다. 대칭 문제에서 고유벡터의 직교성은 협상 대상이 아니므로 이건 치명적이다.

동가라–소렌슨(1987)의 1차 처방은 뭉친 무리마다 그람-슈미트 재직교화를 거는 것이었다. 작동은 하지만 비용이 $O(n^3)$ 으로 튀고, 무리 판정 임계값에 따라 결과가 달라진다.

구–아이젠스탯(Gu–Eisenstat, 1995)의 해법은 훨씬 영리하다. $z$ 를 고치는 것이다. 계산된 근 $\hat\lambda_1,\dots,\hat\lambda_n$ 을 확정한 뒤, 뢰브너(Löwner) 정리를 뒤집어 그 값들이 정확한 고유값이 되게 하는 새 벡터 $\tilde z$ 를 만든다.

$$\tilde z_i^{2} \;=\; \frac{\prod_{j=1}^{n}\left(\hat\lambda_j - d_i\right)}{\prod_{j \ne i}\left(d_j - d_i\right)}$$

부호는 원래 $z_i$ 에서 가져온다. 이제 $\hat\lambda_i$ 는 $D + \rho\,\tilde z\tilde z^{\top}$ 의 정확한 고유값이므로, 그 행렬에 대한 고유벡터 공식은 상쇄 없이 계산되고 결과는 기계 정밀도 수준에서 직교한다. 재직교화가 아예 필요 없다. 곱셈이 그대로 상쇄를 흡수한다는 것이 요점이고, LAPACK dlaed3 가 이 공식을 그대로 담고 있다. 분할 정복이 “빠르지만 위험한 방법”에서 “빠르고 안전한 기본값”으로 승격된 지점이다.

7. 비용·작업공간·언제 쓰나

항목	시프트 QR (dsteqr)	분할 정복 (dstedc)	이분법+역반복 (dstebz+dstein)
고유값만	$O(n^2)$	$O(n^2)$	원하는 개수 $k$ 에 $O(kn)$
고유벡터 전부	$O(n^3)$, 상수 큼	최악 $O(n^3)$, 실측 훨씬 아래	뭉치면 $O(n^3)$
추가 작업공간	$O(n)$	$O(n^2)$	$O(n)$
직교성 보장	회전 누적으로 자동	Gu–Eisenstat 로 보장	재직교화 필요
일부만 뽑기	불가	불가	가능

정리하면 선택 기준은 셋이다.

• 고유값만 필요하다 → 분할 정복의 이점이 거의 없다. QR 도 $O(n^2)$ 이다.
• 고유벡터가 전부 필요하고 $n$ 이 크다 → 분할 정복. 병합 단계의 주 연산이 $Q\,U$ 형태의 행렬-행렬 곱이라 BLAS-3 효율이 그대로 나오는 것도 크다. 반면 QR 은 기븐스 회전을 하나씩 누적해야 해 원리적으로 BLAS-1~2다.
• 몇 개만, 혹은 번호 구간만 필요하다 → 스텀 수열 기반 이분법. 분할 정복은 전부 아니면 전무다.

작업공간은 무시할 수 없는 제약이다. dstedc 는 고유벡터를 요구할 때 실수 배열 $1 + 3n + 2n\lg n + 4n^2$ 개를 요구한다 — 결과 행렬 $Q$ 말고 추가로 $4n^2$ 다. 메모리가 빠듯한 노드에서 dsyevd 가 워크스페이스 조회에 실패하는 사고가 여기서 나온다. 작은 블록에서는 오버헤드가 손해라 LAPACK 은 크기 25 이하 부분문제를 dsteqr 로 넘긴다.

8. 파생과 이웃

같은 착상은 랭크-1 구조가 보이는 곳마다 재활용됐다.

• 특이값 분해. 이중대각행렬에 같은 분할을 걸어 특이값 분해를 구하는 것이 dbdsdc 다(쌍대각화 이후 단계).
• 병렬화. 부분문제들이 완전히 독립이고 병합만 동기화가 필요해, 트리 구조 그대로 분산 메모리에 얹힌다. ScaLAPACK pdsyevd 가 이 구조이며, 1987년 동가라–소렌슨 논문 자체가 병렬 구현 논문이었다.
• MRRR (dstemr). 총 $O(n^2)$ 로 더 싸고 작업공간도 $O(n)$ 이지만, 극단적으로 뭉친 스펙트럼에서 실패 보고가 있어 라이브러리들이 기본값을 분할 정복으로 두는 경우가 많다. 가장 빠른 알고리즘이 항상 기본값이 되지는 않는다는 흔한 사례.
• 밴드행렬·일반화 문제. 밴드 대칭은 삼중대각으로 줄인 뒤 그대로 태우고, $K\phi = \lambda M\phi$ 는 촐레스키 분해로 표준형으로 옮긴 뒤 태운다(모드 해석의 표준 경로).

9. 관련 문서

• 삼중대각화 · 고유값 문제 · QR 알고리즘 · 윌킨슨 시프트
• 스텀 수열 · 역반복법 · 기븐스 회전 · 그람-슈미트
• 특이값 분해 · 쌍대각화 · 촐레스키 분해
• LAPACK · 조건수 · 후진 오차 해석 · 병렬 컴퓨팅
• 모드 해석 · 란초스 알고리즘

10. Footnotes

1. Cuppen, J. J. M. (1981). “A divide and conquer method for the symmetric tridiagonal eigenproblem”. Numerische Mathematik 36, 177–195. 네덜란드계 성이라 한글 표기가 “카우펀/쿠펀/커펜”으로 문헌마다 제각각이다. 논문 한 편으로 40년째 라이브러리 기본값을 차지하고 있는 사람치고는 이름이 알려져 있지 않은 편. ↩

2. 세큘러(secular)는 라틴어 saeculum(세대·긴 시간)에서 온 말로, 천체의 궤도 요소가 아주 긴 주기로 표류하는 현상을 가리키던 천문학 용어다. 그 섭동 계산이 마침 같은 꼴의 방정식이었던 탓에 이름만 선형대수로 이민 왔다. 고유값 계산과 행성의 장주기 표류 사이에 물리적 관계는 없다. ↩

3. 그래서 분할 정복 논문의 성능 표에는 “임의 대칭행렬”과 함께 “1-2-1 행렬”, “윌킨슨 행렬” 같은 인공 행렬이 줄줄이 등장한다. 디플레이션이 잘 되는 행렬을 고르면 $O(n^2)$ 처럼 보이고, 고유값이 고르게 퍼진 행렬을 고르면 $O(n^3)$ 에 가까워진다. 벤치마크를 볼 때 행렬 목록부터 확인하는 습관은 여기서 생긴다. ↩