  This paper deals with finding an $n$-dimensional solution $x$ to a system of
quadratic equations of the form $y_i=|\langle{a}_i,x\rangle|^2$ for $1\le i \le
m$, which is also known as phase retrieval and is NP-hard in general. We put
forth a novel procedure for minimizing the amplitude-based least-squares
empirical loss, that starts with a weighted maximal correlation initialization
obtainable with a few power or Lanczos iterations, followed by successive
refinements based upon a sequence of iteratively reweighted (generalized)
gradient iterations. The two (both the initialization and gradient flow) stages
distinguish themselves from prior contributions by the inclusion of a fresh
(re)weighting regularization technique. The overall algorithm is conceptually
simple, numerically scalable, and easy-to-implement. For certain random
measurement models, the novel procedure is shown capable of finding the true
solution $x$ in time proportional to reading the data $\{(a_i;y_i)\}_{1\le i
\le m}$. This holds with high probability and without extra assumption on the
signal $x$ to be recovered, provided that the number $m$ of equations is some
constant $c>0$ times the number $n$ of unknowns in the signal vector, namely,
$m>cn$. Empirically, the upshots of this contribution are: i) (almost) $100\%$
perfect signal recovery in the high-dimensional (say e.g., $n\ge 2,000$) regime
given only an information-theoretic limit number of noiseless equations,
namely, $m=2n-1$ in the real-valued Gaussian case; and, ii) (nearly) optimal
statistical accuracy in the presence of additive noise of bounded support.
Finally, substantial numerical tests using both synthetic data and real images
corroborate markedly improved signal recovery performance and computational
efficiency of our novel procedure relative to state-of-the-art approaches.
