行列の固有値は、行列が「ベクトルをどれだけ伸ばす向きがあるか」を表す数です。小さな行列なら方程式を解いて求められますが、大きな行列では、適当なベクトルに行列を何度も掛けるだけで、いちばん強く伸びる向きが浮かび上がってきます。この記事では、この単純な繰り返しがなぜ固有値と固有ベクトルに行き着くのか、どれくらいの速さで近づくのか、そして掛ける前に行列を少しずらす工夫でどれだけ速くなるのかを、数値例とともに確かめます。
- 冪乗法は xk+1=Axk/∥Axk∥ を繰り返すだけの算法です。対角化可能で ∣λ1∣>∣λ2∣ なら、xk と固有ベクトル v1 のなす角は O(∣λ2/λ1∣k) で 0 に近づきます。
- Hermite 行列では、Rayleigh 商の誤差は角の正弦の 2 乗に比例します(RA(x)−λ=sin2θ(RA(w)−λ))。このため固有値は ∣λ2/λ1∣2k の速さで求まります。
- シフトつき逆反復は (A−μI)−1 に対する冪乗法で、収束率は ∣λJ−μ∣/∣λK−μ∣(λJ は μ に最も近い固有値、λK は 2 番目に近い固有値)です。シフトを固有値に近づけるほど速くなります。
- シフトを毎回 Rayleigh 商に取り替える Rayleigh 商反復は、Hermite 行列では角の正接が tanθk+1≤(2Δ/δ)tan3θk を満たし、局所 3 次収束(誤差が 1 回ごとにおよそ 3 乗になる速さの収束)します。
- p 本のベクトルを同時に掛けて直交化する同時反復は、優勢な p 次元の不変部分空間(行列を掛けても外へ出ない部分空間)に O(∣λp+1/λp∣k) で近づきます。これが次の章の QR 法の収束の理論の土台になります。
- 計算した近似固有対は、残差 r=Ax−RA(x)x を使って事後に評価できます(残差による評価・Temple の不等式・Bauer–Fike の定理。いずれも、残差の大きさから固有値の誤差を押さえる不等式です)。
研究で使う固有値の計算のほとんどは、numpy.linalg.eigh や scipy.sparse.linalg.eigsh を呼ぶだけで済みます。それでも「何回反復すれば何桁合うのか」「固有値が近いと何が起きるのか」「求めた固有値はどこまで信用できるのか」を答えるには、中で動いている反復の理論が要ります。
固有値問題は多項式の根を求める問題を含むので、5 次以上では有限回の四則と冪根で解く公式がありません。したがって固有値の算法は必ず反復法です。その最も単純な形が冪乗法で、QR 法・Lanczos 法・Arnoldi 法といった実用の算法は、どれも冪乗法を賢く組み合わせたものとして理解できます。
本章では、固有ベクトルを含む方向へベクトルを寄せていく反復を、次の 4 段で扱います。
- 冪乗法:Akx0 の方向は、絶対値最大の固有値の固有ベクトルへ向かう。
- 逆反復:(A−μI)−1 に冪乗法を使えば、μ に最も近い固有値が「最大」になる。
- Rayleigh 商反復:シフト μ を毎回いまの近似固有値に更新すると、収束は 3 次になる。
- 同時反復:1 本のベクトルではなく部分空間を反復すると、上位 p 個の固有値に対応する不変部分空間が求まる。
証明はどれも「固有ベクトルの基底で展開し、各成分が固有値の冪で伸び縮みする様子を見る」という同じ考え方です。
この章では Cn の標準内積 y∗x とそのノルム ∥x∥=(x∗x)1/2 だけを使います。行列 B∈Cm×n のノルムは、このノルムから誘導される作用素ノルム
∥B∥=x=0sup∥x∥∥Bx∥
です。B の第 i 行を bi∗ とすると、Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∥Bx∥2=∑i∣bi∗x∣2≤(∑i∥bi∥2)∥x∥2 なので、この上限は有限です。定義から直ちに ∥Bx∥≤∥B∥∥x∥、∥BC∥≤∥B∥∥C∥ が従います。対角行列 D=diag(d1,…,dn) については ∥Dx∥2=∑i∣di∣2∣xi∣2≤maxi∣di∣2∥x∥2 で、最大の ∣di∣ を与える座標ベクトルで等号が成り立つので ∥D∥=maxi∣di∣ です。
以下の定義は、反復が「固有ベクトルの方向(1 本のときは直線、同時反復のときは部分空間)にどれだけ近づいたか」を測るための物差しです。1 本のベクトルには角を、部分空間どうしには、片方のどのベクトルももう一方に近いかどうかを表す隔たりを使います。
定義 2.1(ベクトルどうし・ベクトルと部分空間のなす角)↗
x,v∈Cn を零でないベクトルとする。v への直交射影を Pvx=v∗vv∗xv とし、
cos∠(x,v)=∥x∥∥v∥∣v∗x∣,sin∠(x,v)=∥x∥∥x−Pvx∥と定める。cos∠=0 のとき tan∠=sin∠/cos∠ とする。角 ∠(x,v)∈[0,π/2] はこの正弦・余弦をもつ角である。
部分空間 S={0} への直交射影を PS とするとき(定理 7.1[内積空間とグラム・シュミット直交化] により存在する)、sin∠(x,S)=∥x−PSx∥/∥x∥ とおき、部分空間 U から S への隔たりを
d(U,S)=x∈U, x=0supsin∠(x,S)で定める。
d(U,S) が小さいとは、「U のどの単位ベクトルも S のすぐ近くにある」ということです。
補題 2.2(角の基本性質)✓
x,v を零でないベクトルとする。
- cos2∠(x,v)+sin2∠(x,v)=1。
- sin∠(x,v)=minα∈C∥x−αv∥/∥x∥。部分空間 S についても sin∠(x,S)=mins∈S∥x−s∥/∥x∥。
- 零でない複素数 a,b について ∠(ax,bv)=∠(x,v)。
- x=av+z、v∗z=0、a=0 と書けるとき tan∠(x,v)=∥z∥/(∣a∣∥v∥)。
証明(補題 2.2)
p=Pvx、z=x−p とおくと v∗z=v∗x−v∗x=0 です。
- p は v の定数倍なので p∗z=0 で、∥x∥2=∥p∥2+∥z∥2 です。∥p∥=∣v∗x∣/∥v∥ なので、両辺を ∥x∥2 で割れば主張になります。
- 任意の α について x−αv=z+(p−αv) で、第 2 項は v の定数倍なので z と直交します。よって ∥x−αv∥2=∥z∥2+∥p−αv∥2≥∥z∥2 で、等号は αv=p で成り立ちます。部分空間の場合も、x−PSx が S に直交すること(定理 7.1[内積空間とグラム・シュミット直交化])から同じ計算で示せます。
- 余弦の定義の分子・分母がともに ∣a∣∣b∣ 倍になるので余弦は変わらず、1 により正弦も変わりません。
- 分解 x=p+z は一意です(v の倍数と v に直交するベクトルへの分解)。よって p=av で、tan∠=∥z∥/∥p∥=∥z∥/(∣a∣∥v∥) です。
∎
固有ベクトル v の方向と、それに直交する方向 w への分解。x = cos θ · v + sin θ · w と書くと、本章の評価はすべて θ の漸化式になる。
定義 2.3(Rayleigh 商と残差)↗
A∈Cn×n、x=0 について
RA(x)=x∗xx∗Axを x の Rayleigh 商、r=Ax−RA(x)x を残差という。
A が Hermite 行列なら RA(x) は実数で、最小の固有値と最大の固有値の間にあります(定理 7.1[スペクトル定理]。そこでの R(x) と同じものです)。
定義 3.1(冪乗法)↗
A∈Cn×n と x0=0 に対し、Axk=0 である限り
xk+1=∥Axk∥Axk,μk+1=RA(xk+1)で xk と固有値の近似 μk を作る反復を冪乗法という。
帰納法により xk は Akx0 の正の定数倍です。正規化は桁あふれを防ぐためだけのもので、方向は Akx0 と同じです。1 回の反復の主な手間は行列とベクトルの積 1 回なので、疎行列なら非零成分の個数に比例する手間で済みます。
下の動画は 2 次の対称行列(固有ベクトル v1 の向きは共通)で xk の向きが v1 に寄っていく様子です。比 ∣λ2/λ1∣ が大きい左は遅く、小さい右は速く収束します。
冪乗法で単位ベクトル x_k(赤い矢印)が固有ベクトル v₁ の方向(緑の破線)へ近づく様子。左は比 0.8、右は比 0.4。薄い点はこれまでの x_j。
定理 3.2(冪乗法の収束)✓
A∈Cn×n は対角化可能で、固有値を重複を込めて ∣λ1∣>∣λ2∣≥∣λ3∣≥⋯≥∣λn∣ と並べられるとする。対応する固有ベクトル v1,…,vn(∥vi∥=1)は基底をなすとし、x0=∑i=1ncivi で c1=0 とする。このとき冪乗法はすべての k で定義され、
r=λ1λ2,S=i=2∑n∣c1∣∣ci∣とおくと、Srk<1 を満たす k について
sin∠(xk,v1)≤1−SrkSrkが成り立つ。とくに sin∠(xk,v1)=O(rk) である。
証明(定理 3.2)
対角化可能なので固有ベクトルの基底がとれます(定理 3.2[対角化とジョルダン標準形])。Akvi=λikvi なので
Akx0=i=1∑nciλikvi.∣λ1∣>∣λ2∣≥0 より λ1=0 で、v1 の係数 c1λ1k は 0 でありません。vi は一次独立なので Akx0=0 で、冪乗法は止まりません。
yk=Akx0/(c1λ1k) とおくと
yk=v1+wk,wk=i=2∑nc1ci(λ1λi)kviで、∣λi/λ1∣≤r(i≥2)と ∥vi∥=1 から三角不等式で ∥wk∥≤Srk です。xk は yk の零でない定数倍なので、補題 2.2 の 3 と 2 から
sin∠(xk,v1)=sin∠(yk,v1)≤∥yk∥∥yk−v1∥=∥v1+wk∥∥wk∥.∥v1+wk∥≥1−∥wk∥≥1−Srk>0 なので主張が従います。
∎
収束率 r=∣λ2/λ1∣ は行列だけで決まり、初期ベクトルは定数 S にしか効きません。r=0.99 なら 1 桁合わせるのに約 230 回かかるので、冪乗法そのものは「最大固有値が他から十分離れている」ときにしか実用になりません。この遅さを直すのが §4 以降のシフトです。
方向が分かれば固有値も分かります。次の系は一般の(Hermite とは限らない)行列での評価です。
系 3.3(Rayleigh 商による固有値の評価(一般の行列))✓
Av=λv、v=0 とする。任意の x=0 について
∣RA(x)−λ∣≤∥A−λI∥sin∠(x,v).とくに 定理 3.2 の仮定のもとで ∣μk−λ1∣=O(rk) である。
証明(系 3.3)
補題 2.2 の 3 により ∥x∥=1 としてよく、Rayleigh 商も x の定数倍で変わりません。d=x−Pvx とおくと ∥d∥=sin∠(x,v) で、(A−λI)Pvx=0 です。したがって
RA(x)−λ=x∗(A−λI)x=x∗(A−λI)dで、Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∣RA(x)−λ∣≤∥x∥∥(A−λI)d∥≤∥A−λI∥∥d∥ です。後半は 定理 3.2 の評価を代入すれば従います。
∎
Hermite 行列では Rayleigh 商の精度が 1 段上がります。誤差が角の正弦の 1 乗ではなく 2 乗で効くのです。たとえば方向が 0.01 ずれていても、Rayleigh 商のずれは 0.0001 に比例する程度で済みます。次の補題は、ベクトルを固有ベクトルの成分とそれに直交する成分に分け、前者が固有値を正確に与え、後者のずれだけが誤差になることを示します。
補題 3.4(Rayleigh 商の 2 次の精度)✓
A を Hermite 行列とし、その最大・最小の固有値を λmax, λmin とする。Av=λv、∥v∥=1 とし、単位ベクトル x を x=αv+βw(∥w∥=1、v∗w=0、∣β∣=sin∠(x,v))と分解する。このとき
RA(x)−λ=sin2∠(x,v)(RA(w)−λ)であり、とくに ∣RA(x)−λ∣≤(λmax−λmin)sin2∠(x,v) である。(sin∠(x,v)=0 のときは w を v に直交する任意の単位ベクトルとする。)
証明(補題 3.4)
分解は αv=Pvx、βw=x−Pvx として作れ、補題 2.2 の 1 から ∣α∣2+∣β∣2=1 です。A は Hermite で λ は実数なので v∗Aw=(Av)∗w=λv∗w=0、同様に w∗Av=0 です。よって交差項が消えて
RA(x)=x∗Ax=∣α∣2λ+∣β∣2w∗Aw=(1−∣β∣2)λ+∣β∣2RA(w)となり、λ を移項すれば等式を得ます。RA(w) と λ はともに [λmin,λmax] に入る(定理 7.1[スペクトル定理] の 1)ので、∣RA(w)−λ∣≤λmax−λmin です。
∎
系 3.5(Hermite 行列の冪乗法)✓
A を Hermite 行列とし、固有値を ∣λ1∣>∣λ2∣≥⋯≥∣λn∣ と並べ、対応する正規直交な固有ベクトルを u1,…,un とする。u1∗x0=0 ならば、冪乗法の xk と μk=RA(xk) は、r=∣λ2/λ1∣ として
tan∠(xk,u1)≤rktan∠(x0,u1),∣μk−λ1∣≤(λmax−λmin)r2ktan2∠(x0,u1)を満たす。
証明(系 3.5)
正規直交な固有ベクトルの基底はスペクトル定理(定理 4.2[スペクトル定理])で存在します。ci=ui∗x0 とおくと Akx0=∑iciλikui で、z=∑i≥2ciλikui は u1 に直交します。補題 2.2 の 4 と 3(xk は Akx0 の正の定数倍)と、正規直交基底での三平方の定理から
tan∠(xk,u1)=∣c1∣∣λ1∣k(∑i≥2∣ci∣2∣λi∣2k)1/2≤rk∣c1∣(∑i≥2∣ci∣2)1/2=rktan∠(x0,u1)です。sin2≤tan2 と 補題 3.4 から後半が従います。
∎
一般の行列では方向の誤差と固有値の誤差が同じ rk でしたが、Hermite 行列では固有値の誤差が r2k になります。反復の回数で言えば、固有値は固有ベクトルの半分の回数で同じ桁数に達します。
# Python 3.12.10 / NumPy 2.5.2(x86-64, OpenBLAS)/乱数種 20261006
rng = np.random.default_rng(20261006)
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 の列
x = rng.standard_normal(n)
print(" k sin(角) |rho-l1| 10 回の比(角) 比(rho)")
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}")
print("予測:", f"{0.8**10:.3f}", f"{0.8**20:.4f}")
k sin(角) |rho-l1| 10 回の比(角) 比(rho)
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
検算. (a) 固有値は構成に使った値 10,8 と [−6,6] の等間隔点なので、r=0.8 が既知です。(b) 角は構成に使った固有ベクトル Q[:, 0] との間で測っています。(c) k≥30 で、10 回あたりの比は角で 0.107≈0.810、Rayleigh 商で 0.0114≈0.820 となり、定理 3.2 と 系 3.5 が予言する傾きに一致します。最初の 10 回の比がそれより小さいのは、初期ベクトルの誤差の大半が ∣λi∣≤6 の成分(10 回で 0.610≈0.006 倍に縮む)にあったためです。k=50 での評価式の右辺は (λmax−λmin)r100tan2∠(x0,u1)≈16×2.0×10−10×4.5≈1.5×10−8 で、実際の誤差 1.9×10−13 はその内側にあります。
上の表の誤差を片対数で描いたもの(横軸は反復回数 k、縦軸は誤差の常用対数)。角(細線)の傾きは log₁₀ 0.8、Rayleigh 商(太線)の傾きはその 2 倍になる。
仮定 ∣λ1∣>∣λ2∣ と対角化可能性を外すと何が壊れるかを、2 つの例で見ます。
例 3.6(絶対値の等しい固有値:回転)
A=(01−10)(90° の回転)の固有値は ±i で、絶対値が等しいので 定理 3.2 の仮定を満たしません。実際、x0=(1,0)T から始めると xk は (1,0)T,(0,1)T,(−1,0)T,(0,−1)T を周期 4 で回り、どの方向にも収束しません。Rayleigh 商はつねに 0 で、どちらの固有値にも近づきません。実行列の冪乗法は、絶対値最大の固有値が複素共役の対になると必ずこうなります。
例 3.7(Jordan 細胞:1/k の遅い収束)
J=(λ01λ)(λ=0)は対角化できません(Jordan 標準形 定理 6.2[対角化とジョルダン標準形] の 2 次の細胞そのもの)。Jk=λk(10k/λ1) なので、x0=(0,1)T に対し Jkx0=λk−1(k,λ)T です。唯一の固有ベクトルの方向 e1 とのなす角は tan∠=∣λ∣/k で、確かに 0 に近づきますが、その速さは幾何級数的ではなく 1/k です。1 桁合わせるのに反復回数を 10 倍にしなければなりません。
冪乗法が見つけるのは絶対値最大の固有値だけです。μ が固有値でないとき、(A−μI)−1 の固有値は 1/(λi−μ) で、固有ベクトルは A と同じです。μ に最も近い固有値が、(A−μI)−1 では絶対値最大の固有値になります。
定義 4.1(シフトつき逆反復)↗
μ∈C が A の固有値でないとき、x0=0 から
(A−μI)yk+1=xk,xk+1=∥yk+1∥yk+1を繰り返す反復を、シフト μ の逆反復という。
逆行列は作らず、A−μI を 1 回 LU 分解しておいて、毎回は前進・後退代入(O(n2))だけを行います。
定理 4.2(シフトつき逆反復の収束)✓
A は 定理 3.2 と同じく対角化可能で、固有値 λ1,…,λn(重複を込める)と単位固有ベクトル v1,…,vn が基底をなすとする。μ は固有値でなく、ある添字 J について
∣λJ−μ∣<∣λK−μ∣≤∣λi−μ∣(i=J)を満たす K=J があるとする(λJ は μ に最も近い固有値で単純、λK は 2 番目に近い固有値)。x0=∑icivi で cJ=0 ならば、
r=∣λK−μ∣∣λJ−μ∣,S=i=J∑∣cJ∣∣ci∣として、Srk<1 のとき sin∠(xk,vJ)≤Srk/(1−Srk) が成り立つ。A が Hermite 行列で μ が実数なら、さらに tan∠(xk,vJ)≤rktan∠(x0,vJ) である(vi は正規直交にとる)。
証明(定理 4.2)
B=(A−μI)−1 とおくと (A−μI)vi=(λi−μ)vi の両辺に B を掛けて (λi−μ)−1 で割れば Bvi=(λi−μ)−1vi です。よって B は同じ固有ベクトルの基底で対角化でき、固有値 νi=(λi−μ)−1 は仮定から ∣νJ∣>∣νK∣≥∣νi∣(i=J)を満たします。逆反復は B に対する冪乗法そのものなので、添字を J が先頭、K が 2 番目になるよう並べ替えて 定理 3.2 を使えば、比は ∣νK/νJ∣=r で前半が従います。A が Hermite で μ が実数なら B も Hermite(B∗=((A−μI)∗)−1=B)なので、系 3.5 の前半を B に使えば後半が従います。
∎
μ を λJ に近づけるほど r は 0 に近づきます。固有値のよい近似(たとえば二分法や QR 法で求めた値)があれば、逆反復は 2〜3 回で固有ベクトルを機械精度まで求めます。これが LAPACK で三重対角行列(対角とその両隣にしか非零成分がない行列)の固有ベクトルを求める古典的な方法でした。
# Python 3.12.10 / NumPy 2.5.2(x86-64, OpenBLAS)/乱数種 20261006
rng = np.random.default_rng(20261006)
lam = np.concatenate(([10.0, 8.0], np.linspace(-6, 6, n - 2)))
Q, _ = np.linalg.qr(rng.standard_normal((n, n)))
v = Q[:, 1] # 固有値 8 の固有ベクトル
x0 = rng.standard_normal(n)
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]:
x, ts = x0 / np.linalg.norm(x0), []
x = np.linalg.solve(M, 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}")
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) 固有値 8 に最も近い他の固有値は 10 と 6(どちらも距離 2)で、μ>8 のとき 2 番目に近いのは 10 なので、予測の比は ∣8−μ∣/∣10−μ∣ です。(b) 角は構成に使った固有ベクトルとの間で測っています。(c) 観測した t3/t2 は予測の比と 2 桁目まで合い、μ を 8 に 10 倍近づけるごとに比も約 10 分の 1 になります。観測値が予測よりわずかに小さいのは、3 番目以降に近い固有値の成分がまだ残っていて、それがさらに速く縮むためです。
逆反復のシフトを固定せず、毎回いまの近似ベクトルの Rayleigh 商に取り替えると、収束の次数そのものが上がります。
定義 5.1(Rayleigh 商反復)↗
単位ベクトル x0 から、μk=RA(xk) とおき、A−μkI が正則である限り
(A−μkI)yk+1=xk,xk+1=∥yk+1∥yk+1を繰り返す反復を Rayleigh 商反復という。μk が固有値になったら止める(固有値が求まっている)。
毎回シフトが変わるので、1 回ごとに A−μkI を分解し直す必要があります(密行列なら O(n3)、三重対角なら O(n))。その代わり、次の定理のとおり収束は 3 次です。
定理に出てくる δ は、注目する固有値 λ から最も近い別の固有値までの距離、Δ は固有値全体の幅です。仮定 Δsin2θ0<δ/2 は、Rayleigh 商の誤差の上限(補題 3.4 の評価)が、他の固有値との距離の半分より小さい、つまり μk が目的の固有値 λ に最も近いままである、ということを保証するためのものです。
定理 5.2(Rayleigh 商反復の局所 3 次収束)✓
A を n 次 Hermite 行列(n≥2)、λ をその単純な固有値、v を対応する単位固有ベクトルとし、
δ=min{∣λ′−λ∣:λ′ は λ と異なる A の固有値},Δ=λmax−λminとおく。Rayleigh 商反復の θk=∠(xk,v) が Δsin2θ0<δ/2 を満たすならば、反復が続く限り次が成り立つ。
- ∣μk−λ∣≤Δsin2θk であり、μk は λ 以外の固有値に一致しない。
- tanθk+1≤δ−Δsin2θkΔsin2θktanθk≤δ2Δtan3θk。とくに θk+1≤θk で、Δsin2θk+1<δ/2 も保たれる。
したがって反復が有限回で止まらなければ tanθk→0 であり、その収束は 3 次である。
証明(定理 5.2)
δ は [λmin,λmax] の 2 点の距離なので δ≤Δ で、仮定から sin2θ0<1/2、とくに cosθ0>0 です。k について帰納的に、Δsin2θk<δ/2 のもとで 1 と 2 を示します。
スペクトル定理(定理 4.2[スペクトル定理])で正規直交な固有ベクトル u1,…,un をとります。λ は単純なので、固有値が λ であるのはただ 1 つ uj だけで、v は uj の絶対値 1 の定数倍です(v を基底で展開すると、Av=λv から固有値が λ と異なる成分は 0)。よって v に直交するベクトル全体は {ui:i=j} で張られます。
反復は xk を絶対値 1 の複素数倍しても、μk も以後の方向も変わらないので、v∗xk=cosθk>0 としてかまいません。sinθk=0 なら xk は固有ベクトルで μk=λ となり反復は止まるので、sinθk>0 とし、xk=cosθkv+sinθkw(∥w∥=1、v∗w=0)と書きます。
1 の証明. 補題 3.4 から ∣μk−λ∣≤Δsin2θk<δ/2 です。λ と異なる固有値 λ′ については ∣λ′−μk∣≥∣λ′−λ∣−∣λ−μk∣>δ−δ/2>0 なので、μk=λ′ です。
2 の証明. 反復が続くので μk=λ で、1 と合わせて A−μkI は正則です。w=∑i=jaiui と書くと
(A−μkI)−1xk=λ−μkcosθkv+sinθki=j∑λi−μkaiuiで、第 2 項は v に直交し、そのノルムは sinθk∥w∥/mini=j∣λi−μk∣≤sinθk/(δ−∣λ−μk∣) 以下です(1 の証明と同じ三角不等式)。補題 2.2 の 4 と 3 から
tanθk+1≤cosθk/∣λ−μk∣sinθk/(δ−∣λ−μk∣)=δ−∣λ−μk∣∣λ−μk∣tanθk.関数 t↦t/(δ−t) は [0,δ) で増加なので、∣λ−μk∣≤Δsin2θk を代入して 2 の最初の不等式を得ます。Δsin2θk<δ/2 から分母は δ/2 より大きく、分子の係数は 1 未満なので tanθk+1<tanθk、すなわち θk+1<θk で、sin2θk+1<sin2θk から帰納法の仮定が保たれます。最後に sin2θk≤tan2θk と分母 >δ/2 から (2Δ/δ)tan3θk の評価を得ます。
最後に収束を示します。ck=Δsin2θk/(δ−Δsin2θk) とおくと、2 から tanθk+1≤cktanθk で、θk は減少するので ck≤c0<1 です(c0<1 は Δsin2θ0<δ/2 と同値です)。よって tanθk≤c0ktanθ0→0 であり、各段で tanθk+1≤(2Δ/δ)tan3θk が成り立つので収束は 3 次です。
∎
正接が 1 回ごとに 3 乗になるので、正しい桁数は 1 回ごとにほぼ 3 倍になります。単根の近くの Newton 法が桁数を 1 回ごとにほぼ 2 倍にする(非線形方程式の反復解法)のと比べて、1 段高い次数です。固有値の誤差は 1 から sin2θk の程度なので、こちらの桁数も 1 回ごとにほぼ 3 倍です。
次のコードは、対角が 2・副対角が −1 の 6 次の三重対角行列(固有値は 2−2cos(jπ/7)、固有ベクトルは (sin(ijπ/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)
A = 2 * mp.eye(n) # 対角 2、副対角 -1 の三重対角行列
A[i, i + 1] = A[i + 1, i] = -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|")
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)
An, xd = np.array(A.tolist(), dtype=float), np.ones(n) / np.sqrt(n)
print("倍精度の |mu_k - 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)
検算. (a) 固有値と固有ベクトルは閉じた形の値と比べています。(b) 倍精度の ∣μk−λ∣ は 100 桁の値と表示した桁で一致し、3 回目で倍精度の丸めの水準に達します。(c) 正接の指数は −1.5→−4.5→−13.4→−40.2 とほぼ 3 倍ずつ、固有値の誤差の指数は −2.8→−8.8→−26.7→−80.3 とやはりほぼ 3 倍ずつ進み、比 tanθk+1/tan3θk は 1.11,1.02,1.00,1.00 です。この行列では Δ=4cos(π/7)≈3.60、δ=2cos(π/7)−2cos(2π/7)≈0.555 で、定理の係数 2Δ/δ≈13 は実際の比より大きく、評価として正しい側にあります。なお k=0 では Δsin2θ0≈0.31>δ/2 で定理の仮定を満たしませんが、k=1 以降は満たしています(定理は十分条件です)。
1 本のベクトルの代わりに p 本の正規直交なベクトルを同時に反復し、毎回直交化し直すと、上位 p 個の固有値に対応する不変部分空間が求まります。直交化しないと、p 本すべてが同じ v1 の方向に寄ってしまい、数値的に一次従属になります。
次の定理は記号が多いので、先に役割を整理します。X は固有ベクトルを並べた行列で、その前半 X1(上位 p 個の固有ベクトル)と後半 X2 に分けます。X−1 の対応する行を W1∗,W2∗ と書くと、W1∗Z0 は初期の p 本が上位 p 個の固有ベクトルの成分をもつかどうかを表す p 次の行列で、これが正則であることが c1=0 の一般化です。証明の方針は冪乗法と同じで、AkZ0 を上位の固有ベクトルの部分と残りの部分に分け、残りが比 ∣λp+1/λp∣k で小さくなることを示します。
定義 6.1(同時反復)↗
Z0∈Cn×p を正規直交な列をもつ行列とする。Yk=AZk−1 の列が一次独立である限り、その列に Gram–Schmidt の直交化(定理 6.1[内積空間とグラム・シュミット直交化])を施して
Yk=ZkRk(Zk は正規直交な列をもち、Rk は対角成分が正の上三角行列)と書き、Zk を次の反復ベクトルとする。これを同時反復(部分空間反復、直交反復)という。p=1 のときは冪乗法である。
Gram–Schmidt の定理の 2(先頭の j 本が張る空間が一致する)と 3(Rk の対角成分 zi∗yi が正)により、この分解は一意で、Zk の先頭 j 列の張る空間は Yk の先頭 j 列の張る空間に等しくなります。実際の計算では Householder 変換(鏡映を順に掛けて列を直交化する方法)による QR 分解を使いますが、張る空間は同じです。
定理 6.2(同時反復の収束)✓
A=XΛX−1 を対角化可能な行列、Λ=diag(λ1,…,λn) とし、ある 1≤p<n について
∣λ1∣≥⋯≥∣λp∣>∣λp+1∣≥⋯≥∣λn∣とする。X=(X1 X2)(X1 は n×p)と分け、X−1 を (W1∗W2∗)(W1 は n×p)と分ける。Xp=spanX1(上位 p 個の固有値の不変部分空間)とし、G=W1∗Z0 が正則であると仮定する。このとき同時反復はすべての k で定義され、spanZk=spanAkZ0 であり、
r=λpλp+1,C=∥W1∥∥X2∥∥W2∗Z0G−1∥とおくと、Crk<1 のとき
d(spanZk,Xp)≤1−CrkCrkが成り立つ。
証明(定理 6.2)
準備. X−1X=I を区分けすると W1∗X1=Ip、W1∗X2=O で、X−1A=ΛX−1 を区分けすると W1∗A=Λ1W1∗ です。ここで Λ1=diag(λ1,…,λp)、Λ2=diag(λp+1,…,λn) とおきました。∣λp∣>∣λp+1∣≥0 から Λ1 は正則です。
反復が定義されること. k についての帰納法で、Zk=AkZ0Mk(Mk は正則な p 次行列)と書けることを示します。k=0 は M0=I です。Zk−1=Ak−1Z0Mk−1 なら Yk=AkZ0Mk−1 で、
W1∗Yk=W1∗AkZ0Mk−1=Λ1kW1∗Z0Mk−1=Λ1kGMk−1は正則です。Ykc=0 なら W1∗Ykc=0 から c=0 なので、Yk の列は一次独立で、Gram–Schmidt の直交化が実行できます。Yk=ZkRk の Rk は正則なので Zk=AkZ0Mk−1Rk−1 で、Mk=Mk−1Rk−1 とすればよい。とくに spanZk=spanAkZ0 です。
隔たりの評価. Ak=XΛkX−1 より
AkZ0=X1Λ1kW1∗Z0+X2Λ2kW2∗Z0で、右から正則な G−1Λ1−k を掛けると
AkZ0G−1Λ1−k=X1+Ek,Ek=X2Λ2k(W2∗Z0G−1)Λ1−kです(G−1Λ1−k は Λ1kG の逆行列で、X1Λ1kGG−1Λ1−k=X1)。対角行列のノルムは対角成分の絶対値の最大なので ∥Λ2k∥=∣λp+1∣k、∥Λ1−k∥=∣λp∣−k で、
∥Ek∥≤∥X2∥∣λp+1∣k∥W2∗Z0G−1∥∣λp∣−k,∥W1∥∥Ek∥≤Crkです。spanZk=span(X1+Ek) なので、その零でない元は y=(X1+Ek)c(c=0)と書けます。W1∗X1=Ip から ∥c∥=∥W1∗X1c∥≤∥W1∥∥X1c∥ なので(ここで ∥W1∗∥≤∥W1∥ を使いました。任意の u について Cauchy–Schwarz の不等式から ∥W1∗u∥2=u∗W1(W1∗u)≤∥u∥∥W1∥∥W1∗u∥ で、W1∗u=0 なら両辺を ∥W1∗u∥ で割れば従います)
∥y∥≥∥X1c∥−∥Ekc∥≥(∥W1∥1−∥Ek∥)∥c∥>0です(∥W1∥∥Ek∥≤Crk<1)。X1c∈Xp なので 補題 2.2 の 2 から
sin∠(y,Xp)≤∥y∥∥y−X1c∥=∥y∥∥Ekc∥≤1/∥W1∥−∥Ek∥∥Ek∥=1−∥W1∥∥Ek∥∥W1∥∥Ek∥≤1−CrkCrkで、y について上限をとれば主張を得ます(t↦t/(1−t) は増加)。
∎
p=1 とすると、W1∗ は左固有ベクトル(X−1 の第 1 行)で、G=0 は「初期ベクトルが v1 の成分をもつ」という 定理 3.2 の条件 c1=0 と同じです。
同時反復の要は、Gram–Schmidt の直交化が先頭の列から順に張る空間を保つことです。このため、先頭の j 列だけを見ると、それ自体が j 本の同時反復になっています。
系 6.3(先頭の列の収束)✓
A=XΛX−1 を 定理 6.2 と同じく対角化可能とし、固有値は ∣λ1∣≥⋯≥∣λn∣ と並んでいるとする。Z0∈Cn×p(1≤p≤n)は正規直交な列をもち、同時反復がすべての k で定義されるとする(たとえば A が正則なら AZk−1 の列はつねに一次独立)。1≤j≤p、j<n について ∣λj∣>∣λj+1∣ を仮定し、Z0 の先頭 j 列を Z0(j)、X−1 の先頭 j 行を (W1(j))∗ と書いて、(W1(j))∗Z0(j) が正則とする。このとき Zk の先頭 j 列 Zk(j) は、Z0(j) から始めた同時反復の k 回目と一致し、その張る空間は Xj=span{x1,…,xj}(X の先頭 j 列)に O(∣λj+1/λj∣k) で近づく。
証明(系 6.3)
k についての帰納法で、Zk(j) が Z0(j) から始めた同時反復の k 回目 Z~k に一致することを示します。Zk−1(j)=Z~k−1 なら、Yk=AZk−1 の先頭 j 列は Y~k=AZ~k−1 です。Gram–Schmidt の直交化の漸化式(定理 6.1[内積空間とグラム・シュミット直交化])は第 i 列を作るのに第 i 列までしか使わないので、Yk を直交化した結果の先頭 j 列は Y~k を直交化した結果に一致し、Zk(j)=Z~k です(Y~k は Yk の先頭 j 列なので、その列は一次独立です)。あとは p を j に置き換えた 定理 6.2 を Z~k に使えば、隔たりは O(∣λj+1/λj∣k) です。
∎
したがって ∣λ1∣>∣λ2∣>⋯>∣λp∣>∣λp+1∣ なら、Zk の第 1 列は x1 の方向に、先頭 2 列は span{x1,x2} に、… と、入れ子になった不変部分空間の列がそれぞれの比の速さで同時に求まります。p=n としたものを行列の相似変換の言葉で書き直したのが QR 法で、その収束の解析は QR 法の章でこの系をもとに行います。
反復を止めたとき、手元にあるのは近似固有対 (μ,x) です。真の固有ベクトルは分からないので、角 θ そのものは測れません。測れるのは残差です。
命題 7.1(Rayleigh 商は残差を最小にする)✓
A∈Cn×n、∥x∥=1 とし、r=Ax−RA(x)x とおく。このとき x∗r=0、∥r∥2=∥Ax∥2−∣RA(x)∣2 であり、任意の μ∈C について
∥Ax−μx∥2=∥r∥2+∣RA(x)−μ∣2.とくに μ↦∥Ax−μx∥ は μ=RA(x) でだけ最小になる。
証明(命題 7.1)
R=RA(x)=x∗Ax とおくと x∗r=x∗Ax−Rx∗x=0 です。Ax−μx=r+(R−μ)x の 2 項は直交するので、ノルムの 2 乗は ∥r∥2+∣R−μ∣2 です。μ=0 とおけば ∥Ax∥2=∥r∥2+∣R∣2 を得ます。
∎
系 7.2(実対称行列の残差評価)✓
A を実対称行列、x∈Rn を単位ベクトル、r=Ax−RA(x)x とする。このとき ∣λ−RA(x)∣≤∥r∥ を満たす A の固有値 λ がある。∥r∥ は、すべての実数 μ についての残差 ∥Ax−μx∥ の中で最小である。
証明(系 7.2)
A は実対称なので RA(x) は実数で、残差による固有値の評価(定理 6.1[行列計算の認証と最適化の証明書])を μ=RA(x) に使えば前半です。後半は 命題 7.1 です。
∎
一般の行列では、残差が小さくても固有値が近いとは限りません。差を測るのは固有ベクトル行列の条件数です。
命題 7.3(対角化可能な行列の残差評価)✓
A=VΛV−1 を対角化可能な行列とし、κ(V)=∥V∥∥V−1∥ とおく。∥x∥=1、μ∈C、r=Ax−μx ならば、mini∣λi−μ∣≤κ(V)∥r∥ である。
証明(命題 7.3)
E=−rx∗ とおくと (A+E)x=Ax−r(x∗x)=μx なので、μ は A+E の固有値です。任意の y について Cauchy–Schwarz の不等式(定理 4.2[内積空間とグラム・シュミット直交化])から ∥Ey∥=∣x∗y∣∥r∥≤∥r∥∥y∥ なので ∥E∥≤∥r∥ です。Bauer–Fike の定理(定理 6.2[行列計算の認証と最適化の証明書])から mini∣μ−λi∣≤κ(V)∥E∥≤κ(V)∥r∥ です。
∎
Hermite 行列では V をユニタリにとれて κ(V)=1 なので、命題 7.3 は 系 7.2 の複素版になっています。∥r∥ の評価は 1 次ですが、固有値の隔たりについての情報が加われば 2 次に改善できます。その代表が Temple の不等式です。
例 7.4(Rayleigh 商反復 1 回のあとの包み込み)
§5 のコードの 6 次の三重対角行列で、x0=(1,…,1)T/6 から Rayleigh 商反復を 1 回行うと(mpmath 30 桁で計算)、η=RA(x1)≈0.199503、∥r∥≈0.0457 です。
- 系 7.2 からは、区間 [η−∥r∥,η+∥r∥]≈[0.154,0.245] に固有値があることが分かります。
- 最小の固有値を λ1、2 番目に小さい固有値を λ2 とします。λ2≥0.75 が分かっていれば(ここでは閉じた形 2−2cos(2π/7)≈0.75302 から。実際の計算では Sturm 列の数え上げなどで確かめます)、Temple の不等式(定理 4.1[無限次元の固有値の包含])を ν=0.75 で使えます。命題 7.1 から Temple の不等式の σ2=∥Ax∥2−η2 は ∥r∥2≈0.00209 に等しく、λ1≥η−∥r∥2/(ν−η)≈0.195712 です。
- 上からは Rayleigh 商そのものが λ1≤η を与えます(定理 7.1[スペクトル定理])。
こうして λ1∈[0.19571,0.19951] と、残差だけの評価より約 24 倍狭い区間に包めます。真の値は 2−2cos(π/7)≈0.198062 です。Temple の評価の幅 ∥r∥2/(ν−η) は残差の 2 乗に比例するので、反復をもう 1 回進めれば幅は桁違いに縮みます。
演習 8.1易
実対称行列で、絶対値最大の固有値と 2 番目の固有値の比が ∣λ2/λ1∣=0.9、初期ベクトルの角が tan∠(x0,u1)=1 であるとする。冪乗法の Rayleigh 商の誤差を (λmax−λmin)×10−10 以下にするには、系 3.5 の評価で何回の反復が十分か。固有ベクトルの角の正接を 10−10 以下にする回数とも比べよ。
解答
系 3.5 より ∣μk−λ1∣≤(λmax−λmin)0.81k なので、0.81k≤10−10、すなわち k≥10/log10(1/0.81)≈10/0.0915≈109.3 から k=110 回で十分です。角の正接は tan∠(xk,u1)≤0.9k なので k≥10/log10(1/0.9)≈218.5、つまり 219 回が必要です。固有値は固有ベクトルの約半分の回数で同じ桁数に達します。
演習 8.2標準
実対称行列の固有値が 1,2,…,10 であるとする。A−σI(σ は実数)に冪乗法を使って固有値 10 を求めるとき、収束率 maxi≤9∣i−σ∣/∣10−σ∣ を最小にする σ と、そのときの収束率を求めよ。シフトなし(σ=0)の場合と比べて、同じ精度に要する反復回数は何倍になるか。
解答
σ<10 の範囲で考えます(σ が 10 に近すぎると固有値 1 の側より 10 の側が小さくなり、絶対値最大でなくなる)。i≤9 について ∣i−σ∣ の最大は max(∣9−σ∣,∣1−σ∣) で、これは σ=5 で最小値 4 をとります。収束率 f(σ)=max(∣9−σ∣,∣1−σ∣)/(10−σ) は、σ≤5 では (9−σ)/(10−σ)=1−1/(10−σ) で σ について減少、5≤σ<10 では (σ−1)/(10−σ) で増加なので、最小は σ=5 の 4/5=0.8 です。シフトなしでは 0.9 なので、反復回数の比は log0.9/log0.8≈0.472、つまり約半分で済みます。
演習 8.3標準
A=diag(a,b)(a=b は実数)と単位ベクトル x=(cosθ,sinθ)T(0<θ<π/2)について、Rayleigh 商反復を 1 回行ったあとのベクトルと e1 のなす角 θ′ が tanθ′=tan3θ を満たすことを示せ(RA(x) が a,b に一致しない場合)。定理 5.2 の評価の係数 2Δ/δ と比べると、この場合の係数は何か。
解答
μ=acos2θ+bsin2θ なので a−μ=(a−b)sin2θ、b−μ=−(a−b)cos2θ です。よって
(A−μI)−1x=a−b1(sin2θcosθ, −cos2θsinθ)Tで、補題 2.2 の 4 から tanθ′=cosθ/sin2θsinθ/cos2θ=cos3θsin3θ=tan3θ です。係数は 1 で、定理の係数 2Δ/δ=2(2 次では Δ=δ=∣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 章。
執筆: Claude Code(Opus 5.5)、Claude Code(Sonnet 5.5)・監査 2 回(Claude Code)・最終 2026-10-12