[ML] Bayesian Optimization이란

Introduction (도입)

보통 우리가 실험을 하거나 입력을 넣고 결과가 나오는 현상에 대해서 이해를 하기 위해서 입력과 결과 사이의 관계를 파악하고 그리고 관계를 이용하여 공식을 도출하고 싶어 합니다. 너무나도 잘 알려진 공식 F=ma 또한 그렇게 나온 것이니...

이때 입력과 결과 사이에 관계를 표현하는 요소를 보통 Black Box 또는 Black box function이라고 표현합니다.

그리고 Black Box의 특징은 다음과 같습니다.

• closed form으로 표현되지 않으며 
• Non-Linear 또는 Non-Convex 하여 수렴이 잘 되지 않는다.
• 또한 상당히 Complex, Noisy 하다.
• 마지막으로 데이터를 한번 얻기 위해 또는 evaluate를 하기 위해서 수많은 비용과 시간이 드는 특징이 있어서 여러 번 수행하여 trial and error와 같은 접근 방법으로는 부적절하다.  
  (즉 y=x라는 함수였다면 수백만 번 테스트가 가능하겠지만 실제 화학물질 제조와 관련된 function이라면 검증을 위해서 수많은 cost가 들기 때문이다) 

이러한 function이 주어진 환경에서 우리가 원하는 답, 즉 global optimum을 찾기 위해서는 보통 2가지 방법이 존재합니다.

1. Grid Search: 관심 있는 하이퍼 파라미터들을 대상으로 가능한 모든 조합을 시도하는 것 

2. Random Search: Grid의 접근법이 처음부터 하나하나 다 시도해보는 것이라면 Random 하게 시도를 해보는 것

첫 번째 방법은 시간이 많이 소요가 될 가능성이 크며 Random은 운이 좋으면 시간을 단축시킬 수도 있겠지만 시간이 확실하게 줄어들 수 없는 가능성이 현저하게 존재합니다.

이럴 때 우리는 Bayesian Optimization을 사용합니다. (줄여서 BO라고 표현)

BO의 목표는 우리가 모르는 Objective (Black-box) function에 대해 알아내는, 또는 근사하는 것이고 근사가 되었을 때 함수에서 최대 (Maximum)또는 최소 (Minimum), 즉 원하는 optimal solution을 보여줄 거 같은 x 지점을 얻어내는 것이 목표인 알고리즘입니다.

활용분야는 다음과 같습니다. 

• Material Discovery & Drug Discovery
• Computer graphics
• Robotics
• Sensor networks
• Architecture configuration
• Particle physics 
• .... 

BO의 근간이 되는 중요한 베이즈 정리부터 설명드리겠습니다. 

 

Bayes theorem (베이즈 정리) 

Bayesian Optimization은 Black box의 정체를 알아내기 위해서 우선 Bayes Theorem을 사용합니다. 

Bayes Theorem은 조건부 확률을 계산하기 위한 접근법 중에 하나입니다. 공식은 아래 사진과 같으며 두 사건 A와 B가 주어졌을 때, 조건부 확률과 각각의 확률 정보를 사용하여 순서가 뒤바뀐 조건부 확률을 쉽게 구할 수 있는 정리입니다. 

즉 우리가 구하고자 하는 것은 좌변에 위치한 P(A|B)이며 우리가 알고 있는 정보들은 우항에 있는 3개의 요소들입니다. 

출처: https://science.howstuffworks.com/math-concepts/bayes-theorem.htm

좀더 더 들어가 봅시다. 

P(A|B): 사건 B라는 증거에 대한 사후 확률 (posterior probability) - 즉 사건 B가 일어났다는 것을 알고 있을 때 그것이 사건 A로부터 일어난 것이라고 생각되는 조건부 확률 

P(B|A): 사건 A가 주어졌을 때 B의 조건부 확률 (Likelihood) - 알려진 결과 (관찰된 결과)에 기초한 어떤 가설에 대한 가능성 (즉, 이 가설을 지지하는 정도) 

P(A): A의 사전 확률 (prior probability) - 과거의 경험 

P(B): B의 사전 확률 (evidence) - 현재의 증거 

여기서 사전 확률은 확률 실험 이전에 사건 발생에 대해 이미 알고 있는 지식을 의미하며 사후 확률은 어떤 사건을 인지한 후에 이들이 어떤 원인에 의해서 출현한 것이라고 생각되는 조건부 확률 지식을 의미합니다. 

우리는 Bayes 공식을 위에서 배운 개념을 토대로 다음과 같이 표현이 가능합니다. 

Posterior = likelihood * prior 

해당 식을 통해서 domain안에서 주어진 samples들을 이용하여 unknown 한 목표 함수에 대한 믿음을 정량화할 수가 있으며 우리는 목표 함수를 통해서 평가를 할 수 있는 중요한 framework를 만들 수가 있으며 Posterior Probability는 곧 Bayesian Optimization에서 중요한 요소 중 하나인 objective function의 surrogate를 만들 수가 있는 기반이 됩니다. 

즉 현재까지 얻은 모델 (Prior)과 추가적인 실험 정보 (Likelihood)를 통해 데이터가 주어졌을 때의 모델 (Posterior)을 추정해 나가는 방식입니다. 

그리고 BO에서는 각각의 요소를 다음과 같이 표현합니다.

Prior = Surrogate Function

Likelihood = Acquistion Function

각각에 대해서 설명드려보겠습니다. 

 

BO Component (1) - Surrogate Function (대체 함수)

Surrogate Function은 지금까지 얻어진 결과들을 토대로 미지의 함수, Objective Function에 대해 대략적인 값을 추정해주는 함수를 의미합니다. 

함수 생성을 위해서 보통 GP (Gaussian Process)를 많이 사용하며 생성 후에는 아래 그림에 나타내어져 있는 점선이 Surrogate Function입니다. 

*파란색 선: 목적함수, 파란 영역: 목적함수가 존재할만한 Confidence Bound (Function의 Variance) 

GP에 대한 구체적인 설명과 코드는 아래의 링크를 참조하시면 더 많은 정보를 찾을 수 있습니다.

1.7. Gaussian Processes

*다른 방법들
- BNN (Bayesian Neural Network)
- Random forest 

 

BO Component (2) - Acquisition Function (획득 함수) 

Acquisition Function의 역할은 predicted 된 평균과 uncertainty를 보여주는 분산 값을 기반으로 확률적으로 다음 추정 지점, 즉 최적의 값을 찾기 위해 알아봐야 할 후보를 제안해주는 함수입니다.

좀 더 정확히는 (2-1) 전체 지점에 대해서 Acquisition Function을 거쳐서 확률 또는 값으로 도출이 됩니다. 도출까지가 획득 함수의 역할이고 우리는 (2-2) 여기서 제일 높은 값이 또는 제일 낮은 값이 어디인지를 max와 같은 함수를 이용하여 다음 탐색해야 할 위치를 찾게 되는데 이러한 기능은 획득 함수와 함께하는 Utility 함수의 역할입니다. 

즉 (2-1)을 통해서 얻은 결과는 전체 지점에 대해서 Acquisition Function값이 얼마인지를 보여주는 초록색 영역이며 (2-2) 과정을 통해서 탐색해야 할 다음 지점이 최고 값을 가지는 별표 위치인 것을 알 수가 있습니다.

그리고 Acquisition Function에서는 중요한 2가지 요소가 있습니다. 

1. Exploration: 탐색이라는 의미이며 단어 의미 그 자체로 우리가 알지 못하는 영역에 대해 파악하는 것을 의미하며 불확실성이 높은 곳에서 지금까지 나온 값들 보다 더 큰 값 또는 작은 값이 나올 것이라는 방향을 가진 전략입니다.  

2. Exploitation: 이용이라는 의미이며 현재까지 가지고 있는 또는 관측된 값들 근처에서 더 큰 값이 나올 것이라는 방향을 가진 탐색 전략입니다.

무조건 아는 것만 이용하는 전략을 지향해서도 안되며 계속 모르는 지역만 탐색하면 안 되기에 Acquisition Function은 해당 2가지 개념을 적절하게 포함할 수 있는 식을 가지고 있으며 각 개념의 영향력을 조절할 수 있는 Hyperparameter가 곱해져 있기에 parameter값을 수정해나가면서 BO의 성능을 개선할 수도 있습니다. 

그러면 Acquisition Function 내부에 어떤 식으로 식이 구성되어 있으며 어떤 입력과 어떤 출력을 내는지, 그리고 종류가 몇 가지가 되는지에 대해서 궁금증이 생길 수가 있기에 밑에서 나열해보겠습니다. 

 

Acquisition Function 종류 

(1). GP-UCB (Gaussian Process - Upper Confidence Bound) 

UCB는 GP regression을 통해서 나온 Uncertainty에서 Upper bound의 최댓값을 다음 관측치로 결정하는 방법입니다. 식은 Mean + kappa * standard deviation을 의미하며 Mean과 Standard deviation은 Surrogate Function에서 온 값이며, 즉 관찰된 데이터들이 보여주는 평균과 분산입니다.

Kappa는 Hyperparameter로 Exploration과 Exploitation의 정도를 조절하는 역할을 하게 됩니다. 즉 kappa값이 클수록 deviation term에 영향력이 커지기에 Exploration을 더 하겠다는 전략이 되며 kappa값이 작아질수록 기존에 가지고 있는 정보를 더 사용하는 방향으로 전략이 바뀌게 됩니다.   

그림에서 보라 색선은 Surrogate Function에 의해 만들어진 결과이며 파란색 선이 UCB Acquisition Function 결과를 통해서 도출된 결과입니다. 여기서 최고치에 해당하는 값인 초록색 지점이 다음 탐색 추천 지점이 되겠습니다.

"""UCB 예시 함수 코드"""
def _ucb(x, gp, kappa):
 with warnings.catch_warnings(): 
       warnings.simplefilter("ignore") 
       mean, std = gp.predict(x, return_std=True) 
   return mean + kappa * std

 

(2). PI (Probability of Improvement) 

PI는 새로운 관측지가 가질 수 있는 mean이 현재 mean보다 더 클 확률에 대해 계산하는 방법입니다.

누적 분포 함수 (cumulative distribution function) 내부에 f_max는 현재까지 찾은 Solution 즉 최대 평균값이며 mean(x)와 sigma(x)는 위의 UCB에서의 요소와 동일하게 Surrogate Function가 가진 평균과 Deviation을 가지고 있고 Xi 또는 크사이는 Hyperparameter로써 또한 Exploration과 Exploitation의 비율을 조정하는 역할을 수행합니다. 

위의 사진을 보면 Xi의 크기에 따라서 탐색하는 (역삼각형) 위치가 다른 것을 볼 수가 있습니다. Xi가 커질수록 Exploration의 영향력이 커져서 Uncertainty가 큰 쪽으로 이동합니다. 

"""PI 예시 함수"""
def _poi(x, gp, y_max, xi):
   with warnings.catch_warnings(): 
      warnings.simplefilter("ignore") 
      mean, std = gp.predict(x, return_std=True) 
   z = (mean - y_max - xi)/std 
   return norm.cdf(z)

 

 

(3). EI (Expected Improvement) 

EI는 확률적으로 얼마나 개선될 것인가에 대한 가능성을 계산해주는 방법이며 현재의 max 보다 더 큰 x를 선택할 큰 함숫값을 도출할 확률과 그 함숫값과 최대 함수값 간의 차이 값을 이용하여 x의 유용성을 반환합니다. 

그림에서 현재까지 조사된 max 위치가 x(+)일 때 다음 후보 입력값 위치에 대해 조사를 하는데 x1, x2 x3에 대해서 조사하는 그림입니다. 

오른쪽 초록색 음영은 최대 함숫값 f(x(+))보다 큰 함숫값을 도출할 확률 PI(x3)를 보여주며 영역이 가장 크기 때문에 x3를 채택하는 것으로 결론이 내어진 그림입니다. 

f(x3)에 대한 평균과 f(x3)-f(x+)를 가중하여 EI 계산이 완료가 되는데 해당 과정의 의미는 max값보다 큰 함숫값을 도출할 가능성뿐만 아니라 실제로 f(x3)가 max f보다 얼마나 더 큰 값인지에 대해서도 반영을 하는 방법입니다. 

*물론 Xi라는 Hyperparameter가 동일하게 들어가 있으며 Exploration과 Exploitation강도를 조절하는 역할을 수행합니다. 

"""EI 예제 함수 코드"""
def _ei(x, gp, y_max, xi): 
   with warnings.catch_warnings(): 
      warnings.simplefilter("ignore")
      mean, std = gp.predict(x, return_std=True) 
      z = (mean - y_max - xi)/std 
   return (mean - y_max - xi) * norm.cdf(z) + std * norm.pdf(z)

 

---

여기까지가 BO를 구성하는 2개의 요소 (Surrogate, Acquisition)와 2가지 중요한 개념 (Exploration, Exploitation) 그리고 개념을 포함하는 Acquisition 함수들의 종류까지 나열해보았습니다. 

 

BO Algorithm Step (알고리즘 순서) 

(1). 초기 데이터 준비 (Observation Data) - Prior 
*만약 초기 데이터가 없을 때는 Random Sampling으로 임의의 점을 생성합니다.

(2). Surrogate Function 생성 - Posterior 
-Gaussian Process 또는 Bayesian Neural Network 등의 방법을 이용하여 현재 알고 있는 Sample들의 Surrogate Model을 생성합니다.

(3). Acquisition Function 계산 + argmax 함수 수행
-Acquisition Function은 Posterior 함수이며 계산을 수행한 다음에 argmax함수를 통해 다음 sample point를 제시해줍니다. 

(4). Evaluation
-제시된 점을 토대로 실험을 수행하고 해당 데이터를 기존 데이터 셋에 추가합니다. 증강된 데이터 셋은 다시 루프를 돌았을 때 Surrogate Model 생성에 사용됩니다. 

(5). Check the Stopping criteria 
-요구된 목표 또는 조건을 만족하였는지 체크하고 더 이상 실험을 멈출 것인지에 대한 방향을 제시해줍니다. 

(6). Restart from (1) - (5) 단계에서 멈추지 않았다면 다시 (1)에서 재 시작합니다. 

 

 

BO Code Example (실습)

import numpy as np
from matplotlib import pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

# added part
from math import sin
from math import pi
from numpy import arange
from numpy import vstack
from numpy import argmax
from numpy import asarray
from numpy.random import normal
from numpy.random import random
from scipy.stats import norm
from warnings import catch_warnings
from warnings import simplefilter



np.random.seed(1)

# objective function 
def f(x):
    """The function to predict."""
    return x * np.sin(x)


# ----------------------------------------------------------------------
#  First the noiseless case (1D Observation data)
#obs_x = np.atleast_2d([1.23, 3.45, 5.1, 6.0, 7.63, 8.2]).T
obs_x = np.array([1.23, 3.45, 5.1, 6.0, 7.63, 8.2]).T
obs_x = np.reshape(obs_x,(len(obs_x),1))


# Observations (1D)
obs_y = f(obs_x).ravel() 
obs_y = np.reshape(obs_y,(len(obs_y),1))


# Mesh the input space for evaluations of the real function, the prediction and
# its MSE
#eval_x = np.atleast_2d(np.linspace(0, 10, 1000)).T
eval_x = np.array(np.linspace(0, 10, 1000)).T
eval_x = np.reshape(eval_x,(len(eval_x),1))

# Instantiate a Gaussian Process model
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)

# Fit to data using Maximum Likelihood Estimation of the parameters
gp.fit(obs_x, obs_y)

# Make the prediction on the meshed x-axis (ask for MSE as well)
pred_y, sigma = gp.predict(eval_x, return_std=True)
pred_y = np.array(pred_y)
sigma = np.array(sigma)


# optimize the acquisition function. 
def opt_acquisition(opt_obs_x, opt_obs_y, model):
    # random search, generate random samples
    random_x = random(100) * 10
    random_x = random_x.reshape(len(random_x), 1)
    #random_x = np.atleast_2d(random_x)
    
    # calculate the acquisition function for each sample
    scores = acquisition(opt_obs_x, random_x, model)
    
    # locate the index of the largest scores
    ix = argmax(scores)
    return random_x[ix, 0]


# probability of improvement acquisition function
def acquisition(acq_obs_x, acq_random_x, model):
    # calculate the best surrogate score found so far
    y_hat, _ = surrogate(model, acq_obs_x)
    best = max(y_hat)
    
    # calculate mean and stdev via surrogate function
    mu, std = surrogate(model, acq_random_x)
    #print(mu)
    
    mu = mu[:, 0]
    
    # calculate the probability of improvement
    probs = norm.cdf((mu - best) / (std+1E-9))
    return probs

# surrogate or approximation for the objective function
def surrogate(model, X):
    # catch any warning generated when making a prediction
    with catch_warnings():
        # ignore generated warnings
        simplefilter("ignore")
        return model.predict(X, return_std=True)


plt.figure()
plt.plot(eval_x, f(eval_x), "r:", label=r"$f(x) = x\,\sin(x)$")
plt.plot(obs_x, obs_y, "b.", markersize=10, label="Observations")
plt.plot(eval_x, pred_y, "b-", label="Prediction")
plt.fill_between(eval_x.ravel(), pred_y.ravel()-1.96*sigma, pred_y.ravel()+1.96*sigma, color='red', alpha=0.5, label='95% Confidence Interval')
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.ylim(-10, 20)
plt.legend(loc="upper left")
    

for i in range(3):
    #1개 
    next_point_x = opt_acquisition(obs_x, obs_y, gp) 
    
    # sample the point
    next_actual_y = f(next_point_x)
    
    # summarize the finding
    est, _ = surrogate(gp, [[next_point_x]])
    print('>%d, x=%.3f, f()=%3f, actual=%.3f' %(i+1, next_point_x, est, next_actual_y))
    
    # add the data to the dataset
    obs_x = vstack((obs_x, [[next_point_x]]))
    obs_y = vstack((obs_y, [[next_actual_y]]))
    
    # update the model
    gp.fit(obs_x, obs_y)

    
pred_y, sigma = gp.predict(eval_x, return_std=True)
pred_y = np.array(pred_y)
sigma = np.array(sigma)

# Plot the function, the prediction and the 95% confidence interval based on
# the MSE

plt.figure()
plt.plot(eval_x, f(eval_x), "r:", label=r"$f(x) = x\,\sin(x)$")
plt.plot(obs_x, obs_y, "b.", markersize=10, label="Observations")
plt.plot(eval_x, pred_y, "b-", label="Prediction")
plt.fill_between(eval_x.ravel(), pred_y.ravel()-1.96*sigma, pred_y.ravel()+1.96*sigma, color='red', alpha=0.5, label='95% Confidence Interval')
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.ylim(-10, 20)
plt.legend(loc="upper left")

초기 Observation으로 만든 Surrogate Function	3번 iteration 후 결과

 

BO Opensource (오픈 소스)

BO는 상당히 역사가 오래된 알고리즘입니다. 때문에 찾을 수 있는 정보들도 많은데 대표적인 코드들만 나열하겠습니다.

Phoenics https://github.com/Atinary-technologies/phoenics

SMAC http://www.cs.ubc.ca/labs/beta/Projects/SMAC/

Hyperopt http://jaberg.github.io/hyperopt/

Bayesopt http://rmcantin.bitbucket.org/ 

Bayeso https://github.com/jungtaekkim/bayeso

 

Reference 

[1] Greenhill, Stewart, et al. "Bayesian optimization for adaptive experimental design: A review." IEEE access 8 (2020): 13937-13948.