Element-wise maximum of two sparse matrices

Question

Is there an easy/build-in way to get the element-wise maximum of two (or ideally more) sparse matrices? I.e. a sparse equivalent of np.maximum.

Answer 1

This did the trick:

def maximum (A, B):
    BisBigger = A-B
    BisBigger.data = np.where(BisBigger.data < 0, 1, 0)
    return A - A.multiply(BisBigger) + B.multiply(BisBigger)

Answer 2

No, there's no built-in way to do this in scipy.sparse. The easy solution is

np.maximum(X.A, Y.A)

but this is obviously going to be very memory-intensive when the matrices have large dimensions and it might crash your machine. A memory-efficient (but by no means fast) solution is

# convert to COO, if necessary
X = X.tocoo()
Y = Y.tocoo()

Xdict = dict(((i, j), v) for i, j, v in zip(X.row, X.col, X.data))
Ydict = dict(((i, j), v) for i, j, v in zip(Y.row, Y.col, Y.data))

keys = list(set(Xdict.iterkeys()).union(Ydict.iterkeys()))

XmaxY = [max(Xdict.get((i, j), 0), Ydict.get((i, j), 0)) for i, j in keys]
XmaxY = coo_matrix((XmaxY, zip(*keys)))

Note that this uses pure Python instead of vectorized idioms. You can try shaving some of the running time off by vectorizing parts of it.

Answer 3

Here's another memory-efficient solution that should be a bit quicker than larsmans'. It's based on finding the set of unique indices for the nonzero elements in the two arrays using code from Jaime's excellent answer here.

import numpy as np
from scipy import sparse

def sparsemax(X, Y):

    # the indices of all non-zero elements in both arrays
    idx = np.hstack((X.nonzero(), Y.nonzero()))

    # find the set of unique non-zero indices
    idx = tuple(unique_rows(idx.T).T)

    # take the element-wise max over only these indices
    X[idx] = np.maximum(X[idx].A, Y[idx].A)

    return X

def unique_rows(a):
    void_type = np.dtype((np.void, a.dtype.itemsize * a.shape[1]))
    b = np.ascontiguousarray(a).view(void_type)
    idx = np.unique(b, return_index=True)[1]
    return a[idx]

Testing:

def setup(n=1000, fmt='csr'):
    return sparse.rand(n, n, format=fmt), sparse.rand(n, n, format=fmt)

X, Y = setup()
Z = sparsemax(X, Y)
print np.all(Z.A == np.maximum(X.A, Y.A))
# True

%%timeit X, Y = setup()
sparsemax(X, Y)
# 100 loops, best of 3: 4.92 ms per loop

Answer 4

The latest scipy (13.0) defines element-wise booleans for sparse matricies. So:

BisBigger = B>A
A - A.multiply(BisBigger) + B.multiply(BisBigger)

np.maximum does not (yet) work because it uses np.where, which is still trying to get the truth value of an array.

Curiously B>A returns a boolean dtype, while B>=A is float64.

Answer 5

Here is a function that returns a sparse matrix that is element-wise maximum of two sparse matrices. It implements the answer by hpaulj:

def sparse_max(A, B):
    """
    Return the element-wise maximum of sparse matrices `A` and `B`.
    """
    AgtB = (A > B).astype(int)
    M = AgtB.multiply(A - B) + B
    return M

Testing:

A = sparse.csr_matrix(np.random.randint(-9,10, 25).reshape((5,5))) 
B = sparse.csr_matrix(np.random.randint(-9,10, 25).reshape((5,5)))

M = sparse_max(A, B)
M2 = sparse_max(B, A)

# Test symmetry:
print((M.A == M2.A).all())
# Test that M is larger or equal to A and B, element-wise:
print((M.A >= A.A).all())
print((M.A >= B.A).all())

Answer 6

from scipy import sparse
from numpy import array
I = array([0,3,1,0])
J = array([0,3,1,2])
V = array([4,5,7,9])
A = sparse.coo_matrix((V,(I,J)),shape=(4,4))

A.data.max()
9

If you haven't already, you should try out ipython, you could have saved your self time my making your spare matrix A then simply typing A. then tab, this will print a list of methods that you can call on A. From this you would see A.data gives you the non-zero entries as an array and hence you just want the maximum of this.

Answer 7

On the current SciPy you may use the object method maximum():

mM = mA.maximum(mB)