[SciPy-User] Sparse matrices and dot product

Nathaniel Smith njs@pobox....
Sun Nov 28 20:32:39 CST 2010


On Sun, Nov 28, 2010 at 7:16 AM, Pauli Virtanen <pav@iki.fi> wrote:
> However, I believe 'dot' should be left to be there. `ndarrays` recently
> gained the same method for matrix products, so it makes sense to leave it
> be also for sparse matrices. This has also the advantage that it becomes
> possible to write generic code that works both on sparse matrices and on
> ndarrays.

If it's been decided that ndarray's should have a dot method, then I
agree that sparse matrices should too -- for compatibility. But it
doesn't actually solve the problem of writing generic code. If A is
dense and B is sparse, then A.dot(B) still won't work.

I just spent a few minutes trying to work out if this is fixable by
defining a protocol -- you need like an  __rdot__ or something? -- but
didn't come up with anything I'd want to actually recommend.

(In fact, there are lots of other problems with writing generic code,
like the behavior of __mul__ and the way sparse matrices like to turn
all dense results into np.matrix's instead of np.ndarray's. The API
seems designed on the assumption that everyone will use np.matrix
instead of np.ndarray for everything, which I guess is fine, but since
I personally never touch np.matrix my generic code ends up pretty
ugly. I don't see how to do better without serious API breakage *and*
a lot more cooperation from numpy. The only full solution might be to
add sparse matrix support to numpy, and eventually deprecate
scipy.sparse?)

In the mean time, maybe it would be a good deed to add this to scipy.sparse?:

def spdot(A, B):
  "The same as np.dot(A, B), except it works even if A or B or both
might be sparse."
  if issparse(A) and issparse(B):
    return A * B
  elif issparse(A) and not issparse(B):
    return (A * B).view(type=B.__class__)
  elif not issparse(A) and issparse(B):
    return (B.T * A.T).T.view(type=A.__class__)
  else:
    return np.dot(A, B)

-- Nathaniel


More information about the SciPy-User mailing list