
I'm trying to improve integration between SciPy's sparse matrices and NumPy's dense array/matrix objects. One problem I'm facing is that NumPy matrices redefine the * operator to call NumPy's dot() function. Since dot() has no awareness of SciPy's sparse matrix objects, this doesn't work for the operation 'dense * sparse'. (It does work for sparse * dense, which calls sparse.__mul__ first.) I'd like to propose the addition of a basic sparse matrix object to NumPy. This wouldn't need to provide any actual functionality, but could instead provide a skeletal interface or base class from which sparse matrices in other packages (SciPy and potentially PySparse) could derive. The spmatrix object in SciPy would be a good starting point. The benefit of this would be that hooks for proper handling of sparse matrices would be easy to provide for functions like dot(), where(), and var(). There may be other ways to make 'dense * sparse' work in SciPy, but I haven't been able to come up with any. This solution would at least be quite flexible and quite straightforward. -- Ed

On Mon, 27 Feb 2006, Ed Schofield wrote:
I'm trying to improve integration between SciPy's sparse matrices and NumPy's dense array/matrix objects. One problem I'm facing is that NumPy matrices redefine the * operator to call NumPy's dot() function. Since dot() has no awareness of SciPy's sparse matrix objects, this doesn't work for the operation 'dense * sparse'. (It does work for sparse * dense, which calls sparse.__mul__ first.)
Have you tried defining sparse.__rmul__? dense.__mul__ should raise an exception when it does not know about the rhs operant and then Python calls <rhs operant>.__rmul__. Pearu

Pearu Peterson wrote:
On Mon, 27 Feb 2006, Ed Schofield wrote:
I'm trying to improve integration between SciPy's sparse matrices and NumPy's dense array/matrix objects. One problem I'm facing is that NumPy matrices redefine the * operator to call NumPy's dot() function. Since dot() has no awareness of SciPy's sparse matrix objects, this doesn't work for the operation 'dense * sparse'. (It does work for sparse * dense, which calls sparse.__mul__ first.)
Have you tried defining sparse.__rmul__? dense.__mul__ should raise an exception when it does not know about the rhs operant and then Python calls <rhs operant>.__rmul__.
Yes, we've defined __rmul__, and this works fine for dense arrays, whose __mul__ raises an exception. The problem is that matrix.__mul__ calls dot(), which doesn't raise an exception, but rather creates an oddball object array: matrix([[ (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0]], dtype=object) We could potentially modify the __mul__ function of numpy's matrix objects to make a guess about whether an array constructed out of the argument will somehow be sane or whether, like here, it should raise an exception. But this would be difficult to get right, since the sparse matrix formats are quite varied (some supporting the map/sequence protocols, some not, etc.). But being able to test isinstance(arg, spmatrix) would make this easy. -- Ed

On Mon, 27 Feb 2006, Ed Schofield wrote:
Pearu Peterson wrote:
On Mon, 27 Feb 2006, Ed Schofield wrote:
I'm trying to improve integration between SciPy's sparse matrices and NumPy's dense array/matrix objects. One problem I'm facing is that NumPy matrices redefine the * operator to call NumPy's dot() function. Since dot() has no awareness of SciPy's sparse matrix objects, this doesn't work for the operation 'dense * sparse'. (It does work for sparse * dense, which calls sparse.__mul__ first.)
Have you tried defining sparse.__rmul__? dense.__mul__ should raise an exception when it does not know about the rhs operant and then Python calls <rhs operant>.__rmul__.
Yes, we've defined __rmul__, and this works fine for dense arrays, whose __mul__ raises an exception. The problem is that matrix.__mul__ calls dot(), which doesn't raise an exception, but rather creates an oddball object array:
matrix([[ (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0]], dtype=object)
We could potentially modify the __mul__ function of numpy's matrix objects to make a guess about whether an array constructed out of the argument will somehow be sane or whether, like here, it should raise an exception. But this would be difficult to get right, since the sparse matrix formats are quite varied (some supporting the map/sequence protocols, some not, etc.). But being able to test isinstance(arg, spmatrix) would make this easy.
Sure, isinstance(arg,spmatrix) would work but it is not a general solution for performing binary operations with matrices and such user defined objects that numpy is not aware of. But these objects may be aware of numpy matrices or arrays. Sparse matrix is one example. Other example is defining a symbolic matrix. Etc. So, IMHO matrix.__mul__ (or dot) should be fixed instead. Pearu

Pearu Peterson wrote:
On Mon, 27 Feb 2006, Ed Schofield wrote:
Yes, we've defined __rmul__, and this works fine for dense arrays, whose __mul__ raises an exception. The problem is that matrix.__mul__ calls dot(), which doesn't raise an exception, but rather creates an oddball object array:
matrix([[ (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0, (1, 0) 0.0 (2, 1) 0.0 (3, 0) 0.0]], dtype=object)
We could potentially modify the __mul__ function of numpy's matrix objects to make a guess about whether an array constructed out of the argument will somehow be sane or whether, like here, it should raise an exception. But this would be difficult to get right, since the sparse matrix formats are quite varied (some supporting the map/sequence protocols, some not, etc.). But being able to test isinstance(arg, spmatrix) would make this easy.
Sure, isinstance(arg,spmatrix) would work but it is not a general solution for performing binary operations with matrices and such user defined objects that numpy is not aware of. But these objects may be aware of numpy matrices or arrays. Sparse matrix is one example. Other example is defining a symbolic matrix. Etc. So, IMHO matrix.__mul__ (or dot) should be fixed instead.
Ah, yes, this could be the simplest solution (at least to the __mul__ problem). We could redefine matrix.__mul__ as def __mul__(self, other): if isinstance(other, N.ndarray) or not hasattr(other, '__rmul__') \ or N.isscalar(other): return N.dot(self, other) else: return NotImplemented This seems to fix multiplication. I may make a case later for sparse matrix hooks for other functions, but I don't see a pressing need right now. ;) -- Ed
participants (2)
-
Ed Schofield
-
Pearu Peterson