30.30 QR 알고리즘의 수렴성 분석과 시프트 전략

30.30 QR 알고리즘의 수렴성 분석과 시프트 전략

1. 기본 QR 알고리즘의 수렴성

기본 QR 반복 A_k = Q_k R_k, A_{k+1} = R_k Q_k에서 A_k의 하삼각 성분이 0에 수렴하는 속도를 분석한다.

1.1 수렴 정리

정리. n \times n 행렬 A의 고유값이 \lvert \lambda_1 \rvert > \lvert \lambda_2 \rvert > \cdots > \lvert \lambda_n \rvert > 0을 만족하고, A가 대각화 가능하면, 기본 QR 반복에서

(A_k)_{ij} = O\left(\left\lvert \frac{\lambda_i}{\lambda_j} \right\rvert^k\right), \quad i > j

이다. 즉, A_k(i, j) 하삼각 성분 (i > j)은 \lvert \lambda_i / \lambda_j \rvert^k의 속도로 0에 수렴한다.

특히, 부대각 성분 (A_k)_{i+1,i}의 수렴 속도는 \lvert \lambda_{i+1} / \lambda_i \rvert^k에 비례한다. 인접 고유값의 비율 \lvert \lambda_{i+1} / \lambda_i \rvert이 작을수록 해당 성분이 빠르게 소거되고, 1에 가까울수록 수렴이 느려진다.

1.2 수렴의 한계

다음의 경우 기본 QR 알고리즘의 수렴이 느리거나 실패한다:

  • 인접 고유값의 절댓값이 근접: \lvert \lambda_{i+1} \rvert \approx \lvert \lambda_i \rvert이면 수렴 비율이 1에 가까워 실용적 수렴에 과도한 반복이 필요하다.
  • 절댓값이 같은 고유값 쌍: \lvert \lambda_i \rvert = \lvert \lambda_{i+1} \rvert (예: 복소 켤레 쌍 또는 \lambda_i = -\lambda_{i+1})이면 해당 블록이 수렴하지 않는다.
  • 영 고유값: \lambda_n = 0이면 A가 특이이고 QR 분해가 정의되지 않을 수 있다.

이러한 한계를 극복하기 위하여 시프트 전략이 도입된다.

2. 시프트된 QR 알고리즘의 수렴성

시프트된 QR 반복 A_k - \sigma_k I = Q_k R_k, A_{k+1} = R_k Q_k + \sigma_k I에서 시프트 \sigma_k를 적절히 선택하면 수렴이 극적으로 가속된다.

2.1 수렴 메커니즘

시프트된 QR 반복에서 A_k(n, n-1) 성분의 수렴 속도는

(A_{k+1})_{n,n-1} = O\left(\left\lvert \frac{\lambda_n - \sigma_k}{\lambda_{n-1} - \sigma_k} \right\rvert\right) (A_k)_{n,n-1}

에 의하여 결정된다. \sigma_k\lambda_n에 가깝게 선택하면 \lvert \lambda_n - \sigma_k \rvert / \lvert \lambda_{n-1} - \sigma_k \rvert \ll 1이 되어 수렴이 가속된다.

이 메커니즘은 시프트된 역멱승법과 동일한 원리이다. 실제로, 시프트된 QR 반복의 마지막 행의 거동은 시프트 \sigma_k를 사용한 역멱승법과 수학적으로 동치이다.

3. 시프트 전략

3.1 단일 시프트 (Rayleigh Quotient Shift)

\sigma_k = (A_k)_{nn}

A_k(n, n) 성분을 시프트로 사용한다. 이 전략의 수렴 특성:

  • 대칭 행렬: 3차 수렴(cubic convergence). \lvert (A_{k+1})_{n,n-1} \rvert = O(\lvert (A_k)_{n,n-1} \rvert^3).
  • 비대칭 행렬: 2차 수렴(quadratic convergence).

대칭 행렬에서의 3차 수렴은 (A_k)_{nn}이 레일리 몫 반복(Rayleigh quotient iteration)의 시프트와 동일하기 때문이다.

그러나 Rayleigh quotient 시프트는 전역 수렴(global convergence)이 보장되지 않는 경우가 있다. 예외적인 초기 조건에서 수렴이 실패할 수 있다.

3.2 Wilkinson 시프트

A_k의 우하단 2 \times 2 블록

B = \begin{pmatrix} a_{n-1,n-1} & a_{n-1,n} \\ a_{n,n-1} & a_{nn} \end{pmatrix}

의 두 고유값 \mu_1, \mu_2a_{nn}에 더 가까운 것을 시프트로 선택한다:

\sigma_k = \mu_j, \quad \text{where } \lvert \mu_j - a_{nn} \rvert = \min(\lvert \mu_1 - a_{nn} \rvert, \lvert \mu_2 - a_{nn} \rvert)

B의 고유값은 이차 공식으로 계산된다:

\mu_{1,2} = \frac{(a_{n-1,n-1} + a_{nn}) \pm \sqrt{(a_{n-1,n-1} - a_{nn})^2 + 4 a_{n-1,n} a_{n,n-1}}}{2}

Wilkinson 시프트의 장점:

  • 대칭 행렬에서 전역 수렴이 보장된다. Wilkinson (1968)에 의하여 증명되었다.
  • 수렴 속도는 최소 2차(quadratic)이며, 대부분의 경우 3차에 가깝다.
  • Rayleigh quotient 시프트의 전역 수렴 실패 문제를 해결한다.

3.3 Francis의 이중 시프트 (Double Shift)

비대칭 실수 행렬이 복소 켤레 고유값 쌍 \alpha \pm i\beta를 가질 때, 단일 실수 시프트로는 해당 2 \times 2 블록을 디플레이트할 수 없다. Francis의 이중 시프트 전략은 두 개의 복소 켤레 시프트를 동시에 적용하되, 암묵적(implicit) 방식으로 실수 산술만을 사용한다.

\sigma_k = \alpha + i\beta\bar{\sigma}_k = \alpha - i\beta에 대하여 두 번의 시프트된 QR 단계를 결합하면

(A_k - \sigma_k I)(A_k - \bar{\sigma}_k I) = A_k^2 - 2\alpha A_k + (\alpha^2 + \beta^2) I

이 행렬은 실수이다. Francis의 알고리즘은 이 이중 단계를 **암묵적 QR 단계(implicit QR step)**로 구현하여, 복소 산술을 완전히 회피한다.

4. 디플레이션(Deflation) 전략

시프트된 QR 반복에서 부대각 성분 (A_k)_{m, m-1}이 충분히 작아지면 (\lvert (A_k)_{m, m-1} \rvert < \epsilon, 여기서 \epsilon은 수렴 판정 기준), 해당 성분을 0으로 설정하고 행렬을 분할한다.

하단 디플레이션. (A_k)_{n, n-1} \approx 0이면 (A_k)_{nn}이 고유값이며, 크기 (n-1) \times (n-1)의 부분 행렬에 대하여 QR 반복을 계속한다.

블록 디플레이션. (A_k)_{n-1, n-2} \approx 0이면 우하단 2 \times 2 블록의 두 고유값 (복소 켤레 쌍일 수 있음)을 추출하고, 크기 (n-2) \times (n-2)의 부분 행렬에 대하여 계속한다.

중간 디플레이션. 중간 위치의 부대각 성분 (A_k)_{m, m-1} \approx 0이면 행렬을 두 블록으로 분할하여 독립적으로 처리한다.

디플레이션에 의하여 문제의 크기가 점진적으로 줄어들며, 최종적으로 모든 고유값이 추출된다.

5. 수렴 판정 기준

부대각 성분의 수렴 판정 기준으로는 다음이 사용된다:

\lvert (A_k)_{i+1, i} \rvert \leq \epsilon_{\text{mach}} (\lvert (A_k)_{ii} \rvert + \lvert (A_k)_{i+1, i+1} \rvert)

여기서 \epsilon_{\text{mach}}는 기계 정밀도(machine epsilon)이다. 이 기준은 상대적 크기에 기반하여 고유값의 스케일에 적응적이다.

6. 시프트 전략에 따른 수렴 속도 비교

시프트 전략수렴 차수 (대칭)수렴 차수 (비대칭)전역 수렴비고
시프트 없음선형선형조건부r = \lvert \lambda_{i+1}/\lambda_i \rvert
Rayleigh quotient3차2차미보장예외적 실패 가능
Wilkinson2차 이상2차대칭에서 보장실용적 표준
Francis 이중 시프트2차실용적 보장비대칭 행렬의 표준

7. 총 계산 비용 분석

하이젠베르크 환원 후의 시프트된 QR 알고리즘의 총 비용을 분석한다.

하이젠베르크 환원: O(\frac{10}{3} n^3) (하우스홀더 변환, 대칭이면 삼대각 환원 O(\frac{4}{3} n^3)).

QR 반복: 반복당 O(n^2) (하이젠베르크 행렬의 QR 분해), 총 반복 횟수는 적절한 시프트 하에서 O(n) (각 고유값의 디플레이션에 평균 2~3회 반복). 따라서 QR 반복의 총 비용은 O(n^3)이나, 상수가 하이젠베르크 환원보다 작다.

전체 비용: O(n^3). 대칭 삼대각 행렬의 경우 QR 반복의 총 비용은 O(n^2)로 축소된다.

8. 예외적 시프트(Exceptional Shift)

극히 드물게, 시프트된 QR 반복이 수렴하지 않고 정체(stagnation)하는 경우가 있다. 이를 감지하면 **예외적 시프트(exceptional shift)**를 적용한다. 이는 무작위적이거나 비표준적인 시프트 값(예: 임의의 랜덤 시프트, 또는 Wilkinson 시프트의 보정)을 사용하여 정체 상태를 탈출하는 전략이다.

LAPACK의 구현에서는 10~20회의 반복 후에도 디플레이션이 발생하지 않으면 예외적 시프트를 적용한다.

9. 실용적 구현에서의 고려 사항

벌거(bulge) 추적. Francis의 암묵적 이중 시프트 QR 단계는 “벌거 추적(bulge chasing)“으로 구현된다. 하이젠베르크 행렬의 좌상단에 도입된 벌거(하이젠베르크 구조를 위반하는 비영 성분)를 기븐스 회전으로 순차적으로 우하단으로 이동시켜 소거하는 과정이다.

다중 시프트. 현대적 구현에서는 2개 이상의 시프트를 동시에 적용하는 다중 시프트(multi-shift) QR 알고리즘이 사용된다. 이는 BLAS-3 수준의 행렬-행렬 곱을 활용하여 캐시 효율을 극대화한다.

공격적 조기 디플레이션(aggressive early deflation). Braman, Byers, Mathias (2002)가 제안한 기법으로, 우하단 블록에서 디플레이션 가능한 고유값을 조기에 감지하여 반복 횟수를 대폭 줄인다. 이 기법은 LAPACK 최신 버전에 구현되어 있다.