sparse matrix failed with element-wise multiplication using numpy.multiply() (Trac #1042)

[scipy-gitbot]

Original ticket http://projects.scipy.org/scipy/ticket/1042 on 2009-11-04 by trac user dingle, assigned to @wnbell.

If a and b are sparse matrices as follows:

>>> a = array([1,0,2])
>>> b = array([2,3,4])
>>> asp = sparse.lil_matrix(a)
>>> bsp = sparse.lil_matrix(b)
>>> multiply(asp, bsp)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.6/dist-packages/scipy/sparse/base.py", line 287, in __mul__
    raise ValueError('dimension mismatch')
ValueError: dimension mismatch

But a.multiply(b) works fine. I think they should be identical.

[scipy-gitbot]

trac user cowlicks wrote on 2013-04-18

It looks like when np.multiply is passed an ndarray objects it first tries to do matrix multiplication on them, if this fails due to a 'dimension mismatch', it does pointwise multiplication (element by element).

However when np.multiply is passed a spmatrix object it asks spmatrix.__mul__ what to do. Where it tries to do matrix multiplication, but if this has a dimension mismatch fails it just coughs up the value error seen. I'll try to patch this tomorrow by falling back on pointwise multiplication when there is a dimension mismatch, like the ndarrays do.

In [8]: np.multiply(A, B)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-8-ab35b70fd7b3> in <module>()
----> 1 np.multiply(A, B)

/usr/lib/python2.7/dist-packages/scipy/sparse/base.pyc in __mul__(self, other)
    244         if issparse(other):
    245             if self.shape[1] != other.shape[0]:
--> 246                 raise ValueError('dimension mismatch')
    247             return self._mul_sparse_matrix(other)
    248 

ValueError: dimension mismatch

Also, I realize this is an old ticket, but you say a.multiply(b) works. However standard nd.arrays don't have a multiply method... So maybe I'm missing something?

[scipy-gitbot]

trac user izzycecil wrote on 2013-04-18

I think they must have meant asp.multiply(bsp).

I was tinkering around with this bug last night, and found a number of odd behaviors.
With the following environment...

>>> import numpy as np
>>> from scipy import sparse
>>> a = np.array([1,2,3])
>>> b = np.array([1,0,2])
>>> asp = sparse.lil_matrix(a)
>>> bsp = sparse.lil_matrix(b)
>>> c = np.matrix([1,2,3])
>>> d = np.matrix([1,0,2])

We have this known fail case

>>> np.multiply(asp,bsp)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/izzy/.virtualenvs/scipy2.7/lib/python2.7/site-packages/scipy/sparse/base.py", line 262, in __mul__
    raise ValueError('dimension mismatch')
ValueError: dimension mismatch

But then, I looked at multiplying a sparse matrix by an array...

>>> np.multiply(asp,b)
array([ <1x3 sparse matrix of type '<type 'numpy.int64'>'
        with 3 stored elements in LInked List format>,
       <1x3 sparse matrix of type '<type 'numpy.int64'>'
        with 0 stored elements in LInked List format>,
       <1x3 sparse matrix of type '<type 'numpy.int64'>'
        with 3 stored elements in LInked List format>], dtype=object)

Where if a normal matrix was multiplied by an array,

>>> np.multiply(c,a)
matrix([[1, 4, 9]])

I would think we would want sparse matrices to simply work like a normal matrix. Or was this indeed the desired behavior? np.multiply(asp,d) will actually return the same result, with type matrix instead of array.

Then I looked at spmatrix.multiply...

>>> asp.multiply(b)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/izzy/.virtualenvs/scipy2.7/lib/python2.7/site-packages/scipy/sparse/base.py", line 215, in multiply
    return self.tocsr().multiply(other)
  File "/home/izzy/.virtualenvs/scipy2.7/lib/python2.7/site-packages/scipy/sparse/compressed.py", line 245, in multiply
    raise ValueError('inconsistent shapes')
ValueError: inconsistent shapes
>>> asp.multiply(bsp)
<1x3 sparse matrix of type '<type 'numpy.int64'>'
        with 2 stored elements in Compressed Sparse Row format>
>>> asp.multiply(d)
matrix([[1, 0, 6]])

It's not playing nice with Array's. This is because of the way it looks at the shape of "other" --- an easy fix.

Looking through the code, I have a pretty good understanding of why all of these are happening, but I'm confused on what the desired behaviors should be. We want spmatrix to essentially act as a matrix, yes? It also seems like there should be more separation between point-wise and matrix multiplication. Am I being nieve?

[scipy-gitbot]

trac user cowlicks wrote on 2013-04-18

Hello Izzy, by [http://docs.scipy.org/doc/numpy/reference/generated/numpy.multiply.html definition] the behavior of numpy.multiply should be pointwise multiplication. So that should be the expected behavior for this.

However it looks like numpy.multipy calls spmatrix.__mul__ which does matrix multiplication by default. I don't know how to fix this without changing __mul__ to pointwise mult which would break things.

I'm submitting a patch for the sparse.multiply 'inconsistent shapes' bug that you noted. so it should work now. If we could somehow redirect numpy.multiply on spmatrix objects to use sparse.multiply I think that would correct the issue initially reported.

nump.matrix has the correct behavior, so we could probably just copy what it does:

In [10]: Am, Bm
Out[10]: 
(matrix([[4, 1, 1],
        [4, 2, 4],
        [4, 3, 3]]),
 matrix([[4, 4, 1],
        [0, 3, 1],
        [1, 2, 0]]))

In [11]: np.multiply(Am, Bm)
Out[11]: 
matrix([[16,  4,  1],
        [ 0,  6,  4],
        [ 4,  6,  0]])

In [12]: Am*Bm
Out[12]: 
matrix([[17, 21,  5],
        [20, 30,  6],
        [19, 31,  7]])

Otherwise, it looks like there is plenty of work for us!

[scipy-gitbot]

@pv wrote on 2013-04-19

The problem with np.multiply here is that np.array(some_sparse_matrix) does something pretty weird:

>>> import numpy as np
>>> from scipy.sparse import lil_matrix
>>> np.array(lil_matrix([[1, 0], [0, 1]]))
array(<2x2 sparse matrix of type '<type 'numpy.int32'>'
    with 2 stored elements in LInked List format>, dtype=object)

Note that it's an 0-dim object array containing the lil_matrix.

One fix is to add an __array__ method that just does self.todense() to the bottom-most sparse matrix base class. This would make np.multiply(asp, b) first cast asp to a dense matrix before doing anything else.

That's a quick and dirty way to fix up the interoperability a bit. A better way would be to fix Numpy so that it does something sensible, such as optionally falling back to __elmul__ and __relmul__ methods inside np.multiply --- it already does some checking, but I think this has not been well thought out so far (see e.g. numpy/core/src/umath/ufunc_object.c:1908).

[pv]

Fixed in gh-2869