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.