How to ignore implicit zeros using scipy.sparse.csr_matrix.minimum?

Viewed 347

I have two matrices mat1 and mat2 that are sparse (most entries are zero) and I'm not interested in the zero-valued entries: I look at the matrices from a graph-theoretical perspective where a zero means that there is no edge between the nodes.

How can I efficiently get the minimum values between non-zero entries only using scipy.sparse matrices? I.e. an equivalent of mat1.minimum(mat2) that would ignore implicit zeros.

Using dense matrices, it is fairly easy to do:

import numpy as np
nnz = np.where(np.multiply(mat1, mat2))
m = mat1 + mat2
m[nnz] = np.minimum(mat1[nnz], mat2[nnz])

But this would be very inefficient with sparse matrices.

NB: a similar question has been asked before but did not get any relevant answer and there is a related PR on the scipy repo that proposes an implementation of this for (arg)min/max but not for minimum.

EDIT: to specify a bit more the desired behavior would be commutative, i.e. this nonzero-minimum would take all values present in only one of the two matrices and the min of the entries that are present in both matrices

2 Answers

Just in case someone also looks for this, my current implementation is below. However, I'd appreciate any proposal that would either speed this up or reduce the memory footprint.

s = mat1.multiply(mat2)
s.data[:] = 1.

a1 = mat1.copy()
a1.data[:] = 1.
a1 = (a1 - s).maximum(0)

a2 = mat2.copy()
a2.data[:] = 1.
a2 = (a2 - s).maximum(0)

res = mat1.multiply(a1) + mat2.multiply(a2) + \
      mat1.multiply(s).minimum(mat2.multiply(s))

If the sparse nonzeros are positive, an alternate way to use the correct UNION behavior of maximum might be to negate and make positive.

Following your lead of mucking with data explicitly. I found

def sp_min_nz_positive(asp,bsp):  # a and b scipy sparse
    amax = asp.max()
    bmax = bsp.max()
    abmaxplus = max(amax, bmax) # + 1.0 : surprise! not needed.
    # invert the direction, while remaining positive
    arev = asp.copy()
    arev.data[:] = abmaxplus - asp.data[:]
    brev = bsp.copy()
    brev.data[:] = abmaxplus - bsp.data[:]
    out = arev.maximum(brev)  # 
    # revert the direction of these positives
    out.data[:] = abmaxplus - out.data[:]
    return out

there may be inexactness due to roundoff

There was also a suggestion to use sparse internals. A rather generic function is sp.find which returns the nonzero elements of anything.

So you could also try out a minimum that handles negative values too, with something like:

import scipy.sparse as sp
def sp_min_union(a, b):
    assert a.shape == b.shape
    assert sp.issparse(a) and sp.issparse(b)
    (ra,ca,_) = sp.find(a)   # over nonzeros only
    (rb,cb,_) = sp.find(b)   # over nonzeros only
    setab = set(zip(ra,ca)).union(zip(rb,cb)) # row-column union-of-nonzero
    r=[]
    c=[]
    v=[]
    for (rr,cc) in setab:
        r.append(rr)
        c.append(cc)
        anz = a[rr,cc]
        bnz = b[rr,cc]
        assert anz!=0 or bnz!=0   # they came from *some* sp.find
        if anz==0: anz = bnz
        #else:
        #    #if bnz==0: anz = anz
        #    #else:      anz=min(anz,bnz)
        # equiv.
        elif bnz!=0: anz=min(anz,bnz)
        v.append(anz)
    # choose what sparse output format you want, many seem
    # constructible as:
    return sp.csr_matrix((v, (r,c)), shape=a.shape)
Related