TOOLS AND RESOURCES
Unsupervised discovery of temporal
sequences in high-dimensional datasets,
with applications to neuroscience
Emily L Mackevicius1†, Andrew H Bahle1†, Alex H Williams2, Shijie Gu1,3,
Natalia I Denisenko1, Mark S Goldman4,5*, Michale S Fee1*
1McGovern Institute for Brain Research, Department of Brain and Cognitive
Sciences, Massachusetts Institute of Technology, Cambridge, United States;
2Neurosciences Program, Stanford University, Stanford, United States; 3School of
Life Sciences and Technology, ShanghaiTech University, Shanghai, China; 4Center
for Neuroscience, Department of Neurobiology, Physiology and Behavior,
University of California, Davis, Davis, United States; 5Department of Ophthamology
and Vision Science, University of California, Davis, Davis, United States
Abstract Identifying low-dimensional features that describe large-scale neural recordings is a
major challenge in neuroscience. Repeated temporal patterns (sequences) are thought to be a
salient feature of neural dynamics, but are not succinctly captured by traditional dimensionality
reduction techniques. Here, we describe a software toolbox—called seqNMF—with new methods
for extracting informative, non-redundant, sequences from high-dimensional neural data, testing
the significance of these extracted patterns, and assessing the prevalence of sequential structure in
data. We test these methods on simulated data under multiple noise conditions, and on several
real neural and behavioral data sets. In hippocampal data, seqNMF identifies neural sequences that
*For correspondence: match those calculated manually by reference to behavioral events. In songbird data, seqNMF
msgoldman@ucdavis.edu (MSG); discovers neural sequences in untutored birds that lack stereotyped songs. Thus, by identifying
fee@mit.edu (MSF) temporal structure directly from neural data, seqNMF enables dissection of complex neural circuits
†These authors contributed without relying on temporal references from stimuli or behavioral outputs.
equally to this work DOI: https://doi.org/10.7554/eLife.38471.001
Competing interests: The
authors declare that no
competing interests exist. Introduction
Funding: See page 29 The ability to detect and analyze temporal sequences embedded in a complex sensory stream is an
Received: 19 May 2018 essential cognitive function, and as such is a necessary capability of neuronal circuits in the brain
Accepted: 04 January 2019 (Clegg et al., 1998; Janata and Grafton, 2003; Bapi et al., 2005; Hawkins and Ahmad, 2016), as
Published: 05 February 2019 well as artificial intelligence systems (Cui et al., 2016; Sutskever et al., 2014). The detection and
characterization of temporal structure in signals is also useful for the analysis of many forms of physi-
Reviewing editor: Laura Colgin,
cal and biological data. In neuroscience, recent advances in technology for electrophysiological and
The University of Texas at Austin,
United States optical measurements of neural activity have enabled the simultaneous recording of hundreds or
thousands of neurons (Chen et al., 2013; Kim et al., 2016; Scholvin et al., 2016; Jun et al., 2017),
Copyright Mackevicius et al.
in which neuronal dynamics are often structured in sparse sequences (Hahnloser et al., 2002;
This article is distributed under
Harvey et al., 2012; MacDonald et al., 2011; Okubo et al., 2015; Fujisawa et al., 2008;
the terms of the Creative
Commons Attribution License, Pastalkova et al., 2008). Such sequences can be identified by averaging across multiple trials, but
which permits unrestricted use only in cases where an animal receives a temporally precise sensory stimulus, or executes a suffi-
and redistribution provided that ciently stereotyped behavioral task.
the original author and source are Neural sequences have been hypothesized to play crucial roles over a much broader range of nat-
credited. ural settings, including during learning, sleep, or diseased states (Mackevicius and Fee, 2018). In
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 1 of 42
Tools and resources Neuroscience
these applications, it may not be possible to use external timing references, either because behav-
iors are not stereotyped or are entirely absent. Thus, sequences must be extracted directly from the
neuronal data using unsupervised learning methods. Commonly used methods of this type, such as
principal component analysis (PCA) or clustering methods, do not efficiently extract sequences,
because they typically only model synchronous patterns of activity, rather than extended spatio-tem-
poral motifs of firing.
Existing approaches that search for repeating neural patterns require computationally intensive
or statistically challenging analyses (Brody, 1999; Mokeichev et al., 2007; Quaglio et al., 2018;
Brunton et al., 2016). While progress has been made in analyzing non-synchronous sequential pat-
terns using statistical models that capture cross-correlations between pairs of neurons (Russo and
Durstewitz, 2017; Gerstein et al., 2012; Schrader et al., 2008; Torre et al., 2016;
Grossberger et al., 2018; van der Meij and Voytek, 2018), such methods may not have statistical
power to scale to patterns that include many (more than a few dozen) neurons, may require long
periods (105 timebins) of stationary data, and may have challenges in dealing with (non-sequential)
background activity. For a review highlighting features and limitations of these methods see
(Quaglio et al., 2018).
Here, we explore a complementary approach, which uses matrix factorization to reconstruct neu-
ral dynamics using a small set of exemplar sequences. In particular, we build on convolutional non-
negative matrix factorization (convNMF) (Smaragdis, 2004; Smaragdis, 2007) (Figure 1B), which
has been previously applied to identify recurring motifs in audio signals such as speech
(O’Grady and Pearlmutter, 2006; Smaragdis, 2007; Vaz et al., 2016), as well as neural
signals (Peter et al., 2017). ConvNMF identifies exemplar patterns (factors) in conjunction with the
times and amplitudes of pattern occurrences. This strategy eliminates the need to average activity
aligned to any external behavioral references.
While convNMF may produce excellent reconstructions of the data, it does not automatically pro-
duce the minimal number of factors required. Indeed, if the number of factors in the convNMF
model is greater than the true number of sequences, the algorithm returns overly complex and
redundant factorizations. Moreover, in these cases, the sequences extracted by convNMF will often
be inconsistent across optimization runs from different initial conditions, complicating scientific inter-
pretations of the results (Peter et al., 2017; Wu et al., 2016).
To address these concerns, we developed a toolbox of methods, called seqNMF, which includes
two different strategies to resolve the problem of redundant factorizations described above. In addi-
tion, the toolbox includes methods for promoting potentially desirable features such as orthogonal-
ity or sparsity of the spatial and temporal structure of extracted factors, and methods for analyzing
the statistical significance and prevalence of the identified sequential structure. To assess these
tools, we characterize their performance on synthetic data under a variety of noise conditions and
also show that they are able to find sequences in neural data collected from two different animal
species using different behavioral protocols and recording technologies. Applied to extracellular
recordings from rat hippocampus, seqNMF identifies neural sequences that were previously found
by trial-averaging. Applied to functional calcium imaging data recorded in vocal/motor cortex of
untutored songbirds, seqNMF robustly identifies neural sequences active in a biologically atypical
and overlapping fashion. This finding highlights the utility of our approach to extract sequences with-
out reference to external landmarks; untutored bird songs are so variable that aligning neural activity
to song syllables would be difficult and highly subjective.
Results
Matrix factorization framework for unsupervised discovery of features
in neural data
Matrix factorization underlies many well-known unsupervised learning algorithms, including PCA
(Pearson, 1901), non-negative matrix factorization (NMF) (Lee and Seung, 1999), dictionary learn-
ing, and k-means clustering (see Udell et al., 2016 for a review). We start with a data matrix, X, con-
taining the activity of N neurons at T timepoints. If the neurons exhibit a single repeated pattern of
synchronous activity, the entire data matrix can be reconstructed using a column vector w represent-
ing the neural pattern, and a row vector h representing the times and amplitudes at which that
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 2 of 42
Tools and resources Neuroscience
Factorization Strategies ConvNMF Redundancies
A NMF: C
. . . Type 1
=
Time
B convNMF:
Type 2
. . .
=
Time
Type 3
Figure 1. Convolutional NMF factorization. (A) NMF (non-negative matrix factorization) approximates a data
matrix describing the activity of N neurons at T timepoints as a sum of K rank-one matrices. Each matrix is
generated as the outer product of two nonnegative vectors: wk of length N, which stores a neural ensemble, and
hk of length T , which holds the times at which the neural ensemble is active, and the relative amplitudes of this
activity. (B) Convolutional NMF also approximates an N  T data matrix as a sum of K matrices. Each matrix is
generated as the convolution of two components: a non-negative matrix wk of dimension N  L that stores a
sequential pattern of the N neurons at L lags, and a vector of temporal loadings, hk , which holds the times at
which each factor pattern is active in the data, and the relative amplitudes of this activity. (C) Three types of
inefficiencies present in unregularized convNMF: Type 1, in which two factors are used to reconstruct the same
instance of a sequence; Type 2, in which two factors reconstruct a sequence in a piece-wise manner; and Type 3,
in which two factors are used to reconstruct different instances of the same sequence. For each case, the factors
(W and H) are shown, as well as the reconstruction (Xe ¼W
 H ¼ w 
1  h þw 
1 2  h2 þ 
 
 
).
DOI: https://doi.org/10.7554/eLife.38471.002
The following figure supplement is available for figure 1:
Figure supplement 1. Quantifying the effect of different penalties on convNMF.
DOI: https://doi.org/10.7554/eLife.38471.003
pattern occurs (temporal loadings). In this case, the data matrix X is mathematically reconstructed
as the outer product of w and h. If multiple component patterns are present in the data, then each
pattern can be reconstructed by a separate outer product, where the reconstructions are summed
to approximate the entire data matrix (Figure 1A) as follows:
XK
Xnt »X
e
nt ¼ WnkHkt ¼ ðWHÞnt (1)
k¼1
where thXnt is the ðntÞ element of matrix X, that is, the activity of neuron n at time t. Here, in order
to store K different patterns, W is a NK matrix containing the K exemplar patterns, and H is a
KT matrix containing the K timecourses:
2 3
2 3  h1 
j j j 6 7
6 7 6 h2 7
W¼ 4w1 w 6 72 
 
 
 wK 5 H¼
46 .. 57
(2)
j j j .
 hK 
Given a data matrix with unknown patterns, the goal of matrix factorization is to discover a small
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 3 of 42
Neurons Neurons
Tools and resources Neuroscience
set of patterns, W, and a corresponding set of temporal loading vectors, H, that approximate the
data. In the case that the number of patterns, K, is sufficiently small (less than N and T), this corre-
sponds to a dimensionality reduction, whereby the data is expressed in more compact form. PCA
additionally requires that the columns of W and the rows of H are orthogonal. NMF instead requires
that the elements of W and H are nonnegative. The discovery of unknown factors is often accom-
plished by minimizing the following cost function, which measures the element-by-element sum of
all squared errors betwqeffiffieffiffinffiffiffiffiffiaffiffiffiffirffieffi construction X
e ¼WH and the original data matrix X using the Fro-
P
benius norm, k 2MkF ¼ ijMij:
ð  W ;H Þ ¼ arg minkXe  k2X F
W;H
(3)
(Note that other loss functions may be substituted if desired, for example to better reflect the
noise statistics; see (Udell et al., 2016) for a review). The factors W and H that minimize this cost
function produce an optimal reconstruction e  ¼  X W H . Iterative optimization methods such as gra-
dient descent can be used to search for global minima of the cost function; however, it is often pos-
sible for these methods to get caught in local minima. Thus, as described below, it is important to
run multiple rounds of optimization to assess the stability/consistency of each model.
While this general strategy works well for extracting synchronous activity, it is unsuitable for dis-
covering temporally extended patterns—first, because each element in a sequence must be repre-
sented by a different factor, and second, because NMF assumes that the columns of the data matrix
are independent ‘samples’ of the data, so permutations in time have no effect on the factorization of
a given data. It is therefore necessary to adopt a different strategy for temporally extended features.
Convolutional matrix factorization
Convolutional nonnegative matrix factorization (convNMF) (Smaragdis, 2004; Smaragdis, 2007)
extends NMF to provide a framework for extracting temporal patterns, including sequences, from
data. While in classical NMF each factor W is represented by a single vector (Figure 1A), the factors
W in convNMF represent patterns of neural activity over a brief period of time. Each pattern is
stored as an N  L matrix, wk, where each column (indexed by ‘ ¼ 1 to L) indicates the activity of
neurons at different timelags within the pattern (Figure 1B). The times at which this pattern/
sequence occurs are encoded in the row vector h1, as for NMF. The reconstruction is produced by
convolving the N  L pattern with the time series h1 (Figure 1B).
If the data contains multiple patterns, each pattern is captured by a different N  L matrix and a
different associated time series vector h. A collection of K different patterns can be compiled
together into an N  K  L array (also known as a tensor), W and a corresponding K  T time series
matrix H. Analogous to NMF, convNMF generates a reconstruction of the data as a sum of K convo-
lutions between each neural activity pattern (W), and its corresponding temporal loadings (H):
e X
K XL1
Xnt »Xnt ¼ Wnk‘Hkðt‘Þ  ðW
 HÞnt (4)
k¼1 ‘¼0
The tensor/matrix convolution operator 
 (notation summary, Table 1) reduces to matrix multi-
plication in the L¼ 1 case, which is equivalent to standard NMF. The quality of this reconstruction
can be measured using the same cost function shown in Equation 3, and W and H may be found
iteratively using similar multiplicative gradient descent updates to standard NMF (Lee and Seung,
1999; Smaragdis, 2004; Smaragdis, 2007).
While convNMF can perform extremely well at reconstructing sequential structure, it can be chal-
lenging to use when the number of sequences in the data is not known (Peter et al., 2017). In this
case, a reasonable strategy would be to choose K at least as large as the number of sequences that
one might expect in the data. However, if K is greater than the actual number of sequences,
convNMF often identifies more significant factors than are minimally required. This is because each
sequence in the data may be approximated equally well by a single sequential pattern or by a linear
combination of multiple partial patterns. A related problem is that running convNMF from different
random initial conditions produces inconsistent results, finding different combinations of partial
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 4 of 42
Tools and resources Neuroscience
Table 1. Notation for convolutional matrix factorization
Shift operator
‘!
The operator ðHÞ shifts a matrix H in the ! direction by ‘ timebins:
‘!  ‘
ðHÞ
t ¼ H
ðtlÞ and likewise ðHÞ
t ¼ H
ðtþ‘Þ
where 
 indicates all elements along the respective matrix dimension.
The shift operator inserts zeros when ðt  ‘Þ< 0 or ðt þ ‘Þ> T
Tensor convolution operator
Convolutive matrix factorization reconstructs a data matrix X
using a N  K  L tensor W and a K  T matrix H:
e P ‘!X ¼W
H ¼ ‘ W

‘H
Note that each neuron n is reconstructed as the sum of k convolutions:
Xe
P P
nt ¼ k ‘ Wnk‘Hkðt‘Þ  ðW
HÞnt
Transpose tensor convolution operator
The following quantity is useful in several contexts:
> P  ‘
W
 X ¼ ‘ðW

‘Þ
>
X
> P P P
Note that each element ð 
 Þ ¼ ð Þ>W X kt l W
k‘ X
ðtþ‘Þ ¼ l n Wnk‘Xnðtþ‘Þ measures
the overlap (correlation) of factor k with the data at time t
convNMF reconstruction
P
X»Xe ¼ k W
k

 Hk
 ¼W
 H
Note that NMF is a special case of convNMF, where L ¼ 1
L1 entrywise norm excluding diagonal elements
P P
For any K  K matrix C, kCk
1;i 6¼j  k j¼6 k Cjk
Special matrices
1 is a K  K matrix of ones
I is the K  K identity matrix
S is a T  T smoothing matrix: Sij ¼ 1 when ji jj< L and otherwise Sij ¼ 0
DOI: https://doi.org/10.7554/eLife.38471.004
patterns on each run (Peter et al., 2017). These inconsistency errors fall into three main categories
(Figure 1C):
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 5 of 42
Tools and resources Neuroscience
. Type 1: Two or more factors are used to reconstruct the same instances of a sequence.
. Type 2: Two or more factors are used to reconstruct temporally different parts of the same
sequence, for instance the first half and the second half.
. Type 3: Duplicate factors are used to reconstruct different instances of the same sequence.
Together, these inconsistency errors manifest as strong correlations between different redundant
factors, as seen in the similarity of their temporal loadings (H) and/or their exemplar activity patterns
(W).
We next describe two strategies for overcoming the redundancy errors described above. Both
strategies build on previous work that reduces correlations between factors in NMF. The first strat-
egy is based on regularization, a common technique in optimization that allows the incorporation of
constraints or additional information with the goal of improving generalization performance or sim-
plifying solutions to resolve degeneracies (Hastie et al., 2009). A second strategy directly estimates
the number of underlying sequences by minimizing a measure of correlations between factors (stabil-
ity NMF; Wu et al., 2016).
Optimization penalties to reduce redundant factors
To reduce the occurrence of redundant factors (and inconsistent factorizations) in convNMF, we
sought a principled way of penalizing the correlations between factors by introducing a penalty
term, R, into the convNMF cost function:
 
ð W ; H Þ ¼ arg min k
2
Xe XkF þR (5)
W;H
Regularization has previously been used in NMF to address the problem of duplicated factors,
which, similar to Type 1 errors above, present as correlations between the H’s (Choi, 2008;
Chen and Cichocki, 2004). Such correlations are measured by computing the correlation matrix
>
HH , which contains the correlations between the temporal loadings of every pair of factors. The
regularization may be implemented using the penalty term R¼ lk >HH k ;i¼6 j, where the seminorm1
k 
 k ;i 6¼j sums the absolute value of every matrix entry except those along the diagonal (notation sum-1
mary, Table 1) so that correlations between different factors are penalized, while the correlation of
each factor with itself is not. Thus, during the minimization process, similar factors compete, and a
larger amplitude factor drives down the temporal loading of a correlated smaller factor. The param-
eter l controls the magnitude of the penalty term R.
In convNMF, a penalty term based on >HH yields an effective method to prevent errors of Type
1, because it penalizes the associated zero lag correlations. However, it does not prevent errors of
the other types, which exhibit different types of correlations. For example, Type 2 errors result in
correlated temporal loadings that have a small temporal offset and thus are not detected by >HH .
One simple way to address this problem is to smooth the H’s in the penalty term with a square win-
dow of length 2L 1 using the smoothing matrix S (Sij ¼ 1 when ji jj<L and otherwise Sij ¼ 0). The
resulting penalty, R ¼ lk >HSH k, allows factors with small temporal offsets to compete, effectively
preventing errors of Types 1 and 2.
This penalty does not prevent errors of Type 3, in which redundant factors with highly similar pat-
terns in W are used to explain different instances of the same sequence. Such factors have temporal
loadings that are segregated in time, and thus have low correlations, to which the cost term
k >HSH k is insensitive. One way to resolve errors of Type 3 might be to include an additional cost
term that penalizes the similarity of the factor patterns in W. This has the disadvantage of requiring
an extra parameter, namely the l associated with this cost.
Instead we chose an alternative approach to resolve errors of Type 3 that simultaneously detects
correlations in W and H using a single cross-orthogonality cost term. We note that, for Type 3
errors, redundant W patterns have a high degree of overlap with the data at the same times, even
though their temporal loadings are segregated at different times. To introduce competition between
these factors, we first compute, for each pattern in W, its overlap with the data at time t. This quan-
>
tity is captured in symbolic form by W
 X (see Table 1). We then compute the pairwise correlation
between the temporal loading of each factor and the overlap of every other factor with the data.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 6 of 42
Tools and resources Neuroscience
This cross-orthogonality penalty term, which we refer to as ’x-ortho’, sums up these correlations
across all pairs of factors, implemented as follows:
>
R¼ lkð >W
 XÞSH k ;i 6¼j (6)1
When incorporated into the update rules, this causes any factor that has a high overlap with the
data to suppress the temporal loadings (H) of any other factors that have high overlap with the data
at that time (Further analysis, Appendix 2). Thus, factors compete to explain each feature of the
Table 2. Regularized NMF and convNMF: cost functions and algorithms
NMF
L ¼ 1 jjXe >Xjj2 þR W W XH
2 2 e >XH þ qR
qW
>
  W XH H
Xe ¼WH > qRW Xeþ
qH
convNMF
>
1 2 ‘!
L ¼ jjXe Xjj þR
2 2 X HW

‘  W

‘ 
Xe ‘!
>
H þ qR
qW

‘
Xe ¼W
 H >
 W
 XH H >
W
 XeþqR
qH
L1 regularization for H ( L1 for W is analogous)
R ¼ ljjHjj qR
1 ¼ 0qW

‘
qR ¼ l1
qH
Orthogonality cost for H
l
R ¼ jj >jj qRHH
2 1;i¼6 j ¼ 0qW

‘
qR ¼ lð1 IÞH
qH
Smoothed orthogonality cost for H (favors ‘events-based’)
¼ l jj > qRR HSH jj
2 1;i 6¼j ¼ 0qW

‘
qR ¼ lð1 IÞHS
qH
Smoothed orthogonality cost for W (favors ‘parts-based’)
¼ lR jj > qRW
2 flat
Wflat jj ;i6¼j ¼ lW ð1 IÞ1 qW flat

‘
where P qR ¼ 0qH
ðWflatÞnk ¼ ‘ Wnk‘
Smoothed cross-factor orthogonality (x-ortho penalty)
>  ‘
R ¼ ljjðW
 XÞ >SH jj qR >
1;i¼6 j ¼ lXSH ð1 IÞqW

‘
>
qR ¼ lð1 IÞW
 XS
qH
DOI: https://doi.org/10.7554/eLife.38471.005
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 7 of 42
Tools and resources Neuroscience
data, favoring solutions that use a minimal set of factors to give a good reconstruction. The resulting
global cost function is:
 
>
ð W ; H Þ ¼ 2 >arg min kXe XkF þlkðW
 XÞSH k1;i¼6 j (7)
W;H
The update rules for W and H are based on the derivatives of this global cost function, leading
to a simple modification of the standard multiplicative update rules used for NMF and convNMF
(Lee and Seung, 1999; Smaragdis, 2004; Smaragdis, 2007) (Table 2). Note that the addition of
this cross-orthogonality term does not formally constitute regularization, because it also includes a
contribution from the data matrix X, rather than just the model variables W and H. However, at
least for the case that the data is well reconstructed by the sum of all factors, the x-ortho penalty
can be shown to be approximated by a formal regularization (Appendix 2). This formal regularization
contains both a term corresponding to a weighted smoothed orthogonality penalty on W and a
term corresponding to a weighted smoothed) orthogonality penalty on H, consistent with the obser-
vation that the x-ortho penalty simultaneously punishes factor correlations in W and H.
There is an interesting relation between our method for penalizing correlations and other meth-
ods for constraining optimization, namely sparsity. Because of the non-negativity constraint imposed
in NMF, correlations can also be reduced by increasing the sparsity of the representation. Previous
efforts have been made to minimize redundant factors using sparsity constraints; however, this
approach may require penalties on both W and H, necessitating the selection of two hyper-parame-
ters (lw and lh) (Peter et al., 2017). Since the use of multiple penalty terms increases the complexity
of model fitting and selection of parameters, one goal of our work was to design a simple, single
penalty function that could regularize both W and H simultaneously. The x-ortho penalty described
above serves this purpose (Equation 6). As we will describe below, the application of sparsity penal-
ties can be very useful for shaping the factors produced by convNMF, and our code includes options
for applying sparsity penalties on both W and H.
Extracting ground-truth sequences with the x-ortho penalty when the
number of sequences is not known
We next examined the effect of the x-ortho penalty on factorizations of sequences in simulated
data, with a focus on convergence, consistency of factorizations, the ability of the algorithm to dis-
cover the correct number of sequences in the data, and robustness to noise (Figure 2A). We first
assessed the model’s ability to extract three ground-truth sequences lasting 30 timesteps and con-
taining 10 neurons in the absence of noise (Figure 2A). The resulting data matrix had a total dura-
tion of 15,000 timesteps and contained on average 60±6 instances of each sequence. Neural
activation events were represented with an exponential kernel to simulate calcium imaging data. The
algorithm was run with the x-ortho penalty for 1000 iterations andit reliably converged to a root-
mean-squared-error (RMSE) close to zero (Figure 2B). RMSE reached a level within 10% of the
asymptotic value in approximately 100 iterations.
While similar RMSE values were achieved using convNMF with and without the x-ortho penalty;
the addition of this penalty allowed three ground-truth sequences to be robustly extracted into
three separate factors (w1, w2, and w3 in Figure 2A) so long as K was chosen to be larger than the
true number of sequences. In contrast, convNMF with no penalty converged to inconsistent factori-
zations from different random initializations when K was chosen to be too large, due to the ambigui-
ties described in Figure 1. We quantified the consistency of each model (see
Materials and methods), and found that factorizations using the x-ortho penalty demonstrated near
perfect consistency across different optimization runs (Figure 2C).
We next evaluated the performance of convNMF with and without the x-ortho penalty on data-
sets with a larger number of sequences. In particular, we set out to observe the effect of the x-ortho
penalty on the number of statistically significant factors extracted. Statistical significance was deter-
mined based on the overlap of each extracted factor with held out data (see Materials and methods
and code package). With the penalty term, the number of significant sequences closely matched the
number of ground-truth sequences. Without the penalty, all 20 extracted sequences were significant
by our test (Figure 2D).
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 8 of 42
Tools and resources Neuroscience
A h1 B
h2 1
h3 10
w1 w2 w3 0
10
-1
10
-2
10
0
10 1 2
3
10 10 10
Time Iterations
C D
1 20X-ortho
convNMF
0.8
15
0.6
10
0.4
5
0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10
Consistency # Sequences
E F                   
Figure 2. Effect of the x-ortho penalty on the factorization of sequences. (A) A simulated dataset with three sequences. Also shown is a factorization
with x-ortho penalty (K ¼ 20, L ¼ 50, l ¼ 0:003). Each significant factor is shown in a different color. At left are the exemplar patterns (W) and on top
are the timecourses (H). (B) Reconstruction error as a function of iteration number. Factorizations were run on a simulated dataset with three sequences
and 15,000 timebins ( » 60 instances of each sequence). Twenty independent runs are shown. Here, the algorithm converges to within 10% of the
asymptotic error value within » 100 iterations. (C) The x-ortho penalty produces more consistent factorizations than unregularized convNMF across 400
independent fits (K ¼ 20, L ¼ 50, l ¼ 0:003). (D) The number of statistically significant factors (Figure 2—figure supplement 1) vs. the number of
ground-truth sequences for factorizations with and without the x-ortho penalty. Shown for each condition is a vertical histogram representing the
number of significant factors over 20 runs (K ¼ 20, L ¼ 50, l ¼ 0:003). (E) Factorization with x-ortho penalty of two simulated neural sequences with
shared neurons that participate at the same latency. (F) Same as E but for two simulated neural sequences with shared neurons that participate at
different latencies.
DOI: https://doi.org/10.7554/eLife.38471.006
The following figure supplements are available for figure 2:
Figure supplement 1. Outline of the procedure used to assess factor significance.
DOI: https://doi.org/10.7554/eLife.38471.007
Figure supplement 2. Number of significant factors as a function of l for datasets containing between 1 and 10 sequences.
DOI: https://doi.org/10.7554/eLife.38471.008
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 9 of 42
Probability Neurons
# Significant factors
RMSE
Tools and resources Neuroscience
We next considered how the x-ortho penalty performs on sequences with more complex structure
than the sparse uniform sequences of activity ediscussed above. We further examined the case in
which a population of neurons is active in multiple different sequences. Such neurons that are shared
across different sequences have been observed in several neuronal datasets (Okubo et al., 2015;
Pastalkova et al., 2008; Harvey et al., 2012). For one test, we constructed two sequences in which
shared neurons were active at a common pattern of latencies in both sequences; in another test,
shared neurons were active in a different pattern of latencies in each sequence. In both tests, factori-
zations using the x-ortho penalty achieved near-perfect reconstruction error, and consistency was
A Probabilistic Participation B Additive Noise C Temporal Jitter D Time Warping
1 1 1 1
0.8 0.8 0.8 0.8
0.6 0.6 0.6 0.6
0.4 0.4 0.4 0.4
0.2 0.2 0.2 0.2
0 0 0 0
0 10 20 30 40 50 60 70 80 90 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 5 10 15 20 25 30 35 40 45 1 2 3 4 5 6 7 8 9 10
% Dropout % Additive noise Jitter SD Warping condition
10 40 80 0.5 2 4 5 20 40 2 5 9
Figure 3. Testing factorization performance on sequences contaminated with noise. Performance of the x-ortho penalty was tested under four different
noise conditions: (A) probabilistic participation, (B) additive noise, (C) temporal jitter, and (D) sequence warping. For each noise type, we show: (top)
examples of synthetic data at three different noise levels; (middle) similarity of extracted factors to ground-truth patterns across a range of noise levels
(20 fits for each level); and (bottom) examples of extracted factors W’s for one of the ground-truth patterns. Examples are shown at the same three
noise levels illustrated in the top row. In these examples, the algorithm was run with K ¼ 20, L ¼ 50 and l = 2l0 (via the procedure described in
Figure 4). For C, jitter displacements were draw from a discrete guassian distribution with the standard deviation in timesteps shown above For D,
timewarp conditions 1–10 indicate: 0, 66, 133, 200, 266, 333, 400, 466, 533 and 600 max % stretching respectively. For results at different values of l, see
Figure 3—figure supplement 1.
DOI: https://doi.org/10.7554/eLife.38471.009
The following figure supplements are available for figure 3:
Figure supplement 1. Robustness to noise at different values of l.
DOI: https://doi.org/10.7554/eLife.38471.010
Figure supplement 2. Robustness to small dataset size when using the x-ortho penalty.
DOI: https://doi.org/10.7554/eLife.38471.011
Figure supplement 3. Robustness to different types of sequences.
DOI: https://doi.org/10.7554/eLife.38471.012
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 10 of 42
Similarity to ground-truth
Tools and resources Neuroscience
similar to the case with no shared neurons (Figure 2E,F). We also examined other types of complex
structure and have found that the x-ortho penalty performs well in data with large gaps between
activity or with large overlaps of activity between neurons in the sequence. This approach also
worked well in cases in which the duration of the activity or the interval between the activity of neu-
rons varied across the sequence (Figure 3—figure supplement 3).
Robustness to noisy data
The cross-orthogonality penalty performed well in the presence of types of noise commonly found in
neural data. In particular, we considered: participation noise, in which individual neurons participate
probabilistically in instances of a sequence; additive noise, in which neuronal events occur randomly
outside of normal sequence patterns; temporal jitter, in which the timing of individual neurons is
shifted relative to their typical time in a sequence; and finally, temporal warping, in which each
instance of the sequence occurs at a different randomly selected speed. To test the robustness of
the algorithm with the x-ortho penalty to each of these noise conditions, we factorized data contain-
ing three neural sequences at a variety of noise levels (Figure 3, top row). The value of l was chosen
using methods described in the next section. Factorizations with the x-ortho penalty proved rela-
tively robust to all four noise types, with a high probability of returning the correct numbers of signif-
icant factors (Figure 4—figure supplement 5). Furthermore, under low-noise conditions, the
algorithm produced factors that were highly similar to ground-truth, and this similarity declined
gracefully at higher noise levels (Figure 3). Visualization of the extracted factors revealed a good
qualitative match to ground-truth sequences even in the presence of high noise except for the case
of temporal jitter (Figure 3). We also found that the x-ortho penalty allows reliable extraction of
sequences in which the duration of each neuron’s activity exhibits substantial random variation
across different renditions of the sequence, and in which the temporal jitter of neural activity exhibits
systematic variation at different points in the sequences (Figure 3—figure supplement 3).
Finally, we wondered how our approach with the x-ortho penalty performs on datasets with only
a small number of instances of each sequence. We generated data containing different numbers of
repetitions ranging from 1 to 20, of each underlying ground-truth sequence. For intermediate levels
of additive noise, we found that three repetitions of each sequence were sufficient to correctly
extract factors with similarity scores close to those obtained with much larger numbers of repetitions
(Figure 3—figure supplement 2).
Methods for choosing an appropriate value of l
The x-ortho penalty performs best when the strength of the regularization term (determined by the
hyperparameter l) is chosen appropriately. For l too small, the behavior of the algorithm
approaches that of convNMF, producing a large number of redundant factors with high x-ortho
cost. For l too large, all but one of the factors are suppressed to zero amplitude, resulting in a fac-
torization with near-zero x-ortho cost, but with large reconstruction error if multiple sequences are
present in the data. Between these extremes, there exists a region in which increasing l produces a
rapidly increasing reconstruction error and a rapidly decreasing x-ortho cost. Thus, there is a single
point, which we term l0, at which changes in reconstruction cost and changes in x-ortho cost are bal-
anced (Figure 4A). We hypothesized that the optimal choice of l (i.e. the one producing the correct
number of ground-truth factors) would lie near this point.
To test this intuition, we examined the performance of the x-ortho penalty as a function ofl in
noisy synthetic data consisting of three non-overlapping sequences (Figure 4A). Our analysis
revealed that, overall, values of l between 2l0 and 5l0 performed well for these data across all noise
types and levels (Figure 4B,C). In general, near-optimal performance was observed over an order of
magnitude range of l (Figure 1). However, there were systematic variations depending on noise
type: for additive noise, performance was better when l was closer to l0, while with other noise
types, performance was better at somewhat higher values of ls ( » 10l0).
Similar ranges of l appeared to work for datasets with different numbers of ground-truth sequen-
ces—for the datasets used in Figure 2D, a range of l between 0.001 and 0.01 returned the correct
number of sequences at least 90% of the time for datasets containing between 1 and 10 sequences
(Figure 2—figure supplement 2). Furthermore, this method for choosing l also worked on datasets
containing sequences with shared neurons (Figure 4—figure supplement 2).
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 11 of 42
Tools and resources Neuroscience
A G
B H
1
1
0.8 0.8
0.6 0.6
0.4 0.4
Reconstruction cost
0.2 X-ortho cost 0.2
C 0 I 0
1
10 110
0
10
D 1 J
0.5
0
E 1 K
0.5
0
F L
0.5
0
-5 0 5 -5 0 5
10 10 10 10 10 10
4
M 10
increasing noise No noise
3 Participation
10
* Additive
Jitter
2 Warping
10
1
10
0
10  = 0
-1
10
Figure 4. Procedure for choosing l for a new dataset based on finding a balance between reconstruction cost and
x-ortho cost. (A) Simulated data containing three sequences in the presence of participation noise (50%
participation probability). This noise condition is used for the tests in (B–F). (B) Normalized reconstruction cost
>
(jj 2Xe XjjF ) and cross-orthogonality cost (jjðW
 >XÞSH jj ;i 6¼j) as a function of l for 20 fits of these data. The cross-1
over point l0 is marked with a black circle. Note that in this plot the reconstruction cost and cross-orthogonality
Figure 4 continued on next page
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 12 of 42
 / 
0
Composite performance Similarity Fraction correct # Significant factors Cost
Tools and resources Neuroscience
Figure 4 continued
cost are normalized to vary between 0 and 1. (C) The number of significant factors obtained as a function of l; 20
fits, mean plotted in orange. Red arrow at left indicates the correct number of sequences (three). (D) Fraction of
fits returning the correct number of significant factors as a function of l. (E) Similarity of extracted factors to
ground-truth sequences as a function of l. (F) Composite performance, as the product of the curves in (D) and (E)
(smoothed using a three sample boxcar, plotted in orange with a circle marking the peak). Shaded region
indicates the range of l that works well ( half height of composite performance). (G–L) same as (A–F) but for
simulated data containing three noiseless sequences. (M) Summary plot showing the range of values of l (vertical
bars), relative to the cross-over point l0, that work well for each noise condition ( half height points of composite
performance). Circles indicate the value of l at the peak of the smoothed composite performance. For each noise
type, results for all noise levels from Figure 3 are shown (increasing color saturation at high noise levels; Green,
participation: 90, 80, 70, 60, 50, 40, 30, and 20%; Orange, additive noise 0.5, 1, 2, 2.5, 3, 3.5, and 4%; Purple, jitter:
SD of the distribution of random jitter: 5, 10, 15, 20, 25, 30, 35, 40, and 45 timesteps; Grey, timewarp: 66, 133, 200,
266, 333, 400, 466, 533, 600, and 666 max % stretching. Asterisk (*) indicates the noise type and level used in
panels (A–F). Gray band indicates a range between 2l0 and 5l0, a range that tended to perform well across the
different noise conditions. In real data, it may be useful to explore a wider range of l.
DOI: https://doi.org/10.7554/eLife.38471.013
The following figure supplements are available for figure 4:
Figure supplement 1. Analysis of the best range of l.
DOI: https://doi.org/10.7554/eLife.38471.014
Figure supplement 2. Procedure for choosing l applied to data with shared neurons.
DOI: https://doi.org/10.7554/eLife.38471.015
Figure supplement 3. Using cross-validation on held-out (masked) data to choose l.
DOI: https://doi.org/10.7554/eLife.38471.016
Figure supplement 4. Quantifying the effect of L1 sparsity penalties on W and H.
DOI: https://doi.org/10.7554/eLife.38471.017
Figure supplement 5. Comparing the performance of convNMF with an x-ortho or a sparsity penalty.
DOI: https://doi.org/10.7554/eLife.38471.018
The value of l may also be determined by cross-validation (see Materials and methods). Indeed,
the l chosen with the heuristic described above coincided with a minimum or distinctive feature in
the cross-validated test error for all the cases we examined (Figure 4—figure supplement 3). The
seqNMF code package accompanying this paper provides functions to determine l both by cross-
validation or in reference to l0.
Sparsity constraints to reduce redundant factors
One of the advantages of the x-ortho penalty is that it includes only a single term to penalize corre-
lations between different factors, and thus requires only a single hyperparameter l. This contrasts
with the approach of incorporating a sparsity constraint on W and H of the form
lwkWk þ lhkHk (Peter et al., 2017). We have found that the performance of the sparsity1 1
approach depends on the correct choice of both hyperparameters lw and lh (Figure 4—figure sup-
plement 4). Given the optimal choice of these parameters, the L1 sparsity constraint yields an overall
performance approximately as good as the x-ortho penalty (Figure 4—figure supplement 4). How-
ever, there are some consistent differences in the performance of the sparsity and x-ortho
approaches depending on noise type; an analysis at moderately high noise levels reveals that the
x-ortho penalty performs slightly better with warping and participation noise, while the L1 sparsity
penalty performs slightly better on data with jitter and additive noise (Figure 4—figure supplement
5). However, given the added complexity of choosing two hyperparameters for L1 sparsity, we pre-
fer the x-ortho approach.
Direct selection of K to reduce redundant factors
An alternative strategy to minimizing redundant factorizations is to estimate the number of underly-
ing sequences and to select the appropriate value of K. An approach for choosing the number of
factors in regular NMF is to run the algorithm many times with different initial conditions, at different
values of K, and choose the case with the most consistent and uncorrelated factors. This strategy is
called stability NMF (Wu et al., 2016) and is similar to other stability-based metrics that have been
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 13 of 42
Tools and resources Neuroscience
used in clustering models (von Luxburg, 2010). The stability NMF score, diss, is measured between
two factorizations, 1 ¼ f 1F W ; 1H g and 2F ¼ f 2W ; 2H g, run from different initial conditions:
 !
K K
1 X X
1 2
dissðF ;F Þ ¼ 2K max Cjk  max Cjk
2K 1kK 1jK
j¼1 k¼1
where C is the cross-correlation matrix between the columns of the matrix 1W and the the columns
of the matrix 2W . Note that diss is low when there is a one-to-one mapping between factors in 1F
and 2F , which tends to occur at the correct K in NMF (Wu et al., 2016; Ubaru et al., 2017). NMF is
run many times and the diss metric is calculated for all unique pairs. The best value of K is chosen as
that which yields the lowest average diss metric.
To use this approach for convNMF, we needed to slightly modify the stability NMF diss metric.
Unlike in NMF, convNMF factors have a temporal degeneracy; that is, one can shift the elements of
hk by one time step while shifting the elements of wk by one step in the opposite direction with little
change to the model reconstruction. Thus, rather than computing correlations from the factor pat-
terns W or loadings H, we computed the diss metric using correlations between factor reconstruc-
tions (Xe k ¼ wk
 hk).
h i
Tr e TXi Xe j
Cij ¼
kXe eikFkXjkF
P
where Tr½
Š denotes the trace operator, Tr½MŠ ¼ iMii. That is, Cij measures the correlation
between the reconstruction of factor i in 1 and the reconstruction of factor j in 2F F . Here, as for sta-
bility NMF, the approach is to run convNMF many times with different numbers of factors (K) and
choose the K which minimizes the diss metric.
We evaluated the robustness of this approach in synthetic data with the four noise conditions
examined earlier. Synthetic data were constructed with three ground-truth sequences and 20
convNMF factorizations were carried out for each K ranging from 1 to 10. For each K the average
diss metric was computed over all 20 factorizations. In many cases, the average diss metric exhibited
a minimum at the ground-truth K (Figure 5—figure supplement 1). As shown below, this method
also appears to be useful for identifying the number of sequences in real neural data.
Not only does the diss metric identify factorizations that are highly similar to the ground truth
and have the correct number of underlying factors, it also yields factorizations that minimize recon-
struction error in held out data (Figure 5, Figure 5—figure supplement 2), as shown using the
same cross-validation procedure described above (Figure 5—figure supplement 2). For simulated
datasets with participation noise, additive noise, and temporal jitter, there is a clear minimum in the
test error at the K given by diss metric. In other cases, there is a distinguishing feature such as a kink
or a plateau in the test error at this K (Figure 5—figure supplement 2).
Strategies for dealing with ambiguous sequence structure
Some sequences can be interpreted in multiple ways, and these interpretations will correspond to
different factorizations. A common example arises when neurons are shared between different
sequences, as is shown in Figure 6A and B. In this case, there are two ensembles of neurons (1 and
2), that participate in two different types of events. In one event type, ensemble one is active alone,
while in the other event type, ensemble one is coactive with ensemble 2. There are two different rea-
sonable factorizations of these data. In one factorization, the two different ensembles are separated
into two different factors, while in the other factorization the two different event types are separated
into two different factors. We refer to these as ’parts-based’ and ’events-based’ respectively. Note
that these different factorizations may correspond to different intuitions about underlying mecha-
nisms. ‘Parts-based’ factorizations will be particularly useful for clustering neurons into ensembles,
and ‘events-based’ factorizations will be particularly useful for correlating neural events with
behavior.
Here, we show that the addition of penalties on either W or H correlations can be used to shape
the factorizations of convNMF, with or without the x-ortho penalty, to produce ‘parts-based’ or
‘events-based’ factorization. Without this additional control, factorizations may be either ‘parts-
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 14 of 42
Tools and resources Neuroscience
A B
0.6 0.6
0.4 0.4
0.2 0.2
0 0
2 4 6 8 10 2 4 6 8 10
K K
C D
0.6 0.6
0.4 0.4
0.2 0.2
0 0
2 4 6 8 10 2 4 6 8 10
K K
Figure 5. Direct selection of K using the diss metric, a measure of the dissimilarity between different
factorizations. Panels show the distribution of diss as a function of K for several different noise conditions. Lower
values of diss indicate greater consistency or stability of the factorizations, an indication of low factor redundancy.
(A) probabilistic participation (60%), (B) additive noise (2.5% bins), (C) timing jitter (SD = 20 bins), and (D) sequence
warping (max warping = 266%). For each noise type, we show: (top) examples of synthetic data; (bottom) the diss
metric for 20 fits of convNMF for K from 1 to 10; the black line shows the median of the diss metric and the dotted
red line shows the true number of factors.
DOI: https://doi.org/10.7554/eLife.38471.019
The following figure supplements are available for figure 5:
Figure supplement 1. Direct selection of K using the diss metric for all noise conditions.
DOI: https://doi.org/10.7554/eLife.38471.020
Figure supplement 2. Estimating the number of sequences in a dataset using cross-validation on randomly
masked held-out datapoints.
DOI: https://doi.org/10.7554/eLife.38471.021
based’, or ‘events-based’ depending on initial conditions and the structure of shared neurons activi-
ties. This approach works because, in ‘events-based’ factorization, the H’s are orthogonal (uncorre-
lated) while the W’s have high overlap; conversely, in the ‘parts-based’ factorization, the W’s are
orthogonal while the H’s are strongly correlated. Note that these correlations in W or H are
unavoidable in the presence of shared neurons and such correlations do not indicate a redundant
factorization. Update rules to implement penalties on correlations in W or H are provided in Table 2
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 15 of 42
Dissimilarity Dissimilarity
Dissimilarity Dissimilarity
Tools and resources Neuroscience
with derivations in Appendix 1. Figure 9—figure supplement 2 shows examples of using these pen-
alties on the songbird dataset described in Figure 9.
L1 regularization is a widely used strategy for achieving sparse model parameters (Zhang et al.,
2016), and has been incorporated into convNMF in the past (O’Grady and Pearlmutter, 2006;
Ramanarayanan et al., 2013). In some of our datasets, we found it useful to include L1 regulariza-
tion for sparsity. The multiplicative update rules in the presence of L1 regularization are included in
Table 2, and as part of our code package. Sparsity on the matrices W and H may be particularly
useful in cases when sequences are repeated rhythmically (Figure 6—figure supplement 1A). For
example, the addition of a sparsity regularizer on the W update will bias the W exemplars to
include only a single repetition of the repeated sequence, while the addition of a sparsity regularizer
on H will bias the W exemplars to include multiple repetitions of the repeated sequence. Like the
ambiguities described above, these are both valid interpretations of the data, but each may be
more useful in different contexts.
Quantifying the prevalence of sequential structure in a dataset
While sequences may be found in a variety of neural datasets, their importance and prevalence is still
a matter of debate and investigation. To address this, we developed a metric to assess how much of
the explanatory power of a seqNMF factorization was due to synchronous vs. asynchronous neural
firing events. Since convNMF can fit both synchronous and sequential events in a dataset, recon-
struction error is not, by itself, diagnostic of the ‘sequenciness’ of neural activity. Our approach is
guided by the observation that in a data matrix with only synchronous temporal structure (i.e. pat-
terns of rank 1), the columns can be permuted without sacrificing convNMF reconstruction error. In
contrast, permuting the columns eliminates the ability of convNMF to model data that contains
Events-based Events-based
A C
Parts-based Parts-based
B D
Figure 6. Using penalties to bias toward events-based and parts-based factorizations. Datasets that have neurons
shared between multiple sequences can be factorized in different ways, emphasizing discrete temporal events
(events-based) or component neuronal ensembles (parts-based), by using orthogonality penalties on H or W to
penalize factor correlations (see Table 2). (Left) A dataset with two different ensembles of neurons that participate
in two different types of events, with (A) events-based factorization obtained using an orthogonality penalty on H
and (B) parts-based factorizations obtained using an orthogonality penalty on W. (Right) A dataset with six
different ensembles of neurons that participate in three different types of events, with (C) events-based and (D)
parts-based factorizations obtained as in (A) and (B).
DOI: https://doi.org/10.7554/eLife.38471.022
The following figure supplement is available for figure 6:
Figure supplement 1. Biasing factorizations between sparsity in W or H.
DOI: https://doi.org/10.7554/eLife.38471.023
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 16 of 42
Tools and resources Neuroscience
A
Only synchronous events 50% sequence (vs sync.) events Only sequence events
TIme
B
100 C 1.2
1
80 0.8
0.6
60
0.4
0.2
40
0
20 -0.2
0 0.5 1 0 0.5 1
Prob. sequence (vs sync) event Prob. sequence (vs sync) event
Figure 7. Using seqNMF to assess the prevalence of sequences in noisy data. (A) Example simulated datasets. Each dataset contains 10 neurons, with
varying amounts of additive noise, and varying proportions of synchronous events versus asynchronous sequences. For the purposes of this figure,
’sequence’ refers to a sequential pattern with no synchrony between different neurons in the pattern. The duration of each dataset used below is 3000
times, and here 300 timebins are shown. (B) Median percent power explained by convNMF (L = 12; K = 2; l=0) for each type of dataset (100 examples
of each dataset type). Different colors indicate the three different levels of additive noise shown in A. Solid lines and filled circles indicate results on
unshuffled datasets. Note that performance is flat for each noise level, regardless of the probability of sequences vs synchronous events. Dotted lines
and open circles indicate results on column-shuffled datasets. When no sequences are present, convNMF performs the same on column-shuffled data.
However, when sequences are present, convNMF performs worse on column-shuffled data. (C) For datasets with patterns ranging from exclusively
synchronous events to exclusively asynchronous sequences, convNMF was used to generate a ‘Sequenciness’ score. Colors correspond to different
noise levels shown in A. Asterisks denote cases where the power explained exceeds the Bonferroni-corrected significance threshold generated from
column-shuffled datasets. Open circles denote cases that do not achieve significance. Note that this significance test is fairly sensitive, detecting even
relatively low presence of sequences, and that the ‘Sequenciness’ score distinguishes between cases where more or less of the dataset consists of
sequences.
DOI: https://doi.org/10.7554/eLife.38471.024
sparse temporal sequences (i.e. high rank patterns) but no synchronous structure. We thus compute
a ‘sequenciness’ metric, ranging from 0 to 1, that compares the performance of convNMF on col-
umn-shuffled versus non-shuffled data matrices (see Materials and methods), and quantify the per-
formance of this metric in simulated datasets containing synchronous and sequential events with
varying prevalence (Figure 7C). We found that this metric varies approximately linearly with the
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 17 of 42
Noise
Neurons
% power explained
‘Sequenciness’
Tools and resources Neuroscience
degree to which sequences are present in a dataset. Below, we apply this method to real experimen-
tal data and obtain high ‘sequenciness’ scores, suggesting that convolutional matrix factorization is a
well-suited tool for summarizing neural dynamics in these datasets.
Application of seqNMF to hippocampal sequences
To test the ability of seqNMF to discover patterns in electrophysiological data, we analyzed multi-
electrode recordings from rat hippocampus (https://crcns.org/data-sets/hc), which were previously
shown to contain sequential patterns of neural firing (Pastalkova et al., 2015). Specifically, rats were
trained to alternate between left and right turns in a T-maze to earn a water reward. Between alter-
nations, the rats ran on a running wheel during an imposed delay period lasting either 10 or 20 sec-
onds. By averaging spiking activity during the delay period, the authors reported long temporal
sequences of neural activity spanning the delay. In some rats, the same sequence occurred on left
and right trials, while in other rats, different sequences were active in the delay period during each
trial types.
Without reference to the behavioral landmarks, seqNMF was able to extract sequences in both
datasets. In Rat 1, seqNMF extracted a single factor, corresponding to a sequence active throughout
the running wheel delay period and immediately after, when the rat ran up the stem of the maze
(Figure 8A); for 10 fits of K ranging from 1 to 10, the average diss metric reached a minimum at 1
and with l ¼ 2l0, most runs using the x-ortho penalty extracted a single significant factor
(Figure 8C–E). Factorizations of thes data with one factor captured 40.8% of the power in the data-
set on average, and had a ‘sequenciness’ score of 0.49. Some runs using the x-ortho penalty
extracted two factors (Figure 8E), splitting the delay period sequence and the maze stem sequence;
this is a reasonable interpretation of the data, and likely results from variability in the relative timing
of running wheel and maze stem traversal. At somewhat lower values of l, factorizations more often
split these sequences into two factors. At even lower values of l, factorizations had even more signif-
icant factors. Such higher granularity factorizations may correspond to real variants of the sequences,
as they generalize to held-out data or may reflect time warping in the data (Figure 5—figure sup-
plement 2J). However, a single sequence may be a better description of the data because the diss
metric displayed a clear minimum at K ¼ 1 (Figure 8C). In Rat 2, seqNMF typically identified three
factors (Figure 8B). The first two correspond to distinct sequences active for the duration of the
delay period on alternating left and right trials. A third sequence was active immediately following
each of the alternating sequences, corresponding to the time at which the animal exits the wheel
and runs up the stem of the maze. For 10 fits of K ranging from 1 to 10, the average diss metric
reached a minimum at three and with l ¼ 1:5l0, most runs with the x-ortho penalty extracted
between 2 and 4 factors (Figure 8F–H). Factorizations of these data with three factors captured
52.6% of the power in the dataset on average, and had a pattern ‘sequenciness’ score of 0.85. Taken
together, these results suggest that seqNMF can detect multiple neural sequences without the use
of behavioral landmarks.
Application of seqNMF to abnormal sequence development in avian
motor cortex
We applied seqNMF methods to analyze functional calcium imaging data recorded in the songbird
premotor cortical nucleus HVC during singing. Normal adult birds sing a highly stereotyped song,
making it possible to detect sequences by averaging neural activity aligned to the song. Using this
approach, it has been shown that HVC neurons generate precisely timed sequences that tile each
song syllable (Hahnloser et al., 2002; Picardo et al., 2016; Lynch et al., 2016). Songbirds learn
their song by imitation and must hear a tutor to develop normal adult vocalizations. Birds isolated
from a tutor sing highly variable and abnormal songs as adults (Fehér et al., 2009). Such ‘isolate’
birds provide an opportunity to study how the absence of normal auditory experience leads to path-
ological vocal/motor development. However, the high variability of pathological ‘isolate’ song makes
it difficult to identify neural sequences using the standard approach of aligning neural activity to
vocal output.
Using seqNMF, we were able to identify repeating neural sequences in isolate songbirds
(Figure 9A). At the chosen l (Figure 9B), x-ortho penalized factorizations typically extracted three
significant sequences (Figure 9C). Similarly, the diss measure has a local minimum at K ¼ 3
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 18 of 42
Tools and resources Neuroscience
reconstruction cost
h X-ortho cost
A 1 C Dw 0.51
1
0.4
0.5
0.3 0
E 10-530 s delay period lambda
0.2 0.6
0.4
h 0.1B 1w 0.21
0 0
1 3 5 7 9 0 2 4 6 8
K # significant factors
h2 F Gw2 0.5
1
0.4
0.5
0.3 0
-5
h H 103
w 0.6 lambda3 0.2
0.4
0.1
0.2
0 0
1 3 5 7 9 0 2 4 6 8
30 s K # significant factors
Figure 8. Application of seqNMF to extract hippocampal sequences from two rats. (A) Firing rates of 110 neurons recorded in the hippocampus of Rat
1 during an alternating left-right task with a delay period (Pastalkova et al., 2015). The single significant extracted x-ortho penalized factor. Both an
x-ortho penalized reconstruction of each factor (left) and raw data (right) are shown. Neurons are sorted according to the latency of their peak
activation within the factor. The red line shows the onset and offset of the forced delay periods, during which the animal ran on a treadmill. (B) Firing
rates of 43 hippocampal neurons recorded in Rat 2 during the same task (Mizuseki et al., 2013). Neurons are sorted according to the latency of their
peak activation within each of the three significant extracted sequences. The first two factors correspond to left and right trials, and the third
corresponds to running along the stem of the maze. (C) The diss metric as a function of K for Rat 1. Black line represents the median of the black
points. Notice the minimum at K = 1. (D) (Left) Reconstruction (red) and correlation (blue) costs as a function of l for Rat 1. Arrow indicates
l ¼ 58 10 , used for the x-ortho penalized factorization shown in (A). (E) Histogram of the number of significant factors across 30 runs of x-ortho
penalized convNMF. (D) The diss metric as a function of K for Rat 2. Notice the minimum at K = 3. (G–H) Same as in (D–E) but for Rat 2. Arrow indicates
l ¼ 8 510 , used for the factorization shown in (B).
DOI: https://doi.org/10.7554/eLife.38471.025
(Figure 9—figure supplement 1B). The three-sequence factorization explained 41% of the total
power in the dataset, with a sequenciness score of 0.7 andhe extracted sequences included sequen-
ces deployed during syllables of abnormally long and variable durations (Figure 9D–F, Figure 9—
figure supplement 1A).
In addition, the extracted sequences exhibit properties not observed in normal adult birds. We
see an example of two distinct sequences that sometimes, but not always, co-occur (Figure 9). We
observe that a shorter sequence (green) occurs alone on some syllable renditions while a second,
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 19 of 42
Neuron # (sorted)
Dissimilarity Dissimilarity
fraction of fits Cost fraction of fits Cost
Tools and resources Neuroscience
A B
Reconstruction cost
X-ortho cost
10-4 10-3 10-2 10-1 100
Lambda
C
50
40
30
20
1 sec
10
0
0 2 4 6 8 10
# significant factors
D E
w1 w2 w3
h1
F h2
h3
w1 w2 w3
Figure 9. SeqNMF applied to calcium imaging data from a singing isolate bird reveals abnormal sequence deployment. (A) Functional calcium signals
recorded from 75 neurons, unsorted, in a singing isolate bird. (B) Reconstruction and cross-orthogonality cost as a function of l. The arrow at l ¼ 0:005
indicates the value selected for the rest of the analysis. (C) Number of significant factors for 100 runs with the x-ortho penalty with K ¼ 10, l ¼ 0:005.
Arrow indicates three is the most common number of significant factors. (D) X-ortho factor exemplars (W’s). Neurons are grouped according to the
Figure 9 continued on next page
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 20 of 42
Neuron # (unsorted)
Neuron # (sorted) Neuron # (sorted)
% seqNMF runs
Cost (au)
Tools and resources Neuroscience
Figure 9 continued
factor in which they have peak activation, and within each group neurons are sorted by the latency of their peak activation within the factor. (E) The
same data shown in (A), after sorting neurons by their latency within each factor as in (D). A spectrogram of the bird’s song is shown at top, with a
purple ‘*’ denoting syllable variants correlated with w2. (F) Same as (E), but showing reconstructed data rather than calcium signals. Shown at top are
the temporal loadings (H) of each factor.
DOI: https://doi.org/10.7554/eLife.38471.026
The following figure supplements are available for figure 9:
Figure supplement 1. Further analysis of sequences.
DOI: https://doi.org/10.7554/eLife.38471.027
Figure supplement 2. Events-based and parts-based factorizations of songbird data.
DOI: https://doi.org/10.7554/eLife.38471.028
longer sequence (purple) occurs simultaneously on other syllable renditions. We found that biasing
x-ortho penalized convNMF towards ’parts-based’ or ’events-based’ factorizations gives a useful
tool to visualize this feature of the data (Figure 9—figure supplement 2). This probabilistic overlap
of different sequences is highly atypical in normal adult birds (Hahnloser et al., 2002; Long et al.,
2010; Picardo et al., 2016; Lynch et al., 2016) and is associated with abnormal variations in syllable
structure—in this case resulting in a longer variant of the syllable when both sequences co-occur.
This acoustic variation is a characteristic pathology of isolate song (Fehér et al., 2009).
Thus, even though we observe HVC generating sequences in the absence of a tutor, it appears
that these sequences are deployed in a highly abnormal fashion.
Application of seqNMF to a behavioral dataset: song spectrograms
Although we have focused on the application of seqNMF to neural activity data, these methods nat-
urally extend to other types of high-dimensional datasets, including behavioral data with applications
to neuroscience. The neural mechanisms underlying song production and learning in songbirds is an
area of active research. However, the identification and labeling of song syllables in acoustic record-
ings is challenging, particularly in young birds in which song syllables are highly variable. Because
automatic segmentation and clustering often fail, song syllables are still routinely labelled by hand
(Okubo et al., 2015). We tested whether seqNMF, applied to a spectrographic representation of
zebra finch vocalizations, is able to extract meaningful features in behavioral data. Using the x-ortho
penalty, factorizations correctly identified repeated acoustic patterns in juvenile songs, placing each
distinct syllable type into a different factor (Figure 10). The resulting classifications agree with previ-
ously published hand-labeled syllable types (Okubo et al., 2015). A similar approach could be
applied to other behavioral data, for example movement data or human speech, and could facilitate
the study of neural mechanisms underlying even earlier and more variable stages of learning. Indeed,
convNMF was originally developed for application to spectrograms (Smaragdis, 2004); notably it
has been suggested that auditory cortex may use similar computations to represent and parse natu-
ral statistics (Młynarski and McDermott, 2018).
Discussion
As neuroscientists strive to record larger datasets, there is a need for rigorous tools to reveal under-
lying structure in high-dimensional data (Gao and Ganguli, 2015; Sejnowski et al., 2014;
Churchland and Abbott, 2016; Bzdok and Yeo, 2017). In particular, sequential structure is increas-
ingly regarded as a fundamental property of neuronal circuits (Hahnloser et al., 2002;
Harvey et al., 2012; Okubo et al., 2015; Pastalkova et al., 2008), but standardized statistical
approaches for extracting such structure have not been widely adopted or agreed upon. Extracting
sequences is particularly challenging when animal behaviors are variable (e.g. during learning) or
absent entirely (e.g. during sleep).
Here, we explored a simple matrix factorization-based approach to identify neural sequences
without reference to animal behavior. The convNMF model elegantly captures sequential structure
in an unsupervised manner (Smaragdis, 2004; Smaragdis, 2007; Peter et al., 2017). However, in
datasets where the number of sequences is not known, convNMF may return inefficient and inconsis-
tent factorizations. To address these challenges, we introduced a new regularization term to penalize
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 21 of 42
Tools and resources Neuroscience
a b c a b c a b c a b
A
B
X-ortho cost
Reconstruction cost
10-5 10-4 10-3 10-2
Lambda
C
D
100 ms
Figure 10. SeqNMF applied to song spectrograms. (A) Spectrogram of juvenile song, with hand-labeled syllable
types (Okubo et al., 2015). (B) Reconstruction cost and x-ortho cost for these data as a function of l. Arrow
denotes l ¼ 0:0003, which was used to run convNMF with the x-ortho penalty (C) W’s for this song, fit with K ¼ 8,
L ¼ 200ms, l ¼ 0:0003. Note that there are three non-empty factors, corresponding to the three hand-labeled
syllables a, b, and c. (D) X-ortho penalized H’s (for the three non-empty factors) and reconstruction of the song
shown in (A) using these factors.
DOI: https://doi.org/10.7554/eLife.38471.029
correlated factorizations, and developed a new dissimilarity measure to assess model stability. Both
proposed methods can be used to infer the number of sequences in neural data and are highly
robust to noise. For example, even when (synthetic) neurons participate probabilistically in sequen-
ces at a rate of 50%, the model typically identifies factors with greater than 80% similarity to the
ground truth (Figure 3A). Additionally, these methods perform well even with very limited amounts
of data: for example successfully extracting sequences that only appear a handful of times in a noisy
data stream (Figure 3—figure supplement 2).
The x-ortho penalty developed in this paper may represent a useful improvement over traditional
orthogonality regularizations or suggest how traditional regularization penalties may be usefully
modified. First, it simultaneously provides a penalty on correlations in both W and H, thus simplify-
ing analyses by having only one penalty term. Second, although the x-ortho penalty does not for-
mally constitute regularization due to its inclusion of the data X, we have described how the penalty
can be approximated by a data-free regularization with potentially useful properties (Appendix 2).
Specifically, the data-free regularization contains terms corresponding to weighted orthogonality in
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 22 of 42
Cost (au)
Tools and resources Neuroscience
(smoothed) H and W, where the weights focus the orthogonality penalty preferentially on those fac-
tors that contribute the most power to the reconstruction. This concept of using power-weighted
regularization penalties may be applicable more generally to matrix factorization techniques.
As in many data analysis scenarios, a variety of statistical approaches may be brought to bear on
finding sequences in neural data. A classic method is to construct cross-correlogram plots, showing
spike time correlations between pairs of neurons at various time lags. However, other forms of spike
rate covariation, such as trial-to-trial gain modulation, can produce spurious peaks in this measure
(Brody, 1999); recent work has developed statistical corrections for these effects (Russo and Dur-
stewitz, 2017). After significant pairwise correlations are identified, one can heuristically piece
together pairs of neurons with significant interactions into a sequence. This bottom-up approach
may be better than seqNMF at detecting sequences involving small numbers of neurons, since
seqNMF specifically targets sequences that explain large amounts of variance in the data. On the
other hand, bottom-up approaches to sequence extraction may fail to identify long sequences with
high participation noise or jitter in each neuron (Quaglio et al., 2018). One can think of seqNMF as
a complementary top-down approach, which performs very well in the high-noise regime since it
learns a template sequence at the level of the full population that is robust to noise at the level of
individual units.
Statistical models with a dynamical component, such as Hidden Markov Models (HMMs)
(Maboudi et al., 2018), linear dynamical systems (Kao et al., 2015), and models with switching
dynamics (Linderman et al., 2017), can also capture sequential firing patterns. These methods will
typically require many hidden states or latent dimensions to capture sequences, similar to PCA and
NMF which require many components to recover sequences. Nevertheless, visualizing the transition
matrix of an HMM can provide insight into the order in which hidden states of the model are visited,
mapping onto different sequences that manifest in population activity (Maboudi et al., 2018). One
advantage of this approach is that it can model sequences that occasionally end prematurely, while
convNMF will always reconstruct the full sequence. On the other hand, this pattern completion prop-
erty makes convNMF robust to participation noise and jitter. In contrast, a standard HMM must pass
through each hidden state to model a sequence, and therefore may have trouble if many of these
states are skipped. Thus, we expect HMMs and related models to exhibit complementary strengths
and weaknesses when compared to convNMF.
Another strength of convNMF is its ability to accommodate sequences with shared neurons, as
has been observed during song learning (Okubo et al., 2015). Sequences with shared neurons can
be interpreted either in terms of ‘parts-based’ or ‘events-based’ factorizations (Figure 9—figure
supplement 2). This capacity for a combinatorial description of overlapping sequences distinguishes
convNMF from many other methods, which assume that neural patterns/sequences do not co-occur
in time. For example, a vanilla HMM can only model each time step with a single hidden state and
thus cannot express parts-based representations of neural sequences. Likewise, simple clustering
models would assign each time interval to a single cluster label. Adding hierarchical and factorial
structure to these models could allow them to test for overlapping neural sequences (see e.g.
Ghahramani and Jordan, 1997); however, we believe seqNMF provides a simpler and more direct
framework to explore this possibility.
Finally, as demonstrated by our development of new regularization terms and stability measures,
convolutional matrix factorization is a flexible and extensible framework for sequence extraction. For
example, one can tune the overall sparsity in the model by introducing additional L1 regularization
terms. The loss function may also be modified, for example substituting in KL divergence or more
general b-divergence (Villasana et al., 2018). Both L1 regularization and b-divergence losses are
included in the seqNMF code package so that the model can be tuned to the particular needs of
future analyses. Future development could incorporate outlier detection into the objective function
(Netrapalli et al., 2014), or online optimization methods for large datasets (Wang et al., 2013).
Other extensions to NMF, for example, Union of Intersections NMF Cluster (Ubaru et al., 2017),
have yielded increased robustness and consistency of NMF factorizations, and could potentially also
be modified for application to convNMF. Thus, adding convolutional structure to factorization-based
models of neural data represents a rich opportunity for statistical neuroscience.
Despite limiting ourselves to a relatively simple model for the purposes of this paper, we
extracted biological insights that would have been difficult to otherwise achieve. For example, we
identified neural sequences in isolated songbirds without aligning to song syllables, enabling new
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 23 of 42
Tools and resources Neuroscience
research into songbird learning on two fronts. First, since isolated and juvenile birds sing highly vari-
able songs that are not easily segmented into stereotyped syllables, it is difficult and highly subjec-
tive to identify sequences by aligning to human-labeled syllables. SeqNMF enables the discovery
and future characterization of neural sequences in these cases. Second, while behaviorally aligned
sequences exist in tutored birds, it is that possible many neural sequences—for example, in different
brain areas or stages of development—are not closely locked to song syllables. Thus, even in cases
where stereotyped song syllables exist, behavioral alignment may overlook relevant sequences and
structure in the data. These lessons apply broadly to many neural systems, and demonstrate the
importance of general-purpose methods that extract sequences without reference to behavior. Our
results show that convolutional matrix factorization models are an attractive framework to meet this
need.
Materials and methods
Key resources table
Reagent type Source or Additional
(species) or resource Designation reference Identifiers information
Software, seqNMF this paper https://github.com/FeeLab/seqNMF start with demo.m
algorithm
Software, convNMF Smaragdis, 2004; https://github.com/colinvaz/nmf-toolbox
algorithm Smaragdis, 2007
Software, sparse convNMF O’Grady and Pearlmutter, 2006; https://github.com/colinvaz/nmf-toolbox
algorithm Ramanarayanan et al., 2013
Software, NMF orthogonality penalties Choi, 2008;
algorithm Chen and Cichocki, 2004
Software, other NMF extensions Cichocki et al., 2009
algorithm
Software, NMF Lee and Seung, 1999
algorithm
Software, CNMF_E (cell extraction) Zhou et al., 2018 https://github.com/zhoupc/CNMF_E
algorithm
Software, MATLAB MathWorks www.mathworks.com,
algorithm RRID:SCR_001622
Strain, strain AAV9.CAG.GCaMP6f. Chen et al., 2013 Addgene viral prep # 100836-AAV9,
background WPRE.SV40 http://n2t.net/addgene:100836,
(adeno-associated virus) RRID:Addgene_100836
Commercial Miniature Inscopix https://www.inscopix.com/nvista
assay or kit microscope
Contact for resource sharing
Further requests should be directed to Michale Fee (fee@mit.edu).
Software and data availability
The seqNMF MATLAB code is publicly available as a github repository, which also includes our song-
bird data (Figure 9) for demonstration (Mackevicius et al., 2018; copy archived at https://github.
com/elifesciences-publications/seqNMF).
The repository includes the seqNMF function, as well as helper functions for selecting l, testing
the significance of factors, plotting, and other functions. It also includes a demo script with an exam-
ple of how to select l for a new dataset, test for significance of factors, plot the seqNMF factoriza-
tion, switch between parts-based and events-based factorizations, and calculate cross-validated
performance on a masked test set.
Generating simulated data
We simulated neural sequences containing between 1 and 10 distinct neural sequences in the pres-
ence of various noise conditions. Each neural sequence was made up of 10 consecutively active
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 24 of 42
Tools and resources Neuroscience
neurons, each separated by three timebins. The binary activity matrix was convolved with an expo-
nential kernel (t ¼ 10 timebins) to resemble neural calcium imaging activity.
SeqNMF algorithm details
The x-ortho penalized convNMF algorithm is a direct extension of the multiplicative update
convNMF algorithm (Smaragdis, 2004), and draws on previous work regularizing NMF to encourage
factor orthogonality (Chen and Cichocki, 2004).
The uniqueness and consistency of traditional NMF has been better studied than convNMF. In
special cases, NMF has a unique solution comprised of sparse, ‘parts-based’ features that can be
consistently identified by known algorithms (Donoho and Stodden, 2004; Arora et al., 2011). How-
ever, this ideal scenario does not hold in many practical settings. In these cases, NMF is sensitive to
initialization, resulting in potentially inconsistent features. This problem can be addressed by intro-
ducing additional constraints or regularization terms that encourage the model to extract particular,
e.g. sparse or approximately orthogonal features (Huang et al., 2014; Kim and Park, 2008). Both
theoretical work and empirical observations suggest that these modifications result in more consis-
tently identified features (Theis et al., 2005; Kim and Park, 2008).
For x-ortho penalized seqNMF, we added to the convNMF cost function a term that promotes
competition between overlapping factors, resulting in the following cost function:
 
ðWf
>
;He Þ ¼ 2 >arg min jjXe XjjF þljjðW
 XÞSH jj1;i¼6 j (8)
W;H
We derived the following multiplicative update rules for W and H (Appendix 1):
 
‘! >
X H
W

‘ W

‘  > (9)
e ‘!  ‘ >X H þlXSH ð1 IÞ
>
W
 X
H H (10)
> >
W
 Xe þlð1 IÞðW
 XSÞ
‘!
where the division and  are element-wise. The operator ð
Þ shifts a matrix in the ! direction by ‘
 ‘
timebins, that is a delay by ‘ timebins, and ð
Þ shifts a matrix in the  direction by ‘ timebins (nota-
tion summary, Table 1). Note that multiplication with the KK matrix ð1 IÞ effectively implements
factor competition because it places in the kth row a sum across all other factors. These update rules
are derived in Appendix 1 by taking the derivative of the cost function in Equation 8 and choosing
an appropriate learning rate for each element.
In addition to the multiplicative updates outlined in Table 2, we also renormalize so rows of H
have unit norm; shift factors to be centered in time such that the center of mass of each W pattern
occurs in the middle; and in the final iteration run one additional step of unregularized convNMF to
prioritize the cost of reconstruction error over the regularization (Algorithm 1). This final step is done
to correct a minor suppression in the amplitude of some peaks in H that may occur within 2L time-
bins of neighboring sequences.
Testing the significance of each factor on held-out data
In order to test whether a factor is significantly present in held-out data, we measured the distribu-
tion across timebins of the overlaps of the factor with the held-out data, and compared the skewness
>
of this distribution to the null case (Figure 1). Overlap with the data is measured as W
 X,
a quantity which will be high at timepoints when the sequence occurs, producing a distribution of
>
W
 X with high skew. In contrast, a distribution of overlaps exhibiting low skew indicates a
sequence is not present in the data, since there are few timepoints of particularly high overlap. We
estimated what skew levels would appear by chance by constructing null factors where temporal
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 25 of 42
Tools and resources Neuroscience
relationships between neurons have been eliminated. To create such null factors, we start from the
real factors then circularly shift the timecourse of each neuron by a random amount between 0 and
L. We measure the skew of the overlap distributions for each null factor, and ask whether the skew
we measured for the real factor is significant at p-value a, that is, if it exceeds the Bonferroni cor-
rected ðð  aÞ  Þth1 100 percentile of the null skews (see Figure 2—figure supplement 1).
K
Algorithm 1: SeqNMF (x-ortho algorithm)
Input: Data matrix X, number of factors K, factor duration L, regularization
strength l
Output: Factor exemplars W, factor timecourses H
1 Initialize W and H randomly
2 Iter = 1
3 While (Iter < maxIter) and (D cost > tolerance) do
4 Update H using multiplicative update from Table 2
5 Shift W and H to center W’s in time
6 Renormalize W and H so rows of H have unit norm
7 Update W using multiplicative update from Table 2
8 Iter = Iter + 1
9 Do one final unregularized convNMF update of W and H
10 return
Note that if l is set too small, seqNMF will produce multiple redundant factors to explain one
sequence in the data. In this case, each redundant candidate sequence will pass the significance test
outlined here. We will address below a procedure for choosing l and methods for determining the
number of sequences.
Calculating the percent power explained by a factorization
In assessing the relevance of sequences in a dataset, it can be useful to calculate what percentage of
the total power in the dataset is explained by the factorization (Xe ). The total power in the data is
P P P P
2
X (abbreviating N Tn¼ t¼ x
2
nt to X
2). The power unexplained by the factorization is
P 1 1
ð  2X Xe Þ . Thus, the percent of the total power explained by the factorization is:
P P P 2
2
X P ðXX
e Þ2 2PXX
e Xe
¼ (11)
2 2
X X
‘Sequenciness’ score
The ‘sequenciness’ score was developed to distinguish between datasets with exclusively synchro-
nous patterns, and datasets with temporally extended sequential patterns. This score relies on the
observation that synchronous patterns are not disrupted by shuffling the columns of the data matrix.
The ‘sequenciness’ score is calculated by first computing the difference between the power
explained by seqNMF in the actual and column-shuffled data. This quantity is then divided by the
power explained in the actual data minus the power explained in data where each neuron is time-
shuffled by a different random permutation.
Choosing appropriate parameters for a new dataset
The choice of appropriate parameters (l, K and L) will depend on the data type (sequence length,
number, and density; amount of noise; etc.).
In practice, we found that results were relatively robust to the choice of parameters. When K or L
is set larger than necessary, seqNMF tends to simply leave the unnecessary factors or times empty.
For choosing l, the goal is to find the ‘sweet spot’ (Figure 4) to explain as much data as possible
while still producing sensible factorizations, that is, minimally correlated factors, with low values of
>
jjð >W
 XÞSH jj ;i 6¼j. Our software package includes demo code for determining the best parameters1
for a new type of data, using the following strategy:
1. Start with K slightly larger than the number of sequences anticipated in the data
2. Start with L slightly longer than the maximum expected factor length
3. Run seqNMF for a range of l’s, and for each l measure the reconstruction error
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 26 of 42
Tools and resources Neuroscience
   >
jj  
 jj2X W H and the correlation cost term jjðW
 >F XÞSH jj1;i¼6 j
4. Choose a l slightly above the crossover point l0
5. Decrease K if desired, as otherwise some factors will be consistently empty
6. Decrease L if desired, as otherwise some times will consistently be empty
In some applications, achieving the desired accuracy may depend on choosing a l that allows
some inconsistency. It is possible to deal with this remaining inconsistency by comparing factors pro-
duced by different random initializations, and only considering factors that arise from several differ-
ent initializations, a strategy that has been previously applied to standard convNMF on neural data
(Peter et al., 2017).
During validation of our procedure for choosing l, we compared factorizations to ground truth
sequences as shown in Figure 4. To find the optimal value of l, we used the product of two curves.
The first curve was obtained by calculating the fraction of fits in which the true number of sequences
was recovered as a function of l. The second curve was obtained by calculating similarity to ground
truth as a function of l (see Materials and methods section ‘Measuring performance on noisy fits by
comparing seqNMF sequence to ground-truth sequences’). The product of these two curves was
smoothed using a three-sample boxcar sliding window, and the width was found as the values of l
on either side of the peak value that correspond most closely to the half-maximum points of the
curve.
Preprocessing
While seqNMF is generally quite robust to noisy data, and different types of sequential patterns,
proper preprocessing of the data can be important to obtaining reasonable factorizations on real
neural data. A key principle is that, in minimizing the reconstruction error, seqNMF is most strongly
influenced by parts of the data that exhibit high variance. This can be problematic if the regions of
interest in the data have relatively low amplitude. For example, high firing rate neurons may be pri-
oritized over those with lower firing rate. As an alternative to subtracting the mean firing rate of
each neuron, which would introduce negative values, neurons could be normalized divisively or by
subtracting off a NMF reconstruction fit using a method that forces a non-negative residual
(Kim and Smaragdis, 2014). Additionally, variations in behavioral state may lead to seqNMF factori-
zations that prioritize regions of the data with high variance and neglect other regions. It may be
possible to mitigate these effects by normalizing data, or by restricting analysis to particular subsets
of the data, either by time or by neuron.
Measuring performance on noisy data by comparing seqNMF sequences to
ground-truth sequences
We wanted to measure the ability of seqNMF to recover ground-truth sequences even when the
sequences are obstructed by noise. Our noisy data consisted of three ground-truth sequences,
obstructed by a variety of noise types. For each ground-truth sequence, we found its best match
among the seqNMF factors. This was performed in a greedy manner. Specifically, we first computed
a reconstruction for one of the ground-truth factors. We then measured the correlation between this
reconstruction and reconstructions generated from each of the extracted factors, and chose the best
match (highest correlation). Next, we matched a second ground-truth sequence with its best match
(highest correlation between reconstructions) among the remaining seqNMF factors, and finally we
found the best match for the third ground-truth sequence. The mean of these three correlations was
used as a measure of similarity between the seqNMF factorization and the ground-truth (noiseless)
sequences.
Testing generalization of factorization to randomly held-out (masked) data
entries
The data matrix X was divided into training data and test data by randomly selecting 5 or 10% of
matrix entries to hold out. Specifically, the objective function (Equation 5, in the Results section) was
modified to:
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 27 of 42
Tools and resources Neuroscience
2
arg min jjMðW
 HXÞjjF þR (12)
W;H
where  indicates elementwise multiplication (Hadamard product) and M is a binary matrix with 5 or
10% of the entries randomly selected to be zero (held-out test set) and the remaining 95 or 90% set
to one (training set). To search for a solution, we reformulate this optimization problem as:
2
arg min jjW
 HZjjF þR
W;H;Z (13)
subjecttoMZ¼MX
where we have introduced a new optimization variable Z, which can be thought of as a surrogate
dataset that is equal to the ground truth data only on the training set. The goal is now to minimize
the difference between the model estimate, Xe ¼W
 H, and the surrogate, Z, while constraining Z
to equal X at unmasked elements (where mij ¼ 1) and allowing Z to be freely chosen at masked ele-
ments (where mij ¼ 0). Clearly, at masked elements, the best choice is to make Z equal to the current
model estimate Xe as this minimizes the cost function without violating the constraint. This leads to
the following update rules which are applied cyclically to update Z, W, and H.

Xnt if Mnt ¼ 1
Znt (14)
ðW
 HÞnt if Mnt ¼ 0
 
‘! >
Z H
W

‘ W

‘  
‘! >
(15)
Xe
 ‘
H þl >ZSH ð1 IÞ
>
W
 Z
H H (16)
> >
W
 Xe þlð1 IÞðW
 ZSÞ
The measure used for testing generalization performance was root mean squared error (RMSE).
For the testing phase, RMSE was computed from the difference between Xe and the data matrix X
only for held-out entries.
Hippocampus data
The hippocampal data was collected in the Buzsaki lab (Pastalkova et al., 2015; Mizuseki et al.,
2013), and is publicly available on the Collaborative Research in Computational Neuroscience
(CRCNS) Data sharing website. The dataset we refer to as ‘Rat 1’ is in the hc-5 dataset, and the data-
set we refer to as ‘Rat 2’ is in the hc-3 dataset. Before running seqNMF, we processed the data by
convolving the raw spike trains with a gaussian kernel of standard deviation 100 ms.
Animal care and use
We used male zebra finches (Taeniopygia guttata) from the MIT zebra finch breeding facility (Cam-
bridge, MA). Animal care and experiments were carried out in accordance with NIH guidelines, and
reviewed and approved by the Massachusetts Institute of Technology Committee on Animal Care
(protocol 0715-071-18).
In order to prevent exposure to a tutor song, birds were foster-raised by female birds, which do
not sing, starting on or before post-hatch day 15. For experiments, birds were housed singly in cus-
tom-made sound isolation chambers.
Data acquisition and preprocessing
The calcium indicator GCaMP6f was expressed in HVC by intracranial injection of the viral vector
AAV9.CAG.GCaMP6f.WPRE.SV40 (Chen et al., 2013) into HVC. In the same surgery, a cranial win-
dow was made using a GRIN (gradient index) lens (1 mm diamenter, 4 mm length, Inscopix). After at
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 28 of 42
Tools and resources Neuroscience
least one week, in order to allow for sufficient viral expression, recordings were made using the
Inscopix nVista miniature fluorescent microscope.
Neuronal activity traces were extracted from raw fluorescence movies using the CNMF_E algo-
rithm, a constrained non-negative matrix factorization algorithm specialized for microendoscope
data by including a local background model to remove activity from out-of-focus cells (Zhou et al.,
2018).
We performed several preprocessing steps before applying seqNMF to functional calcium traces
extracted by CNMF_E. First, we estimated burst times from the raw traces by deconvolving the
traces using an AR-2 process. The deconvolution parameters (time constants and noise floor) were
estimated for each neuron using the CNMF_E code package (Zhou et al., 2018). Some neurons
exhibited larger peaks than others, likely due to different expression levels of the calcium indicator.
Since seqNMF would prioritize the neurons with the most power, we renormalized by dividing the
signal from each neuron by the sum of the maximum value of that row and the th95 percentile of the
signal across all neurons. In this way, neurons with larger peaks were given some priority, but not
much more than that of neurons with weaker signals.
Acknowledgements
This work was supported by a grant from the Simons Collaboration for the Global Brain, the National
Institutes of Health (NIH) [grant number R01 DC009183] and the G Harold and Leila Y Mathers Chari-
table Foundation. ELM received support through the NDSEG Fellowship program. AHB received
support through NIH training grant 5T32EB019940-03. MSG received support from the NIH [grant
number U19NS104648]. AHW received support from the U.S. Department of Energy Computational
Science Graduate Fellowship (CSGF) program. Thanks to Pengcheng Zhou for advice on his
CNMF_E calcium data cell extraction algorithm. Thanks to Wiktor Młynarski for helpful convNMF dis-
cussions. Thanks to Michael Stetner, Galen Lynch, Nhat Le, Dezhe Jin, Edward Nieh, Adam Charles,
Jane Van Velden and Yiheng Wang for comments on the manuscript and on our code package.
Thanks to our reviewers for wonderful suggestions, including the use of diss to select K, and using
seqNMF to measure ‘sequenciness’. Special thanks to the 2017 Methods in Computational Neurosci-
ence course [supported by NIH grant R25 MH062204 and Simons Foundation] at the Woods Hole
Marine Biology Lab, where this collaboration was started.
Additional information
Funding
Funder Grant reference number Author
Simons Foundation Simons Collaboration for Mark S Goldman
the Global Brain Michale S Fee
National Institute on Deafness R01 DC009183 Michale S Fee
and Other Communication
Disorders
G Harold and Leila Y. Mathers Michale S Fee
Foundation
U.S. Department of Defense NDSEG Fellowship program Emily L Mackevicius
Department of Energy, Labor Computational Science Alex H Williams
and Economic Growth Graduate Fellowship (CSGF)
NIH Office of the Director Training Grant Andrew H Bahle
5T32EB019940-03
National Institute of Neurolo- U19 NS104648 Mark S Goldman
gical Disorders and Stroke
National Institute of Mental R25 MH062204 Mark S Goldman
Health Michale S Fee
The funders had no role in study design, data collection and interpretation, or the
decision to submit the work for publication.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 29 of 42
Tools and resources Neuroscience
Author contributions
Emily L Mackevicius, Andrew H Bahle, Conceptualization, Data curation, Software, Formal analysis,
Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and
editing; Alex H Williams, Conceptualization, Formal analysis, Validation, Methodology, Writing—
review and editing; Shijie Gu, Data curation; Natalia I Denisenko, Investigation; Mark S Goldman,
Conceptualization, Formal analysis, Supervision, Funding acquisition, Methodology, Writing—
original draft, Project administration, Writing—review and editing; Michale S Fee, Conceptualization,
Resources, Formal analysis, Supervision, Funding acquisition, Validation, Visualization, Methodology,
Writing—original draft, Project administration, Writing—review and editing
Author ORCIDs
Emily L Mackevicius http://orcid.org/0000-0001-6593-4398
Andrew H Bahle http://orcid.org/0000-0003-0567-7195
Alex H Williams https://orcid.org/0000-0001-5853-103X
Shijie Gu http://orcid.org/0000-0001-6257-5756
Mark S Goldman http://orcid.org/0000-0002-8257-2314
Michale S Fee http://orcid.org/0000-0001-7539-1745
Ethics
Animal experimentation: Animal care and experiments were carried out in accordance with NIH
guidelines, and reviewed and approved by the Massachusetts Institute of Technology Committee on
Animal Care (protocol 0715-071-18).
Decision letter and Author response
Decision letter https://doi.org/10.7554/eLife.38471.040
Author response https://doi.org/10.7554/eLife.38471.041
Additional files
Supplementary files
. Transparent reporting form
DOI: https://doi.org/10.7554/eLife.38471.031
Data availability
Code and songbird dataset is publicly available on github: https://github.com/FeeLab/seqNMF
(copy archived at https://github.com/elifesciences-publications/seqNMF). Rat datasets were col-
lected in the Buzsaki lab, and are publicly available on CRCNS (https://crcns.org/data-sets/hc); users
must first create a free account (https://crcns.org/register) before they can download the datasets
from the site.
The following previously published datasets were used:
Database and
Author(s) Year Dataset title Dataset URL Identifier
Mizuseki K, Sirota 2013 Multiple single unit recordings from https://crcns.org/data- Collaborative
A, Pastalkova E, different rat hippocampal and sets/hc/hc-3 Research in
Diba K entorhinal regions while the Computational
animals were performing multiple Neuroscience, 10.60
behavioral tasks. CRCNS.org. 80/K09G5JRZ
Pastalkova E, Wang 2015 Simultaneous extracellular https://crcns.org/data- Collaborative
Y recordings from left and right sets/hc/hc-5 Research in
hippocampal areas CA1 and right Computational
entorhinal cortex from a rat Neuroscience, 10.60
performing a left / right alternation 80/K0KS6PHF
task and other behaviors. CRCNS.
org.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 30 of 42
Tools and resources Neuroscience
References
Arora S, Ge R, Kannan R, Moitra A. 2011. Computing a nonnegative matrix factorization – provably. arXiv.
https://arxiv.org/abs/1111.0952.
Bapi RS, Pammi VSC, Miyapuram KP, Ahmed. 2005. Investigation of sequence processing: a cognitive and
computational neuroscience perspective. Current Science 89:1690–1698.
Bro R, Kjeldahl K, Smilde AK, Kiers HA. 2008. Cross-validation of component models: a critical look at current
methods. Analytical and Bioanalytical Chemistry 390:1241–1251. DOI: https://doi.org/10.1007/s00216-007-
1790-1, PMID: 18214448
Brody CD. 1999. Correlations without synchrony. Neural Computation 11:1537–1551. DOI: https://doi.org/10.
1162/089976699300016133, PMID: 10490937
Brunton BW, Johnson LA, Ojemann JG, Kutz JN. 2016. Extracting spatial-temporal coherent patterns in large-
scale neural recordings using dynamic mode decomposition. Journal of Neuroscience Methods 258:1–15.
DOI: https://doi.org/10.1016/j.jneumeth.2015.10.010, PMID: 26529367
Bzdok D, Yeo BTT. 2017. Inference in the age of big data: Future perspectives on neuroscience. NeuroImage
155:549–564. DOI: https://doi.org/10.1016/j.neuroimage.2017.04.061, PMID: 28456584
Chen TW, Wardill TJ, Sun Y, Pulver SR, Renninger SL, Baohan A, Schreiter ER, Kerr RA, Orger MB, Jayaraman V,
Looger LL, Svoboda K, Kim DS. 2013. Ultrasensitive fluorescent proteins for imaging neuronal activity. Nature
499:295–300. DOI: https://doi.org/10.1038/nature12354, PMID: 23868258
Chen Z, Cichocki A. 2004. Nonnegative matrix factorization with temporal smoothness and/or spatial
decorrelation constraints. Signal Processing 11.
Choi S. 2008. Algorithms for orthogonal nonnegative matrix factorization. IEEE International Joint Conference on
Neural Networks (IEEE World Congress on Computational Intelligence) 1828–1832.
Churchland AK, Abbott LF. 2016. Conceptual and technical advances define a key moment for theoretical
neuroscience. Nature Neuroscience 19:348–349. DOI: https://doi.org/10.1038/nn.4255, PMID: 26906500
Cichocki A, Zdunek R, Phan AH, Amari S-I. 2009. Nonnegative Matrix and Tensor Factorizations: Applications to
Exploratory Multi-Way Data Analysis and Blind Source Separation. Wiley. ISBN: 9780470747278 DOI: https://
doi.org/10.1002/9780470747278
Clegg BA, Digirolamo GJ, Keele SW. 1998. Sequence learning. Trends in Cognitive Sciences 2:275–281.
DOI: https://doi.org/10.1016/S1364-6613(98)01202-9, PMID: 21227209
Cui Y, Ahmad S, Hawkins J. 2016. Continuous online sequence learning with an unsupervised neural network
model. Neural Computation 28:2474–2504. DOI: https://doi.org/10.1162/NECO_a_00893
Donoho D, Stodden V. 2004. When does non-negative matrix factorization give a correct decomposition into
parts? In: Thrun S, Saul L. K, Schölkopf B (Eds). Advances in Neural Information Processing Systems. 16 MIT
Press. p. 1141–1148.
Fehér O, Wang H, Saar S, Mitra PP, Tchernichovski O. 2009. De novo establishment of wild-type song culture in
the zebra finch. Nature 459:564–568. DOI: https://doi.org/10.1038/nature07994, PMID: 19412161
Fujisawa S, Amarasingham A, Harrison MT, Buzsáki G. 2008. Behavior-dependent short-term assembly dynamics
in the medial prefrontal cortex. Nature Neuroscience 11:823–833. DOI: https://doi.org/10.1038/nn.2134,
PMID: 18516033
Gao P, Ganguli S. 2015. On simplicity and complexity in the brave new world of large-scale neuroscience.
Current Opinion in Neurobiology 32:148–155. DOI: https://doi.org/10.1016/j.conb.2015.04.003, PMID: 25932
978
Gerstein GL, Williams ER, Diesmann M, Grün S, Trengove C. 2012. Detecting synfire chains in parallel spike data.
Journal of Neuroscience Methods 206:54–64. DOI: https://doi.org/10.1016/j.jneumeth.2012.02.003,
PMID: 22361572
Ghahramani Z, Jordan MI. 1997. Factorial hidden markov models. Machine Learning 29:245–273. DOI: https://
doi.org/10.1023/A:1007425814087
Grossberger L, Battaglia FP, Vinck M. 2018. Unsupervised clustering of temporal patterns in high-dimensional
neuronal ensembles using a novel dissimilarity measure. PLOS Computational Biology 14:e1006283.
DOI: https://doi.org/10.1371/journal.pcbi.1006283, PMID: 29979681
Hahnloser RH, Kozhevnikov AA, Fee MS. 2002. An ultra-sparse code underlies the generation of neural
sequences in a songbird. Nature 419:65–70. DOI: https://doi.org/10.1038/nature00974, PMID: 12214232
Harvey CD, Coen P, Tank DW. 2012. Choice-specific sequences in parietal cortex during a virtual-navigation
decision task. Nature 484:62–68. DOI: https://doi.org/10.1038/nature10918, PMID: 22419153
Hastie T, Tibshirani R, Friedman JH. 2009. The Elements of Statistical Learning: Data Mining, Inference, and
Prediction. Springer. DOI: https://doi.org/10.1007/978-0-387-84858-7
Hawkins J, Ahmad S. 2016. Why neurons have thousands of synapses, a theory of sequence memory in
neocortex. Frontiers in Neural Circuits 10:23. DOI: https://doi.org/10.3389/fncir.2016.00023, PMID: 27065813
Huang K, Sidiropoulos ND, Swami A. 2014. Non-negative matrix factorization revisited: uniqueness and
algorithm for symmetric decomposition. IEEE Transactions on Signal Processing 62:211–224. DOI: https://doi.
org/10.1109/TSP.2013.2285514
Janata P, Grafton ST. 2003. Swinging in the brain: shared neural substrates for behaviors related to sequencing
and music. Nature Neuroscience 6:682–687. DOI: https://doi.org/10.1038/nn1081, PMID: 12830159
Jun JJ, Steinmetz NA, Siegle JH, Denman DJ, Bauza M, Barbarits B, Lee AK, Anastassiou CA, Andrei A, Aydın Ç,
Barbic M, Blanche TJ, Bonin V, Couto J, Dutta B, Gratiy SL, Gutnisky DA, Häusser M, Karsh B, Ledochowitsch P,
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 31 of 42
Tools and resources Neuroscience
et al. 2017. Fully integrated silicon probes for high-density recording of neural activity. Nature 551:232–236.
DOI: https://doi.org/10.1038/nature24636, PMID: 29120427
Kao JC, Nuyujukian P, Ryu SI, Churchland MM, Cunningham JP, Shenoy KV. 2015. Single-trial dynamics of motor
cortex and their applications to brain-machine interfaces. Nature Communications 6:7759. DOI: https://doi.org/
10.1038/ncomms8759, PMID: 26220660
Kim TH, Zhang Y, Lecoq J, Jung JC, Li J, Zeng H, Niell CM, Schnitzer MJ. 2016. Long-term optical access to an
estimated one million neurons in the live mouse cortex. Cell Reports 17:3385–3394. DOI: https://doi.org/10.
1016/j.celrep.2016.12.004, PMID: 28009304
Kim J, Park H. 2008. Sparse nonnegative matrix factorization for clustering: Georgia Institute of Technology.
Technical Report GT-CSE.
Kim M, Smaragdis P. 2014. Efficient model selection for speech enhancement using a deflation method for
nonnegative matrix factorization. 2014 IEEE Global Conference on Signal and Information Processing
(GlobalSIP). DOI: https://doi.org/10.1109/GlobalSIP.2014.7032175
Lee DD, Seung HS. 1999. Learning the parts of objects by non-negative matrix factorization. Nature 401:788–
791. DOI: https://doi.org/10.1038/44565, PMID: 10548103
Lee DD, Seung HS. 2001. Algorithms for non-negative matrix factorization. In: Leen T. K, Dietterich T. G, Tresp V
(Eds). Advances in Neural Information Processing Systems. 13 MIT Press. p. 556–562.
Linderman S, Johnson M, Miller A, Adams R, Blei D, Paninski L. 2017. Bayesian learning and inference in
recurrent switching linear dynamical systems. Proceedings of the 20th International Conference on Artificial
Intelligence and Statistics.
Long MA, Jin DZ, Fee MS. 2010. Support for a synaptic chain model of neuronal sequence generation. Nature
468:394–399. DOI: https://doi.org/10.1038/nature09514, PMID: 20972420
Lynch GF, Okubo TS, Hanuschkin A, Hahnloser RH, Fee MS. 2016. Rhythmic continuous-time coding in the
songbird analog of vocal motor cortex. Neuron 90:877–892. DOI: https://doi.org/10.1016/j.neuron.2016.04.
021, PMID: 27196977
Maboudi K, Ackermann E, de Jong LW, Pfeiffer BE, Foster D, Diba K, Kemere C. 2018. Uncovering temporal
structure in hippocampal output patterns. eLife 7:e34467. DOI: https://doi.org/10.7554/eLife.34467, PMID: 2
9869611
MacDonald CJ, Lepage KQ, Eden UT, Eichenbaum H. 2011. Hippocampal "time cells" bridge the gap in memory
for discontiguous events. Neuron 71:737–749. DOI: https://doi.org/10.1016/j.neuron.2011.07.012, PMID: 21
867888
Mackevicius EL, Bahle AH, Williams AH. 2018. seqNMF. GitHub. 25df0d6. https://github.com/FeeLab/seqNMF
Mackevicius EL, Fee MS. 2018. Building a state space for song learning. Current Opinion in Neurobiology 49:59–
68. DOI: https://doi.org/10.1016/j.conb.2017.12.001, PMID: 29268193
Mizuseki, Sirota, Pastalkova, Diba, Buzsáki G. 2013. Multiple single unit recordings from different rat
hippocampal and entorhinal regions while the animals were performing multiple behavioral tasks. CRCNS.
https://crcns.org/data-sets/hc/hc-3c.
Młynarski W, McDermott JH. 2018. Learning midlevel auditory codes from natural sound statistics. Neural
Computation 30:631–669. DOI: https://doi.org/10.1162/neco_a_01048, PMID: 29220308
Mokeichev A, Okun M, Barak O, Katz Y, Ben-Shahar O, Lampl I. 2007. Stochastic emergence of repeating
cortical motifs in spontaneous membrane potential fluctuations in vivo. Neuron 53:413–425. DOI: https://doi.
org/10.1016/j.neuron.2007.01.017, PMID: 17270737
Netrapalli P, Niranjan UN, Sanghavi S, Anandkumar A, Jain P. 2014. Non-convex robust PCA. In: Ghahramani Z,
Welling M, Cortes C, Lawrence ND, Weinberger KQ (Eds). Advances in Neural Information Processing Systems.
27 Curran Associates, Inc. p. 1107–1115.
Okubo TS, Mackevicius EL, Payne HL, Lynch GF, Fee MS. 2015. Growth and splitting of neural sequences in
songbird vocal development. Nature 528:352–357. DOI: https://doi.org/10.1038/nature15741, PMID: 2661
8871
O’Grady PD, Pearlmutter BA. 2006. Convolutive non-negative matrix factorisation with a sparseness constraint.
16th IEEE Signal Processing Society Workshop on Machine Learning for Signal Processing. DOI: https://doi.
org/10.1109/MLSP.2006.275588
Pastalkova E, Itskov V, Amarasingham A, Buzsáki G. 2008. Internally generated cell assembly sequences in the
rat hippocampus. Science 321:1322–1327. DOI: https://doi.org/10.1126/science.1159775, PMID: 18772431
Pastalkova W, Wang Y, Mizuseki K, Buzsáki G. 2015. Simultaneous extracellular recordings from left and right
hippocampal areas CA1 and right entorhinal cortex from a rat performing a left / right alternation task and other
behaviors. CRCNS. https://crcns.org/data-sets/hc/hc-5.
Pearson K. 1901. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and
Dublin Philosophical Magazine and Journal of Science 2:559–572. DOI: https://doi.org/10.1080/
14786440109462720
Peter S, Kirschbaum E, Both M, Campbell L, Harvey B, Heins C, Durstewitz D. 2017. Sparse convolutional coding
for neuronal assembly detection. In: Guyon I, Luxburg U. V, Bengio S, Wallach H, Fergus R, Vishwanathan S,
Garnett R (Eds). Advances in Neural Information Processing Systems. 30 Curran Associates, Inc. p. 3675–3685.
Picardo MA, Merel J, Katlowitz KA, Vallentin D, Okobi DE, Benezra SE, Clary RC, Pnevmatikakis EA, Paninski L,
Long MA. 2016. Population-level representation of a temporal sequence underlying song production in the
zebra finch. Neuron 90:866–876. DOI: https://doi.org/10.1016/j.neuron.2016.02.016, PMID: 27196976
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 32 of 42
Tools and resources Neuroscience
Quaglio P, Rostami V, Torre E, Grün S. 2018. Methods for identification of spike patterns in massively parallel
spike trains. Biological Cybernetics 112:57–80. DOI: https://doi.org/10.1007/s00422-018-0755-0, PMID: 296515
82
Ramanarayanan V, Goldstein L, Narayanan SS. 2013. Spatio-temporal articulatory movement primitives during
speech production: extraction, interpretation, and validation. The Journal of the Acoustical Society of America
134:1378–1394. DOI: https://doi.org/10.1121/1.4812765, PMID: 23927134
Russo E, Durstewitz D. 2017. Cell assemblies at multiple time scales with arbitrary lag constellations. eLife 6:
e19428. DOI: https://doi.org/10.7554/eLife.19428, PMID: 28074777
Scholvin J, Kinney JP, Bernstein JG, Moore-Kochlacs C, Kopell N, Fonstad CG, Boyden ES. 2016. Close-packed
silicon microelectrodes for scalable spatially oversampled neural recording. IEEE Transactions on Biomedical
Engineering 63:120–130. DOI: https://doi.org/10.1109/TBME.2015.2406113, PMID: 26699649
Schrader S, Grün S, Diesmann M, Gerstein GL. 2008. Detecting synfire chain activity using massively parallel
spike train recording. Journal of Neurophysiology 100:2165–2176. DOI: https://doi.org/10.1152/jn.01245.2007,
PMID: 18632888
Sejnowski TJ, Churchland PS, Movshon JA. 2014. Putting big data to good use in neuroscience. Nature
Neuroscience 17:1440–1441. DOI: https://doi.org/10.1038/nn.3839, PMID: 25349909
Smaragdis P. 2004. Non-negative matrix factor deconvolution; extraction of multiple sound sources from
monophonic inputs. Berlin, Heidelberg: Springer. DOI: https://doi.org/10.1007/978-3-540-30110-3_63
Smaragdis P. 2007. Convolutive speech bases and their application to supervised speech separation. IEEE
Transactions on Audio, Speech and Language Processing 15:1–12. DOI: https://doi.org/10.1109/TASL.2006.
876726
Sutskever I, Vinyals O, Le QV. 2014. Sequence to sequence learning with neural networks. In: Ghahramani Z,
Welling M, Cortes C, Weinberger KQ, Lawrence ND (Eds). Advances in Neural Information Processing Systems.
27 MIT Press. p. 3104–3112.
Theis FJ, Stadlthanner K, Tanaka T. 2005. First results on uniqueness of sparse non-negative matrix factorization.
13th European Signal Processing Conference.
Torre E, Canova C, Denker M, Gerstein G, Helias M, Grün S. 2016. ASSET: analysis of sequences of synchronous
events in massively parallel spike trains. PLOS Computational Biology 12:e1004939. DOI: https://doi.org/10.
1371/journal.pcbi.1004939, PMID: 27420734
Ubaru S, Wu K, Bouchard KE. 2017. UoI-NMF cluster: a robust nonnegative matrix factorization algorithm for
improved parts-based decomposition and reconstruction of noisy data. 2017 16th IEEE International
Conference on Machine Learning and Applications (ICMLA).
Udell M, Horn C, Zadeh R, Boyd S. 2016. Generalized low rank models. Foundations and Trends in Machine
Learning 9:1–118. DOI: https://doi.org/10.1561/2200000055
van der Meij R, Voytek B. 2018. Uncovering neuronal networks defined by consistent between-neuron spike
timing from neuronal spike recordings. Eneuro 5:ENEURO.0379-17.2018. DOI: https://doi.org/10.1523/
ENEURO.0379-17.2018, PMID: 29789811
Vaz C, Toutios A, Narayanan S. 2016. Convex hull convolutive non-negative matrix factorization for uncovering
temporal patterns in multivariate time-series data. Interspeech 2016.
Villasana TPJ, Gorlow S, Hariraman AT. 2018. Multiplicative updates for convolutional NMF under b-Divergence.
arVix.
von Luxburg U. 2010. Clustering stability: an overview. Foundations and Trends in Machine Learning 2:235–274.
Wang D, Vipperla R, Evans N, Zheng TF. 2013. Online non-negative convolutive pattern learning for speech
signals. IEEE Transactions on Signal Processing 61:44–56. DOI: https://doi.org/10.1109/TSP.2012.2222381
Wold S. 1978. Cross-validatory estimation of the number of components in factor and principal components
models. Technometrics 20:397–405. DOI: https://doi.org/10.1080/00401706.1978.10489693
Wu S, Joseph A, Hammonds AS, Celniker SE, Yu B, Frise E. 2016. Stability-driven nonnegative matrix
factorization to interpret spatial gene expression and build local gene networks. PNAS 113:4290–4295.
DOI: https://doi.org/10.1073/pnas.1521171113, PMID: 27071099
Zhang Z, Xu Y, Yang J, Li X, Zhang D. 2016. A survey of sparse representation: algorithms and applications.
arXiv. https://arxiv.org/abs/1602.07017.
Zhou P, Resendez SL, Rodriguez-Romaguera J, Jimenez JC, Neufeld SQ, Giovannucci A, Friedrich J,
Pnevmatikakis EA, Stuber GD, Hen R, Kheirbek MA, Sabatini BL, Kass RE, Paninski L. 2018. Efficient and
accurate extraction of in vivo calcium signals from microendoscopic video data. eLife 7:e28728. DOI: https://
doi.org/10.7554/eLife.28728, PMID: 29469809
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 33 of 42
Tools and resources Neuroscience
Appendix 1
DOI: https://doi.org/10.7554/eLife.38471.032
Deriving multiplicative update rules
Standard gradient descent methods for minimizing a cost function must be adapted when
solutions are constrained to be non-negative, since gradient descent steps may result in
negative values. Lee and Seung invented an elegant and widely-used algorithm for non-
negative gradient descent that avoids negative values by performing multiplicative updates
(Lee and Seung, 2001; Lee and Seung, 1999). They derived these multiplicative updates by
choosing an adaptive learning rate that makes additive terms cancel from standard gradient
descent on the cost function. We will reproduce their derivation here, and detail how to
extend it to the convolutional case (Smaragdis, 2004) and apply several forms of
regularization (O’Grady and Pearlmutter, 2006; Ramanarayanan et al., 2013; Chen and
Cichocki, 2004). See Table 2 for a compilation of cost functions, derivatives and multiplicative
updates for NMF and convNMF under several different regularization conditions.
Standard NMF
NMF performs the factorization X »Xe ¼WH. NMF factorizations seek to solve the following
problem:
ðWf ;He Þ ¼ argminLðW;HÞ (17)
W;H
1
ð ; Þ ¼ jj 2L W H Xe XjjF (18)
2
Wf ;He  0 (19)
This problem is convex in W and H separately, not together, so a local minimum is found
by alternating W and H updates. Note that:
q
ð ; Þ ¼ e > >L W H XH XH (20)
qW
q
ð ; Þ ¼ > >L W H W Xe W X (21)
qH
Thus, gradient descent steps for W and H are:
W Wh > >WðXeH XH Þ (22)
 > >H HhHðW Xe W XÞ (23)
To arrive at multiplicative updates, Lee and Seung (Lee and Seung, 2001) set:
W
hW ¼ > (24)
WHH
H
hH ¼ > (25)
W WH
Thus, the gradient descent updates become multiplicative:
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 34 of 42
Tools and resources Neuroscience
> >
XH XH
W  W > ¼W e (26)WHH >XH
> >
W X W X
H  H > ¼ H
W WH >W Xe
(27)
where the division and  are element-wise.
Standard convNMF
P ‘!
Convolutional NMF factorizes data X Xe» ¼ ‘ W

‘H ¼W
 H. ConvNMF factorizations seek
to solve the following problem:
ðWf ;He Þ ¼ argminLðW;HÞ (28)
W;H
1
LðW;HÞ ¼ jjXe Xjj2F (29)
2
Wf ;He  0 (30)
The derivation above for standard NMF can be applied for each ‘, yielding the following
update rules for convNMF (Smaragdis, 2004):
‘!>
XH
W

‘ W

‘ (31)
e ‘!
>
XH
P  ‘ >>
W X W
‘ 

‘  X
H H ¼ H (32)
P  ‘ >>
W e W
 Xe‘ 

‘X
Where the operator ‘! shifts a matrix in the! direction by ‘ timebins, that is a delay by ‘
timebins, and ‘ shifts a matrix in the direction by ‘ timebins (Table 1). Note that NMF is
a special case of convNMF where L ¼ 1.
Incorporating regularization terms
Suppose we want to regularize by adding a new term, R to the cost function:
ðWf ;He Þ ¼ argminLðW;HÞ (33)
W;H
1
LðW;HÞ ¼ jjXe  2XjjF þR (34)
2
Wf ;He  0 (35)
Using a similar trick to Lee and Seung, we choose a hW;hH to arrive at a simple
multiplicative update. Below is the standard NMF case, which generalizes trivially to the
convNMF case.
Note that:
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 35 of 42
Tools and resources Neuroscience
qL e > qR¼ XH  >XH þ (36)
qW qW
qL
¼ > e  > qRW X W Xþ (37)
qH qH
We set:
W
hW ¼ e (38)>þ qRXH
qW
H
hH ¼ >
W Xe
(39)
þ qR
qH
Thus, the gradient descent updates become multiplicative:
>
qL XH
W WhW ¼W e (40)qW >XH þ qR
qW
>
qL W X
H HhH ¼ H (41)
qH > e þ qRW X
qH
where the division and  are element-wise.
The above formulation enables flexible incorporation of different types of regularization or
penalty terms into the multiplicative NMF update algorithm. This framework also extends
naturally to the convolutional case. See Table 2 for examples of several regularization terms,
including L1 sparsity (O’Grady and Pearlmutter, 2006; Ramanarayanan et al., 2013) and
spatial decorrelation (Chen and Cichocki, 2004), as well as the terms we introduce here to
combat the types of inefficiencies and cross correlations we identified in convolutional NMF,
namely, smoothed orthogonality for H and W, and the x-ortho penalty term. For the x-ortho
>
penalty term, ljjð 
 Þ >W X SH jj ;i 6¼j, the multiplicative update rules are:1
 
‘! >
X H
W

‘ W

‘   (42)
‘! >
Xe
 ‘
H þl >XSH ð1 IÞ
>
W
 X
H H (43)
> >
W
 Xe þlð1 IÞðW
 XSÞ
where the division and  are element-wise. Note that multiplication with the K  K matrix
ð1 IÞ effectively implements factor competition because it places in the kth row a sum across
all other factors.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 36 of 42
Tools and resources Neuroscience
Appendix 2
DOI: https://doi.org/10.7554/eLife.38471.032
Relation of the x-ortho penalty to traditional
regularizations
As noted in the main text, the x-ortho penalty term is not formally a regularization because it
includes the data X. In this Appendix, we show how this penalty can be approximated by a
data-free regularization. The resulting regularization contains three terms corresponding to a
weighted orthogonality penalty on pairs of H factors, a weighted orthogonality penalty on
pairs of W factors, and a term that penalizes interactions among triplets of factors. We
analyze each term in both the time domain (Equation 50) and in the frequency domain
(Equations 50 and 69).
Time domain analysis
We consider the cross-orthogonality penalty term:
>
R¼ kðW
 Þ >X SH k ;i 6¼j (44)1
>
and define, R ¼ ðW
 >XÞSH , which is a K  K matrix. Each element Rij is a positive number
describing the overlap or correlation between factor i and factor j in the model. Each element
of R can be written explicitly as:
XXX X
Rij ¼ Wni‘Xnðtþ‘Þ St Hj (45)t t
t n ‘ t
Where the index variables t and t range from 1 to T, n ranges from 1 to N, and ‘ ranges
from 1 to L.
Our goal here is to find a close approximation to this penalty term that does not contain
the data X. This can readily be done if X is well-approximated by the convNMF
decomposition:
XX
Xnt 
» ðW  HÞnt ¼ Wnk‘Hkðt‘Þ (46)
k ‘
P
Substituting this expression into Equation 45 and defining the smoothed matrix St Ht jt
t
as smoothHjt gives:
XXXXX
smooth
Rij » Wni‘Wnk‘0HktHjt (47)
t n ‘ k ‘0
Making the substitution u ¼ ‘ ‘0 gives:
XX XL1 XX
smooth
Rij » Wnið‘0þuÞWnk‘0HkðtþuÞHjt (48)
t n u¼ðL1Þ k ‘0
where in the above expression we have taken u ¼ ‘ ‘0 to extend over the full range from
ðL 1Þ to ðL 1Þ under the implicit assumption that W and H are zero padded such that
values of W for lag indices outside the range 0 to L 1 and values of H for time indices
outside the range 1 to T are taken to be zero.
Relabeling ‘0 as ‘ and gathering terms together yields
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 37 of 42
Tools and resources Neuroscience
 ! !
X XL1 XX X
smooth
Rij » Wnið‘þuÞWnk‘ HkðtþuÞHjt (49)
k u¼ðL1Þ n ‘ t
We note that the above expression contains terms that resemble penalties on orthogonality
between two W factors (first parenthetical) or two H factors (one of which is smoothed,
second parenthetical), but in this case allowing for different time lags u between the factors.
To understand this formula better, we decompose the above sum over k into three
contributions corresponding to k ¼ i, k ¼ j, and k ¼6 i; j
 ! !
XL1 XX X
smooth
Rij » Wnið‘þuÞWni‘ HiðtþuÞHjt þ
u¼ðL1Þ n t ‘ ! !
XL1 XX X
smooth
Wnið‘þuÞWnj‘ HjðtþuÞHjt þ (50)
u¼ðL1Þ n t ‘ ! !
X XL1 XX X
smooth
Wnið‘þuÞWnk‘ HkðtþuÞHjt
k¼6 i;j u¼ðL1Þ n ‘ t
The first term above contains, for u ¼ 0, a simple extension of the ð ; Þthi j element of the H
orthogonality condition >HSH . The extension is that the orthogonality is weighted by the
power, that is the sum of squared elements, in the ith factor of W (the apparent lack of
symmetry in weighing by the ith rather than the jth factor can be removed by simultaneously
considering the term Rji, as shown in the Fourier representation of the following section). This
weighting has the benefit of applying the penalty on H orthogonality most strongly to those
factors whose corresponding W components contain the most power. For u 6¼ 0, this
orthogonality condition is extended to allow for overlap of time-shifted H’s, with weighting at
each time shift by the autocorrelation of the corresponding W factor. Qualitatively, this
enforces that (even in the absence of the smoothing matrix S), H’s that are offset by less than
the width of the autocorrelation of the corresponding W’s will have overlapping convolutions
with these W’s due to the temporal smoothing associated with the convolution operation. We
note that, for sparse sequences as in the examples of Figure 1, there is no time-lagged
component to the autocorrelation, so this term corresponds simply to a smoothed H
orthogonality regularization, weighted by the strength of the corresponding W factors.
The second term above represents a complementary orthogonality condition on the W
components, in which orthogonality in the ith and jth factors are weighted by the (smoothed)
autocorrelation of the H factors. For the case in which the H factors have no time-lagged
autocorrelations, this corresponds to a simple weighting of W orthogonality by the strength
of the corresponding H factors.
Finally, we consider the remaining terms of the cost function, for which k 6¼ i; j. We note
that these terms are only relevant when the factorization contains at least three factors, and
thus their role cannot be visualized from the simple Type 1 to Type 3 examples of Figure 1.
These terms have the form:
 ! !
X XL1 XX X
smooth
Wnið‘þuÞWnk‘ HkðtþuÞHjt (51)
k 6¼i;j u¼ðL1Þ n ‘ t
To understand how this term contributes, we consider each of the expressions in
parentheses. The first expression corresponds, as described above, to the time-lagged cross
correlation of the ith and kth W components. Likewise, the second expression corresponds to
the time-lagged correlation of the (smoothed) jth and kth H components. Thus, this term of Rij
contributes whenever there is a factor (W
k
;Hk
) that overlaps, at the same time lags, with the
ith factor’s W component and the jth factor’s H component. Thus, this term penalizes cases
where, rather than (or in addition to) two factors i and j directly overlapping one another, they
have a common factor k with which they overlap.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 38 of 42
Tools and resources Neuroscience
An example of the contribution of a triplet penalty term, as well as of the paired terms of
Equation 50, is shown in Figure 1 of this Appendix. By inspection, there is a penalty R23 due
to the overlapping values of the pair (h2;h3). Likewise, there is a penalty R13 due to the
overlapping values of the pair (w1;w3). The triplet penalty term contributes to R12 and derives
from the fact that w1 overlaps with w3 at the same time (and with the same, zero time lag) as
h2 overlaps with h3.
In summary, the above analysis shows that for good reconstructions of the data where
X W
»  H, the x-ortho penalty can be well-approximated by the sum of three contributions.
The first corresponds to a penalty on time-lagged (smoothed) H orthogonality weighted at
each time lag by the autocorrelation of the corresponding W factors. The second similarly
corresponds to a penalty on time-lagged W orthogonality weighted at each time lag by the
(smoothed) autocorrelation of the corresponding H factors. For simple cases of sparse
sequences, these contributions reduce to orthogonality in H or W weighted by the power in
the corresponding W or H, respectively, thus focusing the penalties most heavily on those
factors which contribute most heavily to the data reconstruction. The third, triplet contribution
corresponds to the case in which a factor in W and a different factor in H both overlap (at the
same time lag) with a third common factor, and may occur even when the factors W and H
themselves are orthogonal. Further work is needed to determine whether this third
contribution is critical to the x-ortho penalty or is simply a by-product of the x-ortho penalty
procedure’s direct use of the data X.
h1
h2
h3
w1 w2 w3
Appendix 2—figure 1. Example of redundancy with three factors. In addition to the direct
overlap of h2 and h3, and of w1 and w3, there is a ‘triplet’ penalty R12 on factors 1 and 2 that
occurs because each has an overlap (in either W or H) with the 3rd factor (w3;h3). This occurs
even though neither w1 and w2, nor h1 and h2, are themselves overlapping.
DOI: https://doi.org/10.7554/eLife.38471.034
Frequency domain analysis
Additional insight may be obtained by analyzing these three components of R in the Fourier
domain. Before doing so, we below derive the Fourier domain representation of R, and
provide insights suggested by this perspective.
Fourier representation of the x-ortho penalty
As in the time domain analysis, we start with defining:
XXX X
Rij ¼ Wni‘Xnðtþ‘Þ St Hj (52)t t
t n ‘ t
Unpacking the notation above, we note that:
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 39 of 42
Tools and resources Neuroscience
" #
XN ƒ!ðnÞ !ðnÞ !
Rij ¼ CorrðW i ;X Þ 
ConvðHj;
!
s Þ (53)
n¼1
ƒ!ðnÞ !ðnÞ !
where W i is the n
th row of , is the nth !W
i
 X row of X, Hj is Hj
, s is a smoothing vector
corresponding to the entries of each row of the smoothing matrix S, and “
” is a dot product.
For ease of mathematical presentation, in the following, we work with continuous time rather
than the discretely sampled data and extend the W factors, H factors, and data matrix X
through zero-padding on both ends so that:
ƒ Z!ðnÞ !ðnÞ ¥ ƒ!ðnÞ!ðnÞ
CorrðW i ;X ÞðtÞ ¼ W i‘ X‘þtd‘ (54)
¥
and
Z
! ¥!
ConvðHj;
!
s Þ ¼ !Hj st t dt (55)t
¥
Recall that the Fourier transform is defined as:
Z
¥
f̂ ð!Þ  Fðf ðtÞÞ  f ðtÞei!tdt (56)
¥
with inverse Fourier transform:
Z
¥
1
f ðtÞ ¼F1ðf̂ ð!ÞÞ  f̂ ð!Þeþi!td! (57)
2p ¥
Now recall some basic features of Fourier transforms of correlation and convolution
integrals:
FðConvðf ðtÞ;gðtÞÞÞ ¼ f̂ ð!Þĝð!Þ (58)
FðCorrðf ðtÞ;gðtÞÞÞ ¼ f̂ ð!Þĝð!Þ (59)
Z
¥  
f ðtÞ 
 gðtÞ ¼ f ðtÞgðtÞdt¼Corrt¼0ðf ;gÞ ¼ F
1 f̂ ð!Þĝð!Þd!
t¼0
¥ Z (60)
¥
1
f ðtÞ 
 gðtÞ ¼ f̂ ð!Þĝð!Þd!
2p ¥
This final identity, known as Parseval’s theorem, says that the inner product (dot product)
between two functions evaluated in the time and frequency domain are equivalent up to a
proportionality constant of 1=ð2pÞ. With the above identities, we can calculate our quantity of
interest:
" #
XN ƒ!ðnÞ !ðnÞ !
Rij ¼ CorrðW i ;X Þ 
Convð
!
Hj; s Þ (61)
n¼1
First, define:
" #
XN ƒ!ðnÞ !ðnÞ
AðtÞ ¼ CorrðW i ;X Þ (62)
n¼1
!
BðtÞ ¼Convð !Hj; s Þ (63)
From Equation 60 (Parseval’s theorem):
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 40 of 42
Tools and resources Neuroscience
Z
¥
1
Rij ¼ Âð!ÞB̂ð!Þd! (64)
2p ¥
Then, from Equations 58 and 59, we have:
1 XN Z ¥ ðnÞ
Rij ¼ d!Ŵi ð!Þ
ðnÞ
X̂ ð!ÞĤjð!Þ̂sð!Þ (65)
2p
n¼1 ¥
The above formula shows that:
1. Viewed in the frequency domain, the x-ortho penalty reduces to a (sum over neurons and fre-
quencies of a) simple product of Fourier transforms of the four matrices involved in the
penalty.
2. The smoothing can equally well be applied to H or W or X. (For X, note that for symmetric
smoothing function sðtÞ ¼ sðtÞ, we also have ŝð!Þ ¼ ŝð!Þ.)
3. One can view this operation as either of the below:
a. First correlate W and X by summing correlations of each row, and then calculate the
overlap with the h smoothed iH, as described in the mainP R
1 N ¥ ðnÞ
 
text: ¼ d! ð!Þ ðnÞRij n¼  Ŵi X̂ ð!Þ Ĥjð!Þ̂sð!Þ2p 1 ¥
b. Correlate H with each row of X and then calculate the overlap of this correlation with the
corresponding smoothed row h of W. i Then sum over allP R  
¥ ðnÞ
rows:Rij ¼
1 N ðnÞ
n¼  d! Ĥjð!ÞX̂ ð!Þ Ŵi ð!Þ̂sð!Þ2p 1 ¥
Fourier representation of the traditional regularization approximation of
the x-ortho penalty
We now proceed to show how the x-ortho penalty can be approximated by a traditional (data-
free) regularization, expressing the results in the frequency domain. As in the time domain
analysis, we consider the approximation in which the data X are nearly perfectly reconstructed
by the convNMF decomposition (X »W
 H).
Noting that this decomposition is a sum over factors of row-by-row ordinary convolutions,
we can write the Fourier analog for each row of X as:
X
ðnÞ ðnÞ
X̂ ð!Þ» Ŵk ð!ÞĤkð!Þ (66)
k
Thus, substituting X with the reconstruction, W
 H in Equation 65, we have:
Z " #X ¥ N1 X ðnÞ ðnÞ  
Rij » d! Ŵi ð!ÞŴk ð!Þ Ĥkð!ÞĤjð!Þŝð!Þ (67)
2p
k ¥ n¼1
As in the time domain analysis, we separate the sum over k into three cases: k ¼ i, k ¼ j,
and k ¼6 i; j. Recall that for real numbers, f̂ ð!Þ ¼ f̂ ð!Þ, and f̂ ð!Þf̂ ð!Þ ¼ jf̂ ð!Þj2. Thus,
separating the sum over k into the three cases, we have:
Z " #
¥ XN  1  ðnÞ 2  
Rij ¼ d! Ŵi ð!Þ Ĥið!ÞĤjð!Þŝð!Þ þ
2p ¥ n¼1
Z " # (68)
¥ N
1   X2 ðnÞ ðnÞĤjð!Þ Ŵi ð!ÞŴj ð!Þŝð!Þ þY
2p ¥ n¼1
where Y represents the remaining terms for which k 6¼ i; j.
We can obtain a more symmetric form of this equation by summing the contributions of
factors i and j, Rij þRji. For symmetric smoothing functions sðtÞ ¼ sðtÞ, for which
ŝð!Þ ¼ ŝð!Þ, we obtain:
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 41 of 42
Tools and resources Neuroscience
Z "
¥ XN 
#
1    2ðnÞ  2ðnÞ   
RijþRji ¼ d! Ŵ ð!Þ þŴ ð!Þ Ĥið!ÞĤjð!Þŝð!Þ þ
2p i j¥ n¼1
Z h " # (69)¥
1  2  
i XN
2 ðnÞ ðnÞ
Ĥið!Þ þ Ĥjð!Þ Ŵi ð!ÞŴj ð!Þŝð!Þ þY
2p ¥ n¼1
As in the time domain analysis, the first two terms above have a simple interpretation in
comparison to traditional orthogonality regularizations: The first term resembles a traditional
regularization of orthogonality in (smoothed) H, but now weighted frequency-by-frequency by
the summed power at that frequency in the corresponding W factors. For sparse (delta-
function-like) sequences, the power in W at each frequency is a constant and can be taken
outside the integral. In this case, the regularization corresponds precisely to orthogonality in
(smoothed) H, weighted by the summed power in the corresponding W’s. Likewise, the
second term above corresponds to a traditional regularization of orthogonality in (smoothed)
W, weighted by the summed power at each component frequency in the corresponding H
factors.
Altogether, we see that these terms represent a Fourier-power weighted extension of
(smoothed) traditional orthogonality regularizations in W and H. This weighting may be
beneficial relative to traditional orthogonality penalties, since it makes the regularization focus
most heavily on the factors and frequencies that contribute most to the data reconstruction.
Finally, we consider the remaining terms in the cost function, for which k 6¼ i; j. As noted
previously, these terms are only relevant when the factorization contains at least three terms,
so cannot be seen in the simple Type 1, 2 and 3 cases illustrated in Figure 1. These terms
have the form:
XZ
" #
¥
1 XN ðnÞ ðnÞ  
d! Ŵi ð!ÞŴk ð!Þ Ĥkð!ÞĤjð!Þŝð!Þ (70)
2p
k 6¼i;j ¥ n¼1
To understand how this term contributes, we consider each of the expressions in
parentheses. The first expression contains each frequency component of the correlation of the
ith and kth factors’ W’s. The second expression likewise contains each frequency component of
the correlation of the jth and kth factors’ H’s. Thus, analogous to the time domain analysis, this
term of Rij contributes whenever there is a factor (W
k
;Hk
) that overlaps at any frequency
with the ith factor’s component and the jthW factor’s H component. In this manner, this three-
factor interaction term effectively enforces competition between factors i and j even if they are
not correlated themselves, as demonstrated in Figure 1 of this Appendix.
Mackevicius et al. eLife 2019;8:e38471. DOI: https://doi.org/10.7554/eLife.38471 42 of 42