Populate numpy matrix from the difference of two vectors

Question

Is it possible to construct a numpy matrix from a function? In this case specifically the function is the absolute difference of two vectors: S[i,j] = abs(A[i] - B[j]). A minimal working example that uses regular python:

import numpy as np

A = np.array([1,3,6])
B = np.array([2,4,6])
S = np.zeros((3,3))

for i,x in enumerate(A):
    for j,y in enumerate(B):
        S[i,j] = abs(x-y)

Giving:

[[ 1.  3.  5.]
 [ 1.  1.  3.]
 [ 4.  2.  0.]]

It would be nice to have a construction that looks something like:

def build_matrix(shape, input_function, *args)

where I can pass an input function with it's arguments and retain the speed advantage of numpy.

Answer 1

In addition to what @JoshAdel has suggested, you can also use the outer method of any numpy ufunc to do the broadcasting in the case of two arrays.

In this case, you just want np.subtract.outer(A, B) (Or, rather, the absolute value of it).

While either one is fairly readable for this example, in some cases broadcasting is more useful, while in others using ufunc methods is cleaner.

Either way, it's useful to know both tricks.

E.g.

import numpy as np

A = np.array([1,3,6])
B = np.array([2,4,6])

diff = np.subtract.outer(A, B)
result = np.abs(diff)

Basically, you can use outer, accumulate, reduce, and reduceat with any numpy ufunc such as subtract, multiply, divide, or even things like logical_and, etc.

For example, np.cumsum is equivalent to np.add.accumulate. This means you could implement something like a cumdiv by np.divide.accumulate if you even needed to.

Answer 2

I recommend taking a look into numpy's broadcasting capabilities:

In [6]: np.abs(A[:,np.newaxis] - B)
Out[6]: 
array([[1, 3, 5],
       [1, 1, 3],
       [4, 2, 0]])

http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html

Then you could simply write your function as:

In [7]: def build_matrix(func,args):
   ...:     return func(*args)
   ...: 

In [8]: def f1(A,B):
   ...:     return np.abs(A[:,np.newaxis] - B)
   ...: 

In [9]: build_matrix(f1,(A,B))
Out[9]: 
array([[1, 3, 5],
       [1, 1, 3],
       [4, 2, 0]])

This should also be considerably faster than your solution for larger arrays.