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 ³A″ state of SH₂, which facilitates the SH + H ↔ S(³P) + H₂ 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 × 10²) of ab initio points, but it needs substantially more points (∼1 × 10³) to converge reaction probabilities. The implications of these observations for construction of potential energy surfaces are discussed.