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)存储解决方案

• 使用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

由于您追求的是矩阵与其转置的乘积，因此 at 的值[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格式存储，则indices相应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乘积贡献非零元素。当您对所有行求和时，您必然会对某些元素进行多次计数，因此这也是一个上限。

矩阵非零元素的行索引属于indices稀疏 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