枕中稀疏矩阵上的Log-sum-exp技巧

Question

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

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

更新

更好的做法是，我要找的是执行日志和-exp技巧的东西，简单地按exp‘m顺序执行，在scipy.sparse中对行进行求和并执行日志处理是非常简单的。较少琐碎的是以一种干净的方式计算沿行的最大值并减去它，因为稀疏矩阵行中的每个元素被减去相应的最大向量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。

UPDATE Ok，我误解了这个问题: 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并应用上述方法.