首页 文章

使用HDF5进行大型阵列存储（而不是平面二进制文件）是否存在分析速度或内存使用优势？

提问于
浏览 次
66

我正在处理大型3D阵列，我经常需要以各种方式进行切片以进行各种数据分析 . 典型的“立方体”可以是~100GB（将来可能会变大）

似乎python中大型数据集的典型推荐文件格式是使用HDF5（h5py或pytables） . 我的问题是：使用HDF5存储和分析这些立方体而不是将它们存储在简单的平面二进制文件中是否有任何速度或内存使用效益？ HDF5是否更适合表格数据，而不像我正在使用的大型数组？我看到HDF5可以提供很好的压缩，但我对处理速度和处理内存溢出更感兴趣 .

我经常只想分析立方体的一个大的子集 . pytables和h5py的一个缺点是，当我拿一个数组时，我总是得到一个numpy数组，用尽内存 . 但是，如果我切片平面二进制文件的numpy memmap，我可以得到一个视图，它将数据保存在磁盘上 . 因此，似乎我可以更轻松地分析我的数据的特定部分而不会超出我的记忆 .

我已经探讨了pytables和h5py，到目前为止还没有看到我的目的的好处 .

python numpy hdf5 pytables h5py

1 回答

  • 122

    HDF5优点：组织，灵活性，互操作性

    HDF5的一些主要优点是其分层结构（类似于文件夹/文件），与每个项目一起存储的可选任意元数据，以及其灵活性（例如压缩） . 这种组织结构和元数据存储可能听起来微不足道，但在实践中它非常有用 .

    HDF的另一个优点是数据集可以是固定大小或灵活大小 . 因此，无需创建整个新副本，就可以轻松地将数据附加到大型数据集 .

    此外，HDF5是一种标准化格式，其库可用于几乎所有语言，因此使用HDF可以轻松地在Matlab，Fortran，R，C和Python之间共享磁盘数据 . （公平地说，只要你知道C和F的顺序并知道存储数组的形状，类型等，它对于大二进制数组也不会太难 . ）

    HDF对大型阵列的优势：任意切片的I / O更快

    Just as the TL/DR: 对于~8GB的3D阵列，沿着任意轴读取"full"切片，使用分块的HDF5数据集需要大约20秒，对于相同数据的memmapped阵列，需要0.3秒（最佳情况）到超过3小时（最差情况） .

    除了上面列出的内容之外，还有一个很大的优势就是“分块”*磁盘数据格式，如HDF5：读取任意切片（强调任意）通常要快得多，因为磁盘上的数据更加连续平均 .

    * （HDF5不需要它 . 事实上，如果我没记错的话，在 h5py 中创建数据集的默认设置不是块 . ）

    基本上，对于给定切片的数据集，最佳情况下的磁盘读取速度和最坏情况下的磁盘读取速度将与分块的HDF数据集相当接近（假设您选择了合理的块大小或让库为您选择一个） . 使用简单的二进制数组，最好的情况是更快，但最坏的情况更糟 .

    有一点需要注意，如果您有SSD，您可能不会注意到读/写速度的巨大差异 . 但是，对于常规硬盘驱动器，顺序读取比随机读取快得多 . （即常规硬盘驱动器有很长的时间 . ）HDF在SSD上仍然具有优势，但更多的是由于其他功能（例如元数据，组织等）而不是原始速度 .

    首先，要清除混淆，访问 h5py 数据集会返回一个与numpy数组非常相似的对象，但在切片之前不会将数据加载到内存中 . （与memmap类似，但不完全相同 . ）有关详细信息，请查看h5py introduction .

    切片数据集会将数据的一个子集加载到内存中，但可能你想用它做一些事情，此时你无论如何都需要它在内存中 .

    如果您确实想要进行核外计算，则可以使用 pandas 或 pytables 轻松地获取表格数据 . 有可能使用 h5py （对于大型N-D阵列更好），但是您需要下拉到较低级别并自己处理迭代 .

    然而，类似numpy的核外计算的未来是Blaze . Have a look at it如果你真的想要那个路线 .

    “没有问题”的案例

    首先，考虑写入磁盘的3D C排序数组（我将通过调用 arr.ravel() 并打印结果来模拟它，以使事物更加明显）：

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

    这些值将按顺序存储在磁盘上，如下面第4行所示 . （暂时忽略文件系统细节和碎片 . ）

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

    在最好的情况下，让我们沿第一轴切片 . 请注意，这些只是数组的前36个值 . 这将是一个非常快速的阅读！ （一个寻求，一个阅读）

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

    类似地，沿第一轴的下一个切片将仅是接下来的36个值 . 要沿此轴读取完整切片，我们只需要一个 seek 操作 . 如果我们要读的是沿着这个轴的各种切片，那么这就是完美的文件结构 .

    但是，让我们考虑最坏情况：沿最后一个轴的切片 .

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

    要读取此片段，我们需要36次搜索和36次读取，因为所有值都在磁盘上分开 . 它们都不相邻！

    这可能看起来很小，但随着我们越来越大的数组， seek 操作的数量和大小迅速增长 . 对于以这种方式存储并通过 memmap 读取的大型（~10Gb）3D阵列，即使使用现代硬件，沿"worst"轴读取整个切片也可能需要数十分钟 . 同时，沿最佳轴的切片可能需要不到一秒钟 . 为简单起见，我只在单个轴上显示"full"个切片，但是对于任何数据子集的任意切片都会发生完全相同的事情 .

    顺便提一下，有几种文件格式可以利用这一点，并且基本上在磁盘上存储三个巨大的3D阵列副本：一个是C顺序，一个是F顺序，一个是中间的两个 . （一个例子是Geoprobe 's D3D format, though I'我不确定它是否记录在任何地方 . ）谁在乎最终文件大小是4TB，存储是否便宜！关于这一点的疯狂之处在于，因为主要用例是在每个方向上提取单个子切片，所以要进行的读取非常非常快 . 它运作得很好！

    简单的“分块”案例

    假设我们将3D阵列的2x2x2“块”存储为磁盘上的连续块 . 换句话说，像：

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

    所以磁盘上的数据看起来像 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])

    只是为了表明它们是 arr 的2x2x2块，请注意这些是 chunked 的前8个值：

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

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

    要读入沿轴的任何切片，我们将读取6或9个连续的块（数据量是我们需要的两倍），然后只保留我们想要的部分 . 这是最糟糕的情况下最多9次寻求，而最多只寻求36次寻找非分块版本 . （但对于memmapped数组，最好的情况仍然是6次搜索与1次 . ）因为顺序读取与搜索相比非常快，这大大减少了将任意子集读入内存所需的时间 . 再次，对于更大的阵列，这种效果会变得更大 .

    HDF5更进一步 . 块不必连续存储，并且它们由B树索引 . 此外，它们不必在磁盘上具有相同的大小，因此可以对每个块应用压缩 .

    带有h5py的分块数组

    默认情况下， h5py 不会在磁盘上创建分块的HDF文件（相反，我认为 pytables ） . 但是，如果在创建数据集时指定 chunks=True ，则会在磁盘上获得一个分块数组 .

    作为一个快速，最小的例子：

    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?'

    请注意 chunks=True 告诉 h5py 为我们自动选择一个块大小 . 如果您对最常见的用例有更多了解，可以通过指定形状元组来优化块大小/形状（例如，在上面的简单示例中为 (2,2,2) ） . 这使您可以更有效地沿特定轴读取或优化特定大小的读/写 .

    I / O性能比较

    为了强调这一点，让我们比较分块中的读数，从分块的HDF5数据集和包含相同精确数据的大型（~8GB）Fortran排序的3D阵列 .

    我在每次运行之间都有cleared all OS caches，所以我们看到了"cold"的性能 .

    对于每种文件类型，我们将测试沿第一轴的“完整”x切片和沿最后一轴的“完整”z-slize的读数 . 对于Fortran排序的memmapped数组，“x”切片是最坏的情况，“z”切片是最好的情况 .

    使用的代码是in a gist（包括创建 hdf 文件） . 我不能轻松共享此处使用的数据，但您可以通过相同形状的零数组（ 621, 4991, 2600) 和类型 np.uint8 ）来模拟它 .

    chunked_hdf.py 看起来像这样：

    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 是类似的，但触摸更复杂，以确保切片实际上加载到内存中（默认情况下，将返回另一个 memmapped 数组，这不是一个苹果对苹果的比较） .

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

    我们先来看看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

    “完整”x切片和“完整”z切片花费大约相同的时间量（~20秒） . 考虑到这是一个8GB阵列，这还不错 . 大多数时候

    如果我们将它与memmapped数组时间进行比较（它是Fortran有序的：“z-slice”是最好的情况，“x-slice”是最坏的情况 . ）：

    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

    是的，你没有看错 . 一个切片方向为0.3秒，另一个为3.5小时 .

    在"x"方向切片的时间远远长于将整个8GB阵列加载到内存并选择我们想要的切片所需的时间！ （同样，这是一个Fortran排序的数组 . 相反的x / z切片时序就是C有序数组的情况 . ）

    但是，如果我们总是希望沿着最好的方向进行切片，那么磁盘上的大二进制数组非常好 . （~0.3秒！）

    使用memmapped数组，你会遇到这种I / O差异（或者各向异性可能是一个更好的术语） . 但是，对于分块的HDF数据集，您可以选择chunksize，以使访问相等或针对特定用例进行优化 . 它为您提供了更多的灵活性 .

    总结

    希望无论如何，这有助于澄清问题的一部分 . 与“原始”memmaps相比，HDF5还有许多其他优点，但我没有足够的空间在这里展开所有这些 . 压缩可以加速一些事情（我使用的数据不会从压缩中受益很多，所以我很少使用它），并且操作系统级缓存通常比HDF5文件更好地使用“原始”memmaps . 除此之外，HDF5是一种非常出色的容器格式 . 它为您提供了很多管理数据的灵活性，并且可以或多或少地使用任何编程语言 .

    总的来说，尝试一下，看看它是否适用于您的用例 . 我想你可能会感到惊讶 .

    回复于

相关问题