# 対称固有値問題と特異値分解の計算：三重対角化・Sturm 列・Jacobi 法・精度の保証 > 対称行列の三重対角化、Sturm 列による二分法、Jacobi 法の収束を証明し、Weyl の不等式による固有値の誤差、残差による固有ベクトルの sin θ 評価、Jordan–Wielandt 行列による特異値の摂動、二重対角化と Golub–Kahan の三重対角行列を導いて、eigh・svd の精度を数値で確かめる。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/symmetric-eigenproblem-and-svd ## 0. この記事の要点 対称な行列の固有値（行列の「伸び縮みの大きさ」）や、一般の行列の特異値（行列がどれだけ空間をゆがめるかを表す数）を、コンピュータはどう計算し、その答えはどこまで信用できるのか。この章はその 2 つの問いに答えます。まず行列を扱いやすい形（帯状の形）に直し、その形に合った計算法で固有値を求めます。そのうえで、丸め誤差があっても答えの誤差がどのくらいに収まるかを、証明と数値実験で確かめます。読み終えると、`eigh` や `svd` の結果がなぜ信頼できるのか、固有ベクトルだけは信頼できない場合があるのはなぜか、を説明できるようになります。 - 実対称行列は [Householder 鏡映](/mathematics/numerical-linear-algebra/qr-factorization#def-householder)を使った変換で対称な三重対角行列（対角とそのすぐ隣にしか成分がない行列）に直せます（演算量 $4n^3/3$）。三重対角行列の固有値は、[QR 法](/mathematics/numerical-linear-algebra/qr-algorithm#def-qr-iteration)（1 回 $O(n)$）、Sturm 列（固有値が $x$ より小さいものの個数を数える漸化式）による二分法（$k$ 番目の固有値を指定して求められる）、分割統治法・MRRR（どちらも §7 で紹介する、全固有ベクトルを速く求める方法）などで求めます。 - **Jacobi 法**は、2 行 2 列の回転で非対角成分を 1 つずつ消していく算法です。1 回の回転で非対角成分の 2 乗和はちょうど $2a_{pq}^2$ 減り、最大の成分を選ぶ古典的 Jacobi 法は $\mathrm{off}(A_{k+1})^2 \le (1 - 2/(n(n-1)))\,\mathrm{off}(A_k)^2$ で収束します。 - 対称行列の固有値計算が後退安定（出てきた答えが、元の行列を少しだけ動かした行列の厳密な答えになっていること。つまり計算した固有値が $\|E\|_2 \le cnu\|A\|_2$ の $A + E$ の厳密な固有値）なら、Weyl の不等式により各固有値の絶対誤差は $cnu\|A\|_2$ 以下です。固有ベクトルの誤差は固有値の隔たり $\delta$ で割った $\|\boldsymbol{r}\|/\delta$ で抑えられ（sin θ 評価）、近い固有値の固有ベクトルは不安定です。 - 特異値は Jordan–Wielandt 行列 $\begin{pmatrix} O & A \\ A^{\mathsf{T}} & O \end{pmatrix}$ の固有値 $\pm\sigma_i$ として対称固有値問題に帰着し、$|\sigma_i(A + E) - \sigma_i(A)| \le \|E\|_2$ が成り立ちます。$A^{\mathsf{T}}A$ の固有値から計算すると、小さい特異値は桁を失います。 - SVD の計算は、Householder 変換による二重対角化 $A = U\begin{pmatrix} B \\ O \end{pmatrix}V^{\mathsf{T}}$ と、二重対角行列 $B$ を対角が $0$ の対称三重対角行列（Golub–Kahan の行列）の固有値問題に読み替えることで行います。 ## 1. 動機 研究で最も多く呼ばれる固有値の関数は、対称（[Hermite 行列](/mathematics/linear-algebra/spectral-theorem#def-hermitian-unitary)）行列の `eigh`・`eigvalsh` と、特異値分解の `svd` です。対称行列の固有値問題は、一般の行列より 3 つの点で恵まれています。固有値は実数で、固有ベクトルは直交にとれ（）、固有値の摂動は行列の摂動のノルムで抑えられる（条件数 1）ことです。 本章では、[QR 法の章](/mathematics/numerical-linear-algebra/qr-algorithm) の算法を対称行列に特化し、さらに対称行列だけで使える算法（Sturm 列の二分法、Jacobi 法）を加えます。そのうえで「計算した固有値・固有ベクトルはどこまで正しいか」を、[行列のノルムと条件数の章](/mathematics/numerical-linear-algebra/norms-and-conditioning) の Weyl の不等式の系と、残差による固有ベクトルの評価で答えます。最後に特異値分解の計算を対称固有値問題に帰着させます。 大規模な疎行列の一部の固有値を求める Lanczos 法は、[Krylov 部分空間法 I の章](/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient) で扱います。本章は密行列の全固有値を対象にします。 この章では行列は実とし、$u$ は単位丸め、ノルムは 2 ノルム $\|\cdot\|_2$ と [Frobenius ノルム](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-frobenius) $\|\cdot\|_F$ です。対称行列の固有値は、断らない限り大きい順に $\lambda_1(\cdot) \ge \cdots \ge \lambda_n(\cdot)$ と並べます。 ## 2. 三重対角化 実対称行列 $A \in \mathbb{R}^{n\times n}$（$n \ge 3$）について、$n - 2$ 個の Householder 鏡映（または単位行列）の積である直交行列 $Q$ があって、$T = Q^{\mathsf{T}}AQ$ は対称な三重対角行列になる。対称性を使うと演算量は約 $4n^3/3$ である。 この定理のねらいは、固有値を求める前に、成分の多くが $0$ の扱いやすい形へ固有値を変えずに直しておくことです。証明の前半は一般の場合の Hessenberg 化の結果を使い、後半は、鏡映を両側から掛けた結果が「階数 1 の行列 2 つの和」の形に整理できることを示して手間を数えます。 存在は、Hessenberg 化（）を対称行列に使えば、対称な Hessenberg 行列として三重対角行列が得られることから従います。 演算量: 第 $k$ 段で右下の $m\times m$ の対称ブロック $A'$（$m = n - k$）に鏡映 $P = I - \beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}$ を両側から掛けます。$\boldsymbol{p} = \beta A'\boldsymbol{v}$、$\boldsymbol{w} = \boldsymbol{p} - \tfrac{\beta}{2}(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{p})\boldsymbol{v}$ とおくと、展開して $$ PA'P = A' - \beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}A' - \beta A'\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}} + \beta^2(\boldsymbol{v}^{\mathsf{T}}A'\boldsymbol{v})\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}} = A' - \boldsymbol{v}\boldsymbol{w}^{\mathsf{T}} - \boldsymbol{w}\boldsymbol{v}^{\mathsf{T}} $$ です（$\boldsymbol{v}\boldsymbol{p}^{\mathsf{T}} + \boldsymbol{p}\boldsymbol{v}^{\mathsf{T}} - \beta(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{p})\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}} = \boldsymbol{v}\boldsymbol{w}^{\mathsf{T}} + \boldsymbol{w}\boldsymbol{v}^{\mathsf{T}}$）。対称な行列とベクトルの積 $A'\boldsymbol{v}$ が約 $2m^2$、対称な階数 2 の更新が下三角だけで約 $2m^2$ なので、1 段は約 $4m^2$、合計 $\sum_m 4m^2 \approx 4n^3/3$ です。 一般の Hessenberg 化の $10n^3/3$ に比べ、対称性で手間が半分以下になります。固有ベクトルも要るときは $Q$ を累積する手間（約 $4n^3/3$）が加わります。 ## 3. Sturm 列による二分法 対称な三重対角行列 $T$ の対角を $\alpha_1, \ldots, \alpha_n$、副対角を $\beta_1, \ldots, \beta_{n-1}$ とし、$N(x)$ を「$x$ より小さい $T$ の固有値の個数（重複度込み）」とします。計算機援用証明の章の Sturm 列の定理（）により、漸化式 $$ d_1(x) = \alpha_1 - x,\qquad d_k(x) = \alpha_k - x - \frac{\beta_{k-1}^2}{d_{k-1}(x)} $$ の値がすべて定義されて $0$ でない $x$ では、$N(x) = \#\{k : d_k(x) < 0\}$ が $O(n)$ の演算で求まります。 次の命題の考え方は、数当てゲームの「半分に絞る」方法と同じです。「$x$ より小さい固有値は何個か」が安く数えられるので、$x$ を動かして個数が $k - 1$ 個から $k$ 個に変わる場所を挟み、区間を半分ずつ狭めて $k$ 番目の固有値に近づけます。 $T$ の固有値を小さい順に $\mu_1 \le \cdots \le \mu_n$ とし、$1 \le k \le n$ を固定する。$a_0 < b_0$ が $N(a_0) \le k - 1$、$N(b_0) \ge k$ を満たすとする（たとえば Gershgorin の円板 から $a_0 = \min_i(\alpha_i - |\beta_{i-1}| - |\beta_i|)$ とし、$b_0$ をその最大より大きくとる。$\beta_0 = \beta_n = 0$）。$j = 0, 1, \ldots$ について $c_j = (a_j + b_j)/2$ とし、$N(c_j) \ge k$ なら $(a_{j+1}, b_{j+1}) = (a_j, c_j)$、そうでなければ $(a_{j+1}, b_{j+1}) = (c_j, b_j)$ とする（各 $c_j$ で Sturm 列は定義されているとする）。このとき 1. すべての $j$ で $a_j \le \mu_k < b_j$、$b_j - a_j = (b_0 - a_0)/2^j$ であり、$|c_j - \mu_k| \le (b_0 - a_0)/2^{j+1}$。 2. 幅を $\tau > 0$ 以下にするのに要する反復は $\lceil\log_2((b_0 - a_0)/\tau)\rceil$ 回で、全体の演算量は $O(n\log_2((b_0 - a_0)/\tau))$ である。 $N(x) \le k - 1$ は「$x$ より小さい固有値は $k - 1$ 個以下」なので $\mu_k \ge x$ と同値で、$N(x) \ge k$ は $\mu_k < x$ と同値です。初期条件から $a_0 \le \mu_k < b_0$ で、更新の規則は $N(a_j) \le k - 1$ と $N(b_j) \ge k$ を保つので、帰納法で $a_j \le \mu_k < b_j$ です。区間は毎回半分になるので幅は $(b_0 - a_0)/2^j$ で、中点 $c_j$ と区間内の点 $\mu_k$ の距離は幅の半分以下です。2 は $(b_0 - a_0)/2^j \le \tau$ を解いたもので、1 回の $N(c_j)$ の計算は Sturm 列の定理（）により漸化式 $d_1, \ldots, d_n$ の $O(n)$ の演算です。Gershgorin の円板の定理からすべての固有値は $[a_0, \max_i(\alpha_i + |\beta_{i-1}| + |\beta_i|)]$ に入るので、例の $a_0$, $b_0$ は初期条件を満たします（$a_0$ の左に固有値はなく、$b_0$ の左にすべての固有値がある）。 二分法の利点は、$k$ を指定して「下から $k$ 番目の固有値」だけを求められること、区間の中の固有値だけを求められること、そして互いに独立なので並列に計算できることです。固有値が求まった後の固有ベクトルは、[冪乗法・逆反復の章](/mathematics/numerical-linear-algebra/power-and-inverse-iteration) の逆反復で 1〜2 回の反復で求まります。SciPy では `scipy.linalg.eigh_tridiagonal(..., lapack_driver="stebz")` がこの方法（LAPACK の stebz）です。 ## 4. Jacobi 法 Jacobi 法の発想は単純で、行列の「対角以外の成分」が全部 $0$ なら固有値は対角成分そのものなので、対角以外を回転で少しずつ $0$ に近づけていく、というものです。次の定義と補題は、回転の形と、1 回の回転で対角以外の成分の 2 乗和がどれだけ減るかを調べます。 正方行列 $A$ について $\mathrm{off}(A) = \bigl(\sum_{i\ne j}a_{ij}^2\bigr)^{1/2}$ とする。$p < q$ と $c = \cos\phi$, $s = \sin\phi$ について、単位行列の $(p,p), (p,q), (q,p), (q,q)$ 成分を $c, s, -s, c$ に置き換えた直交行列 $J = J(p, q, \phi)$ を **Jacobi 回転**という。 $A$ を実対称行列、$p < q$、$a_{pq} \ne 0$ とする。$\phi$ を $$ (\cos2\phi,\ \sin2\phi) = \frac{(a_{pp} - a_{qq},\ -2a_{pq})}{\sqrt{(a_{pp} - a_{qq})^2 + 4a_{pq}^2}} $$ と選ぶと、$B = J^{\mathsf{T}}AJ$（$J = J(p, q, \phi)$）は対称で $b_{pq} = 0$ であり、 $$ \mathrm{off}(B)^2 = \mathrm{off}(A)^2 - 2a_{pq}^2 $$ が成り立つ。 $AJ$ の第 $p$ 列は $c\,\boldsymbol{a}_p - s\,\boldsymbol{a}_q$、第 $q$ 列は $s\,\boldsymbol{a}_p + c\,\boldsymbol{a}_q$（$\boldsymbol{a}_j$ は $A$ の第 $j$ 列）で、$J^{\mathsf{T}}$ の第 $p$ 行は第 $p$ 成分 $c$・第 $q$ 成分 $-s$ だけをもちます。よって $$ b_{pq} = c(sa_{pp} + ca_{pq}) - s(sa_{qp} + ca_{qq}) = (c^2 - s^2)a_{pq} + cs(a_{pp} - a_{qq}) = a_{pq}\cos2\phi + \tfrac12(a_{pp} - a_{qq})\sin2\phi $$ で、$\phi$ の選び方を代入すると $a_{pq}(a_{pp} - a_{qq}) - a_{pq}(a_{pp} - a_{qq})$ に正の定数を掛けたもので $0$ です。 直交行列を掛けても Frobenius ノルムは変わらないので $\|B\|_F = \|A\|_F$ です。$J$ は第 $p$, $q$ 行・列しか混ぜないので、$B$ の対角成分のうち変わるのは $b_{pp}$, $b_{qq}$ だけです。さらに、$\{p, q\}$ の 2 行 2 列の部分行列について $\begin{pmatrix} b_{pp} & b_{pq} \\ b_{qp} & b_{qq} \end{pmatrix} = R^{\mathsf{T}}\begin{pmatrix} a_{pp} & a_{pq} \\ a_{qp} & a_{qq} \end{pmatrix}R$（$R = \begin{pmatrix} c & s \\ -s & c \end{pmatrix}$ は直交）なので、この部分の Frobenius ノルムも等しく、$b_{pq} = b_{qp} = 0$ から $b_{pp}^2 + b_{qq}^2 = a_{pp}^2 + a_{qq}^2 + 2a_{pq}^2$ です。対角成分の 2 乗和は $2a_{pq}^2$ 増え、全体の 2 乗和は変わらないので、非対角の 2 乗和は $2a_{pq}^2$ 減ります。 $A_0 = A$ を $n$ 次実対称行列（$n \ge 2$）とし、$A_k$ の非対角成分のうち絶対値最大の $(p, q)$（$p < q$）を選んで の回転で $A_{k+1} = J_k^{\mathsf{T}}A_kJ_k$ とする（$\mathrm{off}(A_k) = 0$ なら止める）。このとき $$ \mathrm{off}(A_{k+1})^2 \le \Bigl(1 - \frac{2}{n(n-1)}\Bigr)\mathrm{off}(A_k)^2 $$ であり、$A_k$ の対角成分を大きい順に並べた $d_1^{(k)} \ge \cdots \ge d_n^{(k)}$ は $|d_i^{(k)} - \lambda_i(A)| \le \mathrm{off}(A_k)$ を満たす。とくに対角成分は固有値に収束する。 $\mathrm{off}(A_k)^2$ は $n(n-1)$ 個の非対角成分の 2 乗和で、選んだ $a_{pq}^2$ はその最大なので $a_{pq}^2 \ge \mathrm{off}(A_k)^2/(n(n-1))$ です。 から $\mathrm{off}(A_{k+1})^2 = \mathrm{off}(A_k)^2 - 2a_{pq}^2$ で、評価を得ます。 $A_k$ は $A$ と直交相似なので固有値が等しく（積の定理 から特性多項式が等しい）、$D_k$ を $A_k$ の対角部分とすると $A_k = D_k + (A_k - D_k)$ です。対称行列 $A_k - D_k$ の 2 ノルムは Frobenius ノルム $\mathrm{off}(A_k)$ 以下なので、エルミート行列の固有値の摂動（）から、$D_k$ の固有値（＝対角成分を大きい順に並べたもの）と $A_k$ の固有値の差は $\mathrm{off}(A_k)$ 以下です。 1 掃引（$n(n-1)/2$ 回の回転）あたりの縮小率は評価からは $(1 - 2/(n(n-1)))^{n(n-1)/4} \approx e^{-1/2}$ にすぎませんが、実際には掃引を重ねると収束は 2 次になります（）。最大の成分を探す手間を省き、$(p, q)$ を決まった順に回る巡回 Jacobi 法が実用の形です。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64）／乱数種 20261006 import numpy as np rng = np.random.default_rng(20261006) n = 8 G = rng.standard_normal((n, n)) A = (G + G.T) / 2 ref = np.linalg.eigvalsh(A) def offd(B): # 非対角部分 return B - np.diag(np.diag(B)) N = n * (n - 1) // 2 # 1 掃引 = N 回の回転 print("掃引 off(A) 比 評価の比 (1-2/(n(n-1)))^(N/2)") B, prev = A.copy(), np.linalg.norm(offd(A)) print(f" 0 {prev:.1e}") for sweep in range(1, 4): for _ in range(N): p, q = np.unravel_index(np.argmax(abs(offd(B))), B.shape) th = 0.5 * np.arctan2(-2 * B[p, q], B[p, p] - B[q, q]) J = np.eye(n) J[p, p] = J[q, q] = np.cos(th) J[p, q], J[q, p] = np.sin(th), -np.sin(th) B = J.T @ B @ J o = np.linalg.norm(offd(B)) print(f"{sweep:3d} {o:.1e} {o / prev:.1e}" f" {(1 - 2 / (n * (n - 1))) ** (N / 2):.2f}") prev = o d = np.sort(np.diag(B)) print("固有値の誤差 < 1e-13:", np.abs(d - ref).max() < 1e-13) ``` ```text 掃引 off(A) 比 評価の比 (1-2/(n(n-1)))^(N/2) 0 5.8e+00 1 2.8e-01 4.8e-02 0.60 2 1.5e-03 5.5e-03 0.60 3 4.2e-09 2.8e-06 0.60 固有値の誤差 < 1e-13: True ``` **検算.** (a) 3 掃引後の対角成分は `numpy.linalg.eigvalsh` の固有値と $10^{-13}$ 以内で一致します。 の保証は $\mathrm{off}(A_3) \approx 4\times10^{-9}$ までですが、実際の誤差はずっと小さく、これは対称行列の固有値の誤差が非対角部分の 2 次で効くためです（ と同じ現象）。(b) 回転はコードの `arctan2` で の $\phi$ をそのまま実装しています。(c) 1 掃引ごとの比は評価の $0.60$ をはるかに下回り、$\mathrm{off}$ は $5.8 \to 0.28 \to 1.5\times10^{-3} \to 4.2\times10^{-9}$ と減り、$\mathrm{off}(A_0)$ に対する比で見ると $4.8\times10^{-2} \to 2.6\times10^{-4} \to 7.3\times10^{-10}$ と、1 掃引で桁数が 2 倍以上になっています（漸近 2 次収束）。 巡回 Jacobi 法は、固有値が相異なれば漸近的に 2 次収束することが知られています（掃引の順序に関する精密な解析が要る）。また、Jacobi 法は小さい固有値も**相対**誤差で正確に求められる場合があり（たとえば $A = DHD$、$D$ が対角で $H$ がよい条件の正定値行列のとき）、QR 法にはない利点です。いずれも本章では使いません。Golub–Van Loan『Matrix Computations』4th ed. §8.5、J. Demmel, K. Veselić, "Jacobi's method is more accurate than QR", *SIAM J. Matrix Anal. Appl.* 13 (1992) を参照してください。 ## 5. 計算した固有値・固有ベクトルの精度 $A$ を実対称行列とし、計算した値 $\hat{\lambda}_1 \ge \cdots \ge \hat{\lambda}_n$ が、ある対称行列 $E$（$\|E\|_2 \le \eta\|A\|_2$）について $A + E$ の厳密な固有値であるとする。このとき $|\hat{\lambda}_i - \lambda_i(A)| \le \eta\|A\|_2$（すべての $i$）である。 $\hat{\lambda}_i = \lambda_i(A + E)$ なので、エルミート行列の固有値の摂動（）から $|\hat{\lambda}_i - \lambda_i(A)| \le \|E\|_2 \le \eta\|A\|_2$ です。 三重対角化・対称 QR 法・Jacobi 法はいずれも直交変換を積み重ねる算法で、後退誤差は $\eta = cnu$（$c$ は小さな定数）の程度です（[QR 分解の章](/mathematics/numerical-linear-algebra/qr-factorization) の Householder 変換の後退安定性と同じ議論。デフレーションの誤差は ）。したがって**どの固有値も絶対誤差 $cnu\|A\|_2$ の程度で求まります**。ただし相対誤差ではないので、$|\lambda_i| \ll \|A\|_2$ の小さい固有値は桁を失いえます。 固有ベクトルは事情が違います。次の評価が示すように、誤差は固有値の隔たりに反比例します。 次の定理の記号は、$\mathcal{V}$ が固有値 $\lambda$ の固有空間（その固有値の固有ベクトル全体が作る [部分空間](/mathematics/linear-algebra/vector-spaces#def-subspace)）、$\delta$ が $\theta$ から見て最も近い「別の固有値」までの距離、$\boldsymbol{r}$ が「$\boldsymbol{x}$ が固有ベクトルからどれだけずれているか」を測る残差です。言いたいことは、残差が小さくても、別の固有値が近いと角度は大きくずれうる、ということです。 $A$ を Hermite 行列、$\boldsymbol{x}$ を単位ベクトル、$\theta = \boldsymbol{x}^{*}A\boldsymbol{x}$、$\boldsymbol{r} = A\boldsymbol{x} - \theta\boldsymbol{x}$ とする。$\lambda$ を $A$ の固有値、$\mathcal{V}$ をその固有空間とし、 $$ \delta = \min\{|\lambda' - \theta| : \lambda' \text{ は } \lambda \text{ と異なる } A \text{ の固有値}\} > 0 $$ とする。このとき $\sin\angle(\boldsymbol{x}, \mathcal{V}) \le \|\boldsymbol{r}\|_2/\delta$ である。とくに $\lambda$ が単純で $\boldsymbol{v}$ がその固有ベクトルなら $\sin\angle(\boldsymbol{x}, \boldsymbol{v}) \le \|\boldsymbol{r}\|_2/\delta$。 正規直交な固有ベクトル $\boldsymbol{u}_1, \ldots, \boldsymbol{u}_n$（固有値 $\lambda_j$）で $\boldsymbol{x} = \sum_jc_j\boldsymbol{u}_j$ と展開します（）。$\mathcal{V}$ は $\lambda_j = \lambda$ となる $\boldsymbol{u}_j$ で張られるので（$\mathcal{V}$ の元を展開すると、$\lambda_j \ne \lambda$ の成分は $0$）、$\mathcal{V}$ への直交射影は $\sum_{\lambda_j=\lambda}c_j\boldsymbol{u}_j$ で、 $$ \sin^2\angle(\boldsymbol{x}, \mathcal{V}) = \sum_{\lambda_j\ne\lambda}|c_j|^2 $$ です（）。一方 $\boldsymbol{r} = \sum_j(\lambda_j - \theta)c_j\boldsymbol{u}_j$ なので $$ \|\boldsymbol{r}\|_2^2 = \sum_j(\lambda_j - \theta)^2|c_j|^2 \ge \sum_{\lambda_j\ne\lambda}(\lambda_j - \theta)^2|c_j|^2 \ge \delta^2\sin^2\angle(\boldsymbol{x}, \mathcal{V}) $$ です。単純な場合は $\mathcal{V} = \operatorname{span}\{\boldsymbol{v}\}$ です。 これは Davis–Kahan の sin θ 定理の 1 次元の場合です。計算した固有ベクトル $\hat{\boldsymbol{x}}$ が $A + E$ の厳密な固有ベクトルなら、$A\hat{\boldsymbol{x}} - \hat\theta\hat{\boldsymbol{x}}$ の大きさは $\|E\|_2$ の程度（ から $\|\boldsymbol{r}\|_2 \le \|(A - \hat\lambda I)\hat{\boldsymbol{x}}\|_2 = \|E\hat{\boldsymbol{x}}\|_2$）なので、固有ベクトルの角の誤差はおよそ $\|E\|_2/\delta$ です。固有値が近い（$\delta$ が小さい）と、固有ベクトルは大きく動きえます。 $\varepsilon > 0$ として $A = \begin{pmatrix} 1 + \varepsilon & 0 \\ 0 & 1 - \varepsilon \end{pmatrix}$、$E = \begin{pmatrix} 0 & \varepsilon \\ \varepsilon & 0 \end{pmatrix}$ とします。$A + E = I + \varepsilon\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$ の固有値は $1 \pm \sqrt{2}\varepsilon$ で、$A$ の固有値 $1 \pm \varepsilon$ からの変化は $(\sqrt{2} - 1)\varepsilon \le \|E\|_2 = \varepsilon$ と小さいままです。ところが $A + E$ の固有ベクトルは $(\cos(\pi/8), \sin(\pi/8))$ などで、$A$ の固有ベクトル $\boldsymbol{e}_1$, $\boldsymbol{e}_2$ から $\varepsilon$ によらず $22.5°$ 回転しています。固有値の隔たり $2\varepsilon$ と摂動の大きさ $\varepsilon$ が同じ程度なので、sin θ 評価は $\sin\angle \le O(1)$ しか保証せず、実際に $O(1)$ の回転が起きています。固有値が重なると固有ベクトルは一意でなくなるので、近い固有値の固有ベクトルを個別に求めるのではなく、それらが張る**部分空間**を求めるのが正しい問題設定です（ の $\mathcal{V}$ を固有値の塊の不変部分空間に替えた評価が Davis–Kahan の定理です）。 次の 2 つのコードで確かめます。1 つ目は `scipy.linalg.eigh` の 3 つの駆動ルーチン（後の表）の残差と直交性、2 つ目は固有値の間隔 $g$ を変えたときの固有ベクトルの誤差です。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64, OpenBLAS）／乱数種 20261006 import numpy as np from scipy.linalg import eigh rng = np.random.default_rng(20261006) u = 2.0**-53 def band(x): # 末尾の桁は BLAS で変わるので桁の帯で示す return "< 1" if x < 1 else "< 100" if x < 100 else ">= 100" print(" n driver 残差/(n u) 直交性/(n u)") for n in (100, 400): G = rng.standard_normal((n, n)) A = (G + G.T) / 2 nA = np.linalg.norm(A, 2) for drv in ("ev", "evd", "evr"): w, V = eigh(A, driver=drv) res = np.linalg.norm(A @ V - V * w, 2) / nA orth = np.linalg.norm(V.T @ V - np.eye(n), 2) print(f"{n:4d} {drv:4s} {band(res / (n * u)):>9s}" f" {band(orth / (n * u)):>9s}") n = 400 # 三重対角 Toeplitz（2, -1） T = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1) k = np.arange(1, n + 1) exact = np.sort(4 * np.sin(k * np.pi / (2 * (n + 1))) ** 2) err = np.abs(eigh(T, eigvals_only=True) - exact).max() print("Toeplitz: 最大誤差/(n u ||T||) < 1:", err / (n * u * 4) < 1) ``` ```text n driver 残差/(n u) 直交性/(n u) 100 ev < 1 < 1 100 evd < 1 < 1 100 evr < 1 < 100 400 ev < 1 < 1 400 evd < 1 < 1 400 evr < 1 < 100 Toeplitz: 最大誤差/(n u ||T||) < 1: True ``` **検算.** (a) 三重対角 Toeplitz 行列の固有値は閉じた形 $4\sin^2(k\pi/(2(n+1)))$（[ノルムの章](/mathematics/numerical-linear-algebra/norms-and-conditioning) の三重対角 Toeplitz 行列の補題）で、誤差は $nu\|T\|_2$ 未満（$\|T\|_2 < 4$）と の予言の範囲にあります。(b) 3 つの駆動ルーチンはいずれも残差 $\|AV - V\Lambda\|_2/\|A\|_2$ が $nu$ 未満で、後退安定です。直交性は ev・evd が $nu$ 未満、evr（MRRR）は $nu$ の数倍〜十数倍で、MRRR が直交性を固有ベクトルの計算の精度で確保する算法であることの反映です（）。(c) $n$ を 4 倍にしても $nu$ で割った値は同じ帯に入り、誤差が $n u$ に比例する程度の評価と整合します。値の末尾は BLAS のスレッド数で変わるので、帯で示しました。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）／乱数種 20261006 import numpy as np rng = np.random.default_rng(20261006) n = 50 Q, _ = np.linalg.qr(rng.standard_normal((n, n))) G = rng.standard_normal((n, n)) E = (G + G.T) / 2 E *= 1e-8 / np.linalg.norm(E, 2) # ||E||_2 = 1e-8 print(" 間隔 g sin(角) ||r||/delta ||E||/g") for g in (1e-2, 1e-4, 1e-6, 1e-7): lam = np.concatenate(([1.0, 1.0 + g], np.linspace(2, 3, n - 2))) A = (Q * lam) @ Q.T w, W = np.linalg.eigh(A + E) x = W[:, 0] # A + E の最小の固有ベクトル v = Q[:, 0] # A の最小の固有ベクトル s = np.linalg.norm(x - (v @ x) * v) th = x @ A @ x r = np.linalg.norm(A @ x - th * x) delta = np.min(np.abs(lam[1:] - th)) # 他の固有値との隔たり print(f" {g:.0e} {s:.1e} {r / delta:.1e} {1e-8 / g:.0e}") ``` ```text 間隔 g sin(角) ||r||/delta ||E||/g 1e-02 2.0e-08 5.4e-07 1e-06 1e-04 1.9e-06 5.4e-05 1e-04 1e-06 1.9e-04 5.4e-03 1e-02 1e-07 1.9e-03 5.4e-02 1e-01 ``` **検算.** (a) $A$ の固有値と固有ベクトルは構成に使った値で、角はその固有ベクトルとの間で測っています。(b) どの行でも $\sin\angle \le \|\boldsymbol{r}\|_2/\delta$ が成り立ち、 の評価の内側にあります（$\|\boldsymbol{r}\|_2/\delta$ は $\|E\|_2/g$ と同じ程度）。(c) 間隔 $g$ を 100 分の 1 にすると、角の誤差も評価もちょうど 100 倍になり、$1/g$ に比例します。摂動 $\|E\|_2 = 10^{-8}$ は一定なので、固有値の誤差は $10^{-8}$ 以下のまま、固有ベクトルの誤差だけが隔たりに反比例して増えています。 ## 6. 特異値分解の計算 ここまでは対称行列の話でした。特異値は一般の長方形の行列にも定まりますが、$A$ とその転置を右上と左下に置いた対称な大きい行列を作れば、対称行列の固有値の話に持ち込めます。次の定理は、その大きい行列の固有値が特異値とその符号違い（と、$m > n$ のときの余りの $0$）であることを示します。 $A \in \mathbb{C}^{m\times n}$（$m \ge n$）の特異値を $\sigma_1 \ge \cdots \ge \sigma_n \ge 0$（$\operatorname{rank}A$ より先は $0$ とする）とする。$(m+n)$ 次の Hermite 行列 $$ C(A) = \begin{pmatrix} O & A \\ A^{*} & O \end{pmatrix} $$ の固有値は $\sigma_1, \ldots, \sigma_n, -\sigma_1, \ldots, -\sigma_n$ と $m - n$ 個の $0$ である。とくに大きい順に $\lambda_i(C(A)) = \sigma_i$（$i \le n$）である。 特異値分解 $A = U\Sigma V^{*}$（）の $U$ の列を $\boldsymbol{u}_1, \ldots, \boldsymbol{u}_m$、$V$ の列を $\boldsymbol{v}_1, \ldots, \boldsymbol{v}_n$ とします。$AV = U\Sigma$、$A^{*}U = V\Sigma^{\mathsf{T}}$ から、$i \le n$ で $A\boldsymbol{v}_i = \sigma_i\boldsymbol{u}_i$、$A^{*}\boldsymbol{u}_i = \sigma_i\boldsymbol{v}_i$、$i > n$ で $A^{*}\boldsymbol{u}_i = \boldsymbol{0}$ です（$\Sigma^{\mathsf{T}}$ の第 $i$ 列は $i > n$ で $0$）。よって $i \le n$ について $$ C(A)\begin{pmatrix} \boldsymbol{u}_i \\ \pm\boldsymbol{v}_i \end{pmatrix} = \begin{pmatrix} \pm A\boldsymbol{v}_i \\ A^{*}\boldsymbol{u}_i \end{pmatrix} = \pm\sigma_i\begin{pmatrix} \boldsymbol{u}_i \\ \pm\boldsymbol{v}_i \end{pmatrix} $$ （複号同順）で、$i > n$ では $C(A)(\boldsymbol{u}_i; \boldsymbol{0}) = \boldsymbol{0}$ です。これら $2n + (m - n) = m + n$ 個のベクトルを $\sqrt{2}$ や $1$ で正規化すると、$\boldsymbol{u}_i$ たちと $\boldsymbol{v}_i$ たちがそれぞれ正規直交であることから互いに正規直交です（たとえば $(\boldsymbol{u}_i; \boldsymbol{v}_i)$ と $(\boldsymbol{u}_i; -\boldsymbol{v}_i)$ の内積は $1 - 1 = 0$）。$m + n$ 次元の正規直交な固有ベクトルの基底が得られたので、固有値は上のもので尽きます。$\sigma_n \ge 0$ なので、大きい順の先頭 $n$ 個は $\sigma_1, \ldots, \sigma_n$ です。 $A, E \in \mathbb{C}^{m\times n}$（$m \ge n$）について、すべての $i \le n$ で $|\sigma_i(A + E) - \sigma_i(A)| \le \|E\|_2$ である。 $C(A + E) = C(A) + C(E)$ です。Hermite 行列 $H$ の 2 ノルムは固有値の絶対値の最大です（固有ベクトルの正規直交基底で展開すると $\|H\boldsymbol{x}\|_2^2 = \sum_j\lambda_j^2|c_j|^2$）。 から $C(E)$ の固有値の絶対値の最大は $\sigma_1(E)$ で、これは $\|E\|_2$ に等しい（ の 3）ので $\|C(E)\|_2 = \|E\|_2$ です。 の後半とエルミート行列の固有値の摂動（）から $|\sigma_i(A + E) - \sigma_i(A)| = |\lambda_i(C(A + E)) - \lambda_i(C(A))| \le \|C(E)\|_2 = \|E\|_2$ です。 特異値も固有値と同じく条件数 1 の問題で、後退安定な算法なら各特異値の絶対誤差は $cnu\|A\|_2 = cnu\,\sigma_1$ の程度です。ところが「$A^{\mathsf{T}}A$ の固有値の平方根」として計算すると、こうはなりません。 $A^{\mathsf{T}}A$ を浮動小数で作って固有値を計算すると、 の誤差は $\|A^{\mathsf{T}}A\|_2 = \sigma_1^2$ に対して $u$ の程度、つまり固有値 $\sigma_i^2$ の絶対誤差が $u\sigma_1^2$ の程度です。$\sigma_i$ に戻すと誤差は約 $u\sigma_1^2/(2\sigma_i)$ で、$\sigma_i \lesssim \sqrt{u}\,\sigma_1 \approx 10^{-8}\sigma_1$ の特異値は 1 桁も正しくありません。直接 SVD を計算すれば誤差は $u\sigma_1$ の程度なので、$\sigma_i \gg u\sigma_1$ である限り何桁かは正しく求まります。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）／乱数種 20261006 import numpy as np rng = np.random.default_rng(20261006) m, n = 60, 13 sig = 10.0 ** -np.arange(n) # 1, 1e-1, ..., 1e-12 U, _ = np.linalg.qr(rng.standard_normal((m, n))) V, _ = np.linalg.qr(rng.standard_normal((n, n))) A = (U * sig) @ V.T s1 = np.linalg.svd(A, compute_uv=False) ev = np.linalg.eigvalsh(A.T @ A)[::-1] s2 = np.sqrt(np.maximum(ev, 0.0)) print(" sigma svd の相対誤差 sqrt(eig(A^T A)) の相対誤差") for i in (0, 3, 6, 8, 10, 12): e1 = abs(s1[i] - sig[i]) / sig[i] e2 = abs(s2[i] - sig[i]) / sig[i] t1 = "< 1e-6" if e1 < 1e-6 else f"{e1:.0e}" t2 = "< 1e-6" if e2 < 1e-6 else f"{e2:.0e}" print(f" 1e-{i:<2d} {t1:>8s} {t2:>8s}") ``` ```text sigma svd の相対誤差 sqrt(eig(A^T A)) の相対誤差 1e-0 < 1e-6 < 1e-6 1e-3 < 1e-6 < 1e-6 1e-6 < 1e-6 < 1e-6 1e-8 < 1e-6 7e-04 1e-10 < 1e-6 2e+01 1e-12 2e-06 1e+00 ``` **検算.** (a) 参照値は構成に使った特異値 $10^{-i}$ です（格納した $A$ 自体が丸めで $u\sigma_1$ 程度ずれるので、$10^{-12}$ では SVD でも相対誤差 $10^{-6}$ の程度が残ります。これは の $\|E\|_2 \approx u$ による $u/\sigma_i \approx 10^{-4}$ 以下という予言の範囲です）。(b) 同じ行列について、SVD（gesdd）と $A^{\mathsf{T}}A$ の固有値という 2 つの方法を比べています。(c) $A^{\mathsf{T}}A$ 経由の誤差は $\sigma = 10^{-8} \approx \sqrt{u}$ で目に見えるようになり、それより小さい特異値では $1$ 以上（桁が残らない）になります。 の $u\sigma_1^2/\sigma_i^2$ の相対誤差の見積もりは上界で、実際の値はそれ以下です。 ### 6.1. 二重対角化と Golub–Kahan の三重対角行列 SVD の計算は、$A^{\mathsf{T}}A$ を作らずに、$A$ に左右から別々の直交変換を掛けて行います。 $A \in \mathbb{R}^{m\times n}$（$m \ge n \ge 2$）について、Householder 鏡映（または単位行列）の積である直交行列 $U \in \mathbb{R}^{m\times m}$、$V \in \mathbb{R}^{n\times n}$ があって $$ U^{\mathsf{T}}AV = \begin{pmatrix} B \\ O \end{pmatrix},\qquad B = \begin{pmatrix} d_1 & e_1 & & \\ & d_2 & \ddots & \\ & & \ddots & e_{n-1} \\ & & & d_n \end{pmatrix} $$ （$B$ は上二重対角）となる。$A$ と $B$ の特異値は等しい。演算量は約 $4mn^2 - 4n^3/3$ である。 $k = 1, \ldots, n$ の順に、左から鏡映 $U_k$ を掛けて第 $k$ 列の第 $k+1$ 行以下を $0$ にし、$k \le n - 2$ なら右から鏡映 $V_k$ を掛けて第 $k$ 行の第 $k+2$ 列以降を $0$ にします。$U_k = \operatorname{diag}(I_{k-1}, \tilde{U}_k)$ は第 $k$ 列の第 $k$ 行以下のベクトルの Householder 鏡映（。ベクトルが $\boldsymbol{0}$ なら単位行列）、$V_k = \operatorname{diag}(I_k, \tilde{V}_k)$ は第 $k$ 行の第 $k+1$ 列以降のベクトル（の転置）の Householder 鏡映です。 $U_k$ は第 $k$ 行以下だけを混ぜるので、すでに作った第 $1, \ldots, k-1$ 行（第 $k-1$ 行は第 $k-1$, $k$ 列以外が $0$）と、第 $1, \ldots, k-1$ 列の第 $k$ 行以下の $0$ を壊しません。$V_k$ は第 $k+1$ 列以降だけを混ぜるので、第 $1, \ldots, k$ 列（とくに今作った第 $k$ 列の $0$）を壊しません。帰納法で、最後に得られる行列は上二重対角の $B$ の下に $O$ を置いた形です。$U = U_1\cdots U_n$、$V = V_1\cdots V_{n-2}$ とします。 特異値が等しいのは、$A = U\begin{pmatrix} B \\ O \end{pmatrix}V^{\mathsf{T}}$ に $B$ の特異値分解 $B = U_B\Sigma_BV_B^{\mathsf{T}}$ を代入すると、$A = \bigl(U\operatorname{diag}(U_B, I)\bigr)\begin{pmatrix} \Sigma_B \\ O \end{pmatrix}(VV_B)^{\mathsf{T}}$ が $A$ の特異値分解になり、$\Sigma$ の対角は特異値分解で一意に定まる（$A^{\mathsf{T}}A$ の固有値の平方根、）からです。演算量は、$p = m - k$, $q = n - k$ として左の鏡映が約 $4pq$、右の鏡映が約 $4pq$ で、$\sum_k 8(m - k)(n - k) \approx 8(mn^2/2 - n^3/6) = 4mn^2 - 4n^3/3$ です。 $m \gg n$ のときは、先に QR 分解 $A = QR$ をして $n$ 次の $R$ を二重対角化する方が速くなります（R-二重対角化）。 $B$ を の $n$ 次上二重対角行列とする。$2n$ 次の対称三重対角行列 $$ T_{\mathrm{GK}} = \begin{pmatrix} 0 & d_1 & & & & \\ d_1 & 0 & e_1 & & & \\ & e_1 & 0 & d_2 & & \\ & & d_2 & 0 & \ddots & \\ & & & \ddots & \ddots & d_n \\ & & & & d_n & 0 \end{pmatrix} $$ （対角は $0$、副対角は $d_1, e_1, d_2, e_2, \ldots, e_{n-1}, d_n$）は、Jordan–Wielandt 行列 $C(B)$ を置換行列で並べ替えたもので、その固有値は $\pm\sigma_1(B), \ldots, \pm\sigma_n(B)$ である。 $C(B)\begin{pmatrix} \boldsymbol{x} \\ \boldsymbol{y} \end{pmatrix} = \begin{pmatrix} B\boldsymbol{y} \\ B^{\mathsf{T}}\boldsymbol{x} \end{pmatrix}$ の成分は $$ (B\boldsymbol{y})_i = d_iy_i + e_iy_{i+1},\qquad (B^{\mathsf{T}}\boldsymbol{x})_j = e_{j-1}x_{j-1} + d_jx_j $$ です（$e_0 = e_n = 0$）。未知数を $\boldsymbol{z} = (y_1, x_1, y_2, x_2, \ldots, y_n, x_n)$ の順に並べ替える置換行列を $P$ とすると、$y_j$（第 $2j - 1$ 位置）の式は $x_{j-1}$（第 $2j - 2$ 位置）と $x_j$（第 $2j$ 位置）だけを係数 $e_{j-1}$, $d_j$ で含み、$x_i$（第 $2i$ 位置）の式は $y_i$（第 $2i - 1$ 位置）と $y_{i+1}$（第 $2i + 1$ 位置）だけを係数 $d_i$, $e_i$ で含みます。よって $P^{\mathsf{T}}C(B)P = T_{\mathrm{GK}}$ です。置換行列は直交なので固有値は $C(B)$ と同じで、$m = n$ の から $\pm\sigma_i(B)$ です。 したがって、二重対角行列の特異値は、対角が $0$ の対称三重対角行列の固有値問題として、§3 の二分法でも、[QR 法の章](/mathematics/numerical-linear-algebra/qr-algorithm) のシフトつき QR 法でも求められます。実際の実装（Golub–Kahan–Reinsch の算法）は $T_{\mathrm{GK}}$ を作らず、$B$ に左右から Givens 回転を掛ける「暗黙のシフトつき QR 法」を $B$ の上で直接行います。これは $B^{\mathsf{T}}B$ に対するシフトつき QR 法と同値で、陰的 Q の定理（）により正当化されます。 ## 7. LAPACK と NumPy・SciPy の対応 | 関数 | LAPACK の駆動ルーチン | 算法 | |---|---|---| | `numpy.linalg.eigh`・`eigvalsh` | syevd（heevd） | 三重対角化＋分割統治法 | | `scipy.linalg.eigh(..., driver="ev")` | syev | 三重対角化＋暗黙のシフトつき QR 法 | | `scipy.linalg.eigh(..., driver="evd")` | syevd | 三重対角化＋分割統治法 | | `scipy.linalg.eigh(..., driver="evr")`（標準問題の既定） | syevr | 三重対角化＋MRRR（固有値の一部も指定可） | | `scipy.linalg.eigh(..., driver="evx")` | syevx | 三重対角化＋二分法と逆反復（固有値の一部） | | `scipy.linalg.eigh_tridiagonal(..., lapack_driver="stebz")` | stebz | Sturm 列の二分法（§3） | | `numpy.linalg.svd`、`scipy.linalg.svd`（既定） | gesdd | 二重対角化＋分割統治法 | | `scipy.linalg.svd(..., lapack_driver="gesvd")` | gesvd | 二重対角化＋暗黙のシフトつき QR 法 | どの駆動ルーチンも三重対角化（二重対角化）で後退安定なので、固有値・特異値の絶対誤差は $nu\|A\|_2$ の程度で、違いは速さと固有ベクトルの直交性の細部です。全部の固有ベクトルが要る大きな問題では分割統治法・MRRR が速く、少数の固有値だけなら evr・evx が手間を省けます。大きな疎行列の一部の固有値には、Krylov 部分空間法の章の `scipy.sparse.linalg.eigsh` を使います。 syev の QR 法は Wilkinson シフト（[QR 法の章](/mathematics/numerical-linear-algebra/qr-algorithm) の注意）を使い、すべての対称三重対角行列で大域的に収束します。その証明は Parlett『The Symmetric Eigenvalue Problem』第 8 章にあり、本章では使いません。局所的な 3 次収束は で示しました。 **分割統治法**（Cuppen）は、三重対角行列を 2 つに分けて階数 1 の修正で結び、修正後の固有値を「永年方程式」$1 + \rho\sum_i z_i^2/(d_i - \lambda) = 0$ の根として求めます。固有ベクトルの直交性を保つ工夫（Gu–Eisenstat）を含めて、全固有ベクトルを最も速く求める方法の 1 つです。**MRRR**（多重相対的にロバストな表現、Dhillon–Parlett）は、シフトした $LDL^{\mathsf{T}}$ 分解の表現の上で各固有値を高い相対精度で求め、逆反復の変種で固有ベクトルを互いに独立に $O(n)$ ずつ計算して、全体を $O(n^2)$ で行います。どちらも本章の定理には使いません。J. J. M. Cuppen, "A divide and conquer method for the symmetric tridiagonal eigenproblem", *Numer. Math.* 36 (1981)、I. S. Dhillon, B. N. Parlett, "Multiple representations to compute orthogonal eigenvectors of symmetric tridiagonal matrices", *Linear Algebra Appl.* 387 (2004) を参照してください。 大きな行列の上位 $k$ 個の特異値だけが要るとき、ガウス乱数の行列 $\Omega$（$n\times(k + p)$）を掛けた $Y = A\Omega$ の列空間の正規直交基底 $Q$ を作り、小さな行列 $Q^{\mathsf{T}}A$ の SVD から近似する**乱択 SVD** がよく使われます。その誤差が最良の階数 $k$ 近似の誤差の小さな定数倍に収まることは確率論の評価（ガウス行列の擬似逆のモーメント）で示され、本章の範囲を超えます。N. Halko, P. G. Martinsson, J. A. Tropp, "Finding structure with randomness", *SIAM Review* 53 (2011) を参照してください。 ## 8. 演習 $n = 1000$ 次の対称三重対角行列の固有値がすべて $[-4, 4]$ にあるとき、 で 1 つの固有値を幅 $10^{-12}$ の区間に閉じ込めるのに要する反復回数と、おおよその演算量（Sturm 列 1 回を $3n$ 回の演算とする）を求めよ。全固有値を求める場合はどうか。 $\lceil\log_2(8/10^{-12})\rceil = \lceil\log_2 8 + 12\log_2 10\rceil = \lceil 3 + 39.86\rceil = 43$ 回です。1 回が約 $3000$ 演算なので、固有値 1 つあたり約 $1.3\times10^{5}$、全固有値で約 $1.3\times10^{8}$ 演算です（区間の分割を共有すると減らせます）。三重対角化の約 $4n^3/3 \approx 1.3\times10^{9}$ よりは小さく、固有値だけなら主な手間は三重対角化です。 $A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$ に の回転を 1 回施し、$\phi$ と $B = J^{\mathsf{T}}AJ$ を求めよ。$\mathrm{off}(B)^2 = \mathrm{off}(A)^2 - 2a_{12}^2$ を確かめよ。 $a_{11} - a_{22} = 0$、$-2a_{12} = -2$ なので $(\cos2\phi, \sin2\phi) = (0, -1)$、$\phi = -\pi/4$、$c = 1/\sqrt{2}$、$s = -1/\sqrt{2}$ です。$J = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix}$ で、 $$ AJ = \frac{1}{\sqrt{2}}\begin{pmatrix} 3 & -1 \\ 3 & 1 \end{pmatrix},\qquad B = J^{\mathsf{T}}AJ = \frac12\begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}\begin{pmatrix} 3 & -1 \\ 3 & 1 \end{pmatrix} = \begin{pmatrix} 3 & 0 \\ 0 & 1 \end{pmatrix} $$ です。$\mathrm{off}(A)^2 = 2$、$\mathrm{off}(B)^2 = 0 = 2 - 2\cdot1^2$ で成り立ちます。対角成分の 2 乗和は $8$ から $10$ に $2a_{12}^2 = 2$ 増えています。 $\sigma_1 = 1$、$\sigma_n = 10^{-6}$ の行列について、後退誤差 $\|E\|_2 \le 10u$ の SVD と、$A^{\mathsf{T}}A$ の固有値（後退誤差 $\|E'\|_2 \le 10u\|A^{\mathsf{T}}A\|_2$）の平方根で求めた $\sigma_n$ の相対誤差の上界を、 と から見積もれ（$u \approx 1.1\times10^{-16}$）。 SVD: $|\hat{\sigma}_n - \sigma_n| \le 10u \approx 1.1\times10^{-15}$ で、相対誤差は $1.1\times10^{-9}$ 以下です。$A^{\mathsf{T}}A$: 固有値 $\sigma_n^2 = 10^{-12}$ の絶対誤差は $10u\cdot1 \approx 1.1\times10^{-15}$ 以下で、相対誤差は $1.1\times10^{-3}$ 以下、平方根をとると相対誤差は約半分の $5.5\times10^{-4}$ 以下です。SVD は約 9 桁、$A^{\mathsf{T}}A$ 経由は約 3 桁しか保証されません。 ## 参考文献 - B. N. Parlett『The Symmetric Eigenvalue Problem』SIAM, 1998 — 第 3 章（Sturm 列）、第 7 章（三重対角化）、第 8 章（QR 法）、第 9 章（Jacobi 法）、第 10・11 章（固有値・固有ベクトルの誤差評価）。 - G. H. Golub, C. F. Van Loan『Matrix Computations』4th ed., Johns Hopkins University Press, 2013 — §8.1（対称固有値問題の摂動）、§8.3〜8.5（対称 QR 法・分割統治・Jacobi 法）、§8.6（SVD の計算）。 - J. W. Demmel『Applied Numerical Linear Algebra』SIAM, 1997 — 第 5 章（対称固有値問題と SVD）。 - G. H. Golub, W. Kahan, "Calculating the singular values and pseudo-inverse of a matrix", *SIAM J. Numer. Anal.* 2 (1965), 205–224. - E. Anderson ほか『LAPACK Users' Guide』3rd ed., SIAM, 1999 — 駆動ルーチンの一覧。