Abstract Upcoming million-star spectroscopic surveys have the potential to revolutionize our view of the formation and chemical evolution of the Milky Way. Realizing this potential requires automated approaches to optimize estimates of stellar properties, such as chemical element abundances, from the spectra. The sheer volume and quality of the observations strongly motivate that these approaches should be driven by the data. With this in mind, we introduce SSSpaNG: a data-driven non-Gaussian Process model of stellar spectra. We demonstrate the capabilities of SSSpaNG using a sample of APOGEE red clump stars, whose model parameters we infer using Gibbs sampling. By pooling information between stars to infer their covariance, we permit clear identification of the correlations between spectral pixels. Harnessing this correlation structure, we infer the true spectrum of each red clump star, inpainting missing regions and denoising by a factor of at least two for stars with signal-to-noise of ∼20. As we marginalize over the covariance matrix of the spectra, the effective prior on these true spectra is non-Gaussian and sparsifying, favouring typically small but occasionally large excursions from the mean. The high-fidelity inferred spectra produced with our approach will enable improved chemical elemental abundance estimates for individual stars. Our model also allows us to quantify the information gained by observing portions of a star’s spectrum, and thereby define the most mutually informative spectral regions. Using 25 windows centred on elemental absorption lines, we demonstrate that the iron-peak and alpha-process elements are particularly mutually informative for these spectra, and that the majority of information about a target window is contained in the 10-or-so most informative windows. Such mutual information estimates have the potential to inform models of nucleosynthetic yields and the design of future observations. Our code is made publicly available at https://github.com/sfeeney/ddspectra.