Nonlinear filtering: Grid-based filter

The process of Bayes filtering requires to solve integrals, that are in general intractable. One approach to circumvent this problem is the use of grid-based filtering. In this article, we will derive this method directly from the recursive equations of the Bayes filter. This marks the first part of the nonlinear filtering series. If you haven’t read the introduction, I would recommend to read it first. Before we dive into the derivation, let’s try to state the main idea behind grid-based filtering.

The grid-based filter approximates the Bayes filter by restricting and discretizing the state, input and observation space, to obtain finite domains.

Derivation

We will start the derivation directly from the recursive equations of the Bayes filter with the prediction step

$$p(x_{t+1}|y_{0:t},u_{0:t}) = \int_{x_{t}} p(x_{t+1}|x_{t}, u_{t})p(x_{t}|y_{0:t},u_{0:t-1}) dx_{t}$$

and the update step

$$p(x_t|y_{0:t},u_{0:t-1}) = \frac{p(y_t|x_t)p(x_t|y_{0:t-1},u_{0:t-1})}{\int_{x_t}p(y_t|x_t)p(x_t|y_{0:t-1},u_{0:t-1}) \,dx_t} .$$

The first and most important step in order to arrive at the equations of the grid-based filter is to discretize our system and observation model. Given a set of discretization points \((x^1,\ldots,x^X)\), \((y^1,\ldots,y^Y)\) and \((u^1,\ldots,u^U)\) the result of the discretization can be written as

$$\hat{p}(x_{t+1}|x_{t}, u_{t}) = \frac{1}{Z}p(x_{t+1}|x_{t}, u_{t})\sum_{k=1}^X\sum_{l=1}^X \sum_{m=1}^U \delta_{x^k}(x_{t+1})\delta_{x^l}(x_t)\delta_{u^m}(u_t)$$

for the system model and

$$\hat{p}(y_t|x_t) = \frac{1}{Z}p(y_{t}|x_{t})\sum_{k=1}^Y\sum_{l=1}^X \delta_{y^k}(y_{t})\delta_{x^l}(x_t)$$

for our observation model, where \(Z\) is a normalization factor.

These equations can look confusing and random. Especially if you are not familiar with the Dirac delta function. This box will give a short introduction to the Dirac delta function and provide the tools you will need. Warning: The Dirac delta function is actually not a function, but a distribution or a measure!

The Dirac delta function \(\delta(x)\) can be imagined as a function, that is infinite at \(x=0\) and zero at \(x \neq 0\) with the property

$$\int_x \delta(x)\,dx = 1.$$

It is not possible to plot the Dirac delta function, but we can visualize it schematically by using an arrow pointing upwards, as shown in the figure below.

The length of the line represents the area under the function. The expression \(\delta(x-y)\) is a shifted version of the Dirac delta function with its peak at \(y\).

To unclutter the notation, we will use the shorthand \(\delta_y(x)\) for \(\delta(x-y)\).

The sifting property

$$\int_x f(x)\delta_y(x)\,dx = f(y)$$

will play an important role in our derivation. Intuitively, it is sifting out the function value of \(f(x)\) at \(y\).

The multivariate Dirac delta function, which can be loosely thought of as infinite at \((x_1,\ldots,x_n)=(0,\ldots,0)\) and zero at \((x_1,\ldots,x_n) \neq (0,\ldots,0)\) is defined as the product of several Dirac delta functions:

$$\delta(x_1,\ldots,x_n) = \delta(x_1)* \ldots *\delta(x_n).$$

Several Dirac delta functions can also be combined to form a comb or grid, which is shown in the following figure.

To obtain such a grid, the Dirac delta functions are simply added together

$$\sum_{i=1}^N\delta_{x_i}(x).$$

Finally, we can multiply the grid with a function \(f(x)\) :

$$\sum_{i=1}^Nf(x)\delta_{x_i}(x).$$

The particular Dirac delta functions are weighted by the functional value at the corresponding point. This can be visualized as following.

Now that we have a better understanding of the Dirac delta function, let’s look again at our system and observation model. In our system model, we find a multivariate Dirac delta function

$$\delta_{x^k}(x_{t+1})\delta_{x^l}(x_t)\delta_{u^m}(u_t),$$

where the particular Dirac delta functions could be multivariate as well. We also note, that the multivariate Dirac delta function is embedded in a triple sum over the corresponding discrete spaces. It describes, therefore, a grid of Dirac delta functions. Finally, we multiply this grid with our system model. To obtain a proper probability distribution, we normalize our distribution with \(Z\). The observation model can be treated likewise.

If we start the recursion with a discretized prior distribution

$$\hat{p}(x_0) = \sum_{k=1}^X w^k \delta_{x^k}(x_0)$$

all of our belief distributions will be discretized as well. Therefore, the posterior will have the form

$$\hat{p}(x_{t}|\cdot) = \sum_{k=1}^X w^k \delta_{x^k}(x_t).$$

Prediction step

Let’s see what happens, if we plug the discretized version of our system model into the equations of the prediction step:

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \int_{x_{t}} \frac{1}{Z}p(x_{t+1}|x_{t}, u_{t})\sum_{k=1}^X\sum_{l=1}^X\sum_{m=1}^U \delta_{x^k}(x_{t+1})\delta_{x^l}(x_t)\delta_{u^m}(u_t)\sum_{n=1}^X w^n \delta_{x^n}(x_t) dx_{t}$$

First of all, we know our current input is \(u^i\). To evaluate \(\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) \) at \(u^i\), we simply have to compute \(\int_{u_t}\hat{p}(x_{t+1}|y_{0:t},u_{0:t})\delta_{u^i}(u_t)\,d_t \). As a result, we will obtain

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \frac{1}{Z}\int_{x_{t}}p(x_{t+1}|x_{t}, u_t=u^m) \sum_{k=1}^X\sum_{l=1}^X \delta_{x^k}(x_{t+1})\delta_{x^l}(x_t)\sum_{n=1}^X w^n \delta_{x^n}(x_t) dx_{t}$$

The following expression

$$\sum_{l=1}^X \delta_{x^l}(x_t)\sum_{n=1}^X w^n \delta_{x^n}(x_t)$$

is not zero only if \(l = n\). Therefore, we can replace this expression with

$$\sum_{n=1}^X w^n \delta_{x^n}(x_t)$$

and obtain

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \frac{1}{Z}\int_{x_{t}}p(x_{t+1}|x_{t}, u_t=u^m) \sum_{k=1}^X\sum_{n=1}^X \delta_{x^k}(x_{t+1}) w^n \delta_{x^n}(x_t) dx_{t} .$$

By rearranging the integral and sums we obtain

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \frac{1}{Z}\sum_{k=1}^X\delta_{x^k}(x_{t+1})\int_{x_{t}}p(x_{t+1}|x_{t}, u_t=u^m) \sum_{n=1}^X w^n \delta_{x^n}(x_t) dx_{t}$$

Using the sifting property we will arrive at

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \frac{1}{Z}\sum_{k=1}^X\delta_{x^k}(x_{t+1})\sum_{n=1}^Xp(x_{t+1}|x_{t} = x^n, u_t=u^m) w^n$$

Update step

Let’s start our treatment of the update step and by looking at the numerator. Again, we plug in our discretized distributions to obtain

$$\hat{p}(y_t|x_t)\hat{p}(x_t|y_{0:t-1},u_{0:t-1}) = \frac{1}{Z}p(y_{t}|x_{t})\sum_{k=1}^Y\sum_{l=1}^X \delta_{y^k}(y_{t})\delta_{x^l}(x_t)\sum_{m=1}^X w^m \delta_{x^m}(x_t).$$

We know, that we have an observation \(y^k\), which we can plug in in the same way as the input above:

$$\hat{p}(y_t|x_t)\hat{p}(x_t|y_{0:t-1},u_{0:t-1}) = \frac{1}{Z}p(y_{t}=y^k|x_{t})\sum_{m=1}^Xw^m\delta_{x^m}(x_t)$$

The denominator is simply the integral over \(x_t\) of the expression above:

$$\int_{x_t}\hat{p}(y_t|x_t)\hat{p}(x_t|y_{0:t-1},u_{0:t-1}) \,dx_t = \frac{1}{Z}\int_{x_t}p(y_{t}=y^k|x_{t})\sum_{m=1}^Xw^m\delta_{x^m}(x_t) \,dx_t.$$

With help of the sifting property, we obtain the final expression

$$\int_{x_t}\hat{p}(y_t|x_t)\hat{p}(x_t|y_{0:t-1},u_{0:t-1}) \,dx_t = \frac{1}{Z}\sum_{m=1}^Xw^mp(y_{t}=y^k|x_{t}=x^k)$$

for the denominator.

The full update step is then defined as

$$\hat{p}(x_t|y_{0:t},u_{0:t-1}) = \frac{p(y_{t}=y^k|x_{t})\sum_{m=1}^Xw^m\delta_{x^m}(x_t)}{\sum_{l=1}^Xw^mp(y_{t}=y^k|x_{t}=x^k)} .$$

Let’s summarize our results!

Grid-based filter

The recursive formula for the grid-based filter in state space models consists of the prediction step

$$\hat{p}(x_{t+1}|y_{0:t},u_{0:t}) = \frac{1}{Z}\sum_{k=1}^X\delta_{x^k}(x_{t+1})\sum_{n=1}^Xp(x_{t+1}|x_{t} = x^n, u_t=u^m) w^n$$

and the update step

$$\hat{p}(x_t|y_{0:t},u_{0:t-1}) = \frac{p(y_{t}=y^k|x_{t})\sum_{m=1}^Xw^m\delta_{x^m}(x_t)}{\sum_{l=1}^Xw^mp(y_{t}=y^k|x_{t}=x^k)} .$$

The recursion is started with a discrete prior distribution over the initial state \(\hat{p}(x_0)\).

Discrete probability distribution

Our discretized system and observation model is non-zero only at the points on the grid. Nonetheless, the model itself is still continuous: we can evaluate the function at points, which are not part of the grid. Then again, the posterior will always remain on the grid during the filtering process. Thus, we can represent our system as a discrete probability distribution

$$P(x_{t+1}=k|x_t=l,u=m) = \frac{1}{Z(l,m)} p(x_{t+1}=x^k|x_{t}=x^l, u_t=u^m)$$

with the normalization factor \(Z(l,m) = \sum_{k=1}^X p(x_{t+1}=x^k|x_{t}=x^l, u_t=u^m) \).

Similarly, the observation model can be represented by the discrete probability distribution

$$P(y_{t}=k|x_t=l) = \frac{1}{Z(l)} p(y_{t}=y^k|x_{t}=x^l)$$

with the normalization factor \(Z(l) = \sum_{k=1}^X p(y_{t}=y^k|x_{t}=x^l) \). With this representation, we arrive at the Bayes filter for discrete probability distributions.

Discrete Bayes filter

The recursive formula for the discrete Bayes filter in state space models consists of the prediction step

$$P(x_{t+1}|y_{0:t},u_{0:t}) = \sum_{x_{t}} p(x_{t+1}|x_{t}, u_{t})p(x_{t}|y_{0:t},u_{0:t-1})$$

and the update step

$$P(x_t|y_{0:t},u_{0:t-1}) = \frac{P(y_t|x_t)P(x_t|y_{0:t-1},u_{0:t-1})}{\sum_{x_t} P(y_t|x_t)P(x_t|y_{0:t-1},u_{0:t-1})} .$$

The recursion is started with a discrete prior distribution over the initial state \(P(x_0)\).

Example

Enough of the dry theory! Let’s play around with the grid-based filter in our race track example.

Backward

No action

Forward

Predict step

Observe

Update step

On the outside of the race track, you will notice a blue colored strip. This strip represents our current posterior of the current position of the race car. At the beginning, we have no knowledge about the position of the race car and assign uniform probability over all positions. By pressing the OBSERVE button two things will happen: first, we will take a measurement of the distance to the tree and second, we will display the likelihood for this observed distance on the brown strip inside the race track. By pressing the UPDATE STEP button, we will perform our update step and show the resulting posterior at the outer strip. You will note, that both strips will have the same form after the update. The reason is simple: we just multiplied our likelihood with a constant vector and normalized afterward. Now we are ready for the next time step. Take an action, by pressing the corresponding button below the race track. After the step is performed, you have to update your posterior by pressing the PREDICT STEP button. You will see that the outer strip will change accordingly. Now we finished one full cycle of the filtering process and are ready to start a new cycle by taking a measurement.

What if our distance meter is not working anymore? By either clicking on the tree or pressing the W button on your keyboard, you can turn off your measurement device. Therefore, the observation and update step will be skipped. The tree will become opaque, if your measurement device is turned off.

With the slider below the race track, you can choose a grid size of the discrete probability models. If you want to reset the environment, just press the reset button in the bottom left corner or press the R button on your keyboard. As before you can control the car by using your keyboard: A (Backward), S (Stop), D (Forward) or the buttons below the race track.

Still feeling hungry for Bayes filters? Then you should definitely check out the next part of the nonlinear filtering series covering the derivation of the extended Kalman filter. See you there!

Acknowledgement

The vector graphics of the car were created by Freepik.