Iterated Extended Kalman Smoother-based Variable Splitting for $L_1$-Regularized State Estimation

TL;DR

This paper introduces a novel framework combining iterated extended Kalman smoothers with variable splitting techniques like ADMM to efficiently solve large-scale $L_1$-regularized state estimation problems, offering low complexity and promising results.

Contribution

It develops a new algorithmic framework that integrates Kalman smoothing with variable splitting for $L_1$-regularized state estimation, enabling efficient large-scale problem solving.

Findings

The proposed methods achieve significant speed-ups.

The algorithms demonstrate convergence guarantees.

Experiments show effective handling of high-dimensional problems.

Abstract

In this paper, we propose a new framework for solving state estimation problems with an additional sparsity-promoting $L_1$-regularizer term. We first formulate such problems as minimization of the sum of linear or nonlinear quadratic error terms and an extra regularizer, and then present novel algorithms which solve the linear and nonlinear cases. The methods are based on a combination of the iterated extended Kalman smoother and variable splitting techniques such as alternating direction method of multipliers (ADMM). We present a general algorithmic framework for variable splitting methods, where the iterative steps involving minimization of the nonlinear quadratic terms can be computed efficiently by iterated smoothing. Due to the use of state estimation algorithms, the proposed framework has a low per-iteration time complexity, which makes it suitable for solving a large-scale or high-dimensional state estimation problem. We also provide convergence results for the proposed algorithms. The experiments show the promising performance and speed-ups provided by the methods.

Figures

Tables

Table 1. TABLE I: Different choices of IEKS-based variable splitting algorithms

Method	Related Quadratic Term	${𝚯}_t$	${𝚫}_t$	${𝚺}_t$
PRS	$\frac{\rho }{2}{\Vert {𝐰}_t-{𝛀}_t{𝐱}_t+{𝜼}_t/\rho \Vert }^2$	${𝛀}_t$	${𝐰}_t+{𝜼}_t/\rho$	$𝐈/\rho$
SBM	$\frac{\rho }{2}{\Vert {𝐰}_t-{𝛀}_t{𝐱}_t+{𝜼}_t\Vert }^2$	${𝛀}_t$	${𝜼}_t+{𝐰}_t$	$\rho 𝐈$
FOPD	$\frac{1}{2\rho }{\Vert {𝐱}_t-\left({𝐱}_t^{\left(k\right)}-{𝛀}_t^{\top }{𝐰}_t\right)\Vert }^2$	$𝐈$	${𝐱}_t^{\left(k\right)}-{𝛀}_t^{\top }{𝐰}_t$	$\rho 𝐈$
ADMM	$\frac{\rho }{2}{\Vert {𝐰}_t-{𝛀}_t{𝐱}_t+{𝜼}_t/\rho \Vert }^2$	${𝛀}_t$	${𝐰}_t+{𝜼}_t/\rho$	$𝐈/\rho$

Table 2. TABLE II: Comparison of average CPU time (seconds) with different λ 𝜆 \lambda in the linear case.

$\lambda$	$T$	PRS	SBM	FOPD	ADMM	KS-PRS	KS-SBM	KS-FOPD	KS-ADMM
0.1	${10}^3$	$6.05$	25.42	14.06	6.12	0.23	0.16	0.31	0.16
	${10}^4$	841	6210	2379	851	0.54	0.60	1.55	0.53
	${10}^5$	-	-	-	-	28.1	30.1	11.2	9.7
	${10}^6$	-	-	-	-	39.4	52.1	87.8	38.5
	${10}^7$	-	-	-	-	402	922	341	337
	${10}^8$	-	-	-	-	3330	4121	3510	3378
0.5	${10}^3$	6.07	24.61	14.04	6.06	0.22	0.17	0.32	0.14
	${10}^4$	832	6178	2366	837	0.52	0.60	$1.53$	0.50
	${10}^5$	-	-	-	-	27.9	30.0	10.9	9.4
	${10}^6$	-	-	-	-	38.9	51.3	87.2	38.1
	${10}^7$	-	-	-	-	391	912	337	312
	${10}^8$	-	-	-	-	3121	4003	3421	3115
1	${10}^3$	6.06	$23.65$	$14.04$	$6.04$	$0.22$	$0.15$	$0.31$	$0.14$
	${10}^4$	$\mathrm{𝟖𝟏𝟒}$	$\mathrm{𝟓𝟗𝟒𝟑}$	$\mathrm{𝟐𝟏𝟖𝟗}$	$\mathrm{𝟖𝟐𝟒}$	$0.52$	$0.59$	1.56	$0.47$
	${10}^5$	-	-	-	-	$27.4$	$29.8$	$10.7$	$9.2$
	${10}^6$	-	-	-	-	$37.4$	$48.8$	$83.1$	$37.3$
	${10}^7$	-	-	-	-	$\mathrm{𝟑𝟖𝟑}$	$\mathrm{𝟖𝟖𝟗}$	$\mathrm{𝟑𝟐𝟔}$	$\mathrm{𝟐𝟗𝟖}$
	${10}^8$	-	-	-	-	$\mathrm{𝟐𝟗𝟑𝟔}$	$\mathrm{𝟑𝟗𝟔𝟏}$	$\mathrm{𝟑𝟑𝟐𝟖}$	$\mathrm{𝟑𝟎𝟗𝟔}$
2	${10}^3$	6.12	28.64	16.02	6.16	$0.21$	0.18	0.32	0.14
	${10}^4$	874	6219	2415	868	0.59	0.62	1.54	0.51
	${10}^5$	-	-	-	-	30.0	31.9	12.1	9.2
	${10}^6$	-	-	-	-	40	53	86	38
	${10}^7$	-	-	-	-	417	1002	345	316
	${10}^8$	-	-	-	-	3148	4021	3510	3278

Table 3. TABLE III: Comparison of average CPU time (seconds) with different λ 𝜆 \lambda in the nonlinear case.

$\lambda$	$T$	GN-PRS	GN-SBM	GN-FOPD	GN-ADMM	IEKS-PRS	IEKS-SBM	IEKS-FOPD	IEKS-ADMM
0.01	${10}^2$	$0.26$	$1.15$	0.24	$0.38$	0.07	0.39	0.08	0.05
	${10}^3$	$\mathrm{𝟗𝟗}$	$\mathrm{𝟐𝟑𝟏}$	86	94	0.22	1.32	0.25	0.19
	${10}^4$	1504	5624	4951	1396	1.31	7.27	1.16	2.01
	${10}^5$	-	-	-	-	26.43	146.94	38.12	19.01
	${10}^6$	-	-	-	-	275	1687	389	201
	${10}^7$	-	-	-	-	2621	11542	3460	1963
0.1	${10}^2$	0.27	1.17	$0.23$	0.39	$0.06$	$0.37$	$0.07$	$0.05$
	${10}^3$	102	233	$\mathrm{𝟖𝟓}$	$\mathrm{𝟗𝟑}$	$0.22$	$1.31$	$0.24$	$0.17$
	${10}^4$	$\mathrm{𝟏𝟑𝟎𝟑}$	$5212.7$	$\mathrm{𝟒𝟔𝟒𝟎}$	$1252.6$	$1.29$	$7.26$	$1.16$	$1.99$
	${10}^5$	-	-	-	-	$26.21$	$146.81$	$37.45$	$18.72$
	${10}^6$	-	-	-	-	$\mathrm{𝟐𝟔𝟑}$	$\mathrm{𝟏𝟒𝟕𝟑}$	$\mathrm{𝟑𝟓𝟔}$	$\mathrm{𝟏𝟖𝟕}$
	${10}^7$	-	-	-	-	$\mathrm{𝟐𝟒𝟎𝟐}$	$\mathrm{𝟏𝟏𝟏𝟓𝟓}$	$\mathrm{𝟑𝟐𝟔𝟎}$	$\mathrm{𝟏𝟕𝟏𝟔}$
1	${10}^2$	0.32	1.28	0.26	0.43	0.06	0.39	0.07	0.05
	${10}^3$	105	241	89	102	0.23	1.38	0.31	0.26
	${10}^4$	1597	5711	5001	1450	1.32	7.31	1.19	2.03
	${10}^5$	-	-	-	-	27.21	152.1	39.45	20.12
	${10}^6$	-	-	-	-	291	1492	361	202
	${10}^7$	-	-	-	-	2645	11784	3519	2001
2	${10}^2$	0.34	1.28	0.27	0.43	0.06	0.43	0.07	0.07
	${10}^3$	111	249	95	121	0.26	1.41	0.34	0.27
	${10}^4$	1612	5821	5294	1510	1.48	7.48	1.23	2.45
	${10}^5$	-	-	-	-	28.45	167	41.24	22.12
	${10}^6$	-	-	-	-	312	1625	389	241
	${10}^7$	-	-	-	-	2741	12016	3645	2268

Table 4. TABLE IV: Average CPU time (seconds) for the outer iteration in the tomographic recostruction.

Dataset	$T$	$n_x$	$n_y$	ADMM	KF-ADMM
Emoji	$33$	16384	13020	3657.58	771.49
		16384	6510	1422.38	387.26
		4096	13020	284.83	97.93
		4096	6510	77.79	20.10
Lung	$10$	16384	13020	3584.62	764.78
		16384	6510	1378.41	367.13
		10404	13020	1475.84	487.25
		10404	6510	187.87	61.24

Equations

$$\begin{array}{r}\begin{array}{rl}{\mathbf{x}}_t& ={\mathbf{a}}_t({\mathbf{x}}_{t-1})+{\mathbf{q}}_t,\\ {\mathbf{y}}_t& ={\mathbf{h}}_t({\mathbf{x}}_t)+{\mathbf{r}}_t,\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{rl}& {\mathbf{x}}_{1:T}^{\star }={\mathrm{argmin}}_{{\mathbf{x}}_{1:T}}\frac{1}{2}\sum _{t=1}^T\Vert {\mathbf{y}}_t-{\mathbf{h}}_t({\mathbf{x}}_t){\Vert }_{{\mathbf{R}}_t^{-1}}^2+\lambda \sum _{t=1}^T{\Vert {\mathbf{\Omega }}_t{\mathbf{x}}_t\Vert }_1\\ & +\frac{1}{2}\sum _{t=2}^T\Vert {\mathbf{x}}_t-{\mathbf{a}}_t({\mathbf{x}}_{t-1}){\Vert }_{{\mathbf{Q}}_t^{-1}}^2+\frac{1}{2}\Vert {\mathbf{x}}_1-{\mathbf{m}}_1{\Vert }_{{\mathbf{P}}_1^{-1}}^2,\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{r}\underset{\mathbf{x}}{\min }{\theta }_1(\mathbf{x})+{\theta }_2(\mathbf{\Omega x}),\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{r}\underset{\mathbf{x},\mathbf{w}}{\min }{\theta }_1(\mathbf{x})+{\theta }_2(\mathbf{w}),s.t.\mathbf{w}=\mathbf{\Omega x}.\end{array}\end{array}$$

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

$$+\frac{\rho }{2}\Vert {\mathbf{w}}^{(k)}-\mathbf{\Omega x}{\Vert }^2,$$

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

$$+\frac{\rho }{2}\Vert \mathbf{w}-\mathbf{\Omega }{\mathbf{x}}^{(k+1)}{\Vert }^2,$$

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

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

$$+\frac{\rho }{2}\Vert {\mathbf{w}}^{(k)}-\mathbf{\Omega x}{\Vert }^2,$$

$${\mathbf{\eta }}^{(k+\frac{1}{2})}$$

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

$$+\frac{\rho }{2}\Vert \mathbf{w}-\mathbf{\Omega }{\mathbf{x}}^{(k+1)}{\Vert }^2,$$

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

$${\mathbf{x}}^{(k+1)}={\mathrm{argmin}}_{\mathbf{x}}{\theta }_1(\mathbf{x})+\frac{\rho }{2}\Vert {\mathbf{w}}^{(k)}-\mathbf{\Omega x}+{\mathbf{\eta }}^{(k)}{\Vert }_2^2,$$

$${\mathbf{w}}^{(k+1)}={\mathrm{argmin}}_{\mathbf{w}}{\theta }_2(\mathbf{w})+\frac{\rho }{2}\Vert \mathbf{w}-\mathbf{\Omega }{\mathbf{x}}^{(k+1)}+{\mathbf{\eta }}^{(k)}{\Vert }_2^2,$$

$$\begin{array}{r}\begin{array}{r}{\mathbf{\eta }}^{(k+1)}={\mathbf{\eta }}^{(k)}+({\mathbf{w}}^{(k+1)}-\mathbf{\Omega }{\mathbf{x}}^{(k+1)}).\end{array}\end{array}$$

$${\mathbf{w}}^{(k+1)}={\mathrm{argmin}}_{\mathbf{w}}{\theta }_2^{\ast }(\mathbf{w})+\frac{1}{2\gamma }\Vert \mathbf{w}-({\mathbf{w}}^{(k)}+\gamma \mathbf{\Omega }{\hat{\mathbf{x}}}^{(k)}){\Vert }^2,$$

$${\mathbf{x}}^{(k+1)}={\mathrm{argmin}}_{\mathbf{x}}{\theta }_1(\mathbf{x})+\frac{1}{2\rho }\Vert \mathbf{x}-({\mathbf{x}}^{(k)}-\rho {\mathbf{\Omega }}^{\top }{\mathbf{w}}^{(k+1)}){\Vert }^2,$$

$${\hat{\mathbf{x}}}^{(k+1)}={\mathbf{x}}^{(k+1)}+\tau ({\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}),$$

$$\begin{array}{r}\begin{array}{rl}& {\mathbf{x}}_{1:T}^{\star }={\mathrm{argmin}}_{{\mathbf{x}}_{1:T}}\frac{1}{2}\sum _{t=1}^T\Vert {\mathbf{y}}_t-{\mathbf{h}}_t({\mathbf{x}}_t){\Vert }_{{\mathbf{R}}_t^{-1}}^2\\ & +\frac{1}{2}\sum _{t=2}^T\Vert {\mathbf{x}}_t-{\mathbf{a}}_t({\mathbf{x}}_{t-1}){\Vert }_{{\mathbf{Q}}_t^{-1}}^2+\frac{1}{2}\Vert {\mathbf{x}}_1-{\mathbf{m}}_1{\Vert }_{{\mathbf{P}}_1^{-1}}^2.\end{array}\end{array}$$

$${\mathbf{a}}_t({\mathbf{x}}_{t-1})$$

$${\mathbf{h}}_t({\mathbf{x}}_t)$$

$$\begin{array}{r}\begin{array}{rl}& {\mathbf{x}}_{1:T}^{(i+1)}=\\ & {\mathrm{argmin}}_{{\mathbf{x}}_{1:T}}\frac{1}{2}\sum _{t=1}^T\Vert {\mathbf{y}}_t-{\mathbf{h}}_t({\mathbf{x}}_t^{(i)})-{\mathbf{J}}_{h_t}({\mathbf{x}}_t^{(i)})({\mathbf{x}}_t-{\mathbf{x}}_t^{(i)}){\Vert }_{{\mathbf{R}}_t^{-1}}^2\\ & +\frac{1}{2}\sum _{t=2}^T\Vert {\mathbf{x}}_t-{\mathbf{a}}_t({\mathbf{x}}_{t-1}^{(i)})-{\mathbf{J}}_{a_t}({\mathbf{x}}_{t-1}^{(i)})({\mathbf{x}}_t-{\mathbf{x}}_t^{(i)}){\Vert }_{{\mathbf{Q}}_t^{-1}}^2\\ & +\frac{1}{2}\Vert {\mathbf{x}}_1-{\mathbf{m}}_1{\Vert }_{{\mathbf{P}}_1^{-1}}^2.\end{array}\end{array}$$

$${\mathbf{m}}_t^{-}$$

$${\mathbf{P}}_t^{-}$$

$${\mathbf{S}}_t={\mathbf{J}}_{h_t}({\mathbf{x}}_t^{(i)}){\mathbf{P}}_t^{-}[{\mathbf{J}}_{h_t}({\mathbf{x}}_t^{(i)}){]}^{\top }+{\mathbf{R}}_t,$$

$${\mathbf{K}}_t={\mathbf{P}}_t^{-}[{\mathbf{J}}_{h_t}({\mathbf{x}}_t^{(i)}){]}^{\top }[{\mathbf{S}}_t{]}^{-1},$$

$$\displaystyle\mathbf{m}_{t}=\mathbf{m}_{t}^{-}+\mathbf{K}_{t}\Big{(}\mathbf{y}_{t}-\mathbf{h}_{t}(\mathbf{x}_{t}^{(i)})-\mathbf{J}_{h_{t}}(\mathbf{x}_{t}^{(i)})(\mathbf{m}_{t}^{-}-\mathbf{x}_{t}^{(i)})\Big{)},$$

$${\mathbf{P}}_t={\mathbf{P}}_t^{-}-{\mathbf{K}}_t{\mathbf{S}}_t[{\mathbf{K}}_t{]}^{\top },$$

$${\mathbf{G}}_t$$

$${\mathbf{m}}_t$$

$${\mathbf{P}}_t^s$$

$$\begin{array}{r}\begin{array}{rl}{\mathbf{a}}_t({\mathbf{x}}_{t-1})& ={\mathbf{A}}_t{\mathbf{x}}_{t-1},\\ {\mathbf{h}}_t({\mathbf{x}}_t)& ={\mathbf{H}}_t{\mathbf{x}}_t,\end{array}\end{array}$$

$$\mathbf{x}$$

$$\mathbf{y}$$

$$\mathbf{m}$$

$$\mathbf{H}$$

$$\mathbf{Q}$$

$$\mathbf{R}$$

$$\mathbf{\Omega }$$

$$\mathbf{\Psi }$$

Full text

I Introduction

STATE estimation problems naturally arise in many signal processing applications including target tracking, smart grids, and robotics [1, 2, 3]. In conventional Bayesian approaches, the estimation task is cast as a statistical inverse problem for restoring the original time series from imperfect measurements, based on a statistical model for the measurements given the signal together with a statistical model for the signal. In linear Gaussian models, this problem admits a closed-form solution, which can be efficiently implemented by the Kalman (or Rauch–Tung–Striebel) smoother (KS) [4, 2]. For nonlinear Gaussian models we can use various linearization and sigma-point-based methods [2] for approximate inference. In particular, here we use the so-called iterated extended Kalman smoother (IEKS) [5], which is based on analytical linearisation of the nonlinear functions. Although the aforementioned smoothers are often used to estimate dynamic signals, they lack a mechanism to promote sparsity in the signals.

One approach for promoting sparsity is to add an $L_{1}$-term to the cost function formulation of the state estimation problem. This approach imposes sparsity on the state estimate, which is either based on a synthesis sparse or an analysis sparse signal model. A synthesis sparse model assumes that the signal can be represented as a linear combination of basis vectors, where the coefficients are subject to, for example, an $L_{1}$-penalty, thus promoting sparsity. In the past decade, the use of synthesis sparsity for estimating dynamic signals has drawn a lot of attention [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. For example, a pseudo-measurement technique was used in the Kalman update equations for encouraging sparse solutions [7]. A method based on sparsity was applied compressive sensing to update Kalman innovations or filtering errors [8]. Based on synthesis sparsity, the estimation problem has been formulated as an $L_{1}$-regularized least square problem in [14]. Nevertheless, the previously mentioned methods only consider synthesis sparsity of the signal and assume a linear dynamic system.

On the other hand, analysis sparsity, also called cosparsity, assumes that the signal is not sparse itself, but rather the outcome is sparse or compressible in some transform domain, which leads to the flexibility in the modeling of signals [16, 17, 18, 19, 20]. Analysis sparse models involving an analysis operator – a popular choice being total variation (TV) – have been very successful in image processing. For example, several algorithms [18, 19] have been developed to train an analysis operator and the trained operators have been used for image denoising. In [21] the authors proposed to use the TV regularizer to improve the quality of image reconstruction. However, these approaches are not ideally suited for reconstructing dynamic signals. In state estimation problems, the available methods for analysis sparse priors are still limited. The main goal of this paper is to introduce these kinds of methods for dynamic state estimation.

Formulating a state estimation problem using synthesis and analysis sparsity leads to a general class of optimization problems, which require minimization of composite functions such as an analysis-$L_{1}$-regularized least-square problems. The difficulties arise from the appearance of the nonsmooth regularizer. There are various batch optimization methods such as proximal gradient method [22], Douglas-Rachford splitting (DRS) [23, 24], Peaceman-Rachford splitting (PRS) [25, 26], the split Bregman method (SBM) [27], the alternating method of multipliers (ADMM) [28, 29], and the first-order primal-dual (FOPD) method [30] for addressing this problem. However, these general methods do not take the inherent temporal nature of the optimization problem into account, which leads to bad computational and memory scaling in large-scale or high-dimensional data. This often renders the existing methods intractable due to their extensive memory and computational requirements.

As a consequence, we propose to combine a Kalman smoother with variable splitting optimization methods, which allows us to account for the temporal nature of the data in order to speed up the computations. In this paper, we derive novel methods for efficiently estimating dynamic signals with an extra (analysis) $L_{1}$-regularized term. The developed algorithms are based on using computationally efficient KS and IEKS for solving the subproblems arising within the steps of the optimization methods. Our experiments demonstrate promising performance of the methods in practical applications. The main contributions are as follows:

1. i)

We formulate the state estimation problem as an optimization problem that is based upon a general sparse model containing analysis or synthesis prior. The formulation accommodates a large class of popular sparsifying regularizers (e.g., synthesis $L_{1}$-norm, analysis $L_{1}$-norm, TV norm) for state estimation. 2. ii)

We present novel practical optimization methods, KS-ADMM and IEKS-ADMM, which are based on combining ADMM with KS and IEKS, respectively. 3. iii)

We also prove the convergence of the KS-ADMM method as well as the local convergence of the IEKS-ADMM method. 4. iv)

We generalize our smoother-based approaches to a general class of variable splitting techniques.

The advantage of the proposed approach is that the computational cost per iteration is much less than in the conventional batch solutions. Our approach is computationally superior to the state-of-the-art in a large-scale or high-dimensional state estimation applications.

The rest of the paper is organized as follows. We conclude this section by reviewing variable splitting methods and IEKS. Section II first develops the batch optimization by a classical ADMM method, and then presents a new KS-ADMM method for solving a linear dynamic estimation problem. Furthermore, for the nonlinear case, we present an IEKS-ADMM method in Section III and establish its local convergence properties. Section IV introduces a more general smoother-based variable splitting algorithmic framework. In particular, a general IEKS-based optimization method is formulated. Various experimental results in Section V demonstrate the effectiveness and accuracy in simulated linear and nonlinear state estimation problem. The performance of the algorithm is also illustrated in real-world tomographic reconstruction.

The notation of the paper is as follows. Vectors $\mathbf{x}$ and matrices $\mathbf{X}$ are indicated in boldface. $(\cdot)^{\top}$ stands for transposition, and $(\cdot)^{-1}$ is the matrix inversion. $\mathbf{x}_{1:T}$ stands for the time series from $\mathbf{x}_{1}$ to $\mathbf{x}_{T}$, and $\mathbf{x}^{(k)}$ denotes the value of $\mathbf{x}$ at $k$:th iteration. $\langle\mathbf{x},\mathbf{y}\rangle$ represents the vector inner product $\mathbf{x}^{\top}\mathbf{y}$. We denote by $\mathbb{R}^{n}$ the usual $n$ dimensional Euclidean space. The vector norm $\left\|\cdot\right\|_{p}$ for $p\geq 1$ is the standard $\ell_{p}$-norm. The $\mathbf{R}$-weighted Euclidean norm of a vector $\mathbf{x}$ is denoted by $\|\mathbf{x}\|_{\mathbf{R}}=\sqrt{\mathbf{x}^{\top}\mathbf{R}\mathbf{x}}$. $\theta^{*}$ is the conjugate of a convex function $\theta$, defined as $\theta^{*}(\mathbf{p})=\mathop{\sup}_{\mathbf{x}}\langle\mathbf{x},\mathbf{p}\rangle-\theta(\mathbf{x})$. $\operatorname{sgn}$ represents the signum function. $\operatorname{vec}(\cdot)$ is the vectorization operator, $\operatorname{blkdiag}(\cdot)$ is a block diagonal matrix operator with the elements in its argument on the diagonal, $\partial\phi(\mathbf{x})$ denotes a subgradient of $\phi$ at $\mathbf{x}$, $\mathbf{J}_{\phi}$ is the Jacobian of $\phi(\mathbf{x})$ and $\nabla\phi(\mathbf{x})$ and $\nabla^{2}\phi(\mathbf{x})$ are the gradient and Hessian of the function $\phi(\mathbf{x})$.

I-A Problem Formulation

Consider the dynamic state-space model [1, 2]

[TABLE]

where $t=1,\ldots,T$, $\mathbf{x}_{t}=\begin{bmatrix}x_{1,t}&x_{2,t}&\ldots&x_{n_{x},t}\end{bmatrix}^{\top}\in\mathbb{R}^{n_{x}}$ denotes an $n_{x}$-dimensional state of the system at the time step $t$, and $\mathbf{y}_{t}=\begin{bmatrix}y_{1,t}&y_{2,t}&\ldots&y_{n_{y},t}\end{bmatrix}^{\top}\in\mathbb{R}^{n_{y}}$ is an $n_{y}$-dimensional noisy measurement signal, $\mathbf{h}_{t}:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{y}}$ is a measurement function (typically with $n_{y}\leq n_{x}$), and $\mathbf{a}_{t}:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{x}}$ is a state transition function at time step $t$. The initial state $\mathbf{x}_{1}$ is assumed to have mean $\mathbf{m}_{1}$ and covariance $\mathbf{P}_{1}$. The errors $\mathbf{q}_{t}$ and $\mathbf{r}_{t}$ are assumed to be mutually independent random variables with known positive definite covariance matrices $\mathbf{Q}_{t}$ and $\mathbf{R}_{t}$, respectively. The goal is to estimate the state sequence $\mathbf{x}_{1:{T}}=\left\{\mathbf{x}_{1},\ldots,\mathbf{x}_{T}\right\}$ from the noisy measurement sequence $\mathbf{y}_{1:{T}}=\{\mathbf{y}_{1},\ldots,\mathbf{y}_{T}\}$. In this paper, we focus on estimating $\mathbf{x}_{1:{T}}$ by minimizing the sum of quadratic error terms and an $L_{1}$ sparsity-promoting penalty.

For the sparsity assumption, we add an extra $L_{1}$-penalty for the state $\mathbf{x}_{t}$, and then formulate the optimization problem as

[TABLE]

where $\mathbf{x}^{\star}_{1:{T}}$ is the optimal state sequence, $\mathbf{\Omega}_{t}$ is a linear operator, and $\lambda$ is a penalty parameter, which describes a trade-off between the data fidelity term and the regularizing penalty term. The formulation (2) encompasses two particular cases: by setting $\mathbf{\Omega}_{t}$ to a diagonal matrix (e.g., identity matrix $\mathbf{\Omega}_{t}=\mathbf{I}$), a synthesis sparse model is obtained, which assumes that $\mathbf{x}_{1:T}$ are sparse. Such a case arises frequently in state estimation applications [15, 10, 11, 31, 12]. Correspondingly, an analysis sparse model is obtained when a more general $\mathbf{\Omega}_{t}$ is used. For example, the TV regularization, which is common in tomographic reconstruction, can be obtained by using a finite-difference matrix as $\mathbf{\Omega}_{t}$.

More generally, $\mathbf{\Omega}_{t}$ can be a fixed matrix [32, 20, 16] or a learned matrix [17, 18, 19]. It should be noted that, if the $L_{1}$ term is not used (i.e., when $\lambda=0$) in (2), the objective can be solved by using a linear or non-linear KS [4, 5, 2]. However, when $\lambda>0$, the smoothing is no longer applicable, and the cost function is non-differentiable.

Since $\|\mathbf{\Omega}_{t}\mathbf{x}_{t}\|_{1}$ does not have a closed-form proximal operator in general, we employ variable splitting technique for solving the resulting optimization problem. As mentioned above, many variable splitting methods can be used to solve (2), such as PRS [26], SBM [27], ADMM [28], DRS [23], and FOPD [30]. Especially, ADMM is a popular member of this class. Therefore, we start by presenting algorithms based on ADMM and then extend them to more general variable splitting methods. In the following, we review variable splitting and IEKS methods, before presenting our approach in detail.

I-B Variable Splitting

The methods we develop in this paper are based on variable splitting [33, 34]. Consider an unconstrained optimization problem in which the objective function is the sum of two functions

[TABLE]

where $\theta_{2}(\cdot)=\|\cdot\|_{1}$, and $\mathbf{\Omega}$ is a matrix. Variable splitting refers to the process of introducing an auxiliary constrained variable $\mathbf{w}$ to separate the components in the cost function. More specifically, we impose the constraint $\mathbf{w}=\mathbf{\Omega}\mathbf{x}$, which transforms the original minimization problem (3) into an equivalent constrained minimization problem, given by

[TABLE]

The minimization problem (4) can be solved efficiently by classical constrained optimization methods [35]. The rationale of variable splitting is that it may be easier to solve the constrained problem (4) than the unconstrained one (3). PRS, SBM, FOPD, ADMM, and their variants [36] are a few well-known variable splitting methods – see also [37, 38] for a recent historical overview.

ADMM [28] is one of the most popular algorithms for solving (4). ADMM defines an augmented Lagrangian function, and then alternates between the updates of the split variables. Given $\mathbf{x}^{(0)}$, $\mathbf{w}^{(0)}$, and $\bm{\eta}^{(0)}$, its iterative steps are:

[TABLE]

where $\bm{\eta}$ is a Lagrange multiplier and $\rho$ is a parameter.

The PRS method [25, 26] is similar to ADMM except that it updates the Lagrange multiplier twice. The typical iterative steps for (3) are

[TABLE]

where $\alpha\in(0,1)$.

In SBM [27], we iterate the steps

[TABLE]

$M$ times, and update the extra variable by

[TABLE]

When $M=1$, this is equivalent to ADMM.

There are also other variable splitting methods which alternate proximal steps for the primal and dual variables. One example is FOPD [30], where the $(k+1)$:th iteration consists of the following

[TABLE]

where $\tau$ and $\gamma$ are parameters.

All these variable splitting algorithms provide simple ways to construct efficient iterative algorithms that offer simpler inner subproblems. However, the subproblems such as (5), (6), (7a) and (9b) remain computationally expensive, as they involve large matrix-vector products when the dimensionality of $\mathbf{x}$ is large. We circumvent this problem by combining variable splitting with KS and IEKS.

I-C The Iterated Extended Kalman Smoother

IEKS [5] is an approximative algorithm for solving non-linear optimal smoothing problems. However, it can also be seen as an efficient implementation of the Gauss–Newton algorithm for solving the problem

[TABLE]

That is, it produces the maximum a posteriori (MAP) estimate of the trajectory. The IEKS method works by alternating between linearisation of $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ around a previous estimate $\mathbf{x}^{(i)}_{1:{T}}$, as follows:

[TABLE]

and solving the linearized problem

[TABLE]

The solution of (12) can in turn be efficiently obtained by the Rauch–Tung–Striebel (RTS) smoother [4], which first computes the filtering mean and covariances $\mathbf{m}_{1:T}$ and $\mathbf{P}_{1:T}$, by alternating between prediction

[TABLE]

and update

[TABLE]

where $\mathbf{S}_{t}$ and $\mathbf{K}_{t}$ are the innovation covariance matrix and the Kalman gain at the time step $t$, respectively. The filtering means $\mathbf{m}_{t}$ and covariances $\mathbf{P}_{t}$ are then corrected in a backwards (smoothing) pass

[TABLE]

Now setting $\mathbf{x}^{(i+1)}_{t}=\mathbf{m}_{t}^{s}$ gives the solution to (12). When the functions $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ are linear, the above iteration converges in a single step. This algorithm is the classical RTS smoother or more briefly KS [4].

In this paper, we use the KS and IEKS algorithms as efficient methods for solving generalized versions of the optimization problems given in (10), which arise within the steps of variable splitting.

II Linear State Estimation by KS-ADMM

In this section, we present the KS-ADMM algorithm which is a novel algorithm for solving $L_{1}$-regularized linear Gaussian state estimation problems. In particular, Section II-A describes the batch solution by ADMM. Then, by defining an artificial measurement noise and a pseudo-measurement, we formulate the KS-ADMM algorithm to solve the primal variable update in Section II-B.

II-A Batch Optimization

Let us assume that the state transition function $\mathbf{a}_{t}$ and the measurement function $\mathbf{h}_{t}$ are linear, denoted by

[TABLE]

where $\mathbf{A}_{t}$ and $\mathbf{H}_{t}$ are the transition matrix and the measurement matrix, respectively. In order to reduce this problem to (3), we stack the entire state sequence into a vector, which transforms the objective into a batch optimization problem. Thus, we define the following variables

[TABLE]

The optimization problem introduced in Section I-A can now be reformulated as the following batch optimization problem

[TABLE]

which in turn can be seen to be a special case of (3). Here, our algorithm for solving (18) builds upon the batch ADMM [28].

To derive an ADMM algorithm for (18), we introduce an auxiliary variable $\mathbf{w}=\operatorname{vec}(\mathbf{w}_{1},\ldots,\mathbf{w}_{T})$ and a linear equality constraint $\mathbf{w}=\mathbf{\Omega}\mathbf{x}$. The resulting equality-constrained problem is formulated mathematically as

[TABLE]

The main objective here is to find a stationary point $(\mathbf{x}^{\star},\mathbf{w}^{\star},\bm{\eta}^{\star})$ of the augmented Lagrangian function associated with (19) as the function

[TABLE]

where $\bm{\eta}\in\mathbb{R}^{TP}$ is the dual variable and $\rho$ is a penalty parameter. As described in Section I-B, at each iteration of ADMM we perform the updates

[TABLE]

The update for the primal sequence $\mathbf{x}$ is equivalent to the quadratic optimization problem given by

[TABLE]

which has the closed-form solution

[TABLE]

For the dual sequence $\mathbf{w}$, the iteration in (21b) can be solved by [39]

[TABLE]

where $\mathbf{e}^{(k)}=\mathbf{\Omega}\mathbf{x}^{(k+1)}+\bm{\eta}^{(k)}/\rho$.

While the optimization problem (22) can be solved in closed-form, direct solution is computationally demanding, especially when the number of time points or the dimensionality of the state is large. However, the problem can be recognized to be a special case of optimization problems where the iterations can be solved by KS (see Section I-C) provided that we add pseudo-measurements to the problem. In the following, we present the resulting algorithm.

II-B The KS-ADMM Solver

The proposed KS-ADMM solver is described in Algorithm 1. To extend the batch ADMM to KS-ADMM, we first define an artificial measurement noise $\mathbf{\Sigma}_{t}=\mathbf{I}/\rho$ and a pseudo-measurement $\mathbf{z}_{t}=\mathbf{w}_{t}+\bm{\eta}_{t}/\rho$, and then rewrite (22) as

[TABLE]

The solution to (25) can then be computed by running KS on the state estimation problem

[TABLE]

Here, $\bm{\sigma}_{t}$ is an independent random variable with covariance $\bm{\Sigma}_{t}$. The KS-based solution can be described as a four stage recursive process: prediction, $\mathbf{y}_{t}$-update, $\mathbf{z}_{t}$-update, and a RTS smoother which should be performed for $t=1,\ldots,T$. First, the prediction step is given by

[TABLE]

where $\mathbf{m}_{t}^{-}$ and $\mathbf{P}_{t}^{-}$ are the predicted mean and covariance at the time $t$. Secondly, the update steps for $\mathbf{y}_{t}$ are given by

[TABLE]

Thirdly, the update steps for $\mathbf{z}_{t}$ are

[TABLE]

Here, $\mathbf{S}^{y}_{t}$ and $\mathbf{S}^{z}_{t}$, $\mathbf{K}^{y}_{t}$ and $\mathbf{K}^{z}_{t}$, $\mathbf{m}^{y}_{t}$ and $\mathbf{m}_{t}$, $\mathbf{P}^{y}_{t}$ and $\mathbf{P}_{t}$ are the innovation covariances, gain matrices, means, and covariances for the variables $\mathbf{y}_{t}$ and $\mathbf{z}_{t}$ at the time step $t$, respectively. Finally, we run a RTS smoother [4] for $t$ = $T-1,\ldots,1$, which has the steps

[TABLE]

where $\mathbf{m}_{T}^{s}=\mathbf{m}_{T}$ and $\mathbf{P}_{T}^{s}=\mathbf{P}_{T}$ (see [2] for more details). This gives the update for $\mathbf{x}_{1:T}$ as:

[TABLE]

The remaining updates for $t=1,\ldots,T$ are given as

[TABLE]

where $\mathbf{e}_{t}^{(k)}=\mathbf{\Omega}_{t}\mathbf{x}_{t}^{(k+1)}+\bm{\eta}_{t}^{(k)}/\rho$, and

[TABLE]

It is useful to note that in Algorithm 1, the covariances and gains are independent of the iteration number and thus can be pre-computed outside the ADMM iterations. Furthermore, when the model is time-independent, we can often use stationary Kalman filters and smoothers instead of their general counterparts which can further be used to speed up the computations.

II-C Convergence of KS-ADMM

In this section, we discuss the convergence of KS-ADMM. If the system (26) is detectable [40], then the objective function (20) is convex. The traditional convergence results for ADMM such as in [28, 41] then ensure that the objective globally converges to the stationary (optimal) point $(\mathbf{x}_{1:T}^{\star},\mathbf{w}_{1:T}^{\star},\bm{\eta}_{1:T}^{\star})$. The result is given in the following.

Theorem 1** (Convergence of KS-ADMM).**

Assume that the system (26) is detectable [40]. Then, for a constant $\rho$, the sequence $\{\mathbf{x}_{1:T}^{(k)},\mathbf{w}_{1:T}^{(k)},\bm{\eta}_{1:T}^{(k)}\}$ generated by Algorithm 1 from any starting point $\{\mathbf{x}_{1:T}^{(0)},\mathbf{w}_{1:T}^{(0)},\bm{\eta}_{1:T}^{(0)}\}$ converges to the stationary point $\{\mathbf{x}_{1:T}^{\star},\mathbf{w}_{1:T}^{\star},\bm{\eta}_{1:T}^{\star}\}$ of (20).

Proof.

Due to the detectability assumption, the objective function is convex, and thus the result follows from the classical ADMM convergence proof [28, 41].

III Nonlinear State Estimation by IEKS-ADMM

When $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ are nonlinear, the $\mathbf{x}$ subproblem arising in the ADMM iteration cannot be solved in closed form. In the following, we first present a batch solution of the nonlinear case based on a Gauss–Newton (GN) iteration and then show how it can be efficiently implemented by IEKS.

III-A Batch Optimization

Let us now consider the case where the state transition function $\mathbf{a}_{t}$ and the measurement function $\mathbf{h}_{t}$ in (1) are nonlinear. We now proceed to rewrite the optimization (2) in batch form by defining the following variables

[TABLE]

Note that the variables $\mathbf{x}$, $\mathbf{y}$, $\mathbf{m}$, $\mathbf{Q}$, $\mathbf{R}$ and $\mathbf{\Omega}$ have the same definitions as (17). Using these variables, the $\mathbf{x}$ subproblem can be naturally transformed into

[TABLE]

which is also a special case of (3), similarly to the linear case.

Following the ADMM, we define the augmented Lagrangian function associated with (35) as:

[TABLE]

Since the nonlinear batch solution is based on ADMM, the iteration steps of $\mathbf{w}$ and $\bm{\eta}$ are the same with the linear case (see Equations (LABEL:eq:w_compute_2) and (21c)). Here, we focus on introducing the solution of the primal variable $\mathbf{x}$.

When updating $\mathbf{x}$, the objective is no longer a quadratic function. However, the optimization problem can be solved with GN [42]. Here, the $\mathbf{x}$ subproblem is rewritten as

[TABLE]

where

[TABLE]

Then, the gradient of $f(\mathbf{x})$ is given by

[TABLE]

where

[TABLE]

and the Hessian is $\nabla^{2}f(\mathbf{x})=\mathbf{J}^{\top}\mathbf{J}(\mathbf{x})+\mathbf{H}(\mathbf{x})$, where

[TABLE]

In GN, avoiding the trouble of computing the residual $[\mathbf{H}(\mathbf{x})]_{ij}$, we use the approximation $\nabla^{2}f(\mathbf{x})\approx\mathbf{J}^{\top}\mathbf{J}(\mathbf{x})$ to replace $\nabla^{2}f(\mathbf{x})$, which means $[\mathbf{H}(\mathbf{x})]_{ij}$ is assumed to be small enough. Thus, the primal variable in $\mathbf{x}$ iteration is updated by:

[TABLE]

The iterations can stop after a maximum number of iterations $i_{\max}$ or if the condition $\|\mathbf{x}^{(i+1)}-\mathbf{x}^{(i)}\|_{2}\leq\varepsilon$ is satisfied, where $\varepsilon$ is an error tolerance. If $\varepsilon$ is small enough, then it means that the above algorithm has (almost) converged. The rest of the ADMM updates can be implemented similarly to the linear Gaussian case.

III-B The IEKS-ADMM Solver

We now move on to consider the IEKS-ADMM solver. As discussed in Section I-C, IEKS can be seen as an efficient implementation of the GN method, which inspires us to derive an efficient implementation of the batch ADMM.

Now, we rewrite the $\mathbf{x}$ subproblem (37) as

[TABLE]

In a modest scale (e.g., $T\approx 10^{3}$), $\mathbf{x}_{1:T}$ can be directly computed by (40) although its computations scale as $\mathcal{O}(n_{x}^{3}\times T^{3})$. When $T$ is large, the batch ADMM will have high memory and computational requirements. In this case, the use of IEKS becomes beneficial due to its linear computational scaling. In this paper, the proposed method incorporates IEKS into ADMM to design the IEKS-ADMM algorithm for solving the nonlinear case.

In the IEKS algorithm, the Gaussian smoother is run several times with $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ and their Jacobians are evaluated at the previous (inner loop) iteration. The detailed iteration steps of IEKS-ADMM are described in Algorithm 2. In particular, following the prediction steps (13) in Section I-C, the update steps for $\mathbf{y}_{t}$ are given by

[TABLE]

and for the pseudo-measurement $\mathbf{z}_{t}$, the update steps are the same as in the linear case. They are given in (29).

Additionally, the RTS smoother steps are also described in Section I-C. We can then obtain the solution as $\mathbf{x}_{1:T}^{(k+1)}=\mathbf{m}^{s}_{1:{T}}$. Note that the updates on $\mathbf{w}_{1:T}$ and $\bm{\eta}_{1:T}$ can be implemented by (LABEL:eq:wt_compute) and (33), respectively.

III-C Convergence of IEKS-ADMM

In this section, our aim is to prove the convergence of the IEKS-ADMM algorithm. Although we can rely much on existing convergence results, unfortunately, when $\mathbf{a}(\mathbf{x})$ and $\mathbf{h}(\mathbf{x})$ are nonlinear, the traditional convergence analysis [28, 43, 44, 45] does not work as such. In particular, Jacobian matrices $\mathbf{J}_{a}$, $\mathbf{J}_{h}$ and linear operator $\mathbf{\Omega}$ in this paper are possibly rank-deficient, which is not covered by the existing convergence results. In the following, we will establish the convergence analysis which also covers this case.

For notational convenience, we define $\theta_{1}(\mathbf{x})=\frac{1}{2}\|\mathbf{y}-\mathbf{h}(\mathbf{x})\|_{\mathbf{R}^{-1}}^{2}+\frac{1}{2}\left\|\mathbf{m}-\mathbf{a}(\mathbf{x})\right\|_{\mathbf{Q}^{-1}}^{2}$ and $\theta_{2}(\mathbf{w})=\lambda\|\mathbf{w}\|_{1}$. The variables $\mathbf{x}$ and $\mathbf{w}$ are two sets of time series, and $\theta_{1}(\mathbf{x})$ is a non-quadratic, possibly nonconvex function. The corresponding augmented Lagrangian function can be rewritten as

[TABLE]

where $\mathbf{\Omega}$ can be full-row rank or full-column rank. Thus, the convergence is analyzed in two different cases. We make the following assumptions.

Assumption 1**.**

The gradient $\nabla\theta_{1}(\mathbf{x})$ is Lipschitz continuous with constant $L_{\theta_{1}}$, that is,

[TABLE]

Assumption 2**.**

Function $\theta_{1}(\mathbf{x})+\theta_{2}(\mathbf{w})$ is lower bounded and coercive over the feasible set $\{(\mathbf{x},\mathbf{w}):\mathbf{w}=\mathbf{\Omega}\mathbf{x}\}$.

First, we prove that the sequence $\mathcal{L}(\mathbf{x}^{(k)},\mathbf{w}^{(k)};\bm{\eta}^{(k)})$ is monotonically non-increasing in the following lemma.

Lemma 1** (Nonincreasing sequence).**

Let Assumptions 1 and 2 be satisfied and $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ be the iterative sequence generated by ADMM. Assume that one of two cases is satisfied:

Case (a): There exists $\rho_{0}$ such that when $\rho>\rho_{0}$, $\mathbf{x}\mapsto\mathcal{L}(\mathbf{x},\mathbf{w};\bm{\eta})$ is $\mu_{x}$-strongly convex, that is, $\mathcal{L}(\mathbf{x},\mathbf{w};\bm{\eta})$ satisfies

[TABLE]

Furthermore, assume that $\rho>\max\left(\frac{2L_{\theta_{1}}^{2}}{\kappa_{a}^{2}\mu_{x}},\rho_{0}\right)$ and that $\mathbf{\Omega}$ has full row rank with

[TABLE]

Case (b): $\rho>\frac{L_{\theta_{1}}}{\kappa_{b}^{2}}$, and $\mathbf{\Omega}$ has full-column rank with

[TABLE]

Then, sequence $\mathcal{L}(\mathbf{x}^{(k)},\mathbf{w}^{(k)};\bm{\eta}^{(k)})$ is nonincreasing in $k$.

Proof.

See Appendix A.

Next we prove the convergence of Algorithm 2. For that we need a couple of lemmas which are presented in the following.

Lemma 2** (Convergence of GN).**

Let $\nabla f(\mathbf{x})$ be Lipschitz continuous with constant $L_{f}$, $\mathbf{H}(\mathbf{x})$ be bounded by a constant, that is, $\|\mathbf{H}(\mathbf{x})\|\leq\epsilon_{h}$. If $\mathbf{J}^{\top}\mathbf{J}(\mathbf{x}^{(i)})\succeq\mu^{2}\mathbf{I}$ where $\mu>0$ is a constant, and $\epsilon_{h}<\mu^{2}$, then the sequence $\mathbf{x}^{(i)}$ converges to a local minimum $\mathbf{x}^{\star}$. In particular, the convergence is quadratic when $\epsilon_{h}\to 0$, and (at least) linear convergence is obtained when $\epsilon_{h}<\mu^{2}$.

Proof.

See Appendix B.

Based on Lemmas 1 and 2, we can now establish the convergence by GN-ADMM in the following lemma.

Lemma 3** (Convergence of GN-ADMM).**

Let assumptions of Lemmas 1 and 2 be satisfied. Then, the sequence $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ generated by GN-ADMM algorithm converges to a local minimum $(\mathbf{x}^{\star},\mathbf{w}^{\star},\bm{\eta}^{\star}).$

Proof.

By Lemma 1, the sequence $\mathcal{L}(\mathbf{x}^{(k)},\mathbf{w}^{(k)};\bm{\eta}^{(k)})$ is nonincreasing in $k$. By Assumption 2, the sequence $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ is bounded, because $\mathcal{L}(\mathbf{x}^{(k)},\mathbf{w}^{(k)};\bm{\eta}^{(k)})$ is upper bounded by $\mathcal{L}(\mathbf{x}^{(0)},\mathbf{w}^{(0)};\bm{\eta}^{(0)})$ and nonincreasing. It is also lower bounded by

[TABLE]

By Lemma 2, there exists a local minimum $\mathbf{x}^{\star}$ such that the sequence ${\mathbf{x}^{(i)}}$ converges to $\mathbf{x}^{\star}$, which is a local minimum of $\mathbf{x}$ subproblem. The $\mathbf{w}$ subproblem is convex [46] and thus there exists a unique minimum $\mathbf{w}^{\star}$. We then deduce that the iterative sequence $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ generated by GN-ADMM converges to $(\mathbf{x}^{\star},\mathbf{w}^{\star},\bm{\eta}^{\star})$.

The equivalence of GN and IEKS can now be used to show that IEKS-ADMM converges to a local minimum $(\mathbf{x}^{\star}_{1:T},\mathbf{w}^{\star}_{1:T},\bm{\eta}^{\star}_{1:T})$.

Theorem 2** (Convergence of IEKS-ADMM).**

If the sequence $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ generated by GN-ADMM algorithm converges to a local minimum $(\mathbf{x}^{\star},\mathbf{w}^{\star},\bm{\eta}^{\star})$, then the sequence $\{\mathbf{x}^{(k)}_{1:T},\mathbf{w}^{(k)}_{1:T},\bm{\eta}^{(k)}_{1:T}\}$ generated by IEKS-ADMM algorithm converges to the local minimum $(\mathbf{x}^{\star}_{1:T},\mathbf{w}^{\star}_{1:T},\bm{\eta}^{\star}_{1:T}).$

Proof.

According to[47], the sequence $\{\mathbf{x}^{(k)},\mathbf{w}^{(k)},\bm{\eta}^{(k)}\}$ generated by the GN method and the sequence $\{\mathbf{x}^{(k)}_{1:T},\mathbf{w}^{(k)}_{1:T},\bm{\eta}^{(k)}_{1:T}\}$ generated by IEKS are identical. Based on Lemma 3, we deduce the iterative sequence $\{\mathbf{x}^{(k)}_{1:T},\mathbf{w}^{(k)}_{1:T},\bm{\eta}^{(k)}_{1:T}\}$ generated by IEKS-ADMM is locally convergent to $(\mathbf{x}^{\star}_{1:T},\mathbf{w}^{\star}_{1:T},\bm{\eta}^{\star}_{1:T})$.

IV Extension to general algorithmic Framework

IV-A The Proposed Framework

In this subsection, we present a general algorithmic framework based on the combination of the extended Kalman smoother and variable splitting. As the smoother solution only applies to the $\mathbf{x}_{1:T}$-subproblem, here we only formulate the corresponding $\mathbf{x}_{1:T}$-subproblem which can be solved with IEKS. The different variants in the proposed framework are distinguished by three different choices: the pseudo-measurement, the pseudo-measurement covariance, and the pseudo-measurement model matrix.

When $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ are linear functions, we have the following general objective function for the $\mathbf{x}_{1:T}$-subproblem:

[TABLE]

and when $\mathbf{a}_{t}$ and $\mathbf{h}_{t}$ are nonlinear functions, we have

[TABLE]

In the above objective functions, $\mathbf{\Delta}_{t}$ is the pseudo-measurement, $\mathbf{\Sigma}_{t}$ is the pseudo-measurement covariance, and $\mathbf{\Theta}_{t}$ is the pseudo-measurement model matrix.

As mentioned in Section I-B, various variable splitting such as PRS, SBM, and FOPD can be used to solve the problems (22) and (37). Their KS / IEKS-based counterparts can be obtained by selecting the aforementioned pseudo-measurement model parameters as shown in Table I. The algorithms for solving the optimization problems are then the same as Algorithms 1 and 2 except that the pseudo-measurement updates in (29) are replaced with

[TABLE]

and the updates of the other variables are performed using the appropriate algorithm (see Section I-B).

IV-B Computational Complexity

This section investigates the computational complexity of the KS / IEKS based variable splitting methods. The proposed methods are iterative, in that we use several numbers of iterations to compute the minimal points also for the primal variable $\mathbf{x}^{\star}_{1:T}$ in (35). However, we can always use a bounded number of iterations and thus we only need to determine the complexity of a single iteration to determine the complexity of the whole method. In our case, the computational burden of the auxiliary variable and the dual variable is low compared with the matrix inversions in the primal variable update. Asymptotically, the computational complexities of (LABEL:eq:wt_compute) and (33) are both $\mathcal{O}(n_{x}^{2}\,T)$.

In our method, we compute the primal variable update using KS and IEKS, instead of computing matrix inversions explicitly. The time complexity of (iteration of) KS and IEKS is $\mathcal{O}(n_{x}^{3}\,T)$ [5, 2, 48] (assuming $n_{y}\leq n_{x}$), while the dominating computation in batch variable splitting methods is the matrix inversion with $\mathcal{O}(n_{x}^{3}\,T^{3})$ complexity [28]. Because of the total $\mathcal{O}(n_{x}^{3}\,T)$ computational complexity, the proposed method is especially applicable to large-scale dataset.

V Numerical Experiments

In the following, we demonstrate the KS and IEKS based variable splitting methods in numerical experiments. We first provide several simulated results to study the performance with varying regularization parameter. Then, we turn our attention to the behavior of the proposed methods with respect to the convergence curve and the computational efficiency. Finally, we report the results for large-scale signal estimation and demonstrate the effectiveness of the methodology in a tomographic reconstruction task.

V-A Linear Gaussian Simulation Experiment

Consider a four-dimensional linear tracking model (see, e.g., [2]) where the state contains rectangular coordinates $x_{1}$ and $x_{2}$, and velocity variables $x_{3}$ and $x_{4}$. The state of the system at time step $t$ is $\mathbf{x}_{t}=\begin{bmatrix}x_{1,t}&x_{2,t}&x_{3,t}&x_{4,t}\end{bmatrix}^{\top}$. The transition and measurement model matrices are

[TABLE]

The matrix $\mathbf{\Omega}_{t}$ and the covariance for the transition are

[TABLE]

with $q_{c}=1/2$, $\Delta t=0.1$, the measurement noise covariance $\mathbf{R}_{t}=\sigma^{2}\mathbf{I}$ with $\sigma=0.2$, and $T=100$ (small scale). The relative error is calculated by

[TABLE]

where $\mathbf{x}_{t}^{\text{true}}$ is the ground truth at time step $t$ and $\mathbf{x}_{t}^{(k)}$ is the $k$:th iterate at time step $t$. The goal here is to estimate dynamic signals from the noisy measurements $\mathbf{y}_{1:T}$.

In this experiment, we first illustrate the computations for $1000$ values of the regularizing penalty parameter $\lambda$ in the interval $[0.01,\,10]$. We remark that other parameters, for example, the parameter $\rho$ in ADMM and KS-ADMM, are chosen according to the existing guidelines [28], with no aim at further optimizing the convergence performance. The CPU times of various KS-based variable splitting methods are listed in Fig. 1. The plotted result is an average over $30$ experiments. For $1000$ values of parameter $\lambda$, the total number of ADMM iterations required is less than $20$, which takes around $0.02$ seconds in total. Thus, in the small-scale dataset, the parameter $\lambda$ has a less effect on the computational complexity when $\lambda$ is varying. Fig. 2 shows the relative error as a function of regularization parameter $\lambda$, using KS-PRS, KS-SBM, KS-FOPD, and KS-ADMM. As expected, we observe that the relative error is dependent on a proper choice of $\lambda$. We test different methods for the parameter $\lambda$, and find empirically that the lowest relative errors are achieved with $\lambda=1$.

Additionally, we compare the convergence speed (relative error versus iteration number) and the running time (average CPU time versus iteration number) generated by batch versions of PRS [26], SBM [27], FOPD [30], ADMM [28], to the proposed KS-based variable splitting methods. We also evaluate the performance without adding an extra analysis-$L_{1}$-regularization term (i.e., $\lambda=0$) in which case the optimization problem (25) can be solved by KS. Fig. 3 (a) shows the number of iterations required to solve the estimation problem. In tens of iteration numbers, all the methods have roughly the same relative errors. Fig. 3 (b) shows the average CPU time as function of number of iterations. Note that in order to speed up the KS-based variable splitting methods, we compute the gains $\mathbf{K}_{t}^{y}$, $\mathbf{K}_{t}^{z}$ and $\mathbf{G}_{t}$ only at the first iteration, and use the pre-computated matrices in the following iterations. Although PRS and KS-PRS, SBM and KS-SBM, FOPD and KS-FOPD, ADMM and KS-ADMM have the same convergence speed, not surprisingly, KS-SBM, KS-FOPD and KS-ADMM have a lower CPU time. When $T=100$, KS and the KS-based variable splitting methods have similar CPU time, but KS has a worse relative error. As shown in Fig. 3 (b), running the PRS, SBM, FOPD and ADMM solvers is time-consuming. The KF-based variable splitting methods take about 0.03 seconds to reach 10 iterations, while the classical optimization approaches such as FOPD and ADMM take $7$ times longer.

The benefit of our approach is highlighted by the fact that the methods can efficiently solve a large-scale dynamic signal estimation problem with an extra $L_{1}$-regularized term. Next, we enlarge the time step count from $10^{3}$ to $10^{8}$. All the results reported in Fig. 4 are obtained with $\lambda=1$, which gives the smallest relative error for all the methods. We use $10$ iterations for each method, which in practice is enough for convergence. It can be seen that the proposed method significantly outperforms other batch solutions with respect to the consumed CPU time. In particular, the PRS, SBM, FOPD and ADMM solvers with $10^{4}$ time steps take more time than the proposed methods with $T=10^{8}$. When $T=10^{5},10^{6},10^{7},10^{8}$, the computing operations on PRS, SBM, FOPD, and ADMM run out of memory, and the related results cannot be reported. This is mainly because the KS-based variable splitting methods deal with the objective using recursive computations, which significantly reduces the computational and memory burden, while the batch optimization methods explicitly deal with large vectors and matrices.

Table II summarizes the average CPU time with different regularization parameters $\lambda$ when the number of time steps $T$ is varying from $10^{2}$ to $10^{8}$. The table reports the average CPU time in seconds required by each solver with $10$ iterations. Not surprisingly, the KF-based variable splitting methods (i.e., KF-PRS, KF-SBM, KF-FOPD, and KF-ADMM) are much faster than the batch variable splitting methods when $T$ grows.

V-B Nonlinear Simulation Experiment

We consider a five-dimensional nonlinear coordinated turn model [1]. We set the measurement noise covariance $\mathbf{R}_{t}=\sigma^{2}\mathbf{I}$ with $\sigma=0.3$, $q_{c}=0.01$, $\Delta t=0.2$, and run $k_{\max}=20$ iterations of all the optimization methods. Before moving on, we provide some empirical evidence to support the choice of the regularization parameter $\lambda$ in the IEKS-based variable splitting methods. Similarly to the linear case in Section V-A, we plot the relative errors obtained by varying $\lambda$ in Fig. 5. It can be seen that the IEKS-based variable splitting methods have similar relative errors with $\lambda$ varying in the range $[0.1,\,1]$. Next, we select $\lambda=0.1$ for the following experiments.

Then, we compare IEKS[5], GN-PRS and IEKS-PRS, GN-SBM and IEKS-SBM, GN-FOPD, and IEKS-FOPD, GN-ADMM, and IEKS-ADMM by plotting the relative error and the CPU time as functions of the iteration number. Fig. 6 demonstrates the efficiency of our IEKS-based variable splitting methods against GN-PRS, GN-SBM, GN-FOPD and GN-ADMM in the same experiment. The horizontal axis in Fig. 6 (a) describes the total iteration number, and the vertical axis gives the relative errors. As can be seen, all the methods give the fast convergence in around $5$ iterations.

We are also interested in the performance without adding an extra analysis-$L_{1}$-regularized term (i.e., $\lambda=0$). Since there is no outer iteration in IEKS, we plot the relative error of IEKS after the inner iteration (dashed black line in Fig. 6 (b)). In contrast with the estimation results, we observe a performance gap between the variable splitting methods and IEKS. This gap reveals the benefit of the extra regularization term that is used in these methods. In the average CPU time, IEKS-PRS, IEKS-SBM, IEKS-FOPD, and IEKS-ADMM are clearly superior to batch variable splitting (see Fig. 6 (b)). IEKS-FOPD is the fastest convergent method, needing only $3$ iterations. Although the state estimation problem is relatively small scale, GN-PRS, GN-SBM, GN-FOPD and GN-ADMM are still time-consuming.

Like with linear tracking model, also in this simulation, we enlarge the time step count $T$ from $10^{2}$ to $10^{7}$, and plot the results in Fig. 7. When increasing the time step count, our proposed methods significantly outperform the batch methods with respect to CPU time. In particular, the PRS, SBM, FOPD and ADMM solvers with $10^{4}$ time steps take more time than our proposed methods with $10^{7}$ time steps, when $\lambda\in[0.1,0.5,1,2]$. In Table III, we further list the CPU times with different $\lambda$ when $T$ is varying. The performance benefit of the proposed methods is evident.

V-C Tomographic Reconstruction

In this section, we consider the application of the methodology to X-ray computed tomography (CT) imaging [49, 50]. First, we evaluate the performance of the proposed methods on real tomographic X-ray data of an emoji phantom measured at the University of Helsinki [51]. The dataset consists of 33-point time series of the X-ray sinogram of an emoji made of small squared ceramic stones. In the sequence, the emoji transforms from a face with closed eyes and a straight mouth to a face with smiling eyes and mouth. Typically, we have a sequence of square X-ray images of size $s\times s$ with $s=64,128$, which we are interested in reconstructing from low-dose observations taken from a limited number of angles. These low-dose observations can be modeled by the measurement matrix $\mathbf{H}_{t}$ which describes line integrals through the object (i.e., Radon transform).

Fig. 8 shows the CT reconstruction results for the emoji motion dataset, obtained by KS-ADMM. We set the parameters to $\lambda=10$, $\rho=1$, $k_{\max}=20$, and $n_{x}=4096$. The analysis operator $\mathbf{\Omega}_{t}$ consists of all the vertical and horizontal gradients (one step differences), which corresponds to so called TV regularization [21]. The number of measurements that correspond to $60$ or $30$ projections are $n_{y}=13020$ and $n_{y}=6510$, respectively. Although there is no ground truth to compare the qualitative results, we can observe the visual results from different numbers of projections. When the number of projections is $60$, the method provides good reconstruction results with $20$ iterations. We see that the 30-projection results suffer from the block artifacts as a consequence of the reduction in dose.

Furthermore, we validate the effectiveness of the proposed method on the real inhalation (iBH-CT) and exhalation (eBH-CT) breath-hold CT images, which was acquired as part of the National Heart Lung Blood Institute COPDgene study [52]. The dataset consists of 10 expiratory phase images of the segmented lung voxels. In detail, the parameters are $\lambda=1$, $\rho=0.1$, $k_{\max}=15$, $T=10$, $n_{x}=16384$, and the numbers of measurements are $n_{y}=13020$ and $n_{y}=6510$, corresponding to $60$ and $20$ projections. The ground truth and the reconstruction results are shown in Fig. 9. By visually comparing the results, we observe that moving from $60$ to $20$ projection provides much more drastic change. For example, some additional artifacts exist, but the result with the setting $n_{x}=16384$ and $n_{y}=6510$ is still very much acceptable (see the third column in Fig. 9). The results show that our methods still successfully preserve temporal information when the number of projections is $20$.

In the two experiments, we used a stationary Kalman filter and smoother to implement the optimization. We pre-computed all the gains before the iteration, which significantly speeded up the computations in tomographic reconstruction. We report CPU time (seconds) in Table IV. Table IV shows that KS-ADMM achieves significantly lower CPU time than the batch ADMM although the visual quality of all the reconstructions is equal. For example, in emoji motion dataset, when $n_{x}=16384$ , ADMM takes three time longer than our proposed method. In the lung dataset, when $n_{x}=16384$ and $n_{y}=6510$, KS-ADMM seems to be promising to provide computationally efficient reconstruction.

VI Conclusion

In this paper, we have presented two new classes of methods for solving state estimation problems. The estimation problem has been formulated as an (analysis) $L_{1}$-regularized optimization problem and the resulting problem has been solved by using the combinations of (iterated extended) Kalman smoother and variable splitting methods such as ADMM. The proposed approaches replace the batch solution for the state-update by using the smoother, which has a lower time-complexity than the batch solution. Furthermore, we have extended the proposed methods to a more general algorithmic framework, where the state-update is computed with the smoother. We have also established (local) convergence results for the novel KS-ADMM and IEKS-ADMM methods. In two different linear and nonlinear simulated cases, we have presented experimental results which show the efficiency of the smoother-based variable splitting optimization methods, especially when applied to large-scale or high-dimensional $L_{1}$-regularized state estimation problems. We also applied the methodology to a real-life tomographic reconstruction problem arising in X-ray-based computed tomography. Further work may explore a proper choice of the dual parameters using in KS / IEKS-based variable splitting methods, and discuss the convergence in the adaptive parameter settings.

Acknowledgements

The authors are grateful for the help of Zenith Purisha in preparing the computed tomography experiment and Zheng Zhao for useful comments on the manuscript.

Appendix A Proof of Lemma 1

We first prove for Case (a). By the first-order optimality condition of $\mathbf{x}$ subproblem, we have

[TABLE]

which implies that

[TABLE]

It follows that

[TABLE]

Then, if we assume that $\mathbf{\Omega}$ is full-row rank with $\mathbf{\Omega}\mathbf{\Omega}^{\top}\succeq\kappa_{a}^{2}\mathbf{I}$, we have

[TABLE]

By combining (54) and (55), we get

[TABLE]

Thus, for the $\bm{\eta}$-subproblem, we can use the primal variable to bound as follows:

[TABLE]

Since the $\mathbf{x}$-subproblem is $\mu_{x}$-strongly convex we have

[TABLE]

Similarly, since the $\mathbf{w}$-subproblem is convex, we have the following inequality:

[TABLE]

Thus, by combining (57), (LABEL:eq:strong_convex), and (59), we obtain:

[TABLE]

which will be negative provided by $\rho>{2L_{\theta_{1}}^{2}}/{\kappa_{a}^{2}\mu_{x}}$. Thus, when $\rho>\max\left(\frac{2L_{\theta_{1}}^{2}}{\kappa_{a}^{2}\mu_{x}},\rho_{0}\right)$, the result follows.

Next, we prove Case (b). In this case, we do not assume convexity of the $\mathbf{x}$-subproblem. For the $\bm{\eta}$-subproblem, we obtain

[TABLE]

Let $\mathbf{\Omega}$ be full-column rank with $\mathbf{\Omega}^{\top}\mathbf{\Omega}\succeq\kappa_{b}^{2}\mathbf{I}$, which gives

[TABLE]

For the $\mathbf{x}$-subproblem we get

[TABLE]

and by combining (59), (61), and (63), we get

[TABLE]

which will be nonnegative provided that $\rho>\frac{L_{\theta_{1}}}{\kappa_{b}^{2}}$.

Appendix B Proof of Lemma 2

The error between the local iterate $\mathbf{x}^{(i+1)}$ in the update and the minimizer $\mathbf{x}^{\star}$ satisfies the following recursion:

[TABLE]

Let $\mathbf{H}(\mathbf{x})$ be bounded by $\epsilon_{h}$, that is, $\|\mathbf{H}(\mathbf{x})\|\leq\epsilon_{h}$. We conclude that when $\epsilon_{h}\to 0$, the convergence is quadratic. Now let $\mathbf{J}^{\top}\mathbf{J}(\mathbf{x}^{(i)})\succeq\mu^{2}\mathbf{I}$. Then, linear convergence is obtained when the following condition is satisfied:

[TABLE]