2005/10/11 数値解析学と演習 福井 先生 # 前回 ( 先々週 ) に引続き、固有値問題のヤコビ法 固有値問題 # 弦でいうと、両端を固定して振動させた時に出てくる周波数 A x = \lambda x \lambda_i ( i = 1 - n ) 固有値 x_i 固有ベクトル 固有値の性質 ( 実対称行列の場合は、全ての根が異るので.. ) (x_i x_j) = 0 ( i \ne j ) 二次元の場合 ( 楕円の軸 ) 固有ベクトルは、長軸、短軸の方向を表す。 固有値は、長軸、短軸の長さを表す。 「固有値問題が解けた」ということは、与えられた行列 A に対して、次の変形ができたということ a_1' A => ( a_2' ) 対角化する ( 但し、固有値を変えない状況で..) ... a_ n' 固有値を変えない変換 => 相似変換 A => PAP^t A の要素の一つを 0 にできる 二つの次元からなる平面上の回転変換になる 回転には向きがあるので、そのような変換は二つある。 => 数値計算的には、回転量が少い方を選ぶ。 # 固有値問題は、代数方程式と関係がある # 二次元の場合は、二次方程式と同じ # 四次元までは、巧く行く ( 方程式に公式がある ) # 五次元以上は、fill-in ( 一旦 0 になったものがまた、0 以外になる.. ) がおきて、巧く行かない(対角要素以外を 0 に出来無い)可能性がある.. 相似変換 二つの行,列に関する回転 i j 1 1 i cos x sin x 1 1 j -sin x cos x 1 1 # cf. text p.169 二つの行、列を同時に変更するので、一つの要素を 0 にした時に、他の要素が影響を受けて fill-in がおきる。 # 三重対角まではできる !! ## このままだと解けない ( 対角以外を 0 にできない ) のだが.. 「数値的」には.. 相似変換を行うと対角要素の絶対値が大きなる ( 証明可能 !! ) 対角要素以外は 0 に近付く ( 数値的には、0 に近いければ 0 とみなしてよい ) # どの位の速度で 0 に近付くかは、行列による もし、絶対値最大の非対角要素を選んで、0 にすれば、非対角要素は、速く 0 に収束する ( 反復回数が少くなる ) のだが... => 最大値の要素を探すのが大変 ( O(\frac{n^2}{2}) ) => 探すのは程々に.. 閾値 ( 10^{-2} 程度.. ) を決めてそれより大きい物を選択する ( 0 にする ) => 閾値ヤコビ方法 # 閾値をいくつにするかは、問題に依存する # 計算機では、演算 ( 四則 ) の他に、比較も計算量に関係する ヤコビ法では、( 数学的には、回転角度が必要なのだが.. ) 計算には、直接に、角度が現れない ( 間接的に、sin/cos などの三角比の形で現れる ) 角度 \thita = \frac{1}{2}\arctan{\frac{2a_{pq}}{a_{pp}-a_{qq}}} の計算には、\arctan があるので計算するのは (数値計算的には.. ) 嫌 # 実際に必要なのは角度ではなく三角比なので、.. # => 角度を求めずに三角比を求める wa = \frac{a_{pp}+a_{qq}}{2} sa = \frac{a_{pp}-a_{qq}}{2} r = \sqrt{sa^2+a_{pq}^2} \cos{\thita} = \frac{1}{\sqrt{2}}\sqrt{1\pm \frac{sa}{r}} \pm の所で、角度が二つ出ることになるが、絶対値の小さい方を選ぶ \sin{\thita} = \frac{a_{pq}}{1r \cos{\thita}} # ヤコビ法では、角度を経由せずに、三角比を求める点がポイント # ウィルキンソンの手法 ( 誤差を小くするための工夫がされている ) ## ウィルキンソン ## ぶ厚い本をかいている ( 数値計算の大家 ) ## => 理解しにくい.. ( まともにやったら、もっと厚くなる.. ? ) ## 途中に、突然式が現れて、良く解らない ### 文献検索で「ウィルキンソン」とすれば必ず出てくる # 数学で「適当な値」といったら、「その値を定める良い方法がない ( ので、その場、その場で决める.. )」ということを意味する。 # => これを决めるのが一番面倒 !! # ヤコビ法で、どの要素を 0 にするかは、難しい問題 # => どの閾値 ? どこから探す ? # => プログラムを作る段階で困る # => 問題領域が狭ければ、ある程度、决められるのだが.. # => 汎用の場合は色々と考える必要がある # => EISPACK ( インターネットから得られる ) # 本もある # これを手で入力してみた ..? # a/p/q という表現が.. ? # a/(p*q) だとなんと OverFlow になってしまう !!! ヤコビ法で、座標変換を繰り返すと、対角要素以外は、ほとんど 0 になる。 => その場合は、対角要素 ( n 個 ) が固有値になる 固有ベクトルが欲しいなら.. => 単位行列に相似変換で使った行列を同時にかければ良い。 固有ベクトルにおける注意 固有ベクトルは「向き」だけを表現している 絶対値 ( 長さと向き ) は決らない # 長さに相当するのは、固有値 x が固有ベクトルなら cx も固有ベクトル 固有値をライブラリで求めると.. 固有ベクトルの長さがどうなるかわからない ( 長さが 1 とは限らない ) => ノーマライズした方がよいかも # 複数のライブラリを使うと、固有値は同じだが、固有ベクトルの結果が (見かけ上) 異ることがある。 相似変換 \thita = 0 何もしない変換 \thita = \frac{\pi}{2} 交換になる # 計算機が出たころが、この方法が、良く使われた # その後、他の方法が出て、一時期廃れたが、今は、見直されている # => 意外に素直な方法で、「嫌な問題」も素直に解ける == cf. 戸川先生 「大型行列解法」 ヤコビ法 a_11 a_12 .. a_1n a_11' a_21 a_22 .. a_2n => a_22' .. a_33' a_n1 a_n2 .. a_nn a_44' # いきなり、対角要素意外を全て 0 にしようとする 取りあえず、三重対角行列までの形にする # 一旦 0 にしたものを 0 のままで変換可能 ## fill-in をさける 三重対角行列 ( 中継点数 ) を経由 要素数が n^2 から 3n まで減る ( 対称なら 2n 程度 ) デターミナント ( 行列式 ) を 0 にする値が、固有値 # 踈行列なので、データーミナントの計算が楽 => 二分方法で解くことが可能 バイセクション法 # 普通は「二分法」だが.. QR 法 ( QR 分解 ) # 密行列でも、同じことができるが、それだと演算量が多いので大変 # (三重対角行列にして..) 踈行列にしたので、これらの手法に意味がでてきた。 三重対角化する方法 ギブンス方法 ( ヤコビ法を 0 にする順を固定したもの ) ハウスフォルダー法 ( text p.280 ) # 良くみると似たような方法なのだが.. 三重対角化した行列から、固有値を求める バイセクション QR 法 # 組み合わせて使える # その他にも「ムラタ法」などもあるが.. # これらの方法の問題点 # => 固有値だけしか出ない # => 固有ベクトルが欲しけでば、別に計算する必要がある # 両方ほしければ、ヤコビ法も Okey == [演習] [課題番号] 10/11 [締切] 10/17 [プログラム] 問題 p.275 [入力データ] (戸川先生の本から拝借) <<問題>> A =[ 5 4 1 1 ] 4 5 1 1 1 1 4 2 1 1 2 4 <<答>> \lambda_1 = 10 [ 2 ] 2 1 1 \lambda_2 = 5 [ -1 ] -1 2 2 \lambda_3 = 2 [ 0 ] 0 -1 1 \lambda_4 = 1 [ -1 ] 1 0 0 # 注意、この答は、固有ベクトルの要素の符号のパターンが異る # => 二つの固有ベクトルの内積が 0 になるので、どこかで要素が正負が異る必要がある。