Bayesian Optimization for Turning XGBoost Hyperparameter

R Machine Learning Series

About Xgboost

Xgboost là một thuật toán Machine Learning được phát triển và hình thành từ một dự án của University of Washington bởi Tianqi Chen Carlos Guestrin. Ra đời khá muộn nhưng thuật toán này thường đạt thứ hạng cao trong các cuộc thi về phân tích dữ liệu trên Kaggle. Hiện tại thuật toán này đang có khoảng gần 400 contributors. Lịch sử ra đời của thuật toán này có thể tham khảo từ chính tác giả Tianqi Chen tại đây.

Là một open source có thể chạy trên cả hai hệ hiều hành Windown và Linux cho các ngôn ngữ phổ biến như C++, Java, Python, R, Julia, và Perl và có thể áp dụng cho một loạt các bài toán từ Regression, Classification đến các bài toán dự báo được định nghĩa riêng đặc trưng cho từng dự án (user-defined prediction) và cũng có thể thực thi trên các nền tảng cloud computing như AWS, Azure. Chi tiết về thuật toán này có thể tham khảo ở đây.

Conduct Xgboost in R

Thuật toán này có thể được thực hiện bằng một số thư viện như caret - một tool kit cho Machine Learning tương tự như scikit learn của Python hoặc trực tiếp sử dụng thư viện xgboost linh hoạt hơn và cho phép tinh chỉnh đầy đủ các tham số của thuật toán này và đây cũng là thư viện được sử dụng trong post này theo trình tự sau:

1. Chuẩn bị dữ liệu cho Xgboost.
2. Thực hiện huấn luyện Xgboost mặc định.
3. Sử dụng Bayesian Optimization tinh chỉnh và tìm kiếm tham số tối ưu.

Dữ liệu sử dụng là hmeq.csv ở textbook Wiley and SAS Business: Credit Risk Analytics: Measurement Techniques, Applications, and Examples in SAS.

Data Preparation

Missing Data and Imputation

Các thuật toán Machine Learning đòi hỏi cả features lẫn target phải là đầy đủ. Với dữ liệu trống (missing) thì có thể có nhiều phương pháp thay thế dữ liệu trống (imputation) cho từng loại feature khác nhau. Trước hết chúng ta đánh giá mức độ thiếu của hai nhóm dữ liệu:

#=================================
#       Data Pre-processing
#=================================

# Load some packages for data manipulation: 
library(tidyverse)

# Clear workspace: 
rm(list = ls())

# Import data: 
hmeq <- read_csv("http://www.creditriskanalytics.net/uploads/1/9/5/1/19511601/hmeq.csv")

# Glimpse original data: 
library(knitr)

hmeq %>% 
  head() %>% 
  kable()

BAD	LOAN	MORTDUE	VALUE	REASON	JOB	YOJ	DEROG	DELINQ	CLAGE	NINQ	CLNO	DEBTINC
1	1100	25860	39025	HomeImp	Other	10.5	0	0	94.36667	1	9	NA
1	1300	70053	68400	HomeImp	Other	7.0	0	2	121.83333	0	14	NA
1	1500	13500	16700	HomeImp	Other	4.0	0	0	149.46667	1	10	NA
1	1500	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
0	1700	97800	112000	HomeImp	Office	3.0	0	0	93.33333	0	14	NA
1	1700	30548	40320	HomeImp	Other	9.0	0	0	101.46600	1	8	37.11361

# Check missing data for categorical features: 
sapply(hmeq %>% select_if(is.character), function(x) {table(x, useNA = "ifany")})

## $REASON
## x
## DebtCon HomeImp    <NA> 
##    3928    1780     252 
## 
## $JOB
## x
##     Mgr  Office   Other ProfExe   Sales    Self    <NA> 
##     767     948    2388    1276     109     193     279

# Check missing rate for numeric features: 
sapply(hmeq %>% select_if(is.numeric), function(x) {sum(is.na(x)) / length(x)})

##        BAD       LOAN    MORTDUE      VALUE        YOJ      DEROG     DELINQ 
## 0.00000000 0.00000000 0.08691275 0.01879195 0.08640940 0.11879195 0.09731544 
##      CLAGE       NINQ       CLNO    DEBTINC 
## 0.05167785 0.08557047 0.03724832 0.21258389

Biến REASON có 252 điểm dữ liệu missing. Còn với biến định lượng thì, ví dụ, biến DEBTINC thì tỉ lệ missing là 21.26%. Với biến định lượng (numeric) thì đơn giản nhất là thay thế bằng mean hoặc median. Với biến định tính (categorical) thì có thể thay thế các quan sát missing bằng các nhãn còn lại sao cho tỉ lệ (hay distribution) của cá nhãn thuộc biến được thay thế không thay đổi hoặc xấp xỉ tỉ lệ của các nhãn quan sát được trước khi thay thế. Phương pháp thay thế dữ liệu trống này là đơn giản nhất và sẽ không được sử dụng trong tình huống của bộ dữ liệu hmeq. Thay vì đó chúng ta thay thế dữ liệu trống bằng thuật toán Random Forest của thư viện missRanger như sau:

#---------------------------------------------------
#  Stage 1: Imputate missing data by Random Forest
#---------------------------------------------------

# Impute missing data: 
library(missRanger)

missRanger(hmeq, seed = 29) -> hmeq_imputed

## 
## Missing value imputation by random forests
## 
##   Variables to impute:       MORTDUE, VALUE, REASON, JOB, YOJ, DEROG, DELINQ, CLAGE, NINQ, CLNO, DEBTINC
##   Variables used to impute:  BAD, LOAN, MORTDUE, VALUE, REASON, JOB, YOJ, DEROG, DELINQ, CLAGE, NINQ, CLNO, DEBTINC
## iter 1:  ...........
## iter 2:  ...........
## iter 3:  ...........

Chúng ta có thể kiểm tra để xác nhận rằng các biến có missing data point đã được lấp đầy. Ví dụ với nhóm biến định tính:

# Check data set after imputing: 
sapply(hmeq_imputed %>% select_if(is.character), function(x) {table(x, useNA = "ifany")})

## $REASON
## x
## DebtCon HomeImp 
##    4101    1859 
## 
## $JOB
## x
##     Mgr  Office   Other ProfExe   Sales    Self 
##     784     990    2537    1337     119     193

Feature Engineering for Categorical Variables

Các thuận toán ML hầu hết (nếu không muốn nói là tất cả) đều đòi hỏi rằng các biến số sử dụng để huấn luyện mô hình phải là numeric nếu sử dụng Python. Với R thì yêu cầu này có thể không đòi hỏi. Tuy nhiên thuật toán Xgboost vẫn yêu cầu các biến phải là numeric cho cả R lẫn Python. Do vậy với biến categorical chúng ta có thể chuyển về numeric bằng một kĩ thuật gọi là one-hot encoding. Chúng ta có thể thực hiện kĩ thuật Feature Engineering này với hàm dummyVars() của thư viện caret như sau:

#---------------------------------------------
#  Stage 2: Conduct one-hot encoding process
#---------------------------------------------

library(caret)

# Conduct one-hot encoding for categorical features: 
dummies <- dummyVars(~ REASON + JOB, data = hmeq_imputed)
predict(dummies, hmeq_imputed) %>% as.data.frame() -> features_oneHot

# Final data for modelling: 
hmeq_imputed %>% 
  select_if(is.numeric) %>% 
  bind_cols(features_oneHot) -> df_final

So sánh với dữ liệu ban đầu:

# Compare with original data: 
df_final %>% 
  head() %>% 
  select(BAD, LOAN, contains("JOB")) %>% 
  kable()

BAD	LOAN	JOBMgr	JOBOffice	JOBOther	JOBProfExe	JOBSales	JOBSelf
1	1100	0	0	1	0	0	0
1	1300	0	0	1	0	0	0
1	1500	0	0	1	0	0	0
1	1500	0	0	1	0	0	0
0	1700	0	1	0	0	0	0
1	1700	0	0	1	0	0	0

Đến đây chúng ta có thể sử dụng dữ liệu đã xử lí để huấn luyện Xgboost.

Train Xgboost

Thuật toán Xgboost đòi hỏi rằng dữ liệu không phải ở dạng data frame như thường thấy mà nên là DMatrix form bằng hàm xgb.DMatrix(). Trước hết phân chia dữ liệu ban đầu thành tỉ lệ 70 - 30 rồi chuyển về dạng DMatrix cho train data như sau:

#============================
#   Modelling with XGBoost
#============================

# Prepare data: 

set.seed(1)
id <- createDataPartition(df_final$BAD, p = 0.7, list = FALSE)
train <- df_final[id, ]
test <- df_final[-id, ]


# Convert features to DMatrix form: 

X_train <- train %>% 
  select(-BAD) %>% 
  as.matrix()

Y_train <- train %>% 
  pull(BAD)

X_test <- test %>% 
  select(-BAD) %>% 
  as.matrix()

Y_test <- test %>% 
  pull(BAD)

#------------------------------------------
#   Train XGBoost with default parameters
#------------------------------------------
library(xgboost)

# Convert to DMatrix form for train data: 
dtrain <- xgb.DMatrix(data = X_train, label = Y_train)

Với train data đã ở dạng phù hợp với đòi hỏi của Xgboost chúng ta có thể huấn luyện thuật toán này mà không tinh chỉnh bất cứ tham số nào (default parameters) bằng hàm xgboost:

# Train a default XGBoost: 
xgb1 <- xgboost(data = dtrain, 
                objective = "binary:logistic", 
                verbose = 0, 
                nround = 30)

Để đánh giá chất lượng của mô hình trước hết chúng ta viết một hàm tính toán confusion matrix tương ứng với ngưỡng được chọn:

# Function for calculates confusion matrix: 

my_cm <- function(xgb_obj, test_data, nguong) {
  
  prediction <- predict(xgb_obj, test_data)
  
  pred <- case_when(prediction >= nguong ~ "Bad", TRUE ~ "Good") %>% as.factor()
  
  thuc_te <- case_when(Y_test == 1 ~ "Bad", Y_test == 0 ~ "Good") %>% as.factor()
  
  confusionMatrix(pred, thuc_te, positive = "Bad")
}

Khả năng dự báo của mô hình khi ngưỡng là 0.5 (đây là ngưỡng mặc định được lựa chọn cho phân loại cả ở R lẫn Python):

# Confusion matrix with cutoff = 0.5: 
my_cm(xgb1, X_test, 0.5)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  Bad Good
##       Bad   249   17
##       Good  103 1419
##                                          
##                Accuracy : 0.9329         
##                  95% CI : (0.9203, 0.944)
##     No Information Rate : 0.8031         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.7662         
##                                          
##  Mcnemar's Test P-Value : 8.533e-15      
##                                          
##             Sensitivity : 0.7074         
##             Specificity : 0.9882         
##          Pos Pred Value : 0.9361         
##          Neg Pred Value : 0.9323         
##              Prevalence : 0.1969         
##          Detection Rate : 0.1393         
##    Detection Prevalence : 0.1488         
##       Balanced Accuracy : 0.8478         
##                                          
##        'Positive' Class : Bad            
## 

Cũng có thể đánh giá vai trò - mức độ quan trọng của biến số trong việc phân loại (feature importance) bằng hàm xgb.importance() rồi hình ảnh hóa:

# Extract feature importances in the XGBoost trained: 

importance <- xgb.importance(feature_names = colnames(X_train), model = xgb1)

# Plot feature importance by Gain: 


my_colors <- c("#014d64", "#01a2d9")

importance %>% 
  as.data.frame() %>% 
  mutate(Feature = factor(Feature, Feature)) %>% 
  ggplot(aes(Feature, Gain)) + 
  geom_col(fill = my_colors[2]) + 
  coord_flip() + 
  theme_minimal() + 
  labs(x = NULL, y = NULL, 
       title = "Figure 1: Feature Importance by Gain for XGBoost", 
       subtitle = "Data Source: http://www.creditriskanalytics.net")

Hình dáng của đường ROC và giá trị AUC cũng là một tiêu chí được quan tâm hàng đầu của mô hình phân loại. Chúng ta có thể đánh chúng trên bộ dữ liệu test như sau:

# Use Xgboost for predicting: 
pd <- predict(xgb1, X_test)

# Calculate AUC and plot ROC curve: 
library(pROC)

my_auc <- roc(Y_test, pd)
sen_spec_df <- data_frame(TPR = my_auc$sensitivities, 
                          FPR = 1 - my_auc$specificities, 
                          total = my_auc$sensitivities + my_auc$specificities, 
                          cutoff = my_auc$thresholds)


sen_spec_df %>% 
  ggplot(aes(x = FPR, ymin = 0, ymax = TPR))+
  geom_polygon(aes(y = TPR), fill = "red", alpha = 0.3) +
  geom_path(aes(y = TPR), col = "firebrick", size = 1.2) +
  geom_abline(intercept = 0, slope = 1, color = "gray37", size = 1, linetype = "dashed") + 
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  theme_bw() +
  coord_equal() +
  labs(x = "FPR (1 - Specificity)", 
       y = "TPR (Sensitivity)", 
       title = "Figure 2: Model Performance for Xgboost based on Test Data", 
       subtitle = paste0("AUC Value: ", my_auc$auc %>% round(4)))

Bayesian Optimization

Không tinh chỉnh bất cứ tham số nào cho Xgboost chúng ta cũng có được một mô hình phân loại với AUC = 0.9536. Khả năng phân loại của mô hình có thể sẽ tốt hơn nữa nếu chúng ta tinh chỉnh - tìm kiếm tham số tối ưu cho Xgboost. Hai cách thức tinh chỉnh phổ biến là:

1. Full Grid Search. Phương pháp này liệt kê tất cả những sự kết hợp khác nhau của tham số rồi thực hiện tinh chỉnh tương ứng. Cách thức này tốn rất nhiều thời gian và không hiệu quả - ít nhất là với trường hợp dữ liệu là lớn, nhiều biến số.

2. Random Search. Phương pháp này là một giải pháp cải tiến hơn so với Full Grid Search và sẽ mất ít thời gian hơn đi kèm với cái giá phải trả là phương án tốt nhất (hay bộ tham số tốt nhất) có thể bị bỏ qua. Về cơ chiến thuật này vẫn phải liệt kê ra tất cả các sự kết hợp của tham số và, ví dụ, quá trình tìm kiếm sẽ dừng lại nếu sau một số bước nào đó, model performance tính theo một tiêu chí chọn trước (như ROC-AUC chẳng hạn) không được cải thiện.

3. Bayesian Optimization. Phương pháp tinh chỉnh - tìm kiếm tham số tối ưu này hiệu quả hơn cả vì việc lụa chọn tham số cho bước kết tiếp sẽ được dựa vào thông tin về chất lượng của mô hình trước đó. Chi tiết về phương pháp này có thể tham khảo từ nhiều nguồn, ví dụ như ở đây.

Với R chúng ta có thể sử dụng thư viện rBayesianOptimization để thực hiện tinh chỉnh tham số cho Xgboost theo phương pháp này như sau:

1. Thiết lập domain space - là không gian các tham số sẽ được tinh chỉnh. Bước này khá tương tự với thiết lập Full Grid Search.
2. Viết hàm mục tiêu trả về giá trị của một tiêu chuẩn đánh giá mô hình mà chúng ta muốn tối ưu.
3. Thực hiện tinh chỉnh tham số tối ưu bằng hàm BayesianOptimization().

Các bước chuẩn bị để tinh chỉnh - tìm kiếm tham số tối ưu cũng sẽ tương tự như vậy nếu thực hiện trong bất kì một ngôn ngữ nào khác. Nếu sử dụng Python thì có thể tham khảo ở đây.

Xgboost là mô hình có rất nhiều tham số để tinh chỉnh và có thể mất nhiều thời gian. Dưới đây là R codes thực hiện tinh chỉnh cho một số tham số được lựa chọn là max.depth, min_child_weight và subsample với các tham số còn lại thì sẽ chọn hoặc để mặc định:

# Objective function: 

my_fun <- function(max.depth, min_child_weight, subsample) {
  
  xgb <- xgboost(data = dtrain, 
                 max_depth = max.depth,
                 min_child_weight = min_child_weight,
                 subsample = subsample, 
                 objective = "binary:logistic", 
                 verbose = 0, 
                 nround = 30)
  
  pd <- predict(xgb, X_test)
  my_auc <- roc(Y_test, pd)$auc %>% as.numeric()
  
  list(Score = my_auc, Pred = NULL)
  
  
  
}


# Domain space: 

my_hypers <- list(max.depth = c(2L, 6L),
                  min_child_weight = c(1L, 10L),
                  subsample = c(0.5, 0.8))



# Search optimal hyperparameter by Bayesian Optimization: 
library(rBayesianOptimization)

system.time(OPT_Res <- BayesianOptimization(my_fun, 
                                bounds = my_hypers, 
                                init_points = 10, 
                                n_iter = 20,
                                acq = "ucb", 
                                kappa = 2.576, 
                                eps = 0.0,
                                verbose = FALSE))

## 
##  Best Parameters Found: 
## Round = 23   max.depth = 6.0000  min_child_weight = 1.0000   subsample = 0.8000  Value = 0.9549

##    user  system elapsed 
## 244.222   0.104 233.110

Mất chừng khoảng 200s trên con máy Dell Workstation T5610 chạy hai CPU Intel Xeon E5-2680 V2 để tinh chỉnh và tìm kiếm tham số tối ưu. Các tham số tối ưu và AUC tương ứng trên test data:

# Best parameters: 
OPT_Res$Best_Par

##        max.depth min_child_weight        subsample 
##              6.0              1.0              0.8

# Best AUC on train data: 
OPT_Res$Best_Value

## [1] 0.9548976

AUC trên train data tương ứng với các tham số tối ưu tìm được là 0.9575 cao hơn so với AUC = 0.9536 của Xgboost mặc định. Chúng ta vẫn còn có thể nâng cao chất lượng phân loại của mô hình bằng: (1) tỉnh chỉnh thêm một số tham số, (2) điều chỉnh lại domain space, hoặc (3) sự kết hợp nào đó của (1) và (2).

Brief Summary

Post này trình bày chi tiết các bước từ chuẩn bị dữ liệu (từ xử lí missing data + one-hot encoding cho categorical feature) cho đến bước tinh chỉnh tham số tối ưu theo Bayesian Optimization. Lưu ý rằng việc tinh chỉnh - tìm kiếm tham số tối ưu không hẳn chỉ là tìm tham số để tối đa hóa AUC như chúng ta vừa thực hiện mà có thể là các mục tiêu khác. Ví dụ: với một tổ chức hoạt động vì lợi nhuận như Ngân Hàng thì lợi nhuận là mục tiêu hàng đầu. Hơạc cũng chính Ngân Hàng đó nhưng trong bối cảnh bất ổn về kinh tế thì tiêu tham số tối ưu có thể phải hướng vào việc phân loại chính xác cao nhất các khách hàng xấu khi nộp hồ sơ xin vay hay cấp tín dụng (tức là Recall). Tinh chỉnh theo tham số hướng đến tối ưu các mục tiêu kinh tế này sẽ được trình bày trong một bài viết khác.

Về mặt kĩ thuật, thư viện rBayesianOptimization không hỗ trợ tính toán song song nên thời gian tinh chỉnh có thể lâu. Khắc phục nhược điểm này chúng ta có thể sử dụng thư viện ParBayesianOptimization để tăng tốc độ tìm kiếm tham số tối ưu.

References

1. https://www.r-bloggers.com/bayesian-optimization-of-machine-learning-models/
2. https://bearloga.github.io/bayesopt-tutorial-r/
3. https://xgboost.readthedocs.io/en/latest/index.html
4. https://papers.nips.cc/paper/4522-practical-bayesian-optimization-of-machine-learning-algorithms.pdf
5. Snoek, J., Larochelle, H., & Adams, R. P. (2012). Practical bayesian optimization of machine learning algorithms. In Advances in neural information processing systems (pp. 2951-2959).
6. Bergstra, J., Komer, B., Eliasmith, C., Yamins, D., & Cox, D. D. (2015). Hyperopt: a python library for model selection and hyperparameter optimization. Computational Science & Discovery, 8(1), 014008.