Why can't I assign data to part of sparse matrix in the first "try:"?

Question

I want to assign a value to part of a crs sparse matrix (I know it's expensive but it doesn't matter in my project). I tried to assign a float variable to part of the sparse matrix but it doesn't work the first time. However, if I do the exact same thing in the "except" it will work flawlessly.

I then tried to check the dtype of the sparse matrix and part of it and they are different for some reason. The datatype of the whole matrix is float16 as I assigned, but part of the matrix has a float32 dtype.

Here's a small example for both issues:

from scipy.sparse import csr_matrix
import numpy as np

frame = csr_matrix((10, 10),dtype=np.float16)

print "================\n================ Part 1\n================"
print "Let's assign a value to part of the sparse matrix:"
try:
    frame[0:3,0:3] = np.float16(0.6)
    print "The first attempt worked!"
except:
    print "The first attempt didn't work"

print "let's try again :"

try:
    frame[0:3,0:3] = np.float16(0.6)
    print "The second attempt worked!"
except:
    print "The second attempt didn't work"

print "================\n================ Part 2\n================"
print "Let's check the datatype:"
print "Frame dtype is:",; print frame.dtype
print "Part-of-frame dtype is",; print frame[0:3,0:3].dtype

and here's the result:

================
================ Part 1
================
Let's assign a value to part of the sparse matrix:
The first attempt didn't work
let's try again :
The second attempt worked!
================
================ Part 2
================
Let's check the datatype:
Frame dtype is: float16
Part-of-frame dtype is float32

Still I tried to assign a np.float32 to part of the sparse matrix, and I got the same behaviour. Can someone explain what is happening?

Answer 1

Short answer - the problem is with the np.float16; stick with the usually 32 or 64 floats.

===============================

First a working case (in 0.17)

In [334]: M=sparse.csr_matrix((5,5),dtype=np.float)
In [335]: M[:3,:3]=0.6
/usr/lib/python3/dist-packages/scipy/sparse/compressed.py:730: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)
In [336]: M.A
Out[336]: 
array([[ 0.6,  0.6,  0.6,  0. ,  0. ],
       [ 0.6,  0.6,  0.6,  0. ,  0. ],
       [ 0.6,  0.6,  0.6,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ]])
In [338]: M.data
Out[338]: array([ 0.6,  0.6,  0.6,  0.6,  0.6,  0.6,  0.6,  0.6,  0.6])

If I assign again I don't get the sparsity warning.

Now if I use your dtype, I get the sparsity warning, but also a ValueError (your generic except hid that):

In [339]: M=sparse.csr_matrix((5,5),dtype=np.float16)
In [340]: M[:3,:3]=np.float16(0.6)
/usr/lib/python3/dist-packages/scipy/sparse/compressed.py:730: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-340-aaecba748069> in <module>()
----> 1 M[:3,:3]=np.float16(0.6)

/usr/lib/python3/dist-packages/scipy/sparse/compressed.py in __setitem__(self, index, x)
    654             return
    655         i, j = self._swap((i.ravel(), j.ravel()))
--> 656         self._set_many(i, j, x.ravel())
    657 
    658     def _setdiag(self, values, k):

/usr/lib/python3/dist-packages/scipy/sparse/compressed.py in _set_many(self, i, j, x)
    738             j = j[mask]
    739             j[j < 0] += N
--> 740             self._insert_many(i, j, x[mask])
    741 
    742     def _insert_many(self, i, j, x):

/usr/lib/python3/dist-packages/scipy/sparse/compressed.py in _insert_many(self, i, j, x)
    805             # TODO: only sort where necessary
    806             self.has_sorted_indices = False
--> 807             self.sort_indices()
    808 
    809         self.check_format(full_check=False)

/usr/lib/python3/dist-packages/scipy/sparse/compressed.py in sort_indices(self)
   1039         if not self.has_sorted_indices:
   1040             fn = _sparsetools.csr_sort_indices
-> 1041             fn(len(self.indptr) - 1, self.indptr, self.indices, self.data)
   1042             self.has_sorted_indices = True
   1043 

ValueError: Output dtype not compatible with inputs.

I get the same error if I just attempt to convert the csr to dense:

In [346]: M=sparse.csr_matrix((5,5),dtype=np.float16)
In [347]: M.A
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-347-bdd665fbe1b0> in <module>()
----> 1 M.A

/usr/lib/python3/dist-packages/scipy/sparse/base.py in __getattr__(self, attr)
    511     def __getattr__(self, attr):
    512         if attr == 'A':
--> 513             return self.toarray()
    514         elif attr == 'T':
    515             return self.transpose()

/usr/lib/python3/dist-packages/scipy/sparse/compressed.py in toarray(self, order, out)
    938     def toarray(self, order=None, out=None):
    939         """See the docstring for `spmatrix.toarray`."""
--> 940         return self.tocoo(copy=False).toarray(order=order, out=out)
    941 
    942     ...
ValueError: Output dtype not compatible with inputs.

So the problem isn't with the assignment, but with the handling of that dtype. It looks like a bug, but you can avoid it by not using this dtype.

And lil matrix also has problems with this dtype:

In [348]: M=sparse.lil_matrix((5,5),dtype=np.float16)
In [349]: M[:3,:3]=np.float16(0.6)
...
KeyError: (dtype('int32'), dtype('float16'))

At some point in the conversions, the specified dtype isn't being preserved.

But why use float16 in the first place? float32 and float64 are the normal numpy floats.

Even when I successfully create a np.float16 matrix, that dtype is not preserved during most sparse operations:

In [374]: M=sparse.csr_matrix(np.arange(9).reshape(3,3), dtype=np.float16)
In [375]: M.data
Out[375]: array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.], dtype=float16)
In [376]: 
In [376]: M
Out[376]: 
<3x3 sparse matrix of type '<class 'numpy.float16'>'
    with 8 stored elements in Compressed Sparse Row format>

In [377]: M.A   # same error converting to dense
...
ValueError: Output dtype not compatible with inputs.

In [378]: M.T     # dtype kept during transpose
Out[378]: 
<3x3 sparse matrix of type '<class 'numpy.float16'>'
    with 8 stored elements in Compressed Sparse Column format>

dtype is not preserved with multiplication or indexing. Most likely the compiled matrix multiplication code is written for 32 and 64 (regular C floats and doubles), and not 16.

In [379]: M*M
Out[379]: 
<3x3 sparse matrix of type '<class 'numpy.float32'>'
    with 9 stored elements in Compressed Sparse Row format>
In [380]: M[0,:]
Out[380]: 
<1x3 sparse matrix of type '<class 'numpy.float32'>'
    with 2 stored elements in Compressed Sparse Row format>

I don't know if the documentation warns about float16, but I think it's nearly useless.