Multiway array decomposition has been successful in providing a better understanding of the structure underlying data and in discovering potentially hidden feature dependences serving high-performance decoder applications. However, the computational cost of multiway algorithms can become prohibitive, especially when considering large datasets, rendering them unsuitable for time-critical applications. We propose a multiway regression model for large-scale tensors with optimized performance in terms of time complexity, called fast higher order partial least squares (fHOPLS). We compare fHOPLS with its native version, higher order partial least squares (HOPLS), the state-of-the-art in multilinear regression, under different noise conditions and tensor dimensionalities using synthetic data. We also compare their performance when used for predicting scalp-recorded electroencephalography signals from invasively recorded electrocorticography signals in an oddball experiment. For the sake of exposition, we evaluated the performance of standard unfolded partial least squares (PLS) and linear regression. Our results show that fHOPLS is significantly faster than HOPLS, in particular for big data. In addition, the regression performances of fHOPLS and HOPLS are comparable and outperform both unfolded PLS and linear regression. Another interesting result is that multiway array decoding yields more accurate results than epoch-based averaging procedures traditionally used in the brain computer interfacing community.