numpy矩阵乘法到三角形/稀疏存储？

Question

我正在处理一个非常大的稀疏矩阵乘法 (matmul) 问题。举个例子吧：

• A 是一个二进制 (75 x 200,000) 矩阵。它很稀疏，所以我使用 csc 进行存储。我需要做以下matmul操作：

• B = A.transpose() * A

• 输出将是一个大小为 200Kx200K 的稀疏对称矩阵。

不幸的是，B 将 变大以存储在我笔记本电脑的 RAM（或“内核中”）中。另一方面，我很幸运，因为 B 有一些属性可以解决这个问题。

由于 B 将沿对角线对称且稀疏，我可以使用三角矩阵（上/下）来存储 matmul 运算的结果，而稀疏矩阵存储格式可以进一步减小大小。

我的问题是...能否提前告知 numpy 或 scipy 输出存储要求是什么样的，以便我可以使用 numpy 选择存储解决方案并避免“矩阵太大”运行时计算几分钟（小时）后出错？

换句话说，矩阵乘法的存储需求是否可以通过使用近似计数算法分析两个输入矩阵的内容来近似？

• https://en.wikipedia.org/wiki/Approximate_counting_algorithm

如果没有，我正在寻找蛮力解决方案。来自以下 Web 链接的涉及 map/reduce、核外存储或 matmul 细分解决方案（strassens 算法）的内容：

一对 Map/Reduce 问题细分解决方案

• http://www.norstad.org/matrix-multiply/index.html
• http://bpgergo.blogspot.com/2011/08/matrix-multiplication-in-python.html

核外 (PyTables) 存储解决方案

• Very large matrices using Python and NumPy

matmul 细分解：

• https://en.wikipedia.org/wiki/Strassen_algorithm
• http://facultyfp.salisbury.edu/taanastasio/COSC490/Fall03/Lectures/FoxMM/example.pdf
• http://eli.thegreenplace.net/2012/01/16/python-parallelizing-cpu-bound-tasks-with-multiprocessing/

提前感谢您的任何建议、cmets 或指导！

Answer 1

由于您需要矩阵与其转置的乘积，因此[m, n] 处的值基本上将是原始矩阵中m 和n 列的点积。

我将使用以下矩阵作为玩具示例

a = np.array([[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
              [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
              [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1]])
>>> np.dot(a.T, a)
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2]])

它的形状为(3, 12)，有7 个非零条目。它与它的转置结果当然是(12, 12) 的形状，并且有 16 个非零条目，其中 6 个在对角线上，因此它只需要存储 11 个元素。

您可以通过以下两种方式之一很好地了解输出矩阵的大小：

企业社会责任格式

如果您的原始矩阵有C 非零列，那么您的新矩阵将最多有C**2 非零条目，其中C 在对角线上，并且保证不为零，并且其余条目中您只需要保留一半，因此最多为(C**2 + C) / 2 非零元素。当然，其中许多也将是零，所以这可能是一个严重的高估。

如果您的矩阵以csr 格式存储，则相应scipy 对象的indices 属性具有一个包含所有非零元素的列索引的数组，因此您可以轻松地将上述估计计算为：

>>> a_csr = scipy.sparse.csr_matrix(a)
>>> a_csr.indices
array([ 2, 11,  1,  7, 10,  4, 11])
>>> np.unique(a_csr.indices).shape[0]
6

因此，有 6 列具有非零条目，因此估计最多有 36 个非零条目，比实际的 16 个多得多。

CSC 格式

如果我们有行索引而不是非零元素的列索引，我们实际上可以进行更好的估计。为了使两列的点积不为零，它们必须在同一行中有一个非零元素。如果给定行中有R 非零元素，它们将为产品贡献R**2 非零元素。当你对所有行求和时，你肯定会多次计算一些元素，所以这也是一个上限。

矩阵的非零元素的行索引位于稀疏 csc 矩阵的 indices 属性中，因此该估计值可以按如下方式计算：

>>> a_csc = scipy.sparse.csc_matrix(a)
>>> a_csc.indices
array([1, 0, 2, 1, 1, 0, 2])
>>> rows, where = np.unique(a_csc.indices, return_inverse=True)
>>> where = np.bincount(where)
>>> rows
array([0, 1, 2])
>>> where
array([2, 3, 2])
>>> np.sum(where**2)
17

这与真正的 16 非常接近！而且这个估计实际上与以下内容相同并非巧合：

>>> np.sum(np.dot(a.T,a),axis=None)
17

无论如何，下面的代码应该可以让你看到估计是相当不错的：

def estimate(a) :
    a_csc = scipy.sparse.csc_matrix(a)
    _, where = np.unique(a_csc.indices, return_inverse=True)
    where = np.bincount(where)
    return np.sum(where**2)

def test(shape=(10,1000), count=100) :
    a = np.zeros(np.prod(shape), dtype=int)
    a[np.random.randint(np.prod(shape), size=count)] = 1
    print 'a non-zero = {0}'.format(np.sum(a))
    a = a.reshape(shape)
    print 'a.T * a non-zero = {0}'.format(np.flatnonzero(np.dot(a.T,
                                                                a)).shape[0])
    print 'csc estimate = {0}'.format(estimate(a))

>>> test(count=100)
a non-zero = 100
a.T * a non-zero = 1065
csc estimate = 1072
>>> test(count=200)
a non-zero = 199
a.T * a non-zero = 4056
csc estimate = 4079
>>> test(count=50)
a non-zero = 50
a.T * a non-zero = 293
csc estimate = 294