Matrix factorization underlies many well-known unsupervised learning algorithms, including PCA (Pearson, 1901), non-negative matrix factorization (NMF) (Lee and Seung, 1999), dictionary learning, and k-means clustering (see Udell et al., 2016 for a review). We start with a data matrix, $\mathbf{𝐗}$, containing 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 $\mathbf{𝐰}$ representing the neural pattern, and a row vector $\mathbf{𝐡}$ representing the times and amplitudes at which that pattern occurs (temporal loadings). In this case, the data matrix $\mathbf{𝐗}$ is mathematically reconstructed as the outer product of $\mathbf{𝐰}$ and $\mathbf{𝐡}$. 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:

where ${\mathbf{𝐗}}_{nt}$ is the ${(nt)}^{th}$ element of matrix $\mathbf{𝐗}$, that is, the activity of neuron $n$ at time $t$. Here, in order to store $K$ different patterns, $\mathbf{𝐖}$ is a $N\times K$ matrix containing the $K$ exemplar patterns, and $\mathbf{𝐇}$ is a $K\times T$ matrix containing the $K$ timecourses:

Given a data matrix with unknown patterns, the goal of matrix factorization is to discover a small set of patterns, $\mathbf{𝐖}$, and a corresponding set of temporal loading vectors, $\mathbf{𝐇}$, that approximate the data. In the case that the number of patterns, $K$, is sufficiently small (less than $N$ and $T$), this corresponds to a dimensionality reduction, whereby the data is expressed in more compact form. PCA additionally requires that the columns of $\mathbf{𝐖}$ and the rows of $\mathbf{𝐇}$ are orthogonal. NMF instead requires that the elements of $\mathbf{𝐖}$ and $\mathbf{𝐇}$ are nonnegative. The discovery of unknown factors is often accomplished by minimizing the following cost function, which measures the element-by-element sum of all squared errors between a reconstruction $\widetilde{\mathbf{𝐗}}=\mathbf{𝐖𝐇}$ and the original data matrix $\mathbf{𝐗}$ using the Frobenius norm, ${\parallel \mathbf{𝐌}\parallel }_F=\sqrt{{\sum }_{ij}{\mathbf{𝐌}}_{ij}^2}$:

(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 ${\mathbf{𝐖}}^{*}$ and ${\mathbf{𝐇}}^{*}$ that minimize this cost function produce an optimal reconstruction ${\widetilde{\mathbf{𝐗}}}^{*}={\mathbf{𝐖}}^{*}{\mathbf{𝐇}}^{*}$. Iterative optimization methods such as gradient descent can be used to search for global minima of the cost function; however, it is often possible 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 discovering temporally extended patterns—first, because each element in a sequence must be represented 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 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 $\mathbf{𝐖}$ is represented by a single vector (Figure 1A), the factors $\mathbf{𝐖}$ in convNMF represent patterns of neural activity over a brief period of time. Each pattern is stored as an $N\times L$ matrix, ${\mathbf{𝐰}}_{\mathbf{𝐤}}$, where each column (indexed by $\mathrm{\ell }=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 ${\mathbf{𝐡}}_{\mathrm{𝟏}}$, as for NMF. The reconstruction is produced by convolving the $N\times L$ pattern with the time series ${\mathbf{𝐡}}_{\mathrm{𝟏}}$ (Figure 1B).

If the data contains multiple patterns, each pattern is captured by a different $N\times L$ matrix and a different associated time series vector $\mathbf{𝐡}$. A collection of $K$ different patterns can be compiled together into an $N\times K\times L$ array (also known as a tensor), $\mathbf{𝐖}$ and a corresponding $K\times T$ time series matrix $\mathbf{𝐇}$. Analogous to NMF, convNMF generates a reconstruction of the data as a sum of $K$ convolutions between each neural activity pattern ($\mathbf{𝐖}$), and its corresponding temporal loadings ($\mathbf{𝐇}$):

The tensor/matrix convolution operator ${\displaystyle ⊛}$ (notation summary, Table 1) reduces to matrix multiplication 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 $\mathbf{𝐖}$ and $\mathbf{𝐇}$ 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 challenging 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 patterns on each run (Peter et al., 2017). These inconsistency errors fall into three main categories (Figure 1C):

• 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 ($\mathbf{𝐇}$) and/or their exemplar activity patterns ($\mathbf{𝐖}$).

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 strategy 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 simplifying 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 (stability NMF; Wu et al., 2016).

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, $ℛ$, into the convNMF cost function:

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 $\mathbf{𝐇}$’s (Choi, 2008; Chen and Cichocki, 2004). Such correlations are measured by computing the correlation matrix ${\mathbf{𝐇𝐇}}^{\top }$, which contains the correlations between the temporal loadings of every pair of factors. The regularization may be implemented using the penalty term $ℛ=\lambda {\parallel {\mathbf{𝐇𝐇}}^{\top }\parallel }_{1,i\ne j}$, where the seminorm $\parallel \cdot {\parallel }_{1,i\ne j}$ sums the absolute value of every matrix entry except those along the diagonal (notation summary, 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 parameter $\lambda$ controls the magnitude of the penalty term $ℛ$.

In convNMF, a penalty term based on ${\mathbf{𝐇𝐇}}^{\top }$ 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 ${\mathbf{𝐇𝐇}}^{\top }$. One simple way to address this problem is to smooth the $\mathbf{𝐇}$’s in the penalty term with a square window of length $2L-1$ using the smoothing matrix $\mathbf{𝐒}$ (${\mathbf{𝐒}}_{ij}=1$ when $|i-j|<L$ and otherwise ${\mathbf{𝐒}}_{ij}=0$). The resulting penalty, $ℛ=\lambda \parallel {\mathbf{𝐇𝐒𝐇}}^{\top }\parallel$, 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 patterns in $\mathbf{𝐖}$ 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 $\parallel {\mathbf{𝐇𝐒𝐇}}^{\top }\parallel$ 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 $\mathbf{𝐖}$. This has the disadvantage of requiring an extra parameter, namely the $\lambda$ associated with this cost.

Instead we chose an alternative approach to resolve errors of Type 3 that simultaneously detects correlations in $\mathbf{𝐖}$ and $\mathbf{𝐇}$ using a single cross-orthogonality cost term. We note that, for Type 3 errors, redundant $\mathbf{𝐖}$ 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 $\mathbf{𝐖}$, its overlap with the data at time $t$. This quantity is captured in symbolic form by $\mathbf{𝐖}\stackrel{\top }{⊛}\mathbf{𝐗}$ (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. This cross-orthogonality penalty term, which we refer to as 'x-ortho’, sums up these correlations across all pairs of factors, implemented as follows:

When incorporated into the update rules, this causes any factor that has a high overlap with the data to suppress the temporal loadings ($\mathbf{𝐇}$) 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 data, favoring solutions that use a minimal set of factors to give a good reconstruction. The resulting global cost function is:

The update rules for $\mathbf{𝐖}$ and $\mathbf{𝐇}$ 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 $\mathbf{𝐗}$, rather than just the model variables $\mathbf{𝐖}$ and $\mathbf{𝐇}$. 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 $\mathbf{𝐖}$ and a term corresponding to a weighted smoothed) orthogonality penalty on $\mathbf{𝐇}$, consistent with the observation that the x-ortho penalty simultaneously punishes factor correlations in $\mathbf{𝐖}$ and $\mathbf{𝐇}$.

There is an interesting relation between our method for penalizing correlations and other methods 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 $\mathbf{𝐖}$ and $\mathbf{𝐇}$, necessitating the selection of two hyper-parameters (${\lambda }_w$ and ${\lambda }_h$) (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 $\mathbf{𝐖}$ and $\mathbf{𝐇}$ simultaneously. The x-ortho penalty described above serves this purpose (Equation 6). As we will describe below, the application of sparsity penalties can be very useful for shaping the factors produced by convNMF, and our code includes options for applying sparsity penalties on both $\mathbf{𝐖}$ and $\mathbf{𝐇}$.

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 discover 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 containing 10 neurons in the absence of noise (Figure 2A). The resulting data matrix had a total duration 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.

Figure 2 with 2 supplements see all

Download asset Open asset

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 (${\mathbf{𝐰}}_1$, ${\mathbf{𝐰}}_2$, and ${\mathbf{𝐰}}_3$ 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 factorizations from different random initializations when $K$ was chosen to be too large, due to the ambiguities 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 datasets 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 determined 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).

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, factorizations using the x-ortho penalty achieved near-perfect reconstruction error, and consistency was 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 neurons varied across the sequence (Figure 3—figure supplement 3).

Figure 3 with 3 supplements see all

Download asset Open asset

Figure 4 with 5 supplements see all

Download asset Open asset

An alternative strategy to minimizing redundant factorizations is to estimate the number of underlying 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 used in clustering models (von Luxburg, 2010). The stability NMF score, diss, is measured between two factorizations, ${\mathbf{𝐅}}^1=\{{\mathbf{𝐖}}^1,{\mathbf{𝐇}}^1\}$ and ${\mathbf{𝐅}}^2=\{{\mathbf{𝐖}}^2,{\mathbf{𝐇}}^2\}$, run from different initial conditions:

where $\mathbf{𝐂}$ is the cross-correlation matrix between the columns of the matrix ${\mathbf{𝐖}}^1$ and the the columns of the matrix ${\mathbf{𝐖}}^2$. Note that diss is low when there is a one-to-one mapping between factors in ${\mathbf{𝐅}}^1$ and ${\mathbf{𝐅}}^2$, 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 ${\mathbf{𝐡}}_k$ by one time step while shifting the elements of ${\mathbf{𝐰}}_k$ by one step in the opposite direction with little change to the model reconstruction. Thus, rather than computing correlations from the factor patterns $\mathbf{𝐖}$ or loadings $\mathbf{𝐇}$, we computed the diss metric using correlations between factor reconstructions (${\widetilde{\mathbf{𝐗}}}_k={\mathbf{𝐰}}_{\mathbf{𝐤}}⊛{\mathbf{𝐡}}_{\mathbf{𝐤}}$).

where $\text{Tr}[\cdot ]$ denotes the trace operator, $\text{Tr}[\mathbf{𝐌}]={\sum }_i{\mathbf{𝐌}}_{ii}$. That is, ${\mathbf{𝐂}}_{ij}$ measures the correlation between the reconstruction of factor i in ${\mathbf{𝐅}}^1$ and the reconstruction of factor j in ${\mathbf{𝐅}}^2$. Here, as for stability 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 reconstruction 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).

Figure 5 with 2 supplements see all

Download asset Open asset

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 reasonable 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 mechanisms. ‘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.

Figure 6 with 1 supplement see all

Download asset Open asset

Here, we show that the addition of penalties on either $\mathbf{𝐖}$ or $\mathbf{𝐇}$ 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-based’, or ‘events-based’ depending on initial conditions and the structure of shared neurons activities. This approach works because, in ‘events-based’ factorization, the $\mathbf{𝐇}$’s are orthogonal (uncorrelated) while the $\mathbf{𝐖}$’s have high overlap; conversely, in the ‘parts-based’ factorization, the $\mathbf{𝐖}$’s are orthogonal while the $\mathbf{𝐇}$’s are strongly correlated. Note that these correlations in $\mathbf{𝐖}$ or $\mathbf{𝐇}$ 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 $\mathbf{𝐖}$ or $\mathbf{𝐇}$ are provided in Table 2 with derivations in Appendix 1. Figure 9—figure supplement 2 shows examples of using these penalties 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$ regularization 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 $\mathbf{𝐖}$ and $\mathbf{𝐇}$ 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 $\mathbf{𝐖}$ update will bias the $\mathbf{𝐖}$ exemplars to include only a single repetition of the repeated sequence, while the addition of a sparsity regularizer on $\mathbf{𝐇}$ will bias the $\mathbf{𝐖}$ 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.

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, reconstruction 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. patterns 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 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 column-shuffled versus non-shuffled data matrices (see Materials and methods), and quantify the performance 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 degree to which sequences are present in a dataset. Below, we apply this method to real experimental data and obtain high ‘sequenciness’ scores, suggesting that convolutional matrix factorization is a well-suited tool for summarizing neural dynamics in these datasets.

Figure 7

Download asset Open asset

To test the ability of seqNMF to discover patterns in electrophysiological data, we analyzed multielectrode 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 alternations, the rats ran on a running wheel during an imposed delay period lasting either 10 or 20 seconds. 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 $\lambda =2{\lambda }_0$, 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 dataset 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 $\lambda$, factorizations more often split these sequences into two factors. At even lower values of $\lambda$, factorizations had even more significant 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 supplement 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 $\lambda =1.5{\lambda }_0$, 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.

Figure 8

Download asset Open asset

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 pathological 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 $\lambda$ (Figure 9B), x-ortho penalized factorizations typically extracted three significant sequences (Figure 9C). Similarly, the diss measure has a local minimum at $K=3$ (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 sequences deployed during syllables of abnormally long and variable durations (Figure 9D–F, Figure 9—figure supplement 1A).

Figure 9 with 2 supplements see all

Download asset Open asset

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, 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.

Although we have focused on the application of seqNMF to neural activity data, these methods naturally 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 recordings 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 previously 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 natural statistics (Młynarski and McDermott, 2018).

Figure 10

Download asset Open asset

Further requests should be directed to Michale Fee (fee@mit.edu).

The seqNMF MATLAB code is publicly available as a github repository, which also includes our songbird 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 $\lambda$, testing the significance of factors, plotting, and other functions. It also includes a demo script with an example of how to select $\lambda$ for a new dataset, test for significance of factors, plot the seqNMF factorization, switch between parts-based and events-based factorizations, and calculate cross-validated performance on a masked test set.

We simulated neural sequences containing between 1 and 10 distinct neural sequences in the presence of various noise conditions. Each neural sequence was made up of 10 consecutively active neurons, each separated by three timebins. The binary activity matrix was convolved with an exponential kernel ($\tau =10$ timebins) to resemble neural calcium imaging activity.

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). However, 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 introducing 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 consistently 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:

We derived the following multiplicative update rules for $\mathbf{𝐖}$ and $\mathbf{𝐇}$ (Appendix 1):

where the division and $\times$ are element-wise. The operator $\stackrel{\mathrm{\ell }\to }{(\cdot )}$ shifts a matrix in the $\to$ direction by $\mathrm{\ell }$ timebins, that is a delay by $\mathrm{\ell }$ timebins, and $\stackrel{\leftarrow \mathrm{\ell }}{(\cdot )}$ shifts a matrix in the $\leftarrow$ direction by $\mathrm{\ell }$ timebins (notation summary, Table 1). Note that multiplication with the $K\times K$ matrix $(\mathrm{𝟏}-\mathbf{𝐈})$ effectively implements factor competition because it places in the $k$th 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 $\mathbf{𝐇}$ have unit norm; shift factors to be centered in time such that the center of mass of each $\mathbf{𝐖}$ 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 $\mathbf{𝐇}$ that may occur within $2L$ timebins of neighboring sequences.

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 dataset 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.

We used male zebra finches (Taeniopygia guttata) from the MIT zebra finch breeding facility (Cambridge, 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 custom-made sound isolation chambers.

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 window was made using a GRIN (gradient index) lens (1 mm diamenter, 4 mm length, Inscopix). After at 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 algorithm, 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 ${95}^{th}$ 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.

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.

As noted in the main text, the x-ortho penalty term is not formally a regularization because it includes the data $\mathbf{𝐗}$. 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 $\mathbf{𝐇}$ factors, a weighted orthogonality penalty on pairs of $\mathbf{𝐖}$ 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:

(44) $ℛ={\parallel (\mathbf{𝐖}\stackrel{\top }{⊛}\mathbf{𝐗}){\mathbf{𝐒𝐇}}^{\top }\parallel }_{1,i\ne j}$

and define, $\mathbf{𝐑}=(\mathbf{𝐖}\stackrel{\top }{⊛}\mathbf{𝐗}){\mathbf{𝐒𝐇}}^{\top }$, which is a $K\times K$ matrix. Each element ${\mathbf{𝐑}}_{ij}$ is a positive number describing the overlap or correlation between factor $i$ and factor $j$ in the model. Each element of $\mathbf{𝐑}$ can be written explicitly as:

(45) ${\mathbf{𝐑}}_{ij}=\sum _t\sum _n\sum _{\mathrm{\ell }}{\mathbf{𝐖}}_{ni\mathrm{\ell }}{\mathbf{𝐗}}_{n(t+\mathrm{\ell })}\sum _{\tau }{\mathbf{𝐒}}_{t\tau }{\mathbf{𝐇}}_{j\tau }$

Where the index variables $t$ and $\tau$ range from $1$ to $T$, $n$ ranges from $1$ to $N$, and $\mathrm{\ell }$ ranges from $1$ to $L$.

Our goal here is to find a close approximation to this penalty term that does not contain the data $\mathbf{𝐗}$. This can readily be done if $\mathbf{𝐗}$ is well-approximated by the convNMF decomposition:

(46) ${\mathbf{𝐗}}_{nt}\approx {(\mathbf{𝐖}⊛\mathbf{𝐇})}_{nt}=\sum _k\sum _{\mathrm{\ell }}{\mathbf{𝐖}}_{nk\mathrm{\ell }}{\mathbf{𝐇}}_{k(t-\mathrm{\ell })}$

Substituting this expression into Equation 45 and defining the smoothed matrix ${\sum }_{\tau }{\mathbf{𝐒}}_{t\tau }{\mathbf{𝐇}}_{j\tau }$ as ${\mathbf{𝐇}}_{jt}^{\text{smooth}}$ gives:

(47) ${\mathbf{𝐑}}_{ij}\approx \sum _t\sum _n\sum _{\mathrm{\ell }}\sum _k\sum _{{\mathrm{\ell }}^{\prime }}{\mathbf{𝐖}}_{ni\mathrm{\ell }}{\mathbf{𝐖}}_{nk{\mathrm{\ell }}^{\prime }}{\mathbf{𝐇}}_{kt}{\mathbf{𝐇}}_{jt}^{\text{smooth}}$

Making the substitution $u=\mathrm{\ell }-{\mathrm{\ell }}^{\prime }$ gives:

(48) ${\mathbf{𝐑}}_{ij}\approx \sum _t\sum _n\sum _{u=-(L-1)}^{L-1}\sum _k\sum _{{\mathrm{\ell }}^{\prime }}{\mathbf{𝐖}}_{ni({\mathrm{\ell }}^{\prime }+u)}{\mathbf{𝐖}}_{nk{\mathrm{\ell }}^{\prime }}{\mathbf{𝐇}}_{k(t+u)}{\mathbf{𝐇}}_{jt}^{\text{smooth}}$

where in the above expression we have taken $u=\mathrm{\ell }-{\mathrm{\ell }}^{\prime }$ to extend over the full range from $-(L-1)$ to $(L-1)$ under the implicit assumption that $\mathbf{𝐖}$ and $\mathbf{𝐇}$ are zero padded such that values of $\mathbf{𝐖}$ for lag indices outside the range $0$ to $L-1$ and values of $\mathbf{𝐇}$ for time indices outside the range $1$ to $T$ are taken to be zero.

Relabeling ${\mathrm{\ell }}^{\prime }$ as $\mathrm{\ell }$ and gathering terms together yields

(49) ${\mathbf{𝐑}}_{ij}\approx \sum _k\sum _{u=-(L-1)}^{L-1}\left(\sum _n\sum _{\mathrm{\ell }}{\mathbf{𝐖}}_{ni(\mathrm{\ell }+u)}{\mathbf{𝐖}}_{nk\mathrm{\ell }}\right)\left(\sum _t{\mathbf{𝐇}}_{k(t+u)}{\mathbf{𝐇}}_{jt}^{\text{smooth}}\right)$

We note that the above expression contains terms that resemble penalties on orthogonality between two $\mathbf{𝐖}$ factors (first parenthetical) or two $\mathbf{𝐇}$ 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 ${\displaystyle k}$ into three contributions corresponding to $k=i$, $k=j$, and $k\ne i,j$

(50) ${\displaystyle \begin{array}{rl}{\displaystyle {\mathbf{R}}_{ij}\approx }& {\displaystyle \sum _{u=-(L-1)}^{L-1}\left(\sum _n\sum _{\ell }{\mathbf{W}}_{ni(\ell +u)}{\mathbf{W}}_{ni\ell }\right)\left(\sum _t{\mathbf{H}}_{i(t+u)}{\mathbf{H}}_{jt}^{\text{smooth}}\right)+}\\ & {\displaystyle \sum _{u=-(L-1)}^{L-1}\left(\sum _n\sum _{\ell }{\mathbf{W}}_{ni(\ell +u)}{\mathbf{W}}_{nj\ell }\right)\left(\sum _t{\mathbf{H}}_{j(t+u)}{\mathbf{H}}_{jt}^{\text{smooth}}\right)+}\\ & {\displaystyle \sum _{k\ne i,j}\sum _{u=-(L-1)}^{L-1}\left(\sum _n\sum _{\ell }{\mathbf{W}}_{ni(\ell +u)}{\mathbf{W}}_{nk\ell }\right)\left(\sum _t{\mathbf{H}}_{k(t+u)}{\mathbf{H}}_{jt}^{\text{smooth}}\right)}\end{array}}$

The first term above contains, for $u=0$, a simple extension of the ${(i,j)}^{th}$ element of the $\mathbf{𝐇}$ orthogonality condition ${\mathbf{𝐇𝐒𝐇}}^{\top }$. The extension is that the orthogonality is weighted by the power, that is the sum of squared elements, in the $i^{th}$ factor of $\mathbf{𝐖}$ (the apparent lack of symmetry in weighing by the $i^{th}$ rather than the $j^{th}$ factor can be removed by simultaneously considering the term $R_{ji}$, as shown in the Fourier representation of the following section). This weighting has the benefit of applying the penalty on $\mathbf{𝐇}$ orthogonality most strongly to those factors whose corresponding $\mathbf{𝐖}$ components contain the most power. For $u\ne 0$, this orthogonality condition is extended to allow for overlap of time-shifted $\mathbf{𝐇}$’s, with weighting at each time shift by the autocorrelation of the corresponding $\mathbf{𝐖}$ factor. Qualitatively, this enforces that (even in the absence of the smoothing matrix $\mathbf{𝐒}$), $\mathbf{𝐇}$’s that are offset by less than the width of the autocorrelation of the corresponding $\mathbf{𝐖}$’s will have overlapping convolutions with these $\mathbf{𝐖}$’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 $\mathbf{𝐇}$ orthogonality regularization, weighted by the strength of the corresponding $\mathbf{𝐖}$ factors.

The second term above represents a complementary orthogonality condition on the $\mathbf{𝐖}$ components, in which orthogonality in the $i^{th}$ and $j^{th}\mathbf{𝐖}$ factors are weighted by the (smoothed) autocorrelation of the $\mathbf{𝐇}$ factors. For the case in which the $\mathbf{𝐇}$ factors have no time-lagged autocorrelations, this corresponds to a simple weighting of $\mathbf{𝐖}$ orthogonality by the strength of the corresponding $\mathbf{𝐇}$ factors.

Finally, we consider the remaining terms of the cost function, for which $k\ne 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:

(51) $\sum _{k\ne i,j}\sum _{u=-(L-1)}^{L-1}\left(\sum _n\sum _{\mathrm{\ell }}{\mathbf{𝐖}}_{ni(\mathrm{\ell }+u)}{\mathbf{𝐖}}_{nk\mathrm{\ell }}\right)\left(\sum _t{\mathbf{𝐇}}_{k(t+u)}{\mathbf{𝐇}}_{jt}^{\text{smooth}}\right)$

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 $i^{th}$ and ${\displaystyle k^{th}\mathbf{W}}$ components. Likewise, the second expression corresponds to the time-lagged correlation of the (smoothed) $j^{th}$ and ${\displaystyle k^{th}\mathbf{H}}$ components. Thus, this term of ${\mathbf{𝐑}}_{ij}$ contributes whenever there is a factor (${\mathbf{𝐖}}_{\cdot k\cdot },{\mathbf{𝐇}}_{k\cdot }$) that overlaps, at the same time lags, with the $i^{th}$ factor’s $\mathbf{𝐖}$ component and the $j^{th}$ factor’s $\mathbf{𝐇}$ 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.

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 ${\mathbf{𝐑}}_{23}$ due to the overlapping values of the pair (${\mathbf{𝐡}}_{\mathrm{𝟐}},{\mathbf{𝐡}}_{\mathrm{𝟑}}$). Likewise, there is a penalty ${\mathbf{𝐑}}_{13}$ due to the overlapping values of the pair (${\mathbf{𝐰}}_{\mathrm{𝟏}},{\mathbf{𝐰}}_{\mathrm{𝟑}}$). The triplet penalty term contributes to ${\mathbf{𝐑}}_{12}$ and derives from the fact that ${\mathbf{𝐰}}_{\mathrm{𝟏}}$ overlaps with ${\mathbf{𝐰}}_{\mathrm{𝟑}}$ at the same time (and with the same, zero time lag) as ${\mathbf{𝐡}}_{\mathrm{𝟐}}$ overlaps with ${\mathbf{𝐡}}_{\mathrm{𝟑}}$.

In summary, the above analysis shows that for good reconstructions of the data where $\mathbf{𝐗}\approx \mathbf{𝐖}⊛\mathbf{𝐇}$, the x-ortho penalty can be well-approximated by the sum of three contributions. The first corresponds to a penalty on time-lagged (smoothed) $\mathbf{𝐇}$ orthogonality weighted at each time lag by the autocorrelation of the corresponding $\mathbf{𝐖}$ factors. The second similarly corresponds to a penalty on time-lagged $\mathbf{𝐖}$ orthogonality weighted at each time lag by the (smoothed) autocorrelation of the corresponding $\mathbf{𝐇}$ factors. For simple cases of sparse sequences, these contributions reduce to orthogonality in $\mathbf{𝐇}$ or $\mathbf{𝐖}$ weighted by the power in the corresponding $\mathbf{𝐖}$ or $\mathbf{𝐇}$, 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 $\mathbf{𝐖}$ and a different factor in $\mathbf{𝐇}$ both overlap (at the same time lag) with a third common factor, and may occur even when the factors $\mathbf{𝐖}$ and $\mathbf{𝐇}$ 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 $\mathbf{𝐗}$.

Appendix 2—figure 1

Download asset Open asset