I looked into what Scipy actually does here. There are several reasons it's faster than either Matlab or Julia in this case. First, the transpose of a CSR or CSC matrix in Scipy is nearly free, as it just returns the same exact same index/value data but as the opposite type. Second, sparse matrix-matrix multiplication in Scipy doesn't guarantee sorted indices in the output. So Scipy has implementations for every possible input combination of different sparse matrix formats, and whether or not the operation requires the inputs to have sorted indices. It's a lot of bookkeeping and code, but leads to good performance (I think the calculations are all defined in this templated header https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h) since they always use the right operation for the given inputs.
At the moment Julia only has CSC format, and expects sorted indices everywhere, so some operations are slower - sometimes significantly. Hooking up to a library like PETSc could potentially alleviate this. To get a fair comparison in Scipy for the operations Julia is doing today, you need to modify your Scipy code (which also ran for me in 5.8 seconds, so our laptops are about the same speed). Changing it to y = (x*x.transpose()).sort_indices() slows Scipy down to 22 seconds. Going another step to y = (x*x.transpose().tocsc()).sort_indices() didn't make much difference, x is a lot sparser than y here. Have a look at section 2.8 of Tim Davis' book "Direct Methods for Sparse Linear Systems" if you get a chance. At the end of the section he notes that when C has many more nonzeros than A or B, it's actually better to compute C = (B^T * A^T)^T than it is to compute C = ((A*B)^T)^T which is what Julia's doing now to sort the indices. Making this modification speeds Julia up a little bit from 23.9 seconds to 20.3 seconds on my laptop with your test case. Not quite Matlab speeds (about 12 sec, might be going to MKL here?), but slightly faster than SciPy when you're comparing apples-to-apples on truly the same operations. I'd expect the output to frequently be denser than the inputs in sparse matmul, but there are probably some cases where the current double-transpose is better. For operations that look something like x*x.', I wonder whether Julia's parser would be able to recognize that the two inputs are the same as one another, so it's known in advance that the output will be symmetric and we can skip half the work? -Tony On Friday, February 21, 2014 10:38:53 PM UTC-8, Tony Kelman wrote: > > Is the right way to rectify the global scope problem by wrapping the > benchmark code in a function and running that? > > If so, I'm still able to replicate this performance comparison with > similar results. Matlab's about 2x slower than Scipy which surprised me a > little, but Scipy probably has a specialized operator for sparse A_mul_Bt > whereas Matlab's parser might not be that smart. Doesn't look like anyone > has written sparse matmul with transposes in linalg/sparse.jl yet, but I > would find that useful enough that I'll write a few of them myself at some > point if nobody else does. > > About 65% of the time in your test appears to be spent in the > double-transpose at the end of sparse matmul (lines 175 and 176 of > linalg/sparse.jl). > > -Tony > > On Friday, February 21, 2014 6:48:44 PM UTC-8, John Myles White wrote: >> >> Are you timing these in the global scope? That will cause a substantial >> performance loss. >> >> — John >> >> On Feb 21, 2014, at 6:18 PM, Michael Schnall-Levin <[email protected]> >> wrote: >> >> > I've been doing some benchmarking of Julia vs Scipy for sparse matrix >> multiplication and I'm finding that julia is significantly (~4X - 5X) >> faster in some instances. >> > >> > I'm wondering if I'm doing something wrong, or if this is really true. >> Below are some code snippets for Julia and python. Any help would be very >> appreciated! >> > >> > ----- Julia: >> > Elapsed Time on my laptop: 24.9 seconds ----- >> > x_inds = Int[] >> > y_inds = Int[] >> > vals = Int[] >> > >> > for n = 1:10000 >> > inds = rand(1:2000,10,1) >> > for ind in inds >> > push!(x_inds, ind) >> > push!(y_inds, n) >> > push!(vals,1) >> > end >> > end >> > >> > x = sparse(x_inds, y_inds, vals, 2000, 10000) >> > >> > t = time() >> > for j = 1:250 >> > y = x*transpose(x) >> > end >> > print(string(time() - t, "\n")) >> > ----- >> > >> > ---- Python Elapsed Time on my laptop: 5.8 seconds ----- >> > import numpy >> > import scipy.sparse >> > import time >> > >> > x_inds = [] >> > y_inds = [] >> > vals = [] >> > for n in xrange(10000): >> > inds = numpy.random.randint(0, 2000,10) >> > >> > for ind in inds: >> > x_inds.append(ind) >> > y_inds.append(n) >> > vals.append(1) >> > >> > x_inds = numpy.array(x_inds) >> > y_inds = numpy.array(y_inds) >> > vals = numpy.array(vals) >> > >> > x = scipy.sparse.csc_matrix((vals, (x_inds, y_inds)), shape=(2000, >> 10000)) >> > >> > t = time.time() >> > for j in xrange(250): >> > y = x*x.transpose() >> > print time.time() - t >> > >> >>
