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 [ ]: