Learning the Morphology of Brain Signals Using Alpha-Stable Convolutional Sparse Coding

by   Mainak Jas, et al.
Télécom ParisTech

Neural time-series data contain a wide variety of prototypical signal waveforms (atoms) that are of significant importance in clinical and cognitive research. One of the goals for analyzing such data is hence to extract such 'shift-invariant' atoms. Even though some success has been reported with existing algorithms, they are limited in applicability due to their heuristic nature. Moreover, they are often vulnerable to artifacts and impulsive noise, which are typically present in raw neural recordings. In this study, we address these issues and propose a novel probabilistic convolutional sparse coding (CSC) model for learning shift-invariant atoms from raw neural signals containing potentially severe artifacts. In the core of our model, which we call αCSC, lies a family of heavy-tailed distributions called α-stable distributions. We develop a novel, computationally efficient Monte Carlo expectation-maximization algorithm for inference. The maximization step boils down to a weighted CSC problem, for which we develop a computationally efficient optimization algorithm. Our results show that the proposed algorithm achieves state-of-the-art convergence speeds. Besides, αCSC is significantly more robust to artifacts when compared to three competing algorithms: it can extract spike bursts, oscillations, and even reveal more subtle phenomena such as cross-frequency coupling when applied to noisy neural time series.



page 1

page 2

page 3

page 4


Parameter Estimation of Heavy-Tailed AR Model with Missing Data via Stochastic EM

The autoregressive (AR) model is a widely used model to understand time ...

General Convolutional Sparse Coding with Unknown Noise

Convolutional sparse coding (CSC) can learn representative shift-invaria...

Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals

Frequency-specific patterns of neural activity are traditionally interpr...

Deep Residual Auto-Encoders for Expectation Maximization-based Dictionary Learning

Convolutional dictionary learning (CDL) has become a popular method for ...

A Nonparametric Frequency Domain EM Algorithm for Time Series Classification with Applications to Spike Sorting and Macro-Economics

I propose a frequency domain adaptation of the Expectation Maximization ...

A Robust Score-Driven Filter for Multivariate Time Series

A novel multivariate score-driven model is proposed to extract signals f...

Variational inference for rare variant detection in deep, heterogeneous next-generation sequencing data

The detection of rare variants is important for understanding the geneti...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Neural time series data, either non-invasive such as electroencephalograhy (EEG) or invasive such as electrocorticography (ECoG) and local field potentials (LFP), are fundamental to modern experimental neuroscience. Such recordings contain a wide variety of ‘prototypical signals’ that range from beta rhythms (12–30 Hz) in motor imagery tasks and alpha oscillations (8–12 Hz) involved in attention mechanisms, to spindles in sleep studies, and the classical P300 event related potential, a biomarker for surprise. These prototypical waveforms are considered critical in clinical and cognitive research [1], thereby motivating the development of computational tools for learning such signals from data.

Despite the underlying complexity in the morphology of neural signals, the majority of the computational tools in the community are based on representing the signals with rather simple, predefined bases, such as the Fourier or wavelet bases [2]. While such bases lead to computationally efficient algorithms, they often fall short at capturing the precise morphology of signal waveforms, as demonstrated by a number of recent studies [3, 4]. An example of such a failure is the disambiguation of the alpha rhythm from the mu rhythm [5], both of which have a component around  Hz but with different morphologies that cannot be captured by Fourier- or wavelet-based representations.

Recently, there have been several attempts for extracting more realistic and precise morphologies directly from unfiltered electrophysiology signals, via dictionary learning approaches [6, 7, 8, 9]. These methods all aim to extract certain shift-invariant prototypical waveforms (called ‘atoms’ in this context) to better capture the temporal structure of the signals. As opposed to using generic bases that have predefined shapes, such as the Fourier or the wavelet bases, these atoms provide a more meaningful representation of the data and are not restricted to narrow frequency bands.

In this line of research, Jost et al. [6]

proposed the MoTIF algorithm, which uses an iterative strategy based on generalized eigenvalue decompositions, where the atoms are assumed to be orthogonal to each other and learnt one by one in a greedy way. More recently, the ‘sliding window matching’ (SWM) algorithm

[9] was proposed for learning time-varying atoms by using a correlation-based approach that aims to identify the recurring patterns. Even though some success has been reported with these algorithms, they have several limitations: SWM uses a slow stochastic search inspired by simulated annealing and MoTIF poorly handles correlated atoms, simultaneously activated, or having varying amplitudes; some cases which often occur in practical applications.

A natural way to cast the problem of learning a dictionary of shift-invariant atoms into an optimization problem is a convolutional sparse coding (CSC) approach [10]

. This approach has gained popularity in computer vision 

[11, 12, 13, 14, 15], biomedical imaging [16] and audio signal processing [10, 17], due to its ability to obtain compact representations of the signals and to incorporate the temporal structure of the signals via convolution. In the neuroscience context, Barthélemy et al. [18] used an extension of the K-SVD algorithm using convolutions on EEG data. In a similar spirit, Brockmeier and Príncipe [7] used the matching pursuit algorithm combined with a rather heuristic dictionary update, which is similar to the MoTIF algorithm. In a very recent study, Hitziger et al. [8] proposed the AWL algorithm, which presents a mathematically more principled CSC approach for modeling neural signals. Yet, as opposed to classical CSC approaches, the AWL algorithm imposes additional combinatorial constraints, which limit its scope to certain data that contain spike-like atoms. Also, since these constraints increase the complexity of the optimization problem, the authors had to resort to dataset-specific initializations and many heuristics in their inference procedure.

While the current state-of-the-art CSC methods have a strong potential for modeling neural signals, they might also be limited as they consider an reconstruction error, which corresponds to assuming an additive Gaussian noise distribution. While this assumption could be reasonable for several signal processing tasks, it turns out to be very restrictive for neural signals, which often contain heavy noise bursts and have low signal-to-noise ratio.

In this study, we aim to address the aforementioned concerns and propose a novel probabilistic CSC model called CSC, which is better-suited for neural signals. CSC is based on a family of heavy-tailed distributions called -stable distributions [19] whose rich structure covers a broad range of noise distributions. The heavy-tailed nature of the -stable distributions renders our model robust to impulsive observations. We develop a Monte Carlo expectation maximization (MCEM) algorithm for inference, with a weighted CSC model for the maximization step. We propose efficient optimization strategies that are specifically designed for neural time series. We illustrate the benefits of the proposed approach on both synthetic and real datasets.

2 Preliminaries


For a vector

we denote the norm by . The convolution of two vectors and is denoted by . We denote by the observed signals, the temporal atoms, and the sparse vector of activations. The symbols , , , denote the univariate uniform, exponential, Gaussian, and -stable distributions, respectively.

Convolutional sparse coding: The CSC problem formulation adopted in this work follows the Shift Invariant Sparse Coding (SISC) model from [10]. It is defined as follows:


where denotes one of the observed segments of signals, also referred to as a trials in this paper. We denote by as the length of a trial, and the number of atoms. The aim in this model is to approximate the signals by the convolution of certain atoms and their respective activations, which are sparse. Here, denotes the th atom of the dictionary , and denotes the activation of the th atom in the th trial. We denote by .

The objective function (1) has two terms, an data fitting term that corresponds to assuming an additive Gaussian noise model, and a regularization term that promotes sparsity with an norm. The regularization parameter is called . Two constraints are also imposed. First, we ensure that lies within the unit sphere, which prevents the scale ambiguity between and . Second, a positivity constraint on is imposed to be able to obtain physically meaningful activations and to avoid sign ambiguities between and . This positivity constraint is not present in the original SISC model [10].

Figure 1: (a) PDFs of -stable distributions. (b) Illustration of two trials from the striatal LFP data, which contain severe artifacts. The artifacts are illustrated with dashed rectangles.

-Stable distributions: The -stable distributions have become increasingly popular in modeling signals that might incur large variations [20, 21, 22]

and have a particular importance in statistics since they appear as the limiting distributions in the generalized central limit theorem

[19]. They are characterized by four parameters: , , , and : (i) is the characteristic exponent and determines the tail thickness of the distribution: the distribution will be heavier-tailed as gets smaller. (ii) is the skewness parameter. If , the distribution is symmetric. (iii) is the scale

parameter and measures the spread of the random variable around its mode. Finally, (iv)

is the location parameter.

The probability density function of an

-stable distribution cannot be written in closed-form except for certain special cases; however, the characteristic function can be written as follows:

where for , for , and . As an important special case of the

-stable distributions, we obtain the Gaussian distribution when

and , i.e. . In Fig. 1, we illustrate the (approximately computed) probability density functions (PDF) of the -stable distribution for different values of and . The distribution becomes heavier-tailed as we decrease , whereas the tails vanish quickly when .

The moments of the

-stable distributions can only be defined up to the order , i.e. if and only if

, which implies the distribution has infinite variance when

. Furthermore, despite the fact that the PDFs of -stable distributions do not admit an analytical form, it is straightforward to draw random samples from -stable distributions [23].

3 Alpha-Stable Convolutional Sparse Coding

3.1 The Model

From a probabilistic perspective, the CSC problem can be also formulated as a maximum a-posteriori (MAP) estimation problem on the following probabilistic generative model:


Here, denotes the th element of . We use the same notations for and . It is easy to verify that the MAP estimate for this probabilistic model, i.e. , is identical to the original optimization problem defined in (1).

It has been long known that, due to their light-tailed nature, Gaussian models often fail at handling noisy high amplitude observations or outliers 

[24]. As a result, the ‘vanilla’ CSC model turns out to be highly sensitive to outliers and impulsive noise that frequently occur in electrophysiological recordings, as illustrated in Fig. 1. Possible origins of such artifacts are movement, muscle contractions, ocular blinks or electrode contact losses.

In this study, we aim at developing a probabilistic CSC model that would be capable of modeling challenging electrophysiological signals. We propose an extension of the original CSC model defined in (2) by replacing the light-tailed Gaussian likelihood with heavy-tailed -stable distributions. We define the proposed probabilistic model (CSC) as follows:


where denotes the -stable distribution. While still being able to capture the temporal structure of the observed signals via convolution, the proposed model has a richer structure and would allow large variations and outliers, thanks to the heavy-tailed -stable distributions. Note that the vanilla CSC defined in (2) appears as a special case of CSC, as the -stable distribution coincides with the Gaussian distribution when .

3.2 Maximum A-Posteriori Inference

Given the observed signals , we are interested in the MAP estimates, defined as follows:


As opposed to the Gaussian case, unfortunately, this optimization problem is not amenable to classical optimization tools, since the PDF of the -stable distributions does not admit an analytical expression. As a remedy, we use the product property of the symmetric -stable densities [19, 25] and re-express the CSC model as conditionally Gaussian. It leads to:


where is called the impulse variable that is drawn from a positive -stable distribution (i.e. ), whose PDF is illustrated in Fig. 1. It can be easily shown that both formulations of the

CSC model are identical by marginalizing the joint distribution

over .

The impulsive structure of the CSC model becomes more prominent in this formulation: the variances of the Gaussian observations are modulated by stable random variables with infinite variance, where the impulsiveness depends on the value of . It is also worth noting that when , becomes deterministic and we can again verify that CSC coincides with the vanilla CSC.

The conditionally Gaussian structure of the augmented model has a crucial practical implication: if the impulse variable were to be known, then the MAP estimation problem over and in this model would turn into a ‘weighted’ CSC problem, which is a much easier task compared to the original problem. In order to be able to exploit this property, we propose an expectation-maximization (EM) algorithm, which iteratively maximizes a lower bound of the log-posterior , and algorithmically boils down to computing the following steps in an iterative manner:


where denotes the expectation of a function under the distribution , denotes the iterations, and is a lower bound to and it is tight at the current iterates , .

The E-Step: In the first step of our algorithm, we need to compute the EM lower bound that has the following form:


where denotes equality up to additive constants, denotes the Hadamard (element-wise) product, and the square-root operator is also defined element-wise. Here, are the weights that are defined as follows: . As the variables are expected to be large when cannot explain the observation – typically due to a corruption or a high noise – the weights will accordingly suppress the importance of the particular point . Therefore, the overall approach will be more robust to corrupted data than the Gaussian models where all weights would be deterministic and equal to .

0:  Regularization: , Num. atoms: , Atom length: , Num. iterations: , ,
1:  for  to  do
2:      /* E-step: */
3:     for  to  do
4:        Draw via MCMC (9)
5:     end for
7:      /* M-step: */
8:     for  to  do
9:         = L-BFGS-B on (10)
10:         = L-BFGS-B on the dual of (11)
11:     end for
12:  end for
13:  return  , ,
Algorithm 1 -stable Convolutional Sparse Coding

Unfortunately, the weights

cannot be computed analytically, therefore we need to resort to approximate methods. In this study, we develop a Markov chain Monte Carlo (MCMC) method to approximately compute the weights, where we approximate the intractable expectations with a finite sample average, given as follows:

, where are some samples that are ideally drawn from the posterior distribution . Unfortunately, directly drawing samples from the posterior distribution of is not tractable either, and therefore, we develop a Metropolis-Hastings algorithm [26], that asymptotically generates samples from the target distribution in two steps. In the -th iteration of this algorithm, we first draw a random sample for each and from the prior distribution (cf. (5)), i.e., . We then compute an acceptance probability for each that is defined as follows:


where denotes the iteration number of the MCMC algorithm. Finally, we draw a uniform random number for each and . If , we accept the sample and set ; otherwise we reject the sample and set . This procedure forms a Markov chain that leaves the target distribution invariant, where under mild ergodicity conditions, it can be shown that the finite-sample averages converge to their true values when goes to infinity [27]. More detailed explanation of this procedure is given in the supplementary document.

The M-Step: Given the weights that are estimated during the E-step, the objective of the M-step (7) is to solve a weighted CSC problem, which is much easier when compared to our original problem. This objective function is not jointly convex in and , yet it is convex if one fix either or . Here, similarly to the vanilla CSC approaches [9, 10], we develop a block coordinate descent strategy, where we solve the problem in (7) for either or , by keeping respectively and fixed. We first focus on solving the problem for while keeping fixed, given as follows:


Here, we expressed the convolution of and

as the inner product of the zero-padded activations

, with a Toeplitz matrix , that is constructed from . The matrices are never constructed in practice, and all operations are carried out using convolutions. This problem can be solved by various constrained optimization algorithms. Here, we choose the quasi-Newton L-BFGS-B algorithm [28] with a box constraint: . This approach only requires the simple computation of the gradient of the objective function with respect to (cf. supplementary material). Note that, since each trial is independent from each other, we can solve this problem for each in parallel.

We then solve the problem for the atoms while keeping fixed. This optimization problem turns out to be a constrained weighted least-squares problem. In the non-weighted case, this problem can be solved either in the time domain or in the Fourier domain [10, 11, 12]

. The Fourier transform simplifies the convolutions that appear in least-squares problem, but it also induces several difficulties, such as that the atom

have to be in a finite support , an important issue ignored in the seminal work of [10] and addressed with an ADMM solver in[11, 12]. In the weighted case, it is not clear how to solve this problem in the Fourier domain. We thus perform all the computations in the time domain.

Following the traditional filter identification approach [29], we need to embed the one-dimensional signals into a matrix of delayed signals , where if and elsewhere. Equation (1) then becomes:


Due to the constraint, we must resort to an iterative approach. The options are to use (accelerated) projected gradient methods such as FISTA [30] applied to (11), or to solve a dual problem as done in [10]. The dual is also a smooth constraint problem yet with a simpler positivity box constraint (cf. supplementary material). The dual can therefore be optimized with L-BFGS-B. Using such a quasi-Newton solver turned out to be more efficient than any accelerated first order method in either the primal or the dual (cf. benchmarks in supplementary material)

Our entire EM approach can be summarized in the Algorithm 1. Note that during the alternating minimization, thanks to convexity we can warm start the update and the update using the solution from the previous update. This significantly speeds up the convergence of the L-BFGS-B algorithm, particularly in the later iterations of the overall algorithm.

(a) , .
(b) Time to reach a relative precision of 0.01.
Figure 2: Comparison of state-of-the-art methods with our approach. (a) Convergence plot with the objective function relative to the obtained minimum, as a function of computational time. (b) Time taken to reach a relative precision of , for different settings of and .

4 Experiments

In order to evaluate our approach, we conduct several experiments on both synthetic and real data. First, we show that our proposed optimization scheme for the M-step provides significant improvements in terms of convergence speed over the state-of-the-art CSC methods. Then, we provide empirical evidence that our algorithm is more robust to artifacts and outliers than three competing CSC methods [6, 7, 12]. Finally, we consider LFP data, where we illustrate that our algorithm can reveal interesting properties in electrophysiological signals without supervision, even in the presence of severe artifacts.

Synthetic simulation setup: In our synthetic data experiments, we simulate trials of length by first generating zero mean and unit norm atoms of length

. The activation instants are integers drawn from a uniform distribution in

. The amplitude of the activations are drawn from a uniform distribution in . Atoms are activated only once per trial and are allowed to overlap. The activations are then convolved with the generated atoms and summed up as in (1).

M-step performance: In our first set of synthetic experiments, we illustrate the benefits of our M-step optimization approach over state-of-the-art CSC solvers. We set , and , and use different values for and . To be comparable, we set

and add Gaussian noise to the synthesized signals, where the standard deviation is set to

. In this setting, we have for all , , which reduces the problem to a standard CSC setup. We monitor the convergence of ADMM-based methods by Heide et al. [11] and Wohlberg [12] against our M-step algorithm, using both a single-threaded and a parallel version for the -update. As the problem is non-convex, even if two algorithms start from the same point, they are not guaranteed to reach the same local minimum. Hence, for a fair comparison, we use a multiple restart strategy with averaging across random seeds.

During our experiments we have observed that the ADMM-based methods do not guarantee the feasibility of the iterates. In other words, the norms of the estimated atoms might be greater than during the iterations. To keep the algorithms comparable, when computing the objective value, we project the atoms to the unit ball and scale the activations accordingly. To be strictly comparable, we also imposed a positivity constraint on these algorithms. This is easily done by modifying the soft-thresholding operator to be a rectified linear function. In the benchmarks, all algorithms use a single thread, except “M-step - 4 parallel” which uses 4 threads during the update.

In Fig. 2, we illustrate the convergence behaviors of the different methods. Note that the y-axis is the precision relative to the objective value obtained upon convergence. In other words, each curve is relative to its own local minimum (see supplementary document for details). In the right subplot, we show how long it takes for the algorithms to reach a relative precision of for different settings (cf. supplementary material for more benchmarks). Our method consistently performs better and the difference is even more striking for more challenging setups. This speed improvement on the M-step is crucial for us as this step will be repeatedly executed.

(a) No corruption.
(b) 10% corruption.
(c) 20% corruption
Figure 3: Simulation to compare state-of-the-art methods against CSC.

Robustness to corrupted data: In our second synthetic data experiment, we illustrate the robustness of CSC in the presence of corrupted observations. In order to simulate the likely presence of high amplitude artifacts, one way would be to directly simulate the generative model in (3). However, this would give us an unfair advantage, since CSC is specifically designed for such data. Here, we take an alternative approach, where we corrupt a randomly chosen fraction of the trials ( or ) with strong Gaussian noise of standard deviation , i.e. one order of magnitude higher than in a regular trial. We used a regularization parameter of . In these experiments, by CSC we refer to CSC with , that resembles using only the M-step of our algorithm with deterministic weights for all , . We used a simpler setup where we set , , and . We used atoms, as shown in dashed lines in Fig. 3.

For CSC, we set the number of outer iterations , the number of iterations of the M-step to , and the number of iterations of the MCMC algorithm to . We discard the first samples of the MCMC algorithm as burn-in. To enable a fair comparison, we run the standard CSC algorithm for iterations, i.e. the total number of M-step iterations in CSC. We also compared CSC against competing state-of-art methods previously applied to neural time series: Brockmeier and Príncipe [7] and MoTIF [6]. Starting from multiple random initializations, the estimated atoms with the smallest distance with the true atoms are shown in Fig. 3.

(a) LFP spike data from [8] (b) Estimated atoms
Figure 4: Atoms learnt by CSC on LFP data containing epileptiform spikes with .

In the artifact-free scenario, all algorithms perform equally well, except for MoTIF that suffers from the presence of activations with varying amplitudes. This is because it aligns the data using correlations before performing the eigenvalue decomposition, without taking into account the strength of activations in each trial. The performance of Brockmeier and Príncipe [7] and CSC degrades as the level of corruption increases. On the other hand, CSC is clearly more robust to the increasing level of corruption and recovers reasonable atoms even when 20% of the trials are corrupted.

Results on LFP data In our last set of experiments, we consider real neural data from two different datasets. We first applied CSC on an LFP dataset previously used in [8] and containing epileptiform spikes as shown in Fig. 4(a). The data was recorded in the rat cortex, and is free of artifact. Therefore, we used the standard CSC with our optimization scheme, (i.e. CSC with ). As a standard preprocessing procedure, we applied a high-pass filter at  Hz in order to remove drifts in the signal, and then applied a tapered cosine window to down-weight the samples near the edges. We set , , , , and . The recovered atoms by our algorithm are shown in Fig. 4(b). We can observe that the estimated atoms resemble the spikes in Fig. 4(a). These results show that, without using any heuristics, our approach can recover similar atoms to the ones reported in [8], even though it does not make any assumptions on the shapes of the waveforms, or initializes the atoms with template spikes in order to ease the optimization.

The second dataset is an LFP channel in a rodent striatum from [31]. We segmented the data into trials of length samples, windowed each trial with a tapered cosine function, and detrended the data with a high-pass filter at  Hz. We set , initialized the weights to the inverse of the variance of the trial

. Atoms are in all experiments initialized with Gaussian white noise.

As opposed to the first LFP dataset, this dataset contains strong artifacts, as shown in Fig. 1. In order to be able to illustrate the potential of CSC on this data, we first manually identified and removed the trials that were corrupted by artifacts. In Fig. 5(a), we illustrate the estimated atoms with CSC on the manually-cleaned data. We observe that the estimated atoms correspond to canonical waveforms found in the signal. In particular, the high frequency oscillations around  Hz are modulated in amplitude by the low-frequency oscillation around  Hz, a phenomenon known as cross-frequency coupling (CFC) [32]. We can observe this by computing a comodulogram [33] on the entire signal (Fig. 5(b)). This measures the correlation between the amplitude of the high frequency band and the phase of the low frequency band.

Even though CSC is able to provide these excellent results on the cleaned data set, its performance heavily relies on the manual removal of the artifacts. Finally, we repeated the previous experiment on the full data, without removing the artifacts and compared CSC with CSC, where we set . The results are shown in the middle and the right sub-figures of Fig. 5(a). It can be observed that in the presence of strong artifacts, CSC is not able to recover the atoms anymore. On the contrary, we observe that CSC can still recover atoms as observed in the artifact-free regime. In particular, the cross-frequency coupling phenomenon is still visible.

(a) Atoms learnt by: CSC (clean data), CSC (full data), CSC (full data)
(b) Comodulogram.
Figure 5: (a) Three atoms learnt from a rodent striatal LFP channel, using CSC on cleaned data, and both CSC and CSC on the full data. The atoms capture the cross-frequency coupling of the data (dashed rectangle). (b) Comodulogram presents the cross-frequency coupling intensity computed between pairs of frequency bands on the entire cleaned signal, following [33].

5 Conclusion

This work addresses the present need in the neuroscience community to better capture the complex morphology of brain waves [1, 9]. Our data-driven approach to this problem is a probabilistic formulation of a convolutional sparse coding model [10]. We propose an inference strategy based on Monte Carlo EM to deal efficiently with heavy tailed noise and take into account the polarity of neural activations with a positivity constraint. Our problem formulation allows the use of fast quasi-Newton methods for the M-step which outperform previously proposed state-of-the-art ADMM-based algorithms [11, 12, 34, 35], even when not making use of our parallel implementation. Results on LFP data demonstrate that such algorithms can be robust to the presence of transient artifacts in data and reveal insights on neural time-series without supervision. We will make our code publicly available.

6 Acknowledgement

The work was supported by the French National Research Agency grants ANR-14-NEUC-0002-01 and ANR-16-CE23-0014 (FBIMATRIX), as well as the ERC Starting Grant SLAB ERC-YStG-676943.


  • Cole and Voytek [2017] S. R. Cole and B. Voytek. Brain oscillations and the importance of waveform shape. Trends Cogn. Sci., 2017.
  • Cohen [2014] M. X. Cohen. Analyzing neural time series data: Theory and practice. MIT Press, 2014. ISBN 9780262319560.
  • Jones [2016] S. R. Jones. When brain rhythms aren’t ‘rhythmic’: implication for their mechanisms and meaning. Curr. Opin. Neurobiol., 40:72–80, 2016.
  • Mazaheri and Jensen [2008] A. Mazaheri and O. Jensen. Asymmetric amplitude modulations of brain oscillations generate slow evoked responses. The Journal of Neuroscience, 28(31):7781–7787, 2008.
  • Hari and Puce [2017] R. Hari and A. Puce. MEG-EEG Primer. Oxford University Press, 2017.
  • Jost et al. [2006] P. Jost, P. Vandergheynst, S. Lesage, and R. Gribonval. MoTIF: an efficient algorithm for learning translation invariant dictionaries. In Acoustics, Speech and Signal Processing, ICASSP, volume 5. IEEE, 2006.
  • Brockmeier and Príncipe [2016] A. J. Brockmeier and J. C. Príncipe. Learning recurrent waveforms within EEGs. IEEE Transactions on Biomedical Engineering, 63(1):43–54, 2016.
  • Hitziger et al. [2017] S. Hitziger, M. Clerc, S. Saillet, C. Benar, and T. Papadopoulo. Adaptive Waveform Learning: A Framework for Modeling Variability in Neurophysiological Signals. IEEE Transactions on Signal Processing, 2017.
  • Gips et al. [2017] B. Gips, A. Bahramisharif, E. Lowet, M. Roberts, P. de Weerd, O. Jensen, and J. van der Eerden. Discovering recurring patterns in electrophysiological recordings. J. Neurosci. Methods, 275:66–79, 2017.
  • Grosse et al. [2007] R. Grosse, R. Raina, H. Kwong, and A. Y. Ng. Shift-invariant sparse coding for audio classification. In

    23rd Conference on Uncertainty in Artificial Intelligence

    , UAI’07, pages 149–158. AUAI Press, 2007.
    ISBN 0-9749039-3-0.
  • Heide et al. [2015] F. Heide, W. Heidrich, and G. Wetzstein. Fast and flexible convolutional sparse coding. In

    Computer Vision and Pattern Recognition (CVPR)

    , pages 5135–5143. IEEE, 2015.
  • Wohlberg [2016] B. Wohlberg. Efficient algorithms for convolutional sparse representations. Image Processing, IEEE Transactions on, 25(1):301–315, 2016.
  • Zeiler et al. [2010] M. D. Zeiler, D. Krishnan, G.W. Taylor, and R. Fergus. Deconvolutional networks. In Computer Vision and Pattern Recognition (CVPR), pages 2528–2535. IEEE, 2010.
  • Šorel and Šroubek [2016] M. Šorel and F. Šroubek. Fast convolutional sparse coding using matrix inversion lemma. Digital Signal Processing, 2016.
  • Kavukcuoglu et al. [2010] K. Kavukcuoglu, P. Sermanet, Y-L. Boureau, K. Gregor, M. Mathieu, and Y. Cun. Learning convolutional feature hierarchies for visual recognition. In Advances in Neural Information Processing Systems (NIPS), pages 1090–1098, 2010.
  • Pachitariu et al. [2013] M. Pachitariu, A. M Packer, N. Pettit, H. Dalgleish, M. Hausser, and M. Sahani. Extracting regions of interest from biological images with convolutional sparse block coding. In Advances in Neural Information Processing Systems (NIPS), pages 1745–1753, 2013.
  • Mailhé et al. [2008] B. Mailhé, S. Lesage, R. Gribonval, F. Bimbot, and P. Vandergheynst. Shift-invariant dictionary learning for sparse representations: extending K-SVD. In 16th Eur. Signal Process. Conf., pages 1–5. IEEE, 2008.
  • Barthélemy et al. [2013] Q. Barthélemy, C. Gouy-Pailler, Y. Isaac, A. Souloumiac, A. Larue, and J. I. Mars. Multivariate temporal dictionary learning for EEG. J. Neurosci. Methods, 215(1):19–28, 2013.
  • Samorodnitsky and Taqqu [1994] G. Samorodnitsky and M. S. Taqqu. Stable non-Gaussian random processes: stochastic models with infinite variance, volume 1. CRC press, 1994.
  • Kuruoglu [1999] E. E. Kuruoglu. Signal processing in -stable noise environments: a least Lp-norm approach. PhD thesis, University of Cambridge, 1999.
  • Mandelbrot [2013] B. B. Mandelbrot. Fractals and scaling in finance: Discontinuity, concentration, risk. Selecta volume E. Springer Science & Business Media, 2013.
  • Wang et al. [2016] Y. Wang, Y. Qi, Y. Wang, Z. Lei, X. Zheng, and G. Pan. Delving into -stable distribution in noise suppression for seizure detection from scalp EEG. J. Neural. Eng., 13(5):056009, 2016.
  • Chambers et al. [1976] J. M. Chambers, C. L. Mallows, and B. W. Stuck. A method for simulating stable random variables. Journal of the american statistical association, 71(354):340–344, 1976.
  • Huber [1981] P. J. Huber. Robust Statistics. Wiley, 1981.
  • Godsill and Kuruoglu [1999] S. Godsill and E. Kuruoglu. Bayesian inference for time series with heavy-tailed symmetric -stable noise processes. Proc. Applications of heavy tailed distributions in economics, eng. and stat., 1999.
  • Chib and Greenberg [1995] S. Chib and E. Greenberg. Understanding the Metropolis-Hastings algorithm. The American Statistician, 49(4):327–335, 1995.
  • Liu [2008] J.S. Liu. Monte Carlo strategies in scientific computing. Springer, 2008.
  • Byrd et al. [1995] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16(5):1190–1208, 1995.
  • Moulines et al. [1995] E. Moulines, P. Duhamel, J-F. Cardoso, and S. Mayrargue. Subspace methods for the blind identification of multichannel FIR filters. IEEE Transactions on signal processing, 43(2):516–525, 1995.
  • Beck and Teboulle [2009] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183–202, 2009.
  • Dallérac et al. [2017] G. Dallérac, M. Graupner, J. Knippenberg, R. C. R. Martinez, T. F. Tavares, L. Tallot, N. El Massioui, A. Verschueren, S. Höhn, J.B. Bertolus, et al. Updating temporal expectancy of an aversive event engages striatal plasticity under amygdala control. Nature Communications, 8:13920, 2017.
  • Jensen and Colgin [2007] O. Jensen and L. L. Colgin.

    Cross-frequency coupling between neuronal oscillations.

    Trends in cognitive sciences, 11(7):267–269, 2007.
  • Tort et al. [2010] A. BL. Tort, R. Komorowski, H. Eichenbaum, and N. Kopell. Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol., 104(2):1195–1210, 2010.
  • Wohlberg [2014] B. Wohlberg. Efficient convolutional sparse coding. In Acoustics, Speech and Signal Processing, ICASSP, pages 7173–7177. IEEE, 2014.
  • Bristow et al. [2013] H. Bristow, A. Eriksson, and S. Lucey. Fast convolutional sparse coding. In Computer Vision and Pattern Recognition (CVPR), pages 391–398, 2013.

Appendix A Details of the E-Step

Computing the weights that are required in the M-step requires us to compute the expectation of under the posterior distribution , which is not analytically available.

Monte Carlo methods are numerical techniques that can be used to approximately compute the expectations of the form:


where are some samples drawn from and in our case. However, in our case, sampling directly from is also unfortunately intractable.

MCMC methods generate samples from the target distribution by forming a Markov chain, whose stationary distribution is , so that , where denotes the transition kernel of the Markov chain.

In this study, we develop a Metropolis-Hastings (MH) algorithm, that implicitly forms a transition kernel. The MH algorithm generates samples from a target distribution in two steps. First, it generates a random sample from a proposal distribution , then computes an acceptance probability and draws a uniform random number . If , it accepts the sample and sets ; otherwise it rejects the sample and sets . The acceptance probability is given as follows


where the last equality is obtained by applying the Bayes rule on .

The acceptance probability requires the prior distribution of to be evaluated. Unfortunately, this is intractable in our case since this prior distribution is chosen to be a positive -stable distribution whose PDF does not have an analytical form. As a remedy, we choose the prior distribution of as the proposal distribution, such . This enables us to simplify the acceptance probability. Accordingly, for each , we have the following acceptance probability:


Thanks to the simplification, this probability is tractable and can be easily computed.

Appendix B Details of the M-Step

b.1 Solving for the activations

In the M-step, we optimize (10) to find the activations of each trial independently. To keep the notation simple, we will drop the index for the iteration number of the EM algorithm.

First, this equation can be rewritten by concatenating the Toeplitz matrices for the atoms into a big matrix and the activations for different atoms into a single vector where denotes the transposition operation. Recall that is a zero-padded version of . This leads to a simpler formulation and the objective function :


where is a vector of ones.

The derivative w.r.t. now reads:


In practice, this big matrix is never assembled and all operations are carried out using convolutions. Note also that we do not update the zeros from the padding in . Now that we have the gradient, the activations can be estimated using a efficient quasi-Newton solver such as L-BFGS-B, taking into account the box posititivy constraint .

For each trial, one iteration costs .

b.2 Solving for the atoms

In the M-step, we optimize (11) to find the atoms . As when solving for the activations , we can remove the summation over the atoms by concatenating the delayed matrices into and . This leads to the simpler formulation:


The Lagrangian of this problem is given by:


where are the dual variables. Therefore, the dual problem is:


where , the primal optimal, is given by:


with with . The gradient for the dual variable is given by:


with computed from (20). We can solve this iteratively using again L-BFGS-B taking into account the positivity constraint for all . What we have described so far solves for all the atoms simultaneously. However, it is also possible to estimate the atoms sequentially one at a time using a block coordinate descent (BCD) approach, as in the work of mairal2010online. In each iteration of the BCD algorithm, a residual is computed as given by:


and correspondingly subproblem 17 becomes:


which is solved in the same way as subproblem 17. Now, in the simultaneous case, we construct one linear problem in and one iteration costs . However, in the BCD strategy, we construct linear problems in and one iteration costs only . Interestingly, when the weights are all identical, we can use the fact that for one atom , the matrix is Toeplitz. In this case, we can construct linear problems in only and one iteration costs only .

For the benefit of the reader, we summarize the complexity of the M-step in Table 1. We note and the number of iterations in the L-BFGS-B methods for the activations update and atoms update.

Method Complexity
Solving activations
Solving atoms
Solving atoms (sequential)
Table 1: Complexity analysis of the M-step, where and are the number of iterations in the L-BFGS-B solvers for the activations and atoms updates.

Appendix C Additional Experiments: M-step speed benchmark

c.1 Comparison with state of the art

Here, we compare convergence plots of our algorithm against a number of state-of-art methods. The details of the experimental setup are described in Section 4. Fig. 6 demonstrates on a variety of setups the computational advantage of our quasi-Newton approach to solve the M-step. Note that Fig. 2b is in fact a summary of Fig. 6. Indeed, we can verify that ADMM methods converge quickly to a modest accuracy, but take much longer to converge to a high accuracy.boyd2011distributed

(a) , .
(b) , .
(c) , .
Figure 6: Convergence speed of the relative objective function. The y-axis shows the objective function relative to the obtained minimum for each run:

. Each curve is the geometrical mean over 24 different random initializations.

Next, in Fig. 7, we show more traditional convergence plots. In contrast to Fig. 2 or 6 where the relative objective function is shown, here we plot the absolute value of the objective function. We can now verify that each of the methods have indeed converged to their respective local minimum. Of course, owing to the non-convex nature of the problem, they do not necessarily converge to the same local minimum.

(a) , .
(b) , .
(c) , .
Figure 7: Convergence of the objective function as a function of time. The y-axis shows the absolute objective function . Each curve is the mean over 24 different random initializations.

c.2 Comparison of solver for the activations subproblem

Finally, we compare convergence plots of our algorithm using different solvers for the -update: ISTA, FISTA, and L-BFGS-B. The rationale for choosing a quasi-Newton solver for the -update becomes clear in Fig. 8 as the L-BFGS-B solver turns out to be computationally advantageous on a variety of setups.

(a) , .
(b) , .
(c) , .
Figure 8: Convergence speed of the relative objective function. The y-axis shows the objective function relative to the obtained minimum for each run: . Each curve is the geometrical mean over 24 different random initializations.