python - numpy:广播到多个内积和逆

Question

我有数组e，(形状q by l)f(形状n by l) 和 w (shape 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 更困难，因此第一步是使用列表理解循环遍历 的 0 维，表示 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]])

Answer 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逆

如果是这样

 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。