Universal nonlinear filtering using Feynman path integrals II: the continuous-continuous model with additive noise

Abstract

In this paper, the Feynman path integral formulation of the continuous-continuous filtering problem, a fundamental problem of applied science, is investigated for the case when the noise in the signal and measurement model is Gaussian and additive. It is shown that it leads to an independent and self-contained analysis and solution of the problem. A consequence of this analysis is the configuration space Feynman path integral formula for the conditional probability density that manifests the underlying physics of the problem. A corollary of the path integral formula is the Yau algorithm that has been shown to have excellent numerical properties. The Feynman path integral formulation is shown to lead to practical and implementable algorithms. In particular, the solution of the Yau partial differential equation is reduced to one of function computation and integration.

PACS Codes:02.50.Ey, 02.50.Fz, 05.10.Gg, 89.90.+n, 93E10, 93E11

1 Introduction

1.1 Motivation

The fundamental dynamical laws of physics, both classical and quantum mechanical, are described in terms of variables continuous in time. The continuous nature of the dynamical variables has been verified at all length scales probed so far, even though the relevant dynamical variables, and the fundamental laws of physics, are very different in the microscopic and macroscopic realms. In practical situations, one often deals with macroscopic objects whose state is phenomenologically well-described by classical deterministic laws modified by external disturbances that can be modelled as random noise, or Langevin equations. Even when there is no underlying fundamental dynamical law, the Langevin equation provides an effective description of the state variables in many applications. It is therefore natural to consider the problem of the evolution of a state of a signal of interest described by a Langevin equation called the state process.

When the state model noise is Gaussian (or more generally multiplicatively Gaussian) the state process is a Markov process. Since the process is stochastic, the state process is completely characterized by a probability density function. The Fokker-Planck-Kolmogorov foward equation (FPKfe) describes the evolution of this probability density function (or equivalently, the transition probability density function) and is the complete solution of the state evolution problem.

However, in many applications the signal, or state variables, cannot be directly observed. Instead, what is measured is a nonlinearly related stochastic process called the measurement process. The measurement process can often be modelled by yet another continuous stochastic dynamical system called the measurement model. In other words, the observations, or measurements, are discrete-time samples drawn from a different Langevin equation called the measurement process.

The conditional probability density function of the state variables, given the observations, is the complete solution of the filtering problem. This is because it contains all the probabilistic information about the state process that is in the measurements and in the initial condition [1]. This is the Bayesian approach, i.e., the a priori initial data about the signal process contained in the initial probability distribution of the state is incorporated into the solution. Given the conditional probability density, optimality may be defined under various criterion. Usually, the conditional mean, which is the least mean-squares estimate, is studied due to its richness in results and mathematical elegance. The solution of the optimal nonlinear filtering problem is termed universal, if the initial distribution can be arbitrary.

1.2 Fundamental Sochastic Filtering Results

When the state and measurement processes are linear, the linear filtering problem was solved by Kalman and Bucy [2, 3]. The celebrated Kalman filter has been successfully applied to a large number of problems in many different areas.

Nevertheless, the Kalman filter suffers from some major limitations. The Kalman filter is not optimal even for the linear filtering case if the initial distribution is not Gaussian. It may still be optimal for a linear system under certain criteria, such as minimum mean square error, but not a general criterion. In other words, the Kalman filter is not a universal optimal filter, even when the filtering problem is linear. Secondly, the Kalman filter cannot be an optimal solution for the general nonlinear filtering problem since it assumes that the signal and measurement models are linear. The extended Kalman filter (EKF), obtained by applying the Kalman filter to a linearized model, cannot be a reliable solution, in general. Thirdly, even when the EKF estimates the state well in some cases, it gives no reliable indication of the accuracy of the state estimate, i.e., the conditional variance is unreliable. Finally, the Kalman filter assumes that the conditional probability distribution is Gaussian, which is a very restrictive assumption; for instance, it rules out the possibility of a multi-modal conditional probability distribution.

The continuous-continous nonlinear filtering problem (i.e., continuous-time state and measurement stochastic processes) was studied in [4–6] and [7]. This led to a stochastic differential equation, called the Kushner equation, for the conditional probability density in the continuous-continuous filtering problem. It was noted in [8, 9], and [10] that the stochastic differential equation satisfied by the unnormalized conditional probability density, called the Duncan-Mortensen-Zakai (DMZ) equation, is linear and hence considerably simpler than the Kushner equation. The robust DMZ equation, a partial differential equation (PDE) that follows from the DMZ equation via a gauge transformation, was derived in [11] and [12].

A disadvantage of the robust DMZ equation is that the coefficients depend on the measurements. Thus, one does not know the PDE to solve prior to the measurements. As a result, real-time solution is impossible. A fundamental advance was made in tackling the general nonlinear filtering problem by S-T. Yau and Stephen Yau. In [13], it was proved that the robust DMZ equation is equivalent to a partial differential equation that is independent of the measurements, which is referred to as the Yau Equation (YYe) in this paper. Specifically, the measurements only enter as initial condition at each measurement step. Thus, no on-line solution of a PDE is needed; all PDE computations can be done off-line.

However, numerical solution of partial differential equations presents several challenges. A naïve discretization may not be convergent, i.e., the approximation error may not vanish as the grid size is reduced. Alternatively, when the discretization spacing is decreased, it may tend to a different equation, i.e., be inconsistent. Furthermore, the numerical method may be unstable. Finally, since the solution of the YYe is a probability density, it must be positive which may not be guaranteed by the discretization.

A different approach to solving the PDE was taken in [14] and [15]. An explicit expression for the fundamental solution of the YYe as an ordinary integral was derived. It was shown that the formal solution to the YYe may be written down as an ordinary, but somewhat complicated, multi-dimensional integral, with an infinite series as the integrand. In addition, an estimate of the time needed for the solution to converge to the true solution was presented.

1.3 Outline of the Paper

In this paper, the (Euclidean) Feynman path integral (FPI) formulation is employed to tackle the continuous-continuous nonlinear filtering problem. Specifically, phrasing the stochastic filtering problem in a language common in physics, the solution of the stochastic filtering problem is presented. In particular, no other result in filtering theory (such as the DMZ equation, the robust DMZ equation, etc.) is used. The path integral formulation leads to a path integral formula for the transition probability density for the general additive noise case. A corollary of the FPI formalism is the path integral formula for the fundamental solution of the YYe and the Yau algorithm – a fundamental result of nonlinear filtering theory. It is noted that this paper provides a detailed derivation of results that were used in [16].

The following point needs to be emphasized to readers familiar with the discussion of standard filtering theory – the FPI is different from the Feynman-Kǎc path integral. In filtering theory literature, it is the Feynman-Kǎc formalism that is often used. The Feynman-Kǎc formulation is a rigorous formulation and has led to several rigorous results in filtering theory. However, in spite of considerable effort it has not been proven to be directly useful in the development of reliable practical algorithms with desirable numerical properties. It also obscures the physics of the problem.

In contrast, it is shown that the FPI leads to formulas that are eminently suitable for numerical implementation. It also provides a simple and clear physical picture. Many path integral manipulations have no counterpart in the Feynman-Kǎc approach. Finally, the theoretical insights provided by the FPI are highly valuable, as evidenced by numerous examples in modern theoretical physics (see, for instance, [17]), and shall be employed in subsequent papers.

The outline of this paper is as follows. In the following section, the filtering problem is reformulated in a language common in physics. In Section 3, the path integral formula for the transition probability density is derived for the general additive noise case. The Yau algorithm is then derived in the following section. In Sections 5 and 6 some conceptual remarks and numerical examples are presented. The conclusions are presented in Section 7. In Appendix 1, aspects of continuous-continuous filtering are reviewed.

For more details on the path integral methods, see any modern text on quantum field theory, such as [17], and especially [18] which discusses application of FPI to the study of stochastic processes.

2 The continuous filtering problem: a physical reformulation

In this section, the filtering problem is stated in a language commonly used in theoretical physics.

2.1 Langevin Equation

Consider an ensemble of systems with state variables described by the Langevin equation:

$$\begin{array}{ccc}\dot{x}(t)=f(x(t),t)+e(x(t),t)\nu (t),& x(0)=x_0,& Q(t)\in {ℝ}^{p\times p}\end{array}.$$

(1)

Here, x(t) ∈ ℝⁿ, the drift f(x(t), t) ∈ ℝⁿ, the diffusion vielbein e(x(t), t) ∈ ℝ^n×p, and ν(t) is δ-correlated with covariance matrix Q(t) ∈ ℝ^p×p. When the diffusion vielbein is independent of the state x(t), the noise is termed additive. It is the additive noise that is studied here, since it enables the use of functional methods common in quantum field theory.

Due to the random noise, each system leads to a different vector x(t) that depends on time. Although only one realization of the stochastic process is ever observed, it is meaningful to speak about an ensemble average. For fixed times t = t _i , i = 1, 2, ..., r, the probability density of finding the random vector x(t) in the (n-dimensional) interval x _i ≤ x(t _i ) ≤ x _i + dx _i (1 ≤ i ≤ r) is given by

$$W_r(t_r,x_r;\cdots ;t_1,x_1)=〈{\displaystyle \prod _{i=1}^r{\delta }^n(x(t_i)-x_i)}〉,$$

(2)

where x _i is an n-dimensional column vector and ⟨·⟩ denotes averaging with respect to the signal model noise ν(t). The complete information on the random vector x(t) is contained in the infinite hierarchy of such probability densities. The quantity of interest here is the conditional probability density

$$\begin{array}{llll}p(t_r,x_r|t_{r-1},x_{r-1},\mathrm{...};t_1,x_1)\hfill & =\hfill & 〈{\delta }^n(x(t_r))〉{|}_{x(t_{r-1})=x_{r-1,\mathrm{...},x(t_1)=x_1}},\hfill & x(t_r)\equiv x_r\hfill \\ =\hfill & \frac{W_r(t_r,x_r;\mathrm{...};t_1,x_1)}{{\displaystyle \int W_n(t_r,x_r;\mathrm{...};t_1,x_1)\left\{d^n{x}_r\right\}}}.\hfill \end{array}$$

(3)

Now the process described by the Langevin equation with δ-correlated Langevin force is a Markov process, i.e., the conditional probability density depends only on the value at the immediate previous time:

p(t _n , x _n |t_n-1, x_n-1; ...; t₁, x₁) = p(t _n , x _n |t_n-1, x_n-1).

It can be shown that the transition probability satisfies the Fokker-Planck-Kolmogorov forward equation (FPKfe) (see, for instance, [19])

$$\begin{array}{lll}\frac{\partial p}{\partial t}(t,x)\hfill & =\hfill & -{\displaystyle \sum _{i=1}^n\frac{\partial }{\partial x_i}[f_i(x(t),t)p(t,x)]+\frac{1}{2}}{\displaystyle \sum _{i,j=1}^n{\displaystyle \sum _{a=1}^p\frac{{\partial }^2}{\partial x_i\partial x_j}[{(e(x(t),t)Q(t))}_{ia}{e}_{aj}^T}(x(t),t)p(t,x))]},\hfill \\ \equiv \hfill & {\mathcal{L}}_p(t,x).\hfill \end{array}$$

(5)

Finally, the Gaussian noise process ν(t) can be represented by the following path integral measure

$$\begin{array}{cc}[d\rho (\nu (t))]=[\mathcal{D}\nu (t)]\exp \left[-\frac{1}{2}{\displaystyle \sum _{a,b=1}^p{\displaystyle \int {\nu }_a(t){[Q^{-1}(t)]}_{ab}{\nu }_b(t)dt}}\right],& \nu \in {ℝ}^p.\end{array}$$

(6)

where ν(t) is a real vector for each t. This leads to a configuration space FPI formula for the fundamental solution of the FPKfe (see, for instance, [18]). The path integral formula for the fundamental solution for the FPkfe is applied to the continuous-discrete filtering problem with additive (state model) noise in [20, 21].

2.2 The Continuous-Continuous Filtering Problem

Similarly, the continuous-continuous model can be written as follows:

$$\{\begin{array}{cccc}\dot{x}(t)& =f(x(t),t)+e(x(t),t)\nu (t),& x(0)=x_0,& Q(t)\in {ℝ}^{p\times p}\\ \dot{y}(t)& =h(x(t),t)+n(x(t),t)\mu (t),& y(0)=y_0,& R(t)\in {ℝ}^{q\times q}\end{array}$$

(7)

Here, y(t) ∈ ℝ^m, the measurement model drift h(x(t), t) ∈ ℝ^m, the diffusion vielbein n(x(t), t) ∈ ℝ^m×q, and μ(t) is δ-correlated with covariance matrix W(t) ∈ ℝ^q×q.

Thus, in continuous-continuous filtering, the continuous-time measurement stochastic process needs to be incorporated as well. Consider another ensemble of systems with state variables whose time evolution is governed by the measurement process. The measurement noise means that each system in the ensemble leads to a different time-dependent vector y(t). Thus, even though only one realization of the measurement stochastic process is observed, it is still meaningful to talk about an ensemble average of the measurement process (in addition to one over the state process). Thus, the quantity of interest in continuous-continuous filtering is the conditional probability density

$$P(t_2,x_2;y_2|t_1,x_1;y_1)={〈〈{\delta }^n(x(t_2)-x_2){\delta }^m(y(t_2)-y_2)〉〉}_{\mu }{|}_{x(t_1)=x_1,y(t_1)=y_1},$$

(8)

where ⟨·⟩ _μ denotes averaging with respect to the measurement noise μ(t). A crucial difference between the state and measurement stochastic process is that, unlike the state, the measurement samples are known.

Note that the conditional transition probability density is the complete solution to the continuous-continuous filtering problem, since if the initial distribution is u(t_i-1, x'|Y_i-1), where Y_i-1 is the set of all measurements prior to t_i-1, then the evolved conditional probability distribution is

u(t, x|Y _i ) = ∫ P(t, x; y _i |t_i-1, x'; y_i-1) u (t_i-1, x'|Y_i-1) {dⁿx'}.

In the following sections, the path integral formulas for P(t₂, x₂; y₂|t₁, x₁; y₁) are derived. In Section 4 it shall be shown that it leads to the Yau algorithm. It shall be shown that the YYe plays the same role here that the FPKfe does in continuous-discrete filtering.

3 Path integral formula for the conditional transition probability density

The path integral formula for the conditional transition probability density shall now be derived using functional methods. Note that implicit in the use of these formal functional methods is the use of the Feynman convention, or symmetric discretization for the drift.

As noted in Section 2, the transition probability density is computed by averaging over the signal and measurement ensembles, i.e.,

$$\begin{array}{lll}P(t,x;y|t_0,x_0;y_0)\hfill & =\hfill & {\displaystyle \int [d\rho (\nu (t,t_0))][d\rho (\mu (t,t_0))]}\hfill \\ \times [\mathcal{D}\dot{x}(t)]J{\delta }^n[\dot{x}(t)-f(x(t),t)-e(t)\nu (t)]\hfill \\ \times [\mathcal{D}\dot{y}(t)]J_y{\delta }^m[\dot{y}(t)-h(x(t),t)-n(t)\mu (t)]\hfill \\ \times {\delta }^n[x(t)-x]{|}_{x(t_0)=x_0}{\delta }^m[y(t)-y]{|}_{y(t_0)=y_0}.\hfill \end{array}$$

(10)

From the assumptions of the signal and measurement noise processes, it is evident that

$$[d\rho (\nu (t,t_0))]=[\mathcal{D}\nu (t,t_0)]\exp \left(-\frac{1}{2}{\displaystyle \sum _{a,b=1}^p{\displaystyle {\int }_{t_0}^t{\nu }_a(t){[Q^{-1}(t)]}_{ab}{\nu }_b(t)dt}}\right),$$

(11)

$$[d\rho (\mu (t,t_0))]=[\mathcal{D}\mu (t,t_0)]\exp \left(-\frac{1}{2}{\displaystyle \sum _{c,d=1}^p{\displaystyle {\int }_{t_0}^t{\mu }_c(t){[W^{-1}(t)]}_{cd}{\mu }_d(t)dt}}\right).$$

(12)

The Jacobian J follows from the functional derivative of the Langevin equation:

$$\begin{array}{lll}\frac{\delta }{\delta x_j(t^{\prime })}\left[{\dot{x}}_i(t)-f_i(x(t),t)-{\displaystyle \sum _{a=1}^p{e}_{ia}(t){\nu }_a(t)}\right]\hfill & =\hfill & \left[{\delta }_{ij}\frac{d}{dt}-\frac{\partial f_i}{\partial x_j}(x(t),t)\right]\delta (t-t^{\prime }),\hfill \\ =\hfill & -\frac{d}{d{t}^{\prime }}\left[{\delta }_{ij}\delta (t-t^{\prime })-\theta (t-t^{\prime })\frac{\partial f_i}{\partial x_j}(x(t),t)\right].\hfill \end{array}$$

(13)

Hence,

$$\begin{array}{c}J=\det \left(-\frac{d}{d{t}^{\prime }}\right)\det \left({\delta }_{ij}\delta (t-t^{\prime })-\theta (t-t^{\prime })\frac{\partial f_i}{\partial x_j}(x(t),t)\right),\\ =\mathcal{N}\det \left({\delta }_{ij}\delta (t-t^{\prime })-\theta (t-t^{\prime })\frac{\partial f_i}{\partial x_j}(x(t),t))\right),\end{array}$$

(14)

where $\mathcal{N}$ is an irrelevant constant, or,

$$\begin{array}{c}\ln J=\ln \det \left[{\delta }_{ij}\delta (t-t^{\prime })-\theta (t-t^{\prime })\frac{\partial f_i}{\partial x_j}(x(t),t)\right],\\ =-\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \int \frac{\partial f_i}{\partial x_i}(x(t),t)dt}}.\end{array}$$

(15)

The Jacobian J _y is trivial (as the measurement model drift is y-independent) and can be absorbed into the measure.

It is noteworthy that J is not trivial. In quantum field theory, nontrivial Jacobians usually imply that there is an anomaly, as in the case of chiral anomalies in gauge theories. However, there is no reason for an anomaly here; after all, this is not even a quantum field theoretical system. The puzzle is resolved by noting that path integral anomalies in quantum field theory arise from the "multiplicative" part in the change of variables (i.e., ψ (x) → U(x) ψ (x)). In contrast, the nontrivial Jacobian term here arises from the additive term; the multiplicative term does not contribute to the Jacobian, in accordance with expectations.

Thus, so far,

$$\begin{array}{lll}P(t,x;y|t_0,x_0;y_0)\hfill & =\hfill & {\displaystyle {\int }_{y(t_0)=y_0}^{y(t)=y}{\displaystyle {\int }_{x(t_0)=x_0}^{x(t)=x}[\mathcal{D}x(t)][\mathcal{D}y(t)][d\rho (\nu (t,t_0))][d\rho (\mu (t,t_0))]}}\hfill \\ \times {\delta }^n\left({\dot{x}}_i(t)-f_i(x(t),t)-{\displaystyle \sum _{a=1}^p{e}_{ia}(t){\nu }_a(t)}\right)\exp \left(-\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle {\int }_{t_0}^t\frac{\partial f_i}{\partial x_i}(x(t),t)dt}}\right)\hfill \\ \times {\delta }^m\left({\dot{y}}_k(t)-h_k(x(t),t)-{\displaystyle \sum _{c=1}^q{n}_{kc}(t){\mu }_c(t)}\right).\hfill \end{array}$$

(16)

Using the Fourier integral version of the delta function

$$\begin{array}{lll}P(t,x;y|t_0,x_0;y_0)\hfill & =\hfill & {\displaystyle {\int }_{y(t_0)=y_0}^{y(t)=y}{\displaystyle {\int }_{x(t_0)=x_0}^{x(t)=x}[\mathcal{D}x(t)][\mathcal{D}y(t)][\mathcal{D}\nu (t,t_0)][\mathcal{D}\mu (t,t_0)]}}\hfill \\ [D\lambda (t,t_0)]\exp \left(i{\displaystyle \sum _{i=1}^n{\displaystyle {\int }_{t_0}^t{\lambda }_i(t)\left[{\dot{x}}_i(t)-f_i(x(t),t)-{\displaystyle \sum _{a=1}^p{e}_{ia}(t){\nu }_a(t)}\right]dt}}\right)\hfill \\ [D\kappa (t,t_0)]\exp \left(i{\displaystyle \sum _{k=1}^n{\displaystyle {\int }_{t_0}^t{\kappa }_k(t)\left[{\dot{y}}_k(t)-h_k(x(t),t)-{\displaystyle \sum _{c=1}^q{n}_{kc}(t){\mu }_c(t)}\right]dt}}\right)\hfill \\ \times \exp \left(-\frac{1}{2}{\displaystyle \sum _{a,b=1}^p{\displaystyle {\int }_{t_0}^t{\nu }_a(t){(Q^{-1}(t))}_{ab}{\nu }_b(t)dt}}\right)\times \exp \left(-\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle {\int }_{t_0}^t\frac{\partial f_i}{\partial x_i}(x(t),t)dt}}\right)\hfill \\ \times \exp \left(-\frac{1}{2}{\displaystyle \sum _{c,d=1}^q{\displaystyle {\int }_{t_0}^t{\mu }_c(t){(W^{-1}(t))}_{cd}{\mu }_d(t)dt}}\right).\hfill \end{array}$$

(17)

Integrating over ν(t, t₀), and μ(t, t₀) leads to

$$\begin{array}{c}P(t,x;y|t_0,x_0;y_0)={\displaystyle {\int }_{y(t_0)=y_0}^{y(t)=y}{\displaystyle {\int }_{x(t_0)=x_0}^{x(t)=x}[\mathcal{D}x(t)][\mathcal{D}y(t)][\mathcal{D}\lambda (t,t_0)][\mathcal{D}\kappa (t,t_0)]}}\exp \left(-\frac{1}{2}{\displaystyle \sum _{i=1}^n{\displaystyle \int \frac{\partial f_i}{\partial x_i}(x(t),t)dt}}\right)\\ \times \exp \left(-{\displaystyle {\int }_{t_0}^t\left[{\displaystyle \sum _{i,j=1}^n{\displaystyle \sum _{a,b=1}^p{\lambda }_i(t)e_{ia}(t)Q_{ab}(t)e_{bj}^T(t){\lambda }_j(t)}-i{\displaystyle \sum _{i=1}^n{\lambda }_i(t)({\dot{x}}_i(t)-f_i(x(t),t))}}\right]dt}\right)\\ \times \exp \left(-{\displaystyle {\int }_{t_0}^t\left[{\displaystyle \sum _{k,l=1}^m{\displaystyle \sum _{c,d=1}^q{\kappa }_k(t)n_{kc}(t)W_{cd}(t)n_{dl}^T(t){\kappa }_l(t)}-i{\displaystyle \sum _{k=1}^m{\kappa }_k(t)({\dot{y}}_k(t)-h_k(x(t),t))}}\right]dt}\right)\end{array}$$

(18)

Integrating over λ(t, t₀) and κ(t, t₀), it is clear that

$$P(t,x;y|t_0,x_0;y_0)={\displaystyle {\int }_{y(t_0)=y_0}^{y(t)=y}{\displaystyle {\int }_{x(t_0)=x_0}^{x(t)=y}[\mathcal{D}x(t)][\mathcal{D}y(t)]\exp [-S(t_0,t)]}},$$

(19)

where the action is given by

$$\begin{array}{l}S(t_0,t)=\frac{1}{2}{\displaystyle \sum _{i,j=1}^n{\displaystyle \sum _{a,b=1}^p{\displaystyle {\int }_{t_0}^t[{\dot{x}}_i(t)-f_i(x(t),t)]{(e_{ia}(t)Q_{ab}(t)e_{bj}^T(t))}^{-1}[{\dot{x}}_j(t)-f_j(x(t),t)]dt}}}\hfill \\ +\frac{1}{2}{\displaystyle {\int }_{t_0}^t\left[{\displaystyle \sum _{k,l=1}^n{\displaystyle \sum _{c,d=1}^q[{\dot{y}}_k(t)-h_k(x(t),t)}]{(n_{kc}(t)W_{cd}(t)n_{dl}^T(t))}^{-1}[{\dot{y}}_l(t)-h_l(x(t),t)]+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x(t),t)}}\right]dt.}\hfill \end{array}$$

(20)

4 Derivation of the Yau algorithm

Observe that the path integral formula derived in the previous section is over both the state and measurement variables. In this section, we shall show that in some cases it is possible to reduce the result to a path integral over the state variables only. This has the advantage of being pre-computable since it is independent of measurements. It shall be shown to lead to the Yau algorithm.

4.1 Sampled Continuous Measurements

Suppose measurements are available at time t_i-1 and t _i , and that the conditional transition probability density at time t _i is desired. Further, assume that there are no measurements available between t_i-1 and t _i . The general path integral formula (Equation 19) cannot be simplified unless some additional assumptions are made.

First, consider the case where e(t)Q(t)e^T(t) and n(t)W(t)n^T(t) are ħ _ν I_n×nand ħ _μ I_m×mand h(x(t)) is not explicitly time dependent. Then, the contribution to the action due to the measurement process is

$$-\frac{1}{2{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{\displaystyle {\int }_{t_i-1}^{t_i}dt{[{\dot{y}}_k-h_k(x(t))]}^2}}=-\frac{1}{2{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{\displaystyle {\int }_{t_{i-1}}^{t_i}dt[{\dot{y}}_k^2(t)+h_k^2(x(t))-2h_k(x(t)){\dot{y}}_i(t)]}}.$$

(21)

The quantity of interest is the state and we seek to integrate over the measurement variables.

Now the first term is independent of the state variables. The second term can be added to the action term that is independent of y(t). It remains to investigate the contribution of the third term:

$$\frac{1}{{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{\displaystyle {\int }_{t_{i-1}}^{t_i}{h}_k(x(t)){\dot{y}}_k(t)dt}}.$$

(22)

There are two issues in this evaluation. Firstly, this can be evaluated via the usual integration by parts, but it is important to note that it is valid only for symmetric discretization. Secondly, since the measurements are sampled, there are two possible interpretations when t _i - t_i-1 = ϵ, where ϵ is an infinitesimal:

$$\frac{1}{{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{\displaystyle {\int }_{t_{i-1}}^{t_i}{h}_k(x(t)){\dot{y}}_k(t)dt}}=\{\begin{array}{l}\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{\displaystyle {\int }_{t_{i-1}}^{t_i}{h}_k(x(t)){\dot{y}}_k(t)dt}},\hfill \\ \frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{\displaystyle {\int }_{t_{i-1}-ϵ}^{t_i-ϵ}{h}_k(x(t)){\dot{y}}_k(t)dt}}.\hfill \end{array}$$

(23)

This leads to two possibilities:

$$\frac{1}{{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{\displaystyle {\int }_{t_{i-1}}^{t_i}{h}_k(x(t)){\dot{y}}_k(t)dt}}=\{\begin{array}{l}\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_i)+x(t_{i-1})\right]\right)\left[y_k(t_i)-y_k(t_{i-1})\right]},\hfill \\ \frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_{i-1})+x(t_{i-2})\right]\right)\left[y_k(t_{i-1})-y_k(t_{i-2})\right]}.\hfill \end{array}$$

(24)

Finally, the residual Gaussian path integral over y(t) can be ignored as it is independent of the state. Therefore, the path integral formula simplifies to

$$\begin{array}{lll}P(t_i,x_i;y_i|t_{i-1},x_{i-1};y_{i-1})\hfill & =\hfill & \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\hfill \\ \times \{\begin{array}{l}\exp \left(\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_i)+x(t_{i-1})\right]\right)\left[y_k(t_i)-y_k(t_{i-1})\right]}\right),\hfill \\ \exp \left(\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_{i-1})+x(t_{i-2})\right]\right)\left[y_k(t_{i-1})-y_k(t_{i-2})\right]}\right),\hfill \end{array}\hfill \end{array}$$

(25)

where

$$\tilde{P}(t_i,x_i|t_{i-1},x_{i-1})={\displaystyle {\int }_{x(t_{i-1})=x_{i-1}}^{x(t_i)=x_i}[\mathcal{D}x(t)]\exp (-S(t_{i-1},t_i)),}$$

(26)

and

$$S(t_{i-1},t_i)=\frac{1}{2}{\displaystyle {\int }_{t_{i-1}}^{t_i}dt}\left[\frac{1}{{\hslash }_{\nu }}{\displaystyle \sum _{i=1}^n{\left[{\dot{x}}_i(t)-f_i(x(t)\right]}^2}+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x(t))+\frac{1}{{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{h}_k^2(x(t))}}\right].$$

(27)

Secondly, if the conditions above are relaxed to allowing explicit time dependence of the drift term in the state model, i.e., f(x(t), t), and and C = (n(t)W(t)n^T(t)), then the path integral formula becomes

$$\begin{array}{lll}P(t_i,x_i;y_i|t_{i-1},x_{i-1};y_{i-1})\hfill & =\hfill & \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\hfill \\ \times \{\begin{array}{l}\exp (\left(\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k,l=1}^m{h}_k\left(\frac{1}{2}\left[x(t_i)+x(t_{i-1})\right]\right){(C^{-1})}_{kl}\left[y_l(t_i)-y_l(t_{i-1})\right]}\right),\hfill \\ \exp (\left(\frac{1}{{\hslash }_{\mu }}{\displaystyle {\sum }_{k,l=1}^m{h}_k\left(\frac{1}{2}\left[x(t_{i-1})+x(t_{i-2})\right]\right){(C^{-1})}_{kl}\left[y_l(t_{i-1})-y_l(t_{i-2})\right]}\right),\hfill \end{array}\hfill \end{array}$$

(28)

where

$$\tilde{P}(t_i,x_i|t_{i-1},x_{i-1})={\displaystyle {\int }_{x(t_{i-1})=x_{i-1}}^{x(t_i)=x_i}[Dx(t)]\exp (-S(t_{i-1},t_i)),}$$

(29)

and the action is given by

$$S(t_{i-1},t_i)=\frac{1}{2}{\displaystyle {\int }_{t_{i-1}}^{t_i}dt}\left[\frac{1}{{\hslash }_{\nu }}{\displaystyle \sum _{i=1}^n{\left[{\dot{x}}_i(t)-f_i(x(t),t)\right]}^2}+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x(t),t)+{\displaystyle \sum _{k,l=1}^m{h}_k(x(t)){(C^{-1})}_{kl}{h}_l(x(t))}}\right].$$

(30)

Finally, consider the case when e(t)Q(t)e(t) is time-independent and invertible, but otherwise arbitrary, and C = (n(t)W(t)n^T(t)). Note that C is a constant, symmetric matrix and assumed to be invertible. The path integral formula is given by Equations 28 and 29 and with the action S(t_i-1, t _i ) is given by

$$\begin{array}{c}S(t_{i-1},t_i)=\frac{1}{2}{\displaystyle {\int }_{t_{i-1}}^{t_i}dt}{\displaystyle \sum _{i,j=1}^n{\displaystyle \sum _{a,b=1}^p\left[{\dot{x}}_i(t)-f_i(x(t),t)\right]{\left[e_{ia}(t)Q_{ab}(t)e_{bj}^T(t)\right]}^{-1}\left[{\dot{x}}_j(t)-f_j(x(t),t)\right]}}\\ +\frac{1}{2}{\displaystyle {\int }_{t_{i-1}}^{t_i}dt}\left[{\displaystyle \sum _{k,l=1}^m{h}_k(x(t)){(C^{-1})}_{kl}{h}_l(x(t))}+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x(t),t)}\right].\end{array}$$

(31)

4.2 The Yau Algorithm

In Section 2 it was noted that if v_i-1 (t_i-1, x_i-1) is the conditional probability density at time t_i-1, then the conditional probability density at time t _i is given by

v _i (t _i , x _i ) = ∫ P (t _i , x _i ; y _i |t_i-1, x_i-1; y_i-1) v_i-1 (t_i-1, x_i-1) {dⁿx_i-1}.

For simplicity, let us first consider the case e(t)Q(t)e^T(t) and n(t)W(t)n^T(t) are ħ _ν I_n×nand ħ _μ I_m×m, and h(x(t)) is not explicitly time dependent. Then

$$\begin{array}{l}{v}_i(t_i,x_i)=\\ \{\begin{array}{l}{\displaystyle \int \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\exp \left(-\frac{1}{{\hslash }_{\nu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_i)+x(t_{i-1})\right]\right)\left[y_k(t_i)-y_k(t_{i-1})\right]}\right)v_{i-1}(t_{i-1},x_{i-1})\{d^n{x}_{i-1}\},}\hfill \\ {\displaystyle \int \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\exp \left(-\frac{1}{{\hslash }_{\nu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(\frac{1}{2}\left[x(t_{i-1})+x(t_{i-2})\right]\right)\left[y_k(t_{i-1})-y_k(t_{i-2})\right]}\right)v_{i-1}(t_{i-1},x_{i-1})\{d^n{x}_{i-1}\}.}\hfill \end{array}\end{array}$$

(33)

When t _i - t_i-1 is small

$$h_k\left(\frac{1}{2}[x(t_i)+x(t_{i-1})]\right)\approx h_k(x(t_{i-1}))+O(\sqrt{ϵ}),$$

(34)

and Equation 33 becomes

$$v_i=(t_i,x_i)=\{\begin{array}{l}{\displaystyle \int \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\exp \left(-\frac{1}{{\hslash }_{\nu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(x(t_{i-1})\right)\left[y_k(t_i)-y_k(t_{i-1})\right]}\right)v_{i-1}(t_{i-1},x_{i-1})\{d^n{x}_{i-1}\},}\hfill \\ {\displaystyle \int \tilde{P}(t_i,x_i|t_{i-1},x_{i-1})\exp \left(-\frac{1}{{\hslash }_{\nu }}{\displaystyle {\sum }_{k=1}^m{h}_k\left(x(t_{i-1})\right)\left[y_k(t_{i-1})-y_k(t_{i-2})\right]}\right)v_{i-1}(t_{i-1},x_{i-1})\{d^n{x}_{i-1}\}.}\hfill \end{array}$$

(35)

Following the method originally used by Feynman, the partial differential equation satisfied by $\tilde{P}$ (t, x|t₀, x₀) may be derived (see, for instance, [20]). In particular, $\tilde{P}$ (t _i , x _i |t_i-1, x_i-1) is the fundamental solution of the Yau Equation(YYe):

$$\{\begin{array}{l}\frac{\partial \stackrel{\sim }{P}}{\partial t}(t,x|t_0,x_0)=\frac{{\hslash }_{\nu }}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2\stackrel{\sim }{P}}{\partial x_i^2}(t,x|t_0,x_0)-{\displaystyle \sum _{i=1}^n\frac{\partial }{\partial x_i}\left[f_i(x)\stackrel{\sim }{P}(t,x|t_0,x_0)\right]}}-\frac{1}{2{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{h}_k^2(x)\stackrel{\sim }{P}(t,x|t_0,x_0)},\hfill \\ \stackrel{\sim }{P}(t_0,x|t_0,x_0)={\delta }^n(x-x_0).\hfill \end{array}$$

(36)

This implies that v _i (t _i , x) is the solution at t _i of

$$\{\begin{array}{l}\frac{\partial v_i}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{l=1}^n\frac{{\partial }^2{v}_i}{\partial x_l^2}(t,x)-{\displaystyle \sum _{l=1}^n{f}_l(x)\frac{\partial v_i}{\partial x_l}(t,x)}}-\left({\displaystyle \sum _{l=1}^n\frac{\partial f_l}{\partial x_l}(x)+\frac{1}{2}{\displaystyle \sum _{l=1}^m{h}_l^2(x)}}\right)v_i(t,x),\hfill \\ v_i(t_{i-1},x)=v_{i-1}(t_{i-1},x)\{\begin{array}{l}\exp \left({\displaystyle {\sum }_{j=1}^m{h}_j(x)[y_j(t_i)-y_j(t_{i-1})]}\right),\hfill \\ \exp \left({\displaystyle {\sum }_{j=1}^m{h}_j(x)[y_j(t_{i-1})-y_j(t_{i-2})]}\right).\hfill \end{array}\hfill \end{array}$$

(37)

This is precisely the Yau algorithm.

Likewise, it is straightforward to see that for the general case studied in Section 4.1, the Yau algorithm is extended to this case as follows: v _i (t _i , x _i ) is the solution at t _i of the PDE

$$\{\begin{array}{c}\frac{\partial v_i}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{l,j=1}^n{\displaystyle \sum _{a,b=1}^p\left(e_{la}(t)Q_{ab}(t)e_{bj}^T(t)\right)\frac{{\partial }^2{v}_i}{\partial x_l\partial x_j}(t,x)}}\\ -{\displaystyle \sum _{l=1}^n\frac{\partial }{\partial x_l}[f_l(x,t)v_i(t,x)]-\frac{1}{2}{\displaystyle \sum _{k,l=1}^m{h}_k(x){(C^{-1})}_{kl}{h}_l(x(t))v_i(t,x)}},\\ v_i(t_{i-1},x)=v_{i-1}(t_{i-1},x)\{\begin{array}{l}\exp \left({\displaystyle {\sum }_{j,k=1}^m{h}_j(x){(C^{-1})}_{jk}[y_k(t_i)-y_k(t_{i-1})]}\right),\hfill \\ \exp \left({\displaystyle {\sum }_{j,k=1}^m{h}_j(x){(C^{-1})}_{jk}[y_k(t_{i-1})-y_k(t_{i-2})]}\right).\hfill \end{array}\end{array}$$

(38)

This is a straightforward generalization of the YYe to the state model with explicit time dependence.

5 Some remarks

Following are some remarks on the FPI solution of the filtering problem:

• Note that the FPI formulation has given a complete and self-contained solution of the continuous-continuous filtering problem. For instance, the DMZ equation or its variants were not used as an input. On the contrary, the FPI formula naturally led us to the YYe and the Yau algorithm.

• Note also that the DMZ equation (and variants) cannot be solved reliably in real-time as the various approximations assume that drift is bounded and require solution of a stochastic PDE that depend on measurements. In contrast, the general FPI formula presented in Section 3 can potentially be used for an efficient and reliable real-time solution. This is because the measurement time interval is usually small so that the simplest approximation of the path integral (termed the Dirac-Feynman approximation, see Section 6) is adequate. In contrast, in quantum mechanics and quantum field theory one is interested in the large time case.

• Unlike the Yau algorithm, the PI formula is valid even for the general time-dependent case for large measurement time interval. In other words, one can compute the conditional transition probablity density using the conventional methods (see, for instance, [22]).

• The YYe can be viewed as a local expression of the path integral formula. That is, a path integral is a global object, while the PDE is a local one.

• In this paper, the signal noise and measurement noise are assumed to be additive. This is a stronger condition than the orthogonalilty of the diffusion vielbein assumed in [13]. The solution for the general case has been presented in [16].

• It is noted that other algorithms can also be solved using the FPI formulas with obvious changes. Usually, they require the solution of the FPKfe, which corresponds to the case h(x) = 0. Note that the FPKfe arises naturally in the solution of the continuous-discrete filtering problem [20]. Of course, as noted in Appendix 1, it would be unnecessary since the Yau algorithm has the best numerical properties. What is interesting to note is that the path integral formula naturally leads to the algorithm with the best numerical properties.

• The Yau algorithm also has the form of the "prediction" and "correction" part, as in continuous-discrete filtering [20]. Specifically, the prediction part is the solution of the YYe, whereas the correction part is the multiplicative factor in the initial condition. However, it is crucial to note that the the prediction part contains the measurement model drift. In contrast, the measurement model plays no role in the prediction part in continuous-discrete filtering.

6 Examples

6.1 The Dirac-Feynman Approximation

From the discussion in the previous sections, it is evident that computing the path integral requires computing the Lagrangian, L, defined by S = ∫ dt L(x, $\dot{x}$, t). The simplest (and crudest) approximation (for the case e(t)Q(t)e^T(t) = ħ _ν I_n×nand n(t) W (t)n^T(t) = ħ _μ I_m×m) is to use the following approximation that is valid for infinitesimal time step:

$$\begin{array}{l}\tilde{P}(t^{″},x^{″}|t^{\prime },x^{\prime })=\hfill \\ \exp \left(-\frac{t^{″}-t^{\prime }}{2}\left\{\frac{1}{{\hslash }_{\nu }}{\displaystyle \sum _{i=1}^n{\left[\frac{{x^{″}}_i-{x^{\prime }}_i}{t^{″}-t^{\prime }}-f_i\left(\frac{x^{″}+x^{\prime }}{2}\right)\right]}^2+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}\left(\frac{x^{″}+x^{\prime }}{2}\right)+\frac{1}{{\hslash }_{\mu }}{\displaystyle \sum _{k=1}^m{h}_k^2}\left(\frac{x^{″}+x^{\prime }}{2}\right)}}\right\}\right).\hfill \end{array}$$

(39)

This is the Dirac-Feynman (DF) approximation. The algorithm that follows from applying the DF approximation to the Yau algorithm is the Dirac-Feynman-Yau (DFY) algorithm.

6.2 Example 1

As an example, consider the following continuous-continuous filtering model that has been studied in [23] (and [24])

$$\begin{array}{l}\dot{x}(t)=a\cos (bx(t))+{\sigma }_x\nu (t),\hfill \\ \dot{y}(t)={\tan }^{-1}x(t)+{\sigma }_y\mu (t).\hfill \end{array}$$

(40)

The Lagrangian for this model is given by

$$\frac{1}{2{\sigma }_x^2}{(\dot{x}(t)-a\cos (bx(t)))}^2-\frac{ab}{2}\sin (bx(t))+\frac{1}{2{\sigma }_y^2}{({\tan }^{-1}x(t))}^2.$$

(41)

The parameters chosen were as in the reference. Specifically, a = 1.2, b = 3, σ _x = 0.3, spatial grid spacing Δx = 0.01 and extent [-1.5, 1.5], temporal grid spacing Δt = 0.01 with 200 time steps.

In the first set, the measurement noise was set as σ _y = 0.05. Figure 1 shows a sample of state and measurement processes. The conditional mean, computed using the DFY algorithm, is plotted in Figure 2. Since there was negligible difference in performance between the pre-measurement and post-measurement forms, only the former was employed. Also plotted are 2σ limits. The conditional mean and variance were computed from the computed conditional probability density. The fact that the target was mostly within the 2σ limits of the conditional mean shows that the tracking performance of this algorithm is good.

Figure 1

A sample of state and measurement processes given by Equation 40.

Figure 2

Conditional mean and 2 σ limits computed using the DFY algorithm.

In the next set, the measurement noise was set as σ _y = 0.0125. For this "small noise" case, most of the algorithms studied in [23] failed. In Figure 3 is plotted a sample of signal and measurement processes. The conditional mean and the 2σ limits computed using the DFY algorithm for this instance is plotted in Figure 4.

Figure 3

A sample of state and measurement processes given by Equation 40 – small measurement noise case.

Figure 4

Conditional mean and 2 σ limits computed using the DFY algorithm with small measurement noise.

It is seen that good tracking performance is maintained for the small noise case even when the crudest path integral approximation is used in the Yau algorithm.

6.3 Example 2: Cubic Sensor Problem

The cubic sensor problem is defined by the following signal and measurement model:

$$\begin{array}{l}\dot{x}(t)={\sigma }_x\nu (t),\hfill \\ \dot{y}(t)=x^3(t)+{\sigma }_y\mu (t).\hfill \end{array}$$

(42)

It is a well-studied nonlinear filtering problem because it is one of the simplest examples of a filtering problem that is not finite dimensional (see, for instance, [25] and references therein).

For the simulation of the cubic sensor problem, the following model parameters were chosen (as in [25])

$$\begin{array}{c}\Delta t=0.001,\\ \Delta x=0.01,\\ x_0=N(0,0.01),\\ {\sigma }_x=1,\\ {\sigma }_y=1,\\ T=[0,20].\end{array}$$

(43)

The EKF is a sub-optimal filter which approximates the conditional probability by a Gaussian. For the cubic sensor problem the EKF is given by

$$\begin{array}{c}d\widehat{x}=3{\widehat{x}}^{2}P(dy-{\widehat{x}}^{3}dt),\\ dP=(1-9{\widehat{x}}^4{P}^2)dt.\end{array}$$

(44)

The Yau Lagrangian for the cubic sensor problem is

$$L=\frac{1}{2{\sigma }_x^2}({\dot{x}}^2)+\frac{1}{2{\sigma }_y^2}{x}^6.$$

(45)

The DF approximation the Yau Lagrangian is

$$\tilde{P}(t^{″},x^{″}|t^{\prime },x^{\prime })=\exp \left(-\frac{(t^{″}-t^{\prime })}{2{\sigma }_x^2}{\left(\frac{x^{″}-x^{\prime }}{t^{″}-t^{\prime }}\right)}^2-\left(\frac{t^{″}-t^{\prime }}{2{\sigma }_y^2}\right){\left(\frac{x^{″}+x^{\prime }}{2}\right)}^6\right).$$

(46)

Figure 5 shows the performance of the DFY algorithm. Specifically, the conditional mean along with the two standard deviation bounds computed using the computed conditional probability density is plotted. Observe that the EKF fails completely in this case. As noted in [25], this is because the EKF considers only the first two moments (which vanish here); it is the fourth central moment that plays a crucial role in this example (for the chosen initial condition). Also, note that the state is within the 2σ region for most of the time. This shows that, unlike the EKF, the path integral filter has a reliable error analysis.

Figure 5

A simulation result for the cubic sensor problem.

After an initial period, the performance of the path integral approximation is seen to be excellent and comparable to that obtained using PDE methods in [25]. However, the crucial point is that the path integral solution is equally simple for the higher-dimensional case with more complicated models (e.g., colored noise), whereas a PDE solution would be significantly harder, if not impossible, to implement in real-time.

6.4 Comments

It is remarkable to note that very good performance is obtained using the crudest approximation. Of course, when the time step is large, it will fail (unlike the path integral formula itself). However, the practical situation is that the time steps are often small. Therefore, the DF approximation may be adequate in most cases.

The implementation of this method is trivial. The contrast with other methods, such as those studied in [23], is striking. For instance, many of those methods require off-line computation of complicated partial differential equations with uncertain numerical properties.

The results obtained in this paper used single time-step. More accuracy can be obtained quite simply using multiple time steps. Also, the computation of the transition probability density can be done off-line, but the on-line computation was not an onerous burden for the examples studied here.

It is important to note also that the transition probability density matrix (or tensor in the general case) is sparse (sparsity determined by ħ _μ , ħ _ν ). This is of great importance in higher dimensional filtering problems because

• Sparse matrix storage requirements are small,

• The relevant transition probability density matrix elements can be computed based on the conditional density in the previous step, and

• Sparse matrix computations are very fast.

Note that unlike some other approximation techniques studied in [23], the conditional probability density is obviously always positive (provided, of course, that the initial distribution is positive).

Finally, a comment on measures of performance. Note that a good tracker is one that furnishes not only a good estimate of a state but also provides a reliable measure of the quality of the estimate. For the linear, Gaussian case, the conditional mean and the variance are Gaussian and the Kalman estimates are optimal and provide a complete description. However, for the general nonlinear case, such a concise description is not possible. For instance, the conditional probability density may be highly skewed. It may be multi-modal, in which case the conditional mean is not a very meaningful quantity. A more general measure is to indicate domains of "significant" probability mass; a good filtering solution is then one that guarantees that the state is in the region of significant probability mass with a very high confidence. For the purposes of the paper, the conditional variance was chosen.

7 Conclusion

In this paper, the formal path integral solution to the continuous-continuous nonlinear filtering problem has been presented. The solution is universal, i.e., the initial distribution may be arbitrary. Since the path integral measure is manifestly positive, positivity is maintained if the initial distribution is positive.

A path integral formulation has several advantages. It is well known that Feynman path integrals have led to theoretical insights in other areas including quantum mechanics, quantum field theory and even pure mathematics. It is demonstrated in this paper that it is possible to express the fundamental solution of the YYe in terms of Feynman path integrals. Finally, Feynman path integrals are very suitable for numerical implementation. Practical path integral filtering techniques, especially for solving large dimensional problems, will be presented in subsequent papers.

Appendix 1 continuous-continuous filtering and the Yau equation

In this section, the main results of (continuous-continuous) nonlinear filtering theory are summarized. For the general case (e.g., not the finite-dimensional filter case) and from a practical point of view the most important results are the YYe and the Yau algorithm.

A.1 The Continuous-Continuous Model

The signal and observation model in continuous-continuous filtering is the following:

$$\{\begin{array}{ll}dx(t)=f(x(t),t)dt+e(x(t),t)dv(t),\hfill & x(0)=x_0,\hfill \\ dy(t)=h(x(t),t)dt+n(t)dw(t),\hfill & y(0)=0.\hfill \end{array}$$

(47)

Here x, v, y, and w are ℝⁿ, ℝ^p, ℝ^mand ℝ^qvalued stochastic processes, respectively, and e(x(t), t) ∈ ℝ^n×pand n(t) ∈ ℝ^m×q. These are defined in the Itô sense. The Brownian processes v and w are assumed to be independent with p×p and q ×q covariance matrices Q(t) and W(t), respectively. We denote n(t)W (t)n^T(t) by R(t), a m × m matrix. Also, f is referred to as the drift, e as the diffusion vielbein, and eQe^Tas the diffusion matrix.

In this section, some of the relevant work on continuous-continuous filtering is summarized. Hence, it is assumed that n = p, m = q, f and h are vector-valued C^∞ smooth functions, e(x, t) is an orthogonal matrix-valued C^∞ smooth function, Q(t) is a n × n identity matrix, and n(t) and R(t) are m × m identity matrices. No explicit time dependence is assumed in the model.

A.2 The DMZ Stochastic Differential Equation

The unnormalized conditional probability density, σ (t, x) of the state given the observations {Y (s): 0 ≤ s ≤ t} satisfies the DMZ equation:

$$\begin{array}{cc}d\sigma (t,x)={\mathcal{L}}_0\sigma (t,x)dt+{\displaystyle \sum _{i=1}^n{\mathcal{L}}_i\sigma (t,x)d{y}_i(t),}& \sigma (0,x)={\sigma }_0.\end{array}$$

(48)

Here

$$\begin{array}{c}{\mathcal{L}}_0=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2}{\partial x_i^2}-{\displaystyle \sum _{i=1}^n{f}_i(x)\frac{\partial }{\partial x_i}}}-{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)-\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}},\\ =-{\displaystyle \sum _{i=1}^n\frac{\partial }{\partial x_i}(f_i(x)\cdot )+\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2}{\partial x_i^2}}-\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}},\\ =\frac{1}{2}\left({\displaystyle \sum _{i=1}^n{D}_i^2-\eta (x)}\right),\end{array}$$

(49)

where ${\mathcal{L}}_i$ is the zero-degree differential operator of multiplication by h _i (x), i = 1, ..., m, σ₀ is the probability density of the initial time t₀, and

$$\begin{array}{c}{D}_i\equiv \frac{\partial }{\partial x_i}-f_i(x),\\ \eta (x)\equiv {\displaystyle \sum _{i=1}^n{f}_i^2(x)}+{\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)}+{\displaystyle \sum _{i=1}^m{h}_i^2(x)}.\end{array}$$

(50)

The DMZ equation is to be interpreted in the Stratanovich sense. Note that

$$\begin{array}{c}{D}_i=\frac{\partial }{\partial x_i}-f_i(x),\end{array}$$

(51)

$$\begin{array}{ccc}=e^{F_i}\frac{\partial }{\partial x_i}{e}^{-F_i},& \text{where}& F_i=F_i(x)={\displaystyle {\int }_0^x{f}_i(t)dt}.\end{array}$$

(52)

Hence,

$$D_i^2=e^{F_i}\frac{{\partial }^2}{\partial x_i^2}{e}^{-F_i},$$

(53)

and

$${\mathcal{L}}_0=\frac{1}{2}\left({\displaystyle \sum _{i=1}^n{e}^{F_i}\frac{{\partial }^2}{\partial x_i^2}{e}^{-F_i}-\eta (x)}\right).$$

(54)

A.3 The Robust DMZ Partial Differential Equation

Following [11] and [12] introduce a new unnormalized density

$$u(t,x)=\exp \left(-{\displaystyle \sum _{i=1}^m{h}_i(x)y_i(t)}\right)\sigma (t,x).$$

(55)

Under this transformation, the DMZ SDE is transformed into the following time-varying PDE

$$\{\begin{array}{l}\frac{\partial u}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^{2}u}{\partial x_i^2}(t,u)}+{\displaystyle \sum _{i=1}^n\left(-f_i(x)+{\displaystyle \sum _{j=1}^m{y}_j(t)\frac{\partial h_j}{\partial x_i}}(x)\right)\frac{\partial u}{\partial x_i}(t,x)}\hfill \\ -({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}-\frac{1}{2}{\displaystyle \sum _{i=1}^m{y}_i(t)\Delta h_i(x)+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{y}_i(t)f_j(x)\frac{\partial h_i}{\partial x_j}(x)}}}}\hfill \\ -\frac{1}{2}{\displaystyle \sum _{i,j=1}^m{\displaystyle \sum _{k=1}^n{y}_i(t)y_j(t)\frac{\partial h_i}{\partial x_k}(x)}})u(t,x),\hfill \\ u(0,x)={\sigma }_0(x).\hfill \end{array}$$

(56)

This is called the robust DMZ equation. Here Δ is the Laplacian. The solution of a PDE when the initial condition is a delta function is called the fundamental solution.

A.4 The Yau Equation

Recently, it was proved that the real-time solution of the general nonlinear filtering problem can be obtained reliably [13, 26]. Let $\mathcal{P}$ = {τ₀ <τ₁ < ⋯ <τ _k = τ} be a partition of the time interval [τ₀, τ], and let the norm of the partition ${\mathcal{P}}_k$ be defined as |${\mathcal{P}}_k$| = sup_1≤i≤k{|τ _i - τ_i-1|}. If u _l (t, x) satisfies the equation

$$\{\begin{array}{l}\frac{\partial u_l}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{u}_l}{\partial x_i^2}(t,x)}+{\displaystyle \sum _{i=1}^n\left(-f_i(x)+{\displaystyle \sum _{j=1}^m{y}_j({\tau }_l)\frac{\partial h_j}{\partial x_i}(x)}\right)\frac{\partial u_l}{\partial x_i}(t,x)}\hfill \\ -({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)-\frac{1}{2}{\displaystyle \sum _{i=1}^m{y}_i({\tau }_l)\Delta h_i(x)}}}+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{y}_i({\tau }_l)f_j(x)\frac{\partial h_i}{\partial x_j}(x)}}\\ -\frac{1}{2}{\displaystyle \sum _{i,j=1}^m{\displaystyle \sum _{k=1}^n{y}_i({\tau }_l)y_j({\tau }_l)}\frac{\partial h_i}{\partial x_k}(x)\frac{\partial h_j}{\partial x_k}(x)})u_l(t,x),\\ u_l({\tau }_{l-1},x)=u_{l-1}({\tau }_{l-1},x),\hfill \end{array}$$

(57)

in the time interval τ_l-1 ≤ t ≤ τ _l , then the function ${\tilde{u}}_l$(t, x) defined as

$${\tilde{u}}_l(t,x)=\exp \left({\displaystyle \sum _{i=1}^m{y}_i({\tau }_l)h_i(x)}\right)u_l(t,x)$$

(58)

satisfies the parabolic partial differential equation

$$\begin{array}{c}\frac{\partial {\tilde{u}}_l}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{\tilde{u}}_l}{\partial x_i^2}(t,x)}-{\displaystyle \sum _{i=1}^n{f}_i(x)\frac{\partial {\tilde{u}}_l}{\partial x_i}}(t,x)-\left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}}\right){\tilde{u}}_l(t,x),\\ =-{\displaystyle \sum _{i=1}^n\frac{\partial }{\partial x_i}[f_i(x){\tilde{u}}_l(t,x)]+\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{\tilde{u}}_l}{\partial x_i^2}(t,x)-\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x){\tilde{u}}_l(t,x)},}}\end{array}$$

(59)

in the same time interval. The converse of the statement is also true. In [27], it was also noted that it is sufficient to use the previous observation, i.e., u _l (t, x) satisfies Equation 57 if and only if ${\tilde{u}}_l$(t, x) defined as

$${\tilde{u}}_l(t,x)=\exp \left({\displaystyle \sum _{i=1}^m{y}_i({\tau }_{l-1}})h_i(x)\right)u_l(t,x)$$

(60)

satisfies Equation 59 in the time interval τ_l-1 ≤ t ≤ τ _l . We refer to Equations 58 (60) and Equation 59 as the post-measurement (pre-measurement) forms of the YYe.

Observe that Equation 57 is obtained by setting y(t) to y(τ _l ) in Equation 56. It was proved that the solution of Equation 57 approximates very well the solution of the robust DMZ equation (Equation 56), i.e., it converges to u(t, x) in both pointwise sense and L² sense. Thus, solving Equation 56 is equivalent to solving Equation 59. Finally,

$$\sigma (\tau ,x)=\underset{\left|{\mathcal{P}}_k\right|\to 0}{\lim }{\tilde{u}}_k({\tau }_k,x).$$

(61)

Thus, the solution of the YYe (as |${\mathcal{P}}_k$| → 0) is the desired unnormalized conditional probability density.

Observe that when h(x) = 0, it is simply the FPKfe. However, unlike the FPKfe, the YYe does not satisfy the current conservation condition, i.e., the right-hand term is not a total divergence. This means that it does not conserve probability. This fundamental difference is traced to the fact that the FPKfe evolves the normalized probability density (and preserves the normalization), while the YYe evolves the unnormalized conditional probability density. Therefore, this distinction is made between the two equations in this paper.

A.5 The Yau Algorithm

We may summarize the real-time algorithm, based on both the pre- and post-measurement forms of the YYe, of Yau as follows. Suppose measurements are available at times

⋯ <τ₀ <τ₁ <τ₂ < ⋯ <τ _k = τ.

We seek the solution u _i (t, x), which is the solution of the robust DMZ equation. Let the initial distribution be u(τ₀, x) = σ₀ (x). Then, evolve the initial distribution to the first measurement instant, τ₁, using the YYe:

$$\{\begin{array}{c}\frac{\partial {\tilde{u}}_1}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{\tilde{u}}_1}{\partial x_i^2}(t,x)}-{\displaystyle \sum _{i=1}^n{f}_i(x)\frac{\partial {\tilde{u}}_1}{\partial x_i}}(t,x)-\left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}}\right){\tilde{u}}_1(t,x),\\ {\tilde{u}}_1({\tau }_0,x)={\sigma }_0(x)\{\begin{array}{ll}\exp \left({\displaystyle {\sum }_{j=1}^m[y_j({\tau }_0)-y_j({\tau }_{-1})]h_j(x)}\right)\hfill & (\text{Pre-Measurement})\hfill \\ \exp \left({\displaystyle {\sum }_{j=1}^m[y_j({\tau }_1)-y_j({\tau }_0)]h_j(x)}\right)\hfill & (\text{Post-Measurement}).\hfill \end{array}\end{array}$$

(63)

The solution of equation 63 at time τ₁ is ${\tilde{u}}_1$(τ₁, x). Note that u₁ (τ₁, x) is given by

$$u_1({\tau }_1,x)={\tilde{u}}_1({\tau }_1,x)\{\begin{array}{ll}\exp \left(-{\displaystyle {\sum }_{i=1}^m{y}_i({\tau }_0)h_i(x)}\right)\hfill & (\text{Pre-Measurement})\hfill \\ \exp \left(-{\displaystyle {\sum }_{i=1}^m{y}_i({\tau }_1)h_i(x)}\right)\hfill & (\text{Post-Measurement}).\hfill \end{array}$$

(64)

Next, solve the YYe to the next measurement instant τ₂ with initial condition ${\tilde{u}}_2$(τ₁, x), i.e.,

$$\{\begin{array}{c}\frac{\partial {\tilde{u}}_2}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{\tilde{u}}_2}{\partial x_i^2}(t,x)}-{\displaystyle \sum _{i=1}^n{f}_i(x)\frac{\partial {\tilde{u}}_2}{\partial x_i}}(t,x)-\left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i=1}^m{h}_i^2(x)}}\right){\tilde{u}}_2(t,x),\\ {\tilde{u}}_2({\tau }_1,x)={\tilde{u}}_1({\tau }_1,x)\{\begin{array}{ll}\exp \left({\displaystyle {\sum }_{j=1}^m(y_j({\tau }_1)-y_j({\tau }_0))h_j(x)}\right)\hfill & (\text{Pre-Measurement})\hfill \\ \exp \left({\displaystyle {\sum }_{j=1}^m(y_j({\tau }_2)-y_j({\tau }_1))h_j(x)}\right)\hfill & (\text{Post-Measurement}),\hfill \end{array}\end{array}$$

(65)

to obtain ${\tilde{u}}_2$(τ₂, x). In fact, for i ≥ 2, u _i (τ _i , x) can be computed from ${\tilde{u}}_i$(τ _i , x), where ${\tilde{u}}_i$(t, x) satisfies the equation

$$\{\begin{array}{c}\frac{\partial {\tilde{u}}_i}{\partial t}(t,x)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2\tilde{u}}{\partial x_i^2}(t,x)}-{\displaystyle \sum _{j=1}^n{f}_j(x)\frac{\partial \tilde{u}}{\partial x_j}}(t,x)-\left({\displaystyle \sum _{j=1}^n\frac{\partial f_j}{\partial x_j}(x)+\frac{1}{2}{\displaystyle \sum _{j=1}^m{h}_j^2(x)}}\right){\tilde{u}}_i(t,x),\\ {\tilde{u}}_i({\tau }_{i-1},x)={\tilde{u}}_{i-1}({\tau }_{i-1},x)\{\begin{array}{ll}\exp \left({\displaystyle {\sum }_{j=1}^m(y_j({\tau }_{i-1})-y_j({\tau }_{i-2}))h_j(x)}\right)\hfill & (\text{Pre-Measurement})\hfill \\ \exp \left({\displaystyle {\sum }_{j=1}^m(y_j({\tau }_i)-y_j({\tau }_{i-1}))h_j(x)}\right)\hfill & (\text{Post-Measurement}).\hfill \end{array}\end{array}$$

(66)

The initial condition in Equation 66 follows from noting that (since u _i (τ_i-1, x) = u_i-1 (τ_i-1, x))

$$\begin{array}{c}{\stackrel{\sim }{u}}_i({\tau }_{i-1},x)=u_i({\tau }_{i-1},x)\exp \left({\displaystyle \sum _{j=1}^m{y}_j({\tau }_i)h_j(x)}\right),\\ \begin{array}{cc}=\exp \left(-{\displaystyle \sum _{j=1}^m{y}_j({\tau }_{i-1})h_j(x)}\right){\stackrel{\sim }{u}}_{i-1}({\tau }_{i-1},x)\exp \left({\displaystyle \sum _{j=1}^m{y}_j({\tau }_i)h_j(x)}\right),& (\text{Pre-Measurement})\end{array}\end{array}$$

(67)

$$\begin{array}{c}{\stackrel{\sim }{u}}_i({\tau }_{i-1},x)=u_i({\tau }_{i-1},x)\exp \left({\displaystyle \sum _{j=1}^m{y}_j({\tau }_i)h_j(x)}\right),\\ \begin{array}{cc}=\exp \left(-{\displaystyle \sum _{j=1}^m{y}_j({\tau }_{i-1})h_j(x)}\right){\stackrel{\sim }{u}}_{i-1}({\tau }_{i-1},x)\exp \left({\displaystyle \sum _{j=1}^m{y}_j({\tau }_i)h_j(x)}\right),& (\text{Post-Measurement}).\end{array}\end{array}$$

(68)

Thus the solution of the robust DMZ equation at the ith step is given by

$$u_i({\tau }_i,x)={\tilde{u}}_i({\tau }_i,x)\{\begin{array}{ll}\exp \left(-{\displaystyle {\sum }_{j=1}^m{y}_j({\tau }_{i-1})h_j(x)}\right)\hfill & (\text{Pre-Measurement})\hfill \\ \exp \left(-{\displaystyle {\sum }_{j=1}^m{y}_j({\tau }_i)h_j(x)}\right)\hfill & (\text{Post-Measurement}).\hfill \end{array}$$

(69)

Note that the Yau algorithm is a recursive algorithm as it does not need any previous observation data. Furthermore, the YYe is independent of data and so can be computed off-line, and that the YYe is much simpler than the robust DMZ equation. Finally, note that the output of the Yau algorithm is the desired unnormalized conditional probability density.

Note that while there are several other possible real-time solutions to the continuous-continuous nonlinear filtering problem, all of them assume that the signal and measurement model drifts are bounded (see remarks in [26] and [28]). As a result, those algorithms cannot provide a reliable solution to many real-world problems.