Simulations of colloidal suspensions consisting of mesoscopic particles and smaller species such as ions or depletants are computationally challenging as different length and time scales are involved. Here, 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 we fit 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 effectively state-independent and can be used in direct-coexistence simulations. We show that our ML method reduces the computational cost by several 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.