連立1次方程式：Krylov部分空間法

---

---

クリロフ部分空間†

共役勾配法の計算手順において，残差を用いた． この残差はどこから来たのかを考える． まず，は

とかける．これを反復式とすると，

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

よって，

同様にに関しても，

これらの式から，，は の線形結合で表されることがわかる．これを式にすると，

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

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

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

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

• 双共役勾配法(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法のアルゴリズムを以下に示す．

任意のベクトルを設定(ただし)
for(){
　　
　　
　　
　　もし，なら反復終了
　　
}

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

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

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

を構成する．

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

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

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

上記の式はこのように書ける．これをm次元に一般化する． の行列を使うと，

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

となる．

ちなみに，Aとの固有値は同じとなるので， この方法は固有値計算にも用いられる．

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

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

任意のベクトルを設定(ただし)
for(j = 1,2,...,m){
　　
　　for(i = 1,2,...,j){
　　　　
　　　　
　　}
　　
　　if() 反復終了
　　
}

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

Projection法†

線形システム(Aは)

のステップ目の近似解を初期値から求める． いま，2つの部分空間を考える． そして，

• 
• 

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

ここで，とすると，

よって，

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

Krylov部分空間によるProjection法†

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

• : FOMなど
• : Lanczos, BiCG, CGS, QMRなど
• : GMRESなど

FOM†

Projection法でとする．

ここで，である． この条件を用いる方法としてFOMについて説明する．

近似解を部分空間から探索する． 探索するときの条件はから，

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

とすると，は以下のように表される．

のときの残差は，

が解ならばこの式は0となるので，

両辺にを掛ける．

Arnoldi法より実行列の場合，なので，

をこの式から計算し，に代入することで，近似解が得られる． とし，上式をさらに変形する．

まとめると，

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

FOMの手順†

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

を計算
の行列の各要素を0で初期化
for(j = 1,2,...,m){
　　
　　for(i = 1,2,...,j)\{
　　　　
　　　　
　　}
　　
　　if() m=jとしてループを抜ける
　　
}

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

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

まず，残差ベクトルは，

と表される．いま，であり， また，Arnoldi法より，なので，

より，

よって，残差の大きさは以下となる．

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

FOMの拡張†

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

FOM(m)†

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

1. m=1と設定
2. からFOMでを計算
3. として，2に戻る．

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

IOM†

グラム・シュミット法において，の反復を

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

DIOM†

IOMにおいて問題となるのは，の計算ですべてのが必要になることである． そこで，をからでなく，から計算するように修正したのがDIOM(Direct IOM)である．

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

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

DIOMではからを求める形にする． FOMのに関する式は，

であり，をLU分解した行列で置き換える．

いま，

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

，それぞれの中身について考えてみる．

• について
  の式はとなり，の各列をとすると，が上三角行列であることから， はガウスの消去法の後退代入を使って以下のように計算される．
  
  よって，はm-1までのを使って算出される．ここでを以下のように書くことにする．
  
• について
  についてもであり，例えば，m=5の場合，
  
  よって，
  
  ただし，である． このようにもm-1以下のから計算される．こちらもの時と同様に以下のように書く．
  

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

ここで，はなので，

参考文献†

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