Does numpy provide a generalized inner product?

Question

Most array-oriented languages such as APL or J have some form of generalized inner product, which can act like standard matrix multiplication, but supports arbitrary operations in place of the standard ones. For example, in J +/ . * is the standard multiply-then-sum, but you can also do e.g. <./ . + to get an add-then-min operation (say for incrementally updating the length of the shortest path through a graph).

In slow and 2D-only Python, this would be something like:

import numpy as np

def general_inner(f, g, x, y):
    return np.array([[f(g(x1, y1)) for y1 in y.T] 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 slicing. We can reshape the two arguments so that an operation will get broadcast instead of performed elementwise, and then perform a reducing operation along the 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 hesitate to call it a "generalized inner product", because inner products have a specific mathematical structure that this new version lacks.

The way this basically works is that any numpy.newaxis has a length of 1 and gets broadcast, so:

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

Gives us:

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

Or if it helps you understand (I find broadcasting to be a bit confusing at times),

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,y.T), utilizing broadcasting, 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,y.T)
Out[137]: 
array([[ 0.9,  0.2],
       [ 2.7,  0.4]])
In [138]: np.sum(np.multiply(x,y.T),axis=1)
Out[138]: array([ 1.1,  3.1])

similarly for the maximum of the sums:

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

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

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

There was an old enhancement request to implement this sort of thing within np.einsum. That implements a sum of products calculation with a high degree of control over the indexing. So conceptually it could perform max of sums with the same indexing control. But as far as I know no one has tried to implement it.

That generalized inner product was a cute feature of APL (it was my primary language several decades ago). But apparently not so useful that it migrated out of that family of languages. MATLAB doesn't have anything like it.

Is there anything that APL & J can do with this construct, that can't be done in numpy with kind of broadcasting that we have demonstrated?

---

With more general x and y shapes, I need to add the newaxis as given in the other 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, while using the matmul idea of last dim/2nd last of 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

So while numpy (and MATLAB) doesn't have a special syntax like APL, the idea of an expansion (outter) operation followed by a reduction is common enough.

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)