  We consider the phase retrieval problem, in which the observer wishes to
recover a $n$-dimensional real or complex signal $\mathbf{X}^\star$ from the
(possibly noisy) observation of $|\mathbf{\Phi} \mathbf{X}^\star|$, in which
$\mathbf{\Phi}$ is a matrix of size $m \times n$. We consider a
\emph{high-dimensional} setting where $n,m \to \infty$ with $m/n =
\mathcal{O}(1)$, and a large class of (possibly correlated) random matrices
$\mathbf{\Phi}$ and observation channels. Spectral methods are a powerful tool
to obtain approximate observations of the signal $\mathbf{X}^\star$ which can
be then used as initialization for a subsequent algorithm, at a low
computational cost. In this paper, we extend and unify previous results and
approaches on spectral methods for the phase retrieval problem. More precisely,
we combine the linearization of message-passing algorithms and the analysis of
the \emph{Bethe Hessian}, a classical tool of statistical physics. Using this
toolbox, we show how to derive optimal spectral methods for arbitrary channel
noise and right-unitarily invariant matrix $\mathbf{\Phi}$, in an automated
manner (i.e. with no optimization over any hyperparameter or preprocessing
function).
