Tipo eficiente de memoria de matriz numpy masiva en Python

Necesito ordenar un conjunto de datos genómico MUY grande usando numpy. Tengo una matriz de 2.6 billones de flotadores, dimensiones = (868940742, 3) que ocupa unos 20 GB de memoria en mi máquina una vez cargada y simplemente sentada allí. Tengo una MacBook Pro de 13 ‘a principios de 2015 con 16 GB de RAM, 500 GB de estado sólido de alta definición y un procesador Intel i7 de 3.1 GHz. Solo cargando la matriz se desborda a la memoria virtual pero no al punto en que mi máquina sufre o tengo que detener todo lo demás que estoy haciendo.

Construyo esta matriz MUY grande paso a paso a partir de 22 subarrays más pequeños (N, 2) .

La función FUN_1 genera 2 nuevas matrices (N, 1) utilizando cada una de las 22 subarrays que llamo sub_arr .

La primera salida de FUN_1 se genera mediante la interpolación de los valores de sub_arr[:,0] en la matriz b = array([X, F(X)]) y la segunda salida se genera al colocar sub_arr[:, 0] en bandejas utilizando la matriz r = array([X, BIN(X)]) . Llamo a estas salidas b_arr y rate_arr , respectivamente. La función devuelve una tupla de 3 (N, 1) arreglos:

 import numpy as np def FUN_1(sub_arr): """interpolate b values and rates based on position in sub_arr""" b = np.load(bfile) r = np.load(rfile) b_arr = np.interp(sub_arr[:,0], b[:,0], b[:,1]) rate_arr = np.searchsorted(r[:,0], sub_arr[:,0]) # HUGE efficiency gain over np.digitize... return r[rate_r, 1], b_arr, sub_arr[:,1] 

Llamo a la función 22 veces en un bucle for y relleno una matriz de ceros pre-asignados full_arr = numpy.zeros([868940742, 3]) con los valores:

 full_arr[:,0], full_arr[:,1], full_arr[:,2] = FUN_1 

En cuanto a ahorrar memoria en este paso, creo que esto es lo mejor que puedo hacer, pero estoy abierto a sugerencias. De cualquier manera, no tengo problemas hasta este punto y solo toma alrededor de 2 minutos.

Aquí está la rutina de clasificación (hay dos clases consecutivas)

 for idx in range(2): sort_idx = numpy.argsort(full_arr[:,idx]) full_arr = full_arr[sort_idx] # ... #  

Ahora este tipo ha estado funcionando, aunque lentamente (toma alrededor de 10 minutos). Sin embargo, recientemente comencé a usar una tabla de resolución más fina y fina de [X, F(X)] para el paso de interpolación anterior en FUN_1 que devuelve b_arr y ahora el SORT realmente disminuye, aunque todo lo demás sigue siendo el mismo.

Curiosamente, ni siquiera estoy clasificando en los valores interpolados en el paso donde la clasificación ahora está retrasada. Aquí hay algunos fragmentos de los diferentes archivos de interpolación: el más pequeño es aproximadamente un 30% más pequeño en cada caso y mucho más uniforme en términos de valores en la segunda columna; el más lento tiene una resolución más alta y muchos más valores únicos, por lo que los resultados de la interpolación son probablemente más únicos, pero no estoy seguro de si esto debería tener algún tipo de efecto …

archivo más grande, más lento:

 17399307 99.4 17493652 98.8 17570460 98.2 17575180 97.6 17577127 97 17578255 96.4 17580576 95.8 17583028 95.2 17583699 94.6 17584172 94 

Archivo regular más pequeño, más uniforme:

 1 24 1001 24 2001 24 3001 24 4001 24 5001 24 6001 24 7001 24 

No estoy seguro de lo que podría estar causando este problema y me interesaría cualquier sugerencia o simplemente información general sobre la clasificación en este tipo de caso de límite de memoria.

En el momento en que cada llamada a np.argsort está generando una (868940742, 1) de índices int64, que ocupará ~ 7 GB solo. Además, cuando utiliza estos índices para ordenar las columnas de full_arr , está generando otra matriz de flotantes (868940742, 1) , ya que la indexación de fantasía siempre devuelve una copia en lugar de una vista .

Una mejora bastante obvia sería ordenar full_arr en su lugar utilizando su método .sort() . Desafortunadamente, .sort() no le permite especificar directamente una fila o columna para ordenar. Sin embargo, puede especificar un campo para clasificar por una matriz estructurada. Por lo tanto, puede forzar una ordenación in situ sobre una de las tres columnas al obtener una view en su matriz como una matriz estructurada con tres campos flotantes, y luego ordenar por uno de estos campos:

 full_arr.view('f8, f8, f8').sort(order=['f0'], axis=0) 

En este caso, estoy ordenando full_arr en su lugar por el campo 0, que corresponde a la primera columna. Tenga en cuenta que asumí que hay tres columnas float64 ( 'f8' ); debe cambiar esto en consecuencia si su tipo de dtype es diferente. Esto también requiere que su matriz sea contigua y en formato de fila principal, es decir, full_arr.flags.C_CONTIGUOUS == True .

El crédito por este método debe ir a Joe Kington por su respuesta aquí .


Aunque requiere menos memoria, desafortunadamente, la clasificación de una matriz estructurada por campo es mucho más lenta en comparación con el uso de np.argsort para generar una matriz de índice, como mencionó en los comentarios a continuación (consulte esta pregunta anterior ). Si usa np.argsort para obtener un conjunto de índices para clasificar, puede ver una ganancia de rendimiento modesta al usar np.take lugar de la indexación directa para obtener la matriz ordenada:

  %%timeit -n 1 -r 100 x = np.random.randn(10000, 2); idx = x[:, 0].argsort() x[idx] # 1 loops, best of 100: 148 µs per loop %%timeit -n 1 -r 100 x = np.random.randn(10000, 2); idx = x[:, 0].argsort() np.take(x, idx, axis=0) # 1 loops, best of 100: 42.9 µs per loop 

Sin embargo, no esperaría ver ninguna diferencia en términos de uso de memoria, ya que ambos métodos generarán una copia.


Con respecto a su pregunta acerca de por qué la clasificación de la segunda matriz es más rápida, sí, debe esperar que cualquier algoritmo de clasificación razonable sea más rápido cuando hay menos valores únicos en la matriz porque, en promedio, hay menos trabajo que hacer. Supongamos que tengo una secuencia aleatoria de dígitos entre 1 y 10:

 5 1 4 8 10 2 6 9 7 3 

Hay 10! = 3628800 formas posibles de organizar estos dígitos, pero solo una en la que están en orden ascendente. Ahora supongamos que hay solo 5 dígitos únicos:

 4 4 3 2 3 1 2 5 1 5 

Ahora hay 2⁵ = 32 formas de organizar estos dígitos en orden ascendente, ya que podría intercambiar cualquier par de dígitos idénticos en el vector ordenado sin interrumpir el orden.

De forma predeterminada, np.ndarray.sort() utiliza Quicksort . La variante qsort de este algoritmo funciona seleccionando recursivamente un elemento ‘pivote’ en la matriz, luego reordenando la matriz de tal manera que todos los elementos menos que el valor pivote se coloquen delante de él, y todos los elementos mayores que el valor pivote se coloquen después de. Los valores que son iguales al pivote ya están ordenados. Tener menos valores únicos significa que, en promedio, más valores serán iguales al valor de pivote en cualquier barrido dado, y por lo tanto se necesitan menos barridos para ordenar completamente la matriz.

Por ejemplo:

 %%timeit -n 1 -r 100 x = np.random.random_integers(0, 10, 100000) x.sort() # 1 loops, best of 100: 2.3 ms per loop %%timeit -n 1 -r 100 x = np.random.random_integers(0, 1000, 100000) x.sort() # 1 loops, best of 100: 4.62 ms per loop 

En este ejemplo, los tipos de datos de las dos matrices son los mismos. Si su matriz más pequeña tiene un tamaño de elemento más pequeño en comparación con la matriz más grande, el costo de copiarlo debido a la elegante indización también será menor.

EDITAR: En caso de que alguien nuevo en progtwigción y numpy llegue a esta publicación, quiero señalar la importancia de considerar el np.dtype que está usando. En mi caso, en realidad pude evitar el uso del punto flotante de media precisión, es decir, np.float16 , que redujo un objeto de 20GB en la memoria a 5GB y hizo que la clasificación fuera mucho más manejable. El valor predeterminado utilizado por numpy es np.float64 , que es una gran precisión que puede que no necesite. Consulte el documento aquí, que describe la capacidad de los diferentes tipos de datos. Gracias a @ali_m por señalar esto en los comentarios.

Hice un mal trabajo al explicar esta pregunta, pero descubrí algunas soluciones útiles que creo que serían útiles para compartir para cualquiera que necesite clasificar una matriz numpy verdaderamente masiva.

Estoy creando una matriz numpy muy grande a partir de 22 ” numpy ” de datos del genoma humano que contienen los elementos [position, value] . En última instancia, la matriz final debe ordenarse numéricamente “en su lugar” según los valores de una columna en particular y sin mezclar los valores dentro de las filas.

Las dimensiones de la sub-matriz siguen la forma:

 arr1.shape = (N1, 2) ... arr22.shape = (N22, 2) 

sum([N1..N2]) = 868940742 es decir, hay cerca de las posiciones 1BN para clasificar.

Primero, proceso las 22 subarreglas con la función process_sub_arrs , que devuelve una tupla de 3D de la misma longitud que la entrada. Apilo las matrices 1D en una nueva matriz (N, 3) y las inserto en una matriz np.zeros inicializada para el conjunto de datos completo:

  full_arr = np.zeros([868940742, 3]) i, j = 0, 0 for arr in list(arr1..arr22): # indices (i, j) incremented at each loop based on sub-array size j += len(arr) full_arr[i:j, :] = np.column_stack( process_sub_arrs(arr) ) i = j return full_arr 

EDITAR: Como me di cuenta de que mi conjunto de datos podía representarse con flotadores de media precisión, ahora inicializo full_arr siguiente manera: full_arr = np.zeros([868940742, 3], dtype=np.float16) , que es solo 1/4 del tamaño Y mucho más fácil de ordenar.

El resultado es una matriz masiva de 20 GB:

 full_arr.nbytes = 20854577808 

Como @ali_m señaló en su post detallado, mi rutina anterior fue ineficaz:

 sort_idx = np.argsort(full_arr[:,idx]) full_arr = full_arr[sort_idx] 

la matriz sort_idx , que es un 33% del tamaño de full_arr , se cuelga y desperdicia memoria después de clasificar full_arr . Supuestamente, este tipo genera una copia de full_arr debido a la indexación “elegante”, lo que podría full_arr uso de la memoria al 233% de lo que ya se usa para mantener la gran cantidad de elementos. Este es el paso lento, que dura unos diez minutos y depende en gran medida de la memoria virtual.

Sin embargo, no estoy seguro de que el tipo “elegante” haga una copia persistente. Al observar el uso de la memoria en mi máquina, parece que full_arr = full_arr[sort_idx] elimina la referencia al original sin clasificar, porque después de aproximadamente 1 segundo, todo lo que queda es la memoria utilizada por la matriz ordenada y el índice, incluso si hay una copia transitoria.

Un uso más compacto de argsort() para ahorrar memoria es este:

  full_arr = full_arr[full_arr[:,idx].argsort()] 

Esto sigue causando un pico en el momento de la asignación, donde se crean una matriz de índice transitorio y una copia transitoria, pero la memoria se libera casi instantáneamente de nuevo.

@ali_m señaló un buen truco (acreditado a Joe Kington) para generar una matriz estructurada de facto con una view en full_arr . El beneficio es que estos pueden ordenarse “en el lugar”, manteniendo un orden de filas estable:

 full_arr.view('f8, f8, f8').sort(order=['f0'], axis=0) 

Las vistas funcionan muy bien para realizar operaciones de matrices matemáticas, pero para clasificarlas es demasiado ineficaz incluso para una única sub-matriz de mi conjunto de datos. En general, las matrices estructuradas no parecen escalar muy bien a pesar de que tienen propiedades realmente útiles. Si alguien tiene alguna idea de por qué esto es lo que me interesaría saber.

Una buena opción para minimizar el consumo de memoria y mejorar el rendimiento con arreglos muy grandes es construir una tubería de funciones pequeñas y simples. Las funciones borran las variables locales una vez que se han completado, así que si las estructuras de datos intermedios se están acumulando y agotando la memoria, esta puede ser una buena solución.

Este es un bosquejo de la tubería que he usado para acelerar el ordenamiento masivo de matrices:

 def process_sub_arrs(arr): """process a sub-array and return a 3-tuple of 1D values arrays""" return values1, values2, values3 def build_arr(): """build the initial array by joining processed sub-arrays""" full_arr = np.zeros([868940742, 3]) i, j = 0, 0 for arr in list(arr1..arr22): # indices (i, j) incremented at each loop based on sub-array size j += len(arr) full_arr[i:j, :] = np.column_stack( process_sub_arrs(arr) ) i = j return full_arr def sort_arr(): """return full_arr and sort_idx""" full_arr = build_arr() sort_idx = np.argsort(full_arr[:, index]) return full_arr[sort_idx] def get_sorted_arr(): """call through nested functions to return the sorted array""" sorted_arr = sort_arr()  return statistics 

stack de llamadas: get_sorted_arr -> sort_arr -> build_arr -> process_sub_arrs

Una vez que se completa cada función interna, get_sorted_arr() finalmente solo mantiene la matriz ordenada y luego devuelve una pequeña matriz de estadísticas.

EDITAR: También vale la pena señalar aquí que incluso si es capaz de usar un tipo de dtype más compacto para representar su enorme matriz, deseará usar una mayor precisión para los cálculos de resumen. Por ejemplo, desde full_arr.dtype = np.float16 , el comando np.mean(full_arr[:,idx]) intenta calcular la media en coma flotante de media precisión, pero esto se desborda rápidamente al sumr una matriz masiva. El uso de np.mean(full_arr[:,idx], dtype=np.float64) evitará el desbordamiento.

Publiqué esta pregunta inicialmente porque me sorprendió el hecho de que un conjunto de datos de tamaño idéntico de repente comenzó a ahogar la memoria de mi sistema, aunque hubo una gran diferencia en la proporción de valores únicos en el nuevo conjunto “lento”. @ali_m señaló que, de hecho, es más fácil ordenar datos más uniformes con menos valores únicos:

La variante qsort de Quicksort funciona mediante la selección recursiva de un elemento ‘pivote’ en la matriz, y luego reordena la matriz de tal manera que todos los elementos menores que el valor pivote se colocan antes de ella, y todos los elementos mayores que el valor pivote se colocan después eso. Los valores que son iguales al pivote ya están ordenados, por lo que, de manera intuitiva, cuanto menos valores únicos haya en la matriz, menor será el número de swaps que deben realizarse.

En esa nota, el cambio final que terminé haciendo para intentar resolver este problema fue redondear el conjunto de datos más reciente por adelantado, ya que había un nivel innecesariamente alto de precisión decimal restante de un paso de interpolación. En última instancia, esto tuvo un efecto aún mayor que los otros pasos de ahorro de memoria, lo que demuestra que el algoritmo de clasificación en sí fue el factor limitante en este caso.

Espero otros comentarios o sugerencias que cualquiera pueda tener sobre este tema, y ​​es casi seguro que me he equivocado al hablar sobre algunos problemas técnicos, por lo que me gustaría escucharlos 🙂