Face super-resolution (FSR) is an effective way to solve the problem of inferring a high-resolution (HR) image from one or more low-resolution (LR) face images. In this paper, we propose a new FSR method to reconstruct the HR face images, which utilizes covariation-guided orthonormalized partial least squares to learn the coherent subspace features from LR and high-frequency (HF) facial image patch data. The proposed method models the nonlinear association between different LR and HF facial patch features. This makes it capable to leverage the rich information of features for effective FSR reconstruction. Experimental results on the public CAS-PEAL-R1 and CelebA face image datasets show that our proposed method is able to achieve better FSR performance than the related super-resolution methods.