Phase diversity is a technique for obtaining estimates of both the object and the phase, by exploiting the simultaneous collection of two (or more) short-exposure optical images, one of which has been formed by further blurring the conventional image in some known fashion. This paper concerns a fast computational algorithm based upon a regularized variant of the Gauss-Newton optimization method for phase diversity-based estimation when a Gaussian likelihood fit-to-data criterion is applied. Simulation studies are provided to demonstrate that the method is remarkably robust and numerically efficient.