2004/10/12 解析学とコンピュータ 前回、前々回に、固有値問題について学んだ。 行列 A に対して A x = \lambda x を満たす \lambda のことを A の固有値と呼ぶ 対称な行列に関しては、ヤコビ法で Okey 対角行列の固有値は、自明 # 対角要素の値が固有値そのもの 固有値問題 ( 与えられた行列 A から、その固有値 \lambda を求めること ) を解くためには.. A から、固有値をかえない変換を行い、対角行列に変形すること => 対角化 # 固有値をかえない変換 => 相似変換 A が対称ならば 三重対角要素の下の要素を 0 にできる 対称行列ならば自動的に、上の要素も 0 になる 非対称行列の場合は、上の要素も残ってしまう # 保存則が成立する物理現象ならば、自動的に対称になる # 対称行列が扱えれば、ほぼ Okey 更に、対角化しようと思うと、fill in が生じてしまう。 => ヤコビ法ならば、全ての固有値が得られる 地震の解析 ( 振動問題 ) => 低い周波数 ( <=> 大きい固有値 ) が重要 実用上では、大きい固有値が ( 幾つか.. ) 判れば十分 ヤコビ法は、全ての固有値を求めるので、( 大きい固有値が何個か必要な場合は、 計算の.. ) 無駄がある => 羃乗法 ( 本日の main ) は、大きい固有値を数個求める方法 # 羃乗法、非対称行列を対象にするが、もちろん、対称行列でも使える # ただし、対称性を利用しないことになるので、無駄が生じる == 羃乗法 [仮定] A が、実数の固有値を持つとする # 対称行列の場合は、実数の固有値を持つことが証明できるので don,t care # 非対称行列の場合に必要な仮定 [方法] 適当なベクトル v_0 を初期値とし、 v_i = A v_{i-1} という計算をするだけ.. # これをすると何が起るのだろうか ? 今、(実数の) 固有値を絶対値の大きい方から並べる (重複しないとする) |\lambda_1| > |\lambda_2| > .. > |\lambda_n| 実は A は次の形になる A = \Sum_{k=1}^{n} \lambda_k x_k ( where x_k は \lambda_k の固有ベクトル ) # x_k は、互いに直交するように取れる !! 今、 v_0 = \Sum_{i=1}^n \alpha_i x_i とすれば、 v_k = A^k v_0 = A^k \Sum_{i=1}^n \alpha_i x_i = \Sum_{i=1}^n \alpha_i \lambda_i^k x_i = \lambda_i^k ( \Sum_{i=1}^n \alpha_i \right(\frac{\lambda_i}{\lambda_1}\left)^k x_i ) となる。ここで、仮定より、 |\lambda_1| > |\lambda_i| なので \right|\frac{\lambda_i}{\lambda_1}\left < 1 これより、 k \arrow \infty ならば、 right(\frac{\lambda_i}{\lambda_1}\left)^k \arrow 0 なので、 となる。 つまり、 v_k \arrow \lambda_1^k x_1 なので、(大きさは、ともかく..) 固有値ベクトル x_1 が得られる。 更に、 \frac{v_k|}{v_{k-1}} \arrow \lambda_1 なので、これで、\lambda_1 が得られる。 # 最大固有値 \lambda_1 がえられる。 最大固有値が得られたらどうするか ? A’= A - \lambda_1 x_1 x_1^t # この操作をデフレーションと呼ぶ とすると、A' は、\lambda_1 以外の固有値が A と同じだが \lambda_1 を 固有値を持たない ( 固有値が 0 ? ) ことになる A' に羃乗法を適用すれば、\lambda_2 が得られる。 # もし、これを繰り返せば、全ての固有値が求められるように思われるが.. # 無限桁の計算をしていれば、Yes 。実際には、有限桁の計算なので、後半の # \lambda を計算するときには、誤差が混入して、値が正しくなってしまう.. ## 誤差の伝播を考えると、せいぜい、大きい方から、数個だけ.. 絶対値が小さい方がほしければどうするか ? A の固有値が \lambda_1 < > \lambda_2 < .. ならば A^{-1} の固有値が \frac{1}{\lambda_1} < \frac{1}{\lambda_2} < .. なので、A^{-1} に羃乗法を使えば、A の最小の場合が得られる。 # 実際は、逆行列を求めずに、連立方程式を解く形にする ## v_{i+1} = A^{-1} v_i となるので、 ## A v_{i+1} = v_i という連立方程式を解く どうしても全部の \lambda が欲しいのであれば ? 全ての計算を倍精度で上手く行えば、誤差の拡散を防ぐことができる 羃乗法は、0 ベクトル以外の v_0 を与える # 乱数にするとか all-1 にするとか .. がそれを何にするかによって、収束速度が異る。 \lambda_1 の計算 ( レーリー商 [ Rayleigh Method ] ) # ベクトルの大きさの比を計算するのは大変なので、内積で計算する \lambda_1 = \frac{(v_k v_{k+1})}{(v_k v_k)} # \lambda_1^n の要素が含まれるので、オーバーフローがおきることがある == # 時間に余裕があるので、前回の補足 ヤコビ法を直接適用すると、結構、時間が掛る # 相似変換は、2 x 2 の操作を行うため.. fill-in がおきる => 三重対角までなら、fill-in がおきない => 変換の回数が固定 ( 予測できる ) # 三重対角から対角までが問題 # 非対称の場合は、三重対角の下の部分だけ 0 にできる # => 上の部分は 0 にならない # ヘッセンベルグ行列 ギブンス法 与えられた行列 A をヘッセンベルグ行列に変換すること # 三重対角要素を 0 にしないことを除けば、ヤコビ法と同じ作業 (相似変換) # 村田健郎先生の方法 ( 以下、Text の説明 ) ハウスホルダ法 n x n の行列 A が与えられた時に、 T = P^t A P が、ヘッセンベルグ行列になるような P を求める方法 ここで、v が、A の固有ベクトルならば、 P = I_n - 2\frac{v v^t}{v^t v} が成立する。 # ハウスホルダ法はギブンス法の、計算を効率良く行っているだけで同じ手法 # # 計算順序が異る # ギブンス法では、要素を一つずつ消す # 計算が独立なので、並列化向き # ハウスホルダ法、まとめて計算する # 計算が逐次的なので、ベクトル化向き [まとめ] 対称行列 ヤコビ法 回数が不定 ギブソン / ハウスホルダー 三重対角 バイセクション/QR 法 非対称行列 ベッセンベルグ行列 QR 方法 密行列 / バンド行列 / ランダムースパース行列 実数要素 / 複素数要素 ライブラリ EISPACK 研究者 村田 健郎 先生 ( 既に、リタイア ? ) 別府 先生 電通大学の 片桐 先生 == 連立方程式を解いていて、解けなかったことがある ピボットが 0 になると困る # 正則行列ではない場合など.. 固有値問題でも、同様な問題が生じることがある 普通 ( 物理現象をモデル化したら.. ) は、ヘッセンベルグ形までは Okey 上手く行かないような場合の條件は ? # 戸川先生の分類 \lambda_1 = \lambda_2 の場合 ( A 型 ) x_1 = x_2 の場合 ( B 型 ) 両方が生じる場合 ( AB 型 ) # 物理的には普通、おきない ( 数学的にはあり ) # 非常に稀な状況 このような特別な性質を持つ要素が一つあると 5 重対角になる # cf 伊理先生の本 ( 教育出版 ? ) # アルゴリズム的には、zero div がおきる == 今週の演習 羃乗法 p.68 # 演習的には、山を越した == 応用 普通は、線型方程式を解く # 偏微分方程式の差分化 非線形方程式を解くことは稀 光ファイバー ファイバー内部の光の反射を解くのは 4 次方程式 問題をおきたときにテキストを読み直すためのポイントの紹介を行っている