Pythonで数値計算１

研究をしていると巨大データに対して数値計算を行うという場面が多々あります。数値計算と言えばMatlabが有名ですが、Pythonでも数値計算はできます。ただ、純粋なPythonコードを書くとFortran, C, C++などのコンパイラ言語と比較して許容できないくらい遅くなってしまうので、CやC++で書かれた拡張モジュールを用いることになります。

Pythonで数値計算をするためのモジュール

Pythonで数値計算をするためのモジュールにはNumeric、numarray、numpyの３つがあります。numpyが最も新しく、他２つの特徴を全て含んでいるのでnumpyを使うことをお勧めしたいのですが、世の中にはscitoolsという物がありまして、このscitoolsをインストールするとSciPy, ScientificPython, Gnuplotなどのもっと広い意味で研究に使いそうなモジュールが全部ついてくるというお得なものなのです。当然numpyも入っています。なので、僕はscitoolsで一括して全部入れてしまいました。Pythonを使って演算を行うのならscitoolsを入れておけばいいと思います。

numpyで線形代数

ここからはnumpyを使った話。numpyをインポートするには

from numpy import *

scitoolsを使っている場合はこれらでもOK

from scitools.basics import *     # numpy, scipy, scitools.numpytoolsをインポート
from scitools.numpytools import * # 数値計算に便利な関数がたくさん入ってるらしい
from scitools.std import *        # basicsとeasyviz(結果を表示するためのパッケージ)をインポート

以下、numpyがインポートされているとします。

行列作成

numpyにはndarrayクラスがあります。行列のようなものだと考えればいいと思います(ただし、行列演算を行うにはmatrixクラスのインスタンスに変換する必要があります)。ndarrayを作るにはarray関数を使います。

>>> a = array([1,2,3])
>>> type(a)
<type 'numpy.ndarray'>
>>> a
array([1, 2, 3])
>>> print a
[1 2 3]
>>> a = array([1,2,3], float) # 各要素をfloat型にすることもできる
>>> print a
[ 1.  2.  3.]
>>> array([1., 2., 3.]) # これでもOK
array([ 1., 2., 3.])

nuarrayはtolistメソッドを使うとPythonのListオブジェクトが返されます。

>>> a.tolist()
[1.0, 2.0, 3.0]

行列の形を変えるのも簡単です

>>> a = array([[1,2,3],[4,5,6],[7,8,9]])
>>> print a
[[1 2 3]
 [4 5 6]
 [7 8 9]]
>>> a = a.ravel() # ravelを使うと一行になる
>>> a
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a.shape = (9, 1) # 9x1行列に変換
>>> a
array([[1],
       [2],
       [3],
       [4],
       [5],
       [6],
       [7],
       [8],
       [9]])
>>> a = a.reshape(3,3)　# こういう書き方もできる
>>> a 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

単位行列を作るにはeye関数を使います

>>> eye(3)
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

『-1から+1までを0.5刻みで』みたいなndarrayはscitools.numpyutilsにあるseq関数を使えばOK

>>> from scitools.numpyutils import *
>>> a = seq(-1, 1, 0.5) # -1から+1までを0.5刻みで
>>> a
array([-1. , -0.5,  0. ,  0.5,  1. ])

numpyだけを使うならlinspaceを使って似たようなことができます。linspaceは第3引数で範囲の分割数を指定します。

>>> b = linspace(-1, 1, 5) # -1から1の範囲を5分割
>>> b
array([-1. , -0.5,  0. ,  0.5,  1. ])

fromfunction関数を用いれば関数を指定して行列を生成することもできます。例えば

のような3x3行列は

>>> m = fromfunction(lambda i, j: i**2 + j*3, (3,3))
>>> m  
array([[  0.,   3.,   6.],
       [  1.,   4.,   7.],
       [  4.,   7.,  10.]])

行列要素へのアクセス

nbarrayはPythonのListと同じように要素にアクセスできます。

>>> m # このようなmがあったとすると
array([[  0.,   3.,   6.],
       [  1.,   4.,   7.],
       [  4.,   7.,  10.]])
>>> m[1]
array([ 1.,  4.,  7.])
>>> m[0,1] # m[0][1]と同じ
3.0
>>> m[:,2]
array([  6.,   7.,  10.])
>>> m.max()
10.0
>>> m.min()
0.0
>>> m.argmax() # 一行だったときのm.max()のindex。同様にargmin()も存在する
8
>>> m.ravel()[m.argmax()] # m.max()と同値
10.0

その他の操作

最小値(最大値)を設定して、それ未満(それより大きい)要素を指定した最小値(最大値)に丸め込みたい場合はclip関数を使うと簡単にできます

>>> a = arange(0, 10) # arangeはnbarray版のrange関数
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a.clip(min=2, max=7) # 2以下の数は2に、7以上の数は7に丸め込む
array([2, 2, 2, 3, 4, 5, 6, 7, 7, 7])

全ての要素を+1したい場合はそのまま+1します(ここら辺がmatrixオブジェクトとの違い)

>>> m + 1
array([[  1.,   4.,   7.],
       [  2.,   5.,   8.],
       [  5.,   8.,  11.]])

まとめ

今回はnumpyのnbarrayオブジェクトの簡単な使い方について説明しました。
次回はmatrixオブジェクトを使った行列計算などについて(気が向いたら)説明しようと思います。