2021데이터 크리에이터 캠프 대회를 준비하는 팀원들끼리 스터디를 진행하면서 bayesian optimization 주제를 내가 맞게 되었다. Bayesian optimization을 찾아보던 중 gaussian process에 기반한 방법이라는 것을 알게 되었고, gaussian process를 학습하였다. 이 글의 마지막 부분에서 아주 짧게 gaussian process가 bayesian optimization으로 어떻게 연결되는지 언급하겠다. 본 글의 내용은 작성자가 Cornell 대학교의 Kilian Q. Weinberger 교수님의 Machine Learning Lecture 26을 듣고 이해한 만큼의 내용만 담고 있다.

1. Back to the linear regression

linear regression에서 회귀계수를 추정하는 문제를 생각해보자.
회귀 계수를 추정할 때는 아래의 가정에서 크게 다음과 같이 두 가지 접근법이 있을 수 있다.

1. MLE(Maximum Likelihood Estimation)

• 회귀계수 w를 가정했을 때, training data(=$(x_1, y_1), (x_2, y_2), ... (x_n, y_n))$

• $w = argmax_wP(y_1,x_1,...,y_n,x_n|w) \newline \ \ \ \ =argmax_w∏_{i=1}^nP(y_i,x_i|w) \quad\ \quad \quad \quad (1)$

2. MAP(Maximum A Posteriori)

• 데이터가 주어졌을 때, 나올 likelihood가 가장 높은 w를 찾는 방법. 이때, w의 prior distribution은 정규분포를 따른다고 가정한다.
• $w = argmax_wP(w|y_1,x_1,...,y_n,x_n)\quad\ \quad \quad \quad (2)$

Data가 주어지면 -> w를 추정해서 -> 최종적으로 y를 예측한다.

여기서 아래와 같은 방식으로 생각을 전환할 수 있다.

Data가 주어졌을 때, w의 분포를 추정한 다음 그 분포를 이용해 y를 예측하지 말고, 애초에 test point(=y)의 분포를 모델링 하면 안될까?

2. Estimate the posterior predictive distribution directly

사실 우리의 최종 목표는 회귀계수인 w를 추정하는 것 보다는 test point의 feature vector가 주어졌을 때, y값을 예측하는 데에 있다. 그래서 위와 같은 접근 방식의 전환은 합리적이라 생각된다.

그런데, 위와 같은 접근 방식이 합리적으로 보이기는 하지만, 어떻게 과연 y를 바로 모델링할 수 있을까?

여기에서 바로 bayesian 접근 방식이 활용된다. MAP에서 가정했던 것과 마찬가지로 먼저 w의 prior distribution을 normal distribution으로 가정한다. normal distribution은 normal distribution의 conjugate prior이기 때문에 이와 같은 가정은 유용하다. conjugate priror는 같은 유형의 posterior distribution을 갖기 때문이다.


$P(Y∣D,X)=∫_wP(Y,w∣D,X)dw\newline\qquad=∫_wP(Y∣w,D,X)P(w∣D)dw\quad\ \quad \quad \quad (3)$

위의 식에서 D는 data의 약자로 주어진 데이터의 일체를 의미한다. 즉, training data의 y라벨 뿐만 아니라, test data의 feature vector까지 포함한 집합을 의미한다.

위의 (3)번 방정식을 보고 우리는 the posteior predictive distribution $Y∣D,X$

(3)번 방정식의 가장 우변을 살펴보자. $P(Y∣w,D,X)$

다음으로, 우변의 $P(w∣D)$

$P(w∣D) = \frac{P(D|W)P(W)}{z} \quad , z:normalizing \ factor \quad \quad \quad (4)$

여기서 $P(D|W) = ∏_{i=1}^nP(y_i|x_i;w)$

참고 : (3)번 방정식의 가장 우변에 있는 $∫_wP(Y∣w,D,X)P(w∣D)dw$

3. Fit the parameters $\mu$ and $\sum$

$Y∣D,X$

$Y|X,D \sim N(\mu, \sum)\quad \quad \quad (5)$

여기서 $\mu$는 간결함을 위해 0으로 가정하자. 사실 $\mu$는 큰 관심사가 아니다. $\mu$를 0으로 만들기 위해 단순히 데이터의 label의 평균을 빼버리는 것과 같은 방식을 쉽게 취해줄 수 있다. 다음과 같이 좀 더 명시적으로 표현하겠다.

$P(\begin{bmatrix} y_{1} \\ \vdots \\ y_{n} \\ y_{t} \end{bmatrix} | [x_1, \dots, x_n, x_t]) \sim N(0,\sum)\quad(6), \newline \ n: the\ number\ of\ training\ data\ ,\newline \ x_t : the\ feautre\ vector\ of\ test\ point$

공분산 행렬 $\sum$가 가져야할 성질들을 생각해보자. 우선 $\sum$ 행렬은 positive semi definite이어야 한다. 또한, $\sum$ 행렬의 성분은 feature vector가 서로 가까운 값을 가질 수록 값은 1에 가까워지고 거리가 먼 독립적인 feature vector들 간에는 0의 값을 반환해야 한다.

여기서 공분산 행렬의 성분으로 kernel이 쓰일 여지가 생기게 된다. 예를 들어, rbf커널의 경우 위에서 말한 성질을 만족하며 간단한 형태이므로 아래에 등장하는 kernel(=$K(x_i,x_j)$

$Σ_{ij}=e^{\frac{−∥x_i−x_j∥}{2σ^2}}\qquad(7)$

아래에서 활용될 기호를 정의하겠다.

$K:training\ set\ covariance,\newline K_*:training-test\ set\ covariance\newline K_{**}:test\ set\ covariance$

따라서, $\sum= \begin{pmatrix} K & K_{*} \\ K_{*}^T & K_{**} \end{pmatrix}$

이렇게 표현하고 나면 (6)에서 도출한 결과를 통해 아래의 (8)을 유도할 수 있다.

$P(y_t|y_1,...,y_n,x_1,...,x_n) \sim N(K_*^TK^{-1}y,K_{**}-K_{*}^TK^{-1}K_{*})\qquad\qquad(8)$

multivariate normal distribution을 conditional normal distribution으로 바꿔주는 아래와 같은 공식이 성립한다고 알려져 있다.

최종적으로 우리는 (8)번에 나타나있듯 the posterior predictive distribution의 확률분포를 모두 구하였다. 한 가지 흥미로운 점은 (8)번의 정규분포의 평균이 kernerl regression의 예측치와 동일하다는 점이다.

또한, Kilian Q. Weinberger교수님은 gaussian process regression의 가장 큰 장점은 uncertainty를 얘기할 수 있다는 점이라고 하였다. 우리는 단순히 예측값만 아는 것이 아니라 그 예측값의 확률 분포를 알기 때문에 uncertainty를 측정할 수 있다. 또한, 개인적으로는 이와 같이 uncertainty를 알 수 있다는 게 후에 bayesian optimization과 gaussian process가 연결되는 지점이라고 생각한다. bayesian optimization은 uncertainty가 더 많은 영역의 hyperparameter를 exploration하도록 설계될 수 있기 때문이다.

4. Connection to Bayesian Optimization

위의 그림(출처:https://ax.dev/versions/0.1.3/docs/bayesopt.html)은 Gaussian process와 bayesian optimization에 관해 학습을 하다 보면 자주 마주하게 되는 그림이다.

먼저 위에서 설명한 gaussian process의 연장선상에서 위 그림을 이해하기 위해서는 위에서 언급한 x들을 hyperparameter로 생각하고 y를 해당 hyper parameter에서의 모델 성능으로 생각할 필요가 있다.

예를 들어, RandomForest 모델에서 hyper parameter를 tuning하는 상황이라고 할 때, $x_i =(n\_estimators_i, max\_depth_i)\ for\ i=1,2,3...$

Bayesian Optimization을 이해하는데 필수적인 개념인 Acquisition function에 대해 알아보자.

Acqusition function의 정의는 다음과 같다.

Acquisition functions are mathematical techniques that guide how the parameter space should be explored during Bayesian optimization.
They use the predicted mean and predicted variance generated by the Gaussian process model. 

출처 : https://tune.tidymodels.org/articles/acquisition_functions.html

즉, Acqusition function은 hyper parameter space를 어떻게 탐색해나갈지 알려주는 technique인데, 그 과정에서 gaussian process에 의해 예측된 평균과 분산을 활용한다는 것이다.

hyper parameter를 탐색해 나갈 때는 아래와 같이 두 가지 전략이 고려될 수 있다. 아래 두 가지 전략은 서로 trade-off 관계에 있다.

<Bayesian Optimization에서 hyper parameter search시 고려되는 두가지 방향성>

• exploitation : 현재 가장 좋은 성능이 확인된 parameter 주변을 search하는데 초점을 맞춘다.

• exploration : uncertainty가 많이 남아있는 지역들을 search한다.

이 두가지 전략은 모두 나름대로 타당한 전략인데, exploitation의 경우 큰 함수값을 갖는 좋은 point들은 현재까지 확인된 point들 중에 가장 좋은값이 나온 곳 주변에 있을 가능성이 있기 때문이다. exploration의 경우 현재 정보가 없는 영역을 우리에게 알려준다. 따라서, 현재까지 관측한 point들보다 훨씬 좋은 point가 exploration한 곳에 있을 수 있다.

따라서, acquisition function은 이 두가지를 고려할 수 있게 고안되어야 하는데, 대표적인 acquisition function인 Expected Improvement에 대해 알아본다. Expected Improvement는 현재까지 관측된 hyper parameter의 조합 중 가장 좋은 성능을 보여주는 hyper parameter보다 더 가장 큰 성능의 개선을 가져다 줄 것으로 기대되는 hyper parameter를 찾아나서는 방법이다.

Expected Improvement의 공식은 다음과 같다.

$EI(x)=E[max(y−m_{opt},0)]$

$m_{opt} : From\ previous\ data,\ the\ best \ performance$

이와 같은 expected improvement 산식은 closed form으로 부분적분을 이용해 풀면 아래와 같은 식이 된다는 것이 알려져있다.

$\theta : tuning\ parameter$

$σ(θ) : the\ predicted\ standard\ error\ for\ the\ parameter\ combination\ theta$

$Φ(⋅):the\ cumulative\ standard\ normal$

위의 Expected Improvement의 공식에서 $δ(θ)$항은 exploitation에 해당하는 항이고 $\sigma(θ)$항은 exploration에 해당하는 항이다. 다시 말해, $δ(θ)$항은 새로 탐색할 parameter의 성능 평균이 현재까지의 가장 최고의 성능을 보여준 hyperparameter보다 얼마나 더 클지를 고려하는 항이다. 또한, $\sigma(θ)$항은 새로 탐색할 parameter가 가지고 있는 uncertainty를 고려하는 항이다. 위에서 말했듯, exploitation과 exploration은 모두 합리적인 전략이라고 할 수 있는데, Expected Improvement의 경우 이 두 가지 전략이 고려되고 있음을 알 수 있다. 바로 위 GIF는 acquisition function을 Expected Improvement로 하고 있는데, GP estimate이 높은 곳과 uncertainty가 큰 곳을 위주로 탐색이 먼저 이루어지고 있음을 알 수 있다.

5. Practice with Bayesian Optimization

아래의 코드는 xgboost의 hyper parameter를 tuning할 때, randomized search를 했을 경우와 Bayesian Optimization을 적용했을 경우에 CV score가 어떻게 다르고, 몇번의 iteration만에 그 수준을 달성하는지를 비교할 수 있는 코드(출처:http://krasserm.github.io/2018/03/21/bayesian-optimization/)이다. 아래의 그래프를 통해 알 수 있듯이 Bayesian Optimization 방법이 훨씬 적은 iteration으로 더 높은 모델 성능을 보여주게 하는 hyper parameter를 찾는다는 점을 알 수 있다.

from sklearn import datasets
from sklearn.model_selection import RandomizedSearchCV, cross_val_score

from scipy.stats import uniform
from xgboost import XGBRegressor

# Load the diabetes dataset (for regression)
X, Y = datasets.load_diabetes(return_X_y=True)

# Instantiate an XGBRegressor with default hyperparameter settings
xgb = XGBRegressor()

# and compute a baseline to beat with hyperparameter optimization 
baseline = cross_val_score(xgb, X, Y, scoring='neg_mean_squared_error').mean()


param_dist = {"learning_rate": uniform(0, 1),
              "gamma": uniform(0, 5),
              "max_depth": range(1,50),
              "n_estimators": range(1,300),
              "min_child_weight": range(1,10)}

rs = RandomizedSearchCV(xgb, param_distributions=param_dist, 
                        scoring='neg_mean_squared_error', n_iter=25)

# Run random search for 25 iterations
rs.fit(X, Y);



bds = [{'name': 'learning_rate', 'type': 'continuous', 'domain': (0, 1)},
        {'name': 'gamma', 'type': 'continuous', 'domain': (0, 5)},
        {'name': 'max_depth', 'type': 'discrete', 'domain': (1, 50)},
        {'name': 'n_estimators', 'type': 'discrete', 'domain': (1, 300)},
        {'name': 'min_child_weight', 'type': 'discrete', 'domain': (1, 10)}]

# Optimization objective 
def cv_score(parameters):
    parameters = parameters[0]
    score = cross_val_score(
                XGBRegressor(learning_rate=parameters[0],
                              gamma=int(parameters[1]),
                              max_depth=int(parameters[2]),
                              n_estimators=int(parameters[3]),
                              min_child_weight = parameters[4]), 
                X, Y, scoring='neg_mean_squared_error').mean()
    score = np.array(score)
    return score

optimizer = BayesianOptimization(f=cv_score, 
                                 domain=bds,
                                 model_type='GP',
                                 acquisition_type ='EI',
                                 acquisition_jitter = 0.05,
                                 exact_feval=True, 
                                 maximize=True)

# Only 20 iterations because we have 5 initial random points
optimizer.run_optimization(max_iter=20)


y_rs = np.maximum.accumulate(rs.cv_results_['mean_test_score'])
y_bo = np.maximum.accumulate(-optimizer.Y).ravel()

print(f'Baseline neg. MSE = {baseline:.2f}')
print(f'Random search neg. MSE = {y_rs[-1]:.2f}')
print(f'Bayesian optimization neg. MSE = {y_bo[-1]:.2f}')

plt.plot(y_rs, 'ro-', label='Random search')
plt.plot(y_bo, 'bo-', label='Bayesian optimization')
plt.xlabel('Iteration')
plt.ylabel('Neg. MSE')
plt.ylim(-5000, -3000)
plt.title('Value of the best sampled CV score');
plt.legend();

6. 참고 자료

• https://www.cs.cornell.edu/courses/cs4780/2018fa/lectures/lecturenote15.html
• http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.363.7289&rank=88&q=stochastic
• https://arxiv.org/pdf/1505.02965.pdf
• https://arxiv.org/pdf/1807.02811.pdf
• https://greeksharifa.github.io/bayesian_statistics/2020/07/12/Gaussian-Process/
• https://ax.dev/versions/0.1.3/docs/bayesopt.html
• https://tune.tidymodels.org/articles/acquisition_functions.html
• http://krasserm.github.io/2018/03/21/bayesian-optimization/
• https://www.johndcook.com/blog/2012/10/29/product-of-normal-pdfs/