11. Mixture models and the EM algorithm

11.1. Latent variable models

고차원 결합확률분포의 변수간 의존 관계는 그래프 모델로 모델링할 수 있지만, 관측된 변수들이 숨겨진 공통 원인으로부터 도출됨으로써 상호 연관된다는 가정을 할 수도 있다. 이를 잠재 변수 모델(LVM)이라 한다. 이는 피팅하기 힘들지만, 가측공간에서 상호 연관을 나타내는 매개변수가 더 적은 경우가 많고, 잠재 변수들이 데이터의 표현을 압축하는 병목 역할을 해 줄 수 있다.

11.2. Mixture models

잠재 변수 모델의 가장 간단한 형태는 이산잠재상태 z_{i} \in \{1, \cdots, K\} 가 있을 때 이산적 사전분포 p(z_{i}) = \mathrm{Cat}(\mathbf{\pi}), 가능도 p(\mathbf{x}_{i} | z_{i} = k) = p_{k}(\mathbf{x}_{i})를 두는 것이다. 이 때 p_{k}는 관측의 기반 분포라 한다. 전체 모델은 혼합 모델로서 p(\mathbf{x}_{i} | \mathbf{\theta}) = \sum_{k=1}^{K} \pi_{k} p_{k}(\mathbf{x}_{i} | \mathbf{\theta})의 형태인 혼합가중치를 적용한 볼록합이다.

11.2.1. Mixtures of Gaussians

가장 흔히 사용되는 혼합 모델은 가우시안 혼합, 또는 가우시안 혼합 모델(GMM)로 p(\mathbf{x}_{i} | \mathbf{\theta}) = \sum_{k=1}^{K} \pi_{k} \mathcal{N}(\mathbf{x}_{i} | \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k})이 된다.

가우시안 분포의 혼합. 각 개별 분포의 컨투어, 전체 표면 분포.

11.2.2. Mixture of multinoullis

데이터가 D차원 비트 벡터로 이루어져 있으면 베르누이 분포의 곱으로 조건부 밀도함수를 나타낼 수 있으며, 그 형태는 p(\mathbf{x}_{i} | z_{i} = k, \mathbf{\theta}) = \prod_{j=1}^{D} \mu_{jk}^{x_{ij}}(1 - \mu_{jk})^{1 - x_{ij}}이 된다. 이 때 \mathbf{\Sigma}_{k} = \mathrm{diag}(\mu_{jk}(1-\mu_{jk}))로 뒀을 때 \mathbb{E}[\mathbf{x}] = \sum_{k} \pi_{k} \mathbf{\mu}_{k}, \mathrm{cov}[\mathbf{x}] = \sum_{k} \pi_{k} [\mathbf{\Sigma}_{k} + \mathbf{\mu}_{k}\mathbf{\mu}_{k}^{T}] - \mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{T}이 됨을 알 수 있는데, 이로부터 개별 분포가 인수분해 가능하더라도 결합 분포는 그렇지 않을 수 있다는 것을 알 수 있다. 따라서 혼합 분포는 변수간의 상관 관계도 모델링할 수 있다.

11.2.3. Using mixture models for clustering

혼합 모델의 첫 번째 용도는 이를 블랙박스 분포 모델로 사용하는 것이다. 두 번째 용도는 클러스터링이다. 이는 우선 혼합 모델을 피팅한 뒤 각각 클래스에 속할 책임도를 나타내는 r_{ik} = p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta}) = \frac{p(z_{i} = k | \mathbf{\theta}) p(\mathbf{x}_{i} | z_{i} = k, \mathbf{\theta})}{\sum_{k' = 1}^{K}  p(z_{i} = k' | \mathbf{\theta}) p(\mathbf{x}_{i} | z_{i} = k', \mathbf{\theta}) }을 계산하는 것이다. 이를 소프트 클러스터링이라 하며, 생성적 분류기를 이용해 계산하는 것과 같다. 차이가 있다면 학습 시점에서 혼합 모델은 z_{i}를 관측할 수 없지만 생성적 분류기에서는 y_{i}를 관측하게 된다는 것이 차이이다.

클러스터 할당에 대한 불확실성을 1 - \max_{k} r_{ik}라 하면 최대사후확률추정 z_{i}^{\ast} = \mathrm{argmax}_{k}r_{ik} = \mathrm{argmax}_{k} (\log p(\mathbf{x}_{i} | z_{i} = k, \mathbf{\theta}) + \log p(\mathbf{z}_{i} = k | \mathbf{\theta}))을 계산하는 하드 클러스터링도 가능하다. 이 때 클러스터를 이름짓는 순서는 상관이 없다(라벨 스위칭). 다른 예는 유전자 타입을 무게중심이나 프로토타입을 중심으로 클러스터링할 수 있다.

K-평균 알고리즘을 통한 효모의 유전자형 클러스터링.

다른 예는 MNIST 데이터셋에 대한 클러스터링이다. 아래의 분류 결과에는 오류가 있는데 원인을 추측해보자면

  • 모델이 너무 간단하거나
  • 어떤 숫자의 경우 쓰는 방법이 여러 가지라서 클러스터의 배분이 적절하지 않거나
  • 가능도함수가 볼록하지 않거나

등등이 있을 수 있다.

이진 MNIST 데이터의 분류 예.

11.2.4. Mixtures of experts

혼합 모델로 분류와 회귀에 쓰이는 구별적 모델도 만들 수 있다. 예를 들어 여러 개의 다른 선형 회귀 모델을 가중치합한다고 하면 p(y_{i} | \mathbf{x}_{i}, z_{i} = k, \mathbf{\theta}) = \mathcal{N}(y_{i} | \mathbf{w}_{k}^{T} \mathbf{x}_{i}, \sigma_{k}^{2}), p(z_{i} | \mathbf{x}_{i}, \mathbf{\theta}) = \mathrm{Cat}(z_{i} | \mathcal{S}(\mathbf{V}^{T} \mathbf{x}_{i})) 로 나타낼 수 있는데, 이런 형태의 모델을 전문가의 혼합이라고 한다. 이 때 p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta})는 어떤 전문가 모델을 선택할지를 결정하는 개통 함수가 된다.

3개의 선형 회귀 전문가 모델의 혼합 예.

각각의 전문가 모델에는 어떤 모델을 끼워넣든 상관이 없다. 신경망을 끼워넣으면 혼합 신경망이 된다. 각각의 전문가 모델을 혼합 전문가 모델로 하는 것도 가능하다. 이를 계층 전문가 혼합 모델이라 한다.

11.2.4.1. Application to inverse problems

전문가 혼합 모델은 역함수가 유일하지 않을 때 이를 구하는 역함수 문제에 유용하다. 다른 예는 전진 모델을 통해 영상에서 사람의 움직임 추적을 하는 데 쓸 수도 있다. 이 때에는 선형 전문가 모델을 사용하는데, 사후평균이 아닌 사후최빈값을 사용하는 것이 더 좋다.

11.3. Parameter estimation for mixture models

혼합 모델의 매개변수 추정은 어떻게 하는가? 데이터가 완전할 경우에는 사후분포가 인수분해되므로 게산이 간단하지만 은닉 변수가 있거나 데이터가 완전하지 못하면 이렇게 할 수가 없다. 변수간 조건부독립이 성립하지 않게 되기 때문이다.

11.3.1. Unidentifiability

잠재 변수 모델에서 p(\mathbf{\theta} | \mathcal{D})을 계산하는 데 있어 가장 큰 어려운 점은 사후분포의 최빈값이 유일하지 않을 수 있다는 것이다. 은닉 변수가 존재할 때 각각의 은닉 변수를 채워넣는 데에는 서로 다른 단일 최빈 가능도함수가 존재하므로 은닉 변수를 적분해 빼내면 사후분포는 최빈값이 여러 개가 되므로 최대가능도추정이 유일하지 않아 매개변수가 식별 가능하지 않다. 문제는 최빈값의 갯수도 확실하게 답할 수 없다는 것이다. 최대 K!개일 수 있지만 이 중 몇 개의 최고점은 합쳐질 수 있다. 가우시안 혼합 모델의 경우 최대가능도추정의 정확한 분포를 찾는 것은 NP-난해 문제이다.

최고점이 대칭을 이뤄서 매개변수가 식별 불가능해지는 혼합 모델 분포.

이런 불식별성은 베이지안 추론에 있어 문제가 된다. 일반적 접근법은 하나의 국소적 최빈값을 택하는 것이다. 표본 크기가 클 경우 이는 꽤 잘 먹히는데 매개변수의 수가 표본수보다 훨씬 작아지므로 매개변수의 사후 불확실성이 잠재변수의 사후 불확실성보다 훨씬 작아지기 때문이다. 그러므로 p(\mathbf{\theta} | \mathcal{D}) 대신 p(z_{i} | \mathbf{x}_{i}, \mathbf{\theta})을 계산해도 크게 무리가 없다.

11.3.2. Computing a MAP estimate is non-convex

게다가 최대사후확률추정은 볼록함수가 아니다. 결합확률분포 p(\mathbf{z}, \mathbf{x} | \mathbf{\theta}) = \frac{1}{Z(\mathbf{\theta})} e^{\mathbf{\theta}^{T} \mathbf{\phi}(\mathbf{z}, \mathbf{x})}이 지수족이라면 완전 데이터 로그가능도는 다음과 같다.

l_{c}(\mathbf{\theta}) = \mathbf{\theta}^{T}(\sum_{i} \mathbf{\phi}(\mathbf{x}_{i}, \mathbf{z}_{i})) - NZ(\mathbf{\theta})

이는 오목함수이므로 유일한 최대값을 가지지만, 데이터가 완전하지 않을 때에는 관측 데이터 로그가능도는 다음과 같아진다.

l(\mathbf{\theta}) = \sum_{i} \log [\sum_{\mathbf{z}_{i}}e^{ \mathbf{\theta}^{T}(\sum_{i} \mathbf{\phi}(\mathbf{x}_{i}, \mathbf{z}_{i})) }] - N \log Z(\mathbf{\theta})

이 함수는 일반적으로 볼록/오목함수가 아니므로 국소적 최적점만을 갖는다. 이에 대처하기 위해 보통은 국소적 최적점을 찾되 복수의 무작위 재시작을 통해 보다 좋은 국소적 최적점을 찾는다. 가우시안의 혼합일 경우에는 볼록함수적 해법이 존재하는데, 이는 희소 커널 로지스틱 회귀의 비지도학습 버전이라고 봐도 되지만 이 때 쓰이는 l_{1} 패널티 함수가 희소성을 나타내는 데 있어 좋은 방법이 아니기 때문에 만능이라고 할 수는 없다.

11.4. The EM algorithm

데이터 일부가 유실되었거나 잠재 변수가 존재한다면 최대가능도/최대사후확률추정은 어려워진다. 이 때에는 음의 로그 가능도(NLL) \mathrm{NLL}(\mathbf{\theta}) = -\frac{1}{N} \log p(\mathcal{D} | \mathbf{\theta})을 경사 하강법으로 국소적 최소값을 구하는 것이 하나의 방법이 될 수 있지만, 대개의 경우엔 더 간단한 기대값 최대화(EM) 알고리즘을 사용하는 것이 방법이 될 수 있다. 이는 데이터가 완전히 관측된 경우엔 최대가능도/최대사후확률추정이 계산이 쉽다는 점을 이용하는데, E 단계에서는 매개변수가 주어진 상황에서 유실된 변수를 추론하고 M 단계에서는 이 추론한 변수들로 채워넣은 데이터를 기반으로 매개변수를 최적화한다. EM 알고리즘은 E 단계와 M 단계를 번갈아 반복하는 알고리즘이다.

11.4.1. Basic idea

\mathbf{x}_{i}가 관측된 변수고, \mathbf{z}_{i}가 은닉 변수일 경우엔 목표는 관측된 데이터의 로그 가능도인 l(\mathbf{\theta}) = \sum_{i=1}^{N} \log [\sum_{\mathbf{z}_{i}} p(\mathbf{x}_{i}, \mathbf{z}_{i} | \mathbf{\theta})]를 최대화하는 것이다. 합의 로그를 계산하는 것은 쉽지 않기 때문에 최적화하기 어렵다는 문제가 있다.

이를 위해 EM 알고리즘은 완전 데이터 로그 가능도 l_{c}(\mathbf{\theta}) = \sum_{i=1}^{N} \log p(\mathbf{x}_{i}, \mathbf{z}_{i} | \mathbf{\theta})를 정의해 보조 함수 Q로 정의되는 기대 완전 데이터 로그 가능도 Q(\mathbf{\theta}, \mathbf{\theta}^{t-1}) = \mathbb{E}[l_{c}(\mathbf{\theta}) | \mathcal{D}, \mathbf{\theta}^{t-1}]를 계산한다. 이것이 기대 충족 통계량을 계산하는 E 단계이다.

M 단계에서는 Q 함수를 \mathbf{\theta}에 대해 최적화하여 \mathbf{\theta}^{t} = \mathrm{argmax}_{\mathbf{\theta}} Q(\mathbf{\theta}, \mathbf{\theta}_{t-1}) 을 계산한다. 최대사후확률추정의 경우에는 식에 \log p(\mathbf{\theta})의 항이 붙고, E 단계에서는 차이가 없다.

EM 알고리즘은 관측 데이터 로그 가능도를 단조증가시킨다.

11.4.2. EM for GMMs

11.4.2.1. Auxillary function

Q(\mathbf{\theta}, \mathbf{\theta}^{t-1}) = \mathbb{E}[l_{c}(\mathbf{\theta}) | \mathcal{D}, \mathbf{\theta}_{t-1}] = \sum_{i} \mathbb{E}[\log [\prod_{k=1}^{K} (\pi_{k} p(\mathbf{x}_{i} | \mathbf{\theta}_{k})^{\mathbf{1}_{\mathbf{z}_{i} = k}})]]   = \sum_{i} \sum_{k} \mathbb{E}[\mathbf{1}_{\mathbf{z}_{i} = k}] \log [\pi_{k} p(\mathbf{x}_{i} | \mathbf{\theta}_{k})] = \sum_{i} \sum_{k} p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta}^{t-1}) \log [\pi_{k} p(\mathbf{x}_{i} | \mathbf{\theta}_{k}]

r_{ik} = p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta}_{(t-1)})은 i번째 데이터에 대한 k번째 클러스터의 책임도라 불린다. 이는 E 단계에서 계산된다.

11.4.2.2. E step

r_{ik} = \frac{\pi_{k} p(\mathbf{x}_{i} | \mathbf{\theta}_{k, (t-1)}}{\sum_{k'} \pi_{k'} p(\mathbf{x}_{i} | \mathbf{\theta}_{k', (t-1)}} 이며, 이는 임의의 혼합 모델에 대해 성립한다.

11.4.2.3. M step

\hat{\pi}_{k} = \frac{\sum_{i} r_{ik}}{N} = \frac{r_{k}}{N} 이다.

l(\mathbf{\mu}_{k}, \mathbf{\Sigma}_{k}) = \sum_{k} \sum_{i} r_{ik} \log p(\mathbf{x}_{i} | \mathbf{\theta}_{k}) = -\frac{1}{2} \sum_{i} r_{ik} [\log \lvert \mathbf{\Sigma}_{k} \rvert + (\mathbf{x}_{i} - \mathbf{\mu}_{k})^{T}\mathbf{\Sigma}_{k}^{-1}(\mathbf{x}_{i} - \mathbf{\mu}_{k})]이므로,

\hat{\mathbf{\mu}}_{k} = \frac{\sum_{i} r_{ik} \mathbf{x}_{i}}{r_{k}}

\hat{\mathbf{\Sigma}}_{k} = \frac{\sum_{i} r_{ik} \mathbf{x}_{i}\mathbf{x}_{i}^{T}}{r_{k}} - \mathbf{\mu}_{k}\mathbf{\mu}_{k}^{T} 가 된다.

이제 새 매개변수 추정값으로 다시 E 단계를 시작한다.

11.4.2.4. Example

EM 알고리즘의 용례.

11.4.2.5. K-means algorithm

가우시안 혼합 모델에 대한 EM 알고리즘의 변형으로는 K-평균 알고리즘이 있다. 가우시안 혼합 모델에 대해 \mathbf{\Sigma}_{k} = \sigma^{2} \mathbf{I}_{D}, \pi_{k} = 1/K 이 고정이라 하자. E 단계에서 z_{i}^{\ast} = \mathrm{argmax}_{k} p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta})에 대해 사후분포가 p(z_{i} = k | \mathbf{x}_{i}, \mathbf{\theta}) \simeq \mathbf{1}_{k = z_{i}^{\ast}}로 근사된다고 하자. (즉, 사후확률을 최대로 만드는 클러스터로 강제 지정). 이를 강한 EM이라 부른다. 공분산이 전부 공유되며 구형이므로, 사후확률을 최대로 만드는 클러스터는 z_{i}^{\ast} = \mathrm{argmin}_{k} \lVert \mathbf{x}_{i} - \mathbf{\mu}_{k} \rVert_{2}^{2} 로 구해진다. E 단계마다 이는 기본적으로 O(NKD)의 시간복잡도를 갖지만 이는 더 최적화할 수가 있다.

M 단계에서는 각각의 클러스터 중심이 \mathbf{\mu}_{k} = \frac{1}{N_{k}} \sum_{i : z_{i} = k} \mathbf{x}_{i} 로 업데이트된다.

11.4.2.6. Vector quantization

실변수 벡터 \mathbf{x}_{i}를 K개의 코드북 벡터 \mathbf{\mu}_{k}로 분류하는 것을 벡터 양자화(VQ)라 한다. K-평균 알고리즘과 똑같이 \mathrm{encode}(\mathbf{x}_{i}) = \mathrm{argmin}_{k} \lVert \mathbf{x}_{i} - \mathbf{\mu}_{k} \rVert_{2}^{2} 로 정의하고 \mathrm{decode}(k) = \mathbf{\mu}_{k}로 정의하면 코드북의 정확도는 재구성 오차 또는 왜곡J(\mathbf{\mu}, \mathbf{z} | K, \mathbf{X}) = \frac{1}{N} \sum_{i=1}^{N} \lVert \mathbf{x}_{i} - \mathbf{\mu}_{z_{i}} \rVert_{2}^{2}가 되며 양자화는 이것을 최소화시키는 것이다. 모든 벡터에 대해 코드북을 할당하면 왜곡은 0이 되겠지만 그것은 메모리의 문제 때문에 불가능하며 의미도 없다. 그 대신의 각각의 벡터를 한 번 저장한 뒤 그에 대한 코드워드를 포인터를 통해 할당할 수 있다. 이렇게 하면 인코딩의 비율O(\lg K)가 되어 기존의 O(DC)에 비해 대폭 줄어든다. 용례는 이미지 압축 등이 있다.

벡터 양자화를 이용한 이미지 압축의 예.

11.4.2.7. Initialization and avoiding local minima

K-평균을 초기화할 때 다음과 같은 응용이 가능하다: 첫 점은 무작위로 잡고, 이후의 점들은 지금까지 잡은 마지막 점과의 거리의 제곱에 비례하도록 잡는 것이다. 이런 가장 먼 점 클러스터링K-평균++이라 한다. 이는 왜곡값이 최적값보다 O(\log K) 이상 크게 되지 않음을 보장한다.

클러스터 개수를 지정하는 휴리스틱으로서는 혼합 가중치에 기반해 최초에 클러스터에 점수를 할당한 뒤, 알고리즘의 각 반복마다 최고 점수 클러스터를 쪼개는 것이다. 이를 목표 클러스터 개수에 도달할 때까지 반복한다.

11.4.2.8. MAP estimation

최대가능도추정은 과적합될 수 있으며 이는 가우시안 혼합 모델에서 더 심각해진다. 분산이 적어질 수록 가능도가 한 점에 몰리기 때문이다. 이를 피해서 최대사후확률추정을 사용해 Q'(\mathbf{\theta}, \mathbf{\theta}^{\mathrm{old}}) = Q(\mathbf{\theta}, \mathbf{\theta}^{\mathrm{old}}) + \log p(\mathbf{\pi}) + \sum_{k} \log p(\mathbf{\theta}_{k}) 을 대신 최적화할 수 있다. E 단계는 바뀌지 않고 M 단계만 바뀜을 유의하라.

혼합 가중치에 대해서는 디리클레 사전분포를 사용하면 최대사후확률추정은 \hat{\pi}_{k} = \frac{r_{k} + \alpha_{k} - 1}{N + \sum_{k} \alpha_{k} - K}이 된다. 균등사전분포를 사용하면 최대가능도추정과 같다.

클래스 조건분포의 매개변수 \mathbf{\theta}_{k}에 대해서는 조건분포의 형태에 따라 최대사후확률추정의 형태도 달라진다. 간단한 예로는 NIW 사전분포 p(\mathbf{\mu}_{k}, \mathbf{\Sigma}_{k}) = \mathrm{NIW}(\mathbf{\mu}_{k}, \mathbf{\Sigma}_{k} | \mathbf{m}_{0}, \kappa_{0}, \nu_{0}, \mathbf{S}_{0})를 사용했을 때, \bar{\mathbf{x}}_{k} = \frac{\sum_{i} r_{ik} \mathbf{x}_{i}}{r_{k}}, $\latex \mathbf{S}_{k} = \sum_{i} r_{ik} (\mathbf{x}_{i} – \bar{\mathbf{x}}_{k})(\mathbf{x}_{i} – \bar{\mathbf{x}}_{k})^{T} &s=1$이라 하면 최대사후확률추정은 다음과 같아진다.

\hat{\mathbf{\mu}}_{k} = \frac{r_{k} \bar{\mathbf{x}}_{k} + \kappa_{0} \mathbf{m}_{0}   }{r_{k} + \kappa_{0}}

\hat{\mathbf{\Sigma}}_{k} = \frac{\mathbf{S}_{0} + \mathbf{S}_{k} + \frac{\kappa_{0}r_{k}}{\kappa_{0} + r_{k}}(\bar{\mathbf{x}}_{k} - \mathbf{m}_{0})(\bar{\mathbf{x}}_{k} - \mathbf{m}_{0})^{T}}{\nu_{0} + r_{k} + D + 2}

이 때 \mathbf{S}_{0} = \frac{1}{K^{\frac{1}{D}}} \mathrm{diag}(s_{1}^{2}, \cdots, s_{D}^{2})을 많이 사용한다.

가우시안 혼합 모델의 특이점 발생의 예.
가우시안 혼합 모델에서 최대사후확률추정과 최대가능도추정의 비교.

11.4.3. EM for mixture of experts

전문가 혼합 모델에 대한 EM 알고리즘은 간단히 정의할 수 있다. 기대완전데이터 로그가능도는 \pi_{ik} = \mathcal{S}(\mathbf{V}^{T} \mathbf{x}_{i})_{k}라 할 때 Q(\mathbf{\theta}, \mathbf{\theta}^{\mathrm{old}}) = \sum_{i=1}^{N} \sum_{k=1}^{K} r_{ik} \log[\pi_{ik} \mathcal{N}(y_{i} | \mathbf{w}_{k}^{T} \mathbf{x}_{i}, \sigma_{k}^{2})]이 되고, 이 때 r_{ik} \propto \pi_{ik}^{\mathrm{old}} \mathcal{N}(y_{i} | (\mathbf{w}_{k}^{\mathrm{old}})^{T} \mathbf{x}_{i}, (\sigma_{k}^{\mathrm{old}})^{2})이 된다. E 단계는 일반적인 혼합 모델과 똑같이 계산한다.

M 단계에서는, 선형 회귀 인자에 대해 목적 함수가 Q(\mathbf{\theta}_{k}, \mathbf{\theta}^{\mathrm{old}}) = \sum_{i=1}^{N} r_{ik} [-\frac{1}{\sigma_{k}^{2}} (y_{i} - \mathbf{w}_{k}^{T} \mathbf{x}_{i})]이 된다. 이는 가중치 최소 제곱 문제와 같다. 최대가능도추정은 다음과 같다:

\hat{\mathbf{w}}_{k} = (\mathbf{X}^{T} \mathbf{R}_{k} \mathbf{X})^{-1} \mathbf{X}^{T} \mathbf{R}_{k} \mathbf{y}

\hat{\sigma}_{k}^{2} = \frac{\sum_{i=1}^{N} r_{ik}(y_{i} - \mathbf{w}_{k}^{T} \mathbf{x}_{i})^{2}}{\sum_{i=1}^{N} r_{ik}}

혼합 가중치에 대한 추정은 개통 함수 매개변수 \mathbf{V}에 대한 추정으로 대체할 수 있다. 목적 함수 l(\mathbf{V}) = \sum_{i} \sum_{k} r_{ik} \log \pi_{ik}이 되는데, 이는 다항 로지스틱 회귀를 최적화할 때와 같은 방식으로 최적화하면 된다.

11.4.4. EM for DGMs with hidden variables

위의 아이디어를 은닉 변수가 있는 방향그래프 모델에 대한 EM 알고리즘으로 확장할 수 있다. E 단계에서는 은닉 변수를 추정하고 M 단계에서는 이 변수들로 최대가능도추정을 하는 것이다.

계산의 편의를 위해 조건부 확률분포가 이산적인 경우를 살펴보자.

p(x_{it} | \mathbf{x}_{i, \mathrm{pa}(t)}, \mathbf{\theta}_{t}) = \prod_{c=1}^{K_{\mathrm{pa}(t)}} \prod_{k=1}^{K_{t}} \theta_{tck}^{\mathbf{1}_{x_{it} = 1, \mathbf{x}_{i, \mathrm{pa}(t)} = c}} 이 된다. 그러므로 분류에 대한 실측 개수 \bar{N}_{tck} = \sum_{i} p(x_{it} = k, \mathbf{x}_{i, \mathrm{pa}(t)} = c | \mathcal{D}_{i})에 대해 기대완전데이터 로그가능도는 \mathbb{E}[\log p(\mathcal{D} | \mathbf{\theta})] = \sum_{t} \sum_{c} \sum_{k} \bar{N}_{tck} \log \theta_{tck}이 된다.

주변값 p(x_{it}, \mathbf{x}_{i, \mathrm{pa}(t)} | \mathcal{D}_{i}, \mathbf{\theta})은 그래픽 모델 추론 알고리즘으로 계산될 수 있다. 이가 주어졌을 때 M 단계에 의한 매개변수의 추정값은 \hat{\theta}_{tck} = \frac{\bar{N}_{tck}}{\sum_{k'}\bar{N}_{tck'}}이 된다.

11.4.5. EM for the Student distribution

가우시안 분포는 이상치에 취약하므로 스튜던트 t 분포를 쓸 수 있다. 그러나 최대가능도추정이 닫힌 형태로 주어지지 않기 때문에 반복적 최적화 알고리즘을 써야 한다. 이를 위해서 스튜던트 분포를 다음의 가우시안 비례 혼합 형태인 \mathcal{T}(\mathbf{x}_{i} | \mathbf{\mu}, \mathbf{\Sigma}, \nu) = \int \mathcal{N}(\mathbf{x}_{i} | \mathbf{\mu}, \mathbf{\Sigma}/z_{i}) \mathrm{Ga}(z_{i} | \frac{\nu}{2}, \frac{\nu}{2}) d z_{i}로 나타낼 수 있다는 특성을 이용한다.

이를 이용하면 마할라노비스 거리 \delta_{i} = (\mathbf{x}_{i} - \mathbf{\mu})^{T} \mathbf{\Sigma}^{-1}(\mathbf{x}_{i} - \mathbf{\mu}) 이라 할 때 완전데이터 로그가능도는 l_{c}(\mathbf{\theta}) = -\frac{1}{2} N \log \lvert \mathbf{\Sigma} \rvert - \frac{1}{2} \sum_{i=1}^{N} z_{i} \delta_{i} - N \log \Gamma(\frac{\nu}{2}) + \frac{1}{2} N \nu \log \frac{\nu}{2} + \frac{\nu}{2} \sum_{i=1}^{N} (\log z_{i} - z_{i})이 된다.

11.4.5.1. EM with \nu known

\nu를 알고 있다면 E 단계는 다음과 같다.

\bar{z}_{i, t} = \frac{\nu_{t} + D}{\nu_{t} + \delta_{i, t}}

M 단계는 다음과 같다.

\hat{\mathbf{\mu}}_{t+1} = \frac{\sum_{i} \bar{z}_{i, t} \mathbf{x}_{i}}{\sum_{i} \bar{z}_{i, t}}

\hat{\mathbf{\Sigma}}_{t+1} = \frac{1}{N} \sum_{i} \bar{z}_{i, t} (\mathbf{x}_{i} - \hat{\mathbf{\mu}}_{t+1})(\mathbf{x}_{i} - \hat{\mathbf{\mu}}_{t+1})^{T}

11.4.5.2. EM with \nu unknown

\Psi(x)가 디감마 함수일 때 \bar{l}_{i, t} = \log(\bar{z}_{i, t}) + \Psi(\frac{\nu_{t} + D}{2}) - \log(\frac{\nu_{t} + D}{2}) 이 된다. 이는 1D 제한 최적화 문제로 유일한 해를 갖는다. M 단계에서 경사 하강법에 기반한 최적화를 하는 것은 일반화된 EM 알고리즘의 예시이다. EM 알고리즘은 M 단계에서 매개변수에 부분적 개선을 시행하더라도 국소적 최적값으로 수렴한다.

11.4.5.3. Mixtures of Student distributions

위의 방법은 스튜던트 분포의 혼합을 피팅하는 데 쉽게 확장할 수 있다.

은행 파산에 관련된 데이터셋을 스튜던트 혼합 분포/가우시안 혼합 분포로 피팅했을 때. 스튜던트 혼합 분포가 더 잘 맞는 것을 보여준다.

11.4.6. EM for probit regression

프로빗 회귀에 EM 알고리즘을 적용해 보자. 사전분포 p(\mathbf{w}) = \mathcal{N}(\mathbf{w} | \mathbf{0}, \mathbf{V}_{0})에 대해 완전 데이터 로그가능도는 l(\mathbf{z}, \mathbf{w} | \mathbf{V}_{0}) \propto \sum_{i} \log p(y_{i} | z_{i}) - \frac{1}{2} (\mathbf{z} - \mathbf{X}\mathbf{w})^{T} (\mathbf{z} - \mathbf{X}\mathbf{w}) - \frac{1}{2} \mathbf{w}^{T} \mathbf{V}_{0}^{-1} \mathbf{w} 이 된다.

E 단계에서의 사후 분포는 \mathrm{sg}(y_{i}) = 1 if y_{i} = 1, -1 if y_{i} = 0로 정의할 때 절사 가우시안 p(z_{i} | y_{i}, \mathbf{x}_{i}, \mathbf{w}) = \mathcal{N}(z_{i} | \mathbf{w}^{T} \mathbf{x}_{i}, 1) \mathbf{1}_{z_{i} = \mathrm{sg}(y_{i})}이다. 이의 기대값은 \mathbb{E}[z_{i} | \mathbf{w}, \mathbf{x}_{i}] = \mu_{i} + \mathrm{sg}{y_{i}} \frac{\pi(\mu_{i})}{(1-y_{i}) + \mathrm{sg} \Pi(\mu_{i})}이다 (\mu_{i}  = \mathbf{w}^{T} \mathbf{x}_{i}).

M 단계에서는 \mathbf{w}가 능선 회귀로 추정되어 \hat{\mathbf{w}} = (\mathbf{V}_{0}^{-1} + \mathbf{X}^{T} \mathbf{X})^{-1} \mathbf{X}^{T} \mathbf{\mu} 이 된다.

프로빗 회귀에 대한 EM 알고리즘은 간단하지만 직접적 경사 하강법에 비해 수렴 속도가 훨씬 느릴 수 있다. 이는 E 단계에서 z의 부호만 관찰할 수 있기 때문에 사후 엔트로피가 높을 수 있기 때문이다. 정규화 계수를 강하게 걸거나 데이터 증강을 하면 수렴 속도를 빠르게 할 수 있다.

프로빗 회귀에서 EM 알고리즘과 유사뉴턴법간 수렴속도 비교. 유사뉴턴법이 월등히 빠른 속도를 보여주나, 수렴값의 성능은 떨어진다.

11.4.7. Theoretical basis for EM

EM 알고리즘이 반복시마다 관측 데이터 로그가능도를 단조증가시킴을 증명한다.

11.4.7.1. Expected complete data log likelihood is a lower bound

젠센 부등식에 의해 관측 데이터 로그가능도는 다음의 하한을 갖는다.

l(\mathbf{\theta}) \geq \sum_{i} \sum_{\mathbf{z}_{i}} q_{i}(\mathbf{z}_{i}) \log \frac{p(\mathbf{x}_{i}, \mathbf{z}_{i} | \mathbf{\theta})}{q_{i}(\mathbf{z}_{i})} = \sum_{i} \mathbb{E}_{q_{i}} [\log p(\mathbf{x}_{i}, \mathbf{z}_{i} | \mathbf{\theta})] + \mathbb{H}(q_{i}) =

이는 임의의 분포 q_{i}에 대해 성립하기 때문에 가장 강한 하한을 갖도록 q를 선택하면 된다.

\sum_{\mathbf{z}_{i}} q_{i}(\mathbf{z}_{i}) \log \frac{p(\mathbf{x}_{i}, \mathbf{z}_{i} | \mathbf{\theta})}{q_{i}(\mathbf{z}_{i})} = -\mathbb{KL}(q_{i}(\mathbf{z}_{i}) \lVert p(\mathbf{z}_{i} | \mathbf{x}_{i}, \mathbf{\theta})) + \log p(\mathbf{x}_{i} | \mathbf{\theta}) 이므로, 이를 최대화시키면 \log p(\mathbf{x}_{i} | \mathbf{\theta})이 된다. 이것이 최대화되는 조건은 KL 발산이 0인 경우이므로 q_{i, t}(\mathbf{z}_{i}) = p(\mathbf{z}_{i} | \mathbf{x}_{i}, \mathbf{\theta}_{t})인 경우이다. 이를 보조 함수에 대입하면 완전 데이터 로그가능도를 얻는다.

즉 EM 알고리즘에서 E 단계에서 얻는 관측 데이터 로그가능도에 대한 하한은 가장 강한 하한이다. 이에 의해 M 단계에서의 매개변수 업데이트는 관측 데이터 로그가능도를 증가시킴이 보장된다.

유계 최적화 문제로서의 EM 알고리즘.

11.4.7.2. EM monotonically increases the observed data log likelihood

l(\mathbf{\theta}_{t+1}) \geq Q(\mathbf{\theta}_{t+1}, \mathbf{\theta}_{t}) \geq Q(\mathbf{\theta}_{t}, \mathbf{\theta}_{t}) = l(\mathbf{\theta}_{t})이므로 EM 알고리즘은 관측 데이터 로그가능도를 단조증가시킨다.

11.4.8. Online EM

데이터가 스트리밍될 때에는 온라인 학습이 필요하며 EM 알고리즘에도 예외가 아니다. 이를 온라인 EM이라 하며 첫째로는 가능도의 하한에 대해 q_{i}를 한 번에 하나씩 최적화시키는 증가 EM이 있고 둘째로는 추계적 근사 이론을 적용하는 계단식 EM이 있다.

11.4.8.1. Batch EM review

먼저 배치 EM 알고리즘을 리뷰해 보자. 각각의 데이터 배치에 대한 충족 통계량 벡터를 \mathbf{\phi}(\mathbf{x}_{i}, \mathbf{z})이라 할 때 배치 EM 알고리즘은 매개변수의 평균값 (초기값은 \mathbf{0})을 배치 각각에 대해 \mathbf{s}_{i} =  \sum_{\mathbf{z}} p(\mathbf{z} | \mathbf{x}_{i}, \mathbf{\theta}(\mathbf{\mu})) \mathbf{\phi}(\mathbf{x}_{i}, \mathbf{z})만큼 업데이트시킨다. 이를 수렴할 때까지 반복하는 것이 배치 EM 알고리즘이다.

11.4.8.2. Incremental EM

증가 EM에서는 \mathbf{s}_{i, \mathrm{new}}를 계산하는 것은 같으나 무작위로 선택된 샘플 각각에 대해 매개변수 평균값 (초기값은 \sum_{i} \mathbf{s}_{i}를 기존 값과의 차 \mathbf{s}_{i, \mathrm{new}} - \mathbf{s}_{i, \mathrm{old}}만큼 더한 뒤 뒤 \mathbf{s}_{i}를 재차 업데이트한다. 이를 수렴할 때까지 반복하는 것이 증가 EM 알고리즘이다.

11.4.8.3. Stepwise EM

계단식 EM에서는 \mathbf{s}_{i}를 계산할 때 무작위로 선택된 샘플 각각에 대해 \mathbf{\mu}를 어떤 가중치 \eta_{k}만큼 이동시켜 \mathbf{\mu} = (1 - \eta_{k}) \mathbf{\mu} + \eta_{k} \mathbf{s}_{i}의 형태로 업데이트한다. 이를 수렴할 때까지 반복하는 것이 계단식 EM 알고리즘이다.

11.4.9. Other EM variants

  • 담금질된 EM결정론적 담금질이라는 방법을 사용하는데, 사후분포를 특정 지점까지 상승시킨 뒤 점진적으로 하강시켜 전역 최대값을 찾는 방법이다.
  • 변량 EM은 E 단계에서 추정 계산을 닫힌 형태로 할 수 없을 때 변량 추론 방법으로 근사하는 것이다.
  • 몬테 카를로 EM(MCEM)은 은닉 변수에 대해 사후분포에서 몬테 카를로 샘플링을 한 뒤 그를 이용해 충족 통계랑을 계산한 뒤 평균내는 것이다. 샘플링을 1개만 할 경우 추계적 EM이라 한다. 몬테 카를로 샘플링을 할 때 추계적 근사를 해서 수렴 속도를 빠르게 하는 방법을 추계적 근사 EM이라 한다.
  • 일반화된 EM은 E 단계는 정확히 계산할 수 있으나 M 단계는 정확히 계산할 수 없을 때 사용하며, 경사 하강법을 반복해 부분적 M 단계를 수행해 기대완전데이터 로그가능도를 증가시킨다.
  • 기대 조건분포 최대화(ECM 알고리즘)는 매개변수간 의존 관계가 존재할 때 M 단계를 순차적으로 진행하는 것이다. ECME 알고리즘은 E 단계를 무시하고 M 단계의 관측 데이터 로그가능도를 직접 최대화하는 것이다.
  • 풀린 EM은 유실된 데이터가 많을 때 EM 알고리즘의 수렴 속도가 느리다는 맹점을 보완하며 \mathbf{\theta}_{t+1} = \mathbf{\theta}_{t} + \eta(M(\mathbf{\theta}_{t}) - \mathbf{\theta}_{t})의 업데이트를 수행한다. \eta = 1이면 이것은 일반적인 EM이 되지만 더 큰 값을 쓸 경우 수렴 속도가 훨씬 빨라진다. 하지만 아예 수렴을 실패하게 할 수도 있다.

사실 EM은 한계 최적화 내지는 부차화-최적화(MM) 알고리즘들 중 하나일 뿐이다.

변량 EM의 가능한 패턴. 우측의 경우 하한은 단조증가하나 가능도는 감소할 수 있다. 이 경우 사후분포의 근사값과 참값간 차이가 작아지는 것이며 이것이 정규화 효과가 될 수 있다.
5개의 가우시안을 혼합한 모델에 대한 풀린 EM 알고리즘의 적용례.

11.5. Model selection for latent variable models

잠재 변수 모델을 쓸 때에는 모델 복잡도를 조절하기 위해 변수의 수를 정해야 한다.

11.5.1. Model selection for probabilistic models

베이지안 최적화는 주변가능도를 최대화시키는 K^{\ast} = \mathrm{argmax}_{k} p(\mathcal{D} | K)을 선택한다. 이는 두 가지 문제점이 있다. 하나는 잠재 변수 모델에 대해서는 주변가능도를 계산하기 어렵다는 점이고 하나는 탐색해야 하는 모델의 범위가 넓다는 것이다. 그 대안으로서 모델 공간에서 추계적 샘플링을 해서 수렴 속도를 크게 가속시킬 수 있다. 이는 각각의 K값에 대해 모델을 피팅하는 것은 느리기 때문이다.

11.5.2. Model selection for non-probabilistic models

확률적 모델이 아닌 경우에는 가능도 대신 재구축 오차 E(\mathcal{D}, K) = \frac{1}{\lvert \mathcal{D} \rvert} \sum_{i \in \mathcal{D}} \lVert \mathbf{x}_{i} - \hat{\mathbf{x}}_{i} \rVert^{2} 을 사용하는 것이다. 이는 K-평균 알고리즘에 유용히 쓸 수 있다.

비지도학습에서는 교차검증을 통한 모델 선택이 불가능하므로 K값에 따른 재구축 오차를 찍어본 뒤 튀어나온 지점을 관찰해 택하는 수밖에 없다. 이는 틈 통계량의 사용으로 인해 자동화될 수 있다.

1D k-평균 알고리즘에서의 모델 선택. 데이터의 분포, 모델에 따른 평균제곱오차와 음의 로그가능도의 변화.

11.6. Fitting models with missing data

결합분포 모델을 최대가능도추정으로 피팅하고 싶은데 데이터에 유실이 존재한다고 하자. \mathbf{X}_{v}를 관측된 데이터, \mathbf{X}_{h}를 유실된 데이터, \mathbf{O}를 데이터의 관측 여부를 나타내는 행렬이라 하면 \hat{\mathbf{\theta}} = \mathrm{argmax}_{\mathbf{\theta}} p(\mathbf{X}_{v} | \mathbf{\theta}, \mathbf{O})을 계산하는 것이 목표가 된다.

로그 가능도는 \log p(\mathbf{X}_{v} | \mathbf{\theta}) = \sum_{i} \log [\sum_{\mathbf{x}_{ih}} p(\mathbf{x}_{iv}, \mathbf{x}_{ih} | \mathbf{\theta})]이 되는데 이는 최적화하기 힘드므로 EM 알고리즘을 이용해 근사할 수 있다.

11.6.1. EM for the MLE of an MVN with missing data

다변수 가우시안의 예를 들어 보자.

11.6.1.1. Getting started

먼저 완전히 관측된 데이터 열을 통해 최대가능도추정을 수행해야 한다. 그런 열이 없다면 휴리스틱으로 데이터 전가를 할 수밖에 없다.

11.6.1.2. E step

기대완전데이터 로그가능도는 다음과 같다.

Q(\mathbf{\theta}, \mathbf{\theta}_{t-1}) = -\frac{N}{2} \log \lvert \mathbf{\Sigma} \rvert - \frac{ND}{2} \log (2 \pi) - \frac{1}{2} \mathrm{tr} (\mathbf{\Sigma}^{-1} \mathbb{E}[(\mathbf{x}_{i} - \mathbf{\mu})(\mathbf{x}_{i} - \mathbf{\mu})^{T}])

이를 계산해보면 p(\mathbf{x}_{ih} | \mathbf{x}_{iv}, \mathbf{\theta}) \sim \mathcal{N}(\mathbf{m}_{i}, \mathbf{V}_{i})에 대해 다음을 얻는다.

\mathbf{m}_{i} = \mathbf{\mu}_{h} + \mathbf{\Sigma}_{hv}\mathbf{\Sigma}_{vv}^{-1}(\mathbf{x}_{iv} - \mathbf{\mu}_{v})

\mathbf{V}_{i} = \mathbf{\Sigma}_{hh} - \mathbf{\Sigma}_{hv}\mathbf{\Sigma}_{vv}^{-1} \mathbf{\Sigma}_{vh}

따라서 \mathbb{E}[\mathbf{x}_{i}] = (\mathbf{m}_{i} ; \mathbf{x}_{iv})이다.

11.6.1.3. M step

\nabla Q(\mathbf{\theta}, \mathbf{\theta}_{t-1}) = \mathbf{0}을 풀면 다음을 얻는다.

\hat{\mathbf{\mu}}_{t} = \frac{1}{N} \sum_{i} \mathbb{E}[\mathbf{x}_{i}]

\hat{\mathbf{\Sigma}}_{t} = \frac{1}{N} \sum_{i} \mathbb{E}[\mathbf{x}_{i} \mathbf{x}_{i}^{T}] - \hat{\mathbf{\mu}}_{t} \hat{\mathbf{\mu}}_{t}^{T}

즉 EM은 최대가능도추정에서 변수들을 단순히 기대값으로 치환하는 것이 아니다. 그렇게 하면 사후분산이 무시되기 때문이다. 충족 통계량의 기대값을 계산하고 최대가능도추정에 대입해야만 맞는 결과가 나온다. 이 알고리즘을 개선해 최대사후확률추정을 수행할 수 있다.

11.6.1.4. Example

데이터 전가 문제는 EM을 통해 풀 수 있다. 당연하게도, 성능은 관측된 데이터가 많을 수록 좋아진다.

데이터 전가 문제.

11.6.1.5. Extension to the GMM case

이는 가우시안의 혼합 모델에 대해 쉽게 확장될 수 있다.

요점 정리

  • 잠재 변수 모델: 은닉 변수가 있는 모델로, 매개변수가 더 적으며 은닉 변수가 모델 정보의 병목 역할을 한다는 이점이 있음.
  • 혼합 모델: 기반 분포를 가중치합해서 얻는 모델. 이산적 잠재 상태를 모델링할 때 쓰임.
  • 혼합 모델의 매개변수 추정을 할 때에는 로그 가능도가 볼록하지 않기 때문에 EM 등으로 국소해를 구해야 함.
  • 기대값 최대화(EM) 알고리즘: 기대완전데이터 로그가능도를 계산하고 (E 단계), 이를 최적화하는 매개변수를 추정해 업데이트하는 (M 단계) 과정을 수렴할 때까지 반복해 로그가능도의 국소해를 찾는 과정.
  • 잠재 변수 모델을 선택할 때에는 주변가능도 또는 재구축 오차를 최적화함.
  • 데이터가 유실된 모델을 피팅할 때 EM 알고리즘을 쓸 수 있음.

답글 남기기

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

WordPress.com 로고

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

Google photo

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

Twitter 사진

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

Facebook 사진

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

%s에 연결하는 중