In [1]:

    

import numpy as np
from scipy import sparse
from numba import jit


    

In [2]:

    

n, m = 5, 5
A = sparse.rand(n, m, density=0.4, format='csr')


    

In [3]:

    

print A.toarray()
print A.data
print 'indeces =', A.indices
print 'indptr  =', A.indptr


    

    




[[ 0.          0.95992723  0.          0.          0.        ]
 [ 0.12104877  0.          0.          0.54592141  0.59000769]
 [ 0.          0.          0.79035446  0.6407521   0.        ]
 [ 0.          0.          0.          0.02540113  0.93112072]
 [ 0.75999549  0.          0.          0.          0.59062903]]
[ 0.95992723  0.12104877  0.54592141  0.59000769  0.79035446  0.6407521
  0.02540113  0.93112072  0.75999549  0.59062903]
indeces = [1 0 3 4 2 3 3 4 0 4]
indptr  = [ 0  1  4  6  8 10]

In [4]:

    

@jit(nopython=True)
def _csr_vecmat(num_rows, x, data, indices, indptr, out):
    for i in range(num_rows):
        for k in range(indptr[i], indptr[i+1]):
            out[indices[k]] += x[i] * data[k]


    

In [5]:

    

@jit(nopython=True)
def _csr_matvec(num_rows, data, indices, indptr, x, out):
    for i in range(num_rows):
        dot_prod = 0
        for k in range(indptr[i], indptr[i+1]):
            dot_prod += data[k] * x[indices[k]]
        out[i] = dot_prod


    

In [6]:

    

@jit(nopython=True)
def _csc_vecmat(num_cols, x, data, indices, indptr, out):
    for j in range(num_cols):
        dot_prod = 0
        for k in range(indptr[j], indptr[j+1]):
            dot_prod += data[k] * x[indices[k]]
        out[j] = dot_prod


    

In [7]:

    

x = np.arange(n)


    

In [8]:

    

y = np.empty(m)
A_csc = sparse.csc_matrix(A)
_csc_vecmat(m, x, A_csc.data, A_csc.indices, A_csc.indptr, y)
print y


    

    




[ 3.16103073  0.          1.58070893  1.90362899  5.74588597]

In [9]:

    

x.dot(A.toarray())


    

    
Out[9]:





array([ 3.16103073,  0.        ,  1.58070893,  1.90362899,  5.74588597])

In [10]:

    

y = np.empty(m)
_csr_vecmat(n, x, A.data, A.indices, A.indptr, y)
print y


    

    




[ 3.16103073  0.          1.58070893  1.90362899  5.74588597]

In [11]:

    

y = np.empty(m)
A_csc = sparse.csc_matrix(A)
%timeit _csc_vecmat(m, x, A_csc.data, A_csc.indices, A_csc.indptr, y)

y = np.empty(m)
%timeit _csr_vecmat(n, x, A.data, A.indices, A.indptr, y)


    

    




1000000 loops, best of 3: 611 ns per loop
1000000 loops, best of 3: 627 ns per loop

In [12]:

    

def csr_vecmat(x, A):
    n, m = A.shape
    if len(x) != n:
        raise ValueError('dimensions do not match')
    y = np.zeros(m)
    
    A_csc = sparse.csc_matrix(A)
    _csc_vecmat(m, x, A_csc.data, A_csc.indices, A_csc.indptr, y)
    return y


    

In [13]:

    

csr_vecmat(x, A)


    

    
Out[13]:





array([ 3.16103073,  0.        ,  1.58070893,  1.90362899,  5.74588597])

In [14]:

    

A.T.dot(x)


    

    
Out[14]:





array([ 3.16103073,  0.        ,  1.58070893,  1.90362899,  5.74588597])

In [15]:

    

x.dot(A.toarray())


    

    
Out[15]:





array([ 3.16103073,  0.        ,  1.58070893,  1.90362899,  5.74588597])

In [16]:

    

%timeit csr_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




10000 loops, best of 3: 123 µs per loop
10000 loops, best of 3: 65 µs per loop
10000 loops, best of 3: 84.1 µs per loop

In [17]:

    

n = 10**2
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)
%timeit csr_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




10000 loops, best of 3: 148 µs per loop
10000 loops, best of 3: 68.2 µs per loop
10000 loops, best of 3: 123 µs per loop

In [18]:

    

n = 10**3
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)


    

In [19]:

    

%timeit csr_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




100 loops, best of 3: 7.61 ms per loop
1000 loops, best of 3: 470 µs per loop
100 loops, best of 3: 6.42 ms per loop

In [20]:

    

n = 10**4
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)


    

In [21]:

    

%timeit csr_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




1 loops, best of 3: 1.51 s per loop
10 loops, best of 3: 49.4 ms per loop
1 loops, best of 3: 946 ms per loop

In [22]:

    

def csc_vecmat(x, A):
    n, m = A.shape
    if len(x) != n:
        raise ValueError('dimensions do not match')
    y = np.zeros(m)
    
    _csc_vecmat(m, x, A.data, A.indices, A.indptr, y)
    return y


    

In [23]:

    

n = 10**3
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)
%timeit csc_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




1000 loops, best of 3: 762 µs per loop
1000 loops, best of 3: 454 µs per loop
100 loops, best of 3: 6.11 ms per loop

In [24]:

    

n = 10**4
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)
%timeit csc_vecmat(x, A)
%timeit A.T.dot(x)
%timeit x.dot(A.toarray())


    

    




10 loops, best of 3: 84.8 ms per loop
10 loops, best of 3: 42.9 ms per loop
1 loops, best of 3: 935 ms per loop

In [25]:

    

n = 10**3
m = n
A = sparse.rand(n, m, density=0.4, format='csr')
x = np.ones(n)

y = np.empty(m)
A_csc = sparse.csc_matrix(A)
%timeit _csc_vecmat(m, x, A_csc.data, A_csc.indices, A_csc.indptr, y)

y = np.empty(m)
%timeit _csr_vecmat(n, x, A.data, A.indices, A.indptr, y)


    

    




1000 loops, best of 3: 751 µs per loop
1000 loops, best of 3: 804 µs per loop

In [25]: