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

如果不是这样，我寻找到一个蛮力解决方案。一些涉及到的map/reduce，外的核心存储，或从下面的网页链接MATMUL细分的解决方案（施特拉森算法）：

一对夫妇的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/

预先感谢任何建议，意见或指导！

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 FORMAT

如果你原来的矩阵具有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**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