行列計算における高速アルゴリズム― ポストペタ時代に向けた線形計算アルゴリズムの課題と研究動向 ―2017年6月29日

計算科学技術特論A
電気通信大学 情報理工学研究科 情報・ネットワーク工学専攻
山本 有作

はじめに： スーパーコンピュータの性能動向
性能向上： 10年で1000倍

2018年にはエクサフロップスマシンが登場の見込み

2

本発表の目的• 本発表の目的– エクサフロップス時代に向けた数値計算アルゴリズムの課題を，線形計算に焦点を当てて考える– 最近の研究動向について，簡単なサーベイを行う• 目次– エクサフロップスマシンのハードウェア特性– 線形計算アルゴリズムの技術課題– 最近の研究動向

32018年のスーパーコンピュータ： 4つのタイプ• 汎用型「京」の後継– メモリ容量・帯域・演算性能をバランス良く向上– 「京」のように汎用的に様々な問題に適用可能• 容量・帯域重視型

ベクトル計算機タイプ– 汎用型から演算性能を落として，メモリ性能により多くの資源を割く• 演算重視型

GPU・メニーコアタイプ– メモリ性能を落とし，演算性能により多くの資源を割く• メモリ容量削減型– メモリ容量を極限まで削減し，オンチップメモリですべての計算を完結
出典： HPCI技術ロードマップ白書，2012年3月

4各タイプの性能諸元（予測値）• 条件– 「京」と同程度の消費電力（20MW）– 「京」と同程度の設置面積（2000～3000m2）

最もエクサフロップスに近いのは演算重視型
出典： HPCI技術ロードマップ白書，2012年3月

5エクサフロップスマシンの特徴（1）• 多階層の並列性– 109 個程度の演算コア（動作周波数はGHzオーダー）– コアレベル，チップレベル，ノードレベル，システムレベルの並列性• 複雑なメモリ階層– コア毎のレジスタ– コア毎のキャッシュメモリ– オンチップメモリ– ノード内共有メモリ– 他ノードのメモリ

10チップ / ノード

1000コア / チップ共有メモリ

メモリ

メモリ

105ノード / システム6

エクサフロップスマシンの特徴（2）• データ移動コストの増大： スループット– 主メモリのByte/Flop（データ転送性能と演算性能の比）による比較
総演算性能

総メモリ帯域比演算重視型

1000～2000PFLOPS5～10PB/s

0.005Byte/Flop

京

10PFLOPS

5PB/s

0.5Byte/Flop

主メモリアクセスの相対的コストは，京の場合の100倍
レジスタ演算コア
データ転送速度大キャッシュ
データ転送速度小主メモリ
7

エクサフロップスマシンの特徴（3）• データ移動コストの増大： レイテンシ（遅延時間）– コア間同期・通信レイテンシ： 100サイクル– ノード間通信レイテンシ： 80～200サイクル

AllReduce のような全ノード同期型通信で特に影響大80～200サイクル
8エクサフロップスマシンの特徴（4）• 電力消費の問題– 発熱抑制と節電の両面から，電力消費の抑制が重要– オフチップのデータアクセスが大きな電力を消費• 倍精度演算： 1.1pJ/FLOP• オンチップデータアクセス： 2.1pJ/Word• オフチップデータアクセス： 205pJ/Word100倍• 部品数増加に伴う故障確率上昇– ハードエラー（熱などによる部品の故障）– ソフトエラー（α線などによるビット反転）
ハード／システムソフトレベルでのエラー耐性強化は可能だが，部品数・速度・消費電力の面で大きなコストが必要
出典： HPCI技術ロードマップ白書，2012年3月

9アプリケーション並列化における要求の変化• これまでの並列化研究： 弱スケーリングに重点– ノード当たりの問題サイズが一定という条件下で並列化– ノード数を増やすとともに，問題サイズもどんどん大きくする– 比較的，並列性能を出しやすい• エクサフロップス時代における問題点– すべてのアプリで問題サイズ拡大が求められているわけではない• 問題サイズは固定し，モデルを精密化したい場合• 問題サイズは固定し，時間ステップを増やしたい場合（分子動力学など）– 弱スケーリングでは，ノード数増加につれて計算時間も増加• 多くの場合，演算量は問題サイズに対して線形より速く増加• エクサフロップスマシンの性能を引き出せる問題サイズでは，計算時間が長くなり過ぎて実用的でなくなる場合も10

線形計算アルゴリズムの課題(1)• 109個のコアを活用できる並列性• データ移動の削減（I）： データ移動量の削減– データ移動量（階層間／ノード間）に比例してかかるコストの影響を削減
データ再利用性の向上
データが上位のメモリにある間にできるだけ集中して演算を行う
消費電力節減の面からも有効

演算コア
キャッシュ
主メモリ

• データ移動の削減（II）： データ移動回数の削減– データ移動／同期1回ごとにかかるコストの影響を削減
計算粒度
計算粒度の大きいアルゴリズム
同期やデータ移動をできるだけまとめて行う
通信と通信の間で，できるだけ多くの演算を行う
11

線形計算アルゴリズムの課題(2)• アルゴリズムレベルでの耐故障性– 結果不正のノード，結果を返さないノードがあっても，計算が破綻せずに進行• 計算量のオーダーの低減– ある程度の（確率的）誤差を許容することでオーダーを低減
テンソルの近似に基づくアルゴリズム
確率的アルゴリズム• 強スケーリングの意味で効率的なアルゴリズム– 問題サイズを固定したとき，ノード数を増やすほど実行時間が短縮

強スケーリングの条件下では，通信・同期時間が支配的
大粒度のアルゴリズムはこの場合にも有効12

109個のコアを活用できる並列性： 密行列の場合• 正方行列に対するアルゴリズム– LU分解，QR分解，3重対角化，Hessenberg化
並列性は十分– 全体の計算量： O(n3)，並列性： O(n2)– ネックとなりうる部分： ピボット生成，通信

これらを隠蔽するためのスケジューリングが重要•

DAGスケジューリング

• 帯行列・縦長行列に対するアルゴリズム

bb– 帯行列のLU分解，帯行列の3重対角化，n– 縦長行列のQR分解，ベクトル逐次添加型の直交化– 全体の計算量： O(n2b) or O(nb2) ，並列性： O(b2), O(nb) or O(n)新たな並列性を導入できるアルゴリズムが必要••

帯行列に対する分割統治型のLU分解
Compact WY 表現に基づくベクトル逐次添加型直交化法
13109個のコアを活用できる並列性： 疎行列の場合• 部分空間への射影に基づく反復解法– 連立1次方程式の解法： Krylov部分空間法全般– 固有値問題の解法： Lanczos法，Arnoldi法，Jacobi-Davidson法
部分空間の拡大操作は本質的に逐次的並列性： 行列ベクトル積の並列度 O(nz)• 並列性を向上させる手法– ブロックKrylov部分空間法• 複数の初期ベクトルから出発し，複数の部分空間を同時に生成• 部分空間の数だけ並列性が向上– 数値積分を用いたフィルタで部分空間を生成する手法• 櫻井・杉浦法• 積分点ごとの並列性を新たに利用可能
14データ移動の削減： 密行列の場合• 想定する状況と目標

コア– コア：1個，キャッシュ： M ワード，行列サイズ： n– キャッシュと主メモリとの間のデータ移動を最小化• ノード間通信最小化でも同様の手法を適用可能• ブロック化（タイル）アルゴリズム

レジスタ

キャッシュ
データ移動最小化
主メモリ(M/3)1/2– 行列を，サイズ (M/3)1/2×(M/3)1/2 のブロックに分割• ブロック3個がキャッシュに入る– 各ブロックを要素と見て行列計算を行う• 行列乗算• LU分解，コレスキー分解，QR分解• ブロック3重対角化，ブロック Hessenberg 化

n– 演算の主要部分はブロックどうしの乗算であり，これはキャッシュ上で実行可能
15

ブロック化アルゴリズムの例• ブロック化コレスキー分解– ブロックサイズを 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 dodo J=K+1, n/L
do I=J, n/L
AIJ := AIK AJKT
end doend doend do

対角ブロックのコレスキー分解
ブロックピボット列の作成

ブロックどうしの行列乗算（演算の主要部）– LAPACKで採用されているアルゴリズム
16

ブロック化アルゴリズムの最適性• 定理（Ballard, Demmel, Holtz and Schwartz, 2009）– ブロックサイズを (M/3)1/2 としたブロック化コレスキー分解は，前記の仮定の下で，キャッシュと主メモリの間のデータ移動量をオーダーの意味で最小化する• 証明– 行列乗算についてはデータ移動量の下界がわかっていることに着目– コレスキー分解を用いて行列乗算を計算するアルゴリズムを構築– これより，定数倍の差を除いて，コレスキー分解におけるデータ移動量の下界≧行列乗算におけるデータ移動量の下界が言える– ブロック化コレスキー分解が右辺の下界を達成することを示す
17

Communication optimal なアルゴリズム• データ移動回数の最小化– 1回のデータ移動では，主メモリの連続した領域しか持ってこられないと仮定

連続でないため，1回で持ってこられない– このとき，通常の行列格納形式では，ブロック化コレスキーは移動回数最小にならない– ブロック単位の格納順を使うことで，データ移動回数も最小化

最適でない

最適

Communication optimal なアルゴリズム• データ移動量・回数の下界が知られているアルゴリズムの例– 行列乗算（O(n3)／Strassen）– LU分解，コレスキー分解– QR分解，最小2乗法– 固有値分解，特異値分解両方について下界を達成する（Communication optimal な）アルゴリズムの開発が，活発な研究テーマ
18データ移動の削減： 疎行列の場合• Krylov部分空間法– Kk(A; b) = span{b, Ab, A2b, …, Ak–1b} の中で近似解を求めていく解法– 演算の主要部分は疎行列ベクトル積 y = Ax
– 1ステップ中に複数回の内積（またはノルム計算）が存在• データ移動の観点からの問題点– 行列ベクトル積は行列データの再利用性が低い• y = Ax において，A の各要素は1回しか計算に利用しない– 内積は全ノードでの AllReduce を必要とし，レイテンシの影響大• 課題– 行列ベクトル積におけるデータ再利用性の向上– 複数の内積をまとめて同期回数を削減

以下ではGMRES法を例として手法を紹介
19

GMRES法• 原理– A をかける操作と直交化を繰り返し，部分空間を広げながら正規直交基底 {q1, q2, q3, …}を生成– 各ステップにおいて，||r(n)||2 = ||Ax(n) – b||2 が最小になるよう解を更新• アルゴリズム（正規直交基底の生成部分）新たなベクトルの生成（空間の拡張）
直交化
正規化
20

手法（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本でも高速化できる場合もある
21

手法（II）： k-step GMRES 法• アイディア
– 行列ベクトル積 Ar(m), A2r(m) …, Akr(m) を一度に行って Krylov 部分空間を一度に k 次元拡大し，その後に正規直交基底を生成する• 普通の GMRES 法との比較（1ステップあたり）
GMRES

k-step GMRES

演算量（行列ベクトル積）11+a行列データのアクセス回数

11/k

同期回数*

11/k

• 効果

* 直交化に compact WY 表現またはCGS法を用いた場合

再利用性
k 倍向上
同期回数1/k– 再利用性向上と同期の削減により，大きな性能向上が期待できる– ただし，（直交化前の）基底が線形従属に近い場合は，収束性が悪化
線形独立性を高めるため，A の単項式の代わりに直交多項式を使い，3項間漸化式で計算するなどの工夫が検討されている

22

アルゴリズムレベルでの耐故障性• 想定する状況と目標– 複数のノードが通信を行いつつ協調して計算– そのうち1個のが不正な結果を返すか，あるいは結果を返さなくても，計算は破綻せずに進行– その場合，精度劣化，収束性劣化は許容– 計算の一部を，高信頼モード／高信頼ハードウェア（ただし計算コスト大）で行ってもよい×• 考察– 複数のノードの結果を集めて部分空間を改良するタイプのアルゴリ
ズムは，耐故障性と親和性が高い– 大粒度並列性は，耐故障性にとっても有利• 高信頼モードの使用頻度を削減できる
23

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本のリスタートベクトルの生成ArnoldiArnoldiArnoldi

部分空間合成P本のリスタートベクトルの生成ArnoldiArnoldiArnoldi

高信頼モード24耐故障性を持つその他のアルゴリズム• Fault-Tolerant GMRES法（Hoemmen & Heroux, 2011）– Flexible GMRES法の枠組みを利用– 結果不正を，不適切な前処理と解釈して計算を続行• そのステップでは無駄に部分空間が大きくなるが，精度には悪影響なし
Im z• 櫻井・杉浦法– ある積分点を担当するノードが結果を返さない場合，その点を使わない積分公式を新たに構築して使用

O×Re z• 密行列向けの解法– 直接解法では，1箇所の計算間違いで結果に壊滅的な影響– アルゴリズムレベルでの耐故障性の実現は，原理的に難しそう
25確率的アルゴリズムによる計算量のオーダー低減• 特異値分解： 行列の最良の低ランク近似を与える分解– 画像処理，信号処理，情報検索など広い応用– 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 を選ぶことに帰着
26

Leverage score に基づく確率的アルゴリズム• Leverage score（Drineas et al., 2008）– A の打ち切り型特異値分解を A = Uk Sk VkT とする– このとき，次の量を，VkT の第 j 列の Leverage score と呼ぶ• サンプリングと誤差評価– 確率 pj に従い，A の列をサンプリングする– このとき，確率 0.9 以上で次の誤差評価が成り立つ• 問題点

サンプリング 相対評価
回数に依存

– Leverage score の計算には，VkT が必要（このままでは非実用的）
27

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) で可能に28

強スケーリングの意味で効率的なアルゴリズム• 実対称行列の全固有値・固有ベクトル計算– 分子軌道法をはじめ，様々な問題で重要– n=10,000 程度の中規模問題に対する需要も多い• 想定する状況と目標– 行列サイズは n=10,000 に固定– ノードは何個使ってもよいから，できるだけ高速に解きたい• 分子軌道法では，多電子積分の計算量が O(n4) 以上で並列性も高い• この部分を並列化するため，1万ノード以上を確保している例もある• 標準的な行列計算ライブラリ（ScaLAPACK）での結果– 「京」上では，ノード数を400以上にしても，計算が加速しない• 実行時間の70%以上を占める通信オーバーヘッドのため– 確保した1万ノードの大部分は，固有値計算部分では遊んでしまう29

強スケーリングの意味で効率的なアルゴリズム• ブロックヤコビ法による全固有値・固有ベクトル計算– 行列をブロックに分割し，直交変換により複数の非対角ブロックを並列に消去– これを繰り返すことで，行列を対角行列に近づける=• アルゴリズムの特性– 演算量は3重対角化に基づく解法の10倍以上– Pノードで2次元分割を行った場合，通信回数は O(P1/2 log P * Iter)n=P=10,000 ならば3重対角化の通信回数 O(N log P) に比べて小さい
通信OHの大きい状況では，3重対角化法より高速となる可能性
30

強スケーリングの意味で効率的なアルゴリズム• ScaLAPACK とブロックヤコビ法の性能比較（Kudo, 2013）– 「京」上で1万ノードまでを使用して実行時間を測定– 行列サイズは n=10,000 に固定
35

強スケーリングの条件下ではブロックヤコビ法が優位
消去順序の最適化，前処理により，更に高速化の可能性

ScaLAPACK
Block Jacobi

30

実
行 25
時 20間 15（秒 10）

8.3s

55.5s

0100

400

625

1600

2500

10000

ノード数

• 強スケーリングの条件下でのアルゴリズム設計– 通信・同期オーバーヘッドの削減が重要– そのためならば，演算量をある程度増やしてもよい• 中務氏の新しい固有値・特異値計算アルゴリズム
31

大粒度並列性を持つQR分解手法• A ∈ Rm×n （m≧n）のQR分解– full QR分解 A = QR
A=QA= QR~~

– thin QR分解 A = QR~~

R

• QR分解に対する2つの見方~– Q に着目 → A の列ベクトルの正規直交化– R に着目 → 直交行列による A の上三角化

様々な応用• QR分解のアルゴリズム– ハウスホルダー法– 古典／修正グラム・シュミット法

高性能計算の観点からは課題が多い
32

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).

33

「京」上でのCholeskyQR2法の性能
ScaLAPACK (Householder QR)(sec.)

m=4,194,304 n=16x2.0 speedup
over TSQR

P: #processes (= #nodes)(sec.)

TSQR

m=4,194,304 n=64x2.9 speedup
over TSQR

P: #processes (= #nodes)

CholeskyQR2(sec.)

m=4,194,304 n=256x7.6 speedup
over TSQR

P: #processes (= #nodes)• TSQRとCholeskyQR2はScaLAPACK（ハウスホルダーQR）より高速（n=256でPが大きい場合を除く）• CholeskyQR2はTSQRよりも高速（Pが大きい場合だけでなく，Pが小さい場合でも）コレスキーQR分解を2回繰り返してもまだTSQRより高速！
34

CholeskyQR2法の応用• 第一原理分子動力学法における波動関数の時間発展• 時間離散化を行った場合– Xt = [ψ 1(r, t), …, ψ N(r, t)] とおくと，~

Xt+Dt = Xt +CDt~
Xt+Dt = Xt+Dt R

時間発展
直交化（QR分解）~– Xt は直交化されているから，Dt が小さければ， Xt+Dt も直交に近い– すなわち，条件数は小さい（列が直交する場合，条件数は1）
CholeskyQR2法が適用可能のはず
35

おわりに
• 本発表のまとめ– エクサフロップスマシンでは，多階層の超並列性，データ移動コストの増大，故障確率の上昇などが課題となる．また，強スケーリングでの並列化効率がより重視される．
– そのため，線形計算アルゴリズムの研究では次の点が重要になる．
• 109レベルの並列性• データ移動量・移動回数の削減• アルゴリズムレベルでの耐故障性• （確率的）近似による計算量のオーダーの削減• 強スケーリング性– これらの方向に沿った最近の研究をいくつか紹介した．

36

SIAM Conference on Applied Linear Algebra• 会議概要– 線形代数全般に関する（たぶん）世界最大の国際会議– Communication-avoiding など，HPCに関連する発表も多い– 本発表でも，この会議で聞いた様々な講演を参考にした• 出張報告– http://www.na.scitec.kobe-u.ac.jp/~yamamoto/conference.html

37