numpy 是否提供广义内积？

Question

大多数面向数组的语言（例如 APL 或 J）都有某种形式的广义内积，它可以像标准矩阵乘法一样工作，但支持任意操作来代替标准操作。例如，在 J +/ . * 中是标准的乘然后和，但您也可以这样做，例如<./ . + 得到一个 add-then-min 操作（比如增量更新通过图的最短路径的长度）。

在慢速且仅 2D 的 Python 中，这将类似于：

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))

numpy 是否提供任何东西来直接支持这种模式？

Answer 1

你可以通过切片到达那里。我们可以重塑这两个参数，以便一个操作将被广播而不是按元素执行，然后沿着不需要的轴执行减少操作。

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))

我不愿称其为“广义内积”，因为内积具有新版本所缺乏的特定数学结构。

这基本上工作的方式是任何numpy.newaxis 的长度为 1 并被广播，所以：

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

给我们：

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

或者如果它可以帮助您理解（我发现广播有时有点令人困惑），

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

Answer 2

不那么慢的numpy 等价物是g(x,y.T)，利用广播，然后是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])

类似地求和的最大值：

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])

np.add、np.multiply、np.maximum 是 ufunc，而 np.sum、np.prod、np.max 不是，但使用 axis 参数和 @987654334 @。 （编辑：np.add.reduce 是 ufunc 等价于 np.sum。）

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

在np.einsum 中有一个旧的增强请求来实现这种事情。这实现了sum of products 计算，并对索引进行了高度控制。所以从概念上讲，它可以使用相同的索引控制执行max of sums。但据我所知，没有人尝试过实现它。

广义内积是 APL 的一个可爱特性（几十年前它是我的主要语言）。但显然它没有那么有用，以至于它从那个语言家族中迁移出来了。 MATLAB 没有类似的东西。

APL & J 可以用这个结构做些什么，而 numpy 不能用我们展示的那种广播来做？

---

对于更一般的x 和y 形状，我需要添加其他答案中给出的newaxis

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]])

推广到 3d，同时使用 matmul 的最后一个暗淡/第二个最后一个 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

因此，虽然 numpy（和 MATLAB）没有像 APL 这样的特殊语法，但扩展（外部）操作后跟缩减的想法很常见。

测试其他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)