breezing NumPy-bration - FANCOMI Ad-Tech Blog

　k_oomoriであります！これまで数値計算といえばPythonのNumPy一択かなと思っていたのですが、Scalaでも似たようなライブラリはないものかと調べてみたところ、ありました、breezeが。ここではベクトルや行列の演算について性能を調べてみたいと思います。なお、使用したサーバはAWSのc3.large EC2インスタンス、OSはUbuntu 14.04、ソフトウェアバージョンはPython 2.7.6, NumPy 1.8.2, Scala 2.10.4, breeze 0.10となります。

ベクトルの内積

　まずはベクトルの内積 $\vec{v}\cdot\vec{w}=\sum v_iw_i$ から。NumPyでは例えばランダムな値を持つ3次元ベクトルは

import numpy as np
v = np.random.random(size=(3))
v
array([ 0.8716417 ,  0.01981258,  0.25848169])

で生成できます。自分自身との内積は

np.dot(v,v) # NumPyの関数
v.dot(v) # ndarrayのメソッド
(v * v).sum() # 成分ごとに積を取ってから和を取る
(v ** 2).sum() # 成分ごとに2乗してから和を取る

などといった方法で計算できるので、いろいろな次元のベクトルに対して速度を計測してみます。なお、IPythonの%timeitを使うと計測が捗ります。

%timeit np.dot(v,v)
1000000 loops, best of 3: 672 ns per loop

これは1回あたり数百ナノ秒（1ナノ秒は10億分の1秒）で処理が完了するため、100万回試行してみたら上位3回の平均が672ナノ秒でした！という意味のようです。
　対してbreezeでは

import breeze.linalg._
val v = DenseVector.rand[Double](3)
v dot v //DenseVectorのメソッド
v.t * v //縦ベクトルvを転置して得られた横ベクトルとvとの行列としての積
sum(v:*v) //成分ごとに積を取ってから和を取る

のようになります。では早速結果を。

ベクトルの次元		10	100	1000	10000	100000
NumPy	np.dot(v,v)	557	666	1,560	10,300	98,000
	v.dot(v)	918	1,020	1,920	10,700	98,700
	(v*v).sum()	5,910	6,220	8,140	23,000	171,000
	(v**2).sum()	6,450	6,910	8,790	23,800	172,000
breeze	v dot v	106	188	1,102	10,297	102,300
	v.t * v	116	198	1,104	10,210	101,400
	sum(v:*v)	68	416	4,373	42,110	401,500

時間の単位はナノ秒です。これから読み取れるのは

NumPyでは成分ごとの積を取るよりdotに任せた方が速い
メソッド版よりnp.dotの方が若干速い
2乗を求めるにはv**2よりv*vの方が若干速い
breezeではベクトルの次元が小さい時は成分ごとに掛けて和を取った方が速いが、数十次元以上になる場合はdotや*の方が速い
v dot vとv.t*vはあまり差がなく、次元によって優位性が異なるらしい
NumPy最速対 breeze最速の比較では、次元が10000くらいまでではbreezeが速いが、それ以上になるとNumPyが逆転する

行列の積

　行列の生成とその積を求める方法は、NumPyでは

m = np.mat(np.random.random(size=(5, 3)))
m.T*m

breezeでは

val m = DenseMatrix.rand(5, 3)
m.t*m

とします。（転置を取っているのは、そうしないと非正方行列の積が定義できないからです）

行列サイズ	10x10	100x100	1000x1000	60000x500
NumPy	5	893	752,000	19,500,000
breeze	2	901	743,600	15,800,000

時間単位はマイクロ秒（ $\mu$ s, 100万分の1秒）です。サイズが大きくなったときにbreezeが少し速いかなといった程度で、それほど大きな差は見られません。60000x500という行列のサイズは、たまたま手元にあったデータがそのくらいの大きさでこの検証を始めたことに由来します。

行列の要素ごとの積

　2つの $n\times m$ 行列 $S,T$ に対し要素ごとの積 $S_{ij}T_{ij}\ (1\le i\le n, 1\le j\le m)$ を考えます。こんな積を使うことがあるのか？と思うかもしれませんが、非負値行列因子分解で使うのです。2つの行列 $m,n$ を受け取り、要素ごとの積を取った結果を行列として返すためには、NumPyでは

np.mat(np.array(m) * np.array(n))

breezeでは

m:*n

とします。

行列サイズ	10x10	100x100	1000x1000	60000x500
NumPy	7.44	29	35,300	988,000
breeze	0.43	26	2,800	114,000

単位はマイクロ秒です。np.arrayへの変換等がない分有利なのか、こちらはbreezeの方がかなり速いですね。

行列ノルム（の2乗）

　最後に行列の各成分の2乗和であるフロベニウスノルムを考えます。
$\displaystyle ||M||_F^2=\sum_{i=1}^n\sum_{j=1}^m M_{ij}^2 \tag{1}$
これを2つの行列の差に対して適用すれば、それらがどれだけ離れているかの基準となります。breezeではシンプルにsum(m:*m)として計算できますが、内積計算では成分ごとに計算するよりdotで計算した方が早かったことを思い出すと、行列 $M$ を $n$ 次元の縦ベクトルが横に $m$ 個並んだものとして計算してみたくなります。
$\displaystyle M=(\vec{v}_1 \ \cdots \ \vec{v}_m), \quad ||M||_F^2 = \sum_{j=1}^m \vec{v}_j\cdot \vec{v}_j \tag{2}$
これを実現するScalaコードは

(for (i <- 0 until m.cols) yield {val v = m(::, i); v dot v}).sum

と書けます。逆に $m$ 次元横ベクトルが縦に $n$ 個並んだものとみなすこともでき、その時のコードは

(for (i <- 0 until m.rows) yield {val v = m.t(::, i); v dot v}).sum

となります。
　NumPyの成分版はma = np.array(m); (ma*ma).sum()、ベクトルのdot積版はScalaのように一息で書く方法が思いつかなかったので

def column_vectors(m):
  norm = 0
  for i in xrange(m.shape[1]):
    v = np.array(m[:, i])
    norm += np.dot(v.T,v)[0][0]
  return norm

def row_vectors(m):
  norm = 0
  for i in xrange(m.shape[0]):
    v = np.array(m[i, :])
    norm += np.dot(v,v.T)[0][0]
  return norm

としました。結果は

	行列サイズ	3x2	100x10	1000x100	60000x500
NumPy	ma = np.array(m); (ma*ma).sum()	8	10	378	631,000
	column_vectors(m)	27	128	1,730	360,000
	row_vectors(m)	40	1,240	12,700	838,000
breeze	sum(m:*m)	1	4	369	143,593
	(for (i <- 0 until m.cols) yield {val v = m(::, i); v dot v}).sum	1	3	126	36,900
	(for (i <- 0 until m.rows) yield {val v = m.t(::, i); v dot v}).sum	1	16	284	177,360

単位はマイクロ秒です。これは圧倒的にbreezeの方が速いですね。また、あえて縦に長い行列を用いて実験したのですが、その場合は縦ベクトル、つまり長い方を一つのベクトルとみなしてベクトル演算(dot)した方が効率的であることがわかりました。

結論

私はコレで、Pythonやめました。（昭和を生きた人向けのネタです）

おまけ：インストール方法

breeze→ こちらを参照してください。
NumPy→ python-numpyをyumやapt-getでinstallしてください。RPMパッケージを作成したい方はこちら。