Numpy中的卷积比Matlab中的卷积慢吗?

Question

Matlab中的卷积似乎比Numpy中的卷积快两倍。

Python代码(在我的机器上需要19秒):

import numpy as np
from scipy import ndimage
import time
img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125
start_time = time.time()
ndimage.convolve(img,kernel,mode='constant')
print "Numpy execution took ", (time.time() - start_time), "seconds"

Matlab代码(在我的机器上需要8.7秒):

img = ones(512,512,512);
kernel = ones(5,5,5) / 125;
tic
convn(img, kernel, 'same');
toc

两者给出相同的结果。

是否有办法提高Numpy匹配或击败Matlab的性能在这里?

有趣的是，在许多输入大小下，运行时的这个因子或~2的差异是一致的。

Answer 1

你到底在做什么操作?如果你不需要一般的N-d卷积，ndimage提供了许多优化。

例如，您当前的操作:

img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125
result = ndimage.convolve(img, kernel)

等价于:

img = np.ones((512,512,512))
result = ndimage.uniform_filter(img, 5)

但是，在执行速度上有很大的差异:

In [22]: %timeit ndimage.convolve(img, kernel)
1 loops, best of 3: 25.9 s per loop
In [23]: %timeit ndimage.uniform_filter(img, 5)
1 loops, best of 3: 8.69 s per loop

差异是由于uniform_filter沿每个轴进行一维卷积，而不是一般的三维卷积造成的。

在内核是对称的情况下，您可以进行这些简化并获得显著的速度提升。

我不确定convn，但如果您的输入数据符合某些标准，matlab的函数通常会在幕后进行这些优化。Scipy更常见地使用每个函数一个算法，并期望用户知道在哪种情况下选择哪个算法。

---

你提到了一个"拉普拉斯高斯"滤波器。我总是把术语搞混。

我认为你想要的ndimage函数是scipy.ndimage.gaussian_laplace或scipy.ndimage.gaussian_filter与order=2(通过高斯函数的二阶导数进行过滤)。

无论如何，两者都将操作简化为每个轴上的1-d卷积，这应该会提供显着的(2-3x)加速。

Answer 2

没有答案;只是一个注释:

在比较性能之前，需要确保两个函数返回相同的结果:

如果Matlab的convn返回与Octave的convn相同的结果，则convn不同于ndimage.convolve:

octave> convn(ones(3,3), ones(2,2))
ans =
   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1
In [99]: ndimage.convolve(np.ones((3,3)), np.ones((2,2)))
Out[99]: 
array([[ 4.,  4.,  4.],
       [ 4.,  4.,  4.],
       [ 4.,  4.,  4.]])

ndimage.convolve有其他模式，'reflect','constant','nearest','mirror', 'wrap'，但这些都不匹配convn的默认("full")行为。

---

对于2D数组，scipy.signal.convolve2d比scipy.signal.convolve更快。

对于3D数组，scipy.signal.convolve似乎具有与convn(A,B)相同的行为:

octave> x = convn(ones(3,3,3), ones(2,2,2))
x =
ans(:,:,1) =
   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1
ans(:,:,2) =
   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2
ans(:,:,3) =
   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2
ans(:,:,4) =
   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1
In [109]: signal.convolve(np.ones((3,3,3), dtype='uint8'), np.ones((2,2,2), dtype='uint8'))
Out[109]: 
array([[[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]],
       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],
       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],
       [[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]]], dtype=uint8)

---

注意np.ones((n,m,p))默认创建一个float数组。Matlab ones(n,m,p)似乎创建了一个整数数组。为了进行良好的比较，您应该尝试将numpy数组的dtype与Matlab矩阵的类型相匹配。