EM 알고리즘

EM 알고리즘(expectation-maximization algorithm)은 확률 모델에 관측 불가능한 변수들이 포함되어 있는 경우 최대우도나 최대사후확률 우도를 구하는 방법이다. EM 알고리즘은 기존의 우도를 기반으로 하여 더 좋은 우도를 찾는 계산을 반복하는 구조로 이루어져 있다.

방식

관측할 수 있는 확률변수 $\mathbf{X}$와 관측할 수 없는 확률변수 $\mathbf{Z}$, 그리고 모수 $\boldsymbol{\theta}$가 있을 때, $(\mathbf{X}, \mathbf{Z})$에 대한 확률분포는 $L (\boldsymbol\theta;\mathbf{X},\mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}|\boldsymbol{\theta})$로 주어져 있다. 이 때 최대우도를 계산하고 싶은 경우, 최대화하려는 우도 함수는 $L(\boldsymbol{\theta}; \mathbf{X}) = p(\mathbf{X} |\boldsymbol{\theta}) = \sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z}|\boldsymbol{\theta})$로 정의할 수 있다. 하지만 이 우도함수에는 극대점이 여럿 있을 수 있으며, 그러한 값들을 모두 계산하는 것은 계산적 비용이 높다.

EM 알고리즘은 어떠한 모수 $\boldsymbol{\theta}^{(t)}$를 입력으로 받아서 새로운 모수 $\boldsymbol{\theta}^{(t+1)}$를 찾는다. 이 방법은 크게 다음과 같이 E 단계와 M 단계로 나누어 부른다.

E 단계에서는 $\boldsymbol{\theta}^{(t)}$가 주어졌을 때 새로운 $\boldsymbol{\theta}$를 사용할 때 우도의 기대값 $Q$를 정의한다.

$Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) = \operatorname{E}_{\mathbf{Z}|\mathbf{X},\boldsymbol\theta^{(t)}}\left[ \log L (\boldsymbol\theta;\mathbf{X},\mathbf{Z}) \right] = \sum_{\mathbf{Z}} p(\mathbf{Z}|\mathbf{X}, \boldsymbol{\theta}^{(t)}) \log L(\boldsymbol\theta;\mathbf{X},\mathbf{Z})$

M 단계에서는 $Q$를 최대화하는 새로운 모수 $\boldsymbol\theta^{(t+1)}$를 계산한다.

$\boldsymbol\theta^{(t+1)} = \underset{\boldsymbol\theta}{\operatorname{arg\,max}} \ Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) \,$

실제 EM 알고리즘을 사용할 때에는 모수 $\boldsymbol{\theta}^{(0)}$를 적당한 방법(임의값 등)으로 초기화한 다음, $\boldsymbol{\theta}^{(t)}$를 연속적으로 계산하면서 값이 충분히 수렴될때 멈추는 방식을 사용한다.

EM 알고리즘은 극대점을 찾기 때문에, 구한 값이 실제 최대점은 아닐 수 있다. 이를 피하기 위해 여러 개의 임의 초기값에 대해 계산하는 기법이나 담금질 기법 등이 같이 사용될 수 있다.

예제: 정규분포 혼합 모델 (Gaussian Mixture Model)

두 개의 $d$차원 다변수 정규분포가 혼합되어 있는 경우, 보통 이 모델에서 $d$차원 값을 수집하는 것은 가능하지만 그 값이 두 분포 중 어느 분포에서 나온 것인지는 알 수 없는 경우가 많다. 이 문단에서는 이러한 경우에 대한 EM 알고리즘을 유도한다.

전체 $n$개의 점을 수집하였을 때, $i$번째 수집되는 값은 $(x_i, z_i)$로 나타낼 수 있다. 여기에서 $x_i$는 $d$차원 점을 가리키는 관측 가능한 변수, 그리고 $z_i$는 두 개의 정규분포 중 어느 분포에서 수집되었는지를 가리키는 관측 불가능한 변수이다.

$i$번째 값이 두 분포 중 어느 분포에서 수집되는지에 대한 확률은 $i$와 독립적이기 때문에, $P(z_i)$는 $i$에 관계없는 매개변수로 표현할 수 있다. 이 값을 $P(z_i=1) = \tau_1, P(z_i=2) = \tau_2$인 벡터 $\boldsymbol\tau = (\tau_1, \tau_2)$로 표기하도록 한다($\tau_1 + \tau_2 = 1$).

마지막으로, 두 다변수 정규분포를 각각 $\mathcal{N}_d(\boldsymbol\mu_1, \Sigma_1)$, $\mathcal{N}_d(\boldsymbol\mu_2, \Sigma_2)$로 표기하기로 한다. 그렇다면 이 모델에서 모수는 $\boldsymbol{\theta} = (\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \Sigma_1, \Sigma_2, \boldsymbol\tau)$가 된다.

이제 우도함수 $L(\boldsymbol\theta;\mathbf{x},\mathbf{z})$를 다음과 같이 표현할 수 있다.

$L(\boldsymbol\theta;\mathbf{x},\mathbf{z}) = P(\mathbf{x},\mathbf{z} \vert \boldsymbol\theta) = \prod_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \ \tau_j \ f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j)$

여기에서 $\mathbb{I}$는 표시함수이고, $f$는 다변수 정규분포의 확률 밀도 함수이다. 여기에서 $f$를 풀어서 쓰면 다음과 같이 정리된다.

$L(\boldsymbol\theta;\mathbf{x},\mathbf{z}) = \exp \left\{ \sum_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{n}{2} \log(2\pi) \big] \right\}$

모수 $\boldsymbol\theta^{(t)}$에 대하여 $z_i$에 대한 확률분포 $\operatorname{P}(Z_i=j | X_i=\mathbf{x}_i ;\boldsymbol\theta^{(t)})$는 다음과 같이 계산한다.

$T_{j,i}^{(t)} := \operatorname{P}(Z_i=j | X_i=\mathbf{x}_i ;\boldsymbol\theta^{(t)}) = \frac{\tau_j^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)})}{\tau_1^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_1^{(t)},\Sigma_1^{(t)}) + \tau_2^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_2^{(t)},\Sigma_2^{(t)})}$

그러면 E 단계의 $Q$ 함수는 다음과 같이 계산된다.

$\begin{align}Q(\boldsymbol\theta|\boldsymbol\theta^{(t)}) &= \operatorname{E} [\log L(\boldsymbol\theta;\mathbf{x},\mathbf{Z}) ] \\ &= \sum_{i=1}^n \sum_{j=1}^2 T_{j,i}^{(t)} \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{n}{2} \log(2\pi) \big] \end{align}$

이제 M 단계를 계산한다. $Q$ 함수는 $\tau_j$에 대한 식과 $(\boldsymbol{\mu}_j,\Sigma_j)$에 대한 식의 선형 결합이기 때문에, $\tau_j$와 $(\boldsymbol{\mu}_j,\Sigma_j)$에 대해 독립적으로 최대화를 계산하는 것이 가능하다.

먼저 모수 $\boldsymbol{\tau}^{(t+1)}$는 다음과 같이 유도된다.

$\begin{align}\boldsymbol{\tau}^{(t+1)} &= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}}\ Q(\theta | \theta^{(t)} ) \\ &= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}} \ \left\{ \left[ \sum_{i=1}^n T_{1,i}^{(t)} \right] \log \tau_1 + \left[ \sum_{i=1}^n T_{2,i}^{(t)} \right] \log \tau_2 \right\} \end{align}$

이 식은 이항분포에서의 최대우도를 계산하는 경우와 동일하며, 따라서 다음과 같이 계산할 수 있다.

$\tau^{(t+1)}_j = \frac{\sum_{i=1}^n T_{j,i}^{(t)}}{\sum_{i=1}^n (T_{1,i}^{(t)} + T_{2,i}^{(t)} ) } = \frac{1}{n} \sum_{i=1}^n T_{j,i}^{(t)}$

다음으로 $(\boldsymbol{\mu}_1^{(t+1)}, \Sigma_1^{(t+1)})$에 대한 식은 다음과 같다.

$\begin{align}(\boldsymbol{\mu}_1^{(t+1)},\Sigma_1^{(t+1)}) &= \underset{\boldsymbol{\mu}_1,\Sigma_1} {\operatorname{arg\,max}}\ Q(\theta | \theta^{(t)} ) \\ &= \underset{\boldsymbol{\mu}_1,\Sigma_1} {\operatorname{arg\,max}}\ \sum_{i=1}^n T_{1,i}^{(t)} \left\{ -\tfrac{1}{2} \log |\Sigma_1| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_1)^\top\Sigma_1^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_1) \right\} \end{align}$

이 식은 다변수 정규분포에서의 최대우도를 계산하는 식과 유사하며, 다음과 같이 계산할 수 있다.

$\boldsymbol{\mu}_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{1,i}^{(t)}}$

$\Sigma_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)})^\top }{\sum_{i=1}^n T_{1,i}^{(t)}}$

$(\boldsymbol{\mu}_2^{(t+1)}, \Sigma_2^{(t+1)})$도 같은 방식으로 구할 수 있다.