\[ \def\mcC{\mathcal C} \def\mcF{\mathcal F} \def\bsA{\boldsymbol A} \def\bsC{\boldsymbol C} \def\bsN{\boldsymbol N} \def\bsQ{\boldsymbol Q} \def\bsR{\boldsymbol R} \def\bsS{\boldsymbol S} \def\bsZ{\boldsymbol Z} \]

多項式行列の単因子標準形

2020/12/01
2020/12/06 リニューアルとプログラム追加

単因子標準形

以下に多項式行列の単因子標準形に関する諸概念と諸定理を簡単に紹介しておく。

多項式行列

多項式行列とは、多項式を要素とする行列のことである。
体 $K$ の要素を係数とする $x$ の多項式を、習慣に従って $K[x]$ と表すことにする。
$K[x]$ の要素から成る行列のことを、簡単に $K[x]$ 上の行列と言う。

多項式行列の和や積はまた多項式行列である。任意の $K$ の要素を多項式行列に掛けても多項式行列となる。 $K$ 上の行列も多項式行列である。

$E$ を(適当なサイズの)単位行列とする。
$K$ 上の行列を $M$ とすると $xE-M$ は多項式行列となる。これを $M$ の特性 $x$-行列(あるいは単に特性行列)と言う。

行列の基本変形

多項式行列の行基本変形とは
(a) 行 $i$ と行 $j$ を交換する
(b) 行 $i$ に $K$ の元 $c \;(\not=0)$ を掛ける
(c) 行 $i$ の $f \in K[x]$ 倍を行 $j$ に加える
の3つの操作を言う。同様に列基本変形が定義される。
両者を合わせて多項式行列の基本変形と言う。

基本変形は可逆な操作である。
基本変形は行列で表現できる。

2つの多項式行列 $A$ と $B$ が基本変形の繰り返しによって互いに移るとき、$A$ と $B$ は対等であると言い $A\sim B$ で表す。
この場合 $B=SAT$ で表現できる。
$S$ は単位行列 $E$ に行基本変形を繰り返して生成される行列である。
$T$ は単位行列 $E$ に列基本変形を繰り返して生成される行列である。

単因子定理

多項式行列 $M$ に行の基本変形と列の基本変形を繰り返して次の形の行列 $\hat M$ に持っていける。
(a) $\hat M$ は対角行列である。
(b) $\hat M$ の $0$ でない対角要素は $x$ の多項式であるが、これらの最高次数の係数は $1$ である¹。
(c) $\hat M$ の対角要素を $e_1,e_2,...$ とするとき $e_{k+1}$ は $e_k$ で割り切れる。(ただし $0$ は $0$ で割り切れると考える²)
これを単因子定理と言う。また $\hat M$ を $M$ の単因子標準形と言う。

1. $x$ を含まない多項式の最高次数の係数とは定数項である。
2. つまり対角要素は $e_1,e_2,...,e_l,0,0,...$ の形になるのである。$0,0,...$ の部分は発生するとは限らない。

最小多項式

$M$ を正方行列とする。$M$ に対して或る多項式 $f(x)$ が存在して $f(M)=0$ となる。
ただし $f(x)$ の定数項には単位行列 $E$ を添えておく。
$f(x)$ の中で最小次数のものを $M$ の最小多項式と言う。
$\hat M$ を $xE - M$ の単因子標準形とすると $\hat M$ の対角要素のうち、最大次数の多項式は $M$ の最小多項式となっている。

$xE-M$ の行列式 $D(x)$ は $M$ の特性多項式と呼ばれる。最小多項式は特性多項式を割り切る。従って $D(M)=0$ となる。(ケイリー・ハミルトンの定理)

単因子標準形から Jordan 行列

或る行列 $M$ の特性行列の単因子標準形の対角要素は(逆順に書いて)
$$(x-3)^3 (x-2) ,\quad x-2 ,\quad x-2 ,\quad 1 ,\quad 1 ,\quad 1$$
であったとする。その場合、行列 $M$ の Jordan 標準形 $J$ は次のように求まる。

単因子標準形の対角要素に現れる $(x-\alpha)^\nu$ は1つの Jordan 細胞 $J_\nu(\alpha)$ に対応する。
この場合、Jordan 細胞 $J_1(2)=\begin{pmatrix}2\end{pmatrix}$ が $3$ つ、$J_3(3)=\begin{pmatrix}3&1&0\\0&3&1\\0&0&3\end{pmatrix}$ が $1$ つ存在するので
$$J=\begin{pmatrix}2&0&0&0&0&0 \\ 0&2&0&0&0&0 \\
0&0&2&0&0&0 \\ 0&0&0&3&1&0 \\
0&0&0&0&3&1 \\ 0&0&0&0&0&3 \\ \end{pmatrix}
$$
となる。

Jordan 行列との対応を求める場合、単因子標準形のすべての対角要素は一次式の積に分解されている必要がある。$K$ として複素数体 $\bsC$ が設定されていれば(理論的には)可能であるが、コンピュータにおいては無理数の厳密な扱いは極めて困難である。行列の次元が大きい場合には Jordan 行列は求められないと割り切った方が賢明であろう。

補注: マリツェフ(1) p.131 には「拡張された Jordan 細胞」の定義と説明がある。これは単因子標準形の対角要素が一次式の積に分解されていない場合に使える。分解されている場合には通常の Jordan 細胞に還元する。

行列多項式

行列を係数とする多項式を行列多項式と言う。
多項式行列は行列多項式で表現できる。逆も然り。例えば
$$\begin{pmatrix}x^2-2\,x-1&2\,x+3 \\ x+2&2\,x^2-3\,x \\
\end{pmatrix} = \begin{pmatrix}
1&0 \\
0&2
\end{pmatrix}
x^2 + \begin{pmatrix}
-2&2 \\
1&-3
\end{pmatrix}
x + \begin{pmatrix}
-1&3 \\
2&0
\end{pmatrix}$$
である。
行列多項式において $x$ の最高次数の行列が正則であるとき、正則な行列多項式と言う。
多項式行列 $M$ が逆元を持つとき、$M$ を可逆な多項式行列と言う。
単位行列に行基本変形を繰り返して得られる行列は可逆な多項式行列である。
単位行列に列基本変形を繰り返して得られる行列は可逆な多項式行列である。

行列多項式の除算:
$P,D$ を各々行列多項式とする。
$D$ が正則な行列多項式であれば、$P$ を $D$ で割った商 $Q$ と剰余 $R$ が定義できる。
$P=QD+R$ となる $Q,R$ を各々右整商、右剰余と言う。
$P=DQ+R$ となる $Q,R$ を各々左整商、左剰余と言う。

注意: 数体 $K$ 上の行列では「正則」と「可逆」は同じ意味である。しかるに多項式行列の場合に、「正則」を避けたのは、行列多項式の「正則」と混乱するからである。多項式行列と行列多項式は表している実態は同じであるに関わらず、「可逆」と「正則」の性質はむしろ相反している。
「可逆」は斎藤が多項式行列に対して使っている。「正則」はマリツェフ(1)が行列多項式に対して使っている。
多項式行列と行列多項式に関してはマリツェフ(1)が詳しい。

相似な行列の変換行列

$K$ 上の行列 $M_1$ と $M_2$ が相似であるとは、$P^{-1}M_1 P= M_2$ となる $K$ 上の正則行列 $P$ が存在することである。この $P$ の求め方を示す。

行列 $M_1$ と $M_2$ が相似であることの必要十分条件は、それらの特性行列(すなわち $xE-M_1$ と $xE-M_2$)の単因子標準形が一致していることである。これは単因子標準形に関する非常に重要な定理であるが、その証明は難しい。マリツェフ『線型代数学(1)』 p.132 にその証明と $P$ の求め方がある。あるいは斎藤『線形代数入門』 p.181 にもある。

従って $M_1$ と $M_2$ が相似であれば $U(x)M_1 V(x)=S(s)M_2 T(x)$ となる $K[x]$ の可逆行列 $U(x),V(x),S(x),T(x)$ が存在する。そこで $T(x)V(x)^{-1}$ を $xE-M_1$ で割った右剰余を $R$ とする。$R$ には $x$ は含まれていない。そして$R^{-1}M_2 R = M_1$ の関係が成立している。すなわち、この $R$ が変換行列 $P$ である。

従って $xE-M$ の単因子標準形の対角要素が $x$ の一次式の積にまで分解されていれば、$M$ の Jordan 標準形 $J$ が判り、$P^{-1}MP=J$ となる $P$ を得る。

Maxima のプログラム

有理数体を $\bsQ$ で表す。
ここでは $\bsQ[x]$ 上の行列の単因子標準形を求める Maxima のプログラムを紹介する。
整数行列の単因子標準形を求めるプログラムは http://ar.nyx.link/ide/ で既に公開済みである。多項式行列の場合は基本的な考え方は共通しているもののコーディングは遥かに難しい。

Maxima には膨大な関数が準備されており、その中には数値行列 $M$ の Jordan 標準形 $J$ を求める関数も含まれている。従って特殊な多項式行列である、行列 $M$ の特性行列 $xE - M$ ($E$ は単位行列) の単因子標準形は Jordan 標準形が得られれば直ちに得られるが、一般的な多項式行列の単因子標準形を求める関数は準備されていないようである。ここではそれを求めるプログラムを公開する。

補注: Jordan 標準形は行列の次元が小さいならば得られるが、少し大きくなると一般的には得られる保証はない。
Jordan 標準形から単因子標準形を求める考え方は筋が悪すぎる!

ダウンロード http:elmdiv.mac

$M$ を $x$ の多項式行列、$\hat M$ をその単因子標準形とすると、
\[ SMT=\hat M \] となる可逆な多項式行列 $S,T$ が存在する。このプログラム elmdiv.mac では $\hat M$ と共に $S,T$ をも求めている。なお $S,T$ は一意には定まらないことに注意しておく。

$\hat M,S,T$ を求めるアルゴリズムの基本的な考え方は整数行列の単因子標準形を求める場合と何ら変わらない。従ってここでは改めて解説しない。整数行列の場合には http://ar.nyx.link/ide/ で解説されている。

実行例

Maxima の会話的実行は既知とする。
Maxima を立ち上げ

directory("");

で示されるディレクトリに "elmdiv.mac" を置き

batchload("elmdiv.mac");

を実行する。$M$ が
\[ M=\begin{pmatrix} 1-x & x^2 & x \\ x & x & -x \\ 1+x^2 & x^2 & -x^2 \end{pmatrix} \] の場合(この $M$ はマリツェフ(1) p.114 から借用した)、

(%i2) M:matrix([1-x,x^2,x],[x,x,-x],[1+x^2,x^2,-x^2]);

(%o2) matrix([1-x,x^2,x],[x,x,-x],[x^2+1,x^2,-x^2])
(%i3) elmdiv(M);

step0 0
              [matrix([1,0,0],[0,1,0],[0,0,1]),
               matrix([1-x,x^2,x],[x,x,-x],[x^2+1,x^2,-x^2]),
               matrix([1,0,0],[0,1,0],[0,0,1])]
step1 1
     [matrix([-1,0,0],[0,1,0],[0,0,1]),
      matrix([x-1,-x^2,-x],[x,x,-x],[x^2+1,x^2,-x^2]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
step2 1
     [matrix([-1,0,0],[1,1,0],[x+1,0,1]),
      matrix([x-1,-x^2,-x],[1,x^2+x,0],[2,x^3+2*x^2,x]),
      matrix([1,0,0],[0,1,0],[0,0,1])]
中略
step1 2
     [matrix([1,1,0],[x,x-1,0],[-1,(-x)-1,1]),
      matrix([1,0,0],[0,x,0],[0,0,(-x^2)-x]),
      matrix([1,0,(-x^2)-x],[0,0,1],[0,1,(-x^2)-x+1])]
step1 3
     [matrix([1,1,0],[x,x-1,0],[1,x+1,-1]),
      matrix([1,0,0],[0,x,0],[0,0,x^2+x]),
      matrix([1,0,(-x^2)-x],[0,0,1],[0,1,(-x^2)-x+1])]
(%o3) [matrix([1,1,0],[x,x-1,0],[1,x+1,-1]),
       matrix([1,0,0],[0,x,0],[0,0,x*(x+1)]),
       matrix([1,0,-x*(x+1)],[0,0,1],[0,1,-(x^2+x-1)])]
(%i4)

となる。計算の途中結果も出力されている。
この結果から
\[ S=\begin{pmatrix} 1 & 1 & 0\\ x & x-1 & 0\\ 1 & x+1 & -1 \end{pmatrix} ,\quad \hat M=\begin{pmatrix} 1 & 0 & 0\\ 0 & x & 0\\ 0 & 0 & x(x+1) \end{pmatrix} ,\quad T=\begin{pmatrix} 1 & 0 & -x(x+1)\\ 0 & 0 & 1 \\ 0 & 1 & -(x^2+x-1) \end{pmatrix} \] を得る。

プログラム

elmdiv.mac の主関数だけを紹介すると

elmdiv(M):=block([B,n,k,x,r],
  n:min(length(M),length(M[1])),
  B:tab_gen(M),
  dprint("step0",0,B),
  /* don't use 'for' loop */
  k:1,
  /* don't replace B[2] with A, which behaves differently! */
  while k <= n and chk_zero(B[2],k)=false do block(/* it seems the 'block' is required. why? */
    X:step1(B,k),
    if X=false then (k:n+1, return(false)),
    B:X,
    dprint("step1",k,B),
  	loop, /* we should avoid go(loop) if possible. this is merely an example for go(loop) */
    /* don't replace B[2] with A, which behaves differently! */
    unless chk_col(B[2],k) do (
      B:step2(B,k),
      dprint("step2",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    unless chk_row(B[2],k) do (
      B:step3(B,k),
      dprint("step3",k,B),
      B:step1(B,k),
      dprint("step1",k,B)
    ),
    if chk_col(B[2],k)=false then go(loop),
    /* we may assume here both chk_col(B[2],k) and chk_row(B[2],k) are true */
    B:tab_fix(B,k),
    if chk_mod(B[2],k) then k:k+1
    else (
      X:step4(B,k),
      if X = false then k:n+1
      else (
        B:X,
        dprint("step4",k,B)
      )
    )
  ),
  r:chk_final(B,M),
  if r # 0 then print("# ERROR: bug!",r,B,M),
  return(factor(B))
);

この部分は lde.mac の中の elmdiv(M)と殆ど違わない。基本的なアルゴリズムは同じだからである。

最小多項式

$M$ を数値行列とすると、$xE-M$ の単因子標準形
$$\begin{pmatrix}1&0 \\ 0&x^2-9\,x-1 \\ \end{pmatrix}$$
から直ちに $M$ の最小多項式が求まる。最小多項式は単因子標準形の $x$ の最大次数の対角要素 $\hat d(x)$ なのである。elmdiv.mac には、これを取り出す関数 minpoly(M) が付属している。

M:matrix([2,3],[5,7]);
minpoly(M);

この出力は

x^2-9*x-1

で、これから最小多項式 $ x^2-9x -1 $ を得る。
これを Maxima で確認するには

M^^2 - 9*M - 1*ident(2)

と行列のベキ乗に変換する必要がある。その下で

matrix([0,0],[0,0])

を得る。
なお、最小多項式は Maxima でもサポートされている。ただし導き方が違う。こちらの方は Jordan 標準形を基に導く:

load("diag");
M:matrix([2,3],[5,7]);
minimalPoly(jordan(M));

この結果は

(x+(sqrt(85)-9)/2)*(x-(sqrt(85)+9)/2)

であり。これを

expand((x+(sqrt(85)-9)/2)*(x-(sqrt(85)+9)/2));

として展開してやると同じ結果が得られる。

行列多項式

多項式行列 M から

mp_gen(M);

で M の行列多項式が得られる。
実行例:

M:matrix([x^2-2*x-1,2*x+3],[x+2,2*x^2-3*x]);
mp_gen(M);

で、出力

[matrix([1,0],[0,2]),matrix([-2,2],[1,-3]),matrix([-1,3],[2,0])]

を得る。これは
\[ \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix} x^2 + \begin{pmatrix} -2 & 2 \\ 1 & -3 \end{pmatrix} x + \begin{pmatrix} -1 & 3 \\ 2 & 0 \end{pmatrix} \] を意味する。
$x$ の係数しか表示されない。係数のリストが行列多項式と同一視されている。

逆に $x$ の係数のリスト Mp から多項式行列を求めるには

mp_matrix(Mp);

を実行すればよい。

行列多項式の除算

行列多項式 P を行列多項式 D で右から割った商と剰余を求めるには

L:mp_rdivmod(P,D)
L[1]; /* 右整商 */
L[2]; /* 右剰余 */

である。
同様に左から割った商と剰余を求めるには

L:mp_ldivmod(P,D)
L[1]; /* 左整商 */
L[2]; /* 左剰余 */

である。

例:

\[
P=\begin{pmatrix}2&3 \\ -1&2 \\ \end{pmatrix} x^4 +
\begin{pmatrix}3&1 \\ 4&5 \\ \end{pmatrix} x^3 +
\begin{pmatrix}5&2 \\ 1&3 \\ \end{pmatrix} x^2 +
\begin{pmatrix}2&1 \\ 1&7 \\ \end{pmatrix} x +
\begin{pmatrix}5&1 \\ 2&3 \\ \end{pmatrix}\\
D=\begin{pmatrix}2&5 \\ 1&3 \\ \end{pmatrix} x^2 +
\begin{pmatrix}7&2 \\ 3&5 \\ \end{pmatrix} x +
\begin{pmatrix}-2&1 \\ 3&2 \\ \end{pmatrix}
\]
の下で、$ P=QD+R$ となる $Q,R$ を求めるには

P:[matrix([2,3],[-1,2]),matrix([3,1],[4,5]),matrix([5,2],[1,3]),matrix([2,1],[1,7]),matrix([5,1],[2,3])];
D:[matrix([2,5],[1,3]),matrix([7,2],[3,5]),matrix([-2,1],[3,2])];
L:mp_rdivmod(P,D);

を実行すればよい。結果は

(%i112) Q:L[1];
(%o112) [matrix([3,-4],[-5,9]),matrix([-33,60],[66,-120]),
         matrix([449,-824],[-872,1606])]
(%i113) R:L[2];
(%o113) [matrix([-915,3136],[1779,-6105]),matrix([3375,1200],[-6560,-2337])]

となり、
$$ Q=\begin{pmatrix}3&-4 \\ -5&9 \\ \end{pmatrix} x^2 +
\begin{pmatrix}-33&60 \\ 66&-120 \\ \end{pmatrix} x +
\begin{pmatrix}449&-824 \\ -872&1606 \\ \end{pmatrix}$$
と
$$ R=\begin{pmatrix}-915&3136 \\ 1779&-6105 \\ \end{pmatrix} x +
\begin{pmatrix}3375&1200 \\ -6560&-2337 \\ \end{pmatrix} $$
を得る。

相似な行列の変換行列

行列 $M_1$ と $M_2$ が相似、すなわち $P^{-1}M_1 P= M_2$ となる正則行列 $P$ が存在する場合に、$P$ を求めるプログラム mat_trans(M1,M2) が elmdiv.mac に含まれている。
例えば
$$M_1=\begin{pmatrix}0&-1&-1&0 \\ -1&1&0&1 \\ 2&1&2&-1
\\ -1&-1&-1&1 \\ \end{pmatrix},\qquad
M_2=\begin{pmatrix}1&1&0&0 \\ 0&1&0&0 \\ 0&0&1&1 \\ 0
&0&0&1 \\ \end{pmatrix}$$
とする。

M1:matrix([0,-1,-1,0],[-1,1,0,1],[2,1,2,-1],[-1,-1,-1,1]);
M2:matrix([1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1]);
mat_trans(M1,M2);

より

matrix([1,0,0,0],[-1,-1,-1,0],[0,0,1,0],[1,-1,0,-1])

を得る。すなわち
$$P=\begin{pmatrix}1&0&0&0 \\ -1&-1&-1&0 \\ 0&0&1&0
\\ 1&-1&0&-1 \\ \end{pmatrix}$$
である。
$M_1$ と $M_2$ が相似ではない場合には mat_trans(M1,M2) は false を返す。

Maxima の扱い難いところ

名前の衝突

気がついた点だけを書く。
名前の衝突を避けるのが非常に難しい。
Python の import 文に相当する処理ができない(らしい)。
この問題は深刻で、プログラムを組む側で丁寧に確認しなくてはならない。
それでも完全ではないであろう。今は筆者の環境では問題はなくても、環境が異なれば保証の限りでは無い。
例えば lde.mac と elmdiv.mac は同名の関数を共有している。従ってロードする Maxima を分けざるを得ない。これをやりやすくするために筆者は maxima コマンドを準備している。

Maxima のプログラムが特定の集団だけから配布され、利用者はそれを利用するだけであれば、それでも良いのであるが、利用者がプログラム作成に参加するのであれば、現在の仕様は破綻するであろう。

お節介

Maxima で

matrix([2]).matrix([3]);

を実行すると 6 が表示される。matrix([6]) ではない。
従って $1\times 1$ 行列の演算ではトラブルの原因になる。
このように、結果が $1\times 1$ 行列になったときに、結果を行列要素に置き換える気持ちはどこにあるかと言えば、内積の処理のために必要以上のことを行なっているのである。

もっと酷いのは

matrix([1,2]).matrix([3,4]);

がエラーにならずに

を出すのである。「きっと間違えたのでしょうね...　あんたの欲しいのはこれでしょ」と、忖度して返事するような設計は大嫌いである。怪我の元!

環境改善

プログラミングは、修正と実行の繰り返しである。
それがやりやすいように環境を改善しておく。

elmdiv.mac には

display2d:false;
ap:apropos; /* alias for apropos */
ld():=batchload("elmdiv.mac");

が既に含まれている。

Mac で maxima コマンド

Mac の Terninal.app からコマンドで Maxima を立ち上げられたら便利である。
ついでにヒストリ機能と入力行の編集機能を持たせたらさらに便利である。

#!/bin/sh
exec rlwrap /Applications/Maxima.app/Contents/Resources/maxima.sh "$1" "$2" "$3" "$3" "$4" "$5" "$6" "$7" "$8" "$9"

実行ファイル maxima

rlwrap をインストールしておくと良い。入力行の編集ができる。

SymPy のプログラム

有理数体を $\bsQ$ で表す。
ここでは $\bsQ[x]$ 上の行列の単因子標準形を求める SymPy のプログラムを紹介する。
整数行列の単因子標準形を求めるプログラムは http://ar.nyx.link/ide/ で既に公開済みである。多項式行列の場合は基本的な考え方は共通しているもののコーディングは遥かに難しい。

SymPy には膨大な関数が準備されており、その中には数値行列 $M$ の Jordan 標準形 $J$ を求める関数も含まれている。従って特殊な多項式行列である、行列 $M$ の特性行列 $xE - M$ ($E$ は単位行列) の単因子標準形は Jordan 標準形から直ちに得られるが、一般的な多項式行列の単因子標準形を求める関数は準備されていないようである。ここではそれを求めるプログラムを公開する。

ダウンロード http:melmdiv.py

$M$ を $x$ の多項式行列、$\hat M$ をその単因子標準形とすると、
\[ SMT=\hat M \] となる可逆な多項式行列 $S,T$ が存在する。このプログラム melmdiv.py では $\hat M$ と共に $S,T$ をも求めている。なお $S,T$ は一意には定まらないことに注意しておく。

実行例

SymPy の会話的実行は既知とする。
SymPy を立ち上げ

import os
os.getcwd()

で示されるディレクトリに "melmdiv.py" を置き

import melmdiv
from melmdiv import *

を実行する。$M$ が
\[ M=\begin{pmatrix} 1-x & x^2 & x \\ x & x & -x \\ 1+x^2 & x^2 & -x^2 \end{pmatrix} \] の場合(この $M$ はマリツェフ(1) p.114 から借用した)、

>>> M=Mat([1-x,x**2,x],[x,x,-x],[1+x**2,x*2,-x**2])
>>> elmdiv(M)
step0 0 [Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]), Matrix([
[   1 - x, x**2,     x],
[       x,    x,    -x],
[x**2 + 1,  2*x, -x**2]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
step1 1 [Matrix([
[-1, 0, 0],
[ 0, 1, 0],
[ 0, 0, 1]]), Matrix([
[   x - 1, -x**2,    -x],
[       x,     x,    -x],
[x**2 + 1,   2*x, -x**2]]), Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])]
中略
step1 3 [Matrix([
[  1,         1,    0],
[  x,     x - 1,    0],
[1/2, x/2 + 1/2, -1/2]]), Matrix([
[1, 0,          0],
[0, x,          0],
[0, 0, x**2 - x/2]]), Matrix([
[1, 0,     -x**2 - x],
[0, 0,             1],
[0, 1, -x**2 - x + 1]])]
[Matrix([
[  1,         1,    0],
[  x,     x - 1,    0],
[1/2, (x + 1)/2, -1/2]]), Matrix([
[1, 0,             0],
[0, x,             0],
[0, 0, x*(2*x - 1)/2]]), Matrix([
[1, 0,    -x*(x + 1)],
[0, 0,             1],
[0, 1, -x**2 - x + 1]])]
>>>

プログラム

melmdiv.py の主関数だけを紹介すると

def elmdiv(M):
  n = min(M.shape)
  B = tab_gen(M)
  dprint("step0",-1,B)
  k = 0
  # don't replace B[1] with A, which behaves differently!
  while k < n and chk_zero(B[1],k)==False:
    X = step1(B,k)
    if X == False: break
    B = X
    dprint("step1",k,B)

    while True:
      while chk_col(B[1],k) == False:
        B = step2(B,k)
        dprint("step2",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      while chk_row(B[1],k) == False:
        B = step3(B,k)
        dprint("step3",k,B)
        B = step1(B,k)
        dprint("step1",k,B)
      if chk_col(B[1],k) == True: break

    # we may assume here both chk_col(B[1],k) and chk_row(B[1],k) are true
    if chk_mod(B[1],k): k += 1
    else: # BUG! ned this one. you must also fix mlde.py. examine also lde.mac!
      X = step4(B,k)
      if X==False: break
      B = X
      dprint("step4",k,B)
  r = chk_final(B,M)
  if r != 0 : error("# ERROR: bug!",r)
  return mat_factor(B)

この部分は mlde.py の中の elmdiv(M)と殆ど違わない。基本的なアルゴリズムは同じだからである。

最小多項式

$M$ を数値行列とすると、$xE-M$ の単因子標準形
$$\begin{pmatrix}1&0 \\ 0&x^2-9\,x-1 \\ \end{pmatrix}$$
から直ちに $M$ の最小多項式が求まる。最小多項式は単因子標準形の $x$ の最大次数の対角要素 $\hat d(x)$ なのである。melmdiv.py には、これを取り出す関数 minpoly(M) が付属している。

M=Mat([2,3],[5,7])
minpoly(M)

この出力は

x**2-9*x-1

で、これから最小多項式 $ x^2-9x -1 $ を得る。
これを SymPy で確認するには

M**2 - 9*M - 1*eye(2)

を実行すればよい。その下で

Matrix([
[0, 0],
[0, 0]])

を得る。

行列多項式

多項式行列 M から

mp_gen(M)

で M の行列多項式が得られる。
実行例:

M=Mat([x**2-2*x-1,2*x+3],[x+2,2*x**2-3*x])
mp_gen(M)

で、出力

[Matrix([
[1, 0],
[0, 2]]), Matrix([
[-2,  2],
[ 1, -3]]), Matrix([
[-1, 3],
[ 2, 0]])]

逆に $x$ の係数のリスト Mp から多項式行列を求めるには

mp_matrix(Mp)

を実行すればよい。

行列多項式の除算

行列多項式 P を行列多項式 D で右から割った商と剰余を求めるには

L=mp_rdivmod(P,D)
L[0] # 右整商
L[1] # 右剰余

である。
同様に左から割った商と剰余を求めるには

L=mp_ldivmod(P,D)
L[0] # 左整商
L[1] # 左剰余

である。

例:

P=[Mat([2,3],[-1,2]),Mat([3,1],[4,5]),Mat([5,2],[1,3]),Mat([2,1],[1,7]),Mat([5,1],[2,3])]
D=[Mat([2,5],[1,3]),Mat([7,2],[3,5]),Mat([-2,1],[3,2])]
L=mp_rdivmod(P,D)

を実行すればよい。結果は

>>> L
([Matrix([
[ 3, -4],
[-5,  9]]), Matrix([
[-33,   60],
[ 66, -120]]), Matrix([
[ 449, -824],
[-872, 1606]])], [Matrix([
[-915,  3136],
[1779, -6105]]), Matrix([
[ 3375,  1200],
[-6560, -2337]])])

相似な行列の変換行列

行列 $M_1$ と $M_2$ が相似、すなわち $P^{-1}M_1 P= M_2$ となる正則行列 $P$ が存在する場合に、$P$ を求めるプログラム mat_trans(M1,M2) が melmdiv.py に含まれている。
例えば
$$M_1=\begin{pmatrix}0&-1&-1&0 \\ -1&1&0&1 \\ 2&1&2&-1
\\ -1&-1&-1&1 \\ \end{pmatrix},\qquad
M_2=\begin{pmatrix}1&1&0&0 \\ 0&1&0&0 \\ 0&0&1&1 \\ 0
&0&0&1 \\ \end{pmatrix}$$
とする。

M1=Mat([0,-1,-1,0],[-1,1,0,1],[2,1,2,-1],[-1,-1,-1,1])
M2=Mat([1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1])
mat_trans(M1,M2)

より

Matrix([
[ 1,  0,  0,  0],
[-1, -1, -1,  0],
[ 0,  0,  1,  0],
[ 1, -1,  0, -1]])

を得る。すなわち
$$P=\begin{pmatrix}1&0&0&0 \\ -1&-1&-1&0 \\ 0&0&1&0
\\ 1&-1&0&-1 \\ \end{pmatrix}$$
である。
$M_1$ と $M_2$ が相似ではない場合には mat_trans(M1,M2) は False を返す。

Maxima と SymPy との比較

Maxima も SymPy のプログラムも今回が初めてである。
適当な課題を定めてプログラムをしながら言語を少しづつマスターするのが筆者のやり方である。

根本的に違う世界の言語である。細かなことは違いが多すぎて表にはできない。基本的なことだけを取り上げる。

関数型言語と手続き型言語

Python は手続き型言語であるのに対して、 Maxima は手続き型言語のように見せてはいるが、本質的には関数型言語である。
このことは、Maxima では命令の区切りがあたかも "," であるように見せているが、命令のように見えても関数引数に値を与えているのである。そのために関数引数の区切り記号 "," が使われている。
Maxima の block はブロック関数とでも言うべきもので、lambda と同様に特殊な無名関数なのである。
Maxima の return(...) は常にブロック関数からの戻りを意味する。

Maxima では繰り返しを表現するときには、本質的には関数であることを頭に入れておかないと訳のわからない動作に悩まされる。

名前の参照規則

手続き型言語は名前の参照先はコーディング時に決まる。ところが Maxima は実行時に決まる。
局所変数でない変数は親に何段にも遡って参照される。関数が参照する名前空間が階層構造になっているのである。これは都合の良いことも多いが、戸惑うこともある。

Python のモジュールに相当するものは無いようである。これはとても困る。

言語の生い立ち

Maxima は数式処理に特化した言語であり、元になった Macsyma を含めれば言語仕様は1982年にまで遡る[1]。言語仕様はかなり癖が強い。行列の乗算記号が "*" ではなく "."なのはオブシェクト化ができない代償である。

他方 SymPy は汎用言語の Python のモジュールの一つに過ぎない。従って言語仕様は Python そのものである。Python に慣れている者にとっては親和性が高い。最初の版は2007年であるから[2]、まだ使えるようになってからの歴史はまた浅く、発展途上にあると考えた方が良い。将来性は SymPy に軍配が上がるのではないかと思える。

[1] https://ja.wikipedia.org/wiki/Maxima
[2] https://en.wikipedia.org/wiki/SymPy