Representation of multidimensional global potential energy surfaces suitable for spectral and dynamical calculations from high-level ab initio calculations remains a challenge. Here, we present a detailed study on constructing potential energy surfaces using a machine learning method, namely, Gaussian process regression. Tests for the 3A″ state of SH2, which facilitates the SH + H ↔ S(3P) + H2 abstraction reaction and the SH + H' ↔ SH' + H exchange reaction, suggest that the Gaussian process is capable of providing a reasonable potential energy surface with a small number (∼1 × 102) of ab initio points, but it needs substantially more points (∼1 × 103) to converge reaction probabilities. The implications of these observations for construction of potential energy surfaces are discussed.