数理解析研究所講究録
第 1666 巻 2009 年 49-56

49

行列のスペクトル分解アルゴリズムについて
田島慎一
SHINICHI TAJIMA

新潟大学工学部
FACULTY OF ENGINEERING, NIIGATA UNIVERSITY

飯塚由貴恵
YUKIE IIZUKA

新潟大学工学部
FACULTY OF ENGINEERING, NIIGATA UNIVERSITY

1

はじめに行列のスペクトル分解は重要な概念であるが, 計算機代数の分野ではそれを正確に求めるアルゴリズムの研究は十分になされていないため, 行列のスペクトル分解アルゴリズムは既存の数式処理システムには実装されていない. 理論的には, レゾルベントの留数を用いることでスペクトル分解を表現する式を導けることが知られている. しかし, そのことをそのまま利用したプログラムでは逆行列の計算に時間がかかるため数式処理のプログラムとして実用性に乏しい.

本研究では, 最小多項式を基にスペクトル分解を計算する方法にっいて考察した. 行列多項式の計算を行うことで, 逆行列の計算を避けスペクトル分解の固有値の多項式表現を求めるアルゴリズムを導出し, さらに数式処理システム Risa/Asir に実装した.

2スペクトル分解

2.1

概念$A$を $n$ 次正方行列とし, 行列 $A$ の最小多項式をえる.

$f(\lambda)=0$の相異なる解を $\alpha_{j}(j=1,2,$$l$$\ldots,$$f(\lambda)$で表す. 以下$f(\lambda)$が重複因子を持たない場合を考ただし $l\leq n)$ とおくと, 各$\alpha_{j}$

に対して

(1)

$\{\begin{array}{l}P_{1}+P_{2}+\cdots+P_{l}=E\alpha_{1}P_{1}+\alpha_{2}P_{2}+\cdots+\alpha\iota P_{l}=A\end{array}$を満たす $n$ 次正方行列」らが存在する. ここで, $E$ は $n$ 次単位行列である.
ようにレゾルベント $(\lambda E-A)^{-1}$ の極$\alpha_{j}$$P_{j}$

は射影行列と呼ばれ, 以下のにおける留数として表現できる.

$P_{j}= \frac{1}{2\pi i}\oint_{\alpha_{j}}(\lambda E-A)^{-1}d\lambda$

ここで $(\lambda E-A)^{-1}$をレゾルベントという. 式 (1) が行列 $A$ のスペクトル分解である.

(2)

50

22

計算方法

221

準備

最小多項式

$f(\lambda)$

に対して$f(x)-f(y)=q(x, y)(x-y)$(3)を考える. 式 (3) に $x=A,$ $y=\lambda E$ を代入すると$f(A)-f(\lambda E)=q(A, \lambda E)(A-\lambda E)$を得る.

ここで $f(A)=0,$$f(\lambda E)=f(\lambda)E$(4)より$f(\lambda)E=q(A, \lambda E)(\lambda E-A)$(5)となる. 式 (5) の両辺に右から $(\lambda E-A)^{-1}$ を掛けると$f(\lambda)(\lambda E-A)^{-1}=q(A, \lambda E)$であり, さらに, 式 (6) の両辺を$f(\lambda)$(6)で割ると$( \lambda E-A)^{-1}=\frac{1}{f(\lambda)}q(A, \lambda E)$(7)となるので, 式 (2) に式 (7) を代入することで$P_{j}= \frac{1}{2\pi i}\oint_{\alpha_{j}}\frac{1}{f(\lambda)}q(A, \lambda E)d\lambda$(8)を得る ([1, 4]).

222

$q(A, \lambda E)$

最小多項式の計算$f(\lambda)$

を$f(\lambda)=\lambda^{l}+a_{l-1}\lambda^{l-1}+\cdots+a_{1}\lambda+a_{0}$で表されるモニックな$q(x, y)$

$l$

(9)

次多項式とすると, $q(x, y)$ は$=$$\frac{f(x)-f(y)}{x-y}$$=$$\frac{(x^{l}+a_{l-1}x^{l-1}+\cdots+a_{1}x+a_{0})-(y^{l}+a\iota-1y^{l-1}+\cdots+a_{1}y+a_{0})}{x-y}$$=$$\frac{(x^{l}-y^{l})+a_{l-1}(x^{l-1}-y^{l-1})+\cdots+a_{2}(x^{2}-y^{2})+a_{1}(x-y)}{x-y}$$=$$(x^{l-1}+x^{l-2}y+\cdots+xy^{l-2}+y^{l-1})$$+a_{1-1}(x^{l-2}+x^{l-3}y+\cdots+xy^{\iota-3}+y^{l-2})+a_{2}(x+y)+a_{1}$$=$$y^{l-1}+(x+a_{1-1})y^{l-2}+(x^{2}+a_{1-1}x+a_{l-2})y^{l-3}$$+\cdots+(x^{l-2}+a_{l-1}x^{l-3}+\cdots+a_{2})y+(x^{l-1}+a_{1-1}x^{l-2}+\cdots+a_{2}x+a_{1})$となる. 式 (10) に $x=A,$$q(A, \lambda E)$$=$$y=\lambda E$を代入すると, $q(A,$$\lambda E)$(10)は$E\lambda^{t-1}+(A+a_{l-1}E)\lambda^{l-2}+(A^{2}+a_{l-1}A+a_{l-2}E)\lambda^{l-3}$$+\cdots+(A^{l-2}+a_{1-1}A^{l-3}+\cdots+a_{2}E)\lambda+(A^{l-1}+a_{l-1}A^{l-2}+\cdots+a_{2}A+a_{1}E)(11)$51となる. さらに, 式 (11) の変数$q(A, \lambda E)$$\lambda$の次数ごとの係数行列に注目すると$E\lambda^{l-1}+(A+a_{l-1}E)\lambda^{l-2}+((A+a_{l-1}E)A+a_{l-2}E)\lambda^{l-3}$$=$$+\cdots+(\cdots(((A+a_{l-1}E)A+a_{l-2}E)A+al-3E)A+\cdots+a_{2}E)\lambda$$+((\cdots(((A+a_{l-1}E)A+a_{l-2}E)A+a_{l-3}E)A+\cdots+a_{2}E)A+a1E)$と変形でき,

$k$

(12)

次の係数行列は $k+1$ 次の計算結果を用いて表現できることがわかる. 従って, ホーナー法を用いることで係数行列を効率的に計算できる ([2, 6]).

223

導出

$f(\lambda)$は重複因子を持たないので,$f(\lambda)$と $f’(\lambda)$ は互いに素であることから(13)$a(\lambda)f(\lambda)+b(\lambda)f’(\lambda)=1$を満たす多項式$a(\lambda),$$b(\lambda)$が存在する. よって, 式 (8) は$P_{j}$$=$$\frac{1}{2\tau_{I}i}\oint_{\alpha_{j}}\frac{a(\lambda)f(\lambda)+b(\lambda)f’(\lambda)}{f(\lambda)}q(A, \lambda E)d\lambda$$=$$\frac{1}{2\pi i}\oint_{\alpha_{j}}(a(\lambda)q(A, \lambda E)+\frac{f’(\lambda)}{f(\lambda)}b(\lambda)q(A, \lambda E))d\lambda$(14)と書くことができる.

次に, 式 (14) の被積分関数中の特異性に注目する. 第一項は$\lambda=\alpha j$

で正則であり, 第二項は$\lambda=\alpha j$に極を持っことから(15)$P_{j} \equiv\frac{1}{2\pi i}\oint_{\alpha_{j}}(\frac{f’(\lambda)}{f(\lambda)}b(\lambda)q(A, \lambda E))d\lambda$を得る. さらに, 留数計算の公式

を用いると, 射影行列

$P_{j}$$\frac{1}{2\pi i}\oint_{\alpha}\frac{f’(x)}{f(x)}h(x)dx=h(\alpha)$(16)$P_{j}=b(\alpha_{j})q(A, \alpha_{j}E)$(17)

はで与えられることがわかる (cf. [3, 7]). このように, 各々の固有値を不定文字うに共通の表現となることから, 固有値式 (17) の右辺のをに置き換えたものに対して$\lambda$$\alpha j$$f(\lambda)$$\lambda$$\alpha j$$l$

次最小多項式$f(\lambda)$に対して$b(\lambda)$に対して射影行列乃は式 (17) のよに置き換えることにする.

modと $q(A, \lambda E)$ はともに $l-1$ 次多項式であることから $b(\lambda)q(A, \lambda E)$

きる.

作成したプログラム

.

3.1

(18)

$f(\lambda)$

は $2l-2$ 次多項式となるが, $f(\lambda)$ で割った余りを考えることで $P(A, \lambda)$ を

3

ここで $f(\alpha j)=0$ より,で割った余りを考え, これを $P(A, \lambda)$ とすると$P(A, \lambda)=b(\lambda)q(A, \lambda E)$となる.

$\alpha j$

プログラムの仕様
数式処理システム Risa/Asir ([5]) で動作するプログラム$\lambda$の高々 $l-1$ 次多項式で表現で

52

.
.
.
.

機能 : 最小多項式が重複因子を持たない場合の行列のスペクトル分解を計算
呼出し形式 : 関数名 (Mat, Poly)

引数 :Mat$arrow$

戻り値 :[$b(\lambda)$$P(A, \lambda)$

正方行列, Polyの共通分母,$arrow$$[\lambda^{1-1}$特性多項式の係数行列. . .

$\lambda$の係数行列定数項の係数行列]]が戻り値であり,$P(A, \lambda)=\frac{1}{b(\lambda)\text{の}\#_{\backslash }\text{通分母}}$(($\lambda^{l-1}$の係数行列)$\lambda\iota$-1 $+\cdot\cdot\cdot+$( の係数行列)$\lambda$のように戻り値の要素と対応している. 従って, 式 (19) の不定文宇

値$\alpha_{j}$32が定める固有ベクトル空間への射影行列$P_{j}$$\lambda$$\lambda+$

を固有値$\alpha j$(定数項の係数行列))(19)に置き換えることで, 固有が exact にわかる.

具体例$A=(\begin{array}{lll}-4 5 -4-3 3 3-1 -2 3\end{array}),$$E=(001$$001$$001^{\cdot}),$$f(\lambda)=\lambda^{3}-2\lambda^{2}+2\lambda+66$(20)に対してスペクトル分解を計算する. このとき
$f’(\lambda)=3\lambda^{2}-4\lambda+2$(21)$b( \lambda)=\frac{1}{60134}(2\lambda^{2}-303\lambda+202)$(22)$q(A, \lambda E)=(\begin{array}{lll}l 0 00 1 00 0 1\end{array})\lambda^{2}+(\begin{array}{lll}-6 5 -4-3 1 3-l -2 1\end{array})\lambda+(\begin{array}{lll}15 -7 276 -16 249 -13 3\end{array})$(23)であり, $P(A, \lambda)1$ ま$P(A, \lambda)$$=$$\frac{1}{60134}((\begin{array}{lll}1424 -1509 l250909 -731 -849317 572 -693\end{array}) \lambda^{2}+(\begin{array}{lll}-5267 31l1 -8973-2412 5512 -6678-2925 3543 -245\end{array}) \lambda+(\begin{array}{lll}23556 -2074 59821608 l6370 44521950 -2362 20208\end{array}))$(24)となる. 式 (24) に固有値$\lambda=\alpha_{1},$$\alpha_{2},$$\alpha_{3}$をそれぞれ代入すると射影行列$P_{1}$,$P_{2},$ $P_{3}$がわかり, スペクトル分解が求まる.

33

入出力例

32 節の行列に対する入出力例を示す. あらかじめ関数 spect(仮の関数名) をロードする必要がある.

53[14] $A=newmat$ (3, 3, $[[-4,5, -4],$ $[-3 , 3, 31 , [-1, -2,3]])$ ;$[0]$$[-45-4]$$[-333]$$[-1-23]$[15]$E=newmat(3,3, [[1,0,0], [0,1,0], [0,0,1]])$ ;$[100]$$[010]$$[001]$[16] $F=\det(x*E-A)$ ;$x^{-}3-2*x^{-}2+2*x+66$[17] spect $(A,F)$ ;$[60134,$ $[[1424-15091250]$$[909 -731-849]$$[ -52673111-8973]$$[-24125512-6678]$$[$317 572 $-693$$]$$[ -29253543-245 ]$[ 1608 16370 4452 1[23556-20745982 ][1950-2362 20208 ]$]]$[18]の分母であり, $P(A, \lambda)$ の分母と一致する. ま京都大学の木村欣司氏が作成した特性多項式計
た, 特性多項式を求める際に行列のサイズが大きい場合は,ここで, 関数 spect の出力結果内の 60134 は, 32 節の$b(\lambda)$

算プログラムを用いる.

4アルゴリズム
行列のスペクトル分解を求めるアルゴリズムの流れは以下の通りである.

1. 最小多項式$f(\lambda)$が重複因子を持たないことを確かめる. 重複因子を持っ場合は, エラーメッセージを表示し終了させる.

2.

$b(\lambda)$を求める.

3. $q(A, \lambda E)$ を求める.
4. $P(A, \lambda)=b(\lambda)q(A, \lambda E)$ mod
ここで, $q(A, \lambda E)$ の変数$\lambda$と$b(\lambda)$$f(\lambda)$を求め, $P(A, \lambda)$ を戻り値として出力する.

に対して$\lambda^{k+1}b(\lambda)$

mod$f(\lambda)=\lambda(\lambda^{k}b(\lambda)$

mod$f(\lambda))$mod$f(\lambda)$(25)における計算結果を用いることで, 次数 $k+1$ における計算を効率的に行うことができる. その他, 計算速度を向上させるために以下の工夫を行った.
より, 次数.
.

$k$の変数を扱わずに, 各次数毎の係数行列のみを計算する.

最小多項式$f(\lambda)$多項式の係数が有理数になることに注目し,$b(\lambda)$

式に分けた ($b(\lambda)$$b(\lambda)$

を係数の共通分母となる整数と, 整数係数多項は [整数係数多項式, 共通分母] の形で扱われている).

54.
.

$q(A, \lambda E)$$\lambda^{k}b(\lambda)$の係数行列の計算等, 行列多項式計算を行う場合はホーナー法を利用する (2.2.2 節参照).

mod$f(\lambda)$を先に計算し, その結果に $q(A, \lambda E)$ の係数行列を掛けることで, 乗算の回数を少なくした.

5

計算時間の測定・考察
行列のスペクトル分解を求めるプログラムの CPU 時間, GC 時間, 経過時間を測定した.

5.1.
.

測定環境
使用ソフト

実験で使用したコンピュータ
-OS-

CPU

Windows XP$\ldots$$\ldots$

一メモリ

52.
.
.
.

Risa/Asir

数式処理システムIntel(R) Pentium(R) 4 $(2.52GHz)$

2OOGB

測定方法
行列の要素は整数であり,-512$\sim$

511 (3 桁以内), -65536 65535 (5 桁以内), -8589934592 8589934591$\sim$(10 桁以内) の 3 パターンで実験$n$次正方行列に対して 10 回ずつ測定し, それぞれの平均をとる

CPU 時間, GC 時間, 経過時間は Asir の time

Asir の設定はデフォルトのまま

$()$

で取得

$\sim$

5553

測定結果
表 1: スペクトル分解の計算にかかる時間 (行列の要素 :3 桁以内)

行列の要素が 3 桁以内のときは $n=53,5$ 桁以内のときは $n=43,10$ 桁以内のときは $n=42$ を超えると“Too many heap sections” と表示され, 計算ができなかった.

54

考察

行列の要素の桁数を変えて実験を行ったところ, $n=40$ において要素が 10 桁以内のときの経過時間は, 5桁以内のときと比較して約 34 倍かかった. 時間こそかかるものの, 計算可能な行列のサイズは 5 桁以内のときとそれほど差が無い.

56

次に経過時間に占める CPU 時間の割合に注目する. 3 つの表を比べると, 行列のサイズが小さいときは行列の要素の桁数が少ないときの方が割合が高いが, サイズが大きくなると桁数の多い場合に CPU 時間の割合が高くなる傾向がある.

6

まとめ逆行列の計算を避けることで, スペクトル分解を求めることが可能になった. また, 今回作成したプログ
ラムは固有値の多項式として射影行列を表現する式を求めており, 固有値の値を代入することでスペクトル分解を正確に求めることが可能となっている.
2009 年 1 月, 金沢大学の小原功任氏との共同研究により計算効率の良い新たなアルゴリズムを導出し,Risa/Asir に実装した. また 2009 年 3 月には, 最小多項式が因数分解できる場合 (重複因子を持たない), そのことを利用してスペクトル分解を求めるアルゴリズムを実装した. 以上 2 つのアルゴリズムについては別の機会に詳しく述べることにする.

7

今後の課題
本アルゴリズムは最小多項式が重複因子を持たない場合に限定していたので, 今後最小多項式が重複因子を持つ場合のスペクトル分解の計算法について研究をする予定である. この場合はレゾルベントが 2 位以上の極を持っため本アルゴリズムをそのままの形で適用することはできない. したがって, 本アルゴリズムを拡張していく必要がある.

本研究は, 平成 20 年度内田エネルギー科学振興財団の研究助成を受けている.

参考文献[1] F. シャトラン (伊理正夫伊理由美訳) : 行列の固有値新装版最新の解法と応用, シュプリンガー
フェアラーク東京 (2003).
[2] F. R. Gantmacher : The Theory of Matrices, vol.1, 2, Chelsea (1960).
[3] 加藤涼香, 田島慎一 : 有理関数のローラン展開アルゴリズムと代数的局所コホモロジー, 京都大学数理解
析研究所講究録 1395 「Computer Algebra-Design of Algorithms, Implementations and Applications」(2004), pp.50-56.

[4] T. Kato: A Short Introduction to Perturbation Theory for Linear Operators, Springer-Verlag (1982).
[5] M. Noro and T. Takeshima, Risa/Asir-a computer algebra system, ISSAC 1992 (ed. P. S. Wang),ACM (1992), pp.387-396.

[6] 庄司卓夢, 田島慎一 : 多項式剰余公式の計算アルゴリズム, 京都大学数理解析研究所講究録 1514
「Computer Algebra- Design of Algorithms, Implementations and Applications」 (2006), pp.108-114.
[7] 田島慎一 : 多変数補間問題とホロノミック $D$ -加群, 千葉大学数学セミナーノート No 3 「代数解析学の諸問題」 (1999), pp.73-94.