有效地乘法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之间变化。作为优化过程的一部分，这个方程将被使用近1亿次，因此必须处理得非常快。

Y是随机初始化的，C是一个非常稀疏的矩阵，在300k的对角线上只有几个数不同于0。由于Numpy的对角线函数创建了稠密矩阵，所以我创建了C作为稀疏csr矩阵。但是当你试图解方程的第一部分时：

r = dot(C, Y)

由于内存限制，计算机崩溃。我决定将Y转换成csr_matrix并执行相同的操作：

r = dot(C, Ysparse)

该方法耗时1.38ms。但是这个解决方案有点“棘手”，因为我使用稀疏矩阵来存储稠密矩阵，我想知道这到底有多有效。

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

谢谢你的帮助！

Answer 1

尝试：

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