Does numpy have a common generic product?

Question

Most array-oriented languages, such as APL or J, have some kind of generalized generalized product that can act as a standard matrix multiplication, but supports arbitrary operations instead of standard ones. For example, in J +/ . * +/ . * There is a standard multiplied amount, but you can also do this, for example. <./ . + <./ . + to get the add-then-min operation (say, to gradually update the length of the shortest path through the graph).

In slow and two-dimensional Python, this will be something like:

 import numpy as np def general_inner(f, g, x, y): return np.array([[f(g(x1, y1)) for y1 in yT] for x1 in x]) x = np.arange(1, 5, dtype="float").reshape((2, 2)) y = np.array([[0.9], [0.1]]) assert(x.dot(y) == general_inner(np.sum, np.multiply, x, y)) 

Does numpy provide anything to directly support this pattern?

Answer 1

You can get there with a cut. We can change the two arguments so that the operation receives a broadcast, rather than being performed differently, and then performs a decrease operation along an unwanted axis.

 import numpy a = numpy.array([[1, 2, 3], [4, 5, 6]]) b = numpy.array([[7, 8], [9, 10], [11, 12]]) # Ordinary matrix multiplication print(a @ b) # Matrix multiplication using broadcasting print(numpy.sum(a[:,:,numpy.newaxis] * b[numpy.newaxis,:,:], axis=1)) # Our "generalized" version print(numpy.min(a[:,:,numpy.newaxis] + b[numpy.newaxis,:,:], axis=1)) 

I would not call it a “generalized internal product” because internal products have a certain mathematical structure that this new version does not have.

The principle of operation basically consists in the fact that any numpy.newaxis has a length of 1 and receives broadcast transmission, therefore:

 a[:,:,numpy.newaxis] * b[numpy.newaxis,:,:] 

Gives us:

 result[i,j,k] = a[i,j] * b[j,k] 

Or, if that helps you understand (I find that broadcasting is sometimes a bit confusing),

 aa = a[:,:,numpy.newaxis] bb = b[numpy.newaxis,:,:] result[i,j,k] = aa[i,j,0] * bb[0,j,k] 

Answer 2

The not so slow numpy equivalent is g(x,yT) , using broadcast, followed by f(..., axis=1) .

 In [136]: general_inner(np.sum, np.multiply, x, y) Out[136]: array([[ 1.1], [ 3.1]]) In [137]: np.multiply(x,yT) Out[137]: array([[ 0.9, 0.2], [ 2.7, 0.4]]) In [138]: np.sum(np.multiply(x,yT),axis=1) Out[138]: array([ 1.1, 3.1]) 

similarly for maximum amounts:

 In [145]: general_inner(np.max, np.add, x, y) Out[145]: array([[ 2.1], [ 4.1]]) In [146]: np.max(np.add(x,yT), axis=1) Out[146]: array([ 2.1, 4.1]) 

There is a potential confusion that np.add , np.multiply , np.maximum are ufunc , and np.sum , np.prod , np.max are not, but accept the axis and keepdims . (edit: np.add.reduce is the equivalent of ufunc np.sum .)

 In [152]: np.max(np.add(x,yT), axis=1, keepdims=True) Out[152]: array([[ 2.1], [ 4.1]]) 

There was an old improvement request for implementing this kind of thing at np.einsum . This implements sum of products calculation with a high degree of index control. Thus, conceptually, he could perform max of sums with the same indexing control. But as far as I know, no one tried to implement it.

This generic internal product was a nice feature of APL (it was my main language decades ago). But, apparently, it is not so useful that he migrated from this family of languages. MATLAB has nothing of the kind.

Is there anything APL and J can do with this construct that cannot be done in numpy with the kind of broadcast we demonstrated?

---

With more general forms of x and y I need to add newaxis as indicated in another answer

 In [176]: x = np.arange(3*4).reshape(4,3) In [177]: y = np.arange(3*2).reshape(3,2) In [178]: np.sum(np.multiply(x[...,None],y[None,...]),axis=1) Out[178]: array([[10, 13], [28, 40], [46, 67], [64, 94]]) In [179]: np.max(np.add(x[...,None],y[None,...]),axis=1) Out[179]: array([[ 6, 7], [ 9, 10], [12, 13], [15, 16]]) 

generalizing to 3d, using the idea of matmul last dull / second last from matmul :

 In [195]: x = np.arange(2*4*5).reshape(2,4,5) In [196]: y = np.arange(2*5*3).reshape(2,5,3) In [197]: np.einsum('ijk,ikm->ijm', x, y).shape Out[197]: (2, 4, 3) In [203]: np.add.reduce(np.multiply(x[...,None], y[...,None,:,:]), axis=-2).shape Out[203]: (2, 4, 3) # shapes broadcast: (2,4,5,n) * (2,n,5,3) => (2,4,5,3); sum on the 5 

Thus, although numpy (and MATLAB) does not have a special syntax, such as APL , the idea of ​​an extension operation (outter) followed by a reduction is quite common.

testing other ufunc :

 In [205]: np.maximum.reduce(np.add(x[...,None], y[...,None,:,:]), axis=-2).shape Out[205]: (2, 4, 3) In [208]: np.logical_or.reduce(np.greater(x[...,None], y[...,None,:,:]), axis=-2).shape Out[208]: (2, 4, 3)