Обучение логистической регрессии

Привет, читатель. Сегодня у нас пост про обучение логистической регрессии с L1 и L2 регуляризациями с помощью метода Stochastic Gradient Descent (SGD).

Перед тем как приступить к статье и коду, беглым шагом пробежимся по основным понятиям L1 и L2 регуляризации, логистической регрессии и стахостического градиентного спуска (Stochastic Gradient Descent — SGD).

Итак, самое время поставить перед собой цели, которые к концу статьи должны быть достигнуты:

Реализовать обучение логистической регрессии с L1 и L2 регуляризациями с помощью метода стахостического градиентного спуска;
Для отладки работы алгоритма, реализовать возможность сохранения или вывода ошибки модели, после очередной итерации;
Запустить алгоритм на синтетических данных;
Вывести полученные веса и нарисовать разделяющую границу между классами;
Показать сходимость метода: изобразить графики зависимости значения функции потерь (по всей выборке) после очередной итерации/эпохи (выбрать одно) для разных alpha.

Поехали.

Реализация модели с L1 и L2 регуляризацией с методом SGD

import pandas as pd import numpy as np import matplotlib.pyplot as plt from sklearn.base import BaseEstimator, ClassifierMixin import numpy as np import math %matplotlib inline #Предположим, что в выборке всегда 2 класса class MySGDClassifier(BaseEstimator, ClassifierMixin): def __sigma(self, t): t = np.clip(t, -10, 10) return 1 / (1 + np.exp(-t)) #L1 регуляризация def __l1_penalty(self): return 1/self.C * np.sign(self.theta) #L2 регуляризация def __l2_penalty(self): return 1/self.C * self.theta def __none_penalty(self): return 0 def __init__(self, **kwargs): self.C = kwargs.get('C', 1) #Коэф. регуляризации self.alpha = kwargs.get('alpha', 1) # Скорость спуска self.max_epoch = kwargs.get('max_epoch', 10) #Максимальное количество эпох при обучении модели self.chunk_size = kwargs.get('chunk_size', 5) #Количество элементов в чанке (от 1 до общего числа точек) self.min_err = kwargs.get('min_err', 0.1) #Пороговое значение ошибки penalty = kwargs.get('penalty', 'none') #Способ расчета регуляризации if penalty == 'l1': self.penalty = self.__l1_penalty elif penalty == 'l2': self.penalty = self.__l2_penalty else: self.penalty = self.__none_penalty self.theta = None self.errors = None #Обучение модели def fit(self, X, y=None): import time np.random.seed(int(time.time())) eps = 1e-15 n, m = X.shape #n - число строк, m - столбцов sigma = self.__sigma chunk_size = self.chunk_size errors = np.empty(self.max_epoch) errors[:] = np.nan X_b = np.c_[X, np.ones(n)] self.theta = np.random.randn(m + 1, 1) Y = np.vstack(y) for epoch in range(self.max_epoch): lst = list(range(n)) np.random.shuffle(lst) chunks = [lst[i : i+chunk_size] for i in range(0, n, chunk_size)] for chunk in chunks: xi = X_b[chunk] yi = Y[chunk] gradients = xi.T.dot(sigma(xi.dot(self.theta)) - yi) self.theta = self.theta - self.alpha * (gradients + self.penalty()) p = sigma(X_b.dot(self.theta)) p = np.clip(p, eps, 1-eps) epoch_error = 1 / n * np.sum(-(Y * np.log(p) + (1 - Y)*np.log(1 - p))) errors[epoch] = epoch_error if epoch_error <= self.min_err: break self.errors = errors return self #Возвращение метки класса def predict(self, X): y_hat_proba = self.predict_proba(X) y_hat = np.where(y_hat_proba[0] >= 0.5, 0, 1) return y_hat #Возвращение вероятности каждого из классов def predict_proba(self, X): if self.theta is None: raise Exception("Model is not fitted yet. Use method 'fit'") n, m = X.shape X_b = np.c_[X, np.ones(n)] y1 = self.__sigma(X_b.dot(self.theta)) y0 = 1 - y1 y_hat_proba = np.c_[y0, y1] return y_hat_proba

Запуск алгоритма на синтетических данных и вывод полученных весов с нарисованной разделяющей границей между классами

max_epoch = 20 C1 = np.array([[0., -0.8], [1.5, 0.8]]) C2 = np.array([[1., -0.7], [2., 0.7]]) gauss1 = np.dot(np.random.randn(200, 2) + np.array([5, 3]), C1) gauss2 = np.dot(np.random.randn(200, 2) + np.array([1.5, 0]), C2) X = np.vstack([gauss1, gauss2]) y = np.r_[np.ones(200), np.zeros(200)] fig, ax = plt.subplots() fig.set_size_inches(20,10) ax.scatter(X[:,0], X[:,1], c=y) x_min, x_max = X[:, 0].min(), X[:, 0].max() y_min, y_max = X[:, 1].min(), X[:, 1].max() model = MySGDClassifier(alpha=.01, max_epoch=max_epoch, penalty = 'none', C = 100) model.fit(X, y) print("theta = ", model.theta) t0 = model.theta.item(0) t1 = model.theta.item(1) t2 = model.theta.item(2) x_ = np.array([x_min, x_max]) y_ = -(x_ * t0 + t2) / t1 ax.plot(x_, y_) ax.set_xlim(x_min - 1, x_max + 1) ax.set_ylim(y_min - 1, y_max + 1) plt.show()

Демонстрация сходимости метода с графиками зависимости значения функции потерь после очередной итерации/эпохи для разных alpha

alphas = [0.1, 0.01, 0.001] penalties = ['None', 'l1', 'l2'] penalty_coef = 10 x_min, x_max = X[:, 0].min(), X[:, 0].max() y_min, y_max = X[:, 1].min(), X[:, 1].max() for penalty in penalties: fig, ax = plt.subplots(1, 2, figsize=(20, 7)) ax[1].scatter(X[:,0], X[:,1], c=y) ax[0].set_xlim(0, max_epoch) ax[0].set_ylim(0.2, 1.5) ax[0].set_title('Метод штрафной функции {}'.format(penalty)) for alpha in alphas: model = MySGDClassifier(alpha=alpha, max_epoch=max_epoch, penalty = penalty, C = penalty_coef, min_err = 0.05) model.fit(X, y) ax[0].plot(range(max_epoch), model.errors, label=r'$\alpha$ = {}'.format(alpha)) t0 = model.theta.item(0) t1 = model.theta.item(1) t2 = model.theta.item(2) x_ = np.array([x_min, x_max]) y_ = -(x_ * t0 + t2) / t1 ax[1].plot(x_, y_, label=r'$\alpha$ = {}'.format(alpha)) ax[1].set_xlim(x_min - 1, x_max + 1) ax[1].set_ylim(y_min - 1, y_max + 1) ax[0].legend(loc='upper right') ax[1].legend(loc='lower right') plt.show()

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

Всем знаний!

Подписывайтесь на «Нейрон» в Telegram (@neurondata) ― там свежие статьи и новости из мира науки о данных появляются каждую неделю. Спасибо всем, кто помогает с полезными ссылками и разработкой, особенно Игорю Мариарти и Михаилу Чумашеву.