  We develop a data-driven spectral model for identifying and characterizing
spatially unresolved multiple-star systems and apply it to APOGEE DR13 spectra
of main-sequence stars. Binaries and triples are identified as targets whose
spectra can be significantly better fit by a superposition of two or three
model spectra, drawn from the same isochrone, than any single-star model. From
an initial sample of $\sim$20,000 main-sequence targets, we identify
$\sim$2,500 binaries in which both the primary and secondary star contribute
detectably to the spectrum, simultaneously fitting for the velocities and
stellar parameters of both components. We additionally identify and fit
$\sim$200 triple systems, as well as $\sim$700 velocity-variable systems in
which the secondary does not contribute detectably to the spectrum. Our model
simplifies the process of simultaneously fitting single- or multi-epoch spectra
with composite models and does not depend on a velocity offset between the two
components of a binary, making it sensitive to traditionally undetectable
systems with periods of hundreds or thousands of years. In agreement with
conventional expectations, almost all the spectrally-identified binaries with
measured parallaxes fall above the main sequence in the color-magnitude
diagram. We find excellent agreement between spectrally and dynamically
inferred mass ratios for the $\sim$600 binaries in which a dynamical mass ratio
can be measured from multi-epoch radial velocities. We obtain full orbital
solutions for 64 systems, including 14 close binaries within hierarchical
triples. We make available catalogs of stellar parameters, abundances, mass
ratios, and orbital parameters.
