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

ゆとりデータサイエンティストの諸々所感

データ分析会社で研究開発をしている、ゆとり世代データサイエンティストが学んだ内容や最新トピックについて諸々語る予定

Python高速化 Numba入門 その3

今回もNumbaのドキュメントを読んで行きます。

Numba — numba 0.15.1 documentation

と思ったんですが、読み進めて行くと以外に紹介する内容が少ないことに気づきました。

シンプルなのは良いことなので、今回はUFuncを紹介して、
次回にGPUについて紹介して入門編は一旦終わりたいと思います。

その後に、Numbaを勉強する動機となったものを中心に応用編を話す予定です。

UFunc

UFuncs — numba 0.15.1 documentation

UFunc(ユニバーサル関数)とは、要するにnumpy.ndarrayに対して要素ごとに演算してくれる関数です。

例えばnumpy.logに配列を食わせると、こんな感じです。

>>> import numpy
>>> hoge = numpy.arange(1, 5)
>>> hoge
array([1, 2, 3, 4])
>>> numpy.log(hoge)
array([ 0.        ,  0.69314718,  1.09861229,  1.38629436])

numpy自体にもUFuncの書き方はあるみたいですが、Numbaだとこれもデコレータ一発で行けるようなので見ていきます。

基本は以下の構文。@vectorizeデコレータにリストで型を指定する。

from numba import vectorize
import math
import numpy

@vectorize(['float64(float64, float64)'])
def my_ufunc(x, y):
    return x+y+math.sqrt(x*math.cos(y))

a = numpy.arange(1.0, 10.0)
b = numpy.arange(1.0, 10.0)
# Calls compiled version of my_ufunc for each element of a and b
print(my_ufunc(a, b))

結果は

[  2.73505259          nan          nan          nan  11.1909286
  14.40021285  16.29724091          nan          nan]

マイナスのルート取ろうてして・・・nanが発生。これが例でいいのか・・・

複数の型を受け取りたい時はデコレータの引数のリストにひたすら足していく。

from numba import vectorize
import math
import numpy

@vectorize(['int32(int32, int32)',
            'float64(float64, float64)'])
def my_ufunc(x, y):
    return x+y+math.sqrt(x*math.cos(y))

a = numpy.arange(1.0, 10.0, dtype='f8')
b = numpy.arange(1.0, 10.0, dtype='f8')
print(my_ufunc(a, b))

a = numpy.arange(1, 10, dtype='i4')
b = numpy.arange(1, 10, dtype='i4')
print(my_ufunc(a, b))
[  2.73505259          nan          nan          nan  11.1909286
  14.40021285  16.29724091          nan          nan]
C:\Python27\Scripts\ipython-script.py:16: RuntimeWarning: invalid value encountered in ufunc
[          2 -2147483648 -2147483648 -2147483648          11          14
          16 -2147483648 -2147483648]

@jitみたいに@vectorizeにも型推定ないかなと思い実行。

from numba import vectorize
import math
import numpy

@vectorize
def my_ufunc(x, y):
    return x+y+math.sqrt(x*math.cos(y))

a = numpy.arange(1.0, 10.0)
b = numpy.arange(1.0, 10.0)
# Calls compiled version of my_ufunc for each element of a and b
print(my_ufunc(a, b))
TypeError: wrap() takes exactly 1 argument (2 given)

@vectorizeはしっかり型を指定する必要があるらしい。

既存の関数を後から、UFuncBuilderでUFuncにする方法もある。

from numba.npyufunc.ufuncbuilder import UFuncBuilder
from numba import f8, i4

import math
import numpy

def my_ufunc(x, y):
    return x+y+math.sqrt(x*math.cos(y))

builder = UFuncBuilder(my_ufunc)
builder.add(restype=i4, argtypes=[i4, i4])
builder.add(restype=f8, argtypes=[f8, f8])

compiled_ufunc = builder.build_ufunc()

a = numpy.arange(1.0, 10.0, dtype='f8')
b = numpy.arange(1.0, 10.0, dtype='f8')
print(compiled_ufunc(a, b))

Generalized Ufuncs

こちらは私もよく分かっていませんが、わかっている範囲でご紹介します。

通常のUFuncは要素ごとの処理を行っていましたが、Generalized Ufuncsでは内部の配列ベースで処理を行ってくれるそうです。

C言語でも単に要素のforを回すよりも、配列としてBLASに投げた方が速いので動機としては理解できます。

import math
import numpy
from numba import guvectorize

@guvectorize(['void(int32[:,:], int32[:,:], int32[:,:])',
              'void(float64[:,:], float64[:,:], float64[:,:])'],
              '(x, y),(x, y)->(x, y)')
def my_gufunc(a, b, c):
    for i in range(c.shape[0]):
        for j in range(c.shape[1]):
            c[i, j] = a[i, j] + b[i, j]

a = numpy.arange(1.0, 10.0, dtype='f8').reshape(3,3)
b = numpy.arange(1.0, 10.0, dtype='f8').reshape(3,3)
# Calls compiled version of my_gufunc for each row of a and b
print(my_gufunc(a, b))

型を指定するところは@vectorizeと同じですが、デコレータの第2引数と、本体の関数に戻り値となる引数が増える様です。

デコレータの引数は内部の配列をどう扱うかを指定しています。詳しい話は以下を参考にして下さい。

Generalized Universal Function API — NumPy v1.9 Manual

ちょっと特殊ですが慣れれば行けそうです。

実際にやってみた

前回の@jitと@guvectorizeはどちらが速いか検証します。

# coding: utf-8
import numpy
from numba import double
from numba.decorators import jit
from numba import guvectorize
import time
import traceback

def pairwise_python(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)
    return D

@jit
def pairwise_numba(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)


@jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True)
def pairwise_numba2(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)

    return D


@guvectorize(['void(f8[:, :], f8[:, :])'], '(x, y)->(x, x)')
def pairwise_numba3(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)



if __name__ == '__main__':

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_python(X, D)
    print "python           :", time.time() - t

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba(X, D)
    print "numba_jit        :", time.time() - t

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba2(X, D)
    print "numba_jit_type   :", time.time() - t

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba3(X)
    print "numba_guvectorize:", time.time() - t

結果は・・・

python           : 4.63400006294
numba_jit        : 0.131000041962
numba_jit_type   : 0.010999917984
numba_guvectorize: 0.0119998455048

特に@guvectorizeだからと言って速いわけではないようだ、
私の書き方が正しくない気もするが・・・

次回はGPUの話をする予定
GPU積んでるのがWindows機なので若干手間取るかもしれませんが、頑張ります。

Python高速化 Numba入門 その2

今回は、QuickStartを読んでいきます。

Quick Start — numba 0.15.1 documentation

とりあえず、前回の@jitデコレータだけで動くのは理解した。

from numba import jit

@jit
def sum(x, y):
    return x + y

引数と戻り値の型が陽にわかっている場合には、@jitの引数に『戻り値の型』(引数1の型, 引数2の型, ...)で指定できる。恐らく早くなるのだろう

@jit('f8(f8,f8)')
def sum(x, y):
    return x + y

numpyの配列もサポート。多次元配列はこんな感じf8[:, :, :]

@jit('f8(f8[:])')
def sum1d(array):
    sum = 0.0
    for i in range(array.shape[0]):
        sum += array[i]
    return sum

サポートされている型は以下。主要なものは大体ある。

Type Name Alias Result Type
boolean b1 uint8 (char)
bool_ b1 uint8 (char)
byte u1 unsigned char
uint8 u1 uint8 (char)
uint16 u2 uint16
uint32 u4 uint32
uint64 u8 uint64
char i1 signed char
int8 i1 int8 (char)
int16 i2 int16
int32 i4 int32
int64 i8 int64
float_ f4 float32
float32 f4 float32
double f8 float64
float64 f8 float64
complex64 c8 float complex
complex128 c16 double complex

型は文字列で指定しなくても、importして指定もできる。ハードコーディングにならないのでこっちのが嬉しい。

from numba import jit, f8

@jit(f8(f8[:]))
def sum1d(array):
    ...

型がわからない場合は、通常のPythonで処理するので異常に遅くなってしまう。それを避けるためにPythonでの処理を禁止することもできる。
どうしても型がわからない場合はエラー吐くのかな?

@jit(nopython=True)
def sum1d(array):
    ...

Numbaの型推論の結果を取得するには、inspect_typesメソッドが用意されている。これを参考にデコレータの引数を決めるのも良さそう

sum1d.inspect_types()

実際にやってみる

とりあえず、型がわからない時に、nopython=Trueするとどうなるか

from numba.decorators import jit
import numpy

class DummyClass(object):
    def __init__(self):
        hoge = 0
        huga = 2
        hogo = "hogehoge"

@jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True)
def pairwise_numba3(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)

    return DummyClass()

結果は、やはり推論出来ない型でnopython=Trueはエラーらしい

TypingError: Failed at nopython frontend
Untyped global name 'DummyClass'

型指定したほうが速いのか?

検証のために以下を実行

# coding: utf-8
import numpy
from numba import double
from numba.decorators import jit
import time
import traceback

@jit
def pairwise_numba(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)

    return D

@jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True)
def pairwise_numba2(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)

    return D

if __name__ == '__main__':

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba(X, D)
    print "numba:", time.time() - t

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba2(X, D)
    print "numba2:", time.time() - t

結果は・・・

numba: 0.134999990463
numba2: 0.0120000839233

10倍高速化!

前回、普通のPythonから33倍高速になったので、型指定まですると330倍高速化したことに。

Numba結構イケるやん!

次回は、ユーザーガイドの続きか、First Steps with numbaを読む予定

Python高速化 Numba入門 その1

みなさん、こんにちは
今日からPython高速化 Numbaに入門したいと思います。

入門資料を探しに来た皆様すみませんが、 本記事は私がこれから入門する内容になります。

結果として入門資料に慣れば幸いですが、過度な期待は御無用でお願いします。

基本的には以下を読み進めて行きます。

http://numba.pydata.org/

Numbaとは

JIT(just-in-time)コンパイラを使ってPythonを高速化しよう!』というPythonモジュールです。

LLVMコンパイラを使っており、これはJuliaが高速な理由でもあるので期待大です。

学生時代はCythonを使って高速化をよくしていましたが、以下の理由により今回はNumbaを学びます

候補 今回諦めた理由
Cython cdefとか結構手を入れるのでPythonに戻すのが面倒。pyximportも面倒
C拡張 C言語は極力触りたくない
PyPy numpyサポートが怪しい。いつか触るかも
Julia なるべくPythonでやりたい。これでダメなら入門するかも。

Numbaはどうやらデコレータ一発で一応動くらしい。Cythonよりは使いやすいことを期待したい。

とりあえず通常Pythonと比較

以下のサイトのコードを参考に速度を計測

Numba vs Cython

ちょっと手を加えて、実行したのが以下のコードで比較

# coding: utf-8
import numpy
from numba import double
from numba.decorators import jit
import time

@jit
def pairwise_numba(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)


def pairwise_python(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = numpy.sqrt(d)


if __name__ == '__main__':
    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_python(X, D)
    print "python:", time.time() - t

    t = time.time()
    X = numpy.random.random((1000, 3))
    D = numpy.empty((1000, 1000))
    pairwise_numba(X, D)
    print "numba:", time.time() - t

結果は

python: 4.64800000191
numba: 0.137000083923

34倍の高速化!イケるやん!

次回はQuick Startを読む予定

最適化超入門

SlideShareだけでなく、ブログの記事にもすることに

先日、TokyoWebMning #40にて最適化について熱く語ってきました。

個人的にも結構やりきった感があり、網羅的に最適化手法を紹介出来たと思います。

その後飲んだ研究室の同期には『難しすぎる』と言われましたが、どうなんでしょう・・・

一応はてブが400を超えており、大丈夫だったと信じたい。

というか、最適化の話題ではてブこんなに頂けるのは予想外でした。ありがとうございます!

はてなブックマーク - 最適化超入門

PyData Tokyoの方にもお声掛け頂いたのでまたどこかでお話出来ればと思ってます。
次はセクシー女優みたいにしょーもない話にする予定ですw  
 
 

今回は修士卒の人間が最適化の入門資料を作る事のは、おこがましいと戦々恐々していたのですが、 Twitter上では概ね好評であり、安心しております。

NIIの前原先生からもお褒め頂いたので、最適化チートシート流行らそうと思いますw

f:id:tkm2261:20141206141012p:plain

あと、阪大の梅谷先生からご指摘があり、CPLEXは出来る子とのことです。

CPLEXはコア単位課金でお高いため、企業ではほぼGurobi一択な現状なので誇大に書きすぎたようです。
あとSOS制約の入った問題が高速だったため、うちの会社ではCPLEXからGurobiに移行しました。

研究で使われるかたはどちらもフリーなので、mpsファイルでも吐いて検証してみることをオススメします。

ちなみに、Gurobiは先日バージョン6.0がリリースされ、イベントに参加してきました。

Company - News - Highlights of Gurobi Optimizer 6.0

目玉としては、分散MIPと区分線形関数でした。
MIPも分散する時代になったか・・・といった感じです

区分線形関数もSOS制約でいちいち書かなくて済むので、応用事例が増えることが期待出来そうです。

合わせてBixby氏からNFLの試合日程最適化の話があり、非常に興味深い話でした。 2014年から木曜日に試合を行うことなり、モデルが異常に難しくなったようです。

新日鉄住金ソリューションズの山本さんからもJリーグの日程最適化の話があり、
うちの会社でもスポーツ系やりたい( ´Д`) 

長くなりましたが、益々最適化が流行ればいいな!と思う次第です。

Word2vecで大谷翔平の二刀流論争に終止符を打つ!

皆様、お久しぶりです。

VOYAGE GROUPさん主催の14' Data Scientist MeetUpでLTした時の資料 3時間で作った資料なので超絶適当です。

リファクタしてから上げようかと思いましたが、忙しくて断念

ネタ的にお蔵入りは勿体無かったので、そのまま上げています。

気が向いたらなんかやるかも

Scikit-learnで学ぶ機械学習入門

 今回はこの前勉強会で話してきたこの話

 

『Scikit-learnで学ぶ機械学習入門 』

ついに、このブログでもデータ解析っぽいことを話せて感無量です。

 

詳細な勉強会の模様は天丼丸さんのページをご参照ください。

機械学習勉強会 #2 | /home/by-natures/*

 

 

この勉強会では、機械学習が初めての方も多かったので、

機械学習の全体的な流れを重視して話して来ました。

 

特にp.11の、機械学習の流れは汎用的に使えるかと思います

精度が上がらない時は、このフローに抜け漏れが無いかを確認することが重要です。

f:id:tkm2261:20141002010817p:plain

 

scikit-leanの使い方については、ほとんど解説してませんが、

APIが全て統一されているので、以下だけでほぼ全て使えます。

  • 使用したい手法のインスタンスの生成
  • パラメータのセット
  • fit()関数による学習
  • predict()関数による予測

githubにサンプルコードも載せておいたので、必要でしたら見て下さい。

anaguma2261/scikit-learn-sample · GitHub

 

あとは、各種の手法の利点欠点などをまとめておいたので、

これから機械学習やろうとしている方の一助になればと思います

 

dplyrとR界隈について

Pythonianの私としてはブログの記事として、

Pythonよりも先にRについて書くのは若干気になるものの、

同僚のR使いに面白い共有を頂いたので備忘のため記事にします。

※同僚の受け売りなので、私はあんまりRに詳しくないです

 

近年、R界隈に彗星のごとく現れた、Hadley神ことHadley Wickham氏

had.co.nz

hadley (Hadley Wickham) · GitHub

 

同僚によると、彼の登場により、R界隈に大きい変化の波が来るとのこと

日本語の資料もないので、日本ではほとんど話題になってないが、時間の問題らしい

 

彼の書いた、Advanced R (http://adv-r.had.co.nz/)は

Rのパッケージの作成方法から、

テストといったコードの品質保証まで網羅的に書かれており、

Rをあまり使わない私としても非常に参考になる

 

これにより近年のR(笑)の流れに歯止めがかかり、

Tokyo.RでもRがメイントピックになる日が来る! かな?

 

ウチの会社でも、分析のみの案件では良いものの、システム化が絡むとRは

カワサキバイクのコピペ並に、エンジニアの顔が歪みます

       (それよりも分析官が、テストを書いていないことに起因することが多いですが・・・)

 

ともあれ、分析官にR使いは多いので、

Rがシステム化に耐えられる様になるのは、なにより喜ばしいこと

この流れはこの辺ウォッチしとけば追えるのかな


R-bloggers | R news & tutorials from the web

 

 

前置きが長くなりましたが、dplyrについて語ります。

といっても先人たちがすでに良記事をたくさん書いているので

使い方はそっちを見てください。この記事では思想のみ語ります。

 

日本語資料

英語資料

 

まだまだR-pubs漁れば、たくさん出てきます。

 

dplyrの特徴としては、

  • (Rなのに)速い!
  • filter()、 arrange()、 select()、 mutate()、summarise()の5つの動詞に集約された洗練された設計
  • DB上のテーブルも同じ操作で扱える
  • 近年流行りのパイプ(『%>%』←こいつ)でコードの見渡しが非常に良い。
    SQLライク(筆者感)にもコードが書ける。

 

特にDB上のテーブルを簡単に扱えるのは強力かと

 

Pythonでも Pandas + SQLAlchemy で似たことはできるが、

SQLAlchemyがフルスタックのORマッパなので、

高機能過ぎて学習に時間がかかったり、手順がどうしても多くなる

 

その点、Rはデータ分析に特化しているので、シンプルに機能を絞り込めているので、

使い勝手が良さそう。

 

久しぶりにRで心を動かされたので記事にしてみました。

私もRしっかり勉強するかな・・・