Perform cross-difference in numpy

Question

I would like to perform a 'cross-difference' as follows:

cross_diff ( [a,b], [c,d] ) = [ [ a - c, a - d], [ b - c, b - d] ]

I have a routine to do this in python as follows:

def crossdiff(a,b):
    c = []
    for a1 in range(len(a)):
        for b1 in range(len(b)):
            c.append (a[a1]-b[b1])
    x = numpy.array(c)
    x.reshape(len(a),len(b))
    return x

The problem is that I have to create a python array and stuff the results into it, then convert back to a numpy array. I would like to be able to take numpy vectors a and b and get a numpy array c which contains all the differences, as the performance of the code above is poor for large vector sizes.

Is there away to perform the above calculation as "pure" numpy operations?

EDIT FOR TESTING RESULTS:

I ran all four implementations listed in this thread through the Python profiler for comparison. I had to run them on a workstation since the initial implementation uses ~4GB of RAM with 10k elements.

import numpy
import cProfile

def cross_diff(A, B):
    return A[:,None] - B[None,:]

def crossdiff2 (a,b):
    ap = numpy.tile (a, (numpy.shape(b)[0],1))
    bp = numpy.tile (b, (numpy.shape(a)[0],1))

    return ap - bp.transpose()

def crossdiff(a,b):
    c = []
    for a1 in range(len(a)):
        for b1 in range(len(b)):
            c.append (a[a1]-b[b1])
    x = numpy.array(c)
    x.reshape(len(a),len(b))
    return x

a = numpy.array(range(10000))
b = numpy.array(range(10000))

cProfile.run('crossdiff (a,b)')
cProfile.run('crossdiff2 (a,b)')
cProfile.run('cross_diff (a,b)')
cProfile.run('numpy.subtract.outer (a,b)')

Results: Original python is 74.147 seconds, my version is 1.656 seconds, 3rd implementation 0.296 and 4th 0.288.

Answer 1

You can get the same result without the explicit dimension adding for broadcasting by using the .outer method of ufuncs. For example:

>>> np.subtract.outer([1, 2], [3, 4, 5])
array([[-2, -3, -4],
       [-1, -2, -3]])

It has the added benefit of accepting any iterable as input, you don't need to first convert them to arrays.