Level-3 BLAS に基づく二重対角化アルゴリズムとその性能評価

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 QR法分割統治法 MR3アルゴリズム I-SVDアルゴリズム Bvi =σi xi BTxi =σi yi 二重対角行列の特異値・特異ﾍﾞｸﾄﾙ計算 Bの特異値 {σi}，　特異ベクトル {xi }{yi } vi = V0 yi ui = U0 xi 逆変換逆変換 Aの特異ベクトル {ui }, {vi }

各部分の演算量 • m組の特異値・特異ベクトルを求める場合の演算量 • 二重対角化： (8/3) n3 • 特異値の計算： O(mn) ～ O(n2) • 特異ベクトルの計算： O(mn) ～ O(m2n) • 逆変換： O(mn2) • m ≪ nの場合は演算量の大部分を二重対角化が占める。 • ハウスホルダー法の高速化が必要

２.　従来の二重対角化アルゴリズム ベクトル • ハウスホルダー法による二重対角化 • 鏡像変換 H = I – a wwTによる列の消去 • Hは対称な直交行列 • 与えられたベクトルの第１成分以外をゼロにする。 • 第 k ステップでの処理左からH を乗算 0 0 0 0 0 0 右からHkR を乗算左からHkL を乗算 k ゼロにしたい部分影響を受ける部分非ゼロ要素

ハウスホルダー法のアルゴリズム [Step 1]k = 1 から n –1 まで以下の[Step 2] ～ [Step 8]を繰り返す。 [Step 2]A(k, k+1:n) を 0 にする鏡像変換 HkR= I – akRwkR(wkR)Tの計算 [Step 3]行列ベクトル積：　p := akRA(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 – akLwkL(wkL)Tの計算 [Step 6]行列ベクトル積：　qT := akL (wkL)TA(k+1:n, k:n) [Step 7]行列のrank-1更新：　 A(k +1:n, k:n) := A(k +1:n, k:n) – wkLqT A(k, k+1:n) 0 0 A(k:n, k:n) 0 0 A(k +2:n, k) A(k +1:n, k :n) k

ハウスホルダー法の特徴と問題点 • 特徴 • 全演算量のほとんどが，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 に変換 • 二重対角化を２段階で行うことの利点 • 下三角帯行列への変換は， level-3BLAS のみを使って実行可能 • 下三角帯行列から二重対角行列への変換の演算量は O(n2L) であり，前半部に比べてずっと小さい。次数 n 下三角帯行列化村田法の拡張 0 0 O(n2L) 約 (8/3)n3 0 0 A C B 帯幅 L

（参考）対称行列の三重対角化の場合 • Bischof のアルゴリズム（Bishof et al., 1993, 1994） • 対称密行列 A をまず半帯幅 L の対称帯行列 C に変換 • 次にこの帯行列を三重対角行列 T に変換 • Bischof のアルゴリズムの性能・精度（山本, 2005） • 固有値の誤差は，LAPACKで使われる Dongarra のアルゴリズムに比べ，多くの場合2～3倍程度 • 速度は2～3倍高速 • 様々なマイクロプロセッサ上で，ピークの50～70%の性能を達成次数 n 半帯幅 L 0 帯行列化村田法 0 O(n2L) 約 (4/3)n3 0 0 A C T

下三角帯行列化のアルゴリズム ブロックベクトルブロック鏡像変換によるブロック列の消去 • ブロック鏡像変換 H = I – WαWT • Hは直交行列 • 与えられたブロックベクトルを上三角行列（正確には右上三角部分のみ非零でそれ以外が零の行列）に変形第 Kステップでの処理左からH を乗算 0 0 0 0 0 0 左からHKL を乗算右からHKR を乗算ゼロにしたい部分影響を受ける部分非ゼロ要素

下三角帯行列化のアルゴリズム（続き） [Step 1]K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 [Step 2]A(K, K:N) を上三角行列に変形する鏡像変換 HKR= I – WKRaKR (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 – WKLaKL (WKL)Tの計算 [Step 6]　行列・ﾌﾞﾛｯｸﾍﾞｸﾄﾙ積：　QT := aKL (WKL)TA(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 の第 k ステップ • 第2列に対する bulge-chasing の第 k–2ステップ • 第3列に対する bulge-chasing の第 k–4ステップ　・・・　　が同時に実行可能第2列のbulge-chasing における，右側からの第1の直交変換で更新される要素第1列のbulge-chasing における，右側からの第3の直交変換で更新される要素第1列による二重対角化は，今後　より右の要素にのみ影響を及ぼす。第1列の計算が右下まで行くのを待たずに，第2列の計算を開始できる。

４.　性能評価 • 評価環境 • 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上での性能 L=100 ピークの 80% Performance (GFLOPS) Matrix size n = 12500 のとき，本アルゴリズムはピークの80%の性能を達成

Opteron (1.8GHz) 4CPU上での性能 L=100 ピークの 67% 286s Performance (MFLOPS) 3.3倍 L=100 Matrix size n = 12500 のとき，本アルゴリズムは4CPUで3.3倍の加速率ピークの67%の性能を達成

各部分の演算時間（Opteron, n=12500） Execution time (s) 大規模問題に対しては，村田法の占める時間の割合は小さい。村田法の加速率は，L=100 のとき4CPUで3.3倍

PowerPC G5 (2.0GHz) 1CPU上での性能 L=100 ピークの 57% Performance (GFLOPS) Matrix size n = 12500 のとき，本アルゴリズムはピークの57%の性能を達成

PowerPC G5 (2.0GHz) 2CPU上での性能 L=100 ピークの 42% 421s L=100 Performance (MFLOPS) 1.5倍 Matrix size n = 12500 のとき，本アルゴリズムは2CPUで1.5倍の加速率ピークの42%の性能を達成

各部分の演算時間（PowerPC G5, n=12500） Execution time (s) 大規模問題に対しては，村田法の占める時間の割合は小さい。村田法の加速率は，L=100 のとき2CPUで1.4倍

５.　特異ベクトルの計算手法 n • 方法１：二重対角行列の特異ベクトルを計算して２回逆変換 • 長所 • 二重対角行列の特異値・特異ベクトルを求める任意の手法が適用可能 • 短所 • 逆変換の演算量が 4mn2（ただし，mが小さければ影響は少ない） • 村田法の変換をすべて記憶するため，n2の記憶領域が余計に必要 0 L 0 特異値 {σi} 0 0 QR法 DC法 MR3 I-SVD A C B 2mn2 2mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } B の特異ﾍﾞｸﾄﾙ {xi }{yi }

５.　特異ベクトルの計算手法 n • 方法２：下三角帯行列の特異ベクトルを直接計算 • 長所 • 村田法の逆変換のための演算量（2mn2）・記憶領域（n2）が不要 • 短所 • 帯行列用逆反復法（O(mnL2)）と特異ベクトルの直交化（O(m2n)）が必要 • 特異値σi は Cの高精度な固有値でないため，ツイスト分解は使用不可 0 L 0 0 0 QR法 dqds法 mdLVs法二分法帯行列用逆反復法 A C B 2mn2 O(mnL2+ m2n) A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi } 特異値 {σi}

５.　特異ベクトルの計算手法 n 帯行列版 mdLVs法 • 方法３：村田法を使わず帯行列 Cの特異値・特異ベクトルを直接計算 • 長所 • 村田法の逆変換のための演算量（2mn2）・記憶領域（n2）が不要 • うまく行けば，直交化が不要なアルゴリズムを作れる可能性あり • 問題点 • そもそも帯行列版の mdLVs 法，ツイスト分解は構成できるか？ • 演算量，演算の安定性は？ 0 L 特異値 {σi} 0 A C 帯行列版ﾂｲｽﾄ分解 2mn2 A の特異ﾍﾞｸﾄﾙ {ui }{vi } C の特異ﾍﾞｸﾄﾙ {zi }{wi }

Level-3 BLAS に基づく二重対角化アルゴリズムとその性能評価