  Mean-field spin glasses are families of random energy functions
(Hamiltonians) on high-dimensional product spaces. In this paper we consider
the case of Ising mixed $p$-spin models, namely Hamiltonians $H_N:\Sigma_N\to
{\mathbb R}$ on the Hamming hypercube $\Sigma_N = \{\pm 1\}^N$, which are
defined by the property that $\{H_N({\boldsymbol \sigma})\}_{{\boldsymbol
\sigma}\in \Sigma_N}$ is a centered Gaussian process with covariance ${\mathbb
E}\{H_N({\boldsymbol \sigma}_1)H_N({\boldsymbol \sigma}_2)\}$ depending only on
the scalar product $\langle {\boldsymbol \sigma}_1,{\boldsymbol
\sigma}_2\rangle$. The asymptotic value of the optimum $\max_{{\boldsymbol
\sigma}\in \Sigma_N}H_N({\boldsymbol \sigma})$ was characterized in terms of a
variational principle known as the Parisi formula, first proved by Talagrand
and, in a more general setting, by Panchenko. The structure of superlevel sets
is extremely rich and has been studied by a number of authors. Here we ask
whether a near optimal configuration ${\boldsymbol \sigma}$ can be computed in
polynomial time. We develop a message passing algorithm whose complexity
per-iteration is of the same order as the complexity of evaluating the gradient
of $H_N$, and characterize the typical energy value it achieves. When the
$p$-spin model $H_N$ satisfies a certain no-overlap gap assumption, for any
$\varepsilon>0$, the algorithm outputs ${\boldsymbol \sigma}\in\Sigma_N$ such
that $H_N({\boldsymbol \sigma})\ge (1-\varepsilon)\max_{{\boldsymbol \sigma}'}
H_N({\boldsymbol \sigma}')$, with high probability. The number of iterations is
bounded in $N$ and depends uniquely on $\varepsilon$. More generally,
regardless of whether the no-overlap gap assumption holds, the energy achieved
is given by an extended variational principle, which generalizes the Parisi
formula.
