Part of the problem here is that you're discarding the mask information by forcing the result into a standard array. Your code is equivalent in result to incr = incr + np.array(x*y).

Fair enough. This should be documented down somewhere tho

This is a problem faced by in-place operations more generally, that there exists no __iradd__ and hence the object being added cannot influence how the operation proceeds. This would be solved once we implement MaskedArray using __array_ufunc__ since in that case it can influence the result (more specifically, it can cause it to fail, which is probably appropriate here: the mask should not just be discarded without the user explicitly doing that). Anyway, agreed it would be good to document... PR most welcome...

This problem came up in scipy sparse as a result of a __array_ufunc__ change, didn't it? What was the resolution there regarding correct implementation of __array_ufunc__?

Edit: that was #9019.

Ah, so that was for a change in behaviour of __array_priority__. With __array_ufunc__, the sparse matrix could just return NotImplemented itself and things should work. Anyway, I think we're getting a bit off-topic here, except that I think we could solve the issue in a future MaskedArray that relied on __array_ufunc__.

But this already did the right thing for scipy.sparse before __array_ufunc__, right? How was that achieved? Is it because it was a duck type?

The way __array_priority__ was handled was changed a bit as part of introducing __array_ufunc__ - since it was attempted to sanitize the code allowing overrides more generally - and this is what tripped up scipy.sparse.