Computing Probabilistic Information States

The probabilistic I-states can be quite complicated in practice because each el- ement of prob is a probability distribution or density function. Therefore, sub- stantial effort has been invested in developing efficient techniques for computing probabilistic I-states efficiently. This section can be considered as a continua- tion of the presentations in Sections 11.2.3 (and part of Section 11.4, for the case of continuous state spaces). Section 11.6.1 covers Kalman filtering, which pro- vides elegant computations of probabilistic I-states. It is designed for problems in which the state transitions and sensor mapping are linear, and all acts of nature are modeled by multivariate Gaussian densities. Section 11.6.2 covers a general sampling-based planning approach, which is approximate but applies to a broader class of problems. One of these methods, called particle filtering, has become very popular in recent years for mobile robot localization.

I

• • 1. Kalman Filtering

This section covers the most successful and widely used example of a derived I- space that dramatically collapses the history I-space. In the special case in which both f and h are linear functions, and p(θ), p(ψ), and p(x1) are Gaussian, all probabilistic I-states become Gaussian. This means that the probabilistic I-space,

Iprob, does not need to represent every conceivable probability density function. The probabilistic I-state is always trapped in the subspace of Iprob that corresponds only to Gaussians. The subspace is denoted as Igauss. This implies that an I-map,

κmom : prob gauss, can be applied without any loss of information.

I → I

The model is called linear-Gaussian (or LG). Each Gaussian density on Rn is fully specified by its n-dimensional mean vector µ and an n n symmetric covariance matrix, Σ. Therefore, gauss can be considered as a subset of Rm in which m = 2n + (n). For example, if X = R2, then gauss R5, because two independent parameters specify the mean and three independent parameters

2

I ⊂

I

×

specify the covariance matrix (not four, because of symmetry). It was mentioned in Section 11.4.3 that moment-based approximations can be used in general; however,

for an LG model it is important to remember that Igauss is an exact representation of Iprob.

In addition to the fact that the prob collapses nicely, κmom is a sufficient I-map, and convenient expressions exist for incrementally updating the derived I-states entirely in terms of the computed means and covariance. This implies that we can work directly with gauss, without any regard for the original histories or even the general formulas for the probabilistic I-states from Section 11.4.1. The update expressions are given here without the full explanation, which is lengthy but not difficult and can be found in virtually any textbook on stochastic control (e.g., [95, 564]).

I

I

For Kalman filtering, all of the required spaces are Euclidean, but they may have different dimensions. Therefore, let X = Rn, U = Θ = Rm, and Y = Ψ = Rr. Since Kalman filtering relies on linear models, everything can be expressed in terms of matrix transformations. Let Ak, Bk, Ck, Gk, and Hk each denote a matrix with constant real-valued entries and which may or may not be singular. The dimensions of the matrices will be inferred from the equations in which they will appear (the dimensions have to be defined correctly to make the multiplications work out right). The k subscript is used to indicate that a different matrix may

be used in each stage. In many applications, the matrices will be the same in each stage, in which case they can be denoted by A, B, C, G, and H. Since Kalman filtering can handle the more general case, the subscripts are included (even though they slightly complicate the expressions).

In general, the state transition equation, xk+1 = fk(xk, uk, θk), is defined as

xk+1 = Ak_x_k + Bk_u_k + Gkθk, (11.75) in which the matrices Ak, Bk, and Gk are of appropriate dimensions. The notation

fk is used instead of f, because the Kalman filter works even if f is different in every stage.

Example 11.25 (Linear-Gaussian Example) For a simple example of (11.75), suppose X = R3 and U = Θ = R2. A particular instance is

xk+1 =

0 √2 1

1 1 4

− 

 

2 0 1

xk +

1 0

1 1

0 1

uk +

1 1 

0 1

0 −1 

θk. (11.76)

．

The general form of the sensor mapping yk = hk(xk, ψk) is

yk = Ck_x_k + Hkψk, (11.77)

in which the matrices Ck and Hk are of appropriate dimension. Once again, hk is used instead of h because a different sensor mapping can be used in every stage.

So far the linear part of the model has been given. The next step is to specify the Gaussian part. In each stage, both nature actions θk and ψk are modeled with zero-mean Gaussians. Thus, each has an associated covariance matrix, denoted by Σθ and Σψ , respectively. Using the model given so far and starting with an initial Gaussian density over X, all resulting probabilistic I-states will be Gaussian [564].

Every derived I-state in Igauss can be represented by a mean and covariance.

Let µk and Σk denote the mean and covariance of P(xk ηk). The expressions given in the remainder of this section define a derived information transition equation that computes µk+1 and Σk+1, given µk, Σk, uk, and yk+1. The process starts by computing µ1 and Σ1 from the initial conditions.

|

Assume that an initial condition is given that represents a Gaussian density over Rn. Let this be denoted by µ0, and Σ0. The first I-state, which incorporates

the first observation y1, is computed as µ1 = µ0 + L1(y1 − C1µ0) and

Σ1 = (I − L1C1)Σ0, (11.78)

in which I is the identity matrix and

L = Σ CT （C Σ CT + H Σ H ）−1. (11.79)

1

0

1

1

0

1

1

ψ

1

Although the expression for L1 is complicated, note that all matrices have been specified as part of the model. The only unfortunate part is that a matrix inversion is required, which sometimes leads to numerical instability in practice; see [564] or other sources for an alternative formulation that alleviates this problem.

Now that µ1 and Σ1 have been expressed, the base case is completed. The next part is to give the iterative updates from stage k to stage k + 1. Using µk, the mean at the next stage is computed as

µk+1 = Akµk + Bk_u_k + Lk+1(yk+1 − Ck+1(Akµk + Bk_u_k)), (11.80)

in which Lk+1 will be defined shortly. The covariance is computed in two steps; one is based on applying uk, and the other arises from considering yk+1. Thus, after uk is applied, the covariance becomes

Σ′k+1 = AkΣk_A_T + GkΣθ GT . (11.81)

k k

After yk+1 is received, the covariance Σk+1 is computed from Σ′k+1 as

Σk+1 = (I − Lk+1Ck+1)Σ′k+1. (11.82)

The expression for Lk is

L = Σ′ CT （C Σ′ CT + H Σ H ）−1. (11.83)

k

k

k

k

k

k

k

ψ

k

To obtain Lk+1, substitute k + 1 for k in (11.83). Note that to compute µk+1 using (11.80), Σ′k+1 must first be computed because (11.80) depends on Lk+1, which in turn depends on Σ′k+1.

The most common use of the Kalman filter is to provide reliable estimates of the state xk by using µk. It turns out that the optimal expected-cost feedback plan for a cost functional that is a quadratic form can be obtained for LG systems in a closed-from expression; see Section 15.2.2. This model is often called LQG, to reflect the fact that it is linear, quadratic-cost, and Gaussian. The optimal feedback plan can even be expressed directly in terms of µk, without requiring Σk. This indicates that the I-space may be collapsed down to X; however, the corresponding I-map is not sufficient. The covariances are still needed to compute the means, as is evident from (11.80) and (11.83). Thus, an optimal plan can be specified as π : X U, but the derived I-states in gauss need to be represented for the I-map to be sufficient.

→ I

The Kalman filter provides a beautiful solution to the class of linear Gaussian models. It is even successfully applied quite often in practice for problems that do not even satisfy these conditions. This is called the extended Kalman filter. The success may be explained by recalling that the probabilistic I-space may be approx- imated by mean and covariance in a second-order moment-based approximation. In general, such an approximation may be inappropriate, but it is nevertheless widely used in practice.

• • 1. Sampling-Based Approaches

Since probabilistic I-space computations over continuous spaces involve the eval- uation of complicated, possibly high-dimensional integrals, there is strong mo- tivation for using sampling-based approaches. If a problem is nonlinear and/or non-Gaussian, such approaches may provide the only practical way to compute probabilistic I-states. Two approaches are considered here: grid-based sampling and particle filtering. One of the most common applications of the techniques described here is mobile robot localization, which is covered in Section 12.2.

A grid-based approach Perhaps the most straightforward way to numerically compute probabilistic I-states is to approximate probability density functions over a grid and use numerical integration to evaluate the integrals in (11.57) and (11.58). A grid can be used to compute a discrete probability distribution that approx- imates the continuous probability density function. Consider, for example, using the Sukharev grid shown in Figure 5.5a, or a similar grid adapted to the state space. Consider approximating some probability density function p(x) using a finite set, S X. The Voronoi region surrounding each point can be considered as a “bucket” that holds probability mass. A probability is associated with each sample and is defined as the integral of p(x) over the Voronoi region associated

⊂

with the point. In this way, the samples S and their discrete probability distribu- tion, P(s) for all s ∈ S approximate p(x) over X. Let P(sk) denote the probability

distribution over Sk, the set of grid samples at stage k.

In the initial step, P(s) is computed from p(x) by numerically evaluating the integrals of p(x1) over the Voronoi region of each sample. This can alternatively be estimated by drawing random samples from the density p(x1) and then recording the number of samples that fall into each bucket (Voronoi region). Normalizing the counts for the buckets yields a probability distribution, P(s1). Buckets that have little or no points can be eliminated from future computations, depending on the desired accuracy. Let S1 denote the samples for which nonzero probabilities are associated.

Now suppose that P(sk|ηk) has been computed over Sk and the task is to compute P(sk+1|ηk+1) given uk and yk+1. A discrete approximation, P(sk+1|sk, uk),

to p(xk+1 xk, uk) can be computed using a grid and buckets in the manner described above. At this point the densities needed for (11.57) have been approximated by discrete distributions. In this case, (11.38) can be applied over Sk to obtain a

|

grid-based distribution over Sk+1 (again, any buckets that do not contain enough probability mass can be discarded). The resulting distribution is P(sk+1|ηk, uk),

and the next step is to consider yk+1. Once again, a discrete distribution can be computed; in this case, p(xk+1 yk+1) is approximated by P(sk+1 yk+1) by using the grid samples. This enables (11.58) to be replaced by the discrete counterpart (11.39), which is applied to the samples. The resulting distribution, P(sk+1 ηk+1), represents the approximate probabilistic I-state.

|

| |

Particle filtering As mentioned so far, the discrete distributions can be esti- mated by using samples. In fact, it turns out that the Voronoi regions over the samples do not even need to be carefully considered. One can work directly with a collection of samples drawn randomly from the initial probability density, p(x1). The general method is referred to as particle filtering and has yielded good per- formance in applications to experimental mobile robotics. Recall Figure 1.7 and see Section 12.2.3.

Let S X denote a finite collection of samples. A probability distribution is defined over S. The collection of samples, together with its probability distribu- tion, is considered as an approximation of a probability density over X. Since S is

⊂

used to represent probabilistic I-states, let Pk denote the probability distribution over Sk, which is computed at stage k using the history I-state ηk. Thus, at every stage, there is a new sample set, Sk, and probability distribution, Pk.

The general method to compute the probabilistic I-state update proceeds as follows. For some large number, m, of iterations, perform the following:

• • • 1. Select a state xk ∈ Sk according to the distribution Pk.
      2. Generate a new sample, xk+1, for Sk+1 by generating a single sample accord- ing to the density p(xk+1|xk, uk).
      3. Assign the weight, w(xk+1) = p(yk+1|xk+1).

After the m iterations have completed, the weights over Sk+1 are normalized to obtain a valid probability distribution, Pk+1. It turns out that this method provides an approximation that converges to the true probabilistic I-states as m tends to infinity. Other methods exist, which provide faster convergence [536]. One of the main difficulties with using particle filtering is that for some problems it is difficult to ensure that a sufficient concentration of samples exists in the places where they are needed the most. This is a general issue that plagues many sampling-based algorithms, including the motion planning algorithms of Chapter 5.