QZ 알고리즘

1. 개요

QZ 알고리즘 QZ algorithm
푸는 문제	일반화 고유값 문제 $Ax=\lambda Bx$ (밀집)
결과	일반화 슈어 형 $Q^HAZ=T$, $Q^HBZ=S$ (상삼각)
고유값	쌍 $(t_{ii}, s_{ii})$, $\lambda_i=t_{ii}/s_{ii}$ — $s_{ii}=0$ 이면 무한대
제안	모울러 · 스튜어트 (1973)
비용	$O(n^3)$, 같은 크기 QR 알고리즘의 대략 세 배
구현	LAPACK xGGHRD · xHGEQZ · xGGES · xGGEV

$B^{-1}A$ 를 만들지 않는 것, 그게 알고리즘의 전부다.

QZ 알고리즘(QZ algorithm)은 밀집 행렬의 일반화 고유값 문제 $Ax = \lambda Bx$ 를 푸는 표준 알고리즘이다. 두 유니터리(실수면 직교) 행렬 $Q$, $Z$ 를 반복적으로 갱신해 $A$ 와 $B$ 를 동시에 상삼각 형태로 만든다:

$$Q^H A Z = T, \qquad Q^H B Z = S.$$

이것이 슈어 분해 문서가 한 단락으로 소개한 일반화 슈어 분해이고, 고유값은 대각 원소 쌍의 비 $\lambda_i = t_{ii}/s_{ii}$ 로 읽는다. 모울러와 스튜어트가 1973년에 발표했다. 이름의 Q 와 Z 는 왼쪽·오른쪽 변환 행렬 이름에서 왔다 — QR 알고리즘의 R 처럼 무슨 분해를 뜻하지는 않는다.

핵심 착상은 **「$B^{-1}A$ 에 QR 알고리즘을 돌리는 것과 수학적으로 같은 일을, $B^{-1}$ 을 한 번도 만들지 않고 한다」**는 것이다. $B$ 가 특이하거나 조건이 나쁘면 $B^{-1}A$ 를 만드는 순간 정보가 날아가는데, QZ 는 유니터리 변환만 쓰므로 후진 안정하고, $B$ 가 정확히 특이해도 아무 문제 없이 「무한대 고유값」을 내놓는다. MATLAB 의 eig(A,B)(일반 행렬), SciPy 의 scipy.linalg.eig(a, b), Julia 의 eigen(A, B) 가 모두 LAPACK 의 QZ 루틴으로 내려간다.

2. 왜 $B^{-1}A$ 를 만들면 안 되나

$B$ 가 가역이면 $Ax = \lambda Bx$ 는 $B^{-1}Ax = \lambda x$ 와 같다. 그러니 $B^{-1}A$ 를 만들어 보통 고유값 문제로 풀면 되지 않나? 실험해 보면 답이 나온다. $6\times6$ 무작위 $A$ 와, 특이값이 $(1,1,1,1,1,10^{-13})$ 인 무작위 $B$(조건수 $10^{13}$)로 펜슬을 20개 만들고, 50자리 기준해와 비교해 크기가 보통인 고유값 5개의 최대 상대오차를 쟀다:

방법	중앙값	최선	최악
eigvals(solve(B, A))	$1.6\times10^{-3}$	$2.0\times10^{-4}$	$1.7\times10^{-2}$
QZ (scipy.linalg.eigvals(A, B))	$1.3\times10^{-15}$	$4.2\times10^{-16}$	$2.3\times10^{-14}$

나머지 하나의 고유값은 $10^{12}$ 근처로 튀어 나가는데, 이건 $B$ 가 「거의 특이」하다는 정보를 정직하게 반영한 것이다. 문제는 그 거대한 고유값이 아니라, 멀쩡한 다섯 개까지 12자리를 잃는다는 데 있다. $B^{-1}A$ 를 만드는 단계에서 $\kappa(B)\,u \approx 10^{-3}$ 크기의 오차가 행렬 전체에 퍼지고, 그 뒤 아무리 좋은 고유값 솔버를 써도 되돌릴 수 없다. QZ 는 원래 펜슬 $(A, B)$ 에 대해 후진 안정이라 멀쩡한 고유값은 멀쩡하게 나온다.

$B$ 가 대칭 양정치라면 촐레스키 분해로 $L^{-1}AL^{-T}$ 를 만드는 길이 있고, 대칭 구조를 지킬 수 있어 더 빠르다(LAPACK xSYGV). 하지만 이것도 $B$ 의 조건수가 나쁘면 같은 종류의 손실을 겪는다. 비대칭이거나 $B$ 가 특이할 수 있으면 QZ 가 기본이다.

3. 고유값은 수가 아니라 쌍이다

QZ 의 출력은 $\lambda_i$ 가 아니라 쌍 $(\alpha_i, \beta_i) = (t_{ii}, s_{ii})$ 다. LAPACK xGGEV 도 ALPHAR, ALPHAI, BETA 를 따로 돌려준다. 이 표현 덕분에 세 가지 상황을 같은 틀로 다룬다.

• $\beta_i \ne 0$ — 보통의 유한 고유값 $\alpha_i/\beta_i$.
• $\beta_i = 0$, $\alpha_i \ne 0$ — 무한대 고유값. $B$ 가 특이하다는 뜻이다. 앞의 실험에서 $B$ 의 마지막 특이값을 정확히 0 으로 바꾸면 QZ 는 $s_{66} = 0$ 을 내놓고 SciPy 는 그 고유값을 inf 로 표시한다. 나머지 다섯 개는 아무 일 없다는 듯 그대로 나온다.
• $\alpha_i = \beta_i = 0$ — 특이 펜슬. 모든 $\lambda$ 에 대해 $\det(A - \lambda B) = 0$ 이라 고유값이 정의되지 않는다. 이때 다른 쌍들의 비율은 반올림에 따라 아무 값이나 될 수 있으므로, 실무에서는 이 경우를 먼저 감지해 구조(크로네커 표준형)를 따로 분석해야 한다.

쌍 표현은 $\lambda$ 와 $1/\lambda$ 를 대칭적으로 다룬다는 장점도 있다. $Ax = \lambda Bx$ 와 $Bx = \mu Ax$ 는 $\mu = 1/\lambda$ 인 같은 문제이고, 쌍 $(\alpha, \beta)$ 는 그냥 $(\beta, \alpha)$ 로 뒤집힐 뿐이다. 그래서 정확도도 「상대오차」 대신 현 거리(chordal metric) $\lvert\alpha\hat\beta - \beta\hat\alpha\rvert / \big(\sqrt{\lvert\alpha\rvert^2+\lvert\beta\rvert^2}\sqrt{\lvert\hat\alpha\rvert^2+\lvert\hat\beta\rvert^2}\big)$ 로 재는 것이 스튜어트 이래의 관례다 — 리만 구 위의 거리라 무한대 근처의 고유값도 공평하게 잴 수 있다.

4. 알고리즘 — 세 단계

1단계: 헤센베르크-삼각 축약. 먼저 $B$ 를 QR 분해해 $Q_0^H B$ 를 상삼각으로 만들고 같은 $Q_0^H$ 를 $A$ 에도 곱한다. 그다음 $A$ 를 상 헤센베르크 형태로 만드는데, 왼쪽에서 기븐스 회전으로 $A$ 의 아래쪽 원소를 지우면 $B$ 의 삼각 구조에 부채꼴 원소가 하나 생기고, 그것을 오른쪽 기븐스 회전으로 다시 지운다. 이 「지우고 → 생기고 → 지우는」 짝을 열마다 반복하면 $O(n^3)$ 에 $(A, B)$ 가 (헤센베르크, 상삼각) 쌍이 된다(LAPACK xGGHRD). 보통 고유값 문제의 헤센베르크 축약과 같은 역할이다.

2단계: 암시적 QZ 스텝. 이제 $AB^{-1}$ 이 헤센베르크라는 사실(헤센베르크 × 상삼각의 역 = 헤센베르크)을 이용해, $AB^{-1}$ 에 대한 암시적 이중 시프트 QR 스텝을 $A$ 와 $B$ 위에서 직접 흉내 낸다. 시프트는 $AB^{-1}$ 오른쪽 아래 $2\times2$ 블록의 고유값으로 고르고, 그 블록에 필요한 성분은 $A$ 와 $B$ 의 마지막 몇 행·열만으로 계산된다. 첫 열에 하우스홀더 반사를 걸어 「불룩한 부분(bulge)」을 만든 뒤, 왼쪽 반사로 $A$ 의 헤센베르크 구조를, 오른쪽 반사로 $B$ 의 삼각 구조를 번갈아 복구하며 불룩함을 대각선 아래로 쫓아낸다(bulge chasing). 실수 행렬이면 켤레 복소 고유값 쌍이 $2\times2$ 블록으로 남아 준삼각 형태가 된다 — 실 슈어 형과 똑같은 사정이다.

3단계: 축소(deflation). $A$ 의 부대각 원소가 충분히 작아지면 문제를 둘로 쪼갠다. 여기에 QZ 고유의 경우가 하나 더 있다. $B$ 의 대각 원소가 0(에 가까움)이 되면 무한대 고유값이 있다는 신호이고, 모울러·스튜어트는 회전으로 그 0 을 대각선 맨 아래(또는 위)까지 밀어 내려 무한대 고유값을 통째로 떼어 내는 절차를 제시했다. 이 처리 덕분에 $B$ 가 정확히 특이해도 반복이 멈추거나 나눗셈이 터지지 않는다.

비용은 골럽·반 론 기준으로 고유값만 구하면 약 $30n^3$ flops, $Q$·$Z$ 까지 쌓으면 약 $66n^3$ 로, 같은 크기 비대칭 QR 알고리즘(약 $10n^3$ / $25n^3$)의 대략 세 배다. 현대 LAPACK 은 여러 시프트를 한꺼번에 쫓는 다중 시프트와 반복 도중 이미 수렴한 부분을 일찍 떼어 내는 공격적 조기 축소(aggressive early deflation)를 QZ 에도 적용해 큰 $n$ 에서 몇 배 빨라졌다(코그스트룀·크레스너 2006).¹

5. 후진 안정성과 정렬

QZ 가 계산한 $(T, S)$ 는 어떤 근처 펜슬 $(A + E, B + F)$, $\lVert E\rVert \approx u\lVert A\rVert$, $\lVert F\rVert \approx u\lVert B\rVert$ 의 정확한 일반화 슈어 형이다(후진 오차 해석). 그러니 고유값의 실제 오차는 고유값 조건수의 펜슬판 — 왼·오른 고유벡터 $y$, $x$ 로 $\lVert x\rVert\lVert y\rVert\sqrt{\lVert A\rVert^2 + \lVert B\rVert^2}\big/\sqrt{\lvert y^HAx\rvert^2 + \lvert y^HBx\rvert^2}$ 꼴 — 이 정한다. 「QZ 결과가 이상하다」는 보고의 대부분은 알고리즘이 아니라 이 조건수 문제다.

슈어 형의 대각 순서는 마음대로 바꿀 수 있다(LAPACK xTGSEN, xGGES 의 SELCTG 인자). 원하는 고유값(예: 안정한 것, 단위원 안쪽인 것)을 왼쪽 위로 모으면 $Z$ 의 앞 열들이 그 고유값들의 변형 부분공간(deflating subspace)의 직교 기저가 된다. 고유벡터 하나하나를 구하는 것보다 이 부분공간이 훨씬 잘 조건화돼 있어서, 응용 대부분은 고유벡터가 아니라 정렬된 슈어 형을 쓴다.

6. 어디서 쓰나

• 리카티 방정식. 반 두런(1981)은 연속·이산 대수 리카티 방정식을 확장 펜슬의 정렬된 QZ 로 풀었다. 이산시간에서는 $A$ 의 역행렬이 필요 없어지고, 펜슬이 무한대 고유값을 가질 수 있는 경우($A$ 가 특이)도 그대로 처리된다. 제어 툴박스의 리카티 솔버들이 이 「확장 펜슬 + 정렬된 일반화 슈어 분해」 계열을 쓴다.
• 미분대수방정식과 기술자 시스템. $E\dot x = Ax$ 에서 $E$ 가 특이하면 펜슬 $(A, E)$ 의 무한대 고유값이 대수적 구속에, 유한 고유값이 진짜 동역학에 대응한다. 회로 시뮬레이션(MNA)과 다물체 동역학의 선형화 안정성 해석이 이 꼴이다.
• 구조 동역학. 질량이 없는 자유도(회전 자유도, 라그랑주 승수)가 있으면 질량행렬이 특이해 무한대 고유값이 생긴다. 작은 문제라면 QZ 로 그대로 풀고, 큰 문제는 시프트-역변환 크릴로프로 가되 그 무한대 고유값이 스펙트럼 변환 뒤 0 으로 가는 것을 감안한다. 모드 해석·강성행렬 참조.
• 2차 고유값 문제의 선형화. $(\lambda^2M + \lambda C + K)x = 0$ 을 크기 $2n$ 의 펜슬로 바꾸면 그 펜슬을 푸는 것이 QZ 다. 선형화 전에 계수 스케일링을 해야 한다는 이야기는 그 문서에 있다.
• 거시경제의 합리적 기대 모형. 선형화된 DSGE 모형을 푸는 클라인(2000)의 방법은 정렬된 일반화 슈어 분해로 안정·불안정 고유값을 가른다. 경제학 대학원생이 「QZ 분해」라는 단어를 처음 만나는 곳이 대개 여기다.²

7. 균형화와 스케일링

QZ 도 입력 스케일에 민감하다. 보통 고유값 문제에 행렬 균형화가 있듯, 펜슬에는 워드(1981)의 균형화가 있다 — 왼쪽·오른쪽 대각 행렬 $D_1$, $D_2$ 로 $(D_1AD_2, D_1BD_2)$ 의 원소 크기를 고르게 맞춘다(LAPACK xGGBAL). 다만 펜슬 균형화는 보통 균형화보다 효과가 들쭉날쭉해, LAPACK 드라이버들은 기본으로 「치환만」(고립된 고유값을 떼어 내는 것) 하는 쪽을 택한다. 단위계가 섞인 구조·회로 모델이라면 물리적으로 의미 있는 스케일링을 사용자가 먼저 하는 것이 낫다.

8. 관련 문서

• 슈어 분해 · QR 알고리즘 · 고유값 문제 · 2차 고유값 문제 · 비선형 고유값 문제
• 기븐스 회전 · 하우스홀더 변환 · 후진 오차 해석 · 고유값 조건수 · 행렬 균형화
• 리카티 방정식 · 미분대수방정식 · 크로네커 표준형 · LAPACK

9. Footnotes

1. 「다 끝난 줄 알았던 알고리즘을 30년 뒤에 몇 배 빠르게 만들기」는 수치선형대수의 국룰 장르다. QR 알고리즘도 2002년 브라이어스·바이어스·매슈스가 같은 짓을 했다. 원리를 바꾼 게 아니라 메모리 계층에 맞게 갈아 넣은 결과다. ↩

2. 그래서 경제학 쪽 문헌에서는 QZ 가 「블랜처드-칸 조건을 확인하는 그 분해」로 통한다. 수치해석 사람과 경제학 사람이 같은 LAPACK 루틴을 부르면서 서로 다른 이름으로 부르는 드문 사례. ↩