xgboost. Теория и практика
- 1. Установка пакета xgboost
- 2. Градиентный бустинг: принцип метода и особенности реализации в xgboost
- 3. Использование xgboost совместно с caret
- 4. xgb.train: формат исходных данных и параметры модели
- 5. xgboost и Random Forest
- 6. Пользовательские функции потерь и метрики качества
- 7. Производительность
- Литература
1. Установка пакета xgboost
Библиотека xgboost написана на C++ и может использоваться как автономно (при помощи интерфейса командой строки), так и при помощи библиотек-интерфейсов для R, Python, Julia и Scala. В этом сообщении рассматривается R-версия, работа с другими языками выглядит аналогично.Устанавливать можно с CRAN, но лучше использовать наиболее свежую версию версию с GitHub (актуально для всех ОС):
install.packages("drat", repos = "https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("xgboost", repos = "http://dmlc.ml/drat/", type = "source")2. Градиентный бустинг: принцип метода и особенности реализации в xgboost
Идея градиентного бустинга состоит в построении ансамбля последовательно уточняющих друг друга элементарных моделей. n-ная элементарная модель обучается на “ошибках” ансамбля из n-1 моделей, ответы моделей взвешенно суммируются. “Ошибки” здесь в кавычках, поскольку на самом деле каждая последующая модель приближает антиградиент функции потерь, который не обязательно равен разности фактических и предсказанных значений (т.е. ошибке в буквальном смысле). Различные функции потерь имеют разные производные, но для среднеквадратичной функции потерь, заданной как \(\frac{1}{2}\left [ y_{i} - f\left ( y_{i} \right ) \right ]^{2}\), антиградиент (производная с обратным знаком) представляет собой именно разность между фактическими и предсказанными значениями: \(y_{i} - f\left ( y_{i} \right )\). Это наиболее интуитивный вариант, далее в качестве более сложного примера будет рассмотрена логистическая функция потерь.Элементарная модель здесь названа элементарной вовсе не потому, что она обязана быть очень простой, такой как неглубокое решающее дерево. Под элементарностью подразумевается возможность присутствия модели в качестве составляющей части модели более высокого порядка, в данном случае - модели градиентного бустинга. Бустить можно практически какие угодно модели - общие линейные, обобщенные линейные, деревья решений, K-ближайших соседей и другие, например, на StackOverflow в одном из ответов упоминался бустинг нейросетей с 5 скрытыми слоями. То есть как такового ограничения на сложность модели нет, а есть лишь общая для машинного обучения дилемма смещения-дисперсии (bias-variance tradeoff): модель должна быть достаточно гибкой, чтобы восстанавливать искомую зависимость, но при этом по возможности не должна переобучаться. Элементарная модель в градиентном бустинге известна также как weak learner (непереводимый и довольно абстрактный термин).
К особенностям реализации бустингового алгоритма в xgboost можно отнести использование помимо первой еще и второй производной от функции потерь, наличие встроенной регуляризации, а также возможность задавать пользовательские функции потерь и метрики качества. Первая сугубо техническая особенность повышает эффективность алгоритма. Встроенная регуляризация помогает бороться с переобучением: на очередной итерации решающее дерево не будет строиться до максимальной глубины, если это слишком незначительно улучшает качество модели ценой значительного её усложнения. Наконец, в качестве функции потерь пользователь может задать любую функцию, у которой есть непрерывная первая и вторая производная. Отметим, что в более раннем пакете gbm представлен широкий ассортимент функций потерь в зависимости от распределения целевой переменной:
Currently available options are "gaussian" (squared error), "laplace" (absolute loss), "tdist" (t-distribution loss), "bernoulli" (logistic regression for 0-1 outcomes), "huberized" (huberized hinge loss for 0-1 outcomes), "multinomial" (classification when there are more than 2 classes), "adaboost" (the AdaBoost exponential loss for 0-1 outcomes), "poisson" (count outcomes), "coxph" (right censored observations), "quantile", or "pairwise" (ranking measure using the LambdaMart algorithm).3. Использование xgboost совместно с caret
xgboost имеет множество гиперпараметров, которые нужно как-то настраивать. К счастью, пакет caret в числе прочих моделей поддерживает и xgboost. К сожалению, перебрать по сетке сколько-нибудь значительное количество комбинаций всех гиперпараметров не представляется возможным, поэтому при использовании бустинга решающих деревьев (method = "xgbLinear" в caret::train = booster = "gblinear" в xgboost::xgb.train) подбираются только значения nrounds (количество итераций), max_depth (максимальная глубина дерева), eta (скорость обучения), gamma (минимальное уменьшения значения функции потерь), colsample_bytree (доля переменных, используемых при построении каждого дерева) и min_child_weight (для регрессии - минимальное количество наблюдений в листе дерева).В качестве примера для решения задачи регрессии возьмем набор данных о растворимости химических соединений, рассматриваемый в книге Applied Predictive Modeling. Будем иcпользовать трансформированные предикторы
solTrainXtrans.library(AppliedPredictiveModeling)
library(caret)## Loading required package: lattice## Loading required package: ggplot2library(xgboost)
data(solubility)
ls(pattern = "^solT")## [1] "solTestX" "solTestXtrans" "solTestY" "solTrainX"
## [5] "solTrainXtrans" "solTrainY"dim(solTrainX)## [1] 951 228xgbGrid <- expand.grid(nrounds = 200,
max_depth = c(3, 4, 5),
eta = seq(0.02, 0.1, by = 0.02),
gamma = 0,
colsample_bytree = c(0.6, 0.8),
min_child_weight = 1)set.seed(100)
fitControl <- trainControl(method = "cv",
number = 3)
fit1 <- train(x = solTrainXtrans,
y = solTrainY,
method = "xgbTree",
trControl = fitControl,
metric = "RMSE",
tuneGrid = xgbGrid)## Loading required package: plyrfit1## eXtreme Gradient Boosting
##
## 951 samples
## 228 predictors
##
## No pre-processing
## Resampling: Cross-Validated (3 fold)
## Summary of sample sizes: 634, 635, 633
## Resampling results across tuning parameters:
##
## eta max_depth colsample_bytree RMSE Rsquared
## 0.02 3 0.6 0.7755540 0.8662638
## 0.02 3 0.8 0.7746227 0.8665563
## 0.02 4 0.6 0.7216388 0.8827107
## 0.02 4 0.8 0.7245405 0.8808891
## 0.02 5 0.6 0.6977385 0.8893740
## 0.02 5 0.8 0.7015094 0.8877743
## 0.04 3 0.6 0.6934022 0.8870933
## 0.04 3 0.8 0.7007086 0.8845728
## 0.04 4 0.6 0.6629013 0.8963019
## 0.04 4 0.8 0.6646522 0.8956216
## 0.04 5 0.6 0.6607801 0.8966548
## 0.04 5 0.8 0.6588977 0.8972954
## 0.06 3 0.6 0.6686043 0.8939618
## 0.06 3 0.8 0.6690120 0.8939003
## 0.06 4 0.6 0.6438700 0.9013057
## 0.06 4 0.8 0.6523668 0.8987333
## 0.06 5 0.6 0.6487391 0.9003941
## 0.06 5 0.8 0.6550886 0.8978224
## 0.08 3 0.6 0.6599967 0.8960561
## 0.08 3 0.8 0.6557474 0.8977300
## 0.08 4 0.6 0.6474405 0.9003465
## 0.08 4 0.8 0.6500267 0.8994681
## 0.08 5 0.6 0.6470897 0.9004721
## 0.08 5 0.8 0.6497202 0.8997081
## 0.10 3 0.6 0.6531061 0.8982074
## 0.10 3 0.8 0.6481044 0.8995590
## 0.10 4 0.6 0.6509935 0.8987497
## 0.10 4 0.8 0.6416907 0.9019764
## 0.10 5 0.6 0.6594815 0.8963834
## 0.10 5 0.8 0.6677941 0.8932863
##
## Tuning parameter 'nrounds' was held constant at a value of 200
##
## Tuning parameter 'gamma' was held constant at a value of 0
##
## Tuning parameter 'min_child_weight' was held constant at a value of 1
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nrounds = 200, max_depth = 4,
## eta = 0.1, gamma = 0, colsample_bytree = 0.8 and min_child_weight = 1.pred <- predict(fit1, solTestXtrans)
sqrt(mean((solTestY - pred) ^ 2))## [1] 0.6114814. xgb.train: формат исходных данных и параметры модели
Основной функцией xgboost для создания моделей являетсяxgb.train. Также доступна функция xgboost с более простым интерфейсом. Данные для xgb.train должны быть в формате xgb.DMatrix, в то время как xgboost также принимает на вход обычные матрицы matrix и разреженные матрицы dgCMatrix. Создадим объекты класса xgb.DMatrix с обучающей и тестовой выборками. Необходимо задать матрицу признаков и вектор ответов, также можно задать веса наблюдений (параметр weight).xgb_train <- xgb.DMatrix(as.matrix(solTrainXtrans),
label = solTrainY)
xgb_test <- xgb.DMatrix(as.matrix(solTestXtrans),
label = solTestY)- общих параметров (главный из них - тип бустера);
- параметров, специфичных для каждого бустера;
- параметров решаемой задачи.
booster = "gbtree" - вариант по умолчанию) и его модификацию DART (booster = "dart"), а также бустинг на основе линейных моделей (booster = "gblinear"). Мы рассмотрим первые два, поскольку бустинг линейных моделей обычно менее результативен.Имеется краткое руководство по настройке параметров: параметры
max_depth, min_child_weight и gamma непосредственно ограничивают сложность модели, subsample и colsample_bytree делают её более устойчивой к шуму за счет добавления случайного выбора наблюдений и предикторов, eta влияет на скорость обучения. Авторские рекомендации по настройке параметров, не претендующие на истину в последней инстанции:max.depth- максимальная глубина дерева. Критически важный параметр: выбор слишком глубоких деревьев приводит к переобучению, а слишком маленькие деревья не позволят эффективно восстановить искомую зависимость. Значение по умолчанию равно 6, имеет смысл попробовать значения чуть больше и чуть меньше.eta- скорость обучения модели. Критически важный параметр, контролирующий, с каким весом предсказания каждой следующей модели суммируются с предсказаниями ансамбля. Значение по умолчанию (0.3) является слишком большим, обычно хорошо работают значения меньше 0.1. Слишком маленьким его сделать тяжело, уменьшениеetaкомпенсируется увеличением количества итераций (ценой большего времени обучения, разумеется).gamma- минимальное уменьшения значения функции потерь. Можно оставить значение 0 или сделать равным какому-то небольшому числу, например 0.01.subsample- доля объектов обучающей выборки, используемых на каждой итерации. Еще один критически важный параметр, позволяющий значительно улучшить модель. По умолчанию равен 1, значение 0.5 может оказаться хорошим выбором (но нужно принимать во внимание размер выборки).colsample_bytree- доля переменных, используемых на каждой итерации. Вместе с предыдущим параметром реализует преимущества алгоритма Random Forest (на каждой итерации используется только часть наблюдений и предикторов). Обычно хорошо работают значения между 0.5 и 1, типичные для Random Forest значения скорее всего окажутся слишком маленькими.min_child_weight- минимальное количество наблюдений в листе дерева. Имеет смысл проверить несколько значений больше 1 (например, 3, 5 и более).
param <- list(booster = "gbtree",
max.depth = 4,
eta = 0.07,
gamma = 0,
subsample = 0.5,
colsample_bytree = 0.8,
min_child_weight = 3,
objective = "reg:linear",
eval_metric = "rmse")set.seed(100)
fit2 <- xgb.train(data = xgb_train,
nrounds = 400,
params = param,
verbose = 1,
print_every_n = 50,
watchlist = list(train = xgb_train, test = xgb_test))## [1] train-rmse:3.576950 test-rmse:3.657048
## [51] train-rmse:0.555615 test-rmse:0.690631
## [101] train-rmse:0.414470 test-rmse:0.630226
## [151] train-rmse:0.351295 test-rmse:0.605188
## [201] train-rmse:0.310165 test-rmse:0.594962
## [251] train-rmse:0.278760 test-rmse:0.588556
## [301] train-rmse:0.253982 test-rmse:0.585576
## [351] train-rmse:0.234476 test-rmse:0.586793
## [400] train-rmse:0.217970 test-rmse:0.583388verbose = 0 должен подавлять вывод лога в консоль, однако в действительности заодно отключает и его сохранение в модельном объекте. Поэтому указываем print_every_n = 50, чтобы сократить вывод результатов. Можно посмотреть или нарисовать весь лог обучения:fit2$evaluation_log## iter train_rmse test_rmse
## 1: 1 3.576950 3.657048
## 2: 2 3.353263 3.429678
## 3: 3 3.143332 3.211181
## 4: 4 2.948522 3.015424
## 5: 5 2.764133 2.829077
## ---
## 396: 396 0.219104 0.583646
## 397: 397 0.218913 0.583890
## 398: 398 0.218404 0.583828
## 399: 399 0.218265 0.583748
## 400: 400 0.217970 0.583388plot(fit2$evaluation_log$test_rmse, main = "test_rmse")
Проверка на тестовой выборке (результаты соответствуют последней строке приведенного выше лога):
pred <- predict(fit2, xgb_test)
sqrt(mean((solTestY - pred) ^ 2))## [1] 0.583388importance_frame <- xgb.importance(colnames(xgb_train), model = fit2)
xgb.plot.importance(importance_frame[1:10, ]) 
Для сравнения построим модель с бустером DART и незначительно измененными параметрами. Результат слегка улучшился:
param <- list(booster = "dart",
max.depth = 4,
eta = 0.07,
gamma = 0.01,
subsample = 0.6,
colsample_bytree = 0.7,
min_child_weight = 4,
objective = "reg:linear",
eval_metric = "rmse")
set.seed(100)
fit3 <- xgb.train(data = xgb_train,
nrounds = 400,
params = param,
verbose = 1,
print_every_n = 50,
watchlist = list(train = xgb_train, test = xgb_test))## [1] train-rmse:3.576602 test-rmse:3.659974
## [51] train-rmse:0.547767 test-rmse:0.673526
## [101] train-rmse:0.409768 test-rmse:0.612489
## [151] train-rmse:0.352392 test-rmse:0.592270
## [201] train-rmse:0.308753 test-rmse:0.581850
## [251] train-rmse:0.276780 test-rmse:0.569760
## [301] train-rmse:0.251984 test-rmse:0.567371
## [351] train-rmse:0.234161 test-rmse:0.565282
## [400] train-rmse:0.215848 test-rmse:0.564427xgb.cv для перекрестной проверки. Она может быть полезной, но пакет caret предлагает гораздо больше возможностей для оценки качества моделей.5. xgboost и Random Forest
Благодаря экспериментальному параметруnum_parallel_tree можно задать количество одновременно создаваемых деревьев и представить Random Forest как частный случай бустинговой модели с одной итерацией. А если использовать больше одной итерации, то получится бустинг “случайных лесов”, когда каждый “случайный лес” выступает в качестве элементарной модели.param <- list(booster = "gbtree",
max.depth = 25,
eta = 1,
gamma = 0,
subsample = 1,
colsample_bytree = 0.4,
min_child_weight = 1,
num_parallel_tree = 200,
objective = "reg:linear", eval_metric = "rmse")
set.seed(100)
fit2 <- xgb.train(data = xgb_train,
nrounds = 5,
params = param,
verbose = 1,
watchlist = list(train = xgb_train, test = xgb_test))## [1] train-rmse:0.625530 test-rmse:0.783342
## [2] train-rmse:0.221222 test-rmse:0.651380
## [3] train-rmse:0.114044 test-rmse:0.640410
## [4] train-rmse:0.076963 test-rmse:0.639983
## [5] train-rmse:0.063586 test-rmse:0.6397286. Пользовательские функции потерь и метрики качества
Рассмотрим логистическую функцию потерь из примера. Чтобы задать целевую функцию в xgboost, нужно найти первую и вторую производную функции потерь. Они имеют видgrad <- preds - labels # первая производная
hess <- preds * (1 - preds) # вторая производнаяЛогистическая функция потерь имеет вид
\[log(1+e^{-yP})\]
где \(P\) является логарифмом отношения шансов, а \(y\) является бинарной переменной отклика (0 или 1).
С другой стороны, логистическая функция потерь равна логарифмической функции правдоподобия, взятой с противоположным знаком. Работать будем непосредственно с этой функцией, знак поменяем в самом конце.
Логарифмическая функция правдоподобия:
\[L=y_{i}\cdot log(p_{i})+(1-y_{i})\cdot log(1-p_{i})\]
\(p_{i}\) - значение логистической функции (вероятность): \(p_{i}=\frac{1}{1+e^{-\hat{y}_{i}}}\), где \(\hat{y}_{i}\) - предсказанные значения до логистической трансформации (т.е. логарифм отношения шансов):
\[L=y_{i}\cdot log\left(\frac{1}{1+e^{-\hat{y}_{i}}}\right)+(1-y_{i})\cdot log\left(\frac{e^{-\hat{y}_{i}}}{1+e^{-\hat{y}_{i}}}\right)\]
Первая производная, полученная с помощью Wolfram Alpha:
\[{L}'=\frac{y_{i}-(1-y_{i})\cdot e^{\hat{y}_{i}}}{1+e^{\hat{y}_{i}}}\]
После умножения на \(\frac{e^{-\hat{y}_{i}}}{e^{-\hat{y}_{i}}}\):
\[{L}'=\frac{y_{i}\cdot e^{-\hat{y}_{i}}+y_{i}-1}{1+e^{-\hat{y}_{i}}}= \frac{y_{i}\cdot (1+e^{-\hat{y}_{i}})}{1+e^{-\hat{y}_{i}}}-\frac{1}{1+e^{-\hat{y}_{i}}}=y_{i}-p_{i}\]
Меняем знак и получаем выражение для градиента логистической функции потерь:
\[p_{i}-y_{i}\]
Также мы можем записать первую производную (уже непосредственно функции потерь) как
\[{L}'=\frac{1}{1+e^{-\hat{y}_{i}}}-y_{i}\]
и взять производную от производной (снова с помощью Wolfram Alpha):
\[{L}''=\frac{e^{\hat{y}_{i}}}{\left( 1 + e^{\hat{y}_{i}} \right)^{2}} = \frac{e^{\hat{y}_{i}}}{1 + e^{\hat{y}_{i}}}\cdot\frac{1}{1 + e^{\hat{y}_{i}}}\]
Умножив каждый из двух множителей на \(\frac{e^{-\hat{y}_{i}}}{e^{-\hat{y}_{i}}}\), получим
\[{L}''= \frac{1}{1 + e^{-\hat{y}_{i}}}\cdot\frac{e^{-\hat{y}_{i}}}{1 + e^{-\hat{y}_{i}}}=\frac{1}{1 + e^{-\hat{y}_{i}}}\cdot\frac{1+e^{-\hat{y}_{i}}-1}{1 + e^{-\hat{y}_{i}}}=p(1-p)\]
Мы получили вторую производную \(p(1-p)\). Это в точности те формулы, которые используются в рассматриваемом примере:
logregobj <- function(preds, dtrain) {
labels <- getinfo(dtrain, "label")
# Получаем вероятности из логарифма отношения шансов:
preds <- 1/(1 + exp(-preds))
grad <- preds - labels
hess <- preds * (1 - preds)
return(list(grad = grad, hess = hess))
}
# Соответствующая метрика качества:
evalerror <- function(preds, dtrain) {
labels <- getinfo(dtrain, "label")
err <- as.numeric(sum(labels != (preds > 0)))/length(labels)
return(list(metric = "custom_error", value = err))
}data(agaricus.train, package='xgboost')
data(agaricus.test, package='xgboost')
dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
watchlist <- list(eval = dtest, train = dtrain)
num_round <- 3
param <- list(max_depth = 2,
eta = 1,
nthread = 2,
silent = 1,
objective = logregobj,
eval_metric = evalerror)
xgb.train(param, dtrain, num_round, watchlist)## [1] eval-custom_error:0.042831 train-custom_error:0.046522
## [2] eval-custom_error:0.021726 train-custom_error:0.022263
## [3] eval-custom_error:0.006207 train-custom_error:0.007063## ##### xgb.Booster
## raw: 1.5 Kb
## call:
## xgb.train(params = param, data = dtrain, nrounds = num_round,
## watchlist = watchlist)
## params (as set within xgb.train):
## max_depth = "2", eta = "1", nthread = "2", silent = "1", silent = "1"
## xgb.attributes:
## niter
## callbacks:
## cb.print.evaluation(period = print_every_n)
## cb.evaluation.log()
## niter: 3
## evaluation_log:
## iter eval_custom_error train_custom_error
## 1 0.042830540 0.046522340
## 2 0.021725636 0.022263166
## 3 0.006207325 0.007062797# Воспроизводит результат:
param <- list(max_depth = 2,
eta = 1,
nthread = 2,
silent = 1,
objective = "binary:logistic",
eval_metric = "error")
xgb.train(param, dtrain, num_round, watchlist)## [1] eval-error:0.042831 train-error:0.046522
## [2] eval-error:0.021726 train-error:0.022263
## [3] eval-error:0.006207 train-error:0.007063## ##### xgb.Booster
## raw: 1.5 Kb
## call:
## xgb.train(params = param, data = dtrain, nrounds = num_round,
## watchlist = watchlist)
## params (as set within xgb.train):
## max_depth = "2", eta = "1", nthread = "2", silent = "1", objective = "binary:logistic", eval_metric = "error", silent = "1"
## xgb.attributes:
## niter
## callbacks:
## cb.print.evaluation(period = print_every_n)
## cb.evaluation.log()
## niter: 3
## evaluation_log:
## iter eval_error train_error
## 1 0.042831 0.046522
## 2 0.021726 0.022263
## 3 0.006207 0.007063watchlist <- list(train = dtrain)
num_round <- 3
param <- list(max_depth = 2,
eta = 1,
nthread = 2,
silent = 1,
objective = logregobj,
eval_metric = "logloss",
eval_metric = "error")
xgb.train(param, dtrain, num_round, watchlist)## [1] train-logloss:1.713946 train-error:0.046522
## [2] train-logloss:0.820205 train-error:0.022263
## [3] train-logloss:0.270822 train-error:0.015200## ##### xgb.Booster
## raw: 1.5 Kb
## call:
## xgb.train(params = param, data = dtrain, nrounds = num_round,
## watchlist = watchlist)
## params (as set within xgb.train):
## max_depth = "2", eta = "1", nthread = "2", silent = "1", eval_metric = "logloss", eval_metric = "error", silent = "1"
## xgb.attributes:
## niter
## callbacks:
## cb.print.evaluation(period = print_every_n)
## cb.evaluation.log()
## niter: 3
## evaluation_log:
## iter train_logloss train_error
## 1 1.713946 0.046522
## 2 0.820205 0.022263
## 3 0.270822 0.0152007. Производительность
Высокая производительность xgboost обеспечивается в том числе за счет многопоточности:Parallelization is automatically enabled if OpenMP is present. Number of threads can also be manually specified via nthread parameter.Для работы с большими выборками пригодится возможность загрузки данных с диска. А для работы с очень большими данными можно запускать xgboost поверх систем распределенных вычислений.
Андрей, здравствуйте!
ОтветитьУдалитьспасибо за интересную и актуальную статью.
в 4 пункте, у вас, вероятно, опечатка: subsample по умолчанию равен 1 (а не 0).
Спасибо, исправил.
УдалитьПривет, можно вопрос по R в лс?
ОтветитьУдалитьDon't ask to ask и все такое. Спрашивайте сразу, без предисловий:)
ОтветитьУдалитьАндрей, а почему caret не позволяет вставить регуляризацию: альфу и лямбду, не интересовались?
ОтветитьУдалитьКомбинаторный взрыв: если все сразу тюнить, количество комбинаций гиперпараметров будет слишком баснословным. Поэтому выбраны самые важные (их и так очень много по сравнению с тем же "случайным лесом"). Но можно написать свою модель для caret, которая будет перебирать и эти параметры. Или pipelearner использовать для тех же целей, что даже проще.
Удалить