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

Level-3 BLASに基づく特異値分解
アルゴリズムのSMP上での性能
名古屋大学 計算理工学専攻
山本有作
日本応用数理学会年会
2006年9月18日
目次
1. はじめに
2. 従来の特異値分解アルゴリズムとその問題点
3. Level-3 BLAS に基づく特異値分解アルゴリズム
4. 性能評価
5. まとめと今後の課題
1. はじめに
• 本研究で対象とする問題
– 実正方行列 A の特異値分解 A = US VT
• A : n×n 密行列
• S : n×n 対角行列
• U,V : n×n 直交行列
• 応用分野
–
–
–
–
統計計算 (主成分分析,最小2乗法)
信号処理 (独立成分分析など)
画像処理 (圧縮,ノイズ除去)
電子状態計算
本研究の目的
• 共有メモリ型並列計算機(SMP)上で高性能を達成できる
特異値分解ソルバを作成し,評価
• 背景
– 問題の大規模化
– CPUのマルチコア化などによる SMP 環境の普及
デュアルコア Xeon
Cell プロセッサ
(1+8 コア)
2. 従来の特異値分解アルゴリズムとその問題点
密行列 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アルゴリズム
逆変換
各部分の演算量と実行時間
密行列 A
二重対角化
演算量
実行時間(全特異ベクトル)
(8/3) n3
n = 5000,Xeon 2.8GHz(1~
4PU)
LAPACK での実行時間(秒)
600
逆変換
分割統治法
二重対角化
500
400
O(n2)
~ O(n3)
二重対角行列の
特異値・特異ベクトル計算
300
200
100
逆変換
4mn2
(左右 m 本ずつの
特異ベクトル)
Aの特異ベクトル {ui }, {vi }
0
1
2
4
・二重対角化が実行時間の
大部分を占める
・速度向上率が低い
二重対角化の性能が出ない原因
• 二重対角化の演算パターン
– 左右からのハウスホルダー変換による行・列の消去
• A(k) := (I – a w wT ) A(k)
• 演算は level-2 BLAS(行列ベクトル積と rank-1 更新)
– ただしブロック化により半分は level-3 BLAS にすることが可能
0
0
A(k)
0
右から
の変換
0
0
左から
の変換
0
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
k
• 演算パターンに関する問題点
– Level-2 BLAS はデータ再利用性が低い。
キャッシュの有効利用が困難であり,単体性能が出にくい。
プロセッサ間のアクセス競合により,並列性能向上も困難
3. Level-3 BLAS に基づく特異値分解アルゴリズム
• 2段階の二重対角化アルゴリズム(Bischof et al., ’93)
– 密行列 A をまず帯幅 L の下三角帯行列 C に変換
– 次にこの帯行列を下二重対角行列 B に変換
次数 n
下三角
帯行列化
約 (8/3)n3
A
0
村田法
の拡張
O(n2L)
0
C 帯幅 L
0
0
B
• 二重対角化を2段階で行うことの利点
– 前半の変換は,level-3 BLAS (行列乗算)のみを使って実行可能
キャッシュの有効利用が可能
– 後半の変換は level-2 BLAS が中心だが,演算量は O(n2L)
• 前半部に比べてずっと小さい。
下三角帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 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におけるメモリ競合の影響を低減可能
本アルゴリズムの長所と短所
• 長所
– Level-3 BLAS の利用により,二重対角化の性能を向上可能
• 同様のアイディアに基づく三重対角化アルゴリズムでは,高い単体
性能・並列性能を確認済み
• 短所
– 特異ベクトル計算のための計算量・記憶領域が増大
• 2段階の逆変換あるいは帯行列の特異値分解が必要
• 詳しくは次のスライドで説明
– 二重対角化の高速化効果が大きければ,計算量増大を考慮し
ても全体としては高速化できると予想
– 特に,求める特異ベクトルが少ない場合は効果が大きいはず。
特異ベクトルの計算手法
• 方法1: 下三角帯行列の特異ベクトルを直接計算
n
L
0
0
0
0
A
A の特異ベクトル
{ui }{vi }
C
4mn2
C の特異ベクトル
{zi }{wi }
帯行列用
逆反復法
O(mnL2+ m2n)
B
特異値
{σi }
QR法
dqds法
mdLVs法
二分法
• 長所
– 固有ベクトルの逆変換は1段階のみ
– 逆変換の演算量は 4mn2 (従来法と同じ)
• 短所
– 特異ベクトル計算のための実用的な手法は帯行列用逆反復法のみ
– 直交化が必要であり,演算量は O(mnL2+m2n)
特異ベクトルの計算手法(続き)
• 方法2: 二重対角行列の特異ベクトルを計算して2回逆変換
n
L
0
0
0
0
A
C
A の特異ベクトル
{ui }{vi }
4mn2
C の特異ベクトル
{zi }{wi }
B
4mn2
B の特異ベクトル
{xi }{yi }
特異値
{σi }
QR法
DC法
MR3
I-SVD
• 長所
– 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能
• 短所
– 逆変換の演算量が 8mn2 (従来法の2倍)。ただし level-3 化可能
– 村田法の変換をすべて記憶するため,n2 の記憶領域が余計に必要
SMP 上での level-3 BLAS の高速性を鑑み,方法2を採用
アルゴリズムの全体像
• 2段階の二重対角化と2段階の逆変換
• 二重対角行列の特異値分解には分割統治法を使用
分割統治法
(LAPACK DBDSDC)
0
(8/3)n3
level-3
A
A の特異ベクトル
{ui }{vi }
O(n2L)
level-2
0
C
4mn2
level-3
C の特異ベクトル
{zi }{wi }
0
0
B
4mn2
level-3
O(n2) ~ O(n3)
level-3
B の特異値 {σi }
特異ベクトル {xi }{yi }
• 特徴
– 演算量が O(n3) となる部分はすべて level-3 BLAS で実行可能
– SMP 向け並列化は,基本的に並列版 level-3 BLAS の使用により実現
– 村田法は OpenMP によるパイプライン型の並列化
村田法の並列化
• パイプライン型の並列化
– 第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. 性能評価
• 評価環境
– Xeon (2.8GHz), 1~4PU
• Linux + Intel Fortran ver. 8.1
• BLAS: Intel Math Kernel Library
• LAPACK: Intel Math Kernel Library
• ピーク性能: 5.6GFLOPS/CPU
– 富士通 PrimePower HPC2500 (2.0GHz), 1~32PU
• 富士通 Fortran
• BLAS: 富士通並列化版 BLAS
• LAPACK: 富士通並列化版 LAPACK
• ピーク性能: 8GFLOPS/CPU
• 評価対象・条件
– Level-3 BLAS に基づくアルゴリズムと LAPACK の性能を比較
– n = 1200 ~ 20000 の乱数行列の特異値分解(全特異ベクトルを計算)
• Level-3 アルゴリズムにとっては一番不利な条件
– Level-3 アルゴリズムの L(半帯幅)は各 n ごとに最適値を使用
Xeon での実行時間
• プロセッサ数を変えたときの実行時間
16
n = 1200
120
level-3
LAPACK
14
12
8
6
4
800
level-3
LAPACK
100
実
行
時
間
(
秒
)
10
n = 2500
0
600
80
500
60
400
300
40
200
100
0
1
2
4
level-3
LAPACK
700
20
2
n = 5000
0
1
2
4
PU数
1
2
4
• 結果
– Level-3 アルゴリズムでは PU 数に応じて実行時間が順調に減少
• 4PUでの加速率は 3~3.2 倍
– 4PU の場合は level-3 アルゴリズムが従来法より高速
HPC2500 での実行時間
• プロセッサ数を変えたときの実行時間
800
n = 5000
5000
n = 10000
14000
4500
700
level-3
LAPACK
600
500
400
300
200
実
行
時
間
(
秒
)
12000
level-3
LAPACK
4000
n = 20000
3500
10000
3000
8000
level-3
LAPACK
2500
3.5倍
6000
2000
1500
4000
1000
100
500
0
0
1
2
4
8
16
32
2000
0
1
2
4
8
16
32
PU数
1
2
4
8
16
32
• 結果
– Level-3 アルゴリズムは従来法に比べて最大3.5倍高速
– プロセッサ数が多いとき加速率が鈍るのは,非並列化部分(ブロッ
ク鏡像変換の作成など)の影響と思われる。
両手法の実行時間の内訳
• Xeon,n=5000の場合
800
600
逆変換
分割統治法
二重対角化
500
逆変換2
逆変換1
分割統治法
村田法
帯行列化
700
600
400
500
300
400
300
200
200
100
100
0
0
1
2
4
1
2
4
• 考察
– Level-3 アルゴリズムでは,どの部分の実行時間も順調に減少
– 逆変換1(村田法の逆変換)の占める時間が大きい。
この部分について,さらに高速化が必要
– 必要な特異ベクトルの本数が少ない場合,level-3 アルゴリズムは
さらに有利
両手法の実行時間の内訳
• HPC2500,n=10,000の場合
5000
5000
4500
4500
4000
4000
逆変換
分割統治法
二重対角化
3500
3000
逆変換2
逆変換1
分割統治法
村田法
帯行列化
3500
3000
2500
2500
2000
2000
1500
1500
1000
1000
500
500
0
0
1
2
4
8
16
32
1
2
4
8
16
32
• 考察
– Level-3 アルゴリズムでは,どの部分の実行時間も順調に減少
– 従来法は,二重対角化の部分の加速が鈍い。
• ただし,32PUで6倍程度は加速
• メモリバンド幅が大きいためと思われる。