コンテンツにスキップ

固有値の計算 I:冪乗法・逆反復・Rayleigh 商反復と同時反復

対象読者: 学部の線形代数の対角化とジョルダン標準形、実対称行列のスペクトル定理を学んだ人。固有値を数値的に求める算法は仮定せず、ここから導入します。

前提となる記事:スペクトル定理:対称行列はなぜ「回転だけ」で対角化できるのか、対角化とジョルダン標準形:変換が一番簡単に見える座標を探す

生 Markdown

行列の固有値は、行列が「ベクトルをどれだけ伸ばす向きがあるか」を表す数です。小さな行列なら方程式を解いて求められますが、大きな行列では、適当なベクトルに行列を何度も掛けるだけで、いちばん強く伸びる向きが浮かび上がってきます。この記事では、この単純な繰り返しがなぜ固有値と固有ベクトルに行き着くのか、どれくらいの速さで近づくのか、そして掛ける前に行列を少しずらす工夫でどれだけ速くなるのかを、数値例とともに確かめます。

  • 冪乗法は xk+1=Axk/∥Axk∥\boldsymbol{x}_{k+1} = A\boldsymbol{x}_k/\|A\boldsymbol{x}_k\| を繰り返すだけの算法です。対角化可能で ∣λ1∣>∣λ2∣|\lambda_1| > |\lambda_2| なら、xk\boldsymbol{x}_k と固有ベクトル v1\boldsymbol{v}_1 のなす角は O(∣λ2/λ1∣k)O(|\lambda_2/\lambda_1|^k) で 00 に近づきます。
  • Hermite 行列では、Rayleigh 商の誤差は角の正弦の 2 乗に比例します(RA(x)−λ=sin⁡2θ (RA(w)−λ)R_A(\boldsymbol{x}) - \lambda = \sin^2\theta\,(R_A(\boldsymbol{w}) - \lambda))。このため固有値は ∣λ2/λ1∣2k|\lambda_2/\lambda_1|^{2k} の速さで求まります。
  • シフトつき逆反復は (A−μI)−1(A - \mu I)^{-1} に対する冪乗法で、収束率は ∣λJ−μ∣/∣λK−μ∣|\lambda_J - \mu|/|\lambda_K - \mu|(λJ\lambda_J は μ\mu に最も近い固有値、λK\lambda_K は 2 番目に近い固有値)です。シフトを固有値に近づけるほど速くなります。
  • シフトを毎回 Rayleigh 商に取り替える Rayleigh 商反復は、Hermite 行列では角の正接が tan⁡θk+1≤(2Δ/δ)tan⁡3θk\tan\theta_{k+1} \le (2\Delta/\delta)\tan^3\theta_k を満たし、局所 3 次収束(誤差が 1 回ごとにおよそ 3 乗になる速さの収束)します。
  • pp 本のベクトルを同時に掛けて直交化する同時反復は、優勢な pp 次元の不変部分空間(行列を掛けても外へ出ない部分空間)に O(∣λp+1/λp∣k)O(|\lambda_{p+1}/\lambda_p|^k) で近づきます。これが次の章の QR 法の収束の理論の土台になります。
  • 計算した近似固有対は、残差 r=Ax−RA(x)x\boldsymbol{r} = A\boldsymbol{x} - R_A(\boldsymbol{x})\boldsymbol{x} を使って事後に評価できます(残差による評価・Temple の不等式・Bauer–Fike の定理。いずれも、残差の大きさから固有値の誤差を押さえる不等式です)。

研究で使う固有値の計算のほとんどは、numpy.linalg.eigh や scipy.sparse.linalg.eigsh を呼ぶだけで済みます。それでも「何回反復すれば何桁合うのか」「固有値が近いと何が起きるのか」「求めた固有値はどこまで信用できるのか」を答えるには、中で動いている反復の理論が要ります。

固有値問題は多項式の根を求める問題を含むので、5 次以上では有限回の四則と冪根で解く公式がありません。したがって固有値の算法は必ず反復法です。その最も単純な形が冪乗法で、QR 法・Lanczos 法・Arnoldi 法といった実用の算法は、どれも冪乗法を賢く組み合わせたものとして理解できます。

本章では、固有ベクトルを含む方向へベクトルを寄せていく反復を、次の 4 段で扱います。

  1. 冪乗法:Akx0A^k\boldsymbol{x}_0 の方向は、絶対値最大の固有値の固有ベクトルへ向かう。
  2. 逆反復:(A−μI)−1(A - \mu I)^{-1} に冪乗法を使えば、μ\mu に最も近い固有値が「最大」になる。
  3. Rayleigh 商反復:シフト μ\mu を毎回いまの近似固有値に更新すると、収束は 3 次になる。
  4. 同時反復:1 本のベクトルではなく部分空間を反復すると、上位 pp 個の固有値に対応する不変部分空間が求まる。

証明はどれも「固有ベクトルの基底で展開し、各成分が固有値の冪で伸び縮みする様子を見る」という同じ考え方です。

2. 準備:角・Rayleigh 商・2 ノルム

Section titled “2. 準備:角・Rayleigh 商・2 ノルム”

この章では Cn\mathbb{C}^n の標準内積 y∗x\boldsymbol{y}^{*}\boldsymbol{x} とそのノルム ∥x∥=(x∗x)1/2\|\boldsymbol{x}\| = (\boldsymbol{x}^{*}\boldsymbol{x})^{1/2} だけを使います。行列 B∈Cm×nB \in \mathbb{C}^{m\times n} のノルムは、このノルムから誘導される作用素ノルム

∥B∥=sup⁡x≠0∥Bx∥∥x∥\|B\| = \sup_{\boldsymbol{x}\ne\boldsymbol{0}} \frac{\|B\boldsymbol{x}\|}{\|\boldsymbol{x}\|}

です。BB の第 ii 行を bi∗\boldsymbol{b}_i^{*} とすると、Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∥Bx∥2=∑i∣bi∗x∣2≤(∑i∥bi∥2)∥x∥2\|B\boldsymbol{x}\|^2 = \sum_i |\boldsymbol{b}_i^{*}\boldsymbol{x}|^2 \le (\sum_i \|\boldsymbol{b}_i\|^2)\|\boldsymbol{x}\|^2 なので、この上限は有限です。定義から直ちに ∥Bx∥≤∥B∥∥x∥\|B\boldsymbol{x}\| \le \|B\|\|\boldsymbol{x}\|、∥BC∥≤∥B∥∥C∥\|BC\| \le \|B\|\|C\| が従います。対角行列 D=diag⁡(d1,…,dn)D = \operatorname{diag}(d_1, \ldots, d_n) については ∥Dx∥2=∑i∣di∣2∣xi∣2≤max⁡i∣di∣2∥x∥2\|D\boldsymbol{x}\|^2 = \sum_i |d_i|^2|x_i|^2 \le \max_i |d_i|^2 \|\boldsymbol{x}\|^2 で、最大の ∣di∣|d_i| を与える座標ベクトルで等号が成り立つので ∥D∥=max⁡i∣di∣\|D\| = \max_i |d_i| です。

以下の定義は、反復が「固有ベクトルの方向(1 本のときは直線、同時反復のときは部分空間)にどれだけ近づいたか」を測るための物差しです。1 本のベクトルには角を、部分空間どうしには、片方のどのベクトルももう一方に近いかどうかを表す隔たりを使います。

定義 2.1(ベクトルどうし・ベクトルと部分空間のなす角)↗

x,v∈Cn\boldsymbol{x}, \boldsymbol{v} \in \mathbb{C}^n を零でないベクトルとする。v\boldsymbol{v} への直交射影を Pvx=v∗xv∗vvP_{\boldsymbol{v}}\boldsymbol{x} = \dfrac{\boldsymbol{v}^{*}\boldsymbol{x}}{\boldsymbol{v}^{*}\boldsymbol{v}}\boldsymbol{v} とし、

cos⁡∠(x,v)=∣v∗x∣∥x∥∥v∥,sin⁡∠(x,v)=∥x−Pvx∥∥x∥\cos\angle(\boldsymbol{x}, \boldsymbol{v}) = \frac{|\boldsymbol{v}^{*}\boldsymbol{x}|}{\|\boldsymbol{x}\|\|\boldsymbol{v}\|},\qquad \sin\angle(\boldsymbol{x}, \boldsymbol{v}) = \frac{\|\boldsymbol{x} - P_{\boldsymbol{v}}\boldsymbol{x}\|}{\|\boldsymbol{x}\|}

と定める。cos⁡∠≠0\cos\angle \ne 0 のとき tan⁡∠=sin⁡∠/cos⁡∠\tan\angle = \sin\angle/\cos\angle とする。角 ∠(x,v)∈[0,π/2]\angle(\boldsymbol{x}, \boldsymbol{v}) \in [0, \pi/2] はこの正弦・余弦をもつ角である。

部分空間 S≠{0}\mathcal{S} \ne \{\boldsymbol{0}\} への直交射影を PSP_{\mathcal{S}} とするとき(定理 7.1[内積空間とグラム・シュミット直交化] により存在する)、sin⁡∠(x,S)=∥x−PSx∥/∥x∥\sin\angle(\boldsymbol{x}, \mathcal{S}) = \|\boldsymbol{x} - P_{\mathcal{S}}\boldsymbol{x}\|/\|\boldsymbol{x}\| とおき、部分空間 U\mathcal{U} から S\mathcal{S} への隔たりを

d(U,S)=sup⁡x∈U, x≠0sin⁡∠(x,S)d(\mathcal{U}, \mathcal{S}) = \sup_{\boldsymbol{x}\in\mathcal{U},\ \boldsymbol{x}\ne\boldsymbol{0}} \sin\angle(\boldsymbol{x}, \mathcal{S})

で定める。

d(U,S)d(\mathcal{U}, \mathcal{S}) が小さいとは、「U\mathcal{U} のどの単位ベクトルも S\mathcal{S} のすぐ近くにある」ということです。

補題 2.2(角の基本性質)✓

x,v\boldsymbol{x}, \boldsymbol{v} を零でないベクトルとする。

  1. cos⁡2∠(x,v)+sin⁡2∠(x,v)=1\cos^2\angle(\boldsymbol{x}, \boldsymbol{v}) + \sin^2\angle(\boldsymbol{x}, \boldsymbol{v}) = 1。
  2. sin⁡∠(x,v)=min⁡α∈C∥x−αv∥/∥x∥\sin\angle(\boldsymbol{x}, \boldsymbol{v}) = \min_{\alpha\in\mathbb{C}} \|\boldsymbol{x} - \alpha\boldsymbol{v}\|/\|\boldsymbol{x}\|。部分空間 S\mathcal{S} についても sin⁡∠(x,S)=min⁡s∈S∥x−s∥/∥x∥\sin\angle(\boldsymbol{x}, \mathcal{S}) = \min_{\boldsymbol{s}\in\mathcal{S}} \|\boldsymbol{x} - \boldsymbol{s}\|/\|\boldsymbol{x}\|。
  3. 零でない複素数 a,ba, b について ∠(ax,bv)=∠(x,v)\angle(a\boldsymbol{x}, b\boldsymbol{v}) = \angle(\boldsymbol{x}, \boldsymbol{v})。
  4. x=av+z\boldsymbol{x} = a\boldsymbol{v} + \boldsymbol{z}、v∗z=0\boldsymbol{v}^{*}\boldsymbol{z} = 0、a≠0a \ne 0 と書けるとき tan⁡∠(x,v)=∥z∥/(∣a∣∥v∥)\tan\angle(\boldsymbol{x}, \boldsymbol{v}) = \|\boldsymbol{z}\|/(|a|\|\boldsymbol{v}\|)。
証明(補題 2.2)

p=Pvx\boldsymbol{p} = P_{\boldsymbol{v}}\boldsymbol{x}、z=x−p\boldsymbol{z} = \boldsymbol{x} - \boldsymbol{p} とおくと v∗z=v∗x−v∗x=0\boldsymbol{v}^{*}\boldsymbol{z} = \boldsymbol{v}^{*}\boldsymbol{x} - \boldsymbol{v}^{*}\boldsymbol{x} = 0 です。

  1. p\boldsymbol{p} は v\boldsymbol{v} の定数倍なので p∗z=0\boldsymbol{p}^{*}\boldsymbol{z} = 0 で、∥x∥2=∥p∥2+∥z∥2\|\boldsymbol{x}\|^2 = \|\boldsymbol{p}\|^2 + \|\boldsymbol{z}\|^2 です。∥p∥=∣v∗x∣/∥v∥\|\boldsymbol{p}\| = |\boldsymbol{v}^{*}\boldsymbol{x}|/\|\boldsymbol{v}\| なので、両辺を ∥x∥2\|\boldsymbol{x}\|^2 で割れば主張になります。
  2. 任意の α\alpha について x−αv=z+(p−αv)\boldsymbol{x} - \alpha\boldsymbol{v} = \boldsymbol{z} + (\boldsymbol{p} - \alpha\boldsymbol{v}) で、第 2 項は v\boldsymbol{v} の定数倍なので z\boldsymbol{z} と直交します。よって ∥x−αv∥2=∥z∥2+∥p−αv∥2≥∥z∥2\|\boldsymbol{x} - \alpha\boldsymbol{v}\|^2 = \|\boldsymbol{z}\|^2 + \|\boldsymbol{p} - \alpha\boldsymbol{v}\|^2 \ge \|\boldsymbol{z}\|^2 で、等号は αv=p\alpha\boldsymbol{v} = \boldsymbol{p} で成り立ちます。部分空間の場合も、x−PSx\boldsymbol{x} - P_{\mathcal{S}}\boldsymbol{x} が S\mathcal{S} に直交すること(定理 7.1[内積空間とグラム・シュミット直交化])から同じ計算で示せます。
  3. 余弦の定義の分子・分母がともに ∣a∣∣b∣|a||b| 倍になるので余弦は変わらず、1 により正弦も変わりません。
  4. 分解 x=p+z\boldsymbol{x} = \boldsymbol{p} + \boldsymbol{z} は一意です(v\boldsymbol{v} の倍数と v\boldsymbol{v} に直交するベクトルへの分解)。よって p=av\boldsymbol{p} = a\boldsymbol{v} で、tan⁡∠=∥z∥/∥p∥=∥z∥/(∣a∣∥v∥)\tan\angle = \|\boldsymbol{z}\|/\|\boldsymbol{p}\| = \|\boldsymbol{z}\|/(|a|\|\boldsymbol{v}\|) です。
vwxθcos θ(固有ベクトルの成分)sin θ
固有ベクトル v の方向と、それに直交する方向 w への分解。x = cos θ · v + sin θ · w と書くと、本章の評価はすべて θ の漸化式になる。

定義 2.3(Rayleigh 商と残差)↗

A∈Cn×nA \in \mathbb{C}^{n\times n}、x≠0\boldsymbol{x} \ne \boldsymbol{0} について

RA(x)=x∗Axx∗xR_A(\boldsymbol{x}) = \frac{\boldsymbol{x}^{*}A\boldsymbol{x}}{\boldsymbol{x}^{*}\boldsymbol{x}}

を x\boldsymbol{x} の Rayleigh 商、r=Ax−RA(x)x\boldsymbol{r} = A\boldsymbol{x} - R_A(\boldsymbol{x})\boldsymbol{x} を残差という。

AA が Hermite 行列なら RA(x)R_A(\boldsymbol{x}) は実数で、最小の固有値と最大の固有値の間にあります(定理 7.1[スペクトル定理]。そこでの R(x)R(\boldsymbol{x}) と同じものです)。

定義 3.1(冪乗法)↗

A∈Cn×nA \in \mathbb{C}^{n\times n} と x0≠0\boldsymbol{x}_0 \ne \boldsymbol{0} に対し、Axk≠0A\boldsymbol{x}_k \ne \boldsymbol{0} である限り

xk+1=Axk∥Axk∥,μk+1=RA(xk+1)\boldsymbol{x}_{k+1} = \frac{A\boldsymbol{x}_k}{\|A\boldsymbol{x}_k\|},\qquad \mu_{k+1} = R_A(\boldsymbol{x}_{k+1})

で xk\boldsymbol{x}_k と固有値の近似 μk\mu_k を作る反復を冪乗法という。

帰納法により xk\boldsymbol{x}_k は Akx0A^k\boldsymbol{x}_0 の正の定数倍です。正規化は桁あふれを防ぐためだけのもので、方向は Akx0A^k\boldsymbol{x}_0 と同じです。1 回の反復の主な手間は行列とベクトルの積 1 回なので、疎行列なら非零成分の個数に比例する手間で済みます。

下の動画は 2 次の対称行列(固有ベクトル v1\boldsymbol{v}_1 の向きは共通)で xk\boldsymbol{x}_k の向きが v1\boldsymbol{v}_1 に寄っていく様子です。比 ∣λ2/λ1∣|\lambda_2/\lambda_1| が大きい左は遅く、小さい右は速く収束します。

単位円上で、冪乗法の反復ごとに矢印が固有ベクトルの方向へ回っていく動画。比 0.8 の左は少しずつ、比 0.4 の右は数回でほぼ一致する。
冪乗法で単位ベクトル x_k(赤い矢印)が固有ベクトル v₁ の方向(緑の破線)へ近づく様子。左は比 0.8、右は比 0.4。薄い点はこれまでの x_j。

定理 3.2(冪乗法の収束)✓

A∈Cn×nA \in \mathbb{C}^{n\times n} は対角化可能で、固有値を重複を込めて ∣λ1∣>∣λ2∣≥∣λ3∣≥⋯≥∣λn∣|\lambda_1| > |\lambda_2| \ge |\lambda_3| \ge \cdots \ge |\lambda_n| と並べられるとする。対応する固有ベクトル v1,…,vn\boldsymbol{v}_1, \ldots, \boldsymbol{v}_n(∥vi∥=1\|\boldsymbol{v}_i\| = 1)は基底をなすとし、x0=∑i=1ncivi\boldsymbol{x}_0 = \sum_{i=1}^{n} c_i\boldsymbol{v}_i で c1≠0c_1 \ne 0 とする。このとき冪乗法はすべての kk で定義され、

r=∣λ2λ1∣,S=∑i=2n∣ci∣∣c1∣r = \left|\frac{\lambda_2}{\lambda_1}\right|,\qquad S = \sum_{i=2}^{n}\frac{|c_i|}{|c_1|}

とおくと、Srk<1Sr^k < 1 を満たす kk について

sin⁡∠(xk,v1)≤Srk1−Srk\sin\angle(\boldsymbol{x}_k, \boldsymbol{v}_1) \le \frac{S r^k}{1 - S r^k}

が成り立つ。とくに sin⁡∠(xk,v1)=O(rk)\sin\angle(\boldsymbol{x}_k, \boldsymbol{v}_1) = O(r^k) である。

証明(定理 3.2)

対角化可能なので固有ベクトルの基底がとれます(定理 3.2[対角化とジョルダン標準形])。Akvi=λikviA^k\boldsymbol{v}_i = \lambda_i^k\boldsymbol{v}_i なので

Akx0=∑i=1nciλikvi.A^k\boldsymbol{x}_0 = \sum_{i=1}^{n} c_i\lambda_i^k\boldsymbol{v}_i .

∣λ1∣>∣λ2∣≥0|\lambda_1| > |\lambda_2| \ge 0 より λ1≠0\lambda_1 \ne 0 で、v1\boldsymbol{v}_1 の係数 c1λ1kc_1\lambda_1^k は 00 でありません。vi\boldsymbol{v}_i は一次独立なので Akx0≠0A^k\boldsymbol{x}_0 \ne \boldsymbol{0} で、冪乗法は止まりません。

yk=Akx0/(c1λ1k)\boldsymbol{y}_k = A^k\boldsymbol{x}_0/(c_1\lambda_1^k) とおくと

yk=v1+wk,wk=∑i=2ncic1(λiλ1)kvi\boldsymbol{y}_k = \boldsymbol{v}_1 + \boldsymbol{w}_k,\qquad \boldsymbol{w}_k = \sum_{i=2}^{n}\frac{c_i}{c_1}\left(\frac{\lambda_i}{\lambda_1}\right)^{k}\boldsymbol{v}_i

で、∣λi/λ1∣≤r|\lambda_i/\lambda_1| \le r(i≥2i \ge 2)と ∥vi∥=1\|\boldsymbol{v}_i\| = 1 から三角不等式で ∥wk∥≤Srk\|\boldsymbol{w}_k\| \le Sr^k です。xk\boldsymbol{x}_k は yk\boldsymbol{y}_k の零でない定数倍なので、補題 2.2 の 3 と 2 から

sin⁡∠(xk,v1)=sin⁡∠(yk,v1)≤∥yk−v1∥∥yk∥=∥wk∥∥v1+wk∥.\sin\angle(\boldsymbol{x}_k, \boldsymbol{v}_1) = \sin\angle(\boldsymbol{y}_k, \boldsymbol{v}_1) \le \frac{\|\boldsymbol{y}_k - \boldsymbol{v}_1\|}{\|\boldsymbol{y}_k\|} = \frac{\|\boldsymbol{w}_k\|}{\|\boldsymbol{v}_1 + \boldsymbol{w}_k\|} .

∥v1+wk∥≥1−∥wk∥≥1−Srk>0\|\boldsymbol{v}_1 + \boldsymbol{w}_k\| \ge 1 - \|\boldsymbol{w}_k\| \ge 1 - Sr^k > 0 なので主張が従います。

収束率 r=∣λ2/λ1∣r = |\lambda_2/\lambda_1| は行列だけで決まり、初期ベクトルは定数 SS にしか効きません。r=0.99r = 0.99 なら 1 桁合わせるのに約 230 回かかるので、冪乗法そのものは「最大固有値が他から十分離れている」ときにしか実用になりません。この遅さを直すのが §4 以降のシフトです。

方向が分かれば固有値も分かります。次の系は一般の(Hermite とは限らない)行列での評価です。

系 3.3(Rayleigh 商による固有値の評価(一般の行列))✓

Av=λvA\boldsymbol{v} = \lambda\boldsymbol{v}、v≠0\boldsymbol{v} \ne \boldsymbol{0} とする。任意の x≠0\boldsymbol{x} \ne \boldsymbol{0} について

∣RA(x)−λ∣≤∥A−λI∥ sin⁡∠(x,v).|R_A(\boldsymbol{x}) - \lambda| \le \|A - \lambda I\|\,\sin\angle(\boldsymbol{x}, \boldsymbol{v}) .

とくに 定理 3.2 の仮定のもとで ∣μk−λ1∣=O(rk)|\mu_k - \lambda_1| = O(r^k) である。

証明(系 3.3)

補題 2.2 の 3 により ∥x∥=1\|\boldsymbol{x}\| = 1 としてよく、Rayleigh 商も x\boldsymbol{x} の定数倍で変わりません。d=x−Pvx\boldsymbol{d} = \boldsymbol{x} - P_{\boldsymbol{v}}\boldsymbol{x} とおくと ∥d∥=sin⁡∠(x,v)\|\boldsymbol{d}\| = \sin\angle(\boldsymbol{x}, \boldsymbol{v}) で、(A−λI)Pvx=0(A - \lambda I)P_{\boldsymbol{v}}\boldsymbol{x} = \boldsymbol{0} です。したがって

RA(x)−λ=x∗(A−λI)x=x∗(A−λI)dR_A(\boldsymbol{x}) - \lambda = \boldsymbol{x}^{*}(A - \lambda I)\boldsymbol{x} = \boldsymbol{x}^{*}(A - \lambda I)\boldsymbol{d}

で、Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∣RA(x)−λ∣≤∥x∥∥(A−λI)d∥≤∥A−λI∥∥d∥|R_A(\boldsymbol{x}) - \lambda| \le \|\boldsymbol{x}\|\|(A - \lambda I)\boldsymbol{d}\| \le \|A - \lambda I\|\|\boldsymbol{d}\| です。後半は 定理 3.2 の評価を代入すれば従います。

Hermite 行列では Rayleigh 商の精度が 1 段上がります。誤差が角の正弦の 1 乗ではなく 2 乗で効くのです。たとえば方向が 0.01 ずれていても、Rayleigh 商のずれは 0.0001 に比例する程度で済みます。次の補題は、ベクトルを固有ベクトルの成分とそれに直交する成分に分け、前者が固有値を正確に与え、後者のずれだけが誤差になることを示します。

補題 3.4(Rayleigh 商の 2 次の精度)✓

AA を Hermite 行列とし、その最大・最小の固有値を λmax⁡\lambda_{\max}, λmin⁡\lambda_{\min} とする。Av=λvA\boldsymbol{v} = \lambda\boldsymbol{v}、∥v∥=1\|\boldsymbol{v}\| = 1 とし、単位ベクトル x\boldsymbol{x} を x=αv+βw\boldsymbol{x} = \alpha\boldsymbol{v} + \beta\boldsymbol{w}(∥w∥=1\|\boldsymbol{w}\| = 1、v∗w=0\boldsymbol{v}^{*}\boldsymbol{w} = 0、∣β∣=sin⁡∠(x,v)|\beta| = \sin\angle(\boldsymbol{x}, \boldsymbol{v}))と分解する。このとき

RA(x)−λ=sin⁡2∠(x,v) (RA(w)−λ)R_A(\boldsymbol{x}) - \lambda = \sin^2\angle(\boldsymbol{x}, \boldsymbol{v})\,\bigl(R_A(\boldsymbol{w}) - \lambda\bigr)

であり、とくに ∣RA(x)−λ∣≤(λmax⁡−λmin⁡)sin⁡2∠(x,v)|R_A(\boldsymbol{x}) - \lambda| \le (\lambda_{\max} - \lambda_{\min})\sin^2\angle(\boldsymbol{x}, \boldsymbol{v}) である。(sin⁡∠(x,v)=0\sin\angle(\boldsymbol{x}, \boldsymbol{v}) = 0 のときは w\boldsymbol{w} を v\boldsymbol{v} に直交する任意の単位ベクトルとする。)

証明(補題 3.4)

分解は αv=Pvx\alpha\boldsymbol{v} = P_{\boldsymbol{v}}\boldsymbol{x}、βw=x−Pvx\beta\boldsymbol{w} = \boldsymbol{x} - P_{\boldsymbol{v}}\boldsymbol{x} として作れ、補題 2.2 の 1 から ∣α∣2+∣β∣2=1|\alpha|^2 + |\beta|^2 = 1 です。AA は Hermite で λ\lambda は実数なので v∗Aw=(Av)∗w=λv∗w=0\boldsymbol{v}^{*}A\boldsymbol{w} = (A\boldsymbol{v})^{*}\boldsymbol{w} = \lambda\boldsymbol{v}^{*}\boldsymbol{w} = 0、同様に w∗Av=0\boldsymbol{w}^{*}A\boldsymbol{v} = 0 です。よって交差項が消えて

RA(x)=x∗Ax=∣α∣2λ+∣β∣2 w∗Aw=(1−∣β∣2)λ+∣β∣2RA(w)R_A(\boldsymbol{x}) = \boldsymbol{x}^{*}A\boldsymbol{x} = |\alpha|^2\lambda + |\beta|^2\,\boldsymbol{w}^{*}A\boldsymbol{w} = (1 - |\beta|^2)\lambda + |\beta|^2 R_A(\boldsymbol{w})

となり、λ\lambda を移項すれば等式を得ます。RA(w)R_A(\boldsymbol{w}) と λ\lambda はともに [λmin⁡,λmax⁡][\lambda_{\min}, \lambda_{\max}] に入る(定理 7.1[スペクトル定理] の 1)ので、∣RA(w)−λ∣≤λmax⁡−λmin⁡|R_A(\boldsymbol{w}) - \lambda| \le \lambda_{\max} - \lambda_{\min} です。

系 3.5(Hermite 行列の冪乗法)✓

AA を Hermite 行列とし、固有値を ∣λ1∣>∣λ2∣≥⋯≥∣λn∣|\lambda_1| > |\lambda_2| \ge \cdots \ge |\lambda_n| と並べ、対応する正規直交な固有ベクトルを u1,…,un\boldsymbol{u}_1, \ldots, \boldsymbol{u}_n とする。u1∗x0≠0\boldsymbol{u}_1^{*}\boldsymbol{x}_0 \ne 0 ならば、冪乗法の xk\boldsymbol{x}_k と μk=RA(xk)\mu_k = R_A(\boldsymbol{x}_k) は、r=∣λ2/λ1∣r = |\lambda_2/\lambda_1| として

tan⁡∠(xk,u1)≤rktan⁡∠(x0,u1),∣μk−λ1∣≤(λmax⁡−λmin⁡) r2ktan⁡2∠(x0,u1)\tan\angle(\boldsymbol{x}_k, \boldsymbol{u}_1) \le r^k\tan\angle(\boldsymbol{x}_0, \boldsymbol{u}_1),\qquad |\mu_k - \lambda_1| \le (\lambda_{\max} - \lambda_{\min})\,r^{2k}\tan^2\angle(\boldsymbol{x}_0, \boldsymbol{u}_1)

を満たす。

証明(系 3.5)

正規直交な固有ベクトルの基底はスペクトル定理(定理 4.2[スペクトル定理])で存在します。ci=ui∗x0c_i = \boldsymbol{u}_i^{*}\boldsymbol{x}_0 とおくと Akx0=∑iciλikuiA^k\boldsymbol{x}_0 = \sum_i c_i\lambda_i^k\boldsymbol{u}_i で、z=∑i≥2ciλikui\boldsymbol{z} = \sum_{i\ge2}c_i\lambda_i^k\boldsymbol{u}_i は u1\boldsymbol{u}_1 に直交します。補題 2.2 の 4 と 3(xk\boldsymbol{x}_k は Akx0A^k\boldsymbol{x}_0 の正の定数倍)と、正規直交基底での三平方の定理から

tan⁡∠(xk,u1)=(∑i≥2∣ci∣2∣λi∣2k)1/2∣c1∣∣λ1∣k≤rk (∑i≥2∣ci∣2)1/2∣c1∣=rktan⁡∠(x0,u1)\tan\angle(\boldsymbol{x}_k, \boldsymbol{u}_1) = \frac{\bigl(\sum_{i\ge2}|c_i|^2|\lambda_i|^{2k}\bigr)^{1/2}}{|c_1||\lambda_1|^k} \le r^k\,\frac{\bigl(\sum_{i\ge2}|c_i|^2\bigr)^{1/2}}{|c_1|} = r^k\tan\angle(\boldsymbol{x}_0, \boldsymbol{u}_1)

です。sin⁡2≤tan⁡2\sin^2 \le \tan^2 と 補題 3.4 から後半が従います。

一般の行列では方向の誤差と固有値の誤差が同じ rkr^k でしたが、Hermite 行列では固有値の誤差が r2kr^{2k} になります。反復の回数で言えば、固有値は固有ベクトルの半分の回数で同じ桁数に達します。

# Python 3.12.10 / NumPy 2.5.2(x86-64, OpenBLAS)/乱数種 20261006
import numpy as np
rng = np.random.default_rng(20261006)
n = 50
lam = np.concatenate(([10.0, 8.0], np.linspace(-6, 6, n - 2)))
Q, _ = np.linalg.qr(rng.standard_normal((n, n)))
A = (Q * lam) @ Q.T # 固有値 lam、固有ベクトル Q の列
v1 = Q[:, 0]
x = rng.standard_normal(n)
x /= np.linalg.norm(x)
print(" k sin(角) |rho-l1| 10 回の比(角) 比(rho)")
prev = None
for k in range(51):
if k % 10 == 0:
s = np.linalg.norm(x - (v1 @ x) * v1)
e = abs(x @ A @ x - 10.0)
r = "" if prev is None else \
f"{s / prev[0]:.2e} {e / prev[1]:.2e}"
print(f"{k:2d} {s:.2e} {e:.2e} {r}")
prev = (s, e)
x = A @ x
x /= np.linalg.norm(x)
print("予測:", f"{0.8**10:.3f}", f"{0.8**20:.4f}")
k sin(角) |rho-l1| 10 回の比(角) 比(rho)
0 9.05e-01 8.11e+00
10 3.71e-03 6.59e-05 4.10e-03 8.14e-06
20 2.53e-04 1.29e-07 6.82e-02 1.95e-03
30 2.71e-05 1.47e-09 1.07e-01 1.14e-02
40 2.91e-06 1.69e-11 1.07e-01 1.15e-02
50 3.12e-07 1.94e-13 1.07e-01 1.14e-02
予測: 0.107 0.0115

検算. (a) 固有値は構成に使った値 10,810, 8 と [−6,6][-6, 6] の等間隔点なので、r=0.8r = 0.8 が既知です。(b) 角は構成に使った固有ベクトル Q[:, 0] との間で測っています。(c) k≥30k \ge 30 で、10 回あたりの比は角で 0.107≈0.8100.107 \approx 0.8^{10}、Rayleigh 商で 0.0114≈0.8200.0114 \approx 0.8^{20} となり、定理 3.2 と 系 3.5 が予言する傾きに一致します。最初の 10 回の比がそれより小さいのは、初期ベクトルの誤差の大半が ∣λi∣≤6|\lambda_i| \le 6 の成分(10 回で 0.610≈0.0060.6^{10} \approx 0.006 倍に縮む)にあったためです。k=50k = 50 での評価式の右辺は (λmax⁡−λmin⁡)r100tan⁡2∠(x0,u1)≈16×2.0×10−10×4.5≈1.5×10−8(\lambda_{\max} - \lambda_{\min})r^{100}\tan^2\angle(\boldsymbol{x}_0, \boldsymbol{u}_1) \approx 16 \times 2.0\times10^{-10} \times 4.5 \approx 1.5\times10^{-8} で、実際の誤差 1.9×10−131.9\times10^{-13} はその内側にあります。

11e-41e-81e-1201020304050sin(角)|Rayleigh 商 − λ₁|
上の表の誤差を片対数で描いたもの(横軸は反復回数 k、縦軸は誤差の常用対数)。角(細線)の傾きは log₁₀ 0.8、Rayleigh 商(太線)の傾きはその 2 倍になる。

仮定 ∣λ1∣>∣λ2∣|\lambda_1| > |\lambda_2| と対角化可能性を外すと何が壊れるかを、2 つの例で見ます。

例 3.6(絶対値の等しい固有値:回転)

A=(0−110)A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}(90° の回転)の固有値は ±i\pm i で、絶対値が等しいので 定理 3.2 の仮定を満たしません。実際、x0=(1,0)T\boldsymbol{x}_0 = (1, 0)^{\mathsf{T}} から始めると xk\boldsymbol{x}_k は (1,0)T,(0,1)T,(−1,0)T,(0,−1)T(1,0)^{\mathsf{T}}, (0,1)^{\mathsf{T}}, (-1,0)^{\mathsf{T}}, (0,-1)^{\mathsf{T}} を周期 4 で回り、どの方向にも収束しません。Rayleigh 商はつねに 00 で、どちらの固有値にも近づきません。実行列の冪乗法は、絶対値最大の固有値が複素共役の対になると必ずこうなります。

例 3.7(Jordan 細胞:1/k の遅い収束)

J=(λ10λ)J = \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix}(λ≠0\lambda \ne 0)は対角化できません(Jordan 標準形 定理 6.2[対角化とジョルダン標準形] の 2 次の細胞そのもの)。Jk=λk(1k/λ01)J^k = \lambda^{k}\begin{pmatrix} 1 & k/\lambda \\ 0 & 1 \end{pmatrix} なので、x0=(0,1)T\boldsymbol{x}_0 = (0, 1)^{\mathsf{T}} に対し Jkx0=λk−1(k,λ)TJ^k\boldsymbol{x}_0 = \lambda^{k-1}(k, \lambda)^{\mathsf{T}} です。唯一の固有ベクトルの方向 e1\boldsymbol{e}_1 とのなす角は tan⁡∠=∣λ∣/k\tan\angle = |\lambda|/k で、確かに 00 に近づきますが、その速さは幾何級数的ではなく 1/k1/k です。1 桁合わせるのに反復回数を 10 倍にしなければなりません。

冪乗法が見つけるのは絶対値最大の固有値だけです。μ\mu が固有値でないとき、(A−μI)−1(A - \mu I)^{-1} の固有値は 1/(λi−μ)1/(\lambda_i - \mu) で、固有ベクトルは AA と同じです。μ\mu に最も近い固有値が、(A−μI)−1(A - \mu I)^{-1} では絶対値最大の固有値になります。

定義 4.1(シフトつき逆反復)↗

μ∈C\mu \in \mathbb{C} が AA の固有値でないとき、x0≠0\boldsymbol{x}_0 \ne \boldsymbol{0} から

(A−μI)yk+1=xk,xk+1=yk+1∥yk+1∥(A - \mu I)\boldsymbol{y}_{k+1} = \boldsymbol{x}_k,\qquad \boldsymbol{x}_{k+1} = \frac{\boldsymbol{y}_{k+1}}{\|\boldsymbol{y}_{k+1}\|}

を繰り返す反復を、シフト μ\mu の逆反復という。

逆行列は作らず、A−μIA - \mu I を 1 回 LU 分解しておいて、毎回は前進・後退代入(O(n2)O(n^2))だけを行います。

定理 4.2(シフトつき逆反復の収束)✓

AA は 定理 3.2 と同じく対角化可能で、固有値 λ1,…,λn\lambda_1, \ldots, \lambda_n(重複を込める)と単位固有ベクトル v1,…,vn\boldsymbol{v}_1, \ldots, \boldsymbol{v}_n が基底をなすとする。μ\mu は固有値でなく、ある添字 JJ について

∣λJ−μ∣<∣λK−μ∣≤∣λi−μ∣(i≠J)|\lambda_J - \mu| < |\lambda_K - \mu| \le |\lambda_i - \mu|\quad (i \ne J)

を満たす K≠JK \ne J があるとする(λJ\lambda_J は μ\mu に最も近い固有値で単純、λK\lambda_K は 2 番目に近い固有値)。x0=∑icivi\boldsymbol{x}_0 = \sum_i c_i\boldsymbol{v}_i で cJ≠0c_J \ne 0 ならば、

r=∣λJ−μ∣∣λK−μ∣,S=∑i≠J∣ci∣∣cJ∣r = \frac{|\lambda_J - \mu|}{|\lambda_K - \mu|},\qquad S = \sum_{i\ne J}\frac{|c_i|}{|c_J|}

として、Srk<1Sr^k < 1 のとき sin⁡∠(xk,vJ)≤Srk/(1−Srk)\sin\angle(\boldsymbol{x}_k, \boldsymbol{v}_J) \le Sr^k/(1 - Sr^k) が成り立つ。AA が Hermite 行列で μ\mu が実数なら、さらに tan⁡∠(xk,vJ)≤rktan⁡∠(x0,vJ)\tan\angle(\boldsymbol{x}_k, \boldsymbol{v}_J) \le r^k\tan\angle(\boldsymbol{x}_0, \boldsymbol{v}_J) である(vi\boldsymbol{v}_i は正規直交にとる)。

証明(定理 4.2)

B=(A−μI)−1B = (A - \mu I)^{-1} とおくと (A−μI)vi=(λi−μ)vi(A - \mu I)\boldsymbol{v}_i = (\lambda_i - \mu)\boldsymbol{v}_i の両辺に BB を掛けて (λi−μ)−1(\lambda_i - \mu)^{-1} で割れば Bvi=(λi−μ)−1viB\boldsymbol{v}_i = (\lambda_i - \mu)^{-1}\boldsymbol{v}_i です。よって BB は同じ固有ベクトルの基底で対角化でき、固有値 νi=(λi−μ)−1\nu_i = (\lambda_i - \mu)^{-1} は仮定から ∣νJ∣>∣νK∣≥∣νi∣|\nu_J| > |\nu_K| \ge |\nu_i|(i≠Ji \ne J)を満たします。逆反復は BB に対する冪乗法そのものなので、添字を JJ が先頭、KK が 2 番目になるよう並べ替えて 定理 3.2 を使えば、比は ∣νK/νJ∣=r|\nu_K/\nu_J| = r で前半が従います。AA が Hermite で μ\mu が実数なら BB も Hermite(B∗=((A−μI)∗)−1=BB^{*} = ((A - \mu I)^{*})^{-1} = B)なので、系 3.5 の前半を BB に使えば後半が従います。

μ\mu を λJ\lambda_J に近づけるほど rr は 00 に近づきます。固有値のよい近似(たとえば二分法や QR 法で求めた値)があれば、逆反復は 2〜3 回で固有ベクトルを機械精度まで求めます。これが LAPACK で三重対角行列(対角とその両隣にしか非零成分がない行列)の固有ベクトルを求める古典的な方法でした。

# Python 3.12.10 / NumPy 2.5.2(x86-64, OpenBLAS)/乱数種 20261006
import numpy as np
rng = np.random.default_rng(20261006)
n = 50
lam = np.concatenate(([10.0, 8.0], np.linspace(-6, 6, n - 2)))
Q, _ = np.linalg.qr(rng.standard_normal((n, n)))
A = (Q * lam) @ Q.T
v = Q[:, 1] # 固有値 8 の固有ベクトル
x0 = rng.standard_normal(n)
def tan_angle(x):
c = v @ x
return np.linalg.norm(x - c * v) / abs(c)
print(" mu 予測の比 t1 t2 t3 t3/t2")
for mu in [8.5, 8.1, 8.01]:
M = A - mu * np.eye(n)
x, ts = x0 / np.linalg.norm(x0), []
for k in range(3):
x = np.linalg.solve(M, x)
x /= np.linalg.norm(x)
ts.append(tan_angle(x))
pred = abs(8 - mu) / abs(10 - mu)
t = " ".join(f"{s:.1e}" for s in ts)
print(f"{mu:5.2f} {pred:.2e} {t} {ts[2] / ts[1]:.2e}")
mu 予測の比 t1 t2 t3 t3/t2
8.50 3.33e-01 1.7e+01 5.2e+00 1.7e+00 3.28e-01
8.10 5.26e-02 3.0e+00 1.4e-01 7.0e-03 5.09e-02
8.01 5.03e-03 2.9e-01 1.3e-03 6.3e-06 4.88e-03

検算. (a) 固有値 88 に最も近い他の固有値は 1010 と 66(どちらも距離 22)で、μ>8\mu > 8 のとき 2 番目に近いのは 1010 なので、予測の比は ∣8−μ∣/∣10−μ∣|8 - \mu|/|10 - \mu| です。(b) 角は構成に使った固有ベクトルとの間で測っています。(c) 観測した t3/t2t_3/t_2 は予測の比と 2 桁目まで合い、μ\mu を 88 に 10 倍近づけるごとに比も約 10 分の 1 になります。観測値が予測よりわずかに小さいのは、3 番目以降に近い固有値の成分がまだ残っていて、それがさらに速く縮むためです。

逆反復のシフトを固定せず、毎回いまの近似ベクトルの Rayleigh 商に取り替えると、収束の次数そのものが上がります。

定義 5.1(Rayleigh 商反復)↗

単位ベクトル x0\boldsymbol{x}_0 から、μk=RA(xk)\mu_k = R_A(\boldsymbol{x}_k) とおき、A−μkIA - \mu_k I が正則である限り

(A−μkI)yk+1=xk,xk+1=yk+1∥yk+1∥(A - \mu_k I)\boldsymbol{y}_{k+1} = \boldsymbol{x}_k,\qquad \boldsymbol{x}_{k+1} = \frac{\boldsymbol{y}_{k+1}}{\|\boldsymbol{y}_{k+1}\|}

を繰り返す反復を Rayleigh 商反復という。μk\mu_k が固有値になったら止める(固有値が求まっている)。

毎回シフトが変わるので、1 回ごとに A−μkIA - \mu_k I を分解し直す必要があります(密行列なら O(n3)O(n^3)、三重対角なら O(n)O(n))。その代わり、次の定理のとおり収束は 3 次です。

定理に出てくる δ\delta は、注目する固有値 λ\lambda から最も近い別の固有値までの距離、Δ\Delta は固有値全体の幅です。仮定 Δsin⁡2θ0<δ/2\Delta\sin^2\theta_0 < \delta/2 は、Rayleigh 商の誤差の上限(補題 3.4 の評価)が、他の固有値との距離の半分より小さい、つまり μk\mu_k が目的の固有値 λ\lambda に最も近いままである、ということを保証するためのものです。

定理 5.2(Rayleigh 商反復の局所 3 次収束)✓

AA を nn 次 Hermite 行列(n≥2n \ge 2)、λ\lambda をその単純な固有値、v\boldsymbol{v} を対応する単位固有ベクトルとし、

δ=min⁡{∣λ′−λ∣:λ′ は λ と異なる A の固有値},Δ=λmax⁡−λmin⁡\delta = \min\{|\lambda' - \lambda| : \lambda' \text{ は } \lambda \text{ と異なる } A \text{ の固有値}\},\qquad \Delta = \lambda_{\max} - \lambda_{\min}

とおく。Rayleigh 商反復の θk=∠(xk,v)\theta_k = \angle(\boldsymbol{x}_k, \boldsymbol{v}) が Δsin⁡2θ0<δ/2\Delta\sin^2\theta_0 < \delta/2 を満たすならば、反復が続く限り次が成り立つ。

  1. ∣μk−λ∣≤Δsin⁡2θk|\mu_k - \lambda| \le \Delta\sin^2\theta_k であり、μk\mu_k は λ\lambda 以外の固有値に一致しない。
  2. tan⁡θk+1≤Δsin⁡2θkδ−Δsin⁡2θktan⁡θk≤2Δδtan⁡3θk\displaystyle \tan\theta_{k+1} \le \frac{\Delta\sin^2\theta_k}{\delta - \Delta\sin^2\theta_k}\tan\theta_k \le \frac{2\Delta}{\delta}\tan^3\theta_k。とくに θk+1≤θk\theta_{k+1} \le \theta_k で、Δsin⁡2θk+1<δ/2\Delta\sin^2\theta_{k+1} < \delta/2 も保たれる。

したがって反復が有限回で止まらなければ tan⁡θk→0\tan\theta_k \to 0 であり、その収束は 3 次である。

証明(定理 5.2)

δ\delta は [λmin⁡,λmax⁡][\lambda_{\min}, \lambda_{\max}] の 2 点の距離なので δ≤Δ\delta \le \Delta で、仮定から sin⁡2θ0<1/2\sin^2\theta_0 < 1/2、とくに cos⁡θ0>0\cos\theta_0 > 0 です。kk について帰納的に、Δsin⁡2θk<δ/2\Delta\sin^2\theta_k < \delta/2 のもとで 1 と 2 を示します。

スペクトル定理(定理 4.2[スペクトル定理])で正規直交な固有ベクトル u1,…,un\boldsymbol{u}_1, \ldots, \boldsymbol{u}_n をとります。λ\lambda は単純なので、固有値が λ\lambda であるのはただ 1 つ uj\boldsymbol{u}_j だけで、v\boldsymbol{v} は uj\boldsymbol{u}_j の絶対値 1 の定数倍です(v\boldsymbol{v} を基底で展開すると、Av=λvA\boldsymbol{v} = \lambda\boldsymbol{v} から固有値が λ\lambda と異なる成分は 0)。よって v\boldsymbol{v} に直交するベクトル全体は {ui:i≠j}\{\boldsymbol{u}_i : i \ne j\} で張られます。

反復は xk\boldsymbol{x}_k を絶対値 1 の複素数倍しても、μk\mu_k も以後の方向も変わらないので、v∗xk=cos⁡θk>0\boldsymbol{v}^{*}\boldsymbol{x}_k = \cos\theta_k > 0 としてかまいません。sin⁡θk=0\sin\theta_k = 0 なら xk\boldsymbol{x}_k は固有ベクトルで μk=λ\mu_k = \lambda となり反復は止まるので、sin⁡θk>0\sin\theta_k > 0 とし、xk=cos⁡θk v+sin⁡θk w\boldsymbol{x}_k = \cos\theta_k\,\boldsymbol{v} + \sin\theta_k\,\boldsymbol{w}(∥w∥=1\|\boldsymbol{w}\| = 1、v∗w=0\boldsymbol{v}^{*}\boldsymbol{w} = 0)と書きます。

1 の証明. 補題 3.4 から ∣μk−λ∣≤Δsin⁡2θk<δ/2|\mu_k - \lambda| \le \Delta\sin^2\theta_k < \delta/2 です。λ\lambda と異なる固有値 λ′\lambda' については ∣λ′−μk∣≥∣λ′−λ∣−∣λ−μk∣>δ−δ/2>0|\lambda' - \mu_k| \ge |\lambda' - \lambda| - |\lambda - \mu_k| > \delta - \delta/2 > 0 なので、μk≠λ′\mu_k \ne \lambda' です。

2 の証明. 反復が続くので μk≠λ\mu_k \ne \lambda で、1 と合わせて A−μkIA - \mu_k I は正則です。w=∑i≠jaiui\boldsymbol{w} = \sum_{i\ne j}a_i\boldsymbol{u}_i と書くと

(A−μkI)−1xk=cos⁡θkλ−μk v+sin⁡θk∑i≠jaiλi−μk ui(A - \mu_k I)^{-1}\boldsymbol{x}_k = \frac{\cos\theta_k}{\lambda - \mu_k}\,\boldsymbol{v} + \sin\theta_k\sum_{i\ne j}\frac{a_i}{\lambda_i - \mu_k}\,\boldsymbol{u}_i

で、第 2 項は v\boldsymbol{v} に直交し、そのノルムは sin⁡θk ∥w∥/min⁡i≠j∣λi−μk∣≤sin⁡θk/(δ−∣λ−μk∣)\sin\theta_k\,\|\boldsymbol{w}\|/\min_{i\ne j}|\lambda_i - \mu_k| \le \sin\theta_k/(\delta - |\lambda - \mu_k|) 以下です(1 の証明と同じ三角不等式)。補題 2.2 の 4 と 3 から

tan⁡θk+1≤sin⁡θk/(δ−∣λ−μk∣)cos⁡θk/∣λ−μk∣=∣λ−μk∣δ−∣λ−μk∣tan⁡θk.\tan\theta_{k+1} \le \frac{\sin\theta_k/(\delta - |\lambda - \mu_k|)}{\cos\theta_k/|\lambda - \mu_k|} = \frac{|\lambda - \mu_k|}{\delta - |\lambda - \mu_k|}\tan\theta_k .

関数 t↦t/(δ−t)t \mapsto t/(\delta - t) は [0,δ)[0, \delta) で増加なので、∣λ−μk∣≤Δsin⁡2θk|\lambda - \mu_k| \le \Delta\sin^2\theta_k を代入して 2 の最初の不等式を得ます。Δsin⁡2θk<δ/2\Delta\sin^2\theta_k < \delta/2 から分母は δ/2\delta/2 より大きく、分子の係数は 11 未満なので tan⁡θk+1<tan⁡θk\tan\theta_{k+1} < \tan\theta_k、すなわち θk+1<θk\theta_{k+1} < \theta_k で、sin⁡2θk+1<sin⁡2θk\sin^2\theta_{k+1} < \sin^2\theta_k から帰納法の仮定が保たれます。最後に sin⁡2θk≤tan⁡2θk\sin^2\theta_k \le \tan^2\theta_k と分母 >δ/2> \delta/2 から (2Δ/δ)tan⁡3θk(2\Delta/\delta)\tan^3\theta_k の評価を得ます。

最後に収束を示します。ck=Δsin⁡2θk/(δ−Δsin⁡2θk)c_k = \Delta\sin^2\theta_k/(\delta - \Delta\sin^2\theta_k) とおくと、2 から tan⁡θk+1≤cktan⁡θk\tan\theta_{k+1} \le c_k\tan\theta_k で、θk\theta_k は減少するので ck≤c0<1c_k \le c_0 < 1 です(c0<1c_0 < 1 は Δsin⁡2θ0<δ/2\Delta\sin^2\theta_0 < \delta/2 と同値です)。よって tan⁡θk≤c0ktan⁡θ0→0\tan\theta_k \le c_0^k\tan\theta_0 \to 0 であり、各段で tan⁡θk+1≤(2Δ/δ)tan⁡3θk\tan\theta_{k+1} \le (2\Delta/\delta)\tan^3\theta_k が成り立つので収束は 3 次です。

正接が 1 回ごとに 3 乗になるので、正しい桁数は 1 回ごとにほぼ 3 倍になります。単根の近くの Newton 法が桁数を 1 回ごとにほぼ 2 倍にする(非線形方程式の反復解法)のと比べて、1 段高い次数です。固有値の誤差は 1 から sin⁡2θk\sin^2\theta_k の程度なので、こちらの桁数も 1 回ごとにほぼ 3 倍です。

次のコードは、対角が 22・副対角が −1-1 の 66 次の三重対角行列(固有値は 2−2cos⁡(jπ/7)2 - 2\cos(j\pi/7)、固有ベクトルは (sin⁡(ijπ/7))i(\sin(ij\pi/7))_i と閉じた形で分かる)で、Rayleigh 商反復を 100 桁で回したものです。100 桁(mp.mp.dps = 100、Python の多倍長計算ライブラリ mpmath の設定)を使うのは、桁数が 3 倍ずつ増えるので、倍精度(約 16 桁)では 3 回で限界に達して様子が見えなくなるためです。

# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0(x86-64)
import mpmath as mp
import numpy as np
mp.mp.dps = 100
n = 6
A = 2 * mp.eye(n) # 対角 2、副対角 -1 の三重対角行列
for i in range(n - 1):
A[i, i + 1] = A[i + 1, i] = -1
th = mp.pi / (n + 1)
lam = 2 - 2 * mp.cos(th) # 最小の固有値
v = mp.matrix([mp.sin((j + 1) * th) for j in range(n)])
v = v / mp.norm(v) # その固有ベクトル
x = mp.matrix([1] * n) / mp.sqrt(n)
print("k tan(角) |mu_k - lam|")
for k in range(5):
mu = (x.T * A * x)[0]
c = abs((v.T * x)[0])
t = mp.sqrt(1 - c**2) / c
print(k, mp.nstr(t, 3), mp.nstr(abs(mu - lam), 3))
y = mp.lu_solve(A - mu * mp.eye(n), x)
x = y / mp.norm(y)
An, xd = np.array(A.tolist(), dtype=float), np.ones(n) / np.sqrt(n)
print("倍精度の |mu_k - lam|")
for k in range(4):
mu = xd @ An @ xd
e = abs(mu - float(lam))
print(k, "< 1e-15" if e < 1e-15 else f"{e:.1e}")
yd = np.linalg.solve(An - mu * np.eye(n), xd)
xd = yd / np.linalg.norm(yd)
k tan(角) |mu_k - lam|
0 0.307 0.135
1 0.0321 0.00144
2 3.38e-5 1.56e-9
3 3.87e-14 2.03e-27
4 5.79e-41 4.56e-81
倍精度の |mu_k - lam|
0 1.4e-01
1 1.4e-03
2 1.6e-09
3 < 1e-15

検算. (a) 固有値と固有ベクトルは閉じた形の値と比べています。(b) 倍精度の ∣μk−λ∣|\mu_k - \lambda| は 100 桁の値と表示した桁で一致し、3 回目で倍精度の丸めの水準に達します。(c) 正接の指数は −1.5→−4.5→−13.4→−40.2-1.5 \to -4.5 \to -13.4 \to -40.2 とほぼ 3 倍ずつ、固有値の誤差の指数は −2.8→−8.8→−26.7→−80.3-2.8 \to -8.8 \to -26.7 \to -80.3 とやはりほぼ 3 倍ずつ進み、比 tan⁡θk+1/tan⁡3θk\tan\theta_{k+1}/\tan^3\theta_k は 1.11,1.02,1.00,1.001.11, 1.02, 1.00, 1.00 です。この行列では Δ=4cos⁡(π/7)≈3.60\Delta = 4\cos(\pi/7) \approx 3.60、δ=2cos⁡(π/7)−2cos⁡(2π/7)≈0.555\delta = 2\cos(\pi/7) - 2\cos(2\pi/7) \approx 0.555 で、定理の係数 2Δ/δ≈132\Delta/\delta \approx 13 は実際の比より大きく、評価として正しい側にあります。なお k=0k = 0 では Δsin⁡2θ0≈0.31>δ/2\Delta\sin^2\theta_0 \approx 0.31 > \delta/2 で定理の仮定を満たしませんが、k=1k = 1 以降は満たしています(定理は十分条件です)。

注意 5.3(一般の行列の Rayleigh 商反復)

Hermite でない行列では 補題 3.4 の 2 次の精度が成り立たず(系 3.3 の 1 次の評価しかない)、Rayleigh 商反復は単純な固有値の近くで局所 2 次収束にとどまります。左右の固有ベクトルを同時に反復する 2 側 Rayleigh 商反復にすると 3 次に戻ります。Ostrowski による一連の解析と、Parlett『The Symmetric Eigenvalue Problem』第 4 章の注を参照してください。

注意 5.4(Rayleigh 商反復の大域的な振る舞い)

定理 5.2 は初期ベクトルが固有ベクトルに十分近いときの局所的な主張です。Hermite 行列では、ほとんどすべての初期ベクトルから Rayleigh 商反復が何らかの固有ベクトルに収束すること(どの固有ベクトルに収束するかは予測しにくい)が Kahan と Parlett によって示されています。証明は残差のノルムの単調減少に基づく数頁の議論で、Parlett『The Symmetric Eigenvalue Problem』第 4 章にあります。本章では使いません。

注意 5.5(ほとんど特異な系を解いても壊れない理由)

逆反復や Rayleigh 商反復では、μ\mu が固有値に近いほど A−μIA - \mu I は特異に近く、(A−μI)y=x(A - \mu I)\boldsymbol{y} = \boldsymbol{x} の解の前進誤差は大きくなります。それでも算法が壊れないのは、後退安定な解法で求めた y^\hat{\boldsymbol{y}} が (A−μI+E)y^=x(A - \mu I + E)\hat{\boldsymbol{y}} = \boldsymbol{x}(∥E∥\|E\| は ∥A∥\|A\| の単位丸め倍程度)の厳密解であり、その誤差の大部分が、まさに求めたい固有ベクトルの方向にそろうからです。方向だけが問題の逆反復にとって、その方向の誤差は無害です。精密な議論は Parlett の前掲書第 4 章、Golub–Van Loan『Matrix Computations』§7.6 を参照してください。

1 本のベクトルの代わりに pp 本の正規直交なベクトルを同時に反復し、毎回直交化し直すと、上位 pp 個の固有値に対応する不変部分空間が求まります。直交化しないと、pp 本すべてが同じ v1\boldsymbol{v}_1 の方向に寄ってしまい、数値的に一次従属になります。

次の定理は記号が多いので、先に役割を整理します。XX は固有ベクトルを並べた行列で、その前半 X1X_1(上位 pp 個の固有ベクトル)と後半 X2X_2 に分けます。X−1X^{-1} の対応する行を W1∗,W2∗W_1^{*}, W_2^{*} と書くと、W1∗Z0W_1^{*}Z_0 は初期の pp 本が上位 pp 個の固有ベクトルの成分をもつかどうかを表す pp 次の行列で、これが正則であることが c1≠0c_1 \ne 0 の一般化です。証明の方針は冪乗法と同じで、AkZ0A^kZ_0 を上位の固有ベクトルの部分と残りの部分に分け、残りが比 ∣λp+1/λp∣k|\lambda_{p+1}/\lambda_p|^k で小さくなることを示します。

定義 6.1(同時反復)↗

Z0∈Cn×pZ_0 \in \mathbb{C}^{n\times p} を正規直交な列をもつ行列とする。Yk=AZk−1Y_k = AZ_{k-1} の列が一次独立である限り、その列に Gram–Schmidt の直交化(定理 6.1[内積空間とグラム・シュミット直交化])を施して

Yk=ZkRk(Zk は正規直交な列をもち、Rk は対角成分が正の上三角行列)Y_k = Z_kR_k\qquad(Z_k \text{ は正規直交な列をもち、} R_k \text{ は対角成分が正の上三角行列})

と書き、ZkZ_k を次の反復ベクトルとする。これを同時反復(部分空間反復、直交反復)という。p=1p = 1 のときは冪乗法である。

Gram–Schmidt の定理の 2(先頭の jj 本が張る空間が一致する)と 3(RkR_k の対角成分 zi∗yi\boldsymbol{z}_i^{*}\boldsymbol{y}_i が正)により、この分解は一意で、ZkZ_k の先頭 jj 列の張る空間は YkY_k の先頭 jj 列の張る空間に等しくなります。実際の計算では Householder 変換(鏡映を順に掛けて列を直交化する方法)による QR 分解を使いますが、張る空間は同じです。

定理 6.2(同時反復の収束)✓

A=XΛX−1A = X\Lambda X^{-1} を対角化可能な行列、Λ=diag⁡(λ1,…,λn)\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_n) とし、ある 1≤p<n1 \le p < n について

∣λ1∣≥⋯≥∣λp∣>∣λp+1∣≥⋯≥∣λn∣|\lambda_1| \ge \cdots \ge |\lambda_p| > |\lambda_{p+1}| \ge \cdots \ge |\lambda_n|

とする。X=(X1 X2)X = (X_1\ X_2)(X1X_1 は n×pn\times p)と分け、X−1X^{-1} を (W1∗W2∗)\begin{pmatrix} W_1^{*} \\ W_2^{*} \end{pmatrix}(W1W_1 は n×pn\times p)と分ける。Xp=span⁡X1\mathcal{X}_p = \operatorname{span}X_1(上位 pp 個の固有値の不変部分空間)とし、G=W1∗Z0G = W_1^{*}Z_0 が正則であると仮定する。このとき同時反復はすべての kk で定義され、span⁡Zk=span⁡AkZ0\operatorname{span}Z_k = \operatorname{span}A^kZ_0 であり、

r=∣λp+1λp∣,C=∥W1∥ ∥X2∥ ∥W2∗Z0G−1∥r = \left|\frac{\lambda_{p+1}}{\lambda_p}\right|,\qquad C = \|W_1\|\,\|X_2\|\,\|W_2^{*}Z_0G^{-1}\|

とおくと、Crk<1Cr^k < 1 のとき

d(span⁡Zk,Xp)≤Crk1−Crkd(\operatorname{span}Z_k, \mathcal{X}_p) \le \frac{Cr^k}{1 - Cr^k}

が成り立つ。

証明(定理 6.2)

準備. X−1X=IX^{-1}X = I を区分けすると W1∗X1=IpW_1^{*}X_1 = I_p、W1∗X2=OW_1^{*}X_2 = O で、X−1A=ΛX−1X^{-1}A = \Lambda X^{-1} を区分けすると W1∗A=Λ1W1∗W_1^{*}A = \Lambda_1W_1^{*} です。ここで Λ1=diag⁡(λ1,…,λp)\Lambda_1 = \operatorname{diag}(\lambda_1, \ldots, \lambda_p)、Λ2=diag⁡(λp+1,…,λn)\Lambda_2 = \operatorname{diag}(\lambda_{p+1}, \ldots, \lambda_n) とおきました。∣λp∣>∣λp+1∣≥0|\lambda_p| > |\lambda_{p+1}| \ge 0 から Λ1\Lambda_1 は正則です。

反復が定義されること. kk についての帰納法で、Zk=AkZ0MkZ_k = A^kZ_0M_k(MkM_k は正則な pp 次行列)と書けることを示します。k=0k = 0 は M0=IM_0 = I です。Zk−1=Ak−1Z0Mk−1Z_{k-1} = A^{k-1}Z_0M_{k-1} なら Yk=AkZ0Mk−1Y_k = A^kZ_0M_{k-1} で、

W1∗Yk=W1∗AkZ0Mk−1=Λ1kW1∗Z0Mk−1=Λ1kGMk−1W_1^{*}Y_k = W_1^{*}A^kZ_0M_{k-1} = \Lambda_1^kW_1^{*}Z_0M_{k-1} = \Lambda_1^kGM_{k-1}

は正則です。Ykc=0Y_k\boldsymbol{c} = \boldsymbol{0} なら W1∗Ykc=0W_1^{*}Y_k\boldsymbol{c} = \boldsymbol{0} から c=0\boldsymbol{c} = \boldsymbol{0} なので、YkY_k の列は一次独立で、Gram–Schmidt の直交化が実行できます。Yk=ZkRkY_k = Z_kR_k の RkR_k は正則なので Zk=AkZ0Mk−1Rk−1Z_k = A^kZ_0M_{k-1}R_k^{-1} で、Mk=Mk−1Rk−1M_k = M_{k-1}R_k^{-1} とすればよい。とくに span⁡Zk=span⁡AkZ0\operatorname{span}Z_k = \operatorname{span}A^kZ_0 です。

隔たりの評価. Ak=XΛkX−1A^k = X\Lambda^kX^{-1} より

AkZ0=X1Λ1kW1∗Z0+X2Λ2kW2∗Z0A^kZ_0 = X_1\Lambda_1^kW_1^{*}Z_0 + X_2\Lambda_2^kW_2^{*}Z_0

で、右から正則な G−1Λ1−kG^{-1}\Lambda_1^{-k} を掛けると

AkZ0G−1Λ1−k=X1+Ek,Ek=X2Λ2k (W2∗Z0G−1) Λ1−kA^kZ_0G^{-1}\Lambda_1^{-k} = X_1 + E_k,\qquad E_k = X_2\Lambda_2^k\,(W_2^{*}Z_0G^{-1})\,\Lambda_1^{-k}

です(G−1Λ1−kG^{-1}\Lambda_1^{-k} は Λ1kG\Lambda_1^kG の逆行列で、X1Λ1kG G−1Λ1−k=X1X_1\Lambda_1^kG\,G^{-1}\Lambda_1^{-k} = X_1)。対角行列のノルムは対角成分の絶対値の最大なので ∥Λ2k∥=∣λp+1∣k\|\Lambda_2^k\| = |\lambda_{p+1}|^k、∥Λ1−k∥=∣λp∣−k\|\Lambda_1^{-k}\| = |\lambda_p|^{-k} で、

∥Ek∥≤∥X2∥ ∣λp+1∣k ∥W2∗Z0G−1∥ ∣λp∣−k,∥W1∥ ∥Ek∥≤Crk\|E_k\| \le \|X_2\|\,|\lambda_{p+1}|^k\,\|W_2^{*}Z_0G^{-1}\|\,|\lambda_p|^{-k}, \qquad \|W_1\|\,\|E_k\| \le Cr^k

です。span⁡Zk=span⁡(X1+Ek)\operatorname{span}Z_k = \operatorname{span}(X_1 + E_k) なので、その零でない元は y=(X1+Ek)c\boldsymbol{y} = (X_1 + E_k)\boldsymbol{c}(c≠0\boldsymbol{c} \ne \boldsymbol{0})と書けます。W1∗X1=IpW_1^{*}X_1 = I_p から ∥c∥=∥W1∗X1c∥≤∥W1∥∥X1c∥\|\boldsymbol{c}\| = \|W_1^{*}X_1\boldsymbol{c}\| \le \|W_1\|\|X_1\boldsymbol{c}\| なので(ここで ∥W1∗∥≤∥W1∥\|W_1^{*}\| \le \|W_1\| を使いました。任意の u\boldsymbol{u} について Cauchy–Schwarz の不等式から ∥W1∗u∥2=u∗W1(W1∗u)≤∥u∥ ∥W1∥ ∥W1∗u∥\|W_1^{*}\boldsymbol{u}\|^2 = \boldsymbol{u}^{*}W_1(W_1^{*}\boldsymbol{u}) \le \|\boldsymbol{u}\|\,\|W_1\|\,\|W_1^{*}\boldsymbol{u}\| で、W1∗u≠0W_1^{*}\boldsymbol{u} \ne \boldsymbol{0} なら両辺を ∥W1∗u∥\|W_1^{*}\boldsymbol{u}\| で割れば従います)

∥y∥≥∥X1c∥−∥Ekc∥≥(1∥W1∥−∥Ek∥)∥c∥>0\|\boldsymbol{y}\| \ge \|X_1\boldsymbol{c}\| - \|E_k\boldsymbol{c}\| \ge \Bigl(\frac{1}{\|W_1\|} - \|E_k\|\Bigr)\|\boldsymbol{c}\| > 0

です(∥W1∥∥Ek∥≤Crk<1\|W_1\|\|E_k\| \le Cr^k < 1)。X1c∈XpX_1\boldsymbol{c} \in \mathcal{X}_p なので 補題 2.2 の 2 から

sin⁡∠(y,Xp)≤∥y−X1c∥∥y∥=∥Ekc∥∥y∥≤∥Ek∥1/∥W1∥−∥Ek∥=∥W1∥∥Ek∥1−∥W1∥∥Ek∥≤Crk1−Crk\sin\angle(\boldsymbol{y}, \mathcal{X}_p) \le \frac{\|\boldsymbol{y} - X_1\boldsymbol{c}\|}{\|\boldsymbol{y}\|} = \frac{\|E_k\boldsymbol{c}\|}{\|\boldsymbol{y}\|} \le \frac{\|E_k\|}{1/\|W_1\| - \|E_k\|} = \frac{\|W_1\|\|E_k\|}{1 - \|W_1\|\|E_k\|} \le \frac{Cr^k}{1 - Cr^k}

で、y\boldsymbol{y} について上限をとれば主張を得ます(t↦t/(1−t)t \mapsto t/(1 - t) は増加)。

p=1p = 1 とすると、W1∗W_1^{*} は左固有ベクトル(X−1X^{-1} の第 1 行)で、G≠0G \ne 0 は「初期ベクトルが v1\boldsymbol{v}_1 の成分をもつ」という 定理 3.2 の条件 c1≠0c_1 \ne 0 と同じです。

同時反復の要は、Gram–Schmidt の直交化が先頭の列から順に張る空間を保つことです。このため、先頭の jj 列だけを見ると、それ自体が jj 本の同時反復になっています。

系 6.3(先頭の列の収束)✓

A=XΛX−1A = X\Lambda X^{-1} を 定理 6.2 と同じく対角化可能とし、固有値は ∣λ1∣≥⋯≥∣λn∣|\lambda_1| \ge \cdots \ge |\lambda_n| と並んでいるとする。Z0∈Cn×pZ_0 \in \mathbb{C}^{n\times p}(1≤p≤n1 \le p \le n)は正規直交な列をもち、同時反復がすべての kk で定義されるとする(たとえば AA が正則なら AZk−1AZ_{k-1} の列はつねに一次独立)。1≤j≤p1 \le j \le p、j<nj < n について ∣λj∣>∣λj+1∣|\lambda_j| > |\lambda_{j+1}| を仮定し、Z0Z_0 の先頭 jj 列を Z0(j)Z_0^{(j)}、X−1X^{-1} の先頭 jj 行を (W1(j))∗(W_1^{(j)})^{*} と書いて、(W1(j))∗Z0(j)(W_1^{(j)})^{*}Z_0^{(j)} が正則とする。このとき ZkZ_k の先頭 jj 列 Zk(j)Z_k^{(j)} は、Z0(j)Z_0^{(j)} から始めた同時反復の kk 回目と一致し、その張る空間は Xj=span⁡{x1,…,xj}\mathcal{X}_j = \operatorname{span}\{\boldsymbol{x}_1, \ldots, \boldsymbol{x}_j\}(XX の先頭 jj 列)に O(∣λj+1/λj∣k)O(|\lambda_{j+1}/\lambda_j|^k) で近づく。

証明(系 6.3)

kk についての帰納法で、Zk(j)Z_k^{(j)} が Z0(j)Z_0^{(j)} から始めた同時反復の kk 回目 Z~k\tilde{Z}_k に一致することを示します。Zk−1(j)=Z~k−1Z_{k-1}^{(j)} = \tilde{Z}_{k-1} なら、Yk=AZk−1Y_k = AZ_{k-1} の先頭 jj 列は Y~k=AZ~k−1\tilde{Y}_k = A\tilde{Z}_{k-1} です。Gram–Schmidt の直交化の漸化式(定理 6.1[内積空間とグラム・シュミット直交化])は第 ii 列を作るのに第 ii 列までしか使わないので、YkY_k を直交化した結果の先頭 jj 列は Y~k\tilde{Y}_k を直交化した結果に一致し、Zk(j)=Z~kZ_k^{(j)} = \tilde{Z}_k です(Y~k\tilde{Y}_k は YkY_k の先頭 jj 列なので、その列は一次独立です)。あとは pp を jj に置き換えた 定理 6.2 を Z~k\tilde{Z}_k に使えば、隔たりは O(∣λj+1/λj∣k)O(|\lambda_{j+1}/\lambda_j|^k) です。

したがって ∣λ1∣>∣λ2∣>⋯>∣λp∣>∣λp+1∣|\lambda_1| > |\lambda_2| > \cdots > |\lambda_p| > |\lambda_{p+1}| なら、ZkZ_k の第 1 列は x1\boldsymbol{x}_1 の方向に、先頭 2 列は span⁡{x1,x2}\operatorname{span}\{\boldsymbol{x}_1, \boldsymbol{x}_2\} に、… と、入れ子になった不変部分空間の列がそれぞれの比の速さで同時に求まります。p=np = n としたものを行列の相似変換の言葉で書き直したのが QR 法で、その収束の解析は QR 法の章でこの系をもとに行います。

反復を止めたとき、手元にあるのは近似固有対 (μ,x)(\mu, \boldsymbol{x}) です。真の固有ベクトルは分からないので、角 θ\theta そのものは測れません。測れるのは残差です。

命題 7.1(Rayleigh 商は残差を最小にする)✓

A∈Cn×nA \in \mathbb{C}^{n\times n}、∥x∥=1\|\boldsymbol{x}\| = 1 とし、r=Ax−RA(x)x\boldsymbol{r} = A\boldsymbol{x} - R_A(\boldsymbol{x})\boldsymbol{x} とおく。このとき x∗r=0\boldsymbol{x}^{*}\boldsymbol{r} = 0、∥r∥2=∥Ax∥2−∣RA(x)∣2\|\boldsymbol{r}\|^2 = \|A\boldsymbol{x}\|^2 - |R_A(\boldsymbol{x})|^2 であり、任意の μ∈C\mu \in \mathbb{C} について

∥Ax−μx∥2=∥r∥2+∣RA(x)−μ∣2.\|A\boldsymbol{x} - \mu\boldsymbol{x}\|^2 = \|\boldsymbol{r}\|^2 + |R_A(\boldsymbol{x}) - \mu|^2 .

とくに μ↦∥Ax−μx∥\mu \mapsto \|A\boldsymbol{x} - \mu\boldsymbol{x}\| は μ=RA(x)\mu = R_A(\boldsymbol{x}) でだけ最小になる。

証明(命題 7.1)

R=RA(x)=x∗AxR = R_A(\boldsymbol{x}) = \boldsymbol{x}^{*}A\boldsymbol{x} とおくと x∗r=x∗Ax−R x∗x=0\boldsymbol{x}^{*}\boldsymbol{r} = \boldsymbol{x}^{*}A\boldsymbol{x} - R\,\boldsymbol{x}^{*}\boldsymbol{x} = 0 です。Ax−μx=r+(R−μ)xA\boldsymbol{x} - \mu\boldsymbol{x} = \boldsymbol{r} + (R - \mu)\boldsymbol{x} の 2 項は直交するので、ノルムの 2 乗は ∥r∥2+∣R−μ∣2\|\boldsymbol{r}\|^2 + |R - \mu|^2 です。μ=0\mu = 0 とおけば ∥Ax∥2=∥r∥2+∣R∣2\|A\boldsymbol{x}\|^2 = \|\boldsymbol{r}\|^2 + |R|^2 を得ます。

系 7.2(実対称行列の残差評価)✓

AA を実対称行列、x∈Rn\boldsymbol{x} \in \mathbb{R}^n を単位ベクトル、r=Ax−RA(x)x\boldsymbol{r} = A\boldsymbol{x} - R_A(\boldsymbol{x})\boldsymbol{x} とする。このとき ∣λ−RA(x)∣≤∥r∥|\lambda - R_A(\boldsymbol{x})| \le \|\boldsymbol{r}\| を満たす AA の固有値 λ\lambda がある。∥r∥\|\boldsymbol{r}\| は、すべての実数 μ\mu についての残差 ∥Ax−μx∥\|A\boldsymbol{x} - \mu\boldsymbol{x}\| の中で最小である。

証明(系 7.2)

AA は実対称なので RA(x)R_A(\boldsymbol{x}) は実数で、残差による固有値の評価(定理 6.1[行列計算の認証と最適化の証明書])を μ=RA(x)\mu = R_A(\boldsymbol{x}) に使えば前半です。後半は 命題 7.1 です。

一般の行列では、残差が小さくても固有値が近いとは限りません。差を測るのは固有ベクトル行列の条件数です。

命題 7.3(対角化可能な行列の残差評価)✓

A=VΛV−1A = V\Lambda V^{-1} を対角化可能な行列とし、κ(V)=∥V∥∥V−1∥\kappa(V) = \|V\|\|V^{-1}\| とおく。∥x∥=1\|\boldsymbol{x}\| = 1、μ∈C\mu \in \mathbb{C}、r=Ax−μx\boldsymbol{r} = A\boldsymbol{x} - \mu\boldsymbol{x} ならば、min⁡i∣λi−μ∣≤κ(V)∥r∥\min_i|\lambda_i - \mu| \le \kappa(V)\|\boldsymbol{r}\| である。

証明(命題 7.3)

E=−rx∗E = -\boldsymbol{r}\boldsymbol{x}^{*} とおくと (A+E)x=Ax−r(x∗x)=μx(A + E)\boldsymbol{x} = A\boldsymbol{x} - \boldsymbol{r}(\boldsymbol{x}^{*}\boldsymbol{x}) = \mu\boldsymbol{x} なので、μ\mu は A+EA + E の固有値です。任意の y\boldsymbol{y} について Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∥Ey∥=∣x∗y∣∥r∥≤∥r∥∥y∥\|E\boldsymbol{y}\| = |\boldsymbol{x}^{*}\boldsymbol{y}|\|\boldsymbol{r}\| \le \|\boldsymbol{r}\|\|\boldsymbol{y}\| なので ∥E∥≤∥r∥\|E\| \le \|\boldsymbol{r}\| です。Bauer–Fike の定理(定理 6.2[行列計算の認証と最適化の証明書])から min⁡i∣μ−λi∣≤κ(V)∥E∥≤κ(V)∥r∥\min_i|\mu - \lambda_i| \le \kappa(V)\|E\| \le \kappa(V)\|\boldsymbol{r}\| です。

Hermite 行列では VV をユニタリにとれて κ(V)=1\kappa(V) = 1 なので、命題 7.3 は 系 7.2 の複素版になっています。∥r∥\|\boldsymbol{r}\| の評価は 1 次ですが、固有値の隔たりについての情報が加われば 2 次に改善できます。その代表が Temple の不等式です。

例 7.4(Rayleigh 商反復 1 回のあとの包み込み)

§5 のコードの 66 次の三重対角行列で、x0=(1,…,1)T/6\boldsymbol{x}_0 = (1, \ldots, 1)^{\mathsf{T}}/\sqrt{6} から Rayleigh 商反復を 1 回行うと(mpmath 30 桁で計算)、η=RA(x1)≈0.199503\eta = R_A(\boldsymbol{x}_1) \approx 0.199503、∥r∥≈0.0457\|\boldsymbol{r}\| \approx 0.0457 です。

  • 系 7.2 からは、区間 [η−∥r∥,η+∥r∥]≈[0.154,0.245][\eta - \|\boldsymbol{r}\|, \eta + \|\boldsymbol{r}\|] \approx [0.154, 0.245] に固有値があることが分かります。
  • 最小の固有値を λ1\lambda_1、2 番目に小さい固有値を λ2\lambda_2 とします。λ2≥0.75\lambda_2 \ge 0.75 が分かっていれば(ここでは閉じた形 2−2cos⁡(2π/7)≈0.753022 - 2\cos(2\pi/7) \approx 0.75302 から。実際の計算では Sturm 列の数え上げなどで確かめます)、Temple の不等式(定理 4.1[無限次元の固有値の包含])を ν=0.75\nu = 0.75 で使えます。命題 7.1 から Temple の不等式の σ2=∥Ax∥2−η2\sigma^2 = \|A\boldsymbol{x}\|^2 - \eta^2 は ∥r∥2≈0.00209\|\boldsymbol{r}\|^2 \approx 0.00209 に等しく、λ1≥η−∥r∥2/(ν−η)≈0.195712\lambda_1 \ge \eta - \|\boldsymbol{r}\|^2/(\nu - \eta) \approx 0.195712 です。
  • 上からは Rayleigh 商そのものが λ1≤η\lambda_1 \le \eta を与えます(定理 7.1[スペクトル定理])。

こうして λ1∈[0.19571,0.19951]\lambda_1 \in [0.19571, 0.19951] と、残差だけの評価より約 24 倍狭い区間に包めます。真の値は 2−2cos⁡(π/7)≈0.1980622 - 2\cos(\pi/7) \approx 0.198062 です。Temple の評価の幅 ∥r∥2/(ν−η)\|\boldsymbol{r}\|^2/(\nu - \eta) は残差の 2 乗に比例するので、反復をもう 1 回進めれば幅は桁違いに縮みます。

注意 7.5(2 番目以降の固有値:デフレーション)

冪乗法で (λ1,v1)(\lambda_1, \boldsymbol{v}_1) が求まったあと次の固有値を求めるには、λ1\lambda_1 を取り除いた行列に冪乗法を使います。Hermite 行列なら A−λ1v1v1∗A - \lambda_1\boldsymbol{v}_1\boldsymbol{v}_1^{*}(Hotelling のデフレーション)の固有値は 0,λ2,…,λn0, \lambda_2, \ldots, \lambda_n です。一般の行列では A−v1w∗A - \boldsymbol{v}_1\boldsymbol{w}^{*}(w∗v1=λ1\boldsymbol{w}^{*}\boldsymbol{v}_1 = \lambda_1)とする Wielandt のデフレーションがあります。ただし v1\boldsymbol{v}_1 の誤差が次の固有値に持ち越されるので、多くの固有値が要るときは本章の同時反復や、Krylov 部分空間法の章の Lanczos 法を使うのが普通です。Saad『Numerical Methods for Large Eigenvalue Problems』第 4 章を参照してください。

演習 8.1易

実対称行列で、絶対値最大の固有値と 2 番目の固有値の比が ∣λ2/λ1∣=0.9|\lambda_2/\lambda_1| = 0.9、初期ベクトルの角が tan⁡∠(x0,u1)=1\tan\angle(\boldsymbol{x}_0, \boldsymbol{u}_1) = 1 であるとする。冪乗法の Rayleigh 商の誤差を (λmax⁡−λmin⁡)×10−10(\lambda_{\max} - \lambda_{\min}) \times 10^{-10} 以下にするには、系 3.5 の評価で何回の反復が十分か。固有ベクトルの角の正接を 10−1010^{-10} 以下にする回数とも比べよ。

解答

系 3.5 より ∣μk−λ1∣≤(λmax⁡−λmin⁡) 0.81k|\mu_k - \lambda_1| \le (\lambda_{\max} - \lambda_{\min})\,0.81^k なので、0.81k≤10−100.81^k \le 10^{-10}、すなわち k≥10/log⁡10(1/0.81)≈10/0.0915≈109.3k \ge 10/\log_{10}(1/0.81) \approx 10/0.0915 \approx 109.3 から k=110k = 110 回で十分です。角の正接は tan⁡∠(xk,u1)≤0.9k\tan\angle(\boldsymbol{x}_k, \boldsymbol{u}_1) \le 0.9^k なので k≥10/log⁡10(1/0.9)≈218.5k \ge 10/\log_{10}(1/0.9) \approx 218.5、つまり 219 回が必要です。固有値は固有ベクトルの約半分の回数で同じ桁数に達します。

演習 8.2標準

実対称行列の固有値が 1,2,…,101, 2, \ldots, 10 であるとする。A−σIA - \sigma I(σ\sigma は実数)に冪乗法を使って固有値 1010 を求めるとき、収束率 max⁡i≤9∣i−σ∣/∣10−σ∣\max_{i\le9}|i - \sigma|/|10 - \sigma| を最小にする σ\sigma と、そのときの収束率を求めよ。シフトなし(σ=0\sigma = 0)の場合と比べて、同じ精度に要する反復回数は何倍になるか。

解答

σ<10\sigma < 10 の範囲で考えます(σ\sigma が 1010 に近すぎると固有値 11 の側より 1010 の側が小さくなり、絶対値最大でなくなる)。i≤9i \le 9 について ∣i−σ∣|i - \sigma| の最大は max⁡(∣9−σ∣,∣1−σ∣)\max(|9 - \sigma|, |1 - \sigma|) で、これは σ=5\sigma = 5 で最小値 44 をとります。収束率 f(σ)=max⁡(∣9−σ∣,∣1−σ∣)/(10−σ)f(\sigma) = \max(|9 - \sigma|, |1 - \sigma|)/(10 - \sigma) は、σ≤5\sigma \le 5 では (9−σ)/(10−σ)=1−1/(10−σ)(9 - \sigma)/(10 - \sigma) = 1 - 1/(10 - \sigma) で σ\sigma について減少、5≤σ<105 \le \sigma < 10 では (σ−1)/(10−σ)(\sigma - 1)/(10 - \sigma) で増加なので、最小は σ=5\sigma = 5 の 4/5=0.84/5 = 0.8 です。シフトなしでは 0.90.9 なので、反復回数の比は log⁡0.9/log⁡0.8≈0.472\log 0.9/\log 0.8 \approx 0.472、つまり約半分で済みます。

演習 8.3標準

A=diag⁡(a,b)A = \operatorname{diag}(a, b)(a≠ba \ne b は実数)と単位ベクトル x=(cos⁡θ,sin⁡θ)T\boldsymbol{x} = (\cos\theta, \sin\theta)^{\mathsf{T}}(0<θ<π/20 < \theta < \pi/2)について、Rayleigh 商反復を 1 回行ったあとのベクトルと e1\boldsymbol{e}_1 のなす角 θ′\theta' が tan⁡θ′=tan⁡3θ\tan\theta' = \tan^3\theta を満たすことを示せ(RA(x)R_A(\boldsymbol{x}) が a,ba, b に一致しない場合)。定理 5.2 の評価の係数 2Δ/δ2\Delta/\delta と比べると、この場合の係数は何か。

解答

μ=acos⁡2θ+bsin⁡2θ\mu = a\cos^2\theta + b\sin^2\theta なので a−μ=(a−b)sin⁡2θa - \mu = (a - b)\sin^2\theta、b−μ=−(a−b)cos⁡2θb - \mu = -(a - b)\cos^2\theta です。よって

(A−μI)−1x=1a−b(cos⁡θsin⁡2θ, −sin⁡θcos⁡2θ)T(A - \mu I)^{-1}\boldsymbol{x} = \frac{1}{a - b}\left(\frac{\cos\theta}{\sin^2\theta},\ -\frac{\sin\theta}{\cos^2\theta}\right)^{\mathsf{T}}

で、補題 2.2 の 4 から tan⁡θ′=sin⁡θ/cos⁡2θcos⁡θ/sin⁡2θ=sin⁡3θcos⁡3θ=tan⁡3θ\tan\theta' = \dfrac{\sin\theta/\cos^2\theta}{\cos\theta/\sin^2\theta} = \dfrac{\sin^3\theta}{\cos^3\theta} = \tan^3\theta です。係数は 11 で、定理の係数 2Δ/δ=22\Delta/\delta = 2(22 次では Δ=δ=∣a−b∣\Delta = \delta = |a - b|)より小さく、定理の評価と矛盾しません。2 次の場合は誤差がちょうど 3 乗になり、正しい桁数がきっかり 3 倍になります。

  • G. H. Golub, C. F. Van Loan『Matrix Computations』4th ed., Johns Hopkins University Press, 2013 — §7.3(冪乗法・同時反復)、§7.6(逆反復)、§8.2(対称行列の冪乗法と Rayleigh 商反復)。
  • B. N. Parlett『The Symmetric Eigenvalue Problem』SIAM, 1998 — 第 4 章(冪乗法・逆反復・Rayleigh 商反復)、第 14 章(部分空間反復)。
  • L. N. Trefethen, D. Bau III『Numerical Linear Algebra』SIAM, 1997 — Lecture 27(Rayleigh 商・逆反復)、Lecture 28(同時反復から QR 法へ)。
  • J. W. Demmel『Applied Numerical Linear Algebra』SIAM, 1997 — 第 4 章(非対称固有値問題の反復)、第 5 章(対称固有値問題)。
  • Y. Saad『Numerical Methods for Large Eigenvalue Problems』revised ed., SIAM, 2011 — 第 4 章・第 5 章。

このテーマの 100 本ノック(100 問)

執筆: Claude Code(Opus 5.5)、Claude Code(Sonnet 5.5)・監査 2 回(Claude Code)・最終 2026-10-12

この記事の誤りを報告する ・運営: 夢現技研合同会社 ・料金 ・利用条件 ・特定商取引法に基づく表記

© 2026 夢現技研合同会社 ・コード例は MIT ライセンスです。