卡尔曼滤波与隐马尔可夫模型
Published:
本文关注在动态场景下或者在序列数据上的经典机器学习模型，包括Bayes滤波、Kalman滤波、隐Markov模型的前向-后向等求解算法。
主要参考了 张志华‘强化学习基础 中的相关内容。
Prediction With Hidden State
隐马尔可夫模型（Hidden Markov Model，HMM）是概率图模型的一种，关于概率图模型，感兴趣的读者可以移步至 概率图模型
模型假设具有Markov性质，也即未来仅仅依赖于当前而不依赖于过去。
\[\begin{align} p(s_{t+1} \vert s_{1:t},y_{1:t}) &= p(s_{t+1} \vert s_t) \\ p(y_{t+1} \vert s_{1:t+1},y_{1:t}) &= p(y_{t+1} \vert s_{t+1}) \end{align}\]其中，上述的$s_t$表示时刻$t$之下的状态，而该状态不能被直接观测，为隐变量，而$y_t$为$t$时刻之下的观测数据，为显变量。
对于隐马尔可夫模型的重要任务是利用显变量分布，推断出隐变量的分布，利用Bayes统计的角度理解就是计算一个后验分布。在一个分布的基础上，可以根据最小化某个损失函数的原则，对数据进行预测。
下面介绍最小化$L_1,L_2$范数误差和Dirichlet损失下利用分布产生的预测结果，
定理1 最小化$L_2$范数误差（均方误差）等价于求分布的期望，也即，
\[\begin{align} \min_a E \Vert X - a \Vert_2^2 = E[X] \end{align}\]证明可以使用求导，或者直接证明其最优性，对于 $\forall c$,
\[\begin{align} &E[(X-c)^T(X-c)]- E[(X-EX)^T(X-EX)] \\ =& E[-2X^T c+c^Tc-E[X]^T E[X] +2X^T E[X] ] \\ =&c^T c-2E[X]^T c+ E[X]^T E[X] \\ =& (E[X] - c)^T(E[X]-c) \\ \ge& 0 \end{align}\]本质上，上面结论成立的原因是因为$L_2$范数是一个Bregman散度（Bregman Divergence），其定义为，
\[\begin{align} L(x,y) &= \phi(x) - \phi(y) - \langle \nabla\phi(y),x-y \rangle \\ \text{With } L(x,y) &\ge 0 \text{ And } L(x,y) = 0 \text{ Iff } x=y \end{align}\]该定义和次梯度非常相近，可以证明对于Bregman散度，
\(\min_a E[L(X,a)] = E[X]\) 证明和对于$L_2$范数的情况类似，对于任意其他取值，$\forall c$
\[\begin{align} E[L(X,c)] -E[L(X,a)] &= E[\phi(a) - \phi(c) + \nabla \phi(a)^T(X-a) - \nabla \phi(c)^T(X-c) ] \\ &= E[\phi(a) - \phi(c) - \nabla\phi(c)^T (a-c)] \\ &= \phi(a) - \phi(c) - \nabla\phi(c)^T (a-c) \\ &=L(a,c) \\ &\ge 0 \end{align}\]例如对于经典的回归问题，采用的$L_2$损失是一个Bregman散度，只需验证，
\[\begin{align} L(x,y) &= \phi(x) - \phi(y) - \nabla \phi(y)^T(x-y) \\ &=x^2 - y^2 -2y(x-y), \text{ Let } \phi(x) = x^2 \\ &=(x-y)^2 \\ \end{align}\]而对于分类问题，通常采用基于KL散度的交叉熵损失，而KL散度本质上也是一个Bregman散度，
\(\begin{align} L(x,y) &= \phi(x) - \phi(y) - \nabla \phi(y)^T(x-y) \\ &=\sum_i x_i \log x_i - \sum_i y_i \log y_i - \sum_i (1 + \log y_i)(x_i -y_i) , \text{With } \phi(x)=\sum_i x_i \log x_i ,\sum_i x_i=1\\ &= \sum_i x_i \log x_i - \sum_i x_i \log y_i \\ &= \sum_i x_i \log \frac{x_i}{y_i} \\ &=KL(x \Vert y) \end{align}\) 关于KL散度与熵的内容，可以移步至 熵
定理2 最小化$L_1$范数误差（绝对值误差）等价于求分布的中位数，也即，
\[\min_a E \vert X -a \vert = \text{medium}(X)\]对于离散的情况，证明是简单的，下面仅证明连续随机变量的情况，下面仅对$c > a$的情况证明，$c < a$的情况是类似的，
\[\begin{align} &\forall c>a, \int_{-\infty}^a dF(x) = \int_a^{\infty} dF(x) = \frac{1}{2}, \\ & E \vert X- c \vert - E \vert X - a \vert \\ =& [\int_{-\infty}^a(c-x)dF(x) + \int_{a}^c(c-x)dF(x) + \int_c^{\infty} (x-c)dF(x)] -[\int_{-\infty}^a(a-x)dF(x) + \int_{a}^c(x-a)dF(x) + \int_a^{\infty} (x-a)dF(x)] \\ =& \int_{-\infty}^a(c-a)dF(x) + \int_{a}^c(a+c-2x)dF(x) + \int_c^{\infty} (a-c)dF(x) \\ =& c\int_{-\infty}^a dF(x) + c\int_a^c dF(x) -2\int_a^c xdF(x) - c\int_c^{\infty} dF(x) \\ \ge & c\int_{-\infty}^a dF(x) - c\int_a^c dF(x) - c\int_c^{\infty} dF(x) \\ =& c\int_{-\infty}^a dF(x) - c\int_a^{\infty} dF(x) \\ =& 0 \end{align}\]上述证明尽管看似复杂，但实际上只用到的初等的分段积分方法，较为简单。
定理3 最小化0-1损失（或称Dirichlet损失）等价于求分布的众数（Mode），该结论是显然的，无需证明即可直接得到，
\[\min_a \delta(X-a) = \text{mode}(X)\]从上面的几个结论我们可以发现，对于预测问题的核心实际上在于求解分布，对于观测到部分数据的情况下，就是要求解一个条件分布，或者称为后验分布，因此在后文中我们将关于于对后验分布的求解。
在类似的时移系统中，有如下的术语定义，
- 希望推测过去的后验分布，称为平滑（Smooth）
- 希望推测当下的后验分布，称为滤波（Filter）
- 希望推测未来的后验分布，称为预测（Predict)
Bayesian Filter
上述三个问题中，最常见的是滤波问题，例如在物体跟踪问题中，希望根据当前的观测，来推测当前系统的状态，也就是一个滤波任务。贝叶斯滤波是求解滤波问题的一个框架，基于递推的方法进行。
为了计算， $p(s_{t+1} \vert y_{1:{t+1}}) $ ,可以将其化归为用 $p(s_{t} \vert y_{1:{t}}) $ 表示从而可以递推进行, 使用Beyes公式即可，
\[\begin{align} p(s_{t+1} \vert y_{1:{t+1}}) &= \frac{p(y_{t+1} \vert s_{t+1}) p(s_{t+1} \vert y_{1:t}) }{ p (y_{t+1} \vert y_{1:t})} \\ &= \eta p(y_{t+1} \vert s_{t+1}) p(s_{t+1} \vert y_{1:t}) ,\text{With } \eta = \text{Const}\\ &= \eta p(y_{t+1} \vert s_{t+1}) \int p(s_{t+1} \vert s_t )p(s_{t} \vert y_{1:t}) ds_t \end{align}\]定义所求为 $p(s_{t} \vert y_{1:t}) $为递推中进行传播的信念，则上述公式说明了，
\[\text{belief}(s_{t+1}) = \eta p(y_{t+1} \vert s_{t+1}) \int p(s_{t+1} \vert s_t ) \text{belief}(s_t) ds_t\]上面实际上做了两步操作，第一步是求解得到新的$y_{t+1}$得到的$s_{t+1}$的后验分布，
\[p(s_{t+1} \vert y_{1:{t+1}}) = \eta p(y_{t+1} \vert s_{t+1}) p(s_{t+1} \vert y_{1:t})\]第二步是将$s_{t+1}$的先验分布递推地更新，也即求解$s_t$的分布已知的条件下$s_{t+1}$的条件分布，
\[p(s_{t+1} \vert y_{1:t}) =\int p(s_{t+1} \vert s_t )p(s_{t} \vert y_{1:t}) ds_t\]Kalman Filter
卡尔曼滤波（Kalman FIlter）是基于高斯过程建模的一个贝叶斯滤波，在目标定位追踪等问题中被广泛应用。
其假设状态与观测之间的隐Markov模型可以用高斯过程表示，
\[\begin{align} y_t &= A_t s_t + u_t, u_t \sim \mathcal{N}(0, R_t)\\ s_{t+1} &= B_t s_t + v_t , v_t \sim \mathcal{N}(0, Q_t)\\ \end{align}\]利用Beyes滤波求解，首先求解根据观测样本$y_t$ 得到的$s_t$的后验分布，由于先验分布为正态分布，根据正态分布的均值的共轭分布性质，我们知道其对应的后验分布也是正态分布，下面我们显式求解该后验分布，
假设我们已经得到了 $s_t \sim \mathcal{N}(\mu, \Sigma)$ ,其后验分布，
\(\begin{align} p(s_t \vert y_t ) &= \eta p(y_t \vert s_t) p(s_t) \\ &= K_1 \exp(-\frac{1}{2}(As-y)^TR^{-1} (As-y)) \exp(-\frac{1}{2}(s-\mu)^T \Sigma^{-1} (s-\mu)) \\ &= K_2 \exp(-\frac{1}{2}[s^T (A^T R^{-1} A+ \Sigma^{-1})s- 2 s^T (A^T R^{-1} y +\Sigma^{-1} \mu)]) \end{align}\) 因此可以得到后验分布为，
\[\begin{align} p(s_t \vert y_t) &\sim \mathcal{N}(\tilde \mu, \tilde \Sigma) \\ \text{With } \tilde \mu &= \tilde \Sigma(A^T R^{-1} y+ \Sigma^{-1} \mu) \\ \tilde \Sigma &= (A^T R^{-1} A+ \Sigma^{-1})^{-1} \end{align}\]在计算上，考虑到上式中进行的矩阵求逆操作复杂度较高，而观测变量$y_t$相比起隐状态$s_t$的维度通常小得多，考虑使用Sherman–Morrison–Woodbury公式进行简化运算，
\[(A^{-1} + UB^{-1}V^T)^{-1} = A -AU(B+V^TAU)^{-1} V^T A\]为了减少重复运算定义辅助变量$K$，并且代入上述公式，
\[\begin{align} \tilde \Sigma &= \Sigma - \Sigma A^T (R+A \Sigma A^T)^{-1} A \Sigma \\ &= \Sigma - K A \Sigma \\ \text{With } K &= \Sigma A^T (R+A \Sigma A^T)^{-1} \\ \end{align}\]并且利用关于$K$的等式可以得到,
\[\begin{align} \tilde \mu &= \tilde \Sigma(A^T R^{-1} y+ \Sigma^{-1} \mu) \\ &=(I -KA) \Sigma (A^T R^{-1} y+ \Sigma^{-1} \mu) \\ &= \mu + K(y - A \mu) \\ \text{With } &= (I-KA)\Sigma A^T R^{-1} = K \end{align}\]综上我们得到了求解上述后验分布的动态更新公式，
\[\begin{align} K &= \Sigma A^T (R+A \Sigma A^T)^{-1} \\ \tilde \Sigma &= \Sigma - K A \Sigma \\ \tilde \mu &= \mu + K(y - A \mu) \\ \end{align}\]可见后验分布的均值和方差是在原来的均值和方差的基础上根据$K$矩阵进行调整，$K$也称为Kalman信息增益。
上述方法可也用在递归最小二乘估计中（RLS，Recursive Least Square），也即根据动态获得的数据进行回归，将回归方程写为，
\[y_t = X_t \beta + \epsilon, \epsilon \sim \mathcal{N}(0, \sigma^2)\]即可使用上述算法进行求解。
除了动态更新后验分布，Beyes滤波需要对$s_t$的条件分布进行更新，利用转移等式，
\[s_{t+1} = B_t s_t + v_t , v_t \sim \mathcal{N}(0, Q_t)\]假设我们已经得到了 $s_t \sim \mathcal{N}(\mu, \Sigma)$ ,下一个状态$s_{t+1}$的条件分布可以根据正态分布的线性性质直接得到，
\[s_{t+1} \sim \mathcal{N}(B_t \mu, B_t \Sigma B_t^T + Q_t)\]将所有式子整合起来，总结Kalman滤波的公式，
\[\begin{align} \tilde \mu_{t+1} &= B_t \mu_t, \tilde \Sigma_{t+1} = B_t \Sigma_t B_t^T +Q_t \\ \mu_{t+1} &= \tilde \mu_{t+1} + K(y - A \tilde \mu_{t+1}) ,\Sigma_{t+1} = \tilde \Sigma - K A \tilde \Sigma \\ \end{align}\]Kalman滤波算法和 概率图模型 中的信念传播算法和 EM算法等都有异曲同工之妙。
Hidden Markov Model
本节我们关注于带隐变量的Markov链，关注于在观测数据下每个状态的条件概率，也即，
\[p(s_t \vert y_{1:n}), t <n\]该问题与滤波问题的区别在于滤波问题是动态地更新 $p(s_t \vert y_{1:t})$，而上述问题要求在给定观测序列的前提下计算状态序列的概率，同样的可以基于动态规划的递推方法进行。
定义如下的$\alpha,\beta$, 可以利用局部Markov性质将所求概率拆分为两部分，该方法称为前向-后向算法（Forward-Backward Algorithm）
\[\begin{align} p(s_t \vert y_{1:n}) & = \frac{p(s_t,y_{1:n})}{p(y_{1:n})} \\ &= \frac{p(s_t,y_{1:t}) p(y_{t+1:n} \vert s_t) }{p(y_{1:n})} \\ &= \frac{\alpha(s_t) \beta(s_t)}{\sum_{s_t}\alpha(s_t) \beta(s_t)} \\ \text{With } \alpha(s_t) &= p(s_t,y_{1:t}) , \beta(s_t) = p(y_{t+1:n} \vert s_t) \end{align}\]考虑基于递推方法求解$\alpha,\beta$, 对于$\alpha$可以从前向进行递推，
\(\begin{align} \alpha(s_{t+1}) &= p(s_{t+1},y_{1:t+1}) \\ &= \sum_{s_t} p(s_t,s_{t+1},y_{1:t+1}) \\ &= \sum_{s_t} p(s_t,y_{1:t}) p(s_{t+1} \vert s_t) p(y_{t+1} \vert s_{t+1}) \\ &= \sum_{s_t} \alpha(s_t) p(s_{t+1} \vert s_t) p(y_{t+1} \vert s_{t+1}) \end{align}\) 而对于$\beta$可以从后向进行递推，
\[\begin{align} \beta(s_t) &= p(y_{t+1:n} \vert s_t) \\ &= \sum_{s_{t+1}} p(s_{t+1} ,y_{t+1:n} \vert s_t) \\ &= \sum_{s_{t+1}} p(y_{t+2:n} \vert s_{t+1}) p(y_{t+1} \vert s_{t+1}) p(s_{t+1} \vert s_t) \\ &= \sum_{s_{t+1}} \beta(s_{t+1}) p(y_{t+1} \vert s_{t+1}) p(s_{t+1} \vert s_t) \end{align}\]另外一种改进的方法是在前向算法得到$\alpha$的前提下，在后向算法中直接得到结果，记结果为$\gamma$,满足定义，
\(\gamma(s_t) = p(s_t \vert y_{1:n})\) 下面求解其递推方法，
\[\begin{align} \gamma(s_t) &= p(s_t \vert y_{1:n}) \\ &=\sum_{s_{t+1}} p(s_t,s_{t+1} \vert y_{1:n}) \\ &= \sum_{s_{t+1}} p(s_t \vert s_{t+1},y_{1:n})p(s_{t+1} \vert y_{1:n}) \\ &= \sum_{s_{t+1}} p(s_t \vert s_{t+1},y_{1:t})p(s_{t+1} \vert y_{1:n}) ,\text{By Local Markov Property}\\ &= \sum_{s_{t+1}} \frac{p(s_t ,s_{t+1},y_{1:t})}{p(s_{t+1},y_{1:t})}\gamma(s_{t+1}) \\ &= \sum_{s_{t+1}} \frac{p(s_t,y_{1:t})p(s_{t+1}\vert s_t)}{p(s_{t+1},y_{1:t})}\gamma(s_{t+1}) \\ &= \sum_{s_{t+1}} \frac{\alpha(s_t)p(s_{t+1}\vert s_t)}{\sum_{s_t} \alpha(s_t)p(s_{t+1}\vert s_t)}\gamma(s_{t+1}) ,\text{By Def of } \alpha(s_t) \end{align}\]