`numpy.ufunc.outer` for only certain axes

Viewed 59

I have two numpy 1D arrays of points in space and I want a 2D array of the pairwise difference between the points from each array. It just so happens that if the space were 1D then what I would want is equivalent to

x,y = np.array([1,2,3,4]),np.array([3,4,5,6])
np.difference.outer(x,y)

How do I do this for difference of Nd spaces?

2 Answers
In [118]: np.subtract.outer(x,y)
Out[118]: 
array([[-2, -3, -4, -5],
       [-1, -2, -3, -4],
       [ 0, -1, -2, -3],
       [ 1,  0, -1, -2]])

The outer action can be duplicated with broadcasting:

In [120]: np.subtract(x[:,None],y)    # x[:,None]-y
Out[120]: 
array([[-2, -3, -4, -5],
       [-1, -2, -3, -4],
       [ 0, -1, -2, -3],
       [ 1,  0, -1, -2]])

That can be generalized to more dimensions. just so long as you obey the broadcasting rules.

Because of the power of broadcasting, I haven't tried to work out what the ufunc.outer does in multidimensional array. With a pair of 1d arrays the action is obvious, and simply illustrated in the docs.

Here's a good illustration of how outer combines axes:

In [132]: np.subtract.outer(np.ones((2,3)), np.ones((4,5))).shape
Out[132]: (2, 3, 4, 5)

it doesn't accept an axis parameter, so there's no way to refine the action. Overall it's not that useful of a method.

Here's a way using 3D matrices that works in however many dimensions.

import numpy as np

N = 4 # number of vectors
p = 3 # vector dimensionality
x = np.arange(0, p * N).reshape((N, p))
y = np.arange(2, 2 + p * N).reshape((N, p))

Y = np.zeros((N, N, p))
Y[:] = y # setting rows
X = np.zeros((N, N, p))
X[:] = x
X = np.einsum("ijk->jik", X) # swap rows and columns
diff = np.subtract(X, Y) # reproduces 1D result with p = 0 and print(diff[:,:,0])

There's likely an easier way with indexing to set the columns of the X matrix than setting the rows and then doing the specific transpose.

Related