在涉及线性代数的功能上应用Numpy广播

Question

我想在数学函数上使用numpy广播功能，该功能涉及线性代数（没有分母部分的双变量高斯分布）。 我代码的最小，可再现的示例是：

我有以下功能

import numpy as np
def  gaussian(x):
    mu = np.array([[2],
                   [2]])
    sigma = np.array([[10, 0],
                      [0, 10]])
    xm = x - mu
    result = np.exp((-1/2) * xm.T @ np.linalg.inv(sigma) @ xm)
    return result

该函数假设x是2x1数组。 我的目的是使用该函数生成一个2D数组，其中各个元素是该功能的产品。 我将此功能应用如下：

x, y = np.arange(5), np.arange(5)
xLen, yLen = len(x), len(y)
z = np.zeros((yLen, xLen))

for y_index in range(yLen):
        for x_index in range(xLen):
            element = np.array([[x[x_index]],
                                [y[y_index]]])
            result = gaussian(element)
            z[y_index][x_index] = result

这有效，但是如您所见，我使用两个用于索引。我知道这是不良的做法，当使用更大的阵列时，它非常慢。我想用numpy广播功能解决此问题。我尝试了以下代码：

X, Y = np.meshgrid(x, y, indexing= 'xy')
element = np.array([[X],
                    [Y]])
Z = gaussian(element)

但是我遇到了此错误：valueerror：操作数无法与形状（2,1,5,5）（2,1,1）（2,1）（2,1） for Line xm = x -mu该功能的。我在一定程度上理解了这个错误。

此外，即使我解决了这个问题，我也会遇到另一个错误：value eRror：matmul：输入操作数1在其核心维度0中具有不匹配，带有gufunc签名（n？，k），（k，m？） - ＆gt;（n？，m？）（尺寸5与2不同） result = np.exp（（ - 1/2） * xm.t @ np @ np.linalg.inv（ sigma） @ xm） fuction的行。再次，我明白为什么。 XM不再是2x1数组，并将其乘以Sigma（即2x2）是行不通的。

有人对如何修改我的功能有建议，以便广播实现有效吗？

Answer 1

以下可能有效。
要注意的两件事：

• 我正在使用np.einsum用于矢量 - 矩阵向量乘法。可能会有更快的方法，但这可以很好地处理要播放的其他维度。
• 根据我的经验，对于更大的阵列，使用具有3个以上维度的广播，情况可能不会比简单的嵌套环更快。我还没有挖掘过：也许计算是在错误的维度上完成的（列 - 与行问题），这会减慢速度。因此，也许通过调整或播放尺寸订单，可以加速

设置代码

nx, ny = 5, 5
x, y = np.arange(nx), np.arange(ny)
X, Y = np.meshgrid(x, y, indexing= 'xy')
element = np.array([[X],
                    [Y]])
# Stack X and Y into a nx x ny x 2 array
XY = np.dstack([X, Y])

新功能

def  gaussian(x):
    # Note that I have removed the extra dimension: 
    # mu is a simple array of shape (2,)
    # This is no problem, since we're using einsum
    # for the matrix multiplication
    mu = np.array([2, 2])
    sigma = np.array([[10, 0],
                      [0, 10]])
    # Broadcast xm to x's shape: (nx, ny, 2)
    xm = x - mu[..., :]
    invsigma = np.linalg.inv(sigma)
    # Compute the (double) matrix multiplication
    # Leave the first two dimension (ab) alone
    # The other dimensions will sum up to a single scalar
    # and thus only the ab dimensions are there in the output
    alpha = np.einsum('abi,abj,ji->ab', xm, xm, invsigma)
    result = np.exp((-1/2) * alpha)
    # The shape of result is (nx, ny)
    return result

，然后致电：

gaussian(XY)

显然，请仔细检查。我做了一项简短的检查，这似乎是正确的，但是转录错误可能已交换了尺寸。

Answer 2

因此，（2,1）输入返回a（1,1）结果：

In [83]: gaussian(np.ones((2,1)))
Out[83]: array([[0.90483742]])

添加一些领先维度：

In [84]: gaussian(np.ones((3,4,2,1)))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [84], in <cell line: 1>()
----> 1 gaussian(np.ones((3,4,2,1)))

Input In [80], in gaussian(x)
      4 sigma = np.array([[10, 0],
      5                   [0, 10]])
      6 xm = x - mu
----> 7 result = np.exp((-1/2) * xm.T @ np.linalg.inv(sigma) @ xm)
      8 return result

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

x-mu有效，因为（3,4,2,1）广播（2,1）

该错误发生在（ - 1/2） * xm.t @ np.linalg.inv（sigma）

np.linalg.inv（sigma） is（2,2） ）

xm是（3,4,2,1），因此其转置为（1,2,4,3）。

相反，如果数组为（3,4,1,2） @（2,2） @（3,4,2,1）结果应为（3,4,1,1）。

因此，让我们完善转言：

def  gaussian(x):
    mu = np.array([[2],
                   [2]])
    sigma = np.array([[10, 0],
                      [0, 10]])
    xm = x - mu
    xmt =xm.swapaxes(-2,-1)
    result = np.exp((-1/2) * xmt @ np.linalg.inv(sigma) @ xm)
    return result

现在它适用于原始（2,1）和任何其他（N，M，2,1）形状：

In [87]: gaussian(np.ones((3,4,2,1))).shape
Out[87]: (3, 4, 1, 1)

In [88]: gaussian(np.ones((2,1))).shape
Out[88]: (1, 1)