向量化這個循環在numpy

我想在python與numpy計算矩陣z（定義如下）。向量化這個循環在numpy

這是我目前的解決方案（用於循環使用1）

z = np.zeros((n, k)) 
for i in range(n): 
    v = pi * (1/math.factorial(x[i])) * np.exp(-1 * lamb) * (lamb ** x[i]) 
    numerator = np.sum(v) 
    c = v/numerator 
    z[i, :] = c 
return z

是否有可能完全矢量化這個計算？我需要做數千次迭代的計算，而numpy中的矩陣運算比巨大的循環要快得多。

來源

2017-10-06 Alex Ling

添加最少的樣品情況？ – Divakar

@Divakar感謝評論。可以使用[This gist]（https://gist.github.com/hkalexling/8b97806017cb7cd4ad4937ec1deb157b）作爲示例（python3）。一個EM算法被實現，並且我想向量化E（x）函數。 –

是否「羊肉」不變？ – percusse

這裏是E的矢量化版本。它取代了for-loop和標量運算與NumPy broadcasting和基於陣列的算術：

def alt_E(x): 
    x = x[:, None] 
    z = pi * (np.exp(-lamb) * (lamb**x))/special.factorial(x) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z

我跑em.py獲得意義的x，lamb，pi，n和k的典型大小。在這種規模的數據， alt_E約120倍比E快：

In [32]: %timeit E(x) 
100 loops, best of 3: 11.5 ms per loop 

In [33]: %timeit alt_E(x) 
10000 loops, best of 3: 94.7 µs per loop 

In [34]: 11500/94.7 
Out[34]: 121.43611404435057

這是我使用的基準設置：

import math 
import numpy as np 
import scipy.special as special 

def alt_E(x): 
    x = x[:, None] 
    z = pi * (np.exp(-lamb) * (lamb**x))/special.factorial(x) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z 


def E(x): 
    z = np.zeros((n, k)) 
    for i in range(n): 
     v = pi * (1/math.factorial(x[i])) * \ 
      np.exp(-1 * lamb) * (lamb ** x[i]) 
     numerator = np.sum(v) 
     c = v/numerator 
     z[i, :] = c 
    return z 

n = 576 
k = 2 

x = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5]) 

lamb = np.array([ 0.84835141, 1.04025989]) 
pi = np.array([ 0.5806958, 0.4193042]) 

assert np.allclose(alt_E(x), E(x))

順便說一句，E也可以使用scipy.stats.poisson來計算：

import scipy.stats as stats 
pois = stats.poisson(mu=lamb) 

def alt_E2(x): 
    z = pi * pois.pmf(x[:,None]) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z

但這並不變成是更快，至少對於這個長度的數組：

In [33]: %timeit alt_E(x) 
10000 loops, best of 3: 94.7 µs per loop 

In [102]: %timeit alt_E2(x) 
1000 loops, best of 3: 278 µs per loop

對於較大x，alt_E2更快：

In [104]: x = np.random.random(10000) 

In [106]: %timeit alt_E(x) 
100 loops, best of 3: 2.18 ms per loop 

In [105]: %timeit alt_E2(x) 
1000 loops, best of 3: 643 µs per loop

來源

2017-10-06 15:55:17 unutbu

你是真正的MVP;） –

向量化這個循環在numpy

回答

相關問題