固有値分解 — Hessenberg・Francis QR・対称三重対角・Krylov

概要

固有値分解は線形作用素の本質的な不変量を抽出する操作で、応用範囲は微分方程式の解析、振動系の固有モード、グラフラプラシアン、量子力学のスペクトル、PCA の主成分（実は SVD）など線形代数の中で最も広い。 正方行列 $A \in \mathbb{R}^{n \times n}$ に対し、固有値 $\lambda$ と零でない固有ベクトル $\mathbf{v}$ は

$$A \mathbf{v} = \lambda \mathbf{v}$$

を満たす。 対角化可能であれば $A = V \Lambda V^{-1}$、対称であれば $A = V \Lambda V^\top$ と直交対角化される。 正方とはいえ、計算アルゴリズムは行列の構造（対称、Hermitian、Hessenberg、帯、疎）によって大きく分岐する。

本記事では一般非対称・対称・一般化固有値の標準解法、ベキ乗反復系、Krylov 法、Gershgorin による包含定理、Bauer-Fike 摂動を扱う。 API は LinAlg API — 固有値、Gershgorin 円盤 を参照。

一般行列: Hessenberg + Francis QR

一般非対称行列 $A$ の全固有値計算は、2 段階で行う。

第 1 段: Hessenberg 縮約

Householder 反射を両側相似変換として適用し、$A$ を上 Hessenberg 行列 $H$（subdiagonal のみ非零、それより下は零）に変換する:

$$Q^\top A Q = H, \quad H_{ij} = 0 \text{ for } i > j + 1$$

$Q$ は $n - 2$ 個の Householder の積。フロップ数は $\tfrac{10}{3} n^3$、対称性があれば三重対角まで縮約され $\tfrac{4}{3} n^3$ になる。 この縮約により、後続の QR 反復で各ステップが $O(n^2)$ に落ちる（Hessenberg は QR 反復で形が保たれる）。

第 2 段: Francis double-shift QR

Hessenberg 上で陰的シフト QR 反復を回す。 各ステップでシフト $\sigma$ を選び

$$H - \sigma I = Q R, \quad H' = R Q + \sigma I = Q^\top H Q$$

として相似変換を維持しつつ Hessenberg 構造を保つ。$|H_{n,n-1}|$ が機械精度に達すると最後の固有値が確定し、行列を縮約して同じ手続きを再帰的に適用する。

非対称実行列では複素固有値が共役ペアで現れるため、複素シフト 1 個では実数演算で閉じない。 Francis double-shift は複素共役ペア $(\sigma, \bar{\sigma})$ を一度に適用する 2 ステップを 1 つの実数演算 sweep に統合する:

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

$M$ の第 1 列のみから 1 個の Householder を構成し、続いてバルジ（subdiagonal の下に出現する一時的非零要素）を Givens 連鎖で下方向に追い出す。 このimplicit Q theorem によって毎ステップで陽に QR を計算する必要がなくなり、1 sweep $O(n^2)$ で 2 個の固有値が同時に得られる。

シフト戦略

標準は Wilkinson シフト（$2 \times 2$ 末尾ブロックの 2 つの固有値のうち $h_{nn}$ に近い方）。 ill-conditioned な並んだ固有値や周期的振る舞いに陥った場合は、ad hoc シフトや random shift で脱出する。 QR 反復は典型的に $O(n)$ sweep、トータル $O(n^3)$ で全固有値が確定する。

sangi の eigen は MKL 利用可能時は LAPACK dgeev を呼び、ない場合は内部実装の Hessenberg + Francis QR で同じ結果を得る。

対称行列: 三重対角化と MR3

$A = A^\top$ のとき、Hessenberg 縮約は自動的に対称三重対角に縮約される:

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

$T$ の固有値・固有ベクトルを求めるアルゴリズムは複数あり、それぞれ得手がある:

対称 QR 反復

三重対角に対する陰的 QR 反復で、Wilkinson シフトを使うと 1 sweep $O(n)$ で進む（バルジが $2 \times 2$ なので Givens 1 個で押し出せる）。 全固有値・固有ベクトル取得で $O(n^3)$。

分割統治

$T$ を中央で 2 つの部分三重対角に分け、再帰的に固有値分解した後で secular equation で結合する。 計算量は典型的に $O(n^{2.3})$、$n \gtrsim 100$ で QR 反復より速い。LAPACK DSTEDC。

MR3（Multiple Relatively Robust Representations）

Dhillon-Parlett (2004) によるアルゴリズムで、各固有値クラスタごとに「相対的にロバストな」$LDL^\top$ 表現を選び、 その表現上で $O(n)$ で固有ベクトルを計算する。全体で$O(n^2)$ で全固有値・固有ベクトルが得られる、最も高速な実装。 LAPACK DSTEMR、Eigen の自己内蔵実装も同方針。

計算量比較:

sangi の eigen_symmetric は MKL/OpenBLAS があれば現状 LAPACK DSYEV (三重対角化 + QR 反復) を呼び、ない場合は Jacobi 法にフォールバックする。 大行列向けの DSYEVR (MR3 / DSTEMR) や DSYEVD (分割統治) への切替は今後の課題である。

Hermitian 行列

複素エルミート行列 $A^H = A$ の固有値は実数。 MKL 利用時は LAPACK ZHEEV、ない場合は $n \times n$ エルミートを $2n \times 2n$ 実対称に同型変換して eigen_symmetric で解く実装にフォールバックする。

一般化固有値と QZ

一般化固有値問題 $A \mathbf{x} = \lambda B \mathbf{x}$ は、$B$ が正則であれば $B^{-1} A \mathbf{x} = \lambda \mathbf{x}$ に帰着できるが、 $B$ が悪条件・特異な場合は数値的に破綻する。代わりに $A, B$ をペアとして同時に Schur 形に縮約する QZ アルゴリズムを使う。

$QZ$ は両側ユニタリ相似で $A, B$ を上三角に同時変換する:

$$Q^\top A Z = S, \quad Q^\top B Z = T, \quad S, T \text{ 上三角}$$

固有値ペアは $(\sigma_i, \tau_i) = (S_{ii}, T_{ii})$、固有値は $\lambda_i = \sigma_i / \tau_i$。 $\tau_i = 0$ の場合は無限固有値（$B$ の零空間方向）として扱える。

第 1 段で Hessenberg-三角形ペアに縮約し、第 2 段で QZ 反復（Francis QR の一般化）で対角化する。 LAPACK DGGEV。 $A$ 対称・$B$ 正定値の場合は Cholesky $B = L L^\top$ で $L^{-1} A L^{-\top}$ の対称固有値問題に直接変換できる（sangi eigen_generalized_symmetric）。

ベキ乗法・逆反復・Rayleigh quotient

全固有値が不要で、極端な値（絶対最大・最小・特定値近く）のみ取りたい場合のクラシックな反復系。

ベキ乗法

初期ベクトル $\mathbf{x}_0$ に対して $\mathbf{x}_{k+1} = A \mathbf{x}_k / \|A \mathbf{x}_k\|$ を反復すると、 $\mathbf{x}_k$ は最大絶対固有値 $\lambda_1$ に対応する固有ベクトルに収束し、収束率は $|\lambda_2 / \lambda_1|^k$。 計算は行列-ベクトル積のみで $O(n^2)$/反復、疎行列なら $O(\mathrm{nnz})$/反復。 スペクトルギャップが大きい行列で実用的。Google PageRank の素朴版がこれ。

逆反復

シフト $\mu$ を与え、$(A - \mu I)^{-1} \mathbf{x}_k$ を解いて反復する:

$$\mathbf{x}_{k+1} = (A - \mu I)^{-1} \mathbf{x}_k / \|(A - \mu I)^{-1} \mathbf{x}_k\|$$

$\mu$ に最も近い固有値に対応する固有ベクトルに収束し、収束率は $|(\lambda_i - \mu) / (\lambda_j - \mu)|^k$（$\lambda_i$ が最近、$\lambda_j$ が次に近い）。 近似固有値が得られていればすぐに固有ベクトルを精密化できるため、QR 反復で固有値を出した後の固有ベクトル計算に使う。

Rayleigh quotient iteration

対称行列向け。シフト $\mu_k$ を毎反復で現在の近似 $\mathbf{v}_k$ の Rayleigh quotient

$$\mu_k = \frac{\mathbf{v}_k^\top A \mathbf{v}_k}{\mathbf{v}_k^\top \mathbf{v}_k}$$

に更新する。$\mathbf{v}_k$ が真の固有ベクトルから $O(\varepsilon)$ ずれていると Rayleigh quotient は $O(\varepsilon^2)$ ずれた固有値で、 逆反復がそのシフトで動くため次の $\mathbf{v}_{k+1}$ は $O(\varepsilon^3)$ になる。三次収束で対称固有値の極めて高速な精密化に使われる。 非対称では二次収束に留まる。

Krylov 法: Lanczos と Arnoldi

$n$ が極めて大きく $A$ が疎な場合、Hessenberg 縮約自体が $O(n^3)$ で実行不可能。 Krylov 部分空間 $\mathcal{K}_k(A, \mathbf{v}) = \mathrm{span}(\mathbf{v}, A\mathbf{v}, \ldots, A^{k-1}\mathbf{v})$ への射影で $k \ll n$ 個の極値固有値だけを取り出す方法が必要になる。

Lanczos 法（対称）

対称 $A$ に対し、Krylov 基底を逐次正規直交化すると、過去 2 ベクトルとの直交化のみで次のベクトルが得られる三項漸化式:

$$A \mathbf{q}_j = \beta_{j-1} \mathbf{q}_{j-1} + \alpha_j \mathbf{q}_j + \beta_j \mathbf{q}_{j+1}$$

得られる射影 $T_k = Q_k^\top A Q_k$ は対称三重対角で、$T_k$ の固有値 (Ritz value) が $A$ の極値固有値の近似を与える。 1 反復 $O(n)$（疎行列なら $O(\mathrm{nnz})$）、$k$ 反復で $k$ 個の固有値が得られる。 理論上は厳密だが、丸め誤差で $\mathbf{q}_j$ の直交性が失われる再直交化問題があり、選択的再直交化や Paige の解析に基づく対応が必要。

Arnoldi 法（非対称）

非対称行列では三項漸化に縮約されず、Arnoldi 反復は MGS で過去全ベクトルとの直交化を行い、上 Hessenberg $H_k$ を生成する。 $k$ 反復のコストは $O(nk)$ となる（疎なら $O(\mathrm{nnz} \cdot k)$）。 $H_k$ の固有値が $A$ の Ritz value、極値固有値の近似。GMRES の核でもある。

Shift-and-Invert

絶対最小の固有値や、内部スペクトルの特定値近くだけ欲しい場合は、$(A - \mu I)^{-1}$ に対して Lanczos / Arnoldi を回す shift-and-invert 戦略を使う。固有値マッピング $\lambda \mapsto 1/(\lambda - \mu)$ により $\mu$ 近傍の固有値が極大になり、Krylov 法が速く収束する。 毎反復で $(A - \mu I)$ を係数とする 1 本の連立方程式を解く必要があり、LU 前段がしばしば最大のコスト。

再開戦略

Krylov 基底が大きくなると直交化コストと記憶量が増える。暗黙再開 Arnoldi (IRA)（Sorensen, ARPACK）は 収束済みの Ritz vector を保持しつつ基底を $m$ 次元に絞り込み、新しい $k - m$ ベクトルで拡張する反復で、 実用的な大規模固有値計算の標準。

Gershgorin 円盤定理

Gershgorin の包含定理は、行列の固有値の位置を一切の反復計算なしに矩形領域に束縛する基本道具である:

$$\text{全固有値} \in \bigcup_{i=1}^{n} D_i, \quad D_i = \left\{z \in \mathbb{C} : |z - a_{ii}| \leq \sum_{j \neq i} |a_{ij}|\right\}$$

各 Gershgorin 円盤 $D_i$ は対角要素 $a_{ii}$ を中心、行外の絶対和を半径とする閉円盤。 さらに「$k$ 個の円盤が連結成分を成し、残りと交わらないなら、その連結成分には正確に $k$ 個の固有値が含まれる」という強形が成立する。

応用

• 正定値性の判定: 全ての $D_i$ が複素右半平面に含まれれば全固有値の実部 $> 0$、対称なら正定値。
• 収束の事前判定: Jacobi 法（線形ソルバー）のスペクトル半径 $\rho(B_J) < 1$ の十分条件として、対角優位（$|a_{ii}| > \sum_{j \neq i} |a_{ij}|$）が Gershgorin から従う。
• シフト選択: シフト QR 反復のシフト候補や、shift-and-invert の $\mu$ 選択で固有値範囲の目安として使う。

列に対する Gershgorin も成立し（$A$ と $A^\top$ の固有値が同じため）、行 Gershgorin と列 Gershgorin の交差がより厳しい囲い込みを与える。 sangi は Gershgorin 円盤関数で行/列両方の円盤を返す。

条件数と摂動（Bauer-Fike）

固有値計算の数値安定性は、入力の摂動が固有値にどう伝播するかで決まる。 Bauer-Fike 定理は対角化可能な行列に対する固有値摂動の上界を与える:

$$A = V \Lambda V^{-1}, \quad \min_i |\hat{\lambda} - \lambda_i| \leq \kappa_2(V) \cdot \|E\|_2$$

ここで $\hat{\lambda}$ は摂動行列 $A + E$ の固有値、$\kappa_2(V) = \|V\|_2 \|V^{-1}\|_2$ は固有ベクトル行列の条件数。

結果として:

• 対称行列: $V$ が直交、$\kappa_2(V) = 1$ なので $|\hat{\lambda} - \lambda_i| \leq \|E\|_2$（最良条件、Weyl の不等式）。
• 正規行列: $V$ がユニタリで同上。
• 非対称・非正規: $\kappa_2(V)$ が大きくなる可能性、固有値が摂動に弱い。Wilkinson 行列の極端な例では $\kappa_2(V) = 10^{20}$ オーダー。

個別の固有値ごとには Wilkinson の固有値条件数もあり、左右固有ベクトル $\mathbf{u}_i, \mathbf{v}_i$ に対し $s_i = 1 / |\mathbf{u}_i^\top \mathbf{v}_i|$ で定義され、$s_i \cdot \|E\|_2$ が個別摂動の上界となる。 defective（固有ベクトルが線形独立でない）に近い行列では $s_i$ が発散する。

sangi の選択

参考文献

• Francis, J. G. F. (1961, 1962). "The QR Transformation, I and II". Computer J., 4, 265–271, 332–345.
• Parlett, B. N. (1998). The Symmetric Eigenvalue Problem, SIAM (classic ed.).
• Dhillon, I. S. and Parlett, B. N. (2004). "Orthogonal eigenvectors and relative gaps". SIAM J. Matrix Anal. Appl., 25(3), 858–899.
• Lanczos, C. (1950). "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators". J. Res. Nat. Bur. Standards, 45, 255–282.
• Arnoldi, W. E. (1951). "The principle of minimized iterations in the solution of the matrix eigenvalue problem". Quart. Appl. Math., 9, 17–29.
• Bauer, F. L. and Fike, C. T. (1960). "Norms and exclusion theorems". Numer. Math., 2, 137–141.
• Stewart, G. W. (2001). Matrix Algorithms, Volume II: Eigensystems, SIAM.