Python (SymPy) による行列の計算

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使って行列の計算をしてみたいと思います。Pythonのごく基本的な使い方については以下の記事を参照してください：

pianofisica.hatenablog.com

行列の諸操作
数値的計算

行列の諸操作

行列を入力する

Python (SymPy) で行列を入力するには

import sympy
sympy.init_printing()
sympy.var('m11, m12, n11, n21')
sympy.var('a11, a12, a21, a22, b11, b12, b21, b22') 
M = sympy.Matrix([
[m11,m12]
])
N = sympy.Matrix([
[n11],
[n21]
])
A = sympy.Matrix([
[a11,a12],
[a21,a22]
])
B = sympy.Matrix([
[b11,b12],
[b21,b22]
])

などとすれば良いです。

ここで $M$ は $(1,2)$ -行列、 $N$ は $(2,1)$ -行列、 $A$ 、 $B$ は $(2,2)$ -行列としました：

$\displaystyle{M=\begin{bmatrix} m_{11} & m_{12} \end{bmatrix} \qquad N=\begin{bmatrix} n_{11} \\ n_{21} \end{bmatrix} \qquad A=\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \qquad B=\begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}}$

行列の要素を参照する

行列の要素を参照するには

A[0,1]   #  つまり行列 A の (k ,l) 成分は A[ k-1 , l-1 ] です（ラベルは０から数えます）

B[1,1]

などとします。

行列の演算（積、和、定数倍、ベキ乗）

行列の積の定義から、行列の積 $AB$ 、 $MA$ 、 $BN$ 、 $MBN$ などが許されます：

A*B    #  行列の積はアスタリスク（通常の積と同じ）です

M*A

B*N

M*B*N

などとして計算します。行列のベキ乗は

A**2

です。

同じ型の２つの行列の線形結合は

sympy.var('a, b')
a*A+b*B

となります。

逆行列をとる

逆行列は

A.inv()

または

A**(-1)

によって求められます。

連立１次方程式を解く

求めた逆行列を応用する問題として、連立１次方程式

$~ a_{11}x_1 + a_{12} x_2 + \dots + a_{1n} x_n=y_1 \\ ~ a_{21}x_1 + a_{22} x_2 + \dots + a_{2n} x_n=y_2 \\ ~~ \vdots \\ ~ a_{n1}x_1 + a_{n2} x_2 + \dots + a_{nn} x_n=y_n$

を考えましょう。これは行列を使って

$\displaystyle{~\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

と表すことができます。その解は、行列

$\displaystyle{~A_n=\begin{pmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \dots & \vdots \\ a_{n1} & a_{n2} & \dots & a_{nn} \end{pmatrix}}$

の逆行列 $A_n^{-1}$ （いま逆行列は存在するものと仮定します）を使って

$\displaystyle{~\begin{pmatrix} x_1 \\ x_2 \\ \dots \\ x_n \end{pmatrix}=A_n^{-1}\begin{pmatrix} y_1 \\ y_2 \\ \dots \\ y_n \end{pmatrix}}$

として求めることができます。

例えば３変数 $x$ 、 $y$ 、 $z$ の連立１次方程式

$~ ax+by+cz=j \\ ~ dx+ey+fz=k \\ ~ gx+hy+iz=l$

は行列

$\displaystyle{~A_3=\begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}}$

の逆行列をもとめれば良いわけです：

import sympy
sympy.init_printing()
sympy.var('a, b, c, d, e, f, g, h, i') 
A3 = sympy.Matrix([
[a, b, c],
[d, e, f],
[g, h, i]
])
sympy.cancel(A3**(-1))

により、逆行列が

$\displaystyle{~A_3^{-1}=\frac{1}{\det A_3}\begin{pmatrix} ei-fh & -(bi-ch) & bf-ce \\ -(𝑑𝑖−𝑓𝑔) & 𝑎𝑖−𝑐𝑔 & -(𝑎𝑓−𝑐𝑑) \\ 𝑑ℎ−𝑒𝑔 & -(𝑎ℎ−𝑏𝑔) & 𝑎𝑒−𝑏𝑑 \end{pmatrix}}$

と求まります。ここで

$\displaystyle{~\det A_3=a(ei-fh)-b(di-fg)+c(dh-eg)}$

です。最後はこの逆行列にベクトル

$\displaystyle{~Y=\begin{pmatrix} j \\ k \\ l \end{pmatrix}}$

を掛け合わせれば解が得られます：

$\displaystyle{~\begin{pmatrix} x \\ y \\ z \end{pmatrix} = A_3^{-1}Y= \frac{1}{\det A_3}\begin{pmatrix} j(ei-fh) -k(bi-ch) +l(bf-ce) \\ -j(𝑑𝑖−𝑓𝑔) +k( 𝑎𝑖−𝑐𝑔) -l(𝑎𝑓−𝑐𝑑) \\ j(𝑑ℎ−𝑒𝑔) -k(𝑎ℎ−𝑏𝑔) +l( 𝑎𝑒−𝑏𝑑) \end{pmatrix}}$

これをsolveを使って求めた解と見比べてみましょう：

eq1=sympy.Eq(a*x+b*y+c*z, j)
eq2=sympy.Eq(d*x+e*y+f*z, k)
eq3=sympy.Eq(g*x+h*y+i*z, l)
sympy.solve ([eq1, eq2, eq3], [x, y, z])

より

$\displaystyle{x={{j(ei-fh)-k(bi-ch)+l(bf-ce)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ y={{−𝑗(𝑑𝑖−𝑓𝑔)+𝑘(𝑎𝑖−𝑐𝑔)−𝑙(𝑎𝑓−𝑐𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ z={{𝑗(𝑑ℎ−𝑒𝑔)−𝑘(𝑎ℎ−𝑏𝑔)+𝑙(𝑎𝑒−𝑏𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)} }}$

確かに解になっていることがわかります。

行列式を求める

行列 $A$ の行列式を求めるには

A.det()

と入力します。

行列 $A$ の余因子行列 $\widetilde{A}$ は

A.adjugate()

です：実際

A.inv() == 1/(A.det())*A.adjugate()

より

$\displaystyle{\widetilde{A}=\left(\det A\right) A^{-1}}$

が成り立ちます。

単位行列

単位行列は

sympy.eye(2)

です。引数の数字は行列のサイズを指定します。

たとえば、逆行列になっていることを確認するときなどに

sympy.simplify(A.inv()*A) == sympy.eye(2)

のようにして使うことができます。

トレースをとる

行列 $A$ のトレースをとるには

A.trace()

とします。

転置をとる

行列 $A$ の転置行列

$\displaystyle{A^t=\begin{bmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \end{bmatrix}}$

は

A.transpose()

です。

随伴行列をとる

行列 $A$ の随伴行列

$\displaystyle{A^\dagger=\begin{bmatrix} \overline{a_{11}} & \overline{a_{21}} \\ \overline{a_{12}} & \overline{a_{22}} \end{bmatrix}}$

をとるには

A.adjoint()

とします。

行列を対角化する

行列の固有値を求める

A.eigenvals()

行列の固有ベクトルを求める

A.eigenvects()

対角行列を求める

A.diagonalize()

より詳しくは次の記事を読んでみてください。

pianofisica.hatenablog.com

数値的計算

数値計算ライブラリNumPyでの行列の扱いについては以下の記事にまとめました。あわせて読んでみてください。

pianofisica.hatenablog.com

キーワード：Python、SymPy、行列