Sparse Bayesian Estimation of Parameters in Linear-Gaussian State-Space Models

Let us consider the additive linear-Gaussian state-space model (LGSSM), given by

for $t=1,\dots,T$, where ${\mathbf{x}}_{t}\in\mathbb{R}^{d_{x}}$ is the hidden state with associated observation $\mathbf{y}_{t}\in\mathbb{R}^{d_{y}}$ at time $t$, $\mathbf{A}\in\mathbb{R}^{d_{x}\times d_{x}}$ is the state transition matrix, $\mathbf{H}\in\mathbb{R}^{d_{y}\times d_{x}}$ is the observation matrix, $\mathbf{q}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{Q})$ is the state noise, and $\mathbf{r}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{R})$ is the observation noise. The state prior is $\mathbf{x}_{0}\sim\mathcal{N}(\bar{\mathbf{x}}_{0},\mathbf{P}_{0})$, with $\bar{\mathbf{x}}_{0}$ and $\mathbf{P}_{0}$ known. We assume that the model parameters remain fixed.

A common task in state-space modelling is the estimation of the series of $p({\mathbf{x}}_{t}|{\mathbf{y}}_{1:t})$ for $t\in\{1,\dots,T\}$, also known as the filtering distributions. In the case of the LGSSM, these distributions are obtained exactly via the Kalman filter equations [6, 7]. The linear-Gaussian assumption is not overly restrictive, as many systems can be approximated via linearisation, and for continuous problems Gaussian noises are very common.

Note that, the posterior distribution of any given parameter can be factorised as

where $\bm{\theta}$ is the parameter of interest, $p(\bm{\theta})$ is the prior ascribed to the parameter, and $p(\mathbf{y}_{1:T}|\bm{\theta})$ is extracted from the Kalman filter via the recursion

where $p(\mathbf{y}_{1}|\mathbf{y}_{1:0},\bm{\theta}):=p(\mathbf{y}_{1}|\bm{\theta})$ [6]. This factorisation gives the target distribution for estimating parameters in an LGSSM. In this work we focus on probabilistically estimating ${\mathbf{A}}$, and therefore we are interested in the posterior $p({\mathbf{A}}|\mathbf{y}_{1:T})\propto p(\mathbf{y}_{1:T}|{\mathbf{A}})p({\mathbf{A}})$.

In LGSSMs, ${\mathbf{H}}$ and ${\mathbf{R}}$ are frequently assumed to be known as parameters of the observation instrument, but ${\mathbf{Q}}$ and ${\mathbf{A}}$ are often unknown. For the purposes of this work, we assume that all parameters except ${\mathbf{A}}$ are known, or are suitably estimated, although our method can be extended to all parameters of the linear-Gaussian state-space model, as discussed in Section IV-E.

The estimation of the parameters of a state-space model is, in general, a difficult and computationally intensive task [6, 15]. This difficulty follows from the state dynamics not being directly observed, and stochasticity in the observations.

In this work we focus on Bayesian techniques, as our method is Bayesian, with this focus therefore allowing easier comparison. Frequentist methods for parameter estimation in SSMs are common however, with some relevant references being [29, 30, 31].

There are two main approaches to estimating and summarising parameters in state-space models, which we can broadly classify as point estimation methods and probabilistic methods.

A common method for point estimation in LGSSMs and in general is the expectation-maximisation (EM) algorithm [33], as explicit formulae exist for the conditional MLE of all parameters in the case of the LGSSM, and thus the model parameters can be estimated iteratively. The EM algorithm allows for all model parameters to be estimated simultaneously, but converges much more slowly as the number parameters to estimate increases [6]. Furthermore, this method does not allow for quantification of the uncertainty in the resultant estimates.

Note that there exist probabilistic methods, such as Laplace approximation [34] and variational inference [35], that provide analytical approximations and do not directly give samples. These methods are seldom used in state-space modelling. Distributional methods give more flexibility than point estimates, as they capture distributional behaviour and inherently quantify uncertainty, as well as allowing point estimates to be estimated from the samples, such as the aforementioned MAP estimator being the maximising argument for the posterior likelihood.

When fitting and designing statistical models, the presence of sparsity in parameters is often desirable, as it reduces the number of relevant variables thus easing interpretation and simplifying inference. Furthermore, real systems are often made up of several interacting dense blocks that, when taken as a whole, exhibit complex dynamics [16]. Sparse estimation methods allow for this structure to be recovered, resulting in estimates that can reflect the structure of the underlying system. Sparsity is ubiquitous within signal processing, with signal decomposition into a sparse combination of components being very common [36, 37, 38], which can be parameterised via model parameters. Furthermore, within signal processing, there exist a number of existing sparse Bayesian methods, such as [39, 40, 41, 42], although these do not operate within the paradigm of state-space modelling.

There are several approaches to estimate model parameters such that sparsity may be present. One approach may be to construct many models with unique combinations of sparse and dense elements, fit all of these models, and then select the best model according to some criteria (see [43, 44, 45] for examples). This approach is conceptually sound, but computationally expensive for even a small number of parameters $p$, as $2^{p}$ models must be fitted in order to obtain likelihood estimates, or other goodness-of-fit metrics. Another approach is to estimate the model parameters under a sparsity inducing penalty, with the classic example of such a penalty being the LASSO [46, 47]. This approach, commonly called regularisation, allows for only one model to be fitted rather than many, but increases the computational complexity of fitting the model. This single regularised estimate is, in most cases, more expensive to compute than fitting a non-regularised estimate as required by the previous approach, but this cost is typically far less than the cumulative cost of all required estimates in the previous approach. Regularisation is a common way to obtain sparse estimates, and can be extended to Bayesian modelling and estimation in the form of sparsity inducing priors [47, 48].

In LGSSMs, a sparse estimate of the transition matrix ${\mathbf{A}}$ can be interpreted as the adjacency matrix of a weighed directed graph $G$, with the nodes being the state elements, and the edges the corresponding elements of ${\mathbf{A}}$ [17, 22]. We illustrate this with an example in Fig. 1. We note that there exist a number of graph estimation methods that can be applied to time series data, such as [49, 50, 51, 52], although these methods do not utilise the structure of the state-space model. Furthermore, these methods often employ an acyclicity constraint, which prevents the results from exhibiting cycles, which are a common feature in real world dynamical systems, for example those resulting from discretised systems of ODEs.

The graph $G$ thus encodes the linear, between-step relationships of state elements, simplifying model interpretation. Under this graphical interpretation, $A_{ij}$ being non-zero implies that knowledge of the $j$th element of the state improves the prediction of the $i$th value. The presence of an edge from node $j$ to node $i$ implies a Granger-causal relationship between the state elements, as knowledge of the past values of the $j$th element at time $t$ improves the prediction of the $i$th element at time $t+1$, which is precisely the definition of a Granger-causal relationship [53].

Reversible jump Markov chain Monte Carlo (RJMCMC) was proposed as a method for Bayesian model selection [24], and has since seen use in fields such as ecology [54], Gaussian mixture modelling [27], and hidden Markov modelling [25]. RJMCMC has even been applied within the realm of signal processing, with some relevant references being [55, 56, 57]. However, RJMCMC has not been applied to the estimation the sparsity of model parameters within signal processing.

RJMCMC is an extension of the Metropolis-Hastings algorithm that allows for the sampling of a discrete model space, and thus the inclusion of many models within a single sampling chain. RJMCMC is a hierarchical sampler, with an upper layer sampling the models, and a lower layer sampling the posterior distribution of the parameters within the model. This hierarchy allows the use of standard MCMC methods for the lower layer, with the difficulty coming in designing the upper layer [27].

RJMCMC traverses the model space via transition kernels between pairs of models, with the jumps occurring probabilistically. This lends the model space an interpretation as a directed graph, with nodes representing the models and edges representing the jumps between models. Let $\Theta^{(i)}$ be the parameter space associated with model $M^{(i)}$, and $\bm{\theta}^{(i)}\in\Theta^{(i)}$ an associated realisation of the model parameters. Denote by $\pi_{i,j}$ the probability of jumping from model $M^{(i)}$ to $M^{(j)}$. Note that $\pi_{i,j}$ is zero if and only if $\pi_{j,i}$ is zero. Let $M^{(1)}$ be the current model, and $M^{(2)}$ be a candidate model. In order to construct a Markov kernel for the transition between models, a symmetry constraint is imposed, i.e., if it is possible to jump from $M^{(1)}$ to $M^{(2)}$, it must also be possible to jump from $M^{(2)}$ to $M^{(1)}$ [26]. In general however, the dimension of the parameter spaces is not equal, hence it is not possible to construct an invertible mapping between them, violating the required symmetry. Reversible jump MCMC addresses this by introducing a dimension matching condition [24]; the spaces $\Theta^{(1)}$ and $\Theta^{(2)}$ are augmented with simulated draws from selected distributions such that

where $T_{1,2}$ is a bijection, and $g_{i,j}(\cdot)$ are known distributions.

The parameter mappings and stochastic draws change the equilibrium distribution of the chain, which means that sampling will not be asymptotically correct. This counteracted by modifying the acceptance ratio; in the case of jumping from model $M^{(1)}$ to model $M^{(2)}$, the acceptance ratio is given by

where $p_{i}(\bm{\theta}^{(i)})$ is the density associated, with model $M^{(i)}$ evaluated at $\bm{\theta}^{(i)}$. The proposal ($M^{(2)}$, $\bm{\theta}^{(2)}$) is accepted with probability $\min(\alpha^{(1,2)},1)$, and is otherwise rejected. On rejection, the previous value of the chain is kept, ($M^{(1)}$, $\bm{\theta}^{(1)}$).

In this way, given data ${\mathbf{y}}$, RJMCMC samples the joint posterior $p(\bm{\theta},M|{\mathbf{y}})\propto p(M|{\mathbf{x}})p(\bm{\theta}|M,{\mathbf{y}})$. However, in the case only where only $\bm{\theta}$ is of interest, following standard Monte Carlo rules, we can discard the samples of $M$ to obtain $p(\bm{\theta}|{\mathbf{y}})$ [26].

Reversible jump MCMC incorporates many models into a single chain, so it is simple to compare or average models. However, the parameter mappings and model jump probabilities must be well designed. Poor selection of these parameters will typically lead to poor mixing in the model space [24, 26]. Our method explores only a single overall model, which simplifies the mappings. We impose a pairwise structure on the model space, simplifying the jumps significantly.

In order to use RJMCMC to explore sparsity in the transition matrix of a linear-Gaussian state-space model, we must first construct a set of candidate sub-models of the LGSSM that exhibit various sparsity levels. To this end we introduce the notation in Table I.

Denote by $M_{n}$ the model selected by the algorithm at iteration $n$. This model is uniquely defined by the associated set of indices of dense elements in ${\mathbf{A}}$, which we denote $\mathcal{M}_{n}$. Note that there are thus $2^{d_{x}^{2}}$ models in our model space, precluding parallel evaluation for even small $d_{x}$. It is therefore not possible to use methods such as Bayes factors or marginal likelihood to compare models, as is standard in Bayesian model selection in state-space models [43, 45], as the computational cost is infeasible, due to these requiring all models to be evaluated in order to be compared. The number of elements of $\mathcal{M}_{n}$, denoted by $D_{n}=|\mathcal{M}_{n}|$, is the number of dense elements at iteration $n$, and therefore the number of non-zero elements of ${\mathbf{A}}_{n}$. Denote by $S_{n}=d_{x}^{2}-D_{n}$ the number of sparse elements at iteration $n$. If the true value of a parameter is known, then it will be presented without subscript or superscript. We denote by a superscript $\mathsf{c}$ the complement to a set, and note that $\mathcal{M}^{\mathsf{c}}_{n}$ denotes the set of indices of elements sparse in ${\mathbf{A}}$ at iteration $n$.

Each model has an associated parameter space, which we denote by $\Theta_{n}$. As the parameter space of a sparse parameter is $\{0\}$, we therefore have

where $\bm{\theta}_{(i,j),n}$ is the support for $({\mathbf{A}}_{n})_{ij}$. In the Bayesian paradigm, we can interpret the sparsity of model $M_{n}$ as a prior constraint induced by $p({\mathbf{A}}|M_{n})$, under which the elements indexed by $\mathcal{M}_{n}^{\mathsf{c}}$ are always of value $0$, and the elements indexed by $\mathcal{M}_{n}$ are distributed as $p({\mathbf{A}})$. This is equivalent to fixing the value of the elements of ${\mathbf{A}}$ indexed by $\mathcal{M}_{n}^{\mathsf{c}}$ to $0$, as in Eq. (5).

We denote the $k\times k$ identity matrix by $\mathbf{Id}_{k}$, and the $k$-vector with all elements equal to $1$ by $\mathbf{1}_{k}$. We denote by $\cdot$ any unspecified parameters of a function or distribution. For example, $x\sim g(\cdot)$ means that the parameters of the distribution $g$ are unspecified, typically as they are irrelevant to the discussion.

We provide a table of our most used acronyms in Table II.

In order to explore potential sparsity of ${\mathbf{A}}$ using RJMCMC, we design a model jumping scheme that exploits the structure inherent to the model space.

Since our method applies MCMC to sample the posterior distribution $p({\mathbf{A}}|{\mathbf{y}}_{1:T})$, we must define a parameter proposal routine. Once a model has been proposed, the algorithm proposes a parameter value, ${\mathbf{A}}^{\prime}$, which is constrained to the parameter space of $M^{\prime}$. The process by which we generate the parameter proposal depends on whether or not the proposed model $M^{\prime}$ is the same as the previous model $M_{n-1}$. If $M^{\prime}=M_{n-1}$, the conditional posterior of the transition matrix, $p({\mathbf{A}}|{\mathbf{y}}_{1:T},M^{\prime})$, is sampled. Otherwise, the parameter value is mapped to the parameter space of $M^{\prime}$.

In this section, we first discuss the MH acceptance ratio, and then the way that prior knowledge of the transition matrix is incorporated, and how this relates to the model space. We then discuss the implications of sparsity in the samples.

Our method can be used to obtain sparse estimates of any of the parameters of the LGSSM, although some modifications are required to extend the formulation given above (for the transition matrix). Extending our method to the matrix ${\mathbf{H}}$ requires only for the parameter and model proposal to be changed to reflect the size of ${\mathbf{H}}$.

Extending the method to covariance matrices requires the proposal value to be constrained such that the resulting matrix is positive semi-definite. If the method does not jump models, remaining in model $M_{t-1}$, a possible proposal distribution for the state covariance proposal ${\mathbf{Q}}^{\prime}$ is

where $p$ is larger than $d_{x}$, with $p>30d_{x}$ working well experimentally. This proposal has expectation close to ${\mathbf{Q}}_{t-1}$, and is positive semi-definite. Note that this distribution is not applied to ${\mathbf{Q}}$ if estimating only ${\mathbf{A}}$, e.g., in Section III or in other subsections of Section IV, and applies only to the extension to sampling the covariance parameter. We enforce the sparsity structure of ${\mathbf{Q}}_{t-1}$ in ${\mathbf{Q}}^{\prime}$, by setting elements sparse under $M_{t-1}$ to $0$ after sampling but before the accept-reject step, replacing Step 2.1 in Algorithm 3, where the model does not change.

The model proposal step (Step 2.2 in Algorithm 3) would also need to be modified, with the diagonal being dense at all times, and enforcing indices $(i,j)$ and $(j,i)$ to have the same sparsity. The model space is therefore reduced, and model adjacency is assessed via only the upper triangular. This modification to the model proposal process completes the alterations required to use our method to sparsely sample the state covariance matrix. Note that the covariance parameters cannot be interpreted as encoding state connections, and therefore cannot be interpreted graphically in the same way as in the transition matrix.

The computational cost of our method is very similar to that of regular MCMC methods when applied to state-space models. The most computationally intensive component of the algorithm is the evaluation of the Kalman filtering equations, with this being over $95\%$ of the computational time in our testing. The additional costs compared to a random walk Metropolis-Hastings (RWMH) method are one or two draws from a uniform distribution, zero or one draw from a truncated Poisson distribution (equivalent to a categorical distribution), and some additional array accesses and comparisons. These extra costs are negligible compared to the cost of evaluating the Kalman filtering equations, with the computational cost and complexity being determined by the matrix operations therein, resulting in a complexity of $O(NT(d_{x}^{3}+d_{y}^{3}))$ for our algorithm, where the dominating $d_{x}^{3}$ and $d_{y}^{3}$ terms result from the matrix operations performed by Kalman filter. The computational cost of our method is empirically demonstrated in Section V, and is functionally equivalent to a standard MCMC method that does not explore sparsity. Thus in practice, given the cost of the filtering equations, the sparsity is explored for free.