The thermal properties of β-Ga2O3 can significantly affect the performance and reliability of high-power electronic devices. To date, due to the absence of a reliable interatomic potential, first-principles calculations based on density functional theory (DFT) have been routinely used to probe the thermal properties of β-Ga2O3. DFT calculations can only tackle small-scale systems due to the huge computational cost, while the thermal transport processes are usually associated with large time and length scales. In this work, we develop a machine learning based Gaussian approximation potential (GAP) for accurately describing the lattice dynamics of perfect crystalline β-Ga2O3 and accelerating atomic-scale simulations. The GAP model shows excellent convergence, which can faithfully reproduce the DFT potential energy surface at a training data size of 32 000 local atomic environments. The GAP model is then used to predict ground-state lattice parameters, coefficients of thermal expansion, heat capacity, phonon dispersions at 0 K, and anisotropic thermal conductivity of β-Ga2O3, which are all in excellent agreement with either the DFT results or experiments. The accurate predictions of phonon dispersions and thermal conductivities demonstrate that the GAP model can well describe the harmonic and anharmonic interactions of phonons. Additionally, the successful application of our GAP model to the phonon density of states of a 2500-atom β-Ga2O3 structure at elevated temperature indicates the strength of machine learning potentials to tackle large-scale atomic systems in long molecular simulations, which would be almost impossible to generate with DFT-based molecular simulations at present.