SymPyの使い方11 ～行列の演算 - つれづれなる備忘録

　今回はSymPyを用いて行列演算の方法として逆行列、行列式、固有値、固有ベクトルなどを計算する方法について紹介する。

行列式(Determinant)を得るには.det()を用いる。

from sympy import *
a,b,c,d=symbols('a b c d')
A=Matrix([[a, b], [c, d]])
A.det()
>a*d - b*c

逆行列を求めるには**(-1)または.inv()を用いる。なお**とすると行列の累乗が計算できる。

A**(-1)
>Matrix([
[ d/(a*d - b*c), -b/(a*d - b*c)],
[-c/(a*d - b*c),  a/(a*d - b*c)]])

A.inv()
>Matrix([
[ d/(a*d - b*c), -b/(a*d - b*c)],
[-c/(a*d - b*c),  a/(a*d - b*c)]])

.rref()を用いると階段行列に変形する。連立方程式

$\displaystyle \begin{eqnarray} 5x - 2y &=& 5 \\ x + y &=& 8 \end{eqnarray}$

の左辺の係数と右辺の定数をまとめた拡大係数行列M

$\displaystyle {\bf{M}}=\left( \begin{array}{cc} 5 & -2 & 5\\ 1 & 1 & 8 \end{array} \right)$

に対して.rref()を適用すると、ガウス消去法によって一番右の列に解が得られる。

M = Matrix([[5,-2,5], [1,1,8]])
M.rref()
>(Matrix([
 [1, 0, 3],
 [0, 1, 5]]), (0, 1))

固有値は.eigenvals()を用いる。{固有値:重複度}という辞書形式で出力される。

A.eigenvals()
>{a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1,
 a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1}

固有ベクトルは.eigenvects()を用いる。

A.eigenvects()
>[(a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2, 1, [Matrix([
   [-b/(a/2 - d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2)],
   [                                                   1]])]),
 (a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2, 1, [Matrix([
   [-b/(a/2 - d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2)],
   [                                                   1]])])]

文字式の行列計算は結果が冗長になりやすいのでinit_printing()などで見やすくするとよいかもしれない。

行列の対角化は.diagonalize()を用いる。正則行列Pと対角行列Dが出力され、もとの行列Nに対して

${\bf{N}}={\bf{P}}{\bf{D}}{\bf{P}}^{-1}$

を満たす。

N = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])
P,D=N.diagonalize()
P
>Matrix([
[0, 1, 1,  0],
[1, 1, 1, -1],
[1, 1, 1,  0],
[1, 1, 0,  1]])
D
>Matrix([
[-2, 0, 0, 0],
[ 0, 3, 0, 0],
[ 0, 0, 5, 0],
[ 0, 0, 0, 5]])

P*D*P**-1
>Matrix([
[3, -2,  4, -2],
[5,  3, -3, -2],
[5, -2,  2, -2],
[5, -2, -3,  3]])