対称行列の固有値問題に対する Sylvester 二分法の精度保証
2009SE055 平野 由莉
指導教員：杉浦洋

1

はじめに2

Sylbester の慣性則

修正 Cholesky 分解は，dj = 0 となると式 (1) での割り算が生じて破綻する．また，| dj | が小さいと式 (1) で本研究では，対称三重対角行列式の固有値を Slvester| l | が大きくなり，後で述べる定理により固有値の包囲二分法により求め，後退誤差解析により精度保証を行う． ij
区間幅が大きくなる．これを改善する方法が対角ピボップログラムを作成し，数値実験により，その効果を実証ティングである．//
し，Slvester 二分法による精度保証の特性を分析する．
3.1 対角ピボッティング

定理 2.1(Sylveter[1])対称行列 A を合同変換 A → P AP T = D により，対角行列 D に変換し，D の正の対角成分の数を p，負の対角成分の数を m とする．このとき，p, m は合同変換によらず一定である．//
この p, m を A の符号数という．
この定理からただちに次の定理が導かれる．
系 2.1 対称行列の正の固有値の数，負の固有値の数は，それぞれ A の符号数 p, m と等しい．//3

正則な対角行列 A の修正コレスキー分解は，

A = LDLT ,

10...
0

10
 l21

,
L= .
..
..
.
.
0
 ..
ln1 . . . ln,n−1 1

d1 0 . . . 0.. ..
0 d.
.


2D=.
..
 . ..
.
. 0
.
0 ...
0 dnとすると，A − λI の固有値は λi − λ(1 ≤ i ≤ n) であるから，

p(λ) : λi − λ > 0 となるλi の数.

ck ljk ;k=1djである．目標の固有値を λi とすると
{p(λ) < i =⇒ λ ≥ λip(λ) ≥ i =⇒ λ < λiである．これを用いて，Sylvester 二分法のアルゴリズムにより，任意の精度で λi が包囲できる．すなわち，A の固有値を大きい順に λ1 > λ2 > · · · > λn とし，第 i 番の固有値 λi を求める Sylvester 二分法のアルゴリズムは以下のようになる [3]．

(0) λi を含む初期区間を [a0 , b0 ) とする．
許容誤差 ϵ > 0 を与える． k = 0 とする．

i = j + 1, j + 2, . . . , n の順で，lij =Sylvester 二分法のアルゴリズム

λ1 ≥ λ2 ≥ · · · ≥ λnk=1j−1∑

P AP T = LDLT与えられた λ に対して，A − λI の正の符号数を p(λ) とする．A の固有値をck = dk ljk ;aij −, U = DLTと解釈すれば，A の LU 分解と考えられる．
対角ピボッティングは，Q = P とし，P AQT を A の合同変換とする．これにより，P AP T は対称となり，4である [2]．L が正則なら A → L−1 AL−T = D は合同変換ゆえ，修正 Cholesky 分解で A の符号数がわかる．この L，D を求めるアルゴリズムを示す．
j = 1, 2, 3, . . . , n の順に，
k = 1, 2, . . . , j − 1 の順で，j−1∑

A = LDLT = L(DLT ) = LUと修正 Cholesky 分解される．
対角ピボッティングにより，すべての対角成分が 0 にならないかぎり，0 ピボットが発生せず，アルゴリズムが頑健になる．また，P は置換行列であるため，合同変換は丸め誤差を発生しない．//

修正 Cholesky 分解

dj = aii −0 ピボットを回避するための対角ピボッティングについて説明する．
対称行列 A の修正 Cholesky 分解を

ck ljk;

この計算量は， 16 n3 + O(n2 )flops である．//(1)(1) c =

ak + bkとし，2{p(c) < i =⇒ [ak+1 , bk+1 ) = [c, bk ),p(c) ≥ i =⇒ [ak+1 , bk+1 ) = [ak , c).

(2) bk+1 − ak+1 ≥ ϵ ならば，k → k + 1 として (1) へ．

表 1 Sylvester 二分法による結果(3) 区間 [ak+1 , bk+1 ) ∋ λi を出力．

5Sylvester 二分法の精度保証

5.1

修正 Cholesky 法の後退誤差解析

k12345678910

A − λI に対する修正 Cholesky のアルゴリズムは次と同値
j = 1, 2, 3, . . . , n の順に，
dj = ajj − λ −j−1∑
k=12dk ljk;

i = j + 1, j + 2, . . . , n の順で，
aij −lij =j−1∑

λk4.47 · · · × 101
5.04 · · ·1.87 · · ·1.00 · · ·6.43 · · · × 10−1
4.65 · · · × 10−1
3.66 · · · × 10−1
3.07 · · · × 10−1
2.73 · · · × 10−1
2.55 · · · × 10−1

|mid[λk ]-λk |1.56 × 10−13
1.78 × 10−15
5.33 × 10−15
4.66 × 10−15
2.11 × 10−15
2.16 × 10−15
1.67 × 10−16
1.72 × 10−15
7.22 × 10−16
2.22 × 10−16rad[λk ]3.41 × 10−13
1.84 × 10−13
9.77 × 10−14
8.64 × 10−14
6.92 × 10−14
6.98 × 10−14
6.81 × 10−14
6.83 × 10−14
6.82 × 10−14
6.82 × 10−14

dk lik ljk

k=1dj; //6

これを機械計算した結果を ˜
lij , d˜j と書き，数値実験によく用いられる Frank 行列

A = (aij ) ∈ Rn×n ,L̃ = (˜lij ), D̃ = diag(d˜j ),Ã = L̃D̃L̃T , ∆A = Ã − Aとすると，∥ ∆A ∥∞ ≤∥| L̃ || D̃ || L̃T |∥∞ γn+1 .

(2)

ここで，行列 B = (bij ) のとき | B |= (| bij |) を表す．また，γk = ku/(1 − ku) であり，u = 2−53 は丸め誤差である．式 (2) は，次のように変形でき，右辺の計算量はO(n2 ) である．

∥ ∆A ∥∞ ≤∥| L̃ || D̃ || L̃T | e ∥∞ γn+1 .

数値実験

aij = min{i, j}を取り上げる．また，次数は n = 10 とする．そして，結果を表 1 に示す．表 1 で”k”は固有値番号，”λk ”は固有値，λI は Sylvester 二分法で計算された λk の包囲区間であり，”mid[λI]”はその中心，”rad[λI]”はその半径である．すべての固有値の包囲に成功した．
対角ピボッティングを用いない修正 Cholesky 法でも，すべての固有値の包囲に成功した．しかし，首座小行列の固有値と一致する λ2 , λ4 , λ5 , λ8 では，包囲区間が極端に広くなり，精度保証の品質が劣化した．これは，| L̃ |のノルムが大きくなることが原因である．//(3)

7

ここで，e = (1, 1, · · · , 1)T である．//

おわりに本研究では，対称行列の固有値を求める Slvester 二分法により，固有値を精度保証付きで計算するアルゴリズ以上の結果より次の定理を得る．ここで，A の大きい ムを作成し，Mathematica 上に実装した．
1．修正 Cholesky 分解の後退誤差解析を研究した．
方から i 番目の固有値を λi (A) と書く．
定理 5.1 機械計算による A − λI の桁数を p̃(λ) とする． 2．Sylvester 二分法の出力を 1 により補正し，固有値の厳格な包囲範囲を計算した．
このとき，次の不等式が成立する．
3．修正 Cholesky 法が 0 ピボットのため破綻する条件を明らかにした．
λi (A) ≥ λ − δ(λ),
4．0 ピボットを回避するため，対角ピボッティング付き
λi−1 (A) < λ + δ(λ).
修正 Cholesky 法を作成した．
5．数値実験により，2，4 の効果を確認した．
ここで，今後の課題としては，さらに頑健な修正 Cholesky 分解を
δ(λ) = γn+1 ∥| L̃ || D̃ || L̃T | e ∥∞ .//構成することである．

5.2Slvester 二分法の精度保証

定理 5.2 Sylvester 二分法の機械計算で得られた λi (A)の「包囲区間」を (a, b) とする．このとき，λi (A) ∈ [a − δ(a), b + δ(b)).//

(4)

参考文献
『線形代数学』．東海大学出版会．
[1] 東海大学数学教室：[2] 杉浦洋：『数値計算の基礎と応用—数値解析学への入門—』．サイエンス社，1997．
[3] W.Kahan：ACCURATE EIGENVALUES OF A
SYMMETRIC TRI-DIAGONAL MATRIX，CS41，
Computer Science Dpt，Stanford University，1966．