GBM (Gradient Boosting Machines)에 대한 자세한 설명 (1): Regression

Posted on 2020-09-05 Edited on 2023-05-20 In machine-learning

이번 포스팅은 나무 모형 시리즈의 세 번째 글입니다. 이전 글은 AdaBoost에 대한 자세한 설명과 배깅 (Bagging)과 부스팅 (Boosting)의 원리에서 확인하실 수 있습니다.

이번 포스팅은 나무 모형 시리즈의 세 번째 글입니다. 이전 글은 AdaBoost에 대한 자세한 설명과 배깅 (Bagging)과 부스팅 (Boosting)의 원리에서 확인하실 수 있습니다.
GBM은 LightGBM, CatBoost, XGBoost가 기반하고 있는 알고리즘이기 때문에 해당 원리를 아는 것이 중요합니다.
이 포스팅은 GBM 중 Regression에 초점을 두어 설명했습니다!

Introduction

Gradient Boosting Machines (GBM)은 Boosting의 개념을 Gradient descent라는 최적화 방법으로 이해하는 방법입니다.

부스팅은 이전 포스팅에서도 설명했지만 additive하고 sequential하게 여러 나무를 학습해 결과를 종합하는 방식입니다. (~~"additive하고 sequential하게"라는 말이 "fun하고 cool하게"처럼 느껴질 수도…~~)

"additive하다"는 말은 다음과 같이 여러 나무들의 예측결과를 학습한다는 것을 의미합니다.

F(x_i) = \sum_{t=1}^T f_t(x_i)

여기서 $f_t(x_i)$ 는 $x_i$ 를 입력값으로 넣었을 때 $t$ 번째 약한 학습기가 예측한 결과를 의미하고, 결국 $F(x_i)$ 는 $T$ 개의 $f_t(x_i)$ 를 합친 $y_i$ 의 예측값을 의미합니다.

그리고 "sequential하다"는 말은 이 $T$ 개의 나무를 동시에 학습하는 것이 아니라 한 나무가 학습된 결과를 바탕으로 더 나은 결과를 위해 다시 두 번째 나무를 학습하고, 이 두 나무들보다 나은 결과를 위해 세 번째 나무를 학습하는 등 순서대로 학습한다는 의미를 내포합니다.

AdaBoost vs. GBM

그럼, GBM은 어떻게 약한 학습기로 "더 나은 결과"를 만들까요?

이를 쉽게 이해하기 위해서 AdaBoost와 비교하겠습니다.

GBM은 AdaBoost와 마찬가지로 Boosting 계열의 알고리즘이기 때문에 순차적으로 약한 학습기 (weak learner)를 만들어 가며 이전 학습기의 잔차 (resuidual)를 보완하는 방식입니다.

그러나, 이 둘은 잔차를 보완하는 방법이 다릅니다.

AdaBoost는 이전에 잘못 분류한 데이터에 가중치를 더 많이 주고, 잘 분류한 학습기에 더 많은 가중치를 주어 해결합니다. AdaBoost의 최종 모형인

F(x) = \text{sign} (\sum_{t=1}^{T} \alpha_t f_t(x))

에서 알 수 있듯이 $\alpha_t$ 로 학습기마다 가중치를 다르게 주고, 각 데이터에도 잘못 분류되면 가중치가 커지도록 adaptive하게 잔차를 보완합니다.

이에 비해 GBM은 약한 학습기를 잔차 자체에 적합하고, 이전의 예측값에 예측한 잔차를 더해 주어 예측값을 새로이 업데이트하는 방식입니다.

GBM을 간략하게 설명하면 다음과 같습니다.

$y_i, \ i=1,\ldots,n$ 에 대한 초기 예측값으로 $y_i$ 값들의 평균인 $\overline{y}$ 를 택합니다.
$t=1,\ldots, T$ $t = 1, \dots, T$ 만큼 반복
1. 각 관측치마다 잔차 $y_i - F_t(x)$ 를 구하고, 이 잔차에 첫 번째 약한 학습기를 적합시킵니다.
2. 적합된 결과에 따라 잔차의 예측값 (= 마지막 노드 (leaf)에 속한 잔차들의 평균)을 구하고
3. 잔차 예측값 * 학습률을 이전 예측값과 더해 새롭게 예측값을 업데이트합니다.

이 설명만 보면 gradient descent가 “왜 혹은 어떻게” 모델명에 붙었는지 알 수가 없습니다. 따라서, gradient descent에 대해 먼저 알아보고, 왜 이게 GBM과 관련이 있는지 설명하고자 합니다.

Gradient Descent란?

그럼 gradient descent가 무엇인지 알아봅시다.

Gradient descent를 한국말로 풀면 경사 하강법입니다.

경사하강법은 미분 가능한 convex한 함수, 즉 $y=x^2$ 과 같은 오목 함수에서 $y$ 가 최소값 갖게 하는 $x$ 를 찾기 위한 알고리즘입니다. (위의 그림에서는 $x$ 가 $w$ 로 표시돼 있네요!)
이를 위해 Gradient Descent는 함수의 기울기 혹은 경사 (gradient)를 구하여 기울기가 낮은 쪽으로 $x$ 값을 계속 이동시켜서 $y$ 가 극소값을 갖도록 반복합니다.

참고로

극소값은 local minimum으로 가장 작은 함수값이라 말할 순 없지만 주변에 비해 작은 값
최소값은 global minimum으로 전체 범위에서 가장 작은 함수값

라는 점에서 차이를 보입니다.

구체적인 Gradient Descent 알고리즘은 다음과 같습니다.

초기값 $x_0$ 을 선택하고, 다음의 식을 반복적으로 실행해서 임의의 조건을 만족하면 종료하게 됩니다.

x_t = x_{t-1} -\eta \Delta f(x_{t-1}),\ t=1,2,...

이처럼 $t$ 번째 단계의 $x_{t}$ 는 이전 단계의 값 $x_{t-1}$ 에서 $\eta$ 와 $\Delta f(x_{t-1})$ 를 곱한 값을 빼주어 값을 지속적으로 업데이트합니다.

여기서 $\eta$ 는 학습률 (learning rate), 스텝의 크기 (step size)라 불리우는데, 한 번 업데이트할 때 얼마나 큰 크기로 업데이트할지 조정하는 값으로 사용자에 의해 지정되는 값입니다. 그리고 $\Delta f(x_{t-1})$ 이 바로 기울기 (gradient)에 해당합니다. 기울기는 함수를 1차 미분한 것과 같습니다.

즉 $t$ 번째 값의 $x_t$ 는

기울기 > 0 이면 $x_{t-1}$ 보다 더 작은 값으로 업데이트되고
기울기가 < 0 이면 $x_{t-1}$ 보다 큰 값으로 업데이트됩니다.

기울기 = 0일 때 보통 극값을 갖기 때문에 기울기 > 0이면 $x$ 가 극소값보다 오른쪽에 있다는 뜻이므로 학습률 * 기울기만큼 빼주어 왼쪽으로 이동시켜주는 것이고, <0이면 극소값보다 왼쪽에 있다는 뜻이므로 학습률 * 기울기만큼 더해주어 오른쪽으로 이동시켜주는 원리입니다.

이를 무한 반복하는 것은 아니고 $t$ 번째 값인 $x_t$ 와 $t-1$ 번째 값인 $x_{t-1}$ 번째 값이 거의 같으면 그 값이 수렴했음을 의미하기 때문에 반복을 멈춥니다.

예를 들어 $f(x) = x^4 - 3x^3 +2$ 에서 극값을 찾기 위해 Gradient Descent를 구현해봅시다.

x_t_1 = 0
x_t   = 6 #초기값
eta   = 0.01 #step size
precision = 1e-6 #정밀도

def delta_f(x):
    return 4 * x**3 - 9 * x**2

while abs (x_t - x_t_1) > precision:
    x_t_1 = x_t
    x_t = x_t_1 - eta * delta_f (x_t_1)

print (str(x_t))

여기서 precision부분은 “ $x_{t}$ 와 $t-1$ 번째 값인 $x_{t-1}$ 번째 값이 거의 같으면” 이라는 조건을 수식화해 표현한 것입니다.
이 두 값의 차가 precision값인 0.000001보다 작으면 반복을 멈추게 되어 있기 때문입니다.

Boosting과 Gradient Descent의 관계

자, 그래서 Gradient Descent는 알았는데… Boosting과는 어떻게 연관이 있는거죠?

예측값 $F(x)$ 과 실제 값 $y$ 의 차이를 나타내는 손실 함수 (Loss Function)를 MSE ^{Mean Squared Error}를 2로 나눈 것으로 정의하면, 이 손실 함수의 1차 미분값 (gradient)의 음수를 취한게 바로 잔차 (residual)가 됩니다.

손실 함수는 예측값과 실제값의 차이를 나타내는 함수이기 때문에 값이 작을수록 좋습니다. 따라서 이 손실 함수의 극소값을 찾는 데 gradient descent를 이용한 최적화 문제로 치환할 수 있습니다.

손실 함수를 $L(y,F(x)) = (y - F(x))^2/2$ 로 정의하면 우리가 최소화해야 할 목적 함수 (Objective Function)은 $J = \sum_i L(y_i, F(x_i))$ 가 됩니다.

이 목적 함수에서 극소값을 찾기 위해서 1차 미분을 하면

\begin{aligned} \frac{\partial J}{\partial F(X_i)} &= \frac{\sum_i L(y_i, F(x_i))}{\partial F(x_i)} \\ &= \frac{\partial L(y_i, F(x_i))}{\partial F(x_i)} \\ &= F(x_i) - y_i \end{aligned}

따라서 잔차를 1차미분의 음수값 (negative gradients) 라 해석할 수 있습니다.

y_i - F(x_i) = - \frac{\partial J}{\partial F(x_i)}

그래서요…?

제가 GBM을 공부하면서 가장 막혔던 부분은 이 부분입니다. 손실 함수를 MSE/2로 정의하고, 이에 대한 gradient를 구하면 잔차와 관련 있는건 알겠는데…그래서 어떻게 gradient descent를 이용해서 GBM이 구성되는지는 이해가 되지 않았었습니다.

손실 함수를 최소화하기 위해 gradient descent를 사용하면

\begin{aligned} F_t(x_i) &= F_{t-1}(x_i) - \eta \underbrace{\frac{\partial J}{\partial F_{t-1}(x_i)}}_{F_{t-1} (x_i) - y_i},\ t=1,2,\ldots\\ F_t(x_i) &= F_{t-1}(x_i) + \eta \cdot \underbrace{(y_i - F_{t-1}(x_i))}_{잔차} \end{aligned}

와 같겠죠. 위의 Gradient descent 식에서

$x_t \Rightarrow F_t(x_i)$
$\Delta f(x_{t-1}) = \frac{\partial f(x)}{\partial x}\vert_{x=x_{t-1}} \Rightarrow \frac{\partial J}{\partial F_{t-1}(x_i)}$

로 바뀐 것 뿐입니다. 위에서 설명한 gradient descent 식은 수치적 최적화 (numerical optimization)이라면, 지금의 식은 함수 공간 (functional space)으로 최적화를 확장한 것입니다. 이것은 굳이 알 필요는 없습니다.

결국, $t$ 번째 예측값 $F_t(x_i)$ 는 ( $t-1$ 번째 예측값 + 학습률 * 잔차)로 업데이트됩니다.

자, 그럼 여기서 gradient descent로 예측값 $F_t(x_i)$ 을 업데이트해 MSE/2라는 손실함수를 최소화할 수 있다는 것을 알았습니다.
근데 위의 Introduction에서 분명 잔차에 약한 학습기를 적합시킨다 했는데 왜 굳이 잔차를 바로 구하지 않고 잔차를 예측해서 구할까요?

그 이유는 잔차는 데이터 포인트에서만 구할 수 있는 값인 반면, 약한 학습기에 잔차를 예측하게 되면 데이터 포인트 외의 모든 실수 범위에서 값을 구할 수 있기 때문입니다. 즉, $n$ 개의 데이터만 있으면 잔차인 $y_i - F_{t-1}(x_i)$ 도 $n$ 개만 존재하지만, 예측을 해서 구하면 어떤 $x_i$ 값을 넣든 간에 예측된 잔차값을 알 수 있어 "generalization"이 가능해집니다.

GBM, 자세히 뜯어보자

그래서 Boosting의 개념은 어디로 간거죠?

gradient descent로 예측값을 업데이트할 수 있고, 잔차를 약한 학습기에 적합시키는 것까지 이해했습니다.
근데 sequential하게 나무를 적합하고 additive하게 그 예측결과를 합치는 Boosting의 개념은 아직 나오지 않았습니다.

네, 그래서 이제 위에서 간략히 설명한 GBM의 설명인

$y_i, \ i=1,\ldots,n$ 에 대한 초기 예측값으로 $y_i$ 값들의 평균인 $\overline{y}$ 를 택합니다.
$t=1,\ldots, T$ $t = 1, \dots, T$ 만큼 반복
1. 각 관측치마다 잔차 $y_i - F_t(x)$ 를 구하고, 이 잔차에 첫 번째 약한 학습기를 적합시킵니다.
2. 적합된 결과에 따라 잔차의 예측값 (= 마지막 노드 (leaf)에 속한 잔차들의 평균)을 구하고
3. 잔차 예측값 * 학습률을 이전 예측값과 더해 새롭게 예측값을 업데이트합니다.

을 다시 차근차근 수식과 함께 (^^) 알아보고자 합니다.

모형 초기화: $F_0(x) = \overline{y}$
$t=1,\ldots, T$ $t = 1, \dots, T$ 만큼 반복:
1. $i$ 번째 관측치의 $t$ 번째 잔차 계산
  $r_{it} = - \left[ \frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}\right]_{F(x) = F_{t-1} (x)}$
2. $r_{it},\ i=1,\ldots,n$ 에 대해 나무모형을 적합하고 마지막 노드인 $R_{jt}$ (= terminal nodes, leaves)들에 대해서 잔차 예측값 $\gamma_{jt},\ j=1,\ldots, J_t$ 계산
3. $t$ 번째 예측값 계산
  $F_t (x) = F_{t-1} (x) + \eta \sum_{j=1}^{J_{t}} \gamma_{jt} I(x \in R_{jt})$

여기서 $R_{jt}$ 는 $t$ 번째 나무의 $j$ 번째 마지막 노드를 의미합니다. $t$ 번째 나무의 마지막 노드의 개수는 $J_t$ 개입니다. 그리고 $\gamma_{jt}$ 는 이 마지막 노드들의 예측값을 의미합니다. 저희는 손실 함수를 MSE/2꼴로 두었기 때문에 이 예측값들은 해당 노드에 속한 잔차들의 평균이 됩니다.

마지막으로 $y$ 의 $t$ 번째 예측값인 $F_t(x)$ 는 $t-1$ 번째 예측값 $F_{t-1} (x)$ 에 학습률 * 잔차의 예측값을 더한 것과 같습니다.

손코딩으로 알아보는 GBM

쉬운 예시를 통해서 이 알고리즘을 이해해보겠습니다.
몸무게를 예측하기 위한 변수로 키, 좋아하는 색, 성별을 두었다 가정합시다. 좋아하는 색, 성별은 범주형이기 때문에 get_dummies를 통해 수치형으로 변환해줬습니다.
참고로 이 데이터는 [link]에서 착안해서 비슷하게 만들었는데요. 해당 링크에 굉장히 자세하게 GBM을 설명하고 있습니다. 강추! 합니다. (~~자꾸 BAAM! 하고 노래부르고 그러지만… 나름 교수님이시고 가스펠도 하시는 능력자님…~~)

import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor

df = pd.DataFrame({
    "키": [1.6,1.6,1.5,1.8, 1.3, 1.6],
    "좋아하는 색":["파랑","초록","파랑","초록","분홍","파랑"],
    "성":["M","F","F","M","M","F"]
    ,"몸무게":[88,76,56,80,30,52]})
df = df.get_dummies(df)

첫 번째 업데이트입니다. 먼저 $y$ 값들의 평균으로 모형 초기화를 합니다.

F_0 (x) = 63.67

x = df.loc[:,df.columns!='몸무게']
y = df['몸무게']
F0 = np.mean(y)
print(F0)
> 63.666666666666664

첫 번째 잔차 $r_1 = y - F_0(x)$ 를 구하고

r1 = y - F0
print(r1)

> 0    24.333333
> 1    12.333333
> 2    -7.666667
> 3    16.333333
> 4   -33.666667
> 5   -11.666667

이 잔차에 깊이가 2인 나무모형을 적합해 잔차 예측값인 $\gamma$ 를 구합니다.

tree1  = DecisionTreeRegressor (random_state = 433, max_depth = 2)
r1_fit = tree1.fit(x, r1)
gamma1  = r1_fit.predict(x)
print ('잔차 예측값:{}'.format(np.unique(gamma1)))

> 잔차 예측값:[-33.66666667  -2.33333333  20.33333333]

이 예측값을 시각화하면 다음과 같습니다.

여기서 value에 해당하는 게 $\gamma_{j1}, j=1,2,3$ 값입니다.

마지막으로 예측값을 업데이트합니다.

eta = 0.1
F1 = F0 + eta * gamma1 
print ('첫번째 예측값:{}'.format(F1))
첫번째 예측값:[65.7        63.43333333 63.43333333 65.7        60.3        63.43333333]

여기서 eta는 학습률, gamma1는 첫번째 잔차를 예측한 값, F1는 새로 예측된 $y$ 값을 의미합니다. 예측값을 보면 예측 실력이 형편없죠. 이렇게 약한 학습기는 하나만으로는 좋은 예측 성능을 내지 못합니다.

이 과정을 내가 정한 반복수 $T$ 만큼 반복하면 다음과 같이 GBM을 구현할 수 있습니다.

F_new = np.mean(y) #초기화 모형
tree = DecisionTreeRegressor(random_state=433, max_depth =2) #약한 학습기
eta = 0.1 #학습률

for t in range(100): #T = 100
    print(str(t+1)+'번째 업데이트')
    F_old = F_new
    r     = y - F_old        # 잔차
    r_fit = tree.fit(x,r)    # 잔차 적합
    gamma = r_fit.predict(x)
    print('잔차 예측값:{}'.format(np.round(np.unique(gamma),2)))
    
    F_new = F_old + eta * gamma # 새로운 예측값
    
    print('예측값:{}'.format(np.round(F_new,2)))

이렇게 초반의 업데이트 값들은 잔차도 매우 크게 예측되고, 예측값도 실제 몸무게 값과 많이 다르지만 뒤로 갈수록 잔차 예측값은 0에 가깝고, 예측값도 실제 몸무게 값과 비슷한 결과를 냅니다.

GBM in Python

당연히 GBM도 sklearn패키지에 구현되어 있습니다.
다음은

위의 몸무게 데이터를
GradientBoostingRegressor로 모형을 적합하고,
10-fold 교차 검증을 통해 MSE를 계산한 코드입니다.

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold

# evaluate the model
model = GradientBoostingRegressor()
cv = RepeatedKFold(n_splits=3, n_repeats=3, random_state=1)
n_scores = cross_val_score(model, x, y, scoring='neg_mean_squared_error', cv=cv, n_jobs=-1, error_score='raise')
print('MSE: %.3f (%.3f)' % (-mean(n_scores), std(n_scores)))

> MSE: 885.179 (530.951)

GBM에서 튜닝을 위한 중요한 하이퍼 파라미터는 다음과 같습니다.

loss: 경사 하강법에서 사용할 손실 함수를 지정합니다. default 는 'deviance’입니다.
learning_rate: GBM이 학습을 진행할 때마다 적용하는 학습률입니다. default는 0.1입니다. learning_rate는 너무 작은 값이나 큰 값이면 반복이 끝나도 최소 잔차값을 찾아내지 못할 수 있습니다. 따라서 learning_rate는 n_estimators와 상호보완적으로 조합해 사용합니다. learning_rate를 작게 하고, n_estimators를 크게하면 더 이상 성능이 좋아지지 않는 한계점까지는 예측성능이 조금씩 좋아질 수 있습니다. 그러나 수행 시간이 길어지고, 예측 성능도 역시 현격하게 좋아지지는 않습니다.
n_estimators: 약한 학습기의 개수입니다. 약한 학습기가 순차적으로 오류를 보정하므로 개수가 많을수록 예측 성능이 좋아집니다. default는 100입니다.
subsample: weak learner가 학습에 사용하는 데이터 샘플링 비율입니다. default은 1이며, 이는 전체 학습 데이터를 기반으로 학습한다는 의미입니다. 과적합이 염려되는 경우 subsample을 조정해 1보다 작은 값으로 설정합니다.

이러한 하이퍼 파라미터는 GridSearchCV를 통해 최적화할 수 있습니다.

from sklearn.model_selection import GridSearchCV

params = {
    'n_estimators':[100,200,300]
    ,'learning_rate': [.05, .1]
}

model = GradientBoostingRegressor(random_state = 433)
grid_cv = GridSearchCV(model, param_grid = params, cv = 2, verbose=1)#, scoring='neg_mean_squared_error')
grid_cv.fit(x,y)
print('최적 하이퍼 파라미터:\n',grid_cv.best_params_)
print('최고 예측 정확도: {0:.4f}'.format(grid_cv.best_score_))

> 최적 하이퍼 파라미터:
>  {'learning_rate': 0.1, 'n_estimators': 100}
> 최고 예측 정확도: -1.5024

이렇게 GradientBoostingRegressor를 이용해 교차검증과 그리드 서치를 하는 방법까지 알아보았습니다.

다음 포스팅에서는 GBM을 분류에 어떻게 이용하는지나 XGBoost에 대해 더 알아보겠습니다 !

References

Deep Play님의 Gradient Boosting Algorithm의 직관적 이해 [link]
모두를 위한 컨벡스 최적화 [link]
StatQuest with Josh Starmer YouTube “Gradient Boost Part 2: Regression Details” [link]