Fast Convolutional Dictionary Learning off the Grid

TL;DR

This paper introduces a novel convolutional dictionary learning framework that accurately estimates continuous-time events from sampled signals, using interpolated dictionaries and efficient algorithms, significantly improving off-the-grid event detection and template learning.

Contribution

The paper proposes a new CDL framework with interpolated dictionaries and a fast CSC algorithm, COMP-INTERP, enabling accurate off-the-grid event detection and improved template learning.

Findings

COMP-INTERP achieves similar accuracy to state-of-the-art CBP

COMP-INTERP is two orders of magnitude faster

Dictionary update with overcomplete dictionaries yields more accurate templates

Abstract

Given a continuous-time signal that can be modeled as the superposition of localized, time-shifted events from multiple sources, the goal of Convolutional Dictionary Learning (CDL) is to identify the location of the events--by Convolutional Sparse Coding (CSC)--and learn the template for each source--by Convolutional Dictionary Update (CDU). In practice, because we observe samples of the continuous-time signal on a uniformly-sampled grid in discrete time, classical CSC methods can only produce estimates of the times when the events occur on this grid, which degrades the performance of the CDU. We introduce a CDL framework that significantly reduces the errors arising from performing the estimation in discrete time. Specifically, we construct an expanded dictionary that comprises, not only discrete-time shifts of the templates, but also interpolated variants, obtained by bandlimited interpolation, that account for continuous-time shifts. For CSC, we develop a novel computationally efficient CSC algorithm, termed Convolutional Orthogonal Matching Pursuit with interpolated dictionary (COMP-INTERP). We benchmarked COMP-INTERP to Contiunuous Basis Pursuit (CBP), the state-of-the-art CSC algorithm for estimating off-the-grid events, and demonstrate, on simulated data, that 1) COMP-INTERP achieves a similar level of accuracy, and 2) is two orders of magnitude faster. For CDU, we derive a novel procedure to update the templates given sparse codes that can occur both on and off the discrete-time grid. We also show that 3) dictionary update with the overcomplete dictionary yields more accurate templates. Finally, we apply the algorithms to the spike sorting problem on electrophysiology recording and show their competitive performance.

Tables

Table 1. TABLE I : Notations

Symbol	Description
$𝐇$	Matrix
$𝐡$	Vector
$𝒮$	Set
${𝐇}_i$	$i^{\text{th}}$ column from $𝐇$
${𝐇}^c$	$c^{\text{th}}$ block column from $𝐇$
$𝐡\left[j\right]$	$j^{\text{th}}$ entry from $𝐡$
${𝒮}_i$	$i^{\text{th}}$ element from set $𝒮$
${𝒮}^j$	$j^{\text{th}}$ set
${𝐈}_{L\times L}$	Identity matrix of size $L\times L$
${𝐫}^{\left(t\right)}$	$𝐫$ at $t^{\text{th}}$ iteration
${\mathrm{𝟎}}_L$	a length-$L$ vector with all entries equal to $0$
$n_{j,i}^c$	$i^{\text{th}}$ event from source $c$ in $j^{\text{th}}$ window
$\ast$ , $\star$	convolution / cross-correlation
${\parallel \cdot \parallel }_p$	${\ell }_p$ norm

Table 2. TABLE II : Convolutional greedy methods

	Efficient selection	Efficient projection	${\Delta }_K$
COMP	Yes	Yes	$\Delta$
CMP	Yes	$\cdot$	$\Delta$
COMP-slow	Yes	No	$\Delta$
COMP-INTERP	Yes	Yes	$\Delta /K$

Table 3. TABLE III : Computation time (sec) as a function of T 𝑇 T

	0.5	1	2	3	4	5
COMP	0.041	0.066	0.114	0.175	0.223	0.271
CMP	0.027	0.049	0.089	0.131	0.180	0.237
COMP-slow	0.079	0.176	0.350	0.545	0.759	1.002
COMP-INTERP	0.243	0.445	0.781	1.161	1.565	2.038
CBP	4.24	13.55	51.75	145.27	$\cdot$	$\cdot$

Table 4. TABLE IV : Computation time (sec) as a function of N 1 + N 2 subscript 𝑁 1 subscript 𝑁 2 N_{1}+N_{2}

	10	20	30	40	50
COMP	0.096	0.204	0.267	0.340	0.438
CMP	0.082	0.168	0.218	0.276	0.345
COMP-slow	0.194	0.666	1.189	2.111	3.043
COMP-INTERP	0.744	1.416	1.811	2.377	3.112
CBP	152.44	134.25	135.10	132.60	133.53

Equations

$$y(t)=\sum _{c=1}^C\sum _{i=1}^{N_c}{x}_i^c{h}_c(t-{\tau }_i^c)+\epsilon (t),$$

$$\mathbf{y}[n]=\sum_{c=1}^{C}\big{(}\mathbf{x}_{i}^{c}\ast\mathbf{h}_{c}\big{)}[n]+\boldsymbol{\varepsilon}[n],$$

$$\mathbf{y}=\mathbf{Hx}+\mathbf{\epsilon },$$

$$\begin{split}&\min_{\{\mathbf{h}_{c}\}_{c=1}^{C},\{\mathbf{X}_{j}\}_{j=1}^{J}}\sum_{j=1}^{J}\Big{\lVert}\mathbf{Y}_{j}-\mathbf{H}\mathbf{X}_{j}\Big{\rVert}_{2}^{2}\text{ s.t. }\lVert\mathbf{X}_{j}\rVert_{0}\leq\beta_{1},\\ \end{split}$$

$$\begin{split}&\min_{\{\vec{\mathbf{x}}^{c}_{i}\}}\Big{\lVert}\mathbf{y}-\sum_{c=1}^{C}\sum_{p=1}^{P}\mathbf{\Phi}^{c,p}\mathbf{x}^{c}_{p}\Big{\rVert}_{2}^{2}+\lambda\sum_{c=1}^{C}\lVert\mathbf{x}^{c}_{1}\rVert_{1}\\ &s.t.\,\,\big{[}\mathbf{x}^{c}_{1}[i],\cdots,\mathbf{x}^{c}_{P}[i]\big{]}\in\mathcal{H}\text{ for }i=1,\cdots,N-L+1,\\ \end{split}$$

$$\mathbf{y}[n]=\sum_{c=1}^{C}\sum_{k=0}^{K-1}\big{(}\mathbf{x^{c^{k}}}\ast\mathbf{h}_{c^{k}}\big{)}[n]+\boldsymbol{\varepsilon}[n],$$

$$\mathbf{\widetilde{H}}=\begin{bmatrix}\mathbf{\widetilde{H}}^{1^{0}}\Big{|}\cdots\Big{|}\mathbf{\widetilde{H}}^{1^{K-1}}\Big{|}\cdots\Big{|}\mathbf{\widetilde{H}}^{C^{0}}\Big{|}\cdots\Big{|}\mathbf{\widetilde{H}}^{C^{K-1}}\end{bmatrix},$$

$$\begin{split}\tilde{h}_{c}(t)&=\mathbf{h}_{c}*\Big{(}\frac{\sin\pi t}{\pi t}\Big{)}\\ &=\sum_{m=-\infty}^{m=\infty}\mathbf{h}_{c}[m]\frac{\sin\pi(t-m\Delta)}{\pi(t-m\Delta)}.\end{split}$$

$$\begin{split}\mathbf{h}_{c^{k}}[n]&=h_{c}(t-k\Delta_{K})\Big{|}_{t=n\Delta}\\ &=\sum_{m=-\infty}^{m=\infty}\mathbf{h}_{c}[m]\frac{\sin\pi((n-m)\Delta-k\Delta_{K})}{\pi((n-m)\Delta-k\Delta_{K})}\\ &=\mathbf{h}_{c}*\underbrace{\Bigg{(}\frac{\sin\pi(n\Delta-k\Delta_{K})}{\pi(n\Delta-k\Delta_{K})}\Bigg{)}}_{\mathbf{f}^{k}}.\\ \end{split}$$

$$\mathbf{F}^{k}=\begin{pmatrix}\mathbf{f}^{k}[0]&\mathbf{f}^{k}[-1]&\cdots&\mathbf{f}^{k}[-\frac{L-1}{2}]&\mathbf{0}_{\big{(}\frac{L-1}{2}\big{)}}^{T}\\ \mathbf{f}^{k}[1]&\mathbf{f}^{k}[0]&\cdots&\mathbf{f}^{k}[-\frac{L-1}{2}]&\mathbf{0}_{\big{(}\frac{L-3}{2}\big{)}}^{T}\\ &&\vdots&&\\ \mathbf{f}^{k}[\frac{L-1}{2}]&&\cdots&&\mathbf{f}^{k}[-\frac{L-1}{2}]\\ &&\vdots&&\\ \mathbf{0}_{\big{(}\frac{L-1}{2}\big{)}}^{T}&\mathbf{f}^{k}[\frac{L-1}{2}]&\cdots&&\mathbf{f}^{k}[0]\\ \end{pmatrix}.$$

$$\begin{split}\mathbf{H}^{\text{T}}\mathbf{r}^{(t^{\prime})}=&\big{[}(\mathbf{h}_{1}\star\mathbf{r}^{(t^{\prime})})[1],\cdots,(\mathbf{h}_{1}\star\mathbf{r}^{(t^{\prime})})[N-L+1],\cdots,\\ &(\mathbf{h}_{C}\star\mathbf{r}^{(t^{\prime})})[1],\cdots,(\mathbf{h}_{C}\star\mathbf{r}^{(t^{\prime})})[N-L+1]\big{]}^{\text{T}},\end{split}$$

$$\begin{split}&\mathbf{L}^{(t^{\prime})}=\begin{pmatrix}\mathbf{L}^{(t^{\prime}-1)}&0\\ \mathbf{w}^{\text{T}}&\sqrt{1-\lVert\mathbf{w}\rVert_{2}^{2}}\\ \end{pmatrix}\\ &\text{where }\mathbf{L}^{(t^{\prime}-1)}\mathbf{w}=\mathbf{H}\big{|}_{t^{\prime}-1}^{\text{T}}\mathbf{h}^{(t^{\prime})},\\ \end{split}$$

$$\mathbf{H}{\mathbf{X}}_j=\sum _{c=1}^C\sum _{i=1}^{N_c}{x}_{j,i}^c{\mathbf{S}}_{j,i}^c{\mathbf{h}}_c,$$

$${\mathbf{S}}_{j,i}^c=\left(\begin{array}{c}{0}_{n_{j,i}^c\times L}\\ {\mathbf{I}}_{L\times L}\\ 0_{(W-L-n_{j,i}^c)\times L}\end{array}\right).$$

$$\min_{\{\mathbf{h}_{c}\}_{c=1}^{C}}\sum_{j=1}^{J}\Big{\lVert}\mathbf{Y}_{j}-\sum_{c=1}^{C}\sum_{i=1}^{N_{c}}x_{j,i}^{c}\mathbf{S}_{j,i}^{c}\mathbf{h}_{c}\Big{\rVert}_{2}^{2}.$$

$$\begin{split}&\min_{\{\mathbf{h}_{c}\}_{c=1}^{C}}\sum_{j=1}^{J}\Big{\lVert}\mathbf{Y}_{j}-\sum_{c=1}^{C}\sum_{k=1}^{K}\sum_{i=1}^{N_{c^{k}}}x_{j,i}^{c^{k}}\mathbf{S}_{j,i}^{c^{k}}\mathbf{h}_{c^{k}}\Big{\rVert}_{2}^{2}\\ &\text{s.t. }\mathbf{h}_{c^{k}}=\mathbf{F}^{k}\mathbf{h}_{c},\quad\forall k,c,\\ \end{split}$$

$$\begin{split}&\widehat{\mathbf{h}}_{\overline{c}}=\arg\min_{\mathbf{h}_{\overline{c}}}\sum_{j\in\Omega_{\overline{c}}}\Bigg{\lVert}\mathbf{E}_{j}-\sum_{i=1}^{N_{\overline{c}}}x_{j,i}^{\overline{c}}\mathbf{S}_{j,i}^{\overline{c}}\mathbf{h}_{\overline{c}}\Bigg{\rVert}_{2}^{2},\\ &\text{where }\mathbf{E}_{j}=\mathbf{Y}_{j}-\sum_{c\neq\overline{c}}^{C}\sum_{i=1}^{N_{c}}x_{j,i}^{c}\mathbf{S}_{j,i}^{c}\mathbf{h}_{c},\\ \end{split}$$

$$\begin{split}\widehat{\mathbf{h}}_{\overline{c}}&=\Bigg{(}\sum_{j\in\Omega_{\overline{c}}}\sum_{i=1}^{N_{\overline{c}}}\sum_{m=1}^{N_{\overline{c}}}x_{j,i}^{\overline{c}}\Big{(}\mathbf{S}_{j,i}^{\overline{c}}\Big{)}^{\text{T}}\mathbf{S}_{j,m}^{\overline{c}}x_{j,m}^{\overline{c}}\Bigg{)}^{-1}\\ &\times\Bigg{(}\sum_{j\in\Omega_{\overline{c}}}\sum_{i=1}^{N_{\overline{c}}}x_{j,i}^{\overline{c}}\Big{(}\mathbf{S}_{j,i}^{\overline{c}}\Big{)}^{\text{T}}\mathbf{E}_{j}\Bigg{)}.\\ \end{split}$$

$$\begin{array}{r}{\mathbf{S}}_{j,i}^{c^k}{\mathbf{h}}_{c^k}={\mathbf{S}}_{j,i}^{c^k}{\mathbf{F}}^k{\mathbf{h}}_c={\tilde{\mathbf{S}}}_{j,i}^{c^k}{\mathbf{h}}_c,\end{array}$$

$$\min_{\{\mathbf{h}_{c}\}_{c=1}^{C}}\sum_{j=1}^{J}\Big{\lVert}\mathbf{Y}_{j}-\sum_{c=1}^{C}\sum_{k=1}^{K}\sum_{i=1}^{N_{c^{k}}}x_{j,i}^{c^{k}}\widetilde{\mathbf{S}}_{j,i}^{c^{k}}\mathbf{h}_{c}\Big{\rVert}_{2}^{2}.$$

$$\begin{split}&\widehat{\mathbf{h}}_{\overline{c}}=\arg\min_{\mathbf{h}_{\overline{c}}}\sum_{j\in\Omega_{\overline{c}}}\Bigg{\lVert}\mathbf{E}_{j}-\sum_{k=1}^{K}\sum_{i=1}^{N_{\overline{c}^{k}}}x_{j,i}^{\overline{c}^{k}}\widetilde{\mathbf{S}}_{j,i}^{\overline{c}^{k}}\mathbf{h}_{\overline{c}}\Bigg{\rVert}_{2}^{2}\\ &\text{where }\mathbf{E}_{j}=\mathbf{Y}_{j}-\sum_{c\neq\overline{c}}^{C}\sum_{k=1}^{K}\sum_{i=1}^{N_{c^{k}}}x_{j,i}^{c^{k}}\widetilde{\mathbf{S}}_{j,i}^{c^{k}}\mathbf{h}_{c},\\ \end{split}$$

$$\begin{split}\widehat{\mathbf{h}}_{\overline{c}}&=\Bigg{(}\sum_{j\in\widetilde{\Omega}_{\overline{c}}}\sum_{k=1}^{K}\sum_{i=1}^{N_{\overline{c}}}\sum_{m=1}^{N_{\overline{c}}}x_{j,i}^{\overline{c}^{k}}\Big{(}\widetilde{\mathbf{S}}_{j,i}^{\overline{c}^{k}}\Big{)}^{\text{T}}\widetilde{\mathbf{S}}_{j,m}^{\overline{c}^{k}}x_{j,m}^{\overline{c}^{k}}\Bigg{)}^{-1}\\ &\times\Bigg{(}\sum_{j\in\widetilde{\Omega}_{\overline{c}}}\sum_{k=1}^{K}\sum_{i=1}^{N_{\overline{c}}}x_{j,i}^{\overline{c}^{k}}\Big{(}\widetilde{\mathbf{S}}_{j,i}^{\overline{c}^{k}}\Big{)}^{\text{T}}\mathbf{E}_{j}\Bigg{)}.\\ \end{split}$$

$$\begin{split}&h_{1}(t)\propto\big{(}10^{3}t\big{)}\exp\Big{(}-\big{(}10^{3}t\big{)}^{2}\Big{)}\cos\Big{(}\frac{\pi}{2}\big{(}10^{3}t\big{)}\Big{)}\\ &h_{2}(t)\propto\big{(}10^{3}t\big{)}\exp\Big{(}-\big{(}10^{3}t\big{)}^{2}\Big{)}.\\ \end{split}$$

$$\text{SNR}=10\log \left(\frac{\Vert \mathbf{Hx}{\Vert }_2^2}{\Vert \epsilon {\Vert }_2^2}\right)\text{(dB)},$$

$$\text{err}({\hat{\mathbf{h}}}_c,{\tilde{\mathbf{h}}}_c)=\sqrt{1-\frac{⟨{\hat{\mathbf{h}}}_c,{\tilde{\mathbf{h}}}_c{⟩}^2}{\Vert {\hat{\mathbf{h}}}_c{\Vert }_2^2\Vert {\tilde{\mathbf{h}}}_c{\Vert }_2^2}}.$$

Full text

I Introduction

In recent years, the problem of decomposing an observed signal into a sparse linear combination of elements drawn from a known dictionary, often referred to as sparse approximation [1], has been of great interest to the signal processing community. Specifically, representing the signal as the superposition of time-shifted (or shift-invariant) templates with local support has received special attention [2]. This is due to the observation that many examples of real-world signals can be modeled in this manner. For instance, signals arising from electrophysiological recordings of neural activity can be modeled as the sum of distinct action potentials produced by the neurons near recording electrodes [3]. In studies involving electroencephalography (EEG) or magnetoencephalography (MEG), there is growing evidence that the signal should be studied in terms of the aggregation of transient events with specific templates [4].

A generative model for these signals is the convolution in continuous time between templates, a collection of which is referred to as dictionary, and a set of scaled and time-shifted delta functions, referred to as codes. The amplitude and the location of each code correspond to the magnitude and the time when an event occurs, respectively. Given an observed signal, the goal of Convolutional Dictionary Learning (CDL) frameworks is to estimate the templates and the codes under the generative model, and with sparsity constraints on the codes. These frameworks typically alternate between two steps, a Convolutional Sparse Coding (CSC) step to estimate the codes, and a Convolutional Dictionary Update (CDU) step to estimate the dictionary.

One of the drawbacks of existing CDL frameworks is the assumption that the signal of interest occurs in discrete time, when in fact the underlying signal occurs in continuous time. The discrete-time approximation of the generative model introduces errors known as time-quantization errors [5]. Specifically, if an event in the continuous-time model were to occur at a time that does not coincide with any point on the discrete-time sampling grid, the CSC step would inaccurately identify the event as occurring at a time on the grid. Increasing the sampling rate to obtain a finer grid is a viable solution, but is not always possible for two reasons. First, practically, the sampling rate is limited by hardware specifications and the amount of memory available for storage. Second, the large number of samples will increase the computation time of the CSC and CDU steps.

Continuous Basis Pursuit (CBP) [6], a convex sparse regression framework, was introduced to address such limitations of the conventional approaches for CSC. The CSC step of CBP begins with the derivation, from the original dictionary, of an alternate dictionary that approximates the subspace of continuous-time-shifted copies of the templates (the continuous-time shifts are not necessarily multiples of the sampling interval). Then, using the new dictionary, it solves a $\ell_{1}$-regularized convex regression problem. Despite the method being more accurate, compared to discrete-time approaches, in estimating the times when events occur, it does not scale well with the size of modern datasets. Moreover, because CBP does not include a CDU step, it is not a CDL framework.

With the increasing ability to record larger datasets, there is a need for an efficient CDL framework that addresses the time-quantization issue. As our primary motivation is efficiency, we focus on extending CDL frameworks based on greedy methods, which are known to be less computationally demanding than basis pursuit and $\ell_{1}$-regularized methods [7], [8], [9]. To address the time-quantization issue, we build an overcomplete and interpolated dictionary that accounts for events that occur off the discrete-time sampling grid. Unlike in CBP, this approach also allows us to learn the templates while incorporating the amplitudes and the times of events that occur off the grid.

Our contributions are the following:

• •

A fast convolutional greedy pursuit algorithm We introduce an efficient algorithm for convolutional greedy pursuit under a discrete-time generative model. When the events from the continuous-time model occur on the sampling grid, we term this algorithm Convolutional Orthogonal Matching Pursuit (COMP). COMP is much faster than the alternative based on the basis-pursuit algorithm. Compared to Convolutional Matching Pursuit (CMP), the computational overhead from COMP is negligible (Table III, IV).

• •

A CSC framework that handles non-integer delays We introduce a discrete-time generative model that accounts for events from the continuous-time model that do not occur on the discrete-time sampling grid. This model is inspired by the concept of bandlimited interpolation in digital signal processing. To perform the CSC, we extend COMP and call the resulting algorithm COMP-INTERP. Compared to conventional CSC frameworks, COMP-INTERP is more accurate in identifying the times when of off-the-grid events occur and achieves an accuracy similar to CBP (Figure 5). Moreover, COMP-INTERP is much faster than CBP (Table III, IV).

• •

A dictionary update framework that handles non-integer delays For the first time, we introduce a CDU algorithm that accounts for estimates of the sparse codes from the CSC step that correspond to events off the grid. When compared to conventional dictionary update algorithms, the templates learned from our approach are more accurate (Figure 6). In addition, our algorithm is more robust to varying noise levels (Figure 7).

The rest of our treatment begins in Section II, where we introduce the generative model of interest and formulate the CDL objective functions. In Section III, we review existing work relevant to our CDL framework. In Section IV, we introduce an efficient framework, termed COMP-INTERP, for performing CSC while accounting for events that occur off the grid. In Section V, we introduce a CDU step that can handle events off the grid. We used simulated and real datasets to demonstrate the performance of our algorithms in Section VI. We conclude in Section VII.

II Generative Model and Problem setup

II-A Notations

Table I summarizes our notational conventions. We introduce additional notation as necessary, at the beginning of the section that uses it first. We use the expressions event off the grid and event with a non-integer delay interchangeably.

II-B Continuous and discrete-time generative models

Let $y(t)$ be an observed continuous signal in interval $(0,T]$ and $\{h_{c}(t)\}_{c=1}^{C}$ be templates (filters) from $C$ sources. We assume that the templates each have the same length and are localized in time. The shift-invariant continuous generative model expresses $y(t)$ as follows

[TABLE]

where $N_{c}$ denotes the number of events from source $c$, $\tau_{i}^{c}$ and $x_{i}^{c}$ denote the position and the amplitude of the $i^{\text{th}}$ event from source $c$, respectively. The variable $\varepsilon(t)$ denotes i.i.d. white noise.

To formulate a discrete-time analogue of Eq. 1, let $\Delta$ denote the length of a sampling interval and $f_{s}=\frac{1}{\Delta}$ the associated sampling frequency. The number of intervals of size $\Delta$ in $(0,T]$ is $N=\lfloor\frac{T}{\Delta}\rfloor$. Further let $n=1,\cdots,N\in\mathbb{N}^{+}$ be the discrete-time index, and $n_{i}^{c}$ denote the discrete-time approximation of $\tau_{i}^{c}$, which is such that $n_{i}^{c}\Delta\leq\tau_{i}^{c}<(n_{i}^{c}+1)\Delta$. Finally, we denote by $\mathbf{h}_{c}\in\mathbb{R}^{L}$ the discrete time analogue of $h_{c}(t)$, which we assume is normalized such that $\lVert\mathbf{h}_{c}\rVert_{2}=1,\forall c$. Using this notation, we can obtain discrete-time samples $y[n]=y(n\Delta)$ of $y(t)$ that satisfy

[TABLE]

where we refer to $\mathbf{x}^{c}[n]=\sum_{i=1}^{N_{c}}x_{i}^{c}\delta[n-n_{i}^{c}]$, for $n=1,\cdots,N-L+1$, as the code, and \mathbf{x}^{c}=\big{[}\mathbf{x}^{c}[1],\cdots,\mathbf{x}^{c}[N-L+1]\big{]}^{\text{T}} as the code vector.

We can express Eq. 2 in linear-algebraic form as follows

[TABLE]

where \mathbf{H}=\begin{bmatrix}\mathbf{H}^{1}\Big{|}\cdots\Big{|}\mathbf{H}^{C}\end{bmatrix}\in\mathbb{R}^{N\times C(N-L+1)} is a block-Toeplitz matrix with $c^{\text{th}}$ block $\mathbf{H}^{c}\in\mathbb{R}^{N\times(N-L+1)}$ for $c=1,\cdots,C$ and $\mathbf{x}=[(\mathbf{x}^{1})^{\text{T}},\cdots,(\mathbf{x}^{C})^{\text{T}}]^{\text{T}}\in\mathbb{R}^{C(N-L+1)}$. The columns of the Toeplitz matrix $\mathbf{H}^{c}$ represent delayed versions (time-shifts) of $\mathbf{h}_{c}$, with integer delay between [math] and $N-L$, that have been zero-padded to have length $N$. For each $c$, the non-zero entries of $\mathbf{x}^{c}$ represent the discrete-time indices $\{n_{i}^{c}\}_{i=1}^{N_{c}}$ when source $c$ appears in the signal $y[n]$. Fig. 1 illustrates Eq. 3.

In practice, we divide the signal $\mathbf{y}$ into $J$ non-overlapping windows, each of length $W$ and such that $N=WJ$. We assume that $L<<W<<N$, so that the filters from each source are localized within the signal $y[n]$ and within a window. We denote by $\mathbf{Y}\in\mathbb{R}^{W\times J}$ the matrix whose $j^{\text{th}}$ column is \mathbf{Y}_{j}=\big{[}\mathbf{y}[(j-1)W+1],\cdots,\mathbf{y}[jW]\big{]}^{\text{T}}\in\mathbb{R}^{W}, namely the $j^{\text{th}}$ window from $y[n]$. Similarly, we denote by $\mathbf{X}\in\mathbb{R}^{C(W-L+1)\times J}$ the coefficient matrix whose $j^{\text{th}}$ column $\mathbf{X}_{j}=[(\mathbf{x}_{j}^{1})^{T},\cdots,(\mathbf{x}_{j}^{C})^{\text{T}}]^{\text{T}}\in\mathbb{R}^{C(W-L+1)}$ is the code vector associated with window $j$.

II-C Optimization objective

The goal of Convolutional Dictionary Learning (CDL) is to estimate $\{\mathbf{h}_{c}\}_{c=1}^{C}$ and $\{\mathbf{X}_{j}\}_{j=1}^{J}$ that minimize the error of reconstructing the $\mathbf{Y}_{j}$ in each window using its linear approximation $\mathbf{H}\mathbf{X}_{j}$. We impose a sparsity constraint on the total number of nonzero elements of $\{\mathbf{X}_{j}\}_{j=1}^{J}$ for two reasons. First, without additional constraints, the problem as posed leads to an under-determined system of equations that does not have a unique solution. Second, in many applications, the rate of occurrence of events from the the sources of interest is small compared to $T$, implying that each block $\mathbf{x}_{j}^{c}$ of the vector $\mathbf{X}_{j}$ is sparse. For example, in electrophysiological recordings of neural activity, we expect a sparse number the action potentials from neurons due to their biophysical properties [10]. Following [11], we use the $\ell_{0}$ quasi-norm $\lVert\mathbf{X}_{j}\rVert_{0}$, which counts the number of non-zero elements of a vector, to express the sparsity constraint. This lets us express the problem mathematically as follows

[TABLE]

where $\beta_{1}$ is a pre-defined sparsity threshold. We refer to this optimization as the CDL problem. One limitation of this approach comes from approximating the continuous-domain generative model (Eq. 1) with the discrete-domain generative model (Eq. 2). This approximation results in time-quantization errors, which manifest themselves in two ways: 1) a mismatch between the time when the event occurs in continuous time and its approximation in discrete time, $\tau^{c}_{j,i}\neq n^{c}_{j,i}\Delta$ and 2) template mismatch, $\mathbf{h}_{c}[m-n^{c}_{j,i}]=h_{c}(m\Delta-n^{c}_{j,i}\Delta)\neq h_{c}(m\Delta-\tau^{c}_{j,i})$. One of our contributions is to introduce, in Section IV, a discrete-time generative model that mitigates the effects of time-quantization errors.

CDL is a nonconvex optimization problem, due to the simultaneous optimization over $\{\mathbf{h}_{c}\}_{c=1}^{C}$ and $\{\mathbf{X}_{j}\}_{j=1}^{J}$, as well as the $\ell_{0}$ penalty. A popular approach is to alternatively minimize the objective over one set of variables while the other is fixed, until convergence. At iteration $t+1$, $\mathbf{X}^{(t+1)}$ is computed based on $\mathbf{H}^{(t)}$ through a sparse coding step, after which $\mathbf{H}^{(t+1)}$ is computed using $\mathbf{X}^{(t+1)}$ through a dictionary update step. If $\mathbf{H}$ is a convolutional matrix, we refer to these steps as Convolutional Sparse Coding (CSC) and Convolutional Dictionary Update (CDU), respectively. Existing CSC approaches fall into two categories, which differ with respect to how they enforce the sparsity constraint. One class of approaches, which is the one we follow in this work, uses greedy methods to solve the original problem with the $\ell_{0}$ quasi-norm. Another class of approaches relaxes the $\ell_{0}$ quasi-norm to the $\ell_{1}$ norm, which converts the CSC objective into a convex optimization problem [12], [13], [14]. The advantage of greedy approaches is that they are more efficient computationally [7]. Existing CSC frameworks that address time-quantization errors use the $\ell_{1}$ norm to enforce sparsity [6]. In the next section, we review both classes of CSC approaches, as well as approaches to solve the CDU step.

III Background for CDL

For notational simplicity, we use $\mathbf{x}$ instead of $\mathbf{X}_{j}$ and $\mathbf{y}$ instead of $\mathbf{Y}_{j}$.

III-A Convolutional Greedy Pursuit

Matching Pursuit (MP) [15] and Orthogonal Matching Pursuit (OMP) [16] are greedy methods to solve the CSC step. We introduce the methods first when $\mathbf{H}$ is an arbitrary matrix, and then discuss the convolutional case.

III-A1 Classical greedy pursuit - MP and OMP

Both MP and OMP iteratively select columns from $\mathbf{H}$ to produce an approximation $\mathbf{H}\mathbf{x}$ of $\mathbf{y}$. At iteration step $t^{\prime}+1$, the column of $\mathbf{H}$ with the maximal absolute inner product with the residual $\mathbf{r}^{(t^{\prime})}$ is selected and added to the active filter index set $\mathcal{S}^{(t^{\prime})}$. The initial conditions are $\mathbf{r}^{(0)}=\mathbf{y}$ and $\mathcal{S}^{(0)}=\{\emptyset\}$. The two methods differ in how the coefficients of the chosen columns and the residuals are computed. Let $\mathbf{h}^{(t^{\prime})}$ denote the template chosen at iteration $t^{\prime}$.

• •

MP The coefficient associated with $\mathbf{h}^{(t^{\prime})}$ and the residual are given, respectively, by $\langle\mathbf{h}^{(t^{\prime})},\mathbf{r}^{(t^{\prime})}\rangle$ and $\mathbf{r}^{(t^{\prime}+1)}=\mathbf{r}^{(t^{\prime})}-\langle\mathbf{h}^{(t^{\prime})},\mathbf{r}^{(t^{\prime})}\rangle\mathbf{h}^{(t^{\prime})}$.

• •

OMP The coefficients associated with $\mathbf{h}^{(1)},\cdots\mathbf{h}^{(t^{\prime})}$ are those that minimize the squared error between $\mathbf{y}$ and its linear reconstruction using the columns, i.e. those obtained by least-squares. This is equivalent to projecting $\mathbf{y}$ onto the span of $\mathbf{h}^{(1)},\cdots,\mathbf{h}^{(t^{\prime})}$, and is called a projection step.

The projection step implies two key differences between MP and OMP. First, OMP is slower than MP, due to the matrix inversion and multiplication required in the former. Second, as the residual $\mathbf{r}^{(t^{\prime})}$ in OMP is orthogonal to the span of previously selected columns, a different column of $\mathbf{H}$ is selected at every iteration. This is not the case for MP, which means that the same column can be selected multiple times.

III-A2 Convolutional extensions

CMP has enjoyed popularity in biomedical applications [17, 18] and image recognition [19]. Convolutional OMP, however, has been less popular, primarily because the projection step with a convolutional matrix is computationally expensive. A previous work [20] suggests a way to improve efficiency, at the expense of accuracy, by projecting the residual onto a subset of the active columns rather than all of them.

III-B CSC for estimating continuous time-shifts

Continuous Basis Pursuit (CBP) CBP is a CSC framework that uses the $\ell_{1}$ norm to enforce sparsity, and addresses time-quantization errors [6]. Let $\tau_{i}^{c}\in[(m-\frac{1}{2})\Delta,(m+\frac{1}{2})\Delta)$ with $m\in\mathbb{N}$, and h_{c}\big{(}t-\tau_{i}^{c}\big{)} be a continuous-time shift of $h_{c}(t)$. CBP first uses a local basis, in combination with integer shifts $\{h_{c}(t-m^{\prime}\Delta)\}_{m^{\prime}}$ of $h_{c}(t)$ around $\tau_{i}^{c}$, to define a new set of functions $\{\boldsymbol{\phi}^{c}_{p}\}_{p=1}^{P}\in\mathbb{R}^{L}$ whose linear combination approximates the continuous-time shift. Both $P$ and the number of integer shifts $m^{\prime}$ depend on the choice of basis. A popular option is the Polar basis, where the set of approximating functions are derived from trigonometric splines, with $P=3$ and $m^{\prime}\in[m,m\pm 1]$. Let $\{\mathbf{\Phi}^{c,p}\}_{p=1}^{P}\in\mathbb{R}^{N\times(N-L+1)}$ be the convolutional dictionary constructed from the basis and $\mathbf{x}^{c}_{p}\in\mathbb{R}^{N-L+1}$ a sparse code vector associated with $\mathbf{\Phi}^{c,p}$. Then, it solves the convex $\ell_{1}$-regularized optimization problem

[TABLE]

where $\lVert\mathbf{x}^{c}_{1}\rVert_{1}$ is the $\ell_{1}$ norm of the amplitude of the first basis element and $\mathcal{H}$ is convex constraint set that depends on $\{\boldsymbol{\phi}_{p}^{c}\}_{c,p=1}^{C,P}$. Lastly, the code $\{\mathbf{x}^{c}_{p}\}_{c,p=1}^{C,P}$ is converted to coefficients for the original basis $\{\mathbf{h}_{c}\}_{c=1}^{C}$ through an inverse mapping, which yield the times when the events occur and their amplitudes. Fig. 2 shows an application of CBP, BP, COMP and COMP-INTERP (COMP with interpolated dictionary) to the estimation of the continuous-time shift and amplitude of a single event from one filter ($C=1$). CBP and COMP-INTERP, the approach we propose in Section IV, are able to estimate the continuous-time shift accurately. As we demonstrate in Section VI, COMP-INTERP is orders of magnitude faster than CBP. BP and COMP cannot capture the continuous-time shift.

Continuous OMP In contrast with CBP that uses the $\ell_{1}$ norm, Continuous OMP [21] is a greedy method to solve Eq. 5. It uses an extension of OMP to deal with the $\ell_{0}$ norm. Being greedy, continuous OMP is faster than CBP. Its downside is that it requires the solution to a large number of optimization problems both in the selection and the projection step.

III-C CDU frameworks

The majority of existing CDU frameworks estimate the templates $\{\mathbf{h}_{c}\}_{c=1}^{C}$ by minimizing the error of reconstructing $\mathbf{y}$ using its linear approximation $\mathbf{H}\mathbf{x}$. The key differences between existing approaches are the constraints imposed on the templates and the optimization methods used, as detailed in a recent survey [12]. To the best of our knowledge, existing CDU approaches do not address the problem of learning the templates in the presence of time-quantization errors.

IV Convolutional Orthogonal Matching Pursuit with Interpolated Dictionary

For the CSC step in the alternating-minimization approach to CDL, we introduce an algorithm for off-the-grid sparse coding called Convolutional OMP (COMP) with interpolated dictionary (COMP-INTERP). This is a convolutional greedy pursuit method that minimizes time-quantization errors. From a computational perspective, COMP-INTERP is an efficient alternative to the $\ell_{1}$-based CSC frameworks, such as CBP.

We use $\mathbf{x}$ and $\mathbf{y}$, instead of $\mathbf{X}_{j}$ and $\mathbf{Y}_{j}$, for notational simplicity. We use $t$ and $t^{\prime}$ to denote alternating-minimization iteration and COMP inner iteration, respectively. Since the discussion in this section involves a single iteration of alternating-minimization procedure, we drop $t$.

IV-A Non-integer delay through smooth interpolation

The discrete-time model from Eq. 2 is restrictive because events from the continuous-time generative model of Eq. 1 do not necessarily occur at multiples of the sampling interval $\Delta$. We address this limitation by partitioning $\Delta$ into finer intervals of length $\Delta_{K}:=\frac{1}{K}\Delta$, and modifying $\mathbf{H}$ and $\mathbf{x}$ accordingly. The resulting CSC framework, with finer resolution $\Delta_{K}$, can approximate the time of occurrence of the event $\tau_{i}^{c}$ with $m\Delta+k\Delta_{K}$, where $m\in\mathbb{N}$ and $k=0,\cdots,K-1$, rather than with $m\Delta$, which leads to reduction of time quantization error. That is, $|\tau_{i}^{c}-m\Delta|\geq|\tau_{i}^{c}-(m\Delta+k\Delta_{K})|$. By definition, each template $\mathbf{h}_{c}$ corresponds to discrete-time samples of $h_{c}(t)$ with resolution $\Delta$. Our challenge is to modify $\mathbf{H}$ to account for versions of $h_{c}(t)$ delayed by a non-integer amount $k\Delta_{K}$, that is not an integer multiple of $\Delta$, and sampled at resolution $\Delta$. Let $\mathbf{h}_{c^{k}}\in\mathbb{R}^{L}$ denote $h_{c}(t)$ delayed by a non-integer amount $k\Delta_{K}$ and sampled at resolution $\Delta$. This definition motivates us to reformulate Eq. 2 to account for non-integer shifts of the templates at a finer scale $\Delta_{K}$

[TABLE]

where $\mathbf{x}^{c^{k}}$ denotes the code vector corresponding to $\mathbf{h}_{c^{k}}$. For notational simplicity, we let $\mathbf{h}_{c^{0}}=\mathbf{h}_{c}$. Note that $\mathbf{h}_{c}\neq\mathbf{h}_{c^{k}}$ for $k\neq 0$, as illustrated in Fig. 3. The systematic method of obtaining $\mathbf{h}_{c^{k}}$ from $\mathbf{h}_{c}$ will be discussed in the next section. We use $\{\mathbf{h}_{c^{k}}\}_{c,k=1}^{C,K}$ to construct the interpolated convolutional dictionary $\mathbf{\widetilde{H}}\in\mathbb{R}^{W\times CK(W-L+1)}$

[TABLE]

where $\mathbf{\widetilde{H}}^{c^{k}}$ is the Toeplitz matrix whose columns consists of all integer shifts of $\mathbf{h}_{c^{k}}$. Note that when $K=1$, we get the original convolutional dictionary, i.e. $\widetilde{\mathbf{H}}=\mathbf{H}$. In linear-algebraic form, we can write the generative model as $\mathbf{Y}=\mathbf{\widetilde{H}}\mathbf{X}+\boldsymbol{\varepsilon}$, where $\mathbf{X}\in\mathbb{R}^{CK(W-L+1)\times J}$.

IV-B Smooth interpolation of $\mathbf{h}_{c}$

We use the concept of continuous-time operations on discrete-time signals from digital signal processing theory [22] to obtain $\mathbf{h}_{c^{k}}$ from $\mathbf{h}_{c}$. The process consists of three steps: interpolation, shifting by a non-integer amount, and resampling. These steps, illustrated in Fig. 3, perform the following operations

Interpolation Interpolate $\mathbf{h}_{c}$ with a smooth interpolator to obtain $\tilde{h}_{c}(t)$. 2. 2.

Non-integer shift Shift $\tilde{h}_{c}(t)$ to obtain $\tilde{h}_{c}(t-k\Delta_{K})$. 3. 3.

Resampling Resample $\tilde{h}_{c}(t-k\Delta_{K})$ with resolution $\Delta$ to obtain $\mathbf{h}_{c^{k}}$.

The interpolated template $\tilde{h}_{c}(t)$ depends on the choice of interpolator and does not necessarily coincide with the continuous-time template $h_{c}(t)$. In the next section, we show that a sinc interpolator yields $\tilde{h}_{c}(t)=h_{c}(t)$.

IV-C Bandlimited interpolation of $\mathbf{h}_{c}$

Suppose $h_{c}(t)$ is bandlimited and that the sampling frequency $f_{s}$ is above its Nyquist rate. The sinc interpolator then guarantees that $\tilde{h}_{c}(t)=h_{c}(t)$, namely that we can recover the continuous-time filter from its discrete-time samples. Approximating a continuous-time signal with discrete-time samples in this manner is commonly referred to as bandlimited interpolation [22]. The intuition behind bandlimited interpolation is two-fold: 1) convolving a signal with a sinc interpolator is equivalent to multiplying the signal with an ideal low-pass filter with unit gain in the frequency domain, and thus the signal is not distorted upon interpolation, and 2) when bandlimited $h_{c}(t)$ is sampled at $f_{s}$, there is no aliasing of high frequency content into low frequency bands. The bandlimited assumption translates to the smoothness on $h_{c}(t)$, which requires that the template does not change abruptly between consecutive samples on the sampling grid. Interpolating $\mathbf{h}_{c}$ with a sinc interpolator $\text{sinc}(t)=\frac{\sin\pi t}{\pi t}$ yields

[TABLE]

To obtain $\mathbf{h}_{c^{k}}$, we shift $\tilde{h}_{c}(t)=h_{c}(t)$ by $k\Delta_{K}$ and resample with resolution $\Delta$, where $\mathbf{f}^{k}$ is defined as the discrete sinc interpolator shifted by $k\Delta_{K}$,

[TABLE]

The three-step procedure described previously, and illustrated in Fig. 3, is equivalent to convolving $\mathbf{h}_{c}$ with $\mathbf{f}^{k}$. In practice, since the sinc interpolator has infinite support, we use a truncated sinc interpolator to ensure that $\mathbf{h}_{c^{k}}$ is of finite length. The truncation provides a good approximation to the infinite discrete sinc interpolator [22]. Let $\mathbf{f}^{k}\in\mathbb{R}^{L}$ be the sinc interpolator restricted to the support \big{[}-\frac{L-1}{2},\frac{L-1}{2}\big{]}. We obtain the truncated $\mathbf{h}_{c^{k}}$ by expressing the convolution as a matrix multiplication: $\mathbf{h}_{c^{k}}=\mathbf{F}^{k}\mathbf{h}_{c}$, where $\mathbf{F}^{k}\in\mathbb{R}^{L\times L}$ is defined as

[TABLE]

The matrix $\mathbf{F}^{k}$ is the Toeplitz (convolutional) matrix associated with $\mathbf{f}^{k}$, truncated so that it is of $L\times L$. The truncation ensures that the resulting $\mathbf{h}_{c^{k}}$ is of the same length as $\mathbf{h}_{c}$.

IV-D Efficient algorithm for COMP and COMP-INTERP

Matrix operations involving convolutional dictionaries $\mathbf{H}$ or $\widetilde{\mathbf{H}}$ are expensive both in terms of computations and storage requirements. This is because typical recordings can last on the order of minutes, if not hours, and sampling rates can be on the order of $\sim 10^{4}(Hz)$ for electrophysiology and $10^{3}(Hz)$ for EEG, to name a few examples. Existing greedy algorithms for CSC can handle high-dimensional data in the selection step because it involves cross-correlation operations for which efficient and scalable implementations exist. They fall short, however, of addressing the inefficiency of the projection step in high-dimensional settings. This has led to the wide adoption of CMP and not COMP.

We explore efficient implementations of the selection and the projection step for a convolutional matrix $\mathbf{H}$ (or $\widetilde{\mathbf{H}}$), which lead to efficient implementations of COMP and COMP-INTERP. An efficient implementation of the selection step, which we restate for completeness, was proposed in [20]. For notational simplicity, we focus on $\mathbf{H}$, noting that the same argument holds for $\mathbf{\widetilde{H}}$.

Selection step The selection step requires the inner product between time-shifted $\mathbf{h}_{c}$ and $\mathbf{r}^{(t^{\prime})}$, expressed as $\mathbf{H}^{\text{T}}\mathbf{r}^{(t^{\prime})}$. For large $\mathbf{H}$, explicit computation of the inner product is expensive. However, the convolutional structure of $\mathbf{H}$ lets compute $C$ cross-correlations instead [23]

[TABLE]

where the cross-correlation index $i$ indicates the time offset of $\mathbf{h}_{c}[1]$ relative to the $\mathbf{r}^{(t^{\prime})}[1]$. This formulation has two benefits:

We do not need to construct the convolutional matrix $\mathbf{H}$ explicitly. We only require $O(N)$ memory to store $\{\mathbf{h}_{c}\}_{c=1}^{C}$ and $\mathbf{r}^{(t^{\prime})}$, as opposed to $O(CN^{2})$ memory to store the matrix. 2. 2.

We can compute the $C$ cross-correlation operations using the FFT, which is much more efficient than computing them by multiplication of $\mathbf{H}^{\text{T}}$ and $\mathbf{r}^{(t^{\prime})}$.

Projection step In this step, we project the residual onto the span of \mathbf{H}\big{|}_{t^{\prime}}, which requires the inversion of \mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{H}\big{|}_{t^{\prime}}. The matrix \mathbf{H}\big{|}_{t^{\prime}}\in\mathbb{R}^{N\times t^{\prime}} refers to a convolutional dictionary restricted to columns that have been selected by COMP up to $t^{\prime\text{th}}$ iteration. Consequently, the code $\mathbf{x}|_{t^{\prime}}\in\mathbb{R}^{t^{\prime}}$ refers to the nonzero coefficients from $\mathbf{x}\in\mathbb{R}^{N-L+1}$ corresponding to the columns \mathbf{H}\big{|}_{t^{\prime}}. To avoid the computational cost of inversion in the projection step,  [24] suggested an efficient method for computing the Cholesky factor $\mathbf{L}^{(t^{\prime})}$, which is a lower triangular matrix such that \mathbf{L}^{(t^{\prime})}\big{(}\mathbf{L}^{(t^{\prime})}\big{)}^{\text{T}}=\mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{H}\big{|}_{t^{\prime}}\in\mathbb{R}^{t^{\prime}\times t^{\prime}}. The key idea is that for OMP, \mathbf{H}\big{|}_{t^{\prime}-1} and \mathbf{H}\big{|}_{t^{\prime}} differ only by one column, which is the column selected by OMP at step $t^{\prime}$, and therefore $\mathbf{L}^{(t^{\prime})}$ can be easily computed from $\mathbf{L}^{(t^{\prime}-1)}$ as

[TABLE]

where $\mathbf{h}^{(t^{\prime})}\in\mathbb{R}^{N}$ denotes the column of $\mathbf{H}$ selected at iteration $t^{\prime}$ of COMP. The code \mathbf{x}\big{|}_{t^{\prime}} is obtained as the solution to \mathbf{L}^{(t^{\prime})}\big{(}\mathbf{L}^{(t^{\prime})}\big{)}^{\text{T}}\mathbf{x}|_{t^{\prime}}=\mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{y}, which can be solved more efficiently than the linear system \mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{H}\big{|}_{t^{\prime}}\mathbf{x}|_{t^{\prime}}=\mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{y}.

We extend this idea to the convolutional case, noting that Eq. 12 still requires us to construct \mathbf{H}\big{|}_{t^{\prime}} and to perform multiplications that are expensive in terms of memory and computation. We replace the multiplication operation involving \mathbf{H}\big{|}_{t^{\prime}}^{\text{T}} with a cross-correlation operation, as outlined in Algorithm 1. To keep track of the selected filters, we utilize two sets: 1) The set $\mathcal{S}$ of active template indices chosen by COMP and defined in Section III, and 2) the set $\mathcal{I}$ of times when events associated with each of the templates from $\mathcal{S}$ occur. At COMP iteration $t^{\prime}$, $\mathbf{h}_{\mathcal{S}_{i}^{(t^{\prime})}}$ refers to the template selected at $i^{\text{th}}$ iteration, where $i\leq t$, and $\mathcal{I}^{(t^{\prime})}_{i}$ refers to the time of occurrence of the corresponding template.

With $\mathcal{S}$ and $\mathcal{I}$, neither the convolutional matrix \mathbf{H}\big{|}_{t^{\prime}}, nor the zero-padded filters are required, which makes COMP very efficient. We emphasize that the efficiency gain from the modified projection step outweighs that of the selection step, as we demonstrate in Section VI. Specifically, the efficient projection step essentially removes the huge performance gap between the CMP and the COMP.

Algorithm We summarize the implementation of the efficient COMP-INTERP in Algorithm 2. When $K=1$, COMP-INTERP is equivalent to COMP with the original, non-interpolated, dictionary. The INTERPOLATE function refers to the process of obtaining interpolated templates. The CHOLESKY function refers to the efficient projection step. The superscript $(\ast)$ denotes the quantities at convergence. As in OMP, the convergence criterion can either be when the residual falls below a certain threshold or when a certain sparsity level is reached.

IV-E Comparison with the continuous basis approach

The continuous-basis approximation methods, CBP and continuous OMP, and COMP-INTERP all minimize the effect of time quantization errors due to events off the grid. All three methods assume that the continuous-time templates are smooth. In the continuous-basis approximation approach, the derivation of the mapping that relates the local basis and the templates relies on the smoothness assumption. In COMP-INTERP, the bandlimited assumption implies the smoothness of the templates.

The two approaches differ in how the templates are used, as illustrated in Fig. 4. As detailed in Section III, the continuous-basis approach approximates continuous-time shifts of the templates by using the linear combination of a set of functions indexed by points on the discrete-time sampling grid. COMP-INTERP takes an approach that is simpler and more direct: it uses bandlimited interpolation to approximate continuous-time shifts of the templates by multiples of a arbitrarily fine interval $\Delta_{K}$. This leads to an overcomplete, interpolated dictionary, comprising the original templates and their interpolated variants. In fact, the dictionary $\mathbf{\widetilde{H}}$, can be thought of as a discrete approximation to the infinite-dimensional dictionary that encodes all possible continuous-time shifts  [25].

The advantages of our approach are its simplicity and its speed. It is simpler as it requires neither a local basis nor a mapping between said basis and the templates. It is faster due to its greedy nature, and because of the efficient implementations of operations that involve the convolutional dictionary.

V Convolutional Dictionary Update with interpolated dictionary

COMP-INTERP is an algorithm to solve the CSC step of CDL using a interpolated dictionary that can approximate continuous-time shifts. In this section, we develop an algorithm to solve the CDU step using the interpolated dictionary.

The CDU step involves an optimization problem with respect to $\{\mathbf{h}_{c}\}_{c=1}^{C}$. To simplify it, we first re-write $\mathbf{H}\mathbf{X}_{j}$ as

[TABLE]

where we introduce the matrix representation $\mathbf{S}_{j,i}^{c}\in\mathbb{R}^{W\times L}$ of the linear operator that shifts $\mathbf{h}_{c}$ by $n_{j,i}^{c}$ samples. $\mathbf{S}_{j,i}^{c}$ is a zero-padded identity matrix defined as follows

[TABLE]

Eq. 13 is a result of the commutativity of convolution operation. It allows us to re-write the optimization problem in Eq. 4, with respect to $\{\mathbf{h}_{c}\}_{c=1}^{C}$, as follows

[TABLE]

Compared to Eq. 4, Eq. 15 is simpler because $\mathbf{h}_{c}$ appears as a vector, as opposed to a matrix. Using the interpolated dictionary $\widetilde{\mathbf{H}}$ to account for non-integer delays, we can write the objective similarly

[TABLE]

where $\mathbf{F}^{k}$ is the sinc interpolator matrix defined previously. The constraint enforces the fact we obtain $\mathbf{h}_{c^{k}}$ from $\mathbf{h}_{c}$ by interpolation, as discussed in Section IV. As is customary in the CDL literature, we assume the templates have unit norm, a constraint that we enforce by normalizing the solutions of Eq. 15 or Eq. 16.

V-A CDU with basic dictionary

To solve Eq. 15, we can use any least-squares based algorithm for dictionary. We focus on KSVD [11], specifically shift-invariant KSVD [26]. Two key ideas from KSVD are worth emphasizing: 1) the templates $\{\mathbf{h}_{c}\}_{c=1}^{C}$ are updated sequentially, i.e. one at a time, and 2) when updating a given template, only windows of the data where the template occurs need to be considered.

Suppose we want to update template $\mathbf{h}_{\overline{c}}$. Let $\Omega_{\overline{c}}\subset\{1,\cdots,J\}$ denote the set of indices of windows from which COMP selects at least one occurrence of $\mathbf{h}_{\overline{c}}$. For each window $j\in\Omega_{\overline{c}}$, we split the sum from Eq. 15 into two parts, namely one that involves $\mathbf{h}_{\overline{c}}$ and another that involves the remaining templates. The new estimate of $\mathbf{h}_{\overline{c}}$, denoted by $\widehat{\mathbf{h}}_{\overline{c}}$ is given by

[TABLE]

is the residual from approximating $\mathbf{Y}_{j}$ with templates other than $\mathbf{h}_{\overline{c}}$. Eq.  17 is a least-squares problem whose solution $\widehat{\mathbf{h}}_{\overline{c}}$ is given by

[TABLE]

We can interpret Eq. 18 as the weighted average of segments $\mathbf{E}_{j}$, or \big{(}\mathbf{S}_{j,i}^{\overline{c}}\big{)}^{\text{T}}\mathbf{E}_{j}\in\mathbb{R}^{L} in which template $\mathbf{h}_{\overline{c}}$ occurs. The average is normalized by a factor that accounts for occurrences of $\mathbf{h}_{\overline{c}}$ that overlap: the term \big{(}\mathbf{S}_{j,i}^{\overline{c}}\big{)}^{\text{T}}\mathbf{S}_{j,m}^{\overline{c}}\in\mathbb{R}^{L\times L} is a matrix that is non-zero only if the offset $|n_{j,i}^{\overline{c}}-n_{j,m}^{\overline{c}}|$ between occurrences of $\mathbf{h}_{\overline{c}}$ is less than the template length $L$.

V-B CDU with interpolated dictionary (CDU-INTERP)

To solve Eq. 16, we follow an approach similar to that used to solve Eq. 15. The constraint from Eq. 16 implies that

[TABLE]

where $\widetilde{\mathbf{S}}_{j,i}^{c^{k}}:=\mathbf{S}_{j,i}^{c^{k}}\mathbf{F}^{k}$ is a modified shift operator (with $\mathbf{F}^{k}$, instead of $\mathbf{I}_{L\times L}$, placed in the block starting from row $n_{j,i}^{c^{k}}$). Eq. 19 allows us to rewrite Eq. 16 as the unconstrained optimization problem

[TABLE]

Suppose we want to update template $\mathbf{h}_{\overline{c}}$. Let $\widetilde{\Omega}_{\overline{c}}=\cup_{k=1}^{K}\Omega_{\overline{c}^{k}}$ be the set of indices of windows from which COMP selects at least one occurrence of a template from the set $\{\mathbf{h}_{\overline{c}^{k}}\}_{k=1}^{K}$. Re-arranging Eq. 20 yields the estimate $\widehat{\mathbf{h}}_{\overline{c}}$ of $\mathbf{h}_{\overline{c}}$

[TABLE]

the solution of which is given by

[TABLE]

Similar to Eq. 18, we can interpret Eq. 22 as the weighted average of the interpolated segments $\mathbf{E}_{j}$ in which template any template from the set $\{\mathbf{h}_{\overline{c}^{k}}\}_{k=1}^{K}$ occurs.

V-C CDL algorithm

We summarize the alternating-minimization procedure for CDL in Algorithm 3. COMP-INTERP refers to the Algorithm 2 and CDU refers to the CDU step. The $\ast$ notation from $\mathcal{S}^{j,(\ast)}$ and $\mathcal{I}^{j,(\ast)}$ refer to the index sets for window at COMP convergence.

COMP-INTERP is parallelizable across $J$ windows, and therefore amenable to implementation on GPU and parallel computer architectures. The CDU step is not parallelizable because it needs to aggregate the occurrences of the templates across all $J$ windows.

VI Experiments

We apply the proposed CDL framework to simulated data and the real electrophysiological data from the brain. We use two criteria to evaluate performance: 1) the accuracy and speed of the CSC step, and 2) the accuracy of the CDU step.

For the CSC step, we compare CBP to the following convolutional greedy methods: COMP, CMP, COMP-slow, and COMP-INTERP. CMP does not have a projection step. The projection step from COMP-slow inverts \mathbf{H}\big{|}_{t^{\prime}}^{\text{T}}\mathbf{H}\big{|}_{t^{\prime}} directly without Cholesky factorization. Table II summarizes the similarities and differences of the greedy methods.

VI-A Simulated Dataset

We simulated a signal according to the continuous-time generative model of Eq. 1. We used two $10$-ms-long gamma-tone templates [6] defined for $t$ in the interval from $-5$ to $5$ ms

[TABLE]

We assumed the same number of occurrences $N_{1}=N_{2}$ of the templates. We chose the times when events occur uniformly at random, i.e. $\tau_{i}^{c}\sim\text{Uniform}[0,T]$ for $c=1,2$ and $i=1,\cdots,N_{1}$. We chose the amplitude of each event uniformly at random, i.e. $x_{i}^{c}\sim\text{Uniform}[1,2]$. As explained in subsequent sections, we used a range of values for the variables $T$ and $N_{1}$. We used a sampling rate $f_{s}=10^{4}$ Hz and obtained the discrete-time signal $\mathbf{y}$ by sampling $y(t)$ at every $\Delta=10^{-4}$ seconds. We added white Gaussian noise $\boldsymbol{\varepsilon}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I})$ in discrete time, where $\sigma$ was set according to a desired Signal-to-noise ratio (SNR). We defined the SNR as follows

[TABLE]

where the variance estimate $\widehat{\sigma}^{2}$ is $\widehat{\sigma}^{2}=\frac{\lVert\varepsilon\rVert_{2}^{2}}{N}$, with $N=\lfloor\frac{T}{\Delta}\rfloor$. We obtained $\mathbf{h}_{1},\mathbf{h}_{2}\in\mathbb{R}^{100}$ by acquiring 100 samples from $h_{1}(t)$ and $h_{2}(t)$ and normalizing the resulting vectors to have unit length: $\lVert\mathbf{h}_{1}\rVert_{2}=\lVert\mathbf{h}_{2}\rVert_{2}=1$.

VI-B Results from simulations: CSC step

We set the sparsity level for greedy methods to be the number of events, and fine-tuned the regularization parameter $\lambda$ for CBP (with polar basis) to match the same sparsity level. We use the true templates $\mathbf{h}_{1},\mathbf{h}_{2}$ as dictionary elements.

Sparse-coding computation time We computed the duration of the CSC step, using CBP and the above-mentioned greedy methods, as a function of data length $T$ and total number of events $N_{1}+N_{2}$. Specifically, we ran two sets of experiments: 1) $T\in[0.5,1,2,3,4,5]\,\text{s}$ with fixed $N_{1}+N_{2}=30$ and 2) $N_{1}+N_{2}\in[10,20,30,40,50]$ with fixed $T=3$ s. For each experiment, we report durations averaged over 50 independent repeats (trials). Tables III and IV show the duration of the CSC step for the various methods, respectively as a function of data length and number of occurrences of events. We draw three conclusions from the results

COMP is computationally much faster than CBP, with two possible explanations. In terms of implementation, CBP constructs the full convolutional dictionary for $\ell_{1}$ regularization, whereas COMP does not. Moreover, it is well-known that greedy methods are faster than ones based on basis-pursuit like methods [7], [9]. 2. 2.

COMP is as fast as CMP. This is true even for large $T$ or $N_{1}+N_{2}$, which involve a computationally-demanding projection step. This, along with a comparison of COMP to COMP-slow highlights the importance of making the projection step efficient. COMP reduces the computation time of COMP-slow by 48% $\sim$ 85% on average. 3. 3.

The difference between COMP and COMP-INTERP stems from the fact that the selection step of the former requires $C$ cross-correlation operations, while that of the latter requires $KC$ such operations. The computation time of the projection step is the same for both.

Sparse coding accuracy We computed the average hit error [21] for COMP and CBP, and that for COMP-INTERP as a function of the discretization $\Delta_{K}$. The average hit error measures how far, in terms of absolute displacement, the recovered sparse codes are relative to the true sparse codes. We simulated 50 trials of data with $T=1$ and $N_{1}=N_{2}=10$ and computed the median average hit error across trials for each method. The average hit error of both COMP and CBP is independent of $\Delta_{K}$: COMP operates at the sampling resolution $\Delta$, while CBP operates in continuous time.

Fig. 5 shows that the average hit error for COMP-INTERP is lower than that for COMP. The finer the discretization, the greater the precision with which COMP-INTERP identifies the sparse codes, resulting in the lower average hit error. The reduction in average hit error is marginal below a certain threshold ($\Delta_{10}=\frac{1}{10}\Delta$). The finer the discretization, the more correlated the interpolated templates and, in turn, the columns of the dictionary are. Therefore, as discretization becomes smaller, the returns, in terms of average hit error, from expanding the dictionary diminish. With no discretization, COMP and COMP-INTERP achieve the same average hit errors since the original and the interpolated dictionary coincide. Fig. 5 also shows that there is a small gap, between CBP and COMP-INTERP. We observed in simulations that the difference comes from the fact that COMP-INTERP is not as accurate as CBP in the presence of significant template overlaps. Due to its greedy nature, OMP has worse performance in resolving overlapping templates than approaches based on basis pursuit [13], [27]. This can be mitigated by heuristics [28], which we leave it as a future work.

VI-C Results from simulations: CDU step

Dictionary learning accuracy To assess the accuracy of the CDU algorithms, we use the following error distance metric  [29] between two templates $\widehat{\mathbf{h}}_{c}$ and $\widetilde{\mathbf{h}}_{c}$

[TABLE]

The lower the metric, the closer $\widehat{\mathbf{h}}_{c}$ and $\widetilde{\mathbf{h}}_{c}$ are. If $\widehat{\mathbf{h}}_{c}=\mathbf{h}_{c}^{(t)}$ and $\widetilde{\mathbf{h}}_{c}=\mathbf{h}_{c}$, the metric measures how close the learned template at iteration $t$ is to the true template. If $\widehat{\mathbf{h}}_{c}=\mathbf{h}_{c}^{(0)}$ and $\widetilde{\mathbf{h}}_{c}=\mathbf{h}_{c}^{(\ast)}$, the metric measures how much the template at convergence has changed from the initial template.

We compared the accuracy of the CDU step with the interpolated and the original dictionary. We simulated 25 trials of data, with $T=5$ s and $N_{1}=N_{2}=200$, for several levels of SNR ranging from 8 to 26 dB. We obtained the initial templates $\mathbf{h}_{1}^{(0)}$ and $\mathbf{h}_{2}^{(0)}$ by perturbing the original templates $\mathbf{h}_{1}$ and $\mathbf{h}_{2}$ with Gaussian additive noise, such that $\text{err}(\mathbf{h}_{c}^{(0)},\mathbf{h}_{c})\geq 0.5$ for $c=1,2$. We performed 15 iterations of the alternating-minimization algorithm for CDL to obtain $\mathbf{h}_{1}^{(\ast)}$ and $\mathbf{h}_{2}^{(\ast)}$. We compared the following two methods: 1) COMP & CDU with discretization $\Delta$ (CDL with $\Delta$), which corresponds to conventional CDU methods that use $\mathbf{H}$ instead of $\widetilde{\mathbf{H}}$, and 2) COMP-INTERP & CDU-INTERP with discretization $\Delta_{10}=\frac{1}{10}\Delta$ (CDL with $\Delta_{10}$).

Fig. 6 shows the result of applying the two CDL algorithms to the simulated data. The figure shows that CDL with $\Delta_{10}$ is able to learn templates more accurately than CDL with $\Delta$. This is supported both visual inspection of the learned templates (Fig. 6(a) and (b)), and the trajectory of $\text{err}(\mathbf{h}_{1}^{(t)},\mathbf{h}_{1})$ as a function of alternating minimization iteration $t$, for one of the simulated trials. Although not shown here, we observed that CDL with $\Delta_{10}$ learned $\mathbf{h}_{2}^{(\ast)}$ more accurately than CDL with $\Delta$. CDL with $\Delta$, which does not account account for the non-integer shift, converges to an average of the non-integer shifts of each of the true templates. Fig. 7 shows that both CDL approaches are robust across a range of SNRs.

VI-D Real dataset: Application to spike sorting

We applied our framework to spike sorting. Given a recording of extracellular voltage, the goal of spike sorting is to learn the action potentials (templates) from neurons (sources) near the electrode, and the times when the action potentials occur [10]. We used a dataset that consists of an extracellular recording from the rat hippocampus, along with a simultaneous intracellular recording [30] from one neuron. The intracellular recording provides the ground truth data, as it provides the unequivocal occurrence of an action potential in a single neuron, and thus enables us to evaluate and compare the accuracy of the CDL frameworks. For this dataset, 620 events occurred from the neuron that was recorded intracellularly. The sampling rate of the extracellular data, which comprise 4 channels, is $f_{s}=10^{4}$ Hz. We used $T=150$ seconds of data from channel 1 and preprocessed them following the procedures that are standard in the spike-sorting literature [31]. Specifically, we applied to the data a high-pass filter with cut-off frequency 400 Hz, and whitened it. In addition, we identified peaks from the data that crossed a pre-defined threshold [32] and extracted a segment of length $81$ samples centered around each peak. The resulting collection $\mathbf{Y}\in\mathbb{R}^{5000\times 81}$ of $5,000$ segments is the input to our analyses of the real data.

Method setup We assumed $C=3$, namely that the extracellular recording can detect activity from 3 neurons. We used templates $\{\mathbf{h}_{c}\}_{c=1}^{3}\in\mathbb{R}^{41}$, each of length 4 ms. We applied CDL, with COMP and COMP-INTERP in the CSC step, to the extracellular data. We used the following procedure to initialize the templates to \big{\{}\mathbf{h}_{c}^{(0)}\big{\}}_{c=1}^{3}. Following segment extraction, we first performed PCA on $\mathbf{Y}$ for dimensionality reduction, and then K-means clustering with three clusters in the lower dimension space. We used the centroids of the clusters to obtain \big{\{}\mathbf{h}_{c}^{(0)}\big{\}}_{c=1}^{3}. We used an estimate of the variance of the background noise in the recording as the termination criterion for the CSC step of COMP and COMP-INTERP. We computed this estimate by extracting data from a segment that remained below a pre-defined threshold for more than 500 ms. We ran 15 iterations of the CDL algorithm to obtain $\{\mathbf{h}_{c}^{(\ast)}\}_{c=1}^{3}$.

We compared CBP with $\{\mathbf{h}_{c}^{(0)}\}_{c=1}^{3}$ to the CSC step of COMP and COMP-INTERP using the filters, $\{\mathbf{h}_{c}^{(\ast)}\}_{c=1}^{3}$, learned by the respective CDU algorithms. For a given true spike event from the intracellular data, we associate an event identified using the extracellular data as a true positive if the event is within 30 samples (3 ms) of a true event from the intracellular data. Among the templates learned using the extracellular data, we associate the template with the highest true positive count with the neuron from the intracellular data and refer to it as $\mathbf{h}_{1}$. As we do not have access to the true $\mathbf{h}_{1}$, we treat $\mathbf{h}_{1}^{(\ast)}$ as the best estimate of the true spike template.

VI-E Results from the real dataset

Detection error curve We used two statistics to evaluate the ability of the methods to perform spike sorting on the extracellular data. Following CSC, each method must set a threshold to identify the times when action potentials (spikes) occur. A true miss is a true spike from the intracellular data within 3 ms of which no threshold-crossing event occurs in the extracellular data. The true miss rate is the ratio of the number of true misses to the number of intracellular events. A false alarm is a threshold-crossing event from the extracellular data that is not a true spike, i.e. within 3 ms of which no events occur in the intracellular data. The false alarm rate is the ratio of the number of false alarms to the number of threshold crossing events. Varying the threshold leads to a trade-off between true misses and false alarms. A threshold with high amplitude typically leads to a low number of false alarms and a large number of true misses, and vice versa for a low-amplitude threshold.

Fig. 8 shows the result of sorting spikes associated with $\mathbf{h}_{1}$ using CBP, COMP and COMP-INTERP. The figure shows that the greedy approaches rival CBP, and are better in the low true miss rate regime. The true miss rate for CBP does not decrease below 5%, even with thresholds of low amplitude. This indicates that CBP is not able to identify a subset of true events that COMP and COMP-INTERP identify correctly.

Difference between $\ell_{0}$ and $\ell_{1}$ The discrepancy in the true miss rate of CBP and COMP/COMP-INTERP motivated us to further examine segments for which the number of errors from CBP and the greedy methods differ. Fig. 9 shows examples of such segments. Fig. 9(a) shows that CBP fails to capture the true spike event (red dot), resulting in a true miss event, whereas COMP-INTERP (Fig. 9(b)) uses $\mathbf{h}_{1}$ to correctly identify the event. The failure of CBP and the success of COMP-INTERP point to an important difference between the $\ell_{1}$ and $\ell_{0}$-based methods for the CSC step. To minimize its objective function, $\ell_{1}$-based CBP must strike a balance between the reconstruction error, which it can reduce by using additional templates, and the $\ell_{1}$ penalty, which can be reduced by using fewer templates or ones with lower amplitude. In Fig. 9(a), the choice of regularization parameter $\lambda$ is such that CBP chooses to use one template with large amplitude, thereby missing the true event. Although a smaller regularization parameter $\lambda$ can be used to reduce the effect of the $\ell_{1}$ penalty, thereby identifying the true event and lowering the true miss rate, this would result in spurious events that would increase the number of false alarms. This discussion points to a limitation of $\ell_{1}$-based methods, namely the need to tune $\lambda$ carefully. COMP-INTERP (and COMP), on the other hand, can select as many events as needed to make the reconstruction error below the estimate of the variance of the background noise. In Fig. 9(b), COMP-INTERP first selects $\mathbf{h}_{3}$, and then $\mathbf{h}_{1}$, which corresponds to the true event.

Example of a non-integer shift Fig. 9(c) is an example of a segment where COMP raises a false alarm event, but COMP-INTERP does not (Fig. 9(d)). COMP is forced to use the true spike template, $\mathbf{h}_{1}$, whereas COMP-INTERP uses $\mathbf{h}_{3}$ to select the secondary peak. This likely happens because the spike event does not occur on the discrete-time sampling grid, and further highlights the benefits of using the interpolated dictionary. That being said, we observe from Fig. 8 that the two versions of COMP have similar performance, with COMP-INTERP slightly outperforming COMP in the low true miss regime.

Learned templates Fig. 10 shows the templates that were learned by CDL using the extracellular data, and COMP-INTERP in the CSC step. The fact that the shapes of the learned templates are not significantly different from those of the initial templates suggests that we initialized the templates well. To determine how different the learned templates are from the initial ones, we computed $\text{err}(\mathbf{h}_{c}^{(0)},\mathbf{h}_{c}^{(\ast)})$ for $c=1,2,3$, and then took the maximum of the three values. We found that the maximum equaled 0.32, which indicates that, although not obvious visually, the CDL algorithm did learn new templates.

VII Discussion

We have introduced novel Convolutional Sparse Coding (CSC) and Convolutional Dictionary Update (CDU) algorithms for a class of continuous-time signals that consist of time-shifted copies from multiple sources, each with its own template. The algorithms operate in discrete-time and construct an expanded, overcomplete, dictionary that accounts for the fact that the templates from the continuous-time model do not necessarily occur on the sampling grid of the discrete-time model. Specifically, the elements of the expanded dictionary consist of the original templates, along with their non-integer shifts, obtained by bandlimited interpolation.

The use of an expanded dictionary increases the computational demands of the algorithms, the CSC step in particular. To mitigate this, we chose to focus on greedy pursuit methods and proposed an efficient implementation of convolutional OMP (COMP) and COMP with interpolated dictionary (COMP-INTERP). The efficient COMP algorithm exploits the locality of the templates and is much faster than CBP, which is the state-of-the-art algorithm for CSC in continuous-time. Compared to convolutional MP, the overhead from the efficient COMP algorithm is negligible. We also demonstrated a novel algorithm to perform dictionary update with an expanded convolutional dictionary that accounts for non-integer delays. We showed empirically that more accurate templates were learned with the interpolated dictionary and that the results are robust across a wide range of SNRs. We conclude that our approach is a simple yet efficient paradigm for convolutional dictionary learning, that faithfully accounts for the continuous-time nature of the signal of interest.

We believe that the current work can be extended in several promising directions. First, we plan to extend this framework to the multivariate case. This will be useful for performing dictionary learning using recordings from sensor arrays, a setting in which the spatial information becomes as important as the temporal information of the signal. Then, we would like to mitigate COMP sparse coding errors that are due to template overlaps. Several interesting heuristics, such as bounded local optimization [28], have been proposed for this. These can be incorporated into our framework in efficient way. Next, since the CSC step is embarrassingly parallelizable across windows, we plan to implement the algorithm on GPU.