特異値計算アルゴリズム，dqds 法および mdLVs 法の収束理論
66201

相島健助

指導教員

杉原正顯 教授
仮定 (A). B の対角，副対角成分はすべて正で ak >0 (k = 1, . . . , m),

1 はじめにbk > 0 (k = 1, . . . , m − 1)．

このとき特異値はすべて異なり σ1 > · · · > σm > 0 と行列の特異値計算では，まず与えられた行列を適切おける．

な直交変換により上２重対角化し，この上２重対角行列に対して，ある種の反復計算を行うことで特異値を求めるのが標準的である．
上２重対角行列の特異値計算アルゴリズムとして，近年，Fernando–Parlett により提案された dqds (differ-ential quotient difference with shifts) 法 [1] が注目されている．この dqds 法は数値安定性に優れ，高速に特2 dqds 法について
2.1 dqds 法に対応する行列変形
dqds 法は，B(0) = B と初期設定し，上２重対角行列
B(n) に対する反復式(B(n+1) )T B(n+1) = B(n) (B(n) )T − s(n) I異値計算できるので，線形計算用ライブラリ LAPACKにおいて DLASQ として実装されている．一方 2000
年頃から，中村らのグループにより mdLVs (modifieddiscrete Lotka-Volterra with shifts) 法と呼ばれるアルゴを計算することに相当する．dqds 法の具体的なアルゴリズムは以下のとおりである．

リズムが提案されている [2]．こちらも dqds 法と並ぶアルゴリズム 2.1 dqds 法優れた性能が報告されている．両アルゴリズムでは，初期設定: qk(0)(0)

いたわけではなく，収束速度に関しても，複雑な仮定= ak 2 (k = 1, 2, . . . , m); ek
1, 2, . . . , m − 1)
1: for n := 0, 1, · · · do
2:
シフト量 s(n) (≥ 0) の設定
(n+1)
(n)
3:
d1
:= q1 − s(n)4:for k := 1, · · · , m − 1 do
(n+1)
(n+1)
(n)
5:
qk
:= dk+ ekのもとで，局所収束次数が 2 次であると主張されてい6:

るものの，どのような場合に仮定が満たされるか明ら

7:

ある変数を初期設定して反復計算を行い，その変数の収束先から特異値を求めるが，収束加速のため各反復でシフトというものを適切に設定する必要がある．しかし，そもそも収束自体が完全な形で理論保証されて8:かではない．

9:

本論文では，上記の dqds 法と mdLVs 法の理論的収(1)(n+1)

ek

(n+1)

(n) (n)= bk 2 (k =(n+1):= ek qk+1 /qk(n+1) (n)(n+1)
qk+1 /qk− s(n)

dk+1 := dk
end for(n+1)(n+1)
qm:= dm

10: end for

束性を明らかにした．まず両アルゴリズムの厳密な収束証明を与えた．これにより，シフトに関する一般的な仮定の下，任意の初期行列で収束が保証されることになる．そして，シフト設定として具体的に Johnson
シフトを用いる場合，反復の終盤で常に 1.5 次収束することを示した．
以下，特異値を求めたい m × m の上２重対角行列 Bの各成分を

a10B= .
 ..
0b1a2..
.
...

0..
..

.

.
0
bm−1 amアルゴリズム中の変数と (1) の B(n) は

√(n)
q1

 0(n)
B = ..
 .
0

√
(n)
e
√ 1(n)
q2..
.

..

.

..

.

...

00√(n) e√ m−1 (2)(n)qmという関係式で結ばれる．このアルゴリズムによる特(N)異値の求め方を簡単に述べると，ある N で em−1 ≈ 0となったときに収束したと判定し，一旦アルゴリズ(N)(l)ムを止める．するとこのとき σm 2 ≈ qm + ∑N−1
l=0 sが成り立つので，最小特異値 σm の近似値が求まる．

とおく．ここで [1] に従って一般性を失わず，次のこそしてデフレーションを行うと，すべての特異値をとを仮定する．

σm , σm−1 , . . . , σ1 の順に計算できる．

である．式中の δ (n) は正の範囲で自由に選べるパラ

2.2 主結果：dqds 法の収束定理(n)メータである．ここで B(n) の副対角成分の２乗を ek ，(n)(n)
これまでの議論で，em−1 が 0 に収束すれば特異値(n)が求まることを説明した．本研究では，em−1 を含めてdqds 法の各変数の収束が理論的に保証されることを示(n)した．次の定理中の σmin は B(n) の最小特異値である．
定理 2.1 (dqds 法の大域的収束性). 行列 B が仮定 (A)を満たすとする．dqds 法におけるシフト量 s(n) を，(n)0 ≤ s(n) < (σmin )2 とすれば，(n)lim e = 0n→∞ k(k = 1, 2, . . . , m − 1),(n)lim q = σk 2 −n→∞ k∞∑ s(n)(k = 1, 2, . . . , m)

3.2 主結果：mdLVs 法の収束定理
mdLVs 法においても，特異値を求めるには dqds 法(n)と同様に em−1 を 0 に収束させたい．本研究ではこの収束が理論的に保証されることを示した．
定理 3.1 (mdLVs 法の大域的収束性). 行列 B が仮定 (A)を満たすとする．mdLVs 法アルゴリズムにおいて，パラメータ δ (n) > 0 を，ある非負実数 D0 ≥ 0 に対して(n)ト量 s(n) を 0 ≤ s(n) < (σmin )2 となるようにとれば，次に，Johnson 下界 [3] を利用してシフトを{√(√√ )}1(n)(n)(n)λ = minqk −ek−1 + ek,
k=1,...,m2()2s(n) = max{λ (n) , 0}(3)(n)(n)と定めることにする．このとき em−1 は 1.5 次収束することを示した．
定理 2.2 (1.5 次収束). 行列 B が仮定 (A) を満たすとする．このときシフト量を上記 (3) で決定する dqds 法において，em−1(n)(em−1 )3/2(n)lim e = 0 (k = 1, 2, . . . , m − 1),n→∞ k∞(n)lim qk = σk 2 −s(n) (k = 1, 2, . . . , m)n→∞n=0∑が成立する．
(n)

さらに，具体的にシフト設定に (3) を用いると，em−1は 1.5 次収束することを示した．
定理 3.2 (1.5 次収束). 行列 B が仮定 (A) を満たし，パラメータ δ (n) は，ある正の実数 D1 > 0 に対して(n+1)

limn→∞ δ (n) = +∞ かつ limn→∞ δ (n) em−1 = D1 を満たすように設定する．このときシフト量を (3) で決定す(n+1)

lim(n)を σmin とする．

limn→∞ 1/δ (n) = D0 を満たすようにとり，さらにシフn=0が成立する．

n→∞対角成分の２乗を qk とする．また B(n) の最小特異値1=√,
σm−1 2 − σm 2(n+1)
qm1
lim (n)=√n→∞σm−1 2 − σm 2(qm )3/2が成り立つ．

る mdLVs 法において，(n+1)

limn→∞

em−1(n)(em−1 )3/2(n+1)

lim

qm(n)
n→∞(qm )3/2=√

1σm−1 2 − σm 2,1=√σm−1 2 − σm 2が成り立つ．

3 mdLVs 法について
3.1 mdLVs 法に対応する行列変形
mdLVs 法 は 可 積 分 系 に 基 づ く 特 異 値 計 算 法 で ある [2]．dqds 法に相当する反復式を行列の形で (1) と表現できたが，mdLVs に相当する反復式を行列の形で書くと，

1I,δ (n)1(B̂(n) )T B̂(n) = B̄(n) (B̄(n) )T − (n) I,
δ(n+1) T (n+1)(n) T (n)(B) B= (B̂ ) B̂ − s(n) I(B̄(n) )T B̄(n) = (B(n) )T B(n) +参考文献[1] K. V. Fernando, B. N. Parlett: Accurate Singular Values and Differential qd Algorithms, Numerische Mathematik, vol. 67 (1994), pp. 191–229.
[2] M. Iwasaki and Y. Nakamura: Accurate Computation ofSingular Values in Terms of Shifted Integrable Schemes,
Japan Journal of Industrial and Applied Mathematics,
vol. 23 (2006), pp. 239–259.
[3] C. R. Johnson: A Gersgorin-type Lower Bound for theSmallest Singular Value, Linear Algebra and Its Applications, vol. 112 (1989), pp. 1–7.