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 that the neutrino is a Majorana particle, demonstrate that lepton number conservation is violated in nature, and give 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. We will show preliminary results for the $\pi^{-} \rightarrow \pi^{+} e^{-} e^{-}$ transition amplitude computed on a $16^{3} \times 32$ domain wall fermion lattice with a pion mass of 420 MeV, and discuss improved methods applicable to general lattice calculations of $0 \nu \beta \beta$ decay amplitudes.