삼중대각화

1. 개요

삼중대각화 Tridiagonalization
대상	대칭(실) · 에르미트(복소) 행렬
결과	QTAQ = T, T 는 삼중대각 (고유값 보존)
표준 방법	하우스홀더 반사 $n-2$ 회
비용	(4⁄3)n3 flops (Q 누적 시 +(4⁄3)n3)
비대칭 대응	상헤센베르크 축약 (10⁄3)n3
LAPACK	dsytrd → dsteqr/dstedc/dstemr

대칭 고유값 문제는 두 단계다. 유한 번에 삼중대각까지 가는 앞부분, 그리고 영원히 반복하는 뒷부분.

삼중대각화(tridiagonalization)는 대칭행렬 $A$ 를 직교 상사변환으로 대각과 그 바로 위·아래 부대각만 남은 삼중대각행렬 $T$ 로 옮기는 축약이다.

$$Q^{\top}A\,Q = T = \begin{bmatrix} \alpha_1 & \beta_1 & & \\ \beta_1 & \alpha_2 & \ddots & \\ & \ddots & \ddots & \beta_{n-1} \\ & & \beta_{n-1} & \alpha_n \end{bmatrix}$$

상사변환이므로 고유값이 그대로 보존되고, $Q$ 가 직교이므로 조건수가 1이라 오차를 증폭하지 않는다. 결정적으로 이 축약은 유한 단계로 끝난다 — 하우스홀더 변환 반사 정확히 $n-2$ 번이면 도착한다.¹

밀집 대칭 고유값 solver가 예외 없이 이 형태를 취하는 이유는 단순하다. 저장이 $O(n)$ 으로 줄고, 삼중대각 전용 알고리즘들이 $O(n^2)$ 에 돌아간다. $O(n^3)$ 을 한 번 내고 나면 나머지 반복이 공짜에 가까워지는 구조다.

2. 왜 대각까지 못 가는가

“반사를 더 쓰면 대각까지 갈 수 있지 않나?” 두 겹의 이유로 안 된다.

기계적 이유. 상사변환은 양쪽에서 곱해야 한다. $k$ 번째 단계에서 $A(k+2\!:\!n,\,k)$ 를 0으로 만드는 반사 $H_k$ 는 행 $k+1$ 이하에만 작용하므로 $H_k$ 는 $e_1,\dots,e_k$ 를 고정하고, 따라서 오른쪽에서 $H_k$ 를 곱해도 방금 청소한 $k$ 열이 되살아나지 않는다. 그런데 부대각까지 지우려고 반사를 행 $k$ 이하로 넓히면 오른쪽 곱이 지운 자리를 그대로 다시 채운다. 왼손이 지우고 오른손이 복구한다. 유한 단계로 도달 가능한 한계가 정확히 삼중대각인 이유다.

원리적 이유. 대각까지 유한 번의 사칙연산과 제곱근으로 갈 수 있다면 임의의 $n$ 차 다항식의 근을 유한 번에 구할 수 있다는 뜻이 된다(동반행렬을 생각하면 된다). 아벨-루피니 정리에 정면으로 위배된다. 모든 고유값 알고리즘은 반드시 반복법이어야 한다는 것이 정리 수준의 사실이고, 축약은 그 반복이 시작되기 전에 공짜로 벌 수 있는 만큼을 벌어 두는 작업이다.²

3. 하우스홀더 축약의 실제

$k$ 번째 단계에서 $x = A(k+1\!:\!n,\,k)$ 를 $\mp\|x\|_2 e_1$ 으로 보내는 반사 $H = I - \tau vv^{\top}$ 를 만들고 $A \leftarrow HAH$ 를 적용한다. 여기서 대칭성을 쓰면 비용이 반토막 난다. $p = \tau Av$, $w = p - \tfrac{\tau}{2}(p^{\top}v)\,v$ 로 두면

$$HAH = A - v w^{\top} - w v^{\top}$$

즉 대칭 랭크-2 갱신 한 번이다. 하삼각만 건드리면 되고, $H$ 를 명시적으로 만들 일은 없다. 전체 비용은

$$\tfrac{4}{3}n^3 \ \text{flops}$$

이고, $Q$ 를 실제로 조립하려면 같은 $\tfrac{4}{3}n^3$ 이 더 든다. 반사 벡터 $v$ 들은 소거된 자리에 그대로 저장하고 스칼라 $\tau$ 만 별도 배열에 담으므로 추가 메모리는 사실상 0이다(LAPACK dsytrd).

후진 안정성도 확보된다. 계산된 $\hat T$ 는 어떤 정확한 직교 $\tilde Q$ 와 $\|E\|_F \le c(n)\,u\,\|A\|_F$ 인 $E$ 에 대해 $\hat T = \tilde Q^{\top}(A+E)\tilde Q$ 를 정확히 만족한다. 바일 정리를 얹으면 각 고유값이 최대 $\|E\|_2$ 만큼만 움직이므로, 대칭 문제에서 축약 단계는 오차 걱정의 대상이 아니다.

4. 그런데 이 단계가 병목이다

flops만 보면 완벽한데 실측 성능이 나쁘다. 랭크-2 갱신 자체는 BLAS-3로 블록화되지만, 각 패널에서 $p = \tau Av$ 를 만드는 대칭 행렬-벡터 곱(dsymv)은 원리적으로 BLAS-2다. 그리고 이 부분이 전체 연산량의 약 절반을 차지한다. 결과적으로 dsytrd 는 이론 최고 성능의 몇십 퍼센트에서 놀고, 큰 $n$ 에서는 고유값을 실제로 뽑는 뒷단계보다 준비운동인 축약이 더 오래 걸리는 역전이 일어난다.³

해법이 2단 축약(two-stage, SBR — Bischof–Lang–Sun 계열)이다.

1. 조밀 → 밴드. 대역폭 $b$ 의 밴드행렬까지만 줄인다. 전 과정이 행렬-행렬 곱이라 BLAS-3 효율이 그대로 나온다.
2. 밴드 → 삼중대각. 기븐스 회전이나 작은 반사로 밴드 밖으로 튀어나온 성분(bulge)을 대각선을 따라 끝까지 쫓아 보낸다(bulge chasing, Rutishauser–Schwarz). 연산량은 $O(n^2 b)$ 로 적지만 국소성이 나빠 파이프라인 병렬화가 필수다.

총 flops는 오히려 늘어나는데 실측 시간은 몇 배 빨라진다. ELPA·PLASMA·MAGMA 같은 현대 대규모 라이브러리가 전부 이 구조다. “연산량을 줄이는 것과 시간을 줄이는 것은 다른 문제”라는 이 바닥의 격언이 가장 노골적으로 드러나는 사례.

5. 삼중대각을 얻은 다음

$T$ 를 손에 쥐면 선택지가 넷이다. 무엇을 원하느냐로 갈린다.

가장 먼저 하는 일은 디플레이션이다. 어떤 $\beta_i$ 가 $|\beta_i| \le u(|\alpha_i| + |\alpha_{i+1}|)$ 를 만족하면 그 자리에서 행렬을 두 조각으로 쪼갠다. 문제 크기가 자동으로 줄어든다.

• 암시적 시프트 QL/QR(dsteqr). 우하단 $2\times2$ 에서 윌킨슨 시프트를 뽑아 기븐스 회전으로 bulge를 쫓는다. 대칭 삼중대각에서는 수렴이 보통 3차라 고유값당 반복이 두세 번이면 끝나고, 고유값만 원하면 총 $O(n^2)$ 이다. 고유벡터까지 누적하면 $O(n^3)$ 으로 뛴다. 작은 $n$ 의 기본값.
• 분할 정복(Cuppen 1981, dstedc). $T$ 를 반으로 갈라 $T = \mathrm{diag}(T_1,T_2) + \rho\,vv^{\top}$ 로 쓰고, 두 부분 문제를 재귀로 푼 뒤 랭크-1 수정을 세큘러 방정식 $1 + \rho\sum_i \zeta_i^2/(d_i - \lambda) = 0$ 을 유리함수 반복으로 풀어 병합한다. 고유값이 서로 가깝거나 $\zeta_i$ 가 작으면 그 성분은 계산 없이 통과시키는 디플레이션이 대량으로 일어나 실측이 최악 경우보다 훨씬 빠르다.⁴ 고유벡터 전부가 필요한 큰 $n$ 에서 사실상 최강이며 dsyevd 의 기본 경로다. 자세한 내용은 분할 정복 고유값 알고리즘 참고.
• MRRR(Dhillon–Parlett, dstemr/dsyevr). 뭉친 고유값 무리마다 시프트를 달리한 $LDL^{\top}$ 표현(relatively robust representation)을 따로 두고 뒤틀린 인수분해로 고유벡터를 뽑는다. 재직교화 없이 직교 고유벡터가 나오는 것이 핵심으로, 총 비용이 $O(n^2)$ — 유일한 $O(n^2)$ 선택지다.
• 이분법 + 역반복법(dstebz+dstein). $T - \sigma I = LDL^{\top}$ 에서 $D$ 의 음수 성분 개수가 $\sigma$ 보다 작은 고유값의 개수와 같다(스텀 수열 성질). 이 계수 함수만 있으면 이분법으로 임의의 구간·임의의 번호대 고유값만 골라 완전한 정밀도로 격리할 수 있다. “500번째부터 520번째까지만 필요하다”는 요구에 답하는 유일한 방법. 약점은 뭉친 고유값의 고유벡터가 직교하지 않아 그람-슈미트 재직교화가 필요하다는 것이고, 그 결함을 정면으로 없애려고 나온 것이 위의 MRRR이다.

6. 반복법 대응물과 비대칭 대비

$A$ 가 거대한 희소행렬이면 $O(n^3)$ 축약 자체가 불가능하다. 이때의 대응물이 란초스 알고리즘이다. 3항 점화식으로 크릴로프 부분공간의 정규직교기저를 만들면 그 위에서 $A$ 의 작용이 자동으로 삼중대각 $T_m$ 이 된다 — 행렬-벡터 곱만으로 부분적인 삼중대각화를 수행하는 셈이다. 대가는 유한 정밀도에서의 직교성 붕괴와 그로 인한 유령 고유값이며, 하우스홀더 축약이 무조건 안정한 것과 정확히 대비된다. 유한 단계·완전 안정·$O(n^3)$ 대 무한 반복·불안정·행렬-벡터 곱, 이것이 밀집과 희소의 갈림길이다.

비대칭 행렬로 가면 같은 도구로 한 칸 덜 간다. $Q^{\top}AQ = H$ 상헤센베르크가 한계이고 비용은 $\tfrac{10}{3}n^3$($Q$ 누적 시 $\tfrac{14}{3}n^3$)이다. 대칭성이 없으면 왼쪽에서 지운 성분을 오른쪽 곱이 되살리는 범위가 넓어지기 때문. 여기서 QR 알고리즘을 돌려 슈어 분해로 간다. 직교성을 포기하면 비대칭 행렬도 삼중대각으로 보낼 수 있지만(양측 란초스), 변환행렬의 조건수에 상한이 없고 붕괴(breakdown)가 실제로 일어난다 — 안정성을 팔아 구조를 사는 거래다.

두 가지 이웃도 짚어 두자. 대칭 정부호 일반화 고유값 문제 $K\phi = \lambda M\phi$ 는 $M = LL^{\top}$ 로 촐레스키 분해한 뒤 $C = L^{-1}KL^{-\top}$ 를 삼중대각화하는 것이 표준 경로이며(모드 해석의 기본기), 특이값 분해 쪽의 짝은 좌우에서 서로 다른 반사를 거는 쌍대각화다.

7. 관련 문서

• 하우스홀더 변환 · 기븐스 회전 · 그람-슈미트 · QR 분해
• 고유값 문제 · 슈어 분해 · 특이값 분해 · 역반복법
• 란초스 알고리즘 · 아놀디 알고리즘 · 크리로프 부분공간법
• LAPACK · 조건수 · 후진 오차 해석 · 모드 해석

8. Footnotes

1. $n-2$ 인 이유는 마지막 열에는 지울 것이 남아 있지 않기 때문이다. $n-1$ 로 잡아도 마지막 반사가 아무것도 하지 않으니 답은 맞지만, 이 오프바이원은 for 문의 상한과 반사 벡터 저장 위치를 동시에 어긋나게 만드는 재주가 있어서 축약 코드를 손으로 짜는 사람의 첫 번째 통과의례로 유명하다. ↩

2. 그래서 “고유값 solver의 반복 횟수를 0으로 만들 수 있느냐”는 질문은 수치해석 문제가 아니라 갈루아 이론 문제다. 반대로 $n \le 4$ 에서는 근의 공식이 존재하므로 원리적으로 유한 알고리즘이 가능한데, 아무도 그걸 쓰지 않는다. 4차 근의 공식은 수치적으로 재앙이라서 그냥 반복 돌리는 편이 정확하다. ↩

3. 고유값 계산의 병목이 “고유값을 구하는 부분”이 아니라 “구하기 전 준비운동”이라는 사실은 처음 프로파일러를 켜 본 사람을 대부분 당황시킨다. dsyev 를 프로파일링하면 시간의 절반 이상이 dsytrd 안에 있고, 그 안의 절반은 dsymv 안에 있다. 메모리 대역폭이 유죄. ↩

4. 이 디플레이션 비율이 행렬마다 달라서, 분할 정복의 실측 복잡도를 논문에 쓸 때 “최악 $O(n^3)$, 실제로는 대체로 그보다 한참 아래”라는 애매한 문장이 등장한다. 성능이 데이터에 의존한다는 뜻이고, 벤치마크 행렬을 고르는 사람이 결론을 고를 수 있다는 뜻이기도 하다. ↩