おっぱいそん！

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

外積

numpy

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を使っとけばよさげ。

*1:参考： http://d.hatena.ne.jp/se-kichi/20110214/1297647898 http://www.kamishima.net/mlmpyja/nbayes2/shape.html http://www.kamishima.net/mlmpyja/nbayes2/broadcasting.html

モジュールのimportの仕方

pythonでモジュールをインポートする書き方/仕方はいろいろあるので、まとめておく。

自分でfoofile.pyのなかに以下のような関数を作っていたとする。

#foofile.py
def func(x):
    return x

これをrun.pyで読み込んで使うには、

#run.py
import foofile
print foofile.func(3)

from foofile import func
print func(3)

from foofile import func as fun1
print fun1(3)

import foofile as foo
print foo.func(3)

など。
ただし、モジュールfoofile.pyとrun.pyは同じフォルダに置く。

別のフォルダにモジュールを置きたい場合

testフォルダの中にrun.pyと
folder1/foofile.pyを置く場合、

#run.py
import sys
sys.path.append('/home/name/Desktop/test/folder1')
import foofile as foo
print foo.func(2)

とする。

http://kannokanno.hatenablog.com/entry/20130503/1367571825
http://ameblo.jp/purigen/entry-10145562204.html

[2015.9.29 追記]
フォルダにモジュールを置く場合は__init__.pyを置かないとダメ。
pathを追加しなくても、いいみたい。
下の記事が分かりやすかった。
http://python.matrix.jp/pages/tips/import.html

pythonで特異値分解

numpy scipy

pythonで特異値分解(singular value decomposition,SVD)をする時のメモ。

一般の密行列のSVD

あまり選択肢がないみたい。とりあえず、Numpy or ScipyのSVDを使っとけば間違いなさそう。

  • numpy.linalg.svd(a, full_matrices=1, compute_uv=1)
  • scipy.linalg.svd

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html
http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.svd.html
行列AをSVDしたければ、

import numpy as np
X, Y, Z = np.linalg.svd(A)

import scipy.linalg
X, Y, Z = scipy.linalg.svd(A)

のように使います。
full_matrices=0を指定すれば、上の行列X,Zの次元が小さくなります。
compute_uv=1を指定すれば、X,Yを求めなくなります。
どちらも、LAPACKを使っているみたいです。
この2つに違いがあるのかは知りませんが、僕がテストしてみた範囲では速さや結果は変わらないみたいです。

ちなみに、疎行列用の関数scipy.sparse.linalg.svdsなどでも、密行列のSVDを出来るみたいです。

  • scipy.sparse.linalg.svds(A, k=Aの次元)

もちろん、疎行列のSVDをするなら、scipy.linalg.svdでなく、scipy.sparse.linalg.svdsや下の記事の関数を使う方が早いです。
一方、密行列のSVDはscipy.sparse.linalg.svdsより、scipy.linalg.svdの方が速いようなのですが、メモリはscipy.sparse.linalg.svdsの方が節約出来るみたいです(下のリンクを参照)。
http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/

疎行列のSVD

いろいろあるみたいですが、僕は使ってないので、下の記事を参考に。
幾つかの関数の速さを比較しています。
http://jakevdp.github.io/blog/2012/12/19/sparse-svds-in-python/
http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/

収束しないエラー

numpy.linalg.svdやscipy.linalg.svdを使っていると、

raise LinAlgError, ’SVD did not converge’
numpy.linalg.linalg.LinAlgError: SVD did not converge SVD

とエラーが出る時があります。
SVDは原理的には任意の行列に対して(エルミート行列の対角化を使って)実行できますが、numpy.linalg.svdでは高速化のために、"ある操作"を繰り返してSVDして求められる行列(に近い行列?)に収束させているようです。
その際、収束に必要な"ある操作"の繰り返しのステップ数が一定数n以上になると、諦めて上のようなエラーを吐くようです。
対策としては

  • numpy.linalg.svdモジュールの中の"ある操作"の繰り返しの上限nを書き換えて、収束が遅くてもエラーが出にくくする。
  • SVDする行列を2乗してからSVDする。
  • 対称（orエルミート）行列の対角化を使う。

2つ目の方法のように、2乗してからSVDすると、なぜエラーがでなくなるかは知りません。
（そもそも、モジュールの中でどういうアルゴリズムでSVDしているのか知らないので…とりあえず、使う分にはエラーさえ出なければ困らない）
1つ目の方法は、モジュールの中身を書き換えられない環境の時には使えない。
2つ目の方法は以下のようにする。

import numpy as np
try:
    X, Y, Z = np.linalg.svd(A)
except:
    A2 = np.dot(A.T, A)
    X2, Y2, Z = np.linalg.svd(A2)
    Y = np.sqrt(Y2)
    X = np.dot(A, Z.T); X = np.dot(X, np.linalg.inv(np.diag(Y)))

ここで、何回もSVDする場合には、毎回2乗していると遅くなるので(エラーはたまにしかでないと想定している)、普通にSVDを試みてから、エラーが出たら2乗してSVDするようにしている。
1つ目の方法と2つ目の方法のどちらが高速化は分からない。

3つ目の方法はエルミート行列np.dot(A.T ,A)の対角化に落とす方法です。
（任意の行列Aに対して、特異値分解ができることの証明としてよく使われる方法をそのままコードにした）

try:
	X, Y, Z = np.linalg.svd(A)
except:
	w,v = np.linalg.eigh(np.dot(A.T ,A))
	w = w[::-1]; v = v[:,::-1]
	Y = np.sqrt(w)
	X = np.dot(A,v); X = np.dot(X,np.diag(Y**(-1))); Z = v.T

ただし、Aの形が正方行列出ない時には、形に応じて、上のnp.dot(A.T ,A) or np.dot(A ,A.T)を使わなくてはいけません。適当に書き換えてください。

2つ目の方法より、3つ目の方法の方が安定するようです。

Fancy Indexing：配列の一部を条件をつけて取り出す

numpy

aを適当なN次元配列とする。
a[a_1, a_2, a_3, ]とすると、配列の1つの成分を取り出せる。
a[リスト]とすると、1次元配列のリスト番目の成分のみ取り出した配列を返す。
(上の用に配列の１つの成分だけ取り出したい時にはtuple(リスト)のようにtupleにして渡せば良い）
a[リスト1,リスト2,リスト3, ]とすると、1次元目をリスト1で2次元目をリスト2で、といったように取り出した配列を返す。
np.ix_をつかってa[np.ix_([行],[列])]とすると行列の一部の行と列のみ取り出した新しい行列を作れる。
※ここまでのリストは配列にしても同じ。
同じことをbool値を使って、a[[行のbool],[列のbool],[3次元目のbool], ]とも出来る(これはリストにすると、False=0, True=1と解釈されるので注意)。

import numpy as np
a = np.array([[1,2,3],[4,5,6],[7,8,9]])
print a[1,2]
print a[[1,2]]
print a[[1,2],[0,2]]
print a[np.ix_([1,2],[0,2])]
print a[np.array([False,True,True]),np.array([True,False,True])]

[[1 2 3]
[4 5 6]
[7 8 9]]
6
[[4 5 6]
[7 8 9]]
[4 9]
[[4 6]
[7 9]]
[4 9]

http://ibisforest.org/index.php?python%2Fnumpy#l9d00f8f

Numpyの配列（array）の一部を条件をつけて取り出す方法
import numpy as np
A = np.arange(10,0,-1)
print A>5
print np.where(A>5)
print A[A>5]**2
print A[np.where(A>5)]**2

[ True True True True True False False False False False]
(array([0, 1, 2, 3, 4]),)
[100 81 64 49 36]
[100 81 64 49 36]

配列に対して、論理式（>とか）を使うと、bool型（True or False）の配列を作ってくれる。
それを配列の引数として入れることで、配列の条件を満たす部分だけ取り出せる。
np.whereは基本的にタプルで条件が真の成分の場所を返す（返す値を別のものに指定も出来る）。
np.where(A>5)[0]とすれば、成分の場所を表すリストのみ渡してくれる。

【2014.10.23 追記】
2次元配列の例を追加しておく。

A = np.arange(9).reshape(3,3)
print A
print A>4
print np.where(A>4)
print A[A>4] # or A[np.where(A>4)]

[[0 1 2]
[3 4 5]
[6 7 8]]
[[False False False]
[False False True]
[ True True True]]
(array([1, 2, 2, 2]), array([2, 0, 1, 2]))
[5 6 7 8]

やはり、論理式は同じ形の配列で、成分がbool型のものを返す。
一方、np.whereはTrueの成分の場所を返す。今の場合は、Trueが4つあり、2次元配列なので、1つの成分の場所を指定するのに2つの数字がいるので、合計で8つの数字が並んでいる。
元の配列に論理式 or np.whereを入れると、元の配列が2次元でも、1次元配列として返す。

参考：
http://seesaawiki.jp/met-python/d/array
http://rest-term.com/archives/2999/
上の2つのサイトは配列の扱いに詳しい。
ひと通り読むといいかも。
http://d.hatena.ne.jp/nishiohirokazu/20111118/1321611727
http://www.geopacific.org/opensourcegis/python/python_raster

numpyのテンソル（配列）関係

numpy
numpy.transpose(a, axes=None)

配列の足の順番を入れ替える。
a：入れ替えたい配列
axes：順番の指定（指定なしなら、逆順になる）

import numpy as np

a = np.arange(24).reshape((2,3,4))
print a.shape
print np.transpose(a).shape
print a.transpose().shape
print np.transpose(a,(1,0,2)).shape

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

np.transpose(a)でもa.transpose()でも同じ役割
axesを指定しなければa[a_1,a_2,...,a_N]をa[a_N,...,a_2,a_1]に入れ替える。
http://docs.scipy.org/doc/numpy/reference/generated/numpy.transpose.html


掛け算

a*b = np.multiply(a,b) は要素積
np.dot(a,b)は行列積(a,bが高次のテンソルの場合はaの最後の足とbの最後から2番目の足について掛ける)。
他にもいろいろな種類の掛け算が用意されてる。

multiply or *は基本的にa,bの配列の次元は同じでなければならない。
c[i,j,k,l] = a[i,j,k,l]b[i,j,k,l]
同じでない場合は、片方がスカラー
c[i,j,k,l] = ab[i,j,k,l]
もしくは、1次元配列
c[i,j,k,l] = a[l]b[i,j,k,l] #1次元配列と一方の配列の最後の要素の掛け算になる
でないとだめっぽい。
また、a*bとb*aは同じものを返すっぽい。

import numpy as np

A = np.arange(24).reshape((2,3,4))
B = np.array([1,2,3,4])
print A.shape
print B.shape
print (A*B).shape
print (B*A).shape
print np.multiply(A, B).shape

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

http://ibisforest.org/index.php?python%2Fnumpy
http://docs.scipy.org/doc/numpy/reference/generated/numpy.multiply.html

【追記】任意の軸に関する要素積について、別に記事を書いたので、必要なら参照：
http://oppython.hatenablog.com/entry/2014/10/23/024831