keisukeのブログ

***乱雑です！自分用のメモです！*** 統計や機械学習の勉強と、読み物を書く練習と、備忘録用のブログ

NumPyから最高のパフォーマンスを得る方法

numpy Python

2018/12/27 追記

この記事はかなり昔の情報が混じっている可能性があります。 元の記事にもいくつか修正が入っているようですが、この翻訳記事では特に修正を実施していません。ご注意ください。 正確な情報は元記事を参照してください。

記事

この記事の和訳です。
色々間違っている可能性があるのでご注意ください。

NumPyはPythonの科学計算ソフトフェア群の基礎となるものです。 NumPyはndarrayというベクトル計算に最適化された特別なデータ構造を提供します。 このオブジェクトは、科学数値計算の中の多くのアルゴリズムの核となっています。

特に計算がひとつの命令で多くのデータを操作する (SIMD) パラダイムに沿っている時、Numpy array (配列)を使うことでネイティブなPythonよりもかなりのパフォーマンスの高速化が達成できます。 しかし、最適化されていないNumPyのコードを何気なく書いてしまうこともありえます。

この記事では、あなたが最適化されたNumPyのコードを書くことを助けるいくつかのトリックを紹介します。 時間とメモリを節約するために、不必要な配列のコピーを避ける方法から見ていきましょう。 そのために、NumPyの内部に分け入る必要があるでしょう。

不必要な配列のコピーを避ける方法を学ぶ

NumPyの配列を使った計算はメモリのかたまりを内部でコピーする処理を含むかもしれません。 このコピーは、避けるべきならば必ずしも必要なものではありません。 あなたのコードを最適化するいくつかのヒントを示しましょう。

import numpy as np

配列のメモリアドレスを調べる

無言でおこなわれる配列のコピーを見つけるための最初のステップは、配列がメモリのどこにいるのか見つけることです。 そのために、次の関数を用意しました：

def id(x):
    # この関数は配列のメモリブロックアドレスを返します
    return x.__array_interface__['data'][0]

時々、配列のコピーを作る必要があります。 例えば、もとの配列をメモリに保存したまま、そのコピーを操作したいときなどです。

a = np.zeros(10); aid = id(a); aid

# 71211328

b = a.copy(); id(b) == aid

# False

同じアドレス(idが返す値)にあるふたつの配列は、内部のデータバッファを共有しています。 しかし、データバッファを共有しているふたつの配列が同じアドレスにあるかというと、それはそれらが同じオフセット(最初の要素が同じかどうか)を持っている時だけです。 次に示すように、異なるオフセットをもつ共有配列はメモリの場所がちょっとだけ変わります：

id(a), id(a[1:])

# (71211328, 71211336)

この記事では、この方法を配列が同じオフセットかどうかを確かめるのに使います。 ふたつの配列が共有されたデータを持っているかどうかを調べるもっと信頼できる方法は：

def get_data_base(arr):
    """与えられたNumPyの配列から、本当のデータを
    「持っている」ベース配列を探す"""
    base = arr
    while isinstance(base.base, np.ndarray):
        base = base.base
    return base

def arrays_share_data(x, y):
    return get_data_base(x) is get_data_base(y)

print(arrays_share_data(a,a.copy()),
      arrays_share_data(a,a[1:]))

# False True

これを指摘してこの代替法を提案してくれたMichael Droettboomに感謝します。

in-placeかつ暗黙のコピーがおこなわれる操作

配列の計算はin-placeな操作(最初の例：配列が変更される)を含んだり、暗黙のコピーを生成する操作(2番目の例：新しい配列が生成される)を含んだりする可能性があります：

a *= 2; id(a) == aid

# True

c = a * 2; id(c) == aid

# False

本当に必要な操作を気をつけて選ばなければなりません。 次に示すように、暗黙にコピーを生成する操作はひどく遅いです：

%%timeit a = np.zeros(10000000)
a *= 2

# 10 loops, best of 3: 19.2 ms per loop

%%timeit a = np.zeros(10000000)
b = a * 2

# 10 loops, best of 3: 42.6 ms per loop

配列の形を変更することはコピーを作るかもしれませんし、作らないかもしれません。 その理由は次のように説明できます。 例えば、2次元行列の形を変更することは、転置(もしくはより一般的に不連続)を取らない限り、コピーを生成しません：

a = np.zeros((10, 10)); aid = id(a); aid

# 53423728

配列の順番を保ったままおこなわれる形の変更はコピーを生成しません。

b = a.reshape((1, -1)); id(b) == aid

# True

転置は配列の順番を変更するので、形の変更はコピーを生成します。

c = a.T.reshape((1, -1)); id(c) == aid

# False

よって、後者のコードは前者のコードよりもかなり遅いことになります。

配列のflattenとravelメソッドは、1次元ベクトルへと形を変更します(flattened array)。 前者のメソッドは必ずコピーを返しますが、後者は必要なときだけコピーを返します(ですから特に大きな配列の場合後者がとても速いです)。

d = a.flatten(); id(d) == aid

# False

e = a.ravel(); id(e) == aid

# True

%timeit a.flatten()

# 1000000 loops, best of 3: 881 ns per loop

%timeit a.ravel()

# 1000000 loops, best of 3: 294 ns per loop

ブロードキャスティングルール

ブロードキャスティングルールは異なる形の配列に対して計算を行える方法です。 言い換えると、必ずしも配列の形を変えたりtileしたりして形を合わせなくても良いのです。 次の例は2つの配列の外積の計算です： 最初の方法はtileを用いていますが、2番目の方法はブロードキャスティングを用いています。 2番目の方法のほうがかなり速いです。

n = 1000

a = np.arange(n)
ac = a[:, np.newaxis]
ar = a[np.newaxis, :]

%timeit np.tile(ac, (1, n)) * np.tile(ar, (n, 1))

# 100 loops, best of 3: 10 ms per loop

%timeit ar * ac

# 100 loops, best of 3: 2.36 ms per loop

NumPyで配列の効果的な抽出をする方法

NunmPyは配列のスライス方法をいくつか提供しています。 配列のviewはもとのデータバッファを参照していますが、オフセットや形、データの間引き方が違います。 NumPyはまた配列のある軸にそった任意の抽出をおこなう特別な関数を提供しています。 最後に、もっとも一般的な抽出方法であるfancy indexingを見ますが、この記事の中でもっとも遅い方法となっています。 より速い代替法が使われるべきです。

たくさんの列を持った配列を作りましょう。 その配列の最初の次元にそったスライスを抽出します。

n, d = 100000, 100

a = np.random.random_sample((n, d)); aid = id(a)

配列のviewとfancy indexing

異なる2つの方法(viewとfancy indexing)を使って10列ごとに抽出をおこないましょう。

b1 = a[::10]
b2 = a[np.arange(0, n, 10)]

np.array_equal(b1, b2)

# True

viewはもとのデータバッファを参照していますが、fancy indexingはコピーを生成します。

id(b1) == aid, id(b2) == aid

# (True, False)

両方のパフォーマンスを比較しましょう。

%timeit a[::10]

# 1000000 loops, best of 3: 804 ns per loop

%timeit a[np.arange(0, n, 10)]

# 100 loops, best of 3: 14.1 ms per loop

fancy indexingは大きな配列のコピーを生成するので、桁外れに遅いです。

fancy indexingの代替法：インデックスのリスト

ある次元に沿った間引きでない抽出方法が必要となった時、viewは選択肢になりません。 しかし、この場合でもfancy indexingの代替法があります。 NumPyの関数takeはインデックスのリストを与えられると、ある次元に沿った抽出をおこないます。

i = np.arange(0, n, 10)

b1 = a[i]
b2 = np.take(a, i, axis=0)

np.array_equal(b1, b2)

# True

2番目の方法のほうが速いです：

%timeit a[i]

# 100 loops, best of 3: 13 ms per loop

%timeit np.take(a, i, axis=0)

# 100 loops, best of 3: 4.87 ms per loop

fancy indexingの代替法：ブーリアンマスク

ある次元に沿った抽出がブーリアンマスクによって決定されている時、compressがfancy indexingの代替となります。

i = np.random.random_sample(n) < .5

fancy indexingまたはcompressを使って抽出ができます。

b1 = a[i]
b2 = np.compress(i, a, axis=0)

np.array_equal(b1, b2)

# True

%timeit a[i]

# 10 loops, best of 3: 59.8 ms per loop

%timeit np.compress(i, a, axis=0)

# 10 loops, best of 3: 24.1 ms per loop

二番目の方法のほうがfancy indexingよりもかなり速いです。

fancy indexingは配列から任意の自由な抽出をおこなうもっとも一般的な方法です。 しかし、大抵もっと用途の限定された速いメソッドが存在していて、可能ならばそちらを使うほうが望ましいです。

配列のviewは間引きによる抽出がおこなわれるときはいつでも使われるべきですが、viewは元のデータバッファを参照しているという事実に注意する必要があります。

どうやって動くの？

このセクションでは、NumPyを使っている時に内部で何が起きているのか、そしてこの知識がどうやってこの記事のトリックを理解するのに役立つのかを紹介します。

どうしてNumPyの配列は効率的なのか？

NumPyの配列は基本的にメタデータ(次元数、形、データ型など)と実際のデータで成り立っています。 それらのデータはシステムメモリの特定のアドレス(ランダムアクセスメモリ、RAM)に同一かつ隙間なく並べられています。 このメモリのかたまりはデータバッファと呼ばれます。 listのような純粋なPythonのデータ構造と主に異なっているところは、アイテムがシステムメモリに散らばっていることです。 この側面がNumPyの配列をこんなに効率的にしているクリティカルな特徴です。

なぜこれがそんなに大切なのでしょうか？ 主な理由がいくつかあります：

  1. 配列の計算はCのような低レベル言語でとても効率的に書かれている (そして実際NumPyのほとんどの部分はCで書かれている) 例えば、メモリブロックとデータ型のアドレスを知ることは、全てのアイテムをループする単純な方法です。 同じことをPythonのlistでやろうとすると大きなオーバーヘッドが発生します。
  2. メモリアクセスパターンの特別な局所性はCPUのキャッシュのおかげで、かなりのパフォーマンスの改善につながります。 実際、キャッシュはRAMからCPUのレジスタにバイトのチャンクで読み込まれます。 隣り合ったアイテムはとても効率的に読み込むことが出来ます (シーケンシャル局所性、または参照の局所性)。
  3. データの要素はメモリに連続して保管されるので、NumPyは現代CPUのベクトル化演算、例えばインテルのSSEやAVX、AMDのXOPなどの恩恵を受けることができます。 例えば、複数の連続した浮動小数点数はCPUの演算として実装されているベクトル演算のために128，256，または512ビットのレジスタにロードされます。

加えて、例えばIntel Math Kernel Library (MKL)を通して、NumPyがBLASやLAPACKのような高度に最適化された線形代数ライブラリとリンクされている事実を伝えましょう。 いくつかの特別な行列計算は現代のマルチコアプロセッサの力を得て並列化できているでしょう。

結論として、メモリブロックにデータを連続して配置することで、メモリアクセスパターン、CPUキャッシュ、ベクトル演算などの意味で、現代CPUのアーキテクチャを適切に用いています。

in-placeと暗黙のコピーを生成する操作の違いは何？

配列のすべての値を二倍するa *= 2のような表現はin-place操作に対応しています。 一方、a = a * 2はa * 2の値を持つ新しい配列が作られることを意味していて、その変数はその新しい配列を指すようになります。 古い配列はもう参照されていないのでガベージコレクターによって回収されます。 最初の場合だとメモリアロケーションはおきませんが、2番目の例だとそうではありません。

もっと一般的には、a[i:j]はある配列の一部へのviewです：データを格納しているメモリバッファを指します。 in-placeな操作でそれらを変えることは、もとの配列も変更することになります。 よって、a[:] = a * 2はin-placeな操作となり、a = a * 2はそうではありません。

この微妙なNumPyの動きを知ることはデバッグに役立ったり(viewに対して操作をおこなったせいで配列が暗黙のうちに意図せず変更される)、不必要なコピーを避ける事でメモリや計算速度の最適化に役立ちます。

コピーせずに形を変えられない配列があるのはなぜ？

上で説明したように、転置された二次元行列はコピーなしでは1次元ベクトルにできません。 二次元行列はふたつのインデックス(行と列)を含みますが、内部的には1次元の連続したメモリブロックに配置されていて、ひとつのインデックスでアクセスできます。 行列の要素を1次元ベクトルで配置する方法はひとつではありません： 1番目の列を配置して、次に二番目の列を配置して、...という方法と、1番目の行を配置して、次に2番めの行を配置して、...という方法とがあります。 最初の方法はrow-major orderと呼ばれ、後者はcolumn-major orderと呼ばれます。 どちらを選ぶのかは慣習でしかありません： NumPyはCと同じrow-major orderですが、FORTRANはcolumn-majorです。 f:id:kaisk:20150219223007p:plain

より一般的には、NumPyは高次元行列のインデックスとその1次元ベクトルのメモリ格納法のあいだで、間引き法をしています。 array[i1, i2]の指す要素と、内部データの対応するバイトアドレスは次のようになります：

offset = array.strides[0] * i1 + array.strides[1] * i2

配列の形を変える時、NumPyはstrides属性の値を変更することで形の変更が可能なら、コピーをおこないません。 たとえば、行列を転置するとき、stridesの順番は反転しますが、内部のデータは変わりません。 しかし、転置行列を1次元ベクトルにしようとすると、strides属性の変更だけでは不可能(試してみてください)なので、コピーが必要となります(ハーバードの Chris Beaumont のおかげでこの段落が簡単になりました)。

Recipe 4.6 (NumPyの間引きトリックを使う)にもっと発展した議論を載せています。 また、recipe 4.7 (間引きトリックで効率的な移動平均アルゴリズムを実装する)は特別な配列計算を高速化する間引きトリックについてを含みます。

配列の内部構造は、よく似たNumPyの式の間での予期しないパフォーマンスの違いについての説明をすることが出来ます。 簡単な練習として、次のベンチマークのを説明できますか？

a = np.random.rand(5000, 5000)
%timeit a[0,:].sum()
%timeit a[:,0].sum()

# 100000 loops, best of 3: 9.57 µs per loop
# 10000 loops, best of 3: 68.3 µs per loop