Activity of populations of sensory neurons carries stimulus information in both the temporal and the spatial dimensions. This poses the question of how to compactly represent all the information that the population codes carry across all these dimensions. Here, we developed an analytical method to factorize a large number of retinal ganglion cells' spike trains into a robust low-dimensional representation that captures efficiently both their spatial and temporal information. In particular, we extended previously used single-trial space-by-time tensor decomposition based on non-negative matrix factorization to efficiently discount pre-stimulus baseline activity. On data recorded from retinal ganglion cells with strong pre-stimulus baseline, we showed that in situations where the stimulus elicits a strong change in firing rate, our extensions yield a boost in stimulus decoding performance. Our results thus suggest that taking into account the baseline can be important for finding a compact information-rich representation of neural activity.
This study proposes the use of non-negative matrix factorization accounting for baseline by subtracting the pre-stimulus baseline from each trial and subsequently decompose the data using a 3-way factorization thereby identifying spatial and temporal modules as well as their signed activation. The method is used on data recorded from mouse and pig retinal ganglion cells of time binned spike trains providing improved performance over non-baseline corrected data. 
Pros:
The paper is well written, the analysis interesting and the application of the Tucker2 framework sound. Removing baseline is a reasonable step and the paper includes analysis of several spike-train datasets. The analysis of the approaches in terms of their ability to decode is also sound and interesting.
Cons:
I find the novelty of the paper limited: 
The authors extend the work by (Onken et al. 2016) to subtract baseline (a rather marginal innovation) of this approach. To use a semi-NMF type of update rule (as proposed by Ding et al .2010) and apply the approach to new spike-train datasets evaluating performance by their decoding ability (decoding also considered in Onken et al. 2016).
Multiplicative update-rules are known to suffer from slow-convergence and I would suspect this also to be an issue for the semi-NMF update rules. It would therefore be relevant and quite easy to consider other approaches such as active set or column wise updating also denoted HALS which admit negative values in the optimization, see also the review by N. Giles
https://arxiv.org/abs/1401.5226
as well as for instance:
Nielsen, Søren Føns Vind, and Morten Mørup. "Non-negative tensor factorization with missing data for the modeling of gene expressions in the human brain." Machine Learning for Signal Processing (MLSP), 2014 IEEE International Workshop on. IEEE, 2014.
It would improve the paper to also discuss that the non-negativity constrained Tucker2 model may be subject to local minima solutions and have issues of non-uniqueness (i.e. rotational ambiguity). At least local minima issues could be assessed using multiple random initializations.
The results are in general only marginally improved by the baseline corrected non-negativity constrained approach. For comparison the existing methods ICA, Tucker2 should also be evaluated for the baseline corrected data, to see if it is the constrained representation or the preprocessing influencing the performance. Finally, how performance is influenced by dimensionality P and L should also be clarified.
It seems that it would be naturally to model the baseline by including mean values in the model rather than treating the baseline as a preprocessing step. This would bridge the entire framework as one model and make it potentially possible to avoid structure well represented by the Tucker2 representation to be removed by the preprocessing.
Minor: 
The approach corresponds to a Tucker2 decomposition with non-negativity constrained factor matrices and unconstrained core - please clarify this as you also compare to Tucker2 in the paper with orthogonal factor matrices.
Ding et al. in their semi-NMF work provide elaborate derivation with convergence guarantees.  In the present paper these details are omitted and it is unclear how the update rules are derived from the KKT conditions and the Lagrange multiplier and how they differ from standard semi-NMF, this should be better clarified.