Surface phase diagrams are necessary for understanding surface chemistry in electrochemical catalysis, where a range of adsorbates and coverages exist at varying applied potentials. These diagrams are typically constructed using intuition, which risks missing complex coverages and configurations at potentials of interest. More accurate cluster expansion methods are often difficult to implement quickly for new surfaces. We adopt a machine learning approach to rectify both issues. Using a Gaussian process regression model, the free energy of all possible adsorbate coverages for surfaces is predicted for a finite number of adsorption sites. Our result demonstrates a rational, simple, and systematic approach for generating accurate free-energy diagrams with reduced computational resources. The Pourbaix diagram for the IrO2(110) surface (with nine coverages from fully hydrogenated to fully oxygenated surfaces) is reconstructed using just 20 electronic structure relaxations, compared to approximately 90 using typical search methods. Similar efficiency is demonstrated for the MoS2 surface.