scipy中的向量化逐项logsumexp答案

【问题标题】：Vectorized entry-wise logsumexp in scipyscipy中的向量化逐项logsumexp
【发布时间】：2018-07-01 22:00:57
【问题描述】：

给定一个二维数组 numpy 数组 A 和一个一维数组 c，我想计算二维数组 B 和条目 B[i, j] = scipy.special.logsumexp(np.append(c, A[i, j])).

我可以用矢量化的方式而不是使用双 for 循环吗？

【问题讨论】：

【解决方案1】：

你为什么不直接使用np.log(np.sum(np.exp(c)) + np.exp(A))？

[根据 Paul Panzer 的评论更新]

_c = np.broadcast_to(c, (*A.shape, *c.shape))
B = scipy.special.logsumexp(np.append(_c, A[...,np.newaxis], axis=-1), axis=-1)

【讨论】：

logsumexp 的全部意义在于解决exp 中的溢出问题。
我明白了。非常感谢。
元组拆包重新打包有点可爱。

【解决方案2】：

要模仿logsumexp 的行为，您只需在取exp, sum, log 之前减去其参数的max（检查source code），然后在最后重新添加它。因此，您可以执行以下操作：

>>> import numpy as np
>>> from scipy import special
>>> 
>>> A = np.random.uniform(900, 1100, (4, 4))
>>> c = np.random.uniform(950, 1050, (7,))
>>> 
>>> cm = np.max(c)
>>> mask = A > cm
>>> B = np.empty_like(A)
>>> B[mask] = A[mask] + np.log(np.exp(np.subtract.outer(cm, A[mask])).sum(axis=-1) + 1)
>>> B[~mask] = cm + np.log(np.exp(c - cm).sum() + np.exp(A[~mask] - cm))
>>> 
# compute via logsumexp for reference
>>> cA = np.empty((8, 4, 4))
>>> cA[:-1] = c[:, None, None]
>>> cA[-1] = A
>>> special.logsumexp(cA, axis=0)
array([[ 1048.88855012,  1048.88854955,  1069.83524808,  1048.88854955],
       [ 1048.88854955,  1048.88854955,  1048.88877212,  1048.93142975],
       [ 1048.88854955,  1067.59166572,  1048.88854955,  1069.78737913],
       [ 1048.88854955,  1048.88854955,  1098.61910373,  1072.76058998]])
>>> B
array([[ 1048.88855012,  1048.88854955,  1069.83524809,  1048.88854955],
       [ 1048.88854955,  1048.88854955,  1048.88877212,  1048.93142975],
       [ 1048.88854955,  1067.59166572,  1048.88854955,  1069.78737914],
       [ 1048.88854955,  1048.88854955,  1098.61910374,  1072.76058999]])

【讨论】：

感谢您的回答！实际上，我想知道是否可以将A 广播到形状为(A.shape[0], A.shape[1], len(c)) 的3D 数组Ac（在我的情况下len(c) 非常小），然后简单地应用logsumexp(Ac, axis=-1)？虽然我无法弄清楚如何实现这种广播......
@希尔伯特。当然Ac = np.empty(A.shape + (len(c)+1,), dtype=A.dtype)Ac[..., 0] = AAc[..., 1:] = c。没有测试它，但认为它应该可以工作。