\documentclass[aps,pre,reprint,longbibliography,floatfix]{revtex4-2}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb,latexsym,graphicx}
\usepackage{newtxtext,newtxmath}
\usepackage{bbold,anyfontsize}
\usepackage[dvipsnames]{xcolor}

\begin{document}

\title{
  Conditioning the complexity of random landscapes on marginal optima
}

\author{Jaron Kent-Dobias}
\affiliation{Istituto Nazionale di Fisica Nucleare, Sezione di Roma I, Rome, Italy 00184}

\begin{abstract}
  Marginal optima are minima or maxima of a function with many asymptotically
  flat directions. In settings with many competing optima, marginal ones tend
  to attract algorithms and physical dynamics. Often, the important family of
  marginal attractors are a vanishing minority compared with nonmarginal optima
  and other unstable stationary points. We introduce a generic technique for
  conditioning the statistics of stationary points on their marginality, and
  apply it in three isotropic settings with different typical forms for the
  Hessian at optima: in the spherical spin-glasses, where the Hessian is GOE;
  in a multispherical spin glasses, which are Gaussian but non-GOE; and in a
  model of random nonlinear sum of squares, which is non-Gaussian. In these
  problems we are able to fully characterize the distribution of marginal
  optima in the landscape, including when they are in the minority.
\end{abstract}

\maketitle

\section{Introduction}

Systems with rugged landscapes are important across many disciplines, from the
physics to glasses and spin-glasses to the statistical inference problems. The
behavior of these systems is best understood when equilibrium or optimal
solutions are studied and averages can be taken statically over all possible
configurations. However, such systems are also infamous for their tendency to
defy equilibrium and optimal expectations in practice, due to the presence of
dynamic transitions or crossovers that leave physical or algorithmic dynamics
stuck exploring only a subset of configurations.

In some simple models of such landscapes, it was recently found that marginal
minima are significant as the attractors of gradient descent dynamics
\cite{Folena_2020_Rethinking, Folena_2023_On}. This extends to more novel
algorithms, like message passing \cite{} \textbf{Find out if this is true}.
\textbf{Think of other examples.}
While it is still not known how to predict which marginal minima will be
attractors, this ubiquity of behavior suggests that cartography of marginal
minima is a useful step in bounding out-of-equilibrium dynamical behavior.

In the traditional methods for analyzing the geometric structure of rugged
landscapes, it is not necessarily straightforward to condition an analysis on
the marginality of minima. Using the method of a Legendre transformation of the
Parisi parameter corresponding to a set of real replicas, one can force the
result to be marginal by restricting the value of that parameter, but this
results in only the marginal minima at the energy level at which they are the
majority of stationary points \cite{Monasson_1995_Structural}. It is now
understood that out-of-equilibrium dynamics usually goes to marginal minima at
other energy levels \cite{Folena_2023_On}.

The alternative, used to great success in the spherical models, is to start by
making a detailing understanding of the Hessian matrix at stationary points.
Then, one can condition the analysis on whatever properties of the Hessian are
necessary to lead to marginal minima. This strategy is so successful in the
spherical models because it is very straightforward to implement: a natural
parameter in the analysis of these models linearly shifts the spectrum of the
Hessian, and so fixing this parameter by whatever means naturally allows one to
require that the Hessian spectrum have a pseudogap.
Unfortunately this strategy is less straightforward to generalize. Many models
of interest, especially in inference problems, have Hessian statistics that are
poorly understood.

Here, we introduce a generic method for conditioning the statistics of
stationary points on their marginality. The technique makes use of a novel way
to condition an integral over parameters to select only those that result in a
certain value of the smallest eigenvalue of a matrix that is a function of
those parameters. By requiring that the smallest eigenvalue of the Hessian at
stationary points be zero, we restrict to marginal minima, either those with a
pseudogap in their bulk spectrum or those with outlying eigenvectors. We
provide a heuristic to distinguish these two cases. We demonstrate the method
on the spherical models, where it is unnecessary but instructive, and on
extensions of the spherical models with non-GOE Hessians where the technique is
more useful.

\section{Conditioning on the smallest eigenvalue}



An arbitrary function $g$ of the minimum eigenvalue of a matrix $A$ can be
expressed as
\begin{equation} \label{eq:λmax}
  g(\lambda_\textrm{max}(A))
  =\lim_{\beta\to\infty}\int
    \frac{d\mathbf s\,\delta(N-\mathbf s^T\mathbf s)e^{\beta\mathbf s^TA\mathbf s}}
      {\int d\mathbf s'\,\delta(N-\mathbf s'^T\mathbf s')e^{\beta\mathbf s'^TA\mathbf s'}}
    g\left(\frac{\mathbf s^TA\mathbf s}N\right)
\end{equation}
Assuming
\begin{equation}
  \begin{aligned}
    &\lim_{\beta\to\infty}\int\frac{d\mathbf s\,\delta(N-\mathbf s^T\mathbf s)e^{\beta\mathbf s^TA\mathbf s}}{\int d\mathbf s'\,\delta(N-\mathbf s'^T\mathbf s')e^{\beta\mathbf s'^TA\mathbf s'}}g\left(\frac{\mathbf s^TA\mathbf s}N\right) \\
    &=\int\frac{d\mathbf s\,\delta(N-\mathbf s^T\mathbf s)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{max}(A)I)}(\mathbf s)}{\int d\mathbf s'\,\delta(N-\mathbf s'^T\mathbf s')\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{max}(A)I)}(\mathbf s')}g\left(\frac{\mathbf s^TA\mathbf s}N\right) \\
    &=g(\lambda_\mathrm{max}(A))
    \frac{\int d\mathbf s\,\delta(N-\mathbf s^T\mathbf s)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{max}(A)I)}(\mathbf s)}{\int d\mathbf s'\,\delta(N-\mathbf s'^T\mathbf s')\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{max}(A)I)}(\mathbf s')} \\
    &=g(\lambda_\mathrm{max}(A))
  \end{aligned}
\end{equation}
The first relation extends a technique first introduced in
\cite{Ikeda_2023_Bose-Einstein-like} and used in
\cite{Kent-Dobias_2024_Arrangement}. A Boltzmann distribution is introduced
over a spherical model whose Hamiltonian is quadratic with interaction matrix
given by $A$. In the limit of zero temperature, the measure will concentrate on
the ground states of the model, which correspond with the eigenspace of $A$
associated with its minimum eigenvalue $\lambda_\mathrm{min}$. The second
relation uses the fact that, once restricted to the sphere $\mathbf s^T\mathbf
s=N$ and the minimum eigenspace, $\mathbf s^TA\mathbf s=N\lambda_\mathrm{min}(A)$.

The relationship is formal, but we can make use of the fact that the integral
expression with a Gibbs distribution can be manipulated with replica
techniques, averaged over, and in general treated with a physicist's toolkit.
In particular, we have specific interest in using
$g(\lambda_\mathrm{min}(A))=\delta(\lambda_\mathrm{min}(A))$, a Dirac
delta-function, which can be inserted into averages over ensembles of matrices
$A$ (or indeed more complicated averages) in order to condition that the
minimum eigenvalue is zero.

\subsection{Simple example: shifted GOE}

We demonstrate the efficacy of the technique by rederiving a well-known result:
the large-deviation function for pulling an eigenvalue from the bulk of the
GOE spectrum.
Consider an ensemble of $N\times N$ matrices $A=B+\mu I$ for $B$ drawn from the GOE ensemble with entries
whose variance is $\sigma^2/N$. We know that the bulk spectrum of $A$ is a
Wigner semicircle with radius $2\sigma$ shifted by a constant $\mu$.
Therefore, for $\mu=2\sigma$, the minimum eigenvalue will typically be zero,
while for $\mu>2\sigma$ the minimum eigenvalue would need to be a large
deviation from the typical spectrum and its likelihood will be exponentially
suppressed with $N$. For $\mu<2\sigma$, the bulk of the typical spectrum contains
zero and therefore a larger $N^2$ deviation, moving an extensive number of
eigenvalues, would be necessary. This final case cannot be quantified by this
method, but instead the nonexistence of a large deviation linear in $N$ appears
as the emergence of an imaginary part in the function.

As an example, we compute
\begin{equation} \label{eq:large.dev}
  e^{NG_\lambda^*(\mu)}
  =P_{\lambda_\mathrm{max}(B-\mu I)=\lambda^*}
  =\overline{\delta\big(N\lambda^*-N\lambda_\mathrm{max}(B-\mu I)\big)}
\end{equation}
where the overline is the average over $B$, and we have defined the large
deviation function $G_\sigma(\mu)$.
Using the representation of $\lambda_\mathrm{max}$ defined in \eqref{eq:λmax}, we have
\begin{widetext}
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\overline{
    \lim_{\beta\to\infty}\int\frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)e^{\beta\mathbf x^T(B-\mu I)\mathbf x}}
    {\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')e^{\beta\mathbf x'^T(B-\mu I)\mathbf x'}}\,\delta\big(N\lambda^*-\mathbf x^T(B-\mu I)\mathbf x\big)
  }
\end{equation}
Using replicas to treat the denominator ($x^{-1}=\lim_{n\to0}x^{n-1}$)
and transforming the $\delta$-function to its Fourier
representation, we have
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\overline{\lim_{\beta\to\infty}\lim_{n\to0}\int d\hat\lambda\prod_{a=1}^n\left[d\mathbf x_a\,\delta(N-\mathbf x_a^T\mathbf x_a)\right]
  \exp\left\{\beta\sum_{a=1}^n\mathbf x_a^T(B-\mu I)\mathbf x_a+\hat\lambda\left[N\lambda^*-\mathbf x_1^T(B-\mu I)\mathbf x_1\right]\right\}}
\end{equation}
having introduced the parameter $\hat\lambda$ in the Fourier representation of the $\delta$-function.
The whole expression, so transformed, is a simple exponential integral linear in the matrix $B$.
Taking the average over $B$, we have
\begin{equation}
  \begin{aligned}
  &e^{NG_{\lambda^*}(\mu)}
  =\lim_{\beta\to\infty}\lim_{n\to0}\int d\hat\lambda\prod_{a=1}^n\left[d\mathbf x_a\,\delta(N-\mathbf x_a^T\mathbf x_a)\right] \\
  &\hspace{10em}\exp\left\{N\left[\hat\lambda(\mu+\lambda^*)-n\beta\mu\right]+\frac{\sigma^2}{N}\left[\beta^2\sum_{ab}^n(\mathbf x_a^T\mathbf x_b)^2
    -2\beta\hat\lambda\sum_a^n(\mathbf x_a^T\mathbf x_1)^2
    +\hat\lambda^2N^2
  \right]\right\}
  \end{aligned}
\end{equation}
We make the Hubbard--Stratonovich transformation to the matrix field $Q_{ab}=\frac1N\mathbf x_a^T\mathbf x_b$. This gives
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\lim_{\beta\to\infty}\lim_{n\to0}\int d\hat\lambda\,dQ\,
  \exp N\left\{
    \hat\lambda(\mu+\lambda^*)-n\beta\mu+\sigma^2\left[\beta^2\sum_{ab}^nQ_{ab}^2
      -2\beta\hat\lambda\sum_a^nQ_{1a}^2
    +\hat\lambda^2
  \right]+\frac12\log\det Q\right\}
\end{equation}
\end{widetext}
where $Q_{aa}=1$ because of the spherical constraint. We can evaluate this
integral using the saddle point method. We make a replica symmetric ansatz for
$Q$, because this is a 2-spin model, but with the first row singled out because
of its unique coupling with $\hat\lambda$. This gives
\begin{equation}
  Q=\begin{bmatrix}
    1&\tilde q_0&\tilde q_0&\cdots&\tilde q_0\\
    \tilde q_0&1&q_0&\cdots&q_0\\
    \tilde q_0&q_0&1&\ddots&q_0\\
    \vdots&\vdots&\ddots&\ddots&\vdots\\
    \tilde q_0&q_0&q_0&\cdots&q_0
  \end{bmatrix}
\end{equation}
with $\sum_{ab}Q_{ab}^2=n+2(n-1)\tilde q_0^2+(n-1)(n-2)q_0^2$, $\sum_aQ_{1a}^2=1+(n-1)\tilde q_0^2$,
and
\begin{equation}
  \log\det Q=(n-2)\log(1-q_0)+\log(1+(n-2)q_0-(n-1)\tilde q_0^2)
\end{equation}
Inserting these expressions and taking the limit of $n$ to zero, we find
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}=\lim_{\beta\to\infty}\int d\hat\lambda\,dq_0\,d\tilde q_0\,e^{N\mathcal S_\beta(q_0,\tilde q_0,\hat\lambda)}
\end{equation}
with the effective action
\begin{equation}
  \begin{aligned}
    &\mathcal S_\beta(q_0,\tilde q_0,\hat\lambda) \\
    &\quad=\hat\lambda(\mu+\lambda^*)+\sigma^2\left[
      2\beta^2(q_0^2-\tilde q_0^2)-2\beta\hat\lambda(1-\tilde q_0^2)+\hat\lambda^2
    \right] \\
    &\qquad-\log(1-q_0)+\frac12\log(1-2q_0+\tilde q_0^2)
  \end{aligned}
\end{equation}
We need to evaluate the integral above using the saddle point method, but in the limit of $\beta\to\infty$.
We expect the overlaps to concentrate on one as $\beta$ goes to infinity. We therefore take
\begin{align}
  q_0&=1-y\beta^{-1}-z\beta^{-2}+O(\beta^{-3})
  \\
  \tilde q_0&=1-\tilde y\beta^{-1}-(z-\Delta z)\beta^{-2}+O(\beta^{-3})
\end{align}
However, taking the limit with $y\neq\tilde y$ results in an expression for the
action that diverges with $\beta$. To cure this, we must take $\tilde y=y$. The result is
\begin{equation}
  \begin{aligned}
    \mathcal S_\infty(y,\Delta z,\hat\lambda)
    &=\hat\lambda(\mu+\lambda^*)
    +\sigma^2\big[
      \hat\lambda^2-4(y+\Delta z)
    \big] \\
    &\qquad+\frac12\log\left(1+\frac{2\Delta z}{y^2}\right)
  \end{aligned}
\end{equation}
Extremizing this action over the new parameters $y$, $\Delta z=z-\tilde z$, and $\hat\lambda$, we have
\begin{align}
  \hat\lambda=-\frac1\sigma\sqrt{\frac{(\mu+\lambda^*)^2}{(2\sigma)^2}-1}
  \\
  y=\frac1{2\sigma}\left(\frac{\mu+\lambda^*}{2\sigma}-\sqrt{\frac{(\mu+\lambda^*)^2}{(2\sigma)^2}-1}\right)
  &\\
  \Delta z=\frac1{4\sigma^2}\left(1-\frac{\mu+\lambda^*}{2\sigma}\left(\frac{\mu+\lambda^*}{2\sigma}-\sqrt{\frac{(\mu+\lambda^*)^2}{(2\sigma)^2}-1}\right)\right)
\end{align}
Inserting this solution into $\mathcal S_\infty$ we find
\begin{equation}
  \begin{aligned}
    &G_{\lambda^*}(\mu)
    =\mathop{\textrm{extremum}}_{y,\Delta z,\hat\lambda}\mathcal S_\infty(y,\Delta z,\hat\lambda) \\
    &=-\tfrac{\mu+\lambda^*}{2\sigma}\sqrt{\Big(\tfrac{\mu+\lambda^*}{2\sigma}\Big)^2-1}
    +\log\left(
      \tfrac{\mu+\lambda^*}{2\sigma}+\sqrt{\Big(\tfrac{\mu+\lambda^*}{2\sigma}\Big)^2-1}
    \right)
  \end{aligned}
\end{equation}
This function is plotted in Fig.~\ref{fig:large.dev}. For $\mu<2\sigma$ $G_{\lambda^*}(\mu)$ has an
imaginary part, which makes any additional integral over $\mu$ highly
oscillatory. This indicates that the existence of a marginal minimum for this
parameter value corresponds with a large deviation that grows faster than $N$,
rather like $N^2$, since in this regime the bulk of the typical spectrum is
over zero and therefore extensively many eigenvalues have to have large
deviations in order for the smallest eigenvalue to be zero. For
$\mu\geq2\sigma$ this function gives the large deviation function for the
probability of seeing a zero eigenvalue given the shift $\mu$.
$\mu=2\sigma$ is the maximum of the function with a real value, and
corresponds to the intersection of the average spectrum with zero, i.e., a
pseudogap.

\begin{figure}
  \includegraphics[width=\columnwidth]{figs/large_deviation.pdf}
  \caption{
    The large deviation function $G_\sigma(\mu)$ defined in
    \eqref{eq:large.dev} as a function of the shift $\mu$ to the
    GOE diagonal. As expected, $G_\sigma(2\sigma)=0$, while for
    $\mu>2\sigma$ it is negative and for $\mu<2\sigma$ it gains an
    imaginary part.
  } \label{fig:large.dev}
\end{figure}

Marginal spectra with a pseudogap and those with simple isolated eigenvalues
are qualitatively different, and more attention may be focused on the former.
Here, we see what appears to be a general heuristic for identifying the saddle
parameters for which the spectrum is pseudogapped: the equivalent of this
large-deviation functions will lie on the singular boundary between a purely
real and complex value.

\subsection{Conditioning on a pseudogap}

We have seen that this method effectively conditions a random matrix ensemble
on its lowest eigenvalue being zero. However, this does not correspond on its
own to marginal minima. In the previous example, most values of $\mu$ where
the calculation was valid correspond to matrices with a single isolated
eigenvalue. However, the marginal minima we are concerned with have
pseudogapped spectra, where the continuous part of the spectral density has a
lower bound at zero.

Fortunately, our calculation can be modified to ensure that we consider only
pseudogapped spectra. First, we insert a shift $\mu$ by hand into the `natural'
spectrum of the problem at hand, conditioning the trace to have a specific
value $\mu=\operatorname{Tr}A$. Then, we choose this artificial shift so that
the resulting conditioned spectra are pseudogapped. As seen the previous
subsection, this can be done by starting from a sufficiently large $\mu$ and
decreasing it until the calculation develops an imaginary part, signaling the
breakdown of the large-deviation principle at order $N$.

In isotropic or zero-signal landscapes, there is another way to condition on a
pseudogap. In such landscapes, the typical spectrum does not have an isolated
eigenvalue. Therefore, the condition associated with the bulk of the spectrum
teaching zero, i.e., the pseudogap, will always correspond to the most common
configuration. We can therefore choose $\mu=\mu_\textrm m$ such that
\begin{equation}
  0=\frac\partial{\partial\lambda^*}G_{\lambda^*}(\mu_\mathrm m)\bigg|_{\lambda^*=0}
\end{equation}
In the previous problem, this corresponds precisely to $\mu_\mathrm m=2\sigma$,
the correct marginal shift. Note that when we treat the Dirac $\delta$ function
using its Fourier representation with auxiliary parameter $\hat\lambda$, as in
the previous subsection, this condition corresponds with choosing $\mu$ such
that $\hat\lambda=0$.


\section{Marginal complexity in random landscapes}

The situation in the study of random landscapes is often as follows: an
ensemble of smooth functions $H:\mathbb R^N\to\mathbb R$ define random
landscapes, often with their configuration space subject to one or more
constraints of the form $g(\mathbf s)=0$ for $\mathbf s\in\mathbb R^N$. The
geometry of such landscapes is studied by their complexity, or the average
logarithm of the number of stationary points with certain properties, e.g., of
marginal minima at a given energy.

Such problems can be studied using the method of Lagrange multipliers, with one introduced for every constraint. If the configuration space is defined by $r$ constraints, then the problem is to extremize
\begin{equation}
  H(\mathbf s)+\sum_{i=1}^r\omega_ig_i(\mathbf s)
\end{equation}
with respect to $\mathbf s$ and $\pmb\omega=\{\omega_1,\ldots,\omega_r\}$. The corresponding gradient and Hessian for the problem are
\begin{align}
  \nabla H(\mathbf s,\pmb\omega)=\partial H(\mathbf s)+\sum_{i=1}^r\omega_i\partial g_i(\mathbf s)
  \\
  \operatorname{Hess}H(\mathbf s,\pmb\omega)=\partial\partial H(\mathbf s)+\sum_{i=1}^r\omega_i\partial\partial g_i(\mathbf s)
\end{align}
The number of stationary points in a landscape for a particular realization $H$ is found by integrating over the Kac--Rice measure
\begin{equation}
  d\mu_H(\mathbf s,\pmb\omega)=d\mathbf s\,d\pmb\omega\,\delta\big(\nabla H(\mathbf s,\pmb\omega)\big)\,\delta\big(\mathbf g(\mathbf s)\big)\,\big|\det\operatorname{Hess}H(\mathbf s,\pmb\omega)\big|
\end{equation}
with a $\delta$-function of the gradient and the constraints ensuring that we
count valid stationary points, and the Hessian entering in the determinant of
the Jacobian of the argument to the $\delta$-function. It is usually more
interesting to condition the count on interesting properties of the stationary
points, like the energy,
\begin{equation}
  \begin{aligned}
    &d\mu_H(\mathbf s,\pmb\omega\mid E,\mu) \\
    &\quad=d\mu_H(\mathbf s,\pmb\omega)\,
    \delta\big(NE-H(\mathbf s)\big)
    \,\delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf x,\pmb\omega)\big)
  \end{aligned}
\end{equation}
In this paper we in particular want to exploit our method to condition
complexity on the marginality of stationary points. We therefore define the
number of marginal points in a particular instantiation $H$ as
\begin{widetext}
\begin{equation}
  \begin{aligned}
    &\mathcal N_H(E,\mu,\lambda^*)
    =\int d\mu_H(\mathbf s,\pmb\omega\mid E,\mu)\,\delta\big(N\lambda^*-\lambda_\mathrm{max}(\operatorname{Hess}H(\mathbf s,\pmb\omega))\big) \\
    &=\lim_{\beta\to\infty}\int d\mu_H(\mathbf s,\pmb\omega\mid E,\mu)
    \frac{d\mathbf x\,\delta(N-\mathbf x^T\mathbf x)\delta(\mathbf x^T\partial\mathbf g(\mathbf s))e^{\beta\mathbf x^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x}}
    {\int d\mathbf x'\,\delta(N-\mathbf x'^T\mathbf x')\delta(\mathbf x'^T\partial\mathbf g(\mathbf s))e^{\beta\mathbf x'^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x'}}
    \delta\big(N\lambda^*-\mathbf x^T\operatorname{Hess}H(\mathbf s,\pmb\omega)\mathbf x\big)
  \end{aligned}
\end{equation}
where the $\delta$-functions
\begin{equation}
  \delta(\mathbf x^T\partial\mathbf g(\mathbf s))
  =\prod_{s=1}^r\delta(\mathbf x^T\partial g_i(\mathbf s))
\end{equation}
ensure that the integrals are constrained to the tangent space of the
configuration manifold at the point $\mathbf s$. This likewise allows us to
define the complexity of points with a specific energy, stability, and maximum eigenvalue as
\begin{equation}
  \Sigma_{\lambda^*}(E,\mu)
  =\frac1N\overline{\log\mathcal N_H(E,\mu,\lambda^*)}
\end{equation}
In practice, this can be computed by introducing replicas to treat the
logarithm ($\log x=\lim_{n\to0}\frac\partial{\partial n}x^n$) and replicating
again to treat each of the normalizations in the numerator. This leads to the expression
\begin{equation}
  \begin{aligned}
    \Sigma_{\lambda^*}(E,\mu)
    &=\lim_{\beta\to\infty}\lim_{n\to0}\frac1N\frac\partial{\partial n}\int\prod_{a=1}^n\Bigg[d\mu_H(\mathbf s_a,\pmb\omega_a\mid E,\mu)\,\delta\big(N\lambda^*-(\mathbf x_a^1)^T\operatorname{Hess}H(\mathbf s_a,\pmb\omega_a)\mathbf x_a^1\big)\\
    &\hspace{12em}\times\lim_{m_a\to0}
      \left(\prod_{b=1}^{m_a} d\mathbf x_a^b\,\delta(N-(\mathbf x_a^b)^T\mathbf x_a^b)\delta((\mathbf x_a^b)^T\partial\mathbf g(\mathbf s_a))e^{\beta(\mathbf x_a^b)^T\operatorname{Hess}H(\mathbf s_a,\pmb\omega_a)\mathbf x_a^b}\right)\Bigg]
  \end{aligned}
\end{equation}
\end{widetext}
Finally, the \emph{marginal} complexity is given by fixing $\mu=\mu_\text{m}$ so that the complexity is stationary with respect to changes in the value of the maximum eigenvalue, or
\begin{equation}
  0=\frac\partial{\partial\lambda^*}\Sigma_{\lambda^*}(E,\mu_\text{m}(E))\bigg|_{\lambda^*=0}
\end{equation}
Finally, the marginal complexity is defined by evaluating the complexity conditioned on $\lambda_{\text{min}}=0$ at $\mu_\text{m}$,
\begin{equation}
  \Sigma_\text{m}(E)
  =\Sigma_0(E,\mu_\text m(E))
\end{equation}

\section{Examples}

\subsection{Spherical spin glasses}

\begin{align}
  C_{ab}=\frac1N\mathbf s_a\cdot\mathbf s_b
  &&
  R_{ab}=-i\frac1N\mathbf s_a\cdot\hat{\mathbf s}_b
  &&
  D_{ab}=\frac1N\hat{\mathbf s}_a\cdot\hat{\mathbf s}_b
  \\
  A_{ab}^{cd}=\frac1N\mathbf x_a^c\cdot\mathbf x_b^d
  &&
  X^c_{ab}=\frac1N\mathbf s_a\cdot\mathbf x_b^c
  &&
  \hat X^c_{ab}=\frac1N\hat{\mathbf s}_a\cdot\mathbf x_b^c
\end{align}

\begin{widetext}
\begin{equation}
  \begin{aligned}
    &\sum_{ab}^n\left[\beta\omega A_{aa}^{bb}+\hat x\omega A_{aa}^{11}+\beta^2f''(1)\sum_{cd}^m(A_{ab}^{cd})^2+\hat x^2f''(1)(A_{ab}^{11})^2+\beta\hat xf''(1)\sum_c^m A_{ab}^{1c}\right]\\
    &+\hat\beta^2f(C_{ab})+(2\hat\beta(R_{ab}-F_{ab})-D_{ab})f'(C_{ab})+(R_{ab}^2-F_{ab}^2)f''(C_{ab})
    +\log\det\begin{bmatrix}C&iR\\iR&D\end{bmatrix}-\log\det F
  \end{aligned}
\end{equation}
\end{widetext}

$X^a$ is $n\times m_a$, and $A^{ab}$ is $m_a\times m_b$.
\begin{equation}
  \begin{bmatrix}
    C&iR&X^1&\cdots&X^n \\
    iR&D&i\hat X^1&\cdots&i\hat X^m\\
    (X^1)^T&i(\hat X^1)^T&A^{11}&\cdots&A^{1n}\\
    \vdots&\vdots&\vdots&\ddots&\vdots\\
    (X^n)^T&i(\hat X^n)^T&A^{n1}&\cdots&A^{nn}
  \end{bmatrix}
\end{equation}
$X_{ab}^c$ will be nonzero if the lowest eigenvector of the hessian at the
point $\mathbf s_c$ are correlated with the direction of the point $\mathbf
s_a$. Since the eigenvector problem is always expected to be replica symmetric,
we expect no $b$-dependence of this matrix. $A^{aa}$ is the usual
replica-symmetric overlap matrix of the spherical two-spin problem. $A^{ab}$
describes overlaps between eigenvectors at different stationary points and should be a constant $m_a\times m_b$ matrix.

We will discuss at the end of this paper when these order parameters can be expected to be nonzero, but in this and most isotropic problems all of the $X$s, $\hat X$s, and $A^{ab}$ for $a\neq b$ are zero.

\begin{equation}
  \Sigma_\textrm{marginal}(E)
  =\operatorname{max}_\omega\big[\Sigma(E,\omega)+G_{\sqrt{f''(1)}}(\omega)\big]
\end{equation}
where the maximum over $\omega$ needs to lie at a real value.

\subsection{Twin spherical spin glasses}

$\Omega=S^{N-1}\times S^{N-1}$
\begin{equation}
  H(\pmb\sigma)=H_1(\pmb\sigma^{(1)})+H_2(\pmb\sigma^{(2)})+\epsilon\pmb\sigma^{(1)}\cdot\pmb\sigma^{(2)}
\end{equation}
\begin{equation}
  \overline{H_s(\pmb\sigma_1)H_s(\pmb\sigma_2)}
  =Nf_s\left(\frac{\pmb\sigma_1\cdot\pmb\sigma_2}N\right)
\end{equation}

\begin{widetext}
\begin{equation}
  \mathcal S(C,R,D,W,\hat\beta,\omega)
  =\frac12\frac1n
    \sum_{ab}\left(
      \hat\beta^2f(C_{ab})+(2\hat\beta R_{ab}-D_{ab})f'(C_{ab})+(R_{ab}^2-W_{ab}^2)f''(C_{ab})
    \right)
\end{equation}

\begin{equation}
  \begin{aligned}
    &\mathcal S(C^{11},R^{11},D^{11},W^{11},\hat\beta)+\mathcal S(C^{22},R^{22},D^{22},W^{22},\hat\beta)
    -\epsilon(r_{12}+r_{21})-\omega_1(r^{11}_d-w^{11}_d)-\omega_2(r^{22}_d-w^{22}_d)+\hat\beta E \\
    &+\frac12\log\det\begin{bmatrix}C^{11}&iR^{11}\\iR^{11}&D^{11}\end{bmatrix}
    +\frac12\log\det\left(
    \begin{bmatrix}C^{22}-q_{12}^2C^{11}&iR^{22}\\iR^{22}&D^{22}\end{bmatrix}
    \right)
    -\log\det(W^{11}W^{22}+W^{12}W^{21})
  \end{aligned}
\end{equation}

\begin{equation}
  \begin{aligned}
    &\sum_a^n\left[\hat q_a(Q^{11}_{aa}+Q^{22}_{aa}-1)-\beta(\omega_1Q^{11}_{aa}+\omega_2Q^{22}_{aa}+2\epsilon Q^{12}_{aa})\right]
    +\lambda(\omega_1Q^{11}_{11}+\omega_2Q^{22}_{11}+2\epsilon Q^{12}_{11}) \\
    &+\sum_{i=1,2}f_i''(1)\left[\beta^2\sum_{ab}^n(Q^{ii}_{ab})^2-2\beta\lambda\sum_a^n(Q^{ii}_{1a})^2+\lambda^2(Q^{ii}_{11})^2\right]
    +\frac12\log\det\begin{bmatrix}
      Q^{11}&Q^{12}\\
      Q^{12}&Q^{22}
    \end{bmatrix}
  \end{aligned}
\end{equation}
\end{widetext}
\begin{equation}
  \log\det\begin{bmatrix}
    Q^{11}&Q^{12}\\
    Q^{12}&Q^{22}
  \end{bmatrix}
  +\log\det(Q^{11}Q^{22}-Q^{12}Q^{12})
\end{equation}

\subsection{Random nonlinear least squares}

In this subsection we consider perhaps the simplest example of a non-Gaussian
landscape: the problem of random nonlinear least squares optimization. Though,
for reasons we will see it is easier to make predictions for random nonlinear
\emph{most} squares, i.e., the problem of maximizing the sum of squared terms.
We also take a spherical problem with $\mathbf x\in S^{N-1}$, and consider a set
of $M$ random functions $V_k:\mathbf S^{N-1}\to\mathbb R$ that are centered Gaussians with covariance
\begin{equation}
  \overline{V_i(\mathbf x)V_j(\mathbf x')}=\delta_{ij}f\left(\frac{\mathbf x^T\mathbf x'}N\right)
\end{equation}
The energy or cost function is the sum of squares of the $V_k$, or
\begin{equation}
  H(\mathbf x)=\frac12\sum_{k=1}^MV_k(\mathbf x)^2
\end{equation}
The landscape complexity and large deviations of the ground state for this problem were recently studied in a linear context, with $f(q)=\sigma^2+aq$ \cite{Fyodorov_2020_Counting, Fyodorov_2022_Optimization}. Some results on the ground state of the general nonlinear problem can also be found in \cite{Tublin_2022_A}. In particular, that work indicates that the low-lying minima of the problem tend to be either replica symmetric or full replica symmetry breaking. This is not good news for our analysis or marginal states, because in the former case the problem is typically easy to solve, and in the latter the analysis becomes much more technically challenging.

Fortunately, the \emph{maxima} of this problem have a more amenable structure
for study, as they are typically described by 1-RSB like structure. There is a
heuristic intuition for this: in the limit of $M\to1$, this problem is just the
square of a spherical spin glass landscape. The distribution and properties of
stationary points low and high in the spherical spin glass are not changed,
except that their energies are stretched and minima are transformed into
maxima. This is why the top of the landscape doesn't qualitatively change. The
bottom, however, consists of the zero-energy level set in the spherical spin
glass. This level set is well-connected, and so the ground states should also
be well connected and flat.

Focusing on the top of the landscape and therefore dealing with a 1-RSB like
problem is good for our analysis. First, algorithms will tend to be stuck in
the ways they are for hard optimization problems, and second  we will be able
to explicitly predict where. Therefore, we will study the most squares problem
rather than the least squares one. We calculate the complexity of maxima under a replica symmetric ansatz (which covers 1-RSB like problems) for arbitrary covariance $f$, and then the marginal complexity.

Applying the Lagrange multiplier method detailed above to enforce the spherical constraint, the gradient and Hessian are
\begin{align}
  \nabla H(\mathbf x,\omega)=\sum_k^MV_k(\mathbf x)\partial V_k(\mathbf x)+\omega\mathbf x
  \\
  \operatorname{Hess}H(\mathbf x,\omega)=\partial V_k(\mathbf x)\partial V_k(\mathbf x)+V_k(\mathbf x)\partial\partial V_k(\mathbf x)+\omega I
\end{align}
\begin{widetext}
The number of stationary points in a circumstance where the determinants add constructively is
\begin{equation}
  \begin{aligned}
    &\mathcal N(E,\mu)^n
    =\int\prod_{a=1}^nd\mathbf x_a\frac{d\hat{\mathbf x}_a}{(2\pi)^N}d\omega_a\,d\hat\beta_a\,\hat\mu_a\,d\bar\eta_a\,d\eta_a\,\exp\bigg\{
      i\hat{\mathbf x}_a^T(V^k(\mathbf x_a)\partial V^k(\mathbf x_a)+\omega\mathbf x_a)
      +\hat\beta(NE-\frac12V^k(\mathbf x_a)V^k(\mathbf x_a)) \\
    &  +\bar\eta_a^T(\partial V^k(\mathbf x_a)\partial V^k(\mathbf x_a)^T+V^k(\mathbf x_a)\partial\partial V^k(\mathbf x_a)+\omega_a I)\eta_a
    +\hat\mu_a(N\mu-\partial V^k(\mathbf x_a)^T\partial V^k(\mathbf x_a)-V^k(\mathbf x_a)\operatorname{Tr}\partial\partial V^k(\mathbf x_a)-N\omega_a)
    \bigg\}
  \end{aligned}
\end{equation}
To linearize the argument of the exponential with respect to $V$, we define the following new fields: $w^k_a=V^k(\mathbf x_a)$ and $\mathbf v^k_a=\partial V^k(\mathbf x_1)$. Inserting these in $\delta$ functions, we have
\begin{equation}
  \begin{aligned}
    &\mathcal N(E,\mu)^n
    =\int\prod_{a=1}^nd\mathbf x_a\frac{d\hat{\mathbf x}_a}{(2\pi)^N}d\omega_a\,d\hat\beta_a\,\hat\mu_a\,d\bar\eta_a\,d\eta_a\,\exp\bigg\{
      i\hat{\mathbf x}_a^T(w^k_a\mathbf v^k_a+\omega\mathbf x_a)
      +\hat\beta(NE-\frac12w^k_aw^k_a) \\
    &  +\bar\eta_a^T(\mathbf v^k_a(\mathbf v^k_a)^T+w^k_a\partial\partial V^k(\mathbf x_a)+\omega_a I)\eta_a
    +\hat\mu_a(N\mu-(\mathbf v^k_a)^T\mathbf v^k_a-w^k_a\operatorname{Tr}\partial\partial V^k(\mathbf x_a)-N\omega_a) \\
    & +i\hat w^k_a(w^k_a-V^k(\mathbf x_a))
    +i(\hat{\mathbf v}^k_a)^T(\mathbf v^k_a-\partial V^k(\mathbf x_a))
    \bigg\}
  \end{aligned}
\end{equation}
which is now linear in $V$. Averaging over $V$ yields, from only the terms that depend on it and to highest order in $N$,
\begin{equation}
  -\frac12\left(
    f(C_{ab})\hat w^k_a\hat w^k_b
    +2f'(C_{ab})\hat w^k_a\frac{\mathbf x^T_a\hat{\mathbf v}^k_b}N
    +f'(C_{ab})\frac{(\hat{\mathbf v}^k_a)^T\hat{\mathbf v}^k_b}N
    +f''(C_{ab})\left(\frac{\mathbf x_a^T\hat{\mathbf v}^k_b}N\right)^2
    +f''(C_{ab})w^k_aw^k_bG_{ab}^2
  \right)
\end{equation}
The resulting integrand is Gaussian in the $w$, $\hat w$, $\mathbf y$, and $\hat{\mathbf y}$, with
\begin{equation}
  \exp\left\{
    -\frac12\sum_{k=1}^M\sum_{ab}^n\begin{bmatrix}w_a^k\\\mathbf v_a^k\\\hat w_a^k\\\hat{\mathbf v}_a^k\end{bmatrix}^T
    \begin{bmatrix}
      \hat\beta_a\delta_{ab}+G_{ab}^2f''(C_{ab}) & -i\hat{\mathbf x}_a^T\delta_{ab} & -i\delta_{ab} & 0 \\
      -i\hat{\mathbf x}_a\delta_{ab} & 2(\hat\mu_a I-\bar\eta_a\eta_a^T)\delta_{ab} & 0 & -i\delta_{ab}I\\
      -i\delta_{ab} & 0 & f(C_{ab}) & \frac1Nf'(C_{ab})\mathbf x_a^T \\
      0 & -i\delta_{ab}I & \frac1Nf'(C_{ab})\mathbf x_b & \frac1Nf'(C_{ab})+\frac1{N^2}f''(C_{ab})\mathbf x_a\mathbf x_b^T
    \end{bmatrix}
    \begin{bmatrix}w_b^k\\\mathbf v_b^k\\\hat w_b^k\\\hat{\mathbf v}_b^k\end{bmatrix}
  \right\}
\end{equation}
which produces
\begin{equation}
  \exp\left\{
    \frac M2\log\det\left(
      I+\begin{bmatrix}
      \hat\beta_a\delta_{ab}+G_{ab}^2f''(C_{ab}) & -i\hat{\mathbf x}_a^T\delta_{ab} \\
      -i\hat{\mathbf x}_a\delta_{ab} & 2(\hat\mu_a I-\bar\eta_a\eta_a^T)\delta_{ab}
    \end{bmatrix}
    \begin{bmatrix}
      f(C_{ab})&\frac1Nf'(C_{ab})\mathbf x_a^T \\
      \frac1Nf'(C_{ab})\mathbf x_b & \frac1Nf'(C_{ab})+\frac1{N^2}f''(C_{ab})\mathbf x_a\mathbf x_b^T
    \end{bmatrix}
    \right)
  \right\}
\end{equation}
\begin{equation}
  \begin{bmatrix}
    (\hat\beta_a\delta_{ac}+G_{ac}^2f''(C_{ac}))f(C_{cb}) + R_{ab}f'(C_{ab})
    &
    \frac1N\left[(\hat\beta_a\delta_{ac}+G_{ac}^2f''(C_{ac}))f'(C_{cb})+R_{ab}f''(C_{ab})\right]\mathbf x_b^T-\frac1Nf'(C_{ab})\hat{\mathbf x}_a^T
    \\
    -i\hat{\mathbf x}_af(C_{ab})+\frac1N\hat\mu f'(C_{ab})\mathbf x_b
  \end{bmatrix}
\end{equation}

The condition fixing the maximum eigenvalue adds to the integrand
\begin{equation}
  \frac12\beta\sum_b^{m_a}\mathbf s^T_b(\mathbf v^k_a(\mathbf v^k_a)^T+w^k_a\partial\partial V^k(\mathbf x_a)+\omega I)\mathbf s_b
  +\frac12\hat\lambda\mathbf s_1^T(\mathbf v^k_a(\mathbf v^k_a)^T+w^k_a\partial\partial V^k(\mathbf x_a)+\omega I)\mathbf s_1
\end{equation}
\end{widetext}


\bibliography{marginal}

\end{document}