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

Моя цель - предложение широкого ассортимента товаров и услуг на постоянно высоком качестве обслуживания по самым выгодным ценам.
В этой статье я расскажу как написать свою очень простую машину опорных векторов без scikit-learn или других библиотек с готовой реализацией всего в 30 строчек на Python. Если вам хотелось разобраться в алгоритме SMO, но он показался слишком сложным, то эта статья может быть вам полезна.



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



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

Содержание



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



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



Машина опорных векторов — метод машинного обучения (обучение с учителем) для решения задач классификации, регрессии, детектирования аномалий и т.д. Мы рассмотрим ее на примере задачи бинарной классификации. Наша обучающая выборка — набор векторов фич $\boldsymbol{x}_i$, отнесенных к одному из двух классов $y_i = \pm 1$. Запрос на классификацию — вектор $\boldsymbol{x}$, которому мы должны приписать класс $+1$ или $-1$.

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

Как выбрать лучшую прямую? Интуитивно хочется, чтобы прямая проходила посередине между классами. Для этого записывают уравнение прямой как $\boldsymbol{x} \cdot \boldsymbol{w} + b = 0$ и масштабируют параметры так, чтобы ближайшие к прямой точки датасета удовлетворяли $\boldsymbol{x} \cdot \boldsymbol{w} + b = \pm 1$ (плюс или минус в зависимости от класса) — эти точки и называют опорными векторами.

В таком случае расстояние между граничными точками классов равно $2/|\boldsymbol{w}|$. Очевидно, мы хотим максимизировать эту величину, чтобы как можно более качественно отделить классы. Последнее эквивалентно минимизации $\frac{1}{2} |\boldsymbol{w}|^2$, полностью задача оптимизации записывается

$ \begin{aligned} &\min \frac{1}{2} |\boldsymbol{w}|^2 \\ \text{subject to: } &y_i \left(\boldsymbol{x}_i \cdot \boldsymbol{w} + b\right) - 1 \geq 0. \end{aligned} $

Если её решить, то классификация по запросу $\boldsymbol{x}$ производится так

$ \text{class}(\boldsymbol{x}) = \text{sign}\left(\boldsymbol{x} \cdot \boldsymbol{w} + b\right). $

Это и есть простейшая машина опорных векторов.

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

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

$ \begin{aligned} &\min\left( \frac{1}{2} |\boldsymbol{w}|^2 + C \sum_i \xi_i \right)\\ \text{subject to: }&\xi_i + y_i \left(\boldsymbol{x}_i \cdot \boldsymbol{w} + b\right) - 1 \geq 0,\\ &\xi_i \geq 0, \end{aligned} $

а процедура классификации будет производиться как прежде. Здесь гиперпараметр $C$ отвечает за силу регуляризации, то есть определяет, насколько строго мы требуем от точек соблюдать границу: чем больше $C$ — тем больше $\xi_i$ будет обращаться в ноль и тем меньше точек будут нарушать границу. Опорными векторами в таком случае называют точки, для которых $\xi_i > 0$.

А что если обучающая выборка напоминает логотип группы The Who и точки ни за что нельзя разделить прямой?

Здесь нам поможет остроумная техника — трюк с ядром [4]. Однако, чтобы ее применить, нужно перейти к так называемой двойственной (или дуальной) задаче Лагранжа. Детальное ее описание можно посмотреть в Википедии [5] или в шестой лекции курса [3]. Переменные, в которых решается новая задача, называют дуальными или множителями Лагранжа. Дуальная задача часто проще изначальной и обладает хорошими дополнительными свойствами, например, она вогнута даже если изначальная задача невыпуклая. Хотя ее решение не всегда совпадает с решением изначальной задачи (разрыв двойственности), но есть ряд теорем, которые при определённых условиях гарантируют такое совпадение (сильная двойственность). И это как раз наш случай, так что можно смело перейти к двойственной задаче

$ \begin{aligned} &\max_{\lambda} \sum_{i=1}^n \lambda_i - \frac12 \sum_{i=1}^n \sum_{j=1}^n y_i y_j (\boldsymbol{x}_i \cdot \boldsymbol{x}_j) \lambda_i \lambda_j,\\ \text{subject to: } &0 \leq \lambda_i \leq C, \quad \mbox{ for } i=1, 2, \ldots, n,\\ &\sum_{i=1}^n y_i \lambda_i = 0, \end{aligned} $

где $\lambda_i$ — дуальные переменные. После решения задачи максимизации требуется ещё посчитать параметр $b$, который не вошёл в двойственную задачу, но нужен для классификатора

$ b = \mathbb{E}_{k,\xi_k \neq 0}\left[y_k - \sum_i \lambda_i y_i (\boldsymbol{x}_i \cdot \boldsymbol{x}_k)\right]. $

Классификатор можно (и нужно) переписать в терминах дуальных переменных

$ \text{class}(\boldsymbol{x}) = \text{sign}(f(\boldsymbol{x})),\quad f(\boldsymbol{x}) = \sum_i \lambda_i y_i (\boldsymbol{x}_i \cdot \boldsymbol{x}) + b. $

В чём преимущество этой записи? Обратите внимание, что все векторы из обучающей выборки входят сюда исключительно в виде скалярных произведений $(\boldsymbol{x}_i \cdot \boldsymbol{x}_j)$. Можно сначала отобразить точки на поверхность в пространстве большей размерности, и только затем вычислить скалярное произведение образов в новом пространстве. Зачем это делать видно из рисунка.

При удачном отображении образы точек разделяются гиперплоскостью! На самом деле, всё ещё лучше: отображать-то и не нужно, ведь нас интересует только скалярное произведение, а не конкретные координаты точек. Так что всю процедуру можно эмулировать, заменив скалярное произведение функцией $K(\boldsymbol{x}_i; \boldsymbol{x}_j)$, которую называют ядром. Конечно, быть ядром может не любая функция — должно хотя бы гипотетически существовать отображение $\varphi$, такое что $K(\boldsymbol{x}_i; \boldsymbol{x}_j)=(\varphi(\boldsymbol{x}_i) \cdot \varphi(\boldsymbol{x}_j))$. Необходимые условия определяет теорема Мерсера [6]. В реализации на Python будут представлены линейное ($K(\boldsymbol{x}_i; \boldsymbol{x}_j) = \boldsymbol{x}_i^T \boldsymbol{x}_j$), полиномиальное ($K(\boldsymbol{x}_i; \boldsymbol{x}_j) = (\boldsymbol{x}_i^T \boldsymbol{x}_j)^d$) ядра и ядро радиальных базисных функций ($K(\boldsymbol{x}_i; \boldsymbol{x}_j) = e^{-\gamma |\boldsymbol{x}_i - \boldsymbol{x}_j|^2}$). Как видно из примеров, ядра могут привносить свои специфические гиперпараметры в алгоритм, что тоже будет влиять на его работу.

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


Алгоритм SMO



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

$ \def\M{{\color{red}M}} \def\L{{\color{blue}L}} \begin{aligned} &\max_{\lambda} \sum_{i=1}^n \lambda_i - \frac12 \sum_{i=1}^n \sum_{j=1}^n y_i y_j K(\boldsymbol{x}_i; \boldsymbol{x}_j) \lambda_i \lambda_j,\\ \text{subject to: } &0 \leq \lambda_i \leq C, \quad \mbox{ for } i=1, 2, \ldots, n,\\ &\sum_{i=1}^n y_i \lambda_i = 0, \end{aligned} $

после чего найти параметр

$ b = \mathbb{E}_{k,\xi_k \neq 0}[y_k - \sum_i \lambda_i y_i K(\boldsymbol{x}_i; \boldsymbol{x}_k)], \tag{1} $

а классификатор примет следующий вид

$ \text{class}(\boldsymbol{x}) = \text{sign}(f(\boldsymbol{x})),\quad f(\boldsymbol{x}) = \sum_i \lambda_i y_i K(\boldsymbol{x}_i; \boldsymbol{x}) + b. \tag{2} $

Алгоритм SMO (Sequential minimal optimization, [8]) решения дуальной задачи заключается в следующем. В цикле при помощи сложной эвристики ([9]) выбирается пара дуальных переменных $\lambda_\M$ и $\lambda_\L$, а затем по ним минимизируется целевая функция, с условием постоянства суммы $inline$y_\M\lambda_\M + y_\L\lambda_\L$inline$ и ограничений $0 \leq \lambda_\M \leq C$, $0 \leq \lambda_\L \leq C$ (настройка жёсткости границы). Условие на сумму сохраняет сумму всех $y_i\lambda_i$ неизменной (ведь остальные лямбды мы не трогали). Алгоритм останавливается, когда обнаруживает достаточно хорошее соблюдение так называемых условий ККТ (Каруша-Куна-Такера [10]).

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

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

$ \boldsymbol{K} = \begin{pmatrix} y_1 y_1 K(\boldsymbol{x}_1; \boldsymbol{x}_1) & y_1 y_2 K(\boldsymbol{x}_1; \boldsymbol{x}_2) & \dots & y_1 y_N K(\boldsymbol{x}_1; \boldsymbol{x}_N) \\ y_2 y_1 K(\boldsymbol{x}_2; \boldsymbol{x}_1) & y_2 y_2 K(\boldsymbol{x}_2; \boldsymbol{x}_2) & \dots & y_2 y_N K(\boldsymbol{x}_2; \boldsymbol{x}_N) \\ \cdots & \cdots & \cdots & \cdots \\ y_N y_1 K(\boldsymbol{x}_N; \boldsymbol{x}_1) & y_N y_2 K(\boldsymbol{x}_N; \boldsymbol{x}_2) & \dots & y_N y_N K(\boldsymbol{x}_N; \boldsymbol{x}_N) \\ \end{pmatrix}. \tag{3} $

В дальнейшем я буду использовать обозначение с двумя индексами ($\boldsymbol{K}_{i,j}$), чтобы обратиться к элементу матрицы и с одним индексом ($\boldsymbol{K}_{k}$) для обозначения вектора-столбца матрицы. Дуальные переменные соберём в вектор-столбец $\boldsymbol{\lambda}$. Нас интересует

$ \max_{\boldsymbol{\lambda}} \underbrace{ \sum_{i=1}^n \lambda_i - \frac12 \boldsymbol{\lambda}^T \boldsymbol{K} \boldsymbol{\lambda} }_{\mathscr{L}}. $

Допустим, на текущей итерации мы хотим максимизировать целевую функцию по индексам $\L$ и $\M$. Мы будем брать производные, поэтому удобно выделить слагаемые, содержащие индексы $\L$ и $\M$. Это просто сделать в части с суммой $\lambda_i$, а вот квадратичная форма потребует несколько преобразований.

При расчёте $\boldsymbol{\lambda}^T \boldsymbol{K} \boldsymbol{\lambda}$ суммирование производится по двум индексам, пускай $i$ и $j$. Выделим цветом пары индексов, содержащие $\L$ или $\M$.



Перепишем задачу, объединив всё, что не содержит $\lambda_\L$ или $\lambda_\M$. Чтобы было легче следить за индексами, обратите внимание на $\boldsymbol{K}$ на изображении.

$$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$$

где $\text{const}$ обозначает слагаемые, не зависящие от $\lambda_\L$ или $\lambda_\M$. В последней строке я использовал обозначения

$$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$$

Обратите внимание, что $\boldsymbol{k}_0 + \boldsymbol{Q} \boldsymbol{v}_0$ не зависит ни от $\lambda_\L$, ни от $\lambda_\M$

$$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$$

Ядро — симметрично, поэтому $\boldsymbol{Q}^T = \boldsymbol{Q}$ и можно записать

$ \mathscr{L} = (\boldsymbol{k}_0 + \boldsymbol{Q} \boldsymbol{v}_0 - \boldsymbol{Q} \boldsymbol{v}_0)^T \boldsymbol{v}_0 + \frac{1}{2}\boldsymbol{v}^{\,T}_0 \, \boldsymbol{Q} \, \boldsymbol{v}_0 + \text{const} = (\boldsymbol{k}_0 + \boldsymbol{Q} \boldsymbol{v}_0)^T \boldsymbol{v}_0 - \frac{1}{2} \boldsymbol{v}^{\,T}_0 \, \boldsymbol{Q} \, \boldsymbol{v}_0 + \text{const} $

Мы хотим выполнить максимизацию так, чтобы $inline$y_\L\lambda_\L + y_\M\lambda_\M$inline$ осталось постоянным. Для этого новые значения должны лежать на прямой

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

Несложно убедиться, что для любого $t$

$$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$$

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

$ \mathscr{L}(t) = (\boldsymbol{k}_0 + \boldsymbol{Q} \boldsymbol{v}_0)^T \boldsymbol{v}(t) - \frac12 \boldsymbol{v}^{\,T}(t) \, \boldsymbol{Q} \, \boldsymbol{v}(t) + \text{const}, $

что легко сделать взяв производную

$ \begin{align} \frac{d \mathscr{L}(t)}{d t} = (\boldsymbol{k}_0 + \boldsymbol{Q} \boldsymbol{v}_0)^T \frac{d \boldsymbol{v}}{dt} &- \frac12 \left(\frac{d(\boldsymbol{v}^{\,T} \, \boldsymbol{Q} \, \boldsymbol{v})}{d \boldsymbol{v}}\right)^T \frac{d\boldsymbol{v}}{d t} =\\ &= \boldsymbol{k}_0^T \boldsymbol{u} + \underbrace{\boldsymbol{v}_0^T \boldsymbol{Q}^T \boldsymbol{u} - \boldsymbol{v}^T \boldsymbol{Q}^T \boldsymbol{u}}_{(\boldsymbol{v}_0^T - \boldsymbol{v}^T)\boldsymbol{Q} \boldsymbol{u}} = \boldsymbol{k}_0^T \boldsymbol{u} - t \boldsymbol{u}^T \boldsymbol{Q} \boldsymbol{u}. \end{align} $

Приравнивая производную к нулю, получим

$ t_* = \frac{\boldsymbol{k}^T_0 \boldsymbol{u}}{\boldsymbol{u}^{\,T} \boldsymbol{Q} \boldsymbol{u}}. \tag{5} $

И ещё одно: возможно мы заберёмся дальше чем нужно и окажемся вне квадрата как на картинке. Тогда нужно сделать шаг назад и вернуться на его границу

$$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, классы должны быть указаны как $0$ и $1$ (внутри конвертируются в $-1$ и $+1$, сделано для большей совместимости с sklearn), размерность вектора фич допускается произвольной. Для классификации используется функция predict.


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



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

Сравнение с sklearn.svm.SVC на простом двумерном датасете 
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 ядра)

А это результат работы промышленной версии машины опорных векторов.

Если размерность $2$ кажется слишком маленькой, то можно ещё протестировать на MNIST

Сравнение с sklearn.svm.SVC на 2-х классах из MNIST 
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
Источник: https://habr.com/ru/post/544282/

Вернуться к списку

Интересные статьи

Интересные статьи

Выгрузка пользователей из 1C ЗУП в Битрикс24 или правдивая история о том как настроить интеграцию 1С-Битрикс24 c 1С ЗУП
Выгрузка пользователей из 1C ЗУП в Битрикс24 или правдивая история о том как настроить интеграцию 1С-Битрикс24 с ЗУП без 1С-ника В жизни так бывает, причём бывает чаще чем хотелось б...
Доверяй, но проверяй: контроль неотправленных писем в Битриксе с уведомлением админу
Предыстория Когда-то у меня возникла необходимость проверять наличие неотправленных сообщений в «1С-Битрикс: Управление сайтом» (далее Битрикс) и получать уведомления об этом. Пробле...
1С-Битрикс и попытка его внедрения
Привет, друзья! Меня зовут Петр, я представитель малого белорусского бизнеса со штатом чуть более 20 сотрудников. В данной статье хочу поделиться негативным опытом покупки 1С-Битрикс. ...
Не уникальные описания товаров в интернет-магазине 1С-Битрикс
Не уникальные описания товаров в интернет-магазине 1С-Битрикс
Те, кто собираются открывать интернет-магазин, предварительно начитавшись в интернете о важности уникального контента, о фильтрах, накладываемых поисковиками за копирование материалов с других ресурсо...
Интеграция Zapier и Битрикс24
Интеграция Zapier и Битрикс24
Битрикс24 — популярная в малом бизнесе CRM c большими возможностями даже на бесплатном тарифе. Благодаря API Битрикс24 (даже в облачной редакции) можно легко интегрировать с другими системами.