CODEWITHSUNDEEP
pythonnumpysparse-matrix
nunodsousa
nunodsousa

Reputation: 2785

Numpy matrix product - sparse matrices

Let us consider a matrix A as a diagonal matrix and a matrix B a random matrix, both with size N x N. We want to use the sparse properties of the matrix A to optimize the dot product, i.e. dot(B,A).

However if we compute the product using the sparcity properties of the matrix A, we cannot see any advantage (and it is much slower).

import numpy as np
from scipy.sparse import csr_matrix
# Matrix sizes
N = 1000

#-- matrices generation --
A = np.zeros((N,N), dtype=complex)
for i in range(N):
    A[i][i] = np.random.rand()
B = np.random.rand(N,N)

#product
%time csr_matrix(B).dot(A)
%time np.dot(B,A)

Results:

CPU times: user 3.51 s, sys: 8 ms, total: 3.52 s Wall time: 3.74 s

CPU times: user 348 ms, sys: 0 ns, total: 348 ms Wall time: 230 ms

How to do it correctly?

Upvotes: 1

Views: 1240

Answers (2)

hpaulj
hpaulj

Reputation: 231385

Done right, sparse dot is faster - if the matrices are indeed sparse. But you can't just throw the arrays into the csr_matrix.dot function. 

In [68]: N=1000
In [69]: from scipy import sparse
In [70]: A=np.eye(N)         # the diagonal is more interesting than all zeros
In [71]: B=np.random.rand(N,N)

Base case - dense matrix product

In [72]: timeit np.dot(B,A)
10 loops, best of 3: 98.8 ms per loop

This time is the same for all arrays of the same size (e.g dot(B,B), dot(A,A)).

Make sparse matrix from both. As has lots of zeros, Bs has none, but it is in sparse format

In [73]: As=sparse.csr_matrix(A)
In [74]: Bs=sparse.csr_matrix(B)

Note the conversion times; they are not trivial

In [101]: timeit sparse.csr_matrix(A)
100 loops, best of 3: 13.8 ms per loop
In [102]: timeit sparse.csr_matrix(B)
10 loops, best of 3: 50.1 ms per loop

Matrix product with the csr matrices can be faster. I'll use the Bs.dot(As) form because it's clearer. Bs*As and np.dot(Bs,As) are equivalent. But don't try np.dot(Bs,A) 

In [107]: timeit Bs.dot(As)
100 loops, best of 3: 19 ms per loop

In [112]: timeit sparse.csr_matrix(B).dot(sparse.csr_matrix(A)).A
10 loops, best of 3: 94.1 ms per loop

Noticeably better than the dense version, but marginally better if we include the conversion times.

But note that times vary widely depending on the sparsity of the matrices

In [108]: timeit As.dot(Bs)
100 loops, best of 3: 10 ms per loop
In [109]: timeit As.dot(B)
100 loops, best of 3: 5.82 ms per loop
In [110]: timeit As.dot(As)
1000 loops, best of 3: 215 µs per loop
In [111]: timeit Bs.dot(Bs)
1 loop, best of 3: 3.83 s per loop

Upvotes: 3

jotasi
jotasi

Reputation: 5177

The difference stems from the fact that you convert B to a sparse matrix during the timing (minor effect) and even worse, that dot is not aware of the fact, that A is sparse. If you were to do the conversion before the dot product the sparse dot product is actually faster:

import numpy as np
from scipy.sparse import csr_matrix
# Matrix sizes
N = 1000

#-- matrices generation --
A = np.zeros((N,N), dtype=complex)
for i in range(N):
    A[i][i] = np.random.rand()
B = np.random.rand(N,N)

Asparse = csr_matrix(A)
Bsparse = csr_matrix(B)

%timeit np.dot(B, A)
%timeit csr_matrix(B).dot(A)
%timeit Bsparse.dot(A)
%timeit csr_matrix.dot(B, Asparse)
%timeit csr_matrix.dot(Bsparse, Asparse)

Gives:
np.dot(B, A): 1 loop, best of 3: 414 ms per loop
csr_matrix(B).dot(A): 1 loop, best of 3: 2.22 s per loop
Bsparse.dot(A): 1 loop, best of 3: 2.17 s per loop
csr_matrix.dot(B, Asparse): 10 loops, best of 3: 32.5 ms per loop
csr_matrix.dot(Bsparse, Asparse): 10 loops, best of 3: 28 ms per loop

As you can see the sparse dot product is much faster in all cases where A is in a sparse matrix format which is making dot aware of the fact, that A is sparse. In your timing, the function is actually doing a conversion of B to a sparse format and then a dot product with a dense matrix A.

Upvotes: 2

Related Questions