大内存映射数组的高效点积

Question

我正在使用一些相当大的,密集的numpy浮点数组,这些数组目前驻留在PyTables CArray的磁盘上.我需要能够执行使用这些阵列效率点的产物,例如C = A.dot(B),其中A是一个巨大的(〜1E4 X 3E5 FLOAT32)存储器映射阵列,以及B和C较小numpy的数组,其驻留在核心存储器.

我现在正在做的是使用数据将数据复制到内存映射的numpy数组中np.memmap,然后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.memmap有任何可以提高此类操作效率的花哨缓存？

如果没有,是否有一些预先存在的库函数可以执行I/O高效的点积,或者我应该自己尝试实现它？

更新

我已经完成了一些基准测试,其中一个手动实现np.dot在输入数组的块上运行,这些块被明确地读入核心内存.这些数据至少部分解决了我原来的问题,所以我将其作为答案发布.

Answer 1

我已经实现了一个函数,用于应用np.dot从内存映射数组中显式读入核心内存的块:

import numpy as np

def _block_slices(dim_size, block_size):
    """Generator that yields slice objects for indexing into 
    sequential blocks of an array along a particular axis
    """
    count = 0
    while True:
        yield slice(count, count + block_size, 1)
        count += block_size
        if count > dim_size:
            raise StopIteration

def blockwise_dot(A, B, max_elements=int(2**27), out=None):
    """
    Computes the dot product of two matrices in a block-wise fashion. 
    Only blocks of `A` with a maximum size of `max_elements` will be 
    processed simultaneously.
    """

    m,  n = A.shape
    n1, o = B.shape

    if n1 != n:
        raise ValueError('matrices are not aligned')

    if A.flags.f_contiguous:
        # prioritize processing as many columns of A as possible
        max_cols = max(1, max_elements / m)
        max_rows =  max_elements / max_cols

    else:
        # prioritize processing as many rows of A as possible
        max_rows = max(1, max_elements / n)
        max_cols =  max_elements / max_rows

    if out is None:
        out = np.empty((m, o), dtype=np.result_type(A, B))
    elif out.shape != (m, o):
        raise ValueError('output array has incorrect dimensions')

    for mm in _block_slices(m, max_rows):
        out[mm, :] = 0
        for nn in _block_slices(n, max_cols):
            A_block = A[mm, nn].copy()  # copy to force a read
            out[mm, :] += np.dot(A_block, B[nn, :])
            del A_block

    return out

然后我做了一些基准测试,将我的blockwise_dot函数与np.dot直接应用于内存映射数组的普通函数进行比较(参见下面的基准测试脚本).我正在使用numpy 1.9.0.dev-205598b链接到OpenBLAS v0.2.9.rc1(从源代码编译).这台机器是运行Ubuntu 13.10的四核笔记本电脑,配备8GB RAM和SSD,我已经禁用了交换文件.

结果

正如@Bi Rico所预测的那样,计算点积所花费的时间相对于尺寸而言是漂亮的O(n)A.在缓存块上操作A比仅调用np.dot整个内存映射数组上的普通函数提供了巨大的性能提升:

它对正在处理的块的大小非常不敏感 - 以1GB,2GB或4GB的块处理阵列所需的时间差别很小.我得出结论,无论缓存np.memmap数组本身实现什么,它对计算点产品似乎都是非常不理想的.

进一步的问题

由于我的代码可能必须在具有不同物理内存量和可能不同的操作系统的机器上运行,因此必须手动实现此缓存策略仍然有点痛苦.出于这个原因,我仍然对是否有办法控制内存映射数组的缓存行为以提高性能感兴趣np.dot.

当我运行基准测试时,我注意到一些奇怪的内存处理行为 - 当我调用np.dot整个时,A我从未看到我的Python进程的驻留集大小超过大约3.8GB,即使我有大约7.5GB的RAM空闲.这让我怀疑对np.memmap阵列允许占用的物理内存量有一些限制- 我之前曾假设它会使用操作系统允许它抓取的任何RAM.在我的情况下,能够增加此限制可能是非常有益的.

有没有人对np.memmap数组的缓存行为有任何进一步的了解,这有助于解释这一点？

基准测试脚本

def generate_random_mmarray(shape, fp, max_elements):
    A = np.memmap(fp, dtype=np.float32, mode='w+', shape=shape)
    max_rows = max(1, max_elements / shape[1])
    max_cols =  max_elements / max_rows
    for rr in _block_slices(shape[0], max_rows):
        for cc in _block_slices(shape[1], max_cols):
            A[rr, cc] = np.random.randn(*A[rr, cc].shape)
    return A

def run_bench(n_gigabytes=np.array([16]), max_block_gigabytes=6, reps=3,
              fpath='temp_array'):
    """
    time C = A * B, where A is a big (n, n) memory-mapped array, and B and C are
    (n, o) arrays resident in core memory
    """

    standard_times = []
    blockwise_times = []
    differences = []
    nbytes = n_gigabytes * 2 ** 30
    o = 64

    # float32 elements
    max_elements = int((max_block_gigabytes * 2 ** 30) / 4)

    for nb in nbytes:

        # float32 elements
        n = int(np.sqrt(nb / 4))

        with open(fpath, 'w+') as f:
            A = generate_random_mmarray((n, n), f, (max_elements / 2))
            B = np.random.randn(n, o).astype(np.float32)

            print "\n" + "-"*60
            print "A: %s\t(%i bytes)" %(A.shape, A.nbytes)
            print "B: %s\t\t(%i bytes)" %(B.shape, B.nbytes)

            best = np.inf
            for _ in xrange(reps):
                tic = time.time()
                res1 = np.dot(A, B)
                t = time.time() - tic
                best = min(best, t)
            print "Normal dot:\t%imin %.2fsec" %divmod(best, 60)
            standard_times.append(best)

            best = np.inf
            for _ in xrange(reps):
                tic = time.time()
                res2 = blockwise_dot(A, B, max_elements=max_elements)
                t = time.time() - tic
                best = min(best, t)
            print "Block-wise dot:\t%imin %.2fsec" %divmod(best, 60)
            blockwise_times.append(best)

            diff = np.linalg.norm(res1 - res2)
            print "L2 norm of difference:\t%g" %diff
            differences.append(diff)

        del A, B
        del res1, res2
        os.remove(fpath)

    return (np.array(standard_times), np.array(blockwise_times), 
            np.array(differences))

if __name__ == '__main__':
    n = np.logspace(2,5,4,base=2)
    standard_times, blockwise_times, differences = run_bench(
                                                    n_gigabytes=n,
                                                    max_block_gigabytes=4)

    np.savez('bench_results', standard_times=standard_times, 
             blockwise_times=blockwise_times, differences=differences)

Answer 2

我不认为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点乘积函数,it1并且it2是输入矩阵上的迭代器.

话虽这么说,看起来你的代码可能已经做了正确的事情.在这种情况下,最佳性能实际上比O(n ^ 3/sprt(M))好得多,您可以将IO限制为仅从磁盘读取A的每个项目,或者O(n).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

Answer 3

我建议您使用PyTables而不是numpy.memmap。也请阅读他们关于压缩的介绍，这对我来说听起来很奇怪，但是似乎序列“ 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没有限制。