Kahan の方法による三重対角行列の固有値の精度保証
2010SE001 ： 赤地祐暉
指導教員：杉浦洋

1

はじめに

固有値は線形代数の基本的な概念で，数理科学において多くの応用をもつ．それらの固有値の多くは数値計算によって求められる．固有値の分布が解析対象の特性を決定付けることも多いので，その精度保証は非常に重要である．渡部による先行研究 [1] では，Kahan の符号数計算アルゴリズムにおける 0 割り算の回避が大きな問題となった. 渡部は近似固有値のシフトによる 0 割り算回避を行っているが，最小限のシフトで回避するために複雑なプログラムと大きな計算量を要した．今回はその改良を目指す.

2

Sylvester の慣性則0..

.

0

dnに変換し d1 , · · · , dn の中で正数の数を p 個, 負数の数をm 個とする．このとき，p, m は合同変換によらず一定である．//
上の (p, m) を A の符号数と呼ぶ．

[系] 対称行列 A の符号定数を (p,m) とすると，p は Aの正の固有値の数，m は A の負の固有値の数である．
(証明)A の固有値を λ1 , λ2 , · · · , λn 、対応する固有ベクト
ルを u1 , u2 , · · · , un とする．
λi を対角成分とする対角行列 Λ, ui を列ベクトルとする行列を U = (u1 , u2 , · · · , un ) とすると，とかける．これは対角行列 Λ から A への合同変換である．ゆえに Sylvester の慣性則より，正の固有値の数は p,負の固有値の数は m である．//

..
..

.

.
bn−1
bn−1 anのとき L は二重対角行列1l1L=..
..

d1,D = .
.

..

.
ln−1d2..

.

1dnd1 = a 1 ，bi−12, di = ai − di−1 li−1
(i = 2, 3, · · · , n)li−1 =di−1が得られる．
この di の符号を調べることにより A の符号数がわかる．
符号数の計算には l は不要ゆえに di のみを計算すると)2(b2bi−1= ai − i−1(2)di = ai − di−1di−1di−1となる．これが符号数計算のための Kahan の符号数計算法である．//

4Kahan の方法与えられた λ に対して，A − λI の符号数 (p(λ)，m(λ)) とする．A の固有値をとすると，A − λI の固有値は λi − λ(1 ≤ i ≤ n) であるから，p(λ) : λi − λ > 0 となるλi の数.
m(λ) : λi − λ < 0 となるλi の数.
である．符号数は Kahan の計算法で計算する.

Kahan の符号数計算法
修正 Cholesky 分解は n 次対称行列 A を

A = LDLTa2..
.

λ1 ≤ λ2 ≤ · · · ≤ λn

A = U ΛU T3

b 1A=

b1である．
(1) の両辺を成分ごとに比較すると
a1 = d1 ，2bi−1 = li−1 di−1 , ai = di−1 li−1+ di (2 ≤ i ≤ n)を得る.
これより L, D の成分を計算するアルゴリズム[定理] A を合同変換により対角行列
d1
D = diag(d1 , · · · , dn ) = 

a1(1)と単位下三角行列 L と対角行列 D および LT の積に分解することである．
A を対称三重対角行列以下のグラフ (図 1) より，目標の固有値を λi とすると
{m(λ) < i =⇒ λ ≤ λim(λ) ≥ i =⇒ λ > λiである．これを用いて，二分法のアルゴリズムにより，任意の精度で λi が包囲できる．すなわち，A の固有値を小さい順にλ1 < λ2 < · · · < λn とし，第 i 番の固有値 λi を求める Kahanの方法は以下のようになる．

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

として計算を実行する.
このとき,Kahan のアルゴリズムは

2. c = (ak + bk )/2 とし，{m(c) < i =⇒ [ak+1 , bk+1 ) = [c, bk ),m(c) ≥ i =⇒ [ak+1 , bk+1 ) = [ak , c).
である.

3. bk+1 − ak+1 ≥ ϵ ならば，k → k + 1 として 2. へもどるd1A‘ = A + D, D = diag(δ1 , δ2 · · · , δn ) = Kahan のアルゴリズムの後退誤差解析

Kahan のアルゴリズムの後退誤差解析を行う．符号数を計算する修正 Cholesky 分解のアルゴリズムは，
d1 = a1 − λ(3).

dnに適用されたことになる. ゆえに，{m(λ) < i =⇒ λ ≤ λi (A′ )m(λ) ≥ i =⇒ λ > λi (A′ )(9)D の固有値は δ1 , δ2 · · · , δn であるから，である. さて，次にノルムの評価を行う．
∆αi = α̃i − αi , ∆βi = β̃i − βidi = αi − λ −..

04. 区間 [ak+1 , bk+1 ) ∋ λi を出力．終了．

5

δ2||D||2 = δ = max |δi |(10)

1≤i≤n

2βi−1
di − 10であるから，定理 3.4 より，(4)|λk (A) − λ(A′ )| ≤ δ(1 ≤ k ≤ n)である．ここで(11)である.(2) と (4) より，{{m(λ) < i =⇒ λ − δ ≤ λi (A)m(λ) ≥ i =⇒ λ + δ > λi (A)δ(λ) = max |α1 |u + |λ|u + |β1 |µ3 ,{}max|αi |γ2 + |λ|γ2 + (|βi−1 | + |βi |)µ3 ,2≤i≤n−1}|βn−1 |µ3 + |αn |γ2 + |λ|γ2であり,∥∆A∥∞ ≤ δ(λ) が成立する．
よって A の固有値で小さい方から i 番目を λi (A) と書くと，|λi (Ã) − λi (A)| ≤ ∥∆A∥∞ ≤ δ(λ)

を得るので，二分法の区間縮小を{[ak , c + δ)(m(c) ≥ i)
[ak+1 , bk+1 ) =[c − δ, bk )(m(c) < i)7である．また，δ(λ) は |λ| の単調増加関数であることは明らかだから，|λ| = max |a|, |b|として，厳格な不等式

a − δ(λ) ≤ λm (A) < b + δ(λ)(6)を得る．

6Mathematica による数値実験
テスト行列として
2−1
 −1 2
−1
A=

a − δ(a) ≤ λm (A)である．同じく，λ = b のときの後退誤差解析の結果は，式 (5)より
λm (A) < b + δ(b)(13)に変更する. c + δ, c − δ は精度保障付で計算する.

(5)が成立する.
後退誤差解析の結果 (5) を用いて，4.4 節のアルゴリズムの結果に精度保証を与えることが出来る．4.4 節のアルゴリズムで得られた小さい順で i 番目の固有値 λi の包囲区間を (a, b) とする．λ = a のときの後退誤差解析の結果は，式 (5) より(12)0−12..
.

0−1..
.
−1..

.
2−1 ∈ Rn×n−1 2を用い,n = 2048 で全固有値の包囲に成功した．
n = 128 の場合，最終区間幅の最大値は Max Width= 4.88498×10−15 で 2012 年の渡部優斗の卒業論文 [1] における結果 Max
width= 4.88498 × 10−15 と同じ結果が得られた．
n = 2048 の場合も同等に，最終区間幅の最大値は n = 128 のときと同じで，4.88498 × 10−15 と大次元にも関わらず極めて精密に精度保証ができた．

8

考察

Kahan の方法が 0 割り算で破綻するのを回避する単純で計算量の少ないアルゴリズムを提案した. 結果は渡部 [2] と同等で，非常に精密な精度保証に成功した.

0 割り算回避法

Kahan のアルゴリズムの (2) で,di = 0 となると，次の段階で 0 割り算のため計算が破綻する.
渡部は，λ を変更し，(2) の計算を最初からやり直した.,λ の変 9
おわりに更を最小限にするため，二分法を用いて 0 割り算を起こさない，Kahan の符号数計算アルゴリズムにおける零割り算の回避がλ を探索した. この方法は，計算量が大きくなり，制御も複雑で大きな問題となっていたが，Kahan の方法が渡部による先行研ある.
究 [1] よりも単純で計算量の少ないアルゴリズムが提案でき，す
今回提案するのは，αi を αi + δi (δi ̸= 0) に変更する方法であべての固有値に対して非常に精密な精度保証を求めることに成る. di ̸= 0 ならば，δi = 0,功した.
di = 0 ならば，
δi = u min{1, |ai |}(7)[1] 渡部優斗，Kahan による三重対角行列の固有値の精度保証，
南山大学情報理工学部システム創成学科 2012 年度卒業論文(2013).とし,(5.2) を

di = δi (= αi + δi − λ −2
βi−1)
di − 1

参考文献(8)