# How to initialize observation Matrix in Kalman Filter when there is no clear relationship between measurement and state?

I am try to use Linear Kalman to do time series prediction. I understand that I have to define a model process matrix which indicate how system state evolve, and a measurement matrix H which convert state variable to measurement space. However in my problem, there is no clear relationship between my measurement and state, so there is no way to just give 1 to some element in H to pick up the same physical meaning in state variable X, and 0 to other irrelevant component. Under this situation, how do I design H matrix? I have read a paper, it propose method below which use some data to train all the matrix, including H. But I don't understand how to deduce? can some one shed light on this?

• Can you express your measurement as a function of the state? I.e., if $y_k = x_{1,k} + 4 x_{2, x}$, that's a function. Ditto if $y_k = \sin x_{1, k} + \cos x_{2, k}$. If so, please edit your question with this function. Mar 4 at 3:25
• Well, to use the formula from the paper you need a recorded time series of your state vector $x$ and the corresponding measurements $y$. With those you can then construct Matrix $H$ with which you would have to multiply your measurements to get the state vector. Can you provide a link to the paper? Mar 4 at 20:24
• @MatthiasLa Here is the paper: mitpressjournals.org/doi/abs/10.1162/089976606774841585, titled as 'Bayesian Population Decoding of Motor Cortical Activity Using a Kalman Filter'. There is a lot domain discussion there, but not too relerant to the deduction. So you only need to focus on the supplementary section at the end of the paper. Mar 5 at 9:44

In order to get a better understanding I used the method proposed in the paper to build the measurement function $$y=h(x)$$ using a Matrix $$H$$ and the state transition function $$x_{n+1}=f(x_n)$$ using a Matrix $$A$$. I did this for a practical example to get a feeling how it works.

To be clear in the paper you cited they call the measurements $$z$$ instead of $$y$$. So their measurement function is:

Their state transfer function is:

If the physical relationship between the measurement $$z$$ and the state $$x$$ is unknown one can use the formulas from the paper to calculate $$H$$ and $$A$$ depending on past states and measrements.

The example I tried is from Plane with constant speed. A radar system measures the distance to a flying plane that travels with constant speed. Starting position is 30 km away. Starting speed is 40 m/s. The state Vector has two values: position and speed of the airplane. So the inital state would be: $$x_1=[30000, 40]$$. Because the radar only measures the distance, the measurement equation simply is: $$z_n = \begin{pmatrix} x_{1,n} \\ 0 \end{pmatrix}$$. Thats because only the distance ($$x_{n,1}$$) gets measured but not the velocity ($$x_{n,2}$$).

The state transition equation is $$\begin{pmatrix} x_{1,n+1} \\ x_{2,n+1}\end{pmatrix} = \begin{pmatrix} x_{1,n} + \Delta t \cdot x_{2,n} \\x_{2,n}\end{pmatrix}$$. The velocity $$x_{n,2}$$ stays the same between step $$n$$ and $$n+1$$ but the distance changes by the multiplication of the time step $$\Delta t$$ and the speed during that time step $$x_{2,n}$$.

Using those equations and the measurement values provided in the example I get the following:

Next we try to achieve the same without writing the functions $$f$$ and $$h$$ explicitly but we try to use the formulas from the paper you cited.

For that also the recorded states $$x_{n}$$ and $$x_{n+1}$$ are needed. Luckily they are also provided in the example:

Using the Formulas for $$H$$ and $$A$$ and the recorded traces it is possible to calculate them.

$$H$$ in the case of our example becomes a $$1\times2$$ Matrix. Thats because $$z$$ is a $$1\times10$$ vector and $$x_n$$ is a $$2\times10$$ vector. As the latter is transposed $$H$$ ends up being $$1\times2$$.

Quite the same way $$A$$ ends up being a $$2\times2$$ matrix. $$x_{n+1}$$ and $$x_n$$ are both a $$2\times10$$ matrix. As the latter is transposed by the multiplication one ends up with a $$2\times 2$$ matrix.

Now we can rewrite the equations $$f$$ and $$h$$.

$$x_{n+1}=f(x_n): x_{n+1}=A\cdot x_n$$ and $$z_n=h(x_n): z_n=H\cdot x_n$$

If I put in the recorded values I get $$H=\begin{pmatrix} 0.956 & 35.38 \end{pmatrix}$$ and $$A=\begin{pmatrix} 1.00 & 4.98 \\ 0.00 & 1.00 \end{pmatrix}$$.

With that I rerun the Kalman filter and I get the following:

This looks very similar to the former result which was derived using a physical modell. So in this case the method from the paper seems to work quite fine.

I did this using Matlab. So in case you are interested in the code:

close all

initialStateGuess = [30000 40];
n=length(initialStateGuess);
dt=5;

x_state=[30182  30351.4 30573.3 30769.5 31001.5 31176.4 31333.2 31529.4 31764.3 31952.9;38.2    36  40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
x_state2=[30373 30531.6 30774.3 30968.1 31216.8 31371.5 31509.2 31715.4 31974.8 32152.4; 38.2   36.0    40.2    39.7    43.1    39  35.2    37.2    42.1    39.9];
yMeas = [30110 30265 30740 30750 31135 31015 31180 31610 31960 31865];% Measurement Values
H=(yMeas*x_state2')/(x_state2*x_state2');
A=(x_state2*x_state')/(x_state*x_state');

% h=@(x)[x(1)];% Measurement function
% f=@(x)[x(1)+dt*x(2);x(2)];% State transition function

h=@(x)[H(1)*x(1)+H(2)*x(2)];
f=@(x)[A(1,1)*x(1)+A(1,2)*x(2);A(2,1)*x(1)+A(2,2)*x(2)];

ukf = unscentedKalmanFilter(...
f,... % State transition function
h,... % Measurement function
initialStateGuess,...

yTrue=linspace(30000,32000,10);
R = var(yMeas-yTrue);%0.2; % Variance of the measurement noise v[k]

Q=var(x_state-yTrue);
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 30;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
% Let k denote the current time.
%
% Residuals (or innovations): Measured output - Predicted output
e(k) = yMeas(k) - h(ukf.State); % ukf.State is x[k|k-1] at this point
[xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
predict(ukf);
end

plot([5:5:50],xCorrectedUKF(:,1))
hold on
plot([5:5:50],yMeas)
plot([5:5:50],yTrue)
legend('Kalman FIlter state','Measurement Values','True Value')
grid on

• Remarkable work. Thank you very much. You just clear my doubt of such matrix initialization method. Now you just verify it works. I gauss another thing I want to add is how to deduce the formula. In the formula, A is actually the overdeterminent equation answer under the least square error sense, and W is the residual after the fitting of A using pseudo-inverse division. The same procedure applys to the calculation of H and Q. It's actually quite intuitive, it's strange why a lot of book haven't mentioned this method. Mar 5 at 17:50
• @XiaolongWu I was also surprised how well it works in this case. Actually if you look at the matrices they almost exactly replicate the former equations. $H$ results in $~1\cdot x_{n,1}$ and $~36\cdot x_{n,2}$. The latter is almost $40$ which is the constant speed. With $A$ it is similar. It's $1 \cdot x_{n,1} + ~5 \cdot x_{n,2}$. As $\Lambda t$ ist 5 it's exactly the formula from beforehand. However this of course only works for linear problems and only if some measruements for both $z$ and $x$ are prerecorded. Mar 6 at 11:35