クリロフ部分空間

共役勾配法の計算手順において,残差ls_krylov.eq1.gifを用いた. この残差はどこから来たのかを考える. まず,ls_krylov.eq2.gif

ls_krylov.eq3.gif

とかける.これを反復式とすると,

ls_krylov.eq4.gif

これで残差ベクトルがでてきた. さらに残差ベクトル間の関係を調べる.

ls_krylov.eq5.gif

よって,

ls_krylov.eq6.gif

同様にls_krylov.eq7.gifに関しても,

ls_krylov.eq8.gif

これらの式から,ls_krylov.eq9.gifls_krylov.eq10.gifls_krylov.eq11.gif の線形結合で表されることがわかる.これを式にすると,

ls_krylov.eq12.gif
ls_krylov.eq13.gif

ここで,ls_krylov.eq14.gifはベクトルls_krylov.eq15.gifの線形結合の集合で表される部分空間であり, 上式のような部分空間をクリロフ部分空間と呼ぶ.

ls_krylov.eq16.gif

ls_krylov.eq17.gifの次元nは近似解を求めるための反復ごとに増えていく. そして,クリロフ部分空間内の任意の点はls_krylov.eq18.gifと書ける. ここでls_krylov.eq19.gifは次元がm-1以下の多項式を表している. つまり,クリロフ部分空間内の点はAに関するm-1次以下の多項式とls_krylov.eq20.gifの積の形で書き表せる.

ls_krylov.eq21.gifからクリロフ部分空間ls_krylov.eq17.gifの中を探索することで解を得る非定常な反復解法のことを クリロフ部分空間法と呼ぶ.共役勾配法もクリロフ部分空間法のひとつである (ヤコビ反復やガウス・ザイデルは定常な反復解法).

クリロフ部分空間法としては他に,

  • 双共役勾配法(Bi-Conjugate Gradient method : BiCG法)
  • 安定化双共役勾配法(Bi-Conjugate Gradient STABilized method : BiCGSTAB法)
  • 自乗共役勾配法(Conjugate Gradiate Squared method : CGS法)
  • 共役残差法(Conjugate Residual method : CR法)
  • 一般化共役残差法(Generalized Conjugate Residual method : GCR法)
  • 一般化最小残差法(Generalized Minimal RESidual method : GMRES法) などが提案されている.

アーノルディ法

アーノルディ法(Arnoldi's Method)は非対称行列のKrylov部分空間における正規直交基底を求める方法である. ちなみに対称行列に限定したLanczos法もある.

Arnoldi法のアルゴリズムを以下に示す.

任意のベクトルls_arnoldi.eq1.gifを設定(ただしls_arnoldi.eq2.gif)
for(ls_arnoldi.eq3.gif){
  ls_arnoldi.eq4.gif
  ls_arnoldi.eq5.gif
  ls_arnoldi.eq6.gif
  もし,ls_arnoldi.eq7.gifなら反復終了
  ls_arnoldi.eq8.gif
}

ここでAは対象となる非対称行列である.

wの計算手順をまとめて書くと,

ls_arnoldi.eq9.gif

となり,ls_arnoldi.eq10.gifなので, ls_arnoldi.eq11.gifグラム・シュミットの直交化法でベクトルls_arnoldi.eq12.gifからls_arnoldi.eq13.gifに 直交するベクトルを求めていることになる. そのため, Arnoldi法でj=mまで計算された場合, ls_arnoldi.eq14.gifはKrylov部分空間の正規直交基底

ls_arnoldi.eq15.gif

を構成する.

さて,アルゴリズムより,

ls_arnoldi.eq16.gif

であり,これを変形すると,

ls_arnoldi.eq17.gif

ここでls_arnoldi.eq18.gifである. この式がどのような形になっているのかを確かめるために, m=3の場合で,ls_arnoldi.eq19.gifを列とする行列を使って書き下してみる.

ls_arnoldi.eq20.gif

上記の式はこのように書ける.これをm次元に一般化する. ls_arnoldi.eq21.gifの行列ls_arnoldi.eq22.gifを使うと,

ls_arnoldi.eq23.gif

ls_arnoldi.eq24.gifは基本ベクトルである. ls_arnoldi.eq25.gifはヘッセンベルグ標準形と呼ばれる形式になっており, 行列AをKrylov部分空間に射影した行列となっている. ls_arnoldi.eq26.gifは直交行列なので,この式の両辺にls_arnoldi.eq27.gifを掛けると,

ls_arnoldi.eq28.gif

となる.

ちなみに,Aとls_arnoldi.eq25.gifの固有値は同じとなるので, この方法は固有値計算にも用いられる.

修正グラム・シュミット法を用いたArnoldi法

修正グラムシュミット法(グラム・シュミットの直交化法参照)を使ったArnoldi法のアルゴリズムを以下に示す.

任意のベクトルls_arnoldi.eq1.gifを設定(ただしls_arnoldi.eq2.gif)
for(j = 1,2,...,m){
  ls_arnoldi.eq29.gif
  for(i = 1,2,...,j){
    ls_arnoldi.eq30.gif
    ls_arnoldi.eq31.gif
  }
  ls_arnoldi.eq32.gif
  if(ls_arnoldi.eq33.gif) 反復終了
  ls_arnoldi.eq8.gif
}

丸め誤差がなければ通常のArnoldi法と結果は同じとなる.

Projection法

線形システム(Aはls_arnoldi2_projection.eq1.gif)

ls_arnoldi2_projection.eq2.gif

ls_arnoldi2_projection.eq3.gifステップ目の近似解ls_arnoldi2_projection.eq4.gifを初期値ls_arnoldi2_projection.eq5.gifから求める. いま,2つの部分空間ls_arnoldi2_projection.eq6.gifを考える. そして,

  • ls_arnoldi2_projection.eq7.gif
  • ls_arnoldi2_projection.eq8.gif

として,解を探索するのがProjection法である. Kは解候補が含まれるsearch subspaceで, Lは近似解を得るための制約空間(subspace of constraints)である. Projection法では残差ベクトルls_arnoldi2_projection.eq9.gifがLに対して垂直になるようにする. この条件をPetrov-Galerkin conditionという.

ここで,ls_arnoldi2_projection.eq10.gifとすると,

ls_arnoldi2_projection.eq11.gif

よって,

ls_arnoldi2_projection.eq12.gif

2つめの条件により,ls_arnoldi2_projection.eq13.gifはLに垂直なベクトルとなる. これが基本的なProjection法のステップである. また,K,Lはステップごとに適切なものを選ぶ必要がある.

Krylov部分空間によるProjection法

共役勾配法を代表的なものとして,多くの方法が Projection法において,K,LにKrylov部分空間を用いたものになっている. Krylov部分空間をls_arnoldi2_projection.eq14.gifとすると, これらの方法は大まかに以下のように分類される.

  • ls_arnoldi2_projection.eq15.gif : FOMなど
  • ls_arnoldi2_projection.eq16.gif : Lanczos, BiCG, CGS, QMRなど
  • ls_arnoldi2_projection.eq17.gif : GMRESなど

FOM

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

参考文献

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

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-07-12 (木) 10:56:23