第 28 回数値流体力学シンポジウム
講演番号 E06-1

一般化残差切除法の収束特性の検討
A study on convergence characteristics of the generalized residual cutting method○ 阿部敏彦, ASI 総研, 東京都千代田区内神田 1-18-14, E-mail: toshihiko-abe@asiri.co.jp
関根義人, ASI 総研
菊地一雄, ASI 総研
Toshihiko ABE, ASIRI, 1-18-14, Uchi-Kanda, Chiyoda-ku, Tokyo
Yoshihit SEKINE, ASIRI
Kazuo KIKUCHI, ASIRIThe residual cutting method is known to have superior characteristics of convergence for elliptic partial differentialproblems. In this paper, we propose the generalized version of the residual cutting method for general sparse matrices
and investigate its performance by numerical experiments.

1. はじめに
著者らは楕円型偏微分方程式の数値解法として優れた収束性を有す残差切除法(1)を、内部ソルバの代わりに係数行列と残差ベクトルから生成したベクトルを用いることにより Coupled Perturbed(CP)方程式に適用し、優れた収束性を示すことを報告している(2)。
本報告では、この手法を一般化残差切除法(GRC 法)として、一般の疎行列を対象とした線形ソルバとしての性能を評価する。
2. 残差切除法
2.1 基本アルゴリズム
偏微分方程式を適切に離散化して得られる線形化連立方程式を
𝐇𝐔 = 𝐛(2.1)
により表す。ここで𝐇は係数行列、𝐔は解ベクトルを表す。𝐔の近似解𝐔𝑚 に対する残差 𝐫 𝑚 を
𝐫 𝑚 = 𝐛 − 𝐇 𝐔𝑚(2.2)
により定義し、式(2.1)の解𝐔を近似解𝐔𝑚 と摂動量𝛟で次のように表す.
𝐔 = 𝐔𝑚 + 𝛟(2.3)
以上の式から、摂動量の残差方程式
𝐇𝛟 = 𝐫 𝑚(2.4)
が得られる。
残差切除法により式(2.4)を解くには、まず内部ソルバにより、暫定近似解𝚿𝑚 を求める。内部ソルバとしては、緩和法や ADI 法などの反復解法を用いる。次に残差の L2 ノルムを最小とする合成摂動量𝛟𝑚 と新しい近似解𝐔𝑚+1 を次のように定義する。
𝛟𝑚 = 𝛼1 𝚿𝑚 + ∑𝐿𝑙=2 𝛼𝑙 𝛟𝑚−𝑙+1(2.5)
𝐔 𝑚+1 = 𝐔𝑚 + 𝛟𝑚(2.6)
𝛼𝑙 (𝑙 = 1,2, ⋯ , 𝐿)は次に述べる残差の最小化により定まる定数である．残差の L2 ノルムをつねに小さくするため次のように残差の最小化を行う．すなわち、新しい近似解𝐔𝑚+1 が式(2.6)で表されるとき、これに対する残差𝐫 𝑚+1 は式(2.2)から次のように与えられる。
𝐫 𝑚+1 = 𝐛 − 𝐇𝐔𝑚+1(2.7)
𝐿𝐫 𝑚+1 = 𝐫 𝑚 − 𝛼𝑙 𝐇𝚿𝑚 − ∑𝛼𝑙 𝐇𝛟𝑚−𝑙+1𝑙=2
従って、近似解𝐔𝑚+1 に対する残差𝐫 𝑚+1 の L2 ノルム
𝑟 𝑚+1 = ‖𝐫 𝑚+1 ‖を最小化するためには、
𝜕‖𝑟 𝑚+1 ‖𝟐𝜕𝛼𝑙(2.8)=0(2.9)を解けば良い。式(2.9)は𝛼𝑙 の𝐿元連立方程式となり、これを数値
Copyright © 2014 by JSFM1的に解くことによってα1 ⋯ 𝛼𝑙 が定まる。実際には、𝛼𝑙 の数は数個程度で十分である。このように最小二乗法でα1 ⋯ 𝛼𝑙 が定まれば式(2.5)から𝛟𝑚 が、また式(2.6)から新しい近似解𝐔𝑚+1 が求められる。式(2.2)に戻ってこのプロセスを繰返すことによって、残差のL2 ノルム式(2.8)を零または最小にする解𝐔 ∞ に収束する。
計算の手順は次のようになる。
(1) 初期値𝐔0 を与える。
(2) 残差𝐫 𝑚 = 𝐛 − 𝐇𝐔𝑚 を求める。
(3) 𝐇𝛟 = 𝐫 𝑚 を、内部ソルバによる最小単位の反復で暫定近似解𝚿𝑚 を求める。
(4) 式(2.9)により残差最小化係数𝛼𝑙 (𝑙 = 1,2, ⋯ , 𝐿)を計算する。
(5) 式(2.5),(2.6)により、合成摂動量𝛟𝑚 と新しい近似解
𝐔𝑚+1 を得る。
𝛟𝑚 = 𝛼1 𝚿𝑚 + ∑𝐿𝑙=2 𝛼𝑙 𝛟𝑚−𝑙+1 , 𝐔𝑚+1 = 𝐔𝑚 + 𝛟𝑚(2.10)
ただし𝑚 − 𝐿 + 1 < 0の場合は和の中で𝛟の添え字が負になるものがあるので、その場合は零ベクトルとする。すなわち
𝛟𝑛 = 0 for all 𝑛 < 0(2.11)(2)から(5)の手順を繰り返すことにより、残差を減少させていく。
2.2 一般化残差切除法 (GRC 法)
従来法では𝚿𝑚 は内部ソルバによる𝐇𝛟 = 𝐫 𝑚 の暫定近似解であり、通常 SOR などの緩和法を用いる。緩和法を用いた場合（ここでは簡潔のため、Gauss-Seidel 法の例を示す）は、繰り返しごとの演算は
𝑥 𝑖 = (𝐻𝐿 + 𝐻𝐷 )−1 (−𝐻𝑈 𝑥 𝑖−1 + 𝑟 𝑚 ) ≡ 𝑓(𝑥 𝑖−1 )のように書ける。ただし𝐻𝐿 、𝐻𝑈 、𝐻𝐷 は係数行列をそれぞれ左三角、上三角、対角成分に分割したものである。内部ソルバではこの演算を n 回適用するので、暫定近似解は
𝛹 𝑚 = 𝑓 𝑛 (𝑥0 ) = 𝐹(𝐻, 𝑟 𝑚 )となる。これに対し、提案法では一般の表現として
𝛹 𝑚 = 𝐺 (𝜙 𝑚−1 , 𝐻, 𝑟 𝑚 )とする。ただし G は線形演算子で、ここでは𝐺(𝜙 𝑚−1 ) = (1 − 𝐻)𝜙 𝑚−1 + 𝑟 𝑚を用いる。つまり式(2.4)の暫定近似解を求める代わりに、
𝚿𝑚 = (𝟏 − 𝐇)𝛟𝑚−1 + 𝐔𝑚
(2.12)とする。この演算の目的は、近似解を得ることではなく、Krylov部分空間を構成することである。つまり、(2.9)により残差最小化により近似解を得るために、個々のベクトルは必ずしも近似解である必要はない。

第 28 回数値流体力学シンポジウム
講演番号 E06-1
3. 数値実験
3.1 実験条件
提案した GRC 法を、オリジナルの残差切除法(RC 法)、BiCGSTAB、GMRES 法と比較した。BiCGSTAB および GMRES 法には LIS ライブラリ(3)を用いた。また、できるだけ同じ条件で比較するために、GRC 法および RC 法も LIS のサブルーチンとして実装した。
式(2.5)における𝛼𝑙 (𝑙 = 1,2, ⋯ , 𝐿)の個数については、GRC 法、RC 法ともに𝐿 = 5とした。また、RC 法の内部ソルバとして SOR を用い、緩和係数 1.9 および繰り返し回数は 50 回とした。また GMRES のリスタートまでの回数は、LIS でのデフォルトである 40 回とした。
実験用の係数行列として、(1) 非等間隔格子上でPoisson方程式を離散化して得られる係数行列（フルノイマン条件、格子点数 81×53×55＝約24 万点）(2) 偏微分方程式
uxx + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 + 1000𝑢𝑥 = 𝐹を離散化することにより得られる係数行列(4)（125000 次元）(3) University of Florida sparse matrix collection(5)のデータベースに登録されている係数行列の一つである、‘raefsky2’の 3 つの行列を用いた。
3.2 実験結果
図 3.1、3.2、3.3 に(1)(2)(3)のそれぞれの結果を示す。上段は繰り返しステップ数、下段は時間の関数としての収束特性である。
(1)の場合は、RC 法が最も収束が早い。これは RC 法の内部ソルバとしての緩和法が、差分化された Poisson 方程式に対して効率的に収束しているためと考えられる。それと比較して GRC 法は、前半は遅く、後半は次第に加速していくような収束曲線を示している。BiCGSTAB 法は、残差が単調減少しないため独特の振動を伴うが、収束曲線の概形としてはGRC 法に似ている。
GMRES法は途中で収束が止まっている。
(2)の場合は、RC 法は残差ノルムが全く減少していない。
これは内部ソルバの緩和法が発散しており、残差切除法のアルゴリズムでは原理的に残差ノルムは増加することはないため、一定に保たれている結果であることを確認している。GRC法、GMRES 法とも似たような収束曲線を示している。それよりも遅れて、BiCGSTAB が収束している。
(3)の場合では、RC 法は(2)と同じく内部ソルバが発散しているため、途中までは残差ノルムが変化しないが、その後、発散している。これは内部ソルバの発散の程度が著しく、残差ノルムが 10100 程度のオーダーに達しているため、残差最小化が正確に行われないためであることを確認している。
BiCGSTAB が顕著に早い収束を示した。
かなり遅れて GMRESが収束し、さらに遅れて GRC 法が収束している。

Copyright © 2014 by JSFM2

Fig. 3.1

Results from the coefficient matrix (1)Fig. 3.2

Results from the coefficient matrix (2)第 28 回数値流体力学シンポジウム
講演番号 E06-1
4. まとめ
4 つの手法の中で、まず RC 法については、実験結果(1)のように、内部ソルバとしの緩和法が有効な係数行列に対しては、最も速く収束した。また GRC 法、GMRES 法、BiCGSTAB 法については、それぞれ異なる収束結果を示し、今回の実験結果からはどれか一つが特別に優れているようなことは見られなかった。ただし GRC 法（𝐿 = 5）は GMRES 法（40 ステップ毎リスタート）に比べて保持するベクトル本数が 5/40 となり、そのためメモリ使用量が少なくて済むので、同じ計算環境でより大きな係数行列に適用可能であると考えられる。

参考文献(1) A. Tamura, K. Kikuchi and T. Takahashi, “Residual cutting method for
elliptic boundary value problems: application to Poisson's equation,” J.
Comp. Phys., 137, pp. 247-264, 1997.
(2) T. Abe, Y. Sekine and F. Sato, “Solving a coupled perturbed equationby the residual cutting method,” J. Chem. Phys., 557, pp. 176-181, 2013.
(3) http://www.ssisc.org/lis/index.ja.html(4) G. L.G. Sleijpen and R. Fokkema, “BICGSTAB(L) for linear equations
involving unsymmetric matrices with complex spectrum,” ElectronicTransactions on Numerical Analysis. vol 1, pp. 11-32, 1993.
(5) http://www.cise.ufl.edu/research/sparse/matrices/index.html

Fig. 3.3

Copyright © 2014 by JSFM3

Results from the coefficient matrix (3)