Сгруппируйте похожие изображения с помощью модели смеси Гаусса (алгоритм EM)

Внедрить GMM с нуля для решения проблемы кластеризации изображений

Кластеризация - одна из самых популярных задач машинного обучения без учителя. Мы уже знакомы с алгоритмом кластеризации k-средних, но подождите, в этом есть проблема:

1. Сильно зависит от начального значения центроидов: если мы изменим его инициализацию, а затем положение кластера, то весьма вероятно, что последний кластер изменит свое положение.
2. Кластеры разного размера и плотности: кластеры, имеющие различную форму (разброс f точек данных) и размер (количество точек данных), не могут быть обработаны с помощью этого простого метода.
3. Точки выбросов. Может случиться, что точка данных находится очень далеко от любого из кластеров, и она может выступать в качестве отдельного кластера. Чтобы избежать таких проблем, мы обычно удаляли эти точки.

Если вы не хотите разбираться в этой мерзкой математике, прокрутите вниз и перейдите к разделу Внедрение GMM, где я расскажу о реализации на Python.

Чтобы понять модель гауссовой смеси, вам необходимо понять несколько терминов.

1. Оценка плотности: оценка плотности вероятности - это построение приближенной на основе наблюдаемых данных ненаблюдаемой основной функции плотности вероятности. Этот шаг включает в себя выбор функции распределения вероятностей и параметра этой функции, который описывает наилучшее совместное распределение вероятностей наблюдаемых данных.
2. Оценка максимального правдоподобия (MLE): это метод оценки параметров распределения вероятностей путем максимизации функции логарифмического правдоподобия.
3. Скрытые переменные: они противоположны наблюдаемым переменным, то есть не наблюдаются непосредственно, а скорее выводятся из математической модели. Об этом поговорим позже.

Теперь оценка максимального правдоподобия не очень хорошо работала при наличии скрытых переменных. Алгоритм максимизации ожидания - это один из методов, с помощью которого можно найти подходящие параметры модели при наличии скрытых переменных.

EM - это итерационный метод для поиска максимальной вероятности или максимальных апостериорных оценок параметров, где модель зависит от ненаблюдаемых скрытых переменных.

Позвольте немного упростить это.

Учитывая N наблюдение

и вероятностная модель p (X | θ)

Наша цель - найти значение θ, которое максимизирует p (X | θ). Это называется оценкой максимального правдоподобия.

ЕСЛИ модель представляет собой простое распределение Гаусса, тогда параметр будет иметь вид

где

Если это так, то наиболее распространенным подходом к решению этой проблемы в ML является градиентный спуск путем минимизации логарифмической функции правдоподобия, которая равна –log (p (X | θ)) как функция потерь (помните, что мы фактически максимизируем функцию вероятности при параметре

это простой прием для обеспечения числовой стабильности, поскольку он превращает произведение в сумму, и минимизация отрицательного значения этой функции - то же самое, что максимизация этой функции).

Но градиентный спуск - это метод решения нелинейных задач, обычно для минимума некоторого многомерного распределения. Вычисление среднего (или других моментов) гауссовой смеси не является нелинейной задачей и поэтому не требует методов, разработанных для нелинейных задач. Так что есть способ лучше решить эту проблему.

А вот и алгоритм EM. Давайте посмотрим, как алгоритм EM используется в модели гауссовой смеси.

Оценка максимального правдоподобия (MLE) может быть упрощена путем введения скрытой переменной. Модель скрытых переменных предполагает, что наблюдение xi вызвано некоторой лежащей в основе скрытой переменной. Хммм… все еще не ясно, хорошо, рассмотрим это изображение

Этот тип вероятностной модели называется моделью гауссовой смеси (GMM), взвешенная сумма C гауссовых компонентов в этом примере C = 3 равна

πc и Σc - весовой вектор, средний вектор и ковариационная матрица компонента смеси Cc, соответственно. Веса неотрицательны и в сумме равны 1, т. Е.

Вектор параметра

обозначает набор всех параметров модели. Если мы введем дискретную скрытую переменную 't', которая определяет присвоение наблюдений компонентам смеси, мы можем определить совместное распределение по наблюдаемым и скрытым переменным p (x, t | θ) в терминах условного распределения p (x | t, θ) и априорное распределение p (t | θ)

Где

Значения t закодированы в горячем режиме. Например, t2 = 1 относится ко второму компоненту смеси, что означает t = (0, 1, 0), если всего имеется C = 3 компонента. Маргинальное распределение p (x | θ) получается суммированием по всем возможным состояниям t.

Для каждого наблюдения xi у нас есть одна скрытая переменная ti, которую также называют ответственностью.

Рассматривая X как все наборы наблюдений, набор всех скрытых переменных T, мы можем легко максимизировать полную логарифмическую правдоподобие данных p (X, T | θ). Поскольку после нахождения логарифма правдоподобия мы затем узнаем кластерные присвоения точек, поэтому мы найдем предельное логарифмическое правдоподобие или неполное логарифмическое правдоподобие данных p (X | θ). Таким образом, математически,

Внедрение GMM

Давайте посмотрим шаг за шагом, как наше изображение кластеризуется с помощью модели гауссовой смеси. Я использую здесь python для реализации модели GMM:

Требуется внешняя библиотека Python:

1. imageio: для получения функций RGB из изображения
2. pandas: для обработки набора данных
3. numpy: для математических операций

Шаг 1:

Давайте начнем с нашего набора данных, мы указали путь к папке, в которой находится все изображение. Мы создадим метод, из которого мы сможем извлечь среднее значение вектора R, G B для каждого изображения. Для одного изображения функция должна быть:

######## loading external package dependency ####################
Import pandas as pd
Import numpy as np
from scipy.stats import multivariate_normal
Import imageio
From functools import reduce
def get_image_feature(path):
    Im = imageio.imread(os.path.join(path), pilmode=’RGB’)
    temp = Im/255. # divide by 255 to get in fraction
    mn = temp.sum(axis=0).sum(axis=0)/(temp.shape[0]*temp.shape[1])
    return mn/np.linalg.norm(mn, ord=None) # taking 2nd norm to scale vector

Пусть предоставленный путь изображения - «(some_path) /Image_SunSet.jpg», эта функция возвращает [Red, Green, Blue] интенсивность пикселей изображения как [0,9144867979, 0,3184891297, 0,2495567485], мы можем ясно видеть, что красный цвет имеет большую интенсивность.

Точно так же мы можем сделать это для всех наших изображений и создать фрейм данных pandas, содержащий детали изображения:

• file_name: имя файла
• path: путь к файлу
• extension: расширение изображения, например .jpg, .png и т. д.
• ImageBase64: кодирование изображения Base64
• красный: Красный цвет Средняя интенсивность
• зеленый: зеленый цвет Средняя интенсивность
• синий: Синий цвет Средняя интенсивность

Теперь у нас есть набор данных, содержащий функции (т.е. [красный, зеленый, синий] столбец набора данных, не все остальные данные используются только для просмотра). Если вы хотите использовать большее количество k (количество кластеров), вы можете создать свою собственную функцию, извлеченную из изображения, но представить более широкий диапазон цветового кода, а не только r, g и b. но для простоты я взял эти 3 основных цветовых вектора как функции.

Шаг 2:

Генерация модели (реализация алгоритма EM для моделей гауссовой смеси):

Поскольку в GMM мы предполагаем, что все точки данных вместе взяты из k гауссовского распределения, где k - не что иное, как количество кластеров. Хммм… немного сложновато, позвольте немного упростить.

Рассмотрим это одномерное распределение вероятностей. Если мы предположим, что все точки данных поступают из двух распределений Гаусса, это означает, что мы предполагаем, что k = 2 (количество кластеров), тогда, если мы хотим знать, какие данные были получены из какого распределения Гаусса, все, что нам нужно сделать, это вычислить среднее значение и дисперсия как для гауссовского распределения. Но подождите, мы еще не знаем. Но на секунду, если мы предположим, что знаем эти параметры (то есть среднее значение и дисперсию k-го гауссовского распределения), тогда мы можем узнать вероятность каждой точки данных, из которой оно происходит по гауссовскому закону (эта вероятность относится к ответственности). Ого, но мы этого тоже не знаем. Итак, наш вывод

Если мы знаем параметр (то есть среднее значение и дисперсию по Гауссу), мы можем рассчитать обязанности, а если мы знаем ответственность, мы можем рассчитать параметры.

Вот тут-то и появляется алгоритм EM, где он начинается с размещения «k» гауссианов случайным образом (генерирует случайное среднее значение и дисперсию для гауссианов) в картинку, где он выполняет итерацию между шагами E и M, пока не сойдется, где:

Шаг E: распределите обязанности кластера с учетом текущих параметров

Шаг M: обновление параметров с учетом текущих обязанностей кластера

Поговорим о его реализации подробнее.

Вспомним, для чего используется GMM. В отличие от k- означает, что он не выполняет жестких назначений кластеру, скажем, например, если изображение

Если мы работаем с k = 4, т.е. 4 кластерной моделью, то есть с закатом, лесом, облаками и небом, то в алгоритме k mean эта картинка принадлежит только облаку, это называется жестким назначением, но в GMM оно принадлежит 97% облаков и 3% лес это называется мягким назначением.

Шаг E:

На этом этапе мы вычисляем обязанности кластера. Пусть rik представляет ответственность кластера k за точку данных i. Вы можете принять на себя ответственность как вероятность того, что точка данных «i» принадлежит кластеру «k». Поскольку это вероятность, это означает

Чтобы узнать, насколько кластер отвечает за данную точку данных, мы вычисляем вероятность точки данных при конкретном назначении кластера, умноженную на вес кластера. Для точки данных i и кластера k математически:

Где N (xi | μk, Σk) - гауссово распределение для кластера k (со средним значением μk и ковариацией Σk).

Мы использовали ∝, потому что величина N (xi | μk, Σk) еще не является той ответственностью, которую мы хотим. Чтобы гарантировать, что все обязанности по каждой точке данных в сумме составляют 1, мы добавляем константу нормализации в знаменатель:

Но наши данные - это трехмерный вектор, то есть [красный, зеленый, синий]. Scipy предоставляет удобную функцию для вычисления многомерного нормального распределения. "Проверь это"

то есть multivariate_normal.pdf ([точка данных], среднее значение = [вектор среднего], cov = [матрица ковариаций])

Позже мы поговорим о том, как вычислить средний вектор и ковариационную матрицу на шаге M.

# data (numpy array) : array of observations
# weights (numpy array) : numpy array of weight of each clusters of size (1, n_clusters)
#means (numpy array) : numpy array of means of each clusters of size (n_cluster, dimension)
#covariances(numpy array) : numpy array of covariance metrix of size (n_clusters, dimension, dimension)
def get_responsibilities( data, weights, means, covariances):
    n_data = len(data)
    n_clusters = len(means)
    resp = np.zeros((n_data, n_clusters))
    for i in range(n_data):
       for k in range(n_clusters):
          resp[i, k] = weights[k]*   multivariate_normal.pdf(data[i],means[k],covariances[k],allow_singular=True)
        # Add up responsibilities over each data point and normalize
    row_sums = resp.sum(axis=1)[:, np.newaxis]
    resp = resp / row_sums
    return resp

Помните, что эта ответственность - наша скрытая переменная.

M-Step:

Теперь рассчитана ответственность кластера, мы должны обновить параметры кластера, то есть (веса (πk), средние (μk) и ковариации (Σk)), связанные с каждым кластером.

Обновить веса:

Вес кластера показывает, насколько каждый кластер представляет все точки данных. Математически это определяется как:

код Python

# resp(numpy array) : responsibility numpy array size (n_sample, n_clusters)
def get_soft_counts(resp):
    return np.sum(resp, axis=0)
# counts (numpy array) : count list of sum of soft counts for all clusters of size (n_cluster)
def get_weights(counts):
    n_clusters = len(counts)
    sum_count = np.sum(counts)
    weights = np.array(list(map(lambda k : counts[k]/sum_count, range(n_clusters))))
    return weights

Обновления означают:

Среднее значение каждого кластера равно средневзвешенному значению всех точек данных, взвешенных по обязанностям кластера. Математически для каждой точки данных xi и ответственности за k-й кластер rik is Среднее значение k-го кластера можно определить как:

код Python

# data (numpy array): array of observation points
# resp (numpy array) : responsibility numpy array size (n_sample, n_clusters)
# counts (numpy array) : count list of sum of soft counts for all clusters of size (n_cluster) 
def _get_means( data, resp, counts):
    n_clusters = len(counts)
    n_data = len(data)
    means = np.zeros((n_clusters, len(data[0])))
    
    for k in range(n_clusters):
        weighted_sum = reduce(lambda x,i : x + resp[i,k]*data[i],  range(n_data), 0.0)
        means[k] = weighted_sum/counts[k]
    return means

обновить ковариацию:

Ковариация каждого кластера равна средневзвешенному значению всех внешних продуктов, взвешенных с учетом обязанностей кластера. Математически это определяется как:

Где

это внешний продукт. Давайте рассмотрим простой пример того, как выглядит внешний продукт.

Позволять

- два вектора, то внешнее произведение a и b определяется как:

Транспонирование предназначено только для того, чтобы преобразовать вектор-столбец в вектор-строку. Поскольку по соглашению вектор представлен как вектор-столбец в машинном обучении.

код Python

# data (numpy array) : array of observation points
# resp(numpy array) : responsibility numpy array size (n_sample, n_clusters)
# counts (numpy array) : count list of sum of soft counts for all clusters of size (n_cluster)
# means (numpy array) : numpy array of means of each clusters of size (n_cluster, dimension)
def _get_covariances( data, resp, counts, means):
    n_clusters = len(counts)
    dimension = len(data[0]) # to get dimention of data
    n_data = len(data)
    covariances = np.zeros((n_clusters, dimension, dimension))
    
    for k in range(n_clusters):
       weighted_sum = reduce (lambda x, i :x + resp[i,k] *  np.outer((data[i]-means[k]), (data[i]- means[k]).T), range(n_data), np.zeros((dimension, dimension)))
       
       covariances[k] = weighted_sum /counts[k] # normalize by total sum of counts
    return covariances

Но как мы будем измерять нашу смесь по гауссовской модели. Для этого мы вычислим логарифмическую вероятность гауссовой смеси. Он количественно определяет вероятность наблюдения заданного набора данных при определенных параметрах, то есть средних значениях, ковариации и весах. Мы продолжим повторять алгоритмы EM с разными наборами параметров и проверять их сходимость, то есть логарифмическое правдоподобие не меняется слишком сильно или до тех пор, пока не будет зафиксировано количество итераций.

Код Python

# Z (numpy array) : numpy array of size (1, n_clusters)
def _ln_sum_exp( Z):
    return np.log(np.sum(np.exp(Z)))
# data (numpy array) : array of observation points
# weights (numpy array) : numpy array of weight of each clusters ofsize (1, n_clusters)
# means (numpy array) : numpy array of means of each clusters of size (n_cluster, dimension)
# covs (numpy array) : numpy array of covariance metrix of size (n_clusters, dimension, dimension)
def get_log_likelihood( data, weights, means, covs):
    n_clusters = len(means)
    dimension = len(data[0])
    sum_ln_exp = 0
    
    for d in data:
      
      Z = np.zeros(n_clusters)
      for k in range(n_clusters):
        # compute exponential term in multivariate_normal
        delta = np.array(d) — means[k]
        
        inv = np.linalg.inv(covs[k])
        exponent_term = np.dot (delta.T, np.dot (inv, delta))
        # Compute loglikelihood contribution for this data point and this cluster
        Z[k] += np.log (weights[k])
        det = np.linalg.det(covs[k])
        
        Z[k] -= 1/2. * (dimension * np.log (2*np.pi) + np.log (det) + exponent_term)
        # Increment loglikelihood contribution of this data point across all clusters
      
      sum_ln_exp += _ln_sum_exp(Z)
    return sum_ln_exp

Пришло время собственно основного метода реализации алгоритмов ЭМ.

Код Python

# data (numpy array) : array of observation points
# init_weights (numpy array) : numpy array of initial weight of each clusters ofsize (1, n_clusters)
# init_means (numpy array) : numpy array of initial means of each clusters of size (n_cluster, dimension)
# init_covariances (numpy array) : numpy array of initial covariance metrix of size (n_clusters, dimension, dimension)
# maxiter (int) : maximum iteration to rum EM (optional)
# threshold (float) : maximum threshold to stop EM (optional)
def em_from_parameter(data, init_means, init_covariances, init_weights, maxiter=1000, thresh=1e-4):
    # Make copies of initial parameters
    means = init_means[:]
    covariances = init_covariances[:]
    weights = init_weights[:]
    
    # Infer length of dataset
    n_data = len(data)
    
    # Infer number of cluster
    n_clusters = len(means)
    
    # Initialize some useful variables
    resp = np.zeros((n_data, n_clusters), dtype=np.float64)
    
    l1 = get_log_likelihood(data, weights, means, covariances)
    l1_list = [l1]
    
    for it in range(maxiter):
      
      # E-step: calculate responsibilities
      
      resp = get_responsibilities(data, weights, means, covariances)
      
      # M-step calculate cluster parameter
      counts = get_soft_counts(resp)
      weights = _get_weights(counts)
      means = _get_means(data, resp, counts)
      covariances = _get_covariances(data, resp, counts, means)
      
      l1_new = get_log_likelihood(data, weights, means, covariances)
      
      l1_list.append(l1_new)
      
      if abs(l1_new – l1) < thresh :
        break
      
      l1 = l1_new
    
    param = {'weights': weights, 'means': means, 'covariances': covariances, 'loglikelihood': l1_list, 'responsibility': resp, 'Iterations': it}
    return param

Итак, все настроено, но какой будет наш начальный параметр. Мы можем использовать множество подходов.

Один из способов - инициализировать веса всех кластеров одинаково, то есть (1 / k) для каждого кластера, где k - общее количество кластеров.

Инициализация ковариационной матрицы для 1 кластера представляет собой квадратную диагональную матрицу размерности, равной размеру данных * размерности данных (в данном случае размерность данных равна 3, т.е. [r, g, b]) математически.

Где σxy = - ковариация между x и y

Для средней инициализации это немного сложно. Помните алгоритм k-среднего. Мы используем выходные данные этих алгоритмов для инициализации нашего среднего параметра EM. Поскольку, используя алгоритм k-mean, мы продвигаем координату среднего, чтобы достичь ближайшей точки сходимости. Так что, если мы начнем с этого момента, алгоритм EM займет меньше времени, чтобы сойтись. Но есть проблема с тем, как мы инициализируем начальные центроиды для k-среднего.

Хм. Что ж, для этой проблемы мы воспользуемся стратегией инициализации k mean ++, в которой она выбирает k центроидов, наиболее удаленных друг от друга. Теперь у нас есть все три параметра для начала.

ПРИМЕЧАНИЕ. Для полного решения посмотрите мою ссылку на git-hub PYTHONAPI, вы найдете полный код Flask для этой проблемы кластеризации от начала до конца и может быть использован любой технологией на стороне клиента. Я выбираю Angular 8 как клиентскую сторону, код которой доступен на ML_algorithmUI.

Для тех, кто хочет проверить только код, связанный с EM, ссылка эта содержит две папки Expectation_Maximization и K_Mean, которые содержат относительный код отдельно. Просмотрите код и не стесняйтесь спрашивать меня, если у вас возникнут проблемы.

Теперь вернемся к коду, так что код драйвера этой программы

Код Python

# k: number of clusters
# data: observations coming from distribution
# return final weights, means, covariances, list of loglikelihood till convergence and final responsibility of each observation and number of Iterations it takes till convergence
def em(k, data) :
    dimension = len(data[0])
    # mean co-ordinate of all clusters not recommended at all since it is a very bad method to initialized the mean points alternatively use the k-mean initializing method as discussed early
    means = np.random.rand(k, dimension)
    #### alternate better method to initialize mean points by kmean ####
    ## from KM import KM ## self created KM class available in my github repos
    ## check K_Mean folder in listed github repos path
    km = KM()
    means, dont_require = km.kmeans(data, k) # return centroids and cluster assignments
    
    cov = np.diag(np.var(data, axis=0)) # get covariance of data initialize by diagonal metrix of covariance of data
    
    covariances = np.array([cov] * k) # initiate covariance metrix with diagonal element is covariance of data
    weights = np.array([1/ k] * k) # initialize equal weight to all clusters
    return em_from_parameter(data, means, covariances, weights)

images_path = ['../image2.jpg', '../image2.jpg',…] # список путей изображения ……… .. …………………………………………… …………………………… .. (1)

########### Data generation ################
data = []
For path in images_path:
data.append(get_image_feature(path))
data = np.array(data)

Теперь наши данные готовы, и значение k, то есть номер кластера, является личным выбором. Для простоты я беру k = 4, что означает, что мы предполагаем, что распределение данных происходит от 4 гауссовских значений. Поскольку наши данные и значение k готовы, мы готовы создавать модель.

####### create model #######
parameters = em(k, data) # store returned parameters of model along with responsibilities

. . . . . . ………………………………………………………………………….(2)

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

Код Python

# data (numpy array) : array of observation points
# weights (numpy array) : numpy array of weight of each clusters of size (1, n_clusters)
# means (numpy array) : numpy array of means of each clusters of size (n_cluster, dimension)
# covs (numpy array) : numpy array of covariance metrix of size (n_clusters, dimension, dimension)
def get_responsibilities(self, data, weights, means, covariances):
    n_data = len(data)
    n_clusters = len(means)
    
    resp = np.zeros((n_data, n_clusters))
for i in range(n_data):
for k in range(n_clusters):
resp[i, k] = weights[k]* multivariate_normal.pdf(data[i],means[k],covariances[k])
# Add up responsibilities over each data point and normalize
row_sums = resp.sum(axis=1)[:, np.newaxis]
resp = resp / row_sums
return resp

создайте свой массив наблюдений и вызовите этот метод с параметром, возвращаемым em, как показано в (2), и передайте его этой функции, вы получите мягкое назначение каждого кластера. Например, пусть у вас есть изображение с путем «../../test_example.jpg»

сначала выберите вектор [r, g, b] изображения, как показано в (1), вы получите что-то вроде [0.9144867979, 0.3184891297, 0.2495567485], передайте эти данные в функцию get_responsabilities вместе с возвращаемыми параметрами кластеров. Вы получите это в переменной параметров, т.е.

примечание: помните, что эта функция принимает массив наблюдений, и в этом случае у нас есть только одно изображение, то есть одно наблюдение, поэтому сначала преобразуйте это одиночное наблюдение в массив, а затем вызовите функцию.

data = np.array([ np.array([0.9144867979, 0.3184891297, 0.2495567485 ])])
mean = parameters[‘means’]
cov = parameters[‘covariances’]
weights = parameters[‘weights’]
resp = get_responsibilities(data, weights, mean, cov)

Ура, вы получили свои кластерные задания. Я сделал веб-визуализацию результатов некоторых скриншотов.

Давайте поговорим о том, как он устраняет все недостатки k mean:

1. Лучшая стратегия инициализации: поскольку мы используем стратегию инициализации K mean ++, очень маловероятно, что мы получим неправильные кластерные образования.
2. Различная плотность кластера: в отличие от k означает, что он имеет изменяющуюся плотность кластера, что дает нам вероятностную оценку наблюдений с их соответствующими кластерами.
3. Пункты выброса: выбросы принимаются во внимание, поскольку он содержит очень мало ответственности перед ближайшим кластером.

Заключение: Подходя к завершению, мы можем сказать, что любое изображение с основным набором цветов можно кластеризовать с помощью этой техники. Мы можем использовать этот подход для кластеризации любого количества чисел в Image. Здесь я просто взял простой пример с основным цветом, но при необходимости генерация функций может быть изменена с использованием нескольких оттенков цветов. GMM лучше всего подходит для тех случаев, когда вам нужна вероятность назначения изображения определенному кластеру.

Пришло время закрыть мое сообщение в блоге. Если вам действительно интересен этот код, перейдите по моей ссылке на git-hub и не стесняйтесь вносить свой вклад и вносить предложения по улучшениям. Удачного обучения….

Ссылки и справочная информация:

Мартин Крассер Модель латентных переменных, часть 1, модель гауссовской смеси и алгоритм EM (21 ноября 2019 г.): https://krasserm.github.io/2019/11/21/latent-variable-models-part-1/

ethen8181: http://ethen8181.github.io/machine-learning/clustering/GMM/GMM.html

Вики-ссылка: https://en.wikipedia.org/wiki/Expectation–maximization_algorithm

Ресурсы для разработчиков Google: https://developers.google.com/machine-learning/clustering/algorithm/advantages-disadvantages