--- title: 第9回 配信講義　計算科学技術特論A （2023） tags: #計算科学 #高性能計算技術 #アルゴリズム author: [R-CCS　計算科学研究推進室](https://www.docswell.com/user/2300203199) site: [Docswell](https://www.docswell.com/) thumbnail: https://bcdn.docswell.com/page/LELZZ3V2JR.jpg?width=480 description: 第9回 6月15日 行列計算における高速アルゴリズム2 本講義では、科学技術計算の基盤となる行列計算について、基本的なアルゴリズムの原理と高性能計算技術を紹介する。 第2回ではポストペタ計算機に向けたアルゴリズムの最適化技法を紹介する。また、量子アニーリングマシンの行列計算への適用可能性についても簡単に触れる。 published: June 12, 23 canonical: https://www.docswell.com/s/2300203199/ZRXG71-2023-06-12-130947 --- # Page. 1 ![Page Image](https://bcdn.docswell.com/page/LELZZ3V2JR.jpg) 行列計算における高速アルゴリズム ― エクサフロップス時代における 線形計算アルゴリズムの課題と研究動向 ― 2023年6月15日 計算科学技術特論A 電気通信大学 情報理工学研究科 情報・ネットワーク工学専攻 山本 有作 # Page. 2 ![Page Image](https://bcdn.docswell.com/page/4JMVVDG9JW.jpg) 本講義の構成 • 第1回（6月8日） 「大規模連立1次方程式の反復解法」 – 代表的な行列計算である連立1次方程式の解法について，現在最も 広く利用されているクリロフ部分空間法の基礎を紹介する • 第2回（6月15日） 「エクサフロップス時代における線形計算アルゴリズムの課 題と研究動向」 – 様々な行列計算について，計算機上での実装における課題と，その 解決のための取り組みを紹介する – また，量子アニーリングマシンの行列計算への適用可能性について も簡単に触れる 2 # Page. 3 ![Page Image](https://bcdn.docswell.com/page/PJR224Z9E9.jpg) はじめに： スーパーコンピュータの性能動向 https://samhatfield.co.uk/2019/04/10/supercomputer-trend-data/ 3 # Page. 4 ![Page Image](https://bcdn.docswell.com/page/PEXDD2937X.jpg) 本講義の目的 • 本講義の目的 – エクサフロップス時代における数値計算アルゴリズムの課題を， 線形計算に焦点を当てて考える – 最近の研究動向について，簡単なサーベイを行う • 目次 – エクサフロップスマシンのハードウェア特性 – 線形計算アルゴリズムの技術課題 – 最近の研究動向 4 # Page. 5 ![Page Image](https://bcdn.docswell.com/page/3EKXXMPN7D.jpg) ポストペタ時代のスーパーコンピュータ： 4つのタイプ • 汎用型 「京」の路線 – メモリ容量・帯域・演算性能をバランス良く向上 – 「京」のように汎用的に様々な問題に適用可能 • 容量・帯域重視型 ベクトル計算機タイプ – 汎用型から演算性能を落として，メモリ性能により多くの資源を割く • 演算重視型 メニーコアタイプ – メモリ性能を落とし，演算性能により多くの資源を割く • メモリ容量削減型 – メモリ容量を極限まで削減し，オンチップメモリですべての計算を完結 出典： HPCI技術ロードマップ白書 5 # Page. 6 ![Page Image](https://bcdn.docswell.com/page/L7355YXZE5.jpg) 各タイプの性能諸元（予測値） • 条件 – 「京」と同程度の消費電力（20MW） – 「京」と同程度の設置面積（2000～3000m2） 最も総演算性能が高いのは演算重視型 「富岳」は汎用型 出典： HPCI技術ロードマップ白書 6 # Page. 7 ![Page Image](https://bcdn.docswell.com/page/87DDD5247G.jpg) エクサフロップスマシンの特徴（1） • 多階層の並列性 – 107 ～109 個程度の演算コア（動作周波数はGHzオーダー） – コアレベル，チップレベル，ノードレベル，システムレベルの並列性 • 複雑なメモリ階層 – – – – – コア毎のレジスタ コア毎のキャッシュメモリ オンチップメモリ ノード内共有メモリ 他ノードのメモリ 1～10チップ / ノード 共有メモリ 102～103コア / チップ メモリ メモリ 105～106ノード / システム 7 # Page. 8 ![Page Image](https://bcdn.docswell.com/page/VJPXXG1V78.jpg) エクサフロップスマシンの特徴（2） • データ移動コストの増大： スループット – 主メモリのByte/Flop（データ転送性能と演算性能の比）による比較 総演算性能 総メモリ帯域 比 演算重視型 1000～2000PFLOPS 5～10PB/s 0.005Byte/Flop 富岳*1) 400PFLOPS 160PB/s 0.37Byte/Flop 京 10PFLOPS 5PB/s 0.5Byte/Flop 主メモリアクセスの相対的コストは徐々に増加 レジスタ 演算コア データ転送速度大 キャッシュ データ転送速度小 主メモリ *1) https://www.v-t.co.jp/product/primehpc-a64fx/ 8 # Page. 9 ![Page Image](https://bcdn.docswell.com/page/2EVLL6YRJQ.jpg) エクサフロップスマシンの特徴（3） • データ移動コストの増大： レイテンシ（遅延時間） – コア間同期・通信レイテンシ： 100サイクル – ノード間通信レイテンシ： 80～200サイクル AllReduce のような全ノード同期型通信で特に影響大 80～200サイクル 9 # Page. 10 ![Page Image](https://bcdn.docswell.com/page/57GXXN96EL.jpg) エクサフロップスマシンの特徴（4） • 電力消費の問題 – 発熱抑制と節電の両面から，電力消費の抑制が重要 – オフチップのデータアクセスが大きな電力を消費 • 倍精度演算： 1.1pJ/FLOP • オンチップデータアクセス： 2.1pJ/Word • オフチップデータアクセス： 205pJ/Word 100倍 • 部品数増加に伴う故障確率上昇 – ハードエラー（熱などによる部品の故障） – ソフトエラー（α線などによるビット反転） ハード／システムソフトレベルでのエラー耐性強化は可能だが， 部品数・速度・消費電力の面で大きなコストが必要 出典： HPCI技術ロードマップ白書 10 # Page. 11 ![Page Image](https://bcdn.docswell.com/page/4EQ88QW2JP.jpg) アプリケーション並列化における要求の変化 • これまでの並列化研究： 弱スケーリングに重点 – ノード当たりの問題サイズが一定という条件下で並列化 – ノード数を増やすとともに，問題サイズもどんどん大きくする – 比較的，並列性能を出しやすい • エクサフロップス時代における問題点 – すべてのアプリで問題サイズ拡大が求められているわけではない • 問題サイズは固定し，モデルを精密化したい場合 • 問題サイズは固定し，時間ステップを増やしたい場合（分子動力学など） – 弱スケーリングでは，ノード数増加につれて計算時間も増加 • 多くの場合，演算量は問題サイズに対して線形より速く増加 • エクサフロップスマシンの性能を引き出せる問題サイズでは，計算時間 が長くなり過ぎて実用的でなくなる場合も 11 # Page. 12 ![Page Image](https://bcdn.docswell.com/page/KJ422NKPE1.jpg) 線形計算アルゴリズムの課題(1) • 107～109個のコアを活用できる並列性 • データ移動の削減（I）： データ移動量の削減 – データ移動量（階層間／ノード間）に比例してかかるコストの影響を削減 データ再利用性の向上 データが上位のメモリにある間に できるだけ集中して演算を行う 消費電力節減の面からも有効 演算コア キャッシュ 主メモリ • データ移動の削減（II）： データ移動回数の削減 – データ移動／同期1回ごとにかかるコストの影響を削減 計算粒度 計算粒度の大きいアルゴリズム 同期やデータ移動をできるだけまとめて行う 通信と通信の間で，できるだけ多くの演算を行う 12 # Page. 13 ![Page Image](https://bcdn.docswell.com/page/LE1MM5XX7G.jpg) 線形計算アルゴリズムの課題(2) • アルゴリズムレベルでの耐故障性 – 結果不正のノード，結果を返さないノードがあっても，計算が破綻せ ずに進行 • 計算量のオーダーの低減 – ある程度の（確率的）誤差を許容することでオーダーを低減 テンソルの近似に基づくアルゴリズム 確率的アルゴリズム • 強スケーリングの意味で効率的なアルゴリズム – 問題サイズを固定したとき，ノード数を増やすほど実行時間が短縮 強スケーリングの条件下では，通信・同期時間が支配的 大粒度のアルゴリズムはこの場合にも有効 13 # Page. 14 ![Page Image](https://bcdn.docswell.com/page/GEW335NKE2.jpg) 109個のコアを活用できる並列性： 密行列の場合 • 正方行列に対するアルゴリズム – LU分解，QR分解，3重対角化，Hessenberg化 並列性は十分 – 全体の計算量： O(n3)，並列性： O(n2) – ネックとなりうる部分： ピボット生成，通信 これらを隠蔽するためのスケジューリングが重要 • DAGスケジューリング • 帯行列・縦長行列に対するアルゴリズム b b – 帯行列のLU分解，帯行列の3重対角化， n – 縦長行列のQR分解，ベクトル逐次添加型の直交化 – 全体の計算量： O(n2b) or O(nb2) ，並列性： O(b2), O(nb) or O(n) 新たな並列性を導入できるアルゴリズムが必要 • • 帯行列に対する分割統治型のLU分解 Compact WY 表現に基づくベクトル逐次添加型直交化法 14 # Page. 15 ![Page Image](https://bcdn.docswell.com/page/47ZVVN3NE3.jpg) 109個のコアを活用できる並列性： 疎行列の場合 • 部分空間への射影に基づく反復解法 – 連立1次方程式の解法： Krylov部分空間法全般 – 固有値問題の解法： Lanczos法，Arnoldi法，Jacobi-Davidson法 部分空間の拡大操作は本質的に逐次的 並列性： 行列ベクトル積の並列度 O(nz) • 並列性を向上させる手法 – ブロックKrylov部分空間法 • 複数の初期ベクトルから出発し，複数の部分空間を同時に生成 • 部分空間の数だけ並列性が向上 – 数値積分を用いたフィルタで部分空間を生成する手法 • 櫻井・杉浦法 • 積分点ごとの並列性を新たに利用可能 15 # Page. 16 ![Page Image](https://bcdn.docswell.com/page/YJ6558X97V.jpg) データ移動の削減： 密行列の場合 • 想定する状況と目標 コア – コア：1個，キャッシュ： M ワード，行列サイズ： n – キャッシュと主メモリとの間のデータ移動を最小化 • ノード間通信最小化でも同様の手法を適用可能 • ブロック化（タイル）アルゴリズム レジスタ キャッシュ データ移動 最小化 主メモリ (M/3)1/2 – 行列を，サイズ (M/3)1/2×(M/3)1/2 のブロックに分割 • ブロック3個がキャッシュに入る – 各ブロックを要素と見て行列計算を行う • 行列乗算 • LU分解，コレスキー分解，QR分解 • ブロック3重対角化，ブロック Hessenberg 化 n – 演算の主要部分はブロックどうしの乗算であり，これはキャッシュ上 で実行可能 16 # Page. 17 ![Page Image](https://bcdn.docswell.com/page/GJ5YY9XDJ4.jpg) ブロック化アルゴリズムの例 • ブロック化コレスキー分解 – ブロックサイズを L= (M/3)1/2，第 (I, J) ブロックを AIJ とする do K=1, n/L AKK := Cholesky(AKK) do I=K+1, n/L AIK := AIKAKK-T end do do J=K+1, n/L do I=J, n/L AIJ := AIK AJKT end do end do end do 対角ブロックのコレスキー分解 ブロックピボット列の作成 ブロックどうしの行列乗算 （演算の主要部） – LAPACKで採用されているアルゴリズム 17 # Page. 18 ![Page Image](https://bcdn.docswell.com/page/LE355Y4QJ5.jpg) ブロック化アルゴリズムの最適性 • 定理（Ballard, Demmel, Holtz and Schwartz, 2009） – ブロックサイズを (M/3)1/2 としたブロック化コレスキー分解は，前記の 仮定の下で，キャッシュと主メモリの間のデータ移動量をオーダーの 意味で最小化する • 証明 – 行列乗算についてはデータ移動量の下界がわかっていることに着目 – コレスキー分解を用いて行列乗算を計算するアルゴリズムを構築 – これより，定数倍の差を除いて， コレスキー分解における データ移動量の下界 ≧ 行列乗算における データ移動量の下界 が言える – ブロック化コレスキー分解が右辺の下界を達成することを示す 18 # Page. 19 ![Page Image](https://bcdn.docswell.com/page/8EDDD5QWEG.jpg) Communication optimal なアルゴリズム • データ移動回数の最小化 – 1回のデータ移動では，主メモリの連続した 領域しか持ってこられないと仮定 連続でないため，1回で 持ってこられない – このとき，通常の行列格納形式では，ブロッ ク化コレスキーは移動回数最小にならない – ブロック単位の格納順を使うことで，データ 移動回数も最小化 最適でない 最適 Communication optimal なアルゴリズム • データ移動量・回数の下界が知られているアルゴリズムの例 – – – – 行列乗算（O(n3)／Strassen） LU分解，コレスキー分解 QR分解，最小2乗法 固有値分解，特異値分解 両方について下界を達成する （Communication optimal な） アルゴリズムの開発が，現在 でも活発な研究テーマ 19 # Page. 20 ![Page Image](https://bcdn.docswell.com/page/V7PXXGLXE8.jpg) データ移動の削減： 疎行列の場合 • Krylov部分空間法 – Kk(A; b) = span{b, Ab, A2b, …, Ak–1b} の中で近似解を求めていく解法 – 演算の主要部分は疎行列ベクトル積 y = Ax – 1ステップ中に複数回の内積（またはノルム計算）が存在 • データ移動の観点からの問題点 – 行列ベクトル積は行列データの再利用性が低い • y = Ax において，A の各要素は1回しか計算に利用しない – 内積は全ノードでの AllReduce を必要とし，レイテンシの影響大 • 課題 – 行列ベクトル積におけるデータ再利用性の向上 – 複数の内積をまとめて同期回数を削減 以下ではGMRES法を例として手法を紹介 20 # Page. 21 ![Page Image](https://bcdn.docswell.com/page/2JVLL6Q37Q.jpg) GMRES法 • 原理 – A をかける操作と直交化を繰り返し，部分空間を広げながら正規直 交基底 {q1, q2, q3, …}を生成 – 各ステップにおいて，||r(n)||2 = ||Ax(n) – b||2 が最小になるよう解を更新 • アルゴリズム（正規直交基底の生成部分） 新たなベクトルの生成（空間の拡張） 直交化 正規化 21 # Page. 22 ![Page Image](https://bcdn.docswell.com/page/5EGXXNWYJL.jpg) 手法（I）： ブロック GMRES 法 • アイディア – 複数（b本）の初期ベクトル R(0) = [r1(0), …, rb(0)] から出発し，ブロック Krylov 部分空間 Kb(m)(A; R(0)) 内で解を求める • 普通の GMRES 法との比較（1ステップあたり） GMRES ブロックGMRES 演算量（行列ベクトル積） 1（y=Ax） b （Y=AX） 行列データのアクセス回数 1（y=Ax） 1 （Y=AX） 1 1 同期回数* • 効果 * 直交化に compact WY 表現またはCGS法を用いた場合 再利用性 b 倍向上 同期回数 実質 1 / b – 行列アクセス回数・同期回数は同じため，実行時間増加は b 倍以下 – ブロック Krylov 部分空間を使うことによる収束性向上の効果も存在 右辺が複数本の場合は有利．1本でも高速化できる場合もある 22 # Page. 23 ![Page Image](https://bcdn.docswell.com/page/4JQ88Q367P.jpg) 手法（II）： k-step GMRES 法 • アイディア – 行列ベクトル積 Ar(m), A2r(m) …, Akr(m) を一度に行って Krylov 部分空間 を一度に k 次元拡大し，その後に正規直交基底を生成する • 普通の GMRES 法との比較（1ステップあたり） GMRES k-step GMRES 演算量（行列ベクトル積） 1 1+a 行列データのアクセス回数 1 1/k 同期回数* 1 1/k • 効果 * 直交化に compact WY 表現またはCGS法を用いた場合 再利用性 k 倍向上 同期回数 1/k – 再利用性向上と同期の削減により，大きな性能向上が期待できる – ただし，（直交化前の）基底が線形従属に近い場合は，収束性が悪化 線形独立性を高めるため，A の単項式の代わりに直交多項式を 使い，3項間漸化式で計算するなどの工夫が検討されている 23 # Page. 24 ![Page Image](https://bcdn.docswell.com/page/K7422N1MJ1.jpg) アルゴリズムレベルでの耐故障性 • 想定する状況と目標 – 複数のノードが通信を行いつつ協調して計算 – そのうち1個のが不正な結果を返すか，あるいは 結果を返さなくても，計算は破綻せずに進行 – その場合，精度劣化，収束性劣化は許容 – 計算の一部を，高信頼モード／高信頼ハード ウェア（ただし計算コスト大）で行ってもよい × • 考察 – 複数のノードの結果を集めて部分空間を改良するタイプのアルゴリ ズムは，耐故障性と親和性が高い – 大粒度並列性は，耐故障性にとっても有利 • 高信頼モードの使用頻度を削減できる 24 # Page. 25 ![Page Image](https://bcdn.docswell.com/page/LJ1MM5GYEG.jpg) MERAM（Multiple Explicitly Restarted Arnoldi Method） • アルゴリズム（Emad，Petiton & Edjlali, 2005） (1) P 個のノードで，異なる初期ベクトルを用いて独立にArnoldi法を実行 (2) 各ノードで作ったKrylov部分空間を合わせ，大きな部分空間を生成 (3) その中で P 本の初期ベクトルを新たに生成し，Arnoldi法をリスタート • 耐故障性 – 上記(1)を低信頼モード，(2)，(3)を高信頼モードで実行 • 計算量の大部分は(1)のステップ – (1) において，結果不正のノード／結果を返さないノードがあっても， 収束性が落ちるだけで，計算は破綻しない ノード0 Arnoldi ノード1 Arnoldi ノード2 Arnoldi 低信頼モード 部分空間合成 P本のリスタート ベクトルの生成 Arnoldi Arnoldi Arnoldi 部分空間合成 P本のリスタート ベクトルの生成 Arnoldi Arnoldi Arnoldi 高信頼モード 25 # Page. 26 ![Page Image](https://bcdn.docswell.com/page/GJW335K1J2.jpg) 耐故障性を持つその他のアルゴリズム • Fault-Tolerant GMRES法（Hoemmen & Heroux, 2011） – Flexible GMRES法の枠組みを利用 – 結果不正を，不適切な前処理と解釈して計算を続行 • そのステップでは無駄に部分空間が大きくなるが，精度には悪影響なし Im z • 櫻井・杉浦法 – ある積分点を担当するノードが結果を返さない場合， その点を使わない積分公式を新たに構築して使用 O × Re z • 密行列向けの解法 – 直接解法では，1箇所の計算間違いで結果に壊滅的な影響 – アルゴリズムレベルでの耐故障性の実現は，原理的に難しそう 26 # Page. 27 ![Page Image](https://bcdn.docswell.com/page/4EZVVNZXJ3.jpg) 確率的アルゴリズムによる計算量のオーダー低減 • 特異値分解： 行列の最良の低ランク近似を与える分解 – 画像処理，信号処理，情報検索など広い応用 – m×n 行列（m≧n）の場合の演算量は O(mn2) と大きい • CX分解： 特異値分解の代替手法 – C： A の列ベクトルを k 本（1≦ k ≦ n）選んでできる m×k 行列 – X： 適当な k×n 行列 – このとき，∥A – CX∥F をできるだけ小さくする C と X を求める A の特徴を最もよく表す s 本の列ベクトルを選ぶことに相当 Fノルムに対しては，X = C+A と選ぶのが最適 したがって，問題は C を選ぶことに帰着 27 # Page. 28 ![Page Image](https://bcdn.docswell.com/page/Y76558ZPEV.jpg) Leverage score に基づく確率的アルゴリズム • Leverage score（Drineas et al., 2008） – A の打ち切り型特異値分解を Ak = Uk Sk VkT とする – このとき，次の量を，VkT の第 j 列の Leverage score と呼ぶ • サンプリングと誤差評価 – 確率 pj に従い，A の列をサンプリングする – このとき，確率 0.9 以上で次の誤差評価が成り立つ • 問題点 サンプリング 相対評価 回数に依存 – Leverage score の計算には，VkT が必要（このままでは非実用的） 28 # Page. 29 ![Page Image](https://bcdn.docswell.com/page/G75YY9WP74.jpg) Leverage score の近似 • 近似アルゴリズム（Clarkson et al., 2012） – O(mn log m) の計算量で，Leverage score を近似的に計算するアルゴ リズムが存在する • アイディア – A に左から適当な直交行列をかけることで，一様な Leverage score を 持つ行列 A’ に変換 • Johnson & Lindenstrauss の補題と高速 Walsh-Hadamard 変換を使う – A’ に対して一様なサンプリングを行い，CX 分解を求める – A’ のCX 分解から，A のCX 分解を求める 相対誤差の意味での（確率的）誤差評価を持つCX分解の計算が， O(mn log m) で可能に 29 # Page. 30 ![Page Image](https://bcdn.docswell.com/page/9J2GG2QZJR.jpg) 強スケーリングの意味で効率的なアルゴリズム • 実対称行列の全固有値・固有ベクトル計算 – 分子軌道法をはじめ，様々な問題で重要 – n=10,000 程度の中規模問題に対する需要も多い • 想定する状況と目標 – 行列サイズは n=10,000 に固定 – ノードは何個使ってもよいから，できるだけ高速に解きたい • 分子軌道法では，多電子積分の計算量が O(n4) 以上で並列性も高い • この部分を並列化するため，1万ノード以上を確保している例もある • 標準的な行列計算ライブラリ（ScaLAPACK）での結果 – 「京」上では，ノード数を400以上にしても，計算が加速しない • 実行時間の70%以上を占める通信オーバーヘッドのため – 確保した1万ノードの大部分は，固有値計算部分では遊んでしまう 30 # Page. 31 ![Page Image](https://bcdn.docswell.com/page/DEYVVYWNEM.jpg) 強スケーリングの意味で効率的なアルゴリズム • ブロックヤコビ法による全固有値・固有ベクトル計算 – 行列をブロックに分割し，直交変換により複数の非対角ブロックを並 列に消去 – これを繰り返すことで，行列を対角行列に近づける = • アルゴリズムの特性 – 演算量は3重対角化に基づく解法の10倍以上 – Pノードで2次元分割を行った場合，通信回数は O(P1/2 log P * Iter) n=P=10,000 ならば3重対角化の通信回数 O(N log P) に比べて小さい 通信OHの大きい状況では，3重対角化法より高速となる可能性 31 # Page. 32 ![Page Image](https://bcdn.docswell.com/page/VJNMMP9V78.jpg) 強スケーリングの意味で効率的なアルゴリズム • ScaLAPACK とブロックヤコビ法の性能比較 – 「京」上で1万ノードまでを使用して実行時間を測定 – 行列サイズは n=10,000 に固定 35 強スケーリングの条件下では ブロックヤコビ法が優位 消去順序の最適化，前処理に より，更に高速化の可能性 実 行 時 間 （ 秒 ） ScaLAPACK Block Jacobi 30 25 20 15 10 8.3s 5 5.5s 0 100 400 625 1600 2500 10000 ノード数 • 強スケーリングの条件下でのアルゴリズム設計 – 通信・同期オーバーヘッドの削減が重要 – そのためならば，演算量をある程度増やしてもよい • 中務氏の新しい固有値・特異値計算アルゴリズム 32 # Page. 33 ![Page Image](https://bcdn.docswell.com/page/YE9NN32VE3.jpg) 大粒度並列性を持つQR分解手法 • A ∈ Rm×n （m≧n）のQR分解 – full QR分解 A = QR A = Q A = Q R ~~ – thin QR分解 A = QR ~ ~ R • QR分解に対する2つの見方 ~ – Q に着目 → A の列ベクトルの正規直交化 – R に着目 → 直交行列による A の上三角化 様々な応用 • QR分解のアルゴリズム – ハウスホルダー法 – 古典／修正グラム・シュミット法 高性能計算の観点からは課題が多い 33 # Page. 34 ![Page Image](https://bcdn.docswell.com/page/GE8KKM547D.jpg) CholeskyQR2法による大粒度並列QR分解 • Cholesky QR 法による行列 X∈Rm×n （m≧n）のQR分解 – アルゴリズム （コレスキー分解： A = RTR） – 長所： 大粒度並列，すべての計算がブロック化可能 – 短所： X の条件数が大きいとき，直交行列 Y の直交性が悪化 • CholeskyQR2 法*1 – Cholesky QR 法で得られた Y を再直交化 • Y を X と見て，上記のアルゴリズムをもう一度繰り返す – X の条件数が 108 以下ならば，極めて精度の良い結果が得られる • 直交性 ||QTQ – I|| も残差 ||X – YR|| も丸め誤差レベル *1 Y. Yamamoto et al.: “Roundoff Error Analysis of the CholeskyQR2 Algorithm”, Electronic Transactions on Numerical Analysis, Vol. 44, pp. 306-326 (2015). 34 # Page. 35 ![Page Image](https://bcdn.docswell.com/page/LELZZ3YVJR.jpg) 「京」上でのCholeskyQR2法の性能 ScaLAPACK (Householder QR) (sec.) m=4,194,304 n=16 x2.0 speedup over TSQR P: #processes (= #nodes) (sec.) TSQR m=4,194,304 n=64 x2.9 speedup over TSQR P: #processes (= #nodes) CholeskyQR2 (sec.) m=4,194,304 n=256 x7.6 speedup over TSQR P: #processes (= #nodes) • TSQRとCholeskyQR2はScaLAPACK（ハウスホルダーQR）より高速 （n=256でPが大きい場合を除く） • CholeskyQR2はTSQRよりも高速 （Pが大きい場合だけでなく，Pが小さい場合でも） コレスキーQR分解を2回繰り返してもまだTSQRより高速！ 35 # Page. 36 ![Page Image](https://bcdn.docswell.com/page/4JMVVDNNJW.jpg) CholeskyQR2法の応用 • 第一原理分子動力学法における波動関数の時間発展 • 時間離散化を行った場合 – Xt = [ψ 1(r, t), …, ψ N(r, t)] とおくと， ~ Xt+Dt = Xt +CDt ~ Xt+Dt = Xt+Dt R 時間発展 直交化（QR分解） ~ – Xt は直交化されているから，Dt が小さければ， Xt+Dt も直交に近い – すなわち，条件数は小さい（列が直交する場合，条件数は1） CholeskyQR2法が適用可能のはず 36 # Page. 37 ![Page Image](https://bcdn.docswell.com/page/PJR224D4E9.jpg) 悪条件の行列のための shifted CholeskyQR3法 • Shifted CholeskyQR3 法*1 – 行列 X が悪条件の場合，Cholesky QR2 法において A = X TX が特異 に近くなり，コレスキー分解が破綻することがある – A に微小なシフト sI を加えることで，条件数 1016 程度まで対応可能 – 大粒度並列性など，HPC面から見た利点は Cholesky QR2 法と同じ *1 T. Fukaya et al.: “Shifted CholeskyQR for computing the QR factorization of ill-conditioned matrices”, SIAM Journal on Scientific Computing, Vol. 42, , No. 1, pp. A477-A503 (2020). 37 # Page. 38 ![Page Image](https://bcdn.docswell.com/page/PEXDD21Z7X.jpg) 不完全LU分解型前処理の超並列化 • クリロフ部分空間法と前処理（前回の復習） – 連立1次方程式 K–1Ax = K–1b を考える – いま，行列 K が A の近似行列であるとする – このとき，方程式 K–1Ax = K–1b を考えると，その係数行列 K–1A は A と比べて単位行列に近くなっていると考えられる．これにより，クリロ フ部分空間法の収束はより速くなると期待できる – このような方程式の変形を前処理（左前処理）と呼ぶ • 不完全LU分解型前処理（ILU(0)前処理） – 行列 A のLU分解の際に，フィルインを無視して近似的に A≒LU と 分解し，K = LU とおく．これを ILU(0) 前処理と呼ぶ – ILU(0) 前処理は，幅広いクラスの行列に対して効果があるが，一般 に並列化は困難 • K–1 をかける計算で，前進消去と後退代入が必要となるため 38 # Page. 39 ![Page Image](https://bcdn.docswell.com/page/3EKXXM2V7D.jpg) ブロック赤黒順序付けによる並列化（1） • 概要 – 偏微分方程式を格子点上で離散化して得られる行列を考える – この場合，格子点の番号の付け替え（リオーダリング）により，ILU(0) 分解の並列性を抽出することが可能 – 代表的なリオーダリング手法として，赤黒ブロック順序付けがある 4×4 格子の ブロック分割数 x 方向: 2， y 方向: 2 の順序付けの例 ・ 青字： 自然な格子点の番号 ・ 赤黒字： リオーダリング後の格子点 の番号 岩下武史, 島崎眞昭: 同期点の少ない並列化 ICCG 法のためのブロック化赤-黒 順序付け 情報処理学会論文誌, Vol.43, No.4，pp.893-904，2002． 39 # Page. 40 ![Page Image](https://bcdn.docswell.com/page/L7355YYQE5.jpg) ブロック赤黒順序付けによる並列化（2） ~ • リオーダリングは係数行列 A の行と列の置換 A = PAPT に相当する ~ • ブロック赤黒順序付けによって得られる行列 A は下図のようになる • 同色のブロック間に依存関係がないため，各プロセスに1ブロックを担当 させることで，ILU(0) 分解と前進消去・後退代入が並列化できる 40 # Page. 41 ![Page Image](https://bcdn.docswell.com/page/87DDD55W7G.jpg) ブロック赤黒順序付けによる並列化（3） • ブロック数を増やした場合の収束性 – 120×120×60 格子でのポアソン方程式の求解 – 並列数 = 各色のブロック数 = 全ブロック数 / 2 全 ブロック数が10000程度までは，あまり反復回数が増えない 超並列向き 41 # Page. 42 ![Page Image](https://bcdn.docswell.com/page/VJPXXGGX78.jpg) 低精度計算の利用による高性能化（1） • 近年のCPU/GPUの低精度演算性能 富士通 A64FX（富岳） NVIDIA Quadro GV100 倍精度 3.072TFLOPS 7.40TFLOPS 単精度 6.144TFLOPS 14.80TFLOPS 半精度 12.288TFLOPS 29.60TFLOPS 低精度演算の活用により高速化を達成できる可能性あり • 前処理付きクリロフ部分空間法と低精度演算 – 前処理付きクリロフ部分空間法では，残差の正確性（r(n) = Ax(n) – b が どの程度正確に成り立つか）が最も重要 – この関係さえ保持しておけば，一部で低精度演算を使ったとしても， 解の精度は保たれる • ただし，収束性は低下する可能性あり 42 # Page. 43 ![Page Image](https://bcdn.docswell.com/page/2EVLL663JQ.jpg) 低精度計算の利用による高性能化（2） • 右前処理に低精度計算を使った BiCGSTAB 法1) 右前処理： AM y = b, x = M y 低精度計算による uk, sk の誤差 低精度演算 低精度演算 低精度計算による xk+1 の誤差 低精度計算による rk+1 の誤差 xk+1, rk+1 はそれぞれ誤差を含むが，関係式 rk+1 = Axk+1– b は保たれている 1) H. Tadano et al., “On Single Precision Preconditioners for Krylov Subspace Iterative Methods”, Proceedings of LSSC 2007, LNCS 4818, pp. 721-728, 2008. 43 # Page. 44 ![Page Image](https://bcdn.docswell.com/page/57GXXNNYEL.jpg) 低精度計算の利用による高性能化（3） • 右前処理に低精度計算を使った BiCGSTAB 法の数値結果 – 前処理としてブロック赤黒順序付けのMILU(0) 分解を使用 – ベクトルに前処理行列を適用する部分を単精度で計算 – 倍精度とほぼ同じ結果が得られている A. Shioya and Y. Yamamoto: “Block red–black MILU(0) preconditioner with relaxation on GPU”, Parallel Computing, Vol. 103, No. 1, pp. 102760-102772 (2021). 44 # Page. 45 ![Page Image](https://bcdn.docswell.com/page/4EQ88QQ6JP.jpg) ブロック赤黒順序付けに基づくソルバのGPU性能 • ソルバの概要 – ブロック赤黒順序付けによる並列 MILU(0) 分解 + BiCGSTAB法 – コアレスアクセスのためにデータ格納形式を最適化 – 前処理部分で単精度計算を利用 • NVIDIA Tesla K40t 上での性能 – 例題： 3次元ポアソン方程式 59×59×29 119×119×59 239×239×119 45 # Page. 46 ![Page Image](https://bcdn.docswell.com/page/KJ422NNME1.jpg) おわりに • 本発表のまとめ – エクサフロップスマシンでは，多階層の超並列性，データ移動コスト の増大，故障確率の上昇などが課題となる．また，強スケーリングで の並列化効率がより重視される． – そのため，線形計算アルゴリズムの研究では次の点が重要になる． • 109レベルの並列性 • データ移動量・移動回数の削減 • アルゴリズムレベルでの耐故障性 • （確率的）近似による計算量のオーダーの削減 • 強スケーリング性 – これらの方向に沿った最近の研究をいくつか紹介した． 46