レイリー商とは？固有値を最大化問題として捉える

「固有値を求める」と聞くと、多くの人は特性方程式 $\det(\bm{A} – \lambda \bm{I}) = 0$ を思い浮かべます。これは定義に忠実で、$2 \times 2$ や $3 \times 3$ なら手で解けます。しかし $n = 1000$ の行列で 1000 次方程式を解こうとする人はいません。実際の数値計算ライブラリは、特性多項式を一度も作らずに固有値を求めています。では何をしているのでしょうか。

答えのひとつが「固有値を、方程式の根ではなく最適化問題の最適値として捉え直す」という発想の転換です。対称行列 $\bm{A}$ に対して

$$ R(\bm{x}) = \frac{\bm{x}^\top \bm{A} \bm{x}}{\bm{x}^\top \bm{x}} $$

という量（レイリー商）を考えると、この関数を最大化した値がそのまま最大固有値になります。方程式を解く問題が、丘の頂上を探す問題に化けるのです。しかも「探す」ほうは勾配法でも、反復法でも、幾何的な直感でも攻められます。

左の図では、単位円（灰色）が行列によって細長い楕円に写されています。右上がりの向きは 6 倍に伸びる一方、それと直交する向きは長さが変わりません。行列は「一様に拡大する」のではなく、方向ごとに伸ばし方が違うのです。右の図はその「方向あたりの伸び率」を中心からの距離で表したもので、ひょうたん型のふくらみの先端（赤い点）が最大固有値、くびれた部分（緑の点）が最小固有値に対応します。固有値を求めるとは、この曲線の山頂と谷底を探すことに他なりません。

この見方は理論の飾りではなく、実用の現場で毎日使われています。たとえば主成分分析（PCA）で「データの分散が最大になる方向」を探す作業は、共分散行列のレイリー商を最大化する作業そのものです。構造物の振動解析で「この橋が最も揺れやすい振動数はいくつか」を問うとき、答えは剛性行列と質量行列に関する一般化レイリー商の最小値として現れます。他にも、行列の 2 ノルム、グラフを 2 分割する最適カットの緩和、量子力学の変分法による基底状態エネルギーの推定 — いずれもレイリー商が主役です。

本記事の内容

• レイリー商の直感的な意味と定義、スケール不変性
• 勾配ゼロの条件から $\bm{A}\bm{x} = \lambda \bm{x}$ が出ることの導出（省略なし）
• 固有基底で展開してレイリー商が固有値の凸結合になることの証明
• 第 $k$ 固有値を与えるクーラン・フィッシャーのミニマックス原理
• PCA の「既存成分に直交する方向で分散最大」との完全な対応
• レイリー商反復が 3 次収束することの数値確認

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 固有値分解の理論と応用を完全解説する
• 【対称行列の性質まとめ】定義や固有値、行列式や逆行列を解説
• 【線形代数】行列とベクトルの2次形式や微分の公式を解説

特に「実対称行列は必ず実固有値を持ち、固有ベクトルを正規直交基底に取れる」というスペクトル定理は、本記事の証明の骨格そのものです。ここが曖昧なら、先に対称行列の記事に戻ることをおすすめします。

レイリー商とは — 「方向あたりの伸び率」を測る物差し

バネのついた台車を想像してください。台車を右に 1 cm 引っ張るのと、斜め 45 度に 1 cm 引っ張るのとでは、蓄えられるエネルギーが違います。ばね定数が方向によって違うからです。同じ「1 cm」という移動量でも、方向によってエネルギーの大きさが変わる — この「方向ごとの硬さ」を測りたい、というのがレイリー商の出発点です。

一般に、対称行列 $\bm{A}$ が定める二次形式 $\bm{x}^\top \bm{A} \bm{x}$ は「$\bm{x}$ という変位に対するエネルギー」のような量です。ところがこの値は $\bm{x}$ を 2 倍すると 4 倍になってしまいます。長さの効果と方向の効果が混ざっていて、方向だけの性質を取り出せません。そこで長さの効果を割り算で打ち消します。$\bm{x}^\top \bm{x} = \|\bm{x}\|^2$ で割れば、分子も分母も $\bm{x}$ の 2 次なので、スケールの影響がきれいに消えます。

こうして得られるのがレイリー商です。

$$ \begin{equation} R(\bm{x}) = \frac{\bm{x}^\top \bm{A} \bm{x}}{\bm{x}^\top \bm{x}}, \qquad \bm{x} \neq \bm{0} \end{equation} $$

ここで $\bm{A} \in \mathbb{R}^{n \times n}$ は実対称行列（$\bm{A}^\top = \bm{A}$）、$\bm{x} \in \mathbb{R}^n$ はゼロでないベクトルです。名前はイギリスの物理学者レイリー卿（Lord Rayleigh）に由来します。彼が音響学で振動の固有振動数を近似する際にこの量を使ったのが起源です。

スケール不変性 — 見るべきは方向だけ

まず確認すべき性質は、任意の実数 $c \neq 0$ について

$$ R(c\bm{x}) = \frac{(c\bm{x})^\top \bm{A} (c\bm{x})}{(c\bm{x})^\top (c\bm{x})} = \frac{c^2 \, \bm{x}^\top \bm{A} \bm{x}}{c^2 \, \bm{x}^\top \bm{x}} = R(\bm{x}) $$

が成り立つことです。分子・分母がともに $c^2$ 倍されて相殺します。つまりレイリー商は $\bm{x}$ の長さを一切見ておらず、$\bm{x}$ が指す「方向」だけの関数です。

右の棒グラフを見ると、長さを 0.5 倍・1 倍・2 倍と変えるにつれて分子 $\bm{x}^\top\bm{A}\bm{x}$ は $1.5 \to 6.0 \to 24.0$ と 4 倍ずつ増え、分母 $\bm{x}^\top\bm{x}$ も $0.25 \to 1.0 \to 4.0$ と同じ 4 倍ずつ増えています。両者の増え方が完全に揃っているため、赤い棒（商）はどの長さでも $6.0$ で動きません。分子だけを見ていては「方向の性質」が読み取れないのに対し、割り算をひとつ挟むだけで長さの情報がきれいに落ちる、というのがレイリー商の設計思想です。

この事実は実務上とても便利です。$R$ の最大化を考えるとき、$\bm{x}$ の動く範囲を $\mathbb{R}^n$ 全体（有界でない）ではなく、単位球面 $\|\bm{x}\| = 1$（有界閉集合）に制限してよいことになります。連続関数は有界閉集合上で必ず最大値と最小値を取る（ワイエルシュトラスの最大値定理）ので、この制限のおかげで「最大値が存在するか？」という厄介な問いが一瞬で片付きます。

単位球面に制限すると分母が 1 になるので、レイリー商の最大化は次と同じです。

$$ \begin{equation} \max_{\|\bm{x}\| = 1} \bm{x}^\top \bm{A} \bm{x} \end{equation} $$

「単位球の上を歩き回って、二次形式の値が一番高い場所を探す」— これがレイリー商最大化の絵です。$n = 2$ なら単位円の上を一周するだけ、$n = 3$ なら球面上を動き回るだけ。とても素朴な問題に見えますが、この頂上こそが最大固有値です。

左は二次形式 $\bm{x}^\top\bm{A}\bm{x}$ の曲面で、原点から離れるほど値が跳ね上がるお椀型をしています。この曲面をそのまま最大化しようとしても無限大に発散してしまいますが、単位円という「輪」に沿って切り取ると（赤い曲線）、値は有限の範囲に閉じ込められます。右はその赤い曲線を角度に沿って平らに伸ばしたもので、上端が $\lambda_1 = 6$、下端が $\lambda_2 = 1$ の水平線にぴたりと接しています。青く塗った帯がレイリー商の値域であり、この帯の上下端が固有値だ、というのがこれから証明する主張です。

まず固有ベクトルを代入してみる

証明に入る前に、答えの当たりをつけましょう。$\bm{v}$ が $\bm{A}$ の固有ベクトルで、固有値が $\lambda$ だとします。つまり $\bm{A}\bm{v} = \lambda \bm{v}$ です。これをレイリー商に入れると

$$ R(\bm{v}) = \frac{\bm{v}^\top \bm{A} \bm{v}}{\bm{v}^\top \bm{v}} = \frac{\bm{v}^\top (\lambda \bm{v})}{\bm{v}^\top \bm{v}} = \lambda \frac{\bm{v}^\top \bm{v}}{\bm{v}^\top \bm{v}} = \lambda $$

となります。固有ベクトルを入れると、レイリー商はぴったり対応する固有値を返すわけです。しかもスケール不変性のおかげで、$\bm{v}$ を正規化していなくても構いません。

ここまでで分かったのは「レイリー商の値域には全ての固有値が含まれる」ことだけです。しかし本当に知りたいのはその逆、つまり「レイリー商が取りうる値は固有値の範囲をはみ出さないのか」「最大値は最大固有値ちょうどなのか」です。次の 2 つのセクションで、これを 2 通りの方法で証明します。ひとつは微分（勾配）を使う解析的な方法、もうひとつは固有基底で展開する代数的な方法です。前者は「なぜ固有方程式が出てくるのか」を説明し、後者は「なぜ最大値がそこで止まるのか」を説明します。

勾配ゼロの条件から固有方程式が現れる

レイリー商を最大化したいなら、まず微分して勾配がゼロになる点（臨界点）を探すのが定石です。ここで驚くべきことが起きます。臨界点の条件式が、そのまま固有方程式 $\bm{A}\bm{x} = \lambda\bm{x}$ になるのです。固有値問題を一度も持ち出していないのに、最適化の一階条件から勝手に湧いて出てきます。この事実こそが「固有値 = 最適化問題の解」という見方の核心です。

準備: 二次形式とノルムの勾配

$\bm{A}$ が対称のとき、二次形式 $f(\bm{x}) = \bm{x}^\top \bm{A} \bm{x}$ の勾配は

$$ \nabla f(\bm{x}) = (\bm{A} + \bm{A}^\top)\bm{x} = 2\bm{A}\bm{x} $$

です。成分で確かめておきましょう。$f = \sum_{i}\sum_{j} a_{ij} x_i x_j$ を $x_k$ で偏微分すると、$x_k$ を含む項は「$i=k$ の項」と「$j=k$ の項」の 2 種類あるので

$$ \frac{\partial f}{\partial x_k} = \sum_{j} a_{kj} x_j + \sum_{i} a_{ik} x_i = (\bm{A}\bm{x})_k + (\bm{A}^\top \bm{x})_k $$

となり、対称性 $\bm{A}^\top = \bm{A}$ から $2(\bm{A}\bm{x})_k$ です。同様に $g(\bm{x}) = \bm{x}^\top \bm{x}$ については $\bm{A} = \bm{I}$ の特別な場合なので $\nabla g = 2\bm{x}$ です。

商の微分でレイリー商の勾配を出す

$R = f/g$ なので、商の微分公式（ベクトル版でも形は同じ）を使います。

$$ \nabla R(\bm{x}) = \frac{g(\bm{x}) \nabla f(\bm{x}) – f(\bm{x}) \nabla g(\bm{x})}{g(\bm{x})^2} $$

上で求めた $\nabla f = 2\bm{A}\bm{x}$、$\nabla g = 2\bm{x}$ を代入すると

$$ \nabla R(\bm{x}) = \frac{(\bm{x}^\top \bm{x}) \cdot 2\bm{A}\bm{x} – (\bm{x}^\top \bm{A} \bm{x}) \cdot 2\bm{x}}{(\bm{x}^\top \bm{x})^2} $$

分子と分母から共通因子 $\bm{x}^\top \bm{x}$ を 1 つずつ約分し、$\bm{x}^\top \bm{A}\bm{x}/\bm{x}^\top\bm{x} = R(\bm{x})$ であることを使って第 2 項をまとめると、非常に見通しの良い形になります。

$$ \begin{equation} \nabla R(\bm{x}) = \frac{2}{\bm{x}^\top \bm{x}} \left( \bm{A}\bm{x} – R(\bm{x}) \, \bm{x} \right) \end{equation} $$

この式は一度見たら忘れられません。括弧の中は「$\bm{A}\bm{x}$ と、$\bm{x}$ を $R(\bm{x})$ 倍したもの、の差」です。レイリー商の勾配は、固有方程式の残差そのものなのです。

臨界点 ⟺ 固有ベクトル

$\bm{x} \neq \bm{0}$ なので $\bm{x}^\top\bm{x} > 0$ であり、$\nabla R(\bm{x}) = \bm{0}$ となるのは括弧の中がゼロのとき、つまり

$$ \begin{equation} \bm{A}\bm{x} = R(\bm{x}) \, \bm{x} \end{equation} $$

のときに限ります。これは「$\bm{x}$ が固有値 $R(\bm{x})$ に属する固有ベクトルである」と言っているのと同じです。逆向きも成り立ちます。$\bm{x}$ が固有値 $\lambda$ の固有ベクトルなら、前節で見たとおり $R(\bm{x}) = \lambda$ なので、$\bm{A}\bm{x} – R(\bm{x})\bm{x} = \lambda\bm{x} – \lambda\bm{x} = \bm{0}$ となり勾配はゼロです。

まとめると次の同値が得られました。

$\bm{x} \neq \bm{0}$ がレイリー商 $R$ の臨界点である $\iff$ $\bm{x}$ が $\bm{A}$ の固有ベクトルである。そのとき $R(\bm{x})$ はその固有値に等しい。

言い換えると、単位球面上でレイリー商という「高さ関数」を描いたとき、その峠・山頂・谷底といった平らな点は、すべて固有ベクトルの方向にしか存在しません。固有値問題は「高さ関数の臨界点を全部探せ」という問題と完全に等価だったわけです。

左図の青い矢印が各点での勾配 $\nabla R$ です。矢印は必ず円に接する向きを向いており（$\bm{x}^\top(\bm{A}\bm{x} – R\bm{x}) = 0$ なので、勾配は動径方向の成分を持ちません）、赤い点（山頂）へ向かって流れ込み、緑の点（谷底）から湧き出しています。そして赤・緑の 4 点では矢印が完全に消えています。右図はこれを角度の関数として重ねたもので、勾配の大きさ $\|\nabla R\|$ がゼロに落ちる角度（$26.6^\circ, 116.6^\circ, 206.6^\circ, 296.6^\circ$）が、レイリー商 $R(\theta)$ の山頂・谷底の角度と 1 つの例外もなく一致しています。臨界点は固有ベクトル方向にしかない、という同値関係が目で確認できます。

ラグランジュ未定乗数法で見ても同じ

同じ結論は、制約付き最適化として書いても得られます。$\|\bm{x}\|^2 = 1$ という制約の下で $\bm{x}^\top\bm{A}\bm{x}$ を最大化する問題のラグランジアンは

$$ L(\bm{x}, \lambda) = \bm{x}^\top \bm{A} \bm{x} – \lambda (\bm{x}^\top \bm{x} – 1) $$

です。$\bm{x}$ で偏微分してゼロと置くと、上で計算した勾配を使って

$$ \nabla_{\bm{x}} L = 2\bm{A}\bm{x} – 2\lambda \bm{x} = \bm{0} \quad \Longrightarrow \quad \bm{A}\bm{x} = \lambda\bm{x} $$

となり、やはり固有方程式が出ます。しかも未定乗数 $\lambda$ がそのまま固有値の役を担っているのが面白いところです。制約 $\bm{x}^\top\bm{x} = 1$ を左から $\bm{x}^\top$ を掛けて使えば $\lambda = \bm{x}^\top \bm{A} \bm{x} = R(\bm{x})$ となり、乗数の正体は最適値そのものだと分かります。ラグランジュ乗数が「制約を 1 単位緩めたときの目的関数の増分（影の価格）」であることを思い出すと、固有値は「単位球の半径を少し伸ばしたときにエネルギーがどれだけ増えるか」を表す量だと解釈できます。

ここまでで「候補は固有ベクトルだけ」と分かりました。しかし候補が $n$ 個ある中で、どれが最大でどれが最小かはまだ決まっていません。$n=3$ の球面上に山頂・鞍点・谷底が 1 つずつあったとして、どれが最大固有値に対応するのか。それを一撃で解決するのが、次に見る固有基底での展開です。

固有基底で展開する — レイリー商は固有値の凸結合

スペクトル定理から、実対称行列 $\bm{A}$ には正規直交な固有ベクトルの組 $\bm{v}_1, \dots, \bm{v}_n$ が存在します。対応する固有値を大きい順に並べて

$$ \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n $$

とします。$\{\bm{v}_i\}$ は $\mathbb{R}^n$ の正規直交基底なので、任意の $\bm{x}$ をこの基底で展開できます。

$$ \bm{x} = \sum_{i=1}^{n} c_i \bm{v}_i, \qquad c_i = \bm{v}_i^\top \bm{x} $$

この座標 $c_i$ は「$\bm{x}$ が第 $i$ 固有方向にどれだけ寄っているか」を表す量です。ここからの計算は、正規直交性 $\bm{v}_i^\top \bm{v}_j = \delta_{ij}$ が全ての交差項を消してくれるので驚くほど簡単に進みます。

分母の計算

まず分母から片付けます。展開を代入して 2 重和にすると

$$ \bm{x}^\top \bm{x} = \left(\sum_{i} c_i \bm{v}_i\right)^\top \left(\sum_{j} c_j \bm{v}_j\right) = \sum_{i}\sum_{j} c_i c_j \, \bm{v}_i^\top \bm{v}_j $$

ここで $\bm{v}_i^\top \bm{v}_j$ は $i = j$ のとき 1、それ以外は 0 なので、2 重和のうち対角項だけが生き残ります。

$$ \bm{x}^\top \bm{x} = \sum_{i=1}^{n} c_i^2 $$

これはパーセバルの等式（正規直交基底では座標の二乗和がノルムの二乗に等しい）そのものです。

分子の計算

分子も同じ要領ですが、途中で固有方程式を使うのがポイントです。まず $\bm{A}\bm{x}$ を計算します。$\bm{A}$ は線形なので和の外に出せて、$\bm{A}\bm{v}_j = \lambda_j \bm{v}_j$ を代入すると

$$ \bm{A}\bm{x} = \sum_{j} c_j \bm{A} \bm{v}_j = \sum_{j} c_j \lambda_j \bm{v}_j $$

固有ベクトルに $\bm{A}$ を掛ける操作が、単なるスカラー倍に化けたのが効いています。これに左から $\bm{x}^\top = \sum_i c_i \bm{v}_i^\top$ を掛けると

$$ \bm{x}^\top \bm{A} \bm{x} = \sum_{i}\sum_{j} c_i c_j \lambda_j \, \bm{v}_i^\top \bm{v}_j $$

となり、再び正規直交性で $i = j$ の項だけが残ります。

$$ \bm{x}^\top \bm{A} \bm{x} = \sum_{i=1}^{n} \lambda_i c_i^2 $$

レイリー商 = 重み付き平均

以上を合わせると、レイリー商は次の形になります。

$$ \begin{equation} R(\bm{x}) = \frac{\sum_{i=1}^{n} \lambda_i c_i^2}{\sum_{i=1}^{n} c_i^2} \end{equation} $$

ここで

$$ w_i = \frac{c_i^2}{\sum_{k} c_k^2} $$

と置くと、$w_i \geq 0$ かつ $\sum_i w_i = 1$ です。すなわち $\bm{w} = (w_1, \dots, w_n)$ は確率分布（重みベクトル）の条件を満たしており、

$$ \begin{equation} R(\bm{x}) = \sum_{i=1}^{n} w_i \lambda_i \end{equation} $$

つまり レイリー商とは、固有値たちの凸結合（重み付き平均）である という結論に到達します。重み $w_i$ は「ベクトル $\bm{x}$ のエネルギーのうち、第 $i$ 固有方向に配分されている割合」です。$\bm{x}$ が完全に $\bm{v}_1$ を向いていれば $\bm{w} = (1, 0, \dots, 0)$ となり $R = \lambda_1$、方向がばらけていれば複数の固有値の平均になります。

固有値を $\{7, 4, 2, 1\}$ とした 4 次元の例です。左の積み上げ棒は各ベクトルのエネルギー配分 $w_i$ で、どの棒も高さの合計がちょうど 1 になっています。$\bm{x} = \bm{v}_1$ なら重みが第 1 固有方向に 100% 集中し、$(\bm{v}_1+\bm{v}_2)/\sqrt{2}$ なら 50% ずつ、ランダムなベクトルなら 4 方向にばらけます。右はその結果得られるレイリー商の位置で、$7.000, 5.500, 5.756, 1.000$ といずれも $\lambda_4 = 1$ と $\lambda_1 = 7$ の帯の内側に収まっています。重みが 1 点に集中したときだけ端点（固有値そのもの）に到達する、という凸結合の性質がそのまま現れています。

最大・最小の決定

凸結合の値は、材料となる数の最大値を超えられず、最小値を下回れません。形式的に書けば、$\lambda_n \leq \lambda_i \leq \lambda_1$ を重み付きで足し上げて

$$ \lambda_n = \sum_i w_i \lambda_n \leq \sum_i w_i \lambda_i \leq \sum_i w_i \lambda_1 = \lambda_1 $$

となります。したがって

$$ \begin{equation} \lambda_n \leq R(\bm{x}) \leq \lambda_1 \qquad (\forall \bm{x} \neq \bm{0}) \end{equation} $$

が示せました。しかも $\bm{x} = \bm{v}_1$ とすれば $R = \lambda_1$、$\bm{x} = \bm{v}_n$ とすれば $R = \lambda_n$ なので、この上下限は達成されます。以上をまとめると、レイリー商の変分特徴づけ（Rayleigh の原理）が得られます。

$$ \begin{equation} \lambda_1 = \max_{\bm{x} \neq \bm{0}} R(\bm{x}) = \max_{\|\bm{x}\|=1} \bm{x}^\top \bm{A} \bm{x}, \qquad \lambda_n = \min_{\bm{x} \neq \bm{0}} R(\bm{x}) = \min_{\|\bm{x}\|=1} \bm{x}^\top \bm{A} \bm{x} \end{equation} $$

等号が成立するのはどんなときでしょうか。$\lambda_1$ が単純固有値（重複度 1）なら、$R(\bm{x}) = \lambda_1$ となるのは $w_1 = 1$、すなわち $\bm{x}$ が $\bm{v}_1$ の定数倍のときだけです。$\lambda_1$ が重複しているなら、その固有空間の中のどのベクトルでも最大値を達成します。この「重複すると最適解が一意でなくなる」性質は、PCA で第 1 主成分が定まらない（分散が縮退している）状況として実務にも顔を出します。

副産物 — 二次形式の符号と条件数

この不等式からいくつかの重要な事実がタダで手に入ります。まず、$\bm{A}$ が正定値（すべての $\bm{x} \neq \bm{0}$ で $\bm{x}^\top\bm{A}\bm{x} > 0$）であることは $\lambda_n > 0$ と同値です。実際、$\bm{x}^\top\bm{A}\bm{x} = R(\bm{x})\|\bm{x}\|^2$ なので、二次形式の符号はレイリー商の符号と同じであり、レイリー商の最小値が $\lambda_n$ だからです。半正定値なら $\lambda_n \geq 0$、不定符号なら $\lambda_n < 0 < \lambda_1$ です。

また対称行列の場合、行列の 2 ノルムは $\|\bm{A}\|_2 = \max_i |\lambda_i|$ となり、これも $|R(\bm{x})| \leq \max_i|\lambda_i|$ から従います。一般の（対称でない）行列 $\bm{M}$ についても、$\bm{M}^\top\bm{M}$ は対称半正定値なので

$$ \|\bm{M}\|_2 = \max_{\|\bm{x}\|=1}\|\bm{M}\bm{x}\| = \sqrt{\max_{\|\bm{x}\|=1} \bm{x}^\top \bm{M}^\top\bm{M}\bm{x}} = \sqrt{\lambda_{\max}(\bm{M}^\top\bm{M})} = \sigma_{\max}(\bm{M}) $$

と、最大特異値に等しいことがレイリー商から直ちに出ます。同様に最小特異値は最小値側から出るので、条件数 $\kappa_2(\bm{M}) = \sigma_{\max}/\sigma_{\min}$ も「レイリー商の値域の広がり」として理解できます。数値計算で条件数が悪い行列が嫌われる理由が、この「方向によって伸び率が極端に違う」という幾何に還元されるわけです。

ここまでの議論は、抽象的な $n$ 次元での話でした。次は $2 \times 2$ の具体例で、レイリー商が単位円の上でどんな形をしているかを完全に手計算して、目で確かめましょう。

具体例: 2×2 行列で単位円上のレイリー商を手計算する

次の対称行列を考えます。

$$ \bm{A} = \begin{pmatrix} 5 & 2 \\ 2 & 2 \end{pmatrix} $$

固有値・固有ベクトルを求める

特性方程式は

$$ \det(\bm{A} – \lambda \bm{I}) = (5-\lambda)(2-\lambda) – 4 = \lambda^2 – 7\lambda + 6 = 0 $$

なので $(\lambda – 6)(\lambda – 1) = 0$ より $\lambda_1 = 6$、$\lambda_2 = 1$ です。トレース $5+2 = 7$ が固有値の和、行列式 $10 – 4 = 6$ が固有値の積になっていることも確認できます。

固有ベクトルは、$\lambda_1 = 6$ のとき $(\bm{A} – 6\bm{I})\bm{v} = \bm{0}$ すなわち $-v_1 + 2v_2 = 0$ から $\bm{v}_1 \propto (2, 1)^\top$、$\lambda_2 = 1$ のとき $4v_1 + 2v_2 = 0$ から $\bm{v}_2 \propto (1, -2)^\top$ です。内積は $2 \cdot 1 + 1 \cdot (-2) = 0$ なので確かに直交しています。正規化すると $\bm{v}_1 = (2,1)^\top/\sqrt{5}$、$\bm{v}_2 = (1,-2)^\top/\sqrt{5}$ です。$\bm{v}_1$ の方向は $x$ 軸から $\arctan(1/2) \approx 26.57^\circ$ の向きになります。

レイリー商を角度の関数として書く

単位円上の点を $\bm{x}(\theta) = (\cos\theta, \sin\theta)^\top$ とパラメータ表示します。分母は 1 なので、レイリー商は二次形式そのものです。

$$ R(\theta) = 5\cos^2\theta + 2 \cdot 2 \cos\theta\sin\theta + 2\sin^2\theta $$

（交差項 $a_{12}x_1x_2$ と $a_{21}x_2x_1$ の 2 つがあるので係数が $2 a_{12} = 4$ になります。）ここから三角関数の倍角公式で整理していきます。まず $2\cos\theta\sin\theta = \sin 2\theta$ を使うと第 2 項は $2\sin 2\theta$ です。第 1・第 3 項は $\sin^2\theta = 1 – \cos^2\theta$ を代入してまとめると

$$ 5\cos^2\theta + 2\sin^2\theta = 2(\cos^2\theta + \sin^2\theta) + 3\cos^2\theta = 2 + 3\cos^2\theta $$

さらに半角公式 $\cos^2\theta = (1 + \cos 2\theta)/2$ を代入すると

$$ 2 + \frac{3}{2}(1 + \cos 2\theta) = \frac{7}{2} + \frac{3}{2}\cos 2\theta $$

を得ます。したがって

$$ R(\theta) = \frac{7}{2} + \frac{3}{2}\cos 2\theta + 2 \sin 2\theta $$

最後に、$a\cos\phi + b\sin\phi = \sqrt{a^2+b^2}\cos(\phi – \varphi)$（$\tan\varphi = b/a$）という合成公式を $\phi = 2\theta$、$a = 3/2$、$b = 2$ に適用します。振幅は $\sqrt{(3/2)^2 + 2^2} = \sqrt{9/4 + 4} = \sqrt{25/4} = 5/2$ なので

$$ \begin{equation} R(\theta) = \frac{7}{2} + \frac{5}{2}\cos(2\theta – \varphi), \qquad \tan\varphi = \frac{2}{3/2} = \frac{4}{3} \end{equation} $$

という、極めて明快な形になりました。

左は分解の内訳です。灰色の水平線（定数 $7/2$）に、オレンジの $\frac{3}{2}\cos 2\theta$ と水色の $2\sin 2\theta$ を足すと、青い太線 $R(\theta)$ が再現されます。振動する 2 項はどちらも周期 $180^\circ$ なので、和も同じ周期になります。右は合成後の姿で、中心線 $7/2$ がトレースの半分、紫の両矢印で示した振幅 $5/2$ が固有値の差の半分に一致し、その結果として最大値が $6 = \lambda_1$、最小値が $1 = \lambda_2$ にぴたりと収まっています。$\theta^* = 26.57^\circ$ の赤い縦線が第 1 固有ベクトルの向きです。

読み取れること

この式から一気に多くのことが分かります。

第一に、$R(\theta)$ は $\theta$ の周期 $\pi$ の単純な正弦波です。$2\theta$ で振動するので、単位円を一周する間に山と谷がそれぞれ 2 回来ます。これは $R(-\bm{x}) = R(\bm{x})$（原点対称な 2 点は同じ値）という性質の現れです。

第二に、中心の高さ $7/2$ はトレースの半分、つまり固有値の平均 $(\lambda_1 + \lambda_2)/2 = (6+1)/2 = 3.5$ です。振幅 $5/2$ は固有値の差の半分 $(\lambda_1 – \lambda_2)/2 = 2.5$ です。したがって最大値は $3.5 + 2.5 = 6 = \lambda_1$、最小値は $3.5 – 2.5 = 1 = \lambda_2$ となり、変分特徴づけが具体例で確認できました。

第三に、最大値を取る角度は $2\theta – \varphi = 0$ すなわち $\theta = \varphi/2$ です。$\tan\varphi = 4/3$ より $\varphi = \arctan(4/3) \approx 53.13^\circ$ なので $\theta^* \approx 26.57^\circ$。これは先ほど求めた固有ベクトル $\bm{v}_1 = (2,1)^\top/\sqrt 5$ の方向と完全に一致します。最小値は $2\theta – \varphi = \pi$ すなわち $\theta \approx 116.57^\circ$ で、これは $\theta^*$ より $90^\circ$ ずれた方向、つまり $\bm{v}_2$ の方向です。

第四に、山と谷が必ず直交する理由もこの式が説明します。$\cos(2\theta – \varphi)$ の最大と最小は $2\theta$ が $\pi$ だけ離れたところ、つまり $\theta$ が $\pi/2$ だけ離れたところで起きます。角度が $2\theta$ で回るという事実が、そのまま「対称行列の異なる固有値に属する固有ベクトルは直交する」という定理の $2$ 次元版になっているのです。

数値でも確かめてみましょう。単位円を 20 万点に分割してレイリー商を評価すると、最大値 $5.9999999994$（角度 $26.5644^\circ$）、最小値 $1.0000000006$（角度 $116.5644^\circ$）が得られます。理論値 $6$ と $1$、角度 $26.565^\circ$ と $116.565^\circ$ に、グリッド分解能の範囲でぴたりと一致しています。

$2$ 次元では「最大」と「最小」しかないので、これで全部です。しかし $n$ 次元では $\lambda_2, \lambda_3, \dots$ という中間の固有値があります。これらも変分的に表せるのでしょうか。答えはイエスで、その表現がクーラン・フィッシャーのミニマックス原理です。

第k固有値へ — クーラン・フィッシャーのミニマックス原理

$\lambda_1$ は「単位球全体での最大値」でした。では $\lambda_2$ は何でしょうか。素朴に考えれば「1 番高い山を除いた後の、2 番目に高い山」です。この「除く」を数学的に書くと「$\bm{v}_1$ に直交する平面に限定する」となります。

素朴版: 固有ベクトルへの直交制約

$\bm{x}$ を $\bm{v}_1$ に直交させると、展開係数のうち $c_1 = \bm{v}_1^\top\bm{x} = 0$ になります。すると凸結合の重み $w_1$ がゼロになり

$$ R(\bm{x}) = \sum_{i=2}^{n} w_i \lambda_i \leq \lambda_2 $$

が成り立ちます。$\lambda_2$ 以下の固有値だけの凸結合になるので、その最大は $\lambda_2$ です。しかも $\bm{x} = \bm{v}_2$ は制約を満たし（$\bm{v}_2 \perp \bm{v}_1$）、そのとき $R = \lambda_2$ です。同じ議論を繰り返せば、一般に

$$ \begin{equation} \lambda_k = \max \left\{ R(\bm{x}) \; : \; \bm{x} \neq \bm{0}, \; \bm{x} \perp \bm{v}_1, \dots, \bm{x} \perp \bm{v}_{k-1} \right\} \end{equation} $$

が得られます。これを逐次的な変分特徴づけと呼びます。

しかしこの表現には弱点があります。$\lambda_k$ を特徴づけるのに $\bm{v}_1, \dots, \bm{v}_{k-1}$ という固有ベクトルを既に知っている必要があるのです。「固有値を知らずに固有値を特徴づけたい」という当初の目的からすると、循環している感じが否めません。理論的な応用（たとえば 2 つの行列の固有値を比較する）でも、それぞれの行列の固有ベクトルが違うので、この形では扱いにくい。

クーラン・フィッシャー: 固有ベクトルを使わない表現

この弱点を解消するのがクーラン・フィッシャー（Courant–Fischer）のミニマックス原理です。固有ベクトルを一切使わず、部分空間の次元だけで第 $k$ 固有値を特徴づけます。

$$ \begin{equation} \lambda_k = \max_{\substack{\mathcal{S} \subseteq \mathbb{R}^n \\ \dim \mathcal{S} = k}} \; \min_{\substack{\bm{x} \in \mathcal{S} \\ \bm{x} \neq \bm{0}}} R(\bm{x}) \end{equation} $$

同じことを裏返した表現もあります。

$$ \begin{equation} \lambda_k = \min_{\substack{\mathcal{S} \subseteq \mathbb{R}^n \\ \dim \mathcal{S} = n-k+1}} \; \max_{\substack{\bm{x} \in \mathcal{S} \\ \bm{x} \neq \bm{0}}} R(\bm{x}) \end{equation} $$

前者（max-min 形）を読み下すと、「$k$ 次元の部分空間を 1 つ選ぶ。その部分空間の中でレイリー商が一番低くなる方向を見て、その値を部分空間のスコアとする。あらゆる $k$ 次元部分空間の中でスコアが最も高いものを探すと、そのスコアが $\lambda_k$ になる」という意味です。「一番弱いところで評価して、その中で最強のチームを選ぶ」というミニマックス的な構造です。

証明

（1）$\lambda_k$ が達成可能であること（$\geq$ 方向の下界）

$\mathcal{S}_0 = \mathrm{span}\{\bm{v}_1, \dots, \bm{v}_k\}$ という特別な $k$ 次元部分空間を取ります。この中の任意の $\bm{x} \neq \bm{0}$ は $\bm{x} = \sum_{i=1}^{k} c_i \bm{v}_i$ と書けるので、凸結合表現の重みは $i \leq k$ にしか乗りません。よって

$$ R(\bm{x}) = \sum_{i=1}^{k} w_i \lambda_i \geq \lambda_k \sum_{i=1}^k w_i = \lambda_k $$

となり、$\mathcal{S}_0$ 上の最小値は $\lambda_k$ 以上です。実際 $\bm{x} = \bm{v}_k$ で等号なので、$\min_{\bm{x}\in\mathcal{S}_0} R(\bm{x}) = \lambda_k$ ちょうどです。したがって「あらゆる $k$ 次元部分空間の中で最大」は少なくとも $\lambda_k$ 以上です。

（2）$\lambda_k$ を超えられないこと（$\leq$ 方向の上界）

任意の $k$ 次元部分空間 $\mathcal{S}$ を取ります。ここで次元の数え上げを使います。$\mathcal{T} = \mathrm{span}\{\bm{v}_k, \bm{v}_{k+1}, \dots, \bm{v}_n\}$ は $n – k + 1$ 次元です。2 つの部分空間の次元の和は

$$ \dim \mathcal{S} + \dim \mathcal{T} = k + (n – k + 1) = n + 1 > n $$

なので、これらは必ずゼロでない共通部分を持ちます（もし $\mathcal{S} \cap \mathcal{T} = \{\bm{0}\}$ なら $\dim(\mathcal{S} + \mathcal{T}) = n+1 > n$ となって矛盾するからです）。そこで $\bm{y} \in \mathcal{S} \cap \mathcal{T}$、$\bm{y} \neq \bm{0}$ を取ります。$\bm{y} \in \mathcal{T}$ なので、その展開には $\bm{v}_k$ 以降しか現れません。ゆえに

$$ R(\bm{y}) = \sum_{i=k}^{n} w_i \lambda_i \leq \lambda_k $$

です。この $\bm{y}$ は $\mathcal{S}$ の要素でもあるので、$\min_{\bm{x} \in \mathcal{S}} R(\bm{x}) \leq R(\bm{y}) \leq \lambda_k$ となります。つまりどんな $k$ 次元部分空間を持ってきてもスコアは $\lambda_k$ を超えられません。

（1）と（2）を合わせて、最大値がちょうど $\lambda_k$ で達成されることが示せました。min-max 形も、$\bm{A}$ を $-\bm{A}$ に置き換えて固有値の順序を反転させれば同じ議論で出ます。

数値で確かめる

固有値が $\{7, 4, 2, 1\}$（降順）となる $4 \times 4$ 対称行列を作り、ランダムな $k$ 次元部分空間を 5 万個生成して「部分空間上の最小レイリー商」の最大値を探すと、次の結果が得られます（部分空間の最小レイリー商は、正規直交基底 $\bm{U}$ に対する $\bm{U}^\top\bm{A}\bm{U}$ の最小固有値として計算できます）。

import numpy as np

rng = np.random.default_rng(0)
Q, _ = np.linalg.qr(rng.standard_normal((4, 4)))
lams = np.array([7.0, 4.0, 2.0, 1.0])
A = Q @ np.diag(lams) @ Q.T
A = (A + A.T) / 2

rs = np.random.default_rng(1)
for k in range(1, 5):
    best = -np.inf
    for _ in range(50000):
        U, _ = np.linalg.qr(rs.standard_normal((4, k)))   # ランダムなk次元部分空間
        score = np.linalg.eigvalsh(U.T @ A @ U)[0]        # 部分空間上の最小レイリー商
        best = max(best, score)
    print("k=%d  最大スコア=%.7f  理論値=%g  差=%.1e"
          % (k, best, lams[k - 1], lams[k - 1] - best))

$k$ が小さいほど理論値に届きにくいことが読み取れます。$k=1$（直線）の場合、最適な向きは 4 次元空間の中でただ 1 本しかなく、$k$ 次元部分空間の自由度（グラスマン多様体の次元）は $k(n-k)$ なので $k=1$ では 3 自由度をすべて当てにいく必要があり、ランダムでは詰めきれません。逆に $k$ が大きくなるほど、当てるべき自由度が減ります。とくに $k=3$ が $10^{-9}$ まで届くのには明確な理由があります。$k = n-1 = 3$ の部分空間は法線ベクトル $\bm{u}$ 1 本で決まり、$\bm{u} = \sum_i a_i \bm{v}_i$ と展開したとき、スコアの $\lambda_3$ からの不足量は（$\bm{u}$ が主に $\bm{v}_4$ 方向を向いていれば）

$$ \lambda_3 – \text{スコア} \;\simeq\; \frac{a_3^2}{a_4^2 – a_1^2/5 – a_2^2/2} $$

と、$a_3$ の 2 乗に比例します。つまり「$\bm{u}$ の $\bm{v}_3$ 成分をゼロに近づける」という実質 1 次元の条件を満たせばよく、5 万回引けば $|a_3| \sim 10^{-5}$ 程度は容易に出るので、不足量は $10^{-10}$ 台まで小さくなります。$k=4$ に至っては部分空間が空間全体しかなく、最小値は必ず $\lambda_4 = 1$ になるので完全一致するのは当然です。

左のヒストグラムは、各 $k$ についてランダムな部分空間 5 万個のスコア分布です。どの分布も右端が対応する $\lambda_k$（同色の破線）でぴたりと切れており、$\lambda_k$ を超えるスコアは 1 つも出現しません。これが証明の（2）「どんな $k$ 次元部分空間でもスコアは $\lambda_k$ を超えられない」の実測版です。右のグラフは試行回数に対する最良スコアの推移で、$k=3, 4$ は最初の数試行で理論値に到達する一方、$k=1$ は 5 万回かけても $7$ にわずかに届いていません。上限が破られないことと、その上限が達成可能であること — ミニマックス原理の 2 つの主張が 1 枚の図に収まっています。

何が嬉しいのか — 摂動と比較の道具になる

クーラン・フィッシャー表現の威力は、2 つの行列の固有値を比較できる点にあります。たとえば $\bm{A} \preceq \bm{B}$（$\bm{B} – \bm{A}$ が半正定値）なら、すべての $\bm{x}$ で $R_{\bm{A}}(\bm{x}) \leq R_{\bm{B}}(\bm{x})$ なので、ミニマックスの中身が各点で大小関係を持ち、そのまま $\lambda_k(\bm{A}) \leq \lambda_k(\bm{B})$ が全ての $k$ について従います。固有ベクトルが両者で全く違っていても構いません。この単調性から、ワイルの不等式

$$ \lambda_k(\bm{A} + \bm{E}) \leq \lambda_k(\bm{A}) + \lambda_1(\bm{E}) $$

のような摂動評価が導かれ、「行列を少し揺らしたとき固有値がどれだけ動くか」を評価できます。センサ較正の誤差が共分散行列の固有値をどれだけ動かすか、といった実問題の見積もりに直結する道具です。

さて、この「既に選んだ方向に直交する部分空間の中で、レイリー商を最大化する」という手続きは、実はデータ解析の教科書でおなじみの操作と同じものです。次はそれを見ます。

PCAとの対応 — 「分散最大の方向」はレイリー商最大化そのもの

主成分分析の説明はたいてい「データの分散が最大になる方向を第 1 主成分とする」で始まり、次に「第 1 主成分に直交する方向のうち分散が最大のものを第 2 主成分とする」と続きます。この 2 文が、前節までの変分特徴づけと一字一句対応していることを確認しましょう。

射影の分散はレイリー商

$N$ 個の $d$ 次元データ $\bm{x}_1, \dots, \bm{x}_N$ が中心化されている（平均がゼロ）とします。単位ベクトル $\bm{u}$（$\|\bm{u}\| = 1$）方向へ射影した値は $z_i = \bm{u}^\top \bm{x}_i$ です。この射影値の標本分散は

$$ \mathrm{Var}[z] = \frac{1}{N-1}\sum_{i=1}^{N} (\bm{u}^\top\bm{x}_i)^2 = \frac{1}{N-1}\sum_{i=1}^{N} \bm{u}^\top \bm{x}_i \bm{x}_i^\top \bm{u} = \bm{u}^\top \left( \frac{1}{N-1}\sum_i \bm{x}_i\bm{x}_i^\top \right) \bm{u} $$

と変形できます。括弧の中は標本共分散行列 $\bm{S}$ そのものです。よって

$$ \begin{equation} \mathrm{Var}[z] = \bm{u}^\top \bm{S} \bm{u} = R_{\bm{S}}(\bm{u}) \quad (\|\bm{u}\|=1) \end{equation} $$

つまり 「方向 $\bm{u}$ に射影したときのデータの分散」は、共分散行列のレイリー商にほかなりません。$\bm{u}$ を単位ベクトルに制限しなければ、$\bm{u}^\top\bm{S}\bm{u}/\bm{u}^\top\bm{u}$ という商の形が「長さの効果を除いた、方向あたりの分散」を与えます。

PCA の各主張の翻訳

これで辞書ができました。

とくに最後の行が重要です。「$k$ 個の主成分で張る部分空間が、射影後の分散の総和を最大にする」という PCA の最適性は、レイリー商を部分空間に一般化したトレース最大化

$$ \max_{\bm{U}^\top \bm{U} = \bm{I}_k} \mathrm{tr}(\bm{U}^\top \bm{S} \bm{U}) = \lambda_1 + \cdots + \lambda_k $$

（Ky Fan の定理）として書けます。1 方向ずつ貪欲に選んでも、$k$ 次元をまとめて最適化しても、同じ答えになる — これが PCA が「貪欲なのに最適」である理由です。

数値実験での確認

$2$ 次元の相関のあるデータを 300 点生成し、標本共分散行列を計算すると

$$ \bm{S} = \begin{pmatrix} 2.5962 & 1.4161 \\ 1.4161 & 1.1304 \end{pmatrix} $$

となりました。固有値は $\lambda_1 = 3.4578$、$\lambda_2 = 0.2689$ で、第 1 主成分の方向は $x$ 軸から $31.32^\circ$ です。一方、単位円を 18 万点に区切って射影分散 $\bm{u}^\top\bm{S}\bm{u}$ を総当たりで最大化すると、最大値 $3.4578$、角度 $31.318^\circ$ が得られます。総当たり探索と固有値分解が同じ答えを返すことが確認できました。

さらに、この最適方向へ実際にデータを射影して標本分散を計算すると $3.45780$ となり、固有値 $\lambda_1 = 3.45780$ と 5 桁以上一致します。固有値は抽象的な代数量ではなく、「その方向に射影したときのデータのばらつきの大きさ」という具体的な統計量なのです。

左の散布図では、データ雲が最も長く伸びた向きに赤い第 1 主成分軸が沿い、オレンジの点はデータをその軸へ落とした射影です。射影点の広がりが最も大きくなる向きを探すのが PCA の作業です。右のグラフは射影方向を $0^\circ$ から $180^\circ$ まで回したときの分散で、$31.32^\circ$ で最大固有値 $3.4578$ の水平線に接し、そこから $90^\circ$ 離れた $121.32^\circ$ で最小固有値 $0.2688$ に接しています。この曲線が共分散行列のレイリー商そのものであり、山頂を探す作業と第 1 主成分を求める作業が同一だと分かります。

寄与率は $\lambda_1/(\lambda_1+\lambda_2) = 3.4578/3.7266 = 0.9279$、つまり第 1 主成分だけで全分散の $92.8\%$ を説明します。トレース $\mathrm{tr}(\bm{S}) = 2.5962 + 1.1304 = 3.7266$ が固有値の和 $3.4578 + 0.2689 = 3.7266$ と一致することも、この文脈では「軸を回転させても全分散は変わらない」という意味を持ちます。

デフレーション — 第 2 主成分を「引き算」で取り出す

第 2 主成分を求める素直な方法は、$\bm{v}_1$ 方向の成分を取り除いた行列

$$ \bm{S}’ = (\bm{I} – \bm{v}_1\bm{v}_1^\top) \, \bm{S} \, (\bm{I} – \bm{v}_1\bm{v}_1^\top) $$

を作って、また最大固有値を求めることです。この操作をデフレーションと呼びます。実際に計算すると $\bm{S}’$ の固有値は $\{2.66 \times 10^{-17}, 0.26883\}$ となり、$\lambda_1$ に対応する成分が数値誤差レベルまで消え、$\lambda_2 = 0.26883$ が最大固有値として残ります。その固有ベクトルの向きは $121.32^\circ$ で、第 1 主成分の $31.32^\circ$ からちょうど $90^\circ$ ずれています。

左の棒グラフは対数目盛で、デフレーション後に $\lambda_1 = 3.46$ に対応する固有値だけが $2.66 \times 10^{-17}$ という機械精度レベルまで落ち、$\lambda_2 = 0.269$ はそのまま残っていることを示します。ねらった 1 方向だけをピンポイントで消せているわけです。右のグラフでは、元の $\bm{S}$ の射影分散（青）にあった $31.32^\circ$ の大きな山が消え、紫の曲線では $121.32^\circ$ が新しい山頂になっています。谷だった場所が山になる、という入れ替わりが、「第 1 主成分に直交する中で分散最大」という逐次最大化の操作そのものです。

デフレーションと逐次最大化が同じ結果を返すのは偶然ではありません。射影 $\bm{I} – \bm{v}_1\bm{v}_1^\top$ をかけることは「$\bm{v}_1 \perp$ という制約を、制約なし問題に埋め込む」操作だからです。ただし数値計算の実務では、デフレーションを何度も繰り返すと丸め誤差が蓄積して直交性が崩れるという問題があり、大規模問題では Lanczos 法のような部分空間法が使われます。

理論と対応が付いたところで、次は「では実際にどうやって計算するのか」という数値アルゴリズムの話に移ります。レイリー商は評価するだけの道具ではなく、アルゴリズムを加速する部品でもあるのです。

レイリー商反復 — 良い近似値をシフトに使うと3次収束する

べき乗法の復習と限界

最大固有ベクトルを求める最も単純な方法はべき乗法です。適当な初期ベクトル $\bm{x}_0$ に $\bm{A}$ を繰り返し掛けて正規化するだけ。

$$ \bm{x}_{k+1} = \frac{\bm{A}\bm{x}_k}{\|\bm{A}\bm{x}_k\|} $$

固有基底で展開すれば、$k$ 回掛けたあとの係数は $c_i \lambda_i^k$ になるので、最大固有値の成分が相対的に一番速く育ちます。誤差（固有ベクトルとの角度の正弦）は $|\lambda_2/\lambda_1|^k$ のオーダーで減り、線形収束です。

問題は、この比が 1 に近いと絶望的に遅いことです。$\lambda_1 = 10$、$\lambda_2 = 8$ なら比は $0.8$。1 桁精度を上げるのに $\log 10 / \log(1/0.8) \approx 10$ 回かかります。

レイリー商は「二次の精度」を持つ

ここでレイリー商の隠れた性質が効いてきます。前に「固有ベクトルは $R$ の臨界点」と示しました。臨界点では 1 次の項が消えるので、近似固有ベクトルの誤差が $\varepsilon$ なら、レイリー商の誤差は $O(\varepsilon^2)$ に落ちるのです。テイラー展開で書けば、$\bm{x} = \bm{v}_1 + \varepsilon \bm{e}$（$\bm{e} \perp \bm{v}_1$、$\|\bm{e}\|=1$）として

$$ R(\bm{v}_1 + \varepsilon\bm{e}) = \frac{(\bm{v}_1 + \varepsilon\bm{e})^\top\bm{A}(\bm{v}_1+\varepsilon\bm{e})}{1 + \varepsilon^2} = \frac{\lambda_1 + 2\varepsilon\, \bm{e}^\top\bm{A}\bm{v}_1 + \varepsilon^2 \bm{e}^\top\bm{A}\bm{e}}{1+\varepsilon^2} $$

ここで $\bm{A}\bm{v}_1 = \lambda_1\bm{v}_1$ と $\bm{e} \perp \bm{v}_1$ から、1 次の項 $2\varepsilon\,\bm{e}^\top\bm{A}\bm{v}_1 = 2\varepsilon\lambda_1 \bm{e}^\top\bm{v}_1 = 0$ が消えます。分母を $1 – \varepsilon^2 + O(\varepsilon^4)$ と展開して掛け合わせると

$$ R(\bm{v}_1 + \varepsilon\bm{e}) = \lambda_1 + \varepsilon^2 \left( \bm{e}^\top\bm{A}\bm{e} – \lambda_1 \right) + O(\varepsilon^4) $$

となり、確かに誤差は $\varepsilon^2$ のオーダーです。ベクトルの精度が 1 桁なら、固有値の精度は 2 桁得られる — これがレイリー商を「エネルギー推定の道具」として使う理由でもあります（量子化学の変分法で、波動関数がそこそこでもエネルギーが精度良く出るのはこの性質のおかげです）。

逆反復とシフト

もうひとつの部品が逆反復（inverse iteration）です。シフト $\mu$ を使って $(\bm{A} – \mu\bm{I})^{-1}$ を掛ける反復を行うと、この行列の固有値は $1/(\lambda_i – \mu)$ になります。$\mu$ が $\lambda_j$ に近いほど $1/(\lambda_j – \mu)$ は巨大になるので、$\mu$ に最も近い固有値の成分が爆発的に増幅されます。収束率は

$$ \max_{i \neq j} \left| \frac{\lambda_j – \mu}{\lambda_i – \mu} \right| $$

なので、$\mu$ を $\lambda_j$ に近づけるほど速くなります。しかし当然、$\lambda_j$ を知らないから計算しているのであって、良いシフトを最初から用意することはできません。

2つを組み合わせる — レイリー商反復

そこで「毎回、現在のベクトルのレイリー商を次のシフトに使う」という発想が生まれます。これがレイリー商反復（Rayleigh Quotient Iteration, RQI）です。

1. 初期ベクトル $\bm{x}_0$（$\|\bm{x}_0\| = 1$）を選び、$\mu_0 = R(\bm{x}_0)$ とする
2. $(\bm{A} – \mu_k \bm{I})\bm{y} = \bm{x}_k$ を $\bm{y}$ について解く
3. $\bm{x}_{k+1} = \bm{y}/\|\bm{y}\|$ と正規化する
4. $\mu_{k+1} = R(\bm{x}_{k+1})$ に更新して 2 に戻る

なぜ速いのか、誤差の見積もりで説明できます。ベクトルの誤差（固有ベクトルとの角度の正弦）を $\varepsilon_k$ とすると、

• レイリー商の二次精度から、シフトの誤差は $|\mu_k – \lambda_j| = O(\varepsilon_k^2)$
• 逆反復 1 回の収束率は $O(|\lambda_j – \mu_k|)$ のオーダー（分母のギャップは定数程度）

なので、次のステップの誤差は

$$ \varepsilon_{k+1} = O(\varepsilon_k) \times O(\varepsilon_k^2) = O(\varepsilon_k^3) $$

3 次収束です。1 反復ごとに正しい桁数が 3 倍になる、という猛烈な速さです（対称行列の場合。非対称行列では 2 次収束）。

数値実験

固有値が $\{10, 8, 5, 3, 1.5, 0.5\}$ の $6 \times 6$ 対称行列で試します。

べき乗法は、レイリー商の誤差 $|R(\bm{x}_k) – 10|$ が次のように減ります。反復ごとの誤差比を見ると、$0.650, 0.696, 0.695, 0.684, \dots$ と推移し、$k=20$ 以降は $0.6400$ に張り付きます。これは理論値 $(\lambda_2/\lambda_1)^2 = (8/10)^2 = 0.64$ と完全に一致します（固有値の誤差はベクトルの誤差の 2 乗なので、比も 2 乗になります）。精度 $10^{-12}$ に到達するまでに 65 反復必要でした。

レイリー商反復は、べき乗法を 3 回だけ回して初期値としたあと、次のように進みます。

わずか 4 反復で機械精度に到達しました。べき乗法の 65 反復と比べて桁違いです。

左のグラフは縦軸が対数目盛です。べき乗法（青）はきれいな直線を描いており、これは「1 反復ごとに誤差が一定の割合で減る」線形収束の証拠です。一方レイリー商反復（赤）は直線ではなく、下向きに折れ曲がりながら急降下し、7 反復目（べき乗法によるウォームスタート 3 回を含む）で機械精度の底に達しています。右のグラフはべき乗法の誤差比 $e_{k+1}/e_k$ で、10 反復あたりから理論値 $(\lambda_2/\lambda_1)^2 = 0.64$ の赤い破線に張り付きます。固有値の誤差はベクトルの誤差の 2 乗なので、比も $(\lambda_2/\lambda_1)$ ではなくその 2 乗になる点に注意してください。

3 次収束であることは、角度誤差の列で確かめられます。$\varepsilon_{k+1}/\varepsilon_k^3$ を計算すると

$$ \frac{2.582\times10^{-1}}{(5.375\times10^{-1})^3} = 1.663, \quad \frac{1.913\times10^{-2}}{(2.582\times10^{-1})^3} = 1.111, \quad \frac{7.010\times10^{-6}}{(1.913\times10^{-2})^3} = 1.001 $$

と、ほぼ定数 $1.0$ に収束しています。もし 2 次収束なら $\varepsilon_{k+1}/\varepsilon_k^2$ が定数になるはずですが、そちらを計算すると $0.894 \to 0.287 \to 0.019$ とゼロに落ちていき、収束次数が 2 より高いことが分かります。3 で割ったときだけ定数になる — これが 3 次収束の証拠です。

左の棒グラフを見ると、灰色（2 乗で割った比）は $0.894 \to 0.287 \to 0.019$ と反復のたびに小さくなっていくのに対し、赤（3 乗で割った比）は $1.663 \to 1.111 \to 1.001$ と $1.0$ の点線へ収束しています。ある次数で割ったとき比が定数に落ち着けば、それが真の収束次数です。右のグラフは同じ行列・同じアルゴリズムで初期ベクトルだけを変えた比較で、ウォームスタートありは $\lambda_1 = 10$ に、なしは $\lambda_2 = 8$ に吸い込まれます。どちらも 2〜3 反復で収束しており「速さ」は同じですが、行き先が初期値だけで決まってしまうのがレイリー商反復の性格です。

落とし穴 — どの固有値に収束するかは選べない

RQI には重要な注意点があります。収束先を指定できないのです。同じ行列で、べき乗法による前処理をせず初期ベクトルをそのまま使うと、$\mu_0 = 6.766$ から出発して

$$ 6.766 \to 7.618 \to 7.994 \to 8.0000000033 \to 8.0 $$

と、最大固有値 $10$ ではなく $\lambda_2 = 8$ に収束しました。RQI は「初期のレイリー商に近い固有値」を掴みに行くので、最大固有値がほしいならべき乗法や部分空間反復で大まかな方向を作ってから RQI に渡す、という併用が実用的です。

もうひとつのコストの問題もあります。べき乗法は 1 反復が行列ベクトル積 $O(n^2)$ だけですが、RQI は毎回シフトが変わるので $(\bm{A}-\mu_k\bm{I})$ の分解をやり直す必要があり、密行列なら 1 反復 $O(n^3)$ です。反復回数が 10 分の 1 以下でも、$n$ が大きいと総コストで負けることがあります。実際の固有値ソルバ（LAPACK の対称固有値ルーチンなど）は、まず Householder 変換で三重対角化してから QR 法を回し、その内部でウィルキンソン・シフト（レイリー商シフトの改良版）を使う、という構成になっています。レイリー商は最終的なアルゴリズムの中でも、シフト戦略という形で生き続けているわけです。

なお、$\mu_k$ が固有値に近づくと $(\bm{A} – \mu_k\bm{I})$ は特異に近づき、線形方程式の解は数値的に不安定になります。直感的には破綻しそうですが、実は誤差が増幅される方向がまさに求めたい固有ベクトルの方向なので、正規化した後のベクトルは正しいままです。「ill-conditioned だが答えは合っている」という、数値線形代数でも屈指の美しい現象です。

ここまでの理論を、実際に手を動かして確認しましょう。

Pythonでの実装

以下では 3 つのことを確認します。(1) $2\times2$ 対称行列のレイリー商を単位円上で描き、極値が固有ベクトル方向で起きること。(2) べき乗法と RQI の収束回数の比較。(3) 共分散行列に適用して PCA の第 1 主成分と一致すること。

実験1: 単位円上のレイリー商

まず $\bm{A} = \begin{pmatrix} 5 & 2 \\ 2 & 2\end{pmatrix}$ について、単位円上の角度 $\theta$ を動かしながらレイリー商を評価します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

A = np.array([[5.0, 2.0], [2.0, 2.0]])
w, V = np.linalg.eigh(A)          # 昇順に固有値を返す
print("固有値:", w)               # [1. 6.]
print("固有ベクトル:\n", V)

# 単位円上のレイリー商（分母は1なので二次形式そのもの）
theta = np.linspace(0, 2 * np.pi, 200001)
X = np.vstack([np.cos(theta), np.sin(theta)])       # 2 x M
R = np.einsum("ij,jk,ki->i", X.T, A, X)             # 各列に対する x^T A x

print("最大値: %.10f  角度: %.4f 度" % (R.max(), np.degrees(theta[R.argmax()])))
print("最小値: %.10f  角度: %.4f 度" % (R.min(), np.degrees(theta[R.argmin()])))
# 解析解 R(θ) = 3.5 + 2.5 cos(2θ - φ), tanφ = 4/3 との比較
phi = np.arctan2(2.0, 1.5)
R_analytic = 3.5 + 2.5 * np.cos(2 * theta - phi)
print("解析解との最大差:", np.abs(R - R_analytic).max())

固有値は $\{1, 6\}$、単位円上の最大値は $5.9999999994$（角度 $26.5644^\circ$）、最小値は $1.0000000006$（角度 $116.5644^\circ$）となります。手計算した解析解 $R(\theta) = 3.5 + 2.5\cos(2\theta – \varphi)$ との差は機械精度レベル（$10^{-15}$ 程度）で、倍角公式による変形が正しかったことが確認できます。最大・最小の角度がちょうど $90^\circ$ 離れている点も、固有ベクトルの直交性の現れです。

次に、この関数を「角度に対するグラフ」と「単位円上の高さ」の 2 通りで描きます。

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# 左: 角度 vs レイリー商
ax = axes[0]
ax.plot(np.degrees(theta), R, lw=2, color="tab:blue", label="レイリー商 $R(\\theta)$")
ax.axhline(6, ls="--", color="tab:red", label="最大固有値 $\\lambda_1=6$")
ax.axhline(1, ls="--", color="tab:green", label="最小固有値 $\\lambda_2=1$")
ax.axvline(26.565, ls=":", color="tab:red")
ax.axvline(116.565, ls=":", color="tab:green")
ax.set_xlabel("方向の角度 [度]"); ax.set_ylabel("レイリー商の値")
ax.set_title("単位円を一周したときのレイリー商（周期180度）")
ax.legend(); ax.grid(alpha=0.3)

# 右: 単位円上に R の値を色で描く
ax = axes[1]
sc = ax.scatter(X[0, ::200], X[1, ::200], c=R[::200], cmap="viridis", s=18)
for i, (lab, col) in enumerate([("最小方向 $v_2$", "tab:green"), ("最大方向 $v_1$", "tab:red")]):
    v = V[:, i]
    v = v * np.sign(v[0] if v[0] != 0 else v[1])      # 符号を揃えて右向きにする
    ax.annotate("", xy=1.35 * v, xytext=(0, 0),
                arrowprops=dict(arrowstyle="->", lw=2.5, color=col))
    ax.text(1.45 * v[0], 1.45 * v[1], lab, color=col, ha="center")
plt.colorbar(sc, ax=ax, label="レイリー商の値")
ax.set_aspect("equal"); ax.set_xlim(-1.9, 1.9); ax.set_ylim(-1.9, 1.9)
ax.set_title("単位円上のレイリー商と固有ベクトル方向"); ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()

左のグラフは、レイリー商が角度に対して周期 $180^\circ$ のきれいな正弦波になることを示しています。中心線が固有値の平均 $3.5$、振幅が固有値の差の半分 $2.5$ です。右の図では、色が最も明るい（値が大きい）方向に赤い矢印（$\bm{v}_1$）が、最も暗い方向に緑の矢印（$\bm{v}_2$）が刺さっており、極値がぴったり固有ベクトル方向で起きることが視覚的に分かります。2 本の矢印が直交していることも一目瞭然です。

実験2: べき乗法とレイリー商反復の収束比較

次に $6\times6$ の対称行列で、2 つの反復法の収束速度を比べます。固有値を $\{10, 8, 5, 3, 1.5, 0.5\}$ に設定し、最大固有値の比が $0.8$ とやや厳しい条件にします。

import numpy as np

rng = np.random.default_rng(0)
n = 6
Q, _ = np.linalg.qr(rng.standard_normal((n, n)))      # ランダムな直交行列
lams = np.array([10.0, 8.0, 5.0, 3.0, 1.5, 0.5])
A = Q @ np.diag(lams) @ Q.T
A = (A + A.T) / 2                                     # 対称性を数値的に保証

w, V = np.linalg.eigh(A)
v_top = V[:, -1]                                      # 最大固有値の固有ベクトル

x0 = rng.standard_normal(n); x0 /= np.linalg.norm(x0)

def rayleigh(A, x):
    return x @ A @ x / (x @ x)

# --- べき乗法 ---
x = x0.copy(); err_pm = []
for k in range(80):
    x = A @ x
    x /= np.linalg.norm(x)
    err_pm.append(abs(rayleigh(A, x) - 10.0))

print("べき乗法: 誤差比 (k=20付近) =", err_pm[21] / err_pm[20])   # → 0.64
print("理論値 (lam2/lam1)^2 =", (8 / 10) ** 2)
for tol in [1e-4, 1e-8, 1e-12]:
    idx = [k for k, e in enumerate(err_pm) if e < tol]
    print("  精度 %.0e に必要な反復数: %d" % (tol, idx[0] + 1))

べき乗法の誤差比は $k=20$ 付近で $0.6400$ となり、理論値 $(\lambda_2/\lambda_1)^2 = 0.64$ と一致します。必要な反復数は精度 $10^{-4}$ で 24 回、$10^{-8}$ で 44 回、$10^{-12}$ で 65 回。桁を 4 つ上げるごとに約 20 反復という、線形収束そのものの挙動です。

続いてレイリー商反復です。最大固有値を狙うため、べき乗法を 3 回だけ回して初期ベクトルを作ります。

import numpy as np

def rqi(A, x, n_iter=5):
    """レイリー商反復。各反復のシフト値と固有ベクトル角度誤差を返す"""
    x = x / np.linalg.norm(x)
    mu = x @ A @ x
    hist = [(mu, x.copy())]
    I = np.eye(A.shape[0])
    for _ in range(n_iter):
        try:
            y = np.linalg.solve(A - mu * I, x)       # 逆反復（シフト付き）
        except np.linalg.LinAlgError:
            break                                     # 完全に特異 = 収束済み
        x = y / np.linalg.norm(y)
        mu = x @ A @ x                                # シフトをレイリー商に更新
        hist.append((mu, x.copy()))
    return hist

def sin_angle(x, v):
    c = min(abs(x @ v / (np.linalg.norm(x) * np.linalg.norm(v))), 1.0)
    return np.sqrt(max(0.0, 1 - c ** 2))

# べき乗法3回でウォームスタート
x = x0.copy()
for _ in range(3):
    x = A @ x; x /= np.linalg.norm(x)

hist = rqi(A, x, n_iter=5)
eps = []
for k, (mu, xk) in enumerate(hist):
    s = sin_angle(xk, v_top)
    eps.append(s)
    print("k=%d  mu=%.15f  |mu-10|=%.3e  sin=%.3e" % (k, mu, abs(mu - 10.0), s))

print("\n3次収束の確認 eps[k+1]/eps[k]^3:")
for i in range(3):
    print("  %.4f" % (eps[i + 1] / eps[i] ** 3))

出力は本文の表と同じで、$\mu$ は $9.4055 \to 9.8663 \to 9.99927 \to 9.9999999999 \to 10.0$ と進み、4 反復で機械精度に達します。角度誤差の比 $\varepsilon_{k+1}/\varepsilon_k^3$ は $1.663, 1.111, 1.001$ と定数に収束しており、3 次収束が数値的に裏付けられました。べき乗法が 65 反復かかった精度に、RQI は前処理込みでも合計 7 反復で到達しています。

このアルゴリズムの「弱点」も再現できます。ウォームスタートなしで生の初期ベクトルから始めると、$\mu_0 = 6.766$ から出発して $\lambda_2 = 8$ のほうに吸い込まれます。速いが行き先を選べない、という RQI の性格が数値で確認できます。

import numpy as np

hist_raw = rqi(A, x0, n_iter=5)
print("ウォームスタートなし（初期レイリー商 %.3f）:" % (x0 @ A @ x0))
for k, (mu, _) in enumerate(hist_raw):
    print("  k=%d  mu=%.12f" % (k, mu))
# → 6.765828 → 7.618400 → 7.994054 → 8.000000003 → 8.000000000

出力のとおり、収束先は $8$ です。初期ベクトルのレイリー商 $6.766$ は $10$ よりも $8$ に近く、RQI は「最も近い固有値」を掴みます。最大固有値がほしければ、べき乗法や部分空間反復で先に大まかな方向を作る必要がある、という設計指針がここから読み取れます。

実験3: 共分散行列に適用して PCA と一致することを確認

最後に、レイリー商最大化が PCA の第 1 主成分と同じ答えを返すことを確かめます。

import numpy as np

rng = np.random.default_rng(42)
N = 300
mean = np.array([2.0, 1.0])
C_true = np.array([[3.0, 1.6], [1.6, 1.2]])
L = np.linalg.cholesky(C_true)
X = (L @ rng.standard_normal((2, N))).T + mean        # 相関のある2次元データ

Xc = X - X.mean(axis=0)                               # 中心化
S = Xc.T @ Xc / (N - 1)                               # 標本共分散行列
print("標本共分散行列:\n", S)

w, V = np.linalg.eigh(S)
print("固有値:", w)                                    # [0.26883 3.45780]
# 固有ベクトルの符号は不定なので、角度を0〜180度に畳んで表示する
ang_pc1 = np.degrees(np.arctan2(V[1, 1], V[0, 1])) % 180
print("第1主成分の角度: %.3f 度" % ang_pc1)

# 単位円を総当たりして射影分散（=レイリー商）を最大化
th = np.linspace(0, np.pi, 180001)
U = np.vstack([np.cos(th), np.sin(th)])
R = np.einsum("ij,jk,ki->i", U.T, S, U)
u_best = U[:, R.argmax()]
print("総当たり最大値: %.6f  角度: %.3f 度" % (R.max(), np.degrees(th[R.argmax()])))
print("固有値による予測: %.6f" % w[-1])
print("最適方向へ射影した実データの分散: %.6f" % np.var(Xc @ u_best, ddof=1))
print("寄与率: %.4f" % (w[-1] / w.sum()))

総当たり探索の最大値 $3.457796$ は最大固有値 $3.457796$ と一致し、角度も $31.318^\circ$（固有ベクトルは $31.318^\circ$）で一致します。さらに、その方向へ実データを射影した標本分散も $3.457796$ となり、「固有値 = 射影分散」という対応が数値でも成立しています。寄与率は $0.9279$ なので、第 1 主成分だけで全分散の $92.8\%$ を説明できるデータでした。

散布図の上に主成分軸を描き、方向を変えたときの射影分散の変化も見てみましょう。

import numpy as np
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

ax = axes[0]
ax.scatter(Xc[:, 0], Xc[:, 1], s=14, alpha=0.5, label="中心化したデータ")
for i, (col, lab) in enumerate([("tab:green", "第2主成分"), ("tab:red", "第1主成分")]):
    v = V[:, i] * np.sqrt(w[i]) * 2.5
    ax.annotate("", xy=v, xytext=-v,
                arrowprops=dict(arrowstyle="<->", lw=2.5, color=col))
    ax.text(v[0] * 1.15, v[1] * 1.15, lab, color=col)
ax.set_aspect("equal"); ax.grid(alpha=0.3)
ax.set_xlabel("$x_1$"); ax.set_ylabel("$x_2$")
ax.set_title("データ散布図と主成分軸（軸長は標準偏差に比例）"); ax.legend()

ax = axes[1]
ax.plot(np.degrees(th), R, lw=2, label="射影分散 $u^T S u$")
ax.axhline(w[-1], ls="--", color="tab:red", label="最大固有値 $\\lambda_1$")
ax.axhline(w[0], ls="--", color="tab:green", label="最小固有値 $\\lambda_2$")
ax.axvline(np.degrees(th[R.argmax()]), ls=":", color="tab:red")
ax.set_xlabel("射影方向の角度 [度]"); ax.set_ylabel("射影後の分散")
ax.set_title("方向を変えたときの分散の変化"); ax.legend(); ax.grid(alpha=0.3)

plt.tight_layout(); plt.show()

左の散布図では、データ雲が最も長く伸びた向きに赤い矢印（第 1 主成分）が沿っています。右のグラフは、射影方向を $0^\circ$ から $180^\circ$ まで回したときの分散の変化で、$31.3^\circ$ で最大固有値の水平線に接し、そこから $90^\circ$ 離れた $121.3^\circ$ で最小固有値に接します。「分散最大の方向を探す」という PCA の言葉が、この曲線の頂上を探すことに等しいと一目で分かります。

発展 — 一般化レイリー商とその応用

レイリー商の考え方は、分母を単なるノルムから別の二次形式に取り替えることで大きく広がります。一般化レイリー商は、対称行列 $\bm{A}$ と対称正定値行列 $\bm{B}$ に対して

$$ \begin{equation} R_{\bm{A},\bm{B}}(\bm{x}) = \frac{\bm{x}^\top \bm{A} \bm{x}}{\bm{x}^\top \bm{B} \bm{x}} \end{equation} $$

と定義されます。同じようにラグランジュ未定乗数法で臨界条件を求めると、$2\bm{A}\bm{x} – 2\lambda\bm{B}\bm{x} = \bm{0}$ すなわち一般化固有値問題

$$ \bm{A}\bm{x} = \lambda \bm{B}\bm{x} $$

が現れます。$\bm{B} = \bm{I}$ が通常の固有値問題です。この一般化版は次のような場面で使われます。

構造振動: 剛性行列 $\bm{K}$ と質量行列 $\bm{M}$ について $\bm{K}\bm{\phi} = \omega^2 \bm{M}\bm{\phi}$ が振動モードを与え、レイリー商 $\bm{\phi}^\top\bm{K}\bm{\phi}/\bm{\phi}^\top\bm{M}\bm{\phi}$ は「ポテンシャルエネルギー / 運動エネルギー」の比、すなわち固有角振動数の 2 乗です。レイリー卿がこの量を導入した本来の文脈がこれで、モード形状をおおまかに推定するだけで固有振動数がかなり正確に出る（二次精度）という性質が、実務での近似計算に重宝されてきました。

線形判別分析（LDA）: クラス間分散 $\bm{S}_B$ とクラス内分散 $\bm{S}_W$ の比 $\bm{w}^\top\bm{S}_B\bm{w}/\bm{w}^\top\bm{S}_W\bm{w}$ を最大化すると、クラスが最もよく分かれる射影方向が得られます。PCA が「分散最大」だったのに対し、LDA は「分離度最大」を一般化レイリー商で表現しています。

グラフ分割: グラフラプラシアン $\bm{L}$ に対する $\bm{x}^\top\bm{L}\bm{x}/\bm{x}^\top\bm{x}$ は「隣接ノード間の値の差の二乗和」を測ります。定数ベクトル（自明な最小固有ベクトル）に直交する制約下でこれを最小化すると、第 2 固有ベクトル（フィードラーベクトル）が得られ、その符号でグラフを 2 分割するのがスペクトラルクラスタリングです。ここでも「$\bm{v}_1$ に直交する中で最適化」という逐次的変分特徴づけがそのまま使われています。

部分空間法（Rayleigh–Ritz）: 大規模行列では、$k$ 次元の部分空間 $\mathcal{K}$ に問題を射影し、小さな $k \times k$ 行列 $\bm{U}^\top\bm{A}\bm{U}$ の固有値で元の固有値を近似します。この近似値をリッツ値と呼び、レイリー商の部分空間版そのものです。クリロフ部分空間を使えば Lanczos 法になります。

これらはすべて「制約付きでレイリー商を最適化する」という同じ骨格を持っています。個別のアルゴリズムとして覚えるより、レイリー商という 1 つの型として理解しておくほうが、応用も転用も効きます。

まとめ

本記事では、レイリー商 $R(\bm{x}) = \bm{x}^\top\bm{A}\bm{x}/\bm{x}^\top\bm{x}$ を通じて、固有値を最大化問題として捉える視点を解説しました。

• レイリー商はスケール不変 — 分子・分母がともに 2 次なので方向だけの関数になり、最適化を単位球面上の問題に落とせる
• 勾配が固有方程式の残差 — $\nabla R = \frac{2}{\bm{x}^\top\bm{x}}(\bm{A}\bm{x} – R(\bm{x})\bm{x})$ なので、臨界点であることと固有ベクトルであることが完全に同値
• 固有基底で展開すると凸結合 — $R(\bm{x}) = \sum_i w_i \lambda_i$（$w_i = c_i^2/\sum c_k^2$）なので $\lambda_n \leq R \leq \lambda_1$ が即座に従い、等号は固有ベクトルで達成される
• クーラン・フィッシャー — 第 $k$ 固有値は「$k$ 次元部分空間上の最小レイリー商」の最大値。固有ベクトルを知らなくても書けるので、行列同士の固有値比較や摂動評価の基礎になる
• PCA は共分散行列のレイリー商最大化 — 射影分散 $\bm{u}^\top\bm{S}\bm{u}$ がレイリー商そのもの。「既存成分に直交する方向で分散最大」は逐次的変分特徴づけと一字一句対応する
• レイリー商反復は 3 次収束 — レイリー商の二次精度と逆反復の組み合わせで、数値実験でも角度誤差の比 $\varepsilon_{k+1}/\varepsilon_k^3$ が $1.0$ に収束。べき乗法の 65 反復に対して 4 反復で機械精度に到達する

「固有値は方程式の根である」という代数的な見方と、「固有値は最適化問題の最適値である」という変分的な見方 — 後者を身につけると、PCA・スペクトラルクラスタリング・構造振動・量子化学の変分法・行列ノルムといった、一見バラバラな話題が一本の線でつながって見えてきます。次に固有値が出てくる場面に遭遇したら、「これは何を最大化しているのか」と問い直してみてください。

次のステップとして、以下の記事も参考にしてください。

• 固有値分解の理論と応用を完全解説する
• 主成分分析（PCA）の理論と導出
• 【Python】固有値と固有ベクトルの幾何学的特徴を可視化して理解する
• 【線形代数】行列とベクトルの2次形式や微分の公式を解説