18. State space models

18.1. Introduction

상태 공간 모델(SMM)은 은닉 상태가 연속적인 은닉 마르코프 모델이다. 해당 모델은 다음의 일반적 형태로 나타낼 수 있다:

\mathbf{z}_{t} = g(\mathbf{u}_{t}, \mathbf{z}_{t-1}, \mathbf{\epsilon}_{t})

\mathbf{y}_{t} = h(\mathbf{z}_{t}, \mathbf{u}_{t}, \mathbf{\delta}_{t})

여기서 각각의 변수는 은닉 상태 \mathbf{z}_{t}, 제어 신호 \mathbf{u}_{t}, 관측 \mathbf{y}_{t}, 전이 모델 g, 관측 모델 h, \mathbf{\epsilon}_{t}는 시스템 노이즈, \mathbf{\delta}_{t}는 관측 노이즈를 나타낸다. 모델 매개변수 \mathbf{\theta}는 알려져 있다고 가정하며, 모르는 매개변수는 은닉 상태로 포함시킬 수 있다.

상태 공간 모델의 주 목적은 믿음 상태 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t}, \mathbf{\theta})를 재귀적으로 추정하는 것이다. 또한 사후예측분포 p(\mathbf{y}_{t + 1} | \mathbf{y}_{1 : t})를 통해 믿음 상태를 바탕으로 미래 관측값을 예측할 수도 있다.

상태 공간 모델의 중요한 특수 케이스는 모든 조건부확률분포가 선형 가우시안인 경우이다. 다시 말해,

\mathbf{z}_{t} = \mathbf{A}_{t} \mathbf{z}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t} + \mathbf{\epsilon}_{t}

\mathbf{y}_{t} = \mathbf{C}_{t} \mathbf{z}_{t} + \mathbf{D}_{t} \mathbf{u}_{t} + \mathbf{\delta}_{t}

\mathbf{\epsilon}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_{t})

\mathbf{\delta}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_{t})

인 경우를 말한다. 이를 선형-가우시안 상태 공간 모델 (LG-SSM) 또는 선형 동역학계 (LDS) 라고도 한다. 만약 모델 매개변수 \mathbf{\theta}_{t} = (\mathbf{A}_{t}, \mathbf{B}_{t}, \mathbf{C}_{t}, \mathbf{D}_{t}, \mathbf{Q}_{t}, \mathbf{R}_{t})이 시간에 독립적이라면, 이 모델은 정적이라 한다. 이 케이스가 중요한 이유는 정확한 추론이 가능하기 때문인데, 최초의 믿음 상태가 가우시안 p(\mathbf{z}_{1}) = \mathcal{N}(\mathbf{\mu}_{1 | 0}, \mathbf{\Sigma}_{1 |0})이라면 모든 믿음 상태 또한 가우시안 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}) = \mathcal{N}(\mathbf{\mu}_{t | t}, \mathbf{\Sigma}_{t | t})이 된다. 이는 칼만 필터로 계산할 수 있다.

18.2. Applications of SSMs

상태 공간 모델은 여러 용례가 있다. 여기서는 선형 가우시안 상태 공간 모델을 가정한다.

18.2.1. SSMs for object tracking

칼만 거르기의 최초 응용 중 하나는 물체를 추적하는 것이다. 2D 평면 위에서 움직이는 물체의 위치를 z_{1t}, z_{2t}, 속도를 \dot{z}_{1t}, \dot{z}_{2t}라 하자. 물체는 등속으로 운동하되 무작위 가우시안 노이즈 \mathbf{\epsilon}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})이 있다고 하자. 그러면 해당 계의 동역학을 \mathbf{z}_{t} = \mathbf{A}_{t} \mathbf{z}_{t-1} + \mathbf{\epsilon}_{t}로 모델링할 수 있다. (\mathbf{A}_{t}는 우측 상단의 부분행렬이 표본 주기 \Delta이고 나머지는 단위행렬인 행렬). 이를 무작위 가속 모델이라 한다.

이제 물체의 속도는 측정할 수 없고 위치만 측정할 수 있다고 하자. \mathbf{y}^{t}를 해당 물체의 위치라 하면 \mathbf{y}_{t} = \mathbf{C}_{t} \mathbf{z}_{t} + \mathbf{\delta}_{t}로 놓을 수 있다. (\mathbf{C}_{t}는 왼쪽 부분이 단위행렬인 2 x 4 행렬, \mathbf{\delta}_{t} \sim \mathcal{N}(\mathbf{0}, \mathbf{R})은 측정 노이즈)

이제 물체의 초기 상태 p(\mathbf{z}_{1}) = \mathcal{N}(\mathbf{z}_{1} | \mathbf{\mu}_{1 | 0}, \mathbf{\Sigma}_{1 | 0})을 가우시안으로 정의하면 p(\mathbf{z}_{t} | \mathbf{y}_{1:t})을 칼만 거르기로 계산해 연속적 베이지안 업데이트를 수행할 수 있다.

칼만 필터링과 스무딩의 예

18.2.2. Robotic SLAM

2D 평면을 움직이는 로봇에 대해 고려해 보자. 이 때 로봇은 지도에 대해 학습하고 지도 상의 자신의 위치를 계속 추적해야 하는데 이를 동시 위치 추적 및 매핑(SLAM) 이라 한다. 지도를 고정된 K개의 랜드마크 L_{1}, \cdots, L_{K}라고 할 때, \mathbf{x}_{t}를 로봇의 위치라고 하면 상태 공간을 \mathbf{z}_{t} = (\mathbf{x}_{t}, \mathbf{L}^{1 : K})로 둘 수 있다. \mathbf{y}_{t}를 가장 가까운 랜드마크와의 거리라 할 때 로봇은 랜드마크 좌표의 추정값을 로봇의 시야를 통해 업데이트할 수 있다. 관측 모델 p(\mathbf{y}_{t} | \mathbf{z}_{t}, \mathbf{L})을 선형-가우시안으로 가정하고, p(\mathbf{x}_{t} | \mathbf{x}_{t-1}, \mathbf{u}_{t})을 가우시안 움직임 모델로 가정할 때, 칼만 거르기를 통해 로봇의 위치와 랜드마크의 위치를 측정할 수 있다.

시간에 따라 로봇의 위치에 대한 불확실성은 증가하는데, 로봇이 익숙한 위치로 돌아가게 된다면 이 불확실성은 다시 감소한다. 이를 루프 닫힘이라 한다. 믿음 상태가 가우시안이므로, 사후공분산행렬 \mathbf{\Sigma}_{t}를 시각화할 수 있다. 보통은 그 역행렬 \mathbf{\Lambda}_{t} = \mathbf{\Sigma}_{t}^{-1}을 시각화하고는 하는데, 이 역행렬은 희소한 특성을 가지기 때문이다. 사후정밀도행렬이 희소한 이유는 이 행렬의 성분이 각각 해당되는 비방향 가우시안 그래프 모형의 간선 여부에 대응하기 때문이다. 최초에는 모든 랜드마크가 서로 연결되어 있지 않으므로 \mathbf{\Lambda}_{t}는 대각행렬인데, 로봇이 움직임에 따라 근처 랜드마크간 상호 연관이 생기고 사후정밀도행렬이 조밀해져 예측에 O(K^{3}) 시간이 걸리게 된다. 이를 꼬임 문제라 하며 이 때문에 큰 지도에 알고리즘을 적용하는 데 있어 문제가 생긴다.

이에 대한 해결책은 크게 두 가지가 있는데 첫째는 상호 연관 패턴이 로봇의 위치를 따라 이동한다는 것에 착안해 시간에 따라 약해지는 연관성들을 가지치기하는 방법이고, 둘째는 로봇의 경로 \mathbf{x}_{1:t}에 조건을 걸고 랜드마크의 위치에는 조건을 걸지 않음으로써 p(\mathbf{L} | \mathbf{x}_{1:t}, \mathbf{y}_{1:t}) = \prod_{k=1}^{K} p(\mathbf{L}_{k} | \mathbf{x}_{1:t}, \mathbf{y}_{1:t})의 관계임을 이용하는 것이다.

18.2.3. Online parameter learning using recursive least squares

상태 공간 모델을 이용해 통계적 모델의 매개변수에 대한 온라인 베이지안 추론을 할 수 있다. 여기서는 선형 회귀에 대해 이를 수행해 보자. 기본 아이디어는 회귀 매개변수를 은닉 상태로 놓고 현 데이터를 관측 상태로 놓는 것이다. 즉, 사전분포를 p(\mathbf{\theta}) = \mathcal{N}(\mathbf{\theta} | \mathbf{\theta}_{0}, \mathbf{\Sigma}_{0})로 놓고, 은닉 상태 \mathbf{z}_{t} = \mathbf{\theta}가 바뀌지 않는다고 가정할 때 \mathbf{A}_{t} = \mathbf{I}, \mathbf{Q}_{t} = 0\mathbf{I}가 되므로 p(\mathbf{\theta}_{t} | \mathbf{\theta}_{t-1}) = \mathcal{N}(\mathbf{\theta}_{t} | \mathbf{\theta}_{t-1}, 0 \mathbf{I}) = \delta_{\mathbf{\theta}_{t-1}}(\mathbf{\theta}_{t})로 놓을 수 있다. (매개변수, 즉 은닉 상태가 시간에 따라 바뀔 경우 이를 동적 선형 모델이라 한다.) \mathbf{C}_{t} = \mathbf{x}_{t}^{T}, \mathbf{R}_{t} = \sigma^{2}로 둬서 관측 모델을 \mathcal{N}(\mathbf{y}_{t} | \mathbf{C}_{t}\mathbf{z}_{t}, \mathbf{R}_{t}) = \mathcal{N}(\mathbf{y}_{t} | \mathbf{x}_{t}^{T}\mathbf{\theta}_{t}, \sigma^{2})로 두자. 여기에 칼만 거르기를 수행하면 데이터가 스트리밍됨에 따라 사후상태를 업데이트할 수 있다. 이를 재귀적 최소 제곱법(RLS)이라 한다.

정확히는 사후평균의 칼만 업데이트는 칼만 증가 행렬 \mathbf{K}_{t} = \mathbf{\Sigma}^{t} \mathbf{C}_{t}^{T} \mathbf{R}_{t}^{-1}에 대해 \mathbf{\mu}_{t} = \mathbf{A}_{t} \mathbf{\mu}_{t-1} + \mathbf{K}_{t} (\mathbf{y}_{t} - \mathbf{C}_{t} \mathbf{A}_{t} \mathbf{\mu}_{t-1})로 놓을 수 있으므로, 매개변수의 업데이트는 \hat{\mathbf{\theta}}_{t} = \hat{\mathbf{\theta}}_{t-1} + \frac{1}{\sigma^{2}} \mathbf{\Sigma}_{t | t} (y_{t} - \mathbf{x}_{t}^{T} \hat{\mathbf{\theta}}_{t-1}) \mathbf{x}_{t}로 구할 수 있다. 이 때 \frac{1}{\sigma^{2}} \mathbf{\Sigma}_{t | t - 1} \simeq \eta_{t} \mathbf{I}로 근사하면 최소 평균 제곱(LMS) 알고리즘이 된다. 최소 평균 제곱 알고리즘에 대해서는 매개변수를 업데이트하는 비율 \eta_{t}를 특정해줘야 하고, 데이터에 대한 에포크가 여러 번일 수도 있다. 하지만 재귀적 평균 제곱법에서는 계단 크기를 자동적으로 조정해 주며, 데이터에 대한 한 번의 에포크로 최적의 사후분포로 수렴시킬 수 있다.

재귀적 선형 회귀의 예.

18.2.4. SSM for time series forecasting

상태 공간 모델은 시계열 예측에 매우 특화되어 있다. 상태 공간 모델은 은닉 상태를 추정하기 때문에 미래 상태 예측에는 도움이 안 될거라고 보일 수도 있지만 상태 공간 모델을 이용하면 데이터에 대한 생성적 모델을 잠재 과정에 대해 만들 수 있다. 이후 은닉 변수를 적분해 빼내면 관측 데이터의 사후예측분포를 구성할 수 있다. 모델이 선형-가우시안인 경우 이 과정들을 종합해 구조적 시계열 모델을 구성할 수 있다.

18.2.4.1. Local level model

가장 간단한 잠재 과정은 지역 층 모델로 다음의 형태를 갖는다:

y_{t} = a_{t} + \epsilon_{t}^{y}, \epsilon_{t}^{y} \sim \mathcal{N}(0, R)

a_{t} = a_{t-1} + \epsilon_{t}^{a}, \epsilon_{t-1}^{a} \sim \mathcal{N}(0, Q)

이 때 은닉 상태는 \mathbf{z}_{t} = a_{t}로 잡으며, 관측 상태 y_{t}가 은닉 상태 a_{t}에 분산이 R인 관측 노이즈에 더해진 것으로 생각한다. 또한, 은닉 상태층 a_{t}는 분산이 Q인 시스템 노이즈에 따라 변화해 나가는 것으로 생각한다.

18.2.4.2. Local linear trend

많은 시계열 분포는 위나 아래 방향으로의 선형 트렌드를 최소한 국소적으로는 가진다. 이는 은닉 층 a_{t}가 각 단계마다 b_{t}만큼 변하는 것으로 모델링할 수 있다.

y_{t} = a_{t} + \epsilon_{t}^{y}, \epsilon_{t}^{y} \sim \mathcal{N}(0, R)

a_{t} = a_{t-1} + b_{t-1} + \epsilon_{t}^{a}, \epsilon_{t}^{a} \sim \mathcal{N}(0, Q_{a})

b_{t} = b_{t-1} + \epsilon_{t}^{b}, \epsilon_{t}^{b} \sim \mathcal{N}(0, Q_{b})

Q_{b} = 0인 경우에는 b_{t} = b_{0}이므로 선의 기울기가 상수가 된다. Q_{a} = Q_{b} = 0인 경우에는 a_{t} = a_{t-1} + b_{0}t, 즉 \mathbb{E}[y_{t} | \mathbf{y}_{1 : t-1}] = a_{0} + tb_{0}이 된다. 이는 선형 트렌드 모델의 일반화이다.

18.2.4.3. Seasonality

많은 시계열은 주기적으로 진동하는데, 이는 S 단계의 사이클 동안 합해서 0이 되는 일련의 오프셋 항 c_{t}로 구성되는 잠재 과정들을 추가해 모델링할 수 있다:

c_{t} = -\sum_{s=1}^{S-1} c_{t-s} + \epsilon_{t}^{c}, \epsilon_{t}^{c} \sim \mathcal{N}(0, Q_{c})

시계열에 대한 상태 공간 모델.

18.2.4.4. ARMA models

고전적 시계열 예측은 자동 회귀 이동 평균법(ARMA) 모델에 기반한다. 이는 다음의 형태를 갖는다:

x_{t} = \sum_{i=1}^{p} \alpha_{i} x_{t-i} + \sum_{j=1}^{q} \beta_{j} w_{t-j} + v_{t}

이 때 v_{t}, w_{t} \sim \mathcal{N}(0, 1)은 독립적인 가우시안 노이즈 항이다. q = 0으로 놓으면 i < t - p에 대해 x_{t} \perp x_{i} | x_{t - 1 : t - p}인 순수 자동 회귀 모델이 된다. p = 0으로 놓으면 i < t - q에 대해 x_{t} \perp x_{i}인 순수 이동 평균 모델이 된다. 이 때 w_{t}는 인접한 시간 단계간 연관성을 유도하는 노이즈 변수들로 단기 상호 연관을 모델링한다.

자동 회귀 이동 평균법은 상태 공간 모델로 표현할 수도 있지만, 시계열에 대해서는 자동 회귀 이동 평균법보다 구조적 접근법이 이해하기 더 쉽고, 구조적 접근법은 매개변수가 시간에 따라 변화할 수 있도록 허용하기 때문에 모델이 비정적인 상황에 적응적일 수 있도록 한다.

18.3. Inference in LG-SSM

이 절에서는 선형-가우시안 상태 공간 모델의 추론을 온라인-오프라인 순서대로 다룬다.

18.3.1. The Kalman filtering algorithm

칼만 거르기는 선형 가우시안 상태 공간 모델에 대한 정확한 베이지안 거르기 알고리즘이다. 시점 t에서의 주변사후분포를 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t}) = \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t}, \mathbf{\Sigma}_{t})로 나타내면 모든 것이 가우시안이므로 예측과 업데이트를 모두 닫힌 형태로 표현할 수 있다. 이는 은닉 마르코프 모델과 비슷하다.

18.3.1.1. Prediction step

예측 단계는 다음과 같다.

p(\mathbf{z}_{t} | \mathbf{y}_{1 : t-1}, \mathbf{u}_{1 : t}) = \int \mathcal{N}(\mathbf{z}_{t} | \mathbf{A}_{t} \mathbf{z}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t}, \mathbf{Q}_{t}) \mathcal{N}(\mathbf{z}_{t-1} | \mathbf{\mu}_{t-1}, \mathbf{\Sigma}_{t-1}) d \mathbf{z}_{t-1} = \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t | t - 1}, \mathbf{\Sigma}_{t | t - 1})

\mathbf{\mu}_{t | t - 1} = \mathbf{A}_{t} \mathbf{\mu}_{t-1} + \mathbf{B}_{t} \mathbf{u}_{t}

\mathbf{\Sigma}_{t | t - 1} = \mathbf{A}_{t} \mathbf{\Sigma}_{t - 1} \mathbf{A}_{t}^{T} + \mathbf{Q}_{t}

18.3.1.2. Measurement step

측정 스텝은 베이즈 룰을 통해 다음과 같이 계산할 수 있다.

p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t}) \propto p(\mathbf{y}_{t} | \mathbf{z}_{t}, \mathbf{u}_{t}) p(\mathbf{z}_{t} | \mathbf{y}_{1 : t-1}, \mathbf{u}_{1 : t}) = \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t}, \mathbf{\Sigma}_{t})

\mathbf{\mu}_{t} = \mathbf{\mu}_{t | t - 1} + \mathbf{K}_{t} \mathbf{r}_{t}

\mathbf{\Sigma}_{t} = (\mathbf{I} - \mathbf{K}_{t} \mathbf{C}_{t}) \mathbf{\Sigma}_{t | t-1}

여기서 \hat{\mathbf{y}}_{t} = \mathbb{E}[\mathbf{y}_{t} | \mathbf{y}_{1 : t-1}, \mathbf{u}_{1 : t}] = \mathbf{C}_{t} \mathbf{\mu}_{t | t - 1} + \mathbf{D}_{t} \mathbf{u}_{t}에 대해 \mathbf{r}_{t} = \mathbf{y}_{t} - \hat{\mathbf{y}}_{t}잔여 오차이며, 조건공분산 \mathbf{S}_{t} = \mathrm{cov}[\mathbf{r}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}] = \mathbf{C}_{t} \mathbf{\Sigma}_{t | t-1} \mathbf{C}_{t}^{T} + \mathbf{R}_{t}에 대해 \mathbf{K}_{t} = \mathbf{\Sigma}_{t | t - 1} \mathbf{C}_{t}^{T} \mathbf{S}_{t}^{-1}칼만 증가 행렬이다. 이는 \mathbf{K}_{t} = (\mathbf{\Sigma}_{t | t - 1}^{-1} + \mathbf{C}^{T} \mathbf{R} \mathbf{C})^{-1} \mathbf{C}^{T} \mathbf{R}^{-1}로 나타낼 수 있다.

이 식의 의미를 생각해보자. 먼저 평균 업데이트가 \mathbf{\mu}_{t} = \mathbf{\mu}_{t | t - 1} + \mathbf{K}_{t} \mathbf{r}_{t}로 이루어진다는 점을 생각하자. 즉 새 평균은 이전 평균에다가 오차 신호에 가중치 역할을 하는 칼만 증가 행렬을 곱한 보정 항을 더한 것이다. \mathbf{C}_{t} = \mathbf{I}인 경우 \mathbf{K}_{t} = \mathbf{\Sigma}_{t | t - 1} \mathbf{S}_{t}^{-1}이 되어 사전분포의 공분산과 측정오차의 공분산 간의 비율이 된다. 사전분포가 강하거나 측정 도구의 노이즈가 심하다면 \lvert \mathbf{K}_{t} \rvert은 작을 것이고, 보정 항도 작을 것이다. 반대로 사전분포가 약하거나 측정 도구의 노이즈가 작다면 \lvert \mathbf{K}_{t} \rvert은 클 것이고, 보정 항도 커질 것이다.

18.3.1.3. Marginal likelihood

로그주변가능도도 계산할 수 있다.

\log p(\mathbf{y}_{1 : T} | \mathbf{u}_{1 : T}) = \sum_{t} \log p(\mathbf{y}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}) = \sum_{t} \log \mathcal{N}(\mathbf{y}_{t} | \mathbf{C}_{t} \mathbf{\mu}_{t | t- 1}, \mathbf{S}_{t})

18.3.1.4. Posterior predictive

한 단계 전 사후예측분포는 다음과 같으며, 시계열 예측에 유용하다.

p(\mathbf{y}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}) = \mathcal{N}(\mathbf{y}_{t} | \mathbf{C} \mathbf{\mu}_{t | t - 1}, \mathbf{C} \mathbf{\Sigma}_{t | t - 1}\mathbf{C}^{T} + \mathbf{R})

18.3.1.5. Computational issues

칼만 거르기 알고리즘의 병목은 첫째로는 칼만 증가 행렬을 계산할 때의 행렬 역연산 알고리즘으로, O(\lvert \mathbf{y}_{t} \rvert^{3}) 시간이 걸린다. 둘째로는 공분산 행렬을 계산하는 행렬곱으로 O(\lvert \mathbf{z}_{t} \rvert^{2}) 시간이 걸린다.

\lvert \mathbf{y}_{t} \rvert가 훨씬 큰 경우에는 \mathbf{K}_{t}를 미리 계산해 놓는 방법을 쓰는데, 이는 실제 관측 \mathbf{y}_{ 1 : t}에 의존하지 않기 때문이다. 이 때 \mathbf{\Sigma}_{t}를 재귀적으로 업데이트하는 식은 리카티 방정식이라 한다. 시간 불변 계 (\mathbf{\theta}_{t} = \mathbf{\theta})에서는 이는 고정점으로 수렴한다.

실제로는 정보량 거르기라는 알고리즘으로 \mathbf{\Lambda}_{t} = \mathbf{\Sigma}_{t}^{-1}\mathbf{\eta}_{t} = \mathbf{\Lambda}_{t} \mathbf{\mu}_{t}를 대신 업데이트하는 방법도 있다. 또 다른 방법은 제곱근 거르기로 촐스키 분해 \mathbf{\Sigma}_{t} = \mathbf{U}_{t} \mathbf{D}_{t} \mathbf{U}_{t}에 대해 분해 행렬들을 업데이트하는 방법도 있다.

18.3.1.6. Derivation

칼만 거르기 식들은 행렬 역변환 식으로 구할 수 있다.

18.3.2. The Kalman smoothing algorithm

칼만 거르기 알고리즘은 각 t에 대해 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t})을 계산한다. 이는 위치 추적 같은 온라인 추론 문제에 유용하다. 오프라인 추론 문제에서는 그 대신 p(\mathbf{z}_{t} | \mathbf{y}_{1 : T})을 계산한다. 과거와 미래 데이터 모두에 대해 조건을 걸게 되므로 불확실성이 대폭 감소한다. 이는 칼만 다듬질 또는 RTS 다듬질로 알려져 있다. 이는 은닉 마르코프 모델에서의 전방-후방 알고리즘과 비슷하다.

18.3.2.1. Algorithm

칼만 거르기 알고리즘은 그래프에서 메시지를 왼쪽에서 오른쪽으로 전달하는 것으로 볼 수 있다. 끝까지 가면 p(\mathbf{z}_{T} | \mathbf{y}_{1 : T})를 계산한 것이다. 그와 반대로 오른쪽에서 왼쪽으로 메시지를 전달한 뒤 그를 결합하는 것이 칼만 다듬질이다. 알고리즘은 다음과 같다.

p(\mathbf{z}_{t} | \mathbf{y}_{1 : T}) = \mathcal{N}(\mathbf{\mu}_{t | T}, \mathbf{\Sigma}_{t | T})

\mathbf{\mu}_{t | T} = \mathbf{\mu}_{t | t} + \mathbf{J}_{t}(\mathbf{\mu}_{t + 1 | T} - \mathbf{\mu}_{t + 1 | t})

\mathbf{\Sigma}_{t | T} = \mathbf{\Sigma}_{t | t} + \mathbf{J}_{t}(\mathbf{\Sigma}_{t + 1 | T} - \mathbf{\Sigma}_{t + 1 | t}) \mathbf{\Sigma}_{t}^{T}

\mathbf{J}_{t} = \mathbf{\Sigma}_{t | t} \mathbf{A}_{t + 1}^{T} \mathbf{\Sigma}_{t + 1 | t}^{-1})

이 때 \mathbf{J}_{t}는 후방 칼만 증가 행렬이다. 이는 칼만 거르기 알고리즘에서 구한 \mathbf{\mu}_{T | T}, \mathbf{\Sigma}_{T | T}로 초기화할 수 있다. 후방 알고리즘은 계산 도중 데이터 \mathbf{y}_{1 : T}를 알 필요가 없으므로 메모리를 절약할 수 있다는 사실을 눈여겨보자.

18.3.2.2. Derivation

식의 유도는 마르코프 특성과 가우시안 결합 분포의 성질을 사용한다.

18.3.2.3. Comparison to the forwards-backwards algorithm for HMMs

선형 동역학계에서의 전방-후방 알고리즘은 항상 표준화된 데이터에만 사용하고 후방 알고리즘은 전방 알고리즘의 결과에 의존한다. 그와 반대로 은닉 마르코프 모델의 경우 후방 알고리즘이 전방 알고리즘에 의존하지 않는다. 이에 착안해 칼만 다듬질을 다음과 같이 조금 더 쉽게 계산할 수 있다.

\frac{p(\mathbf{z}_{t + 1} | \mathbf{y}_{1 : T})}{p(\mathbf{z}_{t + 1} | \mathbf{y}_{1 : t})} \propto p(\mathbf{y}_{t + 1 : T} | \mathbf{z}_{t + 1})

그러므로 은닉 마르코프 모델처럼 후방 메시지는 전방 메시지에 상관없이 미래 데이터에 대한 조건부 가능도만으로 계산할 수 있다. 하지만 이 접근법도 단점이 있는데, 첫째로는 원본 관측 서열 \mathbf{y}_{t + 1 : T}를 알아야 한다는 것이고 둘째로는 이것은 사후분포가 아닌 가능도이기 때문에 공분산이 양의 정칙인 가우시안으로 나타내는 것이 항상 가능한 것은 아니다. 셋째로는 정확한 추론이 불가능할 때에는 이러한 후방 가능도 항보다는 다듬질된 분포를 근사하는 방법이 더 낫다는 것이다.

다른 변형으로는 이중 거르기 다듬질로 전방 알고리즘에서는 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t})을 계산하고 후방 알고리즘에서는 걸러진 사후분포 p(\mathbf{z}_{t} | \mathbf{y}_{t + 1 : T})을 계산한 뒤 이를 조합해 p(\mathbf{z}_{t} | \mathbf{y}_{1 : T})을 계산하는 것이다.

18.4. Learning for LG-SSM

선형 가우시안 상태 공간 모델의 매개변수를 추정하는 문제는 계 식별로 불리기도 한다. 이 때 관측 행렬 \mathbf{C}와 전이 행렬 \mathbf{A}는 알려져 있고 고정되어 있으며, 노이즈 공분산 \mathbf{Q}\mathbf{R}만 알면 된다. 이들은 오프라인으로 추정할 수도 있지만, 재귀적으로 사후분포 p(\mathbf{z}_{t}, \mathbf{R}, \mathbf{Q} | \mathbf{y}_{ 1 : t})을 정확히 계산해내는 것도 가능하다.

18.4.1. Identifiability and numerical stability

더 일반적인 상황에서는 \mathbf{A}\mathbf{C}도 학습해야 한다. 이 경우 일반성을 잃지 않고 \mathbf{Q} = \mathbf{I}로 둘 수 있으며 (노이즈를 전이 행렬에 끼워넣으면 되므로) \mathbf{R}도 대각 행렬로 둘 수 있다. 이는 매개변수의 자유도를 낮춰 수치적 안정성을 증가시킨다.

또 다른 유용한 제한 조건은 \mathbf{A}의 고유값에 제한을 두는 것이다. 계 노이즈가 없는 경우 은닉 상태 \mathbf{z}_{t} = \mathbf{A}^{t} \mathbf{z}_{1} = \mathbf{U} \mathbf{\Lambda}_{t} \mathbf{U}^{-1} \mathbf{z}_{1}가 되므로 (고유값 분해를 사용) \mathbf{\lambda}_{i} > 1이 되면 시간에 따라 \mathbf{z}_{t}의 크기가 폭발하게 된다. 따라서 모든 고유값의 크기를 1 미만으로 제한하는 것이 중요하다. 물론 이 경우엔 큰 t에 대해 \mathbb{E}[\mathbf{z}_{t}] = \mathbf{0}이 된다는 문제가 있지만 노이즈를 더하면 상태가 0이 아니게 되므로 이는 해결 가능하다. 아래에서는 이 매개변수들을 어떻게 추정하는지를 다룬다.

18.4.2. Training with fully observed data

은닉 상태 서열을 관측하면 \mathbf{z}_{t - 1} \to \mathbf{z}_{t}, \mathbf{z}_{t} \to \mathbf{y}_{t}의 다변수 선형 회귀 문제를 품으로써 최대가능도추정을 구할 수 있다. 즉, J(\mathbf{A}) = \sum_{t=1}^{2} (\mathbf{z}_{t} - \mathbf{A} \mathbf{z}_{t-1})^{2}를 품으로써 \mathbf{A}와 시스템 노이즈 공분산 \mathbf{Q}를 추정할 수 있고, 비슷한 방법으로 \mathbf{C}\mathbf{R}도 추정할 수 있다.

18.4.3. EM for LG-SSM

출력 서열만 관측할 수 있다면, 매개변수에 대한 최대가능도추정이나 최대사후확률추정은 기대값 최대화 알고리즘으로 이루어진다. 이는 은닉 마르코프 모델에서의 바움-웰치 알고리즘과 비슷하지만, E 단계에서 전방-후방 알고리즘 대신 칼만 다듬질 알고리즘을 사용하므로 M 단계에서의 계산도 달라진다는 점이 다르다.

18.4.4. Subspace methods

기대값 최대화 알고리즘은 초기 매개변수에 의존하기 때문에 항상 만족스러운 결과를 낳지는 않는다. 부분공간법이란 접근법을 사용해서 이를 피할 수 있다. 일단 관측 노이즈와 시스템 노이즈가 없음을 가정하면 \mathbf{z}_{t} = \mathbf{A} \mathbf{z}_{t-1}이고 \mathbf{y}_{t} = \mathbf{C} \mathbf{z}_{t} = \mathbf{C} \mathbf{A}^{t-1} \mathbf{z}_{1}이 되므로 모든 관측은 \mathrm{dim}(\mathbf{z}_{t})차원의 선형 부분공간으로부터 생성된다. 이 부분공간은 주성분 분석으로 알아낼 수 있다. \mathbf{z}_{t}를 추정한 뒤에는 그것이 완전히 관측된 것처럼 모델을 피팅할 수 있다. 이후에는 이것을 기대값 최대화 알고리즘에 쓰던가 하면 된다.

18.4.5. Bayesian methods for “fitting” LG-SSMs

기대값 최대화 알고리즘의 여러 오프라인 베이지안 버전 등이 존재하며 그 예로는 변량 베이즈 기대값 최대화, 막힌 깁스 샘플링 등이 있다. 베이지안 법은 온라인 학습에도 쓸 수 있다. 하지만, 상태 공간 모델 매개변수를 상태 공간에 더하면 일반적으로 그 때부터 모델은 선형 가우시안이 아니게 된다. 따라서 근사적인 온라인 추정 방법을 써야만 한다.

18.5. Approximate online inference for non-linear, non-Gaussian SSMs

모델은 비선형으로 움직일 수도 있다. 선형으로 움직이더라도 모르는 매개변수를 상태 공간에 넣으면 모델은 비선형이 된다. 또한, 이상치 등이나 매개변수의 비선형성 등으로 인해 가우시안이 아닌 노이즈는 매우 흔하다. 이를 위해서는 근사적인 추론을 수행한다.

X가 가우시안이고 f가 비선형 함수이고 모델이 Y = f(X)의 형태로 나타난다면 p(Y)를 가우시안으로 근사하는 방법이 크게 두 가지 있다. 첫째는 f의 일차 근사를 사용하는 것이고, 둘째는 정확한 f를 사용하되 f(X)를 가우시안의 함수공간에 모멘트 매칭으로 사영시키는 것이다.

18.5.1. Extended Kalman filter (EKF)

노이즈가 가우시안이고 전이 모델 g와 관측 모델 h는 비선형이되 미분가능한 다음의 비선형 모델을 가정하자.

\mathbf{z}_{t} = g(\mathbf{u}_{t}, \mathbf{z}_{t - 1}) + \mathcal{N}(\mathbf{0}, \mathbf{Q}_{t})

\mathbf{y}_{t} = h(\mathbf{z}_{t}) + \mathcal{N}(\mathbf{0}, \mathbf{R}_{t})

여기서는 사후분포를 단일 가우시안으로 근사할 수 있는 경우에 집중해보도록 하자. 확장된 칼만 거르기(EKF)는 이러한 형태의 비선형 가우시안 동역학계에 쓰일 수 있다. 기본 발상은 g와 h를 일차 테일러 급수 전개 근사를 이용해 이전 상태 추정값에 대한 선형 함수로 근사하는 것이다. 이후 표준 칼만 거르기 식들을 적용하면 된다. (노이즈 \mathbf{Q}, \mathbf{R}은 변하지 않는다고 가정한다.) 즉, 정적인 비선형 동역학계를 비-정적인 선형 동역학계로 근사하는 것이다.

\mathbf{H}_{t}\mathbf{h}의 자코비안을 사전최빈값에서 구한 값이라 할 때, 측정 모델은 다음과 같이 근사된다:

p(\mathbf{y}_{t} | \mathbf{z}_{t}) \simeq \mathcal{N}(\mathbf{y}_{t} | \mathbf{h}(\mathbf{\mu}_{t | t - 1}) + \mathbf{H}_{t} (\mathbf{y}_{t} - \mathbf{\mu}_{t | t - 1}), \mathbf{R}_{t})

\mathbf{H}_{t} = \mathbf{H} |_{\mathbf{z} = \mathbf{\mu}_{t | t - 1}}

\mathbf{G}_{t}\mathbf{g}의 자코비안을 사전최빈값에서 구한 값이라 할 때, 계 모델은 다음과 같이 근사된다:

p(\mathbf{z}_{t} | \mathbf{z}_{t - 1}, \mathbf{u}_{t}) \simeq \mathcal{N}(\mathbf{z}_{t} | \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1}) + \mathbf{G}_{t}(\mathbf{z}_{t - 1} - \mathbf{\mu}_{t - 1}), \mathbf{Q}_{t})

\mathbf{G}_{t} = \mathbf{G}(\mathbf{u}_{t}) |_{\mathbf{z} = \mathbf{\mu}_{t | t - 1}}

이를 통해 칼만 거르기를 적용하면 사후분포를 다음과 같이 구할 수 있다.

\mathbf{\mu}_{t | t - 1} = \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1})

\mathbf{V}_{t | t - 1} = \mathbf{G}_{t} \mathbf{V}_{t - 1} \mathbf{G}_{t}^{T} + \mathbf{Q}_{t}

\mathbf{K}_{t} = \mathbf{V}_{t | t - 1}\mathbf{H}_{t}^{T} (\mathbf{H}_{t} \mathbf{V}_{t | t - 1}\mathbf{H}_{t}^{T} + \mathbf{R}_{t})^{-1}

\mathbf{\mu}_{t} = \mathbf{\mu}_{t | t - 1} + \mathbf{K}_{t}(\mathbf{y}_{t} - \mathbf{h}(\mathbf{\mu}_{t | t - 1}))

\mathbf{V}_{t} = (\mathbf{I} - \mathbf{K}_{t} \mathbf{H}_{t}) \mathbf{V}_{t | t - 1}

일반 칼만 거르기와 다른 점은 상태를 예측할 때 \mathbf{A}_{t} \mathbf{\mu}_{t - 1} + \mathbf{B}_{t} \mathbf{\mu}_{t} 대신 \mathbf{g}(\mathbf{u}_{t}, \mathbf{\mu}_{t - 1})을 사용한다는 것이고 측정을 업데이트할 때 \mathbf{C}_{t} \mathbf{\mu}_{t | t - 1} 대신 h(\mathbf{\mu}_{t | t - 1})을 사용한다는 것이다. \mathbf{\mu}_{t | t - 1} 대신 \mathbf{\mu}_{t}를 중심으로 선형 근사를 수행해 성능을 향상시킬 수 있는데 이를 반복 확장된 칼만 거르기라고 한다. 더 느리지만 더 좋은 성능을 낸다.

확장된 칼만 거르기가 잘 동작하지 않는 두 가지 경우가 있는데 첫째는 사전 공분산이 클 경우 사전분포가 넓어지기 때문에 선형 근사가 평균에서 많이 떨어진 곳으로 잡히게 된다. 둘째는 현재 평균 근처에서 함수가 심하게 비선형적일 경우 당연히 선형 근사가 잘 동작하지 않는다. 이는 무취 칼만 거르기라는 알고리즘으로 보정할 수 있다.

18.5.2. Unscented Kalman filter (UKF)

무취 칼만 거르기(UKF)는 확장 칼만 거르기를 개선한 것이다. 기본 발상은 함수를 근사하는 것보다 가우시안을 근사하는 게 더 쉬우므로 선형 근사를 수행하는 대신 선택된 점들인 시그마 점들에 가우시안을 전달해 변환된 점들에 대해 가우시안을 피팅하는 것이다. 이를 무취 변환이라 한다.

무취 칼만 거르기는 기본적으로 무취 변환을 2번 수행한다. 한 번은 계 모델 \mathbf{g}를 통해 근사하는 것이고 하나는 측정 모델 \mathbf{h}를 통해 근사하는 것이다. 무취 칼만 거르기와 확장 칼만 거르기는 모두 잠재 상태 공간의 크기가 d일 때 단계당 O(d^{3})의 시간복잡도를 갖는다. 하지만 확장 칼만 거르기는 일차 근사인 데 비해 무취 칼만 거르기는 최소 2차 이상의 정확도를 갖는다. (두 방법 모두 더 고차에 대해 해석적으로 확장될 수 있기는 하다) 또한, 무취 변환은 도함수나 자코비안들에 대한 해석적 계산을 필요로 하지 않는다. (도함수 불요 거르기) 그래서 더 구현하기 쉽고 더 넓은 범위에 적용할 수 있다.

18.5.2.1. The unscented transform

무취 변환을 설명하면 다음과 같다. p(\mathbf{x}) = \mathcal{N}(\mathbf{x} | \mathbf{\mu}, \mathbf{\Sigma})이라 가정하고 \mathbf{y} = \mathbf{f}(\mathbf{x})일 때 p(\mathbf{y})를 추정하는 것을 생각해 보자. 무취 변환은 다음과 같이 수행된다. 우선 2d + 1개의 시그마 점들을 \mathbf{x} = (\mathbf{\mu}, \{ \mathbf{\mu} + \sqrt{(d + \lambda) \mathbf{\Sigma}}_{: i}\}_{i = 1}^{d}, \{ \mathbf{\mu} - \sqrt{(d + \lambda) \mathbf{\Sigma}}_{: i}\}_{i = 1}^{d}) 로 잡는다. (\lambda = \alpha^{2} (d + \kappa) - d은 변환 인자이고 \mathbf{M}_{:i}\mathbf{M}의 i번째 열)

이 시그마 점들은 비선형 함수를 통해 전파되어 \mathbf{y}_{i} =f(\mathbf{x}_{i})로 변환되고 평균과 공분산은 다음과 같이 가중치와 함께 계산된다.

\mathbf{\mu}_{y} = \sum_{i=0}^{2d} w_{m, i} \mathbf{y}_{i}

\mathbf{\Sigma}_{y} = \sum_{i=0}^{2d} w_{c, i} (\mathbf{y}_{i} - \mathbf{\mu}_{y})(\mathbf{y}_{i} - \mathbf{\mu}_{y})^{T}

w_{0, m} = \frac{\lambda}{d + \lambda}

w_{0, c} = \frac{d}{d + \lambda} + 1 - \alpha^{2} + \beta

w_{i, m} = w_{i, c} = \frac{1}{2(d + \lambda)}

이 때 \alpha, \beta, \kappa의 최적의 값은 문제에 따라 달라진다. d = 1일 때는 \alpha = 1, \beta = 0, \kappa = 2이다.

18.5.2.2. The UKF algorithm

무취 칼만 거르기 알고리즘은 무취 변환을 두 번 적용한 것이다. 한 번은 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t})를 계산할 때이고 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t})을 계산할 때이다.

첫 번쨰 단계는 사후예측분포 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t - 1}, \mathbf{u}_{1 : t}) \simeq \mathcal{N}(\mathbf{z}_{t} | \bar{\mathbf{\mu}}_{t}, \bar{\mathbf{\Sigma}}_{t})를 이전 믿음 상태 \mathcal{N}(\mathbf{z}_{t-1}  \mathbf{\mu}_{t - 1}, \mathbf{\Sigma}_{t - 1})를 계 모델 \mathbf{g}를 통해 다음과 같이 전파해 근사하는 것이다:

\mathbf{z}_{0, t - 1} = (\mathbf{\mu}_{t - 1}, \{ \mathbf{\mu}_{t - 1} + \sqrt{ (d + \lambda)\mathbf{\Sigma}_{t-1}}_{: i}\}_{i = 1}^{d}, \{ \mathbf{\mu}_{t - 1} - \sqrt{ (d + \lambda)\mathbf{\Sigma}_{t-1}}_{: i}\}_{i = 1}^{d})

\bar{\mathbf{z}}_{t, \ast i} = \mathbf{g}(\mathbf{u}_{t}, \mathbf{z}_{t-1, 0i})

\bar{\mathbf{\mu}}_{t} = \sum_{i=0}^{2d} w_{m, i} \bar{\mathbf{z}}_{t, \ast i}

\bar{\mathbf{\Sigma}}_{t} = \sum_{i=0}^{2d} w_{c, i} (\bar{\mathbf{z}}_{t, \ast i} - \bar{\mathbf{\mu}}_{t})(\bar{\mathbf{z}}_{t, \ast i} - \bar{\mathbf{\mu}}_{t})^{T} + \mathbf{Q}_{t}

두 번째 단계는 가능도 p(\mathbf{y}_{t} | \mathbf{z}_{t}) \simeq \mathcal{N}(\mathbf{y}_{t} | \hat{\mathbf{y}}_{t}, \mathbf{S}_{t})을 사전분포 \mathcal{N}(\mathbf{z}_{t} | \bar{\mathbf{\mu}}_{t}, \bar{\mathbf{\Sigma}}_{t})을 관측 모델 \mathbf{h}을 통해 다음과 같이 전파해 근사하는 것이다:

\bar{\mathbf{z}}_{0, t} = (\bar{\mathbf{\mu}}_{t}, \{ \bar{\mathbf{\mu}}_{t} + \sqrt{ (d + \lambda)\bar{\mathbf{\Sigma}}_{t}}_{: i}\}_{i = 1}^{d}, \{ \bar{\mathbf{\mu}}_{t} - \sqrt{ (d + \lambda)\bar{\mathbf{\Sigma}}_{t}}_{: i}\}_{i = 1}^{d})

\bar{\mathbf{y}}_{t, \ast i} = \mathbf{h}(\bar{\mathbf{z}}_{t, 0i})

\hat{\mathbf{y}}_{t} = \sum_{i=0}^{2d} w_{m, i} \bar{\mathbf{y}}_{t, \ast i}

\mathbf{S}_{t} = \sum_{i=0}^{2d} w_{c, i} (\bar{\mathbf{y}}_{t, \ast i} - \hat{\mathbf{y}}_{t})(\bar{\mathbf{y}}_{t, \ast i} - \hat{\mathbf{y}}_{t})^{T} + \mathbf{R}_{t}

이제 베이즈 정리를 통해 사후분포 p(\mathbf{z}_{t} | \mathbf{y}_{1 : t}, \mathbf{u}_{1 : t}) \simeq \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t}, \mathbf{\Sigma}_{t})를 근사할 수 있다.

\bar{\mathbf{\Sigma}}_{t, z, y} = \sum_{i=0}^{2d} w_{c, i}(\bar{\mathbf{z}}_{t, \ast i} - \bar{\mathbf{\mu}}_{t})(\bar{\mathbf{y}}_{t, \ast i} - \hat{\mathbf{y}}_{t})^{T}

\mathbf{K}_{t} = \bar{\mathbf{\Sigma}}_{t, z, y} \mathbf{S}_{t}^{-1}

\mathbf{\mu}_{t} = \bar{\mathbf{\mu}}_{t} + \mathbf{K}_{t} (\mathbf{y}_{t} - \hat{\mathbf{y}}_{t})

\mathbf{\Sigma}_{t} = \bar{\mathbf{\Sigma}}_{t} - \mathbf{K}_{t} \mathbf{S}_{t} \mathbf{K}_{t}^{T}

18.5.3. Assumed density filtering (ADF)

업데이트는 정확히 하더라도 사후분포를 가우시안과 같은 익숙한 분포로 근사하고 싶을 때가 있다. 추론하고 싶은 미지 변수를 \mathbf{\theta}_{t}라 두고, \mathcal{Q}를 추적 가능한 분포의 집합이라 하자. 이 때 사전분포를 q_{t-1} \in \mathcal{Q}에 대해 q_{t-1}(\mathbf{\theta}_{t-1}) \simeq p(\mathbf{\theta}_{t-1} | \mathbf{y}_{1 : t - 1})로 근사할 수 있을 때 이를 통해 사후분포를 다음과 같이 근사할 수 있다:

\hat{p}(\mathbf{\theta})_{t} = \frac{1}{Z_{t}} p(\mathbf{y}_{t} | \mathbf{\theta}_{t}) q_{t | t-1} (\mathbf{\theta}_{t})

Z_{t} = \int p(\mathbf{y}_{t} | \mathbf{\theta}_{t}) q_{t | t-1} (\mathbf{\theta}_{t}) d \mathbf{\theta}_{t}

q_{t | t-1} (\mathbf{\theta}_{t}) = \int p(\mathbf{\theta}_{t} | \mathbf{\theta}_{t - 1})q_{t-1}(\mathbf{\theta}_{t-1}) d \mathbf{\theta}_{t-1}

사전분포에 적절한 제한을 걸면 한 단계 업데이트 식은 대개 계산 가능하다. 하지만, 결과 사후분포는 계산 가능하지 않은 경우가 많다. 따라서 업데이트 이후에는 계산 가능한 함수로 이를 근사하는데, q(\mathbf{\theta}_{t}) = \mathrm{argmin}_{q \in \mathcal{Q}} \mathbb{KL} (\hat{p}(\mathbf{\theta}_{t}) || q(\mathbf{\theta}_{t}) )을 계산한다. 이는 \hat{p}를 계산 가능한 함수의 공간에 사영시키는 것으로 생각할 수 있다. 이를 예측-업데이트-사영 사이클이라 하며, 가정된 밀도 거르기(ADF)이라 불린다. q가 지수족이면 KL 분산의 최소화는 모멘트 짝맞춤으로 수행될 수 있다.

18.5.3.1. Boyen-Koller algorithm for online inference in DBNs

이산 상태의 동적 베이즈 망에 대해 추론을 한다면, \theta_{tj}가 j번째 은닉 변수의 시간 t에서의 값이라고 할 때 정확한 사후분포 p(\mathbf{\theta}_{t})은 꼬임 문제 때문에 정확하게 계산할 수 없어진다. 이 때 p에 대한 근사가 q(\mathbf{\theta}_{t}) = \prod_{j=1}^{D} \mathrm{Cat}(\theta_{t, j} | \mathbf{\pi}_{t, j}), \pi_{tjk} = q(\theta_{t, j} = k)로 완전히 인수분해된다고 하자. 이 경우 모멘트 짝맞춤 연산은 \pi_{tjk} = \hat{p}(\theta_{t, j} = k)이 된다. 이는 인수분해된 사전분포를 사용하기 전 예측-업데이트 단계를 수행하고 사후주변분포를 계산함으로써 계산할 수 있다. 이는 보옌-콜러 알고리즘이라 한다.

18.5.3.2. Gaussian approximation for online inference in GLMs

이제 q(\mathbf{\theta}_{t}) = \prod_{j=1}^{D} \mathcal{N}(\theta_{t, j} | \mu_{t, j}, \tau_{t, j})이라 하자. 이 때 사후분포의 계산 가능한 근사 분포에 대한 매개변수의 최적화된 값은 \mu_{t, j} = \mathbb{E}_{\hat{p}} [\theta_{t, j}], \tau_{t, j} = \mathrm{var}_{\hat{p}}[\theta_{t, j}]로 주어진다. 이는 온라인 추론에 유용하다.

이진 로지스틱 회귀 문제에 위를 적용해 보자.

p(y_{t} | \mathbf{x}_{t}, \mathbf{\theta}_{t}) = \mathrm{Ber}(y_{t} | \mathrm{sigm}(\mathbf{x}_{t}^{T} \mathbf{\theta}_{t}))

p(\mathbf{\theta}_{t} | \mathbf{\theta}_{t-1}) = \mathcal{N}(\mathbf{\theta}_{t} | \mathbf{\theta}_{t-1}, \sigma^{2} \mathbf{I})

모델은 위의 형태를 가진다. (\sigma^{2}는 매개변수가 시간에 따라 변화하도록 하는 프로세스 노이즈로, 0이 될 수 있다) 이제 계산 가능한 사전분포의 인수분해를 q_{t-1}(\mathbf{\theta}_{t-1}) = \prod_{j} \mathcal{N}(\theta_{t-1, j} | \mu_{t-1, j}, \tau_{t-1, j})라 하자.

이 때 한 단계 사후예측분포 q_{t | t - 1}(\mathbf{\theta}_{t})은 선형 가우시안 업데이트를 통해 구할 수 있다. 그러므로 남은 것은 측정 업데이트 부분이다.

s_{t} = \mathbf{\theta}_{t}^{T} \mathbf{x}_{t}로 두고 q_{t | t - 1}(\mathbf{\theta}_{t}) = \prod_{j} \mathcal{N}(\theta_{t, j} | \mu_{t | t-1, j}, \tau_{t | t-1, j})로 둔다면, s_{t}의 사후예측분포는 다음과 같이 계산할 수 있다.

q_{t | t-1}(s_{t}) = \mathcal{N}(s_{t} | m_{t | t-1}, v_{t | t-1})

m_{t | t-1} = \sum_{j} x_{t, j} \mu_{t | t - 1, j}

v_{t | t-1} = \sum_{j} x_{t, j}^{2} \tau_{t | t - 1, j}

사후분포는 다음과 같다.

q_{t}(s_{t}) = \mathcal{N}(s_{t} | m_{t}, v_{t})

m_{t} = \int s_{t} \frac{1}{Z_{t}} p(y_{t} | s_{t}) q_{t | t-1}(s_{t}) ds_{t}

v_{t} = \int s_{t}^{2} \frac{1}{Z_{t}} p(y_{t} | s_{t}) q_{t | t-1}(s_{t}) ds_{t} - m_{t}^{2}

Z_{t} = \int p(y_{t} | s_{t}) q_{t | t-1}(s_{t}) ds_{t}

이 때 p(y_{t} | s_{t}) = \mathrm{Ber}(y_{t} | s_{t})이 된다. 이 적분은 1차원이므로 가우시안 구적법으로 계산될 수 있다. 이는 무취 칼만 거르기 알고리즘에서의 한 단계와 같다.

q(s_{t})을 추론한 뒤에는 q(\mathbf{\theta} | s_{t})를 계산해야 한다. \delta_{m} = m_{t} - m_{t | t-1}, \delta_{v} = v_{t} - v_{t | t-1}로 정의하면 모델 변수에 대한 인수분해된 사후분포는 다음과 같음을 보일 수 있다.

q(\theta_{t, j}) = \mathcal{N}(\theta_{t, j} | \mu_{t, j}, \tau_{t, j})

\mu_{t, j} = \mu_{t | t-1, j} + a_{j} \delta_{m}

\tau_{t, j} = \tau_{t | t-1, j} + a_{j}^{2} \delta_{v}

a_{j} = \frac{x_{t, j} \tau_{t | t-1, j}}{\sum_{k} x_{t, k} \tau_{t | t-1, k}}

그러므로 큰 크기를 가졌거나 (\lvert x_{t, j} \rvert이 큰) 불확실성이 큰 (\tau_{t | t-1, j}이 큰) 입력에 대응하는 매개변수가 가장 큰 빈도로 업데이트됨을 알 수 있다.

프로빗 가능도를 이용한 변형은 측정 업데이트는 닫힌 형태로 가능하다. 이 경우 알고리즘은 단계별로 O(D)의 시간밖에 들지 않으므로 많은 매개변수에 대해서도 수행 가능하다. 또한 온라인 알고리즘이므로 큰 데이터셋을 다룰 수도 있다.

18.6. Hybrid discrete/continuous SSMs

많은 계는 이산적/연속적 은닉 변수를 모두 갖고 있으며, 이를 혼성 계라 불린다. 특수한 예로 은닉 마르코프 모델과 선형 상태 공간 모델을 결합한 변환 선형 동역학계(SLDS), 점프 마르코프 선형 계(JMLS), 변환 상태 공간 모델(SSSM) 등이 있다. 구체적으로, 이산적인 잠재 변수 q_{t} \in \{1, \cdots, K\}, 연속적인 잠재 변수 \mathbf{z}_{t} \in \mathbb{R}^{L}, 연속적인 관측 변수 \mathbf{y}_{t} \in \mathbf{R}^{D}, 연속적인 관측 입력/제어 \mathbf{u}_{t} \in \mathbf{R}^{U}이 있다고 하자. 이 때 연속적인 변수들은 선형 가우시안 조건부 확률밀도함수 분포를 따르며 이산적인 상태에 조건이 걸린다고 가정한다.

p(q_{t} = k | q_{t-1} = j, \mathbf{\theta}) = A_{ij}

p(\mathbf{z}_{t} | \mathbf{z}_{t-1}, q_{t} = k, \mathbf{u}_{t}, \mathbf{\theta}) = \mathcal{N}(\mathbf{z}_{t} | \mathbf{A}_{k} \mathbf{z}_{t-1} + \mathbf{B}_{k} \mathbf{u}_{t}, \mathbf{Q}_{k})

p(\mathbf{y}_{t} | \mathbf{z}_{t}, q_{t} = k, \mathbf{u}_{t}, \mathbf{\theta}) = \mathcal{N}(\mathbf{y}_{t} | \mathbf{C}_{k} \mathbf{z}_{t} + \mathbf{D}_{k} \mathbf{u}_{t}, \mathbf{R}_{k})

18.6.1. Inference

혼성 모델을 통한 추론은 계산 가능하지 않다. 왜 그럴까? q_{t}가 이진 변수라고 가정하면 p(\mathbf{z}_{1} | \mathbf{y}_{1})은 2개의 가우시안의 혼합이다. p(\mathbf{z}_{2} | \mathbf{y}_{1})은 4개의 가우시안의 혼합이 될 것이다. 계속해나가다 보면 8개, 16개, …로 늘어서 최빈값의 수가 지수적으로 증가해버릴 것이다.

그래서 이 모델에 대해서는 여러 근사적인 추론 법이 존재한다.

  • 이산 트리에서 확률이 낮은 부분을 가지치기한다. 이는 다중 가설 추적법의 기본이다.
  • 몬테 카를로법을 쓴다. 이산적인 경로들을 샘플링한 뒤 이에 조건을 걸고 연속적인 변수에 해석적인 거르기를 적용한다.
  • 가정된 밀도 거르기법을 쓴다. 지수함수적으로 개수가 큰 가우시안의 혼합을 더 작은 가우시안의 혼합으로 근사한다.

18.6.1.1. A Gaussian sum filter for switching SSMs

가우시안 합 거르기 방법은 각 단계에서의 믿음 상태를 K개의 가우시안의 혼합으로 근사한다. 이는 K개의 칼만 거르기를 병렬로 수행해서 구현할 수 있으며, 변환 상태 공간 모델에 잘 맞는 방법이다. 이 방법 중 하나로 이차 일반화된 유사 베이즈 거르기 (GPB2)가 있다. 상세는 다음과 같다.

사전믿음상태 b_{t - 1}을 K개의 가우시안의 혼합

b_{t-1, i} = p(\mathbf{z}_{t - 1}, q_{t - 1} = i | \mathbf{y}_{1 : t - 1}) = \pi_{t - 1, i} \mathcal{N}(\mathbf{z}_{t - 1} | \mathbf{\mu}_{t - 1, i}, \mathbf{\Sigma}_{t - 1, i})

으로 가정한다. 그리고 이를 K개의 다른 선형 모델로 전파시켜 다음을 얻는다:

b_{t, ij} = p(\mathbf{z}_{t}, q_{t - 1} = i, q_{t} = j | \mathbf{y}_{1 : t}) = \pi_{t, ij} \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t, ij}, \mathbf{\Sigma}_{t, ij}), \pi_{t, ij} = \pi_{t-1, i} p(q_{t} = j | q_{t - 1} = i).

이 때, j의 각 값에 대해 K개의 가우시안 혼합을 단일 혼합으로 근사시켜 다음을 얻을 수 있다:

b_{t, j} = p(\mathbf{z}_{t}, q_{t} = j | \mathbf{y}_{1 : t}) = \pi_{t, j} \mathcal{N}(\mathbf{z}_{t} | \mathbf{\mu}_{t, j}, \mathbf{\Sigma}_{t, j})

가우시안 혼합을 단일 가우시안으로 근사하는 최적의 법은 p(\mathbf{z}) = \sum_{k} \pi_{k} \mathcal{N}(\mathbf{z} | \mathbf{\mu}_{k}, \mathbf{\Sigma}_{k}), q(\mathbf{z}) = \mathcal{N}(\mathbf{z} | \mathbf{mu}, \mathbf{\Sigma})에 대해 q = \mathrm{argmin}_{q} \mathbb{KL} (q || p)을 구하는 것이다. 이는 모멘트 짝맞춤으로 할 수 있다. 즉,

\mathbf{\mu} = \mathbb{E}[\mathbf{z}] = \sum_{k} \pi_{k} \mathbf{\mu}_{k}

\mathbf{\Sigma} = \mathrm{cov}[\mathbf{z}] = \sum_{k} \pi_{k} (\mathbf{\Sigma}_{k} + (\mathbf{\mu}_{k} - \mathbf{\mu})(\mathbf{\mu}_{k} - \mathbf{\mu})^{T})

그래프 모델에서 이는 약한 주변화라 한다. 이 식들을 모델에 적용시키면 b_{t, ij}에서 b_{t, j}를 다음과 같이 구할 수 있다.

\pi_{t, j} = \sum_{i} \pi_{t, ij}

\pi_{t, j | i} = \frac{\pi_{t, ij}}{\sum_{k} \pi_{t, ik}}

\mathbf{\mu}_{t, j} = \sum_{i} \pi_{t, j | i} \mathbf{\mu}_{t, ij}

\mathbf{\Sigma}_{t, j} = \sum_{i} \pi_{t, j | i} (\mathbf{\Sigma}_{t, ij} + (\mathbf{\mu}_{t, ij} - \mathbf{\mu}_{t, j})(\mathbf{\mu}_{t, ij} - \mathbf{\mu}_{t, j})^{T})

이는 각 단계마다 K^{2}개의 거르기를 수행하는 것을 필요로 한다. 더 복잡도가 적게 드는 대안은 각 믿음 상태에 대해 각 단계마다 이산적 상태를 적분해 빼내어 단일 가우시안으로 표현하는 것이다. 이는 가정된 밀도 거르기법의 응용이다. 이의 오프라인 버전은 기대값 보정이라 한다.

다른 휴리스틱한 접근법은 상호 다중 모델(IMM)로, 사전분포를 모멘트 짝맞춤을 통해 단일 가우시안으로 근사한 뒤 K개의 다른 칼만 거르기를 q_{t}의 값별로 업데이트하는 것이다.

18.6.2. Application: data association and multi-target tracking

K개의 물체를 추적할 때 K^{\prime}개의 감지 이벤트가 발생한다고 하자. 이 사이의 대응을 알아내는 문제를 데이터 연관 문제라 한다. 거리에 기반해 데이터를 연관하는 방법을 최근접 근방 데이터 연관 휴리스틱이라 한다. 이는 K \times K^{\prime} 가중치 행렬을 낳는데, 이를 N = \max(K, K^{\prime})에 대해 N \times N 행렬로 확장한 뒤 O(N^{3})헝가리안 알고리즘을 이용해 최대 가중치 이분 짝짓기를 찾을 수 있다. 이를 통해 칼만 거르기 업데이트를 수행할 수 있다. 이의 응용으로 다중 대상 추적법이 쓰인다.

18.6.3. Application: fault diagnosis

근사적 추론을 통해 모델의 결함 판별을 수행할 수 있다.

18.6.4. Application: econometric forecasting

변환 선형 가우시안 상태 공간 모델은 계량 경제 예측국면 전환 모형에 쓰이곤 한다.

요점 정리

  • 상태 공간 모델 : 은닉 상태가 연속적인 은닉 마르코프 모델.
  • 상태 공간 모델은 위치 추적, 시계열 예측 등에 쓰임.
  • 선형 가우시안 상태 공간 모델의 추론은 칼만 필터링 알고리즘으로 이루어짐.
  • 선형 가우시안 상태 공간 모델의 학습은 기대값 최대화법이나 부분공간법을 사용함. 은닉 상태를 전부 알 수 있다면 선형 회귀식으로 풀림.
  • 비선형 비가우시안 상태 공간 모델은 근사적으로 추론할 수밖에 없음. 확장된 칼만 필터, 무취 칼만 필터 등을 사용함.
  • 이산적/연속적 은닉 변수가 혼합된 상태 공간 모델은 가우시안 합 필터법을 이용함.

답글 남기기

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

WordPress.com 로고

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

Google photo

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

Twitter 사진

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

Facebook 사진

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

%s에 연결하는 중