先週は、「ガウスの消去法」を学んだ。 (線型)連立一次方程式 Ax = b を解く。 行列の要素を、(方程式の意味を変えないように..) 変形し、 解きやすい形に変形する方法 ( = ガウスの消去法 ) [今週] 同じ係数行列で、複数の方程式を解く場合がある A x1 = b1 A x2 = b2 .. A xn = bn この場合はどうか ? 数学的な方法 A^{-1} を計算し xi = A^{-1}bi とする ... 利用しない => 計算量も多い ( 2/3 n^3 ) し、安定しない # cf 「ガウスの消去方」は 1/3 n^3 + 1/2 n^2 # もちろん、A^{-1} の係数が欲しいときは別だが.. # => 連立方程式を解くためには利用しないのが常識 もちろん、個々に、毎回ガウスの消去法を適用するのは、無駄 => A の係数の変形は毎回同じなので、そこを省きたい もし、bi が互いに無関係で、予め判っていれば.. A | b1 .. bn の形で、拡大行列を作り、一個の場合と同様に、(拡張された..) ガウスの消去法を利用すればよい。 # しかし... b_{i+1} = A^{-1}b_i のような形で、b_{i+1} が予めはわかってないと.. => 上記の方法は利用できない ( 困った.. ) # (復習)なぜ、A^{-1} を計算しないか ? # 係数を全部 0 にすると、無駄 # => 対角化だけなら、それほど無駄じゃない [LU 分解 方法] A = LU L は、(左)下三角行列 [ ただし、対角要素は 1 ] U は、(右)上三角行列 L = \Ary{ 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0 \\ l_{41} & l_{42} & l_{43} & 1 } U = \Ary{ u_{11} & u_{12} & u_{13} & u_{14} \\ 0 & u_{22} & u_{23} & u_{24} \\ 0 & 0 & u_{33} & u_{34} \\ 0 & 0 & 0 & u_{44} } # 要するに、A を(右)上三角行列化する部分を抜き出したい # なぜ、このような形に分解できるとうれしいか ? 形式的に Ly = b を考える。L は、下三角行列なので、簡単に解ける 1/2 n^2 一方、 A x = (LU) x = L(U x) = L y = b つまり、 U x = y なので、これも U が、上三角行列なので、簡単に解ける 1/2 n^2 ( 後退代入と同じ ) つまり、LU 分解ができていれば、後は、方程式毎に n^2 ( = 1/2 n^2 + 1/2 n^2 ) だけで、方程式は解けてしまう。 # LU 分解するのには 1/3 n^3 かかるので.. 合計で 1/3 n^3 + m ( n^2 ) ( m は方程式の set 数 ) # では、どうやって、A を L と U に分解するか ? # => 基本は、「ガウスの消去法」と同じ # 「ガウスの消去法」で、A を変形すると U と同じ # => L を計算する方法だけが問題。 L : A を U にする為の手順が記載されている行列 => A の係数を消すために行に掛けた数が L の要素 # A | E を拡張「ガウスの消去法」と同じ方法をとれば、 # U | L^{-1} が計算される # Ly = b の計算も y = L^{-1} b ですぐに求めることができる # => L の計算にも 1/3 n^3 かかる # bi の i の数が n 以下なら、LU, そうでなければ L^{-1} の方が # よいかも.. ? 配列の要素領域の再利用 U の計算の結果、どうせ、左下の要素は 0 が入っている。 # ちなみに L の対角要素も 1 => 覚えておく価値はない => そこに L の要素 ( 消すための係数 ) を保存 #! メモリの使用量を減らすための工夫だが、今の計算機で必要か ? # 解りやすさからいえば、別の行列にした方よい # しかし、この関係は、要するに A の次元と、L + U の次元 # が共に n^2 であることを表している。 # 線型空間は次元が同じならば、同形なので、同形対応が、 # 存在し、 A と LU (直交空間の積)が対応することわかる # ので、数学的な背景にもなっている。 # 本来なら、LU 分解の演習もするが、今回は、時間数が足りないのでパス == 今後の予定 07/01 本日 ヤコビ方 07/08 次回 SR 法 07/15 まとめ ( 30 分程度 ) + 前期テスト == LU 法は、本質的にガウスの消去法と同じ ( 全然、変っていない ) # 数値計算の世界はまともな手法だけが残っている 直接法は、ガウスの消去法に付きる # 近年、誤差解析をしたウイルソンくらいしかない # 教科書にのってる、プログラムを全部試した # => どれも同じ == # 方程式系のよさ # f(x) = 0 # => 解を確認できる # => 間違っていたら修正すればよい !! # 直接法以外の方法の可能性 ( 反復方法 ) ## cf. 微分方程式 ## => 解を次々と利用してゆく ## 確認の方法がない.. 反復方法 ( ヤコビ、ガウス=ザイデル ) 解を推定し、誤っていたら修正する方法 1) 修正する方法が違えば、手法も異なる 2) 初期値をどうするかも問題 2-a) ゼロ行列にする => 時々解がでないことがある 2-b) b を利用する # なぜ ( 直接法があるのに.. ) 間接法が必要か ? # => 計算量が問題 # A) そもそも、1/3 n^3 がつらい # n = 10^8 の問題が解きたい !! # B) 微分方程式を解きたい # 係数行列がスパースになる # => この性質を利用して更に、サボりたい # 微分方程式の差分化 # \frac{d^2}{dx^2} # = \frac{ u(x+\delta x) - 2u(x) + u(x-\delta x) ) }{(delta x)^2} # \frac{d^2(x,y)}{dx^2} +\frac{d^2(x,y)}{dy^2} + # = .. # 行列がバンド行列の形になる # ガウスの消去法を行うステップで 0 の要素が多いのに # 0 ではなくなってしまう # => その部分の記憶容量が足りない # cf. 地球シミュレータ # 40T の速度, 10T のメモリ # ( (10^8)^3 = 10^24 ) には足りない ) # 一般に n のサイズの問題では、 # ガウスの消去法 # n^2 # 反復法 # Ax だけが必要 ( 0 でない要素だけを覚えればよい ) # 1 次元の微分方程式 3n # 2 次元の微分方程式 5n # 3 次元の微分方程式 7n # なので、n が大きくなると大変 !! # ( メモリだけでなく、計算量も馬鹿にならない ) ## [注意] 反復法は、収束するかどうかで、計算回数(量) が全然 ## 違うので、「計算全体」としては、量が小さいとはいえない ## ただ、一度の計算ではスパース性が効いてくる ### 仮に収束性が悪くても、アーキテクチャの影響もありうる ### もし、直接法がメモリにはいらず ( ディスク 10^{-2}を利用する ) ### 間接法がメモりに入る ( メモリのみ 10^{-7} ) だと、 ### その差 ( 10^5 違う ) の係数がかかる.. # メモリの観点から、巨大問題は、間接法しか利用できない # でも、精度がわるい ( そもそも、収束するかどうかも.. ) # ハイブリッド ( 両方を組合せる ) という手も # 考えられている ## 本当は、ガウスの消去法がよい、でも、大きな次元スパースな問題は、 ## の場合に限っては、間接法をつかわざるをえない == # 偏微分方程式系から出る係数行列は、対角要素の絶対値が大きい # それを利用して.. [ヤコビ法] A x = b ( A = { a_ij } ) に対して x_i = \frac{1}{a_{ii}}( b_i - ( ai x - a_ii xi ) ) 式を変更し、以下、 x_i^{(0)} = 0 と初期値を定め、 x_i^{(j+1)} = \frac{1}{a_{ii}}( b_i - ( ai x^{(j)} - a_ii xi^{(j)} ) ) として、反復する。 [ガウスザイデル法] x_i^{(j+1)} = \frac{1}{a_{ii}}( b_i - ( ai x^{(k)} - a_ii xi^{(k)} ) ) j+1 i < j k = j i > j # 式は複雑だが、プログラムは簡単 # 計算した結果をとっておき、あとでまとめて更新 # => ヤコビ方 # 計算した結果を直に利用する ( 古い値は捨てる ) # => ガウス=ザイデル法 # => メモリも ( わずかながら. ) 良い ## ヤコビ法のプログラムを間違えると自動的にガウス=ザイデル法に ## なることがある。 # 一般に、ガウス=ザイデル法の方が、収束性がよく精度もよい # ( 1.5 倍位.. ) # but.. 並列化が難しい ( ヤコビ法が復活.. ) ## リソースを無駄にせず、やるか、それとも無駄をしても速度を重 ## 視するか ? ( 並列計算機時代の新しい課題 ) == [本日の演習] ガウス=ザイデル # 余裕のある人は、ヤコビ法もやる # => 収束性の差を確認しよう。 係数行列は text p.58 のものを利用する # 直接法は、適当な係数でも、答ができる # 間接法は、適当な係数だと、収束しないことがある !! == 7/15 はテスト 山をかけてください 予め、出す問題を理解して来て欲しいから # 講義に出ていて、ノートをみれば解るでしょう 先輩に聞いて問題を聞こう !! ( 7, 8 割は当るはず.. ) => 答は聞くな !! # 答が正しいという保障はない # 例年、間違った答を覚えて来る人が沢山いて、 # 共倒れ現象が起る テストでは、プログラムを出すことはない cf. 「何々法で、重要なポイントは ?」