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

  • Tutorial
В этой статье я расскажу как написать свою очень простую машину опорных векторов без 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.

Стоит обратить внимание, что главный потребитель данных из матрицы $\boldsymbol{K}$ — скалярные произведения с $\boldsymbol{\lambda}$ (кроме них на каждой итерации используется ещё только 4 значения для формирования матрицы $\boldsymbol{Q}$). Это означает, что эффективно мы используем только те элементы $\boldsymbol{K}$, что соответствуют ненулевым $\lambda_i$ — остальные будут буквально умножены на ноль. Это важное замечание для больших размеров выборки: если количество переменных прямой задачи соответствовало размерности пространства (числу фич), то для дуальной задачи оно уже равно числу точек (размеру датасета) и квадратичная сложность по памяти может оказаться проблемой. К счастью, лишь небольшое число векторов станут опорными, а значит из всей огромной матрицы $\boldsymbol{K}$ нам понадобятся для расчётов всего несколько элементов, которые можно или пересчитывать каждый раз, или же использовать ленивые вычисления.


Сравнение с 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. https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-98-14.pdf
  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

Средняя зарплата в IT

120 000 ₽/мес.
Средняя зарплата по всем IT-специализациям на основании 6 277 анкет, за 1-ое пол. 2021 года Узнать свою зарплату
Реклама
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее

Комментарии 9

    0
    норм
      –1
      Объясните мне, пожалуйста, в каких случаях может потребоваться разбираться в этом ворохе формул и использовать «упрощенную» версию SVM? Какая практическая польза этого поста?
        +3

        В случае детального изучения машины опорных векторов, конечно же: лучший способ что-то понять — сделать своими руками. А так как промышленная версия SVM весьма сложна, то вполне разумно делать некоторую упрощённую версию, где сохранены только самые важные идеи. Здесь предлагается версия буквально на 30 строк Python — вполне подъёмно даже за полпары сделать.


        К образованию есть разные подходы: от "вызываем эту функцию и чёрная магия делает всё за нас" до "детально разбираемся что и как работает". Я сторонник второго подхода и эта статья написана для студентов, которые тоже придерживаются мнения, что хорошему программисту нужна хорошая база и понимание что и как работает, а не просто лепить одну библиотеку к другой, авось взлетит. Конечно, каждый сам решает, насколько ему нужно фундаментальное образование — можно и ML заниматься "скачал чужую модель — запустил — не подошла — выбросил". Но мне как-то грустно, когда люди по 2 часа не могут развернуть связный список или конвертируют набор битов в число преобразуя его сначала в строку, а то и вообще матрицу на вектор с ошибками множат (всё реальные случаи).


        По поводу "вороха формул". Я так понимаю, имеется в виду раздел "алгоритм SMO". Во-первых, это оригинальный вывод — такого больше нет нигде, ну разве что в моей же заметке на ResearchGate. Так что его в любом случае нужно было привести. Во-вторых, можете посмотреть как это всё выводят другие http://fourier.eng.hmc.edu/e176/lectures/ch9/node9.html и ворох покажется не таким уже и ворохом. Ну и в-третьих, формулы действительно важны. Сравните 30 строк имплементации по формулам в статье и как это имплементируют по стандартным формулам https://github.com/LasseRegin/SVM-w-SMO/blob/master/SVM.py (первый репозитарий, что нагуглился)

        0
        SMO очень элегантный алгоритм с быстрой сходимостью, если пользоваться эвристиками. Раз вы от них уже отказались, то можно вообще написать SVM буквально итерируя две строчки до сходимости (но нужен numpy). www.cs.unm.edu/~ismav/papers/sdm-musik.pdf (5.18)
          0

          Спасибо за ссылку. Не знал об этой работе, обязательно просмотрю детально.


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


          Да и не такая уже и плохая та сходимость, если честно. Вот в приведённой вами статье один из тестов — цифра "2" против всех остальных классов. Признаться, я поленился доставать USPS и в точности повторять всю предобработку, вместо этого просто посмотрел "2" против всех остальных на MNIST, немного модифицировав свой тестовый код из статьи. Ясно, что с алгоритмом из статьи я так сравнивать не имею права, но вот с sklearn можно. При тестовой выборке 899 элементов, sklearn делает от 0 до 2-х ошибок от запуска к запуску (чаще 0). Если считать за итерацию один проход по дуальным переменным (цикл до max_iter у меня в коде), то приведённый в статье алгоритм делает 5 ошибок при 1 итерации, при 3-х итерациях от 0 до 3-х ошибок (чаще 2-3), а при 5 итерациях уже где-то как sklearn. Демо в начале статьи ограничено 60-ю итерациями, что также не сильно много. Я бы не сказал, что это настолько плохая сходимость, что аж срочно надо добавлять эвристики.

            +1
            Я же как раз с вами и соглашаюсь, что SMO это элегантный алгоритм. У Джона Платта вся идея изложена всего на 2х страницах (6 и 7 [9]). Кстати, для полноты картины, можно было бы дать гиперссылку на статью в ссылке [9] выше www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-98-14.pdf Там и псевдокод есть в статье.
            Из интересных особенностей алгоритма — отсутсвие нужды хранить матрицу Грама в памяти, рассчитывая парные ядра по ходу вычисления. Так как SVM это разряженная задача, после первого же прохода большая часть множителей Лагранжа обратится в ноль и пере-рассчитывать ядра почти и не придётся.

            Недопонимание скорее всего возникло, что я не дал пояснения выше, а просто увидел, что вы в коде считаете всю матрицу Грама сразу и подумал, что тогда уж проще сделать код в две строчки с мультипликативными апдейтами (на что и сослался). Кстати тоже простой в имплементации алгоритм не требующий солвера для задачи квадратичного программирования и как таковой имеет дидактическую ценность. А ещё алгоритм Осуны, на который Платт ссылается, тоже довольно поучителен для обучения. Вот рассказываете студентам принцип опорных векторов в пространстве данных, где у вас всего $d$ параметров (по размерности вектора данных), и вдруг переходим в двойственное пространство и мало того, что параметров становится миллионы, если данных много, так еще же это все миллионы в квадрате (задача-то квадратичная). Красиво, конечно, — трюк с ядром теперь можно применять, но студент должен заинтересоваться а стоило ли того. И будет прав — во многих практических случаях не стоило и проще применять линейную SVM. Но если показать студенту алгоритм Осуны и потом его логическое завершение SMO. Потом уже студенты могут сами переходить к fastfood, но это уже совсем другая история ;)
              0

              Спасибо за хорошие замечания! Я поменял ссылку [9] — теперь все ссылки живые. Также очень хорошее замечание о размере матрицы. Я решил добавить про это один абзац в статью (конец раздела об имплементации)


              Стоит обратить внимание, что главный потребитель данных из матрицы K — скалярные произведения с \lambda (кроме них на каждой итерации используется ещё только 4 значения для формирования матрицы Q). Это означает, что эффективно мы используем только те элементы K, что соответствуют ненулевым \lambda_i — остальные будут буквально умножены на ноль. Это важное замечание для больших размеров выборки: если количество переменных прямой задачи соответствовало размерности пространства (числу фич), то для дуальной задачи оно уже равно числу точек (размеру датасета) и квадратичная сложность по памяти может оказаться проблемой. К счастью, лишь небольшое число векторов станут опорными, а значит из всей огромной матрицы K нам понадобятся для расчётов всего несколько элементов, которые можно или пересчитывать каждый раз, или же использовать ленивые вычисления.
          0
          Невозможно объяснить, что такое Матрица машина опорных векторов… Ты должен увидеть это сам.
          Error during service worker registration: s@https://codesandbox.io/static/js/sandbox.743f3294a.js:1:25435

          Невозможно объяснить, почему на скриптах никогда ничего нигде не работает.
            0

            К сожалению, я недостаточно в этом разбираюсь, чтоб исправить (если это проблема скрипта). Всё что я мог — протестировать перед публикацией в доступных мне браузерах: Chrome, Firefox (Windows 10, Ubuntu 14), Edge (Windows 10) и дефолтном браузере на телефоне Honor. Нигде проблем не обнаружил, рассчитывал, что их и не будет.

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

          Самое читаемое