数理解析研究所講究録
第 1791 巻 2012 年 31-36

31

降順に展開した漸化式を使用しない GBiCGSafe 法の考案
A proposal of GBiCGSafe method
without reverse-ordered recurrence
藤野 清次 (九州大学情報基盤研究開発センター)
Seiji FUJINO (Research Instimte for Information Technology, Kyushu University)

Abstract:
We consider to solve a linear system of equations $Ax=b$ by iterative method. Producttype of iterative methods such as GPBiCG, BiCGSafe and GPBiCGAR method constimte apolynomial list by multiplying Lanczos polynomials by acceleration polynomial. In this paper,
we focus on the order of the development of recunence for GPBiCG method and propose aGPBiCGSafe method by improving GPBiCG AR method. Through numerical experiments,
we make clear that the GPBiCGSafe method has excellem convergence rate and robusmess.

1

はじめに

自然科学や工学などの様々な解析分野において，数値シミュレーションは重要な役割を果たしてきた．それらの様々な解析分野から生じる問題は，偏微分方程式で記述されることが多い．偏微分方程式は有限要素法などにより離散化することで，連立一次方程式に帰着される．この連立一次方程式の解法として反復法が提案されている．GPBiCG 法やBiCGSafe 法，GPBiCG-AR 法などの積型反復法は，通常のランチョス多項式に加速多項式を掛けて，積の形の多項式列を構成することで反復法の収束を加速させるという発想に基づく反復法である．最近，阿部らにより GPBiCG 法の変形版が提案された．そこでは，残差ベクトルの更新で用いられる加速多項式に交代漸化式の代わりに三項漸化式が適用された．以下では、 それらを AS-GPB $iCG_{-}v1,$ $AS_{-}GPBiCG_{-}v2$ 法と呼ぶ．
本論文では，GPBiCG 法において，交代漸化式から三項漸化式に変形するときに用いられる降順な漸化式の存在に着目し，それらを消すことによって GPBiCG-AR 法の収束性の改良を図ったので，その結果を検証する．
本論文の構成は以下の通りである．第 2 節で，GPBiCG 法の概要と漸化式の展開について述べ，GPBiCG-AR 法を改良した GPBiCGSafe 法を導く．第 3 節で，数値実験により，GPBiCGSafe 法の収束性を検証する．第 5 節で，まとめと今後の課題を述べる．

2 積型反復法
2.1 GPBiCG 法の算法の中の降順展開した漸化式
解くべき連立一次方程式を$Ax$とする．ただし，係数行列$A$(2.1)$b$$=$は大きさ $n\cross n$ の実非対称行列，$x$ とルと右辺ベクトルと各々する．初期近似解を$b$は次数$n$の解ベクトとし，初期残差ベクトルを $r_{0}=b-Ax_{0}$とする．このとき，BiCG(Bi-Conjugate Gradient) 法において，反復 回目の残差ベクトル$x_{0}$$n$$r_{n}^{BiCG}$と補助ベクトル$p_{n}^{BiCG}$

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

$r_{n}^{BiCG}$$=$$R_{n}(A)r_{0},$$p_{n}^{BiCG}=P_{n}(A)r_{0}$

.

(2.2)

32
多項式瑞 (A) と $P_{n}(A)$ は次の交代漸化式を満たす．
$R_{C}(\lambda)$$=$

1,$R_{n}(\lambda)$$=$$R_{n-1}(\lambda)-\alpha_{n-1}\lambda P_{n-1}(\lambda)$$=$$R_{m}(\lambda)+\beta_{n-1}P_{n-1}(\lambda),$$P_{n}(\lambda)$$P_{0}(\lambda)=1$

,

(2.3)

$BiCG$ 法の残差列を $\{r_{0}^{BiCG}, r_{1}^{BiCG}, \ldots, r_{n}^{BiCG}\}$

,

(2.4)

$n=1,2,$.

$\ldots$(2.5)とする．このとき，適当な手続きにより$\{H_{0}(A)r_{0}^{BiCG}, H_{1}(A)r_{1}^{BiCG}, \ldots, H_{n}(A)r_{n}^{BiCG}\}$ を構
多項式列 $\{H_{0}, H_{1}, \ldots, H_{n}\}$ を生成し，$BiCG$ 法の残差列の収束の加速を図る．このように，残差が $BiCG$ 法の残成することで，差と加速多項式との積で定義された解法は積型反復法と呼ばれる．すなわち，積型反復法の残差ベクトル$r_{n}$

は(2.6)$=H_{n}(A)r_{n}^{BiCG}=H_{n}(A)R_{\eta}(A)r_{0}$$r_{n}$と表される．次に，固有値 を使って多項式列$\lambda$$\{G_{n}(\lambda)\}$と$\{H_{n}(\lambda)\}$を次の (交代) 漸化式で構成する．

,$H_{0}(\lambda)$

$=$

1,

$H_{n}(\lambda)$$=$$H_{n-1}(\lambda)-\lambda G_{n-1}(\lambda)$$G_{n}(\lambda)$$=$$\zeta_{n}H_{n}(\lambda)+\eta_{n}G_{n-1}(\lambda),$$G_{0}(\lambda)=\zeta_{0}$

(2.7)

,

(2.8)

$n=1,2,$$\ldots$.

(2.9)$H_{1}arrow G_{1}arrow H_{2}arrow G_{2}arrow H_{3}arrow G_{3}arrow\cdots$と通常，式 (2.8) と式 (2.9) の漸化式
は，と式(2.9)(2.8)の順番で更新される．すなわち，式は昇順に展開される．一方，GPBiCG$G_{n}$$H_{n}$

法では，多項式列は昇順に展開するのではなく，次のように降順に展開して使用の多項式を次の $(n+1)$ ステップ目のステップ目のされた．すなわち，$\{G_{n}(\lambda)\}$$H_{n+1}(\lambda)$$G_{n}(\lambda)$$n$を使って降順に展開した漸化式を使用した．
$G_{n}(\lambda)$$=$(2.10)$-(H_{n+1}(\lambda)-H_{n}(\lambda))/\lambda$

このように降順に展開した漸化式を使うと，多項式列 $\{H_{n}(A)\}$ は次の 3 項から成る漸化式が得られるという利点がある．一方，数学的には同値でも，降順展開された漸化式は，まるめ誤差の影響を受け易くなる可能性がある．

ここで，$\zeta_{n},$$\eta_{n}$,$H_{0}(\lambda)$

$=$

1,

$H_{n+1}(\lambda)$$=$$(1+\eta_{n}-\zeta_{n}\lambda)H_{n}(\lambda)-\eta_{n}H_{n-1}(\lambda),$$H_{1}(\lambda)=(1-\zeta_{0}\lambda)H_{0}(\lambda)$(2.11)$n=1,2,$は新規のパラメータである．生成された多項式列$\ldots$$\{H_{n}(\lambda)\}$(2.12)

は，任意の$n$に対して $H_{n}(0)=1$ を満たすので，$H_{n+1}(0)-H_{n}(0)=0$ である．ここで，漸化式 (2.12) は式 (29) を式 (28) に代入し，それを式 (2.10) を使って書き直すと容易に得られる．以下でをと略記す．
は，多項式$H_{n},$ $G_{n}$$H_{n}(\lambda),$ $G_{n}(\lambda)$

ここでは，算法を構成する過程で降順な漸化式 (2.10) が使われる様子を，中間ベクト
ル $y_{n}$ と $u_{n}$を取り上げ細かく観察する．すなわち，それらのベクトルの定義式と式変形は次のように各々書き表せる．下線部分が降順の漸化式 (2.10) の部分である．
$\lambda G_{n}R_{\eta+2}$$\lambda G_{n}P_{n}$$=$$(H_{n}-H_{n+1})R_{\eta+2}$$=$$H_{n}(R_{\eta\eta+1}-\alpha_{n+1}P_{n+1})-H_{n+1}(R_{n+1}-\alpha_{n+1}P_{n+1})$$=$$\lambda P_{n}(\zeta_{n}H_{n}+\eta_{n}G_{n-1})$$=$$\zeta$$=$$\zeta_{n}\lambda H_{n}P_{n}+\eta_{n}((H_{n-1}-H_{n})R_{\eta}+\beta_{n-1}\lambda G_{n-1}P_{n-1})$n$\lambda$

HnPnn$+\eta$$(\lambda$

Gn ll ち$\underline$(2.13)$+\beta_{n-1}\lambda G_{n-1}P_{n-1})$(2.14)

33
2.2 GPBiCG-AR 法
GPBiCG 法の算法中の式 (23) において，降順な漸化式を使わないようにした算法が GPBiCGAR 法である．前小節で考察したように，降順な漸化式を使わないようにするためには，中間ベクトル$y_{n}$を使用しなければよいが，パラメター$\zeta_{n}$と$\eta_{n}$を決めることができなくなる．そこで，GPBiCG 法のように残差ベクトルの 2 ノルムの最小化からではな
く，准残差ベクトルの 2 ノルムの最小化からパラメータの値を定めた算法が GPBiCG-AR法である．この施策により，中間ベクトル yn(および ) を使う必要がなくなったが，中間ベクトルはそのまま残ることになった．
$w_{n}$$u_{n}$

2.3 GPBiCGSafe 法
前述のように，GPBiCG AR 法において，中間ベクトル の更新において降順な漸化式をベクトルを用いて次のように書が 1 度だけ使用されている．そこで，ベクトル(下線を付けた部分)により降順な漸化式(2.10)の使用が消滅した．
き表す．この式変形$u_{n}$$u_{n}$$\lambda G_{n}P_{n}$$z_{n}$$=$$\lambda P_{n}(\zeta_{n}H_{n}+\eta_{n}G_{n-1})$$=$$\zeta_{n}\lambda H_{n}P_{n}+\eta_{n}(\lambda G_{n-1}R_{n}+\beta_{n-1}\lambda G_{n-1}P_{n-1})$(2.15)

このように，降順な漸化式 (2. 10) を使わないよう改良した算法を GPBiCGSafe 法と呼ぶ．
また，下線部分に該当するベクトルの計算 $(Az_{n-1})$ は従来の GPBiCG-AR 法でも使用されていたため計算量は増えない．

3

数値実験

3.1

計算機環境と計算条件

計算はすべて倍精度浮動小数点演算で行った．計算機は Nehalem(CPU: Intel Xenon $X5570$ ,クロック周波数 :2.93GHz, メモリ : $24Gbytes$ , OS :RedHat Entelprise Linux 5.6) を用いた．プログラムは Fortran90 を用いて実装し，最適化オプションは”-O3” を使用した．右$b=A\hat{x}$辺項ベクトル は厳密解を $\hat{x}=(1,1, \ldots, 1)^{T}$ とし，で作成した．収束判定条件は相対残差の 2 ノルム : $||r_{k+1}||_{2}/||r_{0}||_{2}\leq 10^{-10}$ とした．初期近似解はすべて とした．
$b$

$0$

$x_{0}$

行列は予め対角スケーリングによって対角項をすべて 1.0 に正規化した．最大反復回数は 10000 回とした．調べた反復法は GPBiCG 法，AS-GPBiCG-vl 法，AS-GPB $iCG_{-}v2$ 法，
GPBiCGAR 法，GPBiCGSafe 法，BiCGSafe 法の計 6 種類とし，前処理としてフィルインを考慮しない ILU(0) 分解を選んだ．初期シャドウ残差$r_{0}^{\star}$には初期残差$r_{0}$を代入した．

3.2 実験結果
表 1-2 に 5 種類の ILU(0) 前処理っき積型反復法の収束性を示す．表中の $\max$ は最大反復回数までで収束しなかったことを意味し，“TRR” は真の相対残差 (True Relative Residual)の常用対数 $\log_{10}(||b-Ax_{k+1}||_{2}/||b-Ax_{0}||_{2})$ の値を意味する．最大反復回数に達するまでに収束しかつ近似解に対する真の相対残差の 2 ノルムの大きさ TRR が $10^{-9.60}$ 以上の場合は偽収束 (表中の TRR の数字に括弧をつけた) と判定した．表 3 に 5 種類の解法における収束 (未収束) ケース，偽収束ケースと最速ケースを示す．表 4 に 6 種類の反復解法の速”度順位とその総合順位を示す．ただし，成績表の集計について，計算時間が同じ場合は同34一順位にした．また収束しなかった場合は順位を 6 位にした．表の観察から以下の知見が得られる．
$\bullet$

GPBiCGAR 法と GPBiCGSafe 法，BiCGSafe 法は全ての行列で収束したが，残り 3種類の解法は反復回数が最大回数に達したり，収束しても偽収束現象が発生したときがあった．

$\bullet$$\bullet$Dc の 2 ケース，AS $-GPBiCG_{-}v2$ 法では行列 $sme3Db$ の 1 ケースで偽収束が発生した．また GPBiCG 法で行列 $sme3Db$ の 1 ケースで偽収束が発生した．

AS-GPBiCG-vl 法では行列 bcircuit,$sme3$

17 種類のテスト行列において，最も速く収束したケースが多い解法は GPBiCGSafe法の 8 ケースである．

$\bullet$

GPBiCG 法の最も速く収束したケースが零に対し，AS-GPBiCG-vl,$-2$

法は 1 ケースと 2 ケースであるが，GPBiCG 法よりも計算時間が長いケースが多い．
$\bullet$

行列 bcircuit において，GPBiCGSafe 法の計算時間は GPBiCG-AR 法の計算時間より約 10% ほど速い．

Table 1: 6 種類の LU(0) 前処理つき積型反復法の収束性

4

まとめ

本論文では，積型反復法で降順な漸化式で定義される多項式 $\{G_{n}(A)\}$ の存在に着目し，

GPBiCG-AR 法において降順な漸化式を使用しない GPBiCGSafe 法を提案した．数値実験の結果，提案法が収束性において優れていることがわかった．

35
Table 2: 6 種類の ILU(O) 前処理つき積型反復法の収束性 (cont’d)

Table 3: 6 種類の解法における収束 (未収束) ケース，偽収束と最速ケース

36
Table 4: 6 種類の反復解法の速度順位とその総合順位

References[1] S.-L. Zhang, GPBi-CG: Generalized product-type methods preconditionings based on Bi-CG for solving nonsymmetric linear systems, SIAM J. Sci. Comput., pp.537-551, 1997.

[2]

藤野清次，藤原牧，吉田正浩，準残差の最小化に基づく BiCGSafe 法の収束性について，Transactions of JSCES, Vol.2005, 20050028, 2005.

[3] Moethuthu, S. Fujino: Stability of GPBiCGAR method based on minimization of asso-ciate residual, Joumal of ASCM, pp. 108-120, 2008.

[4] K. Abe, G.L.G. Sleijpen, Solving linear equations with a stabilized GPBiCG method, $B$

本応用数理学会第 14 回環瀬戸内応用数理研究部会，岡山理科大学，Jan., 2011.