В этой главе мы рассмотрим ещё одно семейство моделей машинного обучения — решающие деревья (decision trees).
Решающее дерево предсказывает значение целевой переменной с помощью применения последовательности простых решающих правил (которые называются предикатами). Этот процесс в некотором смысле согласуется с естественным для человека процессом принятия решений.
Хотя обобщающая способность решающих деревьев невысока, их предсказания вычисляются довольно просто, из-за чего решающие деревья часто используют как кирпичики для построения ансамблей — моделей, делающих предсказания на основе агрегации предсказаний других моделей. О них мы поговорим в следующей главе.
Начнём с небольшого примера. На картинке ниже изображено дерево, построенное для задачи классификации на пять классов:
Объекты в этом примере имеют два признака с вещественными значениями: $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
.
В каждой из прямоугольных областей предсказывается одна и та же стоимость.
Разобравшись с приведёнными выше примерами, мы можем дать определение решающего дерева. Пусть задано бинарное дерево, в котором:
В ходе предсказания осуществляется проход по этому дереву к некоторому листу. Для каждого объекта выборки $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-полной задачей, то есть человечеству пока неизвестны способы решить её за полиномиальное время.
Как быть? Идеального ответа на этот вопрос нет, но до некоторой степени ситуацию можно улучшить двумя не исключающими друг друга способами:
Эти две идеи мы и будем обсуждать в дальнейшем. Сначала попытаемся подробно разобраться с первой — использованием жадного алгоритма.
Пусть $X$ — исходное множество объектов обучающей выборки, а $X_m$ — множество объектов, попавших в текущий лист (в самом начале $X_m = X$). Тогда жадный алгоритм можно верхнеуровнево описать следующим образом:
Данный алгоритм содержит в себе несколько вспомогательных функций, которые надо выбрать так, чтобы итоговое дерево было способно минимизировать $L$:
Это может быть какое-то тривиальное правило: например, остановиться только в тот момент, когда объекты в листе получились достаточно однородными и/или их не слишком много. Более детально мы поговорим о критериях остановки в параграфе про регуляризацию деревьев.
При этом строгой теории, которая бы связывала оптимальность выбора разных вариантов этих функций и разных метрик классификации и регрессии, в общем случае не существует. Однако есть набор интуитивных и хорошо себя зарекомендовавших соображений, с которыми мы вас сейчас познакомим.
Давайте теперь по очереди посмотрим на популярные постановки задач 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)}\]Получившаяся величина неотрицательна: ведь, разделив объекты на две кучки и подобрав ответ для каждой, мы точно не сделаем хуже. Кроме того, она тем больше, чем лучше предлагаемый сплит.
Теперь посмотрим, какими будут критерии ветвления для типичных задач.
Посмотрим на простой пример — регрессию с минимизацией среднеквадратичной ошибки:
\[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 констатного предсказателя (мы это обсуждали в главе про линейные модели).
В качестве информативности выступает абсолютное отклонение от медианы:
\[H(X_m) = \sum\limits_{(x_i, y_i) \in X_m}\frac{|y_i - MEDIAN(Y)|}{|X_m|}\]Пусть в нашей задаче $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. Но этим примером проблемы не исчерпываются. Скажем, бывают ситуации, когда оптимальное с точки зрения выбранной метрики дерево вы получите с критерием ветвления, построенным по другой метрике (например, 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.
Замечание. Если речь идёт о категориальном признаке, может оказаться хорошей идеей ввести дополнительное значение «пропущено» для категориального признака и дальше работать с пропусками, как с обычным значением. Особенно это актуально в ситуациях, когда пропуски имеют системный характер и их наличие несёт в себе определённую информацию.
Мы уже упоминали выше, что деревья легко переобучаются и процесс ветвления надо в какой-то момент останавливать.
Для этого есть разные критерии, обычно используются все сразу:
Делать это можно на разных этапах работы алгоритма, что не меняет сути, но имеет разные устоявшиеся названия:
Теперь временно снимем шапочку 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$ подотрезков, длины которых выбираются так, чтобы в каждый попадало примерно равное число обучающих примеров. После чего значения признака заменяются на номера отрезков, на которые они попали.
Аналогичная процедура проводится для всех признаков выборки. Полная сложность предобработки — $O(DN\log N)$ (сортировка за $O(N\log N)$ для каждого из $D$ признаков).
Теперь в процедуре динамического алгоритма поиска оптимального сплита нам надо перебирать не все $N$ объектов выборки, а всего лишь $b$ подготовленных заранее границ подотрезков. Частичные суммы статистик тоже придётся поддерживать не для исходного массива данных, а для списка из $b$ возможных сплитов. А для того чтобы делать это эффективно, необходим объект, называемый гистограммой: упорядоченный словарь, сопоставляющий каждому значению дискретизированного признака сумму необходимой статистики от таргета на отрезке [B[i-1], B[i]].
Финальный вид алгоритма таков:
root
.build_tree_recursive(root, data)
.Функция build_tree_recursive
выглядит следующим образом:
Общая сложность: $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\]Если вам действительно хочется построить оптимальное (или хотя бы очень близкое к оптимальному) дерево, то на сегодня для решения этой проблемы не нужно придумывать кучу эвристик самостоятельно, а можно воспользоваться специальными солверами, которые решают NP-полные задачи приближённо, но всё-таки почти точно. Так что единственной (и вполне решаемой) проблемой будет представить исходную задачу в понятном для солвера виде. Пример построения оптимального дерева с помощью решения задачи целочисленного программирования.
Как вы, может быть, уже заметили, решающие деревья — это одна большая эвристика для решения NP-полной задачи, практически лишённая какой-либо стройной теоретической подоплёки. В 1970–1990-e годы интерес к ним был весьма велик как в индустрии, где был полезен хорошо интерпретируемый классификатор, так и в науке, где учёные интересовались способами приближённого решения NP-полных задач.
В связи с этим сложилось много хорошо работающих наборов эвристик, у которых даже были имена: например, ID3 был первой реализацией дерева, минимизирующего энтропию, а CART — первым деревом для регрессии. Некоторые из них были запатентованы и распространялись коммерчески.
На сегодня это всё потеряло актуальность в связи с тем, что существуют хорошо написанные библиотеки (например, sklearn, в которой реализована оптимизированная версия CART).