<figure><img src="image_1.png" /></figure>

HAL
open science

# Convolutional Sparse Coding for Time Series Via a l0 Penalty: An Efficient Algorithm With Statistical Guarantees

Charles Truong, Thomas Moreau

► To cite this version:

Charles Truong, Thomas Moreau. Convolutional Sparse Coding for Time Series Via a l0 Penalty: An Efficient Algorithm With Statistical Guarantees. Statistical Analysis and Data Mining, 2024, 17 (6), pp.e70000. 10.1002/sam.70000. hal-04856979

HAL Id: hal-04856979

https://hal.science/hal-04856979v1

Submitted on 27 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RAPID PUBLICATION

OPEN ACCESS

# Convolutional Sparse Coding for Time Series Via a $l_0$ Penalty: An Efficient Algorithm With Statistical Guarantees

Charles Truong¹

Thomas Moreau²

¹Centre Borelli, Université Paris-Saclay, CNRS, ENS Paris-Saclay, Gif-sur-Yvette, France | ²Université Paris-Saclay, Inria, CEA, Palaiseau, France

**Correspondence:** Charles Truong (charles.truong@ens-paris-saclay.fr)

**Received:** 23 May 2024 |

**Revised:** 24 October 2024 |

**Accepted:** 28 November 2024

**Funding:** This work was supported by PhLAMES chair (ENS Paris-Saclay).

**Keywords:** compressive statistical learning | convolutional sparse coding | L0 penalty | signal processing

## ABSTRACT

Identifying characteristic patterns in time series, such as heartbeats or brain responses to a stimulus, is critical to understanding the physical or physiological phenomena monitored with sensors. Convolutional sparse coding (CSC) methods, which aim to approximate signals by a sparse combination of short signal templates (also called atoms), are well-suited for this task. However, enforcing sparsity leads to non-convex and untractable optimization problems. This article proposes finding the optimal solution to the original and non-convex CSC problem when the atoms do not overlap. Specifically, we show that the reconstruction error satisfies a simple recursive relationship in this setting, which leads to an efficient detection algorithm. We prove that our method correctly estimates the number of patterns and their localization, up to a detection margin that depends on a certain measure of the signal-to-noise ratio. In a thorough empirical study, with simulated and real-world physiological data sets, our method is shown to be more accurate than existing algorithms at detecting the patterns' onsets.

## 1 | Introduction

Analyzing continuous temporal data often involves identifying characteristic temporal patterns, which reveal underlying physical phenomena [1–3]. For physiological signals such as ECG [4] or gait recordings [5], detecting anomalous patterns such as irregular heartbeats can help diagnose various pathologies. In neuroscience, characteristic brain responses to stimulus are used to link brain areas and cognitive functions [6]. When prototypical patterns are known, they can be used to detect the occurrence of specific events. This task is often challenging due to varying pattern amplitude or length, pattern similarity, and noise, which produce spurious detections. Extra constraints are considered to improve the detection rate and filter out spurious events to

produce satisfactory results. A typical constraint for physiological signals is to impose that the detected patterns do not overlap. While this extra assumption helps reduce spurious detections, this also makes the algorithmic task harder.

To detect characteristic patterns, various algorithms have been developed. A classical method is the *matched filtering*, which consists in constructing a bank of templates and comparing them individually with the signal to detect related events. It is, for instance, used to detect synaptic events [3], gravitational waves [2], or seismic events [7]. This simple and efficient method is based on cross-correlation and thresholding and performs well. Yet, imposing additional constraints requires tedious post-processing that can yield suboptimal results [8].

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). Statistical Analysis and Data Mining published by Wiley Periodicals LLC.
Other approaches to detect patterns are based on convolutional sparse coding (CSC) [1, 9, 10]. These methods aim to find the best signal reconstruction using only a few atoms or templates. In practice, this reconstruction is computed by solving a convex optimization problem regularized with the $\ell_1$-norm, which is easily adapted to account for additional constraints. This method is solved efficiently due to its convexity, using proximal gradient descent [11] or coordinate descent algorithms [12, 13]. However, the relaxation can create spurious detections that must be post-processed with peak detection algorithms. To improve performances, $\ell_0$-based extensions of CSC have been proposed [14, 15]. These methods allow for sharper detection but make the optimization problem harder due to the non-convex and combinatorial nature of the $\ell_0$-norm [16]. Efficient greedy heuristics have been proposed and give good results in practice. However, only a few works [17, 18] have studied their statistical properties.

In this work, we focus on time series that contains non-overlapping patterns. This simplifying assumption is justified in numerous contexts. In many physiological signals, the observed phenomena happen one after the other without overlap: heartbeats in electrocardiograms (ECG) [4, 19], footsteps in human gait analysis [20, 21], breathing cycles in photoplethysmogram data (PPG) [22] or whole-body plethysmogram data [23, 24], spikes in Local Field Potential (LFP) time series [25], etc. This is also verified in certain industrial settings [26, 27]. In the data mining community, the task of motif discovery assumes that patterns do not overlap Germain et al. [28], Linardi et al. [29], Gao and Lin [30, 31], Ceccarello and Gamper [32], Alae et al. [33].

## 1.1 | Contributions

In this paper, we develop an efficient algorithm for non-overlapping pattern detection in time series based on an $\ell_0$-based CSC approach. Our algorithm is the first that is guaranteed to find the optimal solution to the CSC problem under a no-overlap constraint. Our method is easy to implement and computationally efficient with linear complexity. In a theoretical analysis, we prove that our method correctly estimates the number of patterns and their localization up to a detection margin that depends on a certain measure of the signal-to-noise ratio. The benefits of computing the optimal activations, as opposed to the classical convex relaxation and the greedy approaches, are shown in several biomedical applications.

# 2 | Problem Statement

## 2.1 | Notations

Let $\mathbf{d}_j \in \mathbb{R}^L$, $j = 0, \dots, p-1$, be $p$ atoms (or filters) of length $L$; atoms are assumed to have unit norm, that is, $\|d_j\|^2 = 1$. The matrix $\mathbf{D} := [\mathbf{d}_0 \mathbf{d}_2 \dots \mathbf{d}_{p-1}] \in \mathbb{R}^{L \times p}$ whose columns are the $\mathbf{d}_j$ is called the local dictionary. For two time series $\mathbf{x} \in \mathbb{R}^{T_x}$ and $\mathbf{y} \in \mathbb{R}^{T_y}$, the convolution of $\mathbf{x}$ and $\mathbf{y}$ is denoted by $(\mathbf{x} * \mathbf{y}) \in \mathbb{R}^{T_y+T_x-1}$. For any matrix $\mathbf{M}$, $\mathbf{M}[t]$ will denote the block of $p$ columns of $\mathbf{M}$ from $tp$ to $(t+1)p-1$. For a vector $\mathbf{Z}$, the $t$th block $\mathbf{Z}[t]$ contains the $p$ coefficients from $tp$ to $(t+1)p-1$. Define

$\|\mathbf{Z}\|_{0,1} := \sum_j, 1(\mathbf{Z}[t] \neq 0)$. Denote by $\mathcal{T}(\mathbf{Z})$, or simply $\mathcal{T}$, the temporal support of $\mathbf{Z}: \mathcal{T}(\mathbf{Z}) := \{t \text{ s.t. } \mathbf{Z}[t] \neq 0\}$. Note that $|\mathcal{T}| = \|\mathbf{Z}\|_{0,1}$. Finally, let $(a..b)$ denote the integer interval $\{a, a+1, \dots, b-1\}$ for two integers $a < b$; $\mathbf{y}_{a..b}$ is the sub-signal $\{y_a, y_{a+1}, \dots, y_{b-1}\}$.

## 2.2 | Model and Task

Assume that we observe a real-valued signal $\mathbf{y} = (y_0, y_2, \dots, y_{T-1})$ with $T$ samples, generated according to the linear convolutional model [34]:

$$ \mathbf{y} = \sum_{j=0}^{p-1} \mathbf{d}_j * \mathbf{z}_j^* + \sigma\varepsilon \quad (1) $$

where the $\mathbf{z}_j^*$ are $p$ unknown sparse activation signals of length $\tilde{T} := T - L + 1$, $*$ is the convolution operator, $\varepsilon \in \mathbb{R}^T$ is a noise signal and $\sigma > 0$ is the noise level.

Model (1) can be rewritten in linear form. Indeed, for a given set of activation signals $\{\mathbf{z}_0, \mathbf{z}_1, \dots, \mathbf{z}_{p-1}\} \subset \mathbb{R}^{\tilde{T}}$, define their vectorized form $\mathbf{Z} \in \mathbb{R}^{\tilde{T}^p}$ obtained by concatenating the activations at each timestep, that is, its $t$th block is equal to $\mathbf{Z}[t] := [z_{0,t}, z_{1,t}, \dots, z_{p-1,t}]^T \in \mathbb{R}^p$. Since $\mathbf{Z} \mapsto \sum_{j=0}^{p-1} \mathbf{d}_j * \mathbf{z}_j$ is a linear map from $\mathbb{R}^{\tilde{T}^p}$ to $\mathbb{R}^T$, there exists a matrix $\mathbf{\Phi} \in \mathbb{R}^{(\tilde{T}^p) \times T}$ such that $\mathbf{\Phi}\mathbf{Z} = \sum_{j=0}^{p-1} \mathbf{d}_j * \mathbf{z}_j$. The matrix $\mathbf{\Phi}$ is called the global dictionary (while $\mathbf{D}$ is the local dictionary). A few algebraic manipulations show that $\mathbf{\Phi}$ is a block diagonal matrix given by concatenating all shifted versions of the local dictionary $\mathbf{D}$, as in Figure 1 (see Section A in the Supporting Information: Appendix for a formal definition). Model (1) can be rewritten as follows:

$$ \mathbf{y} = \mathbf{\Phi}\mathbf{Z}^* + \sigma\varepsilon \quad (2) $$

where $\mathbf{Z}^* \in \mathbb{R}^{\tilde{T}^p}$ is such that its $t$th block is equal to $\mathbf{Z}^*[t] := [z_{0,t}^*, z_{1,t}^*, \dots, z_{p-1,t}^*]^T \in \mathbb{R}^p$.

The objective of CSC is to estimate the sparse activation signals $\mathbf{z}_j^*$ from the noisy observations $\mathbf{y}$. In this work, we further assume that the true patterns do not overlap. Formally, introduce $\mathcal{N}\mathcal{O}$ denote the set of all possible temporal supports that satisfy the no-overlap constraint, meaning that $\mathcal{T} = \{t_1, t_2, \dots\} \in \mathcal{N}\mathcal{O}$ iff $0 \le t_1 < t_2 < \dots \le \tilde{T}$ and $t_{k+1} - t_k > L$. Let $\mathcal{Z} \subset \mathbb{R}^{\tilde{T}^p}$ be the set of activation matrices that satisfy the no-overlap constraint, that is, $\mathcal{Z} := \{\mathbf{Z} \in \mathbb{R}^{\tilde{T}^p \text{s.t. } \mathcal{T}(\mathbf{Z}) \in \mathcal{N}\mathcal{O}}\}$. In the following, we, therefore, assume that $\mathbf{Z}^* \in \mathcal{Z}$.

## 2.3 | Related Statistical Guarantees

The CSC problem can be seen as a block sparse coding task [13]. In this setting, several works have proved conditions for the exact recovery of the support $\mathcal{T}^*$ of $\mathbf{Z}^*$ for greedy algorithms, in the noiseless case [35] or in the noisy case [36, 37]. Those conditions are based on (i) the mutual coherence of the local dictionary $\mathbf{D}$, denoted $\mu_0$ and defined by $\mu_0 := \max_{i,j,i \neq j} |\langle \mathbf{d}_i | \mathbf{d}_j \rangle|$, and (ii) the block coherence of the global dictionary $\mathbf{\Phi}$, which is the largest spectral norm of the $\mathbf{\Phi}[i]^T \mathbf{\Phi}[j]$, for $i \neq j$. For the assumptions of those works to hold, both coherences must be bounded away from 1. This is not satisfied in the CSC context—except in limited
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Schematic view of the local dictionary **D**, the global dictionary **Φ**, its *t*th block **Φ**[*t*] and the shifted dictionaries **D**(*s*) and **D**(*s*−*) for a positive *s* (introduced in Section 3.6).</figcaption></figure>

situations, for example, *p* = 1—because of the specific shape of the convolutional global dictionary. Indeed, consecutive blocks **Φ**[*t*] and **Φ**[*t* + 1] are almost identical, therefore the block coherence is close to 1.

To cope with this issue, the authors of Papyan et al. [18] have introduced adapted coherence measures: the shifted mutual coherence which is equal to ||**Φ**[*t*]†**Φ**[*t* + *s*]||∞ for a shift *s* and the stripe coherence, which is more complex and depends on both the sparsity pattern **T**\* and the shifted mutual coherence. For a greedy algorithm, they provided exact support recovery conditions in the noiseless setting, extended to the noisy setting in Fu et al. [17]. The authors of Papyan et al. [18] argue that such conditions are easily satisfied when the mutual coherence μ₀ is larger than all shifted coherence values. While this might be true in some contexts—for instance, when the dictionary’s entries are i.i.d. Gaussian—, we empirically found that it is not verified when the atoms are smooth signals representing physiological patterns. (See, for instance, the illustrative example in Section 3.6.)

**3 | Method**

CSC approaches based on the ℓ₀-norm aim to approximate the signal **y** by a sparse combination of the atoms **d**ᵢ convoluted with the activations **z**ᵢ, that is, solve

$$ \hat{\mathbf{Z}} := \arg \min_{\mathbf{Z} \in \mathcal{Z}} \left[ \frac{1}{2} \|\mathbf{y} - \mathbf{\Phi}\mathbf{Z}\|^2 + \lambda \|\mathbf{Z}\|_{0,1} \right] \quad (3) $$

where λ > 0 controls the trade-off between data fidelity and the sparsity of the solution. In the following, we derive our method for univariate signals to keep simple notations. The adaptation required to apply it to multivariate signals is detailed in the Supporting Information.

A key observation for our method is that the space of all non-overlapping activation matrices can be partitioned as follows:

$$ \{ \mathbf{Z} \in \mathbb{R}^{\hat{T}_p} \text{ s.t. } \mathcal{T}(\mathbf{Z}) \in \mathcal{N}\mathcal{O} \} = \bigcup_{\mathcal{T} \in \mathcal{N}\mathcal{O}} \{ \mathbf{Z} \in \mathbb{R}^{\hat{T}_p} \text{ s.t. } \mathcal{T}(\mathbf{Z}) = \mathcal{T} \}. \quad (4) $$

The optimization problem (3) can thus be divided into two steps: searching for the best localization of the patterns **T** and, given this support, selecting the activation values that best reconstruct the signal. This simple observation allows us to replace the operator min<sub>Z∈Z</sub> in (3) by the two minimizations min<sub>**T**∈**N**<sup>0</sup></sub> min<sub>Z, **τ**(Z)=**T**</sub>. Formally, we propose to solve the equivalent optimization problem:

$$ \min_{\mathcal{T} \in \mathcal{N}\mathcal{O}} \left[ \mathbf{z}, \min_{\mathcal{T}(\mathbf{Z})=\mathcal{T}} \|\mathbf{y} - \mathbf{\Phi}\mathbf{Z}\|^2 + \lambda |\mathcal{T}| \right]. \quad (5) $$

The minimization is first performed on the activation amplitudes for a fixed temporal support, then over all possible temporal supports. The following shows that the objective function in (5) can be written as a sum of simple independent problems. This problem structure makes it easy to solve using dynamic programming, and we derive an efficient algorithm able to minimize (5) over all possible supports and attain the global optima of the original problem.

**3.1 | Fixed Activation Support**

When we consider a fixed temporal support **T** ∈ **N**<sup>0</sup> such that **T** = {t₁, t₂, ..., tₖ} with 0 ≤ t₁ ≤ ... ≤ tₖ ≤ ̂T and an activation matrix **Z** with support **T**(Z) = **T**, we can rewrite the reconstruction error:

$$ \|\mathbf{y} - \mathbf{\Phi}\mathbf{Z}\|^2 = \|\mathbf{y}_{0..t_1}\|^2 + \|\mathbf{y}_{t_1..t_1+L} - \mathbf{DZ}[t_1]\|^2 \\ + \|\mathbf{y}_{t_1+L..t_2}\|^2 + \|\mathbf{y}_{t_2..t_2+L} - \mathbf{DZ}[t_2]\|^2 \\ + \dots + \|\mathbf{y}_{t_K..t_K+L} - \mathbf{DZ}[t_K]\|^2 + \|\mathbf{y}_{t_K+L..T}\|^2. $$

Indeed, since there is no overlap between activated atoms, the convoluted signal **Φ**Z = Σᵢ **d**ᵢ * **z**ᵢ is zero outside the intervals (tₖ, tₖ + L). Simple algebraic manipulations and minimizing over the activations **Z** yield the following lemma.

**Lemma 1.** For a signal **y**, and a fixed temporal support **T** = {t₁, ..., tₖ} ∈ **N**<sup>0</sup>, the optimal reconstruction satisfies the following identity

$$ \min_{\mathbf{z}, \mathcal{T}(\mathbf{Z})=\mathcal{T}} \|\mathbf{y} - \mathbf{\Phi}\mathbf{Z}\|^2 + \lambda |\mathcal{T}| = \|\mathbf{y}\|^2 + \sum_{k=1}^{K} c(y_{t_{k}, t_{k+1}}) \quad (6) $$
where by convention, $t_{K+1} := T$ is a dummy index, and for $(s, t)$ such that $0 \le s < t \le T$, $c(\mathbf{y}_{s,t})$ is the cost function of sub-signal $\mathbf{y}_{s,t}$ defined by

$$c(\mathbf{y}_{s,t}) := \lambda + \min_{\mathbf{u} \in \mathbb{R}^p} [\|\mathbf{Du}\|^2 - 2\mathbf{y}_{s,s+L}^{\top}\mathbf{Du}] \quad (7)$$

if $t-s \ge L$ and $\lambda$ otherwise.

Lemma 1 shows that the reconstruction error (6) can be written as a sum of independent segment costs. Note that the quantity $\|\mathbf{y}\|^2$ in (6) is a constant with respect to $\mathbf{Z}$ and $\mathcal{T}$, and will be omitted in subsequent analysis. This sum structure allows us to devise an efficient sequential algorithm to explore the space of all (non-overlapping) activations and thus solve the CSC problem exactly. Our method is described in the following section and a closed-form expression for the cost function $c(\cdot)$ is given in Section 3.3.

## 3.2 | Finding the Optimal Support

This section describes how to solve the outer minimization in (5), that is, minimization over all possible non-overlapping temporal supports. To that end, define $V_t$ the optimal reconstruction cost for the signal $\mathbf{y}_{0,t}$ for $t \le T$:

$$V_t := \min_{K, 0 \le t_1 < \dots < t_K \le t} \sum_{k=1}^{K} c(\mathbf{y}_{t_k, \cdot, t_{k+1}}).$$

By convention, $V_0 := 0$. The optimal reconstruction error for the complete signal is $V_T$ and can be computed recursively, thanks to the following lemma.

**Lemma 2.** The sequence of approximation errors $(V_t)_t$ has an optimal substructure, that is, $V_t = \min_{s=0, \dots, t-1} [V_s + c(\mathbf{y}_{s,t})]$. In other words, the values of $V_t$ can be computed recursively by keeping in memory the values of the subproblems $V_s$, $s < t$.

This optimal substructure is often found in the change-point detection literature, where dedicated optimization algorithms have a complexity of the order of $\mathcal{O}(T^2)$ ($T$ is the signal length) Truong et al. [38]. A similar approach is described in Jewell et al. [39], where the authors try to fit spikes with a parametric shape and undefined length on a univariate signal. However, in our situation, the local nature of the atoms allows us to design a faster procedure. The following proposition builds on Lemma 2 to give an even more straightforward recursive relationship to compute the values of $V_t$.

**Proposition 1.** The optimal approximation error $V_t$ satisfies the following recursion formula: $V_0 = V_1 = \dots = V_{L-1} = 0$ and $\forall t = L, L+1, \dots, T$,

$$V_t = \min [V_{t-1}, V_{t-L} + c(y_{t-L,t})]. \quad (8)$$

Thanks to Proposition 1, the optimal reconstruction error can be efficiently computed. To recover the temporal support of the activations, one only needs to keep track of the argmin in (8).

## 3.3 | Closed-Form Expressions for $c(\cdot)$

To apply our recursive method, the costs $c(\mathbf{y}_{t,t+L})$ must be computed beforehand. We give here the closed-form expressions for two settings: the NoConstraint setting and the AtMostOne-Activation setting.

**Setting NoConstraint.** Under the NoConstraint setting, activations are not constrained (except for the temporal sparsity). Therefore, the cost function $c(\cdot)$ (7) is straightforwardly $c(\mathbf{y}_{s,t}) = \lambda - \|\mathbf{Py}_{s,s+L}\|^2$ if $t-s \ge L$, and $\lambda$ otherwise, where $\mathbf{P} := \mathbf{D}(\mathbf{D}^{\top}\mathbf{D})^{-1}\mathbf{D}^{\top}$. The penalty value $\lambda$ has an easy interpretation in this situation. It acts as a threshold on the norm of the local reconstruction: if the local reconstruction has a norm larger than $\lambda$, it creates a negative cost and decreases the approximation error. Therefore, only large reconstructions are accepted by our algorithm. However, a simple thresholding of the reconstruction norms is not guaranteed to find the optimal temporal support, contrary to our method.

**Setting AtMostOneActivation.** In many applications, at most one atom is activated at each timestamp. This constraint is formally

$$\text{AtMostOneActivation: } \forall t, \|\mathbf{Z}[t]\|_0 \le 1.$$

Note that, under this constraint, the optimal approximation error for the segment $\mathbf{y}_{t,t+L}$ is

$$\min_{\mathbf{u} \in \mathbb{R}^p, \|\mathbf{u}\|_0 \le 1} \|\mathbf{y}_{t,t+L} - \mathbf{Du}\|^2 = \|\mathbf{y}_{t,t+L}\|^2 - \|\mathbf{D}^{\top}\mathbf{y}_{t,t+L}\|_{\infty}^2$$

with $\|\mathbf{x}\|_{\infty} := \max_i |x_i|$. For $(s, t)$ such that $0 \le s < t \le T$, under the AtMostOneActivation setting, the cost function of sub-signal $\mathbf{y}_{s,t}$ is, therefore, $c(\mathbf{y}_{s,t}) = \lambda - \|\mathbf{D}^{\top}\mathbf{y}_{s,s+L}\|_{\infty}^2$ if $t-s \ge L$, and $\lambda$ otherwise. Then Proposition main also holds for this cost function.

## 3.4 | Final Algorithm

The final algorithm, named CSC-$\ell_0$, consists in first computing all the cost values $c(\mathbf{y}_{t,t+L})$ then finding the optimal support using Proposition 1. Under the NoConstraint setting, the complexity of CSC-$\ell_0$ is $\mathcal{O}(LT(\log T + p))$, and under the AtMostOneActivation, it is $\mathcal{O}(pT \log T)$. Details on the complexities, the pseudo-code, and a multivariate extension can be found in the Supporting Information.

Finally, note that we assume a known dictionary in our setting. When this is not verified, several heuristics exist to learn it from the data in an unsupervised way [40, 41]. We describe one in the Supporting Information: Appendix.

## 3.5 | Complexity Analysis

To compute the complexity of CSC-$\ell_0$, we make the standard assumption that the signal length $T$ is much larger than the filter size $L$ and the number $p$ of filters, that is, $T \gg L$ and $T \gg p$. Note that, for all settings, once all the quantities $c(\mathbf{y}_{t,t+L})$ for all
$t$ are precomputed, the complexity of finding the optimal support is the same, namely $\mathcal{O}(T)$.

Under the NoConstraint setting, the complexity of CSC-$\ell_0$ is $\mathcal{O}(LT(\log T + p))$. To simplify, in this setting only, the number $p$ of filters is assumed to satisfy $p \ll L$, meaning that a small number of atoms approximates the signal. First, the complexity of computing the hat matrix $\mathbf{P}$ is $\mathcal{O}(pL^2)$. Second, computing the $\mathbf{Py}_{t,t+L}$ for all $t$ amounts to $L$ sliding correlations along the signals, which results in a complexity of $\mathcal{O}(LT \log T)$, using the Fast Fourier Transform (FFT) algorithm.

Contrary to NoConstraint, in AtMostOneActivation, the number $p$ of filters can be smaller, equal, or even larger than $L$, if the dictionary $\mathbf{D}$ is overcomplete [42] for instance. The complexity is straightforwardly $\mathcal{O}(pT \log T)$, using FFT to compute the correlations with the filters $\mathbf{d}_t$.

## 3.6 | Statistical Analysis

In this section, we show that, in the NoConstraint setting, our method correctly estimates with high probability the number of patterns as well as their localization, under assumptions on the true signal and the noise.

The first ingredient of our analysis is a coherence measure: for two $L \times p$ matrices $\mathbf{A}$ and $\mathbf{B}$ define the subspace coherence $\tilde{\mu}(\mathbf{A}, \mathbf{B})$ between $\mathbf{A}$ and $\mathbf{B}$ by

$$ \tilde{\mu}(\mathbf{A}, \mathbf{B}) := \sup \left\{ \frac{\mathbf{a}^\mathrm{T}\mathbf{b}}{\|\mathbf{a}\|\|\mathbf{b}\|} \mid \mathbf{a} \in \mathcal{R}(\mathbf{A}), \mathbf{b} \in \mathcal{R}(\mathbf{B}) \right\} $$

where $\mathcal{R}(\mathbf{A})$ is the range of $\mathbf{A}$. This quantity, which is equal to the cosine of the smallest principal angle between $\mathcal{R}(\mathbf{A})$ and $\mathcal{R}(\mathbf{B})$ Deutsch [43], is 0 if $\mathcal{R}(\mathbf{A})$ and $\mathcal{R}(\mathbf{B})$ are orthogonal, and 1 if $\dim(\mathcal{R}(\mathbf{A}) \cap \mathcal{R}(\mathbf{B})) \ge 1$. It is used to quantify the “distance” between subspaces Dyer et al. [44], You and Vidal [45]. In addition, for a shift $s = 1, \dots, L-1$, define the shifted dictionary $\mathbf{D}^{(s)} \in \mathbb{R}^{L \times p}$ which is the concatenation of $s$ rows of zeros and the $L-s$ first rows of $\mathbf{D}$. Similarly $\mathbf{D}^{(-s)} \in \mathbb{R}^{L \times p}$ is the concatenation of the last $L-s$ rows of $\mathbf{D}^{(s)}$ and $s$ rows of zeros. See Figure 1 for an illustration. The shifted subspace coherence is then $\tilde{\mu}_s := \tilde{\mu}(\mathbf{D}, \mathbf{D}^{(s)})$ for $s = -L+1, \dots, L-1$. This quantity measures how much a pattern—a linear combination of atoms—can be reconstructed by overlapping patterns. We can now state our assumptions.

**Assumption 1.** The errors $\varepsilon_t$ are i.i.d. sub-Gaussian variables that satisfy $\mathbb{E}\varepsilon_t = 0$ and $\mathbb{E}\varepsilon_t^2 = 1$.

**Assumption 2.** There exists a constant $C_Z > 0$ such that $\|\mathbf{Z}^*[\mathbf{I}]\|^2 \ge C_Z T$ for all $t \in \mathrm{supp}(\mathbf{Z}^*)$.

**Assumption 3.** (i) There exists a positive constant $C_{\mu_0}$ such that the mutual coherence $\mu_0$ of the dictionary $\mathbf{D}$ satisfies $1 - (p-1)\mu_0 \ge C_{\mu_0}$.

(ii) There exist constants $0 < \delta < 1/2$ and $C_{\tilde{\mu}} > 0$ such that the shifted subspace coherence $\tilde{\mu}_s$ satisfies $1 - 2\tilde{\mu}_s^2 \ge C_{\tilde{\mu}}$ for any $|s| > \delta L$.

(iii) The following holds: $C_Z C_{\mu_0} C_{\tilde{\mu}} > (4\sigma)^2$.

Assumption 1 is used to upper-bound the tail of the noise [46]. Short-range dependence could also be assumed, as in Wang and Samworth [47]. Recall that the noise process in Model 1 is $\sigma\varepsilon$, so the noise standard deviation is $\sigma$ since $\mathbb{E}\varepsilon_t^2 = 1$. Assumption 2 is a lower bound on the activation energy, normalized by the number of samples. The last assumption is novel and is part of our contribution. Assumption 3-(i) ensures that an atom cannot be approximated by a linear combination of other atoms. Assumption 3-(ii) introduces $\delta$ which can be seen as a detection margin. It is the smallest shift such that the shifted subspace coherence is bounded away from one. Intuitively, Assumption 3-(ii) ensures that an atom cannot be approximated by a linear combination of two atoms (including itself) that are shifted by more than $\delta L$. In the main result, patterns that are shifted by less than $\delta L$ samples are not distinguishable. In Assumption 3-(iii), the quantity $C_Z C_{\mu_0} C_{\tilde{\mu}} / (4\sigma)^2$ can be seen as a signal-to-noise ratio. There is an interplay between the three parts of Assumption 3; for instance, a large noise level $\sigma$ implies that $C_{\tilde{\mu}}$ must be large and therefore that the detection margin $\delta$ must be large. However, if $\sigma$ is too large, no $\delta$ can satisfy the hypothesis.

To evaluate the pattern detection, we use the Hausdorff distance between sets of events [47]. Denote by $H(A|B)$ the one-sided Hausdorff distance between two sets $A$ and $B$: $H(A|B) := \sup_{b \in B} \inf_{a \in A} |a-b|$. The Hausdorff distance is defined by $H(A, B) := \max[H(A|B), H(B|A)]$. We can now state our main result.

**Theorem 1.** (Non-asymptotic bound). Assume that Assumption 1, Assumption 2 and Assumption 3 hold true. Define the activation map estimate $\hat{\mathbf{Z}}$ given by (3) with the penalty value set to $\lambda = C_{\lambda}T$ where $4\sigma^2 < C_{\lambda} < C_Z C_{\mu_0} C_{\tilde{\mu}}/2$. Denote $\hat{\mathbf{T}} := \mathcal{T}(\hat{\mathbf{Z}})$ and $\hat{\mathbf{K}} := |\hat{\mathbf{T}}|$. Then

$$ \mathbb{P}(H(\mathbf{T}^*, \hat{\mathbf{T}})/L \ge \delta) \le 6e^{-CT \left[ \frac{\sqrt{C_Z C_{\mu_0} C_{\tilde{\mu}}}}{4\sigma} - 1 \right]^2} $$

where the margin $\delta$ is introduced in Assumption 3.

**Corollary 1.** (Asymptotic bound). In the setting of Theorem 1, the following holds true: $\lim_{T \to \infty} \mathbb{P}(\hat{\mathbf{K}} = K^*) = 1$ and $\lim_{T \to \infty} \mathbb{P}(H(\mathbf{T}^*, \hat{\mathbf{T}})/L \ge \delta) = 0$.

Theorem 1 states that with high probability the approximation error of the patterns’ localization is less than the detection margin $\delta$. In addition, the number of patterns is well estimated.

The assumptions and main result are illustrated on a synthetic example shown on Figure 2. We generated two Gabor atoms and computed their shifted subspace coherence (Figure 2b): it is almost 1 for small shifts and 0 for shifts close to $\pm L$. Their mutual coherence is close to 0 (around $10^{-5}$), therefore $C_{\mu_0} \approx 1$. Then, assuming that the minimum activation energy is $C_Z = 1$, we compute the detection margin $\delta$ for several noise levels (Figure 2c). We observe that the margin is around 0.2 (20% of the atoms’ length) for very low noise levels and remains around this value as $\sigma$ goes to 0.10. Then it gradually increases. Interestingly, for noise levels that are similar to the atoms’ amplitude (around 0.08), the detection margin remains close to 20%. For $\sigma$ above 0.25, there is no valid detection margin and no guarantee
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Illustrative example of the shifted subspace coherence and detection margin. (a) Example dictionary with two Gabor atoms. (b) Associated shifted subspace coherence. For a shift *s*, its normalized version is *s*/*L* where *L* is the atoms' length. (c) Evolution of optimal detection margin *δ* w.r.t. noise level *σ*.</figcaption></figure>

can be made about the accuracy of our procedure. The fact that the detection margin does not decrease to 0 as the noise level goes to 0 is surprising. However, our theoretical analysis only provides an upper bound on the minimum detection margin, meaning that the true detection margin could very well be a lot lower (so better). We leave it for future work to investigate if this is a proof artifact or a real phenomenon.

**4 | Experiments**

This section describes several experiments to evaluate the performance of our method for the event detection task, on simulated data and a real-world bio-medical application. Details on the calibration procedures and the dictionary learning method (when it is unknown) can be found in the *Supporting Information* as well as one additional bio-medical application (spike detection in Local Field Potential data).

**4.1 | Baseline Methods**

We compare our method CSC-$\ell_0$ to two common approaches to solving the convolutional sparse coding problem: CSC-$\ell_1$ and CSC-MP. First, CSC-$\ell_1$ solves the same optimization problem as CSC-$\ell_0$ (3) where the penalization $||Z||_{0,1}$ is replaced by $||Z||_{1,1} = \sum_t ||Z[t]||_{1,1}$. This yields a convex minimization problem, which is solved using proximal gradient descent [11] with the implementation of Dupré La Tour et al. [1]. Because this convex relaxation produces many small but non-zero activations, we only keep the local maxima coefficients in absolute value. To prevent overlap, we enforce the local maxima to be located more than *L* samples away from each other, where *L* is the atoms' length. Second, CSC-MP solves (3) in a greedy fashion. It uses the matching pursuit (MP) iterative approach [14]: (i) find the activation vector Z with only one non-zero element ($||Z||_{0,1} = 1$) that best approximates the signal, (ii) remove the reconstruction from the signal and repeat (i) and (ii) on the residual. The algorithm stops when the reconstruction error is decreased by less than a threshold $\lambda$. Again, we impose that two activations are located more than *L* samples away from each other. For all baselines, the atom length *L* is manually chosen, using prior knowledge: it must be as large as the length of the pattern that we are trying to detect, and not too long so that the no-overlap assumption remains satisfied.

**4.2 | Metrics**

The event detection methods are compared using the F1 score. Denote $\hat{t}_1, \dots, \hat{t}_K$ the $\hat{K}$ starts of estimated events and $t_1^*, \dots, t_K^*$

the $K^*$ starts of true events. A true event $t_k^*$ is declared “detected” (or positive) if at least one estimated event is within a user-defined margin. Formally, precision and recall are defined as follows: precision $:= |TP|/\hat{K}$ and recall $:= |TP|/K^*$ where, for a given margin *M*, the true positives TP are true events for which there is an estimated one at less than *M* samples, that is, $TP := \{t_k^* | \exists \hat{t}_l \text{ s.t. } |\hat{t}_l - t_k^*| \le M \}$. We also add the constraint that a single estimated event cannot detect two true events. For the simulated signal, the margin is *M* = 20 samples; for all real-world applications, the margin is *M* = 0.2s.

**4.3 | Semi-synthetic Data Set**

**4.3.1 | Data**

We evaluate the influence of noise and motif overlap on a simulated data set on the performance of CSC-$\ell_0$. The signals are simulated as follows. We choose atoms from a publicly available data set [5]. Those data were collected by monitoring subjects with foot-worn accelerometers while they were walking. Each atom corresponds to a footstep. A synthetic signal is a sequence of 10 randomly chosen atoms from *d* = 14 randomly selected atoms from the data set (atom length *L* = 70 samples). For a signal, the distance between atoms' onsets is constant and chosen between 50 samples and 150 samples. Note that if the distance between atoms is less than *L* = 70 samples, the atoms overlap. We add white noise to the signals, with the signal-to-noise ratio (SNR) varying from -20 to 10 dB. We only evaluate here the accuracy of the sparse coding, therefore, the true dictionary is known by the methods. For a given distance between events and SNR, a training set of 10 signals is simulated to learn the best penalty value $\lambda$ with a grid search. The F1 score is then computed on 100 other simulated signals.

**4.3.2 | Results**

As can be seen on Figure 3a, CSC-$\ell_0$ has a higher F1-score for all SNR, for a fixed distance between events' onsets (80 samples). However, when the SNR is high enough (above -5 dB), all methods perform similarly and have an F1-score close to 100%. Interestingly, Figure 3b shows that performance drops when the distance between events is below 70 samples, that is, when there is overlap. This phenomenon happens for CSC-$\ell_0$, CSC-$\ell_1$ and CSC-MP, as the no-overlap constraint is hard-coded into the algorithms. As soon as the distance between events is above 70 samples, the no-overlap assumption is satisfied, and CSC-$\ell_0$ is close to 100%. CSC-MP has a lower F1 score just after 70 samples but
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 3** | Results on the synthetic data set. Shaded bands show the 5%–95% range. For (a), the distance between events is 80 samples. For (b), the SNR is 0 dB. If the distance between events is less than 70 samples, the atoms overlap.</figcaption></figure>

quickly reaches the level of CSC-$\ell_0$. CSC-$\ell_1$ remains below the two other methods. Overall, the results indicate that CSC-$\ell_0$ is more robust to noise and that all methods strongly rely on the no-overlap assumption.

## 4.4 | Heartbeat Detection in ECGs

### 4.4.1 | Data

Automatic detection of P-wave onsets in electrocardiogram (ECG) signals is crucial for diagnosing cardiovascular disorders, such as various cardiac arrhythmias. Because P-waves are low-voltage phenomena and their shapes are highly variable [4], their detection is challenging and an open research question. We apply CSC-$\ell_0$ on ECGs from the MIT-BIH Arrhythmia Database (MITDB) [19, 48]. For several signals, manual P-Wave annotations have been available. We report the scores on the non-pathological subjects, as in Saclova et al. [4]. Overall, there are five signals with 650,000 samples and 2043 events per signal on average. As in the step detection task, the dictionary is learned using unsupervised dictionary learning. For comparison, a dedicated P-wave detection algorithm, based on wavelets and denoted WT-Detect [49, 50], is also applied to the ECGs.

### 4.4.2 | Calibration

As in the step detection task, the dictionary is learned in an unsupervised way by alternating between sparse coding and

<table><caption>**TABLE 1** | Performance on the heartbeat detection task. Note that the WT-Detect algorithm is specifically designed for this task. Best score out of CSC-$\ell_0$, CSC-MP and CSC-$\ell_1$ is in bold.</caption><thead><tr><th></th><th>CSC-$\ell_0$ (ours)</th><th>CSC-MP</th><th>CSC-$\ell_1$</th><th>WT-Detect</th></tr></thead><tbody><tr><td>Prec.</td><td><strong>96.6 ± 5.5</strong></td><td>85.7 ± 20.3</td><td>83.4 ± 32.5</td><td>99.8 ± 0.2</td></tr><tr><td>Rec.</td><td><strong>100.0 ± 0.0</strong></td><td>98.9 ± 2.1</td><td>83.5 ± 32.5</td><td>99.9 ± 0.1</td></tr><tr><td>F1</td><td><strong>98.2 ± 3.0</strong></td><td>90.5 ± 14.2</td><td>83.4 ± 32.5</td><td>99.8 ± 0.1</td></tr><tr><td>Time</td><td>1.1 s</td><td>1.0 s</td><td>3.8 s</td><td>17.1 s</td></tr></tbody></table>

<table><caption>**TABLE 2** | Performance on the Human Cyclic Activity signals. Best score is in bold.</caption><thead><tr><th></th><th>CSC-$\ell_0$ (ours)</th><th>CSC-MP</th><th>CSC-$\ell_1$</th></tr></thead><tbody><tr><td>AMI</td><td><strong>0.76 ± 0.23</strong></td><td>0.52 ± 0.34</td><td>0.21 ± 0.27</td></tr><tr><td>Time</td><td>0.6 s</td><td>0.001 s</td><td>21.8 s</td></tr></tbody></table>

dictionary learning on the first 30 s of each signal (atom length: 0.5 s). The hyperparameters are found using cross-validation on the five 30-s-long signals for each method. The performance metrics are computed on the remaining portions of signals. For comparison, a dedicated P-wave detection algorithm, based on wavelets and denoted WT-Detect [49, 50], is also applied to the ECGs.

### 4.4.3 | Results

All results are reported in Table 1. CSC-$\ell_0$ has the best performance by a certain margin compared to CSC-$\ell_1$ and CSC-MP. The precision is smaller than the recall, meaning that CSC-$\ell_0$ detects more events than necessary. This phenomenon is more notable for CSC-MP, which causes its F1 score to be 8 points below CSC-$\ell_0$'s. As for CSC-$\ell_1$, both precision and recall are low. Note that ECGs are not necessarily flat between P-waves, which are the motifs of largest variance. CSC-MP and CSC-$\ell_1$ tend to detect the smaller motifs, while CSC-$\ell_0$ focus first on the largest ones. This highlights the importance of computing the optimal activations and not approximating them. CSC-$\ell_0$ has only a slightly worse F1 score than WT-Detect, a dedicated method with close to perfect performance while being unsupervised and faster to execute.

## 4.5 | Human Gait Classification

Human gait is a complex phenomenon that many neurological disorders can alter. Consequently, medical researchers try to objectively measure gait characteristics to analyze patients' state [20]. Gait can be described as a succession of footsteps, and accurately identifying these events is critical. In this section, we investigate the ability of our method to associate a relevant atom to each pattern it detects. This could be used, for instance, to classify footsteps in long times series collected in free-living environments and characterize abnormal gaits.
<figure><img src="image_7.png" /><figcaption>**FIGURE 4** | Signal example from the Human Gait (Joint Angles) data set. Here, the subject does three activities: Walking backward, walking on an inclined floor, and normal walking. The angles between 12 joint pairs are shown in this plot. The repeated patterns are the footsteps. For each algorithm, we report the time positions of the detected patterns: CSC-$\ell_0$ (●), CSC-MP (▲), and CSC-$\ell_1$ (x). The markers have different colors depending on what atom is activated. On this signal, CSC-$\ell_0$ uses a different atom for each activity; CSC-MP uses the same atom for "inclined walking" and "normal walking"; CSC-$\ell_1$ uses the same atom for all activities and misses a few steps at the beginning.</figcaption></figure>

### 4.5.1 | Data

Protocol and data are described in [21]: participants performed cyclic motions, including different forms of walking, while being monitored with sensors, notably markers on the lower limbs, to capture 3D positions. We focus on kinematic data, specifically, joint angles derived from these marker positions. For each leg, measurements include three hip angles (flexion/extension, adduction/abduction, rotation), one knee angle (flexion/extension), and two ankle angles (dorsiflexion/plantarflexion, eversion/inversion), resulting in $\mathbb{R}^{12}$-valued time series.

We create signals by concatenating three walking activities: normal, backward, and on an inclined floor, with variations in speed and direction (up or down the inclined floor). All time series span 60 s, sampled at 200 Hz.

### 4.5.2 | Calibration

For each signal, the dictionary is composed of three atoms, of length 0.8 s. The atoms are randomly taken from the signal but spread evenly so that there is one atom for each activity. We want to measure the ability of algorithms to attribute a different atom to each activity. The algorithm could be used to classify activity depending on the type of pattern it finds. The predicted class of a detected pattern is the index of the atom it was matched with, while the true class is the activity type of the segment it belongs to. We use the Adjusted Mutual Information (AMI) score to compare the two.

### 4.5.3 | Results

As shown on Table 2, our method has the best AMI, followed by CSC-MP, meaning that it can better detect the activity type.

The execution time (0.6 s on average) is larger than CSC-MP's (0.001 s), but far lower than CSC-$\ell_1$ (21.8 s). An illustrating example is shown on Figure 4. On this signal, CSC-$\ell_0$ finds all patterns (footsteps) and associates a different atom with each activity. CSC-MP find most patterns but attributes the same atom to footsteps of normal walking and walking on an inclined plan. CSC-$\ell_1$ merges all types of footsteps into a single atom. On this signal, CSC-$\ell_0$ has an AMI equal to 1, CSC-MP has an AMI of 0.72 and CSC-$\ell_1$ has an AMI of 0.

## 5 | Conclusion

In this work, we have described an algorithm to solve exactly the convolutional sparse coding with a $\ell_0$ penalty for univariate and multivariate signals under a no-overlap assumption. Our procedure is based on a simple recursive relationship and is easy to implement, with linear complexity. We show that the patterns' onsets are detected with high probability, provided that a certain measure of the signal-to-noise ratio is above a threshold. However, the detection is only up to a detection margin, depending on the signal parameters and noise level. Results on simulated data and real physiological time series show that the $\ell_0$-penalized solution is better adapted to detecting characteristic patterns' onsets, compared to classical greedy and $\ell_1$-penalized approaches.

### Author Contributions

T.M. and C.T. designed methods and experiments and wrote the main manuscript. C.T. conducted the theoretical developments, the simulation experiments and provided the software.

### Conflicts of Interest

The authors declare no conflicts of interest.
# Data Availability Statement

The data that support the findings of this study are openly available in A Human Lower-Limb Biomechanics and Wearable Sensors Dataset at https://doi.org/10.35090/gatech/70296.

# References

1. T. Dupré La Tour, T. Moreau, M. Jas, and A. Gramfort, "Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals," in *Advances in Neural Information Processing Systems (Neurips)* (New York, US: Curran Associates, Inc., 2018), 3292–3302.

2. B. J. Owen and B. S. Sathyaprakash, "Matched Filtering of Gravitational Waves From Inspiraling Compact Binaries: Computational Cost and Template Placement," *Physical Review D* 60, no. 2 (1999): 22002.

3. Y. Shi, Z. Nenadic, and X. Xu, "Novel Use of Matched Filtering for Synaptic Event Detection and Extraction," *PLoS One* 5, no. 11 (2010): e15517.

4. L. Saclova, A. Nemcova, R. Smisek, L. Smital, M. Vitek, and M. Ronzhina, "Reliable P Wave Detection in Pathological ECG Signals," *Scientific Reports* 12, no. 1 (2022): 1–14.

5. C. Truong, R. Barrois-Müller, T. Moreau, et al., "A Data Set for the Study of Human Locomotion With Inertial Measurements Units," *Image Processing on Line (IPOL)* 9 (2019): 381–390.

6. M. X. Cohen, *Analyzing Neural Time Series Data: Theory and Practice. Issues in Clinical and Cognitive Neuropsychology* (Cambridge, Massachusetts: MIT Press, 2014).

7. M. Almadani, U. bin Waheed, M. Masood, and Y. Chen, "Dictionary Learning With Convolutional Structure for Seismic Data Denoising and Interpolation," *Geophysics* 86, no. 5 (2021): V361–V374, https://doi.org/10.1190/geo2019-0689.1.

8. C. R. Rojas and B. Wahlberg, "How to Monitor and Mitigate Stair-Casing in L1 Trend Filtering," in *Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, Brisbane, QLD, Australia (pp. 3946–3950, 2015).

9. R. Grosse, R. Raina, H. Kwong, and A. Y. Ng, "Shift-Invariant Sparse Coding for Audio Classification," in *Conference on Uncertainty in Artificial Intelligence (UIAI)* (Vol. 23, pp. 149–158, 2009).

10. Y. Zhu and S. Lucey, "Convolutional Sparse Coding for Trajectory Reconstruction," *IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI)* 37, no. 3 (2015): 529–540, https://doi.org/10.1109/TPAMI.2013.2295311.

11. B. Wohlberg, "Efficient Algorithms for Convolutional Sparse Representations," *IEEE Transactions on Image Processing* 25, no. 1 (2016): 301–315.

12. T. Moreau, L. Oudre, and N. Vayatis, "DICOD: Distributed Convolutional Sparse Coding," in *International Conference on Machine Learning (ICML)*, vol. 80 (Stockohlm, Sweden: PMLR, 2018), 3626–3634.

13. E. Zisselman, J. Sulam, and M. Elad, "A Local Block Coordinate Descent Algorithm for the Convolutional Sparse Coding Model," in *Proceedings of the EEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, Long Beach, CA, USA (pp. 8200–8209, 2018).

14. E. Plaut and R. Giryes, "Matching Pursuit Based Convolutional Sparse Coding," in *Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, IEEE (pp. 6847–6851, 2018).

15. E. Plaut and R. Giryes, "A Greedy Approach to $\ell_0$-Based Convolutional Sparse Coding," *SIAM Journal on Imaging Sciences* 12, no. 1 (2019): 186–210.

16. J. Kim, H. Zhu, X. Wang, and K. A. Do, "Scalable Network Estimation With L0 Penalty," *Statistical Analysis and Data Mining* 14, no. 1 (2021): 18–30, https://doi.org/10.1002/sam.11483.

17. Y. Fu, T. Zhu, Y. Xiang, Z. Chen, and L. Cai, "Stability Analysis of $\ell_0$-Norm Based Convolutional Sparse Coding Using Stripe Coherence," *IEEE Transactions on Signal Processing* 68 (2020): 5810–5823, https://doi.org/10.1109/TSP.2020.3025728.

18. V. Papyan, J. Sulam, and M. Elad, "Working Locally Thinking Globally: Theoretical Guarantees for Convolutional Sparse Coding," *IEEE Transactions on Signal Processing* 65, no. 21 (2017): 5687–5701, https://doi.org/10.1109/TSP.2017.2733447.

19. L. Maršánová, A. Němcová, R. Smíšek, T. Goldmann, M. Vítek, and L. Smital, "Automatic Detection of P Wave in ECG During Ventricular Extrасystoles," in *World Congress on Medical Physics and Biomedical Engineering* (Czech Republic: Prague, 2019), 381–385.

20. L. Oudre, R. Barrois-Müller, T. Moreau, et al., "Template-Based Step Detection With Inertial Measurement Units," *Sensors* 18, no. 11 (2018): 4033.

21. K. Scherpereel, D. Molinaro, O. Inan, M. Shepherd, and A. Young, "A Human Lower-Limb Biomechanics and Wearable Sensors Dataset During Cyclic and Non-cyclic Activities," *Scientific Data* 10 (2023): 1–12, https://doi.org/10.1038/s41597-023-02840-6.

22. P. Mehrgardt, M. Khushi, S. Poon, and A. Withana, "Deep Learning Fused Wearable Pressure and Ppg Data for Accurate Heart Rate Monitoring," *IEEE Sensors Journal* 21 (2021): 27106–27115, https://doi.org/10.1109/JSEN.2021.3123243.

23. T. Germain, S. Gruffaz, C. Truong, A. O. Durmus, and L. Oudre, "Shape Analysis for Time Series," in *Proceedings of Advances in Neural Information Processing System (NeurIPS)* (2024a).

24. T. Germain, C. Truong, L. Oudre, and E. Krejci, "Unsupervised Classification of Plethysmography Signals With Advanced Visual Representations," *Frontiers in Physiology* 14 (2023): 1154328, https://doi.org/10.3389/fphys.2023.1154328.

25. S. Hitziger, M. Clerc, S. Saillet, C. Bénar, and T. Papadopoulos, "Adaptive Waveform Learning: A Framework for Modeling Variability in Neurophysiological Signals," *IEEE Transactions on Signal Processing* 65, no. 16 (2017): 4324–4338.

26. D. Murray, L. Stankovic, and V. Stankovic, "An Electrical Load Measurements Dataset of United Kingdom Households From a Two-Year Longitudinal Study," *Scientific Data* 4 (2017): 160122, https://doi.org/10.1038/sdata.2016.122.

27. A. Petralia, P. Charpentier, P. Boniol, and T. Palpanas, "Appliance Detection Using Very Low-Frequency Smart Meter Time Series," in *Proceedings of the ACM International Conference on Future Energy Systems* (2023).

28. T. Germain, C. Truong, and L. Oudre, "Persistence-Based Motif Discovery in Time Series," *IEEE Transactions on Knowledge and Data Engineering (TKDE)* 36 (2024): 6814–6827.

29. M. Linardi, Y. Zhu, T. Palpanas, and E. Keogh, "Matrix Profile x: Valmod-Scalable Discovery of Variable-Length Motifs in Data Series," in *Proceedings of the ACM SIGMOD International Conference on Management of Data* (2018), 1053–1066.

30. Y. Gao and J. Lin, "Exploring Variable-Length Time Series Motifs in One Hundred Million Length Scale," *Data Mining and Knowledge Discovery* 32 (2018): 1200–1228, https://doi.org/10.1007/s10618-018-0570-1.

31. Y. Gao and J. Lin, "Hime: Discovering Variable-Length Motifs in Large-Scale Time Series," *Knowledge and Information Systems* 61 (2019): 513–542, https://doi.org/10.1007/s10115-018-1279-6.

32. M. Ceccarello and J. Gamper, "Fast and Scalable Mining of Time Series Motifs With Probabilistic Guarantees," *Proceedings of the VLDB Endowment* 15 (2022): 3841–3853.

33. S. Alaee, K. Kamgar, and E. Keogh, "Matrix Profile XXii: Exact Discovery of Time Series Motifs Under Dtw," in *Proceedings of the IEEE International Conference on Data Mining (ICDM)* (Vol. 2020, Nov. 2020).
34. J. B. Malagnoux and M. Kowalski, "From Convolutional Sparse Coding to *-nmf Factorization of Time-Frequency Coefficients," in *Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* (IEEE, 2024): 5530-5534.

35. Y. C. Eldar, P. Kuppinger, and H. Bölcskei, "Blocksparse Signals: Uncertainty Relations and Efficient Recovery," *IEEE Transactions on Signal Processing* 58, no. 6 (2010): 3042-3054, https://doi.org/10.1109/TSP.2010.2044837.

36. Z. Ben-Haim and Y. C. Eldar, "Near-Oracle Performance of Greedy Block-Sparse Estimation Techniques From Noisy Measurements," *IEEE Journal on Selected Topics in Signal Processing* 5, no. 5 (2011): 1032-1047, https://doi.org/10.1109/JSTSP.2011.2160250.

37. J. Wen, H. Chen, and Z. Zhou, "An Optimal Condition for the Block Orthogonal Matching Pursuit Algorithm," *IEEE Access* 6 (2018): 38179-38185, https://doi.org/10.1109/ACCESS.2018.2853158.

38. C. Truong, L. Oudre, and N. Vayatis, "Selective Review of Offline Change Point Detection Methods," *Signal Processing* 167 (2020): 107299.

39. S. W. Jewell, T. D. Hocking, P. Fearnhead, and D. M. Witten, "Fast Nonconvex Deconvolution of Calcium Imaging Data," *Biostatistics* 21, no. 4 (2020): 709-726.

40. T. Moreau and M. Kowalski, "Understanding Approximate and Unrolled Dictionary Learning for Pattern Recovery," in *Proceedings of the International Conference on Learning Representations (ICLR)* (No. 1, pp. 1-21, 2022).

41. O. Permiakova and T. Burger, "Sketched Stochastic Dictionary Learning for Large-Scale Data and Application to High-Throughput Mass Spectrometry," *Statistical Analysis and Data Mining* 15, no. 1 (2022): 43-56, https://doi.org/10.1002/sam.11542.

42. E. Candès, Y. C. Eldar, D. Needell, and P. Randall, "Compressed Sensing With Coherent and Redundant Dictionaries," *Applied and Computational Harmonic Analysis* 31, no. 1 (2011): 59-73.

43. F. Deutsch, *The Angle Between Subspaces of a Hilbert Space* (Netherlands: Springer, 1995).

44. E. L. Dyer, A. C. Sankaranarayanan, and R. G. Baraniuk, "Greedy Feature Selection for Subspace Clustering," *Journal of Machine Learning Research (JMLR)* 14 (2013): 2487-2517.

45. C. You and R. Vidal, "Geometric Conditions for Subspace-Sparse Recovery," in *Proceedings of the International Conference on Machine Learning (ICML)*, Lille, France (Vol. 2, pp. 1585-1593, 2015).

46. R. Vershynin, *High-Dimensional Probability: An Introduction With Applications in Data Science* (Cambridge, USA: Cambridge University Press, 2018).

47. T. Wang and R. J. Samworth, "High Dimensional Change Point Estimation via Sparse Projection," *Journal of the Royal Statistical Society, Series B: Statistical Methodology* 80, no. 1 (2018): 57-83.

48. A. Goldberger, L. Amaral, L. Glass, et al., "PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals," *Circulation* 101, no. 23 (2000): E215-E220.

49. D. Makowski, T. Pham, Z. J. Lau, et al., "NeuroKit2: A Python Toolbox for Neuro-Physiological Signal Processing," *Behavior Research Methods* 53, no. 4 (2021): 1689-1696.

50. J. P. Martinez, R. Almeida, S. Olmos, A. P. Rocha, and P. Laguna, "A Wavelet-Based ECG Delineator: Evaluation on Standard Databases," *IEEE Transactions on Biomedical Engineering* 51, no. 4 (2004): 570-581.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section.