numpy - 稀疏矩阵上的 einsum

Question

似乎 numpy 的einsum函数不适用于scipy.sparse矩阵。einsum有没有其他方法可以做稀疏矩阵可以做的事情？

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

Answer 1

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.

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]])

Answer 2

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

一对 2D 矩阵上的所有einsum操作都非常容易编写，无需使用einsum和逐点操作，前提是结果不超过二维。dottranspose

因此，如果您需要对多个稀疏矩阵进行特定操作，则很可能无需einsum.

更新：一个具体的实现方式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，您将获得您的特定结果。