対称行列に対する Householder 三重対角化の精度保証
2011SE217 小澤 大地
指導教員：杉浦 洋

1すなわち, i = 1, 2, · · · , n − 2 の順に適切にはじめに対称行列 A の固有値問題における，Householder 三重対角化法 [1] は，Householder 変換により A を相似変換して，三重対角行列 T を作り，問題を T の固有値問題に帰着させる．これを数値的に実行すると，T の近似 T̃ が求まる．本論の目的は T̃ の固有値と A の固有値の差を評価することである．

2(H(q)i =(3)を定め

区間演算

QAQT = T,
Q = H(q)n−2 , · · · , H(q)2 H(q)区間解析では，区間を数の拡張と考える．ここで，区間とは，閉区間[x，x] = {x ∈ R|x ≤ x ≤ x}である．区間を [x] = [x, x] と表すこともある．
関数 f : D ⊂ R → R を領域 D で連続な関数とする．関
数 f を区間関数λi (T̃ ) − λi (A) ≤ w (1 ≤ i ≤ n)

により IR 上の関数に拡張できる．Mathematica には基本的な実数関数の区間関数が実装されている．
区間を要素とする行列を区間行列と言う，Mathematica の区間演算で行列積 ABを含む区間[C] = [A][B] ∋ AB

基本的な定理精度保証に用いる基本的な定理を紹介する.

[定理 1] 対称行列 A ∈ Rn×n と行列 Q ∈ Rn×n について，が計算できる．

QQT = I + E, ∥E∥2 ≤ ϵ < 1Householder 三重対角化

Householder 行列 H(q) ∈ R√∥q∥2 = 2 により ，n×nn(1)で定義される対称直交行列である.
任意のベクトル a = (a1 , a2 , · · · an )T ∈ Rn に対し,Householder 行列 H(q) ∈ Rn×n が存在し,H(q)a = be1 , b = ±∥a∥2(6)のとき，B = QAQT とするとは，ベクトル q ∈ R ，H(q) = I − qq T

(5)を満たす w を計算することである.
ここで, λi (A) は行列 A の小さい方から i 番目の固有値である.

4[A] = ([aij ])(4)とする.
数値的に Householder 三重対角化を行うと, T の近似
行列 T̃ が求まる．Householder 三重対角化の精度保証とは,

f ([x]) = {f (x)|x ∈ [x]} ([x] ⊂ D)

3)
O∈ Rn×nH(qi )IiO|λi (A) − λi (B)| ≤(7)

[定理 2] ベクトル q ∈ Rn について，H(q) = I − qq T とすると∥H(q)H(q)T − I∥2 ≤ (∥q∥22 − 2)∥q∥22 .
[定理 3] 行列 H(q) ∈ Rn×n(2)とすることができる. 符号 ± は任意に選べるが, 計算上で桁落ちの少ない方を選ぶ.
Householder 三重対角化 [1] は, Householder 行列を用いて n 次対称行列 A を相似変換し, 対称三重対角行列 Tを作る.

ϵλi (B).
1−ϵ(8)(1 ∈ i ∈ k) に対して

Q = H(q)n−2 · · · H(q)2 H(q)1 ,(9)∥H(q)i H(q)Ti − I∥2 ≤ δi(10)(1 ≤ i ≤ k)

なら∥QQT − I∥ ≤k∏
i=1(1 + δi ) − 1 =: ϵ2(11)5Householder 三重対角化に対する今回の精度保証アルゴリズム

表 1 計算結果 (乱数行列)

与えられた対称行列 A ∈ Fn×n について，A(1) = A とする．(4.13) の Householder 変換 H1 (q1 ) を数値的に求めたものを()
10T
H1 (q̃1 ) =, q̃1 ∈ F(n−1)
0H(q̃1 )n48163264128とする．さらに，

w
3.47 × 10−14
2.01 × 10−13
1.09 × 10−12
7.26 × 10−12
5.18 × 10−11
4.35 × 10−10wold
5.43 × 10−14
6.98 × 10−12
4.33 × 10−8
11.5601.16 × 1015
2.53 × 1043

W
1.10 × 10−15
2.62 × 10−15
1.32 × 10−15
2.84 × 10−15
2.79 × 10−15
4.84 × 10−15(1)|q̃1T q̃1 − 2| ∥ q̃1 ∥2 ≤ ε1(1)と な る ． ε1を 区 間 演 算 で 求 め る ．こ こ で ，である．これを (4) に代入して
H1 (q̃1 )T A(1) H1 (q̃1 ) を区間演算したものをn−2∑ ε(1)(2)k|λi (T ) − λi (A)| ≤α + ε1(1) k(2)T(1)[A ] = H1 ([q̃1 ]) [A ]H1 ([q̃1 ])k=1 1 − εkを得る．これが今回のアルゴリズムである．
a11 [b1 ] [bT2 ]=  [b1 ]
6 数値実験
[bT2 ]
7 章のアルゴリズムを Mathematica で実現し，数値実とする．そして験を行った．区間演算は Mathematica の区間演算機能を用いた．n × n 対称行列 A はその成分を [−1, 1] 区間の一a11b 1 0T様乱数で生成した．
A(2) =  b1

今回のアルゴリズムで求めた三重対角行列を T ，誤差の0上界をωとする．とする．ここで，b1 = mid[b1 ], A2 = mid[A2 ] である．
このとき|λi (T ) − λi (A)| ≤ ω (1 ≤ i ≤ n)

[A2]A2(1)ε1(2)||A(1) ||2 + ε1|λi (A(2) ) − λi (A(1) )| ≤(1)
1 − ε1(12)を得る．A(2) は A(1) の第 1 列の第 3〜n 行が消去された対称行列であり，以上の操作は Householder 変換の代替となる．しかも，不等式 (3) が得られる．これを修正Householder 変換と呼ぶ．
Householder 三重対角化で，Householder 変換を，修正Householder 変換で置き換えると，行列の列

A(1) , A(2) , · · · , A(n−1)と 1 ≤ k ≤ n − 2 における不等式が保証される．青山 [1] のアルゴリズムで求めた誤差の上界を ωold とする．また

W = max |λi (T ) − λi (A)|1≤i≤nを Mathematica の多倍長機能を用い，10 進 100 桁計算したものを用い，仮に真値とした．
n = 2m (2 ≤ m ≤ 8) で行った実験の結果を表 1 に示す．
青山のアルゴリズムは n ≥ 30 では有効な左上を作れない．また n ≤ 20 でも一貫して今回のアルゴリズムより劣る．
今回のアルゴリズムは，n = 128 でも有効な上昇を計算することに成功している．

(1)|λi (A(k+1) ) − λi (A(k) )| ≤εk(2)(1)

1 − εk||A(k) ||2 + εkを得る．T = A(n−1) は三重対角行列であり，ゆえに，α1 =||A(1) ||∞ とし

αk =

1
(1)
1 − εk−1(2)

αk−1 + εk−1

(k = 2, 3, · · · , n − 2)(13)とすれば||A(k) ||2 ≤ αk(k = 1, 2, · · · , n − 2)7

おわりにHouseholder 三重対角化の精度保証の新しいアルゴリズムを構成した．今回のアルゴリズムは昨年の青山のアルゴリズムより精度が高く，大きな n でも有効であった．

参考文献[1] 青山 亨：対称行列に対する Householder 三重対角化の精度保証，南山大学 (2014)．