Transcript Level-3 BLASに基づく特異値分解アルゴリズムのSMP上での性能

Level-3 BLASに基づく特異値分解
アルゴリズムのSMP上での性能
名古屋大学 計算理工学専攻
山本有作
日本応用数理学会年会
2006年9月18日
目次
１. はじめに
２. 従来の特異値分解アルゴリズムとその問題点
３. Level-3 BLAS に基づく特異値分解アルゴリズム
４. 性能評価
５. まとめと今後の課題
１. はじめに
• 本研究で対象とする問題
– 実正方行列 A の特異値分解 A = US VT
• A ： n×n 密行列
• S ： n×n 対角行列
• U，V ： n×n 直交行列
• 応用分野
–
–
–
–
統計計算 （主成分分析，最小2乗法）
信号処理 （独立成分分析など）
画像処理 （圧縮，ノイズ除去）
電子状態計算
本研究の目的
• 共有メモリ型並列計算機（SMP）上で高性能を達成できる
特異値分解ソルバを作成し，評価
• 背景
– 問題の大規模化
– CPUのマルチコア化などによる SMP 環境の普及
デュアルコア Xeon
Cell プロセッサ
（1+8 コア）
２. 従来の特異値分解アルゴリズムとその問題点
密行列 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 はデータ再利用性が低い。
キャッシュの有効利用が困難であり，単体性能が出にくい。
プロセッサ間のアクセス競合により，並列性能向上も困難
３. Level-3 BLAS に基づく特異値分解アルゴリズム
• ２段階の二重対角化アルゴリズム（Bischof et al., ’93）
– 密行列 A をまず帯幅 L の下三角帯行列 C に変換
– 次にこの帯行列を下二重対角行列 B に変換
次数 n
下三角
帯行列化
約 (8/3)n3
A
0
村田法
の拡張
O(n2L)
0
C 帯幅 L
0
0
B
• 二重対角化を２段階で行うことの利点
– 前半の変換は，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 の利用により，二重対角化の性能を向上可能
• 同様のアイディアに基づく三重対角化アルゴリズムでは，高い単体
性能・並列性能を確認済み
• 短所
– 特異ベクトル計算のための計算量・記憶領域が増大
• ２段階の逆変換あるいは帯行列の特異値分解が必要
• 詳しくは次のスライドで説明
– 二重対角化の高速化効果が大きければ，計算量増大を考慮し
ても全体としては高速化できると予想
– 特に，求める特異ベクトルが少ない場合は効果が大きいはず。
特異ベクトルの計算手法
• 方法１： 下三角帯行列の特異ベクトルを直接計算
n
L
0
0
0
0
A
A の特異ﾍﾞｸﾄﾙ
{ui }{vi }
C
4mn2
C の特異ﾍﾞｸﾄﾙ
{zi }{wi }
帯行列用
逆反復法
O(mnL2+ m2n)
B
特異値
{σi }
QR法
dqds法
mdLVs法
二分法
• 長所
– 固有ベクトルの逆変換は１段階のみ
– 逆変換の演算量は 4mn2 （従来法と同じ）
• 短所
– 特異ベクトル計算のための実用的な手法は帯行列用逆反復法のみ
– 直交化が必要であり，演算量は O(mnL2+m2n)
特異ベクトルの計算手法（続き）
• 方法２： 二重対角行列の特異ベクトルを計算して２回逆変換
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 の高速性を鑑み，方法２を採用
アルゴリズムの全体像
• ２段階の二重対角化と２段階の逆変換
• 二重対角行列の特異値分解には分割統治法を使用
分割統治法
（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 ステップ ・・・
が同時に実行可能
４. 性能評価
• 評価環境
– 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倍程度は加速
• メモリバンド幅が大きいためと思われる。