読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

MATLAB, EigenユーザのためのPython行列演算入門

Python

プログラミングのための線形代数

プログラミングのための線形代数

目次

 

Pythonにおける行列演算

先日紹介したmatplotlibという

グラフ描画ツールを上手く使うと、

PythonをMATLABのように

アルゴリズムのプロトタイプ作成

に使いやすくなります。

myenigma.hatenablog.com

 

また、PythonもMATLABに比べると計算演算も早いので、

そこまで計算コストが大きく無い場合、

C++でアルゴリズムを実装する代わりに

Pythonを使うことができます。

 

様々なアルゴリズムを実装する際に、

重要なのは行列演算機能です。

myenigma.hatenablog.com

 

行列演算の四則演算や、逆行列・固有値の計算を

通常のプログラミング言語で実装するは骨が折れます。

しかし、

MATLABには行列演算機能がデフォルトでbuild-inされていますし、

C++にはEigenという非常に使いやすいライブラリがあります。

myenigma.hatenablog.com

 

実はPythonにおいても、

強力な行列演算ツールがあるため、

それらを使いこなせば、

PythonをMATLABやEigen代わりに使えると思います。

 

今回の記事では、MATLABユーザ向けに

Pythonにおける行列演算ツールと使い方について

概要を説明したいと思います。

 

NumpyとScipy

Pythonで行列演算をする場合、

NumpyとScipyというライブラリを使うと便利です。

 

NumpyはNumeric Pythonの略で数値演算用のPythonライブラリです。

Numpyのメイン機能は、多次元配列(行列)の計算です。

これらのメソッドは、内部でC言語で実装しているため、

大規模な行列演算を実施しても、比較的高速というメリットがあります。

 

ScipyはScientific Pythonの略で科学計算用のライブラリです。

基本的にはNumpyをベースとして開発されており、

Numpyの高速演算メソッドを使って、

統計、最適化、積分、線形代数、フーリエ変換、

信号・イメージ処理、遺伝的アルゴリズム、常微分方程式ソルバなどの

科学計算用の機能が実装されています。

 

Pythonで行列演算を実施したい場合は、

基本的にnumpyを使うと十分かと思います。

 

Pythonで行列演算のサンプルプログラム

下記がNumpyの行列演算のサンプルプログラムです。

 

行列やベクトルの初期化

下記のように、

様々な形の代表的なベクトルや行列は

非常にシンプルな形で初期化して生成することができます。

 

print "====initialize===="
x=np.array([[1],[2],[3]])
print x
#  [[1]
#  [2]
#  [3]]

# ベクトルの繰り返し
print x.repeat(3)
#  [1 1 1 2 2 2 3 3 3]

# ベクトルの繰り返し(繰り返し方向の指定)
print x.repeat(3,axis=1)
#  [[1 1 1]
#  [2 2 2]
#  [3 3 3]]

# ベクトルそのものの繰り返しベクトル
print np.tile([0,1,2,3,4], 2)
#  [0 1 2 3 4 0 1 2 3 4]

#格子行列
a,b = np.meshgrid([1,2,3], [4,5,6,7])
print a
#  [[1 2 3]
#  [1 2 3]
#  [1 2 3]
#  [1 2 3]]
print b
#  [[4 4 4]
#  [5 5 5]
#  [6 6 6]
#  [7 7 7]]

#ゼロベクトル
print np.zeros(5)
#  [ 0.  0.  0.  0.  0.]

#ゼロ行列
print np.zeros([2,2])
#  [[ 0.  0.]
#  [ 0.  0.]]

#全要素が1の行列
print np.ones([2,3])
#  [[ 1.  1.  1.]
#  [ 1.  1.  1.]]

#単位行列 
print np.eye(2,2)
#  [[ 1.  0.]
#  [ 0.  1.]]

#対角行列
print np.tri(3)
#  [[ 1.  0.  0.]
#  [ 1.  1.  0.]
#  [ 1.  1.  1.]]

#1-5までの連続値のベクトル
print np.arange(1,5)
#  [1 2 3 4]

#三つ目の引数はベクトル間の値になる
print np.arange(1, 10, 3)
#  [1 4 7]

#-1で逆方向のベクトルも作れる
print np.arange(10, 1, -1)
#  [10  9  8  7  6  5  4  3  2]

#対角成分の抽出
print np.diag(np.array([[0,1,2], [3,4,5], [6,7,8]]))
#  [0 4 8]

#乱数の生成 (最小値,最大値,個数)
print np.random.randint(0,100,10)
#  [47 43 37  8 96 83 69 10 57 72]

 

配列の情報取得

配列やベクトルの

サイズやデータ量などの情報は

下記のように取得できます。

print "====Data Info===="
a = np.array([[1,2,3],[4,5,6],[7,8,9]])
print a
#  [[1 2 3]
#  [4 5 6]
#  [7 8 9]]

print a.ndim ## 次元数
#  2

print a.size ## 要素数
#  9

print a.shape ## 各次元の要素数 (行数, 列数)
#  (3, 3)

print a.itemsize ## 1要素のバイト数
#  8

print a.strides ## 次の行、と次の列までのバイト数
#  (24, 8)

print a.nbytes ## 配列全体のバイト数 (size*itemsize)
#  72

print a.dtype ## 要素のデータ型(後述)
#  int64

統計量の計算

下記のように、行列やベクトルの各要素を使って、

基本的な統計量も簡単に計算できます。

print "====Statics===="
a = np.array([[1],[2],[3]])
print a
#  [[1]
#  [2]
#  [3]]

print a.mean()     # 平均 np.mean(a)
#  2.0

print a.var()      # 分散 np.var(a)
#  0.666666666667

print np.median(a) # 中央値 
#  2.0

print a.std()      # 標準偏差 np.std(a)        
#  0.816496580928

print a.ptp()      # 値の範囲(最大値-最小値) np.ptp(a)        
#  2

print a.sum()      # 総和 np.sum(a)
#  6

print a.min()      # 最小値 np.min(a)
#  1

print a.max()      # 最大値 np.max(a)
#  3

#相関係数
a = np.random.randint(0,500,20)
b = np.random.randint(0,500,20)
print np.corrcoef(a,b)
#  [[ 1.         -0.16742262]
#  [-0.16742262  1.        ]]

#ヒストグラム (頻度とビン)
print np.histogram(a)
#  (array([4, 1, 0, 1, 1, 4, 0, 2, 3, 4]),
#  array([   4. ,   50.2,   96.4,  142.6,  188.8,  235. ,  281.2,  327.4,
#  373.6,  419.8,  466. ]))

#累積和
print np.cumsum(a)
#  [ 432  827  843  891 1079 1266 1319 1419 1840 1895 1978 2024 2164 2563 2817
#  3316 3434 3856 4087 4449]

#累積
print np.cumprod(a)
#  [                 432               170640              2730240
#  131051520          24637685760        4607247237120
#  244184103567360    24418410356736000 -8166593313523695616
#  -6440774474774020096   371296731333328896 -1367094432376422400
#  -6925779795603619840  3625472610588426240 -1467160596017315840
#  5756625535741460480 -3247717509761073152 -5477727664666050560
#  7470250548101382144 -7440680422603751424]


 

基本的な行列演算

下記のサンプルのように、

固有値や逆行列など、

基本的な行列演算をすることができます。

 

print "======Matrix===="

a = np.array([[1.,0.], [0.,2.]])
print a
#  [[ 1.  0.]
#  [ 0.  2.]]

la, v = np.linalg.eig(a)    # 行列Aの固有値・固有ベクトル
print la    #固有値
#  [ 1.  2.]

print v     #固有ベクトル
#  [[ 1.  0.]
#  [ 0.  1.]]


a = np.matrix([[1,2],[3,4]])
print a
#  [[1 2]
#  [3 4]]

print a.T   #転置
#  [[1 3]
#  [2 4]]

print a.I   #逆行列
#  [[-2.   1. ]
#  [ 1.5 -0.5]]

#行列式
print np.linalg.det(a)
#  -2.0

v1 = np.array((1,0,0))
v2 = np.array((0,1,0))

#内積
print np.dot(v1,v2)
#  0

#外積
print np.cross(v1,v2)
#  [0 0 1]

#ベクトルの長さ
print np.linalg.norm(a)
#  5.47722557505

 

行列の連結

MATLABで非常に便利な機能として、

行列に次々とベクトルを連結して、

時系列のデータを格納する機能があると思います。

 

Numpyでは、

下記のサンプルコードのように、

行方向に連結する場合は、c_ (coloumのc)

列方向に連結する場合は、r_ (rowのr)

という関数で連結できます。 

# 行列の連結
b = np.array([[1],[2],[3]])
print b
#[ 0.  0.  0.]

a = np.empty(3)
for i in range(5):
    a = np.c_[a,b] #連結

print a
#[[ 0.  1.  1.  1.  1.  1.]
# [ 0.  2.  2.  2.  2.  2.]
# [ 0.  3.  3.  3.  3.  3.]]

CSVファイル出力

下記のように、

numpyの行列データは

簡単にCSVファイルに書き込んだり、

逆にCSVファイルから読み込んだりすることができます。

 
## CSVファイル出力

a = np.array([[1],[2],[3]])

#csv形式で保存
np.savetxt("array.csv", a ,  delimiter=",")

#csv形式のファイルを読み込み
arr = np.loadtxt("array.csv", delimiter=",")

 

MATLABユーザ向けNumpy/Scipyの使い方と注意点

下記のScipyのページに、

MATLABユーザ向けのNumpy行列演算の注意点がまとまっています。

 

上記の説明を元に、MATLABユーザ向けの

Numpyの行列演算の注意点をまとめておきたいと思います。

注意点1: 行列データの型

MATLABでは、すべての行列データは、

doubleの多次元配列として表されますが、

Numpyでは、様々な行列データの型があります。

 

Numpyの行列では、浮動小数点型以外にも

整数型やbool型なども格納できます。

 

注意点2: numpy arrayにおける行列の掛け算

numpyのarrayを

通常の掛け算オペレータ*(アスタリスク)で計算すると、

いわゆる行列同士の掛け算ではなく、

要素同士の掛け算になります。

a = np.array([[1,2,3],[4,5,6],[7,8,9]])
print a
#[[1 2 3]
#[4 5 6]
#[7 8 9]]

b = np.array([[1],[2],[3]])
print b
#[[1]
# [2]
# [3]]

print a*b
#  要素毎の掛け算
#  [[ 1  2  3]
#  [ 8 10 12]
#  [21 24 27]]

oppython.hatenablog.com

 

もし行列の掛け算にしたい場合は、

dotメソッドを使う必要があります。

print a.dot(b)
#  行列の掛け算
#  [[14]
#  [32]
#  [50]]

注意点3: インデックスが0始まり

MATLABはインデックスの始まりが1からですが

Pythonは多くのプログラミング言語の配列インデックスと

同様に0始まりです。

基本的ですが、デバックしずらいバグが

入る可能性があるので注意しましょう。

 

注意点4: 行列をスライスしたものはポインタで表される

MATLABでは、

行列からある一列を抽出した場合、

その一列はデータがコピーされ別の変数になりますが、

Numpyの場合、前述のスライスを使って

新しい変数を作っても、そのデータはポインタで

元のデータを参照していることが多いです。

 

ですので、下記のサンプルでは、

atの値を変更するとaの値も変更されてしまいます。

MATLABのようにデータのコピーをして欲しい場合は、

copyメソッドを使いましょう。

at = a[1:3,:]
print at
#  実はポインタ
#  [[4 5 6]
#  [7 8 9]]

atcopy = a[1:3,:].copy()
print atcopy
#  データのコピー
#  [[4 5 6]
#  [7 8 9]]

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

プログラミングのための線形代数

プログラミングのための線形代数