Find two pairs of pairs that sum to the same value

algorithm, numpy, performance, python

Solution

Here is a 'lazy' approach, that scales up to n=10000, using 'only' 4gb of memory, and completing in 10s per call or so. Worst case complexity is O(n^3), but for random data, expected performance is O(n^2). At first sight, it seems like youd need O(n^3) ops; each row combination needs to be produced and inspected at least once. But we need not look at the entire row. Rather, we can perform an early exit strategy on the comparison of rowpairs, once it is clear they are of no use to us; and for random data, we may draw this conclusion typically long before we have considered all columns in a row.

import numpy as np

n = 10
#also works for non-square A
A = np.random.randint(2, size=(n*2,n)).astype(np.int8)
#force the inclusion of some hits, to keep our algorithm on its toes
##A[0] = A[1]


def base_pack_lazy(a, base, dtype=np.uint64):
    """
    pack the last axis of an array as minimal base representation
    lazily yields packed columns of the original matrix
    """
    a = np.ascontiguousarray( np.rollaxis(a, -1))
    init = np.zeros(a.shape[1:], dtype)
    packing = int(np.dtype(dtype).itemsize * 8 / (float(base) / 2))
    for columns in np.array_split(a, (len(a)-1)//packing+1):
        yield reduce(
            lambda acc,inc: acc*base+inc,
            columns,
            init)

def unique_count(a):
    """returns counts of unique elements"""
    unique, inverse = np.unique(a, return_inverse=True)
    count = np.zeros(len(unique), np.int)
    np.add.at(count, inverse, 1)        #note; this scatter operation requires numpy 1.8; use a sparse matrix otherwise!
    return unique, count, inverse

def has_identical_row_sums_lazy(A, combinations_index):
    """
    compute the existence of combinations of rows summing to the same vector,
    given an nxm matrix A and an index matrix specifying all combinations

    naively, we need to compute the sum of each row combination at least once, giving n^3 computations
    however, this isnt strictly required; we can lazily consider the columns, giving an early exit opportunity
    all nicely vectorized of course
    """

    multiplicity, combinations = combinations_index.shape
    #list of indices into combinations_index, denoting possibly interacting combinations
    active_combinations = np.arange(combinations, dtype=np.uint32)

    for packed_column in base_pack_lazy(A, base=multiplicity+1):       #loop over packed cols
        #compute rowsums only for a fixed number of columns at a time.
        #this is O(n^2) rather than O(n^3), and after considering the first column,
        #we can typically already exclude almost all rowpairs
        partial_rowsums = sum(packed_column[I[active_combinations]] for I in combinations_index)
        #find duplicates in this column
        unique, count, inverse = unique_count(partial_rowsums)
        #prune those pairs which we can exclude as having different sums, based on columns inspected thus far
        active_combinations = active_combinations[count[inverse] > 1]
        #early exit; no pairs
        if len(active_combinations)==0:
            return False
    return True

def has_identical_triple_row_sums(A):
    n = len(A)
    idx = np.array( [(i,j,k)
        for i in xrange(n)
            for j in xrange(n)
                for k in xrange(n)
                    if i<j and j<k], dtype=np.uint16)
    idx = np.ascontiguousarray( idx.T)
    return has_identical_row_sums_lazy(A, idx)

def has_identical_double_row_sums(A):
    n = len(A)
    idx = np.array(np.tril_indices(n,-1), dtype=np.int32)
    return has_identical_row_sums_lazy(A, idx)


from time import clock
t = clock()
for i in xrange(10):
    print has_identical_double_row_sums(A)
    print has_identical_triple_row_sums(A)
print clock()-t

Extended to include the calculation over sums of triplets of rows, as you asked above. For n=100, this still takes only about 0.2s

Edit: some cleanup; edit2: some more cleanup

Problem

Original source

Related problems

• how to do scatter/gather operations in numpy