在 numpy 中乘以对数概率矩阵的数值稳定方法

我需要获取包含对数概率的两个 NumPy 矩阵（或其他二维数组）的矩阵乘积。天真的方式np.log(np.dot(np.exp(a), np.exp(b)))由于明显的原因而不是首选。

Using

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) 

可以工作，但运行速度比np.log(np.dot(np.exp(a), np.exp(b)))

Using

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 扩展来计算这个，但我当然宁愿避免这种情况。是否有一种既定的方法可以做到这一点，或者有人知道执行此计算的内存密集程度较低的方法吗？

EDIT:感谢 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 + max_b
    return c

将此方法与上面发布的方法进行快速比较（logdot_old）使用 iPython 的魔力%timeit函数产生以下结果：

In  [1] a = np.log(np.random.rand(1000,2000))

In  [2] b = np.log(np.random.rand(2000,1500))

In  [3] x = logdot(a, b)

In  [4] y = logdot_old(a, b) # this takes a while

In  [5] np.any(np.abs(x-y) > 1e-14)
Out [5] False

In  [6] %timeit logdot_old(a, b)
1 loops, best of 3: 1min 18s per loop

In  [6] %timeit logdot(a, b)
1 loops, best of 3: 264 ms per loop

显然拉斯曼的方法抹杀了我的方法！

logsumexp通过计算等式右侧的值来工作

log(∑ exp[a]) = max(a) + log(∑ exp[a - max(a)])

即，它在开始求和之前取出最大值，以防止溢出exp。在进行向量点积之前也可以应用同样的方法：

log(exp[a] ⋅ exp[b])
 = log(∑ exp[a] × exp[b])
 = log(∑ exp[a + b])
 = max(a + b) + log(∑ exp[a + b - max(a + b)])     { this is logsumexp(a + b) }

但通过在推导中采取不同的方式，我们得到

log(∑ exp[a] × exp[b])
 = max(a) + max(b) + log(∑ exp[a - max(a)] × exp[b - max(b)])
 = max(a) + max(b) + log(exp[a - max(a)] ⋅ exp[b - max(b)])

最终形式的内部有一个向量点积。它也很容易扩展到矩阵乘法，所以我们得到了算法

def logdotexp(A, B):
    max_A = np.max(A)
    max_B = np.max(B)
    C = np.dot(np.exp(A - max_A), np.exp(B - max_B))
    np.log(C, out=C)
    C += max_A + max_B
    return C

这创建了两个A大小的临时体和两个B- 大小的，但每一个都可以通过以下方式消除

exp_A = A - max_A
np.exp(exp_A, out=exp_A)

类似地对于B。 （如果函数可以修改输入矩阵，则可以消除所有临时矩阵。）