[GP-2] GP 회귀 (GP Regression)

가우시안 프로세스 \(f(\mathbf{x})\) 의 관측값에는 노이즈가 포함되어 있다고 가정하는 것이 보다 실제적이다.

 

 

노이즈를 평균이 \(0\) 이고 분산이 \(\sigma_n^2\) 인 가우시안으로 모델링한다면 GP(Gaussian process) 측정 모델은 다음과 같다.

 

\[ \begin{align} & y=f(\mathbf{x})+ \epsilon \tag{1} \\ \\ & \ \ \ \ \ \epsilon \sim \mathcal{N} (0, \sigma_n^2 ), \\ \\ & \ \ \ \ \ f(\mathbf{x}) \sim \mathcal{GP}( \mu(\mathbf{x}), k(\mathbf{x}, \mathbf{x}' ) ) \end{align} \]

 

노이즈가 가우시안 프로세스와 독립이라고 가정하면, 가우시안 프로세스 \(y\) 의 평균과 공분산(covariance)은 다음과 같이 된다.

 

\[ \begin{align} & \mathbb{E} \left[ y \right] = \mathbb{E} \left[ f(\mathbf{x})+ \epsilon \right] = \mu (\mathbf{x}) \tag{2} \\ \\ & \mathbb{E} \left[ ( y-\mu (\mathbf{x}))(y'-\mu (\mathbf{x}')) \right] \\ \\ & \ \ \ \ \ = \mathbb{E} \left[ (f(\mathbf{x})+ \epsilon- \mu (\mathbf{x}))(f(\mathbf{x}')+\epsilon '- \mu (\mathbf{x}')) \right] \\ \\ & \ \ \ \ \ = \mathbb{E} \left[ ((f(\mathbf{x}) - \mu (\mathbf{x}))(f(\mathbf{x}' )- \mu (\mathbf{x}' )) \right] \\ \\ & \ \ \ \ \ \ \ \ \ \ + \mathbb{E} \left[ (f(\mathbf{x})- \mu (\mathbf{x})) \epsilon ' \right] \\ \\ & \ \ \ \ \ \ \ \ \ \ + \mathbb{E} \left[ \epsilon (f(\mathbf{x}' )- \mu(\mathbf{x}' )) \right] + \mathbb{E} \left[ \epsilon \epsilon ') \right] \\ \\ & \ \ \ \ \ = k( \mathbf{x}, \mathbf{x}' ) + \sigma_n^2 \delta_{ii'} \end{align} \]

 

여기서 \(\delta_{ii'}\) 는 크로넥커델타 (Kronecker's delta) 함수로서 \(\mathbf{x}=\mathbf{x}'\) 인 경우만 \(1\) 이다. 위 식을 간단히 표기하면 다음과 같다.

 

\[ y \sim \mathcal{GP} (\mu \left( \mathbf{x}), \ k(\mathbf{x}, \mathbf{x}' ) + \sigma_n^2 \delta_{ii'} \right) \tag{3} \]

 

이제 데이터셋 \(\mathcal{D}= \{(\mathbf{x}_i, y_i ), \ i=1, ... , m \}\) 이 주어졌을 때 입력 \(\{ \mathbf{x}_{m+1}, ... , \mathbf{x}_{m+p} \}\) 에 대응하는 출력 \(\mathbf{y}_* = \begin{bmatrix} y_{m+1} & ... & y_{m+p} \end{bmatrix}^T\) 의 조건부 확률밀도함수 \(p( \mathbf{y}_* \vert \mathbf{y}_{1:m} )\) 를 구해보자. 관련 랜덤벡터의 결합(joint) 확률분포는 다음과 같다.

 

\[ \begin{bmatrix} \mathbf{y}_{1:m} \\ \mathbf{y}_* \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mu \\ \mu_* \end{bmatrix} , \ \begin{bmatrix} K + \sigma_n^2 I & K_* \\ K_*^T & K_{**} \end{bmatrix} \right) \tag{4} \]

 

여기서 \( \mu = \mu (\mathbf{x}_{1:m} ), \ \mu_*= \mu (\mathbf{x}_{m+1:m+p} )\) 이다.

가우시안 랜덤벡터 특성에 의하면, \(\mathbf{y}_{1:m}\) 과 \(\mathbf{y}_*\) 가 결합 가우시안 분포를 가지면, 랜덤벡터 \(\mathbf{y}_*\) 의 조건부 확률밀도함수도 가우시안이며, 평균과 공분산은 다음과 같이 주어진다.

 

\[ p(\mathbf{y}_* \vert \mathbf{y}_{1:m)} = \mathcal{N} (\mu_{pos}, \ \Sigma ) \tag{5} \]

 

여기서

 

\[ \begin{align} & \Sigma = K_{**}-K_*^T \left[ K+ \sigma_n^2 I \right]^{-1} K_* \\ \\ & \mu_{pos} = \mu_*+ K_*^T \left[ K+ \sigma_n^2 I \right]^{-1} \left( \mathbf{y}_{1:m}-\mu \right) \end{align} \]

 

이다.

위 식은 주어진 데이터셋 \(\mathcal{D}= \{ (\mathbf{x}_i, y_i ), \ i=1,...,m \}\) 을 이용하여 추정값 \(\mathbf{y}_*\) 와 추정 확률을 계산해 주는 식이다.

예제로서 함수 \(g(x)=\cos⁡(x)\) 를 가우시안 프로세스 \(f(x)\) 로 추정해보도록 하겠다. 추정 대상인 \(g(x)\) 는 미지의 함수로 가정한다. 측정 노이즈는 평균이 \(0\), 분산이 \(\sigma_n^2=10^{-4}\) 인 가우시안으로 가정한다. 데이터셋은 \(\mathcal{D}= \{(x_i, y_i ),\ i=1,...,5 \}\) 로서 \(x_1=-4\), \(x_2=-3\), \(x_3=-2\), \(x_4=-1\), \(x_5=4\) 에서

 

\[ y_i= \cos⁡(x_i )+ \epsilon_i, \ \ \epsilon_i \sim \mathcal{N} (0, 10^{-4} ) \tag{6} \]

 

로 주어진다고 가정한다.

 

# true function
f = lambda x: np.cos(x).flatten()

 

s = 0.01       # noise std

# training points (given)
X = np.array([ [-4], [-3], [-2], [-1], [4] ])
m = X.shape[0]  # number of training points
y = f(X) + s*np.random.randn(m)  #(m,)

 

 

함수 \(g(x)\) 에 대한 정보는 전혀 없다고 가정하여 가우시안 프로세스의 평균함수는 \(\mu(x)=0\) 으로 놓고, 커널은 다음 식으로 한다.

 

\[ k(x, x' )= \exp \left( -\frac{1}{2 \lambda^2 } (x-x' )^2 \right) \tag{7} \]

 

여기서 \(\lambda^2 =1\) 로 놓는다.

 

# kernel
def kernel(a, b):
    lam2 = 1
    sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a, b.T)
    return np.exp(-.5 * sqdist / lam2)

 

그러면 식 (4)의 공분산 \(K\) 는 다음과 같이 계산할 수 있다.

 

K = kernel(X, X)

 

 

 

테스트 입력 \(\mathbf{x}_*= \{x_1, ... , x_p \}\) 는 \( [-5, \ 5 ]\) 범위에서 일정 간격으로 50개 지정한다.

 

p = 50         # number of test points

# points for prediction
Xstar = np.linspace(-5, 5, p).reshape(-1,1)

 

그러면 식 (4)의 공분산 \(K_*\) 와 \(K_{**}\) 는 다음과 같이 계산된다.

 

Kstar = kernel(X, Xstar)
K2star = kernel(Xstar, Xstar)

 

이제 식 (5)의 조건부 평균과 공분산을 계산하면 된다.

 

\[ \begin{align} & \Sigma = K_{**}-K_*^T \left[ K+ \sigma_n^2 I \right]^{-1} K_* \tag{8} \\ \\ & \mu_{pos} = K_*^T \left[ K+ \sigma_n^2 I \right]^{-1} \mathbf{y}_{1:m} \end{align} \]

 

식 (8)은 역행렬 계산이 포함되므로, 연산량과 수치오차를 고려하여 다음과 같이 촐레스키 분해(Cholesky decomposition)를 이용한다.

 

\[ K+ \sigma_n^2 I = L L^T \tag{9} \]

 

그러면 식 (8)의 공분산 식은 다음과 같이 된다.

 

\[ \Sigma = K_{**}-K_*^T L^{-T} L^{-1} K_* \tag{10} \]

 

\(L_*\) 를 \(K_*=LL_*\) 가 되도록 계산한다면 위 식은 다음과 같이 된다.

 

\[ \begin{align} \Sigma &= K_{**}-L_*^T L^T L^{-T} L^{-1} LL_* \tag{11} \\ \\ &=K_{**}-L_*^T L_* \end{align} \]

 

 

L = np.linalg.cholesky(K + s**2*np.eye(m))
Lstar = np.linalg.solve(L, Kstar)

Sig = K2star - np.dot(Lstar.T, Lstar)

 

또한 식 (8)의 평균식도 다음과 같이 된다.

 

\[ \begin{align} \mu_{pos} &= L_*^T L^T L^{-T} L^{-1} \mathbf{y}_{1:m} \tag{12} \\ \\ &= L_*^T L^{-1} \mathbf{y}_{1:m} \end{align} \]

 

여기서 \(L^{-1} \mathbf{y}_{1:m}=\mathbf{v}\) 로 놓으면, 위 식은 다음과 같이 된다.

 

\[ \begin{align} & \mu_{pos} = L_*^T \mathbf{v} \tag{13} \\ \\ & L \mathbf{v} = \mathbf{y}_{1:m} \end{align} \]

 

 

mu_pos = np.dot(Lstar.T, np.linalg.solve(L, y))

 

식 (8)로 계산한 테스트 입력 \(\mathbf{x}_*= \{x_1, ... , x_p \}\) 에서의 \(\mathbf{y}_*\) 의 평균값과 3-표준편차 (\(3 \sigma\)) 를 그리면 다음과 같다.

 

 

데이터가 밀집된 영역 (\(-4\) 에서 \(-1\) 사이)의 표준편차가 데이터가 없는 영역 (\(0\) 에서 \(4\) 사이) 보다 더 작은 것을 알 수 있다.

 

 

다음 그림은 평균함수가 \(\mu (x)=0\) 이고 공분산이 식 (7)인 가우시안 프로세스에서 10개의 샘플함수를 추출하여 그린 것이다. 처음에 가정했던 가우시안 프로세스 확률 정보를 이용한 것이다. 이를 사전 프로세스 (GP prior) 라고 한다.

 

 

다음 그림은 식 (8)로 계산된 평균함수와 공분산을 갖는 가우시안 프로세스에서 10개의 샘플함수를 추출하여 그린 것이다. 데이터셋을 이용하여 GP prior를 업데이트 하여 얻은 프로세스이므로 사후 프로세스 (GP posterior) 라고 한다.

 

 

다음 그림은 커널 식 (7)에서 \(\lambda^2=0.1\) 로 놓고 GP regression을 수행한 것이다.

 

 

다음 그림은 \(\lambda^2=10\) 일 때 GP regression을 수행한 것이다.

 

 

\(\lambda^2\) 이 작을 수록 데이터셋의 관련성이 작게 반영되므로 분산값이 매우 커지는 것을 알 수 있다. 반대로 이 값이 클수록 데이터셋의 관련성을 과대하게 반영하므로 분산값이 작지만 추정값에 큰 바이어스가 생기는 것을 알 수 있다. 이와 같이 GP regression 의 성능은 커널 파라미터 값에 큰 영향을 받는다. 따라서 데이터셋을 이용하여 커널 파라미터를 최적의 값으로 학습하는 방안을 강구하는 것이 필요하다.

다음은 본 예제에서 사용한 GP regression 의 전체 코드다.

 

# GP regression example
# by st.watermelon

import numpy as np
import matplotlib.pyplot as plt

# kernel
def kernel(a, b):
    lam2 = 1
    sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a, b.T)
    return np.exp(-.5 * sqdist / lam2)

# true function
f = lambda x: np.cos(x).flatten()

# parameters
p = 50         # number of test points
s = 0.01       # noise std

# training points (given)
X = np.array([ [-4], [-3], [-2], [-1], [4] ])
m = X.shape[0]  # number of training points
y = f(X) + s*np.random.randn(m)  #(m,)

K = kernel(X, X)
L = np.linalg.cholesky(K + s**2*np.eye(m))

# points for prediction
Xstar = np.linspace(-5, 5, p).reshape(-1,1)

# posterior mean (mu)
Kstar = kernel(X, Xstar)
Lstar = np.linalg.solve(L, Kstar)
mu_pos = np.dot(Lstar.T, np.linalg.solve(L, y))

# posterior covariance
K2star = kernel(Xstar, Xstar)
Sig = K2star - np.dot(Lstar.T, Lstar)

s2 = np.diag(Sig)
s = np.sqrt(s2)

# plotting
plt.figure(1)
plt.clf()
plt.plot(X, y, 'r+', ms=20)
plt.plot(Xstar, f(Xstar), 'b-')
plt.gca().fill_between(Xstar.flat, mu_pos-3*s, mu_pos+3*s, color="#dddddd")
plt.plot(Xstar, mu_pos, 'r--', lw=2)
plt.title('Mean predictions plus 3 std')
plt.axis([-5, 5, -4, 4])


# samples from the prior
L = np.linalg.cholesky(K2star + 1e-6*np.eye(p))
f_prior = np.dot(L, np.random.normal(size=(p,10)))
plt.figure(2)
plt.clf()
plt.plot(Xstar, f_prior)
plt.title('Ten samples from the GP prior')
plt.axis([-5, 5, -4, 4])


# samples from the posterior
L = np.linalg.cholesky(K2star + 1e-6*np.eye(p) - np.dot(Lstar.T, Lstar))
f_post = mu_pos.reshape(-1,1) + np.dot(L, np.random.normal(size=(p,10)))
plt.figure(3)
plt.clf()
plt.plot(Xstar, f_post)
plt.title('Ten samples from the GP posterior')
plt.axis([-5, 5, -4, 4])

plt.show()