Открытый курс машинного обучения. Тема 10. Градиентный бустинг. Часть 1


    Всем привет! Настало время пополнить наш с вами алгоритмический арсенал.


    Сегодня мы основательно разберем один из наиболее популярных и применяемых на практике алгоритмов машинного обучения — градиентный бустинг. Наша задача — основательно разобраться в бустинге, поэтому статья разбита на 2 части: сегодня мы разберем основную теорию алгоритма, а через 2 недели — практику.


    О том, откуда у бустинга растут корни и что на самом деле творится под капотом алгоритма — в нашем красочном путешествии в мир бустинга под катом. Рванули!



    План этой статьи:


    1. Введение и история бустинга
    2. GBM алгоритм
    3. Функции потерь
    4. Итог про теорию GBM
    5. Домашнее задание №10
    6. Полезные ссылки

    1. Введение и история появления бустинга


    Большинство людей, причастных к анализу данных, хоть раз слышали про бустинг. Этот алгоритм входит в повседневный "джентльменский набор" моделей, которые стоит попробовать в очередной задаче. Xgboost и вовсе часто ассоциируется со стандартным рецептом для победы в ML соревнованиях, породив мем про "стакать xgboost-ы". А еще бустинг является важной частью большинства поисковых систем, иногда выступая еще и их визитной карточкой. Давайте для общего развития посмотрим, как бустинг появился и развивался.


    История появления бустинга


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


    Утвердительный математический ответ на этот вопрос нашелся довольно быстро, что само по себе было важным теоретическим результатом (редкость в ML). Однако, потребовалось несколько лет до появления работоспособных алгоритмов и Adaboost. Их общий подход заключался в жадном построении линейной комбинации простых моделей (базовых алгоритмов) путем перевзвешивания входных данных. Каждая последующая модель (как правило, дерево решений) строилась таким образом, чтобы придавать больший вес и предпочтение ранее некорректно предсказанным наблюдениям.


    Еще немного про Adaboost

    По-хорошему, нам следовало бы последовать примеру большинства остальных курсов по машинному обучению, и перед градиентным бустингом тщательно разобрать его предтечу — Adaboost. Однако, мы решили сразу перейти к самому интересному, так как Adaboost все равно в итоге слился с GBM, когда стало понятно, что это просто его частная вариация.


    Сам алгоритм имеет очень наглядную визуальную интерпретацию стоящей за ним интуиции взвешивания наблюдений. Рассмотрим игрушечный пример задачи классификации, в которой мы будем пробовать на каждой итерации Adaboost разделить данные деревом глубины 1 (так называемым "пнем"). На первых двух итерациях мы увидим следующую картинку:



    Размер точки соответствует полученному ей весу за ошибочное предсказание. И мы видим, как на каждой итерации эти веса растут — пни не могут в одиночку справиться с такой задачей. Однако, когда мы произведем взвешенное голосование ранее построенных пней, мы получим искомое нами разделение:



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



    Adaboost работал хорошо, но из-за того, что обоснований работы алгоритма с его надстройками было мало, вокруг них возник полный спектр спекуляций: кто-то считал его сверх-алгоритмом и волшебной пулей, кто-то был скептичен и разделял мнение, что это малоприменимый подход с жесткой переподгонкой (overfitting). Особенно сильно это касалось применимости на данных с мощными выбросами, к которым Adaboost оказался неустойчив. К счастью, когда за дело взялась профессура Стэнфордской кафедры статистики, уже принесшая миру Lasso, Elastic Net и Random Forest, в 1999 году от Jerome Friedman появилось обобщение наработок алгоритмов бустинга — градиентный бустинг, он же Gradient Boosting (Machine), он же GBM. Этой работой Friedman сразу задал статистическую базу для создания многих алгоритмов, предоставив общий подход бустинга как оптимизации в функциональном пространстве.


    Наследие Стэнфорда

    Вообще, команда Стэнфордской кафедры статистики причастна и к CART, и к bootstrap, и еще ко многим вещам, заранее внеся свои имена в будущие учебники статистики. По большому счету значительная часть нашего повседневного инструментария появилась именно там, и кто знает, что еще появится. Или уже появилось, но еще не нашло достаточного распространения (как, например, glinternet).


    С самим Friedman-ом не так много видеозаписей. Однако, с ним есть очень интересное интервью про создание CART, и вообще про то, как решали статистические задачки (которые мы бы отнесли к анализу данных и data science) 40+ лет назад:



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



    По сути, произошел переход от инженерно-алгоритмических изысканий в построении алгоритмов (так свойственных в ML) к полноценной методологии, как такие алгоритмы строить и изучать. С точки зрения математической начинки, на первый взгляд изменилось не так много: мы всё также добавляем (бустим) слабые алгоритмы, наращивая наш ансамбль постепенными улучшениями тех участков данных, где предыдущие модели "не доработали". Но при построении следующей простой модели, она строится не просто на перевзвешенных наблюдениях, а так, чтобы лучшим образом приближать общий градиент целевой функции. На концептуальном уровне это дало большой простор для фантазии и расширений.


    История становления GBM


    Свое место в "джентльменском наборе" градиентный бустинг занял не сразу — на это потребовалось больше 10 лет с момента появления. Во-первых, у базового GBM появилось много расширений под разные статистические задачи: GLMboost и GAMboost как усиление уже имеющихся GAM моделей, CoxBoost для кривых дожития, RankBoost и LambdaMART для ранжирования. Во-вторых, появилось много реализаций того же GBM под разными названиями и разных платформах: Stochastic GBM, GBDT (Gradient Boosted Decision Trees), GBRT (Gradient Boosted Regression Trees), MART (Multiple Additive Regression Trees), GBM как Generalised Boosting Machines и прочие. К тому же, сообщества machine learner-ов были достаточно разобщены и занимались всем подряд, из-за этого отследить успехи бустинга достаточно сложно.


    В то же время бустинг начали активно применять в задачах ранжирования выдачи поисковых систем. Эту задачу выписали с точки зрения функции потерь, которая штрафует за ошибки в порядке выдачи, так что стало удобно просто вставить ее в GBM. Одними из первых внедрили бустинг в ранжирование AltaVista, а вскоре за ними последовали Yahoo, Yandex, Bing и другие. Причем, говоря о внедрении, речь шла о том, что бустинг на годы вперед становился основным алгоритмом внутри работающих движков, а не еще одной взаимо-заменяемой исследовательской поделкой, живущей в рамках пары научных статей.


    Главную роль в популяризации бустинга сыграли ML соревнования, в особенности kaggle. Исследователям давно не хватало общей платформы, где было бы достаточно участников и задач, чтобы в открытой борьбе за state of the art соревновались люди с их алгоритмами и подходами. Сумрачным немецким гениям, вырастившим у себя в гараже очередной чудо-алгоритм стало нельзя все списывать на закрытые данные, а реальные прорывные библиотеки наоборот получали отличную площадку для развития. Именно это и произошло с бустингом, который прижился на kaggle почти сразу (стоит искать GBM в интервью победителей с 2011 года), а xgboost как библиотека быстро завоевал популярность вскоре после своего появления. При этом xgboost — это не какой-то новый уникальный алгоритм, а просто крайне эффективная реализация классического GBM с некоторыми дополнительными эвристиками.


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


    2. GBM алгоритм


    Постановка ML задачи


    Мы будем решать задачу восстановления функции в общем контексте обучения с учителем. У нас будет набор пар признаков $\large x$ и целевых переменных $\large y$, $\large \left\{ (x_i, y_i) \right\}_{i=1, \ldots,n}$, на котором мы будем восстанавливать зависимость вида $\large y = f(x)$. Восстанавливать будем приближением $\large \hat{f}(x)$, а для понимания, какое приближение лучше, у нас также будет функция потерь $\large L(y,f)$, которую мы будем минимизировать:


    $\large y \approx \hat{f}(x), \\ \large \hat{f}(x) = \underset{f(x)}{\arg\min} \ L(y,f(x))$



    Пока что мы не делаем каких-либо предположений ни о типе зависимости $\large f(x)$, ни о модели наших приближений $\large \hat{f}(x)$, ни о распределении целевой переменной $\large y$. Разве что функция $\large L(y,f)$ должна быть дифференцируемой. Так как задачу надо решать не на всех данных в мире, а только на имеющихся в нашем распоряжении, перепишем все в терминах матожиданий. А именно, будем искать наши приближения $\large \hat{f}(x)$ так, чтобы в среднем минимизировать функцию потерь на тех данных, что есть:


    $\large \hat{f}(x) = \underset{f(x)}{\arg\min} \ \mathbb {E} _{x,y}[L(y,f(x))]$


    К сожалению, функций $\large f(x)$ в мире не просто много — само их функциональное пространство бесконечномерно. Поэтому чтобы хоть как-то решить задачу, в машинном обучении обычно ограничивают пространство поиска каким-нибудь конкретным параметризованным семейством функций $\large f(x, \theta), \theta \in \mathbb{R}^d $. Это сильно упрощает задачу, так как она сводится к уже вполне решаемой оптимизации значений параметров:


    $\large \hat{f}(x) = f(x, \hat{\theta}), \\ \large \hat{\theta} = \underset{\theta}{\arg\min} \ \mathbb {E} _{x,y}[L(y,f(x,\theta))]$


    Аналитические решения для получения оптимальных параметров $\large \hat{\theta}$ в одну строку существуют достаточно редко, поэтому параметры обычно приближают итеративно. Сначала нам надо выписать эмпирическую функцию потерь $\large L_{\theta}(\hat{\theta})$, показывающую, насколько хорошо мы их оценили по имеющимся у нас данных. Также выпишем наше приближение $\large \hat{\theta}$ за $\large M$ итераций в виде суммы (и для наглядности, и чтобы начать привыкать к бустингу):


    $\large \hat{\theta} = \sum_{i = 1}^M \hat{\theta_i}, \\ \large L_{\theta}(\hat{\theta}) = \sum_{i = 1}^N L(y_i,f(x_i, \hat{\theta})) $


    Дело за малым — осталось только взять подходящий итеративный алгоритм, которым мы будем минимизировать $\large L_{\theta}(\hat{\theta})$. Самый простой и часто используемый вариант — градиентный спуск. Для него нужно выписать градиент $\large \nabla L_{\theta}(\hat{\theta})$ и добавлять наши итеративные оценки $\large \hat{\theta_i}$ вдоль него (со знаком минус — мы же хотим уменьшить ошибку, а не нарастить). И все, надо только как-то инициализировать наше первое приближение $\large \hat{\theta_0}$ и выбрать, сколько итераций $\large M$ мы эту процедуру будем продолжать. В нашем неэффективном по памяти виде хранения приближений $\large \hat{\theta}$ наивный алгоритм будет выглядеть следующим образом:


    1. Инициализировать начальное приближение параметров $\large \hat{\theta} = \hat{\theta_0}$
    2. Для каждой итерации $\large t = 1, \dots, M$ повторять:
      1. Посчитать градиент функции потерь $\large \nabla L_{\theta}(\hat{\theta})$ при текущем приближении $\large \hat{\theta}$
        $\large \nabla L_{\theta}(\hat{\theta}) = \left[\frac{\partial L(y, f(x, \theta))}{\partial \theta}\right]_{\theta = \hat{\theta}}$
      2. Задать текущее итеративное приближение $\large \hat{\theta_t}$ на основе посчитанного градиента
        $\large \hat{\theta_t} \leftarrow −\nabla L_{\theta}(\hat{\theta})$
      3. Обновить приближение параметров $\large \hat{\theta}$:
        $\large \hat{\theta} \leftarrow \hat{\theta} + \hat{\theta_t} = \sum_{i = 0}^t \hat{\theta_i} $
    3. Сохранить итоговое приближение $\large \hat{\theta}$
      $\large \hat{\theta} = \sum_{i = 0}^M \hat{\theta_i} $
    4. Пользоваться найденной функцией $\large \hat{f}(x) = f(x, \hat{\theta})$ по назначению


    Функциональный градиентный спуск


    Расширим сознание: представим на секунду, что мы можем проводить оптимизацию в функциональном пространстве и итеративно искать приближения $\large \hat{f}(x)$ в виде самих функций. Выпишем наше приближение в виде суммы инкрементальных улучшений, каждое из которых является функцией. Для удобства сразу будем считать эту сумму, начиная с начального приближения $\large \hat{f_0}(x)$:


    $\large \hat{f}(x) = \sum_{i = 0}^M \hat{f_i}(x)$


    Магии пока не случилось, мы просто решили, что будем искать наше приближение $\large \hat{f}(x)$ не в виде одной большой модели с кучей параметров (как, например, нейросеть), а в виде суммы функций, делая вид, что таким образом мы двигаемся в функциональном пространстве.


    Чтобы решить задачу, нам все равно придется ограничить свой поиск каким-то семейством функций $\large \hat{f}(x) = h(x, \theta)$. Но, во-первых, сумма моделей может быть сложнее чем любая модель из этого семейства (сумму двух деревьев-пней глубины 1 уже не приблизить одним пнем). Во-вторых, общая задача все еще происходит в функциональном пространстве. Сразу учтем, на каждом шаге для функций нам понадобится подбирать оптимальный коэффициент $\large \rho \in \mathbb{R}$. Для шага $\large t$ задача выглядит следующим образом:


    $\large \hat{f}(x) = \sum_{i = 0}^{t-1} \hat{f_i}(x), \\ \large (\rho_t,\theta_t) = \underset{\rho,\theta}{\arg\min} \ \mathbb {E} _{x,y}[L(y,\hat{f}(x) + \rho \cdot h(x, \theta))], \\ \large \hat{f_t}(x) = \rho_t \cdot h(x, \theta_t)$


    А вот теперь время магии. Мы выписывали все наши задачи в общем виде, словно мы можем просто так брать и обучать какие угодно модели $\large h(x, \theta)$ относительно каких угодно функций потерь $\large L(y, f(x, \theta))$. На практике это крайне сложно, поэтому был придуман простой способ свести задачу к чему-то решаемому.


    Зная выражение градиента функции потерь, мы можем посчитать его значения на наших данных. Так давайте обучать модели так, чтобы наши предсказания были наиболее скоррелированными с этим градиентом (со знаком минус). То есть будем решать задачу МНК-регрессии, пытаясь выправлять предсказания по этим остаткам. И для классификации, и для регрессии, и для ранжирования под капотом мы все время будем минимизировать квадрат разности между псевдо-остатками $\large r$ и нашими предсказаниями. Для шага $\large t$ итоговая задача выглядит следующим образом:


    $ \large \hat{f}(x) = \sum_{i = 0}^{t-1} \hat{f_i}(x), \\ \large r_{it} = -\left[\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x)=\hat{f}(x)}, \quad \mbox{for } i=1,\ldots,n ,\\ \large \theta_t = \underset{\theta}{\arg\min} \ \sum_{i = 1}^{n} (r_{it} - h(x_i, \theta))^2, \\ \large \rho_t = \underset{\rho}{\arg\min} \ \sum_{i = 1}^{n} L(y_i, \hat{f}(x_i) + \rho \cdot h(x_i, \theta_t))$



    Классический GBM алгоритм Friedman-а


    Теперь у нас есть все необходимое, чтобы наконец выписать GBM алгоритм, предложенный Jerome Friedman в 1999 году. Мы все так же решаем общую задачу обучения с учителем. На вход алгоритма нужно собрать несколько составляющих:


    • набор данных $\large \left\{ (x_i, y_i) \right\}_{i=1, \ldots,n}$;
    • число итераций $\large M$;
    • выбор функции потерь $\large L(y, f)$ с выписанным градиентом;
    • выбор семейства функций базовых алгоритмов $\large h(x, \theta)$, с процедурой их обучения;
    • дополнительные гиперпараметры $\large h(x, \theta)$, например, глубина дерева у деревьев решений;

    Единственный момент, который остался без внимания — начальное приближение $\large f_0(x)$. Для простоты, в качестве инициализации используют просто константное значение $\large \gamma$. Его, а также оптимальный коэффициент $\large \rho $ находят бинарным поиском, или другим line search алгоритмом относительно исходной функции потерь (а не градиента). Итак, GBM алгоритм:


    1. Инициализировать GBM константным значением $\large \hat{f}(x) = \hat{f}_0, \hat{f}_0 = \gamma, \gamma \in \mathbb{R}$
      $\large \hat{f}_0 = \underset{\gamma}{\arg\min} \ \sum_{i = 1}^{n} L(y_i, \gamma)$
    2. Для каждой итерации $\large t = 1, \dots, M$ повторять:
      1. Посчитать псевдо-остатки $\large r_t$
        $\large r_{it} = -\left[\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x)=\hat{f}(x)}, \quad \mbox{for } i=1,\ldots,n$
      2. Построить новый базовый алгоритм $\large h_t(x)$ как регрессию на псевдо-остатках $\large \left\{ (x_i, r_{it}) \right\}_{i=1, \ldots,n}$
      3. Найти оптимальный коэффициент $\large \rho_t $ при $\large h_t(x)$ относительно исходной функции потерь
        $\large \rho_t = \underset{\rho}{\arg\min} \ \sum_{i = 1}^{n} L(y_i, \hat{f}(x_i) + \rho \cdot h(x_i, \theta))$
      4. Сохранить $\large \hat{f_t}(x) = \rho_t \cdot h_t(x)$
      5. Обновить текущее приближение $\large \hat{f}(x)$
        $\large \hat{f}(x) \leftarrow \hat{f}(x) + \hat{f_t}(x) = \sum_{i = 0}^{t} \hat{f_i}(x)$
    3. Скомпоновать итоговую GBM модель $\large \hat{f}(x)$
      $\large \hat{f}(x) = \sum_{i = 0}^M \hat{f_i}(x) $
    4. Покорять kaggle и мир обученной моделью (предсказания там делать, сами разберетесь)

    Пошаговый пример работы GBM


    Попробуем на игрушечном примере разобраться, как работает GBM. Будем с его помощью восстанавливать зашумленную функцию $\large y = cos(x) + \epsilon, \epsilon \sim \mathcal{N}(0, \frac{1}{5}), x \in [-5,5]$.



    Это задача регрессии на вещественную целевую переменную, поэтому воспользуемся средне-квадратичной функцией потерь. Самих пар наблюдений мы сгенерируем 300 штук, а приближать будем деревьями решений глубины 2. Соберем вместе все, что нам нужно для применения GBM:


    • Игрушечные данные $\large \left\{ (x_i, y_i) \right\}_{i=1, \ldots,300}$
    • Число итераций $\large M = 3$ ✓;
    • Среднеквадратичная функция потерь $\large L(y, f) = (y-f)^2$
      Градиент $\large L(y, f) = L_2$ loss это просто остатки $\large r = (y - f)$ ✓;
    • Деревья решений в качестве базовых алгоритмов $\large h(x)$ ✓;
    • Гиперпараметры деревьев решений: глубина деревьев равна 2 ✓;

    У среднеквадратичной ошибки все просто и с инициализацией $\large \gamma$ и с коэффициентами $\large \rho_t$. А именно, инициализировать GBM мы будем средним значением $\large \gamma = \frac{1}{n} \cdot \sum_{i = 1}^n y_i$, а все $\large \rho_t$ равны 1.


    Запустим GBM и будем рисовать два типа графиков: актуальное приближение $\large \hat{f}(x)$ (синий график), а также каждое построенное дерево $\large \hat{f_t}(x)$ на своих псевдо-остатках (зеленый график). Номер графика соответствует номеру итерации:



    Заметим, что ко второй итерации наши деревья повторили основную форму функции. Однако, на первой итерации мы видим, что алгоритм построил только "левую ветвь" функции ($\large x \in [-5, -4]$). Так получилось банально потому, что нашим деревьям не хватило глубины, чтобы построить симметричную ветку сразу, а ошибка на левой ветви была больше. Так что, правая ветвь "проросла" только на второй итерации.


    В остальном процесс прошел так как мы и ожидали: на каждом шаге наши псевдо-остатки уменьшались, а GBM все ближе приближал исходный косинус. Однако, деревья по построению не могут приблизить непрерывную функцию, поэтому на этом примере GBM полезен, но не идеален. Чтобы самим поиграться с тем, как GBM приближает функции, в блоге Brilliantly wrong есть офигенная интерактивная демка:


    http://arogozhnikov.github.io/2016/06/24/gradient_boosting_explained.html

    3. Функции потерь


    Что делать, если мы хотим решать не обычную среднеквадратичную регрессию, а, скажем, задачу бинарной классификации? Нет проблем, надо только выбрать соответствующую задаче и целевой переменной $\large y$ функцию потерь $\large L(y, f)$. Это самый важный верхнеуровневый момент, определяющий что именно мы будем оптимизировать, и какие свойства ожидать от нашей итоговой модели.


    Как правило, самим нам ничего придумывать и выписывать не надо — исследователи уже все сделали за нас. Сегодня мы разберем функции потерь двух самых часто встречающихся задач: регрессии $\large y \in \mathbb{R}$ и бинарной классификации $\large y \in \left\{-1, 1\right\}$. Что делать с многоклассовой классификаций, рангами, а также всякими промежуточными случаями вроде целочисленной регрессии, расскажем в другой раз.


    Функции потерь регрессии


    Сначала разберемся с регрессией $\large y \in \mathbb{R}$. Выбирая функцию потерь в этом случае, мы прежде всего решаем, какое именно свойство условного распределения $\large (y|x)$ мы хотим восстановить. Наиболее частые варианты:


    • $\large L(y, f) = (y - f)^2$, оно же $\large L_2$ loss, оно же Gaussian loss. Это классическое условное среднее, самый частый и простой вариант. Если нет никакой дополнительной информации или требований к устойчивости (робастности) модели — используйте его.
    • $\large L(y, f) = |y - f|$, оно же $\large L_1$ loss, оно же Laplacian loss. Эта, на первый взгляд, не очень дифференцируемая вещь, на самом деле определяет условную медиану. Медиана, как мы знаем, более устойчива к выбросам, поэтому в некоторых задачах эта функция потерь предпочтительнее, так как она не так сильно штрафует большие отклонения, нежели квадратичная функция.
    • $inline$ \large \begin{equation} L(y, f) =\left\{ \begin{array}{@{}ll@{}} (1 - \alpha) \cdot |y - f|, & \text{if}\ |y-f| \leq 0 \\ \alpha \cdot |y - f|, & \text{if}\ |y-f| >0 \end{array}\right. \end{equation}, \alpha \in (0,1) $inline$, оно же $\large L_q$ loss, оно же Quantile loss. Если бы мы, допустим, захотели не условную медиану с $\large L_1$, а условную 75%-квантиль, мы бы воспользовались этим вариантом с $\large \alpha = 0.75$. Можно видеть, что эта функция ассиметрична и больше штрафует наблюдения, оказывающиеся по нужную нам сторону квантили.


    Давайте попробуем воспользоваться $\large L_q$ функцией потерь на наших игрушечных данных, пытаясь восстановить условную 75%-квантиль косинуса. Соберем все воедино:


    • Игрушечные данные $\large \left\{ (x_i, y_i) \right\}_{i=1, \ldots,300}$
    • Число итераций $\large M = 3$ ✓;
    • Квантильная функция потерь $inline$ \large \begin{equation} L_{0.75}(y, f) =\left\{ \begin{array}{@{}ll@{}} 0.25 \cdot |y - f|, & \text{if}\ |y-f| \leq 0 \\ 0.75 \cdot |y - f|, & \text{if}\ |y-f| >0 \end{array}\right. \end{equation} $inline$ ✓;
    • Градиент $\large L_{0.75}(y, f)$ — индикаторная функция, взвешенная нашим $\large \alpha = 0.75$. Мы будем обучать деревья чему-то, очень похожему на классификатор:
      $\large r_{i} = -\left[\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x)=\hat{f}(x)} = $
      $\large = \alpha I(y_i > \hat{f}(x_i) ) - (1 - \alpha)I(y_i \leq \hat{f}(x_i) ), \quad \mbox{for } i=1,\ldots,300$ ✓;
    • Деревья решений в качестве базовых алгоритмов $\large h(x)$ ✓;
    • Гиперпараметры деревьев решений: глубина деревьев равна 2 ✓;

    У нас есть очевидное начальное приближение — просто взять нужную нам квантиль $\large y$. Однако, про оптимальные коэффициенты $\large \rho_t$ нам ничего не известно, так что воспользуемся стандартным line search. Посмотрим, что у нас получилось:



    Непривычно видеть, что по факту мы обучаем что-то очень непохожее на обычные остатки — на каждой итерации $\large r_{i} $ принимают только два возможных значения. Однако, результат работы GBM достаточно похож на нашу исходную функцию.


    Если оставить алгоритм обучаться на этом игрушечном примере, мы получим почти такой же результат, что и с квадратичной функцией потерь, смещенный на $\large \approx 0.135$. Но если бы мы искали квантили выше 90%, могли бы возникнуть вычислительные трудности. А именно, если соотношение числа точек выше нужной квантили будет слишком мало (как несбалансированные классы), модель не сможет качественно обучиться. Про такие нюансы стоит задумываться, решая нетипичные задачи.


    Еще чуть-чуть про регрессионные функции потерь

    Для задачи регрессии разработано достаточно много функций потерь, в том числе с дополнительными свойствами робастности. Один такой пример — функция потерь Губера, она же Huber loss. Суть функции в том, что на небольших отклонениях она работает как $\large L_2$, а с заранее заданного порога, начинает работать как $\large L_1$. Это позволяет уменьшить вклад выбросов и следующих за ними квадратично-больших ошибок на общий вид функции, при этом не акцентируя внимание на мелких неточностях и отклонениях.


    Можно посмотреть, как работает эта функция потерь на следующем игрушечном примере. За основу возьмем игрушечные данные функции $\large y = \frac{sin(x)}{x}$, к которым был добавлен специальный шум: смесь из Гауссовского распределения и распределения Бернулли, выступающего в роли одностороннего генератора выбросов. Сами функции потерь приведены на графиках A-D, а соответствующие им GBM — на графиках F-H (на графике E — исходная функция):


    И в крупном разрешении.

    В этом примере в качестве базовых алгоритмов для визуальной наглядности были использованы сплайны. Мы ведь уже говорили, что бустить можно не только деревья? О том, какие базовые алгоритмы чаще используются в бустинге и какие у них свойства, мы расскажем в следующей статье про GBM (часть 2). Мы все равно находимся в бонусной секции курса, куда заглянут только самые стойкие, так что легкие спойлеры позволительны.


    По результатам примера, из-за искусственно созданной проблемы с шумом разница между $\large L_2$, $\large L_1$ и Huber loss достаточно заметна. При грамотном подборе параметра Huber loss мы даже получим наилучшую аппроксимацию функции среди наших вариантов. А еще на этом примере хорошо видна разница в условных квантилях (10%, 50% и 90% в нашем случае).


    К сожалению, функция Huber loss реализована не во всех современных библиотеках (в h2o реализована, в xgboost еще нет). Это же касается и других интересных функций потерь, включая и условные квантили, и такие экзотические вещи как условные экспектили. Но в целом, достаточно полезно знать о том, что такие варианты существуют и ими можно пользоваться.


    Функции потерь классификации


    Теперь разберем бинарную классификацию, когда $\large y \in \left\{-1, 1\right\}$. Мы уже видели, что с помощью GBM можно оптимизировать даже не очень дифференцируемые функции потерь. И вообще, можно было бы, не задумываясь, попытаться решить этот случай как еще одну задачу регрессии с каким-нибудь $\large L_2$ loss, но это будет не очень правильно (хотя и возможно).


    Из-за принципиально другой природы распределения целевой переменной, будем предсказывать и оптимизировать не сами метки классов, а их log-правдоподобие. Для этого переформулируем функции потерь над перемноженными предсказаниями и истинными метками $\large y \cdot f$ (не спроста же мы выбрали метки разных знаков). Наиболее известные варианты таких классификационных функций потерь:


    • $\large L(y, f) = log(1 + exp(-2yf))$, она же Logistic loss, она же Bernoulli loss. Интересное свойство заключается в том, что мы штрафуем даже корректно предсказанные метки классов. Нет, это не баг. Наоборот, оптимизируя эту функцию потерь, мы можем продолжать "раздвигать" классы и улучшать классификатор даже если все наблюдения предсказаны верно. Это самая стандартная и часто используемая функция потерь в бинарной классификации.
    • $\large L(y, f) = exp(-yf)$, оно же Adaboost loss. Так получилось, что классический Adaboost алгоритм эквивалентен GBM с этой функцией потерь. Концептуально эта функция потерь очень похожа на Logistic loss, но имеет более жесткий экспоненциальный штраф на ошибки классификации и используется реже.


    Сгенерируем новые игрушечные данные для задачи классификации. За основу возьмем наш зашумленный косинус, а в качестве классов целевой переменной будем использовать функцию sign. Новые данные выглядят следующим образом (jitter-шум добавлен для наглядности):



    Воспользуемся Logistic loss, чтобы посмотреть, что же мы на самом деле бустим. Как и прежде, соберем воедино то, что будем решать:


    • Игрушечные данные $\large \left\{ (x_i, y_i) \right\}_{i=1, \ldots,300}, y_i \in \left\{-1, 1\right\}$
    • Число итераций $\large M = 3$ ✓;
    • В качестве функции потерь — Logistic loss, ее градиент считается так:
      $\large r_{i} = \frac{2 \cdot y_i}{1 + exp(2 \cdot y_i \cdot \hat{f}(x_i)) }, \quad \mbox{for } i=1,\ldots,300$ ✓;
    • Деревья решений в качестве базовых алгоритмов $\large h(x)$ ✓;
    • Гиперпараметры деревьев решений: глубина деревьев равна 2 ✓;

    В этот раз с инициализацией алгоритма все немного сложнее. Во-первых, наши классы несбалансированы и разделены в пропорции примерно 63% на 37%. Во-вторых, аналитической формулы для инициализации для нашей функции потерь неизвестно. Так что будем искать $\large \hat{f_0} = \gamma$ поиском:



    Оптимальное начальное приближение нашлось в районе -0.273. Можно было догадаться, что оно будет отрицательным (нам выгоднее предсказывать всех наиболее популярным классом), но формулы точного значения, как мы уже сказали, нет. А теперь давайте наконец уже запустим GBM и посмотрим, что же на самом деле происходит под его капотом:



    Алгоритм отработал успешно, восстановив разделение наших классов. Можно видеть, как отделяются "нижние" области, в которых деревья больше уверены в корректном предсказании отрицательного класса, и как формируются две ступеньки, где классы были перемешаны. На псевдо-остатках видно, что у нас есть достаточно много корректно классифицированных наблюдений, и какое-то количество наблюдений с большими ошибками, которые появились из-за шума в данных. Как-то выглядит то, что на самом деле предсказывает GBM в задаче классификации (регрессия на псевдо-остатках логистической функции потерь).


    Веса


    Иногда возникает ситуация, когда для задачи хочется придумать более специфичную функцию потерь. Например, в предсказании финансовых рядов мы можем захотеть придавать больший вес крупным движениям временного ряда, а в задаче предсказания клиентского оттока — лучше предсказывать отток у клиентов с высоким LTV (lifetime value, сколько денег клиент нам принесет в будущем).



    Истинный путь статистического воина — придумать свою функцию потерь, выписать для нее производную (а для более эффективного обучения, еще и Гессиан), и тщательно проверить, удовлетворяет ли эта функция требуемым свойствам. Однако, высока вероятность где-то ошибиться, столкнуться с вычислительными трудностями, да и в целом потратить непозволительно много времени на исследования.


    Вместо этого был придуман очень простой инструмент, о котором редко вспоминают на практике — взвешивание наблюдений и задание весовых функций. Простейший пример такого взвешивания — задание весов для балансировки классов. В общем случае, если мы знаем, что какое-то подмножество данных, как во входных переменных $\large x$, так и в целевой переменной $\large y$ имеет большую значимость для нашей модели, мы просто задаем им больший вес $\large w(x,y)$. Главное — выполнить общие требования разумности весов:


    $ \large w_i \in \mathbb{R}, \\ \large w_i \geq 0 \quad \mbox{for } i=1,\ldots,n, \\ \large \sum_{i = 1}^n w_i > 0 $


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


    $ \large L_{w}(y,f) = w \cdot L(y,f), \\ \large r_{it} = - w_i \cdot \left[\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}\right]_{f(x)=\hat{f}(x)}, \quad \mbox{for } i=1,\ldots,n$


    Понятно, что для произвольных весов мы не знаем никаких красивых статистических свойств нашей модели. В общем случае, привязывая веса к значениям $\large y$, мы можем прострелить себе колено. Например, использование весов, пропорциональных $\large |y|$ в $\large L_1$ функции потерь — не эквивалентно $\large L_2$ loss, так как градиент не будет учитывать значения самих предсказаний $\large \hat{f}(x)$.


    Мы обсуждаем все это, чтобы лучше понимать наши возможности. Давайте придумаем какой-нибудь очень экзотический пример весов на наших игрушечных данных. Зададим сильно асимметричную весовую функцию следующим образом:


    $$display$$ \large \begin{equation} w(x) =\left\{ \begin{array}{@{}ll@{}} 0.1, & \text{if}\ x \leq 0 \\ 0.1 + |cos(x)|, & \text{if}\ x >0 \end{array}\right. \end{equation} $$display$$



    С помощью таких весов мы ожидаем увидеть два свойства: меньшую детализацию на отрицательных значениях $\large x$, а также форму функции, в большей степени похожую на исходный косинус. Все остальные настройки GBM мы берем из нашего предыдущего примера с классификацией, включая line search для оптимальных коэффициентов. Посмотрим, что у нас получилось:



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


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


    4. Итог про теорию GBM


    Сегодня мы разобрали основную теорию про градиентный бустинг. GBM — это не просто какой-то конкретный алгоритм, а общая методология, как строить ансамбли моделей. Причем методология достаточно гибкая и расширяемая: можно обучать большое количество моделей с учетом различных функций потерь и при этом еще и навешивать на них разного рода весовые функции.


    Как показывает и практика, и опыт соревнований по машинному обучению, в стандартных задачах (все кроме картинок и аудио, а также сильно разреженных данных), GBM очень часто является самым эффективным алгоритмом (не считая stacking-а и верхнеуровневых ансамблей, где GBM почти всегда является неотъемлемой их частью). При этом существуют адаптации GBM под Reinforcement Learning (Minecraft, ICML 2016), а алгоритм Виола-Джонса, до сих пор используемый в компьютерном зрении, основан на Adaboost.


    В этой статье мы специально опустили вопросы, связанные с регуляризацией GBM, стохастичностью, а также связанными с ними гиперпараметрами алгоритма. Неспроста мы везде выбирали маленькое число итераций алгоритма, $\large M = 3$. Если бы мы выбрали не 3, а 30 деревьев, и запустили GBM как описали выше, результат вышел бы не таким предсказуемым:



    О том, что делать в подобных ситуациях, и как взаимосвязаны друг с другом параметры регуляризации GBM, а также гиперпараметры базовых алгоритмов, мы расскажем в следующей статье. В ней же мы разберем современные пакеты — xgboost, lightgbm и h2o, а также попрактикуемся в их правильной настройке. А пока мы предлагаем вам поиграться с настройками GBM в еще одной очень крутой интерактивной демке Brliiantly wrong:


    http://arogozhnikov.github.io/2016/07/05/gradient_boosting_playground.html

    5. Домашнее задание №10


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


    Максимум за задание – 8 баллов. Дедлайн – 14 ноября 23:59 UTC +3.


    6. Полезные ссылки



    Особая благодарность yorko (Юрию Кашницкому) за ценные комментарии, а также bauchgefuehl (Анастасии Манохиной) за помощь с редактированием.

    Open Data Science 224,89
    Крупнейшее русскоязычное Data Science сообщество
    Поделиться публикацией
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама
    Комментарии 12
    • +2
      По-моему, замечательно. Просьба, Алексей, когда будете писать вторую часть, затронуть такую интересную штуку, как partial dependency plots, особенности интерпретации таких графиков. И сразу вопрос, планируете затронуть настройки класса H2OGradientBoostingEstimator для scikit-learn?
      • +3
        Спасибо :)

        Да, PDP планируются. И да, настройки H2O GBM тоже разберем. По-крайней мере, основные. H2O просто сильно угарели, и запилили штук 30 дополнительных настроек и небольших твиков. Все настройки перебирать смысла нет, обычно достаточно дергать 3-4 ручки (как именно как раз и будем разбираться).
      • +1

        Алексей, понимаю что это мелочи, но все же, у Вас в тексте статьи не отобразились пару формул. Статья хорошая, не все ещё понял, буду разбираться. Спасибо.

        • +1
          Спасибо :)

          Мы попробовали пофиксить этот косяк, но к сожалению, на мобильной версии Хабра ломаются \begin{equation} формулы. Сперва думали, что дело в переносах строк, но это не помогло.
        • +1

          Участникам курса


          Продублирую здесь информацию по окончанию курса. Что еще нас ожидает:


          • Домашняя работа №10 (~ 22 мая появится, времени на нее – неделя)
          • Тьюториалы, которые вы сами можете написать. Подробности тут (до 31 мая)
          • Соревнования Kaggle Inclass (до 29 мая)

          Встреча в московском офисе Mail.ru Group, посвященная финалу курса, пройдет 10 июня (а не в мае, как раньше предполагалось)

          • +1

            И еще ~ 1 июня выйдет 2 часть по бустингу, уже без домашнего задания либо с заданием по желанию.

            • +2

              Финальный рейтинг складывается из:


              • текущего рейтинга
              • 10 домашки (макс. 10 баллов)
              • тьюториалов (макс. 40 баллов)
              • соревнований (макс. по 40 баллов за каждое из двух)

              В веб-форме 10 домашки спросим, кто готов посетить митап 10 июня в московском офисе Mail.ru. Отвечать "да" стоит, если Вам это интересно, есть возможность приехать в Москву, а также если Вы попадаете в топ 200 рейтинга. Приглашены будут топ 100 по рейтингу, но очевидно, кто-то откажется, поэтому ставьте "да", если попадаете в топ 200.
              Программа митапа 10 июня будет оглашена позже, но там точно будет лекция и практика по нейронным сетям (CNN и RNN), будут рассмотрены приложения к задачам текстовой аналитики.

            • +1

              А зачем в logistic loss двойка?


              log(1 + exp(-2yf))

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

              • 0

                Двойка действительно мало на что влияет, кроме исторически сложившегося описания через odds ratio. Даже, вообще говоря, двух обоснований: и через биномиальное лог правдоподобие, и через кросс-энтропию (да, оптимизируем то же самое).


                Я хотел подчистую скопипастить оригинальный вывод из Elements of Statistical Learning, там ровно одна страница 346 (365 в pdf), "10.5 Why Exponential Loss?". Но чтото у меня ТеХ в комментах не работает :(


                А с инженерной точки зрения, эта двойка ни на что не влияет, кроме сравнения разных реализаций (чтобы числа одинаковые в обучении получались).

              • 0

                Участникам курса
                10 домашка по бустингу тут. Дедлайн — 30 мая 23.59 UTC +3.
                В конце веб-формы к этой домашке вопрос про участие в митапе 10 июня. Желательно, чтоб как можно больше участников курса ответили.

                • 0

                  Финишная прямая курса:


                  • Дедлайны по соревнованиям Kaggle Inclass – 29 мая. Надо переименовать свою команду из одного человека в точном соответствии с тем, как оно указано в рейтинге. В течение 2 дней после этого (до 23.59 среды 31 мая) надо загрузить свои решения – по Хабру сюда и по веб-сессиям Элис сюда. Решение должно быть воспроизводимым Python-скриптом (можно, например, тетрадку экспортировать в скрипт), по которому получается нужный файл посылки – обязательно сами прогоните перед тем как отправлять. Чтобы решение воспроизводилось быстро, закомментруйте подбор параметров. Например, если у вас есть блок с GridSearchCV, подбирающим глубину дерева (скажем, лучшее значение – 5), закомментируйте строку, где этот объект обучается (fit) и далее напишете, например, best_depth=5. Напомню, баллы за соревнования начисляются только тем, кто побил бенчмарки на приватном рейтинге (который как раз 29-го откроется), переименовал себя как надо и прислал воспроизводимое решение.
                  • Во вторник дедлайн по 10 домашке. Обратите внимание, что там справшивается про посещение митапа в Мэйле 10 июня. Желательно, чтоб как можно больше людей ответили.
                  • В среду последний день публикации тьюториалов и дедлайн по загрузке скриптов своих kaggle-решений
                  • До пятницы 2 июня учитываются голоса за тьюториалы (отдельно про это напишу, соберем все тьюториалы вместе). Также подводится финальный рейтинг.
                  • примерно в начале июня выйдет еще 2 часть про бустинг, уже без задания либо с опциональным заданием
                  • к 5 июня определяемся с итоговым списком посещения митапа в Мэйле 10 июня (про него еще отдельно сообщение будет)
                  • 0

                    Комментарии к домашнему заданию №10 (были неточности формулировки):


                    • используйте везде random_state=17
                    • log_loss в классификации считайте ненормированным (то есть не надо делить на число объектов) – как в формулах, предложенных во 2 вопросе задания. Лучше тетрадку с заданием смотреть в nbviewer-версии

                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                    Самое читаемое