The shape of earthquake source spectra, traditionally fit by physics-based models, contains important parameters to constrain rupture dimension, duration, and geometry. Here we apply machine learning (ML) to derive a data-driven model of source spectra from a global set of 3675 Mw>5.5 earthquakes. The ML model, in the same degree of freedom as the Brune model, improves the goodness of fit by 8.5%. Specifically, the ML model fits the data without systematic bias, whereas Brune model tends to underestimate at intermediate frequencies and overestimate at high frequencies. The latter descrepancy cannot be modelled by intrinsic attenuation effect, nor by increasing the fall-off exponent in Brune-type or Boatwright-type models. Instead, it could be matched by a Haskell-type double-corner-frequency model with a high-frequency falloff of f ^-2.7. Our results demonstrate that unsupervised machine learning can extract hidden global characteristics of high-dimensional data, which provide hints to amend existing physical models.