Легкая реализация многоклассовой модели Метода опорных векторов (SVM) с нуля на языке Python

Легкое создание многоклассовой модели Метода опорных векторов (SVM) с нуля на языке Python

Поставляется с бесплатным глубоким обзором SVM

В этом материале мы реализуем алгоритм машинного обучения метода опорных векторов в общей форме с использованием мягкого края и ядерного подхода. Мы начнем с краткого обзора SVM и его уравнений для обучения и вывода, а затем переведем их в код для разработки модели SVM. Затем мы расширим нашу реализацию для работы с многоклассовыми сценариями и закончим тестированием модели с использованием Sci-kit Learn.

Таким образом, к концу этого материала:

  • Вы получите четкое представление о различных важных понятиях SVM.
  • Вы сможете реализовать модель SVM с нуля для бинарных и многоклассовых случаев с глубоким пониманием.
Ночь Ван Гога со Строкой Симметрии и Звездами — Создано автором с использованием DALLE 2

Содержание

· Краткий обзорМетод опорных векторов с жестким краемМетод опорных векторов с мягким краемЯдерный метод опорных векторов с мягким краем· РеализацияОсновные импортыОпределение ядер и гиперпараметров SVMОпределение метода прогнозированияОпределение метода прогнозированияТестирование реализацииОбобщение метода Fit на многоклассовый случайОбобщение метода Predict на многоклассовый случайТестирование реализации

Краткий обзор

Метод опорных векторов с жестким краем

Цель метода SVM – подобрать гиперплоскость, которая максимизировала бы отступ (расстояние от ближайших точек из двух классов). Можно показать, и это интуитивно понятно, что такая гиперплоскость (A) имеет лучшие свойства обобщения и более устойчива к шуму, чем гиперплоскость, которая не максимизирует отступ (B).

Рисунок от Ennepetaler86 на Wikimedia. CC BY-SA 3.0 Unported.

Для достижения этого SVM находит гиперплоскость W и b, решая следующую оптимизационную задачу:

Он пытается найти значения W и b, которые максимизируют расстояние до ближайшей точки и правильно классифицируют все (как в ограничении, где y принимает значения ±1). Это можно показать, что эквивалентно следующей оптимизационной задаче:

Для этого можно записать эквивалентную двойственную оптимизационную задачу

Решение этой задачи дает множитель Лагранжа для каждой точки в наборе данных, который мы предполагаем имеет размер m: (α, α₂, …, α_N). Целевая функция явно квадратична по α, и ограничения линейны, что означает, что ее можно легко решить с помощью квадратичного программирования. После нахождения решения следует из вывода двойственной задачи:

(xₛ,yₛ) - любая точка с α>0

Обратите внимание, что только точки с α>0 определяют гиперплоскость (вносят вклад в сумму). Они называются опорными векторами.

Таким образом, уравнение для прогнозирования, которое при заданном новом примере x возвращает его прогноз y=±1, имеет следующий вид:

Включает подстановку, а затем алгебраическое упрощение

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

SVM с мягкой границей

Fit by Soft Margin SVM by Mangat et al on Research Gate. CC BY-SA 4.0 International

Чтобы обобщить SVM с жесткими границами, SVM с мягкой границей адаптирует оптимизационную задачу, вводя константу C (гиперпараметр, задаваемый пользователем), которая контролирует, насколько “жесткой” она должна быть. В частности, она изменяет первоначальную оптимизационную задачу следующим образом:

изменения синим цветом

что позволяет каждой точке совершать некоторое нарушение ϵₙ (например, находиться на неправильной стороне гиперплоскости), но все же стремится в целом уменьшить их путем взвешивания их суммы в целевой функции с помощью C. Приближаясь к бесконечности (обычно до этого), она становится эквивалентной жесткому отступу. Тем временем, меньшее значение C позволит допустить больше нарушений (в обмен на более широкий отступ; то есть, меньший wᵗw).

Довольно удивительно, но можно показать, что эквивалентная двойственная проблема только меняется путем ограничения α для каждой точки ≤C.

Так как допускаются нарушения, опорные векторы (точки с α>0) уже не находятся все на краю отступа. Можно показать, что любой опорный вектор, совершивший нарушение, будет иметь α=C, и что неразделяющие векторы (α=0) не могут совершать нарушения. Мы называем опорные векторы, потенциально совершившие нарушения (α=C), “опорными векторами с неразделяющим отступом”, а остальные (такие, которые не совершали нарушений; находятся на краю) – “опорными векторами с отступом” (0<α<C).

Можно показать, что уравнение вывода не меняется:

Однако, теперь (xₛ,yₛ) должно быть опорным вектором, который не совершил нарушений, так как уравнение предполагает, что он находится на краю отступа (раньше можно было использовать любой опорный вектор).

Ядро Soft Margin SVM

Soft Margin SVM расширяет Hard Margin SVM для обработки шума, но часто данные не разделимы гиперплоскостью из-за факторов, выходящих за пределы шума, таких как естественная нелинейность. Soft Margin SVM можно использовать в таких случаях, но оптимальное решение, скорее всего, будет включать гиперплоскость, которая позволяет гораздо больше ошибок, чем допустимо на практике.

Рисунок: Machine Learner на Wikimedia. CC BY-SA 4.0 International

Ядро Soft Margin SVM обобщает Soft Margin SVM для работы с ситуациями, где данные изначально нелинейны. Например, на показанном слева примере не существует линейной гиперплоскости, которую бы Soft Margin SVM могла найти, вне зависимости от установки C, чтобы прилично разделить данные.

Однако, возможно каждую точку x в наборе данных сопоставить с более высокой размерностью с помощью некоторой функции преобразования z=Φ(x), чтобы сделать данные более линейными (или полностью линейными) в новом пространстве более высокой размерности. Это эквивалентно замене x на z в двойственном варианте:

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

Степень полинома (Q) и RBF γ - гиперпараметры (определяются пользователем)

С этим связана двойная задача оптимизации:

и интуитивно уравнение вывода становится (после алгебраических преобразований):

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

Фото Скотта Грэма на Unsplash

Реализация

Для реализации мы будем использовать

Базовые импорты

Начнем с импорта некоторых базовых библиотек:

import numpy as np                  # для базовых операций над массивами
from scipy.spatial import distance  # для вычисления гауссовского ядра
import cvxopt                       # для решения двойной оптимизационной задачи
import copy                         # для копирования массивов numpy

Определение ядер и гиперпараметров SVM

Начнем с определения трех ядер с использованием соответствующих функций:

Степень полинома (Q) и RBF γ являются гиперпараметрами (определенными пользователем)
class SVM:
    linear = lambda x, xࠤ , c=0: x @ xࠤ.T
    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q
    rbf = lambda x, xࠤ, γ=10: np.exp(-γ*distance.cdist(x, xࠤ,'sqeuclidean'))
    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}

Для согласованности с другими ядрами линейное ядро принимает дополнительный бесполезный гиперпараметр. Как понятно, kernel_funs принимает строку для ядра и возвращает соответствующую функцию ядра.

Теперь продолжим с определением конструктора:

class SVM:
    linear = lambda x, xࠤ , c=0: x @ xࠤ.T
    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q
    rbf = lambda x, xࠤ, γ=10: np.exp(-γ*distance.cdist(x, xࠤ,'sqeuclidean'))
    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}
    
    def __init__(self, kernel='rbf', C=1, k=2):
        # устанавливаем гиперпараметры
        self.kernel_str = kernel
        self.kernel = SVM.kernel_funs[kernel]
        self.C = C                  # параметр регуляризации
        self.k = k                  # параметр ядра
        
        # обучающие данные и опорные векторы (устанавливаются позже)
        self.X, y = None, None
        self.αs = None
        
        # для многоклассовой классификации (устанавливаются позже)
        self.multiclass = False
        self.clfs = []

У SVM есть три основных гиперпараметра: ядро (мы сохраняем данную строку и соответствующую функцию ядра), параметр регуляризации C и гиперпараметр ядра (который передается в функцию ядра); он представляет собой Q для полиномиального ядра и γ для RBF ядра.

Определите метод Fit

Чтобы расширить этот класс с функциями «fit» и «predict» в отдельных ячейках, определим следующую функцию и используем ее в качестве декоратора позже:

SVMClass = lambda func: setattr(SVM, func.__name__, func) or func

Напомним, что подгонка SVM соответствует нахождению опорного вектора α для каждой точки путем решения двойственной оптимизационной задачи:

Пустьα будет переменным столбцовым вектором (α α₂ … α_N)ᵗ, а y будет постоянным столбцовым вектором для меток (y y₂ … y_N)ᵗ, и пусть K будет постоянной матрицей, где K[n, m] вычисляет ядро в (xₙ, xₘ). Напомним следующие эквивалентности на основе индексов для скалярного произведения, внешнего произведения и квадратичной формы соответственно:

чтобы иметь возможность записать двойственную оптимизационную задачу в матричной форме:

Зная, что это квадратичная программа, как мы уже подсказывали, мы можем изучить документацию Quadratic Programming в CVXOPT:

Из документации CVXOPT. GNU General Public License

Квадратные скобки говорят вам, что вы можете вызвать это только с параметрами (P, q), или (P, q, G, h), или (P, q, G, h, A, b) и так далее (любое не указанное значение будет установлено по умолчанию, например, 1).

Чтобы узнать значения для (P, q, G, h, A, b) в нашем случае, мы делаем следующее сравнение:

Упростим сравнение, переписав первое следующим образом:

Обратите внимание, что мы изменили max на min, умножив функцию на -1

Теперь очевидно, что (обратите внимание, что 0≤α эквивалентно -α≤0):

Исходя из этого, мы можем написать следующую функцию fit:

@SVMClassdef fit(self, X, y, eval_train=False):    # если больше двух уникальных меток, вызовите версию для множества классов    if len(np.unique(y)) > 2:        self.multiclass = True        return self.multi_fit(X, y, eval_train)        # если метки заданы в {0,1}, измените их на {-1,1}    if set(np.unique(y)) == {0, 1}: y[y == 0] = -1    # убедитесь, что y - это столбец Nx1 (необходимо для CVXOPT)    self.y = y.reshape(-1, 1).astype(np.double) # Должен быть столбцовым вектором    self.X = X    N = X.shape[0]  # Количество точек        # вычисление ядра для всех возможных пар (x, x') в данных    # с помощью векторизации Numpy это дает матрицу K    self.K = self.kernel(X, X, self.k)        ### Задаем параметры оптимизации    # Для 1/2 x^T P x + q^T x    P = cvxopt.matrix(self.y @ self.y.T * self.K)    q = cvxopt.matrix(-np.ones((N, 1)))        # Для Ax = b    A = cvxopt.matrix(self.y.T)    b = cvxopt.matrix(np.zeros(1))    # Для Gx <= h    G = cvxopt.matrix(np.vstack((-np.identity(N),                                 np.identity(N))))    h = cvxopt.matrix(np.vstack((np.zeros((N,1)),                                 np.ones((N,1)) * self.C)))    # Решение        cvxopt.solvers.options['show_progress'] = False    sol = cvxopt.solvers.qp(P, q, G, h, A, b)    self.αs = np.array(sol["x"])            # наше решение            # Булевый массив, отмечающий точки, которые являются опорными векторами    self.is_sv = ((self.αs-1e-3 > 0)&(self.αs <= self.C)).squeeze()    # индекс некоторого опорного вектора маржи    self.margin_sv = np.argmax((0 < self.αs-1e-3)&(self.αs < self.C-1e-3))        if eval_train:        print(f"Завершена обучение с точностью {self.evaluate(X, y)}")

Мы убеждаемся, что это двоичная задача, а двоичные метки устанавливаются, как предполагает SVM (±1), и что y – это столбцовый вектор размером (N,1). Затем мы решаем оптимизационную задачу для нахождения (α α₂ … α_N)ᵗ.

Мы используем (α α₂ … α_N)ᵗ для получения массива флагов, который равен 1 в любом индексе, соответствующем опорному вектору, чтобы затем применить уравнение прогноза, складывая только опорные векторы и индекс для опорного вектора относительно отступа (xₛ,yₛ). Обратите внимание, что в проверках мы предполагаем, что не-опорные векторы могут не иметь α=0 точно, если это α≤1e-3, то это примерно ноль (мы знаем, что результаты CVXOPT могут быть несколько неточными). Также мы предполагаем, что не-опорные опорные векторы могут не иметь α=C точно.

Определение метода предсказания

Вспомним уравнение прогноза:

@SVMClassdef predict(self, X_t):    if self.multiclass: return self.multi_predict(X_t)    # вычисляем (xₛ, yₛ)    xₛ, yₛ = self.X[self.margin_sv, np.newaxis], self.y[self.margin_sv]    # находим опорные векторы    αs, y, X= self.αs[self.is_sv], self.y[self.is_sv], self.X[self.is_sv]    # вычисляем второй член    b = yₛ - np.sum(αs * y * self.kernel(X, xₛ, self.k), axis=0)    # вычисляем оценку    score = np.sum(αs * y * self.kernel(X, X_t, self.k), axis=0) + b    return np.sign(score).astype(int), score

Вот и все. Мы также можем реализовать метод evaluate для вычисления точности (используется в функции fit выше).

@SVMClassdef evaluate(self, X,y):      outputs, _ = self.predict(X)    accuracy = np.sum(outputs == y) / len(y)    return round(accuracy, 2)

Тестирование реализации

from sklearn.datasets import make_classificationimport numpy as np# Загрузка набора данныхnp.random.seed(1)X, y = make_classification(n_samples=2500, n_features=5,                            n_redundant=0, n_informative=5,                            n_classes=2,  class_sep=0.3)# Тестирование реализованного SVMsvm = SVM(kernel='rbf', k=1)svm.fit(X, y, eval_train=True)y_pred, _ = svm.predict(X)print(f"Точность: {np.sum(y==y_pred)/y.shape[0]}")  #0.9108# Тестирование с помощью Scikitfrom sklearn.svm import SVCclf = SVC(kernel='rbf', C=1, gamma=1)clf.fit(X, y)y_pred = clf.predict(X)print(f"Точность: {sum(y==y_pred)/y.shape[0]}")    #0.9108

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

Обобщение функции fit для многоклассовой классификации

@SVMClassdef multi_fit(self, X, y, eval_train=False):    self.k = len(np.unique(y))      # количество классов    # для каждой пары классов    for i in range(self.k):        # получаем данные для пары        Xs, Ys = X, copy.copy(y)        # изменяем метки на -1 и 1        Ys[Ys!=i], Ys[Ys==i] = -1, +1        # обучаем классификатор        clf = SVM(kernel=self.kernel_str, C=self.C, k=self.k)        clf.fit(Xs, Ys)        # сохраняем классификатор        self.clfs.append(clf)    if eval_train:          print(f"Завершено обучение с точностью {self.evaluate(X, y)}")

Чтобы обобщить модель для многоклассовой классификации с k классами, обучаем бинарный классификатор SVM для каждого класса, где мы проходим по каждому классу и переобозначаем точки, принадлежащие ему, как +1, а точки из всех остальных классов, как -1.

Результат обучения – k классификаторов, когда предоставлено k классов, где i-й классификатор обучен на данных, где i-й класс помечен как +1, а все остальные помечены как -1.

Обобщение predict для многоклассовой классификации

Затем для выполнения прогнозирования для нового примера мы выбираем класс, для которого соответствующий классификатор наиболее уверен (имеет наивысший балл)

@SVMClassdef multi_predict(self, X): # получить прогнозы от всех классификаторов    N = X.shape[0]    preds = np.zeros((N, self.k))    for i, clf in enumerate(self.clfs):        _, preds[:, i] = clf.predict(X)        # получить argmax и соответствующий балл    return np.argmax(preds, axis=1), np.max(preds, axis=1)

Тестирование реализации

from sklearn.datasets import make_classificationimport numpy as np# Загрузить набор данныхnp.random.seed(1)X, y = make_classification(n_samples=500, n_features=2,                            n_redundant=0, n_informative=2,                            n_classes=4, n_clusters_per_class=1,                             class_sep=0.3)# Тестирование SVMsvm = SVM(kernel='rbf', k=4)svm.fit(X, y, eval_train=True)y_pred = svm.predict(X)print(f"Точность: {np.sum(y==y_pred)/y.shape[0]}") # 0.65# Тестирование с использованием Scikitfrom sklearn.multiclass import OneVsRestClassifierfrom sklearn.svm import SVCclf = OneVsRestClassifier(SVC(kernel='rbf', C=1, gamma=4)).fit(X, y)y_pred = clf.predict(X)print(f"Точность: {sum(y==y_pred)/y.shape[0]}")    # 0.65

Построение областей решений для каждого дополнительно приводит к следующему графику:

Иллюстрация автора

Обратите внимание, что, хотя SVM Sci-kit Learn поддерживает OVR по умолчанию (без явного вызова, как показано выше), это потенциально также имеет дополнительные оптимизации, специфичные для SVM.

Полный код

import numpy as np                  # для основных операций с массивамиfrom scipy.spatial import distance  # для вычисления гауссовского ядраimport cvxopt                       # для решения двойственной оптимизационной задачиimport copy                         # для копирования массивов numpy class SVM:    linear = lambda x, xࠤ , c=0: x @ xࠤ .T    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q    rbf = lambda x, xࠤ , γ=10: np.exp(-γ * distance.cdist(x, xࠤ,'sqeuclidean'))    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}    def __init__(self, kernel='rbf', C=1, k=2):        # установить гиперпараметры        self.kernel_str = kernel        self.kernel = SVM.kernel_funs[kernel]        self.C = C                  # параметр регуляризации        self.k = k                  # параметр ядра        # обучающие данные и опорные векторы        self.X, y = None, None        self.αs = None        # для многоклассовой классификации        self.multiclass = False        self.clfs = []      # Это бесполезно здесь (только для блокнота)SVMClass = lambda func: setattr(SVM, func.__name__, func) or func@SVMClassdef fit(self, X, y, eval_train=False):    if len(np.unique(y)) > 2:        self.multiclass = True        return self.multi_fit(X, y, eval_train)    # переименование, если необходимо    if set(np.unique(y)) == {0, 1}: y[y == 0] = -1    # убедиться, что y имеет размерность Nx1    self.y = y.reshape(-1, 1).astype(np.double) # Должен быть столбцовым вектором    self.X = X    N = X.shape[0]    # вычислить ядро для всех возможных пар (x, x') в данных    self.K = self.kernel(X, X, self.k)    # Для 1/2 x^T P x + q^T x    P = cvxopt.matrix(self.y @ self.y.T * self.K)    q = cvxopt.matrix(-np.ones((N, 1)))        # Для Ax = b    A = cvxopt.matrix(self.y.T)    b = cvxopt.matrix(np.zeros(1))    # Для Gx <= h    G = cvxopt.matrix(np.vstack((-np.identity(N),                                 np.identity(N))))    h = cvxopt.matrix(np.vstack((np.zeros((N,1)),                                 np.ones((N,1)) * self.C)))    # Решение        cvxopt.solvers.options['show_progress'] = False    sol = cvxopt.solvers.qp(P, q, G, h, A, b)    self.αs = np.array(sol["x"])            # отображает в опорные векторы    self.is_sv = ((self.αs > 1e-3) & (self.αs <= self.C)).squeeze()    self.margin_sv = np.argmax((1e-3 < self.αs) & (self.αs < self.C - 1e-3))        if eval_train:        print(f"Обучение завершено с точностью {self.evaluate(X, y)}")@SVMClassdef multi_fit(self, X, y, eval_train=False):    self.k = len(np.unique(y))      # количество классов    # для каждой пары классов    for i in range(self.k):        # получить данные для пары        Xs, Ys = X, copy.copy(y)        # изменить метки на -1 и 1        Ys[Ys!=i], Ys[Ys==i] = -1, +1        # обучить классификатор        clf = SVM(kernel=self.kernel_str, C=self.C, k=self.k)        clf.fit(Xs, Ys)        # сохранить классификатор        self.clfs.append(clf)    if eval_train:        print(f"Обучение завершено с точностью {self.evaluate(X, y)}")@SVMClassdef predict(self, X_t):    if self.multiclass: return self.multi_predict(X_t)    xₛ, yₛ = self.X[self.margin_sv, np.newaxis], self.y[self.margin_sv]    αs, y, X= self.αs[self.is_sv], self.y[self.is_sv], self.X[self.is_sv]    b = yₛ - np.sum(αs * y * self.kernel(X, xₛ, self.k), axis=0)    score = np.sum(αs * y * self.kernel(X, X_t, self.k), axis=0) + b    return np.sign(score).astype(int), score@SVMClassdef multi_predict(self, X):    # получить прогнозы от всех классификаторов    preds = np.zeros((X.shape[0], self.k))    for i, clf in enumerate(self.clfs):        _, preds[:, i] = clf.predict(X)    # получить argmax и соответствующий балл    return np.argmax(preds, axis=1)@SVMClassdef evaluate(self, X,y):    outputs, _ = self.predict(X)    accuracy = np.sum(outputs == y) / len(y)    return round(accuracy, 2)from sklearn.datasets import make_classificationimport numpy as np# Загрузить набор данныхnp.random.seed(1)X, y = make_classification(n_samples=500, n_features=2,           n_redundant=0, n_informative=2, n_classes=4,           n_clusters_per_class=1,  class_sep=0.3)# Тестирование SVMsvm = SVM(kernel='rbf
Фото от Нэйтан Ван Эгмонд на Unsplash

В кратком изложении, мы реализовали алгоритм обучения методом опорных векторов (SVM), рассмотрели его общую мягкую границу и ядерную форму. Мы представили обзор SVM, разработали модель в коде, расширили ее для многоклассовых сценариев и проверили правильность нашей реализации с помощью Sci-kit Learn.

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

Ресурсы:

Этот код в основном адаптирован из примера, представленного здесь (MIT License): Перссон, Аладдин. «SVM from Scratch — Machine Learning Python (Support Vector Machine).» YouTube.