Abstract Quantum computing algorithms have been shown to produce performant quantum kernels for machine-learning classification problems. Here, we examine the performance of quantum kernels for regression problems of practical interest. For an unbiased benchmarking of quantum kernels, it is necessary to construct the most optimal functional form of the classical kernels and the most optimal quantum kernels for each given data set. We develop an algorithm that uses an analog of the Bayesian information criterion to optimize the sequence of quantum gates used to estimate quantum kernels for Gaussian process models. The algorithm increases the complexity of the quantum circuits incrementally, while improving the performance of the resulting kernels, and is shown to yield much higher model accuracy with fewer quantum gates than a fixed quantum circuit ansatz. We demonstrate that quantum kernels thus obtained can be used to build accurate models of global potential energy surfaces (PES) for polyatomic molecules. The average interpolation error of the six-dimensional PES obtained with a random distribution of 2000 energy points is 16 cm−1 for H3O+, 15 cm−1 for H2CO and 88 cm−1 for HNO2. We show that a compositional optimization of classical kernels for Gaussian process regression converges to the same errors. This indicates that quantum kernels can achieve the same, though not better, expressivity as classical kernels for regression problems.