読者です 読者をやめる 読者になる 読者になる

F@N Ad-Tech Blog

株式会社ファンコミュニケーションズ nend・nex8のエンジニア・技術ブログ

breezing NumPy-bration

 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

時間単位はマイクロ秒( \mus, 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パッケージを作成したい方はこちら