Легкая реализация многоклассовой модели Метода опорных векторов (SVM) с нуля на языке Python
Легкое создание многоклассовой модели Метода опорных векторов (SVM) с нуля на языке Python
Поставляется с бесплатным глубоким обзором SVM
В этом материале мы реализуем алгоритм машинного обучения метода опорных векторов в общей форме с использованием мягкого края и ядерного подхода. Мы начнем с краткого обзора SVM и его уравнений для обучения и вывода, а затем переведем их в код для разработки модели SVM. Затем мы расширим нашу реализацию для работы с многоклассовыми сценариями и закончим тестированием модели с использованием Sci-kit Learn.
Таким образом, к концу этого материала:
Вы получите четкое представление о различных важных понятиях SVM.
Вы сможете реализовать модель SVM с нуля для бинарных и многоклассовых случаев с глубоким пониманием.
Ночь Ван Гога со Строкой Симметрии и Звездами — Создано автором с использованием DALLE 2
Цель метода SVM – подобрать гиперплоскость, которая максимизировала бы отступ (расстояние от ближайших точек из двух классов). Можно показать, и это интуитивно понятно, что такая гиперплоскость (A) имеет лучшие свойства обобщения и более устойчива к шуму, чем гиперплоскость, которая не максимизирует отступ (B).
Для достижения этого SVM находит гиперплоскость W и b, решая следующую оптимизационную задачу:
Он пытается найти значения W и b, которые максимизируют расстояние до ближайшей точки и правильно классифицируют все (как в ограничении, где y принимает значения ±1). Это можно показать, что эквивалентно следующей оптимизационной задаче:
Для этого можно записать эквивалентную двойственную оптимизационную задачу
Решение этой задачи дает множитель Лагранжа для каждой точки в наборе данных, который мы предполагаем имеет размер m: (α₁, α₂, …, α_N). Целевая функция явно квадратична по α, и ограничения линейны, что означает, что ее можно легко решить с помощью квадратичного программирования. После нахождения решения следует из вывода двойственной задачи:
(xₛ,yₛ) – любая точка с α>0
Обратите внимание, что только точки с α>0 определяют гиперплоскость (вносят вклад в сумму). Они называются опорными векторами.
Таким образом, уравнение для прогнозирования, которое при заданном новом примере x возвращает его прогноз y=±1, имеет следующий вид:
Включает подстановку, а затем алгебраическое упрощение
Эта основная форма SVM называется SVM с жесткими границами, потому что оптимизационная задача, которую она решает (как определено выше), требует, чтобы все точки в обучающем наборе были классифицированы правильно. В практических сценариях может быть некоторый шум, который препятствует или ограничивает существование гиперплоскости, идеально разделяющей данные, в таком случае оптимизационная задача может не вернуть решение или вернуть плохое решение.
Чтобы обобщить 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 можно использовать в таких случаях, но оптимальное решение, скорее всего, будет включать гиперплоскость, которая позволяет гораздо больше ошибок, чем допустимо на практике.
Ядро Soft Margin SVM обобщает Soft Margin SVM для работы с ситуациями, где данные изначально нелинейны. Например, на показанном слева примере не существует линейной гиперплоскости, которую бы Soft Margin SVM могла найти, вне зависимости от установки C, чтобы прилично разделить данные.
Однако, возможно каждую точку x в наборе данных сопоставить с более высокой размерностью с помощью некоторой функции преобразования z=Φ(x), чтобы сделать данные более линейными (или полностью линейными) в новом пространстве более высокой размерности. Это эквивалентно замене x на z в двойственном варианте:
На практике, особенно когда Φ преобразует в очень высокую размерность, вычисление z может занимать очень много времени. Это решается с помощью ядерного трюка, который заменяет zᵗz на эквивалентное вычисление математической функции (называемой ядерной функцией), и которое намного быстрее (например, алгебраическое упрощение zᵗz). Вот некоторые популярные ядерные функции (каждая из которых соответствует некоторому преобразованию Φ в пространство более высокой размерности):
Степень полинома (Q) и RBF γ – гиперпараметры (определяются пользователем)
С этим связана двойная задача оптимизации:
и интуитивно уравнение вывода становится (после алгебраических преобразований):
Полная выводка всех уравнений выше, предполагая, что у вас есть математическое образование, можно найти здесь.
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:
Квадратные скобки говорят вам, что вы можете вызвать это только с параметрами (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 выше).
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.
Надеюсь, вы найдете то, что вы узнали из этой истории, полезным для вашей работы. До скорой встречи, ау-ревуар.