Ускорение алгоритма ЭВМ Гуассиана

Question

Мой код python выглядит следующим образом: он берет навсегда. Могу ли я использовать несколько трюков? Картина я анализирую крошечная и в оттенках серого ...

def gaussian_probability(x,mean,standard_dev): 
     termA = 1.0/(standard_dev*np.sqrt(2.0*np.pi)) 
     termB = np.exp(-((x - mean)**2.0)/(2.0*(standard_dev**2.0))) 
     g = (termA*termB) 
     return g 
def sum_of_gaussians(x): 
     return sum([self.mixing_coefficients[i] * 
        gaussian_probability(x, self.means[i], self.variances[i]**0.5) 
        for i in range(self.num_components)]) 
def expectation(): 
     dim = self.image_matrix.shape 
     rows, cols = dim[0], dim[1] 
     responsibilities = [] 
     for i in range(self.num_components): 
      gamma_k = np.zeros([rows, cols]) 
      for j in range(rows): 
       for k in range(cols): 
        p = (self.mixing_coefficients[i] *  
         gaussian_probability(self.image_matrix[j,k], 
               self.means[i], 
               self.variances[i]**0.5)) 
        gamma_k[j,k] = p/sum_of_gaussians(self.image_matrix[j,k]) 
      responsibilities.append(gamma_k) 
     return responsibilities 

я включил только шаг ожидания, потому что, в то время как шаг максимизации перебирает каждый элемент массива ответственности матриц, кажется, идут относительно быстро (так может быть узким местом является все расчеты gaussian_probability?)

Answer 1

Вы можете значительно ускорить свой расчет, выполнив две вещи:

• не вычислить нормализации в каждом цикле! Как и в настоящее время, для изображения NxN с компонентами M вы вычисляете каждый соответствующий расчет N * N * M раз, что приводит к алгоритму O[N^4 M^2]! Вместо этого вы должны вычислить все элементы один раз, а затем разделить на сумму, которая будет O[N^2 M].

• использование многократной векторизации, а не явных циклов. Это можно сделать очень просто, как вы настроили код.

По сути, ваша expectation функция должна выглядеть следующим образом:

def expectation(self): 
    responsibilities = (self.mixing_coefficients[:, None, None] * 
         gaussian_probability(self.image_matrix, 
              self.means[:, None, None], 
              self.variances[:, None, None] ** 0.5)) 
    return responsibilities/responsibilities.sum(0) 

Вы не обеспечили полный пример, так что мне пришлось импровизировать немного, чтобы проверить и тест, но вот быстрое взятие:

import numpy as np 

def gaussian_probability(x,mean,standard_dev): 
    termA = 1.0/(standard_dev*np.sqrt(2.0*np.pi)) 
    termB = np.exp(-((x - mean)**2.0)/(2.0*(standard_dev**2.0))) 
    return termA * termB 

class EM(object): 
    def __init__(self, N=5): 
     self.image_matrix = np.random.rand(20, 20) 
     self.num_components = N 
     self.mixing_coefficients = 1 + np.random.rand(N) 
     self.means = 10 * np.random.rand(N) 
     self.variances = np.ones(N) 

    def sum_of_gaussians(self, x): 
     return sum([self.mixing_coefficients[i] * 
        gaussian_probability(x, self.means[i], self.variances[i]**0.5) 
        for i in range(self.num_components)]) 

    def expectation(self): 
     dim = self.image_matrix.shape 
     rows, cols = dim[0], dim[1] 
     responsibilities = [] 
     for i in range(self.num_components): 
      gamma_k = np.zeros([rows, cols]) 
      for j in range(rows): 
       for k in range(cols): 
        p = (self.mixing_coefficients[i] *  
         gaussian_probability(self.image_matrix[j,k], 
               self.means[i], 
               self.variances[i]**0.5)) 
        gamma_k[j,k] = p/self.sum_of_gaussians(self.image_matrix[j,k]) 
      responsibilities.append(gamma_k) 
     return responsibilities 

    def expectation_fast(self): 
     responsibilities = (self.mixing_coefficients[:, None, None] * 
          gaussian_probability(self.image_matrix, 
               self.means[:, None, None], 
               self.variances[:, None, None] ** 0.5)) 
     return responsibilities/responsibilities.sum(0) 

Теперь мы можем создать экземпляр объекта и сравнить две реализации на этапе ожидания:

em = EM(5) 
np.allclose(em.expectation(), 
      em.expectation_fast()) 
# True 

Глядя на таймингах, мы про фактор 1000 быстрее для 20x20 изображений с 5 компонентами:

%timeit em.expectation() 
10 loops, best of 3: 65.9 ms per loop 

%timeit em.expectation_fast() 
10000 loops, best of 3: 74.5 µs per loop 

Это улучшение будет расти как размер изображения и количество увеличения компонентов. Удачи!