複素対称連立一次方程式に対する Block Krylov
部分空間反復法の開発と並列固有値解法への応用

多田野 寛人筑波大学システム情報系
A01 高橋班
新学術領域研究「コンピューティクスによる物質デザイン」!
研究会!

東京大学本郷キャンパス!
年3月月17日日
2012年

§1はじめに!

大規模計算環境に適したアルゴリズムの必要性
背景§1 はじめに物性物理，分子軌道計算などの様々な分野で大規模!
固有値問題が現れ，高速に解くことが望まれている．!
!
近年「京コンピュータ」などの超大規模並列環境が登場超大規模並列環境に適した固有値解法!
アルゴリズムが必要！
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日日

 2 !

並列固有値計算アルゴリズム：Sakurai-Sugiura法
法§1 はじめに一般化固有値問題： Fv = λSv
F = FT , S = ST ∈ R n×n, V ∈ R n×L

Im!

周回積分による行列の定義

1Gk ≡2πi�Γkz (zS − F) SVdz
−1

N点点台形則による近似

1 N−1k−1G̃ k =wz(zS−F)
SVj jj∑N j=0z1λm+1λ1O!

Γλ2 . . .

λmz0 λm+2z N−1 Re!
zj●：固有値，●：積分点．!

・複数本の右辺ベクトルを持つ連立一次方程式の求解が必要．
・係数行列 z j S − F は複素対称行列となる．
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日日

 3 !

本研究の目的§1 はじめに問題点!

複数右辺ベクトルをもつ連立一次方程式の求解が!
法の演算時間の大半を占める!
SS法

目的複数右辺ベクトルをもつ複素対称連立一次方程式!
を高速に解くための Block Krylov 部分空間反復法!
を構築する．!
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年年3月月17日
日

 4 !

§2 複素対称行列に対する!
! !Block Krylov 部分空間反復法

複数の右辺ベクトルをもつ複素対称連立一次方程式§2 複素対称行列に対するBlock Krylov部
部分空間反復法

解きたい連立一次方程式
AX = B

ここで，A ∈ C n×n：n 次複素対称行列（ A = AT � AH ），!
�� (1) (2)�
� (1) (2)(L)(L)
, B = b ,b ,...,b
X = x ,x ,...,x

右辺が1本本の場合の Krylov 部分空間反復法!COCG法法!van der Vorst, et al., (1990)!

法!CSQMR法!Freund, (1992)!

法!COCR法!CSYM法法!!

!Sogabe, et al., (2007)!

!Bunse-Gerstner, et al., (1999)!

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日日 6 !

Block Krylov 部分空間反復法§2 複素対称行列に対するBlock Krylov部
部分空間反復法

Block Krylov 部分空間反復法の特長

1. !L 本の方程式をまとめて解くことができる．!
本ずつ解いた場合よりも!
2. !Krylov 部分空間法で 1本!少少ない反復回数で残差が収束することがある．!

複素対称行列用の Block Krylov 部分空間反復法で!
SS法
法の計算時間を短縮できる可能性がある
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日日 7 !

COCG法法のBlock版版アルゴリズム§2 複素対称行列に対するBlock Krylov部
部分空間反復法

COCG法
法

Block COCG法
法

x0 ∈ C n is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 ,For k = 0, 1, . . . , until �r k �2 ≤ ε�b�2 do:
α k = rTk r k / pTk Apk ,X0 ∈ C n×L is an initial guess,
Compute R0 = B − AX0 ,Set P0 = R0 ,For k = 0, 1, . . . , until �R k �F ≤ ε�B�F do:
Solve (PTk AP k )α k = RTk R k for α k ,pk+1 = r k+1 + β k pk ,End forP k+1 = R k+1 + P k β k ,End forx k+1 = x k + α k pk ,r k+1 = r k − α k Apk ,β k = rTk+1 r k+1 /rTk r k ,X k+1 = X k + P k α k ,R k+1 = R k − AP k α k ,Solve (RTk R k )β k = RTk+1 R k+1 for β k ,COCG法
法のBlock 版への拡張は容易にできる!
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日日

 8 !

Block Krylov部
部分空間反復法の問題点§2 複素対称行列に対するBlock Krylov部
部分空間反復法

Block Krylov 部分空間反復法の特長

1. !L 本の方程式をまとめて解くことができる．!
本ずつ解いた場合よりも!
2. !Krylov 部分空間法で 1本!少少ない反復回数で残差が収束することがある．!
!

3. !右右辺ベクトル数 L が多い場合は，残差が収束せず!
!発発散することがある．!
!

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月
月17日日

 9 !

Block COCG法
法の数値的不安定性の原因§2 複素対称行列に対するBlock Krylov部
部分空間反復法X0 ∈ C n×L is an initial guess,
Compute R0 = B − AX0 ,Set P0 = R0 ,For k = 0, 1, . . . , until �R k �F ≤ ε�B�F do:
Solve (PTk AP k )α k = RTk R k for α k ,X k+1 = X k + P k α k ,R k+1 = R k − AP k α k ,Solve (RTk R k )β k = RTk+1 R k+1 for β k ,P k+1 = R k+1 + P k β k ,End forα k , β k を求めるために，!

小規模連立一次方程式を!
解く必要がある．

不安定性の原因!

行列 P k , R k の線形独立性が!

失われ，小規模連立一次方程式!
が悪条件問題になってしまう！!

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  10 !

数値的不安定性の改善§2 複素対称行列に対するBlock Krylov部
部分空間反復法

Hermite行
行列用のBlock CG法
法では修正方向行列 Pk，または!
残差行列 Rk を直交化することで数値的不安定性を改善 [1]．
．

Block COCG法
法も同様の手法で安定化できないか？!
[1] !A. A. Dubrulle, Retooling the method of block conjugate gradients, !
!Electronic Transactions on Numerical Analysis, 12 (2001), pp. 216—233.
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  11 !

Block COCG法
法に直交化を施したアルゴリズム§2 複素対称行列に対するBlock Krylov部
部分空間反復法

Block COCGdQ法
法

Block COCGrQ法
法

X0 ∈ C n×L is an initial guess,
Compute R0 = B − AX0 ,Compute S0 ξ0 = R0 ,For k = 0, 1, . . . , until �R k �F ≤ ε�B�F do:
TSolve (STk AS k )α k = ξ−TR
R for α k ,k kk

X0 ∈ C n×L is an initial guess,
Compute Q0 ξ0 = B − AX0 ,Set S0 = Q0 ,For k = 0, 1, . . . , until �ξ k �F ≤ ε�B�F do:
Solve (STk AS k )α k = QTk Q k for α k ,X k+1 = X k + S k α k ,R k+1 = R k − AS k α k ,TSolve [(RTk R k )ξ−1]β=RRfor β k ,kk
k+1 k+1S k+1 ξ k+1 = R k+1 + S k β k ,End for

修正方向行列直交化版!

X k+1 = X k + S k α k ξ k ,Q k+1 τ k+1 = Q k − AS k α k ,ξ k+1 = τ k+1 ξ k ,Solve (QTk Q k )β k = τTk+1 QTk+1 Q k+1 for β k ,S k+1 = Q k+1 + S k β k ,End for

残差行列直交化版!

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  12 !

§3数値実験!

解く問題とパラメータ設定
解く問題
DNA8塩塩基対の分子軌道計算で現れる!
一般化固有値問題 Fv = λSv をSS法
法で!
解く際に現れる連立一次方程式：§3 数値実験(z0 S − F)Y = SV.

ここで，n = 1,980, nnz(F) = 2,752,198,!
nnz(S) = 693,550.

(a) 行列 F.

(b) 行列 S.

Fig. 1. 非零構造図．

パラメータ設定・右辺ベクトル数 L は L = 1, 4, 8, 16, 32, 64 の6ケケースとした．!
・係数行列のパラメータ z0 は z0 = -1.0 + 1.0 " 10-3 i．
．!
・右辺項 SV は乱数で与えた．!
・反復法の前処理として近似係数行列の完全分解前処理を使用．!
Block COCG法
法, Block COCGdQ法
法, Block COCGrQ法
法の性能を比較した．!
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  14 !

近似係数行列の完全分解前処理§3 数値実験

Step 1. 近似係数行列 AC の生成
A = [ai j ] −→ AC = [ ãi j ]
|ai j | or i = j ai j , |ai j | ≥ δ · maxi, j
ここで，ãi j = 
 0,otherwise! : カットオフパラメータ ( 小さいほど AC は A に近い )．
．

Step 2. 近似係数行列の完全分解
行列 AC を完全Cholesky分分解し，前処理行列として採用! を小さくすると少ない反復回数で解が得られる．!
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  15 !

各解法の相対残差履歴の比較（その１）(a) Block COCG法
法.

(b) Block COCGdQ法
法.

L を増やすと全ての解法で!
反復回数が減少した．!
Fig. 2. ! = 1.0 " 10-3 での相対残差履歴．!
■：L = 1, !■■：L = 4, !■■：L = 8，，!
■：L = 16, !■■：L = 32, !■■：L = 64.

(c) Block COCGrQ法
法.

各解法の相対残差履歴の比較（その２）(a) Block COCG法
法.

(b) Block COCGdQ法
法.

Block COCGrQ法
法では L を!
大きくしても残差が収束．!
Fig. 3. ! = 1.0 " 10-2 での相対残差履歴．!
■：L = 1, !■■：L = 4, !■■：L = 8，，!
■：L = 16, !■■：L = 32, !■■：L = 64.

(c) Block COCGrQ法
法.

SS法
法に対する Block Krylov 部分空間反復法の適用§3 数値実験

解く問題

上皮成長因子受容体2量量体の分子軌道!
計算で現れる一般化固有値問題：!

Fv = λSvn = 96,234, nnz(F) = 456,807,648,!
nnz(S) = 53,146,470.

(a) 行列 F.
Fig. 4.

(b) 行列 S.

非零構造図．

実験条件
・中心 - 0.450415，，半径 1.85815"10-3（虚軸方向は半径を0.1倍倍）!の楕円内に存在する固有値・固有ベクトルを計算．!
・円周上の積分点数は N = 32 （方程式数は16本本）とし，!
右辺ベクトル数は L = 32 とした．!

法，Block COCGrQ法
法で連立一次方程式を解いた．!
・COCG法反復の停止条件は相対残差が 10-10 以下になった場合．

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  18 !

SS法法に対する Block Krylov 部分空間反復法の適用§3 数値実験

Table 1. 実験環境（T2K-Tsukuba 8ノノード）．
CPUMemory
CompilerCompile option
MPI library
Direct solver

AMD Opteron Quad Core 8356 2.3GHz × 4
32.0GBytes
Intel Fortran ver. 11.1-O3 -xO
MVAPICH2 ver. 1.7
SuperLU DIST ver. 3.0

・128コ
コアに対して Flat MPI で並列化!

つの連立一次方程式に対して 2ソソケット を割り当て計算!
・1つ・SuperLU_DIST は反復法の前処理として使用!

近似係数行列の完全分解前処理（! = 1.0 " 10-5）を用いた．

新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  19 !

Krylov 部分空間法と Block Krylov 部分空間法の性能比較§3 数値実験

積分点毎の反復回数!
ばらつきが大きい

反復回数はほぼ一定!

Fig. 5. COCG法
法とBlock COCGrQ法
法の反復回数．!
●：1本本の方程式を解くのに要したCOCG法
法の平均反復回数，!
!●!●●：32本本解くのに要したBlock COCGrQ法
法の反復回数．!
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  20 !

Krylov 部分空間法と Block Krylov 部分空間法の性能比較§3 数値実験■：行列データ分散!
■：前処理行列生成!
■：行列ベクトル積!
■：前進・後退代入!
■：その他!
□：待ち時間!
■：固有対計算

Fig. 6. COCG法
法を用いた場合のSS法
法の計算時間．
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  21 !

Krylov 部分空間法と Block Krylov 部分空間法の性能比較§3 数値実験■：行列データ分散!
■：前処理行列生成!
■：行列ベクトル積!
■：前進・後退代入!
■：その他!
□：待ち時間!
■：固有対計算

Fig. 7. Block COCGrQ法
法を用いた場合のSS法
法の計算時間．
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  22 !

Krylov 部分空間法と Block Krylov 部分空間法の性能比較§3 数値実験(a) COCG法
法．

(b) Block COCGrQ法
法．

Fig. 8. SS法
法の計算時間比較．
新学術領域研究「コンピューティクスによる物質デザイン」研究会 2012年
年3月月17日
日  23 !

§4

まとめと今後の課題

まとめ
1.! 複素対称行列に対するBlock Krylov部
部分空間反復法を!

!構構築した．前処理がより不完全な場合は残差直交化を!

!施施した Block COCGrQ法
法が有効であった．!

2.! SS法
法にBlock Krylov 部分空間反復法を適用することで!
!高高速化ができ，計算時間のばらつきを抑えられた．!

今後の課題
1. !残
残差直交化版である Block COCGrQ法
法がなぜ他の解法!
!よよりも安定性が高いのかを調べる．!

更なる安定化は可能かどうかを検証する．!
2. !更