首页 文章

有效地乘以Numpy / Scipy稀疏矩阵和密集矩阵

提问于
浏览 次
14

我正在努力实现以下等式：

X =(Y.T * Y + Y.T * C * Y) ^ -1

Y是（n×f）矩阵，C是（n×n）对角线; n约为300k，f将在100到200之间变化 . 作为优化过程的一部分，该等式将使用近1亿次，因此必须非常快速地处理 .

Y是随机初始化的，C是一个非常稀疏的矩阵，对角线上300k中只有少数数字与0不同 . 由于Numpy的对角线函数创建了密集矩阵，我创建了C作为稀疏csr矩阵 . 但是当试图解决方程式的第一部分时：

r = dot(C, Y)

由于内存限制，计算机崩溃 . 我决定然后尝试将Y转换为csr_matrix并进行相同的操作：

r = dot(C, Ysparse)

这种方法需要 1.38 ms . 但是这个解决方案有点"tricky"因为我使用稀疏矩阵来存储密集的矩阵，我想知道它真的有多高效 .

所以我的问题是，是否有一些方法可以将稀疏的C和密集的Y相乘，而不必将Y变为稀疏并提高性能？如果不知何故C可以表示为对角线密集而不消耗大量内存，这可能会导致非常有效的性能，但我不知道这是否可行 .

我感谢您的帮助！

python performance numpy scipy sparse-matrix

3 回答

  • 7

    当计算r = dot（C，Y）时，点积遇到内存问题的原因是因为numpy的点函数不具有处理稀疏矩阵的本机支持 . 正在发生的是numpy认为稀疏矩阵C是一个python对象，而不是一个numpy数组 . 如果您进行小规模检查，您可以直接看到问题：

    >>> from numpy import dot, array
>>> from scipy import sparse
>>> Y = array([[1,2],[3,4]])
>>> C = sparse.csr_matrix(array([[1,0], [0,2]]))
>>> dot(C,Y)
array([[  (0, 0)    1
  (1, 1)    2,   (0, 0) 2
  (1, 1)    4],
  [  (0, 0) 3
  (1, 1)    6,   (0, 0) 4
  (1, 1)    8]], dtype=object)

    显然，上述内容并不是您感兴趣的结果 . 相反，您要做的是使用scipy的sparse.csr_matrix.dot函数进行计算：

    r = sparse.csr_matrix.dot(C, Y)

    或者更紧凑

    r = C.dot(Y)
    回复于
  • 2

    尝试：

    import numpy as np
from scipy import sparse

f = 100
n = 300000

Y = np.random.rand(n, f)
Cdiag = np.random.rand(n) # diagonal of C
Cdiag[np.random.rand(n) < 0.99] = 0

# Compute Y.T * C * Y, skipping zero elements
mask = np.flatnonzero(Cdiag)
Cskip = Cdiag[mask]

def ytcy_fast(Y):
    Yskip = Y[mask,:]
    CY = Cskip[:,None] * Yskip  # broadcasting
    return Yskip.T.dot(CY)

%timeit ytcy_fast(Y)

# For comparison: all-sparse matrices
C_sparse = sparse.spdiags([Cdiag], [0], n, n)
Y_sparse = sparse.csr_matrix(Y)
%timeit Y_sparse.T.dot(C_sparse * Y_sparse)

    我的时间：

    In [59]: %timeit ytcy_fast(Y)
100 loops, best of 3: 16.1 ms per loop

In [18]: %timeit Y_sparse.T.dot(C_sparse * Y_sparse)
1 loops, best of 3: 282 ms per loop
    回复于
  • 21

    首先，您确定需要在问题中执行完整矩阵求逆吗？大多数情况下，只需要计算x = A ^ -1 y，这是一个更容易解决的问题 .

    如果确实如此，我会考虑计算逆矩阵的近似值而不是全矩阵求逆 . 因为矩阵求逆真的很昂贵 . 例如，参见Lanczos algorithm，以获得逆矩阵的有效近似 . 近似值可以稀疏地存储为奖励 . 此外，它只需要矩阵向量运算，因此您甚至不必将完整矩阵存储为逆矩阵 .

    作为替代方案，使用pyoperator，您还可以使用.todense方法来使用有效的矩阵向量运算来计算矩阵以进行逆运算 . 对角矩阵有一个特殊的稀疏容器 .

    有关Lanczos算法的实现，您可以查看pyoperators（免责声明：我是此软件的共同作者之一） .

    回复于

相关问题