Bayesian Filtering for Nonlinear Stochastic Systems Using
Holonomic Gradient Method with Integral Transform

In most engineering applications, knowledge of the state of a dynamical system is required for monitoring or control via feedback. The optimal filtering theory was developed to estimate the current state of a system using the history of its outputs that are observable but contaminated by random noise. Following the success of the Kalman filter (KF) [1, 2] for linear systems with Gaussian noise, optimal filtering theory has been extended to applications involving nonlinear systems as well as non-Gaussian noise [3, 4, 5, 6, 7].

Bayesian filtering is a general framework of optimal filtering. Although numerous problems can be formulated as the Bayesian filtering problems, they are difficult to implement as algorithms on computers, particularly for the real-time applications in systems and control. The Monte-Carlo scheme, which was first introduced for the particle filter (PF) [8, 9] and subsequently used in other filtering algorithms, such as [7], is a powerful tool to compute integrations accompanied by statistical operations, such as marginalizations and expectations in nonlinear Bayesian filtering problems. Moreover, the nonlinearity of the system dynamics can be exactly accounted for by updating the particles according to the dynamics. To reap these benefits of the Monte-Carlo scheme, however, the number of particles must be sufficiently large, which requires heavy computational resources and may be unacceptable for applications with short sampling intervals. Another approach involves extending the KF to nonlinear cases such as the extended KF (EKF) [10], unscented KF (UKF) [11], cubature KF (CKF) [12], and Gauss-Hermite KF (GHKF) [13]. In all these methods, the posterior probability density function (PDF) is assumed to be or approximated as a Gaussian PDF. Under this Gaussian approximation, the marginalizations and expectations of the nonlinear functions, which are required to update the posterior mean and variance, can be efficiently computed using the Gauss quadrature rules. However, in exchange for efficiency, the Gauss quadrature rules cannot provide exact values of the integrals, which implies that the nonlinearity of the system dynamics can only be considered in the approximate sense. Thus, there is still room to extend the boundary of the trade-off between the exact consideration of nonlinearity and reduction of computational cost.

In nonlinear cases, the integrals in the posterior PDF are ones of general nonlinear functions that depend on parameters such as the observed output, which are both analytically and numerically intractable to evaluate. In our previous work [14], a symbolic-numeric method called the holonomic gradient method (HGM) [15] was used to evaluate the integrals efficiently. By utilizing the symbolic computation in terms of differential operators, we can perform integrations offline while considering the nonlinearity of the system exactly. The state estimation is then reduced to solving a set of initial-value problems (IVPs) of nonlinear ordinary differential equations (ODEs), which can be efficiently accomplished online with numerical integration methods such as the fourth order Runge-Kutta (RK4) or Adams-Bashforth-Moulton predictor-corrector (ABM4) methods. However, our previous method requires solving an IVP for each component of the posterior mean and variance in the one-step estimation. This incurs a high computational cost for the online computation and is unacceptable for short sampling intervals.

To overcome the computational limitations of our previous method, we use an integral transform. Integral transforms are useful tools to compute the moments of a probability distribution efficiently. For example, in [16], Idan and Speyer used an integral transform called the characteristic function of the multivariate Cauchy distribution. For linear systems under additive Cauchy noise, the characteristic function can be analytically propagated in time and can be used to derive the explicit forms of the posterior mean and variance. In this paper, we introduce an integral transform, which is different from, but similar to, the moment generating function. This similarity enables us to compute all the components of the posterior mean and variance efficiently from the integral transform and its derivatives. Moreover, the integral transform and its derivatives can be simultaneously evaluated by using the HGM only once. Hence, in contrast to our previous method, only one IVP must be solved, which leads to the reduction of the online computational cost.

In the Bayesian filtering approach for discrete-time nonlinear systems, the posterior PDF of the state $x_{k}\in{\boldsymbol{\mathrm{R}}}^{n}$ at time step $k$ conditional on the output history $y_{[0:k]}\subset{\boldsymbol{\mathrm{R}}}^{r}$ from time steps $0$ to $k$ is recursively updated via the following equations [17].

We assume that the system dynamics and observation process are defined by the following state and observation equations.

where $w_{k}\in{\boldsymbol{\mathrm{R}}}^{n}$ and $v_{k}\in{\boldsymbol{\mathrm{R}}}^{r}$ denote the system and observation noises that are assumed to be independent and identically distributed with the PDFs $p_{w}(w_{k})$ and $p_{v}(v_{k})$, respectively. In addition, $u_{k}\in{\boldsymbol{\mathrm{R}}}^{m}$ denotes a given input of the system. Using the rule of transformation of PDFs [10], the conditional PDFs $p(x_{k}\mid x_{k-1})$ and $p(y_{k}\mid x_{k})$ in (1) can be rewritten using functions $f$, $h$, $p_{w}$, and $p_{v}$ as follows:

where the former PDF is conditional on $x_{k-1}$ and parametrized by $u_{k}$ owing to (2). By substituting (4) and (5) into (1), we can describe the update law of Bayesian filtering using $f$, $h$, $p_{w}$, and $p_{v}$.

Since Bayesian filtering is a recursive algorithm, we can focus on the one-step update of the posterior PDF and omit the subscript $k$. Moreover, the past history of outputs $y_{[0:k-1]}$ has nothing to do with the update at time step $k$; it only appears in the previous posterior PDF $p(x_{k-1}\mid y_{[0:k-1]})$, so we also omit the past history $y_{[0:k-1]}$. Thus, (1) can be summarized as

where $x$, $y$, and $u$ denote the current state, output, and input, respectively; $x^{-}$ denotes the previous state, and $p_{xy}(x,y;u)$ denotes the joint PDF of $x$ and $y$ defined as the second line in (1).

The update law (6) can be viewed as a functional that maps the previous posterior PDF $p(x^{-})$ to the current one ${p(x\mid y;u)}$. In the linear case with Gaussian noise, this functional can be reduced to simple arithmetic operations of the mean and variance of the previous posterior PDF to yield the prediction and update steps of the Kalman filter. However, in the nonlinear cases, the Gaussian property of the posterior PDF can no longer be guaranteed, which makes it extremely difficult to compute the functional. Therefore, as the first step to overcome this difficulty, we approximate the posterior PDF using a Gaussian.

For the following discussion, $\tilde{p}$ denotes an approximation of the PDF $p$ with the same stochastic variables. Suppose we have a Gaussian $\tilde{p}(x^{-};\mu^{-},\Sigma^{-})=\mathcal{N}(x^{-};\mu^{-},\Sigma^{-})$ approximating $p(x^{-})$, where $\mathcal{N}$ is defined as

By replacing $p(x^{-})$ with $\tilde{p}(x^{-};\mu^{-},\Sigma^{-})$ in the definition of $p_{xy}(x,y;u)$, it can be approximated as

and $p(x\mid y;u)$ can be approximated as

Note that the approximating PDF on the left-hand side is not Gaussian since we consider the nonlinearities of functions $f$ and $h$ without any approximations. Hence, we further approximate $\tilde{p}(x\mid y;u,\mu^{-},\Sigma^{-})$ with a Gaussian $\mathcal{N}(x;\mu,\Sigma)$ of the same mean $\mu$ and variance $\Sigma$ as those of $\tilde{p}(x\mid y;u,\mu^{-},\Sigma^{-})$, that is, $\tilde{p}(x\mid y;u,\mu^{-},\Sigma^{-})\approx\mathcal{N}(x;\mu,\Sigma)$ where

and

Definitions (9) and (10) provide a finite number of functions that map the parameters $y$, $u$, $\mu^{-}$, and $\Sigma^{-}$ to the current estimates $\mu$ and $\Sigma$. Hence, by evaluating the right-hand sides of (9) and (10) recursively, we can perform Bayesian filtering while exactly considering the nonlinearities of $f$ and $h$. For simplicity, we define $\Phi\left[h(x)\right](y,u,\mu^{-},\Sigma^{-})$ for a scalar-, vector-, or matrix-valued function $h(x)$ as

which is an integral over $x$ depending on the parameters $y$, $u$, $\mu^{-}$, and $\Sigma^{-}$. The computations of (9) and (10) are then reduced to evaluations of three integrals, namely $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$, because $\mu$ and $\Sigma$ can be written as $\mu=\Phi[x]/\Phi[1]$ and $\Sigma=\Phi[xx^{\top}]/\Phi[1]-\mu\mu^{\top}$, respectively.

For nonlinear KFs, the posterior PDF is usually approximated as Gaussian [12, 18]. The integrals $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$ are also approximated without considering the effects of higher order moments induced by the nonlinearity of $f$ and $h$. More specifically, the nonlinear KFs first approximate the integral in (7) by computing the mean and variance of the nonlinear function $f(x^{-},u)$ and further approximate the integral in the denominator of (8) in the same way using $h(x)$. By contrast, the proposed method considers the exact integrals using the HGM.

On the other hand, the PF approximates the integrals appearing in (1) using the Monte-Carlo scheme instead of the Gaussian approximation. Therefore, for a sufficiently large number of particles, the PF must exhibit a better accuracy than the proposed method. However, in the numerical example provided in Section V, we observe that the computational cost for the PF is more than that of the proposed method for the same level of accuracy.

As described in the previous section, the update law of the mean and variance can be formulated as the evaluation of three integrals $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$. For linear systems with Gaussian noise, these integrals can be computed analytically using the Gaussian integral. However, for nonlinear systems, we need a numerical method because it is an extremely complex process to compute the integrals analytically. To achieve the trade-off between the exact consideration of nonlinearity and reduction of computational cost, we use a symbolic-numeric method for evaluating the complex integrals, called the HGM. We only provide the outline of the HGM in a specific situation due to limited space.

Let us consider the evaluation of the following integral.

where $\beta$ and $\gamma$ are smooth scalar-valued functions in $X\in{\boldsymbol{\mathrm{R}}}^{n}$ and $Y\in{\boldsymbol{\mathrm{R}}}^{m}$. For example, $\Phi[h(x)](y,u,\mu^{-},\Sigma^{-})$ can be regarded as integral (11) with $\beta=h$, $\gamma=\tilde{p}_{xy}$, $Y=x$, and $X$ corresponds to $(y,u,\mu^{-},\Sigma^{-})$.

We assume that we know $\mathcal{M}$ and $A_{X_{i}}$ in Definition 1 for the holonomic function $\alpha(X)$. Moreover, we assume that the vector $Q(X_{\mathrm{init}})$ at a certain point $X_{\mathrm{init}}$ is known. For example, the vector $Q(X_{\mathrm{init}})$ can be directly computed from (11) if we use sufficient computational resources. Once $Q(X_{\mathrm{init}})$ is obtained, by virtue of (13), we can efficiently evaluate $\alpha(\hat{X})$ at a point $\hat{X}\neq X_{\mathrm{init}}$ without directly evaluating the multiple integral (11). The target point $\hat{X}$ can be arbitrarily selected considering that the following process is successfully performed.

Let $X(s)$ be a smooth curve $[0,1]\ni s\mapsto X(s)\in{\boldsymbol{\mathrm{R}}}^{n}$ with $X(0)=X_{\mathrm{init}}$ and $X(1)=\hat{X}$. Considering the Pfaffian system (13) on this curve, we obtain an IVP

This IVP can be solved to obtain $Q(\hat{X})$ using numerical solution methods such as the RK4 method if none of the denominators in $A_{X_{i}}$ vanish at any $s\in[0,1]$. Consequently, we can evaluate $\alpha(\hat{X})$ as the first component of $Q(\hat{X})$.

The finite set $\mathcal{M}$ and matrices $A_{X_{i}}$, which are supposed to be known, can be computed from a finite set of differential operators called a basis of a zero-dimensional ideal in $\mathcal{R}_{n}$. The following theorem provides a way to find such a basis for a given function and to determine whether the function is holonomic or not (We have skipped the details for simplicity).

It may be difficult to find a basis of a zero-dimensional ideal for $\alpha$ directly from the expression (11). Using symbolic computation, we can compute a basis of a zero-dimensional ideal for $\alpha$ from those for $\beta$ and $\gamma$ [20], which implies the following corollary [20, 19].

Corollary 1, combined with the following assumption, ensures that $\tilde{p}_{xy}(x,y;u,\mu^{-},\Sigma^{-})$ is holonomic.

Under Assumption 1, all the components of the integrals $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$ are holonomic functions of $u$, $y$, $\mu^{-}$, and $\Sigma^{-}$. Therefore, they can be efficiently evaluated using the HGM if the corresponding Pfaffian systems (13) and initial vectors (12) are computed in advance. The authors propose a method to compute (9) and (10) by directly using the HGM [14]. Although the previous method showed better performance than other existing methods for a numerical example, the IVP (14) has to be solved for each component of the three integrals. If we can decrease the number of IVPs, then it reduces the online computational cost. Hereafter, $z\in{\boldsymbol{\mathrm{R}}}^{N}$ denotes a vector consisting of all the independent components of $y\in{\boldsymbol{\mathrm{R}}}^{r}$, $u\in{\boldsymbol{\mathrm{R}}}^{m}$, $\mu^{-}\in{\boldsymbol{\mathrm{R}}}^{n}$, and $\Sigma^{-}\in\mathrm{PD}(n)$, where $N=r+m+n+n(n+1)/2$. With this notation, for example, $\tilde{p}_{xy}(x,z)$ denotes $\tilde{p}_{xy}(x,y;u,\mu^{-},\Sigma^{-})$ as a function of $x$, $y$, $u$, $\mu^{-}$, and $\Sigma^{-}$.

In this work, we use an integral transform of $\tilde{p}_{xy}$ to reduce the number of IVPs solved in the one-step estimation. Consider an integral transform $\mathcal{T}[\tilde{p}_{xy}]$ defined as follows.

Note that although this definition is similar to that of the moment generating function, it is different because $\tilde{p}_{xy}$ is the PDF of the $x$ and $y$ pair rather than $x$. We assume that there exists a compact subset of ${\boldsymbol{\mathrm{R}}}^{n}$ including the origin such that for all $\xi$ in the subset, (16) can be defined and smooth.

The integrals $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$ are then obtained from the derivatives of $\mathcal{T}[\tilde{p}_{xy}]$ at $\xi=0$ as follows:

Here, the approximating PDF $\tilde{p}_{xy}$ is a holonomic function under Assumption 1, and the kernel $\exp(\xi^{\top}x)$ is also holonomic. This implies that $\mathcal{T}[\tilde{p}_{xy}]$ is also a holonomic function by Corollary 1. Therefore, $\mathcal{M}$ exists such that the vector-valued function $Q(X)$ defined by $\mathcal{M}$ and $\mathcal{T}[\tilde{p}_{xy}]$ satisfies the Pfaffian system (13). Using the Pfaffian system, we can explicitly express $\mathcal{T}[\tilde{p}_{xy}]$ and its derivatives (17) as well as (18) as a linear combination of the components of $Q$ with coefficients in ${\boldsymbol{\mathrm{R}}}(\xi,z)$.

First, $\mathcal{T}[\tilde{p}_{xy}](\xi,z)$ itself is the first component of $Q(\xi,z)$, thus satisfying the following:

where $C^{(0)}$ is a constant coefficient vector $[1\ 0\ \cdots\ 0]\in{\boldsymbol{\mathrm{R}}}^{q}$. Next, the first derivatives $\partial_{\xi_{i}}\mathcal{T}[\tilde{p}_{xy}]\ (i=1,\dots,n)$ are obtained using the first components of the left-hand sides of (13). Hence, the following identities hold.

where the coefficient vector $C^{(1)}_{i}(\xi,z)$ is the first row vector of $A_{\xi_{i}}(\xi,z)\in{\boldsymbol{\mathrm{R}}}(\xi,z)^{q\times q}$ in (13). Finally, the second derivatives $\partial_{\xi_{i}}\partial_{\xi_{j}}\mathcal{T}[\tilde{p}_{xy}]\ (i,j=1,\dots,n)$ are obtained by differentiating both sides of (13); $\partial_{\xi_{i}}\bullet\partial_{\xi_{j}}Q=\partial_{\xi_{i}}\bullet(A_{\xi_{j}}Q)=\partial_{\xi_{i}}A_{\xi_{j}}Q+A_{\xi_{j}}\partial_{\xi_{i}}Q=(\partial_{\xi_{i}}A_{\xi_{j}}+A_{\xi_{j}}A_{\xi_{i}})Q$, and hence

where the coefficient vector $C^{(2)}_{ij}(\xi,z)$ is the first row vector of $\partial_{\xi_{j}}A_{\xi_{i}}+A_{\xi_{i}}A_{\xi_{j}}$.

We assume the explicit expressions of $C^{(1)}_{i}(\xi,z)$ and $C^{(2)}_{ij}(\xi,z)$ are obtained. Using the expressions, we can compute the values of $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$ at $z=\hat{z}$ from (19)–(21) if we have the vector $Q(0,\hat{z})$ for a given $\hat{z}$. The computation of $Q(0,\hat{z})$ can be accomplished by numerically integrating (14) using numerical integration methods, such as the RK4 method. Therefore, the current mean $\mu$ and variance $\Sigma$ for a given data $\hat{z}$ can be computed from (9) and (10) as

In the previous method [14], IVP (14) has to be solved for each of the three integrals $\Phi[1]$, $\Phi[x]$, and $\Phi[xx^{\top}]$ even for the one-dimensional case, that is, three IVPs have to be solved in total. By contrast, only a single IVP for $\mathcal{T}[\tilde{p}_{xy}]$ needs to be solved in the proposed method. This decrease in the number of IVPs leads to a reduction in the computational cost, which can be seen in Section V.

Now, we summarize the offline and online processes of the proposed method as Algorithms 1 and 2. The details of each step in Algorithm 1 can be found in [20, 14].

This section presents a numerical example to show the efficiency of the proposed method. In the following demonstration, we use Risa/Asir [21] for the symbolic computation, Maple for the offline numerical computation, and Python for the online numerical computation on a PC (Intel(R) Core(TM) i7-1065G7 CPU @ 1.30 GHz; RAM: 16 GB).