2.4. コードの最適化

Donald Knuth

「早すぎる最適化は諸悪の根源である」

著者: Gaël Varoquaux

この章では Python コードを速く実行するための戦略について扱います。

事前準備

2.4.1. 最適化ワークフロー

  1. 動くようにしましょう: コードを単純で 判読可能に 書きましょう。

  2. 期待通りに動くようにしましょう: 自動化したテストケースを書きましょう、アルゴリズムが正しいことを確かめて、壊したときにテストが壊れたことをつかめるようにしましょう。

  3. 最適かします、単純なユーザケースをプロファイルしてボトルネックをみつけ、そのボトルネックを速くしましょう、もっとよいアルゴリズムや実装を見つけましょう。 現実的な例でのプロファイリングすることには、単純さや実行スピードのトレードオフがあることを念頭に入れましょう。効率のためには、プロファイリングは10秒程度で終了することが望ましいです。

2.4.2. Python コードのプロファイル

測定せずに最適化してはいけません!

  • 測定: プロファイリング、時間測定

  • 驚くことに: 最速のコードはあなたが考えるものと違います

2.4.2.1. Timeit

IPython では timeit (https://docs.python.org/library/timeit.html) を時間計測の基本操作として利用します:

In [1]: import numpy as np



In [2]: a = np.arange(1000)



In [3]: %timeit a ** 2

100000 loops, best of 3: 5.73 us per loop



In [4]: %timeit a ** 2.1

1000 loops, best of 3: 154 us per loop



In [5]: %timeit a * a

100000 loops, best of 3: 5.56 us per loop

これを利用して最適化戦略を選ぶ助けが得られます。

注釈

長い時間かかる呼び出しには %time を %timeit のかわりに利用しましょう; より不正確な結果になりますが、高速です。

2.4.2.2. プロファイラ

大きなプログラムをプロファイルするときに便利です、例えば following file のように:

# For this example to run, you also need the 'ica.py' file



import numpy as np

from scipy import linalg



from ica import fastica





def test():

    data = np.random.random((5000, 100))

    u, s, v = linalg.svd(data)

    pca = np.dot(u[:, :10].T, data)

    results = fastica(pca.T, whiten=False)



if __name__ == '__main__':

    test()

注釈

これは2つの教師無し学習の組み合わせで、主成分解析 (PCA) と独立成分解析 (ICA) を組み合わせたものです。PCA は次元削減のテクニックで、つまり、データの分散を少ない次元を使って説明するアルゴリズムです。ICA は入力を分離するテクニックで、例えば複数のセンサーで記録された複数の信号をよりわけます。センサーの数よりも信号の方が多い場合には、PCA を最初に実行してその後に ICA することが有用です。詳細については : the FastICA example from scikits-learn を参照して下さい。

これを実行するには ica モジュール も同様にダウンロードする必要があります。IPython では以下のようにしてスクリプトを時間計測できます。

In [1]: %run -t demo.py



IPython CPU timings (estimated):

    User  :    14.3929 s.

    System:   0.256016 s.

プロファイルはこうします:

In [2]: %run -p demo.py



      916 function calls in 14.551 CPU seconds



Ordered by: internal time



ncalls  tottime  percall  cumtime  percall filename:lineno (function)

     1   14.457   14.457   14.479   14.479 decomp.py:849 (svd)

     1    0.054    0.054    0.054    0.054 {method 'random_sample' of 'mtrand.RandomState' objects}

     1    0.017    0.017    0.021    0.021 function_base.py:645 (asarray_chkfinite)

    54    0.011    0.000    0.011    0.000 {numpy.core._dotblas.dot}

     2    0.005    0.002    0.005    0.002 {method 'any' of 'numpy.ndarray' objects}

     6    0.001    0.000    0.001    0.000 ica.py:195 (gprime)

     6    0.001    0.000    0.001    0.000 ica.py:192 (g)

    14    0.001    0.000    0.001    0.000 {numpy.linalg.lapack_lite.dsyevd}

    19    0.001    0.000    0.001    0.000 twodim_base.py:204 (diag)

     1    0.001    0.001    0.008    0.008 ica.py:69 (_ica_par)

     1    0.001    0.001   14.551   14.551 {execfile}

   107    0.000    0.000    0.001    0.000 defmatrix.py:239 (__array_finalize__)

     7    0.000    0.000    0.004    0.001 ica.py:58 (_sym_decorrelation)

     7    0.000    0.000    0.002    0.000 linalg.py:841 (eigh)

   172    0.000    0.000    0.000    0.000 {isinstance}

     1    0.000    0.000   14.551   14.551 demo.py:1 (<module>)

    29    0.000    0.000    0.000    0.000 numeric.py:180 (asarray)

    35    0.000    0.000    0.000    0.000 defmatrix.py:193 (__new__)

    35    0.000    0.000    0.001    0.000 defmatrix.py:43 (asmatrix)

    21    0.000    0.000    0.001    0.000 defmatrix.py:287 (__mul__)

    41    0.000    0.000    0.000    0.000 {numpy.core.multiarray.zeros}

    28    0.000    0.000    0.000    0.000 {method 'transpose' of 'numpy.ndarray' objects}

     1    0.000    0.000    0.008    0.008 ica.py:97 (fastica)

     ...

明らかに svd (decamp.py 内の) が一番時間がかかっています、つまりこれがボトルネックです。このステップをより速くするか、あるいはこのステップを避ける(アルゴリズム的な最適化)方法を探す必要があります。コードの他の部分に時間をかけるのは無意味です。

IPython の外でプロファイルするには ``cProfile`` を実行します

同様に IPython の外でプロファイルすることもできます、組み込みの Python プロファイラ や cProfile そして profile を呼び出します。

$  python -m cProfile -o demo.prof demo.py

-o スイッチを使うことでプロファイラの結果を demo.prof ファイルに出力し、外部ツールで見ることができます。プロファイラの結果を可視化ツールで処理したいときに便利に使えます。

2.4.2.3. 行毎のプロファイラ

プロファイラはどの関数が多くの時間かかったかを教えてくれますが、どこで呼び出されたかは教えてくれません。

そのため line_profiler を利用します: ソースファイル内で調べたいいくつかの関数を @profile (import する必要はありません) でデコレートします。

@profile

def test():

    data = np.random.random((5000, 100))

    u, s, v = linalg.svd(data)

    pca = np.dot(u[:, :10], data)

    results = fastica(pca.T, whiten=False)

次に kernprof.py プログラムを利用してスクリプトを実行します、スイッチとして -l, --line-by-line と -v, --view を利用して行毎のプロファイルと結果を保存しつつ閲覧もできるようにします:

$ kernprof.py -l -v demo.py



Wrote profile results to demo.py.lprof

Timer unit: 1e-06 s



File: demo.py

Function: test at line 5

Total time: 14.2793 s



Line #      Hits         Time  Per Hit   % Time  Line Contents

=========== ============ ===== ========= ======= ==== ========

    5                                           @profile

    6                                           def test():

    7         1        19015  19015.0      0.1      data = np.random.random((5000, 100))

    8         1     14242163 14242163.0   99.7      u, s, v = linalg.svd(data)

    9         1        10282  10282.0      0.1      pca = np.dot(u[:10, :], data)

   10         1         7799   7799.0      0.1      results = fastica(pca.T, whiten=False)

SVD が時間の全てを費やしています この行を最適化する必要があります。

2.4.3. コードの高速化

ボトルネックがわかったら、対応するコードを速くする必要があります。

2.4.3.1. アルゴリズム的最適化

最初にやることはアルゴリズム的最適化です: より少ない計算やもっとよい計算法はないでしょうか?

問題をより高度な視点から見ると、アルゴリズムの背後にある数学への十分な理解は問題解決を助けてくれます。またそれにより 計算やメモリ割り当てを for ループの外に出す ような大きな利益をもたらす小さな変更を見つけることは珍しくありません。

2.4.3.1.1. SVD の例

上の例の両方で SVD - Singular Value Decomposition - が時間のほとんどを費していました。このアルゴリズムの計算コストはおおざっぱに入力行列のサイズに対して \(n^3\) です。

しかし、これらの例では SVD の出力全てを使っておらず、最初に返された引数の最初の数行しか利用していません。scipy の svd 実装を利用する場合、不完全なバージョンの SVD を利用するように要請することができます。ここで scipy の線形代数の実装は numpy のものよりも充実していて、より望ましいと、言及しておきます。

In [3]: %timeit np.linalg.svd(data)

1 loops, best of 3: 14.5 s per loop



In [4]: from scipy import linalg



In [5]: %timeit linalg.svd(data)

1 loops, best of 3: 14.2 s per loop



In [6]: %timeit linalg.svd(data, full_matrices=False)

1 loops, best of 3: 295 ms per loop



In [7]: %timeit np.linalg.svd(data, full_matrices=False)

1 loops, best of 3: 293 ms per loop

この見識を 前のコードの最適化 に利用してみましょう:

def test():

    data = np.random.random((5000, 100))

    u, s, v = linalg.svd(data, full_matrices=False)

    pca = np.dot(u[:, :10].T, data)

    results = fastica(pca.T, whiten=False)

In [1]: import demo



In [2]: %timeit demo.

demo.fastica   demo.np        demo.prof.pdf  demo.py        demo.pyc

demo.linalg    demo.prof      demo.prof.png  demo.py.lprof  demo.test



In [2]: %timeit demo.test()

ica.py:65: RuntimeWarning: invalid value encountered in sqrt

  W = (u * np.diag(1.0/np.sqrt(s)) * u.T) * W  # W = (W * W.T) ^{-1/2} * W

1 loops, best of 3: 17.5 s per loop



In [3]: import demo_opt



In [4]: %timeit demo_opt.test()

1 loops, best of 3: 208 ms per loop

例えば固有ベクトルの内最初の10を計算するような、実数の不完全な SVD については、arpack で計算することができ scipy.sparse.linalg.eigsh で利用できます。

線形代数の計算

特定のアルゴリズムについては、ボトルネックは線形代数の計算であることが多いです。今回のケースでは、問題を解くのに適した関数を使うことが鍵となりました。例えば対象行列の固有値の問題は一般的な行列よりも解くのが簡単です。同様に、多くの場合、逆行列を求めるのを避けて、よりコストの低い(そして、数値的に安全な)演算を利用することができます。

線形代数の計算について知りましょう。本当かどうか疑うなら scipy.linalg を探して %timeit を利用してデータに対して別の代替手段を試してみましょう。

2.4.4. 高速な数値計算コードを書く

numpy の先進的な利用方法に対する一通りの議論が Numpy の先進的な機能 や van der Walt 達による The NumPy array: a structure for efficient numerical computation があります。ここではコードを、いくつかのよく遭遇する高速化のための工夫について議論することにします。

  • ループのベクトル化

    numpy の配列に対する for ループを避ける方法を知りましょう。そのためにマスクやインデックス配列が便利に使えます。

  • ブロードキャスト

    broadcasting を利用して、配列を結合させる前にする演算を最小限にしましょう。

  • インプレースな演算

    In [1]: a = np.zeros(1e7)

    

    In [2]: %timeit global a ; a = 0*a

    10 loops, best of 3: 111 ms per loop

    

    In [3]: %timeit global a ; a *= 0

    10 loops, best of 3: 48.4 ms per loop

    注意: 動作させるために global a が timeit の中に必要です、それが a に代入されて a はローカル変数として扱われます。

  • メモリに優しく: コピーではなくビューを使いましょう

    大きな配列をコピーするのは、単純な数値演算するのと同じくらい、コストがかかります:

    In [1]: a = np.zeros(1e7)

    

    In [2]: %timeit a.copy()

    10 loops, best of 3: 124 ms per loop

    

    In [3]: %timeit a + 1

    10 loops, best of 3: 112 ms per loop

  • キャッシュ効果に気をつけましょう

    グループ化されていればメモリアクセスは安価に済みます: 大きな配列に対する連続アクセスはランダムアクセスに比べて高速です。これはつまり 小さなストライドがより高速 ということです (CPU キャッシュの効果 を参照してください):

    In [1]: c = np.zeros((1e4, 1e4), order='C')

    

    In [2]: %timeit c.sum(axis=0)

    1 loops, best of 3: 3.89 s per loop

    

    In [3]: %timeit c.sum(axis=1)

    1 loops, best of 3: 188 ms per loop

    

    In [4]: c.strides

    Out[4]: (80000, 8)

    これが Fortran 順序か C 順序かで演算操作に大きな違いがある理由です:

    In [5]: a = np.random.rand(20, 2**18)

    

    In [6]: b = np.random.rand(20, 2**18)

    

    In [7]: %timeit np.dot(b, a.T)

    1 loops, best of 3: 194 ms per loop

    

    In [8]: c = np.ascontiguousarray(a.T)

    

    In [9]: %timeit np.dot(b, c)

    10 loops, best of 3: 84.2 ms per loop

    データコピーして対処するのはおそらく効果がないことを注意しておきます:

    In [10]: %timeit c = np.ascontiguousarray(a.T)

    10 loops, best of 3: 106 ms per loop

    numexpr を利用するとこれらの効果に対する最適化を自動でしてくれます。

  • コンパイルされたコードの利用

    最後の手段は、高レベルな最適化全てを実施したなら、ホットスポット、つまり、最も時間のかかる数行または関数、をコンパイルしたコードにすることです。コンパイルされたコードにするには Cython を利用するのが望ましいでしょう: 既存の Python コードをコンパイルされたコードにするのは簡単にでき numpy support をうまく利用することで numpy 配列を効率よく扱えます、例えばループの展開など。

警告

以上のように: 選択肢をプロファイル、時間計測しましょう。理論的な検討に基づいて最適化してはいけません。