Level-3 BLASに基づく二重対角化 アルゴリズムとその性能評価 名古屋大学　計算理工学専攻 山本有作 2006年1月5日-7日 第3回計算数学研究会 日立製作所の山本有作です。 「～」について発表いたします。

目次 １. はじめに ２. 従来の二重対角化アルゴリズム ３. Level-3 BLAS に基づく二重対角化アルゴリズム ４. 性能評価 １.　はじめに ２.　従来の二重対角化アルゴリズム ３.　Level-3 BLAS に基づく二重対角化アルゴリズム ４.　性能評価 ５.　特異ベクトルの計算手法 ６.　まとめと今後の課題 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１. はじめに 本研究で対象とする問題 応用分野 実正方行列 A の下二重対角行列への変換 B = U0TAV0 A ： n×n 密行列 １.　はじめに 本研究で対象とする問題 実正方行列 A の下二重対角行列への変換 B = U0TAV0 A ： n×n 密行列 B ： n×n 下二重対角行列 U0，V0 ： n×n 直交行列 応用分野 行列 A の特異値分解に対する前処理 画像処理の分野では，n = 10,000 以上の大規模行列に対する特異値分解の需要あり 必要な特異ベクトルは数本程度の場合が多い 前処理を数分程度で行えることが望ましい ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

特異値分解の流れ 密行列 A 計算内容 計算手法 U0TAV0 = B (U0, V0: 直交行列) ハウスホルダー法 二重対角化 QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム 二重対角行列の 特異値・特異ﾍﾞｸﾄﾙ計算 Bvi =σi xi BTxi =σi yi Bの特異値 {σi }， 　 特異ベクトル {xi }{yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 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 の場合は演算量の大部分を二重対角化が占める。 ハウスホルダー法の高速化が必要

２. 従来の二重対角化アルゴリズム ハウスホルダー法による二重対角化 第 k ステップでの処理 ２.　従来の二重対角化アルゴリズム ハウスホルダー法による二重対角化 鏡像変換 H = I – a wwT による列の消去 Hは対称な直交行列 与えられたベクトルの第１成分以外をゼロにする。 第 k ステップでの処理 ベクトル 左からH を乗算 右からHkR を乗算 左から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 A(k, k+1:n) A(k:n, k:n) 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-3 BLAS のみを使って実行可能 下三角帯行列から二重対角行列への変換の演算量は O(n2L) であり，前半部に比べてずっと小さい。 次数 n 下三角 帯行列化 村田法 の拡張 約 (8/3)n3 O(n2L) 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 帯行列化 村田法 約 (4/3)n3 O(n2L) A C T

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

下三角帯行列化のアルゴリズム（続き） 本アルゴリズムの特徴 [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におけるメモリ競合の影響を低減可能

下三角帯行列の二重対角化 対称帯行列の三重対角行列への相似変換（村田法） 村田法のアイディアの二重対角化への適用（Lang, 1996） 第 k 列の三重対角化 第 k 列の副対角要素より下を 0 にする鏡像変換を左と右からかける。 Bulge-chasing 上記の結果，帯より下に非零要素がはみ出すので，相似変換を繰り返して，はみ出した分を右下に移動し，右下隅から追い出す。 村田法のアイディアの二重対角化への適用（Lang, 1996） 第 k 列の二重対角化 第 k 列の副対角要素より下を 0 にする鏡像変換を左からかける。 上記の結果，上三角部分に非零要素がはみ出すので，左右から直交変換を繰り返し掛けることにより，はみ出した分を右下に移動し，右下隅から追い出す。

第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 ステップ　・・・ 　　が同時に実行可能 第1列のbulge-chasing における，右側からの 第3の直交変換で更新 される要素 第2列のbulge-chasing における，右側からの 第1の直交変換で更新 される要素 第1列による二重対角化は，今後 　より右の要素にのみ影響を及ぼす。 第1列の計算が右下まで行くのを待たずに，第2列の計算を開始できる。

４. 性能評価 評価環境 評価対象 評価方法 Opteron (1.8GHz), 1～4CPU ４.　性能評価 評価環境 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倍

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

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

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