We present an approach to generate machine-learned force fields (MLFF) with beyond density functional theory (DFT) accuracy. Our approach combines on-the-fly active learning and ∆-machine learning in order to generate an MLFF for zirconia based on the random phase approximation (RPA). Specifically, an MLFF trained on-the-fly during DFT based molecular dynamics simulations is corrected by another MLFF that is trained on the differences between RPA and DFT calculated energies, forces and stress tensors. Thanks to the relatively smooth nature of the differences, the expensive RPA calculations are performed only on a small number of representative structures of small unit cells. These structures are determined by a singular value decomposition rank compression of the kernel matrix with low spatial resolution. This dramatically reduces the computational cost and allows us to generate an MLFF fully capable of reproducing high-level quantum-mechanical calculations beyond DFT. We carefully validate our approach and demonstrate its success in studying the phase transitions of zirconia.Machine learning based regression techniques have become a prominent tool to construct accurate interatomic potentials for materials modeling and simulations [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Machinelearned force fields (MLFF), however, are generally constructed by fitting the energies, forces, and stress tensors derived by density functional theory (DFT) calculations, and therefore the accuracy of the resulting MLFFs is largely limited by DFT. It is not surprising, then, that these MLFFs would fail in problems where DFT is inaccurate, such as in systems where long-range electronic correlation effects play an important role. This implies that pursuing an MLFF beyond DFT is highly desirable.