jupyter_russian/topic10_boosting/lesson10_part1_boosting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 2\n",
    "Автор материала: программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Тема 10. Бустинг\n",
    "## <center>Часть 1. Бустинг, градиентный бустинг и Xgboost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Бустинг представляет собой жадный алгоритм построения композиции алгоритмов. Основная идея заключается в том, чтобы, имея множество относительно слабых алгоритмов обучения, построить их хорошую линейную комбинацию. Он похож на бэггинг тем, что базовый алгоритм обучения фиксирован. Отличие состоит в том, что обучение базовых алгоритмов для композиции происходит итеративно, и каждый следующий алгоритм стремится компенсировать недостатки композиции всех предыдущих алгоритмов.\n",
    "\n",
    "На примере бустинга стало ясно, что хорошим качеством могут обладать сколь угодно сложные композиции классификаторов, при условии, что они правильно настраиваются. Это развеяло существовавшие долгое время представления о том, что для повышения обобщающей способности необходимо ограничивать сложность алгоритмов. \n",
    "\n",
    "Впоследствии этот феномен бустинга получил теоретическое обоснование. Оказалось, что взвешенное голосование не увеличивает эффективную сложность алгоритма, а лишь сглаживает ответы базовых алгоритмов. Эффективность бустинга объясняется тем, что по мере добавления базовых алгоритмов увеличиваются отступы обучающих объектов. Причём бустинг продолжает раздвигать классы даже после достижения безошибочной классификации обучающей выборки.\n",
    "\n",
    "Общая схема бустинга:\n",
    "- Искомый ансамбль алгоритмов имеет вид $a(x) = \\mbox{sign}(\\sum_{t = 1}^T \\alpha_t b_t(x))$, где $b_t$ - базовые алгоритмы.\n",
    "- Ансамбль строится итеративно, оптимизируя на каждом шаге функционал $Q_t$, равный количеству ошибок текущей композиции на обучающей выборке.\n",
    "- При добавлении слагаемого $\\alpha_t b_t(x)$ в сумму, функционал $Q_t$ оптимизируется только по базовому алгоритму $b_t(x)$ и коэффициенту $\\alpha_t$ при нём, все предыдущие слагаемые считаются фиксированными.\n",
    "- Функционал $Q_t$ имеет вид суммы по объектам обучающей выборки пороговых функций вида $[y_i \\sum_{j = 1}^t \\alpha_j b_j(x_i) < 0]$, имеющих смысл \"текущая композиция ошибается на объекте с номером $i$\". Каждое такое слагаемое имеет вид \"ступеньки\" и является разрывной функцией. Для упрощения решения задачи оптимизации такая пороговая функция заменяется на непрерывно дифференцируемую оценку сверху. В итоге получается новый функционал $\\hat{Q}_t \\geqslant Q_t$, минимизация которого приводит к минимизации исходного функционала $Q_t$.\n",
    "\n",
    "Используя различные аппроксимации для пороговой функции потерь $[z < 0]$, будем получать различные виды бустинга. Примеры:\n",
    "- $e^{-z}$ - AdaBoost\n",
    "- $\\log_2(1 + e^{-z})$ - LogitBoost\n",
    "- $(1 - z)^2$ - GentleBoost\n",
    "- $e^{-cz(z+a)}$ - BrownBoost\n",
    "- другие"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlopen\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import xgboost as xgb\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.datasets import load_boston, load_digits, load_iris\n",
    "from sklearn.ensemble import AdaBoostClassifier, GradientBoostingRegressor\n",
    "from sklearn.metrics import confusion_matrix, mean_squared_error\n",
    "from sklearn.model_selection import GridSearchCV, KFold, train_test_split\n",
    "from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-2, 2, 500)\n",
    "\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.plot(x, x < 0, lw=2, label=\"Threshold function: $[M < 0$]\")\n",
    "plt.plot(x, np.exp(-x), lw=2, label=\"AdaBoost\")\n",
    "plt.plot(x, np.log2(1 + np.exp(-x)), lw=2, label=\"LogitBoost\")\n",
    "plt.plot(x, (1 - x) ** 2, lw=2, label=\"GentleBoost\")\n",
    "plt.plot(x, np.exp(-x * (x + 2)), lw=2, label=\"BrownBoost\")\n",
    "plt.title(\"Various approximations of the threshold function\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Margin\")\n",
    "plt.ylabel(\"Loss\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Алгоритм AdaBoost\n",
    "\n",
    "Как было сказано, алгоритм AdaBoost получается из описанной схемы, при аппроксимации пороговой функции потерь с помощью функции $e^{-z}$. Cуществует теорема (Freund, Schapire, 1996), дающая для достаточно богатых семейств базовых классификаторов явные формулы для базового алгоритма $b_t(x)$ и коэффициента $\\alpha_t$ при нём, на которых достигается минимум функционала $\\hat{Q}_t$. \n",
    "\n",
    "Сам алгоритм выглядит следующим образом:\n",
    "- Инициализировать веса объектов $\\Large w_i^{(0)} = \\frac{1}{l}, i = 1, \\dots, l$.\n",
    "- Для всех $t = 1, \\dots, T$\n",
    "    * Обучить базовый алгоритм $\\Large b_t$, пусть $\\epsilon_t$ – его ошибка на обучающей выборке.\n",
    "    * $\\Large \\alpha_t = \\frac{1}{2}ln\\frac{1 - \\epsilon_t}{\\epsilon_t}$.\n",
    "    * Обновить веса объектов: $\\Large w_i^{(t)} = w_i^{(t-1)} e^{-\\alpha_t y_i b_t(x_i)}, i = 1, \\dots, l$.\n",
    "    * Нормировать веса объектов: $\\Large w_0^{(t)} = \\sum_{j = 1}^k w_j^{(t)}, w_i^{(t)} = \\frac{w_i^{(t)}}{w_0^{(t)}}, i = 1, \\dots, l$.\n",
    "- Вернуть $\\sum_t^{T}\\alpha_tb_t$\n",
    "\n",
    "Таким образом, вновь добавляемый алгоритм обучается путём минимизации взвешенной частоты ошибок на обучающей выборке, а не стандартного функционала, равного частоте ошибок. Вес объекта увеличивается в $e^{\\alpha_t}$ раз, когда $b_t$ допускает на нём ошибку, и уменьшается во столько же раз, когда $b_t$ правильно классифицирует этот объект. Таким образом, непосредственно перед настройкой базового алгоритма наибольший вес накапливается у тех объектов, которые чаще оказывались трудными для классификации предыдущими алгоритмами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Пример для иллюстрации \n",
    "Пусть базовые алгоритмы – всего лишь \"пеньки\", то есть деревья решений глубины 1.\n",
    "<img src='../../img/adaboost_toy_step1.png'>\n",
    "Веса объектов, на которых базовый алгоритм ошибается, увеличиваются (кружки увеличиваются в размере).\n",
    "<img src='../../img/adaboost_toy_step2.png'>\n",
    "В конце базовые алгоритмы \"голосуют\", их веса определялись $\\alpha_t$ в процессе построения.\n",
    "<img src='../../img/adaboost_toy_step3.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../../img/boosting_overfitting.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сначала было обнаружено отстутствие переобучения бустинга вплоть до 1000 базовых классификаторов, позже это было теоретически обосновано. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Градиентный бустинг\n",
    "Метод градиентного бустинга в некотором смысле является обобщением остальных методов бустинга, поскольку он позволяет оптимизировать произвольную дифференцируемую функцию потерь. Данный алгоритм похож на метод градиентного спуска, применяемый для решения задач оптимизации. Основная идея заключается в том, что каждый следующий добавляемый в композицию алгоритм настраивается на остатки предыдущих алгоритмов.\n",
    "\n",
    "Пусть дана функция потерь дифференцируемая $\\Large L(F(x), y)$. Сам алгоритм выглядит следующим образом:\n",
    "- Инициализация композиции константным значением $\\Large F_0(x) = \\arg\\min_{\\alpha} \\sum_{i=1}^n L(\\alpha, y_i)$.\n",
    "- Для всех $\\Large t = 1, \\dots, T$:\n",
    "    * Вычислить остатки предыдущей композиции: $\\Large r_{it} = -[\\nabla_{F(x)} L(F(x_i), y_i)]_{F(x) = F_{t-1}(x)}, i = 1, \\dots, n$.\n",
    "    * Настроить базовый алгоритм $\\Large b_t(x)$ на полученные остатки, т.е. обучить его по выборке $\\Large \\{(x_i, r_{it}), i = 1, \\dots, n\\}$.\n",
    "    * Вычислить коэффициент $\\alpha_t$ перед базовым алгоритмом $\\Large b_t(x)$ как решение следующей одномерной задачи оптимизации:\n",
    "    $\\Large \\alpha_t = \\arg\\min_\\alpha \\sum_{i=1}^n L(F_{t-1}(x_i) + \\alpha b_t(x_i), y_i)$.\n",
    "    * Добавить полученное слагаемое в композицию: $\\Large F_t(x) = F_{t-1}(x) + \\alpha_t b_t(x)$.\n",
    "    \n",
    "Одной из возможных модификаций данного алгоритма является стохастический градиентный бустинг (SGB), который заключается в том, чтобы вычислять суммы вида $\\sum_{i=1}^n$ не по всей обучающей выборке, а только по некоторой её случайной подвыборке. Такой подход является одним из способов регуляризации данного алгоритма и позволяет улучшить качество композиции, сходимость алгоритма и время обучения. \n",
    "\n",
    "Другой способ регуляризации - это введение параметра $\\gamma$, называемого темпом обучения. При добавлении нового слагаемого в композицию, будем добавлять его, умноженное на этот коэффициент. Как правило, чем меньше темп обучения, тем лучше качество итоговой композиции.\n",
    "\n",
    "Для задач регресси обычно использую квадратичную функцию потерь $L(x, y) = (x - y)^2$ или модуль отклонения $L(x, y) = |x - y|$.\n",
    "В задаче классификации используется логистическая функция потерь, которая позволяет возвращать вероятности принадлежности объектов к классам.\n",
    "\n",
    "Одним из наиболее популярных семейств базовых алгоритмов являются решающие деревья. Именно такой вариант градиентного бустинга <a href=\"http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html\">реализован</a> в sklearn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xgboost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проблема многих алгоримтов построения деревьев в том, что в них не уделяется должного внимания регуляризации. \n",
    "В классическом градиентном бустинге применяются такие меры:\n",
    "- ограничение на структуру дерева: максимальная глубина (max_depth), минимальное число объектов в листе (min_samples_leaf)\n",
    "- контролирование темпа обучения (learning_rate)\n",
    "- увеличение \"непохожести\" деревьев за счет рандомизации, как в случайном лесе\n",
    "\n",
    "[Xgboost](https://github.com/dmlc/xgboost) использует еще больше параметров для регуляризации базовых деревьев.\n",
    "\n",
    "Целевая функция для оптимизации в Xgboost состоит из двух слагаемых: специфичной функции потерь и регуляризатора $\\Omega (f_k)$ для каждого из $K$ деревьев, где $f_k$ - прогноз $k$-ого дерева.\n",
    "\n",
    "\n",
    "$$\n",
    "obj(\\theta) = \\sum_{i}^{\\ell} l(y_i - \\hat{y_i}) +  \\sum_{k=1}^{K} \\Omega (f_k)\n",
    "$$\n",
    "\n",
    "Функция потерь зависит от решаемой задачи (Xgboost адаптирован под задачи классификации, регрессии и ранжирования, (подробней хорошо описано в [документации](http://xgboost.readthedocs.io/) Xgboost), а регуляризатор выглядит следующим образом:\n",
    "\n",
    "$$\n",
    "\\Omega(f) = \\gamma T + \\frac{1}{2} \\lambda \\sum_{j=1}^{T}w_j^2\n",
    "$$\n",
    "\n",
    "Первое слагаемое ($\\gamma T$) штрафует модель за большое число листьев $T$, а второе ($\\frac{1}{2} \\lambda \\sum_{j=1}^{T}w_j^2$) контролирует сумму весов модели в листьях. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Примеры\n",
    "\n",
    "В sklearn доступны алгоритмы AdaBoost и GradientBoosting для задач классификации и регрессии.\n",
    "В качестве примера рассмотрим решение задачи восстановления одномерной регрессии с помощью <a href=\"http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html\">GradientBoostingRegressor</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_train = 150\n",
    "n_test = 1000\n",
    "noise = 0.1\n",
    "\n",
    "# Generate data\n",
    "def f(x):\n",
    "    x = x.ravel()\n",
    "\n",
    "    return np.exp(-(x ** 2)) + 1.5 * np.exp(-((x - 2) ** 2))\n",
    "\n",
    "\n",
    "def generate(n_samples, noise):\n",
    "    X = np.random.rand(n_samples) * 10 - 5\n",
    "    X = np.sort(X).ravel()\n",
    "    y = (\n",
    "        np.exp(-(X ** 2))\n",
    "        + 1.5 * np.exp(-((X - 2) ** 2))\n",
    "        + np.random.normal(0.0, noise, n_samples)\n",
    "    )\n",
    "    X = X.reshape((n_samples, 1))\n",
    "\n",
    "    return X, y\n",
    "\n",
    "\n",
    "X_train, y_train = generate(n_samples=n_train, noise=noise)\n",
    "X_test, y_test = generate(n_samples=n_test, noise=noise)\n",
    "# One decision tree regressor\n",
    "dtree = DecisionTreeRegressor(random_state=42)\n",
    "\n",
    "dtree.fit(X_train, y_train)\n",
    "d_predict = dtree.predict(X_test)\n",
    "\n",
    "with plt.xkcd():\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    plt.plot(X_test, f(X_test), \"b\")\n",
    "    plt.scatter(X_train, y_train, c=\"b\", s=20)\n",
    "    plt.plot(X_test, d_predict, \"g\", lw=2)\n",
    "    plt.xlim([-5, 5])\n",
    "    plt.title(\"Decision tree regressor, MSE = %.2f\" % np.sum((y_test - d_predict) ** 2))\n",
    "\n",
    "    gbtree = GradientBoostingRegressor(n_estimators=100, subsample=0.5, random_state=42)\n",
    "    gbtree.fit(X_train, y_train)\n",
    "    gb_predict = gbtree.predict(X_test)\n",
    "\n",
    "    plt.figure(figsize=(10, 6))\n",
    "    plt.plot(X_test, f(X_test), \"b\")\n",
    "    plt.scatter(X_train, y_train, c=\"b\", s=20)\n",
    "    plt.plot(X_test, gb_predict, \"r\", lw=2)\n",
    "    plt.xlim([-5, 5])\n",
    "    plt.title(\n",
    "        \"Boosted ensemble of decision tree regressors, MSE = %.2f\"\n",
    "        % np.sum((y_test - gb_predict) ** 2)\n",
    "    );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим пример использования AdaBoostClassifier с деревьями решений единичной глубины (decision stumps) в качестве базовых алгоритмов для решения задачи классификации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading Pima Indians Diabetes data from UCI Machine learning repository\n",
    "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data\"\n",
    "raw_data = urlopen(url)\n",
    "data = np.loadtxt(raw_data, delimiter=\",\")\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    data[:, :8], data[:, 8], random_state=0\n",
    ")\n",
    "\n",
    "dt = DecisionTreeClassifier(random_state=42).fit(X_train, y_train)\n",
    "\n",
    "# AdaBoost\n",
    "ada = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1), random_state=42).fit(\n",
    "    X_train, y_train\n",
    ")\n",
    "\n",
    "print(\"Decision tree accuracy: %.2f\" % dt.score(X_test, y_test))\n",
    "print(\"AdaBoost accuracy: %.2f\" % ada.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пример использования Xgboost для классификации на данных Iris.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = load_iris()\n",
    "X = iris[\"data\"]\n",
    "y = iris[\"target\"]\n",
    "kf = KFold(n_splits=5, shuffle=True, random_state=13)\n",
    "for train_index, test_index in kf.split(y):\n",
    "    xgb_model = xgb.XGBClassifier().fit(X[train_index], y[train_index])\n",
    "    predictions = xgb_model.predict(X[test_index])\n",
    "    actuals = y[test_index]\n",
    "    print(confusion_matrix(actuals, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пример восстановления регрессии с Xgboost на данных boston.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = load_boston()\n",
    "y = boston[\"target\"]\n",
    "X = boston[\"data\"]\n",
    "kf = KFold(n_splits=5, shuffle=True, random_state=17)\n",
    "for train_index, test_index in kf.split(y):\n",
    "    xgb_model = xgb.XGBRegressor().fit(X[train_index], y[train_index])\n",
    "    predictions = xgb_model.predict(X[test_index])\n",
    "    actuals = y[test_index]\n",
    "    print(mean_squared_error(actuals, predictions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Простой пример подбора параметров с GridSearchCV.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = boston[\"data\"]\n",
    "y = boston[\"target\"]\n",
    "\n",
    "xgb_model = xgb.XGBRegressor()\n",
    "xgb_grid = GridSearchCV(xgb_model, {\"max_depth\": [2, 4, 6]}, verbose=1)\n",
    "xgb_grid.fit(X, y)\n",
    "print(xgb_grid.best_score_)\n",
    "print(xgb_grid.best_params_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ранняя остановка используется для того, чтобы прекратить обучение модели (градиентный спуск), если ошибка за несколько итераций не уменьшилась.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = load_digits()\n",
    "\n",
    "X = digits[\"data\"]\n",
    "y = digits[\"target\"]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=17)\n",
    "clf = xgb.XGBClassifier()\n",
    "clf.fit(\n",
    "    X_train,\n",
    "    y_train,\n",
    "    early_stopping_rounds=10,\n",
    "    eval_metric=\"merror\",\n",
    "    eval_set=[(X_test, y_test)],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errors_by_iter = clf.evals_result()[\"validation_0\"][\"merror\"]\n",
    "plt.plot(range(1, len(errors_by_iter) + 1), errors_by_iter)\n",
    "plt.xlabel(\"iter\")\n",
    "plt.ylabel(\"error\")\n",
    "plt.ylim(0, 0.2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Полезные ссылки\n",
    "- <a href=\"https://en.wikipedia.org/wiki/Boosting_(machine_learning)\">Boosting</a>\n",
    "- [Gradient boosting](https://en.wikipedia.org/wiki/Gradient_boosting)\n",
    "- [Лекция](http://www.machinelearning.ru/wiki/images/c/cd/Voron-ML-Compositions-slides.pdf) К.В. Воронцова по композиционным методам классификации\n",
    "- <a href=\"https://github.com/dmlc/xgboost\">Xgboost</a>\n",
    "- <a href=\"https://github.com/ChenglongChen/Kaggle_CrowdFlower\">Обзор</a> решения победителя соревнования Kaggle \"CrowdFlower\" по предсказанию релевантности выдачи поисковика товаров. Решение на основе Xgboost"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "name": "2_1_7_regul.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}