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

Answer 1

HDF5的优势：组织性、灵活性、互操作性

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

HDF的另一个优点是数据集可以是固定大小的或灵活大小的。因此，很容易将数据追加到大型数据集，而无需创建整个新副本。

此外，HDF5是一种标准格式，几乎任何语言都有可用的库，因此，在诸如Matlab、Fortran、R、C和Python之间共享磁盘数据非常容易。（公平地说，只要您知道C与F的顺序，知道存储数组的形状、数据类型等，使用大型二进制数组也不难。）

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

正如TL/DR:对于约8GB的3D数组，使用分块HDF5数据集沿任意轴读取“完整”切片需要约20秒，对于相同数据的memmapped数组，从0.3秒（最佳情况）到需要3小时（最坏情况）。

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

*（HDF5不一定是分块数据格式。它支持分块，但不需要它。事实上，在h5py中创建数据集的默认方法是不分块（如果我没记错的话）

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

一个警告，如果你有一个SSD，你可能不会注意到读/写速度的巨大差异。不过，对于普通硬盘，顺序读取要比随机读取快得多。（例如，普通硬盘有很长的seek时间。）HDF在SSD上仍然具有优势，但它更多地是由于其其他特性（如元数据、组织等）而不是由于原始速度。

---

首先，为了消除混淆，访问h5py数据集会返回一个行为相当类似于numpy数组的对象，但在对数据进行切片之前不会将其加载到内存中。（类似于memmap，但不完全相同）查看^{} introduction获取更多信息。

对数据集进行切片将把数据的一个子集加载到内存中，但您可能希望对其执行某些操作，在这一点上，您仍然需要将其放到内存中。

如果您确实想进行非核心计算，您可以非常容易地使用pandas或pytables进行表格数据。有了h5py（对于大的N-D数组来说更好）这是可能的，但是您需要降到一个较低的级别并自己处理迭代。

然而，类似于numpy的核外计算的未来是光明的。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个seeks和36个reads，因为所有的值在磁盘上都是分开的。它们都不相邻！

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

顺便说一下，有几种文件格式利用了这一点，它们基本上在磁盘上存储了三个大型3D数组的副本：一个是C顺序的，一个是F顺序的，还有一个在两者之间的中间位置。（这方面的一个例子是Geoprobe的D3D格式，虽然我不确定它在任何地方都有文档记录。）谁在乎最终的文件大小是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个seeks，而非分块版本最多36个seeks。（但最好的情况仍然是6个seeks与memmapped数组的1个seeks）因为顺序读取比seeks快得多，这大大减少了将任意子集读入内存所需的时间。再一次，这个效果随着数组的增大而变大。

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，所以我们看到了“冷”的表现。

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

使用的代码是in a gist（包括创建hdf文件）。I c公司很难共享这里使用的数据，但是可以通过一个由相同形状的零组成的数组（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.5hours。

在“x”方向切片的时间比将整个8GB数组加载到内存并选择所需切片所需的时间长得多！（同样，这是一个Fortran有序数组。相反的x/z切片计时适用于C顺序数组。）

但是，如果我们总是想沿着最佳情况的方向分一杯羹，那么磁盘上的大二进制数组就非常好。（~0.3秒！）

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

总结

希望这有助于澄清你的一部分问题，无论如何。与“原始”内存映射相比，HDF5还有许多其他优势，但我在这里没有足够的空间来扩展它们。压缩可以加快一些速度（我使用的数据从压缩中没有太多好处，所以我很少使用它），操作系统级缓存通常比“原始”memmaps更好地处理HDF5文件。除此之外，HDF5是一种非常棒的容器格式。它在管理数据方面给了您很大的灵活性，并且可以或多或少地从任何编程语言中使用。

总的来说，试试看它是否适合您的用例。我想你可能会感到惊讶。