4 계산 — 정리를 알고리즘으로

노트북 4절의 MVN 클래스는 아래 세 알고리즘을 그대로 옮긴 것이다. 원칙은 하나다. 역행렬을 만들지 않는다. $\Sigma^{-1}$ 을 직접 계산하면 비용이 들고, $\Sigma$ 가 특이에 가까울 때 수치 오차가 커진다. Cholesky 인수 $L$ 로 삼각 연립방정식을 푸는 편이 빠르고 안정적이다.

알고리즘 4.1 (표본 추출) 입력 $\bmu$ , $\Sigma\succ0$ , 표본 수 $N$ . 출력 $\x^{(1)},\allowbreak \dots,\allowbreak \x^{(N)}\iid\N_n(\bmu,\allowbreak \Sigma)$ .

1.
Cholesky 분해 $\Sigma=LL\T$ (한 번, 비용 $O(n^3)$ ).
2.
각 $k$ 마다 $\z^{(k)}\sim\N_n(\mathbf0,I)$ 를 뽑고 $\x^{(k)}=\bmu+L\z^{(k)}$ (각 $O(n^2)$ ).

정당성: 정리 3.5에서 $A=L$ .

알고리즘 4.2 (로그밀도) 입력 $\x$ , $\bmu$ , $L$ . 출력 $\log f(\x)$ .

1.
전진대입으로 $L\mathbf u=\x-\bmu$ 를 푼다( $O(n^2)$ ).
2.
$Q=\mathbf u\T\mathbf u$ . (이유: $\Sigma^{-1}=L^{-\prime}L^{-1}$ 이므로 $Q=\|L^{-1}(\x-\bmu)\|^2$ .)
3.
$\log|\Sigma|=2\sum_i\log L_{ii}$ . (이유: $|\Sigma|=|L|^2=\prod L_{ii}^2$ .)
4.
$\log f(\x)=-\tfrac12\bigl(n\log2\pi+\log|\Sigma|+Q\bigr)$ .

알고리즘 4.3 (조건부분포) 입력 블록 $\Sigma_{11},\allowbreak \Sigma_{12},\allowbreak \Sigma_{22}$ , $\bmu$ , 관측 $\x_2$ . 출력 $\bmu_{1\cdot2}$ , $\Sigma_{11\cdot2}$ .

1.
연립방정식 $\Sigma_{22}\,G=\Sigma_{21}$ 을 풀어 $B=G\T=\Sigma_{12}\Sigma_{22}^{-1}$ 을 얻는다.
2.
$\bmu_{1\cdot2}=\bmu_1+B(\x_2-\bmu_2)$ , $\Sigma_{11\cdot2}=\Sigma_{11}-B\Sigma_{21}$ .
3.
반올림으로 생긴 비대칭을 없애려고 $\Sigma_{11\cdot2}\leftarrow\tfrac12(\Sigma_{11\cdot2}+\Sigma_{11\cdot2}\T)$ .

반정치인 경우 Cholesky는 $\Sigma\succ0$ 에서만 작동한다. $\Sigma$ 가 특이하거나 수치적으로 거의 특이하면(6.1절의 GP 공분산이 그렇다) 고유분해 $\Sigma=V\Lambda V\T$ 로 $A=V\Lambda^{1/2}$ 를 만들어 알고리즘 4.1을 그대로 쓴다. 이때 밀도는 없으므로 알고리즘 4.2는 쓸 수 없다(명제 5.6).

연습문제

연습 4.1 (예제 3.16을 코드로 ★) 노트북의 MVN.conditional로 예제 3.16의 $\bmu_{1\cdot2}$ 와 $\Sigma_{11\cdot2}$ 를 $(x_2,x_3)=(1,0)$ 에서 계산하고, $10^6$ 개 표본에서 $|X_2-1|<0.05$ , $|X_3|<0.05$ 인 점들의 $X_1$ 평균·분산과 비교하라.

풀이

이론값은 평균 $5/3=1.667$ , 분산 $4/3=1.333$ 이다. $10^6$ 개 가운데 상자에 들어오는 점은 600여 개라서 Monte Carlo 평균의 표준오차가 $\sqrt{1.33/600}\approx0.05$ 정도다. 노트북의 한 실행에서는 $1.688\pm0.089$ (2 표준오차)였다.

연습 4.2 (KL 발산 구현 ★★) MVN에 kl(self, other)를 추가하라. 식은 명제 6.2이다. 역행렬 없이 $\tr(\Sigma_1^{-1}\Sigma_0)$ 와 이차형식을 계산하는 방법을 설명하고, Monte Carlo 추정 $\tfrac1N\sum_k[\log f_0(\x^{(k)})-\log f_1(\x^{(k)})]$ 과 대조하라.

풀이

$\Sigma_1=L_1L_1\T$ , $\Sigma_0=L_0L_0\T$ 로 쓰면

\begin{gather*} \tr(\Sigma_1^{-1}\Sigma_0)=\bigl\|L_1^{-1}L_0\bigr\|_F^2,\qquad (\bmu_1-\bmu_0)\T\Sigma_1^{-1}(\bmu_1-\bmu_0)=\bigl\|L_1^{-1}(\bmu_1-\bmu_0)\bigr\|^2,\\ \log|\Sigma_k|=2\textstyle\sum_i\log(L_k)_{ii}. \end{gather*}

$L_1^{-1}L_0$ 은 삼각방정식 $L_1M=L_0$ 을 풀어 얻는다. 노트북 4.6절의 두 분포에서 정답은 $0.8621$ 이고, Monte Carlo 추정( $2\times10^5$ 개)은 $0.860$ 이었다.

연습 4.3 (아들과 딸 ★★) GaltonFamilies에서 아들과 딸을 나눠 midparentHeight에 대한 회귀를 각각 구하라. 기울기와 잔차 표준편차를 비교해, Galton의 ×1.08 손질이 평균뿐 아니라 회귀 구조까지 맞추는지 논하라.

풀이

아들: 기울기 $0.713$ , 잔차 표준편차 $2.30$ . 딸: 기울기 $0.661$ , 잔차 표준편차 $2.02$ . 딸에 1.08을 곱하면 기울기 $0.661\times1.08=0.714$ , 잔차 표준편차 $2.18$ 로 아들과 거의 같아진다. 1.08은 평균만이 아니라 회귀 구조까지 꽤 잘 맞춘다. 잔차 분산은 아직 조금 다르다.

이 절을 PDF로 보기 (35쪽)노트북 4절

이 회차에서 읽은 절 0개 · 회차 차례

4 계산 — 정리를 알고리즘으로 ​

4 계산 — 정리를 알고리즘으로