# 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, 2021 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, 2021 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, 2021 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, 2021 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, 2021 at 11:35