密行列固有値解法の最近の発展（II） ー マルチシフトQR法 ー 名古屋大学　計算理工学専攻 山本有作 2006年1月28日 LA研究会 日立製作所の山本有作です。 「～」について発表いたします。

目次 １. はじめに ２. 非対称行列向けのマルチシフトQR法 ３. 対称行列向けのマルチシフトQR法 １.　はじめに ２.　非対称行列向けのマルチシフトQR法 ３.　対称行列向けのマルチシフトQR法 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１. はじめに 本研究で対象とする問題 応用分野 標準固有値問題 Ax = lx A ： n×n 密行列（実対称または実非対称） 対称行列 １.　はじめに 本研究で対象とする問題 標準固有値問題　Ax = lx A ： n×n 密行列（実対称または実非対称） 応用分野 対称行列 分子計算（分子軌道法，FMO-MO法） 統計計算 非対称行列 MHD，化学工学，量子化学，制御理論 Cf. Bai and Demmel: A test matrix collection for non-Hermitian eigenvalue problems. ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

固有値計算の流れ 密行列 A 計算内容 計算手法 QTAQ = T (Q: 直交行列) ハウスホルダー法 固有値の計算 QR法 分割統治法 三重対角化/Hessenberg化 三重対角行列 T / Hessenberg行列 H 固有値の計算 高い信頼性 QR法 分割統治法 MR3アルゴリズム I-SVDアルゴリズム A の固有値 {li } 固有ベクトルの計算 Tui =li ui 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 T or Hの固有ﾍﾞｸﾄﾙ {ui } 逆変換 vi = Q ui 逆変換 Aの固有ベクトル {vi }

各部分の実行時間（非対称行列） 全固有値を求める場合の演算量 Hessenberg 化： (10/3) n3 QR法： 10n3 （経験値） Execution time (min) Origin2000 1PU （R12000，400MHz）上での実行時間 K. Braman et al: “The Multishift QR Algorithm I” より抜粋 Hessenberg 化の時間は推定値

各部分の実行時間（対称行列） 全固有値を求める場合の演算量 三重対角化： (4/3) n3 QR法： O(n2) しかし，三重対角化部分は並列化が容易 QR法も並列化できることが望ましい n = 100,000の行列の固有値計算 5555 min. 8 min. 三重 対角化 10 min. 10 min. QR法 QR法は PrimePower HPC2500 上での実測値 三重対角化の時間は推定値 （1CPUでの性能4GFLOPS，並列化効率70%と仮定） 1CPU 三重対角化部分 のみ1000CPU

QR法の特徴と高性能計算 QR法の特徴 高性能計算の条件 計算の逐次性 （bulge-chasing 型演算） 計算の並列性 低いデータ再利用性 高いデータ再利用性 　　（キャッシュの有効利用） オリジナルのQR法のままでは高性能化が困難 高性能計算の条件を満たす新しいアルゴリズムが必要 マルチシフトQR法

２. 非対称行列向けのマルチシフトQR法 QR法（Francis, 1961 & Kublanovskaya, 1961） QR法の収束 A1 = R1Q1 (= Q1–1A0 Q1 ) A1 = Q2R2 A2 = R2Q2 (= Q2–1A1 Q2 = Q2–1Q1–1A0 Q1Q2 ) QR法の収束 適当な条件の下で，Ak は（ブロック）上三角行列に収束 A0 の固有値を絶対値の大きい順に l1, l2, … , ln とすると，対角要素 aij は li に１次収束 非対角要素 aij （j < i）は収束率 rij ≡ |li| / |lj| で 0 に１次収束

ダブルシフトQR法 シフトの導入 ダブルシフトQR法 Ak から固有値 li の近似値 s を引いた行列に対して QR 法を適用 Ak – s I = Qk Rk Ak+1 = Rk Qk + s I (= Qk–1Ak Qk ) 　　　　　　　　rij ≡ |li – s| / |lj – s| ≒ 0 より収束が加速 ダブルシフトQR法 共役複素数の固有値ペアに対し，シフト s, s による2反復をまとめて行う (Ak – s I)(Ak – s I) = Qk Rk Ak+2 = Qk–1Ak Qk シフトは右下隅の 2×2 行列の固有値に取る 　　　　→　局所的に2次収束 複素数の固有値を持つ場合でも，実数演算だけで QR 法を実行可能

Hessenberg 形と Implicit Q 定理 Ak が Hessenberg 形のとき， QR法の1ステップは O(n2) で実行可能 Ak+1 も Hessenberg 形 Implicit Q 定理の利用により，QR 分解を陽に行う 　 ことなく1ステップの計算を実行可能 　　　　　A0 を Hessenberg 形に相似変換してからQR法を適用 Implicit Q 定理 U，V が直交行列で，G = UtAU， H = VtAV が共に既約な Hessenberg 行列であるとする。 このとき，もし U と V の第1列が等しいならば，±1 を要素に持つ対角行列 D が存在して，V = UD，H = DGD が成り立つ すなわち，行列 A を既約な Hessenberg 行列に相似変換する直交行列は，第1列目だけが与えられれば（実質的に）一意に定まる Hessenberg 形 aij = 0 if i > j+1

Implicit shift QR法 1ステップの計算 この変換の性質 ある直交行列 H による相似変換 (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホル 　　ダー変換H0 を求める Ak’ = H0tAk H0 直交行列による相似変換を繰り返すことにより， Ak’ を再び 　　Hessenberg 行列に変形する（bulge-chasing） この変換の性質 ある直交行列 H による相似変換 H の第1列は，Qk の第1列に等しい 上記第3のステップは第2行・第2列以降のみに影響 H は Ak を Hessenberg 形に変換 Implicit Q 定理より，Ak+1 = H tAk H はQR法の1ステップと等価 Ak H0 AkH0 Ak+1

演算パターンとその特徴 Bulge-chasing の演算パターン 演算パターンの特徴 並列性は O(n) 並列アルゴリズム：　R. Suda et al.（1999），G. Henry et al. (2002) QR法の1ステップで，各行列要素は3回のみ更新 　　データ再利用性が低く，キャッシュの有効利用が困難 左から Hl を乗算 右から Hl を乗算 0にしたい要素 更新される要素

マルチシフトQR法 基本的なアイディア（Bai & Demmel, 1989） m 個のシフトに対するQR法のステップを同時に行う (Ak – smI) ・・・ (Ak – s2I)(Ak – s1I) = Qk Rk Ak+m = Qk–1Ak Qk シフト s1, s2, …, sm は，の右下隅の m×m 行列の固有値に取る マルチシフトQR法の収束（Watkins & Elsner, 1991） 上記のようにシフトを取ったとき，マルチシフトQR法は局所的に2次収束 特に，Q1, Qm+1, Q2m+1, … , Qk の積を Q(k) とするとき，Q(k) の最後の m 列の張る部分空間は，A1のある不変部分空間に2次収束する

マルチシフトQR法の演算パターン 1ステップの計算（implicit 型） 演算パターンの特徴 並列粒度は m/2 倍増加 (Ak – s2I) ・・・ (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホルダー変換H0 を求める Ak’ = H0tAk H0 直交行列による相似変換を繰り返すことにより， Ak’ を再びHessenberg 行列に変形する（bulge-chasing） 演算パターンの特徴 並列粒度は m/2 倍増加 QR法の1ステップで，各行列要素は m+1回更新 データ再利用性の向上により，キャッシュの有効利用が可能 　　（WY representation を用いた BLAS3化）

K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 一方，BLAS3 で性能を出すには m を数十程度に取ることが必要 　　　　このマルチシフトQR法では，計算機性能を引き出せない K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 収束性悪化の原因は？

収束性悪化の解析 シフトの伝播に関する理論的解析（Watkins, 1996） 実験からの考察 bulge を形作る (m+1)×(m+1) 行列を Bl とするとき，シフト 　　s1, s2, …, sm は，行列束 (Bk, N) の m 個の有限固有値 ただし N は固有値 0 の m+1 次ジョルダン細胞 実験からの考察 (Bl, N) の固有値は，m が大きいとき，摂動に対して極めて敏感になる Bl m = 4 の場合 Bulge chasing の過程でシフトが (Bk, N) の固有値として伝わる際に，丸め誤差による摂動でシフトが劣化し，2次収束を阻害するのではないか？

多数の小さい bulge を使ったマルチシフトQR法 基本的なアイディア（K. Braman et al., 2002） 従来法と同様にしてシフト s1, s2, …, sm を計算 シフトを2つずつ組にしてダブルシフトQR法を m/2 回実行 この結果得られる Ak+1 は，無限精度演算では従来法と同じ m/2 個の bulge をできるだけ密に詰めて chasing を行うことで，データ参照の局所性を向上 3行離れていれば，互いの干渉なしに計算が可能 m = 6 の場合 大きな bulge を避けることで，丸め誤差の影響を抑えつつデータの再利用性を高め，性能を向上

K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 収束性に関する実験結果 シフトの数と総演算量との関係 DHSEQR：　従来のマルチシフトQR法（LAPACK） TTQR：　新手法 K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 新手法では，m を数十まで増やしても収束性が低下しない

Level-3 BLAS の利用 更新処理の分割 r の決め方 m/2 個の bulge をそれぞれ r 行 chasing する際， 　 まず対角ブロックのみを更新 非対角ブロックは，更新に使ったハウスホルダー変換を1個の行列に蓄積し，後でまとめて更新 　　非対角ブロックの更新で level-3 BLAS が使用可能 r の決め方 演算量の点から，r ～ 3m とするのが最適 このとき，演算量は従来のマルチシフトQR法の約2倍（非ゼロ構造を利用した場合） Bulge（3×3） 最初に更新 まとめてBLAS3で更新 蓄積されたハウスホルダー変換の非ゼロ構造

従来のマルチシフトQR法との組み合わせ 基本的なアイディア（Kressner, 2005） 3×3の bulge を m/2 個同時に chasing する代わりに，(q+1)×(q+1) の bulge を m/q 個同時にchasing する q ～ 6 であれば，収束性は悪化しない 　　　　　　データ再利用性をさらに向上 Bulge（(q+1)×(q+1)） 最初に更新 まとめてBLAS3で更新

D. Kressner: “On the Use of Larger Bulges in the QR Algorithm”より引用 各手法の性能 比較する対象の手法 DHSEQR：　従来のマルチシフトQR法（LAPACK） MTTQR：　多数の小さいシフトを使ったマルチシフトQR法 　　　　　　　　　（ns： bulge 1個あたりのシフト数） D. Kressner: “On the Use of Larger Bulges in the QR Algorithm”より引用 ほとんどの場合，ns = 4 or 6 とした新手法が最高速

共有メモリ向けの並列化手法 Level-3 BLAS 内部での並列化 より効率的な並列化 非対角ブロックの更新において，次の対角ブロックの更新に必要な部分のみを先に更新 　　並列化困難な対角ブロックの更新を，非対角ブロックの更新と並列に実行可能 次の対角ブロックの 更新で必要 Bulge（3×3） 最初に更新 まとめてBLAS3で更新

シフト数の最適決定 m の最適化 r （1回の bulge chasing の行数）の最適化 m，r に対する自動チューニングの必要性 演算量最小化のためには r ～ 3m が最適 BLAS3 の実行時間は必ずしも演算量に比例しないため，実行時間の最適値はこれとは異なる 　　　　m，r に対する自動チューニングの必要性 n m の最適値 1000～1999 60 2000～2499 116 2500～3999 150 4000～ 180 実験的に求めた m の最適値 （Origin2000上，Braman et al. より引用）