行列関数計算のためのクリロフ部分空間法
数理第 2 研究室 佐藤 暁史

指導教員 室田 一雄 教授

2009 年 2 月 2 日
最後に，以上に述べたアルゴリズムにおいて，赤1 はじめに池のアルゴリズムを実際に実装し，解析結果を確か

本研究では，行列 A ∈ Cn×nとベクトル b ∈ C ，nめる．

そして f (A) が定義されるような関数 f : C ⊂ D →2 行列関数

C が与えられているときに，
f (A)b

Λ(A) = {λ1 , λ2 , . . . , λk } を，行列 A ∈ Cn×n(1)の相違なる k 個の固有値とする．複素関数 f が与えられているとき，行列 f (A) は，f (r) (λj ) がを近似的に評価する方法を扱う．この計算の応用に

r = 0, 1, . . . , nj − 1; j = 1, 2, . . . , k において存在す

おいては，多くの場合行列 A は大規模な疎行列である場合に定義される. この場合，qf,A ∈ PK−1 がたり，なんらかの構造を有している．この場合，最初

かだか K − 1 次の一意に定まる多項式であるとしに f (A) の計算をすることによって式 (1) を評価すて f (A) := qf,A (A) と表すことができる．ここで，るのは不可能であることが多く，よって残りの計算(r)

qf,A (A) ∈ PK−1 は，qf,A (λj ) = f (r) (λj ),をするために利用可能な様々なアルゴリズムを使う

r=0, 1, . . . , nj − 1, j = 1, 2, , . . . , k を満たす一意に定ことができない．f (A) の計算を行わず，式 (1) をまる多項式とする．

直接的に得るための標準的なアプローチは，初期ベクトル b による A のクリロフ部分空間に基づいた3 再出発クリロフ部分空間近似方法である．この方法による長所は，行列とベクトAルの掛け算にのみ A が必要となるということと，f∈Cn×n と 0

̸=b∈

Cn で 与 え

ら れ る ク リ ロ フ 部 分 空 間 を Km (A, b)が滑らかな関数ならば超一次収束をするということ

m−1

span{b, Ab, . . . , Aである．しかし，クリロフ部分空間による近似には:=b} と定義する．f (A)b の計算のための再出発クリロフ部分空間法は以下のよう一つの欠点がある．それは，次元数だけのベクトルになる．

をすべて保存する必要があるということである．ゆえに，しばしばメモリの制限によって，解ける問題

Given: A, b, fのサイズが限定されてしまう．クリロフ部分空間近f (0) := f, f (0) := 0, b(0) := b, γ (0) := ∥b∥

似による非エルミートの線形方程式系の一般的な解for k = 1, 2, . . . until convergence do(k)Km (A, b(k−1) ) の Arnoldi 分 解 AVm法は，クリロフ空間がある最大の次元に達するたび=(k) (k)(k)(k)Vm Hm + µm+1,m v m+1 e⊤m を計算する．
(k)にアルゴリズムを再出発して保存領域を制限することである．Eiermann ら [3] は，有限個のベクトル近似 fを記憶するだけで計算できる再出発クリロフ部分空end for間近似を提案した．

を更新する．

4 再出発長 1 の場合とその性質

Eiermann らのアルゴリズムのうちの特殊な場合
は Afanasjew ら [1] によって最急降下法の一般化にA を正則な n × n 実数行列，x, b を n 次元 実数なっていることが示されている．Afanasjew らの解ベクトル，x を変数とする連立一次方程式 Ax = b

析が赤池 [2] による解析の一般化になっていることの解を x∗ とする．Ax = b を解く最急降下法を以も述べる．

下の手順により定義する：

1

Step 0 x0 を適当な n 次元ベクトルとし，k = 0 と5 数値計算する．

Step 1 Aϵk = Axk − b とする．

以下の値を用い，線型方程式を解くための最急降⊤下法について，赤池の性質を確かめる数値計算を

Step 2 ζk = A Aϵk とする．
Step 3ϵ⊤ A⊤ Aζ
γk = ζk⊤ A⊤ Aζk とする．
kk行った．

 √Step 4 xk+1 = xk − γk ζk とする．

0.00268704
 √
0.01581310
√
0.08234830• A=√

0.17590130
√
0.25946632
√0.49823436• b = 0, x0 = (1, 1, 1, 1, 1, 1)⊤

Step 5 終了条件を満たせば終了．



0Step 6 k の値を 1 増やして Step 1 へと戻る．
これに関して，赤池によって，以下の定理が示され0ている．
定理 4.1.

1. ϵk は A⊤ A の 2 つの固有ベクトルの線形結合

結果の一部を図 1 に示す．に近づいていく．

60

2. ϵk は漸近的に 2 つの特定の方向へ近づく．

50

この性質をみると，ζk と ζk−2 がほとんど同じ
Gamma_k40

方向を向いているとき，xk−2 − xk の方向を用いて xk+1 を計算すると良いことがわかる．すなわち，γ̂k =⊤ϵ⊤k A A(xk−2 −xk )(xk−2 −xk )⊤ A⊤ A(xk−2 −xk )30

2010とし，xk+1 =xk − γ̂k (xk−2 − xk ) によって，xk+1 を計算する．

00510152025
k30

35404550

これを加速手続きとする．
図1

再出発長 1 の場合，再 出 発クリロフ部分空間

γk の値

法は

Given: A, b

参考文献

σ1 := ∥b∥, v1 := b/σ1

[1] M. Afanasjew, M. Eiermann, O. G. Ernst, andfor k = 1, 2, . . . doS. Güttel. A generalization of the steepest de-w := Av k ; ρk := v Hk w

scent method for matrix functions. Electronicw := w − ρk v k ; σk+1 := ∥w∥

Transactions on Numerical Analysis, Vol. 28,

v k+1 := w/σk+1

pp. 206–222, 2008.

end for[2] H. Akaike. On a successive transformation ofとなる．これに関して，次の定理が示された．

probability distribution and its application to定理 4.2. A：をエルミート行列，λmin を最小固有

the analysis of the optimum gradient method.

値，λmax を最大固有値とし，b は λmax や λmin にAnnals of the Institute of Statistical Mathe-対応する固有ベクトルを含むとする．さらに θ ∈matics, Vol. 11, pp. 1–17, 1959.

(0, 1) であるとする．この時以下が成立する．

[3] M. Eiermann and O. G. Ernst. A restarted

lim ρ2k−1 =θλmin + (1 − θ)λmax =: ρ∗1 ,krylov subspace method for the evaluation oflim ρ2k =(1 − θ)λmin + θλmax =: ρ∗2 ,k→∞√lim σk+1 = θ(1 − θ)(λmax − λmin ) =: σ ∗ .

matrix functions.

k→∞

SIAM J. Numer. Anal.,Vol. 44, No. 6, pp. 2481–2504, 2006.

k→∞

これは赤池の結果の一般化である．

2