K - 計算基盤研究室

Download Report

Transcript K - 計算基盤研究室 

応用数理工学特論
線形計算と
ハイパフォーマンスコンピューティング
2008年5月29日
計算理工学専攻 張研究室
山本有作
目次
１. 固有値計算の概要
２. 対称行列の固有値計算のための前処理
・ ハウスホルダー法による三重対角化
・ 三重対角化処理の並列化
・ 三重対角化処理のキャッシュ向け最適化
３. 非対称行列の固有値計算のための前処理
・ ハウスホルダー法によるヘッセンベルグ化
１. 固有値計算の概要
• 固有値計算とその応用
• 固有値・固有ベクトル計算の流れ
• 各部分の演算量
• 三重対角行列 T に対する固有値・固有ベクトル解法
• ヘッセンベルグ行列 H に対する固有値・固有ベクトル解法
１. 固有値計算の概要
• 対象とする問題
– 標準固有値問題 Ax = λx
– λ： 固有値， x ： 固有ベクトル
– A が n×n 行列の場合，固有値は全部で n 個存在
• 固有値問題の分類
– A が対称（エルミート）行列 ⇔ A が非対称（非エルミート）行列
– A が密行列 ⇔ A が疎行列
– 全固有値・固有ベクトルを計算 ⇔ 一部のみを計算
固有値問題の応用
分野
固有値
固有ベクトル
電子状態計算
分子軌道法
エネルギーレベル
電子軌道
構造解析
固有振動数
振動モード
主成分分析
各成分の寄与度
主成分
理論流体力学
（乱流のｽｹｰﾘﾝｸﾞ指
数）
スケーリング指数
相関関数
固有値計算の難しさ
• 固有値の特徴付け
– λが A の固有値 ⇔ det(A –λI) = 0
– λは n 次代数方程式の解
• 固有値問題と代数方程式
– n 次代数方程式が解ければ，固有値λは求まる（数学的には）。
– 逆に，任意の n 次代数方程式は， n×n 行列の固有値問題として
定式化できる。
• 固有値問題の難しさ
– 任意の n 次代数方程式を有限回の演算で解くアルゴリズムは存
在しない。
– よって，任意の固有値問題を有限回の演算で解くアルゴリズムは
存在しない。
求解には反復法が必須
密行列の固有値計算の基本的な考え方
• 密行列に対して直接に反復法を適用すると，演算コストが大きい。
• 有限回の演算により，密行列をまず中間的な行列に相似変換
– 相似変換により，固有値は不変
– 実際には，相似変換の中でも，特に直交変換を利用。
– 中間的な行列とは，三重対角行列あるいはヘッセンベルグ行列
• 中間的な行列に対して反復法を適用し，固有値・固有ベクトルを計算
• 直交変換を逆に施すことにより，固有ベクトルを元の行列の固有ベク
トルに変換
対称行列の固有値・固有ベクトル計算の流れ
実対称行列 A
計算内容
計算手法
三重対角化
Q*AQ = T
(Q: 直交行列)
ハウスホルダー法
三重対角行列 T
三重対角行列の
固有値・固有ﾍﾞｸﾄﾙ計算
Tyi =λi yi
Tの固有値 {λi }，
固有ベクトル {yi }
逆変換
Aの固有ベクトル {xi }
xi = Qyi
QR法
DC法
二分法・逆反復法
MR3アルゴリズム
逆変換
各部分の演算量
• 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)
逆反復法の改良
（dqds，RRR，twisted
dqds = differential qd分解などの利用）
with shift
MR3
アルゴリズム
O(NM)
+ O(NM)
容易
RRR = Relatively Robust Representation; T をLLt の形で表現して計算を進める手法
可
可
非対称行列の固有値・固有ベクトル計算の流れ
計算内容
計算手法
Q*AQ = H
(Q: 直交行列)
ハウスホルダー法
実非対称行列 A
ヘッセンベルグ化
ﾍｯｾﾝﾍﾞﾙｸﾞ行列 H
ヘッセンベルグ行列の
固有値・固有ﾍﾞｸﾄﾙ計算
Hyi =λi yi
Hの固有値 {λi }，
固有ベクトル {yi }
逆変換
Aの固有ベクトル {xi }
xi = Qyi
QR法
DC法
ホモトピー法
複素周回積分法
逆変換
各部分の演算量
• N×N の実対称行列 A に対し，M 組の固有値・固有ベク
トルを求める場合
– ハウスホルダー法によるヘッセンベルグ化： (10/3)N3
– QR法によるヘッセンベルグ行列の固有値・固有ベクトルの計算：
10N3 程度
– 逆変換： 2N2 M
ﾍｯｾﾝﾍﾞﾙｸﾞ行列 H に対する固有値・固有ベクトル解法
解法
概要
計算量
並列化 一部の固有値
・固有ﾍﾞｸﾄﾙ
のみの計算
QR法
相似変換によりHを対
角行列に変換
10N3
程度
可能
不可
複素周回
積分法
ｺｰｼｰの定理を利用し
1 / det(A–λI) の極を
求める
O(N2M)
容易
可
２. 対称行列の固有値計算のための前処理
• ハウスホルダー法による三重対角化
• 三重対角化処理の並列化
– 共有メモリ向け
– 分散メモリ向け
• 三重対角化処理のキャッシュ向け最適化
ハウスホルダー法による三重対角化
• 鏡像変換による列の消去
ベクトル
– 鏡像変換 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)
• 双方向サイクリック 分割による並列化
（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台数が増えるほど， 双方向サイクリック分割が有利
サイクリック列分割と双方向サイクリック分割の性能比較
10
5
性能予測モデルによる推定性能
双方向分割
Mflops
10
4
ｻｲｸﾘｯｸ列分割
10
10
3
2
10
0
10
1
10
2
10
3
10
4
Num. of Proc.
超並列では双方向サイクリック分割が優位
SMPクラスタ上での性能評価
双方向サイクリック 分割による並列化
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
三重対角化処理のキャッシュ向け最適化
• 従来のハウスホルダー法の問題点
• Dongarra のアルゴリズム
• Bischof のアルゴリズム
• 性能・精度評価
従来のハウスホルダー法の問題点
• 演算量
– 行列サイズが N のときの演算量： 約 (4/3)N3
• 行列ベクトル積の演算量：
約 (2/3)N3
• rank-2更新の演算量：
約 (2/3)N3
• データの再利用性
– 演算はほとんどBLAS2のため，行列データの再利用性なし
– キャッシュマシンでは高い性能が期待できない。
ブロック化など，アルゴリズム上の工夫が必要
Dongarra のアルゴリズム
• 原理（Dongarra et al., 1992）
– 複数本のピボット列・行を溜めておき，後でまとめて更新
– 全演算量（(4/3)n3）の半分を行列乗算（BLAS3）にでき，SMPで有効
– LAPACK，ScaLAPACK 等で採用され，広く使われる。
0
0
=
–
×
–
×
複数段分の更新（BLAS3使用）
• 問題点
– 演算の残り半分は，依然として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 を適切に選ぶ必要
がある。
– 逆変換も２段階となるため，逆変換の演算量が 4N2M に倍増
性能評価
• Dongarra のアルゴリズム（LAPACK版）と Bischof のアルゴリズムの
性能を，以下の条件で比較。
– 計算機環境
• Xeon (2.8GHz)
• 富士通 PrimePower HPC2500
– プロセッサ台数
• 1 ～ 32 （計算機により異なる）
– 問題サイズ
• n = 1500, 3000, 6000, 12000, 24000
• 全固有値・固有ベクトルを計算する場合と，固有値のみ計算の場合
– アルゴリズム中のパラメータ
• 両アルゴリズムとも，実験により得られた最適値を採用
– 並列化
• 並列版 BLAS を利用
Xeon (2.8GHz) での性能
実行時間
(sec.)
n = 1500
12
Bischof
Dongarra
10
8
6
4
2
0
1
2
n = 3000
80
70
60
50
40
30
20
10
0
4
n = 6000
600
500
400
300
200
100
プロセッサ数
0
1
2
1
4
2
4
全固有値・固有ベクトルを計算する場合
実行時間
(sec.)
3
n = 1500
2.5
2
1.5
1
Bischof
Dongarra
0.5
0
1
2
n = 3000
20
18
16
14
12
10
8
6
4
2
0
4
1
2
n = 6000
180
160
140
120
100
80
60
40
20
0
4
プロセッサ数
1
2
4
全固有値のみを計算する場合
固有値のみを計算する場合は，Bischof のアルゴリズムが有利
富士通 PrimePower HPC2500 での性能
実行時間
(sec.)
n = 6000
450
400
350
300
250
200
150
100
50
0
n = 12000
3500
Bischof
Dongarra
3000
25000
2500
20000
2000
15000
1500
10000
1000
1
2
4
500
5000
0
0
8 16 32
n = 24000
30000
1
2
4
8 16 32
プロセッサ数
1
2
4
8 16 32
全固有値・固有ベクトルを計算する場合
実行時間
(sec.)
n = 6000
250
Bischof
Dongarra
200
n = 12000
1600
1400
12000
1200
10000
1000
150
8000
800
100
6000
600
50
0
1
2
4
8 16 32
n = 24000
14000
400
4000
200
2000
0
0
1
2
4
8 16 32
プロセッサ数
1
2
4
8 16 32
全固有値のみを計算する場合
プロセッサ数が4以上ならば，全固有値・固有ベクトルを計算する場合でも
Bischof のアルゴリズムが有利
プロセッサ数と必要な固有ベクトル数に応じた
最適アルゴリズム
• Xeon (2.8GHz)
必要な固有
ベクトルの
割合
100%
80%
60%
n = 1500
100%
Dongarra
40%
20%
Bischof
0%
1
2
4
n = 3000
100%
80%
80%
60%
60%
40%
40%
20%
20%
0%
0%
1
2
4
n = 6000
プロセッサ数
1
2
4
• 富士通 PrimePower HPC2500
必要な
固有
ベクトル
の割合
100%
80%
n = 6000
100%
n = 12000
100%
80%
80%
60%
60%
40%
40%
20%
20%
20%
0%
0%
0%
Dongarra
60%
40%
Bischof
1 2 4 8 16 32
1 2 4 8 16 32
n = 24000
1 2 4 8 16 32
プロセッサ
数
精度評価
• テスト例題
–
–
–
–
フランク行列 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 の方法に比べやや大きい。
しかし，増大の程度は１桁以内に収まっている。
キャッシュ向け最適化に関するまとめ
• キャッシュ向け最適化を行った三重対角化手法として Dongarra と
Bischof のアルゴリズムを取り上げ，計算機環境と問題サイズを変え
て性能比較を行った。
• プロセッサ数が増えるに連れ，また，必要な固有ベクトルの割合が
少なくなるに連れ，Bischof のアルゴリズムが有利となる傾向が見ら
れた。
• Bischof のアルゴリズムによる固有値の相対誤差は，Dongarraのア
ルゴリズムに比べてやや大きい。しかし，増大の程度は多くの場合２
～３倍以内，最大でも１桁程度である。
より詳細なデータについては，
http://www.na.cse.nagoya-u.ac.jp/~yamamoto/work.html を参照
グループ発表について
• 概要
– 最後の２回の授業を使って行う。
– ２～３人のグループで１つの発表を行う。
– 講義で勉強したことをもとに，自分たちで数値実験を行った結果
を発表。あるいは論文を読んでまとめて発表。
• スケジュール
– グループ決定： 6/19 まで
– 発表テーマ決定： 6/26 まで
– 発表： 7/10，7/17 （各グループ15～20分程度）
発表テーマの例
• 自分の取り組んでいるシミュレーションのプログラムを，
高性能化あるいは並列化した結果の報告
• 講義で学んだ行列計算アルゴリズムの性能評価
– 分散メモリ型計算機上でのLU分解の並列化
– スパースソルバの性能評価
– 固有値計算アルゴリズムの性能評価，など
• 高性能計算に関する論文を読み，まとめて報告
• その他