Stationary Kalman Filter

Question

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.

So, here am I asking:

• 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?

Answer 1

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.