Kahan の方法による三重対角行列の固有値の精度保証
2009SE179 渡部優斗
指導教員：杉浦洋

1

はじめに3Kahan のアルゴリズム

本研究では，Kahan の二分法による対称三重対角行列与えられた λ に対して，A−λI の符号数を (p(λ)，m(λ))
式の近似固有値の精度保証を行う．後退誤差解析により， とする．A の固有値を λ1 ≤ λ2 ≤ · · · ≤ λn とすると，
Kahan のアルゴリズム [1] の丸め誤差特性を詳しく分析 A − λI の固有値は λi − λ(1 ≤ i ≤ n) であるから，し，精度保証アルゴリズムの構成する．Mathematica 上p(λ) : λi − λ > 0 となるλi の数.
でプログラムを作成し，数値実験によりその有効性を検m(λ) : λi − λ < 0 となるλi の数.
証する．

22 次形式と Sylvester の慣性則

2.12 次形式である．目標の固有値を λi とすると
{m(λ) < i =⇒ λ ≤ λim(λ) ≥ i =⇒ λ > λi

A = (aij ) を n 次実対称行列，x = (x1 x2 · · · xn ) を n次元列ベクトルとするとき
Tf = xT Ax =n∑

aij xi xji,j=1である．これを用いて，二分法のアルゴリズムにより，任意の精度で λi が包囲できる．すなわち，A の固有値を小さい順に λ1 < λ2 < · · · < λn とし，第 i 番の固有値 λi を求める Kahan のアルゴリズムは以下のようになる．

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

を変数 x1 , x2 , · · · , xn に関する 2 次形式という．

2.22. c = (ak + bk )/2 とし，{m(c) < i =⇒ [ak+1 , bk+1 ) = [c, bk ),変数変換

正則行列 P を用いて変数変換 x = P y を行うと 2 次形式 f は f = y T P T AP y と書き換えられる．これにともなって定義行列 A に起こる変換

m(c) ≥ i =⇒ [ak+1 , bk+1 ) = [ak , c).
である.

A → P T AP

3. bk+1 − ak+1 ≥ ϵ ならば，k → k + 1 として 2. へもどる．

を合同変換という．
定理 (Sylvester の慣性則) A を合同変換により対角行列

d1
D = diag(d1 , · · · , dn ) = 00..

.



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

4

修正 Cholesky 分解による符号数計算
修正 Cholesky 分解は n 次対称行列 A を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 = U ΛU Tとかける．これは対角行列 Λ から A への合同変換である．ゆえに Sylvester の慣性則より，正の固有値の数は p,負の固有値の数は m である．//

A = LDLT(1)と単位下三角行列 L と対角行列 D および LT の積に分解することである．[2] D の符号数を数えることにより A の符号数がわかる．A が対称三重対角行列

α1 β1

..
β.
 1 α2
A=

..
..

.
.
βn−1 
βn−1αnのとき，符号数の計算には L は不要ゆえに D = diag(di )のみを計算すると

d1 = α1 ,di = αi −2βi−1di−1(2 ≤ i ≤ n)となる．これが Kahan のアルゴリズムである．//(2)5後退誤差解析

7A − λI に対してアルゴリズム (2) の di を機械計算で求めた結果を d˜i とする．後退誤差解析により，d˜i は行列α̃1 β̃1..
 β̃.
 1 α̃2Ã − λI = 
 − λI..
..
.
.
β̃
n−1 
β̃n−1α̃nに対してアルゴリズム (2) を正確に計算した結果と等しく∥Ã − 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テスト行列として

2 −1
 −12

−1
A=
08である．これより A の固有値で小さい方から i 番目をλi (A) と書くと，|λi (Ã) − λi (A)| ≤ ∥∆A∥∞ ≤ δ(λ)(3)

0−12..
.

−1..
.
−1


 ∈ Rn×n..
.

2 −1 −12を用い，全固有値の包囲区間を計算した．
n = 128 の場合，全ての固有値の包囲に成功した．また，固有値の包囲区間幅の最大値 Max Width= 4.88 · · ·×10−15 で 2011 年の岡山亮士の卒業論文 [4] における結果
Max width= 3.28×10−14 よりも精密に精度保証ができた．
次に n = 2048 の場合，全ての固有値の包囲に成功した．また，固有値の包囲区間幅の最大値は n = 128 のときと同じで，4.88 · · · × 10−15 と大次元にも関わらず極めて精密に精度保証ができた．

である．また

u = 2−53 = 1.11 · · · × 10−16 ,ku|θk | ≤ γk ≡,1 − kuγ√3µ3 =1 + 1 − γ3Mathematica による数値実験

考察

結果として，修正 Cholesky 分解で 0 割り算になってしまった場合のコントロールがやや複雑になってしまったが，高精度であった岡山亮士の精度保証よりもさらに精密にできた．また今回は問題を n = 128, n = 2048 と二通りに設定して実験を行ったが，n の数を大きいものに置き換えても高精度な精度保証ができたことには驚きであった．これは，符号数が 2048 項の数列 dj を延々と漸化式計算した結果であることを考えると驚異的な精度である．

が成立する．

69精度保証付 Kahan の二分法

おわりに本研究では，Kahan の二分法を用いて対称三重対角行列の固有値を精度保証付きで計算することを目指した．後退誤差解析により，Strum 二分法より丸め誤差の影響が小さいことがわかった．実際に Mathematica のプログラ
ムを組み数値実験を行うと，修正 Cholesky 分解で 0 割り算が発生した場合のコントロールが難しくプログラムが複雑になってしまったが，すべての固有値を精度保証付
a − δ(a) ≤ λi (A)きで求めることに成功した．また，Sturm 二分法 [3] よりである．同じく，λ = b のときの後退誤差解析の結果は， 精密な精度保証が得られた．
式 (3) よりλi (A) < b + δ(b)後退誤差解析の結果 (3) を用いて，5 節のアルゴリズムの結果に精度保証を与えることが出来る．5 節のアルゴリズムで得られた小さい方から i 番目の固有値 λi の包囲区間を (a, b) とする．λ = a のときの後退誤差解析の結果は，式 (3) よりである．また，δ(λ) は |λ| の単調増加関数であることは明らかだから，|λ| = max{|a|, |b|}
として，厳格な不等式

a − δ(λ) ≤ λi (A) < b + δ(λ)

(4)

参考文献

[1] W. Kahan: Accurate eigenvalues of a symmetrictri-diagonal matrix, Technical Report of Stanford
University Stanford, CA, USA, 1966.
『数値計算の基礎と応用—数値解析学への入[2] 杉浦洋：門—』．サイエンス社，1997．

[3] 岡山亮士：『卒業論文』．南山大学，を得る．

2011．