White dwarfs (WDs) are powerful tools to study the evolutionary history of stars and binaries in the Galaxy. But do we understand their multiplicity from a theoretical point of view? This can be tested by a comparison with the sample of WDs within 20 pc, which is minimally affected by selection biases. From the literature, we compile the available information of the local WD sample with a particular emphasis on their multiplicity, and compare this to synthetic models of WD formation in single stars and binaries. As part of our population synthesis approach, we also study the effect of different assumptions concerning the star formation history, binary evolution, and the initial distributions of binary parameters. We find that the observed space densities of single and binary WDs are well reproduced by the models. The space densities of the most common WD systems (single WDs and unresolved WD-MS binaries) are consistent within a factor two with the observed value. We find a discrepancy only for the space density of resolved double WDs. We exclude that observational selection effects, fast stellar winds, or dynamical interactions with other objects in the Milky Way explain this discrepancy. We find that either the initial mass ratio distribution in the Solar neighbourhood is biased towards low mass-ratios, or more than ten resolved DWDs have been missed observationally in the 20 pc sample. Furthermore, we show that the low binary fraction of WD systems (~25%) compared to Solartype MS-MS binaries (~50%) is consistent with theory, and is mainly caused by mergers in binary systems, and to a lesser degree by WDs hiding in the glare of their companion stars. Lastly, Gaia will dramatically increase the size of the volume-limited WD sample, detecting the coolest and oldest WDs out to 50 pc. We provide a detailed estimate of the number of single and binary WDs in the Gaia sample.