  Partial Least Squares (PLS) regression emerged as an alternative to ordinary
least squares for addressing multicollinearity in a wide range of scientific
applications. As multidimensional tensor data is becoming more widespread,
tensor adaptations of PLS have been developed. Our investigations reveal that
the previously established asymptotic result of the PLS estimator for a tensor
response breaks down as the tensor dimensions and the number of features
increase relative to the sample size. To address this, we propose Sparse Higher
Order Partial Least Squares (SHOPS) regression and an accompanying algorithm.
SHOPS simultaneously accommodates variable selection, dimension reduction, and
tensor association denoising. We establish the asymptotic accuracy of the SHOPS
algorithm under a high-dimensional regime and verify these results through
comprehensive simulation experiments, and applications to two contemporary
high-dimensional biological data analysis.
