【SymPy】線形代数演算
SymPy パッケージ の行列オブジェクト (MutableDenseMatrix) に対して実行できる 線形代数演算 (行列演算) について解説します。以下のコードは Jupyter Notebook で実行することを想定しています。他の IDE で実行される場合は適宜 print() 関数などを補ってください。
行列オブジェクトの生成
行列の各行をリストで sympy.Matrix() に渡すと、行列オブジェクトが生成されます。
# SYMPY_MATRIX_01
# In[1]
import sympy
# すべての小文字アルファベットを記号として定義
sympy.var("a:z")
# 3次正方行列を作成
M = sympy.Matrix([[a, b, c],
[d, e, f],
[g, h, i]])
M
# Matrix([
# [a, b, c],
# [d, e, f],
# [g, h, i]])
shape属性にアクセスすると行列の形状 (行数と列数) を取得できます。
# In[2] # 行列の形状 M.shape # (3, 3)
転置行列は T属性に記録されています。
# In[3] # 転置行列 M.T # Matrix([ # [a, d, g], # [b, e, h], # [c, f, i]])
行列の要素へのアクセス
行列の特定要素にアクセスするときは行, 列の順に指定します。
# In[4] # 2行3列の値を参照 M[1, 2] # f
row()メソッドで行にアクセスできます。
# In[5] # 第1行を参照 M.row(0) # Matrix([[a, b, c]])
col()メソッドで列にアクセスできます。
# In[6] # 第3列にアクセス M.col(2) # Matrix([ # [c], # [f], # [i]])
行と列の削除・挿入
M.row_del(), M.col_del() は、それぞれ指定した行、列を削除します。行列を変更するメソッドなので、扱うときには注意してください。
# SYMPY_MATRIX_02
# In[1]
import sympy
sympy.var("a:z")
# 3次正方行列を作成
M = sympy.Matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
# 2行目を削除
M.row_del(1)
# Mを表示
# display()はJupyterの関数です
# 他の環境ではprint()を使ってください
display(M)
# 1列目を削除
M.col_del(0)
# Mを表示
display(M)
# Matrix([
# [1, 2, 3],
# [7, 8, 9]])
# Matrix([
# [2, 3],
# [8, 9]])
M.row_insert(), M.col_insert() は、それぞれ指定した行、列に別の行列を挿入します。このメソッドは M を変更しません。
# In[2]
M1 = sympy.Matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
A = sympy.Matrix([[0, 0, 0]])
B = sympy.Matrix([[1],
[1],
[1],
[1]])
# Mの２行目にAを挿入
M2 = M1.row_insert(1, A)
# M2を表示
display(M2)
# M2の4列目にBを挿入
M3 = M2.col_insert(3, B)
# M3を表示
display(M3)
'''
Matrix([
[1, 2, 3],
[0, 0, 0],
[4, 5, 6],
[7, 8, 9]])
Matrix([
[1, 2, 3, 1],
[0, 0, 0, 1],
[4, 5, 6, 1],
[7, 8, 9, 1]])
'''
特殊行列の生成
sympy.eye() は単位行列を生成します。
# SYMPY_MATRIX_03 # In[1] import sympy # 単位行列を生成 sympy.eye(3) # Matrix([ # [1, 0, 0], # [0, 1, 0], # [0, 0, 1]])
sympy.zeros() は零行列を生成します。
# In[2] # 2行3列のゼロ行列を生成 sympy.zeros(2, 3) Matrix([ [0, 0, 0], [0, 0, 0]])
sympy.ones() は 1 で埋め尽くされた行列を生成します。
# In[3] # 全要素が1の行列を生成 sympy.ones(3) # Matrix([ # [1, 1, 1], # [1, 1, 1], # [1, 1, 1]])
sympy.diag() で対角行列を生成できます。
# In[4] # 対角行列を生成 sympy.diag(1, 2, 3) # Matrix([ # [1, 0, 0], # [0, 2, 0], # [0, 0, 3]])
行列オブジェクトのメソッド
2 次正方行列 A, B を生成しておきます。
# SYMPY_MATRIX_04
# In[1]
import sympy
sympy.var("a:z")
# 2次正方行列を生成
A = sympy.Matrix([[a, b],
[c, d]])
B = sympy.Matrix([[e, f],
[g, h]])
M.inv() は M の逆行列を返します。
# In[2] # Aの逆行列 A.inv() # Matrix([ # [ d/(a*d - b*c), -b/(a*d - b*c)], # [-c/(a*d - b*c), a/(a*d - b*c)]])
M.det() は M の行列式の値を計算します。
# In[3] # Aの行列式 A.det() # a*d - b*c
M.transpose() は M の転置行列を返します。
# In[4]
M = sympy.Matrix([[a, b, c],
[d, e, f],
[g, h, i]])
# Mの転置行列
M.transpose()
# Matrix([
# [a, d, g],
# [b, e, h],
# [c, f, i]])
M.adjugate() は M の余因子行列を返します。
# In[5] # Aの余因子行列 A.adjugate() # Matrix([ # [ d, -b], # [-c, a]])
M.trace() は M のトレースを返します。
# In[6] # Aのトレース A.trace() # a + d
M.rank() は M のランク (階数) を返します。
# In[7] # Aのランク A.rank() # 2
行列オブジェクトの演算
演算子を使って行列同士の和、差、積、累乗などを計算できます。
# SYMPY_MATRIX_05
# In[1]
import sympy
sympy.var("a:z")
# 2次正方行列を生成
A = sympy.Matrix([[a, b],
[c, d]])
B = sympy.Matrix([[e, f],
[g, h]])
# 行列の和
A + B
# Matrix([
# [a + e, b + f],
# [c + g, d + h]])
# In[2]
# 行列の差
A - B
# Matrix([
# [a - e, b - f],
# [c - g, d - h]])
# In[3]
# AとBの行列積
A * B
# Matrix([
# [a*e + b*g, a*f + b*h],
# [c*e + d*g, c*f + d*h]])
# In[4]
# Aのスカラー倍
k * A
# Matrix([
# [a*k, b*k],
# [c*k, d*k]])
# In[5]
# Aのべき乗
sympy.simplify(A ** 2)
# Matrix([
# [a**2 + b*c, b*(a + d)],
# [ c*(a + d), b*c + d**2]])
**(-1) は逆行列を取得します。
# In[6] # Aの逆行列 A**(-1) # Matrix([ # [ d/(a*d - b*c), -b/(a*d - b*c)], # [-c/(a*d - b*c), a/(a*d - b*c)]])
ベクトルの生成
sympy.Matrix() に k 個の要素をもつ 1 次元リストを渡すと、k 次元の縦ベクトルが生成されます。すなわち、ベクトルは k × 1 の行列オブジェクトとして扱われます (1 × k の横ベクトルではないことに注意してください)。
# SYMPY_MATRIX_06
# In[1]
import sympy
sympy.var("a:z")
# 3次元縦ベクトルを作成
# 3次元縦ベクトルを作成
V = sympy.Matrix([p, q, r])
W = sympy.Matrix([s, t, u])
display(V, W)
'''
Matrix([
[p],
[q],
[r]])
Matrix([
[s],
[t],
[u]])
'''
V.norm() で V の大きさ (長さ) を計算できます。
# In[2] V.norm() # sqrt(Abs(p)**2 + Abs(q)**2 + Abs(r)**2)
dot()メソッドでベクトルの内積を計算できます。
# In[3] # VとWの内積 V.dot(W) # p*s + q*t + r*u
行列とベクトルの積は * 演算子で実行できます。
# In[4]
# 3次正方行列を作成
M = sympy.Matrix([[a, b, c],
[d, e, f],
[g, h, i]])
# MによるVの線形変換
M*V
# Matrix([
# [a*p + b*q + c*r],
# [d*p + e*q + f*r],
# [g*p + h*q + i*r]])
一般表記行列生成関数
線形代数の学習の際には
\[\begin{bmatrix}a_{11} & a_{12} & a_{13}\\a_{21} & a_{22} & a_{23}\\a_{31} & a_{32} & a_{33}\end{bmatrix}\]
のように各成分に添字のついた行列を使いたいこともあると思います。行数と列数、成分表記に使う文字の種類を指定して、一般的な行列を生成する general_matrix() 関数のコードを載せておきます。
# In[1]
from sympy import Symbol, var, init_printing
from sympy.matrices import Matrix
init_printing()
# 添字付一般表記行列生成関数
def general_matrix(m, n, s):
Symbol(s)
elements = lambda i, j : var('{}{}{}'.format(s, i+1, j+1))
return Matrix(m, n, elements)
3 行 5 列、各成分が記号 $a$ で表される一般行列を生成してみます。
# In[2] # 3行5列の一般表記行列を生成 general_matrix(3, 5, "a")
ランダム係数連立一次方程式
SymPy の行列計算の応用として、ランダム係数をもつ連立一次方程式を生成する関数を定義してみます。
# SYMPY_MATRIX_07
# In[1]
# SymPyをインポート
import sympy
# 小文字アルファベットを記号として定義
sympy.var("a:z")
# 連立一次方程式生成関数
def rand_leq(s, min=-9, max=9):
A = sympy.randMatrix(len(s), min=min, max=max)
x = sympy.randMatrix(len(s), 1, min=min, max=max)
b = A*x
v = sympy.Matrix(s)
return A*v, b, x
第 1 引数 s には未知数として使用する文字 (小文字アルファベット) をリストで渡します。たとえば、[x, y, z] を渡せば、三元一次連立方程式が生成されます。キーワード引数 min, max には係数の最小値と最大値を渡します。sympy.randMatrix() は指定範囲内の乱数を要素にもつ行列を生成する関数です。
rand_leq()関数を使って三元一次連立方程式を生成してみます。
# In[2] # 三元一次連立方程式を生成 # 係数の最小値-9,最大値9 eq = rand_leq([x, y, z], min=-9, max=9) display(eq[0], eq[1], eq[2]) ''' Matrix([ [3*x + 8*y - 4*z], [8*x - 4*y - 5*z], [ -5*x - y - 5*z]]) Matrix([ [ 50], [-58], [ -8]]) Matrix([ [-2], [ 8], [ 2]]) '''
rand_leq() の戻り値は方程式の (左辺、右辺、解) のタプルです。いずれも一次元行列オブジェクト (ベクトル) として格納されます。
コメント
下記は誤植と思われますので、ご確認ください。
MATRIX_02 In[2] プログラムで、M2 = M.row_insert(1, A) → M2 = M1.row_insert(1, A)
実行結果は変わりませんが、M3 = M2.col_insert(4, B) → M3 = M2.col_insert(3, B)
MATRIX_04 In[4] プログラムで、最初にあらためて 3 次正方行列 M の生成が必要。
M = sympy.Matrix([[a, b, c],
[d, e, f],
[g, h, i]])
MATRIX_06 In[4] プログラムで、M*V1 → M*V
直しました。
ありがとうございます。m(_ _)m