Jitter-Adaptive Dictionary Learning - Application to Multi-Trial Neuroelectric Signals
Authors: Sebastian Hitziger, Maureen Clerc, Alexandre Gramfort, Sandrine Saillet, Christian Bénar, Théodore Papadopoulo
ICLR 2013

Abstract: Dictionary Learning has proven to be a powerful tool for many image processing tasks, where atoms are typically defined on small image patches. As a drawback, the dictionary only encodes basic structures. In addition, this approach treats patches of different locations in one single set, which means a loss of information when features are well-aligned across signals. This is the case, for instance, in multi-trial magneto- or electroencephalography (M/EEG). Learning the dictionary on the entire signals could make use of the alignement and reveal higher-level features. In this case, however, small missalignements or phase variations of features would not be compensated for. In this paper, we propose an extension to the common dictionary learning framework to overcome these limitations by allowing atoms to adapt their position across signals. The method is validated on simulated and real neuroelectric data.

arXiv:1301.3611v4  [stat.ML]  24 Jun 2013
Jitter-Adaptive Dictionary Learning - Application to
Multi-Trial Neuroelectric Signals
Sebastian Hitziger
Project-Team Athena
INRIA Sophia Antipolis, France
sebastian.hitziger@inria.fr
Maureen Clerc
Project-Team Athena
INRIA Sophia Antipolis, France
maureen.clerc@inria.fr
Alexandre Gramfort
Institut Mines-Telecom, Telecom ParisTech, CNRS LTCI
alexandre.gramfort@telecom-paristech.fr
Sandrine Saillet
Institut de Neurosciences des Syst` emes
UMR 1106 INSERM
Aix-Marseille Universit´ e
Facult´ e de M´ edecine La Timone
Marseille, France
ssaillet53@gmail.com
Christian B´enar
Institut de Neurosciences des Syst` emes
UMR 1106 INSERM
Aix-Marseille Universit´ e
Facult´ e de M´ edecine La Timone
Marseille, France
christian.benar@univmed.fr
Th´eodore Papadopoulo
Project-Team Athena
INRIA Sophia Antipolis, France
theodore.papadopoulo@inria.fr
Abstract
Dictionary Learning has proven to be a powerful tool for many image processing
tasks, where atoms are typically deﬁned on small image patch es. As a drawback,
the dictionary only encodes basic structures. In addition, this approach treats
patches of different locations in one single set, which mean s a loss of informa-
tion when features are well-aligned across signals. This is the case, for instance,
in multi-trial magneto- or electroencephalography (M/EEG). Learning the dictio-
nary on the entire signals could make use of the alignment andreveal higher-level
features. In this case, however, small misalignments or pha se variations of fea-
tures would not be compensated for. In this paper, we propose an extension to the
common dictionary learning framework to overcome these limitations by allowing
atoms to adapt their position across signals. The method is validated on simulated
and real neuroelectric data.
1 Introduction
The analysis of electromagnetic signals induced by brain activity requires sophisticated tools capable
of efﬁciently treating redundant, multivariate datasets. Redundancy originates for example from the
spatial dimension as in multi-channel magneto- or electroe ncephalography (M/EEG). It can also
1
result from repetitions of the same phenomenon, for instance when a stimulus is presented multiple
times to a subject (typically between 50 and 500 times) and th e neuronal response is recorded; we
will refer to this case as multi-trial analysis.
In the case of multi-channel data, principal component anal ysis (PCA) (Pearson, 1901) and inde-
pendent component analysis (ICA) (Comon, 1994) have been su ccessfully used to decompose the
data into a few waveforms, providing insight into the underl ying neuronal activity and allowing to
enhance the often poor signal-to-noise-ratio (Lagerlund et al., 1997; Makeig et al., 1996). They
both use the fact that the data of all channels are recorded sy nchronically such that features appear
well-aligned and phase-locked.
This condition typically does not hold for multi-trial anal ysis though. In (Woody, 1967) a method
is provided to compensate for temporal jitter across signal s, but it assumes a single underlying
waveform. Matching pursuit (MP) algorithms (Mallat & Zhang , 1993; Durka & Blinowska, 1995)
in turn allow to extract several different features and have recently been adapted to deal with multi-
channel (Durka et al., 2005; Gribonval, 2003) as well as multi-trial (B´ enaret al., 2009) M/EEG data
by compensating for different types of variability. Howeve r, these methods only allow to extract
waveforms that have previously been deﬁned in a dictionary.
In the ﬁeld of image processing, learning dictionaries directly from the data has shown to yield state-
of-the-art results in several applications (Elad & Aharon,2006; Mairal et al., 2008). Typically these
dictionaries are learned on small patches and represent thebasic structures of the images, e.g. edges.
When using this technique for neuroelectric multi-trial an alysis though, the framework should be
carefully adapted to the properties of the data: (i) the wave forms of interest occur approximately
at the same time across trials; (ii) however, they may have slightly different time delays (of a small
fraction of the signal lenth) or phase variations; (iii) thenoise-level is high, partially due to neuronal
background activity which is non-Gaussian; (iv) data is often limited to a few hundred signals. The
latter two properties make the analysis a difﬁcult problem; it is therefore necessary to incorporate
all prior information about the data, in particular (i) and ( ii). We note that similar properties can be
found in many other signal processing applications, such as in other bioelectric or biomagnetic data
(e.g. ECG, EMG).
We suggest that atoms should be learned on the entirety of the signals to provide global high-level
features. The common dictionary learning formulation as a m atrix factorization problem, however,
cannot compensate for the time delays (ii). Existing extens ions known as convolutional or shift-
invariant sparse coding (SISC) (Smith & Lewicki, 2005; Blumensath & Davies, 2006; Grosse et al.,
2007; Ekanadham et al., 2011) learn atoms that are typically smaller than the signal and can occur at
arbitrary and possibly multiple positions per signal. This framework is very general for our purpose
and does not make use of property (i), the approximate alignm ent of waveforms. In addition, the
SISC framework leads to a highly complex algorithm since all shifts of all atoms are considered
in the optimization. The sparse coding step is therefore oft en handled by heuristic preselection of
the active atoms, as described in (Blumensath & Davies, 2006 ). But the update of the dictionary
elements is also a difﬁcult task, as it involves solving a convolutional problem (Grosse et al., 2007).
In this paper, we present a novel dictionary learning framew ork that is designed speciﬁcally to
compensate for small temporal jitter of waveforms across tr ials, leading to the name jitter-adaptive
dictionary learning (JADL). In contrast to SISC, atoms lear ned by JADL are deﬁned on the entire
signal domain and are allowed to shift only up to a small fract ion of the signal length. The most
important difference, however, is a constraint for atoms to occur at most once (i.e. in one position)
per signal, see section 3.2.1. On the one hand, this constrai nt is reasonable since we do not want
to encode a waveform with multiple slightly shifted copies o f one atom. On the other hand, it
signiﬁcantly reduces complexity compared to SISC.
An important difference to previous dictionary learning fr ameworks is the size of the dictionary;
while for image processing dictionaries are often overcomplete, JADL aims at learning only a small
number of atoms. This is not only desired for easy interpretability, but also because of the difﬁculties
introduced by (iii) and (iv), that make it infeasible to learn a large number of atoms. The “unrolled”
version of the dictionary, i.e. the set of all allowed shifts of all atoms may still be large; it therefore
makes sense to speak of sparse solutions with respect to this unrolled dictionary. However, JADL
enforces sparsity to a major part by the explicit constraint mentioned; sparse regularization only
plays a minor role.
2
We begin by brieﬂy stating the common dictionary learning problem after which we will present the
theory and implementation details of JADL. Finally, JADL isevaluated on synthetic and experimen-
tal data.
2 Dictionary learning: prior art
A dictionary consists of a matrix D ∈RN×K that contains for columns the N-dimensional column
vectors {di}K
i=1, its atoms. For a set of signals {xj ∈RN }M
j=1 the problem of ﬁnding a sparse code
aj ∈RK for each xj can be formulated as the following minimization problem:
aj = argmin
aj ∈RK
1
2 ∥xj −Daj ∥2
2 + λ∥aj ∥1 , (1)
where ∥·∥1 denotes the l1-norm and λ> 0 is a regularization parameter. This problem is known as
Lasso (Tibshirani, 1996) and can be solved efﬁciently with algorithms such as least angle regression
(LARS) (Efron et al. , 2004) or the proximal methods ISTA (Combettes & Wajs, 2005) and its
accelerated version FISTA (Beck & Teboulle, 2009).
The case where D is not known beforehand but shall be estimated given the sign als {xj }, leads to
the dictionary learning problem. It is a minimization problem over both the dictionary and the sparse
code, which reads:
min
D, aj
1
2
M∑
j=1
(
∥xj −Daj ∥2
2 + λ∥aj ∥1
)
,
s.t. ∥di∥2 = 1 , i = 1 ,...,K,
(2)
where the latter constraint prevents atoms from growing arbitrarily large.
Most algorithms tackle this non-convex problem iteratively by alternating between the convex sub-
problems: (i) the sparse coding (with D ﬁxed) and (ii) the dictionary update ( {aj }M
j=1 ﬁxed). The
ﬁrst such algorithm was provided in the pioneer work on dictionary learning in (Olshausen & Field,
1997). Many alternatives have been proposed, such as the met hod of optimal directions (MOD) in
(Engan et al., 1999), K-SVD (Aharon et al., 2006), or more recently an online version (Mairalet al.,
2010) to handle large datasets.
3 Jitter-adaptive dictionary learning
This section introduces the main contribution of this paper: a novel technique designed to overcome
the limitations of purely linear signal decomposition meth ods such as PCA, ICA, and dictionary
learning.
We suppose that atoms present in a signal can suffer from unknown time delays, which we will refer
to as jitter. This type of variability addresses the issue of varying lat encies of transient features as
well as oscillations with different phases across trials.
This issue is very important for interpretation of M/EEG data, to answer fundamental questions such
as the link between evoked responses and oscillatory activity (Hanslmayer et al., 2007; Mazaheri &
Jensen, 2008), the correlation between single-trial activ ity and behavioral data (Jung et al., 2001),
the variability of evoked potentials such as the P300 (Holm et al., 2006). Cross-trial variability
is also a precious source of information for simultaneous EE G-fMRI interpretation (B´ enaret al.,
2007).
We therefore provide a dictionary learning framework in whi ch atoms may adapt their position
across trials. This framework can be generalized to address other kinds of variability; the shift
operator introduced below can simply be replaced by the desi red operators. The entire framework
remains the same, only the formula for the update of the dictionary needs to be adapted as described
in Section 3.2.2.
3
3.1 Model and problem statement
Our model is based on the hypothesis that the set of signals of interest {xj }M
j=1 can be generated by
a dictionary D = {di}K
i=1 with few atoms K in the following way: Given a set of shift operations
∆ 1 , for every jthere exist coefﬁcients aij ∈R and shift operators δij ∈∆ , such that
xj =
K∑
i=1
aijδij (di) . (3)
We assume that ∆ contains only small shifts relative to the size of the time wi ndow. Now, we can
formulate the jitter-adaptive dictionary learning problem
min
di,a ij ,δ ij
M∑
j=1

1
2




xj −
K∑
i=1
aij δij(di)





2
2
+ λ∥aj ∥1

 ,
s.t. ∥di∥2 = 1 , δ ij ∈∆ , i = 1 ,...,K, j = 1 ,...,M.
(4)
Note that for ∆ = {I}this reduces to Eq. (2), the problem is thus also non-convex a nd we solve it
similarly using alternate minimizations.
3.2 Implementation
The algorithm that we propose for solving Eq. (4) is based on a n implementation in (Mairal et al.,
2010) for common dictionary learning, which iteratively al ternates between (i) sparse coding and
(ii) dictionary update. We adapt the algorithm least angle regression (LARS) used for solving (i) to
ﬁnd not only the coefﬁcients {aij}but also the latencies {δij}. For (ii), block coordinate descent is
used. The entire procedure is summarized in Algorithm 1 (not ation and details are explained in the
following).
Algorithm 1 Jitter-Adaptive Dictionary Learning
Require: signals {x1,x2,..., xM },shift operators ∆ ,K ∈N,λ ∈R,
1: Initialize D = {d1,d2,..., dK }
2: repeat
3: Set up “unrolled” dictionary DS from D
4: Sparse coding (solve using modiﬁed LARS)
5: for j = 1 to M do
6:
aS
j ←argmin 1
2

xj −DS aS
j

2
2 + λ

aS
j


1 , s.t.


aS,i
j



0
≤1
7: end for
8: Convert {aS
j }to {aij},{δij}
9: Dictionary update (solve using block coordinate descent)
10:
D ←argmin
{di}K
i=1
M∑
j=1
1
2




xj −
K∑
i=1
aijδij (di)





2
2
, s.t. ∥di∥= 1
11: until convergence
1There are different possibilities to deﬁne these shifts. As our entire framework is valid even for arbitrary
linear transforms, we do not specify the choice at this point . While circular shifts, i.e. δn(d) = dn for n ∈ N
and dn[i] := d[(i − n)modN], result in a slightly simpler formulation of the dictionary update and may have
minor computational advantages, they can introduce unwanted boundary effects. Our implementation actually
uses atoms deﬁned on a slightly larger domain (N + S − 1 sample points) than the signals, this way avoiding
circular shifts. For the sake of simplicity, however, we her e assume atoms to be deﬁned on the same domain
as the signals. Although the right way to handle boundary eff ects can be an important question, it is out of
the scope of this paper to discuss this issue in detail. In our experiments, we found the impact of the concrete
deﬁnition of the δ to be small.
4
3.2.1 Sparse coding
When D is ﬁxed, the minimization Eq. (4) can be solved independently for each signal xj ,
min
aij ,δ ij
1
2




xj −
K∑
i=1
aijδij(di)





2
2
+ λ∥aj ∥1 . (5)
We now rewrite this problem into a form similar to the Lasso, w hich allows us to solve it using
a modiﬁcation of LARS. Let us ﬁrst deﬁne an “unrolled” versio n of the dictionary containing all
possible shifts of all its atoms; this is given by DS = {δ(d) : d ∈ D,δ ∈ ∆ }, a matrix of
dimension N ×KS, where S = |∆ |is the number of allowed shifts. The decomposition in Eq. (5)
can now be rewritten as a linear combination over the unrolled dictionary
K∑
i=1
aijδij(di) = DS aS
j ,
where aS
j ∈RKS denotes the corresponding coefﬁcient vector. This vector i s extremely sparse;
in fact, each subvector aS,i
j of aS
j that contains the coefﬁcients corresponding to the shifts o f atom
di shall maximally have one non-zero entry. If such a non-zero e ntry exists, its position indicates
which shift was used for atom di. Now Eq. (5) can be rewritten as
aS
j ← argmin 1
2

xj −DSaS
j

2
2 + λ

aS
j


1 , (6)
s.t.


aS,i
j



0
≤1, i = 1 ,...,K. (7)
Clearly, Eq. (6) is the Lasso, but the constraint (7) leads to a non-convex problem. Therefore the
modiﬁcation of the LARS that we propose below only guarantees convergence to a local minimum.
The LARS algorithm (Efron et al., 2004) follows a stepwise procedure; in each step the coefﬁc ient
of one atom is selected to change from “inactive” to “active” (i.e. it changes from zero to non-zero)
or vice versa. In order to enforce the constraint (7), we make the following modiﬁcation. When
a coefﬁcient is selected for activation, we determine the in dex i such that this coefﬁcient lies in
aS,i
j . We then block all other coefﬁcients contained in the subvec tor aS,i
j such that they cannot get
activated in a later step. In the same manner, we unblock all entries of aS,i
j when its active coefﬁcient
is deactivated.
As mentioned in the introduction, the constraint (7), which is the main difference to SISC, helps to
reduce the complexity of the optimization. In fact, each tim e an atom is activated, all its translates
are blocked and do not need to be considered in the following steps, which facilitates the calculation.
In addition, maximally K steps have to be performed (given that no atom is deactivated in a later
step, which we observed to occur rarely). As suggested for ex ample in (Grosse et al., 2007) the
inital correlations of the shifted atoms with the signal can be computed using fast convolution via
FFT, which speeds up computation in the case of a large number of shifts S.
3.2.2 Dictionary update
For the dictionary update, its unrolled version cannot be used as this would result in updating differ-
ent shifts of the same atom in different ways. Instead, the shifts have to be explicitly included in the
update process. We use block coordinate descent to iterativ ely solve the constrained minimization
problem
dk = argmin
dk
M∑
j=1
1
2




xj −
K∑
i=1
aijδij (di)





2
2
, s.t. ∥dk∥2 = 1
for each atom dk. This can be solved in two steps, the solution of the unconstr ained problem by
differentiation followed by normalization. This is summarized by
˜dk =
M∑
j=1
akj δ−1
kj

xj −
∑
i̸=k
aijδij (di)

 , (8)
dk =
˜dk
∥˜dk∥2
.
5
Original
Atom 1 Atom 2 Atom 3
PCADLJADL
1 2 3
atom
ρ
1 2 3
atom
¯E
 
 
PCA
DL
JADL
Figure 1: Original dictionary and reconstruc-
tions with PCA, DL, and JADL, respectively;
row 5: similarity ρand average energy across
signals ¯E for each atom, the dashed line
marks the value 1 in both plots.
Signal 1
cleanspuriousnoisyPCADL
0 2 4
JADL
time [s]
Signal 2
0 2 4
time [s]
Signal 3
0 2 4
time [s]
Figure 2: Clean, noisy and denoised signals,
row 1: clean signals; row 2: signals plus spu-
rious events; row 3: signals above plus white
noise, row 4-6: denoised signals using PCA,
DL, and JADL, respectively.
with δ−1 the opposite shift of δ. As in (Mairal et al., 2010), we found that one update loop through
all of the atoms was enough to ensure fast convergence of Algorithm 1.
The only difference of this update compared to common dictio nary learning are the shift opera-
tors δij. In contrast to the sparse coding step, the jitter adaptivit y therefore does not increase the
complexity for this step. When high efﬁciency of the algorit hm is needed, this fact can be used by
employing mini-batch or online techniques (Mairal et al., 2010), which increase the frequency of
dictionary update steps with respect to sparse coding steps.
We note that in Eq. (8) we assumed the shifts to be circular; be ing orthogonal transforms, i.e.
δδt = I, they provide for a simple update formula. In the case of non- circular shifts or other
linear operators, the inverse δ−1
kj in the update Eq. (8) needs to be replaced by the adjoint δt
kj . In
addition, the rescaling function ψ= ( ∑ M
j=1 a2
kj δkj δt
kj )−1 has to be applied to the update term.
3.2.3 Hyperparameters and initial dictionary
As mentioned before, the optimal number K of atoms will typically be lower than for dictionary
learning in image processing frameworks; this is due to the high redundancy of the electromagnetic
brain signals as well as the adaptive property of the atoms. W hen oscillatory activity is sought, ∆
should contain shifts of up to the largest period expected to be found. The choice of λplays a less
important role than in common dictionary learning; most of t he sparsity is induced directly by the
constraint (7).
The choice of the initial dictionary can be crucial for the ou tcome, due to the non-convex nature of
the problem. Especially when the initial dictionary alread y produces a small sparse coding error,
the algorithm may converge to a local minimum that is very clo se to the initialization. While this
property allows us to incorporate a priori knowledge by initializing with a predeﬁned dictionary, this
also provides the risk of preventing the algorithm from learning interesting features if they are “far”
from the initial dictionary.
When no dictionary is given a priori, initializations often found in the literature include random
values as in (Olshausen & Field, 1997) as well as random examp les (Aharon et al., 2006) or linear
combinations (e.g. PCA) of the input signals. In order to kee p the bias towards the initialization as
small as possible, we favor random initializations independent from the signals.
6
4 Experiments
Jitter-adaptive dictionary learning (JADL) is next evaluated on synthetic and on experimental data.
Its performance is compared to results obtained with principal component analysis (PCA)(Pearson,
1901) and dictionary learning (DL); for the latter we used th e open-source software package
SPAMS2 whose implementation is described in (Mairal et al., 2010).
4.1 Synthetic data
Generating the dictionary: First, a dic-
K=3 K=4 K=5 K=6 K=8 K=10 K=12
¯ρ PCA 0.522 0.522 0.522 0.522 0.522 0.522 0.522
¯ρ DL 0.563 0.566 0.598 0.615 0.589 0.595 0.581
¯ρ JADL 0.955 0.954 0.946 0.911 0.931 0.881 0.801
λ DL 0.001 0.005 0.001 0.001 0.005 0.050 0.2
λ JADL 0.001 0.005 0.005 0.01 0.005 0.1 0.4
Table 1: First three rows: average similarity ¯ρof origi-
nal and reconstructed atoms; last two rows: parameter
λused for reconstruction.
tionary D with K = 3 normalized atoms
was deﬁned, as shown in Figure 1 (ﬁrst
row). The ﬁrst atom is a delta-shaped spike
and the other two are oscillatory features.
The length of the time window was chosen
as 4 seconds and the sampling rate as 128
Hz, hence N = 512 .
Generating the signals: From the dictio-
nary, M = 200 signals were generated ac-
cording to the model Eq. (3). The coefﬁ-
cients aij and shifts δij were drawn independently from Gaussian distributions withmean µ= 1 and
standard deviation σ = 0 .3 for the coefﬁcients and µ= 0 ,σ = 0 .2 s for the shifts. Three examples
are shown in the ﬁrst row of Figure 2. These signals were then corrupted by two types of noise: (i) to
every signal between 0 and 3 oscillatory events were added, their amplitudes, frequencies, temporal
supports and positions in time were drawn randomly (Fig. 2, r ow 2); (ii) then white Gaussian noise
with a resulting average SNR (energy of clean signals/energy of noise) of 0.790 was added (row 3).
Reconstructions: Performance of the three methods PCA, DL, and JADL given the n oisy signals
(Fig. 2, row 3) was measured on both their ability of recoveri ng the original dictionary and the
denoising of the signals. We performed reconstruction for d ifferent dictionary sizes K. For DL
and JADL we chose the λthat gave the best results; we noticed relatively small sens itivity of the
methods to the choice of λ, especially for small values of K.
Recovering the dictionary: For PCA, the dictionary atoms were deﬁned as the ﬁrst K principal
components. For JADL, the set ∆ was deﬁned to contain all time shifts of maximally±0.6 seconds,
resulting in S = 128 Hz ·1.2 s ≈154 allowed shifts. For each reconstructed dictionary, the thr ee
atoms that showed the highest similarity ρto the original atoms were used to calculate the average
similarity ¯ρ. ρwas deﬁned as the maximal correlation of the reconstructed a tom and all shifts of
maximally ±0.6 seconds of the original atom. The values ¯ρare shown for PCA, DL, and JADL in
Table 1 for each K (for K >12 we observed decreasing similarity values for all three meth ods).
The similarity for JADL is signiﬁcantly higher than for PCA a nd DL; its optimal value is obtained
for K = 3 atoms, but it remains high for larger K, showing its robustness to overestimation of
original atoms. Note that dictionaries obtained by PCA for d ifferent values of K always have their
ﬁrst atoms in common, hence the constant similarity values.
For each method, the K giv-
K=1 K=2 K=3 K=4 K=5 K=6 K=8 K=10 K=12
ǫ PCA 0.871 0.750 0.638 0.539 0.522 0.508 0.502 0.537 0.570
ǫ DL 0.869 0.747 0.635 0.535 0.515 0.498 0.487 0.505 0.539
ǫ JADL 0.505 0.283 0.214 0.230 0.277 0.284 0.317 0.325 0.330
λ DL 0.05 0.05 0.05 0.05 0.05 0.05 0.1 0.1 0.1
λ JADL 0.05 0.2 0.1 0.2 0.2 0.3 0.4 0.5 0.5
Table 2: Relativel2-error ǫproduced by each method for different
K; the last two rows show the values ofλused for DL and JADL.
ing the highest similarity was
determined and the three atoms
with largest ρ-values of the
corresponding dictionaries are
shown in Figure 1 (row 2 - 4).
The atoms found by PCA and
DL contain only mixtures of
the oscillatory atoms, the spike
does not appear at all. JADL
succeeds in separating the three
atoms; their shapes and average energy across signals ¯Eare very close to the originals, as shown in
the bar plots.
2http://spams-devel.gforge.inria.fr/
7
−4 −2 0 2 4
−0.2
−0.1
0
0.1
Epoch 21
time [s]
−4 −2 0 2 4
Epoch 101
time [s]
−4 −2 0 2 4
Epoch 161
time [s]
−4 −2 0 2 4
Evoked potential
time [s]
Figure 3: Different epochs of the local ﬁeld potential (LFP) showing spikes with decreasing energy;
the evoked potential (last) is the average over all epochs.
−0.2
−0.1
0
0.1
0.2
Atom 1
PCA 
Atom 2 Atom 3 Atom 4 Atom 5
0
0.5
1
Av. energy
−0.2
−0.1
0
0.1
0.2
DL 
0
0.5
1
−2 0 2
−0.2
−0.1
0
0.1
0.2
time [s]
JADL 
−2 0 2
time [s]
−2 0 2
time [s]
−2 0 2
time [s]
−2 0 2
time [s]
1 2 3 4 5
atoms
0
0.5
1
Figure 4: Dictionaries learned on LFP epochs using PCA, DL, and JADL.
Denoising the signals: For PCA, denoising was performed by setting the coefﬁcients of all but the
ﬁrst K components to zero. For DL and JADL, the noisy signals were en coded over the learned
dictionaries according to Eq. (1) and (5), respectively. Table 2 shows the average relativel2-errors ǫ
of all denoised signals with respect to the original ones. For each method, three denoised signals for
optimal K are shown in the last three rows of Figure 2. Despite the large r dictionary size (K = 8 )
in the case of PCA and DL, the spike could not be reconstructed due to its jitter across signals. In
addition, the locations of the oscillatory events in the signals denoised with JADL correspond better
to their true locations than it is the case for PCA and DL, especially in the last signal shown. Finally,
the white noise appears more reduced for JADL, which is due toits smaller dictionary size. All three
methods succeeded in removing the spurious events from the signals.
4.2 Real data
Local ﬁeld potentials (LFP) were recorded with an intra-cranial electrode in a Wistar-Han rat, in an
animal model of epilepsy. Bicuculline (a blocker of inhibit ion) was injected in the cortex in order
to elicit epileptic-like discharges. Discharges were selected visually on the LFP traces, and the data
was segmented in 169 epochs centered around the spikes. To simplify the analysis, the epochs were
scaled with a constant factor such that the maximal l2-energy among epochs was equal to 1.
Three examples of epochs are shown in Figure 3, as well as the e voked potential, measured as the
average over the epochs. The only structure visible to the ey e is a spike with changing shape and
decreasing energy across epochs.
Learning the dictionary: Five normalized atoms were learned on the data with PCA, DL an d
JADL; see Figure 4. They were ordered in descending order of their average energy across epochs.
All three methods produce for their ﬁrst atom a spike that resembles the evoked potential (Figure 3);
also, all methods reveal an oscillatory artifact around 1.2 Hz which is not visible in the evoked
potential. However, while the oscillations in the PCA and DLdictionaries are encoded in two atoms
8
epochs
atoms
Absolute coefficient values
 
 
50 100 150
1
2
3
4
5
0
0.2
0.4
0.6
0.8
(a) Absolute coefﬁcients of atoms
learned with JADL.
Latency distributions
−0.3 −0.2 −0.1 0 0.1 0.2 0.3
latency [s]
(b) Distribution of latencies for
each atom.
0 100 200 3000
10
20
30
40
50
Calculation time
time [s]
S
 
 
K=3
K=5
K=7
K=10
K=15
(c) Computation times for differ-
ent values of S and K.
Figure 5: Code visualization and computation times.
(4 and 5) that differ mostly in phase, they are concentrated in atom 5 for JADL. Additionally, in the
case of JADL the oscillations are almost completely separated from the spike, only a small remainder
is still visible. This shows the need of PCA and DL for severalatoms to compensate for phase shifts
while JADL is able to associate oscillations with differentphases; moreover JADL makes use of the
varying phases to clearly separate transient from oscillat ory events. In addition, we can observe a
smoothing effect in the case of PCA and DL: the oscillations look very much like sinusoids whereas
atom 5 in JADL shows a less smooth, rather spiky structure.
Interpreting the code: We visualized the coefﬁcients and the shifts obtained by dec omposing all
the epochs over the dictionary learned with JADL, see Figure 5. Interestingly, each of the ﬁrst three
spikes in the dictionary is active during a contiguous set of epochs during which the other atoms
are hardly used. This allows to segment the epochs into three sets of different spike shapes. The
fourth atom is only active during a few epochs where the domin ant atom is changing from the ﬁrst
to the second spike. The oscillatory artifact (atom 5) in contrast shows very low but constant activity
across all epochs.
The latency distributions (Fig. 5(b)) can provide further i nsight into the data. The highly peaked
distribution for the ﬁrst atoms gives evidence of the accura te alignment of the spikes. The last
atom shows shifts in the whole range from −0.4 to 0.4 seconds, indicating that the phases of the
oscillations were uniformly distributed across epochs.
Computation times: The dictionary and the code above were obtained for 189 itera tions of JADL
after which the algorithm converged. The time of computatio n on a laptop (Intel Core CPU, 2.4
GHz) was 4.3 seconds. As JADL is designed for datasets of simi lar complexity to the one investi-
gated here, computation time should not be an issue for ofﬂine analysis. However, applications such
as M/EEG based brain computer interfaces (BCI) may require computations in real time. We remark
that our implementation has not yet been optimized and could be speeded up signiﬁcantly by paral-
lelization or online techniques mentioned in section 3.2.2 . If training data is available beforehand,
the dictionary may also be calculated ofﬂine in advance; the sparse encoding of new data over this
dictionary then only takes up a small fraction of the training time and can be performed online.
Figure 5(c) illustrates the effect of changing the number of atoms K and shifts S on computation
time t; for each calculation 200 iterations were performed. We can see that tis linearly correlated
with S but increases over-linearly with K: while both, S and K affect the size of the unrolled dic-
tionary DS , an increase of S is handled more efﬁciently by using calculation advantages described
in section 3.2.1; e.g., the non-smooth behavior of the curves at S = 40 results from the fact that for
S >40 an FFT-based convolution is used. In addition, the dictionary update step only uses the com-
pact dictionary D whose size does not increase with S. If the additional factor in computation time
due to S is still not acceptable, ∆ may be subsampled by introducing a minimal distance between
shifts; the tradeoff is a less exact description of the atoms latencies.
9
5 Conclusion
In this paper, a novel method for the analysis and processing of multi-trial neuroelectric signals
was introduced. The method was derived by extending the dictionary learning framework, allowing
atoms to adapt to jitter across trials; hence the name jitter-adaptive dictionary learning (JADL). It was
shown how the resulting non-convex minimization problem ca n be tackled by modifying existing
algorithms used in common dictionary learning.
The method was validated on synthetic and experimental data , both containing variability in laten-
cies of transient and in the phases of oscillatory events. Th e results obtained showed to be superior
to those of common dictionary learning and PCA, both in recov ering the underlying dictionary and
in denoising the signals. The evaluation furthermore demonstrated the usefulness of JADL as a data
exploration tool, capable of extracting global, high-leve l features of the signals and giving insight
into their distributions.
Acknowledgments
This work was supported by the doctoral grant of the region Pr ovence-Alpes-Cˆ ote d’Azur and the
ANR grants CoAdapt (09-EMER-002-01) and MultiModel (2010 BLAN 0309 04).
References
Aharon, M., Elad, M., & Bruckstein, A. 2006. K-SVD: An Algori thm for Designing Overcom-
plete Dictionaries for Sparse Representation. IEEE Transactions on Signal Processing, 54(11),
4311–4322.
Beck, A., & Teboulle, M. 2009. A fast iterative shrinkage-th resholding algorithm for linear inverse
problems. SIAM Journal on Imaging Sciences, 2(1), 183–202.
B´ enar, C.G., Sch¨ on, D., Grimault, S., Nazarian, B., Burle, B., Roth, M., Badier, J.M., Marquis,
P., Liegeois-Chauvel, C., & Anton, J.L. 2007. Single-trial analysis of oddball event-related
potentials in simultaneous EEG-fMRI. Human Brain Mapping, 28, 602613.
B´ enar, C.G., Papadopoulo, T., Torr´ esani, B., & Clerc, M. 2009. Consensus matching pursuit for
multi-trial EEG signals. Journal of neuroscience methods, 180(1), 161–170.
Blumensath, Thomas, & Davies, Mike. 2006. Sparse and shift- invariant representations of music.
Audio, Speech, and Language Processing, IEEE Transactionson, 14(1), 50–57.
Combettes, P. L., & Wajs, V . R. 2005. Signal recovery by proxi mal forward-backward splitting.
Multiscale Modeling and Simulation, 4(4), 1168–1200.
Comon, P. 1994. Independent component analysis, a new conce pt? Signal processing, 36(3), 287–
314.
Durka, PJ, & Blinowska, KJ. 1995. Analysis of EEG transients by means of matching pursuit.
Annals of biomedical engineering, 23(5), 608–611.
Durka, P.J., Matysiak, A., Montes, E.M., Sosa, P.V ., & Blinowska, K.J. 2005. Multichannel match-
ing pursuit and EEG inverse solutions. Journal of neuroscience methods, 148(1), 49–59.
Efron, B., Hastie, T., Johnstone, I., & Tibshirani, R. 2004. Least angle regression. The Annals of
statistics, 32(2), 407–499.
Ekanadham, Chaitanya, Tranchina, Daniel, & Simoncelli, Ee ro P. 2011. Sparse decomposition of
transformation-invariant signals with continuous basis pursuit. Pages 4060–4063 of: Acoustics,
Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on. IEEE.
Elad, M., & Aharon, M. 2006. Image denoising via sparse and re dundant representations over
learned dictionaries. Image Processing, IEEE Transactions on, 15(12), 3736–3745.
Engan, K., Aase, S.O., & Hakon Husoy, J. 1999. Method of optim al directions for frame design.
Pages 2443–2446 of: Acoustics, Speech, and Signal Processi ng, 1999. Proceedings., 1999
IEEE International Conference on, vol. 5. IEEE.
Gribonval, R. 2003. Piecewise linear source separation. Pages 297–310 of: Optical Science and
Technology, SPIE’s 48th Annual Meeting. International Society for Optics and Photonics.
10
Grosse, R, Raina, R, Kwong, H, & Ng, AY . 2007. Shift-invariant sparse coding for audio classiﬁca-
tion. In: Proceedings of the Twenty-third Conference on Uncertai nty in Artiﬁcial Intelligence
(UAI07).
Hanslmayer, S., Klimesch, W., Sauseng, P., Gruber, W., Dopplemayr, M., Freunberger, R., Pecher-
storfer, T., & Birbaumer, N. 2007. Alpha phase reset contrib utes to the generation of ERPs.
Cerebral Cortex, 17, 1–8.
Holm, A., Ranta-aho, P., Sallinen, M., Karjalainen, P., & Ml ler, K. 2006. Relationship of P300
single-trial responses with reaction time and preceding stimulus sequence. Int. J. Psychophys-
iol.
Jung, T.P., Makeig, S., Westerﬁeld, M., Townsend, J., Courc hesne, E., & Sejnowski, T.J. 2001.
Analysis and Visualization of Single-Trial Event-Related Potentials. Human Brain Mapping,
14, 166–185.
Lagerlund, T.D., Sharbrough, F.W., & Busacker, N.E. 1997. S patial ﬁltering of multichannel elec-
troencephalographic recordings through principal compon ent analysis by singular value de-
composition. Journal of Clinical Neurophysiology, 14(1), 73–82.
Mairal, J., Bach, F., Ponce, J., Sapiro, G., & Zisserman, A. 2 008. Supervised dictionary learning.
arXiv preprint arXiv:0809.3083.
Mairal, J., Bach, F., Ponce, J., & Sapiro, G. 2010. Online Lea rning for Matrix Factorization and
Sparse Coding. Journal of Machine Learning Research, 11(Aug.), 19–60.
Makeig, S., Bell, A.J., Jung, T.P., Sejnowski, T.J., et al. 1996. Independent component analysis of
electroencephalographic data. Advances in neural information processing systems, 145–151.
Mallat, S.G., & Zhang, Z. 1993. Matching pursuits with time- frequency dictionaries. Signal Pro-
cessing, IEEE Transactions on, 41(12), 3397–3415.
Mazaheri, A., & Jensen, O. 2008. Amplitude modulations of brain oscillations generate slow evoked
responses. Journal of Neuroscience, 28, 7781–7787.
Olshausen, B.A., & Field, D.J. 1997. Sparse coding with an ov ercomplete basis set: a strategy
employed by V1? Vision research, 37(23), 3311–25.
Pearson, K. 1901. LIII. On lines and planes of closest ﬁt to systems of points in space. The London,
Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559–572.
Smith, Evan, & Lewicki, Michael S. 2005. Efﬁcient coding of t ime-relative structure using spikes.
Neural Computation, 17(1), 19–45.
Tibshirani, R. 1996. Regression shrinkage and selection via the lasso.Journal of the Royal Statistical
Society. Series B (Methodological), 267–288.
Woody, C.D. 1967. Characterization of an adaptive ﬁlter for the analysis of variable latency neuro-
electric signals. Medical and Biological Engineering and Computing, 5(6), 539–554.
11


Peer review:
[Review]:
We thank the reviewers for their constructive comments. We submitted a new version of the paper to arXiv, which should be made available on Wednesday, March 13. As one major change we now point out the similarity to convolutional/shift-invariant sparse coding (SISC)*, but also mention the differences mainly introduced by the l_0 sparsity constraint. A new contribution is an analysis of the algorithm's complexity as well as possibilities for speed ups – although the computation time was already low for the conducted experiments, this could become an issue for real-time analysis. The changes in detail:

                  
[1] Smith, Evan, and Michael S. Lewicki. 'Efficient coding of time-relative structure using spikes.' Neural Computation 17.1 (2005): 19-45.
[2] Blumensath, Thomas, and Davies, Mike. 'Sparse and shift-invariant representations of music.' Audio, Speech, and Language Processing, IEEE Transactions on 14.1 (2006): 50-57.
[3] R. Grosse, R. Raina, H. Kwong, and AY Ng, 'Shift-invariant sparse coding for audio classification,' in Proceedings of the Twenty-third Conference on Uncertainty in Artificial Intelligence (UAI'07), 2007
[4] Ekanadham, Chaitanya, Daniel Tranchina, and Eero P. Simoncelli. 'Sparse decomposition of transformation-invariant signals with continuous basis pursuit.' Acoustics, Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on. IEEE, 2011.

Introduction
The second part of the introduction has been rewritten. Shift-invariant sparse coding (SISC) is introduced and its differences to JADL are pointed out. Most significant is the constraint in JADL, that only one shifted version of each atom shall be active per signal. As a consequence, JADL leads to a less complex algorithm (in both, sparse coding and dictionary update step), which in contrast to SISC does not need heuristic preselection of active atoms. In addition, we remarked that JADL is designed to learn only few atoms, in contrast to most dictionary learning applications. Hence, the term “sparsity” only makes sense with respect to the “unrolled” dictionary. However, for the most part this sparsity is achieved by the explicit constraint, not by sparse regularization.

Section 3, JADL
In the sparse coding section, the computational advantages of the modified LARS have been pointed out and contrasted with SISC.

In the dictionary update section, we noted that the JADL formulation leads to an update of same complexity as regular DL; this does not apply to SISC. This fact could be used by changing the ratio of sparse coding and dictionary update steps in favor to the latter, i.e. by employing mini-batch or online techniques.

Section 4, Experiments
For the real data, the computation time has been investigated for different K (number of atoms) and S (number of allowed shifts). The curve shows linear correlation with S; when increasing K, however, the computation time increases more than linear. This is due to the following: While both S and K affect the size of the unrolled dictionary, an increase in S is handled efficiently in the JADL formulation, as mentioned in the sections sparse coding and dictionary update above. We also mentioned that the computation time for the conducted experiments was very small (4.3 seconds for the real data), hence computational complexity should not become an issue for offline analysis. Employing the proposed speed ups and a further optimization of the code, on the other hand, could even allow for real time analysis (this could be desirable especially for M/EEG-based brain computer interfaces (BCI)) .

Detailed responses to the reviewers' comments
-------------------------------------------------------------

Anonymous 8ed7
-------------------------

Cons
1. [Computational requirements] As mentioned above, we investigated computation time empirically for different values of S and K. Even for K=15 and S=200 (i.e. 3000 elements in the unrolled dictionary ) the computation time remained less than 1 minute for 200 iterations. Therefore, the increased complexity should only matter for real-time applications, for which we proposed several speed-ups in the formulation of the algorithm.

2. (a) [Examples for shifts] We added a comment (footnote) on possible definitions of the shifts in the problem statement, there are different ways how to define the shifts at the borders. As the JADL framework is general enough to allow even arbitrary linear transforms, we do not specify a certain definition at this point. A detailed discussion of the right way to handle boundary effects would be out of the scope of the paper; however, we found that in our experiments the choice of the shifts did not affect the outcome significantly.  

(b) [Examples for types of data that can be analyzed with JADL] In the introduction, we now mention similar properties to those of neuroelectric data in different bioelectric or biomagnetic signals, such as ECG, EMG. [explain “features well-aligned across signals”] We changed this formulation to “each waveform of interest occurs approximately at the same time in every trial”. Is this sufficiently clear?
(c) [Improve explanation how to enforce constraint (7)] changed as suggested

3. [On the importance of lambda] We agree that if the number of atoms is large, the parameter lambda has an important role to ensure sparsity. However, JADL is designed to learn only a small number K of atoms. This is due to several reasons: (i) the jitter-adaptivity ensures a compact representation, as a waveform that is shifted throughout signals can be encoded in a single atom; (ii) for the applications JADL is aimed at, it is either not desired or not feasible to learn many atoms, since (a) the dictionary should be easily interpretable and reveal the main activity across signals (b) the number of training examples M is limited and K<<M must hold to prevent from overfitting, which is a particularily critical aspect due to often high noise level. A similar comment on the different use of “sparsity” in JADL has been added to the introduction to get this point clear.

4. [Comparison to common DL with large K and large lambda] We agree that a comparison to DL with a similar K used as in JADL would not be fair. For the simulated data, different values of K have therefore been used, and lambda has been optimized individually for each K to yield the smallest error/highest similarity w.r.t. the ground truth. The table stops at K=12; for larger K performance becomes worse for all three methods. This is now pointed out in the new version.

Minor Comments

2. We changed the claim from “the biggest challenge” to “an important challenge” and provided references.
3. changed as suggested
4. changed as suggested


Anonymous 5e7a
-------------------------

We agree that a comparison to previous work on convolutional/shift-invariant sparse coding is necessary and hope that the changes made as described above make the similarities and differences to between SISC and JADL sufficiently clear. We found that most papers on SISC do not address the problem of the dictionary update but only focus on sparse coding. An exception was [3], from which it becomes clear that for SISC the dictionary update is a non-trivial problem with increased computational complexity.

Anonymous 8b9c
-------------------------

See comment above for comparison to SISC.

In fact, the LFP data does not contain much more structure than the spikes. Hence, the learned dictionaries look quite redundant and its analysis provides only limited insight in the data. However, we think that the visualization of the code reveals that although the dictionary looks redundant, important differences have been picked up, leading to contiguous sets of epochs dominated by the same atom.