格子QCDのための線形計算 (連立一次方程式、固有値問題)について

Slides:



Advertisements
Similar presentations
1 広島大学 理学研究科 尾崎 裕介 石川 健一. 1. Graphic Processing Unit (GPU) とは? 2. Nvidia CUDA programming model 3. GPU の高速化 4. QCD with CUDA 5. 結果 6. まとめ 2.
Advertisements

Absolute Orientation. Absolute Orientation の問題 二つの座標系の間における剛体 (rigid body) 変換を復元す る問題である。 例えば: 2 台のステレオカメラから得られた3次元情報の間の関 係を推定する問題。 2 台のステレオカメラから得られた3次元情報の間の関.
大規模な三角 Toeplitz 線形方 程式の高速解法とその応用 ○ 安村 修一(法政大学 4 年) 李 磊(法政大学) 日本応用数理学会「行列・固有値の解法とその応用」研究部会 第6回研究会.
QCD Sum rule による中性子電気双極子 モーメントの再評価 永田 夏海(名古屋大学) 2012 年 3 月 27 日 日本物理学会第 67 回年次大会 共同研究者:久野純治,李廷容,清水康弘 関西学院大学.
第 5 章 2 次元モデル Chapter 5 2-dimensional model. Contents 1.2 次元モデル 2-dimensional model 2. 弱形式 Weak form 3.FEM 近似 FEM approximation 4. まとめ Summary.
1. 補間多項式 n 次の多項式 とは、単項式 の 線形結合の事である。 Definitions: ( 区間, 連続関数, abscissas (データ点、格子点、差分点), 多項 式 ) Theorem. (補間多項式の存在と一意性) 各 i = 0, …, n について、 をみたす、次数が高々 n.
Determining Optical Flow. はじめに オプティカルフローとは画像内の明る さのパターンの動きの見かけの速さの 分布 オプティカルフローは物体の動きの よって変化するため、オプティカルフ ローより速度に関する情報を得ること ができる.
有限差分法による 時間発展問題の解法の基礎
CPUとGPUの 性能比較 -行列計算およびN体問題を用いて-
計算理工学基礎 「ハイパフォーマンスコンピューティングの基礎」
Fill-in LevelつきIC分解による 前処理について
GPU上の大規模格子QCDシミュレーション に向けて
一般化Bi-CGSTAB(s, L) (=一般化IDR(s, L))
A Q R QR分解とは? → × ◆QR分解 QTQ = I (単位行列) ◆応用例 ◆主な計算方法 n m 今回はこの方法に注目
格子QCDシミュレーションにおける固有値問題
Finger patternのブロック化による 陰的wavelet近似逆行列前処理の 高速化
剛体の物理シミュレーション は難しい? 佐藤研助手 長谷川晶一.
スペクトル法による数値計算の原理 -一次元線形・非線形移流問題の場合-
AllReduce アルゴリズムによる QR 分解の精度について
重力3体問題の数値積分Integration of 3-body encounter.
「データ学習アルゴリズム」 第3章 複雑な学習モデル 3.1 関数近似モデル ….. … 3層パーセプトロン
P,Q比が変更可能なScaLAPACKの コスト見積もり関数の開発
PCクラスタ上での 連立一次方程式の解の精度保証
BRST対称性とカイラル対称性の破れの 格子QCDシミュレーションによる検証 scirqcd 帝京大学理工学部 古井貞隆
理学部情報科学科 金田研究室 指導教官 金田 康正 工藤 誠
最短路問題のための LMS(Levelwise Mesh Sparsification)
スパコンとJLDG HEPの計算環境 HEPnet-J
ML 演習 第 7 回 新井淳也、中村宇佑、前田俊行 2011/05/31.
シミュレーション演習 G. 総合演習 (Mathematica演習) システム創成情報工学科
正方行列向け特異値分解の CUDAによる高速化
スペクトル・時系列データの前処理方法 ~平滑化 (スムージング) と微分~
文献名 “Performance Tuning of a CFD Code on the Earth Simulator”
高次元データの解析 -平均ベクトルに関する検定統計量の 漸近分布に対する共分散構造の影響-
現実の有限密度QCDの定性的な振る舞いに
領域分割手法について 2008年2月26日 中島研吾.
非エルミート 量子力学と局在現象 羽田野 直道 D.R. Nelson (Harvard)
Graphic Card を使った 高性能計算
格子QCDにおけるGPU計算 広大理 尾崎裕介 共同研究者 石川健一.
格子QCD計算のGPUを用いた加速の可能性について
スペクトル法の一部の基礎の初歩への はじめの一歩
人工知能特論 9.パーセプトロン 北陸先端科学技術大学院大学 鶴岡 慶雅.
Jh NAHI 横田 理央 (東京工業大学) Hierarchical low-rank approximation methods on distributed memory and GPUs 背景  H行列、H2行列、HSS行列などの階層的低ランク近似法はO(N2)の要素を持つ密行列をO(N)の要素を持つ行列に圧縮することができる。圧縮された行列を用いることで、行列積、LU分解、固有値計算をO(NlogN)で行うことができるため、従来密行列の解法が用いられてきた分野では階層的低ランク近似法
6. ラプラス変換.
変換されても変換されない頑固ベクトル どうしたら頑固になれるか 頑固なベクトルは何に使える?
連続領域におけるファジィ制約充足問題の 反復改善アルゴリズムによる解法 Solving by heuristic repair Algorithm of the Fuzzy Constraint Satisfaction Problems with Continuous Domains 北海道大学.
階層的境界ボリュームを用いた 陰関数曲面の高速なレイトレーシング法
パターン認識特論 担当:和田 俊和 部屋 A513 主成分分析
GPUを用いた疎行列の格納形式による行列ベクトル積の評価
pp-wave上の共変的超弦の場 における低エネルギー作用
知識科学研究科 知識システム構築論講座 林研究室 佛明 智
円柱座標系の基底関数系を用いたSCF法による 円盤銀河のシミュレーション
ベイジアンネットワーク概説 Loopy Belief Propagation 茨城大学工学部 佐々木稔
わかりやすいパターン認識 第7章:部分空間法  7.1 部分空間法の基本  7.2 CLAFIC法                  6月13日(金)                  大城 亜里沙.
曲がった時空上の場の理論の熱的な性質と二次元CFT
? 格子QCDシミュレーションによる南部-ゴールドストン粒子の 質量生成機構の研究 質量の起源 ドメインウォールフェルミオン作用
短い部分文字列の ミスマッチトレランスを 高速計算するアルゴリズム
格子ゲージ理論によるダークマターの研究 ダークマター(DM)とは ダークマターの正体を探れ!
計画研究代表者 京都大学基礎物理学研究所 大野木 哲也
ポッツスピン型隠れ変数による画像領域分割
原子核物理学 第7講 殻模型.
ガウス分布における ベーテ近似の理論解析 東京工業大学総合理工学研究科 知能システム科学専攻 渡辺研究室    西山 悠, 渡辺澄夫.
格子QCDの理論的進展と フレーバー物理への応用II
Jh NAHI 横田 理央 (東京工業大学) Hierarchical low-rank approximation methods on distributed memory and GPUs 背景  H行列、H2行列、HSS行列などの階層的低ランク近似法はO(N2)の要素を持つ密行列をO(N)の要素を持つ行列に圧縮することができる。圧縮された行列を用いることで、行列積、LU分解、固有値計算をO(Nlog2N)で行うことができるため、従来密行列の解法が用いられてきた分野では階層的低ランク近似
メモリ使用量の少ないGCR法の提案 東京大学理学部情報科学科 工藤 誠 東京大学情報基盤センター 黒田 久泰
密行列固有値解法の最近の発展 (I) - Multiple Relatively Robust Representation アルゴリズム - 2004年11月26日 名古屋大学 計算理工学専攻 山本有作 日立製作所の山本有作です。 「~」について発表いたします。
シミュレーション演習 G. 総合演習 (Mathematica演習) システム創成情報工学科
媒質中でのカイラル摂動論を用いた カイラル凝縮の解析
目次 はじめに 収束性理論解析 数値実験 まとめ 特異値計算のための dqds 法 シフトによる収束の加速
? リー・ヤンの零点 これまでの格子QCD計算の結果 今年度の計画 リー・ヤンの零点分布から探る有限密度QCDにおける相構造の研究
8.数値微分・積分・微分方程式 工学的問題においては 解析的に微分値や積分値を求めたり, 微分方程式を解くことが難しいケースも多い。
Presentation transcript:

格子QCDのための線形計算 (連立一次方程式、固有値問題)について 石川健一(広島大理) アルゴリズムによる計算科学の融合と発展 @CCS, 筑波大学, 2009年4月22日・23日

1. 目次 2.格子QCDについて 3.ハイブリッドモンテカルロ法 4.格子化クォークのタイプ 5.クォーク伝播関数ソルバー 6.固有値問題 7.まとめ

時空を格子化した有限自由度の格子QCDという方法を用いる。 2. 格子QCDについて QCDを解くためには 時空を格子化した有限自由度の格子QCDという方法を用いる。 K.G.Wilson (1974) 連続時空上の場の変数 → 格子上の場の変数 格子化した時空 連続時空

2. 格子QCDについて 分配関数(ユークリッド化された経路積分) 物理量の期待値 統計力学で用いられてきた方法 作用 有効作用 有効作用を重みとする多重積分 => モンテカルロ積分 統計力学で用いられてきた方法

2. 格子QCDについて 特徴、 自由度がとても大きい 16^4 ~ 32^4 格子 系統誤差 自由度がとても大きい  16^4 ~ 32^4 格子 グルーオン:8x4x(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 典型的な計算時間1~2年(論文、予算、、、、) どのようにしてモンテカルロ積分のための配位を生成するか?

3. ハイブリッドモンテカルロ(HMC)法 格子QCD分配関数(Nf=2:クォークが2種類) ユークリッド化されているので、統計力学の分配関数と同等の計算になる。

モンテカルロ重み=行列式=補助場による分布 D[U]の逆を含む:非局所的 3.ハイブリッドモンテカルロ(HMC)法 クォーク部分の計算が大変 モンテカルロ重み=行列式=補助場による分布 D[U]の逆を含む:非局所的 φで積分すると 真空からクォーク-反クオーク対が対生成、対消滅している効果を表している。真空偏極。

Hybrid Monte Carlo(HMC)法 Exp(-H) の分布に従うU,P,φ を HMC法で生成 Uのサンプル{ U(1),U(2),U(3),…..,U(N)} 物理量の期待値(物理量は U のみの関数 O[U])

HMC法 3.ハイブリッドモンテカルロ(HMC)法 に対して、マルコフチェーンモンテカルロ法の一種HMC法を適用    に対して、マルコフチェーンモンテカルロ法の一種HMC法を適用 分子動力学(Molecular Dynamics) (MD) メトロポリス法(Metropolis) 分子動力学ではエネルギーは差分化のため保存しない 分布がほしい分布        からずれる で補正する 分子動力学法についての発展は今回は省略

Hybrid Monte Carlo(HMC)法 計算が困難な点 分子動力学法の部分で の大規模連立1次方程式を何回も解く必要がある。 幸いD[U]の要素はほとんどゼロ(大規模疎行列)であるので、反復法 を用いて連立方程式を解く。 しかし特にクォークの質量が小さいと、D[U]がゼロに近い固有値を持つため連立方程式を解くのが難しくなる。 クォークの質量が小さいと分子動力学の力が大きくなる。分子動力学の時間刻みを小さくする必要がある。

問題をまとめると、 3.ハイブリッドモンテカルロ(HMC)法 クォークのタイプによって D[U] がいろいろある モンテカルロ法による経路積分の評価で困難が生じるのはクォークの寄与を取り入れる部分。 Uの生成には、 D[U] x = b の連立一次方程式を大量に解く必要がある これまでさまざまな前処理方法が開発されてきた 最近の発展について述べたい。 クォークソルバー

4.格子化クォークのタイプ 離散化には任意性がある 良い性質を持つものほど計算コストがかかる クォークの持つ対称性 ゲージ対称性 :  格子上でも死守すべき対称性 ポアンカレ対称性: 一部破れるけどなるべく保持したい対称性、ユークリッド化後:離散並進、離散回転、反転対称性 カイラル対称性:            格子上で厳密に実現することは不可能(ニールセン-二宮の定理) カイラル対称性の保持の仕方による分類 Kogut-Susskind 型:4つの同じ質量をもつクォーク(4フレーバー)を同時に扱う。U(1)xU(1)カイラル対称性を保持 Wilson-Dirac型:1つのクォークから扱える。カイラル対称性は無い Overlap型、Domain-Wall 型:1つのクォークから扱える。カイラル対称性は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて10-100倍)

4.格子化クォークのタイプ Kogut-Susskind(KS) 型:4つの同じ質量をもつクォーク(4フレーバー)を同時に扱う。U(1)xU(1)カイラル対称性を保持 最も簡単なD[U]の構造を持つ。現実世界のクォークは6種類で質量もばらばら。軽いのはアップとダウンの2種類のクォーク = 4つ同時では扱えない。スキップします。

4.格子化クォークのタイプ Wilson-Dirac型:1つのクォークから扱える。カイラル対称性は無い。 D[U]は一階の差分演算子 クォーク質量 m ∝ 1/κ カラー(3),スピン(4)の自由度が各格子点にある ブロックサイズ12のブロックの帯が並んだ疎行列 M の固有値分布 U=1,4^4 Witzel,Takeda,Wolff, hep-lat[arXiv:0709.4648]

4.格子化クォークのタイプ Overlap型、Domain-Wall 型:1つのクォークから扱える。カイラル対称性は変形された対称性を保持。計算コストがとても高い(Wilson型に比べて10-100倍) 質量 m Wilson-Diracクォークの行列 DW のサイン関数を含む Dov は密な行列 対称性: 質量m=0 でGinparg-Wilson関係を満たす はユニタリー行列 Taken from Talk slid by R.Edwards

4.格子化クォークのタイプ まとめ、 Wilson-Dirac型 Dx = b は反復法で解く 係数行列 D の型 1階差分型:さまざまな前処理が適用されてきた Red-brack(2色), ILU,SSOR, 2色ブロック … etc. Dx = b を直接解く: Non-Hermitian sovler, BiCGStab, GMRES などが使われてきた (D†D)x=D†b を解く: positive Hermitian solver,           CG が使われる。条件数悪し。 γ5Dx=γ5D を解く: indefinite Hermitian solver, MINRESなど、、 前処理つき BiCGStab, GMRES 系統が良く使われる

4.格子化クォークのタイプ Overlap 型 Dをベクトルにかけるときに sign 関数の計算が必要 Dx=b 自体は Hermite化後 CG などで解く Sign関数の計算 Arnoldi / Lanczos 分解による方法   Pade / Rational :部分分数展開近似による方法 部分分数/連分数展開+補助場の導入による方法                          => 5次元化 ソルバーの入れ子を解消 HW や DWの固有値範囲の情報を用い sign関数の近似の精度をコントロールする必要がある。 Overlap 型についてはこれ以上次官の都合上言及しない 日本では JLQCD collaboration でくわしく研究されている KEK 松古、金児、橋本…..

5.クォーク伝播関数ソルバー Dx=bを速く解きたい 格子サイズ: 32^3x64 (現状 O(10)TFlops)          => 256倍: 128^3x256 (O(2)PFlops) 次元: 12x[格子サイズ]= 2500万 => 64億 (1)前処理: Wilson型に関しては構造が簡単(Red-Black, ILD, Domain-Decomposition etc.), Overlap型については難しい⇒5D化で構造が簡単に [省略] (2)固有値固有ベクトルの情報:  Deflation/Multigrid (3)単精度加速: (4)アクセラレータ:

(2) Deflation technique 5.クォーク伝播関数ソルバー (2) Deflation technique LQCD では何回も Dx=b を解く 右辺ベクトルを取り替えながら、係数行列は一定。または、係数行列が少しづつ変化しながら連鎖的に解く必要がある Quark propagator Solver in HMC trajectory ブロック化ソルバー:係数行列固定、右辺ベクトル多数 Deflation technique: 係数行列の固有ベクトルの情報により前処理 Deflation remove/suppress small eigenspace of D.

(2) Deflation technique (Solve Ax=b ・・・(1) ) 5.クォーク伝播関数ソルバー (2) Deflation technique (Solve Ax=b ・・・(1) ) 行列A : 小さい固有値に属する p-次元部分空間を持つとする c,u : この部分空間に属する以下の条件を満たすベクトルとする c から射影演算子を定義 Pを使って問題を前処理 式(1)の解 x は式(2) の解 y を用いて以下のように書かれる “Cp”はどうやって作る? 反復には P の掛け算が伴う(overhead)

5.クォーク伝播関数ソルバー (2) Deflation technique (cont’d) Many works by 正確な不変部分空間の計算はとてもコストが高い     低モード密度 ∝ 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) [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.

(a)固有空間の計算と連立一次方程式を同時に解く. 5.クォーク伝播関数ソルバー (a)固有空間の計算と連立一次方程式を同時に解く. Very effective for few Near zero modes / negative eigen modes case. Near zero modes case 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] [Wilcox, LAT2007]

(b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. 5.クォーク伝播関数ソルバー (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. 基底ベクトルの構成方法が重要 Low mode ベクトルはブロックベクトルで近似できる [Luscher, JHEP07(2007)081]

(b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. 5.クォーク伝播関数ソルバー [Luscher, JHEP07(2007)081] (b)ブロック化された基底ベクトルで元の行列の逆を取り近似固有空間を構成する. B : ブロック化された部分空間でのD[U]の表現(=サブグリッド上に射影された) 基底ベクトル(C)の構成 ランダムベクトルに1/D を数回かけて低モードを増幅 そのベクトルをサブブロックに分割 分割されたベクトルを正規直交化 Deflation の射影演算子 (B^-1)を含む.

(2) MultiGrid Solver 5.クォーク伝播関数ソルバー MultiGrid solver also removes critical slowing down. Choice of subspace basis is important. (Prolongator) 射影のための基底ベクトルの選び方が重要(Luscher の基底と同じ) QCD, 16^3x32 [Brannick,Brower,Clark,Osborn,Rebbi, PRL100(2008)]

(3) 単精度加速 単精度の利点: 有効バンド幅、キャッシュサイズ、(レジスターの数)などが倍になる=効率(実性能/理論性能)が高い 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反復のテクニックで可能 一般の反復法にも組み入れることが出来る(可変前処理) 前処理として単精度ソルバーを使う 残差は倍精度で正しくなるように組む FGMRES, GCR, CG, BiCGStab….. (with flexible prec.) [Buttari,Dongarra,Kurzak,Luszczek,Tomov, AMC Trans.Math.Soft]

Nested-BiCGStab (櫻井-多田野) 5.クォーク伝播関数ソルバー (3) 単精度加速 Nested-BiCGStab (櫻井-多田野) Outer solver : BiCGStab (D.P.) Inner solver: BiCGStab (S.P.)= Outer solverに対する可変前処理 Intel 64, SSE3 使用+リンク再構築で全て倍精度の時より速度が倍に(演算量は増えているが) [PACS-CS collab.]

5.クォーク伝播関数ソルバー (4) アクセラレーター 単精度加速と組み合わせるとさらに有効なもの 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)

5.クォーク伝播関数ソルバー (4) アクセラレーター 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倍以上速くなるのは大きい

5.クォーク伝播関数ソルバー (4) アクセラレーター 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 (ロード)と整数への丸め(ストア)

5.クォーク伝播関数ソルバー (4) アクセラレーター NVIDIA CUDAを使った例 [M.Clark et al. Work shop@CCS, 10-12 March, 2009] 半精度(16bit int)でさらに2割性能上昇 そのほかの工夫でさらに加速 ゲージ固定 SU(3)行列を実8パラメータで保存 GPGPU:演算が圧倒的に速いので、バンド幅、キャッシュの有効利用のほうが重要

[H.Baier at al. ,Lat2008, arXive:0801.1559 5.クォーク伝播関数ソルバー [Spary,Hill,Trew hep-lat/0804.3654; S.Motoki & A. Nakamura Lat2007; F.Belletti et al. LAT2007] (4) アクセラレーター Cell B.E.での計算 大規模 Cell B.E.クラスタ: QPACE (QCD Parallel computing on Cell B.E.) project (EU) 2048 PowerXCell8i, + カスタムネットワーク3Dトーラス(FPGA) 2009春完成? GPUクラスター: 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 ? [H.Baier at al. ,Lat2008, arXive:0801.1559

領域分割前処理 (Domain–decomposition)再び 5.クォーク伝播関数ソルバー 領域分割前処理 (Domain–decomposition)再び 問題の微分方程式離散化 問題空間をいくつかに分割(重複も可) 流れ(基本的に反復改良法) 1.初期解と、初期残差 2.残差に対して問題の式を分割された空間で解く 3.解いた結果を使って近似解を構成 4.残差を計算=>2.にもどる 解の更新の順番とか、境界の扱い方とか、部分空間に制限する方法とか、重複部分をどう扱うかとか、、、、 さまざまなバリエーション この方法だけでは収束しないのでこの反復をKrylov部分空間法の前処理として採用する。 Ω1 Ω2 Γ2 Γ1

Luescher の導入した領域分割前処理 : 5.クォーク伝播関数ソルバー Luescher の導入した領域分割前処理 : Schwartz alternating method (SAP): Two no-overlapping domain = block Schur complement (Luscher) = Multiplicative Schwarz Method 小領域ソルバーをアクセラレータで加速 並列度、ノードに2色入れる必要がある。Overlap 無しのためノードの担当する格子点数は少ない C.f. SAP=Multiplicative Schwarz Solve in Even domain Solve in Odd domain

Multiplicative Schwarz (MS) vs Additve Schwarz (AS) 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

RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない: 5.クォーク伝播関数ソルバー RAS: 領域を重ね、領域間の依存性を無視、重なり領域は戻さない: 並列度が上がる。(1ノード1領域) 小領域を大きくしてアクセラレータの効率を上げる。 領域の重なりにより前処理性能を上げる。 Return only original region Overlapped region Depth 2 Solution vector field x x = x + \sum_i x_i

Restricted (Overlapping) Additive Schwarz (RAS) method 5.クォーク伝播関数ソルバー Restricted (Overlapping) Additive Schwarz (RAS) method Iterate Overlapping improves performance. But task increased. Prolongation is not overlapping (Restricted). This also improves the performance. RASのみでは収束が遅い、MG やDeflation が必要 Each blocks are independent. solved in each block in parallel. GPGPU! C.f. SAP=Multiplicative Schwarz Solve in Even domain Solve in Odd domain

Test on small lattice (16^3x32), timing comparison. 5.クォーク伝播関数ソルバー Test on small lattice (16^3x32), timing comparison. PC Cluster: 16 nodes. Block size: SAP: 8^4 RAS: (8+2d)^3x(16+2d) Deflate10 small eigenvalues. Best case comparison. 計算時間の9割は単精度計算 SAP+Defl RAS(d=1)+Defl Fast Slow SAP with Deflation is the best. RAS(d=1) approaches SAP w/o deflation. RAS(d=2,3) reduce iteration count by 1/2-1/3. But the task in each node is rapidly increasing by overlapping reagion. AS は MS にはやはり勝てない? GPUで比較する予定

6. 固有値問題 D[U]の固有値固有ベクトルが分かるといろいろ便利なことがある。 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) などなど、、、

固有値固有ベクトルを求めるのは Dx=b を一回だけ解くのに比べてはるかに困難 Lattice QCD で使われている手法 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)が使われている。 大体のLQCD論文で ARPACKが言及されている そのほかのアルゴリズムは? Jacobi-Davidoson 法は? [Lehoucq, Sorensen, Yang, Maschhoff]

6. 固有値問題 Arnoldi/Lanczos系での困難 原点に近い固有値固有ベクトルを求めたいとき、逆反復にする(+shiftもするかも) Arnoldi(Lanczos)反復で問題を小さな問題に変換 1. 適当な初期ベクトルから直交基底ベクトルを構成していく(m反復後) 2. Hm の固有値固有ベクトルから 1/D の大きいほうの固有値固有ベクトルを抽出=Dの小さい固有値固有ベクトルを抽出 3. ほしいベクトルを何本か残して式の形が変化しないようにユニタリー変換して縮小してリスタート(Implicit restart) 4. ほしいベクトルが収束するまでリスタートを繰り返す 1/Dの計算は反復法で行う 計算コスト: 大体 (m x [リスタート回数]) 回、反復法で Dx=b を解く 得られる固有ベクトルの精度は 1/D の計算精度で決まってしまう。 Dx=b に対する反復法も良い精度で解く必要がある。 

6. 固有値問題 Arnoldi/Lanczos系での困難 固有ベクトルの情報を何かの計算の加速に使おうと思うと、 連立一次方程式 Dx=b を精度よく何回も解かないといけないので本末転倒しているような気がする。 私がいろいろ調べた結果 Jacobi-Davidson法はその困難がない? Deflation付のソルバーを使わないといけない様である。 固有値固有ベクトルペア(Ritz pair)を1個づつ求めていく様である。 ターゲットの Ritz pair の選び方とか、うまくリスタートする実装の仕方がまだ良く理解できていないのでテストしていません。 RESidual Arnoldi Method (RESAM) 部分空間反復法とArnoldi 法の合成みたいな方法 Jacobi-Davidson法にも似ているが Deflation 付ソルバがいらない D論なのでくわしく書いてある 売りは、Dx=b は低精度で解くので十分なので IRAMよりも速い [C.R.Lee & G.W.Stewart, Tech. Report TR-2007-45; C.R.Lee, Ph.D. Thesis] IRAM(KSRAM) と RESAM を QCD で比較してみました。

6. 固有値問題 RESRAM Step. 4 は精度を求めないので単精度でも OK 1. 初期部分空間を Arnoldi法で作っておく 2. その部分空間から Ritz pair を用いてターゲットとなる近似固有値固有ベクトルを求める 3. その近似固有値固有ベクトルの残差を計算する。収束していたら Deflation に加える、次のターゲットを探すべく 2. へ戻る 4. 残差に対して 1/D の近似をかける(Ritz値へshift しても良い) 5. それを以前の部分空間に対して正規直交化して部分空間に付け加えて部分空間を拡張する 6.部分空間の次元がある最大値に達したら Schur型を経由して次元を縮める。その際収束している部分空間が変化しないようにうまくする。 Step. 4 は精度を求めないので単精度でも OK 最終結果に倍精度を求めるなら、それ以外の計算、特に固有方程式の残差の計算は倍精度で行う必要がある。

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 6700@2.66GHz 原点付近の固有値を20個求める

6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 原点付近の固有値を20個求める 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 欲しい固有値固有ベクトルの精度よりも線形ソルバーの精度を落としたときの固有値固有ベクトルの残差はどうなるか? IRAM:線形ソルバーの精度より良くなることはない RESAM: 常に必要な精度が得られる

6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 原点付近の固有値を20個求める 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 倍精度で解を求める場合の計算時間の比較 IRAM,KSRAM:線形ソルバーは倍精度で説き続ける RESAM:線形ソルバーは1/10の精度で求めればよい。 全体の計算時間のほとんどは線形ソルバーの時間に費やされる。 固有値ソルバーの反復回数は IRAM/KSRAM << RESAM 積算された時間は RESRAM << IRAM/KSRAM

6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 原点付近の固有値を20個求める 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] RESAMで倍精度で解を求める場合の固有ペア残差のヒストリ 1本ずつ原点に近いものから求めていっている。 線形ソルバーの停止条件は残差が1/10 を切ったとき。 線形ソルバーの停止条件をきつくすると、REAAMの反復回数は減っていくが、線形ソルバーの時間の増加を補うほどではないので帰って全体が遅くなる。 Best Time は tol_solver = 1/10の時

6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 原点付近の固有値を20個求める 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 加速率の比較(Best同士の比較) IRAM,KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=1/10 高精度(誤差10^-4以下)の固有ペアを求めるときRESRAM は IRAM/KSRAM より速い。 倍精度ぎりぎりのときは特に倍速い 本格的に大きな格子ではどうか?(並列化必要)

6. 固有値問題 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 原点付近の固有値を20個求める 比較テスト 小問題 [4^4 格子 beta=5.0 クエンチ Matrix market より] 加速率の比較(Best同士の比較) IRAM,KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=1/10 高精度(誤差10^-4以下)の固有ペアを求めるときRESRAM は IRAM/KSRAM より速い。 倍精度ぎりぎりのときは特に倍速い 本格的に大きな格子ではどうか?(並列化必要)

6. 固有値問題 比較テスト 中問題 [16^3x32 格子] 原点付近の固有値を20個求める 線形ソルバー KSRAM: Nested BiCGStab(単精度部分: BiCGStab/GCRO-DR) RESRAM: 単精度 BiCGStab/GCRO-DR 計算時間の比較(Best同士の比較) (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=0.5 すべての必要精度で、RESRAMが速いようだ。

6. 固有値問題 比較テスト 中問題 [16^3x32 格子] 原点付近の固有値を20個求める 加速率の比較(Best同士の比較) KSRAM: (求める固有ペアの精度) = (線形ソルバーの停止精度)/2 RESAM: 線形ソルバーの停止精度=0.5 すべての必要精度で、RESRAMが速いようだ。大体倍ぐらい速い。 QCD(Wilson-Dirac型)では確かにRESRAMは速そうだ。 Hermit 版はどうか? Jacobi-Davidson法も同様に速いのかも? 線形ソルバーを低精度に出来る =>アクセラレータでさらに加速?

7.まとめ 格子QCD: クォーク行列の逆、行列式などの評価の高速化が重要 Wilson-Dirac型クォーク: さまざまな改良: 係数行列の形は簡単=プログラム上の最適化では限界がすぐ来る=>アルゴリズムの改良は重要 さまざまな改良: 前処理:領域分割、Deflation、MultiGrid、 精度を落として加速? 並列度を上げる? アクセラレータ? 固有ペアの計算 これらの改良の組み合わせ 紹介し切れなかったこと 行列関数(Overlap型クォークなど)の評価法 分子動力学の改良 マルコフチェーンモンテカルロ法の改良?