Working with Hamiltonians from chiral effective field theory, we develop a novel framework for describing arbitrary deformed medium-mass nuclei by combining the in-medium similarity renormalization group with the generator coordinate method. The approach leverages the ability of the first method to capture dynamic correlations and the second to include collective correlations without violating symmetries. We use our scheme to compute the matrix element that governs the neutrinoless double beta decay of $^{48}$Ca to $^{48}$Ti, and find it to have the value $0.61$, near or below the predictions of most phenomenological methods. The result opens the door to ab initio calculations of the matrix elements for the decay of heavier nuclei such as $^{76}$Ge, $^{130}$Te, and $^{136}$Xe.