稀疏矩阵上的einsum

Question

numpy的einsum函数似乎不适用于scipy.sparse矩阵。有没有einsum可以用稀疏矩阵做的那些事情的替代方案？

针对@eikenberg的回答：我想要的特定einsum是numpy.einsum("ki,kj->ij",A,A)——行的外积的和。

Answer 1

scipy.sparse矩阵的一个限制是，它们表示线性算子，因此保持二维，这就引出了一个问题：你想做哪个运算？

在没有使用dot、transpose和逐点运算的einsum的情况下，只要结果不超过二维，对一对2D矩阵的所有einsum运算都非常容易写入。

因此，如果你需要对一些稀疏矩阵进行特定的运算，那么很可能你可以在没有einsum的情况下编写它。

UPDATE：实现np.einsum("ki, kj -> ij", A, A)的具体方法是A.T.dot(A)。为了说服自己，请尝试以下示例：

import numpy as np
rng = np.random.RandomState(42)
a = rng.randn(3, 3)
b = rng.randn(3, 3)
the_einsum_ab = np.einsum("ki, kj -> ij", a, b)
the_a_transpose_times_b = a.T.dot(b)
# We write a test in order to assert equality
from numpy.testing import assert_array_equal
assert_array_equal(the_einsum_ab, the_a_transpose_times_b)  # This passes, so equality

这个结果稍微更一般一些。现在，如果您使用b = a，您将获得特定的结果。

Answer 2

einsum使用np.nditer的C版本将索引字符串转换为计算。http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html是对CCD_ 16的一个很好的介绍。请特别注意末尾的Cython示例。

https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py是CCD_ 18的Python模拟。

scipy.sparse有自己的代码（最终用C）来执行基本运算、求和、矩阵乘法等。稀疏矩阵有自己的数据结构。它们可以是列表、字典或一组numpy数组。可以使用Numpy表示法，因为sparse具有适当的__xxx__方法。

稀疏矩阵是一个matrix，一个2d数组对象。可以编写稀疏的einsum，但它最终将使用稀疏矩阵乘法，而不是nditer。所以，这充其量只是一种符号上的便利。

稀疏csr_matrix.dot为：

def dot(self, other):
    """Ordinary dot product
    ...
    """
    return self * other
A=sparse.csr_matrix([[1,2],[3,4]])
A.dot(A.T).A
(A*A.T).A
A.__rmul__(A.T).A
A.__mul__(A.T).A
np.einsum('ij,kj',A.A,A.A)
# array([[ 5, 11],
#        [11, 25]])