目次
- 目次
- Pythonにおける行列演算
- NumpyとScipy
- Pythonで行列演算のサンプルプログラム
- MATLABユーザ向けNumpy/Scipyの使い方と注意点
- 参考資料
- MyEnigma Supporters
Pythonにおける行列演算
先日紹介したmatplotlibという
グラフ描画ツールを上手く使うと、
PythonをMATLABのように
アルゴリズムのプロトタイプ作成
に使いやすくなります。
また、PythonもMATLABに比べると計算演算も早いので、
そこまで計算コストが大きく無い場合、
C++でアルゴリズムを実装する代わりに
Pythonを使うことができます。
様々なアルゴリズムを実装する際に、
重要なのは行列演算機能です。
行列演算の四則演算や、逆行列・固有値の計算を
通常のプログラミング言語で実装するは骨が折れます。
しかし、
MATLABには行列演算機能がデフォルトでbuild-inされていますし、
C++にはEigenという非常に使いやすいライブラリがあります。
実は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.]]
また、hstack, vstackという関数を使うと同じように連結ができます。
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]]
もし行列の掛け算にしたい場合は、
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 Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。
