Summary: We study the problem of solving a quadratic system of equations, i.e., recovering a vector signal \(\mathbf x\in \mathbb{R}^n\) from its magnitude measurements \(y_i=|\langle \mathbf a_i, \mathbf x\rangle|\), \(i=1,\dots, m\). We develop a gradient descent algorithm (referred to as RWF for reshaped Wirtinger flow) by minimizing the quadratic loss of the magnitude measurements. Comparing with Wirtinger flow (WF) [E. J. Candès et al., IEEE Trans. Inf. Theory 61, No. 4, 1985–2007 (2015; Zbl 1359.94069)], the loss function of RWF is nonconvex and nonsmooth, but better resembles the least-squares loss when the phase information is also available. We show that for random Gaussian measurements, RWF enjoys linear convergence to the true signal as long as the number of measurements is \(\mathcal{O}(n)\). This improves the sample complexity of WF (\(\mathcal{O}(n\log n)\)), and achieves the same sample complexity as truncated Wirtinger flow (TWF) [Yuxin Chen and E. J. Candès, Solving random quadratic systems of equations is nearly as easy as solving linear systems. In: Advances in Neural Information Processing Systems (NIPS) Vol. 1, 739–747 (2015)], but without any sophisticated truncation in the gradient loop. Furthermore, RWF costs less computationally than WF, and runs faster numerically than both WF and TWF. We further develop an incremental (stochastic) version of RWF (IRWF) and connect it with the randomized Kaczmarz method for phase retrieval. We demonstrate that IRWF outperforms existing incremental as well as batch algorithms with experiments.