"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
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)
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倍的加速效果:
下面,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
进行基准测试显示,这两个代码片段运行时间基本相同。
cython
。时间上的差异是“预期的”,因为MATLAB有JIT,而CPython没有。如果所有代码都是单个numpy调用,则时间将类似,但你所看到的可能是解释开销。使用cython编写扩展非常容易,通过在正确位置为变量添加一些类型,您可以获得巨大的收益。 - Bakuriux
作为整个矩阵处理,则最好在小维度上循环而不是大维度上循环。但这可能需要一些创意。 - hpaulj