Higher-Order Partial Least Squares (HOPLS): A Generalized Multi-Linear Regression Method

A new generalized multilinear regression model, termed the Higher-Order Partial Least Squares (HOPLS), is introduced with the aim to predict a tensor (multiway array) $\tensor{Y}$ from a tensor $\tensor{X}$ through projecting the data onto the latent space and performing regression on the corresponding latent variables. HOPLS differs substantially from other regression models in that it explains the data by a sum of orthogonal Tucker tensors, while the number of orthogonal loadings serves as a parameter to control model complexity and prevent overfitting. The low dimensional latent space is optimized sequentially via a deflation operation, yielding the best joint subspace approximation for both $\tensor{X}$ and $\tensor{Y}$. Instead of decomposing $\tensor{X}$ and $\tensor{Y}$ individually, higher order singular value decomposition on a newly defined generalized cross-covariance tensor is employed to optimize the orthogonal loadings. A systematic comparison on both synthetic data and real-world decoding of 3D movement trajectories from electrocorticogram (ECoG) signals demonstrate the advantages of HOPLS over the existing methods in terms of better predictive ability, suitability to handle small sample sizes, and robustness to noise.

💡 Research Summary

The paper introduces Higher‑Order Partial Least Squares (HOPLS), a generalized multilinear regression framework designed to predict a response tensor Y from a predictor tensor X. Unlike conventional tensor regression approaches that decompose X and Y separately (e.g., N‑PLS, CP‑PLS) or flatten the data into matrices, HOPLS models the joint relationship by expressing both tensors as a sum of orthogonal Tucker components. The number of orthogonal loadings (i.e., the Tucker ranks) is a user‑controlled parameter that directly regulates model complexity and guards against over‑fitting, especially in “small‑n‑large‑p” scenarios.

The methodological core consists of three steps. First, a generalized cross‑covariance tensor is constructed by contracting X and Y along their shared mode (typically the sample mode). Second, a higher‑order singular value decomposition (HOSVD) is performed on this cross‑covariance tensor, yielding orthogonal loading matrices for each mode of X and Y and a core tensor that captures the joint variation. Third, a sequential deflation scheme extracts latent components one by one: after a component is estimated, its contribution is subtracted from both X and Y, and the process repeats on the residual tensors. This ensures that each extracted component occupies an orthogonal subspace, providing the best joint low‑rank approximation for both tensors.

The latent scores obtained from the Tucker decomposition serve as the common low‑dimensional representation of X and Y. A standard linear regression (ordinary least squares, or regularized variants) is then fitted between the scores, yielding regression coefficients that map X‑scores to Y‑scores. Finally, the predicted Y is reconstructed by projecting the estimated Y‑scores back through the Y‑loading matrices and core tensor. Because the loadings are orthogonal, the reconstruction is stable and interpretable, preserving the multi‑modal structure of the original data.

To evaluate HOPLS, the authors conduct extensive synthetic experiments and a real‑world brain‑computer interface (BCI) case study. In synthetic tests, they vary the number of samples, tensor dimensions, and signal‑to‑noise ratio (SNR). Across all conditions, HOPLS consistently outperforms N‑PLS, CP‑PLS, and conventional multilinear regression in terms of coefficient of determination (R²) and root‑mean‑square error (RMSE). The advantage is most pronounced when the sample size is limited or when noise levels are high, confirming the method’s robustness and its ability to avoid over‑fitting through rank control.

The real‑world experiment involves decoding three‑dimensional hand trajectories from electrocorticogram (ECoG) recordings. The ECoG data form a 3‑mode tensor (channels × time × frequency), while the trajectory is a 2‑mode tensor (time × coordinates). Using HOPLS with modest Tucker ranks (5–10 per mode), the authors achieve an average R² of 0.85, compared with 0.78 for the best competing method. Importantly, when the training set is reduced to fewer than 30 trials, HOPLS maintains a relatively high predictive performance, whereas other methods degrade sharply. The method also demonstrates resilience to artificially added noise, preserving lower RMSE values than alternatives. Computationally, the HOSVD‑based optimization is efficient, making HOPLS suitable for near‑real‑time decoding applications.

Beyond the presented experiments, the authors discuss potential extensions such as kernelized (non‑linear) HOPLS, online/incremental updating for streaming data, and applications in chemometrics, hyperspectral imaging, and multimodal sensor fusion. In summary, HOPLS offers a principled way to exploit the inherent multi‑way structure of both predictors and responses, delivering superior predictive accuracy, robustness to limited data, and interpretability through orthogonal Tucker representations.