The images recorded by optical telescopes are often degraded by aberrations that induce phase variations in the pupil plane. Several wavefront sensing techniques have been proposed to estimate aberrated phases. One of them is phase diversity, for which the joint least-square approach introduced by Gonsalves et al. is a reference method to estimate phase coefficients from the recorded images. In this paper, we rely on the asymptotic theory of Toeplitz matrices to show that Gonsalves' technique provides a consistent phase estimator as the size of the images grows. No comparable result is yielded by the classical joint maximum likelihood interpretation (e.g., as found in the work by Paxman et al.). Finally, our theoretical analysis is illustrated through simulated problems.