CatBoost: от формул к реализации на Python и в Excel
CatBoost: от формул к реализации на Python и в Excel

О чем эта статья

Привет! Меня зовут Федор Тарасов, я занимаюсь машинным обучением в команде collection стрима моделирования розничного бизнеса в Департаменте анализа данных и моделирования банка ВТБ.

Наша команда занимается разработкой ML-моделей для различных процессов взыскания — модели предсказывают вероятность выхода клиента в просрочку, прогнозируют последовательное углубление просрочки при невнесении или неполном внесении очередных платежей по кредиту, а также иные исходы, включая банкротство клиента. Каждая такая ситуация может быть сформулирована как задача классификации, и я на простых примерах хочу показать, как обучается классификатор CatBoost — отечественная библиотека градиентного бустинга, которую мы активно используем в работе. Задача реализовать этот алгоритм вручную интересовала меня с самого начала моей карьеры. Долгое время я шёл к ней окольными путями: цель осложнялась внутренними оптимизациями библиотеки и сложностями в чтении кода на C++, на котором реализованы её вычисления.

В итоге я написал учебную реализацию классификатора CatBoost на Python, подробно иллюстрирующую всю логику его работы, и выложил код на GitHub. Математический блок во многом опирается на лекцию YouTube-канала StatQuest with Josh Starmer, которую я адаптировал под специфику CatBoost. Для упрощения мы не будем использовать категориальные переменные, отключим случайные подвыборки (bootstrap) и механизм внесения случайности в скор сплитов (random_strength). Статья рассчитана на читателей, уже знакомых с основами градиентного бустинга: его назначением и общей идеей.

Может возникнуть вопрос: а зачем вообще разбирать алгоритм по косточкам, если он уже реализован и отлажен? Эта статья предназначена для тех, кто хочет разобраться во всех деталях работы CatBoost и не боится большого количества математических формул.

В этой статье мы:

  • Выведем целевую функцию обучения из логарифмического правдоподобия (Logloss).

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

  • Поймём, зачем нужна предварительная сетка квантизации признаков и как уникальная архитектура симметричных деревьев CatBoost спасает от переобучения.

  • И самое главное — вооружившись Excel-таблицей и Python-скриптом, шаг за шагом воспроизведём на небольшом примере логику вычислений индустриального фреймворка и сойдёмся с его ответом.

Начнём с постановки задачи.

Математика классификатора градиентного бустинга

Представим, что мы решаем задачу предсказания факта банкротства физического лица, являющегося заемщиком по кредиту. Для каждого наблюдения нам нужно получить не просто бинарный ответ «обанкротится / не обанкротится», а вероятность p ∈ [0, 1]. Мы хотим, чтобы эта вероятность была близка к 1 для банкротов (класс 1) и близка к 0 для платёжеспособных клиентов (класс 0). Это позволит не просто получить прогноз, но отранжировать клиентов по степени риска. В терминах машинного обучения рассматриваемая проблема формулируется как задача бинарной классификации.

В машинном обучении эта задача решается через оптимизацию функции потерь. В CatBoost для бинарной классификации по умолчанию используется логарифмическая функция правдоподобия (Log-Likelihood). Её максимизация (эквивалентная минимизации логарифмической функции потерь, LogLoss) будет приближать нас к озвученной выше цели.

LL = \sum_{i=1}^{N} y_i \times \log(p_i) + (1-y_i) \times \log(1-p_i)

где p — это pred модели, y — это 0 или 1 класс таргета. Следует отметить, что таких функций придумано много для разных задач, но наиболее популярна именно эта. Давайте убедимся, что эту функцию действительно нужно максимизировать, и посмотрим, как именно p в ней связана с нашим таргетом.

Рассмотрим для позитивного класса (конечно, это не означает, что мы рассматриваем событие банкротства клиента как позитивное и тем более желательное, таковы лишь особенности устоявшейся терминологии машинного обучения) y=1, правая часть функции при домножении на 1-y, при y=1 зануляется. Для левой части график функции 1*\ln(p) максимизируется при p → 1

Рисунок 1. Логарифмическое правдоподобие для класса 1
Рисунок 1. Логарифмическое правдоподобие для класса 1

Рассмотрим для негативного класса y=0, тогда левая часть функции при домножении на y, при y=0 зануляется. Для правой части график функции (1 - 0)*\ln(1 - p) максимизируется при p → 0.

Рисунок 2. Логарифмическое правдоподобие для класса 0
Рисунок 2. Логарифмическое правдоподобие для класса 0

В реальности у нас будет комбинация 1 и 0, которая разложится на сумму слагаемых, где каждое слагаемое, если мы правильно предскажем класс и его вероятность, будет максимизировать нашу функцию LL.

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

LL = \sum_{i=1}^{N} y_i \times \log(p_i) + (1-y_i) \times \log(1-p_i)

В машинном обучении принято решать задачи оптимизации через поиск минимума. Алгоритму нужна функция потерь (Loss Function, штраф за ошибку), по поверхности которой он будет «спускаться», последовательно уменьшая ошибку. Чтобы превратить задачу максимизации в задачу минимизации, правдоподобие просто умножают на - 1. Так появляется Отрицательное Логарифмическое Правдоподобие (Negative Log-Likelihood, NLL) или Logloss (параметр Loss по умолчанию в CatBoostClassifier):

LogLoss = - LL

Градиентный бустинг не предсказывает вероятность от 0 до 1. Он предсказывает неограниченные сырые значения — логиты (log-odds), обозначим их z_{i} \in ( - \infty, + \infty). Связь между логитом z_{i} и вероятностью p_{i} для i-го объекта выборки задается функцией сигмоиды и обратным ей логарифмическим преобразованием:

p_i = \frac{1}{1+e^{-z_i}} = \frac{e^{z_i}}{e^{z_i}+1} \Leftrightarrow z_i = \ln\left(\frac{p_i}{1-p_i}\right)

Выразим LL через логиты. Подставим это в логарифмы:

\ln(p_i)=\ln(e^{z_i})-\ln(1+e^{z_i})=z_i-\ln(1+e^{z_i})\ln(1-p_i)=\ln(1)-\ln(1+e^{z_i})=0-\ln(1+e^{z_i})

Теперь подставим это в {LL}_{i}

{LL}_{i} = y_{i}\left( z_{i} - \ln(1 + e^{z_{i}}) \right) + (1 - y_{i})\left( - \ln(1 + e^{z_{i}}) \right)

Раскроем скобки:

{LL}_{i} = y_{i}z_{i} - y_{i}\ln(1 + e^{z_{i}}) - \ln(1 + e^{z_{i}}) + y_{i}\ln(1 + e^{z_{i}})

Слагаемые y_{i}\ln(1 + e^{z_{i}}) взаимно уничтожаются, и мы получаем LL через логиты:

LL_i = y_i z_i - \ln(1+e^{z_i}), \qquad -LL_i=-y_i z_i+\ln(1+e^{z_i})

Обозначим индивидуальную потерю L(y_i,z_i)=-LL_i. Заранее найдём первую (gradient — g) и вторую (hessian — h) производные от - {LL}_{i}:

g_{i} = \frac{\partial L}{\partial z_{i}} = - y_{i} + \frac{1}{1 + e^{z_{i}}} \cdot e^{z_{i}} = - y_{i} + p_{i}h_{i} = \frac{\partial^{2}L}{\partial z_{i}^{2}} = p_{i}(1 - p_{i})

Градиентный бустинг обучается пошагово, каждое последующее дерево корректирует предсказанный результат предыдущих. Предположим, у нас есть для наблюдения i предсказание z_{i}^{t - 1} (для нулевого шага p=0.5, z_{i} = 0). Мы обучим дерево t, которое минимизирует общую функцию потерь. Дерево t выдаст для наблюдения i значение O_{value}^{t} (сырое логит значение листа, куда попадет наблюдение). Значит, новое предсказание будет равно: z_{i}^{t} = z_{i}^{t - 1} + O_{value}^{t}. Мы хотим найти такое O_{value}^{t}, которое минимизирует ошибку для всех объектов, попавших в лист. Но чтобы модель не переобучалась (не выдавала в листьях бесконечно огромные числа, которые превращают вероятности строго в 1 или 0), мы добавляем L2-регуляризацию (штраф) за слишком большое значение листа.

Именно так мы получаем итоговую целевую функцию обучения для одного листа дерева t, складывающуюся из суммы n элементов этого листа. Пусть {LL}_{t} — суммарная ошибка всех объектов, которые попали в лист j дерева t, тогда

{LL}_{t} = \sum_{i = 1}^{n}L(y_{i},z_{i}^{t - 1} + O_{value}^{t}) + \frac{1}{2}\lambda{O_{value}^{t}}^{2}

L(y_{i},z_{i}) — функция ошибки от логита z и таргета y. Для оптимизации используем общую формулу Тейлора: f(x + \Delta x) \approx f(x) + f'(x)\Delta x + \frac{1}{2}f''(x)\Delta x^{2}. Применяем её к нашей функции потерь, где роль малого приращения \Delta x играет предсказание нового листа O_{value}^{t}:

\begin{aligned}L(y_i,z_i^{t-1}+O_{value}^t)&\approx L(y_i,z_i^{t-1})\\&\quad+\underbrace{\left[\frac{d}{dz_i}L(y_i,z_i)\right]}_{g_i}O_{value}^t\\&\quad+\frac12\underbrace{\left[\frac{d^2}{dz_i^2}L(y_i,z_i)\right]}_{h_i}(O_{value}^t)^2\end{aligned}

Распишем сумму {LL}_{t} для всех n наблюдений, попавших в конкретный лист:

\begin{aligned}LL_t&=L(y_1,z_1^{t-1}+O_{value}^t)\\&\quad+L(y_2,z_2^{t-1}+O_{value}^t)+\ldots\\&\quad+L(y_n,z_n^{t-1}+O_{value}^t)+\frac12\lambda(O_{value}^t)^2\end{aligned}

Подставляем в неё полученную аппроксимацию Тейлора:

\begin{aligned}LL_t&=L(y_1,z_1^{t-1})+g_1O_{value}^t+\frac12h_1(O_{value}^t)^2\\&\quad+L(y_2,z_2^{t-1})+g_2O_{value}^t+\frac12h_2(O_{value}^t)^2\\&\quad+\ldots+L(y_n,z_n^{t-1})+g_nO_{value}^t\\&\quad+\frac12h_n(O_{value}^t)^2+\frac12\lambda(O_{value}^t)^2\end{aligned}

Мы хотим найти такое выходное значение O_{value}^{t}, которое минимизирует функцию потерь. Для поиска минимума берем производную LL по O_{value} (для упрощения не буду дальше писать t) и приравниваем её к нулю, \underset{0}{\underbrace{\text{зануляя}}} константы:

\begin{aligned}\frac{dLL_t}{dO_{value}}&=\underbrace{\frac{d\left[\sum_{i=1}^{n}L(y_i,z_i^{t-1})\right]}{dO_{value}}}_{0}\\&\quad+\frac{d}{dO_{value}}\left[\begin{aligned}&g_1O_{value}+\frac12h_1O_{value}^2\\&+g_2O_{value}+\frac12h_2O_{value}^2+\ldots\\&+g_nO_{value}+\frac12h_nO_{value}^2+\frac12\lambda O_{value}^2\end{aligned}\right]=0\end{aligned}

Группируем слагаемые с O_{value} и O_{value}^{2}:

\frac{{dLL}_{t}}{{dO}_{value}} = \frac{d\left\lbrack O_{value}(g_{1} + g_{2} + \ldots + g_{n}) + \frac{1}{2}O_{value}^{2}(h_{1} + h_{2} + \ldots + h_{n} + \lambda) \right\rbrack}{{dO}_{value}} = 0

Вычисляем саму производную:

(g_{1} + g_{2} + \ldots + g_{n}) + O_{value}(h_{1} + h_{2} + \ldots + h_{n} + \lambda) = 0

Выражаем O_{value}:

O_{value}^{*}=\frac{-(g_1+g_2+\dots+g_n)}{h_1+h_2+\dots+h_n+\lambda}\quad \left\{\begin{array}{l}g_{i} = (p_{i} - y_{i}) \\h_{i} = p_{i}(1 - p_{i})\end{array} \right.\

Для обучения дерева используем антиградиент -g_{i} = y_{i} - p_{i}: минус перед числителем дроби можно внести под знак суммы.

О знаках и обозначениях. В выводе через ряд Тейлора и в следующих формулах ответа листа и Gain обозначение g_i = p_i - y_i означает градиент функции потерь; направление уменьшения ошибки задаёт антиградиент -g_i = y_i - p_i. Ниже, в записи L2-скора из документации и в формулах косинусного скора, буквой g_i обозначается уже антиградиент. Это смена соглашения об обозначениях: дополнительный минус в этих формулах не нужен.

Если переписать ответ листа в терминах остатков (Residuals) и вероятностей:

O_{value}^{*}=\frac{\sum_i \mathrm{Residual}_i}{\sum_i[\mathrm{PreviousProbability}_i(1-\mathrm{PreviousProbability}_i)]+\lambda}

Итак, мы вывели O_{value}^{*} — идеальное число, которое нужно поместить в лист дерева, чтобы минимизировать ошибку для попавших туда объектов. Но возникает вопрос: как алгоритм изначально понимает, по какому признаку разделить данные на эти листья? (Например, выбрать правило “Возраст > 30” или “Доход < 50k”?). Для этого алгоритму нужна метрика, чтобы оценивать «полезность» каждого возможного разреза. Эта метрика называется Score (или Gain / Выигрыш). Логика её получения: давайте возьмем найденный нами идеальный ответ O_{value}^{*} и подставим его обратно в аппроксимацию Тейлора. Мы хотим узнать: насколько сильно уменьшится наша общая ошибка, если мы создадим этот лист? Чтобы это понять, давайте вернемся к нашей развернутой функции потерь {LL}_{t} (аппроксимированной по Тейлору) и просто сгруппируем в ней слагаемые.

Сумма старых ошибок \sum L(y_{i},z_{i}^{t - 1}) для текущего шага является константой (ошибку предыдущих деревьев мы изменить уже не можем). А вот слагаемые с O_{value} и O_{value}^{2} можно сгруппировать, вынеся эти множители за скобки:

{LL}_{t} = \underset{\text{Ошибка предыдущих деревьев}}{\underbrace{\sum_{i = 1}^{n}L(y_{i},z_{i}^{t - 1})}} + \underset{\text{Изменение ошибки благодаря новому листу}}{\underbrace{\left( \sum_{i=1}^{n} g_i \right)O_{value} + \frac{1}{2}\left( \sum_{i=1}^{n} h_i + \lambda \right)O_{value}^{2}}}\,

Чтобы оценить качество сплита, алгоритму не нужна полная ошибка, ему достаточно знать её изменение. Отбросим константу и оставим только ту часть функции, которая зависит от нашего предсказания. Обозначим это изменение как \Delta LL:

\Delta LL = \left( \sum_i g_i \right)O_{value} + \frac{1}{2}\left( \sum_i h_i + \lambda \right)O_{value}^{2}

Мы хотим узнать: какое минимальное значение примет эта функция (насколько сильно мы сможем уменьшить ошибку), если поместим в лист наш идеальный ответ O_{value}^{*}? Для этого подставим выведенный O_{value}^{*} обратно в \Delta LL. Полученное оптимальное (минимальное) изменение ошибки для листа в литературе по бустингу традиционно обозначают как \mathcal{L}_{\text{leaf}}^{*} (где звездочка указывает на то, что значение вычислено в точке оптимума). Для краткости обозначим суммы градиентов и гессианов в листе как S_{g} = \sum_i g_i и S_{h} = \sum_i h_i:

\mathcal{L}_{\text{leaf}}^{*} = S_{g}\left( - \frac{S_{g}}{S_{h} + \lambda} \right) + \frac{1}{2}(S_{h} + \lambda)\left( - \frac{S_{g}}{S_{h} + \lambda} \right)^{2}

Давайте раскроем скобки и упростим выражение. Возведем третий множитель во втором слагаемом в квадрат (минус при этом уйдет):

\mathcal{L}_{\text{leaf}}^{*} = - \frac{S_{g}^{2}}{S_{h} + \lambda} + \frac{1}{2}(S_{h} + \lambda)\frac{S_{g}^{2}}{(S_{h} + \lambda)^{2}}

Сократим (S_{h} + \lambda) во втором слагаемом:

\mathcal{L}_{\text{leaf}}^{*} = - \frac{S_{g}^{2}}{S_{h} + \lambda} + \frac{1}{2}\frac{S_{g}^{2}}{S_{h} + \lambda}

Обратите внимание: мы получили выражение вида «- 1 + 0.5». Приводим подобные слагаемые и получаем итоговое минимальное значение изменения функции потерь для листа:

\mathcal{L}_{\text{leaf}}^{*} = - \frac{1}{2}\frac{S_{g}^{2}}{S_{h} + \lambda}

При

S_h+\lambda>0 это значение неположительно; если S_g=0, оно равно нулю. Оно показывает изменение квадратичной аппроксимации, а не точное изменение исходного Logloss.

При построении дерева алгоритм перебирает разные варианты сплитов (разбиений) и для каждого оценивает выигрыш. Так как алгоритму удобнее максимизировать положительную величину, знак минус отбрасывают (как часто опускают и константу \frac{1}{2}, которая одинакова для всех сплитов и не влияет на их ранжирование). Так получается ньютоновская оценка полезности листа, с точностью до постоянного множителя:

\mathrm{Score}_{\mathrm{leaf}}=\frac{(\sum_i g_i)^2}{\sum_i h_i+\lambda}

В отличие от классических алгоритмов (XGBoost, LightGBM), которые разделяют каждый узел независимо, CatBoost строит симметричные (oblivious) деревья. Это значит, что один признак и один порог применяются глобально — сразу ко всему уровню дерева. Поэтому алгоритм считает выигрыш (Gain) не локально для пары «левый-правый лист», а суммарно для всего уровня.

Пусть на текущем уровне K листьев. Тогда суммарный выигрыш вычисляется как:

{Gain}_{level} = \sum_{j = 1}^{K}{}\left( \frac{(\sum_{i \in \text{leaf}_j} g_i)^{2}}{\sum_{i \in \text{leaf}_j} h_i + \lambda} \right)

(При поиске лучшего сплита множитель \frac{1}{2} отбрасывается, так как он является константой и не влияет на то, какой вариант разбиения окажется максимальным).

Дерево перебирает сотни признаков и порогов, вычисляет эту метрику и выбирает тот сплит, где {Gain}_{level} максимален. Но давайте посмотрим на эту формулу под другим углом. Вспомним выведенное нами ранее идеальное значение предсказания. Обозначим его для конкретного j-го листа как O_{\mathrm{value},j}^{*}\,:

O_{\mathrm{value},j}^{*}\, = \frac{\sum_{i \in \text{leaf}_{j}} - g_{i}}{\sum_{i \in \text{leaf}_j} h_i + \lambda}

Если вынести O_{\mathrm{value},j} за скобки в формуле {Gain}_{level} (уже без \frac{1}{2}), мы увидим закономерность:

{Gain}_{level} = \sum_{j = 1}^{K}{}\left( O_{\mathrm{value},j}^{*}\, \cdot \sum_{i \in \text{leaf}_{j}} - g_{i} \right)

Если развернуть эту сумму до уровня отдельных наблюдений, окажется, что это не что иное, как скалярное произведение вектора антиградиентов \mathbf{G} (где каждый элемент вектора — это значение антиградиента - g_{i} = y_{i} - p_{i}, то есть остаток для i-го наблюдения) и вектора предсказаний дерева \mathbf{A} (где каждый элемент вектора — это ответ O_{\mathrm{value},j} того листа, куда попал объект). Таким образом (с точностью до константы ½, которую мы отбросили при сравнении сплитов):

{Gain}_{level} \propto \mathbf{G} \cdot \mathbf{A}

Кандидат для скора и итоговый ответ листа

Выше через разложение Тейлора мы получили ньютоновскую оценку полезности листа и его ответ O_{\mathrm{value},j}^{*}. Эта оценка относится к квадратичной аппроксимации функции потерь; найденный ответ задаёт один шаг Ньютона. Для поиска сплитов с L2 или Cosine в рассматриваемом CPU-режиме Plain CatBoost использует другую величину. Ниже рассматриваем числовые признаки, единичные веса объектов, симметричные деревья и один ньютоновский шаг окончательной оценки листьев.

Чтобы не смешивать соглашения о знаке, здесь обозначим остаток через r_i=y_i-p_i. В выводе Тейлора выше это антиградиент r_i=-g_i, поскольку там g_i=p_i-y_i. В записи L2-скора из документации и в формулах Cosine ниже буква g_i уже обозначает сам антиградиент, то есть тот же r_i. Дополнительный минус там не нужен.

Пусть n_j — число объектов в листе j. Кандидат для расчёта L2/Cosine-скора равен:

a_j=\frac{\sum_{i\in\mathrm{leaf}_j}r_i}{n_j+\lambda}

После выбора структуры дерева окончательный ньютоновский ответ листа вычисляется с настоящими гессианами Logloss:

O_{\mathrm{value},j}^{*}=\frac{\sum_{i\in\mathrm{leaf}_j}r_i}{\sum_{i\in\mathrm{leaf}_j}p_i(1-p_i)+\lambda}

Таким образом, кандидат a_j и итоговый ответ O_{\mathrm{value},j}^{*} могут различаться. В первом случае в знаменателе число объектов, во втором — сумма гессианов. Для L2-скора в этом режиме реализация суммирует произведения кандидата и суммы антиградиентов:

\mathrm{Score}_{L2}=\sum_j a_j\sum_{i\in\mathrm{leaf}_j}r_i

Именно знаменатель n_j+\lambda используется в основной функции calculate_score учебного ноутбука. Сумма гессианов нужна в update_predictions, после выбора структуры дерева. Это соответствует исходникам CatBoost 1.2.8: score_calcers.cpp и short_vector_ops.h. В разделе «Упрощение формулы для вычислений под капотом» вектор \mathbf A будет составлен из кандидатов a_j, повторённых для объектов соответствующего листа. Приведённые выше формулы с O_{\mathrm{value},j}^{*} описывают ньютоновскую оценку, а не штатный L2/Cosine-скор этого режима.

Как читать следующий эксперимент. В разделе «Геометрический эксперимент: чувствительность к масштабу» и соответствующем опыте ноутбука кандидат намеренно заменён ньютоновской дробью с суммой гессианов в знаменателе. Названия L2 и Cosine в этом опыте относятся к ненормированному скалярному произведению и его косинусной нормировке для изменённого кандидата. Обсуждаемое там «взрывное» поведение относится к этой модификации учебного кода. В проверке штатного CatBoost 1.2.8 с теми же начальными логитами оба скора выбрали X1 при обоих значениях регуляризации. Отдельной тестовой выборки в опыте нет, поэтому он показывает влияние масштаба кандидатов, но не доказывает преимущество одного критерия на новых данных.

В документации CatBoost для ‘score_function’ L2 = - \sum_{i} w_{i} \cdot (a_{i} - g_{i})^{2}. Здесь g_i обозначает антиградиент. Для оптимальных ответов листьев без регуляризации этот скор после раскрытия квадрата и отбрасывания константы тоже выражается через скалярное произведение. Условия сопоставления с формулой выше разобраны в спойлере.

Вывод эквивалентности L2-скора из документации

Начнем с формулы из документации: Score = - \sum_{i} w_{i}(a_{i} - g_{i})^{2}. Раскроем квадрат и отбросим константу \sum_{i} w_{i}g_{i}^{2} (она не зависит от выбора сплита):

Score \propto \sum_{i}(2w_{i}a_{i}g_{i} - w_{i}a_{i}^{2})

Сгруппируем по листьям j, где предсказание a_{j} постоянно:

Score \propto \sum_{j = 1}^{K}{}\left( 2a_{j}G_{j} - a_{j}^{2}W_{j} \right)

(где G_{j} = \sum_{i \in j} w_{i}g_{i} — сумма антиградиентов, W_{j} = \sum_{i \in j} w_{i}\, — сумма весов). Найдем оптимальное значение листа a_{j}^{*}, максимизирующее это выражение (приравняв производную к нулю):

a_{j}^{*} = \frac{G_{j}}{W_{j}}

Условия этого вывода. Здесь мы оптимизируем L2-скор без регуляризации. Полученный ответ листа совпадает с ранее выведенным ньютоновским ответом при h_i = 1 и \lambda = 0, если в обеих формулах одинаково учтены веса объектов. Для приведённой выше невзвешенной формулы это означает w_i = 1. В численном примере ниже задано reg_lambda = 3.0, поэтому напрямую переносить на него это равенство нельзя: регуляризацию нужно учитывать отдельно.

Подставим a_{j}^{*} обратно в формулу скора:

{Score}^{*} \propto \sum_{j = 1}^{K}{}\left( 2\frac{G_{j}}{W_{j}}G_{j} - \left( \frac{G_{j}}{W_{j}} \right)^{2}W_{j} \right) = \sum_{j = 1}^{K}\frac{G_{j}^{2}}{W_{j}}

Раскроем скалярную структуру. Заметим, что \frac{G_{j}^{2}}{W_{j}} = \left( \frac{G_{j}}{W_{j}} \right) \cdot G_{j} = a_{j}^{*} \cdot G_{j}. Следовательно:

{Score}^{*} \propto \sum_{j = 1}^{K}a_{j}^{*}G_{j} = \mathbf{A} \cdot \mathbf{G}

Мы получили, что выигрыш от сплита пропорционален скалярному произведению вектора оптимальных предсказаний \mathbf{A} и вектора антиградиентов \mathbf{G}.

Геометрический эксперимент: чувствительность к масштабу

Скалярное произведение зависит и от угла между векторами, и от их длин:

\mathbf G\cdot\mathbf A=\|\mathbf G\|\,\|\mathbf A\|\cos\theta

Чтобы отдельно показать влияние масштаба, проведём учебный эксперимент. В разделе ноутбука «L2 Score VS Cosine» мы намеренно заменяем кандидата a_j на ньютоновскую дробь с суммой гессианов в знаменателе. Сравним ненормированный критерий \mathbf G\cdot\mathbf A и его косинусную нормировку. Это изменение формулы учебного алгоритма, а не воспроизведение штатного L2-скора CatBoost.

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

# В отдельном экспериментальном датасете последние три строки — объекты класса 1.
df_experiment['pred_prob'] = 0.5
df_experiment.iloc[-3:, df_experiment.columns.get_loc('pred_prob')] = 0.001

При l2_reg_experiment = 0.001 изменённый ненормированный критерий выбирает X2 ≤ 4,37637, изолируя три объекта. Косинусный критерий выбирает X1 ≤ −0,20326, разделяя основные кластеры:

Рисунок 3. Сравнение разбиений L2 и Cosine при экстремальных вероятностях
Рисунок 3. Сравнение разбиений L2 и Cosine при экстремальных вероятностях

Открыть рисунок 3 в полном размере

При l2_reg_experiment = 0.1 в том же эксперименте оба критерия выбирают X1 ≤ −0,20326:

Рисунок 4. Сравнение разбиений L2 и Cosine при l2_reg_experiment = 0.1
Рисунок 4. Сравнение разбиений L2 и Cosine при l2_reg_experiment = 0.1

Открыть рисунок 4 в полном размере

Этот опыт показывает чувствительность выбранного критерия к масштабу кандидатов. Он не доказывает, что один метод лучше обобщает на новые данные: отдельной тестовой выборки здесь нет. При проверке настоящего CatBoost 1.2.8 с теми же начальными логитами оба штатных скора, L2 и Cosine, выбрали X1 при обоих значениях регуляризации.

Геометрический переход к Cosine

Вернёмся к обычному кандидату a_j=\sum_i r_i/(n_j+\lambda) для поиска сплитов. Косинусное сходство нормирует скалярное произведение на длины векторов:

\cos(\mathbf G,\mathbf A)=\frac{\mathbf G\cdot\mathbf A}{\|\mathbf G\|\,\|\mathbf A\|}

Умножение всего вектора \mathbf A на положительную константу не меняет этот косинус. Изменение только одного листа может менять и направление вектора, поэтому косинусная нормировка сама по себе не гарантирует защиту от выбросов или переобучения.

Упрощение формулы для вычислений под капотом

Давайте адаптируем эту формулу для расчетов в коде. Представим все объекты выборки как точки в N-мерном пространстве (где N — число строк в датасете):

Вектор антиградиентов \mathbf{G} = g_{1},g_{2},\ldots,g_{N}, где в этом разделе g_i = y_i - p_i. Это то направление, куда нам нужно двигаться для минимизации ошибки.

Вектор предсказаний \mathbf{A} = a_{1},a_{2},\ldots,a_{N}. Здесь a_i — кандидат a_j для листа, в который попал объект: a_j=\sum_{i\in\mathrm{leaf}_j}g_i/(n_j+\lambda). Напомним, что в этом разделе g_i означает антиградиент.

Косинус угла между этими векторами вычисляется как:

\cos(\mathbf{G},\mathbf{A})=\frac{\sum_{i=1}^{N}g_i a_i}{\sqrt{\sum_{i=1}^{N}g_i^2}\cdot\sqrt{\sum_{i=1}^{N}a_i^2}}

При поиске лучшего сплита внутри одной итерации алгоритм проводит упрощения:

Отбрасывание константы. Знаменатель содержит длину вектора антиградиентов \sqrt{\sum_{i = 1}^{N}g_{i}^{2}}. Эта величина зависит только от ошибок предыдущих деревьев и абсолютно неизменна для всех проверяемых вариантов текущего сплита. При поиске максимума ее можно смело отбросить.

Группировка числителя по листьям. Числитель \sum_{i = 1}^{N}g_{i}a_{i} — это скалярное произведение. Поскольку значения a_{i} одинаковы для всех объектов внутри одного листа j (и равны a_j), мы можем перегруппировать сумму из N объектов в сумму по K листьям:

\sum_{i = 1}^{N}g_{i}a_{i} = \sum_{j = 1}^{K}{}\left( a_j \cdot \sum_{i \in \text{leaf}_j} g_i \right)

Это L2-скор рассматриваемого CPU/Plain режима. Он имеет ту же форму скалярного произведения, что и ньютоновская оценка выше, но использует другие кандидаты листьев.

Группировка знаменателя по листьям. Оставшаяся часть знаменателя — это норма вектора предсказаний. Мы также группируем ее по листьям. Если n_{j} — количество объектов в листе j, то сумма квадратов одинаковых предсказаний внутри листа сворачивается в простое умножение:

\sqrt{\sum_{i=1}^{N}a_i^2}=\sqrt{\sum_{j=1}^{K}n_j\cdot(a_j)^2}

Сводя всё воедино, мы получаем итоговую, вычислительно эффективную метрику для используемого CPU/Plain режима CatBoost. Она отличается от полного косинуса множителем, одинаковым для всех кандидатов текущего уровня:

\text{Score}_{\text{Cosine}}=\frac{\sum_{j=1}^{K}\left(a_j\cdot\sum_{i\in\text{leaf}_j}g_i\right)}{\sqrt{\sum_{j=1}^{K}n_j\cdot(a_j)^2}}

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

Одна из фишек CatBoost — это симметричные деревья. Это значит, что выбранное правило сплита применяется глобально и одновременно ко всем листьям текущего яруса дерева.

  • В корне 1 лист режется на 2.

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

  • На следующем уровне новое правило разрежет каждый из 4 листьев, создав 8. И так далее.

Структура дерева выбрана. Теперь вычисляем ньютоновские ответы O_{value} для его листьев. Казалось бы, пора обновить предсказание: z_{i}^{t} = z_{i}^{t - 1} + O_{value}^{t}.

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

Бустинг — это процесс осторожного приближения к цели. Поэтому вводится параметр Learning Rate (скорость обучения, lr) — небольшое число (например, 0.03 или 0.1). Например, 0.03 и 0.1 соответствуют 3% и 10% рассчитанного обновления:

z_{i}^{t} = z_{i}^{t - 1} + (lr \times O_{value}^{t})

(На диаграммах и в сохранённой модели CatBoost 1.2.8 значения листьев уже учитывают learning_rate. В наших формулах O_value обозначает ответ до этого умножения. При суммировании готовых значений листьев библиотеки повторно умножать их на learning_rate не нужно.)

Как модель делает финальное предсказание

Давайте соберем весь пазл: как делается предсказание для новых данных?

Модель начинает с базового логита. В нашем примере boost_from_average=False, поэтому z_0=0, что соответствует вероятности 0,5. При другом способе инициализации базовый логит может отличаться.

Информация о новом наблюдении “прокатывается” через все M обученных деревьев. В каждом дереве наблюдение падает в какой-то лист.

Мы суммируем предсказания всех этих деревьев:

Z_{final} = z_{0} + (lr \times O_{value}^{1}) + (lr \times O_{value}^{2}) + \ldots + (lr \times O_{value}^{M})

В результате мы получаем итоговый сырой логит Z_{final} (например, 1.85).

Последний шаг — пропускаем логит через сигмоиду (с которой начинали статью):

p = \frac{1}{1 + e^{- Z_{final}}}

Построение CatBoost на игрушечном примере

Проиллюстрируем это на примере. У CatBoost есть встроенный функционал для визуализации деревьев. Для отображения результата plot_tree нужны Python-пакет graphviz и установленная системная программа dot.

model_toy.plot_tree(
    tree_idx=0,
    pool=pool_toy
)
Рисунок 5. Визуализация дерева CatBoost на игрушечном датасете
Рисунок 5. Визуализация дерева CatBoost на игрушечном датасете

Открыть рисунок 5 в полном размере

Далее мы постараемся вручную без кода построить классификатор и выйти на предсказания, которые дает нам CatBoost из коробки. У нас есть игрушечный датасет (полностью выдуманный):

person_income

loan_percent_income

Target

Catboost_pred

100000

0.10

1

0.675124

80000

0.15

1

0.659438

60000

0.05

1

0.657338

120000

0.18

1

0.677169

50000

0.12

1

0.657338

150000

0.08

1

0.675124

30000

0.40

0

0.340116

45000

0.35

0

0.358294

80000

0.50

0

0.340116

25000

0.45

0

0.340116

70000

0.55

0

0.340116

90000

0.25

1

0.677169

110000

0.28

1

0.613931

35000

0.26

0

0.395619

40000

0.29

0

0.358294

Датасет состоит из 15 клиентов, у которых 1 — это возврат займа, а 0 — выход в просрочку (обратите внимание: в отличие от примера с банкротством из начала статьи, здесь класс 1 — благополучный исход; смысл таргета противоположный). Person Income и Loan Percent Income — наши объясняющие переменные. Catboost_pred — это прогнозы реальной модели, обученной кодом:

target = 'loan_status'
features = ['loan_percent_income', 'person_income'] #, 'person_age'
cb_params = {
  'iterations': 2,
  'depth': 4,
  'learning_rate': 0.8,
  'random_seed': 1,
  'reg_lambda': 3.0,
  'bootstrap_type': 'No',
  'penalties_coefficient': 0,
  'langevin': False,
  'diffusion_temperature': 0,
  'random_strength': 0,
  'logging_level': 'Silent',
  'score_function': 'Cosine',
  'leaf_estimation_method': 'Newton'
}
pool = Pool(df[features], df[target])
pool.quantize()
pool.save_quantization_borders('border1.txt')
bins_df = pd.read_csv('border1.txt', sep='\t', header=None,index_col=0)
log_buffer = io.StringIO()
model = CatBoostClassifier(**cb_params).fit(pool,log_cout=log_buffer,logging_level='Debug')

Заглянем внутрь процесса обучения модели:

print(log_buffer.getvalue())

Вернет:

loan_percent_income, bin=6 score 1.695285187
person_income, bin=8 score 1.896735143
loan_percent_income, bin=10 score 1.9139964
loan_percent_income, bin=6 score 1.9139964
tensor 0 is redundant, remove it and stop
0: learn: 0.5341393 total: 235us remaining: 235us

loan_percent_income, bin=8 score 1.391268638
person_income, bin=3 score 1.557835675
loan_percent_income, bin=3 score 1.582888584
person_income, bin=3 score 1.582888584
tensor 1 is redundant, remove it and stop
1: learn: 0.4241715 total: 336us remaining: 0us

Можно увидеть выбранные CatBoost фичи и их сплиты и соответствующие этим фичам и сплитам скоры.

Рисунок 6. Границы квантизации, деревья и журнал обучения CatBoost
Рисунок 6. Границы квантизации, деревья и журнал обучения CatBoost

Открыть рисунок 6 в полном размере

Слева на картинке приведены значения из файла border1.txt — рассчитанные CatBoost границы квантизации для каждой переменной (loan_percent_income имеет индекс 0, person_income — 1). Справа диаграммы деревьев и относящиеся к ним вывод журнала отладки.

Что такое квантизация и откуда берутся границы сплитов

Обратите внимание на вывод лога: CatBoost пишет не конкретное значение признака, а bin=6 или bin=10. А в файле border1.txt лежат значения вроде 0.2549999952, близкие к 0.255.

Если бы у нас были миллионы наблюдений, алгоритму пришлось бы на каждом шаге перебирать абсолютно каждое уникальное значение признака в поисках идеального разреза. Это вычислительно неэффективно. Поэтому CatBoost (как и другие современные библиотеки бустинга) на этапе инициализации применяет квантизацию (binning) — он заранее разбивает непрерывные значения признаков на диапазоны (бины).

Когда мы вызвали метод pool.quantize() и сохранили border1.txt, алгоритм проанализировал наши данные и создал сетку границ квантизации (borders; в нашем запуске — метод GreedyLogSum), по которым алгоритм будет перебирать возможные сплиты. Оптимальный сплит затем выбирается из этих кандидатов. Например, для переменной loan_percent_income граница с индексом 6 составила 0.2549999952 (приблизительно 0.255). Если посмотреть в наш исходный датасет, это значение лежит аккурат между клиентами со значениями 0.250 и 0.260. Для person_income граница между 80 000 и 90 000 прошла ровно посередине — 85 000.

Перебирая при построении дерева не все сырые данные подряд, а только эту заранее вычисленную сетку границ, библиотека сокращает число рассматриваемых порогов. Индекс bin=6 в логах означает, что лучший Score показала граница с индексом 6, то есть седьмая: индексы начинаются с нуля. Поэтому в наших ручных расчетах в Excel мы будем искать оптимальный сплит не по самим данным, а проверять именно эти заранее подготовленные пороги.

Обучение модели в Excel

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

Рассчитаем скор кандидатов по loan_percent_income на первом уровне первого дерева. Оптимальный сплит будет в точке loan_percent_income≈0.255 со скором 1.695 (максимальным среди всех возможных сплитов как фичи loan_percent_income, так и person_income). В Excel показаны промежуточные расчёты; полный перебор порогов обоих признаков выполняет учебный код.

Рисунок 7. Первое дерево: расчёт скора первого сплита в Excel
Рисунок 7. Первое дерево: расчёт скора первого сплита в Excel

Открыть рисунок 7 в полном размере

Напомним различие этапов: в нашем CPU/Plain режиме для L2/Cosine-скора кандидат листа вычисляется со знаменателем n_j+\lambda. После выбора структуры окончательный ответ методом Ньютона использует \sum_{i\in\mathrm{leaf}_j}p_i(1-p_i)+\lambda. Это не универсальная замена гессиана на единицу для всех режимов CatBoost.

Найдем второй оптимальный сплит. Учитывая то, что у нас уже есть сплит loan_percent_income≈0.255, все данные разделяются согласно этому сплиту.

Рисунок 8. Первое дерево: расчёт скора второго сплита в Excel
Рисунок 8. Первое дерево: расчёт скора второго сплита в Excel

Открыть рисунок 8 в полном размере

Найдём третий оптимальный сплит.

Рисунок 9. Первое дерево: расчёт скора третьего сплита в Excel
Рисунок 9. Первое дерево: расчёт скора третьего сплита в Excel

Открыть рисунок 9 в полном размере

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

В этом примере следующий выбранный порог уже присутствует в дереве и не создаёт нового разбиения. CatBoost удаляет избыточный уровень и останавливает рост дерева; это видно по сообщению tensor ... is redundant. Глубина также ограничена параметром depth, а возможность дальнейшего разбиения зависит от данных и критерия.

Рассчитаем вероятность нашего первого дерева, обращаясь к формулам, выведенным ранее:

z_{1} = z_{0} + (lr \times O_{value}^{1})new\_ pred\_ prob = \frac{1}{1 + e^{- z_{1}}}
Рисунок 10. Первое дерево: значения листьев и прогноз вероятности
Рисунок 10. Первое дерево: значения листьев и прогноз вероятности

Открыть рисунок 10 в полном размере

Обучим второе дерево.

На первом уровне второго дерева оптимальный сплит будет в точке loan_percent_income≈0.285 со скором 1.391

Рисунок 11. Второе дерево: расчёт скора первого сплита в Excel
Рисунок 11. Второе дерево: расчёт скора первого сплита в Excel

Открыть рисунок 11 в полном размере

На втором уровне с учётом разбиения loan_percent_income≈0.285, максимальный скор 1.558 достигается при person_income=42 500.

Рисунок 12. Второе дерево: расчёт скора второго сплита в Excel
Рисунок 12. Второе дерево: расчёт скора второго сплита в Excel

Открыть рисунок 12 в полном размере

На третьем уровне с учётом предыдущих разбиений loan_percent_income≈0.285 и person_income=42 500 наибольший скор 1.5829 достигается при loan_percent_income=0.135

Рассчитаем вероятность нашего второго дерева, обращаясь к формулам, выведенным ранее:

z_{2} = z_{1} + (lr \times O_{value}^{2})new\_ pred\_ prob = \frac{1}{1 + e^{- z_{2}}}
Рисунок 13. Второе дерево: итоговые вероятности и сравнение с CatBoost
Рисунок 13. Второе дерево: итоговые вероятности и сравнение с CatBoost

Открыть рисунок 13 в полном размере

Для всех 15 объектов вероятности Excel совпадают с CatBoost 1.2.8 с максимальной абсолютной разницей около 2.2304\times10^{-9}. В таблице выше результаты округлены до шести знаков.

Обучение модели учебным кодом

Для иллюстрации всей логики обучения обратимся к ноутбуку с кодом.

Результат ниже проверен на CatBoost 1.2.8, NumPy 1.26.4 и pandas 2.2.3 под Windows, на CPU. В фактических параметрах модели: boosting_type=Plain, grow_policy=SymmetricTree, boost_from_average=False, leaf_estimation_iterations=1. Зададим параметры CatBoost модели.

cb_params_toy = {
    'iterations': 2,
    'depth': 4,
    'learning_rate': 0.8,
    'random_seed': 1,
    'reg_lambda': 3.0,
    'bootstrap_type': 'No',
    'penalties_coefficient': 0,
    'langevin': False,
    'diffusion_temperature': 0,
    'random_strength': 0,
    'logging_level': 'Silent',
    'score_function': 'Cosine',
    'leaf_estimation_method': 'Newton'
}

Основные параметры для нашего эксперимента:

  • score_function=‘Cosine’ — функция расчета скора. Учебный код поддерживает также ‘L2’ функцию.

  • bootstrap_type=No отключает bootstrap. Подвыборка признаков регулируется отдельно параметром rsm; в нашем запуске rsm=1, то есть участвуют все признаки.

  • depth — глубина деревьев. Определяет, сколько уровней сплитов будет в каждом дереве.

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

  • learning_rate — коэффициент, контролирующий величину шага нашего обучения на каждой итерации.

  • reg_lambda — коэффициент L2-регуляризации.

Квантируем наши данные. CatBoost произведет биннинг переменных, и далее мы будем работать с разбиненными значениями, разделенными на группы.

pool_toy = Pool(df_toy[features_toy], df_toy[target_toy])
pool_toy.quantize()
pool_toy.save_quantization_borders('border_toy.txt')
borders_dict_toy = parse_borders('border_toy.txt')

Обучим модель кодом CatBoost и нашим кодом

print("⏳ Обучение оригинального CatBoost...")
model_toy = CatBoostClassifier(**cb_params_toy).fit(pool_toy)
l2_reg_toy = model_toy.get_all_params()['l2_leaf_reg']

print("\n🚀 Запуск ручной реализации (Toy Dataset)...")
custom_pred_toy = fit_custom_catboost(
df_toy, features_toy, target_toy, borders_dict_toy,
cb_params_toy['iterations'], cb_params_toy['depth'],
cb_params_toy['learning_rate'], l2_reg_toy,score_function=cb_params_toy['score_function']
)

cb_pred_toy = model_toy.predict_proba(pool_toy)[:, 1]
max_diff_toy = (cb_pred_toy - custom_pred_toy).abs().max()

print("\n" + "="*80)
print(f"🎯 Максимальная разница вероятностей (Toy): {max_diff_toy:.12f}")
if max_diff_toy < 1e-7:
    print("✅ SUCCESS! Идеальное математическое совпадение.")
else:
    print("❌ Найдены расхождения.")
print("="*80)

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

Заключение

Градиентный бустинг часто воспринимается начинающими специалистами как некий магический «черный ящик»: мы просто скармливаем алгоритму датасет, вызываем метод .fit(), подбираем параметры по сетке и радуемся росту метрик. Нашей целью было открыть этот «черный ящик» и показать, что внутри нет никакой магии — лишь строгая, логичная и очень красивая математика.

Конечно, чтобы свести эту задачу к решаемой вручную, нам пришлось пойти на ряд осознанных упрощений. Мы отключили случайные подвыборки (bootstrap) и механизм внесения случайности в скор сплитов (random_strength). За кадром осталась и главная гордость библиотеки — работа с категориальными признаками на лету. Весь код, сопровождающий эту статью, можно найти в репозитории на GitHub. Успешных вам сплитов и высоких метрик!