Optimal Linear Filtering for
Discrete-Time Systems with
Infinite-Dimensional Measurements

Infinite-dimensional representations are widely used in the study of distributed-parameter systems - see e.g. [4, 5, 6]. In such systems, the states and measurements are represented as functions both of time and a continuous-valued index taking values in some domain $\mathcal{D}$. The models may be stochastic [7, 8] or deterministic [4], in either discrete or continuous-time [9]. A critical high-level survey of the approaches, claims, and results in this area has been carried out in [5], and a broader treatment of the estimation and control of distributed-parameter systems can be found in [4, 6, 10, 11]. An overview of various observer design methodologies for linear distributed-parameter systems is also given in the survey [12]. Although this survey does not discuss stochastic systems, the references therein are classified as either early lumping models, where the system is discretized before observer design, and late lumping models, where the system is discretized after observer design. In all the works discussed in [12], both the state and measurement spaces are considered infinite-dimensional.

In contrast with the rich body of work above, in this article we examine systems with infinite-dimensional measurements but finite-dimensional states. complementary problem setup occurs in [13], although in that work the system is assumed to have an infinite-dimensional state space and a finite-dimensional measurement space, which sidesteps the issue of calculating the inverse of an infinite-dimensional object, an issue we will discuss further in Section 6.2. The work in [13] uses an orthogonal projection operator to reduce the state to a finite-dimensional subspace and proceeds to establish upper bounds on the resulting discretization error by way of sensitivity analysis of a Ricatti difference equation. A more common approach is the aforementioned early lumping method, where discretization is performed at the start to produce linear models with high-but-finite-dimensional states and measurements. This is the approach taken by [14], where a Kalman filter is designed for a finite-dimensional subset of the infinite-dimensional measurement, corrupted by the additive Gaussian noise that is spatially white. In this article, we do not perform any such lumping, early or late, and our analysis remains firmly in an infinite-dimensional setting.

See Table 1 for a concise summary of early derivations of linear distributed-parameter filters, as well as more recent work. The papers therein employ a number of different approaches, including Green’s function [15], orthogonal projections [7], least squares arguments [16], the Wiener-Hopf method [8], Bayesian approaches [17], and $\mathcal{H}_{2}$ estimation [18], to name a representative handful. The derivation in this article is closest to the orthogonal projection argument given in [7], although here we deal with a discrete-time system. It should also be noted that [7] assumes that the observation noise is both temporally and spatially white, whereas here we do not require the latter property. We do, however, introduce the assumption that the observation noise is stationary, which allows the novel derivation of a closed form expression of the optimal gain in Section 4.3. This form has not been presented in any of these previous works.

The model that we propose is relevant to systems where high-dimensional measurements, for instance from vision, radar, or lidar sensors, are used to estimate a low-dimensional state, such as the position and velocity of an autonomous agent or moving target. To the best of the authors’ knowledge, such systems have not been the focus of previous work in distributed parameter filtering.

Our key contribution is the derivation of an exact solution for the optimal linear filter for discrete-time systems with measurement noise and process noise modeled by random fields and random vectors respectively (Theorem 4.1, Procedure 1). The assumptions of finite state dimension, and that the measurement noise is stationary on domain $\mathcal{D}$, enable us to derive the optimal gain functional explicitly in terms of a multi-dimensional Fourier transform on the measurement index set $\mathcal{D}$.

This is in contrast to previous works in which the optimal gain is given in terms of a distributed-parameter inverse operator that is either not defined or defined implicitly [8], [17]. Moreover, this inverse can lead to implementation difficulties and is not always guaranteed to exist: these issues are discussed in further detail in Section 6.2. In comparison, our approach allows us to give sufficient conditions for the system to have a well-defined optimal gain function (Section 4): in rough terms, the effective bandwidth of the Fourier transform of the measurement kernel should be smaller than that of the noise covariance kernel.¹¹1The term bandwidth here refers to the effective range of “spatial” frequencies present in a measurement, not temporal frequencies.

An advantage of this continuum limit approach is that it enables us to make precise statements about the filters behaviour in the limit when measurement resolution approaches infinity. This paves the way for further analysis, as various finite-point approximation schemes may be formulated and the error compared with the asymptotic error of the optimal infinite-resolution filter in a principled and quantitative manner. Further discussion and analysis of these finite-point approximation schemes will be addressed in future work.

This article contains five significant extensions in comparison to our preliminary work in the conference paper [21].

The first major development is that this work, by way of a Hilbert space framework, presents necessary and sufficient conditions which guarantee the existence and uniqueness of the optimal gain function. Our previous work [21] employed a calculus of variations approach to yield only necessary conditions on the optimal gain function.

The second development is the extension of our proposed algorithm to a $d$-dimensional, not just scalar, measurement domain $\mathcal{D}$, allowing much richer sensor data to be modeled without obscuring the underlying spatial correlations. Examples include camera images ($d=2$), as well as high-resolution lidar and magnetic resonance imaging (MRI) ($d=3$). This extension entails the use of random fields and Fourier transforms with $d$-dimensional frequency, as opposed to standard random processes and Fourier transforms on a scalar axis. This has also allowed us to present simulation results that involve two-dimensional images, as opposed to the simpler scalar functions previously presented in [21].

The third development is a weakening of the assumptions on the observation noise. The filter provided in [21] was restricted to systems with Gaussian observation noise, whereas the formulation presented here loosens this assumption and instead only requires that the noise be stationary.

The fourth development is that this work clearly defines explicit conditions such that the presented arguments are rigorous, whereas the derivations in [21] were of a more preliminary and exploratory nature. This work provides a list of assumptions such that our derived filter is well-defined, and that each step in the derivation is valid.

The fifth and final development is that we explicitly impose stabilizability and detectability conditions on the system (Section 2.4). The detectability condition in particular involves the measurement noise covariance function, as well as the deterministic dynamical matrix and measurement function. This strongly contrasts with other works in infinite-dimensional systems theory, which impose conditions of strong or weak observability [22, 23] that are agnostic to noise.

These extensions yield a clearer picture of the relationship between the system and the optimal filter, and provide a more generalized solution for the optimal gain.

This article is organized as follows. Section 2 outlines the relevant definitions, theorems, and conditions relating to the Hilbert space of random variables and random fields, as well as the multi-dimensional Fourier transform and notions of system stabilizability and detectability. Section 3 gives our model definition and discusses some concrete examples which our model could be feasibly used to represent. This section also summarizes the key assumptions needed to ensure validity of our analysis, and the motivation behind each assumption. Section 4 derives necessary and sufficient conditions for the optimal gain of the proposed filter (4.1), before presenting a unique gain function which satisfies these conditions 4.2. Section 4.3 will present the full filter equations and stability of this filter will be discussed in Section 5. Section 6 will give a procedure to implement this filter, with Section 6.1 discussing computational complexity and Section 6.2 drawing comparisons to existing results. Finally, Section 7 presents a simulated system and analyzes the empirical performance of the filter compared to its expected performance.

In this work, we denote by $\mathcal{H}$, the Hilbert space of all real-valued random variables with finite second moments [24, Ch. 20]. The inner product of any two elements $u,v\in\mathcal{H}$ is defined by $\langle u,v\rangle=E[uv].$ An inner product operating on two random vectors $u,v$ in the Hilbert space $\mathcal{H}^{n}$ of random $n$-vectors is defined by $\langle u,v\rangle_{\mathcal{H}}=E[u^{\top}v]$. The derivation used in this article involves defining the state of our observed system as a vector residing in a Hilbert space and then selecting an error-minimizing estimate that resides in a subspace. The following theorem [25, p.51, Th. 2] is fundamental to this approach.

This article models measurement noise as an additive random field and some care is required when manipulating these objects. Rigorous and extensive analysis of these fields can be found in a number of textbooks [26], [27], [28]. A brief summary of relevant definitions is given below.

Let $(\Theta,\mathcal{F},P)$ be a probability space and $\mathbb{R}^{d}$ a Euclidian space endowed with Lebesgue measure. A random field is a function $v:\mathbb{R}^{d}\times\Theta\rightarrow\mathbb{R}^{m}$ which is measurable with respect to the product measure on $\mathbb{R}^{d}\times\Theta$.²²2For concision, the dependence on $\theta$ will usually be suppressed. A centered w.s.s. (wide-sense stationary) random field has the properties that i) $E[v(i)]=0$ for all $i$; ii) the random vector $v(i)\in\mathcal{H}$ for all $i$; and iii) the covariance function $R$ only depends on the difference between index points, so that for any $i_{1},i_{2}\in\mathbb{R}^{d},\ E[v(i_{1})v^{\top}(i_{2})]=R(i_{1}-i_{2})$. Such a field that is stationary over a domain $\mathcal{D}$ is often referred to as $\mathcal{D}-$stationary noise. In the special case where the covariance is also independent of direction, the noise is said to be isotropic. The stationarity condition arises in situations where the noise statistics possess translational invariance, in other words, the noise statistics look the same everywhere throughout the domain. Some sources of noise may also be locally stationary, in that they are closely approximated by stationary processes over finite windows. Spatial stationarity is a common assumption in a number of diverse contexts such as image denoising [3, 29, 30], radio signal propagation [31], and geostatistics [32, Ch. 74]. Figure 1 gives three example realizations of a stationary random field with a squared exponential covariance function and varying length scales.

A critical step in this article involves the transformation of an integral equation to an algebraic equation, facilitated by the multi-dimensional Fourier transform, which we now define.

This article will employ notions of stabilizability and detectability, which we define as follows:

In this section, we will prove the following key theorem:

Before we provide a proof of Theorem 4.1, we must first define some important sets, and projections onto those sets.
Consider the set of all possible elements of $\mathcal{H}$ which can be expressed as a bounded linear integral functional operating over all previous measurements:

	$\displaystyle M_{k-1}\triangleq\bigg{\{}m_{k-1}\in\mathcal{H}:m_{k-1}=\sum_{j=1}^{k-1}\int_{\mathbb{R}^{d}}\phi_{j}(i)z_{j}(i)di\bigg{\}},$			(8)
	$\displaystyle\text{ for some }\phi_{j}\in L_{1}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m})\cap L_{2}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m})$	.

Note that $M_{k-1}$ is a subspace of $\mathcal{H}$. We will also assume for the moment that $M_{k-1}$ is a closed subspace of $\mathcal{H}$. Suppose that we are already in possession of the optimal estimate of $x_{k}^{l}$ given the previous measurements. This estimate is the orthogonal projection of $x_{k}^{l}$ onto $M_{k-1}$,

$\displaystyle\hat{x}_{k|k-1}^{l}\triangleq\operatorname{proj}_{M_{k-1}}x_{k}^{l}.$

Given the measurements associated with $M_{k-1},\;\hat{x}_{k|k-1}^{l}$ is the best estimate of $x_{k}^{l}$ in the sense that it is the element of $M_{k-1}$ that is the unique $\mathcal{H}$-norm minimizing vector in relation to the $l^{th}$ component of the true state.

Let $z_{k}^{p}$ be the $p^{th}$ element of ${z}_{k}$ and $\hat{z}_{k}^{p}$ denote the orthogonal projection of $z_{k}^{p}$ onto $M_{k-1}$,

$\displaystyle\hat{z}_{k}^{p}(i)\triangleq\operatorname{proj}_{M_{k-1}}z_{k}^{p}(i)\quad\forall i.$

The measurement innovation vector is $s_{k}(i)\triangleq z_{k}(i)-\hat{z}_{k}(i)$, the $p^{th}$ entry of which is denoted $s^{p}_{k}(i)\in\mathbb{R}.$ Note that for all $m_{k-1}\in M_{k-1}$,

$\displaystyle s^{p}_{k}(i)\perp m_{k-1}\quad\forall i.$

(9)

and let

	$\displaystyle\mathcal{S}_{k}\triangleq\bigg{\{}\tilde{s}_{k}\in\mathcal{H}:\tilde{s}_{k}=\sum_{p=1}^{m}\int_{\mathbb{R}^{d}}\alpha_{k}^{p}(i)s_{k}(i)di\bigg{\}},$		(10)
	$\displaystyle\text{for some }\alpha_{k}^{p}\in L_{1}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m})\cap L_{2}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m}).$

which is a subspace of $\mathcal{H}.$ Here ${\alpha}_{k}\in L_{1}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{n\times m})\cap L_{2}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{n\times m}),$ with the $p^{th}$ row given by $\alpha_{k}^{p}$.

Having shown that the spaces $\mathcal{S}$ and $M_{k-1}$ are orthogonal, we will now find the projection of the true state $x_{k}$ onto the space of all possible elements in $\mathcal{H}$ that can be expressed by a bounded linear integral functional operating over all previous measurements and the current measurement. Let $M^{-}_{k}$ be the set of all possible outputs of linear functionals of any element in $\mathcal{Z}_{k}$,

	$\displaystyle M^{-}_{k}\triangleq\bigg{\{}m^{-}_{k}\in\mathcal{H}:m^{-}_{k}=\int_{\mathbb{R}^{d}}\phi(i)z_{k}(i)di\bigg{\}},$			(12)
	$\displaystyle\text{ for some }\phi\in L_{1}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m})\cap L_{2}(\mathbb{R}^{d},\mathbb{\mathbb{R}}^{1\times m}),\ z_{k}\in\mathcal{Z}_{k}$	.

Note that using this definition, it immediately follows that $M_{k-1}+M^{-}_{k}=M_{k}$, where the addition of sets is given by the Minkowski sum

$\displaystyle U+V=\{u+v|u\in U,v\in V\}.$

We also use the $\oplus$ operator to refer to a direct sum of subspaces, where any element within the direct sum is uniquely represented by a sum of elements within the relevant subspaces [44]. We wish to find the projection of $x_{k}^{l}$ onto the space of $M_{k-1}+M^{-}_{k}$, which is given by

	$\displaystyle\hat{x}^{l}_{k}$	$\displaystyle=\operatorname{proj}_{M_{k-1}+M^{-}_{k}}x_{k}^{l}$
		$\displaystyle=\operatorname{proj}_{M_{k-1}\oplus\mathcal{S}_{k}}x_{k}^{l}$
		$\displaystyle=\operatorname{proj}_{M_{k-1}}x_{k}^{l}+\operatorname{proj}_{\mathcal{S}_{k}}x_{k}^{l}$
		$\displaystyle=\hat{x}_{k|k-1}^{l}+\beta_{k}^{l},\quad\beta_{k}^{l}\triangleq\operatorname{proj}_{\mathcal{S}_{k}}x_{k}^{l}.$		(13)

This sum represents a decomposition of the projected state component given measurements up to time $k$ as a unique sum of a component within $M_{k-1}$ and a component orthogonal to $M_{k-1}.$ Note that as $\beta^{l}_{k}\in{\mathcal{S}_{k}}$, it can be formulated as $\sum_{p=1}^{m}\int_{\mathbb{R}^{d}}\kappa_{k}^{l,p}(i)s_{k}^{p}(i)di$, each $\kappa_{k}^{l,p}(i)$ being the $l^{th}$ row, $p^{th}$ column element of some ${\kappa}_{k}(i)\in\mathbb{R}^{n\times m}$ which must be determined. This matrix-valued function ${{\kappa}_{k}}(i)$ is equivalent to our optimal gain. We are now ready to prove Theorem 4.1, which we will do by formulating necessary and sufficient conditions on the gain such that the optimal post-measurement correction term $\beta_{k}=\int_{\mathbb{R}^{d}}\kappa_{k}(i)s_{k}(i)\,di$ is generated. The proof of Theorem 4.1 is as follows:

An alternative derivation is given in [21] based on the Calculus of Variations approach. The derivation above, which uses the Hilbert Projection Theorem, has the advantage of a guaranteed uniqueness of the correction term $\int_{\mathbb{R}^{d}}\kappa_{k}(i)s_{k}(i)di$ and state estimate $\hat{x}_{k}$, which follows immediately from Theorem 2.1. This approach avoids the use of the functional derivative which has attracted some criticism [8], [5], as well as closely mirroring the methodology of the original Kalman proof [45].

Recall that we initially assumed $M_{k-1}$ is a closed subspace. We will now show that an explicit closed-form solution of the optimal gain function can be derived via $d$-dimensional Fourier based methods, thanks to the w.s.s. assumption of the measurement noise. We will then show that when a unique gain function which satisfies (17) exists, it is indeed the case that $M_{k-1}$ is closed, which is derived in Corollary 4.5.

We may now proceed to find an explicit value for the optimal gain function that satisfies (7). Consider the matrix-valued functional

$\displaystyle F(\kappa_{k})\triangleq\bigg{(}I-\int_{\mathbb{R}^{d}}\kappa_{k}(i)\gamma(i)di\bigg{)}P_{k|k-1}\in\mathbb{R}^{n\times n}$

(18)

and note that,

$\displaystyle\int_{\mathbb{R}^{d}}\kappa_{k}(i)R(i^{\prime}-i)di$

$\displaystyle=F(\kappa_{k})\gamma^{\top}(i^{\prime})$

(19)

The left hand side of (19) is simply a multi-dimensional convolution. Many useful properties of the scalar Fourier transform carry over to the multi-dimensional Fourier transform. We will make use of the convolution property [46],

		$\displaystyle\mathcal{F}\left\{\int_{\mathbb{R}^{d}}f(i_{1}-y_{1},...,i_{d}-y_{d})g(y_{1},...,y_{d})\right\}dy_{1}...dy_{d}$
	$\displaystyle=$	$\displaystyle\mathcal{F}\{(f\ast g)(i_{1},i_{2},...,i_{d})\}=(\mathcal{F}\{f\})(\omega)(\mathcal{F}\{g\})(\omega).$		(20)

This enables us to find that (19) is equivalent to

$\displaystyle(\kappa_{k}\ast R)(i^{\prime})$

$\displaystyle=F(\kappa_{k})\gamma^{\top}(i^{\prime})$

(21)

Note that $\omega$ will be a vector of the same dimension as $i$. We now require three assumptions to ensure existence and invertibility of the Fourier transform of $R$, as well as the existence of the Fourier transform of $\gamma^{\top}$: See A1 See A2 See A3 Applying (21) and (20) to (19) implies

	$\displaystyle\bar{\kappa}_{k}(\omega)\bar{R}(\omega)$	$\displaystyle=F(\kappa_{k})\bar{\gamma}^{\top}(\omega)$
	$\displaystyle\bar{\kappa}_{k}(\omega)$	$\displaystyle=F(\kappa_{k})\bar{\gamma}^{\top}(\omega)\bar{R}(\omega)^{-1}$
	$\displaystyle\kappa_{k}(i)$	$\displaystyle=F(\kappa_{k})f(i),\ f(i)=\mathcal{F}^{-1}\{\bar{\gamma}^{\top}(\omega)\bar{R}(\omega)^{-1}\}.$		(22)

For this method to be applicable we must also make the following assumption See A4 Combining (18) and (22) yields

	$\displaystyle F(\kappa_{k})$	$\displaystyle=\bigg{(}I-F(\kappa_{k})\int_{\mathbb{R}^{d}}f(i)\gamma(i)di\bigg{)}P_{k|k-1}$
	$\displaystyle\implies F(\kappa_{k})$	$\displaystyle=P_{k|k-1}\bigg{(}I+\int_{\mathbb{R}^{d}}f(i)\gamma(i)diP_{k|k-1}\bigg{)}^{-1}$		(23)
	$\displaystyle\implies\kappa_{k}(i)$	$\displaystyle=P_{k|k-1}\bigg{(}I+\int_{\mathbb{R}^{d}}f(i)\gamma(i)diP_{k|k-1}\bigg{)}^{-1}f(i)$

Note that this directly implies that for all $k\in\mathbb{N}$, $\kappa_{k}\in L_{1}(\mathbb{R}^{d},\mathbb{R}^{n\times m})\cap L_{2}(\mathbb{R}^{d},\mathbb{R}^{n\times m})$. For simplicity of notation, recall from (6) that $S=\int_{\mathbb{R}^{d}}f(i)\gamma(i)di\in\mathbb{R}^{n\times n}.$ Having thus found the optimal $\kappa_{k}(i)$ function which generates $\beta_{k}$, we now modify (13) to find for any state component $l$,

	$\displaystyle\operatorname{proj}_{M_{k-1}+M^{-}_{k}}x_{k}^{l}$
	$\displaystyle=\hat{x}_{k|k-1}^{l}+\sum_{p=1}^{m}\int_{\mathbb{R}^{d}}\big{[}P_{k|k-1}(I+SP_{k|k-1})^{-1}f(i)\big{]}^{l,p}s^{p}_{k}(i)di.$

These $l$ equations may be stacked up into a single matrix equation of the form

	$\displaystyle\operatorname{proj}_{M_{k-1}+M^{-}_{k}}x_{k}$
	$\displaystyle=\hat{x}_{k|k-1}+\int_{\mathbb{R}^{d}}P_{k|k-1}(I+SP_{k|k-1})^{-1}f(i)s_{k}(i)di,$		(24)

where the projection operator acting on the state vector is understood to represent the vector of each components projections. We have thus far assumed that $\hat{x}_{k|k-1}$ is known to us, we will complete the chain now by determining this value explicitly as

	$\displaystyle\hat{x}_{k|k-1}$	$\displaystyle=\operatorname{proj}_{M_{k-1}}x_{k}$
		$\displaystyle=\operatorname{proj}_{M_{k-1}}Ax_{k-1}+\operatorname{proj}_{M_{k-1}}w_{k-1}$
		$\displaystyle=A\operatorname{proj}_{M_{k-1}}x_{k-1}$
		$\displaystyle=A\operatorname{proj}_{M_{k-2}+M^{+}_{k-2}}x_{k-1}$
		$\displaystyle=A\hat{x}_{k-1}$		(25)

The form of (24) with this substitution is then

	$\displaystyle\hat{x}_{k}=$	$\displaystyle\operatorname{proj}_{M_{k-1}+M^{-}_{k}}x_{k}$
	$\displaystyle=$	$\displaystyle A{\hat{x}_{k-1}}+P_{k|k-1}(I+SP_{k|k-1})^{-1}\int_{\mathbb{R}^{d}}f(i)s_{k}(i)di$
	$\displaystyle=$	$\displaystyle A{\hat{x}_{k-1}}+$
		$\displaystyle P_{k|k-1}(I+SP_{k|k-1})^{-1}\int_{\mathbb{R}^{d}}\mkern-5.0muf(i)\big{(}z_{k}(i)-\gamma(i)\hat{x}_{k|k-1}\big{)}di,$

where in the last line we substitute $s_{k}(i)=z_{k}(i)-\hat{z}_{k}(i)$, noting that

	$\displaystyle\hat{z}_{k}(i)$	$\displaystyle=\operatorname{proj}_{M_{k-1}}z_{k}(i)$
		$\displaystyle=\operatorname{proj}_{M_{k-1}}\big{(}\gamma(i)x_{k}+v_{k}(i)\big{)}$
		$\displaystyle=\gamma(i)\operatorname{proj}_{M_{k-1}}x_{k}$
		$\displaystyle=\gamma(i)\hat{x}_{k|k-1}.$