分散メモリ型スパースソルバの開発と評価

Level-3 BLASに基づく二重対角化
アルゴリズムとその性能評価
名古屋大学 計算理工学専攻
山本有作
2006年1月5日-7日
第3回計算数学研究会
目次
1. はじめに
2. 従来の二重対角化アルゴリズム
3. Level-3 BLAS に基づく二重対角化アルゴリズム
4. 性能評価
5. 特異ベクトルの計算手法
6. まとめと今後の課題
1. はじめに
• 本研究で対象とする問題
– 実正方行列 A の下二重対角行列への変換 B = U0TAV0
• A : n×n 密行列
• B : n×n 下二重対角行列
• U0,V0 : n×n 直交行列
• 応用分野
– 行列 A の特異値分解に対する前処理
• 画像処理の分野では,n = 10,000 以上の大規模行列に対す
る特異値分解の需要あり
• 必要な特異ベクトルは数本程度の場合が多い
• 前処理を数分程度で行えることが望ましい
特異値分解の流れ
密行列 A
計算内容
計算手法
二重対角化
U0TAV0 = B
ハウスホルダー法
(U0, V0: 直交行列)
二重対角行列 B
二重対角行列の
特異値・特異ベクトル計算
Bvi =σi xi
BTxi =σi yi
Bの特異値 {σi },
特異ベクトル {xi }{yi }
vi = V0 yi
逆変換
ui = U 0 x i
Aの特異ベクトル {ui }, {vi }
QR法
分割統治法
MR3アルゴリズム
I-SVDアルゴリズム
逆変換
各部分の演算量
• m 組の特異値・特異ベクトルを求める場合の演算量
–
–
–
–
二重対角化:
特異値の計算:
特異ベクトルの計算:
逆変換:
(8/3) n3
O(mn) ~ O(n2)
O(mn) ~ O(m2n)
O(mn2)
– m ≪ n の場合は演算量の大部分を二重対角化が占める。
– ハウスホルダー法の高速化が必要
2. 従来の二重対角化アルゴリズム
• ハウスホルダー法による二重対角化
ベクトル
– 鏡像変換 H = I – a wwT による列の消去
• Hは対称な直交行列
• 与えられたベクトルの第1成分以外をゼロにする。
左からH
を乗算
• 第 k ステップでの処理
0
0
0
0
0
右からHkR
を乗算
0
左からHkL
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
ハウスホルダー法のアルゴリズム
[Step 1] k = 1 から n –1 まで以下の[Step 2] ~ [Step 8]を繰り返す。
[Step 2] A(k, k+1:n) を 0 にする鏡像変換 HkR = I – akR wkR (wkR)T の計
算
[Step 3] 行列ベクトル積: p := akR A(k:n, k:n) wkR
[Step 4] 行列のrank-1更新: A(k:n, k:n) := A(k:n, k:n) – p(wkR)T
[Step 5] A(k+2:n, k) を 0 にする鏡像変換 HkL = I – akL wkL (wkL)T の計算
[Step 6] 行列ベクトル積: qT := akL (wkL)T A(k+1:n, k:n)
[Step 7] 行列のrank-1更新: A(k +1:n, k:n) := A(k +1:n, k:n) – wkLqT
0
0
A(k, k+1:n)
A(k:n, k:n)
A(k +2:n, k)
A(k +1:n, k :n)
k
0
0
ハウスホルダー法の特徴と問題点
• 特徴
– 全演算量のほとんどが,level-2 BLAS である行列ベクトル積と
行列のrank-1更新で占められる。
• 全演算量:
約 (8/3)n3
• 行列ベクトル積の演算量:
約 (4/3)n3
• rank-1更新の演算量:
約 (4/3)n3
• 問題点
– level-2 BLASのため,行列データの再利用性なし
• キャッシュミス,SMPでのメモリ競合の影響大
– 最近のマイクロプロセッサ,SMP上での高性能が期待できない。
3. Level-3 BLAS に基づく二重対角化アルゴリズム
• 基本的なアイディア
– 密行列 A をまず帯幅 L の下三角帯行列 C に変換
– 次にこの帯行列を下二重対角行列 B に変換
次数 n
下三角
帯行列化
約 (8/3)n3
A
0
村田法
の拡張
O(n2L)
0
C 帯幅 L
0
0
B
• 二重対角化を2段階で行うことの利点
– 下三角帯行列への変換は, level-3 BLAS のみを使って実行可能
– 下三角帯行列から二重対角行列への変換の演算量は O(n2L) で
あり,前半部に比べてずっと小さい。
(参考)対称行列の三重対角化の場合
• Bischof のアルゴリズム(Bishof et al., 1993, 1994)
– 対称密行列 A をまず半帯幅 L の対称帯行列 C に変換
– 次にこの帯行列を三重対角行列 T に変換
次数 n
半帯幅 L
帯行列化
約 (4/3)n3
A
0
村田法
O(n2L)
0
C
0
0
T
• Bischof のアルゴリズムの性能・精度(山本, 2005)
– 固有値の誤差は,LAPACKで使われる Dongarra のアルゴリズム
に比べ,多くの場合2~3倍程度
– 速度は2~3倍高速
– 様々なマイクロプロセッサ上で,ピークの50~70%の性能を達成
下三角帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – WαWT
• Hは直交行列
• 与えられたブロックベクトルを上三角
左からH
行列(正確には右上三角部分のみ
を乗算
非零でそれ以外が零の行列)に変形
第 K ステップでの処理
0
0
0
0
右からHK
を乗算
非ゼロ要素
R
0
ゼロにしたい部分
左からHKL
を乗算
0
影響を受ける部分
下三角帯行列化のアルゴリズム(続き)
[Step 1] K = 1からN /L–1まで以下の[Step 2] ~ [Step 6]を繰り返す。
[Step 2] A(K, K:N) を上三角行列に変形する鏡像変換
HKR = I – WKR aKR (WKR)T の計算
[Step 3] 行列・ブロックベクトル積: P := A(K:N, K:N) WKR aKR
[Step 4] 行列のrank-L更新: A(K:N, K:N) := A(K:N, K:N) – P(WKR)T
[Step 5] A(K+1:N, K) を上三角行列に変形する鏡像変換
HKL = I – WKL aKL (WKL)T の計算
[Step 6] 行列・ブロックベクトル積: QT := aKL (WKL)T A(K+1:N, K:N)
[Step 7] 行列のrank-L更新:
A(K+1:N, K:N) := A(K+1:N, K:N) – WkLQT
• 本アルゴリズムの特徴
すべて level-3 BLAS(行列乗算)
– 演算が level-3 BLAS 中心のため,キャッシュの有効利用が可能
– SMPにおけるメモリ競合の影響を低減可能
下三角帯行列の二重対角化
• 対称帯行列の三重対角行列への相似変換(村田法)
– 第 k 列の三重対角化
• 第 k 列の副対角要素より下を 0 にする鏡像変換を左と右からかける。
– Bulge-chasing
• 上記の結果,帯より下に非零要素がはみ出すので,相似変換を繰り
返して,はみ出した分を右下に移動し,右下隅から追い出す。
• 村田法のアイディアの二重対角化への適用(Lang, 1996)
– 第 k 列の二重対角化
• 第 k 列の副対角要素より下を 0 にする鏡像変換を左からかける。
– Bulge-chasing
• 上記の結果,上三角部分に非零要素がはみ出すので,左右から直
交変換を繰り返し掛けることにより,はみ出した分を右下に移動し,
右下隅から追い出す。
第1列の二重対角化と bulge-chasing
左側からの
直交変換で
更新された
要素
左側からの
直交変換で
更新された
要素
第2列の二重対角化と bulge-chasing
左側からの
直交変換で
更新された
要素
左側からの
直交変換で
更新された
要素
演算量は 8n2L
SMP向けの並列化
• 下三角帯行列化
– 並列版の Level-3 BLAS を使えば,逐次版のプログラムをその
まま SMP 上で並列化可能
– OpenMP などを用いて手動で並列化すれば,並列化効率を更に
向上可能
SMP向けの並列化(続き)
• 二重対角化向け村田法
– 第1列の二重対角化処理と第2列の二重対角化処理の並列性
第1列のbulge-chasing
における,右側からの
第3の直交変換で更新
される要素
第2列のbulge-chasing
における,右側からの
第1の直交変換で更新
される要素
第1列による二重対角化は,今後 より右の要素にのみ影響を及ぼす。
第1列の計算が右下まで行くのを待たずに,第2列の計算を開始できる。
– 一般の場合の並列性
• 第1列に対する bulge-chasing の第 k ステップ
• 第2列に対する bulge-chasing の第 k–2 ステップ
• 第3列に対する bulge-chasing の第 k–4 ステップ ・・・
が同時に実行可能
4. 性能評価
• 評価環境
– Opteron (1.8GHz), 1~4CPU
• Linux 上で PGI Fortran + GOTO BLAS を使用
• ピーク性能: 3.6GFLOPS/CPU
– PowerPC G5 (2.0GHz), 1~2CPU
• Mac OS X 上で IBM XL Fortran + GOTO BLAS を使用
• ピーク性能: 8GFLOPS/CPU
• 評価対象
– n = 2500 ~ 12500 の乱数行列の二重対角化
– 下三角帯行列化部分は並列版 GOTO BLAS,村田法部分は OpenMP
により並列化
• 評価方法
– 下三角帯行列化 + 村田法の合計時間を測定し,演算量を(8/3)n3 として
MFLOPS 値を算出
– ブロックサイズ L は,各 n に対して最適な値を使用
Opteron (1.8GHz) 1CPU上での性能
Opteron 1CPU
3500
L=100
Performance (GFLOPS)
3000
ピークの
80%
2500
2000
Opteron 1CPU
1500
1000
500
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき,本アルゴリズムはピークの80%の性能を達成
Opteron (1.8GHz) 4CPU上での性能
Performance (MFLOPS)
12000
Opteron 1CPU
Opteron 4CPUs
10000
L=100
ピークの
67%
286s
8000
6000
3.3倍
4000
L=100
2000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき,本アルゴリズムは4CPUで3.3倍の加速率
ピークの67%の性能を達成
各部分の演算時間(Opteron, n=12500)
2500
Execution time (s)
2000
村田法
下三角帯行列化
1500
1000
500
0
L=25
1CPU
L=50
1CPU
L=100
1CPU
L=25
4CPU
L=50
4CPU
L=100
4CPU
大規模問題に対しては,村田法の占める時間の割合は小さい。
村田法の加速率は,L=100 のとき4CPUで3.3倍
PowerPC G5 (2.0GHz) 1CPU上での性能
PowerPC G5 1CPU
8000
Performance (GFLOPS)
7000
PowerPC G5 1CPU
6000
L=100
5000
ピークの
57%
4000
3000
2000
1000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき,本アルゴリズムはピークの57%の性能を達成
PowerPC G5 (2.0GHz) 2CPU上での性能
8000
Performance (MFLOPS)
L=100
PowerPC G5 1CPU
Power PC G5 2CPUs
7000
6000
421s
5000
ピークの
42%
L=100
4000
1.5倍
3000
2000
1000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき,本アルゴリズムは2CPUで1.5倍の加速率
ピークの42%の性能を達成
各部分の演算時間(PowerPC G5, n=12500)
1800
1600
村田法
下三角帯行列化
Execution time (s)
1400
1200
1000
800
600
400
200
0
L=25
1CPU
L=50
1CPU
L=100
1CPU
L=25
2CPU
L=50
2CPU
L=100
2CPU
大規模問題に対しては,村田法の占める時間の割合は小さい。
村田法の加速率は,L=100 のとき2CPUで1.4倍
5. 特異ベクトルの計算手法
• 方法1: 二重対角行列の特異ベクトルを計算して2回逆変換
n
L
0
0
0
0
A
A の特異ベクトル
{ui }{vi }
C
2mn2
C の特異ベクトル
{zi }{wi }
B
2mn2
B の特異ベクトル
{xi }{yi }
特異値
{σi }
QR法
DC法
MR3
I-SVD
• 長所
– 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能
• 短所
– 逆変換の演算量が 4mn2 (ただし,m が小さければ影響は少ない)
– 村田法の変換をすべて記憶するため,n2 の記憶領域が余計に必要
5. 特異ベクトルの計算手法
• 方法2: 下三角帯行列の特異ベクトルを直接計算
n
L
0
0
0
A
A の特異ベクトル
{ui }{vi }
C
2mn2
0
C の特異ベクトル
{zi }{wi }
帯行列用
逆反復法
O(mnL2+ m2n)
B
特異値
{σi }
QR法
dqds法
mdLVs法
二分法
• 長所
– 村田法の逆変換のための演算量(2mn2)・記憶領域(n2)が不要
• 短所
– 帯行列用逆反復法(O(mnL2))と特異ベクトルの直交化(O(m2n))が必要
– 特異値σi は C の高精度な固有値でないため,ツイスト分解は使用不可
5. 特異ベクトルの計算手法
• 方法3: 村田法を使わず帯行列 C の特異値・特異ベクトルを直接計算
n
L
0
帯行列版
mdLVs法
特異値
{σi }
0
A
A の特異ベクトル
{ui }{vi }
C
2mn2
帯行列版
ツイスト分解
C の特異ベクトル
{zi }{wi }
• 長所
– 村田法の逆変換のための演算量(2mn2)・記憶領域(n2)が不要
– うまく行けば,直交化が不要なアルゴリズムを作れる可能性あり
• 問題点
– そもそも帯行列版の mdLVs 法,ツイスト分解は構成できるか?
– 演算量,演算の安定性は?