縮小対角優位行列の高精度固有分解の計算
M2007MM024 太田裕文
指導教員：杉浦洋

1はじめに3縮小対角優位行列は，対角優位行列の左右から対角行列を掛けたものである．本論文は，Barlow と Demmel[1]に従って，縮小対角優位行列の固有分解と特異値分解に関する摂動理論を扱ったものである．
一般の行列の固有値，特異値に関する摂動理論 [2,3] によれば，それらの数値計算誤差は機械精度と行列のノルムに積に比例する誤差限界をもつ．それゆえ，絶対値の小さな固有値と特異値の相対精度が劣化する．
それと対照的に，縮小対角優位行列においては，その固有値と特異値は絶対値の大小にかかわらず，高い相対精度で計算しうる．このことを，摂動理論により示すことが今回の目標である．
実際に，このような高精度で固有値，特異値を計算するためには，特別なアルゴリズムを必要とする．本論文
では，二重対角行列の特異値，正定値対称三重対角行列の固有値に対する高精度計算プログラムを作成した．また，数値実験によりその有効性を確認した．

2基本概念の定義

本論文では，∥ · ∥ は 2− ノルムを表す．行列 A の ij 要
素を Aij で表す．
定義 2.1 A を m × n 行列とする．||| · ||| を行列ノルムとする，A の対角部を D，非対角部を N とするとき(A = D + N )．|||N ||| < mini |Dii | が成立するなら，Aは ||| · ||| に関して対角優位であるという．//
定義 2.2 上の定義で，ある 0 ≤ γ < 1 に対して |||N ||| ≤γ mini |Dii | が成立するとき，A は ||| · ||| に関してγ − 対角優位 (γ − d.d.) であるという．
//
定義 2.3 行 列 H ≡ ∆1 A∆2 で ∆1 ，∆2 が 対 角 行列，A が対角要素 ±1 の γ − d.d. のとき，H を γ −s.d.d.(scaled diagonally dominant，縮小対角優位) と呼ぶ．//
定義 2.4 対称行列 H ≡ ∆A∆ で ∆ が対角，A が対角要
素 ±1 の γ − d.d. のとき，H を対称 γ − s.d.d. と呼ぶ．//
定義 2.5 ペンシル H − λM, H, M ∈ Rn×n が s.d.d.定値であるとは，H と M が対称 s.d.d.であり，M が正定値であることである．
//
定義 2.6 行列 H が単位固有値 λi をもつ対称行列とするとき，gap(λi ) = min |λi − λk | > 0．
//k̸=i(1)

定義 2.7 n 次正方行列 H の固有値を λ1 · · · λn とする
とき，

relgap(λi ) = min|λi − λj |j̸=i |λi λj |1/2．
//(2)

二重対角行列の特異値のための摂動定理

二重対角行列の要素に相対摂動が与えられたときの特異値の相対摂動に関して，次の定理が成立する．
定理 3.1 B は二重対角行列a1 b1

.. ..
.
.

B=.
..





bn−1 an(3)とする．そうでなければ，B は独立な部分行列に分解するので，ai と bi はノンゼロと仮定する．B + δB は，ai の代わりに αi ai そして bi の代わりに βi bi を要素とする摂動二重対角行列とする．そのとき，Bi の特異値 δ1 ≤ · · · ≤ δnと δ̂1 ≤ · · · ≤ δ̂n は gl と gu を以下のように定義し，

gl ≤δ̂i≤ guδi(4)を満たす．また，B の特異値 σ1 ≤ · · · ≤ σn と，B + δBの特異値 σ1′ ≤ · · · ≤ σn′ とする．このとき，{}
βj · · · βkS≡
:1≤j ≤k ≤n−1
αj+1 · · · αk}{αj · · · αk:1≤j≤k≤n(5)∪βj · · · βk−1について，gu = max S ，gl = min S(6)とすると，

gl ≤σi′≤ guσi(7)である．//
行列要素 ai ，bi に対する相対摂動が小さいとき，すな
わち αi ，βi ∼= 1 のとき，gu ，gl ∼= 1 となる．すなわち，′∼σi /σi = 1 が σi の大小にかかわらず成立する．これは，特異値に対する相対摂動が，特異値自身の大小にかかわらず小さいことを意味する．

4固有値の摂動定理

H = ∆A∆ = ∆(E + N )∆ を γ − s.d.d. 対称行列とする．ここで，∆ は対角行列，A の対角要素は ±1，Eは A の対角部，N は A の非対角部で，∥N ∥ = γ < 1である．H の各要素に，ε > 0 以下の相対摂動を与えたものを H + δH とする．|δHij /Hij | ≤ ε である．H +δH = ∆(A + δA)∆ とすると，δA = ∆−1 δH∆−1 について |δAij /Aij | = |δHij /Hij | ≤ ε．条件より，|Aij | ≤ 1 ゆえ，|δAij | ≤ ε を得る．従って，H の要素に対する ε 以下の相対摂動は，A の要素に対する ε 以下の絶対摂動と見なせる．
以下では，H に対する相対摂動を A = ∆−1 H∆−1 への絶対摂動と見なし，解析する．
正定値対称行列の標準固有値問題について，次の定理が成り立つ．
定理 4.1 V を正則とする．A = I +N において，∥N ∥ =γ < 1 とする．このとき，H = V T AV は正定値対称である．δH を H に対する対称摂動とし，∥V −T δHV −1 ∥ ≡η < 1 − γ とする．このとき，H の固有値 λ1 ≤ · · · ≤ λnと H + δH の固有値 λ′1 ≤ · · · ≤ λ′n は1−λ′iηη≤≤1+1−γλi1−γ(8)を満たす．//
この定理で，V を対角行列，N を A の非対角部とすれば，H は γ − s.d.d. となり，ただちに次の系を得る．
系 4.1 H が対称正定値 γ − s.d.d. であり，δH を H に対する対称摂動で ∥∆−1 δH∆−1 ∥ ≡ η < 1 − γ を満たすとする．このとき，

1−ηλ′η
≤ i ≤1+1−γλi1−γ(9)を満たす．//
行列 H ，M が共に正定値行列のとき，ペンシル H −λMの固有値問題について次の定理が成り立つ．

命題 4.1 H = ∆A∆ を 2− ノルムに関して対称 γ−s.d.d.
行列とする．ここで，∆ は対角，A は対角要素が ±1 の対称γ − s.d.d. である．δH を ∥∆−1 δH∆−1 ∥ ≡ η < (1 − γ)/nを満たす H の摂動とする．そのとき，H の固有値を λ1 ≤· · · ≤ λn ，H + δH の固有値を λ′1 ≤ · · · ≤ λ′n とすると，1−λ′nη≤ i ≤1−γλi()−1nη
1−1−γ(12)が成立する．
右辺の因子 n は，次の通りにわずかに命題の条件を強くすることで取り去ることができる．
定理 4.3 H = ∆A∆ を 2-ノルムに関して対称 γ − s.d.d.
行列とする．ここで，∆ は対角，A は対角要素が ±1 の対称 γ − s.d.d. である．δH は ∥∆−1 δH∆−1 ∥ ≡ η 対称摂動，H + ξδH はすべての 0 ≤ ξ ≤ 1 に対して γ − s.d.d.と仮定する．そのとき，λ1 ≤ · · · ≤ λn を H の固有値，λ′1 ≤ · · · ≤ λ′i を H + δH の固有値とすると，(())−ηλ′iηexp≤≤ exp(13)1−γλi1−γが成立する．//
以上の定理，命題，系は，行列 H の要素に対する相対摂動と，固有値の相対摂動の関係を表す．
定値ペンシル H − λM で H が正定値でないとき，次のアルゴリズムにより，標準固有値問題に変換する．
アルゴリズム 4.1 γ − s.d.d. 定値ペンシル H − λM の固有値問題を s.d.d. 行列 Y の標準固有値問題に帰着する．

1/2定理 4.2 行列 NH ，NM は ∥NH ∥ ≤ γ < 1，∥NM ∥ ≤(1)D1 = diag(Mii ) とし，H1 = D1−1 HD1−1 と M1 =γ < 1 を満たし，行列 VH ，VM は正則とする．AH =−1−1D1 M D1 を計算する．このとき，M1 は単位対角を持I + NH ，AM = I + NM に対し，H = VHT AH VH ，ち，そして普通の意味で，対角優位である．
TM = VM
AM VM は正定値行列である．δH と δM は，(2)P は H2 = P H1 P T において，対角要素が (H2 )11 ≤H と M の対称摂動で ∥VH−T δHVH−1 ∥ ≤ η < 1 − γ ，· · · ≤ (H2 )nn と並ぶように選ばれた置換行列とする．
−T−1
∥VM
δM VM
∥ ≤ η < 1 − γ を満たす．このとき，ペM2 = P M1 P T とする．
ンシル H − λM の固有値を λ1 ≤ · · · ≤ λn ，(H + δH) −(3)L は M2 の 下 三 角 コ レ ス キ ー 因 子 と す る ．Y =λ(M + δM ) の固有値を λ′1 ≤ · · · ≤ λ′n とすると，
L−1 H2 L−T とする．このとき，Y と H − λM は同じ固有値を持つ．
1−γ−ηλ′1−γ+η≤ i ≤(10) 命題 4.2 H − λM は γ − s.d.d. 定値ペンシルとする．そ1−γ+ηλi
1−γ−ηして，Y は上の帰着アルゴリズムの出力とする．
を満たす．//
1+γ· [2 + ((2n)1/2 + 1)γ 1/2 ] · γ 1/2 (14)γ ′ ≡ ((2n)1/2 + 1) ·1−γ
この定理で，VH ，VM を対角行列，NH ，NM を A の非対角部とすれば，H も M も γ − s.d.d. となり，ただちと定義する．そのとき，もし，γ が十分小さく，に次の系を得る．
(n + 1)λ′ + γγ̄ ≡≤1(15)
系 4.2 H と M が対称正定値 γ − s.d.d. であり，H と1 − γ′−1−1M の摂動 δH と δM が ∥∆H δH∆H ∥ ≡ η < 1 − γ 及び−1が真であるならば，Y は γ̄ − s.d.d. になる．
∥∆−1M δM ∆M ∥ ≡ η < 1 − γ を満たすとする．このとき，

λ′1−γ+η1−γ−η≤ i ≤1−γ+ηλi1−γ−η(11)を満たす．//
一般の対称行列の標準固有値問題について，定理 4.1 より結果の弱い次の命題と定理が成り立つ．

これにより，定値ペンシルの固有値の摂動問題は，標準固有値問題に帰着され，命題 4.1，定理 4.3 を用いることができる．

5

固有ベクトルの摂動定理

この節では，節 4 と同じ小さい相対摂動の下で対称 s.d.d.の固有ベクトルの感度を議論する．固有ベクトルに対する標準的な摂動理論では，固有値 λi に対応する固有ベクトルの感度を決定するのは，λi と隣接固有値との距離gap(λi ) である．すなわち，y を A + δA の単位固有ベクトルとし，α = y T Ay をレイリー商とする．また，λi をα に近い A の固有値とし，zi を λi の単位固有ベクトルとする．そして，θ(zi , yi ) を y と zi 間の鋭角の角度とす
ると，gap を用いて，標準的な摂動理論の誤差限界は以下のように表される．

sin θ(y, zi ) ≤ 4∥δA∥/gap(λi )[5, p248]． (16)

命題 5.2 (3) の B について
−1/2

B B T DL− I∥ ≤ γL}{aj+1 bjbn−1, max
， (21)= 2 maxrn−1 1≤j≤n−2 rj rj+1−1/2

B T BDR∥DL∥DR−1/2−1/2− I∥ ≤ γR}
aj bj
b1, max． (22)
= 2 maxs1 2≤j≤n−1 sj−1 sj√√
ここで，rj = a2j + b2j ，sj = a2j+1 + b2j である．γL{

1/2
それに対し，対称 s.d.d. 行列の固有ベクトルの感度を決 と γR は両方とも 2 以下である．aj > 3 bj ならば，そ1/2定するのは，隣接固有値との相対距離 relgap(λi ) である． のとき γL < 1，そして aj > 3 bj−1 ならば，そのときγR < 1 である．
定理 5.1 H = ∆A∆ は 2-ノルムに関する γ −s.d.d. 対称ある条件の下で，BB T ，B T B は γ − s.d.d. になる．こ行列とする．行列 E は対称で，∥E∥ = 1 とする．H(ζ) =∆(A+ζE)∆ と定義する．λi (ζ) は H(ζ) の固有値で，λi (0) のとき，B の要素の小さい相対摂動は，特異値の相対 gapは単純とする．したがって十分小さい ζ で単位固有ベク と関係した小さな相対摂動を B の特異ベクトルにもた
トル xi (ζ) は ζ の連続関数になる．そのとき，十分小さ らす．
このような条件で，BB T ，B T B は γ − s.d.d. となり，い ζ に対して
B の特異ベクトルの摂動について，定理 5.1，定理 5.2 を(n − 1)ζ
∥xi (ζ) − xi (0)∥ ≤+ O(ζ 2 ) (17) 用いることができる．
(1 − γ)relgap(λi )6

固有ベクトルの条件数である．
式 (16)，式 (17) の右辺で問題に対する摂動のノル
命題 5.1 H = ∆A∆ は，2-ノルムに関して γ − s.d.d.
対称行列とする．H の固有値を λ1 ，λ2 ，· · · ，λn ，対応する ム，∥δA∥，ζ にかかる係数 4/gap(λi )，(n − 1)/((1 −固有ベクトルによる正規直行系を x1 ，x2 ，· · · ，xn とする． γ)relgap(λi )) を条件数という．
Λ = diag(λ1 ，λ2 ，· · · ，λn )，X = [x1 ，x2 ，· · · ，xn ] とすると，H = XΛX T である． xi (j) は xi の第 j 成分とする． 命題 6.1 H が単純な固有値 λi をもつ対称行列とする．
すなわち，gap(λi ) = mink̸=i |λi − λk | > 0 である．そのそのとき()とき，λi に”対応”する H + δH の固有値が重複であるよ)(1/21/2λiλj1+γ¯· min，．
(18) うな最小の ∥δH∥ は，|xi (j)| ≤ xi (j) ≡1−γλjλimin ∥δH∥ =また，(1+γ1−γ)(∆ii ∆jj，∆jj ∆ii)

gap(λi )2(23)である．”λi に対応する固有値は重複である”というのは，H + ξδH(0 ≤ ξ ≤ 1) の ξ に関して連続な固有値λ (ξ)(λi (0) = λi ) で，λi (1) が重複になる，ということで式 (18)，(19) の右辺の因子 (1 + γ)/(1 − γ) は，Barlow， iある．
3/2
Demmel[1] の {(1 + γ)/(1 − γ)}から改善されたもの
ここで，min は H + δH の固有値 λi (1) が重複であるである．
ものについての最小値である．
定理 5.2 H(ζ) は定理 5.1 と同様とする．xi (ζ)(j) はH(ζ) の第 i 単位固有ベクトルの第 j 成分を示すとする． 今までこの論文で使用されているスケーリングされた
x̄i (j) を命題 5.1 で定義した，xi (j) = xi (0)(j) における 方法で距離を測定し，条件数として 1/relgap(λi ) を用い摂動の上限とする．そのとき，るなら，同様の関係が γ − s.d.d. 対称の行列に成り立つ．
|xi (ζ)(j) − xi (0)(j)|
これは，1/relgap(λi ) が条件数にふさわしい幾何学的な21/2 (n − 1)2}{· x̄i (j) + O(ζ )．(20) 性質をもつ，ということを示した点で興味深い．
≤ζ·(1 − γ) · min relgap(λi )，2−1/2
命題 6.2 H = ∆A∆ は，単純な固有値 λi をもつ γ −最後に，二重対角行列の特異値のための摂動論を考える． s.d.d. 対象行列とする．relgap(λi ) = mink̸=i |λi − λk | ·−1/2> 0 である．さらに，relgap(λi ) ≤ 2−1/2 とB の左特異ベクトルは BB T の固有ベクトルである．ま |λi λk |仮定する．λi とそれの最も近い隣接するものは，因子 5た，B の右特異ベクトルは，B T B の固有ベクトルであと 2 の間で異なることを意味している．そのとき，”対応る．以下の命題により，一定の条件の下で二重対角行列 する”∆(A + δA)∆ の固有値が重複であるような最小のの特異ベクトルの摂動問題は γ − s.d.d. 三重対角行列の ∥δA∥ は，√固有ベクトルの摂動問題に帰着される．
2n(1 + γ)2 relgap(λi )(1 − γ)relgap(λi )|xi (j)| ≤· min(19)√≤ min ∥δA∥ ≤2n + relgap(λi ))(1 − γ)3(24)を満たす．relgap(λi ) << 1 に対して，min ∥δA∥ 上の下限は，relgap(λi )(1−γ)/(21/2 ·3·n)+O((relgap(λi ))2 ) と等しい．言い換えると，上限と下限両方は O(relgap(λi ))である．

8

数値実験

C 言語倍精度で正定値三重対角行列の固有値を求めるプ式 (24) の右辺の因子 (1 + γ)2 /(1 − γ)3 は，Barlow， ログラムを作成した．丸め誤差単位は u = 2−53 = 1.1 ×
Demmel[1] の (1 + γ)4 /(1 − γ)3 から改善されたもので 10−16 である．
ある．
対角要素，副対角要素がそれぞれ Ai,i = 3r2i−2 ，Ai,i+1 = r2i−1 (r = 10−6 ) の 10 次三重対角行列 A の7 アルゴリズム固有値を計算した．また，Mathematica により，計算桁7.1 二重対角行列の特異値
数 200 桁で計算された固有値を真値とした．
二重対角行列の特異値は，0-シフト Givens QR 法に改λ相対誤差良された収束判定条件を用いることにより，高い相対精度で求めることが出来る．
3.00000
2.37 × 10−16−12Givens QR 法は，二重対角行列の特異値を求める標準
2.66667 × 10
2.51 × 10−16−24的なアルゴリズムの一つである．通常，収束を加速する
2.62500 × 10
4.36 × 10−16ため，シフト付き Givens QR 法が用いられる．これは，2.61905 × 10−36
3.73 × 10−16−48
単位行列のシフトパラメータ倍を減じて QR 反復を行う
2.61818 × 10
6.98 × 10−16−60方法である．Givens 回転行列は桁落ちなく高精度で計算
2.61806 × 10
5.33 × 10−16できる．しかし，それを用いた行列の回転は加減算を含
2.61804 × 10−72
4.12 × 10−16−84
み，ここで桁落ちによる精度の劣化が起こる可能性があ
2.61803 × 10
5.37 × 10−16る．これに対し，シフトパラメータを 0 とする 0− シフ
2.61803 × 10−96
5.68 × 10−16−108
ト Givens QR 法を二重対角行列に適用した場合，Givens
2.61803 × 10
7.30 × 10−16
回転は加減算を含まず乗算だけで構成できるので，相対すべての固有値が最大相対誤差 7.3×10−16 で計算するこ精度を損なわない．とが出来た．また，絶対値最小固有値 (絶対値 2.6×10−108 )n 次上三角二重対角行列を B とする．Givens QR 反復で も理論通り高精度で計算することができた．これは従来生成される上三角二重対角行列の列を B = B (0) , B (1) , · · · の誤差理論ではあり得ないことである．
{}とする． B (k) は特異値を対角要素にもつ対角行列に収(k)束する．副対角要素 Bi,j+1 は 0 に収束するが，数値計算では絶対値が十分小さいと判定されたら，収束したと見なし 0 で置き換える．
(k)
通常の判定条件は Bi,i+1 ≤ ∥B∥ u であるが，この条(k)件を満たす Bi,j+1 を 0 で置き換えると，小さな特異値に対して大きな相対摂動を生じる．我々は，新しい摂動理論に基づき判定条件
(k)

(k)

Bi,i+1 ≤ Bi,i u(k)

9

おわりに

以上で述べた摂動定理は，s.d.d. 対称行列の固有値問題が，高い相対精度で解ける可能性があることを示している．実際に，数値実験を行ってみて，高い相対精度を達成するということを実証することができた．
今回のアルゴリズムは二重対角行列の特異値，正定値三重対角行列の固有値の計算に限定されたものである．一般の行列の特異値問題，固有値問題に対する数値アルゴリズムの開発は今後の課題である．

参考文献

を提案する [4]．これを満たす Bi,i+1 を 0 で置き換えた際の特異値の相対摂動が丸め誤差単位 u 程度であることが， [1] Jesse Barlow and James Demmel : Computing accurate eigensystems of scaled diagonally dominant
理論的に保証される．
matrices，SIAM J．Numer．Anal．Vol．27．No．3．
7.2 正値対称三重対角行列の固有値
pp.762-791，June 1990．
n 次正値対称三重対角行列を A のすべての固有値 [2] G．Golub and C．Vanlon：Matrix Computations，
λ1 , λ2 , · · · , λn を高い相対精度で求めるアルゴリズムを
The Johns Hopkins University Press，Baltimore，示す．
MD，1983．
まず，A = B T B と Cholesky 分解し，Cholesky 因子 B
[3] B．Parlett：The Symmetric Eigenproblem，Prenticeの特異値 σ1 , σ2 , · · · , σn を前節の 0− シフト Givens QR 法Hall，Englewood Cliﬀs，NJ，1962．
で求める．これにより，A の固有値を λi = σi2 (1 ≤ i ≤ n)[4] J．Demmel and W．Kahan : Accurate singular valで計算する．
ues of bidiagonal matrices，SIAM J．Sci．Statist．
数値 Cholesky 分解により B の要素に発生する相対摂．
，11 (1990)．
Comput動は 7u + O(u2 ) である．ここで，O(u2 ) は次数 n の多項式である．ゆえに，定理 3.1 より，B の特異値 σi の特異値は高い相対精度で計算され，A の固有値 λi も高い相対精度で計算できる．