The termination in the white dwarf luminosity function is a standard diagnostic tool for measuring the total age of nearby stellar populations. In this paper, an algorithm is presented for inverting the full white dwarf luminosity function to obtain a maximum likelihood estimate of the time-varying star formation rate of the host stellar population. Tests with synthetic data demonstrate that the algorithm converges over a wide class of underlying star formation rate forms. The algorithm successfully estimates the moving average star formation rate as a function of lookback time in the presence of realistic measurement noise, though suffers from degeneracies around discontinuities in the underlying star formation rate. The inversion results are most sensitive to the choice of white dwarf cooling models, with the models produced by different groups giving quite different results. The results are relatively insensitive to the progenitor metallicity, initial mass function, initial-final mass relation and ratio of H/He atmosphere white dwarfs. Application to two independent determinations of the solar neighbourhood white dwarf luminosity function gives similar results. The star formation rate has a bimodal form, with broad peaks at 2-3 and 7-9 Gyr in the past, separated by a significant lull of magnitude 30-90 per cent depending on the choice of cooling models. The onset of star formation occurs around 8-10 Gyr ago. The total integrated star formation rate is ~0.014 stars pc-3 in the solar neighbourhood, for stars more massive than 0.6M.