  We present a novel scheme to accurately predict atomic forces as vector
quantities, rather than sets of scalar components, by Gaussian Process (GP)
Regression. This is based on matrix-valued kernel functions, on which we impose
the requirements that the predicted force rotates with the target configuration
and is independent of any rotations applied to the configuration database
entries. We show that such covariant GP kernels can be obtained by integration
over the elements of the rotation group SO(d) for the relevant dimensionality
d. Remarkably, in specific cases the integration can be carried out
analytically and yields a conservative force field that can be recast into a
pair interaction form. Finally, we show that restricting the integration to a
summation over the elements of a finite point group relevant to the target
system is sufficient to recover an accurate GP. The accuracy of our kernels in
predicting quantum-mechanical forces in real materials is investigated by tests
on pure and defective Ni, Fe and Si crystalline systems.
