Совместимый со Scikit-Learn оценщик KMeans с нуля

В последнее время набирают популярность алгоритмы кластеризации и методы обучения без учителя. Отчасти это связано с экспоненциальным увеличением объема генерируемых данных, а также с тем, что метки для этих данных часто все еще трудно найти. Маркировка данных может занимать много времени и требует человеческих усилий, которые могут быть дорогостоящими. Методы обучения без учителя — это методы машинного обучения, которые можно использовать для извлечения информации из немаркированных данных. Алгоритмы кластеризации специально берут немаркированные точки в наборе данных и пытаются сгруппировать их в «кластеры». Внутри кластеров точки данных очень «похожи» (в некотором смысле, что будет обсуждаться позже), а точки данных между кластерами очень «несходны».

У меня смешанные чувства по поводу кластеризации. Часто бывает трудно количественно оценить, насколько хорошо работает модель, когда у вас нет меры, чтобы определить, что правильно, а что нет. С другой стороны, без помеченных данных зачастую это все, что у нас есть! Несмотря на сложность количественной оценки их эффективности, методы кластеризации могут быть полезны для полу-контролируемого обучения, где имеется небольшой объем размеченных данных и большой объем неразмеченных данных.

В этом посте я расскажу, как написать алгоритм кластеризации k-средних с нуля, используя NumPy. Алгоритм будет объяснен в следующем разделе, и, хотя он кажется простым, его может быть сложно реализовать эффективно! В качестве дополнительного бонуса я расскажу, как реализовать алгоритм таким образом, чтобы он был совместим с Scikit-Learn, чтобы мы могли использовать структуру Scikit-Learn, включая Pipelines и GridSearchCV (что, по общему признанию, не особенно полезно). для этой модели).

Давайте начнем с обсуждения технических деталей алгоритма кластеризации k-средних.

Алгоритм кластеризации k-средних

Алгоритм кластеризации k-средних — это средство разделения набора данных X из n точек и p объектов на k кластеров. Он специально предполагает

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

Основная идея алгоритма состоит в том, чтобы найти центроид для каждого из k кластеров,μ_k, который представляет собой центр кластера. Затем алгоритм назначает точки в наборе данных кластеру, к которому он ближе всего. Это назначение требует, чтобы мы определили метрику d(x_1, x_2), чтобы сообщить алгоритму, насколько близко две точки находятся в нашем пространстве признаков или насколько они похожи. . Чаще всего функция расстояния d принимается за евклидово расстояние.

Алгоритм кластеризации k-средних пытается сформировать кластеры, содержащие точки, сходные по метрике расстояния d. Оказывается, это эквивалентно минимизации дисперсии внутри каждого кластера. Пусть набор всех кластеров равен S = {S_1, …, S_k}, тогда функция стоимости определяется как

Где для каждого x_i у нас есть определение индикаторной функции,

Видно, что сумма,

количество баллов, присвоенное каждому кластеру!

Минимизация функции стоимости является NP-сложной задачей, и для приближенных решений оптимальных центроидов обычно применяется эвристика. Эта эвристика представляет собой жадный итерационный метод, который называется алгоритм Ллойда-Форги и имеет следующие шаги:

1. Инициализируйте центроиды, случайным образом назначив μ_k одной из точек данных.

Потом, пока не сходились,

1. Назначьте каждый x_i кластеру, т. е. найдите r_{i,k}
2. Обновите μ_k, приняв их за среднее значение всех точек данных в кластере,

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

Этот эвристический метод используется для обучения модели k-средних и пытается найти группы данных, которые минимизируют функцию стоимости, которая представляет собой взвешенную сумму квадратов ошибок кластеров. Учитывая, что функция стоимости представляет собой сумму квадратов ошибок, она чувствительна к выбросам, как и линейная регрессия. Кроме того, k-means имеет следующие ограничения:

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

Кроме того, алгоритм чувствителен к начальным условиям, и выбирается ряд различных случайных начальных условий, и в качестве окончательной модели выбирается тот, который дает наилучшие результаты.

Дополнительное примечание о k-means++

k-средних++ — это алгоритм выбора начальных значений (или «начальных значений») для эвристического решения путем указания процедуры инициализации центроидов кластера перед переходом к стандартному итеративному алгоритму k-средних, описанному выше. Алгоритм инициализации состоит из следующих шагов:

• Случайным образом выберите один центр среди точек данных с именем x.
• Для каждой еще не выбранной точки данных X_i вычислите d(x_i, x), расстояние между xи ближайший центр, который уже был выбран.
• Случайным образом выберите одну новую точку данных в качестве нового центра, используя взвешенное распределение вероятностей, где точка x выбирается с вероятностью, пропорциональной d(x_i, x )².
• Повторяйте шаги 2 и 3, пока не будут выбраны k центроидов.
• Теперь, когда начальные центры выбраны, продолжайте использовать стандартную итеративную эвристику.

При инициализации k-средних++ алгоритм гарантированно найдет решение, равное O(log k) и конкурентоспособное с оптимальным решением k-средних.

Написание алгоритма k-средних с помощью NumPy

Сначала давайте начнем с создания простого набора данных в двух измерениях. Мы можем использовать 10 000 точек, принадлежащих 3 кластерам, как показано ниже.

import numpy as np
N = 10000
k = 3

data1 = np.random.randn(N//3,2) + np.array([5,6])
data2 = np.random.randn(N//3,2) + np.array([-5,-6])
data3 = np.random.randn(N//3,2) + np.array([-10,3])
data = np.concatenate((data1, data2, data3))P

Теперь мы давайте нанесем точки,

import matplotlib.pyplot as plt
%matplotlib inline

plt.scatter(data[:,0], data[:,1])

Инициализация центроидов

Чтобы инициализировать центроиды, нам нужно придумать функцию, которая берет набор данных из N точек и выбирает k из них в качестве центроидов. Это наивный процесс инициализации по сравнению с k-means++, но для наших целей он подойдет. Мы определяем эту функцию, используя функцию np.random.choice, чтобы получить записи k разных точек, как показано ниже:

np.random.choice(N, 3,replace=False)
array([1749, 1121, 1310])

Затем мы можем написать функцию для инициализации этих центроидов,

from typing import List

def init_centroids(X: np.array, K: int) -> List[float]:
    N = X.shape[0]
    return X[np.random.choice(N, K,replace=False)]

centroids = init_centroids(data, 3)

Три начальные точки для центроидов:

centroids
array([[-11.23745643,   3.46931906],
       [ -4.91882472,  -4.87376259],
       [  4.78626865,   6.48088   ]])

Мы можем построить их, чтобы увидеть, где они находятся по отношению к остальным точкам,

plt.figure()
plt.scatter(data[:,0], data[:,1])
for c in centroids:
    plt.scatter(*c, color="red")

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

Назначьте каждую точку кластеру

Давайте возьмем пример с первой точкой данных,

point = data[0]
point
array([5.24048539, 6.26693666])

Мы можем напомнить себе, как выглядят центроиды,

centroids
array([[-11.23745643,   3.46931906],
       [ -4.91882472,  -4.87376259],
       [  4.78626865,   6.48088   ]])

Теперь мы хотим найти расстояние от точки до каждого из центроидов и воспользуемся понятием трансляция в NumPy. Мы можем видеть форму point и centroids,

print(f"point.shape = {point.shape}")
print(f"centroids.shape = {centroids.shape}")
point.shape = (2,)
centroids.shape = (3, 2)

Теперь вычитание двух с использованием широковещательной передачи приводит к

point - centroids
array([[16.47794182,  2.79761759],
       [10.15931011, 11.14069925],
       [ 0.45421674, -0.21394334]])

Мы транслировали точку из формы (2,) в (3,2), чтобы она соответствовала форме массива centroids, а затем выполняли поэлементное вычитание.

Теперь мы можем рассчитать расстояние, используя норму из модуля числовой линейной алгебры NumPy. Использование для нормы позволяет нам использовать данные из произвольных измерений!

dists = np.linalg.norm(point - centroids, axis=1)
dists
array([16.71374377, 15.07735924,  0.50208027])

Теперь назначаем точку кластеру, который закрывается, с помощью функции argmin из NumPy.

np.argmin(dists)
2

Это означает, что точка находится ближе всего к кластеру 3 (размер Python имеет индекс 0)! Мы можем видеть это ниже,

plt.figure()
plt.scatter(data[:,0], data[:,1])
for c in centroids:
    plt.scatter(*c, color="red")

plt.scatter(*point, color='purple')

Теперь мы хотим сделать это для каждой точки в наборе данных, мы создаем новый вектор с именем labels, который является кластером, к которому принадлежит каждая точка. Мы используем пустую функцию NumPy, чтобы назначить пустой массив размера N (чтобы память была выделена заранее).

labels = np.empty(data.shape[0])
labels
array([0., 0., 0., ..., 0., 0., 0.])

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

каждая точка в наборе данных,

def update_labels(X: np.array, labels: np.array) -> None:
    for i, point in enumerate(X):
        dists = np.linalg.norm(point - centroids, axis=1) # norm along the rows
        labels[i] = np.argmin(dists)

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

update_labels(data, labels)
labels
array([2., 2., 2., ..., 0., 0., 0.])

Обновление центроидов

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

Мы можем найти все точки, принадлежащие каждому кластеру, используя понятие Маскирование. Чтобы увидеть, как это работает, мы можем найти все помеченные точки, принадлежащие каждому кластеру. Например, мы можем найти, какие точки принадлежат кластеру 0, с помощью следующего.

(labels==0)
array([False, False, False, ...,  True,  True,  True])

Обратите внимание, что это возвращает логический массив, который говорит, является ли каждое значение в массиве 0 или нет. Затем мы можем использовать маскирование, чтобы получить все точки, находящиеся в кластере 0,

data[labels==0]
array([[-10.49346302,   3.01319403],
       [-10.28456084,   2.8061653 ],
       [ -9.34551533,   3.14996967],
       ...,
       [-12.78795698,   2.84267242],
       [-11.20825349,   4.3302826 ],
       [-10.00988014,   4.17066035]])

Мы можем получить значения центроида, взяв среднее значение по каждому столбцу; это можно рассчитать с помощью функции mean из NumPy с axis=0, чтобы показать, что мы суммируем все строки в каждом столбце,

data[labels==0].mean(axis=0)
array([2.72674413, 3.4421267 ])

Написав это как функцию, которая принимает dataframe и labels, мы можем написать цикл for для кластеров, а затем использовать функцию стек для сбора центроидов в виде массива центроидов.

def update_centroids(X: np.array, labels: np.array, K: int) -> None:
    centroids = np.stack([
                        X[labels==i].mean(axis=0) for i in range(K)
    ])
    
    return centroids

И может использовать это для вычисления одной итерации алгоритма k-средних,

update_centroids(data, labels, 3)

plt.figure()
plt.scatter(data[:,0], data[:,1])
for c in centroids:
    plt.scatter(*c, color="red")

Теперь давайте исправим функцию fit, которая будет использовать все функции, которые мы написали до итеративной подгонки под модель k-средних. Последнее, что нам нужно упомянуть, — это сходимость, которая говорит нам, когда наш итерационный метод должен завершиться. Итерация завершается, когда

• Мы достигли некоторого предопределенного максимального количества итераций, называемого max_iter.
• Расстояние между центроидами до и после итерации меньше некоторого предопределенного допуска tol

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

def fit(X: np.array, K: int, max_iters: int = 100, tol:float = 1e-10) -> List[float]:
    centroids = init_centroids(X=X, K=K)
    labels = np.empty(X.shape[0])

    for _ in range(max_iters):

        # label points belonging to clusters
        prev_centroids = centroids

        # update labels
        update_labels(X=X, labels=labels)

        # update centroids
        centroids = update_centroids(X=X, labels=labels, K=K)
        
        if np.linalg.norm(prev_centroids - centroids) < tol:
            break
            
    return centroids

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

centroids = fit(X=data, K=3)

plt.figure()
plt.scatter(data[:,0], data[:,1])
for c in centroids:
    plt.scatter(*c, color="red")

Выглядит довольно хорошо! Теперь давайте поговорим о том, как сделать этот Scikit-Learn совместимым.

Написание оценщика, совместимого со Scikit-Learn

Теперь мы хотим приступить к созданию модели кластеризации k-средних, совместимой со Scikit-Learn, чтобы мы могли использовать встроенные в библиотеку Pipeline и GridSearchCV. Это зависит от нашей способности создать индивидуальный оценщик для Scikit-learn.

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

Одна вещь, которую мы изменили в реализации, заключается в том, что мы позволяем количеству кластеров быть членом класса K, который передается конструктору вместе с максимальным числом итераций max_iter, допуском tol для определения, сходится ли метод, и random_state чтобы увидеть выбор точек данных в качестве начальных центроидов кластера. Другие важные изменения в функциях, написанных выше, в основном либо косметические, либо необходимы для того, чтобы сделать их методами класса. В частности, мы пишем функции как приватные методы класса, и поэтому требуем параметр self и префикс '_' для всех, кроме метода fit, который останется общедоступным. Наконец, мы меняем методы, связанные с центроидами, чтобы они больше не возвращали центроиды, а скорее обновляли член центроидов объектов.

Класс написан ниже,

import numpy as np
from sklearn.base import BaseEstimator, ClusterMixin
from __future__ import annotations

class KMeans(BaseEstimator, ClusterMixin):
    
    def __init__(self, 
                 K: int=2, 
                 max_iter: int=10, 
                 random_state: int = 42,
                 tol: float = 1e-6):
        
        self.K = K
        self.max_iter = max_iter
        self.centroids = None
        self.random_state = random_state
        self.tol = tol
        np.random.seed(random_state)
        
    def _init_centroids(self, X: np.array) -> None:
        N = X.shape[0]
        self.centroids = X[np.random.choice(N,self.K,replace=False)]

    def _update_labels(self, X: np.array, labels: np.array) -> None:
        for i, point in enumerate(X):
            dists = np.linalg.norm(point - self.centroids, axis=1) # sum along the rows
            labels[i] = np.argmin(dists)
            
    def _update_centroids(self, X: np.array, labels: np.array) -> None:
        self.centroids = np.stack([
                            X[labels==i].mean(axis=0) for i in range(self.K)
        ])
    
    def fit(self, X: np.array, y: np.array=None) -> KMeans:
        
        self._init_centroids(X)
        labels = np.empty(X.shape[0])
        
        for _ in range(self.max_iter):
            
            # label points belonging to clusters
            prev_centroids = self.centroids

            # update labels
            self._update_labels(X, labels)
            
            # update centroids
            self._update_centroids(X, labels)
            if np.linalg.norm(prev_centroids - self.centroids) < self.tol:
                break
            
        return self
    
    def predict(self, X: np.array) -> np.array:
        labels = np.empty(X.shape[0])
        self._update_labels(X, labels)
        return labels
    
    def score(self, X: np.array, y: np.array=None) -> np.array:
        if y is None:
            y = self.predict(X)
            
        variance = np.sum([np.linalg.norm(X[y==i] - self.centroids[i], axis=1).sum() 
                       for i in range(self.K)]) / X.shape[0]
        return variance

В дополнение к функции fit у нас также есть общедоступный метод predict, необходимый для классов BaseEstimator и Mixin. Функция предсказания предсказывает, к каким точкам данных кластера принадлежат, путем скрытого вызова метода _update_labels. Тот факт, что мы используем _update_centroids для метода predict, объясняет, почему мы по-прежнему возвращаем массив.

Чтобы наш класс KMeans был совместим с Pipeline и GridSearchCV, нам также нужен метод score для измерения производительности объекта. Для наших целей функция score — это просто сумма квадратов ошибок.

Теперь мы можем создать экземпляр объекта k-средних и подогнать его к набору данных.\

kmeans = KMeans(3)

kmeans = kmeans.fit(data)

Обратите внимание, что функция fit возвращает себя, как это делают оценщики Sciki-learn!

Мы можем видеть центроиды,

kmeans.centroids
array([[-4.99522247, -6.00323426],
       [-9.99131782,  3.00841883],
       [ 4.98813128,  5.98839757]])
plt.figure()
plt.scatter(data[:,0], data[:,1])
for c in kmeans.centroids:
    plt.scatter(*c, color="red")

И предсказать любое количество значений,

kmeans.predict(np.array([10.3242, 5.321]))
array([2., 2.])

Наконец, мы можем найти SSE модели с помощью метода оценки,

kmeans.score(data)
1.2630069857409725

Использование метода локтя и конвейеров

Мы можем попытаться вычислить соответствующее количество кластеров апостериорно, используя метод локтя для различных чисел k.

Идея метода локтя заключается в том, что по мере увеличения количества кластеров сумма внутренней дисперсии или SSE быстро уменьшается, поскольку выборки становятся все более и более однородными. В определенной точке (точка изгиба) кластеры содержат относительно однородные образцы, и уменьшение SSE становится менее значительным по мере увеличения количества кластеров. Считается, что эта точка перегиба возникает из-за того, что мы дольше разделяем законные кластеры, а скорее произвольно разделяем кластеры, и падение SSE должно быть менее значимым.

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

Обратите внимание, что кластеризация чувствительна к масштабированию, поэтому мы используем класс StandardScaler перед применением модели кластеризации, и поэтому нам нужно создать конвейерный объект!

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

pipeline = make_pipeline(StandardScaler(), KMeans(K=k))
model = pipeline.fit(data)

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

K_max = 7
scores = [make_pipeline(StandardScaler(), KMeans(K=k)).fit(data).score(data) 
          for k in range(1,K_max)]

Мы можем построить график зависимости количества кластеров от SSE, чтобы найти оптимальное количество кластеров.

plt.plot(range(1,7), scores)
plt.xlabel("# of clustetrs")
plt.ylabel("SSE", rotation=90)
plt.title("Elbo Curve")

Поскольку 3 кластера так хорошо разделены, колено довольно хорошо определяется на 3.

Давайте возьмем другой пример, где кластеры не так хорошо разделены,

new_data = np.concatenate((
    np.random.randn(N//3,2) + np.array([3,1]),
    np.random.randn(N//3,2) + np.array([-3,0]),
    np.random.randn(N//3,2) + np.array([12,5]),
    np.random.randn(N//3,2) + np.array([8,8]))
)

plt.scatter(new_data[:,0], new_data[:,1])

На этот раз у нас 4 кластера, но две пары несколько перекрываются. Мы можем переобучить модель k-средних для различных значений K и снова построить кривую локтя.

K_max = 10
scores = [make_pipeline(StandardScaler(), KMeans(K=k)).fit(new_data).score(new_data) 
          for k in range(1,K_max)]

plt.plot(range(1,K_max), scores)
plt.xlabel("# of clustetrs")
plt.ylabel("SSE", rotation=90)
plt.title("Elbo Curve")

Оптимальное количество кластеров может быть 2, 3, 4. Очевидно, что разделяющие гиперплоскости не так хорошо определены, как в предыдущем случае!

Наконец, давайте поговорим о том, как использовать GridSearchCV, хотя опять же, это не имеет особого смысла для неконтролируемого обучения, поскольку нет реальной метрики, которую мы пытаемся оптимизировать. Скорее, это просто для того, чтобы показать, как мы создали оценщиков клиентов и как использовать их с остальными инструментами Sciki-learn.

from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import Pipeline

pipeline = Pipeline([
                ('scaler',StandardScaler()), 
                ('model', KMeans())
            ])

Теперь мы можем определить сетку размеров кластеров, а затем подобрать данные,

params = {"model__K": [1,2,3,4,5]}

grid = GridSearchCV(estimator=pipeline,
                    param_grid=params)

results = grid.fit(data)

И просмотрите лучшую подходящую модель,

results.best_estimator_
Pipeline(steps=[('scaler', StandardScaler()), ('model', KMeans(K=1))])

Опять же, результаты не имеют смысла, поскольку поиск по сетке ищет конфигурацию, которая максимизирует SSE, он выбирает 1 кластер. Однако основная идея состоит в том, чтобы показать, как сделать оценщик, совместимый со Scikit-Learn!

Резюме и ссылки

В этом блоге мы рассмотрели, как создать алгоритм кластеризации k-средних, совместимый со Scikit, с использованием NumPy и метода локтя, чтобы попытаться определить, сколько кластеров находится в нашем наборе данных.

Я широко использовал следующие ссылки для написания алгоритма k-средних с нуля.

• https://towardsdatascience.com/k-means-from-scratch-with-numpy-74f79d2b1694
• https://blog.paperspace.com/speed-up-kmeans-numpy-vectorization-broadcasting-profiling/
• https://nicholasvadivelu.com/2021/05/10/fast-k-means/

И собственная документация Scikit-learn для написания Customer Estimator.

Надеюсь, вам понравился этот пост!