Benchmarking different ways of updating a sparse matrix

Problem: We want to build a new sparse matrix by replacing or appending a column/row. How can we do so efficiently? We then want to use these matrices for fast matrix vector multiplication potentially requiring a cast to another sparse matrix type. What is the best overall compromise?



In [1]:

    
import numpy as np
import scipy.sparse



In [2]:

    
sparsity = 0.001
N, M = 10000, 200
rowvec = np.ones(N)
colvec = np.ones(M)
Aorig = np.random.random((N, M)) < sparsity



In [3]:

    
%timeit scipy.sparse.csr_matrix(Aorig)









    



The slowest run took 4.20 times longer than the fastest. This could mean that an intermediate result is being cached.
100 loops, best of 3: 6.38 ms per loop



In [4]:

    
%timeit scipy.sparse.lil_matrix(Aorig)









    



100 loops, best of 3: 18.6 ms per loop



In [5]:

    
Alil = scipy.sparse.lil_matrix(Aorig)
%timeit scipy.sparse.csr_matrix(Alil)









    



100 loops, best of 3: 3.19 ms per loop



In [6]:

    
Adok = scipy.sparse.dok_matrix(Aorig)
%timeit scipy.sparse.csr_matrix(Adok)









    



1000 loops, best of 3: 968 µs per loop

replace row



In [7]:

    
index = 3
newrow = np.random.random(N) < sparsity



In [8]:

    
Arow = Aorig[:, :]
%timeit Arow[:, index] = newrow









    



The slowest run took 4.07 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 32.9 µs per loop



In [9]:

    
def replacerow(A, rowindex, row):
    # see http://stackoverflow.com/questions/12129948/scipy-sparse-set-row-to-zeros
    # zero extant values
    A.data[A.indptr[rowindex]:A.indptr[rowindex+1]] = 0.0
    indices = np.nonzero(newrow)
    A[indices, rowindex] = newrow[indices]



In [10]:

    
Acsrrow = scipy.sparse.csr_matrix(Aorig)
replacerow(Acsrrow, index, newrow)
print(np.all(Acsrrow == Arow))
%timeit replacerow(Acsrrow, index, newrow)









    



True
The slowest run took 4.85 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 91.6 µs per loop



In [11]:

    
Acsrrow = scipy.sparse.csr_matrix(Aorig)
%timeit Acsrrow[:, index] = newrow[:, None]
print(np.all(Acsrrow == Arow))
print(Acsrrow.nnz, np.count_nonzero(Arow))









    



/home/andreas/miniconda2/envs/evolimmune/lib/python2.7/site-packages/scipy/sparse/compressed.py:730: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)






    



The slowest run took 339.45 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 323 µs per loop
True
(12035, 2048)



In [12]:

    
%timeit Acsrrow.eliminate_zeros()
Acsrrow.nnz, np.count_nonzero(Arow)









    



The slowest run took 4.00 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 29.5 µs per loop






    Out[12]:





(2048, 2048)



In [13]:

    
Alilrow = scipy.sparse.lil_matrix(Aorig)
%timeit Alilrow[:, index] = newrow[:, None]
print(np.all(Alilrow == Arow))
print(Alilrow.nnz, np.count_nonzero(Arow))









    



1000 loops, best of 3: 595 µs per loop
True
(2048, 2048)

append row



In [14]:

    
%timeit np.concatenate((Aorig, newrow[:, None]), axis=1)
Aprow = np.concatenate((Aorig, newrow[:, None]), axis=1)









    



1000 loops, best of 3: 332 µs per loop



In [15]:

    
Acsrprow = scipy.sparse.csr_matrix(Aorig)
print(np.all(Aprow==scipy.sparse.hstack((Acsrprow, scipy.sparse.csr_matrix(newrow[:, None])))))
%timeit scipy.sparse.hstack((Acsrprow, scipy.sparse.csr_matrix(newrow[:, None])))









    



True
1000 loops, best of 3: 624 µs per loop



In [16]:

    
def append(A, x):
    A._shape = (A.shape[0]+1, A.shape[1])
    A.indptr = np.hstack((A.indptr, A.indptr[-1]))
    A[-1, :] = x
    return A



In [17]:

    
newcol = np.random.random(M) < sparsity
np.all((append(scipy.sparse.csr_matrix(Aorig), newcol) == scipy.sparse.vstack((scipy.sparse.csr_matrix(Aorig), scipy.sparse.csr_matrix(newcol[None, :])))).todense())









    



/home/andreas/miniconda2/envs/evolimmune/lib/python2.7/site-packages/scipy/sparse/compressed.py:225: SparseEfficiencyWarning: Comparing sparse matrices using == is inefficient, try using != instead.
  " != instead.", SparseEfficiencyWarning)






    Out[17]:





True



In [18]:

    
Acsr = scipy.sparse.csr_matrix(Aorig)
%timeit append(Acsr, newcol)









    



1000 loops, best of 3: 5.06 ms per loop



In [19]:

    
Acsr = scipy.sparse.csr_matrix(Aorig)
%timeit scipy.sparse.vstack((Acsr, scipy.sparse.csr_matrix(newcol[None, :])))









    



1000 loops, best of 3: 268 µs per loop



In [20]:

    
def spmatreplacecol(A, colindex, newcol):
    row, col = A.nonzero()
    A.data[col == colindex] = 0.0
    indices = np.nonzero(newcol)
    A[indices, colindex] = newcol[indices]
    return A



In [36]:

    
colindex = 3
newcol = np.random.random(N) < sparsity



In [32]:

    
Acsr.shape









    Out[32]:





(10000, 200)



In [45]:

    
Acsr = scipy.sparse.csr_matrix(Aorig)
%time Acsr[:, colindex] = newcol[:, None]









    



CPU times: user 104 ms, sys: 0 ns, total: 104 ms
Wall time: 105 ms



In [44]:

    
Acsr2 = scipy.sparse.csr_matrix(Aorig)
%time spmatreplacecol(Acsr2, colindex, newcol)









    



CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 901 µs






    Out[44]:





<10000x200 sparse matrix of type '<type 'numpy.bool_'>'
	with 2060 stored elements in Compressed Sparse Row format>



In [47]:

    
np.all(Acsr2.todense() == Acsr)









    Out[47]:





True



In [ ]: