python - scipy 中的向量化逐项 logsumexp

Question

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

我可以用矢量化的方式而不是使用双 for 循环来做到这一点吗？

Answer 1

要模仿 logsumexp 的行为，您只需在获取 exp 之前减去其参数的 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]])