Stationary Kalman Filter

Ecercise 4.5 from Bayesian Filtering & Smoothing by Simo Särkkä:

Derive the stationary Kalman filter for the Gaussian random walk model. That is, compute the limiting Kalman filter gain when $$k \rightarrow \infty$$ and write down the mean equation of the resulting constant-gain Kalman filter. Plot the frequency response of the resulting time-invariant filter. Which type of digital filter is it?

The Gaussian random walk model comes from Example 4.1:

\begin{align} x_k &= x_{k-1} + q_{k-1}, \qquad q_{k-1} \sim \mathcal N(0,Q)\\ y_k &= x_k + r_k, \quad\;\;\;\qquad r_{k} \sim \mathcal N(0,R) \end{align}

with $$Q=R=1$$.

The Kalman filter prediction and update equations for this particular case are given on page 59 (Example 4.2), or, in more general form on page 57 (Section 4.3).

Relevant posts: How to derive the stationary Kalman filter predictor? This left me with more questions than answers.

My understanding is that for the Kalman filter to be stationary, the filtered mean and variance should not change over time. But I don't understand how is that possible with $$k \rightarrow \infty$$ other than when measurements repeat themselves.

• What is stationarity in state space models?
• Why would the Kalman gain converge to a constant with the number of observations?
• If the Kalman gain is constant, then so must be the variance P, right?

According to the random walk model $$A=H=1$$.

Using this let's derive the equations (assume $$Q$$ and $$R$$ are some constants):

\begin{align} m^-_k &= A_{k-1} m_{k-1} = m_{k-1} \\ P^-_k &= A_{k-1} P_{k-1} A^T_{k-1} + Q_{k-1} = P_{k-1} + Q\\ \\ v_k &= y_k - H_k m^-_k = y_k - m^-_k\\ S_k &= H_{k} P^-_{k} H^T_{k} + R_{k} = P^-_{k} + R\\ \\ K_k &= P^-_{k} H^T_{k} S^{-1}_{k} = \frac{P^-_{k}} {P^-_{k} + R}\\ \\ m_k &= m^-_k + K_k v_k = m^-_k + \frac{P^-_{k}} {P^-_{k} + R} (y_k - m^-_k)\\ P_k &= P^-_k - K_k S_k K^T_k = P^-_k - \frac{P^-_{k}} {P^-_{k} + R} (P^-_{k} + R) \frac{P^-_{k}} {P^-_{k} + R} = P^-_k - \frac{(P^-_{k})^2} {P^-_{k} + R} \end{align}

We can rewrite the last equation like this:

\begin{align} P_k &= P_{k-1} + Q - \frac{(P_{k-1} + Q)^2} {P_{k-1} + Q + R} \end{align}

Now we can see, that $$P_k$$ depends only on itself and on two constants $$Q$$ and $$R$$. At some point it converges to some constant value, let's say $$P^*$$.

It means

\begin{align} P^* &= P^* + Q - \frac{(P^* + Q)^2} {P^* + Q + R} \end{align}

With some elementary math we get a normal quadratic equation:

\begin{align} (P^*)^2 + Q (P^*) - Q R &= 0 \end{align}

By solving this equation you get two roots, one of which is negative and does not fit the restriction for not negative covariance.

So the only solution is:

\begin{align} P^* &= \frac{-Q + \sqrt{Q^2 + 4QR}}{2} \end{align}

By substituting this into expression for $$K_k$$ you get some another constant value for the Kalman gain.

Here is a simple Matlab code to check the equations:

function [] = main()

dt = 0.01;
t=(0:dt:1)';

n = numel(t);

signal_ref = zeros(size(t));

% system noise
Q = 1;

for i=2:numel(t)
signal_ref(i) = signal_ref(i-1) + randn()*sqrt(Q);
end

% measurement noise
R = 1;

signal = signal_ref + randn(size(t)).*sqrt(R);

% state
X = 0;

% covariance matrix
P = 10;

% transition matrix
F = 1;

% observation matrix
H = 1;

% kalman filter output through the whole time
X_arr = zeros(n, 1);
K_arr = zeros(n, 1);
P_arr = zeros(n, 1);

for i = 1:n

y = signal(i);

if (i == 1)
[X] = init_kalman(y); % initialize the state using the 1st measurement
else
[X, P] = prediction(X, P, Q, F); %Prediction

[X, P, K] = update(X, P, y, R, H); %Update

K_arr(i, :) = K;
P_arr(i, :) = P;
end

X_arr(i, :) = X;

end

figure;
subplot(2,1,1);
plot(t, signal_ref, 'LineWidth', 2);
hold on;
plot(t, signal, 'LineWidth', 2);
plot(t, X_arr(:, 1), 'LineWidth', 2);
hold off;
grid on;
legend('Ground Truth', 'Sensor', 'Estimation');

subplot(2,1,2);
plot(t, K_arr, 'LineWidth', 2);
hold on;
plot(t, P_arr, 'LineWidth', 2);
hold off;
grid on;
legend('Kalman gain', 'Covariance P');
end

function [X] = init_kalman(y)
X = y;
end

function [X, P] = prediction(X, P, Q, F)
X = F*X;
P = F*P*F' + Q;
end

function [X, P, K] = update(X, P, y, R, H)
Inn = y - H*X;
S = H*P*H' + R;
K = P*H'/S;

X = X + K*Inn;
P = P - K*H*P;
end


For $$Q = 1$$ and $$R = 2$$ the stationary solution looks like this:

For $$Q = 1$$ and $$R = 1$$ both Kalman gain and covariance converge to the same value:

UPDATE

Regarding the frequency response.

We can rewrite the equation for the state like this:

\begin{align} m_k &= m_{k-1} + K (y_k - m_{k-1}) = (1 - K)m_{k-1} + K y_k\\ \end{align}

This is a discrete equation for a low pass filter (see Wikipedia). The $$K$$ would be the smoothing factor, which relates to the time constant and to the cut-off frequency.

I would say the frequency responce of your filter is the same as for the low pass filter and looks like this (picture taken from Wikipedia):

This is my assumption, I'm not really deep in this stuff, but it seems to me plausible.

• Thank you! Shouldn't you have P = P - K*S*K' there in the update function? You've plotted that to show that P and K converge, right? But, what would "plot the frequency response of the resulting time-invariant filter" mean? Does it mean that I need to use the stationary K to update m and P? – Sandu Ursu Oct 17 '19 at 15:47
• Also, the formula for K_k needs to have S_k^{-1}. – Sandu Ursu Oct 17 '19 at 15:54
• Not obvious to me how can I extend this to the case when I'm dealing with matrices instead of scalars. – Sandu Ursu Oct 17 '19 at 15:59
• I corrected the formula for K_k. Regarding your first comment you can have a look at formula (4.29) in the book. This expression is an alternative formula for (4.21) – Anton Oct 17 '19 at 20:12
• I don't think it would be easy to show the same derivation in matrix form. Scalars give us a better intuition to understand the Kalman filter. I'm sure it works for matrices as well. Ecercise 4.5 deals with scalars, otherwise the R and Q matrices would have been given as identity matrices and not as scalars. Don't you think so? – Anton Oct 17 '19 at 20:18