A Modularized Efficient Framework for Non-Markov Time Series Estimation

TL;DR

This paper introduces a modular, efficient framework for estimating non-Markovian and non-Gaussian time series using a consensus ADMM approach, applicable to various complex signal processing problems.

Contribution

It proposes a novel, modular ADMM-based method for MAP estimation that handles non-Markov dynamics and non-Gaussian measurements with broad applicability.

Findings

Framework successfully estimates complex time series models.

Convergence guaranteed under broad convexity assumptions.

Demonstrated on neuroscience data with non-Gaussian and non-Markovian models.

Abstract

We present a compartmentalized approach to finding the maximum a-posteriori (MAP) estimate of a latent time series that obeys a dynamic stochastic model and is observed through noisy measurements. We specifically consider modern signal processing problems with non-Markov signal dynamics (e.g. group sparsity) and/or non-Gaussian measurement models (e.g. point process observation models used in neuroscience). Through the use of auxiliary variables in the MAP estimation problem, we show that a consensus formulation of the alternating direction method of multipliers (ADMM) enables iteratively computing separate estimates based on the likelihood and prior and subsequently "averaging" them in an appropriate sense using a Kalman smoother. As such, this can be applied to a broad class of problem settings and only requires modular adjustments when interchanging various aspects of the statistical model. Under broad log-concavity assumptions, we show that the separate estimation problems are convex optimization problems and that the iterative algorithm converges to the MAP estimate. As such, this framework can capture non-Markov latent time series models and non-Gaussian measurement models. We provide example applications involving (i) group-sparsity priors, within the context of electrophysiologic specrotemporal estimation, and (ii) non-Gaussian measurement models, within the context of dynamic analyses of learning with neural spiking and behavioral observations.

Tables

Table 1. Table I : Examples of common models. For the multiple modalities case, we define 𝐲 = ( 𝐲 ( 1 ) , … , 𝐲 ( J ) ) 𝐲 superscript 𝐲 1 … superscript 𝐲 𝐽 \mathbf{y}=(\mathbf{y}^{(1)},\dots,\mathbf{y}^{(J)}) to be a J 𝐽 J -tuple of simultaneous and conditionally independent observations, each with its own dimensionality and associated measurement model L ( j ) superscript 𝐿 𝑗 L^{(j)} .

Measurement Models - $L\left(𝐲\mid 𝐱\right)$
Linear Gaussian (LG)	${\sum }_{n=1}^N{\left|\left|{𝐲}_n-{\mathrm{𝐀𝐱}}_n+𝐛\right|\right|}_2^2\triangleq LG\left(𝐱\right)$
Sparse LG	$LG\left(𝐱\right)+{\left|\left|𝐱\right|\right|}_1$
Group sparse LG	$LG\left(𝐱\right)+{\sum }_{k=1}^K{\left({\sum }_{n=1}^N{x}_{k,n}^2\right)}^{\frac{1}{2}}$
Multiple Modalities	${\sum }_{j=1}^J{L}^{\left(j\right)}\left({𝐲}^{\left(j\right)}\mid 𝐱\right)$
System Models - $\varphi \left(𝐰\right)$
LG	${\sum }_{n=1}^N{\left|\left|{\mathrm{𝐂𝐰}}_n-𝐝\right|\right|}_2^2$
Sparse	${\left|\left|𝐰\right|\right|}_1$
Group sparse	${\sum }_{k=1}^K{\left({\sum }_{n=1}^N{w}_{k,n}^2\right)}^{\frac{1}{2}}$

Table 2. Table II : Examples of measurement model/system model pairings in previous works.

	Measurement Model	System Model
Kalman Smoother [4]	LG	LG
State Space Model of Learning [21]	Non-linear/ multiple modalities	Gaussian
Spectrotemporal Pursuit [11]	LG	Group sparse
Lasso-Kalman Smoother [13]	Group sparse	LG
Sparse States and Sparse Innovations [12]	Sparse	Sparse

Table 3. Table III : Performance metrics for the proposed method (ADMM), Fixed-Interval Smoother (FIS), and Sequential Monte Carlo (SMC) averaged over 50 trials with the Gaussian state-space model given by ( 19 ) and the state-space model with sparse variations given by ( 27 ).

	Gaussian State		Sparse Variations
	RMSE	Run Time (s)	RMSE	Run Time (s)
ADMM	$0.165$	$1.8$	$0.141$	7.0
FIS	0.168	2.6	0.181	$5.2$
SMC	0.188	53.5	0.186	105.7

Equations

$$f_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})=\prod _{n=1}^N{f}_{{\mathbf{Y}}_n\mid {\mathbf{X}}_n}({\mathbf{y}}_n\mid {\mathbf{x}}_n)$$

$${\mathbf{W}}_n=\{\begin{array}{ll}{\mathbf{X}}_1& \text{ n=1 }\\ {\mathbf{X}}_n-\mathbf{D}{\mathbf{X}}_{n-1}& \text{ n = 2,…,N }\end{array},$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}-\log f_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})-\log f_{\mathbf{X}}(\mathbf{x})$$

$$\hat{\mathbf{x}}$$

$$=\underset{\mathbf{x}}{\mathrm{argmin}}L(\mathbf{y}\mid \mathbf{x})+\beta \varphi (\mathcal{A}(\mathbf{x}))$$

$$L(\mathbf{y}\mid \mathbf{x})$$

$$\varphi (\mathbf{w})$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\left(\sum _{n=1}^N{L}_n({\mathbf{y}}_n\mid {\mathbf{x}}_n)\right)+\beta \varphi (\mathcal{A}(\mathbf{x}))$$

$$(\hat{\mathbf{x}},\hat{\mathbf{w}})=$$

$$s.t.$$

$$(\hat{\mathbf{x}},\hat{\mathbf{w}},\hat{\mathbf{z}})=$$

$$s.t.$$

$$\mathbf{w}=\mathcal{A}(\mathbf{z}).$$

$${\mathcal{L}}_{\rho }(\mathbf{x},\mathbf{w},\mathbf{z},\mathbf{\lambda },\mathbf{\alpha })=\sum _{n=1}^N{L}_n({\mathbf{y}}_n\mid {\mathbf{x}}_n)+\beta \varphi (\mathbf{w})+⟨\mathbf{\lambda },\mathbf{x}-\mathbf{z}⟩+⟨\mathbf{\alpha },\mathbf{w}-\mathcal{A}(\mathbf{z})⟩+\frac{\rho }{2}||\mathbf{x}-\mathbf{z}|{|}_F^2+\frac{\rho }{2}||\mathbf{w}-\mathcal{A}(\mathbf{z})|{|}_F^2$$

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

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

$${\mathbf{z}}^{(i+1)}$$

$${\mathbf{\lambda }}^{(i+1)}$$

$${\mathbf{\alpha }}^{(i+1)}$$

$${\mathbf{x}}^{(i+1)}=\underset{\mathbf{x}}{\mathrm{argmin}}\sum _{n=1}^N{L}_n({\mathbf{y}}_n\mid {\mathbf{x}}_n)+\frac{\rho }{2}||\mathbf{x}-{\tilde{\mathbf{x}}}^{(i)}|{|}_F^2$$

$${\mathbf{x}}_n^{(i+1)}=\underset{{\mathbf{x}}_n}{\mathrm{argmin}}{L}_n({\mathbf{y}}_n\mid {\mathbf{x}}_n)+\frac{\rho }{2}||{\mathbf{x}}_n-{\tilde{\mathbf{x}}}_n^{(i)}|{|}_2^2$$

$${\mathbf{x}}^{(i+1)}=\underset{\mathbf{x}}{\mathrm{argmin}}\sum _{n=1}^N{L}_n({\mathbf{y}}_n\mid {\mathbf{x}}_n)+\frac{\rho }{2}||{\mathbf{x}}_n-\mathbf{D}{\mathbf{x}}_{n-1}-{\tilde{\mathbf{x}}}_n^{\prime (i)}|{|}_2^2$$

$${\mathbf{w}}^{(i+1)}=\underset{\mathbf{w}}{\mathrm{argmin}}\beta \varphi (\mathbf{w})+\frac{\rho }{2}||\mathbf{w}-{\tilde{\mathbf{w}}}^{(i)}|{|}_F^2$$

$${\mathbf{z}}^{(i+1)}=\underset{\mathbf{z}}{\mathrm{argmin}}||\mathbf{z}-{\tilde{\mathbf{z}}}_{\mathbf{x}}^{(i)}|{|}_F^2+||\mathcal{A}(\mathbf{z})-{\tilde{\mathbf{z}}}_{\mathbf{w}}^{(i)}|{|}_F^2$$

$$\begin{array}{ll}0& =\hat{\mathbf{x}}-\hat{\mathbf{z}}\\ 0& =\hat{\mathbf{w}}-\mathcal{A}(\hat{\mathbf{z}})\end{array}\}PrimalFeasibility$$

$$\begin{array}{ll}0& \in \frac{\partial }{\partial \hat{\mathbf{x}}}L(\mathbf{y}\mid \hat{\mathbf{x}})+\hat{\mathbf{\lambda }}\\ 0& \in \frac{\partial }{\partial \hat{\mathbf{w}}}\beta \varphi (\hat{\mathbf{w}})+\hat{\mathbf{\alpha }}\\ 0& =\hat{\mathbf{\lambda }}+\tilde{\mathcal{A}}(\hat{\mathbf{\alpha }})\end{array}\}DualFeasibility$$

$$\begin{array}{ll}{r}_1^{(i)}& ={\mathbf{x}}^{(i)}-{\mathbf{z}}^{(i)}\\ r_2^{(i)}& ={\mathbf{w}}^{(i)}-\mathcal{A}({\mathbf{z}}^{(i)})\end{array}\}PrimalResiduals$$

$$\begin{array}{ll}{s}_1^{(i)}& =\rho \tilde{\mathcal{A}}({\mathbf{w}}^{(i)}-{\mathbf{w}}^{(i-1)})\\ s_2^{(i)}& =\rho ({\mathbf{z}}^{(i)}-{\mathbf{z}}^{(i-1)})\end{array}\}DualResiduals$$

$${ϵ}_1^{pri}$$

$${ϵ}_2^{pri}$$

$${ϵ}_1^{dual}$$

$${ϵ}_2^{dual}$$

$$X_n=\kappa X_{n-1}+\gamma +V_n$$

$${\color[rgb]{0,0,0}\mathbb{P}}(B_{n}=b_{n}\mid X_{n}=x_{n})=p_{n}^{b_{n}}(1-p_{n})^{1-b_{n}}$$

$$p_{n}=\frac{\exp(\nu+\eta x_{n})}{1{\color[rgb]{0,0,0}+}\exp(\nu+\eta x_{n})}$$

$$R_n\sim \mathcal{N}(\psi +\omega X_n,{\sigma }_R^2)$$

$$\displaystyle-\log{\color[rgb]{0,0,0}\mathbb{P}}(\mathbf{S}_{n}=\mathbf{s}_{n}$$

Full text

I Introduction

We consider the problem of estimating a latent time series based on an underlying dynamic model and noisy measurements. Such a problem appears in a variety settings, including (but certainly not limited to) tracking [1], medical imaging [2], and video denoising [3]. Given the broad applicability of this problem formulation, the underlying models that are used inevitably become increasingly complex.

Certain scenarios are well studied, such as the case of a linear system with Gaussian noise, where it is well known that the maximum a-posteriori (MAP) point estimate can be obtained using a Kalman smoother (KS) [4]. When introducing non-linearities, alternatives include the extended Kalman filter (EKF), which relies on linear approximations, as well as the unscented Kalman filter (UKF) [5] and Particle Filter (PF) [6], which use sample based techniques. While the EKF and UKF are well suited for a broad class of problems, they are not well suited for models with non-Gaussian noise. This is problematic for the increasingly popular problem of incorporating sparsity inducing models to latent signal estimation. These problems include exploiting sparsity in the underlying signal [7, 8, 9, 10] in addition to exploiting sparsity in the signal dynamics [11, 12, 13]. While some of these methods utilize $\ell_{1}$-regularization to enforce sparsity at a local level and enable causal prediction, there is often knowledge of global structures, such as those favored by the group lasso [14], that dictate a need for batch-wise estimation. In such cases, the desired estimation problem deviates from the classical state estimation problem in that the underlying signal is no longer Markov. In such a scenario, there is no clear extension to the EKF, UKF, or PF that may be utilized to address the non-Markovicity of the underlying signal.

The broad scope of the problem in question dictates a need for a systematic approach to latent time series estimation for a variety of measurement models and system models. Furthermore, a solution framework that can compartmentalize these two models facilitates interchangeability and allows new regularization techniques to be easily incorporated to an estimation procedure.

We develop a framework using the alternating direction method of multipliers (ADMM) [15] that, under mild (i.e. log-concavity) assumptions, yields the MAP estimate for problems with non-Markov latent variables and/or nonlinear observations. While ADMM has been utilized to decompose specific dynamic systems into simpler subproblems [16, 13], our approach applies to arbitrary log-concave dynamic models. In particular, we utilize auxiliary variables to enable a solution involving iterative updates to three modules, one that pertains to the measurement model, another that pertains to the prior distribution on the latent signal, and a third that is a Kalman smoother. As such, our framework enables various sparsity models to be easily applied to the signal and/or dynamics with adjustments only required to the corresponding module. We demonstrate implementation of the framework in two distinct applications, namely latent state estimation and spectrotemporal estimation. We show that in the case of state estimation, our method outperforms a fixed interval smoother and particle filter for two state-space models coupled with non-Gaussian observations. In the case of spectrotemporal estimation, we demonstrate the efficacy of our method when using non-Markov priors. The proposed method yields an intuitive approach to latent process estimation with iterative use of a Kalman smoother in tandem with standard convex optimization techniques. We provide a mathematical justification for the intuition by proving that our approach guarantees convergence to the MAP solution under the same relatively mild conditions that apply to general ADMM approaches. Finally, we provide software to enable the reader to reproduce the results of this paper and to easily apply the framework to novel models111An implementation of the proposed framework can be found in the following publicly available GitHub repository: https://github.com/gabeschamberg/nonmarkov-timeseries-estimation/releases/tag/v1.1. This repository includes the iPython [17] notebooks that were used to generate Figs. 2 and 3 and Table III. Our contributions may be summarized as follows:

• •

We present an efficient iterative solution framework for latent time series estimation with a guarantee of convergence to the MAP estimate under mild log-concavity assumptions.

• •

In the presence of non-Linear, non-Gaussian measurement models, our method does not require a Gaussian approximation, unlike KS variants, and is more efficient than Sequential Monte Carlo (SMC) methods.

• •

Our framework accommodates non-Markov signals despite there being no clear method for adapting EKF, UKF, and SMC methods for such a scenario, particularly when the prior applies to highly non-linear functions of the latent process, such as a singular value decomposition.

• •

Through the use of auxiliary variables, the ADMM solution to our reformulated MAP estimation problem is modular, with the observation and system models in disjoint modules that are unified by a Kalman smoother.

The paper is structured as follows: Section II provides the general formulation of the problem we are solving in addition to a brief review of relevant work solving specific instances of the problem. Section III details a novel systematic approach for solving the MAP estimation problem in its general form. Section IV demonstrates the capabilities of the framework through implementation on two existing problems. Section V concludes the paper with a discussion of the results and future work.

II Problem Formulation

II-A Notation

While it is intended that the notation is presented unambiguously, we here present some notational conventions. Bold letters are used to represent vectors and matrices, whereas non-bold letters represent scalars. Subscripts are used for indexing scalar elements of a vector, or columns of a matrix. A double subscript is used to specify scalar elements of a matrix. For example, $x_{n}$ gives the $n$th element of a vector $\mathbf{x}$, $\mathbf{x}_{n}$ gives the $n$th column of a matrix $\mathbf{x}$, and $x_{n,m}$ gives the $m$th row of the $n$th column of a matrix $\mathbf{x}$. Capital/lowercase letter pairs represent either random variable/realization pairs or total count/index pairs. For example, we may have that $\mathbf{x}_{n}$ gives a specific value of the random vector $\mathbf{X}_{n}$, which is the $n$th column of a random matrix $\mathbf{X}$ with $N$ columns in total. We let $f$ and $p$ denote probability density functions (pdfs) and probability mass functions (pmfs), respectively. Various joint and conditional pdfs and pmfs are made clear by their subscripts. For example, the pdf of $X$ given $Y=y$ is $f_{X|Y}(\cdot|y)$. We let $\mathbb{R}$ denote the space of real numbers, $\mathbb{R}_{+}$ denote the non-negative reals, $\mathbb{R}^{A\times B}$ denote the space of $A$ by $B$ real valued matrices, and $\mathbb{R}^{AB}$ denote the space of real valued vectors of length $A$ times $B$.

II-B Problem Setup

Let $\mathcal{X}$ and $\mathcal{Y}$ be measurable spaces and $N$ be the length of time series pertaining to the latent process $\mathbf{X}\in\mathcal{X}^{N}$ and observed process $\mathbf{Y}\in\mathcal{Y}^{N}$. Unless otherwise specified, we assume $\mathcal{X}=\mathbb{R}^{K}$ and $\mathcal{Y}=\mathbb{R}^{P}$ where $K$ is the dimension of the latent process at any time, and $P$ is the dimension of the observation process at any time. As such, $\mathbf{X}\in\mathbb{R}^{K\times N}$ is the latent time series we wish to estimate and $\mathbf{Y}\in\mathbb{R}^{P\times N}$ is the collection of noisy observations. Furthermore, assume that these observations are conditionally independent given the underlying time series:

[TABLE]

where $f_{\mathbf{Y}\mid\mathbf{X}}$ is the likelihood of the entire collection of observations given the entire latent time series and $f_{\mathbf{Y}_{n}\mid\mathbf{X}_{n}}$ is the likelihood of a single observation given the corresponding element of the latent time series.

Next, define the latent signal’s dynamics (or system behavior) in terms of $\mathbf{W}\in\mathbb{R}^{K\times N}$ for which

[TABLE]

where $\mathbf{D}\in\mathbb{R}^{K\times K}$ is a transition matrix and $\mathbf{W}_{n}\in\mathbb{R}^{K}$ and $\mathbf{X}_{n}\in\mathbb{R}^{K}$ represent the $n$th columns of $\mathbf{W}$ and $\mathbf{X}$, respectively. For compactness we write this as $\mathbf{W}=\mathcal{A}(\mathbf{X})$, where $\mathcal{A}$ represents a linear operator that is fully defined by $\mathbf{D}$. We assume that $\mathbf{W}$ is distributed according to a known prior pdf $f_{\mathbf{W}}(\mathbf{w})$. Note that this framework includes, for the special case of $\mathbf{W}_{n}=\mathbf{X}_{n}-\mathbf{X}_{n-1}$ and $\mathbf{W}_{n}\sim\mathcal{N}(\mu_{n},\Sigma_{n})$ are independent Gaussian random vectors for $n=2,\dots,N$, the well-studied scenario in which the underlying time series $\mathbf{X}$ is a Gauss-Markov process.

Here, we consider the problem of finding the maximum a posteriori estimate:

[TABLE]

where $-\log f_{\mathbf{Y}\mid\mathbf{X}}(\mathbf{y}\mid\mathbf{x})$ is the negative log-likelihood and $-\log f_{\mathbf{X}}(\mathbf{x})$ is the negative log-prior. We note that because $\mathbf{W}$ is a linear function of $\mathbf{X}$, we have $f_{\mathbf{X}}(\mathbf{x})\propto f_{\mathbf{W}}(\mathcal{A}(\mathbf{x}))$. This relationship indicates that knowing a prior on either $\mathbf{X}$ or $\mathbf{W}$ induces a prior on the other. Thus, we can equivalently rewrite our problem as:

[TABLE]

with $\beta\in\mathbb{R}_{+}$ and where we define the measurement model $L:\mathbb{R}^{N\times K}\rightarrow\mathbb{R}$ and system model $\phi:\mathbb{R}^{N\times K}\rightarrow\mathbb{R}$ as:

[TABLE]

The inclusion of $\beta$ in (4) is to facilitate the cases when the system model is only known up to a proportionality constant or when $\phi$ is a regularizer used to exploit a desired dynamic characteristic of the latent signal (as opposed to representing the true distribution of $\mathbf{W}$). In either of these cases $\beta$ is interpreted as a tuning parameter used to control the extent to which the system model is weighted (as in $\lambda$ throughout [14]).

Throughout this paper, we will interchangeably use the names log-likelihood/measurement model in reference to $L$, and log-prior/system model/dynamic model in reference to $\phi$. Due to the assumption that observations are conditionally independent given the state variables, the measurement model can be decomposed into a sum over $N$ measurements, each depending on the state variable at a single time instance:

[TABLE]

where $L_{n}(\mathbf{y}_{n}\mid\mathbf{x}_{n})\coloneqq-\log f_{\mathbf{Y}_{n}\mid\mathbf{X}_{n}}(\mathbf{y}_{n}\mid\mathbf{x}_{n})$. It should be noted that the problem presented in (7) is made difficult by the second term. In particular, imposing a prior on the differences of the underlying time series prevents separability across the $N$ time points. Furthermore, by allowing for non-Markov models, it is possible to have models that do not allow the second term to be separated into terms each containing only $\mathbf{x}_{n}$ and $\mathbf{x}_{n-1}$ for each $n=1\dots N$. In the following section, we present a framework for efficiently solving problems in the form of (7) for a broad class of measurement models $L$ and system models $\phi$.

II-C Related Work

Works related to our proposed method include both the investigation of new algorithms for estimating latent time series and the creation/application of new time series models. Notably, the Kalman smoother [4] and its variants [5, 6] provide structured approaches to estimating latent signals in a subset of problems with dynamical system models and noisy measurements. While the Kalman smoother is MAP optimal for the very specific case of a linear system with Gaussian noise, its non-linear variants do not guarantee optimality and do not offer solutions for a comprehensive class of measurement and system models. In particular, there has been growing interest in models exploiting the sparsity of states and/or dynamics of signals [7, 8, 9, 11, 12, 13], which in many cases do not lend themselves to solutions via the existing Kalman smoother variants.

For such sparsity-inducing models, existing causal estimators are often heuristic extensions of the Kalman filter, such as $\ell_{1}$-regularized Kalman filter updates [12] and tracking a belief of the support set [7]. Causal estimation is made particularly challenging for the models that are non-Markov in nature. As such, the aforementioned causal estimators lack performance guarantees. Existing batchwise solutions utilize a Kalman smoother to solve the updates for a particular iterative algorithm, such as IRLS for group sparse dynamics [11] and ADMM for group sparse states [13]. In the latter example, their non-consensus formulation of ADMM is reliant upon the choice of a Gaussian system model.

In addition to the Kalman smoother variants, sample based methods such as Markov chain Monte Carlo (MCMC) and SMC are viable options for latent time series estimation. While these methods can accommodate non-linear and non-Gaussian models [18] and can simultaneously estimate the state and model parameters [19, 20], they are often computationally prohibitive. Furthermore, these methods do not have a straightforward extension to non-Markov and non-linear priors such as the $\ell_{1}/\ell_{2}$ and nuclear norm priors (see Remark 2).

Here we propose a generalized framework for obtaining the MAP estimate in many of the aforementioned problems in a batchwise manner. Tables I and II show the models used in some of these problems and serve to illustrate the primary contribution of our framework, namely that for a given problem, the solution is modular in that the choice of measurement model can be made independently of the system model without requiring a complete rederivation of the solution.

III Modular MAP Estimation Framework

The alternating direction method of multipliers (ADMM) allows large global problems to be decomposed into smaller subproblems whose solutions can be coordinated to achieve the global solution. ADMM offers an iterative solution of the dual problem that has the decomposability of dual descent in addition to the convergence guarantees of the method of multipliers, which hold under fairly mild conditions. While the details of dual optimization and ADMM are omitted here, they can be found in [22] and [15], respectively.

We begin by reformulating (7) to create separability in the objective function by including $\mathbf{w}$ as an optimization variable and introducing a constraint to preserve the relationship between $\mathbf{x}$ and $\mathbf{w}$:

[TABLE]

The optimization problem given by (8) can be solved using ADMM, and would yield a solution that enables the measurement model and penalty function to be addressed in independent subproblems. However, when using the above formulation, the update equations yielded by the ADMM algorithm would require one of the aforementioned approximate or sample-based methods for non-Gaussian measurement models (see Remark 1).

We use a variant of ADMM known as consensus ADMM and construct a modular solution framework shown in Fig. 1 that only requires making local adjustments to the solution when modifying the measurement model ($L$), penalty function ($\phi$), or transition model ($\mathcal{A}$). This is accomplished by introducing an auxiliary variable $\mathbf{z}\in\mathbb{R}^{K\times N}$ to achieve separability (of $\mathbf{x}$ and $\mathbf{w}$) in the constraints as well as the objective function:

[TABLE]

The optimization problem given by (9) is termed the consensus formulation, and $\mathbf{z}$ the consensus variable. By introducing this variable, our iterative updates with respect to the measurement model and penalty function are not only independent of each other, but are also independent of the transition model determined by $\mathcal{A}$.

The first step in solving (9) using ADMM requires generating the augmented Lagrangian:

[TABLE]

where $\boldsymbol{\lambda}\in\mathbb{R}^{K\times N}$ and $\boldsymbol{\alpha}\in\mathbb{R}^{K\times N}$ are Lagrange multipliers, $\langle\cdot{,}\cdot\rangle$ is the Frobenius inner product, $\lvert\lvert\cdot\rvert\rvert_{F}$ is the matrix Frobenius norm, and $\rho\in\mathbb{R}_{+}$ is the penalty parameter for the augmented Lagrangian. Note that in the case where $\rho=0$, the augmented Lagrangian is equivalent to the standard (unaugmented) Lagrangian.

Given the augmented Lagrangian, the ADMM solution is obtained by iteratively alternating between minimization with respect to the primal variables ($\mathbf{x}$, $\mathbf{w}$ and $\mathbf{z}$) and performing gradient ascent on the Lagrange multipliers. These iterations represent a trade off between finding a solution that minimizes the cost function in (9) while ensuring that the Lagrange multipliers are such that the dual function of (9) is increasing in $i$ and thus ensuring the constraints are satisfied. Letting $\mathbf{x}^{(i)}$ represent the estimate of $\mathbf{x}$ after $i$ iterations (similarly for $\mathbf{w}^{(i)}$, $\mathbf{z}^{(i)}$, $\boldsymbol{\lambda}^{(i)}$, and $\boldsymbol{\alpha}^{(i)}$), each iteration of ADMM is composed of the following updates [15, Sec. 3.1]:

[TABLE]

By fixing all but one variable in each update, the objective functions can be simplified by dropping the terms in (10) that do not contain the optimization variable for the corresponding update. As a result, when updating with respect to the measurement model $L$ and the system model $\phi$, we only need to consider the model corresponding to that update and an $\ell_{2}$-norm proximal operator [23] that ensures the update is moving in the appropriate direction to achieve a global consensus. This inclusion of the proximal operators in the augmented Lagrangian enables the use of ADMM with non-smooth objective functions [23, Sec 4.4]. Then, updating of the consensus variable involves “centering” it such that it gives equal representation to our current estimates based on the measurements and our estimates based on the system dynamics. In this sense, our ADMM framework yields a mathematical justification for a very intuitive approach, namely, iteratively finding the best estimate based on measurements, finding the best estimate based on dynamics, and “averaging” the two in the appropriate sense. This viewpoint will be made clearer in the following sections where we detail the specific update equations.

III-A Measurement Model Update

When updating with respect to the measurement model, only terms containing $\mathbf{x}$ in the augmented Lagrangian must be considered. To simplify notation, we will consider the scaled form of the update equations [15, Sec. 3.1.1], which can be obtained by combining the appropriate linear and quadratic terms in (10) by completing the square:

[TABLE]

where $\tilde{\mathbf{x}}^{(i)}\coloneqq\mathbf{z}^{(i)}-\nicefrac{{\boldsymbol{\lambda}^{(i)}}}{{\rho}}$ is fixed within the scope of this update. Details for deriving the scaled form of the update can be found in Appendix A. Given that the squared Frobenius norm can be decomposed to the sum of squared $\ell_{2}$ norms, we note that the measurement model update is separable over $n$, meaning that we can solve for $\mathbf{x}_{n}^{(i+1)}$ for each $n=1,\dots,N$ independently:

[TABLE]

where $\tilde{\mathbf{x}}_{n}^{(i)}\coloneqq\mathbf{z}_{n}^{(i)}-\nicefrac{{\boldsymbol{\lambda}_{n}^{(i)}}}{{\rho}}$.

Remark 1**.**

Note that the ability to separate each of the $N$ updates is a result of the inclusion of the consensus variable. Excluding this variable would require that the dynamics be considered in the update of the measurement model:

[TABLE]

where $\tilde{\mathbf{x}}_{n}^{\prime(i)}\coloneqq\mathbf{w}_{n}^{(i)}-\nicefrac{{\boldsymbol{\gamma}_{n}^{(i)}}}{{\rho}}$, $\mathbf{x}_{0}\coloneqq 0$, and $\boldsymbol{\gamma}$ represents the single Lagrange multiplier that would be required in solving (8) using ADMM. Requiring that the dynamics of the underlying time series be included in the measurement model update prohibits solving for $\mathbf{x}_{n}^{(i)}$ independently across $N$. Thus, using ADMM in this fashion does not offer any simplifications over traditional approaches for non-Gaussian measurement models. As such, incorporation of the consensus variable not only enables faster processing by allowing each update to be parallelized across $N$, but it allows the framework to be applied in a straightforward, non-approximate manner to a broad class of measurement models.

It should be noted that while we assume conditional independence of the observations given the latent time series, one can revert to the update in (12) for the case when the observations are correlated. In this case the ability to parallelize across $n$ is lost, but the ability to ignore system dynamics is preserved (i.e. the optimization problem in (12) still does not depend on $\phi$).

III-B System Model Update

In the system model update, only terms in (10) that contain $\mathbf{w}$ must be included. Again, we consider the scaled form:

[TABLE]

where $\tilde{\mathbf{w}}^{(i)}\coloneqq\mathcal{A}(\mathbf{z}^{(i)})-\nicefrac{{\boldsymbol{\alpha}^{(i)}}}{{\rho}}$. In this form we can clearly interpret the system model update as finding a new collection of latent variable transitions $\mathbf{w}^{(i+1)}$ that is both representative of our system model $\phi$ and proximal to the appropriately scaled current consensus on the transitions $\tilde{\mathbf{w}}^{(i)}$.

The key observation is that this framework imposes no restrictions as to whether or not our underlying signal is Markov. In the case where the signal is indeed Markov, then $\mathbf{w}_{n}^{(i+1)}$ would be updated independently over $n$, but in general we do not assume this is the case. This provides the ability to impose batch-level structures on the dynamics of the signal. Furthermore, we note that the nature of the proximal operator enables closed form solutions when $\phi$ is chosen to be a number of common sparsity inducing priors. In particular, because the proximal operator is not multiplying $\mathbf{w}$ by a non-orthonormal matrix, the $\ell_{1}$, group sparse, and nuclear norm priors all offer soft-thresholding solutions [24]. Furthermore, we note that for a fixed $K$, the complexity of the soft-thresholding solutions for the $\ell_{1}$ and group sparse priors scale linearly with $N$ per iteration. The nuclear norm prior, however, requires a singular value decomposition (SVD), and thus scales quadratically with $N$ per iteration [25]. Similarly, for a fixed $N$, the same scaling factors apply to $K$. It should be noted however, that if increasing $N$ and $K$, the complexity of the SVD will scale quadratically with $\max\{K,N\}$ and cubically with $\min\{K,N\}$.

III-C Consensus Update

Updating the consensus variable depends on neither the measurement model nor the system model. We can think of this step as averaging our current estimates of our signal based on measurements $\mathbf{x}^{(i+1)}$ and based on dynamics $\mathbf{w}^{(i+1)}$:

[TABLE]

where $\tilde{\mathbf{z}}_{\mathbf{x}}^{(i)}\coloneqq\mathbf{x}^{(i+1)}+\nicefrac{{\boldsymbol{\lambda}^{(i)}}}{{\rho}}$ and $\tilde{\mathbf{z}}_{\mathbf{w}}^{(i)}=\mathbf{w}^{(i+1)}+\nicefrac{{\boldsymbol{\alpha}^{(i)}}}{{\rho}}$. Note that given the nature of the linear operator $\mathcal{A}$, (15) can always be solved efficiently using a Kalman smoother.

This step clarifies the notion of “averaging” the current estimates $\mathbf{x}^{(i+1)}$ and $\mathbf{w}^{(i+1)}$. By framing our problem from a consensus ADMM perspective, we can carve out various elements of the model and delegate them to independent updates. Then, given the nature of the relationship between the signal $\mathbf{x}$ and the dynamics $\mathbf{w}$, establishing consensus between the two estimates is a Kalman smoothing problem regardless of the measurement and system models. This is a result of the use of $\ell_{2}$-norms in the augmented Lagrangian, which can be thought of as representing Gaussian noise with identity covariance. In other words, at each iteration $i$, the consensus update is a Kalman smoothing problem where each of our measurements are given by $\tilde{\mathbf{z}}_{\mathbf{x}}^{(i)}$ and each of our predictions are given by $\tilde{\mathbf{z}}_{\mathbf{w}}^{(i)}$. In this sense, the consensus update gives equal weight to the current iterates of our measurement and system estimates. This follows from the fact that the log-likelihood and log-prior have their own uncertainty terms that dictate how far the updates $\mathbf{x}^{(i+1)}$ and $\mathbf{w}^{(i+1)}$ can deviate from the consensus in their respective updates, namely measurement noise and the tuning parameter $\beta$. We note that because both terms in (15) can be thought of as representing Gaussian noise with identity covariance and the transition model $\mathcal{A}$ is invariant over iterations $i$, all matrix inversions required by the Kalman smoother can be precomputed. As a result, each iteration requires on the order of $N$ matrix multiplications.

III-D Convergence

Next we consider the practical and theoretical convergence of the proposed framework. To begin, we present the optimality conditions and the means with which we can in practice implement convergence checks. The derivations are omitted, as they closely follow Section 3.3 of [15]. The optimality conditions for the proposed framework are given by:

[TABLE]

where $\nicefrac{{\partial}}{{\partial\cdot}}$ is the subgradient operator (or gradient when defined, in which case $\in$ becomes an equality) and where $\tilde{\mathcal{A}}(\mathbf{A})_{n}=\mathbf{A}_{n}-\mathbf{D}^{T}\mathbf{A}_{n+1}$ for $n=1,\dots,N-1$ and $\tilde{\mathcal{A}}(\mathbf{A})_{N}=\mathbf{A}_{N}$ for $\mathbf{A}\in\mathbb{R}^{K\times N}$. The primal feasibility conditions ensure that our $\hat{\mathbf{z}}$ preserves the desired relationship between $\hat{\mathbf{x}}$ and $\hat{\mathbf{w}}$, and the dual feasibility conditions serve the purpose of ensuring that the optimal Lagrange multipliers are such that $\hat{\mathbf{x}}$ and $\hat{\mathbf{w}}$ jointly minimize $L$ and $\phi$.

Using these optimality conditions, we can derive the primal and dual residuals:

[TABLE]

where primal feasibility is achieved when $r_{j}^{(i)}=0$ and dual feasibility is achieved when $s_{j}^{(i)}=0$ for all $j\in\{1,2\}$. In practice, we declare the algorithm converged when $\lvert\lvert r_{j}^{(i)}\rvert\rvert_{F}\leq\epsilon_{j}^{pri}$ and $\lvert\lvert s_{j}^{(i)}\rvert\rvert_{F}\leq\epsilon_{j}^{dual}$ for all $j\in\{1,2\}$, with the thresholds given by:

[TABLE]

where $\epsilon^{rel}$ (relative tolerance) and $\epsilon^{abs}$ (absolute tolerance) are small positive parameters.

In general, ADMM does not guarantee convergence for more than two optimization variables [26]. As such, it is not immediately clear that our ADMM framework would guarantee convergence given that it optimizes over $\mathbf{x}$, $\mathbf{w}$, and $\mathbf{z}$. As it turns out, for the particular version of consensus ADMM that we are proposing, we can guarantee convergence under the same mild conditions required in standard ADMM.

Theorem 1**.**

Given an observation $\mathbf{y}$, when $L(\mathbf{y}\mid\cdot)$ and $\phi(\cdot)$ are closed, proper, and convex functions, the ADMM algorithm given by (10) and (11) converges to the solution of (9), i.e. $(\mathbf{x}^{(i)},\mathbf{w}^{(i)},\mathbf{z}^{(i)})\rightarrow(\hat{\mathbf{x}},\hat{\mathbf{w}},\hat{\mathbf{z}})$ as $i\rightarrow\infty$.

The proof of Theorem 1 is based on a consensus ADMM formulation presented in section 5 of [27] and is given in detail in Appendix B.

IV Applications

IV-A State-Space Model of Learning

We begin by demonstrating how the ADMM framework can be applied to a problem with a highly non-linear multimodal measurement model. In the state-space model of learning [28], the system model is a traditional state-space Gauss-Markov process, where the state represents an unobservable cognitive state that represents a subject’s ability to perform a task over time. The corresponding measurement model provides a statistical relationship between the underlying state and the observed task performance for a given trial.

We define $\mathbf{X}\in\mathbb{R}^{1\times N}$ to be the cognitive state (with $K=1$), where $N$ represents the number of trials conducted. The system model is given by:

[TABLE]

where $\kappa\in[0,1]$ is a forgetting factor, $\gamma\in\mathbb{R}_{+}$ is a positive bias that represents a tendency for the cognitive state to increase with time, and $V_{n}\sim\mathcal{N}(0,\sigma_{V}^{2})$ is noise in the system model.

Using the state-space model of learning pertaining with multiple behavioral and neurophysiological measures, we assume that each of the $N$ trials has an associated binary success/failure outcome, a reaction time, and neural spiking behavior. As such, each observation is given by a triplet $\mathbf{Y}_{n}=(B_{n},R_{n},\mathbf{S}_{n})\in\{0,1\}\times\mathbb{R}\times\{0,1\}^{J}$, where $B_{n}$ is a binary random variable indicating whether or not the trial was completed successfully, $R_{n}$ is the log of the subject’s reaction time to complete the task, and $\mathbf{S}_{n}$ is a length $J$ point process that indicates whether or not there was neural spiking activity in each discrete $\Delta t$ time window.

Each of the three observation modalities is associated with an appropriate statistical model. First, the binary success/failure outcomes obey a Bernoulli probability model:

[TABLE]

where $p_{n}$ is given by a logistic function that maps the cognitive state between [math] and $1$:

[TABLE]

where $\nu,\eta\in\mathbb{R}$ are model parameters.

Next, the reaction time obeys a log-normal probability model, with:

[TABLE]

where $\psi\in\mathbb{R}$ is the estimated initial log reaction time, $\omega\in\mathbb{R}_{-}$ is negative to ensure that the reaction time tends to decrease with an increasing cognitive state and $\sigma_{R}^{2}$ represents the level of stochasticity in the relationship between the cognitive state and reaction time.

Lastly, the neural spiking activity is modeled as a point process (as in equation 2.6 of [29]), with the negative log-probability of a given set of spikes given by:

[TABLE]

where $s_{n,j}\in\{0,1\}$ is the $j^{th}$ bit of $\mathbf{s}_{n}$ and $\log\Lambda$ is the conditional intensity function, given by a generalized linear model [30]:

[TABLE]

where $\xi\in\mathbb{R}$ gives a base intensity level, $a\in\mathbb{R}$ determines the effect of the cognitive state on the spiking intensity, and $\mathbf{c}=(c_{1},\dots,c_{M})\in\mathbb{R}^{M}$ accounts for the refractory period in neural spiking, i.e. the fact that it is unlikely to see spiking activity in neighboring bins. The point process model given by (23) represents a discrete approximation of the negative log-likelihood for an inhomogeneous Poisson process where the rate in trial $n$ and time $j$ is $\Lambda_{n,j}$.

Next we adapt the state-space model of learning to the ADMM framework. We begin by considering the negative log-likelihood of the observations given the underlying cognitive state. We note that not only are the observations temporally conditionally independent given a sequence of cognitive states, but each of the three observations within a trial is conditionally independent given the cognitive state corresponding with that trial:

[TABLE]

where the negative log-likelihoods $L_{B_{n}}\coloneqq-\log p_{B_{n}\mid X_{n}}$, $L_{R_{n}}\coloneqq-\log f_{R_{n}\mid X_{n}}$, and $L_{\mathbf{S}_{n}}\coloneqq-\log p_{\mathbf{S}_{n}\mid X_{n}}$ are defined to be the negative log of the appropriate pdf/pmf corresponding with the respective observations. It is important to note that $L$ is indeed convex. Considering this is not immediately obvious, it is shown in Appendix D.

Next we consider the system model. By defining $W_{n}=X_{n}-\kappa X_{n-1}=\gamma+V_{n}$ with $W_{0}=X_{0}$, we get that $W_{n}\sim\mathcal{N}(\gamma,\sigma_{V}^{2})$, i.e. each $W_{n}$ is distributed iid Gaussian. Thus, our negative log-prior is given by:

[TABLE]

where $\mathcal{N}(x;\mu,\sigma^{2})$ gives the value of a normal distribution with mean $\mu$ and variance $\sigma^{2}$ evaluated at $x$. Additionally, under this definition of $\mathbf{W}$ we get that the transition matrix $\mathbf{D}$ is in fact just a scalar, namely $\kappa\in\mathbb{R}$.

Plugging $L$, $\phi$, and $\mathcal{A}$ into equations (12), (14), and (15), we obtain the update equations for solving the state-space model of learning problem. Beginning with the measurement model update, as a result of its separability across trials, each update decomposes into $N$ univariate convex minimization problems. As such, these $N$ problems can be solved in parallel using a convex solver such as CVX [31]. For the system model update, we note that because (14) is separable over $n=1,\dots,N$, the update is reduced to $N$ quadratic minimizations that can be solved in closed form. Given that the density for $\mathbf{W}$ is assumed to be fully known, we set the tuning parameter $\beta=1$. The details of these updates can be found in Appendix C.

We demonstrate the state-space model of learning solution on simulated data with $N=25$, using parameters from section V-A of [21]. The proposed method is compared with the fixed-interval smoother (FIS) detailed in [21] and a sequential Monte Carlo (SMC) method. In particular, we develop a particle smoother using the forward-filtering backward-sampling technique with systematic resampling at each step [32]. For the ADMM method, we set $\rho=30$ and limit the procedure to $25$ iterations, i.e. $\hat{\mathbf{x}}\coloneqq\mathbf{x}^{(25)}$. For the SMC method, we use 100 particles. In Table III we look at the average root-mean-square error (RMSE) and average runtime for each method over 50 trials, where for a given realization $\mathbf{x}$ and a given estimate $\hat{\mathbf{x}}$, $\operatorname{RMSE}(\hat{\mathbf{x}})=\nicefrac{{\lvert\lvert\hat{\mathbf{x}}-\mathbf{x}\rvert\rvert_{2}}}{{\sqrt{N}}}$. We note that the proposed method is both most efficient and most accurate in the RMSE sense. While the SMC method would presumably benefit from a larger number of particles, we see that even with limited samples, it is very computationally intensive. While the difference in RMSE is negligible across all 3 methods, it is worth noting that each method obtains a fundamentally different estimate. To be specific, the proposed method gives the MAP estimate in the limit of large iterations, while the other methods yield conditional expectations of the states given the entire observation sequence. In the case of the FIS, the estimate is the conditional expectation under a Gaussian approximation of the posterior. The SMC method, on the other hand, yields the true conditional expectation in the limit of large particle count.

It should be noted that in the case of a Gaussian state space, the problem formulations given by (8) and (9) are nearly equivalent. In particular, it is possible to omit the consensus variable and modify the constraint such that $\mathbf{W}=\mathbf{X}$. In such a scenario, the measurement model update would remain the same and the system model update would be solvable with a Kalman smoother. Thus, we further demonstrate the utility of our method by considering a second state-space model with sparse variations where such an approach is not possible. We simulate a state-space model with sparse variations by defining $X_{n}=X_{n-1}+V_{n}$ with $V_{n}$ obeying a commonly used sparsity inducing mixture model [33]:

[TABLE]

where $p\in[0,1]$ is a probability, $\sigma\in\mathbb{R}_{+}$ is a positive constant, and we define $U_{n}\sim\chi_{2}^{2}$ as i.i.d. Chi-Squared random variables with two degrees of freedom. This model represents a scenario supported by neurophysiological findings [34, 35] wherein infrequent, discontinuous changes in neural activity arise.

We again conduct 50 trials, setting $N=50$, $p=0.9$, and $\sigma=0.1$, and estimate the state using ADMM, FIS, and SMC approaches. For the ADMM approach, we note that the true system model is no longer log-concave, so we instead use a sparsity inducing $\ell_{1}$ regularizer, i.e. we define $\phi(\mathbf{w})=\beta\left|\left|\mathbf{w}\right|\right|_{1}$. As such, we set $\beta=15$, noting that is no longer determined by the model and must be treated as a tuning parameter. The resulting system model update is given by:

[TABLE]

This problem is known as the LASSO problem and may be efficiently solved by applying a soft threshold operation to $\mathbf{w}^{(i)}$ at each iteration [36].

Given the model mismatch, we observe that the proposed method takes longer to converge on a desirable estimate, and thus increase the maximum number of iterations to 75. For the FIS, given that there is no systematic approach to obtain an estimate with sparse variations, we again utilize a Gaussian approximation, with the noise at each step being modeled by a Gaussian distribution with zero-mean and variance $\operatorname{Var}(U_{n})=4\sigma^{2}$. The SMC method is given the benefit of using the true underlying system model when generating samples on the forward pass. However, when performing the backward pass on sample $x_{n}^{i}$ with respect to a fixed $\hat{x}_{n+1}$, we get that when $x_{n}^{i}>\hat{x}_{n+1}$, the likelihood $f_{X_{n+1},\mathbf{Y}\mid X_{n}}(\hat{x}_{n+1},\mathbf{y}\mid x_{n}^{i})=0$, causing the smoother to continually lower $\hat{x}_{n}$ for $n=N,N-1,\dots,1$ until the smoother fails (i.e. $x_{k}^{i}>\hat{x}_{k+1}$ for all $i$ for some $k\in\{1,\dots,N\}$). As such, we only utilize the forward pass particle filter. Referring to Table III for results, we note that the proposed method again outperforms the other methods in the RMSE sense. From a computational perspective, the 3X increase in iterations causes the ADMM approach to take slightly longer than the FIS, though both remain significantly more efficient than the SMC method.

IV-B Spectrotemporal Pursuit

Next we demonstrate application of the ADMM framework to the method of spectrotemporal pursuit, originally presented in [11]. Spectrotemporal pursuit formulates the problem of estimating time varying frequency coefficients as a compressive sensing problem. We define $\mathbf{Y}\in\mathbb{R}^{P\times N}$ to be a matrix version of an observed time series of length $PN$, where each column of $\mathbf{Y}$ gives a length $P$ window of the time series. Next, we define $\mathbf{X}\in\mathbb{R}^{K\times N}$ to be a matrix of frequency coefficients, with each column $\mathbf{X}_{n}\in\mathbb{R}^{K}$ representing the frequency coefficients corresponding with the time window $\mathbf{Y}_{n}\in\mathbb{R}^{P}$. By defining $\mathbf{X}$ to be real valued, it is implied that the frequency coefficients are in rectangular form, and thus a frequency resolution of $\nicefrac{{K}}{{2}}$ is achieved. Using this representation, we define the quadratic measurement model:

[TABLE]

where $\mathbf{F}_{n}\in\mathbb{R}^{P\times K}$ is an inverse Fourier matrix, i.e. $(\mathbf{F}_{n})_{p,k}\coloneqq\cos(2\pi((n-1)P+p)\frac{k-1}{K}$ and $(\mathbf{F}_{n})_{p,k+\frac{K}{2}}\coloneqq\sin(2\pi((n-1)P+p)\frac{k-1+\nicefrac{{K}}{{2}}}{K}$ for $p=1,\dots,P$ and $k=1,\dots,\nicefrac{{K}}{{2}}$. In this sense we can view the spectrotemporal estimation problem as a traditional linear measurement with Gaussian noise problem. As such, it is well defined when $P\geq K$, which is consistent with the well known fact that the number of frequency coefficients associated with a time series can not exceed the number of samples.

The method of spectrotemporal pursuit removes this constraint by introducing a sparsity inducing prior on the frequency coefficients, paralleling the approaches in compressive sensing used to estimate the coefficients underlying a system with an underdetermined set of observations. In particular, spectrotemporal pursuit imposes a group-sparsity prior on the first differences of the frequency coefficients. Letting $\mathbf{W}_{n}=\mathbf{X}_{n}-\mathbf{X}_{n-1}$ (i.e. $\mathbf{D}$ is the identity matrix), we define the system model:

[TABLE]

We can view this function as the $\ell_{1}$-norm of a vector whose entries are the $\ell_{2}$-norms of the rows of the argument. As such, $\phi(\mathbf{w})$ is small when only a small number of the rows of $\mathbf{w}$ are non-zero. Furthermore, the rows that are non-zero should have a small $\ell_{2}$-norm. Application of this function to the differences of the frequency coefficients over time ensures that throughout a given time series, most frequency coefficients do not vary, and those that do vary are varying smoothly. This time-frequency characterization is known to occur in certain biological time-series. Thus, spectrotemporal pursuit utilizes this knowledge to obtain significantly denoised spectrotemporal estimates while avoiding the time/frequency resolution trade-off without necessitating a sliding window approach. This is again reminiscent of compressive sensing, which makes strong claims regarding the recoverability of a set of coefficients with underdetermined measurements so long as the coefficients are sufficiently sparse.

The spectrotemporal pursuit solution initially proposed in [11] is an iteratively reweighted least squares (IRLS) algorithm. While the IRLS algorithm is also exact and offers convergence guarantees, it requires inversion of $N\times N$ and $K\times K$ matrices $N$ times per iteration of the algorithm. Furthermore, design of the state-covariance matrix obfuscates the problem and requires careful thought when modifying the system model.

The proposed ADMM framework yields a straightforward solution to the spectrotemporal pursuit problem. First, plugging $L$ into equation (12) yields:

[TABLE]

where $\mathbf{C}_{n}\coloneqq(\mathbf{F}_{n}^{T}\mathbf{F}_{n}+\frac{\rho}{2}\mathbf{I})^{-1}$ and $\mathbf{b}_{n}^{(i)}\coloneqq-\frac{1}{2}(\mathbf{F}_{n}^{T}\mathbf{y}_{n}+\rho\tilde{\mathbf{x}}_{n}^{(i)})$. We note that when $P<K$, $\mathbf{F}_{n}^{T}\mathbf{F}_{n}$ is rank deficient and it is our choice of $\rho$ that ensures the update is well formed. Also, it is important to note that each $\mathbf{C}_{n}$ for $n=1,\dots,N$ can be computed once at initialization, as they do not change throughout iterations.

Next, placing the group-sparsity prior in equation (14) shows that the system model update is given by a standard group-lasso problem:

[TABLE]

Furthermore, this special case with an orthonormal regressor matrix (i.e. the identity) yields a closed form solution, namely a row-wise shrinkage operator applied to $\tilde{\mathbf{w}}^{(i)}$ [24]. The shrinkage amount is proportional to the tuning parameter $\beta$, with larger $\beta$ yielding a smaller number of non-zero rows in $\mathbf{w}$.

We demonstrate the ADMM solution for spectrotemporal pursuit on a simulated example recreated from the original paper [11]. Let $\tilde{\mathbf{y}}\in\mathbb{R}^{M}$ be the vectorized version of $\mathbf{y}$ with $M=NP$ and $\mathbf{y}_{n}=[\tilde{y}_{(n-1)P+1},\tilde{y}_{(n-1)P+2},\dots,\tilde{y}_{nP}]^{T}$ for $n=1,\dots,N$. Then, we consider the signal:

[TABLE]

where $f_{0}=0.04$ Hz, $f_{1}=10$ Hz, $f_{0}=11$ Hz, and $v_{m}\sim\mathcal{N}(0,1)$ iid for $m=1,\dots,M$. Letting the sampling frequency be $f_{s}=125$ Hz and $M=7500$ gives a simulated time-series 600 seconds in duration. We note that $\mathbf{y}$ contains a sparse number of active frequency components, and the frequency components that are active are modulated over time in a smooth fashion. Additionally, the active frequency components $f_{1}$ and $f_{2}$ are chosen to be in neighboring frequencies, creating an increased difficulty when trying to distinguish their respective contributions.

The top row of Fig. 3 shows time-frequency estimates of the simulated time-series using traditional methods and spectrotemporal pursuit. First, we observe that the standard spectrogram (Fig. 3A) suffers from significant spectral leakage and is unable to clearly distinguish between the 10 Hz and 11 Hz frequency components. For the spectrotemporal pursuit estimate (Fig. 3B) we select $P<K$, meaning that the number of samples in each time window is less than the number of frequency bins. As such, we are effectively increasing the temporal resolution while still maintaining the spectral resolution without the use of overlapping windows. Because this would in general be an underdetermined problem, the group-sparsity prior is needed to ensure the problem has a unique solution. In addition to increased temporal resolution, we witness that spectrotemporal pursuit enables the contributions from $f_{1}$ and $f_{2}$ to be clearly distinguishable. Further benefits of this approach to spectrotemporal decompositions are given in detail in [11]. Here, we are proposing an algorithm that offers improvements in efficiency, modularity, and interpretability. In particular, we witness a roughly $10\times$ speedup per iteration on the same size data when using the ADMM framework rather than IRLS.

To further illustrate the modularity of the proposed framework, we next demonstrate that we can utilize an entirely different system model with a minor adjustment to a single update. Specifically, we consider a low-rank spectrotemporal decomposition (LRSD) which substitutes the nuclear norm for the group sparsity prior [37]. As such, the LRSD estimate is obtained by substituting the system model update given by (31) with:

[TABLE]

where $\lvert\lvert\cdot\rvert\rvert_{*}$ is the nuclear norm, given by the sum of the singular values of the argument. Conveniently, this update is known as the matrix lasso and yields a straightforward solution via singular value soft thresholding [38]. By making a simple adjustment to the means by which $\mathbf{w}^{(i)}$ is updated, we are able to obtain an entirely different spectrotemporal decomposition.

This point is illustrated by the bottom row of Fig. 3 where we demonstrate the LRSD on human single-channel EEG data using adhesive flexible sensors [39]. The data in question contains a 30-second recording in which the subject’s eyes are closed at the 10 second mark, at which point we would expect to see increased energy in the alpha band (10-12 Hz). The change point nature of the recording suggests that the group sparsity prior on the dynamics, which enforces smoothness across time, is ill-suited for this recording, and the traditional spectrogram (Fig. 3C) suffers significantly from noise. By not explicitly enforcing smoothness in time, the low-rank enforcing nuclear norm prior (Fig. 3D) accommodates the change point and is able to significantly suppress activity outside of the alpha band. Similarly to the spectrotemporal pursuit example, we are able to set $P<K$ and achieve equivalent temporal resolution to the spectrogram without utilizing overlapping windows or sacrificing spectral resolution.

Remark 2**.**

Comparisons with other methods are intentionally omitted in this section given that there is no systematic application to these non-Markov problem formulation. While the original problem proposed in equation (7) does not lend itself to an obvious solution for the discussed non-Markov models, the consensus ADMM formulation given by (11) may be solved in a straightforward manner. In particular, we note that the EKF and UKF have no clear extensions for non-Markov scenarios and the and the use of sampling based methods for such models would require drawing samples of group-sparse or low-rank matrices.

V Discussion

We have presented a unified framework for solving a broad class of dynamic modeling problems. The proposed method can be applied to systems with non-linear measurements and/or non-Markov dynamics. As demonstrated on two applications, our framework can be applied in a straightforward manner to acquire efficient solutions to problems that may otherwise require complex or approximate solutions. Furthermore, we have shown that this algorithm will converge on the true MAP estimate of the latent signal in the limit of large iterations. With this provably accurate algorithm comes a mathematical justification for an intuitive approach to dynamic time-series estimation, namely iteratively computing estimates based on the measurement model and system model and then averaging them in the appropriate sense.

There are a number of extensions to this framework still to be explored. The most glaring shortcomings are the inability to conduct the estimation procedure causally and the necessity to know model parameters a priori. Regarding the former, we note the use of homotopy schemes for causal estimation that gradually incorporate new observations into the solution [40, 41]. Additionally, there has been recent research investigating algorithms for performing ADMM in an online fashion [42, 43] that could potentially be leveraged by our framework. To address the latter, expectation-maximization (EM) techniques can be built into the ADMM iterations in order to estimate model parameters jointly with the desired latent time-series. In that regard, the E-step, which requires sampling from the posterior distribution, is typically the bottleneck. To address that, Langevin based methods and stochastic gradient descent methods can be used to efficiently sample from the posterior distribution [44]. Identifying sufficient conditions on mixing times for generating approximately i.i.d. posterior samples for the M-step could be the subject of future in-depth work. We note that while there exist sample based methods for estimating model parameters [19, Sec. IV], these methods can be computationally prohibitive as witnessed in Table III.

Lastly, we note that there is considerable interest in state-space estimation where the observations or system are subject to noise from heavy-tailed distributions such as the Student’s $t$ or Cauchy distributions [45, 46], which are not log-concave. Recent literature has shown that ADMM can be shown to converge under even milder conditions than those assumed by Theorem 1 [47, 48]. Given that both the Student’s $t$ and Cauchy distributions are log-quasi-concave, continuous, and possess a single local maximum, we could reasonably expect convergence of our framework to the MAP estimate in such a scenario. This topic provides interesting opportunities for future experimental and theoretical work.

Appendix A Derivation of Scaled Form

We will demonstrate the derivation of the scaled form of the measurement model update only, noting that the derivation for the other updates follows almost identical steps. Consider the original measurement model update:

[TABLE]

For ease of notation, the superscript $(i)$ is omitted for the remainder of this appendix. Using the definition of the inner product and Frobenius norm, we can break up the second and third terms across into sums and simplify as follows:

[TABLE]

Defining $\tilde{\mathbf{x}}_{n}=\mathbf{z}_{n}-\frac{\boldsymbol{\lambda}_{n}}{\rho}$ as in Section III-A, we note that $\tilde{\mathbf{x}}_{n}$ does not depend on $\mathbf{x}$, enabling us to complete the square and simplify as follows:

[TABLE]

as was to be shown.

Appendix B Proof of Theorem 1

Consider the problem in its original form:

[TABLE]

The goal is to show that there is an equivalent two-block ADMM problem whose updates match those given by (11). To do so we define the variable $\mathbf{Q}\coloneqq[\mathbf{X}^{T},\mathbf{W}^{T}]^{T}\in\mathbb{R}^{2K\times N}$ ($\mathbf{X},\mathbf{W}\in\mathbb{R}^{K\times N}$) and the function $g(\mathbf{Q})\coloneqq L(\mathbf{y}\mid\mathbf{X})+\beta\phi(\mathbf{W})$. Next, we define $\mathbf{Z}\coloneqq[\mathbf{Z}_{\mathbf{X}}^{T},\mathbf{Z}_{\mathbf{W}}^{T}]^{T}\in\mathbb{R}^{2K\times N}$ and the function:

[TABLE]

Using these newly defined terms, we can write (37) equivalently as:

[TABLE]

Note that if $\mathbf{q}=[\mathbf{x}^{T},\mathbf{w}^{T}]^{T}$ is such that $\mathbf{w}\neq\mathcal{A}(\mathbf{x})$ (the constraints in (37) are not satisfied) and $\mathbf{z}$ is such that $\mathbf{q}-\mathbf{z}=0$ (the constraints (39) are satisfied), then $h(\mathbf{z})=\infty$ and $(\mathbf{q},\mathbf{z})$ are not the minimizers of (39). To solve this problem with ADMM, we first find augmented Lagrangian:

[TABLE]

with Lagrange multiplier $\boldsymbol{\gamma}=[\boldsymbol{\lambda}^{T},\boldsymbol{\alpha}^{T}]^{T}\in\mathbb{R}^{2K\times N}$ ($\boldsymbol{\lambda},\boldsymbol{\alpha}\in\mathbb{R}^{K\times N}$). As a result, we get the following update equations:

[TABLE]

Next we show that the update equations given by (41) are equivalent to those given by (11).

First, consider the $\mathbf{q}$ update:

[TABLE]

where $\mathbf{x}^{(i+1)}$ and $\mathbf{w}^{(i+1)}$ are given by:

[TABLE]

and can be found independently of each other.

Next, consider the $\mathbf{z}$ update:

[TABLE]

where $\mathbf{z}_{\mathbf{x}}^{(i+1)}$ is given as the solution to (15), i.e. the consensus update for our target problem, and the second to last equality follows from the fact that $h(\mathbf{z})$ is infinite if $\mathbf{z}_{\mathbf{w}}\neq\mathcal{A}(\mathbf{z}_{\mathbf{x}})$, so we can treat the problem as a single variable optimization problem.

Next we can substitute these results into the equations for the $\mathbf{q}$ update to obtain:

[TABLE]

which is the (unscaled) update equation (14) for $\mathbf{w}$ in the original formulation, where $\mathbf{z}_{\mathbf{x}}^{(i)}$ in this formulation corresponds with $\mathbf{z}^{(i)}$ in the original formulation. The $\mathbf{x}$ portion of the $\mathbf{q}$ remains unchanged from (43), which is equivalent to the unscaled update equation (12) for $\mathbf{x}$ in the original formulation.

Next, we can decompose the matrix multiplication in the same way as above to show that:

[TABLE]

where $\boldsymbol{\lambda}^{(i+1)}$ and $\boldsymbol{\alpha}^{(i+1)}$ are given by the original updates in (11).

Thus, we have shown that directly solving (39) using ADMM yields the proposed updates detailed in the body of the paper. As such, we will show that the ADMM solution to (39) is convergent. By assumption, $L$ and $\phi$ are closed, proper, and convex, and hence, so is their sum $g$. To show that $h$ is convex, we note that this is an indicator function on the set $H\coloneqq\{(\mathbf{z}_{\mathbf{X}},\mathbf{z}_{\mathbf{W}}):\mathcal{A}(\mathbf{z}_{\mathbf{X}})=\mathbf{z}_{\mathbf{W}}\}\subset\mathbb{R}^{2}$, thus $h$ is convex if and only if $H$ is convex [49, Ch. 2]. Suppose $\mathbf{z}^{1}=[{\mathbf{z}_{\mathbf{x}}^{1}}^{T},{\mathbf{z}_{\mathbf{w}}^{1}}^{T}]^{T}$ and $\mathbf{z}^{2}=[{\mathbf{z}_{\mathbf{x}}^{2}}^{T},{\mathbf{z}_{\mathbf{w}}^{2}}^{T}]^{T}$ are such that $\mathcal{A}(\mathbf{z}_{\mathbf{x}}^{1})=\mathbf{z}_{\mathbf{w}}^{1}$ and $\mathcal{A}(\mathbf{z}_{\mathbf{x}}^{2})=\mathbf{z}_{\mathbf{w}}^{2}$, i.e. $\mathbf{z}^{1},\mathbf{z}^{2}\in H$. Then, if we take a convex combination $\mathbf{z}^{\alpha}\coloneqq\alpha\mathbf{z}^{1}+(1-\alpha)\mathbf{z}^{2}$ for $\alpha\in[0,1]$, we get:

[TABLE]

Thus, we see that $\mathbf{z}^{1},\mathbf{z}^{2}\in H\implies\mathbf{z}^{\alpha}\in H$, i.e. $H$, and therefore $h$, are convex. It then follows from Section 3.2.1 of [15] that the ADMM solution for (39) is convergent, as was to be shown. $\blacksquare$

Appendix C State-Space Model of Learning Updates

We begin by deriving expressions for the negative log-likelihoods for each of the observations:

[TABLE]

[TABLE]

[TABLE]

These expressions can be plugged into equation (25) to obtain the measurement model update equation, which can in turn be solved using Newton’s method.

Next, the system model update can be solved in closed form:

[TABLE]

where $\tilde{w}_{n}^{(i)}\coloneqq z_{n}^{(i)}-\kappa z_{n-1}^{(i)}-\nicefrac{{\alpha_{n}^{(i)}}}{{\rho}}$. Thus, we can solve for each $w_{n}$ separately:

[TABLE]

Finally, given its relatively low dimensionality, we can efficiently solve the consensus update in closed form by posing it as a least squares problem. First, we note that $\mathcal{A}(\mathbf{z})=\mathbf{G}\mathbf{z}$ when we define:

[TABLE]

with $\mathbf{G}\in\mathbb{R}^{N\times N}$. Thus, we have:

[TABLE]

Taking the gradient of the RHS and setting to zero yields:

[TABLE]

Given that $\mathbf{G}$ is known a-priori, we can find $(\mathbf{I}+\mathbf{G}^{T}\mathbf{G})^{-1}$ once and each consensus update becomes a matrix multiplication problem.

Appendix D Convexity State-Space Model of Learning Negative Log-Likelihood

Given that $L$ is the sum of the negative log-likelihoods for each of the observation modalities as in (25), it is sufficient to show that they are each convex in $x_{n}$, which is made easier by use of the simplifications derived in Appendix C. Noting that addition of a constant does not affect convexity, we can assess the final simplification provided in each case. As such, we see that $L_{B_{n}}(b_{n}\mid x_{n})$ is the sum of a term that is linear in $x_{n}$ and a special case of the log sum exponential (LSE) function with an added auxiliary variable constrained to equal zero (giving $e^{0}=1$). Given the convexity of LSE, its sum with a linear term is also convex, and thus $L_{B_{n}}(b_{n}\mid x_{n})$ is convex. Next, $L_{R_{n}}(r_{n}\mid x_{n})$ is quadratic in $x_{n}$ and thus convex. Finally, $L_{\mathbf{S}_{n}}(\mathbf{s}_{n}\mid x_{n})$ is the sum of a term that is linear in $x_{n}$ and a term that is exponential in $x_{n}$, both of which are convex. As a result, $L_{B_{n}}(b_{n}\mid x_{n})$, $L_{R_{n}}(r_{n}\mid x_{n})$, and $L_{\mathbf{S}_{n}}(\mathbf{s}_{n}\mid x_{n})$ are all convex in $x_{n}$ for any $(b_{n},r_{n},\mathbf{s}_{n})\in\{0,1\}\times\mathbb{R}\times\{0,1\}^{J}$, and thus so is their sum $L$.

Acknowledgment

The authors would like to thank M. Wagner for his contributions to the development of the LRSD algorithm and the reviewers for their insightful comments and suggestions for improvements.