The multidimensional nature of neuroscience data has made the use of multi-way statistical analysis suitable in this field. Parallel Factor Analysis (PARAFAC) is a multidimensional generalization of PCA with the advantage of offering unique solutions. However, imposing physiologically acceptable constraints would improve the interpretation of this type of analysis. In this work we propose a new algorithm called Alternating Penalized Least Squares to estimate PARAFAC solutions using different kinds of soft penalization. The algorithm relies on the recent generalization of modified Newton-Raphson techniques to estimate a multiple penalized least squares model. Applied to semi-synthetic and real spontaneous EEG time-varying spectra, we show that a wide range of sparse and smooth solutions can be found separately, as well as with these two properties combined. Smoothness is usually desired in spectra, and different sparse scenarios are observed in the temporal evolution of physiological intermittent phenomena. The degree of constraints can be tuned through the weighting parameters, whose optimal values can be chosen by means of the cross-validation and Corcondia measures.