In this paper, we consider a nonlinear least squares model for the phase retrieval problem. Since the Hessian matrix may not be positive definite and the Gauss–Newton (GN) matrix is singular at any optimal solution, we propose a modified Levenberg–Marquardt (LM) method, where the Hessian is substituted by a summation of the GN matrix and a regularization term. Similar to the well-known Wirtinger flow (WF) algorithm under certain assumptions, we start from an initial point provably close to the set of the global optimal solutions. Global linear convergence and local quadratic convergence to the global solution set are established by estimating the smallest nonzero eigenvalues of the GN matrix, establishing local error bound properties and constructing a modified regularization condition. The computational cost becomes tractable if a preconditioned conjugate gradient method is applied to solve the LM equation inexactly. Specifically, the pre-conditioner is constructed from the expectation of the LM coefficient matrix by assuming the independence between the measurements and iteration point. The initialization can be dropped by slightly modifying the algorithm. Preliminary numerical experiments show that our algorithm is robust and it is often faster than the WF method on both random examples and natural image recovery.