为什么使用 numpy.save 从数组中保存切片会根据切片的方向而变慢？

Question

我有 3D 数据量（x、y、z），我想保存 2D 切片（xy、yz、xz 平面）并保存它们以备将来使用。

我尝试这样做的方法是让一个函数 (slice_data) 获取切片，另一个函数 (save_slices) 调用 slice_data，然后使用 numpy.save 保存切片。

如果我不保存切片，则无论我是否提取 xy、yz、xz 平面，获取切片的时间都是相似的。 但是，如果我保存切片，保存切片的时间取决于切片的方向，并且对于 xy、yz、xz 平面是不同的

这是为什么呢？使用我的完整数据，这种差异从几分钟到几小时......

我在下面写了一个虚拟代码，类似于我在完整数据集中用来演示问题的代码。时间差如下图：

不保存，只是切片

mean time x-slice: 0.00011536836624145507 sec
mean time y-slice: 0.00011417627334594726 sec
mean time z-slice: 0.00011371374130249023 sec

切片和保存：

mean time x-slice: 0.04629791975021362 sec
mean time y-slice: 0.06096100091934204 sec
mean time z-slice: 0.08996494293212891 sec

代码：

import os
import numpy as np
import time
import matplotlib.pyplot as plt

# take a slice of the data
def slice_data(roi):
    dic = {}
    data = np.zeros((512,512,256))
    dic['data'] = np.squeeze( data[roi[0]:roi[1]+1, roi[2]:roi[3]+1, roi[4]:roi[5]+1] )
    return dic


# save slices if the data
def save_slices(roi, save=False):
    var = 'data'
    for i in range(0,6):
                # iterate to simulate a time series of data
        a = slice_data(roi)[var]
        var_dir = 'save_test/'
        if not os.path.exists(var_dir): os.makedirs(var_dir)
        file = var_dir + '{0:04d}{1}'.format(i,'.npy')

        if save is True:
            np.save(file, a)


## define slices
roix=[256, 256, 0, 512, 0, 256] # yz plane slice
roiy=[0, 512, 256, 256, 0, 256] # xz plane slice
roiz=[0, 512, 0, 512, 128, 128] # xy plane slice

## Calculate slices and do not save the results
dtx = []
dty = []
dtz = []
for i in range(100):
    time0 = time.time()
    save_slices(roix)
    time1 = time.time()
    dtx.append(time1-time0)

    time0 = time.time()
    save_slices(roiy)
    time1 = time.time()
    dty.append(time1-time0)


    time0 = time.time()
    save_slices(roiz)
    time1 = time.time()
    dtz.append(time1-time0)

plt.figure(1)
plt.plot(dtx)
plt.plot(dty)
plt.plot(dtz)
plt.title('time to run code without saving data')

print('mean time x-slice: {} sec'.format(np.mean(dtx)))
print('mean time y-slice: {} sec'.format(np.mean(dty)))
print('mean time z-slice: {} sec'.format(np.mean(dtz)))


## Calculate slices and do save the results
dtx = []
dty = []
dtz = []
for i in range(100):
    time0 = time.time()
    save_slices(roix, save=True)
    time1 = time.time()
    dtx.append(time1-time0)

    time0 = time.time()
    save_slices(roiy, save=True)
    time1 = time.time()
    dty.append(time1-time0)


    time0 = time.time()
    save_slices(roiz, save=True)
    time1 = time.time()
    dtz.append(time1-time0)

plt.figure(2)
plt.plot(dtx)
plt.plot(dty)
plt.plot(dtz)
plt.title('time to run code and save data')

print('mean time x-slice: {} sec'.format(np.mean(dtx)))
print('mean time y-slice: {} sec'.format(np.mean(dty)))
print('mean time z-slice: {} sec'.format(np.mean(dtz)))

Answer 1

简答

只有 roix 数组是 c_contiguous。因此从内存到 CPU 的总线传输比非连续数据更快（由于总线以块的形式移动数据并缓存它们）

通过将其设置为 C 连续 np.save(file, np.asarray(a, order='C'))

更多解释

分析

您应该使用 timeit 来计时您的表演，而不是自定义方法。

我已经为你做了一个例子：

在我们得到的一个单元格中：

import os
import numpy as np
import time
import matplotlib.pyplot as plt

# take a slice of the data
def slice_data(roi):
    dic = {}
    data = np.zeros((512,512,256))
    dic['data'] = np.squeeze( data[roi[0]:roi[1]+1, roi[2]:roi[3]+1, roi[4]:roi[5]+1] )
    return dic


# save slices if the data
def save_slices(roi, save=False):
    var = 'data'
    for i in range(0,6):
                # iterate to simulate a time series of data
        a = slice_data(roi)[var]
        var_dir = 'save_test/'
        if not os.path.exists(var_dir): os.makedirs(var_dir)
        file = var_dir + '{0:04d}{1}'.format(i,'.npy')

        if save is True:
            np.save(file, a)


## define slices
roix=[256, 256, 0, 512, 0, 256] # yz plane slice
roiy=[0, 512, 256, 256, 0, 256] # xz plane slice
roiz=[0, 512, 0, 512, 128, 128] # xy plane slice

在其他方面：

%%timeit -n 100
save_slices(roix) # 19.8 ms ± 285 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit -n 100
save_slices(roiy) # 20.5 ms ± 948 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit -n 100
save_slices(roiz) # 20 ms ± 345 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

随着存档

%%timeit -n 10 -r 3
save_slices(roix, True) # 32.7 ms ± 2.31 ms per loop (mean ± std. dev. of 3 runs, 10 loops each)

%%timeit -n 10 -r 3
save_slices(roiy, True) # 101 ms ± 2.61 ms per loop (mean ± std. dev. of 3 runs, 10 loops each)

%%timeit -n 10 -r 3
save_slices(roix, True) # 1.9 s ± 21.1 ms per loop (mean ± std. dev. of 3 runs, 10 loops each)

所以，正如您已经注意到的，没有保存，性能是一样的！ 让我们进入np.save() 方法

np.save 方法

np.save 负责处理 io 流，并调用 write_array 方法。这对于 C_contigous 数组来说真的很快。 （快速存取存储器）

让我们验证这个假设：

np.squeeze( np.zeros((512,512,256))[roix[0]:roix[1]+1, roix[2]:roix[3]+1, roix[4]:roix[5]+1] ).flags.c_contiguous # returns True

np.squeeze( np.zeros((512,512,256))[roiy[0]:roiy[1]+1, roiy[2]:roiy[3]+1, roiy[4]:roiy[5]+1] ).flags.c_contiguous # returns False

np.squeeze( np.zeros((512,512,256))[roiz[0]:roiz[1]+1, roiz[2]:roiz[3]+1, roiz[4]:roiz[5]+1] ).flags.c_contiguous # returns False

所以这可以解释roix 和roiy/roiz 之间的区别。

roiy 和roiz 之间差异的可能解释。数据传输使程序变慢

在那之后，我只能做出假设，roiz 似乎比roiy 更加分散。 write_array 方法需要花费大量时间。

我现在无法自己测试，但是这部分可以使用 linux 中的perf 命令进行验证。 （要查看使用总线的次数，例如缓存未命中的次数）。 如果我不得不做一个疯狂的猜测，我会说缓存未命中率非常高，因为数据不连续。因此，将数据从 RAM 传输到 CPU 确实会减慢处理速度。

其他处理存储的方式

我没有尝试过，但有一个很好的问题，有一些有用的答案：best way to preserve numpy arrays on disk

Answer 2

原因是 Numpy 默认按行主要顺序存储数据。如果你改变了

data = np.zeros((512,512,256))

到

# order F means column major
data = np.zeros((512,512,256), order='F')

您会发现保存 X 切片需要的时间最长。

如果您要保存 XY 平面的多个切片（当您更改为 Z 坐标时），通过转置数组并将其复制到 force a new memory layout.，您将看到更好的性能，这将确保内存布局适合您的访问模式，从而更快地读取（和保存）。下面有更详细的解释。

---

我们以下面的矩阵为例（来自Numpy glossary）：

m = [[1, 2, 3],
     [4, 5, 6]]

如果这在内存中以行主要顺序（numpy 术语中的 C 顺序）表示，则其布局如下：

[1, 2, 3, 4, 5, 6]

如果矩阵在内存中以列主要顺序（或 F 表示 Fortran 顺序）表示，则其布局如下：

[1, 4, 2, 5, 3, 6]

现在，如果您使用m[:,2] 对该数组进行索引，您将获得[3, 6]，而使用m[1,:]，您将获得[4, 5, 6]。如果您回顾内存布局，您会发现值[3, 6] 在列主要表示中是连续的，[4, 5, 6] 在行主要表示中是连续的。

当从数组中读取大量元素时（例如保存一个元素时），连续读取这些值在性能方面要好得多，因为这样可以利用 1-2 的 CPU 缓存数量级faster 比从内存中读取。