我想在numpy中实现以下问题，这是我的代码。
我已经尝试了针对此问题的一个包含for循环的numpy代码。我想知道是否有更高效的方法来进行这个计算？我非常感激您的帮助！
k, d = X.shapem = Y.shape[0]c1 = 2.0*sigma**2c2 = 0.5*np.log(np.pi*c1)c3 = np.log(1.0/k)L_B = np.zeros((m,))for i in xrange(m): if i % 100 == 0: print i L_B[i] = np.log(np.sum(np.exp(np.sum(-np.divide( np.power(X-Y[i,:],2), c1)-c2,1)+c3)))print np.mean(L_B)我想到了通过创建一个3D张量np.expand_dims(X, 2).repeat(Y.shape[0], 2)-Y，这样后续的计算可以通过广播来完成，但当m很大时，这会浪费大量内存。
我还认为np.einsum()只是利用了for循环，所以可能不是那么高效，如果我错了请纠正我。
有什么想法吗？
回答：
优化阶段#1
我的第一级优化是将包含循环的代码直接翻译成基于广播的代码，通过引入一个新轴，但这样做在内存使用上不是那么高效，如下所示 –
p1 = (-((X[:,None] - Y)**2)/c1)-c2p11 = p1.sum(2)p2 = np.exp(p11+c3)out = np.log(p2.sum(0)).mean()优化阶段#2
考虑到我们打算将常数的操作分离出来，我进行了一些优化，最终得到以下结果 –
c10 = -c1c20 = X.shape[1]*c2subs = (X[:,None] - Y)**2p00 = subs.sum(2)p10 = p00/c10p11 = p10-c20p2 = np.exp(p11+c3)out = np.log(p2.sum(0)).mean()优化阶段#3
进一步优化，寻找可以优化操作的地方，我最终使用了Scipy的cdist来替代方差和sum-reduction的繁重工作。这应该非常节省内存，并给出了最终的实现，如下所示 –
from scipy.spatial.distance import cdist# 设置常数c10 = -c1c20 = X.shape[1]*c2c30 = c20-c3c40 = np.exp(c30)c50 = np.log(c40)# 获取与循环操作对应的阶段性操作p1 = cdist(X, Y, 'sqeuclidean')p2 = np.exp(p1/c10).sum(0)out = np.log(p2).mean() - c50运行时间测试
方法 –
def loopy_app(X, Y, sigma): k, d = X.shape m = Y.shape[0] c1 = 2.0*sigma**2 c2 = 0.5*np.log(np.pi*c1) c3 = np.log(1.0/k) L_B = np.zeros((m,)) for i in xrange(m): L_B[i] = np.log(np.sum(np.exp(np.sum(-np.divide( np.power(X-Y[i,:],2), c1)-c2,1)+c3))) return np.mean(L_B)def vectorized_app(X, Y, sigma): # 设置常数 k, d = D_A.shape c1 = 2.0*sigma**2 c2 = 0.5*np.log(np.pi*c1) c3 = np.log(1.0/k) c10 = -c1 c20 = X.shape[1]*c2 c30 = c20-c3 c40 = np.exp(c30) c50 = np.log(c40) # 获取与循环操作对应的阶段性操作 p1 = cdist(X, Y, 'sqeuclidean') p2 = np.exp(p1/c10).sum(0) out = np.log(p2).mean() - c50 return out时间和验证 –
In [294]: # 设置输入，其中m(=D_B.shape[0])是一个大数 ...: X = np.random.randint(0,9,(100,10)) ...: Y = np.random.randint(0,9,(10000,10)) ...: sigma = 2.34 ...: In [295]: np.allclose(loopy_app(X, Y, sigma),vectorized_app(X, Y, sigma))Out[295]: TrueIn [296]: %timeit loopy_app(X, Y, sigma)1 loops, best of 3: 225 ms per loopIn [297]: %timeit vectorized_app(X, Y, sigma)10 loops, best of 3: 23.6 ms per loopIn [298]: # 设置输入，其中m(=Y.shape[0])是一个更大的数 ...: X = np.random.randint(0,9,(100,10)) ...: Y = np.random.randint(0,9,(100000,10)) ...: sigma = 2.34 ...: In [299]: np.allclose(loopy_app(X, Y, sigma),vectorized_app(X, Y, sigma))Out[299]: TrueIn [300]: %timeit loopy_app(X, Y, sigma)1 loops, best of 3: 2.27 s per loopIn [301]: %timeit vectorized_app(X, Y, sigma)1 loops, best of 3: 243 ms per loop大约有10倍的加速！
