  Many problems in contemporary astrophysics---from understanding the formation
of black holes to untangling the chemical evolution of galaxies---rely on
knowledge about binary stars. This, in turn, depends on discovery and
characterization of binary companions for large numbers of different kinds of
stars in different chemical and dynamical environments. Current stellar
spectroscopic surveys observe hundreds of thousands to millions of stars with
(typically) few observational epochs, which allows binary discovery but makes
orbital characterization challenging. We use a custom Monte Carlo sampler (The
Joker) to perform discovery and characterization of binary systems through
radial-velocities, in the regime of sparse, noisy, and poorly sampled
multi-epoch data. We use it to generate posterior samplings in Keplerian
parameters for 232,531 sources released in APOGEE Data Release 16. Our final
catalog contains 19,635 high-confidence close-binary (P < few years, a < few
AU) systems that show interesting relationships between binary occurrence rate
and location in the color-magnitude diagram. We find notable faint companions
at high masses (black-hole candidates), at low masses (substellar candidates),
and at very close separations (mass-transfer candidates). We also use the
posterior samplings in a (toy) hierarchical inference to measure the
long-period binary-star eccentricity distribution. We release the full set of
posterior samplings for the entire parent sample of 232,531 stars. This set of
samplings involves no heuristic "discovery" threshold and therefore can be used
for myriad statistical purposes, including hierarchical inferences about
binary-star populations and sub-threshold searches.
