« 前のひとこと | トップページ | 次のひとこと »

cocolog:93455907

Python の SymPy で行列に関する証明をする方法。 (JRF 9820)

JRF 2022年4月26日 (火)

SymPy は、行列の計算において、直接、行列変数みたいなものは扱えず、行列の各要素を指定したものしか使えないが、例えば…

JRF2022/4/260120

<pre>
from sympy import *
a11, a12, a21, a22 = symbols("a11 a12 a21 a22")
M = Matrix([[a11, a12], [a21, a22]])
print(simplify(Eq(M * (M ** -1), eye(2))))
</pre>

JRF2022/4/261459

…などとすれば、True が表示され、行列変数 M に関する式 M * (M ** -1) == I が証明できる。

任意の n × m 行列について証明を一つ作ることはできなさそうだが、n, m が与えられたとき、それに対応する式の証明を作るようなプログラムを書くことができる場合はあるだろう。

JRF2022/4/262583

ただ、上の例において、det(M) == 0 の場合を考慮していないように、「完璧な証明」が欲しい場合は、SymPy だけでは力不足だろう。

JRF2022/4/269451

……。

追記。

修正 「直接、行列変数みたいなものは扱えず、行列の各要素を指定したものしか使えないが」→「行列の各要素を指定したものを使うとすれば」。

調べ間違ってた。行列変数 MatrixSymbol というものは用意されている。いくつか行列変数を積み上げる BlockMatrix というものも用意されている。

《Matrix Expressions — SymPy 1.10.1 documentation》
https://docs.sympy.org/latest/modules/matrices/expressions.html

JRF2022/4/263656

……。

ところで、上の M は 4 行 4 列になった段階で、今の私のマシンだと M ** -1 の計算が時間がかかり過ぎるものとなって終らなかった。

<pre>
from sympy import *
N = 4
M = Matrix([[Symbol("a" + str(i) + str(j)) for i in range(1, N + 1)]
for j in range(1, N + 1)])
M ** -1
</pre>

…などとすれば終らない。

JRF2022/4/264754

上のようなことをするためには、

<pre>
N = 4
M = MatrixSymbol("M", N, N)
print(simplify(Eq(M * (M ** -1), Identity(N))))
</pre>

…とすればいい。(この場合は N = Symbol('N') でも良い。)

JRF2022/4/260114

ただし、(N = 4 として)

<pre>
print(simplify(Eq(M * (M ** -1), eye(N))))
</pre>

…は True にはならない。Identity(N) と eye(N) はそのままでは同じにならないかららしい。

JRF2022/4/269255

<pre>
print(simplify(Eq((M * (M ** -1)).as_explicit(), eye(N))))
</pre>

…または…

<pre>
print(simplify(M * (M ** -1)).equals(eye(N)))
</pre>

…とすれば良い。

JRF2022/4/267393

« 前のひとこと | トップページ | 次のひとこと »