Transcript Document

Level-3 BLASに基づく二重対角化
アルゴリズムとその性能評価
名古屋大学 計算理工学専攻
山本有作
2006年1月5日-7日
第3回計算数学研究会
目次
１. はじめに
２. 従来の二重対角化アルゴリズム
３. Level-3 BLAS に基づく二重対角化アルゴリズム
４. 性能評価
５. 特異ベクトルの計算手法
６. まとめと今後の課題
１. はじめに
• 本研究で対象とする問題
– 実正方行列 A の下二重対角行列への変換 B = U0TAV0
• A ： n×n 密行列
• B ： n×n 下二重対角行列
• U0，V0 ： n×n 直交行列
• 応用分野
– 行列 A の特異値分解に対する前処理
• 画像処理の分野では，n = 10,000 以上の大規模行列に対す
る特異値分解の需要あり
• 必要な特異ベクトルは数本程度の場合が多い
• 前処理を数分程度で行えることが望ましい
特異値分解の流れ
密行列 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アルゴリズム
逆変換
各部分の演算量
• m 組の特異値・特異ベクトルを求める場合の演算量
–
–
–
–
二重対角化：
特異値の計算：
特異ベクトルの計算：
逆変換：
(8/3) n3
O(mn) ～ O(n2)
O(mn) ～ O(m2n)
O(mn2)
– m ≪ n の場合は演算量の大部分を二重対角化が占める。
– ハウスホルダー法の高速化が必要
２. 従来の二重対角化アルゴリズム
• ハウスホルダー法による二重対角化
ベクトル
– 鏡像変換 H = I – a wwT による列の消去
• Hは対称な直交行列
• 与えられたベクトルの第１成分以外をゼロにする。
左からH
を乗算
• 第 k ステップでの処理
0
0
0
0
0
右からHkR
を乗算
0
左からHkL
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
ハウスホルダー法のアルゴリズム
[Step 1] k = 1 から n –1 まで以下の[Step 2] ～ [Step 8]を繰り返す。
[Step 2] A(k, k+1:n) を 0 にする鏡像変換 HkR = I – akR wkR (wkR)T の計
算
[Step 3] 行列ベクトル積： p := akR A(k:n, k:n) wkR
[Step 4] 行列のrank-1更新： A(k:n, k:n) := A(k:n, k:n) – p(wkR)T
[Step 5] A(k+2:n, k) を 0 にする鏡像変換 HkL = I – akL wkL (wkL)T の計算
[Step 6] 行列ベクトル積： qT := akL (wkL)T A(k+1:n, k:n)
[Step 7] 行列のrank-1更新： A(k +1:n, k:n) := A(k +1:n, k:n) – wkLqT
0
0
A(k, k+1:n)
A(k:n, k:n)
A(k +2:n, k)
A(k +1:n, k :n)
k
0
0
ハウスホルダー法の特徴と問題点
• 特徴
– 全演算量のほとんどが，level-2 BLAS である行列ベクトル積と
行列のrank-1更新で占められる。
• 全演算量：
約 (8/3)n3
• 行列ベクトル積の演算量：
約 (4/3)n3
• rank-1更新の演算量：
約 (4/3)n3
• 問題点
– level-2 BLASのため，行列データの再利用性なし
• キャッシュミス，SMPでのメモリ競合の影響大
– 最近のマイクロプロセッサ，SMP上での高性能が期待できない。
３. Level-3 BLAS に基づく二重対角化アルゴリズム
• 基本的なアイディア
– 密行列 A をまず帯幅 L の下三角帯行列 C に変換
– 次にこの帯行列を下二重対角行列 B に変換
次数 n
下三角
帯行列化
約 (8/3)n3
A
0
村田法
の拡張
O(n2L)
0
C 帯幅 L
0
0
B
• 二重対角化を２段階で行うことの利点
– 下三角帯行列への変換は， level-3 BLAS のみを使って実行可能
– 下三角帯行列から二重対角行列への変換の演算量は O(n2L) で
あり，前半部に比べてずっと小さい。
（参考）対称行列の三重対角化の場合
• Bischof のアルゴリズム（Bishof et al., 1993, 1994）
– 対称密行列 A をまず半帯幅 L の対称帯行列 C に変換
– 次にこの帯行列を三重対角行列 T に変換
次数 n
半帯幅 L
帯行列化
約 (4/3)n3
A
0
村田法
O(n2L)
0
C
0
0
T
• Bischof のアルゴリズムの性能・精度（山本, 2005）
– 固有値の誤差は，LAPACKで使われる Dongarra のアルゴリズム
に比べ，多くの場合2～3倍程度
– 速度は2～3倍高速
– 様々なマイクロプロセッサ上で，ピークの50～70%の性能を達成
下三角帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 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におけるメモリ競合の影響を低減可能
下三角帯行列の二重対角化
• 対称帯行列の三重対角行列への相似変換（村田法）
– 第 k 列の三重対角化
• 第 k 列の副対角要素より下を 0 にする鏡像変換を左と右からかける。
– Bulge-chasing
• 上記の結果，帯より下に非零要素がはみ出すので，相似変換を繰り
返して，はみ出した分を右下に移動し，右下隅から追い出す。
• 村田法のアイディアの二重対角化への適用（Lang, 1996）
– 第 k 列の二重対角化
• 第 k 列の副対角要素より下を 0 にする鏡像変換を左からかける。
– Bulge-chasing
• 上記の結果，上三角部分に非零要素がはみ出すので，左右から直
交変換を繰り返し掛けることにより，はみ出した分を右下に移動し，
右下隅から追い出す。
第1列の二重対角化と bulge-chasing
左側からの
直交変換で
更新された
要素
左側からの
直交変換で
更新された
要素
第2列の二重対角化と bulge-chasing
左側からの
直交変換で
更新された
要素
左側からの
直交変換で
更新された
要素
演算量は 8n2L
SMP向けの並列化
• 下三角帯行列化
– 並列版の Level-3 BLAS を使えば，逐次版のプログラムをその
まま SMP 上で並列化可能
– OpenMP などを用いて手動で並列化すれば，並列化効率を更に
向上可能
SMP向けの並列化（続き）
• 二重対角化向け村田法
– 第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 ステップ ・・・
が同時に実行可能
４. 性能評価
• 評価環境
– Opteron (1.8GHz), 1～4CPU
• Linux 上で PGI Fortran + GOTO BLAS を使用
• ピーク性能： 3.6GFLOPS/CPU
– PowerPC G5 (2.0GHz), 1～2CPU
• Mac OS X 上で IBM XL Fortran + GOTO BLAS を使用
• ピーク性能： 8GFLOPS/CPU
• 評価対象
– n = 2500 ～ 12500 の乱数行列の二重対角化
– 下三角帯行列化部分は並列版 GOTO BLAS，村田法部分は OpenMP
により並列化
• 評価方法
– 下三角帯行列化 + 村田法の合計時間を測定し，演算量を(8/3)n3 として
MFLOPS 値を算出
– ブロックサイズ L は，各 n に対して最適な値を使用
Opteron (1.8GHz) 1CPU上での性能
Opteron 1CPU
3500
L=100
Performance (GFLOPS)
3000
ピークの
80%
2500
2000
Opteron 1CPU
1500
1000
500
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき，本アルゴリズムはピークの80%の性能を達成
Opteron (1.8GHz) 4CPU上での性能
Performance (MFLOPS)
12000
Opteron 1CPU
Opteron 4CPUs
10000
L=100
ピークの
67%
286s
8000
6000
3.3倍
4000
L=100
2000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき，本アルゴリズムは4CPUで3.3倍の加速率
ピークの67%の性能を達成
各部分の演算時間（Opteron, n=12500）
2500
Execution time (s)
2000
村田法
下三角帯行列化
1500
1000
500
0
L=25
1CPU
L=50
1CPU
L=100
1CPU
L=25
4CPU
L=50
4CPU
L=100
4CPU
大規模問題に対しては，村田法の占める時間の割合は小さい。
村田法の加速率は，L=100 のとき4CPUで3.3倍
PowerPC G5 (2.0GHz) 1CPU上での性能
PowerPC G5 1CPU
8000
Performance (GFLOPS)
7000
PowerPC G5 1CPU
6000
L=100
5000
ピークの
57%
4000
3000
2000
1000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき，本アルゴリズムはピークの57%の性能を達成
PowerPC G5 (2.0GHz) 2CPU上での性能
8000
Performance (MFLOPS)
L=100
PowerPC G5 1CPU
Power PC G5 2CPUs
7000
6000
421s
5000
ピークの
42%
L=100
4000
1.5倍
3000
2000
1000
0
2500
5000
7500
10000
12500
Matrix size
n = 12500 のとき，本アルゴリズムは2CPUで1.5倍の加速率
ピークの42%の性能を達成
各部分の演算時間（PowerPC G5, n=12500）
1800
1600
村田法
下三角帯行列化
Execution time (s)
1400
1200
1000
800
600
400
200
0
L=25
1CPU
L=50
1CPU
L=100
1CPU
L=25
2CPU
L=50
2CPU
L=100
2CPU
大規模問題に対しては，村田法の占める時間の割合は小さい。
村田法の加速率は，L=100 のとき2CPUで1.4倍
５. 特異ベクトルの計算手法
• 方法１： 二重対角行列の特異ベクトルを計算して２回逆変換
n
L
0
0
0
0
A
A の特異ﾍﾞｸﾄﾙ
{ui }{vi }
C
2mn2
C の特異ﾍﾞｸﾄﾙ
{zi }{wi }
B
2mn2
B の特異ﾍﾞｸﾄﾙ
{xi }{yi }
特異値
{σi }
QR法
DC法
MR3
I-SVD
• 長所
– 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能
• 短所
– 逆変換の演算量が 4mn2 （ただし，m が小さければ影響は少ない）
– 村田法の変換をすべて記憶するため，n2 の記憶領域が余計に必要
５. 特異ベクトルの計算手法
• 方法２： 下三角帯行列の特異ベクトルを直接計算
n
L
0
0
0
A
A の特異ﾍﾞｸﾄﾙ
{ui }{vi }
C
2mn2
0
C の特異ﾍﾞｸﾄﾙ
{zi }{wi }
帯行列用
逆反復法
O(mnL2+ m2n)
B
特異値
{σi }
QR法
dqds法
mdLVs法
二分法
• 長所
– 村田法の逆変換のための演算量（2mn2）・記憶領域（n2）が不要
• 短所
– 帯行列用逆反復法（O(mnL2)）と特異ベクトルの直交化（O(m2n)）が必要
– 特異値σi は C の高精度な固有値でないため，ツイスト分解は使用不可
５. 特異ベクトルの計算手法
• 方法３： 村田法を使わず帯行列 C の特異値・特異ベクトルを直接計算
n
L
0
帯行列版
mdLVs法
特異値
{σi }
0
A
A の特異ﾍﾞｸﾄﾙ
{ui }{vi }
C
2mn2
帯行列版
ﾂｲｽﾄ分解
C の特異ﾍﾞｸﾄﾙ
{zi }{wi }
• 長所
– 村田法の逆変換のための演算量（2mn2）・記憶領域（n2）が不要
– うまく行けば，直交化が不要なアルゴリズムを作れる可能性あり
• 問題点
– そもそも帯行列版の mdLVs 法，ツイスト分解は構成できるか？
– 演算量，演算の安定性は？