# 固有値の計算 II：QR 法の収束・Hessenberg 化・陰的 Q の定理とシフト > QR 反復が直交相似変換であり同時反復と同値であることを示し、シフトなし QR 法の収束、Householder 変換による Hessenberg 化、陰的 Q の定理、Rayleigh 商シフトと Rayleigh 商反復の関係による局所 3 次収束、デフレーションの後退誤差を証明する。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/qr-algorithm ## 0. この記事の要点 行列の固有値（行列を掛けても向きが変わらないベクトルに対する倍率）を、計算機でまとめて求めるにはどうすればよいか、という話です。ここで扱う QR 法は、「行列を 2 つの行列の積に分け、順序を入れ替えて掛け直す」という操作を繰り返すだけの方法で、続けるうちに対角に固有値が並びます。なぜこの単純な操作で固有値が出るのか、どうすれば速く正確にできるのかを、この記事で順に確かめます。読み終えると、numpy などの固有値計算の中身で何が起きているかが説明できるようになります。 - **QR 法**は、$A_{k-1} = Q_kR_k$ と分解して $A_k = R_kQ_k$ と掛け直すだけの反復です。$A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k$ なので固有値は変わらず、$A_k$ は上三角行列（Schur 形。対角より下がすべて $0$ で、対角に固有値が並ぶ形）に近づいていきます。 - 掛け合わせた直交行列 $\mathbf{Q}_k = Q_1\cdots Q_k$ は $A^k$ の QR 分解の $Q$ で、単位行列から始めた**[同時反復](/mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-subspace-iteration)**そのものです。このため冪乗法の章の同時反復の収束がそのまま QR 法の収束になります。 - $|\lambda_1| > \cdots > |\lambda_n| > 0$ で固有ベクトル行列の逆が LU 分解をもてば、$A_k$ の対角より下は $O(r^k)$（$r = \max_j|\lambda_{j+1}/\lambda_j|$）で $0$ に、対角は固有値に収束します。 - 実用の QR 法は、最初に Householder 変換で **Hessenberg 形**（対角のすぐ下の斜めの列までにしか成分がない形）に直し（$10n^3/3$ 回の演算）、Hessenberg 形を保ったまま 1 回 $O(n^2)$ で反復します。**陰的 Q の定理**により、Hessenberg 形への直交相似は第 1 列だけで決まります。 - シフト $\mu_k = (A_{k-1})_{nn}$ をとると、$\mathbf{Q}_k$ の最後の列は Rayleigh 商反復そのものになり、対称行列では最下行の非対角成分が 1 回ごとにほぼ 3 乗で $0$ に近づきます。$|h_{k+1,k}|$ を単位丸め（計算機の丸め誤差の相対的な大きさ）程度で $0$ とみなすデフレーション（そこで行列を小さい 2 つに切り分けること）は、$\|A\|$ に対して単位丸め程度の後退誤差（答えを、元の問題をほんの少し変えた問題の厳密な答えとみなしたときの、その変化の大きさ）です。 ## 1. 動機 [冪乗法・逆反復の章](/mathematics/numerical-linear-algebra/power-and-inverse-iteration) の方法は、固有値を 1 つずつ、あるいは上位の数個を求めるものでした。密行列のすべての固有値を求める標準的な算法は QR 法で、`numpy.linalg.eig`・`scipy.linalg.eigvals`（LAPACK の geev）も、対称行列の `eigh`（syev 系）も、Hessenberg 化（対称なら三重対角化）のあとに QR 法か、同じ三重対角形の上に立つ後継の算法（分割統治法など）を使っています。 QR 法は 1961 年に Francis と Kublanovskaya が独立に発表しました。算法の記述は 1 行ですが、なぜ収束するのかは一見わかりません。本章の主題は、QR 法が「すべての固有ベクトルについて同時に冪乗法を行い、互いに直交させ続ける」算法であることを示し、そこから収束・計算の工夫（Hessenberg 化）・加速（シフト）を順に導くことです。 この章では行列はすべて実行列とし、直交行列・転置で書きます。[QR 分解の章](/mathematics/numerical-linear-algebra/qr-factorization) の一意性（ の 2。列が一次独立なら、$R$ の対角を正にした QR 分解はただ 1 つ）を何度も使います。複素行列でもユニタリ行列に置き換えれば同じ議論が成り立ちます。ノルムは 2 ノルム $\|\cdot\|_2$ と [Frobenius ノルム](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-frobenius) $\|\cdot\|_F$ を使い、$\|BC\|_F \le \|B\|_2\|C\|_F$（$C$ の各列に $B$ を掛けて評価すればよい）、それを転置した $\|CB\|_F \le \|C\|_F\|B\|_2$（$\|B^{\mathsf{T}}\|_2 = \|B\|_2$、 の 3）と、直交行列 $Q$ について $\|QC\|_F = \|C\|_F$、$\|QC\|_2 = \|C\|_2$ を断りなく使います。 ## 2. QR 反復と同時反復 $A \in \mathbb{R}^{n\times n}$ を正則とする。$A_0 = A$ とし、$k = 1, 2, \ldots$ について、$A_{k-1}$ の QR 分解 $A_{k-1} = Q_kR_k$（$Q_k$ は直交行列、$R_k$ は対角成分が正の上三角行列）をとり、 $$ A_k = R_kQ_k $$ とする。$\mathbf{Q}_k = Q_1Q_2\cdots Q_k$、$\mathbf{R}_k = R_kR_{k-1}\cdots R_1$ とおく（$\mathbf{Q}_0 = \mathbf{R}_0 = I$）。 $A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k$ である。とくに $A_k$ はすべて正則で、$A$ と同じ固有値（重複度を込めて）をもつ。 $R_k = Q_k^{\mathsf{T}}A_{k-1}$ なので $A_k = R_kQ_k = Q_k^{\mathsf{T}}A_{k-1}Q_k$ です。これを繰り返せば $A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k$ で、$\mathbf{Q}_k$ は直交行列の積なので直交行列、とくに $\mathbf{Q}_k^{\mathsf{T}} = \mathbf{Q}_k^{-1}$ です。相似な行列の特性多項式は等しい（積の定理 ）ので固有値も等しく、$\det A_k = \det A \ne 0$ です。$A_k$ が正則なので、次の QR 分解も一意に定まります（ の 2）。 の記号で、すべての $k$ について次が成り立つ。 1. $A^k = \mathbf{Q}_k\mathbf{R}_k$ であり、これは $A^k$ の（$R$ の対角を正にした）QR 分解である。 2. $Z_0 = I$ から始めた同時反復（）の $Z_k$ は $\mathbf{Q}_k$ に等しい。 1. $k$ についての帰納法です。$k = 0$ は明らかです。$A^{k-1} = \mathbf{Q}_{k-1}\mathbf{R}_{k-1}$ とすると、 から $A\mathbf{Q}_{k-1} = \mathbf{Q}_{k-1}A_{k-1}$ なので $$ A^k = A\mathbf{Q}_{k-1}\mathbf{R}_{k-1} = \mathbf{Q}_{k-1}A_{k-1}\mathbf{R}_{k-1} = \mathbf{Q}_{k-1}Q_kR_k\mathbf{R}_{k-1} = \mathbf{Q}_k\mathbf{R}_k $$ です。$\mathbf{R}_k$ は対角成分が正の上三角行列の積なので、対角成分が正の上三角行列です。 2. 帰納法で $Z_{k-1} = \mathbf{Q}_{k-1}$ とします。同時反復は $AZ_{k-1}$ の列を Gram–Schmidt の直交化で正規直交化し、$AZ_{k-1} = Z_kR'_k$（$R'_k$ は対角成分が正の上三角）と書くものでした。一方 $A\mathbf{Q}_{k-1} = \mathbf{Q}_{k-1}A_{k-1} = \mathbf{Q}_{k-1}Q_kR_k = \mathbf{Q}_kR_k$ も同じ行列の、$R$ の対角が正の QR 分解です。$A\mathbf{Q}_{k-1}$ は正則なので、QR 分解の一意性（ の 2）から $Z_k = \mathbf{Q}_k$ です。 QR 法は、$n$ 本の正規直交なベクトル $\mathbf{Q}_{k-1}$ に $A$ を掛けて直交化し直す[同時反復](/mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-subspace-iteration)を、「$A$ を $\mathbf{Q}_k$ の基底で見た行列 $A_k$」の側で記録している、ということです。同時反復では先頭 $j$ 列が優勢な $j$ 次元の不変部分空間に近づきました（）。$\mathbf{Q}_k$ の先頭 $j$ 列が $A$ の不変部分空間をちょうど張れば、$A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k$ の左下の $(n-j)\times j$ ブロックは $0$ になります。これが QR 法で $A_k$ が上三角に近づく理由です。 ## 3. シフトなし QR 法の収束 収束の速さを定量的に示すには、「単位行列に近い行列の QR 分解の $Q$ は単位行列に近い」ことが要ります。次の補題はそのための道具で、$I + F$ の $F$ が小さいほど $Q$ は $I$ に近く、その差が $F$ の大きさのほぼ $\sqrt{n}$ 倍以下で押さえられる、と言っています。証明は、直交行列の各行・各列の 2 乗和が $1$ であることを使い、$Q$ の転置の下三角部分、上三角部分、対角の順に「$I$ との差」を評価していきます。 $F \in \mathbb{R}^{n\times n}$、$s = \|F\|_F$ とし、$(n+1)s^2 < 1$ とする。このとき $I + F$ は正則で、その QR 分解 $I + F = \tilde{Q}\tilde{R}$（$\tilde{R}$ の対角は正）について $$ \|\tilde{Q} - I\|_F \le \sqrt{n}\,s + n s^2 $$ が成り立つ。 $s < 1$ なので、$\|(I + F)\boldsymbol{x}\|_2 \ge \|\boldsymbol{x}\|_2 - \|F\|_2\|\boldsymbol{x}\|_2 \ge (1 - s)\|\boldsymbol{x}\|_2$（$\|F\|_2 \le \|F\|_F$ は Cauchy–Schwarz の不等式 から）で、$I + F$ は正則です。$U = \tilde{Q}^{\mathsf{T}} = (u_{ij})$ とおきます。 **$U$ の狭義下三角部分は小さい.** $\tilde{Q}^{\mathsf{T}}(I + F) = \tilde{R}$ より $U = \tilde{R} - UF$ で、$\tilde{R}$ は上三角なので、$U$ の狭義下三角部分は $-UF$ の狭義下三角部分に等しく、その Frobenius ノルムは $\|UF\|_F = s$ 以下です。 **狭義上三角部分も小さい.** $U$ は直交行列なので、各行・各列の 2 乗和は $1$ です。$1 \le m < n$ について、先頭 $m$ 行の 2 乗和と先頭 $m$ 列の 2 乗和はともに $m$ で、共通の左上ブロック $\sum_{i,j\le m}u_{ij}^2$ を引くと $$ \sum_{i\le mj}u_{ij}^2 \le (n - 1)s^2 $$ です。よって対角を除いた部分を $U - D$（$D$ は $U$ の対角部分）とすると $\|U - D\|_F^2 \le ns^2$ です。 **対角は $1$ に近い.** 第 $j$ 列の 2 乗和が $1$ なので $u_{jj}^2 = 1 - \sum_{i\ne j}u_{ij}^2 \ge 1 - ns^2 > s^2$ です。一方 $U = \tilde{R} - UF$ の対角成分を比べると $u_{jj} = \tilde{r}_{jj} - (UF)_{jj} > -|(UF)_{jj}| \ge -s$ です。$|u_{jj}| > s$ と $u_{jj} > -s$ から $u_{jj} > 0$ で、$0 < u_{jj} \le 1$ です。したがって $1 - u_{jj} \le 1 - u_{jj}^2 = \sum_{i\ne j}u_{ij}^2$ で、$j$ について足して $\|D - I\|_F \le \sum_j(1 - u_{jj}) \le \|U - D\|_F^2 \le ns^2$ です。 以上から $\|\tilde{Q} - I\|_F = \|U - I\|_F \le \|U - D\|_F + \|D - I\|_F \le \sqrt{n}\,s + ns^2$ です。 次の定理が、この節の本題です。$\mathbf{Q}_k$ が固有ベクトルの QR 分解に出てくる直交行列 $Q_X$ に近づくことを、上の補題で見積もり、そこから $A_k$ の対角より下が $0$ に近づく速さを出します。証明に出る記号は、$X$ が固有ベクトルを並べた行列、$S$ が $A$ の上三角形（Schur 形）、$D_k$ が符号をそろえるだけの $\pm1$ の対角行列で、最後の式 $A_k = D_k\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_kD_k$ は「$\tilde{Q}_k \approx I$ なら $A_k \approx D_kSD_k$（$S$ の成分に符号が掛かっただけ）」を言うための形です。 $A \in \mathbb{R}^{n\times n}$ が $A = X\Lambda X^{-1}$、$\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_n)$、 $$ |\lambda_1| > |\lambda_2| > \cdots > |\lambda_n| > 0 $$ と対角化でき、$X^{-1}$ が LU 分解 $X^{-1} = LU$（$L$ は単位下三角、$U$ は対角成分が $0$ でない上三角）をもつとする。$r = \max_{1\le j j),\qquad |(A_k)_{ii} - \lambda_i| \le Cr^k $$ を満たす。とくに $A_k$ の対角より下は $0$ に、対角成分は固有値に、大きさの順に収束する。 **準備.** 実行列の実でない固有値は共役な対で現れ絶対値が等しいので、仮定から固有値はすべて実数です。各固有空間は 1 次元なので、$X$ の各列を零でない複素数倍して実ベクトルにできます。$X' = XD'$（$D'$ は正則な対角行列）とすると $X'^{-1} = D'^{-1}LU = (D'^{-1}LD')(D'^{-1}U)$ も同じ形の LU 分解をもつので、はじめから $X$ は実としてかまいません。このとき $X^{-1}$ は実行列で、その LU 分解の $L$, $U$ も実です（LU 分解をもつので先頭主小行列式は $0$ でなく、 により実数の範囲でも LU 分解をもち、その一意性から両者は一致する）。$X$ の QR 分解を $X = Q_XR_X$（ の 2）とし、$S = R_X\Lambda R_X^{-1}$ とおくと、$S$ は対角成分が $\lambda_1, \ldots, \lambda_n$ の上三角行列で、$Q_X^{\mathsf{T}}AQ_X = R_X\Lambda R_X^{-1} = S$ です（$A$ の実 Schur 形）。 **$A^k$ の分解.** $A^k = X\Lambda^kX^{-1} = X\Lambda^kLU$ で、$\Lambda^kL = (\Lambda^kL\Lambda^{-k})\Lambda^k$ です。$E_k = \Lambda^kL\Lambda^{-k} - I$ は狭義下三角で、その $(i, j)$ 成分は $l_{ij}(\lambda_i/\lambda_j)^k$（$i > j$）なので、$|\lambda_i/\lambda_j| \le r$ から $\|E_k\|_F \le \|L\|_Fr^k$ です。よって $$ A^k = Q_XR_X(I + E_k)\Lambda^kU = Q_X(I + F_k)R_X\Lambda^kU,\qquad F_k = R_XE_kR_X^{-1} $$ で、$\|F_k\|_F \le \|R_X\|_2\|R_X^{-1}\|_2\|L\|_Fr^k =: c_1r^k$ です。 **$\mathbf{Q}_k$ の形.** $(n+1)c_1^2r^{2k} < 1$ となる $k$（$r < 1$ なので十分大きい $k$ はすべて）について、 から $I + F_k = \tilde{Q}_k\tilde{R}_k$、$\|\tilde{Q}_k - I\|_F \le \sqrt{n}c_1r^k + nc_1^2r^{2k}$ です。$\tilde{R}_kR_X\Lambda^kU$ は上三角で、対角成分 $(\tilde{R}_k)_{ii}(R_X)_{ii}\lambda_i^ku_{ii}$ は $0$ でありません。その符号を並べた対角行列を $D_k$（成分は $\pm1$）とすると $$ A^k = (Q_X\tilde{Q}_kD_k)(D_k\tilde{R}_kR_X\Lambda^kU) $$ は、右の因子の対角が正の QR 分解です。 の 1 と QR 分解の一意性から $\mathbf{Q}_k = Q_X\tilde{Q}_kD_k$ です。 **$A_k$ の評価.** から $$ A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k = D_k\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_kD_k $$ です。$\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_k - S = (\tilde{Q}_k - I)^{\mathsf{T}}S\tilde{Q}_k + S(\tilde{Q}_k - I)$ から $$ \|A_k - D_kSD_k\|_F = \|\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_k - S\|_F \le 2\|S\|_2\|\tilde{Q}_k - I\|_F \le 2\|S\|_2(\sqrt{n}c_1 + nc_1^2)r^k $$ です。$D_kSD_k$ は上三角で、対角成分は $\lambda_i$ です。よって十分大きい $k$ で主張の評価が成り立ち、有限個の小さい $k$ は定数 $C$ を大きくとって含めます。 $A_k$ の対角より上の成分は $D_kSD_k$ の成分、すなわち $S$ の成分に符号 $\pm1$ が掛かったものに近づきます。$D_k$ の符号は $\lambda_i^k$ の符号とともに変わりうるので、上三角部分は一般には収束せず符号が振動します。固有値（対角）と、下三角部分が $0$ に近づくことだけが意味をもちます。 の評価は比の最大 $r$ だけで書いていますが、より細かく、第 $j$ 列より左・第 $j$ 行より下のブロックが $O(|\lambda_{j+1}/\lambda_j|^k)$ で、成分 $(i, j)$（$i > j$）が $O(|\lambda_i/\lambda_j|^k)$ で $0$ に近づくことが知られています（ の各 $j$ の評価に対応）。下のコードの表はこの細かい評価と合っています。Watkins『Fundamentals of Matrix Computations』3rd ed., Wiley, 2010 — 第 5・6 章、Golub–Van Loan『Matrix Computations』§7.3 を参照してください。 ```python # Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）／乱数種 20261006 import numpy as np rng = np.random.default_rng(20261006) lam = np.array([6.0, 5.0, 4.0, 3.0, 2.0, 1.0]) X = rng.standard_normal((6, 6)) A = X @ np.diag(lam) @ np.linalg.inv(X) def qr_pos(B): # R の対角を正にそろえた QR 分解 Q, R = np.linalg.qr(B) s = np.sign(np.diag(R)) return Q * s, R * s[:, None] Ak, sub = A.copy(), {} for k in range(1, 101): Q, R = qr_pos(Ak) Ak = R @ Q if k in (40, 60): sub[k] = np.abs(np.diag(Ak, -1)) print(" i |a(i+1,i)|(k=60) 20 回の比 (l(i+1)/l(i))^20") for i in range(5): r = sub[60][i] / sub[40][i] print(f" {i+1} {sub[60][i]:.2e} {r:.3e}" f" {(lam[i+1] / lam[i]) ** 20:.3e}") print("対角（k=100）:", np.round(np.diag(Ak), 10)) ``` ```text i |a(i+1,i)|(k=60) 20 回の比 (l(i+1)/l(i))^20 1 3.71e-06 2.610e-02 2.608e-02 2 2.98e-06 1.153e-02 1.153e-02 3 3.87e-08 3.171e-03 3.171e-03 4 2.21e-11 3.007e-04 3.007e-04 5 6.00e-19 9.537e-07 9.537e-07 対角（k=100）: [6. 5. 4. 3. 2. 1.] ``` **検算.** (a) 固有値は構成に使った $6, 5, \ldots, 1$ で、100 回後の対角は 10 桁に丸めてそれに一致します。(b) 別の算法・別の精度との照合は §4 の最後のコードで行います。(c) 副対角成分 $(i+1, i)$ の 20 回あたりの縮小率は $(\lambda_{i+1}/\lambda_i)^{20}$ と 3〜4 桁一致し、 の成分ごとの速さが見えます。 の $r = 5/6$ は最も遅い成分 $(2, 1)$ の速さです。（$X$ は乱数なので、$X^{-1}$ の LU 分解の存在は確率 1 で成り立ちます。） $A = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$ の固有値は $\pm1$ で絶対値が等しく、 の仮定を満たしません。実際 $A$ は直交行列なので $A = A\cdot I$ が QR 分解で、$A_1 = IA = A$ となり、QR 法は一歩も進みません。 シフト $\mu$（$|\mu| < 1$）を入れて $A - \mu I = QR$、$A_1 = RQ + \mu I$ とします（次節のシフトつき QR 法）。$c = 1 + \mu^2$ とおくと $$ Q = \frac{1}{\sqrt{c}}\begin{pmatrix} -\mu & 1 \\ 1 & \mu \end{pmatrix},\qquad R = \frac{1}{\sqrt{c}}\begin{pmatrix} c & -2\mu \\ 0 & 1 - \mu^2 \end{pmatrix},\qquad A_1 = \frac{1}{c}\begin{pmatrix} -2\mu & 1 - \mu^2 \\ 1 - \mu^2 & 2\mu \end{pmatrix} $$ です（$Q$ の列が正規直交で $QR = A - \mu I$、$R$ の対角が正であることは掛け算で確かめられます）。非対角成分は $1$ から $(1 - \mu^2)/(1 + \mu^2)$ に減ります。たとえば $\mu = 1/2$ なら $0.6$ です。ただし、次節の Rayleigh 商シフト $\mu = a_{22} = 0$ ではこの行列は動きません。Rayleigh 商シフトが失敗するこの種の対称な状況を避けるのが、Wilkinson シフト（）です。 ## 4. Hessenberg 化と陰的 Q の定理 密行列の QR 分解は $O(n^3)$ の演算を要するので、QR 法を $n$ 回以上反復すると全体で $O(n^4)$ になります。そこで最初に一度だけ、QR 反復で保たれる簡単な形に直しておきます。 任意の $A \in \mathbb{R}^{n\times n}$（$n \ge 3$）について、$n - 2$ 個の Householder 鏡映（または単位行列）の積である直交行列 $Q$ があって、$H = Q^{\mathsf{T}}AQ$ は**上 Hessenberg 行列**（$i > j + 1$ で $h_{ij} = 0$）になる。演算量は約 $10n^3/3$ である。$A$ が対称なら $H$ は対称な三重対角行列である。 $k = 1, \ldots, n-2$ の順に、第 $k$ 列の第 $k+2$ 行以下を $0$ にします。第 $k - 1$ 段までで、いまの行列 $A^{(k)}$ の第 $1, \ldots, k-1$ 列は Hessenberg 形になっているとします。$\boldsymbol{x} \in \mathbb{R}^{n-k}$ を第 $k$ 列の第 $k+1$ 行以下とし、$\boldsymbol{x} \ne \boldsymbol{0}$ なら $\boldsymbol{x}$ の Householder 鏡映 $\tilde{P}_k$（）をとります。$\tilde{P}_k\boldsymbol{x} = \alpha\boldsymbol{e}_1$ で、$\tilde{P}_k$ は対称な直交行列です（）。$\boldsymbol{x} = \boldsymbol{0}$ なら $\tilde{P}_k = I$ とします。$P_k = \operatorname{diag}(I_k, \tilde{P}_k)$ とおき、$A^{(k+1)} = P_kA^{(k)}P_k$ とします。 左から $P_k$ を掛けると第 $k+1$ 行以下の行だけが混ざり、第 $1, \ldots, k$ 列のうち第 $k + 1$ 行以下の成分は、第 $1, \ldots, k-1$ 列ではもともと $0$ なので $0$ のまま、第 $k$ 列は $(\alpha, 0, \ldots, 0)^{\mathsf{T}}$ になります。右から $P_k$ を掛けると第 $k+1$ 列以降だけが混ざるので、第 $1, \ldots, k$ 列は変わりません。よって $A^{(k+1)}$ の第 $1, \ldots, k$ 列は Hessenberg 形です。$Q = P_1\cdots P_{n-2}$ とすれば $H = Q^{\mathsf{T}}AQ$ です（$P_k^{\mathsf{T}} = P_k$）。 演算量: $m = n - k$ として、$\tilde{P}_k$ を $m\times(m+1)$ 程度のブロックに左から掛けるのに約 $4m^2$、$n\times m$ のブロックに右から掛けるのに約 $4nm$ の演算が要ります（$\tilde{P}\boldsymbol{y} = \boldsymbol{y} - \beta\boldsymbol{v}(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{y})$ は 1 本あたり約 $4m$）。$\sum_m(4m^2 + 4nm) \approx 4n^3/3 + 2n^3 = 10n^3/3$ です。$A$ が対称なら $H^{\mathsf{T}} = Q^{\mathsf{T}}A^{\mathsf{T}}Q = H$ で、対称な Hessenberg 行列は三重対角です。 $H$ が正則な上 Hessenberg 行列で $H = QR$（$R$ の対角は正）なら、$Q$ と $RQ$ も上 Hessenberg 行列である。この QR 分解は $n - 1$ 個の Givens 回転で $O(n^2)$ の演算で計算でき、$RQ$ の計算も $O(n^2)$ である。$H$ が対称な三重対角行列なら $RQ$ も対称な三重対角行列で、1 回の反復は $O(n)$ の演算で済む。 $Q = HR^{-1}$ で、$R^{-1}$ は上三角です（上三角行列の逆行列は後退代入で列ごとに求まり、上三角になる）。Hessenberg 行列 $H$ の第 $j$ 列は第 $j + 1$ 行までしか成分をもたないので、$HR^{-1}$ の第 $j$ 列（$H$ の第 $1, \ldots, j$ 列の 1 次結合）も第 $j + 1$ 行までしか成分をもちません。よって $Q$ は Hessenberg です。$RQ$ の第 $j$ 列は $R$ の第 $1, \ldots, j+1$ 列の 1 次結合で、上三角の $R$ のそれらの列は第 $j + 1$ 行までしか成分をもたないので、$RQ$ も Hessenberg です。 計算は、$i = 1, \ldots, n-1$ の順に第 $i$・$i+1$ 行に Givens 回転 $G_i$ を掛けて $(i+1, i)$ 成分を $0$ にします（。回転は 2 行だけを変え、それぞれ $O(n)$ の演算）。$G_{n-1}\cdots G_1H = R'$ は上三角で、各行の符号を調整すれば $R$ の対角が正になり、$Q = G_1^{\mathsf{T}}\cdots G_{n-1}^{\mathsf{T}}$（符号の対角行列を含めて）です。$RQ = R'G_1^{\mathsf{T}}\cdots G_{n-1}^{\mathsf{T}}$（符号を含めて）も、右から回転を掛けるごとに 2 列だけが変わるので $O(n^2)$ です。対称な場合は $RQ = Q^{\mathsf{T}}HQ$ が対称かつ Hessenberg なので三重対角で、三重対角行列に回転を掛けるときに変わる成分は定数個なので、各回転が $O(1)$、全体で $O(n)$ です。 Hessenberg 化は一度だけ $O(n^3)$ で、その後の反復は 1 回 $O(n^2)$（対称なら $O(n)$）になります。経験的に固有値 1 つあたり 2〜3 回の反復で済むので、全体の手間は Hessenberg 化と同じ $O(n^3)$ の程度です。 Hessenberg 形に直す直交相似変換は何通りもありそうですが、実は第 1 列を決めると、符号を除いてただ 1 通りしか残りません。次の定理がそれを保証します。「既約」とは、副対角成分 $h_{k+1,k}$ がどれも $0$ でないことです。Hessenberg 形と、対称なときの三重対角形の成分の位置は、次の図のとおりです。
上 Hessenberg 形三重対角形
$A \in \mathbb{R}^{n\times n}$、$Q$, $V$ を直交行列とし、$H = Q^{\mathsf{T}}AQ$ と $G = V^{\mathsf{T}}AV$ がともに上 Hessenberg 行列で、$H$ は**既約**（すべての $k$ で $h_{k+1,k} \ne 0$）とする。$Q$ と $V$ の第 1 列が等しいならば、各成分が $\pm1$ の対角行列 $D$ があって $V = QD$、$G = DHD$ である。 $\boldsymbol{q}_1 = Q\boldsymbol{e}_1 = V\boldsymbol{e}_1$ とし、Krylov 行列 $K = (\boldsymbol{q}_1\ \ A\boldsymbol{q}_1\ \ \cdots\ \ A^{n-1}\boldsymbol{q}_1)$ を考えます。$AQ = QH$ から $A^jQ = QH^j$ なので $$ K = Q\,(\boldsymbol{e}_1\ \ H\boldsymbol{e}_1\ \ \cdots\ \ H^{n-1}\boldsymbol{e}_1) =: QR_H,\qquad K = V\,(\boldsymbol{e}_1\ \ G\boldsymbol{e}_1\ \ \cdots\ \ G^{n-1}\boldsymbol{e}_1) =: VR_G $$ です。Hessenberg 行列を掛けると成分の入る行が 1 つ下に延びるだけなので、帰納法で $H^j\boldsymbol{e}_1$ は第 $j + 1$ 行までしか成分をもたず、その第 $j + 1$ 成分は $h_{21}h_{32}\cdots h_{j+1,j}$ です。よって $R_H$ は上三角で、既約性から対角成分はすべて $0$ でなく、$R_H$ は正則です。同様に $R_G$ も上三角です。$K = QR_H$ は正則なので $R_G = V^{\mathsf{T}}K$ も正則です。 $QR_H = VR_G$ から $V^{\mathsf{T}}Q = R_GR_H^{-1} =: W$ は上三角で、直交行列でもあります。$W^{-1} = W^{\mathsf{T}}$ は上三角行列の逆なので上三角、かつ上三角行列の転置なので下三角で、対角行列です。直交な対角行列の成分は $\pm1$ なので $W = D$、$V^{\mathsf{T}}Q = D$ から $V = QD$（$D^{-1} = D$）です。最後に $G = V^{\mathsf{T}}AV = DQ^{\mathsf{T}}AQD = DHD$ です。 この定理の意味は、「既約な Hessenberg 形への直交相似変換は、第 1 列を決めれば符号を除いて一意」ということです。QR 反復 $H_k = Q_k^{\mathsf{T}}H_{k-1}Q_k$ を行うのに、$Q_k$ を明示的に作らず、第 1 列だけを正しく選んで、あとは Hessenberg 形に戻す（「バルジ追跡」。回転を掛けると Hessenberg 形からはみ出る成分が 1 つできるので、これを次々に斜め下へ送り出して形を回復する操作）だけでよいことになります。形が Hessenberg に戻った時点で、定理により符号を除いて明示的に $Q_k$ を作った場合と同じ行列が得られるからです。Francis の二重シフト QR 法（）や Lanczos 法の陰的再開はこの原理で作られています。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64）／乱数種 20261006 import numpy as np from scipy.linalg import hessenberg rng = np.random.default_rng(20261006) n = 8 lam = np.arange(1.0, n + 1) X = rng.standard_normal((n, n)) cases = {"非対称": X @ np.diag(lam) @ np.linalg.inv(X)} Qs, _ = np.linalg.qr(X) cases["対称"] = Qs @ np.diag(lam) @ Qs.T for name, A in cases.items(): H, P = hessenberg(A, calc_q=True) low = np.abs(np.tril(H, -2)).max() orth = np.linalg.norm(P.T @ P - np.eye(n)) print(f"{name}: 下の帯の外 {low:.0e}, ||P^T P - I|| < 1e-14:", orth < 1e-14) seq = [] for k in range(12): mu = H[-1, -1] # Rayleigh 商シフト Q, R = np.linalg.qr(H - mu * np.eye(n)) H = R @ Q + mu * np.eye(n) seq.append(abs(H[-1, -2])) if seq[-1] < 1e-14: break print(" |h(n,n-1)|:", " ".join(f"{s:.1e}" for s in seq)) ``` ```text 非対称: 下の帯の外 0e+00, ||P^T P - I|| < 1e-14: True |h(n,n-1)|: 4.2e-01 4.8e-01 4.6e-01 2.3e-01 1.8e-02 6.8e-06 2.0e-12 1.6e-25 対称: 下の帯の外 0e+00, ||P^T P - I|| < 1e-14: True |h(n,n-1)|: 7.5e-02 1.0e-04 3.5e-13 1.6e-38 ``` **検算.** (a) `scipy.linalg.hessenberg` の結果は副対角より下が厳密に $0$ で、変換行列の直交性は単位丸めの程度です（）。(b) 下の `mpmath` との照合で、固有値そのものも確かめます。(c) Rayleigh 商シフトの QR 法の $|h_{n,n-1}|$ は、非対称では $1.8\times10^{-2} \to 6.8\times10^{-6} \to 2.0\times10^{-12} \to 1.6\times10^{-25}$ と指数がほぼ 2 倍ずつ（2 次収束）、対称では $1.0\times10^{-4} \to 3.5\times10^{-13} \to 1.6\times10^{-38}$ と指数がほぼ 3 倍ずつ（3 次収束、）進みます。最初の数回は大域的な段階で、収束の次数は現れません。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0／乱数種 20261006 import numpy as np import mpmath as mp from scipy.linalg import eigvals rng = np.random.default_rng(20261006) lam = np.array([6.0, 5.0, 4.0, 3.0, 2.0, 1.0]) X = rng.standard_normal((6, 6)) A = X @ np.diag(lam) @ np.linalg.inv(X) # 丸めで固有値はわずかにずれる Ak = A.copy() # シフトなし QR 法 300 回 for k in range(300): Q, R = np.linalg.qr(Ak) Ak = R @ Q mine = np.sort(np.diag(Ak))[::-1] sp_ = np.sort(eigvals(A).real)[::-1] mp.mp.dps = 40 # 格納された A の厳密な固有値 E, _ = mp.eig(mp.matrix(A.tolist())) ref = sorted([mp.re(e) for e in E], reverse=True) d1 = max(abs(float(r - m)) for r, m in zip(ref, mine)) d2 = max(abs(float(r - s)) for r, s in zip(ref, sp_)) d3 = max(abs(float(r) - l) for r, l in zip(ref, lam)) print("自作 QR 法と mpmath の差の最大 < 1e-12:", d1 < 1e-12) print("eigvals と mpmath の差の最大 < 1e-12:", d2 < 1e-12) print("mpmath と構成値 6,...,1 の差の最大 < 1e-13:", d3 < 1e-13) ``` ```text 自作 QR 法と mpmath の差の最大 < 1e-12: True eigvals と mpmath の差の最大 < 1e-12: True mpmath と構成値 6,...,1 の差の最大 < 1e-13: True ``` **検算.** (a) 構成に使った固有値 $6, \ldots, 1$ と、格納された（丸めを含む）行列 $A$ の 40 桁の固有値の差は $10^{-13}$ 未満です。(b) 自作の QR 法（`numpy.linalg.qr` の符号のままの版。符号は $A_k$ の対角に影響しません）と LAPACK の `eigvals` は、どちらも 40 桁の参照値と $10^{-12}$ 以内で一致します。(c) 300 回の反復では最も遅い比 $(5/6)^{300} \approx 2\times10^{-24}$ まで下三角部分が縮んでおり、残る差は丸めの水準です。 ## 5. シフトと Rayleigh 商反復 シフトなし QR 法の速さ $|\lambda_{j+1}/\lambda_j|$ は行列で決まり、固有値が近いと非常に遅くなります。逆反復と同じく、シフト（行列から $\mu I$ を引いてから分解すること）で改善します。 $A_0 = A$ とし、シフト $\mu_k \in \mathbb{R}$ を選んで、$A_{k-1} - \mu_kI$ が正則なら $A_{k-1} - \mu_kI = Q_kR_k$（$R_k$ の対角は正）と分解し、 $$ A_k = R_kQ_k + \mu_kI $$ とする。$\mu_k = (A_{k-1})_{nn}$ とするものを **Rayleigh 商シフト**という。 $A_k = Q_k^{\mathsf{T}}(A_{k-1} - \mu_kI)Q_k + \mu_kI = Q_k^{\mathsf{T}}A_{k-1}Q_k$ なので、シフトつきでも直交相似変換で、固有値は保たれます。$A_{k-1}$ が Hessenberg なら $A_{k-1} - \mu_kI$ も Hessenberg なので、 により $A_k$ も Hessenberg です。 次の定理で転置が出てくるのは、$B = QR$ の逆行列を取ると $Q$ が $B$ の転置の逆行列を使って書けるためです。$A_{k-1}$ が対称なら $A_{k-1}^{\mathsf{T}} = A_{k-1}$ なので、転置は気にせず普通の逆反復と読んでかまいません。 シフトつき QR 法の 1 段で $A_{k-1} - \mu_kI$ が正則とし、$\boldsymbol{q} = Q_k\boldsymbol{e}_n$（$Q_k$ の最後の列）とする。このとき 1. $\boldsymbol{q} = (R_k)_{nn}\,(A_{k-1} - \mu_kI)^{-\mathsf{T}}\boldsymbol{e}_n$ であり、$(R_k)_{nn} > 0$。すなわち $\boldsymbol{q}$ は、$A_{k-1}^{\mathsf{T}}$ にシフト $\mu_k$ の逆反復を $\boldsymbol{e}_n$ から 1 回行ったベクトルである。 2. $\theta = \boldsymbol{q}^{\mathsf{T}}A_{k-1}\boldsymbol{q}$ とおくと $(A_k)_{nn} = \theta$ であり、$A_k$ の最下行の非対角部分について $$ \Bigl(\sum_{j 1. $B = A_{k-1} - \mu_kI = Q_kR_k$ を転置すると $B^{\mathsf{T}} = R_k^{\mathsf{T}}Q_k^{\mathsf{T}}$ で、その逆行列は $B^{-\mathsf{T}} = Q_k^{-\mathsf{T}}R_k^{-\mathsf{T}} = Q_kR_k^{-\mathsf{T}}$（$Q_k$ は直交）なので $Q_k = B^{-\mathsf{T}}R_k^{\mathsf{T}}$ です。$R_k^{\mathsf{T}}$ は下三角なので $R_k^{\mathsf{T}}\boldsymbol{e}_n = (R_k)_{nn}\boldsymbol{e}_n$ で、$Q_k\boldsymbol{e}_n = (R_k)_{nn}B^{-\mathsf{T}}\boldsymbol{e}_n$ です。$(R_k)_{nn} > 0$ は分解の約束です。 2. $A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k$ の $(n, j)$ 成分は $\boldsymbol{q}^{\mathsf{T}}A_{k-1}\boldsymbol{q}_j = (A_{k-1}^{\mathsf{T}}\boldsymbol{q})^{\mathsf{T}}\boldsymbol{q}_j$（$\boldsymbol{q}_j$ は $Q_k$ の第 $j$ 列）で、$j = n$ では $\theta$ です。$\boldsymbol{s} = A_{k-1}^{\mathsf{T}}\boldsymbol{q} - \theta\boldsymbol{q}$ とおくと、$\theta$ は単位ベクトル $\boldsymbol{q}$ の $A_{k-1}^{\mathsf{T}}$ に関する Rayleigh 商なので $\boldsymbol{s} \perp \boldsymbol{q}$ です（）。$j < n$ では $\boldsymbol{q}^{\mathsf{T}}\boldsymbol{q}_j = 0$ なので $(A_k)_{nj} = \boldsymbol{s}^{\mathsf{T}}\boldsymbol{q}_j$ です。$\boldsymbol{q}_1, \ldots, \boldsymbol{q}_n$ は正規直交基底で $\boldsymbol{s}^{\mathsf{T}}\boldsymbol{q}_n = 0$ なので、$\|\boldsymbol{s}\|_2^2 = \sum_{j つまり、シフトつき QR 法は「最後の列について逆反復を行い、最下行の非対角部分にその残差を記録する」算法です。Rayleigh 商シフト $\mu_k = (A_{k-1})_{nn} = \boldsymbol{e}_n^{\mathsf{T}}A_{k-1}\boldsymbol{e}_n$ は $\boldsymbol{e}_n$ の Rayleigh 商なので、これは [Rayleigh 商反復](/mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-rqi)の 1 段になります。対称行列ではこれを元の座標に戻すと、正確に Rayleigh 商反復の列になります。 次の系は、対称行列では最下行の非対角部分が 1 回ごとに前の 3 乗にほぼ比例して小さくなる、という主張です。式の定数 $32\sqrt{2}\,\Delta^2/\delta^4$ の細かい値より、「次の $\eta_{k+1}$ が今の $\eta_k$ の 3 乗に比例する」ことが要点です。$\delta$, $\Delta$ は固有値の間隔と広がりを表す定数で、定義は冪乗法の章の定理にあります。 $A$ を実対称行列、$\lambda$ をその単純な固有値、$\boldsymbol{v}$ を単位固有ベクトルとし、$\delta$, $\Delta$ を と同じく定める。Rayleigh 商シフトの QR 法で、$\boldsymbol{x}_k = \mathbf{Q}_k\boldsymbol{e}_n$（$\mathbf{Q}_k = Q_1\cdots Q_k$）、$\theta_k = \angle(\boldsymbol{x}_k, \boldsymbol{v})$、$\eta_k = \bigl(\sum_{j）の $k$ 回目に等しく、$\eta_k = \|A\boldsymbol{x}_k - R_A(\boldsymbol{x}_k)\boldsymbol{x}_k\|_2$ である。 2. $\Delta\sin^2\theta_0 < \delta/2$ ならば、反復が続く限り $\tfrac{\delta}{2}\sin\theta_k \le \eta_k \le \Delta\sin\theta_k$ であり、 $$ \eta_{k+1} \le \frac{32\sqrt{2}\,\Delta^2}{\delta^4}\,\eta_k^3 . $$ 1. 帰納法です。$k = 0$ では $\mathbf{Q}_0 = I$ です。$\boldsymbol{x}_{k-1} = \mathbf{Q}_{k-1}\boldsymbol{e}_n$ が Rayleigh 商反復の $k - 1$ 回目とすると、$A_{k-1} = \mathbf{Q}_{k-1}^{\mathsf{T}}A\mathbf{Q}_{k-1}$ から $\mu_k = \boldsymbol{e}_n^{\mathsf{T}}A_{k-1}\boldsymbol{e}_n = R_A(\boldsymbol{x}_{k-1})$ で、これは Rayleigh 商反復のシフトそのものです。$A_{k-1}$ は対称なので、 の 1 から $$ \boldsymbol{x}_k = \mathbf{Q}_{k-1}Q_k\boldsymbol{e}_n = (R_k)_{nn}\mathbf{Q}_{k-1}(A_{k-1} - \mu_kI)^{-1}\boldsymbol{e}_n = (R_k)_{nn}(A - \mu_kI)^{-1}\boldsymbol{x}_{k-1} $$ です（$\mathbf{Q}_{k-1}(\mathbf{Q}_{k-1}^{\mathsf{T}}(A - \mu_kI)\mathbf{Q}_{k-1})^{-1} = (A - \mu_kI)^{-1}\mathbf{Q}_{k-1}$）。$\boldsymbol{x}_k$ は単位ベクトルで $(R_k)_{nn} > 0$ なので、Rayleigh 商反復の正規化と一致します。後半は、 の 2 の右辺が、$\boldsymbol{q} = Q_k\boldsymbol{e}_n$ を $\mathbf{Q}_{k-1}$ で元の座標に戻すと $\|A\boldsymbol{x}_k - R_A(\boldsymbol{x}_k)\boldsymbol{x}_k\|_2$ になること（$\mathbf{Q}_{k-1}$ は直交）から従います。 2. $A$ の正規直交な固有ベクトル $\boldsymbol{u}_i$（固有値 $\lambda_i$。$\boldsymbol{v}$ は $\pm\boldsymbol{u}_{i_0}$）で $\boldsymbol{x}_k = \sum_ic_i\boldsymbol{u}_i$ と展開し、$\rho = R_A(\boldsymbol{x}_k)$ とおきます。 から、反復が続く限り $\Delta\sin^2\theta_k < \delta/2$ と $|\rho - \lambda| \le \Delta\sin^2\theta_k < \delta/2$ が成り立ちます。上からは、$\eta_k$ が残差の最小値であること（）から $$ \eta_k^2 \le \|(A - \lambda I)\boldsymbol{x}_k\|_2^2 = \sum_{i}(\lambda_i - \lambda)^2c_i^2 \le \Delta^2\sum_{i\ne i_0}c_i^2 = \Delta^2\sin^2\theta_k $$ です。下からは、$i \ne i_0$ で $|\lambda_i - \rho| \ge \delta - |\lambda - \rho| > \delta/2$ なので $$ \eta_k^2 = \sum_i(\lambda_i - \rho)^2c_i^2 \ge \frac{\delta^2}{4}\sum_{i\ne i_0}c_i^2 = \frac{\delta^2}{4}\sin^2\theta_k $$ です。最後に、$\sin^2\theta_k < \delta/(2\Delta) \le 1/2$ から $\cos^2\theta_k > 1/2$、$\tan\theta_k < \sqrt{2}\sin\theta_k$ なので、 の 2 と合わせて $$ \eta_{k+1} \le \Delta\sin\theta_{k+1} \le \Delta\tan\theta_{k+1} \le \frac{2\Delta^2}{\delta}\tan^3\theta_k \le \frac{2\Delta^2}{\delta}\bigl(\sqrt{2}\sin\theta_k\bigr)^3 \le \frac{2\Delta^2}{\delta}\cdot 2\sqrt{2}\Bigl(\frac{2\eta_k}{\delta}\Bigr)^3 = \frac{32\sqrt{2}\,\Delta^2}{\delta^4}\eta_k^3 $$ です。 $A$ が対称な三重対角行列なら $\eta_k = |(A_k)_{n,n-1}|$ です。$\eta_k$ が単位丸め程度になったら最下行を切り離し（デフレーション）、$(A_k)_{nn}$ を固有値として、残りの $n - 1$ 次の行列で続けます。非対称な行列では Rayleigh 商反復は 2 次収束なので（）、上のコードのとおり $|h_{n,n-1}|$ も 2 次で収束します。 次の命題は、副対角成分を $0$ とみなして切り離すことで、どれだけ精度を失うかを測ります。結論は、切り離したあとに得られる固有値が、元の行列 $A$ に単位丸め程度の小さな変化を加えた行列の固有値になっている、ということです。 $H = Q^{\mathsf{T}}AQ$（$Q$ は直交）を上 Hessenberg 行列とし、ある $k$ で $$ |h_{k+1,k}| \le u\,(|h_{kk}| + |h_{k+1,k+1}|) $$ とする（$u$ は単位丸め）。$h_{k+1,k}$ を $0$ に置き換えた行列を $H + E$ とすると、 1. $\|E\|_2 = |h_{k+1,k}| \le 2u\|H\|_2 = 2u\|A\|_2$ で、$H + E$ の固有値は左上 $k$ 次と右下 $n - k$ 次の対角ブロックの固有値を合わせたものであり、それらは $A + QEQ^{\mathsf{T}}$（$\|QEQ^{\mathsf{T}}\|_2 = \|E\|_2$）の固有値に等しい。 2. $A$ が対称で、$H$ が対称な三重対角行列のとき、$h_{k+1,k}$ と $h_{k,k+1}$ をともに $0$ に置き換えた行列を $H + E'$ とすると $\|E'\|_2 = |h_{k+1,k}|$ で、大きい順に並べた固有値は $|\lambda_i(H + E') - \lambda_i(A)| \le 2u\|A\|_2$ を満たす。 3. $A = V\Lambda V^{-1}$ が対角化可能なら、$H + E$ の各固有値 $\nu$ について $\min_i|\nu - \lambda_i| \le 2u\,\kappa_2(V)\|A\|_2$ である。 1. $E = -h_{k+1,k}\boldsymbol{e}_{k+1}\boldsymbol{e}_k^{\mathsf{T}}$ で、$\|E\boldsymbol{x}\|_2 = |h_{k+1,k}||x_k| \le |h_{k+1,k}|\|\boldsymbol{x}\|_2$（$\boldsymbol{x} = \boldsymbol{e}_k$ で等号）です。$|h_{ii}| = |\boldsymbol{e}_i^{\mathsf{T}}H\boldsymbol{e}_i| \le \|H\|_2$（Cauchy–Schwarz の不等式 ）から $\|E\|_2 \le 2u\|H\|_2$ で、直交行列を掛けても 2 ノルムは変わらないので $\|H\|_2 = \|A\|_2$ です。$H + E$ はブロック上三角なので、特性多項式はブロック三角行列の行列式（）により 2 つの対角ブロックの特性多項式の積です。$H + E = Q^{\mathsf{T}}(A + QEQ^{\mathsf{T}})Q$ なので、固有値は $A + QEQ^{\mathsf{T}}$ のものと同じです。 2. $E'$ は $(k+1, k)$ と $(k, k+1)$ に $-h_{k+1,k}$ をもつ対称行列で、$\|E'\boldsymbol{x}\|_2^2 = h_{k+1,k}^2(x_k^2 + x_{k+1}^2) \le h_{k+1,k}^2\|\boldsymbol{x}\|_2^2$ から $\|E'\|_2 = |h_{k+1,k}|$ です。$H$ と $A$ は同じ固有値をもち、$H$ と $E'$ は対称なので、エルミート行列の固有値の摂動（）と 1 の評価から従います。 3. $H$ は $A$ と相似で $H = (Q^{\mathsf{T}}V)\Lambda(Q^{\mathsf{T}}V)^{-1}$、$\kappa_2(Q^{\mathsf{T}}V) = \kappa_2(V)$ です。Bauer–Fike の定理（）と 1 から従います。 対称行列では、デフレーションの誤差は固有値の絶対誤差として $2u\|A\|_2$ 以下で、どの固有値にも同じ精度が保証されます。非対称行列では $\kappa_2(V)$ 倍に拡大されうるので、固有値の精度は固有ベクトルの条件で決まります。 のように Rayleigh 商シフトは停滞しうるので、対称三重対角行列では、右下 $2\times2$ ブロックの 2 つの固有値のうち $(A_{k-1})_{nn}$ に近い方をシフトにとる **Wilkinson シフト**を使います。Wilkinson は、このシフトで QR 法がすべての対称三重対角行列について収束し（大域収束）、その速さが少なくとも 2 次、通常は 3 次であることを示しました。証明は数頁の場合分けで、本章では定理として扱いません。J. H. Wilkinson, "Global convergence of tridiagonal QR algorithm with origin shifts", *Linear Algebra Appl.* 1 (1968)、Parlett『The Symmetric Eigenvalue Problem』第 8 章を参照してください。 実の非対称行列は複素共役な固有値の対をもちえ、実のシフトではそれに収束しません。複素共役なシフト $\mu$, $\bar{\mu}$（右下 $2\times2$ ブロックの固有値）で 2 回続けて QR 反復を行うと、結果は実行列 $(A - \mu I)(A - \bar{\mu}I) = A^2 - 2\operatorname{Re}\mu\,A + |\mu|^2I$ の QR 分解だけで決まり、実数演算で実行できます。Francis は、この行列の第 1 列だけを計算し、 に基づいて Hessenberg 形を回復することで、1 回 $O(n^2)$ で実行しました。収束後の $A_k$ は $1\times1$ と $2\times2$ の対角ブロックをもつ準上三角行列（実 Schur 形）になります。J. G. F. Francis, "The QR transformation", *Comput. J.* 4 (1961/62)、Watkins の前掲書第 5・6 章、Golub–Van Loan §7.5 を参照してください。 現代の実装（LAPACK の dhseqr）は、最下行だけでなく右下の大きめのブロック全体について、その部分の Schur 分解から切り離せる固有値を探し、一度に多くの固有値をデフレーションします。また、多数のシフトを小さなバルジとして同時に追跡して、キャッシュを有効に使います。K. Braman, R. Byers, R. Mathias, "The multishift QR algorithm. Part II: Aggressive early deflation", *SIAM J. Matrix Anal. Appl.* 23 (2002) を参照してください。 ## 6. 演習 $A = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}$ に Rayleigh 商シフト $\mu_1 = a_{22} = 1$ の QR 反復を 1 回行い、$A_1$ を求めよ。$A$ の固有値は $(3 \pm \sqrt{5})/2$ である。$\boldsymbol{e}_2$ と小さい方の固有値の固有ベクトルのなす角を $\theta$、$A_1$ についての同じ角を $\theta'$ とするとき、$\tan\theta' = \tan^3\theta$ を確かめよ。 $A - I = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$ の第 1 列を正規化して $\boldsymbol{q}_1 = (1, 1)^{\mathsf{T}}/\sqrt{2}$、$r_{11} = \sqrt{2}$、$r_{12} = \boldsymbol{q}_1^{\mathsf{T}}(1, 0)^{\mathsf{T}} = 1/\sqrt{2}$。$r_{22} > 0$ となるよう $\boldsymbol{q}_2 = (1, -1)^{\mathsf{T}}/\sqrt{2}$ とすると $r_{22} = 1/\sqrt{2}$ です。 $$ RQ = \begin{pmatrix} \sqrt{2} & 1/\sqrt{2} \\ 0 & 1/\sqrt{2} \end{pmatrix}\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \begin{pmatrix} 3/2 & 1/2 \\ 1/2 & -1/2 \end{pmatrix},\qquad A_1 = \begin{pmatrix} 5/2 & 1/2 \\ 1/2 & 1/2 \end{pmatrix} . $$ 跡 $3$ と行列式 $1$ は保たれています。小さい固有値 $\lambda = (3 - \sqrt{5})/2$ の固有ベクトルは、$A$ では $(1, \lambda - 2)$ に比例するので $\tan\theta = 1/(2 - \lambda) = 2/(1 + \sqrt{5}) \approx 0.618$、$A_1$ では $(A_1 - \lambda I)\boldsymbol{v} = \boldsymbol{0}$ の第 1 行から $\boldsymbol{v} \propto (1/2, \lambda - 5/2)$ なので $\tan\theta' = (1/2)/(5/2 - \lambda) = 1/(2 + \sqrt{5}) \approx 0.236$ です。$\tan^3\theta = 8/(1 + \sqrt{5})^3 = 8/(16 + 8\sqrt{5}) = 1/(2 + \sqrt{5})$ なので、$\tan\theta' = \tan^3\theta$ がちょうど成り立ちます（2 次の Rayleigh 商反復の性質。冪乗法の章の演習を参照）。 で $H$ が既約という仮定を外すと結論が成り立たない例を挙げよ。 $A = \operatorname{diag}(1, 2, 3)$ とし、$Q = I$、$V$ を第 2 列と第 3 列を入れ替える置換行列とします。どちらも第 1 列は $\boldsymbol{e}_1$ で、$H = A$ と $G = V^{\mathsf{T}}AV = \operatorname{diag}(1, 3, 2)$ はともに（対角なので）上 Hessenberg 行列です。しかし $V \ne QD$（$D$ は $\pm1$ の対角）です。$H$ の副対角は $0$ なので既約ではありません。証明では $R_H$ が正則であることにこの仮定を使いました。 $n$ 次の行列のすべての固有値を、(a) 密行列のまま QR 法で、(b) Hessenberg 化してから QR 法で求めるとき、固有値 1 つあたり平均 3 回の反復で済むとして、全体の演算量の主要項を比べよ。QR 分解と $RQ$ の積は、密行列で約 $(4/3)n^3 + n^3$、Hessenberg 行列で約 $6n^2$ とする。 (a) 各反復は約 $(7/3)m^3$（$m$ はデフレーションで縮んだ次数）で、$m = n, n-1, \ldots, 1$ について 3 回ずつなので約 $3\cdot(7/3)\sum m^3 \approx (7/4)n^4$ です。(b) Hessenberg 化が約 $(10/3)n^3$、反復は $3\cdot 6\sum m^2 \approx 6n^3$ で、合わせて約 $(28/3)n^3 \approx 9.3n^3$ です（固有ベクトルも求めるなら $Q$ の累積の分が加わる）。$n = 1000$ なら (a) は約 $1.8\times10^{12}$、(b) は約 $9.3\times10^9$ で、約 190 倍の差になります。 ## 参考文献 - G. H. Golub, C. F. Van Loan『Matrix Computations』4th ed., Johns Hopkins University Press, 2013 — §7.3（同時反復と QR 法）、§7.4（Hessenberg 化と陰的 Q の定理）、§7.5（Francis の QR 法）、§8.3（対称 QR 法）。 - D. S. Watkins『Fundamentals of Matrix Computations』3rd ed., Wiley, 2010 — 第 5・6 章。 - B. N. Parlett『The Symmetric Eigenvalue Problem』SIAM, 1998 — 第 8 章（QL 法・QR 法とシフト）。 - L. N. Trefethen, D. Bau III『Numerical Linear Algebra』SIAM, 1997 — Lectures 26（Hessenberg 化）、28–29（QR 法とシフト）。 - J. G. F. Francis, "The QR transformation: a unitary analogue to the LR transformation", *The Computer Journal* 4 (1961), 265–271 および 332–345.