We introduce a computational framework that is able to describe general many-body coarse-grained (CG) interactions of molecules and use it to model the free energy surface of molecular liquids 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 all-atom molecular dynamics (MD) data using Gaussian Approximation Potentials, a type of machine-learning model that employs Gaussian process regression. The resulting CG model is much more accurate than those possible using pair potentials. Though slower than the latter, our model can still be faster than all-atom simulations for solvent-free CG models commonly used in biomolecular simulations.