対象

Python及びNumPy初心者に向けて書いています. 「C言語は使えるけど最近Pythonを始めた」とか「Pythonらしい書き方がよくわからない」に該当する物理系の数値計算を目的とした方には特に有用かもしれません.

また, 自分の不勉強のために間違った記述があるかもしれません. ご容赦ください.

あらまし

Pythonは数値計算系のライブラリが非常に充実しており, かつ使いやすくまとまっています. その多くはFortranで書かれた遺産のラッパーだったりするのですが, C/C++などから呼ぶよりもとても簡単です. 例えばC/C++からLAPACK¹を呼んで固有値問題を計算しようとすると

info = LAPACKE_zheevr(LAPACK_COL_MAJOR, 'V', 'A', 'L', H, (lapack_complex_double*)a, lda, vl, vu, il, iu, abstol, &mm, w, (lapack_complex_double*)z, ldz, isuppz);

のように引数が多すぎて煩雑です. 引数が多いほど計算の自由度が増すとはいえ, あまり初心者向けには見えません. また, C/C++にはネイティブでcomplexに対応していないことから, ラッパーの内部でlapack_complex_double等の型がdefineされており, 混乱を招きやすいです. 一方Pythonでは非常に簡潔です:

w, v = numpy.linalg.eig(H)

LAPACKにあるような引数もちゃんと用意されていますが², デフォルト値が設定されているので, 細かいことを考えずにただ解きたいだけならコレだけでいいのです. シンプルでいいですね. ちなみに各種ソルバーはNumPyよりもSciPyのもののほうが高性能なようです.

こうなるといろいろな数値計算をPythonに任せたくなってくるのですが, この子にも欠点はあります. forループが死ぬほど遅いのです. いわゆる行列演算では3重程度のループは頻出ですし, 数値計算というのは「細かく差分化してループ処理」が基本なので, なかなかに致命的です³. コンパイル言語ではないのでしょうがないとはいえ, 他のLL言語よりも遅い気がします. 行列演算のために

for i in range(N):
    for j in range(N):
        for k in range(N):
            c[i][j] = a[i][k] * b[k][j]

なんて書いたらもうおしまいです.

NumPyの組み込み関数と演算

ではどうするかというと, NumPyの組み込み関数(universal function)を活用します. NumPyの組み込み関数は, NumPy配列ndarrayの各要素に対して作用し, 新しいndarrayを返してくれます. これがとても重要でして, テンソル演算におけるforループを駆逐する可能性を秘めています. 例えばmath.sin(cmath.sin)はスカラーしか渡せませんが, numpy.sinは

>>> theta = numpy.array([numpy.pi/6, numpy.pi/4, numpy.pi])
>>> numpy.sin(theta) 
array([  5.00000000e-01,   7.07106781e-01,   1.22464680e-16])

のように, NumPy配列を投げると各要素にnumpy.sinを作用させてくれます. そして往々にしてこの手の処理はC/C++と遜色ないほど高速です. NumPyの主な実装はCやFortranであり, かつ線形演算ではBLAS⁴/LAPACKが頑張ってくれているからです.

また, ndarrayの定義する四則演算も同様に高速かつ多様です.

その1 : ベクトル

ベクトルに関する線形演算をforループ抜きで書いていきましょう. 内積・外積等の関数がちゃんと用意されているのは想定内なので, 以下では主にndarrayに対する和・積・商などがどのように定義されているのかを見ていきます. ndarrayに対する演算は, Python組み込みのlistに対する演算とは全く異なります
以後import numpy as npのようにエイリアスを設定しています.

ベクトルとスカラー

1次元のndarrayは一種のベクトルとみなしてもよいでしょう.
スカラー(complex, float, int)はndarrayの全ての要素に対して作用します:

>>> a = np.array([1, 2, 3])
>>> a + 2
array([3, 4, 5])
>>> a * (3 + 4j)
array([ 3. +4.j,  6. +8.j,  9.+12.j])

ベクトル同士

ndarrayはインデックスの一致した要素同士が作用します:

>>> a = np.array([1, 2, 3])
>>> b = -np.array([1, 2, 3])
>>> a + b
array([0, 0, 0])
>>> a * b
array([-1, -4, -9])
>>> a / b
array([-1., -1., -1.])

内積でも外積でもないところが面白いですね. 割り算まで定義できます. ある差分化した関数に微分演算を施す場合などによく使います. 次元の異なるndarray同士の和・積は定義されません.

reshapeしたndarray同士

ndarrayにはreshapeと呼ばれるメソッドがあり, 例えば3×1のベクトルを1×3に組み替えるようなことなどができます. 縦ベクトルに対する横ベクトルのようなものですが, 積(*)の定義がベクトル空間のそれとは異なるので単純な対応はしていません. このreshapeしたベクトルとの積はとても面白いです:

>>> c = a.reshape(3, 1)
>>> c
array([[1],
       [2],
       [3]])
>>> c + b
array([[ 0, -1, -2],
       [ 1,  0, -1],
       [ 2,  1,  0]])
>>> c * b
array([[-1, -2, -3],
       [-2, -4, -6],
       [-3, -6, -9]])

(3×1) + or * (1×3) = (3×3)のような規則を持ちつつ, しかし先の演算と同様, この演算は可換です. 初見では不思議に思うかもしれませんが, ちょっと眺めれば前節の規則を用いて説明が可能であることがわかると思います. これを利用すると行列の初期化すらforループを使わずに済む可能性を持っています. 巨大な行列をパラメータを変えるごとに都度初期化して計算を進めるようなタスクでは, 大きな威力を発揮します.

これらの型に縛られない非常に柔軟な仕様は個人的にとても気に入っています⁵. 他の言語よりも直感的な記述ができるように思います.

ベクトルの初期化

演算ではありませんが, そもそもベクトルを用意するのにforループを使っていては芸がありません. ndarrayをつくる関数は幾つもありますが, よく使うものの一部を挙げます:

>>> L = 1
>>> N = 10
>>> np.linspace(-L/2, L/2, N)
array([-0.5       , -0.38888889, -0.27777778, -0.16666667, -0.05555556,
        0.05555556,  0.16666667,  0.27777778,  0.38888889,  0.5       ])

>>> dL = 0.2
>>> np.arange(-L/2, L/2, dL)
array([-0.5, -0.3, -0.1,  0.1,  0.3])

>>> np.logspace(0, 1 ,10, base=np.e)
array([ 1.        ,  1.11751907,  1.24884887,  1.39561243,  1.5596235 ,
        1.742909  ,  1.94773404,  2.17662993,  2.43242545,  2.71828183])

スカラーやベクトルとの演算を組み合わせればforループがなくともなんとかなりそうですね.

その2 : 行列

ベクトルと話の流れはほぼ同じです. 行列積には専用の関数が用意されているので以下も演算の話です. 演算結果はもうだいたい予想がつくかと思います.

行列とスカラー

>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> a
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> a + (2 - 3j)
array([[  3.-3.j,   4.-3.j,   5.-3.j],
       [  6.-3.j,   7.-3.j,   8.-3.j],
       [  9.-3.j,  10.-3.j,  11.-3.j]])
>>> a * 0.2
array([[ 0.2,  0.4,  0.6],
       [ 0.8,  1. ,  1.2],
       [ 1.4,  1.6,  1.8]])

行列とベクトル

>>> b = np.array([1, 2, 3])
>>> a + b
array([[ 2,  4,  6],
       [ 5,  7,  9],
       [ 8, 10, 12]])
>>> a * b
array([[ 1,  4,  9],
       [ 4, 10, 18],
       [ 7, 16, 27]])
>>> a / b
array([[ 1. ,  1. ,  1. ],
       [ 4. ,  2.5,  2. ],
       [ 7. ,  4. ,  3. ]])

行列と行列

>>> b = np.array([[-1, 2, -3], [4, -5, 6], [-7, 8, -9]])
>>> a + b
array([[ 0,  4,  0],
       [ 8,  0, 12],
       [ 0, 16,  0]])
>>> a * b
array([[ -1,   4,  -9],
       [ 16, -25,  36],
       [-49,  64, -81]])

行列の初期化

行列用にも便利な関数が用意されています. 以下はほんの一部です:

>>> np.identity(4)
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

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

>>> c = np.array([1, 2, 3, 4])
>>> np.diag(c)
array([[1, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 3, 0],
       [0, 0, 0, 4]])

おわりに

NumPyの便利な関数や細かい演算の規則などは他の方々がもっと良い記事を書いてくれていると思います. この記事見て, NumPyらしい書き方をすればforループを使わずに様々な計算ができそうであることを感じ取っていただければ幸いです. これらを活用すれば, 理系の大学生がプログラミングの授業で学ぶような数値計算アルゴリズムはC/C++に負けない速度で動いてくれます. もっとも, その手のアルゴリズムは既にSciPyなどに用意されているものも多いです.

少し長くなってしまったので, 具体的な数値計算への応用については別の記事にまとめたいと思います. ありがとうございました.

追記(2016/11/25): 以下に続きます
NumPy・SciPyを用いた数値計算の高速化 : 応用その1
NumPy・SciPyを用いた数値計算の高速化 : 応用その2

LAPACK(Linear Algebra PACKage)はFortran90で書かれた線形代数演算のライブラリです. C/C++用のラッパーではLapackeと呼ばれるものがあります. MKL互換のインターフェイスを持っているので, MKLのマニュアルが参考になります. ↩
詳しくはNumPyやSciPyのマニュアルをどうぞ. ↩
それ以外にも, 関数呼び出しのオーバーヘッドも大きい印象です. プロコンだと「動的計画法で書くと通るのに, メモ化再帰ではTLE」といったことがよく起こります. ↩
LAPACKよりも基本的な線形演算セットを含んだライブラリ. LAPACKの内部ではBLASを呼んでいます. ↩
「型がない」ということを嫌う人もいます. 「ある仕様」を満たすような設計をするプログラムと異なり, 物理における数値計算は出力が正しいかどうかを確認することがそう簡単ではありません. そしてプログラムの実行中に「変数の型が非明示的に変わる」ことが許されるとさらにデバッグは難しくなります. 可観測量はfloatに決まってるので, 静的型付けのほうがバグは生まれにくいという主張は十分に理解できます. でも自分はPythonが好きです. ↩

NumPyによる数値計算の高速化 : 基礎

対象