在 Numpy 中应用不带 for 循环的非平凡矩阵计算

Question

只是想找到为 EM 算法计算这个更新后的协方差矩阵的最佳方法 *

我开发了算法，但使用的是 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 元素 1d 数组）：

total_weight
array([53.51779102, 46.48220898])

关于如何矢量化此计算有什么建议吗？谢谢！

http://www.cse.iitm.ac.in/~vplab/courses/DVP/PDF/gmm.pdf

Answer 1

我们可以利用 NumPy 数组来利用矢量化 ufunc 操作。此外，由于 d 中的列数仅为 2，我们将简单地使用沿该轴的循环（因此只有两次迭代的循环）。因此，我们将使用切片数据而不是在所有方向上扩展数组，这会导致更严重的内存拥塞。我们仍然会在切片数据上利用 broadcasting。最后，我们将使用 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。

在 Numpy 中应用不带 for 循环的非平凡矩阵计算

Applying non-trivial matrix calculating without for-loop in Numpy

python

arrays

numpy

matrix

vectorization