Machine learning methods applied to multi-dimensional surface learning pose some fundamental questions on the importance of the mathematical expression of the feature dimensions of the input space. Moreover, for Gaussian processes which are particularly popular in regression problems, the choice of an appropriate kernel function to construct the desired model is not straightforward. The feature spaces, onto which one projects molecular geometries, have many degrees of freedom and do not necessarily consider the physics of the problem. One then has to consider the consequence of the properties of the feature spaces onto the latent functions of the corresponding Gaussian processes. Moreover, for our particular interest of potential energy surface learning, one usually deals with sparsity in the training data given the high computational cost of ab initio electronic structure calculations. This sparsity creates instability in the models that a Bayesian approach to Gaussian processes creates and the expression of the training data on the feature space is thus important. This thesis covers three aspects of the learning process. Firstly, Gaussian processes that project the input molecular geometries onto spaces that have different properties. Secondly, the optimisation of a Morse projected feature space, which is controlled through its projection parameters. The latter can be selected either with a Bayesian optimisation approach or a separate optimisation. Finally, a different noise kernel that treats to move away from the usual homoscedastic noise treatment. The latter emerges from learning surface from stochastic electronic structure methods which provide, as well as energy data, error estimations of the calculations.