Sparse Bayesian Learning with Dynamic Filtering for Inference of Time-Varying Sparse Signals

TL;DR

This paper introduces SBL-DF, a robust and efficient dynamic filtering algorithm for time-varying sparse signals that outperforms existing methods, especially in challenging dictionary conditions.

Contribution

The paper proposes a novel SBL-DF algorithm that integrates dynamics models into the hyperpriors, improving robustness and efficiency in estimating time-varying sparse signals.

Findings

SBL-DF converges faster and more accurately than standard SBL.

The algorithm is robust to imperfect dynamics models and dictionary coherence.

SBL-DF outperforms state-of-the-art methods in challenging scenarios.

Abstract

Many signal processing applications require estimation of time-varying sparse signals, potentially with the knowledge of an imperfect dynamics model. In this paper, we propose an algorithm for dynamic filtering of time-varying sparse signals based on the sparse Bayesian learning (SBL) framework. The key idea underlying the algorithm, termed SBL-DF, is the incorporation of a signal prediction generated from a dynamics model and estimates of previous time steps into the hyperpriors of the SBL probability model. The proposed algorithm is robust to imperfect dynamics models (due to the propagation of dynamics information through higher-order statistics), robust to certain undesirable dictionary properties such as coherence (due to properties of the SBL framework), allows the use of arbitrary dynamics models, and requires the tuning of fewer parameters than many other dynamic filtering algorithms do. We also extend the fast marginal likelihood SBL inference procedure to the informative hyperprior setting to create a particularly efficient version of the SBL-DF algorithm. Numerical simulations show that SBL-DF converges much faster and to more accurate solutions than standard SBL and other dynamical filtering algorithms. In particular, we show that SBL-DF outperforms state of the art algorithms when the dictionary contains challenging structure such as coherence and column scaling, as is found in many practical applications.

Equations

$${\mathbf{x}}^{(t)}=f_t\left({\mathbf{x}}^{(t-1)}\right)+{\mathbf{n}}^{(t)},$$

$${\mathbf{y}}^{(t)}=\mathbf{\Phi }{\mathbf{x}}^{(t)}+{\mathbf{e}}^{(t)}$$

$$p(\mathbf{y}|\mathbf{x},\lambda )=\mathcal{N}(\mathbf{\Phi }\mathbf{x},\lambda \mathbf{I}).$$

$$p(x_i|{\gamma }_i)={\mathcal{N}}_{x_i}\left(0,{\gamma }_i\right),$$

$$p({\gamma }_i)={\mathcal{I}\mathcal{G}}_{{\gamma }_i}(a_i,b_i)=\frac{b_i^{a_i}}{\Gamma (a_i)}{\gamma }_i^{-a_i-1}{e}^{-b_i{\gamma }_i^{-1}}$$

$$p(\lambda )={\mathcal{I}\mathcal{G}}_{\lambda }(c,d)=\frac{d^c}{\Gamma (c)}{\lambda }^{-c-1}{e}^{-d{\lambda }^{-1}},$$

$$p(\mathbf{x},\mathbf{\gamma },\lambda |\mathbf{y})=p(\mathbf{x}|\mathbf{y},\mathbf{\gamma },\lambda )p(\mathbf{\gamma },\lambda |\mathbf{y}),$$

$$\mathbf{\Sigma }={\left({\mathbf{\Gamma }}^{-1}+{\lambda }^{-1}{\mathbf{\Phi }}^T\mathbf{\Phi }\right)}^{-1}\text{and}\mathbf{\mu }={\lambda }^{-1}\mathbf{\Sigma }{\mathbf{\Phi }}^T\mathbf{y}.$$

$$\mathbf{C}=\lambda \mathbf{I}+\mathbf{\Phi }\mathbf{\Gamma }{\mathbf{\Phi }}^T.$$

$$\ell (\mathbf{\gamma },\lambda )=-\log \left[p(\mathbf{y}|\log \mathbf{\gamma },\log \lambda )\prod _{i=1}^{N}p(\log {\gamma }_i)p(\log \lambda )\right]$$

$$\propto \log \left|\mathbf{C}\right|+{\mathbf{y}}^T{\mathbf{C}}^{-1}\mathbf{y}-2\sum _{i=1}^N\left(a_i\log {\gamma }_i^{-1}-b_i{\gamma }_i^{-1}\right).$$

$${\gamma }_i^{(k+1)}=\frac{{\left(\mathbf{\Sigma }\right)}_{ii}+{\mathbf{\mu }}_i^2+2b_i}{1+2a_i}$$

$${\lambda }^{(k+1)}=\frac{{\Vert \mathbf{y}-\mathbf{\Phi }\mathbf{\mu }\Vert }_2^2+\mathrm{Tr}\left[{\mathbf{\Phi }}^T\mathbf{\Phi }\mathbf{\Sigma }\right]}{M}.$$

$${\ell }_{\mathrm{dyn}}\left(\mathbf{\gamma },\lambda \right)=2\sum _i\left(-a_i\log {\gamma }_i^{-1}+b_i{\gamma }_i^{-1}\right)$$

$${\mathbf{\gamma }}_{\mathrm{dyn}}=\underset{\mathbf{\gamma }}{\mathrm{argmin}}\mathbb{E}\left[{\Vert \hat{\mathbf{x}}-\tilde{\mathbf{x}}\Vert }_2^2\right].$$

$$\frac{\partial }{\partial {\mathbf{\gamma }}^{-1}}$$

$$=$$

$${\gamma }_{\mathrm{dyn},i}={\tilde{x}}_i^2$$

$${\gamma }_i^{\ast }=\frac{b_i}{a_i}.$$

$$a_i=\xi \text{and}{b}_i=\xi {\tilde{x}}_i^2,$$

$$p(x_i|a_i,b_i)$$

$$=\frac{b_i^{a_i}\Gamma (a_i+\frac{1}{2})}{\sqrt{2\pi }\Gamma (a_i)}{\left(b_i+\frac{x_i^2}{2}\right)}^{-\left(a_i+\frac{1}{2}\right)}$$

$$z_i^{(k+1)}$$

$${\gamma }_i^{(k+1)}$$

$${\mathbf{x}}^{(k+1)}$$

$$+2\lambda \sum _{i\in \mathcal{T}}{\left(z_i^{(k+1)}\right)}^{1/2}\sqrt{x_i^2+2b_i}],$$

$$z_i^{(k+1)}$$

$$\mathbf{C}=\lambda \mathbf{I}+\mathbf{\Phi }\mathbf{\Gamma }{\mathbf{\Phi }}^T=\lambda \mathbf{I}+\sum _{i=1}^N{\gamma }_i{\mathbf{\varphi }}_i{\mathbf{\varphi }}_i^T.$$

$$\ell ({\mathbf{\gamma }}_{-j},\lambda )$$

$$-2\sum _{i=1,i\ne j}^N\left(a_i\log {\gamma }_i^{-1}-b_i{\gamma }_i^{-1}\right)$$

$$\ell ({\gamma }_j)$$

$$-\left(2a_j+1\right)\log {\gamma }_j^{-1}+2b_j{\gamma }_j^{-1},$$

$$\frac{\partial \ell ({\gamma }_j)}{\partial {\gamma }_j}$$

Full text

††M. R. O’Shaughnessy, M. A. Davenport, and C. J. Rozell are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332 USA (e-mail: moshaughnessy6@gatech.edu, mdav@gatech.edu, crozell@gatech.edu). This work was supported in part by NSF grants CCF-1350616 and CCF-1409422, James S. McDonnell Foundation grant number 220020399, a gift from the Alfred P. Sloan Foundation, and the National Defense Science and Engineering Graduate (NDSEG) Fellowship.

Note: Single column version of manuscript to appear in IEEE Transactions on Signal Processing111https://doi.org/10.1109/TSP.2019.2961229.

1 Introduction

Many signal processing applications require the reconstruction of sparse and time-varying signals from undersampled and noisy measurements. For example, in medical imaging it is well-known that MRI and CT images can be accurately represented in a wavelet basis by only a few nonzero coefficients that vary slowly over time due to organ movement [1, 2]. In direction of arrival problems, targets may be located in only a few spatial locations, and their movement may be governed by a known model [3]. The signals in each of these applications are both sparse (i.e., can be represented using only a few elements in some basis or dictionary) and vary over time, perhaps according to some imperfectly known dynamics model.

When recovering static sparse signals, incorporating a signal model that exploits this low-dimensional structure can result in dramatic improvements in estimation accuracy, allowing sparse signals to be reconstructed from extremely underdetermined systems [4, 5, 6]. In dynamic settings, we would like to be able to additionally exploit the a priori knowledge that either the signal changes only slowly over time, or according to a known (but potentially imperfect or noisy) dynamics model. Many effective algorithms for improving inference using dynamics information have been developed; a smaller literature has developed algorithms that exploit both dynamics information and the low-dimensional structure of signals. These algorithms are reviewed in detail in Section 2.2.

Most of these methods for exploiting both sparsity and dynamics information do so by extending variants of $\ell_{1}$-minimization algorithms. In this paper, however, we consider a sparsity-based dynamic filtering algorithm that makes use of a sparse recovery procedure based on a hierarchical probabilistic model called sparse Bayesian learning (SBL) [7, 8]. We do so for several reasons. First, in many cases the sparse Bayesian learning algorithm has been shown to produce more accurate estimates than $\ell_{1}$-based sparse recovery methods [9, 10] [11, Sec. I-C]. Second, the SBL procedure has been shown to be robust to undesirable dictionary structure such as the combination of high coherence222We say that a dictionary $\bm{\Phi}$ has high coherence if the normalized inner product of its columns, $\underset{i,j}{\max}\frac{\left|\bm{\Phi}_{i}^{T}\bm{\Phi}_{j}\right|}{\lVert\bm{\Phi}_{i}\rVert\lVert\bm{\Phi}_{j}\rVert}$, is large. Intuitively, this makes reconstruction challenging because it is difficult to differentiate measurements of $x_{i}$ and $x_{j}$, particularly in the presence of measurement noise [12]. and diverse column magnitudes [13], structure that violates sufficient conditions for $\ell_{1}$-based sparse estimation algorithms but that nonetheless commonly appears in real-world applications such as imaging [14]. Third, many applications make use of the SBL framework because its probabilistic framework allows other application-specific priors or constraints to be incorporated into the inference process in a simple and interpretable way (e.g., [15]). Finally, the SBL inference procedure allows the noise variance (sparsity) parameter to be learned automatically, reducing the need for manual parameter tuning compared with similar algorithms.

In this paper we develop an SBL-based algorithm for the online recovery of sparse time-varying signals. Our SBL with dynamic filtering (SBL-DF) algorithm replaces the uninformative hyperpriors typically used by the SBL probability model with informative hyperpriors that have parameters set using a dynamics-based signal prediction. This strategy allows our algorithm to incorporate dynamics information in a manner that is particularly robust to noise or modeling error, utilize arbitrary dynamics models, and inherit the beneficial properties of the SBL framework (particularly its robustness to adversely structured dictionaries).

Mathematically, we consider the problem of recovering a signal $\bm{x}^{(t)}\in\mathbb{R}^{N}$ at times $t=1,\dots,L$ that evolves according to the Markov process

[TABLE]

where $\bm{x}^{(t)}$ is sparse (i.e., $\left\|\bm{x}^{(t)}\right\|_{0}\ll N$), $f_{t}~{}\colon~{}\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is the known dynamics model that describes the evolution of the signal from time $t-1$ to time $t$, and $\bm{n}^{(t)}$ represents the error in the dynamics model (referred to as innovations). Rather than observe the signal directly, we obtain noisy underdetermined measurements $\bm{y}^{(t)}\in\mathbb{R}^{M}$ through the dictionary $\bm{\Phi}\in\mathbb{R}^{M\times N}$ (where typically $M\ll N$) according to the model333Although we consider $\bm{x}$ to be sparse in the canonical domain for clarity of notation, our approach can easily be extended to the case where $\bm{x}$ is sparse in some other basis or dictionary. Similarly, although we consider $\bm{\Phi}$ to be fixed for all $t$ for simplicity, it can easily be extended to be different for each time $t$.

[TABLE]

where $\bm{e}^{(t)}\in\mathbb{R}^{M}$ represents the measurement noise at time $t$.

After reviewing the SBL procedure and existing methods for dynamic filtering of sparse signals in Section 2, we present our main contributions in Section 3: the SBL-DF algorithm in Section 3.1, an investigation of how some properties of static SBL extend to SBL-DF in Section 3.2, and a fast version of the algorithm using an extension of the fast marginal likelihood procedure [16] in Section 3.3. In Section 4, we provide experimental validation of the algorithm and explore how its performance compares to state of the art algorithms in various settings. We conclude in Section 5 with a brief discussion of avenues for future work.

2 Background and Related Work

2.1 Sparse Bayesian learning (SBL)

2.1.1 Probability model

The SBL model makes the usual assumption that measurements are corrupted by i.i.d. Gaussian measurement noise, $\bm{e}\sim\mathcal{N}(0,\lambda\bm{I})$, yielding the likelihood

[TABLE]

A zero-mean Gaussian prior is placed on each element of $\bm{x}$,

[TABLE]

with variances parameterized by hyperparameters $\bm{\gamma}\in\mathbb{R}^{N}$.

This Gaussian prior does not have the high kurtosis of distributions known to encourage sparsity. Instead, the SBL procedure encourages sparse solutions by choosing appropriate values for $\bm{\gamma}$ based on the observed data $\bm{y}$. Intuitively, if a small variance $\gamma_{i}$ is chosen, the probability mass of the prior $p(x_{i}|\gamma_{i})$ is concentrated near zero, so significant evidence from the measurements will be required to make $x_{i}$ nonzero. Conversely, a large variance $\gamma_{i}$ makes the prior very wide, giving the model more flexibility to select nonzero $x_{i}$. The specific method for accomplishing this automatic variable selection is described in the next section.

To complete the probability model, a (conjugate) inverse gamma hyperprior is placed on each variance $\gamma_{i}$ and $\lambda$,

[TABLE]

and

[TABLE]

where $\left\{a_{i}\right\}$ and $c$ are shape parameters, $\left\{b_{i}\right\}$ and $d$ are scale parameters, and $\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx$. In practice, the parameters of these inverse gamma hyperpriors are typically selected to make the hyperpriors either uninformative (i.e., flat) by setting $a_{i}=b_{i}=c=d=0~{}\forall~{}i$ or scale-invariant (i.e., Jeffreys) by setting $a_{i}=c=1,~{}b_{i}=d=0~{}\forall~{}i$ [7], [17, Sec. 2.2.3].

2.1.2 Inference

Computing the estimate $\widehat{\bm{x}}_{\mathrm{SBL}}$ from the SBL probability model requires the calculation of the posterior $p(\bm{x},\bm{\gamma},\lambda|\bm{y})$. Because this quantity cannot be computed analytically, the SBL procedure proceeds by decomposing the full posterior as

[TABLE]

where we call $p(\bm{x}|\bm{y},\bm{\gamma},\lambda)$ the source posterior and $p(\bm{\gamma},\lambda|\bm{y})$ the hyperparameter posterior.

The source posterior $p(\bm{x}|\bm{y},\bm{\gamma},\lambda)$ admits the closed form expression $p(\bm{x}|\bm{y},\bm{\gamma},\lambda)=\mathcal{N}(\bm{\mu},\bm{\Sigma})$, where, denoting $\bm{\Gamma}=\mathrm{diag}(\bm{\gamma})$,

[TABLE]

This expression allows us to calculate the point estimate $\widehat{\bm{x}}_{\mathrm{SBL}}=\mathbb{E}_{\bm{x}}\left[p(\bm{x}|\bm{y},\bm{\gamma},\lambda)\right]=\bm{\mu}$ once the parameters $\bm{\gamma}$ and $\lambda$ have been estimated. These parameters are determined from the observations $\bm{y}$ by approximating the hyperparameter posterior in (3) by its mode (i.e., $p(\bm{\gamma},\lambda|\bm{y})\approx\delta(\bm{\gamma}^{*},\lambda^{*})$), which is determined by maximizing the evidence for the measurements $p(\bm{y}|\bm{\gamma},\lambda)$. Because the evidence is a convolution of Gaussians, it can be computed from the quantities defined in the probability model as $p(\bm{y}|\bm{\gamma},\lambda)=\int p(\bm{y}|\bm{x},\lambda)p(\bm{x}|\bm{\gamma})d\bm{x}=\mathcal{N}_{\bm{y}}(\bm{0},\bm{C})$, where

[TABLE]

The maximization of the evidence $p(\bm{y}|\bm{\gamma},\lambda)$, which is often referred to as type-II maximum likelihood, intuitively represents selecting the hyperparameters that are best supported by the observed data $\bm{y}$. The SBL procedure maximizes the marginal likelihood and hyperparameters over a log scale444Maximizing the hyperparameters over a log scale results in Jeffreys scale invariant hyperpriors; see [7, Appx. A], [17, Sec. 2.2.2—2.2.3], and [18] for more details.. Placing an uninformative prior on the noise variance $\lambda$ (i.e., $c=d=0$) results in the nonconvex negative marginal log-likelihood

[TABLE]

Note that when uninformative priors are also placed on the $\gamma_{i}$ (i.e., $a_{i}=b_{i}=0$ for all $i$), as is the case in the traditional SBL algorithm, the cost function is simply $\ell_{\mathrm{uninf}}(\bm{\gamma},\lambda)=\log\left|\bm{C}\right|+\bm{y}^{T}\bm{C}^{-1}\bm{y}$.

Several methods have been used in the SBL literature to iteratively minimize the cost function (6). The commonly used expectation-maximization procedure can be derived by treating $\bm{x}$ as the hidden variable, resulting in the iterative updates [7, 8]

[TABLE]

and

[TABLE]

During the EM procedure, many variances $\gamma_{i}$ become very small, indicating that the corresponding element $x_{i}$ is zero. To avoid numerical instability, when any $\gamma_{i}$ falls below some small threshold $\tau$ (e.g., $\tau=10^{-12}$), the corresponding $\gamma_{i}$ and dictionary column $\bm{\Phi}_{i}$ are pruned from the model and $x_{i}$ is fixed to zero. We denote the set of unpruned indices by $\mathcal{T}=\left\{i\colon\gamma_{i}\geq\tau\right\}$. An advantageous property of performing this pruning is that each EM iteration (whose cost is dominated by the $O(N^{3})$ matrix inverse required to recompute $\bm{\Sigma}$) becomes significantly faster as more elements are pruned from the model, so in practice larger values of $\tau$ (e.g., $\tau=10^{-4}$) are sometimes used.

2.2 Dynamic filtering of sparse signals

The Kalman filter (KF) [19] is the foundational algorithm for tracking time-varying signals. By assuming the dynamics model is linear, the states have Markovian structure, and the state evolution and observations are corrupted by Gaussian noise, the KF admits optimal closed-form expressions that recursively estimate the state at each time step. Extensions of the KF have also relaxed the linear and Gaussian assumptions by linearizing nonlinear dynamics models [20] or representing the state distribution using carefully chosen points [21, 22]. Further, although the KF framework is causal, Kalman smoothing algorithms (e.g., [20]) refine the estimate of previous time steps using the complete signal estimate.

However, these methods do not exploit the low-dimensional (sparse) structure present in many applications, which can significantly improve estimation accuracy. One body of work exploits sparse structure with the assumption that the signal support is constant or changing only slowly by incorporating a group $\ell_{1}$ penalty into the KF [23], estimating sparse outliers in the measurements and dynamics model jointly with the states [24], or using support estimation as a preprocessing step for the KF [25, 26]. Other methods exploit the sparsity of time-varying signals even when the support is varying rapidly by modifying the KF measurement model to include “pseudo-measurements” that have the effect of encouraging sparsity, inserting a modified orthogonal matching pursuit procedure into the KF procedure [27], extending the Kalman smoothing framework to use a broader class of penalty functions that promote sparsity [28, 29], or modifying the $\ell_{1}$ homotopy algorithm to use the KF framework [30] or sliding window-based processing [31].

Another family of algorithms start from an optimization perspective, augmenting the least-squares objective with $\ell_{1}$ terms that either directly encourage the signal and/or dynamics innovations to be sparse [32] or encourage the signal to be close to a dynamics-based estimate using $\ell_{2}$ [33, Sec. 3] or optimal transport-inspired [34] penalties.

The $\ell_{1}$ methods most similar to our proposed algorithm are methods that model the time-varying signal using a probabilistic framework by setting the priors of future time steps using dynamics information. In [35], a weighted $\ell_{1}$ penalty is used with weights set using the previous signal state and a probabilistic dynamics model. In the RWL1-DF algorithm of [33], dynamics information is incorporated in a robust manner by using a dynamics-based signal prediction to set the hyperpriors of the reweighted $\ell_{1}$ estimator [36, 37]. While these methods use probabilistic models based on the weighted or reweighted $\ell_{1}$ estimator, our SBL-DF algorithm propagates dynamics information through the hyperpriors of the SBL probability model, giving the resulting algorithm the same advantageous properties that SBL has been shown to have.

2.3 Dynamic filtering of sparse signals using SBL

In contrast to the potpourri of algorithms for dynamic filtering of sparse signals using standard $\ell_{1}$-based methods, relatively few algorithms exist for incorporating dynamics information into the SBL framework — despite its widespread use in practical applications.

SBL was initially extended to the offline recovery of dynamic sparse signals with constant support in the MSBL algorithm by assigning a single variance hyperparameter to each row of the matrix $\bm{X}=\left[\bm{x}^{(1)},\bm{x}^{(2)},\dots,\bm{x}^{(L)}\right]$ and jointly recovering the entire sequence of signals [38]. The TSBL algorithm improves this method by incorporating temporal correlation information into the SBL probability model [39]; an efficient online version is derived in [40]. These methods, however, assume that the signal support is constant over the entire time interval.

The constant support requirement is relaxed in [41] by reconstructing the signal over a small sliding window on which the signal support is assumed to be approximately stationary. A similar strategy is used in [42], which reconstructs the signal on each time window using an efficient variational SBL algorithm [43] rather than TSBL. The sliding window approach is improved in [44], which designs windows using a lapped orthogonal transform (similar to [31]) to avoid blocking artifacts and derives an efficient SBL algorithm to add and remove measurements based on the fast marginal likelihood updates of [16]. These methods remove the constant support assumption of MSBL and TSBL, but do not allow the use of a dynamics model and are ill-suited to applications where the signal is changing rapidly.

Echoing similar developments for the $\ell_{1}$ case, [45, 46] incorporate the SBL algorithm into the Kalman framework by replacing the Kalman filter’s Gaussian prior on the dynamics innovations with the Gaussian-inverse Gamma prior of the SBL framework, then performing SBL inference during each time step. Unlike previous methods, this algorithm allows online recovery of time-varying signals using SBL, but in contrast to our method, assumes the dynamics innovations are sparse rather than the signal itself. This sparse-innovations signal model can result in non-sparse states, and is therefore not the appropriate signal model for many applications.

Most similar to our proposed algorithm are methods that incorporate a dynamics-based signal prediction into the SBL probability model or objective function. In [47, 48], an estimate of the signal support is incorporated into the SBL probability model by adding a third-level hyperprior, then performing inference using a variational procedure. Although incorporating the support estimate into hyperpriors in this way makes the model more robust than methods that rigidly enforce (potentially erroneous) support estimates do, the method does not differentiate between small and large nonzero coefficients in the dynamics estimate, and requires an expensive variational procedure to perform inference on the more complex resulting model. In [49], the SBL framework is modified by replacing the hyperprior on the variance parameters with an improper Gaussian parameterized by the previous value of the signal. Inference in this modified model is performed on the signal innovations at each time step, reducing the strength of SBL’s sparsity encouraging mechanism for elements that are large in previous time steps. Unlike the sparse-states signal model we consider, this sparse-innovations model can lead to non-sparse states. Further, none of these methods allow the use of a dynamics model to improve inference.

3 Dynamic Filtering via Sparse Bayesian Learning (SBL-DF)

3.1 The SBL-DF algorithm

In this section, we introduce our algorithm for estimating time-varying sparse signals using SBL. The key insight of our approach is that the estimate of $\bm{x}^{(t)}$ can be improved in a robust manner by injecting information from the estimate of the previous time step and a dynamics model into the hyperparameters of the SBL probability model. A graphical depiction of this strategy is shown in Figure 1.

The dynamics-based prediction of the state at time $t$, $\widetilde{\bm{x}}^{(t)}$, is obtained by propagating $\widehat{\bm{x}}^{(t-1)}$ through the dynamics model $f_{t}\colon\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$. In applications where $f_{t}$ is known, our prediction is simply $\widetilde{\bm{x}}^{(t)}=f_{t}(\bm{x}^{(t-1)})$. Note that, unlike many algorithms based on the Kalman filter, we do not require the dynamics function to be linear. In other applications the dynamics model is unknown, but we expect that $\bm{x}^{(t)}\approx\bm{x}^{(t-1)}$ (i.e., $\bm{x}$ varies only slowly over time). In this case, we simply use either the identity dynamics function $f_{t}(\bm{x}^{(t)})=\bm{x}^{(t-1)}$ or select an $f_{t}$ that reflects a “distribution” on $\bm{x}^{(t)}$ based on $\bm{x}^{(t-1)}$ (for instance, a Gaussian blurring kernel centered on $\bm{x}^{(t-1)}$). In either case, because an important advantage of our method for propagating dynamics information through hyperpriors is the robustness of the resulting model, the use of an inexact dynamics model still provides good results. For example, we show in Section 4 that the identity dynamics function works well even when $\bm{x}$ is slowly time-varying. Note that although our model does not explicitly estimate or exploit temporal correlation structure as batch estimation procedures such as [39] do, if a priori correlation information is known it can be incorporated into the dynamics model $f_{t}$.

We next describe our specific method for mapping the dynamics-based prediction $\widetilde{\bm{x}}^{(t)}=f_{t}\left(\widehat{\bm{x}}^{(t-1)}\right)$ to the hyperparameters of the SBL probability model. For clarity of notation, we drop the superscript $t$, denoting $\widehat{\bm{x}}^{(t)}$ by $\widehat{\bm{x}}$ and $\widetilde{\bm{x}}^{(t)}$ by $\widetilde{\bm{x}}$.

First, we note that the SBL objective (6) can be decomposed as $\ell\left(\bm{\gamma},\lambda\right)=\ell_{\mathrm{uninf}}\left(\bm{\gamma},\lambda\right)+\ell_{\mathrm{dyn}}\left(\bm{\gamma},\lambda\right)$, where

[TABLE]

represents the dynamics-based portion of the objective that we will use to encourage fidelity to the dynamics estimate. The component $\ell_{\mathrm{uninf}}\left(\bm{\gamma},\lambda\right)$, which promote sparsity and fidelity to the measurements, is exactly the standard uninformative SBL objective. Intuitively, the overall magnitudes of $a_{i}$ and $b_{i}$ can be interpreted as controlling how strongly $\ell_{\mathrm{dyn}}\left(\bm{\gamma},\lambda\right)$ is weighed in the objective, and the relative magnitudes of $a_{i}$ and $b_{i}$ can be interpreted as controlling the specific value of $\gamma_{i}$ that $\ell_{\mathrm{dyn}}\left(\bm{\gamma},\lambda\right)$ encourages.

We select the hyperparameter values $\left\{a_{i},b_{i}\right\}$ so that the $\bm{\gamma}$ minimizing $\ell_{\mathrm{dyn}}\left(\bm{\gamma},\lambda\right)$ coincides with the $\bm{\gamma}$ minimizing the expected difference between the SBL-DF signal estimate $\widehat{\bm{x}}$ and the dynamics-based prediction $\widetilde{\bm{x}}$. This selection of $\left\{a_{i},b_{i}\right\}$ is accomplished in two steps. First, we determine the variance parameters as

[TABLE]

A closed-form expression for (10) cannot be found in general, so we make the simplifying assumption that $\bm{\Phi}^{T}\bm{\Phi}$ is diagonal (i.e., $\bm{\Phi}$ has orthogonal columns, but is not necessarily orthonormal). We have found that this approximation works well in practice even when the columns of $\bm{\Phi}$ are not orthogonal. Defining $g_{i}=\bm{\Phi}_{i}^{T}\bm{\Phi}_{i}$ and taking the derivative (with respect to $\bm{\gamma}^{-1}$ for simplicity), we have

[TABLE]

where the first equality follows from substituting the SBL posterior mean in (4) for $\widehat{\bm{x}}$ and then the measurement model (2) for $\bm{y}$. Setting to zero and solving for $\gamma_{i}$ yields

[TABLE]

for each component. This $\bm{\gamma}$ represents the variance parameters that would be optimal if the dynamics estimate $\widetilde{\bm{x}}$ was perfect, and is therefore the value of $\bm{\gamma}$ that we would like the dynamics-based terms of the informative SBL objective $\ell(\bm{\gamma},\lambda)$ to encourage.

We now note that for any fixed $a_{i},b_{i}$, the $\gamma_{i}$ minimizing $\ell_{\mathrm{dyn}}$ is

[TABLE]

Equating the maps (11) and (12) gives the rule $b_{i}/a_{i}=\widetilde{x}_{i}^{2}$. Incorporating a multiplicative “trade-off” parameter gives the dynamics mapping formula

[TABLE]

where the parameter $\xi$ represents how much weight the dynamics-based prediction is assigned in the evidence maximization procedure. Larger values of $\xi$ weigh the dynamics estimate more strongly when selecting $\bm{\gamma}$ and therefore reflect greater confidence in the accuracy of the dynamics model; see Section 4.1 for details.

Note that although the SBL-DF algorithm as presented runs causally using only a single previous estimate $\widehat{\bm{x}}^{(t-1)}$, it can easily be extended to use multiple previous estimates or into a smoothing estimator by setting $a_{i}^{(t)}$ and $b_{i}^{(t)}$ based on the estimate of $\bm{x}$ at other prior or future time steps.

Algorithm 1 summarizes the complete SBL-DF algorithm using expectation-maximization iterations to perform inference in the informative hyperprior SBL model.

3.2 Implications of informative hyperprior SBL

In contrast to the standard uninformative hyperprior SBL model that has enjoyed success in the literature, the SBL-DF algorithm places informative hyperpriors on the variances $\bm{\gamma}$. In this section, we investigate some key properties of the objective function that results from this model and their effects on the SBL-DF algorithm.

3.2.1 “Effective” SBL prior

One method for interpreting the sparsity promoting properties of SBL is to consider the “effective” SBL prior obtained by marginalizing out the first layer Gaussian prior $p(\bm{x}|\bm{\gamma})$. In the SBL probability model, marginalizing out the hyperparameters $\bm{\gamma}$ shows us that each $p(x_{i})$ is the independent non-standardized Student’s-$t$ distribution [7]

[TABLE]

parameterized by location $\mu_{i}=0$, scale $\sigma_{i}=\sqrt{b_{i}/a_{i}}=\left|\widetilde{x}_{i}\right|$, and degrees of freedom $\nu_{i}=2a_{i}=2\xi$. As demonstrated in [7], the high kurtosis of this effective prior intuitively explains why the SBL procedure produces sparse solutions in practice.

From this perspective, we see that SBL-DF’s method for propagating dynamics information into the hyperpriors adheres to our intuition that elements with small predicted magnitude should be given an effective prior $p(x_{i}|a_{i},b_{i})$ tightly peaked about zero, while elements with larger predicted magnitudes should have wider effective priors, which reduces the evidence needed to infer a nonzero value for $x_{i}^{(t+1)}$. Figure 2 shows the effective prior for different values of $\widetilde{x}_{i}$ and $\xi$. We see that, as desired, the probability mass of the effective prior $p(x_{i}|a_{i},b_{i})$ increasingly concentrates near zero as $\widetilde{x}_{i}$ decreases. Meanwhile, as $\xi$ increases the improper hyperprior becomes stronger, increasing the effect of the dynamics estimate on the inferred signal.

3.2.2 Characterization of the informative hyperprior objective

Figure 3 shows the form of $\ell_{\mathrm{dyn}}\left(\gamma\right)$ in one dimension, demonstrating the contrasting effects of the terms $-2a\log\gamma^{-1}$ and $2b\gamma^{-1}$ on the SBL inference procedure. We observe that while the $-2a\log\gamma^{-1}$ term weakly encourages $\gamma$ to be zero-valued (regardless of the value of $\widetilde{x}$), the $2b\gamma^{-1}$ term encourages $\gamma$ to be nonzero when the magnitude of $\widetilde{x}$ is large. The action of $\ell_{\mathrm{dyn}}\left(\gamma\right)$ in SBL-DF can therefore be interpreted as mitigating SBL’s natural sparsity-promoting properties in an amount commensurate with the magnitude of the dynamics-based signal prediction, while more strongly encouraging sparsity when $\widetilde{x}$ is zero-valued.

The minima of the uninformative hyperprior SBL objective $\ell_{\mathrm{uninf}}\left(\bm{\gamma}\right)$ are characterized in [8, Thm. 2], which shows that all local minima $\widehat{\bm{x}}$ are sparse in the sense that $\left\|\widehat{\bm{x}}\right\|_{0}<M$. Figure 3 demonstrates that this property may not always be preserved in the informative hyperprior setting when there are many elements where $x_{i}=0$ but $\widetilde{x}_{i}\neq 0$, since when $\widetilde{x}_{i}$ is nonzero, $2b_{i}\gamma_{i}^{-1}=2\xi\widetilde{x}_{i}\gamma_{i}^{-1}\rightarrow\infty$ as $\gamma_{i}\rightarrow 0$.

3.2.3 Implementation using reweighted $\ell_{1}$ iterations

In [50, 10], Wipf et al. examine the properties of the SBL objective by demonstrating that inference can be performed with a reweighted $\ell_{1}$ minimization procedure (SBL-RWL1) using a modified weight update step. This formulation of the SBL inference procedure elucidates the connection between SBL and $\ell_{1}$ minimization algorithms, and reveals dictionary dependence in the weight updates of SBL-RWL1 that may be responsible for its superior performance in the presence of adversely structured dictionaries [13, 51]. In this section, we extend the results of [10] from the uninformative hyperprior setting to the informative hyperprior setting used by SBL-DF, and show how conclusions drawn in literature from uninformative hyperprior SBL-RWL1 [50, 10, 13, 51, 52] carry over to informative hyperprior SBL-RWL1 and therefore to SBL-DF.

Proposition 1: The informative hyperprior SBL objective (6) can be minimized using the modified reweighted $\ell_{1}$ iterations

[TABLE]

where $\bm{\phi}_{i}$ is the $i^{\mathrm{th}}$ column of $\bm{\Phi}$, the update to $\bm{z}$ is called the majorization step and the update to $\bm{x}$ is called the minimization step. Note that as in the EM procedure a pruning rule is used, fixing $x_{i}$ to zero and removing column/entry $i$ from $\bm{\Phi}$, $\bm{z}$, and $\bm{\gamma}$ when $\gamma_{i}\leq\tau$.

The derivation of this procedure closely follows the method of [50, 10], which takes a majorization-minimization approach to minimizing (6) using a majorizing function cleverly designed to have a form similar to the weighted LASSO objective. Here, we briefly describe the main differences from uninformative hyperprior SBL and their implications; the complete derivation is outlined in Appendix A.

As in [10], the majorizing function is formed as the sum of two upper bounding functions. The concave part of (6), $\log\left|\bm{C}\right|-2\sum_{i}a_{i}\log\gamma_{i}^{-1}$, is bounded using its concave conjugate. The addition of the term $-2\sum_{i}a_{i}\log\gamma_{i}^{-1}$ does not change the linear form of this bound, so the effect of the $\{a_{i}\}$ only appears in the majorization step. The convex part of (6), $\bm{y}^{T}\bm{C}^{-1}\bm{y}+2\sum_{i}b_{i}\gamma_{i}^{-1}$, is bounded by a penalized least-squares term computed by expanding $\bm{y}^{T}\bm{C}^{-1}\bm{y}$ using the Woodbury identity and recognizing the result as the solution to a ridge regression problem; in our case, we add the additional term $2\sum_{i}b_{i}\gamma_{i}^{-1}$ directly into the bound, so the effect of the $\{b_{i}\}$ appears in the minimization step. The resulting majorization-minimization procedure takes the form of the iterations in Proposition 1, and the complete SBL-DF algorithm using this procedure is listed in Algorithm 2.

When $b_{i}=0$ for all $i$, by considering the $z_{i}^{1/2}$ as “weights” in an $\ell_{1}$ minimization problem, (16) is exactly equivalent to the maximization step in the reweighted $\ell_{1}$ algorithm [36]. Further, when $a_{i}=b_{i}=0$ for all $i$, the procedure that iterates between (14) and (16) is equivalent to the uninformative hyperprior SBL-RWL1 procedure of [10].

This procedure facilitates a direct comparison between the SBL procedure and traditional sparsity-encouraging algorithms (see, e.g., [50, 10, 13]). In particular, this allows us to verify that some beneficial properties of standard (uninformative hyperprior) SBL still apply when using informative hyperpriors as in SBL-DF. For instance, consider the weight update (14). It can be shown [51] that (14) can be written as

[TABLE]

where the optimization problem is understood to be performed on the pruned model. The form of (17) allows us to verify that two critical properties of standard (uninformative hyperprior) SBL also hold for informative hyperprior SBL and therefore SBL-DF. First, unlike standard RWL1, the weight update is non-separable, that is, the value of the $i^{\text{th}}$ weight is dependent on all elements in $\bm{x}$ through the latent variables $\bm{\gamma}$. This contrasts with standard reweighted $\ell_{1}$, in which the weight on $x_{i}$ is dependent only on the the single element $x_{i}$ from the previous iteration. Second, the $\ell_{2}$ norm term in (17) reveals that the weights in SBL-RWL1 depend on the coherence structure of $\bm{\Phi}$, despite the inclusion of the nonzero hyperparameter $a_{i}$. This dependence allows SBL to moderate the penalty applied to coefficients that correspond to highly correlated columns of $\bm{\Phi}$ [51], potentially explaining its superior performance with structured dictionaries.

3.3 Improving the speed of SBL-DF using fast marginal likelihood

The EM iterations (7–8) used in the SBL algorithm to minimize $\ell(\bm{\gamma},\lambda)$ require an $O(\widetilde{N}^{3})$ matrix inversion at each iteration to recompute $\bm{\Sigma}$, where $\widetilde{N}=\left|\mathcal{T}=\left\{\gamma_{i}>\tau\right\}\right|$ denotes the number of unpruned elements remaining in the model. In [16], Tipping and Faul propose a more efficient method, the fast marginal likelihood (FML) algorithm, for minimizing the uninformative SBL objective $\ell_{\mathrm{uninf}}(\bm{\gamma},\lambda)$. The key idea of the method in [16] is that the effect of updating an individual coefficient on the state $\bm{\Sigma}^{-1}$, which is required to calculate the state estimate $\bm{\mu}$, can be computed inexpensively by a rank-one update. Therefore, by updating only a single coefficient at each iteration, we can avoid the full recomputation of $\bm{\Sigma}^{-1}$ that each EM iteration requires. Moreover, by strategically choosing which coefficient to update at each iteration, convergence can occur with relatively few of these rank-one updates.

However, the FML updates described in [16] require the hyperparameters of the SBL probability model to be uninformative, so the FML procedure cannot directly be applied to the SBL-DF algorithm. In this section, we derive a generalized FML procedure for the informative hyperprior setting (i.e., to minimize $\ell\left(\bm{\gamma},\lambda\right)$ when $\ell_{\mathrm{dyn}}\left(\bm{\gamma},\lambda\right)\neq 0$) and adapt it to the SBL-DF algorithm. We show that the FML procedure can be extended to the informative hyperprior SBL model at the cost of a more computationally expensive procedure to choose the best coefficient to update at each iteration, an additional cost that grows as $O(N)$.

The critical fact that enables the FML method is that the covariance matrix $\bm{C}$ of the marginal likelihood — and therefore the marginal likelihood $\ell\left(\bm{\gamma},\lambda\right)$ itself — is separable in $\bm{\gamma}$:

[TABLE]

This decomposition allows basis vectors to be efficiently added and removed from the model using cheap rank-one updates to $\bm{C}$. To see this, note that the covariance matrix excluding the $j^{\text{th}}$ element (denoted $\bm{C}_{-j}$) can be calculated as $\bm{C}_{-j}=\bm{C}-\gamma_{j}\bm{\phi}_{j}\bm{\phi}_{j}^{T}$. Substituting these expressions into the objective (6) and then separating the terms involving $\gamma_{j}$ from those involving $\bm{\gamma}_{-j}$ and $\lambda$ shows that the negative log-likelihood (6) can be decomposed as $\ell(\bm{\gamma},\lambda)=\ell(\bm{\gamma}_{-j},\lambda)+\ell(\gamma_{j})$ where

[TABLE]

and

[TABLE]

where, following the notation of [16], $s_{j}=\bm{\phi}_{j}^{T}\bm{C}_{-j}^{-1}\bm{\phi}_{j}$ and $q_{j}=\bm{\phi}_{j}^{T}\bm{C}_{-j}^{-1}\bm{y}$. Critically, this decomposition reveals that finding the value of $\gamma_{j}$ that minimizes $\ell(\bm{\gamma},\lambda)$ is equivalent to finding the value of $\gamma_{j}$ that minimizes $\ell(\gamma_{j})$. We compute this value, which we denote $\widetilde{\gamma}_{j}$, by finding the roots of the derivative

[TABLE]

where $c_{3}=\left(\tfrac{1}{2}+a_{j}\right)s_{j}^{2}$, $c_{2}=\left[\left(\tfrac{1}{2}+2a_{j}\right)s_{j}-\tfrac{1}{2}q_{j}^{2}-b_{j}s_{j}^{2}\right]$, $c_{1}=\left(a_{j}-2b_{j}s_{j}\right)$, and $c_{0}=-b_{j}$. Roots occur when either $\gamma_{j}\rightarrow+\infty$ (from both the $\gamma_{j}$ and $(\gamma_{j}+s_{j}^{-1})^{2}$ terms in the denominator) or when $\gamma_{j}$ is a root of the cubic expression in parentheses; its roots can be computed analytically or by calculating the eigenvalues of its companion matrix. See Appendix B for details.

We immediately remove nonreal fixed points $\widetilde{\gamma}_{j}$, and of the remaining points, keep only the $\widetilde{\gamma}_{j}$ corresponding to local minima of $\ell(\gamma_{j})$ by removing roots where the second derivative,

[TABLE]

is negative. If there are multiple real-valued local minima, we choose the one that results in the smallest value of the negative log likelihood $\ell(\gamma_{j})$. The final value of the root is denoted $\widetilde{\gamma}_{j}$.

With this root-finding procedure for determining $\widetilde{\gamma}_{j}$ in hand, we can now describe the complete FML algorithm. At each iteration, we compute $\widetilde{\gamma}_{i}$ for $i=1,\dots,N$. The calculated value of $\widetilde{\gamma}_{i}$ defines both the action we would take to update that element (either re-estimate, add, or delete) and the change in marginal likelihood $\Delta_{i}$ that performing this action would result in:

• •

If $i\in\mathcal{T}$ and $\widetilde{\gamma}_{i}>\tau$, re-est.: $\gamma_{i}\leftarrow\widetilde{\gamma}_{i}$, $\Delta_{i}=\ell\left(\widetilde{\gamma_{i}}\right)-\ell\left(\gamma_{i}\right)$.

• •

If $i\notin\mathcal{T}$ and $\widetilde{\gamma}_{i}>\tau$, add: $\gamma_{i}\leftarrow\widetilde{\gamma}_{i}$, $\Delta_{i}=\ell\left(\widetilde{\gamma_{i}}\right)$.

• •

If $i\in\mathcal{T}$ and $\widetilde{\gamma}_{i}\leq\tau$, delete: $\gamma_{i},~{}x_{i}\leftarrow 0$, $\Delta_{i}=-\ell\left(\gamma_{i}\right)$.

Once $\Delta_{i}$ has been calculated for $i=1,\dots N$, we perform the action that results in the greatest change in marginal likelihood, updating the element at index $j=\arg\max_{i}\Delta_{i}$. Performing this action can be performed particularly efficiently because by judiciously representing each element by [16]

[TABLE]

the state updates to $\bm{\Sigma}$, $\bm{\mu}$, and $\{S_{j},Q_{j},s_{j},q_{j}\}$ needed to perform the re-estimate, add, and delete actions can be implemented using inexpensive low-rank updates (see Appendix B). The complete SBL-DF algorithm using the FML method for inference is listed in Algorithm 3.

In summary, the FML algorithm initializes an empty model and iteratively re-estimates, adds, or deletes the element that results in the largest improvement of the marginal likelihood. By judiciously representing the state of the model at each iteration using (21–22), each update can be implemented using a few inexpensive rank-one updates. In contrast to the FML algorithm of [16], which is restricted to the uninformative hyperprior SBL probability model, the FML algorithm derived here allows general hyperparamers $a_{i}$ and $b_{i}$ to be used at the cost of a more complicated root-finding procedure to compute $\widetilde{\bm{\gamma}}$. Because this root-finding procedure must be run for $i=1,\dots,N$ at each iteration, the total additional cost introduced by the use of informative hyperpriors grows as $O(N)$.

4 Numerical Simulations

In this section, we evaluate the efficacy of the SBL-DF algorithm with synthetic data. These experiments with known ground truth allow us to characterize the performance of SBL-DF in many different regimes.

First, we demonstrate the performance of SBL-DF’s method for incorporating dynamics information using a single time step, i.e., we recover a single length-$N$ signal using a noisy signal prediction that mimics the prediction we would obtain using a dynamics model. This allows us to explore SBL-DF’s performance as the quality of the dynamics prior information (i.e., accuracy of $\widetilde{\bm{x}}=\left|f_{t}\left(\widehat{\bm{x}}^{(t)}\right)\right|$) changes. We then show that the single time step performance carries over to tracking problems with multiple time steps.

4.1 Performance as dynamics prediction quality degrades

We first show how incorporating a dynamics-based prediction can improve reconstruction accuracy, and how reconstruction accuracy is affected by the quality of the dynamics-based prediction. In this experiment, we assume that our dynamics-based estimate $\widetilde{\bm{x}}=\left|f_{t}\left(\widehat{\bm{x}}^{(t)}\right)\right|$ is the ground truth signal corrupted by one of two types of error: (a) support mismatch, in which the values of $\bar{s}$ nonzero elements are swapped with randomly chosen zero-valued elements; and (b) Gaussian noise, in which i.i.d. mean-zero Gaussian noise with variance $\sigma_{\mathrm{dyn}}^{2}$ is added to the entire signal estimate. We generate data using the measurement model (2) with $N=512$, $s=16$ nonzero elements, $\bm{e}\sim\mathcal{N}\left(\bm{0},\sigma^{2}_{\mathrm{obs}}\bm{I}\right)$ with $\sigma^{2}_{\mathrm{obs}}=10^{-3}$, and construct the dictionary $\bm{\Phi}$ with i.i.d. $\mathcal{N}(0,1/\sqrt{M})$ entries. For recovery, we use FML (Algorithm 3) with threshold $\tau=10^{-1}$ and fixed regularization parameter $\lambda=10^{-3}$. We select the value of the dynamics parameter $\xi$ that minimizes reconstruction error from a fixed grid of $\log_{10}\left(\xi\right)$ from $-2$ to $2$ in increments of $0.1$. As $\xi$ represents how strongly the dynamics-based prediction is considered in the inference procedure, its optimal value is strongly dependent on the quality of the dynamics estimate and therefore a separate value of $\xi$ was selected for each dynamics error level. We vary $M$ for each dynamics noise level, performing $240$ independent trials for each point, and plot the rate of successful recovery (defined as the portion of trials satisfying relative MSE (rMSE) $\left\|\bm{x}-\widehat{\bm{x}}\right\|_{2}^{2}/\left\|\bm{x}\right\|_{2}^{2}<10^{-2}$).

Figure 4 displays the rate of successful recovery when the prediction $\widetilde{\bm{x}}$ is the true signal corrupted with an increasing number of support errors or magnitude of Gaussian noise, respectively. The yellow dotted line represents the performance of the standard uninformative hyperprior SBL algorithm. Both figures demonstrate that SBL-DF’s strategy for including a dynamics-based prediction into the SBL probability model results in vastly improved performance when the dynamics-based prediction is accurate in the sense that the true signal can be recovered accurately with far fewer measurements than traditional SBL requires. Further, even when the prediction is highly flawed, SBL-DF with a reasonable choice of $\xi$ still allows $\widehat{\bm{x}}^{(t+1)}$ to be accurately recovered with fewer measurements than the standard SBL algorithm.

Figure 5 examines the effect of changing the $\xi$ parameter (which controls the weight SBL-DF gives the dynamics-based prediction) using the same experimental setup as Figure 4. In Figure 5 (left), we observe that when the dynamics-based prediction $\widetilde{\bm{x}}$ is accurate, the reconstruction error (rMSE) decreases monotonically with increasing $\xi$. When the dynamics-based prediction is very inaccurate, there is a value of $\xi$ after which increasing $\xi$ (i.e., enforcing the dynamics-based prediction more strongly in inference) decreases reconstruction accuracy. In Figure 5 (right), we observe that when the dynamics-based prediction $\widetilde{\bm{x}}$ is accurate, the reconstruction error (rMSE) decreases monotonically with increasing $\xi$. However, when the dynamics-based prediction is noisy, there is a value of $\xi$ after which increasing $\xi$ decreases reconstruction accuracy.

4.2 Performance with structured dictionaries

Next, we explore the performance of SBL-DF when the dictionary contains challenging structure such as highly correlated columns or column norms of drastically different magnitudes. Although this type of dictionary structure is common in many practical applications (e.g., neural imaging [13]), $\ell_{1}$ minimization-based algorithms provably fail in this setting. In contrast, traditional SBL has been shown to be provably and empirically robust to this type of dictionary structure, and in fact SBL is commonly used in applications with unfavorably structured dictionaries [13, 53, 51]. In this experiment, we explore how this robustness extends to the dynamic case where prior knowledge is incorporated. Specifically, we compare the SBL-DF algorithm with RWL1-DF [33], a state of the art $\ell_{1}$ minimization based tracking algorithm based on reweighted $\ell_{1}$ [36] when this type of undesirable dictionary structure is present.

We consider four dictionary models constructed from $\widetilde{\bm{\Phi}}\in\mathbb{R}^{M\times N}$ containing i.i.d. $\mathcal{N}(0,1/\sqrt{M})$ entries:

• •

i.i.d.: $\bm{\Phi}=\widetilde{\bm{\Phi}}$.

• •

i.i.d. with scaled columns: $\bm{\Phi}=\widetilde{\bm{\Phi}}\bm{D}$, where $\bm{D}\in\mathbb{R}^{N\times N}$ is the diagonal matrix with i.i.d. $U[0,1]$ entries on the diagonal.

• •

Locally coherent: $\bm{\Phi}=\widetilde{\bm{\Phi}}\bm{B}$, where $\bm{B}\in\mathbb{R}^{N\times N}$ is a block diagonal matrix consisting of $4\times 4$ blocks, each with $1$ on its diagonal and $0.8$ elsewhere. This model, which is similar to the setup of the structured dictionary in [51], results in a dictionary with $25$ groups of $4$ highly correlated columns.

• •

Locally coherent with scaled columns: $\bm{\Phi}=\widetilde{\bm{\Phi}}\bm{B}\bm{D}$, where $\bm{B}$ and $\bm{D}$ are defined as above. This model matches the setup of the structured dictionary in [51].

To standardize the effective SNR of each model, we scale each of the resulting dictionaries by the scalar $\lVert\widetilde{\bm{\Phi}}\bm{x}\rVert_{2}/\lVert\bm{\Phi}\bm{x}\rVert_{2}$.

We also consider two models for the signal $\bm{x}$. First, we consider the model where the $s$ nonzero elements of $\bm{x}$ are i.i.d. $\mathcal{N}(0,1)$. This setting, in which SBL has been shown to outperform other sparse estimation methods [9], is frequently considered in the SBL literature. Second, we consider the model where the $s$ nonzero elements of $\bm{x}$ all have value $1$. This setting is more frequently considered in traditional compressive sensing literature.

In this experiment, we set $N=100$, $M=42$, $s=25$, and sweep the observation noise $\sigma^{2}_{\mathrm{obs}}$ from $10^{-8}$ to $10^{-3}$. Note that although this is a relatively small amount of noise in the prototypical i.i.d. dictionary setting, the structured dictionary models we consider here make reconstruction extremely challenging at this noise level. To simulate the dynamics-based prediction $\widetilde{\bm{x}}$, we add support errors to $\bm{x}$ by switching the value of each nonzero element of $\bm{x}$ with a zero-valued element with probability $p=0.1$ and then add white Gaussian noise with variance $\sigma_{\mathrm{dyn}}^{2}=10^{-4}$ to every element of $\widetilde{\bm{x}}$. We perform SBL-DF inference using the EM updates of Algorithm 1 with pruning threshold $\tau=10^{-4}$, select the dynamics parameter $\xi$ by a grid search run independently for each point, and use the automatic learning rule (8) to select the regularization parameter $\lambda$. For RWL1-DF, the three parameters $\lambda$, $\xi$, $\eta$ were selected by a grid search run independently for each point.

Figure 6 shows the rMSE $\left\|\bm{x}-\widehat{\bm{x}}\right\|_{2}^{2}/\left\|\bm{x}\right\|_{2}^{2}$ as the observation noise variance $\sigma_{\mathrm{obs}}^{2}$ increases. In this figure, the main line shows the median rMSE, and error bars show the $25$th and $75$th percentile rMSE to give a sense of the reconstruction error variance. We observe that with an i.i.d. dictionary, both RWL1-DF and SBL-DF are able to reconstruct the signal accurately, with error beginning to increase to unacceptable levels only as $\sigma_{\mathrm{obs}}^{2}$ grows large. However, when the dictionary has the challenging types of structure considered here, SBL-DF tends to perform better than RWL1-DF, recovering the signal with lower error (for i.i.d. with scaled columns and locally coherent with scaled columns models) or lower variance (for the locally coherent model) at high to moderate SNRs. We have observed that SBL-DF similarly outperforms RWL1-DF when using other models of coherent dictionaries, such as dictionaries containing i.i.d. Gaussian entries with nonzero mean. The robustness displayed by SBL-DF aligns with theoretical results in the static SBL literature, which has shown that in the absence of measurement noise, the uninformative hyperprior SBL objective is invariant to column scaling and coherence structure. In contrast, the objective of $\ell_{1}$ based methods are invariant to either column scaling or coherence only in isolation [13].

In addition to superior reconstruction accuracy, SBL-DF also has the significant advantage of only requiring the tuning of a single parameter555Although the pruning threshold $\tau$ could nominally be considered a tuning parameter, it can be typically be set to any reasonable value without a large change in results. Here, we fix $\tau=10^{-4}$ in all experiments using the EM algorithm.. In contrast, RWL1-DF requires tuning three parameters. Since in this experiment we perform a grid search over only one parameter for SBL-DF, and a simultaneous grid search over three parameters for RWL1-DF, the results in this experiment are charitable for RWL1-DF and some degradation in RWL1-DF’s performance would be expected in applications where it is not possible to perform this large grid search.

We further explore SBL-DF’s performance in the presence of challenging dictionary structure in Figure 7. In this experiment, we again set $N=100$, $M=42$, $s=25$, but now fix $\sigma_{\mathrm{obs}}^{2}$ to three different values. As in the previous experiment, the dynamics-based prediction is the ground truth signal corrupted with support errors occurring with probability $p=0.1$ and white Gaussian noise with $\sigma_{\mathrm{dyn}}^{2}=10^{-4}$. The dictionary is constructed using the “locally coherent with scaled columns” model described above, but we now sweep the amount of structure in the dictionary. Specifically, each block of $\bm{B}$ now takes value $1-\tfrac{1}{c}$ off its diagonal and the entries of $\bm{D}$ are i.i.d. $U[\tfrac{1}{c},1]$ so that a larger value of $c$ results in a larger range of column scales. Recall that $\bm{\Phi}$ is scaled to have the same norm as $\widetilde{\bm{\Phi}}$, so the overall scaling of $\bm{D}$ does not affect the SNR. This parameterization results in the dictionary being equivalent to the i.i.d. model described above when $c=1$; as $c$ increases, the dictionary becomes increasingly structured. For each noise level, we sweep the structure parameter $c$, performing $40$ independent trials at each point and plotting the rate of successful recovery (defined as rMSE $\left\|\bm{x}-\widehat{\bm{x}}\right\|_{2}^{2}/\left\|\bm{x}\right\|_{2}^{2}<10^{-2}$). As before we sweep the single parameter $\xi$ for SBL-DF using a grid search for each noise and structure level, and sweep the three parameters, $\lambda$, $\xi$, and $\eta$, for RWL1-DF using a simultaneous grid search for each noise and structure level.

From Figure 7 we observe that, in general, increasing dictionary structure makes it more difficult for both algorithms to accurately recover the signal. However, in contrast to RWL1-DF, in high SNR regimes SBL-DF is able to accurately recover the signal even with very highly structured dictionaries. As the SNR becomes lower and reconstruction becomes increasingly difficult with any amount of structure, both algorithms are able to accurately recover the signal when the dictionary has i.i.d. entries ($c=1$), but not when the amount of structure increases.

4.3 Runtime comparison

Next, we compare the runtime of SBL and SBL-DF using both the EM and FML inference procedures. In this experiment we fix the underdeterminedness ratio to $M/N=1/4$ and the number of nonzero entries to $s=16$, and sweep the signal length from $N=2^{9}$ to $N=2^{14}$. We use the dictionary $\bm{\Phi}\in\mathbb{R}^{M\times N}$ containing i.i.d. $\mathcal{N}(0,1/\sqrt{M})$ entries, set the observation noise variance to $\sigma^{2}_{\mathrm{obs}}=10^{-3}$, and generate the dynamics-based prediction as the true signal corrupted with support errors occurring with probability $p=0.1$ and Gaussian noise with variance $\sigma_{\mathrm{dyn}}^{2}=10^{-4}$. As in the previous experiments, we set the pruning threshold to $\tau=10^{-4}$ for EM and $\tau=10^{-1}$ for FML, fix $\lambda=1.2\times 10^{-3}$ for both EM and FML, and set $\xi=10^{0}$ for SBL-DF. Each algorithm is run until convergence666The FML algorithm maintains and updates values that can be used to efficiently calculate the change in the log likelihood at each iteration, which is used to determine convergence. In contrast, our EM implementation uses $\left\|\bm{\mu}^{\mathrm{(new)}}-\bm{\mu}\right\|_{2}$ as the convergence criteria as it does not maintain the value of the log likelihood. Therefore, to facilitate a fair comparison, in these experiments we use $\left\|\bm{\mu}^{\mathrm{(new)}}-\bm{\mu}\right\|_{2}<\mathrm{tol}$ as the convergence criteria for both algorithms. with a tolerance of $10^{-4}$. This tolerance was selected as it is approximately the value where the algorithm achieves its final error. For this experiment, SBL (EM) and SBL (FML) had similar rMSE, and SBL-DF (EM) and SBL-DF (FML) had similar rMSE. Each trial in this section was run on a single core of an Intel Xeon E5-2680 v3 processor with a base clock rate of 2.5GHz and a shared pool of 100GB memory.

Figure 8 shows the total execution time for each algorithm as the problem size increases, with each point displaying the mean of $24$ independent trials. We first compare the execution times of SBL (EM) and SBL-DF (EM). We see that despite the fact that the initial EM iterations for SBL and SBL-DF have identical $O\left(\widetilde{N}^{3}+\widetilde{N}^{2}M+\widetilde{N}^{2}+\widetilde{N}M+\widetilde{N}\right)$ per-iteration cost, SBL-DF’s incorporation of dynamics information with informative hyperpriors reduces the number of iterations required for convergence by approximately one order of magnitude. Figure 9, which shows the execution time of each iteration, suggests that this improvement is due to SBL-DF inferring and pruning zero-valued elements (and therefore reducing $\widetilde{N}$) much more quickly than SBL.

The relative performance of SBL (FML) and SBL-DF (FML) is slightly more nuanced. The incorporation of informative hyperpriors into the FML method for SBL-DF requires a root-finding procedure that adds a performance penalty to each iteration of SBL-DF (FML) when compared to SBL (FML). However, the cost of this root-finding procedure grows only linearly with $N$, and as the problem size grows this cost is dominated by the $O\left(\widetilde{N}^{2}+MN+N\right)$ matrix operations required for computing the low-rank state updates in each iteration. Additionally, Figure 9 shows that the incorporation of prior information allows SBL-DF to converge with many fewer iterations of the “re-estimate” action. The overall effect of these two phenomena is that SBL-DF (FML) is more expensive than SBL (FML) on small problems, but becomes more efficient than SBL (FML) as $N$ and $M$ increase.

Comparing the execution times of the EM and FML inference procedures, we observe that when using standard (uninformative hyperprior) SBL, the FML method is both significantly more efficient than EM and scales slightly better with the problem size. When using informative hyperpriors (as in SBL-DF), the FML method requires an additional computational burden from the added complexity of root calculation, but has execution time that scales significantly better with problem size than EM. We have observed that, unlike the EM algorithm, the number of iterations required by the FML method is strongly dependent on the sparsity level, so the EM method may still be preferable in large problems where $s/N$ is relatively large. Further, the convergence speed of SBL-DF is dependent on the quality of the dynamics estimate, with SBL-DF converging more quickly when the dynamics estimate is close to the true signal.

4.4 Dynamic filtering

Here, we demonstrate the efficacy of SBL-DF on an $L=30$ time step synthetic tracking problem. The following procedure (similar to that used in [33]) is used to generate synthetic tracking data. First, $\bm{x}^{(1)}\in\mathbb{R}^{N}$ is generated with $s$ nonzero targets drawn from $\mathcal{N}\left(0,1\right)$, but with elements where $\left|x_{i}\right|<0.1$ set to $0.1\times\mathrm{sign}(x_{i})$. To generate the remaining time steps, $\bm{x}^{(2)},\dots,\bm{x}^{(L)}$, each target is assigned a direction (i.e., $-1$ or $+1$), and the linear dynamics model $\bm{F}$ is generated to correspond to each target moving one index in its assigned direction at each time step. However, the actual signal experiences sparse innovations, with each element moving in the opposite direction with probability $p=0.1$ at each time step. We set $N=100$, $M=42$, $s=25$, $\sigma^{2}_{\mathrm{obs}}=10^{-6}$, and generate the dictionary using the “locally coherent with scaled columns” model described above. We compare to static SBL, static RWL1 [36] and RWL1-DF [33]. In addition, we compare to the windowed TMSBL (WTMSBL) algorithm [41], which runs the TMSBL multiple measurement vector SBL algorithm on a series of overlapping windows and is representative of a body of literature that assumes the true signal is approximately stationary over small intervals.

Figure 10 shows the results of the multiple time step dynamic filtering simulation. We observe that both RWL1-DF and SBL-DF outperform their static counterparts that do not take advantage of a dynamics based signal estimate. However, SBL literature has shown that this type of dictionary structure is particularly detrimental to $\ell_{1}$-based algorithms, and here we observe that static SBL outperforms both static RWL1 and RWL1-DF. SBL-DF further improves the performance of SBL, allowing accurate reconstruction even in the presence of challenging dictionary structure and observation noise. We observe that, due to the challenging structure in the dictionary [13], static SBL outperforms RWL1. The windowed TMSBL method performs worse than the standard single measurement vector methods in this setting because the signal violates the slowly time-varying assumption. Both RWL1-DF and SBL-DF are able to improve reconstruction performance over their static counterparts in this setting by using an (imperfect) dynamics model.

5 Discussion

In this paper, we proposed an algorithm for tracking time-varying sparse signals using the sparse Bayesian learning (SBL) framework. We demonstrated that our method for incorporating prior state estimates and an imperfect dynamics model into the SBL probability model allows the algorithm to recover signals more accurately and efficiently than previous methods. Further, we showed that the algorithm outperforms state of the art $\ell_{1}$-based methods when the dictionary contains challenging structure such as coherence and diverse column magnitudes. In addition, we adapted the fast marginal likelihood [16] and reweighted $\ell_{1}$ [10] procedures to the informative hyperprior SBL setting, allowing their use with SBL-DF and allowing us to draw comparisons to static SBL.

In this work, we applied our proposed method to the causal reconstruction of time-varying signals, where the signal estimate is obtained from previously recovered time steps and a dynamics model. However, our general framework for mapping a signal estimate into the SBL probability model developed in Section 3 is applicable to any problem where an a priori estimate of a signal can be used to improve inference.

There are several avenues for future work. One potential improvement could be made by incorporating related work estimating correlation structure or dynamics models from data, allowing the algorithm to benefit from trends automatically discovered in previous time steps. SBL-DF is particularly well-poised to take advantage of research advances in this area due to its use of an arbitrary dynamics model. Additionally, SBL-DF’s probabilistic formulation allows other priors or signal models to be easily incorporated into the algorithm; extensions could add application-specific domain knowledge or signal characteristics such as cluster structure.

A limitation of SBL-DF is inherited from the lack of strong theoretical performance and convergence guarantees in SBL literature, which are not currently as developed as those for $\ell_{1}$-based methods. Theoretical advances in this direction could lead to improved understanding of where and why SBL-DF outperforms similar $\ell_{1}$-based methods, and more rigorously quantify SBL-DF’s performance relative to static SBL and other dynamic filtering algorithms.

Finally, while we have provided a fast algorithm based on the widely used fast marginal likelihood updates, leveraging recent work performing on SBL inference using approximate message passing [54] or variational methods [55] has the potential to admit implementations of SBL-DF that are highly efficient in high-dimensional settings or more amenable to modifications of the SBL probability model.

Appendix A Derivation of Reweighted $\ell_{1}$ Implementation

As in [50, 10], we take a majorization-minimization approach to minimizing (6). We begin by constructing a function that majorizes the SBL cost function (6) as a sum of two upper bounding functions: one that bounds the concave expression $\log\left|\bm{C}\right|-2\sum_{i}a_{i}\log\gamma_{i}^{-1}$, and one that bounds the convex expression $\bm{y}^{T}\bm{C}^{-1}\bm{y}+2\sum_{i}b_{i}\gamma_{i}^{-1}$.

First, we bound $\log\left|\bm{C}\right|-2\sum_{i}a_{i}\log\gamma_{i}^{-1}$ as

[TABLE]

where $g^{*}(\bm{z})=\inf_{\bm{\gamma}}\left[\bm{z}^{T}\bm{\gamma}-g(\bm{\gamma})\right]$ is the concave conjugate of $g(\bm{\gamma})=\log\left|\bm{C}\right|-2\sum_{i}a_{i}\log\gamma_{i}^{-1}$ (see [10] for details). Intuitively, $\bm{z}^{T}\bm{\gamma}-g^{*}(\bm{z})$ is a family of linear upper bounds on $g(\bm{\gamma})$; the tightest bound at a fixed value of $\bm{\gamma}$ is achieved by minimizing over $\bm{z}$. Note that we do not need to explicitly calculate $g^{*}(\bm{z})$ because we will only need to minimize the right hand side of (23) with respect to $\bm{\gamma}$.

Second, we bound the convex expression $\bm{y}^{T}\bm{C}^{-1}\bm{y}+2\sum_{i}b_{i}\gamma_{i}^{-1}$. By applying the matrix inversion lemma and recognizing the form of the result as the solution to the of a Tikhonov-regularized least-squares problem, $\bm{y}^{T}\bm{C}^{-1}\bm{y}$ can be written as [10] $\bm{y}^{T}\bm{C}^{-1}\bm{y}=\underset{\bm{x}}{\arg\min}~{}\frac{1}{\lambda}\left\|\bm{y}-\bm{\Phi}\bm{x}\right\|_{2}^{2}+\left\|\bm{\Gamma}^{-1/2}\bm{x}\right\|_{2}^{2}$, giving the bound

[TABLE]

Combining (23) and (24) gives the majorizing function777We have suppressed the dependence on $\lambda$ in $\ell_{\bm{z}}\left(\bm{\gamma},\bm{x}\right)$ and $\ell\left(\bm{\gamma}\right)$ for clarity.

[TABLE]

We now construct the majorization-minimization iterations. The majorization step updates the majorizing function by calculating the $\bm{z}$ that results in the tightest bound (25) at $\bm{\gamma}$. From duality theory, this value of $\bm{z}$ occurs at the gradient of $g(\bm{z})$ [56], so the majorization step is

[TABLE]

The minimization step consists of solving

[TABLE]

Fixing $\bm{x}$ and minimizing with respect to $\bm{\gamma}$ by setting the derivative to zero gives

[TABLE]

where we have taken the positive component of the radical because as variances, each $\gamma_{i}$ must be positive.

Substituting $\bm{\gamma}^{(k+1)}$ from (27) and minimizing with respect to $\bm{x}$ gives

[TABLE]

The complete majorization-minimization procedure alternates between updating the majorizer using (26) and (27), and minimizing it using (28). (Note that $\bm{C}$ depends on $\bm{\gamma}$ and must be updated at each iteration.) For numerical stability, a pruning rule similar to the one used in static SBL is used in the minimization step: when a $\gamma_{i}$ becomes lower than some threshold $\tau$, the corresponding $x_{i}$ is fixed to zero and dictionary column $\bm{\phi}_{i}$ is removed.

Appendix B Details of Fast Marginal Likelihood

B.1 Calculation of stationary points of $\ell(\gamma_{j})$

In this section, we give expressions for computing the roots of the cubic expression in (18) which are used to determine the optimal index and action for the next fast marginal likelihood iteration.

Defining constants $\Delta_{0}=c_{b}^{2}-3c_{a}c_{c}$, $\Delta_{1}=2c_{b}^{3}-9c_{a}c_{b}c_{c}+27c_{a}^{2}c_{d}$, and determinant $\Delta=18c_{a}c_{b}c_{c}c_{d}-4c_{b}^{3}c_{d}+c_{b}^{2}c_{c}^{2}-4c_{a}c_{c}^{3}-27c_{a}^{2}c_{d}^{2}$, the roots can be shown to be

[TABLE]

where $\zeta=-\frac{\sqrt{3}-1}{2}\sqrt{-1}$ and $C=\sqrt[3]{\frac{1}{2}\left(\Delta_{1}\pm\sqrt{\Delta_{1}^{2}-4\Delta_{0}^{3}}\right)}$. These three cases correspond to when there are non-repeated roots, a double root and a single root, and a triple root, respectively.

We can also find the roots as the eigenvalues of the companion matrix

[TABLE]

We have empirically found that this method is more numerically stable than the analytic expressions. See [57, Sec. 7.4.6] and references therein for more details on the companion matrix method and its stability.

B.2 Efficient updates for fast marginal likelihood

The following expressions efficiently implement the re-estimate, add, and delete actions to element $j$ in the fast marginal likelihood algorithm [16].

B.2.1 Re-estimating an element

Defining $\kappa_{j}=\left[\left(\widetilde{\gamma_{j}}^{-1}-\gamma_{j}^{-1}\right)^{-1}+\left(\bm{\Sigma}\right)_{jj}\right]^{-1}$ and repeatedly applying the matrix inversion lemma yields

[TABLE]

B.2.2 Adding an element

Defining $\omega_{j}=\left(\widetilde{\gamma_{j}}^{-1}+S_{j}\right)^{-1}$, using the block matrix inversion identity [58, Sec. 9.1.3] and applying the matrix inversion lemma yields

[TABLE]

B.2.3 Deleting an element

By setting $\widetilde{\gamma_{j}}=0$, we have $\kappa_{j}=\bm{\Sigma}_{jj}^{-1}$, and the deletion rules follow easily from the re-estimation rules:

[TABLE]