[SciPy-User] Sparse matrices and dot product

[Pauli Virtanen]

Sun, 28 Nov 2010 18:32:39 -0800, Nathaniel Smith wrote:
[clip]
> 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.

Yes, the .dot() method does not fully solve the generic code problem, and 
adding it was not motivated by that. However, in iterative methods you 
often only want to compute matrix-dense-vector products, and for that 
particular case it happens to be enough.

> 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.

The no-brain approach would probably be to ape what Python is doing, 
i.e., __dot__, __rdot__, and returning NotImplemented when the operation 
is not supported by a particular routine. I didn't try to think this 
fully out, though, so I don't know if there are some roadblocks.

> 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?

I would be +1 with having sparse *arrays* in Numpy. But that's expanding 
the scope of Numpy by quite a lot...

> 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)

I would't have objections to adding that.

-- 
Pauli Virtanen