Algorithmic and Machine issues on the

格子ＱＣＤのための線形計算
（連立一次方程式、固有値問題）について
石川健一(広島大理)
アルゴリズムによる計算科学の融合と発展
@CCS, 筑波大学, ２００９年４月２２日・２３日
1
1. 目次
 ２．格子ＱＣＤについて
 ３．ハイブリッドモンテカルロ法
 ４．格子化クォークのタイプ
 ５．クォーク伝播関数ソルバー
 ６．固有値問題
 ７．まとめ
2
２. 格子ＱＣＤについて
 QCDを解くためには
時空を格子化した有限自由度の
格子QCDという方法を用いる。
K.G.Wilson (1974)
連続時空
連続時空上の場の変数 → 格子上の場の変数
格子化した時空
3
２. 格子ＱＣＤについて
 分配関数（ユークリッド化さ
れた経路積分）
Z    dUd e  S (U , )
   dUe  H (U )
作用
 (n)
有効作用
 物理量の期待値
1
O 
Z
 H (U )
dU
O
(
U
)
e

 (n)
U  (n)
有効作用を重みとする多重積分＝＞モンテカルロ積分
統計力学で用いられてきた方法
4
２. 格子ＱＣＤについて
特徴、
自由度がとても大きい
16^4 ～ 32^4 格子
 グルーオン：８ｘ４ｘ（16^4～32^4）=200万～3000万自由度(実数勘
定)
 クォーク：3x4x （16^4～32^4） =160万～1500万自由度(実数勘定)
系統誤差
 格子化に伴う誤差:格子間隔 a = 0.1 fm ～ 0.06 fm → 0 fm
 有限体積による誤差: 核子が大体収まる大きさ L= 3fm～4fm
→∞
 クォークの質量が現実世界と異なるための誤差
m_q=40MeV～100MeV → m_ud = 2MeV～ 10MeV
典型的な計算時間１～２年（論文、予算、、、、）
どのようにしてモンテカルロ積分のための配位を生成する
か？
5
3. ハイブリッドモンテカルロ(HMC)法
 格子QCD分配関数（Nf=2:クォークが２種類）
Z    dPdUd d e
†
 H [ P ,U ,† , ]
,

1
H [ P,U ,  ,  ]   Tr[ P (n) 2 ]  SG [U ]  † D[U ]†
2 n,
†
 D[U ]
1
1

P、  : 補助場、
SG [U ] : U のみの関数、局所的
D[U ] : U に依存する大規模疎行列。
格子点に関する
 (n)
U  (n)
差分で出来ている。
ユークリッド化されているので、統
計力学の分配関数と同等の計算
になる。
6
3.ハイブリッドモンテカルロ(HMC)法
 クォーク部分の計算が大変
1 2 

Z Q [U ]    d d exp  D[U ]  :クォーク部分の分配関数


†


   d† d exp  S Q [† ,  ]
D[U ] : 差分化された Dirac 演算子
φで積分すると
Z Q [U ]  det[ D[U ]]
2
真空からクォーク－反ク
オーク対が対生成、対消
滅している効果を表して
いる。真空偏極。
 モンテカルロ重み＝行列式＝補助場による分布
 D[U]の逆を含む：非局所的
7
3.ハイブリッドモンテカルロ(HMC)法
 Hybrid Monte Carlo（HMC）法
Z    dPdUd d e
†
 H [ P ,U ,† , ]
,
1
H [ P,U ,  ,  ]   Tr[ P (n) 2 ]  Seff [U , † ,  ]
2 n,
†
 Exp(-H) の分布に従うU,P,φ を HMC法で生成
 Uのサンプル{ U(1),U(2),U(3),…..,U(N)}
 物理量の期待値(物理量は U のみの関数 O[U])
†
1
 H [ P ,U ,† , ]
O    dPdUd d O[U ] e
Z
1 N
  O[U ( k ) ] : サンプル平均
N k 1
8
3.ハイブリッドモンテカルロ(HMC)法
 HMC法
Z    dPdUd d e
†
H [ P,U , † ,  ] 
 H [ P ,U ,† , ]
,
1
Tr[ P (n) 2 ]  Seff [U , † ,  ]

2 n,
に対して、マルコフチェーンモンテカルロ法の一種ＨＭＣ法を適用
 分子動力学(Molecular Dynamics) (MD)
H [ P, U , † ,  ]  Tr[ P 2 ] / 2  S eff [U , † ,  ]
の力学はを離散化して分子動力学法を用いて
時間発展させる。（
の微分を差分化）
 メトロポリス法(Metropolis)
 分子動力学ではエネルギーは差分化のため保存しない
 H [ P ,U , ]
 分布がほしい分布 e
からずれる
分子動力学法につ
いての発展は今回
は省略
P[(P,U )  ( P' ,U ' )]  min(1, e H  H ' )
で補正する
9
3.ハイブリッドモンテカルロ(HMC)法
 Hybrid Monte Carlo（HMC）法
計算が困難な点
 分子動力学法の部分で
D[U ]x  b,  x  ( D[U ])1 b
の大規模連立１次方程式を何回も解く必要がある。
幸いＤ［Ｕ］の要素はほとんどゼロ（大規模疎行列）であるので、
反復法
を用いて連立方程式を解く。
 しかし特にクォークの質量が小さいと、D[U]がゼロに近い固有
値を持つため連立方程式を解くのが難しくなる。
 クォークの質量が小さいと分子動力学の力が大きくなる。分子
動力学の時間刻みを小さくする必要がある。
Z
（全体のコスト）


  1  , Z  2～3
 mq 
10
行列式の評価、
HMC法の改良、
分子動力学法に
ついての発展は
今回は省略
3.ハイブリッドモンテカルロ(HMC)法
 問題をまとめると、
クォークのタイプによって D[U] がいろいろある
モンテカルロ法による経路積分の評価で困難が生じるのは
クォークの寄与を取り入れる部分。
Z Q [U ]  det[D[U ]]
   de
† ( D[U ]) 1 
e
 S eff [U ]
Uの生成には、 D[U] x = b の連立一次方程式を大量に解く
必要がある
D[U ]x  b,  x  ( D[U ])1 b
クォークソルバー
これまでさまざまな前処理方法が開発されてきた
最近の発展について述べたい。
11
４.格子化クォークのタイプ
 離散化には任意性がある
 良い性質を持つものほど計算コストがかかる
 クォークの持つ対称性
 ゲージ対称性：格子上でも死守すべき対称性
 ポアンカレ対称性: 一部破れるけどなるべく保持したい対称性、ユーク
リッド化後：離散並進、離散回転、反転対称性
 カイラル対称性： D 5   5 D  0
格子上で厳密に実現することは不可能(ニールセン-二宮の定理)
 カイラル対称性の保持の仕方による分類
 Kogut-Susskind 型：４つの同じ質量をもつクォーク（４フレーバー）を同時
に扱う。U(1)xU(1)カイラル対称性を保持
 Wilson-Dirac型：１つのクォークから扱える。カイラル対称性は無い
 Overlap型、Domain-Wall 型：１つのクォークから扱える。カイラル対称性
は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて
10-100倍)
12
４.格子化クォークのタイプ
 Kogut-Susskind(KS) 型：４つの同じ質量をもつクォーク（４フレーバー）を
同時に扱う。U(1)xU(1)カイラル対称性を保持
 最も簡単なD[U]の構造を持つ。現実世界のクォークは６種類で質量もばらば
ら。軽いのはアップとダウンの２種類のクォーク＝４つ同時では扱えない。
スキップします。
13
４.格子化クォークのタイプ
 Wilson-Dirac型：１つのクォークから扱える。カイラル対称性は無い。
 Ｄ[U]は一階の差分演算子

M (n, m)   1      U  (n)   n ˆ ,m  1      U† (m)   n ˆ ,m
4
 1

DW [U ](n, m)  1 1   n,m   M (n, m)
†

D


(
D
)
,
5
W
5
W

HW   5 DW はエルミート行列
 クォーク質量 m ∝ 1/κ
 カラー(3),スピン(4)の自由度が各格子点にある
 ブロックサイズ１２のブロックの帯が並んだ疎行列
M の固有値分布 U=1,4^4
Witzel,Takeda,Wolff,
hep-lat[arXiv:0709.4648]
14
４.格子化クォークのタイプ
 Overlap型、Domain-Wall 型：１つのクォークから扱える。カイラル対称性
は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて
10-100倍)
DOV  m  1   5 signHW (M 0 )
HW (M 0 )   5 ( DW  M 0 )




質量 m
Wilson-Diracクォークの行列 DW のサイン関数を含む
Dov は密な行列
DOV  5  5 DOV  0
対称性：質量m=0 でGinparg-Wilson関係を満たす   (1  D )
5
OV
5

 5 signHW (M 0 ) はユニタリー行列
 5 signH W ( M 0 )  

 5 HW (M 0 )
HW (M 0 ) 2
DW ( M 0 )
DW ( M 0 )† DW ( M 0 )
Taken from
Talk slid by R.Edwards
15
４.格子化クォークのタイプ
 まとめ、
 Dx = b は反復法で解く
 係数行列 D の型
 Wilson-Dirac型
 １階差分型：さまざまな前処理が適用されてきた
 Red-brack(2色), ILU,SSOR, ２色ブロック … etc.
 Dx = b を直接解く: Non-Hermitian sovler,
BiCGStab, GMRES などが使われてきた
 (D†D)x=D†b を解く: positive Hermitian solver,
が使われる。条件数悪し。
 γ5Dx=γ5D を解く: indefinite Hermitian solver,
MINRESなど、、
 前処理つき BiCGStab, GMRES 系統が良く使われる
CG
16
４.格子化クォークのタイプ
 Overlap 型
 Dをベクトルにかけるときに sign 関数の計算が必要
 Dx=b 自体は Hermite化後 CG などで解く
 Sign関数の計算
T
HW Vk  Vk Tk   k 1vk 1ek ,
 Arnoldi / Lanczos 分解による方法
signHW b  Vk signTk Vk b  ....
†

j
signHW b  H w   2
 j 0 H  
w
j

 部分分数/連分数展開＋補助場の導入による方法
=> ５次元化ソルバーの入れ子を解消
 0 x

1
 Pade / Rational ：部分分数展開近似による方法

 1
1

signx    0 x 
 Schur
2
1 x 



2 x  3
.....

 1 x

  2 n 1 x
 2n




 2n 

 2 n x 
 HW や DW の固有値範囲の情報を用い sign関数の近似の精度をコントロー
ルする必要がある。
Overlap 型についてはこれ以上次官の都合上言及しない
日本では JLQCD collaboration でくわしく研究されている
KEK 松古、金児、橋本…..
17

b


５.クォーク伝播関数ソルバー
 Ｄｘ＝ｂを速く解きたい
 格子サイズ: 32^3x64 (現状 O(10)TFlops)
=> 256倍： 128^3x256 (O(2)PFlops)
 次元: １２ｘ[格子サイズ]= 2500万＝＞ 64億
(1)前処理：
Wilson型に関しては構造が簡単（Red-Black, ILD, DomainDecomposition etc.）, Overlap型については難しい⇒5D化で構造が簡単に [省
略]
(2)固有値固有ベクトルの情報：
(3)単精度加速：
(4)アクセラレータ：
Deflation/Multigrid
18
5.クォーク伝播関数ソルバー
(2) Deflation technique
 LQCD では何回も Dx=b を解く
Dx(i )  b (i ) , i  1,2,3
D (i ) x (i )  b (i ) , i  1,2,3, D (i )  D (i 1) ,
 右辺ベクトルを取り替えながら、係数行列は一定。または、係数行列が
少しづつ変化しながら連鎖的に解く必要がある
 Quark propagator
 Solver in HMC trajectory
 ブロック化ソルバー：係数行列固定、右辺ベクトル多数
 Deflation technique: 係数行列の固有ベクトルの情報により前処理
 Deflation remove/suppress small eigenspace of D.
19
5.クォーク伝播関数ソルバー
(2) Deflation technique
（Solve Ａｘ＝ｂ・・・(1) ）
 行列A ：小さい固有値に属する p次元部分空間を持つとする
c,u : この部分空間に属する以下の
条件を満たすベクトルとする
U p  (u1 , u2 ,, u p )
C p  (c1 , c2 ,, c p )
AU p  C p
 c から射影演算子を定義
Q  I p  U pC A
†
p
 PAy  C pC†pb
 Pb  C pC†pb  b
PA  AQ
 Pを使って問題を前処理
( PA) y  Pb
x  Qy  U pC†pb
 Ax  AQy  AU pC†pb
C†pC p  I p
P  I p  C pC†p
 式(1)の解 x は式(2) の解 y を用い
て以下のように書かれる
(2)
 “Cp”はどうやって作る？
 反復には P の掛け算が伴う
(overhead)
20
5.クォーク伝播関数ソルバー
(2) Deflation technique (cont’d)
Many works by
[Luescher, JHEP07(2007),hep-lat/0710.5417;
A.Stathopoulos, K.Orginos, hep-lat/0707.0131;
W.Wilcox, PoS(LATTICE2007),hep-lat/0710.1813;
A.Abdel-Rehim,R.B.Morgan,W.Wilcox,PoS(LATTICE2007);
R.B.Morgan,W.Wilcox,math-ph/0707.0505,math-ph/0405053;
M.L.Parks, E.De Sturler et al, SIAM J. on Sci.Comp. 28(2006)1651
LATTICE2008: Poster by Abdel-Rehim, Talk by Wilcox]
More details see Wilcox @Lat2007.

正確な不変部分空間の計算はとてもコストが高い
低モード密度 ∝ O(V), モード計算コスト O(V)
=> 全体で O(V^2) の問題
(a) Dx=b を解きながら、同時に不変部分空間を求め使いまわす
 GCRO-DR: Parks & Sturler, used by PACS-CS collab.
 GMRES-DR,GMRES-E..: Wilcox, Morgan & Abdel-Rehim
(b)低モード固有ベクトルは滑らかに局在。ブロック分割可能？
 Luscher’s Domain decomposed subspace blocking with local
coherency. (Used also in HMC algorithm)
21
5.クォーク伝播関数ソルバー
(a)固有空間の計算と連立一次方程式を同時に解く.
Very effective for few Near zero modes / negative eigen modes case.

Near zero modes case
[Wilcox, LAT2007]
 First equation or few equations are solved with GMRES-DR.
Once the subspace converged, change solver with GMRES-proj,
or Deflated solver.
 Normal GMRES stagnates [dot-dot-dashed line]
 Solver with Deflation/Projection converges. [other lines]
22
5.クォーク伝播関数ソルバー
(b)ブロック化された基底ベクトルで元の行列の逆を取り近似固
有空間を構成する.
 基底ベクトルの構成方法が重要 [Luscher, JHEP07(2007)081]
 Low mode ベクトルはブロックベクトルで近似できる
a low mede vector: ( x) 
blocks N



c

 ( j ,  ) j ( x)
j 1
23
5.クォーク伝播関数ソルバー
[Luscher, JHEP07(2007)081]
(b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構
成する.






C   j ( x) : j  1,, N ,   all domain blocks
 基底ベクトル（C）の構成
ランダムベクトルに1/D を数回か
けて低モードを増幅
そのベクトルをサブブロックに分
割
分割されたベクトルを正規直交化
 0 ( x   )
,
 0 ( x   )
 j ( x)  
(i )†   j   ij 
Deflation の射影演算子 (B^-1)を
含む.
P  1  DCB1C† , Q  1  CB 1C† D, B  C† DC,
PD  DQ, P 2  P, Q 2  Q,

B : ブロック化された部分空間でのD[U]の表現（＝サブグリッド上に射影さ
れた） B(i, ; j , )  (  )†  ( D  )
i
j
 B(i, j ,  )  ,   B(i, j , ,  )   ˆ ,  B(i, j , ,  )   ˆ , 
4
 1
24
5.クォーク伝播関数ソルバー
(2) MultiGrid Solver

MultiGrid solver also removes
critical slowing down.
Choice of subspace basis is
important. (Prolongator)

v( x) : a random vector
w  (1  D† D) v v  ( D† D) 1 v :
low mode enhanced,then
w is blockedas


C    ( x) :   all domain blocks
 0 ( x   )
,
 0 ( x   )
  ( x)  
QCD, 16^3x32
[Brannick,Brower,Clark,Osborn,Rebbi,
PRL100(2008)]
(  )†      
C is used for P rolongator.

射影のための基底ベクトルの選び方
が重要（Luscher の基底と同じ）
25
5.クォーク伝播関数ソルバー
(3) 単精度加速
 単精度の利点: 有効バンド幅、キャッシュサイズ、（レジス
ターの数）などが倍になる＝効率（実性能/理論性能）が高い
 Intel 64/AMD 64; Single prec. > Double prec.
 Cell PS3/GPGPU; Single >> Double.
 For Wilson kernel : 3Byte/Flop => 1.5 Byte/Flop
 連立一次方程式の解法を加速するのに使えないか？
 反復改良法/Richardson反復のテクニックで可能
[Buttari,Dongarra,Kurzak,Luszczek,Tomov, AMC Trans.Math.Soft]
 一般の反復法にも組み入れることが出来る（可変前処理）
• 前処理として単精度ソルバーを使う
• 残差は倍精度で正しくなるように組む
• FGMRES, GCR, CG, BiCGStab….. (with flexible prec.)
26
5.クォーク伝播関数ソルバー
(3) 単精度加速
Nested-BiCGStab (櫻井-多田野)
Outer solver : BiCGStab (D.P.)
Inner solver: BiCGStab (S.P.)= Outer solverに対する可変前処理
[PACS-CS collab.]
 Intel 64, SSE3 使用+リンク再構築で全て倍精度の時より速度
が倍に（演算量は増えているが）
27
5.クォーク伝播関数ソルバー
(４) アクセラレーター
 単精度加速と組み合わせるとさらに有効なもの
 GPGPU, CellB.E. は単精度計算が圧倒的に速い
 GRAPE => QCD?
 NVIDIA GT200 (Tesla 10series)
 240 SP (SP cores), 30 DP cores
 ～1,000(or 600)Glops(SP), ～90GFlops(DP)
 We expect > 60 GFlops(SP) for QCD kernel.
(assuming 10% efficiency)
28
5.クォーク伝播関数ソルバー
(４) アクセラレーター
 NVIDIA CUDAを使った例（石川-尾崎）
 CPU: Core2Duo@3GHz
 GPU: GeForce GTX 280
 O(a)-improved Wilson-Dirac
 Red/black site prec’d
 Nested-BiCGStab
 Time
 CPU only: 184sec,
 CPU: 1.9GFlops
 CPU+GPU: 8.6 sec,
 CPU: 1.7GFlops
 GPU: 58GFlops
 GPU D mult.: 102GFlops
この場合 184/8.6 = 21倍の高速化!
10倍以上速くなるのは大きい
29
5.クォーク伝播関数ソルバー
(４) アクセラレーター
 NVIDIA CUDAを使った例（Clark et al. Work shop@CCS, 10-12
March, 2009）
 Nvidia GPU でさらに加速する方法
 単精度＝＞半精度（16bit）化してさらにバンド幅を稼ぐ
 符号付16bit整数をノーコストで[-1.0,1.0] 単精度へ自動変換できる機能
を利用（Texture Fetching）
 有効バンド幅と、有効 Texture Cache サイズが倍、レジスタは32bit 浮
動小数点のものしかないので演算は単精度
 メモリとのやり取りはTexture Fetching (ロード)と整数への丸め(ストア)
30
5.クォーク伝播関数ソルバー
(４) アクセラレーター
 NVIDIA CUDAを使った例
[M.Clark et al. Work shop@CCS, 10-12 March, 2009]
 半精度(16bit int)でさらに２
割性能上昇
 そのほかの工夫でさらに加
速
• ゲージ固定
• ＳＵ（３）行列を実８パラメー
タで保存
 GPGPU:演算が圧倒的に
速いので、バンド幅、キャッ
シュの有効利用のほうが重
要
31
5.クォーク伝播関数ソルバー
(４) アクセラレーター
 Cell B.E.での計算
 大規模 Cell B.E.クラスタ：
[Spary,Hill,Trew hep-lat/0804.3654;
S.Motoki & A. Nakamura Lat2007;
F.Belletti et al. LAT2007]
[H.Baier at al. ,Lat2008, arXive:0801.1559
 QPACE (QCD Parallel computing on Cell B.E.) project (EU)
 2048 PowerXCell8i, + カスタムネットワーク3Dトーラス(FPGA)
 2009春完成？
 ＧＰＵクラスター:
 National Taiwa Univ. 16nodes+16TeslaS1070: 64TFlops
 KEK, 4 nodes + 4 Tesla S1070: 16TFlops
 ….?
 多体系の計算に対する応用がGPGPUでも先行（GRAPE）
 AMD/ATI での計算 [V.Demchik and A. Strelchenko arXiv:0903.0353]
 並列化はどうする？領域分割前処理+Deflation/Multigrid ?
32
5.クォーク伝播関数ソルバー
領域分割前処理 (Domain–decomposition)再び
 問題の微分方程式離散化
 問題空間をいくつかに分割（重複も可）
 流れ（基本的に反復改良法）
 １．初期解と、初期残差
 ２．残差に対して問題の式を分割された空
間で解く
 ３．解いた結果を使って近似解を構成
 ４．残差を計算＝＞２．にもどる
 解の更新の順番とか、境界の扱い方とか、
部分空間に制限する方法とか、重複部分
をどう扱うかとか、、、、
 さまざまなバリエーション
 この方法だけでは収束しないのでこの反
復をKrylov部分空間法の前処理として採
用する。
x0 , r  b  Ax0,
Solve : Ai xi  ri , in i
Update: x  x  f(xi ),
Compute: r  b  Ax,
Ω1
Γ2
Γ１
Ω2
33
5.クォーク伝播関数ソルバー
Luescher の導入した領域分割前処理：
 Schwartz alternating method (SAP): Two no-overlapping domain =
block Schur complement (Luscher) = Multiplicative Schwarz
Method
C.f. SAP=Multiplicative Schwarz
even blocks 1 
v    Dii   r
 j

x  x  v, r  r  Dv
Solve in Even domain
odd blocks 1 
v    Dii   r
 j

x  x  v, r  r  Dv

Solve in Odd domain

: x  x  DEE   DOO   DOO  DOE DEE  r
1
1
1
1
 小領域ソルバーをアクセラレータで加速
 並列度、ノードに２色入れる必要がある。Overlap 無しのためノード
の担当する格子点数は少ない
34
5.クォーク伝播関数ソルバー
 Multiplicative Schwarz (MS) vs Additve Schwarz (AS)
 MS : generalized Block Gauss-Seidel
 AS: generalized Block Jacobi
(MS > AS, factor 2)
 Restricted (Overlapping) Additive Schwarz (RAS) method
[Cai & Sarkis, SIAM J.S.C.21(1999)]
Projection
on a fermion field
Overlapped region
Depth 2
Solve in i-th domain
D_i x_i =r_i
Residual vector field r
35
5.クォーク伝播関数ソルバー
 RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない：
 並列度が上がる。（１ノード１領域）
 小領域を大きくしてアクセラレータの効率を上げる。
 領域の重なりにより前処理性能を上げる。
Return only original
region
Solution vector field x
x = x + \sum_i x_i
Overlapped region
Depth 2
36
5.クォーク伝播関数ソルバー
 Restricted (Overlapping) Additive Schwarz (RAS) method
 Iterate
all blocks 
v    M i r
 j

x  x  v, r  r  Dv

 DR
2d t
i
Mi  R R
0
i
2d
i
Each blocks are independent.
solved in each block in parallel.
GPGPU!
 R
1
,
2d t
i
: Block solver.
Ri2 d : overlapping (depth  2d) block latticeprojectionoperator
C.f. SAP=Multiplicative Schwarz
 Overlapping improves
performance. But task increased.
 Prolongation is not overlapping
(Restricted). This also improves
the performance.
 RASのみでは収束が遅い、MG
やDeflation が必要
even blocks 1 
v    Dii   r
 j

x  x  v, r  r  Dv
Solve in Even domain
odd blocks 1 
v    Dii   r
 j

x  x  v, r  r  Dv

Solve in Odd domain

: x  x  DEE   DOO   DOO  DOE DEE  r
1
1
1
1
37
5.クォーク伝播関数ソルバー
 Test on small lattice (16^3x32), timing comparison.
SAP+Defl
RAS(d=1)+Defl
Fast
SAP with Deflation is the best.
RAS(d=1) approaches SAP w/o deflation.
RAS(d=2,3) reduce iteration count by 1/21/3. But the task in each node is rapidly
increasing by overlapping reagion.
AS は MS にはやはり勝てない？
GPUで比較する予定
PC Cluster: 16 nodes.
Block size:
SAP: 8^4
RAS: (8+2d)^3x(16+2d)
Deflate10 small eigenvalues.
Best case comparison.
計算時間の9割は単精度計算
Slow
38
6. 固有値問題
 Ｄ［Ｕ］の固有値固有ベクトルが分かるといろいろ便利なことがある。
 Low-mode averaging
 物理量にはクォークD[U]の小さい固有値固有ベクトルが重要。
 並進対称性を利用した統計の向上
 Determinant splitting
 行列式の評価で、小さい固有値は正確に計算し残りはモンテカルロで評価する
 Epsilon regime
 D[U]のゼロ付近の固有値分布とカイラル対称性の破れの関係
 Deflation for quark solver
 大量のDx=b を解く時に Dの固有空間が分かるとそれを使って前処理して、高速
化
 Matrix function
 Overlap 型クォークの sign関数、奇数フレーバーのための sqrt(D), 行列式
をTrLog[D]で評価、前処理のための (D†D)(1/n) などなど、、、
39
6. 固有値問題
 固有値固有ベクトルを求めるのは Dx=b を一回だけ解くの
に比べてはるかに困難
 Lattice QCD で使われている手法
 H =γ5D : Hermite 行列 => Lanczos アルゴリズムとその系統
 D : non-Hermitian 行列 => Arnoldi アルゴリズムとその系統
 Lanczos系 Thick-restareted Lanczos
[Wu&Simon]
 PACC-CS collab. で使用中、他のグループでも使われているかも
 Dに対してはほぼ ARPACK (Implicitly Restarted Arnoldi Method)
が使われている。
[Lehoucq, Sorensen, Yang, Maschhoff]
 大体のLQCD論文で ARPACKが言及されている
 そのほかのアルゴリズムは？
 Jacobi-Davidoson 法は？
40
6. 固有値問題
 Arnoldi/Lanczos系での困難
 原点に近い固有値固有ベクトルを求めたいとき、逆反復にする(+shiftも
するかも)
 Arnoldi(Lanczos)反復で問題を小さな問題に変換
 1. 適当な初期ベクトルから直交基底ベクトルを構成していく（m反復後）
D 1Vm  Vm H m  hm 1,m vm 1em ,
T
Vm  (v1 , v2 ,....,vm ) : 正規直交基底( N  m)
H m : m  m 小行列
 2. Hm の固有値固有ベクトルから 1/D の大きいほうの固有値固有ベクトルを抽出
＝Dの小さい固有値固有ベクトルを抽出
 3. ほしいベクトルを何本か残して式の形が変化しないようにユニタリー変換して縮
小してリスタート(Implicit restart)
 4. ほしいベクトルが収束するまでリスタートを繰り返す




1/Dの計算は反復法で行う
計算コスト: 大体 (m x [リスタート回数]) 回、反復法で Dx=b を解く
得られる固有ベクトルの精度は 1/D の計算精度で決まってしまう。
Dx=b に対する反復法も良い精度で解く必要がある。
41
6. 固有値問題
 Arnoldi/Lanczos系での困難
 固有ベクトルの情報を何かの計算の加速に使おうと思うと、
 連立一次方程式 Dx=b を精度よく何回も解かないといけないので本末転
倒しているような気がする。
 私がいろいろ調べた結果
 Jacobi-Davidson法はその困難がない？
 Deflation付のソルバーを使わないといけない様である。
 固有値固有ベクトルペア(Ritz pair)を１個づつ求めていく様である。
 ターゲットの Ritz pair の選び方とか、うまくリスタートする実装の仕方がまだ良く理
解できていないのでテストしていません。
 RESidual Arnoldi Method (RESAM)




[C.R.Lee & G.W.Stewart, Tech. Report TR-2007-45;
C.R.Lee, Ph.D. Thesis]
部分空間反復法とArnoldi 法の合成みたいな方法
Jacobi-Davidson法にも似ているが Deflation 付ソルバがいらない
D論なのでくわしく書いてある
売りは、Dx=b は低精度で解くので十分なので IRAMよりも速い
IRAM(KSRAM) と RESAM を QCD で比較してみました。
42
6. 固有値問題
 RESRAM
 1. 初期部分空間を Arnoldi法で作っておく
 2. その部分空間から Ritz pair を用いてターゲットとなる近似固有値固有ベク
トルを求める
 3. その近似固有値固有ベクトルの残差を計算する。収束していたら
Deflation に加える、次のターゲットを探すべく 2. へ戻る
 4. 残差に対して 1/D の近似をかける(Ritz値へshift しても良い)
 5. それを以前の部分空間に対して正規直交化して部分空間に付け加えて部
分空間を拡張する
 6.部分空間の次元がある最大値に達したら Schur型を経由して次元を縮め
る。その際収束している部分空間が変化しないようにうまくする。
 Step. 4 は精度を求めないので単精度でも OK
 最終結果に倍精度を求めるなら、それ以外の計算、特に固有方程式の
残差の計算は倍精度で行う必要がある。
43
6. 固有値問題
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 IRAM
 ARPACK をそのままコンパイル
 KSRAM
 IRAMを自分で実装できなかったのでIRAMとほぼ同等な Krylov-Schur Restarted Arnoldi Method
[TRLANのnon-Hermitan版] を実装,
 RESAM
原点付近の固有値を２０個求める
 以下を変化させて速度測定
 変化させるパラメータ：
 線形ソルバー(Dx=b)の精度、
 求める固有値固有ベクトルの精度、
 線形ソルバーは共通の GCRO-DR法
 最大部分空間次元は 40
 リスタート時に残す部分空間次元は
 min((すでに求まった部分空間次元＋20),40)
 環境：PC一台 : CPU Core2Duo [email protected]
44
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 欲しい固有値固有ベクトルの
精度よりも線形ソルバーの精
度を落としたときの固有値固
有ベクトルの残差はどうなる
か？
 IRAM:線形ソルバーの精度
より良くなることはない
 RESAM: 常に必要な精度が
得られる
45
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 倍精度で解を求める場合の計算
時間の比較
 IRAM,KSRAM:線形ソルバーは
倍精度で説き続ける
 RESAM:線形ソルバーは1/10の
精度で求めればよい。
 全体の計算時間のほとんどは線
形ソルバーの時間に費やされる。
 固有値ソルバーの反復回数は
 IRAM/KSRAM << RESAM
 積算された時間は
 RESRAM << IRAM/KSRAM
46
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 RESAMで倍精度で解を求める場
合の固有ペア残差のヒストリ
 １本ずつ原点に近いものから求め
ていっている。
 線形ソルバーの停止条件は残差が
1/10 を切ったとき。
 線形ソルバーの停止条件をきつく
すると、REAAMの反復回数は減っ
ていくが、線形ソルバーの時間の増
加を補うほどではないので帰って全
体が遅くなる。
 Best Time は tol_solver = 1/10の
時
47
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 加速率の比較(Best同士の比較)
 IRAM,KSRAM:
 (求める固有ペアの精度)
= (線形ソルバーの停止精度)/2
 RESAM:
 線形ソルバーの停止精度=1/10
 高精度(誤差10^-4以下)の固有ペ
アを求めるときRESRAM は
IRAM/KSRAM より速い。
 倍精度ぎりぎりのときは特に倍速
い
 本格的に大きな格子ではどう
か？（並列化必要）
48
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より]
 加速率の比較(Best同士の比較)
 IRAM,KSRAM:
 (求める固有ペアの精度)
= (線形ソルバーの停止精度)/2
 RESAM:
 線形ソルバーの停止精度=1/10
 高精度(誤差10^-4以下)の固有ペ
アを求めるときRESRAM は
IRAM/KSRAM より速い。
 倍精度ぎりぎりのときは特に倍速
い
 本格的に大きな格子ではどう
か？（並列化必要）
49
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト中問題 [16^3x32 格子]
 線形ソルバー
 KSRAM:
 Nested BiCGStab(単精度部
分: BiCGStab/GCRO-DR)
 RESRAM:
 単精度 BiCGStab/GCRO-DR
 計算時間の比較(Best同士の比
較)
 KSRAM:
 (求める固有ペアの精度)
= (線形ソルバーの停止精度)/2
 RESAM:
 線形ソルバーの停止精度=0.5
 すべての必要精度で、RESRAM
が速いようだ。
50
6. 固有値問題
原点付近の固有値を２０個求める
 比較テスト中問題 [16^3x32 格子]
 加速率の比較(Best同士の比較)
 KSRAM:
 (求める固有ペアの精度)
= (線形ソルバーの停止精度)/2
 RESAM:
 線形ソルバーの停止精度=0.5
 すべての必要精度で、RESRAMが
速いようだ。大体倍ぐらい速い。
 QCD(Wilson-Dirac型)では確かに
RESRAMは速そうだ。
 Hermit 版はどうか？
 Jacobi-Davidson法も同様に速い
のかも？
線形ソルバーを低精度に出来る
＝＞アクセラレータでさらに加速？
51
7．まとめ
 格子QCD: クォーク行列の逆、行列式などの評価の高速化
が重要
 Wilson-Dirac型クォーク：
 係数行列の形は簡単＝プログラム上の最適化では限界がすぐ来る
＝＞アルゴリズムの改良は重要
 さまざまな改良：
 前処理：領域分割、Deflation、MultiGrid、
 精度を落として加速？並列度を上げる？アクセラレータ？
 固有ペアの計算
 これらの改良の組み合わせ
 紹介し切れなかったこと
 行列関数(Overlap型クォークなど)の評価法
 分子動力学の改良
 マルコフチェーンモンテカルロ法の改良？
52

Download Report