Применение нетривиального матричного вычисления без цикла for в Numpy - PullRequest
Новая
       15

Применение нетривиального матричного вычисления без цикла for в Numpy

1 голос
/ 01 апреля 2019

Просто пытаюсь найти лучший способ вычислить эту обновленную ковариационную матрицу для EM-алгоритма enter image description here*

Я разработал алгоритм, но использую цикл 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 - вектор признаков (2 элемента для 100 выборок), где 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])

Любые предложения о том, как векторизовать это вычисление?Спасибо!

1 Ответ

2 голосов
/ 01 апреля 2019

Мы могли бы использовать массивы 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)

Кроме того, поэкспериментируйте с флагом optimize в np.einsum, установив его как True для использования BLAS.

...