修士論文 要旨

双共役勾配法の拡張に関する研究
数理情報学専攻 076218
指導教員

杉原正顯

教授

このような局所的な直交性を繰り返すことにより, ベク

はじめに

1

谷尾 真明

連立一次方程式の反復解法において, Krylov 部分空トル rk , Auk は以下のような大域的な直交性を得る.
命題 1 もし係数 αk , βk の決定において, 破綻が起きなけ間法が優れた枠組みとして知られている. その中でも,れば, rk , Auk ⊥ Kk (A∗ ,er0 ) が成り立つ.

Bi-CGSTAB(s) 法 [2], ML(s)BiCGSTAB 法 [4], IDR(s)2.2

Bi-CGSTAB(L) [1]法 [3] が, 近年注目されている. この 3 つの手法は, 既存

Bi-CG 法の拡張として, Bi-CG 法の残差ベクトルにの手法と比べて理論上の計算量が少ないという点で優L 次の加速多項式を付加して安定化を行うアルゴリズ

れており, 既存の Lanczos 原理に基づく方法では, 高々ムとして, Bi-CGSTAB(L) が知られている. 以下におい2N 回の行列ベクトル積の演算を行えば, 真の解に収束て, その説明を行う. k = mL において, まず L 次の多項する保証があったが, これら 3 つの方法は N + N/s 回の式 pi (i = 1, 2, ..., m) を定義して, 加速多項式 Qk をそれ行列ベクトル積の演算で収束する保証がある. よって,Bらの積:Qk ≡ p1 p2 ...pm と置く. また, ベクトル rB
k , uk−1s ≥ 2 と選ぶことにより, 理論上, より高速に収束するこは, Bi-CG における rk , uk−1 とする. Bi-CGSTAB(L) は,とが分かり, 実際に数多くの数値例において既存手法よ

1 反復において現在のベクトル Qk rBk , Qk uBk−1 から, 新たり高速に解が得られることが確かめられている.

Bなベクトル Qk+L rB
k+L , Qk+L uk+L−1 を生成する. 詳細に関ただ, 3 つの手法に共通する問題点として, 加速多項式の次数が 1 であることが原因で, 係数行列が歪対称に近い問題に対して収束性が悪化することが知られている.
本 発 表 で は, 高 次 の 加 速 多 項 式 か つ 高 速 性(行 列 ベ ク ト ル 積 N + N/s 回) を 兼 ね 備 え た 3
つ の 手 法 GBi-CGSTAB(s, L) 法, Bi-CGSTAB(s, L)
法, ML(s)BiCGSTAB(L) 法 を 提 案 す る こ と に より, こ の 問 題 に 解 決 を 与 え る.

Bi-CGSTAB(s, L)法, ML(s)BiCGSTAB(L) 法 は, そ れ ぞ れ Bi-CG(s),ML(s)BiCG に L 次の加速多項式を付加することにより導いた. GBi-CGSTAB(s, L) は, 高次の shadow resid-ual を持つ Bi-CG を新たに定義した上で (GBi-CG(s)),
L 次の加速多項式を付加を行うことにより導いた. さらに, 3 つの提案手法が generic な条件下で収束まで破綻しないことを示した. また, 数値実験により提案手法の有用性を確認した.

2

既存研究 : Bi-CG 法と Bi-CGSTAB(L) 法

2.1

Bi-CGBi-CG において, 残差ベクトルは, 以下の漸化式によって更新される:
rk+1 = rk − αk Auk ⊥ r̃k , r̃k+1 = r̃k − αk A∗ ũk ,uk+1 = rk+1 − βk+1 uk such that Auk+1 ⊥ r̃k ,ũk+1 = r̃k+1 − βk+1 ũk .

(1)

ただし, rk+1 = rk − αk Auk ⊥ erk は, 更新後の残差 rk+1 がerk と直交するように, 係数 αk を決めることを意味する.

しては [1] を参考にされたい.

本研究

3

予稿においては, 提案する 3 つの手法のうち GBi-CGSTAB(s, L) 法についてのみ触れる.
3.1

GBi-CG( s)ここでは, Bi-CG の shadow residual を, 高次に一般化したものを考える. shadow residual の一般化とは, N × se0 に対して, rk ⊥ Kk (A∗ , Re0 ) を満たのフルランク行列 R
す残差の更新を意味する.
Bi-CG(s) 法 [2] と同様に, shadow residual の高次化に伴って, 補助ベクトルを高次に拡張した上で, 補助ベクトルを反復ごとに入れ替えることにより, 以下のアルゴ
リズム (GBi-CG(s)) が導かれる:
GBi-CG(s) アルゴリズム
1. set x0 and calculate r0 = b − Ax0
2. set s × s matrix U0 = [r0 , Ar0 , ..., A s−1 r0 ]3. k = 0
4. while ||rk || > εek
5. rk+1 = rk − AUk αk ⊥ R
6. xk+1 = xk + Uk αke
7. Uk+1 e1 = rk+1 − Uk β(1)k+1 such that AU k+1 e1 ⊥ Rk
8. v = AUk+1 e19. for j = 2, ..., s
10. set Uk( j)= [rk+1 , AUk+1 e1 , ..., AUk+1 e j−2 , Uk e j , ..., Uk e s ]ek
11. Uk+1 e j = v − U ( j) β( j) such that AUk+1 e j ⊥ Rkk+1

12. v = AUk+1 e j
13. end For
ek = Qk (A∗ )Re0
14. k = k + 1, set R

双共役勾配法の拡張に関する研究 (谷尾15. end while

10真明)

2

Bi-CGSTAB(3)ML(3)Bi-CGSTAB
GBi-CGSTAB(3,3)

また, 上記のアルゴリズムに対して以下の性質が言える:e∗ AUk , σ( j) := Re∗ AU ( j) と定義する. また,命題 2 σk := Rkkkk

Uk0 を N × s の行列で, Uk0 := [rk , AUk e1 , ..., AUk e s−1 ] と定e∗ U 0 を新たに定義する. この時, i < k
義した時, σ0k := Rk k( j)0となる, σi , σi , σi ( j = 2, ..., s) が正則であるならば, (i)アルゴリズムは k 反復まで破綻せず, (ii)rk , AUk e j ( j =e0 ) が満たされる.
1, ..., s) ⊥ Kk (A∗ , R
3.2 GBi-CGSTAB( s, L)

10101010

100-2-4-6-80GBrGBk と行列 U k−1 はそれぞれ, GBi-CG(s) における, 残差

101010ベクトル rk と s 本の補助ベクトル Uk−1 とする.

GBi-CGSTAB(s, L) は, 1 反 復 に お い て 現 在 のGB
GBQk rGB
k , Qk U k−1 から, 次のベクトル Qk+L rk+L と, 行列1010Qk+L UGBk+L−1 を生成する. 詳細は省略する.
3.3

動作保証

今回提案する 3 つのアルゴリズムについて, 以下の命

1500

例 1)

Bi-CGSTAB(L) と全く同様に, k = mL において, 加速多項式を Qk ≡ p1 p2 ...pm と定義する. ここで, ベクトル1000

図 1. s = L = 3 と設定した, 3 つの手法の比較 (数値Bi-CGSTAB(L) の構造を利用することにより, GBiCG(s) 法から GBi-CGSTAB(s, L) が導かれる.

500

102Bi-CGSTAB(3)ML(3)Bi-CGSTAB
GBi-CGSTAB(3,3)

0-2-4-6-80100200300400500600700

図 2. s = L = 3 と設定した, 3 つの手法の比較 (数値例 2)

題が成り立つ :命題 3 GBi-CGSTAB(s, L) 法, Bi-CGSTAB(s, L) 法,効果的であることが分かり, GBi-CGSTAB(3, 3) は Bi-ML(s)BiCGSTAB(L) 法は generic な条件下において, 収CGSTAB(3) と同様の収束性を持つことが分かった. 2 つ束まで破綻しない.

目の数値例としては matrix market より raefsky2 を用いこの命題によって, 初期の shadow residual をランダムに選ぶことでアルゴリズムが確率 1 で収束まで破綻せずに動く保証が得られた.

4数値実験
ここでは, 2 つの数値例を紹介する. 解法としては既存の手法 Bi-CGSTAB(L), ML(s)BiCGSTAB, 提案手法の GBi-CGSTAB(s, L) を用いて, s = L = 3 と設定して実験を行った. 一つ目は 3 次元対流問題で出てくる以下の偏微分方程式について, 境界が [0, 1] × [0, 1] × [0, 1] のディリクレ条件の下で離散化を行った際に出てくる, 以下の連立一次方程式を扱う:

u xx + uyy + uzz + 1000u x = F.
ただし u(x, y, z) = exp(xyz) sin(πx) sin(πy) sin(πz) を満たすように, 関数 F を決定した. 適切な条件下で離散化することで, 125000 × 125000 の連立一次方程式を得た.
実験結果は図 1 のようになった. 実験結果より, この具体例においては, 加速多項式を高次に取ることが非常にた. 右ベクトルは解の要素がすべて 1 になるよう設定した. 結果は図 2 のようになった. 実験結果より, この具体例においては GBi-CGSTAB(3, 3) が ML(3)BiCGSTABが同程度の収束性であり, Bi-CGSTAB(3) がそれより少し劣ることが分かる. この 2 つの結果から, 提案手法の収束性は 2 つの既存手法の, 良い方の収束性に一致することが分かる.

参考文献[1] G. L. G. Sleijpen and D. R. Fokkema, BICGSTAB(L) forequations involving unsymmetric matrices with complex
spectrum, ETNA, 1 (1993), pp. 11–32.
[2] G. L. G. Sleijpen, P. Sonneveld and M. B. van Gijzen, BiCGSTAB as an Induced Dimension Reduction Method,Appl. Math. Anal. REPORT 08-07 (2008), Delft Univ.
Tech.
[3] P. Sonneveld and M. B. van Gijzen, IDR(s): A Family ofSimple and Fast Algorithms for Solving Large Nonsymmetric Systems of Linear Equations, SIAM J. Sci. Comp.,
31 (2008), pp. 1035–1062.
[4] M. Yeung and T. F. Chan, ML(k)BiCGSTAB: A
BiCGSTAB Variant based on multiple Lanczos starting
vectors, SIAM J. Sci. Comp., 29 (1999), pp. 1263–1290.