BUG: sparse in-place operations is not in-place (change data-type and object)

[zerothi]

When doing inplace additions of sparse matrices the resulting data-type is not retained.

I.e. the in-place sparse matrix gets up-casted.

Reproducing code example:

import numpy as np
import scipy.sparse as sp

cc = sp.csr_matrix((10, 10), dtype=np.float64)
cc[0, 0] = 1.

# WRONG
cf = sp.csr_matrix((10, 10), dtype=np.float32)
cf[0, 0] = 1.
print('float32 == ', cf.dtype)
cf += cc
print('float32 == ', cf.dtype)

# WORKS AS EXPECTED
cf = sp.csr_matrix((10, 10), dtype=np.float32)
cf[0, 0] = 1.
print('float32 == ', cf.dtype)
cf[:, :] += cc
print('float32 == ', cf.dtype)

Error message:

There is no error message, it simply upcasts the arrays, also when adding a complex array (say only for real part).

Scipy/Numpy/Python version information:

scipy: 0.19.1
numpy: 1.13.1
Python: 2.7.13 and 3.6.2

[pv]

There's no inplace operation for sparse matrices. `a += b` expands to the default Python operation `a = a + b`.

[zerothi]

So you suggest that it works as expected? I.e. upcasting is an accepted side-effect?

My expectations would be that of numpy (where the above thing works as expected, cf is always np.float32).

[pv]

That there is no in-place operation indeed differs from numpy arrays. This late in the game, with scipy.sparse 10+ years old, I do not really see a way to change this. . The way forward is likely to add new `csr_array` etc classes, which have fully numpy array compatible semantics, and soft-deprecate the `csr_matrix` etc classes.

[matthew-brett]

csr_array etc sounds like an excellent way forwards ...

[zerothi]

I agree. This sounds like a better game-plan.

In my own code I had to implement a replacement for csr_matrix, i.e. dim == 3 where the first 2 are the sparse matrix and the last any number to bypass some of the problems I had with the csr_matrix. However, it has other shortcomings... :(

[denis-bz]

It would be nice if += could do the equivalent of

def inc( X, Y, c=1. ):
""" X += c * Y, X Y sparse or dense
NB not tested on all combinations of: ndarray () (1,) (1,1), numpy matrix, sparse matrix
"""
if (not hasattr( X, "indices" ) # dense += sparse
and hasattr( Y, "indices" )):
# inc an ndarray view, because ndarry += sparse -> matrix --
X = getattr( X, "A", X ).squeeze()
X[Y.indices] += c * Y.data
else:
X += c * Y # sparse + different sparse: SparseEfficiencyWarning
return X

(A matrix of all "X op Y" combinations that work / don't work / upcast
would I think help users and implementers too.)

[zerothi]

I have just refound this bug (time forgets).

I think the proper solution (which may break backward compatibility) would be to raise NotImplementedError when users try to use in-place operations.

I.e. a user may expect:

import numpy as np
from scipy.sparse import csr_matrix, diags
def func(A, b):
    A += b
A = csr_matrix((10, 10))
A[1, 2] = 1
b = diags(np.arange(10))
func(A, b)
print(A)

this will yield:

  (1, 2)	1.0

it should however yield:

 (1, 2)	1.0
  (1, 1)	1.0
  (1, 2)	1.0
  (2, 2)	2.0
  (3, 3)	3.0
  (4, 4)	4.0
  (5, 5)	5.0
  (6, 6)	6.0
  (7, 7)	7.0
  (8, 8)	8.0
  (9, 9)	9.0

The above func will work correctly for numpy arrays.

PS. for others, my current fix is this (which only works for CSC/CSR matrices):

def func(A, b):
    AA = A + b
    # Restore data in the A array
    A.indices = AA.indices
    A.indptr = AA.indptr
    A.data = AA.data
    del AA

[rgommers]

I don't think we want to break backwards compat at this point. At most maybe a warning if a user uses an in-place operation and gets a different dtype back (if this can be checked without costing too much performance).

pydata/sparse#146 adds support for in-place operations to pydata/sparse. I don't see an explicit test for preserving dtype though, and the implementation is similarly a convenience rather than a real in-place operation. So it may suffer from the same issue. @hameerabbasi do you know?

[hameerabbasi]

Hi! There’s no real way that I know of to do in-place addition on Sparse arrays. The reason is that the size of the arrays may change so there’s no real way to do it.

That said, PyData/Sparse does “faux in-place” operations like @pv described, mainly for compatibility with NumPy.

The way it’s laid out currently, we do not preserve the dtype, but preserving that is a few line change at the very most. Since PyData/Sparse prefers NumPy compatibility over all else, doing this is certainly desired for us.

[ev-br]

It's certainly possible to do in-place binops on some formats of sparse arrays. E.g. DOK-type arrays, https://github.com/ev-br/sparr/blob/master/sparr/sp_map.h#L312

[zerothi]

@hameerabbasi I don't see a problem with doing in-place additions since the problem is that the object retaining the data should be the same, see e.g. my func which is basically __iadd__. I see that you cannot ensure the sparsity pattern to be the same, but that is fine.

[hameerabbasi]

Yes, the object remains the same, but all the arrays and memory inside it change. That's why I call them "faux in-place" operations.

As for maintaining the dtype, I've added an issue pydata/sparse#205 and a PR pydata/sparse#206 to fix it and bring behavior in line with NumPy.

[zerothi]

Ok, I see.