稀疏矩阵上的唯一性

Question

看起来numpy的einsum函数不适用于scipy.sparse矩阵。 einsum可以用稀疏矩阵进行各种各样的事情吗？

回复@ eickenberg的回答：我想要的特定义务是numpy.einsum("ki,kj->ij",A,A) - 这些行的外部产品的总和。

Answer 1

scipy.sparse矩阵的一个限制是它们代表线性算子，因此保持二维，这导致了一个问题：您在寻找哪种操作？

在一对2D矩阵的所有einsum操作是非常容易的，而不使用dot，transpose和逐点操作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将索引字符串转换为使用C版本np.nditer的计算。 http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html是一个很好的介绍nditer。最后请注意Cython示例。

https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py是einsum的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]])