The adsorption energies of molecular adsorbates on catalyst surfaces are key descriptors in computational catalysis research. For the relatively large reaction intermediates frequently encountered, e.g., in syngas conversion, a multitude of possible binding motifs leads to complex potential energy surfaces (PES), however. This implies that finding the optimal structure is a difficult global optimization problem, which leads to significant uncertainty about the stability of many intermediates. To tackle this issue, we present a global optimization protocol for surface adsorbate geometries which trains a surrogate machine learning potential on-the-fly. The approach is applicable to arbitrary surface models and adsorbates and minimizes both human intervention and the number of required DFT calculations by iteratively updating the training set with configurations explored by the algorithm. We demonstrate the efficiency of this approach for a diverse set of adsorbates on the Rh(111) and (211) surfaces.