numpy：广播成多个内积和逆答案

【问题标题】：numpy: broadcasting into multiple inner products and inversesnumpy：广播成多个内积和逆
【发布时间】：2016-07-01 03:14:42
【问题描述】：

我有数组e、（形状q by l）f（形状n by l）和w（形状n by l）和我想创建一个数组M，其中M[s,i,j] = np.sum(w[s, :] * e[i, :] * e[j, :])，以及一个数组F，其中F[s,j] = np.sum(w[s, :] * f[s, :] * e[j, :])。

这两者都很容易做到，例如，循环遍历M 的元素，但我想提高效率（我的真实数据有类似 1M 长度为 5k 的条目）。对于F，我可以使用F = np.inner(w * f, e)（我验证它产生与循环相同的答案）。 M 更难，所以第一步是循环通过列表理解的维度零，说M = np.stack([np.inner(r[:] * e, e) for r in w])（我已经验证这也与循环相同）。 np.inner() 不接受任何轴参数，所以我不清楚如何告诉数组只广播 w 的所有行。

最后，我需要使用M和F的元素来创建一个矩阵A，其中A[s,i] = np.sum(np.linalg.inv(M[s, :, :])[i, :] * F[i, :])。这看起来也是内积式的，但是进行大量单个逆运算很耗时，那么有没有一种方法可以计算切片的逆运算，而不是循环？

我的数组中的一些测试值如下：

e = array([[-0.9840087 , -0.17812043],
           [ 0.17812043, -0.9840087 ]])

w = array([[  1.12545297e+01,   1.48690140e+02],
           [  7.30718244e+00,   4.07840612e+02],
           [  2.62753065e+02,   2.27085711e+02],
           [  1.53045364e+01,   5.63025281e+02],
           [  8.00555079e+00,   2.16207407e+02],
           [  1.54070190e+01,   1.87213209e+06],
           [  2.71802081e+01,   1.06392902e+02],
           [  3.46300255e+01,   1.29404438e+03],
           [  7.77638140e+00,   4.18759293e+03],
           [  1.12874849e+01,   5.75023379e+02]])

f = array([[ 0.48907404,  0.06111084],
           [-0.21899297, -0.02207311],
           [ 0.58688524,  0.05156326],
           [ 0.57407751,  0.10004592],
           [ 0.94172351,  0.03895357],
           [-0.7489003 , -0.08911183],
           [-0.7043736 , -0.19014227],
           [ 0.58950925,  0.16587887],
           [-0.35557142, -0.14530267],
           [ 0.24548714,  0.03221844]])

【问题讨论】：

在M[s,i,j] = np.sum(w[s] * e[i] * e[j]) 中的总和是多少？我在两边都看到了i,j,s。是(w[s,k]*e[i,k]*e[j,k]).sum(axis=k)吗？
np.einsum 是一种快速简便的工具，可用于计算此类产品的总和。 s,i,j 方程几乎是自己写的。
w、f 和 e 的行是长度-l 向量，因此三者的元素乘积在其唯一的轴上求和
有时写e[j,:] 比写e[j] 更清晰，即使解释器对两者的处理方式相同。

标签： python numpy matrix-inverse broadcasting inner-product

【解决方案1】：

M[s,i,j] = np.sum(w[s, :] * e[i, :] * e[j, :])

翻译成

M = np.einsum('sk,ik,jk->sij',w,e,e)

和

F[s,j] = np.sum(w[s, :] * f[s, :] * e[j, :])
F = np.einsum('sk,sk,jk->sj', w, f, e)

我没有用你的样本测试过这些，但是翻译很简单。

对于真正的大型数组，您可能必须将表达式分解为多个部分。使用 4 个迭代变量，整个迭代空间可能非常大。但首先看看这些表达式是否适用于中等大小的数组。

至于

A[s,i] = np.sum(np.linalg.inv(M[s, :, :])[i, :] * F[i, :])

我看起来像np.linalg.inv(M) 工作，正在执行s i x i inverses

如果是的话

 IM = np.linalg.inv(M)
 A = np.einsum('skm,ik,im->si', IM, F)

我猜的更多。

同样，维度可能会变得太大，但请先尝试变小。

建议使用典型的线性方程解而不是直接逆解，例如

  A = F/M
  A = np.linalg.solve(M, F)

因为你可能想要A 这样M@A=F（@矩阵乘积）。但我对这些事情有点生疏了。还要检查tensorsolve 和tensorinv。

【讨论】：

A = np.einsum('skm,ik,im->si', IM, F) 的索引太多，但A = np.linalg.solve(M, F) 生成的数组具有正确的维度，成功了一半。
'sij,ij->si` 更好吗？
不，索引调整会导致操作数再次广播错误[original->remapped]: (50,2,2)->(50,2,2) (50,2)->(50,2)。但是，使用np.linalg.solve() 实现，当w 无处不在时，此方法产生与f.dot(e.T) 相同的答案。这是预期的行为，所以我很确定它有效。