データの局所性とキャッシュメモリの考え方
- あるデータが参照されたら、近くのデータも参照されやすいという考え方
- データアクセスがあったら、付近の連続する一定サイズの領域がまとめてキャッシュされる
- キャッシュメモリのサイズは小さいので、キャッシュが溢れたら使っていないデータはキャッシュから取り除かれる
- メインメモリよりキャッシュメモリの方が高速なので、なるべくキャッシュから読んでもらえるようにプログラムを書くのが大事
C言語で大きい配列を扱うプログラムの高速化でよく言われる話ですが、Python使いの方も意識しましょうというお話です。
論外な例
Pythonのリストへのアクセスは超絶遅いです。一つの理由はデータが不連続だから。
これを大規模データ処理のために使っちゃいけません。
例えば、1000×1000の行列Xと、1000次元の（縦）ベクトルyとの積を計算したいとします。
リストベースだとこんな感じ？（Python 2.7の場合）
X = [[i+j for i in xrange(1000)] for j in xrange(1000)]
y = range(1000, 2000)
z = [sum(a*b for a, b in itertools.izip(r, y)) for r in X]
print z[0:10]
# [832333500, 833833000, 835332500, 836832000, 838331500, 839831000, 841330500, 842830000, 844329500, 845829000]
そのスピードはというと…
%timeit [sum(a*b for a, b in itertools.izip(r, y)) for r in X]
# 10 loops, best of 3: 96.4 ms per loop
これを基準として、numpyで書いたものを見てみましょう。
import numpy
X = numpy.arange(1000).reshape((-1,1)) + numpy.arange(1000)
y = numpy.arange(1000, 2000)
z = X.dot(y)
print z[0:10]
# [832333500 833833000 835332500 836832000 838331500 839831000 841330500
# 842830000 844329500 845829000]
これの計算時間を測ります。
%timeit X.dot(y)
# 1000 loops, best of 3: 700 µs per loop
差は歴然ですね。行列サイズが大きくなると差はどんどん広がります。
あまりにデータが大きいと、1回計算するのも一苦労かもしれません。
元のデータがリストであったとしても、計算前にnumpy.arrayで変換しておけば大丈夫です。
X = [[i+j for i in xrange(1000)] for j in xrange(1000)]
y = range(1000, 2000)
# --- ここまで与えられているとする。 ---
X = numpy.array(X)
y = numpy.array(y)
z = X.dot(y)
print z[0:10]
# [832333500 833833000 835332500 836832000 838331500 839831000 841330500
# 842830000 844329500 845829000]
%timeit X.dot(y)
# 1000 loops, best of 3: 699 µs per loop
ここまでは、numpyを使う大きな動機の一つだと思いますので、何をいまさらという感じですが。
問題はこの先です。
行列積と局所性
今度は、行列Xを転置したものにベクトルyを掛けたいと思います。
素直に考えると、こうやるでしょうね。
z = X.T.dot(y)
%timeit X.T.dot(y)
# 100 loops, best of 3: 2.05 ms per loop
遅い、遅すぎる。
Xを転置しないでyを掛けた場合は699マイクロ秒でしたので、3倍近く遅くなっているのです。
この理由を解明するには、numpy.arrayのメモリ上の配置を考える必要があります。
numpy.arrayでM×N行列を作成すると、各成分はメモリ上で
(0, 0), (0, 1), ..., (0, N-1), (1, 0), ..., (M-1, N-2), (M-1, N-1)
のように連続で配置されます。C言語の多次元配列と同じですね。
X.T は転置行列ではありますが、実態としては X を参照しているだけであり、新しいメモリ領域にコピーされているわけではありません。実験してみましょう。
trX = X.T
trX[1, 2] = 10000
print X[2, 1]
# 10000
このように、転置行列を書き換えると、元の行列の対応する成分まで書き換わってしまいます。
行列積ですから、左側の行列の各行に対して、右側のベクトルとの内積計算が行われます。データの局所性を考えると、その際に行単位でメモリ領域が固まっていると都合が良いわけです。ところが、X.T の各行は X の各列ですから、データのアドレスは飛び飛びになってしまうのですね。
これを速くしたいと思えば、X.T を新しいメモリ領域にコピーするのが良いでしょう。
（コピーするにも時間が掛かるので、同じ計算を何度も繰り返す場合という条件はつきますが。以後同様です）
trX = X.T.copy()
z = trX.dot(y)
%timeit trX.dot(y)
# 1000 loops, best of 3: 704 µs per loop
行列同士の積
統計学でいうところの「共分散行列」を求めるときなどに出てくる、$XX^T$ みたいな形の行列積があります。
ここでも局所性を考えてみましょう。
trX = X.T.copy()
%timeit X.dot(trX)
# 1 loop, best of 3: 2.05 s per loop
%timeit X.dot(X.T)
# 1 loop, best of 3: 685 ms per loop
ご覧のように、転置行列をコピーしたものではなく、単純にX.Tとの積を取った方が速いですね。
右側の行列は列単位での処理になるので、列単位でメモリ領域が固まっていたほうが良いわけです。
X.Tの各列はXの各行になっていて、メモリ領域が固まっています。
一方でtrXとの積を取る場合、列ごとに見ると各成分が飛び飛びになるので、遅くなってしまいます。
行列同士の和
%timeit X+X
# 100 loops, best of 3: 3.81 ms per loop
%timeit X+X.T
# 100 loops, best of 3: 6.37 ms per loop
%timeit X.T+X
# 100 loops, best of 3: 6.39 ms per loop
%timeit X.T+X.T
# 100 loops, best of 3: 3.8 ms per loop
片方だけが転置になっていると遅くなるようですね。
転置同士の場合は速いですが、内部処理で計算順序をよろしく調整してくれているのでしょうか？
その他いろいろ
%timeit X.sum(axis=0)
# 1000 loops, best of 3: 421 µs per loop
%timeit X.sum(axis=1)
# 1000 loops, best of 3: 611 µs per loop
行ごとに総和を求める axis=1 の方が速そうに思えましたが、 列ごとに総和を求める axis=0 の方が若干速い模様です。この辺は実装の都合のような気がしています。
また、ある行列の各行、または各列に対して、独立に1次元の畳み込みを行う例も確認します。
w = numpy.array([2, 1])
%timeit scipy.ndimage.convolve1d(X, w, axis=0)
#100 loops, best of 3: 9.34 ms per loop
%timeit scipy.ndimage.convolve1d(X, w, axis=1)
#100 loops, best of 3: 5.64 ms per loop
こちらは axis=1 の方が速いようです。
まとめ
特に行列積を考えるときには局所性を意識しましょう。同じ計算を何度も実行する場合には、計算対象の行列のメモリ配置に気を配りましょう。
numpyやscipyの高度なライブラリ関数だと、データの内部処理も絡んでくるので、必ずしも予想通りの速度差が出るわけではないようですが。