Why does `scipy.sparse.csr_matrix` broadcast multiplication but not subtraction?

Question

I am trying to understand solutions to this question here, and while I can just reuse the code I would prefer to know what is happening before I do.

The question is about how to tile a scipy.sparse.csr_matrix object, and the top answer (by @user3357359) at the time of writing shows how to tile a single row of a matrix across multiple rows as:

from scipy.sparse import csr_matrix
sparse_row = csr_matrix([[0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0]])
repeat_number = 3
repeated_row_matrix = csr_matrix(np.ones([repeat_number,1])) * sparse_row

(I have added the sparse_row and repeat_number initialisation to help make things concrete).

If I now convert this to a dense matrix and print as so:

print(f"repeated_row_matrix.todense() = {repeated_row_matrix.todense()}")

This gives output:

repeated_row_matrix.todense() =
[[0 0 0 0 0 1 0 1 1 0 0 0]
 [0 0 0 0 0 1 0 1 1 0 0 0]
 [0 0 0 0 0 1 0 1 1 0 0 0]]

The operation on the right of the repeated_row_matrix assignment seems to me to be performing broadcasting. The original sparse_row has shape (1,12), the temporary matrix is a (3,1) matrix of ones, and the result is a (3,12) matrix. So far, this is similar behaviour as you would expect from numpy.array. However, if I try the same thing with the subtraction operator:

sparse_row = csr_matrix([[0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0]])
repeat_number = 3
repeated_row_matrix = csr_matrix(np.ones([repeat_number,1])) - sparse_row
print(f"repeated_row_matrix.todense() =\n{repeated_row_matrix.todense()}")

I get an error in the third line:

3 repeated_row_matrix = csr_matrix(np.ones([repeat_number,1])) - sparse_row
...
ValueError: inconsistent shapes

Is this intended behaviour? And if so, why?

I guess that a multiplication between two sparse K-vectors with n1 and n2 non-zeros respectively, would always have less than or equal to min(n1,n2) non-zeros. A subtraction would have in the worst case n1+n2 non-zeros but does this really explain why one behaviour is allowed and one is not.

I wish to perform subtraction of a single row vector from a matrix (for a sparse implementation of K-medoids I am playing with). To perform subtraction, I am creating a temporary sparse array which tiles the original row by using broadcasting with multiplication then I can subtract one array from another. I am sure there should be a better way, but I don't see it.

Also, @"C.J. Jackson" replies in the cmets that a better way to construct the tiling is:

sparse_row[np.zeros(repeat_number),:]

This works, but I have no idea why or what functionality is being employed. Can someone point me to the documentation? If sparse_row were a numpy.array then this does not cause tiling.

Thanks in advance.

Answer 1

With dense arrays, broadcasted multiplication and matrix multiplication can do the same thing for special cases. For example with 2 1d arrays

In [3]: x = np.arange(3); y = np.arange(5)

broadcasted:

In [4]: x[:,None]*y   # (3,1)*(5,) => (3,1)*(1,5) => (3,5)
Out[4]: 
array([[0, 0, 0, 0, 0],
       [0, 1, 2, 3, 4],
       [0, 2, 4, 6, 8]])

dot/matrix multiplication of a (3,1) and (1,5). This is not broadcasting. It is doing sum-of-products on the shared size 1 dimension:

In [5]: x[:,None]@y[None,:]
Out[5]: 
array([[0, 0, 0, 0, 0],
       [0, 1, 2, 3, 4],
       [0, 2, 4, 6, 8]])

Make sparse matrices for these:

In [6]: Mx = sparse.csr_matrix(x);My = sparse.csr_matrix(y)    
In [11]: Mx
Out[11]: 
<1x3 sparse matrix of type '<class 'numpy.intc'>'
    with 2 stored elements in Compressed Sparse Row format>    
In [12]: My
Out[12]: 
<1x5 sparse matrix of type '<class 'numpy.intc'>'
    with 4 stored elements in Compressed Sparse Row format>

Note the shapes (1,3) and (1,5). To do the matrix multiplication, the first needs to be transposed to (3,1):

In [13]: Mx.T@My
Out[13]: 
<3x5 sparse matrix of type '<class 'numpy.intc'>'
    with 8 stored elements in Compressed Sparse Column format>

In [14]: _.A
Out[14]: 
array([[0, 0, 0, 0, 0],
       [0, 1, 2, 3, 4],
       [0, 2, 4, 6, 8]], dtype=int32)

Mx.T*My works the same way, because sparse is modeled on np.matrix (and MATLAB), where * is matrix multiplication.

Element-wise multiplication works in the same way as for dense:

In [20]: Mx.T.multiply(My)
Out[20]: 
<3x5 sparse matrix of type '<class 'numpy.intc'>'
    with 8 stored elements in Compressed Sparse Column format>

I'm a little surprised, it does look at bit like broadcasting, though it doesn't involve any automatic None dimensions (sparse is always 2d). Funny, I can't find a element-wise multiplication for the dense matix.

But as you found Mx.T-My raises the inconsistent shapes error. The sparse developers chose not to implement this kind of subtraction (or addition). In general addition or subtraction of sparse matrices is a problem. It can easily result in a dense matrix, if you add something to all elements, including the "implied" 0s.

In [41]: Mx+1
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Input In [41], in <cell line: 1>()
----> 1 Mx+1

File ~naconda3libsite-packagesscipysparsease.py:410, in spmatrix.__add__(self, other)
    408         return self.copy()
    409     # Now we would add this scalar to every element.
--> 410     raise NotImplementedError('adding a nonzero scalar to a '
    411                               'sparse matrix is not supported')
    412 elif isspmatrix(other):
    413     if other.shape != self.shape:

NotImplementedError: adding a nonzero scalar to a sparse matrix is not supported

To replicate the broadcasted subtraction:

In [54]: x[:,None]-y
Out[54]: 
array([[ 0, -1, -2, -3, -4],
       [ 1,  0, -1, -2, -3],
       [ 2,  1,  0, -1, -2]])

We have to 'tile' the matrices. Your link shows some options (including my answer). Another option is to vstack several instances of the matrices. sparse.vstack actually makes a new matrix, using the coo matrix format:

In [55]: Mxx = sparse.vstack([Mx]*5);Myy = sparse.vstack([My,My,My])    
In [56]: Mxx,Myy
Out[56]: 
(<5x3 sparse matrix of type '<class 'numpy.intc'>'
    with 10 stored elements in Compressed Sparse Row format>,
 <3x5 sparse matrix of type '<class 'numpy.intc'>'
    with 12 stored elements in Compressed Sparse Row format>)

Now two (3,5) matrices can be added or subtracted:

In [57]: Mxx.T-Myy
Out[57]: 
<3x5 sparse matrix of type '<class 'numpy.intc'>'
    with 12 stored elements in Compressed Sparse Column format>

In [58]: _.A
Out[58]: 
array([[ 0, -1, -2, -3, -4],
       [ 1,  0, -1, -2, -3],
       [ 2,  1,  0, -1, -2]], dtype=int32)