frobeniusfg Feb 27 2021 at 16:04

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

16 min

9.1K

Python*Mathematics*Machine learning*Studying in IT

В этой статье я расскажу как написать свою очень простую машину опорных векторов без 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$ , а затем по ним минимизируется целевая функция, с условием постоянства суммы $y_\M\lambda_\M + y_\L\lambda_\L$ и ограничений $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}$ на изображении.

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

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

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

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

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

Ядро — симметрично, поэтому $\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}$

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

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

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

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

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

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

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

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

На этом итерация завершается и выбираются новые индексы.

Реализация

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

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

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

Если размерность $inline$ кажется слишком маленькой, то можно ещё протестировать на 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 строчках — это красиво.

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

Tags:

Hubs: