線性動態系統與卡爾曼濾波

線性動態系統與卡爾曼濾波

一、前言

在之前的文章中我們分別介紹了混合模型和隱馬爾科夫模型,詳見《混合模型與EM演算法》和《隱馬爾科夫模型》兩篇文章。隱馬爾科夫模型是一種時間序列模型,混合模型是一種與時間序列無關的模型。事實上,單看隱馬爾科夫模型的發射概率模型部分,它就是一個混合模型。通過轉移概率模型將一系列的混合模型的隱藏變數連成一個一階馬爾科夫鏈,就構成了隱馬爾科夫模型,因此用隱馬爾科夫模型可以表示比混合模型更豐富的一類模型。然而,根據定義,隱馬爾科夫模型的隱藏變數必須是離散的,而實際中遇到的許多問題需要對連續變數進行建模。因此,我們需要對隱馬爾科夫模型進行擴展,使其能夠對隱藏變數為連續變數的問題建模。在新開一篇文章介紹連續隱變數模型(隱藏變數為連續變數的混合模型)之前,我們先介紹線性動態系統(Linear Dynamic System,LDS),從中可以看到,卡爾曼濾波是LDS的一個子問題。

不論是離散隱變數模型還是連續隱變數模型,它們都是混合模型的一種。將它們的隱藏變數通過轉移概率模型連成一階馬爾科夫鏈,分別得到HMM和LDS。因此,從HMM到LDS,僅僅是隱藏變數從離散變為連續,參數學習和推理演算法的框架是一模一樣的。

由於HMM和LDS的相似性,本文的所有記號均沿用《隱馬爾科夫模型》一文,不再重複說明。本文的大部分推導過程也與《隱馬爾科夫模型》一文一模一樣,只不過因為隱藏變數從離散變為連續了,推導公式中的所有求和要換成求積分。因此如果熟悉了《隱馬爾科夫模型》一文的內容,本文也就不難理解了。

在LDS中,要求發射概率模型和轉移概率模型都是線性高斯模型,因此

p(m{z_n} | m{z_{n-1}}) = N(m{z_n} | Am{z_{n-1},Gamma}) (1)

p(m{x_n} | m{z_n}) = N(m{x_n} | Cm{z_n}, Sigma) (2)

其中, AC 分別是表示線性關係的變換矩陣, GammaSigma 分別是兩個高斯分布的方差。

同理,隱藏變數的初始分布也是一個高斯分布,

p(m{z_1}) = N(m{z_1} | m{mu}_0, P_0) (3)

二、LDS中的推理(Inference)

在這一節中,我們先假設LDS的所有參數已知,介紹其推理過程。關於參數學習的內容在下一節介紹。

推理問題是要求解如下後驗分布

hat{alpha}(m{z}_n) = p(m{z}_n | m{x}_1, m{x}_2, cdots, m{x}_n) (4)

由於在LDS中發射概率模型和轉移概率模型都是線性高斯模型,因此 hat{alpha}(m{z}_n) 也是一個高斯分布,記作

hat{alpha}(m{z}_n) = N(m{z}_n | m{mu_n, V_n}) (5)

因此,我們的目標是求解此高斯分布的均值 m{mu_n} 和方差 V_n

在《隱馬爾科夫模型》一文中,我們推導了 hat{alpha}(m{z}_n) 的遞推式

c_n hat{alpha}(m{z}_n) = p(m{x}_n | m{z}_n)sum_{m{z}_{n-1}}{hat{alpha}(m{z}_{n-1}) p(m{z_n} | m{z}_{n-1})} (53)

只需要將式中的求和符號換成積分符號,即可得到LDS中的遞推式

c_n hat{alpha}(m{z}_n) = p(m{x}_n | m{z}_n)int{hat{alpha}(m{z}_{n-1}) p(m{z_n} | m{z}_{n-1})}dm{z}_{n-1} (6)

將(1)、(2)、(5)式代入(6)式得到高斯分布表示

c_n N(m{z}_n | m{mu_n, V_n}) = N(m{x_n} | Cm{z_n}, Sigma) int{ N(m{z}_{n-1} | m{mu_{n - 1}, V_{n-1}}) N(m{z_n} | Am{z_{n-1},Gamma})}dm{z}_{n-1} (7)

利用《機器學習中的高斯分布》一文中關於邊緣和條件高斯分布的結論,上式被積式是兩個高斯分布的卷積,因此結果也是一個高斯分布,且卷積結果為

int{ N(m{z}_{n-1} | m{mu_{n - 1}, V_{n-1}}) N(m{z_n} | Am{z_{n-1},Gamma})}dm{z}_{n-1} = N(m{z}_n | Am{mu}_{n-1}, P_{n-1}) (8)

其中,

P_{n-1} = AV_{n-1}A^T + Gamma (9)

這樣,(7)式的左右兩邊分別是兩個高斯分布的乘積,形如 p(x|y) p(y) = p(y|x) p(x),其中 p(x)p(y)p(x|y)p(y | x) 都是高斯分布,且 p(x)p(y | x) 已知,要求 p(y)p(x|y)。 再次利用《機器學習中的高斯分布》一文中的結論可以得到

m{mu}_n = A m{mu}_{n-1} + K_n (m{x}_n - CAm{mu}_{n-1}) (10)

V_n = (I - K_n C) P_{n-1} (11)

c_n = N(m{x}_n | CAm{mu}_{n-1}, CP_{n-1}C^T + Sigma) (12)

其中,

K_n = P_{n-1}C^T(CP_{n-1}C^T + Sigma)^{-1} (13)

叫做卡爾曼增益矩陣(Kalman Gain Matrix)。

我相信,當你利用《機器學習中的高斯分布》一文中的相關結論時,可以輕易寫出要求的高斯分布的均值和方差,但你寫出的結果很可能不是(10)~(13)式所示的形式。的確,要變換成(10)~(13)式的形式並不是那麼容易,其中會用到以下矩陣求逆公式

(A + BD^{-1}C) ^{-1} = A^{-1} - A^{-1}B(D + CA^{-1}B)^{-1}CA^{-1} (14)

我們不要求能從你一開始寫出的結果正向推導出形如(10)~(13)式的結論,只需要能夠根據(10)~(13)式的結論反向驗證其正確性即可。

遞推式有了,接下來只需要計算初始條件即可。因為

c_1 hat{alpha}(m{z}_1) = p(m{z}_1) p (m{x}_1 | m{z}_1) (15)

利用同樣的技巧,可以得到

m{mu}_1 = m{mu}_0 + K_1 (m{x}_1 - Cm{mu}_0) (16)

V_1 = (I - K_1C) P_0 (17)

c_1 = N(m{x}_1 | Cm{mu}_0, CP_0C^T + Sigma) (18)

其中,

K_1 = P_0 C^T (C P_0 C^T + Sigma)^{-1} (19)

至此,推理部分就推導完了。其中,(10)~(13)式的遞推式,加上(16)~(19)式的初始條件,即是卡爾曼濾波過程。從後文我們將看到,卡爾曼濾波只是LDS理論的一個子問題,除此之外,LDS還有其他理論問題需要解決。

我們可以回過頭來仔細觀察一下(10)~(13)式的結論,公式是冷冰冰的,因此需要用自然語言的方式解釋一下卡爾曼濾波到底是如何做濾波的。假設我們已經知道到第 n-1 步時隱藏變數的後驗分布的均值 m{mu}_{n - 1} 和方差 V_{n-1} ,怎麼求第 n 步的隱藏變數的後驗分布的均值 m{mu}_{n} 和方差 V_{n} 呢?首先,因為此時還未觀測到第 n 步的觀測變數 m{x}_n ,可以先用轉移概率關係預測 m{z}_n 的均值大約為  A m{mu}_{n-1} ;進一步,利用發射概率關係,可以預測觀測值大約為 CAm{mu}_{n-1}。當觀測到第 n 步的觀測變數 m{x}_n 後,我們發現觀測值和預測值有誤差,因此,我們可以利用此誤差對剛剛的預測進行修正,修正係數正是卡爾曼增益。因此,卡爾曼濾波的過程可以看做不斷對隱藏變數進行預測、並利用觀測值對預測值進行修正的過程。

三、LDS的參數學習(Learning)

同HMM一樣,LDS的參數學習也要用EM演算法,因此,也需要計算兩個關鍵的量,分別是gamma(m{z}_n)xi (m{z}_{n - 1}, m{z}_n)。與HMM不同的是,我們不需要通過先求 hat{eta}(m{z}_n) 間接求gamma(m{z}_n),而是直接求 gamma(m{z}_n) 的遞推式。因為根據定義, hat{eta}(m{z}_n) 是兩個高斯分布的比,它本身並不是一個高斯分布,gamma(m{z}_n) 才是一個高分布,因此,直接求 gamma(m{z}_n) 更容易。將 gamma(m{z}_n) 的高斯分布記為

gamma(m{z}_n) = N(m{z}_n | hat{m{mu}}_n, hat{V}_n)

根據《隱馬爾科夫模型》一文中的結論,

gamma(m{z}_n) = hat{alpha}(m{z}_n)hat{eta}(m{z}_n) (57)

hat{eta}(m{z}_n) 的遞推式為

c_{n+1}hat{eta}(m{z}_n) = sum_{m{z}_{n+1}}{hat{eta}(m{z}_{n+1})p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n+1} | m{z}_n)} (55)

因此,在LDS中, 離散變連續,求和變積分, hat{eta}(m{z}_n) 的遞推式為

c_{n+1}hat{eta}(m{z}_n) = int{hat{eta}(m{z}_{n+1})p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n+1} | m{z}_n)} dm{z}_{n+1} (20)

上式左右兩邊同乘以 hat{alpha}(m{z}_n) 得到

c_{n+1}gamma(m{z}_n) = int{hat{eta}(m{z}_{n+1})p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n+1} | m{z}_n)} hat{alpha}(m{z}_n)dm{z}_{n+1} (21)

重點關註上式被積式中除 hat{eta}(m{z}_{n+1}) 以外的三個高斯分布,我們需要對其進行一些變換。為了記號的簡潔,我們省略條件變數中的 m{x}_1, m{x}_2, cdots, m{x}_n ,從而可以看到以下推導過程:

egin{align} p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n+1} | m{z}_n) hat{alpha}(m{z}_n) & equiv p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n+1} | m{z}_n) p(m{z}_n) \ &= p(m{x}_{n+1} | m{z}_{n+1})p(m{z}_{n} | m{z}_{n+1}) p(m{z}_{n+1}) \ &= p(m{z}_{n+1} | m{x}_{n+1})p(m{z}_{n} | m{z}_{n+1}) p(m{x}_{n+1}) \ & equiv hat{alpha}(m{z}_{n+1}) p(m{z}_{n} | m{z}_{n+1}, m{x}_1, cdots, m{x}_n)c_{n+1} end{align} (22)

其中,第一個「 equiv」 符號之後省略了條件變數中的 m{x}_1, m{x}_2, cdots, m{x}_n ,第二個「 equiv 」符號之後又將省略的條件變數加進來了。將(22)式帶入(21)式從而得到

gamma(m{z}_n) = int{gamma(m{z}_{n+1}) } p(m{z}_{n} | m{z}_{n+1}, m{x}_1, cdots, m{x}_n) dm{z}_{n+1} (23)

現在的目標是要求  p(m{z}_{n} | m{z}_{n+1}, m{x}_1, cdots, m{x}_n) 。觀察(22)式第一個「 = 」之後的推導可以發現,  p(m{z}_{n} | m{z}_{n+1}, m{x}_1, cdots, m{x}_n) 其實是由 p(m{z}_{n+1} | m{z}_n) hat{alpha}(m{z}_n) 這兩個高斯分布經過貝葉斯定理變換之後得到的。記

 p(m{z}_{n} | m{z}_{n+1}, m{x}_1, cdots, m{x}_n) = N(m{z}_n | m{m}_n, M_n) (24)

利用《機器學習中的高斯分布》一文中的結論可得

m{m}_n = M_n (A^T Gamma^{-1} m{z}_{n+1} + V_n^{-1} m{mu}_n) (25)

M_n = (A^T Gamma^{-1} A + V_n^{-1})^{-1} (26)

將(24)~(26)式帶入(23)式,並再次利用《機器學習中的高斯分布》一文中的結論可得

hat{m{mu}}_n = M_n(A^T Gamma^{-1} hat{m{mu}}_{n+1} + V_n^{-1} m{mu}_n) (27)

hat{V}_n = (M_n A^T Gamma^{-1}) hat{V}_{n+1} (M_n A^T Gamma^{-1})^{T} + M_n (28)

這樣就得到了 gamma(m{z}_n) 的均值和協方差矩陣。但是這樣的結論看上去不夠簡潔,既不利於計算機實現(有很多多餘的中間矩陣運算),也無法從式子中直接看出遞推式的解釋。下面對該結論做一些簡化。

首先,根據矩陣求逆公式,

egin{align} M_n & = V_n - V_n A^T (Gamma + A V_n A^T)^{-1} A V_n \ & = V_n - V_n A^T P_n^{-1} A V_n end{align} (29)

那麼,

egin{align} M_n A^T Gamma^{-1}& = (V_n - V_n A^T P_n^{-1} A V_n)A^T Gamma^{-1} \ &=V_n A^T (I - P_n^{-1} A V_n A^T) Gamma^{-1} \ &= V_n A^T (I - P_n^{-1} A V_n A^T - P_n^{-1} Gamma + P_n^{-1} Gamma) Gamma^{-1} \ &= V_n A^T (I - P_n^{-1}P_n + P_n^{-1} Gamma) Gamma^{-1} \ &= V_n A^T P_n^{-1} end{align} (30)

其中用到了 P_n 的定義式(9)。記

J_n = V_n A^T P_n^{-1} (31)

那麼

egin{align} hat{m{mu}}_n &= J_n hat{m{mu}}_{n+1} + (V_n - V_n A^T P_n^{-1} A V_n) V_n^{-1} m{mu}_n \ &= J_n hat{m{mu}}_{n+1} + m{mu}_n - J_n A m{mu}_n \ &= m{mu}_n + J_n (hat{m{mu}}_{n+1} - A m{mu}_n) end{align} (32)

同理,

egin{align} hat{V}_n &= J_n hat{V}_{n+1} J_n^T + V_n - V_n A^T P_n^{-1} A V_n \ &= V_n + J_n (hat{V}_{n+1} - P_n) J_n^T end{align} (33)

其中用到了等式

A V_n= P_n J_n^T (34)

(32)式和(33)式即是最終的結論。

上面求出了 gamma(m{z}_n) 的均值和方差表達式,要用EM演算法學習LDS的參數,還要求 xi(m{z}_{n-1}, m{z}_n) ,它也是個高斯分布。根據定義,

egin{align} xi(m{z}_{n-1}, m{z}_n) &= p(m{z}_{n-1}, m{z}_n | m{x_1}, cdots, m{x_N}) \ &= p(m{z}_n | m{x_1}, cdots, m{x_N}) p(m{z}_{n-1} | m{z}_n, m{x_1}, cdots, m{x_N}) \ &= gamma(m{z}_n) p(m{z}_{n-1} | m{z}_n, m{x_1}, cdots, m{x_{n-1}}) \ &= N(m{z}_n | hat{m{mu}}_n, hat{V}_n) N(m{z}_n | m{m}_n, M_n) end{align} (35)

其中用到了條件獨立等式

p(m{z}_{n-1} | m{z}_n, m{x_1}, cdots, m{x_N}) = p(m{z}_{n-1} | m{z}_n, m{x_1}, cdots, m{x_{n-1}}) (36)

因為根據LDS模型結構,通過概率圖模型條件獨立判斷方法可知,當已知 m{z}_n 時, m{z}_{n-1}m{x}_n , cdots, m{x}_N 條件獨立。

注意到,(35)式的等式右邊正是(23)式右邊的被積式。利用《機器學習中的高斯分布》一文中的結論,高斯分布 xi(m{z}_{n-1}, m{z}_n) 的均值為

left [ hat{m{mu}}_{n-1}, hat{m{mu}}_n 
ight ] (37)

協方差矩陣為

left [ egin{array}{cc} hat{V}_{n-1} & (J_{n-1} hat{V}_n)^T \ J_{n-1} hat{V}_n & hat{V}_{n} end{array} 
ight] (38)

至此,兩個關鍵中間變數 gamma(m{z}_n)xi(m{z}_{n-1}, m{z}_n) 的分布就求出來了。而我們在應用EM演算法進行參數學習時,還需要求得以下幾個期望:

E[z_n] = hat{mu}_n (39)

E[z_n z_{n-1}^T] = hat{V}_n J_{n-1}^T + hat{mu}_n hat{mu}_{n - 1}^T (40)

E[z_n z_{n}^T] = hat{V}_n + hat{mu}_n hat{mu}_{n }^T (41)

至此,用EM演算法進行參數學習的準備工作都做完了。

在LDS中,共有六組參數需要學習,分別是初始概率分布的均值和協方差矩陣 m{mu}_0P_0 ,發射概率模型的變換矩陣和協方差矩陣 CSigma ,以及轉移概率模型的變換矩陣和協方差矩陣 AGamma 。我們將這些參數統一記為 	heta = [m{mu}_0, P_0, C, Sigma, A, Gamma]

根據EM演算法原理,需要最大化全數據的對數似然在(隱藏變數的)後驗分布下的期望,即最大化

Q(	heta, 	heta^{old}) = E_{Z|	heta^{old}}[ln{p(X, Z | 	heta)}] (42)

其中,全數據的對數似然為

egin{array} ln{p(X, Z | 	heta)} = &ln{p(z_1 | mu_0, P_0)} + sum_{n = 1}^N{ln{p(x_n | z_n, C, Sigma)}} \&+ sum_{n = 2} ^N{ln{p(z_n | z_{n-1}, A, Gamma)}} end{array} (43)

首先,考慮參數 m{mu}_0P_0 ,將其它所有參數當做常數。將 {p(z_1 | mu_0, P_0)} 的表達式帶入(43)式可得

egin{array}{} Q(	heta, 	heta^{old}) = &-frac{1}{2}ln{|P_0|} - E_{Z| 	heta^{old}}[frac{1}{2}(z_1 - mu_0)^T P_0^{-1} (z_1 - mu_0)] \ &+ const end{array} (44)

其中, const 表示所有與參數 m{mu}_0P_0 無關的項。最大化上式可求得

mu_o^{new} = E[z_1] (45)

V_0^{new} = E[z_1 z_1^T] - E[z_1]E[z_1^T] (46)

同樣的,考慮參數 AGamma ,將其它所有參數當做常數。將 p(z_n | z_{n-1}, A, Gamma) 的表達式帶入(43)式可得

egin{array}{} Q(	heta, 	heta^{old}) &= - frac{N - 1}{2}ln{Gamma} \ & - E_{Z|	heta^{old}}[frac{1}{2}sum_{n = 2}^{N}{(z_n - A z_{n-1} )^T Gamma^{-1} (z_n - Az_{n-1})}] + const end{array} (47)

其中, const 表示所有與參數 AGamma 無關的項。最大化上式可求得

A^{new} = (sum_{n = 2}^{N}{E[z_n z_{n - 1}^T]}) ({sum_{n = 2}^{N}{E[z_{n-1}z_{n-1}^T]}})^{-1} (48)

egin{array}{} Gamma^{new} = frac{1}{N-1} sum_{n = 2}^{N} { left { E[z_n z_n^T] - A^{new} E[z_{n-1} z_n^T] \ - E[z_n z_{n-1}^T] (A^{new})^T + A^{new} E[z_{n-1}z_{n-1}^T] (A^{new})^T 
ight }} end{array} (49)

最後,考慮參數 CSigma ,將其它所有參數當做常數。將 p(x_n | z_n, C, Sigma) 的表達式帶入(43)式可得

egin{array}{} Q(	heta, 	heta^{old}) &= -frac{N}{2} ln{|Sigma|} \&- E[frac{1}{2} sum_{n = 1}^{N}{(x_n - C z_n)^T Sigma^{-1}(x_n - C z_n)}] + const end{array} (50)

其中, const 表示所有與參數 CSigma 無關的項。最大化上式可求得

C^{new} = (sum_{n = 1}^{N}{x_n E[z_n^T]}) (sum_{n = 1}^{N}{E[z_n z_n^T]})^{-1} (51)

egin{align} Sigma^{new} = frac{1}{N} sum_{n = 1}^{N} left { {x_n x_n^T - C^{new} E[z_n] x_n^T - x_n E[z_n^T](C^{new})^T + C^{new} E[z_nz_n^T]}(C^{new})^T 
ight } end{align} (52)

此文完。


推薦閱讀:

自動控制原理筆記5現代控制理論基礎
自控貓|[問答]連載01-可燃有毒氣體報警控制器能否與DCS合併?
[No Time to Explain] - 知乎最熱鬧的控制群
自動控制總結:第三章:線性系統的時域分析
自動控制原理筆記2

TAG:自動控制 | 卡爾曼濾波KalmanFilter | 機器學習 |