Exact Spike Train Inference Via $\ell_0$ Optimization

TL;DR

This paper introduces an efficient dynamic programming method for exact spike train inference from fluorescence data by replacing an $\\ell_1$ penalty with an $\\ell_0$ penalty, significantly improving accuracy and computational efficiency.

Contribution

The paper demonstrates that an $\\ell_0$ optimization approach for spike inference is computationally feasible and more accurate than previous $\\ell_1$ methods, with practical software implementation.

Findings

The $\\ell_0$ optimization problem can be solved efficiently with dynamic programming.

The proposed method outperforms previous $\\ell_1$ based approaches in simulations and real data.

The algorithm runs in minutes on large datasets with 100,000 timesteps.

Abstract

In recent years, new technologies in neuroscience have made it possible to measure the activities of large numbers of neurons simultaneously in behaving animals. For each neuron, a fluorescence trace is measured; this can be seen as a first-order approximation of the neuron's activity over time. Determining the exact time at which a neuron spikes on the basis of its fluorescence trace is an important open problem in the field of computational neuroscience. Recently, a convex optimization problem involving an $\ell_1$ penalty was proposed for this task. In this paper, we slightly modify that recent proposal by replacing the $\ell_1$ penalty with an $\ell_0$ penalty. In stark contrast to the conventional wisdom that $\ell_0$ optimization problems are computationally intractable, we show that the resulting optimization problem can be efficiently solved for the global optimum using an extremely simple and efficient dynamic programming algorithm. Our R-language implementation of the proposed algorithm runs in a few minutes on fluorescence traces of $100,000$ timesteps. Furthermore, our proposal leads to substantial improvements over the previous $\ell_1$ proposal, in simulations as well as on two calcium imaging data sets. R-language software for our proposal is available on CRAN in the package LZeroSpikeInference. Instructions for running this software in python can be found at https://github.com/jewellsean/LZeroSpikeInference.

Equations

$$y_t={\beta }_0+{\beta }_1{c}_t+{ϵ}_t,{ϵ}_t{\sim }_{\text{ind.}}(0,{\sigma }^2),t=1,\dots ,T;$$

$$c_t=\sum _{i=1}^p{\gamma }_i{c}_{t-i}+s_t,t=p+1,\dots ,T.$$

$$\underset{c_1,\dots ,c_T,s_2,\dots ,s_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda \sum _{t=2}^T|s_t|\right\}\text{mbox}subjectto{s}_t=c_t-\gamma c_{t-1}\ge 0,$$

$$\underset{c_1,\dots ,c_T,s_2,\dots ,s_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda |c_1|+\lambda \sum _{t=2}^T|s_t|\right\}\text{mbox}subjectto{s}_t=c_t-\gamma c_{t-1}\ge 0.$$

$$\underset{c_1,\dots ,c_T,s_2,\dots ,s_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda \sum _{t=2}^T{1}_{\left(s_t\ne 0\right)}\right\}\text{mbox}subjectto{s}_t=c_t-\gamma c_{t-1}\ge 0,$$

$$\underset{c_1,\dots ,c_T,s_2,\dots ,s_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda \sum _{t=2}^T{1}_{\left(s_t\ne 0\right)}\right\}\text{mbox}subjectto{s}_t=c_t-\gamma c_{t-1}.$$

$$\underset{c_1,\dots ,c_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda \sum _{t=2}^T{1}_{\left(c_t-\gamma c_{t-1}\ne 0\right)}\right\}.$$

$$\underset{0={\tau }_0<{\tau }_1<\dots <{\tau }_k<{\tau }_{k+1}=T,k}{\mathrm{minimize}}\left\{\sum _{j=0}^k\mathcal{D}(y_{({\tau }_j+1):{\tau }_{j+1}})+\lambda k\right\},$$

$$\mathcal{D}(y_{a:b})\equiv \underset{c_a,c_t=\gamma c_{t-1},t=a+1,\dots ,b}{\min }\left\{\frac{1}{2}\sum _{t=a}^b{\left(y_t-c_t\right)}^2\right\}.$$

$${\hat{c}}_t$$

$$\mathcal{D}(y_{a:b})$$

$$\mathcal{C}(y_{a:b})$$

$$F(s)$$

$$=\underset{0={\tau }_0<{\tau }_1<\cdots <{\tau }_k<{\tau }_{k+1}=s,k}{\min }\left\{\sum _{j=0}^k\left[\mathcal{D}(y_{({\tau }_j+1):{\tau }_{j+1}})+\lambda \right]-\lambda \right\}$$

$$=\underset{0={\tau }_0<{\tau }_1<\cdots <{\tau }_k<{\tau }_{k+1}=s,k}{\min }\left\{\sum _{j=0}^{k-1}\left[\mathcal{D}(y_{({\tau }_j+1):{\tau }_{j+1}})+\lambda \right]-\lambda +\mathcal{D}(y_{({\tau }_k+1):{\tau }_{k+1}})+\lambda \right\}$$

$$=\underset{0\le {\tau }_k<{\tau }_{k+1}=s}{\min }\left\{\underset{0={\tau }_0<{\tau }_1<\cdots <{\tau }_{k^{\prime }}<{\tau }_{k^{\prime }+1}={\tau }_k,k^{\prime }}{\min }\left\{\sum _{j=0}^{k^{\prime }}\left[\mathcal{D}(y_{({\tau }_j+1):{\tau }_{j+1}})+\lambda \right]-\lambda \right\}+\mathcal{D}(y_{({\tau }_k+1):{\tau }_{k+1}})+\lambda \right\}$$

$$=\underset{0\le \tau <s}{\min }\left\{F(\tau )+\mathcal{D}(y_{(\tau +1):s})+\lambda \right\}.$$

$$F(1)=\underset{0\le \tau <1}{\min }\left\{F(\tau )+\mathcal{D}(y_{(\tau +1):1})+\lambda \right\}=F(0)+\mathcal{D}(y_1)+\lambda ,$$

$$F(2)=\underset{\tau \in \{0,1\}}{\min }\left\{F(\tau )+\mathcal{D}(y_{(\tau +1):2})+\lambda \right\}=\min \left\{F(0)+\mathcal{D}(y_{1:2})+\lambda ,F(1)+\mathcal{D}(y_2)+\lambda \right\}.$$

$${\hat{c}}_t$$

$${\tilde{s}}_t={\hat{s}}_t{1}_{({\hat{s}}_t\ge L)},t=2,\dots ,T,$$

$$\mathrm{MSE}(c,\hat{c})=\frac{1}{T}\sum _{t=1}^T(c_t-{\hat{c}}_t{)}^2.$$

$$\underset{c_1,\dots ,c_T,{\beta }_{01},\dots ,{\beta }_{0T}}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T(y_t-c_t-{\beta }_{0t}{)}^2+\lambda \sum _{t=2}^T{1}_{\left(c_t\ne \gamma c_{t-1},{\beta }_{0t}\ne {\beta }_{0,t-1}\right)}\right\},$$

$$\mathcal{D}(y_{a:b})$$

$$\underset{c_1,\dots ,c_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T(y_t-c_t{)}^2+\lambda \sum _{t=p+1}^T{1}_{\left(c_t\ne {\sum }_{i=1}^p{\gamma }_i{c}_{t-i}\right)}\right\}.$$

$$\mathcal{D}(y_{a:b})$$

$$\underset{c_1,\dots ,c_T}{\mathrm{minimize}}\left\{\frac{1}{2}\sum _{t=1}^T{\left(y_t-c_t\right)}^2+\lambda \sum _{t=2}^T|c_t-c_{t-1}|\right\},$$

$$\mathcal{D}(y_{a:b})\equiv \underset{c_a,c_t=\gamma c_{t-1},t=a+1,\dots ,b}{\min }\left\{\frac{1}{2}\sum _{t=a}^b{\left(y_t-c_t\right)}^2\right\}$$

$$\mathcal{D}(y_{a:b})=\underset{c_a}{\min }\left\{\frac{1}{2}\sum _{t=a}^b{\left(y_t-{\gamma }^{t-a}{c}_a\right)}^2\right\}.$$

$${\hat{c}}_a=\frac{{\sum }_{t=a}^b{y}_t{\gamma }^{t-a}}{{\sum }_{t=a}^b{\gamma }^{2(t-a)}},$$

$$\mathcal{D}(y_{a:b})=\frac{1}{2}\sum _{t=a}^b{\left(y_t-{\gamma }^{t-a}{\hat{c}}_a\right)}^2,$$

$$\hat{\gamma }=\underset{\gamma }{\mathrm{argmin}}\left\{\mathcal{D}\left(y_{a:b}\right)\right\}=\underset{\gamma }{\mathrm{argmin}}\left\{\underset{c_a,c_t=\gamma c_{t-1},t=a+1,\dots ,b}{\min }\left\{\frac{1}{2}\sum _{t=a}^b(y_t-c_t{)}^2\right\}\right\}.$$

Full text

Recently, a convex optimization problem involving an $\ell_{1}$ penalty was proposed for this task. In this paper, we slightly modify that recent proposal by replacing the $\ell_{1}$ penalty with an $\ell_{0}$ penalty. In stark contrast to the conventional wisdom that $\ell_{0}$ optimization problems are computationally intractable, we show that the resulting optimization problem can be efficiently solved for the global optimum using an extremely simple and efficient dynamic programming algorithm. Our R-language implementation of the proposed algorithm runs in a few minutes on fluorescence traces of $100,000$ timesteps. Furthermore, our proposal leads to substantial improvements over the previous $\ell_{1}$ proposal, in simulations as well as on two calcium imaging data sets.

R*-language software for our proposal is available on CRAN in the package

LZeroSpikeInference. Instructions for running this software in python can be found at https://github.com/jewellsean/LZeroSpikeInference.*

*Department of Statistics

University of Washington

Seattle, Washington 98195

\safe@setrefe1@*

*Departments of Biostatistics and Statistics

University of Washington

Seattle, Washington 98195

\safe@setrefe2@*

@afterheading

1 Introduction

When a neuron spikes, calcium floods the cell. In order to quantify intracellular calcium levels, calcium imaging techniques make use of fluorescent calcium indicator molecules (Dombeck et al., 2007; Ahrens et al., 2013; Prevedel et al., 2014). Thus, a neuron’s fluorescence trace can be seen as a first-order approximation of its activity level over time.

However, the fluorescence trace itself is typically not of primary scientific interest: instead, it is of interest to determine the underlying neural activity, that is, the specific times at which the neuron spiked. Inferring the spike times on the basis of a fluorescence trace amounts to a challenging deconvolution problem, which has been the focus of substantial investigation (Grewe et al., 2010; Pnevmatikakis et al., 2013; Theis et al., 2016; Deneux et al., 2016; Sasaki et al., 2008; Vogelstein et al., 2009; Yaksi and Friedrich, 2006; Vogelstein et al., 2010; Holekamp, Turaga and Holy, 2008; Friedrich and Paninski, 2016; Friedrich, Zhou and Paninski, 2017). In this paper, we propose a new approach for this task, which is based upon the following insight: an auto-regressive model for calcium dynamics that has been extensively studied in the neuroscience literature (Vogelstein et al., 2010; Friedrich and Paninski, 2016; Friedrich, Zhou and Paninski, 2017) leads directly to a simple $\ell_{0}$ optimization problem, for which an efficient and exact algorithm is available.

1.1 An Auto-Regressive Model for Calcium Dynamics

In this paper, we will revisit an auto-regressive model for calcium dynamics that has been considered by a number of authors in the recent literature (Vogelstein et al., 2010; Pnevmatikakis et al., 2016; Friedrich and Paninski, 2016; Friedrich, Zhou and Paninski, 2017). We closely follow the notation of Friedrich, Zhou and Paninski (2017). This model posits that $y_{t}$, the fluoresence at the $t$th timestep, is a noisy realization of $c_{t}$, the unobserved underlying calcium concentration at the $t$th timestep. In the absence of a spike at the $t$th timestep ($s_{t}=0$), the calcium concentration decays according to a $p$th-order auto-regressive process. However, if a spike occurs at the $t$th timestep ($s_{t}>0$), then the calcium concentration increases. Thus,

[TABLE]

In (1), the quantities $\gamma_{1},\ldots,\gamma_{p}$ are the parameters in the auto-regressive model. Note that the quantity $y_{t}$ in (1) is observed; all other quantities are unobserved. Since we would like to know whether a spike occurred at the $t$th timestep, the parameter of interest is $s_{t}$. Figure 1(a) displays a small dataset generated according to (1).

In what follows, for ease of exposition, we assume $\beta_{0}=0$ and $\beta_{1}=1$ in (1). This assumption is made without loss of generality, since $\beta_{0}$ and $\beta_{1}$ can be estimated from the data, and the observed fluorescence $y_{1},\ldots,y_{T}$ centered and scaled accordingly. See Section 5 for additional details.

Vogelstein et al. (2010), Friedrich and Paninski (2016), and Friedrich, Zhou and Paninski (2017) seek to interpret $s_{t}$ in (1) as the number of spikes at the $t$th timestep. Thus, in principle it would be desirable to use a count-valued distribution, such as the Poisson distribution, as a prior on $s_{t}$. However, because maximum a posteriori estimation of $s_{t}$ in (1) using a Poisson distribution is computationally intractable, they instead suppose that $s_{t}$ has an exponential distribution (Vogelstein et al., 2010). In the case of the first-order auto-regressive process ($p=1$), this leads Vogelstein et al. (2010) to the optimization problem

[TABLE]

where $\lambda$ is a nonnegative tuning parameter that controls the trade-off between the fit of the estimated calcium to the observed fluorescence, and the sparsity of the estimated spike vector $\hat{s}_{2},\ldots,\hat{s}_{T}$. Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017) instead consider a closely-related problem that results from including an additional $\ell_{1}$ penalty for the initial calcium concentration,

[TABLE]

Both (2) and (3) are convex optimization problems, which can be solved for the global optimum using a well-developed set of optimization algorithms (Boyd and Vandenberghe, 2004; Hastie, Tibshirani and Friedman, 2009; Hastie, Tibshirani and Wainwright, 2015; Bien and Witten, 2016). Because $\hat{s}_{2},\ldots,\hat{s}_{T}$ are not integer-valued, they cannot be directly interpreted as the number of spikes at each timestep; however, informally, a larger value of $\hat{s}_{t}$ can be interpreted as indicating greater certainty that one or more spikes occurred at the $t$th timestep.

In this paper, we re-consider the model (1) that originally motivated the optimization problems (2) and (3) in the recent literature (Vogelstein et al., 2010; Friedrich and Paninski, 2016; Friedrich, Zhou and Paninski, 2017). Rather than interpreting $s_{t}$ in (1) as the number of spikes at the $t$th timestep, we interpret its sign as an indicator for whether or not at least one spike occurred: that is, $s_{t}=0$ indicates no spikes at the $t$th timestep, and $s_{t}>0$ indicates the occurrence of at least one spike. Then, in the case of a first-order auto-regressive model ($p=1$), (1) leads naturally to the optimization problem

[TABLE]

where $1_{(A)}$ is an indicator variable that equals $1$ if the event $A$ holds, and [math] otherwise. In (4), $\lambda$ is a non-negative tuning parameter that controls the trade-off between the fit of the estimated calcium to the observed fluorescence, and the number of timesteps at which a spike is estimated to occur.

Unfortunately, the optimization problem (4) is highly non-convex, due to the presence of the indicator variable. In the statistics literature, this term is known as an $\ell_{0}$ penalty. It is well-known that optimization involving $\ell_{0}$ penalties is typically computationally intractable: in general, no efficient algorithms are available to solve for the global optimum.

In fact, the convex optimization problem (2) considered in Vogelstein et al. (2010), and its close cousin (3) considered in Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017), can be viewed as convex relaxations to the problem (4). That is, if we replace the $\ell_{0}$ penalty in (4) with an $\ell_{1}$ penalty, then that we arrive exactly at the problem (2).

1.2 Contribution of This Paper

In the previous subsection, we established that the optimization problems (2) and (3), studied by Vogelstein et al. (2010), Friedrich and Paninski (2016), and Friedrich, Zhou and Paninski (2017), can be seen as convex relaxations of the $\ell_{0}$ optimization problem (4), which follows directly from the model (1). In fact, under the model (1), (4) is the “right” optimization problem to be solving; (2) and (3) are simply computationally tractable approximations to this problem. (In fact, Friedrich, Zhou and Paninski (2017) allude to this in the “Hard shrinkage and $\ell_{0}$ penalty” section of their paper.)

However, using an $\ell_{1}$ norm to approximate an $\ell_{0}$ norm comes with computational advantages at the expense of substantial performance disadvantages: in particular, the use of an $\ell_{1}$ penalty tends to overshrink the fitted estimates (Zou, 2006). This can be seen quite clearly in Figures 1(b) and 1(c). Retaining only the four spikes in Figure 1(c) associated with the largest increases in calcium leads to an improvement in spike detection (Figure 1(e); this is referred to as the post-thresholding $\ell_{1}$ estimator in what follows), but still one of the four true spikes is missed.

In this paper, we consider a slight modification of (4) that results from removing the positivity constraint,

[TABLE]

In practice, the distinction between the problems (5) and (4) is quite minor: on real data applications, for appropriate choices of the decay rate $\gamma$, the solution to (5) tends to satisfy the constraint in (4), and so the solutions are identical.

Like problem (4), solving problem (5) for the global optimum appears, at a glance, to be computationally intractable — we (the authors) are only aware of a few $\ell_{0}$ optimization problems for which exact solutions can be obtained via efficient algorithms!

However, in this paper, we show that in fact, (5) is a rare $\ell_{0}$ optimization problem that can be exactly solved for the global optimum using an efficient algorithm. This is because (5) can be seen as a changepoint detection problem, for which efficient algorithms that run in no more than $\mathcal{O}(T^{2})$ time, and often closer to $\mathcal{O}(T)$ time, are available. Furthermore, our implementation of the exact algorithm for solving (5) yields excellent results relative to the convex approximation (3) considered by Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017). This vastly improved performance can be seen in Figure 1(d).

The rest of this paper is organized as follows. In Section 2, we present an exact algorithm for solving the $\ell_{0}$ problem (5). In Section 3, we investigate the performance of this algorithm, relative to the algorithm of Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017) for solving the $\ell_{1}$ problem (3), in a simulation study. In Section 4, we investigate the performances of both algorithms for spike train inference on a data set for which the true spike times are known (Chen et al., 2013; GENIE Project, 2015) and on a data set from the Allen Brain Observatory (Allen Institute for Brain Science, 2016; Hawrylycz et al., 2016). In Section 5 we generalize the problem (5) in order to allow for efficient estimation of an intercept term, and to accommodate an auto-regressive model of order $p>1$ (1). Finally, we close with a discussion in Section 6. Technical details and additional results can be found in the Appendix.

2 An Exact Algorithm for Solving Problem (5)

In Section 2.1, we show that problem (5) can be viewed as a changepoint detection problem. In Sections 2.2 and 2.3, we apply existing algorithms for changepoint detection in order to efficiently solve (5) for the global optimum in $\mathcal{O}(T^{2})$ and in substantially fewer than $O(T^{2})$ operations, respectively. Timing results are presented in Section 2.4. We discuss selection of the tuning parameter $\lambda$ and auto-regressive parameter $\gamma$ in (5) in Appendix B.

2.1 Recasting (5) as a Changepoint Detection Problem

Recall that our goal is to solve the $\ell_{0}$ optimization problem (5), or equivalently, to compute $\hat{c}_{1},\ldots,\hat{c}_{T}$ that solve the optimization problem

[TABLE]

We estimate a spike event at the $t$th timestep if $\hat{c}_{t}\neq\gamma\hat{c}_{t-1}$. (We refer to this as a “spike event”, rather than a spike, since $\hat{c}_{t}\neq\gamma\hat{c}_{t-1}$ indicates the presence of at least one spike at the $t$th timepoint, but does not directly provide an estimate of the number of spikes.) We now make two observations about this optimization problem.

Given that a spike event is estimated at the $t$th timestep, the estimated calcium concentration at any time $t_{1}<t$ is independent of the estimated calcium concentration at any time $t_{2}\geq t$. 2. 2.

Given that two spike events are estimated at the $t$th and $t^{\prime}$th timesteps with $t<t^{\prime}$, and no spike events are estimated in between the $t$th and $t^{\prime}$th timesteps, the calcium concentration is estimated to decay exponentially between the $t$th and $t^{\prime}$th timesteps.

This motivates us to consider the relationship between (5) and a changepoint detection problem (Aue and Horváth, 2013; Braun and Muller, 1998; Davis, Lee and Rodriguez-Yam, 2006; Yao, 1988; Lee, 1995; Jackson et al., 2005; Killick, Fearnhead and Eckley, 2012; Maidstone et al., 2016) of the form

[TABLE]

where

[TABLE]

In (6), we are simultaneously minimizing the objective over the times at which the changepoints ($\tau_{1},\ldots,\tau_{k}$) occur and the number of changepoints ($k$); the parameter $\lambda$ controls the relative importance of these two terms.

The following result establishes an equivalence between (6) and (5).

Proposition 1

*There is a one-to-one correspondence between the set of estimated spike events in the solution to (5) and the set of changepoints $0=\tau_{0},\tau_{1},\ldots,\tau_{k},\tau_{k+1}=T$ in the solution to (6), in the sense that $\hat{c}_{t}\neq\gamma\hat{c}_{t-1}$ if and only if $t\in\{\tau_{1}+1,\ldots,\tau_{k}+1\}$. Furthermore, given the set of changepoints, the solution to (5) takes the form *

[TABLE]

for $j\in\left\{0,\ldots,k\right\}$.

Proposition 1 indicates that in order to solve (5), it suffices to solve (6). (We note that due to a slight discrepancy between the conventions used in the changepoint detection literature and the notion of a spike event in this paper, the indexing in Proposition 1 is a little bit awkward, in the sense that the $k$th spike event is estimated to occur at time $\tau_{k}+1$, rather than at time $\tau_{k}$.)

In the next two sections, we will make use of the following result.

Proposition 2

The quantity (7) has a closed-form expression,

[TABLE]

*Furthermore, given ${\mathcal{D}}(y_{a:b})$, we can calculate ${\mathcal{D}}(y_{a:(b+1)})$ in constant time.

Propositions 1 and 2 are proven in Appendix A.

2.2 An Algorithm for Solving (5) in $\mathcal{O}(T^{2})$ Operations

In this section, we apply a dynamic programming algorithm proposed by Jackson et al. (2005) and Auger and Lawrence (1989) in order to solve the changepoint detection problem (6) for the global optimum in $\mathcal{O}(T^{2})$ time. Due to the equivalence between (6) and (5) established in Proposition 1, this algorithm also solves problem (5).

Roughly speaking, this algorithm recasts the very difficult problem of choosing the times of all changepoints simultaneously into the much simpler problem of choosing the time of just the most recent changepoint. In greater detail, consider solving (6) on the first $s$ timesteps. Define $F(0)\equiv-\lambda$, and for $s\geq 1$, define

[TABLE]

In other words, in order to solve (6), we need simply identify the time of the most recent changepoint, and then solve (6) on all earlier times.

This recursion gives a simple recipe for evaluating $F(T)$ efficiently: set $F(0)=-\lambda$, and compute $F(1),F(2),\ldots,F(T)$ based on previously calculated (and stored) values. For example, at $s=1$, calculate and store

[TABLE]

and then at $s=2$ use the previously calculated values $F(0)$ and $F(1)$ to compute the minimum over a finite set with $s$ elements

[TABLE]

Given $F(1),\ldots,F(s-1)$, computing $F(s)$ requires minimizing over a finite set of size $s$, and therefore it has computational cost linear in $s$. The total cost of computing $F(T)$ is quadratic in the total number of timesteps, $T$, since there are $T+1$ subproblems: $\sum_{s=0}^{T}s=\mathcal{O}(T^{2})$.

Full details are provided in Algorithm 1. We note that this algorithm is particularly efficient in light of Proposition 2, which makes it possible to perform a constant-time update to ${\mathcal{D}}(y_{(\tau+1):s})$ in order to compute ${\mathcal{D}}(y_{(\tau+1):(s+1)})$.

2.3 Dramatic Speed-Ups Using Cost-Complexity Pruning

In a recent paper, Killick, Fearnhead and Eckley (2012) considered problems of the form (6) for which an assumption on ${\mathcal{D}}(\cdot)$ holds; this assumption is satisfied by (7).

The main insight of their paper is as follows. Suppose that $s<r$ and $F(s)+{\mathcal{D}}(y_{(s+1):r})>F(r)$. Then for any $q>r$, it is mathematically impossible for the most recent changepoint before the $q$th timestep to have occurred at the $s$th timestep. This allows us to prune the set of candidate changepoints that must be considered in each step of Algorithm 2, leading to drastic speed-ups. Details are provided in Algorithm 2, which solves (5) for the global optimum.

Under several technical assumptions, Killick, Fearnhead and Eckley (2012) show that the expected complexity of this algorithm is $\mathcal{O}(T)$. The main assumption is that the expected number of changepoints in the data increases linearly with the length of the data; this is reasonable in the context of calcium imaging data, in which we expect the number of neuron spike events to be linear in the length of the recording.

2.4 Timing Results for Solving (5)

We simulated data from (1) with $\gamma=0.998$, $\sigma=0.15$, and $s_{t}\sim_{\text{ind.}}\text{Poisson}(\theta)$ with $\theta\in\{0.1,0.01,0.001\}$. We solved (5) with $\lambda=1$, using our R-language implementations of Algorithms 1 and 2.

Timing results, averaged over 50 simulated data sets, are displayed in Figure 2. As expected, the running time of Algorithm 1 scales quadratically in the length of the time series, whereas the running time of Algorithm 2 is upper bounded by that of Algorithm 1. Furthermore, the running time of Algorithm 2 decreases as the firing rate increases. The Chen et al. (2013) dataset explored in Section 4.1 has firing rate on the same order of magnitude as the middle panel, $\theta=0.01$. Using Algorithm 2, we can solve (5) for the global optimum in a few minutes on a 2.5 GHz Intel Core i7 Macbook Pro on fluorescence traces of length $100,000$ with moderate to high firing rates.

We note here that Algorithm 2 for solving (5) is much slower than the algorithm of Friedrich, Zhou and Paninski (2017) for solving (3), which is implemented in Cython and has approximately linear running time. It should be possible to develop a faster algorithm for solving (5) using ideas from Johnson (2013), Maidstone et al. (2016), and Hocking et al. (2017). Furthermore, a much faster implementation of Algorithm 2 would be possible using a language other than R. We leave such improvements to future work.

3 Simulation Study

3.1 Comparison Methods

In this section, we use in silico data to demonstrate the performance advantages of the $\ell_{0}$ approach (5) over two competing approaches:

The $\ell_{1}$ proposal (3) of Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017), which involves a single tuning parameter $\lambda$. 2. 2.

A thresholded version of the $\ell_{1}$ estimator. Letting $\hat{s}_{2},\ldots,\hat{s}_{T}$ denote the solution to (3), we define the post-thresholding estimator as

[TABLE]

for $L$ a positive constant. In other words, the post-thresholding estimator retains only the estimated spikes for which the estimated increase in calcium exceeds a threshold $L$. The post-thresholding estimator involves two tuning parameters: $\lambda$ in (3), as well as the value of $L$ used to perform thresholding.

The post-thresholding estimator is motivated by the fact that the solution to (3) tends to yield many “small” spikes: i.e. $\hat{s}_{t}$ is near zero, but not exactly equal to zero, for many timesteps. In fact, this can be seen in Figure 1(c). As seen in Figure 1(e), the post-thresholding estimator has the potential to improve the performance of the $\ell_{1}$ estimator by removing some of these small spikes. Of course, the post-thresholding estimator with $L=0$ is identical to the $\ell_{1}$ estimator from (3).

3.2 Performance Measures

We measure performance of each method based on two criteria: (i) error in calcium estimation, and (ii) error in spike detection.

We consider the mean of squared differences between the true calcium concentration in (1) and the estimated calcium concentration that solves (5),

[TABLE]

This quantity involves the unobserved calcium concentrations, $c_{1},\ldots,c_{T}$, and thus can only be computed on simulated data. Furthermore, this quantity can be computed for our $\ell_{0}$ proposal (5) and for the $\ell_{1}$ proposal (3), but not for the post-thresholding estimator (9), since the post-thresholding estimator does not lead to an estimate of the underlying calcium concentrations.

We now consider the task of quantifying the error in spike detection. We make use of the Victor-Purpura distance metric (Victor and Purpura, 1996, 1997), which defines the distance between two spike trains as the minimum cost of transforming one spike train to the other through spike insertion, deletion, or translation. We also use the van Rossum distance (van Rossum, 2001), defined as the mean squared difference between two spike trains that have been convolved with an exponential kernel with timescale $\tau=2$.

3.3 Results

We generated 100 simulated data sets according to (1) with parameter settings $\gamma=0.96$, $T=5000$, $\sigma=0.15$, and $s_{t}\sim_{\text{i.i.d.}}\text{Poisson(0.01)}$.

On each simulated data set, we solved (5) and (3) for a range of values of the tuning parameter $\lambda$. Moreover, we post-thresholded the $\ell_{1}$ solution, as in (9), with five different threshold values: $L\in\{0,0.125,0.250,0.375,0.500\}$.

Figure 3(a) displays the error in spike event detection for the van Rossum distance, Figure 3(b) displays the error in spike event detection for the Victor-Purpura distance metric, and Figure 3(c) displays the error in calcium estimation (10), for the $\ell_{0}$ problem (5) and the $\ell_{1}$ problem (3), for a range of values of $\lambda$. Results are averaged over the 50 simulated data sets.

As mentioned earlier, since the calcium concentration is not defined for the post-thresholding estimator (9), the post-thresholding estimator is not displayed in Figure 3(c). In Figures 3(a) and 3(b), five distinct curves are displayed for the post-thresholding operator; each corresponds to a distinct value of $L$. Note that as $L$ increases, the maximum possible number of estimated spikes from the post-thresholding estimator decreases. For example, with $\lambda=0$ and $L=0.5$, no more than approximately 50 spikes are estimated by the post-thresholding estimator. For this reason, some of the curves corresponding to the post-thresholding estimator appear truncated in Figures 3(a) and 3(b).

Figure 3 reveals that the $\ell_{0}$ estimator (5) results in dramatically lower errors in both calcium estimation and spike detection than the $\ell_{1}$ estimator (3) (which is equivalent to the post-thresholding operator with $L=0$). Although post-thresholding with $L>0$ improves upon the unthresholded $\ell_{1}$ estimator, the $\ell_{0}$ estimator still substantially outperforms all competitors in Figures 3(a) and 3(b). Moreover, the $\ell_{0}$ estimator requires just a single tuning parameter $\lambda$ in (5), whereas the post-thresholding procedure involves two tuning parameters, $\lambda$ in (3) and $L$ in (9), leading to challenges in tuning parameter selection.

Furthermore, the $\ell_{0}$ problem (5) achieves the lowest errors in both calcium estimation and spike detection when applied using a value of the tuning parameter $\lambda$ that yields approximately 50 estimated spikes, which is the expected number of spikes in this simulation. This suggests that it should be possible to use a cross-validation scheme to select the tuning parameter $\lambda$ for the $\ell_{0}$ approach; we propose such a scheme in Appendix B. By contrast, in Figure 3(b), the $\ell_{1}$ approach achieves its lowest error in calcium estimation when far more than 50 spikes are estimated. This is a consequence of the fact that the $\ell_{1}$ penalty simultaneously reduces the number of estimated spikes and shrinks the estimated calcium. Therefore, the value of the tuning parameter $\lambda$ in (3) that yields the most accurate estimate of calcium will result in severe over-estimation of the number of spikes. This means that the cross-validation scheme detailed in Appendix B will not perform well for the $\ell_{1}$ approach.

4 Application to Calcium Imaging Data

In this section, we apply our $\ell_{0}$ proposal (5) and the $\ell_{1}$ proposal of Friedrich and Paninski (2016) and Friedrich, Zhou and Paninski (2017) (3), both with and without post-thresholding (9), to two calcium imaging data sets. In the first data set, the true spike times are known (Chen et al., 2013; GENIE Project, 2015), and so we can directly assess the spike detection accuracy of each proposal. In the second data set, the true spike times are unknown (Allen Institute for Brain Science, 2016; Hawrylycz et al., 2016); nonetheless, we are able to make a qualitative comparison of the results of the $\ell_{1}$ and $\ell_{0}$ proposals.

4.1 Application to Chen et al. (2013) Data

We first consider a data set that consists of simultaneous calcium imaging and electrophysiological measurements (Chen et al., 2013; GENIE Project, 2015), obtained from the Collaborative Research in Computational Neuroscience portal (http://crcns.org/data-sets/methods/cai-1/about-cai-1). In what follows, we refer to the spike times inferred from the electrophysiological measurements as the “true” spikes.

The top panel of Figure 4 shows a 40-second recording from cell 2002, which expresses GCaMP6s. The data are measured at 60 Hz, for a total of 2400 timesteps. The raw fluorescence traces are DF/F transformed with a $20\%$ percentile filter as in Figure 3 of Friedrich, Zhou and Paninski (2017). In this 40-second recording, there are a total of 23 true spikes; therefore, we solved the $\ell_{0}$ and $\ell_{1}$ problems with $\gamma=0.9864405$ using values of $\lambda$ in (5) and (3) that yield 23 estimated spikes. Additionally, we solved the $\ell_{1}$ problem with $\lambda=1$, and post-thresholded it according to (9) using $L=0$, $0.1$, and $0.13$; these threshold values yielded $230$, $54$, and $23$ estimated spikes, respectively.

Figure 4 displays the estimated spikes resulting from the $\ell_{0}$ proposal, the estimated spikes resulting from the $\ell_{1}$ proposal, the estimated spikes from post-thresholding the $\ell_{1}$ solutions, and the ground truth spikes. We see that the $\ell_{0}$ proposal has one false negative (i.e. it misses one true spike at around $7$ seconds) and one false positive (i.e. it estimates a spike at around 36 seconds, where there is no true spike). By contrast, the $\ell_{1}$ problem concentrates the 23 estimated spikes at three points in time, and therefore suffers from a substantial number of false positives as well as false negatives. Because the $\ell_{1}$ penalty controls both the number of spikes and the estimated calcium, the $\ell_{1}$ problem tends to put a large number of spikes in a row, each of which is associated with a very modest increase in calcium. This is consistent with the results seen in Figures 1 and 3. Post-thresholding the $\ell_{1}$ estimator does lead to an improvement in results relative to the unthresholded $\ell_{1}$ method; however, the post-thresholded solution with 23 spikes still tends to estimate a number of spikes in short succession when in fact only one true spike is present, and also misses several true spike events.

We note that the $\ell_{0}$ method tends to estimate spike times one or two timesteps ahead of the true spike times. This is due to model misspecification: model (1) with $p=1$ assumes that the calcium concentration increases instantaneously due to a spike event, and subsequently decays; however, we see from Figure 4 that in reality, a spike event is followed by an increase in calcium over the course of a few timesteps, before the onset of exponential decay. We see two possible avenues to address this relatively minor issue: estimated spike times from the $\ell_{0}$ method can be adjusted to account for this empirical observation; or else the optimization problem (5) can be adjusted in order to allow for more realistic calcium dynamics (for example, by solving an $\ell_{0}$ optimization problem corresponding to (1) with $p>1$). We explore the second alternative in Section 5.

In Appendix C, we apply an approach proposed by Friedrich, Zhou and Paninski (2017) to approximate the solution to a non-convex problem using a greedy algorithm. This alternative approach performs quite a bit better than solving the $\ell_{1}$ problem (3); however, it does not achieve the global optimum.

4.2 Application to Allen Brain Observatory Data

We now consider a data set from the Allen Brain Observatory, a large open-source repository of calcium imaging data from the mouse visual cortex (Allen Institute for Brain Science, 2016; Hawrylycz et al., 2016). For this data, the true spike times are not available, and so it is difficult to objectively assess the performances of the $\ell_{1}$, post-thresholded $\ell_{1}$, and $\ell_{0}$ methods. Instead, for each method we present several fits that differ in the number of detected spikes. We argue that the $\ell_{0}$ problem yields results that are qualitatively superior to those of the competitors, in the sense that they are better supported by visual inspection of the data.

For the second ROI in NWB 510221121, we applied the $\ell_{1}$, post-thresholded $\ell_{1}$, and $\ell_{0}$ methods to the first $10,000$ timesteps of the $DF/F$-transformed fluorescence traces. Since the data are measured at 30 Hz, this amounts to the first $333$ seconds of the recording. Figure 5 shows the results obtained with $\gamma=0.981756$. For the $\ell_{0}$ and $\ell_{1}$ estimators, we chose the values of $\lambda$ in (3) and (5) in order to obtain 27, 49, and 128 estimated spikes. For the post-thresholded estimator (9), we set $\lambda=1$, and then selected $L$ to yield 27, 49, and 128 estimated spikes.

As in the previous subsection, we see that when faced with a large increase in fluorescence, the $\ell_{1}$ problem tends to estimate a very large number of spikes in quick succession. For example, when 27 spikes are estimated, the $\ell_{1}$ problem concentrates the estimated spikes at three points in time (Figure 5(a)). Even when 128 spikes are estimated, the $\ell_{1}$ problem still seems to miss all but the largest peaks in the fluorescence data (Figure 5(c)). Post-thresholding the $\ell_{1}$ estimator improves upon this issue somewhat, but spikes corresponding to smaller increases in fluorescence are still missed; this issue can be clearly seen in Figures 5(d)–(f), which zoom in on a smaller time window.

By contrast, the $\ell_{0}$ problem can assign an arbitrarily large increase in calcium to a single spike event. Therefore, it seems to capture most of the visible peaks in the fluorescence data when 49 spikes are estimated (Figures 5(b) and 5(e)), and it captures all of them when 128 spikes are estimated (Figures 5(c) and 5(f)).

Though the true spike times are unknown for this data, based on visual inspection, the results for the $\ell_{0}$ proposal seem superior to those of the $\ell_{1}$ and post-thresholded $\ell_{1}$ proposals.

5 Extensions

We now present two straightforward extensions to the optimization problem (5), for which computationally attractive algorithms along the lines of the one proposed in Section 2 are available.

5.1 Estimation of the Intercept in (1)

The model for calcium dynamics considered in this paper (1) allows for an intercept term, $\beta_{0}$. In order to arrive at (5), we assumed that the intercept was known and (without loss of generality) equal to zero. However, in practice, we might want to fit the model (1) without knowing the value of the intercept $\beta_{0}$. In fact, in many settings, this may be of great practical importance, since the meaning of the model (1) (and, for instance, the rate of exponential decay $\gamma$) is inextricably tied to the value of the intercept.

We now propose a modification to the $\ell_{0}$ optimization problem (5) that allows for estimation of the intercept $\beta_{0}$. So that the resulting problem can be efficiently solved using the ideas laid out in Section 2, we must ensure that given the estimated spike times, the calcium can be estimated separately between each pair of consecutive spikes. Consequently, we must allow for a separate intercept term between each pair of consecutive spikes. This suggests the optimization problem

[TABLE]

where the indicator variable $1_{(A,B)}$ equals 1 if the event $A\cup B$ holds, and equals zero otherwise. Therefore, $1_{\left(c_{t}\neq\gamma c_{t-1},\beta_{0t}\neq\beta_{0,t-1}\right)}$ equals $1$ if the calcium concentration stops decaying or if the intercept changes. Note that in the solution to (11), the intercept is constant between adjacent timesteps, unless there is a spike.

Problem (11) can be recast as a changepoint problem of the form (6) with

[TABLE]

Given ${\mathcal{D}}(y_{a:b})$, ${\mathcal{D}}(y_{a:(b+1)})$ can be updated in constant time. Thus, the algorithms introduced in Section 2 can be easily modified in order to solve (11) for the global optimum.

5.2 An Auto-Regressive Model with $p>1$ in (1)

The model (1) allows for the calcium dynamics to follow a $p$th order auto-regressive process. For simplicity, this paper focused on the case of $p=1$. We now consider developing an $\ell_{0}$ optimization problem for the model (1) with $p>1$.

It is natural to consider the $\ell_{0}$ optimization problem

[TABLE]

However, (12) cannot be expressed in the form (6): the penalty in (12) induces a dependency in the calcium that spans more than two timesteps, so that the calcium at a given timestep may depend on the calcium prior to the most recent spike. As a result, (12) is computationally intractable.

Instead, we consider a changepoint detection problem of the form (6) with cost function

[TABLE]

Thus, the calcium follows a $p$th order auto-regressive model between any pair of spikes; furthermore, once a spike occurs, the calcium concentrations are reset completely. That is, the calcium after a spike is not a function of the calcium before a spike. Consequently, it is straightforward to develop a fast algorithm for solving this changepoint detection problem for the global optimum, using ideas detailed in Section 2.

In particular, a popular model for calcium dynamics assumes that between any pair of spikes, the calcium can be well-approximated by the difference between two exponentially-decaying functions (Brunel and Wang, 2003; Mazzoni et al., 2008; Cavallari, Panzeri and Mazzoni, 2016; Volgushev, Ilin and Stevenson, 2015). This would perhaps be a better model for the data from the Allen Brain Observatory, in which increases in fluorescence due to a spike occur over the course of a few timesteps, rather than instantaneously (see Figure 5). This “difference of exponentials” models falls directly within the framework of (1) with $p=2$, and hence could be handled handled using the changepoint detection problem just described.

6 Discussion

In this paper, we considered solving the seemingly intractable $\ell_{0}$ optimization problem (5) corresponding to the model (1). By recasting (5) as a changepoint detection problem, we were able to derive an algorithm to solve (5) for the global optimum in expected linear time. It should be possible to develop an even more efficient algorithm for solving (5) that exploits recent algorithmic developments for changepoint detection (Johnson, 2013; Maidstone et al., 2016; Hocking et al., 2017); we leave this as an avenue for future work.

We have shown in this paper that solving the $\ell_{0}$ optimization problem (5) leads to more accurate spike event detection than solving the $\ell_{1}$ optimization problem (3) proposed by Friedrich, Zhou and Paninski (2017). Indeed, this finding is intuitive: the $\ell_{1}$ penalty and positivity constraint in (3) serves as a exponential prior on the increase in calcium at any given time point, and thereby effectively limits the amount that calcium can increase in response to a spike event. By contrast, the $\ell_{0}$ penalty in (5) is completely agnostic to the amount by which a spike event increases the level of calcium. Consequently, it can allow for an arbitrarily large (or small) increase in fluorescence as a result of a spike event.

While approximations to the solution to the $\ell_{0}$ problem (5) are possible (de Rooi and Eilers, 2011; de Rooi, Ruckebusch and Eilers, 2014; Hugelier et al., 2016; Scott and Knott, 1974; Olshen et al., 2004; Fryzlewicz et al., 2014; Friedrich, Zhou and Paninski, 2017), there is no guarantee that such approaches will yield an attractive local optimum on a given data set. In this paper, we completely bypass this concern by solving the $\ell_{0}$ problem for the global optimum.

In this paper, we have focused on the empirical benefits of the $\ell_{0}$ problem (5) over the $\ell_{1}$ problem (3). However, it is natural to wonder whether these empirical benefits are backed by statistical theory. Conveniently, both the $\ell_{0}$ and $\ell_{1}$ optimization problems are very closely-related to problems that have been well-studied in the statistical literature from a theoretical standpoint. In particular, in the special case of $\gamma=1$, the $\ell_{0}$ problem (5) was extensively studied in Yao and Au (1989) and Boysen et al. (2009). Furthermore, when $\gamma=1$, the $\ell_{1}$ problem (5) is very closely-related to the fused lasso optimization problem,

[TABLE]

which has also been extremely well-studied (Tibshirani et al., 2005; Mammen et al., 1997; Davies and Kovac, 2001; Rinaldo et al., 2009; Harchaoui and Lévy-Leduc, 2010; Qian and Jia, 2012; Rojas and Wahlberg, 2014; Lin et al., 2016; Dalalyan et al., 2017). However, we leave a formal theoretical analysis of the relative merits of (5) and (3), in terms of $\ell_{2}$ error bounds and spike recovery properties, to future work.

Our R-language software for our proposal is available on CRAN in the package

LZeroSpikeInference. Instructions for running this software in python can be found at https://github.com/jewellsean/LZeroSpikeInference.

Acknowledgments

We are grateful to three reviewers and an associate editor for helpful comments that improved the quality of this work. We thank Michael Buice and Kyle Lepage at the Allen Institute for Brain Science for conversations that motivated this work. Johannes Friedrich at Columbia University provided assistance in using his software to solve the $\ell_{1}$ problem (3). Sean Jewell received funding from the Natural Sciences and Engineering Research Council of Canada. This work was partially supported by NIH Grant DP5OD009145 and NSF CAREER Award DMS-1252624 to Daniela Witten.

A Proof of Propositions

A.1 Proof of Proposition 1

Proof A.1**.**

The first sentence follows by inspection. To establish the second sentence, we observe that the cost

[TABLE]

can be rewritten by direct substitution of the constraint as

[TABLE]

This is a least squares problem and is minimized at

[TABLE]

which implies that

[TABLE]

and furthermore that for $a<t\leq b$ the fitted values are $\hat{c}_{t}=\gamma\hat{c}_{t-1}$. Applying this argument to each segment gives the result stated in Proposition 1.

A.2 Proof of Proposition 2

Proof A.2**.**

The first equation follows by expanding the square for the final form of ${\mathcal{D}}(y_{a:b})$ in the proof of Proposition 1. Given ${\mathcal{D}}(y_{a:b})$ we can calculate ${\mathcal{D}}(y_{a:(b+1)})$ in constant time by storing $\sum_{t=a}^{b}\frac{y_{t}^{2}}{2}$ and $\sum_{t=a}^{b}y_{t}\gamma^{t-a}$, and updating each of these sums for the new data point $y_{b+1}$; we use a closed form expression to calculate $\sum_{t=a}^{b+1}\gamma^{2(t-a)}$. With each of these quantities stored, ${\mathcal{D}}(\cdot)$ and ${\mathcal{C}}(\cdot)$ are updated in constant time.

B Choosing $\lambda$ and $\gamma$

Recall that in (5), the parameters $\lambda$ and $\gamma$ are unknown. The nonnegative parameter $\lambda$ controls the trade-off between the number of estimated spike events and the quality of the estimated calcium fit to the observed fluorescence. The parameter $\gamma,0<\gamma<1$, controls the rate of exponential decay of the calcium. We consider two approaches for choosing $\gamma$ and $\lambda$.

B.1 Approach 1

To estimate $\gamma$, we manually select a segment $y_{a:b}$ that, based on visual inspection, appears to exhibit exponential decay. We then estimate $\gamma$ as

[TABLE]

This can be done via numerical optimization.

Next, given $\gamma$, we select $\lambda$ via cross-validation. For each value of $\lambda$ that we consider, we solve (6) on a training set, and then evaluate the mean squared error (MSE) on a hold-out set. Details are provided in Algorithm 3.

B.2 Approach 2

Pnevmatikakis et al. (2013), Friedrich and Paninski (2016), and Friedrich, Zhou and Paninski (2017) propose to select the exponential decay parameter $\gamma$ based on the autocovariance function, and to choose the tuning parameter $\lambda$ such that $||y-\hat{c}||_{2}\leq\sigma\sqrt{T}$, where the standard deviation $\sigma$ is estimated through the power spectral density of $y$, and $T$ is the number of timepoints. We refer the reader to Friedrich, Zhou and Paninski (2017) and Pnevmatikakis et al. (2016) for additional details.

C A Greedy Approach for Approximating the Solution to a Non-Convex Problem

Friedrich, Zhou and Paninski (2017) consider a variant of the optimization problem (3),

[TABLE]

obtained from (3) by setting $\lambda=0$, and changing the convex positivity constraint to the non-convex constraint that $s_{t}$ lies within a non-convex set. Like (5), (13) is non-convex. Friedrich, Zhou and Paninski (2017) do not attempt to solve (13) for the global optimum; instead, they provide a heuristic modification to their algorithm for solving (5), which is intended to approximate the solution to (13).

Figure 6 illustrates the behavior of this approximate algorithm when applied to the same data as in Figure 4. We set $\gamma=0.9864405$, and considered three values of $s_{min}$. When $s_{min}=10^{-8}$ and $s_{min}=0.1$, in panels (a) to (b), too many spikes are estimated. But when $s_{min}=0.3$, in panel (c), the solution to (13) is very similar to the solution to (5) with $\lambda=0.6$. Both almost perfectly recover the ground truth spikes. Therefore, in this example, the approximate algorithm of Friedrich, Zhou and Paninski (2017) for solving (13) performs quite well.

However, (13) is a non-convex problem, and the approximate algorithm of Friedrich, Zhou and Paninski (2017) is not guaranteed to find the global minimum. In fact, we can see that on the data shown in Figure 6, this approximate algorithm does not find the global optimum. When applied with $s_{min}=0.3$, the approximate algorithm yields an objective value of $8.57$. By contrast, our algorithm for solving (5) yields a solution that is feasible for (13), and which results in a value of $7.86$ for the objective of (13). We emphasize that this is quite remarkable: even though the algorithm proposed in Section 2 solves (5) and not (13), it nonetheless yields a solution that is closer to the global optimum of (13) than does the approximate algorithm of Friedrich, Zhou and Paninski (2017), which is intended to solve (13).

In many cases, the greedy algorithm of Friedrich, Zhou and Paninski (2017) for solving (13) might yield good results that are near the global optimum of (13), and potentially even near the global optimum of (5). However, there is no guarantee that this algorithm will yield a “good” local optimum on any given data set. By contrast, in this paper we have proposed an elegant and efficient algorithm for exactly solving the $\ell_{0}$ problem (5).