スライド 0

応用数理工学特論 第10回
2008年7月3日
計算理工学専攻
山本有作
目次
1. はじめに
2. 共有メモリ型並列計算機での特異値分解の高速化
~ 正方行列の特異値分解 ~
3. SIMD超並列プロセッサでの特異値分解の高速化
~ 長方行列の特異値分解 ~
1
1. はじめに
2
特異値分解とは
• 任意の実正方行列 A は,直交行列 U,対角行列 S,直交
行列 VT の積に分解できる。
A :n×n 実行列
n
=
n
×
×
n
U :n×n 直交行列
V : n×n 直交行列
S : n×n 対角行列
n
• 参考
• A が実対称行列の場合,A=UDUT
• A が任意の実行列の場合,A=SDS–1
• A が任意の実行列の場合,A=URUT
3
特異値と特異ベクトル
• A=USV T において,
– V の各列 vi を右特異ベクトル
– U の各列 ui を左特異ベクトル
– S の要素 si を特異値,と呼ぶ。
• 特異ベクトルの性質
– AV=US より, Avi = si ui
– UTA=S V T より, uiTA = si viT
4
固有値分解(対角化)との関係
• A=USV T より, ATA= VSUTUSV T = VS 2V T
よって,V は ATA を対角化する直交行列
A の右特異ベクトル = ATA の固有ベクトル
• A の特異値 = ATA の固有値の平方根
• 一方, AAT= USV TVSU T = US 2U T
よって,V は ATA を対角化する直交行列
A の左特異ベクトル = AAT の固有ベクトル
5
特異値分解を求めるアルゴリズム
• 二重対角行列への変換
U0TAV0 = B (U0, V0: 直交行列)
• 二重対角行列の特異値分解
U1TBV1 = S
(U1, V1: 直交行列)
• A の特異値分解の計算
U1TU0TAV0V1 = S
よって,U= U0U1 ,V= V0V1 とおくと(逆変換),
A = USV T
6
特異値分解の応用
• 信号処理
• 画像処理
• 電子状態計算
– Filter Diagonalization Method
• 文書検索
– Latent Semantic Indexing
• 各種統計計算
– 主成分分析
– 独立成分分析
– 最小二乗法
7
対象とする計算機
8
共有メモリ型並列計算機
• アーキテクチャ
– 複数のプロセッサ(PU)がバスを
通してメモリを共有
– PUはそれぞれキャッシュを持つ
キャッシュ PU0
PU1
PU2
PU3
バス
• 例1: マルチコアプロセッサ
メモリ
– Intel社 デュアルコアXeon
– 1チップに2個のプロセッサを共有メモリ型で集積
• 例2: 共有メモリ型のスーパーコンピュータ
– 富士通 PrimePower HPC2500
– 最大128個のPUを共有メモリで結合
デュアルコア Xeon
PrimePower HPC2500
9
SIMD超並列プロセッサ
• 1チップに100個単位のプロセッサを搭載
• 全プロセッサが別々のデータに対して同じ命令を実行する
SIMD(Single Instruction Multiple Data)型の並列計算機
– 線形計算に最適
• 汎用プロセッサの10~100倍の性能
– ClearSpeed社CSX600: 96プロセッサ,48GFLOPS(倍精度)
– 東大GRAPE-DR: 512プロセッサ,256GFLOPS(倍精度)
• 演算性能は極めて高いが,データ転送速度は汎用プロセッ
サと同程度
• データ再利用性の向上が,共有メモリ型並列計算機以上に
重要
10
512コアプロセッサ(GRAPE-DR)
11
グラフィックスプロセッサ(GPU)
• GPU(Graphics Processing Unit)の高速化
– CPUを大きく上回るペースで演算性能が向上
– グラフィックスメモリも大容量化・高速化
nVIDIA社 GeForce8800GTX
・ 240個の演算プロセッサ
・ 演算性能: 933GFLOPS(単精度)
90GFLOPS(倍精度)
・ メモリ: 1GB
GPUを汎用の数値計算に使うGPGPU(General-Purpose
GPU)が注目を集める
– 計算機としてのアーキテクチャはSIMD超並列型
12
2. 共有メモリ型並列計算機での
特異値分解の高速化
~ 正方行列の特異値分解 ~
13
正方行列の特異値分解
A :n×n 密行列
n
=
n
×
×
n
U :n×n 直交行列
V : n×n 直交行列
S : n×n 対角行列
n
<応用>
– 各種統計計算
– より応用の多い長方行列の特異値分解は,正方行列の特異値分解に
帰着される。
14
従来の特異値分解アルゴリズムとその問題点
密行列 A
計算内容
計算手法
二重対角化
U0TAV0 = B
ハウスホルダー法
(U0, V0: 直交行列)
二重対角行列 B
二重対角行列の
特異値・特異ベクトル計算
Bvi =σi xi
BTxi =σi yi
Bの特異値 {σi },
特異ベクトル {xi }{yi }
vi = V0 yi
逆変換
ui = U 0 x i
QR法
分割統治法
MR3アルゴリズム
I-SVDアルゴリズム
逆変換
Aの特異ベクトル {ui }, {vi }
15
各部分の演算量と実行時間
密行列 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
・二重対角化が実行時間の
大部分を占める
・速度向上率が低い
16
特異値分解アルゴリズムの高速化
• ハウスホルダー法による二重対角化
ベクトル
– 鏡像変換 H = I – a wwT による列の消去
• Hは対称な直交行列
• 与えられたベクトルの第1成分以外をゼロにする
• HA(k) = A(k) – au (utA(k))
行列ベクトル積
rank-1更新
左からH
を乗算
• 第 k ステップでの処理
0
0
A(k)
0
右からHkR
0
を乗算
0
左からHkL
0
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
17
ハウスホルダー法の特徴と問題点
• 特徴
– 全演算量のほとんどが,BLAS2 である行列ベクトル積と行列の
rank-1更新で占められる。
• 全演算量:
約 (8/3)n3
• 行列ベクトル積の演算量: 約 (4/3)n3
• rank-1更新の演算量:
約 (4/3)n3
• 問題点
– BLAS2のため,行列データの再利用性なし
• キャッシュミス,メモリ競合の影響大
– 共有メモリ型並列計算機上で高性能が出ない理由
18
BLAS3 に基づく二重対角化アルゴリズム
• 基本的なアイディア
– 密行列 A をまず帯幅 L の下三角帯行列 C に変換
– 次にこの帯行列を下二重対角行列 B に変換
次数 n
下三角
帯行列化
約 (8/3)n3
A
0
村田法
O(n2L)
0
C 帯幅 L
0
0
B
• 二重対角化を2段階で行うことの利点
– 下三角帯行列への変換は, BLAS3 のみを使って実行可能
キャッシュの有効利用が可能
– 下三角帯行列から二重対角行列への変換の演算量は O(n2L)
• 前半部に比べてずっと小さい。
19
下三角帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – WαWT
• Hは直交行列
• 与えられたブロックベクトルを上三角
行列(正確には右上三角部分のみ
非零でそれ以外が零の行列)に変形
• HA(k) = A(k) – WαWTA(k)
第 K ステップでの処理
行列積 = BLAS3
0
0
A(k)
左からH
を乗算
0
0
右からHK
を乗算
非ゼロ要素
R
0
ゼロにしたい部分
左からHKL
を乗算
0
影響を受ける部分
20
帯行列の二重対角化(村田法)
左側からの
直交変換で
更新された
要素
左側からの
直交変換で
更新された
要素
・演算量は 8n2L
・並列化も可能
21
本アルゴリズムの長所と短所
• 長所
– BLAS3 の利用により,二重対角化の性能を向上可能
• 同様のアイディアに基づく三重対角化アルゴリズムでは,高い単体性能・並
列性能を確認済み
• 短所
– 特異ベクトル計算のための計算量・記憶領域が増大
• 2段階の逆変換が必要
• 詳しくは次のスライドで説明
– 二重対角化の高速化効果が大きければ,計算量増大を考慮し
ても全体としては高速化できると予想
– 特に,求める特異ベクトルが少ない場合は効果が大きいはず。
22
特異ベクトルの計算手法
• 二重対角行列の特異ベクトルを計算して2回逆変換
n
L
0
0
0
0
A
A の特異ベクトル
{ui }{vi }
C
4mn2
C の特異ベクトル
{zi }{wi }
B
4mn2
B の特異ベクトル
{xi }{yi }
特異値
{σi }
QR法
DC法
MR3
I-SVD
• 長所
– 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能
• 短所
– 逆変換の演算量が 8mn2 (従来法の2倍)。ただし BLAS3 で実行可能
– 村田法の変換をすべて記憶するため,n2 の記憶領域が余計に必要
23
性能評価
24
性能評価
• 評価環境
– 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
• 評価対象・条件
– BLAS3 に基づくアルゴリズムと LAPACK の性能を比較
– n = 1200 ~ 20000 の乱数行列の特異値分解(全特異ベクトルを計算)
• BLAS3 アルゴリズムにとっては一番不利な条件
– BLAS3 アルゴリズムの L(半帯幅)は各 n ごとに最適値を使用
25
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
• 結果
– BLAS3 アルゴリズムでは PU 数に応じて実行時間が順調に減少
• 4PUでの加速率は 3~3.2 倍
– 4PU の場合は BLAS3 アルゴリズムが従来法より高速
26
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
• 結果
– BLAS3 アルゴリズムは従来法に比べて最大3.5倍高速
– プロセッサ数が多いとき加速率が鈍るのは,非並列化部分(ブロッ
ク鏡像変換の作成など)の影響と思われる。
27
両手法の実行時間の内訳
• Xeon,n=5000の場合
800
600
逆変換
分割統治法
二重対角化
500
400
逆変換2
逆変換1
分割統治法
村田法
帯行列化
700
600
500
300
400
300
200
200
100
100
0
0
1
2
4
1
2
4
• 考察
– BLAS3 アルゴリズムでは,どの部分の実行時間も順調に減少
– 逆変換1(村田法の逆変換)の占める時間が大きい。
この部分について,さらに高速化が必要
– 必要な特異ベクトルの本数が少ない場合,BLAS3 アルゴリズム
はさらに有利
28
両手法の実行時間の内訳
• 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
• 考察
– BLAS3 アルゴリズムでは,どの部分の実行時間も順調に減少
– 従来法は,二重対角化の部分の加速が鈍い。
• ただし,32PUで6倍程度は加速
• メモリバンド幅が大きいためと思われる。
29
まとめと今後の課題
30
3. SIMD超並列プロセッサによる
特異値分解の高速化
~ 長方行列の特異値分解 ~
32
長方行列の特異値分解
A :m×n 密行列
m
=
n
×
×
n
U :m×n 列直交行列
V : n×n 直交行列
S : n×n 対角行列
n
<応用>
– 画像処理
– 電子状態計算 (Filter Diagonalization Method)
– 文書検索 (Latent Semantic Indexing)
– 各種統計計算 (主成分分析,最小二乗法)
(例)
10万
5千
33
高い演算能力を持つSIMD超並列プロセッサ
• ClearSpeed CSX600
– 1個の汎用コア + 96個の演算コア
– 48GFLOPS (倍精度)
• GRAPE-DR
– 512個の演算コア
– 512GFLOPS(単精度)
– 256GFLOPS(倍精度)
• GeForce GTX280 (GPU)
– 240個の演算コア
– 933GFLOPS(単精度)
– 90GFLOPS(倍精度)
•多数の演算装置の集積により極めて高いFLOPS値
•メモリアクセス性能の相対的な低下によるメモリネック
34
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)
行列積を効率的に使えれば,一般のアルゴリズムも高速化可能
35
本章の目的
• 長方行列の特異値分解アルゴリズムをCSX600を利用して
高速化する
• 既存のアルゴリズムをそのまま使用せず,行列積をなるべく
効率的に使う形で CSX600 向けに最適化する
• 性能評価を行い,更なる高速化に向けての課題を明らかに
する
36
ClearSpeed CSX600
について
37
CSX600のアーキテクチャと性能
•
CSX600 チップ
– 1個の主プロセッサ
– 96個の演算専用プロセッサ
• 64ビット
• 倍精度2演算/サイクル
• 128B レジスタファイル
• 6KB SRAM
– 250MHz で動作
– ピーク性能: 48GFLOPS
•
ClearSpeed Advance ボード
–
–
–
–
2個の CSX600 プロセッサ
1GB DRAM
PCI-X バスにより PC 本体と接続
ピーク性能: 96GFLOPS
38
CSX600の利用環境
• Software Development Kit
– コンパイラ: Cn 言語によるチップ内並列プログラミング
– デバッガ
– シミュレータ
• CSXL ライブラリ
–
–
–
–
今回はこれを利用
ClearSpeed Advance ボード用の BLAS
行列データを PC の主メモリ上に置いてコール
PC ⇔ ボード間の転送はライブラリ内で実行
公称性能: DGEMM(行列乗算)で 50GFLOPS
• CSFFT ライブラリ
39
性能(MFLOPS)
CSXLのDGEMMの性能
50000
m = k = 450 1000 ≦ n ≦ 6000
50000
A,B:非転置
A:非転置,B:転置
40000
40000
30000
30000
20000
20000
10000
10000
0
0
1000 2000 3000 4000 5000 6000
n
m
C
k = 450 1000 ≦ m = n ≦ 6000
n
1000 2000 3000 4000 5000 6000
n
k
+= A × B
m
C
n,m
k
+= A × B
性能を引き出すにはサイズパラメータのうち少なくとも
2個を大きい値にとる必要がある
40
CSX600向けの
高速化手法
41
長方行列の特異値分解の手順
n
QR分解
A = QR
m
二重対角化
R = U1 B V1T
特異値分解
B = U2 S V2
U’ = U1 U2
逆変換
V = V1 V2
Qの乗算 U =
QU’
A
n
T
n
m
Q
n
R
R = U’ S V T
n
A= US V
n
T
B
42
長方行列の特異値分解の手順
m ≫ n (例:m =100000, n =5000)の場合
QR分解
計算量
2mn2
A = QR
二重対角化
R = U1 B V1T
(8/3)n3
特異値分解
B = U2 S V2T
O(n2) ~ O(n3)
U’ = U1 U2
逆変換
V = V1 V2
R = U’ S V T
2n3 ~ 4n3
Qの乗算 U = QU’
A= US V T
4mn2
実
行
時
間
の
大
部
分
を
占
め
る
43
高速化の方針
§CSX600を利用する部分
QR分解
行列乗算のみ高速化可能
A = QR
Qの乗算 U = QU’
A= US V T
行列乗算中心のアルゴリズム
§CSX600を利用しない(CPUのみで実行する)部分
二重対角化
R = U1 B V1T
U’ = U1 U2
逆変換
V = V1 V2
特異値分解
R = U’ S V T
B = U2 S V2T
LAPACKのDGEBRD
LAPACKのDORMBR
I-SVD(IDBDSLV)
44
ハウスホルダー変換による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
45
複数のハウスホルダー変換の合成
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–
×
×
複数のハウスホルダー変換の作用を,行列乗算として
まとめて実行可能
46
ハウスホルダー変換のブロック化
HL ・・・ H1
=
複数のハウスホルダー変換の合成
(I – Y T YT )
=
(I – Y T YT )
=
行列乗算としてまとめて作用
(I – Y T YT )
=
CSX600で高速化
47
(Ⅰ)ブロックQR分解(LAPACKで使用)
分解
L
I – Y 1 T1 Y 1 T
更新
分解
更新
I – Y 2 T2 Y 2 T
分解
I – Y 3 T3 Y 3 T
・・・行列乗算
H3H2H1= I – Y1T1Y1T
H1
H2
•ブロックの分解(
H3
)は行列・ベクトル積
•演算量: L (ブロック幅) << n のとき 2mn2
•CSX600を使うにはL≧450 ⇒行列・ベクトル積の演算量の増加
48
(Ⅱ)再帰的QR分解
分解
更新
分解
合成
(I – Y1T1Y1T)
×(I – Y2T2Y2T)
= I – Y T YT
I – Y 1 T1 Y 1 T
I – Y 2 T2 Y 2 T
I – Y 1 T1 Y 1
・・・行列乗算
T
•ほぼ全ての演算が行列乗算
I – Y11T11Y11
T
•演算量: 3mn2
•Qの乗算の効率が良い
49
(Ⅲ)再帰的QR分解の拡張
L
更新
分解
分解
I – Y 2 T2 Y 2 T
I – Y 1 T1 Y 1 T
I – Y 1 T1 Y 1
更新
分解
I – Y 3 T3 Y 3 T
・・・行列乗算
T
•ほぼ全ての演算が行列乗算
•演算量:Ⅰ(2mn2)とⅡ(3mn2)の中間
•ブロック幅Lを適切に選ぶことで,他の2方法より高速になる可能性がある
50
Q の乗算
(Ⅰ)ブロックQR分解・(Ⅲ)再帰的QR分解の拡張の場合
Q = ( I – Yn/L Tn/L Yn/LT ) ・・・ (I – Y1 T1 Y1T )
I–
L
×
・・・ I –
×
×
×
(Ⅱ)再帰的QR分解の場合
Q = I – Y T YT
n
I–
×
×
•行列サイズの大きいⅡの方がCSXLの性能を引き出せる
•Ⅰ・Ⅲの方が使用するメモリの量が少ない
51
性能評価
52
評価内容
評価環境
• Xeon 3.2GHz,メモリ8GB
• Intel Fortran -O3 + Intel Math Kernel Library
• ClearSpeed Advance ボード
評価例題
• [-0.5,0.5] の一様乱数を要素とする m×n 乱数行列の特異値分解
• 10000 ≦ m ≦ 100000,1000 ≦ n ≦ 4000
評価項目
• ClearSpeed ボードでの3種のQR分解法の性能比較
• 特異値分解全体の ClearSpeed ボードによる高速化効果
• 精度評価
53
CSX600使用時の各手法の性能比較(1)
m =100000 n =4000
2000
実
行
時
間
(
秒
)
1800
Qの乗算
1600
QR分解
1400
1200
1000
800
600
400
200
0
Ⅰ.ブロック
QR分解
Ⅱ.再帰的
QR分解
500
1000
1500
2000
Ⅲ.再帰的QR分解の拡張(数字はブロック幅)
54
CSX600使用時の各手法の性能比較(2)
m =10000 n =4000
160
実
行 140
時
間 120
(
秒 100
)
Qの乗算
QR分解
80
60
40
20
0
Ⅰ.ブロック Ⅱ.再帰的
QR分解
QR分解
500
1000
1500
2000
Ⅲ.再帰的QR分解の拡張(数字はブロック幅)
55
CSX600による特異値分解全体の高速化
実 25
行
時
間 20
(
秒
) 15
m = 10000 n=1000 (m:n = 10:1)
1.2倍
全体
Qの乗算
Rの特異値分解
QR分解
1.8倍
3000
m = 100000 n=4000 (m:n = 25:1)
1.3倍
2500
2000
3.1倍
1500
10
1000
5
500
0
0
PC
PC+CSX
LAPACKのDGESDD
PC
PC+CSX
再帰的QR分解を利用
PC
PC+CSX
LAPACKのDGESDD
PC
PC+CSX
再帰的QR分解を利用
56
行列サイズ による高速化率の変化(1)
§再帰的QR分解を用いた場合
5
4
高 3
速
化 2
率
100000
50000
1
0
10000
4000
3000
2000
m
1000
n
高速化率=PC単体での実行時間/PC+CSX600での実行時間
57
各部分の高速化率の変化
m = 50000
m = 100000
10
10
9
9
QR分解
8
高
速
化
率
8
Qの乗算
7
7
6
6
5
5
4
4
3
3
2
2
1
1
0
0
1000
2000
3000
n
4000
1000
2000
3000
4000
n
58
行列サイズ による高速化率の変化(2)
§LAPACKのDGESDDを用いた場合
5
4
高 3
速
化 2
率
100000
50000
1
0
10000
4000
3000
2000
m
1000
n
高速化率=PC単体での実行時間/PC+CSX600での実行時間
59
精度評価
左特異ベクトルの直交性
残差
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
60
まとめと今後の課題
61
まとめと今後の課題
§まとめ
• CSX600による長方行列特異値分解の高速化を行った
• CSX600に適したQR分解アルゴリズムの選択を行った
• 適切なQR分解のアルゴリズムを選ぶことで,最大で3倍
程度の高速化ができた
§今後の課題
• よりCSX600に適したQR分解アルゴリズムの開発
• CSX600による Rの特異値分解の高速化
• 他の行列計算へのCSX600の適用
62