Transcript Document

応用数理工学特論
線形計算と
ハイパフォーマンスコンピューティング
第7回
名古屋大学 計算理工学専攻
山本有作
目次
１. 固有値計算の概要
２. 三重対角化処理の並列化
３. 三重対角化処理のキャッシュ向け最適化
１. 固有値計算の概要
• 固有値計算とその応用
• 固有値・固有ベクトル計算の流れ
• 各部分の演算量
• 三重対角行列 T に対する固有値・固有ベクトル解法
固有値計算とその応用
• 対象とする問題
– 標準固有値問題 Ax = λx （A： 実対称またはエルミート行列）
• 対象とする計算機
– 共有メモリ型並列計算機
– 分散メモリ型並列計算機
– SMPクラスタ
• 固有値計算の応用分野
– 分子計算，統計計算，構造解析など
固有値・固有ベクトル計算の流れ
実対称行列 A
計算内容
計算手法
三重対角化
Q*AQ = T
(Q: 直交行列)
ハウスホルダー法
三重対角行列 T
三重対角行列の
固有値・固有ﾍﾞｸﾄﾙ計算
Tyi =λi yi
Tの固有値 {λi }，
固有ベクトル {yi }
逆変換
Aの固有ベクトル {xi }
xi = Qyi
QR法
DC法
二分法・逆反復法
Dhillon のｱﾙｺﾞﾘｽﾞﾑ
逆変換
各部分の演算量
• N×N の実対称行列 A に対し，M 組の固有値・固有ベク
トルを求める場合
– ハウスホルダー法による三重対角化： (4/3)N3
– 三重対角行列の固有値・固有ベクトルの計算： O(NM) ～ O(N3)
（解法および固有値の分布により大きく異なる。）
– 逆変換： 2N2 M
M≪N の場合は，三重対角化の演算量が支配的
本報告ではこの部分に絞って並列化・最適化の手法を述べる。
三重対角行列 T に対する固有値・固有ベクトル解法
解法
概要
計算量
並列化 一部の固有値
・固有ﾍﾞｸﾄﾙ
のみの計算
QR法
相似変換によりTを
対角行列に変換
O(N2)
+ O(N3)
容易
不可
DC法
（分割統治
法）
Tを再帰的に半分の
ｻｲｽﾞの行列に分割
O(N2)
～ O(N3)
容易
不可
二分法・
逆反復法
二分法で固有値，逆
反復法で固有ﾍﾞｸﾄﾙ
を計算
O(NM) 非自明
+ O(NM2)
Dhillon の
アルゴリズム
逆反復法の改良
（dqds，RRR，twisted
dqds = differential qd分解などの利用）
with shift
O(NM)
+ O(NM)
容易
RRR = Relatively Robust Representation; T をLLt の形で表現して計算を進める手法
可
可
２. 三重対角化処理の並列化
• ハウスホルダー法による三重対角化
• 共有メモリ向けの並列化
• サイクリック列分割による並列化
• Scattered Square 分割による並列化
• 性能評価
ハウスホルダー法による三重対角化
• 鏡像変換による列の消去
ベクトル
– 鏡像変換 H = I – a uut
• Hは対称な直交行列
• 与えられたベクトルの第１成分以外をゼロにする。
左からH
を乗算
• 第 k ステップでの処理
0
0
0
0
0
左からH
を乗算
0
右からH
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
各ステップでの処理の詳細
•鏡像変換ベクトルの作成
0
σ= sqrt(d td)
u = (d1 – sgn(d1)σ, d2, … , dN – k )
α= 2 /∥u∥2
0
C
d
•H による両側からの乗算
k
HCH = (I – a uu t)C (I – a uu t)
= C – a uu tC – aC uu t + a2 uu tC uu t
= C – up t – pu t + (1/2) a uu tpu t + (1/2) a up tuu t
( p ≡ aC u )
= C – uq t – qu t
( q ≡ p – (1/2) a (p tu) u )
•実際には，C の対称性を利用し， HCH は下三角部分のみ計算する。
•p ≡ aC u の計算でも，C の上三角部分は下三角部分で代用する。
ハウスホルダー法のアルゴリズム
[Step 1] k = 1からN –2まで以下の[Step 2] ～ [Step 8]を繰り返す。
[Step 2] 内積 σ(k) = sqrt(d (k)t d (k)) の計算
[Step 3] 鏡像変換ベクトル作成：
u(k) = (d (k)1 – sgn(d (k)1)σ(k), d (k)2, … , d (k)N – k )
[Step 4] 規格化定数の計算： α(k) = 2 /∥u(k)∥2
[Step 5] 行列ベクトル積： p(k) =α(k)C (k)u(k)
[Step 6] ベクトルの内積： β(k) =α(k) u(k)tp(k) /2
[Step 7] ベクトルの加算： q(k) = p(k) –β(k)u(k)
[Step 8] 行列のrank-2更新： C (k) = C (k) – u(k)q(k)t – q(k)u(k)t
0
0
d (k )1
C (k )
d (k )
k
ピボット列
ピボット行
ハウスホルダー法の特徴
• 演算量
– 行列サイズが N のときの演算量： 約 (4/3)N3
• 行列ベクトル積の演算量：
約 (2/3)N3
• rank-2更新の演算量：
約 (2/3)N3
• データの再利用性
– 演算はほとんどBLAS2のため，行列データの再利用性なし
– キャッシュマシンでは高い性能が期待できない。
ブロック化など，アルゴリズム上の工夫が必要
• 並列粒度
– BLAS2のため O(N2/p) であり，比較的高い。
共有メモリ向けの並列化
• 基本方針
– 行列ベクトル積 および rank-2更新 の部分について，BLAS2レベ
ルで並列化
– 並列BLASを使えば，並列化は極めて容易
• PU間でのアクセス競合の回避
– BLAS2はデータ再利用性がないため，
バスを通した共有メモリへのアクセスが
頻発
– アクセス競合による性能低下を回避
するには，後に述べるキャッシュ向け
最適化技法が有効
キャッシュ
PU0
PU1
PU2
バス
メモリ
PU3
分散メモリ向けの並列化 (1-1)
• ブロックサイクリック列分割による並列化
（Dongarra et al., 1992）
0
0
– 行列を幅 b の列ブロックに分割し，ブロックを
周期的にプロセッサに割り当てる。
C (k )
鏡像変換
ベクトル
k
• 各ステップの処理
(1) 鏡像変換ベクトル u(k) およびα(k) の作成
・第 k 列を担当するPUが実行
ﾌﾞﾛｰﾄﾞ
ｷｬｽﾄ
(2) u(k) およびα(k) のブロードキャスト
・第 k 列を担当するPUが他の全PUに送信
・二進木を使った場合，通信にかかる時間は O(N log p)
（ベクトル u(k) の長さは平均 N/2）
u (k )
分散メモリ向けの並列化 (1-2)
• 各ステップの処理（続き）
(3) C (k) の下三角部分と u(k) との積
・各PUは自分の行列要素を使って部分和ベクトルを計算
・結果の集計／ブロードキャストの通信時間はO(N log p)
(4) C (k) の上三角部分と u(k) との積
・各PUは自分の行列要素を使って結果ﾍﾞｸﾄﾙの一部を計算
・全対全通信により全PUに結果ﾍﾞｸﾄﾙの全要素を分配する
通信時間はO(N log p)
(5) (3)と(4)の結果の集計，β(k) の計算，q(k) の計算
×
＝
部分和の集計
／ﾌﾞﾛｰﾄﾞｷｬｽﾄ
× ＝
全対全通信により
結果ﾍﾞｸﾄﾙを分配
(6) rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t
・各PUで独立に計算
以上より，ブロックサイクリック分割の通信時間は O(N log p)
分散メモリ向けの並列化 (2-1)
• Scattered Square 分割による並列化
（Choi et al., 1995）
0
0
– 行列を b×b のブロックに分割し，ブロックに対して
大きさ √p× √ p のプロセッサテンプレートを周期的
に割り当てる。
鏡像変換
ベクトル
C (k )
k
• 各ステップの処理
(1) 鏡像変換ベクトル u(k) およびα(k) の作成
・第 k 列を担当するPU群が（通信しながら）実行
(2)
u(k)
およびα(k) のブロードキャスト
u (k )
・第 k 列担当のPUが，まず横方向のPUにﾌﾞﾛｰﾄﾞｷｬｽﾄ
・次に，対角ブロックを担当するPUが縦方向のPUに
ﾌﾞﾛｰﾄﾞｷｬｽﾄ
・通信にかかる時間は O(N log p /√p)
ﾌﾞﾛｰ
ﾄﾞｷｬｽ
ﾄ
ﾌﾞﾛｰﾄﾞ
ｷｬｽﾄ
分散メモリ向けの並列化 (2-2)
• 各ステップの処理（続き）
(3) C (k) の下三角部分と u(k) との積
×
＝
・各PUは自分の行列要素を使って部分和ベクトルを計算
・横方向に集計し，対角PUに集める（O(N log p /√p) ）。
部分和を横方向に集計し，
結果を対角PUに集める。
(4) C (k) の上三角部分と u(k) との積
・各PUは自分の行列要素を使って部分和ベクトルを計算
・縦方向に集計し，対角PUに集める（O(N log p /√p) ）。
(5) (3)と(4)の結果の集計とブロードキャスト
×
＝
部分和を縦方向に集計し，
結果を対角PUに集める。
・対角PUにおいて，(3)と(4)の結果を集計
・結果の部分ベクトルを，横方向のPUにブロードキャスト（O(N log p /√p) ）
・結果の部分ベクトルを，縦方向のPUにブロードキャスト（O(N log p /√p) ）
分散メモリ向けの並列化 (2-3)
• 各ステップの処理（続き）
(6) β(k) の計算，q(k) の計算
(7) rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t
・各PUで独立に計算
• 両並列化方式の比較
– 各PUの計算時間
• 各ステップの計算時間は両方式ともO(N2 / p) （p に反比例して減少）
– 通信時間
• ブロックサイクリック列分割は O(N log p) （pとともに増加）
• ブロックScattered Square 分割は O(N log p /√p) （pとともに減少）
PU台数が増えるほど， ブロックSS分割が有利
サイクリック列分割とSS分割の性能比較
10
5
性能予測モデルによる推定性能
SS分割
Mflops
10
4
ｻｲｸﾘｯｸ列分割
10
10
3
2
10
0
10
1
10
2
10
3
10
4
Num. of Proc.
超並列ではScattered Square分割が優位
SMPクラスタ上での性能評価
ブロックScattered Square 分割による並列化（MATRIX/MPP HZES1MDP使用）
SR8000/F1（SMPクラスタ）の１～16ノードで評価。キャッシュ向け最適化あり。
N=2000 ～ 32000 のエルミート行列での性能
100
90
Gflops
•
•
•
80
70
60
50
40
30
20
10
0
1000
1832[sec]
16-nodes
4-nodes
902[sec]
1-node
447[sec]
10000
Matrix size
100000
４. 三重対角化処理のキャッシュ向け最適化
• 従来のハウスホルダー法の問題点
• Bischof のアルゴリズム
• 性能・精度評価
従来のハウスホルダー法の問題点
• 演算量
– 行列サイズが N のときの演算量： 約 (4/3)N3
• 行列ベクトル積の演算量：
約 (2/3)N3
• rank-2更新の演算量：
約 (2/3)N3
• データの再利用性
– 演算はほとんどBLAS2のため，行列データの再利用性なし
– キャッシュマシンでは高い性能が期待できない。
ブロック化など，アルゴリズム上の工夫が必要
Bischof のアルゴリズム
• 原理（Bishof et al., 1993, 1994）
– 密行列 A をまず半帯幅Lの帯行列 B に変換し，次にこの帯行列
を三重対角行列 T に変換する。
次数 N
半帯幅 L
帯行列化
約 (4/3)N3
A
0
村田法
O(N2L)
0
B
0
0
T
• 三重対角化を２段階で行うことの利点
– 帯行列への変換は，BLAS3のみを使って実行可能。
– 帯行列から三重対角行列への変換の演算量はO(N2L)であり，
前半部に比べてずっと小さい。
Bischof のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – UαUt
• Hは直交行列
• 与えられたブロックベクトルの
第１ブロックを上三角行列にし，
それ以外をゼロにする。
左からH
を乗算
第kステップでの処理
0
0
左からH
を乗算
非ゼロ要素
0
0
0
ゼロにしたい部分
右からH
を乗算
0
影響を受ける部分
Bischof のアルゴリズム
[Step 1] K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。
[Step 2] D(K) を第１ブロックが上三角行列で第２ブロック以下が
ゼロ行列であるブロックベクトルに変換するブロック鏡像変換
I – U(K) (K)U(K)tを求める。
[Step 3] 行列・ブロックベクトル積： P(K) = C (K)U(K)α(K)
[Step 4] ブロックベクトルの内積： β(K) = α(K)tU(K)tP(K) / 2
[Step 5] ブロックベクトルの加算： Q(K) = P(K) –U(K)β(K)
[Step 6] 行列のrank-2L更新： C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t
すべてBLAS3（行列乗算）
Bischof のアルゴリズムの特徴と問題点
• データの再利用性
– 演算のうち，O(N3) の部分はすべて BLAS3 化
– キャッシュサイズに応じて L を選ぶことで，データ再利用性を最
大化することが可能
• 計算精度
– 多様な例題を用いて体系的に評価した結果は，まだ報告されて
いない。
• 演算量
– 後半部（村田法）の演算量は L に比例して増加
– 全体としての実行時間が最小になるよう，L を適切に選ぶ必要
がある。
４. 性能・精度評価
• 性能評価
– テスト例題
• N = 480，960，1920，3840 のフランク行列の三重対角化
• Bischof のアルゴリズムの性能を LAPACK と比較
– 評価マシン
• Opteron (1.6GHz)
• Alpha 21264A (750MHz)
• Power 5 (1.9GHz)
– 評価方法
• 演算量を(4/3)N3 とし，三重対角化に要した時間より各アルゴ
リズムのMFLOPS値を算出して比較
• ブロックサイズ L は，全実行時間が最小になるよう決定
• BLASルーチンとしては，全マシンで ATLAS 3.6.0 を使用
Opteron (1.6GHz) 上での性能
1800
Performance (GFLOPS)
1600
1400
LAPACK
Bishof
L=48
1200
1000
800
L’=32
600
400
200
0
480
960
1920
3840
Matrix size
Wu の方法は LAPACK に比べて約２倍の性能を達成
N = 3840 のとき，Bischof の方法はピークの50%以上の性能を達成
Alpha 21264A (750MHz) 上での性能
800
Performance (GFLOPS)
700
L=24, L’=4
LAPACK
Bishof
600
L=48
500
400
L’=32
300
200
100
0
480
960
1920
3840
Matrix size
Wu の方法は LAPACK に比べて約２倍の性能を達成
N = 3840 のとき，Bischof の方法はピークの50%以上の性能を達成
Power5 (1.9GHz) 上での性能
5.00E+00
Performance (GFLOPS)
4.50E+00
L=48, L’=2
LAPACK
Bishof
4.00E+00
ピークの
60%
L=96
3.50E+00
3.00E+00
2.50E+00
L’=64
2.00E+00
1.50E+00
1.00E+00
5.00E-01
0.00E+00
480
960
1920
3840
7680
11520
Matrix size
Wu の方法は LAPACK に比べて２倍以上の性能を達成
N = 11520 のとき，Bischof の方法はピークの60%の性能を達成
Power5 (1.9GHz) 上でのSMP性能
Performance (GFLOPS)
60
L=48, L’=2
N=480
N=960
N=1920
N=3840
N=7680
50
40
L=24, L’=4
30
L=24, L’=2
20
L=24, L’=2
10
L=12, L’=2
0
1
2
4
8
16
Number of processors
•帯行列化部分のみの性能（村田法は含まず）
•N = 7680 のとき， Bischof の方法は16PEで10倍の加速を達成
精度評価
• テスト例題
–
–
–
–
フランク行列 aij = min(i, j)
２次元ラプラス方程式の５点差分により得られる行列
Harwell-Boeing Sparse Matrix Collection の行列
乱数行列 （要素が [0, 1] の一様乱数である対称行列）
• 評価マシン
– Opteron (1.6GHz)
• 評価方法
– 各解法に対し，計算した固有値λi’の相対誤差の最大値
max 1≦i≦N |λi’ – λi| / |λi| を評価
– 比較対象の固有値λiとしては，解析解あるいはオリジナルのハウ
スホルダー法により求めた固有値を使用
フランク行列に対する計算精度
Matrix size
1.00E-06
480
960
1920
3840
Maximum relative error
1.00E-07
1.00E-08
1.00E-09
Dongarra
Bischof
Wu
1.00E-10
1.00E-11
1.00E-12
Bischof 及び Wu の方法の誤差は，Dongarra の方法と同程度
２次元５点差分の行列に対する計算精度
Matrix size
1.00E-10
480
960
1920
3840
Maximum relative error
1.00E-11
1.00E-12
1.00E-13
Dongarra
Bischof
Wu
1.00E-14
1.00E-15
1.00E-16
Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。
しかし，増大の程度は１桁以内に収まっている。
Harwell-Boeing Collection の行列に対する計算精度
Matrix name
1.00E-03
BCSSTK8
BCSSTK12
BCSSTK13
BCSSTK23
BCSSTK24
Maximum relative error
1.00E-04
1.00E-05
1.00E-06
Dongarra
Bischof
Wu
1.00E-07
1.00E-08
1.00E-09
Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。
しかし，増大の程度は２～３倍以内に収まっている。
乱数行列に対する計算精度
Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。
しかし，増大の程度は１桁以内に収まっている。
キャッシュ向け最適化に関するまとめ
• 性能・精度評価の結論
– キャッシュマシン向け三重対角化アルゴリズムとしては Wu のア
ルゴリズムが最高速であり，大規模問題では Dongarraのアルゴ
リズムの約２倍の性能が達成できる。
– 特に N = 3840の場合，Wu のアルゴリズムでは Opteron，Alpha
21264Aでピークの50%以上の性能が達成できる。
– Bischof / Wuのアルゴリズムによる固有値の相対誤差は，
Dongarraのアルゴリズムに比べてやや大きい。しかし，増大の程
度は多くの場合２～３倍以内，最大でも１桁程度である。
より詳細なデータについては，
http://www.na.cse.nagoya-u.ac.jp/~yamamoto/work.html を参照
５. 今後の展開
• より多様な例題での精度評価
• 固有ベクトル計算部分の実装と性能・精度評価
• 共有／分散メモリ型計算機上での並列化と性能評価
• 性能パラメータ L，L’ の自動最適化