我已经开发了算法,但使用了 for 循环。我正在尝试确定如何利用 Numpy 矢量化。
cov_c = []
for cluster, u, w in zip(r.T, mu_c, total_weight):
s = 0
for n in range(len(d)):
s += cluster[n]*np.outer(d[n] - u, d[n] - u)
cov_c.append(s / w)
cov_c 是一个两元素列表,每个元素都有一个协方差矩阵 (2x2)
[array([[0.19, 0.23],[0.23, 0.39]]),
array([[4.05, -5.01,[-5.018, 6.22]])]
d和r都是二维数组(加权样本)d是特征向量(100个样本的2个特征),其中r是2个高斯每个特征的权重
d.shape
(100, 2)
r.shape
(100, 2)
mu_c 是均值向量的二元素列表
mu_c
[array([ 0.24387682, -0.27793324]), array([ 2.37853451, -1.86454301])]
总重量是一个归一化因子(简单的 2 元素一维数组):
total_weight
array([53.51779102, 46.48220898])
关于如何向量化此计算有什么建议吗?谢谢!
最佳答案
我们可以利用 NumPy 数组来利用矢量化 ufunc 运算。此外,由于 d
中的列数仅为 2
,因此我们只需沿该轴使用循环(因此循环只需两次迭代)。因此,我们将使用切片数据,而不是在各个方向上扩展数组,这会导致更严重的内存拥塞。我们仍然会对切片数据利用广播
。最后,我们将使用np.einsum
来取代外部总和减少,这可能是我们获益最多的地方。
我们最终会得到这样的结果 -
mu_c = np.asarray(mu_c)
total_weight = np.asarray(total_weight)
n = d.shape[1]
out = np.empty((n,2,2))
for i in range(n):
du = d-mu_c[i]
out[i] = np.einsum('i,ij,ik->jk',r[:,i],du,du)
cov_c_out = out/total_weight[:,None,None]
或者,einsum
部分可以替换为矩阵乘法步骤 -
out[i] = (r[:,i,None]*du).T.dot(du)
为了完整性或只是为了好玩,这是一个完全矢量化的解决方案,它是内存密集型的,因此很可能会更慢 -
dmuc = d[:,None,:]-mu_c
out = np.einsum('ij,ijk,ijl->jkl',r,dmuc,dmuc)
此外,还可以通过将 np.einsum
中的 optimize
标志设置为 True
来使用 BLAS。
关于python - 在 Numpy 中应用不使用 for 循环的非平凡矩阵计算,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55460073/