QR 알고리즘

1. 개요[편집]

QR 알고리즘 QR algorithm
고안	J. G. F. Francis (1961–62), V. N. Kublanovskaya (1961)
선조	Rutishauser의 LR 알고리즘 (1955)
기본 반복	$A_k = Q_k R_k,\ \ A_{k+1} = R_k Q_k$
실전 구성	헤센베르크 축약 + 암시적 이중 이동 + 디플레이션
비용	고유값만 $\approx 10n^3$, 슈어 벡터까지 $\approx 25n^3$
LAPACK	dgehrd → dhseqr (드라이버 dgeev)

분해하고, 순서를 뒤집어 다시 곱한다. 이 어처구니없이 단순한 반복이 세상의 고유값을 전부 책임진다.

QR 알고리즘(QR algorithm)은 정사각 행렬을 QR 분해한 뒤 두 인수를 반대 순서로 다시 곱하는 조작을 반복해 고유값을 뽑아내는 반복법이다. 1961~62년 존 프랜시스(J. G. F. Francis)와 베라 쿠블라노프스카야(V. N. Kublanovskaya)가 독립적으로 제시했고, 60여 년이 지난 지금도 조밀 비대칭 행렬의 전체 고유값을 구하는 사실상 유일한 표준이다. numpy.linalg.eig 한 줄 뒤에서 도는 것이 이것이다.

한 줄짜리 정의는 순진하지만 실제 구현은 그렇지 않다. 실무 QR 알고리즘은 (1) 헤센베르크 축약으로 반복당 비용을 $O(n^3)$ 에서 $O(n^2)$ 로 내리고, (2) 이동(shift)으로 수렴을 선형에서 2차 이상으로 끌어올리고, (3) 암시적 이중 이동으로 실수 산술만 쓰면서 복소 고유값을 처리하고, (4) 디플레이션으로 문제 크기를 줄여 나가는 네 겹의 장치로 굴러간다. 이 문서는 그 네 겹을 차례로 벗긴다. 결과물인 슈어 분해 자체의 성질과 응용은 그쪽 문서가 다룬다.¹

2. 순수 QR 반복의 정체 — 부분공간 반복[편집]

이동 없는 순수 반복

$$A_k = Q_k R_k, \qquad A_{k+1} = R_k Q_k = Q_k^{\mathsf H} A_k Q_k$$

는 매 스텝이 유니터리 상사변환이므로 고유값이 보존된다. 문제는 왜 이게 상삼각형으로 수렴하느냐인데, 답은 이 반복이 위장한 거듭제곱법이라는 데 있다.

누적곱 $\hat{Q}_k = Q_1 Q_2 \cdots Q_k$, $\hat{R}_k = R_k \cdots R_2 R_1$ 을 두면 귀납법 한 줄로

$$A^k = \hat{Q}_k \hat{R}_k$$

가 나온다.² 즉 $\hat{Q}_k$ 는 $A^k$ 의 QR 분해에서 나오는 직교 인수다. 이것이 전부다.

• $\hat{Q}_k$ 의 첫 열은 $A^k e_1$ 을 정규화한 것 — 정확히 거듭제곱법이고, 지배 고유벡터로 $|\lambda_2/\lambda_1|^k$ 의 비율로 수렴한다.
• 앞쪽 $j$ 개 열이 펼치는 공간은 $A^k$ 의 앞 $j$ 열이 펼치는 공간 — 즉 직교(부분공간) 반복이며, 지배 $j$ 차원 불변 부분공간으로 $|\lambda_{j+1}/\lambda_j|^k$ 로 수렴한다.

불변 부분공간으로 수렴한다는 말은 $A_{k+1} = \hat{Q}_k^{\mathsf H} A \hat{Q}_k$ 의 좌하단 블록이 0으로 죽는다는 뜻이고, 모든 $j$ 에 대해 이것이 일어나면 상삼각형이 된다. 부대각 성분 $|a^{(k)}_{j+1,j}|$ 는

$$\left|a^{(k)}_{j+1,j}\right| = O\!\left(\left|\frac{\lambda_{j+1}}{\lambda_j}\right|^{k}\right)$$

로 줄어든다. 여기서 세 가지가 즉시 따라 나온다. 고유값의 절댓값이 서로 달라야 하고(같은 크기의 켤레쌍이 있으면 그 자리는 영원히 안 죽는다), 수렴은 선형이며, 비율이 1에 가까우면 실용적으로 멈춘 것과 같다. 뒤의 이동은 정확히 이 세 결함을 겨냥한 장치다.

동시에 반대 방향의 해석도 성립한다. $\hat{Q}_k$ 의 마지막 열은 $A^{-\mathsf{H}}$ 에 대한 역반복법을 $e_n$ 에서 출발해 $k$ 번 돌린 것과 같다. 즉 순수 QR 반복은 위에서 거듭제곱법, 아래에서 역반복을 동시에 굴리고 있다 — 이동이 왜 아래쪽 모서리를 먼저 죽이는지가 여기서 설명된다.

3. 헤센베르크 축약 — 반복당 $O(n^3) \to O(n^2)$[편집]

밀집 행렬에 QR 분해를 한 번 하는 것만 $\tfrac{4}{3}n^3$ flops다. 수십 번 반복하면 답이 없다. 그래서 반복에 들어가기 전에 하우스홀더 변환 $n-2$ 번으로 상헤센베르크 형($i > j+1$ 에서 $h_{ij}=0$, 부대각 바로 아래까지만 채워진 형태)으로 만든다. 유한 단계로 끝나고 $\tfrac{10}{3}n^3$ flops, 직교 인수까지 누적하면 $\tfrac{14}{3}n^3$ 이다.

이 축약이 값을 하는 이유는 헤센베르크 구조가 QR 스텝에서 보존되기 때문이다. $H = QR$ 에서 $H$ 가 헤센베르크면 $Q$ 도 헤센베르크이고, 상삼각 $R$ 을 곱한 $RQ$ 도 헤센베르크다. 그래서 한 번 축약해 두면 이후 모든 반복이 그 안에서 논다. 게다가 헤센베르크 행렬의 QR 분해는 죽일 성분이 부대각 $n-1$ 개뿐이라 기븐스 회전 $n-1$ 개면 끝나므로, 반복 하나가 $O(n^2)$ 으로 떨어진다. 대칭 행렬이면 헤센베르크가 곧 삼중대각이므로 삼중대각화가 같은 역할을 하고, 반복 하나는 $O(n)$ 이 된다.

축약 후에는 반드시 기약(unreduced, 모든 부대각이 0이 아님) 여부를 본다. 부대각이 0이면 그 자리에서 문제가 두 개로 쪼개진다 — 이것이 아래의 디플레이션과 같은 장치다.

4. 이동 — 수렴 차수를 올리는 장치[편집]

이동을 넣은 반복은

$$H_k - \mu_k I = Q_k R_k, \qquad H_{k+1} = R_k Q_k + \mu_k I$$

이고 여전히 상사변환이다. 앞의 부분공간 해석에서 $A$ 가 $A - \mu I$ 로 바뀌었으므로, 아래쪽 부대각의 수렴비는

$$\left|\frac{\lambda_n - \mu}{\lambda_{n-1} - \mu}\right|$$

가 된다. $\mu$ 를 $\lambda_n$ 바로 옆에 갖다 놓으면 이 비가 0에 가까워지고, 매 스텝 $\mu$ 를 갱신하면 비율 자체가 반복마다 줄어들어 선형 수렴이 초선형으로 바뀐다. 이동 QR이 사실상 이동 역반복이라는 앞의 관찰이 여기서 값을 한다.

• 레일리 몫 이동 $\mu = h_{nn}$: 아래쪽 모서리 성분을 그대로 이동으로 쓴다. $h_{nn}$ 이 $e_n$ 방향의 레일리 몫이므로 이 반복은 정확히 레일리 몫 반복(RQI) 이다. 따라서 일반적으로 2차 수렴, 대칭(에르미트)이면 3차 수렴한다.³ 흔히 “이동이 3차 수렴을 만든다”고 뭉뚱그리지만, 3차는 대칭인 경우의 이야기이고 비대칭 일반에서는 2차다.
• 윌킨슨 이동: 우하단 $2\times2$ 블록의 고유값 중 $h_{nn}$ 에 가까운 쪽을 쓴다. $h_{nn}$ 을 그냥 쓰면 실수 산술에서 켤레쌍 앞에서 갇히는데, $2\times2$ 를 보면 켤레쌍을 알아채므로 그 함정을 피한다. 대칭 삼중대각에 대해서는 윌킨슨 이동 QR의 전역 수렴이 증명돼 있고, 실측으로 고유값 하나당 평균 두세 번의 반복이면 끝난다.

주의할 것은 비대칭 일반에서는 전역 수렴이 증명돼 있지 않다는 점이다. 배터슨·스밀리 등이 이동 QR이 갇히는 예를 만들었고, 그래서 실무 코드는 안전장치를 둔다. LAPACK의 dlahqr 은 같은 부분문제에서 10회 반복이 지나면 예외 이동(exceptional shift)을 한 번 끼워 넣어 대칭성을 깨고 빠져나온다. “수렴은 신에게 맡긴다”가 아니라 “10번 봐주고 발로 찬다”에 가깝다.

5. 암시적 이중 이동과 실 슈어 형[편집]

실행렬의 고유값은 켤레복소쌍으로 나올 수 있는데, 좋은 이동은 그 복소값이다. 그러면 복소 산술을 써야 하고 메모리·연산이 두 배가 된다. 프랜시스의 해법이 암시적 이중 이동(implicit double shift)이다.

켤레쌍 $\mu, \bar\mu$ 를 한꺼번에 적용하는 두 스텝을 합치면

$$M = (H - \mu I)(H - \bar\mu I) = H^2 - 2\,\mathrm{Re}(\mu)\,H + |\mu|^2 I$$

이 되고, 계수가 전부 실수이므로 $M$ 은 실행렬이다. 그런데 $M$ 을 통째로 만드는 것은 $O(n^3)$ 이라 그것도 못 쓴다. 여기서 암시적 Q 정리가 구원한다 — 기약 헤센베르크로 축약하는 직교행렬은 첫 열이 정해지면 (열 부호를 빼고) 유일하다. 그러므로 $M$ 의 첫 열만 계산해(헤센베르크라서 성분 세 개뿐이다) 그것을 $e_1$ 로 보내는 작은 하우스홀더 반사를 때린 뒤, 망가진 헤센베르크 구조를 복원하기만 하면 결과는 $M$ 을 명시적으로 QR 분해한 것과 같다.

구조 복원 과정이 유명한 부풀림 쫓기(bulge chasing)다. 첫 반사가 좌상단에 $3\times3$ 짜리 혹을 만들고, 이어지는 작은 반사들이 그 혹을 대각선을 따라 우하단으로 밀어내 결국 행렬 밖으로 떨어뜨린다. 스윕 하나가 $O(n^2)$, 명시적 QR 분해는 한 번도 하지 않는다(그래서 “암시적”이다).

이 방식의 산출물이 실 슈어 형이다. 실수 산술만 썼으니 복소 고유값을 대각에 놓을 수 없고, 대신 켤레쌍마다 $2\times2$ 대각 블록이 남는 준상삼각 행렬 $T$ 가 나온다. 각 $2\times2$ 블록의 고유값을 근의 공식으로 읽으면 켤레쌍이 나온다.

6. 디플레이션과 전체 비용[편집]

부대각 성분 하나가 충분히 작아지면 — 실무 판정은 절대 기준이 아니라

$$|h_{k+1,k}| \le u\,\big(|h_{kk}| + |h_{k+1,k+1}|\big)$$

같은 상대 기준이다($u$ 는 단위 반올림) — 그것을 0으로 놓고 문제를 두 조각으로 쪼갠다. 보통은 우하단 $1\times1$ 또는 $2\times2$ 가 떨어져 나가고, 활성 창이 한두 칸씩 줄면서 반복이 계속된다. 이 디플레이션이 없으면 이미 수렴한 고유값에 계속 연산을 낭비하게 된다.

현대 구현은 여기에 공격적 조기 디플레이션(Braman–Byers–Mathias, 2002)을 얹는다. 우하단 창을 슈어 형으로 따로 풀어 보고 “사실상 수렴한” 고유값들을 한꺼번에 떼어내는 방식인데, 전통적 판정보다 훨씬 이르게 디플레이션이 일어나 큰 행렬에서 수 배의 속도 차이를 낸다. 여기에 다중 이동(한 스윕에 이동 수십 개를 동시에 밀어 넣어 BLAS-3 블록 연산으로 바꾸는 것)을 결합한 것이 dhseqr 의 현재 모습이다.

전체 비용은 고유값만 필요하면 $\approx 10n^3$ flops, 슈어 벡터까지 누적하면 $\approx 25n^3$ 이다. 반복법인데도 비용을 $n^3$ 의 상수배로 적을 수 있는 이유는, 이동이 워낙 잘 들어서 고유값 하나당 평균 반복 수가 사실상 상수이기 때문이다. 원리적으로는 유한 단계 알고리즘이 존재할 수 없는 문제인데도(아벨-루피니) 실무에서는 직접법처럼 굴러가는 것이 QR 알고리즘의 기묘한 지위다.

7. 대칭이면 길이 갈린다[편집]

행렬이 대칭이면 삼중대각화 후의 선택지가 여러 갈래로 벌어지고, QR은 더 이상 유일한 답이 아니다.

방법	LAPACK	고유값만	전체 고유쌍	성격
대칭 QR (윌킨슨 이동)	dsteqr	$O(n^2)$	$O(n^3)$	가장 견고, 전역 수렴 증명됨
분할 정복	dstedc	—	$O(n^3)$, 실측 최속	$O(n^2)$ 작업공간 필요
MRRR	dstemr	$O(n^2)$	$O(nk)$	재직교화 없이 직교 고유벡터
이분법 + 역반복	dstebz/dstein	부분집합	군집에서 재직교화	스텀 수열 기반
야코비 회전	—	느림	느림	상대 정확도 최고

갈림길의 기준은 고유벡터가 필요한가다. 삼중대각 QR은 고유값만 뽑으면 $O(n^2)$ 로 압도적이지만, 고유벡터를 얻으려면 기븐스 회전 전체를 누적해야 해서 $O(n^3)$ 이 되고 여기서 분할 정복에 진다. 실제로 dstedc 는 부분문제 크기가 SMLSIZ(기본 25) 아래로 내려가면 스스로 dsteqr 을 부른다 — 작은 문제에서는 QR이 여전히 이긴다는 뜻이다. MRRR(Multiple Relatively Robust Representations, 딜런·파를렛)은 잘 고른 $LDL^{\mathsf T}$ 표현들 위에서 고유벡터를 뽑아 재직교화를 아예 생략하며, 고유쌍 $k$ 개만 필요할 때 $O(nk)$ 라는 최적 비용을 낸다. 다만 극단적으로 뭉친 군집에서 직교성이 요구 수준에 못 미치는 사례가 보고돼 있어, LAPACK은 지금도 세 경로를 모두 유지한다.

비대칭 쪽에는 이런 사치가 없다. 비대칭 조밀 고유값 문제에는 QR 알고리즘 말고 실용적 경쟁자가 없다. 대안이 등장하려면 QDWH 기반 스펙트럼 분할 정복(행렬 부호 함수 참조)처럼 통신량을 줄이는 완전히 다른 노선으로 가야 한다.

8. 여담[편집]

• 프랜시스는 QR 알고리즘 논문 두 편을 낸 뒤 수치해석계를 떠나 산업계로 갔고, 반세기 동안 자기 알고리즘이 20세기 10대 알고리즘에 뽑힌 것도 모르고 살았다. 2007년에야 진 골룹 등이 그를 찾아냈고, 본인은 자신이 유명한 줄 전혀 몰랐다고 한다.⁴
• 이름은 QR 분해에서 왔지만 실전 구현은 QR 분해를 명시적으로 하지 않는다. 부풀림 쫓기가 그 자리를 대신한다. “QR 알고리즘의 코드에서 QR 분해를 찾아보라”는 것이 수치해석 수업의 고전 낚시다.
• 이동을 껐다 켜 보면 이 알고리즘의 절반이 이동이라는 게 실감 난다. 같은 행렬에서 이동 없이 40~50 스윕 걸리던 것이 윌킨슨 이동 다섯 스윕에 같은 자리에 닿는다.

9. 관련 문서[편집]

• QR 분해 · 슈어 분해 · 고유값 문제
• 하우스홀더 변환 · 기븐스 회전 · 삼중대각화
• 거듭제곱법 · 역반복법 · 레일리 몫
• 분할 정복 고유값 알고리즘 · 스텀 수열 · 야코비 고유값 알고리즘
• 후진 오차 해석 · 조건수 · 비정규 행렬
• LAPACK · 크릴로프-슈어 방법 · 특이값 분해

10. Footnotes[편집]

1. 정확히 말하면 QR 알고리즘은 “고유값을 구하는 알고리즘”이 아니라 “슈어 분해를 계산하는 알고리즘”이고, 고유값은 그 부산물로 대각에서 읽는 것이다. 고유벡터가 필요하면 슈어 형 $T$ 에 대해 삼각계를 풀어(dtrevc) 되돌리는 별도 단계가 붙는다. ↩

2. $A^k = A\,A^{k-1} = A\hat{Q}_{k-1}\hat{R}_{k-1} = \hat{Q}_{k-1}A_k\hat{R}_{k-1} = \hat{Q}_{k-1}Q_kR_k\hat{R}_{k-1} = \hat{Q}_k\hat{R}_k$. 세 번째 등호에 $A_k = \hat{Q}_{k-1}^{\mathsf H}A\hat{Q}_{k-1}$ 을 썼다. 두 줄짜리 귀납법이 알고리즘 전체의 수렴 이론을 낳는다. ↩

3. 레일리 몫 문서에서 보듯 대칭 행렬에서는 고유벡터 오차 $\varepsilon$ 이 고유값 오차 $\mathcal{O}(\varepsilon^2)$ 로 이어진다. 이 “공짜 한 자릿수”가 역반복의 2차 수렴에 곱해져 3차가 된다. 비대칭에서는 좌·우 고유벡터가 달라 이 공짜 자릿수가 없고, 그래서 2차에 머문다. ↩

4. Golub & Uhlig, “The QR algorithm: 50 years later” (IMA J. Numer. Anal., 2009)에 그 추적기가 실려 있다. 참고로 같은 시기 독립적으로 같은 결과에 도달한 쿠블라노프스카야는 소련 쪽 문헌 장벽 때문에 서방에서 한참 늦게 인정받았다. 알고리즘 하나에 냉전이 두 번 얽힌 셈. ↩