python - scipy 中稀疏矩阵的 Log-sum-exp 技巧

Question

将 scipy.misc.logsumexp 之类的东西应用于稀疏矩阵(例如 scipy.sparse.csr_matrix)并指定一个轴的最佳方法是什么？

重点是将零排除在计算之外。

更新

最好指明我正在寻找执行 log-sum-exp trick 的东西，在 scipy.sparse 中简单地连续执行 exp elem-wise，对行求和并按 elem-wise 执行日志是微不足道的。更简单的是以一种干净的方式计算沿行的最大值并减去它，因为稀疏矩阵行中的每个元素都减去相应的最大向量 elem(最后保留一个稀疏矩阵)。

Answer 1

CSR 矩阵 X 的非零项通过以下方式获得

X[i].data

和(排列)实际行的值将通过向其附加 X.shape[1] - len(X[i].data) 零来获得。

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

对于向量a。让我们设置 b = X[i].data 和 k = X.shape[1] - len(X[i].data) 并表示我们之前排列的行X 为

(b, 0ₖ)

使用 0ₖ 表示长度为 k 的零向量，并使用 (⋅, ⋅) 进行连接。然后

logsumexp((b, 0ₖ))
 = max((b, 0ₖ)) + log(∑ exp[(b, 0ₖ) - max((b, 0ₖ))])
 = max(max(b), 0) + log(∑ exp[(b, 0ₖ) - max(max(b), 0)])
 = max(max(b), 0) + log(∑ exp[b - max(max(b), 0)] + ∑ exp[0ₖ - max(max(b), 0)])
 = max(max(b), 0) + log(∑ exp[b - max(max(b), 0)] + k × exp[-max(max(b), 0)])

这样我们就得到了算法

def logsumexp_csr_row(x):
    data = x.data
    mx = max(np.max(data), 0)
    tmp = data - mx
    r = np.exp(tmp, out=tmp).sum()
    k = X.shape[1] - len(data)
    return mx + np.log(r + k * np.exp(-mx))

对于 CSR 行向量。将此算法扩展到完整矩阵可以通过列表理解轻松完成，尽管更有效的形式是使用 indptr 遍历行:

def logsumexp_csr_rows(X):
    result = np.empty(X.shape[0])
    for i in range(X.shape[0]):
        data = X.data[X.indptr[i]:X.indptr[i+1]]
        # fill in from logsumexp_csr_row
        result[i] = mx + np.log(r + k * np.exp(-mx))
    return result

列式版本要复杂得多；转置矩阵并转换回 CSR 可能是最简单的。

---

更新 好的，我误解了这个问题:OP 根本没有兴趣处理零，所以上面的推导是无用的，算法应该是

def logsumexp_row_nonzeros(X):
    result = np.empty(X.shape[0])
    for i in range(X.shape[0]):
        result[i] = logsumexp(X.data[X.indptr[i]:X.indptr[i+1]])
    return result

这只是填写 CSR 矩阵上按行操作的一般方案。对于列式，转置，转换回 CSR 并应用上述内容。