본 글은 2014년에 발표된 생성 모델인 Variational AutoEncoder에 대해 설명하고 이를 코드로 구현하는 내용을 담고 있다. VAE에 대해서 알기 위해서는 Variational Inference (변분 추론)에 대한 사전지식이 필요하다. 이에 대해 알고 싶다면 이 글을 참조하길 바란다.

본 글은 크게 3가지 파트로 구성되어 있다. Chapter1에서는 VAE 논문을 리뷰할 것이다. Chapter2에서는 먼저 논문을 간단히 요약하고, 논문의 부록에 수록되어 있었던 미분 가능한 KL-Divergence에 대한 예시를 소개할 것이다. (요약본에 대해서 먼저 보고 싶다면 2.1을 먼저 보라.) Chapter3에서는 Tensorflow를 통해 VAE를 구현할 것이다.

1. Auto-Encoding Variational Bayes 논문 리뷰

1.1. Introduction

연속형 잠재 변수와 파라미터가 다루기 힘든 사후 분포를 갖는 방향성 확률 모델에 대해 효율적인 근사 추론 및 학습을 수행할 수 있는 방법이 없을까? Variational Bayesian 접근법은 다루기 힘든 사후 분포에 대한 근사의 최적화를 내포한다.

불행히도, 일반적인 Mean-Field 접근법은 근사적 사후 분포에 대해 기댓값의 analytic한 해결법을 요구하는데 이는 보통 굉장히 intractable한 방법이다. 본 논문은 Variational Lower Bound(ELBO)의 Reparameterization이 Lower Bound의 미분 가능한 불편향 estimator를 만드는 방법에 대해 보여줄 것이다. 이 Stochastic Gradient Variational Bayes: SGVB estimator는 연속형 잠재변수나 파라미터를 갖고 있는 대부분의 모델에 대해 효율적인 근사 사후 추론을 가능하게 하며, 표준 Stochastic Gradient Ascent 스킬을 사용하여 최적화하기에 굉장히 편리하다.

iid 데이터셋이고, 데이터 포인트 별로 연속형 잠재변수를 갖고 있는 경우에 대해 본 논문은 Auto-Encoding VB 알고리즘을 제안한다. 이 알고리즘에서는 Simple Ancestral Sampling을 이용하여 근사 사후 추론을 하는 인식 모델을 최적화하기 위해 SGVB estimator를 사용하여 추론과 학습을 효율적으로 해낸다. 이 과정은 MCMC와 같이 데이터포인트 별로 반복적인 추론을 행하여 많은 연산량을 요구하지 않는 장점을 가진다.

학습된 근사 사후 추론 모델은 recognition, denoising, representation, visualization의 목적으로 활용될 수 있다. 본 알고리즘이 인식(recognition) 모델에 사용될 때, 이를 Variational Auto-Encoder라고 부를 것이다.

1.2. Method

본 섹션에서는 연속형 잠재 변수를 내포하는 다양한 방향성 그래픽 모델(Directed Graphical Model)에서 Stochastic 목적 함수인 Lower Bound Estimator를 끌어내는 과정을 설명할 것이다. 데이터포인트 별 잠재변수는 iid한 상황이라는 가정 하에 본 논문에서는 (전역) 파라미터에 대해 Maximum Likelihood와 Maximum A Posteriori 추론을 수행하고 잠재변수에 대해 Variational Inference를 수행할 것이다. 이러한 방법은 온라인 러닝에도 사용될 수 있지만 본 논문에서는 간단히 하기 위해 고정된 데이터셋을 사용할 것이다.

1.2.1. Problem Scenario

N개의 Sample을 가진 $X$ 라는 데이터가 있다고 해보자. 본 논문은 이 데이터가 관측되지 않은 연속형 확률 변수 $\mathbf{z}$ 를 내포하는 어떤 Random Process에 의해 형성되었다고 가정한다.

이 과정은 2가지 단계로 구성된다.

[\mathbf{z}^{(i)} \sim Prior: p_{\theta^*}(\mathbf{z})]

(여기서 $\mathbf{z}$는 원인, $\mathbf{x}$는 결과라고 보면 이해가 쉬울 것이다.)

이 때 우리는 위 2개의 확률이 모두 아래 두 분포의 Parametric Families of Distributions에서 왔다고 가정한다.

이들의 확률밀도함수는 거의 모든 $\theta, \mathbf{z}$에 대해 미분가능하다고 전제한다.

불행히도, 이러한 과정의 많은 부분은 우리가 직접 확인하기 어렵다. True 파라미터인 $\theta^*$ 와 잠재 변수의 값 $\mathbf{z}^{(i)}$ 는 우리에게 알려져 있지 않다.

본 논문은 주변 확률이나 사후 확률에 대한 단순화를 위한 일반적인 가정을 취하지 않고, 아래에서 제시한 상황처럼 분포가 intractable하고 큰 데이터셋을 마주하였을 경우를 위한 효율적인 알고리즘에 대해 이야기하고자 한다.

1) Intractability

(1) Marginal Likelihood $p_{\theta}(\mathbf{x})$ 의 적분은 위 식으로 표현되는데, 이 식이 intractable한 경우가 존재한다. 이 경우는 Evidence가 적분이 불가능한 경우를 의미한다.

(2) True Posterior Density가 intractable한 경우 (EM알고리즘이 사용될 수 없음)

True Posterior Density는 아래와 같다.

(3) 어떠한 합리적인 Mean-Field VB 알고리즘을 위한 적분이 불가능한 경우

이러한 Intractability는 굉장히 흔하며, 복잡한 Likelihood 함수를 갖는 신경망 네트워크에서 발견할 수 있다.

2) A Large Dataset

데이터가 너무 크면 배치 최적화는 연산량이 매우 많다. 우리는 작은 미니배치나 데이터포인트에 대해 파라미터 업데이트를 진행하고 싶은데, Monte Carlo EM과 같은 Sampling Based Solution은 데이터 포인트별로 Sampling Loop를 돌기 때문에 너무 느리다.

위 시나리오에서 설명한 문제들에 대해 본 논문은 아래와 같은 해결책을 제시한다.

첫 번째로, 파라미터 $\theta$ 에 대한 효율적인 근사 ML/MAP 추정을 제안한다. 이 파라미터들은 숨겨진 Random Process를 흉내내고 실제 데이터를 닮은 인공적인 데이터를 생성할 수 있게 해준다.

두 번째로, 파라미터 $\theta$ 의 선택에 따라 관측값 $\mathbf{x}$ 가 주어졌을 때 잠재 변수 $\mathbf{z}$ 에 대한 효율적인 근사 사후 추론을 제안한다.

세 번째로, 변수 $\mathbf{x}$ 에 대한 효율적인 근사 Marginal Inference를 제안한다. 이는 $\mathbf{x}$ 에 대한 prior이 필요한 모든 추론 task를 수행할 수 있게 해준다.

위 문제를 해결하기 위해 아래와 같은 인식 모델이 필요하다.

이 모델은 intractable한 True Posterior의 근사 버전이다.

기억해야 할 것이, Mean-Field Variational Inference에서의 근사 Posterior와는 다르게 위 인식 모델은 꼭 계승적(factorial)일 필요도 없고, 파라미터 $\phi$ 가 닫힌 형식의 기댓값으로 계산될 필요도 없다.

본 논문에서는 인식 모델 파라미터인 $\phi$ 와 생성 모델 파라미터인 $\theta$를 동시에 학습하는 방법에 대해 이야기할 것이다.

코딩 이론의 관점에서 보면, 관측되지 않은 변수 $\mathbf{z}$ 는 잠재 표현 또는 code라고 해석될 수 있다. 본 논문에서는 인식 모델을 Encoder라고 부를 것이다. 왜냐하면 데이터 포인트 $\mathbf{x}$ 가 주어졌을 때 이 Encoder가 데이터 포인트 $\mathbf{x}$ 가 발생할 수 code $\mathbf{z}$의 가능한 값에 대한 분포를 생산하기 때문이다.

비슷한 맥락에서 우리는 생성 모델을 확률적 Decoder라고 명명할 것인데, 왜냐하면 code $\mathbf{z}$ 가 주어졌을 때 이 Decoder가 상응하는 가능한 $\mathbf{x}$ 의 값에 대해 분포를 생성하기 때문이다.

1.2.2. The Variational Bound

Marginal Likelihood는 각 데이터 포인트의 Marginal Likelihood의 합으로 구성된다. 이를 식으로 표현하면 아래와 같다.

[log p_{\theta}(\mathbf{x}^{(1)}, …, \mathbf{x}^{(N)}) = \sum_{i=1}^N log p_{\theta} (\mathbf{x}^{(i)})]

그런데 데이터 포인트 하나에 대한 Marginal Likelihood는 아래와 같이 재표현이 가능하다.

우변의 첫 번째 항은 True Posterior의 근사 KL Divergence이다. 이 KL Divergence는 음수가 아니기 때문에, 두 번째 항인 $\mathcal{L}(\theta, \phi; \mathbf{x}^{(i)})$ 는 i번째 데이터 포인트의 Marginal Likelihood의 Varitaional Lower Bound라고 한다. 변분 추론 글을 읽어보고 왔다면 알겠지만, 이는 Evidence의 하한 값을 뜻하기도 하기 때문에 ELBO라고 부르기도 한다.

부등식으로 나타내면 아래와 같다.

이 식은 또 아래와 같이 표현할 수 있다.

우리는 Lower Bound $L(\theta, \phi; \mathbf{x}^{(i)})$ 를 Variational 파라미터와 생성 파라미터인 $\phi, \theta$ 에 대하여 미분하고 최적화하고 싶은데, 이 Lower Bound의 $\phi$ 에 대한 Gradient는 다소 복잡하다.

이러한 문제에 대한 일반적인 Monte Carlo Gradient Estimator는 아래와 같다.

그런데 이 Gradient Estimator는 굉장히 큰 분산을 갖고 있어서 우리의 목적에 적합하지 않다.

1.2.3. The SGVB estimator and AEVB algorithm

이 섹션에서는 Lower Bound의 실용적인 추정량과 파라미터에 대한 그 추정량의 미분 값을 소개할 것이다.

일단 우리는 위와 같은 근사 Posterior를 가정한다. 다만 x에 대한 조건부를 가정하지 않은 케이스인 $q_{\phi}(\mathbf{z})$ 와 같은 경우에도 같은 테크닉을 적용할 수 있음에 유의하자. Posterior를 추론하는 Fully Variational Bayesian 방법론은 본 논문의 부록에 소개되어 있다.

위에서 설정한 근사 Posterior에 대해 1.2.4 섹션에서 정한 certain mild conditions 하에 우리는 그 근사 Posterior를 따르는 확률 변수 $\tilde{\mathbf{z}}$ 를 Reparameterize 할 수 있는데, 이는 Auxiliary Noise 변수 $\epsilon$ 의 Diffrentiable Transformation($g_{\phi} (\epsilon, \mathbf{x})$)를 통해 이루어진다.

[\tilde{\mathbf{z}} \sim g_{\phi} (\epsilon, \mathbf{x}), \epsilon \sim p(\epsilon)]

다음 Chapter를 보면 적절한 $p(\epsilon)$ 분포와 $g_{\phi} (\epsilon, \mathbf{x})$ 함수를 선정하는 일반적인 방법에 대해 알 수 있다. 이제 이 Chapter 서두에서 언급한 근사 Posterior에 대해 어떤 함수 $f(\mathbf{z})$ 기댓값의 Monte Carlo 추정량을 다음과 같이 쓸 수 있다.

이는 실제로 계산하기 어려운 어떤 함수에 대해 $L$ 개의 Sample를 뽑아 이에 대한 근사치로 추정량을 구하는 방법이다.

이제 이 테크닉을 ELBO에 대해 적용하면 SGVB: Stochastic Gradient Variational Bayes 추정량을 얻을 수 있다.

[\tilde{\mathcal{L}}^A (\theta, \phi ; \mathbf{x}^{(i)}) \simeq \mathcal{L} (\theta, \phi ; \mathbf{x}^{(i)})]

[g_{\phi} (\epsilon^{(i, l)}, \mathbf{x}^{(i)}), \epsilon^{(l)} \sim p(\epsilon)]

잠시 아래 식에서 쿨백-라이블리 발산에 주목해보자.

이 쿨백-라이블리 발산 값은 종종 analytic 하게 적분될 수 있는데 이렇게 되면 오직 Expected Reconstruction Error만이 샘플링에 의한 추정을 필요로하게 된다. Expected Reconstruction Error 항은 아래 식을 의미한다.

쿨백-라이블리 발산 항은 근사 Posterior를 Prior $p_{\theta}(z)$ 에 가깝게 만들어서 $\phi$ 를 규제하는 것으로 해석될 수 있다.

이러한 과정을 SGVB 추정량의 두 번째 버전으로 이어지는데, 이 추정량은 일반적인 추정량에 비해 작은 분산을 갖고 있다.

이 때

[\mathbf{z}^{(i, l)} = g_{\phi} (\epsilon^{(i, l)}, \mathbf{x}^{(i)}), \epsilon^{(l)} \sim p(\epsilon)]

$N$ 개의 데이터 포인트를 갖고 있는 데이터셋 $X$ 에서 복수의 데이터 포인트가 주어졌을 때 우리는 미니배치에 기반하여 전체 데이터셋에 대한 Marginal Likelihood Lower Bound의 추정량을 구성할 수 있다.

[\mathcal{L} (\theta, \phi ; X) \simeq \tilde{\mathcal{L}}^M (\theta, \phi ; X^M) = \frac{N}{M} \Sigma_{i=1}^M \tilde{\mathcal{L}} (\theta, \phi ; \mathbf{x}^{(i)})]

이 때

[X^M = [\mathbf{x}^{i}]_{i=1}^M]

미니배치 $X^M$ 은 전체 데이터셋 $X$ 에서 랜덤하게 뽑힌 샘플들을 의미한다. 본 논문의 실험에서, 미니 배치 사이즈 $M$ 이 (예를 들어 100) 충분히 크면 데이터 포인트 별 sample의 크기인 $L$ 이 1로 설정될 수 있다는 사실을 알아 냈다.

[\triangledown_{\theta, \phi} \tilde{\mathcal{L}} (\theta; X^M)]

위와 같은 derivate가 도출될 수 있고, 이에 따른 Gradients는 SGD나 Adagrad와 같은 확률적 최적화 방법과 연결되어 사용될 수 있다.

다음은 기본적인 Stochastic Gradients를 계산하는 방법에 대한 내용이다.

아래 목적 함수를 보면 Auto-Encoder와의 연결성이 더욱 뚜렷해진다.

Prior로부터 나온 근사 Posterior에 대한 쿨백-라이블리 발산 값인 첫 번째 항은 Regularizer의 역할을 하며, 두 번째 항은 Expected Negative Reconstruction Error의 역할을 하게 된다.

$g_{\phi} (.)$ 라는 함수는 데이터 포인트 $\mathbf{x}^{(i)}$ 와 Random Noise Vector $\epsilon^{(l)}$ 을 데이터 포인트 $\mathbf{z}^{(i, l)}$ 을 위한 근사 Posterior로 부터 추출된 Sample로 매핑하는 역할을 수행한다.

그 후, 이 Sample $\mathbf{z}^{(i, l)}$ 은 아래 함수의 Input이 된다.

이 함수는 $\mathbf{z}^{(i, l)}$ 이 주어졌을 때 생성 모델 하에서 데이터 포인트 $\mathbf{x}^{(i)}$ 의 확률 밀도 함수와 동일하다. 이 항은 Auto-Encoder 용어에서 Negative Reconstruction Error에 해당한다.

1.2.4. The Reparamaterization Trick

위에 제시된 근사 Posterior로부터 Sample을 생성하기 위해 본 논문에서는 다른 방법을 적용하였다. 본질적인 Parameterization 트릭은 굉장히 단순하다. $\mathbf{z}$ 가 연속형 확률 변수이고 아래와 같이 어떤 조건부 분포를 따른다고 하자.

그렇다면 이제 이 확률 변수 $\mathbf{z}$ 를 다음과 같은 Deterministic 변수라고 표현할 수 있다.

[\mathbf{z} = g_{\phi} (\epsilon, \mathbf{x})]

이 때 $\epsilon$ 은 독립적인 Marginal $p(\epsilon)$ 을 가지는 보조 변수이고, $g_{\phi} (.)$ 함수는 $\phi$ 에 의해 parameterized 되는 vector-valued 함수이다.

이 Reparameterization은 매우 유용하다. 왜냐하면 근사 Posterior의 기댓값을 $\phi$에 대해 미분 가능한 기댓값의 Monte Carlo 추정량으로 재표현하는 데에 사용될 수 있기 때문이다.

증명은 다음과 같다. $z = g_{\phi} (\epsilon, \mathbf{x})$ 라는 Deterministic Mapping이 주어졌을 때 우리는 다음 사실을 알 수 있다.

이는 각 Density에 derivatives를 곱한 것이 서로 같다는 것을 의미한다. 참고로 위에서 $dz_i$ 라고 표기한 것은, 무한소에 대한 수식이기 때문이며 사실 $d\mathbf{z} = \prod_i dz_i$ 이다. 따라서 아래와 같은 식을 구성할 수 있다.

마지막 식을 보면 $\mathbf{z}$ 를 Deterministic 하게 표현하여 오로지 $\epsilon, \mathbf{x}$ 로만 식을 구성한 것을 알 수 있다. 이제 미분 가능한 추정량을 구성할 수 있다.

[f(g_{\phi} (\epsilon, \mathbf{x}) ) d\epsilon \simeq \frac{1}{L} \Sigma_{l=1}^L f( g_{\phi} (\mathbf{x}, \epsilon^{(l)}) )]

사실 이는 이전 Chapter에서 ELBO의 미분 가능한 추정량을 얻기 위해 적용하였던 트릭이다.

예를 들어 단변량 정규분포의 케이스를 생각해보자.

이 경우 적절한 Reparameterization은 $z = \mu + \sigma \epsilon$이다. 이 때 $\epsilon$은 $\mathcal{N} (0, 1)$을 따르는 보조 Noise 변수이다. 그러므로 우리는 다음 식을 얻을 수 있다.

[E_{\mathcal{N} (z; \mu, \sigma^2)} [f(z)] = E_{\mathcal{N} (\epsilon; 0, 1)} [f(\mu + \sigma \epsilon)] \simeq \frac{1}{L} \Sigma^L_{l=1} f(\mu + \sigma \epsilon^{(l)}), \epsilon^{(l)} \sim \mathcal{N} (0, 1)]

그렇다면 어떤 근사 Posterior에 따라 미분 가능한 변환 함수 $g_{\phi}(.)$ 와 보조 변수 $\epsilon \sim p(\epsilon)$ 은 어떻게 선택하는가? 기본적인 방법은 다음과 같다.

첫 번째로, Tractable Inverse CDF를 선택한다. $\epsilon \sim \mathcal{U} (\mathbf{0}, \mathbf{I})$ 일 때, $g_{\phi} (\epsilon, \mathbf{x})$ 는 근사 Posterior의 Inverse CDF라고 해보자. 예를 들면, Exponential, Cauchy, Logistic, Rayleigh, Pareto, Weibull, Reciprocal, Gompertz, Gumbel, Erlang 분포가 있다.

두 번째로, 정규 분포 예시와 유사하게 어떠한 location-scale family of distributions에 대해 우리는 location=0, scale=1인 표준 분포를 보조 변수 $\epsilon$ 으로 선택할 수 있다.

그렇다면 $g(.)$ = location + scale * $\epsilon$ 이 될 것이다. 예를 들면, Laplace, Elliptical, T, Logistic, Uniform, Triangular, Gaussian이 있다.

마지막으로 분포를 결합하는 방식을 채택할 수 있다. 종종 확률 변수를 보조 변수의 다른 변환으로 표현할 수 있다. 예를 들어 Log-Normal, Gamma, Dirichlet, Beta, Chi-Squared, F가 있다.

이 3가지 방법이 모두 통하지 않는다면, Inverse CDF의 좋은 근사법은 PDF에 비해 많은 연산량과 소요 시간을 요한다.

1.3. Example: Variational Auto-Encoder

이번 Chapter에서는 확률론적 Encoder (생성 모델의 Posterior의 근사) 와 AEVB 알고리즘을 통해 파라미터 $\phi$ 와 $\theta$ 가 Jointly 최적화되는 신경망에 대한 예시를 다루도록 할 것이다.

잠재 변수에 대한 Prior는 Centered Isotropic Multivariate Gaussian $p_{\theta} (\mathbf{z}) = \mathcal{N} (\mathbf{z}; \mathbf{0}, \mathbf{I})$ 라고 하자. 이 경우에서는 Prior에 파라미터가 없다는 사실을 염두에 두자.

위 분포는 다변량 정규 분포 혹은 베르누이 분포(각각 실수 데이터, 또는 Binary 데이터 일때)로 설정하고, 분포의 파라미터는 Multi Layer Perceptron에서 $\mathbf{z}$ 로 계산된다.

이 경우 위 식의 왼쪽 부분인 True Posterior는 intractable하다. 위 식의 오른쪽 부분인 근사 Posterior에 더 많은 자유가 있기 때문에, 우리는 True(but intractable) Posterior가 Approximately Diagonal Covariance를 가진 근사 정규분포를 따른다고 가정한다. 이 경우 우리는 Variational Approximate Posterior가 Diagonal Covariance 구조를 가진 다변량 정규분포라고 설정할 수 있다.

이 때 근사 Posterior의 평균과 표준편차는 MLP의 Output이다. 이 때 MLP는 Variational Parameter인 $\phi$ 와 데이터 포인트 $\mathbf{x}^{(i)}$ 의 비선형적 함수라고 생각할 수 있다.

1.2.4 에서 설명하였듯이, 우리는 다음과 같이 근사 Posterior에서 Sampling을 한다.

이 때 아래 사실을 이용한다.

[\mathbf{z}^{(i, l)} = g_{\phi} (\mathbf{x}^{(i)} , {\epsilon}^{(l)} ) = {\mu}^{(i)} + {\sigma}^{(i)} \odot {\epsilon}^{(l)}, {\epsilon}^{(l)} \sim \mathcal{N} (\mathbf{0}, \mathbf{I})]

이 모델에서 Prior와 근사 Posterior는 정규 분포를 따른다. 그러므로 우리는 이전 Chapter에서 보았던 추정량을 사용할 수 있는데, 이 때 쿨백-라이블리 발산 값은 추정 없이 계산되고 미분될 수 있다.

이 때

[\mathbf{z}^{(i, l)} = g_{\phi} (\mathbf{x}^{(i)} , {\epsilon}^{(l)} ) = {\mu}^{(i)} + {\sigma}^{(i)} \odot {\epsilon}^{(l)}, {\epsilon}^{(l)} \sim \mathcal{N} (\mathbf{0}, \mathbf{I})]

[log p_{\theta} (\mathbf{x}^{(i)}, \mathbf{z}^{(i, l)})]

바로 위에 있는 Decoding Term은 우리가 모델링하는 데이터의 종류에 따라 베르누이 또는 정규분포 MLP가 된다.

1.4. Related Work

Wake-Sleep 알고리즘은 연속적인 잠재 변수를 갖고 있는 같은 문제를 갖고 있는 상황에 적용할 수 있는 다른 방법론이다. 본 논문에서 제시된 방법과 마찬가지로 이 알고리즘은 True Posterior를 근사하기 위해 인식 모델을 사용한다. 그러나 이 이 알고리즘의 경우 동시에 발생하는 2개의 목적함수를 최적화해야 하기 때문에 Marginal Likelihood의 최적화에 적합하지 않다. 이 알고리즘은 이산적인 잠재변수를 다루는 모델에 적합하다. Wake-Sleep 알고리즘의 데이터 포인트 별 계산 복잡도는 AEVB와 같다.

Stochastic Variational Inference는 주목 받는 알고리즘이다. 이 논문에서는 2.1 장에서 Naive Gradient Estimator의 높은 분산을 감소시키기 위해 Control Variate Scheme을 소개하였고, 이를 Posterior의 Exponential Family Approximation에 적용하였다.

AEVB 알고리즘은 Variational Objective로 학습되는 Directed Probabilistic Model과 Auto-Encoder 사이의 연결성을 밝혀준다. 선형적인 Auto-Encoder와 Generative Linear-Gaussian 모델의 특정 종류 사이의 연결성은 오래 전부터 알려져왔다. 이 논문은 PCA가 아래 조건에 해당하는 Linear-Gaussian Model의 Maximum Likelihood라는 점을 밝히고 있다. (특히 매우 작은 $\epsilon$ 일 때)

Auto-Encoder와 관련된 연구 중에서 이 논문은 unregularized auto-encoder는 input $X$ 와 잠재 표현 $Z$ 의 Mutual Information의 Maximum Lower Bound라는 것을 보였다. Mutual Information을 최대화하는 것은 결국 조건부 엔트로피를 최대화하는 것과 같은데, 이 조건부 엔트로피는 autoencoding 모델 하에서 데이터의 expected loglikelihood에 의해 lower bounded된다. (Negative Reconstruction Error)

그러나 이 Reconstruction 기준은 유용한 표현을 학습하기에는 sufficient하지 않다. denoising, contractive, sparse autoencoder 변인들이 Auto-Encoder로 하여금 더 유용한 표현을 학습하도록 제안되었다.

SGVB 목표함수는 Variational Bound에 의해 좌우되는 규제화 항을 포함하며, 일반적인 불필요한 규제화 하이퍼 파라미터들은 포함하지 않는다.

Predictive Sparse Decomposition과 같은 Encoder-Decoder 구조 역시 본 논문에서 영감을 받은 방법론이다. Generative Stochastic Network를 발표한 이 논문에서는 Noisy autoencoder가 데이터 분포로부터 Sampling을 하는 Markov Chain의 Transition Operator를 학습한다는 내용이 소개되어 있다. 다른 논문에서는 Deep Boltzmann 머신과 함께 인식 모델이 사용되기도 하였다. 이러한 방법들은 비정규화된 모델들에 타겟팅되어 있거나 희소 코딩 모델에 있어 제한적이며, 본 논문에서 제안된 알고리즘은 Directed Probabilisitic Model의 일반적인 경우를 학습할 수 있기 때문에 차별화 된다.

1.5. Experiments

(논문 참조)

1.6. Conclusion

본 논문에서 Variational Lower Bound의 새로운 추정량인 SGVB를 새롭게 소개하였다. 이 알고리즘은 연속적인 잠재 변수들에 대한 효율적인 근사적 추론을 가능하게 해준다. 제안된 추정량은 직접적으로 미분이 가능하고 표준적인 Stochastic Gradient 방법들을 사용하여 최적화될 수 있다.

iid 데이터셋과 연속적인 잠재 변수에 한해 본 논문에서는 효과적인 추론과 학습 방법인 AEVB: Auto-Encoding VB = VAE를 제안하였는데, 이 알고리즘은 SGVB 추정량을 사용하여 근사적인 추론을 행하는 모델이다. 이론적인 장점은 실험 결과에 반영되어 있다.

2. 보충 설명

2.1. VAE 요약

우리가 풀고 싶은 문제는 이것이다. 연속적인 잠재 변수가 존재한다고 할 때 데이터에 기반하여 이에 대한 효과적인 학습과 추론을 행하고 싶다. 그런데 문제가 존재한다.

위와 같은 Evidence, Posterior, Likelihood가 intractable한 경우가 매우 흔하게 존재한다. 그리고 데이터셋이 크다면 이들에 대해 배치 최적화나 Monte-Carlo Sampling을 통한 추론을 행하기에는 시간이 너무 오래 걸린다.

이를 위해 아래와 같은 새로운 인식 모델이 제안된다. 이 모델은 True Posterior를 근사하기 위해 제안되었으며 이러한 방법론을 Variational Inference 혹은 Variational Bayes라고 부른다.

기본적으로 ELBO를 최대화하는 $\theta, \phi$ 를 찾는 것으로 근사가 이루어진다.

[\theta^, \phi^ = argmax_{\theta, \phi} \Sigma_{i=1}^{N=1} \mathcal{L} (\theta, \phi ; \mathbf{x}^{(i)})]

최적화된 파라미터를 찾는 analytic한 전통적인 방법은 Coordinate Ascent Mean-Field Variational Inference이다. 글 서두에서도 밝혔듯이 이에 대한 내용은 이 글에서 확인할 수 있다. 이 방법은 여러 단점이 있는데, 그 중 하나는 factorial한 표현을 전제로 하기 때문에 본 논문에서와 같이 intractable한 Likelihood를 갖는 경우에는 사용이 불가능하다. 따라서 본 논문에서는 SGVB 추정량을 제안하고 있다.

SGVB는 파라미터의 Gradient를 구해 Stochastic하게 업데이트하는 방식을 취하는데, ELBO에서 $\phi$ 의 Gradient를 얻는 것은 다소 복잡하다. 따라서 Monte Carlo 추정량을 얻는 것을 생각할 수 있는데, 이 또한 분산이 너무 커서 직접적으로 사용하기에는 무리가 있다.

그래서 최종적으로 채택한 방법이 근사 Posterior를 따르는 확률 변수 $\tilde{\mathbf{z}}$ 에 대해 Reparameterization을 행하는 것이다. 근사 Posterior에서 직접 $\mathbf{z}$ 를 Sampling 하는 것이 아니라, 보조 Noise 변수 $\epsilon$ 을 사용하여 미분 가능한 분포에서 Deterministic하게 정해지는 것으로 파악하는 것이다.

[\tilde{\mathbf{z}} = g_{\phi} (\epsilon, x), \epsilon \sim p(\epsilon)]

기존의 ELBO는 아래와 같다.

지금부터 2가지 케이스가 존재한다. 만약 첫 번째 항인 쿨백-라이블리 발산이 analytic하게 적분이 되지 않는다면 위 식 전체를 Monte-Carlo 추정을 통해 구해야 한다. 만약 가능하다면, 오직 두 번째 항만을 Monte-Carlo 추정을 통해 구하면 된다. (방금 전 분산이 커져서 사용하기 어렵다고 했던 부분은, $\mathbf{z}$ 를 근사 Posterior로부터 직접 Sampling을 할 때의 이야기이다.)

첫 번째 케이스(A)를 살펴보자. ELBO는 아래와 같이 다시 표현할 수 있다.

그대로 Monte-Carlo 추정을 시행하면 아래와 같은 SGVB 추정량을 얻을 수 있다.

이제 두 번째 케이스(B)를 살펴보자. 쿨백-라이블리 발산이 적분이 가능하다고 하였으므로, 두 번째 항인 Expected Reconstruction Error 만이 Sampling에 의한 Monte-Carlo 추정을 필요로 하게 된다. 참고로 첫 번째 항은 Regualizer의 역할을 수행한다.

이 때

[\mathbf{z}^{(i, l)} = g_{\phi} (\epsilon^{(l)}, \mathbf{x}^{(i)}), \epsilon^{(l)} \sim p(\epsilon)]

$\mathbf{z}^{(i, l)}$ 은 위 목적 함수(SGVB-B)의 두 번째 항의 Input이다. 이 데이터 포인트는 $g_{\phi} (\epsilon^{(i)}, \mathbf{x}^{(i)})$ 에서 Sampling 되는 것이며 이 $g_{\phi} (.)$ 라는 함수는 일반적으로 일변량 정규분포로 설정된다. ($g_{\phi} (\epsilon, x) = \mu + \sigma \epsilon$) 이렇게 하면, 위 목적함수를 최적화하고 역전파를 이용하여 학습을 진행할 수 있게 된다.

2.2. Solution of Negative KL-Divergence

이전 Chapter에서 SGVB-B를 구할 때, 쿨백-라이블리 발산 값이 analytic하게 적분될 수 있는 경우를 가정하였다.

Prior와 근사 Posterior 모두 정규분포로 설정해보자.

$J$ 는 $\mathbf{z}$ 의 차원이라고 할 때, Negative 쿨백-라이블리 발산은 아래와 같이 정리할 수 있다.

왜냐하면,

2.3. MLP’s as Probabilistic Encoders and Decoders

데이터의 종류에 따라 Decoder는 Gaussian Output을 반환할 수도, Bernoulli Output을 반환할 수도 있다. 두 경우 논문의 부록에 잘 설명되어 있다.

베르누이 분포일 경우

이 때 $f_{\sigma}(.)$ 함수는 elementwise sigmoid 활성화 함수이다.

정규 분포일 경우

3. Tensorflow로 VAE 구현

Tensorflow 홈페이지에는 (흔히 그렇듯) MNIST 예제로 VAE를 적용하는 방법에 대해 가이드를 제시하고 있다. 코드도 깔끔하고 설명도 어느 정도 되어 있기 때문에 참고하기를 추천하며, 이번 Chapter 역시 그 가이드에 기반하여 작성되었음을 밝힌다. 다만 Tensorflow 홈페이지에서 제시한 예시는 SGVB-A를 사용하고 있는데, 본문을 읽어보면 단순한 예시를 들기 위해 ELBO 전체에 대해 Monte-Carlo 추정을 시행하였다고 밝히고 있다. 이번 Chapter에서는 SGVB-B를 활용하여 Loss Function을 설계할 것이다.

먼저 모델을 정의해보자.

class CVAE(tf.keras.Model):
    def __init__(self, latent_dim):
        super(CVAE, self).__init__()
        self.latent_dim = latent_dim
        self.encoder = tf.keras.Sequential(
            [
                tf.keras.layers.InputLayer(input_shape=(28, 28, 1)),
                tf.keras.layers.Conv2D(filters=32, kernel_size=3, strides=(2, 2), activation='relu'),
                tf.keras.layers.Conv2D(filters=64, kernel_size=3, strides=(2, 2), activation='relu'),
                tf.keras.layers.Flatten(),
                # No activation
                tf.keras.layers.Dense(latent_dim + latent_dim),
            ])

        self.decoder = tf.keras.Sequential(
            [
                tf.keras.layers.InputLayer(input_shape=(latent_dim, )),
                tf.keras.layers.Dense(units=7*7*32, activation=tf.nn.relu),
                tf.keras.layers.Reshape(target_shape=(7, 7, 32)),
                tf.keras.layers.Conv2DTranspose(filters=64, kernel_size=3, strides=2, padding='same', activation='relu'),
                tf.keras.layers.Conv2DTranspose(filters=32, kernel_size=3, strides=2, padding='same', activation='relu'),
                # No activation
                tf.keras.layers.Conv2DTranspose(filters=1, kernel_size=3, strides=1, padding='same'),
            ])

    @tf.function
    def sample(self, eps=None):
        if eps is None:
            eps = tf.random.normal(shape=(100, self.latent_dim))
        return self.decode(eps, apply_sigmoid=True)

    def encode(self, x):
        # encoder의 Output은 (batch_size, latent_dim * 2) 이다. 
        # 각 mini-batch에서 이를 반으로 쪼갠다.
        # logvar: Linear Layer를 통과한 후 음수의 값을 가질 수도 있기 때문에 이와 같이 표기한다.
        mean, logvar = tf.split(self.encoder(x), num_or_size_splits=2, axis=1)
        stddev = 1e-8 + tf.nn.softplus(logvar)
        return mean, stddev

    def reparameterize(self, mean, stddev):
        # 보조 노이즈 변수: eps
        eps = tf.random.normal(shape=mean.shape)
        z = mean + eps * stddev
        return z

    def decode(self, z, apply_sigmoid=False):
        logits = self.decoder(z)
        if apply_sigmoid:
            probs = tf.sigmoid(logits)
            return probs
        return logits

optimizer = tf.keras.optimizers.Adam(1e-4)

모델의 구조는 아래와 같은 그림으로 이해하면 쉬울 것이다.

Prior와 근사 Posterior가 모두 정규 분포라는 가정 하에 Negative KL-Divergence는 아래와 같다.

그리고 2.3을 참고했을 때 다음과 같이 Log Likelihood를 얻을 수 있다.

자 그럼 이제 SGVB-B를 정확히 구해보자. ($J$ 는 잠재 변수의 차원이다.)

[= \frac{1}{2} \Sigma_{j=1}^J (1 + log( (\sigma_j^{(i)})^2 ) - (\mu_j^{(i)})^2 - (\sigma_j^{(i)})^2 ) + \frac{1}{L}\Sigma_{l=1}^L x_i log y_i + (1-x_i) * log(1 - y_i)]

이를 코드로 구현하면 아래와 같다.

def compute_loss(model, x):
    mean, stddev = model.encode(x)
    z = model.reparameterize(mean, stddev)
    x_logit = model.decode(z, True)
    x_logit = tf.clip_by_value(x_logit, 1e-8, 1-1e-8)

    # Loss
    marginal_likelihood = tf.reduce_sum(x * tf.math.log(x_logit) + (1 - x) * tf.math.log(1 - x_logit), axis=[1, 2])
    loglikelihood = tf.reduce_mean(marginal_likelihood)

    kl_divergence = -0.5 * tf.reduce_sum(1 + tf.math.log(1e-8 + tf.square(stddev)) - tf.square(mean) - tf.square(stddev),
                                         axis=[1])
    kl_divergence = tf.reduce_mean(kl_divergence)

    ELBO = loglikelihood - kl_divergence
    loss = -ELBO

    return loss

학습 및 테스트 코드는 아래와 같다.

@tf.function
def train_step(model, x, optimizer):
    with tf.GradientTape() as tape:
        loss = compute_loss(model, x)
    gradients = tape.gradient(loss, model.trainable_variables)
    optimizer.apply_gradients(zip(gradients, model.trainable_variables))

    return loss


epochs = 50
latent_dim = 2
model = CVAE(latent_dim)

# Train
for epoch in range(1, epochs + 1):
    train_losses = []
    for train_x in train_dataset:
        loss = train_step(model, train_x, optimizer)
        train_losses.append(loss)

    print('Epoch: {}, Loss: {:.2f}'.format(epoch, np.mean(train_losses)))

    # metric = tf.keras.metrics.Mean()
    # for test_x in test_dataset:
    #     metric(compute_loss(model, test_x))
    # elbo = -metric.result()


# Test
def generate_images(model, test_sample, random_sample=False):
    mean, stddev = model.encode(test_sample)
    z = model.reparameterize(mean, stddev)

    if random_sample:
        predictions = model.sample(z)
    else:
        predictions = model.decode(z, True)

    predictions = tf.clip_by_value(predictions, 1e-8, 1 - 1e-8)

    fig = plt.figure(figsize=(4, 4))

    for i in range(predictions.shape[0]):
        plt.subplot(4, 4, i + 1)
        plt.imshow(predictions[i, :, :, 0], cmap='gray')
        plt.axis('off')

    plt.show()


num_examples_to_generate = 16
random_vector_for_generation = tf.random.normal(shape=[num_examples_to_generate, latent_dim])
test_sample = next(iter(test_dataset))[0:num_examples_to_generate, :, :, :]

for i in range(test_sample.shape[0]):
    plt.subplot(4, 4, i + 1)
    plt.imshow(test_sample[i, :, :, 0], cmap='gray')
    plt.axis('off')

plt.show()

generate_images(model, test_sample, False)

실제로 학습을 시켜보면, Loss가 155까지는 빠르게 떨어지고, 그 이후에는 아주 서서히 감소하는 것을 알 수 있다. 이렇게 구현한 Convolutional VAE의 경우 성능은 좀 아쉽다. 아주 미세하게 숫자를 구분해 내지는 못한다. 아래의 이미지는 Epoch 30 이후의 결과이다.

원본

생성본

References

1) https://arxiv.org/abs/1312.6114
2) https://ratsgo.github.io/generative%20model/2018/01/27/VAE/
3) https://www.youtube.com/watch?v=SAfJz_uzaa8
4) https://taeu.github.io/paper/deeplearning-paper-vae/
5) https://dnddnjs.github.io/paper/2018/06/20/vae2/
6) https://www.tensorflow.org/tutorials/generative/cvae

이번 글에서는 Monte Carlo 근사 방법에 대해 간단히 알아볼 것이다. 사실 개념 자체는 간단하니 가볍게 파악해보록 하자. 본 글은 Kevin P.Murphy의 Machine Learning: A Probabilistic Perspective 내용을 번역 및 요약한 것임을 밝힌다.

1. Monte Carlo Approximation

일반적으로 변수 변환 공식을 사용하여 확률 변수의 함수의 분포를 계산하는 것은 굉장히 어렵다. 그러나 몬테카를로 근사법은 이에 대한 확실한 대안이 될 수 있다.

어떤 분포로부터 $S$ 개의 Sample $x_1, …, x_s$ 를 추출하였다고 하자. 이 Sample에 기반하여 Empirical Distribution인 $[f(x_i)]_{i=1}^S$ 를 이용하여 우리는 $f(X)$ 라는 분포를 근사할 수 있다. 이를 몬테카를로 적분이라고 부른다.

확률 변수의 어떠한 함수에 대해서도 기댓값을 구할 수 있는 것이 몬테카를로 근사법이다. Sample을 뽑고 이에 대해 Arithmetic Mean을 취하면 된다.

[E[f(X)] = \int f(x) p(x) dx \approx \frac{1}{S} \Sigma_{i=1}^S f(x_i)]

f 함수를 변환시켜 다음과 같은 여러 특성도 얻을 수 있다.

[\bar{x} = \frac {1}{X} \Sigma_{i=1}^S x_i \to E[X]]

[\Sigma_{i=1}^S (x_i - \bar{x})^2 \to var[X]]

[median {x_1, …, x_s} \to median(X)]

2. Example: estimating $\pi$ by Monte Carlo Integration

몬테카를로 적분하면 떠오르는 가장 기본적인 예제이다. 원의 면적을 구해보자. 원의 면적을 $I$ 라고 해보자. Indicator Function을 사용하여 나타내면 아래와 같다.

[I = \int_{-r}^{r} \int_{-r}^{r} \mathbb{I} (x^2 + y^2 \le r^2) dxdy]

Indicatior Funciton을 $f(x, y)$ 라고 하고, 원 밖에 점이 찍히면 0, 안에 찍히면 1의 값을 갖는 함수라고 하자. $p(x), p(y)$ 는 $-r, r$ 구간 내에 존재하는 균일 분포 함수이다. 따라서 $p(x) = p(y) = \frac{1}{2r}$ 이다.

[I = (2r)(2r) \int \int f(x, y) p(x) p(y) dx dy]

[\approx 4r^2 \frac{1}{S} \Sigma_{i=1}^S f(x_i, y_i)]

3. Accuracy of Monte Carlo Approximation

기본적으로 몬테카를로 근사법은 표본의 크기가 클수록 정확도가 높아진다. 실제 평균을 $\mu = E[f(X)]$ 라고 하자. 그리고 몬테카를로 근사법에 의한 근사 평균을 $\hat{\mu}$ 라고 하자. 독립된 표본이라는 가정 하에

[\hat{\mu} - \mu \to \mathcal{N} (0, \frac{\sigma^2}{S})]

[\sigma^2 = var[f(X)] = E[f(X)^2] - E[f(X)]^2]

위는 중심극한정리에 의한 결과이다. 물론 $\sigma^2$ 는 미지의 값이다. 그러나 이 또한 몬테카를로 적분에 의해 근사될 수 있다.

[\hat{\sigma}^2 = \frac{1}{S} \Sigma_{i=1}^S (f(x_i) - \hat{\mu})^2]

[P(\mu - 1.96 \frac{\hat{\sigma}}{\sqrt{S}} \le \hat{\mu} \le \mu + 1.96 \frac{\hat{\sigma}}{\sqrt{S}}) \approx 0.95]

이 때 $\sqrt{ \frac{\hat{\sigma^2}} {S}}$ 는 Numerical(Empirical) Standard Error 라고 불리며, 이는 $\mu$ 의 추정량에 대한 불확실성에 대한 추정치이다. 만약 우리가 95%의 확률로 $\pm \epsilon$ 내부에 위치할 만큼 정확한 Report를 원한다면 우리는 $1.96 \sqrt{ \frac{\hat{\sigma}^2}{S}} \le \epsilon$ 을 만족시키는 $S$ 개의 Sample이 필요하다.

이번 글에서는 통계학을 비롯하여 여러 분야에서 중요하게 사용되는 Expected Maximization 알고리즘에 대해 설명하도록 할 것이다. 본 글은 Bishop의 Pattern Recognition And Machine Learning Chapter9 부분을 번역하여 정리한 것이다. 더욱 상세한 내용을 원한다면 직접 책을 참고하면 좋을 것이다.

EM 알고리즘은 잠재 변수를 가진 모델에서 MLE를 찾는 일반적인 방법이다. Stochastic Gradient Descent와 같은 함수 최적화 방식으로 MLE를 찾을 수 있지만 쉽게 이해하고 구현할 수 있다는 장점을 바탕으로 여전히 널리 사용되고 있는 알고리즘이다.

1. Mixture Models

1.1. Mixture of Gaussians

Gaussian Mixture Model을 이산적인 잠재 변수를 갖는 모델로서 생각해보자. Gaussian Mixture 분포는 아래와 같이 정규분포의 선형 결합으로 표현될 수 있다.

이 때 $\pi_k$ 는 k라는 cluster에 속할 확률을, $K$ 는 Cluster의 개수를 의미한다.

위와 같이 표현이 가능한 이유에 대해 알아보자. 1-of-K representation(하나의 원소만 1이고, 나머지는 0인)을 갖고 $K$ 차원을 가진 이산 확률 변수 $\mathbf{z}$ 가 있다고 하자.

우리는 이제 결합확률 분포 $p(\mathbf{x}, \mathbf{z})$ 를 아래 두 분포의 관점에서 정의하고자 한다. (베이지안 관점에서 생각해보면 된다.)

$\mathbf{z}$ 에 대한 Marginal 분포는 아래와 같이 mixing 계수 $\pi_k$ 의 관점에서 구체화될 수 있다.

[p(z_k = 1) = \pi_k]

$\mathbf{z}$ 가 1-of-K representation을 갖고 있기 때문에 우리는 다음과 같이 $p(\mathbf{z})$ 를 정의할 수 있다.

[p(\mathbf{z}) = \prod_{k=1}^K \pi_k^{z_k}]

유사하게, 특정 $\mathbf{z}$ 값이 주어졌을 때 $\mathbf{x}$ 의 조건부 분포는 아래와 같이 정규분포이다.

이 사실을 이용하여 다음의 결과를 얻을 수 있다.

최종적으로 이번 Chapter 서두에서 보았던 $p(\mathbf{x})$ 를 얻을 수 있다.

각 $x_i$ 에 대해 이에 상응하는 $z_i$ 가 존재한다는 사실을 알 수 있다. 이제 명시적인 잠재 변수를 포함하는 Gaussian Mixture 공식을 얻긴 하였는데, 사실 우리는 $p(\mathbf{x}, \mathbf{z})$ 를 얻고 싶다. 이는 EM 알고리즘을 통해 수행될 수 있다.

Posterior 또한 매우 중요하다. $\mathbf{x}$ 가 주어졌을 때 $z_k = 1$ 일 조건부 확률을 $\gamma(z_k)$ 라고 하자. 베이즈 정리에 의해 다음 식을 도출할 수 있다.

이제 이 $\gamma(z_k)$ 는 관측값 x를 설명하는 데에 있어 component k의 Responsibility로 파악할 수 있다.

1.2. Maximum Likelihood

우리가 $\mathbf{x}_1, …, \mathbf{x}_N$ 의 관측값을 가진 데이터셋을 갖고 있다고 하자. 그리고 우리는 GMM을 이용하여 이 $(N, D)$ 형상의 데이터 $X$를 모델링하고 싶다. 이 때 n번째 행은 $\mathbf{x}_n^T$ 라고 표현될 수 있을 것이다.

유사하게도 이에 상응하는 잠재 변수들은 $(N, K)$ 형상의 $Z$ 행렬로 표현될 수 있고, 이 때 행은 $\mathbf{z}_n^T$ 라고 표현할 수 있을 것이다. 모든 데이터가 분포로부터 독립적으로 추출되었다면, 우리는 GMM을 다음과 같은 그래프 모델로 표현할 수 있다.

이 때 Log Likelihood 함수는 아래와 같이 표현될 수 있다.

이제 MLE를 통해서 이 함수를 최대화 해야 한다. 그런데 그게 그렇게 녹록치만은 않다. 2가지 이슈가 있다.

1) 특이점 문제: Problem of Singularities
2) 혼합의 식별 문제: Indentifiability of Mixtures

먼저 특이점 문제에 대해 살펴볼 것인데, 공분산 행렬이 $\Sigma_k = \sigma^2_k I$ 인 간단한 Gaussiam Mixture에 대해 생각해보자.

이 때 j번째 분포만을 고려한다면 이 분포의 평균은 당연히 $\mu_j$ 일 것인데, 그런데 만약 이 j번째 분포에 오직 1개의 Sample만이 존재한다고 생각해보자.

그렇다면 $\mu_j = \mathbf{x}_n$ 가 될 것이고 Likelihood 함수는 다음과 같이 표현될 수 있다.

만약 이 때 $\sigma_j \to 0$ 이라면, 이 항이 무한을 향해 나아갈 것이고 또한 Log Likelihood 또한 무한이 될 것이다. (Likelihood 함수는 각각의 분포의 곱으로 표현된다.)

따라서 이러한 형태의 Likelihood 함수는 well-posed된 문제라고 할 수 없다. 왜냐하면 이러한 Singularities는 Gaussian Component들 중 하나라도 특정 데이터 포인트에 쏠리기만 해도 이러한 현상은 언제나 발생하기 때문이다.

단 하나의 Gaussian에서는 이러한 문제가 거의 발생하지 않지만, Gaussian Mixture에서는 우연히 하나의 데이터 포인트가 하나의 분포를 이루는 경우가 종종 존재하기 때문에 문제가 되는 것이다.

이러한 Singularity는 Maximum Likelihood 방법론에서 발생할 수 있는 심각한 과적합 문제에 대한 또다른 예시를 제공해주는 셈이다. 우리가 만약 Bayesian Approach를 도입한다면 이러한 문제를 겪지 않을 수 있다. 적절한 휴리스틱을 도입함으로써 이러한 문제를 해결할 수 있는데, 예를 들어 Gaussian Component가 평균이 특정 데이터 포인트에 collapse할 경우 그 component의 평균을 랜덤하게 선택된 값으로 설정하고, 공분산 역시 어떤 큰 값으로 설정한 후 최적화를 다시 진행할 수 있을 것이다.

이제 두 번째 문제인 식별 문제에 대해 알아보자. K개의 파라미터를 K개의 Component에 할당하는 방법은 K!개이기 때문에, K-component Mixture에는 MLE 결과를 동일하게 만드는 K!개의 Solution이 발생하게 된다.

따라서 임의의 한 Sample에 대해 똑같은 분포를 만들어 낼 수 있는 K!-1 개의 데이터 포인트가 존재하게 된다. 이러한 문제를 식별 문제, Identifiability라고 부르며 모델에 의해 밝혀진 파라미터를 해석하는 데에 있어서 굉장히 중요한 이슈이다.

1.3. EM for Gaussian Mixtures

다음 Chapter에서는 EM 알고리즘이 변분 추론의 한 방법임을 밝히면서 좀 더 일반화된 설명을 진행하겠지만, 일단 이번 Chapter에서는 GMM의 맥락에서 설명하도록 하겠다.

GMM을 EM 알고리즘으로 풀어내기 위해서는 Log Likelihood를 최대화하는 파라미터들인 $\mu, \Sigma, \pi$ 를 찾아내야 한다. 지금부터 차근차근 단계를 밟아나가 보도록 하겠다.

Gaussian Component의 평균 $\mu_k$ 에 대한 위 식의 미분 값을 0이라고 설정하면 우리는 다음 식을 얻을 수 있다.

이 때 중간의 분수 식은 $\gamma(z_{nk})$ 임을 기억하자. $\Sigma_k^{-1}$ 를 곱하고 정리하면 다음을 얻을 수 있다. 이 식을 A 식이라고 하자.

[\mu_k = \frac{1}{N_k} \Sigma_{n=1}^N \gamma(z_{nk}) \mathbf{x}_n]

[N_k = \Sigma_{n=1}^N \gamma(z_{nk})]

이 때 $N_k$ 는 Cluster k에 배정된 데이터 포인트의 효과적인 수로 해석할 수 있다. k번째 Gaussian Component의 평균 $\mu_k$ 는 데이터셋에 있는 모든 데이터 포인트의 가중 평균을 취함으로써 얻을 수 있다.

이 때 데이터 포인트 $\mathbf{x}_n$ 의 가중 Factor는,
Component k가 $\mathbf{x}_n$ 을 생성하는 데 Responsible할 사후 확률을 의미하며 그 확률은 아래와 같다.

[\gamma(z_{nk})]

이번에는 $\Sigma_k$ 에 대한 Log Likelihood의 미분 값을 0이라고 설정하고 위와 같은 유사한 과정을 거치면 single Gaussian에 대한 공분산 행렬의 Maximum Likelihood Solution을 아래와 같이 얻을 수 있다. 이 식을 B 식이라고 하자.

[\Sigma_k = \frac{1}{N_k} \Sigma_{n=1}^N \gamma(z_{nk}) (\mathbf{x}_n - \mu_k) (\mathbf{x}_n - \mu_k)^T]

위 식은 데이터셋에 적합한 single Gaussian의 결과와 같은 형상을 취하지만, 각 데이터 포인트는 상응하는 사후 확률에 의해 가중되어 있고, 각 Component의 데이터 포인트의 효과적인 수로 나눠지고 있다.

이제 마지막으로 이 Log Likelihood를 mixing 계수 $\pi_k$ 에 대해 최대화해보자. 이 때 mixing 계수들의 총합은 1이라는 조건이 만족되어야 하는데, 이는 Lagrange Multiplier를 사용하는 것으로 해결할 수 있다.

그래서 위 식을 최대화하면,

이 때 $\pi_k$ 를 양변에 곱해주면,

[\to 0 = \Sigma_{n=1}^N \gamma(r_{nk}) + \lambda]

[\to 0 = N_k + \lambda]

$\lambda = N_k$ 임을 알았다. $\lambda$ 대신 $N_k$ 를 이용하여 식을 다시 정리하면, 아래와 같은 결과를 얻을 수 있다. 이 식은 C 식이라고 하자.

[\pi_k = \frac{N_k}{N}]

해석해보면, k번째 Component의 mixing 계수는 그 Component가 데이터 포인트들을 설명하는 평균 Responsiblity를 의미한다고 볼 수 있다.

이렇게 얻은 $\mu_k, \Sigma_k, \pi_k$ 는 Closed-form의 해를 갖지 못하느데, 왜냐하면 이 파라미터들은 아래 식에서처럼 $\gamma(z_k)$ 에 영향을 받기 때문이다.

따라서 MLE를 얻기 위해서는 Iterative(반복적인) 방법이 제안된다. 일단 처음에는 이 파라미터들에 대한 초깃값을 설정하고, E step과 M step이라 부르는 2가지 단계를 통해 업데이트를 교대로 수행하면 된다.

먼저 E step에서는 위 $\gamma(z_k)$ 에 대한 식을 통해 파라미터의 현재의 값으로 사후 확률을 평가한다. 그리고 이 확률은 M step에서 파라미터들을 A, B, C 식으로 재추정하는데에 사용된다.

정리해보자.

EM for Gaussian Mixtures
1) 평균, 공분산, mixing 계수를 초기화하고 Log Likelihood의 초깃값을 평가한다.
2) E step: 현재의 파라미터 값에 기반하여 Responsibility를 구한다.

3) M step: 위에서 구한 Responsibility와 A, B, C 식을 바탕으로 새로운 파라미터값을 재추정한다.

[\mu_k = \frac{1}{N_k} \Sigma_{n=1}^N \gamma(z_{nk}) \mathbf{x}_n]

[\Sigma_k = \frac{1}{N_k} \Sigma_{n=1}^N \gamma(z_{nk}) (\mathbf{x}_n - \mu_k) (\mathbf{x}_n - \mu_k)^T]

[\pi_k = \frac{N_k}{N}]

4) Log Likelihood를 평가한다.

그리고 파라미터나 Log Likelihood가 수렴하는지 확인한다. 만약 수렴하지 않는다면, 계속해서 2단계로 돌아가서 반복한다.

2. An Alternative View of EM

이번 Chapter에서는 EM 알고리즘에서 잠재 변수가 갖는 핵심적인 역할에 대해 살펴볼 것이다.
다시 설명하지만, EM 알고리즘의 목표는 잠재 변수를 갖는 모델의 Maximum Likelihood를 찾는 것이다. 파라미터가 주어졌을 때 Log Likelihood는 아래와 같다.

물론 연속적인 잠재 변수가 존재하면 위 식에 적분 기호를 적용하면 된다.

중요한 점은, 잠재 변수에 대한 $\Sigma$ 가 Log 안쪽에 있다는 것이다. 이 부분은 Logarithm이 결합 분포에 직접적으로 영향을 미치는 것을 막아 Maximum Likelihood Solution을 찾는 데 있어서 어려움을 가져다 준다.

잠시 Complete 과 Incomplete의 개념을 알아보자. $\mathbf{Z}, \mathbf{X}$ 를 모두 알고 있다면, 이를 Complete 데이터라고 부른다. 만약 잠재 변수에 대해서는 알지 못하고 $\mathbf{X}$ 에 대해서만 안다면, 이는 Incomplete 데이터이다.

이 때 Complete 데이터의 경우 Likelihood 함수는 단지 위 식과 같이 표현되기 때문에 계산이 어렵지 않다.

그러나 실제로 $\mathbf{Z}$ 를 아는 경우는 매우 드물다. 우리가 잠재 변수에 대해 아는 것은 오직 Posterior일 뿐이다. 우리가 Complete-data Log Likelihood를 쓸 수 없기 때문에, 잠재 변수의 사후 분포 하의 이것의 기댓값을 사용하는 것을 생각해 보아야 하는데, 이는 사실 E step에 해당하는 부분이다.

이어지는 M step에서는 이 기댓값을 최대화 한다. 만약 파라미터에 대한 현재의 추정값을 $\theta^{old}$ 라고 한다면, E, M step 이후의 파라미터는 $\theta^{new}$ 로 수정되어 표현할 수 있겠다.

E step 에서는 $\theta^{old}$를 사용하여 아래와 같은 잠재 변수에 대한 Posterior를 찾게 된다.

그리고 이 Posterior를 사용하여 어떤 일반적인 파라미터 값인 $\theta$ 를 찾기 위해 Complete-data Log Likelihood의 기댓값을 계산하게 된다. 이 기댓값은 $\mathcal{Q} (\theta, \theta^{old})$ 라고 표현하며 아래와 같은 형상을 취한다.

그리고 M step에서는 위 식을 최대화하는 $\theta$ 를 찾아 $\theta^{new}$ 로 설정한다.

[\theta^{new} = \underset{\theta}{argmax} \mathcal{Q} (\theta, \theta^{old})]

$\mathcal{Q} (\theta, \theta^{old})$ 에서 Logarithm은 결합 분포에 직접적으로 영향을 미치기 때문에, M step 에서의 최적화는 tractable하다.

요약하자면 아래와 같다.

The General EM Algorithm
1) $\theta^{old}$ 에 대한 초깃값을 설정한다.
2) E step: 아래 사후 확률을 계산한다.

3) M step: $\theta^{new}$ 를 찾고 평가한다.

[\theta^{new} = \underset{\theta}{argmax} \mathcal{Q} (\theta, \theta^{old})]

4) 파라미터 또는 Log Likelihood가 수렴하는지 확인하고, 그렇지 않다면 $\theta^{old} \leftarrow \theta^{new}$ 로 업데이트를 진행한 후 2단계 부터 다시 진행한다.

파라미터에 대해 Prior $p(\theta)$ 가 정의된 모델에 대해 MAP Solution을 찾는 데에도 EM 알고리즘을 쓰일 수 있다. 이 때 E step은 위와 동일하며 M step의 경우 아래와 같이 약간의 수정이 가미된다.

[\mathcal{Q} (\theta, \theta^{old}) + logp(\theta)]

Prior가 절절히 선택된다면 MAP Solution을 찾는 방식으로 Singularities를 제거할 수 있다.

지금까지 잠재 변수에 대해 대응하는 방식으로써 EM 알고리즘을 설명하였는데, 사실 이 알고리즘은 결측값이 존재할 때 역시 적용될 수 있는 알고리즘이다. 이 부분에 대해서는 다른 참고 자료를 확인하길 바란다.

이번 Chapter의 내용을 바탕으로 GMM을 다시 한 번 설명한 부분은 Reference의 책의 442페이지를 참조하길 바란다.

3. The EM algorithm in General

이 Chapter의 주된 내용은 EM 알고리즘을 변분 추론의 관점에서 설명하는 것이다. 변분 추론에 대해서 알고 싶다면 이 글을 참조하라.

지금까지 계속 보았듯이 우리의 목적은 아래 Likelihood를 최대화하는 것이다.

잠재 변수에 대한 분포 $q(\mathbf{Z})$ 를 설정하자. 어떻게든 이 분포에 대해 알 수 있다면 아래와 같은 식은 성립한다.

위 식은 변분 추론에 대한 사전 지식이 있다면 자주 보았던 식일 것이다. 아래와 같은 확률의 곱셈 법칙을 이용하면,

위에서 기술한 ELBO식을 다시 정리할 수 있을 것이다.

EM 알고리즘은 2단계로 MLE를 찾는 알고리즘이다. 위에서 본 분해 식을 활용하여 Log Likelihood를 최대화해보자.

현재의 파라미터를 $\theta^{old}$ 라고 한다면, E step에서는 이를 고정한 채, ELBO $\mathcal{L} (q, \theta)$ 를 $q(\mathbf{Z})$ 에 대해 최대화할 것이다. 이 최대화 해는 쉽게 찾을 수 있다. 그 이유에 대해 설명하겠다.

일단 앞서 확인한 확률의 곱셈법칙으로 인해 $\theta^{old}$ 가 주어졌을 때 아래 식은 $\mathbf{Z}$ 에 의존하지 않는다는 것을 알 수 있다.

그렇다면 ELBO가 최대화된다는 뜻은 결국 KL-divergence가 최소화된다는 뜻이다. 이 경우는 곧 KL-divergence 식에 존재하는 두 분포가 유사해질 때 발생한다는 것을 알 수 있다. (아래 참조)

이제 M step 에서는 $q(\mathbf{Z})$ 를 위 식과 같이 고정하고 ELBO를 최대화하여 $\theta^{new}$ 를 얻는다. 이는 결국 Log Likelihood를 증가시키는 길이 될 것이다.

조금 더 상세히 알아보자. ELBO는 아래와 같은 식인 것을 확인했다.

자 이제 $q(\mathbf{Z})$ 를 고정해보자.

[= \mathcal{Q} (\theta, \theta^{old}) + Const]

위 상수는 사실 그냥 q 분포에 대한 Negative Entropy이다. 결국 M step 에서 최대화하고 있는 것은 이전 Chapter에서 설명하였던 것처럼, Complete-Data Log Likelihood인 것이다. 물론 만약 아래와 같은 결합 분포가 Exponential Family에 속한다면 Log가 Exponential을 cancel하고 계산은 훨씬 간단해질 것이다.

M step을 그림으로 나타내면 아래와 같다.

다음 그림은 EM 알고리즘의 전반적인 과정을 그림으로 표현한 것이다. 이와 같이 EM 알고리즘은 파라미터의 공간으로 해석할 수 있다. 이러한 측면 덕분에 EM 알고리즘은 Bayesian Optimization에서 Acqusiition Function으로 널리 사용되고 있다.

빨간 선이 우리가 최대화하고 싶은 Complete-Data Log Likelihood이다. 파란 선은 초깃값으로 설정한 ELBO이다. 지속적으로 EM을 수행하면 초록 선과 같은 결과를 얻을 수 있다.

References

Bishop, Pattern Recognition And Machine Learning, Chapter9

본 글에서는 정보 이론과 관련이 있는 Kullback-Leibler Divergence와 이를 기반으로 한 Variational Inference에 대해 간략히 정리해보고자 한다. 시작에 앞서, 변분 추론은 근사 추정의 대표적인 방법이라는 점을 밝히고 싶으며, 본 글에서 소개된 변분 추론 기법은 Vanilla Variational Inference라고 부를 수 있는 CAVI이다. CAVI의 단점을 보완한 다양한 변분 추론 기법들이 연구되었으며 이에 대한 내용은 후속 글에서 다루도록 할 것이다.

1. Kullback-Leibler Divergence

정보이론에서 정보량은 불확실성이 커질수록 많아지는 것으로 정의한다. Shannon Entropy는 확률의 값에 $log$ 를 씌우고 -1을 곱해준 값으로, 모든 사건의 정보량의 Expectation을 의미한다. 확률 분포 P에 대한 섀넌 엔트로피는 아래와 같이 정의할 수 있다.

[H(P) = E_{X \sim P}[-logP(x)]]

Kullback-Leibler Divergence는 분포의 유사한 정도를 나타내는 척도인데, 비대칭적이기 때문에 거리를 나타내는 척도로는 사용될 수 없다. 그보다는 정보량의 차이를 나타내는 것이라고 이해하면 좋다.

어떤 확률 분포 $P(x)$ 가 주어졌을 때, 이 분포와 $Q(x)$ 라는 분포와의 차이를 알고 싶다면 쿨백-라이블리 발산 값을 구하면 된다. 이 값은 Cross-Entropy에서 Self-information을 뺀 값으로 정의되는데, Cross-Entropy는 아래와 같이 정의된다.

[H(P, Q) = E_{X \sim P}[-log(Q(x))] = -\Sigma_x P(x) log Q(x)]

$H(P)$ 로 표현되는 Self-Information은 그 확률변수가 갖고 있는 정보량을 의미하는데,

이산확률변수라면 $-\Sigma_x P(x) logP(x)$로 표현할 수 있다.

이제 쿨백-라이블리 발산 식을 보자.

[= -\Sigma_x P(x) log Q(x) + \Sigma_x P(x) logP(x)]

[= -\Sigma_x [P(x) log Q(x) - P(x)logP(x)]]

[= -\Sigma_x P(x) [logQ(x)-logP(x)]]

[= -\Sigma_x P(x) log(\frac{Q(x)}{P(x)})]

[= E_{X \sim P}[-log \frac{Q(x)}{P(x)}]]

쿨백-라이블리 발산을 최소화하는 것은 기본적으로 두 확률 분포 사이의 차이를 최소화하는 것이다. 일반적으로 P를 우리가 갖고 있는 데이터의 분포, 혹은 계산하기 편리한 쉬운 분포라고 한다면 Q는 모델이 추정한 분포 또는 확률적으로 계산하기 어려운 분포라고 생각할 수 있다.

머신러닝 관점에서 생각해본다면 Self-Information 부분은 학습 과정에 있어서 어떤 변동을 겪지는 않는다. 따라서 쿨백-라이블리 발산을 최소화하는 것은 크로스-엔트로피를 최소화하는 것과 동일한 의미를 지닌다.

두 분포가 동일하면 쿨백-라이블리 발산은 0의 값을 가지고 다를수록 점점 큰 값을 가진다.

2. Variational Inference

2.1. ELBO

$\mathbf{x}$ 란 확률 변수가 있고, 이 변수의 특성의 상당 부분은 잠재 변수인 $\mathbf{z}$ 에 의해 설명된다고 하자. 우리는 이 때 우리는 $\mathbf{x}$ 의 실현 값인 데이터가 존재할 때 $\mathbf{z}$ 의 분포, 즉 Posterior $p(\mathbf{z}|\mathbf{x})$ 를 알고 싶다. 그런데 Posterior는 많은 경우에 Numerical 계산이 불가능하다. 따라서 우리는 이 Posterior를 알기 쉬운 분포 $q(\mathbf{z})$ 로 바꾸고 싶다.

변분추론은 이렇게 Posterior를 다루기 쉬운 분포로 근사하는 방법론을 의미한다.

이 때 $q(z)$ 는 어떤 함수의 집합 $Q$ 의 한 원소라고 생각할 수 있다.

용어를 잠시 정리해보자.

$p(x)$ 는 Marginal Probability 또는 Evidence를, $p(z)$ 는 Prior를 의미한다.
$p(x|z)$ 는 Likelihood를, $p(z|x)$ 는 Posterior를 의미한다.

위에서 확인한 쿨백-라이블리 발산을 이용하여 변분추론을 설명하면, 변분추론은 아래와 같이 쿨백-라이블리 발산 값을 최소화하는 $Q$ 집합 내의 q 함수를 찾는 것이 된다.

이전 Chapter에서 쿨백-라이블리 발산은 분포의 유사도를 측정하는 Index라고 하였다. 이 식을 위 관점을 다시 한 번 분석해보자.

이렇게 사후확률과 다루기 쉬운 분포 사이의 쿨백-라이블리 발산은 총 3가지 항으로 분해할 수 있다.

쿨백-라이블리 발산은 그 정의 때문에 0 이상의 값을 가진다. 즉 Non-negative이다. 따라서 제일 마지막 줄은 아래와 같이 다시 표현할 수 있다. (이 부분은 Jensen의 부등식을 통해서도 추론할 수 있다.)

우항을 ELBO(Evidence Lower BOund)라고 부른다. Evidence의 하한선이라는 의미이다.

Variational Density $q(z)$ 와 Posterior 사이의 쿨백-라이블리 발산 값 부터 다시 표현해보면,

위 식을 보면 Evidence는 ELBO와 쿨백-라이블리 발산의 합으로 구성된다는 것을 알 수 있다. 이는 매우 중요한 수식이다. 일단은 쿨백-라이블리 발산을 최소화하는 것은 곧 ELBO를 최대화하는 것과 의미가 같다는 것은 쉽게 파악할 수 있다. 다만 ELBO와 쿨백-라이블리 발산 모두 q 함수에 의존적이기 때문에 단순히 한 쪽을 최소화하는 q 함수를 찾았다고 해서 이것이 반드시 Evidence의 값을 최소화한다고 말하기는 어렵다는 부분은 잊으면 안된다.

다음 Chapter 부터는 최적화 방법에 대해 설명할 것인데, 전통적인 방법인 CAVI를 통해 설명을 진행하도록 하겠다.

2.2. CAVI: Cooridinate Ascent mean field Variational Inference

(Mean-Field Variational Inference)

그렇다면 q 함수는 대체 어떤 함수인가? 아주 클래식한 방법으로 설명하자면, Mean Field Variational Family를 언급해야 할 것이다. 잠재 변수 $\mathbf{z}$ 가 모두 독립적이라고 할 때, q 함수는 아래와 같이 분해될 수 있다.

[q(\mathbf{z}) = \prod_j q_j(z_j)]

이렇게 Variational 분포를 각각의 곱으로 분해하고 나면, 우리는 각 Factor에 대해 Coordinate Ascent Optimization을 적용할 수 있다. CAVI는 한 쪽을 고정한 채, Mean-Field Variational Density의 각 Factor를 반복적으로 최적화한다. 이 알고리즘은 ELBO의 Local Optimum으로 이끈다.

ELBO를 다시 확인해보자.

[= \int_z q(z) log p(\mathbf{x}, z) - q(z)logq(z) dz]

[= E_{z \sim q(z)} [logp(\mathbf{x}, z)] - E_{z \sim q(z)}[logq(z)]]

CAVI의 핵심 아이디어는 q 함수가 분해될 수 있다는 사실을 이용하는 것이다. $z_j$ 를 j번째 잠재 변수라고 하자. (아래의 - 기호는 그 index를 제외한 것을 의미한다.) 이 때 우리가 알고 싶은 것은, $\mathbf{x}$ 와 $\mathbf{z}_{-j}$ 가 모두 주어졌을 때 $z_j$ 의 Complete 조건부 확률이다. 이는 아래와 같이 표현할 수 있다.

그런데 앞서 했던 가정에 따라 모든 잠재 변수는 독립적이다. 따라서 위 식은 아래와 같다.

[= logp(z_j, \mathbf{z}_{-j}, \mathbf{x})]

지금부터 이 사실을 염두에 두고 위에서 보았던 ELBO 식을 $q_j$ 의 관점에서 풀어쓸 것이다. 아래에서 나오는 $l$ 기호는 $j$ 가 아닌 Index를 의미한다. (j번째 잠재 변수가 아닌 나머지 Variational Factors: $q_l(\mathbf{z}_l)$ )

[ELBO = E_q[logp(\mathbf{x}, z_j, \mathbf{z}{-j})] - E{q_l}[logq_l(\mathbf{z}_l)]]

Iterative Expectation을 이용하면,

첫 항의 안쪽 부분을 보자. 기댓값의 정의에 따라 다음과 같이 식을 전개할 수 있다.

[= \int_{-j} logp(\mathbf{x}, z_j, \mathbf{z}{-j}) q(\mathbf{z}{-j}) dq_{-j}]

[= E_{-j} [logp(\mathbf{x}, z_j, z_{-j})]]

최종적으로 $q_j$ 에 대한 ELBO는 아래와 같다.

[ELBO = E_j[ E_{-j} [logp(\mathbf{x}, z_j, \mathbf{z}_{-j})]] - E_j[logq_j] + Const]

첫 번재 항을 최대로 하는 것이 $q_j$ 에 대한 ELBO를 최대화하는 길이다. 따라서 $q_j$ 에 대한 Optimal Solution은 아래와 같이 표현할 수 있다.

[q^*j{z_j} \propto exp( E{-j} [logp(\mathbf{x}, z_j, \mathbf{z}_{-j})] )]

쿨백-라이블리 발산이 아래와 같이 표현되었다는 사실을 기억해보자.

위 쿨백-라이블리 발산의 개념을 적용해보면, $q_j$ 에 대한 ELBO 식은 $q^{*}_j z_j, q_j(z_j)$ 사이의 Negative 쿨백-라이블리 발산 값을 의미한다.

따라서 이를 해석해보면, j번째 잠재변수의 Variational Density를 j번째 잠재변수의 최적화된 Variational Density와 유사하게 만드는 것이 $q_j$ 의 ELBO를 최대화하는 것이고, 이러한 과정을 모든 j에 대해, ELBO가 수렴할 때까지 반복한다면 우리가 원하는 q 함수를 얻을 수 있다는 의미가 된다.

지금까지 설명한 부분을 정리해보자.

CAVI의 일반적인 절차는 아래와 같다.
1) Variational 분포 q를 설정한다.
2) 각 잠재 변수의 Gradient를 잡아 각 $q_j$ 를 최적화한다.
3) ELBO를 계산한다.
4) ELBO가 수렴할 때 까지 위 과정을 반복한다.

CAVI는 클래식하고 좋은 방법론이지만, Non-Convex 최적화 문제에서 Global Optimum에 도달할 것이라고 보장해주지는 못한다. 즉, 충분히 쿨백-라이블리 발산 값을 최소화하지 못할 수도 있다는 뜻이다. 또한 MCMC와 같은 Posterior Estimation 보다는 (최적화 방법이기에) 속도가 빠르지만 한 쪽을 고정하고 다른 쪽을 교대로 계산하는 방법을 채택하고 있기 때문에 기본적으로 속도가 아주 빠르지는 않다는 단점도 지니고 있다.

이에 대한 보완책으로 여러 방법론이 대두되었는데, Stochastic Gradient Descent, Convex Relaxation, Monte Carlo Sampling 등의 개념을 활용한 알고리즘들이 등장하였다. 글 서두에서도 밝혔듯이 이러한 알고리즘들에 대해서는 후속 글에서 다루도록 하겠다.

Reference

1) 변분추론 설명 블로그1
2) 변분추론 설명 블로그2
3) 패턴인식-머신러닝 책 정리 사이트
4) 변분추론 논문

Gaussian Process에 대해 알아보자!

Gaussian Process는 Random(Stochastic) Process의 한 예이다. 이 때 Random Process는 시간이나 공간으로 인덱싱된 Random Variable의 집합을 의미한다. GP의 정의는 아래와 같다.

Stochastic process such that every finite collection of those random variables has a multivariate normal distribution.

이는 또한 이 Random Variable들의 선형 결합 일변량 정규분포를 따른다는 말과 동일한 설명이고, GP의 일부를 가져와도 이는 항상 다변량 정규분포를 따른다는 것을 의미한다. 또한 다변량 정규분포는 주변 분포와 조건부 분포 역시 모두 정규분포를 따르기 때문에 이를 계산하기 매우 편리하다는 장점을 지닌다.

GP는 일종의 Bayesian Non-parametric method으로 설명되는데, 이 때 Non-parametric이라는 것은 parameter의 부재를 의미하는 것이 아니라 parameter가 무한정 (infinite) 있다는 것을 의미한다.

지금부터는 GP에 대해 이해하기 위해 단계적으로 설명을 진행할 것이다.

1. Basics of Gaussian Process

다변량 정규분포를 생각해보자.

[p(x, y) \sim \mathcal{N}( \begin{bmatrix} \mu_x \ \mu_y \end{bmatrix}, \begin{bmatrix} \Sigma_x \Sigma_{xy} \ \Sigma_{xy}^T \Sigma_y \end{bmatrix})]

앞서 언급하였듯이 이 다변량 정규분포를 이루는 확률변수의 어떠한 부분집합에 대해서도 주변 분포와 조건부 분포 모두 정규분포를 따른다. GP는 여기서 한발 더 나아가서, 이러한 다변량 정규분포를 무한 차원으로 확장시키는 개념으로 생각하면 된다.

이 무한의 벡터를 일종의 함수로 생각할 수도 있을 것이다. 연속형 값을 인풋으로 받아들이는 함수를 가정하면, 이는 본질적으로 input에 의해 인덱싱된 값들을 반환하는 무한한 벡터로 생각할 수 있다. 이 아이디어를 무한 차원의 정규분포에 적용하면 이것이 바로 GP의 개념이 된다.

따라서 Gaussian Process는 함수에 대한 분포라고 표현할 수 있다. 다변량 정규분포가 평균 벡터와 공분산 행렬로 표현되는 것처럼, GP 또한 평균 함수와 공분산 함수를 통해 다음과 같이 정의된다.

[P(X) \sim GP(m(t), k(x, x\prime))]

GP에 있어서 Marginalization Property는 매우 중요한 특성이다. 우리가 관심 없거나 관측되지 않은 수많은 변수에 대해 Marginalize할 수 있다.

GP의 구체적 예시를 다음과 같이 들 수 있을 것이다. 실제로 이는 가장 흔한 설정이다.

[m(x) = 0]

[k(x, x\prime) = \theta_1 exp( - \frac{\theta_2}{2} ( x - x\prime)^2 )]

여기서 공분산 함수로 Squared Exponential을 사용하였는데, $x$와 $x\prime$이 유사한 값을 가질 수록 1에 수렴하는 함수이다. (거리가 멀수록 0에 가까워짐) 평균 함수로는 0을 사용하였는데, 사실 평균 함수로 얻을 수 있는 정보는 별로 없기에 단순한 설정을 하는 것이 가장 편리하다.

유한한 데이터 포인트에 대해 GP는 위에서 설정한 평균과 공분산을 가진 다변량 정규분포가 된다.

다음 Chapter부터는 본격적으로 이론에 대한 부분을 정리하도록 하겠다. 2개의 논문을 정리하였는데, 첫 번째 논문은 Gaussian Process의 가장 기본적이고 중요한 내용을 담은 논문이며, 두 번째 논문은 좀 더 개념을 확장하고 직관적으로 Gaussian Process Regression에 대해 서술한 논문이다.

2. Gaussian Process in Machine Learning

본 논문은 GP가 회귀를 위한 Bayesian 프레임워크를 형성하기 위해 어떻게 사용되는지, Random(Stochastic) Process가 무엇이고 이것이 어떻게 지도학습에 사용되는지를 설명하는 것이 주 목적이다. 또한 공분산 함수의 Hyperparameter 설정에 관한 부분, 그리고 주변 우도와 Automatic Occam’s Razor에 관한 이야기도 포함한다.

(Occam’s Razor: 오캄의 면도날 원칙, 단순함이 최고다.)

2.1. Posterior Gaussian Process

GP는 함수에 대한 분포로 정의되며, 이러한 GP는 베이지안 추론의 Prior로 사용된다. 이 Prior는 학습 데이터에 의존하지 않으며 함수들에 대한 어떤 특성을 구체화한다. Posterior Gaussian Process의 목적은 학습데이터가 주어졌을 때, 이 Prior를 업데이트하는 방법을 도출해내는 것이다. 나아가 새로운 데이터가 주어졌을 때 적절히 예측 값을 반환하는 것이 목표가 될 것이다.

기존의 학습데이터와 새로운 테스트 데이터를 분리하여 다음과 같은 결합 분포를 상정해보자.

[\begin{bmatrix} \mathbf{f} \ \mathbf{f} \end{bmatrix} \sim \mathcal{N}( \begin{bmatrix} \mathbf{\mu} \ \mathbf{\mu_}\end{bmatrix}, \begin{bmatrix} \Sigma, \Sigma_* \ \Sigma_*^T, \Sigma_{**} \end{bmatrix})]

[\mathbf{\mu} = m(x_i), i = 1, … , n]

이제 우리가 알고 싶어하는 f*의 조건부 분포는 아래와 같은 형상을 지녔다. 아래 식은 테스트 데이터에 대한 사후분포에 해당한다.

이와 같은 분포를 얻을 수 있는 이유는 결합 정규분포를 조건화하는 공식인 다음의 결과에 기인한다.

위에서 확인한 사후분포에 기반하여 Posterior Process를 구해보면 아래와 같다.

이 때 $\Sigma(X, x)$ 는 모든 학습 데이터와 $x$ 의 공분산 벡터를 의미한다. 이제 위 식을 자세히 뜯어보자. Posterior Process의 공분산 함수는 Prior의 공분산 함수에서 양의 값을 뺀 것과 같다. 즉 Posterior Process의 공분산 함수는 Prior의 그것보다 언제나 작은 값을 가진다는 의미이다. 이것은 논리적인데, 데이터가 우리에게 정보를 제공하였기 때문에 Posterior의 분산이 감소하는 것이다.

자 이제 학습 데이터의 Noise를 고려해야 한다. 이에 대해서 정규 분포를 설정하는 것이 일반적이다. Noise를 고려한 후 다시 정리하면 아래와 같다.

이제 Posterior Process에서 샘플링을 진행할 수 있다. 이렇게 평균 함수와 공분산 함수를 정의함으로써 학습 데이터가 주어졌을 때 Prior를 Posterior로 업데이트할 수 있게 되었다. 그러나 문제가 찾아왔다. 어떻게 평균 함수와 공분산 함수를 적절히 설정하는가? 그리고 Noise Level( $\sigma_n^2$ )은 어떻게 추정하는가?

2.2. Training a Gaussian Process

사실 일반적인 머신러닝 적용 케이스에서 Prior에 대해 충분한 정보를 갖고 있는 것은 드문 경우이다. 즉, 평균 함수와 공분산 함수를 정의하기에는 갖고 있는 정보가 부족하다는 것이다. 우리는 갖고 있는 학습 데이터에 기반하여 평균, 공분산 함수에 대해 적절한 추론을 행해야 한다.

Hyperparameter에 의해 모수화되는 평균, 공분산 함수를 가진 Hierarchical Prior를 사용해보자.

[f \sim \mathcal{GP}(m, k)]

[m(x) = ax^2 + bx + c]

[k(x, x\prime) = \sigma_{y}^2 exp(- \frac{( x - x\prime)^2}{2l^2} ) + \sigma_n^2 \delta_{ii\prime}]

이제 우리는 $\theta=[a, b, c, \sigma_y, \sigma_n, l]$ 이라는 Hyperparameter 집합을 설정하였다. 이러한 계층적 구체화 방법은 vague한 Prior 정보를 간단히 구체화할 수 있게 해준다.

우리는 데이터가 주어졌을 때 이러한 모든 Hyperparameter에 대해 추론을 행하고 싶다. 이를 위해서는 Hyperparameter가 주어졌을 때 데이터의 확률을 계산해야 한다. 이는 어렵지 않다. 주어진 데이터의 분포는 정규 분포임을 가정했기 때문이다. 이제 Log Marginal Likelihood를 구해보자.

이제 편미분 값을 통해 이 주변 우도를 최적화(여기서는 최소화)하는 Hyperparameter의 값을 찾을 수 있다. 아래에서 $\theta_m$ 와 $\theta_k$ 는 평균과 공분산에 관한 Hyperparameter를 나타내기 위한 parameter이다.

위 값들은 Conjugate Gradients와 같은 Numerical Optimization에 사용된다.

GP는 Non-parametric 모델이기 때문에 Marginal Likelihood의 형상은 Parametric 모델에서 보던 것과는 사뭇 다르다. 사실 만약 우리가 Noise Level인 $\sigma_n^2$ 를 0으로 설정한다면, 모델은 정확히 학습 데이터 포인트와 일치하는 평균 예측 함수를 생성할 것이다. 하지만 이것은 주변 우도를 최적화하는 일반적인 방법이 아니다.

Log Marginal Likelihood 식은 3가지 항으로 구성되어 있는데, 첫 번째는 Complexity Penalty Term으로 모델의 복잡성을 측정하고 이에 대해 페널티를 부과한다. Negative Quadratic인 두 번째 항은 데이터에 적합하는 역할을 수행하며, 오직 이 항만이 학습 데이터의 Output인 $\mathbf{y}$ 에 의존적이다. 세 번째 항은 Log-Normalization Term으로 데이터에 대하여 독립적이며 사실 뭐 그리 중요한 항은 아니다.

GP에서 페널티와 데이터에 대한 적합의 trade-off는 자동적이다. 왜냐하면 Cross Validation과 같은 외부적 방법이 세팅될 Parameter가 존재하지 않기 때문이다. 실제로 이와 같은 특성은 보통의 머신러닝 알고리즘 상에 존재하는 Hyperparameter 튜닝에 소요되는 시간을 절약하게 해주기 때문에 학습을 더욱 간단하게 만드는 장점을 갖게 된다.

2.3. Conclusions and Future Directions

본 논문에서는 GP가 굉장히 변동성이 크고 유연한 비선형적 회귀를 구체화하는 데 편리하게 사용되는 과정에 대해 알아보았다. 본 논문에서는 오직 1가지 종류의 공분산 함수가 사용되었지만, 다른 많은 함수들이 사용될 수 있다. 또한 본 논문에서는 오직 가장 간단한 형태인 정규 분포의 Noise를 가정하였지만, 그렇지 않을 경우 학습은 더욱 복잡해지고 Laplace 근사와 같은 방법이 도입되거나 Sampling이 이루어져야만 non-Gaussian Posterior를 정규분포와 유사하게 만들 수 있다.

또 중요한 문제점은 계산 복잡성이다. 공분산 행렬의 역행렬을 구하기 위해서는 메모리 상에서는 $\mathcal{O}(n^2$ 의 복잡도가, 계산 상에서는 $\mathcal{O}(n^3)$ 의 복잡도가 발생한다. 리소스에 따라 다르지만, 행이 10,000개만 넘어가도 직접적으로 역행렬을 계산하기에는 많은 무리가 따른다. 따라서 근사적인 방법이 요구되는데, 본 논문이 나온 시점이 2006년임을 고려하면, 이후에도 많은 연구가 진행되었음을 짐작할 수 있을 것이다.

한 예로 이 논문이 있는데, 추후에 다루도록 할 것이다.

3. Gaussian Process Regression

본 Chapter에서는 두 번째 논문을 기반으로 좀 더 단계적으로 설명을 해볼 것이다.

논문의 내용을 설명하기 전에 전체적인 구조를 다시 한번 되짚어보도록 하자.

3.1. Overview

비선형 회귀 문제를 생각해보자. 우리는 데이터가 주어졌을 때 이를 표현하는 어떤 함수 f를 학습하고 싶고 이 함수는 확률 모델이기 때문에 신뢰 구간 또는 Error Bar를 갖게 된다.

[Data: \mathcal{D} = [\mathbf{x}, \mathbf{y}]]

Gaussian Process는 평균 함수와 공분산 함수를 통해 이 함수에 대해 분포를 정의한다. 이 함수는 Input Space $\mathcal{X}$ 를 $\mathcal{R}$ 로 mapping하며, 만약 두 공간이 정확히 일치할 경우 이 함수는 infinite dimensional quantity가 된다.

[p(f) = f(x) \sim \mathcal{GP}(m, k)]

그리고 베이즈 정리에 따라 위 확률은 Bayesian Regression에 사용된다.

Posterior를 구하기 위해서는 당연히 Prior와 Likelihood가 필요한데, 이 때 Prior는 Gaussian Process를 따른다고 가정한다. 이제 Likelihood를 구해야 한다.

우리가 수집한 데이터 $\mathcal{D}$ 는 일반적으로 Noise를 포함하고 있다. 따라서 우리의 정확한 목표는 $f(x)$ 를 추정하는 것이 아니라 Noise를 포함한 $y$ 를 추정하는 것이어야 한다. 평균 함수를 0으로 가정하고 $y$ 를 비롯하여 GPR에 필요한 모델들에 대해 정리해보자.

[y = f(x) + \epsilon]

[\epsilon \sim \mathcal{N}(0, \sigma_n^2)]

다음 Chapter에서도 나오겠지만 이 Noise의 분산을 공분산 함수 속으로 집어넣을 수 있다. (자세한 수식은 다음 Chapter를 참조하라) 그러면 사실 아래의 $K$ 는 $K + \sigma_n^2$ 를 의미하게 된다.

[f \sim \mathcal{GP}(0, K)]

f의 Prior는 GP고, Likelihood는 정규분포이므로 f에 대한 Posterior 또한 GP이다. 일단 주어진 데이터에 기반하여 Marginal Likelihood를 구해보자.

[= \mathcal{N}(0, K)]

그런데 이 때 이전 Chapter와 마찬가지로 공분산 함수를 정의할 때 사용되는 Hyperparameter로 $\theta$ 를 정의하게 되면, Marginal Likelihood는 정확히 아래와 같이 표현할 수 있다.

이 식에 Log를 취해서 다시 정리하면 Log Marginal Likelihood가 된다. ( $\theta$ subscript는 생략한다.)

Numeric한 방법으로 위 목적 함수를 최적화(최소화)하는 $\theta$ 를 구하면 이는 공분산 함수의 최적 Hyperparameter가 된다. 이제 예측을 위한 분포를 확인해보자. 새로운 데이터 포인트 $x_*$ 가 주어졌을 때의 예측 값에 관한 사후분포이다.

[= \mathcal{N}( K_K^{-1}\mathbf{y}, K_{**} - K_ K^{-1} K_*^T )]

이제 위 분포를 바탕으로 Sampling을 진행하고, 평균과 분산을 바탕으로 그래프를 그리면 본 글의 가장 서두에서 본 것과 같은 아름다운 그래프를 볼 수 있다.

평균인 $K_*K^{-1}\mathbf{y}$ 는 다음과 같이 $\mathbf{y}$ 에 대한 선형결합으로 표현할 수도 있다.

[K_K^{-1}\mathbf{y} = \Sigma_{i=1}^n \alpha_i k(x_i, x_), \alpha = K^{-1}\mathbf{y}]

지금까지 설명한 내용이 바로 Gaussian Process를 Function Space View로 이해한 것이다.

3.2. Definition of Gaussian Process

지금부터는 논문의 내용을 정리한 것이다. 사실 GP의 기본적인 설명은 끝났다고 봐도 무방하지만, 그럼에도 이 세심한 논문의 설명을 다시 한 번 읽어보지 않을 수가 없다. 정의에 대한 부분은 처음에 설명하였으므로 생략하도록 하겠다.

Chapter1에서 공분산 함수 $ k(x, x\prime) $에 대해서 설명하였는데, 본 논문에 맞추어 Notation을 살짝 변형하도록 하겠다. (이전 Chapter에서는 이 공분산 함수를 가장 단순한 버전인 $\Sigma$ 로 표현하였다.)

[k(x, x\prime) = \sigma_f^2 exp( - \frac{( x - x\prime)^2}{2l^2} )]

정말 기호만 살짝 바뀌었다. $x$가 $x\prime$과 유사할 수록 $f(x)$라 $f(x\prime)$과 상관성(Correlation)을 가진다고 해석할 수 있다. 이것은 좋은 의미이다. 왜냐하면 함수가 smooth해지고 이웃한 데이터 포인트끼리 더욱 유사해지기 때문이다.

만약 그 반대의 경우 2개의 데이터 포인트는 서로 마주칠 수도 없다. 즉 새로운 $x$값이 삽입될 때, 이와 먼 곳에 있는 관측값들은 그다지 큰 영향을 미칠 수 없다. 이러한 분리가 갖는 효과는 사실 length parameter인 $l$에 달렸있는데, 이 때문에 이 공분산 함수는 상당한 유연성을 지니는 식이 된다.

하지만 데이터는 일반적으로 Noise를 포함하고 있다. 이 때문에 언제나 측정 오차는 발생하기 마련이다. 따라서 $y$ 관측값은 $f(x)$에 더불어 Gaussian Noise를 포함하고 있다고 가정하는 것이 옳다.

[y = f(x) + \mathcal{N}(0, \sigma_n^2)]

많이 보았던 회귀식 같아 보인다. 이 Noise를 공분산 함수안에 집어넣으면 아래와 같은 형식을 갖추게 된다.

[k(x, x\prime) = \sigma_f^2 exp(- \frac{( x - x\prime)^2}{2l^2} ) + \sigma_n^2 \delta(x, x\prime)]

여기서 $\delta(x, x\prime)$은 Kronecker Delta Function이다.

많은 이들은 GP를 사용할 때 $\sigma_n$을 공분산 함수와 분리해서 생각하지만 사실 우리의 목적은 y* 를 예측하는 것이지 정확한 f* 를 예측하는 것이 아니기 때문에 위와 같이 설정하는 것이 맞다.

Gaussian Process Regression을 준비하기 위해 모든 존재하는 데이터포인트에 대해 아래와 같은 공분한 함수를 계산하도록 하자.

$K$의 대각 원소는 $\sigma_f^2 + \sigma_n^2$ 이고, 비대각 원소 중 끝에 있는 원소들은 $x$ 가 충분히 큰 domain을 span할수록 0에 가까운 값을 갖게 된다.

3.3. How to Regress using Gaussian Process

GP에서 가장 중요한 가정은 우리의 데이터가 다변량 정규 분포로부터 추출된 Sample로 표현된다는 것이므로 아래와 같이 표현할 수 있다.

[\begin{bmatrix} \mathbf{y} \ y* \end{bmatrix} \sim \mathcal{N}(0, \begin{bmatrix} K, K_^T \ K_, K_{**} \end{bmatrix})]

우리는 물론 조건부 확률에 대해 알고 싶다.

이 확률은 데이터가 주어졌을 때 $y_*$ 에 대한 예측의 확실한 정도를 의미한다.

정규분포이므로, $y_*$ 에 대한 Best Estimate는 평균이 될 것이다.

[\bar{y}* = K*K^{-1}\mathbf{y}]

그리고 분산 또한 아래와 같다.

[var(y_) = K_{**} - K_ K^{-1} K_*^T]

이제 본격적으로 예제를 사용해보자. Noise가 존재하는 데이터에서 다음 포인트 $x_*$ 에서의 예측 값은 얼마일까?

6개의 관측값이 다음과 같이 주어졌다.

x = [-1.5, -1, -0.75, -0.4, -0.25, 0]

Noise의 표준편차 $\sigma_n$ 이 0.3이라고 하자. $\sigma_f$ 와 $l$ 을 적절히 설정하였다면 아래와 같은 행렬 $K$를 얻을 수 있다.

공분산 함수를 통해 아래 사실을 알 수 있다.

[K_{**} = 3]

[K_* = [0.38, 0.79, 1.03, 1.35, 1.46, 1.58]]

[\bar{y}_* = 0.95]

[var(y_*) = 0.21]

[x* = 0.2]

그런데 매번 이렇게 귀찮게 구할 필요는 없다. 엄청나게 많은 데이터 포인트가 존재하더라도 이를 한번에 큰 $K_*$ 과 $K_{**}$ 을 통해 계산해버리면 그만이다.

만약 1000개의 Test Point가 존재한다면 $K_{**}$ 는 (1000, 1000)일 것이다.

95% Confidence Interval은 아래 식으로 구할 수 있고 이를 그래프로 표현하면 아래 그림과 같다.

[\bar{y}* \pm 1.96\sqrt{var(y*)}]

3.4. GPR in the Real World

이전 Chapter에서 보았던 내용이 신뢰를 얻기 위해서는 사실 우리가 얼마나 공분산 함수를 잘 선택하느냐에 달려있다. $\theta = [l, \sigma_f, \sigma_n]$ 라는 Parameter 집합이 적절히 설정되어야만 결과가 합리적일 것이다.

$\theta$ 의 Maximum a Posteriori Estimate는 다음 식이 최댓값을 가질 때 찾을 수 있다.

베이즈 정리에 따라 우리가 $\theta$ 에 대해 거의 아는 것이 없다고 가정할 때 우리는 다음과 같은 식을 최대화해야 한다.

다변량 최적화 알고리즘(예: Conjugate Gradients, Nelder-Mead simplex)을 이용하면 예를 들어 $l=1, \sigma_f=1.27$ 과 같은 좋은 값을 얻을 수 있다.

그런데 이건 그냥 단지 좋은 값 에 불과하다. 수많은 옵션 중에 딱 하나 좋은 답이 있으면 안되는가? 이 질문에 대한 답은 다음 장에서 찾을 수 있다.

좀 더 복잡한 문제에 대해 생각해보자. 아래와 같은 Trend를 갖는 데이터가 있다고 하자.

좀 더 복잡한 공분한 함수가 필요할 것 같다.

[k(x, x\prime) = \sigma_{f_1}^2 exp(- \frac{( x - x\prime)^2}{2l_1^2} ) + \sigma_{f_2}^2 exp(- \frac{( x - x\prime)^2}{2l_2^2} ) + \sigma_n^2 \delta(x, x\prime)]

위 식의 우항에서 첫 번째 부분은 예측변수의 작은 변동을 잡아내기 위함이고, 두 번째 부분은 좀 더 긴 기간 동안의 변동성을 포착하기 위해 설계되었다. ( $l_2 \approx 6l_1$ )

이 공분산 함수는 $K$ 가 positive definite이기만 하면 복잡한 데이터에 적합하게 무한대로 확장할 수 있다.

그런데 이 함수가 정말 시간적 흐름을 포착할 수 있을까? 보완을 위해 새로운 항을 추가해보자.

[k(x, x\prime) = \sigma_{f}^2 exp(- \frac{( x - x\prime)^2}{2l^2} ) + exp( -2sin^2[\nu \pi (x-x\prime)] ) + \sigma_n^2 \delta(x, x\prime)]

우항의 첫 부분은 마찬가지로 장기간의 트렌드를 포착하기 위해 설계된 부분이고, 두 번째 부분은 빈도를 나타내는 $\nu$ 와 함께 periodicity를 반영하게 된다. 위에서 살펴본 그림의 검은 실선이 위 공분산 함수를 이용하여 적합한 것이다.

4. Fitting Gaussian Process with Python

베이지안 방법론을 위한 대표적인 라이브러리로 PyMC3가 있지만 본 글에서는 scikit-learn 라이브러리를 이용하겠다.

회귀 문제에서는 공분산 함수(kernel)를 명시함으로써 GaussianProcessRegressor를 사용할 수 있다. 이 때 적합은 주변 우도의 로그를 취한 값을 최대화하는 과정을 통해 이루어진다. 이 Class는 평균 함수를 명시할 수 없는데, 왜냐하면 평균 함수는 0으로 고정되어 있기 때문이다.

분류 문제에서는 GaussianProcessClassifier를 사용할 수 있을 것이다. 언뜻 생각하면 범주형 데이터를 적합하기 위해 정규 분포를 사용하는 것이 이상하다. 이는 Latent Gaussian Response Variable을 사용한 뒤 이를 unit interval(다중 분류에서는 simplex interval)로 변환하는 작업을 통해 해결할 수 있다. 이 알고리즘의 결과는 일반적인 머신러닝 알고리즘에 비해 부드럽고 확률론적인 분류 결과를 반환한다. (이에 대한 자세한 내용은 Reference에 있는 2번째 논문의 7페이지를 참조하길 바란다.)

GP의 Posterior는 정규분포가 아니기 때문에 Solution을 찾기 위해 주변 우도를 최대화하는 것이 아니라 Laplace 근사를 이용한다.

이제부터 아주 간단한 예를 통해 라이브러리를 사용하는 법에 대해 소개하겠다. 본 내용은 scikit-learn 라이브러리 홈페이지에서 확인할 수 있다.

아래와 같은 함수를 추정하는 것이 우리의 목표이다.

import numpy as np

# X, y는 학습 데이터

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

X = np.linspace(0.1, 9.9, 20)
X = np.atleast_2d(X).T

# Observations and noise
y = f(X).ravel()
dy = 0.5 + 1.0 * np.random.random(y.shape)
noise = np.random.normal(0, dy)
y += noise

공분산 함수(kernel)를 정의하고 GPR 적합을 시작한다. 본 예시에는 kernel을 구성할 때의 Hyperparameter 최적화에 대한 내용은 포함되어 있지 않다.

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

kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel,
                              n_restarts_optimizer=9,
                              optimizer='fmin_l_bfgs_b',
                              random_state=0)
gp.fit(X, y)

아주 드넓은 공간에서 함수 추정을 해보자.

x = np.atleast_2d(np.linspace(0, 10, 1000)).T
y_pred, sigma = gp.predict(x, return_std=True)

plt.figure()
plt.plot(x, f(x), 'r:', label=r'$f(x) = x\,\sin(x)$')
plt.errorbar(X.ravel(), y, dy, fmt='r.', markersize=10, label='Observations')
plt.plot(x, y_pred, 'b-', label='Prediction')
plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate([y_pred - 1.9600 * sigma,
                        (y_pred + 1.9600 * sigma)[::-1]]),
         alpha=.5, fc='b', ec='None', label='95% confidence interval')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.ylim(-10, 20)
plt.legend(loc='upper left')
plt.show()

Reference

1) GP 논문1
2) GP 논문2
3) GP 설명 블로그
4) scikit-learn 홈페이지
5) PyMC3 홈페이지