Transcript 長方行列向け特異値分解の浮動小数点コプロセッサによる高速化

長方行列向け特異値分解の
浮動小数点コプロセッサによる
高速化
2007年1月18日 HPCS2007
深谷 猛，山本 有作 （名古屋大学）
畝山 多加志 （京都大学）
堀 玄，梅野 健 （理化学研究所）
概要
•
•
•
•
•
•
背景と目的
CSX600のアーキテクチャと性能
長方行列の特異値分解アルゴリズム
CSX600向けの高速化手法
性能評価
まとめと今後の課題
1
１. 背景と目的
• 長方行列の特異値分解
A ：m×n 密行列
=
×
×
U ： m×n 列直交行列
V ： n×n 直交行列
S ： n×n 対角行列
• 応用
–
–
–
–
画像処理
電子状態計算 （Filter Diagonalization Method）
文書検索 （Latent Semantic Indexing）
各種統計計算 （主成分分析，最小二乗法）
数万×数千以上の行列の特異値分解
高速化が必要
2
高い浮動小数点能力を持つ専用プロセッサの登場
• Cell
– 1個の汎用CPUコア + 8個の専用コア
– 256GFLOPS （単精度）
• GRAPE-DR
– 512個の演算コア
– 512GFLOPS（単精度）
– 384GFLOPS（倍精度）
• ClearSpeed CSX600
– 1個の汎用コア + 96個の演算コア
– 48GFLOPS （倍精度）
3
専用プロセッサの特徴
• 多数の浮動小数点演算器
– 極めて高いGFLOPS値
• メモリアクセス性能は相対的に低下
– バンド幅はあまり大きくなっていない
– PCIバス経由でコプロセッサとして使う場合は，特に顕著
メモリネックを軽減できるアルゴリズムが性能向上の鍵
4
Level-3 BLASの利用
• 行列乗算 C:=C+AB
C
=
C
+
A
×
B
データ量： O(N 2)
演算量： O(N 3)
– 演算量に比べ，データ量は O(1/N)
– キャッシュをうまく使うことで，メモリネックを軽減可能
※行列ベクトル積（ y := y + Ax）では，データ量・演算量ともにO(N 2)
行列乗算を効率的に使えれば，一般のアルゴリズムも高速化可能
（LAPACKの基本的な考え方）
5
本研究の目的
• 長方行列の特異値分解アルゴリズムをCSX600を利用して
高速化する
• 既存のアルゴリズムをそのまま使用せず，行列乗算をなる
べく効率的に使う形で CSX600 向けに最適化する
• 性能評価を行い，更なる高速化に向けての課題を明らかに
する
6
２. CSX600のアーキテクチャと性能
• CSX600 チップ
– 1個の主プロセッサ
– 96個の演算専用プロセッサ
• 64ビット
• 倍精度2演算/サイクル
• 128B レジスタファイル
• 6KB SRAM
– 250MHz で動作
– ピーク性能： 48GFLOPS
• ClearSpeed Advance ボード
–
–
–
–
2個の CSX600 プロセッサ
1GB DRAM
PCI-X バスにより PC 本体と接続
ピーク性能： 96GFLOPS
7
CSX600の利用環境
• Software Development Kit
– コンパイラ： Cn 言語によるチップ内並列プログラミング
– デバッガ
– シミュレータ
• CSXL ライブラリ
–
–
–
–
–
今回はこれを利用
ClearSpeed Advance ボード用の BLAS
行列データを PC の主メモリ上に置いてコール
PC ⇔ ボード間の転送はライブラリ内で実行
公称性能： DGEMM（行列乗算）で sustained 50GFLOPS
現状では非常に制約が多い
• 使用可能なルーチンは DGEMM と DGETRF のみ
• DGEMM では，M，N，K ≧ 448 でないと動かない
• CSFFT ライブラリ
8
CSXL DGEMMの性能 (1)
• サイズパラメータのうち2個が450の場合
性能（MFLOPS）
50000
40000
N
A, B 非転置
A 非転置，B 転置
30000
C
K
+= A ×
B
M
20000
・ M = K = 450
・ 1000 ≦ N ≦ 6000
10000
0
N
1000 2000 3000 4000 5000 6000
– 性能は最大でも 11GFLOPS 程度
– 他の2個が450の場合でも，ほぼ同じ結果
CSX600の性能が活かせない
9
CSXL DGEMMの性能 (2)
• サイズパラメータのうち1個のみが450の場合
50000
性能（MFLOPS）
N
K
40000
30000
20000
10000
C
A, B 非転置
A 非転置，B 転置
0
1000 2000 3000 4000 5000 6000
+= A ×
B
M
N
・ K = 450
・ 1000 ≦ M = N ≦ 6000
– M = N が大きい場合，性能は 40GFLOPS 以上
– 他の2個が大きい場合でも，ほぼ類似の性能
性能を引き出すには，サイズパラメータのうち少なくとも2個を
大きい値にとる必要あり
10
３. 長方行列の特異値分解アルゴリズム
m×n 密行列 A
QR分解 A=QR
2mn2
二重対角化
(8/3)n3
特異値分解
O(n2) ~ O(n3)
逆変換
2n3 ~ 4n3
Qの乗算
4mn2
n×n 上三角行列 R
二重対角行列 B
Bの特異値・特異ベクトル
（B と R の特異値は同じ）
R の右特異
ベクトル
=
R の左特異
ベクトル
A の左特異
ベクトル
A の右特異
ベクトル
11
M >> N の場合の実行時間
• M = 40000 の場合
250
200
実
行
時
間
（
秒
）
150
環境
・ Xeon 3.2GHz
・ g77 -O3 + Intel MKL
Qの乗算
逆変換
特異値分解
二重対角化
QR分解
各部分のアルゴリズム（詳細は後述）
・ QR分解： ブロック化 + 再帰的QR
・ 二重対角化： Bischof + 村田法
・ 特異値分解： 分割統治法
100
50
0
N=1000
N=2000
QR分解とQの乗算の時間が支配的
12
４. CSX600向けの高速化手法
高速化の方針
• QR分解とQの乗算
– 行列乗算を効率的に利用できるようアルゴリズムを変更
– CSXLで実行
• 二重対角化および逆変換
– 二重対角化はBischof のアルゴリズムと村田法を利用
– CSXLは利用しない（CPUのみで実行）
• 二重対角行列の特異値分解
– 行列乗算中心のアルゴリズムである分割統治法を利用
– LAPACKのDBDSDCを利用（CSXLが使える部分は使う）
13
ハウスホルダー変換によるQR分解
ハウスホルダー変換による列の消去
H1 A = ( I – t1 y1 y1T ) A
= A(1)
level-2 BLAS
CSXLを使えない
上三角化とQR分解
Hn・・・ H2 H1 A = A(n)
A = H1 H2・・・ Hn A(n) = QR
・・・
A
A(1)
A(2)
A(n) = R
14
複数のハウスホルダー変換の合成
Compact WY representation
Hn・・・ H2 H1 = ( I – tn yn ynT ) ・・・ ( I – t2 y2 y2T )( I – t1 y1 y1T )
= I – Yn Tn YnT
ただし， Yn = [ y1 | y2 | ・・・ | yn] （m×n 行列）
Tn： n×n 下三角行列
I–
× ×
・・・ I –
× ×
= I–
×
×
複数のハウスホルダー変換の作用を，行列乗算として
まとめて実行可能
15
QR分解における行列乗算の利用法 Ⅰ
単純なブロック化（LAPACKで使用）
分解
L
更新
分解
I – Y 2 T2 Y 2 T
I – Y 1 T1 Y 1 T
更新
分解
I – Y 3 T3 Y 3 T
H3H2H1= I – Y1T1Y1T
H1
H2
H3
•ブロックのQR分解はlevel-2 BLAS
•演算量： L (ブロック幅) << n のとき 2mn2
•CSXLを使うにはL≧450 ⇒ level-2 BLASの演算量の増加
16
QR分解における行列乗算の利用法 Ⅱ
再帰的QR分解
分解
更新
分解
合成
(I – YL TL YLT )
×(I – YR TR YRT )
= I – Y T YT
I – YL TL YLT
I – YR TR YRT
I – YL TL YLT
•ほぼ全ての演算が行列乗算
I – Y’L T’L Y’LT
•演算量： 3mn2
17
QR分解における行列乗算の利用法 Ⅲ
本研究で採用した方法
L
更新
分解
分解
I – Y 2 T2 Y 2 T
I – Y 1 T1 Y 1 T
更新
分解
I – Y 3 T3 Y 3 T
I – Y 1 T1 Y 1 T
•ほぼ全ての演算が行列乗算
•演算量：Ⅰ(2mn2)とⅡ(3mn2)の中間
•ブロック幅Lを適切に選ぶことで，他の２方法より高速になる可能性がある
18
Q の乗算
Ⅰ（単純なブロック化）・Ⅲ（本研究で採用した方法）の場合
Q = ( I – Yn/L Tn/L Yn/LT ) ・・・ (I – Y1 T1 Y1T )
I–
L
×
・・・ I –
×
×
×
Ⅱ（再帰的QR分解）の場合
Q = I – Y T YT
n
I–
×
×
•行列サイズの大きいⅡの方がCSXLの性能を引き出せる
•Ⅰ・Ⅲの方が使用するメモリの量が少ない
19
５. 性能評価
評価環境
• Xeon 3.2GHz，メモリ2GB
• g77 -O3 + Intel Math Kernel Library
• ClearSpeed Advance ボード
評価例題
• [-0.5,0.5] の一様乱数を要素とする m×n 乱数行列の特異値分解
• 10000 ≦ m ≦ 40000，1000 ≦ n ≦ 3000
評価項目
• ClearSpeed ボードでの3種のQR分解法の性能比較
• 特異値分解全体の ClearSpeed ボードによる高速化効果
• 精度評価
20
ClearSpeed ボードでの3種のQR分解法の性能比較
140
m = 10000，n = 3000
130
実
行
時
間
（
秒
）
280
260
Qの乗算
120
240
QR分解
110
220
100
200
90
180
80
160
70
140
60
120
50
100
40
80
30
60
20
40
10
20
0
0
I
L=500
II
m = 20000，n = 3000
III
III
L=500
L=1000
I
L=500
•方法ⅠはQR分解が極めて遅い
•方法ⅡはQの乗算が最高速
•方法ⅢはLを適切に選ぶことで，全体として最高速になる
II
III
III
L=500
L=1000
21
特異値分解全体の ClearSpeed ボードによる高速化効果
m = 10000
90
80
実
行
時
間
（
秒
）
m = 40000
250
Qの乗算
逆変換
70
60
50
200
特異値分解
二重対角化
150
QR分解(方法Ⅲ)
40
100
30
20
50
10
0
0
MKL
CS
n = 1000
MKL
CS
n = 2000
•QR分解は最大で2倍程度高速化
•Qの乗算は最大で5倍程度高速化
•特異値分解全体としては最大で2.3倍高速化
MKL
n = 1000
CS
MKL
CS
n = 2000
22
m，n による加速率の変化
特異値分解全体の加速率
nによる性能変化（m=30000）
40000
2.5
2
加
速
率
1.5
1
0.5
40000
30000
20000
0
3000
2000
n
10000
m
1000
性能（MFLOPS）
35000
QR分解
Qの乗算
30000
25000
20000
15000
10000
5000
0
1000
2000
3000
n
加速率はm, n が大きくなるほど増大する
•m / n の増加： 特異値分解全体に対してCSXLを使える部分の割合が増加
•n の増加 ： QR分解・Qの乗算の部分におけるCSXLの性能が向上
ただし， CSXLを使わない部分（Rの特異値分解）の演算量も増加
23
精度評価
左特異ベクトルの直交性
残差
1.00E-10
1.00E-10
CS
∥US VT – A∥F
∥UTU – I∥F
MKL
1.00E-11
1.00E-12
1.00E-13
m ： 10000
n：
20000 30000 40000 10000 20000 30000 40000 10000 20000 30000
1000
2000
3000
1.00E-11
1.00E-12
1.00E-13
m ： 10000 20000
n：
30000 40000 10000 20000 30000 40000 10000 20000 30000
1000
2000
3000
•左特異ベクトルの直交性，残差ともCSの方が精度が悪い
•悪化の程度は一桁以下
24
６. まとめと今後の課題
本研究のまとめ
• ClearSpeed Advance ボードと CSXL を用いて，長方行列向け特異値
分解の高速化を行った
• QR分解のアルゴリズムを行列乗算を効果的に使う形に改良することで，
ボードによる高速化の効果を高めた
• 特異値分解全体では， 40000×2000 程度の行列においてボードを使う
ことで2.3倍の高速化が得られた
• 行列の規模が大きくなるほど高速化が期待できる
今後の課題
• より大規模な行列で性能評価
• SDKの使用による高速化
• 他の行列計算アルゴリズムへのCSX600の適用
25