The nuclear matrix elements that govern the rate of neutrinoless double beta decay must be accurately calculated if experiments are to reach their full potential. Theorists have been working on the problem for a long time but have recently stepped up their efforts as ton-scale experiments have begun to look feasible. Here we review past and recent work on the matrix elements in a wide variety of nuclear models and discuss work that will be done in the near future. Ab initio nuclear-structure theory, which is developing rapidly, holds out hope of more accurate matrix elements with quantifiable error bars.