  High-resolution spectroscopic surveys of the Milky Way have entered the Big
Data regime and have opened avenues for solving outstanding questions in
Galactic archaeology. 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 to
disentangle the component of spectral data space intrinsic to the stars from
that due to systematics. Using functional principal component analysis on a
sample of $18,933$ giant spectra from APOGEE, we find that the intrinsic
structure above the level of observational uncertainties requires ${\approx}$10
functional principal components (FPCs). Our FPCs can reduce the dimensionality
of spectra, remove systematics, and impute masked wavelengths, thereby enabling
accurate studies of stellar populations. To demonstrate the applicability of
our FPCs, we use them to infer stellar parameters and abundances of 28 giants
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. By hierarchically combining the inferred abundances, we
limit the spread of the following elements in M67: $\mathrm{Fe} \lesssim 0.02$
dex; $\mathrm{C} \lesssim 0.03$ dex; $\mathrm{O}, \mathrm{Mg}, \mathrm{Si},
\mathrm{Ni} \lesssim 0.04$ dex; $\mathrm{Ca} \lesssim 0.05$ dex; $\mathrm{N},
\mathrm{Al} \lesssim 0.07$ dex (at 68% confidence). Our constraints suggest a
lack of self-pollution by core-collapse supernovae in M67, which has promising
implications for the future of chemical tagging to understand the star
formation history and dynamical evolution of the Milky Way.
