一般化Bi-CGSTAB(s, L) (=一般化IDR(s, L)) 東京大学大学院　情報理工学系研究科 数理情報学専攻 谷尾　真明 杉原　正顯

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

Krylov部分空間法 Bi-CG GMRES CGS FOM Bi-CGSTAB GCR Bi-CGSTAB(L) Orthomin 概要 Krylov部分空間法 Bi-CG GMRES CGS FOM Bi-CGSTAB Bi-CGSTAB(L) GCR Orthomin GP-BiCG Lanczos原理 (短い漸化式) Arnoldi原理 (長い漸化式) IDR(induced dimension reduction)(s)[1](2007) [1]P. Sonneveld and M. B. van Gijzen, IDR(s): A Family of Simple and Fast Algorithms for Solving Large Nonsymmetric Systems of Linear Equations, Reports Depart. Appl. Math. Anal., REPORT 07-07, Delft Univ., Tech., 2007. (短い漸化式)

Lanczos原理に基づくKrylov部分空間法 vs IDR(s)法 概要 一番重たい計算＝行列とベクトルの積 Krylov部分空間法 高々2N回行えば収束(理論) IDR(s)法 高々N+N/s回行えば収束(理論) A u v 行列ベクトル積 s>1の時、既存の方法より少ない １回 実際の数値計算においても 早く収束する！

IDR(s)法の一般化 IDR(s)法 GIDR(s,L)法[2] 1 L 計算量 N+N/s N+N/s 加速多項式 の次数 概要 歪対称に近い係数行列と相性が良くない [2]谷尾真明 杉原正顯, GIDR(s,L): 一般化IDR(s), 応用数理学会2008年度年会, pp. 411—412, 千葉.

既存のKrylov部分空間法と IDR(s), GIDR(s,L)の関係図 概要 IDR(1) = Bi-CGSTAB GIDR(1,L) = Bi-CGSTAB(L) 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s IDR(s) GIDR(s,L) 拡張 ？ 拡張

既存研究 Bi-CG(s)⇒Bi-CGSTAB(s)(=IDR(s)) (Sleijpen, 2008) 概要 加速多項式の次数 1 L 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s Bi-CG(s) IDR(s) =Bi-CGSTAB(s) GIDR(s,L) 拡張

本研究 GBi-CG(s)⇒GBi-CGSTAB(s,1)(=IDR(s)) ⇒GBi-CGSTAB(s,L)(=GIDR(s,L)) 概要 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s Bi-CG(s) GBi-CG(s) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) GIDR(s,L) GBi-CGSTAB(s,L) 拡張

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

IDR定理(1/2) IDR(s)法 定義: :　span{ } と直交する空間 : N次元ベクトル空間 (注意) は非ゼロのスカラー量

IDR定理(2/2) IDR(s)法 定理 (genericな条件の下で) -s次元 -s次元 ・

IDR(s)法の概要 N+N/sMATVECsで0に到達 IDR(s)法 残差ベクトルが入れ子になっている空間列 　　の中を　　　　　　　　　　　　と潜っていくように更新. 残差ベクトル N+N/sMATVECsで0に到達

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS) GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

IDR(s)法の不満点 係数行列が歪対称行列に近いと収束性悪化. IDR定理自身が持つ問題点. GIDR(s,L)法 加速多項式に相当 →高次にしたい！

数値実験による確認 IDR(3) (1次の加速多項式) vs Bi-CGSTAB(3) (3次の加速多項式) GIDR(s,L)法 3次元対流問題 離散化 ・125000×125000 ・歪対称行列に近い Bi-CGSTAB(3) IDR(3)

一般化IDR定理(1/2) GIDR(s,L)法 一般化 定義: :　span{ } : N次元ベクトル空間 高次 　　: 　　　　　　　　　　　　　　　

一般化IDR定理(2/2) GIDR(s,L)法 定理 (genericな条件の下で) -sL次元 -sL次元 ・

GIDR(s,L)のアルゴリズム GIDR(s,L)法 1. repeat until 14. end 2. For i = 0,1,…,L-1 15. Select 3. For j = 1,…,s 16. For i=1,…,L 4. Solve for 17. 5. 18. 19. 6. Compute 7. 20. end 8. end 21. 9. Solve for 22. 10. 23. end repeat 11. 12. 13. (k=0,…,i) 残差の更新の際に 補助ベクトルを導入している

数値実験による確認 IDR(3) (1次の加速多項式) vs GIDR(3,3) (3次の加速多項式) GIDR(s,L)法 3次元対流問題 離散化 ・125000×125000 ・歪対称行列に近い GIDR(3,3) IDR(3)

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

関係図 =IDR(1) =GIDR(1,L) =GBi-CG(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) M = : Mathematically equivalent 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

GBi-CG(s) Bi-CG法 Krylov 部分空間 反復 条件 shadow residual

GBi-CG(s) Bi-CG法の1反復 s. t. s. t. 残差が満たす性質

Bi-CG法の一般化 shadow residualの高次化 GBi-CG(s) Bi-CG法の一般化 shadow residualの高次化 Bi-CG GBi-CG(s) shadow residual N ・・・ N 1 s 条件 span

Bi-CG法の一般化 shadow residualの高次化 GBi-CG(s) Bi-CG法の一般化　 shadow residualの高次化 ML(s)BiCG(Yeung and Chan, 1999) Bi-CG(s) (Sleijpen, 2008) GBi-CG(s) 加速多項式の次数 1 L Bi-CG Bi-CGSTAB Bi-CGSTAB(L) s Bi-CG(s) GBi-CG(s) ML(s)BiCG shadow residual

GBi-CG(s)の1反復 GBi-CG(s) s. t. For j = 1,…,s set s. t. 残差が満たす性質 end

GBi-CG(s)法の計算量 total 2s × N/s = 2N MATVECs 1反復あたりの行列ベクトル積の演算 shadow residual ・・・・・・s 収束するのに必要な反復回数 ×s (genericな条件下で) total 2s × N/s = 2N MATVECs (Bi-CG法, Bi-CGSTAB法などと同じ)

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

拡張のアイデア =IDR(1) =GIDR(1,L) =GBi-CG(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) GBi-CGSTAB(s,L) Bi-CG⇒Bi-CGSTAB(L)と同じ構造を利用することによりGBi-CG(s)⇒GBi-CGSTAB(s,L)を達成！ 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

加速多項式の仕組み (Bi-CG) (1/3) =IDR(1) =GIDR(1,L) =GBi-CG(s) =Bi-CGSTAB(s) GBi-CGSTAB(s,L) 加速多項式の仕組み (Bi-CG) (1/3) Bi-CG⇒Bi-CGSTAB 　　　　　⇒Bi-CGSTAB(L) 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

GBi-CGSTAB(s,L) 加速多項式の仕組み (Bi-CG) (2/3) ・Bi-CG法における任意性 s. t. s. t.

加速多項式の仕組み (Bi-CG) (3/3) 残差 shadow residual GBi-CGSTAB(s,L) Bi-CG Bi-CG + 加速多項式 加速多項式 (fixed)

加速多項式の仕組み (GBi-CG(s)) (1/3) GBi-CGSTAB(s,L) 加速多項式の仕組み (GBi-CG(s)) (1/3) GBi-CG(s)⇒GBi-CGSTAB(s,1)(=IDR(s)) 　　　　　　　　⇒GBi-CGSTAB(s,L)(=GIDR(s,L)) M 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

加速多項式の仕組み (GBi-CG(s)) (2/3) GBi-CGSTAB(s,L) 加速多項式の仕組み (GBi-CG(s)) (2/3) ・GBi-CG(s)における任意性 For i = 1,…,s set s. t. 残差が満たす性質 end 35

加速多項式の仕組み (GBi-CG(s)) (3/3) GBi-CGSTAB(s,L) 加速多項式の仕組み (GBi-CG(s)) (3/3) 残差 shadow residual GBi-CG(s) GBi-CG(s) + 加速多項式 (fixed) 加速多項式

GBi-CGSTAB(s,L)=GIDR(s,L) 加速多項式の次数 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

GIDR(s,L)のアルゴリズム GBi-CGSTAB(s,L) 再掲 GBi-CGSTAB(s,L) 1. repeat until 14. end 2. For i = 0,1,…,L-1 15. Select 3. For j = 1,…,s 16. For i=1,…,L 4. Solve for 17. 5. 18. 19. 6. Compute 7. 20. end 8. end 21. 9. Solve for 22. 10. 23. end repeat 11. 12. 13. (k=0,…,i) 補助ベクトル ⇔ GBi-CG(s)における補助ベクトル

GBi-CGSTAB(s,L)の計算量 2N N+N/s k反復後の k反復後の残差 shadow residual Matvec ks k = [N/s]で収束 GBi-CG(s) + 加速多項式 (fixed) Matvec 0 計算量 Matvec ks + k N+N/s

発表の流れ 概要 IDR(s)法, GIDR(s,L)法の導出 (応用数理年会) GBi-CGSTAB(s,L)法の導出 (RIMS研究集会) GBi-CG(s)法 GBi-CGSTAB(s,L)法(=GIDR(s,L)法) まとめ

結論 =IDR(1) =GIDR(1,L) =GBi-CG(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) GIDR(s,L)法を, Bi-CG法の拡張として解釈出来た(GBi-CGSTAB(s,L)法) 1 L Bi-CG Bi-CGSTAB =IDR(1) Bi-CGSTAB(L) =GIDR(1,L) s ＝Bi-CG(s) =GBi-CG(s) (ML(s)BiCG) IDR(s) =Bi-CGSTAB(s) =GBi-CGSTAB(s,1) (ML(s)BiCGSTAB) GIDR(s,L) GBi-CGSTAB(s,L) M M M M M

今後の課題 前処理を含めたGBi-CGSTAB(s, L)法の提案. パラメータsとLの決定法. 加速多項式の次数 1 L Bi-CG 1 L Bi-CG Bi-CGSTAB Bi-CGSTAB(L) s Bi-CG(s) GBi-CG(s) ML(s)BiCG IDR(s) Bi-CGSTAB(s) GBi-CGSTAB(s,1) ML(s)BiCGSTAB GIDR(s,L) GBi-CGSTAB(s,L) shadow residual

実験環境 MATLAB7.30 で反復終了 加速多項式の次数 1 2 3 4 shadow residual IDR(1) Bi-CGSTAB(1) GBi-CGSTAB(1,1) Bi-CGSTAB(2) Bi-CGSTAB(3) Bi-CGSTAB(4) IDR(2) GBi-CG STAB(2,2) IDR(3) STAB(3,3) IDR(4) STAB(4,4) shadow residual

実験その1　3次元対流問題 離散化 ・125000×125000 ・歪対称に近い

実験その1 Bi-CGSTAB(L) Bi-CGSTAB(1) shadow residual Bi-CGSTAB(2) 加速多項式の次数 shadow residual 1 2 3 4 Bi-CGSTAB(2) Bi-CGSTAB(3) Bi-CGSTAB(4)

実験その1 IDR(s) IDR(1) shadow residual IDR(2) IDR(4) IDR(3) 加速多項式の次数 1 2

実験その1 GBi-CGSTAB(s, L) GBi-CGSTAB(1,1) shadow residual GBi-CGSTAB(2,2) 加速多項式の次数 shadow residual 1 2 3 4 GBi-CGSTAB(2,2) GBi-CGSTAB(3,3) GBi-CGSTAB(4,4)

実験その1 s = L = 3 shadow residual GBi-CGSTAB(3,3) IDR(3) Bi-CGSTAB(3) 加速多項式の次数 shadow residual 1 2 3 4 GBi-CGSTAB(3,3) IDR(3) Bi-CGSTAB(3)

実験その2　　raefsky2 (matrix market より) x b A 3242×3242

実験その2 Bi-CGSTAB(L) Bi-CGSTAB(1) shadow residual shadow residual 加速多項式の次数 shadow residual shadow residual 1 2 3 4 Bi-CGSTAB(2) Bi-CGSTAB(3) Bi-CGSTAB(4)

実験その2 IDR(s) shadow residual IDR(4) IDR(3) IDR(2) IDR(1) 加速多項式の次数 1 2

実験その2 GBi-CGSTAB(s,L) GBi-CGSTAB(1,1) shadow residual GBi-CGSTAB(4,4) 加速多項式の次数 shadow residual 1 2 3 4 GBi-CGSTAB(4,4) GBi-CGSTAB(3,3) GBi-CGSTAB(2,2)

実験その2 s = L = 3 shadow residual IDR(3) GBi-CGSTAB(3,3) Bi-CGSTAB(3) 加速多項式の次数 shadow residual 1 2 3 4 IDR(3) GBi-CGSTAB(3,3) Bi-CGSTAB(3)

IDR定理の空間列 に関する定理(2008, Sleijpen)

一般化IDR定理の空間列 に関する定理

IDR(s)法に対する 2つの解釈 (2008, Sleijpen) 加速多項式を取り払った 空間と残差 IDR定理 -s次元 -s次元

GIDR(s,L)に対する 2つの解釈 加速多項式を取り払った 残差の動き 一般化IDR定理 -sL次元 -sL次元