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 (3)A '' state of SH2, which facilitates the SH + H <-> S(P-3) + H-2 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 (similar to 1 X 10(2)) of ab initio points, but it needs substantially more points (similar to 1 x 10(3)) to converge reaction probabilities. The implications of these observations for construction of potential energy surfaces are discussed.