While neutrino oscillation experiments have demonstrated that neutrinos have small, nonzero masses, much remains unknown about their properties and decay modes. One potential decay mode --- neutrinoless double beta decay ($0 \nu \beta \beta$) --- is a particularly interesting target of experimental searches, since its observation would imply both the violation of lepton number conservation in nature as well as the existence of at least one Majorana neutrino, in addition to giving further constraints on the neutrino masses and mixing angles. Relating experimental constraints on $0 \nu \beta \beta$ decay rates to the neutrino masses, however, requires theoretical input in the form of non-perturbative nuclear matrix elements which remain difficult to calculate reliably. In this talk we will discuss progress towards first-principles calculations of relevant nuclear matrix elements using lattice QCD and effective field theory techniques, assuming neutrinoless double beta decay mediated by a light Majorana neutrino.