EM 알고리즘

Expectation-maximization algorithm; EM 알고리즘, 기댓값 최대화 알고리즘

관측되지 않는 잠재 변수가 있는 확률 모형에서, 잠재 변수의 사후 분포를 계산하는 E 단계와 그 기대 완전 데이터 로그 우도를 최대화하는 M 단계를 번갈아 반복해 최대 우도 추정값을 찾는 반복 알고리즘

EM 알고리즘(Expectation-maximization algorithm, 기댓값 최대화 알고리즘)은 데이터 일부가 관측되지 않았거나 모형에 숨은 변수가 있어서 우도를 직접 최대화하기 어려울 때 쓰는 반복 최적화 방법이다. "숨은 변수를 알았다면 쉬웠을" 문제로 바꿔, 현재 매개변수로 숨은 변수를 추측(E 단계)하고 그 추측을 바탕으로 매개변수를 다시 추정(M 단계)하기를 되풀이한다. 매 반복에서 관측 데이터의 우도가 줄지 않는다는 보장이 있다.

비슷한 계산은 이전에도 특수한 문제마다 따로 쓰였는데(1958년 하틀리의 불완전 데이터 최대 우도 추정, 1970년 바움 등의 은닉 마르코프 모형 추정 등), 1977년 아서 뎀프스터(Arthur Dempster), 낸 레어드(Nan Laird), 도널드 루빈(Donald Rubin)이 일반적인 틀로 정리하고 EM이라는 이름을 붙였다. 이 논문은 우도가 단조 증가함을 보이고 결측값, 그룹화·중도 절단 데이터, 유한 혼합 모형, 분산 성분 추정, 요인 분석 등을 예로 들었다. 기계 학습에서는 가우시안 혼합 모델, 은닉 마르코프 모델, 토픽 모델, 기계 번역의 단어 정렬 등 잠재 변수 모형 학습의 기본 도구다.

문제 설정

관측 데이터 $X$, 잠재 변수 $Z$, 매개변수 $\theta$가 있을 때 최대 우도 추정의 목표는 주변 로그 우도

$\ell (\theta )=\log p(X\mid \theta )=\log \sum _{Z}p(X,Z\mid \theta )$

를 최대화하는 것이다($Z$가 연속이면 합 대신 적분). 로그 안에 합이 있어서 미분해 0으로 두어도 닫힌 해가 잘 나오지 않는다. 반대로 $Z$까지 관측된 완전 데이터의 로그 우도 $\log p(X,Z\mid \theta )$는 지수족 모형이라면 로그가 곱을 합으로 풀어 주어 쉽게 최대화된다. 예를 들어 가우시안 혼합 모델에서 각 점이 어느 성분에서 나왔는지 안다면 성분별 표본 평균과 표본 공분산을 구하면 끝이다. EM은 모르는 $Z$를 그 사후 분포에 대한 기댓값으로 대신한다.

유도

증거 하한

$Z$에 대한 임의의 분포 $q(Z)$를 끼워 넣고 로그가 오목 함수라는 젠센 부등식을 쓰면

$\log p(X\mid \theta )=\log \sum _{Z}q(Z)\frac{p(X,Z\mid \theta )}{q(Z)}\ge \sum _{Z}q(Z)\log \frac{p(X,Z\mid \theta )}{q(Z)}=:ℒ(q,\theta )$

이다. $ℒ$을 증거 하한(evidence lower bound, ELBO)이라 한다. 등호가 어디서 성립하는지는 다음 분해로 정확히 드러난다.

$\log p(X\mid \theta )=ℒ(q,\theta )+\mathrm{K}\mathrm{L}(q(Z)\Vert p(Z\mid X,\theta ))$

쿨백-라이블러 발산은 0 이상이고 두 분포가 같을 때만 0이므로, $q$가 잠재 변수의 사후 분포와 같아지면 하한이 로그 우도에 딱 붙는다.

두 단계

EM은 하한 $ℒ(q,\theta )$를 $q$와 $\theta$에 대해 번갈아 최대화하는 좌표 상승법이다.

• E 단계: ${\theta }^{(t)}$를 고정하고 $q$에 대해 최대화한다. 답은 $q^{(t)}(Z)=p(Z\mid X,{\theta }^{(t)})$이고, 이때 $ℒ(q^{(t)},{\theta }^{(t)})=\ell ({\theta }^{(t)})$가 된다. 실제로는 다음 단계에 필요한 충분 통계량의 기댓값만 계산하면 된다.
• M 단계: $q^{(t)}$를 고정하고 $\theta$에 대해 최대화한다. $q$의 엔트로피 항은 $\theta$와 무관하므로 Q 함수

$Q(\theta \mid {\theta }^{(t)})={𝔼}_{Z\sim p(Z\mid X,{\theta }^{(t)})}[\log p(X,Z\mid \theta )]$

를 최대화하는 것과 같다. ${\theta }^{(t+1)}=\arg \underset{\theta }{\max }Q(\theta \mid {\theta }^{(t)})$.

단조 증가

위 두 단계를 이으면

$\ell ({\theta }^{(t+1)})\ge ℒ(q^{(t)},{\theta }^{(t+1)})\ge ℒ(q^{(t)},{\theta }^{(t)})=\ell ({\theta }^{(t)})$

이다. 첫 부등호는 하한이라서, 둘째는 M 단계가 하한을 키워서, 등호는 E 단계에서 하한을 우도에 붙였기 때문이다. 따라서 관측 데이터 로그 우도는 반복마다 줄지 않는다. M 단계에서 최대화까지 하지 않고 $Q$를 늘리기만 해도 이 성질은 유지된다(일반화 EM, GEM).

의사코드

입력: 관측 X, 초깃값 theta, 허용 오차 tol
반복:
    E 단계: q(Z) <- p(Z | X, theta)              # 또는 필요한 기대 충분 통계량
    M 단계: theta <- argmax_theta E_q[log p(X, Z | theta)]
    ll <- log p(X | theta)
    ll 증가량 < tol 이면 멈춤

예시: 두 동전 문제

앞면이 나올 확률이 ${\theta }_A$, ${\theta }_B$인 동전 두 개 중 하나를 골라 10번 던지는 실험을 5번 했는데, 매번 어느 동전을 골랐는지는 기록하지 않았다고 하자. 앞면 수는 5, 9, 8, 4, 7이다. 동전 정체가 잠재 변수다. 도(Chuong B. Do)와 바초글루(Serafim Batzoglou)가 EM 소개 글에서 쓴 예다.

• E 단계: 묶음 $i$의 앞면 수가 $h_i$일 때 동전 A였을 사후 확률 $w_i=\frac{{\theta }_A^{h_i}(1-{\theta }_A{)}^{10-h_i}}{{\theta }_A^{h_i}(1-{\theta }_A{)}^{10-h_i}+{\theta }_B^{h_i}(1-{\theta }_B{)}^{10-h_i}}$(두 동전을 고를 확률은 같다고 둔다).
• M 단계: 가중 앞면 비율 ${\theta }_A=\frac{\sum _i{w}_i{h}_i}{\sum _{i}10w_i}$, ${\theta }_B$도 $1-w_i$로 같은 방식.

import numpy as np
from scipy.stats import binom

heads = np.array([5, 9, 8, 4, 7])   # 묶음마다 10번 던져 나온 앞면 수
n = 10
theta_A, theta_B = 0.6, 0.5          # 초깃값

for it in range(1, 11):
    # E 단계: 각 묶음이 동전 A에서 나왔을 사후 확률
    la = binom.pmf(heads, n, theta_A)
    lb = binom.pmf(heads, n, theta_B)
    w_A = la / (la + lb)
    w_B = 1 - w_A
    # M 단계: 사후 확률로 가중한 앞면 비율
    theta_A = (w_A * heads).sum() / (w_A * n).sum()
    theta_B = (w_B * heads).sum() / (w_B * n).sum()
    # 관측 데이터 로그 우도(두 동전을 반반 확률로 고른다고 둘 때)
    ll = np.log(0.5 * binom.pmf(heads, n, theta_A) + 0.5 * binom.pmf(heads, n, theta_B)).sum()
    if it in (1, 2, 10):
        print(f"{it:2d}회: theta_A={theta_A:.2f}, theta_B={theta_B:.2f}, 로그 우도={ll:.4f}")

 1회: theta_A=0.71, theta_B=0.58, 로그 우도=-10.0860
 2회: theta_A=0.75, theta_B=0.57, 로그 우도=-9.9498
10회: theta_A=0.80, theta_B=0.52, 로그 우도=-9.7969

첫 반복 값(0.71, 0.58)과 10회 뒤 값(0.80, 0.52)은 원 글의 그림과 같고, 로그 우도는 반복마다 커진다(NumPy 2.1, SciPy 1.15에서 실행). 동전 정체를 알았다면 동전별 앞면 비율을 바로 구했을 것이고, 0 또는 1로 딱 잘라 배정하는 대신 사후 확률로 나눠 세는 것이 EM의 핵심이다.

예시: 가우시안 혼합 모델

$K$개 정규 분포의 혼합 $p(x)=\sum _k{\pi }_k𝒩(x\mid {\mu }_k,{\Sigma }_k)$에서 잠재 변수는 각 점의 성분 번호다.

• E 단계: 책임값(responsibility) ${\gamma }_{ik}=\frac{{\pi }_k𝒩(x_i\mid {\mu }_k,{\Sigma }_k)}{\sum _j{\pi }_j𝒩(x_i\mid {\mu }_j,{\Sigma }_j)}$
• M 단계: $N_k=\sum _i{\gamma }_{ik}$, ${\pi }_k=N_k/N$, ${\mu }_k=\frac{1}{N_k}\sum _i{\gamma }_{ik}{x}_i$, ${\Sigma }_k=\frac{1}{N_k}\sum _i{\gamma }_{ik}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^{\mathrm{\top }}$

책임값을 0 또는 1로 굳히고 공분산을 같은 크기의 구형으로 고정하면 k-평균의 두 단계(가장 가까운 중심에 배정, 중심 다시 계산)가 된다. 구현과 실무 내용은 가우시안 혼합 모델 문서 참고.

수렴 성질

• 정류점 수렴: 우도가 위로 유계면 우도 값의 열은 수렴한다. 다만 매개변수 열이 어디로 가는지는 따로 따져야 한다. 우(C. F. Jeff Wu)는 1983년 완전 데이터가 (곡선) 지수족이고 매개변수 공간이 콤팩트하면 EM 열의 모든 극한점이 우도의 정류점이며, 우도가 봉우리 하나이고 미분 가능성 조건을 만족하면 유일한 최대 우도 추정값으로 수렴함을 보였다. 정류점에는 국소 최댓값뿐 아니라 안장점도 포함된다.
• 국소 최적해: 혼합 모형처럼 우도에 봉우리가 여럿이면 초깃값에 따라 다른 답에 도달한다. 여러 초깃값으로 돌려 우도가 가장 큰 결과를 고르거나, k-평균 결과로 초기화한다.
• 속도: 수렴은 보통 선형(1차)이고, 그 속도는 전체 정보 가운데 결측 정보가 차지하는 비율로 정해진다. 숨은 부분이 클수록(군집이 많이 겹칠수록) 느려진다. 최적해 근처에서 느려지는 문제는 뉴턴 방법이나 준뉴턴 방법으로 마무리해 보완한다.
• 퇴화: 가우시안 혼합에서 한 성분이 데이터 점 하나에 붙어 분산이 0으로 가면 우도가 무한대가 된다. 공분산 대각에 작은 값을 더하거나 사전 분포를 둔다.
• 멈춤 기준: 로그 우도 증가량이나 매개변수 변화량이 허용 오차보다 작아지면 멈춘다. 로그 우도가 줄었다면 구현 오류다. 단조성은 좋은 디버깅 도구다.

변형

변형	바뀌는 부분	쓰는 경우
일반화 EM(GEM)	M 단계에서 $Q$를 늘리기만 한다(경사 한 걸음 등)	M 단계의 닫힌 해가 없을 때
ECM	M 단계를 매개변수 묶음별 조건부 최대화 여러 번으로 나눈다	전체 최대화는 어렵지만 일부를 고정하면 쉬울 때
몬테카를로 EM	E 단계의 기댓값을 표본 평균으로 근사	사후 분포의 기댓값을 계산할 수 없을 때(MCMC와 결합)
온라인(확률적) EM	미니배치마다 충분 통계량을 이동 평균으로 갱신	데이터가 커서 전체를 한 번에 못 볼 때
점진적·희소 EM	데이터 일부의 $q$만 갱신해도 하한이 늘어난다는 관점	대규모 데이터, 책임값 대부분이 0일 때
하드 EM(비터비 학습)	사후 분포 대신 가장 그럴듯한 $Z$ 하나로 M 단계	빠르지만 우도 단조성 보장이 약하다. k-평균이 대표 예
변분 EM	정확한 사후 분포 대신 다루기 쉬운 분포족 안에서 $q$를 최적화	사후 분포가 다루기 어려울 때. 잠재 디리클레 할당의 원 논문 방식
MAP-EM	M 단계 목적에 로그 사전 분포를 더한다	퇴화 방지, 데이터가 적을 때

응용

모형·과제	잠재 변수	비고
가우시안 혼합 모델	각 점의 성분 번호	연성 군집, 밀도 추정
은닉 마르코프 모델	숨은 상태 열	바움-웰치 알고리즘. E 단계를 전방-후방 알고리즘으로 계산한다
결측값이 있는 다변량 정규 모형	빠진 값	결측값을 기대값으로 채우면서 평균·공분산 추정(결측치 처리)
통계적 기계 번역의 단어 정렬	원문 단어와 번역 단어의 대응	IBM 모델 1~5의 학습
준지도 텍스트 분류	레이블 없는 문서의 클래스	적은 레이블과 많은 무레이블 문서로 나이브 베이즈 학습
전문가 혼합	입력마다 담당 전문가	계층적 전문가 혼합을 EM으로 학습 오늘날 전문가 혼합 신경망은 경사 하강으로 학습한다
요인 분석, 확률적 PCA	저차원 잠재 요인	주성분 분석의 확률 모형판
베이지안 네트워크	관측 안 된 노드, 결측값	조건부 확률표를 기대 빈도로 추정

다른 방법과의 관계

• 경사 상승법: 로그 우도를 직접 미분해 올라가도 된다. 실제로 ${\mathrm{\nabla }}_{\theta }\ell (\theta )={\mathrm{\nabla }}_{\theta }Q(\theta \mid {\theta }^{(t)}){|}_{\theta ={\theta }^{(t)}}$이므로 두 방법은 같은 정류점을 찾는다. EM은 학습률이 없고 제약(확률의 합이 1, 공분산이 양의 정부호)을 M 단계 닫힌 해가 자동으로 지킨다는 장점이 있다. 신경망처럼 M 단계 닫힌 해가 없는 모형은 경사법이 자연스럽다.
• 변분 추론과 VAE: 증거 하한 $ℒ(q,\theta )$는 변분 추론과 변분 오토인코더가 최대화하는 바로 그 ELBO다. EM은 $q$를 정확한 사후 분포로 맞출 수 있는 특수한 경우이고, VAE는 데이터마다 $q$를 최적화하는 대신 인코더 신경망이 $q(z\mid x)$를 한 번에 내놓게 하고(상각 추론) $q$와 $\theta$를 경사법으로 함께 갱신한다.
• MCMC: 점추정 대신 매개변수의 사후 분포 전체가 필요하면 깁스 샘플링 같은 마르코프 연쇄 몬테카를로를 쓴다. 깁스 샘플링은 "잠재 변수를 기댓값으로 채우는" 대신 "표본으로 채우는" EM의 확률적 형제로 볼 수 있다.

한계

• 국소 최적해와 초깃값 의존성이 크다.
• 결측 정보가 많으면 수렴이 매우 느리다.
• E 단계에 필요한 사후 분포가 다루기 어려운 모형(복잡한 그래프 모형, 신경망 생성 모형)에는 그대로 쓸 수 없고 근사가 필요하다.
• 점추정만 주므로 매개변수 불확실성은 따로 구해야 한다(관측 정보 행렬, 부트스트랩 등).

같이 보기

• 최대 우도 추정
• 가우시안 혼합 모델
• 은닉 마르코프 모델
• 잠재 디리클레 할당
• 변분 추론
• 변분 오토인코더
• 쿨백-라이블러 발산
• 마르코프 연쇄 몬테카를로
• K-평균 군집화
• 비지도 학습
• 우도

참고 문헌

• Christopher M. Bishop, 《Pattern Recognition and Machine Learning》, Springer, 2006, 9장. 마이크로소프트 리서치 공개본
• Geoffrey J. McLachlan, Thriyambakam Krishnan, 《The EM Algorithm and Extensions》 2판, Wiley, 2008.