Svm

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

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

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

Содержание

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

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

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

В простейшем случае классы обучающей выборки можно разделить проведя всего одну прямую как на рисунке (для большего числа фич это была бы гиперплоскость). Теперь, когда придёт запрос на классификацию некоторой точки $\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$ , но также желательно, чтобы таких нарушителей было как можно меньше. Этого можно достичь с помощью модификации целевой функции дополнительным слагаемым (регуляризация $inline$ ):

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

а процедура классификации будет производиться как прежде. Здесь гиперпараметр $inline$ отвечает за силу регуляризации, то есть определяет, насколько строго мы требуем от точек соблюдать границу: чем больше $inline$ тем больше $\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$ дуальные переменные. После решения задачи максимизации требуется ещё посчитать параметр $inline$ , который не вошёл в двойственную задачу, но нужен для классификатора

$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}$ суммирование производится по двум индексам, пускай $inline$ и $inline$ . Выделим цветом пары индексов, содержащие $\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$$

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

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

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

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

Сравнение с sklearn.svm.SVC на простом двумерном датасете

from sklearn.svm import SVCimport matplotlib.pyplot as pltimport seaborn as sns; sns.set()from sklearn.datasets import make_blobs, make_circlesfrom matplotlib.colors import ListedColormapdef 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 ядра)

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

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

Сравнение с sklearn.svm.SVC на 2-х классах из MNIST

from sklearn import datasets, svmfrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import confusion_matriximport matplotlib.pyplot as pltimport seaborn as snsclass_A = 3class_B = 8digits = 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,1X_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 строчках это красиво.

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

https://fbeilstein.github.io/simplest_smo_ever/
страница на http://www.machinelearning.ru
Начала Машинного Обучения, КАУ
https://en.wikipedia.org/wiki/Kernel_method
https://en.wikipedia.org/wiki/Duality_(optimization)
Статья на http://www.machinelearning.ru
https://www.youtube.com/watch?v=MTY1Kje0yLg
https://en.wikipedia.org/wiki/Sequential_minimal_optimization
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).
https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions
http://cs229.stanford.edu/materials/smo.pdf
https://www.researchgate.net/publication/344460740_Yet_more_simple_SMO_algorithm
http://fourier.eng.hmc.edu/e176/lectures/ch9/node9.html
https://github.com/fbeilstein/simplest_smo_ever

Svm

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