  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 B$_x$N$_y$O$_z$ 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 B$_{x+2}$N$_{x}$O$_{3}$ compounds with $x \geq 3$
(like B$_5$N$_3$O$_3$, B$_6$N$_4$O$_3$, etc.) are potentially superhard and
thermodynamically favorable. Our meta-GGA density functional theory
calculations indicate that these materials are also wide bandgap ($\ge 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.
