고유값의 수치 계산

특성 다항식으로 구하면 안 되는 이유

이것이 2장(조건수와 안정성)에서 본 윌킨슨 다항식입니다. 계수를 배정밀도로 저장만 해도(전개 계산 자체는 정확히 해도) 반올림 오차가 생기고, 그 근을 정확히 구해도 $\lambda = 12, 13, 14$ 근처에서 오차가 $6\times10^{-4}$ 정도가 되고, 근을 배정밀도로 계산하면 오차가 더 커집니다. 윌킨슨이 1959년에 보인 대로 $\lambda^{19}$의 계수 $-210$을 $2^{-23} \approx 1.2\times10^{-7}$만큼만 바꾸면 근 열 개가 복소수 쌍이 되고, $16.73 \pm 2.81i$ 같은 근이 나타납니다. 다항식의 근은 계수에 대해 조건이 나쁜 문제입니다.

그런데 행렬의 고유값 자체는 그렇지 않을 수 있습니다. 문제를 "다항식의 근"으로 바꾸는 순간 원래 문제보다 훨씬 예민한 문제로 바꾼 것입니다. 행렬 쪽의 민감도는 다음 정리가 말해 줍니다.

이제 결론이 분명합니다. 고유값을 구하는 좋은 알고리즘은 다항식을 거치지 않고, 행렬에 직교 닮음 변환 $A \mapsto Q^{\mathsf T}AQ$만 쌓아서 고유값을 읽을 수 있는 꼴(삼각·대각)로 만듭니다. 직교 변환은 노름을 보존하므로 반올림 오차가 커지지 않고, 결과는 "$\|E\| \approx \varepsilon_{\text{mach}}\|A\|$만큼 교란된 행렬의 정확한 고유값"이 됩니다(후방 안정). 그다음 정확도는 바우어–파이크가 말하는 문제 자체의 조건에만 달려 있습니다.

거듭제곱법과 수렴 속도

$A$가 고유벡터 $\mathbf v_1, \dots, \mathbf v_n$으로 대각화되고 $|\lambda_1| \gt |\lambda_2| \ge \cdots \ge |\lambda_n|$이라 합시다. 시작 벡터 $\mathbf x_0 = \sum c_i\mathbf v_i$에 $A$를 $k$번 곱하면

function power_iteration(A, x, tol, maxit)
  x = x / norm(x)
  for k = 1 .. maxit
    y = A x                                // 반복당 행렬–벡터 곱 한 번
    lam = dot(x, y)                        // 레일리 몫 xᵀAx (‖x‖ = 1)
    if norm(y - lam * x) <= tol * abs(lam) then return (lam, x)
    x = y / norm(y)
  return (lam, x)

역반복, 이동, 레일리 몫 반복

$A\mathbf v = \lambda\mathbf v$이면 $(A - \sigma I)^{-1}\mathbf v = \frac{1}{\lambda - \sigma}\mathbf v$입니다. 고유벡터는 그대로이고 고유값만 $1/(\lambda_i - \sigma)$로 바뀝니다. 이동량 $\sigma$에 가장 가까운 고유값이 절댓값이 가장 큰 고유값이 되므로, 이 행렬에 거듭제곱법을 쓰면 그 고유벡터로 갑니다. 역행렬을 만들지 않고 매번 $(A - \sigma I)\mathbf y = \mathbf x$를 풀며, LU 분해는 처음 한 번만 하면 됩니다.

한 걸음 더 나가면 매 단계 $\sigma$를 지금의 가장 좋은 추정, 곧 레일리 몫으로 바꿀 수 있습니다. 이것이 레일리 몫 반복(Rayleigh quotient iteration, RQI)입니다.

function rayleigh_quotient_iteration(A, x, tol)
  x = x / norm(x);  rho = dot(x, A x)
  repeat
    y = solve(A - rho I, y = x)            // 매번 새 행렬 분해: 비싸다
    x = y / norm(y)
    rho = dot(x, A x)                      // 이동량을 갱신
  until norm(A x - rho x) <= tol * norm(A)
  return (rho, x)

게르시고린 원판: 계산 없이 고유값 가두기

13장의 엄격한 대각 우세는 모든 원판이 원점을 포함하지 않는다는 뜻이고, 그래서 대각 우세 행렬은 가역입니다. 대각이 양수인 대칭 대각 우세 행렬은 모든 고유값이 양수, 곧 양정치입니다. 원판이 너무 크면 대각 닮음 변환 $DAD^{-1}$(고유값은 같다)으로 $i$행의 비대각 원소를 $d_i$배, $i$열을 $1/d_i$배 해서 원하는 원판을 줄일 수 있습니다.

수축: 찾은 고유값을 치우기

거듭제곱법으로 $\lambda_1, \mathbf v_1$을 찾았다면 다음 고유값은 어떻게 찾을까요? 찾은 것을 행렬에서 "치워" 두 번째 것이 가장 크게 만드는 과정을 수축(deflation)이라 합니다.

예를 들어 $A = \begin{bmatrix}2 & 1\\ 1 & 2\end{bmatrix}$에서 $\lambda_1 = 3$, $\mathbf v_1 = (1, 1)/\sqrt2$이면 $A_1 = A - 3\mathbf v_1\mathbf v_1^{\mathsf T} = \begin{bmatrix}0.5 & -0.5\\ -0.5 & 0.5\end{bmatrix}$이고, 고유값은 0과 1입니다. 남은 고유값 1을 거듭제곱법으로 바로 얻습니다.

QR 알고리즘: 모든 고유값을 한꺼번에

먼저 기본 성질입니다. $R_k = Q_k^{\mathsf T}A_k$이므로 $A_{k+1} = Q_k^{\mathsf T}A_kQ_k$, 직교 닮음이라 고유값이 보존됩니다. 그리고 $\bar Q_k = Q_0Q_1\cdots Q_{k-1}$의 열들은 $A^k$의 열에 그람–슈미트를 한 결과와 같습니다. 곧 QR 알고리즘(QR algorithm)은 $n$개 벡터에 거듭제곱법을 동시에 하며 서로 직교하게 유지하는 동시 반복입니다. 실용 알고리즘은 여기에 세 가지를 더합니다.

function shifted_qr_eigenvalues(A)
  H = hessenberg(A)                        // 하우스홀더, 한 번만: 10n³/3
  m = n
  while m > 1
    mu = wilkinson_shift(H[m-1..m, m-1..m])  // 오른쪽 아래 2×2의 고유값 중 h_mm에 가까운 것
    Q R = qr(H[1..m,1..m] - mu I)            // 기븐스 m-1개: O(m²)
    H[1..m,1..m] = R Q + mu I
    if |H[m][m-1]| <= eps * (|H[m][m]| + |H[m-1][m-1]|)
      H[m][m-1] = 0;  m = m - 1                // 수축: 고유값 H[m][m] 확정
  return diag(H)

이동이 두 번째 열쇠입니다. $A_k - \mu_kI = Q_kR_k$, $A_{k+1} = R_kQ_k + \mu_kI$는 여전히 닮음 변환이고, 마지막 대각 아래 원소는 대략 $|\lambda_n - \mu|/|\lambda_{n-1} - \mu|$의 비율로 줄어듭니다. 역반복과 똑같은 구조입니다. 실제로 이동 QR의 마지막 행은 $\mathbf e_n$에 이동 역반복을 하는 것과 같습니다. $\mu = h_{nn}$(레일리 이동)을 쓰면 RQI처럼 2차(대칭이면 3차) 수렴합니다. 레일리 이동은 대칭적인 경우에 막힐 수 있어서, 오른쪽 아래 2 × 2 블록의 고유값 중 $h_{nn}$에 가까운 것을 쓰는 윌킨슨 이동(Wilkinson shift)이 표준이며, 대칭 삼중대각 행렬에서는 항상 수렴함이 증명되어 있습니다. 세 번째 열쇠는 수축입니다. $h_{m,m-1}$이 반올림 수준이 되면 잘라 내고 $(m-1)$차로 계속합니다.

표의 마지막 줄이 윌킨슨 이동이 필요한 이유입니다. 반대각 행렬은 $h_{nn} = 0$이고 고유값 $\pm1$이 0에서 같은 거리에 있어, 이동이 없거나 레일리 이동이면 QR 단계가 행렬을 전혀 바꾸지 못합니다. 고유값이 고유값 하나당 평균 두세 단계 만에 확정되므로 전체 비용은 축소 비용이 지배하고, 비대칭 행렬의 고유값 전부가 약 $10n^3$ flop, 대칭이면 약 $\tfrac43n^3$ flop에 나옵니다. 실수 비대칭 행렬의 복소 고유값 쌍은 켤레 이동 두 개를 실수 연산으로 한 번에 적용하는 프랜시스의 이중 이동으로 처리하며, 결과는 대각에 2 × 2 블록이 남는 실수 슈르 꼴입니다.

대칭 행렬과 야코비 회전

대칭 행렬에는 훨씬 오래된, 그리고 지금도 쓰이는 방법이 있습니다. 비대각 원소 하나 $a_{pq}$를 골라, $p, q$ 좌표 평면의 회전 $J$로 $J^{\mathsf T}AJ$의 $(p, q)$ 원소를 0으로 만듭니다. 2 × 2 대칭 행렬을 대각화하는 회전각은 언제나 있습니다.

야코비는 QR보다 몇 배 느리지만, 작은 고유값까지 상대 정확도가 높고(양정치 행렬에서 특히), 회전 하나하나가 독립적이라 병렬화하기 쉽습니다. 고유벡터는 회전들의 곱 $V = J_1J_2\cdots$으로 저절로 나옵니다. 엔진의 MB.mat.eigSym도 이 방법입니다.

특이값 계산과 저계수 근사

$m \times n$ 행렬 $A$의 특이값 $\sigma_i$는 $A^{\mathsf T}A$의 고유값의 제곱근입니다(LinearBook: SVD). 그렇다고 $A^{\mathsf T}A$를 만들어 고유값을 구하면 안 됩니다. 조건수가 제곱되기 때문입니다.

표준 알고리즘(골럽–카한)은 $A$의 왼쪽과 오른쪽에 서로 다른 하우스홀더 반사를 번갈아 곱해 이중대각 행렬 $B = U_1^{\mathsf T}AV_1$을 만들고($O(mn^2)$), $B$에 $B^{\mathsf T}B$를 만들지 않은 채 암묵적 이동 QR을 적용합니다. $A^{\mathsf T}A$에 대한 대칭 QR 알고리즘을 $A$ 자체에서 흉내 내는 것입니다. 한쪽 야코비 회전으로 $A$의 열을 서로 직교하게 만드는 방법도 쓰입니다.

예를 들어 5 × 5 힐베르트 행렬 $h_{ij} = 1/(i + j - 1)$의 특이값은 1.567, 0.2085, 0.01141, 3.06×10⁻⁴, 3.29×10⁻⁶입니다. 계수 2의 근사도 2-노름 오차 0.0114, 곧 $\|A\|_2$의 0.7%로 맞습니다. 특이값이 빨리 줄어드는 행렬은 적은 수의 벡터로 압축되며, 이것이 영상 압축, 추천 시스템, 주성분 분석의 수학입니다. 아주 큰 행렬에서는 무작위 벡터 $k + 10$개 정도에 $A$를 곱해 치역을 근사한 뒤 작은 SVD를 하는 무작위 SVD가 쓰입니다.

응용: 페이지랭크, 진동 모드, 에너지 준위

페이지 $j$에서 나가는 링크가 $d_j$개면 $P_{ij} = 1/d_j$($j \to i$ 링크가 있을 때). 나가는 링크가 없는 페이지는 모든 페이지로 균등하게 간다고 둡니다. 서퍼는 확률 $\alpha$(감쇠 인자, 보통 0.85)로 링크를 따르고 $1 - \alpha$로 아무 페이지로나 순간 이동합니다.

$G$는 꽉 찬 행렬이지만 만들 필요가 없습니다. $G\mathbf r = \alpha P\mathbf r + \frac{1 - \alpha}{n}\mathbf 1$($\sum r_i = 1$일 때)이므로 희소 행렬 $P$와의 곱 한 번이면 됩니다. $\alpha = 0.85$이면 $0.85^k = 10^{-8}$에서 $k \approx 113$번, 웹 전체 규모에서도 같은 횟수입니다.

진동 모드. 질량 $m$인 물체 $n$개를 강성 $k$인 용수철로 일렬로 잇고 양 끝을 고정하면 운동 방정식은 $M\ddot{\mathbf x} = -K\mathbf x$, $K = k\,\mathrm{tridiag}(-1, 2, -1)$이고, $\mathbf x = \mathbf v\cos\omega t$를 넣으면 $K\mathbf v = \omega^2M\mathbf v$, 곧 일반화 고유값 문제입니다. 이 사슬에서는 $\omega_j^2 = \frac{k}{m}\bigl(2 - 2\cos\frac{j\pi}{n+1}\bigr)$를 정확히 알지만, 실제 구조물의 유한 요소 모델은 크기가 수십만이고 필요한 것은 가장 낮은 진동수 몇십 개입니다. 그래서 $\sigma = 0$ 근처의 이동 역반복을 크릴로프 공간과 결합한 이동–역 란초시 방법(ARPACK, MATLAB eigs)이 쓰입니다(LinearBook: 고유값).

주성분 분석. 데이터 행렬의 공분산 $C = X^{\mathsf T}X/(N - 1)$의 가장 큰 고유벡터들이 데이터가 가장 많이 퍼진 방향입니다(LinearBook: PCA). 수치적으로는 예제 14.6의 교훈대로 $C$를 만들지 않고 중심화한 $X$의 SVD를 구합니다.

양자 에너지 준위. 1차원 조화 진동자 $-\tfrac12\psi'' + \tfrac12x^2\psi = E\psi$(단위를 골라 정확한 준위는 $E_j = j + \tfrac12$)를 2배 하여 $-\psi'' + x^2\psi = 2E\psi$로 쓰고, $[-6, 6]$을 $n$개 내부점으로 나눠 2차 중앙 차분을 쓰면 대칭 삼중대각 행렬의 고유값 문제가 됩니다. 2E의 정확한 값 1, 3, 5, 7, …과 비교하면

바닥 상태의 오차는 0.00347, 0.00088, 0.00022, 0.00006으로 $h$가 반이 될 때마다 1/4, 곧 차분의 $O(h^2)$를 그대로 따릅니다. 높은 준위일수록 파동 함수가 빨리 진동해 같은 격자에서 오차가 큽니다. 고유값 계산은 스투름 수열을 쓴 이분법(삼중대각 대칭 행렬에서 $\lambda$보다 작은 고유값의 수를 $O(n)$에 센다)으로 했습니다.

이 도구가 쓰이는 곳