根据位置将函数应用于numpy矩阵答案

【问题标题】：Apply function to numpy matrix dependent on position根据位置将函数应用于numpy矩阵
【发布时间】：2021-11-23 12:55:10
【问题描述】：

给定一个形状为 [m,m] 的二维 numpy 数组 X，我希望应用一个函数并获得一个新的二维 numpy 矩阵 P，它的形状也是 [m,m]，其 [i ,j]th 元素得到如下：

P[i][j] = exp (-|| X[i] - x[j] ||**2)

其中||.|| 表示向量的标准 L-2 范数。有没有比简单的嵌套 for 循环更快的方法？

例如，

X = [[1,1,1],[2,3,4],[5,6,7]]

然后，在对角线条目处，访问的行将是相同的，它们的差异的范数/大小将为 0。因此，

P[0][0] = P[1][1] = P[2][2] = exp (0) = 1.0

还有，

P[0][1] = exp (- || X[0] - X[1] ||**2) = exp (- || [-1,-2,-3] || ** 2) = exp (-14)

等等

使用嵌套for循环的最简单的解决方案如下：

import numpy as np
X = np.array([[1,2,3],[4,5,6],[7,8,9]])
P = np.zeros (shape=[len(X),len(X)])
for i in range (len(X)):
    for j in range (len(X)):
        P[i][j] = np.exp (- np.linalg.norm (X[i]-X[j])**2)
        
print (P)

打印出来：

P = [[1.00000000e+00 1.87952882e-12 1.24794646e-47]
    [1.87952882e-12 1.00000000e+00 1.87952882e-12]
    [1.24794646e-47 1.87952882e-12 1.00000000e+00]]

这里，m 是 5e4 的数量级。

【问题讨论】：

@j1-lee 其实不是，X[i] 代表一维向量
为了清楚起见，给我们一个工作示例，必要时使用嵌套循环。然后花点时间阅读numpybroadcasting。
@hpaulj 我刚刚给出了一个简单实现的工作代码

标签： python numpy matrix-multiplication

【解决方案1】：

如果您提供一个示例数组，这会更容易。您可以创建一个大小为 [m, m, m] 的数组 Q，其中 Q[i, j, k] = X[i, k] - X[j, k] 使用

X[None,:,:] - X[:,None,:]

此时，您正在对第三个轴执行简单的 numpy 操作。

【讨论】：

但是 m 是 5e4 的数量级。我现在提供了一个例子。

【解决方案2】：

In [143]: X = np.array([[1,2,3],[4,5,6],[7,8,9]])
     ...: P = np.zeros (shape=[len(X),len(X)])
     ...: for i in range (len(X)):
     ...:     for j in range (len(X)):
     ...:         P[i][j] = np.exp (- np.linalg.norm (X[i]-X[j]))
     ...: 
In [144]: P
Out[144]: 
array([[1.00000000e+00, 5.53783071e-03, 3.06675690e-05],
       [5.53783071e-03, 1.00000000e+00, 5.53783071e-03],
       [3.06675690e-05, 5.53783071e-03, 1.00000000e+00]])

无循环版本：

In [145]: np.exp(-np.sqrt(((X[:,None,:]-X[None,:,:])**2).sum(axis=2)))
Out[145]: 
array([[1.00000000e+00, 5.53783071e-03, 3.06675690e-05],
       [5.53783071e-03, 1.00000000e+00, 5.53783071e-03],
       [3.06675690e-05, 5.53783071e-03, 1.00000000e+00]])

我不得不删除您的 **2 以匹配值。

将norm 应用于 3d 差分数组：

In [148]: np.exp(-np.linalg.norm(X[:,None,:]-X[None,:,:], axis=2))
Out[148]: 
array([[1.00000000e+00, 5.53783071e-03, 3.06675690e-05],
       [5.53783071e-03, 1.00000000e+00, 5.53783071e-03],
       [3.06675690e-05, 5.53783071e-03, 1.00000000e+00]])

在其中一个scikit 包中（学习？）有一个cdist 可以更快地处理这类事情。

【讨论】：

上述方法的空间复杂度是O(m^3)吗？
即使只有 m = 4e3，我也会收到以下内存错误：numpy.core._exceptions.MemoryError: Unable to allocate 93.5 GiB for an array with shape (4000, 4000, 784) and data type浮动64。还有其他解决方案吗？ :(

【解决方案3】：

正如 hpaulj 提到的，cdist 做得更好。请尝试以下操作。

from scipy.spatial.distance import cdist
import numpy as np

np.exp(-cdist(X,X,'sqeuclidean'))

注意sqeuclidean。这意味着 scipy 不采用平方根，因此您不必像上面那样对范数进行平方。

【讨论】：