We search for new superhard B-N-O compounds with an iterative machine learning (ML) procedure, where ML models are trained using sample crystal structures from evolutionary algorithm. We first use cohesive energy to evaluate the thermodynamic stability of varying BxNyOz compositions, and then gradually focus on compositional regions with high cohesive energy and high hardness. The results converge quickly after a few iterations. Our resulting ML models show that Bx+2NxO3 compounds with x ≥ 3 (like B5N3O3, B6N4O3, etc.) are potentially superhard and thermodynamically favorable. Our meta-GGA density functional theory calculations indicate that these materials are also wide bandgap (≥ 4.4 eV) insulators, with the valence band maximum related to the p-orbitals of nitrogen atoms near vacant sites. This study demonstrates that an iterative method combining ML and ab initio simulations provides a powerful tool for discovering novel materials.