2

我有数组e，（shape qby l）f（shape nby l）和w（shape nby l），我想创建一个数组MwhereM[s,i,j] = np.sum(w[s, :] * e[i, :] * e[j, :])和一个数组F， where 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]])
4

1 回答 1

3 
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)作品，表演sixi 逆

如果是的话

 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。

于 2016-07-01T05:00:56.320 回答