行列の固有値（行列を掛けても向きが変わらないベクトルに対する倍率）を、計算機でまとめて求めるにはどうすればよいか、という話です。ここで扱う QR 法は、「行列を 2 つの行列の積に分け、順序を入れ替えて掛け直す」という操作を繰り返すだけの方法で、続けるうちに対角に固有値が並びます。なぜこの単純な操作で固有値が出るのか、どうすれば速く正確にできるのかを、この記事で順に確かめます。読み終えると、numpy などの固有値計算の中身で何が起きているかが説明できるようになります。
QR 法 は、A k − 1 = Q k R k A_{k-1} = Q_kR_k A k − 1 = Q k R k A k = R k Q k A_k = R_kQ_k A k = R k Q k A k = Q k T A k − 1 Q k A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k A k = Q k T A k − 1 Q k A k A_k A k 0 0 0 掛け合わせた直交行列 Q k = Q 1 ⋯ Q k \mathbf{Q}_k = Q_1\cdots Q_k Q k = Q 1 ⋯ Q k A k A^k A k Q Q Q 同時反復 **そのものです。このため冪乗法の章の同時反復の収束がそのまま QR 法の収束になります。
∣ λ 1 ∣ > ⋯ > ∣ λ n ∣ > 0 |\lambda_1| > \cdots > |\lambda_n| > 0 ∣ λ 1 ∣ > ⋯ > ∣ λ n ∣ > 0 A k A_k A k O ( r k ) O(r^k) O ( r k ) r = max j ∣ λ j + 1 / λ j ∣ r = \max_j|\lambda_{j+1}/\lambda_j| r = max j ∣ λ j + 1 / λ j ∣ 0 0 0 実用の QR 法は、最初に Householder 変換で Hessenberg 形 （対角のすぐ下の斜めの列までにしか成分がない形）に直し（10 n 3 / 3 10n^3/3 10 n 3 /3 O ( n 2 ) O(n^2) O ( n 2 ) 陰的 Q の定理 により、Hessenberg 形への直交相似は第 1 列だけで決まります。
シフト μ k = ( A k − 1 ) n n \mu_k = (A_{k-1})_{nn} μ k = ( A k − 1 ) nn Q k \mathbf{Q}_k Q k 0 0 0 ∣ h k + 1 , k ∣ |h_{k+1,k}| ∣ h k + 1 , k ∣ 0 0 0 ∥ A ∥ \|A\| ∥ A ∥
同時反復 : 複数のベクトルに行列を掛けては直交化し直すことで、大きい固有値に対応する不変部分空間を求める反復。→ 同時反復 Rayleigh 商反復 : Rayleigh 商をシフトに使った逆反復を繰り返し、固有ベクトルに速く近づく方法。→ Rayleigh 商反復 Householder 鏡映 : ベクトルを、ある決めた向きへ写す直交行列（鏡に映す変換）。→ Householder 鏡映 Givens 回転 : 2 つの成分だけを回転させて片方を 0 0 0 Givens 回転 LU 分解 : 行列を単位下三角行列と上三角行列の積に分けること。→ LU 分解 相似 : ある正則行列 P P P P − 1 A P P^{-1}AP P − 1 A P 相似 冪乗法・逆反復の章 の方法は、固有値を 1 つずつ、あるいは上位の数個を求めるものでした。密行列のすべての固有値を求める標準的な算法は QR 法で、numpy.linalg.eig・scipy.linalg.eigvals（LAPACK の geev）も、対称行列の eigh（syev 系）も、Hessenberg 化（対称なら三重対角化）のあとに QR 法か、同じ三重対角形の上に立つ後継の算法（分割統治法など）を使っています。
QR 法は 1961 年に Francis と Kublanovskaya が独立に発表しました。算法の記述は 1 行ですが、なぜ収束するのかは一見わかりません。本章の主題は、QR 法が「すべての固有ベクトルについて同時に冪乗法を行い、互いに直交させ続ける」算法であることを示し、そこから収束・計算の工夫（Hessenberg 化）・加速（シフト）を順に導くことです。
この章では行列はすべて実行列とし、直交行列・転置で書きます。QR 分解の章 の一意性（定理 4.1［QR 分解］ の 2。列が一次独立なら、R R R ∥ ⋅ ∥ 2 \|\cdot\|_2 ∥ ⋅ ∥ 2 Frobenius ノルム ∥ ⋅ ∥ F \|\cdot\|_F ∥ ⋅ ∥ F ∥ B C ∥ F ≤ ∥ B ∥ 2 ∥ C ∥ F \|BC\|_F \le \|B\|_2\|C\|_F ∥ B C ∥ F ≤ ∥ B ∥ 2 ∥ C ∥ F C C C B B B ∥ C B ∥ F ≤ ∥ C ∥ F ∥ B ∥ 2 \|CB\|_F \le \|C\|_F\|B\|_2 ∥ C B ∥ F ≤ ∥ C ∥ F ∥ B ∥ 2 ∥ B T ∥ 2 = ∥ B ∥ 2 \|B^{\mathsf{T}}\|_2 = \|B\|_2 ∥ B T ∥ 2 = ∥ B ∥ 2 命題 4.2［行列のノルムと線形方程式の条件数］ の 3）と、直交行列 Q Q Q ∥ Q C ∥ F = ∥ C ∥ F \|QC\|_F = \|C\|_F ∥ QC ∥ F = ∥ C ∥ F ∥ Q C ∥ 2 = ∥ C ∥ 2 \|QC\|_2 = \|C\|_2 ∥ QC ∥ 2 = ∥ C ∥ 2
定義 2.1 （QR 法（シフトなし） ） ↗
A ∈ R n × n A \in \mathbb{R}^{n\times n} A ∈ R n × n A 0 = A A_0 = A A 0 = A k = 1 , 2 , … k = 1, 2, \ldots k = 1 , 2 , … A k − 1 A_{k-1} A k − 1 A k − 1 = Q k R k A_{k-1} = Q_kR_k A k − 1 = Q k R k Q k Q_k Q k R k R_k R k
A k = R k Q k A_k = R_kQ_k A k = R k Q k とする。Q k = Q 1 Q 2 ⋯ Q k \mathbf{Q}_k = Q_1Q_2\cdots Q_k Q k = Q 1 Q 2 ⋯ Q k R k = R k R k − 1 ⋯ R 1 \mathbf{R}_k = R_kR_{k-1}\cdots R_1 R k = R k R k − 1 ⋯ R 1 Q 0 = R 0 = I \mathbf{Q}_0 = \mathbf{R}_0 = I Q 0 = R 0 = I
命題 2.2 （QR 反復は直交相似変換である ） ✓
A k = Q k T A k − 1 Q k = Q k T A Q k A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k A k = Q k T A k − 1 Q k = Q k T A Q k A k A_k A k A A A
証明（命題 2.2） R k = Q k T A k − 1 R_k = Q_k^{\mathsf{T}}A_{k-1} R k = Q k T A k − 1 A k = R k Q k = Q k T A k − 1 Q k A_k = R_kQ_k = Q_k^{\mathsf{T}}A_{k-1}Q_k A k = R k Q k = Q k T A k − 1 Q k A k = Q k T A Q k A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k A k = Q k T A Q k Q k \mathbf{Q}_k Q k Q k T = Q k − 1 \mathbf{Q}_k^{\mathsf{T}} = \mathbf{Q}_k^{-1} Q k T = Q k − 1 定理 6.1［行列式とその性質］ ）ので固有値も等しく、det A k = det A ≠ 0 \det A_k = \det A \ne 0 det A k = det A = 0 A k A_k A k 定理 4.1［QR 分解］ の 2）。
∎ 定理 2.3 （QR 法と同時反復の同値性 ） ✓
定義 2.1 の記号で、すべての k k k
A k = Q k R k A^k = \mathbf{Q}_k\mathbf{R}_k A k = Q k R k A k A^k A k R R R Z 0 = I Z_0 = I Z 0 = I 定義 6.1［固有値の計算 I］ ）の Z k Z_k Z k Q k \mathbf{Q}_k Q k 証明（定理 2.3）
k k k k = 0 k = 0 k = 0 A k − 1 = Q k − 1 R k − 1 A^{k-1} = \mathbf{Q}_{k-1}\mathbf{R}_{k-1} A k − 1 = Q k − 1 R k − 1 命題 2.2 から A Q k − 1 = Q k − 1 A k − 1 A\mathbf{Q}_{k-1} = \mathbf{Q}_{k-1}A_{k-1} A Q k − 1 = Q k − 1 A k − 1 A k = A Q k − 1 R k − 1 = Q k − 1 A k − 1 R k − 1 = Q k − 1 Q k R k R k − 1 = Q k R k 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 A k = A Q k − 1 R k − 1 = Q k − 1 A k − 1 R k − 1 = Q k − 1 Q k R k R k − 1 = Q k R k です。R k \mathbf{R}_k R k
帰納法で Z k − 1 = Q k − 1 Z_{k-1} = \mathbf{Q}_{k-1} Z k − 1 = Q k − 1 A Z k − 1 AZ_{k-1} A Z k − 1 A Z k − 1 = Z k R k ′ AZ_{k-1} = Z_kR'_k A Z k − 1 = Z k R k ′ R k ′ R'_k R k ′ A Q k − 1 = Q k − 1 A k − 1 = Q k − 1 Q k R k = Q 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 A Q k − 1 = Q k − 1 A k − 1 = Q k − 1 Q k R k = Q k R k R R R A Q k − 1 A\mathbf{Q}_{k-1} A Q k − 1 定理 4.1［QR 分解］ の 2）から Z k = Q k Z_k = \mathbf{Q}_k Z k = Q k
∎ QR 法は、n n n Q k − 1 \mathbf{Q}_{k-1} Q k − 1 A A A 同時反復 を、「A A A Q k \mathbf{Q}_k Q k A k A_k A k j j j j j j 系 6.3［固有値の計算 I］ ）。Q k \mathbf{Q}_k Q k j j j A A A A k = Q k T A Q k A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k A k = Q k T A Q k ( n − j ) × j (n-j)\times j ( n − j ) × j 0 0 0 A k A_k A k
収束の速さを定量的に示すには、「単位行列に近い行列の QR 分解の Q Q Q I + F I + F I + F F F F Q Q Q I I I F F F n \sqrt{n} n 1 1 1 Q Q Q I I I
補題 3.1 （単位行列に近い行列の QR 分解 ） ✓
F ∈ R n × n F \in \mathbb{R}^{n\times n} F ∈ R n × n s = ∥ F ∥ F s = \|F\|_F s = ∥ F ∥ F ( n + 1 ) s 2 < 1 (n+1)s^2 < 1 ( n + 1 ) s 2 < 1 I + F I + F I + F I + F = Q ~ R ~ I + F = \tilde{Q}\tilde{R} I + F = Q ~ R ~ R ~ \tilde{R} R ~
∥ Q ~ − I ∥ F ≤ n s + n s 2 \|\tilde{Q} - I\|_F \le \sqrt{n}\,s + n s^2 ∥ Q ~ − I ∥ F ≤ n s + n s 2 が成り立つ。
証明（補題 3.1） s < 1 s < 1 s < 1 ∥ ( I + F ) x ∥ 2 ≥ ∥ x ∥ 2 − ∥ F ∥ 2 ∥ x ∥ 2 ≥ ( 1 − s ) ∥ x ∥ 2 \|(I + F)\boldsymbol{x}\|_2 \ge \|\boldsymbol{x}\|_2 - \|F\|_2\|\boldsymbol{x}\|_2 \ge (1 - s)\|\boldsymbol{x}\|_2 ∥ ( I + F ) x ∥ 2 ≥ ∥ x ∥ 2 − ∥ F ∥ 2 ∥ x ∥ 2 ≥ ( 1 − s ) ∥ x ∥ 2 ∥ F ∥ 2 ≤ ∥ F ∥ F \|F\|_2 \le \|F\|_F ∥ F ∥ 2 ≤ ∥ F ∥ F 定理 4.2［内積空間とグラム・シュミット直交化］ から）で、I + F I + F I + F U = Q ~ T = ( u i j ) U = \tilde{Q}^{\mathsf{T}} = (u_{ij}) U = Q ~ T = ( u ij )
U U U Q ~ T ( I + F ) = R ~ \tilde{Q}^{\mathsf{T}}(I + F) = \tilde{R} Q ~ T ( I + F ) = R ~ U = R ~ − U F U = \tilde{R} - UF U = R ~ − U F R ~ \tilde{R} R ~ U U U − U F -UF − U F ∥ U F ∥ F = s \|UF\|_F = s ∥ U F ∥ F = s
狭義上三角部分も小さい. U U U 1 1 1 1 ≤ m < n 1 \le m < n 1 ≤ m < n m m m m m m m m m ∑ i , j ≤ m u i j 2 \sum_{i,j\le m}u_{ij}^2 ∑ i , j ≤ m u ij 2
∑ i ≤ m < j u i j 2 = ∑ i ≤ m < j u j i 2 \sum_{i\le m<j}u_{ij}^2 = \sum_{i\le m<j}u_{ji}^2 i ≤ m < j ∑ u ij 2 = i ≤ m < j ∑ u j i 2 です。これを m = 1 , … , n − 1 m = 1, \ldots, n-1 m = 1 , … , n − 1 ( i , j ) (i, j) ( i , j ) i < j i < j i < j m = i , … , j − 1 m = i, \ldots, j - 1 m = i , … , j − 1 j − i j - i j − i
∑ i < j u i j 2 ≤ ∑ i < j ( j − i ) u i j 2 = ∑ i < j ( j − i ) u j i 2 ≤ ( n − 1 ) ∑ i > j u i j 2 ≤ ( n − 1 ) s 2 \sum_{i<j}u_{ij}^2 \le \sum_{i<j}(j - i)u_{ij}^2 = \sum_{i<j}(j - i)u_{ji}^2 \le (n - 1)\sum_{i>j}u_{ij}^2 \le (n - 1)s^2 i < j ∑ u ij 2 ≤ i < j ∑ ( j − i ) u ij 2 = i < j ∑ ( j − i ) u j i 2 ≤ ( n − 1 ) i > j ∑ u ij 2 ≤ ( n − 1 ) s 2 です。よって対角を除いた部分を U − D U - D U − D D D D U U U ∥ U − D ∥ F 2 ≤ n s 2 \|U - D\|_F^2 \le ns^2 ∥ U − D ∥ F 2 ≤ n s 2
対角は 1 1 1 第 j j j 1 1 1 u j j 2 = 1 − ∑ i ≠ j u i j 2 ≥ 1 − n s 2 > s 2 u_{jj}^2 = 1 - \sum_{i\ne j}u_{ij}^2 \ge 1 - ns^2 > s^2 u j j 2 = 1 − ∑ i = j u ij 2 ≥ 1 − n s 2 > s 2 U = R ~ − U F U = \tilde{R} - UF U = R ~ − U F u j j = r ~ j j − ( U F ) j j > − ∣ ( U F ) j j ∣ ≥ − s u_{jj} = \tilde{r}_{jj} - (UF)_{jj} > -|(UF)_{jj}| \ge -s u j j = r ~ j j − ( U F ) j j > − ∣ ( U F ) j j ∣ ≥ − s ∣ u j j ∣ > s |u_{jj}| > s ∣ u j j ∣ > s u j j > − s u_{jj} > -s u j j > − s u j j > 0 u_{jj} > 0 u j j > 0 0 < u j j ≤ 1 0 < u_{jj} \le 1 0 < u j j ≤ 1 1 − u j j ≤ 1 − u j j 2 = ∑ i ≠ j u i j 2 1 - u_{jj} \le 1 - u_{jj}^2 = \sum_{i\ne j}u_{ij}^2 1 − u j j ≤ 1 − u j j 2 = ∑ i = j u ij 2 j j j ∥ D − I ∥ F ≤ ∑ j ( 1 − u j j ) ≤ ∥ U − D ∥ F 2 ≤ n s 2 \|D - I\|_F \le \sum_j(1 - u_{jj}) \le \|U - D\|_F^2 \le ns^2 ∥ D − I ∥ F ≤ ∑ j ( 1 − u j j ) ≤ ∥ U − D ∥ F 2 ≤ n s 2
以上から ∥ Q ~ − I ∥ F = ∥ U − I ∥ F ≤ ∥ U − D ∥ F + ∥ D − I ∥ F ≤ n s + n s 2 \|\tilde{Q} - I\|_F = \|U - I\|_F \le \|U - D\|_F + \|D - I\|_F \le \sqrt{n}\,s + ns^2 ∥ Q ~ − I ∥ F = ∥ U − I ∥ F ≤ ∥ U − D ∥ F + ∥ D − I ∥ F ≤ n s + n s 2
∎ 次の定理が、この節の本題です。Q k \mathbf{Q}_k Q k Q X Q_X Q X A k A_k A k 0 0 0 X X X S S S A A A D k D_k D k ± 1 \pm1 ± 1 A k = D k Q ~ k T S Q ~ k D k A_k = D_k\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_kD_k A k = D k Q ~ k T S Q ~ k D k Q ~ k ≈ I \tilde{Q}_k \approx I Q ~ k ≈ I A k ≈ D k S D k A_k \approx D_kSD_k A k ≈ D k S D k S S S
定理 3.2 （シフトなし QR 法の収束 ） ✓
A ∈ R n × n A \in \mathbb{R}^{n\times n} A ∈ R n × n A = X Λ X − 1 A = X\Lambda X^{-1} A = X Λ X − 1 Λ = diag ( λ 1 , … , λ n ) \Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_n) Λ = diag ( λ 1 , … , λ n )
∣ λ 1 ∣ > ∣ λ 2 ∣ > ⋯ > ∣ λ n ∣ > 0 |\lambda_1| > |\lambda_2| > \cdots > |\lambda_n| > 0 ∣ λ 1 ∣ > ∣ λ 2 ∣ > ⋯ > ∣ λ n ∣ > 0 と対角化でき、X − 1 X^{-1} X − 1 X − 1 = L U X^{-1} = LU X − 1 = LU L L L U U U 0 0 0 r = max 1 ≤ j < n ∣ λ j + 1 / λ j ∣ r = \max_{1\le j<n}|\lambda_{j+1}/\lambda_j| r = max 1 ≤ j < n ∣ λ j + 1 / λ j ∣ C C C A k A_k A k k k k
∣ ( A k ) i j ∣ ≤ C r k ( i > j ) , ∣ ( A k ) i i − λ i ∣ ≤ C r k |(A_k)_{ij}| \le Cr^k\quad(i > j),\qquad |(A_k)_{ii} - \lambda_i| \le Cr^k ∣ ( A k ) ij ∣ ≤ C r k ( i > j ) , ∣ ( A k ) ii − λ i ∣ ≤ C r k を満たす。とくに A k A_k A k 0 0 0
証明（定理 3.2） 準備. 実行列の実でない固有値は共役な対で現れ絶対値が等しいので、仮定から固有値はすべて実数です。各固有空間は 1 次元なので、X X X X ′ = X D ′ X' = XD' X ′ = X D ′ D ′ D' D ′ X ′ − 1 = D ′ − 1 L U = ( D ′ − 1 L D ′ ) ( D ′ − 1 U ) X'^{-1} = D'^{-1}LU = (D'^{-1}LD')(D'^{-1}U) X ′ − 1 = D ′ − 1 LU = ( D ′ − 1 L D ′ ) ( D ′ − 1 U ) X X X X − 1 X^{-1} X − 1 L L L U U U 0 0 0 定理 2.2［Gauss 消去と LU 分解］ により実数の範囲でも LU 分解をもち、その一意性から両者は一致する）。X X X X = Q X R X X = Q_XR_X X = Q X R X 定理 4.1［QR 分解］ の 2）とし、S = R X Λ R X − 1 S = R_X\Lambda R_X^{-1} S = R X Λ R X − 1 S S S λ 1 , … , λ n \lambda_1, \ldots, \lambda_n λ 1 , … , λ n Q X T A Q X = R X Λ R X − 1 = S Q_X^{\mathsf{T}}AQ_X = R_X\Lambda R_X^{-1} = S Q X T A Q X = R X Λ R X − 1 = S A A A
A k A^k A k A k = X Λ k X − 1 = X Λ k L U A^k = X\Lambda^kX^{-1} = X\Lambda^kLU A k = X Λ k X − 1 = X Λ k LU Λ k L = ( Λ k L Λ − k ) Λ k \Lambda^kL = (\Lambda^kL\Lambda^{-k})\Lambda^k Λ k L = ( Λ k L Λ − k ) Λ k E k = Λ k L Λ − k − I E_k = \Lambda^kL\Lambda^{-k} - I E k = Λ k L Λ − k − I ( i , j ) (i, j) ( i , j ) l i j ( λ i / λ j ) k l_{ij}(\lambda_i/\lambda_j)^k l ij ( λ i / λ j ) k i > j i > j i > j ∣ λ i / λ j ∣ ≤ r |\lambda_i/\lambda_j| \le r ∣ λ i / λ j ∣ ≤ r ∥ E k ∥ F ≤ ∥ L ∥ F r k \|E_k\|_F \le \|L\|_Fr^k ∥ E k ∥ F ≤ ∥ L ∥ F r k
A k = Q X R X ( I + E k ) Λ k U = Q X ( I + F k ) R X Λ k U , F k = R X E k R X − 1 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} A k = Q X R X ( I + E k ) Λ k U = Q X ( I + F k ) R X Λ k U , F k = R X E k R X − 1 で、∥ F k ∥ F ≤ ∥ R X ∥ 2 ∥ R X − 1 ∥ 2 ∥ L ∥ F r k = : c 1 r k \|F_k\|_F \le \|R_X\|_2\|R_X^{-1}\|_2\|L\|_Fr^k =: c_1r^k ∥ F k ∥ F ≤ ∥ R X ∥ 2 ∥ R X − 1 ∥ 2 ∥ L ∥ F r k =: c 1 r k
Q k \mathbf{Q}_k Q k ( n + 1 ) c 1 2 r 2 k < 1 (n+1)c_1^2r^{2k} < 1 ( n + 1 ) c 1 2 r 2 k < 1 k k k r < 1 r < 1 r < 1 k k k 補題 3.1 から I + F k = Q ~ k R ~ k I + F_k = \tilde{Q}_k\tilde{R}_k I + F k = Q ~ k R ~ k ∥ Q ~ k − I ∥ F ≤ n c 1 r k + n c 1 2 r 2 k \|\tilde{Q}_k - I\|_F \le \sqrt{n}c_1r^k + nc_1^2r^{2k} ∥ Q ~ k − I ∥ F ≤ n c 1 r k + n c 1 2 r 2 k R ~ k R X Λ k U \tilde{R}_kR_X\Lambda^kU R ~ k R X Λ k U ( R ~ k ) i i ( R X ) i i λ i k u i i (\tilde{R}_k)_{ii}(R_X)_{ii}\lambda_i^ku_{ii} ( R ~ k ) ii ( R X ) ii λ i k u ii 0 0 0 D k D_k D k ± 1 \pm1 ± 1
A k = ( Q X Q ~ k D k ) ( D k R ~ k R X Λ k U ) A^k = (Q_X\tilde{Q}_kD_k)(D_k\tilde{R}_kR_X\Lambda^kU) A k = ( Q X Q ~ k D k ) ( D k R ~ k R X Λ k U ) は、右の因子の対角が正の QR 分解です。定理 2.3 の 1 と QR 分解の一意性から Q k = Q X Q ~ k D k \mathbf{Q}_k = Q_X\tilde{Q}_kD_k Q k = Q X Q ~ k D k
A k A_k A k 命題 2.2 から
A k = Q k T A Q k = D k Q ~ k T S Q ~ k D k A_k = \mathbf{Q}_k^{\mathsf{T}}A\mathbf{Q}_k = D_k\tilde{Q}_k^{\mathsf{T}}S\tilde{Q}_kD_k A k = Q k T A Q k = D k Q ~ k T S Q ~ k D k です。Q ~ k T S Q ~ k − S = ( Q ~ k − I ) T S Q ~ k + S ( Q ~ k − I ) \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) Q ~ k T S Q ~ k − S = ( Q ~ k − I ) T S Q ~ k + S ( Q ~ k − I )
∥ A k − D k S D k ∥ F = ∥ Q ~ k T S Q ~ k − S ∥ F ≤ 2 ∥ S ∥ 2 ∥ Q ~ k − I ∥ F ≤ 2 ∥ S ∥ 2 ( n c 1 + n c 1 2 ) r k \|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 ∥ A k − D k S D k ∥ F = ∥ Q ~ k T S Q ~ k − S ∥ F ≤ 2∥ S ∥ 2 ∥ Q ~ k − I ∥ F ≤ 2∥ S ∥ 2 ( n c 1 + n c 1 2 ) r k です。D k S D k D_kSD_k D k S D k λ i \lambda_i λ i k k k k k k C C C
∎ A k A_k A k D k S D k D_kSD_k D k S D k S S S ± 1 \pm1 ± 1 D k D_k D k λ i k \lambda_i^k λ i k 0 0 0
# Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）／乱数種 20261006
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 分解
return Q * s, R * s[:, None ]
sub[k] = np. abs ( np. diag ( Ak , - 1 ))
print ( " i |a(i+1,i)|(k=60) 20 回の比 (l(i+1)/l(i))^20 " )
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 ))
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 , … , 1 6, 5, \ldots, 1 6 , 5 , … , 1 ( i + 1 , i ) (i+1, i) ( i + 1 , i ) ( λ i + 1 / λ i ) 20 (\lambda_{i+1}/\lambda_i)^{20} ( λ i + 1 / λ i ) 20 注意 3.3 の成分ごとの速さが見えます。定理 3.2 の r = 5 / 6 r = 5/6 r = 5/6 ( 2 , 1 ) (2, 1) ( 2 , 1 ) X X X X − 1 X^{-1} X − 1
例 3.4 （固有値の絶対値が等しいと収束しない、シフトで直る ）
A = ( 0 1 1 0 ) A = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} A = ( 0 1 1 0 ) ± 1 \pm1 ± 1 定理 3.2 の仮定を満たしません。実際 A A A A = A ⋅ I A = A\cdot I A = A ⋅ I A 1 = I A = A A_1 = IA = A A 1 = I A = A
シフト μ \mu μ ∣ μ ∣ < 1 |\mu| < 1 ∣ μ ∣ < 1 A − μ I = Q R A - \mu I = QR A − μ I = QR A 1 = R Q + μ I A_1 = RQ + \mu I A 1 = R Q + μ I c = 1 + μ 2 c = 1 + \mu^2 c = 1 + μ 2
Q = 1 c ( − μ 1 1 μ ) , R = 1 c ( c − 2 μ 0 1 − μ 2 ) , A 1 = 1 c ( − 2 μ 1 − μ 2 1 − μ 2 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 = c 1 ( − μ 1 1 μ ) , R = c 1 ( c 0 − 2 μ 1 − μ 2 ) , A 1 = c 1 ( − 2 μ 1 − μ 2 1 − μ 2 2 μ ) です（Q Q Q Q R = A − μ I QR = A - \mu I QR = A − μ I R R R 1 1 1 ( 1 − μ 2 ) / ( 1 + μ 2 ) (1 - \mu^2)/(1 + \mu^2) ( 1 − μ 2 ) / ( 1 + μ 2 ) μ = 1 / 2 \mu = 1/2 μ = 1/2 0.6 0.6 0.6 μ = a 22 = 0 \mu = a_{22} = 0 μ = a 22 = 0 注意 5.5 ）です。
密行列の QR 分解は O ( n 3 ) O(n^3) O ( n 3 ) n n n O ( n 4 ) O(n^4) O ( n 4 )
定理 4.1 （Householder 変換による Hessenberg 化 ） ✓
任意の A ∈ R n × n A \in \mathbb{R}^{n\times n} A ∈ R n × n n ≥ 3 n \ge 3 n ≥ 3 n − 2 n - 2 n − 2 Q Q Q H = Q T A Q H = Q^{\mathsf{T}}AQ H = Q T A Q 上 Hessenberg 行列 （i > j + 1 i > j + 1 i > j + 1 h i j = 0 h_{ij} = 0 h ij = 0 10 n 3 / 3 10n^3/3 10 n 3 /3 A A A H H H
証明（定理 4.1） k = 1 , … , n − 2 k = 1, \ldots, n-2 k = 1 , … , n − 2 k k k k + 2 k+2 k + 2 0 0 0 k − 1 k - 1 k − 1 A ( k ) A^{(k)} A ( k ) 1 , … , k − 1 1, \ldots, k-1 1 , … , k − 1 x ∈ R n − k \boldsymbol{x} \in \mathbb{R}^{n-k} x ∈ R n − k k k k k + 1 k+1 k + 1 x ≠ 0 \boldsymbol{x} \ne \boldsymbol{0} x = 0 x \boldsymbol{x} x P ~ k \tilde{P}_k P ~ k 定義 3.1［QR 分解］ ）をとります。P ~ k x = α e 1 \tilde{P}_k\boldsymbol{x} = \alpha\boldsymbol{e}_1 P ~ k x = α e 1 P ~ k \tilde{P}_k P ~ k 命題 3.2［QR 分解］ ）。x = 0 \boldsymbol{x} = \boldsymbol{0} x = 0 P ~ k = I \tilde{P}_k = I P ~ k = I P k = diag ( I k , P ~ k ) P_k = \operatorname{diag}(I_k, \tilde{P}_k) P k = diag ( I k , P ~ k ) A ( k + 1 ) = P k A ( k ) P k A^{(k+1)} = P_kA^{(k)}P_k A ( k + 1 ) = P k A ( k ) P k
左から P k P_k P k k + 1 k+1 k + 1 1 , … , k 1, \ldots, k 1 , … , k k + 1 k + 1 k + 1 1 , … , k − 1 1, \ldots, k-1 1 , … , k − 1 0 0 0 0 0 0 k k k ( α , 0 , … , 0 ) T (\alpha, 0, \ldots, 0)^{\mathsf{T}} ( α , 0 , … , 0 ) T P k P_k P k k + 1 k+1 k + 1 1 , … , k 1, \ldots, k 1 , … , k A ( k + 1 ) A^{(k+1)} A ( k + 1 ) 1 , … , k 1, \ldots, k 1 , … , k Q = P 1 ⋯ P n − 2 Q = P_1\cdots P_{n-2} Q = P 1 ⋯ P n − 2 H = Q T A Q H = Q^{\mathsf{T}}AQ H = Q T A Q P k T = P k P_k^{\mathsf{T}} = P_k P k T = P k
演算量: m = n − k m = n - k m = n − k P ~ k \tilde{P}_k P ~ k m × ( m + 1 ) m\times(m+1) m × ( m + 1 ) 4 m 2 4m^2 4 m 2 n × m n\times m n × m 4 n m 4nm 4 nm P ~ y = y − β v ( v T y ) \tilde{P}\boldsymbol{y} = \boldsymbol{y} - \beta\boldsymbol{v}(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{y}) P ~ y = y − β v ( v T y ) 4 m 4m 4 m ∑ m ( 4 m 2 + 4 n m ) ≈ 4 n 3 / 3 + 2 n 3 = 10 n 3 / 3 \sum_m(4m^2 + 4nm) \approx 4n^3/3 + 2n^3 = 10n^3/3 ∑ m ( 4 m 2 + 4 nm ) ≈ 4 n 3 /3 + 2 n 3 = 10 n 3 /3 A A A H T = Q T A T Q = H H^{\mathsf{T}} = Q^{\mathsf{T}}A^{\mathsf{T}}Q = H H T = Q T A T Q = H
∎ 命題 4.2 （QR 反復は Hessenberg 形を保つ ） ✓
H H H H = Q R H = QR H = QR R R R Q Q Q R Q RQ R Q n − 1 n - 1 n − 1 O ( n 2 ) O(n^2) O ( n 2 ) R Q RQ R Q O ( n 2 ) O(n^2) O ( n 2 ) H H H R Q RQ R Q O ( n ) O(n) O ( n )
証明（命題 4.2） Q = H R − 1 Q = HR^{-1} Q = H R − 1 R − 1 R^{-1} R − 1 H H H j j j j + 1 j + 1 j + 1 H R − 1 HR^{-1} H R − 1 j j j H H H 1 , … , j 1, \ldots, j 1 , … , j j + 1 j + 1 j + 1 Q Q Q R Q RQ R Q j j j R R R 1 , … , j + 1 1, \ldots, j+1 1 , … , j + 1 R R R j + 1 j + 1 j + 1 R Q RQ R Q
計算は、i = 1 , … , n − 1 i = 1, \ldots, n-1 i = 1 , … , n − 1 i i i i + 1 i+1 i + 1 G i G_i G i ( i + 1 , i ) (i+1, i) ( i + 1 , i ) 0 0 0 命題 3.4［QR 分解］ 。回転は 2 行だけを変え、それぞれ O ( n ) O(n) O ( n ) G n − 1 ⋯ G 1 H = R ′ G_{n-1}\cdots G_1H = R' G n − 1 ⋯ G 1 H = R ′ R R R Q = G 1 T ⋯ G n − 1 T Q = G_1^{\mathsf{T}}\cdots G_{n-1}^{\mathsf{T}} Q = G 1 T ⋯ G n − 1 T R Q = R ′ G 1 T ⋯ G n − 1 T RQ = R'G_1^{\mathsf{T}}\cdots G_{n-1}^{\mathsf{T}} R Q = R ′ G 1 T ⋯ G n − 1 T O ( n 2 ) O(n^2) O ( n 2 ) R Q = Q T H Q RQ = Q^{\mathsf{T}}HQ R Q = Q T H Q O ( 1 ) O(1) O ( 1 ) O ( n ) O(n) O ( n )
∎ Hessenberg 化は一度だけ O ( n 3 ) O(n^3) O ( n 3 ) O ( n 2 ) O(n^2) O ( n 2 ) O ( n ) O(n) O ( n ) O ( n 3 ) O(n^3) O ( n 3 )
Hessenberg 形に直す直交相似変換は何通りもありそうですが、実は第 1 列を決めると、符号を除いてただ 1 通りしか残りません。次の定理がそれを保証します。「既約」とは、副対角成分 h k + 1 , k h_{k+1,k} h k + 1 , k 0 0 0
上 Hessenberg 形 三重対角形
6 次の行列で、成分が入りうる位置（塗り）。左が上 Hessenberg 形（対角のすぐ下の列まで）、右が三重対角形（対角とそのすぐ上下だけ）。塗っていない位置は 0。 定理 4.3 （陰的 Q の定理 ） ✓
A ∈ R n × n A \in \mathbb{R}^{n\times n} A ∈ R n × n Q Q Q V V V H = Q T A Q H = Q^{\mathsf{T}}AQ H = Q T A Q G = V T A V G = V^{\mathsf{T}}AV G = V T A V H H H 既約 （すべての k k k h k + 1 , k ≠ 0 h_{k+1,k} \ne 0 h k + 1 , k = 0 Q Q Q V V V ± 1 \pm1 ± 1 D D D V = Q D V = QD V = Q D G = D H D G = DHD G = D H D
証明（定理 4.3） q 1 = Q e 1 = V e 1 \boldsymbol{q}_1 = Q\boldsymbol{e}_1 = V\boldsymbol{e}_1 q 1 = Q e 1 = V e 1 K = ( q 1 A q 1 ⋯ A n − 1 q 1 ) K = (\boldsymbol{q}_1\ \ A\boldsymbol{q}_1\ \ \cdots\ \ A^{n-1}\boldsymbol{q}_1) K = ( q 1 A q 1 ⋯ A n − 1 q 1 ) A Q = Q H AQ = QH A Q = Q H A j Q = Q H j A^jQ = QH^j A j Q = Q H j
K = Q ( e 1 H e 1 ⋯ H n − 1 e 1 ) = : Q R H , K = V ( e 1 G e 1 ⋯ G n − 1 e 1 ) = : V R G 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 K = Q ( e 1 H e 1 ⋯ H n − 1 e 1 ) =: Q R H , K = V ( e 1 G e 1 ⋯ G n − 1 e 1 ) =: V R G です。Hessenberg 行列を掛けると成分の入る行が 1 つ下に延びるだけなので、帰納法で H j e 1 H^j\boldsymbol{e}_1 H j e 1 j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 h 21 h 32 ⋯ h j + 1 , j h_{21}h_{32}\cdots h_{j+1,j} h 21 h 32 ⋯ h j + 1 , j R H R_H R H 0 0 0 R H R_H R H R G R_G R G K = Q R H K = QR_H K = Q R H R G = V T K R_G = V^{\mathsf{T}}K R G = V T K
Q R H = V R G QR_H = VR_G Q R H = V R G V T Q = R G R H − 1 = : W V^{\mathsf{T}}Q = R_GR_H^{-1} =: W V T Q = R G R H − 1 =: W W − 1 = W T W^{-1} = W^{\mathsf{T}} W − 1 = W T ± 1 \pm1 ± 1 W = D W = D W = D V T Q = D V^{\mathsf{T}}Q = D V T Q = D V = Q D V = QD V = Q D D − 1 = D D^{-1} = D D − 1 = D G = V T A V = D Q T A Q D = D H D G = V^{\mathsf{T}}AV = DQ^{\mathsf{T}}AQD = DHD G = V T A V = D Q T A Q D = D H D
∎ この定理の意味は、「既約な Hessenberg 形への直交相似変換は、第 1 列を決めれば符号を除いて一意」ということです。QR 反復 H k = Q k T H k − 1 Q k H_k = Q_k^{\mathsf{T}}H_{k-1}Q_k H k = Q k T H k − 1 Q k Q k Q_k Q k Q k Q_k Q k 注意 5.6 ）や Lanczos 法の陰的再開はこの原理で作られています。
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64）／乱数種 20261006
from scipy.linalg import hessenberg
rng = np.random. default_rng ( 20261006 )
lam = np. arange ( 1.0 , n + 1 )
X = rng. standard_normal ( (n, n) )
cases = { " 非対称 " : X @ np. diag ( lam ) @ np.linalg. inv ( 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:" ,
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 ]))
print ( " |h(n,n-1)|: " , " " . join ( f " {s :.1e } " for s in seq ))
非対称: 下の帯の外 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 0 0 定理 4.1 ）。(b) 下の mpmath との照合で、固有値そのものも確かめます。(c) Rayleigh 商シフトの QR 法の ∣ h n , n − 1 ∣ |h_{n,n-1}| ∣ h n , n − 1 ∣ 1.8 × 10 − 2 → 6.8 × 10 − 6 → 2.0 × 10 − 12 → 1.6 × 10 − 25 1.8\times10^{-2} \to 6.8\times10^{-6} \to 2.0\times10^{-12} \to 1.6\times10^{-25} 1.8 × 1 0 − 2 → 6.8 × 1 0 − 6 → 2.0 × 1 0 − 12 → 1.6 × 1 0 − 25 1.0 × 10 − 4 → 3.5 × 10 − 13 → 1.6 × 10 − 38 1.0\times10^{-4} \to 3.5\times10^{-13} \to 1.6\times10^{-38} 1.0 × 1 0 − 4 → 3.5 × 1 0 − 13 → 1.6 × 1 0 − 38 系 5.3 ）進みます。最初の数回は大域的な段階で、収束の次数は現れません。
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0／乱数種 20261006
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 回
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: " ,
自作 QR 法と mpmath の差の最大 < 1e-12: True
eigvals と mpmath の差の最大 < 1e-12: True
mpmath と構成値 6,...,1 の差の最大 < 1e-13: True
検算. (a) 構成に使った固有値 6 , … , 1 6, \ldots, 1 6 , … , 1 A A A 10 − 13 10^{-13} 1 0 − 13 numpy.linalg.qr の符号のままの版。符号は A k A_k A k eigvals は、どちらも 40 桁の参照値と 10 − 12 10^{-12} 1 0 − 12 ( 5 / 6 ) 300 ≈ 2 × 10 − 24 (5/6)^{300} \approx 2\times10^{-24} ( 5/6 ) 300 ≈ 2 × 1 0 − 24
シフトなし QR 法の速さ ∣ λ j + 1 / λ j ∣ |\lambda_{j+1}/\lambda_j| ∣ λ j + 1 / λ j ∣ μ I \mu I μ I
定義 5.1 （シフトつき QR 法 ） ↗
A 0 = A A_0 = A A 0 = A μ k ∈ R \mu_k \in \mathbb{R} μ k ∈ R A k − 1 − μ k I A_{k-1} - \mu_kI A k − 1 − μ k I A k − 1 − μ k I = Q k R k A_{k-1} - \mu_kI = Q_kR_k A k − 1 − μ k I = Q k R k R k R_k R k
A k = R k Q k + μ k I A_k = R_kQ_k + \mu_kI A k = R k Q k + μ k I とする。μ k = ( A k − 1 ) n n \mu_k = (A_{k-1})_{nn} μ k = ( A k − 1 ) nn Rayleigh 商シフト という。
A k = Q k T ( A k − 1 − μ k I ) Q k + μ k I = Q k T A k − 1 Q k 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 = Q k T ( A k − 1 − μ k I ) Q k + μ k I = Q k T A k − 1 Q k A k − 1 A_{k-1} A k − 1 A k − 1 − μ k I A_{k-1} - \mu_kI A k − 1 − μ k I 命題 4.2 により A k A_k A k
次の定理で転置が出てくるのは、B = Q R B = QR B = QR Q Q Q B B B A k − 1 A_{k-1} A k − 1 A k − 1 T = A k − 1 A_{k-1}^{\mathsf{T}} = A_{k-1} A k − 1 T = A k − 1
定理 5.2 （シフトつき QR 法と逆反復 ） ✓
シフトつき QR 法の 1 段で A k − 1 − μ k I A_{k-1} - \mu_kI A k − 1 − μ k I q = Q k e n \boldsymbol{q} = Q_k\boldsymbol{e}_n q = Q k e n Q k Q_k Q k
q = ( R k ) n n ( A k − 1 − μ k I ) − T e n \boldsymbol{q} = (R_k)_{nn}\,(A_{k-1} - \mu_kI)^{-\mathsf{T}}\boldsymbol{e}_n q = ( R k ) nn ( A k − 1 − μ k I ) − T e n ( R k ) n n > 0 (R_k)_{nn} > 0 ( R k ) nn > 0 q \boldsymbol{q} q A k − 1 T A_{k-1}^{\mathsf{T}} A k − 1 T μ k \mu_k μ k e n \boldsymbol{e}_n e n θ = q T A k − 1 q \theta = \boldsymbol{q}^{\mathsf{T}}A_{k-1}\boldsymbol{q} θ = q T A k − 1 q ( A k ) n n = θ (A_k)_{nn} = \theta ( A k ) nn = θ A k A_k A k ( ∑ j < n ( A k ) n j 2 ) 1 / 2 = ∥ A k − 1 T q − θ q ∥ 2 . \Bigl(\sum_{j<n}(A_k)_{nj}^2\Bigr)^{1/2} = \|A_{k-1}^{\mathsf{T}}\boldsymbol{q} - \theta\boldsymbol{q}\|_2 . ( j < n ∑ ( A k ) nj 2 ) 1/2 = ∥ A k − 1 T q − θ q ∥ 2 . 証明（定理 5.2）
B = A k − 1 − μ k I = Q k R k B = A_{k-1} - \mu_kI = Q_kR_k B = A k − 1 − μ k I = Q k R k B T = R k T Q k T B^{\mathsf{T}} = R_k^{\mathsf{T}}Q_k^{\mathsf{T}} B T = R k T Q k T B − T = Q k − T R k − T = Q k R k − T B^{-\mathsf{T}} = Q_k^{-\mathsf{T}}R_k^{-\mathsf{T}} = Q_kR_k^{-\mathsf{T}} B − T = Q k − T R k − T = Q k R k − T Q k Q_k Q k Q k = B − T R k T Q_k = B^{-\mathsf{T}}R_k^{\mathsf{T}} Q k = B − T R k T R k T R_k^{\mathsf{T}} R k T R k T e n = ( R k ) n n e n R_k^{\mathsf{T}}\boldsymbol{e}_n = (R_k)_{nn}\boldsymbol{e}_n R k T e n = ( R k ) nn e n Q k e n = ( R k ) n n B − T e n Q_k\boldsymbol{e}_n = (R_k)_{nn}B^{-\mathsf{T}}\boldsymbol{e}_n Q k e n = ( R k ) nn B − T e n ( R k ) n n > 0 (R_k)_{nn} > 0 ( R k ) nn > 0
A k = Q k T A k − 1 Q k A_k = Q_k^{\mathsf{T}}A_{k-1}Q_k A k = Q k T A k − 1 Q k ( n , j ) (n, j) ( n , j ) q T A k − 1 q j = ( A k − 1 T q ) T q j \boldsymbol{q}^{\mathsf{T}}A_{k-1}\boldsymbol{q}_j = (A_{k-1}^{\mathsf{T}}\boldsymbol{q})^{\mathsf{T}}\boldsymbol{q}_j q T A k − 1 q j = ( A k − 1 T q ) T q j q j \boldsymbol{q}_j q j Q k Q_k Q k j j j j = n j = n j = n θ \theta θ s = A k − 1 T q − θ q \boldsymbol{s} = A_{k-1}^{\mathsf{T}}\boldsymbol{q} - \theta\boldsymbol{q} s = A k − 1 T q − θ q θ \theta θ q \boldsymbol{q} q A k − 1 T A_{k-1}^{\mathsf{T}} A k − 1 T s ⊥ q \boldsymbol{s} \perp \boldsymbol{q} s ⊥ q 命題 7.1［固有値の計算 I］ ）。j < n j < n j < n q T q j = 0 \boldsymbol{q}^{\mathsf{T}}\boldsymbol{q}_j = 0 q T q j = 0 ( A k ) n j = s T q j (A_k)_{nj} = \boldsymbol{s}^{\mathsf{T}}\boldsymbol{q}_j ( A k ) nj = s T q j q 1 , … , q n \boldsymbol{q}_1, \ldots, \boldsymbol{q}_n q 1 , … , q n s T q n = 0 \boldsymbol{s}^{\mathsf{T}}\boldsymbol{q}_n = 0 s T q n = 0 ∥ s ∥ 2 2 = ∑ j < n ( s T q j ) 2 \|\boldsymbol{s}\|_2^2 = \sum_{j<n}(\boldsymbol{s}^{\mathsf{T}}\boldsymbol{q}_j)^2 ∥ s ∥ 2 2 = ∑ j < n ( s T q j ) 2
∎ つまり、シフトつき QR 法は「最後の列について逆反復を行い、最下行の非対角部分にその残差を記録する」算法です。Rayleigh 商シフト μ k = ( A k − 1 ) n n = e n T A k − 1 e n \mu_k = (A_{k-1})_{nn} = \boldsymbol{e}_n^{\mathsf{T}}A_{k-1}\boldsymbol{e}_n μ k = ( A k − 1 ) nn = e n T A k − 1 e n e n \boldsymbol{e}_n e n Rayleigh 商反復 の 1 段になります。対称行列ではこれを元の座標に戻すと、正確に Rayleigh 商反復の列になります。
次の系は、対称行列では最下行の非対角部分が 1 回ごとに前の 3 乗にほぼ比例して小さくなる、という主張です。式の定数 32 2 Δ 2 / δ 4 32\sqrt{2}\,\Delta^2/\delta^4 32 2 Δ 2 / δ 4 η k + 1 \eta_{k+1} η k + 1 η k \eta_k η k δ \delta δ Δ \Delta Δ
系 5.3 （Rayleigh 商シフトの局所 3 次収束（対称行列） ） ✓
A A A λ \lambda λ v \boldsymbol{v} v δ \delta δ Δ \Delta Δ 定理 5.2［固有値の計算 I］ と同じく定める。Rayleigh 商シフトの QR 法で、x k = Q k e n \boldsymbol{x}_k = \mathbf{Q}_k\boldsymbol{e}_n x k = Q k e n Q k = Q 1 ⋯ Q k \mathbf{Q}_k = Q_1\cdots Q_k Q k = Q 1 ⋯ Q k θ k = ∠ ( x k , v ) \theta_k = \angle(\boldsymbol{x}_k, \boldsymbol{v}) θ k = ∠ ( x k , v ) η k = ( ∑ j < n ( A k ) n j 2 ) 1 / 2 \eta_k = \bigl(\sum_{j<n}(A_k)_{nj}^2\bigr)^{1/2} η k = ( ∑ j < n ( A k ) nj 2 ) 1/2
μ k \mu_k μ k A A A x k \boldsymbol{x}_k x k x 0 = e n \boldsymbol{x}_0 = \boldsymbol{e}_n x 0 = e n 定義 5.1［固有値の計算 I］ ）の k k k η k = ∥ A x k − R A ( x k ) x k ∥ 2 \eta_k = \|A\boldsymbol{x}_k - R_A(\boldsymbol{x}_k)\boldsymbol{x}_k\|_2 η k = ∥ A x k − R A ( x k ) x k ∥ 2 Δ sin 2 θ 0 < δ / 2 \Delta\sin^2\theta_0 < \delta/2 Δ sin 2 θ 0 < δ /2 δ 2 sin θ k ≤ η k ≤ Δ sin θ k \tfrac{\delta}{2}\sin\theta_k \le \eta_k \le \Delta\sin\theta_k 2 δ sin θ k ≤ η k ≤ Δ sin θ k η k + 1 ≤ 32 2 Δ 2 δ 4 η k 3 . \eta_{k+1} \le \frac{32\sqrt{2}\,\Delta^2}{\delta^4}\,\eta_k^3 . η k + 1 ≤ δ 4 32 2 Δ 2 η k 3 . 証明（系 5.3）
帰納法です。k = 0 k = 0 k = 0 Q 0 = I \mathbf{Q}_0 = I Q 0 = I x k − 1 = Q k − 1 e n \boldsymbol{x}_{k-1} = \mathbf{Q}_{k-1}\boldsymbol{e}_n x k − 1 = Q k − 1 e n k − 1 k - 1 k − 1 A k − 1 = Q k − 1 T A Q k − 1 A_{k-1} = \mathbf{Q}_{k-1}^{\mathsf{T}}A\mathbf{Q}_{k-1} A k − 1 = Q k − 1 T A Q k − 1 μ k = e n T A k − 1 e n = R A ( x k − 1 ) \mu_k = \boldsymbol{e}_n^{\mathsf{T}}A_{k-1}\boldsymbol{e}_n = R_A(\boldsymbol{x}_{k-1}) μ k = e n T A k − 1 e n = R A ( x k − 1 ) A k − 1 A_{k-1} A k − 1 定理 5.2 の 1 から
x k = Q k − 1 Q k e n = ( R k ) n n Q k − 1 ( A k − 1 − μ k I ) − 1 e n = ( R k ) n n ( A − μ k I ) − 1 x k − 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} x k = Q k − 1 Q k e n = ( R k ) nn Q k − 1 ( A k − 1 − μ k I ) − 1 e n = ( R k ) nn ( A − μ k I ) − 1 x k − 1 です（Q k − 1 ( Q k − 1 T ( A − μ k I ) Q k − 1 ) − 1 = ( A − μ k I ) − 1 Q 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} Q k − 1 ( Q k − 1 T ( A − μ k I ) Q k − 1 ) − 1 = ( A − μ k I ) − 1 Q k − 1 x k \boldsymbol{x}_k x k ( R k ) n n > 0 (R_k)_{nn} > 0 ( R k ) nn > 0 定理 5.2 の 2 の右辺が、q = Q k e n \boldsymbol{q} = Q_k\boldsymbol{e}_n q = Q k e n Q k − 1 \mathbf{Q}_{k-1} Q k − 1 ∥ A x k − R A ( x k ) x k ∥ 2 \|A\boldsymbol{x}_k - R_A(\boldsymbol{x}_k)\boldsymbol{x}_k\|_2 ∥ A x k − R A ( x k ) x k ∥ 2 Q k − 1 \mathbf{Q}_{k-1} Q k − 1
A A A u i \boldsymbol{u}_i u i λ i \lambda_i λ i v \boldsymbol{v} v ± u i 0 \pm\boldsymbol{u}_{i_0} ± u i 0 x k = ∑ i c i u i \boldsymbol{x}_k = \sum_ic_i\boldsymbol{u}_i x k = ∑ i c i u i ρ = R A ( x k ) \rho = R_A(\boldsymbol{x}_k) ρ = R A ( x k ) 定理 5.2［固有値の計算 I］ から、反復が続く限り Δ sin 2 θ k < δ / 2 \Delta\sin^2\theta_k < \delta/2 Δ sin 2 θ k < δ /2 ∣ ρ − λ ∣ ≤ Δ sin 2 θ k < δ / 2 |\rho - \lambda| \le \Delta\sin^2\theta_k < \delta/2 ∣ ρ − λ ∣ ≤ Δ sin 2 θ k < δ /2 η k \eta_k η k 命題 7.1［固有値の計算 I］ ）から η k 2 ≤ ∥ ( A − λ I ) x k ∥ 2 2 = ∑ i ( λ i − λ ) 2 c i 2 ≤ Δ 2 ∑ i ≠ i 0 c i 2 = Δ 2 sin 2 θ 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 η k 2 ≤ ∥ ( A − λ I ) x k ∥ 2 2 = i ∑ ( λ i − λ ) 2 c i 2 ≤ Δ 2 i = i 0 ∑ c i 2 = Δ 2 sin 2 θ k です。下からは、i ≠ i 0 i \ne i_0 i = i 0 ∣ λ i − ρ ∣ ≥ δ − ∣ λ − ρ ∣ > δ / 2 |\lambda_i - \rho| \ge \delta - |\lambda - \rho| > \delta/2 ∣ λ i − ρ ∣ ≥ δ − ∣ λ − ρ ∣ > δ /2
η k 2 = ∑ i ( λ i − ρ ) 2 c i 2 ≥ δ 2 4 ∑ i ≠ i 0 c i 2 = δ 2 4 sin 2 θ k \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 η k 2 = i ∑ ( λ i − ρ ) 2 c i 2 ≥ 4 δ 2 i = i 0 ∑ c i 2 = 4 δ 2 sin 2 θ k です。最後に、sin 2 θ k < δ / ( 2 Δ ) ≤ 1 / 2 \sin^2\theta_k < \delta/(2\Delta) \le 1/2 sin 2 θ k < δ / ( 2Δ ) ≤ 1/2 cos 2 θ k > 1 / 2 \cos^2\theta_k > 1/2 cos 2 θ k > 1/2 tan θ k < 2 sin θ k \tan\theta_k < \sqrt{2}\sin\theta_k tan θ k < 2 sin θ k 定理 5.2［固有値の計算 I］ の 2 と合わせて
η k + 1 ≤ Δ sin θ k + 1 ≤ Δ tan θ k + 1 ≤ 2 Δ 2 δ tan 3 θ k ≤ 2 Δ 2 δ ( 2 sin θ k ) 3 ≤ 2 Δ 2 δ ⋅ 2 2 ( 2 η k δ ) 3 = 32 2 Δ 2 δ 4 η k 3 \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 η k + 1 ≤ Δ sin θ k + 1 ≤ Δ tan θ k + 1 ≤ δ 2 Δ 2 tan 3 θ k ≤ δ 2 Δ 2 ( 2 sin θ k ) 3 ≤ δ 2 Δ 2 ⋅ 2 2 ( δ 2 η k ) 3 = δ 4 32 2 Δ 2 η k 3 です。
∎ A A A η k = ∣ ( A k ) n , n − 1 ∣ \eta_k = |(A_k)_{n,n-1}| η k = ∣ ( A k ) n , n − 1 ∣ η k \eta_k η k ( A k ) n n (A_k)_{nn} ( A k ) nn n − 1 n - 1 n − 1 注意 5.3［固有値の計算 I］ ）、上のコードのとおり ∣ h n , n − 1 ∣ |h_{n,n-1}| ∣ h n , n − 1 ∣
次の命題は、副対角成分を 0 0 0 A A A
命題 5.4 （デフレーションの後退誤差 ） ✓
H = Q T A Q H = Q^{\mathsf{T}}AQ H = Q T A Q Q Q Q k k k
∣ h k + 1 , k ∣ ≤ u ( ∣ h k k ∣ + ∣ h k + 1 , k + 1 ∣ ) |h_{k+1,k}| \le u\,(|h_{kk}| + |h_{k+1,k+1}|) ∣ h k + 1 , k ∣ ≤ u ( ∣ h k k ∣ + ∣ h k + 1 , k + 1 ∣ ) とする（u u u h k + 1 , k h_{k+1,k} h k + 1 , k 0 0 0 H + E H + E H + E
∥ E ∥ 2 = ∣ h k + 1 , k ∣ ≤ 2 u ∥ H ∥ 2 = 2 u ∥ A ∥ 2 \|E\|_2 = |h_{k+1,k}| \le 2u\|H\|_2 = 2u\|A\|_2 ∥ E ∥ 2 = ∣ h k + 1 , k ∣ ≤ 2 u ∥ H ∥ 2 = 2 u ∥ A ∥ 2 H + E H + E H + E k k k n − k n - k n − k A + Q E Q T A + QEQ^{\mathsf{T}} A + QE Q T ∥ Q E Q T ∥ 2 = ∥ E ∥ 2 \|QEQ^{\mathsf{T}}\|_2 = \|E\|_2 ∥ QE Q T ∥ 2 = ∥ E ∥ 2 A A A H H H h k + 1 , k h_{k+1,k} h k + 1 , k h k , k + 1 h_{k,k+1} h k , k + 1 0 0 0 H + E ′ H + E' H + E ′ ∥ E ′ ∥ 2 = ∣ h k + 1 , k ∣ \|E'\|_2 = |h_{k+1,k}| ∥ E ′ ∥ 2 = ∣ h k + 1 , k ∣ ∣ λ i ( H + E ′ ) − λ i ( A ) ∣ ≤ 2 u ∥ A ∥ 2 |\lambda_i(H + E') - \lambda_i(A)| \le 2u\|A\|_2 ∣ λ i ( H + E ′ ) − λ i ( A ) ∣ ≤ 2 u ∥ A ∥ 2 A = V Λ V − 1 A = V\Lambda V^{-1} A = V Λ V − 1 H + E H + E H + E ν \nu ν min i ∣ ν − λ i ∣ ≤ 2 u κ 2 ( V ) ∥ A ∥ 2 \min_i|\nu - \lambda_i| \le 2u\,\kappa_2(V)\|A\|_2 min i ∣ ν − λ i ∣ ≤ 2 u κ 2 ( V ) ∥ A ∥ 2 証明（命題 5.4）
E = − h k + 1 , k e k + 1 e k T E = -h_{k+1,k}\boldsymbol{e}_{k+1}\boldsymbol{e}_k^{\mathsf{T}} E = − h k + 1 , k e k + 1 e k T ∥ E x ∥ 2 = ∣ h k + 1 , k ∣ ∣ x k ∣ ≤ ∣ h k + 1 , k ∣ ∥ x ∥ 2 \|E\boldsymbol{x}\|_2 = |h_{k+1,k}||x_k| \le |h_{k+1,k}|\|\boldsymbol{x}\|_2 ∥ E x ∥ 2 = ∣ h k + 1 , k ∣∣ x k ∣ ≤ ∣ h k + 1 , k ∣∥ x ∥ 2 x = e k \boldsymbol{x} = \boldsymbol{e}_k x = e k ∣ h i i ∣ = ∣ e i T H e i ∣ ≤ ∥ H ∥ 2 |h_{ii}| = |\boldsymbol{e}_i^{\mathsf{T}}H\boldsymbol{e}_i| \le \|H\|_2 ∣ h ii ∣ = ∣ e i T H e i ∣ ≤ ∥ H ∥ 2 定理 4.2［内積空間とグラム・シュミット直交化］ ）から ∥ E ∥ 2 ≤ 2 u ∥ H ∥ 2 \|E\|_2 \le 2u\|H\|_2 ∥ E ∥ 2 ≤ 2 u ∥ H ∥ 2 ∥ H ∥ 2 = ∥ A ∥ 2 \|H\|_2 = \|A\|_2 ∥ H ∥ 2 = ∥ A ∥ 2 H + E H + E H + E 補題 4.1［行列の固有値評価と不等式］ ）により 2 つの対角ブロックの特性多項式の積です。H + E = Q T ( A + Q E Q T ) Q H + E = Q^{\mathsf{T}}(A + QEQ^{\mathsf{T}})Q H + E = Q T ( A + QE Q T ) Q A + Q E Q T A + QEQ^{\mathsf{T}} A + QE Q T
E ′ E' E ′ ( k + 1 , k ) (k+1, k) ( k + 1 , k ) ( k , k + 1 ) (k, k+1) ( k , k + 1 ) − h k + 1 , k -h_{k+1,k} − h k + 1 , k ∥ E ′ x ∥ 2 2 = h k + 1 , k 2 ( x k 2 + x k + 1 2 ) ≤ h k + 1 , k 2 ∥ x ∥ 2 2 \|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 ′ x ∥ 2 2 = h k + 1 , k 2 ( x k 2 + x k + 1 2 ) ≤ h k + 1 , k 2 ∥ x ∥ 2 2 ∥ E ′ ∥ 2 = ∣ h k + 1 , k ∣ \|E'\|_2 = |h_{k+1,k}| ∥ E ′ ∥ 2 = ∣ h k + 1 , k ∣ H H H A A A H H H E ′ E' E ′ 系 8.1［行列のノルムと線形方程式の条件数］ ）と 1 の評価から従います。
H H H A A A H = ( Q T V ) Λ ( Q T V ) − 1 H = (Q^{\mathsf{T}}V)\Lambda(Q^{\mathsf{T}}V)^{-1} H = ( Q T V ) Λ ( Q T V ) − 1 κ 2 ( Q T V ) = κ 2 ( V ) \kappa_2(Q^{\mathsf{T}}V) = \kappa_2(V) κ 2 ( Q T V ) = κ 2 ( V ) 定理 6.2［行列計算の認証と最適化の証明書］ ）と 1 から従います。
∎ 対称行列では、デフレーションの誤差は固有値の絶対誤差として 2 u ∥ A ∥ 2 2u\|A\|_2 2 u ∥ A ∥ 2 κ 2 ( V ) \kappa_2(V) κ 2 ( V )
演習 6.1 標準
A = ( 2 1 1 1 ) A = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix} A = ( 2 1 1 1 ) μ 1 = a 22 = 1 \mu_1 = a_{22} = 1 μ 1 = a 22 = 1 A 1 A_1 A 1 A A A ( 3 ± 5 ) / 2 (3 \pm \sqrt{5})/2 ( 3 ± 5 ) /2 e 2 \boldsymbol{e}_2 e 2 θ \theta θ A 1 A_1 A 1 θ ′ \theta' θ ′ tan θ ′ = tan 3 θ \tan\theta' = \tan^3\theta tan θ ′ = tan 3 θ
解答 A − I = ( 1 1 1 0 ) A - I = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} A − I = ( 1 1 1 0 ) q 1 = ( 1 , 1 ) T / 2 \boldsymbol{q}_1 = (1, 1)^{\mathsf{T}}/\sqrt{2} q 1 = ( 1 , 1 ) T / 2 r 11 = 2 r_{11} = \sqrt{2} r 11 = 2 r 12 = q 1 T ( 1 , 0 ) T = 1 / 2 r_{12} = \boldsymbol{q}_1^{\mathsf{T}}(1, 0)^{\mathsf{T}} = 1/\sqrt{2} r 12 = q 1 T ( 1 , 0 ) T = 1/ 2 r 22 > 0 r_{22} > 0 r 22 > 0 q 2 = ( 1 , − 1 ) T / 2 \boldsymbol{q}_2 = (1, -1)^{\mathsf{T}}/\sqrt{2} q 2 = ( 1 , − 1 ) T / 2 r 22 = 1 / 2 r_{22} = 1/\sqrt{2} r 22 = 1/ 2
R Q = ( 2 1 / 2 0 1 / 2 ) 1 2 ( 1 1 1 − 1 ) = ( 3 / 2 1 / 2 1 / 2 − 1 / 2 ) , A 1 = ( 5 / 2 1 / 2 1 / 2 1 / 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} . R Q = ( 2 0 1/ 2 1/ 2 ) 2 1 ( 1 1 1 − 1 ) = ( 3/2 1/2 1/2 − 1/2 ) , A 1 = ( 5/2 1/2 1/2 1/2 ) . 跡 3 3 3 1 1 1 λ = ( 3 − 5 ) / 2 \lambda = (3 - \sqrt{5})/2 λ = ( 3 − 5 ) /2 A A A ( 1 , λ − 2 ) (1, \lambda - 2) ( 1 , λ − 2 ) tan θ = 1 / ( 2 − λ ) = 2 / ( 1 + 5 ) ≈ 0.618 \tan\theta = 1/(2 - \lambda) = 2/(1 + \sqrt{5}) \approx 0.618 tan θ = 1/ ( 2 − λ ) = 2/ ( 1 + 5 ) ≈ 0.618 A 1 A_1 A 1 ( A 1 − λ I ) v = 0 (A_1 - \lambda I)\boldsymbol{v} = \boldsymbol{0} ( A 1 − λ I ) v = 0 v ∝ ( 1 / 2 , λ − 5 / 2 ) \boldsymbol{v} \propto (1/2, \lambda - 5/2) v ∝ ( 1/2 , λ − 5/2 ) tan θ ′ = ( 1 / 2 ) / ( 5 / 2 − λ ) = 1 / ( 2 + 5 ) ≈ 0.236 \tan\theta' = (1/2)/(5/2 - \lambda) = 1/(2 + \sqrt{5}) \approx 0.236 tan θ ′ = ( 1/2 ) / ( 5/2 − λ ) = 1/ ( 2 + 5 ) ≈ 0.236 tan 3 θ = 8 / ( 1 + 5 ) 3 = 8 / ( 16 + 8 5 ) = 1 / ( 2 + 5 ) \tan^3\theta = 8/(1 + \sqrt{5})^3 = 8/(16 + 8\sqrt{5}) = 1/(2 + \sqrt{5}) tan 3 θ = 8/ ( 1 + 5 ) 3 = 8/ ( 16 + 8 5 ) = 1/ ( 2 + 5 ) tan θ ′ = tan 3 θ \tan\theta' = \tan^3\theta tan θ ′ = tan 3 θ
演習 6.2 易
定理 4.3 で H H H
解答 A = diag ( 1 , 2 , 3 ) A = \operatorname{diag}(1, 2, 3) A = diag ( 1 , 2 , 3 ) Q = I Q = I Q = I V V V e 1 \boldsymbol{e}_1 e 1 H = A H = A H = A G = V T A V = diag ( 1 , 3 , 2 ) G = V^{\mathsf{T}}AV = \operatorname{diag}(1, 3, 2) G = V T A V = diag ( 1 , 3 , 2 ) V ≠ Q D V \ne QD V = Q D D D D ± 1 \pm1 ± 1 H H H 0 0 0 R H R_H R H
演習 6.3 標準
n n n R Q RQ R Q ( 4 / 3 ) n 3 + n 3 (4/3)n^3 + n^3 ( 4/3 ) n 3 + n 3 6 n 2 6n^2 6 n 2
解答 (a) 各反復は約 ( 7 / 3 ) m 3 (7/3)m^3 ( 7/3 ) m 3 m m m m = n , n − 1 , … , 1 m = n, n-1, \ldots, 1 m = n , n − 1 , … , 1 3 ⋅ ( 7 / 3 ) ∑ m 3 ≈ ( 7 / 4 ) n 4 3\cdot(7/3)\sum m^3 \approx (7/4)n^4 3 ⋅ ( 7/3 ) ∑ m 3 ≈ ( 7/4 ) n 4 ( 10 / 3 ) n 3 (10/3)n^3 ( 10/3 ) n 3 3 ⋅ 6 ∑ m 2 ≈ 6 n 3 3\cdot 6\sum m^2 \approx 6n^3 3 ⋅ 6 ∑ m 2 ≈ 6 n 3 ( 28 / 3 ) n 3 ≈ 9.3 n 3 (28/3)n^3 \approx 9.3n^3 ( 28/3 ) n 3 ≈ 9.3 n 3 Q Q Q n = 1000 n = 1000 n = 1000 1.8 × 10 12 1.8\times10^{12} 1.8 × 1 0 12 9.3 × 10 9 9.3\times10^9 9.3 × 1 0 9
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.
執筆: Claude Code（Opus 5.5）、Claude Code（Sonnet 5.5）・監査 2 回（Claude Code）・最終 2026-10-12