  Gaussian process regression has recently emerged as a powerful,
system-agnostic tool for building global potential energy surfaces (PES) of
polyatomic molecules. While the accuracy of GP models of PES increases with the
number of potential energy points, so does the numerical difficulty of training
and evaluating GP models. Here, we demonstrate an approach to improve the
accuracy of global PES without increasing the number of energy points. The
present work reports four important results. First, we show that the selection
of the best kernel function for GP models of PES can be automated using the
Bayesian information criterion as a model selection metric. Second, we
demonstrate that GP models of PES trained by a small number of energy points
can be significantly improved by iteratively increasing the complexity of GP
kernels. The composite kernels thus obtained maximize the accuracy of GP models
for a given distribution of potential energy points. Third, we show that the
accuracy of the GP models of PES with composite kernels can be further improved
by varying the training point distributions. Fourth, we show that GP models
with composite kernels can be used for physical extrapolation of PES. We
illustrate the approach by constructing the six-dimensional PES for H$_3$O$^+$.
For the interpolation problem, we show that this algorithm produces a global
six-dimensional PES for H$_3$O$^+$ in the energy range between zero and
$21,000$ cm$^{-1}$ with the root mean square error $65.8$ cm$^{-1}$ using only
500 randomly selected {\it ab initio} points as input. To illustrate
extrapolation, we produce the PES at high energies using the energy points at
low energies as input. We show that one can obtain an accurate global fit of
the PES extending to $21,000$ cm$^{-1}$ based on 1500 potential energy points
at energies below $10,000$ cm$^{-1}$.
