NumPy和Matlab性能差异

27
我正在为稀疏自编码器计算反向传播算法。我已经使用numpy和matlab分别实现了它的python代码。两份代码几乎相同,但性能却有很大差异。matlab完成任务所需的时间是0.252454秒,而numpy需要0.973672151566秒,这几乎是四倍之多。由于我将在后面的最小化问题中多次调用此代码,因此这种差异会导致实现之间数分钟的延迟。这是正常的行为吗?我该如何改善numpy的性能?
Numpy实现:
Sparse.rho是一个调整参数,sparse.nodes是隐藏层中节点的数量(25),sparse.input(64)是输入层中节点的数量,theta1和theta2分别是第一层和第二层的权重矩阵,维度为25x64和64x25,m等于10000,rhoest的尺寸为(25,),x的尺寸为10000x64,a3的尺寸为10000x64,a2的尺寸为10000x25。
更新:我根据一些回答的想法对代码进行了更改。现在性能为numpy: 0.65 vs matlab: 0.25。
partial_j1 = np.zeros(sparse.theta1.shape)
partial_j2 = np.zeros(sparse.theta2.shape)
partial_b1 = np.zeros(sparse.b1.shape)
partial_b2 = np.zeros(sparse.b2.shape)
t = time.time()

delta3t = (-(x-a3)*a3*(1-a3)).T

for i in range(m):

    delta3 = delta3t[:,i:(i+1)]
    sum1 =  np.dot(sparse.theta2.T,delta3)
    delta2 = ( sum1 + sum2 ) * a2[i:(i+1),:].T* (1 - a2[i:(i+1),:].T)
    partial_j1 += np.dot(delta2, a1[i:(i+1),:])
    partial_j2 += np.dot(delta3, a2[i:(i+1),:])
    partial_b1 += delta2
    partial_b2 += delta3

print "Backprop time:", time.time() -t

Matlab实现:

tic
for i = 1:m

    delta3 = -(data(i,:)-a3(i,:)).*a3(i,:).*(1 - a3(i,:));
    delta3 = delta3.';
    sum1 =  W2.'*delta3;
    sum2 = beta*(-sparsityParam./rhoest + (1 - sparsityParam) ./ (1.0 - rhoest) );
    delta2 = ( sum1 + sum2 ) .* a2(i,:).' .* (1 - a2(i,:).');
    W1grad = W1grad + delta2* a1(i,:);
    W2grad = W2grad + delta3* a2(i,:);
    b1grad = b1grad + delta2;
    b2grad = b2grad + delta3;
end
toc

1
有一个名为mlabwrap的模块。您可以通过导入它将Matlab用作Python库。语法非常简单。您可以在此处找到源代码和详细文档。http://mlabwrap.sourceforge.net/ - Koustav Ghosal
5
看看cython。时间上的差异是“预期的”,因为MATLAB有JIT,而CPython没有。如果所有代码都是单个numpy调用,则时间将类似,但你所看到的可能是解释开销。使用cython编写扩展非常容易,通过在正确位置为变量添加一些类型,您可以获得巨大的收益。 - Bakuriu
“data”的形状是什么?具体来说,“m”与其他维度相比如何? - hpaulj
m = 10000,x 是一个 10000x64 的矩阵,theta1 是一个 25x64 的矩阵,theta2 是一个 64x25 的矩阵。 - pabaldonedo
如果您无法将x作为整个矩阵处理,则最好在小维度上循环而不是大维度上循环。但这可能需要一些创意。 - hpaulj
3个回答

56
"Matlab比NumPy快"或反之亦然是错误的说法。它们的性能通常是可比较的。使用NumPy时,要获得良好的性能,您必须记住NumPy的速度来自于调用用C/C++/Fortran编写的底层函数。当您将这些函数应用于整个数组时,它的表现很好。一般而言,在Python循环中调用这些NumPy函数处理较小的数组或标量时,性能较差。
那么Python循环有什么问题呢?每次通过Python循环迭代都会调用next方法。每次使用[]索引都会调用__getitem__方法。每次+=都会调用__iadd__。每次点属性查找(例如像np.dot)都涉及函数调用。这些函数调用会显著影响速度。这些钩子赋予了Python表达能力——例如,对于字符串进行索引与对于字典进行索引意义不同。相同的语法,不同的含义。通过给对象不同的__getitem__方法来实现这种魔法。"

但是这种表达能力的提升会以速度为代价。因此,当您不需要所有动态表达能力时,为了获得更好的性能,请尽量限制自己只使用整个数组的NumPy函数调用。

因此,去掉for循环;尽可能使用“向量化”方程式。例如,而不是

for i in range(m):
    delta3 = -(x[i,:]-a3[i,:])*a3[i,:]* (1 - a3[i,:])    

你可以一次性计算每个delta3

delta3 = -(x-a3)*a3*(1-a3)

for 循环 中,delta3 是一个向量,而使用向量化的方程式,delta3 就是一个矩阵。


for循环中,部分计算与 i 无关,因此应该在循环外进行。例如,sum2 看起来像一个常数:

sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest) )

这里有一个可运行的示例,其中包含您代码(orig)的替代实现(alt)。我的时间基准测试显示速度提高了6.8倍
In [52]: %timeit orig()
1 loops, best of 3: 495 ms per loop

In [53]: %timeit alt()
10 loops, best of 3: 72.6 ms per loop

import numpy as np


class Bunch(object):
    """ http://code.activestate.com/recipes/52308 """
    def __init__(self, **kwds):
        self.__dict__.update(kwds)

m, n, p = 10 ** 4, 64, 25

sparse = Bunch(
    theta1=np.random.random((p, n)),
    theta2=np.random.random((n, p)),
    b1=np.random.random((p, 1)),
    b2=np.random.random((n, 1)),
)

x = np.random.random((m, n))
a3 = np.random.random((m, n))
a2 = np.random.random((m, p))
a1 = np.random.random((m, n))
sum2 = np.random.random((p, ))
sum2 = sum2[:, np.newaxis]

def orig():
    partial_j1 = np.zeros(sparse.theta1.shape)
    partial_j2 = np.zeros(sparse.theta2.shape)
    partial_b1 = np.zeros(sparse.b1.shape)
    partial_b2 = np.zeros(sparse.b2.shape)
    delta3t = (-(x - a3) * a3 * (1 - a3)).T
    for i in range(m):
        delta3 = delta3t[:, i:(i + 1)]
        sum1 = np.dot(sparse.theta2.T, delta3)
        delta2 = (sum1 + sum2) * a2[i:(i + 1), :].T * (1 - a2[i:(i + 1), :].T)
        partial_j1 += np.dot(delta2, a1[i:(i + 1), :])
        partial_j2 += np.dot(delta3, a2[i:(i + 1), :])
        partial_b1 += delta2
        partial_b2 += delta3
        # delta3: (64, 1)
        # sum1: (25, 1)
        # delta2: (25, 1)
        # a1[i:(i+1),:]: (1, 64)
        # partial_j1: (25, 64)
        # partial_j2: (64, 25)
        # partial_b1: (25, 1)
        # partial_b2: (64, 1)
        # a2[i:(i+1),:]: (1, 25)
    return partial_j1, partial_j2, partial_b1, partial_b2


def alt():
    delta3 = (-(x - a3) * a3 * (1 - a3)).T
    sum1 = np.dot(sparse.theta2.T, delta3)
    delta2 = (sum1 + sum2) * a2.T * (1 - a2.T)
    # delta3: (64, 10000)
    # sum1: (25, 10000)
    # delta2: (25, 10000)
    # a1: (10000, 64)
    # a2: (10000, 25)
    partial_j1 = np.dot(delta2, a1)
    partial_j2 = np.dot(delta3, a2)
    partial_b1 = delta2.sum(axis=1)
    partial_b2 = delta3.sum(axis=1)
    return partial_j1, partial_j2, partial_b1, partial_b2

answer = orig()
result = alt()
for a, r in zip(answer, result):
    try:
        assert np.allclose(np.squeeze(a), r)
    except AssertionError:
        print(a.shape)
        print(r.shape)
        raise

提示:注意我保留了所有中间数组的形状注释。知道数组的形状有助于我理解你的代码在做什么。数组的形状可以帮助指导你使用正确的NumPy函数。或者至少,关注形状可以帮助你知道操作是否合理。例如,当你计算时。

np.dot(A, B)

如果 A.shape = (n, m)B.shape = (m, p),那么 np.dot(A, B) 将是一个形状为 (n, p) 的数组。

如果使用np.dot,构建C_CONTIGUOUS顺序的数组可能会有多达3倍的加速效果:

下面,xxf相同,只是x是C_CONTIGUOUS而xf是F_CONTIGUOUS——yyf也是同样的关系。

import numpy as np

m, n, p = 10 ** 4, 64, 25
x = np.random.random((n, m))
xf = np.asarray(x, order='F')

y = np.random.random((m, n))
yf = np.asarray(y, order='F')

assert np.allclose(x, xf)
assert np.allclose(y, yf)
assert np.allclose(np.dot(x, y), np.dot(xf, y))
assert np.allclose(np.dot(x, y), np.dot(xf, yf))

%timeit 基准测试显示速度差异:

In [50]: %timeit np.dot(x, y)
100 loops, best of 3: 12.9 ms per loop

In [51]: %timeit np.dot(xf, y)
10 loops, best of 3: 27.7 ms per loop

In [56]: %timeit np.dot(x, yf)
10 loops, best of 3: 21.8 ms per loop

In [53]: %timeit np.dot(xf, yf)
10 loops, best of 3: 33.3 ms per loop

关于Python中的基准测试:

仅使用一对time.time()调用的差异来测试Python代码的速度可能会误导。您需要多次重复测量。最好禁用自动垃圾收集器。同时,重要的是测量大时间跨度(例如至少10秒的重复),以避免由于时钟计时器的低分辨率而产生错误,并降低time.time调用开销的重要性。Python提供了timeit模块,可以代替自己编写所有这些代码。我基本上使用它来计算代码的运行时间,只不过为了方便,我通过一个IPython终端进行调用。

我不确定这是否会影响你的基准测试,但请注意它可能会产生差异。根据我提供的链接到的问题,根据time.time,两个代码片段的差异因子为1.7倍,而使用timeit进行基准测试显示,这两个代码片段运行时间基本相同。


预先计算“for-loop”之前的“delta3”,并将“sum2”移到外部有所帮助(我已经更新了问题),但它仍然比Matlab慢两倍以上。令我印象深刻的是,Matlab在for循环中计算“delta3”的时间几乎与numpy访问现在作为矩阵的预计算delta3的行的时间相同。这总是numpy比matlab慢吗? - pabaldonedo
我已经更改了求和的方法,如下所示:sum2 = np.dot(sum2.reshape(-1,1),np.ones((1,sum1.shape[1])))。现在它可以正常工作了,是否有更好的方法来实现呢?非常感谢您的回复。 - pabaldonedo
1
你可以使用 sum2 = sum2[:, np.newaxis]sum2 从形状为 (25,) 的数组转换为形状为 (25,1) 的数组。NumPy广播会负责将其升级到形状为 (25,10000),而不会消耗重复相同值10000次的不必要内存。在我的电脑上,sum2[:, np.newaxis]np.dot(sum2.reshape(-1,1),np.ones((1,sum1.shape[1]))) 快了约4300倍。当然,我们只需要执行一次,因此总速度提升微不足道。但是,这是一个很好的技巧。 - unutbu
我稍微改进了NumPy的结果,将行sum1 = (sparse.theta2.T, delta3)更改为sum1 = scipy.linalg.blas.dgemm(alpha=1.0, a=sparse.theta2.T, b=delta3, trans_b=False),从而将时间缩短到了0.10s,因为这两个矩阵都是Fortran连续的。如果我重新设计它,使得所有操作(或至少大部分操作)都在Fortran连续的矩阵之间或C连续的矩阵之间进行,并应用最方便的点积运算符,那么我能否进一步优化它,或者改进会很微不足道? - pabaldonedo
1
@hpauljпјҡжІЎй”ҷпјҢдҪҶжҳҜpabaldonedoд»ҺдёҖдёӘеҪўзҠ¶дёә(25,)зҡ„ж•°з»„ејҖе§ӢгҖӮд»–йңҖиҰҒдёҖз§Қж–№жі•е°Ҷе…¶йҮҚж–°еЎ‘йҖ дёә(25,1)гҖӮnp.random.random((p,))еҸӘжҳҜжҲ‘еҲ¶дҪңзҡ„дёҖдёӘж•°з»„пјҢз”ЁдҪңд»–зңҹе®һж•°з»„зҡ„жӣҝд»Је“ҒгҖӮ - unutbu
显示剩余6条评论

3
我会从原地操作开始,以避免每次分配新的数组:
partial_j1 += np.dot(delta2, a1[i,:].reshape(1,a1.shape[1]))
partial_j2 += np.dot(delta3, a2[i,:].reshape(1,a2.shape[1]))
partial_b1 += delta2
partial_b2 += delta3

你可以替换这个表达式:
a1[i,:].reshape(1,a1.shape[1])

使用更简单、更快速(感谢Bi Rico)的方法:

a1[i:i+1]

此外,这一行:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest))

每次循环似乎都是相同的,您不需要重新计算。

而且,可能是一个小优化,您可以将所有的x[i,:]替换为x[i]

最后,如果您能够承担m倍的内存分配,可以按照unutbu的建议对循环进行向量化:

for m in range(m):
    delta3 = -(x[i]-a3[i])*a3[i]* (1 - a3[i])

使用:

delta3 = -(x-a3)*a3*(1-a3)

你可以随时使用Numba,无需矢量化(也无需使用更多内存),就能显著提高速度。


1
我已经检查过了,原地操作几乎没有什么区别。 - pabaldonedo
3
a1[i,:].reshape(1,a1.shape[1])可以写成a[i:i+1] - Bi Rico
毕里克,我不这么认为。 - user2304916

2
numpy和matlab之间的性能差异一直让我感到沮丧。它们通常归结于底层的lapack库。据我所知,matlab默认使用完整的atlas lapack,而numpy使用lapack light。Matlab认为人们不关心空间和体积,而numpy认为人们关心。类似问题有一个很好的答案。

3
在这种情况下,我几乎不相信是 LAPACK 的错,因为他们只使用点积。更有可能的是 MATLAB 对循环进行了即时编译以加快速度。 - user2304916
我的经验是,numpy的运行速度与旧版Matlab或Octave大致相同(或最多只慢一半)。但是新版本的Matlab似乎更加积极地进行向量化或编译(jit)。对于有经验的“旧”Matlab用户来说,for i = 1:ma3(i,:)是缓慢的代码标志。 - hpaulj
FWIW,MATLAB 目前已经停止使用 ATLAS,转而使用英特尔 MKL(大约从 v7 开始,也就是十多年前)。你也可以编译 NumPy 以便支持 MKL。Christoph Gohlke 提供了 NumPy-MKL 的 Windows 二进制文件:http://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy - Amro
是的,我同意这更可能是jit的一个特性。引入numpypy可以提高这种速度吗?Matlab的jit非常厉害,能够找到语法相似的Matlab例程并调用预编译的C代码块。如果你像在C中编码一样编写Matlab代码,它的运行速度就和真正的C代码一样快,因为它已经在预编译的C上运行了。 - Philliproso

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接