18

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

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

Y 是 (nxf) 矩阵，C 是 (nxn) 对角矩阵；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。但是这个解决方案有点“棘手”，因为我使用稀疏矩阵来存储密集矩阵，我想知道这到底有多有效。

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

我感谢您的帮助！

4

4 回答 4

29

计算 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)
于 2013-05-25T22:26:38.127 回答
9

尝试：

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
于 2012-11-07T16:48:43.837 回答
2

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

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

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

对于 Lanczos 算法的实现，您可以查看pyoperators（免责声明：我是该软件的合著者之一）。

于 2012-11-07T15:46:38.380 回答
0

我不知道问这个问题时是否有可能；但现在，广播是你的朋友。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
于 2021-01-19T04:20:19.013 回答