Машина опорных векторов в 30 строчек27.02.2021 23:02

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

Узнав, что на Хабре есть возможность вставлять медиаэлементы, я создал небольшое демо (если вдруг не сработает — можно ещё попытать счастья с версией на гитхабе [1]). Поместите на плоскость (пространство двух фич и ) несколько красных и синих точек (это наш датасет) и машина произведёт классификацию (каждая точка фона закрашивается в зависимости от того, куда был бы классифицирован соответствующий запрос). Подвигайте точки, поменяйте ядро (советую попробовать Radial Basis Functions) и твёрдость границы (константа ). Мои извинения за ужасный код на JS — писал на нём всего несколько раз в жизни, чтобы разобраться в алгоритме используйте код на Python далее в статье.

Содержание

• В следующем разделе я бегло опишу математическую постановку задачи обучения машины опорных векторов, к какой задаче оптимизации она сводится, а также некоторые гиперпараметры, позволяющие регулировать работу алгоритма. Этот материал лишь подводящий к нашей конечной цели и нужен чтобы вспомнить ключевые факты, если такое напоминание не требуется — можете его пропустить без ущерба для понимания дальнейших частей. Если же вы ранее не сталкивались с машинами опорных векторов, то понадобиться куда более полное изложение — обратите внимание на соответствующую лекцию уже ставшего классикой курса Воронцова [2] или на десятую лекцию курса [3], в которую, кстати, входит представленный ниже метод.
• В разделе «Алгоритм SMO» я подробно расскажу как решить поставленную задачу минимизации упрощенным методом SMO и в чём, собственно, состоит упрощение. Будут выкладки, но их объём гораздо меньше, чем в тех подходах к SMO, что доводилось видеть мне.
• Наконец, в разделе «Реализация» будет представлен код классификатора на Python и схема обучения в псевдокоде.
• Узнать насколько алгоритм удался можно в разделе «Сравнение с sklearn.svc.svm» — там приведено визуальное сравнение для небольших датасетов в 2D и confusion matrix для двух классов из MNIST.
• А в «Заключении» что-нибудь да заключим.

Машина опорных векторов

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

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

Как выбрать лучшую прямую? Интуитивно хочется, чтобы прямая проходила посередине между классами. Для этого записывают уравнение прямой как  и масштабируют параметры так, чтобы ближайшие к прямой точки датасета удовлетворяли (плюс или минус в зависимости от класса) — эти точки и называют опорными векторами.

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

А что делать в случае когда точки разных классов взаимно проникают как на рисунке?

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

При удачном отображении образы точек разделяются гиперплоскостью! На самом деле, всё ещё лучше: отображать-то и не нужно, ведь нас интересует только скалярное произведение, а не конкретные координаты точек. Так что всю процедуру можно эмулировать, заменив скалярное произведение функцией , которую называют ядром. Конечно, быть ядром может не любая функция — должно хотя бы гипотетически существовать отображение , такое что . Необходимые условия определяет теорема Мерсера [6]. В реализации на Python будут представлены линейное (), полиномиальное () ядра и ядро радиальных базисных функций (). Как видно из примеров, ядра могут привносить свои специфические гиперпараметры в алгоритм, что тоже будет влиять на его работу.

Возможно, вы видели видео, где действие гравитации поясняют на примере натянутой резиновой пленки в форме воронки [7]. Это работает, так как движение точки по искривленной поверхности в пространстве большой размерности эквивалентно движению её образа в пространстве меньшей размерности, если снабдить его нетривиальной метрикой. Фактически, ядро искривляет пространство.

Алгоритм SMO

Итак, мы у цели, осталось решить дуальную задачу, поставленную в предыдущем разделе

Я собираюсь сделать несколько упрощений.

• Откажусь от сложной эвристики выбора индексов (так сделано в курсе Стэнфордского университета [11]) и буду итерировать по одному индексу, а второй — выбирать случайным образом.
• Откажусь от проверки ККТ и буду выполнять наперёд заданное число итераций.
• В самой процедуре оптимизации, в отличии от классической работы [9] или стэнфордского подхода [11], воспользуюсь векторным уравнением прямой. Это существенно упростит выкладки (сравните объём [12] и [13]).

При расчёте суммирование производится по двум индексам, пускай и . Выделим цветом пары индексов, содержащие или .

Перепишем задачу, объединив всё, что не содержит или . Чтобы было легче следить за индексами, обратите внимание на  на изображении.

$$display$$ \begin{aligned} \mathscr{L} &= \lambda_\M + \lambda_\L — \sum_{j} \lambda_\M \lambda_j K_{\M, j} — \sum_{i} \lambda_\L \lambda_i K_{\L, i} + \text{const} + \\ {\text{компенсация}\atop\text{двойного подсчета}} \rightarrow\qquad &+ \frac{1}{2}\lambda_\M^2 K_{\M,\M} + \lambda_\M \lambda_\L K_{\M,\L} + \frac{1}{2}\lambda_\L^2 K_{\L,\L} = \\ &= \lambda_\M \left (1-\sum_{j} \lambda_j K_{\M, j}\right) + \lambda_\L \left (1-\sum_{i} \lambda_i K_{\L, i}\right)+\\ &+\frac{1}{2}\left (\lambda_\M^2 K_{\M,\M} + 2 \lambda_\M \lambda_\L K_{\M,\L}+\lambda_\L^2 K_{\L,\L} \right) + \text{const} = \\ &=\boldsymbol{k}^T_0 \boldsymbol{v}_0 + \frac{1}{2}\boldsymbol{v}^{\, T}_0 \, \boldsymbol{Q} \, \boldsymbol{v}_0 + \text{const}, \end{aligned} $$display$$

$$display$$ \begin{align} \boldsymbol{v}_0 &= (\lambda_\M, \lambda_\L)^T, \tag{4a}\\ \boldsymbol{k}_0 &= \left (1 — \boldsymbol{\lambda}^T\boldsymbol{K}_{\M}, 1 — \boldsymbol{\lambda}^T\boldsymbol{K}_{\L}\right)^T, \tag{4b}\\ \boldsymbol{Q} &= \begin{pmatrix} K_{\M,\M} & K_{\M,\L} \\ K_{\L,\M} & K_{\L,\L} \\ \end{pmatrix},\tag{4c}\\ \boldsymbol{u} &= (-y_\L, y_\M)^T. \tag{4d} \end{align} $$display$$

$$display$$ \boldsymbol{k}_0 = \begin{pmatrix} 1 — \lambda_\M K_{\M,\M} — \lambda_\L K_{\M,\L} — \sum_{i \neq \M,\L} \lambda_i K_{\M, i}\\ 1 — \lambda_\M K_{\L,\M} — \lambda_\L K_{\L,\L} — \sum_{i \neq \M,\L} \lambda_i K_{\L, i}\\ \end{pmatrix} = \begin{pmatrix} 1 — \sum_{i \neq \M,\L} \lambda_i K_{\M, i}\\ 1 — \sum_{i \neq \M,\L} \lambda_i K_{\L, i}\\ \end{pmatrix} — \boldsymbol{Q} \boldsymbol{v}_0. $$display$$

$$display$$ (\lambda_\M^\text{new}, \lambda_\L^\text{new})^T = \boldsymbol{v}(t) = \boldsymbol{v}_0 + t \boldsymbol{u}. $$display$$

$$display$$ y_\M\lambda_\M^\text{new} + y_\L\lambda_\L^\text{new} = y_\M \lambda_\M + y_\L \lambda_\L + t (-y_\M y_\L + y_\L y_\M) = y_\M\lambda_\M + y_\L\lambda_\L. $$display$$

$$display$$ (\lambda_\M^\text{new}, \lambda_\L^\text{new}) = \boldsymbol{v}_0 + t_*^{\text{restr}} \boldsymbol{u}. $$display$$

Реализация

Принципиальную схему обучения упрощённой машины опорных векторов можно записать как

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

import numpy as np

class SVM:
  def __init__(self, kernel='linear', C=10000.0, max_iter=100000, degree=3, gamma=1):
    self.kernel = {'poly'  : lambda x,y: np.dot(x, y.T)**degree,
         'rbf': lambda x,y: np.exp(-gamma*np.sum((y-x[:,np.newaxis])**2,axis=-1)),
         'linear': lambda x,y: np.dot(x, y.T)}[kernel]
    self.C = C
    self.max_iter = max_iter

  # ограничение параметра t, чтобы новые лямбды не покидали границ квадрата
  def restrict_to_square(self, t, v0, u): 
    t = (np.clip(v0 + t*u, 0, self.C) - v0)[1]/u[1]
    return (np.clip(v0 + t*u, 0, self.C) - v0)[0]/u[0]

  def fit(self, X, y):
    self.X = X.copy()
    # преобразование классов 0,1 в -1,+1; для лучшей совместимости с sklearn
    self.y = y * 2 - 1 
    self.lambdas = np.zeros_like(self.y, dtype=float)
    # формула (3)
    self.K = self.kernel(self.X, self.X) * self.y[:,np.newaxis] * self.y
    
    # выполняем self.max_iter итераций
    for _ in range(self.max_iter):
      # проходим по всем лямбда 
      for idxM in range(len(self.lambdas)):                                    
        # idxL выбираем случайно
        idxL = np.random.randint(0, len(self.lambdas))                         
        # формула (4с)
        Q = self.K[[[idxM, idxM], [idxL, idxL]], [[idxM, idxL], [idxM, idxL]]] 
        # формула (4a)
        v0 = self.lambdas[[idxM, idxL]]                                        
        # формула (4b)
        k0 = 1 - np.sum(self.lambdas * self.K[[idxM, idxL]], axis=1)           
        # формула (4d)
        u = np.array([-self.y[idxL], self.y[idxM]])                            
        # регуляризированная формула (5), регуляризация только для idxM = idxL
        t_max = np.dot(k0, u) / (np.dot(np.dot(Q, u), u) + 1E-15) 
        self.lambdas[[idxM, idxL]] = v0 + u * self.restrict_to_square(t_max, v0, u)
    
    # найти индексы опорных векторов
    idx, = np.nonzero(self.lambdas > 1E-15) 
    # формула (1)
    self.b = np.mean((1.0-np.sum(self.K[idx]*self.lambdas, axis=1))*self.y[idx]) 
  
  def decision_function(self, X):
    return np.sum(self.kernel(X, self.X) * self.y * self.lambdas, axis=1) + self.b

  def predict(self, X): 
    # преобразование классов -1,+1 в 0,1; для лучшей совместимости с sklearn
    return (np.sign(self.decision_function(X)) + 1) // 2

При создании объекта класса SVM можно указать гиперпараметры. Обучение производится вызовом функции fit, классы должны быть указаны как  и  (внутри конвертируются в  и , сделано для большей совместимости с sklearn), размерность вектора фич допускается произвольной. Для классификации используется функция predict.

Сравнение с sklearn.svm.SVC

Не то, чтобы данное сравнение имело особый смысл, ведь речь идёт о крайне упрощённом алгоритме, разработанном исключительно в целях обучения студентов, но всё же. Для тестирования (и чтобы посмотреть как этим всем пользоваться) можно выполнить следующее (этот код тоже есть на гитхаб [14])

from sklearn.svm import SVC
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from sklearn.datasets import make_blobs, make_circles
from matplotlib.colors import ListedColormap

def test_plot(X, y, svm_model, axes, title):
  plt.axes(axes)
  xlim = [np.min(X[:, 0]), np.max(X[:, 0])]
  ylim = [np.min(X[:, 1]), np.max(X[:, 1])]
  xx, yy = np.meshgrid(np.linspace(*xlim, num=700), np.linspace(*ylim, num=700))
  rgb=np.array([[210, 0, 0], [0, 0, 150]])/255.0
  
  svm_model.fit(X, y)
  z_model = svm_model.decision_function(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
  
  plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
  plt.contour(xx, yy, z_model, colors='k', levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
  plt.contourf(xx, yy, np.sign(z_model.reshape(xx.shape)), alpha=0.3, levels=2, cmap=ListedColormap(rgb), zorder=1)
  plt.title(title)

X, y = make_circles(100, factor=.1, noise=.1)
fig, axs = plt.subplots(nrows=1,ncols=2,figsize=(12,4))
test_plot(X, y, SVM(kernel='rbf', C=10, max_iter=60, gamma=1), axs[0], 'OUR ALGORITHM')
test_plot(X, y, SVC(kernel='rbf', C=10, gamma=1), axs[1], 'sklearn.svm.SVC')

X, y = make_blobs(n_samples=50, centers=2, random_state=0, cluster_std=1.4)
fig, axs = plt.subplots(nrows=1,ncols=2,figsize=(12,4))
test_plot(X, y, SVM(kernel='linear', C=10, max_iter=60), axs[0], 'OUR ALGORITHM')
test_plot(X, y, SVC(kernel='linear', C=10), axs[1], 'sklearn.svm.SVC')

fig, axs = plt.subplots(nrows=1,ncols=2,figsize=(12,4))
test_plot(X, y, SVM(kernel='poly', C=5, max_iter=60, degree=3), axs[0], 'OUR ALGORITHM')
test_plot(X, y, SVC(kernel='poly', C=5, degree=3), axs[1], 'sklearn.svm.SVC')

После запуска будут сгенерированы картинки, но так как алгоритм рандомизированный, то они будут слегка отличаться для каждого запуска. Вот пример работы упрощённого алгоритма (слева направо: линейное, полиномиальное и rbf ядра)

from sklearn import datasets, svm
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns

class_A = 3
class_B = 8

digits = datasets.load_digits()
mask = (digits.target == class_A) | (digits.target == class_B)
data = digits.images.reshape((len(digits.images), -1))[mask]
target = digits.target[mask] // max([class_A, class_B]) # rescale to 0,1
X_train, X_test, y_train, y_test = train_test_split(data, target, test_size=0.5, shuffle=True)

def plot_confusion(clf):
  clf.fit(X_train, y_train)
  y_fit = clf.predict(X_test)

  mat = confusion_matrix(y_test, y_fit)
  sns.heatmap(mat.T, square=True, annot=True, fmt='d', cbar=False, xticklabels=[class_A,class_B], yticklabels=[class_A,class_B])
  plt.xlabel('true label')
  plt.ylabel('predicted label');
  plt.show()

print('sklearn:')
plot_confusion(svm.SVC(C=1.0, kernel='rbf', gamma=0.001))
print('custom svm:')
plot_confusion(SVM(kernel='rbf', C=1.0, max_iter=60, gamma=0.001))

Для MNIST я попробовал несколько случайных пар классов (упрощённый алгоритм поддерживает только бинарную классификацию), но разницы в работе упрощённого алгоритма и sklearn не обнаружил. Представление о качестве даёт следующая confusion matrix.

Заключение

Спасибо всем, кто дочитал до конца. В этой статье мы рассмотрели упрощённую учебную реализацию машины опорных векторов. Конечно, она не может состязаться с промышленным образцом, но благодаря крайней простоте и компактному коду на Python, а также тому, что все основные идеи SMO были сохранены, эта версия SVM вполне может занять своё место в учебной аудитории. Стоит отметить, что алгоритм проще не только весьма мудрёного алгоритма [9], но даже его упрощённой версии от Стэнфордского университета [11]. Все-таки машина опорных векторов в 30 строчках — это красиво.

Список литературы

1. https://fbeilstein.github.io/simplest_smo_ever/
2. страница на http://www.machinelearning.ru
3. «Начала Машинного Обучения», КАУ
4. https://en.wikipedia.org/wiki/Kernel_method
5. https://en.wikipedia.org/wiki/Duality_(optimization)
6. Статья на http://www.machinelearning.ru
7. https://www.youtube.com/watch? v=MTY1Kje0yLg
8. https://en.wikipedia.org/wiki/Sequential_minimal_optimization
9. Platt, John. Fast Training of Support Vector Machines using Sequential Minimal Optimization, in Advances in Kernel Methods — Support Vector Learning, B. Scholkopf, C. Burges, A. Smola, eds., MIT Press (1998).
10. https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions
11. http://cs229.stanford.edu/materials/smo.pdf
12. https://www.researchgate.net/publication/344460740_Yet_more_simple_SMO_algorithm
13. http://fourier.eng.hmc.edu/e176/lectures/ch9/node9.html
14. https://github.com/fbeilstein/simplest_smo_ever