首页 文章

Numpy:智能矩阵乘法到稀疏结果矩阵

提问于
浏览
4

在numpy的python中,假设我有两个矩阵:

  • S ,稀疏的 x*x 矩阵

  • M ,一个密集的 x*y 矩阵

现在我想做 np.dot(M, M.T) ,它将返回一个密集的 x*x 矩阵 S_ .

但是,我只关心 S 中非零的单元格,这意味着如果我这样做,它对我的应用程序没有任何影响

S_ = S*S_

显然,这将浪费操作,因为我想省去 S 中给出的不相关的单元格 . 请记住,在矩阵乘法中

S_[i,j] = np.sum(M[i,:]*M[:,j])

所以我想只为 i,j 做这个操作 S[i,j]=True .

是否通过在C中运行的numpy实现以某种方式支持它,以便我不需要使用python循环实现它?


EDIT 1 [solved]: 我仍然有这个问题,实际上 M 现在也很稀疏 .

现在,给定 S 的行和列,我实现了这样:

data = np.array([ M[rows[i],:].dot(M[cols[i],:]).data[0] for i in xrange(len(rows)) ])
S_   = csr( (data, (rows,cols)) )

......但它仍然很慢 . 有什么新想法吗?

EDIT 2: jdehesa给出了一个很好的解决方案,但我想节省更多的内存 .

解决方案是执行以下操作:

data = M[rows,:].multiply(M[cols,:]).sum(axis=1)

然后从 rowscolsdata 构建一个新的稀疏矩阵 .

但是,当运行上面的行时,scipy构建一个(连续的)numpy数组,其元素与第一个子矩阵的 nnz 和第二个子矩阵的 nnz 一样多,这可能导致 MemoryError 在我的情况下 .

为了节省更多内存,我想迭代地将每一行与其各自的“伙伴”列相乘,然后求和并丢弃结果向量 . 使用简单的python实现这一点,基本上我回到了非常慢的版本 .

有解决这个问题的快捷方法吗?

1 回答

  • 1

    以下是使用NumPy / SciPy进行密集和稀疏 M 矩阵的方法:

    import numpy as np
    import scipy.sparse as sp
    
    # Coordinates where S is True
    S = np.array([[0, 1],
                  [3, 6],
                  [3, 4],
                  [9, 1],
                  [4, 7]])
    
    # Dense M matrix
    # Random big matrix
    M = np.random.random(size=(1000, 2000))
    # Take relevant rows and compute values
    values = np.sum(M[S[:, 0]] * M[S[:, 1]], axis=1)
    # Make result matrix from values
    result = np.zeros((len(M), len(M)), dtype=values.dtype)
    result[S[:, 0], S[:, 1]] = values
    
    # Sparse M matrix
    # Construct sparse M as COO matrix or any other way
    M = sp.coo_matrix(([10, 20, 30, 40, 50],  # Data
                       ([0, 1, 3, 4, 6],      # Rows
                        [4, 4, 5, 5, 8])),    # Columns
                      shape=(1000, 2000))
    # Convert to CSR for fast row slicing
    M_csr = M.tocsr()
    # Take relevant rows and compute values
    values = M_csr[S[:, 0]].multiply(M_csr[S[:, 1]]).sum(axis=1)
    values = np.squeeze(np.asarray(values))
    # Construct COO sparse matrix from values
    result = sp.coo_matrix((values, (S[:, 0], S[:, 1])), shape=(M.shape[0], M.shape[0]))
    

相关问题