python - Productos de puntos eficientes de grandes matrices mapeadas en memoria

arrays performance

Question

Estoy trabajando con algunas matrices flotantes numpy densas y bastante grandes que actualmente residen en el disco en PyTables CArray s. Necesito poder realizar productos de puntos eficientes usando estas matrices, por ejemplo C = A.dot(B) , donde A es una matriz mapeada de memoria enorme (~ 1E4 x 3E5 float32), y B y C son matrices numpy más pequeñas que son residentes en la memoria del núcleo.

Lo que estoy haciendo en este momento es copiar los datos en matrices numpy mapeadas en memoria usando np.memmap , y luego llamar a np.dot directamente en las matrices mapeadas en memoria. Esto funciona, pero sospecho que el np.dot estándar (o más bien las funciones BLAS subyacentes que llama) probablemente no sea muy eficiente en términos del número de operaciones de E / S requeridas para calcular el resultado.

Encontré un ejemplo interesante en este artículo de revisión . Un producto de puntos ingenuo calculado utilizando 3x bucles anidados, como este:

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

requiere O (n ^ 3) operaciones de E / S para calcular.

Sin embargo, procesando las matrices en bloques de tamaño apropiado:

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

donde M es la cantidad máxima de elementos que cabrán en la memoria del núcleo, el número de operaciones de E / S se reduce a O (n ^ 3 / sqrt (M)) .

¿Qué tan inteligente es np.dot y / o np.memmap ? ¿ np.dot a np.dot realizar un producto de punto de bloqueo eficiente de E / S? ¿ np.memmap hace algún caché elegante que mejore la eficiencia de este tipo de operación?

Si no es así, ¿existe alguna función de biblioteca preexistente que realice productos de E / S eficientes, o debería intentar implementarla yo mismo?

Actualizar

Hice algunas evaluaciones comparativas con una implementación np.dot de np.dot que opera en bloques de la matriz de entrada, que se leen explícitamente en la memoria del núcleo. Esta información al menos parcialmente responde a mi pregunta original, así que la estoy publicando como respuesta.

Answer 1

Implementé una función para aplicar np.dot a bloques que se leen explícitamente en la memoria del núcleo de la matriz mapeada en memoria:

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

Luego hice algunos benchmarking para comparar mi función blockwise_dot con la función np.dot normal aplicada directamente a una matriz mapeada en memoria (ver a continuación el script de benchmarking). Estoy usando numpy 1.9.0.dev-205598b vinculado con OpenBLAS v0.2.9.rc1 (compilado de la fuente). La máquina es una computadora portátil de cuatro núcleos con Ubuntu 13.10, con 8 GB de RAM y una SSD, y he desactivado el archivo de intercambio.

Resultados

Como predijo @Bi Rico, el tiempo necesario para calcular el producto escalar es maravillosamente O (n) con respecto a las dimensiones de A Operar en bloques en caché de A brinda una gran mejora en el rendimiento sobre la simple invocación de la función np.dot normal en toda la matriz mapeada en memoria:

Es sorprendentemente insensible al tamaño de los bloques que se procesan: hay muy poca diferencia entre el tiempo necesario para procesar la matriz en bloques de 1 GB, 2 GB o 4 GB. Concluyo que cualquiera que sea el cache de np.memmap matrices de np.memmap implementadas de forma nativa, parece ser muy poco óptimo para el cálculo de productos de puntos.

Mas preguntas

Todavía es un poco doloroso tener que implementar manualmente esta estrategia de almacenamiento en caché, ya que es probable que mi código tenga que ejecutarse en máquinas con diferentes cantidades de memoria física y sistemas operativos potencialmente diferentes. Por esa razón, todavía estoy interesado en si hay formas de controlar el comportamiento de almacenamiento en caché de las matrices mapeadas en memoria para mejorar el rendimiento de np.dot .

np.dot comportamiento de manejo de memoria extraño mientras ejecutaba los puntos de referencia: cuando llamé a np.dot en todo A , nunca vi que el tamaño del conjunto residente de mi proceso de Python superara los 3.8 GB, a pesar de que tengo aproximadamente 7.5 GB de RAM gratis. Esto me lleva a sospechar que existe un límite impuesto sobre la cantidad de memoria física que puede ocupar una matriz np.memmap : anteriormente había supuesto que usaría la RAM que el sistema operativo le permitiera capturar. En mi caso, podría ser muy beneficioso poder aumentar este límite.

¿Alguien tiene más información sobre el comportamiento de almacenamiento en caché de matrices np.memmap que ayudaría a explicar esto?

Guión de Benchmarking

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

No creo que numpy optimice el producto de puntos para las matrices de memmap, si miras el código de multiplicación de matrices, que obtuve here , verás que la función MatrixProduct2 (como se implementa actualmente) calcula los valores de la matriz de resultados en c orden de memoria:

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);
}

En el código anterior, op es la matriz de retorno, dot es la función de producto de 1d punto e it1 e it2 son iteradores sobre las matrices de entrada.

Dicho esto, parece que su código ya está haciendo lo correcto. En este caso, el rendimiento óptimo es en realidad mucho mejor que O (n ^ 3 / sprt (M)), puede limitar su IO a solo leer cada elemento de A una vez desde el disco, o O (n). Los arreglos de Memmap naturalmente tienen que hacer algo de almacenamiento en caché detrás de la escena y el bucle interno opera en it2 , por lo que si A está en orden C y el caché de memmap es lo suficientemente grande, es posible que su código ya esté funcionando. Puede aplicar el almacenamiento en caché de filas de A de forma explícita haciendo algo como:

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

Te recomiendo que uses PyTables en lugar de numpy.memmap. También leo sus presentaciones sobre la compresión, me parece extraño, pero parece que la secuencia "comprimir-> transferir-> descomprimir" es más rápida que simplemente transferir sin comprimir .

También use np.dot con MKL. Y no sé cómo el numexpr (las tablas también parecen tener algo así ) se puede usar para la multiplicación de matrices, pero por ejemplo para calcular la norma euclidiana es la forma más rápida (en comparación con el numpy).

Intente comparar este código de muestra:

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()

El problema es que no sé cómo ajustar el tamaño del fragmento y la tasa de compresión a la máquina actual, por lo que creo que el rendimiento puede depender de los parámetros.

También tenga en cuenta que todas las matrices en el código de muestra se almacenan en el disco, si algunos de ellos se almacenan en la memoria RAM, creo que será más rápido.

Por cierto estoy usando la máquina x32 y con numpy.memmap tengo algunas limitaciones en el tamaño de la matriz (no estoy seguro, pero parece que el tamaño de la vista puede ser solo ~ 2Gb) y las PyTables no tienen limitaciones.