Simulations of colloidal suspensions consisting of mesoscopic particles and smaller species like ions or depletants are computationally intractable as different length and time scales are involved.
We introduce a machine learning (ML) approach in which the degrees of freedom of the microscopic species are integrated out and the mesoscopic particles interact with effective many-body potentials, which are fitted as a function of all colloid coordinates with a set of symmetry functions. We apply this approach to a colloid-polymer mixture. Remarkably, the ML potentials can be assumed to be state-independent and can be used in direct-coexistence simulations. We show that our ML method reduces the computational cost by up to four orders of magnitude compared to a numerical evaluation, and accurately describes the phase behavior and structure, even for state points where the effective potential is largely determined by many-body contributions.