Transcript I – n

応用数理工学特論 第10回
2008年7月3日
計算理工学専攻
山本有作
目次
１. はじめに
２. 共有メモリ型並列計算機での特異値分解の高速化
～ 正方行列の特異値分解 ～
３. SIMD超並列プロセッサでの特異値分解の高速化
～ 長方行列の特異値分解 ～
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
バス
• 例１： マルチコアプロセッサ
メモリ
– Intel社 デュアルコアXeon
– 1チップに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
２. 共有メモリ型並列計算機での
特異値分解の高速化
～ 正方行列の特異値分解 ～
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は対称な直交行列
• 与えられたベクトルの第１成分以外をゼロにする
• 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
• 二重対角化を２段階で行うことの利点
– 下三角帯行列への変換は， 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 の利用により，二重対角化の性能を向上可能
• 同様のアイディアに基づく三重対角化アルゴリズムでは，高い単体性能・並
列性能を確認済み
• 短所
– 特異ベクトル計算のための計算量・記憶領域が増大
• ２段階の逆変換が必要
• 詳しくは次のスライドで説明
– 二重対角化の高速化効果が大きければ，計算量増大を考慮し
ても全体としては高速化できると予想
– 特に，求める特異ベクトルが少ない場合は効果が大きいはず。
22
特異ベクトルの計算手法
• 二重対角行列の特異ベクトルを計算して２回逆変換
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
３. 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
性能を引き出すにはサイズパラメータのうち少なくとも
２個を大きい値にとる必要がある
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 ）
=
CSX６００で高速化
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を適切に選ぶことで，他の２方法より高速になる可能性がある
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使用時の各手法の性能比較（１）
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使用時の各手法の性能比較（２）
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
行列サイズ による高速化率の変化（１）
§再帰的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
行列サイズ による高速化率の変化（２）
§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
まとめと今後の課題
§まとめ
• CSX６００による長方行列特異値分解の高速化を行った
• CSX６００に適したQR分解アルゴリズムの選択を行った
• 適切なQR分解のアルゴリズムを選ぶことで，最大で３倍
程度の高速化ができた
§今後の課題
• よりCSX600に適したQR分解アルゴリズムの開発
• CSX600による Rの特異値分解の高速化
• 他の行列計算へのCSX６００の適用
62