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 comments 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 ~\anaconda3\lib\site-packages\scipy\sparse\base.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)

In the linear algebra world (esp. finite difference and finite elements) where sparse math was developed, matrix multiplication was important. Other math that operates just on the nonzero elements is fairly easy. But operations that change sparsity are (relatively) expensive. New values, and submatrices are best added to the coo inputs. coo duplicates are added when converted to csr. So addition/subtraction of whole matrices is (intentionally) limited.