Accurate and robust wavefront reconstruction methods for pyramid wavefront sensors are in high demand, as these sensors are planned to be part of many instruments currently under development for ground-based telescopes. The pyramid sensor relates the incoming wavefront and its measurements in a nonlinear way. Nevertheless, almost all existing reconstruction algorithms are based on a linearization of the model. The assumption of a linear pyramid sensor response is justifiable in closed-loop adaptive optics (AO) when the measured phase information is small, but, depending on the system, may not be feasible due to unpreventable errors such as non-common path aberrations. In order to solve the nonlinear inverse problem of wavefront reconstruction from pyramid sensor data, we introduce two new methods based on the nonlinear Landweber and Landweber-Kaczmarz iterations. Using these algorithms, we experience high-quality wavefront estimation, especially for the non-modulated sensor, by still keeping the numerical effort feasible for large-scale AO systems.