MemoryError al crear un producto cartesiano en Numpy

Tengo 3 matrices numpy y necesito formar el producto cartesiano entre ellas. Las dimensiones de las matrices no son fijas, por lo que pueden tomar valores diferentes, un ejemplo podría ser A = (10000, 50), B = (40, 50), C = (10000,50).

Luego, realizo algunos procesos (como a + bc) A continuación se muestra la función que estoy usando para el producto.

def cartesian_2d(arrays, out=None): arrays = [np.asarray(x) for x in arrays] dtype = arrays[0].dtype n = np.prod([x.shape[0] for x in arrays]) if out is None: out = np.empty([n, len(arrays), arrays[0].shape[1]], dtype=dtype) m = n // arrays[0].shape[0] out[:, 0] = np.repeat(arrays[0], m, axis=0) if arrays[1:]: cartesian_2d(arrays[1:], out=out[0:m, 1:, :]) for j in range(1, arrays[0].shape[0]): out[j * m:(j + 1) * m, 1:] = out[0:m, 1:] return out a = [[ 0, -0.02], [1, -0.15]] b = [[0, 0.03]] result = cartesian_2d([a,b,a]) // array([[[ 0. , -0.02], [ 0. , 0.03], [ 0. , -0.02]], [[ 0. , -0.02], [ 0. , 0.03], [ 1. , -0.15]], [[ 1. , -0.15], [ 0. , 0.03], [ 0. , -0.02]], [[ 1. , -0.15], [ 0. , 0.03], [ 1. , -0.15]]]) 

La salida es la misma que con itertools.product . Sin embargo, estoy usando mi función personalizada para aprovechar las numerosas operaciones vectorizadas, que funcionan bien en comparación con itertools.product en mi caso.

Después de esto, hago

 result[:, 0, :] + result[:, 1, :] - result[:, 2, :] //array([[ 0. , 0.03], [-1. , 0.16], [ 1. , -0.1 ], [ 0. , 0.03]]) 

Así que este es el resultado final esperado.

    La función funciona como se espera, siempre y cuando mi matriz encaje en la memoria. Pero mi caso de uso requiere que trabaje con datos enormes y obtengo un MemoryError en la línea np.empty() ya que no puede asignar la memoria requerida. Estoy trabajando con datos de alrededor de 20 GB en este momento y esto podría boost en el futuro.

    Estas matrices representan vectores y deberán almacenarse en float , por lo que no puedo usar int . Además, son matrices densas, por lo que el uso sparse no es una opción.

    Estaré usando estos arreglos para procesamiento adicional e idealmente no me gustaría almacenarlos en archivos en esta etapa. Por memmap h5py formato memmap / h5py puede no ser útil, aunque no estoy seguro de esto.

    Si hay otras formas de formar este producto, eso también estaría bien.

    Como estoy seguro de que hay aplicaciones con conjuntos de datos mucho más grandes que esto, espero que alguien haya encontrado este tipo de problemas antes y quisiera saber cómo manejar este problema. Por favor ayuda.

    Si al menos tu resultado cabe en la memoria.

    Lo siguiente produce su resultado esperado sin depender de un intermediario tres veces el tamaño del resultado. Utiliza la radiodifusión.

    Tenga en cuenta que casi cualquier operación NumPy se puede transmitir de esta manera, por lo que en la práctica probablemente no sea necesario un producto cartesiano explícito:

     #shared dimensions: sh = a.shape[1:] aba = (a[:, None, None] + b[None, :, None] - a[None, None, :]).reshape(-1, *sh) aba #array([[ 0. , 0.03], # [-1. , 0.16], # [ 1. , -0.1 ], # [ 0. , 0.03]]) 

    Direccionando las filas de resultados por ‘ID’

    Puede considerar dejar de lado la reshape . Eso le permitiría direccionar las filas en el resultado mediante un índice combinado. Si sus ID de componente son solo 0,1,2, … como en su ejemplo, esto sería lo mismo que la ID combinada. Por ejemplo, aba [1,0,0] correspondería a la fila obtenida como segunda fila de a + primera fila de b – primera fila de a.

    Un poco de explicacion

    Transmisión: cuando, por ejemplo, agregar dos matrices, sus formas no tienen que ser idénticas, solo son compatibles debido a la transmisión. La difusión es, en cierto sentido, una generalización de la adición de escalares a los arreglos:

      [[2], [[7], [[2], 7 + [3], equiv to [7], + [3], [4]] [7]] [4]] 

    Radiodifusión:

      [[4], [[1, 2, 3], [[4, 4, 4], [[1, 2, 3]] + [5], equiv to [1, 2, 3], + [5, 5, 5], [6]] [1, 2, 3]] [6, 6, 6]] 

    Para que esto funcione, cada dimensión de cada operando debe ser 1 o igual a la dimensión correspondiente en cada otro operando, a menos que sea 1. Si un operando tiene menos dimensiones que los otros, su forma se rellena con unas a la izquierda . Tenga en cuenta que las matrices equivalentes que se muestran en la ilustración no se crean explícitamente.

    Si el resultado tampoco encaja

    En ese caso no veo cómo posiblemente puedas evitar usar el almacenamiento, así que h5py o algo así es.

    Eliminando la primera columna de cada operando

    Esto es sólo una cuestión de rebanar:

     a_no_id = a[:, 1:] 

    etc. Tenga en cuenta que, a diferencia de las listas de Python, las matrices NumPy cuando se dividen no devuelven una copia sino una vista. Por lo tanto, la eficiencia (memoria o tiempo de ejecución) no es un problema aquí.

    Una solución alternativa es crear un producto cartesiano de índices (que es más fácil, ya que existen soluciones para productos cartesianos de arreglos 1D):

     idx = cartesian_product( np.arange(len(a)), np.arange(len(b)) + len(a), np.arange(len(a)) ) 

    Y luego use la indexación de fantasía para crear la matriz de salida:

     x = np.concatenate((a, b)) result = x[idx.ravel(), :].reshape(*idx.shape, -1) 

    Escribir resultados eficientemente en el disco.

    Al principio algunas mentes sobre el tamaño de los datos resultantes.

    Tamaño de los datos de resultados.

     size_in_GB = A.shape[0]**2*A.shape[1]*B.shape[0]*(size_of_datatype)/1e9 

    En su pregunta mencionó A.shape = (10000,50), B = (40,50). Usando float64 su resultado será de aproximadamente 1600 GB. Esto se puede hacer sin problemas si tiene suficiente espacio en el disco, pero debe pensar qué va a hacer con los datos a continuación. Tal vez esto sea solo un resultado intermedio y el procesamiento de datos en bloques es posible.

    Si este no es el caso, aquí hay un ejemplo de cómo manejar 1600 GB de datos de manera eficiente (el uso de RAM será de unos 200 MB). El rendimiento del canal debe ser de alrededor de 200 MB / s en datos realistas.

    El código que calcula los resultados es de @PaulPanzer.

     import numpy as np import tables #register blosc import h5py as h5 import h5py_cache as h5c a=np.arange(500*50).reshape(500, 50) b=np.arange(40*50).reshape(40, 50) # isn't well documented, have a look at https://github.com/Blosc/hdf5-blosc compression_opts=(0, 0, 0, 0, 5, 1, 1) compression_opts[4]=9 #compression level 0...9 compression_opts[5]=1 #shuffle compression_opts[6]=1 #compressor (I guess that's lz4) File_Name_HDF5='Test.h5' f = h5.File(File_Name_HDF5, 'w',chunk_cache_mem_size=1024**2*300) dset = f.create_dataset('Data', shape=(a.shape[0]**2*b.shape[0],a.shape[1]),dtype='d',chunks=(a.shape[0]*b.shape[0],1),compression=32001,compression_opts=(0, 0, 0, 0, 9, 1, 1), shuffle=False) #Write the data for i in range(a.shape[0]): sh = a.shape[1:] aba = (a[i] + b[:, None] - a).reshape(-1, *sh) dset[i*a.shape[0]*b.shape[0]:(i+1)*a.shape[0]*b.shape[0]]=aba f.close() 

    Leyendo los datos

     File_Name_HDF5='Test.h5' f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300) dset=f['Data'] chunks_size=500 for i in range(0,dset.shape[0],chunks_size): #Iterate over the first column data=dset[i:i+chunks_size,:] #avoid excessive calls to the hdf5 library #Do something with the data f.close() f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300) dset=f['Data'] for i in range(dset.shape[1]): # Iterate over the second dimension # fancy indexing eg[:,i] will be much slower # use np.expand_dims or in this case np.squeeze after the read operation from the dset # if you wan't to have the same result than [:,i] (1 dim array) data=dset[:,i:i+1] #Do something with the data f.close() 

    En este ejemplo de prueba, obtengo un rendimiento de escritura de aproximadamente 550 M / s, una lectura a través de aproximadamente (500 M / s primer dim, 1000M / s segundo dim) y una relación de compresión de 50. El memmap de Numpy solo proporcionará una velocidad aceptable si lee o escribe datos a lo largo de la dirección de cambio más rápido (en C la última dimensión), con un formato de datos fragmentados utilizado por HDF5 aquí, esto no es un problema en absoluto. La compresión tampoco es posible con el memmap Numpy, lo que lleva a tamaños de archivo más altos y una velocidad más lenta.

    Tenga en cuenta que el filtro de compresión y la forma del trozo deben configurarse según sus necesidades. Esto depende de cómo no desee leer los datos posteriormente y los datos reales.

    Si hace algo completamente incorrecto, el rendimiento puede ser de 10 a 100 veces más lento en comparación con una forma adecuada de hacerlo (por ejemplo, la forma de fragmentos puede optimizarse para el primer o el segundo ejemplo de lectura).