JP2006172137A - 情報処理装置、情報処理方法 - Google Patents

Abstract

【課題】 係数行列が正則かどうか不明である連立一次方程式を解くこと。
【解決手段】 共役勾配法でもって連立一次方程式を解く計算処理が発散すれば（Ｓ１０７）、初期値をX_MINに設定し、共役残差法でもって連立一次方程式を解く計算処理を行う（Ｓ１０８）。
【選択図】 図１

Description

本発明は、連立一次方程式の解を求める為の計算処理技術に関するものである。

様々な分野で解かれる問題の多くは、離散近似等により連立一次方程式に帰着されることが知られている。

これら連立一次方程式の係数行列は通常正則である。従って、Gaussの消去法によって解を得ることが出来る。また、係数行列が大規模でスパースな場合には、Jacobi法、Gauss−Seidel法等の反復法や、共役勾配法が有効であることが知られている（非特許文献１を参照）。対象とする問題に対してどの解法を選択するかによって、必要とされる計算機のスペックや計算速度が異なることから、以下のような技術が開示されている。

特許文献１では、行列がスパースであるかどうか、及び電源の周波数がパラメタ以下かどうかによって、連立一次方程式の解法として、直接法と反復法を選択する方法を提供している。これによってメモリ容量に応じた最適な解法が選択できるとしている。

また特許文献２では、求解時間を予測する工程を設け、その結果を用いて、行列に対する複数の三角分解解法のうちから最適な解法を選択する方法を提供している。これによって、与えられた問題を最速に解くことができるとしている。

しかし、対象とする問題によっては係数行列が特異になる場合がある。このときには、直接法としてのGauss法や、Jacobi法に基づく一連の反復法は破綻する。

一方、共役勾配法(以下CG法と略記)は、係数行列が特異であってもある条件が成り立つ場合には、解を与えることが知られている（非特許文献２を参照）。また、CG法が破綻するような場合でも、共役残差法(以下CR法と略記)は破綻せず最適な解に収束することが知られている（非特許文献３）。
特開平05-073527号公報 特開平11-242664号公報 森、杉原、室田："数値計算、"岩波講座 応用数学[方法２]、岩波書店(1994)。 EF. Kaasschieter : "Preconditioned conjugate gradients for solving singular systems," Journal of Computational and Applied Mathematics, 24, pp.265-275 (1988). 阿部、緒方、杉原、張、三井："特異な係数行列をもつ連立一次方程式に対するCR法の収束性、"日本応用数理学会論文誌、Vol.9, No.1, pp.1-13 (1999).

以上の説明をまとめると、係数行列が正則であれば、直接法、反復法、CG法、CR法の全ての方法が有効である。係数行列が特異である場合、ある条件を満たせば、CG法、CR法が有効である。それ以外の場合は、CR法のみが有効である。

しかしこのように係数行列が正則であるか、特異であるかに応じて、解法が異なる故、係数行列が正則かどうか不明であるような連立一次方程式を解く場合には、事前に解法を決めることができない。

本発明は以上の問題に鑑みて成されたものであり、係数行列が正則かどうか不明である連立一次方程式を解くための計算技術を提供することを目的とする。

本発明の目的を達成するために、例えば本発明の情報処理装置は以下の構成を備える。

即ち、連立一次方程式の解データを求める為の第１の計算方法に従った処理を行う第１の計算手段と、
前記第１の計算手段が計算処理を反復する毎に得られる計算結果のデータに基づいて、当該計算処理が発散しているか否かを判断する判断手段と、
前記判断手段が発散していると判断した場合には、前記第１の計算手段による計算処理の代わりに、前記連立一次方程式の解データを求める為の第２の計算方法に従った処理を行う第２の計算手段と
を備えることを特徴とする。

本発明の目的を達成するために、例えば本発明の情報処理方法は以下の構成を備える。

即ち、連立一次方程式の解データを求める為の第１の計算方法に従った処理を行う第１の計算工程と、
前記第１の計算工程で計算処理を反復する毎に得られる計算結果のデータに基づいて、当該計算処理が発散しているか否かを判断する判断工程と、
前記判断工程で発散していると判断した場合には、前記第１の計算工程による計算処理の代わりに、前記連立一次方程式の解データを求める為の第２の計算方法に従った処理を行う第２の計算工程と
を備えることを特徴とする。

本発明の構成により、係数行列が正則かどうか不明である連立一次方程式を解くことができる。

以下添付図面を参照して、本発明を好適な実施形態に従って詳細に説明する。

図２は、本実施形態に係る情報処理装置として機能するコンピュータの基本構成を示すブロック図である。

２０１はＣＰＵで、ＲＡＭ２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。

２０２は表示装置で、ＣＲＴや液晶画面などにより構成されており、ＣＰＵ２０１による処理結果を文字が画像などでもって表示することができる。

２０３は入力装置で、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。

２０４はＲＡＭで、ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。

２０５はハードディスク装置で、ここにはＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されており、これらの一部もしくは全部はＣＰＵ２０１による制御に従ってＲＡＭ２０４にロードされる。

２０６はＮＩＣ（ネットワークインターフェース）で、外部装置とのデータ通信を行う際のインターフェースとして機能するものであり、本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。

２０７は上述の各部を繋ぐバスである。

次に、上記構成を備える。コンピュータが行う、連立一次方程式の解を求める為の数値計算処理について説明する。

図１は、本実施形態に係るコンピュータが行う、連立一次方程式の解を求めるための数値計算処理のフローチャートである。なお、同図のフローチャートに従った処理をＣＰＵ２０１に実行させるためのプログラムやデータはハードディスク装置２０５に保存されており、ＣＰＵ２０１の制御に従ってＲＡＭ２０４にロードされ、これを用いてＣＰＵ２０１が処理を行うことで、本実施形態に係るコンピュータは以下説明する処理を行うことになる。

先ず、以下の処理で用いる各変数の為の記憶領域をＲＡＭ２０４内に確保すると共に、解く対象となる連立一次方程式に係るデータ（Ｍ×Ｎの係数行列をＡ、Ｍ×１の変数行列（変数ベクトル）をＸ、Ｍ×１の右辺を示す行列をｂとすると、連立一次方程式はＡＸ＝ｂの形でもって表現できるので、連立一次方程式に係るデータは係数行列Ａのデータ、右辺の行列Ｂのデータを示す）や以下の説明で既知のパラメータとして説明するものをハードディスク装置２０５からＲＡＭ２０４にロードし、さらには、以下の処理で用いる各種の変数を初期化するのであるが、例えば（反復回数を代入するための変数）変数ｔを０に初期化する（ステップＳ１０１）。

次に、共役勾配法に従った処理を１回（１ステップ）実行する（ステップＳ１０２）。そして残差R_tの最小値R_MINを更新する（ステップＳ１０３）。t=0のときは、R_MIN=R₀とする。

次に、R_MINに対する変数ベクトルX_MIN、及び反復回数t_MINのそれぞれが保持する値をＲＡＭ２０４に一時的に記憶する（ステップＳ１０４）。

次に、累積減少量CRD_t、累積増加量CRI_tを求める（ステップＳ１０５）。ステップＳ１０５における処理については詳しくは図３〜６を用いて後述する。

次に、共役勾配法の収束判定条件に基づいて、計算が収束したかどうかを判定し（ステップＳ１０６）、収束していれば本処理を終了する。一方、収束していなければ処理をステップＳ１０７に進める。ここで、収束しているのか否かの判定処理については特に限定するものではないが、例えば反復回数、残差ベクトルのノルム等によって予め決定された方法を用いるようにしても良い。

ステップＳ１０６において収束していないと判断された場合には、以下の式（１）に従った大小比較処理を行う（ステップＳ１０７）。

εCRD_t＜CRI_t （１）
ここでεは経験的に決定される正の実数である。CRD_t、CRI_tについては後述する。そしてこのような大小比較処理を行った結果、εCRD_t≧CRI_tである場合には処理をステップＳ１０２に戻し、変数ｔに１を加算して以降の処理を繰り返す。

一方、εCRD_t＜CRI_tの場合、即ち、共役勾配法による計算処理が発散したと判断される場合には、処理をステップＳ１０８に進め、変数ベクトルの初期値にX_MINを設定し、共役残差法に従って、連立一次方程式の解を求める処理を行う（ステップＳ１０８）。

次に、上記ステップＳ１０５における処理の詳細について説明する。上記CRD_t、CRI_tは以下の式に従って求める。

CRD_t ← CRD_t-1 + max[(R₀-R_t),0] （２）
CRI_t ← CRI_t-1 + max[(R_t-R_MIN),0] （３）
ここで、ｍａｘ（ｘ、ｙ）はｘ、ｙの大きい方の値を出力する関数である。図４は、共役勾配法、共役残差法それぞれにおける残差と反復回数との関係を示す図である。図４に示す如く、縦軸は残差、横軸は反復回数を示す。また、同図において実線の折れ線４０３は、共役勾配法に従った計算処理を行っているときの残差R_tの反復回数に対する変化を示し、二点鎖線４０４は、共役残差法に従った計算処理を行っているときの残差R_tの反復回数に対する変化を示す。また４０５はR_MINであり、この時点での反復回数がt_MIN、変数ベクトルがX_MINである。

４０６は初期残差R₀である。４０７は４０８が現れるまでの一時的なR_MINであり、４０９は４１０が現れるまでの一時的なR_MINである。４１１は共役勾配法に従った計算処理が発散したと判断されたときの反復回数である。即ち、図４は、上記式（１）を満たしている（εCRD_t＜CRI_t）場合のものである。

図３は、ステップＳ１０５における処理のフローチャートである。先ず、残差R_t＜R_MINであるのか否かを判断し（ステップＳ３０１）、残差R_t＜R_MINの場合には処理をステップＳ３０２に進め、以下の式（４），（５）に従った処理を行う（ステップＳ３０２）。

R_MIN = R_t （４）
X_MIN = X_t （５）
これは図４において４０６から４０７の間、４０８から４０９の間、４１０から４０５の間で成り立つ。ステップＳ３０２における処理が完了すると、処理をステップＳ３０４に進める。

一方、ステップＳ３０１における判断処理で残差R_t≧R_MINの場合には処理をステップＳ３０３に進め、以下の式（６）に従った計算処理を行い、CRI_tを更新する（ステップＳ３０３）。

CRI_t ← CRI_t-1 + (R_t - R_MIN) （６）
ステップＳ３０２における処理が完了すると、処理をステップＳ３０４に進める。以上の処理が完了すると次に、R_t＜R₀であるのか否かを判断する（ステップＳ３０４）。R_t≧R₀の場合には本処理、即ちステップＳ１０５における処理を完了し、処理をステップＳ１０６に進める。一方、R_t＜R₀の場合には処理をステップＳ３０５に進め、以下の式（７）に従った計算処理を行い、CRD_tを更新する（ステップＳ３０５）。

CRD_t ← CRD_t-1 + (R₀ - R_t) （７）
なお、CRD_t、CRI_tの更新はそれぞれ式（６）、（７）に従って更新することに限定するものではなく、例えば以下の式（８）、（９）に従って更新するようにしても良い。

CRD_t ← CRD_t-1 + max[(R₀-R_t),0] t ≦ t_MIN （８）
CRI_t ← CRI_t-1 + max[(R_t-R_MIN),0] t_MIN ＜ t （９）
以下では、CRD_t、CRI_tの更新を式（８）、（９）に従って行う場合について説明する。

図６は、共役勾配法、共役残差法それぞれにおける残差と反復回数との関係を示す図である。図６に示す如く、縦軸は残差、横軸は反復回数を示す。また、同図において実線の折れ線６０３は、共役勾配法に従った計算処理を行っているときの残差R_tの反復回数に対する変化を示し、二点鎖線６０４は、共役残差法に従った計算処理を行っているときの残差R_tの反復回数に対する変化を示す。また６０５はR_MINであり、この時点での反復回数がt_MIN、変数ベクトルがX_MINである。

６０６は初期残差R₀である。６０７は共役勾配法に従った計算処理が発散したと判断されたときの反復回数である。即ち、図６は、上記式（１）を満たしている（εCRD_t＜CRI_t）場合のものである。

図５は、ステップＳ１０５における処理のフローチャートである。先ず、残差R_t＜R_MINであるのか否かを判断し（ステップＳ５０１）、残差R_t＜R_MINの場合には処理をステップＳ５０２に進め、上記式（４），（５）に従った処理を行う（ステップＳ５０２）。

そして以下の式（１０）に従ってCRD_tを更新する処理を行う（ステップＳ５０３）。

CRD_t ← CRD_t-1 + CRD2 （１０）
ただし、CRD₀、CRD2はステップＳ１０１において０に初期化されているものとする。そして処理をステップＳ５０４に進め、CRI_t、CRD2を以下の式（１１）、（１２）に示す如く、０に初期化する（ステップＳ５０４）。

CRD2 = 0 （１１）
CRI_t = 0 （１２）
そして本フローチャートに従った処理を完了させ、処理をステップＳ１０６に戻す。

一方、ステップＳ５０１において、残差R_t≧R_MINの場合には処理をステップＳ５０５に進め、以下の式（１３）に従ってCRD2を更新する処理を行う（ステップＳ５０５）。

CRD2 ← CRD2 + (R₀ - R_t) （１３）
そして次に、以下の式（１４）に従ってCRI_tを更新する処理を行う（ステップＳ５０６）。

CRI_t ← CRI_t-1 + (R_t - R_MIN) （１４）
以上の処理が完了すると、処理をステップＳ１０６に戻す。

次に、上記共役勾配法について説明する。図７は、共役勾配法に従った計算処理（ステップＳ１０２における処理）のフローチャートである。なお、同図に示したフローチャートに従った処理は１ステップ分のものであるので実際にはこのような１ステップ分の処理を複数回繰り返して行う。

先ず、以下の式（１５）に従って、更新係数α_tの値を求める（ステップＳ７０１）。

α_t+1 = (R_t,p_t)／(p_t,Ap_t) （１５）
ここで、Ａは上記係数行列である。

次に、変数ベクトルＸ_ｔを以下の式（１６）に従って更新する（ステップＳ７０２）。

X_t+1 = X_t + α_tp_t （１６）
次に、残差ベクトルR_tを以下の式（１７）に従って更新する（ステップＳ７０３）。

R_t+1 = R_t -α_tAp_t （１７）
次に、探索方向ベクトルp_tの更新係数β_t+1を以下の式（１８）に従って求める（ステップＳ７０４）。

β_t+1 = - (R_t+1,Ap_t)／(Ap_t,Ap_t) （１８）
次に、探索方向ベクトルp_tを以下の式（１９）に従って更新する（ステップＳ７０５）。

p_t+1 = R_t+1 + β_t+1p_t （１９）
次に、上記共役残差法について説明する。図８は、共役残差法に従った計算処理（ステップＳ１０８における処理）のフローチャートである。なお、同図に示したフローチャートに従った処理は１ステップ分のものであるので実際にはこのような１ステップ分の処理を複数回繰り返して行う。

先ず、以下の処理で用いる各変数の為の記憶領域をＲＡＭ２０４内に確保すると共に、それぞれの変数の初期化処理を行う（ステップＳ８０１）。このような変数には、変数ベクトルX_k、変位拘束条件B、評価関数L₁が含まれている。

また、変数の初期化については例えば、変数ベクトルX、残差ベクトルR、探索方向ベクトルpを初期化する。それぞれの初期値をX₀、R₀、p₀と表す。ここでX₀＝X_MINとする。

また、R₀については上記X₀を用いて以下の式（２０）に従って求める。

R_k = b - AX_k （２０）
ここで、ｂは上記右辺の行列（ベクトル）である。また、p₀についてはR₀と同じ値とする。また本明細書における説明において、ｋ回目の繰り返し計算で得られたものについては、それを表す変数に添え字kを付す。

また、反復回数を示す変数ｋを０に初期化する。

図８に戻って、次に、状態ベクトルの更新係数α_kを以下の式（２１）に従って求める（ステップＳ８０２）。

α_k = (R_k-1,Ap_k-1)／(Ap_k-1,Ap_k-1) （２１）
次に、変数ベクトルX_kを以下の式（２２）に従って更新する（ステップＳ８０３）。

X_k = X_k-1 + α_kp_k-1 （２２）
次に、残差ベクトルR_kを以下の式（２３）に従って更新する（ステップＳ８０４）。

R_k = R_k-1 -α_kAp_k-1 （２４）
次に、探索方向ベクトルpの更新係数β_kを以下の式（２４）に従って求める（ステップＳ８０５）。

β_k = - (AR_k,Ap_k-1)／(Ap_k-1,Ap_k-1) （２４）
次に、探索方向ベクトルp_kを以下の式（２５）に従って更新する（ステップＳ８０６）。

p_k = R_k + β_kp_k-1 （２５）
そして次に、以上の計算の結果、計算が収束しているのか否かを判断し（ステップＳ８０７）、収束していれば本処理を終了するのであるが、収束していない場合には処理をステップＳ８０２に戻し、以降の処理を行う。ここで、収束しているのか否かの判定処理については特に限定するものではないが、例えば反復回数、残差ベクトルのノルム等によって予め決定された方法を用いるようにしても良い。

なお、上記ステップＳ１０２では共役勾配法の派生法を用いるようにしても良いし、同様に、ステップＳ１０８では共役残差法の派生法を用いるようにしても良い。

なお、共役残差法およびその様々な派生法は、例えば、速見：“特異な系に対するGCR(k)法の収束性について、”京都大学数理解析研究所講究録１２６５（２００２年）、阿部、緒方、杉原、張、三井：“特異な係数行列をもつ連立一次方程式に対するCR法の収束性、”日本応用数理学会論文誌、Vol.9, No.1, pp.1-13 (1999年)、K. Hayami: ”On the Behaviour of the Conjugate Residual Method for Singular Systems,” NII Technical Report, NII-2001-002E (2002年)、S-L. Zhang, Y. Oyanagi, M. Sugihara: “Necessary and sufficient conditions for the convergence of Orthomin(k) on singular and inconsistent systems,” Numerische Mathemarik, 87, pp.391-405 (2000年)等に開示されている。

また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体（または記憶媒体）を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記録媒体に格納されたプログラムコードを読み出し実行することによっても、達成されることは言うまでもない。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録した記録媒体は本発明を構成することになる。

また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

さらに、記録媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書込まれた後、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。

本発明の実施形態に係るコンピュータが行う、連立一次方程式の解を求めるための数値計算処理のフローチャートである。 本発明の実施形態に係る情報処理装置として機能するコンピュータの基本構成を示すブロック図である。 ステップＳ１０５における処理のフローチャートである。 共役勾配法、共役残差法それぞれにおける残差と反復回数との関係を示す図である。 ステップＳ１０５における処理のフローチャートである。 共役勾配法、共役残差法それぞれにおける残差と反復回数との関係を示す図である。 共役勾配法に従った計算処理（ステップＳ１０２における処理）のフローチャートである。 共役残差法に従った計算処理（ステップＳ１０８における処理）のフローチャートである。

Claims (7)

1. 連立一次方程式の解データを求める為の第１の計算方法に従った処理を行う第１の計算手段と、
   前記第１の計算手段が計算処理を反復する毎に得られる計算結果のデータに基づいて、当該計算処理が発散しているか否かを判断する判断手段と、
   前記判断手段が発散していると判断した場合には、前記第１の計算手段による計算処理の代わりに、前記連立一次方程式の解データを求める為の第２の計算方法に従った処理を行う第２の計算手段と
   を備えることを特徴とする情報処理装置。
2. 前記第１の計算方法は共役勾配法及びその派生法であり、前記第２の計算方法は共役残差法及びその派生法であることを特徴とする請求項１に記載の情報処理装置。
3. 前記第２の計算手段は、残差の最小値R_MINを更新し、当該最小値における状態変数の値X_MINを記憶し、残差の初期残差からの減少量の累積値CR-と，R_MINからの増加量の累積値CR+をそれぞれ計算し、CR-とCR+との値が予め設定された条件を満足したときに、前記第１の計算手段による計算処理の代わりに、前記値X_MINを初期値として前記第２の計算方法に従った処理を行うことを特徴とする請求項１又は２に記載の情報処理装置。
4. 前記累積値CR-は、前記最小値R_MINが出現する反復回数までの累積値であり、前記累積値CR+は、前記最小値R_MINが出現する反復回数からの累積値であることを特徴とする請求項３に記載の情報処理装置。
5. 連立一次方程式の解データを求める為の第１の計算方法に従った処理を行う第１の計算工程と、
   前記第１の計算工程で計算処理を反復する毎に得られる計算結果のデータに基づいて、当該計算処理が発散しているか否かを判断する判断工程と、
   前記判断工程で発散していると判断した場合には、前記第１の計算工程による計算処理の代わりに、前記連立一次方程式の解データを求める為の第２の計算方法に従った処理を行う第２の計算工程と
   を備えることを特徴とする情報処理方法。
6. コンピュータに、
   連立一次方程式の解データを求める為の第１の計算方法に従った処理を行う第１の計算工程と、
   前記第１の計算工程で計算処理を反復する毎に得られる計算結果のデータに基づいて、当該計算処理が発散しているか否かを判断する判断工程と、
   前記判断工程で発散していると判断した場合には、前記第１の計算工程による計算処理の代わりに、前記連立一次方程式の解データを求める為の第２の計算方法に従った処理を行う第２の計算工程と
   を実行させることを特徴とするプログラム。
7. 請求項６に記載のプログラムを格納することを特徴とする、コンピュータ読み取り可能な記憶媒体。

Images