We present the virial equation of state of low-density nuclear matter composed of neutrons, protons and alpha particles. The virial equation of state is model-independent, and therefore sets a benchmark for all nuclear equations of state at low densities. We calculate the second virial coefficients for nucleon-nucleon, nucleon-alpha and alpha-alpha interactions directly from the relevant binding energies and scattering phase shifts. The virial approach systematically takes into account contributions from bound nuclei and the resonant continuum, and consequently provides a framework to include strong-interaction corrections to nuclear statistical equilibrium models. The virial coefficients are used to make model-independent predictions for a variety of properties of nuclear matter over a range of densities, temperatures and compositions. Our results provide constraints on the physics of the neutrinosphere in supernovae. The resulting alpha particle concentration differs from all equations of state currently used in supernova simulations. Finally, the virial equation of state greatly improves our conceptual understanding of low-density nuclear matter.