NunodeSousa NunodeSousa -4 years ago 330
Python Question

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?

Answer Source

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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download