MATLABのプログラムをPythonに移植する際のTIPS・注意点

  • 34
    いいね
  • 0
    コメント

年に1度ぐらいの頻度で必要に迫られてMATLABのプログラムをPythonに移植せざるを得ないことがあるのだが、その際に気づいたTIPSや注意点をまとめておく。随時追加予定。
(以下の内容には、MATLABとの相違点に関するNumPy/SciPyの公式ドキュメントNumpy for MATLAB usersの内容と重複するものが多く含まれる。)

転置

行列Mの転置は、MATLABでは M' 、NumPyでは M.T 。

配列A, Bに対する [A B]

MATLABでの配列A,Bに対する [A B] という操作は、NumPyでの numpy.hscack([A, B]) に相当する。(AとBは、最後の次元以外ではサイズが同じである必要がある。)

size()

配列Aの各次元のサイズはMATLABでは size(A) で求めるが、NumPyでは numpy.shape(A) で求める。numpy.size(A) とすると、配列Aに含まれる要素数になってしまうので注意。

配列のindex

MATLABでは配列の要素のindexが1から始まり、NumPyでは0から始まるというのは誰もが移植の際に注意していることであろうが、配列の次元に対するindexの違いにも注意する必要がある。すなわち、MATLABでの sum(A, 2) は配列Aの2番目の次元に沿って和をとるが、これはNumPyでは numpy.sum(A, axis=1) としなければならない。

特定の次元に沿って演算を行う関数(sum(), mean(), std(), var(), max(), min(), ...)

MATLABでの sum(A) は配列Aの最初の次元に沿った和を取るが、NumPyでは(axisを指定しなければ)配列Aのすべての要素の和をとる。よって、MATLABでの sum(A) はNumPyでの numpy.sum(A, axis=0) に相当することになる。この挙動は mean(), std(), var(), max(), min() 等の、特定の次元に沿って演算を行う他の関数でも同様。

関数適用後の配列のsize/shape

上記に関連して、sum(), mean(), std(), var() 等で特定の次元に沿って演算を行う際、結果の配列がどのような次元を持つかについてMATLABとNumPyで違いがあり、注意する必要がある。例えば、「3次元配列Aの2番目の次元に沿って平均を取り、それを元の配列から引く」という計算をしたい場合、MATLABでは A - mean(A, 2) で求められるが、NumPyで A - numpy.mean(A, 1) とするとエラーが出る。これはMATLABでの mean(A, 2) の結果は3次元配列である(平均を取った次元のサイズは1になる)のに対し、NumPyでの numpy.mean(A, 1) の結果は2次元配列である(平均を取った次元が消去される)ため、3次元配列との減算が(一般的には)できないからである。この場合、NumPyでは A - numpy.mean(A, 1)[:, numpy.newaxis, :] として、消去された次元を numpy.newaxis で明示的に補ってやる必要がある。

reshape()

配列の要素を変えずに形だけを変えたいとき reshape() を使うが、例えば配列Aを [1, 2, 3, 4] としてこれを2x2の2次元配列に並べ替える際、MATLABで reshape(A, 2, 2) とすると [1 3; 2 4] が生成されるのに対し、NumPyで numpy.reshape(A, (2, 2)) とすると [[1, 2], [3, 4]] が生成される。つまり、MATLABではまず最初の次元にそって要素を埋めていくが、NumPyではまず最後の次元に沿って要素を埋めていく。NumPyでMATLABのような順番で要素を埋めたい場合、numpy.reshape() のキーワード引数 order に 'F' を与える必要がある。

配列Aに対する A(:)

上記とも関連するが、MATLABでの多次元配列Aに対する A(:) という操作は、Aの全要素を1次元に並べなおした列ベクトルを返すが、この際、MATLABは最初の要素からまず最初の次元に沿って要素を取り出して並べる。これに対し、NumPyの numpy.ravel(A) も同様に多次元配列を1次元配列に並べなおすが、この際、NumPyは最初の要素からまず最後の次元に沿って要素を取り出して並べるため、MATLABの A(:) の結果とは要素の順番が異なってしまう。NumPyでMATLABのような順番で要素を並べたい場合、numpy.ravel() のキーワード引数 order に 'F' を与える必要がある。

zeros(), ones()

MATLABでの zeros(n) はサイズが n x n で、すべての要素の値が0の2次元配列を生成するが、NumPyの numpy.zeros(n) はサイズ n の1次元配列を生成する。よって、MATLABでの zeros(n) はNumPyでの numpy.zeros((n, n)) に相当することになる(生成される配列の shape をタプルで与える必要があることに注意)。引数が2つ以上の場合はMATLABでもNumPyでも同様の挙動をする。つまり、MATLABでの zeros(m, n) もNumPyの numpy.zeros((m, n)) もサイズが m x n の2次元配列を生成する。この挙動は ones() でも同様。

rank()

行列MのランクはMATLABでは rank(M) で求めるが、NumPyでは numpy.linalg.matrix_rank(M) で求める。numpy.rank(M) とすると、配列Mの次元数(つまり M.ndim すなわち len(M.shape))になってしまうので注意。

eig()

行列Mの固有値はMATLABでは eig(M) で求めるが、NumPyでは numpy.linalg.eigvals(M) で求める。numpy.eig(M) とすると、固有値を要素に持つベクトルに加え、固有ベクトルを列に持つ行列も同時に受け取ることになるので注意(w, V = numpy.eig(M) で、ベクトルwが固有値、行列Vが固有ベクトル)。MATLABでも、eig(M) から戻り値を2つ受け取るようにすると固有値と固有ベクトルが得られるが、この場合、固有ベクトルを列に持つ行列と、固有値を対角成分に持つ対角行列、という形で受け取ることになるので注意([V,D] = eig(M) で、行列Vが固有ベクトル、行列Dが固有値)。

svd()

行列Mの特異値分解にはMATLABでは svd(M) を、NumPyでは numpy.linalg.svd(M) を使う。MATLABで戻り値を3つ受け取るようにすると行列 U, S, V が得られ、Sが特異値を対角要素に持つ行列、U, V がそれぞれ左・右特異ベクトルを列に持つ行列である。NumPyのsvd(M)からは配列 u, d, v が得られ、u が左特異ベクトルを列に持つ行列、d が特異値を要素に持つ1次元配列 、v が右特異ベクトルを行に持つ行列である。つまり、MATLABの U とNumPyの u は一致するが他は一致せず、S は numpy.diag(d) に、V は v の転置(複素行列の場合はエルミート転置)に相当するので注意が必要。

cov()

多変数のデータ行列から変数間の共分散を計算する関数だが、MATLABとNumPyではデータ行列の行と列の解釈が逆になっている。すなわち、MATLABでは各が各変数に対応するが、NumPyでは各が各変数に対応する。よって、MATLABの cov(X) はNumPyの numpy.cov(X.T) に相当することになる。

linspace(a, b, n)

MATLABとNumPyでほとんど同じ挙動をするが、nが1の際に、MATLABではbを返すが、NumPyではaを返すという違いがある。

構造体

MATLABの構造体をPythonで表現するにはいくつかの方法が考えられるが、NumPyのstructured arrayを用いるのが便利であるケースが多いと思う。ただし、MATLABの構造体は動的にフィールドを追加できるのに対し、NumPyのstructured arrayは宣言の際にすべてのフィールドを定義しておかなければならないことに注意。

構造体に対する[・]

MATLABでの構造体Sのフィールドfに対する [S.f] という操作は、NumPyのstructured array Sのフィールドfに対する S['f'] に相当する。ただし、S['f']が多次元配列の場合は、numpy.hstack(S['f']) とする必要がある(つまり、最初の次元に関して配列を分解し、最後の次元の方向で結合しなおす)。