Решающие деревья · Учебник по ML от ШАД

Учебник по ML от ШАД решающие деревья

Решающие деревья

Автор:  Филипп Синицин

Введение

В этой главе мы рассмотрим ещё одно семейство моделей машинного обучения — решающие деревья (decision trees).

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

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

Пример решающего дерева

Начнём с небольшого примера. На картинке ниже изображено дерево, построенное для задачи классификации на пять классов:

< {: style=”width:35vw” }

Объекты в этом примере имеют два признака с вещественными значениями: $X$ и $Y$. Решение о том, к какому классу будет отнесён текущий объект выборки, будет приниматься с помощью прохода от корня дерева к некоторому листу.

В каждом узле этого дерева находится предикат. Если предикат верен для текущего примера из выборки, мы переходим в правого потомка, если нет — в левого. В данном примере все предикаты — это просто взятие порога по значению какого-то признака:

\[B(x, j, t) = [ x_j \le t ]\]

В листьях записаны предсказания (например, метки классов). Как только мы дошли до листа, мы присваиваем объекту ответ, записанный в вершине.

На картинке ниже визуализирован процесс построения решающих поверхностей, порождаемых деревом (правая часть картинки):

<

Каждый предикат порождает разделение текущего подмножества пространства признаков на две части. На первом этапе, когда происходило деление по $[ X \le X_1 ]$, вся плоскость была поделена на две соответствующие части. На следующем уровне часть плоскости, для которой выполняется $X \le X_1$, была поделена на две части по значению второго признака $Y \le Y_1$ — так образовались области 1 и 2. То же самое повторяется для правой части дерева — и так далее до листовых вершин: получится пять областей на плоскости. Теперь любому объекту выборки будет присваиваться один из пяти классов в зависимости от того, в какую из образовавшихся областей он попадает.

Этот пример хорошо демонстрирует, в частности, то, что дерево осуществляет кусочно-постоянную аппроксимацию целевой зависимости. Ниже приведён пример визуализации решающей поверхности, которая соответствует дереву глубины 4, построенному для объектов данных из Ames Hounsing Dataset, где из всех признаков, описывающих объекты недвижимости, были выбраны ширина фасада (Lot_Frontage) и площадь (Lot_Area), а предсказать нужно стоимость. Для более понятной визуализации перед построением дерева из датасета были выкинуты объекты с Lot_Frontage > 150 и с Lot_Area > 20000.

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

Определение решающего дерева

Разобравшись с приведёнными выше примерами, мы можем дать определение решающего дерева. Пусть задано бинарное дерево, в котором:

  • каждой внутренней вершине $v$ приписан предикат \(B_v: \mathbb{X} \to \{ 0, 1 \}\);
  • каждой листовой вершине $v$ приписан прогноз $c_v \in \mathbb{Y}$, где $\mathbb{Y}$ — область значений целевой переменной (в случае классификации листу может быть также приписан вектор вероятностей классов).

В ходе предсказания осуществляется проход по этому дереву к некоторому листу. Для каждого объекта выборки $x$ движение начинается из корня. В очередной внутренней вершине $v$ проход продолжится вправо, если $B_v(x) = 1$, и влево, если $B_v(x) = 0$. Проход продолжается до момента, пока не будет достигнут некоторый лист, и ответом алгоритма на объекте $x$ считается прогноз $c_v$, приписанный этому листу.

Предикат $B_v$ может иметь, вообще говоря, произвольную структуру, но, как правило, на практике используют просто сравнение с порогом $t \in \mathbb{R}$ по какому-то $j$-му признаку:

\[B_v(x, j, t) = [ x_j \le t ]\]

При проходе через узел дерева с данным предикатом объекты будут отправлены в правое поддерево, если значение $j$-го признака у них меньше либо равно $t$, и в левое — если больше. В дальнейшем рассказе мы будем по умолчанию использовать именно такие предикаты.

Из структуры дерева решений следует несколько интересных свойств:

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

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

Выученная закономерность очень сложная; иногда она выделяет прямоугольником одну точку. И действительно, как мы увидим дальше, деревья умеют идеально подстраиваться под обучающую выборку. Это чемпионы переобучения. Помешать этому может лишь ограничение на высоту дерева. Вот как будет выглядеть дерево высоты 3, построенное на том же датасете:

Не всё идеально, но уже гораздо менее безумный классификатор.

Теперь обратимся к задаче регрессии. Обучим решающее дерево, не ограничивая его высоту.

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

Почему построение оптимального решающего дерева — сложная задача?

Пусть, как обычно, у нас задан датасет $(X, y)$, где \(y=\{y_i\}_{i=1}^N \subset \mathbb{R}^N\) — вектор таргетов, а \(X=\{x_i\}_{i = 1}^N \in \mathbb{R}^{N \times D}, x_i \in \mathbb{R}^D\) — матрица признаков, в которой $i$-я строка — это вектор признаков $i$-го объекта выборки. Пусть у нас также задана функция потерь $L(f, X, y)$, которую мы бы хотели минимизировать.

Наша задача — построить решающее дерево, наилучшим образом предсказывающее целевую зависимость. Однако, как уже было замечено выше, оптимизировать структуру дерева с помощью градиентного спуска не представляется возможным. Как ещё можно было бы решить эту задачу? Давайте начнём с простого — научимся строить решающие пни, то есть решающие деревья глубины 1.

Как и раньше, мы будем рассматривать только самые простые предикаты:

\[B_{j, t}(x_i) = [ x_{ij} \le t ]\]

Ясно, что задачу можно решить полным перебором: существует не более $(N - 1) D$ предикатов такого вида. Действительно, индекс $j$ (номер признака) пробегает значения от $1$ до $D$, а всего значений порога $t$, при которых меняется значение предиката, может быть не более $N - 1$:

Решение, которое мы ищем, будет иметь вид:

\[(j_{opt}, t_{opt}) = \arg \min_{j, t} L(B_{j, t}, X, y)\]

Для каждого из предикатов $B_{j, t}$ нам нужно посчитать значение функции потерь на всей выборке, что, в свою очередь, тоже занимает $O(N)$. Следовательно, полный алгоритм выглядит так:

min_loss = inf
optimal_border = None

for j in range(D):
    for t in X[:, j]:     # Можно брать сами значения признаков в качестве порогов
        loss = calculate_loss(t, j, X, y)
        if loss < min_loss:
           min_loss, optimal_border = loss, (j, t)

Сложность алгоритма — $O(N^2 D)$. Это не заоблачная сложность, хотя, конечно, не идеальная. Но это была схема возможного алгоритма поиска оптимального дерева высоты 1.

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

Как быть? Идеального ответа на этот вопрос нет, но до некоторой степени ситуацию можно улучшить двумя не исключающими друг друга способами:

  1. Разрешить себе искать не оптимальное решение, а просто достаточно хорошее. Начать можно с того, чтобы строить дерево с помощью жадного алгоритма, то есть не искать всю структуру сразу, а строить дерево этаж за этажом. Тогда в каждой внутренней вершине дерева будет решаться задача, схожая с задачей построения решающего пня. Для того чтобы этот подход хоть как-то работал, его придётся прокачать внушительным набором эвристик.
  2. Заняться оптимизацией с точки зрения computer science — наивную версию алгоритма (перебор наборов возможных предикатов и порогов) можно ускорить и асимптотически, и в константу раз.

Эти две идеи мы и будем обсуждать в дальнейшем. Сначала попытаемся подробно разобраться с первой — использованием жадного алгоритма.

Жадный алгоритм построения решающего дерева

Пусть $X$ — исходное множество объектов обучающей выборки, а $X_m$ — множество объектов, попавших в текущий лист (в самом начале $X_m = X$). Тогда жадный алгоритм можно верхнеуровнево описать следующим образом:

  1. Создаём вершину $v$.
  2. Если выполнен критерий остановки $Stop(X_m)$, то останавливаемся, объявляем эту вершину листом и ставим ей в соответствие ответ $Ans(X_m)$, после чего возвращаем её.
  3. Иначе: находим предикат (иногда ещё говорят сплит) $B_{j, t}$, который определит наилучшее разбиение текущего множества объектов $X_m$ на две подвыборки $X_{\ell}$ и $X_r$, максимизируя критерий ветвления $Branch(X_m, j, t)$.
  4. Для $X_\ell$ и $X_r$ рекурсивно повторим процедуру.

Данный алгоритм содержит в себе несколько вспомогательных функций, которые надо выбрать так, чтобы итоговое дерево было способно минимизировать $L$:

  • $Ans(X_m)$, вычисляющая ответ для листа по попавшим в него объектам из обучающей выборки, может быть, например:
    • в случае задачи классификации — меткой самого частого класса или оценкой дискретного распределения вероятностей классов для объектов, попавших в этот лист;
    • в случае задачи регрессии — средним, медианой или другой статистикой;
    • простой моделью. К примеру, листы в дереве, задающем регрессию, могут быть линейными функциями или синусоидами, обученными на данных, попавших в лист. Впрочем, везде ниже мы будем предполагать, что в каждом листе просто предсказывается константа.
  • Критерий остановки $Stop(X_m)$ — функция, которая решает, нужно ли продолжать ветвление или пора остановиться.

Это может быть какое-то тривиальное правило: например, остановиться только в тот момент, когда объекты в листе получились достаточно однородными и/или их не слишком много. Более детально мы поговорим о критериях остановки в параграфе про регуляризацию деревьев.

  • Критерий ветвления $Branch(X_m, feature, value)$ — пожалуй, самая интересная компонента алгоритма. Это функция, измеряющая, насколько хорош предлагаемый сплит. Чаще всего эта функция оценивает, насколько улучшится некоторая финальная метрика качества дерева в случае, если получившиеся два листа будут терминальными, по сравнению с ситуацией, когда сама исходная вершина является листом. Выбирается такой сплит, который даёт наиболее существенное улучшение. Впрочем, есть и другие подходы.

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

Критерии ветвления: общая идея

Давайте теперь по очереди посмотрим на популярные постановки задач ML и под каждую подберём свой критерий.

Ответы дерева будем кодировать так: \(c \in \mathbb{R}\) — для ответов регрессии и меток класса; для случаев, когда надо ответить дискретным распределением на классах, $c \in \mathbb{R}^K$ будет вектором вероятностей:

\[c = (c_1, \ldots, c_K), \quad \sum_{i = 1}^K c_i = 1\]

Предположим также, что задана некоторая функция потерь $L(y_i, c)$. О том, что это может быть за функция, мы поговорим ниже.

В момент, когда мы ищем оптимальный сплит $X_m = X_l\sqcup X_r$, мы можем вычислить для объектов из $X_m$ тот константный таргет $c$, которые предсказало бы дерево, будь текущая вершина терминальной, и связанное с ними значение исходного функционала качества $L$. А именно — константа $c$ должна минимизировать среднее значение функции потерь:

\[\frac{1}{\vert X_m\vert}\sum\limits_{(x_i, y_i) \in X_m} L(y_i, c)\]

Оптимальное значение этой величины

\[H(X_m) = \min\limits_{c \in Y} \frac{1}{\vert X_m\vert}\sum\limits_{(x_i, y_i) \in X_m} L(y_i, c)\]

обычно называют информативностью, или impurity. Чем она ниже, тем лучше объекты в листе можно приблизить константным значением.

Похожим образом можно определить информативность решающего пня. Пусть, как и выше, $X_l$ — множество объектов, попавших в левую вершину, а $X_r$ — в правую; пусть также $c_l$ и $c_r$ — константы, которые предсказываются в этих вершинах. Тогда функция потерь для всего пня в целом будет равна

\[\frac1{|X_m|}\left(\sum_{x_i\in X_l}L(y_i, c_l) + \sum_{x_i\in X_r}L(y_i, c_r)\right) =\]

Вопрос на подумать. Как информативность решающего пня связана с информативностью его двух листьев?

Ответ (не открывайте сразу; сначала подумайте сами!)

Преобразуем выражение:

\[\frac1{|X_m|}\left(\sum_{x_i\in X_l}L(y_i, c_l) + \sum_{x_i\in X_r}L(y_i, c_r)\right) =\] \[= \frac1{|X_m|}\left(|X_l|\cdot\frac{1}{|X_l|}\sum_{x_i\in X_l}L(y_i, c_l) + |X_r|\cdot\frac1{|X_r|}\sum_{x_i\in X_r}L(y_i, c_r)\right)\]

Эта сумма будет минимальна при оптимальном выборе констант $c_l$ и $c_r$, и информативность пня будет равна

\[\frac{|X_l|}{|X_m|}H(X_l) + \frac{|X_r|}{|X_m|}H(X_r)\]

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

Разность информативности исходной вершины и решающего пня равна

\[H(X_m) - \frac{|X_l|}{|X_m|}H(X_l) - \frac{|X_r|}{|X_m|}H(X_r)\]

Для симметрии её приятно умножить на $\vert X_m\vert$; тогда получится следующий критерий ветвления:

\[\color{#348FEA}{Branch (X_m, j, t) = |X_m| \cdot H(X_m) - |X_l| \cdot H(X_l) - |X_r| \cdot H(X_r)}\]

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

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

Информативность в задаче регрессии: MSE

Посмотрим на простой пример — регрессию с минимизацией среднеквадратичной ошибки:

\[L(y_i, c) = (y_i -c)^2\]

Информативность листа будет выглядеть следующим образом:

\[H(X_m) = \frac{1}{\vert X_m\vert} \min\limits_{c \in Y} \sum\limits_{(x_i, y_i) \in X_m} (y_i - c)^2\]

Как мы уже знаем, оптимальным предсказанием константного классификатора для задачи минимизации MSE является среднее значение, то есть

\[c = \frac{\sum y_i}{|X_m|}\]

Подставив в формулу информативности сплита, получаем:

\[\color{#348FEA}{H(X_m) = \sum\limits_{(x_i, y_i) \in X_m}\frac{\left(y_i - \overline{y} \right)^2}{|X_m|}, ~ \text{где} ~ \overline{y} = \frac{1}{\vert X_m \vert} \sum_i y_i}\]

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

Информативность в задаче регрессии: MAE

\[L(y_i, c) = |y_i -c|\]

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

В качестве информативности выступает абсолютное отклонение от медианы:

\[H(X_m) = \sum\limits_{(x_i, y_i) \in X_m}\frac{|y_i - MEDIAN(Y)|}{|X_m|}\]

Критерий информативности в задаче классификации: misclassification error

Пусть в нашей задаче $K$ классов, а $p_k$ — доля объектов класса $k$ в текущей вершине $X_m$:

\[p_k = \frac{1}{|X_m|} \sum_{(x_i, y_i) \in X_m} \mathbb{I}[ y_i = k ]\]

Допустим, мы заботимся о доле верно угаданных классов, то есть функция потерь — это индикатор ошибки:

\[L(y_i, c) = \mathbb{I}[y_i \ne c]\]

Пусть также предсказание модели в листе — один какой-то класс. Информативность для такой функции потерь выглядит так:

\[H(X_m) = \min_{c \in Y} \frac{1}{|X_m|} \sum_{(x_i, y_i) \in X_m} \mathbb{I}[y_i \ne c]\]

Ясно, что оптимальным предсказанием в листе будет наиболее частотный класс $k_{\ast}$, а выражение для информативности упростится следующим образом:

\[H(X_m) = \frac{1}{|X_m|} \sum_{(x_i, y_i) \in X_m} \mathbb{I}[y_i \ne k_{\ast}] = 1 - p_{k_{\ast}}\]

Информативность в задаче классификации: энтропия

Если же мы собрались предсказывать вероятностное распределение классов $(c_1, \ldots, c_K)$, то к этому вопросу можно подойти так же, как мы поступали при выводе логистической регрессии: через максимизацию логарифма правдоподобия (= минимизацию минус логарифма) распределения Бернулли. А именно, пусть в вершине дерева предсказывается фиксированное распределение $c$ (не зависящее от $x_i$), тогда правдоподобие имеет вид

\[P(y\vert x, c) = P(y\vert c) = \prod\limits_{(x_i, y_i) \in X_m}P(y_i\vert c) = \prod\limits_{(x_i, y_i) \in X_m} \prod\limits_{k = 1}^K c_k^{\mathbb{I}[ y_i = k ]},\]

откуда

\[H(X_m) = \min\limits_{\sum\limits_{k} c_k = 1} \left( -\frac{1}{|X_m|} \sum\limits_{(x_i, y_i) \in X_m} \sum\limits_{k = 1}^K \mathbb{I}[ y_i = k ] \log c_k \right)\]

То, что оценка вероятностей в листе $c_k$, минимизирующая $H(X_m)$, должна быть равна $p_k$, то есть доле попавших в лист объектов этого класса, до некоторой степени очевидно, но это можно вывести и строго.

Доказательство для любопытных

Из-за наличия условия на $\sum_k c_k = 1$ нам придётся минимизировать лагранжиан

\[L(c, \lambda) = \min_{c, \lambda} \left( \left( -\frac{1}{|X_m|} \sum\limits_{(x_i, y_i) \in X_m} \sum\limits_{k = 1}^K \mathbb{I}[ y_i = k ] \log c_k \right) + \lambda \sum\limits_{k = 1}^K c_k \right)\]

Как обычно, возьмём частную производную и решим уравнение:

\[0 = \frac{\partial }{\partial c_j}L(c, \lambda) =\] \[= \frac{\partial }{\partial c_j} \left( \left( -\frac{1}{|X_m|} \sum\limits_{(x_i, y_i) \in X_m} \sum\limits_{k = 1}^K \mathbb{I}[ y_i = k ] \log c_k \right) + \lambda \sum\limits_{k = 1}^K c_k \right) =\] \[= \left( \left( -\frac{1}{|X_m|} \sum\limits_{(x_i, y_i) \in X_m} \mathbb{I}[ y_i = j ] \frac{1}{c_j} \right) + \lambda \right) = - \frac{p_j}{c_j} + \lambda,\]

откуда выражаем $c_j = \frac{p_j}{\lambda}$. Суммируя эти равенства, получим:

\[1 = \sum_{k = 1}^K c_k = \frac{1}{\lambda} \sum_{k = 1}^K p_k = \frac{1}{\lambda},\]

откуда $\lambda = 1$, а значит, $c_k = p_k$.

Подставляя вектор \(c = (p_1,\ldots,p_K)\) в выражение выше, мы в качестве информативности получим энтропию распределения классов:

\[\color{#348FEA}{H(X_m) = -\sum_{k = 1}^K p_k \log p_k}\]
Немного подробнее об энтропии

Величина

\[H(X_m) = -\sum_k p_k \log p_k\]

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

Так как $p_k\in[0,1]$, энтропия неотрицательна. Если случайная величина принимает только одно значение, то она абсолютно предсказуема и её энтропия равна $ - 1 \log(1) = 0$.

Наибольшего значения энтропия достигает для равномерно распределённой случайной величины — и это отражает тот факт, что среди всех величин с данной областью значений она наиболее «непредсказуема». Для равномерно распределённой на множестве $\{1,\ldots,K\}$ случайной величины значение энтропии будет равно:

\[-\sum_{k = 1}^K \frac{1}{K}\log\frac{1}{K} = \log K\]

На следующем графике приведены три дискретных распределения на множестве \(\{0, 1, \ldots, 20\}\) с их энтропиями. Как и указано выше, максимальную энтропию будет иметь равномерное распределение; у двух других проявляются пики разной степени остроты — и тем самым реализации этих величин обладают меньшей неопределённостью: мы можем с большей уверенностью говорить, что будет сгенерировано.

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

В изначальном датасете по 25 точек каждого класса; энтропия состояния равна

\[S_0 = -\frac{25}{50}\log_2{\frac{25}{50}}-\frac{25}{50}\log_2{\frac{25}{50}} = 1\]

Для первого разбиения, по $[X_1 \leqslant 3]$ в левую часть попадают $25$ точек класса 0 и $12$ точек класса 1, а в правую — $0$ точек класса 0 и $13$ точек класса 1. Энтропия левой группы равна

\[S_l = -\frac{25}{37}\log_2{\frac{25}{37}}-\frac{12}{37}\log_2{\frac{12}{37}}\approx 0.9\]

Энтропия правой группы, строго говоря, не определена (под логарифмом ноль), но если заменить несуществующее значение на $\lim_{t\rightarrow 0+}t\log_2{t} = 0$, то получится

\[S_r = - 0 - 1\log_2{1}=0\]

Что, в принципе, логично: с вероятностью 1 выпадает единица, мы всегда уверены в результате, и энтропия у такого, вырожденного распределения тоже минимальная, равная нулю.

Как можно заметить, энтропия в обеих группах уменьшилась по сравнению с изначальной. Получается, что, разделив шарики по значению координаты, равному 12, мы уменьшили общую неопределённость системы. Но какое из разбиений даёт самое радикальное уменьшение? После несложных вычислений, мы получаем для нарисованных выше разбиений:

\[Branch^{(1)} \approx 16.4\] \[Branch^{(2)} \approx 30.5\] \[Branch^{(3)} \approx 7\]

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

Кстати, Шеннон изначально собирался назвать информационную энтропию или «информацией», или «неопределённостью», но в итоге выбрал название «энтропия», потому что концепция со схожим смыслом в статистической механике уже была названа энтропией. Употребление термина из другой научной области выглядело убедительным преимуществом при ведении научных споров.

Информативность в задаче классификации: критерий Джини

Пусть предсказание модели — это распределение вероятностей классов $(c_1, \ldots, c_k)$. Вместо логарифма правдоподобия в качестве критерия можно выбрать, например, метрику Бриера (за которой стоит всего лишь идея посчитать MSE от вероятностей). Тогда информативность получится равной

\[H(X_m) = \min_{\sum_k c_k = 1} \frac{1}{|X_m|} \sum_{(x_i, y_i) \in X_m} \sum_{k = 1}^K (c_k - \mathbb{I}[ y_i = k ] )^2\]

Можно показать, что оптимальное значение этой метрики, как и в случае энтропии, достигается на векторе $c$, состоящем из выборочных оценок частот классов $(p_1, \ldots, p_k)$, $p_i = \frac{1}{\vert X_m\vert}\sum_i \mathbb{I}[ y_i = k ]$. Если подставить $(p_1, \ldots, p_k)$ в выражение выше и упростить его, получится критерий Джини:

\[\color{#348FEA}{H(X_m) = \sum_{k = 1}^K p_k (1 - p_k)}\]

Критерий Джини допускает и следующую интерпретацию: $H(X_m)$ равно математическому ожиданию числа неправильно классифицированных объектов в случае, если мы будем приписывать им случайные метки из дискретного распределения, заданного вероятностями $(p_1, \ldots, p_k)$.

Неоптимальность полученных критериев

Казалось бы, мы вывели критерии информативности для всех популярных задач, и они довольно логично следуют из их постановок, но получилось ли у нас обмануть NP-полноту и научиться строить оптимальные деревья легко и быстро?

Конечно, нет. Простейший пример — решение задачи XOR с помощью жадного алгоритма и любого критерия, который мы построили выше:

XOR

(источник картинки)

Вне зависимости от того, что вы оптимизируете, жадный алгоритм не даст оптимального решения задачи XOR. Но этим примером проблемы не исчерпываются. Скажем, бывают ситуации, когда оптимальное с точки зрения выбранной метрики дерево вы получите с критерием ветвления, построенным по другой метрике (например, MSE-критерий для MAE-задачи или Джини для misclassification error).

Особенности данных

Категориальные признаки

На первый взгляд, деревья прекрасно могут работать с категориальными переменными. А именно, если признак $x^i$ принимает значения из множества $C = \{c_1,\ldots,c_M\}$, то при очередном разбиении мы можем рассматривать по этому признаку произвольные сплиты вида $C = C_l\sqcup C_r$ (предикат будет иметь вид $[x^i\in C_r]$). Это очень логично и естественно; проблема в том, что при больших $M$ у нас будет $2^{M-1}-1$ сплитов, и перебирать их будет слишком долго. Было бы здорово уметь каким-то образом упорядочивать значения $c_m$, чтобы работать с ними так же, как с обычными числами: разделяя на значения, «не превосходящие» и «большие» определённого порога.

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

Так, для задачи бинарной классификации значения $c_m$ можно упорядочить по неубыванию доли объектов класса 1 с $x^i = c_m$, после чего работать с ними, как со значениями вещественного признака. Показано, что в случае, если мы выбираем таким образом сплит, оптимальный с точки зрения энтропийного критерия или критерия Джини, то он будет оптимальным среди всех $2^{M-1}-1$ сплитов.

Для задачи регрессии с функцией потерь MSE значения $c_m$ можно упорядочивать по среднему значению таргета на подмножестве $\{X\mid x^i = c_m\}$. Полученный таким образом сплит тоже будет оптимальным.

Работа с пропусками

Одна из приятных особенностей деревьев — это способность обрабатывать пропуски в данных. Разберёмся, что при этом происходит на этапе обучения и на этапе применения дерева.

Пусть у нас есть некоторый признак $x^i$, значение которого пропущено у некоторых объектов. Как обычно, обозначим через $X_m$ множество объектов, пришедших в рассматриваемую вершину, а через $V_m$ — подмножество $X_m$, состоящее из объектов с пропущенным значением $x^i$. В момент выбора сплитов по этому признаку мы будет просто игнорировать объекты из $V_m$, а когда сплит выбран, мы отправим их в оба поддерева. При этом логично присвоить им веса: $\frac{\vert X_l\vert}{\vert X_m\vert}$ для левого поддерева и $\frac{\vert X_r\vert}{\vert X_m\vert}$ для правого. Веса будут учитываться как коэффициенты при $L(y_i, c)$ в формуле информативности.

Вопрос на подумать. Во всех критериях ветвления участвуют мощности множеств $X_m$, $X_l$ и $X_r$. Нужно ли уменьшение размера выборки учитывать в формулах для информативности? Если нужно, то как?

Ответ (не открывайте сразу; сначала подумайте сами!)

Чтобы ответить на этот вопрос, вспомним, как работают критерии ветвления. Их суть в сравнении информативности $H(X_m)$ для исходной вершины с информативностью для решающего пня, порождённого рассматриваемым сплитом. Величина $H(X_m)$ никак не меняется; надо понять, что будет с пнём. Для него можно посчитать информативность так, словно объектов с пропущенными значениями нет:

\[\frac1{|X_m\setminus V_m|}\left(\sum_{x_i\in X_l\setminus V_m}L(y_i, c_l) + \sum_{x_i\in X_r\setminus V_m}L(y_i, c_r)\right)\]

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

\[Branch = |X_m\setminus V_m|\cdot H(X_m) - |X_l\setminus V_m|\cdot H(X_l\setminus V_m) - |X_r\setminus V_m|\cdot H(X_r\setminus V_m)\]

Если же мы предполагаем, что объекты из $V_m$ отправляются в новые листья с весами, как описано выше, то формула оказывается другой, и коэффициенты можно не менять:

\[\frac1{|X_m|}\left(\underbrace{\sum_{x_i\in X_l\setminus V_m}L(y_i, c_l) + \frac{|X_l|}{|X_m|}\sum_{x_i\in V_m}L(y_i, c_l)}_{\text{Левый лист}} + \underbrace{\sum_{x_i\in X_r\setminus V_m}L(y_i, c_r) + \frac{|X_r|}{|X_m|}\sum_{x_i\in V_m}L(y_i, c_r)}_{\text{Правый лист}}\right)\]

 

Теперь рассмотрим этап применения дерева. Допустим, в вершину, где сплит идёт по $i$-му признаку, пришёл объект $x_0$ с пропущенным значением этого признака. Предлагается отправить его в каждую из дальнейших веток и получить по ним предсказания $\widehat{y}_l$ и $\widehat{y}_r$. Эти предсказания мы усредним с весами $\frac{\vert X_l\vert}{\vert X_m\vert}$ и $\frac{\vert X_r\vert}{\vert X_m\vert}$ (которые мы запомнили в ходе обучения):

\[\widehat{y} = \frac{\vert X_l\vert}{\vert X_m\vert}\widehat{y}_l + \frac{\vert X_r\vert}{\vert X_m\vert}\widehat{y}_r\]

Для задачи регрессии это сразу даст нам таргет, а в задаче бинарной классификации — оценку вероятности класса 1.

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

Методы регуляризации решающих деревьев

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

Для этого есть разные критерии, обычно используются все сразу:

  • ограничение по максимальной глубине дерева;
  • ограничение на минимальное количество объектов в листе;
  • ограничение на максимальное количество листьев в дереве;
  • требование, чтобы функционал качества $Branch$ при делении текущей подвыборки на две улучшался не менее чем на $s$ процентов.

Делать это можно на разных этапах работы алгоритма, что не меняет сути, но имеет разные устоявшиеся названия:

  • можно проверять критерии прямо во время построения дерева, такой способ называется pre-pruning или early stopping;
  • а можно построить дерево жадно без ограничений, а затем провести стрижку (pruning), то есть удалить некоторые вершины из дерева так, чтобы итоговое качество упало не сильно, но дерево начало подходить под условия регуляризации. При этом качество стоит измерять на отдельной, отложенной выборке.

Алгоритмические трюки

Теперь временно снимем шапочку ML-аналитика, наденем шапочку разработчика и специалиста по computer science и посмотрим, как можно сделать полученный алгоритм более вычислительно эффективным.

В базовом алгоритме мы в каждой вершине дерева для всех возможных значений сплитов вычисляем информативность. Если в вершину пришло $q$ объектов, то мы рассматриваем $qD$ сплитов и для каждого тратим $O(q)$ операций на подсчёт информативности. Отметим, что в разных вершинах, находящихся в нашем дереве на одном уровне, оказываются разные объекты, то есть сумма этих $q$ по всем вершинам заданного уровня не превосходит $N$, а значит, выбор сплитов во всех вершинах уровня потребует $O(N^2 D)$ операций. Таким образом, общая сложность построения дерева — $O(hN^2 D)$ (где $h$ — высота дерева), и доминирует в ней перебор всех возможных предикатов на каждом уровне построения дерева. Посмотрим, что с этим можно сделать.

Динамическое программирование

Постараемся оптимизировать процесс выбора сплита в одной конкретной вершине.

Вместо того чтобы рассматривать все $O(ND)$ возможных сплитов, для каждого тратя $O(N)$ на вычисление информативности, можно использовать одномерную динамику. Для этого заметим, что если отсортировать объекты по какому-то признаку, то, проходя по отсортированному массиву, можно одновременно и перебирать все значения предикатов, и поддерживать все необходимые статистики для пересчёта значений информативности за $O(1)$ для каждого следующего варианта сплита (против изначальных $O(N)$).

Давайте разберём, как это работает, на примере построения дерева для MSE. Чтобы оценить информативность для листа, нам нужно знать несколько вещей:

  • дисперсию и среднее значение таргета в текущем листе;
  • дисперсию и среднее значение таргета в обоих потомках для каждого потенциального значения сплита.

Дисперсию и среднее текущего листа легко посчитать за $O(n)$.

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

\[\overline{Y} = \frac1N\sum y_i,\] \[\sigma^2 (Y) = \overline{Y^2} - (\overline{Y})^2 = \frac1N\sum y_i^2 - \frac{1}{N^2}(\sum y_i)^2\]

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

Если в вершину дерева пришло $q$ объектов, сложность построения одного сплита складывается из $D$ сортировок каждая по $O(q\log q)$ и одного линейного прохода с динамикой, всего $O(qD\log q + qD) = O(qD\log q)$, что лучше исходного $O(q^2D)$. Итоговая сложность алгоритма построения дерева — $O(hND\log N)$ (где $h$ – высота дерева) против $hN^2D$ в наивной его версии.

Какие именно статистики накапливать (средние, медианы, частоты), зависит от критерия, который вы используете.

Гистограммный метод

Если бы мощность множества значений признаков была ограничена какой-то разумной константой $b \ll N$, то сортировку в предыдущем способе можно было бы заменить сортировкой подсчётом и за счёт этого существенно ускорить алгоритм: ведь сложность такой сортировки — $O(N)$.

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

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

discr

Аналогичная процедура проводится для всех признаков выборки. Полная сложность предобработки — $O(DN\log N)$ (сортировка за $O(N\log N)$ для каждого из $D$ признаков).

Теперь в процедуре динамического алгоритма поиска оптимального сплита нам надо перебирать не все $N$ объектов выборки, а всего лишь $b$ подготовленных заранее границ подотрезков. Частичные суммы статистик тоже придётся поддерживать не для исходного массива данных, а для списка из $b$ возможных сплитов. А для того чтобы делать это эффективно, необходим объект, называемый гистограммой: упорядоченный словарь, сопоставляющий каждому значению дискретизированного признака сумму необходимой статистики от таргета на отрезке [B[i-1], B[i]].

Финальный вид алгоритма таков:

  1. Дискретизируем каждый из признаков на $b$ значений. Сложность $O(DN\log N)$.
  2. Создаём корневую вершину root.
  3. Вызываем build_tree_recursive(root, data).

Функция build_tree_recursive выглядит следующим образом:

  1. Проверяем, не пора ли остановиться. Если пора — считаем значение в листе.
  2. Теперь мы снова используем динамический алгоритм, но объекты будем сортировать не по исходным значениям признаков, а по их дискретизированным версиям, упорядочивая их с помощью сортировки подсчётом (для вершины, в которую попало $q$ объектов, сложность будет равна $O(qD)$ против $O(q\log q\cdot D)$ в стандартной динамике).
  3. Находим оптимальный сплит за $O(qD)$.
  4. Делим данные, запускаем процедуру рекурсивно для обоих поддеревьев.

Общая сложность: $O(DN\log N + hND)$

Гистограммный метод на примере

Давайте посчитаем информативность сплита для MSE на данных с одним признаком, по которому мы уже отсортировали данные:

x -3 -2 -0.05 1 1 2 6 8
y 0 0.5 -1 0 1 2 1 4

Дискретизируем на три отрезка с границами $B = [-1, 1.5]$

discretized_x 0 0 1 1 1 2 2 2
y 0 0.5 -1 0 1 2 1 4

Строим гистограмму частичных сумм, сумм квадратов и количества объектов для двух отрезков (для третьего можно посчитать, используя общую сумму):

b 0 1
cnt 2 3
sum_y 0.5 0
sum_y_sq 0.25 2

Всего объектов восемь, сумма таргетов — $\text{sum_total} = 7.5$, сумма квадратов — $\text{sum_sq_total} = 23.25$.

Информативность текущего листа:

\[D[X_d] = \frac{\text{sum_sq_total}}{\text{cnt_total}} - \left(\frac{\text{sum_total}}{\text {cnt_total}}\right)^2 = \\ = \frac{23.25}{8} - \left( \frac{7.5}{8} \right)^{2} = 2.0273\]

Посчитаем значение критерия ветвления для $b=-1$. Информативность левого листа, то есть дисперсия в нём, равна:

\[D[X_l] = E[X_l^2] - E^2[X_l] = \\ = \frac{\text{sum_y_sq}[0]}{\text{cnt}[0]} - \left(\frac{\text{sum_y}}{\text{cnt}[0]}\right)^2 = \\ = 0.125 - 0.25 ^ 2 = 0.0625\]

Информативность правого листа, то есть дисперсия в нём, равна:

\[D[X_r] = E[X_r^2] - E^2[X_r] = \\ = \frac{\text{sum_sq_total} - \text{sum_y_sq}[0]}{\text{cnt_total} - \text{cnt}[0]} - \left(\frac{\text{sum_total} - \text{sum_y}}{\text{cnt_total} - \text{cnt}[0]}\right)^2 \\ \approx 2.47\]

Полное значение критерия для разбиения по $b=-1$:

\[Branch(b = -1) = 8 \cdot D[X_d] - 2 \cdot D[X_l] - 6 \cdot D[X_r] = \\ = 6.555\]

Mixed integer optimization

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

Историческая справка

Как вы, может быть, уже заметили, решающие деревья — это одна большая эвристика для решения NP-полной задачи, практически лишённая какой-либо стройной теоретической подоплёки. В 1970–1990-e годы интерес к ним был весьма велик как в индустрии, где был полезен хорошо интерпретируемый классификатор, так и в науке, где учёные интересовались способами приближённого решения NP-полных задач.

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

На сегодня это всё потеряло актуальность в связи с тем, что существуют хорошо написанные библиотеки (например, sklearn, в которой реализована оптимизированная версия CART).