  Phase retrieval refers to the problem of recovering a signal
$\mathbf{x}_{\star}\in\mathbb{C}^n$ from its phaseless measurements
$y_i=|\mathbf{a}_i^{\mathrm{H}}\mathbf{x}_{\star}|$, where
$\{\mathbf{a}_i\}_{i=1}^m$ are the measurement vectors. Many popular phase
retrieval algorithms are based on the following two-step procedure: (i)
initialize the algorithm based on a spectral method, (ii) refine the initial
estimate by a local search algorithm (e.g., gradient descent). The quality of
the spectral initialization step can have a major impact on the performance of
the overall algorithm. In this paper, we focus on the model where the
measurement matrix $\mathbf{A}=[\mathbf{a}_1,\ldots,\mathbf{a}_m]^{\mathrm{H}}$
has orthonormal columns, and study the spectral initialization under the
asymptotic setting $m,n\to\infty$ with $m/n\to\delta\in(1,\infty)$. We use the
expectation propagation framework to characterize the performance of spectral
initialization for Haar distributed matrices. Our numerical results confirm
that the predictions of the EP method are accurate for not-only Haar
distributed matrices, but also for realistic Fourier based models (e.g. the
coded diffraction model). The main findings of this paper are the following:
  (1) There exists a threshold on $\delta$ (denoted as
$\delta_{\mathrm{weak}}$) below which the spectral method cannot produce a
meaningful estimate. We show that $\delta_{\mathrm{weak}}=2$ for the
column-orthonormal model. In contrast, previous results by Mondelli and
Montanari show that $\delta_{\mathrm{weak}}=1$ for the i.i.d. Gaussian model.
  (2) The optimal design for the spectral method coincides with that for the
i.i.d. Gaussian model, where the latter was recently introduced by Luo,
Alghamdi and Lu.
