Recovering Latent Signals from a Mixture of Measurements using a Gaussian Process Prior

Gaussian processes (GPs) [5] are a nonparametric prior distribution on functions $f:X\mapsto\mathbb{R}$, for any set $X$ where a covariance function can be defined (e.g., a metric space), such that for any finite collection of inputs $\mathbf{x}=[x_{1},x_{2},\ldots,x_{n}]^{\top}\in X^{n}$, the corresponding outputs $f(x_{1}),f(x_{2}),\ldots,f(x_{n})\in\mathbb{R}$ are jointly Gaussian, that is,

where the mean and covariance functions, $\mu(\cdot)$ and $K(\cdot,\cdot)$ respectively, are parametric forms that determine the spatial properties of samples drawn from the GP prior. Training the GP to observed data involves finding the parameters of the mean and covariance functions, then, the posterior distribution of the entire function $f$ conditional to a set of observations is Gaussian.

The GP framework is well suited to the sensing setting in fig. 1, since modelling latent signal as a GP results in the posterior distribution of the latent process (given the mixture of measurements) being Gaussian as well. Previous GP-based models for convolution processes [6, 7, 8, 9] model signals a convolution between a continuous-time filter and a white-noise process, which is unsuitable to represent the latent process in the sensing application where the spatial correlation of the process is fundamental. Conversely, [10] allows GPs as latent functions but to address the multi-output case, where the aim is to perform inference on the outputs rather than the latent processes. Furthermore, these methods consider continuous-time convolution filters, which is computationally demanding and requires, e.g., variational approximations [11, 12]. Consequently, closed-form and computationally-efficient Bayesian reconstruction of the latent process is still an open problem in sensing applications.

Our approach will consider the weights $w_{i,j}$ to be unknown and learn them from the observations, however, let us assume they are known in this section in order to address the problem as a linear system; for ease of notation, we assume there are no common locations across the measurements (i.e., $i\neq i^{\prime},j\neq j^{\prime}\Rightarrow x_{i,j}\neq x_{i^{\prime},j^{\prime}}$). With these assumptions, eq. (4) is an underdetermined linear system: there are $L=MN$ unknowns and $N$ equations, where $N\ll L$, in fact, the general solution to such a system has $L-N$ free parameters or degrees of freedom (neglecting the inconsistent case).

This underdetermined system has infinite solutions, with the minimum-norm solution given by $\hat{\mathbf{m}}=\mathbf{W}^{+}\mathbf{y}$, where $\mathbf{W}^{+}$ is the Moore-Penrose pseudoinverse of $\mathbf{W}$ [13]. Using this solution to recover the latent signal has a number of drawbacks: (i) it requires the weights $w_{i,j}$ to be known, (ii) it only recovers the process at the measured locations without providing any insight about regions not measured, and (iii) it does not provide a measure of uncertainty for the estimates, e.g., in the form of error bars.

As the system in eq. (4) has infinite solutions, a constraint (or regularisation criterion) has to be imposed to reduce the number of solutions, or critically, to find a single solution. The Moore-Penrose pseudoinverse imposes the minimum-norm criterion, we instead assume a probabilistic condition (a prior) on the spatial structure of the solution. Specifically, we (i) place a GP prior on the latent signal to then (ii) find the posterior distribution of the latent signal conditional to the mixture of measurements, even at locations that were not measured. A key property of this approach is that these two steps are performed analytically, since the latent signal and the mixture of measurements are jointly Gaussian. We next present a formal description of the proposed generative model.

We model the latent process $f(\cdot)$ in (2) as a GP over the set of locations $X$ given by

where $\mu_{f}:X\mapsto\mathbb{R}$ and $K_{f}:X\times X\mapsto\mathbb{R}$ are the GP mean and covariance functions respectively. As the linear combination of jointly-Gaussian random variables (RVs) is Gaussian, the observations in eq. (3) are Gaussian RVs indexed by $\mathbf{x}_{i}=[x_{i,1},\ldots,x_{i,M}]^{\top}\in X^{M}$ with mean and covariance respectively given by

where $\sigma^{2}_{\eta},\sigma^{2}_{\epsilon}$ are the variances of the measurement and observation noises respectively. Additionally, note that if measurements are always taken at different locations (which is the case if the set of locations is continuous) we have $\delta(x_{i,j}-x_{i^{\prime},j^{\prime}})=0\ \forall i,i^{\prime},j,j^{\prime}$ and

where ${K_{f}}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})\in\mathbb{R}^{M\times M}$ is the Gram matrix evaluated on $\mathbf{x}_{i}$ and $\mathbf{x}_{i^{\prime}}$, $\mathbf{x}_{i}=[x_{i,1},\ldots,x_{i,M}]^{\top}$ and $\mathbf{w}_{i}=[w_{i,1},\ldots,w_{i,M}]^{\top}$.

The covariance of the observations is a mixture of evaluations of the covariance function of the process $f(\cdot)$ at the locations measured. This implies that a single entry in common between measurement locations $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ is enough for the observations $y_{i}(\mathbf{x}_{i})$ and $y_{j}(\mathbf{x}_{j})$ to be correlated. This mixture-of-kernel structure resembles additive GPs [14] and multikernel learning [15, 16, 17], however, these methods combine different kernels (evaluated on a common input) for expressive kernel design, whereas the proposed model combines evaluations of a single kernel on different locations to find the spatial structure of the latent signal. We refer to the presented model as Gaussian processes mixture of measurements (GPMM) and give a graphical model illustration in fig. 2.

The posterior density of $f$ given the observations $\mathcal{D}$, $p(f|\mathcal{D})$, is Gaussian and determined by (i) the prior mean—assumed to be zero in this case, (ii) the autocovariance of the mixture-of-measurement process $y$—given in eq. (7), and (iii) the covariance between the latent process $f(x)$ and the observation at $\mathbf{x}_{i}=[x_{i,1},\ldots,x_{i,M}]^{\top}$ given by $y(\mathbf{x}_{i})=y_{i}=\sum_{j=1}^{M}w_{i,j}m_{i,j}$, this covariance is

Denoting $\mathbf{y}=[y(\mathbf{x}_{1}),\ldots,y(\mathbf{x}_{N})]^{\top}$, the predictive posterior is:

thus, the reconstruction of the latent process $f$ can be computed in closed-form (we used the notation $K_{fy}^{\top}=K_{yf}$).

Without loss of generality, let us consider that observations $y_{i}$ were taken at a grid $\mathbf{x}_{\text{grid}}=[x_{1},\ldots,x_{M}]$ were only a few weights are nonzero per observation. Furthermore, without evidence for spatial correlation of the latent process $f(\cdot)$, its covariance matrix is the identity multiplied by $\sigma_{f}^{2}$ (signal power)—i.e., $\mathbf{K}(\mathbf{x}_{\text{grid}},\mathbf{x}_{\text{grid}})=\sigma_{f}^{2}\mathbf{I}_{M}$. Consequently, denoting the mixing weights by $\mathbf{W}=[\mathbf{w}_{1},\ldots,\mathbf{w}_{N}]$ and $\mathbf{x}$ the input of the observed process, from eqs. (9)-(10) the covariances are given by

Finally, introducing the above two expressions in the posterior mean in eq. (11) gives the solution to the linear system

The connection between solutions to linear systems and GP models is therefore established: when the noise variances are negligible w.r.t. the signal power $\left(\sigma_{f}^{2}\gg\sigma_{\epsilon}^{2},\sigma_{\eta}^{2}\right)$ the ratios in eq. (15) ${\sigma_{f}^{2}}/{(\sigma_{\epsilon}^{2}+\sigma_{f}^{2})}\rightarrow 1$ and ${\sigma_{\eta}^{2}}/{(\sigma_{\epsilon}^{2}+\sigma_{f}^{2})}\rightarrow 0$, and the Moore-Penrose inverse is obtained. On the contrary, when the noise variance $\sigma_{\epsilon}^{2}$ is large the estimate decays to zero, or reverts to the prior, since the measurements are not reliable. Furthermore, when the observation noise $\sigma^{2}_{\eta}=0$, eq. (15) is equivalent to the ordinary least squares and when $\sigma^{2}_{\eta}>0$ to the regularised least squares (ridge regression).

Finally, we emphasise that, unlike the Moore-Penrose method, the proposed GPMM jointly infers the mixing weights and the complete latent process, while also providing a measure of uncertainty given by the variance in eq. (12).

This first toy example considered a draw from a GP with SE covariance as the latent function and 120 observations with 7 measurements per observation. Fig. 3 shows the GP estimates and their mean square error (MSE) for both the standard GP with SE kernel (termed GP-SE) and the proposed GPMM with SE kernel (termed GPMM-SE). Notice how GP-SE fails to recover all the extrema of the latent process and adjust very tightly to the observations, this is because the convolution performed by the sensor smooths out the extrema in the observations which are then trusted as true values by GP-SE. Conversely, the GPMM-SE was able to recover all the extrema, place appropriate error bars on the latent function and report a lower estimation error.

Instantaneous-frequency estimation is performed averaging over a time window thus motivating the use of the proposed GPMM. We used a heart-rate time-series from the MIT-BIH Database [18] (ecg.mit.edu/time-series) and constructed a lower-resolution version of it composed by a mixture of measurements. Fig. 4 shows the recovery of a heart-rate signal from such observations for both the GP-SE and the proposed GPMM-SE, where the GPMM-SE again outperforms the GP-SE in terms of estimation MSE and uncertainty representation. We used 240 observations with 7 measurements each.

Motivated by edge-detection applications, we then considered the Heaviside step function. To cater for discontinuous signals, we used the neural network (NN) kernel [19] and implemented both the standard GP-NN and the proposed GPMM-NN to recover the latent step function. Observe in 5, how GPMM-NN successfully recovered the discontinuity with low variance. We used 120 observations with 7 measurements each.

Our aim is to recover spatial structure in the latent process, in this sense, Fig. 6 shows the PSDs for the smooth, heart-rate, and step function from top to bottom, observe how the GPMM (blue) was able to better recover the spectrum in all the experiments considered unlike the standard GP (red), where the latter failed to identify the spectral content of the latent process that was removed by the convolution sensor.