Accurate nuclear matrix elements (NMEs) for neutrinoless double beta decays of candidate nuclei are important for the design and interpretation of future experiments. Significant progress has been made in the modeling of these NMEs from first principles. The NME for 48Ca shows a good agreement among three different ab initio calculations starting from the same nuclear interaction constructed within the chiral EFT and the same decay operator. These studies open the door to ab initio calculations of the matrix elements for the decay of heavier nuclei such as 76Ge, 130Te, and 136Xe. The ultimate goal is the computation of NMEs in many-body calculations with controllable approximations, using nuclear interactions and weak transition operators derived consistently from chiral EFT. We are expecting more progress towards this goal in the near future.