矩阵乘法的子集，快速且稀疏

Question

将协同过滤代码转换为使用稀疏矩阵我对以下问题感到困惑：给定两个完整矩阵 X（m x l）和 Theta（n x l），以及一个稀疏矩阵 R（m x n），有没有一种快速计算稀疏内积的方法。大尺寸是 m 和 n（100000 阶），而 l 很小（10 阶）。对于大数据来说，这可能是一个相当常见的操作，因为它出现在大多数线性回归问题的成本函数中，所以我希望 scipy.sparse 中有一个内置的解决方案，但我还没有发现任何明显的东西。

在 python 中执行此操作的简单方法是 R.multiply(XTheta.T)，但这将导致计算完整矩阵 XTheta.T（m x n，阶数 100000 **2) 占用太多内存，然后因为 R 是稀疏的，所以转储大部分条目。

有一个pseudo solution already here on stackoverflow，但是一步不稀疏：

def sparse_mult_notreally(a, b, coords):
    rows, cols = coords
    rows, r_idx = np.unique(rows, return_inverse=True)
    cols, c_idx = np.unique(cols, return_inverse=True)
    C = np.array(np.dot(a[rows, :], b[:, cols])) # this operation is dense
    return sp.coo_matrix( (C[r_idx,c_idx],coords), (a.shape[0],b.shape[1]) )

这对我来说在足够小的数组上运行良好且快速，但在我的大数据集上却出现以下错误：

... in sparse_mult(a, b, coords)
      132     rows, r_idx = np.unique(rows, return_inverse=True)
      133     cols, c_idx = np.unique(cols, return_inverse=True)
  --> 134     C = np.array(np.dot(a[rows, :], b[:, cols])) # this operation is not sparse
      135     return sp.coo_matrix( (C[r_idx,c_idx],coords), (a.shape[0],b.shape[1]) )

ValueError: array is too big.

一个实际上稀疏但非常慢的解决方案是：

def sparse_mult(a, b, coords):
    rows, cols = coords
    n = len(rows)
    C = np.array([ float(a[rows[i],:]*b[:,cols[i]]) for i in range(n) ]) # this is sparse, but VERY slow
    return sp.coo_matrix( (C,coords), (a.shape[0],b.shape[1]) )

有人知道一种快速、完全稀疏的方法吗？

Answer 1

我针对您的问题分析了 4 种不同的解决方案，看起来对于任何大小的数组，numba jit 解决方案都是最好的。紧随其后的是@Alexander 的 cython 解决方案。

这是结果（M 是x 数组中的行数）：

M = 1000
function sparse_dense    took 0.03 sec.
function sparse_loop     took 0.07 sec.
function sparse_numba    took 0.00 sec.
function sparse_cython   took 0.09 sec.
M = 10000
function sparse_dense    took 2.88 sec.
function sparse_loop     took 0.68 sec.
function sparse_numba    took 0.00 sec.
function sparse_cython   took 0.01 sec.
M = 100000
function sparse_dense    ran out of memory
function sparse_loop     took 6.84 sec.
function sparse_numba    took 0.09 sec.
function sparse_cython   took 0.12 sec.

我用来分析这些方法的脚本是：

import numpy as np
from scipy.sparse import coo_matrix
from numba import autojit, jit, float64, int32
import pyximport
pyximport.install(setup_args={"script_args":["--compiler=mingw32"],
                              "include_dirs":np.get_include()},
                  reload_support=True)

def sparse_dense(a,b,c):
    return coo_matrix(c.multiply(np.dot(a,b)))

def sparse_loop(a,b,c):
    """Multiply sparse matrix `c` by np.dot(a,b) by looping over non-zero
    entries in `c` and using `np.dot()` for each entry."""
    N = c.size
    data = np.empty(N,dtype=float)
    for i in range(N):
        data[i] = c.data[i]*np.dot(a[c.row[i],:],b[:,c.col[i]])
    return coo_matrix((data,(c.row,c.col)),shape=(a.shape[0],b.shape[1]))

#@autojit
def _sparse_mult4(a,b,cd,cr,cc):
    N = cd.size
    data = np.empty_like(cd)
    for i in range(N):
        num = 0.0
        for j in range(a.shape[1]):
            num += a[cr[i],j]*b[j,cc[i]]
        data[i] = cd[i]*num
    return data

_fast_sparse_mult4 = \
    jit(float64[:,:](float64[:,:],float64[:,:],float64[:],int32[:],int32[:]))(_sparse_mult4)

def sparse_numba(a,b,c):
    """Multiply sparse matrix `c` by np.dot(a,b) using Numba's jit."""
    assert c.shape == (a.shape[0],b.shape[1])
    data = _fast_sparse_mult4(a,b,c.data,c.row,c.col)
    return coo_matrix((data,(c.row,c.col)),shape=(a.shape[0],b.shape[1]))

def sparse_cython(a, b, c):
    """Computes c.multiply(np.dot(a,b)) using cython."""
    from sparse_mult_c import sparse_mult_c

    data = np.empty_like(c.data)
    sparse_mult_c(a,b,c.data,c.row,c.col,data)
    return coo_matrix((data,(c.row,c.col)),shape=(a.shape[0],b.shape[1]))

def unique_rows(a):
    a = np.ascontiguousarray(a)
    unique_a = np.unique(a.view([('', a.dtype)]*a.shape[1]))
    return unique_a.view(a.dtype).reshape((unique_a.shape[0], a.shape[1]))

if __name__ == '__main__':
    import time

    for M in [1000,10000,100000]:
        print 'M = %i' % M
        N = M + 2
        L = 10

        x = np.random.rand(M,L)
        t = np.random.rand(N,L).T

        # number of non-zero entries in sparse r matrix
        S = M*10

        row = np.random.randint(M,size=S)
        col = np.random.randint(N,size=S)

        # remove duplicate rows and columns       
        row, col = unique_rows(np.dstack((row,col)).squeeze()).T

        data = np.random.rand(row.size)

        r = coo_matrix((data,(row,col)),shape=(M,N))

        a2 = sparse_loop(x,t,r)

        for f in [sparse_dense,sparse_loop,sparse_numba,sparse_cython]:
            t0 = time.time()
            try:
                a = f(x,t,r)
            except MemoryError:
                print 'function %s ran out of memory' % f.__name__
                continue
            elapsed = time.time()-t0
            try:
                diff = abs(a-a2)
                if diff.nnz > 0:
                    assert np.max(abs(a-a2).data) < 1e-5
            except AssertionError:
                print f.__name__
                raise
            print 'function %s took %.2f sec.' % (f.__name__,elapsed)

cython 函数是@Alexander 代码的略微修改版本：

# working from tutorial at: http://docs.cython.org/src/tutorial/numpy.html
cimport numpy as np

# Turn bounds checking back on if there are ANY problems!
cimport cython
@cython.boundscheck(False) # turn of bounds-checking for entire function
def sparse_mult_c(np.ndarray[np.float64_t, ndim=2] a,
                  np.ndarray[np.float64_t, ndim=2] b,
                  np.ndarray[np.float64_t, ndim=1] data,
                  np.ndarray[np.int32_t, ndim=1] rows,
                  np.ndarray[np.int32_t, ndim=1] cols,
                  np.ndarray[np.float64_t, ndim=1] out):

    cdef int n = rows.shape[0]
    cdef int k = a.shape[1]
    cdef int i,j

    cdef double num

    for i in range(n):
        num = 0.0
        for j in range(k):
            num += a[rows[i],j] * b[j,cols[i]]
        out[i] = data[i]*num

Answer 2

根据有关 cmets 的额外信息，我认为让您失望的是对 np.unique 的调用。尝试以下方法：

import numpy as np
import scipy.sparse as sps
from numpy.core.umath_tests import inner1d

n = 100000
x = np.random.rand(n, 10)
theta = np.random.rand(n, 10)
rows = np.arange(n)
cols = np.arange(n)
np.random.shuffle(rows)
np.random.shuffle(cols)


def sparse_multiply(x, theta, rows, cols):
    data = inner1d(x[rows], theta[cols])
    return sps.coo_matrix((data, (rows, cols)),
                          shape=(x.shape[0], theta.shape[0]))

我得到以下时间：

n = 1000
%timeit sparse_multiply(x, theta, rows, cols)
1000 loops, best of 3: 465 us per loop

n = 10000
%timeit sparse_multiply(x, theta, rows, cols)
100 loops, best of 3: 4.29 ms per loop

n = 100000
%timeit sparse_multiply(x, theta, rows, cols)
10 loops, best of 3: 61.5 ms per loop

当然还有n = 100:

>>> np.allclose(sparse_multiply(x, theta, rows, cols).toarray()[rows, cols],
                x.dot(theta.T)[rows, cols])
>>> True

Answer 3

尚未测试 Jaime 的答案（再次感谢！），但我实现了另一个同时使用 cython 的答案。

文件sparse_mult_c.pyx：

# working from tutorial at: http://docs.cython.org/src/tutorial/numpy.html
cimport numpy as np

# Turn bounds checking back on if there are ANY problems!
cimport cython
@cython.boundscheck(False) # turn of bounds-checking for entire function
def sparse_mult_c(np.ndarray[np.float64_t, ndim=2] a,
                  np.ndarray[np.float64_t, ndim=2] b,
                  np.ndarray[np.int32_t, ndim=1] rows,
                  np.ndarray[np.int32_t, ndim=1] cols,
                  np.ndarray[np.float64_t, ndim=1] C ):

    cdef int n = rows.shape[0]
    cdef int k = a.shape[1]
    cdef int i,j

    for i in range(n):
        for j in range(k):
            C[i] += a[rows[i],j] * b[j,cols[i]]

然后按照http://docs.cython.org/src/userguide/tutorial.html进行编译

然后在我的 python 代码中，我包含以下内容：

def sparse_mult(a, b, coords):
    #a,b are np.ndarrays
    from sparse_mult_c import sparse_mult_c
    rows, cols = coords
    C = np.zeros(rows.shape[0])
    sparse_mult_c(a,b,rows,cols,C)
    return sp.coo_matrix( (C,coords), (a.shape[0],b.shape[1]) )

这完全是稀疏的，而且运行速度甚至比原来的（对我来说内存效率低）解决方案还要快。