We present an approach combining a representation of a multivariate function using subdimensional functions with machine learning based representation of component functions: Random sampling high dimensional model representation Gaussian process regression (RS-HDMR-GPR). The use of Gaussian process regressions to represent component functions allows nonparametric (unbiased) representation and the possibility to work only with functions of desired dimensionality, obviating the need to build an expansion over orders of coupling. All component functions are determined from a single set of samples. The method is tested by fitting six- and 15-dimensional potential energy surfaces (PES) of polyatomic molecules as well as by computing vibrational spectra for a four-atomic molecule.