numpy：广播到多个内积和逆

Question

我有数组 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]])

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 个迭代变量，整个迭代 space 可能会非常大。但首先看看这些表达式是否适用于中等大小的数组。

至于

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.

numpy：广播到多个内积和逆

numpy: broadcasting into multiple inner products and inverses

python

numpy

broadcasting

matrix-inverse

inner-product