The development of reliable, yet computationally efﬁcient interatomic forceﬁelds is key to facilitate the modeling of glasses. However, the parameterization of novel forceﬁelds is challenging as the high number of parameters renders traditional optimization methods inefﬁcient or subject to bias. Here, we present a new parameterization method based on machine learning, which combines ab initio molecular dynamics simulations and Bayesian optimization. By taking the example of glassy silica, we show that our method yields a new interatomic forceﬁeld that offers an unprecedented agreement with ab initio simulations. This method offers a new route to efﬁciently parameterize new interatomic forceﬁelds for disordered solids in a non-biased fashion.