数理解析研究所講究録
第 1815 巻 2012 年 21-28

21

最小消去多項式を用いた行列スペクトル分解計算の並列化
小原功任
KATSUYOSHI OHARA

金沢大学理工研究域数物科学系*
田島慎一
SHINICHI TAJIMA

筑波大学大学院数理物質科学研究科\dagger

1

行列スペクトル分解
正方行列 $A$ の固有値を$\lambda_{1},$$\cdots,$$\lambda_{m}$とし，固有値いるように、 このとき互いに可換な行列$P_{1},$$\cdots$,$\lambda_{i}$に関する不変部分空間を$P_{m},$ $D$$V_{i}$とする。 よく知られてが存在して$E=P_{1}+\cdots+P_{m},$$A=\lambda_{1}P_{1}+\cdots+\lambda_{m}P_{m}+D,$$P_{i}$$D$は鷲への射影行列，は巾零行列，となる。 この $\{P_{1}, \cdots, P_{m};D\}$ を行列 $A$ のスペクトル分解という。

本研究の目的は、 $Q$ 上の行列のスペクトル分解をできるだけ高速に exact に計算することである。 そのために、数学的にはレゾルベントに関する留数解析を、 計算技術としては並列計算の手法を用いる。 また、$D$ は固有値上の行列のスペクトル分解であるから、 一般に行列の多項式または有理式を成分と$Q$$P_{i},$$\lambda_{i}$する。 このときに冗長な表現で求めてしまうと、 一見して高速に計算できたように思えても、 スペクトル分解の後処理としての simplification あるいは一意化の手間が馬鹿にならずかえって時間がかかってしまう。

したがって本研究では最小次数の多項式を成分とするような結果を求めるものである。本稿でははじめにスペクトル分解アルゴリズムの概略を、 その後に今回工夫を行った最小消去多項式を利用した高速計算法の概略を述べ、最後に並列計算アルゴリズムについて説明する。 また、 われわれのアルゴリズムは計算代数システム $Risa/Asir$ に実装しているが、 そのタイミングデータについても述べる。

2レゾルベントと射影行列

定義 1. 正方行列 $A$ に対して、 行列値複素関数 $R(z)=(zE-A)^{-1}$ を $A$ のレゾルベントという。
*ohara@air.s.kanazawa-u.ac.jp\dagger tajima@math.tsukuba.ac.jp

22の特性
簡単に分かるように、 レゾルベント $R(z)$ は $C$ 上の行列値有理関数となる。 また $R(z)$ の極はレゾルベントについては次の定理がよく知られすなわち固有値である。の根、多項式 $x(z)=\det(zE-A)$$A$ている。

定理 1. $A$ の固有値を$\lambda_{1},$$\cdots,$$\lambda_{m}$とする。

$\lambda_{i}$におけるローラン展開$R(z)= \cdots+\frac{1}{(z-\lambda_{i})^{2}}D_{i}+\frac{1}{z-\lambda_{i}}P_{i}+\cdots$

は、 $A$ のスペクトル分解 $\{P_{1}, \cdots,P_{m};\sum_{i=1}^{m}D_{i}\}$ を与える。

したがってスペクトル分解を求めるには、 $R(z)$ の固有値のまわりでのローラン展開について調べればよいことになる。特に、コーシーの積分定理から次がしたがう。
定理 2. 固有値

影行列

$P_{i}$$z=\lambda_{i}$のまわりを反時計回りに一周まわる閉曲線をは、周回積分$P_{i}= \frac{1}{2\pi\sqrt{-1}}\int_{c_{:}}R(z)dz$$C_{i}$とする。 このとき、$\lambda_{i}$に対する射によって与えられる。

また、 レゾルベントの有理関数としての構造については次の定理が得られる。

定理 3. 多項式 $f(z)\in C[z]\backslash \{0\}$ が $f(A)=O$ を満たせば、 レゾルベントは$R(z)= \frac{1}{f(z)}\Psi_{f}(zE,A)$,

ただし$\Psi_{f}(x,y)=\frac{f(x)-f(y)}{x-y}\in C[x,y].$と表される。

例行列$A=(\begin{array}{lll}0 4 0-1 4 00 0 3\end{array})$である。 このとき、$A$の最小多項式は $\pi(x)=(x-2)^{2}(x-3)$ であり、 固有値は$\lambda_{1}=2,$のレゾルベントは$R(z)= \frac{1}{(z-2)^{2}}(-6E+5A-A^{2})+\frac{1}{z-2}(-3E+4A-A^{2})+\frac{1}{z-3}(4E-4A+A^{2})$と書かれ、その係数行列$P_{1}=-3E+4A-A^{2}=(\begin{array}{lll}1 0 00 1 00 0 0\end{array}), P_{2}=4E-4A+A^{2}=(\begin{array}{lll}0 0 00 0 00 0 1\end{array}),$$D_{1}=-6E+5A-A^{2}=(\begin{array}{lll}-2 4 0-1 2 00 0 0\end{array}), D_{2}=O$

は $A$ のスペクトル分解 $\{P_{1},P_{2};D_{1}+D_{2}\}$ を与える。

例 2. 行列$A=(\begin{array}{llll}0 2 0 11 0 0 00 0 0 20 0 1 0\end{array})$$\lambda_{2}=3$

23の最小多項式は $\pi(x)=(x^{2}-2)^{2}$ であり、 固有値は$\lambda_{1}=\sqrt{2},$$\lambda_{2}=$-西である。 このとき、$A$のスペク

トル分解 $\{P_{1}, P_{2};D_{1}+D_{2}\}$ は$P_{i}=(\begin{array}{llll}\frac{1}{2} \frac{\lambda}{2} 0 \lambda\vec{8}\lambda\vec{4} \frac{1}{2} -\lambda\vec{16} 00 0 \frac{1}{2} \lambda\vec{2}0 0 r\lambda 4 \frac{1}{2}\end{array}), D_{i}=(\begin{array}{llll}0 0 \lambda\vec{8} 410 0 \frac{1}{8} arrow\lambda 80 0 0 00 0 0 0\end{array}) (i=1,2)$で与えられる。

3レゾルベントの留数解析とネーター作用素
定理 3 より行列 $A$ を消去する多項式があれば、 レゾルベントを具体的に書き下すことができるが、 計算量的観点からは、 次数が小さな多項式を選ぶべきである。 したがってレゾルベントは最小多項式 $\pi(z)$ を用いて $R(z)= \frac{1}{\pi(z)}\Psi_{\pi}(zE, A)$ と表されることになる。 さらに、 定理 2 より、 レゾルベントの周回積分で射影行列は表されるので、 最小多項式で表されたレゾルベントについて、 周回積分を exact に求めることができれば、 射影行列の計算はできることになる。 そのために代数的局所コホモロジー群とネーター作用素の概念を導入する。 ネーター作用素計算の $Risa/$ Asir への実装については、 [4] を参照のこと。本節では、 射影行列

$P_{i}$についてのみ説明するが、 巾零部分

$D_{i}$も同様に計算可能である。

$f(z)\in Q[z]$ を既約多項式とする。集合 $V(f)=\{z|f(z)=0\}$に同値関係$\equiv$に特異点をもつ有理関数の全体 $Q[z, f^{-1}]$を次のように入れる。

$r_{1}(z)\equiv r_{2}(z)\Leftrightarrow r_{1}(z)-r_{2}(z)\in Q[z]$この同値類を$[r_{1}]$で表し、 $V(f)$ に台をもつ代数的局所コホモロジー類という。 また，この同値類の全体を$V(f)$ に台をもつ代数的局所コホモロジー群という。 このとき、 与えられた有理関数 $g(z)/f(z)^{\ell}\in Q[z, f^{-1}]$に対して，次の多項式係数 $\ell-1$ 階微分作用素$[ \frac{g(z)}{f(z)^{\ell}}]=T[\frac{f’(z)}{f(z)}],$$T$を計算することができる ([3], [4])。

$T= \sum_{k=0}^{l-1}(-\partial_{z})^{\ell-k}t_{k}(z)\in Q\langle z,$

この $T$ をネーター作用素という。 より一般には、 $V(fi\cdots f_{m})$$\partial_{z}\rangle.$に台をもつ代数的局所コホモロジー群を考えることができ、 有理関数 $r(z)\in Q[z, f_{1}^{-1}, \cdots , f_{m}^{-1}]$ の代数的局所コホモロジー類に対して、 既約成分ごとのネーター作用素表示$[r(z)]=T_{1}[ \frac{f_{1}’(z)}{f_{1}(z)}]+\cdots+T_{m}[\frac{f_{m}’(z)}{f_{m}(z)}],$$T_{i}\in Q\langle z,$$\partial_{z}\rangle$が求まる。
行列 $A$ の固有値$\lambda_{i}$に対する射影行列$P_{i}$を求めるために、 ネーター作用素が使えることを示そう。 いま、 $A$ の最小多項式を $\pi(z)=fi(z)^{\ell_{m}}\cdots f_{m}(z)^{\ell_{m}}$ とする。 $f_{i}(z)$ は固有値

ネーター作用素$[ \frac{1}{\pi(z)}]=T_{1}[\frac{f_{1}’(z)}{f_{1}(z)}]+\cdots+T_{m}[\frac{f_{m}’(z)}{f_{m}(z)}]$$\lambda_{i}$の最小多項式である。 まず、24を求める。 ネーター作用素$T_{i}= \sum_{k=0}^{\ell_{:-1}}(-\partial_{z})^{\ell_{i}-k}t_{k}(z)$に対して、 随伴作用素$T_{i}^{*}= \sum_{k=0}^{\ell_{i}-1}t_{k}(z)(\partial_{z})^{\ell_{i}-k}$は容易に得られる。部分積分と留数定理により$=$$P_{i}$$\frac{1}{2\pi\sqrt{-1}}\int_{c_{:}}R(z)dz=\frac{1}{2\pi\sqrt{-1}}\int_{C_{i}}\Psi_{\pi}(zE, A)\frac{1}{\pi(z)}dz$$= \frac{1}{2\pi\sqrt{-1}}\int_{C_{1}}\Psi_{\pi}(zE, A)\sum_{j=1}^{m}T_{j}\frac{f_{j}’(z)}{f_{j}(z)}dz=\frac{1}{2\pi\sqrt{-1}}\int_{C_{i}}\Psi_{\pi}(zE, A)T_{i}\frac{f_{i}’(z)}{f_{i}(z)}dz$$= \frac{1}{2\pi\sqrt{-1}}\int_{C_{l}}\{T_{i}^{*}\Psi_{\pi}(zE, A)\}\frac{f_{i}’(z)}{f_{i}(z)}dz=\{T_{i}^{*}\Psi_{\pi}(zE, A)\}|_{z=\lambda_{i}}$このことから 2 変数多項式$(T_{i}^{*}\Psi_{\pi})(z, y)$$\Psi_{\pi}(z, y)$に対して、 微分作用素についての行列多項式昂$(\lambda_{i}E, A)$$\tau_{i}*$を作用させ、 得られた多項式 $F_{i}(z, y)=$が射影行列疏になることが分かった。

4最小消去多項式を用いた高速計算法とその並列化

4.1

最小消去多項式とスペクトル分解アルゴリズムの関係

定義 2. 正方行列$A$とベクトル $v\neq 0$ に対し、 $f(A)v=0$ を満たす次数最小のモニック多項式を、$v$の最小消去多項式という。
定理 3 の多項式 $f(z)$ を最小消去多項式で置き換えることで、 次の定理が得られる。

定理 4. 第$i$

基本ベクトル $ej$ の最小消去多項式を $\pi j(z)$ とする。 このとき、 $A$ のレゾルベント $R(z)$ について$R(z)e_{j}= \frac{1}{\pi_{j}(z)}\Psi_{j}(zE, A)e_{j}$,ただし$\Psi_{j}(x, y)=\frac{\pi_{j}(x)-\pi_{j}(y)}{x-y}\in Q[x, y].$が成り立つ。
$A$の最小多項式 $\pi(z)$ と最小消去多項式 $\pi j(z)$ を比較すると、 $\deg\pi j\leq\deg\pi$ であることは容易に分かる。 つまり最小多項式を最小消去多項式に置き換えることで、 アルゴリズムをを高速化することができる。

鳥を求めるための新しいスペクトル分解アルゴリズムは次のようになる。
アルゴリズム 1.
入力: $Q$ 上の正方行列 $A$

出力: スペクトル分解 $\{P_{i};D_{i}\}$

1.

$A$

2. 各の特性多項式$e_{j}$$\chi(z)$を求める。の最小消去多項式$\pi j(z)=f_{1}^{\ell_{1j}}\cdots f_{m}^{\ell_{mj}}$

3. $[1/\pi j(z)]$ のネーター作用素$T_{1j},$$\cdots,$$T_{mj}$とその既約分解を求める。

を求める。

4. 多項式 $F_{ij}(z, y)=T_{ij}^{*}\Psi_{j}(z, y)mod f_{i}(z)$ を求める。

5. 行列多項式 $Piej=F_{ij}(\lambda_{i}E, A)e_{j}$ を求める。
6. 乃を得る。

25
前節と同様にアルゴリズム 1 における計算量を見積もると、となる。

$\frac{1}{n}\sum_{j=1}^{n}\deg\pi_{j}\leq\deg\pi$$P_{i}$を求めるためには、 $O(n^{2} \deg f_{i}\sum_{j}^{n}=1\deg\pi_{j})$であるから、 アルゴリズム 1 のほうが計算量が小さい。 この差が大きければ大きいほど高速化されたことになる (ブロック対角化されている場合など)。

アルゴリズムの並列化の観点から見ると、 ステップ 3 については、 列番号でそれぞれ独立に計算でき、および列番号 でそれぞれ独立に計算できることが分かる。 さらに後述するようにステップ 2 における最小消去多項式の構成においても並列化技法が使用できる。 したがってアルゴリズム 1 は並列化の面からは非常に扱いやすい。
ステップ 4,5 については、 既約因子番号4.2$i$$j$$i$

最小消去多項式導出の逐次的計算法

本小節では、 第基本ベクトルに関する最小消去多項式の逐次的計算法について概説する。 アルゴリ
ズムの詳細については [8] を参照すること。
$i$$e_{j}$特性多項式の既約分解が $\chi(z)=fi(z)^{\ell_{1}}\cdots f_{m}(z)^{\ell_{m}}$ であれば、 最小消去多項式は$\pi_{j}(z)=f_{1}(z)^{l_{1j}}\cdots f_{m}(z)^{\ell_{mj}},$$(0\leq l_{ij}\leq\ell_{i})$の形になる。つまり、 既約因子は分かるため指数ベクトルの手順で求めることができる。
アルゴリズム 2. $k=1,$

1.

$\cdots,$ $m$$\chi_{k}(z)=\chi(z)/f_{k}(z)^{\ell_{k}}$および $j=1,$

とおき$\grave{}$$\cdots,$ $n$$b_{kj}=\chi_{k}(A)e_{j}$$(\ell_{1j}, \cdots , \ell_{mj})$が問題である。 指数ベクトルは次について次を行う。
をすべて求める。

2. $b_{kj}=0$ ならば、 $l_{kj}=0$ である。 $b_{kj}\neq 0$ ならば、 順に $f_{k}(A)$ を左からかけて、 $f_{k}(A)^{p}b_{kj}=0$ となる最小の $p>0$ が $\ell_{kj}=p$ である。
さて、 この方法で最小消去多項式が計算できることは分かったが、 問題点は計算量が大きいことである。

したがって、確率的アルゴリズムを導入して高速化を図る。確率的アルゴリズムのポイントは、 どんなベク
トル$v$に対しても$\pi_{j}(A)e_{j}=0$ならば$t_{v\pi_{j}(A)e_{j}=0}$が成り立つという事実である。 もちろん逆は成り立たないが、 乱数ベクトル をとり、 $t_{v\pi_{j}(A)e_{j}=0}$ で最小消去多項式を探索すれば、 ほぼ確率 1 で最小消去多項式を発見できることが期待される。 また、 $f_{k}(A)e_{j}$は次元ベクトルだが、 $tvf_{k}(A)e_{j}$ はスカラー量なので、 計算量も $1/n$ となる。
$v$$n$

アルゴリズム 3. 乱数ベクトル$v$を選び、 $k=1,$$\cdots,$$m$について次を行う。

1. $\chi_{k}(z)=\chi(z)/f_{k}(z)^{\ell_{k}}$ とおき、 $t_{v_{k}=}t_{v\chi_{k}(A)}$ をすべて求める。

2. $i=1,$について次を行う。
ならば、 $\ell_{kj}=0$ と置く。

$\cdots,$ $n$$t_{v_{k}e_{j}=0}$$0$

でなければ、$t_{v}$に順に $f_{k}(A)$ を右からかけて、 $tv_{k}f_{k}(A)^{p}$ を求め、 $t_{v_{k}f_{k}(A)^{p}b_{kj}=0}$ となる最小の $p>0$ で $\ell_{kj}=p$ と置く。

3. $j=1,$について次を行う。
推定した指数ベクトル ([最小消去多項式候補」 と呼ぶ) が本当に最小消去多項式になっているか、 実$\cdots,$ $n$

際に $\pi j(A)e_{j}$ を計算して確認する。 最小消去多項式でなければアルゴリズム 2 でやり直し。

26
ステップ 1 では、 すべての$v_{1},$$\cdots,$ $v_{m}$を求めるには、 ベキ乗算

$tu\mapsto tuf_{k}(A)^{\ell_{k}}$が$m^{2}$回必要であるように思えるが、 2 分探索の方法を用いれば、 およそ $m\log_{2}m$ 回で実行可能である。 これを詳しく説明しよう。 いま、 $S=\{1, \cdots, m\}$ の部分集合$I$に対して、 ベクトル

$v_{I}$を$t_{v_{I}=^{t}v\prod_{i\in S\backslash I}f_{i}(A)^{\ell_{i}}}$で定める。すると、 $v_{k}=v_{\{k\}}$ かつ $v=v_{S}$ である。 $S$ を根とし、 singleton $\{k\},$ $(k=1, \cdots, m)$ を葉とする二分木を次の規則で定める。
$I$が親で$J,$$K$がその子であるならば、 $I=J\sqcup K$ (disjoint union) である。

例 3. $m=4$ とする。 このとき、 集合$T=\{\{1,2,3,4\}, \{1,2\}, \{3,4\}, \{1\}, \{2\}, \{3\}, \{4\}\}$

は $S=\{1,2,3,4\}$ を根とし、 $\{k\}$ を葉とする二分木の構造が入る。

さて、$I$と $J$ が親子関係にあるとき、 ベクトルの親子$v_{I}$, v」の間にも関係式$t_{v_{J}=^{t}v_{I}\prod_{i\in I\backslash J}f_{i}(A)^{\ell}}.,$がある。 よって親ノードの計算結果を用いて子ノードのベクトルが計算できる。 二分木を辿ることによって、 すべての葉におけるベクトル$v_{k}$が求まる。

は一度計算しておけば、 どの基本ベクトルの計算にも利用できることにも注意する。 ま$tv\chi_{k}(A)$の計算の方がはるかに計算量が小さい。 よって、た行列を計算するより、 ベクトルを計算しておいて、 再利用を試みるのは計算量的に危険だが、 $tv\chi_{k}(A)$ ではそれほどでもない。
さらに、$v_{k}$$e_{j}$$\chi_{k}(A)$$\chi_{k}(A)$4.3

最小消去多項式導出の並列計算法

ここでは、 アルゴリズム 3 の並列化について議論したい。 ステップ 2,3 は単純なループであり、 ステップ

$2$は各 $j,$$k$について、 ステップ 3 は各 $i$ について計算が独立であるから、 並列化できることは明らかである。 よって問題になるのはステップ 1 である。 前節で説明したように、 ステップ 1 では二分木を辿りながら$v_{1},$$\cdots,v_{m}$を求めることができる。 このとき親子関係にないノードでは、 独立に計算を実行できるが、 親子関係にあるときには、 計算順序の制御が必要である。数式処理システム $Risa/Asir$ における並列計算の実装には、 oh-p パッケージが利用できるので、 これにもとついて二分木を辿る並列計算を説明しよう ([9])。
いま二分木を根の方から辿るような計算をすることとし、 子ノードの計算には親ノードの計算結果が必要であると仮定する。親ノードと子ノードが同時に計算されると、 不都合が生じるので、 同時に計算を行うのは、親子関係にないノードでなければならない。 また最初は根の計算の時には、 個のプロセスがあって$n$も 1 個だけ計算を行い、他のプロセスは待機することになる。とし、 既に計算が終わったジョブの全体、 計算中のジョいま、 二分木の各ノードの全体をジョブ集合とする。つまり、 $T=T_{f}\sqcup T_{b}u$ 窺である。 またブの全体、残りのジョブの全体をそれぞれ、とする。次のような関数
を考える。
$T$

$T_{f},T_{b},T_{u}$$\emptyset\not\in T$$\varphi$$\varphi(T_{u}, T_{f})=\{\begin{array}{l}x ( x\in T_{u} かつ parent (x)\in T_{f})\emptyset (その他)\end{array}$

27
このは iterator として働く。返り値がのときには、 ビジー状態であるとして oh-p.foreach r はが変化するまで待機する。 oh-p では oh-p foreach のオプションで iterator を引数として渡すことができる。 ま$\varphi$$\emptyset$$T_{f}$

た、 親ノードにおける計算結果を利用するために、 oh-p foreach の filter オプションを併用する。

4.4

最小消去多項式導出に関する実験結果

前節で述べた最小消去多項式導出の並列計算法を数式処理システム Risa $/Asir$ に実装し、 実験を行った。

..
.

実験に用いた環境は以下の通りである。

CPU: Xeon 2. $93GHz$ (4-core)$\cross 2$

Memory: $48GB$$OS$: Debian GNU Linux 2.6

また、 160 次整数行列 $A$ を、 以下の条件を満たすように、 コンパニオン行列のブロック対角行列に対し.て、 行列基本変形を繰り返して作成した。
$A$の特性多項式$\chi(z)$と最小多項式 $\pi(z)$ が、$\chi(z)=f_{1}(z)^{4}\cdots f_{8}(z)^{4}, \pi(z)=f_{1}(z)^{3}\cdots f_{8}(z)^{3}$となる既約分解をもつ。
$\circ$

各既約因子の次数は 5 である。

このとき、 160 次整数行列 $A$ を入カデータとして、 次のタイミングデータを得た。

表 1:160 次整数行列に対する最小消去多項式導出.
.

このことから以下のことが観察される。

計算時間は並列数にほぼ比例する。
並列数 8 でおよそ 4 倍の性能である。

参考文献[1] 小原功任，田島慎一: 行列のスペクトル分解・固有ベクトルの分散計算，京都大学数理解析研究所講究
録 1666(2008), 65-68.
[2] K. Ohara and S. Tajima: Spectral Decomposition and Eigenvectors of Matrices by Residue Calculus,Proceedings of the Joint Conference of ASCM 2009 and MACIS 2009, $COE$ Lecture Note 22 (2009),
Kyushu University, 137-140.

[3] 加藤涼香，田島慎一: 有理関数のローラン展開アルゴリズムと代数的局所コホモロジー，京都大学数理解析研究所講究録 1395 (2004), 50-56.

28[4] 庄司卓夢，田島慎一:1 変数代数的局所コホモロジー類に対する Risa$/Asir$ 用パッケージ taji-alc,$Risa/$ Asir Journa12 (2007), 1-32.

[5] 田島慎一，飯塚由貴恵: 行列のスペクトル分解アルゴリズムについて，京都大学数理解析研究所講究録
1666 (2008), 49-56.
[6] 飯塚由貴恵，田島慎一: 行列のスペクトル分解アルゴリズムー最小多項式が重複因子を持つ場合，数式処理 16-2 (2009).

[7] 飯塚由貴恵，田島慎一: 行列のスペクトル分解アルゴリズムー最小多項式が複数の重複因子から成る
場合，京都大学数理解析研究所講究録掲載予定．

[8] 田島慎一，奈良洗平: 最小消去多項式候補とその応用，京都大学数理解析研究所講究録投稿中．

[9] 小原功任: Open $XM$ を用いた $Risa/Asir$ 並列計算フレームワークの開発，数式処理 $1S$ -1 (2011).