密行列固有値解法の最近の発展 (I) －　Multiple Relatively Robust Representation アルゴリズム　－ 2004年11月26日 名古屋大学　計算理工学専攻 山本有作 日立製作所の山本有作です。 「～」について発表いたします。

目次 １. はじめに ２. Multiple Relatively Robust Representation アルゴリズム １.　はじめに ２.　Multiple Relatively Robust Representation アルゴリズム ３.　対称密行列の固有値計算への適用 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１. はじめに 本報告で対象とする問題 応用分野 標準固有値問題 Au = λu A： 実対称またはエルミートの n×n 密行列 １.　はじめに 本報告で対象とする問題 標準固有値問題　Au = λu A： 実対称またはエルミートの n×n 密行列 全部または一部の固有値・固有ベクトルを求める。 応用分野 分子計算，統計計算，構造解析など ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

固有値・固有ベクトル計算の流れ Q*AQ = T (Q: 直交行列) ｜T – λi I｜= 0 Tvi =λi vi ui = Qvi 計算内容 Q*AQ = T (Q: 直交行列) 三重対角化 三重対角行列 T 二分法 ｜T – λi I｜= 0 Tの固有値 {λi }， 逆反復法 Tvi =λi vi 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 Tの固有ベクトル{vi } 逆変換 ui = Qvi Aの固有ベクトル {ui }

逆反復法による固有ベクトル計算 逆反復法の原理 vi (m) := ( T – λ’i I ) –1 vi (m – 1) λ’i : Tの固有値λi の近似値 適当な初期ベクトル vi (0) から出発し，次の反復を行う。 vi に平行な成分が，１反復毎に (λi–λi’) –1 倍に拡大される。 vi (m) := ( T – λ’i I ) –1 vi (m – 1)

逆反復法の長所と短所 長所 短所 一部の固有ベクトルのみの計算が可能 固有値が十分に離れている場合，k 本の固有ベクトルを計算するための計算量は O(kn) 短所 固有値が密集している場合，固有ベクトルの直交化が必要 固有ベクトルを全部直交化する場合，計算量は O(k2n) に増加 大規模問題（n > 1000）ではほとんど常にこの状況 直交化が必要な場合，並列化が困難（不可能ではないが） 　　直交化を行わずに高精度な固有ベクトルを求める方法ができれば， 　　計算量と並列性の面で非常に有利

本報告の目的 直交化を行わずに三重対角行列の高精度な固有ベクトルを計算する方法である Multiple Relatively Robust Representation アルゴリズム （MR3 アルゴリズム，Dhillon (1997)） について，概要を紹介する。 対称密行列の固有値計算に MR3 アルゴリズムを適用する際の課題について考察する。

２. Multiple Relatively Robust Representation アルゴリズム 基本的なアイディア 固有値の相対ギャップが大きい場合 固有値の相対ギャップが小さい場合 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

基本的なアイディア 固有ベクトルに関する sin theorem T を対称な三重対角行列，λ’を固有値の近似値，λをλ’にもっとも近い固有値とする。このとき，長さ１の任意のベクトル x に対して次の不等式が成り立つ。 　　　　　　　sin|∠(x, v)| 　≦　|| Tx – xλ’|| / gap(λ) ここで，gap(λ) = |μ–λ|，μはλ以外で最もλ’に近い固有値。

基本的なアイディア（続き） sin theorem の利用 いま，固有値の近似値λ’と固有ベクトルの近似ベクトル x が次の条件を満たすように求められたとする。 このとき，sin theorem より ここで，relgap(λ) = gap(λ) / |λ|。 　(*)式の成立を保証できれば，固有値の相対ギャップが大きい場合には直交化なしで自動的に精度の高い（したがって直交性も良い）固有ベクトルが求まる。 || Tx – xλ’||　=　O(ne) |λ’|　　　---　(*) sin|∠(x, v)| 　≦　|| Tx – xλ’|| / gap(λ) 　　　　 =　 O(ne) |λ’| / gap(λ) 　　　　～　O(ne) / relgap(λ)

基本的なアイディア（続き） 従来のアルゴリズムの問題点 新しいアルゴリズムの概要（相対ギャップが大きい場合） 従来の二分法・逆反復法では，次の不等式しか成り立たない。 小さい固有値に対しては，相対残差が大きくなる可能性がある。 新しいアルゴリズムの概要（相対ギャップが大きい場合） (1) T +μI が正定値となるようにμを選び，T +μI = LDLT と改訂 　　　 コレスキー分解を行う。 (2) LDLTの固有値の近似値λ’を，相対誤差の意味で高精度に 　　　 計算する （dqds法などを利用）。 (3) twisted 分解を用いて，λ’に対する固有ベクトルを相対残差 　　　 が小さくなるよう高精度に計算する。 || Tx – xλ’||　=　O(ne) ||T||

固有値の相対ギャップが大きい場合 なぜ分解 T +μI = LDLT が必要か 計算した固有値λ’の誤差は，通常，後退誤差解析 + 摂動論により評価する。 しかし，三重対角行列 T に対しては，dλが ||dT|| でしか押さえられない 。 相対誤差の意味で高精度とするには， dλを ||λdT||で押さえたい。 後退誤差解析： λ’はあるdT に対して T+dT の厳密な固有値 摂動論： T → T+dT のとき，固有値はdλだけずれる。 LDLT の形で表現された行列の固有値問題（すなわち LD1/2 の特異値問題）に対しては，dλを ||λd(LD1/2 )||で押さえることが可能（Kahan, 1967）　→　Relatively Robust Representation

LDLT の固有値の高精度計算 特異値分解アルゴリズムの利用 二重対角行列に対しては，その特異値を相対誤差の意味で高精度に計算するアルゴリズムが存在 二分法の改良 （Kahan, 1967） dqdsアルゴリズム （Fernando & Parlett, 1994） これを LD1/2 に適用することにより， LDLT の固有値λを相対誤差の意味で高精度に計算可能

固有ベクトルの高精度計算 Twisted分解 逆反復法の良い初期ベクトルを求めるための手法 近似固有値λ’に対し，LDLT –λ’I を各 k （1 ≦ k ≦ n）に対して次のように分解 （計算にはdqds法を用いる）。 このうち，γkが最小になるような k を求め，(LDLT –λ’I )x = γkek を（上式の右辺を用いて）解く。

|| (LDLT –λ’I ) x || / ||x|| ≦ n |λ–λ’ | ・M / (M – 1) 固有ベクトルの高精度計算（続き） Twisted分解（続き） このとき，得られた解ベクトル x は次の式を満たすことが示せる。（Dhillon, 1997） ただし，Mはある正の定数。 λ’が相対誤差の意味で高精度（ |λ–λ’ | = O(e) |λ’| ）ならば， 　　|| Tx – xλ’||　=　O(ne) |λ’|　が言える。 固有ベクトルの近似値 x は高精度。 || (LDLT –λ’I ) x || / ||x||　≦　n |λ–λ’ | ・M / (M – 1)

固有値の相対ギャップが小さい場合 問題点 行列のシフトの利用 以上のアルゴリズムで言えるのは 　 sin|∠(x, v)| ≦ O(ne) / relgap(λ)まで。 relgap(λ)が大きい場合は，固有ベクトルの高精度性が言えない。 行列のシフトの利用 T の固有ベクトルと T –νI の固有ベクトルは共通。 ν～λと取れば，relgap(λ)は大きくできる。 既約な三重対角行列に重複固有値は存在しない。 上記の変形を行った上で，相対ギャップが大きい場合のアルゴリズムを適用。

固有値の相対ギャップが小さい場合（続き） 課題１ T –νI は一般に正定値行列ではない。 LDLT分解は可能だが，それが Relatively Robust Representation である （固有値を相対誤差の意味で高精度で決定する） とは一般に言えない。 Dhillon (1997) では，「証明はできないが，数値実験の結果では，ほとんどの場合， R3 を与えるνがλの近くに存在」と主張。 課題２ 異なる固有値に属する固有ベクトルの計算には，複数の R3 が 必要（MR3）。これらの間の変形を高精度にできるか？ この変形にも dqdsアルゴリズムを使うことを提案。

３. 対称密行列の固有値計算への適用 MR3アルゴリズムの性能 O(kn) の計算量 高い並列性 分散メモリ型並列計算機上で高い性能 ３.　対称密行列の固有値計算への適用 MR3アルゴリズムの性能 O(kn) の計算量 高い並列性 分散メモリ型並列計算機上で高い性能 三重対角化と逆変換の時間が相対的に増大 Pentium 4クラスタ（16PU）上での性能 （Dhillon, 2004）

三重対角化のための高速アルゴリズム Dongarra のアルゴリズム Bischof / Wu のアルゴリズム ハウスホルダー法におけるrank-2更新を多段化 Level-3 BLAS で書けるのは全演算量の1/2のみ キャッシュマシンではピークの10～25%の性能 通信回数が多い （各ステップで通信） Bischof / Wu のアルゴリズム 行列をいったん帯行列に変換し，村田法により三重対角化 全演算量のほとんどを level-3 BLAS で実行可能 通信回数が少ない （Dongarra のアルゴリズムの 1/L） 半帯幅 L A B T 次数 N 約 (4/3)N3 O(N2L) 帯行列化 村田法

各アルゴリズムの性能（Opteron, 1.6GHz） L=24, L’=4 L=48 Performance (GFLOPS) L’=32 Matrix size Wu の方法は Dongarra の方法に比べて約２倍の性能を達成 N = 3840 のとき，Wu の方法はピークの50%以上の性能を達成

Bischof / Wu のアルゴリズムでの固有ベクトル計算 （従来の逆反復法，直交化が必要ない場合） 計算法１ 三重対角行列に対して逆反復法を行い，得られる固有ベクトルに２段階の逆変換を行う。 計算法２ 三重対角行列の固有値を用いて帯行列に対して逆反復法を行い，１段階の逆変換を行う。 T A B O(kn) {λi } {λi } 2kn2 2kn2 {ui } {wi } {vi } O(kn) T A B O(kn) O(kLn2L2) {λi } {λi } 2kn2 {ui } {wi } L2 ≪ n ならば計算法２のほうが高速

計算法２にMR3アルゴリズムを適用する際の問題点 固有値の相対精度の問題 B → T → {λi } という経路で求めた固有値は，相対誤差の意味で B の高精度な固有値になっていない。 T の高精度な固有値には当然なっている。 三重対角化アルゴリズム（村田法）の問題ではなく，三重対角行列への変形自体が相対精度を破壊すると思われる。 Twisted 分解による固有ベクトル計算アルゴリズム （の拡張） を適用するための前提が成り立たない。

解決策 （案１） 計算法１を用いる。 （案２） MR3アルゴリズム全体を帯行列に拡張 三重対角行列 T の固有値・固有ベクトルをMR3で計算 固有ベクトルを２段階に逆変換 （2kn2 + 2kn2） （案２） MR3アルゴリズム全体を帯行列に拡張 Twisted 分解，dqds法等を帯行列に対して拡張（可能か？） L2 ≪ n かつ k ～ n ならば，案１より高速になると予想される。 帯行列に対して適用することで，dqds 法の収束性を三重対角行列の場合より向上できる可能性　→　更なる高速化