计算矩阵中的每个元素

Question

这里是我需要为Numpy矩阵的每个元素计算的公式：

Mi_j = Sum_v(Av * Xi_v) + Sum_v(Bv * Wj_v) + Sum_v(Gv * Zij_v)

我不太清楚如何用numpy的方式编写它（在python中太长了）：矢量化/切片/C API。

你有什么建议，能给我一个简单的例子吗？我对numpy很新。

@编辑索引

- A、B、G是一维数组[x,x,x] - Xi和Wj也是如此（X是一个矩阵，W是一个矩阵） - Zij是一个一维数组

- Touki

---

“Sum(Av * Xi_v)”的意思是你正在对“v”索引求和吗？如果是这样，那么“Sum(Gu * Zij_v)”又代表什么意思呢？ - unutbu

...并且u也是一个索引吗？ - askewchan

我已经编辑了这条消息。 - Touki

输入的 X 和 W 是否作为完整矩阵，Z 是否作为完整的三维数组，如公式所假设的那样，还是像最后列出的 Xi、Wj 和 Zij 一样作为1d数组？我的答案允许完整的数组。 - askewchan

Answer 1

我个人认为，先确定代数过程，然后使用numpy矩阵作为标准来执行它们，会更易读。如果你的工作涉及到数学，使用numpy的matrix类将更容易将数学转换为代码，反之亦然。

这也将帮助你避免仔细广播。

从以下内容开始：

Mi_j = Sum_v(Av * Xi_v) + Sum_v(Bv * Wj_v) + Sum_v(Gv * Zij_v)

在numpy中，它变成了：

M = X*A + (W*B).T + Z*G

如果您将每个矩阵初始化为np.matrix，适当的代数运算会自动完成。

import numpy as np
N = 5

A = np.asmatrix(np.arange(N)).T
B = np.asmatrix(np.arange(N)).T
G = np.asmatrix(np.arange(N)).T

X = np.asmatrix(np.arange(N*N).reshape(N,N))
W = np.asmatrix(np.arange(N*N).reshape(N,N))

Z = np.asmatrix(np.arange(N**3).reshape(N,N,N))

请注意，我已经对一维矩阵进行了转置，因为一维矩阵默认是行向量。真正的向量应该是列向量。之后你就不需要再担心广播问题了。

M = X*A + (W*B).T + Z*G
print M
[[  90  190  290  390  490]
 [ 390  490  590  690  790]
 [ 690  790  890  990 1090]
 [ 990 1090 1190 1290 1390]
 [1290 1390 1490 1590 1690]]

- askewchan

Answer 2

让我们通过一个简单的例子来进行说明：

如果我们定义：

import numpy as np
N = 5
A = np.arange(N)
X = np.arange(N*N).reshape(N,N)

B = np.arange(N)
W = np.arange(N*N).reshape(N,N)

G = np.arange(N)
Zij = np.arange(N)

那么第一个和，Sum_v(Av * Xi_v) 可以使用 np.dot 计算：

In [54]: X
Out[54]: 
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

In [55]: A
Out[55]: array([0, 1, 2, 3, 4])

In [56]: np.dot(X, A)
Out[56]: array([ 30,  80, 130, 180, 230])

同样地，第二个总和 Sum_v(Bv * Wj_v) 可以计算为：

In [58]: np.dot(W,B)
Out[58]: array([ 30,  80, 130, 180, 230])

然而，我们希望第一项求和结果是一个沿着i索引变化的向量，而我们希望第二项求和结果是一个沿着j索引变化的向量。为了在numpy中实现这一点，使用广播：

In [59]: np.dot(X,A) + np.dot(W,B)[:,None]
Out[59]: 
array([[ 60, 110, 160, 210, 260],
       [110, 160, 210, 260, 310],
       [160, 210, 260, 310, 360],
       [210, 260, 310, 360, 410],
       [260, 310, 360, 410, 460]])

第三个总和是两个一维数组之间的简单点积：

In [60]: np.dot(Zij, G)
Out[60]: 30

把它们全部放在一起，就是这样：

In [61]: M = np.dot(X,A) + np.dot(W,B)[:,None] + np.dot(Zij, G)

In [62]: M
Out[62]: 
array([[ 90, 140, 190, 240, 290],
       [140, 190, 240, 290, 340],
       [190, 240, 290, 340, 390],
       [240, 290, 340, 390, 440],
       [290, 340, 390, 440, 490]])

---

请注意，我可能误解了Zij的含义。虽然您说它是一个一维数组，但也许您的意思是对于每个i,j，它都是一个一维数组。那么Z将是三维的。

为了具体起见，假设Z的前两个轴表示i和j索引，而Z的最后一个轴是您想要求和的轴。

在这种情况下，您需要的最后一个术语是np.dot(Z, G):

In [13]: Z = np.arange(N**3).reshape(N,N,-1)

In [14]: np.dot(X,A) + np.dot(W,B)[:,None] + np.dot(Z, G)
Out[14]: 
array([[  90,  190,  290,  390,  490],
       [ 390,  490,  590,  690,  790],
       [ 690,  790,  890,  990, 1090],
       [ 990, 1090, 1190, 1290, 1390],
       [1290, 1390, 1490, 1590, 1690]])

- unutbu

---

作为一点小小的挑剔，你使用 None 作为 newaxis 的原因是什么呢？（非常好的答案，无论如何加上一票） - ev-br

这假设Z_ij对于所有的i和j都是常数，但一般情况下并非如此。 - askewchan

@askewchan：谢谢，你可能是对的。我已经相应地编辑了我的帖子。顺便说一下，虽然矩阵类使前两个术语看起来更漂亮，但如果Z是三维的，您将无法使用矩阵类进行“Z”。Z = np.asmatrix(...reshape(N,N,N))会抛出一个ValueError。 - unutbu

个人而言，我更喜欢将矩阵保留为“matrix”，因为这样我就不会意外地执行“A*X”操作，因为它会发出警告，因为内部索引的维度不相同。严格的类型有助于防止错误，并允许进行外积等操作，而不需要太多思考（对我来说）。 - askewchan

我更喜欢使用普通的ndarrays而不是矩阵，因为矩阵只适用于二维。虽然这使得一些表达式看起来更像数学，但符号并不具有普适性。我宁愿不编写假定我的输入是二维的函数。 - unutbu