The development of reliable, yet computationally efficient interatomic forcefields is key to facilitate the modeling of glasses. However, the parameterization of novel forcefields is challenging as the high number of parameters renders traditional optimization methods inefficient 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 forcefield that offers an unprecedented agreement with ab initio simulations. This method offers a new route to efficiently parameterize new interatomic forcefields for disordered solids in a non-biased fashion.