使用numpy矢量化方程式

Question

我试图将上述公式实现为矢量化形式. K=3在这里,X是150x4numpy数组.mu是3x4numpy数组.Gamma是一个150x3numpy数组.Sigma是一个kx4x4numpy数组.因此Sigma[k]是一个4x4numpy数组.N=150

N_k = np.sum(Gamma, axis=0)
for k in range(K): # Correct
         x_new = X - mu[k] #Correct
         a = np.dot(x_new.T, x_new) #Incorrect from here I feel
         for i in range(len(data)):
             sigma[k] = Gamma[i][k] * a
         sigma[k]=sigma[k]/N_k #totally incorrect

如何解决这个问题？

Answer 1

一个产品？听起来像是np.einsum的工作:

import numpy as np
N = 150
K = 3
M = 4
x = np.random.random((N,M))
mu = np.random.random((K,M))
gamma = np.random.random((N,K))

xbar = x-mu[:,None,:] # shape (3, 150, 4)
sigma = np.einsum('nk,knm,kno->kmo', gamma, xbar, xbar)
sigma /= gamma.sum(axis=0)[:,None,None]

解码'nk,knm,kno->kmo':

此下标规范在array(->)左侧有三个组件,后面是右侧的一个组件.

左边的三个组成部分与标进行对应gamma,xbar并且xbar,被传递到操作数np.einsum.

gamma有下标nk,就像你发布的公式一样. xbar有形状(3,150,4).您可以将其视为具有下标knm,其中k和n具有与您发布的公式中相同的含义,并且m是表示长度为4的轴的下标,该公式未在您的公式中明确提及,但显然已经给出了对阵列的形状.

现在第三个下标组件是kno.使用o下标是因为它o扮演与m下标相同的角色,但我们不希望求和m.事实上,我们希望m和o下标单独迭代,而不是锁定步骤.因此我们给第三个下标不同的字母.

请注意,它n出现在左侧(nk, knm, kno)的下标中,但不出现在右侧(中kmo).这说明np.einsum总结一下n.

该k在左边的标出现,并在右边.这告诉np.einsum我们希望以k锁定步骤推进下标,但是(因为它出现在右边)我们不想总结k.

由于kmo出现在右侧,这些下标仍保留在结果中.这导致sigma具有形状(K,M,M)(即(3,4,4)).