Thank you for leading the investigation. After all it appears that things are not that bad. Regarding your last question, have a look at line 104 in linalg/matmul.jl. It not on the parser level, but here A_mul_Bt checks if A and B are the same matrix and calls syrk in that case. I guess we could do something similar for sparse matrices.
2014-02-22 11:23 GMT+01:00 Tony Kelman <[email protected]>: > 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 >>> > >>> >>> -- Med venlig hilsen Andreas Noack Jensen
