【问题标题】:Difference on performance between numpy and matlabnumpy 和 matlab 之间的性能差异
【发布时间】:2013-09-02 05:18:13
【问题描述】:

我正在为稀疏自动编码器计算 backpropagation 算法。我已经使用numpymatlab 在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。

UPDATE:我在响应的一些想法之后对代码进行了更改。性能现在是 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

【问题讨论】:

  • 有一个名为 mlabwrap 的模块。您可以通过导入将 matlab 用作 python 库。语法非常简单。您可以在此处找到源代码和详细文档。mlabwrap.sourceforge.net
  • 看看cython。时间上的差异是预期的,因为 MATLAB 有 JIT,而 CPython 没有。如果所有代码都是单个 numpy 调用,那么时间会相似,但您看到的可能是解释开销。使用 cython 编写扩展非常简单,您可能会在正确的位置向变量添加一些类型,从而获得巨大的收益。
  • data 的形状是什么?具体m与其他维度相比如何?
  • m = 10000,x 为 10000x64 矩阵,theta1 为 25x64 矩阵,theta2 为 64x25。
  • 如果您不能将x 作为一个整体矩阵工作,最好在小维度上循环而不是在大维度上循环。但这可能需要一些独创性。

标签: python performance matlab numpy backpropagation


【解决方案1】:

我会从就地操作开始,以避免每次都分配新数组:

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 并显着提高速度,而无需矢量化(也无需使用更多内存)。

【讨论】:

  • 我已经检查过了,inplace 操作几乎没有区别。
  • a1[i,:].reshape(1,a1.shape[1]) 可以写成a[i:i+1]
  • Bi Rico,我不这么认为。
【解决方案2】:

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

【讨论】:

  • 在这种情况下,我很难相信 LAPACK 应该受到责备,因为他们只使用点积。更可能 MATLAB 做了一些 jit 来加速循环。
  • 我的经验是,numpy 的运行速度与旧的 Matlab 或 Octave 的运行速度大致相同(或最差的一半)。但新的 Matlab 版本似乎更积极地进行矢量化或编译(jit)。对于有“旧”Matlab 经验的人,for i = 1:ma3(i,:) 是慢代码标志。
  • fwiw,MATLAB 暂时停止使用 ATLAS 支持 Intel MKL(我认为是从 v7 开始,那是 10 多年前)。您也可以针对 MKL 编译 NumPy。 Christoph Gohlke 提供 NumPy-MKL windows 二进制文件:lfd.uci.edu/~gohlke/pythonlibs/#numpy
  • 是的,这更有可能是我同意的 jit 的一个功能。引入 numpypy 可以提高这个速度吗? Matlabs jit 非常神奇地发现语法相似的 matlab 例程并调用预编译的 C 代码位。如果你在 matlab 中编码,就像你在 C 中编码一样快,因为它已经在运行预编译的 C。
【解决方案3】:

说“Matlab 总是比 NumPy 快”是错误的,反之亦然 反之亦然。他们的表现通常是可比的。使用 NumPy 时,要变得更好 性能你必须牢记 NumPy 的速度来自于调用 用 C/C++/Fortran 编写的底层函数。申请时效果很好 这些函数到整个数组。一般来说,在 Python 循环中对较小的数组或标量调用这些 NumPy 函数时,性能会较差。

你问的 Python 循环有什么问题?通过 Python 循环的每次迭代都是 调用next 方法。每次使用[] 索引都是调用 __getitem__ 方法。每个+= 都是对__iadd__ 的调用。每个虚线属性 查找(例如np.dot)涉及函数调用。那些函数调用 加起来是速度的重大障碍。这些钩子给了 Python 表现力——字符串的索引意味着与索引不同的东西 例如,对于 dicts。相同的语法,不同的含义。魔术是通过为对象提供不同的__getitem__ 方法来实现的。

但这种表现力是以速度为代价的。所以当你不需要所有 那种动态的表现力,为了得到更好的表现,尽量把自己限制在 NumPy 函数调用整个数组。

所以,删除 for 循环;尽可能使用“矢量化”方程。例如,而不是

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

您可以同时为每个i 计算delta3

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) )

这是一个可运行的示例,其中包含您的代码 (orig) 的替代实现 (alt)。

我的 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

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

np.dot(A, B)

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


它可以帮助以 C_CONTIGUOUS 顺序构建数组(至少,如果使用 np.dot)。这样做可能会提高 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 中的基准测试:

It can be misleading 使用 time.time() 调用对的差异来对 Python 中的代码速度进行基准测试。 您需要多次重复测量。最好禁用自动垃圾收集器。测量大的时间跨度(例如至少 10 秒的重复次数)也很重要,以避免由于时钟计时器分辨率差而导致的错误,并减少 time.time 调用开销的重要性。 Python 不是自己编写所有代码,而是为您提供timeit module。我基本上是用它来计时代码片段,只是为了方便起见,我通过IPython terminal 调用它。

我不确定这是否会影响您的基准测试,但请注意它可能会有所作为。在question I linked to 中,根据time.time,两段代码相差1.7 倍,而使用timeit 的基准测试显示这两条代码运行的时间基本相同。

【讨论】:

  • precomputing delta3before the for-loopand take sum2 outside help(我已经更新了问题)但它仍然比matlab慢两倍多。让我印象深刻的是,matlab 在 for 循环内计算 delta3 所需的时间几乎与 numpy 访问预先计算的 delta3 行作为矩阵所需的时间相同。与matlab 相比,这总是numpy 这么慢吗?
  • 感谢您的彻底回复,但操作 sum1+sum2 在我的计算机中崩溃,sum1 的尺寸为 25,10000sum2 的尺寸为(25,)
  • 我已经更改了总和,添加了上一行,如下sum2 = np.dot(sum2.reshape(-1,1),np.ones((1,sum1.shape[1])))。现在它可以工作了,有没有更好的方法呢?非常感谢您的回复。
  • 您可以使用sum2 = sum2[:, np.newaxis]sum2 从形状数组 (25,) 转换为形状数组 (25,1)。 NumPy broadcasting 将负责将其“升级”为形状 (25, 10000) 而不 消耗不必要的内存,重复相同的值 10000 次。 sum2[:, np.newaxis] 在我的电脑上比 np.dot(sum2.reshape(-1,1),np.ones((1,sum1.shape[1]))) 快大约 4300 倍。当然,我们只这样做一次,因此总的速度增益可以忽略不计。不过,这是一个很好的技巧。
  • @hpaulj:没错,但 pabaldonedo 是从一个形状数组 (25, ) 开始的。他需要一种方法将其重塑为(25, 1)np.random.random((p, )) 只是我制作的一个数组,用来代替他的真实数组。