5

If I have a numpy.ndarray A and a scipy.sparse.csc_matrix B, how do I take A dot B? I can do B dot A by saying B.dot(A), but the other way I can only think of this:

B.T.dot(A.T).T

Is there a more direct method to do this?

ali_m
  • 71,714
  • 23
  • 223
  • 298
shaoyl85
  • 1,854
  • 18
  • 30
  • 1
    Why not simply `A.dot(B)` (or `np.dot(A, B)`)? Am I missing something here? – ali_m Apr 27 '15 at 10:45
  • The dense `dot` probably can't handle a sparse matrix. But the sparse version will know about arrays. I wonder about the speed of this v.v `A.dot(B.A)` – hpaulj Apr 27 '15 at 11:07
  • @hpaulj Yes it will, although the result will of course be dense – ali_m Apr 27 '15 at 11:39
  • @hpaulj Ah, I see that it may be a version issue. The `dot` override was added in [this PR](https://github.com/scipy/scipy/pull/2869), and is not present in versions of scipy older than 0.14.0. – ali_m Apr 27 '15 at 11:47

1 Answers1

3

Your question initially confused me, since for my version of scipy, A.dot(B) and np.dot(A, B) both work fine; the .dot method of the sparse matrix simply overrides np.dot. However it seems that this feature was added in this pull request, and is not present in versions of scipy older than 0.14.0. I'm guessing that you have one of these older versions.

Here's some test data:

import numpy as np
from scipy import sparse

A = np.random.randn(1000, 2000)
B = sparse.rand(2000, 3000, format='csr')

For versions of scipy >= 0.14.0, you can simply use:

C = A.dot(B)
C = np.dot(A, B)

For versions < 0.14.0, both of these will raise a ValueError:

In [6]: C = A.dot(B)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-7fbaa337fd94> in <module>()
----> 1 C = A.dot(B)

ValueError: Cannot find a common data type.

Instead, you could use one of:

# your original solution
%timeit B.T.dot(A.T).T
# 10 loops, best of 3: 93.1 ms per loop

# post-multiply A by B
%timeit B.__rmul__(A)
# 10 loops, best of 3: 91.9 ms per loop

As you can see there's basically no performance difference, although I personally think the second version is more readable.

Update:

As @shaoyl85 just pointed out, one can just use the * operator rather than calling the __rmul__() method directly:

# equivalent to B.__rmul__(A)
C = A * B

It seems that matrices have a higher priority when determining the behavior of the * operator than ndarrays. This is a potential gotcha for those of us who are more used to the behavior of ndarrays (where * means elementwise multiplication).

ali_m
  • 71,714
  • 23
  • 223
  • 298
  • Turns out we visited this issue a year ago: http://stackoverflow.com/a/23144255/901925 – hpaulj Apr 27 '15 at 15:41
  • 1
    Interesting. I was indeed using an earlier version. And on seeing your answer I realize that I can simply do `A * B` which effectively calls the `__rmul__`. – shaoyl85 Apr 27 '15 at 16:36
  • Ah, it didn't occur to me to try `*` - I'm much more used to dealing with arrays rather than matrices (where `*` means element-wise multiplication). In fact this is a potential gotcha: it's not at all obvious that it's the matrix in the expression that determines the behavior of the `*` operator, rather than the ndarray. – ali_m Apr 27 '15 at 17:05
  • 1
    I think so. In this case the `__rmul__` is called simply because `A.__mul__(B)` returns `NotImplemented`. Otherwise the behavior may become an element-wise multiplication. – shaoyl85 Apr 27 '15 at 17:08