Temporal Parallelization of Bayesian Smoothers

TL;DR

This paper introduces algorithms that enable the parallel processing of Bayesian smoothers over time, significantly reducing computational complexity from linear to logarithmic by leveraging all-prefix-sums operations.

Contribution

It formulates Bayesian smoothing as all-prefix-sums problems and develops parallel algorithms for filtering and smoothing, especially for linear/Gaussian models.

Findings

Reduces smoothing complexity from linear to logarithmic

Applies parallel scan algorithms to Bayesian smoothing

Specializes algorithms for linear/Gaussian models

Abstract

This paper presents algorithms for temporal parallelization of Bayesian smoothers. We define the elements and the operators to pose these problems as the solutions to all-prefix-sums operations for which efficient parallel scan-algorithms are available. We present the temporal parallelization of the general Bayesian filtering and smoothing equations and specialize them to linear/Gaussian models. The advantage of the proposed algorithms is that they reduce the linear complexity of standard smoothing algorithms with respect to time to logarithmic.

Equations

$$\begin{array}{rl}{x}_k& \sim p(x_k\mid x_{k-1}),\\ y_k& \sim p(y_k\mid x_k).\end{array}$$

$$p(x_k\mid y_{1:k-1})=\int p(x_k\mid x_{k-1})p(x_{k-1}\mid y_{1:k-1})\mathrm{d}{x}_{k-1},$$

$$p(x_k\mid y_{1:k})=\frac{p(y_k\mid x_k)p(x_k\mid y_{1:k-1})}{\int p(y_k\mid x_k)p(x_k\mid y_{1:k-1})\mathrm{d}{x}_k}.$$

$$\begin{array}{rl}& p(x_k\mid y_{1:n})\\ & =p(x_k\mid y_{1:k})\int \frac{p(x_{k+1}\mid x_k)p(x_{k+1}\mid y_{1:n})}{p(x_{k+1}\mid y_{1:k})}\mathrm{d}{x}_{k+1}.\end{array}$$

$$(a_1,a_1\otimes a_2,\dots ,a_1\otimes \cdots \otimes a_n).$$

$$\mathcal{F}=\left\{\left(f,g\right):\int f\left(y\mid z\right)\mathrm{d}y=1\right\},$$

$$\left(f_i,g_i\right)\otimes \left(f_j,g_j\right)$$

$$f_{ij}\left(x\mid z\right)$$

$$g_{ij}\left(z\right)$$

$$f_k\left(x_k\mid x_{k-1}\right)$$

$$g_k\left(x_{k-1}\right)$$

$$a_1\otimes a_2\otimes \cdots \otimes a_k$$

$$\mathcal{S}=\left\{a:\int a\left(x\mid z\right)\mathrm{d}x=1\right\}.$$

$$a_i\otimes a_j$$

$$a_{ij}\left(x\mid z\right)$$

$$a_k\otimes a_{k+1}\otimes \cdots \otimes a_n$$

$$p\left(y_{1:n}\right)$$

$$p\left(y_k\mid y_{1:k-1}\right)$$

$$a_k$$

$$a_k=p\left(x_{lk}\mid y_{1:l\left(k+1\right)-1},x_{l\left(k+1\right)}\right)$$

$$x_k$$

$$y_k$$

$$p\left(x_k\mid x_{k-1}\right)$$

$$p\left(y_k\mid x_k\right)$$

$$f_k\left(x_k\mid x_{k-1}\right)$$

$$g_k\left(x_{k-1}\right)$$

$$\begin{array}{rl}{A}_k& =\left(I_{n_x}-K_k{H}_k\right)F_{k-1},\\ b_k& =u_{k-1}+K_k\left(y_k-H_k{u}_{k-1}-d_k\right),\\ C_k& =\left(I_{n_x}-K_k{H}_k\right)Q_{k-1},\\ K_k& =Q_{k-1}{H}_k^{\top }{S}_k^{-1},\\ S_k& =H_k{Q}_{k-1}{H}_k^{\top }+R_k,\end{array}$$

$$\begin{array}{rl}{m}_1^{-}& =F_0{m}_0+u_0,\\ P_1^{-}& =F_0{P}_0{F}_0^{\top }+Q_0,\\ S_1& =H_1{P}_1^{-}{H}_1^{\top }+R_1,\\ K_1& =P_1^{-}{H}_1^{\top }{S}_1^{-1},\\ A_1& =0,\\ b_1& =m_1^{-}+K_1[y_1-H_1{m}_1^{-}-d_1],\\ C_1& =P_1^{-}-K_1{S}_1{K}_1^{\top }.\end{array}$$

$$\begin{array}{rl}{\eta }_k& =F_{k-1}^{\top }{H}_k^{\top }{S}_k^{-1}\left(y_k-H_k{u}_{k-1}-d_k\right),\\ J_k& =F_{k-1}^{\top }{H}_k^{\top }{S}_k^{-1}{H}_k{F}_{k-1},\end{array}$$

$$f_i\left(y\mid z\right)$$

$$g_i\left(z\right)$$

$$f_j\left(y\mid z\right)$$

$$g_j\left(z\right)$$

$$\left(f_i,g_i\right)\otimes \left(f_j,g_j\right)$$

Full text

{IEEEkeywords}

Bayesian smoothing, Kalman filtering and smoothing, parallel computing, parallel scan, prefix sums

1 Introduction

Parallel computing is rapidly transforming from a scientists’ computational tool to a general purpose computational paradigm. The availability of affordable massively-parallel graphics processing units (GPUs) as well as widely-available parallel grid and cloud computing systems [1, 2, 3] drive this transformation by bringing parallel computing technology to everyday use. This creates a demand for parallel algorithms that can harness the full power of the parallel computing hardware.

Stochastic state-space models allow for modeling of time-behaviour and uncertainties of dynamic systems, and they have long been used in various tracking, automation, communications, and imaging applications [4, 5, 6, 7, 8]. More recently, they have also been used as representations of prior information in machine learning setting (see, e.g., [9]). In all of these applications, the main problem can be mathematically formulated as a state-estimation problem on the stochastic model, where we estimate the unknown phenomenon from a set of noisy measurement data. Given the mathematical problem, the remaining task is to design efficient computational methods for solving the inference problems on large data sets such that they utilize the computational hardware as effectively as possible.

Bayesian filtering and smoothing methods [6] provide the classical [10] solutions to state-estimation problems which are computationally optimal in the sense that their computational complexities are linear with respect to the number of data points. Although these solutions are optimal for single central processing unit (CPU) systems, due to the sequential nature of the algorithms, their complexity remains linear also in parallel multi-CPU systems. However, genuine parallel algorithms are often capable to perform operations in logarithmic number of steps by massive parallelization of the operations. More precisely, their span-complexity [3], that is, the number of computational steps as measured by a wall-clock, is often logarithmic with respect to the number of data points. However, their total number of operations, the work-complexity, is still linear as all the data points need to be processed.

Despite the long history of state-estimation methods, the existing parallelization methods have concentrated on parallelizing the subproblems arising in Bayesian filtering and smoothing methods, but there is a lack of algorithms that are specifically designed for solving state-estimation problems in parallel architectures. There are, however, some existing approaches for parallelizing Kalman type of filters as well as particle filters. One approach was studied in [11] and [12] is to parallelize the corresponding batch formulation, which leads to sub-linear computational methods, because the matrix computations can be parallelized. If the state-space of the Kalman filter is large, it is then possible to speed up the matrix computations via parallelization [13, 14]. Particle filters can also be parallelized over the particles [15, 16] the bottleneck being the resampling step. In some specific cases such as in multiple target tracking [17] it is possible to develop parallelized algorithms by using the structure of the specific problem.

The contribution of this article is to propose a novel general algorithmic framework for parallel computation of Bayesian smoothing solutions for state space models. We also present algorithms for parallelizing the Bayesian filtering solutions, but our focus is in smoothing, because the parallel computation is done off-line in the sense that all the data needs to be available during the parallel computations and it cannot arrive sequentially. Our approach to parallelization differs from the aforementioned approaches in the aspect that we replace the whole Bayesian filtering and smoothing formalism with another, parallelizable formalism. We replace the Bayesian filtering and smoothing equations [4, 6] with another set of equations that can be combined with so-called scan or prefix-sums algorithm [18, 19, 20, 2], which is one of the fundamental algorithm frameworks in parallel computing. The advantage of this is that it allows for reduction of the linear $O(n)$ complexity of batch filtering and smoothing algorithms to logarithmic $O(\log n)$ span-complexity in the number $n$ of data points. Based on the novel formulation we develop parallel algorithms for computing the filtering and smoothing solutions to linear Gaussian systems with the logarithmic span-complexity. As this parallelization is done in temporal direction, the individual steps of the resulting algorithm could further be parallelized in the same way as Kalman filters and particle filters have previously been parallelized [13, 17, 14, 15, 16].

The organization of the article is the following. In Section 2 we review the classical Bayesian filtering and smoothing methodology as well as the parallel scan algorithm for computing prefix sums. In Section 3 we present the general framework for parallelizing Bayesian filtering and smoothing methods. Section 4 is concerned with specializing the general framework to linear Gaussian systems. Numerical example of linear/Gaussian systems is given in Section 5 and finally Section 6 concludes the article along with discussion on various aspects of the methodology.

2 Background

2.1 Bayesian filtering and smoothing

Bayesian filtering and smoothing methods [6] are algorithms for statistical inference in probabilistic state-space models of the form

[TABLE]

with $x_{0}\sim p(x_{0})$. Above, the state $x_{k}\in\mathbb{R}^{n_{x}}$ at time step $k$ evolves as a Markov process with transition density $p(x_{k}\mid x_{k-1})$. State $x_{k}$ is observed by the measurement $y_{k}\in\mathbb{R}^{n_{y}}$ whose density is $p(y_{k}\mid x_{k})$.

The objective of Bayesian filtering is to compute the posterior density $p(x_{k}\mid y_{1:k})$ of the state $x_{k}$ given the measurements $y_{1:k}=(y_{1},\ldots,y_{k})$ up to time step $k$. Given the measurements up to a time step $n$, the objective is smoothing is to compute the density $p(x_{k}\mid y_{1:n})$ of the state $x_{k}$ for $k<n$.

The key insight of Bayesian filters and smoothers is that the computation of the required densities can be done in linear $O(n)$ number of computational steps by using recursive (forward) filtering and (backward) smoothing algorithms. This is significant, because a naive computation of the posterior distribution would typically take at least $O(n^{3})$ computational steps.

The Bayesian filter is a sequential algorithm, which iterates the following prediction and update steps:

[TABLE]

[TABLE]

Given the filtering outputs for $k=1,\ldots,n$, the Bayesian forward-backward smoother consists of the following backward iteration for $k=n-1,\ldots,1$:

[TABLE]

When applied to a batch of data of size $n$, the computational complexity of the filter and smoother is $O(n)$ as they perform $n$ sequential steps in forward and backward directions when looping over the data. The Kalman filter and Rauch–Tung–Striebel (RTS) smoother [21, 22] are the solutions to these recursions when the transition densities are linear and Gaussian. The filtering and smoothing equations can also be analogously solved in closed form for discrete-state models [7]. In this paper, we show how to parallelize the previous recursions using parallel scan-algorithm, which is reviewed next.

2.2 Parallel scan-algorithm

The parallel-scan algorithm [20] is a general parallel computing framework that can be used to convert sequential $O(n)$ algorithms with certain associative property to $O(\log n)$ parallel algorithms. The algorithm was originally developed for computing prefix sums [18], where it uses the associative property of summation. The algorithm has since been generalized to arbitrary associative binary operators and it is used as the basis of multitude of parallel algorithms including sorting, linear programming, and graph algorithms. This kind of algorithms are especially useful in GPU-based computing systems and they are likely to be fundamental algorithms in a many future parallel computing systems.

The problem that the parallel-scan algorithm [20] solves is the all-prefix-sums operation, which is defined next.

Definition 1

Given a sequence of elements $(a_{1},a_{2},\ldots,a_{n})$, where $a_{i}$ belongs to a certain set, along with an associative binary operator $\otimes$ on this set, the all-prefix-sums operation computes the sequence

[TABLE]

For example, if we have $n=4$, $a_{i}=i$, and $\otimes$ denotes the ordinary summation, the all-prefix-sums are $(1,3,6,10)$. If $\otimes$ denotes the subtraction, the all-prefix-sums are $(1,-1,-4,-8)$. It should be noted that the operator is not necessarily commutative so we use a product symbol, as matrix products are not commutative, instead of a summation symbol.

The all prefix-sums operation can be computed sequentially by processing one element after the other. However, this direct sequential iteration inherently takes $O(n)$ time. We can now see the analogy of the iteration to the Bayesian filter discussed in previous section – both of the algorithms have linear $O(n)$ complexity, because they need to loop over all the elements in forward direction. A similar argument applies to the Bayesian smoothing pass.

Fortunately, the all-prefix-sum sums operation can be computed in parallel in $O(\log n)$ span-time by using up-sweep and down-sweep algorithms [20] shown in Fig. 1. These algorithms correspond to up and down traversals in a binary tree which are used for computing partial (generalized) sums of the elements. A final pass is then used to construct the final result. The algorithms can be used for computing all-prefix-sums (5) for an arbitrary associative operator $\otimes$.

3 Parallel Bayesian filtering and smoothing

In this section, we explain how to define the elements and the binary operators to be able to perform Bayesian filtering and smoothing using parallel scan algorithms.

3.1 Bayesian filtering

In order to perform parallel Bayesian filtering, we need to find the suitable element $a_{k}$ and the binary associative operator $\otimes$. As we will see in this section, an element $a$ consists of a pair $\left(f,g\right)\in\mathcal{F}$ where $\mathcal{F}$ is

[TABLE]

and $f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{x}}\rightarrow\left[0,\infty\right)$ represents a conditional density, and $g:\mathbb{R}^{n_{x}}\rightarrow\left[0,\infty\right)$ represents a likelihood.

Definition 2

Given two elements $\left(f_{i},g_{i}\right)\in\mathcal{F}$ and $\left(f_{j},g_{j}\right)\in\mathcal{F}$, the binary associative operator $\otimes$ for Bayesian filtering is

[TABLE]

where

[TABLE]

The proof that $\otimes$ has the associative property is given in Appendix 7.

Theorem 3

Given the element $a_{k}=\left(f_{k},g_{k}\right)\in\mathcal{F}$ where

[TABLE]

$p\left(x_{1}\mid y_{1},x_{0}\right)=p\left(x_{1}\mid y_{1}\right)$, and $p\left(y_{1}\mid x_{0}\right)=p\left(y_{1}\right)$, the $k$-th prefix sum is

[TABLE]

Theorem 3 is proved in Appendix 7. Theorem 3 implies that we can parallelise the computation of all filtering distributions $p\left(x_{k}\mid y_{1:k}\right)$ and the marginal likelihoods $p\left(y_{1:k}\right)$, of which the latter ones can be used for parameter estimation [6].

Remark 4

If we only know $p\left(y_{k}\mid x_{k-1}\right)$ up to a proportionality constant, which means that $g_{k}\left(x_{k-1}\right)\propto p\left(y_{k}\mid x_{k-1}\right)$, we can still recover the filtering density $p\left(x_{k}\mid y_{1:k}\right)$ by the above operations. However, we will not be able to recover the marginal likelihoods $p\left(y_{1:k}\right)$. We can nevertheless recover $p\left(y_{1:k}\right)$ by an additional parallel pass, as will be explained in Section 3.3.

3.2 Bayesian smoothing

The Bayesian smoothing pass requires that the filtering densities have been obtained beforehand. In smoothing, we consider a different type of element $a$ and binary operator $\otimes$ than those used in filtering. As we will see in this section, an element $a$ is a function $a:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{x}}\rightarrow\left[0,\infty\right)$ that belongs to the set

[TABLE]

Definition 5

Given two elements $a_{i}\in\mathcal{S}$ and $a_{j}\in\mathcal{S}$, the binary associative operator $\otimes$ for Bayesian smoothing is

[TABLE]

where

[TABLE]

The proof that $\otimes$ has the associative property is included in Appendix 8.

Theorem 6

Given the element $a_{k}=p\left(x_{k}\mid y_{1:k},x_{k+1}\right)\in\mathcal{S}$, with $a_{n}=p\left(x_{n}\mid y_{1:n}\right)$, we have that

[TABLE]

Theorem 6 is proved in Appendix 8. Theorem 6 implies that we can compute all smoothing distributions in parallel form. However, it should be noted we should apply the parallel scan algorithm with elements in reverse other, that is, with elements $b_{k}=a_{n-k+1}$, so that the prefix-sums $b_{1}\otimes\cdots\otimes b_{k}$ recover the smoothing densities.

3.3 Additional aspects

We proceed to discuss additional aspects of the previous formulation of filtering and smoothing. In Section 3.1, it was indicated that the marginal likelihood $p\left(y_{1:n}\right)$ is directly available from the parallel scan algorithm if $g_{k}\left(x_{k-1}\right)=p\left(y_{k}\mid x_{k-1}\right)$. However, sometimes we only know $p\left(y_{k}\mid x_{k-1}\right)$ up to a proportionality constant so $g_{k}\left(x_{k-1}\right)\propto p\left(y_{k}\mid x_{k-1}\right)$, as will happen in Section 4. Although in this case, the parallel scan Bayesian filtering algorithm provides us with the filtering densities but not the marginal likelihood $p\left(y_{1:n}\right)$, we can still recover the marginal likelihoods as follows. We first run the parallel filtering algorithm to recover all filtering distributions $p\left(x_{k}\mid y_{1:k}\right)$ for $k=1$ to $n$ and then, we perform the following decomposition for $p\left(y_{1:n}\right)$

[TABLE]

where

[TABLE]

Each factor $p\left(y_{k}\mid y_{1:k-1}\right)$ can be computed in parallel using the predictive density $p\left(x_{k}\mid y_{1:k-1}\right)$ and the likelihood $p\left(y_{k}\mid x_{k}\right)$. We can then recover all $p\left(y_{1:k}\right)$ by $O(\log n)$ parallel recursive pairwise multiplications of the adjacent terms.

It is also possible to perform the parallelization at block level instead of at individual element level. When using the parallel scan algorithm, we do not need to assign each single-measurement element to a single computational node, but instead we can perform initial computations in blocks such that a single node processes a block of measurements before combining the results with other blocks. The results of the blocks can then be used as the elements in the parallel-scan algorithm. This kind of procedure corresponds to selecting the elements for the scan algorithm to consist of blocks of length $l$:

[TABLE]

in filtering and

[TABLE]

in smoothing instead of the corresponding terms with $l=1$. A practical advantage of this is that we can more easily distribute the computations to a limited number of computational nodes while still getting the optimal speedup from parallelization.

4 Parallel linear/Gaussian filter and smoother

The parallel linear/Gaussian filter and smoother are obtained by particularising the element $a$ and binary operator $\otimes$ for Bayesian filtering and smoothing explained in the previous section to linear/Gaussian systems. The sequential versions of these algorithms correspond to the Kalman filter and the RTS smoother.

We consider the linear/Gaussian state space model

[TABLE]

where $F_{k-1}\in\mathbb{R}^{n_{x}\times n_{x}}$ and $H_{k}\in\mathbb{R}^{n_{y}\times n_{x}}$ are known matrices, $u_{k-1}\in\mathbb{R}^{n_{x}}$ and $d_{k}\in\mathbb{R}^{n_{y}}$ are known vectors, and $q_{k-1}$ and $r_{k}$ are zero-mean, independent Gaussian noises with covariance matrices $Q_{k-1}\in\mathbb{R}^{n_{x}\times n_{x}}$ and $R_{k}\in\mathbb{R}^{n_{y}\times n_{y}}$. The initial distribution is given as $x_{0}\sim\mathrm{N}(m_{0},P_{0})$. With this model, we have that

[TABLE]

In this section, we use the notation $\mathrm{N}_{I}\left(\cdot;\eta,J\right)$ to denote a Gaussian density parameterised in information form so that $\eta$ is the information vector and $J$ is the information matrix. If a Gaussian distribution has mean $\overline{x}$ and covariance matrix $P$, its parameterisation in information form is $\eta=P^{-1}\overline{x}$ and $J=P^{-1}$. This parametrization corresponds to so-called information form of Kalman filter [23]. We also use $I_{n_{x}}$ to denote an identity matrix of size $n_{x}$.

4.1 Linear/Gaussian filtering

We first describe the representation of an element $a_{k}\in\mathcal{F}$ for filtering in linear and Gaussian systems by the following lemma.

Lemma 7

For linear/Gaussian systems, the element $a_{k}\in\mathcal{F}$ for filtering becomes

[TABLE]

where the parameters of the first term are given for $k>1$ as

[TABLE]

and for $k=1$ as

[TABLE]

The parameters of the second term are given as

[TABLE]

for $k=1,\ldots,n$.

In Lemma 7, densities $p\left(x_{k}\mid y_{k},x_{k-1}\right)$ and $p\left(y_{k}\mid x_{k-1}\right)$ are obtained by applying the Kalman filter update with measurement $y_{k}$, distributed according to (9), applied to the density $p\left(x_{k}\mid x_{k-1}\right)$ in (8) and matching the terms. For the first step we have applied the Kalman filter prediction and update steps starting from $x_{0}\sim\mathcal{N}(m_{0},P_{0})$ and matched the terms.

Therefore, an element $a_{k}$ can be parameterised by $\left(A_{k},b_{k},C_{k},\eta_{k},J_{k}\right)$, which can be computed for each element in parallel. Also, it is relevant to notice that if the system parameters $\left(F_{k},u_{k},Q_{k},H_{k},d_{k},R_{k}\right)$ do not depend on the time step $k$, the only parameters of $a_{k}$ that depend on $k$ are $b_{k}$ and $\eta_{k}$, as they depend on the measurement $y_{k}$.

Lemma 8

Given two elements $\left(f_{i},g_{i}\right)\in\mathcal{F}$ and $\left(f_{j},g_{j}\right)\in\mathcal{F}$, with parameterisations

[TABLE]

the binary operator $\otimes$ for filtering becomes

[TABLE]

where

[TABLE]

with

[TABLE]

The proof is provided in Appendix 9.

4.2 Linear/Gaussian smoothing

We first describe the representation of an element $a_{k}\in\mathcal{S}$ for smoothing in linear and Gaussian systems by the following lemma.

Lemma 9

For linear/Gaussian systems, the element $a_{k}\in\mathcal{S}$ for smoothing becomes

[TABLE]

where for $k<n$

[TABLE]

and for $k=n$ we have

[TABLE]

Above, $\overline{x}_{k}$ and $P_{k}$ are the filtering mean and covariance matrix at time step $k$, such that $p\left(x_{k}\mid y_{1:k}\right)=\mathrm{N}\left(x_{k};\overline{x}_{k},P_{k}\right)$.

Lemma 9 is obtained by performing a Kalman filter update on density $p\left(x_{k}\mid y_{1:k}\right)$ with an observation $x_{k+1}$, whose distribution is given by (8). Element $a_{k}$ for smoothing with linear/Gaussian systems can be parameterised as $a_{k}=\left(E_{k},g_{k},L_{k}\right)$.

Lemma 10

Given two elements $a_{i}\in\mathcal{S}$ and $a_{j}\in\mathcal{S}$ with parameterisation

[TABLE]

the binary operator $\otimes$ for smoothing becomes

[TABLE]

where

[TABLE]

and

[TABLE]

5 Numerical experiment

In order to illustrate the benefit of parallelization we consider a simple tracking model (see, e.g., [5, 6]) with the state $x=\begin{pmatrix}u&v&\dot{u}&\dot{v}\end{pmatrix}^{\top}$, where $(u,v)$ is the 2D position and $(\dot{u},\dot{v})$ is the 2D velocity of the tracked object. From noisy measurements of the position $(u,v)$, we aim to solve the smoothing problem in order to determine the whole trajectory of the target.

The model has the form

[TABLE]

where $q_{k}\sim\mathrm{N}(0,Q)$, $r_{k}\sim\mathrm{N}(0,R)$, and

[TABLE]

along with

[TABLE]

In our simulations we used the parameters $\sigma=0.5$, $\Delta t=0.1$, $q=1$, and started the trajectory from a random Gaussian initial condition with mean $m_{0}=\begin{pmatrix}0&0&1&-1\end{pmatrix}^{\top}$ and covariance $P_{0}=I_{4}$.

Fig. 2 shows a typical trajectory and measurements from the model defined by Eqs. (15) and (16) along with the Kalman filter and RTS smoother solutions. As the parallel algorithms produce exactly the same filter and smoothing solutions as the classic sequential algorithms, this result also illustrates the typical result produced by the proposed algorithms.

We now aim to evaluate the required number of floating point operations (flops) for generating the smoothing solution for this model. In order to do that, we run the sequential filter and smoothing methods (KF and RTS) as well as the proposed parallel algorithms (PKF and PRTS) over simulated data sets of different sizes and evaluate their span and work flops. The span flops here refers to the minimum number of floating point steps when the parallelizable operations in the algorithm are done in parallel – this corresponds to the actual execution time required to do the computations in a parallel computer. The work flops refers to the total number of operations that the parallel computer needs to perform – it measures the total energy required for the computations or equivalently the time required by the algorithm in a single-core computer. As the classic sequential KF and RTS algorithms are not parallelizable, their span and work flops are equal. The flops have been computed by estimating how many flops each of the matrix operations takes (multiplication, summation, LU-factorization) and incrementing the flops counter after every operation in the code.

Fig. 3 shows the flops required by the sequential KF along with the span flops and work flops required by the parallel Kalman filter algorithm. As expected, with small data set sizes the number of span flops required by the parallel KF is larger than that of the sequential KF, but already starting from time step count of around 20, the span flops is lower for the parallel KF. The logarithmic growth of the span flops in the parallel algorithm can be clearly seen while the number of flops for the sequential KF grows linearly. However, the work flops required by the parallel KF is approximately 8 times the flops of the sequential KF. This means that although the execution time for the parallel algorithms is smaller than for the sequential algorithms, they need to perform more floating point operations in total.

The flops required by the sequential RTS smoother along with the span flops and work flops required by the parallel RTS smoother are shown in Fig. 4. In this case, the parallel algorithm reaches the sequential algorithm speed already with data set of size less than 10. Furthermore, the total number of floating point operations required by the parallel algorithm is approximately 4 times the operations required by the sequential algorithm. The ratios of these total (work) operations for both the filter and smoother are shown in Fig. 5.

6 Conclusion and Discussion

In this article we have proposed a novel general algorithmic framework for parallel computation of batch Bayesian filtering and smoothing solutions for state space models. The framework is based on formulating the computations in terms of associative operations between suitably defined elements such that the all-prefix-sums operation computed by a parallel-scan algorithm exactly produces the Bayesian filtering and smoothing solutions. The advantage of the framework is that the parallelization allows for performing the computations in $O(\log n)$ span complexity, where $n$ is the number of data points, while sequential filtering and smoothing algorithms have an $O(n)$ complexity. Parallel versions of Kalman filters and Rauch–Tung–Striebel smoothers were derived as special cases of the framework. The computational advantages of the framework were illustrated in a numerical simulation.

A disadvantage of the proposed methodology is that although the wall-clock time of execution is significantly reduced, the total amount of computations (and hence required energy) is larger than with conventional sequential algorithms. Although the total amount of computations is only increased by a constant factor, in some systems, such as small-scale mobile systems, even if parallelization would be possible, it can be beneficial to use the classic algorithms. However, the speedup gain of the proposed approach is beneficial in applications such as data-assimilation based weather forecasting [24] and other spatio-temporal systems appearing, for example, in tomographic reconstruction [8] or machine learning [25], where the computations take a significant amount of time. In these systems, it is possible to dedicate the required amount of extra computational resources to gain the significant speedup provided by parallelization.

Although we have restricted our consideration to specific types parallel-scan algorithms, it is also possible to use other kinds of algorithms for computing the prefix sums corresponding to the Bayesian filtering and smoothing solutions. We could also select algorithms for given computer or network architectures, for minimizing the communication between the nodes, or for minimizing the energy consumption [26, 27]. The present formulation of the computations in terms of local associative operations is likely to have other applications beyond parallelization. For example, in decentralized systems, it is advantageous to be able to first perform operations locally and then combine them to produce the full state-estimation solution.

The proposed framework is also valid for discrete state spaces as well as for other state spaces provided that we consider the elements with the appropriate domain and replace the Lebesgue integrals by integrals with respect to the corresponding reference measure, e.g., counting measure in the case of discrete-state models.

The framework could be extended to non-linear and non-Gaussian models by replacing the exact Kalman filters and smoothers with iterated extended Kalman filters and smoothers [28, 29] or their sigma-point/numerical-integration versions such as posterior linearization filters and smoothers [30, 31, 32]. Possible future work also includes developing particle filter and smoother methods (see, e.g., [6]) for the present framework along with various other Bayesian filter and smoother approximations proposed in literature.

\appendices

7

In this appendix, we prove the required results for Bayesian filtering: the associative property of the operator in Definition 2 and Theorem 3.

7.1 Associative property

In order to prove the associative property of $\otimes$ for filtering, we need to prove that for three elements $\left(f_{i},g_{i}\right)$, $\left(f_{j},g_{j}\right)$, $\left(f_{k},g_{k}\right)\in\mathcal{F}$, the following relation holds

[TABLE]

We proceed to perform the calculations on both sides of the equation to check that they yield the same result.

7.1.1 Left-hand side

We use Definition 2 in the left-hand side of (18) and obtain

[TABLE]

where

[TABLE]

and

[TABLE]

7.1.2 Right-hand side

We first use operator $\otimes$ to the elements with indices $j$ and $k$ in the right-hand side of (18), see Definition 2,

[TABLE]

where

[TABLE]

Then, the right-hand side of (18) becomes

[TABLE]

where

[TABLE]

and

[TABLE]

In (19), (20), it is met that $f^{\prime}_{ijk}\left(x\mid z\right)=f{}_{ijk}\left(x\mid z\right)$ and $g^{\prime}_{ijk}\left(z\right)=g{}_{ijk}\left(z\right)$, which proves the associative property of $\otimes$ in Definition 2.

7.2 Proof of Theorem 3

In this appendix, we prove Theorem 3. We first prove by induction that

[TABLE]

for $l<k+1$. Relation (23) holds for $l=0$ by definition of $a_{k}$. Then, assuming that

[TABLE]

holds, we need to prove that (23) holds.

We calculate the first element of $a_{k-l}\otimes b_{k-l+1}$, denoted by $f_{ab}$, where $b_{k-l+1}=a_{k-l+1}\otimes\cdots\otimes a_{k-1}\otimes a_{k}$. We have

[TABLE]

Function $f_{ab}$ corresponds to the first element of (23), as required. We further get

[TABLE]

Function $g_{ab}$ corresponds to the second element of (23), as required.

Substituting $l=k+2$ into (23), we obtain

[TABLE]

We now calculate the first element of $a_{1}\otimes\left[a_{2}\otimes\cdots\otimes a_{k}\right]$, denoted as $f_{1k}$, where $a_{1}$ is given in Theorem 3:

[TABLE]

The second element of $a_{1}\otimes\left[a_{2}\otimes\cdots\otimes a_{k}\right]$, denoted as $g_{1k}$ is

[TABLE]

Results (26) and (27) finish the proof of Theorem 3.

8

In this appendix, we prove the required results for Bayesian smoothing: the associative property of the operator in Definition 5 and Theorem 6.

8.1 Associative property

In order to prove the associative property of $\otimes$ for filtering, we need to prove that, for three elements $a_{i}$, $a_{j}$, $a_{k}\in\mathcal{S}$ , the following relation holds:

[TABLE]

We proceed to perform the calculations on both sides of the equation to check that they yield the same result.

8.1.1 Left-hand side

We apply the operator in Definition 5 on the left-hand side of (28) to obtain

[TABLE]

where

[TABLE]

8.1.2 Right-hand side

We first calculate $a_{j}\otimes a_{k}$ using the operator in Definition 5 which gives

[TABLE]

Then, we calculate the right hand side of (28) we have that

[TABLE]

where

[TABLE]

We can see that $a_{ijk}$ in (29) is equal to $a^{\prime}_{ijk}$ in (30), which proves the associative property of the operator in Definition 5.

8.2 Proof of Theorem 6

In this appendix, we prove Theorem 6. We first prove by induction that

[TABLE]

for $l<n-k$. Relation (31) holds for $l=0$ by definition of $a_{k}$. Then, assuming that

[TABLE]

holds, we need to prove that (23) holds.

We use $a_{k+l}$ in Theorem 6 to calculate

[TABLE]

This proves (31).

If $l=n-k-1$ and $a_{n}$ as in Theorem 6, we have

[TABLE]

This result finishes the proof of Theorem 6.

9

In this appendix, we prove Lemma 8. We have the following easily verifiable identities:

[TABLE]

and

[TABLE]

We also have

[TABLE]

By using Definition 2 for $f_{ij}$ and $g_{ij}$ together with parameterizations in Lemma 8, elementary computations lead to (13) and (14).

Acknowledgment

The authors would like to thank Academy of Finland for financial support.