This thesis introduces a framework that is able to describe general many-body coarse-grained interactions. We make use of this to describe the free energy surface as a cluster expansion in terms of monomer, dimer, and trimer terms. The contributions to the free energy due to these terms are inferred from MD results of the underlying all-atom model using Gaussian Approximation Potentials, a type of machine-learning potential based on Gaussian process regression. This provides CG interactions that are much more accurate than is possible with site-based pair potentials. While slower than these, it can still be faster than all-atom simulations for solvent-free CG models of systems with a large amount of solvent, as is common in biomolecular simulations.