Particle Filter

A particle filter carries out Bayesian filtering by representing the distribution of a hidden state with weighted random samples, called particles. At every step the particles are moved by the motion model, reweighted by how well they explain the new measurement, and periodically resampled to concentrate computation where the probability is. Because the samples can take any arrangement, the filter can represent skewed and multimodal beliefs and use any model that can be simulated and evaluated, which made it a standard tool for tracking in clutter and for robot localization. The price is computation, which grows quickly with the size of the problem.

Problem

The Kalman filter solves the Bayesian recursion exactly for linear models with Gaussian noise, and the extended Kalman filter approximates it with a single Gaussian otherwise. Both fail when the belief has several peaks, which is common in visual tracking: a person passes behind one of two pillars, or a contour detector responds to several edges near the predicted outline. A single Gaussian must then commit to one peak, risking a confident wrong answer, or cover them all with a mean in an unlikely region between them. The goal is to approximate the full posterior $p(\mathbf{x}_k \mid \mathbf{z}_{1:k})$, whatever its shape.

Inputs and Outputs

Inputs: a way to sample from the transition model $p(\mathbf{x}_k \mid \mathbf{x}_{k-1})$ or another proposal; a likelihood $p(\mathbf{z}_k \mid \mathbf{x}_k)$ that can be evaluated pointwise, up to a constant; samples from the initial distribution $p(\mathbf{x}_0)$; and the number of particles $N$.

Outputs, at each step: a weighted particle set $\{\mathbf{x}_k^{(i)}, w_k^{(i)}\}_{i=1}^{N}$ approximating the posterior, from which estimates such as the weighted mean or the probability of a region are computed.

Intuition

Imagine releasing many hypothetical copies of the tracked object. At every frame, each copy moves according to the motion model plus a random perturbation, so the cloud spreads as uncertainty grows. Each copy is scored by how well it explains the new measurement, and resampling clones high scorers and discards low ones, so the density of copies approximates the probability that the object is there. If the measurement is ambiguous, the cloud keeps several clusters alive until later measurements rule all but one out.

Algorithm

The bootstrap filter, or sampling importance resampling (SIR) filter, introduced by Gordon, Salmond and Smith (1993), is the basic form:

1. Initialize. Draw $\mathbf{x}_0^{(i)} \sim p(\mathbf{x}_0)$ for $i = 1, \ldots, N$ and set $w_0^{(i)} = 1/N$.
2. For each time step $k = 1, 2, \ldots$:
   1. Propagate. Draw $\mathbf{x}_k^{(i)} \sim p(\mathbf{x}_k \mid \mathbf{x}_{k-1}^{(i)})$ by running the motion model with sampled noise.
   2. Weight. Set $w_k^{(i)} \propto w_{k-1}^{(i)} \, p(\mathbf{z}_k \mid \mathbf{x}_k^{(i)})$ and normalize the weights to sum to 1. With no measurement, skip this step.
   3. Estimate the required quantities from the weighted set.
   4. Resample if the weights have become too uneven: draw $N$ particles from the set with probabilities equal to their weights, and reset the weights to $1/N$.

The original bootstrap filter resamples at every step; resampling only when the effective sample size falls below a threshold is now common.

Mathematical Formulation

Weighted-sample approximation

The particle set approximates the posterior by weighted point masses,

$$p(\mathbf{x}_k \mid \mathbf{z}_{1:k}) \approx \sum_{i=1}^{N} w_k^{(i)} \, \delta\!\left(\mathbf{x}_k - \mathbf{x}_k^{(i)}\right), \qquad \sum_{i=1}^{N} w_k^{(i)} = 1,$$

so the posterior expectation of any function $g$ of the state is approximated by $\sum_i w_k^{(i)} g(\mathbf{x}_k^{(i)})$.

Sequential importance sampling

The posterior cannot be sampled directly, so particles are drawn from a proposal distribution $q(\mathbf{x}_k \mid \mathbf{x}_{k-1}^{(i)}, \mathbf{z}_k)$, which may use the newest measurement, and the mismatch is corrected by importance weights:

$$w_k^{(i)} \propto w_{k-1}^{(i)} \, \frac{p(\mathbf{z}_k \mid \mathbf{x}_k^{(i)}) \, p(\mathbf{x}_k^{(i)} \mid \mathbf{x}_{k-1}^{(i)})}{q(\mathbf{x}_k^{(i)} \mid \mathbf{x}_{k-1}^{(i)}, \mathbf{z}_k)} .$$

This is sequential importance sampling (SIS). The bootstrap filter uses the transition model as the proposal, so the transition terms cancel and the update reduces to multiplying by the likelihood. Sampling from the transition model approximates the Bayes filter’s prediction integral, and weighting by the likelihood applies Bayes’ rule.

Degeneracy and effective sample size

Without resampling, SIS degenerates: the variance of the importance weights can only increase over time, and after a few steps all but one particle carry negligible weight (Arulampalam et al., 2002). Degeneracy is measured by the estimated effective sample size

$$\hat{N}_{\text{eff}} = \frac{1}{\sum_{i=1}^{N} \left(w_k^{(i)}\right)^2},$$

which is $N$ when all weights are equal and 1 when one particle holds all the weight. Resampling is triggered when $\hat{N}_{\text{eff}}$ falls below a threshold $N_{\text{T}}$, often $N/2$.

Resampling

Multinomial resampling draws the $N$ indices independently. Systematic resampling draws a single offset $u \sim \mathcal{U}[0, 1)$ and uses the evenly spaced points

$$u_j = \frac{u + j - 1}{N}, \qquad j = 1, \ldots, N,$$

selecting for each $u_j$ the particle $i$ with $c_{i-1} \le u_j < c_i$, where $c_i = \sum_{l \le i} w_k^{(l)}$ is the cumulative weight and $c_0 = 0$. It runs in $O(N)$ time and adds less Monte Carlo variation than multinomial resampling, which is why Arulampalam et al. (2002) prefer it.

Resampling cures degeneracy but causes sample impoverishment: heavily weighted particles are copied many times and the set loses diversity. With very small process noise the particles can collapse onto a single point within a few steps, and because few distinct particle histories survive, estimates of past states from stored paths degenerate (Arulampalam et al., 2002).

Parameters

• Number of particles $N$. Cost grows linearly with $N$. The number required depends mainly on how concentrated the likelihood is relative to the prediction: a precise sensor and a vague prediction leave few particles with real weight.
• Proposal distribution $q$. The most important design choice. The transition prior ignores the current measurement, so it wastes particles when the likelihood is narrow and is sensitive to outliers. The proposal $p(\mathbf{x}_k \mid \mathbf{x}_{k-1}, \mathbf{z}_k)$ minimizes the variance of the weights but can rarely be sampled exactly; Gaussian approximations of it, built by local linearization or the unscented transform, are common compromises (Arulampalam et al., 2002).
• Resampling threshold $N_{\text{T}}$. Resampling too rarely lets the set degenerate; too often accelerates impoverishment.
• Process noise and likelihood sharpness. Process noise also keeps resampled duplicates from staying identical, and is sometimes inflated for that reason. Appearance-based likelihoods are often heuristic, and making them too peaked lets a few particles take all the weight.

Complexity

Each step costs $O(N)$ evaluations of the motion model and likelihood, plus $O(N)$ for normalization and systematic resampling; memory is $O(N n)$ for states of dimension $n$. In vision the likelihood usually dominates, since scoring a particle means comparing image features. Propagation and weighting parallelize well; resampling needs all the weights.

Implementation

An object moves along a line. A range sensor measures its distance to landmark A at position 0, which cannot tell which side of A the object is on; landmark B at $-25$ is detected only within 12 m. The object starts at $-4$ and drifts towards B. The bootstrap filter is compared with an EKF on the same model, started from a guess on the wrong side and gating each measurement with a 99% chi-squared test, as trackers commonly do.

import numpy as np

rng = np.random.default_rng(0)
T, dt = 50, 1.0
q, sig = 0.002, 0.5              # acceleration variance; range noise (m)
LA, LB, reach_B = 0.0, -25.0, 12.0   # landmark A always seen; B within 12 m

F = np.array([[1, dt], [0, 1]])
G = np.array([dt**2 / 2, dt])

def ranges(x):
    """Distances from position(s) x to landmarks A and B."""
    return np.abs(x - LA), np.abs(x - LB)

# Ground truth: starts at -4 and drifts away from A, towards B.
truth = np.array([-4.0, -0.4])
xs, zs = [], []
for k in range(T):
    truth = F @ truth + G * rng.normal(0, np.sqrt(q))
    rA, rB = ranges(truth[0])
    zA = rA + rng.normal(0, sig)
    zB = rB + rng.normal(0, sig) if rB < reach_B else None   # B out of reach
    xs.append(truth[0]); zs.append((zA, zB))

# --- Bootstrap particle filter -------------------------------------------
N = 2000
p = np.column_stack([rng.uniform(-20, 20, N), rng.normal(0, 1, N)])
w = np.full(N, 1.0 / N)

def systematic_resample(w):
    u = (rng.random() + np.arange(len(w))) / len(w)
    return np.minimum(np.searchsorted(np.cumsum(w), u), len(w) - 1)

def pf_step(p, w, z):
    p = p @ F.T + np.outer(rng.normal(0, np.sqrt(q), len(p)), G)  # propagate
    rA, rB = ranges(p[:, 0])
    loglik = -0.5 * ((z[0] - rA) / sig) ** 2       # Gaussian log-likelihood
    if z[1] is not None:
        loglik -= 0.5 * ((z[1] - rB) / sig) ** 2
    w = w * np.exp(loglik - loglik.max()); w /= w.sum()
    if 1.0 / np.sum(w**2) < len(w) / 2:          # effective sample size
        idx = systematic_resample(w)
        p, w = p[idx], np.full(len(w), 1.0 / len(w))
    return p, w

# --- EKF on the same model, started from a guess on the wrong side -------
x, P = np.array([1.0, 0.0]), np.diag([100.0, 1.0])
Q, R1 = q * np.outer(G, G), sig**2

def ekf_step(x, P, z):
    x, P = F @ x, F @ P @ F.T + Q
    for L, zi in ((LA, z[0]), (LB, z[1])):
        if zi is None:
            continue
        H = np.array([np.sign(x[0] - L), 0.0])    # d|x - L| / dx
        S = H @ P @ H + R1
        y = zi - abs(x[0] - L)
        if y**2 / S > 6.63:                       # 99% chi-square gate
            continue                              # rejected as an outlier
        K = P @ H / S
        x = x + K * y
        P = P - np.outer(K, H @ P)
    return x, P

err_ekf, err_pf = [], []
print(" k  true x   EKF x   PF mean   PF P(x<0)   B seen")
for k in range(T):
    p, w = pf_step(p, w, zs[k])
    x, P = ekf_step(x, P, zs[k])
    err_ekf.append(abs(x[0] - xs[k])); err_pf.append(abs(w @ p[:, 0] - xs[k]))
    if k in (0, 5, 10, 16, 17, 18, 30, 49):
        print(f"{k:2d}  {xs[k]:6.2f}  {x[0]:6.2f}  {w @ p[:, 0]:7.2f}"
              f"   {w[p[:, 0] < 0].sum():8.2f}   {'yes' if zs[k][1] is not None else 'no'}")

k_B = next((k for k in range(T) if zs[k][1] is not None), None)  # B first seen
if k_B is None:
    print("B never came into range; the side remains ambiguous")
else:
    print(f"mean |error| from step {k_B} on: EKF {np.mean(err_ekf[k_B:]):.2f} m,"
          f" PF {np.mean(err_pf[k_B:]):.2f} m")

Output:

 k  true x   EKF x   PF mean   PF P(x<0)   B seen
 0   -4.40    4.32    -0.44       0.55   no
 5   -6.24    6.14    -3.07       0.75   no
10   -8.89    8.97    -1.76       0.60   no
16  -12.49   12.39    -2.87       0.62   no
17  -13.14   13.09   -13.16       1.00   yes
18  -13.79   13.86   -13.71       1.00   yes
30  -21.02   21.54   -21.43       1.00   yes
49  -31.05   31.03   -30.92       1.00   yes
mean |error| from step 17 on: EKF 45.18 m, PF 0.19 m

Until B is detected at step 17, the particle filter keeps both mirror-image hypotheses, giving the left side probability 0.55 to 0.75 in the steps shown. Its weighted mean, between about $-3$ and $0$, falls between the two modes, where the object almost certainly is not. The first measurement of B eliminates the wrong mode.

The EKF, linearizing $|x|$ at its initial guess, confidently tracks the mirror image. When B appears, its measurements disagree with the EKF’s prediction by more than 40 standard deviations, the gate rejects them as outliers, and the filter stays wrong. Without the gate, the large innovations drag the EKF across zero within a few steps and it settles near the truth after an overshoot, but an ungated filter accepts clutter as readily.

With other random seeds B comes into range at a different step, or in some runs not within the 50 steps; in 100 seeds tested, whenever B was detected the particle filter settled on the correct side and the gated EKF did not.

Properties and Behavior

• Asymptotic consistency. As $N \to \infty$, the weighted particle set becomes an equivalent representation of the posterior, and particle estimates approach the optimal Bayesian estimate (Arulampalam et al., 2002; Doucet et al., 2001). For finite $N$ the result is random and varies between runs.
• Curse of dimensionality. Snyder et al. (2008) showed that the number of particles needed to avoid weight collapse grows exponentially with the variance of the observation log-likelihood, which typically grows with the state and measurement dimensions; in their simple Gaussian example, a 200-dimensional state required at least $10^{11}$ particles. Plain particle filters therefore suit low-dimensional states; articulated body pose or many jointly tracked objects strain them.
• Loss of track. If no particle lies near the true state, for example after a motion the model considers nearly impossible, reweighting cannot recover it; inflated process noise, or injecting particles near detector outputs, mitigates this.

Limitations

• Cost. When the belief is close to Gaussian, a Kalman-type filter is usually the better choice at a fraction of the computation.
• Point estimates. The mean of a multimodal belief can be meaningless, as in the example, and the mode of a sample set is noisy; extracting a single answer needs care, such as clustering the particles first.

Variants

The auxiliary particle filter of Pitt and Shephard (1999) looks at the new measurement before propagating: it resamples the previous particles according to how well a point prediction of each explains the measurement, then propagates only the selected ones. This helps when process noise is small; when it is large, a point prediction characterizes a particle poorly and performance can degrade (Arulampalam et al., 2002).

The Rao–Blackwellized particle filter samples only part of the state and computes the rest analytically for each sample, typically with a Kalman filter per particle when the remainder is linear-Gaussian; marginalizing analytically reduces the variance of the estimates and the number of particles needed. FastSLAM (Montemerlo et al., 2002) samples the robot’s path and gives each particle a separate small EKF per landmark, since landmark positions are independent given the path.

The regularized particle filter resamples from a continuous kernel density estimate around the particles, so resampled particles are jittered copies rather than duplicates, counteracting impoverishment at the cost of a bandwidth parameter; Markov chain Monte Carlo moves after resampling serve the same purpose (Arulampalam et al., 2002).

Uses in Computer Vision

• Contour and appearance tracking in clutter. CONDENSATION (Isard & Blake, 1998; first presented in 1996) brought particle filtering to visual tracking, propagating a random sample set with learned dynamical models to follow object outlines through dense clutter, where a unimodal Kalman filter cannot hold alternative hypotheses. Later trackers used color-histogram and other appearance likelihoods.
• Monte Carlo localization. Dellaert et al. (1999) applied particle filtering to mobile robot localization in a known map. The method can localize a robot without knowing its starting position, and the authors reported that it was faster, more accurate, and less memory-intensive than earlier grid-based methods.
• Nonlinear and multi-hypothesis tracking. Bearing-only and range-only tracking, where Gordon et al. (1993) found the bootstrap filter greatly superior to the extended Kalman filter in simulation, and tracking with ambiguous data association.