Python 数値計算入門 Jupyter NoteBook を活用したPython3 プログラミング

[SymPy] 行列とベクトル

[SymPy] 行列とベクトル
目次

SymPy 線型代数演算

　SymPy の行列オブジェクト (MutableDenseMatrix) に対して実行できる線型代数演算 (行列演算) について解説します。以下のコードは Jupyter Notebook で実行することを想定しています。他の IDE で実行される場合は適宜 print()関数などを補ってください。

行列オブジェクトの生成

　行列の各行をリストで sympy.Matrix() に渡すと、行列オブジェクトが生成されます。

# リストMTX01-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属性にアクセスすると行列の形状 (行数と列数) を取得できます。

# リストMTX01-2

# 行列の形状
M.shape
(3, 3)

　転置行列は T属性に記録されています。

# リストMTX01-3

# 転置行列
M.T
Matrix([
[a, d, g],
[b, e, h],
[c, f, i]])

　

行列の要素へのアクセス

　行列の特定要素にアクセスするときは行, 列の順に指定します。

# リストMTX01-4

# 2行3列の値を参照
M[1, 2]
f

　row()メソッドで行にアクセスできます。

# リストMTX01-5

# 第1行を参照
M.row(0)
Matrix([[a, b, c]])

　col()メソッドで列にアクセスできます。

# リストMTX01-6

# 第3列にアクセス
M.col(2)
Matrix([
[c],
[f],
[i]])

　

行と列の削除・挿入

　M.row_del(), M.col_del() は、それぞれ指定した行、列を削除します。M を変更するメソッドなので、扱うときには注意してください。

# リストMTX02

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 を変更しません。

# リストMTX03

import sympy
sympy.var("a:z")

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 = M.row_insert(1, A)

# M2を表示
display(M2)

# M2の4列目にBを挿入
M3 = M2.col_insert(4, 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() は単位行列を生成します。

# リストMTX04-1

import sympy

# 単位行列を生成
sympy.eye(3)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])

　sympy.zeros() は零行列を生成します。

# リストMTX04-2

# 2行3列のゼロ行列を生成
sympy.zeros(2, 3)
Matrix([
[0, 0, 0],
[0, 0, 0]])

　sympy.ones() は 1 で埋め尽くされた行列を生成します。

# リストMTX04-3

# 全要素が1の行列を生成
sympy.ones(3)
Matrix([
[1, 1, 1],
[1, 1, 1],
[1, 1, 1]])

　sympy.diag() で対角行列を生成できます。

# リストMTX04-4

# 対角行列を生成
sympy.diag(1, 2, 3)
Matrix([
[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])

　

行列オブジェクトのメソッド

　2 次正方行列 A, B を生成しておきます。

# リストMTX05-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 の逆行列を返します。

# リストMTX05-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 の行列式の値を計算します。

# リストMTX05-3

# Aの行列式
A.det()
a*d - b*c

　M.transpose() は M の転置行列を返します。

# リストMTX05-4

# Mの転置行列
M.transpose()
Matrix([
[a, d, g],
[b, e, h],
[c, f, i]])

　M.adjugate() は M の余因子行列を返します。

# リストMTX05-5

# Aの余因子行列
A.adjugate()
Matrix([
[ d, -b],
[-c,  a]])

　M.trace() は M のトレースを返します。

# リストMTX05-6

# Aのトレース
A.trace()
a + d

　M.rank() は M のランク (階数) を返します。

# リストMTX05-7

# Aのランク
A.rank()
2

　

行列オブジェクトの演算

　演算子を使って行列同士の和、差、積、累乗などを計算できます。リストMTX05-1 で定義した行列 A, B を使うので、引き続きリストMTX05 シリーズのコードです。

# リストMTX05-6

# 行列の和
A + B
Matrix([
[a + e, b + f],
[c + g, d + h]])

　

# リストMTX05-7

# 行列の差
A - B
Matrix([
[a - e, b - f],
[c - g, d - h]])

　

# リストMTX05-8

# AとBの行列積
A * B
Matrix([
[a*e + b*g, a*f + b*h],
[c*e + d*g, c*f + d*h]])

　

# リストMTX05-9

# Aのスカラー倍
k * A
Matrix([
[a*k, b*k],
[c*k, d*k]])

　

# リストMTX05-10

# Aのべき乗
sympy.simplify(A ** 2)
Matrix([
[a**2 + b*c,  b*(a + d)],
[ c*(a + d), b*c + d**2]])

　
　**(-1) は逆行列を取得します。

# リストMTX05-11

# 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 の横ベクトルではないことに注意してください)。

# リストMTX06-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 の大きさ (長さ) を計算できます。

# リストMTX06-2

V.norm()
sqrt(Abs(p)**2 + Abs(q)**2 + Abs(r)**2)

　dot()メソッドでベクトルの内積を計算できます。

# リストMTX06-2

# VとWの内積
V.dot(W)
p*s + q*t + r*u

　行列とベクトルの積は * 演算子で実行できます。

# リストMTX06-3

# 3次正方行列を作成
M = sympy.Matrix([[a, b, c],
                  [d, e, f],
                  [g, h, i]])

# MによるVの線型変換
M*V1
Matrix([
[a*p + b*q + c*r],
[d*p + e*q + f*r],
[g*p + h*q + i*r]])

　

ランダム係数連立一次方程式

　SymPy の行列計算の応用として、ランダム係数をもつ連立一次方程式を生成する関数を定義してみます。

# リストMTX07-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()関数を使って三元一次連立方程式を生成してみます。

# リストMTX07-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() の戻り値は方程式の (左辺、右辺、解) のタプルです。
　いずれも一次元行列オブジェクト (ベクトル) として格納されます。