使用numpy的

Question

我試圖實現上述公式作爲向量化形式向量化的方程。 K=3這裏，X是150x4 numpy數組。 mu是3x4 numpy數組。 Gamma是一個150x3 numpy數組。 Sigma是一個kx4x4 numpy數組。因此Sigma[k]是一個4x4 numpy數組。 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'：

此下標規範具有三個分量到後跟一個分量，其右側的陣列（->）的左側。

左邊的三個分量對應與標爲gamma，xbar和xbar，被傳遞到np.einsum操作數。

gamma具有下標nk，就像在您發佈的公式中一樣。 xbar有形狀（3,150,4）。你可以認爲它有下標knm，其中k和n與你發佈的公式中的含義相同，而m是一個代表長度軸的下標，這在你的公式中沒有明確提到，但顯然是給定的你對數組形狀的描述。

現在第三個下標成分是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））。

Answer 2

只是在unubtu偉大的答案頂部的幾個性能指針。

np.einsum無法優化具有兩個以上參數的調用。只要有可能，它通常更快的計算手動分成的兩個參數組是，例如：

def unubtu(): 
    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] 
    return sigma 

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

In [50]: %timeit unubtu() 
10000 loops, best of 3: 147 µs per loop 

In [51]: %timeit faster() 
10000 loops, best of 3: 129 µs per loop 

提高了12％並不多，但差異可以得到（多少）有更大的陣列大。

此外，即使np.einsum是一個很好的工具，它使得很容易的事情變得非常困難，但它並不像np.dot那樣優化，如果你的numpy是用一個好的線性代數庫構建的話。在你的情況，考慮到K小，使用np.dot和Python的循環是更快：

def even_faster(): 
    sigma = np.empty((K, M, M)) 
    for k in xrange(K): 
     x_ = x - mu[k] 
     sigma[k] = np.dot((x_ * gamma[:, k, None]).T, x_) 
    sigma /= gamma.sum(axis=0)[:,None,None] 
    return sigma 

In [52]: %timeit even_faster() 
10000 loops, best of 3: 101 µs per loop