[Python] 要素数の少ないnumpy arrayで最大値・最小値を求める際の注意点

[PythonでFXシストレのバックテスト](http://qiita.com/toyolab/items/e8292d2f051a88517cb2)
の記事でちょっと触れていますが、[GitHub](https://github.com/toyolab/MT5IndicatorsPy)に公開したテクニカル指標のうち、パラボリックSAR(iSAR)の関数だけがやけに時間がかかっていました。アルゴリズムが複雑だからしょうがないかと思っていたのですが、実は最大値、最小値の求め方に問題があったことが判明しました。

本記事では、サンプルコードを使って問題点をまとめておきます。

# 例題
適当な時系列データから３サンプルずつ取り出して、その最大値を順次求める問題を考えます。式で書くとこんな感じです。

$$y(n)=\max\\{x(n), x(n-1), x(n-2)\\}$$

時系列データは、次のように乱数列として作成しておきます。

````py
import numpy as np
x = np.random.randint(1000, size=100000)
````

# ４種類のコード
３サンプルの最大値を時系列として求めるコードを４種類試してみます。
## func1
式のままコードにします。Pythonの組み込み関数`max`の引数に`x[i],x[i-1],x[i-2]`を渡します。

````py
def func1(x):
    y = np.empty(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = max(x[i], x[i-1], x[i-2])
    return y
````

## func2
ちょっとPythonらしく、スライスを使って`max`の引数に渡します。

````py
def func2(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = max(x[i-2:i+1])
    return y
````

## func3
numpyにも`max`関数があるから、`max`の代わりに`np.max`を使ってみます。

````py
def func3(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = np.max(x[i-2:i+1])
    return y
````

## func4
`x[i],x[i-1],x[i-2]`の３つの要素をリストにして`np.max`の引数に渡してみます。

````py
def func4(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = np.max([x[i], x[i-1], x[i-2]])
    return y
````

# 実行時間の比較

上記の４つの関数の実行時間を比較してみます。

````py
%timeit y1 = func1(x)
%timeit y2 = func2(x)
%timeit y3 = func3(x)
%timeit y4 = func4(x)
````
````
10 loops, best of 3: 91.6 ms per loop
1 loop, best of 3: 304 ms per loop
1 loop, best of 3: 581 ms per loop
1 loop, best of 3: 1.29 s per loop
````

func1が最も速くて、func2, func3, func4と遅くなります。func4はfunc1の14倍時間がかかっています。

# numbaによる高速化
同じコードで、numbaによる高速化の比較を行ってみます。

````py
from numba import jit

@jit
def func1(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = max(x[i], x[i-1], x[i-2])
    return y

@jit
def func2(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = max(x[i-2:i+1])
    return y

@jit
def func3(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = np.max(x[i-2:i+1])
    return y

@jit
def func4(x):
    y = np.zeros(len(x), dtype=int)
    for i in range(2,len(x)):
        y[i] = np.max([x[i], x[i-1], x[i-2]])
    return y
````

````py
%timeit y1 = func1(x)
%timeit y2 = func2(x)
%timeit y3 = func3(x)
%timeit y4 = func4(x)
````
````
1000 loops, best of 3: 365 µs per loop
1 loop, best of 3: 377 ms per loop
100 loops, best of 3: 4.33 ms per loop
1 loop, best of 3: 1.36 s per loop
````

時間の単位に注意して比べてみると、func1がµsなので、圧倒的に速くなっています。次に速いのはfunc3です。この二つはnumbaによる高速化の効果がはっきりとわかります。それに比べて、func2とfunc4はnumbaの効果がないどころか、むしろ遅くなっています。

その結果、func4とfunc1では、3700倍と差が広がっています。結局、少ない要素数で最大値を求める場合、組み込み関数の`max`に要素を個別に渡した方が最も高速に処理できることがわかりました。

実は、numpyのドキュメントの[amax]((http://docs.scipy.org/doc/numpy/reference/generated/numpy.amax.html))のところには、「maximum(a[0], a[1]) is faster than amax(a, axis=0).」と、それらしいコメントが書いてありました。

# iSARが遅かった理由

話を戻すと、iSARが遅かった理由は、func4のような書き方をしていたためでした。それをfunc1のような書き方にしたところ、劇的に速くなったわけです。numbaの効果もあり、これまで最も遅かったテクニカル指標は、一転、高速なテクニカル指標の一つになりました。

ほんとPythonってどこにボトルネックがあるかわからないですね。