[SciPy-dev] sparsetools - new and (hopefully) improved!

Nathan Bell wnbell at gmail.com
Fri Dec 29 17:33:03 CST 2006

I've rewritten sparse/sparsetools to speed up some of the slower bits
in current code.  The previous implementation of CSR/CSC matrix
addition, multiplication, and conversions were all very slow, often
orders of magnitude slower than MATLAB or PySparse.

In the new implementation:
- CSR<->CSC conversion is cheap (linear time)
- sparse addition is fast (linear time)
- sparse multiplication is fast (noticeably faster than MATLAB in my testing)
- CSR column indices do not need to be sorted (same for CSC row indices)
- the results of all the above operations contain no explicit zeros
- COO->CSR and COO->CSC is linear time and sums duplicate entries

This last point is useful for constructing stiffness/mass matrices.
Matrix multiplication is essentially optimal (modulo Strassen-like
algorithms) and is based on the SMMP algorithm I mentioned a while

I implemented the functions in C++ with SWIG bindings.  The C++
consists of a single header file with templated functions for each
operation.  SWIG makes it easy to instantiate templates for each data
type (float,double,long
double,npy_cfloat,npy_cdouble,npy_clongdouble).  For the complex types
I had some trouble at first, as they do no support arithmetic
operations out of the box, that is, A+B works when A and B are floats
or doubles, but not npy_cfloats or npy_cdoubles.  To remedy this I
wrote C++ operators for the complex types (in another header file) to
implement the necessary operations.

I used the SWIG example in the NumPy codebase (with some alterations),
so strided arrays should be handled correctly.  Also, when two
different types are used together (e.g. a CSR matrix of type float64
is added to a CSR matrix of type float32) SWIG overloading takes care
of the conversion.  It will determine that float32->float64 is a safe
cast and use the float64 version of the function.  Some care is needed
to ensure that the most economical function is chosen, but I think
this can be accomplished just by ordering the functions appropriately
(e.g float32 before float64).  Currently the indices are assumed to be
C ints.  If we wanted both 32bit and 64bit ints it would be possible
to template the index type as well.

One challenge in writing sparse matrix routines is that the size of
the output is not known in advance.  For example, multiplying two
sparse matrices together leads to an unpredicitable number of nonzeros
in the result.  To get around this I used STL vectors as SWIG ARGOUTs.
 This produces a nice pythonic interface for most of the functions.
For functions operating on dense 2D arrays I chose to use an INPLACE
array, since the size is known in advance.

The code is available here:

The sparse directory in the archive is indented to replace the sparse
directory in scipy/Lib/
I've successfully installed the above code with the standard scipy
setup.py and all the unittests pass.

A simple benchmark for sparse matrix multiplication SciPy and MATLAB
is available here:

On my laptop SciPy takes about 0.33 seconds to MATLAB's 0.55, so 40% faster.

The code was all written my myself, aside from the bits pilfered from
the NumPy SWIG examples and SciPy UMFPACK SWIG wrapper, so there
shouldn't be any license issues.

I welcome any comments, questions, and criticism you have.  If you
know a better or more robust way to interface the C++ with Scipy
please speak up!.

Also, who should I contact regarding SVN commit privileges?

Nathan Bell wnbell at gmail.com

More information about the Scipy-dev mailing list