python - una - numpy tutorial español pdf

Subconjunto de una multiplicación de matriz, rápido y escaso

Question

Convertir un código de filtrado colaborativo para usar matrices dispersas Me resulta desconcertante el siguiente problema: dadas dos matrices completas X (m por l) y Theta (n por l), y una matriz dispersa R (m por n), ¿hay alguna forma rápida de calcular el producto interno escaso. Las dimensiones grandes son m y n (orden 100000), mientras que l es pequeño (orden 10). Esta es probablemente una operación bastante común para Big Data, ya que aparece en la función de costo de la mayoría de los problemas de regresión lineal, por lo que esperaría una solución integrada en scipy.sparse, pero no he encontrado nada obvio todavía.

La forma ingenua de hacer esto en python es R.multiply (X Theta.T), pero esto dará como resultado la evaluación de la matriz completa X Theta.T (m por n, orden 100000 ** 2) que ocupa demasiada memoria, luego volcando la mayoría de las entradas ya que R es escaso.

Ya hay una pseudo solución aquí en stackoverflow , pero no es dispersa en un solo paso:

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

Esto funciona bien, y rápido, para mí en matrices lo suficientemente pequeñas, pero muestra errores en mis grandes conjuntos de datos con el siguiente error:

... 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.

Una solución que ES realmente escasa, pero muy lenta, es:

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

¿Alguien sabe de una manera rápida y totalmente escasa de hacer esto?

Answer 1

Aún no he probado la respuesta de Jaime (¡gracias de nuevo!), Pero implementé otra respuesta que funciona mientras tanto usando cython.

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

Luego, compílelo de acuerdo con http://docs.cython.org/src/userguide/tutorial.html

Luego, en mi código python, incluyo lo siguiente:

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

Esto funciona completamente escaso y también se ejecuta más rápido que incluso la solución original (que no me sirve de memoria).

Answer 2

En base a la información adicional sobre los comentarios, creo que lo que te está desanimando es la llamada a np.unique . Pruebe el siguiente enfoque:

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

Obtengo los siguientes tiempos:

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

Y por supuesto, con n = 100 :

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

Answer 3

Hice un perfil de 4 soluciones diferentes para su problema, y ​​parece que para cualquier tamaño de la matriz, la solución de numba jit es la mejor. Un segundo cercano es la solución cython de @ Alexander.

Aquí están los resultados (M es el número de filas en la matriz 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.

El script que utilicé para perfilar estos métodos es:

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)

La función cython es una versión ligeramente modificada del código de @ 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