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

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

 

 

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

 

$$\begin{array}{r}\text{(1)}& & y=f(\mathbf{x})+ϵ\\ & ϵ\sim \mathcal{N}(0,{\sigma }_n^2),\\ \\ & f(\mathbf{x})\sim \mathcal{GP}(\mu (\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime }))\end{array}$$

 

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

 

$$\begin{array}{r}\text{(2)}& & \mathbb{E}\left[y\right]=\mathbb{E}[f(\mathbf{x})+ϵ]=\mu (\mathbf{x})\\ & \mathbb{E}[(y-\mu (\mathbf{x}))(y^{\prime }-\mu ({\mathbf{x}}^{\prime }))]\\ \\ & =\mathbb{E}[(f(\mathbf{x})+ϵ-\mu (\mathbf{x}))(f({\mathbf{x}}^{\prime })+{ϵ}^{\prime }-\mu ({\mathbf{x}}^{\prime }))]\\ \\ & =\mathbb{E}[((f(\mathbf{x})-\mu (\mathbf{x}))(f({\mathbf{x}}^{\prime })-\mu ({\mathbf{x}}^{\prime }))]\\ \\ & +\mathbb{E}[(f(\mathbf{x})-\mu (\mathbf{x})){ϵ}^{\prime }]\\ \\ & +\mathbb{E}[ϵ(f({\mathbf{x}}^{\prime })-\mu ({\mathbf{x}}^{\prime }))]+\mathbb{E}[ϵ{ϵ}^{\prime })]\\ \\ & =k(\mathbf{x},{\mathbf{x}}^{\prime })+{\sigma }_n^2{\delta }_{i{i}^{\prime }}\end{array}$$

 

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

 

$$\begin{array}{}\text{(3)}& y\sim \mathcal{GP}(\mu (\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime })+{\sigma }_n^2{\delta }_{i{i}^{\prime }})\end{array}$$

 

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

 

$$\begin{array}{}\text{(4)}& \left[\begin{array}{c}{\mathbf{y}}_{1:m}\\ {\mathbf{y}}_{\ast }\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}\mu \\ {\mu }_{\ast }\end{array}\right],\left[\begin{array}{cc}K+{\sigma }_n^{2}I& K_{\ast }\\ K_{\ast }^T& K_{\ast \ast }\end{array}\right])\end{array}$$

 

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

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

 

$$\begin{array}{}\text{(5)}& p({\mathbf{y}}_{\ast }|{\mathbf{y}}_{1:m)}=\mathcal{N}({\mu }_{pos},\mathrm{\Sigma })\end{array}$$

 

여기서

 

$$\begin{array}{rl}& \mathrm{\Sigma }=K_{\ast \ast }-K_{\ast }^T{[K+{\sigma }_n^{2}I]}^{-1}{K}_{\ast }\\ \\ & {\mu }_{pos}={\mu }_{\ast }+K_{\ast }^T{[K+{\sigma }_n^{2}I]}^{-1}({\mathbf{y}}_{1:m}-\mu )\end{array}$$

 

이다.

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

예제로서 함수 $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$ 에서

 

$$\begin{array}{}\text{(6)}& y_i=\cos (x_i)+{ϵ}_i,{ϵ}_i\sim \mathcal{N}(0,{10}^{-4})\end{array}$$

 

로 주어진다고 가정한다.

 

# 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$ 으로 놓고, 커널은 다음 식으로 한다.

 

$$\begin{array}{}\text{(7)}& k(x,x^{\prime })=\exp (-\frac{1}{2{\lambda }^2}(x-x^{\prime }{)}^2)\end{array}$$

 

여기서 ${\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}}_{\ast }=\{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_{\ast }$ 와 $K_{\ast \ast }$ 는 다음과 같이 계산된다.

 

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

 

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

 

$$\begin{array}{r}\text{(8)}& & \mathrm{\Sigma }=K_{\ast \ast }-K_{\ast }^T{[K+{\sigma }_n^{2}I]}^{-1}{K}_{\ast }\\ & {\mu }_{pos}=K_{\ast }^T{[K+{\sigma }_n^{2}I]}^{-1}{\mathbf{y}}_{1:m}\end{array}$$

 

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

 

$$\begin{array}{}\text{(9)}& K+{\sigma }_n^{2}I=L{L}^T\end{array}$$

 

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

 

$$\begin{array}{}\text{(10)}& \mathrm{\Sigma }=K_{\ast \ast }-K_{\ast }^T{L}^{-T}{L}^{-1}{K}_{\ast }\end{array}$$

 

$L_{\ast }$ 를 $K_{\ast }=L{L}_{\ast }$ 가 되도록 계산한다면 위 식은 다음과 같이 된다.

 

$$\begin{array}{r}\text{(11)}& \mathrm{\Sigma }& =K_{\ast \ast }-L_{\ast }^T{L}^T{L}^{-T}{L}^{-1}L{L}_{\ast }\\ & =K_{\ast \ast }-L_{\ast }^T{L}_{\ast }\end{array}$$

 

 

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{array}{r}\text{(12)}& {\mu }_{pos}& =L_{\ast }^T{L}^T{L}^{-T}{L}^{-1}{\mathbf{y}}_{1:m}\\ & =L_{\ast }^T{L}^{-1}{\mathbf{y}}_{1:m}\end{array}$$

 

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

 

$$\begin{array}{r}\text{(13)}& & {\mu }_{pos}=L_{\ast }^T\mathbf{v}\\ & L\mathbf{v}={\mathbf{y}}_{1:m}\end{array}$$

 

 

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

 

식 (8)로 계산한 테스트 입력 ${\mathbf{x}}_{\ast }=\{x_1,...,x_p\}$ 에서의 ${\mathbf{y}}_{\ast }$ 의 평균값과 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()