Kalman Filter : MMSE Estimator

常常說 Kalman filter 是 optimal estimate (in linear dynamic model), 如何解釋? 或者更基本，Kalman filter 的理論基礎為何? Kalman filter 和一般常用的 LPF (Butterworth, Cheby filter, etc.) 有什麼本質上的差異?

Kalman filter 和常用 LPF 差異

一般信號處理常用的 filter (FIR, IIR, Butterworth, Cheby) based on Laplace or Fourier tranform , (i) 都是 linear time invariant (LTI) filter;  (ii) filter 的特性 (e.g.頻寛)和輸入信號基本上 independent.

要計算輸出信號，我們通常先 derive (or define) filter 的 transfer function and impulse response, 再用 frequency domain (s/z-domain) multiplication 或 time domain convolution 來獲得。

Kalman filter, on the other hand, 是由 stochastic signal (random sequence 隨機信號) estimation 出發 (i) 輸入信號和相關雜訊的 statistics 決定了 Kalman filter 的特性;  (ii) Kalman filter 基本上是 time varying filter.  (在很多情況下 Kalman filter 會收歛到 steady state filter, 近似 time invariant filter).  以上兩點和常用 LPF 完全不同。

Kalman filter 和常用 LPF 可以類比嗎?

前文 ( Kalman filter for Engineer) 提到 Kalman filter 可以視為一個 feedback controller with time varying Kalman gain (Kg).  因此可以視為一個 time varying IIR filter, 頻寛會隨時間改變。

在 steady state (如果 steady state 存在)時，Kg 基本上近似常數，趨近一個 time invariant filter).  不過 Kg 由輸入信號和相關雜訊決定，和一般 LPF 仍有不同。

Kalman filter 和 Wiener Filter 的差異

可以 reference Binghamton 的 slides.

Wiener filter assumes WSS (wide sense stationary) signal + noise

Kalman filter assumes Dynamic model + observation noise

Wiener filter 是 time invariant filter, 但 filter 特性由輸入信號和雜訊決定，這點和 Kalman filter 相似。

If we let consider only after much time has elapsed: given IIR Wiener case and steady state Kalman filter, then Kalman = Wiener.

Step1: Random Variable MMSE Estimator

先從最簡單的 one variable 開始，X 的 MMSE (Minimum Mean-Square Error) estimator 是 find x* to minimize E(|X-x*|^2).  很容易可証明: x* = E(X),  MSE=E(|X-E(X)|^2)=σ^2.        Any other estimate E(X)+a 的 MSE=σ^2 + a^2 > σ^2.

當有 two variables X, Y, the MMSE (Minimum Mean-Squared Error) estimator:
Let y be observed.   Fnd x* an estimate of the unobserved x to minimize  E(|X-x*|^2).  Hint: E(|X-x*|^2) = Ex Ey (|X-x*|^2|y)

很明顯如果 X, Y are independent or uncorreleated. y is observed or not 對 X estimation 並沒有任何幫助:

x* = E(X) and MSE = σ_x^2.

如果 X and Y are correlated, y is observed 對 X estimation 多少有幫助。可以証明 x* = E(X|Y=y) 是 MMSE estimator. 並且 MSE error 會比沒有 y is observed 小。直覺上 make sense, 多了 y 的資訊對減小 X esimation error 多少有幫助。

注意 x* = E(X|Y=y) 並沒有假設任何 X and Y 的統計特性 (如 Gaussian)!  這是 Baysian estimator 的基礎。

X, Y Gaussian

再來假設 X, Y Gaussian 來進一步說明。

[x, y] ~ N( [ E(X), E(Y)], [σ_x^2, σ_xy; σ_xy, σ_y^2])

x* = E(X|Y=y) = E(X) + σ_xy (y-E(y)) / σ_y^2

MSE = σ_x^2 – |σ_xy|^2/σ_y^2  < σ_x^2

(i) x* = E(X|Y=y) = a + b y  是一個 linear function of y

(ii) MSE is reduced with additional observable y

(iii) If X, Y are uncorrelated, σ_xy = 0.  x* = E(X); MSE = σ_x^2 as expected

X, Y Non-Gaussian

在 X and Y 非 Gaussian 情況下要得到 x* = E(X|Y=y) 就比較困難。上述 Gaussian 結論需要修正:

(i) x* = E(X|Y=y) = g(y) in general in NOT a linear function of y.  不過可以証明在所有的 linear function (a + b y) 中，之前 Gaussian 所得的公式:

x* = E(X) + σ_xy (y-E(y)) / σ_y^2

MSE 是最小的，稱為 LMSE and LMMSE (Linear Minimum Mean Square Error).   同時可以証明 LMMSE 的確是 minimum (at least local minimum) and unbias.

但是可能存在其他非線性的  x* = g(y) 的 MSE 更小 (global minimum), some example are here.

(ii) and (iii) 仍然是正確。

Step2: MMSE with linear measurement

Consider the following special case (Y is linear transformation x plus another Gaussian noise)

y = H x + v,  x~N(E(X), σ_x^2),  v~N(E(V), σ_v^2)

x* = EX+ K (y-EY)   where EY = H EX + EV, 

K = σ_x^2 H’ (H σ_x^2 H’ + σ_v^2) ^-1  from equation of last section.  K 事實上就是 Kalman gain.

MSE = σ_x^2 – K * H σ_x^2 = σ_x^2 (1 – K H)   Therefore, MSE is reduced due to y is observed.

The above equation forms the base of Kalman filter. 

 

Step3: MMSE with sequences of linear measurement

y = Hx + v    given y1, y2, y3, … yk

x*(k) = E(X|y1, y2, y3, ..yk) =  E(X|y1|y2| .. |yk)

-> x*(1) = E(X|y1) = EX + K(1) (y1 – EY)

EY=H EX+EV  MSE(1)=P(1)=sig_x^2 H’(H sig_x^2 H’+ sig_v^2)^-1

K(1) = sig_x^2 H’ (H sig_x^2 H’ + sig_v^2) ^-1

-> x*(2) = E(X|y1|y2) = E(X|y1) + K(2)(y2 – EY(2))

                = x*(1) + K(2)(y2 – EY(2))

E(X|y1)=x*(1)     EY(2)=Hx*(1)+EV   MSE(2)=P(2)=P(1)^2 H’(HP(1)^2H’+sig_v^2)^-1  K(2) = sig_x^2 H’ (H sig_x^2 H’ + sig_v^2) ^-1

-> x*(3) = x*(2) + K(2)(y2 – EY(2))

其實上述公式就是 Kalman filter 的雛型。只要再加上 state  dynamic model  x(k+1) = A x(k) + W  就變成 Kalman filter, 這部份先帶過。  反之 Kalman filter 的 state equation 簡化為 x(k+1) = x(k), 就變回 sequencial MMSE estimator. 

Sequencial MMSE estimation 和 Kalman filter 一個很大的不同是沒有 W (state process noise).  因此愈多的 measurements,  MSE error 會一直減小，最後會接近 0, 也不會有 steady state. 

在 Kalman filter 有 W (state process noise).  在 prediction stage 時，MSE noise 會增加，update stage 時，MSE noise 才會減少。因此 Kalman filter MSE 最後會收歛到一個非零的 steady state. 圖示如下。

 

MMSE with nonlinear measurement

y = f(x) + v

直覺的 approach, 就是用 Taylor expansion 來 approximate f(x).   However, there are two problems:

(i) In general, y is no longer Gaussian 上述的公式不適合。

(ii) For some common nonlinear function f(x) may introduce bias to disturb the estimation.

除了 Taylor expansion 外，另一方法是用 Gaussian 近似 f(x) 的 pdf.  一個方法是 unscented transform (UT).  之後再介紹。