numpy和matlab之间的性能差异

Question

我正在计算backpropagation稀疏自动编码器的算法。我已经使用numpy和在python中实现了它matlab。代码几乎相同，但是性能却大不相同。matlab完成任务所需的时间为0.252454秒，而numpy为0.973672151566，几乎是原来的四倍。在最小化问题中，我将在以后多次调用此代码，因此这种差异会导致实现之间的延迟几分钟。这是正常行为吗？如何提高numpy的性能？

numpy实现：

Sparse.rho是调整参数，sparse.nodes是隐藏层中的节点数（25），sparse.input（64）是输入层中的节点数，theta1和theta2是第一层和第二层的权重矩阵第二层的尺寸分别为25x64和64x25，m等于10000，最圆的尺寸为（25，），x的尺寸为10000x64，a3
10000x64和a2 10000x25。

UPDATE：我已经根据响应的一些想法对代码进行了更改。现在的性能很麻木：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

Answer 1

说“ Matlab总是比NumPy更快”是错误的，反之亦然。通常他们的表现是可比的。使用NumPy时，要获得良好的性能，必须记住NumPy的速度来自调用用C
/ C ++ /
Fortran编写的基础函数。当您将这些函数应用于整个数组时，它的性能很好。通常，在Python循环中的较小数组或标量上调用那些NumPy函数时，性能会变差。

您问的Python循环有什么问题？Python循环中的每次迭代都是对next方法的调用。[]索引的每次使用都是对__getitem__方法的调用
。每个+=都是对的呼唤__iadd__。每个点分属性查找（例如like中的np.dot）都涉及函数调用。这些函数调用加起来大大阻碍了速度。这些钩子赋予Python强大的表达能力-
例如，为字符串建立索引与为dict建立索引的含义有所不同。相同的语法，不同的含义。通过为对象提供不同的__getitem__方法来实现魔术。

但是，这种表达能力的代价是速度。因此，当您不需要所有动态表达时，为了获得更好的性能，请尝试将自己限制为对整个数组进行NumPy函数调用。

因此，删除for循环；尽可能使用“向量化”公式。例如，代替

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

您可以一次delta3为每个计算i：

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

而中的for-loop delta3是向量，使用向量化方程式delta3是矩阵。

运算中的某些运算for-loop不依赖于运算，i因此应提升到循环之外。例如，sum2看起来像一个常量：

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

这是一个可运行的示例，其中包含alt代码（orig）的替代实现（）。

我的timeit基准测试显示 速度提高 了 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)。

它可以帮助以C_CONTIGUOUS顺序构建数组（至少，如果使用np.dot）。这样做最多可以将速度提高3倍：

下面，x是一样的xf，除了x是C_CONTIGUOUS并且 xf是F_CONTIGUOUS -对于同样的关系y和yf。

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的基准测试表明，这两段代码的运行时间基本相同。