일단, 칼만 필터가 무엇인지 정확하게 알아보기 전에, 이것이 쓰인 대표적인 사건을 알아보자.
바로 아폴로 미션이다.
아무것도 없는 우주 공간에서 현재 자신의 위치를 아는 것은 굉장히 중요하다. 이를 기반으로 목적지까지 도달하는 방법이 달라지기 때문이다. 심지어 우주선의 위치는 계속하여 바뀌고, 측정기의 성능이 완벽하지 않으므로 이는 더욱 어려워진다. 외력이 거의 없으므로 예측하기 쉽겠다고 생각하는 사람이 있을지도 모르니, 한 극단적인 예시를 들어보자.
인류 역사상 가장 성공적인 실패라고 불리는 ‘아폴로 13 미션’은 달 탐사를 가던 도중, 산소 탱크가 폭파되어 탐사를 포기하고 지구로 복귀한 사건이다. 우주선의 동력과 내부 산소 등이 모두 부족한 극한 환경에서 대원들은 기적적으로 모두 무사 복귀하였다.
이 아폴로 13 미션에서는, 지구로 복귀하는 과정에서 그 어느때보다 현재 자신의 위치를 알아내는 기법인 localization이 중요했다. 만일 그 위치를 조금이라도 어긋나게 계산한다면 중력권을 벗어나서 영영 우주미아가 되거나 중력의 영향을 너무 세게 받아서 지면과 큰 충돌이 일어나 죽게 될지도 모르는 상황이었다.
심지어는 그 당시에 산소가 새고 있었기 때문에, 우주선은 수시로 랜덤하게 흔들렸고 이는 더욱 localization을 까다롭게 했다.
더욱이, 그 당시의 컴퓨터는 우리의 상상보다 훨씬 느렸고, 전력도 심각하게 부족하여, 정밀한 값을 얻기 위해 큰 비용의 계산을 시킬 수 없는 상황이었다.
그때 이들을 구원한 것은 바로 칼만 필터(Kalman Filter)였다. 이는 지금까지의 측정값을 기반으로 현재 위치를 추정하는 방법으로, 많이 측정할 수록 재귀적으로 정확도가 높아지는 경향이 있는 알고리즘이다. 또한, 기존의 알고리즘에 비해 매우 빠르고 값싸게 수행될 수 있어서 비단 아폴로 13 미션 뿐만 아닌, 다른 아폴로 미션에서도 유용하게 쓰였던 알고리즘이다.
결국, 조금 과장하자면 그 당시 칼만 필터가 없었다면 아마도 아폴로 13 미션은 인류 역사에 있어서 굉장히 끔찍한 참사로 남게 되었을지도 모른다. 그래서, NASA 측에서도 칼만 필터의 최초의 활용 사례로 아폴로 미션을 꼽고 있으며, 칼만 필터의 고안에 투자를 한 투자사도 이를 대대적으로 홍보하고 있다.
localization 예시를 계속 생각해보자. t t t x t x_t x t y t y_t y t
이때, 실제 위치인 x t x_t x t 1↗ y t y_t y t
여기서 다음과 같은 종속관계를 가정하자.
{ x t + 1 = α x t + β + ε x y t = λ x t + ε y \begin{cases}
x_{t + 1} &= \alpha x_{t} + \beta + \varepsilon_{x} \\
y_{t} &= \lambda x_t + \varepsilon_y
\end{cases} { x t + 1 y t = α x t + β + ε x = λ x t + ε y 단, 여기서 α , β ∈ R \alpha, \beta \in \R α , β ∈ R ε x , ε y \varepsilon_x, \varepsilon_y ε x , ε y 0 0 0 x 0 x_0 x 0
당연히 다음 위치는 이전 위치와 밀접한 관계가 있고, 측정값은 실제 위치와 밀접한 관계가 있다고 생각하면 이 가정은 자연스럽다.
이러면, 다음과 같은 graphical model을 생각해볼 수 있다.
x 0 → x 1 → ⋯ → x t → ⋯ → x T ↓ ↓ ↓ ↓ y 0 → y 1 → ⋯ → y t → ⋯ → y T \begin{CD}
x_0 @>>> x_1 @>>> \cdots @>>> x_t @>>> \cdots @>>> x_T \\
@VVV @VVV @. @VVV @. @VVV \\
y_0 @>>> y_1 @>>> \cdots @>>> y_t @>>> \cdots @>>> y_T
\end{CD} x 0 ↓ ⏐ y 0 x 1 ↓ ⏐ y 1 ⋯ ⋯ x t ↓ ⏐ y t ⋯ ⋯ x T ↓ ⏐ y T 즉, joint probability p ( x 0 , ⋯ , x T , y 0 , ⋯ , y T ) p(x_0, \cdots, x_T, y_0 , \cdots, y_T) p ( x 0 , ⋯ , x T , y 0 , ⋯ , y T )
p ( x 0 ) p ( y 0 ∣ x 0 ) p ( x 1 ∣ x 0 ) p ( y 1 ∣ x 1 ) ⋯ p ( x T ∣ x T − 1 ) p ( y T ∣ x T ) p(x_0)p(y_0 | x_0) p(x_1 | x_0)p(y_1 | x_1) \cdots p(x_{T} | x_{T-1}) p(y_T | x_T) p ( x 0 ) p ( y 0 ∣ x 0 ) p ( x 1 ∣ x 0 ) p ( y 1 ∣ x 1 ) ⋯ p ( x T ∣ x T − 1 ) p ( y T ∣ x T ) 여기서 우리의 목적은 기존의 측정값인 y 0 , ⋯ , y T − 1 y_0, \cdots, y_{T-1} y 0 , ⋯ , y T − 1 y T y_T y T x T x_T x T
p ( x T ∣ y 0 , ⋯ , y T ) p(x_T | y_0, \cdots, y_T) p ( x T ∣ y 0 , ⋯ , y T ) 이를 위해서는 x 0 , ⋯ , x T − 1 x_0 , \cdots, x_{T-1} x 0 , ⋯ , x T − 1 T T T y 0 , ⋯ , y T y_0, \cdots, y_T y 0 , ⋯ , y T
gaussian을 marginalize 해도 gaussian!
gaussian을 conditioning 해도 gaussian!
즉, 결과물은 gaussian이므로 그것의 평균과 분산만 알면 되는 것이다.
또한, 다음을 알고 있으므로
μ 1 ∣ 2 = μ 1 + Σ 12 Σ 2 − 1 ( x 2 − μ 2 ) \bm\mu_{1|2} = \bm\mu_1 + \Sigma_{12} \Sigma_{2}^{-1} ({\bf x}_2 - \bm\mu_2) μ 1∣2 = μ 1 + Σ 12 Σ 2 − 1 ( x 2 − μ 2 ) x T ∣ y 0 , ⋯ , y T x_{T} | y_0, \cdots, y_T x T ∣ y 0 , ⋯ , y T
μ x T ∣ y 0 , ⋯ , y T = μ x T + Σ x T , y 0 , ⋯ , y T Σ y 0 , ⋯ , y T − 1 ( ( y 0 , ⋯ , y T ) ⊤ − μ y 0 , ⋯ , y T ) \mu_{x_T | y_0, \cdots, y_T} = \mu_{x_T} + \Sigma_{x_T, y_0, \cdots, y_T} \Sigma _{y_0, \cdots, y_T}^{-1} ((y_0, \cdots, y_T)^\top - \bm\mu_{y_0, \cdots, y_T}) μ x T ∣ y 0 , ⋯ , y T = μ x T + Σ x T , y 0 , ⋯ , y T Σ y 0 , ⋯ , y T − 1 (( y 0 , ⋯ , y T ) ⊤ − μ y 0 , ⋯ , y T ) 와! 근데…
물론 closed form이긴 한데… 조금 생각해보면 이 부분이 마음에 걸린다.
Σ y 0 , ⋯ , y T − 1 \Sigma _{y_0, \cdots, y_T}^{-1} Σ y 0 , ⋯ , y T − 1 이 covariance matrix는 R ( T + 1 ) × ( T + 1 ) \R^{(T+1) \times (T+1)} R ( T + 1 ) × ( T + 1 )
Ω ( ( T + 1 ) 3 ) \Omega((T+1)^3) Ω (( T + 1 ) 3 ) 세제곱의 order를 가지는 알고리즘은 흔히 ‘좋지 않다고 ’ 여겨진다.
심지어, 6~70년대의 컴퓨터로 이정도의 계산을 하려고 하면 정말 오래 걸릴 뿐더러, 아폴로 13 미션 당시의 급박한 상황에서는 섣불리 이러한 계산을 할 엄두가 나지 않았을 것이다.
결국 inverse가 가장 큰 overhead인데, 이를 최소화 하는 방법이 있을까? 그렇다. 없었다면, 아폴로 미션이 존재하지 못했을 것이다.
Option II - Recursive update 한꺼번에 closed form을 구하려고 하지 말고, 앞서 가정한 종속관계에 따라 이를 재귀적으로 update해보자.
현재 시각 T T T t t t
앞으로, 혼동이 없을 때에는 편하게 다음과 같은 표기를 사용한다.
t + 1 ∣ t + 1 : = x t + 1 ∣ y 0 , ⋯ , y t + 1 {t+1 ~|~ t+1} := {x_{t+1} | y_0, \cdots, y_{t+1}} t + 1 ∣ t + 1 := x t + 1 ∣ y 0 , ⋯ , y t + 1 t + 1 , t + 1 ∣ t : = x t + 1 , y t + 1 ∣ y 0 , ⋯ , y t t+1 , t+1 ~|~ t := x_{t+1}, y_{t+1} | y_{0} ,\cdots, y_t t + 1 , t + 1 ∣ t := x t + 1 , y t + 1 ∣ y 0 , ⋯ , y t Observation update: μ t + 1 ∣ t + 1 \mu _{t+1 | t+1} μ t + 1∣ t + 1 먼저, 새로운 측정값 y t + 1 y_{t + 1} y t + 1
conditional mean은 다음과 같다.
μ t + 1 ∣ t + 1 = μ t + 1 ∣ t + σ t + 1 , t + 1 ∣ t σ y t + 1 ∣ t − 2 ( y t + 1 − E [ y t + 1 ∣ t ] ) \mu_{t+1 | t+1} = \mu_{t+1 | t} + \sigma_{t+1, t+1 | t} ~\sigma^{-2}_{y_{t+1} | t} ~(y_{t+1} - \mathbb{E}[y_{t+1} | t]) μ t + 1∣ t + 1 = μ t + 1∣ t + σ t + 1 , t + 1∣ t σ y t + 1 ∣ t − 2 ( y t + 1 − E [ y t + 1 ∣ t ]) 이를 정리하기에 앞서, 이전에 overhead가 되었던 부분인 σ y t + 1 ∣ t \sigma_{y_{t+1} | t} σ y t + 1 ∣ t O ( 1 ) O(1) O ( 1 )
이것이 t = 1 , ⋯ , T t = 1 , \cdots, T t = 1 , ⋯ , T O ( T ) O(T) O ( T ) μ T ∣ T \mu_{T | T} μ T ∣ T
엄청난 발전이다. 하지만 아직 놀라긴 이르다. 위 식을 좀 더 정리하기 위해서 σ t + 1 , t + 1 ∣ t \sigma_{t+1, t+1 | t} σ t + 1 , t + 1∣ t σ y t + 1 ∣ t \sigma_{y_{t+1} | t} σ y t + 1 ∣ t
σ y t + 1 ∣ t 2 = E [ ( λ x t + 1 + ε y ) 2 ∣ t ] − E [ λ x t + 1 + ε y ∣ t ] 2 = λ 2 E [ x t + 1 2 ∣ t ] + λ E [ x t + 1 ε y ∣ t ] + E [ ε y 2 ∣ t ] − λ 2 E [ x t + 1 ∣ t ] 2 \begin{align*}
\sigma^2_{y_{t+1} | t} &= \mathbb{E}[(\lambda x_{t+1} + \varepsilon_y)^2 | t] - \mathbb{E}[\lambda x_{t+1} + \varepsilon_y | t]^2 \\
&= \lambda ^2 \mathbb{E}[x_{t+1}^2 | t] + \lambda \mathbb{E}[x_{t+1}\varepsilon_y | t ] + \mathbb{E}[\varepsilon_y^2 | t] - \lambda^2 \mathbb{E}[x_{t+1} | t]^2
\end{align*} σ y t + 1 ∣ t 2 = E [( λ x t + 1 + ε y ) 2 ∣ t ] − E [ λ x t + 1 + ε y ∣ t ] 2 = λ 2 E [ x t + 1 2 ∣ t ] + λ E [ x t + 1 ε y ∣ t ] + E [ ε y 2 ∣ t ] − λ 2 E [ x t + 1 ∣ t ] 2 여기서 ε y \varepsilon_y ε y y 0 , ⋯ , y t y_0 ,\cdots, y_t y 0 , ⋯ , y t E [ ε y ∣ t ] = E [ ε y ] = 0 \mathbb{E}[\varepsilon _y | t] = \mathbb{E}[\varepsilon_y] = 0 E [ ε y ∣ t ] = E [ ε y ] = 0
x t + 1 x_{t+1} x t + 1 E [ x t + 1 ε y ∣ t ] = E [ x t + 1 ∣ t ] E [ ε y ∣ t ] = 0 \mathbb{E}[x_{t+1} \varepsilon_y | t] = \mathbb{E}[x_{t+1} | t]\mathbb{E}[ \varepsilon_y | t] = 0 E [ x t + 1 ε y ∣ t ] = E [ x t + 1 ∣ t ] E [ ε y ∣ t ] = 0
또한, E [ x t + 1 2 ∣ t ] − E [ x t + 1 ∣ t ] 2 = σ t + 1 ∣ t 2 \mathbb{E}[x_{t+1} ^2 | t] - \mathbb{E}[x_{t+1} | t]^2 = \sigma_{t+1 | t}^2 E [ x t + 1 2 ∣ t ] − E [ x t + 1 ∣ t ] 2 = σ t + 1∣ t 2
σ y t + 1 ∣ t 2 = λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] \sigma_{y_{t+1} | t}^2 = \lambda^2 \sigma_{t+1 | t}^2 + \mathbb{E}[\varepsilon^2_y] σ y t + 1 ∣ t 2 = λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] σ t + 1 , t + 1 ∣ t \sigma_{t+1, t+1 | t} σ t + 1 , t + 1∣ t
σ t + 1 , t + 1 ∣ t = E [ x t + 1 ( λ x t + 1 + ε y ) ∣ t ] − E [ x t + 1 ∣ t ] E [ λ x t + 1 + ε y ∣ t ] = λ E [ x t + 1 2 ∣ t ] − λ E [ x t + 1 ∣ t ] 2 = λ σ t + 1 ∣ t 2 \begin{align*}
\sigma_{t+1, t+1 | t} &= \mathbb{E}[x_{t+1}(\lambda x_{t+1} + \varepsilon_y) | t] - \mathbb{E}[x_{t+1} | t]~\mathbb{E}[\lambda x_{t+1} + \varepsilon_y | t] \\
&= \lambda \mathbb{E}[x_{t+1} ^2 | t ] - \lambda\mathbb{E}[x_{t+1} | t]^2 \\
&= \lambda \sigma_{t+1 | t}^2
\end{align*} σ t + 1 , t + 1∣ t = E [ x t + 1 ( λ x t + 1 + ε y ) ∣ t ] − E [ x t + 1 ∣ t ] E [ λ x t + 1 + ε y ∣ t ] = λ E [ x t + 1 2 ∣ t ] − λ E [ x t + 1 ∣ t ] 2 = λ σ t + 1∣ t 2 따라서, mean은 다음과 같다.
μ t + 1 ∣ t + 1 = μ t + 1 ∣ t + λ σ t + 1 ∣ t 2 λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] ( y t + 1 − λ μ t + 1 ∣ t ) = E [ ε y 2 ] λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] μ t + 1 ∣ t + λ 2 σ t + 1 ∣ t 2 λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] y t + 1 λ \begin{align*}
\mu_{t+1 | t+1} &= \mu_{t+1 | t} + \frac{\lambda \sigma^2_{t+1 | t}}{\lambda^2 \sigma^2_{t+1 | t} + \mathbb{E}[\varepsilon^2_y ]} (y_{t+1} - \lambda \mu_{t+1 | t}) \\
&=\textcolor{crimson}{\frac{\mathbb{E}[\varepsilon_y^2]}{{\lambda^2 \sigma^2_{t+1 | t} + \mathbb{E}[\varepsilon^2_y ]}}} \mu_{t+1 | t} + \textcolor{cornflowerblue}{\frac{{\lambda^2 \sigma_{t+1 | t}^2 }}{{\lambda^2 \sigma^2_{t+1 | t} + \mathbb{E}[\varepsilon^2_y ]}} }\frac{y_{t+1}}{\lambda}
\end{align*} μ t + 1∣ t + 1 = μ t + 1∣ t + λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] λ σ t + 1∣ t 2 ( y t + 1 − λ μ t + 1∣ t ) = λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] E [ ε y 2 ] μ t + 1∣ t + λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] λ 2 σ t + 1∣ t 2 λ y t + 1 Observation update: σ t + 1 ∣ t + 1 2 \sigma^2_{t+1 | t+1} σ t + 1∣ t + 1 2 이제는 분산을 구해보자. conditional covariance는 다음과 같이 구할 수 있다.
σ t + 1 ∣ t + 1 2 = σ t + 1 ∣ t 2 − σ t + 1 , t + 1 ∣ t σ y t + 1 ∣ t − 2 σ t + 1 , t + 1 ∣ t \sigma^2_{t+1 | t+1} = \sigma^2_{t+1 | t} - \sigma_{t+1, t+1 | t} \sigma^{-2}_{y_{t+1} | t} \sigma_{t+1, t+1 | t} σ t + 1∣ t + 1 2 = σ t + 1∣ t 2 − σ t + 1 , t + 1∣ t σ y t + 1 ∣ t − 2 σ t + 1 , t + 1∣ t 여기서 σ t + 1 , t + 1 ∣ t = λ σ t + 1 ∣ t 2 \sigma_{t+1, t+1 | t} = \lambda \sigma_{t+1 | t}^2 σ t + 1 , t + 1∣ t = λ σ t + 1∣ t 2 σ y t + 1 ∣ t 2 = λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] \sigma^{2}_{y_{t+1} | t} = \lambda^2 \sigma_{t+1 | t}^2 + \mathbb{E}[\varepsilon^2_y] σ y t + 1 ∣ t 2 = λ 2 σ t + 1∣ t 2 + E [ ε y 2 ]
σ t + 1 ∣ t + 1 2 = σ t + 1 ∣ t 2 − λ 2 σ t + 1 ∣ t 4 λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] = E [ ε y 2 ] λ 2 σ t + 1 ∣ t 2 + E [ ε y 2 ] σ t + 1 ∣ t 2 \begin{align*}
\sigma^{2}_{t+1 | t+1} &= \sigma^2_{t+1 | t} - \cfrac{\lambda^2 \sigma_{t+1 | t}^4}{\lambda^2 \sigma_{t+1 | t}^2 + \mathbb{E}[\varepsilon^2_y]} \\
&= \textcolor{crimson}{\cfrac{ \mathbb{E}[\varepsilon^2_y]}{\lambda^2 \sigma_{t+1 | t}^2 + \mathbb{E}[\varepsilon^2_y]}} ~ \textcolor{cornflowerblue}{\sigma_{t+1 | t}^2}
\end{align*} σ t + 1∣ t + 1 2 = σ t + 1∣ t 2 − λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] λ 2 σ t + 1∣ t 4 = λ 2 σ t + 1∣ t 2 + E [ ε y 2 ] E [ ε y 2 ] σ t + 1∣ t 2 여기서 붉게 표시된 부분은 1 이하의 양수이고, 파랗게 표시된 부분은 이전의 분산이다.
따라서, 새로운 데이터 y t + 1 y_{t+1} y t + 1
만일 E [ ε y 2 ] \mathbb{E} [\varepsilon^2_y] E [ ε y 2 ]
Time update: μ t + 1 ∣ t \mu_{t+1 | t} μ t + 1∣ t 그럼, 시간이 흐름에 따라 추정하는 경우를 생각해보자. 즉, 이번에는 새로운 데이터 y t + 1 y_{t+1} y t + 1 t t t t + 1 t+1 t + 1
conditional mean은 다음과 같다.
μ t + 1 ∣ t = E [ α x t + β + ε x ∣ t ] = α μ t ∣ t + β \begin{align*}
\mu_{t+1 | t} &= \mathbb{E}[\alpha x_t + \beta + \varepsilon_x | t] \\
&= \alpha \mu_{t | t} + \beta
\end{align*} μ t + 1∣ t = E [ α x t + β + ε x ∣ t ] = α μ t ∣ t + β 즉, 이전 mean을 update한 것이다. 직관적이다.
Time update: σ t + 1 ∣ t 2 \sigma^2_{t+1 | t} σ t + 1∣ t 2 conditional covariance는 다음과 같다.
σ t + 1 ∣ t 2 = E [ x t + 1 2 ∣ t ] − E [ x t + 1 ∣ t ] 2 = E [ ( α x t + β + ε x ) 2 ∣ t ] − E [ α x t + β + ε x ∣ t ] 2 = α 2 σ t ∣ t 2 + E [ ε x 2 ] \begin{align*}
\sigma^2_{t+1|t} &= \mathbb{E}[x_{t+1}^2 | t ] - \mathbb{E}[x_{t+1}|t]^2 \\
&= \mathbb{E}[(\alpha x_t + \beta + \varepsilon_x)^2 | t] - \mathbb{E}[\alpha x_t + \beta + \varepsilon_x | t]^2 \\
&= \alpha^2 \sigma^{2}_{t | t} + \mathbb{E}[\varepsilon_x^2]
\end{align*} σ t + 1∣ t 2 = E [ x t + 1 2 ∣ t ] − E [ x t + 1 ∣ t ] 2 = E [( α x t + β + ε x ) 2 ∣ t ] − E [ α x t + β + ε x ∣ t ] 2 = α 2 σ t ∣ t 2 + E [ ε x 2 ] 이는 observation update와는 다르게, variance를 증가 시키는 term E [ ε x 2 ] \mathbb{E}[\varepsilon^2_x] E [ ε x 2 ]
직관적으로 생각하면, 시간이 지날수록 이전 추정값에 대한 신뢰도가 낮아진다고 생각할 수 있다.
결국 이 time update와 observation update를 반복하면 임의의 시점에서의 위치를 추정할 수 있고, 그 정확도는 점점 높아진다는 것을 확인할 수 있다. 또한, 각 update는 O ( 1 ) O(1) O ( 1 )
이렇게 아폴로 미션은 성공할 수 있었다. 하지만, 사실 위치를 아는 것 만으로는 부족하다. 정확히는 현재 위치와 속도를 알아야 미래의 운동을 예측할 수 있기 때문이다. 따라서, 현실적으로는 위치와 속도를 묶어서 하나의 벡터(혹은 행렬)로 생각하여 kalman filter를 적용한다.
이에 대해서는 다음에 알아보자!
