요약

가우스 혼합 모델의 매개변수를 자료에서 찾는 방법이다. “각 점이 어느 무리에서 왔는지”를 알면 무리마다 평균과 퍼짐을 바로 계산할 수 있고, 반대로 평균과 퍼짐을 알면 각 점이 어느 무리에서 왔을 확률을 계산할 수 있다. 둘 다 모르니 아무렇게나 시작해 두 계산을 번갈아 되풀이한다. 되풀이할 때마다 자료가 나올 가능성은 줄지 않지만, 가장 좋은 답이 아닌 곳에 멈출 수 있어 처음 값이 중요하다.

예시로 보기

GMM의 책임도를 계산하려면 무리별 평균, 분산, 비중이 필요한데, 그것을 구하려면 어느 점이 어느 무리인지 알아야 한다. 닭이 먼저냐 달걀이 먼저냐다.

1차원 점 0, 1, 2, 6, 7, 8을 두 무리로 나눈다. 처음 평균 1과 2, 분산 1, 비중 \(\frac12\)로 시작한다1.

E 단계 (책임도). 무리 1의 책임도를 점마다 계산한다.

점 0 1 2 6 7 8
무리 1 책임도 0.818 0.622 0.378 0.011 0.004 0.002

점 0은 평균 1에 더 가까워 무리 1의 책임이 크고, 6, 7, 8은 평균 2 쪽(무리 2)이 훨씬 그럴듯하다.

M 단계 (다시 맞추기). 책임도를 무게로 삼아 무리마다 가중 평균, 가중 분산, 비중을 다시 계산한다.

  평균 분산 비중
무리 1 0.809 0.885 0.306
무리 2 5.405 7.076 0.694

무리 2의 평균이 2에서 5.4로 크게 움직였다. 이를 되풀이하면 평균 1과 7, 비중 \(\frac12\)씩으로 다가간다.

그림

같은 자료로 E 단계와 M 단계를 15번 되풀이했다. 무리 2의 평균은 2 → 5.4 → … → 7로 올라가고, 무리 1은 1 근처에 자리 잡는다. 오른쪽의 로그가능도는 한 번도 줄지 않고, 평균이 자리를 잡으면 함께 멈춘다2.

검증: 위 두 표, 수렴값, 자료 2,000개에서 매개변수 되찾기, 무작위 자료 100개에서 로그가능도가 한 번도 줄지 않음 — 30_em-gmm_impl.py

정의

GMM의 매개변수를 EM으로 추정하는 순서는 다음과 같다. ① 매개변수를 정한다 ② 군집 배정(책임도)을 계산한다 ③ 새 배정으로 매개변수를 고친다 ④ 수렴할 때까지 되풀이한다.

  • E 단계(기댓값): 지금 매개변수로 각 점의 군집 소속 확률을 계산한다. 어느 가우스 분포가 각 점을 만들었는지를 사후 확률로 나눠 준다.
  • M 단계(최대화): 그 확률을 놓고, 자료를 만들 가능성(가능도)이 가장 커지도록 매개변수를 고친다.
flowchart TD
    S["평균·공분산·비중의 처음 값을 정한다"] --> E["E 단계: 지금 매개변수로 점마다 책임도를 계산한다"]
    E --> M["M 단계: 책임도를 무게로 매개변수를 다시 계산한다"]
    M --> Q{"수렴했는가?"}
    Q -->|"아니다"| E
    Q -->|"그렇다"| F["매개변수와 책임도를 돌려준다"]

E 단계는 매개변수를 고정하고 소속 확률을, M 단계는 소속 확률을 고정하고 매개변수를 고친다. k-평균의 배정·갱신 고리와 모양이 같다.

의사코드

GMM의 EM을 의사코드로 쓰면 이렇다.

μ_α, Σ_α는 무작위로, π_α는 1/k로 정한다
repeat
    E 단계: 모든 i, α에 대해
        [z_i]_α = π_α N(x_i | μ_α, Σ_α) / Σ_l π_l N(x_i | μ_l, Σ_l)
    M 단계: α = 1..k에 대해
        μ_α = Σ_i [z_i]_α x_i / Σ_i [z_i]_α                         # 새 중심
        Σ_α = Σ_i [z_i]_α (x_i − μ_α)(x_i − μ_α)ᵀ / Σ_i [z_i]_α        # 새 모양
        π_α = (1/n) Σ_i [z_i]_α                                      # 새 비중
until 수렴

M 단계는 책임도를 무게로 한 평균과 공분산이다. 책임도가 큰 점일수록 그 군집의 계산에 많이 들어간다.

처음 비중은 \(1/k\)로 둔다. 군집 \(k\)개의 비중을 모두 더하면 1이 되어야 하기 때문이다. M 단계 식의 \(\frac1n\)은 책임도 합을 점 수 \(n\)으로 나누는 것이라 처음 값과 다르다.

정확성: 가능도는 줄지 않는다

주장. 매 반복에서 로그가능도 \(\ell = \sum_i \log\sum_\alpha \pi_\alpha\mathcal N(\mathbf x_i \mid \boldsymbol\mu_\alpha, \Sigma_\alpha)\)는 줄지 않는다1.

증명 스케치
  1. 아래 경계. 아무 확률 \(q_{i\alpha}\)(\(\sum_\alpha q_{i\alpha} = 1\))에 대해 \(\log\sum_\alpha \pi_\alpha\mathcal N_{i\alpha} = \log\sum_\alpha q_{i\alpha}\frac{\pi_\alpha\mathcal N_{i\alpha}}{q_{i\alpha}} \ge \sum_\alpha q_{i\alpha}\log\frac{\pi_\alpha\mathcal N_{i\alpha}}{q_{i\alpha}}\). — 로그는 오목해서 평균의 로그 ≥ 로그의 평균(옌센 부등식)
  2. E 단계는 경계를 딱 붙인다. \(q_{i\alpha}\)를 책임도 \([\mathbf z_i]_\alpha\)로 두면 \(\frac{\pi_\alpha\mathcal N_{i\alpha}}{q_{i\alpha}}\)가 \(\alpha\)에 상관없이 \(\sum_l\pi_l\mathcal N_{il}\)로 같아져, 1단계의 부등호가 등호가 된다. — 대입
  3. M 단계는 경계를 올린다. \(q\)를 고정하고 오른쪽 식을 매개변수에 대해 최대화한다. 그 답이 위 의사코드의 가중 평균·가중 공분산·비중이다. — 미분과 최소화
  4. 합치기. 옛 매개변수에서 \(\ell_{\text{옛}}\) = 경계(2단계) ≤ 새 매개변수의 경계(3단계) ≤ \(\ell_{\text{새}}\)(1단계). ∎

[증명 스케치] 3단계의 최댓값 계산은 생략했다. Bishop, Pattern Recognition and Machine Learning, 9.2~9.4절을 따른다.

스스로 설명해 보기

1. 1단계에 옌센 부등식을 쓰는 이유

로그 안의 합 \(\log\sum\)은 다루기 어렵다. 부등식으로 합을 로그 밖으로 꺼내면 \(\sum\log\)가 되어, 각 군집의 식이 따로 떨어져 M 단계에서 쉽게 최대화된다.

2. 2단계에서 책임도를 고르면 등호가 되는 이유

옌센 부등식은 평균 내는 값들이 모두 같을 때 등호다. 책임도로 나누면 \(\frac{\pi_\alpha\mathcal N_{i\alpha}}{[\mathbf z_i]_\alpha} = \sum_l\pi_l\mathcal N_{il}\)로 모든 \(\alpha\)에서 같은 값이 된다.

3. “줄지 않는다”가 “가장 좋은 답에 간다”는 뜻이 아닌 이유

매 단계가 지금 위치 근처의 경계만 올린다. 가능도에 봉우리가 여럿이면 처음 값이 가까운 봉우리에서 멈춘다(국소 최적). 그래서 처음 값을 바꿔 여러 번 돌린다.

이 증명의 핵심 아이디어는?

어려운 목적 함수 아래에 다루기 쉬운 경계를 깔고(E), 그 경계를 올린다(M). 경계가 현재 위치에서 목적 함수에 붙어 있으므로 경계를 올리면 목적 함수도 올라간다.

이 방법을 쓸 수 있는 다른 상황은?

빠진 값이 있는 자료의 최대가능도 추정, 은닉 마르코프 모델의 학습(바움-웰치), 그리고 하드 배정 버전인 k-평균(배정 = E, 평균 = M).

복잡도

한 번 반복에 점 \(n\)개와 군집 \(k\)개의 책임도를 모두 계산한다. \(d\)차원 공분산을 쓰면 가우스 밀도 계산에 역행렬과 행렬식이 들어 \(O(nkd^2 + kd^3)\)이다1. 차원이 크면 공분산 추정이 불안정해진다.

활용

  • 구현: 30_em-gmm_impl.py(1차원)
  • 처음 값으로 k-평균의 결과를 쓰는 경우가 많다(사이킷런 GaussianMixture의 기본 설정)1.

연결

자주 하는 오해

“EM은 GMM의 참 매개변수를 반드시 찾는다”

아니다. EM은 가능도를 줄이지 않을 뿐, 처음 값에 따라 다른 봉우리에 멈춘다. 또 한 성분이 점 하나에만 붙으면 분산이 0으로 줄며 가능도가 무한대로 커지는 퇴화도 생긴다(구현에서 분산 아래한계를 두는 이유). 그럴듯해 보이는 이유는 예시처럼 잘 떨어진 자료에서는 거의 늘 참값에 가기 때문이다. 확인하는 법: 처음 값을 바꿔 여러 번 돌려 로그가능도를 비교한다.

확인 문제

C1 EM의 E 단계와 M 단계가 각각 계산하는 것을 식과 함께 쓰라.

답: E: 책임도 \([\mathbf z_i]_\alpha = \frac{\pi_\alpha\mathcal N(\mathbf x_i \mid \boldsymbol\mu_\alpha, \Sigma_\alpha)}{\sum_l\pi_l\mathcal N(\mathbf x_i \mid \boldsymbol\mu_l, \Sigma_l)}\). M: \(\boldsymbol\mu_\alpha\)는 책임도 가중 평균, \(\Sigma_\alpha\)는 책임도 가중 공분산, \(\pi_\alpha = \frac1n\sum_i[\mathbf z_i]_\alpha\).

C2 점 0, 1, 2, 6, 7, 8과 처음 값(평균 1, 2, 분산 1, 비중 반반)에서 첫 E 단계 뒤 무리 1의 책임도 합과, M 단계의 무리 1 비중을 구하라.

답: 책임도 0.818 + 0.622 + 0.378 + 0.011 + 0.004 + 0.002 ≈ 1.835. 비중 \(\pi_1 = \frac{1.835}{6} \approx 0.306\).

C3 의사코드의 μ_α = Σ_i [z_i]_α x_i / Σ_i [z_i]_α 줄이 하는 일을 한 문장으로 쓰라.

답: 각 점이 군집 \(\alpha\)에 속할 확률을 무게로 삼아 점들의 가중 평균을 내서, 군집 \(\alpha\)의 새 중심을 정한다.

C4 증명 4단계의 부등식 사슬 \(\ell_{\text{옛}} \le \ell_{\text{새}}\)에서 첫 등호와 둘째·셋째 부등호의 근거를 각각 대라.

답: 첫 등호(\(\ell_{\text{옛}}\) = 경계): E 단계에서 책임도를 고르면 옌센 부등식이 등호가 된다. 둘째(경계 ≤ 새 경계): M 단계가 그 경계를 최대화한다. 셋째(새 경계 ≤ \(\ell_{\text{새}}\)): 옌센 부등식은 아무 \(q\)에서나 맞다.

C5 \(\pi_\alpha\)를 \(1/k\)가 아니라 \(1/n\)(\(n\)은 점 수)으로 시작하면 무엇이 깨지는가? 점 6개, 군집 2개로 보여라.

답: 비중이 \(\frac16, \frac16\)이라 합이 \(\frac13\)이다. \(p(\mathbf x) = \sum\pi_\alpha\mathcal N\)이 확률밀도가 아니게 되어(전체 적분이 1/3) 첫 로그가능도가 틀린다. 책임도는 분자·분모의 공통 배수가 약분되어 우연히 같지만, 비중의 뜻은 깨진다. \(1/k = \frac12\)로 시작해야 한다.

  1. 구현 코드로 확인했다(Dempster, Laird, Rubin, 1977). ↩ ↩2 ↩3 ↩4

  2. 30_em-algorithm_plot.py로 그렸고, 위 두 표의 값(책임도 0.818 … 0.002, 평균 0.809·5.405, 분산 0.885·7.076, 비중 0.306·0.694), 로그가능도가 줄지 않음, 평균 1과 7·비중 \(\frac12\)로 수렴함을 같은 코드로 확인했다. ↩