4.numpyを使いこなそう
Numpy
はNumerical Pythonの略称であり,Pythonにおいて数的データ処理を行うことに特化したライブラリです.Python=Numpyと言っても過言ではないぐらいPythonの可能性を広げてくれる究極のライブラリの一つです.本書の主題もNumpyに置いています.ここでは,Numpyの概念からインストール,多用する関数,データの入出力,関数を活用した実務的なデータ処理といった内容までを一気に紹介します.
4-1.Numpyはなぜ早い?
さて,Pythonを用いる上でしばしば批判されるのが実行速度です.確かに動的型付け言語なので,配列の各要素を参照する際に型の認定が必要となるため速度が遅くなるのは事実です.しかし,Python=遅いとするのはステレオタイプで不勉強な人の発言だと考えて下さい.Numpyを使えば他のプログラミング言語を凌駕する速度を得られる可能性すらあります.その理由は以下のとおりです.
- Numpyは元から配列の型(Int,float)が決まっている
- ブロードキャスティングによる配列内部の一括処理が可能
- Numpyはその中身がコンパイル言語であるC(Cython)で書かれている
- 関数内部の記述も論理的で無駄がない
- 様々な計算がすでに関数化されているのでコードを書く時間も短くなる(c.f.車輪の大発明)
- 正しい手順で導入した場合計算環境が自分のPCに最適化される(例:Intel MKL)
等々,様々な理由があります.ここでは,(1)と(2)の説明を順次行います.
(1)Numpyは元から配列の型(Int,float)が決まっている
コンパイル言語との違いでPythonは動的型付け言語であると説明しました.動的な型付けとは,FortranやC言語のように変数にintやfloatなどの型付けが必要ないというものです.しかし,Python内部での実際の計算では各数値の型を判別した上で計算が行われています.ここが,Pythonが遅いと言われる所以です.何事も経験が一番なので以下のコードを実行してみましょう.私のPC環境では,Python-native codeの実行時間は7.8s,Numpy_codeの実行時間は1.9sでした.つまり,pythonだけで計算するのとNumpyを用いて計算するのとでは4倍の差があります.配列サイズが大きくなる程にこの差は大きくなります.侮ることなかれ,4倍は小さく見えますが,1回のシミュレーションに4日必要なプログラムが1日で終了できると考えれば非常に大きな改善です.ただし,配列のサイズが小さい場合はリストを用いたほうが速いケースもあります.
コンパイル言語との違いでPythonは動的型付け言語であると説明しました.動的な型付けとは,FortranやC言語のように変数にintやfloatなどの型付けが必要ないというものです.しかし,Python内部での実際の計算では各数値の型を判別した上で計算が行われています.ここが,Pythonが遅いと言われる所以です.何事も経験が一番なので以下のコードを実行してみましょう.私のPC環境では,Python-native codeの実行時間は7.8s,Numpy_codeの実行時間は1.9sでした.つまり,pythonだけで計算するのとNumpyを用いて計算するのとでは4倍の差があります.配列サイズが大きくなる程にこの差は大きくなります.侮ることなかれ,4倍は小さく見えますが,1回のシミュレーションに4日必要なプログラムが1日で終了できると考えれば非常に大きな改善です.ただし,配列のサイズが小さい場合はリストを用いたほうが速いケースもあります.
Python-native code
- lstA = [5]*500000000
- lstB = [10]*500000000
- lstC = []
-
- for i in range(len(lstA)):
- X = lstA[i] + lstB[i]
- lstC.append(X)
- print(lstC)
-
[15,15,15・・・]
注意:メモリサイズの小さいPCではフリーズの原因になる.具体的にはサイズを5000にするなどの対応策がある.
Numpy_code
- import numpy as np
- a = np.full(500000000,(5))
- b = np.full(500000000,(10))
- c = np.empty(500000000)
-
- for i in range(len(a)):
- c[i] = a[i] + b[i]
- print(c)
- print(a.dtype,b.dtype,c.dtype)
-
[15. 15. 15. ... 15. 15. 15.] int32 int32 float64
np.full(shape, (fill_value))
指定した fill_valueで満たされた大きさがShapeの配列を生成します.
指定した fill_valueで満たされた大きさがShapeの配列を生成します.
np.empty(shape)
ランダムな値で満たされた大きさがShapeの配列を生成します.
ランダムな値で満たされた大きさがShapeの配列を生成します.
np.darray.dtype
指定したnp.darrayの型を返します.なお,darrayのdは次元数(Dimension)を表す.
指定したnp.darrayの型を返します.なお,darrayのdは次元数(Dimension)を表す.
先程のプログラムでは,長さの同じ配列lstAとlstBの要素を一つずつ足し合わせ,その結果であるlstCを作成するというものでした.今回の配列長は5億ですから,5億回足し算をしていることになります.その度にlstA[i]とlstB[i]の型を検索・同定する処理が行われます.ここに時間がかかっています.
一方,NumpyはNumpy codeの実行結果から分かるように,配列darray全体に対して型が決められています.これにより,a[i]とb[i]における型の検索・同定が不要になります.余計なプロセスが減った結果,同じ計算でも速度が早くなるのです.当然,(3)にあったC言語で実施されているということも速度が早い理由でもあります.
(2)ブロードキャスティングによる配列内部の一括処理が可能
Pythonの有名な機能の一つはブロードキャスティング(Broadcasting)です.ブロードキャスティングは放送という意味がメインですが,広げるという意味もあります.放送というのも情報を拡散するという意味からの派生で生まれました.Python-nativeなリスト[1,2,3]を2倍すると,[1,2,3,1,2,3]となります.普通は[1,2,3]×2を行うと,[2,4,6]になるスカラー倍をイメージするのが普通ではないでしょうか?Numpyはこれが出来ます.スカラーベースなので[1,2,3]+2=[3,4,5]や[2,4,6]/2=[1,2,3]という計算が出来ます。
Pythonの有名な機能の一つはブロードキャスティング(Broadcasting)です.ブロードキャスティングは放送という意味がメインですが,広げるという意味もあります.放送というのも情報を拡散するという意味からの派生で生まれました.Python-nativeなリスト[1,2,3]を2倍すると,[1,2,3,1,2,3]となります.普通は[1,2,3]×2を行うと,[2,4,6]になるスカラー倍をイメージするのが普通ではないでしょうか?Numpyはこれが出来ます.スカラーベースなので[1,2,3]+2=[3,4,5]や[2,4,6]/2=[1,2,3]という計算が出来ます。
Broadcasting.py
- import numpy as np
- a = [1,2,3]
- b = np.full((5,3),(2))
- print(b)
- print(a+b)
-
[[2 2 2] [2 2 2] [2 2 2] [2 2 2] [2 2 2]]
[[3 4 5] [3 4 5] [3 4 5] [3 4 5] [3 4 5]]
しかし,これはNumpyがスカラーを扱えるだけであって,ブロードキャスティングではありません.ブロードキャスティングは以下のように自動的にndarrayを拡張する機能です.Broadcasting.pyの5行目でa+bが行われています.しかし,サイズの異なる配列同士を加減する概念はスカラーにはありません.加減を行うには,配列のサイズを一致させ,行列として対応する部分の加減を行う必要があります.Numpyは自動でそれをしてくれています.ブロードキャスティングの概念は以下のように定式化出来ます.
ブロードキャスティングは異なるサイズの配列の小さい方を大きい方に合わせるように拡大し,加減を可能にします.異なるサイズの配列同士の乗除も同様にブロードキャスティングが適用されます.以下にその例を定式化します.ここで出てくる記号はアダマール除算を意味します.
上記
は割り算へ適用されるものですが,
はアダマール積を意味します.これはPoint-wise calculattionとも呼ばれます.すなわち,2つの行列を比較して同じ位置にあるもの同士を乗除するものです.
また内部計算は行列ではなく,要素ごとの積であるアダマール積が実行されます.配列の位置的に対応する部分が乗除されるので非常に感覚的な計算が出来ます.
Numpyの行列計算は基本的にアダマール演算に従っていると考えましょう.内積をしたければnp.dot(a,b)で可能です.1行N列のベクトルであればnp.dot(a,b) = np.sum(a*b)となります.
Numpy_code_2
Numpy_code_2
- import numpy as np
- a = np.full(500000000,(5))
- b = np.full(500000000,(10))
- c = np.empty(500000000)
-
- for i in range(len(a)):
- c[i] = a[i] + b[i]
- print(c)
-
[15. 15. 15. ... 15. 15. 15.]
Numpy_code_3
- import numpy as np
- a = np.full(500000000,(5))
-
- c = a + 10
- print(c)
-
[15. 15. 15. ... 15. 15. 15.]
Numpy_codeを少し書き換えたNumpy_code_2とNumpy_code_3を例にブロードキャスティングの有無を比較してみます.Numpy_code_2はNumpy_codeと同様に2.0sが必要でしたが,Numpy_code_3は僅か0.2sで終了しました.計算の内容自体はどちらも変わりませんが,速度は10倍です.従って,Python-nativeなコードと比較して40倍の改善になります.さらに,Numpy_code_3は変数の宣言も少なく,スカラーを活用してより感覚的に計算をプログラムできています.
Numpy_code_3では,5億列の配列aに合わせて10がブロードキャストされ,a+10の行列計算が行われています.コンパイル言語では,C = A + Bが忠実に行えるようにAとBのサイズが一致するように予め変数を整えておく必要があります.また,A+Bの結果であるCもAとBのサイズに合うようにint C[500000000]のような宣言を予めしておく必要があります.
では,なぜNumpy_code_3のようにブロードキャスティングを活用して書くと速度が早くなるかといえば,それは配列の拡大や計算結果である配列cの生成をNumpyの内部,つまりはC言語を通して行ってくれるからです.Numpy_code_2では,Numpyで生成した配列の要素を一つ一つ抜き出して計算を行い,最後に予め生成した同じサイズの配列に埋め込むという作業をすべてPython上で行っています.一方,Numpy_code_3では,事前の配列の作成は一つだけである上,計算にfor文が出ていません.しかし,しっかりと計算結果は出力されています.a+10という計算に必要な配列サイズの拡大から行列計算,計算結果の変数出力までをすべてNumpy上で行ってくれているのです.Numpy上で行っている足し算はNumpy_code_2のfor i in range(len(a)): c[i] = a[i] + b[i]と変わりありませんが,これがC言語で書かれているため同じ計算でも速度が向上します.つまり,同じ計算でもNumpyを通すことで速度が顕著に改善されるのです.
まずはExcelができるようになろう]では,for文による計算は速度低下の原因になることを述べました.これは,Python上でfor文を使うような配列計算はNumpyによって書き換えられる可能性が高く,Numpyを通すことでfor文を使うことに変わりはなくとも,それがC言語で書かれているために計算速度が向上するためです.さらに,Numpyを活用することでスカラーの概念に沿ったより感覚的な記述やより少ない変数宣言が可能になり,コード記述速度の改善や可読性の向上にもなります.
コラム5.Ipythonによる実行速度の比較
コマンドプロンプトでpythonを入力すると簡易版pythonが出てくることは学びました.この簡易版pythonをインタラクティブモード呼びます.これを強化したものとしてIpython(Jupyter Notebook)があります.Anacondaなどでパッケージインストールした人はすでにIpythonが導入されているかと思いますが,導入されていない人はpip install ipythonでインストールしておきます.Ipythonは便利な機能が沢山ありますが,今回はその一つである時間計測機能を活用してNumpyの関数実行速度を比較します.Ipythonによるインタラクティブモードを使用するには,コマンドプロンプト上でipythonを入力し,実行します.すると,In [1]:が表示されるはずです.これが出ていれば成功です.
ここからはIpythonを使って関数の実行速度を比較していきます.比較するのはnp.zeros,np.ones,np.full,np.empty,Python-nativeなリストの5種類です.Ipythonで関数の実行速度を測るには関数の前に%timeitを付与するだけです.これはIpythonのマジックコマンドと呼ばれるものの一つです.例えば,np.zerosの実行速度を測る場合,import numpy as npをIn [1]で実行した上で,In [2]にて %timeit np.zeros(1000)を実行します.すると,次の行に複数回実行した時に得られる平均実行速度±標準偏差が出力されます.
3.38 µs ± 1.19 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)
これを同様にnp.ones,np.full,np.emptyについても行います.私のPCにおける計測結果は以下の通りとなりました.
np.zeros(1000) 3.38 µs ± 1.19 µs np.ones(1000) 9.91 µs ± 2.17 µs np.full(1000,(5)) 8.84 µs ± 2.14µs np.empty(1000) 2.03 µs ± 277 ns [5]*1000 6.16 µs ± 1.87 µs
ということで,結果は
np.empty > np.zeros >> Python-native list >> np.ful l> np.ones
となり,np.emptyの高速さが際立つ結果となりました.あるサイズの配列を作る際,pythonはPCのメモリースペースをはじめに確保します.その後,用途に合わせて配列の中身の数値を設定していきます.zerosであれば,そのメモリースペースをすべてゼロに.fullであれば,すべてを指定した数値に変換します.np.emptyはメモリースペースを確保しておしまいです.確保したメモリースペースにおける数値はその時々によって変化します.Ipythonでprint(np.empty(5))の実行結果を見てみましょう.
[4.47685025e-290 4.47802457e-290 4.47661539e-290 4.47739827e-290 4.47747656e-290]
print(np.empty(10))はどうでしょうか.
[4.00276266e-289 1.77922788e-290 4.70809539e-292 4.05614237e-307 4.05646827e-307 4.05679421e-307 4.31378456e-289 4.29502680e-289 4.30041301e-289 4.35994307e-289]
このようにランダム数値が並んでいることが分かるかと思います.しかし,このように確保した配列がNumpy codeのように別の計算結果に書き換えられるのであれば,確保した時の数値は関係がありません.このようなケースではnp.emptyによって配列を確保することがより合理的と言えるでしょう.
この便利な%timeitマジックは関数だけでなく%%timeitとすることでforループ内の速度計測にも使用できます.
この便利な%timeitマジックは関数だけでなく%%timeitとすることでforループ内の速度計測にも使用できます.