Sebastian Hitziger
email: sebastian.hitziger@gmx.de

Maureen Clerc

Sandrine Saillet

Christian Bénar

Théodore Papadopoulo

Adaptive Waveform Learning: A Framework for Modeling Variability in Neurophysiological Signals

Keywords: dictionary learning, epileptiform discharges, lo-24 cal field potential (LFP), sparse representations, signal variability, 25

When analyzing brain activity such as local field 1 potentials (LFP), it is often desired to represent neural events by 2 stereotypic waveforms. Due to the non-deterministic nature of 3 the neural responses, an adequate waveform estimate typically 4 requires to record multiple repetitions of the neural events. It 5 is common practice to segment the recorded signal into event-6 related epochs and calculate their average. This approach suffers 7 from two major drawbacks: (i) epoching can be problematic, 8 especially in the case of overlapping neural events and (ii) vari-9 ability of the neural events across epochs (such as varying onset 10 latencies) is not accounted for, which may lead to a distorted 11 average. 12 In this paper, we propose a novel method called adaptive wave-13 form learning (AWL). It is designed to learn multi-component 14 representations of neural events while explicitly capturing and 15 compensating for waveform variability, such as changing latencies 16 or more general shape variations. Thanks to its generality, it can 17 be applied to both epoched (i.e., segmented) and continuous (i.e., 18 non-epoched) signals by making the corresponding specializations 19 to the algorithm. We evaluate AWL's performance and robustness 20 to noise on simulated data and demonstrate its empirical util-21 ity on an electrophysiological recording containing intracranial 22 epileptiform discharges (epileptic spikes).

23

 [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF]and independent component analysis (ICA) [START_REF] Comon | Independent component analysis, a new concept?[END_REF]. Especially ICA has proved a valuable tool for separating multi-channel electroencephalogry (EEG) recordings into components representing different active brain sources by assuming their statistical independence [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF], [START_REF] Lagerlund | Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition[END_REF]. Multilinear techniques, such as parallel factor analysis (PARAFAC), allow to decompose multi-channel EEG recordings into components with multiple dimensions, such as space, time, and frequency [START_REF] Miwakeichi | Decomposing eeg data into spacetime-frequency components using parallel factor analysis[END_REF]. However, as the techniques described above rely on a linear framework, they require isochronicity of the multi-variate input signal. This is a reasonable assumption if the input components are synchronously acquired signals from different recording channels, but it typically does not hold for different experimental trials. The more recent method, differentially variable component analysis (dVCA), combines features of PCA/ICA and Woody's method by extending the linear multicomponent framework to include latency variability of each waveform component and has been applied to LFP [START_REF] Truccolo | Estimation of single-trial multicomponent erps: Differentially variable component analysis (dVCA)[END_REF] and EEG recordings [START_REF] Knuth | Differentially variable component analysis: identifying multiple evoked components using trial-to-trial variability[END_REF].

Another approach for analyzing brain signals consists in sparse representations calculated by techniques such as matching pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] or least angle regression (LARS) [START_REF] Efron | Least angle regression[END_REF]. These methods allow the detection of features taken from a dictionary, a predefined and often overcomplete set of atoms (i.e., basis waveforms). MP has been applied to EEG data in [START_REF] Durka | Analysis of EEG transients by means of matching pursuit[END_REF], and extensions have been proposed to specifically address multi-channel [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF], [START_REF]Piecewise linear source separation[END_REF] and multi-trial [START_REF] Bénar | Consensus matching pursuit for multi-trial EEG signals[END_REF], [START_REF] Sielużycki | Single-trial evoked brain responses modeled by multivariate matching pursuit[END_REF] data. A drawback of sparse coding techniques is the fact that the optimal dictionary is often unknown a priori and typical choices such as the symmetric Gabor wavelets [START_REF] Bénar | Consensus matching pursuit for multi-trial EEG signals[END_REF] may not well represent the neural events. A remedy consists in learning the dictionary directly from the data [START_REF] Olshausen | Sparse coding with an overcomplete basis set: a strategy employed by V1?[END_REF], a popular technique especially in the image processing community [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], [START_REF] Mairal | Supervised dictionary learning[END_REF]. Dictionary learning has furthermore been extended to translation-invariant settings [START_REF] Blumensath | Sparse and shift-invariant representations of music[END_REF], allowing to explicitly account for variable latencies of neural events [START_REF] Jost | Learning 1195 redundant dictionaries with translation invariance property: the MoTIF 1196 algorithm[END_REF]- [START_REF] Hitziger | Jitter-adaptive dictionary learning-application to multi-trial 1202 neuroelectric signals[END_REF].

In many neurophysiological applications, it is common practice to segment a long continuous recording into eventrelated epochs in a preprocessing step to facilitate further analysis. However, this epoching step can be problematic if the latencies of the neural events are not exactly known a priori. In addition, in the case of overlapping neural events, epoching may lead to significant errors. Optimally, a method should thus be capable of processing the recording as a whole, which requires a model that allows repetitions of the neural events at different latencies.

In this work, we introduce a new framework, called adaptive waveform learning (AWL), to learn single-or multi-

x m = K k=1 a km d k (• -δ m ) + m , (3)
which is the underlying model of dVCA [START_REF] Truccolo | Estimation of single-trial multicomponent erps: Differentially variable component analysis (dVCA)[END_REF]. Note that in the models above, each component waveform d k occurs at most once per epoch x m . In order to include repetitions, we can add another sum over different translations δ p corresponding to the time samples in each x m ,

x m = K k=1 P p=1 a kpm d k (• -δ p ) + m . (4)
The set of all KP translated waveforms may be very large and overcomplete if KP > T . Hence, the coefficients a kpm should be sparse, i.e., a kpm = 0 for most triplets (k, p, m). This sparse model underlies translation-invariant dictionary learning techniques, where the d k are often called kernels or generating functions. The translation-invariant dictionary then contains all shifted versions of these kernels. An example for an application to EEG recordings is reported in [START_REF] Jost | Learning 1195 redundant dictionaries with translation invariance property: the MoTIF 1196 algorithm[END_REF], where the authors introduce the translation-invariant dictionary learning algorithm MoTIF.

B. AWL model

The idea of the technique presented in this paper is the efficient modeling of the neural events through a small set of kernels d k which are sufficiently adaptive to capture the variability across signal epochs. For this purpose, we extend model (4) by including dilations (i.e., linear temporal scaling), which can account, for instance, for changing signal durations and varying frequencies. This yields

x m = K k=1 P p=1 Q q=1 a kpqm 1 √ γ q d k 1 γ q (• -δ p ) + m . (5)
Note that the idea of multi-scale approaches in dictionary learning is not entirely new: in [START_REF] Mairal | Learning multiscale sparse rep-1205 resentations for image and video restoration[END_REF], the authors learn dictionaries of image patches with blocks of different sizes, and the technique presented in [START_REF] Ophir | Multi-scale dictionary learning using 1208 wavelets[END_REF] makes use of predefined scaleinvariant wavelets. However, in neither approach is the learned dictionary itself scale-invariant.

In model [START_REF] Lagerlund | Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition[END_REF], every neural event is represented as an instantiation of a kernel d k at a specific temporal location δ p and with a specific duration γ q and amplitude a kpqm . This is illustrated in Fig. 1. In addition to capturing the waveform amplitudes, the coefficients a kpqm have a crucial role in selecting the relevant waveforms. That is, each nonzero coefficient a kpqm denotes an occurrence of a neural event given through (δ p , γ q , d k ) in an epoch x m . It is therefore essential for the coefficients to be sparse, since we do not expect neural events to occur at every possible time instant (or with every possible dilation parameter). The specific way of imposing this sparsity on the coefficients will be treated in the following sections.

We note that model [START_REF] Lagerlund | Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition[END_REF], has many unknown parameters (i.e., the kernels d k , their parameters δ p , γ q , and their variable amplitudes a kpqm ), which bears the risk of overfitting the problem and makes interpretation difficult. For the applications considered in this paper, we will therefore never address [START_REF] Lagerlund | Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition[END_REF] in its full complexity, but instead specialize it to different settings.

.........

where the operators φ l may represent translations, dilations,

191
and their compositions, but can also describe more general 197

The abstract model [START_REF] Miwakeichi | Decomposing eeg data into spacetime-frequency components using parallel factor analysis[END_REF] has the advantage that the following 198 analysis is not limited to translations and dilations. In fact, we 199 will only require the operators φ l to be (i) linear

d k 2 = 1, (8)
where • 2 denotes the l 2 -norm. Another indeterminacy con-214 sists in the order of the kernels d k , which will be addressed in 1 Linearity means that

φ l K k=1 a k d k = K k=1 a k φ l (d k )
holds for any sets of kernels {d 1 , . . . , d K } and coefficients {a 1 , . . . , a K }. In particular, translations and dilations are linear operations. 2 Any temporal rescaling can be represented by a linear operator φ.

III. MINIMIZATION PROBLEM AND PROPOSED

ALGORITHMS

Based on model ( 6), we formulate a minimization problem in order to learn the kernels d k , as well as their instantiations in the data, given through the coefficient vector a ≡ {a klm } and the selected operators φ l . For this purpose, we first remain in a general setting, resulting in a template algorithm which will then be implemented for an epoched signal setting (E-AWL) and a continuous (i.e., non-epoched) setting (C-AWL).

A. AWL template algorithm

As discussed in the previous section, model [START_REF] Miwakeichi | Decomposing eeg data into spacetime-frequency components using parallel factor analysis[END_REF] is only useful for interpretation if the coefficient vector a ≡ {a klm } is sparse. This sparsity will be induced through the following exclusivity constraint in order to prevent neural events from being detected multiple times by similar instantiated kernels: we impose that an instantiation of a kernel d k may exclude certain other instantiations of d k (or other similar kernels d k ) in the same signal x m , which can be expressed as a klm = 0 ⇒ a k l m = 0 for appropriate index tuples (k, k , l, l ). In Sections III-B and III-C, we will give concrete implementations of this constraint, which we denote as C(a) throughout this section. Since neural events typically occur with the same polarity within the recordings, it is furthermore reasonable to assume non-negativity of the coefficients a ≥ 0. This reduces the parameter space of the optimization problem, and the following algorithms can ensure this constraint without an increase in computational complexity.

AWL problem: First note that the finite set of operators Φ = {φ l } L l=1 is not directly learned, but instead determined a priori, e.g., as the set of permitted translations and dilations. The relevant operators corresponding to the neural events are implicitly selected from this set through the non-zero coefficients a klm . Hence, the unknowns in model ( 6) are only the coefficient vector a ≡ {a klm } and the kernels d k . Taking into account coefficient sparsity and non-negativity as described above, as well as the normalization of the kernels (8), we formulate the minimization problem

min a,{d k } M m=1 x m - K k=1 L l=1 a klm φ l (d k ) 2 2 , (9)
s.t. d k 2 = 1 for all k, ( 10
) a ≥ 0, (11) C(a), (12)
Note that in the case of Φ = {id} this problem is similar to the dictionary learning problem [START_REF] Olshausen | Sparse coding with an overcomplete basis set: a strategy employed by V1?[END_REF], where in our case, sparsity is induced through C(a). This non-convex joint optimization problem is often solved through alternating minimization: starting with an initial set {d k }, the coefficients a and the kernels d k are iteratively updated in separate steps. We adapt this alternating framework to account for the operators φ l and the constraints [START_REF] Durka | Analysis of EEG transients by means of matching pursuit[END_REF], [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF].

Coefficient update:

The coefficient update consists in the minimization of ( 9)-( 12) with respect to the coefficients a, while leaving the kernels d k and the operators φ l fixed. The operators φ l can thus be eliminated from the cost function by applying them to the kernels d k , i.e., by creating the dictionary

D = {d l k } with atoms d l k def = φ l (d k ).
Since D is fixed and there are no dependencies of the coefficients a across different epochs x m , the resulting minimization can be performed separately for each x m . Hence, for each m = 1, . . . , M , we have to solve

argmin {a••m} x m - K k=1 L l=1 a klm d l k 2 2 , (13)
s.t. a ≥ 0, (14) C(a). (15)
While ( 13) is an ordinary least squares problem, many choices for the exclusity constraint C(a) result in a problem ( 13)-( 15) that is non-convex. For its solution, we will make use of sparse coding techniques. In particular, we shall focus on matching pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and least angle regression shrinkage (LARS) [START_REF] Efron | Least angle regression[END_REF]. The advantage of MP and LARS is that both algorithms iteratively select active atoms (i.e., those with nonzero coefficients) that have maximal dot product with the data. After each selection step, MP subtracts the contribution of the activated atom from the current residual signal and performs the following iteration on the updated residual for the remaining atoms. In contrast, LARS never fully subtracts an atom's contribution. Instead, it keeps track of all activated atoms and can deactivate a selected atom in a later step. While in some cases LARS produces better solutions, it also has a higher computational complexity.

As LARS and MP proceed in successive activation steps, constraints [START_REF] Bénar | Consensus matching pursuit for multi-trial EEG signals[END_REF], [START_REF] Sielużycki | Single-trial evoked brain responses modeled by multivariate matching pursuit[END_REF] can easily be ensured: First, we impose coefficient non-negativity by selecting only atoms which have positive dot product with the data. For LARS, this variant is also mentioned in [START_REF] Efron | Least angle regression[END_REF]. Second, after each activation of some coefficient a klm , we exclude from later selection those coefficients a k l m which would violate C(a).

Note that LARS is typically used to solve the Lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] which includes l 1 -regularization. In the following Sections III-B and III-B we will provide concrete implementations for ( 13)- [START_REF] Sielużycki | Single-trial evoked brain responses modeled by multivariate matching pursuit[END_REF] and will describe how LARS can be used without l 1 -regularization.

Kernel update: Minimizing ( 9)-( 12) for the kernels d k while leaving the coefficients and the operators fixed is a convex problem since the constraints [START_REF] Durka | Analysis of EEG transients by means of matching pursuit[END_REF], [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF] only concern the coefficients a. We can efficiently solve it through block coordinate descent, i.e., by performing loops through the index set k ∈ {1, . . . , K} and minimizing for each d k separately.

For each k, we thus have to solve

argmin d k M m=1 x m - K k =1 k =k L l=1 a k lm φ l (d k ) - L l=1 a klm φ l (d k ) 2 2 , (16)
and rewrite [START_REF] Olshausen | Sparse coding with an overcomplete basis set: a strategy employed by V1?[END_REF] as

310 argmin d k M m=1 x m - K k =1 k =k ψ k m (d k ) -ψ km (d k ) 2 2
.

The ψ k m are fixed and known at this point, so we can 311 differentiate the minimization term with respect to d k and 312 write the necessary condition for a minimum. This yields the 313 closed form solution

314 d k ← M m=1 ψ t km ψ km + M m=1 ψ t km (r km ) , (18)
with ψ t km denoting the adjoint operators, (•) + the Moore-315 Penrose pseudoinverse, and

316 r km = x m - K k =1 k =k a k m ψ k m (d k )
the residual of the signal x m after subtraction of all but 317 kernel d k 's contribution. The kernel update [START_REF] Mairal | Supervised dictionary learning[END_REF] thus describes 318 a generalized average over these residuals, with the adjoint 319 operations ψ t km (r km ) performing realignments. Invertibility 320 and numerical conditioning of the operator ψ t km ψ km are 321 discussed later in the concrete applications.

322

After numerical convergence of the iterative updates [START_REF] Mairal | Supervised dictionary learning[END_REF], 323 which we generally observed already after one cycle through 324 the index set {1, . . . , K}, the d k are normalized in order to 325 ensure constraint [START_REF] Efron | Least angle regression[END_REF]. Note that this normalization of the d k 326 requires a corresponding adjustment of the coefficients, which 327 is automatically done in the next coefficient update.

328

Hierarchical learning: Contrary to other dictionary learning 329 applications where dictionaries are typically large, for the 330 applications addressed in this paper, we are interested in 331 learning only a small number K of kernels. This makes 332 it feasible to hierarchically learn representations of growing 333 cardinalities K: First a representation with a single kernel d 1 334 is learned. Then, a second kernel d 2 is initialized and learning 335 is repeated on the set {d 1 , d 2 }. This process is repeated until a 336 maximal representation size K max is reached. The advantage 337 of this approach is that we obtain a set of representations with 338 different cardinalities, whose comparison can give interesting 339 insight. In addition, it allows us to determine the optimal 340 representation size K a posteriori whose choice is often a 341 difficult task. We only need to ensure that K max is chosen 342 sufficiently large, which may depend on the task and the 343 desired interpretation. For all applications shown in this paper, 344 we found K max = 5 to be sufficient.

345

Note that the hierarchical learning approach also provides 346 an ordering of the kernels, where the last kernels are the ones 347 most recently added to the learning process. This avoids the 348 ordering indeterminacy described in Section II-B.

Implementation: The alternating minimization scheme pre-350 sented above provides the bricks for an algorithm solving (9)-351 [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF]. For a concrete implementation, however, it is necessary to make several additional specifications, including (i) the choice 353 of the operators φ l , (ii) the exact formulation of the exclusivity constraint C(a), (iii) the choice between MP and LARS in the 355 coefficient update, and (iv) the initialization of the kernels d k .

356

The choices for (i)-(iv) should be carefully adapted to the 357 specific applications. This is done for two applications in the 358 following sections, leading to the algorithms E-AWL and C-359 AWL.

Algorithm 1 Hierarchical E-AWL Input: {x m } M m=1 , {δ p } P p=-P , K max ∈ N.

1: for K = 1 to K max do 2:
Initialize d K with white Gaussian noise.

3:

loop 4:
{a kpm } ← COEFF UPDATE({x m }, {δ p }, {d k } K k=1 ).

5:

Drop index p, keeping only the non-zero coefficients 6:

a km and their corresponding latencies δ km .

7:

if stopping criterion reached: break. Create a dictionary D = {d p k } with d p k = d k (• -δ p ).

3:

for m = 1 to M do 4:
Solve through LARS-0: ∀ k : (a k(-P )m , . . . , a kP m ) 0 ≤ 1.

8:

end for 9: end procedure, return a

1: procedure KERNEL UPDATE({x m }, {a km }, {δ km }, {d k }) 2:
for k ∈ {1, . . . , K} do 3: for k = 1 to K do 7:

d k ← M m=1 a km r km (• + δ km ), 4: with r km = x m - k =k a k m (d k (• -δ km )).
d k ← d k (• -δk ), with δk = m akmδkm m akm . 8: d k ← d k / d k 2 . 9:
end for 10: end procedure, return {d k }

B. Epoched AWL

The general AWL problem (9)-( 12) includes the possibility 362 of repeating neural events within a single x m by allowing 363 several instantiations of each kernel d k . In this section, we assume each neural event to occur at most once per signal epoch x m . In addition, we limit the variability to translations about δ p ∈ {δ -P , . . . , δ P } of the kernels. 3 The values for δ -P and δ P determine the maximal shifts to the left and right, respectively, and can be used to control the permitted amount of latency variability. The general AWL problem ( 9)-( 12) can now be specialized to

min a,{d k }    M m=1 x m - K k=1 P p=-P a kpm d k (• -δ p ) 2 2    (19) s.t. d k 2 = 1 for all k, ( 20
) a ≥ 0, ( 21
) (a k(-P )m , . . . , a kP m ) 0 ≤ 1 for all k, m. (22)
Note that the exclusivity constraint C(a) is specified by the l 0 -constraint ( 22) which allows at most one instantiation of each kernel

d k per epoch x m .
Now we can use the alternating minimization scheme from the previous section to derive the concrete hierarchical E-AWL algorithm. Its pseudocode is given in Algorithm 1, followed by the routines COEFF UPDATE and KERNEL UPDATE. In the following, we will discuss this algorithm in detail.

Kernel initialization:

In order to learn the kernels d k blindly with the least possible bias, we suggest random initialization of the d k with white Gaussian noise. Alternatively, inital kernels can be extracted from the data or calculated in a preprocessing step, e.g., by performing a PCA or ICA. However, note that due to the non-convexity of the problem, this bears the risk of converging to a local minimum close to the initialization.

Coefficient update: The minimization of ( 19)-( 22) w.r.t. the coefficients is summarized in the routine COEFF UPDATE, which we solve using a modification of the LARS algorithm denoted as LARS-0. Standard LARS [START_REF] Efron | Least angle regression[END_REF] is designed to solve the Lasso problem

argmin {a••m} x m - K k=1 P p=-P a kpm d p k 2 2 + λ a 1 ,
with λ ≥ 0 denoting a regularization parameter and the l 1norm being defined as a

1 def = k,p,m |a kpm |.
In fact, a special feature of LARS is its ability to calculate the full regularization path, that is, the solution a for any parameter λ ≥ 0. For sufficiently large λ , this solution is a ≡ 0. When decreasing λ , certain entries a kpm in the solution vector a will successively become active, that is, change from zero to non-zero. However, once activated, an entry may become deactivated again on the further regularization path. In our modification LARS-0, we exclude (reinclude) after each activation (deactivation) all entries a kp m , corresponding to translates of the activated (deactivated) kernel, from later activation. This ensures the l 0 -constraint (line 7 in COEFF UPDATE) and is illustrated in Fig. 2. The nonnegativity constraint (line 6) is implemented by only activating Fig. 2: The l0-constraint [START_REF] Hitziger | Jitter-adaptive dictionary learning-application to multi-trial 1202 neuroelectric signals[END_REF] can easily be enforced when performing sparse coding with the LARS algorithm: After each activation of a kernel (green circle) in one LARS step, we exclude all other translations of this kernel from later activation (red crosses). In contrast, if an already active kernel becomes deactivated, its previously excluded translates become available again for activation. entries a kpm if these become positive, otherwise keeping them

ψ k m (d k ) = a k m d k (• -δ k m ).
Since translations are orthogonal operators, we have

ψ t k m ψ k m = a 2 k m
• id, and the update formula [START_REF] Mairal | Supervised dictionary learning[END_REF] further-427 more simplifies to

428 d k ← 1 M m=1 a 2 km M m=1
a km r km (• + δ km ) , where

r km = x m - K k =1 k =k a k m (d k (• -δ km )),
d k ← d k (• -δk ), (23)
where

δk def = M m=1 a km δ km M m=1
a km describes the weighted mean of the previously detected shifts δ km . The new position of the kernel d k thus represents 438 the mean latency of its instantiations, ensuring that it can make optimal use of the permitted latencies {δ -P , . . . , δ P }: Suppose, to the contrary, that a kernel d k is shifted to the left in most of its instantiations (resulting, e.g., from random initialization). Then d k would not be able to detect events located on the extreme left, i.e., farther than the maximally allowed shift δ -P . This issue can be improved through the realignment [START_REF] Mairal | Learning multiscale sparse rep-1205 resentations for image and video restoration[END_REF].

As before, the kernel updates conclude with normalization in order to meet constraint [START_REF] Jost | Learning 1195 redundant dictionaries with translation invariance property: the MoTIF 1196 algorithm[END_REF]. The steps above are summarized in the procedure KERNEL UPDATE.

Note that the kernel realignment and normalization need to be compensated by making corresponding adjustments to the latencies δ km and the coefficients a km , respectively. This is automatically done in the next coefficient update in Algorithm 1, which is why the stopping criterion is placed below the routine COEFF UPDATE rather than KERNEL UPDATE.

C. Continuous AWL

In many applications, the epoched trials addressed in the previous section result from the segmentation of a continuous signal. Such an epoching step can be problematic, especially if the latencies of the neural events are not exactly known or if their waveforms overlap. Therefore, we now present an approach for directly processing a single continuous signal that contains repetitions of neural events.

In the present setting, we will consider both translations and dilations. By making the appropriate specializations to the general AWL problem ( 9)-( 12), we obtain min a,{dk}

x -

K k=1 P p=1 Q q=-Q a kpq 1 √ γq d k 1 γq (• -δ p ) 2 2 (24) s.t. d k 2 = 1 for all k, (25)
a kpq < ρ ⇒ a kpq = 0 for all k, p, q, (26)
a kpq = 0 ⇒ a kp q = 0 if |δ p -δ p | < ∆. (27)
We note that x now denotes a single long signal, whereas the kernels d k are defined on shorter domains. Each translation

d k → d k (• -δ p )
is thus implemented by shifting d k to the time point δ p in the signal domain and then zero-padding. Contrary to the previous section where we limited the maximal shifts δ -P and δ P , we now include translations {δ 1 , . . . , δ P } over the entire signal x, allowing kernels to be instantiated at any time sample. We use logarithmically spaced dilations γ q ∈ {γ -Q , . . . , γ Q }, with maximal compression and stretch given by γ -Q and γ Q , respectively.

Note that we replaced coefficient non-negativity by the stronger constraint [START_REF] Hyvärinen | Independent component analysis: algorithms 1214 and applications[END_REF] with a given threshold ρ > 0. This ensures that events are only detected if the corresponding kernels have sufficiently large correlation with the data.

The constraint [START_REF] Chevallier | On the need for metrics 1216 in dictionary learning assessment[END_REF] implements the exclusivity constraint C(a) for this case. It ensures that different instantiations of a kernel do not fully overlap: with l denoting the length of the kernels, the maximally allowed overlap is (l -∆)/l. This overlap limitation is frequently used in translation-invariant dictionary learning (see for instance [START_REF] Blumensath | Sparse and shift-invariant representations of music[END_REF]) and is due to the fact that waveforms which are slightly shifted are similar to themselves, i.e., have high dot product. Without controlling

Algorithm 2 Hierarchical C-AWL Input: x, {δ p } P p=1 , {γ q } Q q=-Q , 0 < α < 1, K max ∈ N. 1: for K = 1 to K max do 2:
Initialize d K with data segment in x.

3:

loop 4:
a ← COEFF UPDATE(x, {δ p }, {γ q }, {d k } K k=1 , α).

5:

if stopping criterion reached: break.

6:

{d k } ← KERNEL UPDATE(x, a, {δ p }, {γ q }, {d k }). Save representation R K ← (a, {δ p }, {γ q }, {d k } K k=1 ).

9: end for Output: R 1 , . . . , R Kmax .

1: procedure COEFF UPDATE(x, {δ p }, {γ q }, {d k }, α) 2: Create D = {d pq k } with d pq k = 1 √ γq d k 1 γq (• -δ p ) . 3: if max k,p,q d pq k , x < 0 : set d pq k = -d pq k , d k = -d k , ∀k, p, q. 4: Set λ = α • max k,p,q d pq k , x . 5:
Initialize a = 0, I = {(k, p, q)} k,p,q , r = x.

6:

while I = ∅ do 7:

( k, p, q) ← argmax (k,p,q)∈I d pq k , r .

8:

if d pq k , r < λ : break.

9:

Refine dilation γ q .

10:

ak pq ← d pq k , r . 11:
r ← r -ak pq d pq k .

12:

I ← I \ {( k, p, q); |δ p -δ p | < ∆, -Q ≤ q ≤ Q}, 13:
end while 14: end procedure, return a (and {d k } if sign changed, line 3) 1: procedure KERNEL UPDATE(x, a, {δ p }, {γ q }, {d k }) 2:

for k = 1 to K do 3: d k ← (ψ t k ψ k ) + (ψ t k (r k ))
, where 4:

r k = x - k =k a k ψ k (d k ),
ψ k (d k ) = P p=1 Q q=-Q a k pq √ γq d k 1 γq (• -δ p ) . 6:
Align d k w.r.t. prominent landmark.

7:

d k ← 1 √ γk d k 1 γk • , where 8: γk = P p=1 Q q=-Q γ |akpq| q 1 p,q
|akpq| . 9: 499 Kernel initialization: In the continuous case, we generally initialize the kernels with predefined templates. This is necessary because the latencies of the neural events are entirely unknown, making their correct detection a more difficult task than in the epoched case. Thus, initializing with Gaussian noise would bear the risk of only detecting random structures in the data. For the processing of the dataset in Section V, we initialized the kernels with epileptiform spikes taken directly from the data.

d k ← d k / d k 2 .
Coefficient update: For a long signal x with high sampling rate, the set of possible latencies {δ 1 , . . . , δ P } is large, making calculation with LARS impractical. In the continuous setting, we thus use MP for the coefficient updates, which we found to yield very good results in the continuous case. This is due to the fact that when processing a long signal, most of the instantiated kernels have mutually non-overlapping support and thus vanishing dot product. In addition, the constraint [START_REF] Chevallier | On the need for metrics 1216 in dictionary learning assessment[END_REF] further limits the overlap between instantiated kernels. Since MP is exact for orthogonal dictionaries, the error commited by MP is thus relatively low in the present setting.

Our MP implementation as described in COEFF UPDATE successively searches for atoms d pq k that have maximal dot product with the current data residual and subtracts their contribution. It stops when the dot product of every remaining atom with the data is less than the threshold ρ from [START_REF] Hyvärinen | Independent component analysis: algorithms 1214 and applications[END_REF]. To facilitate the choice of ρ, we define it as a fraction 0 < α < 1 of the maximal dot product of all atoms d pq k with the signal x. The parameter α should be chosen dependent on the signalto-noise ratio (SNR), in order to avoid noise fitting. Note that constraint [START_REF] Chevallier | On the need for metrics 1216 in dictionary learning assessment[END_REF] is enforced through the index set I in COEFF UPDATE, which controls the indices of the permitted atoms d pq k .

Note that the dilations are more costly to implement than translations (see Section III-D), and directly using a fine resolution γ q+1 /γ q would be computationally infeasible. Hence, we suggest a multi-resolution approach, initially using a coarse resolution of the set {γ -Q , . . . , γ Q }. After each activation of an atom d pq k in line 7 of COEFF UPDATE, we then refine the corresponding dilation factor γ q (noted in line 9).

Kernel update:

The kernel update is performed by block coordinate descent as described in Section III-A. Now, the operator ψ k from ( 17) is given by

ψ k (d k ) = P p=1 Q q=-Q a k pq √ γq d k 1 γq (• -δ p ) ,
which defines a convolution with a "stretchable" kernel. The update formula [START_REF] Mairal | Supervised dictionary learning[END_REF] reduces to

d k ← ψ t k ψ k + (ψ t k (r k ))
, where

r k = x - K k =1 k =k a k ψ k (d k ).
Hence, the kernel updates are given by a sort of deconvolution of the residual signals r k . The conditioning of the operators ψ k strongly depends on the differences between detected latencies, which can be controlled by the parameter ∆ above.

In case of a poor condition number, regularization should be Fig. 3: K = 3 kernels were defined in order to generate signal epochs (see Fig. 4). They represent different types of activity of interest in neurological recordings. Note that neural background activity is not represented by these kernels but is modeled through pink noise.

considered. However, this did not occur in our experiments, as we chose ∆ sufficiently large.

As in Section III-B, we lift an indeterminacy in the model by realigning the kernels. While we previously used the mean latency across the different epochs for the realignment, this approach is not applicable to the continuous setting where we only have one signal x. Instead, we suggest to align kernels with respect to a prominent landmark, such as the absolute peak of a spike, as done in the following applications.

In addition, the learned kernels should represent the mean duration of their instantiations in the data, in order to make optimal use of the permitted dilations {γ -Q , . . . , γ Q }, cf.

comment after [START_REF] Mairal | Learning multiscale sparse rep-1205 resentations for image and video restoration[END_REF]. Hence, the following rescaling is applied: Dilations were implemented by resampling the discrete signals using linear interpolation; in the downsampling cases, we previously applied an anti-aliasing filter. This implementation means significantly higher computational costs than for translations. In order to still maintain a high resolution between different dilations in Φ, we used the multi-resolution approach described in the coefficient update in Section III-C. Fig. 4: Three randomly chosen trials, generated from the kernels in Fig. 3 according to the model underlying ( 19)- [START_REF] Hitziger | Jitter-adaptive dictionary learning-application to multi-trial 1202 neuroelectric signals[END_REF]. Top row: noiseless trials. Middle row: trials plus pink noise (SNR: 5 dB). Bottom row: trials plus pink noise (SNR: -5 dB). See Section IV-A for more details.

d k ← 1 √ γk d k 1 γk • ,
translation-invariant dictionary learning algorithm MoTIF and 591 independent component analysis (ICA). The performance of 592 C-AWL is demonstrated on real data in the next section.

593

A. Data generation 594

We started by defining K = 3 kernels, representing 5-595 second long signals with 100 Hz sampling rate. They include 596 both transient and oscillatory waveforms (see Fig. 3). These 597 kernels were used to create 200 signal epochs (or trials) 598 according to the E-AWL model underlying the minimization 599 problem ( 19)- [START_REF] Hitziger | Jitter-adaptive dictionary learning-application to multi-trial 1202 neuroelectric signals[END_REF]. Amplitudes and latencies were drawn 600 independently for each kernel from Gaussian distributions with 601 respective means 1 and 0, and respective standard deviations 602 σ a and σ δ specified in the following paragraphs. Negative 603 amplitudes were discarded to ensure constraint [START_REF] Barthélemy | Multivariate temporal dictionary learning for EEG[END_REF]. We 604 simulated pink noise with a 1/f -shaped power spectrum, 605 which is typical for neural background activity. We varied the 606 standard deviation σ of the noise throughout the simulations, 607 resulting in different signal-to-noise ratios (SNR) defined as 608 20 log(σ x /σ ) [dB], with σ x the standard deviation of the 609 simulated (noiseless) signals. Fig. 4 shows three examples of 610 generated trials with different levels of pink noise. Note that 611 the level of latency jitter σ δ here is very low and therefore 612 hardly visible.

613

B. Compared methods

614

The 200 generated signals were processed with MoTIF, 615 ICA, and E-AWL to recover the underlying kernels. In order 616 to be able to compare to the original kernels, we considered 617 the number K to be known a priori.

618

Like hierarchical E-AWL, the translation-invariant MoTIF 619 algorithm [START_REF] Jost | Learning 1195 redundant dictionaries with translation invariance property: the MoTIF 1196 algorithm[END_REF] proceeds by incrementally learning the different 620 kernels. In each such step, the new kernel to be learned is con-621 strained to have minimal cross-correlation with all previously 622 learned kernels, to avoid recovering the same kernel multiple 623 times. In contrast to E-AWL, however, the hierarchy in the 624 approach is strict in the sense that after a kernel is calculated, 625 it is not altered anymore while learning the next kernels. This 626 implies a severe drawback of MoTIF: the first learned kernel 627 naturally captures the maximal variance in the data and is 628 therefore susceptible to contain a linear combination of the 629 original kernels, which cannot be corrected in a later step.

630

For the present comparison, we used the original MATLAB 631 implementation provided to us by the authors of [START_REF] Jost | Learning 1195 redundant dictionaries with translation invariance property: the MoTIF 1196 algorithm[END_REF]. In order 632 to avoid edge effects, we employed zero-padding at both ends 633 of each trial.

634

ICA was calculated using the Matlab software package 635 FastICA 4 described in [START_REF] Hyvärinen | Independent component analysis: algorithms 1214 and applications[END_REF]. As suggested by the authors, we 636 performed a PCA prior to the ICA, in order to whiten the data 637 and reduce its dimension. This distance generalizes the one proposed in [START_REF] Chevallier | On the need for metrics 1216 in dictionary learning assessment[END_REF] by replacing 652 the dot product between kernels by the maximal value of 653 their cross-correlation, thus providing for a shift-invariant 654 measure (see also [START_REF] Lesage | Apprentissage de dictionnaires structurés pour la 1220 modélisation parcimonieuse des signaux multicanaux[END_REF]). Note that the use of the absolute

ε({d k }, { dk }) def = min π∈P(K) 1 K K k=1 ε(d k , dπ(k) ) ∈ [0, 1]. (28)
The calculation of ε can thus be described as finding a pairing 664 of the kernels {d k } with the { dk } such that the average 665 distance between all of these pairs is minimal. Note that in 666 our applications, the number K of kernels is small, such that 667 brute-force minimization over P(K) is a feasible task.

668

Due to its invariance properties, ε is only a pseudo-metric 669 since the separability axiom does not hold. For more informa-670 tion on dictionary metrics, see for instance [START_REF] Chevallier | On the need for metrics 1216 in dictionary learning assessment[END_REF] .

671

D. Quantitative comparisons 672

We investigate the effects of varying kernel amplitudes and 673 latencies, as well as different noise levels on the performances 674 4 http://research.ics.aalto.fi/ica/fastica/ of MoTIF, ICA, and E-AWL, measured by the distance ε between calculated and original kernels. The results are shown in Fig. 5 and discussed below. Note that, contrary to MoTIF and E-AWL, ICA does not explicitly account for varying latencies and is typically used to separate events across different recording channels, where latency jitter is only a minor concern. However, we think that comparison to ICA is instructive, as it shows ICA's tolerance w.r.t. increasing latency jitter.

1) Increasing number of kernels: We first measured the methods' performances for an increasing number of kernels K. For this purpose, we simulated trials with only the first, the first two, and all three kernels shown in Fig. 3. Amplitude and latency variability were fixed to σ a = 0.3 and σ δ = 0.01, respectively (cf. Section IV-A). Pink noise was added, resulting in a signal-to-noise ratio (SNR) of 10 dB. As shown in the upper left plot in Fig. 5, all methods succeed well in recovering a single kernel from the trials (low error ε). However, in the case of two and three kernels, only E-AWL shows good results, both for random and ICA initializations. ICA cannot properly separate the different kernels due to the latency jitter. MoTIF performs even worse at separating several kernels, which is due to the strictly hierarchical learning approach described in Section IV-B. In fact, in the presence of several kernels in the data, the first kernel learned by MoTIF is a linear combination of these kernels (see Fig. 6), which cannot be corrected in a later step.

2) Varying amplitudes: In order to investigate the effect of varying amplitudes, we simulated signals using all K = 3 kernels and setting σ δ = 0 (no latency variability). Amplitude variability σ a was increased throughout the simulations from 0 to 1. Again, the SNR was 10 dB. The resulting errors ε are shown in the upper right panel of Fig. 5. ICA yields very good results for increasing σ a . In fact, this case is optimal for ICA, as the amplitudes were drawn independently which allows ICA to separate the kernels. In addition, there is no latency variability to cope with. In contrast, MoTIF cannot make use of the amplitude variability to separate the kernels: even for high σ a , the first kernel is learned to maximize the data variance and thus contains a mixture of all original waveforms, which cannot be corrected in a later step (cf. Section IV-B). E-AWL initialized with Gaussian noise improves with increasing σ a but does not reach the level of ICA. However, when initialized directly with the ICA components, the E-AWL algorithm converges close to this initialization, as it already provides a very good estimate.

3) Varying latencies: Throughout the next simulations, latency variability σ δ was increased from 0 to 5 seconds, and SNR was kept at 10 dB. In order to allow for the detection of the largely shifted kernels, the length of each trial was extended from 5 to 15 seconds and the permitted translations for MoTIF and E-AWL were increased to ±5 seconds (previously ±0.1 seconds). We maintained some amplitude variability (σ a = 0.3) to allow ICA to separate the waveforms. The results in the lower left plot of Fig. 5 show that ICA's initially good performance decreases with latency jitter above σ δ = 0.004 seconds. Around the same value, the kernel error of E-AWL decreases. In fact, E-AWL does not only compensate for the row) shows a linear mixture of the true kernels (Fig. 3). Since MoTIF learns the kernels strictly hierarchically, this kernel is not corrected when learning the next kernels. The second and third kernels learned by MoTIF do not strongly resemble the original kernels. ICA correctly learns the first original kernel but produces two versions with different phases, resulting from its incapacity to compensate for the latency jitter. In addition, it is not able to correctly separate the second and third original kernels but instead produces a mixture. E-AWL correctly separates all three waveforms. It does so even when initialized with the suboptimal ICA components (last row). However the third learned kernel shows some small baseline change, resulting from fitting low frequencies of the pink noise.

In case of high contamination with pink noise (right half of Fig. 6), MoTIF yields similar results as for high SNR, showing it to be robust against noise. For ICA, the first two kernels are similar to those learned in the high SNR setting. The third kernel, however, seems to capture some low-frequency noise from the pink noise contamination. The kernels learned with E-AWL show strong contaminations with low-frequency noise. E-AWL's ability to compensate for varying latencies makes it susceptible to fitting the low frequency components in the pink noise. Even when initialized with the ICA kernels, E-AWL still strongly picks up this noise.

V. APPLICATIONS TO NEUROLOGICAL SIGNALS

In this section, we demonstrate the usefulness of the AWL framework as a data exploration tool capable of producing compact and insightful representations. For this purpose, we apply AWL to a neuroelectrical recording containing epileptiform discharges or spikes using two approaches. The first approach requires a prior segmentation step to produce a set of short signal epochs which are then processed with E-AWL. Such an epoched approach is frequently used in neurological signal processing and allows us to compare to other methods that require multiple input signals, such as ICA and MoTIF. Finally, we demonstrate how the single, continuous (i.e., nonepoched) recording can be directly processed with C-AWL and illustrate the complementary benefits of this second approach.

A. Data acquisition

In an animal model of epilepsy, an electrode was placed in the cortex of a Wistar-Han rat for measuring local field Learned kernels (pink noise, SNR 5 dB) Learned kernels (pink noise, SNR -5 dB)

Fig. 6: Kernels recovered from trials contaminated with pink noise (Fig. 4). For high SNR of 5 dB (left), E-AWL shows best performance and correctly identifies the three original kernels (Fig. 3), both for initialization with random noise (third row) and with ICA (bottom row). MoTIF (top row) and ICA (second row) do not separate all kernels correctly. For low SNR of -5 dB (right), the performances of MoTIF and ICA only slightly worsen w.r.t. the high SNR setting. The kernels learned with E-AWL, however, show strong contamination with low-frequency noise. an inhibition blocker (bicuculline) had been injected into the cortex to provoke epileptiform discharges. This data acquisition was performed simultaneously with other multi-modal recordings and the full experimental protocol can be found in [START_REF] Saillet | Interneurons 1224 contribute to the hemodynamic/metabolic response to epileptiform dis-1225 charges[END_REF].

As can be seen in Fig. 7, the spiking activity changed throughout the recording, with periods of high and low spiking densities and different spike amplitudes. However, the exact spike shapes and their evolution across the dataset cannot be directly seen in this plot. Therefore, the goal of the following analysis was to obtain a compact representation of this LFP dataset, which could provide insight into the spike shapes as well as their variable parameters, such as amplitudes, durations, and spiking rates.

B. Epoched processing

In the first approach to process the LFP signal, 169 spikes with at least 10-second inter-spike intervals (peak-to-peak)

were manually selected and segmented into 10-second time windows, centered around the spikes. The time windows were chosen relatively large w.r.t. the duration of the spikes in order to possibly recover other signal structures in their vicinity. In fact, the identification of an oscillatory artefact as shown in Fig. 9: Kernels learned with E-AWL from 169 spike epochs (Fig. 8). Each of the first five rows corresponds to one set of normalized kernels obtained in the hierarchical learning approach. Note that only the central 3 seconds of the 10-second long kernels are shown. The two bottom rows show the representations learned non-hierarchically, using white Gaussian and ICA initialization, respectively. The last column shows the coefficients of each kernel (row) used in each epoch (column); light colors correspond to large values. The red curve under each kernel shows the distribution of latencies used for this kernel across epochs. Note that these distributions are shifted 0.5 seconds to the right for better visualization, avoiding overlaps with the spikes.

the following paragraphs would not have been possible on 838 short epochs. Fig. 8 shows seven sample epochs, as well as 839 the average over all 169 epochs.

840

Processing with E-AWL: Hierarchical E-AWL (Algo-841 rithm 1) was used to learn kernel representations of increasing 842 cardinalities K = 1, . . . , 5. In order to enable blind learning of 843 interesting signal structures, we initialized each newly added 844 kernel with Gaussian noise. We compared this hierarchical 845 approach to the direct learning of the K = 5 kernels. In order 846 to allow E-AWL to identify waveforms with large jitter or 847 phase variability, we used translations {δ -P , . . . , δ P } from -2 848 to 2 seconds.

849

For each K = 1, . . . , 5, the corresponding kernel represen-850 tation of hierarchical E-AWL is shown in the first five rows 851 of Fig. 9 and consists of the K kernels (black signals), their latency distributions (red curves below the kernels), and their coefficients across the epochs (grey values in the last column, see caption). The two bottom rows show the representations learned with non-hierarchical approaches, where the kernels were initialized with Gaussian noise and with ICA, respectively. Note that only 3 seconds of the 10-second long kernels are shown as the remainders of the kernels did not contain any interesting information.

For K = 1, the resulting kernel is simply a weighted average across epochs and resembles the average spike shown in Fig. 8. In fact, since the prior manual epoching step accurately aligned spikes across epochs, this kernel's latency changes were negligible, as reflected by the sharply peaked latency distribution. Adding a second kernel results in two spike components (second row). Only after learning the third kernel does an entirely new, oscillatory waveform appear.

These oscillations were later identified to be an artefact from the recording device. With the fourth and fifth learned kernels, the spike is further refined into different components. While all kernels representing spike components almost always have zero latency, the oscillatory kernel takes different latencies mainly in a range of about 1 second, corresponding to its period. This dispersed latency distribution indicates the independence of the waveform's phase w.r.t. the positions of the spikes.

Comparing the last three rows of Fig. 9 shows that the kernels produced by E-AWL differ depending on the learning approach (hierarchical vs. non-hierarchical) and the used initialization (random vs. ICA). However, the results are qualitatively similar, each containing four spike components and one periodic waveform.

The learned kernel coefficients (last column of Fig. 9) provide insight into the evolution of the spikes across the recording. For K = 1, the coefficient profile shows decreasing energy across the epochs. Interestingly, for K = 2, the coefficients of the second learned spike component (second row) only take non-zero values after the first 60 epochs, indicating a sudden change in the spike shape. For K > 2, the coefficient profiles reveal even more detailed structural information about the spike's evolution. These profiles may be taken as an indicator for the optimal number of kernels to be learned: while the profiles for K < 4 look relatively smooth, we see more frequent changes in the coefficients for resent the different phases of the oscillations through linear 910 combinations of the differently shifted sinusoidal kernels 4 911 and 5. Second, these sinusoidal functions only capture the 912 fundamental frequency of the oscillations and do not show the 913 distinctive pointy shape of the oscillatory waveform clearly 914 visible in the E-AWL representation. Third, only E-AWL 915 clearly separates the oscillatory waveform from the spike 916 components, while the sinusoidal kernels of ICA contain spike 917 artefacts.

K = 5,

918

Note that the original kernels learned with MoTIF contained 919 the spike components at arbitrary temporal locations, due to 920 the translation-invariance in the approach (we only aligned 921 them here for better visualization).

922

The learned coefficients (last column of Fig. 10) also reveal 923 important differences between the three methods. For MoTIF, 924 we can observe very similar coefficient profiles for most 925 kernels. Only the second kernel shows a very low profile, 926 suggesting that it does not well capture an actual pattern in 927 the data. In fact, the shape of the kernel seems to contain 928 artefacts, possibly resulting from MoTIF's maximal decorre-929 lation constraint (cf. Section IV-B). In the case of ICA, the first 930 three kernels are active together in the first half of the epochs. 931 The coefficient profile of E-AWL appears more contrasted and 932 provides a detailed structuring of the epochs.

933

For a quantitative comparison, we calculated the distances ε 934 as defined in [START_REF] Lesage | Apprentissage de dictionnaires structurés pour la 1220 modélisation parcimonieuse des signaux multicanaux[END_REF] between the kernel sets learned with MoTIF, 935 ICA, and the three different E-AWL approaches (hierarchical 936 and random vs. ICA initialization). These distances are visu-937 alized in the matrix in Fig. 11. The MoTIF kernels have the 938 largest distance to the kernels obtained with the other methods. 939 The smallest distances are found between the different E-AWL 940 approaches.

941

C. Continuous processing 942

The epoched approach above suffers from several draw-943 backs. First, the manual epoching is time-consuming and 944 would not be feasible for a larger set of recordings. Second, 945 this approach requires spikes to be well isolated, which was 946 only the case for a small subset of spikes in the given data. 947 We now demonstrate how the recording can be processed with 948 C-AWL without prior epoching. For the kernels, we used 1.5-second long time windows.

We set the parameter ∆ = 0.2 seconds (cf. Section III-C), resulting in a maximal overlap between kernel instantiations of 87%. This constraint prevented spikes from being detected multiple times. At the same time, ∆ was chosen small enough to still allow the detection of spikes in close succession.

The relative correlation threshold was set to α = 0.1. For the presented dataset, which had a high SNR, the small value for α allowed us to even detect low-amplitude spikes. In cases of lower SNR, α should be chosen larger to avoid noise fitting.

We applied C-AWL in two different ways to explain the spike variability: (1) using a multi-class model with different constant kernels d k and (2) using a single-class model with one kernel d of adaptive duration. Note that both are special cases of the C-AWL model. We found that using different kernels and variable duration in a single approach provided too many parameters to describe the spike variability and led to redundancies in the representation.

In order to verify the performance of both approaches, the spikes were first detected manually (n = 520) and their temporal locations were compared to the ones detected with C-AWL. Note, however, that the main objective of this section is to demonstrate the qualitative advantages of C-AWL compared to the epoched approach from Section V-B. For a more exhaustive quantitative evaluation we refer the reader to Chapter 6 in [START_REF] Hitziger | Modeling the variability of electrical activity in the brain[END_REF], where C-AWL is compared to template matching in terms of detection performance for different noise levels.

1) Multiple kernels of constant durations: We learned hierarchical representations with K = 1, . . . , 5 kernels using Algorithm 2 without dilation-invariance. Here, we only analyze the representation for K = 5, which is illustrated in Upper left: the five learned kernels. Middle left: spike coefficients plotted in time across the recording. Lower left: spike coefficients plotted against the temporal distances to previous spikes, with the dashed line describing the system's maximal spiking potential. Right: each learned kernel is plotted on top of the spikes that it represents in the data. Fig. 12, the colors correspond to the different kernels. We 997 see an overall decrease in spike energy to about one third 998 of the initial energy towards the end of the recording. This 999 is more than the decrease of the spike peaks from about 1000 -700 µV to -400 µV, which we can observe in the original 1001 recording in Fig. 7. This suggests that the decreasing spike 1002 energies result not only from their different amplitudes but 1003 also their different durations. Besides the global decrease, the 1004 coefficient plot shows well-separated clusters corresponding to 1005 the spike classes, which provides an interesting structuring of 1006 the dataset. Around 0, 700, 1250, and 1950 seconds, we can 1007 see clusters of slightly smaller coefficients. These correspond 1008 to the periods of dense spiking activity, which can be directly 1009 observed in the original recording in Fig. 7.

1010

The relationship between the spike coefficients (i.e., their 1011 l 2 -norms) and the spiking density becomes even clearer from 1012 the lower left of Figure 12, where we plotted each coefficient 1013 against the inter-spike delay w.r.t. the preceding spike (log 1014 scale). We can observe that the coefficients are larger for 1015 longer inter-spike intervals, suggesting that the system requires 1016 some time to regain its full spiking potential after each spike. 1017 In fact, the dashed line clearly shows this maximal spiking 1018 potential as a function of the inter-spike intervals.

1019

Comparison with the manually detected spikes showed that 1020 all 518 spikes were true positives. Only 2 spikes were missed 1021 by the C-AWL algorithm, i.e., 100% precision and 99.6% 1022 recall.

1023

Note that the processing of the recording with C-AWL using 1024 different kernels of constant durations showed to result in a 1025 combined spike detection and clustering algorithm. This is 1026 similar to an approach recently proposed in [START_REF] Nonclercq | Cluster-based spike detection algorithm 1230 adapts to interpatient and intrapatient variation in spike morphology[END_REF]. However, 1027 the latter algorithm requires the choice of several correlation 1028 and feature thresholds for spike detection, whereas C-AWL 1029 uses only the correlation threshold α. Another advantage of C-1030 AWL is the possibility to describe the spike variability directly 1031 through a dilation parameter, as demonstrated in the following 1032 epoched (E-AWL) and continuous single-channel recordings (C-AWL). Through the explicit modeling of waveform variability, AWL is capable of capturing variations across the recorded neural events, such as different amplitudes, latencies, and dilations.

The application to recorded local field potentials (LFP) containing epileptiform discharges showed the capability of both E-AWL and C-AWL to learn interesting data representations. In fact, due to their complementary approaches, the two algorithms provided very different insights into the recording: Using previously epoched spike segments, E-AWL produced detailed decompositions of the spikes into several components. In addition, it revealed a hidden oscillatory artefact in the data. In turn, C-AWL did not require the time-consuming epoching step but was able to automatically detect the spikes in the continuous signal. This gave a more complete representation of the dataset since C-AWL detected even close spike occurrences, which had to be omitted in E-AWL. In summary, E-AWL proved capable of revealing the patterns constituting a signal, while C-AWL is better designed to detect given patterns inside a signal. This suggests a combination of both methods for a fully automatic pattern recognition methodology in future works. For all experiments, we were able to use high resolution across translations and dilations, thanks to an efficient implementation using the fast Fourier transform and a multi-resolution approach.

The general AWL framework furthermore allows to implement other types of waveform variability. For instance, we found that dilations (together with varying amplitudes) could account for the majority of the spike variability, but not for all of it. More general temporal rescaling functions could thus be considered, for example, through the use of dynamic time warping. Besides the processing of epileptiform spikes, E-AWL and C-AWL can also be applied to other neurophysiological signal processing tasks, such as: (i) the identification of event-related brain potentials (ERPs) across a set of experimental trials as well as the description of the inter-trial (or inter-subject) variability with E-AWL and (ii) the automatic detection of spontaneously occurring neural events with C-AWL, such as sleep spindles during stage 2 sleep.

As a drawback, we saw that E-AWL's ability to compensate for latency variability makes it susceptible to fitting lowfrequency noise components. In cases of strong low-frequency noise, high-pass filtering, either as a preprocessing step or as a postprocessing of the learned kernels, could therefore be considered. The susceptibility to low-frequency noise furthermore shows that the complexity of the AWL framework (type and amount of permitted variability, number of kernels, additional constraints) should be carefully adapted to each application and the signal-to-noise ratio (SNR).

Currently, AWL is being extended to process multi-channel recordings (e.g., EEG), based on the observation that different channels can be treated similarly to the different trials in E-AWL. However, latencies should only vary across trials, since neural events typically appear across channels without time delay. This is similar to the multi-channel extension of dVCA proposed in [START_REF] Knuth | Differentially variable component analysis: identifying multiple evoked components using trial-to-trial variability[END_REF]. Preliminary work on multi-channel AWL has recently been presented at the International Conference on

Fig. 1 :( 5 a

 15 Fig.1: The AWL framework (5) models each neural event in the signals xm as the instantiation of a kernel d k through a specific latency δp, duration γq, and amplitude a kpqm . Note that a kernel may be used multiple times in the same signal xm to model repeating events (bottom row).

192 morphological deformations. Note that the model parameters 193 in ( 6 )

 1936 that need to be learned are only the coefficients a klm 194 and the kernels d k , whereas the finite set {φ l } L l=1 is defined 195 a priori. As mentioned above, the non-zero coefficients a klm 196 thus fulfil the role of selecting the relevant operators φ l .

  while leaving all d k with k = k fixed. Since we assume the operators φ l to be linear (cf. Section II-B), we can define operators 309

  407 zero (cf. Section III-A). Since the problem in lines 5-7 of 408 COEFF UPDATE does not contain an l 1 -regularization term, 409 we calculate the regularization path until λ = 0. While this 410 use of LARS may seem unconventional, it has two important 411 advantages: (i) While the unconstrained problem in line 5 of 412 COEFF UPDATE could also be solved by an ordinary least 413 squares solver, ensuring the additional constraints in lines 6,7 414 is non-trivial, but can be conveniently handled in LARS's reg-415 ularization path. (ii) When considering only the least squares 416 problem in line 5, following the LARS path until λ = 0 does 417 in fact provide the exact solution, whereas matching pursuit 418 (MP) would only calculate an approximation. 419 Kernel update: As ensured by constraint (22), there is 420 maximally one non-zero coefficient a kpm per epoch x m and 421 kernel d k . We can thus drop the index p, denoting by a km 422 only the non-zero coefficients and by δ km the corresponding 423 latencies. The operators ψ k m defined in the previous section 424 in (17) thus reduce to 425

resulting in realigned averages of the residual signals r km . 429

 429 Note that the absolute temporal positions of the kernels 430 are arbitrary in the sense that we could obtain an equivalent 431 representation by slightly shifting a kernel d k and correspond-432 ingly adjusting the detected latencies δ km . In order to lift 433 this indeterminacy, we fix the absolute position of each kernel 434 through the alignment 435

  and

11 :

 11 end procedure, return {d k } the maximal overlap, a neural event might thus be encoded by 489 several slightly shifted versions of the same kernel. This not 490 only complicates interpretation, it also makes the following 491 algorithm less stable (see kernel updates). In cases where the 492 kernels are very similar, such as the spike classes learned 493 in Section V-C, we also limit the overlap between different 494 kernels by replacing the second index k by k in (27). Problem (24)-(27) can again be solved by implementing 496 the alternate minimization scheme from Section III-A, The 497 pseudocode is given in Algorithm 2 and discussed in the 498 following paragraphs.

  where γk is the geometric mean of the dilations used in the instantiations of d k . As before, the kernel update is concluded by normalizing each d k . The steps above are summarized in the routine KERNEL UPDATE. D. Implementation details Both LARS and MP are based on the dot products between the atoms and the data. In case of translated kernels d k (•-δ p ), this requires the computation of cross-correlations which can be efficiently calculated through the fast Fourier transform. This efficient calculation allows us to use a resolution δ p+1 -δ p equal to the sampling resolution of the signals (both for E-AWL and C-AWL).

  Both E-AWL and C-AWL have been implemented in C++ with MATLAB interface (mex-files). The code for these implementations and the following experiments are freely available at https://github.com/hitziger/AWL.IV. SYNTHETIC EXPERIMENTSWe use simulated data to evaluate the capability of the E-AWL algorithm to identify three kernels from a set of signals in the presence of amplitude and latency variability as well as noise. The results are compared to those obtained by the

638E-

  AWL was implemented according to Algorithm 1. We 639 used two different initializations to compare their impact on 640 the learned kernels: (i) random Gaussian noise and (ii) the 641 kernels obtained with ICA. To distinguish these two initializa-642 tions, approach (ii) is denoted as ICA + E-AWL.643For both MoTIF and E-AWL, we allowed translations 644 {δ -P , . . . , δ P } ranging from -0.1 to 0.1 seconds, with a 645 resolution δ p+1 -δ p equal to the sampling period (0.01 s).646C. Kernel distances 647In order to quantify the methods' performances, we defined 648 a distance between the original kernels and the calculated ones. 649 For this purpose, let ε first denote the distance between two 650 normalized kernels d and d, given by 651 ε(d, d) def = 1 -max t τ d(τ ) d(τ + t) .

655

  value furthermore yields sign-invariance. Both properties are 656 important in the present setting, due to indeterminacies in 657 relative latencies and signs of the calculated kernels. Another 658 indeterminacy consists in the order of the learned kernels, 659 which needs to be accounted for when extending ε to measure 660 the distance between kernel sets. For this purpose, let P(K) 661 denote the set of permutations of {1, . . . , K}. For two sets of 662 normalized kernels {d k } and { dk }, we can now define 663

Fig. 5 :

 5 Fig.5: Performances of MoTIF, ICA, and E-AWL in recovering kernels from trials, measured by the distance ε between original and calculated kernels. Four different settings where studied (see plot titles). For E-AWL, we initialized with both Gaussian noise and ICA components (the latter is denoted by ICA + E-AWL, see legend in first plot). See Section IV-D for a detailed discussion.

744 4 )Fig. 7 :

 47 Fig.7: Recording of local field potentials (LFP) in the cortex of a rat. The vertical lines are epileptiform discharges (spikes), whose density visibly changes throughout the recording. The maximal negative spike amplitudes decrease towards the end of the recording from about -700 µV to -400 µV.

Fig. 8 :

 8 Fig.8: The first seven plots show sample epochs from the 169 epoched data segments (original order in data maintained). The average over all 169 epochs is plotted in the bottom right. It is apparent that spikes are changing throughout the dataset, decreasing in duration and amplitude. Note that for better visualization only 3 seconds of the 10-second long epochs are shown.

Fig. 10 :

 10 Fig. 10: Five kernels learned with MoTIF, ICA, and hierarchical E-AWL, see respective rows. The last column shows the absolute coefficients of the kernels (rows) used across the epochs (columns). Light colors correspond to large absolute values. While ICA and E-AWL both recover an oscillatory artefact, only E-AWL clearly separates it from the spike components and encodes it in a single kernel.

Fig. 11 :

 11 Fig.11: Distance ε between the kernels learned with MoTIF, ICA, and E-AWL. For the latter, we used three different approaches: hierarchical and non-hierarchical learning with both random and ICA initialization of the kernels.

Fig. 12 .Fig. 12 :

 1212 Fig. 12. The five learned kernels are shown in the upper left plot, where they are scaled with the average coefficients of their respective occurrences in the data. Note that the time window was chosen sufficiently large to learn not only the first negative wave but also the slow positive wave following it. We can see that the spike classes represented by the kernels differ mainly in duration and average amplitude. The plot on the right shows the negative waves of the learned kernels plotted on top of the respective spikes they represent in the recording. Note the sharp kink around 0.03 seconds in the first three kernels, which can also be observed in the real spikes and gives evidence of the good time resolution properties of C-AWL. The coefficients of the 518 detected spike occurrences are plotted in time across the recording in the middle left of

Fig. 13 :1033 2 )

 132 Fig. 13: Spike representation learned with C-AWL using a single kernel with adaptive duration. Upper left: the learned kernel, scaled with the mean amplitude of its instantiations. Middle left: spike coefficients plotted in time across the recording. Lower left: dilation factors plotted in time across the recording. Right: the top plot shows the superposed detected spikes and the learned kernel. The four bottom plots illustrate samples of detected spikes with corresponding instantiations of the kernel.

  KERNEL UPDATE({x m }, {a km }, {δ km }, {d k }). Save representation R K ← ({a km }, {δ km }, {d k }) K k=1 .

	8: {d k } ← 9: end loop
	10:
	11: end for
	Output: R 1 , . . . , R Kmax .

1: procedure COEFF UPDATE({x m }, {δ p }, {d k }) 2:

the hierarchical learning approach at the end of Section III-A.

This model already proved sufficiently rich for the following epoched applications, and adding dilation invariance did not provide better results. Dilation invariance is therefore explicitly studied only in the continuous setting in Section III-C.

Basic and Clinical Multimodal Imaging. ACKNOWLEDGMENT This work was supported by a doctoral grant of the region Provence-Alpes-Côte d'Azur and the ANR grants COADAPT (09-EMER-002-01), MULTIMODEL (2010-BLAN-0309-04), and VIBRATIONS (ANR-13-PRTS-0011). The electrophysiological recordings (Section V) were obtained within the MULTIMODEL project.