小编mar*_*art的帖子

我需要获取包含对数概率的两个NumPy矩阵(或其他2d数组)的矩阵乘积.np.log(np.dot(np.exp(a), np.exp(b)))出于显而易见的原因,天真的方式不是优选的.

运用

from scipy.misc import logsumexp
res = np.zeros((a.shape[0], b.shape[1]))
for n in range(b.shape[1]):
    # broadcast b[:,n] over rows of a, sum columns
    res[:, n] = logsumexp(a + b[:, n].T, axis=1) 

工作但运行速度比慢100倍 np.log(np.dot(np.exp(a), np.exp(b)))

运用

logsumexp((tile(a, (b.shape[1],1)) + repeat(b.T, a.shape[0], axis=0)).reshape(b.shape[1],a.shape[0],a.shape[1]), 2).T

或者其他瓦片和重塑的组合也起作用,但是比上面的循环运行得更慢,因为实际大小的输入矩阵需要非常大量的存​​储器.

我目前正在考虑在C中编写一个NumPy扩展来计算它,但当然我宁愿避免这种情况.是否有既定的方法来执行此操作,或者是否有人知道执行此计算的内存密集程度较低的方法？

编辑: 感谢larsmans提供此解决方案(参见下面的推导):

def logdot(a, b):
    max_a, max_b = np.max(a), np.max(b)
    exp_a, exp_b = a - max_a, b - max_b
    np.exp(exp_a, out=exp_a)
    np.exp(exp_b, out=exp_b)
    c = np.dot(exp_a, exp_b)
    np.log(c, out=c)
    c += max_a + …