高效地相乘Numpy/Scipy稀疏矩阵和密集矩阵

Question

我正在努力实现以下方程：

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

Y是一个(n x f)的矩阵，C是一个(n x n)的对角矩阵；n大约为300k，而f将在100到200之间变化。作为优化过程的一部分，这个方程式会被使用近100百万次，因此必须处理得非常快。

Y是随机初始化的，而C是一个非常稀疏的矩阵，只有少数几个数在对角线上与0不同。由于Numpy的对角函数会创建密集矩阵，所以我将C创建为一个稀疏csr矩阵。但在尝试解决方程式的第一部分时：

r = dot(C, Y)

电脑由于内存限制而崩溃。我决定尝试将 Y 转换为 csr_matrix 并执行相同的操作：

r = dot(C, Ysparse)

这种方法花费了1.38毫秒。但是，这种解决方案有点“棘手”，因为我使用稀疏矩阵来存储密集矩阵，我想知道这样的效率如何。

所以我的问题是，如果有一种方法可以将稀疏矩阵C和密集矩阵Y相乘，而不必将Y转换为稀疏矩阵并提高性能吗？ 如果某种方式可以表示C作为对角线元素密集型而不消耗大量内存，那么这可能会导致非常高效的性能，但我不知道是否可能。

感谢您的帮助！

- Willian Fuks

---

只是出于好奇，这是为了一个推荐系统吗？（例如http://bit.ly/1aqCsfs） - Tobias Domhan

@TobiasDomhan 是的，我正在实现Koren关于隐式反馈数据集的论文 =)。 - Willian Fuks

Answer 1

当计算r = dot(C,Y)时，点积遇到内存问题的原因是因为numpy的dot函数没有本地支持处理稀疏矩阵。正在发生的事情是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)

- M.H.

---

既然这是一项校内项目，让我指出，“紧凑性”不应该成为目标。然而，你的第二个例子确实是更面向对象的方法，当然，这也是通常更好的方法。但它并不是关于变得简短。可读性比简洁更重要。 - Will

我主要是发布了一篇帖子，因为当天我也遇到了同样的问题，而且没有看到解释为什么dot（sparse，dense）函数没有返回您预期的结果。只希望给那些遇到这个问题的人提供一些帮助。 - M.H.

Answer 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

- pv.

---

谢谢帮忙！我尝试了一下，时间大致相同。目前最好的解决方案是将密集的Y表示为稀疏矩阵，这样速度更快。也许没有办法改进这个问题。 - Willian Fuks

这是一个很棒的解决方案！ :) - Tobias Domhan

Answer 3

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

如果确实如此，我建议您计算逆矩阵的近似值，而不是完整的矩阵求逆。因为矩阵求逆真的非常昂贵。例如，可以使用Lanczos算法来高效地近似逆矩阵。该近似值可以作为奖励稀疏存储。此外，它仅需要矩阵向量操作，因此您甚至无需存储完整矩阵以进行求逆。

另外，使用pyoperators，您还可以使用.todense方法通过有效的矩阵向量操作计算要求逆的矩阵。有一种特殊的对角线矩阵稀疏容器。

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

- Nicolas Barbey

---

谢谢提供这个信息，我一定会仔细研究！在这个问题中，反演过程的成本并不太高，在这里只需要258微秒（对于100 x 100矩阵，如果f = 100）。如果我能进一步减少这个时间，那就太好了，因为它将重复数百万次。主要瓶颈是乘法中的1.38毫秒，这需要几天才能完成。我会看看pyoperators并尝试一下！谢谢！ - Willian Fuks

我在你的PB上进行了一些测试，看起来最长的矩阵向量运算是Y.T * v。你的对角矩阵是否是正定的？这将允许将问题写成(B.T * B)^-1，其中B = sqrt(C) * Y。此外，在1亿次使用中，C或Y是否会发生变化？ - Nicolas Barbey

好的观点，我不确定C是否确实是一个正定矩阵，我认为不是（我必须进一步研究这个主题，我正在遵循Koren的论文，他没有提到这一点）。C是常数，在此过程中Y也是如此。只有在计算完所有Xu之后，Y才会发生变化，然后Y的更新过程开始。 - Willian Fuks

由于C是对角矩阵，如果所有元素>0，则它是正定的。 - Nicolas Barbey

Answer 4

我不知道在问题被提出时是否可能；但是现在，广播是您的朋友。一个n * n对角矩阵只需要成为对角线元素的数组即可用于矩阵乘积：

>>> n, f = 5, 3
>>> Y = np.random.randint(0, 10, (n, f))
>>> C = np.random.randint(0, 10, (n,))
>>> Y.shape
(5, 3)
>>> C.shape
(5,)
>>> np.all(Y.T @ np.diag(C) @ Y == Y.T*C @ Y)
True

请注意，Y.T*C @ Y 不是结合的：

>>> Y.T*(C @ Y)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (3,5) (3,)

但是 Y.T @ (C[:, np.newaxis]*Y) 将会得到预期的结果:

>>> np.all(Y.T*C @ Y == Y.T@(C[:, np.newaxis]*Y))
True

- Aubergine