python - ¿Existe una ventaja de velocidad de análisis o uso de memoria al usar HDF5 para almacenamiento de matriz grande(en lugar de archivos binarios planos)?

numpy pytables

Question

Estoy procesando grandes matrices 3D, que a menudo necesito cortar de varias maneras para hacer una variedad de análisis de datos. Un "cubo" típico puede ser de ~ 100 GB (y probablemente se hará más grande en el futuro)

Parece que el formato de archivo típico recomendado para grandes conjuntos de datos en Python es usar HDF5 (h5py o pytables). Mi pregunta es: ¿hay algún beneficio de velocidad o uso de memoria al usar HDF5 para almacenar y analizar estos cubos en lugar de almacenarlos en archivos binarios planos simples? ¿HDF5 es más apropiado para datos tabulares, en comparación con grandes matrices como con las que estoy trabajando? Veo que HDF5 puede proporcionar una buena compresión, pero estoy más interesado en la velocidad de procesamiento y en lidiar con el desbordamiento de memoria.

Con frecuencia quiero analizar solo un gran subconjunto del cubo. Parece que un inconveniente de pytables y h5py es que cuando tomo una porción de la matriz, siempre obtengo una matriz numpy, usando memoria. Sin embargo, si divido un mapa de memoria numpy de un archivo binario plano, puedo obtener una vista, que mantiene los datos en el disco. Entonces, parece que puedo analizar más fácilmente sectores específicos de mis datos sin saturar mi memoria.

He explorado pytables y h5py, y no he visto el beneficio de ninguno hasta ahora para mi propósito.

Answer 1

Ventajas de HDF5: organización, flexibilidad, interoperabilidad

Algunas de las principales ventajas de HDF5 son su estructura jerárquica (similar a las carpetas / archivos), los metadatos arbitrarios opcionales almacenados con cada elemento y su flexibilidad (por ejemplo, compresión). Esta estructura organizativa y almacenamiento de metadatos puede parecer trivial, pero es muy útil en la práctica.

Otra ventaja de HDF es que los conjuntos de datos pueden ser de tamaño fijo o flexible. Por lo tanto, es fácil agregar datos a un gran conjunto de datos sin tener que crear una copia completamente nueva.

Además, HDF5 es un formato estandarizado con bibliotecas disponibles para casi cualquier idioma, por lo que compartir sus datos en el disco entre, digamos Matlab, Fortran, R, C y Python es muy fácil con HDF. (Para ser justos, no es demasiado difícil con una gran matriz binaria, siempre que conozca el orden C vs. F y conozca la forma, el tipo, etc. de la matriz almacenada).

Ventajas de HDF para una gran variedad: E / S más rápida de un segmento arbitrario

Al igual que el TL / DR: para una matriz 3D de ~ 8GB, leer un segmento "completo" a lo largo de cualquier eje tomó ~ 20 segundos con un conjunto de datos HDF5 fragmentado, y 0.3 segundos (mejor caso) a más de tres horas (peor caso) durante una matriz memmapped de los mismos datos.

Más allá de las cosas enumeradas anteriormente, hay otra gran ventaja de un formato de datos en disco "fragmentado" * como HDF5: la lectura de un segmento arbitrario (énfasis en arbitrario) generalmente será mucho más rápido, ya que los datos en disco son más contiguos en promedio.

* (HDF5 no tiene que ser un formato de datos fragmentados. Es compatible con la fragmentación, pero no lo requiere. De hecho, el valor predeterminado para crear un conjunto de datos en h5py no es fragmentar, si no recuerdo h5py ).

Básicamente, la velocidad de lectura del disco en el mejor de los casos y la velocidad de lectura del disco en el peor de los casos para una porción determinada de su conjunto de datos estará bastante cerca con un conjunto de datos HDF fragmentado (suponiendo que elija un tamaño de fragmento razonable o deje que una biblioteca elija uno para usted). Con una matriz binaria simple, el mejor de los casos es más rápido, pero el peor de los casos es mucho peor.

Una advertencia, si tiene un SSD, es probable que no note una gran diferencia en la velocidad de lectura / escritura. Sin embargo, con un disco duro normal, las lecturas secuenciales son mucho, mucho más rápidas que las lecturas aleatorias. (es decir, un disco duro normal tiene un tiempo de seek prolongado). HDF todavía tiene una ventaja en un SSD, pero se debe más a sus otras características (por ejemplo, metadatos, organización, etc.) que a la velocidad bruta.

En primer lugar, para aclarar la confusión, acceder a un conjunto de datos h5py devuelve un objeto que se comporta de manera bastante similar a una matriz numpy, pero no carga los datos en la memoria hasta que se corta. (Similar a memmap, pero no idéntico). h5py introducción de h5py para obtener más información.

Al dividir el conjunto de datos se cargará un subconjunto de los datos en la memoria, pero presumiblemente desea hacer algo con él, en cuyo momento lo necesitará en la memoria de todos modos.

Si desea realizar cálculos fuera del núcleo, puede obtener datos tabulares con pandas o pytables con bastante facilidad. Es posible con h5py (mejor para grandes matrices ND), pero debe desplegarse a un nivel inferior táctil y manejar la iteración usted mismo.

Sin embargo, el futuro de los cálculos fuera de núcleo como Numpy es Blaze. Echa un vistazo si realmente quieres tomar esa ruta.

El caso "no fragmentado"

En primer lugar, considere una matriz 3D ordenada en C escrita en el disco (lo arr.ravel() llamando a arr.ravel() e imprimiendo el resultado, para que las cosas sean más visibles):

In [1]: import numpy as np

In [2]: arr = np.arange(4*6*6).reshape(4,6,6)

In [3]: arr
Out[3]:
array([[[  0,   1,   2,   3,   4,   5],
        [  6,   7,   8,   9,  10,  11],
        [ 12,  13,  14,  15,  16,  17],
        [ 18,  19,  20,  21,  22,  23],
        [ 24,  25,  26,  27,  28,  29],
        [ 30,  31,  32,  33,  34,  35]],

       [[ 36,  37,  38,  39,  40,  41],
        [ 42,  43,  44,  45,  46,  47],
        [ 48,  49,  50,  51,  52,  53],
        [ 54,  55,  56,  57,  58,  59],
        [ 60,  61,  62,  63,  64,  65],
        [ 66,  67,  68,  69,  70,  71]],

       [[ 72,  73,  74,  75,  76,  77],
        [ 78,  79,  80,  81,  82,  83],
        [ 84,  85,  86,  87,  88,  89],
        [ 90,  91,  92,  93,  94,  95],
        [ 96,  97,  98,  99, 100, 101],
        [102, 103, 104, 105, 106, 107]],

       [[108, 109, 110, 111, 112, 113],
        [114, 115, 116, 117, 118, 119],
        [120, 121, 122, 123, 124, 125],
        [126, 127, 128, 129, 130, 131],
        [132, 133, 134, 135, 136, 137],
        [138, 139, 140, 141, 142, 143]]])

Los valores se almacenarían en el disco secuencialmente como se muestra en la línea 4 a continuación. (Ignoremos los detalles del sistema de archivos y la fragmentación por el momento).

In [4]: arr.ravel(order=''C'')
Out[4]:
array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
        13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
        26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
        39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
        52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
        65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
        78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
        91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
       104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
       117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
       130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143])

En el mejor de los casos, tomemos un corte a lo largo del primer eje. Tenga en cuenta que estos son solo los primeros 36 valores de la matriz. ¡Esta será una lectura muy rápida! (una búsqueda, una lectura)

In [5]: arr[0,:,:]
Out[5]:
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35]])

Del mismo modo, el siguiente segmento a lo largo del primer eje serán los siguientes 36 valores. Para leer un corte completo a lo largo de este eje, solo necesitamos una operación de seek . Si todo lo que vamos a leer son varios cortes a lo largo de este eje, entonces esta es la estructura de archivo perfecta.

Sin embargo, consideremos el peor de los casos: un corte a lo largo del último eje.

In [6]: arr[:,:,0]
Out[6]:
array([[  0,   6,  12,  18,  24,  30],
       [ 36,  42,  48,  54,  60,  66],
       [ 72,  78,  84,  90,  96, 102],
       [108, 114, 120, 126, 132, 138]])

Para leer este segmento, necesitamos 36 búsquedas y 36 lecturas, ya que todos los valores están separados en el disco. ¡Ninguno de ellos es adyacente!

Esto puede parecer bastante menor, pero a medida que llegamos a matrices cada vez más grandes, el número y el tamaño de las operaciones de seek aumentan rápidamente. Para una matriz 3D de gran tamaño (~ 10 Gb) almacenada de esta manera y leída a través de memmap , leer un segmento completo a lo largo del "peor" eje puede llevar fácilmente decenas de minutos, incluso con hardware moderno. Al mismo tiempo, un corte a lo largo del mejor eje puede tomar menos de un segundo. Para simplificar, solo estoy mostrando secciones "completas" a lo largo de un solo eje, pero exactamente lo mismo sucede con secciones arbitrarias de cualquier subconjunto de datos.

Por cierto, hay varios formatos de archivo que aprovechan esto y básicamente almacenan tres copias de enormes matrices 3D en el disco: una en orden C, una en orden F y una en el intermedio entre las dos. (Un ejemplo de esto es el formato D3D de Geoprobe, aunque no estoy seguro de que esté documentado en ningún lado). ¿A quién le importa si el tamaño final del archivo es de 4 TB, el almacenamiento es barato! Lo loco de eso es que debido a que el caso de uso principal es extraer una sola subcarte en cada dirección, las lecturas que desea hacer son muy, muy rápidas. ¡Funciona muy bien!

El simple caso "fragmentado"

Digamos que almacenamos "fragmentos" 2x2x2 de la matriz 3D como bloques contiguos en el disco. En otras palabras, algo como:

nx, ny, nz = arr.shape
slices = []
for i in range(0, nx, 2):
    for j in range(0, ny, 2):
        for k in range(0, nz, 2):
            slices.append((slice(i, i+2), slice(j, j+2), slice(k, k+2)))

chunked = np.hstack([arr[chunk].ravel() for chunk in slices])

Entonces los datos en el disco se verían chunked :

array([  0,   1,   6,   7,  36,  37,  42,  43,   2,   3,   8,   9,  38,
        39,  44,  45,   4,   5,  10,  11,  40,  41,  46,  47,  12,  13,
        18,  19,  48,  49,  54,  55,  14,  15,  20,  21,  50,  51,  56,
        57,  16,  17,  22,  23,  52,  53,  58,  59,  24,  25,  30,  31,
        60,  61,  66,  67,  26,  27,  32,  33,  62,  63,  68,  69,  28,
        29,  34,  35,  64,  65,  70,  71,  72,  73,  78,  79, 108, 109,
       114, 115,  74,  75,  80,  81, 110, 111, 116, 117,  76,  77,  82,
        83, 112, 113, 118, 119,  84,  85,  90,  91, 120, 121, 126, 127,
        86,  87,  92,  93, 122, 123, 128, 129,  88,  89,  94,  95, 124,
       125, 130, 131,  96,  97, 102, 103, 132, 133, 138, 139,  98,  99,
       104, 105, 134, 135, 140, 141, 100, 101, 106, 107, 136, 137, 142, 143])

Y solo para mostrar que son bloques de arr de 2x2x2, observe que estos son los primeros 8 valores de chunked :

In [9]: arr[:2, :2, :2]
Out[9]:
array([[[ 0,  1],
        [ 6,  7]],

       [[36, 37],
        [42, 43]]])

Para leer en cualquier segmento a lo largo de un eje, leemos en 6 o 9 fragmentos contiguos (el doble de datos que necesitamos) y luego solo conservamos la parte que deseamos. Es un máximo de 9 búsquedas en el peor de los casos frente a un máximo de 36 búsquedas para la versión no fragmentada. (Pero el mejor caso sigue siendo 6 búsquedas vs 1 para la matriz memmapped.) Debido a que las lecturas secuenciales son muy rápidas en comparación con las búsquedas, esto reduce significativamente la cantidad de tiempo que lleva leer un subconjunto arbitrario en la memoria. Una vez más, este efecto se hace más grande con matrices más grandes.

HDF5 lleva esto unos pasos más allá. Los fragmentos no tienen que almacenarse de forma contigua, y están indexados por un B-Tree. Además, no tienen que ser del mismo tamaño en el disco, por lo que se puede aplicar compresión a cada fragmento.

Matrices h5py con h5py

Por defecto, h5py no crea archivos HDF fragmentados en el disco (creo que pytables hace, por el contrario). Sin embargo, si especifica chunks=True al crear el conjunto de datos, obtendrá una matriz fragmentada en el disco.

Como un ejemplo rápido y mínimo:

import numpy as np
import h5py

data = np.random.random((100, 100, 100))

with h5py.File(''test.hdf'', ''w'') as outfile:
    dset = outfile.create_dataset(''a_descriptive_name'', data=data, chunks=True)
    dset.attrs[''some key''] = ''Did you want some metadata?''

Tenga en cuenta que chunks=True le dice a h5py que elija automáticamente un tamaño de fragmento para nosotros. Si sabe más sobre su caso de uso más común, puede optimizar el tamaño / forma del fragmento especificando una tupla de forma (por ejemplo, (2,2,2) en el ejemplo simple anterior). Esto le permite hacer lecturas a lo largo de un eje particular más eficiente u optimizar para lecturas / escrituras de cierto tamaño.

Comparación de rendimiento de E / S

Solo para enfatizar el punto, comparemos la lectura en segmentos de un conjunto de datos HDF5 fragmentado y una gran matriz 3D (~ 8GB) ordenada por Fortran que contiene los mismos datos exactos.

He borrado todos los cachés del sistema operativo entre cada ejecución, por lo que estamos viendo el rendimiento "en frío".

Para cada tipo de archivo, probaremos la lectura en una división x "completa" a lo largo del primer eje y una división z "completa" a lo largo del último eje. Para la matriz memmapped ordenada por Fortran, el segmento "x" es el peor de los casos, y el segmento "z" es el mejor de los casos.

El código utilizado está en una esencia (incluida la creación del archivo hdf ). No puedo compartir fácilmente los datos utilizados aquí, pero podría simularlos con una matriz de ceros de la misma forma ( 621, 4991, 2600) y escribir np.uint8 .

El chunked_hdf.py ve así:

import sys
import h5py

def main():
    data = read()

    if sys.argv[1] == ''x'':
        x_slice(data)
    elif sys.argv[1] == ''z'':
        z_slice(data)

def read():
    f = h5py.File(''/tmp/test.hdf5'', ''r'')
    return f[''seismic_volume'']

def z_slice(data):
    return data[:,:,0]

def x_slice(data):
    return data[0,:,:]

main()

memmapped_array.py es similar, pero tiene un toque más complejo para garantizar que los cortes se carguen realmente en la memoria (de forma predeterminada, se memmapped otra matriz memmapped , que no sería una comparación de manzanas con manzanas).

import numpy as np
import sys

def main():
    data = read()

    if sys.argv[1] == ''x'':
        x_slice(data)
    elif sys.argv[1] == ''z'':
        z_slice(data)

def read():
    big_binary_filename = ''/data/nankai/data/Volumes/kumdep01_flipY.3dv.vol''
    shape = 621, 4991, 2600
    header_len = 3072

    data = np.memmap(filename=big_binary_filename, mode=''r'', offset=header_len,
                     order=''F'', shape=shape, dtype=np.uint8)
    return data

def z_slice(data):
    dat = np.empty(data.shape[:2], dtype=data.dtype)
    dat[:] = data[:,:,0]
    return dat

def x_slice(data):
    dat = np.empty(data.shape[1:], dtype=data.dtype)
    dat[:] = data[0,:,:]
    return dat

main()

Veamos primero el rendimiento de HDF:

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python chunked_hdf.py z
python chunked_hdf.py z  0.64s user 0.28s system 3% cpu 23.800 total

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python chunked_hdf.py x
python chunked_hdf.py x  0.12s user 0.30s system 1% cpu 21.856 total

Una rebanada x "completa" y una rebanada z "completa" toman aproximadamente la misma cantidad de tiempo (~ 20 segundos). Teniendo en cuenta que se trata de una matriz de 8 GB, no está tan mal. La mayor parte del tiempo

Y si comparamos esto con los tiempos de matriz memmapped (está ordenado por Fortran: un "corte en z" es el mejor caso y un "corte en x" es el peor de los casos):

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python memmapped_array.py z
python memmapped_array.py z  0.07s user 0.04s system 28% cpu 0.385 total

jofer at cornbread in ~ 
$ sudo ./clear_cache.sh

jofer at cornbread in ~ 
$ time python memmapped_array.py x
python memmapped_array.py x  2.46s user 37.24s system 0% cpu 3:35:26.85 total

Sí, lo leiste bien. 0.3 segundos para una dirección de corte y ~ 3.5 horas para la otra.

¡El tiempo para cortar en la dirección "x" es mucho más largo que el tiempo que llevaría cargar toda la matriz de 8GB en la memoria y seleccionar la porción que queríamos! (Nuevamente, esta es una matriz ordenada por Fortran. El tiempo de corte x / z opuesto sería el caso de una matriz ordenada en C).

Sin embargo, si siempre queremos tomar un corte en la dirección del mejor de los casos, la gran matriz binaria en el disco es muy buena. (~ 0.3 segundos!)

Con una matriz memmapped, está atrapado con esta discrepancia de E / S (o quizás la anisotropía es un término mejor). Sin embargo, con un conjunto de datos HDF fragmentado, puede elegir el tamaño de fragmento de manera que el acceso sea igual o esté optimizado para un caso de uso particular. Te da mucha más flexibilidad.

En resumen

Con suerte, eso ayudará a aclarar una parte de su pregunta, en cualquier caso. HDF5 tiene muchas otras ventajas sobre los mapas de memoria "en bruto", pero no tengo espacio para expandirlos todos aquí. La compresión puede acelerar algunas cosas (los datos con los que trabajo no se benefician mucho de la compresión, por lo que rara vez la uso), y el almacenamiento en caché a nivel del sistema operativo a menudo se reproduce mejor con los archivos HDF5 que con los memmaps "en bruto". Más allá de eso, HDF5 es un formato contenedor realmente fantástico. Le brinda mucha flexibilidad en la administración de sus datos y puede usarse desde más o menos cualquier lenguaje de programación.

En general, pruébelo y vea si funciona bien para su caso de uso. Creo que te sorprenderás.