We refer readers to recent reviews for a comprehensive list of developments in this rapidly growing field. In this study, we focus on GAP models, which typically require less data to be trained than neural network potentials, and have good scalability and computational efficiency for large-scale molecular dynamics simulations. In a previous study, we demonstrated that GAP models trained with DFT calculations 9 provide accurate estimates of the thermal expansion properties of graphene. Other studies have also shown the success of GAP models for thermal properties of 2DMs, such as graphene, 10 carbon allotropes, 11 monolayer h-BN, 12 h-BN allotropes, 13 silicene, and monolayer MoS 2 . In this work, we propose a streamlined procedure for generating highly accurate GAP models for 2DMs. We demonstrate our approach for graphene, silicene (as a buckled 2DM), and h-XN (X = B, Al, and Ga, as binary 2DMs). In addition to traditional validation using randomly selected data points (excluded from the training set) from ab initio molecular dynamics (AIMD) trajectories, we also use thermal property calculations to validate the accuracy of the MLIPs. The first test is the phonon dispersion curves, which provide information about dynamical stability and two-phonon interactions through 2 nd order interatomic force constants (IFCs). The second test is the iterative solution of lattice thermal conductivity, which also depends on 3 rd order IFCs and provides information about the anharmonicity of the structure. The accuracy of force prediction for 2 nd and 3 rd order force constant calculations is also tested by implementing the HIPHIVE package. The HIPHIVE package is a Python library that aims to generate high-order IFCs with reduced number of force calculations by using advanced optimation and machine learning techniques. Therefore, reference forces (i.e., training data) must be calculated with high (first principles level) accuracy to obtain a reliable output. Finally, the performance of large-scale simulations is confirmed via phonon density of states computed with velocity auto-correlation calculations based on molecular dynamics simulations.