Here we introduce a new approach to compute the finite temperature lattice dynamics from first-principles via the newly developed slave mode expansion. We study PbTe where inelastic neutron scattering (INS) reveals strong signatures of nonlinearity as evidenced by anomalous features which emerge in the phonon spectra at finite temperature. Using our slave mode expansion in the classical limit, we compute the vibrational spectra and show remarkable agreement with temperature dependent INS measurements. Furthermore, we resolve experimental controversy by showing that there are no appreciable local nor global spontaneously broken symmetries at finite temperature and that the anomalous spectral features simply arise from two anharmonic interactions. Our approach should be broadly applicable across the periodic table.Inelastic neutron scattering (INS) is a fundamental probe of materials which has allowed a unique view into the most basic aspects of mechanical behavior [1]. With the advent of the Spallation Neutron Source, massive amounts of data can be accumulated at an unprecedented rate, allowing for an extremely detailed inspection of phenomena throughout reciprocal space. Despite the enormous success of this experimental method, theory has greatly lagged behind in the context of vibrations. At present, one cannot compute the temperature dependence of the phonon spectrum for materials with appreciable phonon interactions from first-principles. Moreover, this scenario cannot even be routinely handled in the classical limit for indirect reasons. This assertion is empirically well illustrated by the strongly interacting phonon material PbTe. INS measurements demonstrated signatures of strong interactions in the temperature dependence of the phonon spectra [2,3], attracting significant attention to this system, yet complementary theoretical predictions have not yet been made. In this paper, we circumvent previous theoretical limitations, and resolve the experimental anomalies in PbTe.In many materials, density functional theory (DFT) is expected to describe the structural energetics to a high degree of accuracy. Therefore, the interatomic potential generated via solving the DFT equations (assuming the Born-Oppenheimer approximation) should be reliable to perform both quantum and classical dynamics of the nuclei. However, the problem is that the scaling of DFT is sufficiently prohibitive to forbid this in many scenarios. Even at the level of classical mechanics (aka. ab initio molecular dynamics), one would need a very large unit cell along with many time steps. In most cases this prevents one from directly being able to compute even the classical vibrational spectrum at finite temperature, as evidenced by the sparse number of such publications in the literature. There are a variety of approaches that attempt to circumvent this problem. One approach is the use of an empirical potential in place of DFT. If an accurate empirical potential exists this is acceptable, but this is generally not the case. Another approach wou...