更快的函数评估调用

Question

我有一个函数需要至少评估 1e6 次。它只需要两个数组作为输入。

import numpy as np
np.random.seed(1)

res = np.random.rand(0,1,1000)
f1 = lambda r, sig: -0.5*(np.sum(2*np.log(sig) + np.log(2*np.pi) + (r/sig)**2))    

目前，

%timeit f1(res,np.sqrt(res)) 

时钟在 26 微秒。有什么方法可以让我有更快的评估调用/执行时间？

任何方法都可以；例如不使用 lambda 函数，使用本机数学而不是 numpy，或者任何涉及使用在 python 中实现的更快的数学运算来操纵方程的技巧。

编辑：

我犯了一个愚蠢的错误，即生成了一个空的 res 数组。所以应该是

res=np.random.rand(1000)

另外，第二个数组np.sqrt(res) 实际上独立于res，通常包含比res 小得多的值，但我使用sqrt(res) 只是为了说明。

所以重写代码更清楚

import numpy as np
np.random.seed(1)

res1 = np.random.rand(1000)
res2 = np.sqrt(res) #not true in general but res2<res1
f1 = lambda r, sig: -0.5*(np.sum(2*np.log(sig) + np.log(2*np.pi) + (r/sig)**2)) 

现在%f1(res1,res2) 在我的机器上以 60 微秒的速度运行。

Answer 1

作为参考，f1 你定义的方式给了我：

10000 次循环，3 次中的最佳：每个循环 20.9 µs

避免每次都计算np.log(2*np.pi)的简单优化：

log2pi = np.log(2*np.pi)
f2 = lambda r, sig: -0.5*(np.sum(2*np.log(sig) + log2pi + (r/sig)**2))

这将我的timeit 输出减少了大约 15%：

100000 次循环，3 次取胜：每个循环 17.8 µs

另外，计算x*x 通常比计算x**2 更快。使用这种方法，您也可以避免2* 操作。

def f3(r, sig):
   sig2 = sig * sig
   return -0.5*(np.sum(np.log(sig2) + log2pi + (r*r/sig2)))

100000 次循环，3 次中的最佳：每个循环 15 µs

此外，根据您使用参数res, np.sqrt(res) 调用此函数的频率，您绝对可以考虑将其简化为：

f4 = lambda a: -0.5*np.sum(np.log(a) + log2pi + a)

您可以检查f4(res) == f1(res, np.sqrt(res))，它的速度也明显加快（大约 44%）。

100000 次循环，3 次中的最佳：每个循环 11.7 µs

Answer 2

我很确定最昂贵的东西就是 logs 和那些你可以保存的东西（即调用一次而不是 sig.size - 在这个例子中 = 1000 - 次）。基础身份是log(a) + log(b) = log(a*b)。让我的速度提高了三倍。

顺便说一句。我将产生空数组的random.rand(0, 1, 1000) 更改为random.rand(1000)。另外我认为sig=sqrt(res) 不是一般假设，所以我没有尝试利用它。

我做了一些其他优化，但它们的效果远不及摆脱 logs 的影响。 1) 预计算 log(2*pi) 2) 使用 np.dot 进行最后一项。

更新：

正如@JPdL 指出的那样，这种方法容易受到上溢和下溢的影响。可以使用分块来交换各种性能以防止出现这种情况，请参阅下面的f4. f5. f6。

import numpy as np
np.random.seed(1)

res = np.random.rand(1000)
f1 = lambda r, sig: -0.5*(np.sum(2*np.log(sig) + np.log(2*np.pi) + (r/sig)**2)) 

def f2(r, sig):
    return -0.5*(np.sum(2*np.log(sig) + np.log(2*np.pi) + (r/sig)**2))

# gung ho
def f3(r, sig, precomp=np.log(2*np.pi)):
    rs = r/sig
    return -np.log(np.prod(sig)) - 0.5*(r.size * precomp + np.dot(rs, rs))

# chunk and hope for the best
def f4(r, sig, chnk=32, precomp=np.log(2*np.pi)):
    rs = r/sig
    sumlogsig = np.log(np.multiply.reduceat(sig, np.arange(0, len(r), chnk))).sum()
    return -sumlogsig - 0.5*(r.size * precomp + np.dot(rs, rs))

# chunk and check for extreme values
def f5(r, sig, chnk=32,
       precomp=np.log(2*np.pi), precomp2=np.exp([-8, 8])):
    rs = r/sig
    bad = np.where((sig<precomp2[0]) | (sig>precomp2[1]))[0]
    sumlogsig = np.log(sig[bad]).sum()
    keep = sig[bad]
    sig[bad] = 1
    sumlogsig += np.log(np.multiply.reduceat(sig, np.arange(0, len(r), chnk))).sum()
    sig[bad] = keep
    return -sumlogsig - 0.5*(r.size * precomp + np.dot(rs, rs))

# chunk and try to be extra clever
def f6(r, sig, chnk=512,
       precomp=np.log(2*np.pi), precomp2=np.exp(np.arange(-18, 19))):
    binidx = np.digitize(sig, precomp2[1::2])<<1
    rs = r/sig
    sig = sig*precomp2[36 - binidx]
    bad = np.where((binidx==0) | (binidx==36))[0]
    sumlogsig = binidx.sum() - 18*r.size + np.log(sig[bad]).sum()
    sig[bad] = 1
    sumlogsig += np.log(np.multiply.reduceat(sig, np.arange(0, len(r), chnk))).sum()
    return -sumlogsig - 0.5*(r.size * precomp + np.dot(rs, rs))

from timeit import timeit
sr = np.sqrt(res)
print(timeit('f1(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f1':f1}))
print(timeit('f2(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f2':f2}))
print(timeit('f3(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f3':f3}))
print(timeit('f4(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f4':f4}))
print(timeit('f5(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f5':f5}))
print(timeit('f6(res,sr)', number=100, globals={'res':res, 'np':np, 'sr':sr, 'f6':f6}))
print(f1(res,np.sqrt(res)))
print(f2(res,np.sqrt(res)))
print(f3(res,np.sqrt(res)))
print(f4(res,np.sqrt(res)))
print(f5(res,np.sqrt(res)))
print(f6(res,np.sqrt(res)))

样本输出：

0.004246247990522534
0.00418912700843066
0.0011273059935774654
0.0013386670034378767
0.0022679700050503016
0.004274581006029621
-662.250886322
-662.250886322
-662.250886322
-662.250886322
-662.250886322
-662.250886322

Answer 3

实际上，您可以跳过np.sum()，因为您在通过应用+ 运算符获得的单个参数上调用它。此外，np.log(2*np.pi) 是常量，不依赖于 lambda 参数，因此您可以在 lambda 之外计算并重用。此外，这可能不会有太大变化，但请检查 a**2 是比 a*a 快还是慢；如果后者更快，计算一次d = r/sig 并使用d*d 而不是(r/sig)**2。您还可以检查0.5*x 是比x/2.0 快还是慢。这并不算多，但它可能会让您加速一到两个使用秒。

现在，一些数学。 根据上述建议，您会得到：

log_of_pi = np.log(2*np.pi)
f1 = lambda r, sig: -0.5*(2*np.log(sig) + log_of_pi + (r/sig)**2))
f1(res,np.sqrt(res)) 

所以基本上，你计算：

c1*(2*log(sqrt(r)) + c2 + (r/sqrt(r))**2)

首先：

log(sqrt(r)) = log(r**(1/2)) = 1/2*log(r)

还有：

(r/sqrt(r))**2 = sqrt(r)**2 = abs(r)

所以基本上，您正在尝试计算：

c1*(2*log(sqrt(r)) + c2 + (r/sqrt(r))**2) = c1*(2*1/2*log(r) + c2 + abs(r)) = c1*(log(r) + c2 + abs(r))

因此，您可以跳过昂贵的除法平方（总和的最后一部分）和计算 r 的平方根（在调用 f1 之前完成）。

当然，前提是您将始终使用 sqrt(r) 作为 sig。

PS。尊重 grovina，他 answered 在我之前 15 秒有同样的想法，并且也提供了测量结果。你愿意接受我的回答吗，接受另一个，因为它更完整。

Answer 4

您正在使用 lambda 进行隐式复制。您可以通过在可能的情况下使用函数的就地变体来提高性能。 http://ipython-books.github.io/featured-01/

def f2(r, sig):
   x = np.log(sig)
   x *= 2  # in-place operation

   x += np.log(2*np.pi)  # in-place operation

   y = r / sig
   y **= 2  # in-place operation

   x += y  # in-place operation

   return -0.5 * x.sum()

另外，检查

的输出

import numpy.distutils.system_info as sysinfo
sysinfo.show_all()

...

openblas_lapack_info:
  libraries openblas not found in ['C:\\Users\\User\\AppData\\Local\\Programs\\Python\\Python35-32\\lib', 'C:\\', 'C:\\Users\\User\\AppData\\Local\\Programs\\Python\\Python35-32\\libs']
  NOT AVAILABLE

如果 blas、lapack 和 atlas 是 NOT AVAILABLE，那么您无法充分利用 numpy。不是由一个长镜头。这些库使 numpy 能够使用SIMD 操作（单指令多数据——一种并行化形式）。