scipy/numpyのBoost.Pythonによる高速化

Boost.PythonとPyUblasに関して, 今まで適当に書いていた日記をまとめたもの. 詳細は過去の日記へのリンクをあちこちに貼ってある.

前書き

Pythonは言わずと知れた非常に遅い言語だ. Pythonで大きなループを含む処理を書くことは即ち死を意味する. 例えば手元のthinkpad Core2Duo2.5GHzでは, Pythonで2つの500×500の行列同士の行列積を計算するのに29秒かかる. これは, C++では僅か0.6秒で終了することを考えると非常に遅い. (ただし, 乱数による行列の生成処理を含み, コンパイルオプションは-O2である)

しかし世間は広く, この遅いPythonで数値解析をしたい奇特な人間がいるらしい. そこで恐らく今もっとも使われているライブラリがscipy(numpy)である.

このscipyには様々な数値解析用の関数などが含まれている. 例えば, 先ほどの500×500の行列積の問題はわずか0.4秒で終了する. これはC++よりも速い! 先ほど試したC++の行列積のプログラムは私が適当に書いたものなので, scipyできちんとチューニングしたものの方が速いだけなのだが.

ともかく, この心強いライブラリであるところのscipyを使えば典型的なユースケース(行列積とか総和とか)では非常に役に立つ. ただ, 自分が行いたい処理がscipyに含まれていない場合はどうすればいいのか?

scipyはarrayという配列・行列用のクラスが含まれている. これは恐らく数値の型と配列サイズを最初に指定することで, データをメモリの連続した場所に配置し高速に処理するためのデータ型である.

これを使ってPythonで処理を書けばよくね? とすると, それは即ち死を意味する. 驚くべきことにarrayの各要素へのアクセスはPythonのリストよりも遅い. ある長さ1000の配列targetの要素を全て足し合わせる処理を1万回行なった場合,

result = 0
for x in target:
    result += x

targetがlistの場合は0.4秒なのに対し, arrayの場合は4.5秒もかかる(もちろんscipyのsum関数を使えば話は別だが). このように, scipyのarrayはループ処理には全く適さない.

この, 「scipyを基本的に使いたいけど, 一部の処理が遅すぎる!」という状況に対して世間では色々な解を用意しているらしい. ただ, 基本的には動的ライブラリ(.so)の形で実装されたPythonの拡張モジュールを作成して, ループを用いる重い処理はそちらで行い, その結果をPythonに返すことで解決するようだ.

拡張モジュールを簡単に作るためにいくつかの手段が提案されている

  • Cython : Pythonに変な文法を混ぜた奴. 簡単に書けるとか言ってるけど覚えたくない
  • SWIG : 昔からあるCのコードをスクリプト言語から呼べるようにする奴. めんどくさそうだけど, C側のコードをあんまり汚さないイメージ
  • Boost.Python : 神

これらと違ってJITコンパイラーを使って速くしようとしているPyPyというプロジェクトもあるが2011年時点ではscipyに対応してないぽい(今はどうなんだろう).

さて前置きが長くなってしまったが, 今回はBoost.Pythonを用いたscipyの高速化を行なった. 簡単なBoost.Pythonの使い方はBoost.Python の機能をざっと紹介してみるがわかりやすい.

題材

実装したのは, Dynamic Time Warpingと呼ばれるベクトルx,y の間の距離を計測する関数である. 例えば, それぞれの要素がある時間の音の大きさを表す2つの音声信号
x = [12 24 6 17 40 40]
y = [12 12 24 6 17 40]
を比較する. この場合, xの1秒目の音の大きさは12, 2秒目の音の大きさは24となる.

xとyの間の距離は通常のユークリッド距離を用いるとEUC(x,y)=33.43となるが, よく見るとこの2つのベクトルは時間的に1秒ずれているだけである. このような時間的なずれを考慮した距離の一種がDTWであり, 詳細は省くがDTW(x,y)=0となる.

この距離関数はscipyに実装されておらず, 二重ループを用いる必要がある. 以下はPythonのlistによる実装である.

def dtw(x,y):
    n = len(x)
    m = len(y)

    cost = [[0 for j in xrange(m)] for i in xrange(n)]

    cost[0][0] = abs(x[0] - y[0])
    for i in xrange(1,n):
        cost[i][0] = cost[i-1][0] + abs(x[i]-y[0])
    for i in xrange(1,m):
        cost[0][i] = cost[0][i-1] + abs(x[0]-y[i])

    for i in xrange(1,n):
        for j in xrange(1,m):
            cost[i][j] = min(cost[i-1][j],cost[i][j-1],cost[i-1][j-1]) + abs(x[i] - y[j])

    return cost[n-1][m-1]

これをBoost.Pythonを用いて高速化する. ちなみにこのcostを2次元のscipy.arrayにすると3〜4倍遅くなる.

仕様

入力として二つのscipy.array x,yを受け取り, 出力としてその二つの間の距離DTW(x,y)を返したい. ここで面倒なのは入力としてscipy.arrayを取る点である. Boost.Python単体ではscipy.arrayを素直に受け取ると不便な形で取得してしまう(Boost.Pythonの使い方とNumpy/Scipy arrayの渡し方).

そこで, PyUblasというライブラリを用いる. これは乱暴に言って, 1次元のscipy.arrayをC++側でBoost::numeric::ublas::vectorとして扱えるようにするライブラリである. なお2次元のscipy.arrayはublas::matrixとして扱う. PyUblasを使う

PyUblasではBoost::numeric::ublas::vectorの代わりに, pyublas::numpy_vectorをarrayへのインターフェイスとして使うことも出来る. これはnumpy.arrayが持つメモリ空間をC++側で共有する. すなわち渡されたnumpy.arrayの一部の値を書き換えることが出来るし, 無駄なメモリコピーも発生しない. これは, ublas::vectorではメモリコピーを行い, 値を書き換えられないこととは対照的である.
しかし, 後に示すが残念ながらこのクラスはarrayの要素へのアクセスが遅い. そのため, 要素へのアクセスが多数発生するようなケースではublas::vectorを用いた方がよい.

実装

細かい仕様はここまでに紹介したリンク先を見てもらうことにして実際のコードを示そう.
C++側 subset_extra_distance.cc

#include <cmath>
#include <boost/numeric/ublas/vector.hpp>
#include <pyublas/numpy.hpp>

using namespace std;

static inline double euc(double x,double y){
  //return sqrt((x-y)*(x-y));
  return std::abs(x-y);
}

static double dtw(const boost::numeric::ublas::vector<double>& x,
                  const boost::numeric::ublas::vector<double>& y){
  using namespace pyublas;
  int n = x.size(),m = y.size();

  //0は番兵
  double costs[n][m];

  costs[0][0] = euc(x[0],y[0]);
  for(uint i = 1;i < x.size();i++){
    costs[i][0] = costs[i-1][0] + euc(x[i],y[0]);
  }

  for(uint i = 1;i < y.size();i++){
    costs[0][i] = costs[0][i-1] + euc(x[0],y[i]);
  }

  for(int i = 1;i < n;i++){
    for(int j = 1;j < m;j++){
      costs[i][j] = min(min(costs[i-1][j],costs[i-1][j-1]),costs[i][j-1]) + euc(x[i],y[j]);
    }
  }

  return costs[n-1][m-1];
}

//BOOST_PYTHON_MODULE()の引数がmodule名
BOOST_PYTHON_MODULE(subset_extra_distance){
  //import_array();
  boost::python::def("dtw",dtw);
}

Python

from scipy import *

import subset_extra_distance
import pyublas

x = array([12,24,6,17,40,40],dtype=float)
y = array([12,12,24,6,17,40],dtype=float)

print subset_extra_distance.dtw(x,y)

このように簡単な記述で拡張モジュールを作成することが出来, 呼び出し側もimport pyublasという記述以外は普通のモジュールのように使うことが出来る.
C++STLvectorに慣れた人なら違和感無く使うことが出来るだろう.

パフォーマンス計測

さて, では拡張モジュールは手間に見合うだけのリターンがあるのか確かめよう. 次の4つのDTWの実装を比較する.

  • PurePython : Pythonのリストで実装
  • ublas_vector : PyUblasのublas::vectorを使用
  • numpy_vector : PyUblasでnumpy_vectorを使用
  • mlpy : mlpyという機械学習ライブラリに含まれるdtw_std関数. DTW関数はCythonを使って実装されているっぽい.

与えるデータは, 2つの長さ16のベクトルを25000回それぞれの手法で比較した. 結果は次のようになった.

手法 時間
PurePython 51.07(s)
ublas_vector 0.80(s)
numpy_vector 2.55(s)
mlpy 2.60(s)

このようにublas::vectorを用いることでPythonと比較して60倍程度高速化した. また前述の通りnumpy_vectorはメモリコピーが生じないにも関わらずublas_vectorよりも遅い.

mlpyよりも3倍程度, ublas_vectorが速い. これはCythonに対するPyUblasの優位を即座に示すものでは無いが, 既存のものより速いのは嬉しい.

余談だがPurePythonにおいて, 計算に用いるテーブルcostをlistでは無くscipy.arrayを用いると16秒程度かかってしまう. 本当にscipy.arrayの要素アクセスは遅い.

また, ublas::vectorは実用上は更に高速化が可能である. mlpyのDTWは二つのarrayを比較するインターフェイスしか持っていないが, 実用上はarrayの二つの集合X, Yの全ての組み合わせの距離を求めたいことが多い. このときX, Yの要素を一つずつ渡すのでは無く, 一度に渡すことでX, Yの要素数が500個ずつだったときに, C++<=>Python間の型変換が50000回から1000回に減少する.
実際に, このように一度にX,Yを渡すようにするとublas::vector版のDTWは1.5倍程度速くなる.

まとめ

ここまでで示したように, 拡張モジュールを作成することは大幅なパフォーマンスの向上をもたらす.
しかもBoost.Pythonを用いることで, 新しい言語を覚える必要もなく, またscipy.arrayとのシームレスな結合も可能である.

また, Python上でのscipy.arrayの要素アクセスの遅さも注目すべき要素である. このことを考えるとPython内では基本的にlistでデータを管理し, scipyに投げるときやsliceをがんがん使いたい時のみarrayにするのもありかもしれない. 用途によるが.

他にもimport pyublasは, おそらくBoost::Pythonの型変換テーブルへの登録を行なっているとか色々書きたいことはあるが長くなったのでこの辺で.