python - 行列 - ドット積 内積




大きなメモリマップ配列の効率的なドット積 (2)

numpy.memmapの代わりにPyTablesを使うことをお勧めします。 また、圧縮についての彼らのプレゼンテーションを読んで、それは私には奇妙に聞こえるが、シーケンス"compress-> transfer-> uncompress"が速く、圧縮されずに転送されるようです。

また、MKLでnp.dotを使用します。 そして、私はnumexpr( pytablesにも同様のものがあるように見える )が行列の乗算にどのように使用できるか分かりませんが、ユークリッドノルムを計算するのに(numpyと比較して)最も速い方法です。

このサンプルコードをベンチマークしてみてください:

import numpy as np
import tables
import time
n_row=1000
n_col=1000
n_batch=100
def test_hdf5_disk():
    rows = n_row
    cols = n_col
    batches = n_batch
    #settings for all hdf5 files
    atom = tables.Float32Atom()
    filters = tables.Filters(complevel=9, complib='blosc') # tune parameters
    Nchunk = 4*1024  # ?
    chunkshape = (Nchunk, Nchunk)
    chunk_multiple = 1
    block_size = chunk_multiple * Nchunk

    fileName_A = 'carray_A.h5'
    shape_A = (n_row*n_batch, n_col)  # predefined size
    h5f_A = tables.open_file(fileName_A, 'w')
    A = h5f_A.create_carray(h5f_A.root, 'CArray', atom, shape_A, chunkshape=chunkshape, filters=filters)
    for i in range(batches):
        data = np.random.rand(n_row, n_col)
        A[i*n_row:(i+1)*n_row]= data[:]
    rows = n_col
    cols = n_row
    batches = n_batch
    fileName_B = 'carray_B.h5'
    shape_B = (rows, cols*batches)  # predefined size
    h5f_B = tables.open_file(fileName_B, 'w')
    B = h5f_B.create_carray(h5f_B.root, 'CArray', atom, shape_B, chunkshape=chunkshape, filters=filters)
    sz= rows/batches
    for i in range(batches):
        data = np.random.rand(sz, cols*batches)
        B[i*sz:(i+1)*sz]= data[:]
    fileName_C = 'CArray_C.h5'
    shape = (A.shape[0], B.shape[1])
    h5f_C = tables.open_file(fileName_C, 'w')
    C = h5f_C.create_carray(h5f_C.root, 'CArray', atom, shape, chunkshape=chunkshape, filters=filters)
    sz= block_size
    t0= time.time()
    for i in range(0, A.shape[0], sz):
        for j in range(0, B.shape[1], sz):
            for k in range(0, A.shape[1], sz):
                C[i:i+sz,j:j+sz] += np.dot(A[i:i+sz,k:k+sz],B[k:k+sz,j:j+sz])
    print (time.time()-t0)
    h5f_A.close()
    h5f_B.close()
    h5f_C.close()

チャンクサイズと圧縮率を現在のマシンに合わせる方法がわからないという問題は、パフォーマンスがパラメータに依存する可能性があると思います。

また、サンプルコードのすべての行列がディスクに保存されていることに注意してください。その中にはRAMに格納されているものがあれば、高速になると思います。

ところで、私はx32マシンとnumpy.memmapを使っていますが、私は行列サイズにいくつかの制限があります(私は確かではありませんが、ビューサイズは〜2Gbしかないと思われます)。PyTablesには制限がありません。

私は現在、PyTables CArrayのディスク上に常駐している、かなり大きくて密度の高い浮動小数点配列を扱っています。 私はAが巨大な(〜1E4 x 3E5 float32)メモリマップされた配列で、 BCがより小さなnumpyの配列である、 C = A.dot(B)ようなこれらの配列を使って効率的なドットプロダクトを実行できるようにする必要があります。コアメモリに常駐しています。

現在私がやっているのは、 np.memmapを使ってメモリマップされたnumpy配列にデータをコピーし、次にメモリマップされた配列にnp.dot直接呼び出すことnp.dot 。 これはうまくいきますが、標準のnp.dot (またはそれが呼び出す基本的なBLAS関数)は、おそらく結果を計算するために必要なI / O操作の数に関して非常に効率的ではないと思われます。

私はこのレビュー記事で興味深い例を見つけました。 次のように、3xネストループを使用して計算された素朴なドット積:

def naive_dot(A, B, C):
    for ii in xrange(n):
        for jj in xrange(n):
            C[ii,jj] = 0
            for kk in xrange(n):
                C[ii,jj] += A[ii,kk]*B[kk,jj]
    return C

計算にはO(n ^ 3)の I / O操作が必要です。

しかし、適切なサイズのブロックで配列を処理することによって:

def block_dot(A, B, C, M):
    b = sqrt(M / 3)
    for ii in xrange(0, n, b):
        for jj in xrange(0, n, b):
            C[ii:ii+b,jj:jj+b] = 0
            for kk in xrange(0, n, b):
                C[ii:ii+b,jj:jj+b] += naive_dot(A[ii:ii+b,kk:kk+b], 
                                                B[kk:kk+b,jj:jj+b],
                                                C[ii:ii+b,jj:jj+b])
    return C

ここで、 Mはコアメモリに収まる要素の最大数、I / O操作の数はO(n ^ 3 / sqrt(M))に減少します。

np.dotおよび/またはnp.memmapはどのくらいスマートですか? np.dotを呼び出すと、I / O効率の高いブロックnp.dotドット製品が実行されますか? np.memmapはこのタイプの操作の効率を向上させるためのnp.memmapキャッシングを行いますか?

そうでない場合は、I / O効率のよいドットプロダクトを実行する既存のライブラリ関数がいくつかありますか、それとも自分で実装するべきですか?

更新

私は、コアメモリに明示的に読み込まれた入力配列のブロックで動作するnp.dotnp.dotで実装するベンチマーキングを行ってきました。 このデータは少なくとも部分的に私の元の質問に対応しているので、回答として投稿しています。


私はnumpyがmemmap配列のドットプロダクトを最適化するとは思っていませんが、行列乗算のコードを見れば、関数MatrixProduct2 (現在実装されている)は結果行列の値をcメモリの順序:

op = PyArray_DATA(ret); os = PyArray_DESCR(ret)->elsize;
axis = PyArray_NDIM(ap1)-1;
it1 = (PyArrayIterObject *)
    PyArray_IterAllButAxis((PyObject *)ap1, &axis);
it2 = (PyArrayIterObject *)
    PyArray_IterAllButAxis((PyObject *)ap2, &matchDim);
NPY_BEGIN_THREADS_DESCR(PyArray_DESCR(ap2));
while (it1->index < it1->size) {
    while (it2->index < it2->size) {
        dot(it1->dataptr, is1, it2->dataptr, is2, op, l, ret);
        op += os;
        PyArray_ITER_NEXT(it2);
    }
    PyArray_ITER_NEXT(it1);
    PyArray_ITER_RESET(it2);
}

上のコードでは、 opは戻り行列、 dotは1d内積関数、 it1it2は入力行列に対する反復子です。

つまり、あなたのコードはすでに正しいことをしているようです。 この場合、最適なパフォーマンスは実際にはO(n ^ 3 / sprt(M))よりもはるかに優れています.IOをディスクから一回だけ読み込むように制限することができます。 Memmap配列は当然シーンの背後でキャッシュを行う必要があり、内部ループはit2上で動作するため、AがCオーダーでmemmapキャッシュが十分に大きければ、コードはすでに動作している可能性があります。 Aの行を明示的にキャッシュするには、次のようにします。

def my_dot(A, B, C):

    for ii in xrange(n):
        A_ii = np.array(A[ii, :])
        C[ii, :] = A_ii.dot(B)

    return C




linear-algebra