text:
\begin{document}
\maketitle
\begin{abstract}
    This paper studies the effect of parametric mismatch  in  minimum mean
    square error (MMSE) estimation. In particular, we consider the problem of estimating the input signal
    from the output of an additive white Gaussian channel whose gain is fixed, but unknown.
    The input distribution is known, and the estimation process consists of two algorithms. First, a channel
    estimator blindly estimates the channel gain using past observations. Second, a mismatched
    MMSE estimator, optimized for the estimated channel gain, estimates the input signal.
    We analyze the {\em regret}, i.e., the additional mean square error, that is raised in this process.
    We derive upper-bounds on both absolute and relative regrets. Bounds are expressed in terms of
    the Fisher information. We also study regret for unbiased, efficient channel estimators, and derive
    a simple trade-off between Fisher information and relative regret.
    This trade-off shows that the product of a certain function of
    relative regret and Fisher information equals  the signal-to-noise ratio, independent of
    the input distribution.
    The trade-off  relation implies
    that higher Fisher information results to smaller~expected~relative~regret.
\end{abstract}

    \section{Introduction}
        Consider an application that you are given the output of a system, and you seek
        to recover the input of the system. You know that the system is noisy, e.g., it adds white Gaussian noise
        to the output. You know the
        distribution of the input, but you do not know the system parameters. Problems
        of this sort arise in different applications in signal processing and communication systems.
        Some examples include blind deconvolution \cite{Donoho:Deconvolution}, dereverberation \cite{Majid:ICASSP:10},
        denoising \cite{Weissman:DUDE:Uncetainty}, and mismatch decoding \cite{Mismatch:Lapidoth}.
        These applications differ in their fundamental models, fidelity criteria,
        and methodologies.  However, they have one thing in common: they all suffer from parametric mismatch in
        recovering the input signals.

        The motivation of this work is blind deconvolution and dereverberation applications. Linear time-invariant
        channels serve as common models in these applications. As the input signal passes through these channels, it convolves
        with the unknown finite-impulse response (FIR) of the channel, and it adds with additive white Gaussian noise (of known variance).
        Recovering the input signals from the noisy output could be impossible even with perfect knowledge about the
        channel response. This is out of our scope. Instead, we aim to study the penalty and performance degradation that
        is specifically caused by the lack of knowledge about the channel response.

        We benchmark performance
        against that of perfect channel knowledge scenario. We are concerned about
        issues such as required sample complexity or training in channel estimation to bring performance of input estimation
        within a desired range. As a counterpart problem in communication systems, one may think of block fading channels
        and the trade-off between accuracy of channel estimation and performance of decoding~\cite{Hassibi:Training}. Note
        that channel estimation in our case is blind as we have no control of the source.

        As a first step to address these problems, in this work, we focus on the most basic system in which
        the unknown channel is just a single gain.  We expect that the results and intuitions
        of this work will shed lights on the analysis of generic FIR channels.\footnote{Analogous to the case between the analysis of flat-fading and the analysis of frequency-selective channels.}
        In treating the problem, we consider an estimation process that consists of two algorithms. First, a {\em channel
        estimator} blindly estimates the channel gain using past output observations. Second, a mismatched
        {\em minimum mean square error} (MMSE) estimator, optimized for the estimated channel gain, estimates the input signal.
        Figure~\ref{fig:system} illustrates the building blocks of this process. Due to estimation error in
        channel estimation, the MMSE estimator that is used in recovering the input signal results in a
        mean square error that is larger than that of the ideal MMSE estimator.
        We call this additional error as {\em regret}, and
        we derive novel upper-bounds on both absolute and relative regrets. The bounds are simple and
        demonstrate interesting connections to  the Fisher information. To this end, one might attempt
        to exploit the results of \cite{Sergio:Mismatched} and
        \cite{Guo:RelativeEntropy:Score} to derive other alternative bounds.

        We also quantify  regret for unbiased, {\em efficient} channel estimators. Since these
        estimators achieve Cramer-Rao bound, they result in  a simple trade-off relation between Fisher information and relative regret.
        This trade-off relation expresses that the product of a certain function of
        relative regret and Fisher information is equivalent to the signal-to-noise ratio, independent of the input distribution.
        Trade-off suggests that higher Fisher information results to smaller expected relative regret.
        Although, intuitively, this may seem expected, simplicity of the trade-off relation makes it worthwhile.

    \section{Setup}
        Consider a linear dynamic system
        \begin{eqnarray}
            Y_n = aX_{n}+ V_n
            \label{eqn:dynamic}
        \end{eqnarray}
        in which $\{V_n\}$ is an independent, identically, distributed (i.i.d.) Gaussian noise
        such that $V_n\sim\mcal N(0,\s^2_v)$. The input $X_{n}$
        is an i.i.d. process whose distribution is known to be $\probdist{}{X}$.
        Parameter $a \in \mbb R^+$ is a fixed,  unknown channel gain. It results to
        a derived parametric family of probability measures $\probdist{a}{X,Y}$,
        the joint distribution of $X$ and $Y$, governing the system dynamic \eqref{eqn:dynamic}.
        The objective is to observe a realization of the output process
        \[Y^n=(Y_1,Y_2, \cdots, Y_n)\] and estimate the realization of the
        underlying input process, i.e.,
        \[X^n=(X_1, X_2, \cdots, X_n).\]

        Let $\mcal X =\mbb R$ and $\mcal Y = \mbb R$ denote the input and output spaces, respectively.
        We consider memoryless input estimators, e.g., $\phi \colon \mcal Y \to \mcal X$ where
        $\phi(Y_n)$ is an estimate for $X_n$.  The {\em mean square error} (MSE) for  $\f$ is defined
        \begin{eqnarray}
             \Exp{(X-\phi(Y))^2} = \int (x-\phi(y))^2 d \probdistn{a}.
            \label{eqn:mse}
        \end{eqnarray}
        In Eq. \eqref{eqn:mse} and henceforth we follow the convention that
        unsubscribed expectations are measured according to $\probdist{a}{X,Y}$.
        Moreover, we use concise notations like $\probdistn{a}=\probdist{a}{X,Y}$
        and $\probdistn{a|y}=\probdist{a}{X|Y=y}$ to denote joint and conditional
        distributions, respectively.

        One seeks to find an estimator that minimizes MSE~\eqref{eqn:mse}. The main challenge, however, is that $a$ and
        $\probdistn{a}$ are unknown. If we had
        oracle knowledge about $a$, the MMSE estimator for $X$
        is defined
        \begin{align}
             \phi_a(y) = \Exp{X|Y=y}.
             \label{eqn:optimal:mmse:estimator}
        \end{align}
        for an observation $Y=y$. Any other estimator $\f$ results to additional
        error that we call it {\em regret}. The motivation for this name is that
        it measures degradation on performance, an impact caused by imprecise knowledge about $a$.

        In this paper, we assess regret for a special class of mismatched estimators. Namely,
        we consider an estimation process that is depicted in Figure \ref{fig:system}. A channel
        estimation works in parallel with an MMSE input estimation as follows. At time instance $n$,
        a channel estimator  finds an  estimate $\hat{a}=\hat{a}_n$ of $a$ using the observed values $Y^{n-1}$.
        Then, it uses the optimal estimator of $\probdist{\hat{a}}{X, Y}$
        to compute
        \begin{align}
            \phi_{\hat{a}}(y_n) = \Expi{\hat{a}}{X_n|Y=y_n}
            \label{eqn:mismatch:estimator}
        \end{align}
        as an estimate for $X_n$. Function  $\phi_{\hat{a}}$ is a mismatch MMSE
        estimator that causes regret when used in place of $\f_a$. In the following sections, we
        study two types of regret: {\em absolute regret} and {\em relative regret}.

            \begin{figure}[t]
                \begin{center}
                    \psfrag{X}[Bc][B1][\figfontsmall][0]{$X_n$}
                    \psfrag{Y}[Bc][B1][\figfontsmall][0]{$Y_n$}
                    \psfrag{Yn}[Bc][B1][\figfontsmall][0]{$Y^{n-1}$}
                    \psfrag{V}[Bc][B1][\figfontsmall][0]{$V_n$}
                    \psfrag{a}[Bc][B1][\figfontsmall][0]{$a$}
                    \psfrag{aX}[Bc][B1][\figfontsmall][0]{$aX_n$}
                    \psfrag{ah}[Bc][B1][\figfontsmall][0]{$\hat{a}$}
                    \psfrag{Xh}[Bc][B1][\figfontsmall][0]{$\hat{X}_n$}
                    \psfrag{C}[Bc][B1][\figfontsmall][0]{channel}
                    \psfrag{E}[Bc][B1][\figfontsmall][0]{estimation}
                    \psfrag{I}[Bc][B1][\figfontsmall][0]{MMSE}
                    \psfrag{e}[Bc][B1][\figfontsmall][0]{estimation process}
                    \includegraphics[width=.95\columnwidth]{system.eps}
                \end{center}
                \caption{Figure depicts the building blocks of the system setup.
                The estimation process consists of two individual algorithms: 1) a channel estimation algorithm
                that blindly estimates $\hat{a}$ as an estimate for $a$, 2) a mismatch MMSE estimation algorithm, optimized for $\hat{a}$, that recovers $X_n$.}
                \label{fig:system}
            \end{figure}

    \section{Absolute Regret}
        \subsection{Deviation Analysis}
        The absolute regret corresponding to $\phi_{\hat{a}}$ is
        \begin{align}
            \regret{\hat{a},a} = \Exp{(X-\phi_{\hat{a}}(Y))^2} - \Exp{(X-\phi_a(Y))^2}.
        \end{align}
        Application of orthogonality principle results to
        \begin{align}
            \regret{\hat{a},a} = \Exp{(\phi_{\hat{a}}(Y)  - \phi_a(Y))^2}.
            \label{eqn:regret}
        \end{align}
        Eq. \eqref{eqn:regret} quantifies the absolute {\em regret} of using $\phi_{\hat{a}}$ instead of $\phi_a$. The following lemma
        states and proves an upper-bound on  \eqref{eqn:regret}.
        \begin{lemma}
            For every $\hat{a}$, the following holds true
            \begin{align}
                 \nonumber \regret{\hat{a},a}  \leq (\hat{a}-a)^2 & \Exp{\left(6 \s^2_x + 8 \frac{Y^2}{{a}^2}\right)  \fisher{X;a||Y}}\\&+o(\hat{a}-a)^2
                \label{eqn:lemma:upperbound:regret}
            \end{align}
            in which the expectation is with respect to $Y$,

            and
            \begin{eqnarray}
                 \fisher{X;a||Y}  \triangleq \Expi{}{\big(\nabla \ln \dfi{a}{X|Y} \big)^2 |Y}
                 \label{eqn:lemma:fisher:random}
            \end{eqnarray}
            is the Fisher information of $X$ relative to $a$, conditioned on $Y$.
            Here, $\dfi{a}{X|Y}$ is  the density of $\probdistn{{a}|Y}$.
            \label{lemma:upperbound:regret}
        \end{lemma}
        \begin{proof}
            Refer to Appendix \ref{app:upperbound:regret}.
        \end{proof}
        Lemma \ref{lemma:upperbound:regret} describes a bound \eqref{eqn:lemma:upperbound:regret} that
        comprises two  multiplicative terms. The first term $(\hat{a}-a)^2$ measures the
        channel estimation error. The second term is the weighted average of conditional
        Fisher information. Intuitively, this term  measures the amount of information
        that an observable random variable $X$ carries about unknown parameter $a$ conditioned on $Y$, assigning
        more weight to larger values of $Y$.

        \begin{corollary}
            If  $|\hat{a}-a| << 1$, and

            if $\fisher{X;a||Y}$ and $Y^2$ are uncorrelated, we obtain the simple
            bound
            \begin{align}
                \regret{\hat{a},a} \leq (\hat{a}-a)^2 (14 \s^2_x + 8 \frac{\s^2_v}{{a}^2}) \fisher{X;a|Y}
                \label{eqn:upperbound:regret:small:deviation:simple}
            \end{align}
            in which
                \begin{eqnarray}
                    \fisher{X;a|Y} \triangleq \Expi{}{\big(\nabla \ln \dfi{a}{X|Y} \big)^2}
                    \label{eqn:fisher:average}
                \end{eqnarray}
            is the average of ${\fisher{X;a||Y}}$ with respect to $Y$.\footnote{Lookout for the subtle notational difference between ${\fisher{X;a||Y}}$, a random variable,
            and $\fisher{X;a|Y}$, a scalar.}
        \label{cor:regret:approximation}
        \end{corollary}

        \subsection{Efficient Channel Estimation}
        Neither Eq. \eqref{eqn:lemma:upperbound:regret} nor Eq. \eqref{eqn:upperbound:regret:small:deviation:simple} depend on the channel estimation
        algorithm that estimates $a$. They simply relate small deviation between $\hat{a}$ and $a$ to
        absolute regret in estimating $X$. To incorporate the effect of channel estimation algorithm, we proceed as follows.

        As mentioned earlier, at time $n$, $\hat{a}$ is obtained through
        observation of $Y^{n-1}=(Y_i)_{i=1}^{n-1}$. In formal terms,
        \[\hat{a} = A_n(Y^{n-1})\]
        where $A=(A_1, A_2, \cdots)$ is a channel estimation algorithm in which $A_n\colon \mcal Y^{n-1} \to \mbb R^+$.
        \begin{lemma}
            Let $\mcal A$ denote the class of all unbiased channel estimation algorithms.
            If $\mcal A$ contains an {\em efficient estimator} \cite[p. 92]{Borovkov:MathematicalStatistics},
            the following holds true
            \begin{align}
                \nonumber \inf_{A\in \mcal A} \Exp{\regret{A_n(Y^{n-1}),  a}} \leq \qquad \qquad \qquad\\ \frac{1}{n-1}  \frac{\Exp{\left(6 \s^2_x + 8 \frac{Y^2}{{a}^2}\right)  \fisher{X;a||Y}}}{\fisher{Y;a}}
                \label{eqn:esterror:aregret:Cramer-Rao}
            \end{align}
            for sufficiently large values of $n$.\footnote{The expectation in the LHS is with respect to $Y^{n-1}$.}
            \label{lemma:regret:Cramer-Rao}
        \end{lemma}
        \begin{proof}
            Refer to Appendix \ref{app:prf:regret:Cramer-Rao}.
        \end{proof}

    \section{Relative Regret}
        \subsection{Deviation Analysis}
        Let
        \begin{align}
            \rregret{\hat{a},a} = \Exp{\frac{(\phi_{\hat{a}}(Y)  - \phi_a(Y))^2}{\Expi{\hat{a}}{X^2|Y} + \Expi{{a}}{X^2|Y}}}
            \label{eqn:relative:regret}
        \end{align}
        denote the relative regret. The following lemma states and proves a simple upper-bound
        on  $ \rregret{\hat{a},a}$.
        \begin{lemma}
            For every $\hat{a}$, we have
            \begin{align}
                \rregret{\hat{a},a}  \leq & {(\hat{a}-a)^2} \;  \fisher{X;a|Y} + o(\hat{a}-a)^2
                \label{eqn:lemma:relative:regret}
            \end{align}
            where $\fisher{X;a|Y}$ is defined as Eq. \eqref{eqn:fisher:average}
            and denotes the conditional Fisher information of $X$ relative to $a$.
            \label{lemma:upperbound:relative:regret}
        \end{lemma}
        \begin{proof}
            See  Appendix \ref{app:prf:lemma:upperbound:relative:regret}.
        \end{proof}
        Eq. \eqref{eqn:lemma:relative:regret}
        results to a simple upper-bound on the relative regret for small deviations between $\hat{a}$
        and $a$.

        \subsection{Efficient Channel Estimation}

        Similar to the case for absolute regret, we now state the following
        result.
        \begin{lemma}
            Let $\mcal A$ denote the class of all unbiased estimation algorithms.
            If $\mcal A$ contains an {\em efficient estimator},
            the following holds true
            \begin{eqnarray}
                \inf_{A\in \mcal A}\Exp{\rregret{A_n(Y^{n-1}), a}} \leq  \frac{1}{n-1} \; \frac{ \fisher{X;a|Y}}{\fisher{Y;a}}
                \label{eqn:rregret:cramer-rao}
            \end{eqnarray}
            for sufficiently large values of $n$.
            \label{lemma:rregret:Cramer-Rao}
        \end{lemma}
        \begin{proof}
            The proof of this lemma is essentially the same as the proof of~Lemma~\ref{lemma:regret:Cramer-Rao}.
        \end{proof}
        Lemma \ref{lemma:rregret:Cramer-Rao} describes a bound on the expected relative regret, should
        an efficient estimator be used. This bound determines the smallest upper-bound on average relative regret, when
        sufficiently good unbiased channel estimators are used.

        \subsection{Regret Scalar}
        The constant value in the RHS of Eq. \eqref{eqn:rregret:cramer-rao} worths attention.
        It does not change with respect to $n$, and as $n \to \infty$, it becomes the sole
        scalar that determines the level of relative regret.
        We define this quantity as the {\em regret scalar} and denote it by
        \begin{eqnarray}
            \r(a)=\frac{\fisher{X;a|Y}}{\fisher{Y;a}}.
            \label{eqn:regret coefficient:fisher}
        \end{eqnarray}

        \begin{lemma}
            For every zero-mean input distribution $\probdist{}{X}$, the following trade-off holds true between
            regret scalar and output Fisher information
            \begin{eqnarray}
                (\r(a)+1) {\fisher{Y;a}} = \frac{\s^2_x}{\s^2_v}.
                \label{eqn:trade-off}
            \end{eqnarray}
            \label{lemma:tradeoff}
        \end{lemma}
        \begin{proof}
            See Appendix \ref{prf:lemma:tradeoff}.
        \end{proof}
        In Eq. \eqref{eqn:trade-off}, the RHS is the signal-to-noise ratio that is independent of $a$. Thus,
        Eq.~\eqref{eqn:trade-off} presents a simple product trade-off
        relationship between $\r(a)$ and $\fisher{Y;a}$. It suggest that the higher the Fisher information, the smaller the regret scalar,
        and vice-versa. The following example explicates this trade-off.

    \begin{example}[Gaussian Input]
        Assume $X_n\sim\mcal N(0,\s^2_x)$ and
        $V_n\sim\mcal N(0,\s^2_v)$ are  i.i.d. implying that
        $Y_n\sim\mcal N(0,a^2\s^2_x+\s^2_v)$ and $Y_n|x_n \sim \mcal N(a x_n, \s^2_v )$.
        With perfect knowledge of $a$, the ideal estimator for $X$ given $Y=y$ is
        \begin{align}
            \phi_a(y) = \frac{a\s^2_x}{a^2\s^2_x+\s^2_v} y.
        \end{align}
        The MMSE error resulting from this estimator is
        \begin{align}
            {\Expi{}{(X-\phi_a(Y))^2}}=\frac{\s^2_x \s^2_v}{a^2 \s^2_x +\s^2_v}.
        \end{align}
        A mismatch estimator for $\hat{a}$ is
        \begin{align}
            \phi_{\hat{a}}(y) = \frac{\hat{a}\s^2_x}{\hat{a}^2\s^2_x+\s^2_v} y.
        \end{align}

        We have
        \begin{eqnarray}
            \fisher{Y;a|X} = \frac{\s^2_x}{\s^2_v}
        \end{eqnarray}
         and
        \begin{align}
        \fisher{Y;a} =  \frac{2a^2\s^4_x}{(a^2\s^2_x+\s^2_v)^2}.
        \end{align}

        Thus,
        \begin{eqnarray}
            \r(a) = \frac{1}{2} \left(\frac{a^2\s^2_x}{\s^2_v} + \frac{\s^2_v}{a^2\s^2_x}\right)
        \end{eqnarray}
        Figure~\ref{fig:regret} depicts the behavior of $\r(a)$ and $\fisher{Y;a}$ with respect to
        $a$. The ${\rm SNR}=\frac{\s^2_x}{\s^2_v}=10 \text{ dB}$ and at $a=.35$, the minimum regret
        scalar coincides with maximum Fisher information.
            \begin{figure}[t]
                \begin{center}
                    \psfrag{S}[Bc][B1][\figfontsmall][0]{${\rm SNR} = \frac{\s^2_x}{\s^2_v}$}
                    \psfrag{Min}[Bc][B1][\figfontsmall][0]{Minimum Regret}
                    \psfrag{Max}[Bc][B1][\figfontsmall][0]{Maximum Fisher}
                    \psfrag{a}[Bc][B1][\figfontlarge][0]{$a$}
                    \includegraphics[width=.9\columnwidth]{regret.eps}
                \end{center}
                \caption{Figure illustrates the multiplicative trade-off between regret scalar and Fisher information.
                Smaller Fisher information results to larger regret scalar and vice versa. The ${\rm SNR}=10 \text{ dB}$ and the minimum
                regret scalar is coincident with maximum Fisher information.}
                \label{fig:regret}
            \end{figure}

    \end{example}

    \section{Recap and Conclusion}
        We considered the problem of estimating the input signal  from the output of an additive white Gaussian
        noise  channel subject to parametric uncertainty. Namely, the channel gain is fixed, but unknown.
        In treating the problem, we considered an estimation process that consists of two algorithms: a blind channel
        estimator and a mismatched MMSE estimator to estimate the input.
        We studied the  regret that is raised as a result of mismatch estimation.
        Simple upper-bounds on both absolute and relative regrets were presented. These bounds
        provide useful tools in assessing deviation in estimating the
        input when there exists a small deviation in channel gain estimation.
        The bounds are simple and expressed in terms of the Fisher information. This makes them
        more intuitive and could potentially bridge to other known results in the literature.

        We also quantified regret for unbiased, { efficient} channel estimators. Using Caramer-Rao bound,
        we derived a simple trade-off between Fisher information and relative regret.
        This trade-off expresses that the  product of a certain function of relative
        regret and the Fisher information is equivalent to the signal-to-noise ratio, independent of the input distribution.
        The trade-off suggests that the higher the Fisher information, the smaller the expected relative regret.

        This work is our initial attempt to shed light on information-theoretic
        limits of blind deconvolution and dereverberation systems. We are currently working on generalization of
        these results to these applications.

\appendix

    \section{Proofs}

    \subsection{Proof of Lemma \ref{lemma:upperbound:regret}}
    \label{app:upperbound:regret}
        To derive an upperbound on absolute regret, we first state
        and prove the following results.
        \begin{proposition}
            For every $\hat{a}$ and $y \in \mcal Y$, we have
            \begin{align}
                \nonumber (\phi_{\hat{a}}(y) & - \phi_a(y))^2  \leq \\ & 2 (\Expi{\hat{a}}{X^2|y} + \Expi{{a}}{X^2|y})  \kullback{\probdistn{\hat{a}|y}}{\probdistn{{a}|y}}.
                \label{eqn:lemma:pointwise:regret}
            \end{align}
            \label{lemma:upperbound:pointwise:regret}
        \end{proposition}
        \begin{proof}
            By definition, we have
            \begin{align*}
                \left(\phi_{\hat{a}}(y)- \phi_a(y)\right)^2  &=  \left(\int x \big( \frac{d P_{\hat{a}|y} }{ d Q} - \frac{d P_{{a}|y} }{ d Q} \big) { d Q}  \right)^2

            \end{align*}
            for every probability measure $Q$ such that $P_{a|y} \ll Q$ and $P_{\hat{a}|y} \ll Q$.
            By Cauchy Schwartz inequality, we obtain
            \begin{align}
            \nonumber (\phi_{\hat{a}}(y)- \phi_a(y))^2  & \leq   \int x^2 \left(\sqrt{\frac{d P_{\hat{a}|y} }{ d Q}} + \sqrt{\frac{d P_{{a}|y} }{ d Q}} \right)^2 dQ \\
                        & .  \int  \left(\sqrt{\frac{d P_{\hat{a}|y} }{ d Q}} - \sqrt{\frac{d P_{{a}|y} }{ d Q}} \right)^2 dQ
                \label{eqn:cauchy}
            \end{align}
            By inequality $(a+b)^2 \leq 2(a^2+b^2)$, one can show that the first term in the RHS of the above inequality is smaller than or equal to
            \[2 (\Expi{a}{X^2|y} + \Expi{\hat{a}}{X^2|y}).\]
            The second term in the RHS of inequality \eqref{eqn:cauchy}  is known as {\em Kakutani-Hellinger} distance  between $\probdist{{a}}{X|y}$ and $\probdist{\hat{a}}{X|y}$,
            denoted by \cite[p. 363]{Shiryaev:Probability}
            \begin{align*}
                r^2(\probdistn{\hat{a}|y}, \probdistn{{a}|y})=  \frac{1}{2}\int  \left(\sqrt{\frac{d P_{\hat{ a}|y} }{ d Q}} - \sqrt{\frac{d P_{{a}|y} }{ d Q}} \right)^2 dQ.

            \end{align*}
            Moreover, we know of the following inequality between Kakutani-Hellinger distance and Kullback-Leibler distance  \cite[p. 369]{Shiryaev:Probability}
            \begin{align*}
                2 r^2(\probdistn{\hat{a}|y}, \probdistn{{a}|y}) \leq \kullback{\probdistn{\hat{a}|y}}{\probdistn{{a}|y}}.

            \end{align*}
            Substituting in \eqref{eqn:cauchy}, we obtain Eq. \eqref{eqn:lemma:pointwise:regret}.
        \end{proof}

        \begin{proposition}
            For every ${a}$ and $y \in \mcal Y $, the following inequality holds true
            \begin{align}
                \Expi{{a}}{X^2|y}
                                    &\leq 3 \s^2_x + 4 \frac{y^2}{{a}^2}.
                \label{eqn:lemma:conditional:second:moment}
           \end{align}
            \label{lemma:upperbound:conditional:second:moment}
        \end{proposition}
        \begin{proof}
        Let $\dfi{{a}}{y|x}$ and $\dfi{{a}}{y}$ denote the conditional and marginal densities for
        $\probdist{{a}}{X,Y}$. Then,
        \begin{align}
            \nonumber \Expi{{a}}{X^2|y} &= \int x^2 \frac{\dfi{{a}}{y|x} }{\dfi{{a}}{y}} \dfi{}{x}dx\\
            \nonumber     & =\int_{x: \dfi{{a}}{y|x}\leq {\dfi{{a}}{y}}} + \int_{x: \dfi{{a}}{y|x}> {\dfi{{a}}{y}}}\\
                                &\leq \Expi{}{X^2} + \int_{x: \dfi{{a}}{y|x}> {\dfi{{a}}{y}}}
            \label{eqn:expectation:break}
        \end{align}
        To simplify the second term, we substitute $x^2$ by the inequality that is derived as follows
        \begin{align*}
            \dfi{{a}}{y|x} > &{\dfi{{a}}{y}} \Rightarrow\\
            (y-{a}x)^2 &< -2 \s^2_v \ln \left(\sqrt{2 \pi} \s_v \dfi{{a}}{y} \right).

        \end{align*}
        Taking the square roots, we obtain
        \begin{align*}

            \abs{y-{a}x} &< \sqrt{-2 \s^2_v \ln \left(\sqrt{2 \pi} \s_v \dfi{{a}}{y} \right)} \Rightarrow\\
            \abs{{a}x} &< \abs{y}+ \sqrt{-2 \s^2_v \ln \left(\sqrt{2 \pi} \s_v \dfi{{a}}{y} \right)}.

        \end{align*}
        Taking the square of both sides of the previous inequality and using the inequality  $(a+b)^2 \leq 2(a^2+b^2)$, we obtain
        \begin{align*}

            x^2 &< 2\frac{y^2}{{a}^2}-4 \frac{\s^2_v}{{a}^2} \ln \left(\sqrt{2 \pi} \s_v \dfi{{a}}{y} \right) \Rightarrow\\
            x^2 &< 2\frac{y^2}{{a}^2}+ 4 \frac{\s^2_v}{{a}^2}  \int \frac{(y-{a}x)^2}{2 \s^2_v} \dfi{}{x} dx  \Rightarrow\\
            x^2 &< 4\frac{y^2}{{a}^2}+ 2 \s^2_x
        \end{align*}
        By substituting for $x^2$ in  the second term of the RHS of Eq. \eqref{eqn:expectation:break},
        we conclude Eq. \eqref{eqn:lemma:conditional:second:moment}.
        \end{proof}

        As a result of Propositions \ref{lemma:upperbound:pointwise:regret} and \ref{lemma:upperbound:conditional:second:moment}, we obtain
        \begin{align}
           \nonumber(\phi_{\hat{a}}(y)- \phi_a(y))^2 & \leq \\ & 2 (6 \s^2_x + 4 \frac{y^2}{\hat{a}^2}+4 \frac{y^2}{{a}^2}) \kullback{\probdistn{\hat{a}|y}}{\probdistn{{a}|y}}.
           \label{eqn:lemma:basic}
        \end{align}
        Moreover, the following equality is known between Kullback-Leibler distance and Fisher information \cite[p.55]{Kullback:InformationTheory}
        \begin{align}
            \nonumber \kullback{\probdistn{\hat{a}|y}}{\probdistn{{a}|y}} = \frac{(\hat{a}-a)^2}{2}& \fisher{X;a||Y=y}\\&+ o(\hat{a}-a)^2,
            \label{eqn:kulback:fisher}
        \end{align}
         where $ \fisher{X;a||Y=y}  \triangleq \Expi{}{\big(\nabla \ln \dfi{a}{X|Y}\big)^2|Y=y} $
        is the Fisher information of $X$ relative to $a$, conditioned on $Y=y$.
        Substitute Eq. \eqref{eqn:kulback:fisher} in Eq. \eqref{eqn:lemma:basic} and
        note that \[\frac{1}{\hat{a}^2}=\frac{1}{a^2}+o(\hat{a}-a)^2\]
        for $|\frac{\hat{a}-a}{a}|<<1$.
        Taking the expectation
        with respect to $Y$, we conclude the proof of Lemma \ref{lemma:upperbound:regret}.

    \subsection{Proof of Lemma \ref{lemma:regret:Cramer-Rao}}
        \label{app:prf:regret:Cramer-Rao}
        We know that $\hat{a}=A_n(Y^{n-1})$. For an unbiased estimator and for sufficiently large values of $n$, $|A_n(Y^{n-1})-a|<<1$ and
        \begin{align}
            \nonumber \regret{A_n(Y^{n-1}),a} & \leq (A_n(Y^{n-1})-a)^2 \\&\Exp{\left(6 \s^2_x + 8 \frac{Y^2}{{a}^2}\right) \fisher{X;a||Y}}.
            \label{eqn:upperbound:regret:small:deviation}
        \end{align}
        holds true with arbitrarily high probability. Taking the expectation
        of both sides of Eq. \eqref{eqn:upperbound:regret:small:deviation} with respect to $Y^{n-1}$, we obtain
        \begin{align*}
            \nonumber  & \Exp{\regret{A_n(Y^{n-1}), a}} \leq  \\
            &  \Exp{(A_n(Y^{n-1})-a)^2} \; \Exp{\left(6 \s^2_x + 8 \frac{Y^2}{{a}^2}\right)  \fisher{X;a||Y}}

        \end{align*}
        Take the infimum of both sides over $\mcal A$ and assume
        $\mcal A$ contains an {\em efficient estimator} \cite[p. 92]{Borovkov:MathematicalStatistics}.
        By definition an efficient estimator achieves the Cramer-Rao bound. This means
        \begin{align*}
            \Exp{(A_n(Y^{n-1})-a)^2} = \frac{1}{\fisher{Y^{n-1};a}}.
        \end{align*}
        Since $Y_n$ is i.i.d., by additivity of Fisher information \[\fisher{Y^{n-1};a} = (n-1) \fisher{Y;a}.\]
        As a result, we obtain
        \begin{align*}
            \nonumber \inf_{A\in \mcal A} \Exp{\regret{A_n(Y^{n-1}),  a}} \leq \qquad \qquad \qquad\\ \frac{1}{n-1}  \frac{\Exp{\left(6 \s^2_x + 8 \frac{Y^2}{{a}^2}\right)  \fisher{X;a||Y}}}{\fisher{Y;a}}.
        \end{align*}

    \subsection{Proof of Lemma \ref{lemma:upperbound:relative:regret}}
    \label{app:prf:lemma:upperbound:relative:regret}
        By Proposition \ref{lemma:upperbound:pointwise:regret}, we have
        \begin{align*}
           \frac{ (\phi_{\hat{a}}(y)- \phi_a(y))^2}{\Expi{\hat{a}}{X^2|y}+\Expi{a}{X^2|y}}   &\leq   2 \kullback{P_{\hat{a}|y}}{P_{a|y}}.
        \end{align*}
        Substituting from Eq. \eqref{eqn:kulback:fisher} and taking the average with respect to $Y$, we obtain
        \begin{align*}
             \rregret{\hat{a}, a}   \leq   (\hat{a}-a)^2  \Expi{}{\big(\nabla\ln \dfi{a}{X|Y}\big)^2}+ o(\hat{a}-a)^2,
        \end{align*}
        and conclude the proof.

   \subsection{Proof of Lemma \ref{lemma:tradeoff}}
        \label{prf:lemma:tradeoff}

        Since $X$ does not depend on $a$, $\fisher{X;a} = 0 $, and hence
        \begin{eqnarray*}
            \r(a)=\frac{\fisher{X;a|Y}}{\fisher{Y;a}} = \frac{\fisher{Y;a|X}}{\fisher{Y;a}}-1.
        \end{eqnarray*}
        Moreover, since the additive noise is Gaussian, the equality
        \begin{eqnarray*}
            \fisher{Y;a|X} = \frac{\s^2_x}{\s^2_v}
        \end{eqnarray*}
        holds true for every distribution $P(X)$ with zero mean. As a result, we
        obtain Eq. \eqref{eqn:trade-off}.

\vspace{10pt}
    \section*{Acknowledgement}
        The author is thankful to the anonymous reviewers whose
        comments and suggestions improved the presentation of the paper.
\vspace{10pt}


text:
\chapter{Fast and Stable Signal Deconvolution via Compressible State-Space Models} \label{chap:css}
\chaptermark{Compressible State-Space Deconvolution}

Common biological measurements are in the form of noisy convolutions of signals of interest with possibly unknown and transient blurring kernels. Examples include EEG and calcium imaging data. Thus, signal deconvolution of these measurements is crucial in understanding the underlying biological processes. The objective of this chapter is to develop fast and stable solutions for signal deconvolution from noisy, blurred and undersampled data, where the signals are in the form of discrete events distributed in time and space. Due to their smoothness properties, in these application Gaussian state-space models exhibit poor performance in recovering the states and the sharp transitions. In this chapter we consider the problem of estimating compressible state space models, where the sparsity lies in the transitions of the states.

We introduce compressible state-space models as a framework to model and estimate such discrete events. These state-space models admit abrupt changes in the states and have a convergent transition matrix, and are coupled with compressive linear measurements. We consider a dynamic compressive sensing optimization problem and develop a fast solution, using two nested Expectation Maximization algorithms, to jointly estimate the states as well as their transition matrices. Under suitable sparsity assumptions on the dynamics, we prove optimal stability guarantees for the recovery of the states and present a method for the identification of the underlying discrete events with precise confidence bounds. We present simulation studies as well as application to calcium deconvolution and sleep spindle detection,  which verify our theoretical results and show significant improvement over existing techniques. Our results show that by explicitly modeling the dynamics of the underlying signals, it is possible to construct signal deconvolution solutions that are scalable, statistically robust, and achieve high temporal resolution.  Our proposed methodology provides a framework for modeling and deconvolution of noisy, blurred, and undersampled measurements in a fast and stable fashion, with potential application to a wide range of biological data.

\section{Introduction}

In many signal processing applications such as estimation of brain activity from magnetoencephalography (MEG) time-series \cite{phillips1997meg}, estimation of time-varying networks \cite{kolar2010estimating},  electroencephalogram (EEG) analysis \cite{nunez1995neocortical}, calcium imaging \cite{vogelstein2010fast}, functional magnetic resonance imaging (fMRI) \cite{chang2010time}, and video compression \cite{jung2010motion}, the signals often exhibit abrupt changes that are blurred through convolution with unknown kernels due to intrinsic measurement constraints. Extracting the underlying signals from blurred and noisy measurements is often referred to as signal deconvolution. Traditionally, state-space models have been used for such signal deconvolution problems, where the states correspond to the unobservable signals. Gaussian state-space models in particular are widely used to model smooth state transitions. Under normality assumptions, posterior mean filters and smoothers are optimal estimators, where the analytical solution is given respectively by the Kalman filter and the fixed interval smoother \cite{anderson1979optimal,haykin2008adaptive}.

When applied to observations from abruptly changing states, Gaussian state-space models exhibit poor performance in recovering sharp transitions of the states due to their underlying smoothing property. Although filtering and smoothing recursions can be obtained in principle for non-Gaussian state-space models, exact calculations are no longer possible. Apart from crude approximations like the extended Kalman filter, several methods have been proposed including numerical methods for low-dimensional states \cite{kitagawa1998self}, Monte Carlo filters \cite{kitagawa1998self,hurzeler1998monte}, posterior mode estimation \cite{fruhwirth1994applied,fruhwirth1994data}, and fully Bayesian smoothing using Markov chain Monte Carlo simulation \cite{knorr1999conditional, shephard1997likelihood}. In order to exploit sparsity, several dynamic compressed sensing (CS) techniques, such as the Kalman filtered CS algorithm, have been proposed that typically assume partial information about the sparse support or estimate it in a greedy and online fashion \cite{vaswani2010ls, vaswani2008kalman, carmi2010methods, ziniel2013dynamic, zhan2015time}. However, little is known about the theoretical performance guarantees of these algorithms.

In this chapter, we consider the problem of estimating state dynamics from noisy and undersampled observations, where the state transitions are governed by autoregressive models with compressible innovations. Motivated by the theory of CS, we employ an objective function formed by the $\ell_1$-norm of the state innovations \cite{ba2012exact}. Unlike the traditional compressed sensing setting, the sparsity is associated with the dynamics and not the states themselves. In the absence of observation noise, the CS recovery guarantees are shown to extend to this problem \cite{ba2012exact}. However, in a realistic setting in the presence of observation noise, it is unclear how the CS recovery guarantees generalize to this estimation problem.

We will present stability guarantees for this estimator under a convergent state transition matrix, which confirm that the CS recovery guarantees can be extended to this problem. The corresponding optimization problem in its Lagrangian form is akin to the MAP estimator of the states in a linear state-space model where the innovations are Laplace distributed. This allows us to integrate methods from Expectation-Maximization (EM) theory and Gaussian state-space estimation to derive efficient algorithms for the estimation of states as well as the state transition matrix, which is usually unknown in practice. To this end, we construct two nested EM algorithms in order to jointly estimate the states and the transition matrix. The outer EM algorithm for state estimation is akin to the fixed interval smoother, and the inner EM algorithm uses the state estimates to update the state transition matrix \cite{shumway1982approach}. The resulting EM algorithm is recursive in time, which makes the computational complexity of our method scale linearly with temporal dimension of the problem. This provides an advantage over existing methods based on convex optimization, which typically scale super-linearly with the temporal dimension.

Our results are related to parallel applications in spectral estimation, source localization, and channel equalization \cite{fevrier1999reduced}, where the measurements are of the form $\mathbf{Y} = \mathbf{A} \mathbf{X}+ \mathbf{N}$, with $\mathbf{Y}$ is the observation matrix, $\mathbf{X}$ denotes the unknown parameters, $\mathbf{A}$ is the measurement matrix, and $\mathbf{N}$ is the additive noise. These problems are referred to as Multiple Measurement Vectors (MMV) \cite{cotter2005sparse} and Multivariate Regression \cite{obozinski2011support}. In these applications, solutions with row sparsity in $\mathbf{X}$ are desired. Recovery of sparse signals with Gaussian innovations is studied in \cite{zhang2011sparse}. Several recovery algorithms including the $\ell_1\!\!\!-\!\!\ell_q$ minimization methods, subspace methods and greedy pursuit algorithms \cite{davies2012rank} have been proposed for support union recovery in this setup. Our contributions are distinct in that we directly model the state innovations as a compressible sequence, for recovery of which we present both sharp theoretical guarantees as well as fast algorithms from state-space estimation.

Finally, we provide simulation results as well as applications to two experimentally-acquired data sets: calcium imaging recordings of neuronal activity, and EEG data during sleep. In the former, the deconvolution problem concerns estimating the location of spikes given the temporally blurred calcium fluorescence, and in the latter, the objective is to detect the occurrence and onset of sleep spindles. Our simulation studies confirm our theoretical predictions on the performance gain obtained by compressible state-space estimation over those obtained by traditional estimators such as the basis pursuit denoising. Our real data analyses reveal that our compressible state-space modeling and estimation framework outperforms two of the commonly-used methods for calcium deconvolution and sleep spindle detection. In the spirit of easing reproducibility, we have made MATLAB implementations of our codes publicly available \cite{code}.

\vspace{-2mm}
\section{Methods}\label{sec:tv_formulation}
In this section we introduce the experimental procedures of recording the data for analysis and establish our problem formulation and notational conventions and present our main theoretical analysis and algorithm development.
\subsection{Experimental Procedures}
\subsubsection{Surgery} 2 hours before surgery, $0.1~\text{cc}$ dexamethasone ($2~\text{mg/ml}$, VetOne) was injected subcutaneously to reduce brain swelling during craniotomy. Anesthesia is induced with $4\%$ isoflurane (Fluriso, VetOne) with a calibrated vaporizer (Matrx VIP 3000). During surgery, isoflurane level was reduced to and maintained at a level of $1.5\%$--$2\%$. Body temperature of the animal is maintained at $36.0$ degrees Celsius during surgery. Hair on top of head of the animal was removed using Hair Remover Face Cream (Nair), after which Betadine (Purdue Products) and $70\%$ ethanol was applied sequentially 3 times to the surface of the skin before removing the skin. Soft tissues and muscles were removed to expose the skull. Then a custom designed 3D printed stainless headplate was mounted over left auditory cortex and secured with C\&B-Metabond (Parkell). A craniotomy with a diameter of around $3.5~\text{mm}$ was then performed over left auditory cortex. A three layered cover slip was used as cranial window, which is made by gluing (NOA71, Norland Products) $2$ pieces of $3~\text{mm}$ coverslips (64-0720 (CS-3R), Warner Instruments) with a $5~\text{mm}$ coverslip ($64$--$0700$ (CS-5R), Warner Instruments). Cranial window was quickly dabbed in kwik-sil (World Precision Instruments) before mounted $3~\text{mm}$ coverslips facing down onto the brain. After kwik-sil cured, Metabond was applied to secure the position of the cranial window. Synthetic Black Iron Oxide (Alpha Chemicals) was then applied to the hardened Metabond surface. $0.05~\text{cc}$ Cefazolin ($1$ gram/vial, West Ward Pharmaceuticals) was injected subcutaneously when entire procedure was finished. After the surgery the animal was kept warm under heat light for $30$ minutes for recovery before returning to home cage. Medicated water (Sulfamethoxazole and Trimethoprim Oral Suspension, USP $200~\text{mg}/40~\text{mg}$ per $5~\text{ml}$, Aurobindo Pharms USA; $6~\text{ml}$ solution diluted in $100~\text{ml}$ water) substitute normal drinking water for $7$ days before any imaging was performed.

\subsubsection{Awake two-photon imaging} Spontaneous activity data of population of layer $2/3$ auditory cortex (A1) neurons is collected from adult (3-month old) Thy1-GCaMP6s female mouse implanted with chronic window following the above procedure, using two-photon imaging. Acquisition is performed using a two-photon microscope (Thorlabs Bscope 2) equipped with a Vision 2 Ti:Sapphire laser (Coherent), equipped with a GaAsP photo detector module (Hamamatsu) and resonant scanners enabling faster high-resoluation scanning at $30$--$60~\text{Hz}$ per frame. The excitation wavelength was $920~\text{nm}$. Regions ($\sim 300~\mu\text{m}^2$) within A1 were scanned at $30~\text{Hz}$ through a $20\text{x}$, $0.95$ NA water-immersion objective (Olympus). During imaging the animal was head-fixed and awake. The microscope was rotated $45$ degrees and placed over the left A1 where window was placed.  An average image of field of view was generated by choosing a time window where minimum movement of the brain was observed and used as reference image for motion correction using TurboReg plugin in ImageJ. GCaMP6s positive cells are selected manually by placing a ring like ROI over each identified cell. Neuropil masks were generated by placing a $20~\mu\text{m}$ radius circular region over each cell yet excluding all cell soma regions. Traces of soma and neuropil were generated by averaging image intensity within respective masks at each time point. A ratio of $0.7$ was used to correct for neuropil contamination.

\subsubsection{Cell-attached patch clamp recordings and two-photon imaging} Recordings were performed in vitro in voltage clamp to simultaneously measure spiking activity and $\Delta F / F$. Thalamocortical slices containing A1 were prepared as previously described \cite{meng2015visual}. The extracellular recording solution consisted of artificial cerebral spinal fluid (ACSF) containing: $130$ NaCl, $3$ KCl, $1.25$ KH2PO4, $20$ NaHCO3, $10$ glucose, $1.3$ MgSO4, $2.5$ CaCl2 (pH $7.35$-$7.4$, in $95\%$ O2 – $5\%$ CO2). Action potentials were recorded extracellularly in loose-seal cell-attached configuration (seal resistance typically $20$--$30~\text{M}\Omega$) in voltage clamp mode. Borosilicate glass patch pipettes were filled with normal ACSF diluted $10\%$, and had a tip resistance of $\sim3$-$5~\text{M}\Omega$ in the bath. Data were acquired with a Multiclamp 700B patch clamp amplifier (Molecular Devices), low-pass filtered at $3$-$6~\text{kHz}$, and digitized at $10~\text{kHz}$ using the MATLAB-based software. Action potentials were stimulated with a bipolar electrode placed in L1 or L23 to stimulate the apical dendrites of pyramidal cells (pulse duration $1$-$5~\text{ms}$). Data were analyzed offline using MATLAB. Imaging was largely performed using a two-photon microscope (Ultima, Prairie Technologies) and a MaiTai DeepSee laser (SpectraPhysics), equipped with a GaAsP photo detector module (Hamamatsu) and resonant scanners enabling faster high-resoluation scanning at $30$-$60~\text{Hz}$ per frame. Excitation was set at $900~\text{nm}$. Regions were scanned at $30~\text{Hz}$ through a $40\text{x}$ water-immersion objective (Olympus). Cells were manually selected as ring-like regions of interest (ROIs) that cover soma but exclude cell nuclei, and pixel intensity within each ROI was averaged to generate fluorescence over time and changes in fluorescence ($\Delta F / F$) were then calculated.

\vspace{-2mm}
\subsection{Problem Formulation and Theoretical Analysis}

We consider the linear compressible state-space model given by
\begin{equation}
\label{eq:tv_lap_state_space}
\mathbf{x}_t = \mathbf{\Theta} \mathbf{x}_{t-1}+ \mathbf{w}_t, \qquad \mathbf{y}_t = \mathbf{A}_t \mathbf{x}_t + \mathbf{v}_t,
\end{equation}
where $(\mathbf{x}_t)_{t=1}^{T} \in \mathbb{R}^p$ denote the sequence of unobservable states, $\mathbf{\Theta}$ is the state transition matrix satisfying $\|\mathbf{\Theta}\|<1$, $\mathbf{w}_t \in \mathbb{R}^p$ is the state innovation sequence, $(\mathbf{y}_t)_{t=1}^T \in \mathbb{R}^{n_t}$ are the linear observations, $\mathbf{A}_t \in \mathbb{R}^{n_t \times p}$ denotes the measurement matrix, and {$\mathbf{v}_t \in \mathbb{R}^{n_t}$} denotes the measurement noise. The main problem is to estimate the unobserved sequence $(\mathbf{x}_t)_{t=1}^T$ (and possibly $\boldsymbol{\Theta}$), given the sequence of observations $(\mathbf{y}_t)_{t=1}^T$. This problem is in general ill-posed, when $n_t < p$, for some $t$. We therefore need to make additional assumptions in order to seek a stable solution.

We assume that the state innovations are sparse (resp. compressible), i.e. $\mathbf{x}_t-\mathbf{\Theta} \mathbf{x}_{t-1}$ is $s_t$-sparse (resp. $(s_t,\xi)$-compressible) with $s_1 \gg s_t$ for $t \in [T]\backslash\{1\}$. Our theoretical analysis pertain to the compressed sensing regime where $1 \ll s_t < n_t \ll p$.  We assume that the rows of $\mathbf{A}_t$ are a subset of the rows of $\mathbf{A}_1$, i.e. $\mathbf{A}_t = (\mathbf{A}_{1})_{n_t}$, and define $\widetilde{\mathbf{A}}_t = \sqrt{\frac{n_1}{n_t}}\mathbf{A}_t$. Other than its technical usefulness, the latter assumption helps avoid prohibitive storage of all the measurement matrices. In order to promote sparsity of the state innovations, we consider the dynamic $\ell_1$-regularization (dynamic CS from now on) problem defined as
\begin{align}
\label{eq:tv_prob_def_primal}
\minimize \limits_{(\mathbf{x}_t)_{t=1}^{T},\mathbf{\Theta}} \quad \displaystyle \sum_{t=1}^T \frac{\|\mathbf{x}_t-\mathbf{\Theta} \mathbf{x}_{t-1}\|_1}{\sqrt{s_t}} \ \ \text{s.t.} \ \ \|\mathbf{y}_t-\mathbf{A}_t\mathbf{x}_t\|_2 \leq \sqrt{\frac{n_t}{n_1}}\varepsilon.
\end{align}
where $\varepsilon$ is an upper bound on the observation noise, i.e., $\| v_t \|_2 \le \varepsilon$ for all $t$. Note that this problem is a variant of the dynamic CS problem introduced in \cite{ba2012exact}. We also consider the modified Lagrangian form of (\ref{eq:tv_prob_def_primal}) given by
\begin{equation}
\label{eq:tv_prob_def_dual}
\minimize_{(\mathbf{x}_t)_{t=1}^{T},\mathbf{\Theta}} \quad \lambda \sum_{t=1}^T \frac{\|\mathbf{x}_t-\mathbf{\Theta} \mathbf{x}_{t-1}\|_1}{\sqrt{s_t}} +  \frac{1}{n_t}\frac{\|\mathbf{y}_t-\mathbf{A}_t\mathbf{x}_t\|_2^2}{2\sigma^2}.
\end{equation}
for some constants $\sigma_2^2$ and $\lambda \ge 0$. Note that if $\mathbf{v}_t \sim \mathcal{N}(\mathbf{0}, n_t \sigma^2 \mathbf{I})$, then Eq. (\ref{eq:tv_prob_def_dual}) is akin to the maximum \emph{a posteriori} (MAP) estimator of the states in (\ref{eq:tv_lap_state_space}), assuming that the \textit{state} innovations were independent Laplace random variables with respective parameters $\lambda / \sqrt{s_t}$. We will later use this analogy to derive fast solutions to the optimization problem in (\ref{eq:tv_prob_def_dual}).

Uniqueness and exact recovery of the sequence $(\mathbf{x}_t)_{t=1}^T$ in the absence of noise was proved in \cite{ba2012exact} for $\mathbf{\Theta} =\mathbf{I}$, by an inductive construction of dual certificates. The special case $\mathbf{\Theta} =\mathbf{I}$ can be considered as a generalization of the total variation (TV) minimization problem \cite{poon2015role}. Our main result on stability of the solution of (\ref{eq:tv_prob_def_primal}) is the following:

\begin{thm}[Stable Recovery of Activity in the Presence of Noise]
\label{thm:tv_main}
Let $(\mathbf{x}_t)_{t=1}^T \in \mathbb{R}^p$ be a sequence of states with a known transition matrix $\mathbf{\Theta}=\theta \mathbf{I}$, where $|\theta|<1$ and  $\widetilde{\mathbf{A}}_t$, $t\geq 1$ satisfies RIP of order ${4s}$ with $\delta_{4s} <1/3$. Suppose that $n_1 > n_2 = n_3 = \cdots = n_T$. Then, the solution $(\widehat{\mathbf{x}}_t)_{t=1}^{T}$ to the dynamic CS problem (\ref{eq:tv_prob_def_primal}) satisfies
\begin{align}
\label{eq:tv_main_stable}
 \frac{1}{T} \sum_{t=1}^T \|\mathbf{x}_t-\widehat{\mathbf{x}}_t\|_2 \le\displaystyle \frac{1-\theta^T}{1-\theta}\left(12.6 \left( 1+\frac{1}{T} \sqrt{\frac{n_1}{n_2}}-\frac{1}{T} \right)\varepsilon + \frac{3}{T} \sum_{t=1}^T \frac{\sigma_{s_t}(\mathbf{x}_t-\mathbf{\Theta} \mathbf{x}_{t-1})}{\sqrt{s_t}}\right).
\end{align}
\end{thm}

\noindent \textit{\textbf{Proof Sketch.}} The proof of Theorem \ref{thm:tv_main} is based on establishing a modified cone and tube constraint for the dynamic CS problem (\ref{eq:tv_prob_def_primal}) and using the boundedness of the Frobenius norm of the inverse first-order differencing operator. Details of the proof are given in Appendix \ref{app:tv_main_proof}.

\noindent { \textbf{Remark 1.}} The first term on the right hand side of Theorem \ref{thm:tv_main} implies that the average reconstruction error of the sequence $(\mathbf{x}_t)_{t=1}^T$ is upper bounded proportional to the noise level $\varepsilon$, which implies the stability of the estimate. The second term is a measure of compressibility of the innovation sequence and vanishes when the sparsity condition is exactly met.

\noindent \textbf{Remark 2.} Similar to the general compressive sensing setup, the results of Theorem \ref{thm:tv_main} are general and hold even when the sparsity assume is not exactly met, namely when the compressibility (sparsity) assumption is violated, the second term on the right hand side could be of order of a constant. However in order for $\widetilde{\mathbf{A}}_t$ to satisfy RIP of order $4s$ one requires $n \sim \mathcal{O}(s \log p)$ which in the absence of sparsity ($s \sim \mathcal{O}(p)$) reduces to the denoising regime.

\noindent \textbf{Remark 3.} The theory of LASSO suggests a choice of the regularization parameter ${\lambda}  \sim \mathcal{O} \left(\sigma \sqrt{\frac{\log p}{n_t}}\right)$ \cite{Negahban}.  In the absence of sparsity (compressibility) $s \sim \mathcal{O}(n)$, one can think of two possible scenarios: In the high-dimensional setup where $n \gg 1$ , $\lambda$ will be very small and (\ref{eq:tv_prob_def_dual}) will be similar to solving the Maximum-Likelihood problem which is known to result in an unbiased estimator. In the low-dimensional case $\lambda$ is chosen via cross-validation which adapts the problem to the sparsity level.

\noindent \textbf{Remark 4.} Finally, the choice of the $\ell_1$-regularized maximum-likelihood estimation in (\ref{eq:tv_prob_def_dual}) is motivated by its superior performance over ML estimation in the compressive sensing literature and the applications of interest in this chapter where exhibit sparse activity patterns in time. It is well-known that Laplace distribution is not a compressible distribution \cite{gribonval2012compressible}. Therefore we have not made the assumption that the innovations $\mathbf{w}_t$ are Laplace distributed, but have made the observation that an $\ell_1$-regularized ML estimator is akin to the MAP estimator for a Laplace state-space model. Theorem \ref{thm:tv_main} suggests that a Laplace state-space model is asymptotically \textit{as good as} any other compressible state-space model up to a constant factor in the error bounds.

\vspace{-2mm}
\subsection{Fast Iterative Solution via the EM Algorithm}
\label{sec:tv_algorithm}
Due to the high dimensional nature of the state estimation problem, algorithms with polynomial complexity exhibit poor scalability. Moreover, when the state transition matrix is not known, the dynamic CS optimization problem (\ref{eq:tv_prob_def_dual}) is not convex in $\left( (\mathbf{x}_t)_{t=1}^T,\mathbf{\Theta} \right)$. Therefore standard convex optimization solvers cannot be directly applied. This problem can be addressed by employing the Expectation-Maximization algorithm \cite{shumway1982approach}. A related existing result considers weighted $\ell_1$-regularization to adaptively capture the state dynamics \cite{charles2013dynamic}. Our approach is distinct in that we derive a fast solution to (\ref{eq:tv_prob_def_dual}) via two nested EM algorithms, in order to jointly estimate the states and their transition matrix. The outer EM algorithm converts the estimation problem to a form suitable for the usage of the traditional Fixed Interval Smoothing (FIS) by invoking the EM interpretation of the Iterative Re-weighted Least Squares (IRLS) algorithms \cite{babadi_IRLS}. The inner EM algorithm performs state and parameter estimation efficiently using the FIS. We refer to our estimated as the Fast Compressible State-Space ({\sf FCSS}) estimator.

\subsubsection*{The outer EM loop of {\sf FCSS}}

In \cite{babadi_IRLS}, the authors established the equivalence of the IRLS algorithm as an instance of the EM algorithm for solving $\ell_1$-minimization problems via the Normal/Independent (N/I) characterization of the Laplace distribution. Consider the $\epsilon$-perturbed $\ell_1$-norm as
\begin{equation}
\label{eq:tv_eps_l1}
\|\mathbf{x}\|_{1,\epsilon} = \sqrt{x_1^2 + \epsilon^2}+\sqrt{x_2^2 + \epsilon^2}+\cdots+\sqrt{x_p^2 + \epsilon^2}.
\end{equation}
Note that for $\epsilon =0$, $\|\mathbf{x}\|_{1,\epsilon}$ coincides with the usual $\ell_1$-norm. We define the $\epsilon$-perturbed version of the dual problem (\ref{eq:tv_prob_def_dual}) by
\begin{equation}
\label{eq:tv_eps_lag}
\minimize_{\left(\mathbf{x}_t\right)_{t=1}^T, \mathbf{\Theta}} \quad \lambda \sum_{t=1}^T \frac{\|\mathbf{x}_t-\mathbf{\Theta} \mathbf{x}_{t-1}\|_{1,\epsilon}}{\sqrt{s_t}} +  \frac{1}{n_t}\frac{\|\mathbf{y}_t-\mathbf{A}_t\mathbf{x}_t\|_2^2}{2\sigma^2}.
\end{equation}
If instead of the $\ell_{1,\epsilon}$-norm, we had the square $\ell_2$ norm, then the above problem could be efficiently solved using the FIS. The outer EM algorithm indeed transforms the problem of Eq. (\ref{eq:tv_eps_lag}) into this form. Note that the $\epsilon$-perturbation only adds a term of the order $\mathcal{O}(\epsilon p)$ to the estimation error bound of Theorem \ref{thm:tv_main}, which is negligible for small enough $\epsilon$ \cite{babadi_IRLS}.

The problem of Eq. (\ref{eq:tv_eps_lag}) can be interpreted as a MAP problem: the first term corresponds to the state-space prior $-\log$ $ p(\mathbf{x}_t | \mathbf{x}_{t-1}, \boldsymbol{\Theta}) = -\log p_{s_t}(\mathbf{x}_t -  \boldsymbol{\Theta} \mathbf{x}_{t-1})$, where $p_{s_t}(\mathbf{x}) \sim \exp$~$(-\lambda \|\mathbf{x}\|_{1,\epsilon} / \sqrt{s_t} )$ denoting the $\epsilon$-perturbed Laplace distribution; the second term is the negative log-likelihood of the data given the state, assuming a zero-mean Gaussian observation noise with covariance $\sigma^2 \mathbf{I}$. Suppose that at the end of the $l^{\sf th}$ iteration, the estimates $(\widehat{\mathbf{x}}_t^{(l)})_{t=1}^T,\widehat{\mathbf{\Theta}}^{(l)}$ are obtained, given the observations $(\mathbf{y}_t)_{t=1}^T$. As it is shown in Appendix \ref{app:tv_EM}, the outer EM algorithm transforms the optimization problem to:
\begin{align}
\label{eq:tv_prob_def_dual_irls}
\minimize_{\left(\mathbf{x}_t\right)_{t=1}^T, \mathbf{\Theta}} \quad &\frac{\lambda}{2} \sum_{j=1}^p\sum_{t=1}^T \frac{\left( \mathbf{x}_{t}-\mathbf{\Theta} \mathbf{x}_{t-1}\right)_j^2+\epsilon^2}{\sqrt{s_t}\sqrt{\left( \widehat{\mathbf{x}}_{t}^{(l)}-\widehat{\mathbf{\Theta}}^{(l)} \widehat{\mathbf{x}}_{t-1}^{(l)}\right)_j^2+\epsilon^2}}  +  \sum_{t=1}^T \frac{1}{n_t}\frac{\|\mathbf{y}_t-\mathbf{A}_t\mathbf{x}_t\|_2^2}{2\sigma^2}.
\end{align}
in order to find $(\widehat{\mathbf{x}}_t^{(l+1)})_{t=1}^T,\widehat{\mathbf{\Theta}}^{(l+1)}$. Under mild conditions, convergence of the solution of (\ref{eq:tv_prob_def_dual_irls}) to that of (\ref{eq:tv_prob_def_dual}) was established in \cite{babadi_IRLS}. The objective function of (\ref{eq:tv_prob_def_dual_irls}) is still not jointly convex in $\left( (\mathbf{x}_t)_{t=1}^T,\mathbf{\Theta} \right)$. Therefore, to carry out the optimization, i.e. the outer M step, we will employ another instance of the EM algorithm, which we call the inner EM algorithm, to alternate between estimating of $\left(\mathbf{x}_t\right)_{t=1}^T$ and $\mathbf{\Theta}$.

\subsubsection*{The inner EM loop of {\sf FCSS}}

Let $\mathbf{W}_t^{(l)}$ be a diagonal matrix such that
\[
(\mathbf{W}_t^{(l)})_{j,j} = {s_t}^{-1/2} \left \{\left( \widehat{\mathbf{x}}_t^{(l)} -{\widehat{\mathbf{\Theta}}}^{(l)} \widehat{\mathbf{x}}_{t-1}^{(l)}\right)_j^2+\epsilon^2\right\}^{-1/2}.
\]
Consider an estimate $\widehat{\mathbf{\Theta}}^{(l,m)}$, corresponding to the $m^{\sf th}$ iteration of the inner EM algorithm within the $l^{\sf th}$ M-step of the outer EM. In this case, Eq. (\ref{eq:tv_prob_def_dual_irls}) can be thought of the MAP estimate of the Gaussian state-space model given by:
\begin{align}
\label{eq:tv_dynamic}
\begin{array}{l}
\mathbf{x}_t = \widehat{\mathbf{\Theta}}^{(l,m)} \mathbf{x}_{t-1}+ \mathbf{w}_t, \quad \mathbf{w}_t\sim \mathcal{N}\left(\mathbf{0}, \frac{1}{\lambda} {\mathbf{W}_t^{(l)}}^{-1} \right) \\
\mathbf{y}_t = \mathbf{A}_t \mathbf{x}_t + \mathbf{v}_t, \quad \mathbf{v}_t\sim \mathcal{N}(\mathbf{0},n_t\sigma^2 I)
\end{array},
\end{align}
In order to obtain the inner E step, one needs to find the density of $(\mathbf{x}_t)_{t=1}^T$ given $(\mathbf{y}_t)_{t=1}^T$ and $\widehat{\mathbf{\Theta}}^{(l,m)} \}$. Given the Gaussian nature of the state-space in Eq. (\ref{eq:tv_dynamic}), this density is a multivariate Gaussian density, whose means and covariances can be efficiently computed using the FIS. For all $t \in [T]$, the FIS performs a forward Kalman filter and a backward smoother to generate \cite{rauch1965maximum,anderson1979optimal}:
\begin{align*}
\mathbf{x}^{(l, m+1)}_{{t|T}} := {\mathbb{E}}\left\{\mathbf{x}_t|(\mathbf{y}_t)_{t=1}^T, \widehat{\mathbf{\Theta}}^{(l,m)}\right\},
\end{align*}
\begin{equation*}
\boldsymbol{\Sigma}^{(l,m+1)}_{t|T} := {\mathbb{E}}\left\{\mathbf{x}_t \mathbf{x}_t'|(\mathbf{y}_t)_{t=1}^T, \widehat{\mathbf{\Theta}}^{(l,m)}\right\},
\end{equation*}
and
\begin{equation*}
\mathbf{\Sigma}^{(l,m+1)}_{t-1,t|T}={\mathbf{\Sigma}}^{(l,m+1)}_{t,t-1|T}={\mathbb{E}}\left\{\mathbf{x}_{t-1} \mathbf{x}_t'|(\mathbf{y}_t)_{t=1}^T,\widehat{\mathbf{\Theta}}^{(l,m)}\right\}.
\end{equation*}

Note that due to the quadratic nature of all the terms involving $(\mathbf{x}_t)_{t=1}^T$, the outputs of the FIS suffice to compute the expectation of the objective function in Eq. (\ref{eq:tv_prob_def_dual_irls}), i.e., the inner E step, which results in:
\begin{align}
\label{eq:tv_Q_calculated}
\notag \maximize_{\mathbf{\Theta}} & -\frac{\lambda}{2}  \displaystyle \left ({\mathbf{\Theta}}\left( \sum_{t=1}^T   \mathbf{W}_t^{(l)} \left (\mathbf{x}^{(l, m+1)}_{t-1|T}{\mathbf{x}'}^{(l, m+1)}_{t-1|T} + \mathbf{\Sigma}_{t-1|T}^{(l,m+1)} \right )\right){\mathbf{\Theta}}^T\right)\\
&  + \frac{\lambda}{2} \operatorname{Tr} \displaystyle { \left ({\mathbf{\Theta}} \left (\sum_{t=1}^T \mathbf{W}_t^{(l)} \left ( \mathbf{x}^{(l, m+1)}_{t-1|T} {\mathbf{x}'}^{(l, m+1)}_{t|T} + \mathbf{x}^{(l, m+1)}_{t|T} {\mathbf{x}'}^{(l, m+1)}_{t-1|T} + 2\mathbf{\Sigma}_{t-1,t|T}^{(l,m+1)} \right ) \right) \right)},
\end{align}
to obtain $\widehat{\boldsymbol{\Theta}}^{(l,m)}$. The solution has a closed-form given by:
\begin{align}
\label{eq:tv_theta_upd}
\notag  \widehat{\mathbf{\Theta}}^{(l, m+1)} = \displaystyle & \left(   \sum_{t=1}^T 2 \mathbf{W}_t^{(l)}\left(\mathbf{x}^{(l, m+1)}_{t-1|T}{\mathbf{x}'}^{(l, m+1)}_{t-1|T} + \mathbf{\Sigma}_{t-1|T}^{(l,m+1)} \right)  \right)^{-1}\\
 &  \displaystyle \left( \sum_{t=1}^T \mathbf{W}_t^{(l)} \left(\mathbf{x}^{(l, m+1)}_{t-1|T} {\mathbf{x}'}^{(l, m+1)}_{t|T}+ \mathbf{x}^{(l, m+1)}_{t|T} {\mathbf{x}'}^{(l, m+1)}_{t-1|T} + 2 \mathbf{\Sigma}_{t-1,t|T}^{(l,m+1)}\right) \right).
\end{align}

\noindent \begin{minipage}{\columnwidth}
\begin{center}
\begin{algorithm}[H]

\caption{\small The Fast Compressible State-Space ({\sf FCSS}) Estimator}
\label{alg:DS3}
\begin{algorithmic}[1]
\Procedure{{\sf FCSS}}{}
\State Initialize: $\widehat{\mathbf{\Theta}}^{(0)}=\mathbf{0}$, $(\widehat{\mathbf{x}}_t^{(0)} = \mathbf{0})_{t=1}^T $, $(\mathbf{W}^{(0)}_t = \mathbf{0})_{t=1}^T$.
\Repeat
\State $l =0$ .

\State \textbf{Outer E-step:}
\State $\mathbf{W}_t^{(l)} = {\sf diag} \Bigg \{ \frac{1}{\sqrt{s_t}\sqrt{\left( \widehat{\mathbf{x}}_t^{(l)} -{\widehat{\mathbf{\Theta}}}^{(l)} \widehat{\mathbf{x}}_{t-1}^{(l)}\right)_j^2+\epsilon^2}} \Bigg \}_{j=1}^p$.
\State \textbf{Outer M-step:}
\Repeat
\State $m =0$ .
\State \textbf{Inner E-Step:} Find the smoothed estimates  $\mathbf{x}^{(l,m+1)}_{{t|T}}$,  $\boldsymbol{\Sigma}^{(l,m+1)}_{{t|T}}$ and $\boldsymbol{\Sigma}^{(l,m+1)}_{{t-1,t|T}}$ using a Fixed Interval Smoother for
\begin{align*}
\left \{
\begin{array}{l}
\mathbf{x}_t = \widehat{\mathbf{\Theta}}^{(l,m)} \mathbf{x}_{t-1}+ \mathbf{w}_t, \quad \mathbf{w}_t\sim \mathcal{N}\left(\mathbf{0}, \frac{1}{\lambda} {\mathbf{W}_t^{(l)}}^{-1} \right) \\
\mathbf{y}_t = \mathbf{A}_t \mathbf{x}_t + \mathbf{v}_t, \quad \mathbf{v}_t\sim \mathcal{N}(\mathbf{0},n_t\sigma^2 I)
\end{array} \right. .
\end{align*}

\State \textbf{Inner M-Step:} Update $\widehat{\mathbf{\Theta}}^{(l, m+1)}$ via Eq. (\ref{eq:tv_theta_upd}).
\State $m \leftarrow m+1$.
\Until{convergence criteria met}
\State Update the estimates:

\begin{equation*}
\widehat{\mathbf{\Theta}}^{(l+1)} \leftarrow \widehat{\mathbf{\Theta}}^{(l,m)}, \quad \left(\widehat{\mathbf{x}}_t^{(l+1)}\right)_{t=1}^{T} \leftarrow \left({\mathbf{x}}_{t|T}^{(l,m)}\right)_{t=1}^{T}.
\end{equation*}
\State $l \leftarrow l+1$.
\Until{convergence criteria met}

$\widehat{\mathbf{\Theta}} \leftarrow \widehat{\mathbf{\Theta}}^{(l)}.  \quad \left(\widehat{\mathbf{x}}_t \right)_{t=1}^{T} \leftarrow  \left(\widehat{\mathbf{x}}_t^{(l)}\right)_{t=1}^{T} $
\EndProcedure
\end{algorithmic}
\end{algorithm}
\end{center}
\end{minipage}

\medskip
\medskip
\medskip

\noindent \textbf{Remark 1.}
The EM algorithm is known to converge linearly in the number of iterations \cite{mclachlan2007algorithm}. In applications of interest in this chapter, the stability assumption results in finite time constants and the EM algorithm can be assumed to converge in finite number of steps. As a result the complexity of the {\sf FCSS} algorithm is equal to complexity per EM iteration.
 By virtue of the FIS procedure, the compexity of the {\sf FCSS} algorithm is \emph{linear} in $T$, i.e., the observation duration. As we will show in Section \ref{sec:tv_sim}, this makes the {\sf FCSS} algorithm scale favorably when applied to long data sets.  Each iteration of the FIS algorithm requires two inversions which is of complexity $\mathcal{O}(p^3)$. Similarly, updating $\mathbf{\Theta}$ requires an inversion which is of complexity $\mathcal{O}(p^3)$. Altogether, complexity of the {\sf FCSS} algorithm amounts to $\mathcal{O}(p^3 T)$. For some applications of interest in this chapter, such as calcium deconvolution in the denoising regime and sleep spindle detection, the inversion is performed on a diagonal matrix and hence the complexity reduces to $\mathcal{O}(T)$.

\noindent \textbf{Remark 2.} In order to update $\mathbf{\Theta}$ in the inner M-step given by E. (\ref{eq:tv_theta_upd}), we have not specifically enforced the condition $\|\mathbf{\Theta}\| <1$ in the maximization step. This condition is required to obtain a convergent state transition matrix which results in the stability of the state dynamics. It is easy to verify that the set of matrices $\mathbf{\Theta}$ satisfying $\|\mathbf{\Theta}\| < 1- \eta$, is a closed convex set for small positive $\eta$, and hence one can perform the maximization in (\ref{eq:tv_theta_upd}) by projection onto this closed convex set. Alternatively, matrix optimization methods with operator norm constraints can be used \cite{mishra2013low}. We have avoided this technicality by first finding the global minimum and examining the largest eigenvalue. In the applications of interest in this chapter which follow next, the largest eigenvalue has always been found to be less than $1$.

\vspace{-2mm}
\section{Results}
\label{sec:tv_sim}
In this section, we study the performance of the {\sf FCSS} estimator on simulated data as well real data from two-photon calcium imaging recordings of neuronal activity and sleep spindle detection from EEG.

\vspace{-2mm}
\subsection{Application to Simulated Data}
We first apply the {\sf FCSS} algorithm to simulated data and compare its performance with the Basis Pursuit Denoising (BPDN) algorithm. The parameters are chosen as $p = 200, T = 200, s_1 =8,s_2 = 4,\epsilon = 10^{-10}$, and $\mathbf{\Theta} = 0.95 \mathbf{I}$. We define the quantity $1- n/p$ as the compression ratio. We refer to the case of $n_t=p$, i.e., no compression, as the denoising setting. The measurement matrix $\mathbf{A}$ is an $n_t \times p$ i.i.d. Gaussian random matrix, where $n_t$ is chosen such that $\frac{s_t}{n_t}$ is a fixed ratio. An initial choice of ${\lambda}  \geq 2\sqrt{2} \sigma \sqrt{\frac{\log p}{n_t}}$ is made inspired by the theory of LASSO \cite{Negahban}, which is further tuned using two-fold cross-validation.

Figures \ref{fig:tv_ds}--(a) and \ref{fig:tv_ds}--(b) show the estimated states as well as the innovations for different compression ratios for one sample component. In the denoising regime, all the innovations (including the two closely-spaced components) are exactly recovered. As we take fewer measurements, the performance of the algorithm degrades as expected. However, the overall structure of the innovation sequence is captured even for highly compressed measurements.

\begin{figure}[htb!]
\vspace{-2mm}
\centering
\includegraphics[width=1\columnwidth]{figures/tv_ds-01}
\vspace{-2mm}
\caption{Reconstruction results of {\sf FCSS} on simulated data vs. compression levels. (a) reconstructed states, (b) reconstructed spikes. The {\sf FCSS} estimates degrade gracefully as the compression level increases.}\label{fig:tv_ds}
\vspace{-1mm}
\end{figure}

\begin{figure}[htb!]
\vspace{-2mm}
\centering
\includegraphics[width=.7\columnwidth]{figures/tv_mse-01}
\vspace{-2mm}
\caption{MSE vs. SNR comparison betwee {\sf FCSS} and BPDN. The {\sf FCSS} significantly outperforms the BPDN, even for moderate SNR levels.}
\label{fig:tv_mse_comp}
\vspace{-2mm}
\end{figure}

\begin{figure}[htb!]
\centering
\includegraphics[width=.7\columnwidth]{figures/tv_vs_bpdn-01}
\vspace{-2mm}
\caption{Example of state reconstruction results for {\sf FCSS} and BPDN. Top: true states, Middle: {\sf FCSS} state estimates, Bottom: BPDN estimates. The {\sf FCSS} reconstruction closely follows the true state evolution, while the BPDN fails to capture the state dynamics.}
\label{fig:tv_dcs_vs_ccs}

\end{figure}

\begin{figure}[htb!]
\centering
\includegraphics[width=0.9\columnwidth]{figures/imagesc_sim}
\caption{Raster plots of the observed and estimated states via {\sf FCSS}. Left: noisy observations, Right: {\sf FCSS} estimates. The {\sf FCSS} significantly denoises the observed states.}
\label{fig:imagesc_sim}
\vspace{-4mm}
\end{figure}

Figure \ref{fig:tv_mse_comp} shows the MSE comparison of the {\sf FCSS} vs. BPDN, where the MSE is defined as $10\log_{10}\frac{1}{T}\sum_{t=1}^T \|\widehat{\mathbf{x}}_t-\mathbf{x}_t\|_2^2$. The {\sf FCSS} algorithm significantly outperforms BPDN, especially at high SNR values. Figure \ref{fig:tv_dcs_vs_ccs} compares the performance of {\sf FCSS} and BPDN on a sample component at a compression level of $n/p = 2/3$, in order to visualize the performance gain implied by Figure \ref{fig:tv_dcs_vs_ccs}.

Finally, Figure \ref{fig:imagesc_sim} shows the comparison of the estimated states for the entire simulated data in the denoising regime. As can be observed from the figure, the sparsity pattern of the states and innovations are captured while significantly denoising the observed states.

\vspace{-3mm}
\subsection{Application to Calcium Signal Deconvolution}

Calcium imaging takes advantage of intracellular calcium flux to directly visualize calcium signaling in living neurons. This is done by using calcium indicators, which are fluorescent molecules that can respond to the binding of calcium ions by changing their fluorescence properties and using a fluorescence or two-photon microscope and a CCD camera to capture the visual patterns \cite{smetters1999detecting, stosiek2003vivo}. Since spikes are believed to be the units of neuronal computation, inferring spiking activity from calcium recordings, referred to as calcium deconvolution, is an important problem in neural data analysis. Several approaches to calcium deconvolution have been proposed in the neuroscience literature, including model-free approaches such as sequential Monte Carlo methods \cite{vogelstein2009spike} and model-based approaches such as non-negative deconvolution methods \cite{vogelstein2010fast, pnevmatikakis2016simultaneous}. These approaches require solving convex optimization problems, which do not scale well with the temporal dimension of the data. In addition, they lack theoretical performance guarantees and do not provide clear measures for assessing the statistical significance of the detected spikes.

In order to construct confidence bounds for our estimates, we employ recent results from high-dimensional statistics \cite{van2014asymptotically}. We first compute the confidence intervals around the outputs of the {\sf FCSS} estimates using the node-wise regression procedure of \cite{van2014asymptotically}, at a confidence level of $1-\frac{\alpha}{2}$. We perform the node-wise regression separately for each time $t$. For an estimate $\widehat{\mathbf{x}}_t$, we obtain $\widehat{\mathbf{x}}_t^{\sf u}$ and $\widehat{\mathbf{x}}_t^{\sf l}$ as the upper and lower confidence bounds, respectively. Next, we partition the estimates into small segments, starting with a local minimum (trough) and ending in a local maximum (peak). For the $i^{\sf th}$ component of the estimate, let $t_{\sf min}$ and $t_{\sf max}$ denote the time index corresponding to two such consecutive troughs and peaks. If the difference $(\widehat{\mathbf{x}}_{t_{\sf max}}^{\sf l})_i - (\widehat{\mathbf{x}}_{t_{\sf min}}^{\sf u})_i$ is positive, the detected innovation component is declared significant (i.e., spike) at a confidence level of $1- \alpha$, otherwise it is discarded (i.e., no spike). We refer to this procedure as Pruned-{\sf FCSS} ({\sf PFCSS}).

\begin{figure*}[t!]
\centering
\includegraphics[width=1\textwidth]{figures/tv_cal_simul_ss_vs_foopsi}
\caption{Ground-truth performance comparison between {\sf PFCSS} and constrained f-oopsi. Top: the observed calcium traces (black) and ground-truth electrophysiology data (blue), Middle: {\sf PFCSS} state estimates (black) with $95\%$ confidence intervals (orange) and the detected spikes (blue), Bottom: the constrained f-oopsi state estimates (black) and the detected spikes (blue). The {\sf FCSS} spike estimates closely match the ground-truth spikes with only a few false detections, while the constrained f-oopsi estimates contain significant clustered false detections.}
\label{fig:tv_cal_simul_ss_vs_foopsi}
\end{figure*}

\begin{figure*}[htb!]
\centering
\includegraphics[width=1\textwidth]{figures/tv_3neuron-01}
\caption{Performance of {\sf FCSS} on large-scale calcium imaging data. Top, middle and bottom rows correspond to three selected neurons labeled as Neuron 1, 2, and 3, respectively. (a) raw calcium traces, (b) {\sf FCSS} reconstruction with no compression, (c) {\sf FCSS} reconstruction with $2/3$ compression ratio. Orange hulls show $95\%$ confidence intervals. The {\sf FCSS} significantly denoises the observed traces in both the uncompressed and compressed settings.\label{fig:tv_cal}}
\end{figure*}

\begin{figure*}[htb!]
\centering
\includegraphics[width=1\textwidth]{figures/spikes_3neuron-01}
\caption{Reconstructed spikes of {\sf PFCSS} and constrained f-oopsi from large-scale calcium imaging data. Top, middle and bottom rows correspond to three selected neurons labeled as Neuron 1, 2, and 3, respectively. (a) constrained f-oopsi spike estimates, (b) {\sf PFCSS} spike estimates with no compression, (c) {\sf PFCSS} spike estimates with $2/3$ compression ratio. The {\sf PFCSS} estimates in both the uncompressed and compressed settings are sparse in time, whereas the constrained f-oopsi estimates are in the form of clustered spikes. \label{fig:tv_cal_spikes}}
\end{figure*}

\begin{figure}[htb!]
\centering
\includegraphics[width=1\columnwidth]{figures/tv_raster-01}
\caption{Raster plot of the estimated spikes from large-scale calcium imaging data. Left: the constrained f-oopsi estimates, Middle: {\sf PFCSS} estimates with no compression, Right: {\sf PFCSS} estimates with a $\frac{1}{3}$ compression ratio. The {\sf PFCSS} estimates are spatiotemporally sparse, whereas the constrained f-oopsi outputs temporally clustered spike estimates.}
\label{fig:tv_raster_cal}
\end{figure}

We first apply the {\sf FCSS} algorithm for calcium deconvolution in a scenario where the ground-truth spiking is recorded \emph{in vitro} through simultaneous electrophysiology (cell-attached patch clamp) and two-photon calcium imaging. The calcium trace as well as the ground-truth spikes are shown for a sample neuron in Figure \ref{fig:tv_cal_simul_ss_vs_foopsi}--(a). The {\sf FCSS} denoised estimate of the states (black) and the detected spikes (blue) using $95\%$ confidence intervals (orange hulls) and the corresponding quantities for the constrained f-oopsi algorithm \cite{pnevmatikakis2016simultaneous} are shown in Figures \ref{fig:tv_cal_simul_ss_vs_foopsi}--(b) and --(c), respectively. Both algorithms detect the large dynamic changes in the data, corresponding to the spikes, which can also be visually captured in this case. However, in doing so, the f-oopsi algorithm incurs a high rate of false positive errors, manifested as clustered spikes around the ground truth events. Similar to f-oopsi, most state-of-the-art model-based methods suffer from high false positive rate, which makes the inferred spike estimates unreliable. Thanks to the aforementioned pruning process based on the confidence bounds, the {\sf PFCSS} is capable of rejecting the insignificant innovations, and hence achieve a lower false positive rate. One factor responsible for this performance gap can be attributed to the underestimation of the calcium decay rate in the transition matrix estimation step of f-oopsi. However, we believe the performance gain achieved by {\sf FCSS} is mainly due to the explicit modeling of the sparse nature of the spiking activity by going beyond the Gaussian state-space modeling paradigm.

Next, we apply the {\sf FCSS} algorithm to large-scale \emph{in vivo} calcium imaging recordings, for which the ground-truth is not available due to measurement constraints. The data used in our analysis was recorded from $219$ spontaneously active neurons in mouse auditory cortex. The two-photon microscope operates at a rate of $30$ frames per second. We chose $T=2000$ samples corresponding to $1$ minute for the analysis. We chose $p=108$ well-separated neurons visually. We estimate the measurement noise variance by appropriate re-scaling of the power spectral density in the high frequency bands where the signal is absent. We chose a value of $\epsilon = 10^{-10}$. It is important to note that estimation of the measurement noise variance is critical, since it affects the width of the confidence intervals and hence the detected spikes. Moreover, we estimate the baseline fluorescence by averaging the signal over values within a factor of $3$ standard deviations of the noise. By inspecting Eq. (\ref{eq:tv_prob_def_dual}), one can see a trade-off between the choice of $\lambda$ and the estimate of the observation noise variance $\sigma^2$. We have done our analysis in both the compression regime, with a compression ratio of $1/3$ ($n/p = 2/3$), and the denoising regime. The measurements in the compressive regime were obtained from applying i.i.d. Gaussian random matrices to the observed calcium traces. The latter is done to motivate the use of compressive imaging, as opposed to full sampling of the field of view.

Figure \ref{fig:tv_cal}--(a) shows the observed traces for four selected neurons. The reconstructed states using {\sf FCSS} in the compressive and denoising regimes are shown in Figures \ref{fig:tv_cal}--(b) and --(c), respectively. The $90\%$ confidence bounds are shown as orange hulls. The {\sf FCSS} state estimates are significantly denoised while preserving the calcium dynamics. Figure \ref{fig:tv_cal_spikes} shows the detected spikes using constrained f-oopsi and {\sf PFCSS} in both the compressive and denoising regimes. Finally, Figure \ref{fig:tv_raster_cal} shows the corresponding raster plots of the reconstructed spikes for the entire ensemble of neurons. Similar to the preceding application on ground-truth date, the f-oopsi algorithm detects clusters of spikes, whereas the {\sf PFCSS} procedure results in sparser spike detection. This results in the detection of seemingly more active neurons in the raster plot. However, motivated by the foregoing ground-truth analysis, we believe that a large fraction of these detected spikes may be due to false positive errors. Strikingly, even with a compression ratio of $1/3$ the performance of the {\sf PFCSS} is similar to the denoising case. The latter observation corroborates the feasibility of compressed two-photon imaging, in which only a random fraction of the field of view is imaged, which in turn can result in higher acquisition rates.

In addition to the foregoing discussion on the comparisons in Figures \ref{fig:tv_cal_simul_ss_vs_foopsi}, \ref{fig:tv_cal}, \ref{fig:tv_cal_spikes}, and \ref{fig:tv_raster_cal}, two remarks are in order. First, the iterative solution at the core of {\sf FCSS} is linear in the observation length and hence significantly faster than the batch-mode optimization procedure used for constrained f-oopsi. Our comparisons suggest that the {\sf FCSS} reconstruction is at least 3 times faster than f-oopsi for moderate data sizes of the order of tens of minutes. Moreover, the vector formulation of {\sf FCSS} allows for easy parallelization {(without the need for GPU implementations)}, which allows simultaneous processing of ROI's without losing speed. {As a numerical example the results of Figure \ref{fig:tv_cal}--(b) took an average of 60-70 seconds to calculate for all $p=108$ ROI's and $T=2000$ frames on an Apple Macintosh desktop computer}. Second, using only about two-thirds of the measurements achieves similar results by {\sf FCSS} as using the full measurements.

\vspace{-2mm}
\subsection{Application to Sleep Spindle Detection}
In this section we use compressible state-space models in order to model and detect sleep spindles. A sleep spindle is a burst of oscillatory brain activity manifested in the EEG that occurs during stage 2 non-rapid eye movement (NREM) sleep. It consists of stereotypical $12$--$14~\text{Hz}$ wave packets that last for at least $0.5$ seconds \cite{de2003sleep}. The spindles occur with a rate of $2$--$5\%$ in time, which makes their generation an appropriate candidate for compressible dynamics. Therefore, we hypothesize that the spindles can be modeled using a combination of few echoes of the response of a second order compressible state-space model. As a result, the spindles can be decomposed as sums of modulated sine waves.

\begin{figure}[b!]
\centering
\includegraphics[width=1\columnwidth]{figures/tv_spindle_sim-01}
\caption{Performance of {\sf FCSS} on simulated spindles. Top: simulated clean data (black) and ground-truth spindle events (red), Middle: simulated noisy data, Bottom: the denoised signal (black) and deconvolved spindle events (red). The {\sf FCSS} estimates are significantly denoised and closely match the ground-truth data shown in the top panel.}
\label{fig:tv_spindle_sim}

\end{figure}

\begin{figure*}[t!]
\centering
\includegraphics[width=1\textwidth]{figures/tv_spindle_dream-01}
\caption{Performance comparison between {\sf FCSS} and band-pass filtered EEG data. Left and right panels correspond to two selected electrodes labeled as 1 and 2, respectively. The orange blocks show the extent of the detected spindles by the expert. Top: raw EEG data, Middle: band-pass filtered EEG data in the $12$--$14~\text{Hz}$ band, Bottom: {\sf FCSS} spindle estimates. The {\sf FCSS} estimates closely match the expert annotations, while the band-pass filtered data contains significant signal components outside of the orange blocks.}
\label{fig:tv_spindle_dream}
\vspace{-3mm}
\end{figure*}

\begin{figure}[t!]
\centering
\includegraphics[width=0.5\columnwidth]{figures/tv_spindle_ROC}
\caption{The ROC curves for {\sf FCSS} (solid red) and bandpass-filtered RMS (dashed black). The {\sf FCSS} outperforms the widely-used band-pass filtered RMS method as indicated by the ROC curves.}
\label{fig:tv_spindle_roc}
\end{figure}

In order to model the oscillatory nature of the spindles, we consider a second order autoregressive (AR) model where the pole locations are given by $a e^{-j2\pi \frac{f}{f_s}}$ and $a e^{+j2\pi \frac{f}{f_s}}$, where $0<a<1$ is a positive constant controlling the effective duration of the impulse response, $f_s$ is the sampling frequency and $f$ is a putative frequency accounting for the dominant spindle frequency. The equivalent state-space model for which the MAP estimation admits the {\sf FCSS} solution is therefore:

\begin{equation}
\label{eq:tv_lap_state_space_spindle}
\begin{array}{ll}
\mathbf{x}_t = 2a\cos\left(\textstyle 2\pi \frac{f}{f_s}\right) \mathbf{x}_{t-1}-a^2 \mathbf{x}_{t-2}+ \mathbf{w}_t,\\
\mathbf{y}_t = \mathbf{A}_t \mathbf{x}_t + \mathbf{v}_t.
\end{array}
\vspace{-1mm}
\end{equation}

Note that impulse response of the state-space dynamics is given by $h_n  = a^n \cos\left(\textstyle 2\pi \frac{f}{f_s} n\right)u_n$, which is a stereotypical decaying sinusoid. By defining the augmented state $\widetilde{\mathbf{x}}_t = [\mathbf{x}_t', \mathbf{x}_{t-1}']'$,
Eq. (\ref{eq:tv_lap_state_space_spindle}) can be expressed in the canonical form:
\begin{equation}
\label{eq:tv_lap_state_space_spindle_canonical}
\begin{array}{l}
\widetilde{\mathbf{x}}_t = \widetilde{\mathbf{\Theta}} \widetilde{\mathbf{x}}_{t-1}+ \widetilde{\mathbf{w}}_t, \qquad \mathbf{y}_t = \widetilde{\mathbf{A}}_t \widetilde{\mathbf{x}}_t + \mathbf{v}_t,
\end{array}
\end{equation}
where $\widetilde{\mathbf{w}}_t := [\mathbf{w}'_t
, \mathbf{0}']'$ , $\widetilde{\mathbf{A}}_t = \begin{array}{c:c}
[\mathbf{A}_t & \mathbf{0}]
\end{array}$ and
\[\widetilde{\mathbf{\Theta}} = \left[\begin{array}{c:c}
2a\cos\left(\textstyle 2\pi \frac{f}{f_s}\right)\mathbf{I} & -a^2 \mathbf{I}\\
\hdashline
\mathbf{I} & \mathbf{0}
\end{array}
\right].
\]

Eq. (\ref{eq:tv_Q_calculated}) can be used to update $\widetilde{\mathbf{\Theta}}$ in the M step of the {\sf FCSS} algorithm. However, $\widetilde{\mathbf{\Theta}}$ has a specific structure in this case, determined by $a$ and $f$, which needs to be taken into account in the optimization step. Let $\phi =: 2a\cos\left(2\pi \frac{f}{f_s}\right)$ and $\psi = a^2$, and let
\[
\sum_{t=1}^T   \widetilde{\mathbf{W}}_t^{(l)} \left ({\widetilde{\mathbf{x}}}^{(l, m+1)}_{t-1|T}{\operatorname{\widetilde{\mathbf{x}}}'}^{(l, m+1)}_{t-1|T} + \widetilde{\mathbf{\Sigma}}_{t-1|T}^{(l,m+1)} \right ) =:  \left [ \begin{tabular}{c:c} $\mathbf{A}$ & $\mathbf{B}$\\\hdashline$\mathbf{C}$ & $\mathbf{D}$ \end{tabular} \right],
\]
\begin{align}
\nonumber & \displaystyle \sum_{t=1}^T \widetilde{\mathbf{W}}_t^{(l)} \left(\widetilde{\mathbf{x}}^{(l, m+1)}_{t-1|T} {\operatorname{\widetilde{\mathbf{x}}}'}^{(l, m+1)}_{t|T}+ \widetilde{\mathbf{x}}^{(l, m+1)}_{t|T} {\operatorname{\widetilde{\mathbf{x}}}'}^{(l, m+1)}_{t-1|T} + 2 \widetilde{\mathbf{\Sigma}}_{t-1,t|T}^{(l,m+1)}\right) =: \left [ \begin{tabular}{c:c} $\mathbf{E}$ & $\mathbf{F}$\\\hdashline$\mathbf{G}$ & $\mathbf{H}$ \end{tabular} \right]\!.
\end{align}
\vspace{-2mm}

\noindent Then, Eq. (\ref{eq:tv_Q_calculated}) is equivalent to maximizing
\begin{align}
\label{eq:tv_estimate_af}
\maximize_{\phi,\psi} \ \ & \frac{\lambda}{2} (\phi^2 \text{Tr}(\mathbf{A}) \!-\! \phi \psi \text{Tr} (\mathbf{B}+\mathbf{C}) + \psi^2 \text{Tr} (\mathbf{D}) + \text{Tr}(\mathbf{A})) - {\textstyle \frac{\lambda}{2}} (\phi \text{Tr}(\mathbf{E}) - \psi \text{Tr}(\mathbf{G}) + \text{Tr}(\mathbf{F})).
\end{align}
subject to $0 \le \psi \le 1$ and $\phi^2 \le 4 \psi$, which can be solved using interior point methods. In our implementation, we have imposed additional priors of $f \sim {\sf Uniform} (12, 14)$ Hz and $a \sim {\sf Uniform} (0.95,0.99)$, which simplifies the constraints on $\phi$ and $\psi$ to
\[
0.95^2 \leq \psi \leq 0.99^2, \ \ \text{and} \ \  4\psi \cos^2\left(\textstyle 2\pi \frac{12}{f_s}\right)\leq \phi^2 \leq 4\psi \cos^2\left(\textstyle 2\pi \frac{14}{f_s}\right).
\]
 Given the convexity of the cost function in (\ref{eq:tv_estimate_af}), one can conclude from the KKT conditions that if the global minimum is not achieved inside the region of interest it must be achieved on the boundaries.

Figure \ref{fig:tv_spindle_sim}, top panel, shows two instances of simulated spindles (black traces) with parameters $f_s = 200$ Hz, $f = 13$ Hz and $a = 0.95$, with the ground truth events generating the wave packets shown in red. The middle panel shows the noisy version of the data with an SNR of $-7.5$ dB. The noise was chosen as white Gaussian noise plus slowly varying (2 Hz) oscillations to resemble the slow oscillations in real EEG data. As can be observed the simulated signal exhibits visual resemblance to real spindles, which verifies the hypothesis that spindles can be decomposed into few combinations of wave packets generated by a second order AR model. The third panel, shows the denoised data using {\sf FCSS}, which not only is successful in detecting the the ground-truth dynamics (red bars), but also significantly denoises the data.

We next apply {\sf FCSS} to real EEG recordings from stage-2 NREM sleep. Manual scoring of sleep spindles can be very time-consuming, and achieving accurate manual scoring on a long-term recording is a highly demanding task with the associated risk of decreased diagnosis. Although automatic spindle detection would be attractive, most available algorithms sensitive to variations in spindle amplitude and frequency that occur between both subjects and derivations, reducing their effectiveness \cite{nonclercq2013sleep, schimicek1994automatic}. Moreover most of these algorithms require significant pre- and post-processing and manual tuning. Examples include algorithms based on Empirical Mode Decomposition (EMD) \cite{yang2006detection, fraiwan2012automated, causa2010automated}, data-driven Bayesian methods \cite{babadi2012diba}, and machine learning approaches \cite{duman2009efficient, ventouras2005sleep}. Unlike our approach, none of the existing methods consider modeling the generative dynamics of spindles, as transient sparse events in time, in the detection procedure.

The data used in our analysis is part of the recordings in the DREAMS project \cite{devuyst2011automatic}, recorded using a 32-channel polysomnograpgh. We have used the EEG channels in our analysis. The data was recorded at a rate of $f_s=200$ Hz for 30 minutes. The data was scored for sleep spindles independently by two experts. We have used expert annotations to separate regions which include spindles for visualization purposes. For comparison purposes, we use a bandpass filtered version of the data within $12$--$14~\text{Hz}$, which is the basis of several spindle detection algorithms \cite{schimicek1994automatic,clemens2005overnight, huupponen2007development}, hallmarked by the widely-used bandpass filtered root-mean-square (RMS) method \cite{schimicek1994automatic}.

Figure \ref{fig:tv_spindle_dream} shows the detection results along with the bandpass filtered version of the data for two of the EEG channels. The red bars show the expert markings of the onset and offset of the spindle events. The {\sf FCSS} simultaneously captures the spindle events and suppressed the activity elsewhere, whereas the bandpass filtered data produces significant activity in the $12$--$14~\text{Hz}$ throughout the observation window, resulting in high false positives. To quantify this observation, we have computed the ROC curves of the {\sf FCSS} and bandpass filtering followed by root mean square (RMS) computation in Figure \ref{fig:tv_spindle_roc}, which confirms the superior performance of the {\sf FCSS} algorithm over the data set. The annotations of one of the experts has been used for as the ground truth benchmark.

\vspace{-2mm}
\section{Discussion}
\label{sec:disc}
In this section, we discuss the implication of our techniques in regard to the application domains as well as existing methods.

\vspace{-3mm}
\subsection{Connection to existing literature in sparse estimation}
Contrary to the traditional compressive sensing, our linear measurement operator does not satisfy the RIP \cite{ba2012exact}, despite the fact that $\mathbf{A}_t$'s satisfy the RIP. Nevertheless, we have extended the near-optimal recovery guarantees of CS to our compressible state-space estimation problem via Theorem \ref{thm:tv_main}. Closely related problems to our setup are the super-resolution and sparse spike deconvolution problems \cite{candes2014towards,duval2015exact}, in which abrupt changes with minimum separation in time are resolved in fine scales using coarse (lowpass filtered) frequency information, which is akin to working in the compressive regime.

Theoretical guarantees of CS require the number of measurements to be roughly proportional to the sparsity level for stable recovery \cite{Negahban}. These results do not readily generalize to the cases where the sparsity lies in the dynamics, not the states per se. Most of the dynamic compressive sensing techniques such as Kalman filtered compressed sensing, assume partial information about the support or estimate them in a greedy and often ad-hoc fashion \cite{vaswani2010ls, vaswani2008kalman, carmi2010methods, ziniel2013dynamic, zhan2015time}. As another example, the problem of recovering discrete signals which are approximately sparse in their gradients using compressive measurements, has been studied in the literature using Total Variation (TV) minimization techniques \cite{needell2013stable,poon2015role}. For one-dimensional signals, since the gradient operator is not orthonormal, the Frobenius operator norm of its inverse grows linearly with the discretization level \cite{needell2013stable}. Therefore, stability results of TV minimization scale poorly with respect to discretization level. In higher dimensions, however, the fast decay of Haar coefficients allow for near-optimal theoretical guarantees \cite{cohen1999nonlinear}. A major difference of our setup with those of CS for TV-minimization is the \emph{structured} and \emph{causal} measurements, which unlike the non-causal measurements in \cite{needell2013stable}, do not result in an overall measurement matrix satisfying RIP. We have considered dynamics with convergent transition matrices, in order to generalize the TV minimization approach. To this end, we showed that using the state-space dynamics one can infer temporally global information from local and causal measurements. Another closely related problem is the fused lasso \cite{tibshirani2005sparsity} in which sparsity is promoted both on the covariates and their differences.

\vspace{-3mm}
\subsection{Application to calcium deconvolution} In addition to scalability and the ability to detect abrupt transitions in the states governed by discrete events in time (i.e., spikes), our method provides several other benefits compared to other spike deconvolution methods based on state-space models, such as the constrained f-oopsi algorithm. First, our sampling-complexity trade-offs are known to be optimal from the theory of compressive sensing, whereas no performance guarantee exists for constrained f-oopsi. Second, we are able to construct precise confidence intervals on the estimated states, whereas constrained f-oopsi does not produce confidence intervals over the detected spikes. A direct consequence of these confidence intervals is estimation of spikes with high fidelity and low false alarm. Third, our comparisons suggest that the {\sf FCSS} reconstruction is at least 3 times faster than f-oopsi for moderate data sizes of the order of tens of minutes. Finally, our results corroborate the possibility of using compressive measurement for reconstruction and denoising of calcium traces. From a practical point of view, a compressive calcium imaging setup can lead to higher scanning rate as well as better reconstructions, which allows monitoring of larger neuronal populations \cite{pnevmatikakis2013sparse}. Due to the structured nature of our sampling and reconstruction schemes, we can avoid prohibitive storage problems and benefit from parallel implementations.

\vspace{-3mm}
\subsection{Application to sleep spindle detection}
Another novel application of our modeling and estimating framework is to case sleep spindle generation as a second-order dynamical system governed by compressive innovations, for which {\sf FCSS} can be efficiently used to denoise and detect the spindle events. Our modeling framework suggest that spectrotemporal spindle dynamics cannot be fully captured by just pure sinusoids via bandpass filtering, as the data consistently contains significant $12$--$14~\text{Hz}$ oscillations almost everywhere (See Figure \ref{fig:tv_spindle_dream}). Therefore, using the bandpass filtered data for further analysis purposes clearly degrades the performance of the resulting spindle detection and scoring algorithms. The {\sf FCSS} provides a robust alternative to bandpass filtering in the form of model-based denoising.

In contrast to state-of-the-art methods for spindle detection, our spindle detection procedure requires minimal pre- and post-processing steps. We expect similar properties for higher order AR dynamics, which form a useful generalization of our methods for deconvolution of other transient neural signals. In particular, K-complexes during the stage 2 NREM sleep form another class of transient signals with high variability. A potential generalization of our method using higher order models can be developed for simultaneous detection of K-complexes and spindles.

\vspace{-2mm}
\section{Concluding Remarks}\label{sec:tv_conc}

In this chapter, we considered estimation of compressible state-space models, where the state innovations consist of compressible discrete events. For dynamics with convergent state transition dynamics, using theory of compressed sensing we provided an optimal error bound and stability guarantees for the dynamic $\ell_1$-regularization algorithm which is akin to the MAP estimator for a Laplace state-space model. We also developed a fast and low-complexity iterative algorithm, namely {\sf FCSS}, for estimation of the states as well as their transition matrix. We further verified the validity of our theoretical results through simulation studies as well as application to spike deconvolution from calcium traces and detection of sleep spindles from EEG data. Our methodology has two unique major advantages: first, we have proven theoretically why our algorithm performs well, and characterized its error performance. Second, we have developed a fast algorithm, with guaranteed convergence to a solution of the deconvolution problem, which for instance, is $\sim 3$ times faster than the widely-used f-oopsi algorithm in calcium deconvolution applications.

While we focused on two specific application domains, our modeling and estimation techniques can be generalized to apply to broader classes of signal deconvolution problems: we have provided a framework to model transient phenomena which are driven by sparse generators in time domain, and whose event onsets are of importance. Examples include heart beat dynamics and rapid changes in the covariance structure of neural data (e.g., epileptic seizures). In the spirit of easing reproducibility, we have made a MATLAB implementation of our algorithm publicly available \cite{code}.

config: arxiv