[SciPy-dev] sparse matrix multiplication - poor performance

[Jonathan Guyer]

On Dec 6, 2006, at 5:16 PM, Nathan Bell wrote:

> There appears to be something terribly wrong with the current
> implementation of sparse matrix multiplication.  Consider the test
> below done in both MATLAB and scipy.
>
> MATLAB is  [660.5872  936.8530  702.5506  650.2876] times faster for
> these four problems.  Changing N = 100,000 puts the ratio in the
> neighborhood of 4000.
>
> Also, doing "B = A*A.T" takes *hundreds* of times longer than "B =
> A*A" or "B = A.T*A".  This is more than the time required for dense
> matrix multiplication (try N=1000).
>
> Can anyone comment on this?  Is the SPARSKIT code really that slow?

I don't know the reason for it, but I can verify the results (the  
scipy side of it, anyway). About a year ago, I posted some  
benchmarking results to <http://www.scipy.org/wikis/featurerequests/ 
SparseSolvers>, the upshot of which was that scipy.sparse was much  
slower than PySparse. I have no idea where that wiki page is since  
the site moved, though.


On my 1.67 GHz PowerBook, scipy gives

--------------------SCIPY OUTPUT--------------------

5.61
1.95
28.74
5.64

whereas the PySparse equivalent is comparable to your MATLAB results:

-------------------PYSPARSE CODE--------------------

import spmatrix
from scipy import *
from Numeric import *
import time
N = 10*1000

A = spmatrix.ll_mat(N, N, N)
ids = arange(N)
A.put(ones(N), ids, ids)
start = time.clock()
B = spmatrix.dot(A, A)
print time.clock() - start
start = time.clock()
B = spmatrix.matrixmultiply(A, A)
print time.clock() - start


A = spmatrix.ll_mat(N, N, 3*N)
ids1 = array((arange(N)-1, arange(N), arange(N)+1))
ids2 = array((arange(N), arange(N), arange(N)))
A.put(array(rand(3,N)).flat, ids1.flat, ids2.flat)
start = time.clock()
B = spmatrix.dot(A, A)
print time.clock() - start
start = time.clock()
B = spmatrix.matrixmultiply(A, A)
print time.clock() - start

------------------PYSPARSE OUTPUT-------------------

0.0
0.0
0.01
0.02