15. Gaussian processes

15.1. Introduction

지도 학습에서는, 입력 \mathbf{x}_{i}와 출력 y_{i}에 대해 y_{i} = f(\mathbf{x}_{i})의 꼴로 어떠한 (노이즈가 낄 수 있는) 함수로 나타내어진다고 가정한다. 최적의 접근법은 이러한 함수들의 데이터에 대한 분포 p(f | \mathbf{X}, \mathbf{y})를 추론해 이를 이용해 새 입력에 대한 예측 p(y_{\ast} | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) = \int p(y_{\ast} | f, \mathbf{x}_{\ast}) p(f | \mathbf{X}, \mathbf{y}) df을 계산한다. 이 때 지금까지는 함수에 대한 추론을 할 때 매개변수에 대한 추론 p(\mathbf{\theta} | \mathcal{D})를 수행했다. 이제부터는 함수 자체에 대한 추론을 수행하도록 한다.

우리의 접근법은 가우시안 과정(GP)에 기반한다. 가우시안 과정은 함수들에 대한 사전분포를 정의한 뒤 관측된 데이터를 기반으로 이것을 함수들에 대한 사후분포로 변환시킨다. 함수에 대한 분포를 표현하기 어려워 보일 수 있겠지만 우리는 유한한 임의의 점들 \mathbf{x}_{1}, \cdots, \mathbf{x}_{n}에 대한 함수의 값들에 대한 분포를 정의할 수만 있으면 충분하다. 가우시안 과정은 p(f(\mathbf{x}_{1}, \cdots, f(\mathbf{x}_{N}))이 어떠한 평균 \mathbf{\mu}(\mathbf{x})와 원소 각각이 양의 정부호 커널로 정의된 공분산 \mathbf{\Sigma}(\mathbf{x})에 대해 결합적으로 가우시안임을 가정한다. 핵심 아이디어는 \mathbf{x}_{i}\mathbf{x}_{j}가 커널에 의해 유사하다고 판명이 되었다면, 이 점들에서의 함수값도 유사할 것이라고 기대하는 것에서 출발한다.

회귀 문제에서는 이 모든 계산은 닫힌 형태로 O(N^{3}) 시간에 될 수 있다. 분류 문제에서는 사후분포가 정확히 가우시안이 되지 않기 때문에 가우시안 근사 같은 근사법을 활용해야 한다.

가우시안 과정은 L1VM, RVM, SVM같은 커널법의 베이지안식 대안으로 볼 수 있다. 이런 커널법들이 더 희소하고 더 빠르긴 하지만, 확률적으로 잘 계기된 결과를 내지는 못한다. 여러 상황에서는 적절히 튜닝된 확률적 출력이 중요하다. 가우시안 과정의 유도법에는 고전적 방법이 있으며 크리깅 법이 있다.

15.2. GPs for regression

회귀시의 사전분포가 가우시안 과정을 따른다고 가정하자. 즉, 적당한 평균 함수 m(\mathbf{x})와 양의 정부호 커널 함수 \kappa(\mathbf{x}, \mathbf{x}^{\prime})에 대해 f(\mathbf{x}) \sim \mathrm{GP}(m(\mathbf{x}), \kappa(\mathbf{x}, \mathbf{x}^{\prime}))으로 나타내어지는 것이다. 이 경우 임의의 유한 개의 점에 대해 이 과정은 결합 가우시안 p(\mathbf{f} | \mathbf{X}) = \mathcal{N}(\mathbf{f} | \mathbf{\mu}, \mathbf{K})을 정의한다. (\mathbf{\mu} = (m(\mathbf{x}_{1}), \cdots, m(\mathbf{x}_{N})), K_{ij} = \kappa(\mathbf{x}_{i}, \mathbf{x}_{j}))

가우시안 과정은 평균 함수에 대해 유연성을 갖기 때문에 m(\mathbf{x}) = 0으로 종종 정의하고는 한다.

15.2.1. Predictions using noise-free observations

f_{i} = f(\mathbf{x}_{i})가 노이즈 없는 관측일 때 학습 데이터셋 \mathcal{D} = \{(\mathbf{x}_{i}, f_{i}) : i = 1, \cdots, N\}을 잡자. 테스트 셋 \mathbf{X}_{\ast}에 대한 출력 \mathbf{f}_{\ast}은 어떻게 예측할 수 있을까?

f(\mathbf{x})을 예측할 때 \mathbf{x}가 이미 관측된 값이라면 우리가 가우시안 과정을 쓸 때는 f(\mathbf{x})가 불확실성 없이 예측될 수 있도록 기대할 것이다. 즉, 관측에 노이즈가 없을 경우 가우시안 과정은 학습 데이터를 통한 외삽으로 기능한다.

가우시안 과정의 정의에 따라 결합분포는 다음의 형태를 가지게 된다. (여기서 \mathbf{K} = \kappa(\mathbf{X}, \mathbf{X}), \mathbf{K}_{\ast} = \kappa(\mathbf{X}, \mathbf{X}_{\ast}), \mathbf{K}_{\ast \ast} = \kappa(\mathbf{X}_{\ast}, \mathbf{X}_{\ast}))

\begin{pmatrix} \mathbf{f} \\ \mathbf{f}_{\ast} \end{pmatrix} \sim \mathcal{N}(\begin{pmatrix} \mathbf{\mu} \\ \mathbf{\mu}_{\ast} \end{pmatrix}, \begin{pmatrix} \mathbf{K} & \mathbf{K}_{\ast} \\ \mathbf{K}_{\ast}^{T} & \mathbf{K}_{\ast \ast} \end{pmatrix})

이 경우 사후분포는 다음과 같다.

p(\mathbf{f}_{\ast} | \mathbf{X}_{\ast}, \mathbf{X}, \mathbf{f}) = \mathcal{N}(\mathbf{f}_{\ast} | \mathbf{\mu}_{\ast}, \mathbf{\Sigma}_{\ast}), \mathbf{\mu}_{\ast} = \mathbf{\mu}(\mathbf{X}_{\ast}) + \mathbf{K}_{\ast}^{T} \mathbf{K}^{-1} (\mathbf{f} - \mathbf{\mu}(\mathbf{X})), \mathbf{\Sigma}_{\ast} = \mathbf{K}_{\ast \ast} - \mathbf{K}_{\ast}^{T} \mathbf{K}^{-1} \mathbf{K}_{\ast}.

당연하지만 어떤 커널을 쓰느냐에 따라 가우시안 과정도 달라지는데 자주 쓰이는 커널 중 하나로는 제곱지수커널(가우시안 커널, RBF 커널) \kappa(x, x^{\prime}) = \sigma_{f}^{2} e^{-\frac{1}{2 l^{2}} (x - x^{\prime})^{2} }이 있다. 여기서 l은 함수의 가로 방향 분산을, \sigma_{f}^{2}는 세로 방향 분산을 제어한다.

노이즈 없는 가우시안 과정 회귀의 여러 매개변수에 대한 예시.

노이즈 없는 가우시안 과정 회귀의 용례는 복잡한 시뮬레이션에 대해 계산 복잡도가 적게 드는 근사가 있다.

15.2.2. Predictions using noisy observations

이제 노이즈가 있는 관측 y = f(\mathbf{x}) + \epsilon, \epsilon \sim \mathcal{N}(0, \sigma_{y}^{2})에 대해 알아보자. 이 경우 모델은 데이터를 정확히 외삽하는 대신에 데이터에 근접하도록 기대된다. 공분산은 \mathrm{cov}[\mathbf{y} | \mathbf{X}] = \mathbf{K} + \sigma_{y}^{2} \mathbf{I}_{N} = \mathbf{K}_{y}가 된다.

이 때 관측 데이터와, 테스트 데이터에 대한 노이즈 없는 잠재 예측함수의 결합 확률밀도함수는 다음과 같다. (편의성을 위해 평균은 0으로 가정한다)

\begin{pmatrix} \mathbf{y} \\ \mathbf{f}_{\ast} \end{pmatrix} \sim \mathcal{N}(\mathbf{0}, \begin{pmatrix} \mathbf{K}_{y} & \mathbf{K}_{\ast} \\ \mathbf{K}_{\ast}^{T} & \mathbf{K}_{\ast \ast} \end{pmatrix})

따라서 사후예측분포는 다음과 같다.

p(\mathbf{f}_{\ast} | \mathbf{X}_{\ast}, \mathbf{X}, \mathbf{y}) = \mathcal{N}(\mathbf{f}_{\ast} | \mathbf{\mu}_{\ast}, \mathbf{\Sigma}_{\ast}), \mathbf{\mu}_{\ast} = \mathbf{K}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{y}, \mathbf{\Sigma}_{\ast} = \mathbf{K}_{\ast \ast} - \mathbf{K}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{K}_{\ast}.

단일 테스트 입력의 경우에는 다음과 같다. (\mathbf{k}_{\ast} = [\kappa(\mathbf{x}_{\ast}, \mathbf{x}_{1}), \cdots, \kappa(\mathbf{x}_{\ast}, \mathbf{x}_{N})], k_{\ast \ast}=\kappa(\mathbf{x}_{\ast}, \mathbf{x}_{\ast}))

p(f_{\ast} | \mathbf{X}_{\ast}, \mathbf{X}, \mathbf{y}) = \mathcal{N}(f_{\ast} | \mathbf{k}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{y}, k_{\ast \ast} - \mathbf{k}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{k}_{\ast})

사후평균의 경우 \mathbf{\alpha} = \mathbf{K}_{y}^{-1} \mathbf{y}로 놓으면 \bar{f}_{\ast} = \sum_{i=1}^{N} \alpha_{i} \kappa(\mathbf{x}_{i}, \mathbf{x}_{\ast})로 쓸 수도 있다.

15.2.3. Effect of the kernel parameters

가우시안 과정의 성능은 커널의 선택에 따라 달라진다. 노이즈 낀 관측에 대해 다음의 제곱지수 커널을 사용한다고 가정하자: \kappa_{y}(x_{p}, x_{q}) = \sigma_{f}^{2} e^{-\frac{1}{2 l^{2}} (x_{p} - x_{q})^{2}} + \sigma_{y}^{2} \mathbf{1}_{p = q} 이 때 l은 함수의 가로 방향 분산을, \sigma_{f}^{2}는 세로 방향 분산을, \sigma_{y}^{2}는 노이즈 분산을 제어한다.

여러 커널 매개변수에 대한 노이즈 낀 1D 가우시안 과정.

제곱지수커널을 다차원으로 다음과 같이 확장할 수 있다: \kappa_{y}(\mathbf{x}_{p}, \mathbf{x}_{q}) = \sigma_{f}^{2} e^{-\frac{1}{2} (\mathbf{x}_{p} - \mathbf{x}_{q})^{T} \mathbf{M} (\mathbf{x}_{p} - \mathbf{x}_{q})} + \sigma_{y}^{2} \mathbf{1}_{p = q} 이 때 행렬 \mathbf{M}은 다양한 방법으로 정의할 수 있는데, 가장 간단한 방법은 등장 행렬 \mathbf{M}_{1} = l^{-2} \mathbf{I}로 정의하는 것이고, 아니면 각각의 차원마다 스케일을 달리 해서 \mathbf{M}_{2} = \mathrm{diag}(\mathbf{l})^{-2}로 잡을 수도 있다. 이 때 스케일이 커지는 차원의 경우 해당 특성은 사실상 무시된다. D \times K 행렬 \mathbf{\Lambda}에 대해 \mathbf{M}_{3} = \mathbf{\Lambda} \mathbf{\Lambda}^{T} + \mathrm{diag}(\mathbf{l})^{-2}로 잡는 방법도 있는데, 이를 인자 분석 거리 함수라 한다.

여러 초매개변수에 대한 2D 가우시안 과정.

15.2.4. Estimating the kernel parameters

커널 매개변수를 추정하기 위해서는 SVM처럼 매개변수 각각의 조합으로 만들어지는 그리드에 대해 검증 손실 함수를 목적함수로 해서 탐색을 할 수도 있는데 이는 매우 느리다. 그 대신에 연속성 있는 최적화법을 제공하는 실측 베이스법을 쓸 수 있다. 즉, 주변가능도 p(\mathbf{y} | \mathbf{X}) = \int p(\mathbf{y} | \mathbf{f}, \mathbf{X}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f})을 최대화시키는 것이다.

p(\mathbf{f} | \mathbf{X}) = \mathcal{N}(\mathbf{f} | \mathbf{0}, \mathbf{K}), p(\mathbf{y} | \mathbf{f}) = \prod_{i} \mathcal{N}(y_{i} | f_{i}, \sigma_{y}^{2})이므로, 주변가능도는 \log p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2} \mathbf{y} \mathbf{K}_{y}^{-1} \mathbf{y} - \frac{1}{2} \log \lvert \mathbf{K}_{y} \rvert - \frac{N}{2} \log (2 \pi)으로 주어진다. 첫 번째 항은 데이터에 대한 피팅을 나타내고, 두 번째 항은 모델의 복잡도를 나타내고, 세 번째 항은 상수항이다. 데이터에 대한 피팅과 모델의 복잡도는 트레이드오프 관계에 있는데, 데이터에 대해 피팅이 잘 되어 첫 번째 항이 작아지면 \mathbf{K}는 거의 대각행렬이 되므로 두 번째 항은 커진다. 이는 모델의 복잡도가 높아짐을 나타낸다. 반대로 모델의 복잡도를 낮추면 \mathbf{K}는 거의 단위행렬이 되므로 두 번째 항은 작아지지만 길이 비례상수가 커지므로 첫 번째 항은 커져서 모델은 잘 피팅되지 않게 된다.

주변가능도를 최대화시키려면 어떻게 해야 할까? 커널 매개변수를 \mathbf{\theta}로 놓으면 \mathbf{\alpha} = \mathbf{K}_{y}^{-1} \mathbf{y}일 때 \frac{\partial}{\partial \theta_{j}} \log p(\mathbf{y} | \mathbf{X}) = \frac{1}{2} \mathbf{y}^{T} \mathbf{K}_{y}^{-1} \frac{\partial \mathbf{K}_{y}}{\partial \theta_{j}} \mathbf{K}_{y}^{-1} \mathbf{y} - \frac{1}{2} \mathrm{tr}(\mathbf{K}_{y}^{-1} \frac{\partial \mathbf{K}_{y}}{\partial \theta_{j}}) = \frac{1}{2} \mathrm{tr}((\mathbf{\alpha} \mathbf{\alpha}^{T} - \mathbf{K}_{y}^{-1}) \frac{\partial \mathbf{K}_{y}}{\partial \theta_{j}} )로 계산된다. \mathbf{K}_{y}^{-1}을 계산하는 데는 O(N^{3}) 시간이, 그라디언트를 계산하기 위한 초매개변수 각각을 계산하는 데는 O(N^{2}) 시간이 걸린다. 첫째 항 \frac{\partial \mathbf{K}_{y}}{\partial \theta_{j}}는 커널의 형태와 미분하는 매개변수에 따라 달린다. 초매개변수에 \sigma_{y}^{2} \geq 0과 같은 제한이 걸릴 때가 많은데, 이 경우 \theta = \log (\sigma_{y}^{2})로 놓고 연쇄법칙을 쓴다. 이 주변가능도는 아무 경사 하강법을 이용해 최적화할 수 있으나 목적함수가 볼록함수가 아니기 때문에 국소적 해에 갇히는 문제가 있다.

15.2.4.1. Example

지수제곱커널을 최적화할 때에는 적은 데이터로는 국소적 해에 갇히기 쉽다.

지수제곱커널을 사용한 가우시안 과정. 2가지의 국소적 해가 존재한다.

15.2.4.2. Bayesian inference for the hyper-parameters

초매개변수에 대한 점 근사를 수행하는 대신 그 사후분포를 계산할 수 있다. \mathbf{\theta}의 차원수가 적다면 가능한 값의 조합에 대한 격자를 만들 수 있다. 이 격자는 최대사후확률추정 \hat{\mathbf{\theta}}을 중심으로 한다. 그 이후 잠재 변수에 대한 사후분포는 \delta_{s}가 격자점 s에 대한 가중치라고 할 때 p(\mathbf{f} | \mathcal{D}) \propto \sum_{s=1}^{S} p(\mathbf{f} | \mathcal{D}, \mathbf{\theta}_{s}) p(\mathbf{\theta}_{s} | \mathcal{D}) \delta_{s}로 근사할 수 있다.

매개변수의 차원수가 커지면 격자법은 차원의 저주에 시달린다. 몬테카를로 법이 대안이 될 수 있지만 느리다는 단점이 있다. 다른 대안은 유사 몬테카를로 법을 쓰는데, 최빈값에 격자점을 놓고 각각 차원마다 최빈값 양옆으로 1개씩의 격자점을 놓아 2 \lvert \mathbf{\theta} \rvert + 1개의 격자점만 쓰는 것이다. 이를 중심 복합 설계라 한다. 이 근사를 조금 더 정확하게 하기 위해서 초매개변수를 로그 변환하기도 한다.

15.2.4.3. Multiple kernel learning

커널 매개변수를 최적화하는 또 다른 방법으로 복수 커널 학습이 있다. 이는 커널을 기본 커널의 가중치합 \kappa(\mathbf{x}, \mathbf{x}^{\prime}) = \sum_{j} w_{j} \kappa_{j}(\mathbf{x}, \mathbf{x}^{\prime})로 나타낸 뒤 커널 자체가 아니라 가중치항 w_{j}를 최적화하는 것이다. 이는 다른 종류의 데이터를 융합시킬 때 유용하다.

15.2.5. Computational and numerical issues

사후예측평균값은 \bar{f}_{\ast} = \mathbf{k}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{y}로 쓸 수 있는데, \mathbf{K}_{y}를 직접 역행렬 변환하면 계산 복복잡도가 높아지므로 촐스키 분해 \mathbf{K}_{y} = \mathbf{L} \mathbf{L}^{T}를 계산하고 사후예측평균과 분산, 로그주변가능도를 계산할 수 있다. 이 때 \mathbf{\alpha} = \mathbf{L}^{-T} \mathbf{L}^{-1} \mathbf{y}, \mathbb{E}[f_{\ast}] = \mathbf{k}_{\ast}^{T} \mathbf{\alpha}, \mathbf{v} = \mathbf{L}^{-1} \mathbf{k}_{\ast}, \mathrm{var}[f_{\ast}] = \kappa(\mathbf{x}_{\ast}, \mathbf{x}_{\ast}) - \mathbf{v}^{T} \mathbf{v}, \log p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2} \mathbf{y}^{T} \mathbf{\alpha} - \sum_{i} \log L_{ii} - \frac{N}{2} \log 2 \pi이 된다. 촐스키 분해를 계산하는 데 O(N^{3}) 시간, \mathbf{\alpha}를 계산하는 데 O(N^{2}) 시간, 각 테스트 케이스마다 평균 계산에 O(N) 시간, 분산 계산에 O(N^{2}) 시간이 걸린다.

다른 대안은 \mathbf{K}_{y} \mathbf{\alpha} = \mathbf{y}을 켤레 그라디언트를 사용해 푸는 것이다. 루프 k번 이후에 멈추면 O(kN^{2}) 시간이 걸린다. 또 다른 방법은 켤레 그라디언트에 필요한 행렬-벡터곱을 빠른 가우시안 변환을 이용해 구할 수 있는데, 이 방법은 고차원 입력으로 확장할 수 없다는 단점이 있다.

15.2.6. Semi-parametric GPs

가우시안 과정의 평균값에 대한 선형 모델 f(\mathbf{x}) = \mathbf{\beta}^{T} \mathbf{\phi}(\mathbf{x}) + r(\mathbf{x})이 유용할 수 있는데, 이 때 r(\mathbf{x}) \sim \mathrm{GP}(0, \kappa(\mathbf{x}, \mathbf{x}^{\prime}))은 나머지 항을 모델링한다. 이는 매개변수 모델과 비매개변수 모델의 접근법 각각을 혼합하는 반매개변수 모델이다. \mathbf{\beta} \sim \mathcal{N}(\mathbf{b}, \mathbf{B})로 놓으면, 전체 매개변수는 새 가우시안 과정으로 모델링된다:

f(\mathbf{x}) \sim \mathrm{GP}(\mathbf{\phi}(\mathbf{x})^{T} \mathbf{b}, \kappa(\mathbf{x}, \mathbf{x}^{\prime}) + \mathbf{\phi}(\mathbf{x})^{T} \mathbf{B} \mathbf{\phi}(\mathbf{x}^{\prime}))

여기서 \mathbf{\beta}를 적분해 빼내면 테스트 입력 \mathbf{X}_{\ast}에 대한 사후예측분포는 다음과 같아진다.

p(\mathbf{f}_{\ast} | \mathbf{X}_{\ast}, \mathbf{X}, \mathbf{y}) = \mathcal{N}(\bar{\mathbf{f}}_{\ast}, \mathrm{cov}[f_{\ast}])

\bar{\mathbf{f}}_{\ast} = \mathbf{\Phi}_{\ast}^{T} \bar{\mathbf{\beta}} + \mathbf{K}_{\ast}^{T} \mathbf{K}_{y}^{-1} (\mathbf{y} - \mathbf{\Phi} \bar{\mathbf{\beta}})

\bar{\mathbf{\beta}} = (\mathbf{Phi}^{T} \mathbf{K}_{y}^{-1} \mathbf{\Phi} + \mathbf{B}^{-1})^{-1} (\mathbf{\Phi} \mathbf{K}_{y}^{-1} \mathbf{y} + \mathbf{B}^{-1} \mathbf{b})

\mathrm{cov}[\mathbf{f}_{\ast}] = \mathbf{K}_{\ast \ast} - \mathbf{K}_{\ast}^{T} \mathbf{K}_{y}^{-1} \mathbf{K}_{\ast} + \mathbf{R}^{T} (\mathbf{B}^{-1} + \mathbf{\Phi} \mathbf{K}_{y}^{-1} \mathbf{\Phi}^{T})^{-1} \mathbf{R}

\mathbf{R} = \mathbf{\Phi}_{\ast} - \mathbf{\Phi} \mathbf{K}_{y}^{-1} \mathbf{\Phi}_{\ast}

이 때 사후평균은 선형 모델의 출력에 가우시안 과정의 보정 항을 더한 것이 된다. 사후 공분산은 가우시안 과정의 공분산에 \mathbf{\beta}의 불확실성으로 인해 발생하는 항을 더한 것이다.

15.3. GPs meet GLMs

가우시안 과정을 일반화된 선형 모델 세팅으로 확장해 분류에 써 보자. 베이지안 로지스틱 회귀와 같이, 가장 어려운 부분은 가우시안 사전분포가 베르누이/다항 베르누이 가능도와 켤레함수가 아니라는 점이다. 여러 방법이 있으나 여기서는 가장 빠르고 간단한 가우시안 근사를 알아보기로 한다.

15.3.1. Binary classification

이진 분류의 경우 모델은 p(y_{i} | \mathbf{x}_{i}) = \sigma(y_{i} f(\mathbf{x}_{i}))의 형태가 된다. 가우시안 과정 회귀와 같이, f \sim \mathrm{GP}(0, \kappa)로 놓는다.

15.3.1.1. Computing the posterior

표준화되지 않은 사후분포의 로그값은 다음과 같이 정의된다.

l(\mathbf{f}) = \log p(\mathbf{y} | \mathbf{f}) + \log p(\mathbf{f} | \mathbf{X}) = \log p(\mathbf{y} | \mathbf{f}) - \frac{1}{2} \mathbf{f}^{T} \mathbf{K}^{-1} \mathbf{f} - \frac{1}{2} \log \lvert \mathbf{K} \rvert - \frac{N}{2} \log 2 \pi.

목적 함수 J(f) = -l(f)로 놓으면 그라디언트와 헤시안은 다음과 같다.

\mathbf{g} = - \nabla \log p(\mathbf{y} | \mathbf{f}) + \mathbf{K}^{-1} \mathbf{f}

\mathbf{H} = -\nabla \nabla \log p(\mathbf{y} | \mathbf{f}) + \mathbf{K}^{-1} = \mathbf{W} + \mathbf{K}^{-1}

이 때 데이터가 동일하고 독립적으로 분포되므로 \mathbf{W} = - \nabla \nabla \log p(\mathbf{y} | \mathbf{f})은 대각행렬이다.

최대사후확률추정은 반복 재가중 최소 제곱법으로 구할 수 있다. 이 때 업데이트는 다음과 같이 이루어진다.

\mathbf{f}^{\mathrm{new}} = \mathbf{f} - \mathbf{H}^{-1} \mathbf{g} = (\mathbf{K}^{-1} + \mathbf{W})^{-1} (\mathbf{W} \mathbf{f} + \nabla \log p(\mathbf{y} | \mathbf{f}))

이가 수렴하면 사후확률의 가우시안 근사는 p(\mathbf{f} | \mathbf{X}, \mathbf{y}) \simeq \mathcal{N}(\hat{\mathbf{f}}, (\mathbf{K}^{-1} + \mathbf{W})^{-1})의 형태가 된다.

15.3.1.2. Computing the posterior predictive

테스트 케이스 \mathbf{x}_{\ast}의 사후예측분포의 평균값은 다음과 같다.

\mathbb{E}[\mathbf{f}_{\ast} | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}] = \int \mathbb{E}[\mathbf{f}_{\ast} | \mathbf{f}, \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}] p(\mathbf{f} | \mathbf{X}, \mathbf{y}) d \mathbf{f} = \int \mathbf{k}_{\ast}^{T} \mathbf{K}^{-1} \mathbf{f} p(\mathbf{f} | \mathbf{X}, \mathbf{y}) d \mathbf{f} = \mathbf{k}_{\ast}^{T} \mathbf{K}^{-1} \mathbb{E}[\mathbf{f} | \mathbf{X}, \mathbf{y}] \simeq \mathbf{k}_{\ast}^{T} \mathbf{K}^{-1} \hat{\mathbf{f}}

사후예측분산은 다음과 같다.

\mathrm{var}[f_{\ast}] = \mathbb{E}[\mathrm{var}[f_{\ast} | \mathbf{f}]] + \mathrm{var}[\mathbb{E}[f_{\ast} | \mathbf{f}]] = k_{\ast \ast} - \mathbf{k}_{\ast}^{T} (\mathbf{K}^{-1} - \mathbf{K}^{-1} \mathrm{cov}[\mathbf{f}] \mathbf{K}^{-1}) \mathbf{k}_{\ast} \simeq k_{\ast \ast} - \mathbf{k}_{\ast}^{T} (\mathbf{K} + \mathbf{W}^{-1})^{-1} \mathbf{k}_{\ast}

따라서 p(f_{\ast} | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) = \mathcal{N}(\mathbb{E}[f_{\ast}], \mathrm{var}[f_{\ast}])이 된다.

이를 이진 분류 세팅에 적용하면 \pi_{\ast} = p(y_{\ast} = 1 | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) \simeq \int \sigma(f_{\ast}) p(f_{\ast} | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) d f_{\ast}이 된다. 이는 베이지안 로지스틱 회귀에 썼던 여러 방법들로 근사될 수 있다.

15.3.1.3. Computing the marginal likelihood

라플라스 근사를 사용해 주변가능도를 근사하면 \log p(\mathbf{y} | \mathbf{X}) \simeq \log p(\mathbf{y} | \hat{\mathbf{f}}) - \frac{1}{2} \hat{\mathbf{f}}^{T} \mathbf{K}^{-1} \hat{\mathbf{f}} - \frac{1}{2} \log \lvert \mathbf{K} \rvert - \frac{1}{2} \log \lvert \mathbf{K}^{-1} + \mathbf{W} \rvert 이 된다. \hat{\mathbf{f}}, \mathbf{W}\mathbf{\theta}에 의존하므로 도함수 계산은 회귀 세팅에 비해 더 복잡하다.

15.3.1.4. Numerically stable computation

수치적으로 안정된 계산을 하기 위해서는 \mathbf{K}\mathbf{W}의 역행렬을 직접 계산하지 않는 것이 좋다. 대신에 \mathbf{B} = \mathbf{I}_{N} + \mathbf{W}^{\frac{1}{2}} \mathbf{K} \mathbf{W}^{\frac{1}{2}}로 정의하면 고유값들은 1 이상이고 1 + \frac{N}{4} \max_{ij} K_{ij} 이하이므로 안전하게 역행렬을 구할 수 있다. 이 때 반복 재가중 최소 제곱법 업데이트는 \mathbf{B} = \mathbf{L} \mathbf{L}^{T}를 촐스키 분해로 잡았을 때 다음과 같다:

\mathbf{f}^{\mathrm{new}} = \mathbf{K} \mathbf{a},

\mathbf{a} = \mathbf{b} - \mathbf{W}^{\frac{1}{2}} \mathbf{L}^{-T} \mathbf{L}^{-1} \mathbf{W}^{\frac{1}{2}} \mathbf{K} \mathbf{b},

\mathbf{b} = \mathbf{W} \mathbf{f} + \nabla \log p(\mathbf{y} | \mathbf{f}).

뉴턴법 반복을 T번 수행할 때 피팅은 O(TN^{3}) 시간 복잡도와 O(N^{2}) 공간 복잡도를 소모한다. 수렴 시에는 \mathbf{a} = \mathbf{K}^{-1} \hat{\mathbf{f}}이 되며 \log p(\mathbf{y} | \mathbf{X}) = \log p(\mathbf{y} | \hat{\mathbf{f}}) - \frac{1}{2} \mathbf{a}^{T} \hat{\mathbf{f}} - \sum_{i} \log L_{ii}로 로그 주변가능도를 계산할 수 있다.

사후예측분포의 경우 평균은 \mathbb{E}[f_{\ast}] = \mathbf{k}_{\ast}^{T} \nabla \log p(\mathbf{y} | \hat{\mathbf{f}}), 분산은 \mathbf{v} = \mathbf{L}^{-1} \mathbf{W}^{\frac{1}{2}} \mathbf{k}_{\ast}일 때 \mathrm{var}[f_{\ast}] = k_{\ast \ast} - \mathbf{v}^{T} \mathbf{v}이 된다. 피팅은 O(N^{3}) 시간이 걸리고, 예측은 테스트 케이스 각각에 대해 O(N^{2}) 시간이 걸린다.

15.3.1.5. Example

초매개변수를 학습시킨 뒤 분류를 하면 모델이 데이터에 대해 더 깐깐한 설명력을 갖추게 됨을 알 수 있다.

가우시안 과정을 통한 학습 전/후 이진 분류.

15.3.2. Multi-class classification

다중 클래스 분류에 대해서는 p(y_{i} | \mathbf{x}_{i}) = \mathrm{Cat}(y_{i} | \mathcal{S}(\mathbf{f}_{i})), \mathbf{f}_{i} = (f_{i1}, \cdots, f_{iC}), f_{.c} \sim \mathrm{GP}(0, \kappa_{c})이 된다. 즉 클래스마다 하나씩 잠재함수가 있고, 이 잠재 함수들은 사전적으로 독립이며, 각각 다른 커널을 쓸 수 있다. 이전과 마찬가지로 사후분포에는 가우시안 근사를 사용한다.

15.3.2.1. Computing the posterior

비표준화된 로그 사후분포는 다음과 같다.

l(\mathbf{f}) = -\frac{1}{2} \mathbf{f}^{T} \mathbf{K}^{-1} \mathbf{f} + \mathbf{y}^{T} \mathbf{f} - \sum_{i=1}^{N} \log(\sum_{c=1}^{C} e^{f_{ic}}) - \frac{1}{2} \log \lvert \mathbf{K} \rvert - \frac{CN}{2} \log (2 \pi)

\mathbf{f} = (f_{11}, \cdots, f_{NC})^{T}, \mathbf{y} = (y_{11}, \cdots, y_{NC})^{T}, \mathbf{K}_{c} = [\kappa_{c}(\mathbf{x}_{i}, \mathbf{x}_{j})]

그라디언트와 헤시안은 다음과 같다.

\nabla l = -\mathbf{K}^{-1} \mathbf{f} + \mathbf{y} - \mathbf{\pi}

\nabla \nabla l = -\mathbf{K}^{-1} - \mathbf{W}

\mathbf{W} = \mathrm{diag}(\mathbf{\pi}) - \mathbf{\Pi} \mathbf{\Pi}^{T}, \mathbf{\Pi} = [\mathrm{diag}(\mathbf{\pi}_{:1}), \cdots, \mathrm{diag}(\mathbf{\pi}_{:C})]^{T}

최빈값은 반복 재가중 최소 제곱법으로 다음과 같이 업데이트된다.

\mathbf{f}^{\mathrm{new}} = (\mathbf{K}^{-1} + \mathbf{W})^{-1} (\mathbf{W}\mathbf{f} + \mathbf{y} - \mathbf{\pi}).

직접 계산은 O(N^{3}C^{3}) 시간이 걸리나 O(CN^{3}) 시간으로 단축시킬 수 있다.

15.3.2.2. Computing the posterior predictive

사후예측분포는 앞서 이진 분류 세팅 때와 비슷하게 구할 수 있다. \mathbf{Q}_{\ast} = \begin{pmatrix} \mathbf{k}_{1}(\mathbf{x}_{\ast}) & \cdots & \mathbf{0} \\ \cdots & \cdots & \cdots \\ \mathbf{0} & \cdots & \mathbf{k}_{C}(\mathbf{x}_{\ast}) \end{pmatrix} 로 놓으면

\mathbb{E}[\mathbf{f}_{\ast}] = \mathbf{Q}_{\ast}^{T} (\mathbf{y} - \hat{\mathbf{\pi}}),

\mathrm{cov}[\mathbf{f}_{\ast}] = \mathrm{diag}(\mathbf{k}(\mathbf{x}_{\ast}, \mathbf{x}_{\ast})) - \mathbf{Q}_{\ast}^{T} (\mathbf{K} + \mathbf{W}^{-1})^{-1} \mathbf{Q}_{\ast}이 된다.

이를 계산하기 위해서는 다음의 근사를 사용한다.

p(y | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) \simeq \int \mathrm{Cat}(y | \mathcal{S}(\mathbf{f}_{\ast})) \mathcal{N}(\mathbf{f}_{\ast} | \mathbb{E}[\mathbf{f}_{\ast}], \mathrm{cov}[\mathbf{f}_{\ast}]) d \mathbf{f}_{\ast}

또는 몬테 카를로 법을 사용할 수도 있다.

15.3.2.3. Computing the marginal likelihood

주변가능도는 다음과 같이 근사되며, 다른 때와 비슷하게 수치적으로 최적화될 수 있다.

\log p(\mathbf{y} | \mathbf{X} ) \simeq -\frac{1}{2} \hat{\mathbf{f}}^{T} \mathbf{K}^{-1} \hat{\mathbf{f}} + \mathbf{y}^{T} \mathbf{f} - \sum_{i=1}^{N} \log (\sum_{c=1}^{C} e^{\hat{f}_{ic}}) - \frac{1}{2} \log \lvert \mathbf{I}_{C_{N}} + \mathbf{W}^{\frac{1}{2}} \mathbf{K} \mathbf{W}^{\frac{1}{2}} \rvert

15.3.2.4. Numerical and computational issues

뉴턴 재가중 반복을 T번 수행할 때 피팅에는 O(TCN^{3}) 시간이, O(CN^{2}) 공간이 소모된다. 예측은 테스트 케이스 N_{\ast}개에 대해 O(CN^{3} + CN^{2} N_{\ast}) 시간이 소모된다.

15.3.3. GPs for Poisson regression

푸아송 회귀에도 가우시안 과정을 쓸 수 있다. 병증 매핑이 그 예인데, 모델의 형태는 상대 리스크 r_{i}에 대해 y_{i} \sim \mathrm{Poi}(e_{i}r_{i})의 꼴이 된다.

병증 매핑에 가우시안 과정 푸아송 회귀를 적용한 예.

15.4. Connection with other methods

다른 기계학습 방법과 가우시안 과정은 어떻게 관련되어 있을까?

15.4.1. Linear models compared to GPs

사전분포가 p(\mathbf{w}) = \mathcal{N}(\mathbf{0}, \mathbf{\Sigma})로 주어진 베이지안 선형 회귀를 고려해 보자. 사후예측분포는 다음과 같다:

p(f_{\ast} | \mathbf{x}_{\ast}, \mathbf{X}, \mathbf{y}) = \mathcal{N}(\mu, \sigma^{2})

\mu = \frac{1}{\sigma_{y}^{2}} \mathbf{x}_{\ast}^{T} \mathbf{A}^{-1} \mathbf{x}^{T} \mathbf{y}

\sigma^{2} = \mathbf{x}_{\ast}^{T} \mathbf{A}^{-1} \mathbf{x}_{\ast}

\mathbf{A} = \sigma_{y}^{-2} \mathbf{X}^{T} \mathbf{X} + \mathbf{\Sigma}^{-1}.

이를 풀면 베이지안 선형 회귀는 공분산 함수 \kappa(\mathbf{x}, \mathbf{x}^{\prime}) = \mathbf{x}^{T} \mathbf{\Sigma}\mathbf{x}^{\prime}를 차용한 가우시안 과정과 동일하다. 하지만, 이는 퇴화된 형태의 공분산 함수인데, 0이 아닌 고유값을 최대 D개까지밖에 가질 수 없기 때문이다. 즉, 표현할 수 있는 함수의 수가 제한되어 있음을 의미한다. 이는 과소적합을 불러일으킬 뿐더러, 사후분포 또한 유연하지 못하게 된다.

15.4.2. Linear smoothers compared to GPs

선형 다듬질가중치 함수 w_{i}(\mathbf{x}_{\ast})에 대해 학습 출력에 대한 선형 함수가 되는 회귀 함수 \hat{f} (\mathbf{x}_{\ast}) = \sum_{i} w_{i} (\mathbf{x}_{\ast}) y_{i}을 뜻한다. 이의 예는 커널 회귀, 국소 가중 회귀, 다듬질 쐐기, 가우시안 과정 회귀 등이 있다. 가우시안 과정 회귀가 선형 다듬질의 일부가 되는 이유는, 사후예측분포는 가중치 함수 w_{i}(\mathbf{x}_{\ast}) [(\mathbf{K} + \sigma_{y}^{2} \mathbf{I}_{N})^{-1} \mathbf{k}_{\ast}]_{i}에 대해 \bar{f}(\mathbf{x}_{\ast}) = \sum_{i=1}^{N} y_{i} w_{i}(\mathbf{x}_{\ast})의 꼴로 나타낼 수 있기 때문이다.

커널 회귀에서는 가중치 함수를 메르세르 커널 대신 다듬질 커널로부터 유도할 수 있다. 이 경우엔 가중치 함수가 국소적 지지집합을 갖게 된다. 가우시안 과정에 대해서는 가중치 함수가 \mathbf{K}^{-1}에 의존하게 되므로 가중치 함수의 지지집합은 명확하지 않아진다. 특정한 가우시안 과정 커널 함수에 대해서는 w_{i}(\mathbf{x})을 해석적으로 유도할 수 있는데, 이를 등가 커널이라 한다. 이 때 \sum_{i=1}^{N} w_{i}(\mathbf{x}_{\ast})=1이 됨을 보일 수 있는데 (가중치가 음수일 수 있으므로 볼록 결합은 아니다), 가우시안 과정에서 쓰였던 원본 커널 함수가 국소적이지 않더라도 각각의 가중치 함수는 국소적이 될 수 있다. 또한, 커널 다듬질을 할 때 대역폭을 수동으로 조절해야 했던 반면에, 가우시안 과정의 등가 커널의 유효 대역폭은 표본 크기가 커짐에 따라 자동으로 작아진다.

15.4.2.1. Degrees of freedom of linear smoothers

선형 다듬질은 왜 ‘다듬질’이라는 이름이 붙는가? 이는 가우시안 과정으로 잘 설명할 수 있다. 학습 데이터에 대한 예측값 \bar{\mathbf{f}} = \mathbf{K}(\mathbf{K} + \sigma_{y}^{2})^{-1} \mathbf{y}을 고려해 보자. 고유값 분해 \mathbf{K} = \sum_{i=1}^{N} \lambda_{i} \mathbf{u}_{i} \mathbf{u}_{i}^{T}를 생각해 보면 \mathbf{K}가 실수 대칭 양의 정부호 행렬이므로 고유값 \lambda_{i}는 음이 아닌 실수가 되며 고유벡터 \mathbf{u}_{i}는 정규직교한다. \gamma_{i} = \mathbf{u}_{i}^{T} \mathbf{y}, \mathbf{y} = \sum_{i=1}^{N} \gamma_{i} \mathbf{u}_{i}로 놓으면 위의 식은 \bar{\mathbf{f}} = \sum_{i=1}^{N} \frac{\gamma_{i} \lambda_{i}}{\lambda_{i} + \sigma_{y}^{2}}\mathbf{u}_{i}로 다시 쓸 수 있다. 즉 \mathbf{y}의 고주파 성분들은 다듬질된다. 이 때 선형 다듬질 모델의 유효 자유도\mathrm{dof} = \mathrm{tr}(\mathbf{K}(\mathbf{K} + \sigma_{y}^{2} \mathbf{I})^{-1}) = \sum_{i=1}^{N} \frac{\lambda_{i}}{\lambda_{i} + \sigma_{y}^{2}}이 된다. 이는 그래프가 얼마나 울퉁불퉁한지를 나타낸다.

15.4.3. SVMs compared to GPs

이진 분류에 보조 벡터 기계를 쓸 때 목적 함수는 J(\mathbf{w}) = \frac{1}{2} \lVert \mathbf{w} \rVert^{2} + C \sum_{i=1}^{N} (1 - y_{i}f_{i})_{+}이 된다. 이 때 최적해는 \mathbf{w} = \sum_{i} \alpha_{i} \mathbf{x}_{i}의 형태가 된다. 이를 커널화하면 \mathbf{f} = \mathbf{K} \mathbf{\alpha}에 대해 \lVert \mathbf{w} \rVert^{2} = \mathbf{f}^{T} \mathbf{K}^{-1}\mathbf{f}를 얻는다. 즉 보조 벡터 기계 목적 함수는 J(\mathbf{w}) = \frac{1}{2} \mathbf{f}^{T} \mathbf{f} + C \sum_{i=1}^{N} (1 - y_{i} f_{i})_{+}로 쓸 수 있고, 이는 가우시안 과정 분류기의 최대사후확률추정 목적함수 J(\mathbf{w}) = \frac{1}{2} \mathbf{f}^{T} \mathbf{f} - \sum_{i=1}^{N} \log p(y_{i} | f_{i})와 닮은 구석이 있다.

보조 벡터 기계를 적절한 가능도함수를 찾아 가우시안 과정으로 끼워맞출 수 있지 않겠냐는 물음을 던질 수 있지만 안타깝게도 그런 가능도함수는 없다. 유사가능도 함수는 존재한다.

15.4.4. L1VM and RVMs compared to GPs

희소 커널 기계들은 선형 모델에 기저 함수 확장 \mathbf{\phi}(\mathbf{x}) = [\kappa(\mathbf{x}, \mathbf{x}_{1}), \cdots, \kappa(\mathbf{x}, \mathbf{x}_{N})]을 적용한 것이다. 이는 다음 커널 \kappa(\mathbf{x}, \mathbf{x}^{\prime}) = \sum_{j=1}^{D} \frac{1}{\alpha_{j}} \phi_{j}(\mathbf{x}) \phi_{j} (\mathbf{x}^{\prime})을 적용하고 사전분포 p(\mathbf{w}) = \mathcal{N}(\mathbf{0}, \mathrm{diag}(\alpha_{j}^{-1}))을 적용한 가우시안 과정과 같다.

이 커널 함수에는 2가지 흥미로운 점이 있는데, 첫째로는 이는 퇴화 커널이어서 0이 아닌 고유값을 최대 N개까지만 가질 수 있어 결합분포가 나타낼 수 있는 함수의 폭이 매우 제한되어 있다는 점이다. 둘째로는 커널이 학습 데이터에 의존하므로 학습 데이터를 기반으로 추론할 때 신뢰도가 과대 평가된다는 점이다.

15.4.5. Neural networks compared to GPs

신경망은 일반화된 선형 모델의 비선형 일반화이다. 이진 분류의 경우 신경망은 로지스틱 회귀 모델을 로지스틱 회귀 모델에 적용한 것이라 볼 수 있다.

p(y | \mathbf{x}, \mathbf{\theta}) = \mathrm{Ber}(y | \mathrm{sigm}(\mathbf{w}^{T} \mathrm{sigm}(\mathbf{V}\mathbf{x})))

은닉층 1개가 있는 회귀 신경망은 p(y | \mathbf{x}, \mathbf{\theta}) = \mathcal{N}(y | f(\mathbf{x}; \mathbf{\theta}), \sigma^{2}), f(\mathbf{x}) = b + \sum_{j=1}^{H} v_{j} g(\mathbf{x} ; \mathbf{u}_{j})로 나타낼 수 있다 (b는 편향, v_{j}는 출력 가중치, \mathbf{u}_{j}는 입력 가중치, g는 시그모이드나 tanh 등의 은닉층 활성 함수).

가중치에 대한 사전분포를 b \sim \mathcal{N}(0, \sigma_{b}^{2}), \mathbf{v} \sim \prod_{j} \mathcal{N}(v_{j} | 0, \sigma_{w}^{2}), \mathbf{u} \sim \prod_{j} p(\mathbf{u}_{j})로 나타낼 때 가중치 전체를 \mathbf{\theta}로 두면 다음을 얻는다:

\mathbb{E}_{\mathbf{\theta}} [f(\mathbf{x})] = 0

\mathbb{E}_{\mathbf{\theta}} [f(\mathbf{x}) f(\mathbf{x}^{\prime})] = \sigma_{b}^{2} + H \sigma_{v}^{2} \mathbb{E}_{\mathbf{u}}[g(\mathbf{x}; \mathbf{u}) g(\mathbf{x}^{\prime}; \mathbf{u})]

g의 상한/하한이 존재한다면 중심극한정리에 의해 은닉 층의 유닛 수가 무한대가 된다면 가우시안 과정으로 수렴함을 알 수 있다.

활성 함수 g(\mathbf{x}; \mathbf{u}) = \mathrm{erf}(u_{0} + \sum_{j=1}^{D} u_{j} x_{j})로 잡고, \mathbf{u} \sim \mathcal{N}(\mathbf{0}, \mathbf{\Sigma})로 두면 공분산 커널은 \tilde{\mathbf{x}} = (1, x_{1}, \cdots, x_{D})일 때 \kappa_{\mathrm{NN}}(\mathbf{x}, \mathbf{x}^{\prime}) = \frac{2}{\pi} \arcsin (\frac{2 \tilde{\mathbf{x}} \mathbf{\Sigma} \tilde{\mathbf{x}}^{\prime} }{\sqrt{(1 + 2 \tilde{\mathbf{x}}^{T} \mathbf{\Sigma} \tilde{\mathbf{x}}) (1 + 2 (\tilde{\mathbf{x}}^{\prime})^{T} \mathbf{\Sigma} \tilde{\mathbf{x}}^{\prime})}})가 된다. 이것이 시그모이드같이 양의 정부호 커널이 아닌 커널과는 다른 “진짜” 신경망 커널이다. 방사기저함수 커널과는 다르게, 이 커널로부터 샘플된 함수는 데이터로부터 멀리 떨어진다고 해서 0에 가까워지지는 않는다. 그 대신에 데이터로부터 멀어지면 같은 값으로 머물게 됨을 알 수 있다.

“진짜” 신경망 커널 함수의 여러 값에 대한 플롯.

신경망에서 은닉층 활성 함수를 g(\mathbf{x} ; \mathbf{u}) = e^{-\frac{ \lvert \mathbf{x} - \mathbf{u} \rvert^{2}}{2\sigma_{g}^{2}}}, \mathbf{u} \sim \mathcal{N}(\mathbf{0}, \sigma_{u}^{2} \mathbf{I})로 잡으면 대응되는 커널은 방사기저함수 커널이 된다.

15.4.6. Smoothing splines compared to GPs

다듬질 쐐기는 데이터를 매끄럽게 외삽하는 비매개변수법이다. 가우시안 과정의 일부이기도 하며, 입력이 1/2차원일 때 주로 쓰인다.

15.4.6.1. Univariate splines

기본 아이디어는 함수 f를 피팅할 때 너무 울퉁불퉁한 함수들을 징벌하는 항을 더해 데이터와의 괴리도를 최소화시키는 것이다. 함수의 m차 도함수를 징벌하면, 목적 함수는 J(f) = \sum_{i=1}^{N} (f(x_{i}) - y_{i})^{2} + \lambda \int (\frac{d^{m}}{dx^{m}} f(x))^{2} dx이 된다. 이 때의 해 f는 2m-1차 다항식을 이어붙인 구간적 다항 함수로, \mathcal{I} = \bigcup_{i=2}^{N} [x_{i-1}, x_{i}]일 때 f(x) = \sum_{j=0}^{m-1} \beta_{j} x^{j} + \mathbf{1}_{\mathcal{I}} (\sum_{i=1}^{N} \alpha_{i} (x - x_{i})_{+}^{2m-1}) + \mathbf{1}_{\mathcal{I}^{c}} (\sum_{i=1}^{N} \alpha_{i} (x - x_{i})_{+}^{m-1})이 된다. m = 2일 경우 3차 다항식을 이어붙인 3차 쐐기를 얻는다. 이는 능선 회귀를 통해 피팅할 수 있으며, \hat{\mathbf{w}} = (\mathbf{\Phi}^{T} \mathbf{\Phi} + \lambda \mathbf{I}_{N})^{-1} \mathbf{\Phi}^{T} \mathbf{y}이 된다. 이는 O(N) 시간에 피팅하는 것도 가능하다.

15.4.6.2. Regression splines

일반적으로, 이 다항식들을 K개의 매듭 \xi_{k}에 대해 고정시킬 수 있다. 이를 회귀 쐐기라 한다. 이는 다음의 기저 함수 확장을 사용하는 매개변수 모델이다:

f(x) = \beta_{0} + \beta_{1} x + \sum_{k=1}^{K} \alpha_{j} (x - \xi_{k})_{+}^{3}

매듭의 수와 위치를 결정하는 것은 보조 벡터 기계에서 하는 방식과 같다. 회귀 계수 \alpha_{j}l_{2} 정규화자를 두면 징벌된 쐐기가 된다.

15.4.6.3. The connection with GPs

3차 쐐기는 사전분포 p(\beta_{j}) \propto 1, 나머지 r(x) \sim \mathrm{GP}(0, \sigma_{f}^{2} \kappa_{\mathrm{sp}}(x, x^{\prime})), \kappa_{\mathrm{sp}}(x, x^{\prime}) = \int_{0}^{1} (x - u)_{+} (x^{\prime} - u)_{+} du일 때 함수 f(x) = \beta_{0} + \beta_{1} x + r(x)에 대한 최대사후확률추정과 같다.

이 때 커널 함수는 부자연스럽고, 유도되는 가우시안 과정에 의한 사후표본도 매끄럽지 않지만, 사후최빈값과 사후평균값은 매끄러운 함수가 된다. 즉, 정규화자를 두는 것이 항상 좋은 사전분포를 만드는 것은 아니다.

15.4.6.4. 2d input (thin-plate splines)

3차 쐐기는 2차원 입력으로 확장할 수 있는데, 다음 형태의 정규화자를 두면 된다:

\int \int [ (\frac{\partial^{2} f(x)}{\partial x_{1}^{2}})^{2} + 2 (\frac{\partial^{2} f(x)}{\partial x_{1} \partial x_{2}})^{2} + (\frac{\partial^{2} f(x)}{\partial x_{2}^{2}})^{2}   ]dx_{1} dx_{2}

이 때의 해는 \mathbf{\phi}_{i}(\mathbf{x}) = \eta(\lVert \mathbf{x} - \mathbf{x}_{i} \rVert), \eta_{z} = z^{2} \log (z^{2})일 때 f(x) = \beta_{0} + \mathbf{\beta}_{1}^{T} \mathbf{x} + \sum_{i=1}^{N} \alpha_{i} \mathbf{\phi}_{i} (\mathbf{x})이 된다. 이는 얇은 판 쐐기법으로 알려져 있다. 이는 가우시안 과정의 최대사후확률추정에 대응된다.

15.4.6.5. Higher-dimensional inputs

더 높은 차원에 대해서는 최적해를 해석적으로 풀기 힘들다. 하지만, 매개변수 회귀 쐐기 세팅에서는 기저 함수를 더 유연하게 정의할 수 있다. 한 방법은 1차원 기저 함수들의 외적으로 정의되는 텐서곱 기저를 이용하는 것이다. 예를 들어, 2차원 입력에 대해서는 f(x_{1}, x_{2}) = \beta_{0} + \sum_{m} \beta_{1m} (x_{1} - \xi_{1m})_{+} + \sum_{m} \beta_{2m} (x_{2} - \xi_{2m})_{+} + \sum_{m} \beta_{12m}(x_{1} - \xi_{1m})_{+} (x_{2} - \xi_{2m})_{+}로 정의할 수 있다.

데이터의 차원이 높아질수록 고차원 상호 작용을 모델링하긴 어려워진다. 피팅해야 할 매개변수가 너무 많아지기 때문이다. 이에 대해서는 다변수 적응적 회귀 쐐기(MARS) 같은 방법을 쓸 수 있다.

15.4.7. RKHS methods compared to GPs

다듬질 쐐기법은 함수의 도함수를 징벌하는 방법이었다. 이를 일반화시킬 수는 없을까?

메르세르의 정리에 의해 임의의 양의 정부호 커널 함수는 고유함수에 대해 \kappa(\mathbf{x}, \mathbf{x}^{\prime}) = \sum_{i=1}^{\infty} \lambda_{i} \phi_{i}(\mathbf{x}) \phi_{i} (\mathbf{x}^{\prime})로 나타낼 수 있다. 이 때 \phi_{i}재생산 커널 힐베르트 공간 (RKHS) \mathcal{H}_{k} = \{f : f(\mathbf{x}) = \sum_{i=1}^{\infty} f_{i} \phi_{i}(\mathbf{x}), \sum_{i=1}^{\infty} \frac{f_{i}^{2}}{\lambda_{i}} < \infty\}의 기저로 잡은 뒤 함수 f(\mathbf{x}) = \sum_{i=1}^{\infty} f_{i} \phi_{i} (\mathbf{x}), g(\mathbf{x}) = \sum_{i=1}^{\infty} g_{i} \phi_{i} (\mathbf{x})의 내적을 <f, g>_{\mathcal{H}} = \sum_{i=1}^{\infty} \frac{f_{i} g_{i}}{\lambda_{i}}로 잡으면, <\kappa(\mathbf{x}_{1}, \cdot), \kappa(\mathbf{x}_{2}, \cdot)>_{\mathcal{H}} = \kappa(\mathbf{x}_{1}, \mathbf{x}_{2})이 된다. 이를 재생산 특성이라 한다.

이 때 함수의 놈 \lVert f \rVert_{H}^{2} = <f, f>_{\mathcal{H}}일 때 다음 형태의 최적화 문제를 고려해 보자:

J(f) = \frac{1}{2 \sigma_{y}^{2}} \sum_{i=1}^{N} (y_{i} - f(\mathbf{x}_{i}))^{2} + \frac{1}{2} \lVert f \rVert_{H}^{2}

직관적으로 커널에 대해 복잡도가 높은 함수는 높은 놈값을 갖게 됨을 알 수 있는데, 표현하기 위한 고유함수의 수가 많아지게 되기 때문이다. 따라서 복잡도가 낮으면서 데이터를 잘 피팅하는 함수가 선호된다. 이 때 해는 f(\mathbf{x}) = \sum_{i=1}^{N} \alpha_{i} \kappa(\mathbf{x}, \mathbf{x}_{i})의 꼴이 됨을 알 수 있다. 이를 표현자 정리라 하며, 다른 볼록 손실함수에 대해서도 성립한다.

이를 대입하면 J(\mathbf{\alpha}) = \frac{1}{2 \sigma_{y}^{2}} \lvert \mathbf{y} - \mathbf{K} \mathbf{\alpha} \rvert^{2} + \frac{1}{2} \mathbf{\alpha}^{T} \mathbf{K} \mathbf{\alpha}을 얻는데, 이를 풀면 \hat{\mathbf{\alpha}} = (\mathbf{K} + \sigma_{y}^{2} \mathbf{I})^{-1}을 얻는다. 따라서 \hat{f}(\mathbf{x}_{\ast}) = \sum_{i} \hat{\alpha}_{i} \kappa(\mathbf{x}_{\ast}, \mathbf{x}_{i}) = \mathbf{k}_{\ast}^{T} (\mathbf{K} + \sigma_{y}^{2} \mathbf{I})^{-1} \mathbf{y}이 된다. 이는 가우시안 과정의 사후예측분포의 형태와 같다. 사실, 가우시안의 경우 평균과 최빈값이 같은 분포이므로, 재생산 커널 힐베르트 공간 정규화자를 적용한 선형 회귀는 가우시안 과정의 최대사후확률추정과 같다. 가우시안 과정 로지스틱 회귀에 대해서도 같은 명제가 성립한다.

15.5. GP latent variable model

또 다른 확률적 주성분 분석에 커널을 결합하는 방법으로 가우시안 과정 잠재 변수 모델(GP-LVM)이 있다. 이를 설명하기 위해서는 먼저 확률적 주성분 분석으로부터 시작한다.

p(\mathbf{z}_{i}) = \mathcal{N}(\mathbf{z}_{i} | \mathbf{0}, \mathbf{I})

p(\mathbf{y}_{i} | \mathbf{z}_{i}, \mathbf{\theta}) = \mathcal{N}(\mathbf{y}_{i} | \mathbf{W}\mathbf{z}_{i}, \sigma^{2} \mathbf{I})

이는 \mathbf{z}_{i}를 적분해 빼내고 \mathbf{W}를 최대화시킴으로써 최대가능도추정으로 피팅할 수 있다. 목적함수는 \mathbf{C} = \mathbf{W}\mathbf{W}^{T} + \sigma^{2} \mathbf{I}일 때 p(\mathbf{Y} | \mathbf{W}, \sigma^{2}) = (2 \pi)^{-\frac{DN}{2}} \lvert \mathbf{C} \rvert^{-\frac{N}{2}} e^{-\mathrm{tr}(\mathbf{C}^{-1} \mathbf{Y}^{T} \mathbf{Y}) }로 주어진다. 이 때의 최대가능도추정은 \mathbf{Y}^{T} \mathbf{Y}의 고유벡터에 대해 나타낼 수 있다.

이의 거울상 문제로 \mathbf{Z}를 최대화시키고 \mathbf{W}를 적분해 뺴내는 문제를 생각해 보자. 사전분포 p(\mathbf{W}) = \prod_{j} \mathcal{N} (\mathbf{w}_{j} | \mathbf{0}, \mathbf{I})를 사용하면 해당하는 가능도는 \mathbf{K}_{z} = \mathbf{Z} \mathbf{Z}^{T} + \sigma^{2} \mathbf{I}일 때 p(\mathbf{Y} | \mathbf{Z}, \sigma^{2}) = \prod_{d=1}^{D} \mathcal{N} (\mathbf{y}_{:, d} | \mathbf{0}, \mathbf{Z} \mathbf{Z}^{T} + \sigma^{2} \mathbf{I}) = (2 \pi)^{-\frac{DN}{2}} \lvert \mathbf{K}_{z} \rvert^{-\frac{D}{2}} e^{-\frac{1}{2} \mathrm{tr}(\mathbf{K}_{z}^{-1} \mathbf{Y} \mathbf{Y}^{T})}이 된다. 이도 고유값 법으로 풀 수 있다.

선형 커널을 사용한다면 주성분 분석을 얻는다. 하지만 보다 일반적인 커널을 사용할 수도 있다: \mathbf{K}_{z} = \mathbf{K} + \sigma^{2} \mathbf{I}를 사용하는 것이다. 이 때 최대가능도추정 \hat{\mathbf{Z}}은 더 이상 고유값 방법으로 직접 구할 수는 없지만 경사 하강법 기반 최적화를 사용할 수는 있다. 목적 함수는 다음과 같다.

l = -\frac{D}{2} \log \lvert \mathbf{K}_{z} \rvert - \frac{1}{2} \mathrm{tr}(\mathbf{K}_{z}^{-1} \mathbf{Y} \mathbf{Y}^{T})

그라디언트는 다음과 같다.

\frac{\partial l}{\partial Z_{ij}} = \frac{\partial l}{\partial \mathbf{K}_{z}} \frac{\partial \mathbf{K}_{z}}{\partial Z_{ij}}, \frac{\partial l}{\partial \mathbf{K}_{z}} = \mathbf{K}_{z}^{-1} \mathbf{Y} \mathbf{Y}^{T} \mathbf{K}_{z}^{-1} - D \mathbf{K}_{z}^{-1}

두번째 항 \frac{\partial \mathbf{K}_{z}}{\partial Z_{ij}}은 당연하지만 사용하는 커널에 따라 달라진다. 이후 일반적인 경사 하강 최적화를 사용하면 된다.

가우시안 과정 잠재 변수 모델과 커널 주성분 분석을 비교하면 어떻게 되는가? 커널 주성분 분석에서는 관측 데이터 공간으로부터 잠재 변수 공간으로의 커널화된 함수를 학습한다. 반면에 가우시안 과정 잠재 변수 모델에서는 잠재 변수 공간으로부터 관측 데이터 공간으로의 커널화된 함수를 학습한다. 가우시안 과정 잠재 변수 모델은 여러 확률적 생성 모델의 이점을 그대로 갖는데, 이 이점들에는 유실된 데이터와 여러 다른 형태의 데이터를 다룰 수 있는 것, 커널 매개변수를 튜닝할 때 격자점 탐색 대신 경사 하강법을 쓸 수 있는 것, 사전분포를 쓸 수 있는 것 등이 있다.

15.6. Approximation methods for large datasets

가우시안 과정의 주된 단점은 사용하는 데 O(N^{3}) 시간이 걸린다는 점이다. 이는 커널 행렬 \mathbf{K}의 역행렬을 구해야 하기 때문이다. 사용자 정의 매개변수 M에 대해 O(M^{2} N)의 시간을 갖는 근사법이 존재한다.

요점 정리

  • 가우시안 과정 : 밀도함수에 대한 베이지안 추정. 함수에 대한 사전분포를 정의한 뒤 관측된 데이터를 기반으로 이것을 함수에 대한 사후분포로 변환시킨다.
  • 결합 가우시안 분포의 성질을 통해 회귀 문제에서 가우시안 과정을 쓸 수 있음.
  • 일반화된 선형 모델과 결합시켜 분류 문제에서 가우시안 과정을 쓸 수 있음.
  • 가우시안 과정은 기계학습에 쓰이는 여러 다른 방법과 밀접한 연관이 있음.
  • 가우시안 과정은 잠재 변수 모델과도 결합될 수 있음.
  • 가우시안 과정은 O(N^{3})의 시간복잡도를 가지지만, 더 짧은 시간이 걸리도록 근사할 수 있음.

답글 남기기

아래 항목을 채우거나 오른쪽 아이콘 중 하나를 클릭하여 로그 인 하세요:

WordPress.com 로고

WordPress.com의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Google photo

Google의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Twitter 사진

Twitter의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

Facebook 사진

Facebook의 계정을 사용하여 댓글을 남깁니다. 로그아웃 /  변경 )

%s에 연결하는 중