おっぱいそん!

pythonを数値計算に使うときのテクニックとかをまとめていきたい。

Pythonの代入の挙動

Pythonで代入をした時に、各変数がどのオブジェクトを指すかをまとめておく。特に、list(リスト)とNdarray(配列)の違いをば。
前回の記事も参照 : http://oppython.hatenablog.com/entry/2015/01/15/014729

例1

前回の記事でも書いたようにPythonでは、代入(=)は変数に値(=オブジェクト)を格納しているのではなく(値渡し)、値(=オブジェクト)への参照を格納している(参照渡し)。

x = np.array([1])
y = x
print 'id(x)_before =', id(x)
print 'id(y)_berore =',  id(y)
x = np.array([2])
print 'id(x)_after =', id(x)
print 'id(y)_after =',  id(y)
print y

id(x)_before = 4300100832
id(y)_berore = 4300100832
id(x)_after = 4300101152
id(y)_after = 4300100832
[1]

y = xとすると、変数yもxと同じオブジェクトnp.array([1])を指すようになる(参照渡し)。
変数xに新しいオブジェクトnp.array([2])を代入すると、xは新しいオブジェクトを参照するようになる。
もちろん、yには何もしていないので、始めのオブジェクトを参照したままである。
上の挙動は、ndarrayだけでなく、タプル、文字列、リスト、数値などでも共通(型が変更可能かどうかに依らない)。

続きを読む

Pythonの型の分類と代入

変数が指すオブジェクトがどれかなどの振る舞いをすぐ忘れてしまうのでメモしておく。

組み込み型の分類

まず、Pythonの組み込み型についてまとめておく。

管理方式

管理方式の違いで分類すると、以下のようになる。

  • コレクション型(=コンテナ型):複数のデータをまとめたもの

+ シーケンス型:要素を順番に並べて管理。整数のindexでデータを管理できる。
例)文字列、リスト、タプル
+ マップ型:キーと値のペアで要素を管理。
例)ディクショナリ

  • 非コレクション

例)数値型(int, float, long, complex)、ブール型、None

文字型がシーケンス型なのは面白い。シーケンスなので、たとえば、

a = "oppai"
print a[4] #i

のようにスライスを使って、要素の指定とかができる。

変更可能性

変更可能性で分類すると、以下のようになる。

  • 変更可能 (mutable)

例)リスト、ディクショナリ

  • 変更不能 (immutable)

例)数値型、文字列、タプル

基本的に、代入なりなんなりをした時の挙動は、変更可能かどうかで決まる。
オブジェクトはメモリに保存されているが、変更可能ならば、そのメモリに書かれているオブジェクトは書き換え可能になる。

続きを読む

numpy.matrixとnumpy.ndarrayの違い

numpy.matrixと2次元のnumpy.ndarrayは同じことができるが、numpy.matrixでできることはすべてnumpy.ndarrayでできる。
numpy.matrixがある理由はたぶん、行列を扱う時にはnumpy.ndarrayより見やすいので。
行列以外の高次のテンソルを扱う時には、混乱を避けるため、numpy.ndarrayだけ使うようにしたほうがよいと思われ。

個人的には、行列以外のテンソルを扱うことも多いし、numpy.matrixを使う理由はなさそう。

Ref:
http://wiki.scipy.org/NumPy_for_Matlab_Users
http://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/linalg.html

固有値問題

Python固有値問題を解く方法についてメモしておく。
メジャーな方法として、以下の3つがある

  1. numpy.linalgの関数を使う。
  2. scipy.linalgの関数を使う。
  3. scipy.sparse.linalgの関数を使う。


numpy.linalgとscipy.linalgには以下の4つの関数がある。

関数名のhはHermitianの略。Scipyだと一般化固有値問題*1もオプションで出来る。
ちなみに、scipy.linalgはnumpy.linalgに含まれる関数はすべて含んでいて、さらに追加で他の関数も含んでいる。また、scipy.linalgが常にBLAS/LAPACKコンパイルされるのに対し、numpy.linalgはそうとは限らないので、scipyの方が早い場合もある。なので、特別な理由(コードが既にnumpy.linalgで書かれていて、書きなおすのがめんどいとか)がない限り、scipy.linalgを使うほうがよい。たぶん。
http://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/linalg.html


scipy.sparse.linalgには以下の2つの関数がある。

関数名のsはsparseの略。一般化固有値問題も出来る。

*1:普通の固有値問題がAx=axに対し、一般化固有値問題はAx=aMxのような固有値問題のこと。

続きを読む

外積

Numpyで外積を計算する方法をまとめておく。

ベクトル(1次元配列)同士の外積、つまり
c[i,j] = a[i]b[j]
はnumpy.outerを使って、
http://docs.scipy.org/doc/numpy/reference/generated/numpy.outer.html

import numpy as np
a = np.ones(3)
b = np.arange(2)
print np.outer(a,b)

とできる。


高次元配列の外積
numpy.multiply.outer もしくは numpy.einsumを使うといい。
numpy.multiply.outerはmultiplyの部分を別の演算子に変更することも出来る。
参考:
http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.outer.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html
http://stackoverflow.com/questions/24839481/python-matrix-outer-product

例:c[i,j,k] = a[i]b[j,k]

a = np.ones(2)
b = np.arange(12).reshape((3,4))

print a.shape, b.shape
print np.multiply.outer(a,b).shape
print np.outer(a,b).reshape(2,3,4).shape
print np.einsum('a,bc->abc',a,b).shape

(2,) (3, 4)
(2, 3, 4)
(2, 3, 4)
(2, 3, 4)

要素積を任意の軸について

要素積 A*B or np.multiply(A,B) をndarray A, B の任意の軸について実行する方法。
(基本的な使い方は、以前の記事http://oppython.hatenablog.com/entry/2014/01/05/004454で紹介したので、そちらを参考に。)

np.multiply(A, B)はn次元配列Aと1次元配列Bを掛けると、Aの最後の軸とBの要素積をする。:
C[i,j,k] = A[i,j,k] B[k]
したがって、
C[i,j,k] = A[i,j,k] B[j]
のような、最後の軸以外の軸と要素積をしたい場合には、そのままでは使えない。

そこで、少し工夫して、最後の軸以外の軸との要素積を実装する。

もっとも、メジャーな方法は、1次元配列にnp.newaxis*1を使って軸を追加する方法。
np.newaxisはあるn次元配列B(例えば B.shape=(4,5,3) )に対し、Bnew = B[:,np.newaxis,:,np.newaxis,:]とすることで、Bnew.shape = (4,1,5,1,3)のように軸を追加できる機能。
例えば、上のC[i,j,k] = A[i,j,k] B[j]の例であれば、

import numpy as np

A = np.arange(24).reshape((2,3,4))
B = np.array([0,2,3])

print A.shape
print B.shape
print (A*B[np.newaxis,:,np.newaxis]).shape  # or (A*B[:,np.newaxis]) 

(2, 3, 4)
(3,)
(2, 3, 4)

とすれば良い。
参考:
http://stackoverflow.com/questions/7096371/in-numpy-what-is-the-fastest-way-to-multiply-the-second-dimension-of-a-3-dimens


他にも、配列の次元に依って、いくつか方法がある。
もっとも、単純な場合である2次元配列A(つまり行列)と1次元配列B(つまりベクトル)について、Aの行とBの要素積をやり方を以下に列挙する(どれも同じ結果を返す)。

A = np.arange(6).reshape((2,3))
B = np.array([0,2])

print A
print B
print A*B[:,None] # or A*B[:,np.newaxis] or B[:,None]*A
print A*np.expand_dims(B,1) 
print (A.T * B).T
print np.dot(np.diag(B), A)
print np.einsum('ij,i->ij',A,B)

[[0 1 2]
[3 4 5]]
[0 2]
[[ 0 0 0]
[ 6 8 10]]
....

None、np.newaxis、np.expand_dimsはどれも軸の追加で、同じもの。
転置とdot+diagは配列の次元が上がったり、要素積したい軸が変わると使えなくなる。
np.einsumはアインシュタインの縮約で、これ以外にも色々使い方があるみたい。

速度は方法に依って異なる。
速度の点では、
np.einsum > 転置 = np.newaxis > dot+diag の順で優れているっぽい。
参考:
http://stackoverflow.com/questions/18522216/multiplying-across-in-a-numpy-array


速度を気にしない時には、分かりやすいnp.newaxisを、速度が気になる時にはnp.einsumを使っとけばよさげ。