Transcript (I) --- Multiple Relatively Robust Representation アルゴリズム

密行列固有値解法の最近の発展 (I)
－ Multiple Relatively Robust Representation
アルゴリズム －
2004年11月26日
名古屋大学 計算理工学専攻
山本有作
目次
１. はじめに
２. Multiple Relatively Robust Representation アルゴリズム
３. 対称密行列の固有値計算への適用
１. はじめに
• 本報告で対象とする問題
– 標準固有値問題 Au = λu
• A： 実対称またはエルミートの n×n 密行列
• 全部または一部の固有値・固有ベクトルを求める。
• 応用分野
– 分子計算，統計計算，構造解析など
固有値・固有ベクトル計算の流れ
実対称行列 A
計算内容
三重対角化
Q*AQ = T
(Q: 直交行列)
三重対角行列 T
｜T – λi I｜= 0
二分法
Tの固有値 {λi }，
逆反復法
Tvi =λi vi
Tの固有ベクトル{vi }
逆変換
Aの固有ベクトル {ui }
ui = Qvi
逆反復法による固有ベクトル計算
• 逆反復法の原理
– T: 三重対角行列
– λ’i : Tの固有値λi の近似値
– 適当な初期ベクトル vi (0) から出発し，次の反復を行う。
vi (m) := ( T – λ’i I ) –1 vi (m – 1)
– vi に平行な成分が，１反復毎に (λi–λi’) –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 が次の
条件を満たすように求められたとする。
|| Tx – xλ’|| = O(ne) |λ’|
--- (*)
– このとき，sin theorem より
sin|∠(x, v)| ≦ || Tx – xλ’|| / gap(λ)
= O(ne) |λ’| / gap(λ)
～ O(ne) / relgap(λ)
– ここで，relgap(λ) = gap(λ) / |λ|。
(*)式の成立を保証できれば，固有値の相対ギャップが大きい場
合には直交化なしで自動的に精度の高い（したがって直交性も
良い）固有ベクトルが求まる。
基本的なアイディア（続き）
• 従来のアルゴリズムの問題点
– 従来の二分法・逆反復法では，次の不等式しか成り立たない。
|| Tx – xλ’|| = O(ne) ||T||
– 小さい固有値に対しては，相対残差が大きくなる可能性がある。
• 新しいアルゴリズムの概要（相対ギャップが大きい場合）
– (1) T +μI が正定値となるようにμを選び，T +μI = LDLT と改訂
コレスキー分解を行う。
– (2) LDLTの固有値の近似値λ’を，相対誤差の意味で高精度に
計算する （dqds法などを利用）。
– (3) twisted 分解を用いて，λ’に対する固有ベクトルを相対残差
が小さくなるよう高精度に計算する。
固有値の相対ギャップが大きい場合
• なぜ分解 T +μI = LDLT が必要か
– 計算した固有値λ’の誤差は，通常，後退誤差解析 + 摂動論により
評価する。
後退誤差解析： λ’はあるdT に対して T+dT の厳密な固有値
摂動論：
T → T+dT のとき，固有値はdλだけずれる。
– しかし，三重対角行列 T に対しては，dλが ||dT|| でしか押さえられ
ない 。
– 相対誤差の意味で高精度とするには， dλを ||λdT||で押さえたい。
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法を用いる）。
LDL
T
 1

1
 L




λ I 








1
k 2
L
1
k 1
L
1
U
1
k

U
k 1

1
  D 1










n 1 
U 

1  
2
D
k 1
D
γk
k 1
D
n 1
D
 1












n 
D  
U
1

1
U
2

1
k
L
U
k 2

1
k 1
L
U
k 1

1
– このうち，γkが最小になるような k を求め，(LDLT –λ’I )x = γkek を
（上式の右辺を用いて）解く。
1
n 1
L













1 
固有ベクトルの高精度計算（続き）
• Twisted分解（続き）
– このとき，得られた解ベクトル x は次の式を満たすことが示せる。
（Dhillon, 1997）
|| (LDLT –λ’I ) x || / ||x|| ≦ n |λ–λ’ | ・M / (M – 1)
– ただし，Mはある正の定数。
– λ’が相対誤差の意味で高精度（ |λ–λ’ | = O(e) |λ’| ）ならば，
|| Tx – xλ’|| = O(ne) |λ’| が言える。
– 固有ベクトルの近似値 x は高精度。
固有値の相対ギャップが小さい場合
• 問題点
– 以上のアルゴリズムで言えるのは
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) の計算量
– 高い並列性
分散メモリ型並列計算機上で
高い性能
三重対角化と逆変換の時間が
相対的に増大
Pentium 4クラスタ（16PU）上での性能
（Dhillon, 2004）
三重対角化のための高速アルゴリズム
• Dongarra のアルゴリズム
–
–
–
–
ハウスホルダー法におけるrank-2更新を多段化
Level-3 BLAS で書けるのは全演算量の1/2のみ
キャッシュマシンではピークの10～25%の性能
通信回数が多い （各ステップで通信）
• Bischof / Wu のアルゴリズム
– 行列をいったん帯行列に変換し，村田法により三重対角化
– 全演算量のほとんどを level-3 BLAS で実行可能
– 通信回数が少ない （Dongarra のアルゴリズムの 1/L）
次数 N
半帯幅 L
帯行列化
約 (4/3)N3
A
0
村田法
O(N2L)
0
B
0
0
T
各アルゴリズムの性能（Opteron, 1.6GHz）
1800
Performance (GFLOPS)
1600
1400
L=24, L’=4
Dongarra
Bishof
Wu
L=48
1200
1000
800
L’=32
600
400
200
0
480
960
1920
3840
Matrix size
Wu の方法は Dongarra の方法に比べて約２倍の性能を達成
N = 3840 のとき，Wu の方法はピークの50%以上の性能を達成
Bischof / Wu のアルゴリズムでの固有ベクトル計算
（従来の逆反復法，直交化が必要ない場合）
• 計算法１
– 三重対角行列に対して逆反復法を行
い，得られる固有ベクトルに２段階の
逆変換を行う。
A
B
0
T
0
0
0
O(kn)
{λi }
{λi }
2kn2
{ui }
2kn2
{wi }
{vi }
O(kn)
• 計算法２
– 三重対角行列の固有値を用いて帯
行列に対して逆反復法を行い，１段
階の逆変換を行う。
A
B
0
0
0
0
O(kLn2L2)
{λi }
{ui }
O(kn)
{λi }
2kn2
L2 ≪ n ならば計算法２のほうが高速
T
{wi }
計算法２にMR3アルゴリズムを適用する際の問題点
• 固有値の相対精度の問題
– B → T → {λi } という経路で求めた固有値は，相対誤差の意味で
B の高精度な固有値になっていない。
• T の高精度な固有値には当然なっている。
• 三重対角化アルゴリズム（村田法）の問題ではなく，三重対角行列
への変形自体が相対精度を破壊すると思われる。
– Twisted 分解による固有ベクトル計算アルゴリズム （の拡張） を
適用するための前提が成り立たない。
解決策
• （案１） 計算法１を用いる。
– 三重対角行列 T の固有値・固有ベクトルをMR3で計算
– 固有ベクトルを２段階に逆変換 （2kn2 + 2kn2）
• （案２） MR3アルゴリズム全体を帯行列に拡張
– Twisted 分解，dqds法等を帯行列に対して拡張（可能か？）
– L2 ≪ n かつ k ～ n ならば，案１より高速になると予想される。
– 帯行列に対して適用することで，dqds 法の収束性を三重対角行
列の場合より向上できる可能性 → 更なる高速化