対称行列の正規直交基底を算出するLanczos法とその線形方程式への適用,そして,共役勾配法のアルゴリズムの導出について説明する.



Lanczos法

Arnoldi法では非対称行列を扱ったが,Aが対称行列だった場合, アルゴリズムを簡素化できる. これがランチョス(Lanczos)法である.

まず,Arnoldi法でAが実対称行列だったとき, ls_lanzcos1.eq1.gifも対称行列になる. さらにls_lanzcos1.eq2.gifはヘッセンベルグ行列であり,ls_lanzcos1.eq3.gif,ただし i > j+1 である. よって,ls_lanzcos1.eq2.gifが対称行列ならば,i < j-1 でもls_lanzcos1.eq3.gifとなる. つまり,以下のような三重対角行列(対角成分とその左右のみに値がある行列)となる.

ls_lanzcos1.eq4.gif

ここで,ls_lanzcos1.eq5.gifとした.

ls_lanzcos1.eq6.gifの代わりにls_lanzcos1.eq7.gifを使って, Arnoldi法(修正グラム・シュミットを用いたもの)を書き換える.

今,i < j-1でls_lanzcos1.eq6.gifが0となることから,i=j-1,jについてのみ考えればよい. i=j-1では,ls_lanzcos1.eq8.gifなので,ls_lanzcos1.eq9.gifと置き換えることができる. ただし,ls_lanzcos1.eq10.gifと置いておく. 次に,i=jでは,ls_lanzcos1.eq11.gifなので,ls_lanzcos1.eq12.gifと置き換えることができる. 最後にls_lanzcos1.eq13.gifと置き換え,ls_lanzcos1.eq14.gifは次の反復においてls_lanzcos1.eq15.gifとして用いられる. これらのことを適用するとアルゴリズムは以下のようになる.

任意のベクトルls_lanzcos1.eq16.gifを設定(ただしls_lanzcos1.eq17.gif)
ls_lanzcos1.eq18.gifを設定
for(j = 1,2,...,m){
  ls_lanzcos1.eq9.gif
  ls_lanzcos1.eq19.gif
  ls_lanzcos1.eq12.gif
  ls_lanzcos1.eq20.gif
  if(ls_lanzcos1.eq21.gif) 反復終了
  ls_lanzcos1.eq22.gif
}

これがLanczos法である.

共役勾配法のアルゴリズムの導出

Lanczos法から共役勾配法を導くことができる. 共役勾配法は,そのため,FOMと同じくls_lanzcos2_ls.eq1.gifとしたProjection法となる. ここでは,Lanczos法,Direct版のLanczos法を説明し,そこから導き出せる関係式,そして, その関係式を使って共役勾配法のアルゴリズムを導出するまでを説明する.

Lanczos法を使った線形システムの解法

線形システムls_lanzcos2_ls.eq2.gif(ただし,Aは対称行列)について,初期近似値ls_lanzcos2_ls.eq3.gifが与えられたとき, mステップ目の近似値ls_lanzcos2_ls.eq4.gifはFOMと同様に以下のようになる.

ls_lanzcos2_ls.eq5.gif

ここで,ls_lanzcos2_ls.eq6.gifが三重対角行列(tridiagonal matrix)になるのでls_lanzcos2_ls.eq7.gifと置き換えている. また,ls_lanzcos2_ls.eq8.gifである.

この式からLanczos法を使った線形システムの解法アルゴリズムは以下となる.

ls_lanzcos2_ls.eq9.gifを計算
ls_lanzcos2_ls.eq10.gifを設定
for(j = 1,2,...,m){
  ls_lanzcos2_ls.eq11.gif
  ls_lanzcos2_ls.eq12.gif
  ls_lanzcos2_ls.eq13.gif
  ls_lanzcos2_ls.eq14.gif
  if(ls_lanzcos2_ls.eq15.gif) 反復終了
  ls_lanzcos2_ls.eq16.gif
}
ls_lanzcos2_ls.eq17.gifで構成される三重対角行列ls_lanzcos2_ls.eq7.gifを設定
ls_lanzcos2_ls.eq18.gif
ls_lanzcos2_ls.eq19.gif

Direct版のLanczos法

FOMからDIOMを導き出したときと同様にDirect版のLanczos法を考える. Lanczos法によるls_lanzcos3_direct.eq1.gifは三重対角行列なので,DIOMにおけるk=2の場合(値がある部分の幅が3)と考えられる. m=5のとき,ls_lanzcos3_direct.eq1.gifのLU分解は以下のようになる.

ls_lanzcos3_direct.eq2.gif

DIOMと同様に,ls_lanzcos3_direct.eq3.gifとおくと,

ls_lanzcos3_direct.eq4.gif

となる.

  • ls_lanzcos3_direct.eq5.gifに関して
    DIOMのls_lanzcos3_direct.eq5.gifの最後の列ls_lanzcos3_direct.eq6.gifに関する式,
    ls_lanzcos3_direct.eq7.gif
    において,k=2なので,i=m-1のときだけ考えればよい.つまり,
    ls_lanzcos3_direct.eq8.gif
  • ls_lanzcos3_direct.eq9.gifに関して
    ls_lanzcos3_direct.eq10.gif
    となるので,
    ls_lanzcos3_direct.eq11.gif
    とする.

最終的にDIOMと同様にls_lanzcos3_direct.eq13.gifからls_lanzcos3_direct.eq14.gifを求める.

ls_lanzcos3_direct.eq15.gif

この式によりls_lanzcos3_direct.eq14.gifを更新していくのがDirect版のLanczos法である.

Direct版でのls_lanzcos3_direct.eq16.gifはガウス消去法のステップから,

ls_lanzcos3_direct.eq17.gif

により求めることができる.

ls_lanzcos3_direct.eq18.gifを計算
ls_lanzcos3_direct.eq19.gifを設定
for(j = 1,2,...){
  ls_lanzcos3_direct.eq20.gif
  ls_lanzcos3_direct.eq21.gif
  if(ls_lanzcos3_direct.eq22.gif) ls_lanzcos3_direct.eq23.gif
  ls_lanzcos3_direct.eq24.gif
  ls_lanzcos3_direct.eq25.gif
  ls_lanzcos3_direct.eq26.gif
  if(収束判定) 反復終了
  ls_lanzcos3_direct.eq27.gif
  ls_lanzcos3_direct.eq28.gif
  ls_lanzcos3_direct.eq29.gif
}

直交・共役関係

DIOMで得られた残差ベクトルls_lanzcos4_direct2.eq1.gifls_lanzcos4_direct2.eq2.gifについて考えてみる. FOMの収束判定で導いたように,

ls_lanzcos4_direct2.eq3.gif

であり,ls_lanzcos4_direct2.eq1.gifls_lanzcos4_direct2.eq4.gifのスカラー倍になっている. ls_lanzcos4_direct2.eq5.gifは直交系であるので,

ls_lanzcos4_direct2.eq6.gif

また,ls_lanzcos4_direct2.eq2.gifに関しては,以下の関係が成り立つ.

ls_lanzcos4_direct2.eq7.gif

これはAに関する共役を意味しており,A-共役 (A-conjugate)と呼ばれる.Aが単位行列ならばrと同じく直交になる.

これに関する証明を以下に示す.

ls_lanzcos4_direct2.eq8.gifであるならば,ls_lanzcos4_direct2.eq9.gifは対角行列になるはずである. ls_lanzcos4_direct2.eq10.gifなので,

ls_lanzcos4_direct2.eq11.gif

Aは対称行列なので,ls_lanzcos4_direct2.eq9.gifも対称行列になる. よって,ls_lanzcos4_direct2.eq12.gifも対称行列である.また,ls_lanzcos4_direct2.eq13.gifは下三角行列(上三角行列の逆行列もまた上三角行列なので)である. 対称な下三角行列とは要するに対角行列である.よって,ls_lanzcos4_direct2.eq9.gifは対角行列であり, ls_lanzcos4_direct2.eq2.gifのA-共役が成り立つ.

共役勾配アルゴリズムの導出

ls_lanzcos5_cg.eq1.gifls_lanzcos5_cg.eq2.gifに関する関係式を使って共役勾配法のアルゴリズムを導き出してみる. まず,基本的なProjection法のステップを思いだそう.

ls_lanzcos5_cg.eq3.gif

ここで,ls_lanzcos5_cg.eq4.gifである. これをDirect版Lanczosアルゴリズムに当てはめる.ls_lanzcos5_cg.eq5.gifとすると,

ls_lanzcos5_cg.eq6.gif

これまで,ls_lanzcos5_cg.eq7.gifから始まっていたが,ここではより一般的なls_lanzcos5_cg.eq8.gifからのスタートにしてあるので注意. また,ls_lanzcos5_cg.eq9.gifと書き換えている. 次にls_lanzcos5_cg.eq2.gifの更新式を考える. Direct版のLanczosではls_lanzcos5_cg.eq10.gifであり, ls_lanzcos5_cg.eq1.gifls_lanzcos5_cg.eq11.gifのスカラー倍になっているので,残差ベクトルls_lanzcos5_cg.eq1.gifと前ステップのls_lanzcos5_cg.eq2.gifを使って,

ls_lanzcos5_cg.eq12.gif

と書ける. なお,式中のls_lanzcos5_cg.eq13.gifls_lanzcos5_cg.eq14.gifの要素として使っていたls_lanzcos5_cg.eq13.gifとは別物なので注意.

前節で述べたls_lanzcos5_cg.eq1.gifls_lanzcos5_cg.eq2.gifの直交,共役関係を使ってls_lanzcos5_cg.eq13.gifを算出することで,共役勾配法のアルゴリズムが得られる.

  • ls_lanzcos5_cg.eq15.gifの導出
    ls_lanzcos5_cg.eq1.gifに関する関係式よりls_lanzcos5_cg.eq16.gifなので,
    ls_lanzcos5_cg.eq17.gif
    これをls_lanzcos5_cg.eq15.gifについて解くと以下となる.
    ls_lanzcos5_cg.eq18.gif
    分母について,ls_lanzcos5_cg.eq2.gifの性質ls_lanzcos5_cg.eq19.gifを用いると,
    ls_lanzcos5_cg.eq20.gif
    となる.よって,
    ls_lanzcos5_cg.eq21.gif
  • ls_lanzcos5_cg.eq22.gifの導出
    ls_lanzcos5_cg.eq23.gifの関係を用いる.まず,
    ls_lanzcos5_cg.eq24.gif
    であり,よって以下のようにls_lanzcos5_cg.eq22.gifに関する式が求められる.
    ls_lanzcos5_cg.eq25.gif
    これをさらに変形する. いま,ls_lanzcos5_cg.eq26.gifから, ls_lanzcos5_cg.eq27.gifなので,
    ls_lanzcos5_cg.eq28.gif
    ls_lanzcos5_cg.eq29.gifより,最終的に,
    ls_lanzcos5_cg.eq30.gif

これらの式から,以下の共役勾配法のアルゴリズムが得られる.

ls_lanzcos5_cg.eq31.gifを計算
for(j = 0,1,...){
  ls_lanzcos5_cg.eq32.gif
  ls_lanzcos5_cg.eq33.gif
  ls_lanzcos5_cg.eq26.gif
  if(収束判定) 反復終了
  ls_lanzcos5_cg.eq34.gif
  ls_lanzcos5_cg.eq35.gif
}

参考文献

  • Yousef Saad, Iterative methods for sparse linear systems 2nd ed., SIAM, 2003.

添付ファイル: filels_lanzcos2_ls.eq1.gif 787件 [詳細]

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2022-11-30 (水) 13:48:13