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]:
This paper introduces a dictionary learning technique that incorporates time delays or shifts on the learned dictionary, called JADL, to better account for this structure in multi-trial neuroelectric signals. The algorithm uses the previous dictionary learning framework and non-convex optimization, but adds a selection step over possible shifts for each atom (for each point), framed as an l-0 optimization. This objective is the main contribution of the paper, which enables better performance for time-delayed data as well as potentially useful temporal structure to be extracted from the data.

The paper introduces a novel objective for addressing the time shift problem (e.g in M/EEG data), but frames a typical coordinate descent approach for solving the resulting non-convex problem.  The main difference in the optimization is (1) ensuring that the coefficients, a, for the dictionary, D, block-wise satisfy the l-0 constraint by disabling updates to all but one coefficient within a block and (2) modifying the gradient update in block coordinate descent on the dictionary, D, which now has a shift operator around D. Taking this obvious solution route leads to a non-convex optimization and potentially lengthy computation (as the delta set size increases). The quality of writing and experiments is high.

Pros
1. The proposed JADL algorithm facilitates application of dictionary learning techniques to M/EEG data, which is an important application. Moreover, as a secondary benefit, it allows time-delay structure to be learned.

2. The writing is mostly clear and the paper is well organized.

3. Experimental results are comprehensive and include important details of their experimental procedure.

Cons
1. The computational requirements of this algorithm are not explored, though the larger dictionary in JADL (due to the addition of delta shifts) could significantly slow learning.

2. For clarity: (a) Include examples of shifts, Delta, in the problem statement (such as the ones used in the experiments). (b) Include examples of the types of data that could benefit from this framework, to better justify the importance of framing the problem with time-shifts and better explain what is meant by 'features [being] well-aligned across signals'.  (c ) The explanation of how to enforce constraint (7) should be improved, e.g. 'block all other coefficients a_j^{S,i}' should probably be 'block all other coefficients in segment a_j^{S,i}', but the meaning is actually significantly different and this was quite confusing.

3. The comment that the parameter, lambda, is no longer important because sparsity is induced by the constraint in (7) suggests that as the size of the set of delta increases, this problem formulation no longer learns a sparse solution over the chosen dictionary. I would suggest that this is not the case, but rather that the datasets used in this paper had a small dictionary and did not require the coefficients to be sparse. The constraint in (7) simply ensures that only one delta is chosen per atom, but does not guarantee that the final solution over the delta-shifted dictionary will be sparse. Therefore, if the number of atoms is large, the regularizer || a_j ||_1 should still be important. It is true that constraint (7) ensures the solution is sparse over all possible delta-shifted dictionaries; this is however an unfair comparison to other dictionary learning techniques which have a much smaller dictionary space to weight over.

4. Con (3) suggests that learning with a very large dictionary (the size of the delta-shifted dictionary set) and setting the lambda parameter large might have more comparable performance to the algorithm suggested in this paper and should be included. Of course, this highly regularized approach on a large dictionary would not explicitly provide the time shift structure in the data as does JADL, but would be an interesting and more fair comparison. However, if the time-shift structure is not actually useful (and is simply used to improve learning), then DL with a large dictionary and large regularization parameter, lambda, could be all that is needed to deal with this problem for EEG data. The authors should clarify this difference and contribution more clearly.


Minor Comments:
1. For a reference on convex solution to the dictionary learning problem, see 'Convex sparse matrix factorizations', F. Bach, J. Mairal and J. Ponce. 2008; and 'Convex Sparse Coding, Subspace Learning and Semi-Supervised Extensions', X. Zhang, Y. Yu, M. White, R. Huang and D. Schuurmans. 2011.
2. There should be citations for the claim: 'This issue is currently the biggest challenge
in M/EEG multi-trial analysis.'
3. Bottom of page 3: 'which allows to solve it' -> 'which allows us to solve it'
4. Page 5: 'property allows to' -> 'property allows us to'