Projection法でls_arnoldi3_fom.eq1.gifとする.

ls_arnoldi3_fom.eq2.gif

ここで,ls_arnoldi3_fom.eq3.gifである. この条件を用いる方法としてFOMについて説明する.

近似解ls_arnoldi3_fom.eq4.gifを部分空間ls_arnoldi3_fom.eq5.gifから探索する. 探索するときの条件はls_arnoldi3_fom.eq6.gifから,

ls_arnoldi3_fom.eq7.gif

いま,ls_arnoldi3_fom.eq8.gifとし,Arnoldi法で得られるKrylov部分空間の正規直交ベクトルを 並べた行列ls_arnoldi3_fom.eq9.gifにより,

ls_arnoldi3_fom.eq10.gif

とすると,ls_arnoldi3_fom.eq4.gifは以下のように表される.

ls_arnoldi3_fom.eq11.gif

ls_arnoldi3_fom.eq4.gifのときの残差は,

ls_arnoldi3_fom.eq12.gif

ls_arnoldi3_fom.eq4.gifが解ならばこの式は0となるので,

ls_arnoldi3_fom.eq13.gif

両辺にls_arnoldi3_fom.eq14.gifを掛ける.

ls_arnoldi3_fom.eq15.gif

Arnoldi法より実行列の場合,ls_arnoldi3_fom.eq16.gifなので,

ls_arnoldi3_fom.eq17.gif

ls_arnoldi3_fom.eq18.gifをこの式から計算し,ls_arnoldi3_fom.eq19.gifに代入することで,近似解が得られる. ls_arnoldi3_fom.eq20.gifとし,上式をさらに変形する.

ls_arnoldi3_fom.eq21.gif

まとめると,

ls_arnoldi3_fom.eq22.gif

この式に基づき近似解を求める方法がFOM(Full Orthogonalization Method)である.

FOMの手順

FOMでls_arnoldi3_fom.eq23.gifの近似解を求めるアルゴリズムを以下に示す.

ls_arnoldi3_fom.eq24.gifを計算
ls_arnoldi3_fom.eq25.gifの行列ls_arnoldi3_fom.eq26.gifの各要素を0で初期化
for(j = 1,2,...,m){
  ls_arnoldi3_fom.eq27.gif
  for(i = 1,2,...,j)\{
    ls_arnoldi3_fom.eq28.gif
    ls_arnoldi3_fom.eq29.gif
  }
  ls_arnoldi3_fom.eq30.gif
  if(ls_arnoldi3_fom.eq31.gif) m=jとしてループを抜ける
  ls_arnoldi3_fom.eq32.gif
}
ls_arnoldi3_fom.eq33.gif
ls_arnoldi3_fom.eq19.gif

ls_arnoldi3_fom.eq34.gifのループの中は修正グラム・シュミット法を用いたArnoldi法(Arnoldi法#pc5d0806参照)そのものである. ここではls_arnoldi3_fom.eq35.gifが0かどうかで反復を抜ける判断をしているが, 実際には近似解に対する残差ls_arnoldi3_fom.eq36.gifの大きさで判断したい. かつなるべく残差を計算するのにコストは掛けたくない. そのため,以下の式を用いる.

ls_arnoldi3_fom.eq37.gif

この式の導出を以下に示す.

まず,残差ベクトルは,

ls_arnoldi3_fom.eq38.gif

と表される.いま,ls_arnoldi3_fom.eq39.gifであり, また,Arnoldi法より,ls_arnoldi3_fom.eq40.gifなので,

ls_arnoldi3_fom.eq41.gif

ls_arnoldi3_fom.eq42.gifより,

ls_arnoldi3_fom.eq43.gif

よって,残差の大きさは以下となる.

ls_arnoldi3_fom.eq44.gif

この式を使って残差を求め,それを反復終了条件とすればよい.

FOMの拡張

FOMの計算コストはls_arnoldi3_fom.eq45.gifである.計算コストを減らすためにいくつかの手法が提案されている. ここでは,FOMの拡張であるFOM(m),IOM,について簡単に紹介する(概要だけ).

FOM(m)

FOMの拡張の一つでRestarted FOMである.FOM(m)のアルゴリズムは,

  1. m=1と設定
  2. ls_arnoldi3_fom.eq46.gifからFOMでls_arnoldi3_fom.eq4.gifを計算
  3. ls_arnoldi3_fom.eq47.gifとして,2に戻る.

mが設定した上限値に達するか,残差が十分小さくなるまで,2,3を繰り返す.

IOM

グラム・シュミット法において,ls_arnoldi3_fom.eq48.gifの反復を

ls_arnoldi3_fom.eq49.gif

と変更する.これをincomplete orthogonalizationといい,これを用いたFOMを IOM(Incomplete Orthogonalization Method)という.

DIOM

IOMにおいて問題となるのは,ls_arnoldi4_diom.eq1.gifの計算ですべてのls_arnoldi4_diom.eq2.gifが必要になることである. そこで,ls_arnoldi4_diom.eq3.gifls_arnoldi4_diom.eq4.gifからでなく,ls_arnoldi4_diom.eq5.gifから計算するように修正したのがDIOM(Direct IOM)である.

IOMでi=j-k+1〜jとしたことで,ls_arnoldi4_diom.eq6.gifは値を持つ部分の幅がk+1な行列になる.以下にm=5, k=3の場合の例を示す.

ls_arnoldi4_diom.eq7.gif

まずls_arnoldi4_diom.eq8.gifとLU分解する. ピボッティングがなければ,分解した行列は以下のようになる.

ls_arnoldi4_diom.eq9.gif

DIOMではls_arnoldi4_diom.eq5.gifからls_arnoldi4_diom.eq10.gifを求める形にする. FOMのls_arnoldi4_diom.eq3.gifに関する式は,

ls_arnoldi4_diom.eq11.gif

であり,ls_arnoldi4_diom.eq6.gifをLU分解した行列で置き換える.

ls_arnoldi4_diom.eq12.gif

いま,

ls_arnoldi4_diom.eq13.gif

とおくと以下のように変形できる.

ls_arnoldi4_diom.eq14.gif

ls_arnoldi4_diom.eq15.gifls_arnoldi4_diom.eq16.gifそれぞれの中身について考えてみる.

  • ls_arnoldi4_diom.eq15.gifについて
    ls_arnoldi4_diom.eq15.gifの式はls_arnoldi4_diom.eq17.gifとなり,ls_arnoldi4_diom.eq15.gifの各列をls_arnoldi4_diom.eq18.gifとすると,ls_arnoldi4_diom.eq19.gifが上三角行列であることから, ls_arnoldi4_diom.eq20.gifはガウスの消去法の後退代入を使って以下のように計算される.
    ls_arnoldi4_diom.eq21.gif
    よって,ls_arnoldi4_diom.eq20.gifはm-1までのls_arnoldi4_diom.eq18.gifを使って算出される.ここでls_arnoldi4_diom.eq15.gifを以下のように書くことにする.
    ls_arnoldi4_diom.eq22.gif
  • ls_arnoldi4_diom.eq16.gifについて
    ls_arnoldi4_diom.eq16.gifについてもls_arnoldi4_diom.eq23.gifであり,例えば,m=5の場合,
    ls_arnoldi4_diom.eq24.gif
    よって,
    ls_arnoldi4_diom.eq25.gif
    ただし,ls_arnoldi4_diom.eq26.gifである. このようにls_arnoldi4_diom.eq27.gifもm-1以下のls_arnoldi4_diom.eq28.gifから計算される.こちらもls_arnoldi4_diom.eq15.gifの時と同様に以下のように書く.
    ls_arnoldi4_diom.eq29.gif

これらのことを使って式を書き換えると,

ls_arnoldi4_diom.eq30.gif

ここで,ls_arnoldi4_diom.eq31.gifls_arnoldi4_diom.eq5.gifなので,

ls_arnoldi4_diom.eq32.gif

添付ファイル: filels_arnoldi3_fom.eq40.gif 539件 [詳細] filels_arnoldi3_fom.eq28.gif 541件 [詳細] filels_arnoldi3_fom.eq24.gif 558件 [詳細] filels_arnoldi3_fom.eq18.gif 584件 [詳細] filels_arnoldi3_fom.eq22.gif 556件 [詳細] filels_arnoldi3_fom.eq25.gif 516件 [詳細] filels_arnoldi3_fom.eq51.gif 543件 [詳細] filels_arnoldi3_fom.eq7.gif 539件 [詳細] filels_arnoldi3_fom.eq31.gif 604件 [詳細] filels_arnoldi3_fom.eq34.gif 538件 [詳細] filels_arnoldi3_fom.eq12.gif 592件 [詳細] filels_arnoldi3_fom.eq47.gif 581件 [詳細] filels_arnoldi3_fom.eq37.gif 516件 [詳細] filels_arnoldi3_fom.eq38.gif 545件 [詳細] filels_arnoldi3_fom.eq45.gif 524件 [詳細] filels_arnoldi3_fom.eq26.gif 503件 [詳細] filels_arnoldi3_fom.eq16.gif 578件 [詳細] filels_arnoldi3_fom.eq43.gif 527件 [詳細] filels_arnoldi3_fom.eq46.gif 531件 [詳細] filels_arnoldi3_fom.eq49.gif 580件 [詳細] filels_arnoldi3_fom.eq42.gif 538件 [詳細] filels_arnoldi3_fom.eq36.gif 543件 [詳細] filels_arnoldi3_fom.eq17.gif 604件 [詳細] filels_arnoldi3_fom.eq20.gif 501件 [詳細] filels_arnoldi3_fom.eq39.gif 573件 [詳細] filels_arnoldi3_fom.eq48.gif 486件 [詳細] filels_arnoldi3_fom.eq2.gif 487件 [詳細] filels_arnoldi3_fom.eq52.gif 544件 [詳細] filels_arnoldi3_fom.eq13.gif 894件 [詳細] filels_arnoldi3_fom.eq19.gif 524件 [詳細] filels_arnoldi3_fom.eq4.gif 526件 [詳細] filels_arnoldi3_fom.eq44.gif 476件 [詳細] filels_arnoldi3_fom.eq5.gif 507件 [詳細] filels_arnoldi3_fom.eq30.gif 531件 [詳細] filels_arnoldi3_fom.eq32.gif 534件 [詳細] filels_arnoldi3_fom.eq15.gif 596件 [詳細] filels_arnoldi3_fom.eq21.gif 497件 [詳細] filels_arnoldi3_fom.eq23.gif 489件 [詳細] filels_arnoldi3_fom.eq9.gif 509件 [詳細] filels_arnoldi3_fom.eq1.gif 564件 [詳細] filels_arnoldi3_fom.eq8.gif 503件 [詳細] filels_arnoldi3_fom.eq33.gif 489件 [詳細] filels_arnoldi3_fom.eq3.gif 566件 [詳細] filels_arnoldi3_fom.eq27.gif 539件 [詳細] filels_arnoldi3_fom.eq50.gif 553件 [詳細] filels_arnoldi3_fom.eq35.gif 473件 [詳細] filels_arnoldi3_fom.eq6.gif 575件 [詳細] filels_arnoldi3_fom.eq10.gif 458件 [詳細] filels_arnoldi3_fom.eq14.gif 540件 [詳細] filels_arnoldi3_fom.eq29.gif 505件 [詳細] filels_arnoldi3_fom.eq41.gif 538件 [詳細] filels_arnoldi3_fom.eq11.gif 540件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-07-12 (木) 15:26:10 (2995d)