数理解析研究所講究録
第 1719 巻 2010 年 7-117

上二重対角行列の最小特異値の下界に関する最近の進展について
山下

京都大学大学院情報学研究科

木村巧，欣司 ，中村佳正(Takumi Yamashita, Kinji Kimura, Yoshimasa Nakamura)

Graduate School of Informatics,
Kyoto University

1

はじめに行列の特異値について，最小特異値のよりシャープな下界を少ない計算量で求めることは理論，応用の両面で有用である．特に上二重対角行列の最小特異値の下界については，その下界の 2 乗をシフト付き特異値計算アルゴリズム [1,2] におけるシフト量として用いることができる．

本講演では，最小特異値の下界に関する著者達による最近の研究について報告する．

2 上二重対角行列の最小特異値の一般化 Newton 下界
正則な $N\cross N(N\geq 2)$ 上二重対角行列 $B$ について，その上副対角成分は全て零でないと仮定する．

このとき， の特異値を$B$

列 $(B^{T}B)^{M}$ の固有値を$tJ_{1},$$\cdots,$$\lambda_{1}^{(M)},$$o_{N}(0_{1}>\cdots>0_{N}>0)$劫

$\cdots,$$\lambda_{N}^{(M)}(\lambda_{1}^{(}$とする．任意の正の整数 $M$ に対し，正定値行$>\cdots>\lambda_{N}^{(M)}>0)$とする．ここに，$\lambda_{i}^{(M)}=0_{i}^{2M}(1\leq i\leq N)$$(B^{T}B)^{M}$ の固有方程式
である．

$f(\lambda^{(M)})\equiv\det((B^{T}B)^{M}-\lambda^{(M)}I)=0$を考える．

$(B^{T}B)^{M}$が正定値行列であることから，(切 $=0$ を起点としてよく知られた Newton 法を$\lambda$一回適用して得られる値$- \frac{f(0)}{f(0)}=(\sum_{i=1}^{N}\frac{1}{\lambda_{i}^{(M)}})^{-1}=(\sum_{i=1}^{N}\frac{1}{o_{i}^{2M}})^{-1}\equiv \mathcal{F}_{M}^{1}$

は $(B^{T}B)^{M}$ の最小固有値$\lambda$

卿の下界を与える．ここに， は，$J_{M}$$J_{M}=Tr(((B^{T}B)^{M})^{-1})=Tr(((BB^{T})^{M})^{-1})$を満たす．これより， の最小特異値$B$$0_{N}$の下界として次の次数 $M(M=1,2, \cdots)$ 一般化 Newton

下界$\Theta_{M}=J_{M}^{-1}m$$M=1$ の場合については，例えば， の主対角成分および副対角成分が全て正である場が得られる．
合に Femando と Parlett [1] か 主対角成分および副対角成分が全て非零である場合に von Matt [3]$M$ に対してがの最小特異値の下界として同様の結果を与えている．一般化 Newton 下界
は，単調増加性を示す．すなわち，$B$$\searrow$$B$$\Theta_{M}$$\Theta_{1}<\Theta_{2}<\cdots<0_{N}$$Marrow\infty$ の極限でである．さらに，$\lim_{Marrow\infty}\Theta_{M}=(J_{N}$である．

$\Theta_{M}$は最小特異値に収束する．すなわち，

83 一般化 Newton 下界を計算するための公式 (減算を含む形式)$(B^{T}B)^{M}$ または $(BB^{T})^{M}$ の逆行列のトレースが分かればよ
一般化 Newton 下界を計算するには，い．そのためには，これらの逆行列の対角成分を求めればよい．次の記号を導入する．上二重対角行
列 の第 行目の対角成分，副対角成分をそれぞれ $b_{j}(1\leq i\leq N)$ , ci $(1 \leq i\leq N-1)$ とする．次のよ$B$$i$うに行列$V^{(m)},$$Y^{(q)}(0\leq m\leq M, 0\leq q\leq M-1)$$W^{m)},$ $X^{(q)},$を定める．

$\{\begin{array}{l}V^{(m)}\equiv((B^{T}B)^{m})^{-1},W^{m)}\equiv((BB^{T})^{m})^{-1},X^{(q)}\equiv(B(B^{T}B)^{q})^{-1}=((BB^{T})^{q}B)^{-1},Y^{(q)}\equiv(X^{(q)})^{T}.\end{array}$$\mathcal{V}^{\langle m)},$$W^{(m)},$ $X^{(q)},$の第 行目の対角成分をそれぞれ$Y^{(q)}$$i$$z_{i}^{(q)}\equiv b_{i}(x_{i}^{(q)}+y_{i}^{(q)})(1\leq i\leq N)$$v_{i}^{(0)}=1$$v_{i}^{(m)},$$w_{i}^{(m)},$を定義する．

$V^{(0)}=W^{(0)}=\Gamma^{1}=I$,$w_{i}^{(0)}=1$である．このとき，求めるべき対角成分,$v_{i}^{(M)},$$x_{i}^{(q)},y_{i}^{(q)}(1\leq i\leq N)$とする．さらに，より，(1)$(1\leq i\leq N)$$w_{i}^{(}$

閉について，次の定理が成り立つ．

定理 1[4]$M$を任意の正の整数とする．求めるべき逆行列$(1\leq i\leq$$((B^{T}B)^{M})^{-1},$ $((BB^{T})^{M})^{-1}$の対角成分$v_{i}^{(M)},$$w_{i}^{(M)}$

めは，式 (1) および次の漸化式$v_{i}^{(J)}= \frac{1}{b_{i}^{2}}(c_{i}^{2}v_{i+1}^{(p)}+z_{i}^{(p-1)}-w_{i}^{(p-1)})$,,$w_{j}^{(p)}= \frac{1}{b_{j}^{2}}(c_{j-1}^{2}w_{j-1}^{(p)}+z_{j}^{(p-1)}-v_{j}^{(p-1)})$$z_{j}^{(q)}=z_{j-1}^{(q)}+2(v_{j}^{(q)}-w_{j-1}^{(q)})$

$(p)$

1

$(p-1)$’$v_{N}=w_{N}\overline{b_{N}^{2}}$

$(p)$

1

,

$(p-1)$$w_{1}$$=v_{1}\overline{b_{1}^{2}}$’$z_{1}^{(q)}=2v_{1}^{(q)}$から有限回の四則演算によって得られる．ここに，$i,$$1\leq p\leq M,$ $0\leq q\leq M-1$

である．

証明は [4] を参照されたい．

$\square$$j,$ $p,$$q$は，整数 1 $\leq i\leq N-1,2\leq j\leq N$ ,

9
注意 2のとき，定理 1 の$p=1$$v_{i}^{(p)},$$w_{i}^{(p)}(1\leq i\leq N)$に関する漸化式は次のように簡略化される．
$(1\leq i\leq N-1)$ ,$v_{i}^{(p)}= \frac{1}{b_{i}^{2}}(c_{i}^{2}v_{i+1}^{(p)}+1)$$(2\leq j\leq N)$$w_{j}^{(p)}= \frac{1}{b_{j}^{2}}(c_{j-1}^{2}w_{j-1}^{(p)}+1)$,

1$(p)$$v_{N}=\overline{b_{N}^{2}}$

1

$(p)$
$w_{1}$’.

$=\square \overline{b_{1}^{2}}$を計算するためのコストは無視できるとすると，定理 1 および注意 2 に基づいて一般化 Newton を計算するために必要な演算回数は $O(MN)$ のオーダーである． は固定された定数であるから，このオーダーは $O(N)$ と見倣せる．ただし，数値計算を行う場合，$M\geq 2$ のときは計算中の変数に大きな相対誤差が入り得る．例を以下に示す． をマシンイプシロンとする．数値計算において，乗算と除算は，加算と減算より先に行われると仮定する．さらに，計算式に括弧がある場合，括$J_{M}$から$J_{M}^{m^{1}}$$M$$\epsilon$

弧内の計算が先に行われると仮定する．

$b_{1}^{2}v_{2}^{(1)}<\epsilon$

おいて， から$b_{1}$$b_{1}^{-2}$

かつである場合を考える．数値計算に$c_{1}^{2}v_{2}^{(1)}<\epsilon$を計算したとき，この値が丸め誤差を含めて$\beta$になるものとする．数値計算により小さおいて，注意 2 の漸化式によって計算されるについて，いとすると，の計算においていわゆる情報落ちによってこれが無視され，最終的に は$c_{1}^{2}v_{2}^{(1)}$$v_{1}^{(1)}$が丸め誤差を含めても$\epsilon$$v_{1}^{(1)}$$c_{1}^{2}v_{2}^{(1)}+1$$\beta$$2=2^{1}$ より，2 をかける乗算では誤差が混入しないと仮定する．このとき，数値計算におけになる．

るは $z_{1}^{(1)}=2v_{1}^{(1)}=2\beta$$z_{1}^{(1)}$となる．

$w_{1}^{(1)}=\beta$

であることを考慮すると，数値計算として，$z_{2}^{(1)}=z_{1}^{(1)}+2(v_{2}^{(1)}-w_{1}^{(1)})=2\beta+2\cross(v_{2}^{(1)}-\beta)$

であるが，括弧内の計算において$v_{2}^{(1)}$は無視されるので，この計算結果は になる．厳密な計算にお$0$

いては.

$z_{2}^{(1)}= \frac{2}{b_{1}^{2}}(c_{1}^{2}v_{2}^{(1)}+1)+2(v_{2}^{(1)}-\frac{1}{b_{1}^{2}})=2(\frac{c_{1}^{2}}{b_{1}^{2}}+1)v_{2}^{(1)}>0$

であるから，この場合の数値計算における$z_{2}^{(1)}$の相対誤差は 1 のオーダーである．この問題を考慮し，定理 1 の漸化式を改良する．

4 一般化 Newton 下界を計算するための公式 (減算を含まない形式)
定理 1 の漸化式から出発してこれと等価で減算を含まない漸化式が導出される．

まず，以下の値$\check{B}_{i}=\frac{1}{b_{i}^{2}}$$\check{B}_{i}$,$(1 \leq i\leq N)$ $F_{i}(1\leq i\leq N-1),\tilde{F}_{i}(2\leq.i\leq N)$$(1\leq i\leq$

ハリ，$F_{i}=c_{i}^{2}\check{B}_{i}$$(1 \leq i\leq N-1)$ ,$\tilde{F}_{i}=c_{i-1}^{2}\check{B}_{i}$$(2\leq l\leq N)$.

を導入する．

10$p=1$ のとき，注意 2 の漸化式は次のように書き直せる．
この記号を用いると，$\{\begin{array}{ll}v_{i}^{(p)}=F_{j}v_{i+1}^{(p)}+B_{i} (1 \leq i\leq N-1),w_{j}^{(p)}=\tilde{F}_{i}w_{j-1}^{(p)}+\check{B}_{i} (2\leq j\leq N),v_{N}^{(p)}=\check{B}_{N}, w_{1}^{(p)}=\check{B}_{1}.

さらに，次 o) 値

(2)\end{array}$$g_{i}^{(r)},\tilde{g}_{i}^{(r)}(1\leq i\leq N, r=1,2, \cdots)$を定義する．

定義 3.
.
.

$g_{i}^{(r)}(1\leq i\leq N, r=1,2, \cdots)$

任意の正の整数$r=1$$r$は次で定義される．

に対し，$g_{N}^{(r)}=0$

.

に対し，

$g_{i}^{(1)}=F_{i}v_{i+1}^{(1)}(1\leq i\leq N-1)$$r=2,3,$$\cdots$.

に対し，$g_{i}^{(r)}=F_{i}g_{i+1}^{(r)}+B_{i+1}g_{i}^{(r-1)}+ \sum_{k=1}^{r-1}g_{i+1}^{(k)}g_{i}^{(r-k)}(1\leq i\leq N-1)$

.

定義 4

.

$\tilde{g}_{i}^{(r)}(1\leq i\leq N, r=1,2, \cdots)$.
.

任意の正の整数$r=1$$r$は次で定義される．

に対し，$\tilde{g}_{1}^{(r)}=0$.
.

に対し，$r=2,3,$$\tilde{g}_{i}^{(1)}=\tilde{F}_{j}w_{i-1}^{(1)}(2\leq i\leq N)$$\cdots$.

に対し，$\tilde{g}^{(r)}=\tilde{F}_{i}\tilde{g}_{i-1}^{(r)}+B_{i-1\tilde{g}_{i}^{(r-1)}+\sum_{k=1}^{r-1}\tilde{g}_{i-1}^{(k)}\tilde{g}_{i}^{(r-k)}}(2\leq i\leq N)$

注意 5$g_{N-1}^{(r)},\tilde{g}_{2}^{(r)}(r=2,3, \cdots)$

はそれぞれ，$g_{N-1}^{(r)}=\check{B}_{Ng_{N-1}^{(r-1)},\tilde{g}_{2}^{(r)}=B_{1\tilde{g}_{2}^{(r-1)}}}$

で計算できる．

これらの定義式を用いて，次の定理が証明される．
定理 6$M\geq 2$のとき，求めるべき逆行列$((B^{T}B)^{M})^{-1},$ $((BB^{T})^{M})^{-1}$

式 (2) および次の漸化式$v_{N}^{(s)}=B_{Nw_{N}^{(s-1)}}$$w_{1}^{(s)}=\check{B}_{1}v_{1}^{(s-1)}$

,,
$s-1$

$v_{i}^{(s)}=F_{j}v_{i+1}^{(s)}+ \check{B}_{i^{W_{i}^{(s-1)}}}+2\sum g_{i}^{(k)}w_{i}^{(s-k)}$

,

$k=1$

$s-1$
$w_{j}^{(s)}= \tilde{F}_{j}w_{j-1}^{(s)}+B_{j}v_{j}^{(s-1)}+2\sum g_{j}^{(k)}v_{j}^{(s-k)}\sim$$k=1$の対角成分$v_{i}^{(M)},$$w_{i}^{(M)}(1\leq i\leq N)$は，11から有限回の四則演算によって得られる．ここに，$i,$$2\leq s\leq M$$J_{M}$

から$j,$ $p,$$q$は，整数 $1\leq i\leq N-1,2\leq i\leq N$,である．
$\square$$J_{M}^{z^{1}\pi}$を計算するためのコストは無視できるとすると，定理 6 に基づいて一般化 Newton$M$ は固定された定数であるから，を計算するために必要な演算回数は $O(M^{2}N)$ のオーダーである．
このオーダーは O(恥と見倣せる．

アンダーフロー，オーバーフローは生じないとの仮定の下で，$M=2,3$ の場合について丸め誤差またはの全ての対の蓄積を見積もる．注意 5 を考慮し，定理 6 に基づいて角成分を求める．このときの丸め誤差の蓄積による相対誤差は，計算中に現れる全変数に対し，マシ$((B^{T}B)^{M})^{-1}$ンイプシロンを$((BB^{T})^{M})^{-1}$として，最悪の場合でも以下のオーダーである．したがって，少なくと$M\geq 4$ の場合の解析
も $M=2,3$ のとき，求めるべき逆行列の対角成分の計算は数値安定性を持つ．
$\epsilon$$O(M^{2}N\epsilon)$

は今後の課題である．

5

まとめ上二重対角行列の最小特異値の下界として一般化 Newton 下界を導入した．この下界を計算するための公式を漸化式中に減算を含む形式と含まない形式の 2 通りの形式で与えた．これらの公式より一般化 Newton 下界を求めるのに必要な逆行列の対角成分が有限回の四則演算で得られる．

参考文献[1] K. V Femando and B. N. Parlett, Accurate singular values and differential qd algorithms, Numer.

Math. 67 (1994) 191-229.
[2] M. Iwasaki and Y. Nakamura, Accurate computation of singular values in terms of shifted inte-grable schemes, Japan J. Indust. Appl. Math. 23 (2006) 239-259.
[3] U. von Matt, The orthogonal qd-algorithm, SIAM J. Sci. Comput. 18 (1997) 1163-1186.

[4] K. Kilnura, T. Yamashita and Y Nakamura, On $O(N)$ formula for the diagonal elements of inverse

powers of symmetric positive definite tridiagonal matrices, 京都大学大学院情報学研究科数理工学専攻テクニカルレポート 2008-010, 2008