Numpy matrix product - sparse matrices

1
votes

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?

2
pythonnumpysparse-matrix
%time csr_matrix(B).dot(A) does first a conversion of B to a csr_matrix. Is that something you want to time as well? - jotasi
But B is not a sparse matrix. csr_matrix should only be applied to sparse matrices. - nunodsousa
I would say go with the sparse matrix method only when the matrix is extremely big and sparse. np.dot is extremely efficient under normal scenarios. - Divakar
Probably I didn't explain my problem correctly. I want to do np.dot(B,A) taking advantage of A being sparse. - nunodsousa
See my answer for a more elaborate explanation. - jotasi

2 Answers

3
votes

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
2
votes

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.