数理解析研究所講究録
第 1791 巻 2012 年 11-20

11近似解の精度を改善する IDRstab 法

A variant of the IDRstab method for improving the accuracy ofapproximate solutions
相原 研輔 *, 阿部 邦美 \dag er , 石渡 恵美子 *****\dag er 岐阜聖徳学園大学経済情報学部，
東京理科大学理学部
東京理科大学大学院，

Kensuke Aihara , Kuniyoshi Abe \dag er and Emiko Ishiwata*
Graduate School of Science, Tokyo University of Science\dag erFaculty of Economics and Information, Gifu Shotoku University
Department of Mathematical Information Science, Tokyo University of Science
*

**

1

**

はじめに

大規模非対称行列を係数に持つ線形方程式 $Ax=b$ を近似的に解く反復法として，帰納的次元縮小$A$ は(Induced Dimension Reduction, IDR) 定理に基づく IDR(s) 法 [13] が提案されている．ここで，$s=1$ のとき，Bi-Conjugate Gradient
行列，右辺項 は 次元ベクトルである．IDR(s) 法は，$s>1$ のとき，BiCGSTAB 法における初期
STABilized (BiCGSTAB) 法 [15] と数学的に等価になり，シャドウ残差を 本に拡張したものであると解釈できる [9].
$b$$n\cross n$

$n$

$s$

$IDR(s)$ 法の反復過程では，BiCGSTAB 法と同様に，1 次の安定化多項式を用いて残差ノルムの最小化を行う．従って，非対称性の強い実行列に対して，IDR $(s)$ 法は不安定な収束性を示す場合がある．この問題を解決するため，IDR(s) 法の安定化多項式を 次に拡張した IDRstab 法 [12] が提案さ$\ell=1$ の$s=1$ のとき，BiCGstab $(l)$ 法 [8] と，れている．IDRstab 法は，パラメータ
に対して，
$(s)$
$(s)$
とき，IDR 法とそれぞれ数学的に等価である．また，IDR 法に 次の安定化多項式を付加した
GBi-CGSTAB $(s, L)$ 法 [14] も提案されている．この解法は，IDRstab 法と数学的に等価であるが，実装法が異なる．本論文は IDRstab 法に限って議論する．
$s>1,$ $\ell>1$ を用いる IDRstab 法は，しばしば IDR $(s)$ 法や BiCGstab $(l)$ 法よりも優れた収束性を示す．しかし，丸め誤差の影響により，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いが$\ell$$s,$$\ell$$L$異なる場合がある．このとき，真の残差ノルムは漸化式から求まる残差ノルムよりも大きい値で停滞し，近似解の精度が劣化する．そこで，我々は残差を更新する漸化式を修正することで，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いを一致させ，近似解の精度の改善を試みる．しかし，
この漸化式の修正によって，行列ベクトル積の計算回数が増加する．一方，IDRstab 法は や が大きくなるにつれて，行列ベクトル積に比べ，ベクトルを更新する AXPY の計算回数が相対的に多くなる．よって，計算時間は AXPY の計算回数に大きく左右される．AXPY とは，スカラー と 次$\ell$

$s$

$a$

$n$

元ベクトル $x,$ $y$ に対して，$ax+y$ の形式の演算を表す．本論文では，残差を更新する漸化式を修正すると同時に，AXPY の計算回数を削減する工夫を行う．結果として，近似解の精度を改善し，かつ計算効率の良い新たな IDRstab 法の実装法を提案する．
$BiCG$ 系統の解法
さらに，従来の $BiCG$ 系統の解法と同様に，IDRstab 法の残差ノルムは振動する．
に対して，残差ノルムの振る舞いを滑らかにする擬似最小残差 (Quasi-Minimal Residual, QMR) アプ
ローチの解法が開発されている．その代表的な解法として，QMR 法 [5] や QMRCGSTAB 法 [2] などがある．これらの解法が $BiCG$ 法 [4] や BiCGSTAB 法の QMR 化によって得られるのと同様に，IDR $(s)$法の残差ノルムの振る舞いを滑らかにした QMRIDR $(s)$ 法 [3] が近年では提案されている．また，文献 [16] では，数値的に安定な基底の生成過程に基づいて，Flexible QMRIDR 法および Multi-Shift
QMRIDR 法が提案されており，いずれも残差ノルムの振る舞いが滑らかである．しかし，高次の安定化多項式を用いる IDRstab 法の QMR 化は，文献 [3, 16] では扱われていない．

12
一方，残差ノルムの振る舞いを滑らかにするための異なる手法として，スムージング [7,17,18] が提$BiCG$ 法に QMR, スムージング [18] を適用すると，QMR 法によって生成される
案されてきた．特に，残差と数学的に等価な残差を生成できることが知られている [18]. 我々は，この QMR, スムージングを本論文で提案する IDRstab 法の変形アルゴリズムに適用する．このとき，我々の提案する IDRstab法の変形アルゴリズムは，残差を更新する漸化式を修正したため，計算量を大きく増加させることなく適用できる．一般に，スムージングによって得られる近似解の精度は，スムージング前の近似解の精度に比べて，向上こそしないものの，同程度の精度が維持されることが知られている [6]. すなわ

ち，我々の提案する IDRstab 法の変形アルゴリズムを利用することで，近似解の精度が改善されるという利点を引き継いだまま，効率良く IDRstab 法の QMR 化が実現できる．数値実験において，近似解の精度が改善され，残差ノルムの振る舞いが滑らかになることを示す．

2IDRstab 法
本節では，IDRstab 法の概略を述べる．IDR 定理に基づく反復法は，次元が縮小する空間列に対して， $j=0,1,2,$ に属する残差と対応する近似解を構築する．ただし，$\mathcal{G}_{j},$$\{\mathcal{G}_{j}\}$$\mathcal{G}0\equiv \mathcal{K}_{n}(A, r_{0})$$\ldots$,

は行列 $A$ と初期残差 $r0\equiv b-Ax0$ で張る 次のであり，
Krylov 部分空間，瑞は $n\cross s$ 行列島の像空間に対する直交補空間， は非零のスカラーを表す．
に属する残差IDRstab 法では， に属する残差 をへと更新する．ただし，整数 は$\mathcal{G}_{j+1}\equiv(I-\omega j+1A)(\mathcal{G}_{j}\cap\tilde{R}_{0}^{\perp})$$n$$\mathcal{K}_{n}(A, r_{0})$$\omega j$$\mathcal{G}_{k}$$\mathcal{G}_{k+\ell}$$r_{k}$$k$$r_{k+\ell}$の倍数である．この残差の更新過程を “1 サイクル” と数える．この 1 サイクルは， に属するベク
トルを口瑞へ写す “IDR step” と，IDR step で得られた 口瑞に属するベクトルに 4 次の多項式を作用させ， へと写す “polynomial step” の 2 つのステップによって構成される．以下に各ス$\ell$$\mathcal{G}_{k}$$\mathcal{G}_{k}$$\mathcal{G}_{k}$$\mathcal{G}_{k+\ell}$

テップの計算過程を示す．

2.1

IDR step

と
, および
行列
と対応する近似解
まず，1 サイクルの始まりにおいて，残差の各列ベクトルは対応する行列砿が得られているとする．ただし，に属するものとする．こ$r_{k}\in \mathcal{G}_{k}$こで，行列$\Pi_{i}^{(j)}$$AU_{k}$$n\cross s$$x_{k}$$AU_{k}$$\mathcal{G}_{k}$を$\Pi_{i}^{(j)}\equiv I-A^{i}U_{k}^{(j-1)}\sigma_{j}^{-1}\tilde{R}_{0}^{*}A^{j-i},$と定義すると，$\tilde{R}_{t}^{*}\Pi_{j}^{(j)}=O$$\sigma_{j}\equiv\tilde{R}_{0}^{*}A^{j}U_{k}^{(j-1)},$$i=0,1,$ $\ldots,j,$ $j=1,2,$$l$$\ldots,$

かつ $\Pi_{i+1}^{(j)}A=A\Pi_{i}^{(j)}$ を満たす．

を用いて残差が 回更新される．この更新に伴う計算過程を 回の “繰り返し” と呼び， 回目 $(j=1,2, \ldots, \ell)$ の繰り返しを添え字 $(j)$ ” で表す．いま，$j-1$ 回の繰り返しにおよびベクトル$\ldots,j-1$ と行列,と対応する近似解$i=1,2,$
よって，残差$i=0,1,2,$ $\ldots,j$ が得られたとする．ただし，
砿である．こIDR step では，$\Pi_{1}^{(j)}$$\ell$$\ell$$j$$A^{i}r_{k}^{(j-1)},$$x_{k}^{(j-)}$$r_{k}^{(j-1)}$$x_{k}^{(0)}\equiv x_{k},$$A^{i}U_{k}^{(j-1)},$$r_{k}^{(0)}\equiv r_{k},$$U_{k}^{(0)}\equiv$のとき， 回目の繰り返しは，以下のように行う．
$i$

まず，ベクトル$A^{i}r_{k}^{0)},$$i=0,1,$ $\ldots,j-1$

は，$\Pi_{i+1}^{(j)}$による写像によって，$A^{i}r_{k}^{(j)}\equiv\Pi_{i+1}^{(j)}A^{i}r_{k}^{(j-1)}=A^{i}r_{k}^{0-1)}-A^{i+1}U_{k}^{(j-1)}\tilde{\alpha}^{(j)},\vec{\alpha}^{0)}\equiv\sigma_{j}^{-1}(\tilde{R}_{0}^{*}A^{j-1}r_{k}^{(j-1)})$と更新される．また，対応する近似解$x_{k}^{(j)}$の漸化式は，$x_{k}^{(j)}=x_{k}^{(j-1)}+U_{k}^{(j-1)}\vec{\alpha}^{(j)}$

となる．

(1)

(2)

13
の各列ベクトルは， による写像と行列ベクトル積を用いて，Krylov 部分空に
間の基底となるように構築される．具体的には，まずベクトル$i=0,1,$を掛けることで，が得られる．そして，の 1 列目は， を用いて，次に，行列$A^{i}U_{k}^{(j)}$$\Pi_{i}^{(j)}$$A^{j-1}r_{k}^{(j)}$$\mathcal{K}_{s}(\Pi_{i}^{(j)}A, \Pi_{i}^{(j)}A^{i}r_{k}^{(j)})$$A^{j}r_{k}^{(j)}$$A^{i}U_{k}^{(j)},$$\ldots,j$$A^{i}U_{k}^{(j)}e_{1}\equiv\Pi_{i}^{(j)}A^{i}r_{k}^{(j)}=A^{i}r_{k}^{(j)}-A^{i}U_{k}^{(j-1)_{\beta_{1}^{(j)}}^{\neg}},$$A$$\Pi_{i}^{(j)}$$\beta_{1}^{(j)}\neg\equiv\sigma_{j}^{-1}\rho_{1}^{\prec j)},$(3)$\rho_{1}^{\prec j)}\equiv\tilde{R}_{0}^{*}A^{j}r_{k}^{(j)}$にと計算される．同様にして，$q<s$ に対して，ベクトルを掛けることで，のを用いて，$q+1$ 列目は，
$i=0,1,$が計算され，$c_{q}^{(j)}\equiv A(A^{j}U_{k}^{(j)}e_{q})$$A^{i}U_{k}^{(j)},$$A^{j+1}U_{k}^{(j)}e_{q}=$$A$$A^{j}U_{k}^{(j)}e_{q}$$\Pi_{i}^{(j)}$$\ldots,j$(4)$A^{i}U_{k}^{(j)}e_{q+1}=\Pi_{i}^{(j)}A^{i+1}U_{k}^{(j)}e_{q}=A^{i+1}U_{k}^{(j)}e_{q}-A^{i}U_{k}^{(j-1)}\vec{\beta}_{q+1}^{(j)},\vec{\beta}_{q+1}^{(j)}\equiv\sigma_{j}^{-1}\rho_{q+1}^{4j)},$ $\rho_{q+1}^{\prec j)}\equiv\tilde{R}_{0}^{*}c_{q}^{(j)}$$i=0,1,$ $\ldots,j-1$ かつとなる．以上の計算過程によって， 回目の繰り返しでは，$A^{i}r_{k}^{(j)}\in \mathcal{G}_{k}\cap\tilde{R}_{C}^{\perp},$$j$$A^{i}U_{k}^{(j)}e_{q}\in \mathcal{G}_{k}\cap\tilde{R}_{0}^{\perp},$2.2$q=1,2,$$\ldots,$$i=1,2,$ $\ldots,j$ を満たす．

$s,$polynomial step

IDR step の$p$

回の繰り返しによって，残差および乏 $+2$ 個の $n\cross s$ 行列$AU_{k}^{(\text{の}}$$i=0,1,$$A^{i}U_{k}^{(\ell)},$

ただし，係数れ$\ldots,$$\ell+1$は，以下のように更新される．

$r_{k+\ell}=r_{k}^{(\ell)}-\gamma_{1,k}Ar_{k}^{(l)}-\cdots-\gamma_{\ell,k}A^{\ell}r_{k}^{(\ell)},$かつと対応する近似解$r_{k}^{(\ell)}$$\gamma_{1,k},$ $\gamma_{2,k},$$AU_{k+\ell}e_{i}\in \mathcal{G}_{k+\ell},$$r_{k}^{(\ell)},$$AU_{k}^{(\ell)}$$\ldots,$$\gamma_{\ell,k}$$i=1,2,$は，ノルム$\ldots,$$s$$x_{k}^{(\ell)}$, ベクトル$A^{i}r_{k}^{(\ell)},$$i=1,2,$$l$$\ldots,$が得られる．polynomial step では，$r_{k}^{(l)}$$AU_{k+\ell}=AU_{k}^{(\ell)}-\gamma_{1,k}A^{2}U_{k}^{(\ell)}-\cdots-\gamma_{\ell,k}A^{\ell+1}U_{k}^{(\ell)}$$\Vert r_{k+\ell}\Vert_{2}$$x_{k+\ell},$.

の最小化から決定される．このとき，$r_{k+\ell}\in \mathcal{G}_{k+l}$および行列に対しても，それぞが得られる．以上，IDRstab 法の 1 サイクを満たす [12]. 近似解に対応する更新を行うことで，と$U_{k+\ell}$$U_{k}^{(\ell)}$$x_{k}^{(\ell)}$

ルを終了する．

3

IDRstab 法の変形アルゴリズム
従来の IDRstab 法では，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いが異なり，近似解の精度が劣化する場合がある．そこで，まず漸化式から求まる残差ノルムと真の残差ノルムの振る舞いを一致させるため，残差を更新する漸化式を修正する．さらに，AXPY の計算回数を削減する工夫を行うことで，効率良く近似解の精度を改善する IDRstab 法の変形アルゴリズムを導出する．

31

残差を更新する漸化式の見直し

漸化式から求まる残差ノルムと真の残差ノルムの振る舞いの違いは，近似解と残差を更新する漸化$BiCG$ 系統の解法
式に深い関わりがある [10, 11]. 我々は，近似解と残差を更新する漸化式に対して，に近い以下の形式を考える．
$x_{m+1}=x_{m}+u_{m},$ $r_{m+1}=r_{m}-Au_{m},$ここで，$u_{m}$$m=0,1,2,$は適当な探索方向ベクトルであり，ベクトル $Au_{m}$ はる．このように$A$$u_{m}$を陽に掛ける操作は，漸化式から求まる残差ノルムの振る舞いを一致させるために重要である．
一方，IDRstab 法の残差を更新する漸化式 (1) では，行列 $AU_{k}^{(j-1)}$ は$\ldots$.

(5)に $A$ を陽に掛けることで得$\Vert r_{m+1}\Vert_{2}$と真の残差ノルム$\Vert b-Ax_{m+1}\Vert_{2}$を陽に掛けるのではなく，$\Pi_{1}^{(j-1)}$$U_{k}^{(j-1)}$に対して (5) のように $A$

による写像で構築される．このため，計算されたベクトル$AU_{k}^{(j-1)}\vec{\alpha}^{(j)}$14$U_{k}^{(0)}$$x_{k}^{(0)}$$U_{k}^{(1)}e_{1}$$x_{k}^{(1)}$$arrow$...

$U_{k}^{(1)}e_{2}$$U_{k}^{(1)}e_{s}$$\Pi_{0,\nearrow}^{(1)}$$\Pi_{0,\nearrow}^{(1)}$$\downarrow A$$\downarrow A$$\downarrow A$

4$r_{k}^{(0)}$図 1: 変形した IDR step の 1 回目の繰り返し．

$U_{k}^{(1)}$$x_{k}^{(1)}$$U_{k}^{(2)}e_{1}$$x_{k}^{(2)}$$arrow$

J.

$\Pi_{0,\nearrow}^{(2)}$$AU_{k}^{(1)}$. ..

$U_{k}^{(2)}e_{2}$$U_{k}^{(2)}e_{s}$$\Pi_{0,\nearrow}^{(2)}$$2^{-1}$$J.\ulcorner 1$$\downarrow_{A}\cdot 1$$\downarrow A$$\downarrow A$$\downarrow A$$A^{2}U_{k}^{(2)}e_{1}$$A^{2}U_{k}^{(2)}e_{2}$$arrow A$$r_{k}^{(1)}$$\Pi_{1,\nearrow}^{(2)}$$\Pi_{1,\nearrow}^{(2)}$$\downarrow A$. ..

$A^{2}U_{k}^{(2)}e_{s}$$\downarrow A$$A^{2}r_{k}^{(2)}$図 2: 変形した IDR step の 2 回目の繰り返し．

は，丸め誤差の影響で正しい値からかけ離れる可能性があり，漸化式から求まる残差ノルムと真の残差ノルムの間に差が生じる場合がある．そこで，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いを一致させるため，(5) の形式に基づく以下の漸化式を導入する．
$r_{k}^{(j)}=r_{k}^{(j-1)}-A(U_{k}^{(j-1)}\vec{\alpha}^{(j)}),$(6) では，探索方向ベクトル$U_{k}^{(j-1)}\vec{\alpha}^{(j)}$に $A$のアプローチを用いる．すなわち，残差$j=1,2,$$\ell$$\ldots,$.

(6)を陽に掛ける．また，polynomial step においても，同様$r_{k+l}$は以下のように更新する．

.

(7)$r_{k+p}=r_{k}^{(p)}-A(\gamma_{1,k}r_{k}^{(p)}-\gamma_{2,k}Ar_{k}^{(l)}-\cdots-\gamma_{\ell,k}A^{p-1}r_{k}^{(\ell)})$

ここで，ベクトル$A( \sum_{i=1}^{p}\gamma_{i,k}A^{i-1}r_{k}^{(p)})$は$\sum_{i=1}^{\ell}\gamma_{i,k}A^{i-1}r_{k}^{(l)}$に $A$

を陽に掛けることで得る．以上の(6) と (7) の導入によって，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いが一致し，結果として近似解の精度の劣化を回避できると考えられる．

32

計算量の削減漸化式 (6), (7) を用いる場合，従来の方法に比べて行列ベクトル積が余分に必要となるため，計算量が増大する．一方，IDRstab 法は や が大きくなるにつれて，行列ベクトル積に比べ，AXPY の計算回数が相対的に多くなる．そこで，全体の計算量を軽減するために，AXPY の計算回数を削減す$s$$\ell$

る工夫を行う．

まず，行列$\rho_{1}^{4)}$および$A^{*}\tilde{R}_{0}$を事前に計算し，保持する．そして，文献 [1] の記述に従い，IDR step の$\rho_{q+1}^{\prec j)}(q<s)$を，それぞれ従来とは異なる式$(A^{*}\tilde{R}_{0})^{*}A^{j-1}U_{k}^{(j-1)},$ $\sigma_{j}^{-1}((A^{*}\tilde{R}_{4})^{*}A^{j-2}r_{k}^{(j-1)})$,,

によって計算する．これらの計算方法の変更によって，IDR stepを用いずに実行できる．
の繰り返しは，行列の各列ベクトルを
IDR step の 1 回目の繰り返しでは，行列による写像で構築し， にを掛けることで，行列を得る．このとき，行列が無ければ，本来は (1) の更新$(A^{*}\tilde{R}_{0})^{*}A^{j-1}r_{k}^{(j)},$$(A^{*}\tilde{R}_{0})^{*}A^{j}U_{k}^{(j)}e_{q}$の $i$ 回目 $(j=1,2, \ldots, \ell)$$A^{j}U_{k}^{0-1)}$$U_{k}^{(1)}$$A$$\sigma j,\vec{\alpha}^{(j)}$$AU_{k}^{(1)}$$\Pi_{0}^{(1)}$$AU_{k}^{(0)}\equiv AU_{k}$$U_{k}^{(1)}$15
アルゴリズム 1. Quasi-Minimal Residual スムージング [18]1. Set $s_{0}=r_{0},$ $y_{0}=x_{0},\hat{u}_{0}=\hat{v}_{0}=0,$
2. While$\tau_{0}=\Vert r_{0}\Vert_{2}$and $m=1$ .

$\Vert s_{m}\Vert_{2}>tol$

3. Compute $p_{m}=x_{m}-x_{m-1}$ and $Ap_{m}$4. $\hat{v}_{m}=\hat{v}_{m-1}+p_{m},\hat{u}_{m}=\hat{u}_{m-1}+Ap_{m}$5.

$\rho_{m}=\Vert s_{m-1}-\hat{u}_{m}\Vert_{2}$6.

$y_{m}=y_{m-1}+(\tau_{m}^{2}/\rho_{m}^{2})\hat{v}_{m},$7.

$\hat{v}_{m}=(1-(\tau_{m}^{2}/\rho_{m}^{2}))\hat{v}_{m},\hat{u}_{m}=(1-(\tau_{m}^{2}/\rho_{m}^{2}))\hat{u}_{m}$and define$\tau_{m}$by $1/\tau_{m}^{2}=1/\tau_{m-1}^{2}+1/\rho_{m}^{2}$$s_{m}=s_{m-1}-(\tau_{m}^{2}/\rho_{m}^{2})\hat{u}_{m}$8. $m=m+1$
9. End while

が行えず，残差 は得られない．しかし，我々は (6) によって残差を更新するため， は不要でに対して，(1) の更新を行い，ベクトルある．2 回目以降の繰り返しでは，$i=1,2,$$r_{k}^{(1)}$$AU_{k}^{(0)}$$\ldots,j-2$$A^{j-2}r_{k}^{(j)}$を掛けることで，を得る．さらに，$i=0,1,$ $\ldots,j-1$ に対して，(3) および (4) の更新をに行う．そして，回目の繰り返しの最後に，を掛け， を得る．各繰り返しにおけるに$A$$A^{j-1}r_{k}^{(j)}$$A^{\ell-1}r_{k}^{(\ell)}$$\ell$$A$$A^{\ell}r_{k}^{(p)}$近似解と残差は，常に (2) と (6) によって更新する．以上の計算過程により，IDR step の 回目の繰$i$

り返しでは，行列$AU_{k+\ell}$$A^{j+1}U_{k}^{(j)}$を構築する必要がなくなり，また polynomial step では，行列$AU_{k}^{(\ell)}$からへの更新が不要となる．この結果，従来の方法に比べ，1 サイクルあたり $\ell(s^{2}+3s)+\frac{1}{2}(\ell-1)-s$

回の AXPY を削減できる．計算量の詳細は，42 節で述べる．

$\ell=2$ とし，表記は文献 [12] に倣う．こ図 1-2 に，変形した IDR step の計算過程を示す．ただし，

こで，$\blacksquare$

で囲まれたベクトルは， に関する行列ベクトル積を陽に行うことで得られる．特に，従$A$

来の IDRstab 法では，

$r_{k}^{(j)}$は$\Pi_{1}^{(j)}$による写像で更新されるのに対し，図 1-2 では，(6) のように行列

ベクトル積を陽に行う．以後，本節で述べた我々が提案する IDRstab 法の変形アルゴリズムを，“変形 1” と呼ぶ．IDRstab 法の変形 1 は，従来の IDRstab 法と数学的に等価である．

4

QMR スムージングによる IDRstab 法の QMR 化
従来の $BiCG$ 系統の解法と同様に，IDRstab 法の残差ノルムは振動する．一方，我々が提案した

IDRstab 法の変形 1 は，計算量を大きく増加させることなくスムージングを適用できる．そこで本節
では，QMR スムージングを IDRstab 法の変形 1 に適用することで，効率良く IDRstab 法の QMR 化を実現し，残差ノルムの振る舞いを滑らかにする．

41

QMR スムージングの適用

文献 [18] で提案された QMR スムージングを，アルゴリズム 1 に示す．いま，$\{r_{m}\}$と $\{x_{m}\}$ を適に対する探当な反復法によって生成される残差列と対応する近似解列とする．このとき，近似解
索方向ベクトルとベクトルに対して，アルゴリズム 1 の 4-6 行の計算を行うことで，スムー$x_{m}$$p_{m}$$Ap_{m}$

ジングされた残差と対応する近似解が得られる．このようにして得られた残差ノルムの振る舞いは，QMR アプローチから得られる残差ノルムのように滑らかである [18].
アルゴリズム 1 におけるベクトル $Ap_{m}$ は，数値的安定性の理由から，探索方向ベクトルに $A$
が (5) の形式で更新される場合は，行列ベクトル積をと残差を陽に掛けることで得る．近似解
陽に行って得られたベクトル $Ap_{m}$ を再利用できるため，行列ベクトル積を増加させることなく QMR$s_{m}$$\Vert s_{m}\Vert_{2}$$y_{m}$$p_{m}$$x_{m}$$r_{m}$

16
表 1:1 サイクルあたりの行列ベクトル積数，ベクトル更新回数，内積数．

スムージングを適用できる．ここで，従来の IDRstab 法は，残差を (5) とは異なる形式で更新する．
を得るために行列ベクトル積を余分に行そのため，アルゴリズム 1 を実装するには，ベクトルう必要があり，計算効率が悪い．そこで本論文では，QMR スムージングを我々の提案した IDRstab法の変形 1 に適用する．IDRstab 法の変形 1 では，IDR step と polynomial step において，近似解と$Ap_{m}$

残差を (5) の形式で更新する．従って，IDR step と polynomial step のそれぞれにアルゴリズム 1 の計算過程を組み込むことで，行列ベクトル積を増加させることなく QMR 化が実現できる．また，
般にスムージングによって得られる近似解の精度は，スムージング前の近似解の精度と同程度となる
ことが知られている [6]. このため，QMR 化した IDRstab 法によって得られる近似解の精度は，変形
1 によって得られる近似解の精度が維持されると考えられる．
以後，IDRstab 法の変形 1 に QMR スムージングを適用したアルゴリズムを “変形 2” と呼び，アルゴ$q\leq s$ に対して，行列リズム 2 に示す．ただし，アルゴリズムの表記は MATLAB コードに倣う．特に，$W=[w_{1}, \ldots, w_{s}]$ の部分行列 $[w_{1}, \ldots, w_{q}]$ および列ベクトルと記述す
は，それぞれ $W_{(:,1:q)},$$W_{(:,q)}$$w_{q}$

る．また，$[W0;\cdots$のとき，$U_{i},$$V_{i},$;$W_{j]}\equiv[W_{0}^{T}, \ldots, W_{j}^{T}]^{T}$$r_{i}$

および$u_{i}$である．さらに，アルゴリズム 2 において，$i=0,1,$$\ldots,j$$U=[U0;\cdots;U_{j}],$ $V=$
は，それぞれ U, V, r, u の成分を表し，$r=[r0;\cdots;r_{j}],$ $u=[u_{0};\ldots;uj]$ である．
アルゴリズム 2 において，下線で記した坦行，盤行は，3 節で導入した漸化式 (6), (7) をそれぞれ$\overline{T}$ , および mT が QMR スムージングの計算，表している．また，数字で記した$[V_{0;}\ldots;V_{j}],$$\vec{\tau},$

すなわちアルゴリズム 1 の計算過程に対応する．本論文では，IDRstab 法の変形 1 のアルゴリズムを具体的に明記しないが，これはアルゴリズム 2 から QMR スムージングの計算を削除し，近似解の漸化式を導入することで得られる．なお，文献 [12] の記述に従い，アルゴリズム 2 の 6-7 行および 25-26の各列ベクトルの正規直交化を行っている．
行において，アーノルディ過程による行列$U_{0},$

42$A^{j}U_{k}^{(j)}$

計算量の比較

本節では，従来の IDRstab 法と提案した変形 1, 変形 2 の計算量を比較する．表 1 に，1 サイクルあ次元ベクトル同士の内積数 (DOTs)
たりの行列ベクトル積数 (MVs), ベクトル更新回数 (AXPYs),$n$を示す．ただし，IDR step におけるアーノルディ過程の計算量は含めない．
表 1 より，次のことが言える．従来の IDRstab 法に比べ，変形 1 は 1 サイクルあたり $\ell+1$ 回の行列ベクトル積が増加する代わりに，AXPY は $\ell(s^{2}+3s)+\frac{1}{2}(P-1)-s$ 回少ない．従って，行列の疎性が高い場合，AXPY の削減によって，計算量の増加を十分に抑えることができる．また，変形 2 は

変形 1 を利用しているため，行列ベクトル積の増加はなく，1 サイクルあたり$\ell+1$

回の DOT のみで QMR スムージングを適用している．

5

数値実験

$\frac{7}{2}(\ell+1)$

回の AXPY と本節では，従来の IDRstab 法と提案した変形 1, 変形 2 の収束性を，数値実験によって比較する．
特に，変形 1 によって，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いが一致し，近似解の精度が改善されること，変形 2 によって，残差ノルムが滑らかに収束することをそれぞれ示す．

17
アルゴリズム 2. IDRstab 法の変形 2

1. Select an initial guess and an $(n\cross s)$ matrix
2. Compute $r_{0}=b-$ Ax, $r=[r_{0}]$$x$$\tilde{R}_{\mathbb{C}}$.

% Generate an initial $(n\cross s)$ matrix $U=[U_{0}]$4. For $q=1,$
5. if $q=1,$ $u_{0}=r_{0}$ , else, $u_{0}=Au_{0}$6. $\vec{\mu}=(U_{0(:,1:q-1)})^{*}u_{0},$ $u_{0}=u_{0}-U_{0(:,1:q-1)}\vec{\mu}$7.
$U_{0(:,q)}=u_{0}$8. End for
9. While$\ldots,$$s$$u_{0}=u_{0}/\Vert u_{0}\Vert_{2},$$\Vert s\Vert_{2}>tol$

10.
%
11.
12.
13.
14.
15.
16.

$\underline{\prod 20}$

For $j=1,$
The $IDR$ step$\ell$$\ldots,$$\sigma=(A^{*}\tilde{R}))^{*}U_{j-1}$

if$j=1,\vec{\alpha}=\sigma^{-1}(\tilde{R}_{0}^{*}r_{0})$, else, $\vec{\alpha}=\sigma^{-}((A’ R)^{*}rj-2)$$r_{0}=r_{0}-A(U_{0}\tilde{\alpha})$For $i=1,$ $\ldots,j-2$$r_{i}=r_{i}-U_{i+1}\alpha^{\neg}$

End for$\hat{v}=(1-\eta)\hat{v},\hat{u}=(1-\eta)\hat{u}$21. if $j>1,$ $r=[r;ArJ-2]$
22. For $q=1,$
23. if $q=$ l, u $=r$ , else, $u=[u_{1};\ldots;uj]$24.
$\ldots,$$s$$\vec{\beta}=\sigma^{-1}((A^{*}\tilde{R}_{0})^{*}u_{j-1}),$

.

$25$

$\vec{\mu}=(V_{j(:,1:q-1)})^{*}uJ,$

26.
27.

End for

.

$U=V$

$2S$

$u=u/\Vert uj\Vert_{2},$$u=u-U\vec{\beta^{\ovalbox{\tt\small REJECT}}},$$u=[u;Au_{i-1}]$$u=u-V_{(:,1:q-1)}\vec{\mu}$$V_{(:,q)}=u$

29. End for
30. $r=[r;Ar_{l-1}]$%The polynomial step31.
32. $r_{0}=r_{0}-A([r_{0}, \ldots, r\ell-1]$ う$\vec{\gamma}=[\gamma_{1};\ldots ; \gamma_{\ell}]=\arg\min_{\vec{\gamma}}\Vert r_{0}-[r_{1}, \ldots, r_{\ell}]\vec{\gamma}\Vert_{2}$$)$$\hat{v}=\hat{v}+[r_{0}, \ldots, r_{\ell-1}]\gamma\neg,\hat{u}=\hat{u}+A([r_{0}, \ldots, r_{\ell-1}]\vec{\gamma})$

図$\rho=\Vert s-\hat{u}\Vert_{2}^{2},$$\underline{\cap 35}y=y+\eta\hat{v},$$\tau=1/(1/\tau+1/\rho),$$s=s-\eta\hat{u}$$\hat{v}=(1-\eta)\hat{v},\hat{u}=(1-\eta)\hat{u}$37.

$U=[U_{0}-\sum_{j=1}^{\ell}\gamma_{j}U_{j}]$38. End while$\eta=\tau/\rho$

18

サイクル数

図 3: $(s, \ell)=(4,2)$ に対する収束履歴．

サイクル数

図 4: $(s, l)=(8,8)$ に対する収束履歴．

数値実験は，HP の PC(Intel Core $i72.67GHz$ CPU) で，Intel $C++11$ . 1.048 コンパイラの倍精度演算で行った．テスト行列は，文献 [12] と同様に，Matrix-Market に収納されている非対称疎行列より，SHERMAN5 を取り上げた．行列の次元数および非零要素数は，それぞれ 3312,20793 である．
右辺項は，真の解あ $=(1,1, \ldots, 1)^{T}$ に対して，と設定した．初期値は , 収束判定条件は$b=A\overline{x}$$tol=10^{-12}\Vert b\Vert_{2}$とした．行列島は，各成分を区間$(0,1)$$(s, l)$ の組み合わせは (4, 2), (4, 4),

交化して与えた．

5.1$0$

上の乱数で生成し，各列ベクトルを正規直(8, 2), (4, 8), (8, 8) とした．

計算結果

およ
$(s, \ell)=(4,2),$ $(8,8)$ に対する各解法の相対残差 2 ノルム
図 3-4 に，び真の相対残差 2 ノルムの振る舞いを示す．グラフの横軸は$(\Vert r_{k}\Vert_{2}/\Vert b\Vert_{2}, \Vert s_{k}\Vert_{2}/\Vert b\Vert_{2})$$(\Vert b-Ax_{k}\Vert_{2}/\Vert b\Vert_{2}, \Vert b-Ay_{k}\Vert_{2}/\Vert b\Vert_{2})$19
表 2: 各解法のサイクル数，行列ベクトル積数，計算時間，真の相対残差 2 ノルム．

IDRstab 法におけるサイクル数，縦軸は相対残差 2 ノルムを表す．また表 2 に，各解法の収束までに要したサイクル数 (Cycles), 行列ベクトル積数 (MVs), 計算時間 (Time $[\sec]$ ), および収束時の真の相対残差 2 ノルム (True res) を示す．

および表 2 より，次のことが言える．まず従来の IDRstab 法では，漸化式から求まる残差
ノルムと真の残差ノルムの振る舞いが異なり，近似解の精度が劣化している．これに対し変形 1 では，漸化式から求まる残差ノルムと真の残差ノルムの振る舞いが一致し，近似解の精度が改善されている．
$(s, \ell)=(4,2),$ $(4,4),$ $(8,2)$ の場合，変形 1 の収束までに要したサイクル数は，従来の IDRstab 法と同程度である．このとき，変形 1 の行列ベクトル積数は，従来の IDRstab 法に比べて増加しているが，$(s, \ell)=(4$ ,計算時間は短い．これは，変形 1 で AXPY が十分に削減されたためと考えられる．一方，8 , (8, 8) の場合，変形 1 の収束までに要したサイクル数は，従来の IDRstab 法よりも多く，収束速度
図 3-4,$)$

が大きい場合に起きる傾向がある．
が低下している．このような収束速度の低下は，パラメータ
また，従来の IDRstab 法および変形 1 では，残差ノルムが振動している．これに対し変形 2 では，残差ノルムが滑らかに収束している．このとき，得られた近似解の精度は，変形 1 の場合と同程度であり，変形 1 に対する変形 2 の計算時間の増分は，平均して約 16%であった．従って，変形 2 は近似解の精度が改善されるという変形 1 の性質を引き継いだまま，計算量を大きく増加させることなく，$s,$$\ell$

残差ノルムの振る舞いが滑らかになったと言える．

6

まとめ

IDRstab 法に対して，近似解の精度を改善するため，漸化式から求まる残差ノルムと真の残差ノル
ムの振る舞いが一致する変形アルゴリズム (変形 1) を提案した．この変形 1 は，従来の IDRstab 法に比べ，行列ベクトル積を余分に必要とするものの，AXPY を大幅に削減することで，計算量の増加を抑えることができる．さらに，変形 1 を利用することで，行列ベクトル積を増加させずに，効率良く
IDRstab 法の QMR 化 (変形 2) を実現した．変形 2 の計算時間の増加はわずか 16%程度であり，得られる近似解の精度は，変形 1 によって得られる近似解の精度と同程度である．数値実験では，提案手

20法によって近似解の精度が改善され，また残差ノルムの振る舞いが滑らかになることを示した．ただし，パラメータのが大きいとき，収束速度が低下する場合があった．この原因の究明と改善策に$s,$$\ell$

ついては，今後の課題としたい．

本研究について，貴重なご意見を頂いた Utrecht 大学の Gerard L. G. Sleijpen 先生，HamburgHarburg 工科大学の Jens-Peter M. Zemke 先生に深く御礼申し上げます．

謝辞

参考文献
[1] K. Aihara, K. Abe and E. Ishiwata, An alternative implementation of the IDRstab method saving vector
updates. JSIAM Letters. 3 (2011), 69-72.

[2] T. F. Chan, E. Gallopoulos, V. Simoncini, T. Szeto and C.H. Tong, A quasi-minimal residual variantof the Bi-CGSTAB algorithm for nonsymmetric systems, SIAM J. Sci. Comput., 15 (1994), 338-347.

[3] L. Du. T. Sogabe and S.-L. Zhang. A variant of the IDR $(s)$ method with the quasi-minimal residual
strategy, J. Comput. A $ppl$ . Math., 236 (2011), 621-630.
[4] R. Fletcher, Conjugate gradient methods for indefinite systems. Lecture Notes in Mathematics, 506(1976), 73-89.
[5] R. W. Fkeund and N. M. Nachtigal. QMR: a quasi minimal residual method for non-Hermitian linear
systems, Numer. Math., 60 (1991), 315-339.
[6] M. H. Gutknecht and M. Rozlo\v{z}nik, Residual smoothing techniques: do they improve the limiting
accuracy of iterative solvers?, BIT. 41 (2001), 86-114.
[7] W. Sch\"onauer, Scientific Computing on Vector Computers, Elsevier, Amsterdam, 1987.

[8] G. L. G. Sleijpen and D. R. Fokkema, BiCGstab for linear equations involving unsymmetric matrices. Numer. Anal., 1 (1993), 11-32.
with complex spectrum, Electron.
$(\ell)$$\mathcal{T}Vans$[9] G. L. G. Sleijpen, P. Sonneveld and M. B. van Gijzen, Bi-CGSTAB as an induced dimension reduction
method, Appl. Numer. Math., 60 (2010). 1100-1114.
[10] G. L. G. Sleijpen and H. A. van der Vorst, Reliable updated residuals in hybrid Bi-CG methods, Comput.,
2 (1996), 141-163.

[11] G. L. G. Sleijpen, H. A. van der Vorst and D. R. Fokkema, BiCGstabmethods, Numer. A lgori thms, 7 (1994), 75-109.
[12] G. L. G. Sleijpen and M. B. van Gijzen, Exploiting BiCGstab
SIAM J. Sci. Comput., 32 (2010), 2687-2709.

$(\ell)$$(\ell)$and other hybrid Bi-CG

strategies to induce dimension reduction,[13] P. Sonneveld and M. B. van Gijzen, IDR $(s)$ : a family of simple and fast algorithms for solving large
nonsymmetric linear systems, SIAM J. Sci. Comput., 31 (2008), 1035-1062.

[14] M. Tanio and M. Sugihara, GBi-CGSTAB $(s, L)$ : IDR $(s)$ with higher-order stabilization polynomials,
Comput. Appl. Math., 235 (2010), 765-784.

$J$.

[15] H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of
nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 13 (1992), 631-644.
[16] M. B. van Gijzen, G. L. G. Sleijpen and J.-P. M. Zemke, Flexible and multi-shift induced dimension
reduction algorithms for solving large sparse linear system, Report 11-06, DIAM, TU Delft & Bericht
156, INS, TU Hamburg-Harburg (2011).
[17] R. Weiss, Convergence behavior of generalized conjugate gradient methods, Ph.D. thesis, Univ. of Karlsruhe, Germany, 1990.
[18] L. Zhou and H. F. Walker, Residual smoothing techniques for iterative methods, SIAM J. Sci. Comput.,
15 (1994), 297-312.