A special type of row by row multiplication by two sparse matrices in Python

What I'm looking for: A way to implement in Python a special multiplication operation for matrices that end up in scipy sparse ( csr ). This is a special kind of multiplication, not matrix multiplication, but Kronecker multiplication and Hadamard aka pointwise multiplication and doesn't seem to have built-in support in scipy.sparse.

Required operation: Each line of the output must contain the results of each product of the elements of the corresponding lines in the two input matrices. Therefore, starting with two equally sized matrices, each with dimensions m by n, the result should have dimensions m by n ^ 2.

It looks like this:

graphic of desired multiplication operation

Python code:

import scipy.sparse
A = scipy.sparse.csr_matrix(np.array([[1,2],[3,4]]))
B = scipy.sparse.csr_matrix(np.array([[0,5],[6,7]]))

# C is the desired product of A and B. It should look like:
C = scipy.sparse.csr_matrix(np.array([[0,5,0,10],[18,21,24,28]]))


What would be a good or efficient way to do this? I've tried looking here on stackoverflow as elsewhere until no luck. So far, I've found it best to do a row-by-row operation in a for loop, but that sounds awful since my input matrices have several million rows and several thousand columns, mostly sparse .


source to share

2 answers

In your example C

, this is the first and last linekron

In [4]: A=np.array([[1,2],[3,4]])
In [5]: B=np.array([[0,5],[6,7]])
In [6]: np.kron(A,B)
array([[ 0,  5,  0, 10],
       [ 6,  7, 12, 14],
       [ 0, 15,  0, 20],
       [18, 21, 24, 28]])
In [7]: np.kron(A,B)[[0,3],:]
array([[ 0,  5,  0, 10],
       [18, 21, 24, 28]])



contains the same values ​​as np.outer

but in a different order.

For large dense arrays einsum

can provide good speed:




does the same as np.kron


As = sparse.csr_matrix(A); Bs ...



is written in Python, so you can probably change it if it does unnecessary calculations.

An iterative solution looks like this:

sparse.vstack([sparse.kron(a,b) for a,b in zip(As,Bs)]).A


Being iterative, I don't expect it to be faster than parrying full kron

. But without going into the logic sparse.kron

, perhaps this is the best I can do.


uses bmat

, so the calculation is:

sparse.bmat([[sparse.kron(a,b)] for a,b in zip(As,Bs)])


But bmat

quite complex, so simplifying it won't be easy.

The solution np.einsum

cannot be easily extended to sparse - does not exist sparse.einsum

, and the intermediate is 3d, which is sparse not processed.


uses a format coo

that is not suitable for working with strings. But working in the spirit of this function, I developed a function that iterates over rows of format matrices csr

. Like kron

and bmat

, I build arrays data

, row

, col

and build a coo_matrix

one of them. This, in turn, can be converted to other formats.

def test_iter(A, B):
    m,n1 = A.shape
    n2 = B.shape[1]
    Cshape = (m, n1*n2)
    data = np.empty((m,),dtype=object)
    col =  np.empty((m,),dtype=object)
    row =  np.empty((m,),dtype=object)
    for i,(a,b) in enumerate(zip(A, B)):
        data[i] = np.outer(a.data, b.data).flatten()
        #col1 = a.indices * np.arange(1,a.nnz+1) # wrong when a isn't dense
        col1 = a.indices * n2   # correction
        col[i] = (col1[:,None]+b.indices).flatten()
        row[i] = np.full((a.nnz*b.nnz,), i)
    data = np.concatenate(data)
    col = np.concatenate(col)
    row = np.concatenate(row)
    return sparse.coo_matrix((data,(row,col)),shape=Cshape)


With these small 2x2 matrices as well as larger ones (e.g. A1=sparse.rand(1000,2000).tocsr()

), this is about 3x faster than the version with bmat

. For large enough matrices, the dense version einsum

(which may have memory errors) is better .



The non-optimal way to do this is to use kron separately for each line:

def my_mult(A, B):
    nrows = A.shape[0]
    prodrows = []
    for i in xrange(0, nrows):
        Arow = A.getrow(i)
        Brow = B.getrow(i)
        prodrow = scipy.sparse.kron(Arow,Brow)
    return scipy.sparse.vstack(prodrows)


This is about 3x worse performance than @hpaulj's solution here , as you can see from the following code:

A=scipy.sparse.rand(20000,1000, density=0.05).tocsr()
B=scipy.sparse.rand(20000,1000, density=0.05).tocsr()

# Check memory
%memit C1 = test_iter(A,B)
%memit C2 = my_mult(A,B)

# Check time
%timeit C1 = test_iter(A,B)
%timeit C2 = my_mult(A,B)

# Last but not least, check correctness!
print (C1 - C2).nnz == 0



Hpaulj method:

peak memory: 1993.93 MiB, increment: 1883.80 MiB
1 loops, best of 3: 6.42 s per loop


this method:

peak memory: 2456.75 MiB, increment: 1558.78 MiB
1 loops, best of 3: 18.9 s per loop




All Articles