Gauge Freedom within the Class of Linear Feedback Particle Filters

TL;DR

This paper explores the gauge freedom in linear feedback particle filters, revealing how different particle transformations can preserve the posterior distribution, thus extending the class of known linear FPFs.

Contribution

It characterizes gauge freedom in linear FPFs for Gaussian filtering, expanding the set of parametrized linear feedback particle filters.

Findings

Identifies gauge transformations that leave the posterior invariant.

Extends known families of linear FPFs through gauge characterization.

Provides a framework for understanding equivalence classes of particle filters.

Abstract

Feedback particle filters (FPFs) are Monte-Carlo approximations of the solution of the filtering problem in continuous time. The samples or particles evolve according to a feedback control law in order to track the posterior distribution. However, it is known that by itself, the requirement to track the posterior does not lead to a unique algorithm. Given a particle filter, another one can be constructed by applying a time-dependent transformation of the particles that keeps the posterior distribution invariant. Here, we characterize this gauge freedom within the class of FPFs for the linear-Gaussian filtering problem, and thereby extend previously known parametrized families of linear FPFs.

Equations

$$d{X}_t$$

$$d{Y}_t$$

$$d{X}_t$$

$$d{Y}_t$$

$$d{\mu }_t$$

$$d{P}_t$$

$$d{S}_t=f(S_t)dt+g(S_t)d{\bar{W}}_t+u_t(S_t)dt+K_t(S_t)\circ d{Y}_t,$$

$$\nabla \cdot (p(x)K^{.j}(x))=-\left(h^j(x)-{\hat{h}}^j\right)p(x)$$

$$d{S}_t=A{S}_{t}dt+P_t{C}^{\top }\left(d{Y}_t-\frac{1}{2}C(S_t+{\mu }_t)dt\right)+Bd{\bar{W}}_t.$$

$$\begin{gathered} d{S}_t=A{S}_{t}dt+P_t{C}^{\top }\left(d{Y}_t-\frac{1}{2}C(S_t+{\mu }_t)dt\right) \\ +{\hat{\Omega }}_t{P}_t^{-1}(S_t-{\mu }_t)dt+\frac{1}{2}B{B}^{\top }{P}_t^{-1}(S_t-{\mu }_t)dt, \end{gathered}$$

$$\hat{\Omega}_{t}P_{t}^{-1}+P_{t}^{-1}\hat{\Omega}_{t}=A^{\top}-A\\ +\tfrac{1}{2}\Big{(}P_{t}^{-1}BB^{\top}-BB^{\top}P_{t}^{-1}\\ +P_{t}C^{\top}C-C^{\top}CP_{t}\Big{)}.$$

$$g_t{P}_t{g}_t^{\top }=P_t,t\ge 0,$$

$$d{g}_t=\left(M_{t}dt+\sum _{i=1}^m{{\rm Y}}_t^{(i)}d{Y}_t^{(i)}\right)P_t^{-1}{g}_t,$$

$$M_t=\frac{1}{2}\left({\dot{P}}_t-g_t{\dot{P}}_t{g}_t^{\top }-\sum _{i=1}^m{{\rm Y}}_t^{(i)}{P}_t^{-1}({{\rm Y}}_t^{(i)}{)}^{\top }\right)+{{\rm Y}}_t^{(0)},$$

$$d{E}_t=\left(G_{t}dt+\sum _{i=1}^m{\Omega }_t^{(i)}d{Y}_t^{(i)}\right)P_t^{-1}{E}_t+\sum _{j=1}^r{H}_t^{(j)}d{\bar{W}}_t^{(j)},$$

$$G_{t}=\frac{1}{2}\Bigg{(}\dot{P}_{t}+\sum_{i=1}^{m}\Omega^{(i)}_{t}P_{t}^{-1}\Omega^{(i)}_{t}\\ -\sum_{j=1}^{r}H^{(j)}_{t}(H^{(j)}_{t})^{\top}\Bigg{)}+\Omega^{(0)}_{t},$$

$$\begin{gathered} d{S}_t=A{S}_{t}dt+P_t{C}^{\top }\left(d{Y}_t-\frac{1}{2}C(S_t+{\mu }_t)dt\right) \\ +{\tilde{\Omega }}_t{P}_t^{-1}(S_t-{\mu }_t)dt+\frac{1}{2}B{B}^{\top }{P}_t^{-1}(S_t-{\mu }_t)dt \\ +\sum _{j=1}^r{H}_t^{(j)}d{\bar{W}}_t^{(j)}-\frac{1}{2}\sum _{j=1}^r{H}_t^{(j)}(H_t^{(j)}{)}^{\top }{P}_t^{-1}(S_t-{\mu }_t)dt, \end{gathered}$$

$${\tilde{\Omega }}_t^{(0)}$$

$$+\frac{1}{2}\sum _{i=1}^m\left({{\rm Y}}_t^{(i)}{P}_t^{-1}{\tilde{\Omega }}_t^{(i)}-{\tilde{\Omega }}_t^{(i)}{P}_t^{-1}{{\rm Y}}_t^{(i)}\right),$$

$${\tilde{\Omega }}_t^{(i)}$$

$${\tilde{H}}_t^{(j)}$$

$$d{g}_t{E}_t$$

$$g_{t}d{E}_t$$

$$+\sum _{j=1}^r{g}_t{H}_t^{(j)}d{\bar{W}}_t^{(j)},$$

$$d{g}_{t}d{E}_t$$

$$\begin{gathered} g_t{G}_t{g}_t^{\top }+M_t=\frac{1}{2}{g}_t{\dot{P}}_t{g}_t^{\top }+\frac{1}{2}\sum _{i=1}^m{g}_t{\Omega }_t^{(i)}{P}_t^{-1}{\Omega }_t^{(i)}{g}_t^{\top } \\ -\frac{1}{2}\sum _{j=1}^r{g}_t{H}_t^{(j)}(g_t{H}_t^{(j)}{)}^{\top }+g_t{\Omega }_t^{(0)}{g}_t^{\top } \\ +\frac{1}{2}{\dot{P}}_t-\frac{1}{2}{g}_t{\dot{P}}_t{g}_t^{\top }-\frac{1}{2}\sum _{i=1}^m{{\rm Y}}_t^{(i)}{P}_t^{-1}({{\rm Y}}_t^{(i)}{)}^{\top }+{{\rm Y}}_t^{(0)}. \end{gathered}$$

$$g_{t}G_{t}g_{t}^{\top}+M_{t}=\frac{1}{2}\Big{(}\dot{P}_{t}+\sum_{i=1}^{m}\tilde{\Omega}^{(i)}_{t}P_{t}^{-1}\tilde{\Omega}^{(i)}_{t}\\ -\sum_{j=1}^{r}\tilde{H}^{(j)}_{t}(\tilde{H}^{(j)}_{t})^{\top}\Big{)}+g_{t}\Omega^{(0)}_{t}g_{t}^{\top}+\Upsilon^{(0)}_{t},$$

$${\mathcal{L}}_{P,r}({\Omega }^{(0)},...,{\Omega }^{(m)},H^{(1)},...,H^{(r)})=\text{tr}G{P}^{-1}{G}^{\top },$$

$${\Omega }_t^{(0)}={\Omega }_t^{\ast }=\underset{\Omega \in \text{Skew}(n)}{\mathrm{argmin}}{\mathcal{L}}_{P_t,0}(\Omega ,0,...,0)$$

$$d{E}_t=\left(A+\frac{1}{2}B{B}^T-\frac{1}{2}{P}_t{C}^{\top }C{P}_t+{\hat{\Omega }}_t\right)P_t^{-1}{E}_{t}dt,$$

$${\Omega }_t^{(0)}=\frac{1}{2}\left(A{P}_t-P_t{A}^{\top }+{\hat{\Omega }}_t\right).$$

$$\text{tr}\left[X\left(P^{-1}{G}^{\top }-G{P}^{-1}\right)\right]=0,\forall X\in \text{Skew}(n).$$

$$\Omega P^{-1}+P^{-1}\Omega=\tfrac{1}{2}\left(P^{-1}\dot{P}-\dot{P}P^{-1}\right)\\ =\tfrac{1}{2}\Big{(}P^{-1}BB^{\top}-BB^{\top}P^{-1}+P^{-1}AP+A^{\top}\\ -A-PA^{\top}P^{-1}+PC^{\top}C-C^{\top}CP\Big{)}.$$

$$\sum _{j=1}^r{H}_t^{(j)}(H_t^{(j)}{)}^{\top }-\sum _{i=1}^m{\Omega }_t^{(i)}{P}_t^{-1}{\Omega }_t^{(i)}={\dot{P}}_t.$$

$$\begin{gathered} 0=\text{tr}\left[K^{(0)}\left(P^{-1}{G}^{\top }-G{P}^{-1}\right)\right] \\ +\sum _{i=1}^m\text{tr}\left[K^{(i)}{P}^{-1}{\Omega }^i\left(P^{-1}(G{)}^{\top }+G{P}^{-1}\right)\right] \\ -\sum _{j=1}^m\text{tr}\left[l^{(j)}(H^{(j)}{)}^{\top }\left(P^{-1}(G{)}^{\top }+G{P}^{-1}\right)\right] \end{gathered}$$

Full text

I Introduction

The filtering problem is the problem of estimating a quantity evolving in time that is accessible only through partial and noisy observations. This is commonly formalized as the problem of finding the conditional distribution of the hidden state at time $t$ given all observations up to that time. Despite its long history and rich developments around its theory (see [1], Section 1.3), the principal challenge of implementing filters in practical applications centers around the lack of closed-form solutions and the resulting necessity to find efficient and robust numerical approximations.

The issue has become even more severe in the era of big data where the dimensionality of the processes is very large. The main approach to approximate the conditional distribution, sequential Monte-Carlo or particle filtering (see [2] for a survey and pointers to the literature), is known to exhibit a curse of dimensionality as the number of dimensions of the observations grows [3]-[6]. The problem can be traced down to the use of importance weights and their increasing degeneracy as time progresses (see [6] and the references therein).

More recently, particle filters without importance weights [7]-[11] have been gaining attention. While lacking the principal vulnerability of weighted particle filters, many theoretical questions remain open [12]-[14]. One question that has recently received some interest is the non-uniqueness of the law of the process approximating the filtering distribution. It has been shown [15]-[16] that at least in the linear-Gaussian case there are many ways to construct particle dynamics – both deterministic and stochastic – that track the exact posterior distribution (given by the Kalman-Bucy filter).

Here, we systematically investigate the degrees of freedom in choosing dynamics of the particles (within some constraints) for the linear-Gaussian filtering problem while keeping their distribution aligned with the exact conditional distribution (assuming that the initial distribution of the filter is Gaussian). We characterize a group of transformations (which we call gauge transformations, taking inspiration by similar transformations appearing in theoretical physics) that preserve the conditional distribution and describe how it acts on a class of linear feedback particle filters that extends the ones in the literature. Moreover, we propose a cost function on the family of all such filters. We identify a known and a new feedback particle filter as a (constrained) optimum of this cost function.

The remainder of the paper is structured as follows: in Section II, we discuss the filtering problem in general and the linear-Gaussian case that is the focus of this paper in particular. We define the notion of particle filter that is the focus of this work and provide an overview of previous work regarding such filters and the history of the non-uniqueness problem. In Section III, we discuss the symmetries (gauge transformations) of the problem, define certain classes of filters on which the symmetries act, and describe this action in detail. In Section III-C we show how existing linear feedback particle filters arise as special cases of the class of linear feedback particle filters defined earlier. In Section IV, we introduce an optimality criterion and optimize it with and without constraints, obtaining two specific types of feedback particle filters. Lastly, in Section V, we discuss the implications of our results and comment on future directions.

II Preliminaries and background

The classical filtering problem is to find the conditional distribution of the hidden state $X_{t}$ given observations $\mathscr{F}^{Y}_{t}$ for the stochastic system given by

[TABLE]

where $X_{t}$ and $Y_{t}$ are valued in $\mathds{R}^{n}$ and $\mathds{R}^{m}$ respectively, $f,g,h$ are (known) vector- or matrix-valued functions satisfying suitable conditions for the well-posedness of the stochastic differential equations (SDEs), and $W_{t}$ and $V_{t}$ are independent $n^{\prime}$ and $m$-dimensional Brownian motions respectively. The distribution of $X_{0}$ is assumed to be known and independent of the Brownian motions.

II-A Linear-Gaussian filtering problem

In this paper, we will focus on the special linear-Gaussian case of the above problem, where $f,g$, and $h$ are chosen such that

[TABLE]

for some $A\in\mathds{R}^{n\times n}$, $B\in\mathds{R}^{n\times n^{\prime}}$, and $C\in\mathds{R}^{m\times n}$, and $X_{0}$ has Gaussian distribution. This filtering problem has an exact solution, called the Kalman-Bucy filter [17]. The conditional distribution of $X_{t}$ given $\mathscr{F}^{Y}_{t}$ is multivariate Gaussian with mean $\mu_{t}$ and covariance matrix $P_{t}$ which jointly evolve as

[TABLE]

where $\mu_{0}$ and $P_{0}$ are chosen according to the distribution of $X_{0}$. For convenience, throughout this article it will be assumed that $P_{0}$ is strictly positive definite. Under this assumption, $P_{t}$ remains strictly positive definite for all $t\geq 0$ (see Proposition 1.1 in [18]).

II-B Particle filters

In the context of this paper, a particle filter is any approximation of the conditional distribution of $X_{t}$ given $\mathscr{F}^{Y}_{t}$ by a set of (unweighted) samples or particles $S_{t}^{(i)}$ for $i=1,...,N$, such that $S_{t}^{(i)}$ are $\mathscr{F}^{Y,Z}_{t}$-adapted processes.111This use of the term ‘particle filter’ differs of the usual one, where the samples are weighted by importance weights. In this paper, we shall only be concerned with unweighted particle filters, which are also known as interacting particle systems or ensemble Kalman filters [19, 7]. Here, $\mathscr{F}^{Y,Z}_{t}$ is the filtration generated by the process $(Y_{t},Z_{t})$, where $Z_{t}$ is some process independent from $\mathscr{F}^{X,Y}_{t}$ (for example, some additional noise in the particle dynamics). In the following, we will only consider ‘symmetric’ particle filters for which all particles have the same conditional distribution given $\mathscr{F}^{Y}_{t}$. We will therefore talk about the dynamics of a single representative particle $S_{t}$, omitting the particle index $i$. The distribution of $S_{t}$ will still depend on the number of particles $N$. Such a particle filter is called asymptotically exact if the distribution of $S_{t}$ given $\mathscr{F}^{Y}_{t}$ converges to the conditional distribution of $X_{t}$ as $N\to\infty$.

II-C Feedback particle filter (FPF)

An asymptotically exact filter has been found in [9] and [10] based on mean-field optimal control. It is usually referred to simply as Feedback Particle Filter, but in order to distinguish it from similar algorithms, we will refer to this specific algorithm as stochastic feedback particle filter (sFPF). Specifically, the sFPF is derived by finding control terms $u$ and $K$ for the particle dynamics222The FPF is more naturally expressed in Stratonovich form (notation: $\circ$).

[TABLE]

such that the particle filter whose particles $S_{t}^{(i)}$ evolve according to the SDE (7) is asymptotically exact. Here, $\bar{W}_{t}$ is an $n^{\prime}$-dimensional Brownian motion that is independent of $\mathscr{F}^{X,Y}_{t}$ and plays the role of $Z_{t}$ in the previous paragraph. The sFPF can be derived by ‘aligning’ the Fokker-Planck equation of the particle filter with the Kushner-Stratonovich filtering equation in an appropriate sense (see [10], Section 2), which leads to a McKean-Vlasov or mean-field equation where $u_{t}$ and $K_{t}$ depend on the distribution of $S_{t}$. As a result, the FPF is exact in the mean-field sense if initialized with the correct initial distribution. However, the computation of the gain term $K$ requires the solution of a linear boundary value problem (BVP) at each instant in time, which accounts for the bulk of the computational cost this algorithm.

Specifically, the $j$’th column of the gain matrix $K$ is any vector-valued function that satisfies a weighted Poisson equation

[TABLE]

with suitable boundary conditions, where $p$ is the density of the current particle distribution and $\hat{h}^{j}=\int h^{j}(y)p(y)dy$. This equation does not have uniqueness of solutions, as any solution $K$ may be perturbed by a divergence-free vector field $K^{\prime}$, i.e. $\nabla\cdot(pK^{\prime})=0$ to produce a new solution $K+K^{\prime}$. Uniqueness can be obtained by requiring $K$ to be a gradient333Under suitable conditions on $p$, this gradient solution can be interpreted as the minimum-energy gain, which is related to the dynamic (Benamou-Brenier) formulation of optimal transport (see Lemma 8.4.2 in [20])., i.e. $K=\nabla\phi$. Whatever $K$ is chosen, the general solution for $u$ can be written as $u=-\tfrac{1}{2}(h+\hat{h})K+u^{\prime}$, where $u^{\prime}$ is divergence-free in the sense above.

In the linear-Gaussian case, the gradient-form solution of the BVP (8) can be obtained in closed form. It is given by the Kalman gain $K_{t}\equiv P_{t}C^{\top}$ (see [10], Section 4.2). By choosing $u_{t}=-\tfrac{1}{2}(h+\hat{h})K_{t}$, the resulting FPF – which will henceforth be referred to as the stochastic linear FPF (slFPF) – is equivalent to the square-root form of the Ensemble Kalman-Bucy filter. It takes the form

[TABLE]

If $S_{0}$ is normally distributed with mean $\mu_{0}$ and covariance matrix $P_{0}$, the conditional distribution of $S_{t}$ given $\mathscr{F}^{Y}_{t}$ has mean and covariance $\mu_{t}$ and $P_{t}$ respectively for all $t\geq 0$.

II-D Non-uniqueness of the control law

The linear case allows us to explore the non-uniqueness inherent in the choice of $K$ and $u$. In the linear case, a larger class of solutions of Eq. (8) can be obtained by adding a linear divergence-free field to each column of $K$, i.e. setting $K^{.j}_{t}(x)=P_{t}(C^{j.}+\Pi^{(j)}_{t}(x-\mu_{t}))$, where $\Pi^{(j)}_{t}$ are skew-symmetric matrices. In addition, one may choose $u_{t}(x)=-\tfrac{1}{2}K_{t}(x)C(x+\mu_{t})+P_{t}\tilde{\Pi}_{t}(x-\mu_{t})$, where $\tilde{\Pi}_{t}$ is yet another skew-symmetric matrix.

Besides the modification of $K$ and $u$ by divergence-free terms, there are other modifications of the slFPF dynamics (9) for which the distribution of $S_{t}$ agrees with the Kalman-Bucy filter. In [15] a deterministic linear FPF (deterministic refers to the fact that there is no independent noise term) was derived from an optimal transport perspective. We will refer to it as optimal transport deterministic linear feedback particle filter (OTdetFPF). Its dynamics are given by

[TABLE]

where $\hat{\Omega}_{t}$ is the unique skew-symmetric matrix satisfying

[TABLE]

The OTdetFPF replaces the Brownian motion term of the sFPF by a repulsive term that drives the particles away from their mean. In addition, it has a skew-symmetric term that arises from the minimization of the transportation cost. Since the skew-symmetric term does not affect the distribution, a parametrized family of deterministic FPFs (detFPF) may be constructed by choosing it arbitrarily, in particular by setting it to zero.

In another paper [16], an approach based on optimal control theory and duality formalisms is used to extend the linear FPF to a two-parameter family of FPFs. In the next section, all of the previously mentioned transformations will be captured by a common parametrization, see in particular Section III-C.

II-E Notations and definitions

We denote by $S$ the process $(S_{t})_{t\geq 0}$ as a random variable with values in a suitable subspace of functions $\mathds{R}\to\mathds{R}^{n}$ (for most cases, the space of continuous functions will suffice), with its law given by a probability measure on the corresponding Borel $\sigma$-algebra. If the conditional distribution of $S_{t}$ given $\mathscr{F}^{Y}_{t}$ agrees with the conditional distribution of $X_{t}$ given $\mathscr{F}^{Y}_{t}$ for all $t\geq 0$, such a process is called particle filter444We only consider exact particle filters here. The idea is that fore finite $N$ the mean-field quantities appearing in the control law are replaced by empirical estimates, yielding an asymptotically exact filter (although the details are nontrivial, see e.g. [21]). If $S$ is adapted to $\mathscr{F}^{Y,Z}_{t}$, where $Z$ consists of $r\geq 0$ Brownian motions independent of $Y$ and if $S$ is the solution to an SDE, it will be called feedback particle filter. If $r=0$ it is called deterministic feedback particle filter. We denote by $\text{Skew}(n)$ the real vector-space of skew-symmetric $n\times n$-matrices, i.e. of matrices $X$ such that $X^{\top}=-X$, where $\cdot^{\top}$ denotes transposition. The notation $\text{GL}(n)$ is used for the multiplicative group of invertible $n\times n$-matrices.

III The class of linear particle filters and their symmetries

In the linear-Gaussian case, a particle filter is exact if the conditional distribution of $S_{t}$ given $\mathscr{F}^{Y}_{t}$ is multivariate Gaussian with mean equal to $\mu_{t}$ and covariance matrix equal to $P_{t}$ (the mean and covariance matrix of the Kalman-Bucy filter). If we split off the mean and define $E_{t}:=S_{t}-\mu_{t}$, the process $E_{t}$ can be any Gaussian process with conditional mean equal to zero and conditional covariance function $k(t,t^{\prime})$ such that $k(t,t)=P_{t}$ for all $t\geq 0$. In this section, we first look at a group of linear time-dependent transformations that preserve this structure, called gauge transformations. Then we consider the subclass $\mathfrak{F}$ of filters for which $E_{t}$ is the solution to an SDE, which we will call the class of linear feedback particle filters, and the subclass $\mathfrak{F}_{\text{det}}\subset\mathfrak{F}$ for which $E_{t}$ does not have an independent source of noise, called deterministic linear feedback particle filters. We then describe the action of the group on $\mathfrak{F}$, proving that the group of gauge transformations acts transitively on $\mathfrak{F}_{\text{det}}$.

III-A Gauge freedom of linear particle filters

From any given particle filter $S_{t}=\mu_{t}+E_{t}$, another particle filter $\tilde{S}_{t}=\mu_{t}+\tilde{E}_{t}$ can be constructed by applying a linear transformation $\tilde{E}_{t}=g_{t}E_{t}$, where $g_{t}\in\text{GL}(n)$ is adapted to $\mathscr{F}^{Y}_{t}$. By construction, the conditional means of $\tilde{S}_{t}$ and $S_{t}$ agree, and if in addition $g_{t}$ satifies

[TABLE]

the conditional variances of $\tilde{S}_{t}$ and $S_{t}$ also agree.

For fixed $t$, the subset consisting of those $g\in\text{GL}(n)$ that satisfy Eq. (12) forms a (random and time-dependent) Lie subgroup of $\text{GL}(n)$, which we denote by $\mathcal{G}_{t}$. Its Lie algebra $\mathfrak{g}_{t}$ consists of all matrices of the form $\Omega P_{t}^{-1}$, where $\Omega\in\text{Skew}(n)$ is arbitrary.

In principle, the choice of the function $g:[0,\infty)\to GL(n)$, $t\mapsto g_{t}$ does not have to be regular. For example, after simulating a sample of particles up to time $t$, a sample for $s\leq t$ can be modified by a transformation $g\in\mathcal{G}_{s}$ without concern for samples at other times as long as only information up to time $s$ is used. However, in the following we will restrict to choices of $g$ with more regularity.

For example, $g$ may be chosen to be a continuously differentiable function governed by an ordinary differential equation (ODE) $\dot{g}_{t}=M_{t}P_{t}^{-1}g_{t}$. This parametrization is convenient because Eq. (12) implies the constraint $M_{t}=\tfrac{1}{2}(\dot{P}_{t}-g_{t}\dot{P}_{t}g_{t}^{\top})+\Upsilon^{(0)}_{t}$, where $\Upsilon^{(0)}$ is a continuous but otherwise arbitrary function with values in $\text{Skew}(n)$. The skew-symmetric component of $M_{t}$ accounts for motion of $g_{t}$ along the group $\mathcal{G}_{t}$, whereas the symmetric component is due to the change in $\mathcal{G}_{t}$ itself.

More generally, $g_{t}$ may be given by the solution of an SDE involving the observations, i.e.

[TABLE]

where $Y^{(i)}_{t}$ denotes the component number $i$ of $Y_{t}$. By differentiating Eq. (12) and matching terms, we find that $\Upsilon^{(i)}_{t}\in\text{Skew}(n)$, $i=1,...,m$ can be chosen arbitrarily, whereas $M_{t}$ has to satisfy the constraint

[TABLE]

where $\Upsilon^{(0)}_{t}\in\text{Skew}(n)$ is again arbitrary. We call any $g$ of the above form (13)-(14) a deterministic gauge transformation.

III-B A general class of linear feedback particle filters

Suppose that $S_{t}=\mu_{t}+E_{t}$, where the dynamics of $E_{t}$ are given by

[TABLE]

where $r\in\{0,1,2,...\}$ (for $r=0$ the sum over $j$ is zero by convention), $\bar{W}^{(j)}_{t}$ are scalar Brownian motions independent of each other and of $\mathscr{F}^{X,Y}_{t}$, and $E_{0}$ has zero mean and covariance matrix equal to $P_{0}$. In addition, all coefficients are assumed to be adapted to $\mathscr{F}^{Y}_{t}$ and bounded. By construction, $S_{t}$ has conditional mean equal to $\mu_{t}$. Moreover, $S_{t}$ has conditional covariance equal to $P_{t}$ for all $t\geq 0$ if and only if $\Omega^{(i)}_{t}\in\text{Skew}(n)$, $i=1,...,m$, and $G_{t}$ satisfies the constraint

[TABLE]

where $\Omega^{(0)}_{t}\in\text{Skew}(n)$ is arbitrary.

Definition III.1

The class of all processes $S_{t}=\mu_{t}+E_{t}$ with $E_{t}$ having dynamics according to Eqs. (15)-(16) where $\Omega^{(i)}_{t}$, $i=0,1,...,m$ are $\mathscr{F}^{Y}_{t}$-adapted $\text{Skew}(n)$-valued processes, is denoted by $\mathfrak{F}$ and called the class of linear feedback particle filters. The subclass $\mathfrak{F}_{\text{det}}\subset\mathfrak{F}$ for which $r=0$ is called deterministic linear feedback particle filters.

III-C Special cases of linear feedback particle filters

First, it is worth pointing out that the filters obtained by adding linear divergence-free vector fields to $K$ and $u$ (as described in the beginning of Section II-D) belong to $\mathscr{F}$, as they can be related to the above parametrization by setting $\Omega^{(i)}_{t}=P_{t}\Pi^{(i)}_{t}P_{t}$ and $\Omega^{(0)}_{t}=\tfrac{1}{2}AP_{t}-\tfrac{1}{2}P_{t}A-P_{t}\tilde{\Pi}_{t}P_{t}-\tfrac{1}{2}\sum_{i=1}^{m}C^{i.}\mu_{t}P_{t}\Pi^{(i)}_{t}P_{t}$.

By adding the mean dynamics ($d\mu_{t}$) to $dS_{t}$ from Eqs. (15)-(16), we obtain the dynamics of particle filters in $\mathfrak{F}$:

[TABLE]

where $\tilde{\Omega}_{t}=\tfrac{1}{2}\Omega^{(0)}_{t}-\tfrac{1}{2}(AP_{t}-P_{t}A^{\top})+\tfrac{1}{2}\sum_{i=1}^{m}\Omega^{(i)}_{t}P_{t}^{-1}\Omega^{(i)}_{t}$.

Note that although $\tilde{\Omega}_{t}$ is comprised of several terms, the presence of the arbitrary skew-symmetric matrix $\Omega^{(0)}_{t}$ means that there are no constraints on its time evolution. We may thus split it as $\tilde{\Omega}_{t}=\hat{\Omega}_{t}+(\tilde{\Omega}_{t}-\hat{\Omega}_{t})$, where $\hat{\Omega}_{t}$ is given by the solution to Eq. (11). Thus, for $r=0$ (or $r>0$ and $H^{(j)}_{t}=0$, $j=1,...,r$) and $\tilde{\Omega}_{t}-\hat{\Omega}_{t}=0$, Eq. (17) takes the form of the OTdetFPF from Eq. (10). If we set $r=n^{\prime}$, $H^{(j)}_{t}=j\text{'th column of }B$, and $\tilde{\Omega}_{t}=0$, the last term in the second line of (17) cancels with the last term in the third line, and we obtain the slFPF from Eq. (9). As a last example, if $r=n^{\prime}+m$, $H^{(j)}_{t}=\gamma_{1}\cdot j\text{'th column of }B$ for $j=1,...,n^{\prime}$, and $H^{(n^{\prime}+j)}_{t}=\gamma_{2}\cdot j\text{'th column of }P_{t}C^{\top}$ for $j=1,...,m$, we obtain the two-parameter family with parameters $\gamma_{1},\gamma_{2}$ from [16].

III-D Action of gauge transformations on the class of linear feedback particle filters

The following two results summarize the details of how gauge transformations act upon the classes $\mathfrak{F}$ and its subclasses defined above.

Proposition III.2

Let $S\in\mathfrak{F}$ and $\tilde{S}_{t}=\mu_{t}+\tilde{E}_{t}$, where $\tilde{E}_{t}=g_{t}E_{t}$ with $g$ a deterministic gauge transformation according to Eqs. (13)-(14). Then $\tilde{S}$ also belongs to $\mathfrak{F}$ with $\tilde{\Omega}^{(i)}_{t}$, $i=0,1,...,m$ and $\tilde{H}^{(j)}_{t}$, given by

[TABLE]

Proof:

We have $d\tilde{E}_{t}=dg_{t}E_{t}+g_{t}dE_{t}+dg_{t}dE_{t}$. By using Eqs. (13) and (15), we obtain

[TABLE]

By rewriting $P_{t}^{-1}E_{t}=g_{t}^{\top}g_{t}^{-\top}P_{t}^{-1}g_{t}^{-1}\tilde{E}_{t}$, noting that $g_{t}^{-\top}P_{t}^{-1}g_{t}^{-1}=P_{t}^{-1}$, and then collecting all the terms, we obtain Eqs. (19) and (20) as well as $\tilde{G}_{t}=g_{t}G_{t}g_{t}^{\top}+M_{t}$. By using Eqs. (14) and (16), we obtain

[TABLE]

By rewriting $P_{t}^{-1}=g_{t}^{\top}P_{t}^{-1}g_{t}$ and then using Eq. (19) to substitute $g_{t}\Omega^{(i)}_{t}g_{t}^{\top}$, after cancelling all terms we obtain

[TABLE]

from which Eq. (III.2) follows. ∎

Corollary III.3

The deterministic gauge transformations act transitively on $\mathfrak{F}_{\text{det}}$, i.e. every deterministic gauge transformation of an $S\in\mathfrak{F}_{\text{det}}$ is in $\mathfrak{F}_{\text{det}}$, and for any pair of $S,\tilde{S}\in\mathfrak{F}_{\text{det}}$ there is a deterministic gauge transformation that maps $S$ to $\tilde{S}$.

IV Optimality criteria

We will now return to the multidimensional case. In order to select a particle filter among all linear feedback particle filters, additional criteria are required. Since every filter $S\in\mathfrak{F}$ is specified by a choice of $\Omega^{(i)}_{t}\in\text{Skew}(n)$ for $i=0,1,...,m$ and a vector $H^{(j)}_{t}\in\mathds{R}^{n}$ for $j=1,...,r$ and for all $t$, an optimality criterion can be formulated as a function of these quantities. Define

[TABLE]

where tr denotes the trace operator and $G$ is given by Eq. (16). This term appears to quadratic order in $dt$ when expanding $\text{tr}\,\mathds{E}[dE_{t}dE_{t}^{\top}|\mathscr{F}^{Y}_{t}]$ and can be associated with the transport cost as in [15]. It is therefore unsurprising that its minimization (under the constraint that the observation and Brownian motion terms are absent) yields the OTdetFPF from [15].

Proposition IV.1

The filter given by $\Omega^{(i)}_{t}=0$, $i=1,..,m$, $H^{(j)}=0$, $j=1,...,r$, and

[TABLE]

is identical to the optimal transport deterministic linear feedback particle filter (OTdetFPF).

Proof:

The OTdetFPF has dynamics555The reference [15] assumed $B=\mathds{1}$. For this reason, a few additional terms involving $B$ appear here.

[TABLE]

where $\hat{\Omega}_{t}$ is the unique solution of Eq. (11). This corresponds to Eq. (15) with $\Omega^{(i)}_{t}=0$, $i=1,..,m$, $H^{(j)}=0$, $j=1,...,r$, and

[TABLE]

Any critical point of the function $\mathscr{L}_{P,0}(\Omega,0,...,0)$ satisfies

[TABLE]

Since $(X,Y)\mapsto\text{tr}\,XY$ defines an inner product on $\text{Skew}(n)$, this implies $P^{-1}G^{\top}-GP^{-1}=0$. Substituting $G=\tfrac{1}{2}\dot{P}+\Omega$, we obtain the equation

[TABLE]

It can be checked that the unique solution is given by $\Omega=\tfrac{1}{2}(AP-PA^{\top}+\hat{\Omega})$. ∎

Whereas the optimal transport formulation in [15] always yields deterministic filters, the cost function (26) may be minimized without the constraints of Proposition IV.1. This yields a new class of filters that has not yet explicitly appeared in the literature.

Proposition IV.2

Let $S$ be a particle filter specified by data $\mathscr{D}_{t}=(\Omega^{(0)}_{t},\Omega^{(1)}_{t},...,\Omega^{(m)}_{t},H^{(1)}_{t},...,H^{(r)}_{t})$ consisting of $\Omega^{(i)}_{t}\in\text{Skew}(n)$ for $i=0,1,...,m$ and a vector $H^{(j)}_{t}\in\mathds{R}^{n}$ for $j=1,...,r$ for some $r\geq 0$. The following properties are equivalent:

1. i)

$\mathscr{L}_{P_{t},r}$* has a critical point at $\mathscr{D}_{t}$,* 2. ii)

$\mathscr{L}_{P_{t},r}$* has a global minimum at $\mathscr{D}_{t}$,* 3. iii)

$\mathscr{D}_{t}$* is given by $\Omega^{(0)}_{t}=0$ and*

[TABLE]

Proof:

By differentiating $\mathscr{L}_{P,r}$, we have i) if and only if

[TABLE]

for all $K^{(i)}\in\text{Skew}(n)$ and all $l^{(j)}\in\mathds{R}^{n}$, $i=0,...,m$ and $j=1,...,r$. Since tr is an inner product, this is equivalent to

[TABLE]

which in turn is equivalent to $G=0$ and hence to iii). This shows i)$\Leftrightarrow$iii) Moreover, $G=0$ when plugged into $\mathscr{L}_{P,r}$ gives a value of [math], which is the global minimum for this function. This shows that iii)$\Rightarrow$ii), which concludes the proof. ∎

Remark IV.3

For $n=1$, since all skew-symmetric $1\times 1$-matrices are zero, condition iii) in Proposition IV.2 reads $\sum_{j=1}^{r}(h^{(j)}_{t})^{2}=b^{2}+2aP_{t}-c^{2}P_{t}^{2}$. This has a solution if and only if $P_{0}\leq P_{t}\leq\frac{a+\sqrt{a^{2}+b^{2}c^{2}}}{c^{2}}$. Otherwise there are no global minimizers of $\mathscr{L}_{P_{t},r}$.

V Discussion

In this paper, we studied the intrinsic freedom in choosing dynamics of a particle filter if the only requirement for it is that the distribution of $S_{t}$ match the exact posterior distribution. We studied this for the example of the linear-Gaussian filtering problem where all calculations can be done explicitly, and the transformations can be assumed to be linear.

As noted in Section II-D, part of this freedom is clearly present in the nonlinear case in terms of the freedom of modifying any given pair of solutions for $K$ and $u$ by divergence-free vector fields. However, since the distribution is not known, such divergence-free fields cannot be easily found. Other forms of freedom, such as the exchange of noise terms for deterministic terms, that are easily accomplished in the linear case, also present more difficulties in the nonlinear case since it is not known, without having detailed knowledge about the distribution, how to compensate a given diffusion term by a deterministic one.

Generally speaking, in the nonlinear case both the gain estimation problem as well as the problem of finding gauge transformations are equally difficult because the distribution is unknown, and cannot be easily reduced to mean and covariance matrix such as in the linear case. However, the presence of gauge transformations suggests alternative formulation of the gain estimation problem. In particular, it raises the question whether the gradient solution for $K$, as it is used in the literature, is the most relevant solution.

The freedom in the choice of the gain can be viewed from yet another angle when considering the filtering problem for a hidden process with values in a smooth manifold $M$. If a feedback particle filter is to be designed on a manifold, there is no canonical riemannian metric and therefore no preferred gradient solution for $K$. Part of the gauge freedom of nonlinear feedback particle filters can therefore be attributed to the choice of a riemannian metric on the manifold. As explained in Section II-C, under certain conditions on the distribution (in the case of compact $M$, for all distributions that can be expressed as a smooth volume form) the choice of a riemannian metric on the state space gives rise to the gradient solution as a minimum-energy solution, which can also be interpreted as the solution to the following dynamical optimal transport problem for the 2-Wasserstein metric induced by the riemannian metric: given a direction in the space of probability distributions, prescribed by the filtering equation, what is the corresponding vector field in the state space such that the incremental transportation cost is minimized?

In future research, the connection between this optimal transport principle (for the gain function) and the optimal transport deterministic feedback particle filter should be further explored. Also, the broader significance of optimal transport principles and the choice of riemannian metric in nonlinear filtering deserves further study. In particular, desiderata such as numerical stability should be connected, if possible, to optimal transport principles in order to enable a principled design of feedback particle filters.

VI Acknowledgements

We thank Jean-Pascal Pfister, Anna Kutschireiter, Amir-hossein Taghvaei, and Prashant Mehta for helpful discussions, as well as the anonymous referees for their suggestions for improving the manuscript.