Arnoldi法(Lanczos法)についてのメモ

動機

久しぶりのブログエントリになります。
休んでる間何もしてなかったかというと、勿論そうではなくて、比重がインプットに偏っていたのが原因でした。
偏りなくできるのが理想ですが、あまり時間がない時はひとつのことに集中するのが、経験上効率的だと思っているので、今度書く時もまた半年ぶりとかになる可能性もありますが、まぁ自分用のメモみたいなものなので、気楽にやっていきたいと思います。

前置きはこの辺にして、本題に入ります。
最近、個人的な興味があって固有値の求め方を調べていました。
固有値は、特異値分解において重要な役割を持ち、そして特異値分解は主成分分析、潜在意味インデキシングといった応用の基本となるアルゴリズムです。特異値分解自体は、それ以外にも非常に多数の応用が考えられますが、詳細については省略します。

実際に固有値を求める方法は、学部時代に固有方程式を手計算で解き、また数値解析の授業でべき乗法とかヤコビ法を授業で聞いたなぁという程度でしたが、自然言語処理等でよく使われる大規模疎行列でも高速に計算可能と言われる方法について、一度はちゃんと勉強しておかなければいけないと思ったのが、今回調査したきっかけとなります。

なお、実用上は大規模疎行列の特異値分解や固有値を求めることが必要な状況であれば、ARPACKやredsvdを利用すれば良いと思いますので、別に興味がなければ知らないままでいいと思います。

事前準備

そもそも固有値は行列 ${ \mathbf{A} }$ に対して、下記の式で成り立つような ${ \lambda }$ 及びベクトル ${ x\neq0 }$ が存在する時に ${ \lambda }$ を固有値と呼びます。

${ \displaystyle \mathbf{A}\mathbf{x}=\lambda\mathbf{x} }$

固有値を考える上では複素ベクトル空間で考えなければならないので、 ${ \mathbf{A} }$ の随伴行列(転置行列の複素共役行列)を ${ \mathbf{A}^{*} }$ と表すとします。
ここで、 ${ \mathbf{A}^{*}\mathbf{A} = \mathbf{I} }$ ( ${ \mathbf{I} }$ は単位ベクトル)が成り立つ行列をユニタリ行列と言います。
ユニタリ行列は、各列ベクトル(または各行ベクトル)を ${v_1..v_m}$ とすると、内積が以下のように記述できるので、各列ベクトルは正規直交基底であると言えます。

正規直交基底の定義
(内積に下記の関係が成り立つベクトル列であること)

${ (v_i, v_j)= \begin{cases} 0 & (i \neq j) \\ 1 & (i = j) \end{cases} }$

こういった正規直交基底は、グラムシュミットの正規直交化法というアルゴリズムで任意の独立なベクトル列から生成することが可能です。

グラムシュミットの正規直交化法



1. Start: Compute  := ||. If  = 0 stop, else  := /

2. Loop: For j = 2, . . . , p do:

    Compute  := ( , ) for i = 1, 2, . . . , j − 1

     :=  - 

     := ||

    If  = 0 then stop, else  := /

固有値を求めるフレームワーク

大規模行列 ${ \mathbf{A} }$ の固有値を効率的に求める際に一般的に、下記のようなフレームワークがあります。

${ \mathbf{A} }$ を部分空間に射影し、新しい行列 ${ \mathbf{H} }$ を得る
${ \mathbf{H} }$ の固有値を求める

ここで、 ${ \mathbf{H}}$ の固有値が ${ \mathbf{A} }$ の固有値の良い近似であることを証明すれば、このフレームワークは十分に機能することとなります。
一般に1の過程はReductionと呼ばれ、よく名前の上がるArnoldi法やLanczos法は、1の過程でKrylov部分空間に射影する手法のことであり、 ${ \mathbf{H} }$ の固有値を求める際はArnoldi法であればQR分解などまた別の手法を用います。

Krylov部分空間

Krylov部分空間は、べき乗法のアイディアを自然に発展させたものと考えることができます。
べき乗法は、初期値を適当なベクトル ${ \mathbf{x} }$ を取ると、
${ \displaystyle \mathbf{x_{n+1}} = \mathbf{A}\mathbf{x_{n}} }$
という計算を繰り返すことにより ${ \mathbf{x_{n}} }$ が絶対値最大の固有値の固有ベクトルに収束する性質を利用した固有値の求め方です。(証明はこちら)

べき乗法の場合、途中の計算結果は捨ててしまいますが、Krylov部分空間を用いた手法は途中の計算結果を利用して下記のようなベクトル列を構成します。

${ \displaystyle \{ \mathbf{x},\mathbf{Ax},\mathbf{A}^{2}\mathbf{x}, \mathbf{A}^{3}\mathbf{x} \dots \} }$

このベクトル列を基底とする空間をKrylov部分空間と言います。
Krylov部分空間は、べき乗法と同様、絶対値最大の固有値をよく近似することが可能であり、またそのままべき乗法を利用する場合は最大固有値しか得られないという問題があるのですが、一度Krylov部分空間に射影してから固有値を求めるという方法を取ることにより、より高速に計算できるようになります。

Arnoldi法

それでは実際に非対称行列 ${ \mathbf{A} }$ をArnoldi法を用いてKrylov部分空間に射影する方法について考えます。
対称行列の場合は、より計算が簡便なLanczos法を用いますが、これはArnoldi法の特殊ケースであると考えることができます。

${ \mathbf{A}}$ をKrylov部分空間に射影するアルゴリズムは下記のように計算ができます。



1. Start: Choose  where 

2. Loop: For j = 1, . . .  do:

    Compute  := ( , ) for i = 1, 2, . . . , j

     :=  - 

     := ||

    If  = 0 then stop, else  :=  /

このアルゴリズムが具体的にどういうことを行っているのかということを説明すると、まずn回目のイテレーションで、Krylov部分空間を構成するベクトル列からなる行列を ${ \mathbf{Q_{n}} = q_1,q_2,q_3,\dots,q_{n} }$ と置き、AをKrylov部分空間に射影した行列を ${ \mathbf{H_n} }$ と置くと、 ${ \mathbf{Q_{n}} }$ は各ベクトルをグラムシュミットの正規直交化法により正規化し、下記のように ${ \mathbf{H_n} }$ を求めていることとなります。

${ \displaystyle \mathbf{H_n} = \mathbf{Q_{n}^{*}}\mathbf{A}\mathbf{Q_n} }$

求められた ${ \mathbf{H_n} }$ は、ヘッセンベルグ標準形と呼ばれる形式となります。

ヘッセンベルグ標準形

${ \displaystyle \mathbf{H_n} = \begin{bmatrix} h_{1,1} & h_{1,2} & h_{1,3} & \cdots & h_{1,n} \\ h_{2,1} & h_{2,2} & h_{2,3} & \cdots & h_{2,n} \\ 0 & h_{3,2} & h_{3,3} & \cdots & h_{3,n} \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & h_{n,n-1} & h_{n,n} \end{bmatrix} }$

ヘッセンベルグ標準形はギブンス回転を用いたQR分解を実施することで固有値を高速に計算可能です。

${ \mathbf{H_n} }$ の固有値 ${ \simeq }$ ${ \mathbf{A} }$ の固有値の証明

では、次に ${ \mathbf{H_n} }$ の固有値が ${ \mathbf{A} }$ の固有値の近似となるかを考えます。
n回目の ${ mathbf{H_n} }$ を考えると、誤差項を ${ \mathbf{E_n} }$ と置くと

${ \displaystyle \mathbf{A}\mathbf{Q_n} = \mathbf{Q_n}\mathbf{H_n} + \mathbf{E_n} }$

が成り立ちます。
ここで ${ H_{n} }$ の固有ベクトルを ${ \mathbf{y} }$ 、固有値を ${ \lambda }$ と置いて、上式の両辺に ${ \mathbf{y} }$ をかけます。

${ \displaystyle \mathbf{A}\mathbf{Q_n}\mathbf{y} = \mathbf{Q_n}\mathbf{H_n}\mathbf{y} + \mathbf{E_n}\mathbf{y} \\ \mathbf{A}\mathbf{Q_n}\mathbf{y} = \lambda\mathbf{Q_n}\mathbf{y} + \mathbf{E_n}\mathbf{y} \\ \mathbf{A}\mathbf{Q_n}\mathbf{y} - \lambda\mathbf{Q_n}\mathbf{y} = \mathbf{E_n}\mathbf{y} \\ \left(\mathbf{A} - \lambda\mathbf{I}\right)\mathbf{u} = \mathbf{E_n}\mathbf{y} }$

ただし、 ${ \mathbf{u} = \mathbf{Q_n}\mathbf{y} }$

Arnoldi法の収束条件はこの誤差項が十分に0に近づく時ですので、

${ \displaystyle (\mathbf{A} - \lambda\mathbf{I})\mathbf{u} \simeq 0 }$

が成り立ち、これは ${ \mathbf{A} }$ の固有値を求める式と一致します。

まとめ

Lanczos法やArnoldi法は、固有値を直接求めるアルゴリズムではなく、Krylov部分空間に射影するアルゴリズム
Krylov部分空間はべき乗法のアイディアを発展させたもの
実用上は既存パッケージを使えばいい

Arnoldi法、Lanczos法はいくつも派生があり、非対称行列にも適用可能なLanczos法なども存在します。
部分空間に射影する方法は、redsvdが実装している方法も同様で、ランダム行列を使って射影を行い、特異値分解を行っているだけですのでアルゴリズム自体は実装も簡単です。ただ、理論的に近似誤差を求めているのですが、丁寧に説明が書かれているにも関わらず、まだよく理解できていません。特異値分解や固有値の求め方は、どちらかというと工学よりも理学(数学)よりですので、理論が難しくなりますね・・・。

また説明の際は簡単のため、各行列, ベクトルの次元を省略しています。
厳密な説明は参考文献を参照してください。

参考文献

数値解析(森正武)
数値計算の基礎
Numerical Methods for Large Eigenvalue Problems
英語版wikipedia(Arnoldi Iteration)
固有値問題ノートの補足
固有値計算アルゴリズムの最近の進展

動機