The quaternion kernel least squares algorithm (QKLS) is introduced as a generic kernel framework for the estimation of multivariate quaternion valued signals. This is achieved based on the concepts of quaternion inner product and quaternion positive definiteness, allowing us to define quaternion kernel regression. Next, the least squares solution is derived using the recently introduced HR calculus. We also show that QKLS is a generic extension of standard kernel least squares, and their equivalence is established for real valued kernels. The superiority of the quaternion-valued linear kernel with respect to its real-valued counterpart is illustrated for both synthetic and real-world prediction applications, in terms of accuracy and robustness to overfitting.