Matlab/OctaveユーザーのためのPython NumPyライブラリ概要

私のブログからの転載です。プログラミングやソフトウェア周りの記事を書いているのでご興味あればのぞいてみてください。

この記事のまとめ

• MATLAB/OctaveとPythonのNumPyライブラリ比較
• NumPyのarray型とmatrix型の比較

背景

最近の私がよく使うプログラミング言語はPythonとMATLABなのですが、基本的にPythonで算術計算をあまりすることはなく、行列計算などを必要とするプログラムについては、MATLABでしかやったことがありませんでした。

機械学習を学んでいると、Pythonで行列計算などの算術計算をする必要がありそうなので、MATLABと対比しながら、Pythonでどのように行列計算を行っていくかをまとめたいと思います。

基本的にはSciPyの公式ページにまとまっている記事があるので、要点だけ抽出してこちらに日本語で記載していきたいと思います。

MATLABとPython NumPyの重要な違い

対比して簡単にまとめると次の通りです。

MATLAB	NumPy
基本のデータ型がdouble精度の浮動小数点型の多次元配列。2次元インスタンスにおける演算は線形代数の行列演算用に設計されている。	NumPyではarray型が多次元配列。通常の演算子による演算は要素ごとの演算となる。ただし、線形代数用にarray型のサブクラスとしてmatrix型もある。matrixクラスでの演算は線形代数演算となる。
配列の要素のインデックスが1から始まる	配列の要素のインデックスが0から始まる

他にもいろいろありますが、基本的な概念のとして、MATLABはプロトタイピングに適した言語であって、非常にわかりやすく、使い易い反面、メモリー消費が激しかったりして、それを商用的に使うような言語ではありません。

それに対して、Pythonは一般的にはスクリプト言語でもあるため、実行速度は単純なコーディングでは大して早くないですが、NumPyなどの科学計算用のライブラリを備えながら、並列計算やGPGPU用などのさまざまなライブラリやコンパイラなど高速化できる要素が出てきていることもあって、いわゆるデータサイエンティストが機械学習用に商用利用しても遜色のない言語になってきています。

そのあたりを前提として、それぞれの違いを考えた方がよいとは思います。

詳細は下記の記事などを読むとよくわかると思います。

• どこまで速くできる？ 達人に学ぶPython超高速データ分析～PyData.Tokyo Meetup #4イベントレポート

array型を使うべきか？matrix型を使うべきか

端的に言うとarray型を使うべき

端的な理由は下記の通りです。
1. 多くのNumPy関数はarray型を返す
2. 要素ごとの演算と線形代数演算が明確に分かれている
3. 配列として扱うことも、行ベクトル、列ベクトルも扱うこともできる

詳細に書くと下記の通りです。

array型	matrix型
*演算子は要素ごとの掛け算。dot()関数が行列積演算として使われる。	*演算子は行列積演算。multiply()関数が要素ごとの掛け算として使われる。
Nx1、1xNのベクトルを扱うこともできるし、N個の1次元配列を扱うことができる。例えば、A[:,1]とした場合、N個の1次元配列として扱われる（Nx1の列ベクトルではない）。	1次元配列は、Nx1、1xNの列ベクトルや行ベクトルには変換されない。A[:,1]とした場合、Nx1の列ベクトルとして扱われる。
2次元より大きいオブジェクトを扱うことができる	常に2次元のオブジェクトのみを扱う
.Tによって転置ができる	.Tで転置だけでなく、.H、.I、.Aによって複素共役転置（随伴行列、エルミート転置）、逆行列、asarray()をできる
arrayコンストラクターはPythonの配列表現で初期化する: array([[1,2,3],[4,5,6]])	matrixコンストラクターは文字列表現で初期化することもできる: matrix("[1 2 3; 4 5 6]")

なお、array型とmatrix型の演算を行うとmatrix型になります。

MATLABユーザーがよく使う行列表記とNumPy表記の比較

上記の通り、NumPyを使う場合にはarray型を使うべきなので、下記の比較はarray型のみを対象にしております。

もし、matrix型を使いたいという方は、NumPyの中にmatlibモジュールがあり、これにはMATLABユーザーにとってより使い慣れた関数表現が可能です。ただし、matlibモジュールの関数はmatrix型の返り値を取りますので、array型を使う方には適していませんのでご注意ください。

MATLAB	NumPy	備考
[1 2 3; 4j 5j 6j]	numpy.array([[1.,2.,3.],[4.j,5.j,6.j]])	2x3の複素行列表記
[a b; c d]	vstack([hstack([a,b]), hstack([c,d])])	配列の結合
A(end)	A[-1]	最後の要素（ただしPythonの場合は1次元目の最後の要素）
A(1,3)	A[0,2]	1行目3列目の要素
A(1,:)	A[0,:]	1行目のすべての要素
A(1:2,2:3)	A[0:2][1:3]	1-2行目の2-3列目の要素
A.'	A.T	転置
A'	A.conj().T	エルミート転置（複素共役転置、随伴行列）
A*B	A.dot(B)	行列積
A.*B	A*B	要素ごとの積
A.^3	A**3	要素ごとのべき乗
numel(A)	numpy.size(A) or A.size	配列の要素数
size(A)	numpy.shape(A) or A.shape	配列のサイズ
size(A,n)	A.shape[n-1]	配列のn次元目のサイズ
find(A>0.5)	nonzero(A>0.5)	0.5以上の要素のインデックス
C=A	C=A.copy()	値渡しによる代入（C=Aとすると参照渡しになります。部分参照の場合でも同様）
D=A(:)	D=A.flatten()	行列のベクトル化
1:10	arange(1.,11.) or r_[1.:11.]	1から10までの連続数の配列（MATLABは行ベクトル、列ベクトルが配列として扱われる）
1:10	arange(1.,11.).reshape(1,10)	1から10までの連続数の行ベクトル
[1:10]'	arange(1.:11.).reshape(10,1) or c_[1.:11.]	1から10までの連続数の列ベクトル
zeros(3,4)	zeros((3,4))	64-bit浮動小数点の0で満たされた3x4の行列
ones(3,4)	ones((3,4))	64-bit浮動小数点の1で満たされた3x4の行列
eye(3)	eye(3)	3x3の単位行列
diag(a)	diag(a)	配列aの要素を対角要素に持つ対角行列
rand(3,4)	random.rand(3,4)	ランダムな要素の3x4の行列
repmat(A,2,3)	tile(A,(2,3))	行列Aを行方向に2個、列方向に3個コピーした行列
max(max(A))	A.max()	行列Aのすべての要素の最大値
max(A)	A.max(0)	1次元目の配列の要素の最大値
inv(A)	linalg.inv(A)	逆行列
pinv(A)	linalg.pinv(A)	疑似逆行列
A\B	linalg.solve(A,B) if A is squared linalg.lstsq(A,B) otherwise	$A^{-1}B$を解く

とりあえず最低限必要なものはこの程度でしょうか。

必要に応じて適宜追加していきたいと思います。

もし、実際のコードを打ちながら試したいという方は本記事の最後にarray型とmatrix型でのコードを示しておきますのでお試しください。

今回は以上です。
最後まで読んでいただき、ありがとうございます。

サンプルコード

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Make all numpy available via shorter 'np' prefix
import numpy as np
# Make all matlib functions accessible at the top level via npm.func()
import numpy.matlib as npm

# 3x3 の複素行列表記
## array型
A = np.array([[1., 2., 3.], [4.j, 5.j, 6.j], [7., 8., 9.]])
print('A: \n{0}'.format(A))
## matrix型
M = np.matrix('[1. 2. 3.; 4.j 5.j 6.j; 7. 8. 9.]')
print('M: \n{0}'.format(M))

# 配列の結合
## array型
a = np.array([1, 1])
b = np.array([2, 2])
c = np.array([3, 3])
d = np.array([4, 4])
e = np.array([5, 5])
f = np.array([6, 6])
B = np.vstack( [np.hstack([a, b]), np.hstack([c, d]), np.hstack([e, f])] )
print('vstack([hstack([a,b]),hstack([c,d]),hstack([e,f])]): \n{0}'.format(B))
## matrix型
m = np.matrix('[1 1]')
n = np.matrix('[2 2]')
o = np.matrix('[3 3]')
p = np.matrix('[4 4]')
q = np.matrix('[5 5]')
r = np.matrix('[6 6]')
N = np.bmat('m n; o p; q r')
print("bmat('m n; o p; q r'): \n{0}".format(N))

# 最後の要素
## array型、matrix型共通
print('A[-1]: \n{0}'.format(A[-1]))

# 3行目1列目の要素
## array型、matrix型共通
print('A[2,0]: \n{0}'.format(A[2,0]))

# 1行目のすべての要素
## array型、matrix型共通
print('A[0,:]: \n{0}'.format(A[0,:]))

# 2-3行目の1-2列目の要素
## array型、matrix型共通
print('A[1:3,0:2]: \n{0}'.format(A[1:3,0:2]))

# 転置
## array型、matrix型共通
print('A.T: \n{0}'.format(A.T))

# エルミート転置（複素共役転置、随伴行列）
## array型
print('A.conj().T: \n{0}'.format(A.conj().T))
## matrix型
print('M.H: \n{0}'.format(M.H))

# 行列積
## array型
print('A.dot(B): \n{0}'.format(A.dot(B)))
## matrix型
print('M*N: \n{0}'.format(M*N))

# 要素ごとの積
## array型
print('A*A: \n{0}'.format(A*A))
## matrix型
print('multiply(M,M): \n{0}'.format(np.multiply(M,M)))

# 要素ごとのべき乗
## array型
print('A**3 \n{0}'.format(A**3))
## matrix型
print('power(M,M): \n{0}'.format(np.power(M,3)))

# 行列の要素数
## array型、matrix型共通
print('np.size(A): {0}'.format(np.size(A)))
print('A.size: {0}'.format(A.size))

# 行列のサイズ
## array型、matrix型共通
print('np.shape(A): {0}'.format(np.shape(A)))
print('A.shape: {0}'.format(A.shape))

# 1次元目の要素数
## array型、matrix型共通
print('A.shape[0]: {0}'.format(A.shape[0]))

# 0.5以上の要素のインデックス
## array型、matrix型共通
print('nonzero(A>0.5): {0}'.format(np.nonzero(A>0.5)))

# 値渡しによる代入
## array型、matrix型共通
C=A.copy()
C[0,0] = 99.
print('A: \n{0}'.format(A))
print('C: \n{0}'.format(C))

# 行列のベクトル化
D=A.flatten()
print('D: \n{0}'.format(D))

# 1から10までの連続数の配列
## array型のみ（matrix型は1次元配列の定義ができない）
print('arange(1.,11.): \n{0}'.format(np.arange(1.,11.)))
print('r_[1.:11.]: \n{0}'.format(np.r_[1.:11.]))

# 1から10までの連続数の行ベクトル
## array型
print('arange(1.,11.).reshape(1,10): \n{0}'.format(np.arange(1.,11.).reshape(1,10)))
## matrix型
print('matrix(arange(1.,11.)): \n{0}'.format(np.matrix(np.arange(1.,11.))))

# 1から10までの連続数の列ベクトル
## array型
print('arange(1.,11.).reshape(10,10): \n{0}'.format(np.arange(1.,11.).reshape(10,1)))
print('c_[1.:11.]: \n{0}'.format(np.c_[1.:11.]))
## matrix型
print('matrix(arange(1.,11.)).T: \n{0}'.format(np.matrix(np.arange(1.,11.)).T))

# 0で満たされた3x4の行列
## array型
print('numpy.zeros((3,4)): \n{0}'.format(np.zeros((3,4))))
## matrix型
print('numpy.matlib.zeros((3,4)): \n{0}'.format(npm.zeros((3,4))))

# 1で満たされた3x4の行列
## array型
print('numpy.ones((3,4)): \n{0}'.format(np.ones((3,4))))
## matrix型
print('numpy.matlib.ones((3,4)): \n{0}'.format(npm.ones((3,4))))

# 3x3の単位行列
## array型
print('numpy.eye(3): \n{0}'.format(np.eye(3)))
## matrix型
print('numpy.matlib.eye(3): \n{0}'.format(npm.eye(3)))

# 配列aの要素を対角要素に持つ対角行列
## array型
print('numpy.diag(a): \n{0}'.format(np.diag(a)))
## matrix型
print('numpy.matlib.diag(a): \n{0}'.format(npm.diag(a)))

# 乱数を要素に持つ3x4の行列
## array型
print('numpy.random.rand(3,4): \n{0}'.format(np.random.rand(3,4)))
## matrix型
print('numpy.matlib.random.rand(3,4): \n{0}'.format(npm.random.rand(3,4)))

# 行列のすべての要素の最大値
## array型、matrix型共通
print('A.max(): {0}'.format(A.max()))

# 1次元目(行)の要素について2次元目(列)以降ごとの最大値
## array型、matrix型共通
print('A.max(0): {0}'.format(A.max(0)))

# 行列のコピー
## array型
print('tile(A,(2,3)): \n{0}'.format(np.tile(A,(2,3))))
## matrix型
print('repmat(): \n{0}'.format(npm.repmat(M,2,3)))

# 逆行列
## array型、matrix型共通
print('numpy.linalg.inv(A): \n{0}'.format(np.linalg.inv(A)))

# 疑似逆行列
## array型、matrix型共通
print('numpy.linalg.pinv(A): \n{0}'.format(np.linalg.pinv(A)))

# A^-1*Bを解く
## array型、matrix型共通
print('numpy.linalg.solve(A,B): \n{0}'.format(np.linalg.solve(A,B)))