We present an efficient phase retrieval approach for imaging systems with high numerical aperture based on the vectorial model of the point spread function. The algorithm is in the class of alternating minimization methods and can be adjusted for applications with either known or unknown amplitude of the field in the pupil. The algorithm outperforms existing solutions for high-numerical-aperture phase retrieval: (1) the generalization of the method of Hanser et al., based on extension of the scalar diffraction theory by representing the out-of-focus diversity applied to the image by a spherical cap, and (2) the method of Braat et al., which assumes through the use of extended Nijboer–Zernike expansion the phase to be smooth. The former is limited in terms of accuracy due to model deviations, while the latter is of high computational complexity and excludes phase retrieval problems where the phase is discontinuous or sparse. Extensive numerical results demonstrate the efficiency, robustness, and practicability of the proposed algorithm in various practically relevant simulations.