べき乗法とは何か？

べき乗法（Power Method）とは、行列 A の絶対値最大の固有値とその固有ベクトルを、行列とベクトルの積を繰り返すことで求める反復法である。初期ベクトルに A を繰り返し掛け、正規化することで、支配的な固有ベクトルの方向に収束する。収束速度は第2固有値と第1固有値の比 |λ₂/λ₁| に依存する。

逆べき乗法（逆反復法）とは何か？

逆べき乗法は、(A − σI) の逆行列にべき乗法を適用する手法である。シフト σ に最も近い固有値に収束し、σ が固有値に近いほど高速に収束する。実装では逆行列を明示的に求めず、LU分解を利用して連立方程式を解くことで効率的に計算する。

レイリー商反復法とは何か？

レイリー商反復法は、逆反復法のシフト σ を毎反復ごとにレイリー商 R(x) = xᵀAx / xᵀx で更新する手法である。対称行列に対して3次収束（cubic convergence）を達成し、通常3〜4回の反復で機械精度に到達する。

第23章: べき乗法と逆反復法

Power Method and Inverse Iteration

1. 導入と動機

$N \times N$ 行列 $A$ の固有値問題 $A\mathbf{v} = \lambda \mathbf{v}$ の数値解法には、大きく分けて2つのアプローチがある。

全固有値を同時に求める手法: QR法（第5章参照）
特定の固有値を選択的に求める反復法: べき乗法とその変形

本章では (2) のアプローチを詳しく解説する。べき乗法は絶対値最大の固有値を求める最も基本的な手法であり、逆反復法はその変形として任意の目標値に最も近い固有値を効率的に求めることができる。

別名

べき乗法は冪乗法、累乗法、乗冪法とも呼ばれる。Richard von Mises（1929年）によって体系化された。

2. べき乗法

2.1 前提条件

べき乗法の前提条件

$A$ は対角化可能である（固有ベクトルが $\mathbb{R}^N$（または $\mathbb{C}^N$）の基底をなす）。
固有値を絶対値の大きい順に並べたとき、優位固有値が唯一: $$|\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \cdots \geq |\lambda_N|$$
初期ベクトル $\mathbf{x}^{(0)}$ を固有ベクトルで展開したとき $$\mathbf{x}^{(0)} = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \cdots + c_N \mathbf{v}_N$$ の係数が $c_1 \neq 0$ を満たす。

$c_1 \neq 0$ の意味

初期ベクトルが $\mathbf{v}_1$ 方向の成分を持たなければ、べき乗法は $\lambda_1$ に収束しない。対称行列では固有ベクトルが直交するため $c_1 = \mathbf{x}^{(0)} \cdot \mathbf{v}_1 / \|\mathbf{v}_1\|^2$ であり、$\mathbf{x}^{(0)} \cdot \mathbf{v}_1 \neq 0$ と同値になる。しかし一般の行列では固有ベクトルは直交しないため、$c_1 \neq 0$ が正確な条件である。実用上は、ランダムな初期ベクトルを選べば $c_1 = 0$ となることはまずないため問題にならない。

2.2 アルゴリズム

べき乗法のアルゴリズム

入力: $N \times N$ 行列 $A$、初期ベクトル $\mathbf{x}^{(0)}$（$\|\mathbf{x}^{(0)}\| = 1$）、収束判定閾値 $\varepsilon$、最大反復回数 $k_{\max}$

$k = 1, 2, \ldots, k_{\max}$ に対して:
1. $\mathbf{y}^{(k)} = A \mathbf{x}^{(k-1)}$（行列ベクトル積: $O(N^2)$）
2. $\mathbf{x}^{(k)} = \mathbf{y}^{(k)} / \|\mathbf{y}^{(k)}\|$（正規化）
3. $\mu^{(k)} = (\mathbf{x}^{(k)})^T A \mathbf{x}^{(k)}$（レイリー商で固有値近似）
4. $|\mu^{(k)} - \mu^{(k-1)}| < \varepsilon$ ならば終了

出力: 固有値の近似 $\mu^{(k)}$、固有ベクトルの近似 $\mathbf{x}^{(k)}$

3. 収束の証明

3.1 基本性質: $A^k \mathbf{v}_n = \lambda_n^k \mathbf{v}_n$

固有値方程式 $A\mathbf{v}_n = \lambda_n \mathbf{v}_n$ から、帰納法により次式が成り立つ。

\begin{equation} A^k \mathbf{v}_n = \lambda_n^k \mathbf{v}_n \quad (k = 1, 2, 3, \ldots) \label{eq:Akvn} \end{equation}

証明: $k$ について帰納法を用いる。$k = 1$ は定義より成立。$k-1$ で成立すると仮定すると、

$$A^k \mathbf{v}_n = A(A^{k-1} \mathbf{v}_n) = A(\lambda_n^{k-1} \mathbf{v}_n) = \lambda_n^{k-1} (A\mathbf{v}_n) = \lambda_n^{k-1} (\lambda_n \mathbf{v}_n) = \lambda_n^k \mathbf{v}_n$$

3.2 固有ベクトル展開による収束の導出

べき乗法の収束定理

$A$ が対角化可能で $|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_N|$ を満たすとする。初期ベクトルの固有ベクトル展開で $c_1 \neq 0$ ならば、正規化付きべき乗法で生成される列 $\{\mathbf{x}^{(k)}\}$ は $\lambda_1$ に対応する固有ベクトルの方向に収束する。

証明: 初期ベクトルを固有ベクトルで展開する。

$$\mathbf{x}^{(0)} = c_1 \mathbf{v}_1 + c_2 \mathbf{v}_2 + \cdots + c_N \mathbf{v}_N, \quad c_1 \neq 0$$

$A^k$ を作用させると、式 $\eqref{eq:Akvn}$ より

\begin{align} A^k \mathbf{x}^{(0)} &= c_1 \lambda_1^k \mathbf{v}_1 + c_2 \lambda_2^k \mathbf{v}_2 + \cdots + c_N \lambda_N^k \mathbf{v}_N \\ &= c_1 \lambda_1^k \left[ \mathbf{v}_1 + \frac{c_2}{c_1}\left(\frac{\lambda_2}{\lambda_1}\right)^k \mathbf{v}_2 + \cdots + \frac{c_N}{c_1}\left(\frac{\lambda_N}{\lambda_1}\right)^k \mathbf{v}_N \right] \label{eq:expansion} \end{align}

$$\lim_{k \to \infty} \left(\frac{\lambda_i}{\lambda_1}\right)^k = 0 \quad (i = 2, 3, \ldots, N)$$

したがって $k \to \infty$ で括弧内は $\mathbf{v}_1$ に収束し、正規化後のベクトルは

$$\lim_{k \to \infty} \frac{A^k \mathbf{x}^{(0)}}{\|A^k \mathbf{x}^{(0)}\|} = \frac{\mathbf{v}_1}{\|\mathbf{v}_1\|}$$

を満たす。$\square$

3.3 固有値の抽出

$\|\mathbf{x}^{(k)}\| = 1$ を保つ正規化により、レイリー商で固有値を近似できる。

$$\mu^{(k)} = (\mathbf{x}^{(k)})^T A \mathbf{x}^{(k)}$$

$\mathbf{x}^{(k)} \to \mathbf{v}_1 / \|\mathbf{v}_1\|$ のとき、$A\mathbf{x}^{(k)} \to \lambda_1 \mathbf{x}^{(k)}$ であるから

$$\lim_{k \to \infty} \mu^{(k)} = \lambda_1 \|\mathbf{x}^{(k)}\|^2 = \lambda_1$$

3.4 正規化の必要性

4. 収束速度の解析

べき乗法の収束速度

固有ベクトルの収束: $\mathbf{x}^{(k)}$ と真の固有ベクトルの誤差は

$$\|\mathbf{x}^{(k)} - (\pm \mathbf{v}_1 / \|\mathbf{v}_1\|)\| = O\!\left(\left|\frac{\lambda_2}{\lambda_1}\right|^k\right)$$

固有値の収束: レイリー商 $\mu^{(k)}$ の誤差は

$$|\mu^{(k)} - \lambda_1| = O\!\left(\left|\frac{\lambda_2}{\lambda_1}\right|^{2k}\right)$$

レイリー商は固有ベクトルの誤差に対して2次の精度を持つため、固有値の収束は固有ベクトルの2倍の速さとなる。

収束比 $r = |\lambda_2/\lambda_1|$ は0と1の間の値を取り、$r$ が小さいほど収束が速い。$k$ 反復後の固有値の相対誤差を $\varepsilon$ 以下にするには

$$k \geq \frac{\log \varepsilon}{2 \log r}$$

回の反復が必要である。

反復回数の見積もり

$\lambda_1 = 10,\ \lambda_2 = 9$ のとき $r = 0.9$ である。固有値を6桁（$\varepsilon = 10^{-6}$）の精度で求めるには

$$k \geq \frac{-6}{2 \log_{10} 0.9} = \frac{6}{2 \times 0.0458} \approx 66 \text{ 回}$$

$|\lambda_2 / \lambda_1|$ が1に近いと多くの反復が必要になる。

5. 逆べき乗法（逆反復法）

5.1 原理

$A$ の固有値が $\lambda_1, \lambda_2, \ldots, \lambda_N$ であるとき、$A^{-1}$ の固有値は $1/\lambda_1, 1/\lambda_2, \ldots, 1/\lambda_N$ であり、固有ベクトルは同じである。したがって $A^{-1}$ にべき乗法を適用すれば、$|1/\lambda_i|$ が最大、すなわち $|\lambda_i|$ が最小の固有値に対応する固有ベクトルに収束する。

逆べき乗法（Inverse Iteration）

$A^{-1}$ にべき乗法を適用する。すなわち $\mathbf{y}^{(k)} = A^{-1} \mathbf{x}^{(k-1)}$ を反復する。実装では逆行列を明示的に計算せず、連立方程式 $A\mathbf{y}^{(k)} = \mathbf{x}^{(k-1)}$ を解く。

5.2 アルゴリズム

逆べき乗法のアルゴリズム

前処理: $A = LU$（LU分解を1回だけ計算: $O(N^3)$）

$k = 1, 2, \ldots$ に対して:
1. $L\mathbf{z} = \mathbf{x}^{(k-1)}$ を前進代入で解く（$O(N^2)$）
2. $U\mathbf{y}^{(k)} = \mathbf{z}$ を後退代入で解く（$O(N^2)$）
3. $\mathbf{x}^{(k)} = \mathbf{y}^{(k)} / \|\mathbf{y}^{(k)}\|$
4. $\mu^{(k)} = (\mathbf{x}^{(k)})^T A \mathbf{x}^{(k)}$（レイリー商）

出力: 絶対値最小の固有値 $\mu^{(k)}$ と対応する固有ベクトル $\mathbf{x}^{(k)}$

計算量

LU分解は1回の $O(N^3)$ で済み、各反復は前進・後退代入の $O(N^2)$ のみである。べき乗法と同じ反復コストで、絶対値最小の固有値を求められる。

6. シフト付き逆べき乗法

6.1 シフトの導入

目標のシフト値 $\sigma$ を導入し、行列 $(A - \sigma I)^{-1}$ にべき乗法を適用する。$(A - \sigma I)$ の固有値は $\lambda_i - \sigma$ であるから、$(A - \sigma I)^{-1}$ の固有値は

$$\frac{1}{\lambda_i - \sigma}$$

である。$\sigma$ に最も近い固有値 $\lambda_j$ に対して $|1/(\lambda_j - \sigma)|$ が最大となるため、べき乗法は $\lambda_j$ の固有ベクトルに収束する。

シフト付き逆べき乗法（Shifted Inverse Iteration）

シフト $\sigma$ に最も近い固有値とその固有ベクトルを求める。反復では連立方程式

$$(A - \sigma I)\mathbf{y}^{(k)} = \mathbf{x}^{(k-1)}$$

を解き、$\mathbf{x}^{(k)} = \mathbf{y}^{(k)} / \|\mathbf{y}^{(k)}\|$ で正規化する。

6.2 アルゴリズム

シフト付き逆べき乗法のアルゴリズム

入力: $A$、シフト $\sigma$、$\mathbf{x}^{(0)}$（$\|\mathbf{x}^{(0)}\| = 1$）、$\varepsilon$

前処理: $A - \sigma I = LU$（LU分解を1回）

$k = 1, 2, \ldots$ に対して:
1. $(A - \sigma I)\mathbf{y}^{(k)} = \mathbf{x}^{(k-1)}$ を前進・後退代入で解く
2. $\mathbf{x}^{(k)} = \mathbf{y}^{(k)} / \|\mathbf{y}^{(k)}\|$
3. $\mu^{(k)} = (\mathbf{x}^{(k)})^T A \mathbf{x}^{(k)}$
4. $|\mu^{(k)} - \mu^{(k-1)}| < \varepsilon$ ならば終了

出力: $\sigma$ に最も近い固有値 $\mu^{(k)}$ と固有ベクトル $\mathbf{x}^{(k)}$

6.3 収束速度

シフト付き逆べき乗法の収束速度

$\lambda_j$ を $\sigma$ に最も近い固有値、$\lambda_{\text{next}}$ を2番目に近い固有値とすると、固有ベクトルの収束比は

$$r = \left|\frac{\lambda_j - \sigma}{\lambda_{\text{next}} - \sigma}\right|$$

$\sigma$ が $\lambda_j$ に近いほど $r$ は小さくなり、収束が加速する。

シフトの選び方

Gershgorin 円定理で固有値のおおよその範囲を推定する
QR法の途中経過で得られる固有値の近似値を使う
三重対角化（Householder 変換）後の対角要素を使う

LAPACK では QR法で固有値を求めた後、シフト付き逆反復で固有ベクトルを計算するのが標準的な戦略である。

Gershgorin 円定理によるシフトの推定

$A = \begin{pmatrix}2&1\\1&3\end{pmatrix}$ の場合、Gershgorin 円は

行1: 中心 $2$、半径 $|1| = 1$ → 円 $[1, 3]$
行2: 中心 $3$、半径 $|1| = 1$ → 円 $[2, 4]$

すべての固有値は $[1, 4]$ に含まれる。第2固有値を狙うなら $\sigma = 1.3$（左端付近）、第1固有値を狙うなら $\sigma = 3.5$（右端付近）のようにシフトを選ぶ。

7. レイリー商反復法

7.1 レイリー商の性質

レイリー商（Rayleigh Quotient）

$$R(\mathbf{x}) = \frac{\mathbf{x}^T A \mathbf{x}}{\mathbf{x}^T \mathbf{x}}$$

$\mathbf{x}$ が固有ベクトル $\mathbf{v}_j$ に一致するとき $R(\mathbf{v}_j) = \lambda_j$ となる。

レイリー商の2次近似性

$\mathbf{x} = \mathbf{v}_j + \boldsymbol{\varepsilon}$ のとき

$$R(\mathbf{x}) = \lambda_j + O(\|\boldsymbol{\varepsilon}\|^2)$$

すなわちレイリー商は固有ベクトルの誤差に対して2次の精度を持つ。

7.2 アルゴリズム

シフト付き逆べき乗法のシフト $\sigma$ を、毎反復ごとにレイリー商で更新する。

レイリー商反復法

$k = 1, 2, \ldots$ に対して:
1. $\sigma^{(k)} = R(\mathbf{x}^{(k-1)}) = (\mathbf{x}^{(k-1)})^T A \mathbf{x}^{(k-1)}$（シフトを更新）
2. $(A - \sigma^{(k)} I)\mathbf{y}^{(k)} = \mathbf{x}^{(k-1)}$ を解く
3. $\mathbf{x}^{(k)} = \mathbf{y}^{(k)} / \|\mathbf{y}^{(k)}\|$
4. $|\sigma^{(k)} - \sigma^{(k-1)}| < \varepsilon$ ならば終了

レイリー商反復法の収束速度

一般の行列: 2次収束（quadratic convergence）
対称行列: 3次収束（cubic convergence）

対称行列の場合、3〜4回の反復で機械精度に到達する。

利点と欠点

利点: 極めて高速な収束（対称行列で3次）
欠点: シフトが毎回変わるため、LU分解を毎反復で再計算する必要がある（$O(N^3)$/反復）
欠点: どの固有値に収束するかは初期ベクトルに依存し、事前に予測できない

実用上は、シフト付き逆反復を数回行って固有ベクトルの近似を得てからレイリー商反復に切り替えるのが効果的である。

8. 数値例

8.1 べき乗法

例1: 2×2 行列のべき乗法

行列

$$A = \begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}$$

の固有値は $\lambda_1 = \frac{5+\sqrt{5}}{2} \approx 3.618$、$\lambda_2 = \frac{5-\sqrt{5}}{2} \approx 1.382$ であり、収束比は $|\lambda_2/\lambda_1| \approx 0.382$ である。

初期ベクトル $\mathbf{x}^{(0)} = (1/\sqrt{2},\ 1/\sqrt{2})^T$ として反復した結果を示す。

$k$	$\mathbf{x}^{(k)}$	$\mu^{(k)}$	$\|\mu^{(k)} - \lambda_1\|$
0	$(0.7071,\ 0.7071)$	$3.5000$	$1.2 \times 10^{-1}$
1	$(0.6000,\ 0.8000)$	$3.6000$	$1.8 \times 10^{-2}$
2	$(0.5547,\ 0.8321)$	$3.6154$	$2.7 \times 10^{-3}$
3	$(0.5369,\ 0.8436)$	$3.6176$	$3.8 \times 10^{-4}$
4	$(0.5300,\ 0.8480)$	$3.6180$	$5.5 \times 10^{-5}$

$\mu^{(k)}$ の誤差は毎反復で約 $0.382^2 \approx 0.146$ 倍に減少しており、理論通りの収束を確認できる。

8.2 シフト付き逆べき乗法

例2: シフト付き逆べき乗法

同じ行列 $A$ に対し、シフト $\sigma = 1.3$ として第2固有値 $\lambda_2 \approx 1.382$ を求める。

$\sigma = 1.3$ は $\lambda_2$ に近いため、収束比は

$$r = \left|\frac{\lambda_2 - \sigma}{\lambda_1 - \sigma}\right| = \left|\frac{1.382 - 1.3}{3.618 - 1.3}\right| = \frac{0.082}{2.318} \approx 0.035$$

と非常に小さくなり、2〜3回の反復で高精度に収束する。べき乗法の $r = 0.382$ と比較すると、シフトの選択による収束加速の効果は劇的である。

8.3 レイリー商反復法

例3: レイリー商反復法の3次収束

同じ行列 $A$ に初期ベクトル $\mathbf{x}^{(0)} = (1/\sqrt{2},\ 1/\sqrt{2})^T$ でレイリー商反復を適用する。

$k$	$\sigma^{(k)}$	$\|\sigma^{(k)} - \lambda_1\|$
0	$3.5000$	$1.2 \times 10^{-1}$
1	$3.6176$	$3.9 \times 10^{-4}$
2	$3.618033988738$	$1.2 \times 10^{-11}$
3	$3.618033988749895$	$\approx 10^{-16}$

3次収束では、正確な桁数が毎反復で約3倍に増える。実際、正確な桁数は $1 \to 3 \to 11 \to 16$（機械精度）と推移しており、3回の反復で機械精度に到達している。

9. 収束の可視化

9.1 べき乗法の収束グラフ

図1: $A = \begin{pmatrix}2&1\\1&3\end{pmatrix}$ に対するべき乗法の収束。反復ごとにレイリー商の誤差が $|\lambda_2/\lambda_1|^2 \approx 0.146$ 倍に減少する。

9.2 3手法の収束比較

図2: 3手法の収束速度の比較（$\mathbf{x}_0 = (1,1)^T/\sqrt{2}$、シフト $\sigma = 3.5$）。べき乗法（青）は線形収束、シフト付き逆反復（緑）は高速な線形収束、レイリー商反復法（赤）は3次収束を示す。3手法とも同じ初期ベクトルから出発し、k=1 ではシフト逆反復とレイリー商反復が同一点に到達する（初回のシフトが一致するため）。

10. 実装上の注意

10.1 初期ベクトルの選択

ランダムベクトルを初期値に使えば、$c_1 = 0$ となることはまずない。実用上はランダムな初期ベクトルで十分であるが、何らかの事前情報がある場合はそれを活用することで反復回数を減らせる。

10.2 収束判定

以下のいずれかを使用する。

固有値の変化: $|\mu^{(k)} - \mu^{(k-1)}| < \varepsilon$
残差: $\|A\mathbf{x}^{(k)} - \mu^{(k)}\mathbf{x}^{(k)}\| < \varepsilon$（残差ノルム）

残差ベースの判定はより信頼性が高い。なお、固有値が非常に小さい場合は絶対誤差では不十分なので、実装では

$$\frac{|\mu^{(k)} - \mu^{(k-1)}|}{\max(1,\, |\mu^{(k)}|)} < \varepsilon$$

のように相対誤差を用いることが多い。

10.3 べき乗法が適用できない場合

収束しないケース

$|\lambda_1| = |\lambda_2|$: 例えば $\lambda_1 = -\lambda_2$（反対符号）や $\lambda_1 = \overline{\lambda_2}$（複素共役）の場合、べき乗法は振動して収束しない。
$A$ が対角化不可能: Jordan ブロックが存在する場合、収束は保証されない。

これらの場合は QR法（第5章）や Lanczos 法を使用する。

10.4 計算量の比較

手法	前処理	1反復あたり	収束速度
べき乗法	なし	$O(N^2)$	線形 $\|\lambda_2/\lambda_1\|^k$
逆べき乗法	LU分解 $O(N^3)$	$O(N^2)$	線形 $\|\lambda_{N-1}/\lambda_N\|^k$
シフト付き逆反復	LU分解 $O(N^3)$	$O(N^2)$	線形（$\sigma$ 依存）
レイリー商反復	なし	$O(N^3)$	2次 / 3次（対称）

10.5 他の手法との関係

QR法: 全固有値を同時に求める。暗黙的にべき乗法の並列実行と見なせる。
Lanczos法: 大規模疎行列に対する固有値計算。べき乗法のKrylov部分空間への拡張。
LAPACK の戦略: QR法で固有値を求め、シフト付き逆反復で固有ベクトルを計算する。

11. まとめ

べき乗法は絶対値最大の固有値を求める最も基本的な手法である。収束速度は $O(|\lambda_2/\lambda_1|^k)$ で、$|\lambda_2/\lambda_1|$ が1に近いと遅い。
逆べき乗法は $A^{-1}$ にべき乗法を適用し、絶対値最小の固有値を求める。
シフト付き逆べき乗法は $(A - \sigma I)^{-1}$ にべき乗法を適用し、$\sigma$ に最も近い固有値を求める。LAPACK の固有ベクトル計算の標準手法である。
レイリー商反復法はシフトをレイリー商で動的更新し、対称行列で3次収束を達成する。

参考文献

二宮, 吉田, 長谷川, 泰野, 杉浦, 櫻井, 細田「数値計算のわざ」pp.87–89
戸川隼人「マトリクスの数値計算」pp.28–29
磯田和男, 大野豊「数値計算ハンドブック」pp.89–99
戸川隼人「科学技術計算ハンドブック」pp.332–338
山本哲朗「数値解析入門増訂版」pp.101–102
Wikipedia: 冪乗法