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

Question

我正在处理一个非常大的稀疏矩阵乘法（matmul）问题。作为一个例子，让我们说：

>

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

B=A.转置（）*A

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

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

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

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

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

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

如果没有，我正在寻找一种蛮力解决方案。从以下网络链接中获取涉及map/duce、核心外存储或matmul细分解决方案（strassens算法）的内容：

一对映射/约简问题细分解决方案

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

核心外（PyTables）存储解决方案

• 使用Python和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/

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

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个元素。

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

CSR格式

如果原始矩阵有非零列，则新矩阵最多将有非零项，其中非零项位于对角线中，并确保不为零，其余项只需保留一半，即最多有非零元素。当然，其中许多也将为零，因此这可能是一个严重的高估。

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

>>> 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矩阵的索引属性中，因此该估计值可按如下方式计算：

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