High-resolution spectroscopic surveys of stars in the Milky Way have entered the Big Data regime and have opened avenues for solving outstanding questions in the history of galaxies. However, exploiting their full potential is limited by complex systematics, whose characterization has not received much attention in modern spectroscopic analyses. In this work, we present a novel method based on Functional Principal Component Analysis that disentangles the component of spectral data space intrinsic to the stars from that due to systematics. This method also reduces the dimensionality of spectra and imputes masked wavelengths, thereby enabling accurate studies of populations of stars. To demonstrate the applicability of our results, we use them to infer chemical compositions of stars in the open cluster M67. We employ Sequential Neural Likelihood, a simulation-based Bayesian inference method that learns likelihood functions using neural density estimators, to incorporate non-Gaussian effects in spectral likelihoods. We discuss this application, and its implications for understanding the formation and evolution of the Milky Way.