4

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

在 python 中执行此操作的最简单的方法是 R.multiply(X Theta.T),但这会导致计算占用太多内存的完整矩阵 X Theta.T(m x n,顺序 100000**2),然后转储大部分条目,因为 R 是稀疏的。

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]) )

有谁知道一种快速、完全稀疏的方法来做到这一点?

4

3 回答 3

3

我针对您的问题分析了 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
于 2013-09-14T01:35:33.030 回答
1

根据评论中的额外信息,我认为让你失望的是对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
于 2013-09-13T20:12:36.870 回答
0

尚未测试 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]) )

这完全是稀疏的,甚至比原来的(对我来说内存效率低)解决方案运行得更快。

于 2013-09-13T20:28:46.927 回答