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つがある
- numpy.linalgの関数を使う。
- scipy.linalgの関数を使う。
- scipy.sparse.linalgの関数を使う。
numpy.linalgとscipy.linalgには以下の4つの関数がある。
- eig:一般の行列の固有値・固有ベクトルを求める。
- eigh:エルミート(or 実対称)行列の固有値・固有ベクトルを求める。
- eigvals:一般の行列の固有値のみを求める。
- eigvalsh:エルミート(or 実対称)行列の固有値のみを求める。
関数名の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の略。一般化固有値問題も出来る。
続きを読む外積
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)
他にも、配列の次元に依って、いくつか方法がある。
もっとも、単純な場合である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を使っとけばよさげ。
モジュールの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で特異値分解
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つ目の方法の方が安定するようです。